Opened 8 years ago

Closed 7 years ago

#13596 closed enhancement (fixed)

Improvements to IntegerMod is_square

Reported by: roed Owned by: AlexGhitza
Priority: major Milestone: sage-6.1
Component: basic arithmetic Keywords:
Cc: nbruin, saraedum, jpflori Merged in:
Authors: David Roe, Peter Bruin Reviewers: Francis Clarke, Peter Bruin, Jean-Pierre Flori
Report Upstream: N/A Work issues:
Branch: u/pbruin/13596-IntegerMod_is_square (Commits) Commit: 6c726cc75a36e8d1b06058efc06942c910bb2515
Dependencies: #15193, #11868 Stopgaps:

Description (last modified by pbruin)

Currently we don't check first to see if the Jacobi symbol is -1 (in which case we can declare non-squareness without factoring the modulus). This patch fixes this. If this test does not prove that the element is not a square, we now call PARI's Zn_issquare function with the factored modulus.

Attachments (2)

13596.patch (3.9 KB) - added by roed 8 years ago.
13596_new.patch (7.4 KB) - added by pbruin 7 years ago.
based on 5.12.beta4 + #15124

Download all attachments as: .zip

Change History (26)

Changed 8 years ago by roed

comment:1 Changed 8 years ago by roed

  • Status changed from new to needs_review

comment:2 Changed 8 years ago by roed

  • Cc nbruin saraedum added

comment:3 Changed 8 years ago by fwclarke

  • Reviewers set to Francis Clarke
  • Status changed from needs_review to positive_review

Looks good. Passes doctests. Much faster when Jacobi symbol is -1 and no slower otherwise.

Positive review.

comment:4 Changed 8 years ago by jdemeyer

  • Status changed from positive_review to needs_work
sage -t --long "devel/sage/sage/rings/finite_rings/integer_mod.pyx"
Exception raised by doctesting framework. Use -verbose for details.
         [34.4 s]

----------------------------------------------------------------------
The following tests failed:


        sage -t --long "devel/sage/sage/rings/finite_rings/integer_mod.pyx" # Exception from doctest framework
Total time for all tests: 34.6 seconds

comment:5 Changed 8 years ago by roed

  • Dependencies set to #12415
  • Status changed from needs_work to needs_review

Sorry about that. I'm using syntax that depends on #12415. I'm going to mark this as depending on #12415, though if that ticket continues to stall I can change the test here instead.

comment:6 Changed 8 years ago by jdemeyer

  • Milestone changed from sage-5.5 to sage-pending

comment:7 Changed 8 years ago by jpflori

  • Cc jpflori added

comment:8 follow-ups: Changed 7 years ago by pbruin

  • Reviewers changed from Francis Clarke to Francis Clarke, Peter Bruin
  • Status changed from needs_review to needs_work

The patch did not apply anymore due to the fact that the pari() method was renamed to _pari_(). I fixed this, and then did the following (inspired by the new doctest):

def test_native():
    for p,q,r in cartesian_product_iterator([[3,5],[11,13],[17,19]]):
        for ep,eq,er in cartesian_product_iterator([[0,1,2,3],[0,1,2,3],[0,1,2,3]]):
            for e2 in [0,1,2,3,4]:
                n = p^ep*q^eq*r^er*2^e2
                for _ in range(2):
                    a = Zmod(n).random_element()
                    b = a.is_square()

def test_pari():
    for p,q,r in cartesian_product_iterator([[3,5],[11,13],[17,19]]):
        for ep,eq,er in cartesian_product_iterator([[0,1,2,3],[0,1,2,3],[0,1,2,3]]):
            for e2 in [0,1,2,3,4]:
                n = p^ep*q^eq*r^er*2^e2
                for _ in range(2):
                    a = Zmod(n).random_element()
                    b = bool(a._pari_().issquare())

I got the following timings:

sage: %timeit test_native()
1 loops, best of 3: 3.11 s per loop
sage: %timeit test_pari()
1 loops, best of 3: 2.98 s per loop

From the perspective of someone who does not want to reinvent the wheel, this suggests that we should just call PARI instead of using a native Sage implementation.

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

Replying to pbruin:

From the perspective of someone who does not want to reinvent the wheel, this suggests that we should just call PARI instead of using a native Sage implementation.

+1

comment:10 in reply to: ↑ 8 ; follow-up: Changed 7 years ago by fwclarke

Replying to pbruin:

I got the following timings:

The following case makes me think that the examples in the doctest are rather too small to rely on for timing information.

sage: print version()
Sage Version 5.10, Release Date: 2013-06-17
sage: a, b, n = 176338465705, 161089773287, 217267613611
sage: n.factor()
341983 * 635317
sage: jacobi_symbol(a, n), jacobi_symbol(b, n)
(1, -1)
sage: %timeit Zmod(n)(a).is_square()
100000 loops, best of 3: 10.2 us per loop
sage: %timeit Zmod(n)(a)._pari_().issquare()
1000 loops, best of 3: 215 us per loop
sage: %timeit Zmod(n)(b).is_square()
100000 loops, best of 3: 8.83 us per loop
sage: %timeit Zmod(n)(b)._pari_().issquare()
100000 loops, best of 3: 14.8 us per loop

It certainly looks as though PARI checks the Jacobi symbol first, but even without the patch the current code wins. Perhaps it is worth reinventing this wheel.

comment:11 Changed 7 years ago by jdemeyer

  • Milestone changed from sage-pending to sage-5.12

comment:12 in reply to: ↑ 10 Changed 7 years ago by pbruin

Replying to fwclarke:

Replying to pbruin:

The following case makes me think that the examples in the doctest are rather too small to rely on for timing information.

It certainly looks as though PARI checks the Jacobi symbol first, but even without the patch the current code wins. Perhaps it is worth reinventing this wheel.

What is behind this is that Sage caches the factorisation of n in Zmod(n), while PARI has to factor n every time. With your values of a and n, the PARI version used by Sage on my system needs 367 us for issquare(Mod(a, n)), while factor(n) takes 360 us. This means that 98% of the cost of issquare() in this case is spent on factoring n.

So for a single call, PARI is probably slightly faster, while the Sage implementation is faster for multiple calls within one Zmod(n) thanks to the cached factorisation.

comment:13 Changed 7 years ago by pbruin

Here are some timings (on Linux x86_64):

a = 176338465705
b = 161089773287
n = 217267613611

issquare(Mod(a, n)) resp. Zmod(n)(a).is_square()
pari 2.5.4: 367 us
pari 2.6.1: 396 us
sage 5.11.rc0: 8.74 us
sage + #13596: 25 us

issquare(Mod(b, n)) resp. Zmod(n)(b).is_square()
pari 2.5.4: 1.18 us
pari 2.6.1: 155 us
sage 5.11.rc0: 15.2 us
sage + #13596: 10.6 us

kronecker(a, n) resp. a.jacobi(n)
pari 2.5.4: 0.960 us
pari 2.6.1: 0.961 us
sage 5.11.rc0: 0.693 us

kronecker(b, n) resp. b.jacobi(n)
pari 2.5.4: 0.970 us
pari 2.6.1: 0.981 us
sage 5.11.rc0: 0.709 us

factorint(n) resp. n.factor()
pari 2.5.4: 360 us
pari 2.6.1: 330 us
sage 5.11.rc0: 157 us

Note:

  • almost all the time for issquare in PARI is spent on factoring n (see above).
  • there is a dramatic slowdown between PARI 2.5.4 and 2.6.1 for issquare(Mod(b, n)). (Edit: this has been fixed in the latest Git version of PARI.)
  • Kronecker symbols and factoring are faster in Sage than in PARI for this value of n; this is because Sage uses the right pre-invented wheel (GMP resp. Flint).
  • this patch slows down the Sage function by a factor of 3 for a mod n, but speeds it up by 30% for b mod n.
Last edited 7 years ago by pbruin (previous) (diff)

comment:14 Changed 7 years ago by pbruin

It turns out that PARI has a function Zn_issquare(d, n) to check whether d is a square modulo n; for n it accepts either an integer or its factorisation matrix. This looks like a nice preinvented wheel to me.

comment:15 Changed 7 years ago by pbruin

After experimenting with various mixtures of Cython and PARI, I believe that the most efficient solution is to use Cython for those checks for non-squareness that don't require factoring, and if that fails, call PARI's Zn_issquare function with the modulus in factored form. I made a patch and am now running doctests.

There seems to be a very slight slowdown with respect to the current implementation if the element is a square (perhaps because the above quick tests are done once in Cython and again in PARI), but a noticeable speedup if it is not.

comment:16 Changed 7 years ago by pbruin

  • Authors changed from David Roe to David Roe, Peter Bruin
  • Dependencies changed from #12415 to #15193
  • Description modified (diff)
  • Status changed from needs_work to needs_review

Changed 7 years ago by pbruin

based on 5.12.beta4 + #15124

comment:17 Changed 7 years ago by pbruin

Again some timings:

a = 176338465705
b = 161089773287
n = 217267613611

def test1():
    for p,q,r in cartesian_product_iterator([[3,5],[11,13],[17,19]]):
        for ep,eq,er in cartesian_product_iterator([[0,1,2,3],[0,1,2,3],[0,1,2,3]]):
            for e2 in [0,1,2,3,4]:
                n = p^ep*q^eq*r^er*2^e2
                for _ in range(2):
                    a = Zmod(n).random_element()
                    b = a.is_square()

def test2(bound):
    for n in xrange(bound):
        for a in xrange(n):
            b = Mod(a, n).is_square()

# without patch
sage: %timeit -c -n 100000 -r 1 Zmod(n)(a).is_square()
100000 loops, best of 1: 11.5 us per loop
sage: %timeit -c -n 100000 -r 1 Zmod(n)(b).is_square()
100000 loops, best of 1: 9.64 us per loop
sage -t --long sage/rings/finite_rings/integer_mod.pyx
cpu time: 5.5 seconds
sage: %timeit -c -r 1 test1()
1 loops, best of 1: 1.27 s per loop
sage: %timeit -c -r 1 test2(1000)
1 loops, best of 1: 7.12 s per loop

# with patch
sage: %timeit -c -n 100000 -r 1 Zmod(n)(a).is_square()
100000 loops, best of 1: 12.1 us per loop
sage: %timeit -c -n 100000 -r 1 Zmod(n)(b).is_square()
100000 loops, best of 1: 6.52 us per loop
sage -t --long sage/rings/finite_rings/integer_mod.pyx
cpu time: 7.8 seconds
sage: %timeit -c -r 1 test1()
1 loops, best of 1: 1.28 s per loop
sage: %timeit -c -r 1 test2(1000)
1 loops, best of 1: 6.84 s per loop

(On x86_64 GNU/Linux, but a different machine than in comment:13. The increase in doctest time is because of the new test in this patch.)

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

comment:18 Changed 7 years ago by jpflori

Hi, I didn't have time to go through this ticket, but FYI I've reimplemented the modular squareroot in cython using mpz's (translating the original python/cython implem) at least when the mod is prime or 4 times a prime as part of an implementation of ECPP. It's much faster than what we have and faster than PARI IIRC. This should get published together with the rest of the ECPP implem when the code is a little more polished, let's say before the end of the year hopefully. (We'll also present it in Paris at some Sage afternoon next week.)

comment:19 Changed 7 years ago by pbruin

Patchbot:

apply 13596_new.patch

comment:20 Changed 7 years ago by pbruin

  • Created changed from 10/13/12 02:05:55 to 10/13/12 02:05:55
  • Dependencies changed from #15193 to #15193, #11868
  • Description modified (diff)
  • Modified changed from 09/23/13 10:30:22 to 09/23/13 10:30:22

comment:21 Changed 7 years ago by pbruin

  • Branch set to u/pbruin/13596-IntegerMod_is_square

comment:22 Changed 7 years ago by git

  • Commit set to 6c726cc75a36e8d1b06058efc06942c910bb2515

Branch pushed to git repo; I updated commit sha1. This was a forced push. Recent commits:

6c726cc* First do checks for non-squareness that do not require factoring. * Use PARI otherwise.

comment:23 Changed 7 years ago by jpflori

  • Reviewers changed from Francis Clarke, Peter Bruin to Francis Clarke, Peter Bruin, Jean-Pierre Flori
  • Status changed from needs_review to positive_review

Looks good to me, nice work. Successfully built and tested on one system as well.

comment:24 Changed 7 years ago by vbraun

  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.