Mercurial > dropbear
annotate bn_mp_div.c @ 200:c5c969ed76f3 libtommath
propagate from branch 'au.asn.ucc.matt.ltm-orig' (head 7fa10cba9535de3461cedb14b877c24858826204)
to branch 'au.asn.ucc.matt.dropbear.ltm' (head fc26f60de0370ab0a281fa41a2d13fb17c9d90a8)
| author | Matt Johnston <matt@ucc.asn.au> |
|---|---|
| date | Wed, 11 May 2005 16:15:27 +0000 |
| parents | d8254fc979e9 |
| children |
| rev | line source |
|---|---|
| 142 | 1 #include <tommath.h> |
| 2 #ifdef BN_MP_DIV_C | |
| 2 | 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, [email protected], http://math.libtomcrypt.org | |
| 16 */ | |
| 142 | 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); | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
52 if (((res = mp_abs(a, &ta)) != MP_OKAY) || |
|
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
53 ((res = mp_abs(b, &tb)) != MP_OKAY) || |
| 142 | 54 ((res = mp_mul_2d(&tb, n, &tb)) != MP_OKAY) || |
| 55 ((res = mp_mul_2d(&tq, n, &tq)) != MP_OKAY)) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
56 goto LBL_ERR; |
| 142 | 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)) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
63 goto LBL_ERR; |
| 142 | 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)) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
68 goto LBL_ERR; |
| 142 | 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); | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
77 c->sign = (mp_iszero(c) == MP_YES) ? MP_ZPOS : n2; |
| 142 | 78 } |
| 79 if (d != NULL) { | |
| 80 mp_exch(d, &ta); | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
81 d->sign = (mp_iszero(d) == MP_YES) ? MP_ZPOS : n; |
| 142 | 82 } |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
83 LBL_ERR: |
| 142 | 84 mp_clear_multi(&ta, &tb, &tq, &q, NULL); |
| 85 return res; | |
| 86 } | |
| 87 | |
| 88 #else | |
| 2 | 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) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
132 goto LBL_Q; |
| 2 | 133 } |
| 134 | |
| 135 if ((res = mp_init (&t2)) != MP_OKAY) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
136 goto LBL_T1; |
| 2 | 137 } |
| 138 | |
| 139 if ((res = mp_init_copy (&x, a)) != MP_OKAY) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
140 goto LBL_T2; |
| 2 | 141 } |
| 142 | |
| 143 if ((res = mp_init_copy (&y, b)) != MP_OKAY) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
144 goto LBL_X; |
| 2 | 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) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
156 goto LBL_Y; |
| 2 | 157 } |
| 158 if ((res = mp_mul_2d (&y, norm, &y)) != MP_OKAY) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
159 goto LBL_Y; |
| 2 | 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} */ | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
171 goto LBL_Y; |
| 2 | 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) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
177 goto LBL_Y; |
| 2 | 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) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
219 goto LBL_Y; |
| 2 | 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) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
231 goto LBL_Y; |
| 2 | 232 } |
| 233 | |
| 234 if ((res = mp_lshd (&t1, i - t - 1)) != MP_OKAY) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
235 goto LBL_Y; |
| 2 | 236 } |
| 237 | |
| 238 if ((res = mp_sub (&x, &t1, &x)) != MP_OKAY) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
239 goto LBL_Y; |
| 2 | 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) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
245 goto LBL_Y; |
| 2 | 246 } |
| 247 if ((res = mp_lshd (&t1, i - t - 1)) != MP_OKAY) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
248 goto LBL_Y; |
| 2 | 249 } |
| 250 if ((res = mp_add (&x, &t1, &x)) != MP_OKAY) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
251 goto LBL_Y; |
| 2 | 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 */ | |
| 142 | 263 x.sign = x.used == 0 ? MP_ZPOS : a->sign; |
| 2 | 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 mp_div_2d (&x, norm, &x, NULL); | |
| 273 mp_exch (&x, d); | |
| 274 } | |
| 275 | |
| 276 res = MP_OKAY; | |
| 277 | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
278 LBL_Y:mp_clear (&y); |
|
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
279 LBL_X:mp_clear (&x); |
|
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
280 LBL_T2:mp_clear (&t2); |
|
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
281 LBL_T1:mp_clear (&t1); |
|
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
282 LBL_Q:mp_clear (&q); |
| 2 | 283 return res; |
| 284 } | |
| 142 | 285 |
| 286 #endif | |
| 287 | |
| 288 #endif |