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 static int bn_clz(BignumInt x)
 536 {
 537     /*
 538      * Count the leading zero bits in x. Equivalently, how far left
 539      * would we need to shift x to make its top bit set?
 540      *
 541      * Precondition: x != 0.
 542      */
 543
 544     /* FIXME: would be nice to put in some compiler intrinsics under
 545      * ifdef here */
 546     int i, ret = 0;
 547     for (i = BIGNUM_INT_BITS / 2; i != 0; i >>= 1) {
 548         if ((x >> (BIGNUM_INT_BITS-i)) == 0) {
 549             x <<= i;
 550             ret += i;
 551         }
 552     }
 553     return ret;
 554 }
 555
 556 static BignumInt reciprocal_word(BignumInt d)
 557 {
 558     BignumInt dshort, recip;
 559     BignumDblInt product;
 560     int corrections;
 561
 562     /*
 563      * Input: a BignumInt value d, with its top bit set.
 564      */
 565     assert(d >> (BIGNUM_INT_BITS-1) == 1);
 566
 567     /*
 568      * Output: a value, shifted to fill a BignumInt, which is strictly
 569      * less than 1/(d+1), i.e. is an *under*-estimate (but by as
 570      * little as possible within the constraints) of the reciprocal of
 571      * any number whose first BIGNUM_INT_BITS bits match d.
 572      *
 573      * Ideally we'd like to _totally_ fill BignumInt, i.e. always
 574      * return a value with the top bit set. Unfortunately we can't
 575      * quite guarantee that for all inputs and also return a fixed
 576      * exponent. So instead we take our reciprocal to be
 577      * 2^(BIGNUM_INT_BITS*2-1) / d, so that it has the top bit clear
 578      * only in the exceptional case where d takes exactly the maximum
 579      * value BIGNUM_INT_MASK; in that case, the top bit is clear and
 580      * the next bit down is set.
 581      */
 582
 583     /*
 584      * Start by computing a half-length version of the answer, by
 585      * straightforward division within a BignumInt.
 586      */
 587     dshort = (d >> (BIGNUM_INT_BITS/2)) + 1;
 588     recip = (BIGNUM_TOP_BIT + dshort - 1) / dshort;
 589     recip <<= BIGNUM_INT_BITS - BIGNUM_INT_BITS/2;
 590
 591     /*
 592      * Newton-Raphson iteration to improve that starting reciprocal
 593      * estimate: take f(x) = d - 1/x, and then the N-R formula gives
 594      * x_new = x - f(x)/f'(x) = x - (d-1/x)/(1/x^2) = x(2-d*x). Or,
 595      * taking our fixed-point representation into account, take f(x)
 596      * to be d - K/x (where K = 2^(BIGNUM_INT_BITS*2-1) as discussed
 597      * above) and then we get (2K - d*x) * x/K.
 598      *
 599      * Newton-Raphson doubles the number of correct bits at every
 600      * iteration, and the initial division above already gave us half
 601      * the output word, so it's only worth doing one iteration.
 602      */
 603     product = MUL_WORD(recip, d);
 604     product += recip;
 605     product = -product;                /* the 2K shifts just off the top */
 606     product &= (((BignumDblInt)BIGNUM_INT_MASK << BIGNUM_INT_BITS) +
 607                 BIGNUM_INT_MASK);
 608     product >>= BIGNUM_INT_BITS;
 609     product = MUL_WORD(product, recip);
 610     product >>= (BIGNUM_INT_BITS-1);
 611     recip = (BignumInt)product;
 612
 613     /*
 614      * Now make sure we have the best possible reciprocal estimate,
 615      * before we return it. We might have been off by a handful either
 616      * way - not enough to bother with any better-thought-out kind of
 617      * correction loop.
 618      */
 619     product = MUL_WORD(recip, d);
 620     product += recip;
 621     corrections = 0;
 622     if (product >= ((BignumDblInt)1 << (2*BIGNUM_INT_BITS-1))) {
 623         do {
 624             product -= d;
 625             recip--;
 626             corrections++;
 627         } while (product >= ((BignumDblInt)1 << (2*BIGNUM_INT_BITS-1)));
 628     } else {
 629         while (product < ((BignumDblInt)1 << (2*BIGNUM_INT_BITS-1)) - d) {
 630             product += d;
 631             recip++;
 632             corrections++;
 633         }
 634     }
 635
 636     return recip;
 637 }
 638
 639 /*
 640  * Compute a = a % m.
 641  * Input in first alen words of a and first mlen words of m.
 642  * Output in first alen words of a
 643  * (of which first alen-mlen words will be zero).
 644  * Quotient is accumulated in the `quotient' array, which is a Bignum
 645  * rather than the internal bigendian format.
 646  *
 647  * 'recip' must be the result of calling reciprocal_word() on the top
 648  * BIGNUM_INT_BITS of the modulus (denoted m0 in comments below), with
 649  * the topmost set bit normalised to the MSB of the input to
 650  * reciprocal_word. 'rshift' is how far left the top nonzero word of
 651  * the modulus had to be shifted to set that top bit.
 652  */
 653 static void internal_mod(BignumInt *a, int alen,
 654                          BignumInt *m, int mlen,
 655                          BignumInt *quot, BignumInt recip, int rshift)
 656 {
 657     int i, k;
 658
 659 #ifdef DIVISION_DEBUG
 660     {
 661         int d;
 662         printf("start division, m=0x");
 663         for (d = 0; d < mlen; d++)
 664             printf("%0*llx", BIGNUM_INT_BITS/4, (unsigned long long)m[d]);
 665         printf(", recip=%#0*llx, rshift=%d\n",
 666                BIGNUM_INT_BITS/4, (unsigned long long)recip, rshift);
 667     }
 668 #endif
 669
 670     /*
 671      * Repeatedly use that reciprocal estimate to get a decent number
 672      * of quotient bits, and subtract off the resulting multiple of m.
 673      *
 674      * Normally we expect to terminate this loop by means of finding
 675      * out q=0 part way through, but one way in which we might not get
 676      * that far in the first place is if the input a is actually zero,
 677      * in which case we'll discard zero words from the front of a
 678      * until we reach the termination condition in the for statement
 679      * here.
 680      */
 681     for (i = 0; i <= alen - mlen ;) {
 682         BignumDblInt product, subtmp, t;
 683         BignumInt aword, q;
 684         int shift, full_bitoffset, bitoffset, wordoffset;
 685
 686 #ifdef DIVISION_DEBUG
 687         {
 688             int d;
 689             printf("main loop, a=0x");
 690             for (d = 0; d < alen; d++)
 691                 printf("%0*llx", BIGNUM_INT_BITS/4, (unsigned long long)a[d]);
 692             printf("\n");
 693         }
 694 #endif
 695
 696         if (a[i] == 0) {
 697 #ifdef DIVISION_DEBUG
 698             printf("zero word at i=%d\n", i);
 699 #endif
 700             i++;
 701             continue;
 702         }
 703
 704         aword = a[i];
 705         shift = bn_clz(aword);
 706         aword <<= shift;
 707         if (shift > 0 && i+1 < alen)
 708             aword |= a[i+1] >> (BIGNUM_INT_BITS - shift);
 709
 710         t = MUL_WORD(recip, aword);
 711         q = (BignumInt)(t >> BIGNUM_INT_BITS);
 712
 713 #ifdef DIVISION_DEBUG
 714         printf("i=%d, aword=%#0*llx, shift=%d, q=%#0*llx\n",
 715                i, BIGNUM_INT_BITS/4, (unsigned long long)aword,
 716                shift, BIGNUM_INT_BITS/4, (unsigned long long)q);
 717 #endif
 718
 719         /*
 720          * Work out the right bit and word offsets to use when
 721          * subtracting q*m from a.
 722          *
 723          * aword was taken from a[i], which means its LSB was at bit
 724          * position (alen-1-i) * BIGNUM_INT_BITS. But then we shifted
 725          * it left by 'shift', so now the low bit of aword corresponds
 726          * to bit position (alen-1-i) * BIGNUM_INT_BITS - shift, i.e.
 727          * aword is approximately equal to a / 2^(that).
 728          *
 729          * m0 comes from the top word of mod, so its LSB is at bit
 730          * position (mlen-1) * BIGNUM_INT_BITS - rshift, i.e. it can
 731          * be considered to be m / 2^(that power). 'recip' is the
 732          * reciprocal of m0, times 2^(BIGNUM_INT_BITS*2-1), i.e. it's
 733          * about 2^((mlen+1) * BIGNUM_INT_BITS - rshift - 1) / m.
 734          *
 735          * Hence, recip * aword is approximately equal to the product
 736          * of those, which simplifies to
 737          *
 738          * a/m * 2^((mlen+2+i-alen)*BIGNUM_INT_BITS + shift - rshift - 1)
 739          *
 740          * But we've also shifted recip*aword down by BIGNUM_INT_BITS
 741          * to form q, so we have
 742          *
 743          * q ~= a/m * 2^((mlen+1+i-alen)*BIGNUM_INT_BITS + shift - rshift - 1)
 744          *
 745          * and hence, when we now compute q*m, it will be about
 746          * a*2^(all that lot), i.e. the negation of that expression is
 747          * how far left we have to shift the product q*m to make it
 748          * approximately equal to a.
 749          */
 750         full_bitoffset = -((mlen+1+i-alen)*BIGNUM_INT_BITS + shift-rshift-1);
 751 #ifdef DIVISION_DEBUG
 752         printf("full_bitoffset=%d\n", full_bitoffset);
 753 #endif
 754
 755         if (full_bitoffset < 0) {
 756             /*
 757              * If we find ourselves needing to shift q*m _right_, that
 758              * means we've reached the bottom of the quotient. Clip q
 759              * so that its right shift becomes zero, and if that means
 760              * q becomes _actually_ zero, this loop is done.
 761              */
 762             if (full_bitoffset <= -BIGNUM_INT_BITS)
 763                 break;
 764             q >>= -full_bitoffset;
 765             full_bitoffset = 0;
 766             if (!q)
 767                 break;
 768 #ifdef DIVISION_DEBUG
 769             printf("now full_bitoffset=%d, q=%#0*llx\n",
 770                    full_bitoffset, BIGNUM_INT_BITS/4, (unsigned long long)q);
 771 #endif
 772         }
 773
 774         wordoffset = full_bitoffset / BIGNUM_INT_BITS;
 775         bitoffset = full_bitoffset % BIGNUM_INT_BITS;
 776 #ifdef DIVISION_DEBUG
 777         printf("wordoffset=%d, bitoffset=%d\n", wordoffset, bitoffset);
 778 #endif
 779
 780         /* wordoffset as computed above is the offset between the LSWs
 781          * of m and a. But in fact m and a are stored MSW-first, so we
 782          * need to adjust it to be the offset between the actual array
 783          * indices, and flip the sign too. */
 784         wordoffset = alen - mlen - wordoffset;
 785
 786         if (bitoffset == 0) {
 787             BignumInt c = 1;
 788             BignumInt prev_hi_word = 0;
 789             for (k = mlen - 1; wordoffset+k >= i; k--) {
 790                 BignumInt mword = k<0 ? 0 : m[k];
 791                 product = MUL_WORD(q, mword);
 792                 product += prev_hi_word;
 793                 prev_hi_word = product >> BIGNUM_INT_BITS;
 794 #ifdef DIVISION_DEBUG
 795                 printf("  aligned sub: product word for m[%d] = %#0*llx\n",
 796                        k, BIGNUM_INT_BITS/4,
 797                        (unsigned long long)(BignumInt)product);
 798 #endif
 799 #ifdef DIVISION_DEBUG
 800                 printf("  aligned sub: subtrahend for a[%d] = %#0*llx\n",
 801                        wordoffset+k, BIGNUM_INT_BITS/4,
 802                        (unsigned long long)(BignumInt)product);
 803 #endif
 804                 subtmp = (BignumDblInt)a[wordoffset+k] +
 805                     ((BignumInt)product ^ BIGNUM_INT_MASK) + c;
 806                 a[wordoffset+k] = (BignumInt)subtmp;
 807                 c = subtmp >> BIGNUM_INT_BITS;
 808             }
 809         } else {
 810             BignumInt add_word = 0;
 811             BignumInt c = 1;
 812             BignumInt prev_hi_word = 0;
 813             for (k = mlen - 1; wordoffset+k >= i; k--) {
 814                 BignumInt mword = k<0 ? 0 : m[k];
 815                 product = MUL_WORD(q, mword);
 816                 product += prev_hi_word;
 817                 prev_hi_word = product >> BIGNUM_INT_BITS;
 818 #ifdef DIVISION_DEBUG
 819                 printf("  unaligned sub: product word for m[%d] = %#0*llx\n",
 820                        k, BIGNUM_INT_BITS/4,
 821                        (unsigned long long)(BignumInt)product);
 822 #endif
 823
 824                 add_word |= (BignumInt)product << bitoffset;
 825
 826 #ifdef DIVISION_DEBUG
 827                 printf("  unaligned sub: subtrahend for a[%d] = %#0*llx\n",
 828                        wordoffset+k,
 829                        BIGNUM_INT_BITS/4, (unsigned long long)add_word);
 830 #endif
 831                 subtmp = (BignumDblInt)a[wordoffset+k] +
 832                     (add_word ^ BIGNUM_INT_MASK) + c;
 833                 a[wordoffset+k] = (BignumInt)subtmp;
 834                 c = subtmp >> BIGNUM_INT_BITS;
 835
 836                 add_word = (BignumInt)product >> (BIGNUM_INT_BITS - bitoffset);
 837             }
 838         }
 839
 840         if (quot) {
 841 #ifdef DIVISION_DEBUG
 842             printf("adding quotient word %#0*llx << %d\n",
 843                    BIGNUM_INT_BITS/4, (unsigned long long)q, full_bitoffset);
 844 #endif
 845             internal_add_shifted(quot, q, full_bitoffset);
 846 #ifdef DIVISION_DEBUG
 847             {
 848                 int d;
 849                 printf("now quot=0x");
 850                 for (d = quot[0]; d > 0; d--)
 851                     printf("%0*llx", BIGNUM_INT_BITS/4,
 852                            (unsigned long long)quot[d]);
 853                 printf("\n");
 854             }
 855 #endif
 856         }
 857     }
 858
 859 #ifdef DIVISION_DEBUG
 860     {
 861         int d;
 862         printf("end main loop, a=0x");
 863         for (d = 0; d < alen; d++)
 864             printf("%0*llx", BIGNUM_INT_BITS/4, (unsigned long long)a[d]);
 865         if (quot) {
 866             printf(", quot=0x");
 867             for (d = quot[0]; d > 0; d--)
 868                 printf("%0*llx", BIGNUM_INT_BITS/4,
 869                        (unsigned long long)quot[d]);
 870         }
 871         printf("\n");
 872     }
 873 #endif
 874
 875     /*
 876      * The above loop should terminate with the remaining value in a
 877      * being strictly less than 2*m (if a >= 2*m then we should always
 878      * have managed to get a nonzero q word), but we can't guarantee
 879      * that it will be strictly less than m: consider a case where the
 880      * remainder is 1, and another where the remainder is m-1. By the
 881      * time a contains a value that's _about m_, you clearly can't
 882      * distinguish those cases by looking at only the top word of a -
 883      * you have to go all the way down to the bottom before you find
 884      * out whether it's just less or just more than m.
 885      *
 886      * Hence, we now do a final fixup in which we subtract one last
 887      * copy of m, or don't, accordingly. We should never have to
 888      * subtract more than one copy of m here.
 889      */
 890     for (i = 0; i < alen; i++) {
 891         /* Compare a with m, word by word, from the MSW down. As soon
 892          * as we encounter a difference, we know whether we need the
 893          * fixup. */
 894         int mindex = mlen-alen+i;
 895         BignumInt mword = mindex < 0 ? 0 : m[mindex];
 896         if (a[i] < mword) {
 897 #ifdef DIVISION_DEBUG
 898             printf("final fixup not needed, a < m\n");
 899 #endif
 900             return;
 901         } else if (a[i] > mword) {
 902 #ifdef DIVISION_DEBUG
 903             printf("final fixup is needed, a > m\n");
 904 #endif
 905             break;
 906         }
 907         /* If neither of those cases happened, the words are the same,
 908          * so keep going and look at the next one. */
 909     }
 910 #ifdef DIVISION_DEBUG
 911     if (i == mlen) /* if we printed neither of the above diagnostics */
 912         printf("final fixup is needed, a == m\n");
 913 #endif
 914
 915     /*
 916      * If we got here without returning, then a >= m, so we must
 917      * subtract m, and increment the quotient.
 918      */
 919     {
 920         BignumInt c = 1;
 921         for (i = alen - 1; i >= 0; i--) {
 922             int mindex = mlen-alen+i;
 923             BignumInt mword = mindex < 0 ? 0 : m[mindex];
 924             BignumDblInt subtmp = (BignumDblInt)a[i] +
 925                 ((BignumInt)mword ^ BIGNUM_INT_MASK) + c;
 926             a[i] = (BignumInt)subtmp;
 927             c = subtmp >> BIGNUM_INT_BITS;
 928         }
 929     }
 930     if (quot)
 931         internal_add_shifted(quot, 1, 0);
 932
 933 #ifdef DIVISION_DEBUG
 934     {
 935         int d;
 936         printf("after final fixup, a=0x");
 937         for (d = 0; d < alen; d++)
 938             printf("%0*llx", BIGNUM_INT_BITS/4, (unsigned long long)a[d]);
 939         if (quot) {
 940             printf(", quot=0x");
 941             for (d = quot[0]; d > 0; d--)
 942                 printf("%0*llx", BIGNUM_INT_BITS/4,
 943                        (unsigned long long)quot[d]);
 944         }
 945         printf("\n");
 946     }
 947 #endif
 948 }
 949
 950 /*
 951  * Compute (base ^ exp) % mod, the pedestrian way.
 952  */
 953 Bignum modpow_simple(Bignum base_in, Bignum exp, Bignum mod)
 954 {
 955     BignumInt *a, *b, *n, *m, *scratch;
 956     BignumInt recip;
 957     int rshift;
 958     int mlen, scratchlen, i, j;
 959     Bignum base, result;
 960
 961     /*
 962      * The most significant word of mod needs to be non-zero. It
 963      * should already be, but let's make sure.
 964      */
 965     assert(mod[mod[0]] != 0);
 966
 967     /*
 968      * Make sure the base is smaller than the modulus, by reducing
 969      * it modulo the modulus if not.
 970      */
 971     base = bigmod(base_in, mod);
 972
 973     /* Allocate m of size mlen, copy mod to m */
 974     /* We use big endian internally */
 975     mlen = mod[0];
 976     m = snewn(mlen, BignumInt);
 977     for (j = 0; j < mlen; j++)
 978         m[j] = mod[mod[0] - j];
 979
 980     /* Allocate n of size mlen, copy base to n */
 981     n = snewn(mlen, BignumInt);
 982     i = mlen - base[0];
 983     for (j = 0; j < i; j++)
 984         n[j] = 0;
 985     for (j = 0; j < (int)base[0]; j++)
 986         n[i + j] = base[base[0] - j];
 987
 988     /* Allocate a and b of size 2*mlen. Set a = 1 */
 989     a = snewn(2 * mlen, BignumInt);
 990     b = snewn(2 * mlen, BignumInt);
 991     for (i = 0; i < 2 * mlen; i++)
 992         a[i] = 0;
 993     a[2 * mlen - 1] = 1;
 994
 995     /* Scratch space for multiplies */
 996     scratchlen = mul_compute_scratch(mlen);
 997     scratch = snewn(scratchlen, BignumInt);
 998
 999     /* Skip leading zero bits of exp. */
1000     i = 0;
1001     j = BIGNUM_INT_BITS-1;
1002     while (i < (int)exp[0] && (exp[exp[0] - i] & ((BignumInt)1 << j)) == 0) {
1003         j--;
1004         if (j < 0) {
1005             i++;
1006             j = BIGNUM_INT_BITS-1;
1007         }
1008     }
1009
1010     /* Compute reciprocal of the top full word of the modulus */
1011     {
1012         BignumInt m0 = m[0];
1013         rshift = bn_clz(m0);
1014         if (rshift) {
1015             m0 <<= rshift;
1016             if (mlen > 1)
1017                 m0 |= m[1] >> (BIGNUM_INT_BITS - rshift);
1018         }
1019         recip = reciprocal_word(m0);
1020     }
1021
1022     /* Main computation */
1023     while (i < (int)exp[0]) {
1024         while (j >= 0) {
1025             internal_mul(a + mlen, a + mlen, b, mlen, scratch);
1026             internal_mod(b, mlen * 2, m, mlen, NULL, recip, rshift);
1027             if ((exp[exp[0] - i] & ((BignumInt)1 << j)) != 0) {
1028                 internal_mul(b + mlen, n, a, mlen, scratch);
1029                 internal_mod(a, mlen * 2, m, mlen, NULL, recip, rshift);
1030             } else {
1031                 BignumInt *t;
1032                 t = a;
1033                 a = b;
1034                 b = t;
1035             }
1036             j--;
1037         }
1038         i++;
1039         j = BIGNUM_INT_BITS-1;
1040     }
1041
1042     /* Copy result to buffer */
1043     result = newbn(mod[0]);
1044     for (i = 0; i < mlen; i++)
1045         result[result[0] - i] = a[i + mlen];
1046     while (result[0] > 1 && result[result[0]] == 0)
1047         result[0]--;
1048
1049     /* Free temporary arrays */
1050     smemclr(a, 2 * mlen * sizeof(*a));
1051     sfree(a);
1052     smemclr(scratch, scratchlen * sizeof(*scratch));
1053     sfree(scratch);
1054     smemclr(b, 2 * mlen * sizeof(*b));
1055     sfree(b);
1056     smemclr(m, mlen * sizeof(*m));
1057     sfree(m);
1058     smemclr(n, mlen * sizeof(*n));
1059     sfree(n);
1060
1061     freebn(base);
1062
1063     return result;
1064 }
1065
1066 /*
1067  * Compute (base ^ exp) % mod. Uses the Montgomery multiplication
1068  * technique where possible, falling back to modpow_simple otherwise.
1069  */
1070 Bignum modpow(Bignum base_in, Bignum exp, Bignum mod)
1071 {
1072     BignumInt *a, *b, *x, *n, *mninv, *scratch;
1073     int len, scratchlen, i, j;
1074     Bignum base, base2, r, rn, inv, result;
1075
1076     /*
1077      * The most significant word of mod needs to be non-zero. It
1078      * should already be, but let's make sure.
1079      */
1080     assert(mod[mod[0]] != 0);
1081
1082     /*
1083      * mod had better be odd, or we can't do Montgomery multiplication
1084      * using a power of two at all.
1085      */
1086     if (!(mod[1] & 1))
1087         return modpow_simple(base_in, exp, mod);
1088
1089     /*
1090      * Make sure the base is smaller than the modulus, by reducing
1091      * it modulo the modulus if not.
1092      */
1093     base = bigmod(base_in, mod);
1094
1095     /*
1096      * Compute the inverse of n mod r, for monty_reduce. (In fact we
1097      * want the inverse of _minus_ n mod r, but we'll sort that out
1098      * below.)
1099      */
1100     len = mod[0];
1101     r = bn_power_2(BIGNUM_INT_BITS * len);
1102     inv = modinv(mod, r);
1103     assert(inv); /* cannot fail, since mod is odd and r is a power of 2 */
1104
1105     /*
1106      * Multiply the base by r mod n, to get it into Montgomery
1107      * representation.
1108      */
1109     base2 = modmul(base, r, mod);
1110     freebn(base);
1111     base = base2;
1112
1113     rn = bigmod(r, mod);               /* r mod n, i.e. Montgomerified 1 */
1114
1115     freebn(r);                         /* won't need this any more */
1116
1117     /*
1118      * Set up internal arrays of the right lengths, in big-endian
1119      * format, containing the base, the modulus, and the modulus's
1120      * inverse.
1121      */
1122     n = snewn(len, BignumInt);
1123     for (j = 0; j < len; j++)
1124         n[len - 1 - j] = mod[j + 1];
1125
1126     mninv = snewn(len, BignumInt);
1127     for (j = 0; j < len; j++)
1128         mninv[len - 1 - j] = (j < (int)inv[0] ? inv[j + 1] : 0);
1129     freebn(inv);         /* we don't need this copy of it any more */
1130     /* Now negate mninv mod r, so it's the inverse of -n rather than +n. */
1131     x = snewn(len, BignumInt);
1132     for (j = 0; j < len; j++)
1133         x[j] = 0;
1134     internal_sub(x, mninv, mninv, len);
1135
1136     /* x = snewn(len, BignumInt); */ /* already done above */
1137     for (j = 0; j < len; j++)
1138         x[len - 1 - j] = (j < (int)base[0] ? base[j + 1] : 0);
1139     freebn(base);        /* we don't need this copy of it any more */
1140
1141     a = snewn(2*len, BignumInt);
1142     b = snewn(2*len, BignumInt);
1143     for (j = 0; j < len; j++)
1144         a[2*len - 1 - j] = (j < (int)rn[0] ? rn[j + 1] : 0);
1145     freebn(rn);
1146
1147     /* Scratch space for multiplies */
1148     scratchlen = 3*len + mul_compute_scratch(len);
1149     scratch = snewn(scratchlen, BignumInt);
1150
1151     /* Skip leading zero bits of exp. */
1152     i = 0;
1153     j = BIGNUM_INT_BITS-1;
1154     while (i < (int)exp[0] && (exp[exp[0] - i] & ((BignumInt)1 << j)) == 0) {
1155         j--;
1156         if (j < 0) {
1157             i++;
1158             j = BIGNUM_INT_BITS-1;
1159         }
1160     }
1161
1162     /* Main computation */
1163     while (i < (int)exp[0]) {
1164         while (j >= 0) {
1165             internal_mul(a + len, a + len, b, len, scratch);
1166             monty_reduce(b, n, mninv, scratch, len);
1167             if ((exp[exp[0] - i] & ((BignumInt)1 << j)) != 0) {
1168                 internal_mul(b + len, x, a, len,  scratch);
1169                 monty_reduce(a, n, mninv, scratch, len);
1170             } else {
1171                 BignumInt *t;
1172                 t = a;
1173                 a = b;
1174                 b = t;
1175             }
1176             j--;
1177         }
1178         i++;
1179         j = BIGNUM_INT_BITS-1;
1180     }
1181
1182     /*
1183      * Final monty_reduce to get back from the adjusted Montgomery
1184      * representation.
1185      */
1186     monty_reduce(a, n, mninv, scratch, len);
1187
1188     /* Copy result to buffer */
1189     result = newbn(mod[0]);
1190     for (i = 0; i < len; i++)
1191         result[result[0] - i] = a[i + len];
1192     while (result[0] > 1 && result[result[0]] == 0)
1193         result[0]--;
1194
1195     /* Free temporary arrays */
1196     smemclr(scratch, scratchlen * sizeof(*scratch));
1197     sfree(scratch);
1198     smemclr(a, 2 * len * sizeof(*a));
1199     sfree(a);
1200     smemclr(b, 2 * len * sizeof(*b));
1201     sfree(b);
1202     smemclr(mninv, len * sizeof(*mninv));
1203     sfree(mninv);
1204     smemclr(n, len * sizeof(*n));
1205     sfree(n);
1206     smemclr(x, len * sizeof(*x));
1207     sfree(x);
1208
1209     return result;
1210 }
1211
1212 /*
1213  * Compute (p * q) % mod.
1214  * The most significant word of mod MUST be non-zero.
1215  * We assume that the result array is the same size as the mod array.
1216  */
1217 Bignum modmul(Bignum p, Bignum q, Bignum mod)
1218 {
1219     BignumInt *a, *n, *m, *o, *scratch;
1220     BignumInt recip;
1221     int rshift, scratchlen;
1222     int pqlen, mlen, rlen, i, j;
1223     Bignum result;
1224
1225     /*
1226      * The most significant word of mod needs to be non-zero. It
1227      * should already be, but let's make sure.
1228      */
1229     assert(mod[mod[0]] != 0);
1230
1231     /* Allocate m of size mlen, copy mod to m */
1232     /* We use big endian internally */
1233     mlen = mod[0];
1234     m = snewn(mlen, BignumInt);
1235     for (j = 0; j < mlen; j++)
1236         m[j] = mod[mod[0] - j];
1237
1238     pqlen = (p[0] > q[0] ? p[0] : q[0]);
1239
1240     /*
1241      * Make sure that we're allowing enough space. The shifting below
1242      * will underflow the vectors we allocate if pqlen is too small.
1243      */
1244     if (2*pqlen <= mlen)
1245         pqlen = mlen/2 + 1;
1246
1247     /* Allocate n of size pqlen, copy p to n */
1248     n = snewn(pqlen, BignumInt);
1249     i = pqlen - p[0];
1250     for (j = 0; j < i; j++)
1251         n[j] = 0;
1252     for (j = 0; j < (int)p[0]; j++)
1253         n[i + j] = p[p[0] - j];
1254
1255     /* Allocate o of size pqlen, copy q to o */
1256     o = snewn(pqlen, BignumInt);
1257     i = pqlen - q[0];
1258     for (j = 0; j < i; j++)
1259         o[j] = 0;
1260     for (j = 0; j < (int)q[0]; j++)
1261         o[i + j] = q[q[0] - j];
1262
1263     /* Allocate a of size 2*pqlen for result */
1264     a = snewn(2 * pqlen, BignumInt);
1265
1266     /* Scratch space for multiplies */
1267     scratchlen = mul_compute_scratch(pqlen);
1268     scratch = snewn(scratchlen, BignumInt);
1269
1270     /* Compute reciprocal of the top full word of the modulus */
1271     {
1272         BignumInt m0 = m[0];
1273         rshift = bn_clz(m0);
1274         if (rshift) {
1275             m0 <<= rshift;
1276             if (mlen > 1)
1277                 m0 |= m[1] >> (BIGNUM_INT_BITS - rshift);
1278         }
1279         recip = reciprocal_word(m0);
1280     }
1281
1282     /* Main computation */
1283     internal_mul(n, o, a, pqlen, scratch);
1284     internal_mod(a, pqlen * 2, m, mlen, NULL, recip, rshift);
1285
1286     /* Copy result to buffer */
1287     rlen = (mlen < pqlen * 2 ? mlen : pqlen * 2);
1288     result = newbn(rlen);
1289     for (i = 0; i < rlen; i++)
1290         result[result[0] - i] = a[i + 2 * pqlen - rlen];
1291     while (result[0] > 1 && result[result[0]] == 0)
1292         result[0]--;
1293
1294     /* Free temporary arrays */
1295     smemclr(scratch, scratchlen * sizeof(*scratch));
1296     sfree(scratch);
1297     smemclr(a, 2 * pqlen * sizeof(*a));
1298     sfree(a);
1299     smemclr(m, mlen * sizeof(*m));
1300     sfree(m);
1301     smemclr(n, pqlen * sizeof(*n));
1302     sfree(n);
1303     smemclr(o, pqlen * sizeof(*o));
1304     sfree(o);
1305
1306     return result;
1307 }
1308
1309 Bignum modsub(const Bignum a, const Bignum b, const Bignum n)
1310 {
1311     Bignum a1, b1, ret;
1312
1313     if (bignum_cmp(a, n) >= 0) a1 = bigmod(a, n);
1314     else a1 = a;
1315     if (bignum_cmp(b, n) >= 0) b1 = bigmod(b, n);
1316     else b1 = b;
1317
1318     if (bignum_cmp(a1, b1) >= 0) /* a >= b */
1319     {
1320         ret = bigsub(a1, b1);
1321     }
1322     else
1323     {
1324         /* Handle going round the corner of the modulus without having
1325          * negative support in Bignum */
1326         Bignum tmp = bigsub(n, b1);
1327         assert(tmp);
1328         ret = bigadd(tmp, a1);
1329         freebn(tmp);
1330     }
1331
1332     if (a != a1) freebn(a1);
1333     if (b != b1) freebn(b1);
1334
1335     return ret;
1336 }
1337
1338 /*
1339  * Compute p % mod.
1340  * The most significant word of mod MUST be non-zero.
1341  * We assume that the result array is the same size as the mod array.
1342  * We optionally write out a quotient if `quotient' is non-NULL.
1343  * We can avoid writing out the result if `result' is NULL.
1344  */
1345 static void bigdivmod(Bignum p, Bignum mod, Bignum result, Bignum quotient)
1346 {
1347     BignumInt *n, *m;
1348     BignumInt recip;
1349     int rshift;
1350     int plen, mlen, i, j;
1351
1352     /*
1353      * The most significant word of mod needs to be non-zero. It
1354      * should already be, but let's make sure.
1355      */
1356     assert(mod[mod[0]] != 0);
1357
1358     /* Allocate m of size mlen, copy mod to m */
1359     /* We use big endian internally */
1360     mlen = mod[0];
1361     m = snewn(mlen, BignumInt);
1362     for (j = 0; j < mlen; j++)
1363         m[j] = mod[mod[0] - j];
1364
1365     plen = p[0];
1366     /* Ensure plen > mlen */
1367     if (plen <= mlen)
1368         plen = mlen + 1;
1369
1370     /* Allocate n of size plen, copy p to n */
1371     n = snewn(plen, BignumInt);
1372     for (j = 0; j < plen; j++)
1373         n[j] = 0;
1374     for (j = 1; j <= (int)p[0]; j++)
1375         n[plen - j] = p[j];
1376
1377     /* Compute reciprocal of the top full word of the modulus */
1378     {
1379         BignumInt m0 = m[0];
1380         rshift = bn_clz(m0);
1381         if (rshift) {
1382             m0 <<= rshift;
1383             if (mlen > 1)
1384                 m0 |= m[1] >> (BIGNUM_INT_BITS - rshift);
1385         }
1386         recip = reciprocal_word(m0);
1387     }
1388
1389     /* Main computation */
1390     internal_mod(n, plen, m, mlen, quotient, recip, rshift);
1391
1392     /* Copy result to buffer */
1393     if (result) {
1394         for (i = 1; i <= (int)result[0]; i++) {
1395             int j = plen - i;
1396             result[i] = j >= 0 ? n[j] : 0;
1397         }
1398     }
1399
1400     /* Free temporary arrays */
1401     smemclr(m, mlen * sizeof(*m));
1402     sfree(m);
1403     smemclr(n, plen * sizeof(*n));
1404     sfree(n);
1405 }
1406
1407 /*
1408  * Decrement a number.
1409  */
1410 void decbn(Bignum bn)
1411 {
1412     int i = 1;
1413     while (i < (int)bn[0] && bn[i] == 0)
1414         bn[i++] = BIGNUM_INT_MASK;
1415     bn[i]--;
1416 }
1417
1418 Bignum bignum_from_bytes(const unsigned char *data, int nbytes)
1419 {
1420     Bignum result;
1421     int w, i;
1422
1423     assert(nbytes >= 0 && nbytes < INT_MAX/8);
1424
1425     w = (nbytes + BIGNUM_INT_BYTES - 1) / BIGNUM_INT_BYTES; /* bytes->words */
1426
1427     result = newbn(w);
1428     for (i = 1; i <= w; i++)
1429         result[i] = 0;
1430     for (i = nbytes; i--;) {
1431         unsigned char byte = *data++;
1432         result[1 + i / BIGNUM_INT_BYTES] |=
1433             (BignumInt)byte << (8*i % BIGNUM_INT_BITS);
1434     }
1435
1436     bn_restore_invariant(result);
1437     return result;
1438 }
1439
1440 Bignum bignum_from_bytes_le(const unsigned char *data, int nbytes)
1441 {
1442     Bignum result;
1443     int w, i;
1444
1445     assert(nbytes >= 0 && nbytes < INT_MAX/8);
1446
1447     w = (nbytes + BIGNUM_INT_BYTES - 1) / BIGNUM_INT_BYTES; /* bytes->words */
1448
1449     result = newbn(w);
1450     for (i = 1; i <= w; i++)
1451         result[i] = 0;
1452     for (i = 0; i < nbytes; ++i) {
1453         unsigned char byte = *data++;
1454         result[1 + i / BIGNUM_INT_BYTES] |=
1455             (BignumInt)byte << (8*i % BIGNUM_INT_BITS);
1456     }
1457
1458     bn_restore_invariant(result);
1459     return result;
1460 }
1461
1462 Bignum bignum_from_decimal(const char *decimal)
1463 {
1464     Bignum result = copybn(Zero);
1465
1466     while (*decimal) {
1467         Bignum tmp, tmp2;
1468
1469         if (!isdigit((unsigned char)*decimal)) {
1470             freebn(result);
1471             return 0;
1472         }
1473
1474         tmp = bigmul(result, Ten);
1475         tmp2 = bignum_from_long(*decimal - '0');
1476         result = bigadd(tmp, tmp2);
1477         freebn(tmp);
1478         freebn(tmp2);
1479
1480         decimal++;
1481     }
1482
1483     return result;
1484 }
1485
1486 Bignum bignum_random_in_range(const Bignum lower, const Bignum upper)
1487 {
1488     Bignum ret = NULL;
1489     unsigned char *bytes;
1490     int upper_len = bignum_bitcount(upper);
1491     int upper_bytes = upper_len / 8;
1492     int upper_bits = upper_len % 8;
1493     if (upper_bits) ++upper_bytes;
1494
1495     bytes = snewn(upper_bytes, unsigned char);
1496     do {
1497         int i;
1498
1499         if (ret) freebn(ret);
1500
1501         for (i = 0; i < upper_bytes; ++i)
1502         {
1503             bytes[i] = (unsigned char)random_byte();
1504         }
1505         /* Mask the top to reduce failure rate to 50/50 */
1506         if (upper_bits)
1507         {
1508             bytes[i - 1] &= 0xFF >> (8 - upper_bits);
1509         }
1510
1511         ret = bignum_from_bytes(bytes, upper_bytes);
1512     } while (bignum_cmp(ret, lower) < 0 || bignum_cmp(ret, upper) > 0);
1513     smemclr(bytes, upper_bytes);
1514     sfree(bytes);
1515
1516     return ret;
1517 }
1518
1519 /*
1520  * Read an SSH-1-format bignum from a data buffer. Return the number
1521  * of bytes consumed, or -1 if there wasn't enough data.
1522  */
1523 int ssh1_read_bignum(const unsigned char *data, int len, Bignum * result)
1524 {
1525     const unsigned char *p = data;
1526     int i;
1527     int w, b;
1528
1529     if (len < 2)
1530         return -1;
1531
1532     w = 0;
1533     for (i = 0; i < 2; i++)
1534         w = (w << 8) + *p++;
1535     b = (w + 7) / 8;                   /* bits -> bytes */
1536
1537     if (len < b+2)
1538         return -1;
1539
1540     if (!result)                       /* just return length */
1541         return b + 2;
1542
1543     *result = bignum_from_bytes(p, b);
1544
1545     return p + b - data;
1546 }
1547
1548 /*
1549  * Return the bit count of a bignum, for SSH-1 encoding.
1550  */
1551 int bignum_bitcount(Bignum bn)
1552 {
1553     int bitcount = bn[0] * BIGNUM_INT_BITS - 1;
1554     while (bitcount >= 0
1555            && (bn[bitcount / BIGNUM_INT_BITS + 1] >> (bitcount % BIGNUM_INT_BITS)) == 0) bitcount--;
1556     return bitcount + 1;
1557 }
1558
1559 /*
1560  * Return the byte length of a bignum when SSH-1 encoded.
1561  */
1562 int ssh1_bignum_length(Bignum bn)
1563 {
1564     return 2 + (bignum_bitcount(bn) + 7) / 8;
1565 }
1566
1567 /*
1568  * Return the byte length of a bignum when SSH-2 encoded.
1569  */
1570 int ssh2_bignum_length(Bignum bn)
1571 {
1572     return 4 + (bignum_bitcount(bn) + 8) / 8;
1573 }
1574
1575 /*
1576  * Return a byte from a bignum; 0 is least significant, etc.
1577  */
1578 int bignum_byte(Bignum bn, int i)
1579 {
1580     if (i < 0 || i >= (int)(BIGNUM_INT_BYTES * bn[0]))
1581         return 0;                      /* beyond the end */
1582     else
1583         return (bn[i / BIGNUM_INT_BYTES + 1] >>
1584                 ((i % BIGNUM_INT_BYTES)*8)) & 0xFF;
1585 }
1586
1587 /*
1588  * Return a bit from a bignum; 0 is least significant, etc.
1589  */
1590 int bignum_bit(Bignum bn, int i)
1591 {
1592     if (i < 0 || i >= (int)(BIGNUM_INT_BITS * bn[0]))
1593         return 0;                      /* beyond the end */
1594     else
1595         return (bn[i / BIGNUM_INT_BITS + 1] >> (i % BIGNUM_INT_BITS)) & 1;
1596 }
1597
1598 /*
1599  * Set a bit in a bignum; 0 is least significant, etc.
1600  */
1601 void bignum_set_bit(Bignum bn, int bitnum, int value)
1602 {
1603     if (bitnum < 0 || bitnum >= (int)(BIGNUM_INT_BITS * bn[0])) {
1604         if (value) abort();                    /* beyond the end */
1605     } else {
1606         int v = bitnum / BIGNUM_INT_BITS + 1;
1607         BignumInt mask = (BignumInt)1 << (bitnum % BIGNUM_INT_BITS);
1608         if (value)
1609             bn[v] |= mask;
1610         else
1611             bn[v] &= ~mask;
1612     }
1613 }
1614
1615 /*
1616  * Write a SSH-1-format bignum into a buffer. It is assumed the
1617  * buffer is big enough. Returns the number of bytes used.
1618  */
1619 int ssh1_write_bignum(void *data, Bignum bn)
1620 {
1621     unsigned char *p = data;
1622     int len = ssh1_bignum_length(bn);
1623     int i;
1624     int bitc = bignum_bitcount(bn);
1625
1626     *p++ = (bitc >> 8) & 0xFF;
1627     *p++ = (bitc) & 0xFF;
1628     for (i = len - 2; i--;)
1629         *p++ = bignum_byte(bn, i);
1630     return len;
1631 }
1632
1633 /*
1634  * Compare two bignums. Returns like strcmp.
1635  */
1636 int bignum_cmp(Bignum a, Bignum b)
1637 {
1638     int amax = a[0], bmax = b[0];
1639     int i;
1640
1641     /* Annoyingly we have two representations of zero */
1642     if (amax == 1 && a[amax] == 0)
1643         amax = 0;
1644     if (bmax == 1 && b[bmax] == 0)
1645         bmax = 0;
1646
1647     assert(amax == 0 || a[amax] != 0);
1648     assert(bmax == 0 || b[bmax] != 0);
1649
1650     i = (amax > bmax ? amax : bmax);
1651     while (i) {
1652         BignumInt aval = (i > amax ? 0 : a[i]);
1653         BignumInt bval = (i > bmax ? 0 : b[i]);
1654         if (aval < bval)
1655             return -1;
1656         if (aval > bval)
1657             return +1;
1658         i--;
1659     }
1660     return 0;
1661 }
1662
1663 /*
1664  * Right-shift one bignum to form another.
1665  */
1666 Bignum bignum_rshift(Bignum a, int shift)
1667 {
1668     Bignum ret;
1669     int i, shiftw, shiftb, shiftbb, bits;
1670     BignumInt ai, ai1;
1671
1672     assert(shift >= 0);
1673
1674     bits = bignum_bitcount(a) - shift;
1675     ret = newbn((bits + BIGNUM_INT_BITS - 1) / BIGNUM_INT_BITS);
1676
1677     if (ret) {
1678         shiftw = shift / BIGNUM_INT_BITS;
1679         shiftb = shift % BIGNUM_INT_BITS;
1680         shiftbb = BIGNUM_INT_BITS - shiftb;
1681
1682         ai1 = a[shiftw + 1];
1683         for (i = 1; i <= (int)ret[0]; i++) {
1684             ai = ai1;
1685             ai1 = (i + shiftw + 1 <= (int)a[0] ? a[i + shiftw + 1] : 0);
1686             ret[i] = ((ai >> shiftb) | (ai1 << shiftbb)) & BIGNUM_INT_MASK;
1687         }
1688     }
1689
1690     return ret;
1691 }
1692
1693 /*
1694  * Left-shift one bignum to form another.
1695  */
1696 Bignum bignum_lshift(Bignum a, int shift)
1697 {
1698     Bignum ret;
1699     int bits, shiftWords, shiftBits;
1700
1701     assert(shift >= 0);
1702
1703     bits = bignum_bitcount(a) + shift;
1704     ret = newbn((bits + BIGNUM_INT_BITS - 1) / BIGNUM_INT_BITS);
1705
1706     shiftWords = shift / BIGNUM_INT_BITS;
1707     shiftBits = shift % BIGNUM_INT_BITS;
1708
1709     if (shiftBits == 0)
1710     {
1711         memcpy(&ret[1 + shiftWords], &a[1], sizeof(BignumInt) * a[0]);
1712     }
1713     else
1714     {
1715         int i;
1716         BignumInt carry = 0;
1717
1718         /* Remember that Bignum[0] is length, so add 1 */
1719         for (i = shiftWords + 1; i < ((int)a[0]) + shiftWords + 1; ++i)
1720         {
1721             BignumInt from = a[i - shiftWords];
1722             ret[i] = (from << shiftBits) | carry;
1723             carry = from >> (BIGNUM_INT_BITS - shiftBits);
1724         }
1725         if (carry) ret[i] = carry;
1726     }
1727
1728     return ret;
1729 }
1730
1731 /*
1732  * Non-modular multiplication and addition.
1733  */
1734 Bignum bigmuladd(Bignum a, Bignum b, Bignum addend)
1735 {
1736     int alen = a[0], blen = b[0];
1737     int mlen = (alen > blen ? alen : blen);
1738     int rlen, i, maxspot;
1739     int wslen;
1740     BignumInt *workspace;
1741     Bignum ret;
1742
1743     /* mlen space for a, mlen space for b, 2*mlen for result,
1744      * plus scratch space for multiplication */
1745     wslen = mlen * 4 + mul_compute_scratch(mlen);
1746     workspace = snewn(wslen, BignumInt);
1747     for (i = 0; i < mlen; i++) {
1748         workspace[0 * mlen + i] = (mlen - i <= (int)a[0] ? a[mlen - i] : 0);
1749         workspace[1 * mlen + i] = (mlen - i <= (int)b[0] ? b[mlen - i] : 0);
1750     }
1751
1752     internal_mul(workspace + 0 * mlen, workspace + 1 * mlen,
1753                  workspace + 2 * mlen, mlen, workspace + 4 * mlen);
1754
1755     /* now just copy the result back */
1756     rlen = alen + blen + 1;
1757     if (addend && rlen <= (int)addend[0])
1758         rlen = addend[0] + 1;
1759     ret = newbn(rlen);
1760     maxspot = 0;
1761     for (i = 1; i <= (int)ret[0]; i++) {
1762         ret[i] = (i <= 2 * mlen ? workspace[4 * mlen - i] : 0);
1763         if (ret[i] != 0)
1764             maxspot = i;
1765     }
1766     ret[0] = maxspot;
1767
1768     /* now add in the addend, if any */
1769     if (addend) {
1770         BignumDblInt carry = 0;
1771         for (i = 1; i <= rlen; i++) {
1772             carry += (i <= (int)ret[0] ? ret[i] : 0);
1773             carry += (i <= (int)addend[0] ? addend[i] : 0);
1774             ret[i] = (BignumInt) carry & BIGNUM_INT_MASK;
1775             carry >>= BIGNUM_INT_BITS;
1776             if (ret[i] != 0 && i > maxspot)
1777                 maxspot = i;
1778         }
1779     }
1780     ret[0] = maxspot;
1781
1782     smemclr(workspace, wslen * sizeof(*workspace));
1783     sfree(workspace);
1784     return ret;
1785 }
1786
1787 /*
1788  * Non-modular multiplication.
1789  */
1790 Bignum bigmul(Bignum a, Bignum b)
1791 {
1792     return bigmuladd(a, b, NULL);
1793 }
1794
1795 /*
1796  * Simple addition.
1797  */
1798 Bignum bigadd(Bignum a, Bignum b)
1799 {
1800     int alen = a[0], blen = b[0];
1801     int rlen = (alen > blen ? alen : blen) + 1;
1802     int i, maxspot;
1803     Bignum ret;
1804     BignumDblInt carry;
1805
1806     ret = newbn(rlen);
1807
1808     carry = 0;
1809     maxspot = 0;
1810     for (i = 1; i <= rlen; i++) {
1811         carry += (i <= (int)a[0] ? a[i] : 0);
1812         carry += (i <= (int)b[0] ? b[i] : 0);
1813         ret[i] = (BignumInt) carry & BIGNUM_INT_MASK;
1814         carry >>= BIGNUM_INT_BITS;
1815         if (ret[i] != 0 && i > maxspot)
1816             maxspot = i;
1817     }
1818     ret[0] = maxspot;
1819
1820     return ret;
1821 }
1822
1823 /*
1824  * Subtraction. Returns a-b, or NULL if the result would come out
1825  * negative (recall that this entire bignum module only handles
1826  * positive numbers).
1827  */
1828 Bignum bigsub(Bignum a, Bignum b)
1829 {
1830     int alen = a[0], blen = b[0];
1831     int rlen = (alen > blen ? alen : blen);
1832     int i, maxspot;
1833     Bignum ret;
1834     BignumDblInt carry;
1835
1836     ret = newbn(rlen);
1837
1838     carry = 1;
1839     maxspot = 0;
1840     for (i = 1; i <= rlen; i++) {
1841         carry += (i <= (int)a[0] ? a[i] : 0);
1842         carry += (i <= (int)b[0] ? b[i] ^ BIGNUM_INT_MASK : BIGNUM_INT_MASK);
1843         ret[i] = (BignumInt) carry & BIGNUM_INT_MASK;
1844         carry >>= BIGNUM_INT_BITS;
1845         if (ret[i] != 0 && i > maxspot)
1846             maxspot = i;
1847     }
1848     ret[0] = maxspot;
1849
1850     if (!carry) {
1851         freebn(ret);
1852         return NULL;
1853     }
1854
1855     return ret;
1856 }
1857
1858 /*
1859  * Create a bignum which is the bitmask covering another one. That
1860  * is, the smallest integer which is >= N and is also one less than
1861  * a power of two.
1862  */
1863 Bignum bignum_bitmask(Bignum n)
1864 {
1865     Bignum ret = copybn(n);
1866     int i;
1867     BignumInt j;
1868
1869     i = ret[0];
1870     while (n[i] == 0 && i > 0)
1871         i--;
1872     if (i <= 0)
1873         return ret;                    /* input was zero */
1874     j = 1;
1875     while (j < n[i])
1876         j = 2 * j + 1;
1877     ret[i] = j;
1878     while (--i > 0)
1879         ret[i] = BIGNUM_INT_MASK;
1880     return ret;
1881 }
1882
1883 /*
1884  * Convert a (max 32-bit) long into a bignum.
1885  */
1886 Bignum bignum_from_long(unsigned long nn)
1887 {
1888     Bignum ret;
1889     BignumDblInt n = nn;
1890
1891     ret = newbn(3);
1892     ret[1] = (BignumInt)(n & BIGNUM_INT_MASK);
1893     ret[2] = (BignumInt)((n >> BIGNUM_INT_BITS) & BIGNUM_INT_MASK);
1894     ret[3] = 0;
1895     ret[0] = (ret[2]  ? 2 : 1);
1896     return ret;
1897 }
1898
1899 /*
1900  * Add a long to a bignum.
1901  */
1902 Bignum bignum_add_long(Bignum number, unsigned long addendx)
1903 {
1904     Bignum ret = newbn(number[0] + 1);
1905     int i, maxspot = 0;
1906     BignumDblInt carry = 0, addend = addendx;
1907
1908     for (i = 1; i <= (int)ret[0]; i++) {
1909         carry += addend & BIGNUM_INT_MASK;
1910         carry += (i <= (int)number[0] ? number[i] : 0);
1911         addend >>= BIGNUM_INT_BITS;
1912         ret[i] = (BignumInt) carry & BIGNUM_INT_MASK;
1913         carry >>= BIGNUM_INT_BITS;
1914         if (ret[i] != 0)
1915             maxspot = i;
1916     }
1917     ret[0] = maxspot;
1918     return ret;
1919 }
1920
1921 /*
1922  * Compute the residue of a bignum, modulo a (max 16-bit) short.
1923  */
1924 unsigned short bignum_mod_short(Bignum number, unsigned short modulus)
1925 {
1926     BignumDblInt mod, r;
1927     int i;
1928
1929     r = 0;
1930     mod = modulus;
1931     for (i = number[0]; i > 0; i--)
1932         r = (r * (BIGNUM_TOP_BIT % mod) * 2 + number[i] % mod) % mod;
1933     return (unsigned short) r;
1934 }
1935
1936 #ifdef DEBUG
1937 void diagbn(char *prefix, Bignum md)
1938 {
1939     int i, nibbles, morenibbles;
1940     static const char hex[] = "0123456789ABCDEF";
1941
1942     debug(("%s0x", prefix ? prefix : ""));
1943
1944     nibbles = (3 + bignum_bitcount(md)) / 4;
1945     if (nibbles < 1)
1946         nibbles = 1;
1947     morenibbles = 4 * md[0] - nibbles;
1948     for (i = 0; i < morenibbles; i++)
1949         debug(("-"));
1950     for (i = nibbles; i--;)
1951         debug(("%c",
1952                hex[(bignum_byte(md, i / 2) >> (4 * (i % 2))) & 0xF]));
1953
1954     if (prefix)
1955         debug(("\n"));
1956 }
1957 #endif
1958
1959 /*
1960  * Simple division.
1961  */
1962 Bignum bigdiv(Bignum a, Bignum b)
1963 {
1964     Bignum q = newbn(a[0]);
1965     bigdivmod(a, b, NULL, q);
1966     while (q[0] > 1 && q[q[0]] == 0)
1967         q[0]--;
1968     return q;
1969 }
1970
1971 /*
1972  * Simple remainder.
1973  */
1974 Bignum bigmod(Bignum a, Bignum b)
1975 {
1976     Bignum r = newbn(b[0]);
1977     bigdivmod(a, b, r, NULL);
1978     while (r[0] > 1 && r[r[0]] == 0)
1979         r[0]--;
1980     return r;
1981 }
1982
1983 /*
1984  * Greatest common divisor.
1985  */
1986 Bignum biggcd(Bignum av, Bignum bv)
1987 {
1988     Bignum a = copybn(av);
1989     Bignum b = copybn(bv);
1990
1991     while (bignum_cmp(b, Zero) != 0) {
1992         Bignum t = newbn(b[0]);
1993         bigdivmod(a, b, t, NULL);
1994         while (t[0] > 1 && t[t[0]] == 0)
1995             t[0]--;
1996         freebn(a);
1997         a = b;
1998         b = t;
1999     }
2000
2001     freebn(b);
2002     return a;
2003 }
2004
2005 /*
2006  * Modular inverse, using Euclid's extended algorithm.
2007  */
2008 Bignum modinv(Bignum number, Bignum modulus)
2009 {
2010     Bignum a = copybn(modulus);
2011     Bignum b = copybn(number);
2012     Bignum xp = copybn(Zero);
2013     Bignum x = copybn(One);
2014     int sign = +1;
2015
2016     assert(number[number[0]] != 0);
2017     assert(modulus[modulus[0]] != 0);
2018
2019     while (bignum_cmp(b, One) != 0) {
2020         Bignum t, q;
2021
2022         if (bignum_cmp(b, Zero) == 0) {
2023             /*
2024              * Found a common factor between the inputs, so we cannot
2025              * return a modular inverse at all.
2026              */
2027             freebn(b);
2028             freebn(a);
2029             freebn(xp);
2030             freebn(x);
2031             return NULL;
2032         }
2033
2034         t = newbn(b[0]);
2035         q = newbn(a[0]);
2036         bigdivmod(a, b, t, q);
2037         while (t[0] > 1 && t[t[0]] == 0)
2038             t[0]--;
2039         while (q[0] > 1 && q[q[0]] == 0)
2040             q[0]--;
2041         freebn(a);
2042         a = b;
2043         b = t;
2044         t = xp;
2045         xp = x;
2046         x = bigmuladd(q, xp, t);
2047         sign = -sign;
2048         freebn(t);
2049         freebn(q);
2050     }
2051
2052     freebn(b);
2053     freebn(a);
2054     freebn(xp);
2055
2056     /* now we know that sign * x == 1, and that x < modulus */
2057     if (sign < 0) {
2058         /* set a new x to be modulus - x */
2059         Bignum newx = newbn(modulus[0]);
2060         BignumInt carry = 0;
2061         int maxspot = 1;
2062         int i;
2063
2064         for (i = 1; i <= (int)newx[0]; i++) {
2065             BignumInt aword = (i <= (int)modulus[0] ? modulus[i] : 0);
2066             BignumInt bword = (i <= (int)x[0] ? x[i] : 0);
2067             newx[i] = aword - bword - carry;
2068             bword = ~bword;
2069             carry = carry ? (newx[i] >= bword) : (newx[i] > bword);
2070             if (newx[i] != 0)
2071                 maxspot = i;
2072         }
2073         newx[0] = maxspot;
2074         freebn(x);
2075         x = newx;
2076     }
2077
2078     /* and return. */
2079     return x;
2080 }
2081
2082 /*
2083  * Render a bignum into decimal. Return a malloced string holding
2084  * the decimal representation.
2085  */
2086 char *bignum_decimal(Bignum x)
2087 {
2088     int ndigits, ndigit;
2089     int i, iszero;
2090     BignumDblInt carry;
2091     char *ret;
2092     BignumInt *workspace;
2093
2094     /*
2095      * First, estimate the number of digits. Since log(10)/log(2)
2096      * is just greater than 93/28 (the joys of continued fraction
2097      * approximations...) we know that for every 93 bits, we need
2098      * at most 28 digits. This will tell us how much to malloc.
2099      *
2100      * Formally: if x has i bits, that means x is strictly less
2101      * than 2^i. Since 2 is less than 10^(28/93), this is less than
2102      * 10^(28i/93). We need an integer power of ten, so we must
2103      * round up (rounding down might make it less than x again).
2104      * Therefore if we multiply the bit count by 28/93, rounding
2105      * up, we will have enough digits.
2106      *
2107      * i=0 (i.e., x=0) is an irritating special case.
2108      */
2109     i = bignum_bitcount(x);
2110     if (!i)
2111         ndigits = 1;                   /* x = 0 */
2112     else
2113         ndigits = (28 * i + 92) / 93;  /* multiply by 28/93 and round up */
2114     ndigits++;                         /* allow for trailing \0 */
2115     ret = snewn(ndigits, char);
2116
2117     /*
2118      * Now allocate some workspace to hold the binary form as we
2119      * repeatedly divide it by ten. Initialise this to the
2120      * big-endian form of the number.
2121      */
2122     workspace = snewn(x[0], BignumInt);
2123     for (i = 0; i < (int)x[0]; i++)
2124         workspace[i] = x[x[0] - i];
2125
2126     /*
2127      * Next, write the decimal number starting with the last digit.
2128      * We use ordinary short division, dividing 10 into the
2129      * workspace.
2130      */
2131     ndigit = ndigits - 1;
2132     ret[ndigit] = '\0';
2133     do {
2134         iszero = 1;
2135         carry = 0;
2136         for (i = 0; i < (int)x[0]; i++) {
2137             carry = (carry << BIGNUM_INT_BITS) + workspace[i];
2138             workspace[i] = (BignumInt) (carry / 10);
2139             if (workspace[i])
2140                 iszero = 0;
2141             carry %= 10;
2142         }
2143         ret[--ndigit] = (char) (carry + '0');
2144     } while (!iszero);
2145
2146     /*
2147      * There's a chance we've fallen short of the start of the
2148      * string. Correct if so.
2149      */
2150     if (ndigit > 0)
2151         memmove(ret, ret + ndigit, ndigits - ndigit);
2152
2153     /*
2154      * Done.
2155      */
2156     smemclr(workspace, x[0] * sizeof(*workspace));
2157     sfree(workspace);
2158     return ret;
2159 }
2160
2161 #ifdef TESTBN
2162
2163 #include <stdio.h>
2164 #include <stdlib.h>
2165 #include <ctype.h>
2166
2167 /*
2168  * gcc -Wall -g -O0 -DTESTBN -o testbn sshbn.c misc.c conf.c tree234.c unix/uxmisc.c -I. -I unix -I charset
2169  *
2170  * Then feed to this program's standard input the output of
2171  * testdata/bignum.py .
2172  */
2173
2174 void modalfatalbox(const char *p, ...)
2175 {
2176     va_list ap;
2177     fprintf(stderr, "FATAL ERROR: ");
2178     va_start(ap, p);
2179     vfprintf(stderr, p, ap);
2180     va_end(ap);
2181     fputc('\n', stderr);
2182     exit(1);
2183 }
2184
2185 int random_byte(void)
2186 {
2187     modalfatalbox("random_byte called in testbn");
2188     return 0;
2189 }
2190
2191 #define fromxdigit(c) ( (c)>'9' ? ((c)&0xDF) - 'A' + 10 : (c) - '0' )
2192
2193 int main(int argc, char **argv)
2194 {
2195     char *buf;
2196     int line = 0;
2197     int passes = 0, fails = 0;
2198
2199     while ((buf = fgetline(stdin)) != NULL) {
2200         int maxlen = strlen(buf);
2201         unsigned char *data = snewn(maxlen, unsigned char);
2202         unsigned char *ptrs[5], *q;
2203         int ptrnum;
2204         char *bufp = buf;
2205
2206         line++;
2207
2208         q = data;
2209         ptrnum = 0;
2210
2211         while (*bufp && !isspace((unsigned char)*bufp))
2212             bufp++;
2213         if (bufp)
2214             *bufp++ = '\0';
2215
2216         while (*bufp) {
2217             char *start, *end;
2218             int i;
2219
2220             while (*bufp && !isxdigit((unsigned char)*bufp))
2221                 bufp++;
2222             start = bufp;
2223
2224             if (!*bufp)
2225                 break;
2226
2227             while (*bufp && isxdigit((unsigned char)*bufp))
2228                 bufp++;
2229             end = bufp;
2230
2231             if (ptrnum >= lenof(ptrs))
2232                 break;
2233             ptrs[ptrnum++] = q;
2234
2235             for (i = -((end - start) & 1); i < end-start; i += 2) {
2236                 unsigned char val = (i < 0 ? 0 : fromxdigit(start[i]));
2237                 val = val * 16 + fromxdigit(start[i+1]);
2238                 *q++ = val;
2239             }
2240
2241             ptrs[ptrnum] = q;
2242         }
2243
2244         if (!strcmp(buf, "mul")) {
2245             Bignum a, b, c, p;
2246
2247             if (ptrnum != 3) {
2248                 printf("%d: mul with %d parameters, expected 3\n", line, ptrnum);
2249                 exit(1);
2250             }
2251             a = bignum_from_bytes(ptrs[0], ptrs[1]-ptrs[0]);
2252             b = bignum_from_bytes(ptrs[1], ptrs[2]-ptrs[1]);
2253             c = bignum_from_bytes(ptrs[2], ptrs[3]-ptrs[2]);
2254             p = bigmul(a, b);
2255
2256             if (bignum_cmp(c, p) == 0) {
2257                 passes++;
2258             } else {
2259                 char *as = bignum_decimal(a);
2260                 char *bs = bignum_decimal(b);
2261                 char *cs = bignum_decimal(c);
2262                 char *ps = bignum_decimal(p);
2263
2264                 printf("%d: fail: %s * %s gave %s expected %s\n",
2265                        line, as, bs, ps, cs);
2266                 fails++;
2267
2268                 sfree(as);
2269                 sfree(bs);
2270                 sfree(cs);
2271                 sfree(ps);
2272             }
2273             freebn(a);
2274             freebn(b);
2275             freebn(c);
2276             freebn(p);
2277         } else if (!strcmp(buf, "modmul")) {
2278             Bignum a, b, m, c, p;
2279
2280             if (ptrnum != 4) {
2281                 printf("%d: modmul with %d parameters, expected 4\n",
2282                        line, ptrnum);
2283                 exit(1);
2284             }
2285             a = bignum_from_bytes(ptrs[0], ptrs[1]-ptrs[0]);
2286             b = bignum_from_bytes(ptrs[1], ptrs[2]-ptrs[1]);
2287             m = bignum_from_bytes(ptrs[2], ptrs[3]-ptrs[2]);
2288             c = bignum_from_bytes(ptrs[3], ptrs[4]-ptrs[3]);
2289             p = modmul(a, b, m);
2290
2291             if (bignum_cmp(c, p) == 0) {
2292                 passes++;
2293             } else {
2294                 char *as = bignum_decimal(a);
2295                 char *bs = bignum_decimal(b);
2296                 char *ms = bignum_decimal(m);
2297                 char *cs = bignum_decimal(c);
2298                 char *ps = bignum_decimal(p);
2299
2300                 printf("%d: fail: %s * %s mod %s gave %s expected %s\n",
2301                        line, as, bs, ms, ps, cs);
2302                 fails++;
2303
2304                 sfree(as);
2305                 sfree(bs);
2306                 sfree(ms);
2307                 sfree(cs);
2308                 sfree(ps);
2309             }
2310             freebn(a);
2311             freebn(b);
2312             freebn(m);
2313             freebn(c);
2314             freebn(p);
2315         } else if (!strcmp(buf, "pow")) {
2316             Bignum base, expt, modulus, expected, answer;
2317
2318             if (ptrnum != 4) {
2319                 printf("%d: pow with %d parameters, expected 4\n", line, ptrnum);
2320                 exit(1);
2321             }
2322
2323             base = bignum_from_bytes(ptrs[0], ptrs[1]-ptrs[0]);
2324             expt = bignum_from_bytes(ptrs[1], ptrs[2]-ptrs[1]);
2325             modulus = bignum_from_bytes(ptrs[2], ptrs[3]-ptrs[2]);
2326             expected = bignum_from_bytes(ptrs[3], ptrs[4]-ptrs[3]);
2327             answer = modpow(base, expt, modulus);
2328
2329             if (bignum_cmp(expected, answer) == 0) {
2330                 passes++;
2331             } else {
2332                 char *as = bignum_decimal(base);
2333                 char *bs = bignum_decimal(expt);
2334                 char *cs = bignum_decimal(modulus);
2335                 char *ds = bignum_decimal(answer);
2336                 char *ps = bignum_decimal(expected);
2337
2338                 printf("%d: fail: %s ^ %s mod %s gave %s expected %s\n",
2339                        line, as, bs, cs, ds, ps);
2340                 fails++;
2341
2342                 sfree(as);
2343                 sfree(bs);
2344                 sfree(cs);
2345                 sfree(ds);
2346                 sfree(ps);
2347             }
2348             freebn(base);
2349             freebn(expt);
2350             freebn(modulus);
2351             freebn(expected);
2352             freebn(answer);
2353         } else if (!strcmp(buf, "divmod")) {
2354             Bignum n, d, expect_q, expect_r, answer_q, answer_r;
2355             int fail;
2356
2357             if (ptrnum != 4) {
2358                 printf("%d: divmod with %d parameters, expected 4\n", line, ptrnum);
2359                 exit(1);
2360             }
2361
2362             n = bignum_from_bytes(ptrs[0], ptrs[1]-ptrs[0]);
2363             d = bignum_from_bytes(ptrs[1], ptrs[2]-ptrs[1]);
2364             expect_q = bignum_from_bytes(ptrs[2], ptrs[3]-ptrs[2]);
2365             expect_r = bignum_from_bytes(ptrs[3], ptrs[4]-ptrs[3]);
2366             answer_q = bigdiv(n, d);
2367             answer_r = bigmod(n, d);
2368
2369             fail = FALSE;
2370             if (bignum_cmp(expect_q, answer_q) != 0) {
2371                 char *as = bignum_decimal(n);
2372                 char *bs = bignum_decimal(d);
2373                 char *cs = bignum_decimal(answer_q);
2374                 char *ds = bignum_decimal(expect_q);
2375
2376                 printf("%d: fail: %s / %s gave %s expected %s\n",
2377                        line, as, bs, cs, ds);
2378                 fail = TRUE;
2379
2380                 sfree(as);
2381                 sfree(bs);
2382                 sfree(cs);
2383                 sfree(ds);
2384             }
2385             if (bignum_cmp(expect_r, answer_r) != 0) {
2386                 char *as = bignum_decimal(n);
2387                 char *bs = bignum_decimal(d);
2388                 char *cs = bignum_decimal(answer_r);
2389                 char *ds = bignum_decimal(expect_r);
2390
2391                 printf("%d: fail: %s mod %s gave %s expected %s\n",
2392                        line, as, bs, cs, ds);
2393                 fail = TRUE;
2394
2395                 sfree(as);
2396                 sfree(bs);
2397                 sfree(cs);
2398                 sfree(ds);
2399             }
2400
2401             freebn(n);
2402             freebn(d);
2403             freebn(expect_q);
2404             freebn(expect_r);
2405             freebn(answer_q);
2406             freebn(answer_r);
2407
2408             if (fail)
2409                 fails++;
2410             else
2411                 passes++;
2412         } else {
2413             printf("%d: unrecognised test keyword: '%s'\n", line, buf);
2414             exit(1);
2415         }
2416
2417         sfree(buf);
2418         sfree(data);
2419     }
2420
2421     printf("passed %d failed %d total %d\n", passes, fails, passes+fails);
2422     return fails != 0;
2423 }
2424
2425 #endif