Opened 8 years ago

Closed 8 years ago

#16308 closed enhancement (fixed)

Improve sums of squares

Reported by: Jeroen Demeyer Owned by:
Priority: major Milestone: sage-6.3
Component: number theory Keywords:
Cc: Nathann Cohen, Karl-Dieter Crisman, Vincent Delecroix Merged in:
Authors: Jeroen Demeyer Reviewers: Peter Bruin, Leif Leonhardy , Vincent Delecroix
Report Upstream: N/A Work issues:
Branch: f071a8d (Commits, GitHub, GitLab) Commit: f071a8db1316a7bf8503f6c6c9227fcaa6b953f7
Dependencies: Stopgaps:

Status badges

Description (last modified by Vincent Delecroix)

This ticket will:

  • add a Sage implementation of writing a number as a sum of 2 squares
  • add a Sage implementation of writing a number as a sum of 3 squares
  • add a Sage implementation of writing a number as a sum of 4 squares
  • add a Sage implementation of writing a number as a sum of k squares for any k

follow up: #16374

Change History (49)

comment:1 Changed 8 years ago by Peter Bruin

Component: basic arithmeticnumber theory
Description: modified (diff)

comment:2 Changed 8 years ago by Jeroen Demeyer

Description: modified (diff)

comment:3 Changed 8 years ago by Nathann Cohen

Cc: Nathann Cohen added

comment:4 Changed 8 years ago by Jeroen Demeyer

Description: modified (diff)

comment:5 in reply to:  2 Changed 8 years ago by Peter Bruin

Replying to jdemeyer:

OK, I should have written a comment to ask what the idea was instead of adding possible interpretations to the description. Given that your name is in the Author field, are you working on this?

comment:6 Changed 8 years ago by Leif Leonhardy

Define "number".

comment:7 in reply to:  6 Changed 8 years ago by Jeroen Demeyer

Replying to leif:

Define "number".

Element of ZZ.

comment:8 Changed 8 years ago by Jeroen Demeyer

Branch: u/jdemeyer/ticket/16308
Created: May 8, 2014, 9:28:38 AMMay 8, 2014, 9:28:38 AM
Modified: May 8, 2014, 3:59:34 PMMay 8, 2014, 3:59:34 PM

comment:9 Changed 8 years ago by Jeroen Demeyer

Commit: 6b94a6576031e6e341143a4fc1fd83263a800ed9
Status: newneeds_review

New commits:

6b94a65Improve sums of squares of integers

comment:10 Changed 8 years ago by Jeroen Demeyer

It is embarrassing to say so, but for many examples, the old GAP code is actually faster. The reason seems to be that GAP beats PARI on factoring small numbers (which is a big surprise to me!). But GAP doesn't support large numbers:

sage: two_squares(3^100+1)
RuntimeError: Gap produced error output
Error, sorry,  cannot factor [ 1109667565348268894775749972955601 ]
type 'return;' to try again with a larger number of trials in
FactorsRho (or use option 'RhoTrials')


   executing TwoSquares(515377520732011331036461129765621272702107522002);

comment:11 Changed 8 years ago by Jeroen Demeyer

GAP is amazing:

sage: timeit('libgap(19180172397815991981).FactorsInt()', number=1000, repeat=10)
1000 loops, best of 10: 103 µs per loop
sage: timeit('pari(19180172397815991981).factor(proof=False)', number=1000, repeat=10)
1000 loops, best of 10: 2.43 ms per loop

comment:12 Changed 8 years ago by Peter Bruin

Interesting... Since this code is probably not the place where Sage should have to make a decision about which factoring algorithm to use, it looks tempting to use Integer.factor() instead of pari.factor() and spend some effort on #1145. (And possibly add support for algorithm=gap in Integer.factor() if it is worthwhile?)

comment:13 Changed 8 years ago by Leif Leonhardy

And I'd still support the algorithm keyword...

Does make returning (a,b) with a (or b) zero any sense?

comment:14 in reply to:  13 ; Changed 8 years ago by Leif Leonhardy

Replying to leif:

Does make returning (a,b) with a (or b) zero any sense?

To be more precise, IMHO both a and b should be < n1/2 (for two_squares()).

I don't think anybody would expect e.g. two_squares(25) to return the trivial solution (0,5) instead of (3,4).

comment:15 Changed 8 years ago by Leif Leonhardy

(GAP apparently only returns (0,sqrt(n)) if there is no other solution, with sqrt(n) integer of course.)

comment:16 Changed 8 years ago by Karl-Dieter Crisman

Cc: Karl-Dieter Crisman added

I don't think anybody would expect e.g. two_squares(25) to return the trivial solution (0,5) instead of (3,4).

Hmm, well that IS a solution. And if you want to use it to approximate the value of pi, you need all of them (though this one in any case wouldn't give negative solutions).

comment:17 in reply to:  16 Changed 8 years ago by Leif Leonhardy

Replying to kcrisman:

I don't think anybody would expect e.g. two_squares(25) to return the trivial solution (0,5) instead of (3,4).

Hmm, well that IS a solution.

Depends on the definition of two_squares(); the previous one in Sage was pretty vague.

comment:18 in reply to:  14 ; Changed 8 years ago by Jeroen Demeyer

Replying to leif:

I don't think anybody would expect e.g. two_squares(25) to return the trivial solution (0,5) instead of (3,4).

From a number-theory point of view, there is nothing wrong with (0,5). In fact, in some cases, the only solution is the "trivial" one, for example for 9.

comment:19 in reply to:  12 Changed 8 years ago by Jeroen Demeyer

Replying to pbruin:

Interesting... Since this code is probably not the place where Sage should have to make a decision about which factoring algorithm to use, it looks tempting to use Integer.factor() instead of pari.factor()

I guess you're right, but I decided to do everything in PARI for speed, also with the idea that the code could be changed to use direct calls to the PARI library if needed.

comment:20 in reply to:  18 ; Changed 8 years ago by Leif Leonhardy

Replying to jdemeyer:

Replying to leif:

I don't think anybody would expect e.g. two_squares(25) to return the trivial solution (0,5) instead of (3,4).

From a number-theory point of view, there is nothing wrong with (0,5). In fact, in some cases, the only solution is the "trivial" one, for example for 9.

As mentioned, depends on how you specify two_squares().

From a user perspective, I think there should be some optional argument no_trivial_solutions, as feeding the function with perfect squares otherwise doesn't make much sense.

comment:21 in reply to:  20 ; Changed 8 years ago by Jeroen Demeyer

Replying to leif:

as feeding the function with perfect squares otherwise doesn't make much sense.

Depends on your application...

Anyway, if people really care about non-trivial solutions, let those people implement such a flag, I'm not going to do it.

comment:22 in reply to:  18 Changed 8 years ago by Peter Bruin

Replying to jdemeyer:

Replying to leif:

I don't think anybody would expect e.g. two_squares(25) to return the trivial solution (0,5) instead of (3,4).

From a number-theory point of view, there is nothing wrong with (0,5). In fact, in some cases, the only solution is the "trivial" one, for example for 9.

I agree. Given that 0 is a square, it makes complete sense to allow it as a term in a sum of squares. That convention makes a lot of things nicer, for example:

  • every square is also a sum of two squares, and more generally every sum of k squares is also a sum of k + 1 squares;
  • n is a sum of 2 squares if and only if n can be factored non-trivially in Z[i];
  • n is a sum of k squares if and only if the lattice Zk has a point of squared distance n from the origin.

There are also nice formulae for the number of ways to write n as a sum of k squares if you interpret this number as the number of points (x1 , ... , xk) in Zk with x12 + ... + xk2 = n.

comment:23 in reply to:  21 Changed 8 years ago by Leif Leonhardy

Replying to jdemeyer:

Anyway, if people really care about non-trivial solutions, let those people implement such a flag, I'm not going to do it.

algorithm="gap"?

We'd need that for backwards compatibility anyway, besides that you said GAP is faster (in factoring at least) for some (which kind of?) inputs.

comment:24 Changed 8 years ago by Vincent Delecroix

Cc: Vincent Delecroix added

comment:25 Changed 8 years ago by Vincent Delecroix

Hi there,

For very small integers (<500) a naive cython implementation is much faster (at least x100)... as this was the initial purpose of Nathann I make the remark here.

def two_squares_pyx(int n):
    cdef int i,ii,j,nn
    if n < 2:
        return True

    i = ii = 0
    while ii < n:
        j = 0
        while j <= i:
            nn = ii + j*j
            if nn >= n:
                break
            j += 1

        if nn == n:
            return (i,j)

        i += 1
        ii = i*i

    raise ValueError("%d is not the sum of 2 squares"%n)

then

sage: %timeit two_squares(37)
10000 loops, best of 3: 108 µs per loop
sage: %timeit two_squares_pyx(37)
1000000 loops, best of 3: 464 ns per loop

sage: %timeit two_squares(97)
10000 loops, best of 3: 185 µs per loop
sage: %timeit two_squares_pyx(97)
1000000 loops, best of 3: 558 ns per loop

sage: %timeit two_squares(212)
10000 loops, best of 3: 130 µs per loop
sage: %timeit two_squares_pyx(212)
1000000 loops, best of 3: 734 ns per loop

sage: %timeit two_squares(500)
10000 loops, best of 3: 131 µs per loop
sage: %timeit two_squares_pyx(500)
1000000 loops, best of 3: 1 µs per loop

What do you think of using it for entries < 500?

Vincent

comment:26 Changed 8 years ago by Vincent Delecroix

(and the Cython version can even be improved)

comment:27 Changed 8 years ago by Karl-Dieter Crisman

I was hoping someone would do this.

I agree with pbruin, though of course this is not the same function as "number of ways to write as two squares" or "all ways to write as sum of two squares". This is what happens when we "unwrap" :-)

comment:28 in reply to:  25 Changed 8 years ago by Leif Leonhardy

Replying to vdelecroix:

def two_squares_pyx(int n):
    cdef int i,ii,j,nn
    if n < 2:
        return True

    i = ii = 0
    while ii < n:
        j = 0
        while j <= i:
            nn = ii + j*j
            if nn >= n:
                break
            j += 1

        if nn == n:
            return (i,j)

        i += 1
        ii = i*i

    raise ValueError("%d is not the sum of 2 squares"%n)

Minor correction: while ii <= n (such that we get "trivial" solutions with j==0 as well). Also, I'd swap i and j, or return (j,i).

comment:29 Changed 8 years ago by Leif Leonhardy

Well, I don't care much about a==0 (a matter of definition, it's ok, and one usually has to treat zeroes differently anyway), but normally one isn't interested in a trivial solution, unless it's the only one, and if there's no way to get (or even get to know whether there is) a non-trivial solution [as well]...

comment:30 Changed 8 years ago by Leif Leonhardy

And feel free to implement a different function that returns a list of [all] solutions...

comment:31 Changed 8 years ago by Vincent Delecroix

Authors: Jeroen DemeyerJeroen Demeyer, Vincent Delecroix
Branch: u/jdemeyer/ticket/16308u/vdelecroix/16308
Commit: 6b94a6576031e6e341143a4fc1fd83263a800ed99b81872a5ff12ba45acbbc32a711c85f2fd58d08

Hello,

I implemented what I said in Cython (i.e. two_squares_pyx, three_squares_pyx and four_squares_pyx). Timings are now much better for small values (but it also influences small values).

Based on the timings below (with the old version), I put the threshold to 500 for the three functions. But if you have strong opinion for another threshold please tell me.

There are many comments, but is there somebody up for a review ?

Vincent

Two squares

sage: %timeit two_squares(101)
10000 loops, best of 3: 94.5 µs per loop
sage: %timeit two_squares_pyx(101)
1000000 loops, best of 3: 780 ns per loop

sage: timeit("try:\n    two_squares(102)\nexcept ValueError:\n    pass")
625 loops, best of 3: 57.4 µs per loop
sage: timeit("try:\n    two_squares_pyx(102)\nexcept ValueError:\n    pass")
625 loops, best of 3: 6.4 µs per loop

sage: %timeit two_squares(200)
10000 loops, best of 3: 87.3 µs per loop
sage: %timeit two_squares_pyx(200)
1000000 loops, best of 3: 788 ns per loop

sage: timeit("try:\n    two_squares(201)\nexcept ValueError:\n    pass")
625 loops, best of 3: 39.6 µs per loop
sage: timeit("try:\n    two_squares_pyx(201)\nexcept ValueError:\n    pass")
625 loops, best of 3: 6.61 µs per loop

sage: %timeit two_squares(305)
10000 loops, best of 3: 158 µs per loop
sage: %timeit two_squares_pyx(305)
1000000 loops, best of 3: 998 ns per loop

sage: timeit("try:\n    two_squares(304)\nexcept ValueError:\n    pass")
625 loops, best of 3: 50.8 µs per loop
sage: timeit("try:\n    two_squares_pyx(304)\nexcept ValueError:\n    pass")
625 loops, best of 3: 6.71 µs per loop

sage: %timeit two_squares(400)
10000 loops, best of 3: 69.4 µs per loop
sage: %timeit two_squares_pyx(400)
1000000 loops, best of 3: 1.01 µs per loop

sage: timeit("try:\n    two_squares(402)\nexcept ValueError:\n    pass")
625 loops, best of 3: 57.5 µs per loop
sage: timeit("try:\n    two_squares_pyx(402)\nexcept ValueError:\n    pass")
625 loops, best of 3: 6.85 µs per loop

sage: %timeit two_squares(500)
10000 loops, best of 3: 119 µs per loop
sage: %timeit two_squares_pyx(500)
1000000 loops, best of 3: 1.18 µs per loop

sage: timeit("try:\n    two_squares(501)\nexcept ValueError:\n    pass")
625 loops, best of 3: 39.5 µs per loop
sage: timeit("try:\n    two_squares_pyx(501)\nexcept ValueError:\n    pass")
625 loops, best of 3: 6.76 µs per loop

Three squares

sage: %timeit three_squares(106)
1000 loops, best of 3: 266 µs per loop
sage: %timeit three_squares_pyx(106)
1000000 loops, best of 3: 1.38 µs per loop

sage: timeit("try:\n    three_squares(103)\nexcept ValueError:\n    pass")
625 loops, best of 3: 52.2 µs per loop
sage: timeit("try:\n    three_squares_pyx(103)\nexcept ValueError:\n    pass")
625 loops, best of 3: 6.9 µs per loop

sage: %timeit three_squares(200)
1000 loops, best of 3: 205 µs per loop
sage: %timeit three_squares_pyx(200)
1000000 loops, best of 3: 1.81 µs per loop

sage: timeit("try:\n    three_squares(303)\nexcept ValueError:\n    pass")
625 loops, best of 3: 52.1 µs per loop
sage: timeit("try:\n    three_squares_pyx(303)\nexcept ValueError:\n    pass")
625 loops, best of 3: 9.24 µs per loop

sage: %timeit three_squares(406)
1000 loops, best of 3: 189 µs per loop
sage: %timeit three_squares_pyx(406)
100000 loops, best of 3: 3.46 µs per loop

sage: timeit("try:\n    three_squares(407)\nexcept ValueError:\n    pass")
625 loops, best of 3: 51.6 µs per loop
sage: timeit("try:\n    three_squares_pyx(407)\nexcept ValueError:\n    pass")
625 loops, best of 3: 10.5 µs per loop

sage: %timeit three_squares(500)
10000 loops, best of 3: 167 µs per loop
sage: %timeit three_squares_pyx(500)
100000 loops, best of 3: 3.89 µs per loop

sage: timeit("try:\n    three_squares(503)\nexcept ValueError:\n    pass")
625 loops, best of 3: 51.2 µs per loop
sage: timeit("try:\n    three_squares_pyx(503)\nexcept ValueError:\n    pass")
625 loops, best of 3: 10.7 µs per loop

Four squares

sage: %timeit four_squares(503)
1000 loops, best of 3: 259 µs per loop
sage: %timeit four_squares_pyx(503)
100000 loops, best of 3: 13.5 µs per loop

New commits:

616c606trac #16308: cython version of x_squares
9b81872trac #16308: add the file arith_pyx.pyx

comment:32 Changed 8 years ago by Vincent Delecroix

Of course, we can revert to the jeroen branch if my contribution seems useless...

Vincent

comment:33 Changed 8 years ago by Vincent Delecroix

Authors: Jeroen Demeyer, Vincent DelecroixJeroen Demeyer
Branch: u/vdelecroix/16308u/jeroen/ticket/16308
Commit: 9b81872a5ff12ba45acbbc32a711c85f2fd58d08

Hi,

Sorry. I was a bit fast. The decision of choosing my branch or not might be the choice of the author of the ticket. So if you want to have a look at my changes, see either u/vdelecroix/16308 or public/16308

Vincent

comment:34 Changed 8 years ago by Vincent Delecroix

Branch: u/jeroen/ticket/16308u/jdemeyer/ticket/16308
Commit: 6b94a6576031e6e341143a4fc1fd83263a800ed9

New commits:

6b94a65Improve sums of squares of integers

comment:35 Changed 8 years ago by Leif Leonhardy

There's even more room for optimizations; a C/Cython version is probably faster up to at least 10000 or even 216 (for two_squares()).

I'd also compare the overall speed over a range.

comment:36 Changed 8 years ago by Leif Leonhardy

P.S.: If a solution exists, k_squares() should IMHO always return a k-tuple.

comment:37 in reply to:  33 ; Changed 8 years ago by Leif Leonhardy

Replying to vdelecroix:

Sorry. I was a bit fast. The decision of choosing my branch or not might be the choice of the author of the ticket. So if you want to have a look at my changes, see either u/vdelecroix/16308 or public/16308

Probably open a separate ticket "Improve sums of squares for small input" and we can later (re)base one of them on the other.

comment:38 Changed 8 years ago by Vincent Delecroix

Description: modified (diff)

comment:39 in reply to:  37 Changed 8 years ago by Vincent Delecroix

Replying to leif:

Replying to vdelecroix:

Sorry. I was a bit fast. The decision of choosing my branch or not might be the choice of the author of the ticket. So if you want to have a look at my changes, see either u/vdelecroix/16308 or public/16308

Probably open a separate ticket "Improve sums of squares for small input" and we can later (re)base one of them on the other.

Good idea. Done in #16374.

comment:40 in reply to:  36 ; Changed 8 years ago by Jeroen Demeyer

Replying to leif:

P.S.: If a solution exists, k_squares() should IMHO always return a k-tuple.

Sure, in which case does it not do that?

comment:41 in reply to:  40 Changed 8 years ago by Leif Leonhardy

Replying to jdemeyer:

Replying to leif:

P.S.: If a solution exists, k_squares() should IMHO always return a k-tuple.

Sure, in which case does it not do that?

Sorry, I had been looking at some diff of Vincent's code against yours, and only saw some return (foo,) (and I swear it was more than one!1!!1 ;-) ), but in your original code, or the one that's again here now, it's of course only in the case k=1 (some of one square), which is correct.

Sorry for the noise.

Version 0, edited 8 years ago by Leif Leonhardy (next)

comment:42 Changed 8 years ago by Vincent Delecroix

Status: needs_reviewneeds_info

Hi,

Related to comment:12, can you justify why using pari? We already have factor and is_pseudoprime and sqrtrem might be the replacement for N.sqrtint.

There are also trailing whitespaces on the lines 5161 and 5246.

Vincent

comment:43 Changed 8 years ago by git

Commit: 6b94a6576031e6e341143a4fc1fd83263a800ed97b33df2464a78e70766fbb031af1750e5ed6239a

Branch pushed to git repo; I updated commit sha1. New commits:

c4f603bUse Integer.factor() for factoring
f6754d0Merge remote-tracking branch 'origin/develop' into ticket/16308
7b33df2Speed up sums of squares

comment:44 Changed 8 years ago by Jeroen Demeyer

Status: needs_infoneeds_review

comment:45 Changed 8 years ago by Vincent Delecroix

Hey,

Looks cool! I added micro optimizations and random tests at u/vdelecroix/16308. If you like my changes, this is good for a positive review!

Vincent

comment:46 Changed 8 years ago by Vincent Delecroix

Reviewers: Peter Bruin, Leif Leonhardy , Vincent Delecroix

comment:47 Changed 8 years ago by Jeroen Demeyer

Branch: u/jdemeyer/ticket/16308u/vdelecroix/16308
Commit: 7b33df2464a78e70766fbb031af1750e5ed6239af071a8db1316a7bf8503f6c6c9227fcaa6b953f7

New commits:

f071a8dsum of squares: micro improvement + random tests

comment:48 Changed 8 years ago by Jeroen Demeyer

Status: needs_reviewpositive_review

comment:49 Changed 8 years ago by Volker Braun

Branch: u/vdelecroix/16308f071a8db1316a7bf8503f6c6c9227fcaa6b953f7
Resolution: fixed
Status: positive_reviewclosed
Note: See TracTickets for help on using tickets.