Faster primality test between 65 and 128 bits

doc/source/ulong_extras.rst

@@ -776,6 +776,33 @@ precomputed quotient for `b`:
 5. compute `h = b \check{a} - q n` (two single-word multiplications)
 6. if `h \ge n` then `q \gets q+1`
 
+Double-limb modular arithmetic
+--------------------------------------------------------------------------------
+
+..  function:: void n_ll_small_2_powmod(nn_ptr res, nn_srcptr exp, nn_srcptr m, nn_srcptr minv)
+               void n_ll_small_powmod_triple(nn_ptr res1, nn_ptr res2, nn_ptr res3, ulong b1, ulong b2, ulong b3, nn_srcptr exp, nn_srcptr m, nn_srcptr minv)
+
+    Computes `b^e \bmod m` where `b` is a small single-limb integer,
+    `e` is a double-limb integer stored in the array ``exp`` and
+    `m` is a double-limb integer accompanied by the two-limb
+    precomputed inverse `\lfloor 2^{3 \cdot \mathtt{FLINT\_BITS}} / \rfloor`
+    stored in ``minv``. The function :func:`n_ll_small_2_powmod` is
+    specialized for base `b = 2`, writing the double-limb result to
+    ``res``, while :func:`n_ll_small_powmod_triple`
+    performs three simultaneous exponentations, writing the results
+    for bases `b_1 \le b_2 \le b_3` to ``res1``, ``res2`` and ``res3``
+    respectively.
+
+    We require `m > 2^{\mathtt{FLINT\_BITS}}`. In addition, both `m` and
+    `b` must be "small" for correctness, i.e. the whole double-limb
+    operand range is not supported.
+    A sufficient criterion for correctness is that
+    `(12 m)^2 < 2^{3 \cdot \mathtt{FLINT\_BITS}}` for the base-2 test
+    and `(6 b_3 m)^2 < 2^{3 \cdot \mathtt{FLINT\_BITS}}` for the general
+    test.
+
+    These are helper functions used by :func:`n_ll_is_prime`.
+
 Divisibility testing
 --------------------------------------------------------------------------------
 
@@ -1090,6 +1117,26 @@ Primality testing
 
     This function is obsolete and currently just wraps :func:`n_is_prime`.
 
+.. function:: int n_ll_is_prime(ulong nhi, ulong nlo)
+
+    Primality test for a double-limb integer `n` represented by
+    low part ``nlo`` and high part ``nhi``. The high part must be nonzero.
+    Returns 1 if `n` is certainly prime, 0 if `n` is certainly composite,
+    and -1 if unknown.
+
+    For `n` up to about 81 bits on a 64-bit machine, this function first does
+    trial division and then performs a strong probable prime test (Miller-Rabin
+    test) with the first 13 primes as witnesses. This has been shown to prove
+    primality for integers in this range [SorWeb2016]_, so the return
+    value in this range is always 0 or 1.
+
+    For larger `n` on a 64-bit machine, this function does trial division
+    and a base-2 test, returning either 0 or -1.
+    A return value of -1 thus indicates that `n` is at least a base
+    strong probable prime. Users may fall back on
+    :func:`fmpz_is_prime` for a proved result in this case.
+
+    On 32-bit machines, this function currently always returns -1.
 
 Chinese remaindering
 --------------------------------------------------------------------------------

src/fmpz/is_prime.c

@@ -25,7 +25,7 @@ static int _fmpz_is_prime(const fmpz_t n, int proved)
    const ulong * primes;
    const double * pinv;
 
-   fmpz_t F1, Fsqr, Fcub, R, t;
+   fmpz_t F1, Fsqr, Fcub, R;
    int res = -1;
 
     if (!COEFF_IS_MPZ(*n))
@@ -34,9 +34,8 @@ static int _fmpz_is_prime(const fmpz_t n, int proved)
        if (v <= 1)
            return 0;
 
-       /* Note: we want n_is_prime rather than n_is_probabprime
-          even when proved == 0, because ns_is_prime handles general
-          input faster. */
+       /* Note: n_is_prime and n_is_probabprime are currently identical,
+          so we ignore the proved flag. */
        return n_is_prime(v);
    }
    else
@@ -57,45 +56,39 @@ static int _fmpz_is_prime(const fmpz_t n, int proved)
       if (d[0] % 2 == 0)
           return 0;
 
-      bits = size * FLINT_BITS + FLINT_BIT_COUNT(d[size-1]);
-      trial_primes = bits;
+      if (size == 2 && FLINT_BITS == 64)
+      {
+         /* n_ll_is_prime does trial division, a base-2 sprp test, and may
+            additionally be able to certify some inputs. Currently the
+            certifications in n_ll_is_prime are faster than a Lucas test,
+            so use it even when !proved. */
+         res = n_ll_is_prime(d[1], d[0]);
+         if (res != -1)
+            return res;
+      }
+      else
+      {
+          bits = size * FLINT_BITS + FLINT_BIT_COUNT(d[size-1]);
+          trial_primes = bits;
 
-      if (flint_mpn_factor_trial(d, size, 1, trial_primes))
-           return 0;
-    }
+          if (flint_mpn_factor_trial(d, size, 1, trial_primes))
+               return 0;
 
-   /* todo: use fmpz_is_perfect_power? */
-   if (fmpz_is_square(n))
-      return 0;
+           /* todo: use fmpz_is_perfect_power? */
+           if (fmpz_is_square(n))
+              return 0;
 
-   if (!proved)
-      return fmpz_is_probabprime_BPSW(n);
-
-   /* Fast deterministic Miller-Rabin test up to about 81 bits. This choice of
-      bases certifies primality for n < 3317044064679887385961981;
-      see https://doi.org/10.1090/mcom/3134 */
-   fmpz_init(t);
-   fmpz_tdiv_q_2exp(t, n, 64);
-   if (fmpz_cmp_ui(t, 179817) < 0)
-   {
-      static const char bases[] = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 0 };
-
-      for (i = 0; bases[i] != 0; i++)
-      {
-         fmpz_set_ui(t, bases[i]);
-         if (!fmpz_is_strong_probabprime(n, t))
-             return 0;  /* no need to clear t since it is small */
+           /* Do a single base-2 test to rule out most composites */
+           fmpz base2 = 2;
+           if (!fmpz_is_strong_probabprime(n, &base2))
+             return 0;
       }
-
-      return 1;
    }
 
-   /* Do a single base-2 test to rule out most composites */
-   fmpz_set_ui(t, 2);
-   if (!fmpz_is_strong_probabprime(n, t))
-     return 0;
-
-   fmpz_clear(t);
+   /* At this point n has no small factor and is at least a base-2 sprp.
+      Adding a Lucas test makes this a BPSW test. */
+   if (!proved)
+     return fmpz_is_probabprime_lucas(n);
 
    logd = fmpz_dlog(n);
    limit = (ulong) (logd*logd*logd/100.0) + 20;

src/ulong_extras.h

@@ -371,6 +371,11 @@ ulong n_invmod(ulong x, ulong y)
 
 ulong n_binvert(ulong a);
 
+void n_ll_small_preinv(nn_ptr minv, nn_srcptr m);
+void n_ll_small_powmod_triple(nn_ptr res1, nn_ptr res2, nn_ptr res3, ulong b1,
+    ulong b2, ulong b3, nn_srcptr exp, nn_srcptr m, nn_srcptr minv);
+void n_ll_small_2_powmod(nn_ptr res, nn_srcptr exp, nn_srcptr m, nn_srcptr minv);
+
 /* Modular multiplication with fixed operand **********************************/
 
 ULONG_EXTRAS_INLINE
@@ -504,6 +509,8 @@ void n_prime_pi_bounds(ulong *lo, ulong *hi, ulong n);
 
 ulong n_nextprime(ulong n, int FLINT_UNUSED(proved));
 
+int n_ll_is_prime(ulong nhi, ulong nlo);
+
 /* Factorisation *************************************************************/
 
 #define FLINT_FACTOR_TRIAL_PRIMES_BEFORE_PRIMALITY_TEST 64