If we start with an integer n>0 with known (or computable) factorization n=\prod_{p\in P} p^{r(p)}, for some finite set of primes P, each prime being one modulo 4, and corresponding powers function p\to r(p), then the algorithm to go straightforward to the result is as follows. Write each p as a product (a+ib)(a-ib) of two primes in \Bbb R=\Bbb Z[i], 0 < a\le b. If p=2 the representation is (1+i)(1-i), and the factors are associated. Else the two factors are different primes. Now from each p^{r(p)} prime power factorization in \Bbb Z we take the corresponding factorization in \Bbb Z, (a+ib)^r(p)(a-ib)^r(p), and distribute it in all (r(p), p\ne2) possibilities in two conjugated factors e(p),f(p) of same norm, p^{r(p)}=e(p)f(p), e(p)=A+iB, f(p)=\overline{e(p)}=A-iB.
We obtain numerous factorizations
\prod e(p)\cdot \prod f(p)\ ,
and writing \prod e(p)=S+iT, \prod f(p)=S-iT the two numbers are relatively prime iff each choice p^{r(p)}=e(p)f(p) has relatively prime factors e(p),f(p).
So we have in essence only two choices for e(p), either
(a+ib)^r(p), or (a-ib)^r(p)$.
The algorithm is thus to start with an n, collect all prime powers, represent them as product of powers over \Bbb Z[i], build for each p the factor e(p) in two ways for p\ne 2, compute \prod e(p)=S+iT, then S^2+T^2=n. Among the many generated solutions only the half is needed in case there is a prime p\ne 2 in P, because the choices \prod e(p)=S+iT, and \prod e'(p)=S-iT, with e'=\bar e lead to the same solution. So for the first prime p\ne 2 we make one choice.
def reps(n):
"""Searching for all different representations n = S² + T², where n > 0 is an integer
"""
sols = [] # and we append solutions sol one by one when found
first_prime, choices, R = True, [], GaussianIntegers()
for p, mul in factor(n):
if p == 2:
choices.append( [(1 + i)^mul] )
elif p % 4 != 1:
return []
else:
pp = R(p).factor()[0][0] # first factor of p seen in ZZ[i]
if first_prime:
first_prime = False
choices.append( [pp^mul] )
else:
choices.append( [pp^mul, pp.conjugate()^mul] )
for choice in cartesian_product(choices):
S, T = R(prod(choice))
S, T = abs(S), abs(T)
if S > T: S, T = T, S
sols.append((S, T))
sols.sort()
return sols
Then with the above code:
sage: N = 10^40 + 1
sage: reps(N)
[(1, 100000000000000000000),
(19999999800000001, 99999998000000020000),
(21175073018521651199, 97732370700092384320),
(21194619068964748801, 97728133730883635680),
(27288575518388109025, 96204644618527795624),
(47041175529588244705, 88244704117552958824),
(47058823529411764705, 88235294117647058824),
(64675591816386108385, 76269704491365444424)]
sage: factor(N)
17 * 5070721 * 5882353 * 19721061166646717498359681
An other check:
sage: N = 2^90 + 1
sage: factor(N)
5^2 * 13 * 37 * 41 * 61 * 109 * 181 * 1321 * 54001 * 29247661
sage: factor(len(reps(N)))
2^9
And we have all 2^9 expected representation, the corresponding list of solutions is of course not shown. One more check:
sage: reps(5^20 * 13^20 * 17^20)
[(1120228314932443982452255202016, 2472109092323644501221860856287),
(1222652640782476248129183940513, 2423087980638188359909730969616),
(1808156181728664217638411829584, 2024056833293480968965675605087),
(1891237057569836360342883758113, 1946652828318445656133130477184)]
sage: factor(1120228314932443982452255202016^2 + 2472109092323644501221860856287^2)
5^20 * 13^20 * 17^20