1/* $NetBSD: bn_mp_div.c,v 1.2 2017/01/28 21:31:47 christos Exp $ */ 2 3#include <tommath.h> 4#ifdef BN_MP_DIV_C 5/* LibTomMath, multiple-precision integer library -- Tom St Denis 6 * 7 * LibTomMath is a library that provides multiple-precision 8 * integer arithmetic as well as number theoretic functionality. 9 * 10 * The library was designed directly after the MPI library by 11 * Michael Fromberger but has been written from scratch with 12 * additional optimizations in place. 13 * 14 * The library is free for all purposes without any express 15 * guarantee it works. 16 * 17 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org 18 */ 19 20#ifdef BN_MP_DIV_SMALL 21 22/* slower bit-bang division... also smaller */ 23int mp_div(mp_int * a, mp_int * b, mp_int * c, mp_int * d) 24{ 25 mp_int ta, tb, tq, q; 26 int res, n, n2; 27 28 /* is divisor zero ? */ 29 if (mp_iszero (b) == 1) { 30 return MP_VAL; 31 } 32 33 /* if a < b then q=0, r = a */ 34 if (mp_cmp_mag (a, b) == MP_LT) { 35 if (d != NULL) { 36 res = mp_copy (a, d); 37 } else { 38 res = MP_OKAY; 39 } 40 if (c != NULL) { 41 mp_zero (c); 42 } 43 return res; 44 } 45 46 /* init our temps */ 47 if ((res = mp_init_multi(&ta, &tb, &tq, &q, NULL) != MP_OKAY)) { 48 return res; 49 } 50 51 52 mp_set(&tq, 1); 53 n = mp_count_bits(a) - mp_count_bits(b); 54 if (((res = mp_abs(a, &ta)) != MP_OKAY) || 55 ((res = mp_abs(b, &tb)) != MP_OKAY) || 56 ((res = mp_mul_2d(&tb, n, &tb)) != MP_OKAY) || 57 ((res = mp_mul_2d(&tq, n, &tq)) != MP_OKAY)) { 58 goto LBL_ERR; 59 } 60 61 while (n-- >= 0) { 62 if (mp_cmp(&tb, &ta) != MP_GT) { 63 if (((res = mp_sub(&ta, &tb, &ta)) != MP_OKAY) || 64 ((res = mp_add(&q, &tq, &q)) != MP_OKAY)) { 65 goto LBL_ERR; 66 } 67 } 68 if (((res = mp_div_2d(&tb, 1, &tb, NULL)) != MP_OKAY) || 69 ((res = mp_div_2d(&tq, 1, &tq, NULL)) != MP_OKAY)) { 70 goto LBL_ERR; 71 } 72 } 73 74 /* now q == quotient and ta == remainder */ 75 n = a->sign; 76 n2 = (a->sign == b->sign ? MP_ZPOS : MP_NEG); 77 if (c != NULL) { 78 mp_exch(c, &q); 79 c->sign = (mp_iszero(c) == MP_YES) ? MP_ZPOS : n2; 80 } 81 if (d != NULL) { 82 mp_exch(d, &ta); 83 d->sign = (mp_iszero(d) == MP_YES) ? MP_ZPOS : n; 84 } 85LBL_ERR: 86 mp_clear_multi(&ta, &tb, &tq, &q, NULL); 87 return res; 88} 89 90#else 91 92/* integer signed division. 93 * c*b + d == a [e.g. a/b, c=quotient, d=remainder] 94 * HAC pp.598 Algorithm 14.20 95 * 96 * Note that the description in HAC is horribly 97 * incomplete. For example, it doesn't consider 98 * the case where digits are removed from 'x' in 99 * the inner loop. It also doesn't consider the 100 * case that y has fewer than three digits, etc.. 101 * 102 * The overall algorithm is as described as 103 * 14.20 from HAC but fixed to treat these cases. 104*/ 105int mp_div (mp_int * a, mp_int * b, mp_int * c, mp_int * d) 106{ 107 mp_int q, x, y, t1, t2; 108 int res, n, t, i, norm, neg; 109 110 /* is divisor zero ? */ 111 if (mp_iszero (b) == 1) { 112 return MP_VAL; 113 } 114 115 /* if a < b then q=0, r = a */ 116 if (mp_cmp_mag (a, b) == MP_LT) { 117 if (d != NULL) { 118 res = mp_copy (a, d); 119 } else { 120 res = MP_OKAY; 121 } 122 if (c != NULL) { 123 mp_zero (c); 124 } 125 return res; 126 } 127 128 if ((res = mp_init_size (&q, a->used + 2)) != MP_OKAY) { 129 return res; 130 } 131 q.used = a->used + 2; 132 133 if ((res = mp_init (&t1)) != MP_OKAY) { 134 goto LBL_Q; 135 } 136 137 if ((res = mp_init (&t2)) != MP_OKAY) { 138 goto LBL_T1; 139 } 140 141 if ((res = mp_init_copy (&x, a)) != MP_OKAY) { 142 goto LBL_T2; 143 } 144 145 if ((res = mp_init_copy (&y, b)) != MP_OKAY) { 146 goto LBL_X; 147 } 148 149 /* fix the sign */ 150 neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG; 151 x.sign = y.sign = MP_ZPOS; 152 153 /* normalize both x and y, ensure that y >= b/2, [b == 2**DIGIT_BIT] */ 154 norm = mp_count_bits(&y) % DIGIT_BIT; 155 if (norm < (int)(DIGIT_BIT-1)) { 156 norm = (DIGIT_BIT-1) - norm; 157 if ((res = mp_mul_2d (&x, norm, &x)) != MP_OKAY) { 158 goto LBL_Y; 159 } 160 if ((res = mp_mul_2d (&y, norm, &y)) != MP_OKAY) { 161 goto LBL_Y; 162 } 163 } else { 164 norm = 0; 165 } 166 167 /* note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4 */ 168 n = x.used - 1; 169 t = y.used - 1; 170 171 /* while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} } */ 172 if ((res = mp_lshd (&y, n - t)) != MP_OKAY) { /* y = y*b**{n-t} */ 173 goto LBL_Y; 174 } 175 176 while (mp_cmp (&x, &y) != MP_LT) { 177 ++(q.dp[n - t]); 178 if ((res = mp_sub (&x, &y, &x)) != MP_OKAY) { 179 goto LBL_Y; 180 } 181 } 182 183 /* reset y by shifting it back down */ 184 mp_rshd (&y, n - t); 185 186 /* step 3. for i from n down to (t + 1) */ 187 for (i = n; i >= (t + 1); i--) { 188 if (i > x.used) { 189 continue; 190 } 191 192 /* step 3.1 if xi == yt then set q{i-t-1} to b-1, 193 * otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */ 194 if (x.dp[i] == y.dp[t]) { 195 q.dp[i - t - 1] = ((((mp_digit)1) << DIGIT_BIT) - 1); 196 } else { 197 mp_word tmp; 198 tmp = ((mp_word) x.dp[i]) << ((mp_word) DIGIT_BIT); 199 tmp |= ((mp_word) x.dp[i - 1]); 200 tmp /= ((mp_word) y.dp[t]); 201 if (tmp > (mp_word) MP_MASK) 202 tmp = MP_MASK; 203 q.dp[i - t - 1] = (mp_digit) (tmp & (mp_word) (MP_MASK)); 204 } 205 206 /* while (q{i-t-1} * (yt * b + y{t-1})) > 207 xi * b**2 + xi-1 * b + xi-2 208 209 do q{i-t-1} -= 1; 210 */ 211 q.dp[i - t - 1] = (q.dp[i - t - 1] + 1) & MP_MASK; 212 do { 213 q.dp[i - t - 1] = (q.dp[i - t - 1] - 1) & MP_MASK; 214 215 /* find left hand */ 216 mp_zero (&t1); 217 t1.dp[0] = (t - 1 < 0) ? 0 : y.dp[t - 1]; 218 t1.dp[1] = y.dp[t]; 219 t1.used = 2; 220 if ((res = mp_mul_d (&t1, q.dp[i - t - 1], &t1)) != MP_OKAY) { 221 goto LBL_Y; 222 } 223 224 /* find right hand */ 225 t2.dp[0] = (i - 2 < 0) ? 0 : x.dp[i - 2]; 226 t2.dp[1] = (i - 1 < 0) ? 0 : x.dp[i - 1]; 227 t2.dp[2] = x.dp[i]; 228 t2.used = 3; 229 } while (mp_cmp_mag(&t1, &t2) == MP_GT); 230 231 /* step 3.3 x = x - q{i-t-1} * y * b**{i-t-1} */ 232 if ((res = mp_mul_d (&y, q.dp[i - t - 1], &t1)) != MP_OKAY) { 233 goto LBL_Y; 234 } 235 236 if ((res = mp_lshd (&t1, i - t - 1)) != MP_OKAY) { 237 goto LBL_Y; 238 } 239 240 if ((res = mp_sub (&x, &t1, &x)) != MP_OKAY) { 241 goto LBL_Y; 242 } 243 244 /* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */ 245 if (x.sign == MP_NEG) { 246 if ((res = mp_copy (&y, &t1)) != MP_OKAY) { 247 goto LBL_Y; 248 } 249 if ((res = mp_lshd (&t1, i - t - 1)) != MP_OKAY) { 250 goto LBL_Y; 251 } 252 if ((res = mp_add (&x, &t1, &x)) != MP_OKAY) { 253 goto LBL_Y; 254 } 255 256 q.dp[i - t - 1] = (q.dp[i - t - 1] - 1UL) & MP_MASK; 257 } 258 } 259 260 /* now q is the quotient and x is the remainder 261 * [which we have to normalize] 262 */ 263 264 /* get sign before writing to c */ 265 x.sign = x.used == 0 ? MP_ZPOS : a->sign; 266 267 if (c != NULL) { 268 mp_clamp (&q); 269 mp_exch (&q, c); 270 c->sign = neg; 271 } 272 273 if (d != NULL) { 274 mp_div_2d (&x, norm, &x, NULL); 275 mp_exch (&x, d); 276 } 277 278 res = MP_OKAY; 279 280LBL_Y:mp_clear (&y); 281LBL_X:mp_clear (&x); 282LBL_T2:mp_clear (&t2); 283LBL_T1:mp_clear (&t1); 284LBL_Q:mp_clear (&q); 285 return res; 286} 287 288#endif 289 290#endif 291 292/* Source: /cvs/libtom/libtommath/bn_mp_div.c,v */ 293/* Revision: 1.4 */ 294/* Date: 2006/12/28 01:25:13 */ 295