sum of 2 squares

edit

Here is a possible approach:

def all_two_squares(n):
    return [(abs(d[0]),abs(d[1])) for d in divisors(GaussianIntegers()(n)) if norm(d)==n]

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