# sum of 2 squares

If the factorization of the integer $n$ is either $n=p_1^{a_1}\cdots p_s^{a_s}$ or $n=2p_1^{a_1}\cdots p_s^{a_s}$, where all $p_i$ are primes of form $4k+1$, then $n$ has $2^{s-1}$ representations as $n=a^2+b^2$, where $\gcd(a,b)=1$ (see for example this post ). The SAGE command sum_of_k_squares(2,n) will give one such representation. Is there a way to find all of them?

edit retag close merge delete

Sort by » oldest newest most voted

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]

more

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

more