1 #include <tommath.h> 2 #ifdef BN_FAST_MP_INVMOD_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, tomstdenis (at) gmail.com, http://math.libtomcrypt.com 16 */ 17 18 /* computes the modular inverse via binary extended euclidean algorithm, 19 * that is c = 1/a mod b 20 * 21 * Based on slow invmod except this is optimized for the case where b is 22 * odd as per HAC Note 14.64 on pp. 610 23 */ 24 int fast_mp_invmod (mp_int * a, mp_int * b, mp_int * c) 25 { 26 mp_int x, y, u, v, B, D; 27 int res, neg; 28 29 /* 2. [modified] b must be odd */ 30 if (mp_iseven (b) == 1) { 31 return MP_VAL; 32 } 33 34 /* init all our temps */ 35 if ((res = mp_init_multi(&x, &y, &u, &v, &B, &D, NULL)) != MP_OKAY) { 36 return res; 37 } 38 39 /* x == modulus, y == value to invert */ 40 if ((res = mp_copy (b, &x)) != MP_OKAY) { 41 goto LBL_ERR; 42 } 43 44 /* we need y = |a| */ 45 if ((res = mp_mod (a, b, &y)) != MP_OKAY) { 46 goto LBL_ERR; 47 } 48 49 /* 3. u=x, v=y, A=1, B=0, C=0,D=1 */ 50 if ((res = mp_copy (&x, &u)) != MP_OKAY) { 51 goto LBL_ERR; 52 } 53 if ((res = mp_copy (&y, &v)) != MP_OKAY) { 54 goto LBL_ERR; 55 } 56 mp_set (&D, 1); 57 58 top: 59 /* 4. while u is even do */ 60 while (mp_iseven (&u) == 1) { 61 /* 4.1 u = u/2 */ 62 if ((res = mp_div_2 (&u, &u)) != MP_OKAY) { 63 goto LBL_ERR; 64 } 65 /* 4.2 if B is odd then */ 66 if (mp_isodd (&B) == 1) { 67 if ((res = mp_sub (&B, &x, &B)) != MP_OKAY) { 68 goto LBL_ERR; 69 } 70 } 71 /* B = B/2 */ 72 if ((res = mp_div_2 (&B, &B)) != MP_OKAY) { 73 goto LBL_ERR; 74 } 75 } 76 77 /* 5. while v is even do */ 78 while (mp_iseven (&v) == 1) { 79 /* 5.1 v = v/2 */ 80 if ((res = mp_div_2 (&v, &v)) != MP_OKAY) { 81 goto LBL_ERR; 82 } 83 /* 5.2 if D is odd then */ 84 if (mp_isodd (&D) == 1) { 85 /* D = (D-x)/2 */ 86 if ((res = mp_sub (&D, &x, &D)) != MP_OKAY) { 87 goto LBL_ERR; 88 } 89 } 90 /* D = D/2 */ 91 if ((res = mp_div_2 (&D, &D)) != MP_OKAY) { 92 goto LBL_ERR; 93 } 94 } 95 96 /* 6. if u >= v then */ 97 if (mp_cmp (&u, &v) != MP_LT) { 98 /* u = u - v, B = B - D */ 99 if ((res = mp_sub (&u, &v, &u)) != MP_OKAY) { 100 goto LBL_ERR; 101 } 102 103 if ((res = mp_sub (&B, &D, &B)) != MP_OKAY) { 104 goto LBL_ERR; 105 } 106 } else { 107 /* v - v - u, D = D - B */ 108 if ((res = mp_sub (&v, &u, &v)) != MP_OKAY) { 109 goto LBL_ERR; 110 } 111 112 if ((res = mp_sub (&D, &B, &D)) != MP_OKAY) { 113 goto LBL_ERR; 114 } 115 } 116 117 /* if not zero goto step 4 */ 118 if (mp_iszero (&u) == 0) { 119 goto top; 120 } 121 122 /* now a = C, b = D, gcd == g*v */ 123 124 /* if v != 1 then there is no inverse */ 125 if (mp_cmp_d (&v, 1) != MP_EQ) { 126 res = MP_VAL; 127 goto LBL_ERR; 128 } 129 130 /* b is now the inverse */ 131 neg = a->sign; 132 while (D.sign == MP_NEG) { 133 if ((res = mp_add (&D, b, &D)) != MP_OKAY) { 134 goto LBL_ERR; 135 } 136 } 137 mp_exch (&D, c); 138 c->sign = neg; 139 res = MP_OKAY; 140 141 LBL_ERR:mp_clear_multi (&x, &y, &u, &v, &B, &D, NULL); 142 return res; 143 } 144 #endif 145 146 /* $Source: /cvs/libtom/libtommath/bn_fast_mp_invmod.c,v $ */ 147 /* $Revision: 1.3 $ */ 148 /* $Date: 2006/03/31 14:18:44 $ */ 149
