1#include <tommath.h>
2#ifdef BN_MP_DIV_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#ifdef BN_MP_DIV_SMALL
19
20/* slower bit-bang division... also smaller */
21int mp_div(mp_int * a, mp_int * b, mp_int * c, mp_int * d)
22{
23   mp_int ta, tb, tq, q;
24   int    res, n, n2;
25
26  /* is divisor zero ? */
27  if (mp_iszero (b) == 1) {
28    return MP_VAL;
29  }
30
31  /* if a < b then q=0, r = a */
32  if (mp_cmp_mag (a, b) == MP_LT) {
33    if (d != NULL) {
34      res = mp_copy (a, d);
35    } else {
36      res = MP_OKAY;
37    }
38    if (c != NULL) {
39      mp_zero (c);
40    }
41    return res;
42  }
43
44  /* init our temps */
45  if ((res = mp_init_multi(&ta, &tb, &tq, &q, NULL) != MP_OKAY)) {
46     return res;
47  }
48
49
50  mp_set(&tq, 1);
51  n = mp_count_bits(a) - mp_count_bits(b);
52  if (((res = mp_abs(a, &ta)) != MP_OKAY) ||
53      ((res = mp_abs(b, &tb)) != MP_OKAY) ||
54      ((res = mp_mul_2d(&tb, n, &tb)) != MP_OKAY) ||
55      ((res = mp_mul_2d(&tq, n, &tq)) != MP_OKAY)) {
56      goto LBL_ERR;
57  }
58
59  while (n-- >= 0) {
60     if (mp_cmp(&tb, &ta) != MP_GT) {
61        if (((res = mp_sub(&ta, &tb, &ta)) != MP_OKAY) ||
62            ((res = mp_add(&q, &tq, &q)) != MP_OKAY)) {
63           goto LBL_ERR;
64        }
65     }
66     if (((res = mp_div_2d(&tb, 1, &tb, NULL)) != MP_OKAY) ||
67         ((res = mp_div_2d(&tq, 1, &tq, NULL)) != MP_OKAY)) {
68           goto LBL_ERR;
69     }
70  }
71
72  /* now q == quotient and ta == remainder */
73  n  = a->sign;
74  n2 = (a->sign == b->sign ? MP_ZPOS : MP_NEG);
75  if (c != NULL) {
76     mp_exch(c, &q);
77     c->sign  = (mp_iszero(c) == MP_YES) ? MP_ZPOS : n2;
78  }
79  if (d != NULL) {
80     mp_exch(d, &ta);
81     d->sign = (mp_iszero(d) == MP_YES) ? MP_ZPOS : n;
82  }
83LBL_ERR:
84   mp_clear_multi(&ta, &tb, &tq, &q, NULL);
85   return res;
86}
87
88#else
89
90/* integer signed division.
91 * c*b + d == a [e.g. a/b, c=quotient, d=remainder]
92 * HAC pp.598 Algorithm 14.20
93 *
94 * Note that the description in HAC is horribly
95 * incomplete.  For example, it doesn't consider
96 * the case where digits are removed from 'x' in
97 * the inner loop.  It also doesn't consider the
98 * case that y has fewer than three digits, etc..
99 *
100 * The overall algorithm is as described as
101 * 14.20 from HAC but fixed to treat these cases.
102*/
103int mp_div (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
104{
105  mp_int  q, x, y, t1, t2;
106  int     res, n, t, i, norm, neg;
107
108  /* is divisor zero ? */
109  if (mp_iszero (b) == 1) {
110    return MP_VAL;
111  }
112
113  /* if a < b then q=0, r = a */
114  if (mp_cmp_mag (a, b) == MP_LT) {
115    if (d != NULL) {
116      res = mp_copy (a, d);
117    } else {
118      res = MP_OKAY;
119    }
120    if (c != NULL) {
121      mp_zero (c);
122    }
123    return res;
124  }
125
126  if ((res = mp_init_size (&q, a->used + 2)) != MP_OKAY) {
127    return res;
128  }
129  q.used = a->used + 2;
130
131  if ((res = mp_init (&t1)) != MP_OKAY) {
132    goto LBL_Q;
133  }
134
135  if ((res = mp_init (&t2)) != MP_OKAY) {
136    goto LBL_T1;
137  }
138
139  if ((res = mp_init_copy (&x, a)) != MP_OKAY) {
140    goto LBL_T2;
141  }
142
143  if ((res = mp_init_copy (&y, b)) != MP_OKAY) {
144    goto LBL_X;
145  }
146
147  /* fix the sign */
148  neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;
149  x.sign = y.sign = MP_ZPOS;
150
151  /* normalize both x and y, ensure that y >= b/2, [b == 2**DIGIT_BIT] */
152  norm = mp_count_bits(&y) % DIGIT_BIT;
153  if (norm < (int)(DIGIT_BIT-1)) {
154     norm = (DIGIT_BIT-1) - norm;
155     if ((res = mp_mul_2d (&x, norm, &x)) != MP_OKAY) {
156       goto LBL_Y;
157     }
158     if ((res = mp_mul_2d (&y, norm, &y)) != MP_OKAY) {
159       goto LBL_Y;
160     }
161  } else {
162     norm = 0;
163  }
164
165  /* note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4 */
166  n = x.used - 1;
167  t = y.used - 1;
168
169  /* while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} } */
170  if ((res = mp_lshd (&y, n - t)) != MP_OKAY) { /* y = y*b**{n-t} */
171    goto LBL_Y;
172  }
173
174  while (mp_cmp (&x, &y) != MP_LT) {
175    ++(q.dp[n - t]);
176    if ((res = mp_sub (&x, &y, &x)) != MP_OKAY) {
177      goto LBL_Y;
178    }
179  }
180
181  /* reset y by shifting it back down */
182  mp_rshd (&y, n - t);
183
184  /* step 3. for i from n down to (t + 1) */
185  for (i = n; i >= (t + 1); i--) {
186    if (i > x.used) {
187      continue;
188    }
189
190    /* step 3.1 if xi == yt then set q{i-t-1} to b-1,
191     * otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */
192    if (x.dp[i] == y.dp[t]) {
193      q.dp[i - t - 1] = ((((mp_digit)1) << DIGIT_BIT) - 1);
194    } else {
195      mp_word tmp;
196      tmp = ((mp_word) x.dp[i]) << ((mp_word) DIGIT_BIT);
197      tmp |= ((mp_word) x.dp[i - 1]);
198      tmp /= ((mp_word) y.dp[t]);
199      if (tmp > (mp_word) MP_MASK)
200        tmp = MP_MASK;
201      q.dp[i - t - 1] = (mp_digit) (tmp & (mp_word) (MP_MASK));
202    }
203
204    /* while (q{i-t-1} * (yt * b + y{t-1})) >
205             xi * b**2 + xi-1 * b + xi-2
206
207       do q{i-t-1} -= 1;
208    */
209    q.dp[i - t - 1] = (q.dp[i - t - 1] + 1) & MP_MASK;
210    do {
211      q.dp[i - t - 1] = (q.dp[i - t - 1] - 1) & MP_MASK;
212
213      /* find left hand */
214      mp_zero (&t1);
215      t1.dp[0] = (t - 1 < 0) ? 0 : y.dp[t - 1];
216      t1.dp[1] = y.dp[t];
217      t1.used = 2;
218      if ((res = mp_mul_d (&t1, q.dp[i - t - 1], &t1)) != MP_OKAY) {
219        goto LBL_Y;
220      }
221
222      /* find right hand */
223      t2.dp[0] = (i - 2 < 0) ? 0 : x.dp[i - 2];
224      t2.dp[1] = (i - 1 < 0) ? 0 : x.dp[i - 1];
225      t2.dp[2] = x.dp[i];
226      t2.used = 3;
227    } while (mp_cmp_mag(&t1, &t2) == MP_GT);
228
229    /* step 3.3 x = x - q{i-t-1} * y * b**{i-t-1} */
230    if ((res = mp_mul_d (&y, q.dp[i - t - 1], &t1)) != MP_OKAY) {
231      goto LBL_Y;
232    }
233
234    if ((res = mp_lshd (&t1, i - t - 1)) != MP_OKAY) {
235      goto LBL_Y;
236    }
237
238    if ((res = mp_sub (&x, &t1, &x)) != MP_OKAY) {
239      goto LBL_Y;
240    }
241
242    /* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */
243    if (x.sign == MP_NEG) {
244      if ((res = mp_copy (&y, &t1)) != MP_OKAY) {
245        goto LBL_Y;
246      }
247      if ((res = mp_lshd (&t1, i - t - 1)) != MP_OKAY) {
248        goto LBL_Y;
249      }
250      if ((res = mp_add (&x, &t1, &x)) != MP_OKAY) {
251        goto LBL_Y;
252      }
253
254      q.dp[i - t - 1] = (q.dp[i - t - 1] - 1UL) & MP_MASK;
255    }
256  }
257
258  /* now q is the quotient and x is the remainder
259   * [which we have to normalize]
260   */
261
262  /* get sign before writing to c */
263  x.sign = x.used == 0 ? MP_ZPOS : a->sign;
264
265  if (c != NULL) {
266    mp_clamp (&q);
267    mp_exch (&q, c);
268    c->sign = neg;
269  }
270
271  if (d != NULL) {
272    mp_div_2d (&x, norm, &x, NULL);
273    mp_exch (&x, d);
274  }
275
276  res = MP_OKAY;
277
278LBL_Y:mp_clear (&y);
279LBL_X:mp_clear (&x);
280LBL_T2:mp_clear (&t2);
281LBL_T1:mp_clear (&t1);
282LBL_Q:mp_clear (&q);
283  return res;
284}
285
286#endif
287
288#endif
289
290/* $Source: /cvsroot/tcl/libtommath/bn_mp_div.c,v $ */
291/* $Revision: 1.4 $ */
292/* $Date: 2006/12/01 19:45:38 $ */
293