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