]> asedeno.scripts.mit.edu Git - PuTTY.git/blobdiff - sshbn.c
Fix copy-and-paste error in testbn main program.
[PuTTY.git] / sshbn.c
diff --git a/sshbn.c b/sshbn.c
index 9d69b38588aba720ae1856dde3b308e9dc8f3e41..afc38a0bb2b5dfd056a50fae8483fa7dd020e6ba 100644 (file)
--- a/sshbn.c
+++ b/sshbn.c
@@ -6,93 +6,12 @@
 #include <assert.h>
 #include <stdlib.h>
 #include <string.h>
+#include <limits.h>
+#include <ctype.h>
 
 #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.
- * <http://msdn2.microsoft.com/en-us/library/bf1dw62z.aspx> */
-#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;
@@ -101,6 +20,7 @@ typedef BignumInt *Bignum;
 
 BignumInt bnZero[1] = { 0 };
 BignumInt bnOne[2] = { 1, 1 };
+BignumInt bnTen[2] = { 1, 10 };
 
 /*
  * The Bignum format is an array of `BignumInt'. The first
@@ -116,13 +36,15 @@ BignumInt bnOne[2] = { 1, 1 };
  * nonzero.
  */
 
-Bignum Zero = bnZero, One = bnOne;
+Bignum Zero = bnZero, One = bnOne, Ten = bnTen;
 
 static Bignum newbn(int length)
 {
-    Bignum b = snewn(length + 1, BignumInt);
-    if (!b)
-       abort();                       /* FIXME */
+    Bignum b;
+
+    assert(length >= 0 && length < INT_MAX / BIGNUM_INT_BITS);
+
+    b = snewn(length + 1, BignumInt);
     memset(b, 0, (length + 1) * sizeof(*b));
     b[0] = length;
     return b;
@@ -154,7 +76,11 @@ void freebn(Bignum b)
 
 Bignum bn_power_2(int n)
 {
-    Bignum ret = newbn(n / BIGNUM_INT_BITS + 1);
+    Bignum ret;
+
+    assert(n >= 0);
+
+    ret = newbn(n / BIGNUM_INT_BITS + 1);
     bignum_set_bit(ret, n, 1);
     return ret;
 }
@@ -589,7 +515,7 @@ 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;
@@ -598,6 +524,7 @@ static void internal_add_shifted(BignumInt *number,
     addend = (BignumDblInt)n << bshift;
 
     while (addend) {
+        assert(word <= number[0]);
        addend += number[word];
        number[word] = (BignumInt) addend & BIGNUM_INT_MASK;
        addend >>= BIGNUM_INT_BITS;
@@ -619,11 +546,11 @@ static void internal_mod(BignumInt *a, int alen,
                         BignumInt *m, int mlen,
                         BignumInt *quot, int qshift)
 {
-    BignumInt m0, m1;
-    unsigned int h;
+    BignumInt m0, m1, h;
     int i, k;
 
     m0 = m[0];
+    assert(m0 >> (BIGNUM_INT_BITS-1) == 1);
     if (mlen > 1)
        m1 = m[1];
     else
@@ -631,7 +558,7 @@ static void internal_mod(BignumInt *a, int alen,
 
     for (i = 0; i <= alen - mlen; i++) {
        BignumDblInt t;
-       unsigned int q, r, c, ai1;
+        BignumInt q, r, c, ai1;
 
        if (i == 0) {
            h = 0;
@@ -686,7 +613,7 @@ static void internal_mod(BignumInt *a, int alen,
        for (k = mlen - 1; k >= 0; k--) {
            t = MUL_WORD(q, m[k]);
            t += c;
-           c = (unsigned)(t >> BIGNUM_INT_BITS);
+           c = (BignumInt)(t >> BIGNUM_INT_BITS);
            if ((BignumInt) t > a[i + k])
                c++;
            a[i + k] -= (BignumInt) t;
@@ -769,7 +696,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++;
@@ -782,7 +709,7 @@ Bignum modpow_simple(Bignum base_in, Bignum exp, Bignum mod)
        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) {
+           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);
            } else {
@@ -815,20 +742,15 @@ Bignum modpow_simple(Bignum base_in, Bignum exp, Bignum mod)
        result[0]--;
 
     /* Free temporary arrays */
-    for (i = 0; i < 2 * mlen; i++)
-       a[i] = 0;
+    smemclr(a, 2 * mlen * sizeof(*a));
     sfree(a);
-    for (i = 0; i < scratchlen; i++)
-       scratch[i] = 0;
+    smemclr(scratch, scratchlen * sizeof(*scratch));
     sfree(scratch);
-    for (i = 0; i < 2 * mlen; i++)
-       b[i] = 0;
+    smemclr(b, 2 * mlen * sizeof(*b));
     sfree(b);
-    for (i = 0; i < mlen; i++)
-       m[i] = 0;
+    smemclr(m, mlen * sizeof(*m));
     sfree(m);
-    for (i = 0; i < mlen; i++)
-       n[i] = 0;
+    smemclr(n, mlen * sizeof(*n));
     sfree(n);
 
     freebn(base);
@@ -873,6 +795,7 @@ Bignum modpow(Bignum base_in, Bignum exp, Bignum mod)
     len = mod[0];
     r = bn_power_2(BIGNUM_INT_BITS * len);
     inv = modinv(mod, r);
+    assert(inv); /* cannot fail, since mod is odd and r is a power of 2 */
 
     /*
      * Multiply the base by r mod n, to get it into Montgomery
@@ -923,7 +846,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++;
@@ -936,7 +859,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 {
@@ -965,23 +888,17 @@ Bignum modpow(Bignum base_in, Bignum exp, Bignum mod)
        result[0]--;
 
     /* Free temporary arrays */
-    for (i = 0; i < scratchlen; i++)
-       scratch[i] = 0;
+    smemclr(scratch, scratchlen * sizeof(*scratch));
     sfree(scratch);
-    for (i = 0; i < 2 * len; i++)
-       a[i] = 0;
+    smemclr(a, 2 * len * sizeof(*a));
     sfree(a);
-    for (i = 0; i < 2 * len; i++)
-       b[i] = 0;
+    smemclr(b, 2 * len * sizeof(*b));
     sfree(b);
-    for (i = 0; i < len; i++)
-       mninv[i] = 0;
+    smemclr(mninv, len * sizeof(*mninv));
     sfree(mninv);
-    for (i = 0; i < len; i++)
-       n[i] = 0;
+    smemclr(n, len * sizeof(*n));
     sfree(n);
-    for (i = 0; i < len; i++)
-       x[i] = 0;
+    smemclr(x, len * sizeof(*x));
     sfree(x);
 
     return result;
@@ -999,6 +916,12 @@ Bignum modmul(Bignum p, Bignum q, Bignum mod)
     int pqlen, mlen, rlen, i, j;
     Bignum result;
 
+    /*
+     * The most significant word of mod needs to be non-zero. It
+     * should already be, but let's make sure.
+     */
+    assert(mod[mod[0]] != 0);
+
     /* Allocate m of size mlen, copy mod to m */
     /* We use big endian internally */
     mlen = mod[0];
@@ -1018,6 +941,13 @@ Bignum modmul(Bignum p, Bignum q, Bignum mod)
 
     pqlen = (p[0] > q[0] ? p[0] : q[0]);
 
+    /*
+     * Make sure that we're allowing enough space. The shifting below
+     * will underflow the vectors we allocate if pqlen is too small.
+     */
+    if (2*pqlen <= mlen)
+        pqlen = mlen/2 + 1;
+
     /* Allocate n of size pqlen, copy p to n */
     n = snewn(pqlen, BignumInt);
     i = pqlen - p[0];
@@ -1064,25 +994,49 @@ Bignum modmul(Bignum p, Bignum q, Bignum mod)
        result[0]--;
 
     /* Free temporary arrays */
-    for (i = 0; i < scratchlen; i++)
-       scratch[i] = 0;
+    smemclr(scratch, scratchlen * sizeof(*scratch));
     sfree(scratch);
-    for (i = 0; i < 2 * pqlen; i++)
-       a[i] = 0;
+    smemclr(a, 2 * pqlen * sizeof(*a));
     sfree(a);
-    for (i = 0; i < mlen; i++)
-       m[i] = 0;
+    smemclr(m, mlen * sizeof(*m));
     sfree(m);
-    for (i = 0; i < pqlen; i++)
-       n[i] = 0;
+    smemclr(n, pqlen * sizeof(*n));
     sfree(n);
-    for (i = 0; i < pqlen; i++)
-       o[i] = 0;
+    smemclr(o, pqlen * sizeof(*o));
     sfree(o);
 
     return result;
 }
 
+Bignum modsub(const Bignum a, const Bignum b, const Bignum n)
+{
+    Bignum a1, b1, ret;
+
+    if (bignum_cmp(a, n) >= 0) a1 = bigmod(a, n);
+    else a1 = a;
+    if (bignum_cmp(b, n) >= 0) b1 = bigmod(b, n);
+    else b1 = b;
+
+    if (bignum_cmp(a1, b1) >= 0) /* a >= b */
+    {
+        ret = bigsub(a1, b1);
+    }
+    else
+    {
+        /* Handle going round the corner of the modulus without having
+         * negative support in Bignum */
+        Bignum tmp = bigsub(n, b1);
+        assert(tmp);
+        ret = bigadd(tmp, a1);
+        freebn(tmp);
+    }
+
+    if (a != a1) freebn(a1);
+    if (b != b1) freebn(b1);
+
+    return ret;
+}
+
 /*
  * Compute p % mod.
  * The most significant word of mod MUST be non-zero.
@@ -1096,6 +1050,12 @@ static void bigdivmod(Bignum p, Bignum mod, Bignum result, Bignum quotient)
     int mshift;
     int plen, mlen, i, j;
 
+    /*
+     * The most significant word of mod needs to be non-zero. It
+     * should already be, but let's make sure.
+     */
+    assert(mod[mod[0]] != 0);
+
     /* Allocate m of size mlen, copy mod to m */
     /* We use big endian internally */
     mlen = mod[0];
@@ -1147,11 +1107,9 @@ static void bigdivmod(Bignum p, Bignum mod, Bignum result, Bignum quotient)
     }
 
     /* Free temporary arrays */
-    for (i = 0; i < mlen; i++)
-       m[i] = 0;
+    smemclr(m, mlen * sizeof(*m));
     sfree(m);
-    for (i = 0; i < plen; i++)
-       n[i] = 0;
+    smemclr(n, plen * sizeof(*n));
     sfree(n);
 }
 
@@ -1171,6 +1129,8 @@ Bignum bignum_from_bytes(const unsigned char *data, int nbytes)
     Bignum result;
     int w, i;
 
+    assert(nbytes >= 0 && nbytes < INT_MAX/8);
+
     w = (nbytes + BIGNUM_INT_BYTES - 1) / BIGNUM_INT_BYTES; /* bytes->words */
 
     result = newbn(w);
@@ -1178,14 +1138,93 @@ 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);
+    }
+
+    bn_restore_invariant(result);
+    return result;
+}
+
+Bignum bignum_from_bytes_le(const unsigned char *data, int nbytes)
+{
+    Bignum result;
+    int w, i;
+
+    assert(nbytes >= 0 && nbytes < INT_MAX/8);
+
+    w = (nbytes + BIGNUM_INT_BYTES - 1) / BIGNUM_INT_BYTES; /* bytes->words */
+
+    result = newbn(w);
+    for (i = 1; i <= w; i++)
+        result[i] = 0;
+    for (i = 0; i < nbytes; ++i) {
+        unsigned char byte = *data++;
+        result[1 + i / BIGNUM_INT_BYTES] |=
+            (BignumInt)byte << (8*i % BIGNUM_INT_BITS);
+    }
+
+    bn_restore_invariant(result);
+    return result;
+}
+
+Bignum bignum_from_decimal(const char *decimal)
+{
+    Bignum result = copybn(Zero);
+
+    while (*decimal) {
+        Bignum tmp, tmp2;
+
+        if (!isdigit((unsigned char)*decimal)) {
+            freebn(result);
+            return 0;
+        }
+
+        tmp = bigmul(result, Ten);
+        tmp2 = bignum_from_long(*decimal - '0');
+        result = bigadd(tmp, tmp2);
+        freebn(tmp);
+        freebn(tmp2);
+
+        decimal++;
     }
 
-    while (result[0] > 1 && result[result[0]] == 0)
-       result[0]--;
     return result;
 }
 
+Bignum bignum_random_in_range(const Bignum lower, const Bignum upper)
+{
+    Bignum ret = NULL;
+    unsigned char *bytes;
+    int upper_len = bignum_bitcount(upper);
+    int upper_bytes = upper_len / 8;
+    int upper_bits = upper_len % 8;
+    if (upper_bits) ++upper_bytes;
+
+    bytes = snewn(upper_bytes, unsigned char);
+    do {
+        int i;
+
+        if (ret) freebn(ret);
+
+        for (i = 0; i < upper_bytes; ++i)
+        {
+            bytes[i] = (unsigned char)random_byte();
+        }
+        /* Mask the top to reduce failure rate to 50/50 */
+        if (upper_bits)
+        {
+            bytes[i - 1] &= 0xFF >> (8 - upper_bits);
+        }
+
+        ret = bignum_from_bytes(bytes, upper_bytes);
+    } while (bignum_cmp(ret, lower) < 0 || bignum_cmp(ret, upper) > 0);
+    smemclr(bytes, upper_bytes);
+    sfree(bytes);
+
+    return ret;
+}
+
 /*
  * Read an SSH-1-format bignum from a data buffer. Return the number
  * of bytes consumed, or -1 if there wasn't enough data.
@@ -1247,7 +1286,7 @@ int ssh2_bignum_length(Bignum bn)
  */
 int bignum_byte(Bignum bn, int i)
 {
-    if (i >= (int)(BIGNUM_INT_BYTES * bn[0]))
+    if (i < 0 || i >= (int)(BIGNUM_INT_BYTES * bn[0]))
        return 0;                      /* beyond the end */
     else
        return (bn[i / BIGNUM_INT_BYTES + 1] >>
@@ -1259,7 +1298,7 @@ int bignum_byte(Bignum bn, int i)
  */
 int bignum_bit(Bignum bn, int i)
 {
-    if (i >= (int)(BIGNUM_INT_BITS * bn[0]))
+    if (i < 0 || i >= (int)(BIGNUM_INT_BITS * bn[0]))
        return 0;                      /* beyond the end */
     else
        return (bn[i / BIGNUM_INT_BITS + 1] >> (i % BIGNUM_INT_BITS)) & 1;
@@ -1270,11 +1309,11 @@ int bignum_bit(Bignum bn, int i)
  */
 void bignum_set_bit(Bignum bn, int bitnum, int value)
 {
-    if (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
@@ -1306,7 +1345,18 @@ int ssh1_write_bignum(void *data, Bignum bn)
 int bignum_cmp(Bignum a, Bignum b)
 {
     int amax = a[0], bmax = b[0];
-    int i = (amax > bmax ? amax : bmax);
+    int i;
+
+    /* Annoyingly we have two representations of zero */
+    if (amax == 1 && a[amax] == 0)
+        amax = 0;
+    if (bmax == 1 && b[bmax] == 0)
+        bmax = 0;
+
+    assert(amax == 0 || a[amax] != 0);
+    assert(bmax == 0 || b[bmax] != 0);
+
+    i = (amax > bmax ? amax : bmax);
     while (i) {
        BignumInt aval = (i > amax ? 0 : a[i]);
        BignumInt bval = (i > bmax ? 0 : b[i]);
@@ -1328,6 +1378,8 @@ Bignum bignum_rshift(Bignum a, int shift)
     int i, shiftw, shiftb, shiftbb, bits;
     BignumInt ai, ai1;
 
+    assert(shift >= 0);
+
     bits = bignum_bitcount(a) - shift;
     ret = newbn((bits + BIGNUM_INT_BITS - 1) / BIGNUM_INT_BITS);
 
@@ -1347,6 +1399,44 @@ Bignum bignum_rshift(Bignum a, int shift)
     return ret;
 }
 
+/*
+ * Left-shift one bignum to form another.
+ */
+Bignum bignum_lshift(Bignum a, int shift)
+{
+    Bignum ret;
+    int bits, shiftWords, shiftBits;
+
+    assert(shift >= 0);
+
+    bits = bignum_bitcount(a) + shift;
+    ret = newbn((bits + BIGNUM_INT_BITS - 1) / BIGNUM_INT_BITS);
+
+    shiftWords = shift / BIGNUM_INT_BITS;
+    shiftBits = shift % BIGNUM_INT_BITS;
+
+    if (shiftBits == 0)
+    {
+        memcpy(&ret[1 + shiftWords], &a[1], sizeof(BignumInt) * a[0]);
+    }
+    else
+    {
+        int i;
+        BignumInt carry = 0;
+
+        /* Remember that Bignum[0] is length, so add 1 */
+        for (i = shiftWords + 1; i < ((int)a[0]) + shiftWords + 1; ++i)
+        {
+            BignumInt from = a[i - shiftWords];
+            ret[i] = (from << shiftBits) | carry;
+            carry = from >> (BIGNUM_INT_BITS - shiftBits);
+        }
+        if (carry) ret[i] = carry;
+    }
+
+    return ret;
+}
+
 /*
  * Non-modular multiplication and addition.
  */
@@ -1398,8 +1488,7 @@ Bignum bigmuladd(Bignum a, Bignum b, Bignum addend)
     }
     ret[0] = maxspot;
 
-    for (i = 0; i < wslen; i++)
-        workspace[i] = 0;
+    smemclr(workspace, wslen * sizeof(*workspace));
     sfree(workspace);
     return ret;
 }
@@ -1583,6 +1672,8 @@ Bignum bigdiv(Bignum a, Bignum b)
 {
     Bignum q = newbn(a[0]);
     bigdivmod(a, b, NULL, q);
+    while (q[0] > 1 && q[q[0]] == 0)
+        q[0]--;
     return q;
 }
 
@@ -1593,6 +1684,8 @@ Bignum bigmod(Bignum a, Bignum b)
 {
     Bignum r = newbn(b[0]);
     bigdivmod(a, b, r, NULL);
+    while (r[0] > 1 && r[r[0]] == 0)
+        r[0]--;
     return r;
 }
 
@@ -1629,12 +1722,31 @@ Bignum modinv(Bignum number, Bignum modulus)
     Bignum x = copybn(One);
     int sign = +1;
 
+    assert(number[number[0]] != 0);
+    assert(modulus[modulus[0]] != 0);
+
     while (bignum_cmp(b, One) != 0) {
-       Bignum t = newbn(b[0]);
-       Bignum q = newbn(a[0]);
+       Bignum t, q;
+
+        if (bignum_cmp(b, Zero) == 0) {
+            /*
+             * Found a common factor between the inputs, so we cannot
+             * return a modular inverse at all.
+             */
+            freebn(b);
+            freebn(a);
+            freebn(xp);
+            freebn(x);
+            return NULL;
+        }
+
+        t = newbn(b[0]);
+       q = newbn(a[0]);
        bigdivmod(a, b, t, q);
        while (t[0] > 1 && t[t[0]] == 0)
            t[0]--;
+       while (q[0] > 1 && q[q[0]] == 0)
+           q[0]--;
        freebn(a);
        a = b;
        b = t;
@@ -1750,6 +1862,7 @@ char *bignum_decimal(Bignum x)
     /*
      * Done.
      */
+    smemclr(workspace, x[0] * sizeof(*workspace));
     sfree(workspace);
     return ret;
 }
@@ -1761,13 +1874,13 @@ char *bignum_decimal(Bignum x)
 #include <ctype.h>
 
 /*
- * gcc -g -O0 -DTESTBN -o testbn sshbn.c misc.c -I unix -I charset
+ * 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(char *p, ...)
+void modalfatalbox(const char *p, ...)
 {
     va_list ap;
     fprintf(stderr, "FATAL ERROR: ");
@@ -1778,6 +1891,12 @@ void modalfatalbox(char *p, ...)
     exit(1);
 }
 
+int random_byte(void)
+{
+    modalfatalbox("random_byte called in testbn");
+    return 0;
+}
+
 #define fromxdigit(c) ( (c)>'9' ? ((c)&0xDF) - 'A' + 10 : (c) - '0' )
 
 int main(int argc, char **argv)
@@ -1864,11 +1983,49 @@ int main(int argc, char **argv)
             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);
+                printf("%d: pow with %d parameters, expected 4\n", line, ptrnum);
                 exit(1);
             }