A299103 Primes p = x^2 + y^2, not of the form z^2 + 1, such that 2^(x^2) == 2^(y^2) == 1 (mod p).

73, 433, 601, 673, 1801, 4513, 18433, 32377, 37633, 54001, 55201, 61681, 63901, 66529, 100801, 115201, 121369, 122921, 168781, 178481, 187417, 203617, 210913, 258721, 286721, 370661, 414721, 588061, 649657, 695701, 737537, 1781921, 3194101, 4674797, 4681801, 5039581, 6433561, 7593961, 7692697

Offset

1,1

Links

Max Alekseyev, Table of n, a(n) for n = 1..100

Example

For 73 = 3^2 + 8^2 we have 2^9 == 1 (mod 73) and 2^64 == 2 (mod 73).

Maple

f:= proc(p) local F, x, y;
  if not isprime(p) then return false fi;
  if issqr(p-1) then return false fi;
  F:= GaussInt:-GIfactors(p)[2];
  x, y:= (Re, Im)(F[1][1]);
  2 &^ (x^2) mod p = 1 or 2 &^ (y^2) mod p = 1
end proc:
select(f, [seq(i, i=5..10^7, 4)]); # Robert Israel, Feb 02 2018

Mathematica

pmax = 8000000;
xmax = Ceiling@ Sqrt@ pmax;
Table[p = x^2 + y^2; If[p <= pmax && PrimeQ[p] && (PowerMod[2, x^2, p] == 1 || PowerMod[2, y^2, p] == 1), p, Nothing], {x, 2, xmax}, {y, 2, x}] // Flatten // Union (* Jean-François Alcover, Oct 16 2020 *)

Prog

(PARI) B=bnfinit(x^2+1); { is_A299103(p) = my(z); if(p%4!=1 || issquare(p-1), return(0)); z=abs(Vec(bnfisintnorm(B, p)[1])); Mod(2, p)^(z[1]^2)==1 || Mod(2, p)^(z[2]^2)==1; } \\ Max Alekseyev, Feb 02 2018

Crossrefs

Cf. A002313, A002496.

Sequence in context: A200908 A031418 A142229 * A201961 A142145 A227340

Adjacent sequences: A299100 A299101 A299102 * A299104 A299105 A299106

Keyword

nonn

Author

Thomas Ordowski, Feb 02 2018

Extensions

Terms a(7) onward from Max Alekseyev, Feb 02 2018

Status

approved