A deterministic polynomial-time algorithm for finding a representation of a prime as a sum of four squares
- Clone this repository.
$ git clone https://github.com/c910335/four-squares
$ cd four-squares
- Compile it.
$ shards build --release --no-debug
Run and enter some primes.
$ bin/four_squares
1
1 is not a prime.
2
1^2 + 1^2 + 0^2 + 0^2 = 2
65539
233^2 + 75^2 + 75^2 + 0^2 = 65539
2147483647
33135^2 + 31941^2 + 4310^2 + 3279^2 = 2147483647
170141183460469231731687303715884105727
9485443033339641291^2 + 7911815779561688934^2 + 4104118515834236449^2 + 852605535991260783^2 = 170141183460469231731687303715884105727
The algorithm is from Sums of Four Squares.
- Given a prime P, define P - 1 = Q2S with Q odd; define the level of a number n to be the smallest k that satisfies n2kQ ≡ 1 (mod P).
- Find the first element whose level is not 1 in the following sequence.
- c0 = P - 1
- ci = ci-1 / 2 if ci-1 is even.
- ci = ci-1 - 1 if ci-1 is odd.
- Let cj be the first element whose level is not 1:
- If the level of cj is greater than 1, then P - 1 is a quadratic residue modulo P, and its square root x satisfies x2 + 1 = kP for some k.
- If the level of cj is 0, then both cj and P - cj-1 are quadratic residues modulo P, and their square roots x and y satisfy x2 + y2 + x2(or 1 if cj-1 is odd) = kP for some k.
- Find the square roots x, y, z that satisfy the equation x2 + y2 + z2 = kP in the previous step by the Tonelli–Shanks algorithm.
- Find the gcd (a + bi + cj + dk) of the two Hurwitz integers (x + yi + zj) and P by the Euclidean algorithm.
- If their gcd is not a Lipschitz integer, left multiply it by (e1 - e2i - e3j - e4k) / 2, where e1 ≡ a * 2, e2 ≡ b * 2, e3 ≡ c * 2, e4 ≡ d * 2 (mod 4), and e1, e2, e3, e4 ∈ {±1}.
- Output a2 + b2 + c2 + d2 = P.
If the input number is not a prime, the algorithm fails with high probability and thus we can know that it is a composite number.
However, it is known that a lot of composite numbers pass this test.
Here are some examples.
51^2 + 26^2 + 0^2 + 0^2 = 3277 # 29 * 113
63^2 + 8^2 + 0^2 + 0^2 = 4033 # 37 * 109
49^2 + 40^2 + 28^2 + 26^2 = 5461 # 43 * 127
65^2 + 64^2 + 0^2 + 0^2 = 8321 # 53 * 157
171^2 + 10^2 + 0^2 + 0^2 = 29341 # 13 * 37 * 61
210^2 + 71^2 + 0^2 + 0^2 = 49141 # 157 * 313
255^2 + 16^2 + 0^2 + 0^2 = 65281 # 97 * 673
- Fork it (https://github.com/c910335/four-squares/fork)
- Create your feature branch (
git checkout -b my-new-feature
) - Commit your changes (
git commit -am 'Add some feature'
) - Push to the branch (
git push origin my-new-feature
) - Create a new Pull Request
- Tatsujin Chin - creator and maintainer