Refinement of the deterministic part of prime_is_prime

Author: czurnieden

mp_prime_extra_strong_lucas.c

@@ -25,7 +25,6 @@
     Jon Grantham's paper is heartily recommended, it gives some useful insights
     into this theme.
 
-    Comments starting with PARI are pari/gp code, just grep for it.
 
     [1] Baillie, "Extra strong Lucas pseudoprimes", OEIS A217719, https://oeis.org/A217719
     [2] Crandall and Pomerance, Prime Numbers: A Computational Perspective, 2nd ed. Springer, 2005
@@ -134,7 +133,5 @@ mp_err mp_prime_extra_strong_lucas(const mp_int *N, bool *result)
    return err;
 }
 #endif
-
-
-
 #endif
+

mp_prime_is_prime.c

@@ -121,8 +121,16 @@ mp_err mp_prime_is_prime(const mp_int *a, int t, bool *result)
    }
 #endif
 
-   /* run at least one Miller-Rabin test with a random base */
+
    if (t == 0) {
+#ifndef LTM_USE_ONLY_MR
+      /* The BPSW version as used here has no counter-example below 2^64 */
+      if (mp_count_bits(a) < 64) {
+         *result = true;
+         goto LBL_B;
+      }
+#endif
+      /* run at least one Miller-Rabin test with a random base if n > 2^64 */
       t = 1;
    }
 
@@ -131,35 +139,82 @@ mp_err mp_prime_is_prime(const mp_int *a, int t, bool *result)
 
       It uses the bases necessary for a deterministic M-R test if the input is
       smaller than  3317044064679887385961981
-      The caller has to check the size.
-      TODO: can be made a bit finer grained but comparing is not free.
+      The caller has to check the maximum size.
    */
    if (t < 0) {
-      int p_max = 0;
+      int p_max = 0, bits;
+      bits = mp_count_bits(a);
 
-      /*
-          Sorenson, Jonathan; Webster, Jonathan (2015).
-           "Strong Pseudoprimes to Twelve Prime Bases".
-       */
-      /* 0x437ae92817f9fc85b7e5 = 318665857834031151167461 */
-      if ((err =   mp_read_radix(&b, "437ae92817f9fc85b7e5", 16)) != MP_OKAY) {
+#ifndef LTM_USE_ONLY_MR
+      if (bits < 64) {
+         /* Just complete the BPSW test */
+         if ((err = mp_prime_extra_strong_lucas(a, &res)) != MP_OKAY) {
+            goto LBL_B;
+         }
+         *result = res;
          goto LBL_B;
       }
+#else
+      /* Base 2 has been done already at this point */
 
-      if (mp_cmp(a, &b) == MP_LT) {
-         p_max = 12;
-      } else {
+      mp_digit bases32 = {7u, 61u};
+      /* 2, 325, 9375, 28178, 450775, 9780504, 1795265022 found by Jim Sinclair 2011 */
+      uint32_t bases64 = {325ull, 9375ull, 28178ull, 450775ull, 9780504ull, 1795265022ull};
+      if (bits < 32) {
+         for (ix = 0; ix < 2; ix++) {
+            mp_set(&b, bases32[ix]);
+            if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
+               goto LBL_B;
+            }
+            if (!res) {
+               goto LBL_B;
+            }
+         }
+         *result = true;
+         goto LBL_B;
+      } else if ((bits >= 32) && (bits < 64)) {
+         for (ix = 0; ix < 6; ix++) {
+            mp_set_u32(&b, bases64[ix]);
+            if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
+               goto LBL_B;
+            }
+            if (!res) {
+               goto LBL_B;
+            }
+         }
+         *result = true;
+         goto LBL_B;
+      }
+#endif
+      /*
+          Sorenson, Jonathan; Webster, Jonathan,  "Strong Pseudoprimes to Twelve Prime Bases". (2015) https://arxiv.org/abs/1509.00864
+          Z. Zhang, "Finding strong pseudoprimes to several bases," Math. Comp., 70:234 (2001) 863--87
+          Z. Zhang, "Finding C3-strong pseudoprimes," Math. Comp., 74:250 (2005) 1009--1024
+          Zhang, Zhenxiang, "Two kinds of strong pseudoprimes up to 10^36," Math. Comp., 76:260 (2007) 2095--2107
+       */
+      /*
+         Comparing bigints is not free, so give the magnitude of "n" a rough check
+         before spending computing time.
+      */
+      if (bits <= 78) {
+         /* 0x437ae92817f9fc85b7e5 = 318665857834031151167461 */
+         if ((err =   mp_read_radix(&b, "437ae92817f9fc85b7e5", 16)) != MP_OKAY) {
+            goto LBL_B;
+         }
+         if (mp_cmp(a, &b) == MP_LT) {
+            p_max = 12;
+         }
+      } else if ((bits > 78) && (bits < 82)) {
          /* 0x2be6951adc5b22410a5fd = 3317044064679887385961981 */
          if ((err = mp_read_radix(&b, "2be6951adc5b22410a5fd", 16)) != MP_OKAY) {
             goto LBL_B;
          }
-
          if (mp_cmp(a, &b) == MP_LT) {
             p_max = 13;
-         } else {
-            err = MP_VAL;
-            goto LBL_B;
          }
+      } else {
+         err = MP_VAL;
+         goto LBL_B;
       }
 
       /* we did bases 2 and 3  already, skip them */