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:  sage6.3 
Component:  number theory  Keywords:  
Cc:  Nathann Cohen, KarlDieter 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: 
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 8 years ago by
Component:  basic arithmetic → number theory 

Description:  modified (diff) 
comment:2 followup: 5 Changed 8 years ago by
Description:  modified (diff) 

comment:3 Changed 8 years ago by
Cc:  Nathann Cohen added 

comment:4 Changed 8 years ago by
Description:  modified (diff) 

comment:5 Changed 8 years ago by
comment:8 Changed 8 years ago by
Branch:  → u/jdemeyer/ticket/16308 

Created:  May 8, 2014, 9:28:38 AM → May 8, 2014, 9:28:38 AM 
Modified:  May 8, 2014, 3:59:34 PM → May 8, 2014, 3:59:34 PM 
comment:9 Changed 8 years ago by
Commit:  → 6b94a6576031e6e341143a4fc1fd83263a800ed9 

Status:  new → needs_review 
New commits:
6b94a65  Improve sums of squares of integers

comment:10 Changed 8 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 8 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 8 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 8 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 followup: 18 Changed 8 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 8 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 8 years ago by
Cc:  KarlDieter 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 Changed 8 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 followups: 20 22 Changed 8 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 Changed 8 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 followup: 21 Changed 8 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 followup: 23 Changed 8 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 Changed 8 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 Changed 8 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 8 years ago by
Cc:  Vincent Delecroix added 

comment:25 followup: 28 Changed 8 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:27 Changed 8 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 Changed 8 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 8 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 8 years ago by
And feel free to implement a different function that returns a list of [all] solutions...
comment:31 Changed 8 years ago by
Authors:  Jeroen Demeyer → Jeroen Demeyer, Vincent Delecroix 

Branch:  u/jdemeyer/ticket/16308 → u/vdelecroix/16308 
Commit:  6b94a6576031e6e341143a4fc1fd83263a800ed9 → 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 8 years ago by
Of course, we can revert to the jeroen branch if my contribution seems useless...
Vincent
comment:33 followup: 37 Changed 8 years ago by
Authors:  Jeroen Demeyer, Vincent Delecroix → Jeroen Demeyer 

Branch:  u/vdelecroix/16308 → u/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
Branch:  u/jeroen/ticket/16308 → u/jdemeyer/ticket/16308 

Commit:  → 6b94a6576031e6e341143a4fc1fd83263a800ed9 
New commits:
6b94a65  Improve sums of squares of integers

comment:35 Changed 8 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 8 years ago by
P.S.: If a solution exists, k_squares()
should IMHO always return a ktuple.
comment:37 followup: 39 Changed 8 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 8 years ago by
Description:  modified (diff) 

comment:39 Changed 8 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 followup: 41 Changed 8 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 Changed 8 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 (some of one square), which is correct.
Sorry for the noise.
comment:42 Changed 8 years ago by
Status:  needs_review → 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 8 years ago by
Commit:  6b94a6576031e6e341143a4fc1fd83263a800ed9 → 7b33df2464a78e70766fbb031af1750e5ed6239a 

comment:44 Changed 8 years ago by
Status:  needs_info → needs_review 

comment:45 Changed 8 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 8 years ago by
Reviewers:  → Peter Bruin, Leif Leonhardy , Vincent Delecroix 

comment:47 Changed 8 years ago by
Branch:  u/jdemeyer/ticket/16308 → u/vdelecroix/16308 

Commit:  7b33df2464a78e70766fbb031af1750e5ed6239a → f071a8db1316a7bf8503f6c6c9227fcaa6b953f7 
New commits:
f071a8d  sum of squares: micro improvement + random tests

comment:48 Changed 8 years ago by
Status:  needs_review → positive_review 

comment:49 Changed 8 years ago by
Branch:  u/vdelecroix/16308 → f071a8db1316a7bf8503f6c6c9227fcaa6b953f7 

Resolution:  → fixed 
Status:  positive_review → 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?