sshbn.c

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