Home | History | Annotate | Download | only in libtommath 
      1 #include <tommath.h>
      2 #ifdef BN_MP_PRIME_MILLER_RABIN_C
      3 /* LibTomMath, multiple-precision integer library -- Tom St Denis
      4  *
      5  * LibTomMath is a library that provides multiple-precision
      6  * integer arithmetic as well as number theoretic functionality.
      7  *
      8  * The library was designed directly after the MPI library by
      9  * Michael Fromberger but has been written from scratch with
     10  * additional optimizations in place.
     11  *
     12  * The library is free for all purposes without any express
     13  * guarantee it works.
     14  *
     15  * Tom St Denis, tomstdenis (at) gmail.com, http://math.libtomcrypt.com
     16  */
     17 
     18 /* Miller-Rabin test of "a" to the base of "b" as described in
     19  * HAC pp. 139 Algorithm 4.24
     20  *
     21  * Sets result to 0 if definitely composite or 1 if probably prime.
     22  * Randomly the chance of error is no more than 1/4 and often
     23  * very much lower.
     24  */
     25 int mp_prime_miller_rabin (mp_int * a, mp_int * b, int *result)
     26 {
     27   mp_int  n1, y, r;
     28   int     s, j, err;
     29 
     30   /* default */
     31   *result = MP_NO;
     32 
     33   /* ensure b > 1 */
     34   if (mp_cmp_d(b, 1) != MP_GT) {
     35      return MP_VAL;
     36   }
     37 
     38   /* get n1 = a - 1 */
     39   if ((err = mp_init_copy (&n1, a)) != MP_OKAY) {
     40     return err;
     41   }
     42   if ((err = mp_sub_d (&n1, 1, &n1)) != MP_OKAY) {
     43     goto LBL_N1;
     44   }
     45 
     46   /* set 2**s * r = n1 */
     47   if ((err = mp_init_copy (&r, &n1)) != MP_OKAY) {
     48     goto LBL_N1;
     49   }
     50 
     51   /* count the number of least significant bits
     52    * which are zero
     53    */
     54   s = mp_cnt_lsb(&r);
     55 
     56   /* now divide n - 1 by 2**s */
     57   if ((err = mp_div_2d (&r, s, &r, NULL)) != MP_OKAY) {
     58     goto LBL_R;
     59   }
     60 
     61   /* compute y = b**r mod a */
     62   if ((err = mp_init (&y)) != MP_OKAY) {
     63     goto LBL_R;
     64   }
     65   if ((err = mp_exptmod (b, &r, a, &y)) != MP_OKAY) {
     66     goto LBL_Y;
     67   }
     68 
     69   /* if y != 1 and y != n1 do */
     70   if (mp_cmp_d (&y, 1) != MP_EQ && mp_cmp (&y, &n1) != MP_EQ) {
     71     j = 1;
     72     /* while j <= s-1 and y != n1 */
     73     while ((j <= (s - 1)) && mp_cmp (&y, &n1) != MP_EQ) {
     74       if ((err = mp_sqrmod (&y, a, &y)) != MP_OKAY) {
     75          goto LBL_Y;
     76       }
     77 
     78       /* if y == 1 then composite */
     79       if (mp_cmp_d (&y, 1) == MP_EQ) {
     80          goto LBL_Y;
     81       }
     82 
     83       ++j;
     84     }
     85 
     86     /* if y != n1 then composite */
     87     if (mp_cmp (&y, &n1) != MP_EQ) {
     88       goto LBL_Y;
     89     }
     90   }
     91 
     92   /* probably prime now */
     93   *result = MP_YES;
     94 LBL_Y:mp_clear (&y);
     95 LBL_R:mp_clear (&r);
     96 LBL_N1:mp_clear (&n1);
     97   return err;
     98 }
     99 #endif
    100 
    101 /* $Source: /cvs/libtom/libtommath/bn_mp_prime_miller_rabin.c,v $ */
    102 /* $Revision: 1.3 $ */
    103 /* $Date: 2006/03/31 14:18:44 $ */
    104