Line data    Source code

       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@gmail.com, http://libtom.org
      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         160 : 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         160 :   if (mp_iseven (b) == 1) {
      31           0 :     return MP_VAL;
      32             :   }
      33             : 
      34             :   /* init all our temps */
      35         160 :   if ((res = mp_init_multi(&x, &y, &u, &v, &B, &D, NULL)) != MP_OKAY) {
      36           0 :      return res;
      37             :   }
      38             : 
      39             :   /* x == modulus, y == value to invert */
      40         160 :   if ((res = mp_copy (b, &x)) != MP_OKAY) {
      41           0 :     goto LBL_ERR;
      42             :   }
      43             : 
      44             :   /* we need y = |a| */
      45         160 :   if ((res = mp_mod (a, b, &y)) != MP_OKAY) {
      46           0 :     goto LBL_ERR;
      47             :   }
      48             : 
      49             :   /* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
      50         160 :   if ((res = mp_copy (&x, &u)) != MP_OKAY) {
      51           0 :     goto LBL_ERR;
      52             :   }
      53         160 :   if ((res = mp_copy (&y, &v)) != MP_OKAY) {
      54           0 :     goto LBL_ERR;
      55             :   }
      56         160 :   mp_set (&D, 1);
      57             : 
      58      346863 : top:
      59             :   /* 4.  while u is even do */
      60     1017435 :   while (mp_iseven (&u) == 1) {
      61             :     /* 4.1 u = u/2 */
      62      346861 :     if ((res = mp_div_2 (&u, &u)) != MP_OKAY) {
      63           0 :       goto LBL_ERR;
      64             :     }
      65             :     /* 4.2 if B is odd then */
      66      346861 :     if (mp_isodd (&B) == 1) {
      67      173135 :       if ((res = mp_sub (&B, &x, &B)) != MP_OKAY) {
      68           0 :         goto LBL_ERR;
      69             :       }
      70             :     }
      71             :     /* B = B/2 */
      72      346861 :     if ((res = mp_div_2 (&B, &B)) != MP_OKAY) {
      73           0 :       goto LBL_ERR;
      74             :     }
      75             :   }
      76             : 
      77             :   /* 5.  while v is even do */
      78     1017745 :   while (mp_iseven (&v) == 1) {
      79             :     /* 5.1 v = v/2 */
      80      347171 :     if ((res = mp_div_2 (&v, &v)) != MP_OKAY) {
      81           0 :       goto LBL_ERR;
      82             :     }
      83             :     /* 5.2 if D is odd then */
      84      347171 :     if (mp_isodd (&D) == 1) {
      85             :       /* D = (D-x)/2 */
      86      173750 :       if ((res = mp_sub (&D, &x, &D)) != MP_OKAY) {
      87           0 :         goto LBL_ERR;
      88             :       }
      89             :     }
      90             :     /* D = D/2 */
      91      347171 :     if ((res = mp_div_2 (&D, &D)) != MP_OKAY) {
      92           0 :       goto LBL_ERR;
      93             :     }
      94             :   }
      95             : 
      96             :   /* 6.  if u >= v then */
      97      346871 :   if (mp_cmp (&u, &v) != MP_LT) {
      98             :     /* u = u - v, B = B - D */
      99      173721 :     if ((res = mp_sub (&u, &v, &u)) != MP_OKAY) {
     100           0 :       goto LBL_ERR;
     101             :     }
     102             : 
     103      173721 :     if ((res = mp_sub (&B, &D, &B)) != MP_OKAY) {
     104           0 :       goto LBL_ERR;
     105             :     }
     106             :   } else {
     107             :     /* v - v - u, D = D - B */
     108      173150 :     if ((res = mp_sub (&v, &u, &v)) != MP_OKAY) {
     109           0 :       goto LBL_ERR;
     110             :     }
     111             : 
     112      173150 :     if ((res = mp_sub (&D, &B, &D)) != MP_OKAY) {
     113           0 :       goto LBL_ERR;
     114             :     }
     115             :   }
     116             : 
     117             :   /* if not zero goto step 4 */
     118      346871 :   if (mp_iszero (&u) == 0) {
     119      323551 :     goto top;
     120             :   }
     121             : 
     122             :   /* now a = C, b = D, gcd == g*v */
     123             : 
     124             :   /* if v != 1 then there is no inverse */
     125         160 :   if (mp_cmp_d (&v, 1) != MP_EQ) {
     126           0 :     res = MP_VAL;
     127           0 :     goto LBL_ERR;
     128             :   }
     129             : 
     130             :   /* b is now the inverse */
     131         160 :   neg = a->sign;
     132         451 :   while (D.sign == MP_NEG) {
     133         131 :     if ((res = mp_add (&D, b, &D)) != MP_OKAY) {
     134           0 :       goto LBL_ERR;
     135             :     }
     136             :   }
     137         160 :   mp_exch (&D, c);
     138         160 :   c->sign = neg;
     139         160 :   res = MP_OKAY;
     140             : 
     141         160 : LBL_ERR:mp_clear_multi (&x, &y, &u, &v, &B, &D, NULL);
     142         160 :   return res;
     143             : }
     144             : #endif
     145             : 
     146             : /* $Source: /cvs/libtom/libtommath/bn_fast_mp_invmod.c,v $ */
     147             : /* $Revision: 1.4 $ */
     148             : /* $Date: 2006/12/28 01:25:13 $ */