Below is the file 'mpi.c' from this revision. You can also download the file.
/* Start: bn_error.c */ /* LibTomMath, multiple-precision integer library -- Tom St Denis * * LibTomMath is a library that provides multiple-precision * integer arithmetic as well as number theoretic functionality. * * The library was designed directly after the MPI library by * 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@iahu.ca, http://math.libtomcrypt.org */ #include <ltc_tommath.h> static const struct { int code; char *msg; } msgs[] = { { MP_OKAY, "Successful" }, { MP_MEM, "Out of heap" }, { MP_VAL, "Value out of range" } }; /* return a char * string for a given code */ char *mp_error_to_string(int code) { int x; /* scan the lookup table for the given message */ for (x = 0; x < (int)(sizeof(msgs) / sizeof(msgs[0])); x++) { if (msgs[x].code == code) { return msgs[x].msg; } } /* generic reply for invalid code */ return "Invalid error code"; } /* End: bn_error.c */ /* Start: bn_fast_mp_invmod.c */ /* LibTomMath, multiple-precision integer library -- Tom St Denis * * LibTomMath is a library that provides multiple-precision * integer arithmetic as well as number theoretic functionality. * * The library was designed directly after the MPI library by * 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@iahu.ca, http://math.libtomcrypt.org */ #include <ltc_tommath.h> /* computes the modular inverse via binary extended euclidean algorithm, * that is c = 1/a mod b * * Based on mp_invmod except this is optimized for the case where b is * odd as per HAC Note 14.64 on pp. 610 */ int fast_mp_invmod (mp_int * a, mp_int * b, mp_int * c) { mp_int x, y, u, v, B, D; int res, neg; /* 2. [modified] b must be odd */ if (mp_iseven (b) == 1) { return MP_VAL; } /* init all our temps */ if ((res = mp_init_multi(&x, &y, &u, &v, &B, &D, NULL)) != MP_OKAY) { return res; } /* x == modulus, y == value to invert */ if ((res = mp_copy (b, &x)) != MP_OKAY) { goto __ERR; } /* we need y = |a| */ if ((res = mp_abs (a, &y)) != MP_OKAY) { goto __ERR; } /* 3. u=x, v=y, A=1, B=0, C=0,D=1 */ if ((res = mp_copy (&x, &u)) != MP_OKAY) { goto __ERR; } if ((res = mp_copy (&y, &v)) != MP_OKAY) { goto __ERR; } mp_set (&D, 1); top: /* 4. while u is even do */ while (mp_iseven (&u) == 1) { /* 4.1 u = u/2 */ if ((res = mp_div_2 (&u, &u)) != MP_OKAY) { goto __ERR; } /* 4.2 if B is odd then */ if (mp_isodd (&B) == 1) { if ((res = mp_sub (&B, &x, &B)) != MP_OKAY) { goto __ERR; } } /* B = B/2 */ if ((res = mp_div_2 (&B, &B)) != MP_OKAY) { goto __ERR; } } /* 5. while v is even do */ while (mp_iseven (&v) == 1) { /* 5.1 v = v/2 */ if ((res = mp_div_2 (&v, &v)) != MP_OKAY) { goto __ERR; } /* 5.2 if D is odd then */ if (mp_isodd (&D) == 1) { /* D = (D-x)/2 */ if ((res = mp_sub (&D, &x, &D)) != MP_OKAY) { goto __ERR; } } /* D = D/2 */ if ((res = mp_div_2 (&D, &D)) != MP_OKAY) { goto __ERR; } } /* 6. if u >= v then */ if (mp_cmp (&u, &v) != MP_LT) { /* u = u - v, B = B - D */ if ((res = mp_sub (&u, &v, &u)) != MP_OKAY) { goto __ERR; } if ((res = mp_sub (&B, &D, &B)) != MP_OKAY) { goto __ERR; } } else { /* v - v - u, D = D - B */ if ((res = mp_sub (&v, &u, &v)) != MP_OKAY) { goto __ERR; } if ((res = mp_sub (&D, &B, &D)) != MP_OKAY) { goto __ERR; } } /* if not zero goto step 4 */ if (mp_iszero (&u) == 0) { goto top; } /* now a = C, b = D, gcd == g*v */ /* if v != 1 then there is no inverse */ if (mp_cmp_d (&v, 1) != MP_EQ) { res = MP_VAL; goto __ERR; } /* b is now the inverse */ neg = a->sign; while (D.sign == MP_NEG) { if ((res = mp_add (&D, b, &D)) != MP_OKAY) { goto __ERR; } } mp_exch (&D, c); c->sign = neg; res = MP_OKAY; __ERR:mp_clear_multi (&x, &y, &u, &v, &B, &D, NULL); return res; } /* End: bn_fast_mp_invmod.c */ /* Start: bn_fast_mp_montgomery_reduce.c */ /* LibTomMath, multiple-precision integer library -- Tom St Denis * * LibTomMath is a library that provides multiple-precision * integer arithmetic as well as number theoretic functionality. * * The library was designed directly after the MPI library by * 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@iahu.ca, http://math.libtomcrypt.org */ #include <ltc_tommath.h> /* computes xR**-1 == x (mod N) via Montgomery Reduction * * This is an optimized implementation of mp_montgomery_reduce * which uses the comba method to quickly calculate the columns of the * reduction. * * Based on Algorithm 14.32 on pp.601 of HAC. */ int fast_mp_montgomery_reduce (mp_int * x, mp_int * n, mp_digit rho) { int ix, res, olduse; mp_word W[MP_WARRAY]; /* get old used count */ olduse = x->used; /* grow a as required */ if (x->alloc < n->used + 1) { if ((res = mp_grow (x, n->used + 1)) != MP_OKAY) { return res; } } /* first we have to get the digits of the input into * an array of double precision words W[...] */ { register mp_word *_W; register mp_digit *tmpx; /* alias for the W[] array */ _W = W; /* alias for the digits of x*/ tmpx = x->dp; /* copy the digits of a into W[0..a->used-1] */ for (ix = 0; ix < x->used; ix++) { *_W++ = *tmpx++; } /* zero the high words of W[a->used..m->used*2] */ for (; ix < n->used * 2 + 1; ix++) { *_W++ = 0; } } /* now we proceed to zero successive digits * from the least significant upwards */ for (ix = 0; ix < n->used; ix++) { /* mu = ai * m' mod b * * We avoid a double precision multiplication (which isn't required) * by casting the value down to a mp_digit. Note this requires * that W[ix-1] have the carry cleared (see after the inner loop) */ register mp_digit mu; mu = (mp_digit) (((W[ix] & MP_MASK) * rho) & MP_MASK); /* a = a + mu * m * b**i * * This is computed in place and on the fly. The multiplication * by b**i is handled by offseting which columns the results * are added to. * * Note the comba method normally doesn't handle carries in the * inner loop In this case we fix the carry from the previous * column since the Montgomery reduction requires digits of the * result (so far) [see above] to work. This is * handled by fixing up one carry after the inner loop. The * carry fixups are done in order so after these loops the * first m->used words of W[] have the carries fixed */ { register int iy; register mp_digit *tmpn; register mp_word *_W; /* alias for the digits of the modulus */ tmpn = n->dp; /* Alias for the columns set by an offset of ix */ _W = W + ix; /* inner loop */ for (iy = 0; iy < n->used; iy++) { *_W++ += ((mp_word)mu) * ((mp_word)*tmpn++); } } /* now fix carry for next digit, W[ix+1] */ W[ix + 1] += W[ix] >> ((mp_word) DIGIT_BIT); } /* now we have to propagate the carries and * shift the words downward [all those least * significant digits we zeroed]. */ { register mp_digit *tmpx; register mp_word *_W, *_W1; /* nox fix rest of carries */ /* alias for current word */ _W1 = W + ix; /* alias for next word, where the carry goes */ _W = W + ++ix; for (; ix <= n->used * 2 + 1; ix++) { *_W++ += *_W1++ >> ((mp_word) DIGIT_BIT); } /* copy out, A = A/b**n * * The result is A/b**n but instead of converting from an * array of mp_word to mp_digit than calling mp_rshd * we just copy them in the right order */ /* alias for destination word */ tmpx = x->dp; /* alias for shifted double precision result */ _W = W + n->used; for (ix = 0; ix < n->used + 1; ix++) { *tmpx++ = (mp_digit)(*_W++ & ((mp_word) MP_MASK)); } /* zero oldused digits, if the input a was larger than * m->used+1 we'll have to clear the digits */ for (; ix < olduse; ix++) { *tmpx++ = 0; } } /* set the max used and clamp */ x->used = n->used + 1; mp_clamp (x); /* if A >= m then A = A - m */ if (mp_cmp_mag (x, n) != MP_LT) { return s_mp_sub (x, n, x); } return MP_OKAY; } /* End: bn_fast_mp_montgomery_reduce.c */ /* Start: bn_fast_s_mp_mul_digs.c */ /* LibTomMath, multiple-precision integer library -- Tom St Denis * * LibTomMath is a library that provides multiple-precision * integer arithmetic as well as number theoretic functionality. * * The library was designed directly after the MPI library by * 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@iahu.ca, http://math.libtomcrypt.org */ #include <ltc_tommath.h> /* Fast (comba) multiplier * * This is the fast column-array [comba] multiplier. It is * designed to compute the columns of the product first * then handle the carries afterwards. This has the effect * of making the nested loops that compute the columns very * simple and schedulable on super-scalar processors. * * This has been modified to produce a variable number of * digits of output so if say only a half-product is required * you don't have to compute the upper half (a feature * required for fast Barrett reduction). * * Based on Algorithm 14.12 on pp.595 of HAC. * */ int fast_s_mp_mul_digs (mp_int * a, mp_int * b, mp_int * c, int digs) { int olduse, res, pa, ix; mp_word W[MP_WARRAY]; /* grow the destination as required */ if (c->alloc < digs) { if ((res = mp_grow (c, digs)) != MP_OKAY) { return res; } } /* clear temp buf (the columns) */ memset (W, 0, sizeof (mp_word) * digs); /* calculate the columns */ pa = a->used; for (ix = 0; ix < pa; ix++) { /* this multiplier has been modified to allow you to * control how many digits of output are produced. * So at most we want to make upto "digs" digits of output. * * this adds products to distinct columns (at ix+iy) of W * note that each step through the loop is not dependent on * the previous which means the compiler can easily unroll * the loop without scheduling problems */ { register mp_digit tmpx, *tmpy; register mp_word *_W; register int iy, pb; /* alias for the the word on the left e.g. A[ix] * A[iy] */ tmpx = a->dp[ix]; /* alias for the right side */ tmpy = b->dp; /* alias for the columns, each step through the loop adds a new term to each column */ _W = W + ix; /* the number of digits is limited by their placement. E.g. we avoid multiplying digits that will end up above the # of digits of precision requested */ pb = MIN (b->used, digs - ix); for (iy = 0; iy < pb; iy++) { *_W++ += ((mp_word)tmpx) * ((mp_word)*tmpy++); } } } /* setup dest */ olduse = c->used; c->used = digs; { register mp_digit *tmpc; /* At this point W[] contains the sums of each column. To get the * correct result we must take the extra bits from each column and * carry them down * * Note that while this adds extra code to the multiplier it * saves time since the carry propagation is removed from the * above nested loop.This has the effect of reducing the work * from N*(N+N*c)==N**2 + c*N**2 to N**2 + N*c where c is the * cost of the shifting. On very small numbers this is slower * but on most cryptographic size numbers it is faster. * * In this particular implementation we feed the carries from * behind which means when the loop terminates we still have one * last digit to copy */ tmpc = c->dp; for (ix = 1; ix < digs; ix++) { /* forward the carry from the previous temp */ W[ix] += (W[ix - 1] >> ((mp_word) DIGIT_BIT)); /* now extract the previous digit [below the carry] */ *tmpc++ = (mp_digit) (W[ix - 1] & ((mp_word) MP_MASK)); } /* fetch the last digit */ *tmpc++ = (mp_digit) (W[digs - 1] & ((mp_word) MP_MASK)); /* clear unused digits [that existed in the old copy of c] */ for (; ix < olduse; ix++) { *tmpc++ = 0; } } mp_clamp (c); return MP_OKAY; } /* End: bn_fast_s_mp_mul_digs.c */ /* Start: bn_fast_s_mp_mul_high_digs.c */ /* LibTomMath, multiple-precision integer library -- Tom St Denis * * LibTomMath is a library that provides multiple-precision * integer arithmetic as well as number theoretic functionality. * * The library was designed directly after the MPI library by * 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@iahu.ca, http://math.libtomcrypt.org */ #include <ltc_tommath.h> /* this is a modified version of fast_s_mp_mul_digs that only produces * output digits *above* digs. See the comments for fast_s_mp_mul_digs * to see how it works. * * This is used in the Barrett reduction since for one of the multiplications * only the higher digits were needed. This essentially halves the work. * * Based on Algorithm 14.12 on pp.595 of HAC. */ int fast_s_mp_mul_high_digs (mp_int * a, mp_int * b, mp_int * c, int digs) { int oldused, newused, res, pa, pb, ix; mp_word W[MP_WARRAY]; /* calculate size of product and allocate more space if required */ newused = a->used + b->used + 1; if (c->alloc < newused) { if ((res = mp_grow (c, newused)) != MP_OKAY) { return res; } } /* like the other comba method we compute the columns first */ pa = a->used; pb = b->used; memset (W + digs, 0, (pa + pb + 1 - digs) * sizeof (mp_word)); for (ix = 0; ix < pa; ix++) { { register mp_digit tmpx, *tmpy; register int iy; register mp_word *_W; /* work todo, that is we only calculate digits that are at "digs" or above */ iy = digs - ix; /* copy of word on the left of A[ix] * B[iy] */ tmpx = a->dp[ix]; /* alias for right side */ tmpy = b->dp + iy; /* alias for the columns of output. Offset to be equal to or above the * smallest digit place requested */ _W = W + digs; /* skip cases below zero where ix > digs */ if (iy < 0) { iy = abs(iy); tmpy += iy; _W += iy; iy = 0; } /* compute column products for digits above the minimum */ for (; iy < pb; iy++) { *_W++ += ((mp_word) tmpx) * ((mp_word)*tmpy++); } } } /* setup dest */ oldused = c->used; c->used = newused; /* now convert the array W downto what we need * * See comments in bn_fast_s_mp_mul_digs.c */ for (ix = digs + 1; ix < newused; ix++) { W[ix] += (W[ix - 1] >> ((mp_word) DIGIT_BIT)); c->dp[ix - 1] = (mp_digit) (W[ix - 1] & ((mp_word) MP_MASK)); } c->dp[newused - 1] = (mp_digit) (W[newused - 1] & ((mp_word) MP_MASK)); for (; ix < oldused; ix++) { c->dp[ix] = 0; } mp_clamp (c); return MP_OKAY; } /* End: bn_fast_s_mp_mul_high_digs.c */ /* Start: bn_fast_s_mp_sqr.c */ /* LibTomMath, multiple-precision integer library -- Tom St Denis * * LibTomMath is a library that provides multiple-precision * integer arithmetic as well as number theoretic functionality. * * The library was designed directly after the MPI library by * 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@iahu.ca, http://math.libtomcrypt.org */ #include <ltc_tommath.h> /* fast squaring * * This is the comba method where the columns of the product * are computed first then the carries are computed. This * has the effect of making a very simple inner loop that * is executed the most * * W2 represents the outer products and W the inner. * * A further optimizations is made because the inner * products are of the form "A * B * 2". The *2 part does * not need to be computed until the end which is good * because 64-bit shifts are slow! * * Based on Algorithm 14.16 on pp.597 of HAC. * */ int fast_s_mp_sqr (mp_int * a, mp_int * b) { int olduse, newused, res, ix, pa; mp_word W2[MP_WARRAY], W[MP_WARRAY]; /* calculate size of product and allocate as required */ pa = a->used; newused = pa + pa + 1; if (b->alloc < newused) { if ((res = mp_grow (b, newused)) != MP_OKAY) { return res; } } /* zero temp buffer (columns) * Note that there are two buffers. Since squaring requires * a outer and inner product and the inner product requires * computing a product and doubling it (a relatively expensive * op to perform n**2 times if you don't have to) the inner and * outer products are computed in different buffers. This way * the inner product can be doubled using n doublings instead of * n**2 */ memset (W, 0, newused * sizeof (mp_word)); memset (W2, 0, newused * sizeof (mp_word)); /* This computes the inner product. To simplify the inner N**2 loop * the multiplication by two is done afterwards in the N loop. */ for (ix = 0; ix < pa; ix++) { /* compute the outer product * * Note that every outer product is computed * for a particular column only once which means that * there is no need todo a double precision addition * into the W2[] array. */ W2[ix + ix] = ((mp_word)a->dp[ix]) * ((mp_word)a->dp[ix]); { register mp_digit tmpx, *tmpy; register mp_word *_W; register int iy; /* copy of left side */ tmpx = a->dp[ix]; /* alias for right side */ tmpy = a->dp + (ix + 1); /* the column to store the result in */ _W = W + (ix + ix + 1); /* inner products */ for (iy = ix + 1; iy < pa; iy++) { *_W++ += ((mp_word)tmpx) * ((mp_word)*tmpy++); } } } /* setup dest */ olduse = b->used; b->used = newused; /* now compute digits * * We have to double the inner product sums, add in the * outer product sums, propagate carries and convert * to single precision. */ { register mp_digit *tmpb; /* double first value, since the inner products are * half of what they should be */ W[0] += W[0] + W2[0]; tmpb = b->dp; for (ix = 1; ix < newused; ix++) { /* double/add next digit */ W[ix] += W[ix] + W2[ix]; /* propagate carry forwards [from the previous digit] */ W[ix] = W[ix] + (W[ix - 1] >> ((mp_word) DIGIT_BIT)); /* store the current digit now that the carry isn't * needed */ *tmpb++ = (mp_digit) (W[ix - 1] & ((mp_word) MP_MASK)); } /* set the last value. Note even if the carry is zero * this is required since the next step will not zero * it if b originally had a value at b->dp[2*a.used] */ *tmpb++ = (mp_digit) (W[(newused) - 1] & ((mp_word) MP_MASK)); /* clear high digits of b if there were any originally */ for (; ix < olduse; ix++) { *tmpb++ = 0; } } mp_clamp (b); return MP_OKAY; } /* End: bn_fast_s_mp_sqr.c */ /* Start: bn_mp_2expt.c */ /* LibTomMath, multiple-precision integer library -- Tom St Denis * * LibTomMath is a library that provides multiple-precision * integer arithmetic as well as number theoretic functionality. * * The library was designed directly after the MPI library by * 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@iahu.ca, http://math.libtomcrypt.org */ #include <ltc_tommath.h> /* computes a = 2**b * * Simple algorithm which zeroes the int, grows it then just sets one bit * as required. */ int mp_2expt (mp_int * a, int b) { int res; /* zero a as per default */ mp_zero (a); /* grow a to accomodate the single bit */ if ((res = mp_grow (a, b / DIGIT_BIT + 1)) != MP_OKAY) { return res; } /* set the used count of where the bit will go */ a->used = b / DIGIT_BIT + 1; /* put the single bit in its place */ a->dp[b / DIGIT_BIT] = 1 << (b % DIGIT_BIT); return MP_OKAY; } /* End: bn_mp_2expt.c */ /* Start: bn_mp_abs.c */ /* LibTomMath, multiple-precision integer library -- Tom St Denis * * LibTomMath is a library that provides multiple-precision * integer arithmetic as well as number theoretic functionality. * * The library was designed directly after the MPI library by * 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@iahu.ca, http://math.libtomcrypt.org */ #include <ltc_tommath.h> /* b = |a| * * Simple function copies the input and fixes the sign to positive */ int mp_abs (mp_int * a, mp_int * b) { int res; /* copy a to b */ if (a != b) { if ((res = mp_copy (a, b)) != MP_OKAY) { return res; } } /* force the sign of b to positive */ b->sign = MP_ZPOS; return MP_OKAY; } /* End: bn_mp_abs.c */ /* Start: bn_mp_add.c */ /* LibTomMath, multiple-precision integer library -- Tom St Denis * * LibTomMath is a library that provides multiple-precision * integer arithmetic as well as number theoretic functionality. * * The library was designed directly after the MPI library by * 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@iahu.ca, http://math.libtomcrypt.org */ #include <ltc_tommath.h> /* high level addition (handles signs) */ int mp_add (mp_int * a, mp_int * b, mp_int * c) { int sa, sb, res; /* get sign of both inputs */ sa = a->sign; sb = b->sign; /* handle two cases, not four */ if (sa == sb) { /* both positive or both negative */ /* add their magnitudes, copy the sign */ c->sign = sa; res = s_mp_add (a, b, c); } else { /* one positive, the other negative */ /* subtract the one with the greater magnitude from */ /* the one of the lesser magnitude. The result gets */ /* the sign of the one with the greater magnitude. */ if (mp_cmp_mag (a, b) == MP_LT) { c->sign = sb; res = s_mp_sub (b, a, c); } else { c->sign = sa; res = s_mp_sub (a, b, c); } } return res; } /* End: bn_mp_add.c */ /* Start: bn_mp_add_d.c */ /* LibTomMath, multiple-precision integer library -- Tom St Denis * * LibTomMath is a library that provides multiple-precision * integer arithmetic as well as number theoretic functionality. * * The library was designed directly after the MPI library by * 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@iahu.ca, http://math.libtomcrypt.org */ #include <ltc_tommath.h> /* single digit addition */ int mp_add_d (mp_int * a, mp_digit b, mp_int * c) { int res, ix, oldused; mp_digit *tmpa, *tmpc, mu; /* grow c as required */ if (c->alloc < a->used + 1) { if ((res = mp_grow(c, a->used + 1)) != MP_OKAY) { return res; } } /* if a is negative and |a| >= b, call c = |a| - b */ if (a->sign == MP_NEG && (a->used > 1 || a->dp[0] >= b)) { /* temporarily fix sign of a */ a->sign = MP_ZPOS; /* c = |a| - b */ res = mp_sub_d(a, b, c); /* fix sign */ a->sign = c->sign = MP_NEG; return res; } /* old number of used digits in c */ oldused = c->used; /* sign always positive */ c->sign = MP_ZPOS; /* source alias */ tmpa = a->dp; /* destination alias */ tmpc = c->dp; /* if a is positive */ if (a->sign == MP_ZPOS) { /* add digit, after this we're propagating * the carry. */ *tmpc = *tmpa++ + b; mu = *tmpc >> DIGIT_BIT; *tmpc++ &= MP_MASK; /* now handle rest of the digits */ for (ix = 1; ix < a->used; ix++) { *tmpc = *tmpa++ + mu; mu = *tmpc >> DIGIT_BIT; *tmpc++ &= MP_MASK; } /* set final carry */ ix++; *tmpc++ = mu; /* setup size */ c->used = a->used + 1; } else { /* a was negative and |a| < b */ c->used = 1; /* the result is a single digit */ if (a->used == 1) { *tmpc++ = b - a->dp[0]; } else { *tmpc++ = b; } /* setup count so the clearing of oldused * can fall through correctly */ ix = 1; } /* now zero to oldused */ while (ix++ < oldused) { *tmpc++ = 0; } mp_clamp(c); return MP_OKAY; } /* End: bn_mp_add_d.c */ /* Start: bn_mp_addmod.c */ /* LibTomMath, multiple-precision integer library -- Tom St Denis * * LibTomMath is a library that provides multiple-precision * integer arithmetic as well as number theoretic functionality. * * The library was designed directly after the MPI library by * 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@iahu.ca, http://math.libtomcrypt.org */ #include <ltc_tommath.h> /* d = a + b (mod c) */ int mp_addmod (mp_int * a, mp_int * b, mp_int * c, mp_int * d) { int res; mp_int t; if ((res = mp_init (&t)) != MP_OKAY) { return res; } if ((res = mp_add (a, b, &t)) != MP_OKAY) { mp_clear (&t); return res; } res = mp_mod (&t, c, d); mp_clear (&t); return res; } /* End: bn_mp_addmod.c */ /* Start: bn_mp_and.c */ /* LibTomMath, multiple-precision integer library -- Tom St Denis * * LibTomMath is a library that provides multiple-precision * integer arithmetic as well as number theoretic functionality. * * The library was designed directly after the MPI library by * 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@iahu.ca, http://math.libtomcrypt.org */ #include <ltc_tommath.h> /* AND two ints together */ int mp_and (mp_int * a, mp_int * b, mp_int * c) { int res, ix, px; mp_int t, *x; if (a->used > b->used) { if ((res = mp_init_copy (&t, a)) != MP_OKAY) { return res; } px = b->used; x = b; } else { if ((res = mp_init_copy (&t, b)) != MP_OKAY) { return res; } px = a->used; x = a; } for (ix = 0; ix < px; ix++) { t.dp[ix] &= x->dp[ix]; } /* zero digits above the last from the smallest mp_int */ for (; ix < t.used; ix++) { t.dp[ix] = 0; } mp_clamp (&t); mp_exch (c, &t); mp_clear (&t); return MP_OKAY; } /* End: bn_mp_and.c */ /* Start: bn_mp_clamp.c */ /* LibTomMath, multiple-precision integer library -- Tom St Denis * * LibTomMath is a library that provides multiple-precision * integer arithmetic as well as number theoretic functionality. * * The library was designed directly after the MPI library by * 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@iahu.ca, http://math.libtomcrypt.org */ #include <ltc_tommath.h> /* trim unused digits * * This is used to ensure that leading zero digits are * trimed and the leading "used" digit will be non-zero * Typically very fast. Also fixes the sign if there * are no more leading digits */ void mp_clamp (mp_int * a) { /* decrease used while the most significant digit is * zero. */ while (a->used > 0 && a->dp[a->used - 1] == 0) { --(a->used); } /* reset the sign flag if used == 0 */ if (a->used == 0) { a->sign = MP_ZPOS; } } /* End: bn_mp_clamp.c */ /* Start: bn_mp_clear.c */ /* LibTomMath, multiple-precision integer library -- Tom St Denis * * LibTomMath is a library that provides multiple-precision * integer arithmetic as well as number theoretic functionality. * * The library was designed directly after the MPI library by * 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@iahu.ca, http://math.libtomcrypt.org */ #include <ltc_tommath.h> /* clear one (frees) */ void mp_clear (mp_int * a) { /* only do anything if a hasn't been freed previously */ if (a->dp != NULL) { /* first zero the digits */ memset (a->dp, 0, sizeof (mp_digit) * a->used); /* free ram */ XFREE(a->dp); /* reset members to make debugging easier */ a->dp = NULL; a->alloc = a->used = 0; a->sign = MP_ZPOS; } } /* End: bn_mp_clear.c */ /* Start: bn_mp_clear_multi.c */ /* LibTomMath, multiple-precision integer library -- Tom St Denis * * LibTomMath is a library that provides multiple-precision * integer arithmetic as well as number theoretic functionality. * * The library was designed directly after the MPI library by * 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@iahu.ca, http://math.libtomcrypt.org */ #include <ltc_tommath.h> #include <stdarg.h> void mp_clear_multi(mp_int *mp, ...) { mp_int* next_mp = mp; va_list args; va_start(args, mp); while (next_mp != NULL) { mp_clear(next_mp); next_mp = va_arg(args, mp_int*); } va_end(args); } /* End: bn_mp_clear_multi.c */ /* Start: bn_mp_cmp.c */ /* LibTomMath, multiple-precision integer library -- Tom St Denis * * LibTomMath is a library that provides multiple-precision * integer arithmetic as well as number theoretic functionality. * * The library was designed directly after the MPI library by * 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@iahu.ca, http://math.libtomcrypt.org */ #include <ltc_tommath.h> /* compare two ints (signed)*/ int mp_cmp (mp_int * a, mp_int * b) { /* compare based on sign */ if (a->sign != b->sign) { if (a->sign == MP_NEG) { return MP_LT; } else { return MP_GT; } } /* compare digits */ if (a->sign == MP_NEG) { /* if negative compare opposite direction */ return mp_cmp_mag(b, a); } else { return mp_cmp_mag(a, b); } } /* End: bn_mp_cmp.c */ /* Start: bn_mp_cmp_d.c */ /* LibTomMath, multiple-precision integer library -- Tom St Denis * * LibTomMath is a library that provides multiple-precision * integer arithmetic as well as number theoretic functionality. * * The library was designed directly after the MPI library by * 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@iahu.ca, http://math.libtomcrypt.org */ #include <ltc_tommath.h> /* compare a digit */ int mp_cmp_d(mp_int * a, mp_digit b) { /* compare based on sign */ if (a->sign == MP_NEG) { return MP_LT; } /* compare based on magnitude */ if (a->used > 1) { return MP_GT; } /* compare the only digit of a to b */ if (a->dp[0] > b) { return MP_GT; } else if (a->dp[0] < b) { return MP_LT; } else { return MP_EQ; } } /* End: bn_mp_cmp_d.c */ /* Start: bn_mp_cmp_mag.c */ /* LibTomMath, multiple-precision integer library -- Tom St Denis * * LibTomMath is a library that provides multiple-precision * integer arithmetic as well as number theoretic functionality. * * The library was designed directly after the MPI library by * 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@iahu.ca, http://math.libtomcrypt.org */ #include <ltc_tommath.h> /* compare maginitude of two ints (unsigned) */ int mp_cmp_mag (mp_int * a, mp_int * b) { int n; mp_digit *tmpa, *tmpb; /* compare based on # of non-zero digits */ if (a->used > b->used) { return MP_GT; } if (a->used < b->used) { return MP_LT; } /* alias for a */ tmpa = a->dp + (a->used - 1); /* alias for b */ tmpb = b->dp + (a->used - 1); /* compare based on digits */ for (n = 0; n < a->used; ++n, --tmpa, --tmpb) { if (*tmpa > *tmpb) { return MP_GT; } if (*tmpa < *tmpb) { return MP_LT; } } return MP_EQ; } /* End: bn_mp_cmp_mag.c */ /* Start: bn_mp_cnt_lsb.c */ /* LibTomMath, multiple-precision integer library -- Tom St Denis * * LibTomMath is a library that provides multiple-precision * integer arithmetic as well as number theoretic functionality. * * The library was designed directly after the MPI library by * 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@iahu.ca, http://math.libtomcrypt.org */ #include <ltc_tommath.h> static const int lnz[16] = { 4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0 }; /* Counts the number of lsbs which are zero before the first zero bit */ int mp_cnt_lsb(mp_int *a) { int x; mp_digit q, qq; /* easy out */ if (mp_iszero(a) == 1) { return 0; } /* scan lower digits until non-zero */ for (x = 0; x < a->used && a->dp[x] == 0; x++); q = a->dp[x]; x *= DIGIT_BIT; /* now scan this digit until a 1 is found */ if ((q & 1) == 0) { do { qq = q & 15; x += lnz[qq]; q >>= 4; } while (qq == 0); } return x; } /* End: bn_mp_cnt_lsb.c */ /* Start: bn_mp_copy.c */ /* LibTomMath, multiple-precision integer library -- Tom St Denis * * LibTomMath is a library that provides multiple-precision * integer arithmetic as well as number theoretic functionality. * * The library was designed directly after the MPI library by * 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@iahu.ca, http://math.libtomcrypt.org */ #include <ltc_tommath.h> /* copy, b = a */ int mp_copy (mp_int * a, mp_int * b) { int res, n; /* if dst == src do nothing */ if (a == b) { return MP_OKAY; } /* grow dest */ if (b->alloc < a->used) { if ((res = mp_grow (b, a->used)) != MP_OKAY) { return res; } } /* zero b and copy the parameters over */ { register mp_digit *tmpa, *tmpb; /* pointer aliases */ /* source */ tmpa = a->dp; /* destination */ tmpb = b->dp; /* copy all the digits */ for (n = 0; n < a->used; n++) { *tmpb++ = *tmpa++; } /* clear high digits */ for (; n < b->used; n++) { *tmpb++ = 0; } } /* copy used count and sign */ b->used = a->used; b->sign = a->sign; return MP_OKAY; } /* End: bn_mp_copy.c */ /* Start: bn_mp_count_bits.c */ /* LibTomMath, multiple-precision integer library -- Tom St Denis * * LibTomMath is a library that provides multiple-precision * integer arithmetic as well as number theoretic functionality. * * The library was designed directly after the MPI library by * 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@iahu.ca, http://math.libtomcrypt.org */ #include <ltc_tommath.h> /* returns the number of bits in an int */ int mp_count_bits (mp_int * a) { int r; mp_digit q; /* shortcut */ if (a->used == 0) { return 0; } /* get number of digits and add that */ r = (a->used - 1) * DIGIT_BIT; /* take the last digit and count the bits in it */ q = a->dp[a->used - 1]; while (q > ((mp_digit) 0)) { ++r; q >>= ((mp_digit) 1); } return r; } /* End: bn_mp_count_bits.c */ /* Start: bn_mp_div.c */ /* LibTomMath, multiple-precision integer library -- Tom St Denis * * LibTomMath is a library that provides multiple-precision * integer arithmetic as well as number theoretic functionality. * * The library was designed directly after the MPI library by * 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@iahu.ca, http://math.libtomcrypt.org */ #include <ltc_tommath.h> /* 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 * case that y has fewer than three digits, etc.. * * 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) { 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 __Q; } if ((res = mp_init (&t2)) != MP_OKAY) { goto __T1; } if ((res = mp_init_copy (&x, a)) != MP_OKAY) { goto __T2; } if ((res = mp_init_copy (&y, b)) != MP_OKAY) { goto __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 __Y; } if ((res = mp_mul_2d (&y, norm, &y)) != MP_OKAY) { goto __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 __Y; } while (mp_cmp (&x, &y) != MP_LT) { ++(q.dp[n - t]); if ((res = mp_sub (&x, &y, &x)) != MP_OKAY) { goto __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 __Y; } /* 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); /* 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 __Y; } if ((res = mp_lshd (&t1, i - t - 1)) != MP_OKAY) { goto __Y; } if ((res = mp_sub (&x, &t1, &x)) != MP_OKAY) { goto __Y; } /* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */ if (x.si