1 /* $NetBSD: bn_mp_prime_miller_rabin.c,v 1.2 2017/01/28 21:31:47 christos 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 */
mp_prime_miller_rabin(mp_int * a,mp_int * b,int * result)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