sshbn.c

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