text: 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
  Asked: 2012-05-08 20:51:59 -0600
  Seen: 422 times
  Last updated: May 09 '12
  OrderedSetPartitions (variant)
  Strings with given frequency
  Row reduction modulo prime powers
  Using `SimplicalComplex` on a set of sets
  Number of (and way to enumerate) head to head schedules?
  symmetric skew Macdonald polynomials
  4-Digit Rounding Arithmetic:System of Equations
  How to draw special trees from a list consisting of tuples with Sage?
  The value of binomial(-1,-1) is?
  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.