|
282
|
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 * |
|
386
|
15 * Tom St Denis, [email protected], http://math.libtomcrypt.com |
|
282
|
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 |
|
386
|
289 |
|
|
290 /* $Source: /cvs/libtom/libtommath/bn_mp_div.c,v $ */ |
|
|
291 /* $Revision: 1.3 $ */ |
|
|
292 /* $Date: 2006/03/31 14:18:44 $ */ |
