diff options
Diffstat (limited to 'libtommath/bn_mp_div.c')
-rw-r--r-- | libtommath/bn_mp_div.c | 483 |
1 files changed, 246 insertions, 237 deletions
diff --git a/libtommath/bn_mp_div.c b/libtommath/bn_mp_div.c index de4ca04..44e3cb9 100644 --- a/libtommath/bn_mp_div.c +++ b/libtommath/bn_mp_div.c @@ -1,4 +1,4 @@ -#include <tommath.h> +#include "tommath_private.h" #ifdef BN_MP_DIV_C /* LibTomMath, multiple-precision integer library -- Tom St Denis * @@ -9,77 +9,74 @@ * Michael Fromberger but has been written from scratch with * additional optimizations in place. * - * The library is free for all purposes without any express - * guarantee it works. - * - * Tom St Denis, tomstdenis@gmail.com, http://math.libtomcrypt.com + * SPDX-License-Identifier: Unlicense */ #ifdef BN_MP_DIV_SMALL /* slower bit-bang division... also smaller */ -int mp_div(mp_int * a, mp_int * b, mp_int * c, mp_int * d) +int mp_div(const mp_int *a, const mp_int *b, mp_int *c, mp_int *d) { mp_int ta, tb, tq, q; int res, n, n2; - /* is divisor zero ? */ - if (mp_iszero (b) == 1) { - return MP_VAL; - } - - /* if a < b then q=0, r = a */ - if (mp_cmp_mag (a, b) == MP_LT) { - if (d != NULL) { - res = mp_copy (a, d); - } else { - res = MP_OKAY; - } - if (c != NULL) { - mp_zero (c); - } - return res; - } - - /* init our temps */ - if ((res = mp_init_multi(&ta, &tb, &tq, &q, NULL) != MP_OKAY)) { - return res; - } - - - mp_set(&tq, 1); - n = mp_count_bits(a) - mp_count_bits(b); - if (((res = mp_abs(a, &ta)) != MP_OKAY) || - ((res = mp_abs(b, &tb)) != MP_OKAY) || - ((res = mp_mul_2d(&tb, n, &tb)) != MP_OKAY) || - ((res = mp_mul_2d(&tq, n, &tq)) != MP_OKAY)) { + /* is divisor zero ? */ + if (mp_iszero(b) == MP_YES) { + return MP_VAL; + } + + /* if a < b then q=0, r = a */ + if (mp_cmp_mag(a, b) == MP_LT) { + if (d != NULL) { + res = mp_copy(a, d); + } else { + res = MP_OKAY; + } + if (c != NULL) { + mp_zero(c); + } + return res; + } + + /* init our temps */ + if ((res = mp_init_multi(&ta, &tb, &tq, &q, NULL)) != MP_OKAY) { + return res; + } + + + mp_set(&tq, 1uL); + n = mp_count_bits(a) - mp_count_bits(b); + if (((res = mp_abs(a, &ta)) != MP_OKAY) || + ((res = mp_abs(b, &tb)) != MP_OKAY) || + ((res = mp_mul_2d(&tb, n, &tb)) != MP_OKAY) || + ((res = mp_mul_2d(&tq, n, &tq)) != MP_OKAY)) { goto LBL_ERR; - } - - while (n-- >= 0) { - if (mp_cmp(&tb, &ta) != MP_GT) { - if (((res = mp_sub(&ta, &tb, &ta)) != MP_OKAY) || - ((res = mp_add(&q, &tq, &q)) != MP_OKAY)) { - goto LBL_ERR; - } - } - if (((res = mp_div_2d(&tb, 1, &tb, NULL)) != MP_OKAY) || - ((res = mp_div_2d(&tq, 1, &tq, NULL)) != MP_OKAY)) { - goto LBL_ERR; - } - } - - /* now q == quotient and ta == remainder */ - n = a->sign; - n2 = (a->sign == b->sign ? MP_ZPOS : MP_NEG); - if (c != NULL) { - mp_exch(c, &q); - c->sign = (mp_iszero(c) == MP_YES) ? MP_ZPOS : n2; - } - if (d != NULL) { - mp_exch(d, &ta); - d->sign = (mp_iszero(d) == MP_YES) ? MP_ZPOS : n; - } + } + + while (n-- >= 0) { + if (mp_cmp(&tb, &ta) != MP_GT) { + if (((res = mp_sub(&ta, &tb, &ta)) != MP_OKAY) || + ((res = mp_add(&q, &tq, &q)) != MP_OKAY)) { + goto LBL_ERR; + } + } + if (((res = mp_div_2d(&tb, 1, &tb, NULL)) != MP_OKAY) || + ((res = mp_div_2d(&tq, 1, &tq, NULL)) != MP_OKAY)) { + goto LBL_ERR; + } + } + + /* now q == quotient and ta == remainder */ + n = a->sign; + n2 = (a->sign == b->sign) ? MP_ZPOS : MP_NEG; + if (c != NULL) { + mp_exch(c, &q); + c->sign = (mp_iszero(c) == MP_YES) ? MP_ZPOS : n2; + } + if (d != NULL) { + mp_exch(d, &ta); + d->sign = (mp_iszero(d) == MP_YES) ? MP_ZPOS : n; + } LBL_ERR: mp_clear_multi(&ta, &tb, &tq, &q, NULL); return res; @@ -87,202 +84,214 @@ LBL_ERR: #else -/* integer signed division. +/* integer signed division. * c*b + d == a [e.g. a/b, c=quotient, d=remainder] * HAC pp.598 Algorithm 14.20 * - * Note that the description in HAC is horribly - * incomplete. For example, it doesn't consider - * the case where digits are removed from 'x' in - * the inner loop. It also doesn't consider the + * Note that the description in HAC is horribly + * incomplete. For example, it doesn't consider + * the case where digits are removed from 'x' in + * the inner loop. It also doesn't consider the * case that y has fewer than three digits, etc.. * - * The overall algorithm is as described as + * The overall algorithm is as described as * 14.20 from HAC but fixed to treat these cases. */ -int mp_div (mp_int * a, mp_int * b, mp_int * c, mp_int * d) +int mp_div(const mp_int *a, const mp_int *b, mp_int *c, mp_int *d) { - mp_int q, x, y, t1, t2; - int res, n, t, i, norm, neg; - - /* is divisor zero ? */ - if (mp_iszero (b) == 1) { - return MP_VAL; - } - - /* if a < b then q=0, r = a */ - if (mp_cmp_mag (a, b) == MP_LT) { - if (d != NULL) { - res = mp_copy (a, d); - } else { - res = MP_OKAY; - } - if (c != NULL) { - mp_zero (c); - } - return res; - } - - if ((res = mp_init_size (&q, a->used + 2)) != MP_OKAY) { - return res; - } - q.used = a->used + 2; - - if ((res = mp_init (&t1)) != MP_OKAY) { - goto LBL_Q; - } - - if ((res = mp_init (&t2)) != MP_OKAY) { - goto LBL_T1; - } - - if ((res = mp_init_copy (&x, a)) != MP_OKAY) { - goto LBL_T2; - } - - if ((res = mp_init_copy (&y, b)) != MP_OKAY) { - goto LBL_X; - } - - /* fix the sign */ - neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG; - x.sign = y.sign = MP_ZPOS; - - /* normalize both x and y, ensure that y >= b/2, [b == 2**DIGIT_BIT] */ - norm = mp_count_bits(&y) % DIGIT_BIT; - if (norm < (int)(DIGIT_BIT-1)) { - norm = (DIGIT_BIT-1) - norm; - if ((res = mp_mul_2d (&x, norm, &x)) != MP_OKAY) { - goto LBL_Y; - } - if ((res = mp_mul_2d (&y, norm, &y)) != MP_OKAY) { - goto LBL_Y; - } - } else { - norm = 0; - } - - /* note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4 */ - n = x.used - 1; - t = y.used - 1; - - /* while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} } */ - if ((res = mp_lshd (&y, n - t)) != MP_OKAY) { /* y = y*b**{n-t} */ - goto LBL_Y; - } - - while (mp_cmp (&x, &y) != MP_LT) { - ++(q.dp[n - t]); - if ((res = mp_sub (&x, &y, &x)) != MP_OKAY) { - goto LBL_Y; - } - } - - /* reset y by shifting it back down */ - mp_rshd (&y, n - t); - - /* step 3. for i from n down to (t + 1) */ - for (i = n; i >= (t + 1); i--) { - if (i > x.used) { - continue; - } - - /* step 3.1 if xi == yt then set q{i-t-1} to b-1, - * otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */ - if (x.dp[i] == y.dp[t]) { - q.dp[i - t - 1] = ((((mp_digit)1) << DIGIT_BIT) - 1); - } else { - mp_word tmp; - tmp = ((mp_word) x.dp[i]) << ((mp_word) DIGIT_BIT); - tmp |= ((mp_word) x.dp[i - 1]); - tmp /= ((mp_word) y.dp[t]); - if (tmp > (mp_word) MP_MASK) - tmp = MP_MASK; - q.dp[i - t - 1] = (mp_digit) (tmp & (mp_word) (MP_MASK)); - } - - /* while (q{i-t-1} * (yt * b + y{t-1})) > - xi * b**2 + xi-1 * b + xi-2 - - do q{i-t-1} -= 1; - */ - q.dp[i - t - 1] = (q.dp[i - t - 1] + 1) & MP_MASK; - do { - q.dp[i - t - 1] = (q.dp[i - t - 1] - 1) & MP_MASK; - - /* find left hand */ - mp_zero (&t1); - t1.dp[0] = (t - 1 < 0) ? 0 : y.dp[t - 1]; - t1.dp[1] = y.dp[t]; - t1.used = 2; - if ((res = mp_mul_d (&t1, q.dp[i - t - 1], &t1)) != MP_OKAY) { - goto LBL_Y; + mp_int q, x, y, t1, t2; + int res, n, t, i, norm, neg; + + /* is divisor zero ? */ + if (mp_iszero(b) == MP_YES) { + return MP_VAL; + } + + /* if a < b then q=0, r = a */ + if (mp_cmp_mag(a, b) == MP_LT) { + if (d != NULL) { + res = mp_copy(a, d); + } else { + res = MP_OKAY; + } + if (c != NULL) { + mp_zero(c); + } + return res; + } + + if ((res = mp_init_size(&q, a->used + 2)) != MP_OKAY) { + return res; + } + q.used = a->used + 2; + + if ((res = mp_init(&t1)) != MP_OKAY) { + goto LBL_Q; + } + + if ((res = mp_init(&t2)) != MP_OKAY) { + goto LBL_T1; + } + + if ((res = mp_init_copy(&x, a)) != MP_OKAY) { + goto LBL_T2; + } + + if ((res = mp_init_copy(&y, b)) != MP_OKAY) { + goto LBL_X; + } + + /* fix the sign */ + neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG; + x.sign = y.sign = MP_ZPOS; + + /* normalize both x and y, ensure that y >= b/2, [b == 2**DIGIT_BIT] */ + norm = mp_count_bits(&y) % DIGIT_BIT; + if (norm < (DIGIT_BIT - 1)) { + norm = (DIGIT_BIT - 1) - norm; + if ((res = mp_mul_2d(&x, norm, &x)) != MP_OKAY) { + goto LBL_Y; } + if ((res = mp_mul_2d(&y, norm, &y)) != MP_OKAY) { + goto LBL_Y; + } + } else { + norm = 0; + } - /* find right hand */ - t2.dp[0] = (i - 2 < 0) ? 0 : x.dp[i - 2]; - t2.dp[1] = (i - 1 < 0) ? 0 : x.dp[i - 1]; - t2.dp[2] = x.dp[i]; - t2.used = 3; - } while (mp_cmp_mag(&t1, &t2) == MP_GT); + /* note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4 */ + n = x.used - 1; + t = y.used - 1; - /* step 3.3 x = x - q{i-t-1} * y * b**{i-t-1} */ - if ((res = mp_mul_d (&y, q.dp[i - t - 1], &t1)) != MP_OKAY) { + /* while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} } */ + if ((res = mp_lshd(&y, n - t)) != MP_OKAY) { /* y = y*b**{n-t} */ goto LBL_Y; - } + } - if ((res = mp_lshd (&t1, i - t - 1)) != MP_OKAY) { - goto LBL_Y; - } + while (mp_cmp(&x, &y) != MP_LT) { + ++(q.dp[n - t]); + if ((res = mp_sub(&x, &y, &x)) != MP_OKAY) { + goto LBL_Y; + } + } - if ((res = mp_sub (&x, &t1, &x)) != MP_OKAY) { - goto LBL_Y; - } + /* reset y by shifting it back down */ + mp_rshd(&y, n - t); - /* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */ - if (x.sign == MP_NEG) { - if ((res = mp_copy (&y, &t1)) != MP_OKAY) { - goto LBL_Y; + /* step 3. for i from n down to (t + 1) */ + for (i = n; i >= (t + 1); i--) { + if (i > x.used) { + continue; } - if ((res = mp_lshd (&t1, i - t - 1)) != MP_OKAY) { - goto LBL_Y; + + /* step 3.1 if xi == yt then set q{i-t-1} to b-1, + * otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */ + if (x.dp[i] == y.dp[t]) { + q.dp[(i - t) - 1] = ((mp_digit)1 << (mp_digit)DIGIT_BIT) - (mp_digit)1; + } else { + mp_word tmp; + tmp = (mp_word)x.dp[i] << (mp_word)DIGIT_BIT; + tmp |= (mp_word)x.dp[i - 1]; + tmp /= (mp_word)y.dp[t]; + if (tmp > (mp_word)MP_MASK) { + tmp = MP_MASK; + } + q.dp[(i - t) - 1] = (mp_digit)(tmp & (mp_word)MP_MASK); + } + + /* while (q{i-t-1} * (yt * b + y{t-1})) > + xi * b**2 + xi-1 * b + xi-2 + + do q{i-t-1} -= 1; + */ + q.dp[(i - t) - 1] = (q.dp[(i - t) - 1] + 1uL) & (mp_digit)MP_MASK; + do { + q.dp[(i - t) - 1] = (q.dp[(i - t) - 1] - 1uL) & (mp_digit)MP_MASK; + + /* find left hand */ + mp_zero(&t1); + t1.dp[0] = ((t - 1) < 0) ? 0u : y.dp[t - 1]; + t1.dp[1] = y.dp[t]; + t1.used = 2; + if ((res = mp_mul_d(&t1, q.dp[(i - t) - 1], &t1)) != MP_OKAY) { + goto LBL_Y; + } + + /* find right hand */ + t2.dp[0] = ((i - 2) < 0) ? 0u : x.dp[i - 2]; + t2.dp[1] = ((i - 1) < 0) ? 0u : x.dp[i - 1]; + t2.dp[2] = x.dp[i]; + t2.used = 3; + } while (mp_cmp_mag(&t1, &t2) == MP_GT); + + /* step 3.3 x = x - q{i-t-1} * y * b**{i-t-1} */ + if ((res = mp_mul_d(&y, q.dp[(i - t) - 1], &t1)) != MP_OKAY) { + goto LBL_Y; + } + + if ((res = mp_lshd(&t1, (i - t) - 1)) != MP_OKAY) { + goto LBL_Y; + } + + if ((res = mp_sub(&x, &t1, &x)) != MP_OKAY) { + goto LBL_Y; } - if ((res = mp_add (&x, &t1, &x)) != MP_OKAY) { - goto LBL_Y; + + /* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */ + if (x.sign == MP_NEG) { + if ((res = mp_copy(&y, &t1)) != MP_OKAY) { + goto LBL_Y; + } + if ((res = mp_lshd(&t1, (i - t) - 1)) != MP_OKAY) { + goto LBL_Y; + } + if ((res = mp_add(&x, &t1, &x)) != MP_OKAY) { + goto LBL_Y; + } + + q.dp[(i - t) - 1] = (q.dp[(i - t) - 1] - 1uL) & MP_MASK; } + } - q.dp[i - t - 1] = (q.dp[i - t - 1] - 1UL) & MP_MASK; - } - } - - /* now q is the quotient and x is the remainder - * [which we have to normalize] - */ - - /* get sign before writing to c */ - x.sign = x.used == 0 ? MP_ZPOS : a->sign; - - if (c != NULL) { - mp_clamp (&q); - mp_exch (&q, c); - c->sign = neg; - } - - if (d != NULL) { - mp_div_2d (&x, norm, &x, NULL); - mp_exch (&x, d); - } - - res = MP_OKAY; - -LBL_Y:mp_clear (&y); -LBL_X:mp_clear (&x); -LBL_T2:mp_clear (&t2); -LBL_T1:mp_clear (&t1); -LBL_Q:mp_clear (&q); - return res; + /* now q is the quotient and x is the remainder + * [which we have to normalize] + */ + + /* get sign before writing to c */ + x.sign = (x.used == 0) ? MP_ZPOS : a->sign; + + if (c != NULL) { + mp_clamp(&q); + mp_exch(&q, c); + c->sign = neg; + } + + if (d != NULL) { + if ((res = mp_div_2d(&x, norm, &x, NULL)) != MP_OKAY) { + goto LBL_Y; + } + mp_exch(&x, d); + } + + res = MP_OKAY; + +LBL_Y: + mp_clear(&y); +LBL_X: + mp_clear(&x); +LBL_T2: + mp_clear(&t2); +LBL_T1: + mp_clear(&t1); +LBL_Q: + mp_clear(&q); + return res; } #endif #endif + +/* ref: $Format:%D$ */ +/* git commit: $Format:%H$ */ +/* commit time: $Format:%ai$ */ |