[flint] 04/24: Rewrite n_is_prime and n_is_probabprime

[felix salfelder]

This is an automated email from the git hooks/post-receive script.

felix-guest pushed a commit to branch master
in repository flint.

commit 12126fb2d24e82930f354ae90c4c9753b8290eeb
Author: Dana Jacobsen <dana at acm.org>
Date:   Tue May 13 01:07:45 2014 -0700

    Rewrite n_is_prime and n_is_probabprime
---
 ulong_extras/is_prime.c       | 18 +++++---
 ulong_extras/is_probabprime.c | 99 +++++++++++++++++++++----------------------
 2 files changed, 59 insertions(+), 58 deletions(-)

diff --git a/ulong_extras/is_prime.c b/ulong_extras/is_prime.c
index 56030b2..126efa4 100644
--- a/ulong_extras/is_prime.c
+++ b/ulong_extras/is_prime.c
@@ -34,12 +34,16 @@ int n_is_prime(mp_limb_t n)
        [1]  http://www.janfeitsma.nl/math/psp2/database
     */
 
-#if !FLINT64
+    if (n < 11) {
+        if (n == 2 || n == 3 || n == 5 || n == 7)   return 1;
+        else                                        return 0;
+    }
+    if (!(n%2) || !(n%3) || !(n%5) || !(n%7))       return 0;
+    if (n <  121) /* 11*11 */                       return 1;
+    if (!(n%11) || !(n%13) || !(n%17) || !(n%19) ||
+        !(n%23) || !(n%29) || !(n%31) || !(n%37) ||
+        !(n%41) || !(n%43) || !(n%47) || !(n%53))   return 0;
+    if (n < 3481) /* 59*59 */                       return 1;
+
     return n_is_probabprime(n);
-#else
-    if (n < UWORD(10000000000000000))
-        return n_is_probabprime(n);
-    else 
-        return n_is_probabprime_BPSW(n);
-#endif
 }
diff --git a/ulong_extras/is_probabprime.c b/ulong_extras/is_probabprime.c
index 4bca967..365abcb 100644
--- a/ulong_extras/is_probabprime.c
+++ b/ulong_extras/is_probabprime.c
@@ -33,79 +33,76 @@
 
     Currently it acts as such all the way up to 2^64.
 */
+
 int n_is_probabprime(mp_limb_t n)
 {
     mp_limb_t d;
     unsigned int norm;
-	mp_limb_t ninv;
+    int isprime;
+    double npre;
 
     if (n <= UWORD(1)) return 0;
     if (n == UWORD(2)) return 1;
     if ((n & UWORD(1)) == 0) return 0;
 
-    if (n_is_perfect_power235(n)) return 0;
+    if (n < FLINT_ODDPRIME_SMALL_CUTOFF)
+        return n_is_oddprime_small(n);
+    if (n < FLINT_PRIMES_TAB_DEFAULT_CUTOFF)
+        return n_is_oddprime_binary(n);
 
 #if FLINT64
-    if (n >= UWORD(10000000000000000)) return n_is_probabprime_BPSW(n);
+    /* Avoid the unnecessary inverse */
+    if (n >= UWORD(1050535501))
+        return n_is_probabprime_BPSW(n);
 #endif
 
+    isprime = 0;
     d = n - 1;
     count_trailing_zeros(norm, d);
     d >>= norm;
+    npre = n_precompute_inverse(n);
 
-#if FLINT64
-    if (n < UWORD(1122004669633))
-#else
-    if (n < UWORD(2147483648))
-#endif  
+#if !FLINT64
+
+    /* For 32-bit, just the 2-base or 3-base Miller-Rabin is enough */
+    if (n < UWORD(9080191))
+    {
+        isprime = n_is_strong_probabprime_precomp(n, npre, UWORD(31), d)
+               && n_is_strong_probabprime_precomp(n, npre, UWORD(73), d);
+    }
+    else
     {
-        double npre;
-        if (n < FLINT_ODDPRIME_SMALL_CUTOFF)
-            return n_is_oddprime_small(n);
-
-        if (n < FLINT_PRIMES_TAB_DEFAULT_CUTOFF)
-            return n_is_oddprime_binary(n);
-      
-        npre = n_precompute_inverse(n);
-
-        if (n < UWORD(9080191)) 
-        {
-            if (n_is_strong_probabprime_precomp(n, npre, UWORD(31), d)
-                && n_is_strong_probabprime_precomp(n, npre, UWORD(73), d)) return 1;
-            else return 0;
-        }
+        isprime = n_is_strong_probabprime_precomp(n, npre, UWORD(2), d)
+               && n_is_strong_probabprime_precomp(n, npre, UWORD(7), d)
+               && n_is_strong_probabprime_precomp(n, npre, UWORD(61), d);
+    }
 
-#if FLINT64
-        if (n < UWORD(4759123141))
-        {
-#endif
-        if (n_is_strong_probabprime_precomp(n, npre, UWORD(2), d) 
-            && n_is_strong_probabprime_precomp(n, npre, UWORD(7), d) 
-            && n_is_strong_probabprime_precomp(n, npre, UWORD(61), d)) return 1;
-        else return 0;
-#if FLINT64
-        }
-
-        if (n_is_strong_probabprime_precomp(n, npre, UWORD(2), d) 
-            && n_is_strong_probabprime_precomp(n, npre, UWORD(13), d) 
-            && n_is_strong_probabprime_precomp(n, npre, UWORD(23), d) 
-            && n_is_strong_probabprime_precomp(n, npre, UWORD(1662803), d))
-            if (n != UWORD(46856248255981)) return 1;
-        return 0;
+#else
+
+    /* For 64-bit, BPSW seems to be a little bit faster than 3 bases. */
+    if (n < UWORD(341531))
+    {
+        isprime = n_is_strong_probabprime_precomp(n, npre, UWORD(9345883071009581737), d);
+    }
+    else if (n < UWORD(1050535501))
+    {
+        isprime = n_is_strong_probabprime_precomp(n, npre, UWORD(336781006125), d)
+               && n_is_strong_probabprime_precomp(n, npre, UWORD(9639812373923155), d);
+    }
+#if 0
+    else if (n < UWORD(350269456337))
+    {
+        isprime = n_is_strong_probabprime_precomp(n, npre, UWORD(4230279247111683200), d)
+               && n_is_strong_probabprime_precomp(n, npre, UWORD(14694767155120705706), d)
+               && n_is_strong_probabprime_precomp(n, npre, UWORD(16641139526367750375), d);
+    }
 #endif
+    else
+    {
+        isprime = n_is_probabprime_BPSW(n);
     }
 
-	ninv = n_preinvert_limb(n);
-
-    if (n_is_strong_probabprime2_preinv(n, ninv, UWORD(2), d) 
-        && n_is_strong_probabprime2_preinv(n, ninv, UWORD(3), d) 
-        && n_is_strong_probabprime2_preinv(n, ninv, UWORD(7), d) 
-        && n_is_strong_probabprime2_preinv(n, ninv, UWORD(61), d) 
-        && n_is_strong_probabprime2_preinv(n, ninv, UWORD(24251), d))
-#if FLINT64
-        if (n != UWORD(46856248255981)) 
 #endif
-        return 1;
 
-    return 0;
+    return isprime;
 }

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/flint.git