Mercurial > hg > Members > anatofuz > MoarVM

comparison 3rdparty/libtommath/bn_mp_div.c @ 0:2cf249471370

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
convert mercurial for git
author Takahiro SHIMIZU <anatofuz@cr.ie.u-ryukyu.ac.jp>
date Tue, 08 May 2018 16:09:12 +0900
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:2cf249471370
1 #include <tommath_private.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, tstdenis82@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) == MP_YES) {
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) == MP_YES) {
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 }
202 q.dp[(i - t) - 1] = (mp_digit) (tmp & (mp_word) (MP_MASK));
203 }
204
205 /* while (q{i-t-1} * (yt * b + y{t-1})) >
206 xi * b**2 + xi-1 * b + xi-2
207
208 do q{i-t-1} -= 1;
209 */
210 q.dp[(i - t) - 1] = (q.dp[(i - t) - 1] + 1) & MP_MASK;
211 do {
212 q.dp[(i - t) - 1] = (q.dp[(i - t) - 1] - 1) & MP_MASK;
213
214 /* find left hand */
215 mp_zero (&t1);
216 t1.dp[0] = ((t - 1) < 0) ? 0 : y.dp[t - 1];
217 t1.dp[1] = y.dp[t];
218 t1.used = 2;
219 if ((res = mp_mul_d (&t1, q.dp[(i - t) - 1], &t1)) != MP_OKAY) {
220 goto LBL_Y;
221 }
222
223 /* find right hand */
224 t2.dp[0] = ((i - 2) < 0) ? 0 : x.dp[i - 2];
225 t2.dp[1] = ((i - 1) < 0) ? 0 : x.dp[i - 1];
226 t2.dp[2] = x.dp[i];
227 t2.used = 3;
228 } while (mp_cmp_mag(&t1, &t2) == MP_GT);
229
230 /* step 3.3 x = x - q{i-t-1} * y * b**{i-t-1} */
231 if ((res = mp_mul_d (&y, q.dp[(i - t) - 1], &t1)) != MP_OKAY) {
232 goto LBL_Y;
233 }
234
235 if ((res = mp_lshd (&t1, (i - t) - 1)) != MP_OKAY) {
236 goto LBL_Y;
237 }
238
239 if ((res = mp_sub (&x, &t1, &x)) != MP_OKAY) {
240 goto LBL_Y;
241 }
242
243 /* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */
244 if (x.sign == MP_NEG) {
245 if ((res = mp_copy (&y, &t1)) != MP_OKAY) {
246 goto LBL_Y;
247 }
248 if ((res = mp_lshd (&t1, (i - t) - 1)) != MP_OKAY) {
249 goto LBL_Y;
250 }
251 if ((res = mp_add (&x, &t1, &x)) != MP_OKAY) {
252 goto LBL_Y;
253 }
254
255 q.dp[(i - t) - 1] = (q.dp[(i - t) - 1] - 1UL) & MP_MASK;
256 }
257 }
258
259 /* now q is the quotient and x is the remainder
260 * [which we have to normalize]
261 */
262
263 /* get sign before writing to c */
264 x.sign = (x.used == 0) ? MP_ZPOS : a->sign;
265
266 if (c != NULL) {
267 mp_clamp (&q);
268 mp_exch (&q, c);
269 c->sign = neg;
270 }
271
272 if (d != NULL) {
273 if ((res = mp_div_2d (&x, norm, &x, NULL)) != MP_OKAY) {
274 goto LBL_Y;
275 }
276 mp_exch (&x, d);
277 }
278
279 res = MP_OKAY;
280
281 LBL_Y:mp_clear (&y);
282 LBL_X:mp_clear (&x);
283 LBL_T2:mp_clear (&t2);
284 LBL_T1:mp_clear (&t1);
285 LBL_Q:mp_clear (&q);
286 return res;
287 }
288
289 #endif
290
291 #endif
292
293 /* $Source$ */
294 /* $Revision$ */
295 /* $Date$ */