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