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