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