| blob | 51cecdf2b0247725803e92e5c2e8affcc08876ec |
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