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