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