rcs> Your solution comes from (3^64 + i)/(1+i).
Whenever X = 2N+1 is odd, (X+i)/(1+i) generates the decomposition (X2 + 1)/2 = 2N2 + 2N + 1 = (N+1)^2 + N^2.
The 128 solutions: Reduce is obviously counting all orders and all sign combinations as different. In this scheme, 5 = 4+1 = 2^2+1^2 has 8 solutions as a2+b2. Right. To just count them: In[16]:= SquaresR[2, (3^128 + 1)/2]
Out[16]= 128
Normally, we'd say that 5 splits as the sum of two squares in 1 way.
Right. To get just those, In[17]:= Length[PowersRepresentations[(3^128 + 1)/2, 2, 2]] Out[17]= 16 --rwg (As always, docfiend NeilB is showing me these things.)
Now 3^128+1 = 2 257 275201 138424618868737 3913786281514524929 153849834853910661121
Ignoring the 2, each of these primes is 256K+1, and hence 4J+1, and splits (uniquely) into two squares. E.g. 257 = 16^2 + 1^2. Each additional 4J+1 prime doubles the number of solutions, as we may choose either sign for the imaginary terms in the product (16 +- i) (524 +- 25i) (10663711 +- 4970904i) (1503205840 +- 1286140927i) (11066340060 +- 5602316711i).
--- Slightly related: I once went looking for Fermatish bases N, so that N+1, N2+1, N4+1, N8+1, &c are all prime. For N=2: 2^2^K + 1 is prime for K=0-4, but composite when K=5. For N=6: 7, 37, and 1297 are prime, but 6^8+1 = 17*98801. I found a few N which worked for K=0-5, but no K=0-6 primes. N = 7072833120, 9736020616, and 12852419340 work for K=0...5, but generate composites when K=6. Standard conjectures assert that any finite set of Ks should have an infinite set of prime-generating bases N.
Gee, and you can't even find six in a row? Tsk.
Rich
------ Quoting Bill Gosper <billgosper@gmail.com>:
Reduce[(3^128 + 1)/2 == a^2 + b^2, {b, a}, Integers] finds 128 solutions, including (b == 1716841910146256242328924544640 && a == 1716841910146256242328924544641) --rwg