libs/libtommath/bn_mp_invmod_slow.c

   1 #include <tommath.h>
   2 #ifdef BN_MP_INVMOD_SLOW_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://math.libtomcrypt.com
  16  */
  17
  18 /* hac 14.61, pp608 */
  19 int mp_invmod_slow (mp_int * a, mp_int * b, mp_int * c)
  20 {
  21   mp_int  x, y, u, v, A, B, C, D;
  22   int     res;
  23
  24   /* b cannot be negative */
  25   if (b->sign == MP_NEG || mp_iszero(b) == 1) {
  26     return MP_VAL;
  27   }
  28
  29   /* init temps */
  30   if ((res = mp_init_multi(&x, &y, &u, &v,
  31                            &A, &B, &C, &D, NULL)) != MP_OKAY) {
  32      return res;
  33   }
  34
  35   /* x = a, y = b */
  36   if ((res = mp_mod(a, b, &x)) != MP_OKAY) {
  37       goto LBL_ERR;
  38   }
  39   if ((res = mp_copy (b, &y)) != MP_OKAY) {
  40     goto LBL_ERR;
  41   }
  42
  43   /* 2. [modified] if x,y are both even then return an error! */
  44   if (mp_iseven (&x) == 1 && mp_iseven (&y) == 1) {
  45     res = MP_VAL;
  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 (&A, 1);
  57   mp_set (&D, 1);
  58
  59 top:
  60   /* 4.  while u is even do */
  61   while (mp_iseven (&u) == 1) {
  62     /* 4.1 u = u/2 */
  63     if ((res = mp_div_2 (&u, &u)) != MP_OKAY) {
  64       goto LBL_ERR;
  65     }
  66     /* 4.2 if A or B is odd then */
  67     if (mp_isodd (&A) == 1 || mp_isodd (&B) == 1) {
  68       /* A = (A+y)/2, B = (B-x)/2 */
  69       if ((res = mp_add (&A, &y, &A)) != MP_OKAY) {
  70          goto LBL_ERR;
  71       }
  72       if ((res = mp_sub (&B, &x, &B)) != MP_OKAY) {
  73          goto LBL_ERR;
  74       }
  75     }
  76     /* A = A/2, B = B/2 */
  77     if ((res = mp_div_2 (&A, &A)) != MP_OKAY) {
  78       goto LBL_ERR;
  79     }
  80     if ((res = mp_div_2 (&B, &B)) != MP_OKAY) {
  81       goto LBL_ERR;
  82     }
  83   }
  84
  85   /* 5.  while v is even do */
  86   while (mp_iseven (&v) == 1) {
  87     /* 5.1 v = v/2 */
  88     if ((res = mp_div_2 (&v, &v)) != MP_OKAY) {
  89       goto LBL_ERR;
  90     }
  91     /* 5.2 if C or D is odd then */
  92     if (mp_isodd (&C) == 1 || mp_isodd (&D) == 1) {
  93       /* C = (C+y)/2, D = (D-x)/2 */
  94       if ((res = mp_add (&C, &y, &C)) != MP_OKAY) {
  95          goto LBL_ERR;
  96       }
  97       if ((res = mp_sub (&D, &x, &D)) != MP_OKAY) {
  98          goto LBL_ERR;
  99       }
 100     }
 101     /* C = C/2, D = D/2 */
 102     if ((res = mp_div_2 (&C, &C)) != MP_OKAY) {
 103       goto LBL_ERR;
 104     }
 105     if ((res = mp_div_2 (&D, &D)) != MP_OKAY) {
 106       goto LBL_ERR;
 107     }
 108   }
 109
 110   /* 6.  if u >= v then */
 111   if (mp_cmp (&u, &v) != MP_LT) {
 112     /* u = u - v, A = A - C, B = B - D */
 113     if ((res = mp_sub (&u, &v, &u)) != MP_OKAY) {
 114       goto LBL_ERR;
 115     }
 116
 117     if ((res = mp_sub (&A, &C, &A)) != MP_OKAY) {
 118       goto LBL_ERR;
 119     }
 120
 121     if ((res = mp_sub (&B, &D, &B)) != MP_OKAY) {
 122       goto LBL_ERR;
 123     }
 124   } else {
 125     /* v - v - u, C = C - A, D = D - B */
 126     if ((res = mp_sub (&v, &u, &v)) != MP_OKAY) {
 127       goto LBL_ERR;
 128     }
 129
 130     if ((res = mp_sub (&C, &A, &C)) != MP_OKAY) {
 131       goto LBL_ERR;
 132     }
 133
 134     if ((res = mp_sub (&D, &B, &D)) != MP_OKAY) {
 135       goto LBL_ERR;
 136     }
 137   }
 138
 139   /* if not zero goto step 4 */
 140   if (mp_iszero (&u) == 0)
 141     goto top;
 142
 143   /* now a = C, b = D, gcd == g*v */
 144
 145   /* if v != 1 then there is no inverse */
 146   if (mp_cmp_d (&v, 1) != MP_EQ) {
 147     res = MP_VAL;
 148     goto LBL_ERR;
 149   }
 150
 151   /* if its too low */
 152   while (mp_cmp_d(&C, 0) == MP_LT) {
 153       if ((res = mp_add(&C, b, &C)) != MP_OKAY) {
 154          goto LBL_ERR;
 155       }
 156   }
 157
 158   /* too big */
 159   while (mp_cmp_mag(&C, b) != MP_LT) {
 160       if ((res = mp_sub(&C, b, &C)) != MP_OKAY) {
 161          goto LBL_ERR;
 162       }
 163   }
 164
 165   /* C is now the inverse */
 166   mp_exch (&C, c);
 167   res = MP_OKAY;
 168 LBL_ERR:mp_clear_multi (&x, &y, &u, &v, &A, &B, &C, &D, NULL);
 169   return res;
 170 }
 171 #endif
 172
 173 /* $Source: /cvs/libtom/libtommath/bn_mp_invmod_slow.c,v $ */
 174 /* $Revision: 1.3 $ */
 175 /* $Date: 2006/03/31 14:18:44 $ */