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