%I #22 Jan 23 2022 20:04:22
%S 4,16,20,24,36,40,56,64,84,100,116,120,124,140,144,156,160,176,180,
%T 184,196,204,224,236,240,256,260,264,276,280,284,296,300,324,340,344,
%U 384,396,400,404,420,436,440,444,464,480,484,496,516,536,540,544,556,576
%N Numbers k such that every prime factor of k^2 + 1 is congruent to 1 (mod 8).
%C From _Robert Israel_, Mar 29 2020: (Start)
%C All terms are divisible by 4.
%C Includes all terms of A005574 that are divisible by 4. (End)
%H Robert Israel, <a href="/A216033/b216033.txt">Table of n, a(n) for n = 1..10000</a>
%e 64 is in the sequence because 64^2 + 1 = 17*241 and {17, 241} == 1 (mod 8).
%p with(numtheory):for n from 1 to 1000 do:x:=factorset(n^2+1):n1:=nops(x):s1:=0:for m from 1 to n1 do: if irem(x[m],8)=1 then s1:=s1+1:else fi:od:if s1=n1 then printf(`%d, `,n):else fi:od:
%p # Alternative:
%p select(n -> numtheory:-factorset(n^2+1) mod 8 = {1}, 4*[$1..1000]); # _Robert Israel_, Mar 29 2020
%t Select[Range[576], Union[Mod[Transpose[FactorInteger[#^2 + 1]][[1]], 8]] == {1} &] (* _T. D. Noe_, Aug 31 2012 *)
%t Select[Range[600],AllTrue[FactorInteger[#^2+1][[All,1]],Mod[#,8]==1&]&] (* Requires Mathematica version 10 or later *) (* _Harvey P. Dale_, Jan 31 2021 *)
%o (Magma) [n: n in [1..600] | forall{PrimeDivisors(n^2+1)[i]: i in [1..#PrimeDivisors(n^2+1)] | IsOne(PrimeDivisors(n^2+1)[i] mod 8)}]; // _Bruno Berselli_, Aug 30 2012
%Y Cf. A005574, A216032.
%K nonn
%O 1,1
%A _Michel Lagneau_, Aug 30 2012
|