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