Year-end appeal: Please make a donation to the OEIS Foundation to support ongoing development and maintenance of the OEIS. We are now in our 61st year, we have over 378,000 sequences, and we’ve reached 11,000 citations (which often say “discovered thanks to the OEIS”).
%I #27 Mar 23 2022 14:40:20
%S 1,2,3,5,6,8,9,11,15,21,30,39,50,63,83,95,99,156,173,350,854,1308,
%T 1769,2903,5250,5345,5639,6195,7239,21368,41669,47684,58619,63515,
%U 69468,70539,133508,134993,187160,493095
%N Numbers n such that 1 + Sum_{i=1..n} 2^(2i-1) is prime.
%C If this sequence is infinite then so is A124401.
%C Equals A127965(n)/2.
%C The sum has the simple closed form 1 + 2/3*(4^n-1). - _Stefan Steinerberger_, Nov 24 2007
%C Terms beyond a(30) correspond to probable primes, cf. A000978. - _M. F. Hasler_, Aug 29 2008
%F a(n) = floor(A000978(n)/2) = ceiling(log(4)(A000979(n))); A000978(n) = 2 a(n) + 1; A000979(n) = (2*4^a(n)+1)/3. - _M. F. Hasler_, Aug 29 2008
%e a(1)=1 because 1 + 2 = 3 is prime;
%e a(2)=2 because 1 + 2 + 2^3 = 11 is prime;
%e a(3)=3 because 1 + 2 + 2^3 + 2^5 = 43 is prime;
%e a(4)=5 because 1 + 2 + 2^3 + 2^5 + 2^7 + 2^9 = 683 is prime;
%e ...
%t a = {}; Do[If[PrimeQ[1 + Sum[2^(2n - 1), {n, 1, x}]], AppendTo[a, x]], {x, 1, 1000}]; a
%t b = {}; Do[c = 1 + Sum[2^(2n - 1), {n, 1, x}]; If[PrimeQ[c], AppendTo[b, c]], {x, 0, 1000}]; a = {}; Do[AppendTo[a, FromDigits[IntegerDigits[b[[x]], 2]]], {x, 1, Length[b]}]; d = {}; Do[AppendTo[d, (1/2)(DigitCount[a[[x]], 10, 0]+DigitCount[a[[x]], 10, 1]]), {x, 1, Length[a]}]; d
%t Position[Accumulate[2^(2*Range[1000]-1)],_?(PrimeQ[#+1]&)]//Flatten (* The program generates the first 21 terms of the sequence. To generate more, increase the Range constant. *) (* _Harvey P. Dale_, Mar 23 2022 *)
%o (PARI) for(n=1,999, ispseudoprime(2^(2*n+1)\3+1) & print1(n",")) \\ _M. F. Hasler_, Aug 29 2008
%o (Haskell)
%o import Data.List (findIndices)
%o a127936 n = a127936_list !! (n-1)
%o a127936_list = findIndices ((== 1) . a010051'') a007583_list
%o -- _Reinhard Zumkeller_, Mar 24 2013
%o (Python)
%o from sympy import isprime
%o A127936 = [i for i in range(1,10**3) if isprime(int('01'*i+'1',2))]
%o # _Chai Wah Wu_, Sep 05 2014
%Y Cf. A127962, A127963, A127964, A127965, A127961, A000979, A000978, A124400, A126614, A127955, A127956, A127957, A127958, A127936, A127936, A124401, A010051, A007583.
%K nonn,more
%O 1,2
%A _Artur Jasinski_, Feb 08 2007, Feb 09 2007
%E Edited by _N. J. A. Sloane_ at the suggestion of _Andrew S. Plewe_, Jun 11 2007
%E 2 more terms from _Stefan Steinerberger_, Nov 24 2007
%E 6 more terms from _Dmitry Kamenetsky_, Jul 12 2008
%E a(30)-a(40) calculated from A000978 by _M. F. Hasler_, Aug 29 2008