Mercurial > dropbear

comparison libtommath/bn_mp_div.c @ 389:5ff8218bcee9

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
propagate from branch 'au.asn.ucc.matt.ltm.dropbear' (head 2af95f00ebd5bb7a28b3817db1218442c935388e) to branch 'au.asn.ucc.matt.dropbear' (head ecd779509ef23a8cdf64888904fc9b31d78aa933)
author Matt Johnston <matt@ucc.asn.au>
date Thu, 11 Jan 2007 03:14:55 +0000
parents 1c7a072000e0
children 60fc6476e044
comparison
equal deleted inserted replaced
388:fb54020f78e1 389:5ff8218bcee9
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, [email protected], http://math.libtomcrypt.com
16 */
17
18 #ifdef BN_MP_DIV_SMALL
19
20 /* slower bit-bang division... also smaller */
21 int 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 }
83 LBL_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 */
103 int 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 if ((res = mp_div_2d (&x, norm, &x, NULL)) != MP_OKAY) {
273 goto LBL_Y;
274 }
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.3 $ */
294 /* $Date: 2006/03/31 14:18:44 $ */