Opened 7 years ago

Closed 7 years ago

#16308 closed enhancement (fixed)

Improve sums of squares

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

Description (last modified by vdelecroix)

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 7 years ago by pbruin

  • Component changed from basic arithmetic to number theory
  • Description modified (diff)

comment:2 follow-up: Changed 7 years ago by jdemeyer

  • Description modified (diff)

comment:3 Changed 7 years ago by ncohen

  • Cc ncohen added

comment:4 Changed 7 years ago by jdemeyer

  • Description modified (diff)

comment:5 in reply to: ↑ 2 Changed 7 years ago by pbruin

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 follow-up: Changed 7 years ago by leif

Define "number".

comment:7 in reply to: ↑ 6 Changed 7 years ago by jdemeyer

Replying to leif:

Define "number".

Element of ZZ.

comment:8 Changed 7 years ago by jdemeyer

  • Branch set to u/jdemeyer/ticket/16308
  • Created changed from 05/08/14 09:28:38 to 05/08/14 09:28:38
  • Modified changed from 05/08/14 15:59:34 to 05/08/14 15:59:34

comment:9 Changed 7 years ago by jdemeyer

  • Commit set to 6b94a6576031e6e341143a4fc1fd83263a800ed9
  • Status changed from new to needs_review

New commits:

6b94a65Improve sums of squares of integers

comment:10 Changed 7 years ago by jdemeyer

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 7 years ago by jdemeyer

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 follow-up: Changed 7 years ago by 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() and spend some effort on #1145. (And possibly add support for algorithm=gap in Integer.factor() if it is worthwhile?)

comment:13 follow-up: Changed 7 years ago by leif

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 ; follow-up: Changed 7 years ago by leif

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 7 years ago by leif

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

comment:16 follow-up: Changed 7 years ago by kcrisman

  • Cc kcrisman 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 7 years ago by leif

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 ; follow-ups: Changed 7 years ago by 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.

comment:19 in reply to: ↑ 12 Changed 7 years ago by jdemeyer

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 ; follow-up: Changed 7 years ago by leif

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 ; follow-up: Changed 7 years ago by jdemeyer

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 7 years ago by pbruin

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 7 years ago by leif

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 7 years ago by vdelecroix

  • Cc vdelecroix added

comment:25 follow-up: Changed 7 years ago by vdelecroix

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 7 years ago by vdelecroix

(and the Cython version can even be improved)

comment:27 Changed 7 years ago by kcrisman

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 7 years ago by leif

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 7 years ago by leif

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 7 years ago by leif

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

comment:31 Changed 7 years ago by vdelecroix

  • Authors changed from Jeroen Demeyer to Jeroen Demeyer, Vincent Delecroix
  • Branch changed from u/jdemeyer/ticket/16308 to u/vdelecroix/16308
  • Commit changed from 6b94a6576031e6e341143a4fc1fd83263a800ed9 to 9b81872a5ff12ba45acbbc32a711c85f2fd58d08

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 7 years ago by vdelecroix

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

Vincent

comment:33 follow-up: Changed 7 years ago by vdelecroix

  • Authors changed from Jeroen Demeyer, Vincent Delecroix to Jeroen Demeyer
  • Branch changed from u/vdelecroix/16308 to u/jeroen/ticket/16308
  • Commit 9b81872a5ff12ba45acbbc32a711c85f2fd58d08 deleted

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 7 years ago by vdelecroix

  • Branch changed from u/jeroen/ticket/16308 to u/jdemeyer/ticket/16308
  • Commit set to 6b94a6576031e6e341143a4fc1fd83263a800ed9

New commits:

6b94a65Improve sums of squares of integers

comment:35 Changed 7 years ago by leif

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 follow-up: Changed 7 years ago by leif

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

comment:37 in reply to: ↑ 33 ; follow-up: Changed 7 years ago by 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.

comment:38 Changed 7 years ago by vdelecroix

  • Description modified (diff)

comment:39 in reply to: ↑ 37 Changed 7 years ago by vdelecroix

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 ; follow-up: Changed 7 years ago by 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?

comment:41 in reply to: ↑ 40 Changed 7 years ago by leif

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 (sum of one square), which is correct.

Sorry for the noise.

Last edited 7 years ago by leif (previous) (diff)

comment:42 Changed 7 years ago by vdelecroix

  • Status changed from needs_review to needs_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 7 years ago by git

  • Commit changed from 6b94a6576031e6e341143a4fc1fd83263a800ed9 to 7b33df2464a78e70766fbb031af1750e5ed6239a

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 7 years ago by jdemeyer

  • Status changed from needs_info to needs_review

comment:45 Changed 7 years ago by vdelecroix

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 7 years ago by vdelecroix

  • Reviewers set to Peter Bruin, Leif Leonhardy , Vincent Delecroix

comment:47 Changed 7 years ago by jdemeyer

  • Branch changed from u/jdemeyer/ticket/16308 to u/vdelecroix/16308
  • Commit changed from 7b33df2464a78e70766fbb031af1750e5ed6239a to f071a8db1316a7bf8503f6c6c9227fcaa6b953f7

New commits:

f071a8dsum of squares: micro improvement + random tests

comment:48 Changed 7 years ago by jdemeyer

  • Status changed from needs_review to positive_review

comment:49 Changed 7 years ago by vbraun

  • Branch changed from u/vdelecroix/16308 to f071a8db1316a7bf8503f6c6c9227fcaa6b953f7
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.