From 6fbbc0e043cb8c7b6cbd96bbe6d0118feb32ec45 Mon Sep 17 00:00:00 2001
From: czurnieden <czurnieden@gmx.de>
Date: Tue, 4 Apr 2023 06:39:39 +0200
Subject: [PATCH] Further refinement of the deterministic part of
 prime_is_prime

---
 mp_prime_is_prime.c | 56 +++++++++++++++++++++++++++++++++------------
 1 file changed, 42 insertions(+), 14 deletions(-)

diff --git a/mp_prime_is_prime.c b/mp_prime_is_prime.c
index 743c379c..cbe30330 100644
--- a/mp_prime_is_prime.c
+++ b/mp_prime_is_prime.c
@@ -16,9 +16,9 @@ static unsigned int s_floor_ilog2(int value)
 mp_err mp_prime_is_prime(const mp_int *a, int t, bool *result)
 {
    mp_int  b;
-   int     ix;
+   int     ix, bits;
    bool    res;
-   mp_err  err;
+   mp_err  err = MP_OKAY;
 
    /* default to no */
    *result = false;
@@ -59,11 +59,19 @@ mp_err mp_prime_is_prime(const mp_int *a, int t, bool *result)
    if ((err = s_mp_prime_is_divisible(a, &res)) != MP_OKAY) {
       return err;
    }
-
    /* return if it was trivially divisible */
    if (res) {
       return MP_OKAY;
    }
+   /* floor(log_2(a)) */
+   bits = mp_count_bits(a) - 1;
+
+   /* If the whole prime table up to p = 1619 has been tested, than  all
+      numbers below 1621^2 = 2,627,641 are prime now. log_2(1621^2) ~ 21.33 */
+   if (bits < 21) {
+      *result = true;
+      return MP_OKAY;
+   }
 
    /*
        Run the Miller-Rabin test with base 2 for the BPSW test.
@@ -78,6 +86,15 @@ mp_err mp_prime_is_prime(const mp_int *a, int t, bool *result)
    if (!res) {
       goto LBL_B;
    }
+   /* If the whole prime table up to p = 1619 and the Miller-Rabin tests to base two
+      has been applied, than all numbers below 4,469,471 (~2^{22.1}) are prime now.
+      With 1659 SPSPs < 2^32 left */
+#if ((defined S_MP_PRIME_IS_DIVISIBLE_C) && (MP_PRIME_TAB_SIZE >= 256))
+   if (bits < 22) {
+      *result = true;
+      goto LBL_B;
+   }
+#endif
    /*
       Rumours have it that Mathematica does a second M-R test with base 3.
       Other rumours have it that their strong L-S test is slightly different.
@@ -91,6 +108,15 @@ mp_err mp_prime_is_prime(const mp_int *a, int t, bool *result)
       goto LBL_B;
    }
 
+   /* If the whole prime table up to p = 1619 and the Miller-Rabin tests to bases
+      two and three have been applied, than all numbers below 11,541,307 (~2^{23.5}) are prime now.
+      With 89 SPSPs < 2^32 left */
+#if ((defined S_MP_PRIME_IS_DIVISIBLE_C) && (MP_PRIME_TAB_SIZE >= 256))
+   if (bits < 23) {
+      *result = true;
+      goto LBL_B;
+   }
+#endif
    /*
     * Both, the Frobenius-Underwood test and the the extra strong Lucas test are quite
     * slow so if speed is an issue, define LTM_USE_ONLY_MR to use M-R tests with
@@ -125,7 +151,7 @@ mp_err mp_prime_is_prime(const mp_int *a, int t, bool *result)
    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) {
+      if (bits < 64) {
          *result = true;
          goto LBL_B;
       }
@@ -142,9 +168,7 @@ mp_err mp_prime_is_prime(const mp_int *a, int t, bool *result)
       The caller has to check the maximum size.
    */
    if (t < 0) {
-      int p_max = 0, bits;
-      bits = mp_count_bits(a);
-
+      int p_max = 0;
 #ifndef LTM_USE_ONLY_MR
       if (bits < 64) {
          /* Just complete the BPSW test */
@@ -155,8 +179,8 @@ mp_err mp_prime_is_prime(const mp_int *a, int t, bool *result)
          goto LBL_B;
       }
 #else
-      /* Base 2 has been done already at this point */
-
+      /* Base 2 has been done already at this point. Also possible: { (2, 3,) 5, 13} to avoid waste */
+      /* 2, 7, 61 found by Gerhard Jaeschke 1993 */
       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};
@@ -196,21 +220,27 @@ mp_err mp_prime_is_prime(const mp_int *a, int t, bool *result)
          Comparing bigints is not free, so give the magnitude of "n" a rough check
          before spending computing time.
       */
-      if (bits <= 78) {
+
+      else if ((bits >= 64) && (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 {
+            p_max = 13;
          }
-      } else if ((bits > 78) && (bits < 82)) {
+      } else if ((bits >= 78) && (bits <= 81)) {
          /* 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;
@@ -232,10 +262,9 @@ mp_err mp_prime_is_prime(const mp_int *a, int t, bool *result)
        Do "t" M-R tests with random bases between 3 and "a".
        See Fips 186.4 p. 126ff
    */
-   else if (t > 0) {
+   if (t > 0) {
       unsigned int mask;
       int size_a;
-
       /*
        * The mp_digit's have a defined bit-size but the size of the
        * array a.dp is a simple 'int' and this library can not assume full
@@ -283,7 +312,6 @@ mp_err mp_prime_is_prime(const mp_int *a, int t, bool *result)
       for (ix = 0; ix < t; ix++) {
          unsigned int fips_rand;
          int len;
-
          /* mp_rand() guarantees the first digit to be non-zero */
          if ((err = mp_rand(&b, 1)) != MP_OKAY) {
             goto LBL_B;