Ask Your Question

All decompositions of a prime as a sum of four squares

asked 2012-05-08 20:51:59 -0600

parzan gravatar image

updated 2012-05-08 20:52:38 -0600

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 flag offensive 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.

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

1 answer

Sort by ยป oldest newest most voted

answered 2012-05-09 21:17:15 -0600

Mike Hansen gravatar image

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()]
    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
sage: 8*(389+1)
edit flag offensive delete link 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.

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

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools


Asked: 2012-05-08 20:51:59 -0600

Seen: 419 times

Last updated: May 09 '12