sshbn.c

   1 /*
   2  * Bignum routines for RSA and DH and stuff.
   3  */
   4
   5 #include <stdio.h>
   6 #include <assert.h>
   7 #include <stdlib.h>
   8 #include <string.h>
   9
  10 #include "misc.h"
  11
  12 #if defined __GNUC__ && defined __i386__
  13 typedef unsigned long BignumInt;
  14 typedef unsigned long long BignumDblInt;
  15 #define BIGNUM_INT_MASK  0xFFFFFFFFUL
  16 #define BIGNUM_TOP_BIT   0x80000000UL
  17 #define BIGNUM_INT_BITS  32
  18 #define MUL_WORD(w1, w2) ((BignumDblInt)w1 * w2)
  19 #define DIVMOD_WORD(q, r, hi, lo, w) \
  20     __asm__("div %2" : \
  21             "=d" (r), "=a" (q) : \
  22             "r" (w), "d" (hi), "a" (lo))
  23 #else
  24 typedef unsigned short BignumInt;
  25 typedef unsigned long BignumDblInt;
  26 #define BIGNUM_INT_MASK  0xFFFFU
  27 #define BIGNUM_TOP_BIT   0x8000U
  28 #define BIGNUM_INT_BITS  16
  29 #define MUL_WORD(w1, w2) ((BignumDblInt)w1 * w2)
  30 #define DIVMOD_WORD(q, r, hi, lo, w) do { \
  31     BignumDblInt n = (((BignumDblInt)hi) << BIGNUM_INT_BITS) | lo; \
  32     q = n / w; \
  33     r = n % w; \
  34 } while (0)
  35 #endif
  36
  37 #define BIGNUM_INT_BYTES (BIGNUM_INT_BITS / 8)
  38
  39 #define BIGNUM_INTERNAL
  40 typedef BignumInt *Bignum;
  41
  42 #include "ssh.h"
  43
  44 BignumInt bnZero[1] = { 0 };
  45 BignumInt bnOne[2] = { 1, 1 };
  46
  47 /*
  48  * The Bignum format is an array of `BignumInt'. The first
  49  * element of the array counts the remaining elements. The
  50  * remaining elements express the actual number, base 2^BIGNUM_INT_BITS, _least_
  51  * significant digit first. (So it's trivial to extract the bit
  52  * with value 2^n for any n.)
  53  *
  54  * All Bignums in this module are positive. Negative numbers must
  55  * be dealt with outside it.
  56  *
  57  * INVARIANT: the most significant word of any Bignum must be
  58  * nonzero.
  59  */
  60
  61 Bignum Zero = bnZero, One = bnOne;
  62
  63 static Bignum newbn(int length)
  64 {
  65     Bignum b = snewn(length + 1, BignumInt);
  66     if (!b)
  67         abort();                       /* FIXME */
  68     memset(b, 0, (length + 1) * sizeof(*b));
  69     b[0] = length;
  70     return b;
  71 }
  72
  73 void bn_restore_invariant(Bignum b)
  74 {
  75     while (b[0] > 1 && b[b[0]] == 0)
  76         b[0]--;
  77 }
  78
  79 Bignum copybn(Bignum orig)
  80 {
  81     Bignum b = snewn(orig[0] + 1, BignumInt);
  82     if (!b)
  83         abort();                       /* FIXME */
  84     memcpy(b, orig, (orig[0] + 1) * sizeof(*b));
  85     return b;
  86 }
  87
  88 void freebn(Bignum b)
  89 {
  90     /*
  91      * Burn the evidence, just in case.
  92      */
  93     memset(b, 0, sizeof(b[0]) * (b[0] + 1));
  94     sfree(b);
  95 }
  96
  97 Bignum bn_power_2(int n)
  98 {
  99     Bignum ret = newbn(n / BIGNUM_INT_BITS + 1);
 100     bignum_set_bit(ret, n, 1);
 101     return ret;
 102 }
 103
 104 /*
 105  * Compute c = a * b.
 106  * Input is in the first len words of a and b.
 107  * Result is returned in the first 2*len words of c.
 108  */
 109 static void internal_mul(BignumInt *a, BignumInt *b,
 110                          BignumInt *c, int len)
 111 {
 112     int i, j;
 113     BignumDblInt t;
 114
 115     for (j = 0; j < 2 * len; j++)
 116         c[j] = 0;
 117
 118     for (i = len - 1; i >= 0; i--) {
 119         t = 0;
 120         for (j = len - 1; j >= 0; j--) {
 121             t += MUL_WORD(a[i], (BignumDblInt) b[j]);
 122             t += (BignumDblInt) c[i + j + 1];
 123             c[i + j + 1] = (BignumInt) t;
 124             t = t >> BIGNUM_INT_BITS;
 125         }
 126         c[i] = (BignumInt) t;
 127     }
 128 }
 129
 130 static void internal_add_shifted(BignumInt *number,
 131                                  unsigned n, int shift)
 132 {
 133     int word = 1 + (shift / BIGNUM_INT_BITS);
 134     int bshift = shift % BIGNUM_INT_BITS;
 135     BignumDblInt addend;
 136
 137     addend = (BignumDblInt)n << bshift;
 138
 139     while (addend) {
 140         addend += number[word];
 141         number[word] = (BignumInt) addend & BIGNUM_INT_MASK;
 142         addend >>= BIGNUM_INT_BITS;
 143         word++;
 144     }
 145 }
 146
 147 /*
 148  * Compute a = a % m.
 149  * Input in first alen words of a and first mlen words of m.
 150  * Output in first alen words of a
 151  * (of which first alen-mlen words will be zero).
 152  * The MSW of m MUST have its high bit set.
 153  * Quotient is accumulated in the `quotient' array, which is a Bignum
 154  * rather than the internal bigendian format. Quotient parts are shifted
 155  * left by `qshift' before adding into quot.
 156  */
 157 static void internal_mod(BignumInt *a, int alen,
 158                          BignumInt *m, int mlen,
 159                          BignumInt *quot, int qshift)
 160 {
 161     BignumInt m0, m1;
 162     unsigned int h;
 163     int i, k;
 164
 165     m0 = m[0];
 166     if (mlen > 1)
 167         m1 = m[1];
 168     else
 169         m1 = 0;
 170
 171     for (i = 0; i <= alen - mlen; i++) {
 172         BignumDblInt t;
 173         unsigned int q, r, c, ai1;
 174
 175         if (i == 0) {
 176             h = 0;
 177         } else {
 178             h = a[i - 1];
 179             a[i - 1] = 0;
 180         }
 181
 182         if (i == alen - 1)
 183             ai1 = 0;
 184         else
 185             ai1 = a[i + 1];
 186
 187         /* Find q = h:a[i] / m0 */
 188         if (h >= m0) {
 189             /*
 190              * Special case.
 191              *
 192              * To illustrate it, suppose a BignumInt is 8 bits, and
 193              * we are dividing (say) A1:23:45:67 by A1:B2:C3. Then
 194              * our initial division will be 0xA123 / 0xA1, which
 195              * will give a quotient of 0x100 and a divide overflow.
 196              * However, the invariants in this division algorithm
 197              * are not violated, since the full number A1:23:... is
 198              * _less_ than the quotient prefix A1:B2:... and so the
 199              * following correction loop would have sorted it out.
 200              *
 201              * In this situation we set q to be the largest
 202              * quotient we _can_ stomach (0xFF, of course).
 203              */
 204             q = BIGNUM_INT_MASK;
 205         } else {
 206             DIVMOD_WORD(q, r, h, a[i], m0);
 207
 208             /* Refine our estimate of q by looking at
 209              h:a[i]:a[i+1] / m0:m1 */
 210             t = MUL_WORD(m1, q);
 211             if (t > ((BignumDblInt) r << BIGNUM_INT_BITS) + ai1) {
 212                 q--;
 213                 t -= m1;
 214                 r = (r + m0) & BIGNUM_INT_MASK;     /* overflow? */
 215                 if (r >= (BignumDblInt) m0 &&
 216                     t > ((BignumDblInt) r << BIGNUM_INT_BITS) + ai1) q--;
 217             }
 218         }
 219
 220         /* Subtract q * m from a[i...] */
 221         c = 0;
 222         for (k = mlen - 1; k >= 0; k--) {
 223             t = MUL_WORD(q, m[k]);
 224             t += c;
 225             c = t >> BIGNUM_INT_BITS;
 226             if ((BignumInt) t > a[i + k])
 227                 c++;
 228             a[i + k] -= (BignumInt) t;
 229         }
 230
 231         /* Add back m in case of borrow */
 232         if (c != h) {
 233             t = 0;
 234             for (k = mlen - 1; k >= 0; k--) {
 235                 t += m[k];
 236                 t += a[i + k];
 237                 a[i + k] = (BignumInt) t;
 238                 t = t >> BIGNUM_INT_BITS;
 239             }
 240             q--;
 241         }
 242         if (quot)
 243             internal_add_shifted(quot, q, qshift + BIGNUM_INT_BITS * (alen - mlen - i));
 244     }
 245 }
 246
 247 /*
 248  * Compute (base ^ exp) % mod.
 249  */
 250 Bignum modpow(Bignum base_in, Bignum exp, Bignum mod)
 251 {
 252     BignumInt *a, *b, *n, *m;
 253     int mshift;
 254     int mlen, i, j;
 255     Bignum base, result;
 256
 257     /*
 258      * The most significant word of mod needs to be non-zero. It
 259      * should already be, but let's make sure.
 260      */
 261     assert(mod[mod[0]] != 0);
 262
 263     /*
 264      * Make sure the base is smaller than the modulus, by reducing
 265      * it modulo the modulus if not.
 266      */
 267     base = bigmod(base_in, mod);
 268
 269     /* Allocate m of size mlen, copy mod to m */
 270     /* We use big endian internally */
 271     mlen = mod[0];
 272     m = snewn(mlen, BignumInt);
 273     for (j = 0; j < mlen; j++)
 274         m[j] = mod[mod[0] - j];
 275
 276     /* Shift m left to make msb bit set */
 277     for (mshift = 0; mshift < BIGNUM_INT_BITS-1; mshift++)
 278         if ((m[0] << mshift) & BIGNUM_TOP_BIT)
 279             break;
 280     if (mshift) {
 281         for (i = 0; i < mlen - 1; i++)
 282             m[i] = (m[i] << mshift) | (m[i + 1] >> (BIGNUM_INT_BITS - mshift));
 283         m[mlen - 1] = m[mlen - 1] << mshift;
 284     }
 285
 286     /* Allocate n of size mlen, copy base to n */
 287     n = snewn(mlen, BignumInt);
 288     i = mlen - base[0];
 289     for (j = 0; j < i; j++)
 290         n[j] = 0;
 291     for (j = 0; j < base[0]; j++)
 292         n[i + j] = base[base[0] - j];
 293
 294     /* Allocate a and b of size 2*mlen. Set a = 1 */
 295     a = snewn(2 * mlen, BignumInt);
 296     b = snewn(2 * mlen, BignumInt);
 297     for (i = 0; i < 2 * mlen; i++)
 298         a[i] = 0;
 299     a[2 * mlen - 1] = 1;
 300
 301     /* Skip leading zero bits of exp. */
 302     i = 0;
 303     j = BIGNUM_INT_BITS-1;
 304     while (i < exp[0] && (exp[exp[0] - i] & (1 << j)) == 0) {
 305         j--;
 306         if (j < 0) {
 307             i++;
 308             j = BIGNUM_INT_BITS-1;
 309         }
 310     }
 311
 312     /* Main computation */
 313     while (i < exp[0]) {
 314         while (j >= 0) {
 315             internal_mul(a + mlen, a + mlen, b, mlen);
 316             internal_mod(b, mlen * 2, m, mlen, NULL, 0);
 317             if ((exp[exp[0] - i] & (1 << j)) != 0) {
 318                 internal_mul(b + mlen, n, a, mlen);
 319                 internal_mod(a, mlen * 2, m, mlen, NULL, 0);
 320             } else {
 321                 BignumInt *t;
 322                 t = a;
 323                 a = b;
 324                 b = t;
 325             }
 326             j--;
 327         }
 328         i++;
 329         j = BIGNUM_INT_BITS-1;
 330     }
 331
 332     /* Fixup result in case the modulus was shifted */
 333     if (mshift) {
 334         for (i = mlen - 1; i < 2 * mlen - 1; i++)
 335             a[i] = (a[i] << mshift) | (a[i + 1] >> (BIGNUM_INT_BITS - mshift));
 336         a[2 * mlen - 1] = a[2 * mlen - 1] << mshift;
 337         internal_mod(a, mlen * 2, m, mlen, NULL, 0);
 338         for (i = 2 * mlen - 1; i >= mlen; i--)
 339             a[i] = (a[i] >> mshift) | (a[i - 1] << (BIGNUM_INT_BITS - mshift));
 340     }
 341
 342     /* Copy result to buffer */
 343     result = newbn(mod[0]);
 344     for (i = 0; i < mlen; i++)
 345         result[result[0] - i] = a[i + mlen];
 346     while (result[0] > 1 && result[result[0]] == 0)
 347         result[0]--;
 348
 349     /* Free temporary arrays */
 350     for (i = 0; i < 2 * mlen; i++)
 351         a[i] = 0;
 352     sfree(a);
 353     for (i = 0; i < 2 * mlen; i++)
 354         b[i] = 0;
 355     sfree(b);
 356     for (i = 0; i < mlen; i++)
 357         m[i] = 0;
 358     sfree(m);
 359     for (i = 0; i < mlen; i++)
 360         n[i] = 0;
 361     sfree(n);
 362
 363     freebn(base);
 364
 365     return result;
 366 }
 367
 368 /*
 369  * Compute (p * q) % mod.
 370  * The most significant word of mod MUST be non-zero.
 371  * We assume that the result array is the same size as the mod array.
 372  */
 373 Bignum modmul(Bignum p, Bignum q, Bignum mod)
 374 {
 375     BignumInt *a, *n, *m, *o;
 376     int mshift;
 377     int pqlen, mlen, rlen, i, j;
 378     Bignum result;
 379
 380     /* Allocate m of size mlen, copy mod to m */
 381     /* We use big endian internally */
 382     mlen = mod[0];
 383     m = snewn(mlen, BignumInt);
 384     for (j = 0; j < mlen; j++)
 385         m[j] = mod[mod[0] - j];
 386
 387     /* Shift m left to make msb bit set */
 388     for (mshift = 0; mshift < BIGNUM_INT_BITS-1; mshift++)
 389         if ((m[0] << mshift) & BIGNUM_TOP_BIT)
 390             break;
 391     if (mshift) {
 392         for (i = 0; i < mlen - 1; i++)
 393             m[i] = (m[i] << mshift) | (m[i + 1] >> (BIGNUM_INT_BITS - mshift));
 394         m[mlen - 1] = m[mlen - 1] << mshift;
 395     }
 396
 397     pqlen = (p[0] > q[0] ? p[0] : q[0]);
 398
 399     /* Allocate n of size pqlen, copy p to n */
 400     n = snewn(pqlen, BignumInt);
 401     i = pqlen - p[0];
 402     for (j = 0; j < i; j++)
 403         n[j] = 0;
 404     for (j = 0; j < p[0]; j++)
 405         n[i + j] = p[p[0] - j];
 406
 407     /* Allocate o of size pqlen, copy q to o */
 408     o = snewn(pqlen, BignumInt);
 409     i = pqlen - q[0];
 410     for (j = 0; j < i; j++)
 411         o[j] = 0;
 412     for (j = 0; j < q[0]; j++)
 413         o[i + j] = q[q[0] - j];
 414
 415     /* Allocate a of size 2*pqlen for result */
 416     a = snewn(2 * pqlen, BignumInt);
 417
 418     /* Main computation */
 419     internal_mul(n, o, a, pqlen);
 420     internal_mod(a, pqlen * 2, m, mlen, NULL, 0);
 421
 422     /* Fixup result in case the modulus was shifted */
 423     if (mshift) {
 424         for (i = 2 * pqlen - mlen - 1; i < 2 * pqlen - 1; i++)
 425             a[i] = (a[i] << mshift) | (a[i + 1] >> (BIGNUM_INT_BITS - mshift));
 426         a[2 * pqlen - 1] = a[2 * pqlen - 1] << mshift;
 427         internal_mod(a, pqlen * 2, m, mlen, NULL, 0);
 428         for (i = 2 * pqlen - 1; i >= 2 * pqlen - mlen; i--)
 429             a[i] = (a[i] >> mshift) | (a[i - 1] << (BIGNUM_INT_BITS - mshift));
 430     }
 431
 432     /* Copy result to buffer */
 433     rlen = (mlen < pqlen * 2 ? mlen : pqlen * 2);
 434     result = newbn(rlen);
 435     for (i = 0; i < rlen; i++)
 436         result[result[0] - i] = a[i + 2 * pqlen - rlen];
 437     while (result[0] > 1 && result[result[0]] == 0)
 438         result[0]--;
 439
 440     /* Free temporary arrays */
 441     for (i = 0; i < 2 * pqlen; i++)
 442         a[i] = 0;
 443     sfree(a);
 444     for (i = 0; i < mlen; i++)
 445         m[i] = 0;
 446     sfree(m);
 447     for (i = 0; i < pqlen; i++)
 448         n[i] = 0;
 449     sfree(n);
 450     for (i = 0; i < pqlen; i++)
 451         o[i] = 0;
 452     sfree(o);
 453
 454     return result;
 455 }
 456
 457 /*
 458  * Compute p % mod.
 459  * The most significant word of mod MUST be non-zero.
 460  * We assume that the result array is the same size as the mod array.
 461  * We optionally write out a quotient if `quotient' is non-NULL.
 462  * We can avoid writing out the result if `result' is NULL.
 463  */
 464 static void bigdivmod(Bignum p, Bignum mod, Bignum result, Bignum quotient)
 465 {
 466     BignumInt *n, *m;
 467     int mshift;
 468     int plen, mlen, i, j;
 469
 470     /* Allocate m of size mlen, copy mod to m */
 471     /* We use big endian internally */
 472     mlen = mod[0];
 473     m = snewn(mlen, BignumInt);
 474     for (j = 0; j < mlen; j++)
 475         m[j] = mod[mod[0] - j];
 476
 477     /* Shift m left to make msb bit set */
 478     for (mshift = 0; mshift < BIGNUM_INT_BITS-1; mshift++)
 479         if ((m[0] << mshift) & BIGNUM_TOP_BIT)
 480             break;
 481     if (mshift) {
 482         for (i = 0; i < mlen - 1; i++)
 483             m[i] = (m[i] << mshift) | (m[i + 1] >> (BIGNUM_INT_BITS - mshift));
 484         m[mlen - 1] = m[mlen - 1] << mshift;
 485     }
 486
 487     plen = p[0];
 488     /* Ensure plen > mlen */
 489     if (plen <= mlen)
 490         plen = mlen + 1;
 491
 492     /* Allocate n of size plen, copy p to n */
 493     n = snewn(plen, BignumInt);
 494     for (j = 0; j < plen; j++)
 495         n[j] = 0;
 496     for (j = 1; j <= p[0]; j++)
 497         n[plen - j] = p[j];
 498
 499     /* Main computation */
 500     internal_mod(n, plen, m, mlen, quotient, mshift);
 501
 502     /* Fixup result in case the modulus was shifted */
 503     if (mshift) {
 504         for (i = plen - mlen - 1; i < plen - 1; i++)
 505             n[i] = (n[i] << mshift) | (n[i + 1] >> (BIGNUM_INT_BITS - mshift));
 506         n[plen - 1] = n[plen - 1] << mshift;
 507         internal_mod(n, plen, m, mlen, quotient, 0);
 508         for (i = plen - 1; i >= plen - mlen; i--)
 509             n[i] = (n[i] >> mshift) | (n[i - 1] << (BIGNUM_INT_BITS - mshift));
 510     }
 511
 512     /* Copy result to buffer */
 513     if (result) {
 514         for (i = 1; i <= result[0]; i++) {
 515             int j = plen - i;
 516             result[i] = j >= 0 ? n[j] : 0;
 517         }
 518     }
 519
 520     /* Free temporary arrays */
 521     for (i = 0; i < mlen; i++)
 522         m[i] = 0;
 523     sfree(m);
 524     for (i = 0; i < plen; i++)
 525         n[i] = 0;
 526     sfree(n);
 527 }
 528
 529 /*
 530  * Decrement a number.
 531  */
 532 void decbn(Bignum bn)
 533 {
 534     int i = 1;
 535     while (i < bn[0] && bn[i] == 0)
 536         bn[i++] = BIGNUM_INT_MASK;
 537     bn[i]--;
 538 }
 539
 540 Bignum bignum_from_bytes(const unsigned char *data, int nbytes)
 541 {
 542     Bignum result;
 543     int w, i;
 544
 545     w = (nbytes + BIGNUM_INT_BYTES - 1) / BIGNUM_INT_BYTES; /* bytes->words */
 546
 547     result = newbn(w);
 548     for (i = 1; i <= w; i++)
 549         result[i] = 0;
 550     for (i = nbytes; i--;) {
 551         unsigned char byte = *data++;
 552         result[1 + i / BIGNUM_INT_BYTES] |= byte << (8*i % BIGNUM_INT_BITS);
 553     }
 554
 555     while (result[0] > 1 && result[result[0]] == 0)
 556         result[0]--;
 557     return result;
 558 }
 559
 560 /*
 561  * Read an ssh1-format bignum from a data buffer. Return the number
 562  * of bytes consumed, or -1 if there wasn't enough data.
 563  */
 564 int ssh1_read_bignum(const unsigned char *data, int len, Bignum * result)
 565 {
 566     const unsigned char *p = data;
 567     int i;
 568     int w, b;
 569
 570     if (len < 2)
 571         return -1;
 572
 573     w = 0;
 574     for (i = 0; i < 2; i++)
 575         w = (w << 8) + *p++;
 576     b = (w + 7) / 8;                   /* bits -> bytes */
 577
 578     if (len < b+2)
 579         return -1;
 580
 581     if (!result)                       /* just return length */
 582         return b + 2;
 583
 584     *result = bignum_from_bytes(p, b);
 585
 586     return p + b - data;
 587 }
 588
 589 /*
 590  * Return the bit count of a bignum, for ssh1 encoding.
 591  */
 592 int bignum_bitcount(Bignum bn)
 593 {
 594     int bitcount = bn[0] * BIGNUM_INT_BITS - 1;
 595     while (bitcount >= 0
 596            && (bn[bitcount / BIGNUM_INT_BITS + 1] >> (bitcount % BIGNUM_INT_BITS)) == 0) bitcount--;
 597     return bitcount + 1;
 598 }
 599
 600 /*
 601  * Return the byte length of a bignum when ssh1 encoded.
 602  */
 603 int ssh1_bignum_length(Bignum bn)
 604 {
 605     return 2 + (bignum_bitcount(bn) + 7) / 8;
 606 }
 607
 608 /*
 609  * Return the byte length of a bignum when ssh2 encoded.
 610  */
 611 int ssh2_bignum_length(Bignum bn)
 612 {
 613     return 4 + (bignum_bitcount(bn) + 8) / 8;
 614 }
 615
 616 /*
 617  * Return a byte from a bignum; 0 is least significant, etc.
 618  */
 619 int bignum_byte(Bignum bn, int i)
 620 {
 621     if (i >= BIGNUM_INT_BYTES * bn[0])
 622         return 0;                      /* beyond the end */
 623     else
 624         return (bn[i / BIGNUM_INT_BYTES + 1] >>
 625                 ((i % BIGNUM_INT_BYTES)*8)) & 0xFF;
 626 }
 627
 628 /*
 629  * Return a bit from a bignum; 0 is least significant, etc.
 630  */
 631 int bignum_bit(Bignum bn, int i)
 632 {
 633     if (i >= BIGNUM_INT_BITS * bn[0])
 634         return 0;                      /* beyond the end */
 635     else
 636         return (bn[i / BIGNUM_INT_BITS + 1] >> (i % BIGNUM_INT_BITS)) & 1;
 637 }
 638
 639 /*
 640  * Set a bit in a bignum; 0 is least significant, etc.
 641  */
 642 void bignum_set_bit(Bignum bn, int bitnum, int value)
 643 {
 644     if (bitnum >= BIGNUM_INT_BITS * bn[0])
 645         abort();                       /* beyond the end */
 646     else {
 647         int v = bitnum / BIGNUM_INT_BITS + 1;
 648         int mask = 1 << (bitnum % BIGNUM_INT_BITS);
 649         if (value)
 650             bn[v] |= mask;
 651         else
 652             bn[v] &= ~mask;
 653     }
 654 }
 655
 656 /*
 657  * Write a ssh1-format bignum into a buffer. It is assumed the
 658  * buffer is big enough. Returns the number of bytes used.
 659  */
 660 int ssh1_write_bignum(void *data, Bignum bn)
 661 {
 662     unsigned char *p = data;
 663     int len = ssh1_bignum_length(bn);
 664     int i;
 665     int bitc = bignum_bitcount(bn);
 666
 667     *p++ = (bitc >> 8) & 0xFF;
 668     *p++ = (bitc) & 0xFF;
 669     for (i = len - 2; i--;)
 670         *p++ = bignum_byte(bn, i);
 671     return len;
 672 }
 673
 674 /*
 675  * Compare two bignums. Returns like strcmp.
 676  */
 677 int bignum_cmp(Bignum a, Bignum b)
 678 {
 679     int amax = a[0], bmax = b[0];
 680     int i = (amax > bmax ? amax : bmax);
 681     while (i) {
 682         BignumInt aval = (i > amax ? 0 : a[i]);
 683         BignumInt bval = (i > bmax ? 0 : b[i]);
 684         if (aval < bval)
 685             return -1;
 686         if (aval > bval)
 687             return +1;
 688         i--;
 689     }
 690     return 0;
 691 }
 692
 693 /*
 694  * Right-shift one bignum to form another.
 695  */
 696 Bignum bignum_rshift(Bignum a, int shift)
 697 {
 698     Bignum ret;
 699     int i, shiftw, shiftb, shiftbb, bits;
 700     BignumInt ai, ai1;
 701
 702     bits = bignum_bitcount(a) - shift;
 703     ret = newbn((bits + BIGNUM_INT_BITS - 1) / BIGNUM_INT_BITS);
 704
 705     if (ret) {
 706         shiftw = shift / BIGNUM_INT_BITS;
 707         shiftb = shift % BIGNUM_INT_BITS;
 708         shiftbb = BIGNUM_INT_BITS - shiftb;
 709
 710         ai1 = a[shiftw + 1];
 711         for (i = 1; i <= ret[0]; i++) {
 712             ai = ai1;
 713             ai1 = (i + shiftw + 1 <= a[0] ? a[i + shiftw + 1] : 0);
 714             ret[i] = ((ai >> shiftb) | (ai1 << shiftbb)) & BIGNUM_INT_MASK;
 715         }
 716     }
 717
 718     return ret;
 719 }
 720
 721 /*
 722  * Non-modular multiplication and addition.
 723  */
 724 Bignum bigmuladd(Bignum a, Bignum b, Bignum addend)
 725 {
 726     int alen = a[0], blen = b[0];
 727     int mlen = (alen > blen ? alen : blen);
 728     int rlen, i, maxspot;
 729     BignumInt *workspace;
 730     Bignum ret;
 731
 732     /* mlen space for a, mlen space for b, 2*mlen for result */
 733     workspace = snewn(mlen * 4, BignumInt);
 734     for (i = 0; i < mlen; i++) {
 735         workspace[0 * mlen + i] = (mlen - i <= a[0] ? a[mlen - i] : 0);
 736         workspace[1 * mlen + i] = (mlen - i <= b[0] ? b[mlen - i] : 0);
 737     }
 738
 739     internal_mul(workspace + 0 * mlen, workspace + 1 * mlen,
 740                  workspace + 2 * mlen, mlen);
 741
 742     /* now just copy the result back */
 743     rlen = alen + blen + 1;
 744     if (addend && rlen <= addend[0])
 745         rlen = addend[0] + 1;
 746     ret = newbn(rlen);
 747     maxspot = 0;
 748     for (i = 1; i <= ret[0]; i++) {
 749         ret[i] = (i <= 2 * mlen ? workspace[4 * mlen - i] : 0);
 750         if (ret[i] != 0)
 751             maxspot = i;
 752     }
 753     ret[0] = maxspot;
 754
 755     /* now add in the addend, if any */
 756     if (addend) {
 757         BignumDblInt carry = 0;
 758         for (i = 1; i <= rlen; i++) {
 759             carry += (i <= ret[0] ? ret[i] : 0);
 760             carry += (i <= addend[0] ? addend[i] : 0);
 761             ret[i] = (BignumInt) carry & BIGNUM_INT_MASK;
 762             carry >>= BIGNUM_INT_BITS;
 763             if (ret[i] != 0 && i > maxspot)
 764                 maxspot = i;
 765         }
 766     }
 767     ret[0] = maxspot;
 768
 769     sfree(workspace);
 770     return ret;
 771 }
 772
 773 /*
 774  * Non-modular multiplication.
 775  */
 776 Bignum bigmul(Bignum a, Bignum b)
 777 {
 778     return bigmuladd(a, b, NULL);
 779 }
 780
 781 /*
 782  * Create a bignum which is the bitmask covering another one. That
 783  * is, the smallest integer which is >= N and is also one less than
 784  * a power of two.
 785  */
 786 Bignum bignum_bitmask(Bignum n)
 787 {
 788     Bignum ret = copybn(n);
 789     int i;
 790     BignumInt j;
 791
 792     i = ret[0];
 793     while (n[i] == 0 && i > 0)
 794         i--;
 795     if (i <= 0)
 796         return ret;                    /* input was zero */
 797     j = 1;
 798     while (j < n[i])
 799         j = 2 * j + 1;
 800     ret[i] = j;
 801     while (--i > 0)
 802         ret[i] = BIGNUM_INT_MASK;
 803     return ret;
 804 }
 805
 806 /*
 807  * Convert a (max 32-bit) long into a bignum.
 808  */
 809 Bignum bignum_from_long(unsigned long nn)
 810 {
 811     Bignum ret;
 812     BignumDblInt n = nn;
 813
 814     ret = newbn(3);
 815     ret[1] = (BignumInt)(n & BIGNUM_INT_MASK);
 816     ret[2] = (BignumInt)((n >> BIGNUM_INT_BITS) & BIGNUM_INT_MASK);
 817     ret[3] = 0;
 818     ret[0] = (ret[2]  ? 2 : 1);
 819     return ret;
 820 }
 821
 822 /*
 823  * Add a long to a bignum.
 824  */
 825 Bignum bignum_add_long(Bignum number, unsigned long addendx)
 826 {
 827     Bignum ret = newbn(number[0] + 1);
 828     int i, maxspot = 0;
 829     BignumDblInt carry = 0, addend = addendx;
 830
 831     for (i = 1; i <= ret[0]; i++) {
 832         carry += addend & BIGNUM_INT_MASK;
 833         carry += (i <= number[0] ? number[i] : 0);
 834         addend >>= BIGNUM_INT_BITS;
 835         ret[i] = (BignumInt) carry & BIGNUM_INT_MASK;
 836         carry >>= BIGNUM_INT_BITS;
 837         if (ret[i] != 0)
 838             maxspot = i;
 839     }
 840     ret[0] = maxspot;
 841     return ret;
 842 }
 843
 844 /*
 845  * Compute the residue of a bignum, modulo a (max 16-bit) short.
 846  */
 847 unsigned short bignum_mod_short(Bignum number, unsigned short modulus)
 848 {
 849     BignumDblInt mod, r;
 850     int i;
 851
 852     r = 0;
 853     mod = modulus;
 854     for (i = number[0]; i > 0; i--)
 855         r = (r * (BIGNUM_TOP_BIT % mod) * 2 + number[i] % mod) % mod;
 856     return (unsigned short) r;
 857 }
 858
 859 #ifdef DEBUG
 860 void diagbn(char *prefix, Bignum md)
 861 {
 862     int i, nibbles, morenibbles;
 863     static const char hex[] = "0123456789ABCDEF";
 864
 865     debug(("%s0x", prefix ? prefix : ""));
 866
 867     nibbles = (3 + bignum_bitcount(md)) / 4;
 868     if (nibbles < 1)
 869         nibbles = 1;
 870     morenibbles = 4 * md[0] - nibbles;
 871     for (i = 0; i < morenibbles; i++)
 872         debug(("-"));
 873     for (i = nibbles; i--;)
 874         debug(("%c",
 875                hex[(bignum_byte(md, i / 2) >> (4 * (i % 2))) & 0xF]));
 876
 877     if (prefix)
 878         debug(("\n"));
 879 }
 880 #endif
 881
 882 /*
 883  * Simple division.
 884  */
 885 Bignum bigdiv(Bignum a, Bignum b)
 886 {
 887     Bignum q = newbn(a[0]);
 888     bigdivmod(a, b, NULL, q);
 889     return q;
 890 }
 891
 892 /*
 893  * Simple remainder.
 894  */
 895 Bignum bigmod(Bignum a, Bignum b)
 896 {
 897     Bignum r = newbn(b[0]);
 898     bigdivmod(a, b, r, NULL);
 899     return r;
 900 }
 901
 902 /*
 903  * Greatest common divisor.
 904  */
 905 Bignum biggcd(Bignum av, Bignum bv)
 906 {
 907     Bignum a = copybn(av);
 908     Bignum b = copybn(bv);
 909
 910     while (bignum_cmp(b, Zero) != 0) {
 911         Bignum t = newbn(b[0]);
 912         bigdivmod(a, b, t, NULL);
 913         while (t[0] > 1 && t[t[0]] == 0)
 914             t[0]--;
 915         freebn(a);
 916         a = b;
 917         b = t;
 918     }
 919
 920     freebn(b);
 921     return a;
 922 }
 923
 924 /*
 925  * Modular inverse, using Euclid's extended algorithm.
 926  */
 927 Bignum modinv(Bignum number, Bignum modulus)
 928 {
 929     Bignum a = copybn(modulus);
 930     Bignum b = copybn(number);
 931     Bignum xp = copybn(Zero);
 932     Bignum x = copybn(One);
 933     int sign = +1;
 934
 935     while (bignum_cmp(b, One) != 0) {
 936         Bignum t = newbn(b[0]);
 937         Bignum q = newbn(a[0]);
 938         bigdivmod(a, b, t, q);
 939         while (t[0] > 1 && t[t[0]] == 0)
 940             t[0]--;
 941         freebn(a);
 942         a = b;
 943         b = t;
 944         t = xp;
 945         xp = x;
 946         x = bigmuladd(q, xp, t);
 947         sign = -sign;
 948         freebn(t);
 949         freebn(q);
 950     }
 951
 952     freebn(b);
 953     freebn(a);
 954     freebn(xp);
 955
 956     /* now we know that sign * x == 1, and that x < modulus */
 957     if (sign < 0) {
 958         /* set a new x to be modulus - x */
 959         Bignum newx = newbn(modulus[0]);
 960         BignumInt carry = 0;
 961         int maxspot = 1;
 962         int i;
 963
 964         for (i = 1; i <= newx[0]; i++) {
 965             BignumInt aword = (i <= modulus[0] ? modulus[i] : 0);
 966             BignumInt bword = (i <= x[0] ? x[i] : 0);
 967             newx[i] = aword - bword - carry;
 968             bword = ~bword;
 969             carry = carry ? (newx[i] >= bword) : (newx[i] > bword);
 970             if (newx[i] != 0)
 971                 maxspot = i;
 972         }
 973         newx[0] = maxspot;
 974         freebn(x);
 975         x = newx;
 976     }
 977
 978     /* and return. */
 979     return x;
 980 }
 981
 982 /*
 983  * Render a bignum into decimal. Return a malloced string holding
 984  * the decimal representation.
 985  */
 986 char *bignum_decimal(Bignum x)
 987 {
 988     int ndigits, ndigit;
 989     int i, iszero;
 990     BignumDblInt carry;
 991     char *ret;
 992     BignumInt *workspace;
 993
 994     /*
 995      * First, estimate the number of digits. Since log(10)/log(2)
 996      * is just greater than 93/28 (the joys of continued fraction
 997      * approximations...) we know that for every 93 bits, we need
 998      * at most 28 digits. This will tell us how much to malloc.
 999      *
1000      * Formally: if x has i bits, that means x is strictly less
1001      * than 2^i. Since 2 is less than 10^(28/93), this is less than
1002      * 10^(28i/93). We need an integer power of ten, so we must
1003      * round up (rounding down might make it less than x again).
1004      * Therefore if we multiply the bit count by 28/93, rounding
1005      * up, we will have enough digits.
1006      */
1007     i = bignum_bitcount(x);
1008     ndigits = (28 * i + 92) / 93;      /* multiply by 28/93 and round up */
1009     ndigits++;                         /* allow for trailing \0 */
1010     ret = snewn(ndigits, char);
1011
1012     /*
1013      * Now allocate some workspace to hold the binary form as we
1014      * repeatedly divide it by ten. Initialise this to the
1015      * big-endian form of the number.
1016      */
1017     workspace = snewn(x[0], BignumInt);
1018     for (i = 0; i < x[0]; i++)
1019         workspace[i] = x[x[0] - i];
1020
1021     /*
1022      * Next, write the decimal number starting with the last digit.
1023      * We use ordinary short division, dividing 10 into the
1024      * workspace.
1025      */
1026     ndigit = ndigits - 1;
1027     ret[ndigit] = '\0';
1028     do {
1029         iszero = 1;
1030         carry = 0;
1031         for (i = 0; i < x[0]; i++) {
1032             carry = (carry << BIGNUM_INT_BITS) + workspace[i];
1033             workspace[i] = (BignumInt) (carry / 10);
1034             if (workspace[i])
1035                 iszero = 0;
1036             carry %= 10;
1037         }
1038         ret[--ndigit] = (char) (carry + '0');
1039     } while (!iszero);
1040
1041     /*
1042      * There's a chance we've fallen short of the start of the
1043      * string. Correct if so.
1044      */
1045     if (ndigit > 0)
1046         memmove(ret, ret + ndigit, ndigits - ndigit);
1047
1048     /*
1049      * Done.
1050      */
1051     sfree(workspace);
1052     return ret;
1053 }