fp_prime_miller_rabin.c

00001 /* TomsFastMath, a fast ISO C bignum library.
00002  * 
00003  * This project is meant to fill in where LibTomMath
00004  * falls short.  That is speed ;-)
00005  *
00006  * This project is public domain and free for all purposes.
00007  * 
00008  * Tom St Denis, tomstdenis@gmail.com
00009  */
00010 #include <tfm.h>
00011 
00012 /* Miller-Rabin test of "a" to the base of "b" as described in 
00013  * HAC pp. 139 Algorithm 4.24
00014  *
00015  * Sets result to 0 if definitely composite or 1 if probably prime.
00016  * Randomly the chance of error is no more than 1/4 and often 
00017  * very much lower.
00018  */
00019 void fp_prime_miller_rabin (fp_int * a, fp_int * b, int *result)
00020 {
00021   fp_int  n1, y, r;
00022   int     s, j;
00023 
00024   /* default */
00025   *result = FP_NO;
00026 
00027   /* ensure b > 1 */
00028   if (fp_cmp_d(b, 1) != FP_GT) {
00029      return;
00030   }     
00031 
00032   /* get n1 = a - 1 */
00033   fp_init_copy(&n1, a);
00034   fp_sub_d(&n1, 1, &n1);
00035 
00036   /* set 2**s * r = n1 */
00037   fp_init_copy(&r, &n1);
00038 
00039   /* count the number of least significant bits
00040    * which are zero
00041    */
00042   s = fp_cnt_lsb(&r);
00043 
00044   /* now divide n - 1 by 2**s */
00045   fp_div_2d (&r, s, &r, NULL);
00046 
00047   /* compute y = b**r mod a */
00048   fp_init(&y);
00049   fp_exptmod(b, &r, a, &y);
00050 
00051   /* if y != 1 and y != n1 do */
00052   if (fp_cmp_d (&y, 1) != FP_EQ && fp_cmp (&y, &n1) != FP_EQ) {
00053     j = 1;
00054     /* while j <= s-1 and y != n1 */
00055     while ((j <= (s - 1)) && fp_cmp (&y, &n1) != FP_EQ) {
00056       fp_sqrmod (&y, a, &y);
00057 
00058       /* if y == 1 then composite */
00059       if (fp_cmp_d (&y, 1) == FP_EQ) {
00060          return;
00061       }
00062       ++j;
00063     }
00064 
00065     /* if y != n1 then composite */
00066     if (fp_cmp (&y, &n1) != FP_EQ) {
00067        return;
00068     }
00069   }
00070 
00071   /* probably prime now */
00072   *result = FP_YES;
00073 }
00074 
00075 /* $Source: /cvs/libtom/tomsfastmath/src/numtheory/fp_prime_miller_rabin.c,v $ */
00076 /* $Revision: 1.1 $ */
00077 /* $Date: 2007/01/24 21:25:19 $ */