X-Git-Url: https://asedeno.scripts.mit.edu/gitweb/?a=blobdiff_plain;f=sshbn.c;h=6768204bc4cd2ec930799bcc3d250407928d63a2;hb=510f49e405e71ba5c97875e7a019364e1ef5fac9;hp=a35bf31444a23e788e1fb09fa16e7d4d7b2861c9;hpb=89da2ddf564a93414ee9ab2df3f053608094e417;p=PuTTY.git diff --git a/sshbn.c b/sshbn.c index a35bf314..6768204b 100644 --- a/sshbn.c +++ b/sshbn.c @@ -11,90 +11,7 @@ #include "misc.h" -/* - * Usage notes: - * * Do not call the DIVMOD_WORD macro with expressions such as array - * subscripts, as some implementations object to this (see below). - * * Note that none of the division methods below will cope if the - * quotient won't fit into BIGNUM_INT_BITS. Callers should be careful - * to avoid this case. - * If this condition occurs, in the case of the x86 DIV instruction, - * an overflow exception will occur, which (according to a correspondent) - * will manifest on Windows as something like - * 0xC0000095: Integer overflow - * The C variant won't give the right answer, either. - */ - -#if defined __GNUC__ && defined __i386__ -typedef unsigned long BignumInt; -typedef unsigned long long BignumDblInt; -#define BIGNUM_INT_MASK 0xFFFFFFFFUL -#define BIGNUM_TOP_BIT 0x80000000UL -#define BIGNUM_INT_BITS 32 -#define MUL_WORD(w1, w2) ((BignumDblInt)w1 * w2) -#define DIVMOD_WORD(q, r, hi, lo, w) \ - __asm__("div %2" : \ - "=d" (r), "=a" (q) : \ - "r" (w), "d" (hi), "a" (lo)) -#elif defined _MSC_VER && defined _M_IX86 -typedef unsigned __int32 BignumInt; -typedef unsigned __int64 BignumDblInt; -#define BIGNUM_INT_MASK 0xFFFFFFFFUL -#define BIGNUM_TOP_BIT 0x80000000UL -#define BIGNUM_INT_BITS 32 -#define MUL_WORD(w1, w2) ((BignumDblInt)w1 * w2) -/* Note: MASM interprets array subscripts in the macro arguments as - * assembler syntax, which gives the wrong answer. Don't supply them. - * */ -#define DIVMOD_WORD(q, r, hi, lo, w) do { \ - __asm mov edx, hi \ - __asm mov eax, lo \ - __asm div w \ - __asm mov r, edx \ - __asm mov q, eax \ -} while(0) -#elif defined _LP64 -/* 64-bit architectures can do 32x32->64 chunks at a time */ -typedef unsigned int BignumInt; -typedef unsigned long BignumDblInt; -#define BIGNUM_INT_MASK 0xFFFFFFFFU -#define BIGNUM_TOP_BIT 0x80000000U -#define BIGNUM_INT_BITS 32 -#define MUL_WORD(w1, w2) ((BignumDblInt)w1 * w2) -#define DIVMOD_WORD(q, r, hi, lo, w) do { \ - BignumDblInt n = (((BignumDblInt)hi) << BIGNUM_INT_BITS) | lo; \ - q = n / w; \ - r = n % w; \ -} while (0) -#elif defined _LLP64 -/* 64-bit architectures in which unsigned long is 32 bits, not 64 */ -typedef unsigned long BignumInt; -typedef unsigned long long BignumDblInt; -#define BIGNUM_INT_MASK 0xFFFFFFFFUL -#define BIGNUM_TOP_BIT 0x80000000UL -#define BIGNUM_INT_BITS 32 -#define MUL_WORD(w1, w2) ((BignumDblInt)w1 * w2) -#define DIVMOD_WORD(q, r, hi, lo, w) do { \ - BignumDblInt n = (((BignumDblInt)hi) << BIGNUM_INT_BITS) | lo; \ - q = n / w; \ - r = n % w; \ -} while (0) -#else -/* Fallback for all other cases */ -typedef unsigned short BignumInt; -typedef unsigned long BignumDblInt; -#define BIGNUM_INT_MASK 0xFFFFU -#define BIGNUM_TOP_BIT 0x8000U -#define BIGNUM_INT_BITS 16 -#define MUL_WORD(w1, w2) ((BignumDblInt)w1 * w2) -#define DIVMOD_WORD(q, r, hi, lo, w) do { \ - BignumDblInt n = (((BignumDblInt)hi) << BIGNUM_INT_BITS) | lo; \ - q = n / w; \ - r = n % w; \ -} while (0) -#endif - -#define BIGNUM_INT_BYTES (BIGNUM_INT_BITS / 8) +#include "sshbn.h" #define BIGNUM_INTERNAL typedef BignumInt *Bignum; @@ -128,8 +45,6 @@ static Bignum newbn(int length) assert(length >= 0 && length < INT_MAX / BIGNUM_INT_BITS); b = snewn(length + 1, BignumInt); - if (!b) - abort(); /* FIXME */ memset(b, 0, (length + 1) * sizeof(*b)); b[0] = length; return b; @@ -172,20 +87,17 @@ Bignum bn_power_2(int n) /* * Internal addition. Sets c = a - b, where 'a', 'b' and 'c' are all - * big-endian arrays of 'len' BignumInts. Returns a BignumInt carried - * off the top. + * big-endian arrays of 'len' BignumInts. Returns the carry off the + * top. */ -static BignumInt internal_add(const BignumInt *a, const BignumInt *b, - BignumInt *c, int len) +static BignumCarry internal_add(const BignumInt *a, const BignumInt *b, + BignumInt *c, int len) { int i; - BignumDblInt carry = 0; + BignumCarry carry = 0; - for (i = len-1; i >= 0; i--) { - carry += (BignumDblInt)a[i] + b[i]; - c[i] = (BignumInt)carry; - carry >>= BIGNUM_INT_BITS; - } + for (i = len-1; i >= 0; i--) + BignumADC(c[i], carry, a[i], b[i], carry); return (BignumInt)carry; } @@ -199,13 +111,10 @@ static void internal_sub(const BignumInt *a, const BignumInt *b, BignumInt *c, int len) { int i; - BignumDblInt carry = 1; + BignumCarry carry = 1; - for (i = len-1; i >= 0; i--) { - carry += (BignumDblInt)a[i] + (b[i] ^ BIGNUM_INT_MASK); - c[i] = (BignumInt)carry; - carry >>= BIGNUM_INT_BITS; - } + for (i = len-1; i >= 0; i--) + BignumADC(c[i], carry, a[i], ~b[i], carry); } /* @@ -269,7 +178,7 @@ static void internal_mul(const BignumInt *a, const BignumInt *b, int toplen = len/2, botlen = len - toplen; /* botlen is the bigger */ int midlen = botlen + 1; - BignumDblInt carry; + BignumCarry carry; #ifdef KARA_DEBUG int i; #endif @@ -398,9 +307,7 @@ static void internal_mul(const BignumInt *a, const BignumInt *b, i = 2*len - botlen - 2*midlen - 1; while (carry) { assert(i >= 0); - carry += c[i]; - c[i] = (BignumInt)carry; - carry >>= BIGNUM_INT_BITS; + BignumADC(c[i], carry, c[i], 0, carry); i--; } #ifdef KARA_DEBUG @@ -414,7 +321,6 @@ static void internal_mul(const BignumInt *a, const BignumInt *b, } else { int i; BignumInt carry; - BignumDblInt t; const BignumInt *ap, *bp; BignumInt *cp, *cps; @@ -427,11 +333,8 @@ static void internal_mul(const BignumInt *a, const BignumInt *b, for (cps = c + 2*len, ap = a + len; ap-- > a; cps--) { carry = 0; - for (cp = cps, bp = b + len; cp--, bp-- > b ;) { - t = (MUL_WORD(*ap, *bp) + carry) + *cp; - *cp = (BignumInt) t; - carry = (BignumInt)(t >> BIGNUM_INT_BITS); - } + for (cp = cps, bp = b + len; cp--, bp-- > b ;) + BignumMULADD2(carry, *cp, *ap, *bp, *cp, carry); *cp = carry; } } @@ -516,7 +419,6 @@ static void internal_mul_low(const BignumInt *a, const BignumInt *b, } else { int i; BignumInt carry; - BignumDblInt t; const BignumInt *ap, *bp; BignumInt *cp, *cps; @@ -529,11 +431,8 @@ static void internal_mul_low(const BignumInt *a, const BignumInt *b, for (cps = c + len, ap = a + len; ap-- > a; cps--) { carry = 0; - for (cp = cps, bp = b + len; bp--, cp-- > c ;) { - t = (MUL_WORD(*ap, *bp) + carry) + *cp; - *cp = (BignumInt) t; - carry = (BignumInt)(t >> BIGNUM_INT_BITS); - } + for (cp = cps, bp = b + len; bp--, cp-- > c ;) + BignumMULADD2(carry, *cp, *ap, *bp, *cp, carry); } } } @@ -600,125 +499,440 @@ static void monty_reduce(BignumInt *x, const BignumInt *n, } static void internal_add_shifted(BignumInt *number, - unsigned n, int shift) + BignumInt n, int shift) { int word = 1 + (shift / BIGNUM_INT_BITS); int bshift = shift % BIGNUM_INT_BITS; - BignumDblInt addend; - - addend = (BignumDblInt)n << bshift; - - while (addend) { + BignumInt addendh, addendl; + BignumCarry carry; + + addendl = n << bshift; + addendh = (bshift == 0 ? 0 : n >> (BIGNUM_INT_BITS - bshift)); + + assert(word <= number[0]); + BignumADC(number[word], carry, number[word], addendl, 0); + word++; + if (!addendh && !carry) + return; + assert(word <= number[0]); + BignumADC(number[word], carry, number[word], addendh, carry); + word++; + while (carry) { assert(word <= number[0]); - addend += number[word]; - number[word] = (BignumInt) addend & BIGNUM_INT_MASK; - addend >>= BIGNUM_INT_BITS; + BignumADC(number[word], carry, number[word], 0, carry); word++; } } +static int bn_clz(BignumInt x) +{ + /* + * Count the leading zero bits in x. Equivalently, how far left + * would we need to shift x to make its top bit set? + * + * Precondition: x != 0. + */ + + /* FIXME: would be nice to put in some compiler intrinsics under + * ifdef here */ + int i, ret = 0; + for (i = BIGNUM_INT_BITS / 2; i != 0; i >>= 1) { + if ((x >> (BIGNUM_INT_BITS-i)) == 0) { + x <<= i; + ret += i; + } + } + return ret; +} + +static BignumInt reciprocal_word(BignumInt d) +{ + BignumInt dshort, recip, prodh, prodl; + int corrections; + + /* + * Input: a BignumInt value d, with its top bit set. + */ + assert(d >> (BIGNUM_INT_BITS-1) == 1); + + /* + * Output: a value, shifted to fill a BignumInt, which is strictly + * less than 1/(d+1), i.e. is an *under*-estimate (but by as + * little as possible within the constraints) of the reciprocal of + * any number whose first BIGNUM_INT_BITS bits match d. + * + * Ideally we'd like to _totally_ fill BignumInt, i.e. always + * return a value with the top bit set. Unfortunately we can't + * quite guarantee that for all inputs and also return a fixed + * exponent. So instead we take our reciprocal to be + * 2^(BIGNUM_INT_BITS*2-1) / d, so that it has the top bit clear + * only in the exceptional case where d takes exactly the maximum + * value BIGNUM_INT_MASK; in that case, the top bit is clear and + * the next bit down is set. + */ + + /* + * Start by computing a half-length version of the answer, by + * straightforward division within a BignumInt. + */ + dshort = (d >> (BIGNUM_INT_BITS/2)) + 1; + recip = (BIGNUM_TOP_BIT + dshort - 1) / dshort; + recip <<= BIGNUM_INT_BITS - BIGNUM_INT_BITS/2; + + /* + * Newton-Raphson iteration to improve that starting reciprocal + * estimate: take f(x) = d - 1/x, and then the N-R formula gives + * x_new = x - f(x)/f'(x) = x - (d-1/x)/(1/x^2) = x(2-d*x). Or, + * taking our fixed-point representation into account, take f(x) + * to be d - K/x (where K = 2^(BIGNUM_INT_BITS*2-1) as discussed + * above) and then we get (2K - d*x) * x/K. + * + * Newton-Raphson doubles the number of correct bits at every + * iteration, and the initial division above already gave us half + * the output word, so it's only worth doing one iteration. + */ + BignumMULADD(prodh, prodl, recip, d, recip); + prodl = ~prodl; + prodh = ~prodh; + { + BignumCarry c; + BignumADC(prodl, c, prodl, 1, 0); + prodh += c; + } + BignumMUL(prodh, prodl, prodh, recip); + recip = (prodh << 1) | (prodl >> (BIGNUM_INT_BITS-1)); + + /* + * Now make sure we have the best possible reciprocal estimate, + * before we return it. We might have been off by a handful either + * way - not enough to bother with any better-thought-out kind of + * correction loop. + */ + BignumMULADD(prodh, prodl, recip, d, recip); + corrections = 0; + if (prodh >= BIGNUM_TOP_BIT) { + do { + BignumCarry c = 1; + BignumADC(prodl, c, prodl, ~d, c); prodh += BIGNUM_INT_MASK + c; + recip--; + corrections++; + } while (prodh >= ((BignumInt)1 << (BIGNUM_INT_BITS-1))); + } else { + while (1) { + BignumInt newprodh, newprodl; + BignumCarry c = 0; + BignumADC(newprodl, c, prodl, d, c); newprodh = prodh + c; + if (newprodh >= BIGNUM_TOP_BIT) + break; + prodh = newprodh; + prodl = newprodl; + recip++; + corrections++; + } + } + + return recip; +} + /* * Compute a = a % m. * Input in first alen words of a and first mlen words of m. * Output in first alen words of a * (of which first alen-mlen words will be zero). - * The MSW of m MUST have its high bit set. * Quotient is accumulated in the `quotient' array, which is a Bignum - * rather than the internal bigendian format. Quotient parts are shifted - * left by `qshift' before adding into quot. + * rather than the internal bigendian format. + * + * 'recip' must be the result of calling reciprocal_word() on the top + * BIGNUM_INT_BITS of the modulus (denoted m0 in comments below), with + * the topmost set bit normalised to the MSB of the input to + * reciprocal_word. 'rshift' is how far left the top nonzero word of + * the modulus had to be shifted to set that top bit. */ static void internal_mod(BignumInt *a, int alen, BignumInt *m, int mlen, - BignumInt *quot, int qshift) + BignumInt *quot, BignumInt recip, int rshift) { - BignumInt m0, m1; - unsigned int h; int i, k; - m0 = m[0]; - assert(m0 >> (BIGNUM_INT_BITS-1) == 1); - if (mlen > 1) - m1 = m[1]; - else - m1 = 0; +#ifdef DIVISION_DEBUG + { + int d; + printf("start division, m=0x"); + for (d = 0; d < mlen; d++) + printf("%0*llx", BIGNUM_INT_BITS/4, (unsigned long long)m[d]); + printf(", recip=%#0*llx, rshift=%d\n", + BIGNUM_INT_BITS/4, (unsigned long long)recip, rshift); + } +#endif - for (i = 0; i <= alen - mlen; i++) { - BignumDblInt t; - unsigned int q, r, c, ai1; + /* + * Repeatedly use that reciprocal estimate to get a decent number + * of quotient bits, and subtract off the resulting multiple of m. + * + * Normally we expect to terminate this loop by means of finding + * out q=0 part way through, but one way in which we might not get + * that far in the first place is if the input a is actually zero, + * in which case we'll discard zero words from the front of a + * until we reach the termination condition in the for statement + * here. + */ + for (i = 0; i <= alen - mlen ;) { + BignumInt product; + BignumInt aword, q; + int shift, full_bitoffset, bitoffset, wordoffset; - if (i == 0) { - h = 0; - } else { - h = a[i - 1]; - a[i - 1] = 0; - } +#ifdef DIVISION_DEBUG + { + int d; + printf("main loop, a=0x"); + for (d = 0; d < alen; d++) + printf("%0*llx", BIGNUM_INT_BITS/4, (unsigned long long)a[d]); + printf("\n"); + } +#endif - if (i == alen - 1) - ai1 = 0; - else - ai1 = a[i + 1]; - - /* Find q = h:a[i] / m0 */ - if (h >= m0) { - /* - * Special case. - * - * To illustrate it, suppose a BignumInt is 8 bits, and - * we are dividing (say) A1:23:45:67 by A1:B2:C3. Then - * our initial division will be 0xA123 / 0xA1, which - * will give a quotient of 0x100 and a divide overflow. - * However, the invariants in this division algorithm - * are not violated, since the full number A1:23:... is - * _less_ than the quotient prefix A1:B2:... and so the - * following correction loop would have sorted it out. - * - * In this situation we set q to be the largest - * quotient we _can_ stomach (0xFF, of course). - */ - q = BIGNUM_INT_MASK; - } else { - /* Macro doesn't want an array subscript expression passed - * into it (see definition), so use a temporary. */ - BignumInt tmplo = a[i]; - DIVMOD_WORD(q, r, h, tmplo, m0); - - /* Refine our estimate of q by looking at - h:a[i]:a[i+1] / m0:m1 */ - t = MUL_WORD(m1, q); - if (t > ((BignumDblInt) r << BIGNUM_INT_BITS) + ai1) { - q--; - t -= m1; - r = (r + m0) & BIGNUM_INT_MASK; /* overflow? */ - if (r >= (BignumDblInt) m0 && - t > ((BignumDblInt) r << BIGNUM_INT_BITS) + ai1) q--; - } - } + if (a[i] == 0) { +#ifdef DIVISION_DEBUG + printf("zero word at i=%d\n", i); +#endif + i++; + continue; + } - /* Subtract q * m from a[i...] */ - c = 0; - for (k = mlen - 1; k >= 0; k--) { - t = MUL_WORD(q, m[k]); - t += c; - c = (unsigned)(t >> BIGNUM_INT_BITS); - if ((BignumInt) t > a[i + k]) - c++; - a[i + k] -= (BignumInt) t; - } + aword = a[i]; + shift = bn_clz(aword); + aword <<= shift; + if (shift > 0 && i+1 < alen) + aword |= a[i+1] >> (BIGNUM_INT_BITS - shift); - /* Add back m in case of borrow */ - if (c != h) { - t = 0; - for (k = mlen - 1; k >= 0; k--) { - t += m[k]; - t += a[i + k]; - a[i + k] = (BignumInt) t; - t = t >> BIGNUM_INT_BITS; - } - q--; - } - if (quot) - internal_add_shifted(quot, q, qshift + BIGNUM_INT_BITS * (alen - mlen - i)); + { + BignumInt unused; + BignumMUL(q, unused, recip, aword); + (void)unused; + } + +#ifdef DIVISION_DEBUG + printf("i=%d, aword=%#0*llx, shift=%d, q=%#0*llx\n", + i, BIGNUM_INT_BITS/4, (unsigned long long)aword, + shift, BIGNUM_INT_BITS/4, (unsigned long long)q); +#endif + + /* + * Work out the right bit and word offsets to use when + * subtracting q*m from a. + * + * aword was taken from a[i], which means its LSB was at bit + * position (alen-1-i) * BIGNUM_INT_BITS. But then we shifted + * it left by 'shift', so now the low bit of aword corresponds + * to bit position (alen-1-i) * BIGNUM_INT_BITS - shift, i.e. + * aword is approximately equal to a / 2^(that). + * + * m0 comes from the top word of mod, so its LSB is at bit + * position (mlen-1) * BIGNUM_INT_BITS - rshift, i.e. it can + * be considered to be m / 2^(that power). 'recip' is the + * reciprocal of m0, times 2^(BIGNUM_INT_BITS*2-1), i.e. it's + * about 2^((mlen+1) * BIGNUM_INT_BITS - rshift - 1) / m. + * + * Hence, recip * aword is approximately equal to the product + * of those, which simplifies to + * + * a/m * 2^((mlen+2+i-alen)*BIGNUM_INT_BITS + shift - rshift - 1) + * + * But we've also shifted recip*aword down by BIGNUM_INT_BITS + * to form q, so we have + * + * q ~= a/m * 2^((mlen+1+i-alen)*BIGNUM_INT_BITS + shift - rshift - 1) + * + * and hence, when we now compute q*m, it will be about + * a*2^(all that lot), i.e. the negation of that expression is + * how far left we have to shift the product q*m to make it + * approximately equal to a. + */ + full_bitoffset = -((mlen+1+i-alen)*BIGNUM_INT_BITS + shift-rshift-1); +#ifdef DIVISION_DEBUG + printf("full_bitoffset=%d\n", full_bitoffset); +#endif + + if (full_bitoffset < 0) { + /* + * If we find ourselves needing to shift q*m _right_, that + * means we've reached the bottom of the quotient. Clip q + * so that its right shift becomes zero, and if that means + * q becomes _actually_ zero, this loop is done. + */ + if (full_bitoffset <= -BIGNUM_INT_BITS) + break; + q >>= -full_bitoffset; + full_bitoffset = 0; + if (!q) + break; +#ifdef DIVISION_DEBUG + printf("now full_bitoffset=%d, q=%#0*llx\n", + full_bitoffset, BIGNUM_INT_BITS/4, (unsigned long long)q); +#endif + } + + wordoffset = full_bitoffset / BIGNUM_INT_BITS; + bitoffset = full_bitoffset % BIGNUM_INT_BITS; +#ifdef DIVISION_DEBUG + printf("wordoffset=%d, bitoffset=%d\n", wordoffset, bitoffset); +#endif + + /* wordoffset as computed above is the offset between the LSWs + * of m and a. But in fact m and a are stored MSW-first, so we + * need to adjust it to be the offset between the actual array + * indices, and flip the sign too. */ + wordoffset = alen - mlen - wordoffset; + + if (bitoffset == 0) { + BignumCarry c = 1; + BignumInt prev_hi_word = 0; + for (k = mlen - 1; wordoffset+k >= i; k--) { + BignumInt mword = k<0 ? 0 : m[k]; + BignumMULADD(prev_hi_word, product, q, mword, prev_hi_word); +#ifdef DIVISION_DEBUG + printf(" aligned sub: product word for m[%d] = %#0*llx\n", + k, BIGNUM_INT_BITS/4, + (unsigned long long)product); +#endif +#ifdef DIVISION_DEBUG + printf(" aligned sub: subtrahend for a[%d] = %#0*llx\n", + wordoffset+k, BIGNUM_INT_BITS/4, + (unsigned long long)product); +#endif + BignumADC(a[wordoffset+k], c, a[wordoffset+k], ~product, c); + } + } else { + BignumInt add_word = 0; + BignumInt c = 1; + BignumInt prev_hi_word = 0; + for (k = mlen - 1; wordoffset+k >= i; k--) { + BignumInt mword = k<0 ? 0 : m[k]; + BignumMULADD(prev_hi_word, product, q, mword, prev_hi_word); +#ifdef DIVISION_DEBUG + printf(" unaligned sub: product word for m[%d] = %#0*llx\n", + k, BIGNUM_INT_BITS/4, + (unsigned long long)product); +#endif + + add_word |= product << bitoffset; + +#ifdef DIVISION_DEBUG + printf(" unaligned sub: subtrahend for a[%d] = %#0*llx\n", + wordoffset+k, + BIGNUM_INT_BITS/4, (unsigned long long)add_word); +#endif + BignumADC(a[wordoffset+k], c, a[wordoffset+k], ~add_word, c); + + add_word = product >> (BIGNUM_INT_BITS - bitoffset); + } + } + + if (quot) { +#ifdef DIVISION_DEBUG + printf("adding quotient word %#0*llx << %d\n", + BIGNUM_INT_BITS/4, (unsigned long long)q, full_bitoffset); +#endif + internal_add_shifted(quot, q, full_bitoffset); +#ifdef DIVISION_DEBUG + { + int d; + printf("now quot=0x"); + for (d = quot[0]; d > 0; d--) + printf("%0*llx", BIGNUM_INT_BITS/4, + (unsigned long long)quot[d]); + printf("\n"); + } +#endif + } + } + +#ifdef DIVISION_DEBUG + { + int d; + printf("end main loop, a=0x"); + for (d = 0; d < alen; d++) + printf("%0*llx", BIGNUM_INT_BITS/4, (unsigned long long)a[d]); + if (quot) { + printf(", quot=0x"); + for (d = quot[0]; d > 0; d--) + printf("%0*llx", BIGNUM_INT_BITS/4, + (unsigned long long)quot[d]); + } + printf("\n"); + } +#endif + + /* + * The above loop should terminate with the remaining value in a + * being strictly less than 2*m (if a >= 2*m then we should always + * have managed to get a nonzero q word), but we can't guarantee + * that it will be strictly less than m: consider a case where the + * remainder is 1, and another where the remainder is m-1. By the + * time a contains a value that's _about m_, you clearly can't + * distinguish those cases by looking at only the top word of a - + * you have to go all the way down to the bottom before you find + * out whether it's just less or just more than m. + * + * Hence, we now do a final fixup in which we subtract one last + * copy of m, or don't, accordingly. We should never have to + * subtract more than one copy of m here. + */ + for (i = 0; i < alen; i++) { + /* Compare a with m, word by word, from the MSW down. As soon + * as we encounter a difference, we know whether we need the + * fixup. */ + int mindex = mlen-alen+i; + BignumInt mword = mindex < 0 ? 0 : m[mindex]; + if (a[i] < mword) { +#ifdef DIVISION_DEBUG + printf("final fixup not needed, a < m\n"); +#endif + return; + } else if (a[i] > mword) { +#ifdef DIVISION_DEBUG + printf("final fixup is needed, a > m\n"); +#endif + break; + } + /* If neither of those cases happened, the words are the same, + * so keep going and look at the next one. */ + } +#ifdef DIVISION_DEBUG + if (i == mlen) /* if we printed neither of the above diagnostics */ + printf("final fixup is needed, a == m\n"); +#endif + + /* + * If we got here without returning, then a >= m, so we must + * subtract m, and increment the quotient. + */ + { + BignumCarry c = 1; + for (i = alen - 1; i >= 0; i--) { + int mindex = mlen-alen+i; + BignumInt mword = mindex < 0 ? 0 : m[mindex]; + BignumADC(a[i], c, a[i], ~mword, c); + } + } + if (quot) + internal_add_shifted(quot, 1, 0); + +#ifdef DIVISION_DEBUG + { + int d; + printf("after final fixup, a=0x"); + for (d = 0; d < alen; d++) + printf("%0*llx", BIGNUM_INT_BITS/4, (unsigned long long)a[d]); + if (quot) { + printf(", quot=0x"); + for (d = quot[0]; d > 0; d--) + printf("%0*llx", BIGNUM_INT_BITS/4, + (unsigned long long)quot[d]); + } + printf("\n"); } +#endif } /* @@ -727,7 +941,8 @@ static void internal_mod(BignumInt *a, int alen, Bignum modpow_simple(Bignum base_in, Bignum exp, Bignum mod) { BignumInt *a, *b, *n, *m, *scratch; - int mshift; + BignumInt recip; + int rshift; int mlen, scratchlen, i, j; Bignum base, result; @@ -750,16 +965,6 @@ Bignum modpow_simple(Bignum base_in, Bignum exp, Bignum mod) for (j = 0; j < mlen; j++) m[j] = mod[mod[0] - j]; - /* Shift m left to make msb bit set */ - for (mshift = 0; mshift < BIGNUM_INT_BITS-1; mshift++) - if ((m[0] << mshift) & BIGNUM_TOP_BIT) - break; - if (mshift) { - for (i = 0; i < mlen - 1; i++) - m[i] = (m[i] << mshift) | (m[i + 1] >> (BIGNUM_INT_BITS - mshift)); - m[mlen - 1] = m[mlen - 1] << mshift; - } - /* Allocate n of size mlen, copy base to n */ n = snewn(mlen, BignumInt); i = mlen - base[0]; @@ -782,7 +987,7 @@ Bignum modpow_simple(Bignum base_in, Bignum exp, Bignum mod) /* Skip leading zero bits of exp. */ i = 0; j = BIGNUM_INT_BITS-1; - while (i < (int)exp[0] && (exp[exp[0] - i] & (1 << j)) == 0) { + while (i < (int)exp[0] && (exp[exp[0] - i] & ((BignumInt)1 << j)) == 0) { j--; if (j < 0) { i++; @@ -790,14 +995,26 @@ Bignum modpow_simple(Bignum base_in, Bignum exp, Bignum mod) } } + /* Compute reciprocal of the top full word of the modulus */ + { + BignumInt m0 = m[0]; + rshift = bn_clz(m0); + if (rshift) { + m0 <<= rshift; + if (mlen > 1) + m0 |= m[1] >> (BIGNUM_INT_BITS - rshift); + } + recip = reciprocal_word(m0); + } + /* Main computation */ while (i < (int)exp[0]) { while (j >= 0) { internal_mul(a + mlen, a + mlen, b, mlen, scratch); - internal_mod(b, mlen * 2, m, mlen, NULL, 0); - if ((exp[exp[0] - i] & (1 << j)) != 0) { + internal_mod(b, mlen * 2, m, mlen, NULL, recip, rshift); + if ((exp[exp[0] - i] & ((BignumInt)1 << j)) != 0) { internal_mul(b + mlen, n, a, mlen, scratch); - internal_mod(a, mlen * 2, m, mlen, NULL, 0); + internal_mod(a, mlen * 2, m, mlen, NULL, recip, rshift); } else { BignumInt *t; t = a; @@ -810,16 +1027,6 @@ Bignum modpow_simple(Bignum base_in, Bignum exp, Bignum mod) j = BIGNUM_INT_BITS-1; } - /* Fixup result in case the modulus was shifted */ - if (mshift) { - for (i = mlen - 1; i < 2 * mlen - 1; i++) - a[i] = (a[i] << mshift) | (a[i + 1] >> (BIGNUM_INT_BITS - mshift)); - a[2 * mlen - 1] = a[2 * mlen - 1] << mshift; - internal_mod(a, mlen * 2, m, mlen, NULL, 0); - for (i = 2 * mlen - 1; i >= mlen; i--) - a[i] = (a[i] >> mshift) | (a[i - 1] << (BIGNUM_INT_BITS - mshift)); - } - /* Copy result to buffer */ result = newbn(mod[0]); for (i = 0; i < mlen; i++) @@ -932,7 +1139,7 @@ Bignum modpow(Bignum base_in, Bignum exp, Bignum mod) /* Skip leading zero bits of exp. */ i = 0; j = BIGNUM_INT_BITS-1; - while (i < (int)exp[0] && (exp[exp[0] - i] & (1 << j)) == 0) { + while (i < (int)exp[0] && (exp[exp[0] - i] & ((BignumInt)1 << j)) == 0) { j--; if (j < 0) { i++; @@ -945,7 +1152,7 @@ Bignum modpow(Bignum base_in, Bignum exp, Bignum mod) while (j >= 0) { internal_mul(a + len, a + len, b, len, scratch); monty_reduce(b, n, mninv, scratch, len); - if ((exp[exp[0] - i] & (1 << j)) != 0) { + if ((exp[exp[0] - i] & ((BignumInt)1 << j)) != 0) { internal_mul(b + len, x, a, len, scratch); monty_reduce(a, n, mninv, scratch, len); } else { @@ -998,7 +1205,8 @@ Bignum modpow(Bignum base_in, Bignum exp, Bignum mod) Bignum modmul(Bignum p, Bignum q, Bignum mod) { BignumInt *a, *n, *m, *o, *scratch; - int mshift, scratchlen; + BignumInt recip; + int rshift, scratchlen; int pqlen, mlen, rlen, i, j; Bignum result; @@ -1015,16 +1223,6 @@ Bignum modmul(Bignum p, Bignum q, Bignum mod) for (j = 0; j < mlen; j++) m[j] = mod[mod[0] - j]; - /* Shift m left to make msb bit set */ - for (mshift = 0; mshift < BIGNUM_INT_BITS-1; mshift++) - if ((m[0] << mshift) & BIGNUM_TOP_BIT) - break; - if (mshift) { - for (i = 0; i < mlen - 1; i++) - m[i] = (m[i] << mshift) | (m[i + 1] >> (BIGNUM_INT_BITS - mshift)); - m[mlen - 1] = m[mlen - 1] << mshift; - } - pqlen = (p[0] > q[0] ? p[0] : q[0]); /* @@ -1057,19 +1255,21 @@ Bignum modmul(Bignum p, Bignum q, Bignum mod) scratchlen = mul_compute_scratch(pqlen); scratch = snewn(scratchlen, BignumInt); + /* Compute reciprocal of the top full word of the modulus */ + { + BignumInt m0 = m[0]; + rshift = bn_clz(m0); + if (rshift) { + m0 <<= rshift; + if (mlen > 1) + m0 |= m[1] >> (BIGNUM_INT_BITS - rshift); + } + recip = reciprocal_word(m0); + } + /* Main computation */ internal_mul(n, o, a, pqlen, scratch); - internal_mod(a, pqlen * 2, m, mlen, NULL, 0); - - /* Fixup result in case the modulus was shifted */ - if (mshift) { - for (i = 2 * pqlen - mlen - 1; i < 2 * pqlen - 1; i++) - a[i] = (a[i] << mshift) | (a[i + 1] >> (BIGNUM_INT_BITS - mshift)); - a[2 * pqlen - 1] = a[2 * pqlen - 1] << mshift; - internal_mod(a, pqlen * 2, m, mlen, NULL, 0); - for (i = 2 * pqlen - 1; i >= 2 * pqlen - mlen; i--) - a[i] = (a[i] >> mshift) | (a[i - 1] << (BIGNUM_INT_BITS - mshift)); - } + internal_mod(a, pqlen * 2, m, mlen, NULL, recip, rshift); /* Copy result to buffer */ rlen = (mlen < pqlen * 2 ? mlen : pqlen * 2); @@ -1112,12 +1312,9 @@ Bignum modsub(const Bignum a, const Bignum b, const Bignum n) /* Handle going round the corner of the modulus without having * negative support in Bignum */ Bignum tmp = bigsub(n, b1); - if (tmp) { - ret = bigadd(tmp, a1); - freebn(tmp); - } else { - ret = NULL; - } + assert(tmp); + ret = bigadd(tmp, a1); + freebn(tmp); } if (a != a1) freebn(a1); @@ -1136,7 +1333,8 @@ Bignum modsub(const Bignum a, const Bignum b, const Bignum n) static void bigdivmod(Bignum p, Bignum mod, Bignum result, Bignum quotient) { BignumInt *n, *m; - int mshift; + BignumInt recip; + int rshift; int plen, mlen, i, j; /* @@ -1152,16 +1350,6 @@ static void bigdivmod(Bignum p, Bignum mod, Bignum result, Bignum quotient) for (j = 0; j < mlen; j++) m[j] = mod[mod[0] - j]; - /* Shift m left to make msb bit set */ - for (mshift = 0; mshift < BIGNUM_INT_BITS-1; mshift++) - if ((m[0] << mshift) & BIGNUM_TOP_BIT) - break; - if (mshift) { - for (i = 0; i < mlen - 1; i++) - m[i] = (m[i] << mshift) | (m[i + 1] >> (BIGNUM_INT_BITS - mshift)); - m[mlen - 1] = m[mlen - 1] << mshift; - } - plen = p[0]; /* Ensure plen > mlen */ if (plen <= mlen) @@ -1174,19 +1362,21 @@ static void bigdivmod(Bignum p, Bignum mod, Bignum result, Bignum quotient) for (j = 1; j <= (int)p[0]; j++) n[plen - j] = p[j]; - /* Main computation */ - internal_mod(n, plen, m, mlen, quotient, mshift); - - /* Fixup result in case the modulus was shifted */ - if (mshift) { - for (i = plen - mlen - 1; i < plen - 1; i++) - n[i] = (n[i] << mshift) | (n[i + 1] >> (BIGNUM_INT_BITS - mshift)); - n[plen - 1] = n[plen - 1] << mshift; - internal_mod(n, plen, m, mlen, quotient, 0); - for (i = plen - 1; i >= plen - mlen; i--) - n[i] = (n[i] >> mshift) | (n[i - 1] << (BIGNUM_INT_BITS - mshift)); + /* Compute reciprocal of the top full word of the modulus */ + { + BignumInt m0 = m[0]; + rshift = bn_clz(m0); + if (rshift) { + m0 <<= rshift; + if (mlen > 1) + m0 |= m[1] >> (BIGNUM_INT_BITS - rshift); + } + recip = reciprocal_word(m0); } + /* Main computation */ + internal_mod(n, plen, m, mlen, quotient, recip, rshift); + /* Copy result to buffer */ if (result) { for (i = 1; i <= (int)result[0]; i++) { @@ -1227,11 +1417,11 @@ Bignum bignum_from_bytes(const unsigned char *data, int nbytes) result[i] = 0; for (i = nbytes; i--;) { unsigned char byte = *data++; - result[1 + i / BIGNUM_INT_BYTES] |= byte << (8*i % BIGNUM_INT_BITS); + result[1 + i / BIGNUM_INT_BYTES] |= + (BignumInt)byte << (8*i % BIGNUM_INT_BITS); } - while (result[0] > 1 && result[result[0]] == 0) - result[0]--; + bn_restore_invariant(result); return result; } @@ -1249,11 +1439,11 @@ Bignum bignum_from_bytes_le(const unsigned char *data, int nbytes) result[i] = 0; for (i = 0; i < nbytes; ++i) { unsigned char byte = *data++; - result[1 + i / BIGNUM_INT_BYTES] |= byte << (8*i % BIGNUM_INT_BITS); + result[1 + i / BIGNUM_INT_BYTES] |= + (BignumInt)byte << (8*i % BIGNUM_INT_BITS); } - while (result[0] > 1 && result[result[0]] == 0) - result[0]--; + bn_restore_invariant(result); return result; } @@ -1271,6 +1461,7 @@ Bignum bignum_from_decimal(const char *decimal) tmp = bigmul(result, Ten); tmp2 = bignum_from_long(*decimal - '0'); + freebn(result); result = bigadd(tmp, tmp2); freebn(tmp); freebn(tmp2); @@ -1398,11 +1589,11 @@ int bignum_bit(Bignum bn, int i) */ void bignum_set_bit(Bignum bn, int bitnum, int value) { - if (bitnum < 0 || bitnum >= (int)(BIGNUM_INT_BITS * bn[0])) - abort(); /* beyond the end */ - else { + if (bitnum < 0 || bitnum >= (int)(BIGNUM_INT_BITS * bn[0])) { + if (value) abort(); /* beyond the end */ + } else { int v = bitnum / BIGNUM_INT_BITS + 1; - int mask = 1 << (bitnum % BIGNUM_INT_BITS); + BignumInt mask = (BignumInt)1 << (bitnum % BIGNUM_INT_BITS); if (value) bn[v] |= mask; else @@ -1500,7 +1691,6 @@ Bignum bignum_lshift(Bignum a, int shift) bits = bignum_bitcount(a) + shift; ret = newbn((bits + BIGNUM_INT_BITS - 1) / BIGNUM_INT_BITS); - if (!ret) return NULL; shiftWords = shift / BIGNUM_INT_BITS; shiftBits = shift % BIGNUM_INT_BITS; @@ -1566,12 +1756,11 @@ Bignum bigmuladd(Bignum a, Bignum b, Bignum addend) /* now add in the addend, if any */ if (addend) { - BignumDblInt carry = 0; + BignumCarry carry = 0; for (i = 1; i <= rlen; i++) { - carry += (i <= (int)ret[0] ? ret[i] : 0); - carry += (i <= (int)addend[0] ? addend[i] : 0); - ret[i] = (BignumInt) carry & BIGNUM_INT_MASK; - carry >>= BIGNUM_INT_BITS; + BignumInt retword = (i <= (int)ret[0] ? ret[i] : 0); + BignumInt addword = (i <= (int)addend[0] ? addend[i] : 0); + BignumADC(ret[i], carry, retword, addword, carry); if (ret[i] != 0 && i > maxspot) maxspot = i; } @@ -1600,17 +1789,16 @@ Bignum bigadd(Bignum a, Bignum b) int rlen = (alen > blen ? alen : blen) + 1; int i, maxspot; Bignum ret; - BignumDblInt carry; + BignumCarry carry; ret = newbn(rlen); carry = 0; maxspot = 0; for (i = 1; i <= rlen; i++) { - carry += (i <= (int)a[0] ? a[i] : 0); - carry += (i <= (int)b[0] ? b[i] : 0); - ret[i] = (BignumInt) carry & BIGNUM_INT_MASK; - carry >>= BIGNUM_INT_BITS; + BignumInt aword = (i <= (int)a[0] ? a[i] : 0); + BignumInt bword = (i <= (int)b[0] ? b[i] : 0); + BignumADC(ret[i], carry, aword, bword, carry); if (ret[i] != 0 && i > maxspot) maxspot = i; } @@ -1630,17 +1818,16 @@ Bignum bigsub(Bignum a, Bignum b) int rlen = (alen > blen ? alen : blen); int i, maxspot; Bignum ret; - BignumDblInt carry; + BignumCarry carry; ret = newbn(rlen); carry = 1; maxspot = 0; for (i = 1; i <= rlen; i++) { - carry += (i <= (int)a[0] ? a[i] : 0); - carry += (i <= (int)b[0] ? b[i] ^ BIGNUM_INT_MASK : BIGNUM_INT_MASK); - ret[i] = (BignumInt) carry & BIGNUM_INT_MASK; - carry >>= BIGNUM_INT_BITS; + BignumInt aword = (i <= (int)a[0] ? a[i] : 0); + BignumInt bword = (i <= (int)b[0] ? b[i] : 0); + BignumADC(ret[i], carry, aword, ~bword, carry); if (ret[i] != 0 && i > maxspot) maxspot = i; } @@ -1680,40 +1867,52 @@ Bignum bignum_bitmask(Bignum n) } /* - * Convert a (max 32-bit) long into a bignum. + * Convert an unsigned long into a bignum. */ -Bignum bignum_from_long(unsigned long nn) +Bignum bignum_from_long(unsigned long n) { + const int maxwords = + (sizeof(unsigned long) + sizeof(BignumInt) - 1) / sizeof(BignumInt); Bignum ret; - BignumDblInt n = nn; + int i; + + ret = newbn(maxwords); + ret[0] = 0; + for (i = 0; i < maxwords; i++) { + ret[i+1] = n >> (i * BIGNUM_INT_BITS); + if (ret[i+1] != 0) + ret[0] = i+1; + } - ret = newbn(3); - ret[1] = (BignumInt)(n & BIGNUM_INT_MASK); - ret[2] = (BignumInt)((n >> BIGNUM_INT_BITS) & BIGNUM_INT_MASK); - ret[3] = 0; - ret[0] = (ret[2] ? 2 : 1); return ret; } /* * Add a long to a bignum. */ -Bignum bignum_add_long(Bignum number, unsigned long addendx) +Bignum bignum_add_long(Bignum number, unsigned long n) { - Bignum ret = newbn(number[0] + 1); - int i, maxspot = 0; - BignumDblInt carry = 0, addend = addendx; + const int maxwords = + (sizeof(unsigned long) + sizeof(BignumInt) - 1) / sizeof(BignumInt); + Bignum ret; + int words, i; + BignumCarry carry; - for (i = 1; i <= (int)ret[0]; i++) { - carry += addend & BIGNUM_INT_MASK; - carry += (i <= (int)number[0] ? number[i] : 0); - addend >>= BIGNUM_INT_BITS; - ret[i] = (BignumInt) carry & BIGNUM_INT_MASK; - carry >>= BIGNUM_INT_BITS; - if (ret[i] != 0) - maxspot = i; + words = number[0]; + if (words < maxwords) + words = maxwords; + words++; + ret = newbn(words); + + carry = 0; + ret[0] = 0; + for (i = 0; i < words; i++) { + BignumInt nword = (i < maxwords ? n >> (i * BIGNUM_INT_BITS) : 0); + BignumInt numword = (i < number[0] ? number[i+1] : 0); + BignumADC(ret[i+1], carry, numword, nword, carry); + if (ret[i+1] != 0) + ret[0] = i+1; } - ret[0] = maxspot; return ret; } @@ -1722,13 +1921,17 @@ Bignum bignum_add_long(Bignum number, unsigned long addendx) */ unsigned short bignum_mod_short(Bignum number, unsigned short modulus) { - BignumDblInt mod, r; + unsigned long mod = modulus, r = 0; + /* Precompute (BIGNUM_INT_MASK+1) % mod */ + unsigned long base_r = (BIGNUM_INT_MASK - modulus + 1) % mod; int i; - r = 0; - mod = modulus; - for (i = number[0]; i > 0; i--) - r = (r * (BIGNUM_TOP_BIT % mod) * 2 + number[i] % mod) % mod; + for (i = number[0]; i > 0; i--) { + /* + * Conceptually, ((r << BIGNUM_INT_BITS) + number[i]) % mod + */ + r = ((r * base_r) + (number[i] % mod)) % mod; + } return (unsigned short) r; } @@ -1886,7 +2089,7 @@ char *bignum_decimal(Bignum x) { int ndigits, ndigit; int i, iszero; - BignumDblInt carry; + BignumInt carry; char *ret; BignumInt *workspace; @@ -1933,11 +2136,33 @@ char *bignum_decimal(Bignum x) iszero = 1; carry = 0; for (i = 0; i < (int)x[0]; i++) { - carry = (carry << BIGNUM_INT_BITS) + workspace[i]; - workspace[i] = (BignumInt) (carry / 10); + /* + * Conceptually, we want to compute + * + * (carry << BIGNUM_INT_BITS) + workspace[i] + * ----------------------------------------- + * 10 + * + * but we don't have an integer type longer than BignumInt + * to work with. So we have to do it in pieces. + */ + + BignumInt q, r; + q = workspace[i] / 10; + r = workspace[i] % 10; + + /* I want (BIGNUM_INT_MASK+1)/10 but can't say so directly! */ + q += carry * ((BIGNUM_INT_MASK-9) / 10 + 1); + r += carry * ((BIGNUM_INT_MASK-9) % 10); + + q += r / 10; + r %= 10; + + workspace[i] = q; + carry = r; + if (workspace[i]) iszero = 0; - carry %= 10; } ret[--ndigit] = (char) (carry + '0'); } while (!iszero); @@ -1956,204 +2181,3 @@ char *bignum_decimal(Bignum x) sfree(workspace); return ret; } - -#ifdef TESTBN - -#include -#include -#include - -/* - * gcc -Wall -g -O0 -DTESTBN -o testbn sshbn.c misc.c conf.c tree234.c unix/uxmisc.c -I. -I unix -I charset - * - * Then feed to this program's standard input the output of - * testdata/bignum.py . - */ - -void modalfatalbox(const char *p, ...) -{ - va_list ap; - fprintf(stderr, "FATAL ERROR: "); - va_start(ap, p); - vfprintf(stderr, p, ap); - va_end(ap); - fputc('\n', stderr); - exit(1); -} - -#define fromxdigit(c) ( (c)>'9' ? ((c)&0xDF) - 'A' + 10 : (c) - '0' ) - -int main(int argc, char **argv) -{ - char *buf; - int line = 0; - int passes = 0, fails = 0; - - while ((buf = fgetline(stdin)) != NULL) { - int maxlen = strlen(buf); - unsigned char *data = snewn(maxlen, unsigned char); - unsigned char *ptrs[5], *q; - int ptrnum; - char *bufp = buf; - - line++; - - q = data; - ptrnum = 0; - - while (*bufp && !isspace((unsigned char)*bufp)) - bufp++; - if (bufp) - *bufp++ = '\0'; - - while (*bufp) { - char *start, *end; - int i; - - while (*bufp && !isxdigit((unsigned char)*bufp)) - bufp++; - start = bufp; - - if (!*bufp) - break; - - while (*bufp && isxdigit((unsigned char)*bufp)) - bufp++; - end = bufp; - - if (ptrnum >= lenof(ptrs)) - break; - ptrs[ptrnum++] = q; - - for (i = -((end - start) & 1); i < end-start; i += 2) { - unsigned char val = (i < 0 ? 0 : fromxdigit(start[i])); - val = val * 16 + fromxdigit(start[i+1]); - *q++ = val; - } - - ptrs[ptrnum] = q; - } - - if (!strcmp(buf, "mul")) { - Bignum a, b, c, p; - - if (ptrnum != 3) { - printf("%d: mul with %d parameters, expected 3\n", line, ptrnum); - exit(1); - } - a = bignum_from_bytes(ptrs[0], ptrs[1]-ptrs[0]); - b = bignum_from_bytes(ptrs[1], ptrs[2]-ptrs[1]); - c = bignum_from_bytes(ptrs[2], ptrs[3]-ptrs[2]); - p = bigmul(a, b); - - if (bignum_cmp(c, p) == 0) { - passes++; - } else { - char *as = bignum_decimal(a); - char *bs = bignum_decimal(b); - char *cs = bignum_decimal(c); - char *ps = bignum_decimal(p); - - printf("%d: fail: %s * %s gave %s expected %s\n", - line, as, bs, ps, cs); - fails++; - - sfree(as); - sfree(bs); - sfree(cs); - sfree(ps); - } - freebn(a); - freebn(b); - freebn(c); - freebn(p); - } else if (!strcmp(buf, "modmul")) { - Bignum a, b, m, c, p; - - if (ptrnum != 4) { - printf("%d: modmul with %d parameters, expected 4\n", - line, ptrnum); - exit(1); - } - a = bignum_from_bytes(ptrs[0], ptrs[1]-ptrs[0]); - b = bignum_from_bytes(ptrs[1], ptrs[2]-ptrs[1]); - m = bignum_from_bytes(ptrs[2], ptrs[3]-ptrs[2]); - c = bignum_from_bytes(ptrs[3], ptrs[4]-ptrs[3]); - p = modmul(a, b, m); - - if (bignum_cmp(c, p) == 0) { - passes++; - } else { - char *as = bignum_decimal(a); - char *bs = bignum_decimal(b); - char *ms = bignum_decimal(m); - char *cs = bignum_decimal(c); - char *ps = bignum_decimal(p); - - printf("%d: fail: %s * %s mod %s gave %s expected %s\n", - line, as, bs, ms, ps, cs); - fails++; - - sfree(as); - sfree(bs); - sfree(ms); - sfree(cs); - sfree(ps); - } - freebn(a); - freebn(b); - freebn(m); - freebn(c); - freebn(p); - } else if (!strcmp(buf, "pow")) { - Bignum base, expt, modulus, expected, answer; - - if (ptrnum != 4) { - printf("%d: mul with %d parameters, expected 4\n", line, ptrnum); - exit(1); - } - - base = bignum_from_bytes(ptrs[0], ptrs[1]-ptrs[0]); - expt = bignum_from_bytes(ptrs[1], ptrs[2]-ptrs[1]); - modulus = bignum_from_bytes(ptrs[2], ptrs[3]-ptrs[2]); - expected = bignum_from_bytes(ptrs[3], ptrs[4]-ptrs[3]); - answer = modpow(base, expt, modulus); - - if (bignum_cmp(expected, answer) == 0) { - passes++; - } else { - char *as = bignum_decimal(base); - char *bs = bignum_decimal(expt); - char *cs = bignum_decimal(modulus); - char *ds = bignum_decimal(answer); - char *ps = bignum_decimal(expected); - - printf("%d: fail: %s ^ %s mod %s gave %s expected %s\n", - line, as, bs, cs, ds, ps); - fails++; - - sfree(as); - sfree(bs); - sfree(cs); - sfree(ds); - sfree(ps); - } - freebn(base); - freebn(expt); - freebn(modulus); - freebn(expected); - freebn(answer); - } else { - printf("%d: unrecognised test keyword: '%s'\n", line, buf); - exit(1); - } - - sfree(buf); - sfree(data); - } - - printf("passed %d failed %d total %d\n", passes, fails, passes+fails); - return fails != 0; -} - -#endif