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