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