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:  sage6.1 
Component:  basic arithmetic  Keywords:  
Cc:  nbruin, saraedum, jpflori  Merged in:  
Authors:  David Roe, Peter Bruin  Reviewers:  Francis Clarke, Peter Bruin, JeanPierre Flori 
Report Upstream:  N/A  Work issues:  
Branch:  u/pbruin/13596IntegerMod_is_square (Commits)  Commit:  6c726cc75a36e8d1b06058efc06942c910bb2515 
Dependencies:  #15193, #11868  Stopgaps: 
Description (last modified by )
Currently we don't check first to see if the Jacobi symbol is 1 (in which case we can declare nonsquareness 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)
Change History (26)
Changed 8 years ago by
comment:1 Changed 8 years ago by
 Status changed from new to needs_review
comment:2 Changed 8 years ago by
 Cc nbruin saraedum added
comment:3 Changed 8 years ago by
 Reviewers set to Francis Clarke
 Status changed from needs_review to positive_review
comment:4 Changed 8 years ago by
 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
 Dependencies set to #12415
 Status changed from needs_work to needs_review
comment:6 Changed 8 years ago by
 Milestone changed from sage5.5 to sagepending
comment:7 Changed 8 years ago by
 Cc jpflori added
comment:8 followups: ↓ 9 ↓ 10 Changed 7 years ago by
 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
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 ; followup: ↓ 12 Changed 7 years ago by
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: 20130617 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
 Milestone changed from sagepending to sage5.12
comment:12 in reply to: ↑ 10 Changed 7 years ago by
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
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 factoringn
(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 preinvented 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.
comment:14 Changed 7 years ago by
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
After experimenting with various mixtures of Cython and PARI, I believe that the most efficient solution is to use Cython for those checks for nonsquareness 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
 Dependencies changed from #12415 to #15193
 Description modified (diff)
 Status changed from needs_work to needs_review
comment:17 Changed 7 years ago by
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.)
comment:18 Changed 7 years ago by
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
Patchbot:
apply 13596_new.patch
comment:20 Changed 7 years ago by
 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
 Branch set to u/pbruin/13596IntegerMod_is_square
comment:22 Changed 7 years ago by
 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 nonsquareness that do not require factoring. * Use PARI otherwise.

comment:23 Changed 7 years ago by
 Reviewers changed from Francis Clarke, Peter Bruin to Francis Clarke, Peter Bruin, JeanPierre 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
 Resolution set to fixed
 Status changed from positive_review to closed
Looks good. Passes doctests. Much faster when Jacobi symbol is 1 and no slower otherwise.
Positive review.