OSDN Git Service

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