Hi everyone!

As I continue working through Project Euler, I want to write a blog about another piece of general mathematical knowledge that is both interesting on its own and might be useful in some problems. Consider the following Diophantine equation:

We assume that we're given a specific number $$$z \in \mathbb Z$$$, and we need to check if there are $$$x, y \in \mathbb Z$$$ for which the identity above holds. Then, if such numbers exist we should find them and report. Example of a problem that might need it:

Timus — 1593. Square Country. Version 2. You're given an integer $$$n$$$. Find minimum $$$k$$$ such that $$$n = a_1^2+\dots+a_k^2$$$.

Tl;dr.

Let $$$z = 2^{k_0} p_1^{k_1} \dots p_n^{k_n} p_{n+1}^{k_{n+1}} \dots p_m^{k_m}$$$, where $$$p_1, \dots, p_n$$$ are different prime numbers with remainder $$$3$$$ modulo $$$4$$$, and $$$p_{n+1}, \dots, p_m$$$ are different prime numbers with remainder $$$1$$$ modulo $$$4$$$. Then there are two cases. If any of $$$k_{1}, \dots, k_n$$$ is odd, there are no solutions. Otherwise there is always a solution $$$z = x^2 + y^2$$$ that looks like

where $$$i^2=-1$$$ and $$$x_k^2+y_k^2 = p_k^2$$$ for $$$k$$$ from $$$n+1$$$ to $$$m$$$. For each $$$p_k$$$, to find such $$$x_k, y_k$$$ we need to find an integer $$$i$$$ such that $$$i^2 \equiv -1 \pmod{p}$$$, then find a minimum $$$x_k = i y_k \bmod p_k$$$ for $$$1 \leq y_k < \sqrt {p_k}$$$. This is doable in $$$O(\log p_k)$$$.

And if we want to count solutions, their number is given by Jakobi's two-square theorem: The number of ways of representing $$$z$$$ as the sum of two squares is $$$4(d_1(z) - d_3(z))$$$, where $$$d_k(z)$$$ is the number of divisors of $$$z$$$ that have remainder $$$k$$$ modulo $$$4$$$.

From exact equation to $$$\bmod z$$$

First of all, let's solve a bit easier problem. One obvious thing we can do is to take remainder modulo $$$z$$$ on both parts:

This relaxed version is equivalent to finding $$$x, y, k \in \mathbb Z$$$ such that

Hmmm... Remainders modulo arbitrary number $$$z$$$ is not the most pleasant thing to work on directly. But remainders modulo prime number $$$p$$$ are usually nice. On the other hand, if $$$p$$$ is some prime factor of $$$z$$$ and there is a solution for $$$z$$$, it means that there will as well be a solution for $$$p$$$ with $$$k=\frac{z}{p}$$$. So, let's assume $$$z$$$ to be prime, for now.

From arbitrary $$$z$$$ to prime $$$p$$$

Now, we have another equation

where $$$p$$$ is a prime number. What's good about prime numbers is that remainders modulo prime numbers form a field (i.e. they work very similarly to rationals, and we can expect similar results to hold). For $$$p=2$$$, there is a non-trivial solution $$$x=y=1$$$. What about odd numbers $$$p$$$? There are two cases to consider, as the remainder of $$$p$$$ modulo $$$4$$$ is either $$$1$$$ or $$$-1$$$.

Fermat's theorem on sums of two squares tells us that for an odd prime $$$p$$$, the solution exists if and only if $$$p$$$ has a remainder $$$1$$$ modulo $$$4$$$. Moreover, the sum of two squares theorem tells us that the number $$$z$$$ is expressible as a sum of two squares if and only if its prime decomposition does not have a term $$$p^k$$$, where $$$p \equiv -1 \pmod 4$$$, and $$$k$$$ is odd. Let's find out why.

$$$p \equiv 1 \pmod 4$$$

Of course, it's not yet clear why these two cases are important. Let's assume that there is an integer $$$i$$$ such that

that is there is a remainder modulo $$$p$$$ which behaves very similarly to imaginary unit from complex numbers. Then

This reduces the initial equation to a union of linear equations

For each $$$y$$$, except $$$y=0$$$, there are $$$2$$$ possible values of $$$x = \pm iy$$$, so there are a total of $$$2p+1$$$ solutions. Noteworthy, it is always possible to find a pair of solutions $$$(x,y)$$$ such that $$$1 \leq x, y < \sqrt p$$$, which means that $$$x^2 + y^2 = p$$$ is satisfied exactly.

How to find it? Find $$$i$$$, and consider the minimum value of $$$ik\bmod p$$$ among $$$1 \leq k < \sqrt p$$$. Due to pigeonhole principle, there will be $$$k_1 \neq k_2$$$ such that $$$i (k_1 - k_2) \bmod p \leq \sqrt p$$$. This is actually very similar to 102354I - From Modular to Rational!

Now, when does such $$$i$$$ exist and how to find it? It is known that remainders modulo $$$p$$$ have a primitive root $$$g$$$ such that its powers from $$$0$$$ to $$$p-2$$$ run through all possible remainders modulo $$$p$$$. Note that for odd $$$p$$$ it always holds that

Then, if such $$$i$$$ exists we should be able to find it from

Well, technically $$$-g^{\frac{p-1}{4}}$$$ also can be used as $$$i$$$, but it's not that important. What's important is that it is possible to do as above only when $$$p-1$$$ is divisible by $$$4$$$. In other words, when $$$p \equiv 1 \pmod 4$$$.

$$$p \equiv -1 \pmod 4$$$

Now, let's think about the other case. If there is no such $$$i$$$, we can introduce it! Now, we can formally consider numbers that look like $$$x+iy$$$, where $$$i$$$ is not a remainder modulo $$$p$$$. Numbers of this kind, if treated formally, also form a field. If you're familiar with field theory, I should mention that it is isomorphic to the Galois field $$$GF(p^2)$$$. If you're not familiar with it, ignore what I just wrote.

The thing now is that we can try to find all solutions in this new, extended field. And it reduces to the same union of equations

so for every $$$y$$$, the only possible solutions are $$$x = \pm iy$$$. The problem is, this time such $$$x$$$ would not be a remainder modulo $$$p$$$, unless $$$y=0$$$. Instead, it will be an "imaginary" solution. So, the only "real" solution is $$$x \equiv y \equiv 0 \pmod p$$$. It means that all solutions to

look like $$$x = px'$$$ and $$$y=py'$$$. Thus,

So, if $$$k$$$ is not divisible by $$$p$$$, there are no solutions. Otherwise $$$k=pk'$$$ reduces it to

after which similar argument could be applied. So, if $$$k$$$ is divisible by an odd power of such $$$p$$$, there are no solutions. We're only one step away from solving the whole $$$x^2+y^2=z$$$ problem now, assuming that we know the factorization of $$$z$$$.

Back to arbitrary $$$z$$$

Now we need to use one more fact from complex numbers. There, we can introduce a norm

Its crucial property for this task is that it is multiplicative, that is

This gives the Brahmagupta–Fibonacci identity

from which it follows that if we can represent $$$z$$$ as a product of several several numbers that are expressible as a sum of two squares, we can use the identity above to also express $$$z$$$ as a sum of two squares. In complex number terms, it means that we will find a complex number $$$x+iy$$$ such that $$$\|x+iy\| = z$$$.

Repeating what's written in tl'dr section, let $$$z = 2^{k_0} p_1^{k_1} \dots p_n^{k_n} p_{n+1}^{k_{n+1}} \dots p_m^{k_m}$$$, where $$$p_1, \dots, p_n$$$ are different prime numbers with remainder $$$3$$$ modulo $$$4$$$, and $$$p_{n+1}, \dots, p_m$$$ are different prime numbers with remainder $$$1$$$ modulo $$$4$$$. Then there are two cases. If any of $$$k_{1}, \dots, k_n$$$ is odd, there are no solutions. Otherwise there is always a solution $$$z = x^2 + y^2$$$ that looks like

where $$$i^2=-1$$$ and $$$x_k^2+y_k^2 = p_k^2$$$ for $$$k$$$ from $$$n+1$$$ to $$$m$$$. This result is tightly connected to the following ones:

Classification of Gaussian primes. A Gaussian integer $$$a+ib$$$ is prime if and only if either of the following is true:

• $$$a=0$$$ or $$$b=0$$$ and the non-zero number is prime with remainder $$$3$$$ modulo $$$4$$$,
• $$$a^2 + b^2$$$ is $$$2$$$ or a prime number with remainder $$$1$$$ modulo $$$4$$$.

The result also allows to factor Gaussian integers by representing their norm as a sum of two squares in the way suggested above.

Jakobi's two-square theorem. The number of ways of representing $$$z$$$ as the sum of two squares is $$$4(d_1(z) - d_3(z))$$$, where $$$d_k(z)$$$ is the number of divisors of $$$z$$$ that have remainder $$$k$$$ modulo $$$4$$$. If there is a prime divisor with remainder $$$3$$$ mod $$$4$$$, this number is equal to $$$0$$$. Otherwise, it is $$$4$$$ times the number of divisors of $$$p_{n+1}^{k_{n+1}} \dots p_m^{k_m}$$$, we may interpret it that the divisor decides how much of multipliers in the product would correspond to $$$x_{k} + iy_k$$$, and how much to $$$y_k + i x_k$$$, after which $$$4$$$ accounts for all the possible ways to multiply $$$x+iy$$$ with Gaussian integers units, i.e. with $$$1$$$, $$$-1$$$, $$$i$$$ or $$$-i$$$, which accounts for all possible combinations of $$$(\pm x, \pm y)$$$.