# What is the probability that a natural number is a sum of two squares?

Some natural numbers can be expressed as a sum of two squares:

$$2=1^2+1^2$$ $$25=3^2+4^2$$ $$50=7^2+1^2$$

If one chooses a random natural number, what would be the probability that that number is a sum of two squares? Is it zero?

I read about Lagrange´s theorem on squares, but it looks it can´t be useful here.

NOTE 1: "Square" means "square of a natural number".

NOTE 2: I am aware that the expression "random natural number" is not a strict math notion. However, as I said in a comment, one can adopt a reasonable strict definition, which is not difficult to devise at all. It is mentioned also in an answer below.

NOTE 3: A related question on SE: How to determine whether a number can be written as a sum of two squares?

• What probability distribution are you using? Are you sure such a distribution exists?
– user223391
Aug 30 '15 at 8:55
• Aug 30 '15 at 9:00
• @avid19 I leave to you and your intuition the interpretation of expression "random natural number". One possible approach would be first to limit choice to finite number of natural numbers (lets say first $N$), use uniform distribution, and see how it behaves for $N$ taking large values. Aug 30 '15 at 9:00
• Perhaps also of some interest: the probability that a natural number is a sum of three squares is 5/6. Sep 9 '15 at 14:16
• @MichaelLugo: Yes, I've conjectured (in the answer below) that it is $5/6$. (Thanks: I made an error) But how can we prove it? Sep 9 '15 at 20:08

Without going into detail what a "random natural number" might be, we could consider the density of such numbers.

The possibility to express $n$ depends on the prime factorization of $n$: It may be divisible by $2$ and by primes $\equiv 1\pmod 4$ as much as it likes, but each prime $\equiv 3\pmod 4$ must occur in an even power.

The prime $3$ "spoils" all numbers divisible by $3$ (that's $\frac13$), except those divisible by $9$ (that's $\frac19$), though it does spoil those divisible by $27$, execept ... All in all the prime $3$ spoils $\frac13-\frac19+\frac1{27}-\frac1{81}\pm\ldots = \frac14$ (geometric series) of all numbers in a large range. In general, a prime $p\equiv 3\pmod 4$ spoils $\frac1{p+1}$ and the effects of distinct primes are independent. Hence the limit density of natural numbers expressible as sum of squares is $$\prod_{p\equiv 3\pmod 4}\left(1-\frac1{p+1}\right)$$ This product however diverges to $0$ because it is well-known that $\sum_p\frac1p$ diverges (where it doesn't matter that we use only half of all primes).

Thus: For sufficiently large $N$, the probability that a number picked uniformly from $\{1,\ldots,N\}$ is the sum of two squares becomes arbitrarily small. It is for any $\let\epsilon\varepsilon\epsilon>0$ we can find $N_0$ such that for all $N>N_0$ said probability is $<\epsilon$.

• Nice answer! It is worth to add that the average number of solutions to the equation $x^2 + y^2 = n$ is positive (something like $\pi/4$). Why is this interesting? Because this is kind of ''paradox'': this ''random variable'' (number of solutions) is almost surely zero, but has a positive expectation! Aug 30 '15 at 9:21
• @AdityaAgarwal, are you familiar with the concept of a measure space? $\mathbb{Q}$ is infinite but has zero measure as a subset of $\mathbb{R}$. It's possible for things to be infinitely big and yet "take up no space at all". Aug 30 '15 at 10:09
• I understood, sorry for the confusion. Aug 30 '15 at 10:18
• You say that half the primes satisfy $p \equiv 3\pmod 4$. That sounds plausible, but has it been proven? Aug 30 '15 at 15:53
• @user7530 That's exactly why they are independent. See, of the numbers $1,\ldots,210000$ exactly $70000$ are divisible by $3$, exactly $30000$ are divisible by $7$ and exactly $10000$ are divisible by both. Hence exactly $\frac17$ of those divisible by $3$ are additionally divisible by $7$, and of those $140000$ not divisible by $3$ exactly $\frac17$ ($20000$) are divisible by $7$. Hence divisibility by $7$ is completely independent of divisibility by $3$. Thus a share $(1-\frac13)(1-\frac17)$ of these $210000$ numbers is divisible by neither $3$ nor $7$. Aug 31 '15 at 20:05

Though this is not an answer to the original question, I'll write it in answer to VividD's question under Hagen von Eitzen post (and I believe it is not totally unrelated to the question).

Let $A_n$ be the number of different pairs $(x,y)$ of non-negative integers solving the equation $x^2 + y^2 = n$. Then it is natural to define the expected number of compositions of a randomly chosen natural number into sum of two squares as
$$A = \lim_{n\to\infty} \frac1n \sum_{k=1}^n A_k.$$ But the latter quantity is $\frac1n$ times the number of points of the form $(x/\sqrt{n},y/\sqrt{n})$, where $x,y$ are integers, inside the non-negative quarter of the unit circle centered at the origin (minus one, but this does not matter). So here we have a Jordan approximation of the area of the quarter. Therefore, $A = \frac\pi4$.

So, quite paradoxically, the expected value of the number of compositions of a randomly chosen non-negative integer into two squares is positive despite this number is zero "with probability 1".

• How can this be extended to the variant of the original question for sum of three and/or four squares? Does these cases correspond to 3 and 4 dmensional space? What would be the expectation in these cases? Aug 30 '15 at 12:14
• @VividD, no, for squares it will not work. The reason why it works in $\mathbb R^2$ is that $1/n = (1/\sqrt{n})^2$, the volume of a square in the Jordan approximation. So the expectation will be infinite for more than three squares. But, for the same reason, it will work for three cubes, four 4th powers etc. I'm too lazy to compute the volumes of those objects, sorry :) Aug 30 '15 at 14:02
• @zhoraster: My answer is based upon your argument and confirms it (apart from a factor $1/2$). Sep 5 '15 at 13:22

Consider the function $r: \mathbb{N} \to \mathbb{C}$: $$r(n) = \begin{cases} 1 &n \text{ is the sum of two squares (0 is a square)} \\ 0 &\text{otherwise}. \end{cases}$$ One way to define the "probability that a random integer is the sum of two squares" would be to consider the distribution on the integers where $n$ selected with probability proportional to $n^x$, for $x > 1$, and then to take the limit as $x \to 1$. That is, we can try to compute: $$\lim_{x \to 1} \frac{\sum_{i=1}^\infty \frac{r(n)}{n^x}}{\sum_{i=1}^\infty \frac{1}{n^x}} \tag{1}.$$

As mentioned here, $r(n)$ is $1$ iff every prime of the form $4k-1$ occurs an even number of times in $n$. It follows that $r(n)$ is multiplicative, and \begin{align*} \frac{\sum_{i=1}^\infty \frac{r(n)}{n^x}}{\sum_{i=1}^\infty \frac{1}{n^x}} &= \frac{\prod_{p \text{ prime}} \sum_{i=0}^\infty \frac{r(p^i)}{p^{ix}}}{\prod_{p \text{ prime}} \sum_{i=0}^\infty \frac{1}{p^{ix}}} \\ &= \prod_{p \text{ prime}} \frac{1 + \frac{r(p)}{p^x} + \frac{r(p^2)}{p^{2x}} + \cdots}{1 + \frac{1}{p^x} + \frac{1}{p^{2x}} + \cdots} \\ &= \prod_{p \equiv 3 \pmod{4}} \frac{1 + 0 + \frac{1}{p^{2x}} + 0 + \frac{1}{p^{4x}} + \cdots}{1 + \frac{1}{p^x} + \frac{1}{p^{2x}} + \cdots} \prod_{p \equiv 1,2 \pmod{4}} \frac{1 + \frac{1}{p^x} + \frac{1}{p^{2x}} + \cdots}{1 + \frac{1}{p^x} + \frac{1}{p^{2x}} + \cdots} \\ &= \prod_{p \equiv 3 \pmod{4}} \frac{1 + \frac{1}{p^{2x}} + \frac{1}{p^{4x}} + \cdots}{1 + \frac{1}{p^x} + \frac{1}{p^{2x}} + \cdots} \\ &= \prod_{p \equiv 3 \pmod{4}} \frac{1 \big/ \left(1 - \tfrac{1}{p^{2x}}\right)}{1 \big/ \left(1 - \tfrac{1}{p^{x}}\right)} \\ &= \prod_{p \equiv 3 \pmod{4}} \frac{1}{1 + p^{-x}}. \\ \end{align*} Now if you plug in $x = 1$, you get the product $$\prod_{p \equiv 3 \pmod{4}} \frac{p}{p + 1} = \prod_{p \equiv 3 \pmod{4}} \left( 1 - \frac{1}{p + 1} \right) = 0,$$ by the reasoning of Hagen von Eitzen. So by this definition of probability, the probability that a random integer is a sum of two squares is zero.

• i repeat .The sum of the reciprocals of all primes is infinite,This does not automatically imply the reciprocals of those congruent to 3 mod 4 sums to infinity. And you have an inconsistency about the function r . As defined,r(n)=0 if n is a power of 4. But according to the sentence after formula (1) it is 1,not 0. Sep 10 '15 at 15:36
• @user254665 Thanks for your criticism. $r(n) = 1$ for all powers of 4, and in general for all perfect squares. $4^n = (2^n)^2 + 0^2$. Your point with primes $\equiv 3 \mod 4$ is true--it's not as obvious as Hagen's answer states. But it's still true. I suspect Hagen or someone else on this site could tell you how the reasoning goes to prove that's infinite. Sometimes, a good answer leaves out the proof for some facts, and instead just cites them.
– 6005
Sep 10 '15 at 15:48
• It is a theorem that asymptotically there are as many primes 4n+1 as 4n-1 and it's certainly not elementary, even assuming the prime number theorem . Sep 19 '15 at 6:34
• @user254665 Do you perhaps have an asymptotic answer to my question here?
– 6005
Sep 19 '15 at 10:18

the density is zero, and one may be quite precise about it: the count of numbers up to some large real $x$ that are the sum of two squares is asymptotic to $$\frac{0.7642... \, x}{\sqrt{\log x}}$$ where the logarithm is base $e,$ and the $0.7642...$ is defined by an infinite product. See the last few pages in LeVeque. This is combined volumes 1 and 2, it is the last few pages in volume 2.

A simple numerical experiment shall confirm the answer given by zhoraster: the latter is only deviant from mine by a factor $2$. We define an initial segment of the naturals with length $N$ and count all sums of two squares in that segment. The accompanying program is in Pascal:

program kwadraat;
procedure test(N : integer);
type
data = record
b : boolean; i,j : integer;
end;
var
i,j,k,t : integer;
rij : array of data;
begin
t := 0;
SetLength(rij,N);
for k := 0 to N-1 do
begin
rij[k].b := false;
end;
i := 0;
while true do
begin
if sqr(i) > N-1 then Break;
j := i;
while true do
begin
k := sqr(i)+sqr(j);
if k > N-1 then Break;
{ Writeln(k,' = ',i,'^2 + ',j,'^2'); }
rij[k].i := i; rij[k].j := j;
rij[k].b := true;
j := j + 1;
t := t + 1;
end;
i := i + 1;
end;
if N < 100 then
for k := 0 to N-1 do
begin
if rij[k].b then
Writeln(k,' = ',rij[k].i,'^2 + ',rij[k].j,'^2');
end;
Writeln(t/N,' ->',Pi/8);
end;
begin
test(10);
test(1000);
test(100000);
test(10000000);
end.

Output (details for $N=10$ only):

0 = 0^2 + 0^2
1 = 0^2 + 1^2
2 = 1^2 + 1^2
4 = 0^2 + 2^2
5 = 1^2 + 2^2
8 = 2^2 + 2^2
9 = 0^2 + 3^2
7.00000000000000E-0001 -> 3.92699081698724E-0001
4.19000000000000E-0001 -> 3.92699081698724E-0001
3.95420000000000E-0001 -> 3.92699081698724E-0001
3.92969900000000E-0001 -> 3.92699081698724E-0001

Note that the results converge to $\;\pi/8$ , quite in agreement with the argument given by zhoraster, provided though that $i^2 + j^2$ and $j^2 + i^2$ give a double count of $\,\pi/4\,$ which must be halved.

EDIT. Question & Answer is related to : Double Think about Numerosity .

BONUS. In one of the comments with the answer by zhoraster, VividD has been asking for a variant of the original question for the sum of three/four squares. Minor modification of the above program gives the following output for the three squares case. It is seen that some numbers can be written as a sum of three squares in more than one way. Therefore two cases shall be distinguished: with or without these duplicates. Details again for $N=10$ :

0 = 0^2 + 0^2 + 0^2
1 = 0^2 + 0^2 + 1^2
4 = 0^2 + 0^2 + 2^2
9 = 0^2 + 0^2 + 3^2
2 = 0^2 + 1^2 + 1^2
5 = 0^2 + 1^2 + 2^2
8 = 0^2 + 2^2 + 2^2
3 = 1^2 + 1^2 + 1^2
6 = 1^2 + 1^2 + 2^2
9 = 1^2 + 2^2 + 2^2
with duplicates = 10/10
without duplicates = 9/10
with duplicates = 3254/1000
without duplicates = 835/1000
with duplicates = 2807201/100000
without duplicates = 83336/100000
with duplicates = 87741031/1000000
without duplicates = 833336/1000000

If the duplicates are counted, then the results are seen to diverge $\to \infty$ .
It is conjectured that, without the duplicates, the results converge to : $5/6$ (Michael Lugo).
Any takers to prove the latter statement?

• I don't see why you say this throws doubt on the other answers; if I understand what you've done here, you've computed the expected number of ways that an integer may be written as a sum of two squares, not the expected number of integers that may be so written. These are different problems, and have different answers. Sep 9 '15 at 13:42
• @SimonRose: If that is so, how can I conjecture then the comment by Michael Lugo? An integer may not be written as the sum of 2 squares, as the sum of 2 squares uniquely, as the sum of 2 squares in two ways: $$25 = 0^2 + 5^2 = 3^2 + 4^2 \\ 50 = 1^2 + 7^2 = 5^2 + 5^2 \\ 65 = 1^2 + 8^2 = 4^2 + 7^2 \\ 85 = 2^2 + 9^2 = 6^2 + 7^2$$ And so on and so forth. So I don't understand what you mean. Sep 9 '15 at 20:25
• You have in fact confirmed the answer by zhoraster, but it seems you did not understand it. Even that answer itself says it is "not an answer to the original question," that is, it comes up with a different value than the other answers because it is answering a different question. It does not contradict the other answers. Sep 9 '15 at 21:25
• @DavidK: Would you please elaborate about the difference between the two questions, because indeed it seems that I do not understand it. For my own comprehension, if the sums of two squares are counted without the duplicates, then I get much lower values. Are these converging to zero perhaps? Sep 10 '15 at 10:29
• Yes, the difference between counting duplicates and not counting duplicates is that without duplicates the ratio converges to zero, and with duplicates it converges to $\pi/8$ or $\pi/4$, depending on whether you distinguish $3^2+4^2$ from $4^2+3^2$. If $X$ is the number of ways that the chosen number can be written as two squares, $X>0$ when the number is the sum of two squares and OP asks for $P(X>0)$. If you count duplicates you compute $E(X)$. Sep 10 '15 at 12:18