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.