# 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?

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.