src/TortoisePlink/SSHBN.C

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