Is there a way to make two_squares faster ? (it uses the interface with GAP)

[Nathann Cohen]

Helloooooooooo everybody !

In the brand new design code that Volker is merging into Sage these days [1], one of the functions calls two_squares very often. two_squares tells you if an integer can be written as the sum of two squares, and because it calls GAP this functions is a bit too slow, and it is actually the bottleneck in some of the computations...

Do you know if :

1) There would be a way to make it quicker to communicate with GAP, as we call it often

2) There is another implementation of two_squares in some other library we ship that we could use instead

3) It can be re-written easily with pure Sage tools ?

Thanks for your help ! :-)

Aaaaaaaaaaaaaaand I hope we will eventually hear some praise for Sage because of its designs code... It is quite cool already :-P

Nathann

[1] sorry for the conflicts Volker, I really try to make it as clean as I can and to fix them quickly, but almost-linearly-ordered tickets is hell to manage :-/

sage: 4367424234**2+325235235**2                  

19180172397815991981

sage: %prun two_squares(19180172397815991981)

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)

        2    0.001    0.000    0.001    0.000 {posix.write}

        6    0.000    0.000    0.001    0.000 pexpect.py:918(expect_list)

      112    0.000    0.000    0.000    0.000 {method 'search' of '_sre.SRE_Pattern' objects}

        1    0.000    0.000    0.001    0.001 gap.py:576(_execute_line)

        1    0.000    0.000    0.000    0.000 {eval}

        1    0.000    0.000    0.001    0.001 arith.py:4910(two_squares)

        1    0.000    0.000    0.001    0.001 expect.py:1154(eval)

       37    0.000    0.000    0.000    0.000 {method 'start' of '_sre.SRE_Match' objects}

        1    0.000    0.000    0.001    0.001 gap.py:510(eval)

        1    0.000    0.000    0.001    0.001 gap.py:664(_eval_line)

        1    0.000    0.000    0.000    0.000 pexpect.py:502(read_nonblocking)

....

[John Cremona]

Try William's course notes here: http://modular.math.washington.edu/edu/2007/spring/ent/ent-html/node75.html where he even gives Sage code of a fast algorithm.

John

--
You received this message because you are subscribed to the Google Groups "sage-devel" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sage-devel+...@googlegroups.com.
To post to this group, send email to sage-...@googlegroups.com.
Visit this group at http://groups.google.com/group/sage-devel.
For more options, visit https://groups.google.com/d/optout.

[John Cremona]

On 7 May 2014 20:55, John Cremona <john.c...@gmail.com> wrote:

Try William's course notes here: http://modular.math.washington.edu/edu/2007/spring/ent/ent-html/node75.html where he even gives Sage code of a fast algorithm.

Sorry, that's only for primes.  I doubt if factoring your number first is the best thing to do.  I expect that lattice basis reduction is.  To write the number as a sum of two rational squares (a different problem, I admit!) has a lot written about it, e.g. by me:  http://www.ams.org/journals/mcom/2003-72-243/S0025-5718-02-01480-1/home.html  but that would also require factorization.

Just some ideas...

John

[John Cremona]

Last word:  this might be efficient and needs only a little work:

sage: n = 19180172397815991981 

sage: R = ZZ[sqrt(-1)] 

sage: R(n).factor() 

(-1) * (-283307*I - 233098) * (283307*I - 233098) * (-3252*I - 2293) * 3^2 * (-3252*I + 2293)  

[Simon King]

Hi Nathann,

On 7 May 2014 20:38, Nathann Cohen <nathan...@gmail.com> wrote:
> Do you know if :
>
> 1) There would be a way to make it quicker to communicate with GAP, as we
> call it often

There is now libgap, and it has less overhead:
sage: libgap(1235^2+2398234^2).TwoSquares().sage()
[888355, 2227634]
sage: _[0]^2+_[1]^2 == 1235^2+2398234^2
True
sage: def Lib2sq(L):
....: for n in L:
....: try:
....: libgap(n).TwoSquares().sage()
....: except ValueError:
....: pass
....:
sage: def Sage2sq(L):
....: for n in L:
....: try:
....: two_squares(n)
....: except ValueError:
....: pass
....:
sage: L = [ZZ.random_element(10^6,10^12) for i in range(1000)]
sage: %timeit Sage2sq(L)
1 loops, best of 3: 742 ms per loop
sage: %timeit Lib2sq(L)
1 loops, best of 3: 416 ms per loop

But this is just about the overhead in talking with GAP, it is the same
underlying algorithm.

Cheers,
Simon

[Vincent Delecroix]

Ho,

There is stil a factor 10 compared to factorisation:

sage:

2014-05-07 22:08 UTC+02:00, Simon King <simon...@uni-jena.de>:

[Vincent Delecroix]

Hi,

Calling factor is 10 times faster

sage: n = 1235^2+2398234^2
sage: %timeit libgap(n).TwoSquares().sage()
1000 loops, best of 3: 680 µs per loop
sage: %timeit n.factor()
10000 loops, best of 3: 42 µs per loop

Moreover, in the context of designs, the numbers are of the size of
1000. So it is worth in that case to factorize the number... I am
pretty sure that given that context, it would be the fastest.

Vincent

[Nathann Cohen]

Hello !

 

There is now libgap, and it has less overhead:

Thanks ! So at the very least Sage's function can be updated, even if the "factor" test is not a good idea in general (for Sage) while it is for our purposes. I will write a patch for that.

Nathann

[John Cremona]

After sending 3 replies yesterday I thought that I should not send any more, but here is a reason why you should not expect to solve this without factorization:  let p, q be primes both 1 mod 4 and n=p*q.    Now n has 2 essentially different ways to be expressed as a sum of two squares, and if you know both of them then you can factor n (write n=a^2+b^2=c^2+d^2 and take the gcd of a+b*i and c+d*i in ZZ[i] using different signs).

The lattice method I mentioned first requires solving the congruence x^2=-1 (mod n) and all known methods to do that require factoring n also.

So there may well be speedups possible compared with calling gap, but don't expect miracles for numbers too large to be factored in reasonable time.  (Another case to think about is n=p*q where p=q=3 (mod 4);  then n is not a sum of 2 squares but n=1 (mod 4) and if you wrongly assumed that n were prime then you would think it possible.)

--

[Jeroen Demeyer]

On 2014-05-07 21:58, John Cremona wrote:
>
> On 7 May 2014 20:55, John Cremona <john.c...@gmail.com

> <mailto:john.c...@gmail.com>> wrote:
>
> Try William's course notes here:
> http://modular.math.washington.edu/edu/2007/spring/ent/ent-html/node75.html
> where he even gives Sage code of a fast algorithm.
>
>
> Sorry, that's only for primes. I doubt if factoring your number first
> is the best thing to do. I expect that lattice basis reduction is. To
> write the number as a sum of two rational squares (a different problem,
> I admit!) has a lot written about it, e.g. by me:
> http://www.ams.org/journals/mcom/2003-72-243/S0025-5718-02-01480-1/home.html
> but that would also require factorization.

If I'm not mistaken, the problem of writing N as a^2 + b^2 is equivalent
to finding a square root of -1 modulo N. I know that computing arbitrary
square roots modulo N is equivalent to factoring, so this seems pretty
close to require factoring.

Jeroen.

[John Cremona]

On 8 May 2014 09:30, Jeroen Demeyer <jdem...@cage.ugent.be> wrote:

 

If I'm not mistaken, the problem of writing N as a^2 + b^2 is equivalent to finding a square root of -1 modulo N. I know that computing arbitrary square roots modulo N is equivalent to factoring, so this seems pretty close to require factoring.

Agreed -- see my post of a couple of minutes ago!

John

[Nathann Cohen]

Hmmmmmmmm guys I have to admit that I do not follow you in your number-theoretic remarks :-P

What I know for sure is that Simon's way of calling gap would make the function faster in all cases.

What I do not know for sure is whether what you all think of would be faster in the general case (besides not knowing how to implement  it).

Thus :

1) If the only way we can improve Sage in the general case is to call GAP differently I can do it.

2) If we shoud rewrite this function differently with Sage tools to make it faster, I cannot do it because I do not understand how.

SOooooooooooooo what should we do ?

Nathann

--
You received this message because you are subscribed to a topic in the Google Groups "sage-devel" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/sage-devel/n6gip5AK4P8/unsubscribe.
To unsubscribe from this group and all its topics, send an email to sage-devel+...@googlegroups.com.

[Jeroen Demeyer]

On 2014-05-08 11:01, Nathann Cohen wrote:
> Hmmmmmmmm guys I have to admit that I do not follow you in your
> number-theoretic remarks :-P
>
> What I know for sure is that Simon's way of calling gap would make the
> function faster in all cases.
>
> What I do not know for sure is whether what you all think of would be
> faster in the general case (besides not knowing how to implement it).
>
> Thus :
>
> 1) If the only way we can improve Sage in the general case is to call
> GAP differently I can do it.
> 2) If we shoud rewrite this function differently with Sage tools to make
> it faster, I cannot do it because I do not understand how.

Is there a ticket already? I would propose a direct implementation in
Sage (possibly using the PARI library).

Jeroen.

[Nathann Cohen]

> Is there a ticket already? I would propose a direct implementation in Sage (possibly using the PARI library).

None.

Nathann

[Jeroen Demeyer]

See
http://trac.sagemath.org/ticket/16308

[Nathann Cohen]

Just added my name in cc:. Thanks ! ;-)

Nathann

On 8 May 2014 11:28, Jeroen Demeyer <jdem...@cage.ugent.be> wrote:
> See
> http://trac.sagemath.org/ticket/16308

[leif]

Vincent Delecroix wrote:
> Moreover, in the context of designs, the numbers are of the size of
> 1000.

Do you mean numbers with 1000 (decimal) digits, or numbers around 10^3
(about 10 bits!)?

In the latter case, presumably even the "naive method", implemented in
Cython/C, would outperform any of the other available, and one could use
a fairly small look-up table as well.


-leif

> So it is worth in that case to factorize the number... I am
> pretty sure that given that context, it would be the fastest.

--
() The ASCII Ribbon Campaign
/\ Help Cure HTML E-Mail

[Nathann Cohen]

Yooooooooo !!

> Do you mean numbers with 1000 (decimal) digits, or numbers around 10^3
> (about 10 bits!)?

Ahahahh. No he means between 0 and 1000. Okay, let's say 5000 if we
add some caching to the functions, but nothing larger than that :-)

> In the latter case, presumably even the "naive method", implemented in
> Cython/C, would outperform any of the other available, and one could use a
> fairly small look-up table as well.

Yepyep, but it would be cool if we didn't have to add anything in our
totally unrelated code. I just hoped that there would be a way, like
libgap, to make this function faster. Right now Jeroen has created a
ticket for that, aaaaaaand so perhaps we could just play it lazy on
the design side and use the general function, even though we could be
a bit faster if we had hmh

[leif]

John Cremona wrote:
> On 8 May 2014 09:30, Jeroen Demeyer <jdem...@cage.ugent.be

> <mailto:jdem...@cage.ugent.be>> wrote:
>
> If I'm not mistaken, the problem of writing N as a^2 + b^2 is
> equivalent to finding a square root of -1 modulo N. I know that
> computing arbitrary square roots modulo N is equivalent to
> factoring, so this seems pretty close to require factoring.
>
>
> Agreed -- see my post of a couple of minutes ago!

Where "pretty close" means?

Obviously, for the two squares problem, one doesn't necessarily have to
fully factor N if no solution exists. And since the relative number of
those for which a solution exists tends to zero, one shouldn't do so in
an efficient implementation... ;-)

[Of course that's more work, perhaps unless one can use lazy factoring,
but since you are ambitious and eager...]

(The currently proposed solution first fully factors positive N. [1])


-leif

[1] http://trac.sagemath.org/ticket/16308

[John Cremona]

On 8 May 2014 20:39, leif <not.r...@online.de> wrote:

John Cremona wrote:

On 8 May 2014 09:30, Jeroen Demeyer <jdem...@cage.ugent.be

<mailto:jdem...@cage.ugent.be>> wrote:

    If I'm not mistaken, the problem of writing N as a^2 + b^2 is
    equivalent to finding a square root of -1 modulo N. I know that
    computing arbitrary square roots modulo N is equivalent to
    factoring, so this seems pretty close to require factoring.


Agreed -- see my post of a couple of minutes ago!

Where "pretty close" means?

I made this rather precise in the case N=p*q in an earlier message.

 

Obviously, for the two squares problem, one doesn't necessarily have to fully factor N if no solution exists.  And since the relative number of those for which a solution exists tends to zero, one shouldn't do so in an efficient implementation... ;-)

I don't think it is possible to tell whether or not a solution exists without factoring N (except using the stupid method of course).

 

[Of course that's more work, perhaps unless one can use lazy factoring, but since you are ambitious and eager...]

(The currently proposed solution first fully factors positive N. [1])


-leif

[1] http://trac.sagemath.org/ticket/16308

--
() The ASCII Ribbon Campaign
/\   Help Cure HTML E-Mail

--
You received this message because you are subscribed to the Google Groups "sage-devel" group.

To unsubscribe from this group and stop receiving emails from it, send an email to sage-devel+unsubscribe@googlegroups.com.