sshbn.c

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