All decompositions of a prime as a sum of four squares

Jacobi showed that a prime $p$ has $8(p+1)$ decompositions as a sum of four squares. sage.rings.arith.four_squares computes one such decomposition. Is there a way to compute all of them?

edit retag close merge delete

Somewhat surprisingly, even in Hardy and Wright I don't seem to find a quaternion decomposition or something else that gives a construction of them... maybe I didn't read it closely enough. It sounds like at least a little new code would be needed, not just a straight-up application of a theorem.

( 2012-05-09 03:02:15 -0600 )edit

Sort by ยป oldest newest most voted

It is straightforward to do it brute-force:

def sum_of_squares(n, number_of_squares):
if number_of_squares <= 1:
if n.is_square():
yield [n.isqrt()]
if n != 0:
yield [-n.isqrt()]
return
for i in range(n.isqrt(), -1, -1):
k = n - i**2
for rest in sum_of_squares(k, number_of_squares - 1):
yield [i] + rest
if i != 0:
yield [-i] + rest


You can check that the theorem holds:

sage: %time len(list(sum_of_squares(389,4)))
CPU times: user 0.03 s, sys: 0.00 s, total: 0.03 s
Wall time: 0.02 s
3120
sage: 8*(389+1)
3120

more

I sort of assumed the poster had already tried this. I do like the yields - think this could now be Cythonized and actually put in Sage? Presumably for numbers larger than 389 this would quickly get slow.

( 2012-05-10 02:19:45 -0600 )edit