By adamant, history, 4 months ago, 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:

$x^2 + y^2 = z.$

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

$x+ iy = (1+i)^{k_0} p_{1}^{k_{1}/2} \dots p_n^{k_n/2} (x_{n+1}+iy_{n+1})^{k_{n+1}} \dots (x_m + iy_m)^{k_m},$

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:

$x^2+y^2 \equiv 0 \pmod z.$

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

$x^2 + y^2 = kz.$

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

$x^2 + y^2 \equiv 0 \pmod{p},$

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

$i^2 \equiv -1 \pmod{p},$

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

$x^2 + y^2 \equiv (x+iy)(x-iy) \pmod p.$

This reduces the initial equation to a union of linear equations

$\left[ \begin{array}{ll} x+iy \equiv 0 \pmod p, \\ x-iy \equiv 0 \pmod p. \end{array} \right .$

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 - От модулей к рациональным!

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

$-1 \equiv g^{\frac{p-1}{2}} \pmod p.$

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

$i^2 \equiv g^{\frac{p-1}{2}} \pmod p \iff i \equiv g^{\frac{p-1}{4}} \pmod p.$

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

$\left[ \begin{array}{ll} x+iy \equiv 0 \pmod p, \\ x-iy \equiv 0 \pmod p, \end{array} \right .$

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

$x^2 + y^2 = kp$

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

$p^2 x'^2 + p^2 y'^2 = kp.$

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

$x'^2+y'^2 = k',$

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

$\|x+iy \| = (x+iy)(x-iy) = x^2+y^2.$

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

$\|(a+ib)(c+id)\| = \| a + ib \| \cdot \| c + id \| .$

This gives the Brahmagupta–Fibonacci identity

$(a^2+b^2)(c^2+d^2)= (ac-bd)^2 + (ad+bc)^2,$

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

$x+ iy = (1+i)^{k_0} p_{1}^{k_{1}/2} \dots p_n^{k_n/2} (x_{n+1}+iy_{n+1})^{k_{n+1}} \dots (x_m + iy_m)^{k_m},$

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$, it's $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 $1$, $-1$, $i$ or $-i$, which accounts for all possible combinations of $(\pm x, \pm y)$.

And, of course, there are similar results for $3$ squares and $4$ squares:  Comments (16)
 » what is i % p ?
•  » » 4 months ago, # ^ | ← Rev. 4 →   Any solution to $i^2 \equiv -1 \pmod p$. It can be found as $i \equiv g^{\frac{p-1}{4}} \pmod p$, where $g$ is a primitive root of $p$, if $p = 4n+1$. Otherwise, $i$ does not exist among remainders modulo $p$.
 » 4 months ago, # | ← Rev. 2 →   I proposed a trivial version of this problem on codechef once.
 » I got stuck on problem 66 "Diophantine equation" for a very long time two years ago, because I thought you could solve it without googling for anything, as almost all of the first project euler problems. Definitely couldn't :/. Good topic.
•  » » I think another blog is more appropriate for problem 66 :)
 » For those interested:Task related to two/three/four-square theorems: BOJ 17633Harder Task: BOJ 17646
•  » » That's very interesting! I wonder what's the construction for 3 and 4 squares is. Would it be sufficient to take one of squares at random until you get a number that factorizes into smaller number of squares?
•  » » » Construction for 1 square is trivial I will go over rest.Minimum number of squares needed is 4 : As we know a number cannot be represented as a sum of 3 squares if $4^a\left(8b+7\right) = n$ for any non-negative integer $a,b$ so a number that can be expressed as $8b+6$ has a way to express it as a sum of three powers. Assuming we can find a solution to $8b+6=x^2+y^2+z^2$ our answer will be $\left(x\cdot 2^a\right)^2+\left(y\cdot 2^a\right)^2+\left(z\cdot 2^a\right)^2+\left(2^a\right)^2$. Our next step will be finding the solution for 3 squares.Minimum number of squares needed is 3 : We can just use randomization like you mentioned above. In my testing generating a random number less than $10^{18}$ and than generating random squares until we find a number such that we can represent our original number minus the square as sum of two squares it took less than 100 random squares generated for $10^{5}$ tests I ran.I am very curious if there is a solution without randomization.
 » I and vipin were discussing about Pythagorean triplets yesterday. I remembered a problem where we had to find the number of integer solutions of $a + b^2 = c^2$ till $10^5$. Does anyone remember where this problem has appeared?
•  » » Not exactly $a + b^2 = c^2$, but $a^2 + b^2 = c^2$ appeared HERE.
•  » » Similar problem: https://codeforces.com/contest/1487/problem/D
 » 4 months ago, # | ← Rev. 2 →   Bro how is your math so good? Did you prepare for some math Olympiad? I am really intimidated by your blogs.
•  » » Please don't be :)I don't have that much of Math preparation, really. I just studied some basics in university, and I enjoy reading and/or writing about some stuff I find interesting from time to time.
 » 4 months ago, # | ← Rev. 2 →   Thanks for another beautiful article!I have a question regarding implementation part: in tl;dr you wrote 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_k}$, then find a minimum $x_k = i y_k \pmod{p_k}$ for $1 \le y_k < \sqrt{p_k}$. This is doable in $\mathcal{O}(\log p_k)$How do you suggest to implement it in the given complexity? I think I have a way to do it, but it is rather different from yours, so I am very interested, what did you mean.1) First of all, I would struggle to find a primitive root. According to wikipedia, picking random number and verifying, is it a primitive root or not, takes $\mathcal{O}(\log \log p \cdot \log p)$ time: $\frac{p}{\phi(p - 1)} = \mathcal{O}(\log \log p)$.On the other side, being primitive root is too strong property: for $p = 4k + 1$ consider any element $x \in \mathbb{Z}_p^*$. If $4 | ord(x)$, then $x^{\frac{ord(x)}{4}} = \pm i$ and $\mathbb{Z}_p^*$ contains at least $\frac{p - 1}{2} = \mathcal{O}(p)$ elements of that type (any odd power of a fixed primitive root fits). To avoid direct calculating of $ord(x)$ one can choose $i = x^{\frac{p - 1}{2^k}}$ for some $k$. Therefore, we can find $i$ with non-deterministic $\mathcal{O}(\log p)$.2) However, I am far more interested in this, second, point. Any way of finding primitive root seems to be fast enough in practice.How do you find such $y_k$ that minimizes $x_k$? I can propose a different way: let $a \in \mathbb{Z}$ be an element from previous item: $a^2 \equiv -1 \pmod{p}$. Then one can find $gcd(a + i, p)$ in gaussian integers $\mathbb{Z}[i]$. It can be shown that $\parallel gcd(a + i, p)\parallel= p$: $gcd(a + i, P) \cdot \overline{gcd(a + i, p)} = gcd(a + i, p) \cdot gcd(a - i, p) = gcd(a^2 + 1, p) = p$and you can find $gcd$ in gaussian integers, since it's a euclidian domain. Moreover, each time you replace $(a, b) \rightsquigarrow (b, a \% b)$ the norm of second term at least halves, so $gcd$ workd in $\mathcal{O}(\log p)$.
•  » » 4 months ago, # ^ | ← Rev. 2 →   For 1), I'd do it in standard way to find primitive root. But I like that your idea avoids computing $\operatorname{ord}(x)$, which requires knowing the factorization of $p-1$. On the other hand, raising random numbers to $\frac{p-1}{4}$ is fine too, as they're likely to be primitive roots.For 2), I expected kind of same solution as to Library Judge — Min of Mod of Linear. There are several ways to do it, you can e.g. find minimum element with a Euclidean-like process, or with continued fractions (see the solution to 102354I - От модулей к рациональным there).But I like your solution more, it is an unexpected result for me that $\|\gcd(a+i, p)\| = p$.
•  » » » Thank you, definitely would read the tutorials you provided.