Ask Your Question
0

Finding representation of a given number as a sum of squares

asked 2023-04-09 20:01:04 +0200

klaaa gravatar image

updated 2023-04-09 20:03:49 +0200

Is there a quick method using Sage to find all representations of a given (possible very large) positive integer B as a sum $x_1^2+....+x_n^2$ of exactly $n$ integers $x_i >0$, where we can assume that $x_i \leq x_{i+1}$ for all $i$?

Thanks for any help

edit retag flag offensive close merge delete

2 Answers

Sort by ยป oldest newest most voted
1

answered 2023-04-11 15:59:23 +0200

dan_fulea gravatar image

Let us implement the function with bare hands.

def my_reps(k, n):
    if k <= 0 or n == 0:
        return []
    if k == 1:
        return [[isqrt(n)]] if n.is_square() else []
    sols = []
    for xn in [isqrt(n/k)..isqrt(n - k + 1)]:
       shorter_sols = [sol for sol in my_reps(k - 1, n - xn^2) if sol[0] <= xn]  
       sols += [[xn] + s for s in shorter_sols]
    return sols

The cheap recursion makes the code less optimal in cases of a "big $k$". We can check by counting, if the function is doing the right job. We have for instance:

sage: my_reps(3, 1001) [[24, 16, 13], [24, 19, 8], [24, 20, 5], [26, 15, 10], [26, 17, 6], [26, 18, 1], [27, 16, 4], [29, 12, 4], [30, 10, 1], [31, 6, 2]] sage: len(_) 10

Some counting generating function for the solutions is $f^k$, where $f=q+q^4+q^9+q^{16}+q^{25}+\dots$ - but order of components matter. We have collected only decreasing solutions. Fortunately, all displayed solutions have three different components, so we expect a $10\cdot 3!$ count in $f^3$ as coefficient of $q^{1001}$:

sage: R. = PowerSeriesRing(QQ, default_prec=1600, sparse=True) sage: R. = PowerSeriesRing(QQ, default_prec=1600) sage: f = sum([q^(n^2) for n in [1..40]]) + O(q^1601)

And indeed, there is a part with 12*q^1000 + 60*q^1001 + 24*q^1002 + 12*q^1003 + 21*q^1004 + 24*q^1005 + 30*q^1006 + 27*q^1009 in the answer. For $n=1588$ there are particularly few solutions, let us check this degree, too.

sage: my_reps(3, 1588)
[[36, 16, 6]]

An other check may be the brute force loop:

sage: for a, b, c in cartesian_product([[1..40], [1..40], [1..40]]):
....:     if a >= b and b >= c and a^2 + b^2 + c^2 == 1001:
....:         print(a,b,c)
....: 
24 16 13
24 19 8
24 20 5
26 15 10
26 17 6
26 18 1
27 16 4
29 12 4
30 10 1
31 6 2

Again, ten solutions. With a bigger $k$ and a bigger $n$.

sage: sols = my_reps(3, 1234567)
sage: len(sols)
0
sage: sols = my_reps(4, 1234567)
len(sols) 
sage: len(sols)
26234

And after a long time, the list is finally built... The recursive lazy implementation above is of course not optimal. Instead, explicitly implementing an iterator that translates the recursion may be optimal. The number of solution grows quickly with $k$. I do not see any reason for having them stored in memory. Why do you want all those $26234$ solutions for $k=4$ and $n=26234$?

If counting is the purpose (well, avoiding zero, restriction to positive and reject of negatives, restriction to in- or decreasing solutions - all this breaks the arithmetics when counting, but ok...) then a generating series may be the better way.

edit flag offensive delete link more
1

answered 2023-04-09 21:36:32 +0200

tmonteil gravatar image

updated 2023-04-09 21:37:29 +0200

Partial answer : for a single representation, you can use the sum_of_k_squares function :

sage: sum_of_k_squares(3, 12345678)
(63, 147, 3510)
sage: 63^2 + 147^2 + 3510^2
12345678

Finding all representations looks costly. Could you please provide an example of B and n you are willing to deal with ?

edit flag offensive delete link more

Comments

Thanks. Numbers less than 100 000 would be ok for a start and small n such as n<=8. The command sum_of_k_squares(4, 2016) shows that it also includes $x_i=0$. Is there a way to fix this?

klaaa gravatar imageklaaa ( 2023-04-09 21:46:45 +0200 )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

1 follower

Stats

Asked: 2023-04-09 20:01:04 +0200

Seen: 652 times

Last updated: Apr 11