Mercurial > dropbear

comparison libtommath/bn_mp_prime_miller_rabin.c @ 330:5488db2e9e4e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
merge of 332f709a4cb39cde4cedab7c3be89e05f3023067 and ca4ca78b82c5d430c69ce01bf794e8886ce81431
author Matt Johnston <matt@ucc.asn.au>
date Sat, 10 Jun 2006 16:39:40 +0000
parents eed26cff980b
children 5ff8218bcee9
comparison
equal deleted inserted replaced
329:8ed0dce45126 330:5488db2e9e4e
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, [email protected], http://math.libtomcrypt.org
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