# A search for integers which can be written as a sum of two squares in multiple ways

As part of a number theory hobby project, I'm looking for a computational way to enumerate all integers $n$ which can be written as a sum of two integer squares in three or more ways. The range of $n$ will likely get quite large ( up to $\approx 2^{64}$) so I'm looking for methods that let me specify a search range for $n$, say $11\cdot10^{13}$ to $12\cdot10^{14}$, and return any values such as

$$1193162282546 = 1043815^2+321889^2 = 1066211^2+237395^2 = 1076911^2+182825^2$$

My first thought of solving this problem is rather brute force:

Enumerating $every$ sum of two squares within a range $n_0..n_1$ by iterating over $n=a^2+b^2$ in a double loop over $1\le a \le\sqrt{n_1-1}$ and $\sqrt{n_0^2-a^2}\le b \le \sqrt{n_1^2-a^2}$ , then sorting the list and scanning it for duplicate values of 3 or more values. I think this will work, but it may be very inefficient, with the vast majority of the computation thrown out.

The puzzle for you: Is there a better computational method?

-
It is interesting to note that the asymptotics can be calculated analytically. In fact, the asymptotic quantity will be almost all integers of the form $a^2+b^2$. – Eric Naslund Jan 21 '13 at 1:40
Isn't it the case that a number can be written as a sum of two squares in multiple ways if and only if it has the form $pq$, where $p$ and $q$ are numbers that can be written as the sum of two squares? And the form of numbers of the latter type is well-known. If so, you can generate the numbers very efficiently, by representing them in prime-factorization form. – MJD Jan 21 '13 at 1:42
According to this, "A number $n$ is a sum of two squares if and only if all prime factors of $n$ of the form $4m+3$ have even exponent in the prime factorization of $n$." I'm not sure if this will make your actual computation faster or slower. – George V. Williams Jan 21 '13 at 1:45
MJD, correct, but not efficient. Factoring is much more effort than enumerating, and this is going to be done for trillions of values. Factoring is useful if we want to test a single $n$ value, though. – Rich Denver Jan 21 '13 at 2:10
@RichDenver I think you misunderstand MJD's comment. Generating prefactored numbers is much more efficient than factoring, and certainly seems competitive with enumeration. By analogy, if you wanted to generate a list of almost-primes, it would make far more sense to start from a list of primes, rather than perform any factorizations. – Erick Wong Jan 21 '13 at 2:17

Since you do not require relatively prime $x,y$ in $x^2 + y^2 = n,$ this has a very simple characterization, using the prime factorization of $n:$

(A) the exponent of 2 is irrelevant

(B) if $q \equiv 3 \pmod 4,$ the exponent must be even, otherwise this prime is also irrelevant as far as the count of $x^2 + y^2$ expressions

(C) there must either be $p^{4 + k}$ with $k \geq 0, \; \; p \equiv 1 \pmod 4,$ or $p_1^{1 + k_1} p_2^{2 + k_2}$ with $k_1, k_2 \geq 0, \; \; p_1, p_2 \equiv 1 \pmod 4,$ where we do not care whether $p_1$ or $p_2$ is larger, or $p_1^{1 + k_1} p_2^{1 + k_2} p_3^{1 + k_3}$ with $k_1, k_2, k_3 \geq 0, \; \; p_1, p_2, p_3 \equiv 1 \pmod 4.$ Here $p_1, p_2, p_3$ are distinct.

Right, so three distinct good primes does it, or one and another squared or better, or one to the fourth power or more.

EEEDDDIITTT: given your anti-factoring stance, the best way to do this is to spray things. Make a list of primes $1 \pmod 4$ up to an appropriate bound. Use that to make a list up to $1.2 \cdot 10^{15}$ of the successful numbers of type (C). Save that, and multiply by arbitrary $q^2$ with $q \equiv 3 \pmod 4$ and arbitrary $2^k,$ and keep the result if it is not too large.

I decided to run a little program, odd numbers up to $1 + 13^4$ that were expressible in at least three ways as $x^2 + y^2$ with $0 \leq x \leq y,$ and were not divisible by any prime $q \equiv 3 \pmod 4.$ Note that multiplying by 2 or by such $q^2$ does not change the number of representations. However, I found the numbers with at least three representations first, then i said only print out iff odd and not divisible by any $q \equiv 3 \pmod 4,$ and finally find the prime factorization for each and print that. So this little table experimentally confirms the patterns in part (C) above; I did not build those in.

325 = 5^2 * 13
425 = 5^2 * 17
625 = 5^4
725 = 5^2 * 29
845 = 5 * 13^2
925 = 5^2 * 37
1025 = 5^2 * 41
1105 = 5 * 13 * 17
1325 = 5^2 * 53
1445 = 5 * 17^2
1525 = 5^2 * 61
1625 = 5^3 * 13
1825 = 5^2 * 73
1885 = 5 * 13 * 29
2125 = 5^3 * 17
2225 = 5^2 * 89
2405 = 5 * 13 * 37
2425 = 5^2 * 97
2465 = 5 * 17 * 29
2525 = 5^2 * 101
2665 = 5 * 13 * 41
2725 = 5^2 * 109
2825 = 5^2 * 113
2873 = 13^2 * 17
3125 = 5^5
3145 = 5 * 17 * 37
3425 = 5^2 * 137
3445 = 5 * 13 * 53
3485 = 5 * 17 * 41
3625 = 5^3 * 29
3725 = 5^2 * 149
3757 = 13 * 17^2
3925 = 5^2 * 157
3965 = 5 * 13 * 61
4205 = 5 * 29^2
4225 = 5^2 * 13^2
4325 = 5^2 * 173
4505 = 5 * 17 * 53
4525 = 5^2 * 181
4625 = 5^3 * 37
4745 = 5 * 13 * 73
4825 = 5^2 * 193
4901 = 13^2 * 29
4925 = 5^2 * 197
5125 = 5^3 * 41
5185 = 5 * 17 * 61
5365 = 5 * 29 * 37
5525 = 5^2 * 13 * 17
5725 = 5^2 * 229
5785 = 5 * 13 * 89
5825 = 5^2 * 233
5945 = 5 * 29 * 41
6025 = 5^2 * 241
6205 = 5 * 17 * 73
6253 = 13^2 * 37
6305 = 5 * 13 * 97
6409 = 13 * 17 * 29
6425 = 5^2 * 257
6565 = 5 * 13 * 101
6625 = 5^3 * 53
6725 = 5^2 * 269
6845 = 5 * 37^2
6925 = 5^2 * 277
6929 = 13^2 * 41
7025 = 5^2 * 281
7085 = 5 * 13 * 109
7225 = 5^2 * 17^2
7325 = 5^2 * 293
7345 = 5 * 13 * 113
7565 = 5 * 17 * 89
7585 = 5 * 37 * 41
7625 = 5^3 * 61
7685 = 5 * 29 * 53
7825 = 5^2 * 313
7925 = 5^2 * 317
8125 = 5^4 * 13
8177 = 13 * 17 * 37
8245 = 5 * 17 * 97
8381 = 17^2 * 29
8405 = 5 * 41^2
8425 = 5^2 * 337
8585 = 5 * 17 * 101
8725 = 5^2 * 349
8825 = 5^2 * 353
8845 = 5 * 29 * 61
8905 = 5 * 13 * 137
8957 = 13^2 * 53
9061 = 13 * 17 * 41
9125 = 5^3 * 73
9265 = 5 * 17 * 109
9325 = 5^2 * 373
9425 = 5^2 * 13 * 29
9605 = 5 * 17 * 113
9685 = 5 * 13 * 149
9725 = 5^2 * 389
9805 = 5 * 37 * 53
9925 = 5^2 * 397
10025 = 5^2 * 401
10205 = 5 * 13 * 157
10225 = 5^2 * 409
10309 = 13^2 * 61
10525 = 5^2 * 421
10585 = 5 * 29 * 73
10625 = 5^4 * 17
10693 = 17^2 * 37
10825 = 5^2 * 433
10865 = 5 * 41 * 53
10933 = 13 * 29^2
10985 = 5 * 13^3
11125 = 5^3 * 89
11225 = 5^2 * 449
11245 = 5 * 13 * 173
11285 = 5 * 37 * 61
11425 = 5^2 * 457
11525 = 5^2 * 461
11645 = 5 * 17 * 137
11713 = 13 * 17 * 53
11765 = 5 * 13 * 181
11849 = 17^2 * 41
12025 = 5^2 * 13 * 37
12125 = 5^3 * 97
12325 = 5^2 * 17 * 29
12337 = 13^2 * 73
12505 = 5 * 41 * 61
12545 = 5 * 13 * 193
12625 = 5^3 * 101
12665 = 5 * 17 * 149
12725 = 5^2 * 509
12805 = 5 * 13 * 197
12905 = 5 * 29 * 89
13025 = 5^2 * 521
13325 = 5^2 * 13 * 41
13345 = 5 * 17 * 157
13481 = 13 * 17 * 61
13505 = 5 * 37 * 73
13525 = 5^2 * 541
13625 = 5^3 * 109
13925 = 5^2 * 557
13949 = 13 * 29 * 37
14045 = 5 * 53^2
14065 = 5 * 29 * 97
14125 = 5^3 * 113
14225 = 5^2 * 569
14297 = 17 * 29^2
14365 = 5 * 13^2 * 17
14425 = 5^2 * 577
14645 = 5 * 29 * 101
14705 = 5 * 17 * 173
14825 = 5^2 * 593
14885 = 5 * 13 * 229
14965 = 5 * 41 * 73
15025 = 5^2 * 601
15041 = 13^2 * 89
15145 = 5 * 13 * 233
15317 = 17^2 * 53
15325 = 5^2 * 613
15385 = 5 * 17 * 181
15425 = 5^2 * 617
15457 = 13 * 29 * 41
15625 = 5^6
15665 = 5 * 13 * 241
15725 = 5^2 * 17 * 37
15805 = 5 * 29 * 109
16025 = 5^2 * 641
16133 = 13 * 17 * 73
16165 = 5 * 53 * 61
16325 = 5^2 * 653
16385 = 5 * 29 * 113
16393 = 13^2 * 97
16405 = 5 * 17 * 193
16465 = 5 * 37 * 89
16525 = 5^2 * 661
16705 = 5 * 13 * 257
16745 = 5 * 17 * 197
16825 = 5^2 * 673
16925 = 5^2 * 677
17069 = 13^2 * 101
17125 = 5^3 * 137
17225 = 5^2 * 13 * 53
17425 = 5^2 * 17 * 41
17485 = 5 * 13 * 269
17525 = 5^2 * 701
17629 = 17^2 * 61
17725 = 5^2 * 709
17797 = 13 * 37^2
17945 = 5 * 37 * 97
18005 = 5 * 13 * 277
18125 = 5^4 * 29
18241 = 17 * 29 * 37
18245 = 5 * 41 * 89
18265 = 5 * 13 * 281
18325 = 5^2 * 733
18421 = 13^2 * 109
18605 = 5 * 61^2
18625 = 5^3 * 149
18685 = 5 * 37 * 101
18785 = 5 * 13 * 17^2
18925 = 5^2 * 757
19025 = 5^2 * 761
19045 = 5 * 13 * 293
19097 = 13^2 * 113
19225 = 5^2 * 769
19325 = 5^2 * 773
19345 = 5 * 53 * 73
19465 = 5 * 17 * 229
19625 = 5^3 * 157
19669 = 13 * 17 * 89
19721 = 13 * 37 * 41
19805 = 5 * 17 * 233
19825 = 5^2 * 13 * 61
19865 = 5 * 29 * 137
19885 = 5 * 41 * 97
19925 = 5^2 * 797
19981 = 13 * 29 * 53
20165 = 5 * 37 * 109
20213 = 17 * 29 * 41
20225 = 5^2 * 809
20345 = 5 * 13 * 313
20485 = 5 * 17 * 241
20525 = 5^2 * 821
20605 = 5 * 13 * 317
20705 = 5 * 41 * 101
20725 = 5^2 * 829
20905 = 5 * 37 * 113
21025 = 5^2 * 29^2
21097 = 17^2 * 73
21125 = 5^3 * 13^2
21325 = 5^2 * 853
21425 = 5^2 * 857
21437 = 13 * 17 * 97
21605 = 5 * 29 * 149
21625 = 5^3 * 173
21845 = 5 * 17 * 257
21853 = 13 * 41^2
21905 = 5 * 13 * 337
21925 = 5^2 * 877
22025 = 5^2 * 881
22265 = 5 * 61 * 73
22321 = 13 * 17 * 101
22345 = 5 * 41 * 109
22525 = 5^2 * 17 * 53
22625 = 5^3 * 181
22685 = 5 * 13 * 349
22765 = 5 * 29 * 157
22865 = 5 * 17 * 269
22945 = 5 * 13 * 353
22997 = 13 * 29 * 61
23125 = 5^4 * 37
23153 = 13^2 * 137
23165 = 5 * 41 * 113
23225 = 5^2 * 929
23273 = 17 * 37^2
23425 = 5^2 * 937
23525 = 5^2 * 941
23545 = 5 * 17 * 277
23585 = 5 * 53 * 89
23725 = 5^2 * 13 * 73
23825 = 5^2 * 953
23885 = 5 * 17 * 281
24089 = 13 * 17 * 109
24125 = 5^3 * 193
24245 = 5 * 13 * 373
24425 = 5^2 * 977
24505 = 5 * 13^2 * 29
24565 = 5 * 17^3
24625 = 5^3 * 197
24905 = 5 * 17 * 293
24925 = 5^2 * 997
24973 = 13 * 17 * 113
25085 = 5 * 29 * 173
25181 = 13^2 * 149
25225 = 5^2 * 1009
25285 = 5 * 13 * 389
25325 = 5^2 * 1013
25345 = 5 * 37 * 137
25493 = 13 * 37 * 53
25525 = 5^2 * 1021
25625 = 5^4 * 41
25705 = 5 * 53 * 97
25721 = 17^2 * 89
25789 = 17 * 37 * 41
25805 = 5 * 13 * 397
25825 = 5^2 * 1033
25925 = 5^2 * 17 * 61
26065 = 5 * 13 * 401
26129 = 17 * 29 * 53
26225 = 5^2 * 1049
26245 = 5 * 29 * 181
26525 = 5^2 * 1061
26533 = 13^2 * 157
26585 = 5 * 13 * 409
26605 = 5 * 17 * 313
26645 = 5 * 73^2
26725 = 5^2 * 1069
26765 = 5 * 53 * 101
26825 = 5^2 * 29 * 37
26945 = 5 * 17 * 317
27145 = 5 * 61 * 89
27325 = 5^2 * 1093
27365 = 5 * 13 * 421
27425 = 5^2 * 1097
27521 = 13 * 29 * 73
27565 = 5 * 37 * 149
27625 = 5^3 * 13 * 17
27725 = 5^2 * 1109
27925 = 5^2 * 1117
27985 = 5 * 29 * 193
28033 = 17^2 * 97
28085 = 5 * 41 * 137
28145 = 5 * 13 * 433
28225 = 5^2 * 1129
28249 = 13 * 41 * 53
28561 = 13^4


# ---------------------------

-
This strategy will work, combined with Eric's point that I don't need to factor each n, I need to GENERATE n progressively with new factors. This is the effectively the "spray" method. Thanks! – Rich Denver Jan 21 '13 at 3:21