putty/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 #elif defined _LP64\r
  55 /* 64-bit architectures can do 32x32->64 chunks at a time */\r
  56 typedef unsigned int BignumInt;\r
  57 typedef unsigned long BignumDblInt;\r
  58 #define BIGNUM_INT_MASK  0xFFFFFFFFU\r
  59 #define BIGNUM_TOP_BIT   0x80000000U\r
  60 #define BIGNUM_INT_BITS  32\r
  61 #define MUL_WORD(w1, w2) ((BignumDblInt)w1 * w2)\r
  62 #define DIVMOD_WORD(q, r, hi, lo, w) do { \\r
  63     BignumDblInt n = (((BignumDblInt)hi) << BIGNUM_INT_BITS) | lo; \\r
  64     q = n / w; \\r
  65     r = n % w; \\r
  66 } while (0)\r
  67 #elif defined _LLP64\r
  68 /* 64-bit architectures in which unsigned long is 32 bits, not 64 */\r
  69 typedef unsigned long BignumInt;\r
  70 typedef unsigned long long BignumDblInt;\r
  71 #define BIGNUM_INT_MASK  0xFFFFFFFFUL\r
  72 #define BIGNUM_TOP_BIT   0x80000000UL\r
  73 #define BIGNUM_INT_BITS  32\r
  74 #define MUL_WORD(w1, w2) ((BignumDblInt)w1 * w2)\r
  75 #define DIVMOD_WORD(q, r, hi, lo, w) do { \\r
  76     BignumDblInt n = (((BignumDblInt)hi) << BIGNUM_INT_BITS) | lo; \\r
  77     q = n / w; \\r
  78     r = n % w; \\r
  79 } while (0)\r
  80 #else\r
  81 /* Fallback for all other cases */\r
  82 typedef unsigned short BignumInt;\r
  83 typedef unsigned long BignumDblInt;\r
  84 #define BIGNUM_INT_MASK  0xFFFFU\r
  85 #define BIGNUM_TOP_BIT   0x8000U\r
  86 #define BIGNUM_INT_BITS  16\r
  87 #define MUL_WORD(w1, w2) ((BignumDblInt)w1 * w2)\r
  88 #define DIVMOD_WORD(q, r, hi, lo, w) do { \\r
  89     BignumDblInt n = (((BignumDblInt)hi) << BIGNUM_INT_BITS) | lo; \\r
  90     q = n / w; \\r
  91     r = n % w; \\r
  92 } while (0)\r
  93 #endif\r
  94 \r
  95 #define BIGNUM_INT_BYTES (BIGNUM_INT_BITS / 8)\r
  96 \r
  97 #define BIGNUM_INTERNAL\r
  98 typedef BignumInt *Bignum;\r
  99 \r
 100 #include "ssh.h"\r
 101 \r
 102 BignumInt bnZero[1] = { 0 };\r
 103 BignumInt bnOne[2] = { 1, 1 };\r
 104 \r
 105 /*\r
 106  * The Bignum format is an array of `BignumInt'. The first\r
 107  * element of the array counts the remaining elements. The\r
 108  * remaining elements express the actual number, base 2^BIGNUM_INT_BITS, _least_\r
 109  * significant digit first. (So it's trivial to extract the bit\r
 110  * with value 2^n for any n.)\r
 111  *\r
 112  * All Bignums in this module are positive. Negative numbers must\r
 113  * be dealt with outside it.\r
 114  *\r
 115  * INVARIANT: the most significant word of any Bignum must be\r
 116  * nonzero.\r
 117  */\r
 118 \r
 119 Bignum Zero = bnZero, One = bnOne;\r
 120 \r
 121 static Bignum newbn(int length)\r
 122 {\r
 123     Bignum b = snewn(length + 1, BignumInt);\r
 124     if (!b)\r
 125         abort();                       /* FIXME */\r
 126     memset(b, 0, (length + 1) * sizeof(*b));\r
 127     b[0] = length;\r
 128     return b;\r
 129 }\r
 130 \r
 131 void bn_restore_invariant(Bignum b)\r
 132 {\r
 133     while (b[0] > 1 && b[b[0]] == 0)\r
 134         b[0]--;\r
 135 }\r
 136 \r
 137 Bignum copybn(Bignum orig)\r
 138 {\r
 139     Bignum b = snewn(orig[0] + 1, BignumInt);\r
 140     if (!b)\r
 141         abort();                       /* FIXME */\r
 142     memcpy(b, orig, (orig[0] + 1) * sizeof(*b));\r
 143     return b;\r
 144 }\r
 145 \r
 146 void freebn(Bignum b)\r
 147 {\r
 148     /*\r
 149      * Burn the evidence, just in case.\r
 150      */\r
 151     memset(b, 0, sizeof(b[0]) * (b[0] + 1));\r
 152     sfree(b);\r
 153 }\r
 154 \r
 155 Bignum bn_power_2(int n)\r
 156 {\r
 157     Bignum ret = newbn(n / BIGNUM_INT_BITS + 1);\r
 158     bignum_set_bit(ret, n, 1);\r
 159     return ret;\r
 160 }\r
 161 \r
 162 /*\r
 163  * Internal addition. Sets c = a - b, where 'a', 'b' and 'c' are all\r
 164  * big-endian arrays of 'len' BignumInts. Returns a BignumInt carried\r
 165  * off the top.\r
 166  */\r
 167 static BignumInt internal_add(const BignumInt *a, const BignumInt *b,\r
 168                               BignumInt *c, int len)\r
 169 {\r
 170     int i;\r
 171     BignumDblInt carry = 0;\r
 172 \r
 173     for (i = len-1; i >= 0; i--) {\r
 174         carry += (BignumDblInt)a[i] + b[i];\r
 175         c[i] = (BignumInt)carry;\r
 176         carry >>= BIGNUM_INT_BITS;\r
 177     }\r
 178 \r
 179     return (BignumInt)carry;\r
 180 }\r
 181 \r
 182 /*\r
 183  * Internal subtraction. Sets c = a - b, where 'a', 'b' and 'c' are\r
 184  * all big-endian arrays of 'len' BignumInts. Any borrow from the top\r
 185  * is ignored.\r
 186  */\r
 187 static void internal_sub(const BignumInt *a, const BignumInt *b,\r
 188                          BignumInt *c, int len)\r
 189 {\r
 190     int i;\r
 191     BignumDblInt carry = 1;\r
 192 \r
 193     for (i = len-1; i >= 0; i--) {\r
 194         carry += (BignumDblInt)a[i] + (b[i] ^ BIGNUM_INT_MASK);\r
 195         c[i] = (BignumInt)carry;\r
 196         carry >>= BIGNUM_INT_BITS;\r
 197     }\r
 198 }\r
 199 \r
 200 /*\r
 201  * Compute c = a * b.\r
 202  * Input is in the first len words of a and b.\r
 203  * Result is returned in the first 2*len words of c.\r
 204  *\r
 205  * 'scratch' must point to an array of BignumInt of size at least\r
 206  * mul_compute_scratch(len). (This covers the needs of internal_mul\r
 207  * and all its recursive calls to itself.)\r
 208  */\r
 209 #define KARATSUBA_THRESHOLD 50\r
 210 static int mul_compute_scratch(int len)\r
 211 {\r
 212     int ret = 0;\r
 213     while (len > KARATSUBA_THRESHOLD) {\r
 214         int toplen = len/2, botlen = len - toplen; /* botlen is the bigger */\r
 215         int midlen = botlen + 1;\r
 216         ret += 4*midlen;\r
 217         len = midlen;\r
 218     }\r
 219     return ret;\r
 220 }\r
 221 static void internal_mul(const BignumInt *a, const BignumInt *b,\r
 222                          BignumInt *c, int len, BignumInt *scratch)\r
 223 {\r
 224     if (len > KARATSUBA_THRESHOLD) {\r
 225         int i;\r
 226 \r
 227         /*\r
 228          * Karatsuba divide-and-conquer algorithm. Cut each input in\r
 229          * half, so that it's expressed as two big 'digits' in a giant\r
 230          * base D:\r
 231          *\r
 232          *   a = a_1 D + a_0\r
 233          *   b = b_1 D + b_0\r
 234          *\r
 235          * Then the product is of course\r
 236          *\r
 237          *  ab = a_1 b_1 D^2 + (a_1 b_0 + a_0 b_1) D + a_0 b_0\r
 238          *\r
 239          * and we compute the three coefficients by recursively\r
 240          * calling ourself to do half-length multiplications.\r
 241          *\r
 242          * The clever bit that makes this worth doing is that we only\r
 243          * need _one_ half-length multiplication for the central\r
 244          * coefficient rather than the two that it obviouly looks\r
 245          * like, because we can use a single multiplication to compute\r
 246          *\r
 247          *   (a_1 + a_0) (b_1 + b_0) = a_1 b_1 + a_1 b_0 + a_0 b_1 + a_0 b_0\r
 248          *\r
 249          * and then we subtract the other two coefficients (a_1 b_1\r
 250          * and a_0 b_0) which we were computing anyway.\r
 251          *\r
 252          * Hence we get to multiply two numbers of length N in about\r
 253          * three times as much work as it takes to multiply numbers of\r
 254          * length N/2, which is obviously better than the four times\r
 255          * as much work it would take if we just did a long\r
 256          * conventional multiply.\r
 257          */\r
 258 \r
 259         int toplen = len/2, botlen = len - toplen; /* botlen is the bigger */\r
 260         int midlen = botlen + 1;\r
 261         BignumDblInt carry;\r
 262 #ifdef KARA_DEBUG\r
 263         int i;\r
 264 #endif\r
 265 \r
 266         /*\r
 267          * The coefficients a_1 b_1 and a_0 b_0 just avoid overlapping\r
 268          * in the output array, so we can compute them immediately in\r
 269          * place.\r
 270          */\r
 271 \r
 272 #ifdef KARA_DEBUG\r
 273         printf("a1,a0 = 0x");\r
 274         for (i = 0; i < len; i++) {\r
 275             if (i == toplen) printf(", 0x");\r
 276             printf("%0*x", BIGNUM_INT_BITS/4, a[i]);\r
 277         }\r
 278         printf("\n");\r
 279         printf("b1,b0 = 0x");\r
 280         for (i = 0; i < len; i++) {\r
 281             if (i == toplen) printf(", 0x");\r
 282             printf("%0*x", BIGNUM_INT_BITS/4, b[i]);\r
 283         }\r
 284         printf("\n");\r
 285 #endif\r
 286 \r
 287         /* a_1 b_1 */\r
 288         internal_mul(a, b, c, toplen, scratch);\r
 289 #ifdef KARA_DEBUG\r
 290         printf("a1b1 = 0x");\r
 291         for (i = 0; i < 2*toplen; i++) {\r
 292             printf("%0*x", BIGNUM_INT_BITS/4, c[i]);\r
 293         }\r
 294         printf("\n");\r
 295 #endif\r
 296 \r
 297         /* a_0 b_0 */\r
 298         internal_mul(a + toplen, b + toplen, c + 2*toplen, botlen, scratch);\r
 299 #ifdef KARA_DEBUG\r
 300         printf("a0b0 = 0x");\r
 301         for (i = 0; i < 2*botlen; i++) {\r
 302             printf("%0*x", BIGNUM_INT_BITS/4, c[2*toplen+i]);\r
 303         }\r
 304         printf("\n");\r
 305 #endif\r
 306 \r
 307         /* Zero padding. midlen exceeds toplen by at most 2, so just\r
 308          * zero the first two words of each input and the rest will be\r
 309          * copied over. */\r
 310         scratch[0] = scratch[1] = scratch[midlen] = scratch[midlen+1] = 0;\r
 311 \r
 312         for (i = 0; i < toplen; i++) {\r
 313             scratch[midlen - toplen + i] = a[i]; /* a_1 */\r
 314             scratch[2*midlen - toplen + i] = b[i]; /* b_1 */\r
 315         }\r
 316 \r
 317         /* compute a_1 + a_0 */\r
 318         scratch[0] = internal_add(scratch+1, a+toplen, scratch+1, botlen);\r
 319 #ifdef KARA_DEBUG\r
 320         printf("a1plusa0 = 0x");\r
 321         for (i = 0; i < midlen; i++) {\r
 322             printf("%0*x", BIGNUM_INT_BITS/4, scratch[i]);\r
 323         }\r
 324         printf("\n");\r
 325 #endif\r
 326         /* compute b_1 + b_0 */\r
 327         scratch[midlen] = internal_add(scratch+midlen+1, b+toplen,\r
 328                                        scratch+midlen+1, botlen);\r
 329 #ifdef KARA_DEBUG\r
 330         printf("b1plusb0 = 0x");\r
 331         for (i = 0; i < midlen; i++) {\r
 332             printf("%0*x", BIGNUM_INT_BITS/4, scratch[midlen+i]);\r
 333         }\r
 334         printf("\n");\r
 335 #endif\r
 336 \r
 337         /*\r
 338          * Now we can do the third multiplication.\r
 339          */\r
 340         internal_mul(scratch, scratch + midlen, scratch + 2*midlen, midlen,\r
 341                      scratch + 4*midlen);\r
 342 #ifdef KARA_DEBUG\r
 343         printf("a1plusa0timesb1plusb0 = 0x");\r
 344         for (i = 0; i < 2*midlen; i++) {\r
 345             printf("%0*x", BIGNUM_INT_BITS/4, scratch[2*midlen+i]);\r
 346         }\r
 347         printf("\n");\r
 348 #endif\r
 349 \r
 350         /*\r
 351          * Now we can reuse the first half of 'scratch' to compute the\r
 352          * sum of the outer two coefficients, to subtract from that\r
 353          * product to obtain the middle one.\r
 354          */\r
 355         scratch[0] = scratch[1] = scratch[2] = scratch[3] = 0;\r
 356         for (i = 0; i < 2*toplen; i++)\r
 357             scratch[2*midlen - 2*toplen + i] = c[i];\r
 358         scratch[1] = internal_add(scratch+2, c + 2*toplen,\r
 359                                   scratch+2, 2*botlen);\r
 360 #ifdef KARA_DEBUG\r
 361         printf("a1b1plusa0b0 = 0x");\r
 362         for (i = 0; i < 2*midlen; i++) {\r
 363             printf("%0*x", BIGNUM_INT_BITS/4, scratch[i]);\r
 364         }\r
 365         printf("\n");\r
 366 #endif\r
 367 \r
 368         internal_sub(scratch + 2*midlen, scratch,\r
 369                      scratch + 2*midlen, 2*midlen);\r
 370 #ifdef KARA_DEBUG\r
 371         printf("a1b0plusa0b1 = 0x");\r
 372         for (i = 0; i < 2*midlen; i++) {\r
 373             printf("%0*x", BIGNUM_INT_BITS/4, scratch[2*midlen+i]);\r
 374         }\r
 375         printf("\n");\r
 376 #endif\r
 377 \r
 378         /*\r
 379          * And now all we need to do is to add that middle coefficient\r
 380          * back into the output. We may have to propagate a carry\r
 381          * further up the output, but we can be sure it won't\r
 382          * propagate right the way off the top.\r
 383          */\r
 384         carry = internal_add(c + 2*len - botlen - 2*midlen,\r
 385                              scratch + 2*midlen,\r
 386                              c + 2*len - botlen - 2*midlen, 2*midlen);\r
 387         i = 2*len - botlen - 2*midlen - 1;\r
 388         while (carry) {\r
 389             assert(i >= 0);\r
 390             carry += c[i];\r
 391             c[i] = (BignumInt)carry;\r
 392             carry >>= BIGNUM_INT_BITS;\r
 393             i--;\r
 394         }\r
 395 #ifdef KARA_DEBUG\r
 396         printf("ab = 0x");\r
 397         for (i = 0; i < 2*len; i++) {\r
 398             printf("%0*x", BIGNUM_INT_BITS/4, c[i]);\r
 399         }\r
 400         printf("\n");\r
 401 #endif\r
 402 \r
 403     } else {\r
 404         int i;\r
 405         BignumInt carry;\r
 406         BignumDblInt t;\r
 407         const BignumInt *ap, *bp;\r
 408         BignumInt *cp, *cps;\r
 409 \r
 410         /*\r
 411          * Multiply in the ordinary O(N^2) way.\r
 412          */\r
 413 \r
 414         for (i = 0; i < 2 * len; i++)\r
 415             c[i] = 0;\r
 416 \r
 417         for (cps = c + 2*len, ap = a + len; ap-- > a; cps--) {\r
 418             carry = 0;\r
 419             for (cp = cps, bp = b + len; cp--, bp-- > b ;) {\r
 420                 t = (MUL_WORD(*ap, *bp) + carry) + *cp;\r
 421                 *cp = (BignumInt) t;\r
 422                 carry = (BignumInt)(t >> BIGNUM_INT_BITS);\r
 423             }\r
 424             *cp = carry;\r
 425         }\r
 426     }\r
 427 }\r
 428 \r
 429 /*\r
 430  * Variant form of internal_mul used for the initial step of\r
 431  * Montgomery reduction. Only bothers outputting 'len' words\r
 432  * (everything above that is thrown away).\r
 433  */\r
 434 static void internal_mul_low(const BignumInt *a, const BignumInt *b,\r
 435                              BignumInt *c, int len, BignumInt *scratch)\r
 436 {\r
 437     if (len > KARATSUBA_THRESHOLD) {\r
 438         int i;\r
 439 \r
 440         /*\r
 441          * Karatsuba-aware version of internal_mul_low. As before, we\r
 442          * express each input value as a shifted combination of two\r
 443          * halves:\r
 444          *\r
 445          *   a = a_1 D + a_0\r
 446          *   b = b_1 D + b_0\r
 447          *\r
 448          * Then the full product is, as before,\r
 449          *\r
 450          *  ab = a_1 b_1 D^2 + (a_1 b_0 + a_0 b_1) D + a_0 b_0\r
 451          *\r
 452          * Provided we choose D on the large side (so that a_0 and b_0\r
 453          * are _at least_ as long as a_1 and b_1), we don't need the\r
 454          * topmost term at all, and we only need half of the middle\r
 455          * term. So there's no point in doing the proper Karatsuba\r
 456          * optimisation which computes the middle term using the top\r
 457          * one, because we'd take as long computing the top one as\r
 458          * just computing the middle one directly.\r
 459          *\r
 460          * So instead, we do a much more obvious thing: we call the\r
 461          * fully optimised internal_mul to compute a_0 b_0, and we\r
 462          * recursively call ourself to compute the _bottom halves_ of\r
 463          * a_1 b_0 and a_0 b_1, each of which we add into the result\r
 464          * in the obvious way.\r
 465          *\r
 466          * In other words, there's no actual Karatsuba _optimisation_\r
 467          * in this function; the only benefit in doing it this way is\r
 468          * that we call internal_mul proper for a large part of the\r
 469          * work, and _that_ can optimise its operation.\r
 470          */\r
 471 \r
 472         int toplen = len/2, botlen = len - toplen; /* botlen is the bigger */\r
 473 \r
 474         /*\r
 475          * Scratch space for the various bits and pieces we're going\r
 476          * to be adding together: we need botlen*2 words for a_0 b_0\r
 477          * (though we may end up throwing away its topmost word), and\r
 478          * toplen words for each of a_1 b_0 and a_0 b_1. That adds up\r
 479          * to exactly 2*len.\r
 480          */\r
 481 \r
 482         /* a_0 b_0 */\r
 483         internal_mul(a + toplen, b + toplen, scratch + 2*toplen, botlen,\r
 484                      scratch + 2*len);\r
 485 \r
 486         /* a_1 b_0 */\r
 487         internal_mul_low(a, b + len - toplen, scratch + toplen, toplen,\r
 488                          scratch + 2*len);\r
 489 \r
 490         /* a_0 b_1 */\r
 491         internal_mul_low(a + len - toplen, b, scratch, toplen,\r
 492                          scratch + 2*len);\r
 493 \r
 494         /* Copy the bottom half of the big coefficient into place */\r
 495         for (i = 0; i < botlen; i++)\r
 496             c[toplen + i] = scratch[2*toplen + botlen + i];\r
 497 \r
 498         /* Add the two small coefficients, throwing away the returned carry */\r
 499         internal_add(scratch, scratch + toplen, scratch, toplen);\r
 500 \r
 501         /* And add that to the large coefficient, leaving the result in c. */\r
 502         internal_add(scratch, scratch + 2*toplen + botlen - toplen,\r
 503                      c, toplen);\r
 504 \r
 505     } else {\r
 506         int i;\r
 507         BignumInt carry;\r
 508         BignumDblInt t;\r
 509         const BignumInt *ap, *bp;\r
 510         BignumInt *cp, *cps;\r
 511 \r
 512         /*\r
 513          * Multiply in the ordinary O(N^2) way.\r
 514          */\r
 515 \r
 516         for (i = 0; i < len; i++)\r
 517             c[i] = 0;\r
 518 \r
 519         for (cps = c + len, ap = a + len; ap-- > a; cps--) {\r
 520             carry = 0;\r
 521             for (cp = cps, bp = b + len; bp--, cp-- > c ;) {\r
 522                 t = (MUL_WORD(*ap, *bp) + carry) + *cp;\r
 523                 *cp = (BignumInt) t;\r
 524                 carry = (BignumInt)(t >> BIGNUM_INT_BITS);\r
 525             }\r
 526         }\r
 527     }\r
 528 }\r
 529 \r
 530 /*\r
 531  * Montgomery reduction. Expects x to be a big-endian array of 2*len\r
 532  * BignumInts whose value satisfies 0 <= x < rn (where r = 2^(len *\r
 533  * BIGNUM_INT_BITS) is the Montgomery base). Returns in the same array\r
 534  * a value x' which is congruent to xr^{-1} mod n, and satisfies 0 <=\r
 535  * x' < n.\r
 536  *\r
 537  * 'n' and 'mninv' should be big-endian arrays of 'len' BignumInts\r
 538  * each, containing respectively n and the multiplicative inverse of\r
 539  * -n mod r.\r
 540  *\r
 541  * 'tmp' is an array of BignumInt used as scratch space, of length at\r
 542  * least 3*len + mul_compute_scratch(len).\r
 543  */\r
 544 static void monty_reduce(BignumInt *x, const BignumInt *n,\r
 545                          const BignumInt *mninv, BignumInt *tmp, int len)\r
 546 {\r
 547     int i;\r
 548     BignumInt carry;\r
 549 \r
 550     /*\r
 551      * Multiply x by (-n)^{-1} mod r. This gives us a value m such\r
 552      * that mn is congruent to -x mod r. Hence, mn+x is an exact\r
 553      * multiple of r, and is also (obviously) congruent to x mod n.\r
 554      */\r
 555     internal_mul_low(x + len, mninv, tmp, len, tmp + 3*len);\r
 556 \r
 557     /*\r
 558      * Compute t = (mn+x)/r in ordinary, non-modular, integer\r
 559      * arithmetic. By construction this is exact, and is congruent mod\r
 560      * n to x * r^{-1}, i.e. the answer we want.\r
 561      *\r
 562      * The following multiply leaves that answer in the _most_\r
 563      * significant half of the 'x' array, so then we must shift it\r
 564      * down.\r
 565      */\r
 566     internal_mul(tmp, n, tmp+len, len, tmp + 3*len);\r
 567     carry = internal_add(x, tmp+len, x, 2*len);\r
 568     for (i = 0; i < len; i++)\r
 569         x[len + i] = x[i], x[i] = 0;\r
 570 \r
 571     /*\r
 572      * Reduce t mod n. This doesn't require a full-on division by n,\r
 573      * but merely a test and single optional subtraction, since we can\r
 574      * show that 0 <= t < 2n.\r
 575      *\r
 576      * Proof:\r
 577      *  + we computed m mod r, so 0 <= m < r.\r
 578      *  + so 0 <= mn < rn, obviously\r
 579      *  + hence we only need 0 <= x < rn to guarantee that 0 <= mn+x < 2rn\r
 580      *  + yielding 0 <= (mn+x)/r < 2n as required.\r
 581      */\r
 582     if (!carry) {\r
 583         for (i = 0; i < len; i++)\r
 584             if (x[len + i] != n[i])\r
 585                 break;\r
 586     }\r
 587     if (carry || i >= len || x[len + i] > n[i])\r
 588         internal_sub(x+len, n, x+len, len);\r
 589 }\r
 590 \r
 591 static void internal_add_shifted(BignumInt *number,\r
 592                                  unsigned n, int shift)\r
 593 {\r
 594     int word = 1 + (shift / BIGNUM_INT_BITS);\r
 595     int bshift = shift % BIGNUM_INT_BITS;\r
 596     BignumDblInt addend;\r
 597 \r
 598     addend = (BignumDblInt)n << bshift;\r
 599 \r
 600     while (addend) {\r
 601         addend += number[word];\r
 602         number[word] = (BignumInt) addend & BIGNUM_INT_MASK;\r
 603         addend >>= BIGNUM_INT_BITS;\r
 604         word++;\r
 605     }\r
 606 }\r
 607 \r
 608 /*\r
 609  * Compute a = a % m.\r
 610  * Input in first alen words of a and first mlen words of m.\r
 611  * Output in first alen words of a\r
 612  * (of which first alen-mlen words will be zero).\r
 613  * The MSW of m MUST have its high bit set.\r
 614  * Quotient is accumulated in the `quotient' array, which is a Bignum\r
 615  * rather than the internal bigendian format. Quotient parts are shifted\r
 616  * left by `qshift' before adding into quot.\r
 617  */\r
 618 static void internal_mod(BignumInt *a, int alen,\r
 619                          BignumInt *m, int mlen,\r
 620                          BignumInt *quot, int qshift)\r
 621 {\r
 622     BignumInt m0, m1;\r
 623     unsigned int h;\r
 624     int i, k;\r
 625 \r
 626     m0 = m[0];\r
 627     if (mlen > 1)\r
 628         m1 = m[1];\r
 629     else\r
 630         m1 = 0;\r
 631 \r
 632     for (i = 0; i <= alen - mlen; i++) {\r
 633         BignumDblInt t;\r
 634         unsigned int q, r, c, ai1;\r
 635 \r
 636         if (i == 0) {\r
 637             h = 0;\r
 638         } else {\r
 639             h = a[i - 1];\r
 640             a[i - 1] = 0;\r
 641         }\r
 642 \r
 643         if (i == alen - 1)\r
 644             ai1 = 0;\r
 645         else\r
 646             ai1 = a[i + 1];\r
 647 \r
 648         /* Find q = h:a[i] / m0 */\r
 649         if (h >= m0) {\r
 650             /*\r
 651              * Special case.\r
 652              * \r
 653              * To illustrate it, suppose a BignumInt is 8 bits, and\r
 654              * we are dividing (say) A1:23:45:67 by A1:B2:C3. Then\r
 655              * our initial division will be 0xA123 / 0xA1, which\r
 656              * will give a quotient of 0x100 and a divide overflow.\r
 657              * However, the invariants in this division algorithm\r
 658              * are not violated, since the full number A1:23:... is\r
 659              * _less_ than the quotient prefix A1:B2:... and so the\r
 660              * following correction loop would have sorted it out.\r
 661              * \r
 662              * In this situation we set q to be the largest\r
 663              * quotient we _can_ stomach (0xFF, of course).\r
 664              */\r
 665             q = BIGNUM_INT_MASK;\r
 666         } else {\r
 667             /* Macro doesn't want an array subscript expression passed\r
 668              * into it (see definition), so use a temporary. */\r
 669             BignumInt tmplo = a[i];\r
 670             DIVMOD_WORD(q, r, h, tmplo, m0);\r
 671 \r
 672             /* Refine our estimate of q by looking at\r
 673              h:a[i]:a[i+1] / m0:m1 */\r
 674             t = MUL_WORD(m1, q);\r
 675             if (t > ((BignumDblInt) r << BIGNUM_INT_BITS) + ai1) {\r
 676                 q--;\r
 677                 t -= m1;\r
 678                 r = (r + m0) & BIGNUM_INT_MASK;     /* overflow? */\r
 679                 if (r >= (BignumDblInt) m0 &&\r
 680                     t > ((BignumDblInt) r << BIGNUM_INT_BITS) + ai1) q--;\r
 681             }\r
 682         }\r
 683 \r
 684         /* Subtract q * m from a[i...] */\r
 685         c = 0;\r
 686         for (k = mlen - 1; k >= 0; k--) {\r
 687             t = MUL_WORD(q, m[k]);\r
 688             t += c;\r
 689             c = (unsigned)(t >> BIGNUM_INT_BITS);\r
 690             if ((BignumInt) t > a[i + k])\r
 691                 c++;\r
 692             a[i + k] -= (BignumInt) t;\r
 693         }\r
 694 \r
 695         /* Add back m in case of borrow */\r
 696         if (c != h) {\r
 697             t = 0;\r
 698             for (k = mlen - 1; k >= 0; k--) {\r
 699                 t += m[k];\r
 700                 t += a[i + k];\r
 701                 a[i + k] = (BignumInt) t;\r
 702                 t = t >> BIGNUM_INT_BITS;\r
 703             }\r
 704             q--;\r
 705         }\r
 706         if (quot)\r
 707             internal_add_shifted(quot, q, qshift + BIGNUM_INT_BITS * (alen - mlen - i));\r
 708     }\r
 709 }\r
 710 \r
 711 /*\r
 712  * Compute (base ^ exp) % mod, the pedestrian way.\r
 713  */\r
 714 Bignum modpow_simple(Bignum base_in, Bignum exp, Bignum mod)\r
 715 {\r
 716     BignumInt *a, *b, *n, *m, *scratch;\r
 717     int mshift;\r
 718     int mlen, scratchlen, i, j;\r
 719     Bignum base, result;\r
 720 \r
 721     /*\r
 722      * The most significant word of mod needs to be non-zero. It\r
 723      * should already be, but let's make sure.\r
 724      */\r
 725     assert(mod[mod[0]] != 0);\r
 726 \r
 727     /*\r
 728      * Make sure the base is smaller than the modulus, by reducing\r
 729      * it modulo the modulus if not.\r
 730      */\r
 731     base = bigmod(base_in, mod);\r
 732 \r
 733     /* Allocate m of size mlen, copy mod to m */\r
 734     /* We use big endian internally */\r
 735     mlen = mod[0];\r
 736     m = snewn(mlen, BignumInt);\r
 737     for (j = 0; j < mlen; j++)\r
 738         m[j] = mod[mod[0] - j];\r
 739 \r
 740     /* Shift m left to make msb bit set */\r
 741     for (mshift = 0; mshift < BIGNUM_INT_BITS-1; mshift++)\r
 742         if ((m[0] << mshift) & BIGNUM_TOP_BIT)\r
 743             break;\r
 744     if (mshift) {\r
 745         for (i = 0; i < mlen - 1; i++)\r
 746             m[i] = (m[i] << mshift) | (m[i + 1] >> (BIGNUM_INT_BITS - mshift));\r
 747         m[mlen - 1] = m[mlen - 1] << mshift;\r
 748     }\r
 749 \r
 750     /* Allocate n of size mlen, copy base to n */\r
 751     n = snewn(mlen, BignumInt);\r
 752     i = mlen - base[0];\r
 753     for (j = 0; j < i; j++)\r
 754         n[j] = 0;\r
 755     for (j = 0; j < (int)base[0]; j++)\r
 756         n[i + j] = base[base[0] - j];\r
 757 \r
 758     /* Allocate a and b of size 2*mlen. Set a = 1 */\r
 759     a = snewn(2 * mlen, BignumInt);\r
 760     b = snewn(2 * mlen, BignumInt);\r
 761     for (i = 0; i < 2 * mlen; i++)\r
 762         a[i] = 0;\r
 763     a[2 * mlen - 1] = 1;\r
 764 \r
 765     /* Scratch space for multiplies */\r
 766     scratchlen = mul_compute_scratch(mlen);\r
 767     scratch = snewn(scratchlen, BignumInt);\r
 768 \r
 769     /* Skip leading zero bits of exp. */\r
 770     i = 0;\r
 771     j = BIGNUM_INT_BITS-1;\r
 772     while (i < (int)exp[0] && (exp[exp[0] - i] & (1 << j)) == 0) {\r
 773         j--;\r
 774         if (j < 0) {\r
 775             i++;\r
 776             j = BIGNUM_INT_BITS-1;\r
 777         }\r
 778     }\r
 779 \r
 780     /* Main computation */\r
 781     while (i < (int)exp[0]) {\r
 782         while (j >= 0) {\r
 783             internal_mul(a + mlen, a + mlen, b, mlen, scratch);\r
 784             internal_mod(b, mlen * 2, m, mlen, NULL, 0);\r
 785             if ((exp[exp[0] - i] & (1 << j)) != 0) {\r
 786                 internal_mul(b + mlen, n, a, mlen, scratch);\r
 787                 internal_mod(a, mlen * 2, m, mlen, NULL, 0);\r
 788             } else {\r
 789                 BignumInt *t;\r
 790                 t = a;\r
 791                 a = b;\r
 792                 b = t;\r
 793             }\r
 794             j--;\r
 795         }\r
 796         i++;\r
 797         j = BIGNUM_INT_BITS-1;\r
 798     }\r
 799 \r
 800     /* Fixup result in case the modulus was shifted */\r
 801     if (mshift) {\r
 802         for (i = mlen - 1; i < 2 * mlen - 1; i++)\r
 803             a[i] = (a[i] << mshift) | (a[i + 1] >> (BIGNUM_INT_BITS - mshift));\r
 804         a[2 * mlen - 1] = a[2 * mlen - 1] << mshift;\r
 805         internal_mod(a, mlen * 2, m, mlen, NULL, 0);\r
 806         for (i = 2 * mlen - 1; i >= mlen; i--)\r
 807             a[i] = (a[i] >> mshift) | (a[i - 1] << (BIGNUM_INT_BITS - mshift));\r
 808     }\r
 809 \r
 810     /* Copy result to buffer */\r
 811     result = newbn(mod[0]);\r
 812     for (i = 0; i < mlen; i++)\r
 813         result[result[0] - i] = a[i + mlen];\r
 814     while (result[0] > 1 && result[result[0]] == 0)\r
 815         result[0]--;\r
 816 \r
 817     /* Free temporary arrays */\r
 818     for (i = 0; i < 2 * mlen; i++)\r
 819         a[i] = 0;\r
 820     sfree(a);\r
 821     for (i = 0; i < scratchlen; i++)\r
 822         scratch[i] = 0;\r
 823     sfree(scratch);\r
 824     for (i = 0; i < 2 * mlen; i++)\r
 825         b[i] = 0;\r
 826     sfree(b);\r
 827     for (i = 0; i < mlen; i++)\r
 828         m[i] = 0;\r
 829     sfree(m);\r
 830     for (i = 0; i < mlen; i++)\r
 831         n[i] = 0;\r
 832     sfree(n);\r
 833 \r
 834     freebn(base);\r
 835 \r
 836     return result;\r
 837 }\r
 838 \r
 839 /*\r
 840  * Compute (base ^ exp) % mod. Uses the Montgomery multiplication\r
 841  * technique where possible, falling back to modpow_simple otherwise.\r
 842  */\r
 843 Bignum modpow(Bignum base_in, Bignum exp, Bignum mod)\r
 844 {\r
 845     BignumInt *a, *b, *x, *n, *mninv, *scratch;\r
 846     int len, scratchlen, i, j;\r
 847     Bignum base, base2, r, rn, inv, result;\r
 848 \r
 849     /*\r
 850      * The most significant word of mod needs to be non-zero. It\r
 851      * should already be, but let's make sure.\r
 852      */\r
 853     assert(mod[mod[0]] != 0);\r
 854 \r
 855     /*\r
 856      * mod had better be odd, or we can't do Montgomery multiplication\r
 857      * using a power of two at all.\r
 858      */\r
 859     if (!(mod[1] & 1))\r
 860         return modpow_simple(base_in, exp, mod);\r
 861 \r
 862     /*\r
 863      * Make sure the base is smaller than the modulus, by reducing\r
 864      * it modulo the modulus if not.\r
 865      */\r
 866     base = bigmod(base_in, mod);\r
 867 \r
 868     /*\r
 869      * Compute the inverse of n mod r, for monty_reduce. (In fact we\r
 870      * want the inverse of _minus_ n mod r, but we'll sort that out\r
 871      * below.)\r
 872      */\r
 873     len = mod[0];\r
 874     r = bn_power_2(BIGNUM_INT_BITS * len);\r
 875     inv = modinv(mod, r);\r
 876 \r
 877     /*\r
 878      * Multiply the base by r mod n, to get it into Montgomery\r
 879      * representation.\r
 880      */\r
 881     base2 = modmul(base, r, mod);\r
 882     freebn(base);\r
 883     base = base2;\r
 884 \r
 885     rn = bigmod(r, mod);               /* r mod n, i.e. Montgomerified 1 */\r
 886 \r
 887     freebn(r);                         /* won't need this any more */\r
 888 \r
 889     /*\r
 890      * Set up internal arrays of the right lengths, in big-endian\r
 891      * format, containing the base, the modulus, and the modulus's\r
 892      * inverse.\r
 893      */\r
 894     n = snewn(len, BignumInt);\r
 895     for (j = 0; j < len; j++)\r
 896         n[len - 1 - j] = mod[j + 1];\r
 897 \r
 898     mninv = snewn(len, BignumInt);\r
 899     for (j = 0; j < len; j++)\r
 900         mninv[len - 1 - j] = (j < (int)inv[0] ? inv[j + 1] : 0);\r
 901     freebn(inv);         /* we don't need this copy of it any more */\r
 902     /* Now negate mninv mod r, so it's the inverse of -n rather than +n. */\r
 903     x = snewn(len, BignumInt);\r
 904     for (j = 0; j < len; j++)\r
 905         x[j] = 0;\r
 906     internal_sub(x, mninv, mninv, len);\r
 907 \r
 908     /* x = snewn(len, BignumInt); */ /* already done above */\r
 909     for (j = 0; j < len; j++)\r
 910         x[len - 1 - j] = (j < (int)base[0] ? base[j + 1] : 0);\r
 911     freebn(base);        /* we don't need this copy of it any more */\r
 912 \r
 913     a = snewn(2*len, BignumInt);\r
 914     b = snewn(2*len, BignumInt);\r
 915     for (j = 0; j < len; j++)\r
 916         a[2*len - 1 - j] = (j < (int)rn[0] ? rn[j + 1] : 0);\r
 917     freebn(rn);\r
 918 \r
 919     /* Scratch space for multiplies */\r
 920     scratchlen = 3*len + mul_compute_scratch(len);\r
 921     scratch = snewn(scratchlen, BignumInt);\r
 922 \r
 923     /* Skip leading zero bits of exp. */\r
 924     i = 0;\r
 925     j = BIGNUM_INT_BITS-1;\r
 926     while (i < (int)exp[0] && (exp[exp[0] - i] & (1 << j)) == 0) {\r
 927         j--;\r
 928         if (j < 0) {\r
 929             i++;\r
 930             j = BIGNUM_INT_BITS-1;\r
 931         }\r
 932     }\r
 933 \r
 934     /* Main computation */\r
 935     while (i < (int)exp[0]) {\r
 936         while (j >= 0) {\r
 937             internal_mul(a + len, a + len, b, len, scratch);\r
 938             monty_reduce(b, n, mninv, scratch, len);\r
 939             if ((exp[exp[0] - i] & (1 << j)) != 0) {\r
 940                 internal_mul(b + len, x, a, len,  scratch);\r
 941                 monty_reduce(a, n, mninv, scratch, len);\r
 942             } else {\r
 943                 BignumInt *t;\r
 944                 t = a;\r
 945                 a = b;\r
 946                 b = t;\r
 947             }\r
 948             j--;\r
 949         }\r
 950         i++;\r
 951         j = BIGNUM_INT_BITS-1;\r
 952     }\r
 953 \r
 954     /*\r
 955      * Final monty_reduce to get back from the adjusted Montgomery\r
 956      * representation.\r
 957      */\r
 958     monty_reduce(a, n, mninv, scratch, len);\r
 959 \r
 960     /* Copy result to buffer */\r
 961     result = newbn(mod[0]);\r
 962     for (i = 0; i < len; i++)\r
 963         result[result[0] - i] = a[i + len];\r
 964     while (result[0] > 1 && result[result[0]] == 0)\r
 965         result[0]--;\r
 966 \r
 967     /* Free temporary arrays */\r
 968     for (i = 0; i < scratchlen; i++)\r
 969         scratch[i] = 0;\r
 970     sfree(scratch);\r
 971     for (i = 0; i < 2 * len; i++)\r
 972         a[i] = 0;\r
 973     sfree(a);\r
 974     for (i = 0; i < 2 * len; i++)\r
 975         b[i] = 0;\r
 976     sfree(b);\r
 977     for (i = 0; i < len; i++)\r
 978         mninv[i] = 0;\r
 979     sfree(mninv);\r
 980     for (i = 0; i < len; i++)\r
 981         n[i] = 0;\r
 982     sfree(n);\r
 983     for (i = 0; i < len; i++)\r
 984         x[i] = 0;\r
 985     sfree(x);\r
 986 \r
 987     return result;\r
 988 }\r
 989 \r
 990 /*\r
 991  * Compute (p * q) % mod.\r
 992  * The most significant word of mod MUST be non-zero.\r
 993  * We assume that the result array is the same size as the mod array.\r
 994  */\r
 995 Bignum modmul(Bignum p, Bignum q, Bignum mod)\r
 996 {\r
 997     BignumInt *a, *n, *m, *o, *scratch;\r
 998     int mshift, scratchlen;\r
 999     int pqlen, mlen, rlen, i, j;\r
1000     Bignum result;\r
1001 \r
1002     /* Allocate m of size mlen, copy mod to m */\r
1003     /* We use big endian internally */\r
1004     mlen = mod[0];\r
1005     m = snewn(mlen, BignumInt);\r
1006     for (j = 0; j < mlen; j++)\r
1007         m[j] = mod[mod[0] - j];\r
1008 \r
1009     /* Shift m left to make msb bit set */\r
1010     for (mshift = 0; mshift < BIGNUM_INT_BITS-1; mshift++)\r
1011         if ((m[0] << mshift) & BIGNUM_TOP_BIT)\r
1012             break;\r
1013     if (mshift) {\r
1014         for (i = 0; i < mlen - 1; i++)\r
1015             m[i] = (m[i] << mshift) | (m[i + 1] >> (BIGNUM_INT_BITS - mshift));\r
1016         m[mlen - 1] = m[mlen - 1] << mshift;\r
1017     }\r
1018 \r
1019     pqlen = (p[0] > q[0] ? p[0] : q[0]);\r
1020 \r
1021     /* Allocate n of size pqlen, copy p to n */\r
1022     n = snewn(pqlen, BignumInt);\r
1023     i = pqlen - p[0];\r
1024     for (j = 0; j < i; j++)\r
1025         n[j] = 0;\r
1026     for (j = 0; j < (int)p[0]; j++)\r
1027         n[i + j] = p[p[0] - j];\r
1028 \r
1029     /* Allocate o of size pqlen, copy q to o */\r
1030     o = snewn(pqlen, BignumInt);\r
1031     i = pqlen - q[0];\r
1032     for (j = 0; j < i; j++)\r
1033         o[j] = 0;\r
1034     for (j = 0; j < (int)q[0]; j++)\r
1035         o[i + j] = q[q[0] - j];\r
1036 \r
1037     /* Allocate a of size 2*pqlen for result */\r
1038     a = snewn(2 * pqlen, BignumInt);\r
1039 \r
1040     /* Scratch space for multiplies */\r
1041     scratchlen = mul_compute_scratch(pqlen);\r
1042     scratch = snewn(scratchlen, BignumInt);\r
1043 \r
1044     /* Main computation */\r
1045     internal_mul(n, o, a, pqlen, scratch);\r
1046     internal_mod(a, pqlen * 2, m, mlen, NULL, 0);\r
1047 \r
1048     /* Fixup result in case the modulus was shifted */\r
1049     if (mshift) {\r
1050         for (i = 2 * pqlen - mlen - 1; i < 2 * pqlen - 1; i++)\r
1051             a[i] = (a[i] << mshift) | (a[i + 1] >> (BIGNUM_INT_BITS - mshift));\r
1052         a[2 * pqlen - 1] = a[2 * pqlen - 1] << mshift;\r
1053         internal_mod(a, pqlen * 2, m, mlen, NULL, 0);\r
1054         for (i = 2 * pqlen - 1; i >= 2 * pqlen - mlen; i--)\r
1055             a[i] = (a[i] >> mshift) | (a[i - 1] << (BIGNUM_INT_BITS - mshift));\r
1056     }\r
1057 \r
1058     /* Copy result to buffer */\r
1059     rlen = (mlen < pqlen * 2 ? mlen : pqlen * 2);\r
1060     result = newbn(rlen);\r
1061     for (i = 0; i < rlen; i++)\r
1062         result[result[0] - i] = a[i + 2 * pqlen - rlen];\r
1063     while (result[0] > 1 && result[result[0]] == 0)\r
1064         result[0]--;\r
1065 \r
1066     /* Free temporary arrays */\r
1067     for (i = 0; i < scratchlen; i++)\r
1068         scratch[i] = 0;\r
1069     sfree(scratch);\r
1070     for (i = 0; i < 2 * pqlen; i++)\r
1071         a[i] = 0;\r
1072     sfree(a);\r
1073     for (i = 0; i < mlen; i++)\r
1074         m[i] = 0;\r
1075     sfree(m);\r
1076     for (i = 0; i < pqlen; i++)\r
1077         n[i] = 0;\r
1078     sfree(n);\r
1079     for (i = 0; i < pqlen; i++)\r
1080         o[i] = 0;\r
1081     sfree(o);\r
1082 \r
1083     return result;\r
1084 }\r
1085 \r
1086 /*\r
1087  * Compute p % mod.\r
1088  * The most significant word of mod MUST be non-zero.\r
1089  * We assume that the result array is the same size as the mod array.\r
1090  * We optionally write out a quotient if `quotient' is non-NULL.\r
1091  * We can avoid writing out the result if `result' is NULL.\r
1092  */\r
1093 static void bigdivmod(Bignum p, Bignum mod, Bignum result, Bignum quotient)\r
1094 {\r
1095     BignumInt *n, *m;\r
1096     int mshift;\r
1097     int plen, mlen, i, j;\r
1098 \r
1099     /* Allocate m of size mlen, copy mod to m */\r
1100     /* We use big endian internally */\r
1101     mlen = mod[0];\r
1102     m = snewn(mlen, BignumInt);\r
1103     for (j = 0; j < mlen; j++)\r
1104         m[j] = mod[mod[0] - j];\r
1105 \r
1106     /* Shift m left to make msb bit set */\r
1107     for (mshift = 0; mshift < BIGNUM_INT_BITS-1; mshift++)\r
1108         if ((m[0] << mshift) & BIGNUM_TOP_BIT)\r
1109             break;\r
1110     if (mshift) {\r
1111         for (i = 0; i < mlen - 1; i++)\r
1112             m[i] = (m[i] << mshift) | (m[i + 1] >> (BIGNUM_INT_BITS - mshift));\r
1113         m[mlen - 1] = m[mlen - 1] << mshift;\r
1114     }\r
1115 \r
1116     plen = p[0];\r
1117     /* Ensure plen > mlen */\r
1118     if (plen <= mlen)\r
1119         plen = mlen + 1;\r
1120 \r
1121     /* Allocate n of size plen, copy p to n */\r
1122     n = snewn(plen, BignumInt);\r
1123     for (j = 0; j < plen; j++)\r
1124         n[j] = 0;\r
1125     for (j = 1; j <= (int)p[0]; j++)\r
1126         n[plen - j] = p[j];\r
1127 \r
1128     /* Main computation */\r
1129     internal_mod(n, plen, m, mlen, quotient, mshift);\r
1130 \r
1131     /* Fixup result in case the modulus was shifted */\r
1132     if (mshift) {\r
1133         for (i = plen - mlen - 1; i < plen - 1; i++)\r
1134             n[i] = (n[i] << mshift) | (n[i + 1] >> (BIGNUM_INT_BITS - mshift));\r
1135         n[plen - 1] = n[plen - 1] << mshift;\r
1136         internal_mod(n, plen, m, mlen, quotient, 0);\r
1137         for (i = plen - 1; i >= plen - mlen; i--)\r
1138             n[i] = (n[i] >> mshift) | (n[i - 1] << (BIGNUM_INT_BITS - mshift));\r
1139     }\r
1140 \r
1141     /* Copy result to buffer */\r
1142     if (result) {\r
1143         for (i = 1; i <= (int)result[0]; i++) {\r
1144             int j = plen - i;\r
1145             result[i] = j >= 0 ? n[j] : 0;\r
1146         }\r
1147     }\r
1148 \r
1149     /* Free temporary arrays */\r
1150     for (i = 0; i < mlen; i++)\r
1151         m[i] = 0;\r
1152     sfree(m);\r
1153     for (i = 0; i < plen; i++)\r
1154         n[i] = 0;\r
1155     sfree(n);\r
1156 }\r
1157 \r
1158 /*\r
1159  * Decrement a number.\r
1160  */\r
1161 void decbn(Bignum bn)\r
1162 {\r
1163     int i = 1;\r
1164     while (i < (int)bn[0] && bn[i] == 0)\r
1165         bn[i++] = BIGNUM_INT_MASK;\r
1166     bn[i]--;\r
1167 }\r
1168 \r
1169 Bignum bignum_from_bytes(const unsigned char *data, int nbytes)\r
1170 {\r
1171     Bignum result;\r
1172     int w, i;\r
1173 \r
1174     w = (nbytes + BIGNUM_INT_BYTES - 1) / BIGNUM_INT_BYTES; /* bytes->words */\r
1175 \r
1176     result = newbn(w);\r
1177     for (i = 1; i <= w; i++)\r
1178         result[i] = 0;\r
1179     for (i = nbytes; i--;) {\r
1180         unsigned char byte = *data++;\r
1181         result[1 + i / BIGNUM_INT_BYTES] |= byte << (8*i % BIGNUM_INT_BITS);\r
1182     }\r
1183 \r
1184     while (result[0] > 1 && result[result[0]] == 0)\r
1185         result[0]--;\r
1186     return result;\r
1187 }\r
1188 \r
1189 /*\r
1190  * Read an SSH-1-format bignum from a data buffer. Return the number\r
1191  * of bytes consumed, or -1 if there wasn't enough data.\r
1192  */\r
1193 int ssh1_read_bignum(const unsigned char *data, int len, Bignum * result)\r
1194 {\r
1195     const unsigned char *p = data;\r
1196     int i;\r
1197     int w, b;\r
1198 \r
1199     if (len < 2)\r
1200         return -1;\r
1201 \r
1202     w = 0;\r
1203     for (i = 0; i < 2; i++)\r
1204         w = (w << 8) + *p++;\r
1205     b = (w + 7) / 8;                   /* bits -> bytes */\r
1206 \r
1207     if (len < b+2)\r
1208         return -1;\r
1209 \r
1210     if (!result)                       /* just return length */\r
1211         return b + 2;\r
1212 \r
1213     *result = bignum_from_bytes(p, b);\r
1214 \r
1215     return p + b - data;\r
1216 }\r
1217 \r
1218 /*\r
1219  * Return the bit count of a bignum, for SSH-1 encoding.\r
1220  */\r
1221 int bignum_bitcount(Bignum bn)\r
1222 {\r
1223     int bitcount = bn[0] * BIGNUM_INT_BITS - 1;\r
1224     while (bitcount >= 0\r
1225            && (bn[bitcount / BIGNUM_INT_BITS + 1] >> (bitcount % BIGNUM_INT_BITS)) == 0) bitcount--;\r
1226     return bitcount + 1;\r
1227 }\r
1228 \r
1229 /*\r
1230  * Return the byte length of a bignum when SSH-1 encoded.\r
1231  */\r
1232 int ssh1_bignum_length(Bignum bn)\r
1233 {\r
1234     return 2 + (bignum_bitcount(bn) + 7) / 8;\r
1235 }\r
1236 \r
1237 /*\r
1238  * Return the byte length of a bignum when SSH-2 encoded.\r
1239  */\r
1240 int ssh2_bignum_length(Bignum bn)\r
1241 {\r
1242     return 4 + (bignum_bitcount(bn) + 8) / 8;\r
1243 }\r
1244 \r
1245 /*\r
1246  * Return a byte from a bignum; 0 is least significant, etc.\r
1247  */\r
1248 int bignum_byte(Bignum bn, int i)\r
1249 {\r
1250     if (i >= (int)(BIGNUM_INT_BYTES * bn[0]))\r
1251         return 0;                      /* beyond the end */\r
1252     else\r
1253         return (bn[i / BIGNUM_INT_BYTES + 1] >>\r
1254                 ((i % BIGNUM_INT_BYTES)*8)) & 0xFF;\r
1255 }\r
1256 \r
1257 /*\r
1258  * Return a bit from a bignum; 0 is least significant, etc.\r
1259  */\r
1260 int bignum_bit(Bignum bn, int i)\r
1261 {\r
1262     if (i >= (int)(BIGNUM_INT_BITS * bn[0]))\r
1263         return 0;                      /* beyond the end */\r
1264     else\r
1265         return (bn[i / BIGNUM_INT_BITS + 1] >> (i % BIGNUM_INT_BITS)) & 1;\r
1266 }\r
1267 \r
1268 /*\r
1269  * Set a bit in a bignum; 0 is least significant, etc.\r
1270  */\r
1271 void bignum_set_bit(Bignum bn, int bitnum, int value)\r
1272 {\r
1273     if (bitnum >= (int)(BIGNUM_INT_BITS * bn[0]))\r
1274         abort();                       /* beyond the end */\r
1275     else {\r
1276         int v = bitnum / BIGNUM_INT_BITS + 1;\r
1277         int mask = 1 << (bitnum % BIGNUM_INT_BITS);\r
1278         if (value)\r
1279             bn[v] |= mask;\r
1280         else\r
1281             bn[v] &= ~mask;\r
1282     }\r
1283 }\r
1284 \r
1285 /*\r
1286  * Write a SSH-1-format bignum into a buffer. It is assumed the\r
1287  * buffer is big enough. Returns the number of bytes used.\r
1288  */\r
1289 int ssh1_write_bignum(void *data, Bignum bn)\r
1290 {\r
1291     unsigned char *p = data;\r
1292     int len = ssh1_bignum_length(bn);\r
1293     int i;\r
1294     int bitc = bignum_bitcount(bn);\r
1295 \r
1296     *p++ = (bitc >> 8) & 0xFF;\r
1297     *p++ = (bitc) & 0xFF;\r
1298     for (i = len - 2; i--;)\r
1299         *p++ = bignum_byte(bn, i);\r
1300     return len;\r
1301 }\r
1302 \r
1303 /*\r
1304  * Compare two bignums. Returns like strcmp.\r
1305  */\r
1306 int bignum_cmp(Bignum a, Bignum b)\r
1307 {\r
1308     int amax = a[0], bmax = b[0];\r
1309     int i = (amax > bmax ? amax : bmax);\r
1310     while (i) {\r
1311         BignumInt aval = (i > amax ? 0 : a[i]);\r
1312         BignumInt bval = (i > bmax ? 0 : b[i]);\r
1313         if (aval < bval)\r
1314             return -1;\r
1315         if (aval > bval)\r
1316             return +1;\r
1317         i--;\r
1318     }\r
1319     return 0;\r
1320 }\r
1321 \r
1322 /*\r
1323  * Right-shift one bignum to form another.\r
1324  */\r
1325 Bignum bignum_rshift(Bignum a, int shift)\r
1326 {\r
1327     Bignum ret;\r
1328     int i, shiftw, shiftb, shiftbb, bits;\r
1329     BignumInt ai, ai1;\r
1330 \r
1331     bits = bignum_bitcount(a) - shift;\r
1332     ret = newbn((bits + BIGNUM_INT_BITS - 1) / BIGNUM_INT_BITS);\r
1333 \r
1334     if (ret) {\r
1335         shiftw = shift / BIGNUM_INT_BITS;\r
1336         shiftb = shift % BIGNUM_INT_BITS;\r
1337         shiftbb = BIGNUM_INT_BITS - shiftb;\r
1338 \r
1339         ai1 = a[shiftw + 1];\r
1340         for (i = 1; i <= (int)ret[0]; i++) {\r
1341             ai = ai1;\r
1342             ai1 = (i + shiftw + 1 <= (int)a[0] ? a[i + shiftw + 1] : 0);\r
1343             ret[i] = ((ai >> shiftb) | (ai1 << shiftbb)) & BIGNUM_INT_MASK;\r
1344         }\r
1345     }\r
1346 \r
1347     return ret;\r
1348 }\r
1349 \r
1350 /*\r
1351  * Non-modular multiplication and addition.\r
1352  */\r
1353 Bignum bigmuladd(Bignum a, Bignum b, Bignum addend)\r
1354 {\r
1355     int alen = a[0], blen = b[0];\r
1356     int mlen = (alen > blen ? alen : blen);\r
1357     int rlen, i, maxspot;\r
1358     int wslen;\r
1359     BignumInt *workspace;\r
1360     Bignum ret;\r
1361 \r
1362     /* mlen space for a, mlen space for b, 2*mlen for result,\r
1363      * plus scratch space for multiplication */\r
1364     wslen = mlen * 4 + mul_compute_scratch(mlen);\r
1365     workspace = snewn(wslen, BignumInt);\r
1366     for (i = 0; i < mlen; i++) {\r
1367         workspace[0 * mlen + i] = (mlen - i <= (int)a[0] ? a[mlen - i] : 0);\r
1368         workspace[1 * mlen + i] = (mlen - i <= (int)b[0] ? b[mlen - i] : 0);\r
1369     }\r
1370 \r
1371     internal_mul(workspace + 0 * mlen, workspace + 1 * mlen,\r
1372                  workspace + 2 * mlen, mlen, workspace + 4 * mlen);\r
1373 \r
1374     /* now just copy the result back */\r
1375     rlen = alen + blen + 1;\r
1376     if (addend && rlen <= (int)addend[0])\r
1377         rlen = addend[0] + 1;\r
1378     ret = newbn(rlen);\r
1379     maxspot = 0;\r
1380     for (i = 1; i <= (int)ret[0]; i++) {\r
1381         ret[i] = (i <= 2 * mlen ? workspace[4 * mlen - i] : 0);\r
1382         if (ret[i] != 0)\r
1383             maxspot = i;\r
1384     }\r
1385     ret[0] = maxspot;\r
1386 \r
1387     /* now add in the addend, if any */\r
1388     if (addend) {\r
1389         BignumDblInt carry = 0;\r
1390         for (i = 1; i <= rlen; i++) {\r
1391             carry += (i <= (int)ret[0] ? ret[i] : 0);\r
1392             carry += (i <= (int)addend[0] ? addend[i] : 0);\r
1393             ret[i] = (BignumInt) carry & BIGNUM_INT_MASK;\r
1394             carry >>= BIGNUM_INT_BITS;\r
1395             if (ret[i] != 0 && i > maxspot)\r
1396                 maxspot = i;\r
1397         }\r
1398     }\r
1399     ret[0] = maxspot;\r
1400 \r
1401     for (i = 0; i < wslen; i++)\r
1402         workspace[i] = 0;\r
1403     sfree(workspace);\r
1404     return ret;\r
1405 }\r
1406 \r
1407 /*\r
1408  * Non-modular multiplication.\r
1409  */\r
1410 Bignum bigmul(Bignum a, Bignum b)\r
1411 {\r
1412     return bigmuladd(a, b, NULL);\r
1413 }\r
1414 \r
1415 /*\r
1416  * Simple addition.\r
1417  */\r
1418 Bignum bigadd(Bignum a, Bignum b)\r
1419 {\r
1420     int alen = a[0], blen = b[0];\r
1421     int rlen = (alen > blen ? alen : blen) + 1;\r
1422     int i, maxspot;\r
1423     Bignum ret;\r
1424     BignumDblInt carry;\r
1425 \r
1426     ret = newbn(rlen);\r
1427 \r
1428     carry = 0;\r
1429     maxspot = 0;\r
1430     for (i = 1; i <= rlen; i++) {\r
1431         carry += (i <= (int)a[0] ? a[i] : 0);\r
1432         carry += (i <= (int)b[0] ? b[i] : 0);\r
1433         ret[i] = (BignumInt) carry & BIGNUM_INT_MASK;\r
1434         carry >>= BIGNUM_INT_BITS;\r
1435         if (ret[i] != 0 && i > maxspot)\r
1436             maxspot = i;\r
1437     }\r
1438     ret[0] = maxspot;\r
1439 \r
1440     return ret;\r
1441 }\r
1442 \r
1443 /*\r
1444  * Subtraction. Returns a-b, or NULL if the result would come out\r
1445  * negative (recall that this entire bignum module only handles\r
1446  * positive numbers).\r
1447  */\r
1448 Bignum bigsub(Bignum a, Bignum b)\r
1449 {\r
1450     int alen = a[0], blen = b[0];\r
1451     int rlen = (alen > blen ? alen : blen);\r
1452     int i, maxspot;\r
1453     Bignum ret;\r
1454     BignumDblInt carry;\r
1455 \r
1456     ret = newbn(rlen);\r
1457 \r
1458     carry = 1;\r
1459     maxspot = 0;\r
1460     for (i = 1; i <= rlen; i++) {\r
1461         carry += (i <= (int)a[0] ? a[i] : 0);\r
1462         carry += (i <= (int)b[0] ? b[i] ^ BIGNUM_INT_MASK : BIGNUM_INT_MASK);\r
1463         ret[i] = (BignumInt) carry & BIGNUM_INT_MASK;\r
1464         carry >>= BIGNUM_INT_BITS;\r
1465         if (ret[i] != 0 && i > maxspot)\r
1466             maxspot = i;\r
1467     }\r
1468     ret[0] = maxspot;\r
1469 \r
1470     if (!carry) {\r
1471         freebn(ret);\r
1472         return NULL;\r
1473     }\r
1474 \r
1475     return ret;\r
1476 }\r
1477 \r
1478 /*\r
1479  * Create a bignum which is the bitmask covering another one. That\r
1480  * is, the smallest integer which is >= N and is also one less than\r
1481  * a power of two.\r
1482  */\r
1483 Bignum bignum_bitmask(Bignum n)\r
1484 {\r
1485     Bignum ret = copybn(n);\r
1486     int i;\r
1487     BignumInt j;\r
1488 \r
1489     i = ret[0];\r
1490     while (n[i] == 0 && i > 0)\r
1491         i--;\r
1492     if (i <= 0)\r
1493         return ret;                    /* input was zero */\r
1494     j = 1;\r
1495     while (j < n[i])\r
1496         j = 2 * j + 1;\r
1497     ret[i] = j;\r
1498     while (--i > 0)\r
1499         ret[i] = BIGNUM_INT_MASK;\r
1500     return ret;\r
1501 }\r
1502 \r
1503 /*\r
1504  * Convert a (max 32-bit) long into a bignum.\r
1505  */\r
1506 Bignum bignum_from_long(unsigned long nn)\r
1507 {\r
1508     Bignum ret;\r
1509     BignumDblInt n = nn;\r
1510 \r
1511     ret = newbn(3);\r
1512     ret[1] = (BignumInt)(n & BIGNUM_INT_MASK);\r
1513     ret[2] = (BignumInt)((n >> BIGNUM_INT_BITS) & BIGNUM_INT_MASK);\r
1514     ret[3] = 0;\r
1515     ret[0] = (ret[2]  ? 2 : 1);\r
1516     return ret;\r
1517 }\r
1518 \r
1519 /*\r
1520  * Add a long to a bignum.\r
1521  */\r
1522 Bignum bignum_add_long(Bignum number, unsigned long addendx)\r
1523 {\r
1524     Bignum ret = newbn(number[0] + 1);\r
1525     int i, maxspot = 0;\r
1526     BignumDblInt carry = 0, addend = addendx;\r
1527 \r
1528     for (i = 1; i <= (int)ret[0]; i++) {\r
1529         carry += addend & BIGNUM_INT_MASK;\r
1530         carry += (i <= (int)number[0] ? number[i] : 0);\r
1531         addend >>= BIGNUM_INT_BITS;\r
1532         ret[i] = (BignumInt) carry & BIGNUM_INT_MASK;\r
1533         carry >>= BIGNUM_INT_BITS;\r
1534         if (ret[i] != 0)\r
1535             maxspot = i;\r
1536     }\r
1537     ret[0] = maxspot;\r
1538     return ret;\r
1539 }\r
1540 \r
1541 /*\r
1542  * Compute the residue of a bignum, modulo a (max 16-bit) short.\r
1543  */\r
1544 unsigned short bignum_mod_short(Bignum number, unsigned short modulus)\r
1545 {\r
1546     BignumDblInt mod, r;\r
1547     int i;\r
1548 \r
1549     r = 0;\r
1550     mod = modulus;\r
1551     for (i = number[0]; i > 0; i--)\r
1552         r = (r * (BIGNUM_TOP_BIT % mod) * 2 + number[i] % mod) % mod;\r
1553     return (unsigned short) r;\r
1554 }\r
1555 \r
1556 #ifdef DEBUG\r
1557 void diagbn(char *prefix, Bignum md)\r
1558 {\r
1559     int i, nibbles, morenibbles;\r
1560     static const char hex[] = "0123456789ABCDEF";\r
1561 \r
1562     debug(("%s0x", prefix ? prefix : ""));\r
1563 \r
1564     nibbles = (3 + bignum_bitcount(md)) / 4;\r
1565     if (nibbles < 1)\r
1566         nibbles = 1;\r
1567     morenibbles = 4 * md[0] - nibbles;\r
1568     for (i = 0; i < morenibbles; i++)\r
1569         debug(("-"));\r
1570     for (i = nibbles; i--;)\r
1571         debug(("%c",\r
1572                hex[(bignum_byte(md, i / 2) >> (4 * (i % 2))) & 0xF]));\r
1573 \r
1574     if (prefix)\r
1575         debug(("\n"));\r
1576 }\r
1577 #endif\r
1578 \r
1579 /*\r
1580  * Simple division.\r
1581  */\r
1582 Bignum bigdiv(Bignum a, Bignum b)\r
1583 {\r
1584     Bignum q = newbn(a[0]);\r
1585     bigdivmod(a, b, NULL, q);\r
1586     return q;\r
1587 }\r
1588 \r
1589 /*\r
1590  * Simple remainder.\r
1591  */\r
1592 Bignum bigmod(Bignum a, Bignum b)\r
1593 {\r
1594     Bignum r = newbn(b[0]);\r
1595     bigdivmod(a, b, r, NULL);\r
1596     return r;\r
1597 }\r
1598 \r
1599 /*\r
1600  * Greatest common divisor.\r
1601  */\r
1602 Bignum biggcd(Bignum av, Bignum bv)\r
1603 {\r
1604     Bignum a = copybn(av);\r
1605     Bignum b = copybn(bv);\r
1606 \r
1607     while (bignum_cmp(b, Zero) != 0) {\r
1608         Bignum t = newbn(b[0]);\r
1609         bigdivmod(a, b, t, NULL);\r
1610         while (t[0] > 1 && t[t[0]] == 0)\r
1611             t[0]--;\r
1612         freebn(a);\r
1613         a = b;\r
1614         b = t;\r
1615     }\r
1616 \r
1617     freebn(b);\r
1618     return a;\r
1619 }\r
1620 \r
1621 /*\r
1622  * Modular inverse, using Euclid's extended algorithm.\r
1623  */\r
1624 Bignum modinv(Bignum number, Bignum modulus)\r
1625 {\r
1626     Bignum a = copybn(modulus);\r
1627     Bignum b = copybn(number);\r
1628     Bignum xp = copybn(Zero);\r
1629     Bignum x = copybn(One);\r
1630     int sign = +1;\r
1631 \r
1632     while (bignum_cmp(b, One) != 0) {\r
1633         Bignum t = newbn(b[0]);\r
1634         Bignum q = newbn(a[0]);\r
1635         bigdivmod(a, b, t, q);\r
1636         while (t[0] > 1 && t[t[0]] == 0)\r
1637             t[0]--;\r
1638         freebn(a);\r
1639         a = b;\r
1640         b = t;\r
1641         t = xp;\r
1642         xp = x;\r
1643         x = bigmuladd(q, xp, t);\r
1644         sign = -sign;\r
1645         freebn(t);\r
1646         freebn(q);\r
1647     }\r
1648 \r
1649     freebn(b);\r
1650     freebn(a);\r
1651     freebn(xp);\r
1652 \r
1653     /* now we know that sign * x == 1, and that x < modulus */\r
1654     if (sign < 0) {\r
1655         /* set a new x to be modulus - x */\r
1656         Bignum newx = newbn(modulus[0]);\r
1657         BignumInt carry = 0;\r
1658         int maxspot = 1;\r
1659         int i;\r
1660 \r
1661         for (i = 1; i <= (int)newx[0]; i++) {\r
1662             BignumInt aword = (i <= (int)modulus[0] ? modulus[i] : 0);\r
1663             BignumInt bword = (i <= (int)x[0] ? x[i] : 0);\r
1664             newx[i] = aword - bword - carry;\r
1665             bword = ~bword;\r
1666             carry = carry ? (newx[i] >= bword) : (newx[i] > bword);\r
1667             if (newx[i] != 0)\r
1668                 maxspot = i;\r
1669         }\r
1670         newx[0] = maxspot;\r
1671         freebn(x);\r
1672         x = newx;\r
1673     }\r
1674 \r
1675     /* and return. */\r
1676     return x;\r
1677 }\r
1678 \r
1679 /*\r
1680  * Render a bignum into decimal. Return a malloced string holding\r
1681  * the decimal representation.\r
1682  */\r
1683 char *bignum_decimal(Bignum x)\r
1684 {\r
1685     int ndigits, ndigit;\r
1686     int i, iszero;\r
1687     BignumDblInt carry;\r
1688     char *ret;\r
1689     BignumInt *workspace;\r
1690 \r
1691     /*\r
1692      * First, estimate the number of digits. Since log(10)/log(2)\r
1693      * is just greater than 93/28 (the joys of continued fraction\r
1694      * approximations...) we know that for every 93 bits, we need\r
1695      * at most 28 digits. This will tell us how much to malloc.\r
1696      *\r
1697      * Formally: if x has i bits, that means x is strictly less\r
1698      * than 2^i. Since 2 is less than 10^(28/93), this is less than\r
1699      * 10^(28i/93). We need an integer power of ten, so we must\r
1700      * round up (rounding down might make it less than x again).\r
1701      * Therefore if we multiply the bit count by 28/93, rounding\r
1702      * up, we will have enough digits.\r
1703      *\r
1704      * i=0 (i.e., x=0) is an irritating special case.\r
1705      */\r
1706     i = bignum_bitcount(x);\r
1707     if (!i)\r
1708         ndigits = 1;                   /* x = 0 */\r
1709     else\r
1710         ndigits = (28 * i + 92) / 93;  /* multiply by 28/93 and round up */\r
1711     ndigits++;                         /* allow for trailing \0 */\r
1712     ret = snewn(ndigits, char);\r
1713 \r
1714     /*\r
1715      * Now allocate some workspace to hold the binary form as we\r
1716      * repeatedly divide it by ten. Initialise this to the\r
1717      * big-endian form of the number.\r
1718      */\r
1719     workspace = snewn(x[0], BignumInt);\r
1720     for (i = 0; i < (int)x[0]; i++)\r
1721         workspace[i] = x[x[0] - i];\r
1722 \r
1723     /*\r
1724      * Next, write the decimal number starting with the last digit.\r
1725      * We use ordinary short division, dividing 10 into the\r
1726      * workspace.\r
1727      */\r
1728     ndigit = ndigits - 1;\r
1729     ret[ndigit] = '\0';\r
1730     do {\r
1731         iszero = 1;\r
1732         carry = 0;\r
1733         for (i = 0; i < (int)x[0]; i++) {\r
1734             carry = (carry << BIGNUM_INT_BITS) + workspace[i];\r
1735             workspace[i] = (BignumInt) (carry / 10);\r
1736             if (workspace[i])\r
1737                 iszero = 0;\r
1738             carry %= 10;\r
1739         }\r
1740         ret[--ndigit] = (char) (carry + '0');\r
1741     } while (!iszero);\r
1742 \r
1743     /*\r
1744      * There's a chance we've fallen short of the start of the\r
1745      * string. Correct if so.\r
1746      */\r
1747     if (ndigit > 0)\r
1748         memmove(ret, ret + ndigit, ndigits - ndigit);\r
1749 \r
1750     /*\r
1751      * Done.\r
1752      */\r
1753     sfree(workspace);\r
1754     return ret;\r
1755 }\r
1756 \r
1757 #ifdef TESTBN\r
1758 \r
1759 #include <stdio.h>\r
1760 #include <stdlib.h>\r
1761 #include <ctype.h>\r
1762 \r
1763 /*\r
1764  * gcc -g -O0 -DTESTBN -o testbn sshbn.c misc.c -I unix -I charset\r
1765  *\r
1766  * Then feed to this program's standard input the output of\r
1767  * testdata/bignum.py .\r
1768  */\r
1769 \r
1770 void modalfatalbox(char *p, ...)\r
1771 {\r
1772     va_list ap;\r
1773     fprintf(stderr, "FATAL ERROR: ");\r
1774     va_start(ap, p);\r
1775     vfprintf(stderr, p, ap);\r
1776     va_end(ap);\r
1777     fputc('\n', stderr);\r
1778     exit(1);\r
1779 }\r
1780 \r
1781 #define fromxdigit(c) ( (c)>'9' ? ((c)&0xDF) - 'A' + 10 : (c) - '0' )\r
1782 \r
1783 int main(int argc, char **argv)\r
1784 {\r
1785     char *buf;\r
1786     int line = 0;\r
1787     int passes = 0, fails = 0;\r
1788 \r
1789     while ((buf = fgetline(stdin)) != NULL) {\r
1790         int maxlen = strlen(buf);\r
1791         unsigned char *data = snewn(maxlen, unsigned char);\r
1792         unsigned char *ptrs[5], *q;\r
1793         int ptrnum;\r
1794         char *bufp = buf;\r
1795 \r
1796         line++;\r
1797 \r
1798         q = data;\r
1799         ptrnum = 0;\r
1800 \r
1801         while (*bufp && !isspace((unsigned char)*bufp))\r
1802             bufp++;\r
1803         if (bufp)\r
1804             *bufp++ = '\0';\r
1805 \r
1806         while (*bufp) {\r
1807             char *start, *end;\r
1808             int i;\r
1809 \r
1810             while (*bufp && !isxdigit((unsigned char)*bufp))\r
1811                 bufp++;\r
1812             start = bufp;\r
1813 \r
1814             if (!*bufp)\r
1815                 break;\r
1816 \r
1817             while (*bufp && isxdigit((unsigned char)*bufp))\r
1818                 bufp++;\r
1819             end = bufp;\r
1820 \r
1821             if (ptrnum >= lenof(ptrs))\r
1822                 break;\r
1823             ptrs[ptrnum++] = q;\r
1824             \r
1825             for (i = -((end - start) & 1); i < end-start; i += 2) {\r
1826                 unsigned char val = (i < 0 ? 0 : fromxdigit(start[i]));\r
1827                 val = val * 16 + fromxdigit(start[i+1]);\r
1828                 *q++ = val;\r
1829             }\r
1830 \r
1831             ptrs[ptrnum] = q;\r
1832         }\r
1833 \r
1834         if (!strcmp(buf, "mul")) {\r
1835             Bignum a, b, c, p;\r
1836 \r
1837             if (ptrnum != 3) {\r
1838                 printf("%d: mul with %d parameters, expected 3\n", line);\r
1839                 exit(1);\r
1840             }\r
1841             a = bignum_from_bytes(ptrs[0], ptrs[1]-ptrs[0]);\r
1842             b = bignum_from_bytes(ptrs[1], ptrs[2]-ptrs[1]);\r
1843             c = bignum_from_bytes(ptrs[2], ptrs[3]-ptrs[2]);\r
1844             p = bigmul(a, b);\r
1845 \r
1846             if (bignum_cmp(c, p) == 0) {\r
1847                 passes++;\r
1848             } else {\r
1849                 char *as = bignum_decimal(a);\r
1850                 char *bs = bignum_decimal(b);\r
1851                 char *cs = bignum_decimal(c);\r
1852                 char *ps = bignum_decimal(p);\r
1853                 \r
1854                 printf("%d: fail: %s * %s gave %s expected %s\n",\r
1855                        line, as, bs, ps, cs);\r
1856                 fails++;\r
1857 \r
1858                 sfree(as);\r
1859                 sfree(bs);\r
1860                 sfree(cs);\r
1861                 sfree(ps);\r
1862             }\r
1863             freebn(a);\r
1864             freebn(b);\r
1865             freebn(c);\r
1866             freebn(p);\r
1867         } else if (!strcmp(buf, "pow")) {\r
1868             Bignum base, expt, modulus, expected, answer;\r
1869 \r
1870             if (ptrnum != 4) {\r
1871                 printf("%d: mul with %d parameters, expected 3\n", line);\r
1872                 exit(1);\r
1873             }\r
1874 \r
1875             base = bignum_from_bytes(ptrs[0], ptrs[1]-ptrs[0]);\r
1876             expt = bignum_from_bytes(ptrs[1], ptrs[2]-ptrs[1]);\r
1877             modulus = bignum_from_bytes(ptrs[2], ptrs[3]-ptrs[2]);\r
1878             expected = bignum_from_bytes(ptrs[3], ptrs[4]-ptrs[3]);\r
1879             answer = modpow(base, expt, modulus);\r
1880 \r
1881             if (bignum_cmp(expected, answer) == 0) {\r
1882                 passes++;\r
1883             } else {\r
1884                 char *as = bignum_decimal(base);\r
1885                 char *bs = bignum_decimal(expt);\r
1886                 char *cs = bignum_decimal(modulus);\r
1887                 char *ds = bignum_decimal(answer);\r
1888                 char *ps = bignum_decimal(expected);\r
1889                 \r
1890                 printf("%d: fail: %s ^ %s mod %s gave %s expected %s\n",\r
1891                        line, as, bs, cs, ds, ps);\r
1892                 fails++;\r
1893 \r
1894                 sfree(as);\r
1895                 sfree(bs);\r
1896                 sfree(cs);\r
1897                 sfree(ds);\r
1898                 sfree(ps);\r
1899             }\r
1900             freebn(base);\r
1901             freebn(expt);\r
1902             freebn(modulus);\r
1903             freebn(expected);\r
1904             freebn(answer);\r
1905         } else {\r
1906             printf("%d: unrecognised test keyword: '%s'\n", line, buf);\r
1907             exit(1);\r
1908         }\r
1909 \r
1910         sfree(buf);\r
1911         sfree(data);\r
1912     }\r
1913 \r
1914     printf("passed %d failed %d total %d\n", passes, fails, passes+fails);\r
1915     return fails != 0;\r
1916 }\r
1917 \r
1918 #endif\r