mp_prime_miller_rabin.c

   1 #include "tommath_private.h"
   2 #ifdef MP_PRIME_MILLER_RABIN_C
   3 /* LibTomMath, multiple-precision integer library -- Tom St Denis */
   4 /* SPDX-License-Identifier: Unlicense */
   5
   6 /* Miller-Rabin test of "a" to the base of "b" as described in
   7  * HAC pp. 139 Algorithm 4.24
   8  *
   9  * Sets result to 0 if definitely composite or 1 if probably prime.
  10  * Randomly the chance of error is no more than 1/4 and often
  11  * very much lower.
  12  */
  13 mp_err mp_prime_miller_rabin(const mp_int *a, const mp_int *b, bool *result)
  14 {
  15    mp_int  n1, y, r;
  16    mp_err  err;
  17    int     s, j;
  18
  19    /* ensure b > 1 */
  20    if (mp_cmp_d(b, 1uL) != MP_GT) {
  21       return MP_VAL;
  22    }
  23
  24    /* get n1 = a - 1 */
  25    if ((err = mp_init_copy(&n1, a)) != MP_OKAY) {
  26       return err;
  27    }
  28    if ((err = mp_sub_d(&n1, 1uL, &n1)) != MP_OKAY) {
  29       goto LBL_ERR1;
  30    }
  31
  32    /* set 2**s * r = n1 */
  33    if ((err = mp_init_copy(&r, &n1)) != MP_OKAY) {
  34       goto LBL_ERR1;
  35    }
  36
  37    /* count the number of least significant bits
  38     * which are zero
  39     */
  40    s = mp_cnt_lsb(&r);
  41
  42    /* now divide n - 1 by 2**s */
  43    if ((err = mp_div_2d(&r, s, &r, NULL)) != MP_OKAY) {
  44       goto LBL_ERR2;
  45    }
  46
  47    /* compute y = b**r mod a */
  48    if ((err = mp_init(&y)) != MP_OKAY) {
  49       goto LBL_ERR2;
  50    }
  51    if ((err = mp_exptmod(b, &r, a, &y)) != MP_OKAY) {
  52       goto LBL_END;
  53    }
  54
  55    /* if y != 1 and y != n1 do */
  56    if ((mp_cmp_d(&y, 1uL) != MP_EQ) && (mp_cmp(&y, &n1) != MP_EQ)) {
  57       j = 1;
  58       /* while j <= s-1 and y != n1 */
  59       while ((j <= (s - 1)) && (mp_cmp(&y, &n1) != MP_EQ)) {
  60          if ((err = mp_sqrmod(&y, a, &y)) != MP_OKAY) {
  61             goto LBL_END;
  62          }
  63
  64          /* if y == 1 then composite */
  65          if (mp_cmp_d(&y, 1uL) == MP_EQ) {
  66             *result = false;
  67             goto LBL_END;
  68          }
  69
  70          ++j;
  71       }
  72
  73       /* if y != n1 then composite */
  74       if (mp_cmp(&y, &n1) != MP_EQ) {
  75          *result = false;
  76          goto LBL_END;
  77       }
  78    }
  79
  80    /* probably prime now */
  81    *result = true;
  82
  83 LBL_END:
  84    mp_clear(&y);
  85 LBL_ERR2:
  86    mp_clear(&r);
  87 LBL_ERR1:
  88    mp_clear(&n1);
  89    return err;
  90 }
  91 #endif