Line data Source code
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://libtom.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 1992 : 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 1992 : if (mp_iszero (b) == 1) {
110 0 : return MP_VAL;
111 : }
112 :
113 : /* if a < b then q=0, r = a */
114 1992 : if (mp_cmp_mag (a, b) == MP_LT) {
115 358 : if (d != NULL) {
116 358 : res = mp_copy (a, d);
117 : } else {
118 0 : res = MP_OKAY;
119 : }
120 358 : if (c != NULL) {
121 0 : mp_zero (c);
122 : }
123 344 : return res;
124 : }
125 :
126 1634 : if ((res = mp_init_size (&q, a->used + 2)) != MP_OKAY) {
127 0 : return res;
128 : }
129 1634 : q.used = a->used + 2;
130 :
131 1634 : if ((res = mp_init (&t1)) != MP_OKAY) {
132 0 : goto LBL_Q;
133 : }
134 :
135 1634 : if ((res = mp_init (&t2)) != MP_OKAY) {
136 0 : goto LBL_T1;
137 : }
138 :
139 1634 : if ((res = mp_init_copy (&x, a)) != MP_OKAY) {
140 0 : goto LBL_T2;
141 : }
142 :
143 1634 : if ((res = mp_init_copy (&y, b)) != MP_OKAY) {
144 0 : goto LBL_X;
145 : }
146 :
147 : /* fix the sign */
148 1634 : neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;
149 1634 : x.sign = y.sign = MP_ZPOS;
150 :
151 : /* normalize both x and y, ensure that y >= b/2, [b == 2**DIGIT_BIT] */
152 1634 : norm = mp_count_bits(&y) % DIGIT_BIT;
153 1634 : if (norm < (int)(DIGIT_BIT-1)) {
154 1634 : norm = (DIGIT_BIT-1) - norm;
155 1634 : if ((res = mp_mul_2d (&x, norm, &x)) != MP_OKAY) {
156 0 : goto LBL_Y;
157 : }
158 1634 : if ((res = mp_mul_2d (&y, norm, &y)) != MP_OKAY) {
159 0 : goto LBL_Y;
160 : }
161 : } else {
162 0 : 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 1634 : n = x.used - 1;
167 1634 : t = y.used - 1;
168 :
169 : /* while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} } */
170 1634 : if ((res = mp_lshd (&y, n - t)) != MP_OKAY) { /* y = y*b**{n-t} */
171 0 : goto LBL_Y;
172 : }
173 :
174 3373 : while (mp_cmp (&x, &y) != MP_LT) {
175 179 : ++(q.dp[n - t]);
176 179 : if ((res = mp_sub (&x, &y, &x)) != MP_OKAY) {
177 0 : goto LBL_Y;
178 : }
179 : }
180 :
181 : /* reset y by shifting it back down */
182 1634 : mp_rshd (&y, n - t);
183 :
184 : /* step 3. for i from n down to (t + 1) */
185 72946 : for (i = n; i >= (t + 1); i--) {
186 71312 : if (i > x.used) {
187 0 : 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 71312 : if (x.dp[i] == y.dp[t]) {
193 0 : q.dp[i - t - 1] = ((((mp_digit)1) << DIGIT_BIT) - 1);
194 : } else {
195 : mp_word tmp;
196 71312 : tmp = ((mp_word) x.dp[i]) << ((mp_word) DIGIT_BIT);
197 71312 : tmp |= ((mp_word) x.dp[i - 1]);
198 71312 : tmp /= ((mp_word) y.dp[t]);
199 71312 : if (tmp > (mp_word) MP_MASK)
200 0 : tmp = MP_MASK;
201 71312 : 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 71312 : q.dp[i - t - 1] = (q.dp[i - t - 1] + 1) & MP_MASK;
210 : do {
211 117043 : q.dp[i - t - 1] = (q.dp[i - t - 1] - 1) & MP_MASK;
212 :
213 : /* find left hand */
214 117043 : mp_zero (&t1);
215 117043 : t1.dp[0] = (t - 1 < 0) ? 0 : y.dp[t - 1];
216 117043 : t1.dp[1] = y.dp[t];
217 117043 : t1.used = 2;
218 117043 : if ((res = mp_mul_d (&t1, q.dp[i - t - 1], &t1)) != MP_OKAY) {
219 0 : goto LBL_Y;
220 : }
221 :
222 : /* find right hand */
223 117043 : t2.dp[0] = (i - 2 < 0) ? 0 : x.dp[i - 2];
224 117043 : t2.dp[1] = (i - 1 < 0) ? 0 : x.dp[i - 1];
225 117043 : t2.dp[2] = x.dp[i];
226 117043 : t2.used = 3;
227 117043 : } 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 71312 : if ((res = mp_mul_d (&y, q.dp[i - t - 1], &t1)) != MP_OKAY) {
231 0 : goto LBL_Y;
232 : }
233 :
234 71312 : if ((res = mp_lshd (&t1, i - t - 1)) != MP_OKAY) {
235 0 : goto LBL_Y;
236 : }
237 :
238 71312 : if ((res = mp_sub (&x, &t1, &x)) != MP_OKAY) {
239 0 : goto LBL_Y;
240 : }
241 :
242 : /* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */
243 71312 : if (x.sign == MP_NEG) {
244 0 : if ((res = mp_copy (&y, &t1)) != MP_OKAY) {
245 0 : goto LBL_Y;
246 : }
247 0 : if ((res = mp_lshd (&t1, i - t - 1)) != MP_OKAY) {
248 0 : goto LBL_Y;
249 : }
250 0 : if ((res = mp_add (&x, &t1, &x)) != MP_OKAY) {
251 0 : goto LBL_Y;
252 : }
253 :
254 0 : 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 1634 : x.sign = x.used == 0 ? MP_ZPOS : a->sign;
264 :
265 1634 : if (c != NULL) {
266 0 : mp_clamp (&q);
267 0 : mp_exch (&q, c);
268 0 : c->sign = neg;
269 : }
270 :
271 1634 : if (d != NULL) {
272 1634 : mp_div_2d (&x, norm, &x, NULL);
273 1634 : mp_exch (&x, d);
274 : }
275 :
276 1560 : res = MP_OKAY;
277 :
278 1634 : LBL_Y:mp_clear (&y);
279 1634 : LBL_X:mp_clear (&x);
280 1634 : LBL_T2:mp_clear (&t2);
281 1634 : LBL_T1:mp_clear (&t1);
282 1634 : LBL_Q:mp_clear (&q);
283 1634 : return res;
284 : }
285 :
286 : #endif
287 :
288 : #endif
289 :
290 : /* $Source: /cvs/libtom/libtommath/bn_mp_div.c,v $ */
291 : /* $Revision: 1.4 $ */
292 : /* $Date: 2006/12/28 01:25:13 $ */
|
