cp's OEIS Frontend

This is a front-end for the Online Encyclopedia of Integer Sequences, made by Christian Perfect. The idea is to provide OEIS entries in non-ancient HTML, and then to think about how they're presented visually. The source code is on GitHub.

A182199 Largest integer N such that a^(2^k) + b^(2^k) for 1 <= k <= N is prime, where p = a^2 + b^2 is the n-th prime of the form 4m+1.

Original entry on oeis.org

4, 2, 3, 2, 2, 2, 2, 1, 2, 2, 1, 1, 1, 1, 2, 2, 2, 1, 2, 1, 1, 1, 1, 1, 2, 2, 1, 2, 2, 2, 1, 1, 2, 1, 2, 1, 2, 2, 2, 2, 1, 1, 1, 2, 1, 1, 2, 2, 1, 2, 2, 1, 1, 1, 1, 1, 2, 1, 1, 1, 2, 2, 2, 1, 2, 2, 2, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 2, 1, 1, 1, 1, 2, 1, 2, 1, 2
Offset: 1

Views

Author

Thomas Ordowski, Apr 20 2012

Keywords

Comments

a(1) corresponds to the first four Fermat primes. - Thomas Ordowski, Apr 22 2012
Schinzel's hypothesis H implies that there are arbitrarily large terms in the sequence. - Thomas Ordowski, Apr 26 2012
First value > 4 is a(102416) = 5 (corresponding to p = 2823521). - Robert Israel, May 28 2015
First value > 5 is a(4250044) = 6 (corresponding to p = 151062433). - Giovanni Resta, Jun 09 2015

Examples

			Let f(p,k) = a^(2^k)+b^(2^k), where f(p,1) = p is a prime of form 4k+1.
f(5,1) = 5, f(5,2) = 17, f(5,3) = 257, f(5,4) = 65537, f(5,5) = 641*6700417. So N = 4. Next prime of form 4k+1 is 13; N = 2. 17; N = 3. etc.
		

Programs

  • Maple
    N:= 10^4: # to get values corresponding to primes <= N
    Primes:= select(isprime,[4*i+1 $ i=1..floor((N-1)/4)]):
    G:= map(p -> [Re,Im](GaussInt:-GIfactors(p)[2][1][1]),Primes):
    f:= proc(ab) local j;
      for j from 2 do if not isprime(ab[1]^(2^j)+ab[2]^(2^j)) then return(j-1) fi od
    end proc:
    map(f,G); # Robert Israel, May 28 2015
  • Mathematica
    nn = 35; pr = {}; Do[p = a^2 + b^2; If[p < nn^2 && PrimeQ[p], AppendTo[pr, {p, a, b}]], {a, nn}, {b, a}]; pr = Sort[pr]; {jnk, a, b} = Transpose[pr]; Table[i = 1; While[PrimeQ[a[[n]]^2^i + b[[n]]^2^i], i++]; i - 1, {n, 2, Length[pr]}] (* T. D. Noe, Apr 24 2012 *)
  • PARI
    f(p)=my(s=lift(sqrt(Mod(-1, p))), x=p, t); if(s>p/2, s=p-s); while(s^2>p, t=s; s=x%s; x=t); s
    forprime(p=5,1e3,if(p%4==1,a=f(p);b=sqrtint(p-a^2);n=1; while(ispseudoprime(a^(2^n)+b^(2^n)),n++);print1(n-1", ")))
    \\ Charles R Greathouse IV, Apr 24 2012