Mercurial > dropbear

comparison libtommath/bn_mp_div.c @ 1439:8d24733026c5 coverity

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
merge
author Matt Johnston <matt@ucc.asn.au>
date Sat, 24 Jun 2017 23:33:16 +0800
parents 60fc6476e044
children 8bba51a55704
comparison
equal deleted inserted replaced
1400:238a439670f5 1439:8d24733026c5
1 #include <tommath.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 *
5 * LibTomMath is a library that provides multiple-precision 5 * LibTomMath is a library that provides multiple-precision
6 * integer arithmetic as well as number theoretic functionality. 6 * integer arithmetic as well as number theoretic functionality.
10 * additional optimizations in place. 10 * additional optimizations in place.
11 * 11 *
12 * The library is free for all purposes without any express 12 * The library is free for all purposes without any express
13 * guarantee it works. 13 * guarantee it works.
14 * 14 *
15 * Tom St Denis, [email protected], http://math.libtomcrypt.com 15 * Tom St Denis, [email protected], http://libtom.org
16 */ 16 */
17 17
18 #ifdef BN_MP_DIV_SMALL 18 #ifdef BN_MP_DIV_SMALL
19 19
20 /* slower bit-bang division... also smaller */ 20 /* slower bit-bang division... also smaller */
22 { 22 {
23 mp_int ta, tb, tq, q; 23 mp_int ta, tb, tq, q;
24 int res, n, n2; 24 int res, n, n2;
25 25
26 /* is divisor zero ? */ 26 /* is divisor zero ? */
27 if (mp_iszero (b) == 1) { 27 if (mp_iszero (b) == MP_YES) {
28 return MP_VAL; 28 return MP_VAL;
29 } 29 }
30 30
31 /* if a < b then q=0, r = a */ 31 /* if a < b then q=0, r = a */
32 if (mp_cmp_mag (a, b) == MP_LT) { 32 if (mp_cmp_mag (a, b) == MP_LT) {
38 if (c != NULL) { 38 if (c != NULL) {
39 mp_zero (c); 39 mp_zero (c);
40 } 40 }
41 return res; 41 return res;
42 } 42 }
43 43
44 /* init our temps */ 44 /* init our temps */
45 if ((res = mp_init_multi(&ta, &tb, &tq, &q, NULL) != MP_OKAY)) { 45 if ((res = mp_init_multi(&ta, &tb, &tq, &q, NULL)) != MP_OKAY) {
46 return res; 46 return res;
47 } 47 }
48 48
49 49
50 mp_set(&tq, 1); 50 mp_set(&tq, 1);
51 n = mp_count_bits(a) - mp_count_bits(b); 51 n = mp_count_bits(a) - mp_count_bits(b);
52 if (((res = mp_abs(a, &ta)) != MP_OKAY) || 52 if (((res = mp_abs(a, &ta)) != MP_OKAY) ||
53 ((res = mp_abs(b, &tb)) != MP_OKAY) || 53 ((res = mp_abs(b, &tb)) != MP_OKAY) ||
54 ((res = mp_mul_2d(&tb, n, &tb)) != MP_OKAY) || 54 ((res = mp_mul_2d(&tb, n, &tb)) != MP_OKAY) ||
55 ((res = mp_mul_2d(&tq, n, &tq)) != MP_OKAY)) { 55 ((res = mp_mul_2d(&tq, n, &tq)) != MP_OKAY)) {
56 goto LBL_ERR; 56 goto LBL_ERR;
57 } 57 }
58 58
69 } 69 }
70 } 70 }
71 71
72 /* now q == quotient and ta == remainder */ 72 /* now q == quotient and ta == remainder */
73 n = a->sign; 73 n = a->sign;
74 n2 = (a->sign == b->sign ? MP_ZPOS : MP_NEG); 74 n2 = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;
75 if (c != NULL) { 75 if (c != NULL) {
76 mp_exch(c, &q); 76 mp_exch(c, &q);
77 c->sign = (mp_iszero(c) == MP_YES) ? MP_ZPOS : n2; 77 c->sign = (mp_iszero(c) == MP_YES) ? MP_ZPOS : n2;
78 } 78 }
79 if (d != NULL) { 79 if (d != NULL) {
85 return res; 85 return res;
86 } 86 }
87 87
88 #else 88 #else
89 89
90 /* integer signed division. 90 /* integer signed division.
91 * c*b + d == a [e.g. a/b, c=quotient, d=remainder] 91 * c*b + d == a [e.g. a/b, c=quotient, d=remainder]
92 * HAC pp.598 Algorithm 14.20 92 * HAC pp.598 Algorithm 14.20
93 * 93 *
94 * Note that the description in HAC is horribly 94 * Note that the description in HAC is horribly
95 * incomplete. For example, it doesn't consider 95 * incomplete. For example, it doesn't consider
96 * the case where digits are removed from 'x' in 96 * the case where digits are removed from 'x' in
97 * the inner loop. It also doesn't consider the 97 * the inner loop. It also doesn't consider the
98 * case that y has fewer than three digits, etc.. 98 * case that y has fewer than three digits, etc..
99 * 99 *
100 * The overall algorithm is as described as 100 * The overall algorithm is as described as
101 * 14.20 from HAC but fixed to treat these cases. 101 * 14.20 from HAC but fixed to treat these cases.
102 */ 102 */
103 int mp_div (mp_int * a, mp_int * b, mp_int * c, mp_int * d) 103 int mp_div (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
104 { 104 {
105 mp_int q, x, y, t1, t2; 105 mp_int q, x, y, t1, t2;
106 int res, n, t, i, norm, neg; 106 int res, n, t, i, norm, neg;
107 107
108 /* is divisor zero ? */ 108 /* is divisor zero ? */
109 if (mp_iszero (b) == 1) { 109 if (mp_iszero (b) == MP_YES) {
110 return MP_VAL; 110 return MP_VAL;
111 } 111 }
112 112
113 /* if a < b then q=0, r = a */ 113 /* if a < b then q=0, r = a */
114 if (mp_cmp_mag (a, b) == MP_LT) { 114 if (mp_cmp_mag (a, b) == MP_LT) {
185 for (i = n; i >= (t + 1); i--) { 185 for (i = n; i >= (t + 1); i--) {
186 if (i > x.used) { 186 if (i > x.used) {
187 continue; 187 continue;
188 } 188 }
189 189
190 /* step 3.1 if xi == yt then set q{i-t-1} to b-1, 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 */ 191 * otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */
192 if (x.dp[i] == y.dp[t]) { 192 if (x.dp[i] == y.dp[t]) {
193 q.dp[i - t - 1] = ((((mp_digit)1) << DIGIT_BIT) - 1); 193 q.dp[(i - t) - 1] = ((((mp_digit)1) << DIGIT_BIT) - 1);
194 } else { 194 } else {
195 mp_word tmp; 195 mp_word tmp;
196 tmp = ((mp_word) x.dp[i]) << ((mp_word) DIGIT_BIT); 196 tmp = ((mp_word) x.dp[i]) << ((mp_word) DIGIT_BIT);
197 tmp |= ((mp_word) x.dp[i - 1]); 197 tmp |= ((mp_word) x.dp[i - 1]);
198 tmp /= ((mp_word) y.dp[t]); 198 tmp /= ((mp_word) y.dp[t]);
199 if (tmp > (mp_word) MP_MASK) 199 if (tmp > (mp_word) MP_MASK) {
200 tmp = MP_MASK; 200 tmp = MP_MASK;
201 q.dp[i - t - 1] = (mp_digit) (tmp & (mp_word) (MP_MASK)); 201 }
202 } 202 q.dp[(i - t) - 1] = (mp_digit) (tmp & (mp_word) (MP_MASK));
203 203 }
204 /* while (q{i-t-1} * (yt * b + y{t-1})) > 204
205 xi * b**2 + xi-1 * b + xi-2 205 /* while (q{i-t-1} * (yt * b + y{t-1})) >
206 206 xi * b**2 + xi-1 * b + xi-2
207 do q{i-t-1} -= 1; 207
208 do q{i-t-1} -= 1;
208 */ 209 */
209 q.dp[i - t - 1] = (q.dp[i - t - 1] + 1) & MP_MASK; 210 q.dp[(i - t) - 1] = (q.dp[(i - t) - 1] + 1) & MP_MASK;
210 do { 211 do {
211 q.dp[i - t - 1] = (q.dp[i - t - 1] - 1) & MP_MASK; 212 q.dp[(i - t) - 1] = (q.dp[(i - t) - 1] - 1) & MP_MASK;
212 213
213 /* find left hand */ 214 /* find left hand */
214 mp_zero (&t1); 215 mp_zero (&t1);
215 t1.dp[0] = (t - 1 < 0) ? 0 : y.dp[t - 1]; 216 t1.dp[0] = ((t - 1) < 0) ? 0 : y.dp[t - 1];
216 t1.dp[1] = y.dp[t]; 217 t1.dp[1] = y.dp[t];
217 t1.used = 2; 218 t1.used = 2;
218 if ((res = mp_mul_d (&t1, q.dp[i - t - 1], &t1)) != MP_OKAY) { 219 if ((res = mp_mul_d (&t1, q.dp[(i - t) - 1], &t1)) != MP_OKAY) {
219 goto LBL_Y; 220 goto LBL_Y;
220 } 221 }
221 222
222 /* find right hand */ 223 /* find right hand */
223 t2.dp[0] = (i - 2 < 0) ? 0 : x.dp[i - 2]; 224 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[1] = ((i - 1) < 0) ? 0 : x.dp[i - 1];
225 t2.dp[2] = x.dp[i]; 226 t2.dp[2] = x.dp[i];
226 t2.used = 3; 227 t2.used = 3;
227 } while (mp_cmp_mag(&t1, &t2) == MP_GT); 228 } while (mp_cmp_mag(&t1, &t2) == MP_GT);
228 229
229 /* step 3.3 x = x - q{i-t-1} * y * b**{i-t-1} */ 230 /* 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 if ((res = mp_mul_d (&y, q.dp[(i - t) - 1], &t1)) != MP_OKAY) {
231 goto LBL_Y; 232 goto LBL_Y;
232 } 233 }
233 234
234 if ((res = mp_lshd (&t1, i - t - 1)) != MP_OKAY) { 235 if ((res = mp_lshd (&t1, (i - t) - 1)) != MP_OKAY) {
235 goto LBL_Y; 236 goto LBL_Y;
236 } 237 }
237 238
238 if ((res = mp_sub (&x, &t1, &x)) != MP_OKAY) { 239 if ((res = mp_sub (&x, &t1, &x)) != MP_OKAY) {
239 goto LBL_Y; 240 goto LBL_Y;
242 /* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */ 243 /* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */
243 if (x.sign == MP_NEG) { 244 if (x.sign == MP_NEG) {
244 if ((res = mp_copy (&y, &t1)) != MP_OKAY) { 245 if ((res = mp_copy (&y, &t1)) != MP_OKAY) {
245 goto LBL_Y; 246 goto LBL_Y;
246 } 247 }
247 if ((res = mp_lshd (&t1, i - t - 1)) != MP_OKAY) { 248 if ((res = mp_lshd (&t1, (i - t) - 1)) != MP_OKAY) {
248 goto LBL_Y; 249 goto LBL_Y;
249 } 250 }
250 if ((res = mp_add (&x, &t1, &x)) != MP_OKAY) { 251 if ((res = mp_add (&x, &t1, &x)) != MP_OKAY) {
251 goto LBL_Y; 252 goto LBL_Y;
252 } 253 }
253 254
254 q.dp[i - t - 1] = (q.dp[i - t - 1] - 1UL) & MP_MASK; 255 q.dp[(i - t) - 1] = (q.dp[(i - t) - 1] - 1UL) & MP_MASK;
255 } 256 }
256 } 257 }
257 258
258 /* now q is the quotient and x is the remainder 259 /* now q is the quotient and x is the remainder
259 * [which we have to normalize] 260 * [which we have to normalize]
260 */ 261 */
261 262
262 /* get sign before writing to c */ 263 /* get sign before writing to c */
263 x.sign = x.used == 0 ? MP_ZPOS : a->sign; 264 x.sign = (x.used == 0) ? MP_ZPOS : a->sign;
264 265
265 if (c != NULL) { 266 if (c != NULL) {
266 mp_clamp (&q); 267 mp_clamp (&q);
267 mp_exch (&q, c); 268 mp_exch (&q, c);
268 c->sign = neg; 269 c->sign = neg;
269 } 270 }
270 271
271 if (d != NULL) { 272 if (d != NULL) {
272 if ((res = mp_div_2d (&x, norm, &x, NULL)) != MP_OKAY) { 273 if ((res = mp_div_2d (&x, norm, &x, NULL)) != MP_OKAY) {
273 goto LBL_Y; 274 goto LBL_Y;
274 } 275 }
275 mp_exch (&x, d); 276 mp_exch (&x, d);
276 } 277 }
277 278
278 res = MP_OKAY; 279 res = MP_OKAY;
279 280
287 288
288 #endif 289 #endif
289 290
290 #endif 291 #endif
291 292
292 /* $Source: /cvs/libtom/libtommath/bn_mp_div.c,v $ */ 293 /* $Source$ */
293 /* $Revision: 1.3 $ */ 294 /* $Revision$ */
294 /* $Date: 2006/03/31 14:18:44 $ */ 295 /* $Date$ */