Improvements to IntegerMod is_square

[roed]

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.

[fwclarke]

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

Positive review.

[jdemeyer]

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

[roed]

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.

[pbruin]

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.

[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

[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.

[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.

[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.

[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.

[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.

[pbruin]

based on 5.12.beta4 + #15124

[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.)

[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.)

[pbruin]

Patchbot:

apply 13596_new.patch

[git]

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.

[jpflori]

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