Processing math: 100%

First time here? Check out the FAQ!

Ask Your Question
0

Finding representation of a given number as a sum of squares

asked 1 year ago

klaaa gravatar image

updated 1 year ago

Is there a quick method using Sage to find all representations of a given (possible very large) positive integer B as a sum x21+....+x2n of exactly n integers xi>0, where we can assume that xixi+1 for all i?

Thanks for any help

Preview: (hide)

2 Answers

Sort by » oldest newest most voted
1

answered 1 year ago

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 fk, where f=q+q4+q9+q16+q25+ - but order of components matter. We have collected only decreasing solutions. Fortunately, all displayed solutions have three different components, so we expect a 103! count in f3 as coefficient of q1001:

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.

Preview: (hide)
link
1

answered 1 year ago

tmonteil gravatar image

updated 1 year ago

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 ?

Preview: (hide)
link

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 xi=0. Is there a way to fix this?

klaaa gravatar imageklaaa ( 1 year ago )

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: 1 year ago

Seen: 936 times

Last updated: Apr 11 '23