Line data    Source code

       1             : #include "tommath_private.h"
       2             : #ifdef BN_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           0 : mp_err mp_prime_miller_rabin(const mp_int *a, const mp_int *b, mp_bool *result)
      14             : {
      15             :    mp_int  n1, y, r;
      16             :    mp_err  err;
      17             :    int     s, j;
      18             : 
      19             :    /* default */
      20           0 :    *result = MP_NO;
      21             : 
      22             :    /* ensure b > 1 */
      23           0 :    if (mp_cmp_d(b, 1uL) != MP_GT) {
      24           0 :       return MP_VAL;
      25             :    }
      26             : 
      27             :    /* get n1 = a - 1 */
      28           0 :    if ((err = mp_init_copy(&n1, a)) != MP_OKAY) {
      29           0 :       return err;
      30             :    }
      31           0 :    if ((err = mp_sub_d(&n1, 1uL, &n1)) != MP_OKAY) {
      32           0 :       goto LBL_N1;
      33             :    }
      34             : 
      35             :    /* set 2**s * r = n1 */
      36           0 :    if ((err = mp_init_copy(&r, &n1)) != MP_OKAY) {
      37           0 :       goto LBL_N1;
      38             :    }
      39             : 
      40             :    /* count the number of least significant bits
      41             :     * which are zero
      42             :     */
      43           0 :    s = mp_cnt_lsb(&r);
      44             : 
      45             :    /* now divide n - 1 by 2**s */
      46           0 :    if ((err = mp_div_2d(&r, s, &r, NULL)) != MP_OKAY) {
      47           0 :       goto LBL_R;
      48             :    }
      49             : 
      50             :    /* compute y = b**r mod a */
      51           0 :    if ((err = mp_init(&y)) != MP_OKAY) {
      52           0 :       goto LBL_R;
      53             :    }
      54           0 :    if ((err = mp_exptmod(b, &r, a, &y)) != MP_OKAY) {
      55           0 :       goto LBL_Y;
      56             :    }
      57             : 
      58             :    /* if y != 1 and y != n1 do */
      59           0 :    if ((mp_cmp_d(&y, 1uL) != MP_EQ) && (mp_cmp(&y, &n1) != MP_EQ)) {
      60           0 :       j = 1;
      61             :       /* while j <= s-1 and y != n1 */
      62           0 :       while ((j <= (s - 1)) && (mp_cmp(&y, &n1) != MP_EQ)) {
      63           0 :          if ((err = mp_sqrmod(&y, a, &y)) != MP_OKAY) {
      64           0 :             goto LBL_Y;
      65             :          }
      66             : 
      67             :          /* if y == 1 then composite */
      68           0 :          if (mp_cmp_d(&y, 1uL) == MP_EQ) {
      69           0 :             goto LBL_Y;
      70             :          }
      71             : 
      72           0 :          ++j;
      73             :       }
      74             : 
      75             :       /* if y != n1 then composite */
      76           0 :       if (mp_cmp(&y, &n1) != MP_EQ) {
      77           0 :          goto LBL_Y;
      78             :       }
      79             :    }
      80             : 
      81             :    /* probably prime now */
      82           0 :    *result = MP_YES;
      83           0 : LBL_Y:
      84           0 :    mp_clear(&y);
      85           0 : LBL_R:
      86           0 :    mp_clear(&r);
      87           0 : LBL_N1:
      88           0 :    mp_clear(&n1);
      89           0 :    return err;
      90             : }
      91             : #endif