# All decompositions of a prime as a sum of four squares

 1 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? asked May 09 '12 parzan 858 ● 3 ● 12 ● 30 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.kcrisman (May 09 '12)

 1 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  posted May 10 '12 Mike Hansen 3840 ● 21 ● 46 ● 84 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.kcrisman (May 10 '12)

[hide preview]

## Stats:

Seen: 59 times

Last updated: May 10 '12