%I #22 Aug 19 2013 13:18:31
%S 6,12,18,36,66,48,132,96,264,288,528,1026,1032,384,2064,1152,2112,
%T 8196,1536,4224,16386,32772,8448,32784,65538,9216,16896,65568,131076,
%U 12288,18432,66048,262176,131328,262272,132096,524352,1048578,1048584
%N a(n) = 2^x+2^y where p(n) is the n-th prime of the form 4*k+1 and x, y is the unique integer solution to p(n) = x^2+y^2.
%C Gauss proved that any prime of the form 4*k+1 (A002144) is equal to the unique sum of two squares. This sequence identifies these summed squares and uniquely maps them into a decimal from which the two squares can be retrieved. The mapping is given by a(n) = 2^x+2^y where p(n) = x^2+y^2.
%H L. A. Butler, <a href="http://www.maths.bris.ac.uk/~malab/PDFs/2ndYearEssay.pdf">A Classification of Gaussian Primes</a>, School of Mathematics, Bristol Univ.
%e p(6) = 41 and 41 = 4^2 + 5^2 hence a(6) = 2^4+2^5 = 48. To retrieve the values 4 and 5 from 48 convert 48 to binary. The 1 bits (there are only ever two) select 2^4 and 2^5. So x, y are 4, 5.
%t next1m4prime[n1_] := (n2=n1+1; While[!PrimeQ[n2]||!Mod[n2, 4]==1, n2++]; n2); getbinmap[m1_] := (m2=m1; m3=Floor[Sqrt[m2]]; While[!IntegerQ@Sqrt[m2-m3^2], m3--]; 2^Sqrt[m2-m3^2] + 2^m3); SetAttributes[getbinmap, Listable]; getbinmap[Table[Nest[next1m4prime, 1, n], {n, 1, 100}]]
%Y Cf. A002144
%K nonn
%O 1,1
%A _Frank M Jackson_, Aug 19 2013