mp_prime_is_prime.c

   1 #include "tommath_private.h"
   2 #ifdef MP_PRIME_IS_PRIME_C
   3 /* LibTomMath, multiple-precision integer library -- Tom St Denis */
   4 /* SPDX-License-Identifier: Unlicense */
   5
   6 /* portable integer log of two with small footprint */
   7 static unsigned int s_floor_ilog2(int value)
   8 {
   9    unsigned int r = 0;
  10    while ((value >>= 1) != 0) {
  11       r++;
  12    }
  13    return r;
  14 }
  15
  16 mp_err mp_prime_is_prime(const mp_int *a, int t, bool *result)
  17 {
  18    mp_int  b;
  19    int     ix;
  20    bool    res;
  21    mp_err  err;
  22
  23    /* default to no */
  24    *result = false;
  25
  26    /* Some shortcuts */
  27    /* N > 3 */
  28    if (a->used == 1) {
  29       if ((a->dp[0] == 0u) || (a->dp[0] == 1u)) {
  30          *result = false;
  31          return MP_OKAY;
  32       }
  33       if (a->dp[0] == 2u) {
  34          *result = true;
  35          return MP_OKAY;
  36       }
  37    }
  38
  39    /* N must be odd */
  40    if (mp_iseven(a)) {
  41       return MP_OKAY;
  42    }
  43    /* N is not a perfect square: floor(sqrt(N))^2 != N */
  44    if ((err = mp_is_square(a, &res)) != MP_OKAY) {
  45       return err;
  46    }
  47    if (res) {
  48       return MP_OKAY;
  49    }
  50
  51    /* is the input equal to one of the primes in the table? */
  52    for (ix = 0; ix < MP_PRIME_TAB_SIZE; ix++) {
  53       if (mp_cmp_d(a, s_mp_prime_tab[ix]) == MP_EQ) {
  54          *result = true;
  55          return MP_OKAY;
  56       }
  57    }
  58    /* first perform trial division */
  59    if ((err = s_mp_prime_is_divisible(a, &res)) != MP_OKAY) {
  60       return err;
  61    }
  62
  63    /* return if it was trivially divisible */
  64    if (res) {
  65       return MP_OKAY;
  66    }
  67
  68    /*
  69        Run the Miller-Rabin test with base 2 for the BPSW test.
  70     */
  71    if ((err = mp_init_set(&b, 2uL)) != MP_OKAY) {
  72       return err;
  73    }
  74
  75    if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
  76       goto LBL_B;
  77    }
  78    if (!res) {
  79       goto LBL_B;
  80    }
  81    /*
  82       Rumours have it that Mathematica does a second M-R test with base 3.
  83       Other rumours have it that their strong L-S test is slightly different.
  84       It does not hurt, though, beside a bit of extra runtime.
  85    */
  86    b.dp[0]++;
  87    if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
  88       goto LBL_B;
  89    }
  90    if (!res) {
  91       goto LBL_B;
  92    }
  93
  94    /*
  95     * Both, the Frobenius-Underwood test and the the Lucas-Selfridge test are quite
  96     * slow so if speed is an issue, define LTM_USE_ONLY_MR to use M-R tests with
  97     * bases 2, 3 and t random bases.
  98     */
  99 #ifndef LTM_USE_ONLY_MR
 100    if (t >= 0) {
 101 #ifdef LTM_USE_FROBENIUS_TEST
 102       err = mp_prime_frobenius_underwood(a, &res);
 103       if ((err != MP_OKAY) && (err != MP_ITER)) {
 104          goto LBL_B;
 105       }
 106       if (!res) {
 107          goto LBL_B;
 108       }
 109 #else
 110       if ((err = mp_prime_strong_lucas_selfridge(a, &res)) != MP_OKAY) {
 111          goto LBL_B;
 112       }
 113       if (!res) {
 114          goto LBL_B;
 115       }
 116 #endif
 117    }
 118 #endif
 119
 120    /* run at least one Miller-Rabin test with a random base */
 121    if (t == 0) {
 122       t = 1;
 123    }
 124
 125    /*
 126       Only recommended if the input range is known to be < 3317044064679887385961981
 127
 128       It uses the bases necessary for a deterministic M-R test if the input is
 129       smaller than  3317044064679887385961981
 130       The caller has to check the size.
 131       TODO: can be made a bit finer grained but comparing is not free.
 132    */
 133    if (t < 0) {
 134       int p_max = 0;
 135
 136       /*
 137           Sorenson, Jonathan; Webster, Jonathan (2015).
 138            "Strong Pseudoprimes to Twelve Prime Bases".
 139        */
 140       /* 0x437ae92817f9fc85b7e5 = 318665857834031151167461 */
 141       if ((err =   mp_read_radix(&b, "437ae92817f9fc85b7e5", 16)) != MP_OKAY) {
 142          goto LBL_B;
 143       }
 144
 145       if (mp_cmp(a, &b) == MP_LT) {
 146          p_max = 12;
 147       } else {
 148          /* 0x2be6951adc5b22410a5fd = 3317044064679887385961981 */
 149          if ((err = mp_read_radix(&b, "2be6951adc5b22410a5fd", 16)) != MP_OKAY) {
 150             goto LBL_B;
 151          }
 152
 153          if (mp_cmp(a, &b) == MP_LT) {
 154             p_max = 13;
 155          } else {
 156             err = MP_VAL;
 157             goto LBL_B;
 158          }
 159       }
 160
 161       /* we did bases 2 and 3  already, skip them */
 162       for (ix = 2; ix < p_max; ix++) {
 163          mp_set(&b, s_mp_prime_tab[ix]);
 164          if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
 165             goto LBL_B;
 166          }
 167          if (!res) {
 168             goto LBL_B;
 169          }
 170       }
 171    }
 172    /*
 173        Do "t" M-R tests with random bases between 3 and "a".
 174        See Fips 186.4 p. 126ff
 175    */
 176    else if (t > 0) {
 177       unsigned int mask;
 178       int size_a;
 179
 180       /*
 181        * The mp_digit's have a defined bit-size but the size of the
 182        * array a.dp is a simple 'int' and this library can not assume full
 183        * compliance to the current C-standard (ISO/IEC 9899:2011) because
 184        * it gets used for small embedded processors, too. Some of those MCUs
 185        * have compilers that one cannot call standard compliant by any means.
 186        * Hence the ugly type-fiddling in the following code.
 187        */
 188       size_a = mp_count_bits(a);
 189       mask = (1u << s_floor_ilog2(size_a)) - 1u;
 190       /*
 191          Assuming the General Rieman hypothesis (never thought to write that in a
 192          comment) the upper bound can be lowered to  2*(log a)^2.
 193          E. Bach, "Explicit bounds for primality testing and related problems,"
 194          Math. Comp. 55 (1990), 355-380.
 195
 196             size_a = (size_a/10) * 7;
 197             len = 2 * (size_a * size_a);
 198
 199          E.g.: a number of size 2^2048 would be reduced to the upper limit
 200
 201             floor(2048/10)*7 = 1428
 202             2 * 1428^2       = 4078368
 203
 204          (would have been ~4030331.9962 with floats and natural log instead)
 205          That number is smaller than 2^28, the default bit-size of mp_digit.
 206       */
 207
 208       /*
 209         How many tests, you might ask? Dana Jacobsen of Math::Prime::Util fame
 210         does exactly 1. In words: one. Look at the end of _GMP_is_prime() in
 211         Math-Prime-Util-GMP-0.50/primality.c if you do not believe it.
 212
 213         The function mp_rand() goes to some length to use a cryptographically
 214         good PRNG. That also means that the chance to always get the same base
 215         in the loop is non-zero, although very low.
 216         If the BPSW test and/or the additional Frobenious test have been
 217         performed instead of just the Miller-Rabin test with the bases 2 and 3,
 218         a single extra test should suffice, so such a very unlikely event
 219         will not do much harm.
 220
 221         To preemptively answer the dangling question: no, a witness does not
 222         need to be prime.
 223       */
 224       for (ix = 0; ix < t; ix++) {
 225          unsigned int fips_rand;
 226          int len;
 227
 228          /* mp_rand() guarantees the first digit to be non-zero */
 229          if ((err = mp_rand(&b, 1)) != MP_OKAY) {
 230             goto LBL_B;
 231          }
 232          /*
 233           * Reduce digit before casting because mp_digit might be bigger than
 234           * an unsigned int and "mask" on the other side is most probably not.
 235           */
 236          fips_rand = (unsigned int)(b.dp[0] & (mp_digit) mask);
 237          if (fips_rand > (unsigned int)(INT_MAX - MP_DIGIT_BIT)) {
 238             len = INT_MAX / MP_DIGIT_BIT;
 239          } else {
 240             len = (((int)fips_rand + MP_DIGIT_BIT) / MP_DIGIT_BIT);
 241          }
 242          /*  Unlikely. */
 243          if (len < 0) {
 244             ix--;
 245             continue;
 246          }
 247          if ((err = mp_rand(&b, len)) != MP_OKAY) {
 248             goto LBL_B;
 249          }
 250          /*
 251           * That number might got too big and the witness has to be
 252           * smaller than "a"
 253           */
 254          len = mp_count_bits(&b);
 255          if (len >= size_a) {
 256             len = (len - size_a) + 1;
 257             if ((err = mp_div_2d(&b, len, &b, NULL)) != MP_OKAY) {
 258                goto LBL_B;
 259             }
 260          }
 261          /* Although the chance for b <= 3 is miniscule, try again. */
 262          if (mp_cmp_d(&b, 3uL) != MP_GT) {
 263             ix--;
 264             continue;
 265          }
 266          if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
 267             goto LBL_B;
 268          }
 269          if (!res) {
 270             goto LBL_B;
 271          }
 272       }
 273    }
 274
 275    /* passed the test */
 276    *result = true;
 277 LBL_B:
 278    mp_clear(&b);
 279    return err;
 280 }
 281
 282 #endif