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