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