%I #13 Oct 09 2014 16:26:08
%S 2,5,13,41,137,149,229,293,397,509,661,677,709,761,809,877,881,1217,
%T 1249,1277,1601,2053,2633,3637,3701,4481,4729,5101,5449,5749,5861,
%U 7121,7237,7517,8009,8089,8117,8377,9661,14869,14897,18229,19609,20369,20441,21493,22349,23917,24781,24977,25717
%N Primes dividing nonzero terms in A003095: the iterates of x^2 + 1 starting at 0.
%C Relative density in the primes is 0, see Jones theorem 5.5.
%H Charles R Greathouse IV, <a href="/A247981/b247981.txt">Table of n, a(n) for n = 1..500</a>
%H Rafe Jones, <a href="http://arxiv.org/abs/math/0612415">The density of prime divisors in the arithmetic dynamics of quadratic polynomials</a>, J. Lond. Math. Soc. (2) 78 (2) (2008), pp. 523-544.
%F a(n) << exp(k^n) for some constant k > 0, see Jones theorem 6.1. In particular this sequence is infinite. - _Charles R Greathouse IV_, Sep 28 2014
%e 2 and 13 are in the sequence since A003095(4) = 26. 3 is not in the sequence since it does not divide any member of A003095.
%t Select[Table[d=0; t=0; Do[t=Mod[t^2+1,Prime[j]]; If[t==0,d=1],{k,1,Prime[j]}]; If[d==1,Prime[j],0],{j,1,1000}],#!=0&] (* _Vaclav Kotesovec_, Oct 04 2014 *)
%o (PARI) is(p)=my(v=List([1]),t=1); while(t,t=(t^2+1)%p; for(i=1,#v, if(v[i]==t, return(0))); listput(v,t)); isprime(p)
%Y Cf. A003095, A248218, A248219.
%K nonn
%O 1,1
%A _Charles R Greathouse IV_, Sep 28 2014