Mercurial > dropbear

comparison libtommath/bn_mp_div.c @ 1739:13d834efc376 fuzz

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
merge from main
author Matt Johnston <matt@ucc.asn.au>
date Thu, 15 Oct 2020 19:55:15 +0800
parents 1051e4eea25a
children
comparison
equal deleted inserted replaced
1562:768ebf737aa0 1739:13d834efc376
1 #include <tommath_private.h> 1 #include "tommath_private.h"
2 #ifdef BN_MP_DIV_C 2 #ifdef BN_MP_DIV_C
3 /* LibTomMath, multiple-precision integer library -- Tom St Denis 3 /* LibTomMath, multiple-precision integer library -- Tom St Denis */
4 * 4 /* SPDX-License-Identifier: Unlicense */
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://libtom.org
16 */
17 5
18 #ifdef BN_MP_DIV_SMALL 6 #ifdef BN_MP_DIV_SMALL
19 7
20 /* slower bit-bang division... also smaller */ 8 /* slower bit-bang division... also smaller */
21 int mp_div(mp_int * a, mp_int * b, mp_int * c, mp_int * d) 9 mp_err mp_div(const mp_int *a, const mp_int *b, mp_int *c, mp_int *d)
22 { 10 {
23 mp_int ta, tb, tq, q; 11 mp_int ta, tb, tq, q;
24 int res, n, n2; 12 int n, n2;
25 13 mp_err err;
26 /* is divisor zero ? */ 14
27 if (mp_iszero (b) == MP_YES) { 15 /* is divisor zero ? */
28 return MP_VAL; 16 if (MP_IS_ZERO(b)) {
29 } 17 return MP_VAL;
30 18 }
31 /* if a < b then q=0, r = a */ 19
32 if (mp_cmp_mag (a, b) == MP_LT) { 20 /* if a < b then q=0, r = a */
33 if (d != NULL) { 21 if (mp_cmp_mag(a, b) == MP_LT) {
34 res = mp_copy (a, d); 22 if (d != NULL) {
35 } else { 23 err = mp_copy(a, d);
36 res = MP_OKAY; 24 } else {
37 } 25 err = MP_OKAY;
38 if (c != NULL) { 26 }
39 mp_zero (c); 27 if (c != NULL) {
40 } 28 mp_zero(c);
41 return res; 29 }
42 } 30 return err;
43 31 }
44 /* init our temps */ 32
45 if ((res = mp_init_multi(&ta, &tb, &tq, &q, NULL)) != MP_OKAY) { 33 /* init our temps */
46 return res; 34 if ((err = mp_init_multi(&ta, &tb, &tq, &q, NULL)) != MP_OKAY) {
47 } 35 return err;
48 36 }
49 37
50 mp_set(&tq, 1); 38
51 n = mp_count_bits(a) - mp_count_bits(b); 39 mp_set(&tq, 1uL);
52 if (((res = mp_abs(a, &ta)) != MP_OKAY) || 40 n = mp_count_bits(a) - mp_count_bits(b);
53 ((res = mp_abs(b, &tb)) != MP_OKAY) || 41 if ((err = mp_abs(a, &ta)) != MP_OKAY) goto LBL_ERR;
54 ((res = mp_mul_2d(&tb, n, &tb)) != MP_OKAY) || 42 if ((err = mp_abs(b, &tb)) != MP_OKAY) goto LBL_ERR;
55 ((res = mp_mul_2d(&tq, n, &tq)) != MP_OKAY)) { 43 if ((err = mp_mul_2d(&tb, n, &tb)) != MP_OKAY) goto LBL_ERR;
56 goto LBL_ERR; 44 if ((err = mp_mul_2d(&tq, n, &tq)) != MP_OKAY) goto LBL_ERR;
57 } 45
58 46 while (n-- >= 0) {
59 while (n-- >= 0) { 47 if (mp_cmp(&tb, &ta) != MP_GT) {
60 if (mp_cmp(&tb, &ta) != MP_GT) { 48 if ((err = mp_sub(&ta, &tb, &ta)) != MP_OKAY) goto LBL_ERR;
61 if (((res = mp_sub(&ta, &tb, &ta)) != MP_OKAY) || 49 if ((err = mp_add(&q, &tq, &q)) != MP_OKAY) goto LBL_ERR;
62 ((res = mp_add(&q, &tq, &q)) != MP_OKAY)) { 50 }
63 goto LBL_ERR; 51 if ((err = mp_div_2d(&tb, 1, &tb, NULL)) != MP_OKAY) goto LBL_ERR;
64 } 52 if ((err = mp_div_2d(&tq, 1, &tq, NULL)) != MP_OKAY) goto LBL_ERR;
65 } 53 }
66 if (((res = mp_div_2d(&tb, 1, &tb, NULL)) != MP_OKAY) || 54
67 ((res = mp_div_2d(&tq, 1, &tq, NULL)) != MP_OKAY)) { 55 /* now q == quotient and ta == remainder */
68 goto LBL_ERR; 56 n = a->sign;
69 } 57 n2 = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;
70 } 58 if (c != NULL) {
71 59 mp_exch(c, &q);
72 /* now q == quotient and ta == remainder */ 60 c->sign = MP_IS_ZERO(c) ? MP_ZPOS : n2;
73 n = a->sign; 61 }
74 n2 = (a->sign == b->sign) ? MP_ZPOS : MP_NEG; 62 if (d != NULL) {
75 if (c != NULL) { 63 mp_exch(d, &ta);
76 mp_exch(c, &q); 64 d->sign = MP_IS_ZERO(d) ? MP_ZPOS : n;
77 c->sign = (mp_iszero(c) == MP_YES) ? MP_ZPOS : n2; 65 }
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: 66 LBL_ERR:
84 mp_clear_multi(&ta, &tb, &tq, &q, NULL); 67 mp_clear_multi(&ta, &tb, &tq, &q, NULL);
85 return res; 68 return err;
86 } 69 }
87 70
88 #else 71 #else
89 72
90 /* integer signed division. 73 /* integer signed division.
98 * case that y has fewer than three digits, etc.. 81 * case that y has fewer than three digits, etc..
99 * 82 *
100 * The overall algorithm is as described as 83 * The overall algorithm is as described as
101 * 14.20 from HAC but fixed to treat these cases. 84 * 14.20 from HAC but fixed to treat these cases.
102 */ 85 */
103 int mp_div (mp_int * a, mp_int * b, mp_int * c, mp_int * d) 86 mp_err mp_div(const mp_int *a, const mp_int *b, mp_int *c, mp_int *d)
104 { 87 {
105 mp_int q, x, y, t1, t2; 88 mp_int q, x, y, t1, t2;
106 int res, n, t, i, norm, neg; 89 int n, t, i, norm;
107 90 mp_sign neg;
108 /* is divisor zero ? */ 91 mp_err err;
109 if (mp_iszero (b) == MP_YES) { 92
110 return MP_VAL; 93 /* is divisor zero ? */
111 } 94 if (MP_IS_ZERO(b)) {
112 95 return MP_VAL;
113 /* if a < b then q=0, r = a */ 96 }
114 if (mp_cmp_mag (a, b) == MP_LT) { 97
115 if (d != NULL) { 98 /* if a < b then q=0, r = a */
116 res = mp_copy (a, d); 99 if (mp_cmp_mag(a, b) == MP_LT) {
117 } else { 100 if (d != NULL) {
118 res = MP_OKAY; 101 err = mp_copy(a, d);
119 } 102 } else {
120 if (c != NULL) { 103 err = MP_OKAY;
121 mp_zero (c); 104 }
122 } 105 if (c != NULL) {
123 return res; 106 mp_zero(c);
124 } 107 }
125 108 return err;
126 if ((res = mp_init_size (&q, a->used + 2)) != MP_OKAY) { 109 }
127 return res; 110
128 } 111 if ((err = mp_init_size(&q, a->used + 2)) != MP_OKAY) {
129 q.used = a->used + 2; 112 return err;
130 113 }
131 if ((res = mp_init (&t1)) != MP_OKAY) { 114 q.used = a->used + 2;
132 goto LBL_Q; 115
133 } 116 if ((err = mp_init(&t1)) != MP_OKAY) goto LBL_Q;
134 117
135 if ((res = mp_init (&t2)) != MP_OKAY) { 118 if ((err = mp_init(&t2)) != MP_OKAY) goto LBL_T1;
136 goto LBL_T1; 119
137 } 120 if ((err = mp_init_copy(&x, a)) != MP_OKAY) goto LBL_T2;
138 121
139 if ((res = mp_init_copy (&x, a)) != MP_OKAY) { 122 if ((err = mp_init_copy(&y, b)) != MP_OKAY) goto LBL_X;
140 goto LBL_T2; 123
141 } 124 /* fix the sign */
142 125 neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;
143 if ((res = mp_init_copy (&y, b)) != MP_OKAY) { 126 x.sign = y.sign = MP_ZPOS;
144 goto LBL_X; 127
145 } 128 /* normalize both x and y, ensure that y >= b/2, [b == 2**MP_DIGIT_BIT] */
146 129 norm = mp_count_bits(&y) % MP_DIGIT_BIT;
147 /* fix the sign */ 130 if (norm < (MP_DIGIT_BIT - 1)) {
148 neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG; 131 norm = (MP_DIGIT_BIT - 1) - norm;
149 x.sign = y.sign = MP_ZPOS; 132 if ((err = mp_mul_2d(&x, norm, &x)) != MP_OKAY) goto LBL_Y;
150 133 if ((err = mp_mul_2d(&y, norm, &y)) != MP_OKAY) goto LBL_Y;
151 /* normalize both x and y, ensure that y >= b/2, [b == 2**DIGIT_BIT] */ 134 } else {
152 norm = mp_count_bits(&y) % DIGIT_BIT; 135 norm = 0;
153 if (norm < (int)(DIGIT_BIT-1)) { 136 }
154 norm = (DIGIT_BIT-1) - norm; 137
155 if ((res = mp_mul_2d (&x, norm, &x)) != MP_OKAY) { 138 /* note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4 */
156 goto LBL_Y; 139 n = x.used - 1;
157 } 140 t = y.used - 1;
158 if ((res = mp_mul_2d (&y, norm, &y)) != MP_OKAY) { 141
159 goto LBL_Y; 142 /* while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} } */
160 } 143 /* y = y*b**{n-t} */
161 } else { 144 if ((err = mp_lshd(&y, n - t)) != MP_OKAY) goto LBL_Y;
162 norm = 0; 145
163 } 146 while (mp_cmp(&x, &y) != MP_LT) {
164 147 ++(q.dp[n - t]);
165 /* note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4 */ 148 if ((err = mp_sub(&x, &y, &x)) != MP_OKAY) goto LBL_Y;
166 n = x.used - 1; 149 }
167 t = y.used - 1; 150
168 151 /* reset y by shifting it back down */
169 /* while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} } */ 152 mp_rshd(&y, n - t);
170 if ((res = mp_lshd (&y, n - t)) != MP_OKAY) { /* y = y*b**{n-t} */ 153
171 goto LBL_Y; 154 /* step 3. for i from n down to (t + 1) */
172 } 155 for (i = n; i >= (t + 1); i--) {
173 156 if (i > x.used) {
174 while (mp_cmp (&x, &y) != MP_LT) { 157 continue;
175 ++(q.dp[n - t]); 158 }
176 if ((res = mp_sub (&x, &y, &x)) != MP_OKAY) { 159
177 goto LBL_Y; 160 /* step 3.1 if xi == yt then set q{i-t-1} to b-1,
178 } 161 * otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */
179 } 162 if (x.dp[i] == y.dp[t]) {
180 163 q.dp[(i - t) - 1] = ((mp_digit)1 << (mp_digit)MP_DIGIT_BIT) - (mp_digit)1;
181 /* reset y by shifting it back down */ 164 } else {
182 mp_rshd (&y, n - t); 165 mp_word tmp;
183 166 tmp = (mp_word)x.dp[i] << (mp_word)MP_DIGIT_BIT;
184 /* step 3. for i from n down to (t + 1) */ 167 tmp |= (mp_word)x.dp[i - 1];
185 for (i = n; i >= (t + 1); i--) { 168 tmp /= (mp_word)y.dp[t];
186 if (i > x.used) { 169 if (tmp > (mp_word)MP_MASK) {
187 continue; 170 tmp = MP_MASK;
188 } 171 }
189 172 q.dp[(i - t) - 1] = (mp_digit)(tmp & (mp_word)MP_MASK);
190 /* step 3.1 if xi == yt then set q{i-t-1} to b-1, 173 }
191 * otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */ 174
192 if (x.dp[i] == y.dp[t]) { 175 /* while (q{i-t-1} * (yt * b + y{t-1})) >
193 q.dp[(i - t) - 1] = ((((mp_digit)1) << DIGIT_BIT) - 1); 176 xi * b**2 + xi-1 * b + xi-2
194 } else { 177
195 mp_word tmp; 178 do q{i-t-1} -= 1;
196 tmp = ((mp_word) x.dp[i]) << ((mp_word) DIGIT_BIT); 179 */
197 tmp |= ((mp_word) x.dp[i - 1]); 180 q.dp[(i - t) - 1] = (q.dp[(i - t) - 1] + 1uL) & (mp_digit)MP_MASK;
198 tmp /= ((mp_word) y.dp[t]); 181 do {
199 if (tmp > (mp_word) MP_MASK) { 182 q.dp[(i - t) - 1] = (q.dp[(i - t) - 1] - 1uL) & (mp_digit)MP_MASK;
200 tmp = MP_MASK; 183
201 } 184 /* find left hand */
202 q.dp[(i - t) - 1] = (mp_digit) (tmp & (mp_word) (MP_MASK)); 185 mp_zero(&t1);
203 } 186 t1.dp[0] = ((t - 1) < 0) ? 0u : y.dp[t - 1];
204 187 t1.dp[1] = y.dp[t];
205 /* while (q{i-t-1} * (yt * b + y{t-1})) > 188 t1.used = 2;
206 xi * b**2 + xi-1 * b + xi-2 189 if ((err = mp_mul_d(&t1, q.dp[(i - t) - 1], &t1)) != MP_OKAY) goto LBL_Y;
207 190
208 do q{i-t-1} -= 1; 191 /* find right hand */
192 t2.dp[0] = ((i - 2) < 0) ? 0u : x.dp[i - 2];
193 t2.dp[1] = x.dp[i - 1]; /* i >= 1 always holds */
194 t2.dp[2] = x.dp[i];
195 t2.used = 3;
196 } while (mp_cmp_mag(&t1, &t2) == MP_GT);
197
198 /* step 3.3 x = x - q{i-t-1} * y * b**{i-t-1} */
199 if ((err = mp_mul_d(&y, q.dp[(i - t) - 1], &t1)) != MP_OKAY) goto LBL_Y;
200
201 if ((err = mp_lshd(&t1, (i - t) - 1)) != MP_OKAY) goto LBL_Y;
202
203 if ((err = mp_sub(&x, &t1, &x)) != MP_OKAY) goto LBL_Y;
204
205 /* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */
206 if (x.sign == MP_NEG) {
207 if ((err = mp_copy(&y, &t1)) != MP_OKAY) goto LBL_Y;
208 if ((err = mp_lshd(&t1, (i - t) - 1)) != MP_OKAY) goto LBL_Y;
209 if ((err = mp_add(&x, &t1, &x)) != MP_OKAY) goto LBL_Y;
210
211 q.dp[(i - t) - 1] = (q.dp[(i - t) - 1] - 1uL) & MP_MASK;
212 }
213 }
214
215 /* now q is the quotient and x is the remainder
216 * [which we have to normalize]
209 */ 217 */
210 q.dp[(i - t) - 1] = (q.dp[(i - t) - 1] + 1) & MP_MASK; 218
211 do { 219 /* get sign before writing to c */
212 q.dp[(i - t) - 1] = (q.dp[(i - t) - 1] - 1) & MP_MASK; 220 x.sign = (x.used == 0) ? MP_ZPOS : a->sign;
213 221
214 /* find left hand */ 222 if (c != NULL) {
215 mp_zero (&t1); 223 mp_clamp(&q);
216 t1.dp[0] = ((t - 1) < 0) ? 0 : y.dp[t - 1]; 224 mp_exch(&q, c);
217 t1.dp[1] = y.dp[t]; 225 c->sign = neg;
218 t1.used = 2; 226 }
219 if ((res = mp_mul_d (&t1, q.dp[(i - t) - 1], &t1)) != MP_OKAY) { 227
220 goto LBL_Y; 228 if (d != NULL) {
221 } 229 if ((err = mp_div_2d(&x, norm, &x, NULL)) != MP_OKAY) goto LBL_Y;
222 230 mp_exch(&x, d);
223 /* find right hand */ 231 }
224 t2.dp[0] = ((i - 2) < 0) ? 0 : x.dp[i - 2]; 232
225 t2.dp[1] = ((i - 1) < 0) ? 0 : x.dp[i - 1]; 233 err = MP_OKAY;
226 t2.dp[2] = x.dp[i]; 234
227 t2.used = 3; 235 LBL_Y:
228 } while (mp_cmp_mag(&t1, &t2) == MP_GT); 236 mp_clear(&y);
229 237 LBL_X:
230 /* step 3.3 x = x - q{i-t-1} * y * b**{i-t-1} */ 238 mp_clear(&x);
231 if ((res = mp_mul_d (&y, q.dp[(i - t) - 1], &t1)) != MP_OKAY) { 239 LBL_T2:
232 goto LBL_Y; 240 mp_clear(&t2);
233 } 241 LBL_T1:
234 242 mp_clear(&t1);
235 if ((res = mp_lshd (&t1, (i - t) - 1)) != MP_OKAY) { 243 LBL_Q:
236 goto LBL_Y; 244 mp_clear(&q);
237 } 245 return err;
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 } 246 }
288 247
289 #endif 248 #endif
290 249
291 #endif 250 #endif
292
293 /* ref: $Format:%D$ */
294 /* git commit: $Format:%H$ */
295 /* commit time: $Format:%ai$ */