Opened 7 years ago
Closed 7 years ago
#16308 closed enhancement (fixed)
Improve sums of squares
Reported by:  jdemeyer  Owned by:  

Priority:  major  Milestone:  sage6.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 )
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 anyk
follow up: #16374
Change History (49)
comment:1 Changed 7 years ago by
 Component changed from basic arithmetic to number theory
 Description modified (diff)
comment:2 followup: ↓ 5 Changed 7 years ago by
 Description modified (diff)
comment:3 Changed 7 years ago by
 Cc ncohen added
comment:4 Changed 7 years ago by
 Description modified (diff)
comment:5 in reply to: ↑ 2 Changed 7 years ago by
comment:6 followup: ↓ 7 Changed 7 years ago by
Define "number".
comment:7 in reply to: ↑ 6 Changed 7 years ago by
comment:8 Changed 7 years ago by
 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
 Commit set to 6b94a6576031e6e341143a4fc1fd83263a800ed9
 Status changed from new to needs_review
New commits:
6b94a65  Improve sums of squares of integers

comment:10 Changed 7 years ago by
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
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 followup: ↓ 19 Changed 7 years ago by
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 followup: ↓ 14 Changed 7 years ago by
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 ; followup: ↓ 18 Changed 7 years ago by
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 < n^{1/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
(GAP apparently only returns (0,sqrt(n)) if there is no other solution, with sqrt(n) integer of course.)
comment:16 followup: ↓ 17 Changed 7 years ago by
 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
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 ; followups: ↓ 20 ↓ 22 Changed 7 years ago by
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 numbertheory 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
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 ofpari.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 ; followup: ↓ 21 Changed 7 years ago by
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 numbertheory 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 ; followup: ↓ 23 Changed 7 years ago by
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 nontrivial 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
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 numbertheory 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 nontrivially in Z[i];
 n is a sum of k squares if and only if the lattice Z^{k} 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 (x_{1} , ... , x_{k}) in Z^{k} with x_{1}^{2} + ... + x_{k}^{2} = n.
comment:23 in reply to: ↑ 21 Changed 7 years ago by
Replying to jdemeyer:
Anyway, if people really care about nontrivial 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
 Cc vdelecroix added
comment:25 followup: ↓ 28 Changed 7 years ago by
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
(and the Cython version can even be improved)
comment:27 Changed 7 years ago by
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
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
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 nontrivial solution [as well]...
comment:30 Changed 7 years ago by
And feel free to implement a different function that returns a list of [all] solutions...
comment:31 Changed 7 years ago by
 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:
616c606  trac #16308: cython version of x_squares

9b81872  trac #16308: add the file arith_pyx.pyx

comment:32 Changed 7 years ago by
Of course, we can revert to the jeroen branch if my contribution seems useless...
Vincent
comment:33 followup: ↓ 37 Changed 7 years ago by
 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
 Branch changed from u/jeroen/ticket/16308 to u/jdemeyer/ticket/16308
 Commit set to 6b94a6576031e6e341143a4fc1fd83263a800ed9
New commits:
6b94a65  Improve sums of squares of integers

comment:35 Changed 7 years ago by
There's even more room for optimizations; a C/Cython version is probably faster up to at least 10000 or even 2^{16} (for two_squares()
).
I'd also compare the overall speed over a range.
comment:36 followup: ↓ 40 Changed 7 years ago by
P.S.: If a solution exists, k_squares()
should IMHO always return a ktuple.
comment:37 in reply to: ↑ 33 ; followup: ↓ 39 Changed 7 years ago by
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
 Description modified (diff)
comment:39 in reply to: ↑ 37 Changed 7 years ago by
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 ; followup: ↓ 41 Changed 7 years ago by
Replying to leif:
P.S.: If a solution exists, k
_squares()
should IMHO always return a ktuple.
Sure, in which case does it not do that?
comment:41 in reply to: ↑ 40 Changed 7 years ago by
Replying to jdemeyer:
Replying to leif:
P.S.: If a solution exists, k
_squares()
should IMHO always return a ktuple.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.
comment:42 Changed 7 years ago by
 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
 Commit changed from 6b94a6576031e6e341143a4fc1fd83263a800ed9 to 7b33df2464a78e70766fbb031af1750e5ed6239a
comment:44 Changed 7 years ago by
 Status changed from needs_info to needs_review
comment:45 Changed 7 years ago by
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
 Reviewers set to Peter Bruin, Leif Leonhardy , Vincent Delecroix
comment:47 Changed 7 years ago by
 Branch changed from u/jdemeyer/ticket/16308 to u/vdelecroix/16308
 Commit changed from 7b33df2464a78e70766fbb031af1750e5ed6239a to f071a8db1316a7bf8503f6c6c9227fcaa6b953f7
New commits:
f071a8d  sum of squares: micro improvement + random tests

comment:48 Changed 7 years ago by
 Status changed from needs_review to positive_review
comment:49 Changed 7 years ago by
 Branch changed from u/vdelecroix/16308 to f071a8db1316a7bf8503f6c6c9227fcaa6b953f7
 Resolution set to fixed
 Status changed from positive_review to closed
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?