Mercurial > dropbear

comparison libtommath/bn_mp_div.c @ 284:eed26cff980b

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 6c790cad5a7fa866ad062cb3a0c279f7ba788583) to branch 'au.asn.ucc.matt.dropbear' (head fff0894a0399405a9410ea1c6d118f342cf2aa64)
author Matt Johnston <matt@ucc.asn.au>
date Wed, 08 Mar 2006 13:23:49 +0000
parents
children 1c7a072000e0
comparison
equal deleted inserted replaced
283:bd240aa12ba7 284:eed26cff980b
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.org
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 mp_div_2d (&x, norm, &x, NULL);
273 mp_exch (&x, d);
274 }
275
276 res = MP_OKAY;
277
278 LBL_Y:mp_clear (&y);
279 LBL_X:mp_clear (&x);
280 LBL_T2:mp_clear (&t2);
281 LBL_T1:mp_clear (&t1);
282 LBL_Q:mp_clear (&q);
283 return res;
284 }
285
286 #endif
287
288 #endif