1 #include <tommath.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, tomstdenis (at) gmail.com, http://math.libtomcrypt.com 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) == 1) { 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) == 1) { 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 q.dp[i - t - 1] = (mp_digit) (tmp & (mp_word) (MP_MASK)); 202 } 203 204 /* while (q{i-t-1} * (yt * b + y{t-1})) > 205 xi * b**2 + xi-1 * b + xi-2 206 207 do q{i-t-1} -= 1; 208 */ 209 q.dp[i - t - 1] = (q.dp[i - t - 1] + 1) & MP_MASK; 210 do { 211 q.dp[i - t - 1] = (q.dp[i - t - 1] - 1) & MP_MASK; 212 213 /* find left hand */ 214 mp_zero (&t1); 215 t1.dp[0] = (t - 1 < 0) ? 0 : y.dp[t - 1]; 216 t1.dp[1] = y.dp[t]; 217 t1.used = 2; 218 if ((res = mp_mul_d (&t1, q.dp[i - t - 1], &t1)) != MP_OKAY) { 219 goto LBL_Y; 220 } 221 222 /* find right hand */ 223 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[2] = x.dp[i]; 226 t2.used = 3; 227 } while (mp_cmp_mag(&t1, &t2) == MP_GT); 228 229 /* 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 goto LBL_Y; 232 } 233 234 if ((res = mp_lshd (&t1, i - t - 1)) != MP_OKAY) { 235 goto LBL_Y; 236 } 237 238 if ((res = mp_sub (&x, &t1, &x)) != MP_OKAY) { 239 goto LBL_Y; 240 } 241 242 /* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */ 243 if (x.sign == MP_NEG) { 244 if ((res = mp_copy (&y, &t1)) != MP_OKAY) { 245 goto LBL_Y; 246 } 247 if ((res = mp_lshd (&t1, i - t - 1)) != MP_OKAY) { 248 goto LBL_Y; 249 } 250 if ((res = mp_add (&x, &t1, &x)) != MP_OKAY) { 251 goto LBL_Y; 252 } 253 254 q.dp[i - t - 1] = (q.dp[i - t - 1] - 1UL) & MP_MASK; 255 } 256 } 257 258 /* now q is the quotient and x is the remainder 259 * [which we have to normalize] 260 */ 261 262 /* get sign before writing to c */ 263 x.sign = x.used == 0 ? MP_ZPOS : a->sign; 264 265 if (c != NULL) { 266 mp_clamp (&q); 267 mp_exch (&q, c); 268 c->sign = neg; 269 } 270 271 if (d != NULL) { 272 if ((res = mp_div_2d (&x, norm, &x, NULL)) != MP_OKAY) { 273 goto LBL_Y; 274 } 275 mp_exch (&x, d); 276 } 277 278 res = MP_OKAY; 279 280 LBL_Y:mp_clear (&y); 281 LBL_X:mp_clear (&x); 282 LBL_T2:mp_clear (&t2); 283 LBL_T1:mp_clear (&t1); 284 LBL_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.3 $ */ 294 /* $Date: 2006/03/31 14:18:44 $ */ 295
