• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1 #include <time.h>
2 
3 #ifdef IOWNANATHLON
4 #include <unistd.h>
5 #define SLEEP sleep(4)
6 #else
7 #define SLEEP
8 #endif
9 
10 #include "tommath.h"
11 
ndraw(mp_int * a,char * name)12 void ndraw(mp_int * a, char *name)
13 {
14    char buf[16000];
15 
16    printf("%s: ", name);
17    mp_toradix(a, buf, 10);
18    printf("%s\n", buf);
19 }
20 
draw(mp_int * a)21 static void draw(mp_int * a)
22 {
23    ndraw(a, "");
24 }
25 
26 
27 unsigned long lfsr = 0xAAAAAAAAUL;
28 
lbit(void)29 int lbit(void)
30 {
31    if (lfsr & 0x80000000UL) {
32       lfsr = ((lfsr << 1) ^ 0x8000001BUL) & 0xFFFFFFFFUL;
33       return 1;
34    } else {
35       lfsr <<= 1;
36       return 0;
37    }
38 }
39 
myrng(unsigned char * dst,int len,void * dat)40 int myrng(unsigned char *dst, int len, void *dat)
41 {
42    int x;
43 
44    for (x = 0; x < len; x++)
45       dst[x] = rand() & 0xFF;
46    return len;
47 }
48 
49 
50 
51 char cmd[4096], buf[4096];
main(void)52 int main(void)
53 {
54    mp_int a, b, c, d, e, f;
55    unsigned long expt_n, add_n, sub_n, mul_n, div_n, sqr_n, mul2d_n, div2d_n,
56       gcd_n, lcm_n, inv_n, div2_n, mul2_n, add_d_n, sub_d_n, t;
57    unsigned rr;
58    int i, n, err, cnt, ix, old_kara_m, old_kara_s;
59    mp_digit mp;
60 
61 
62    mp_init(&a);
63    mp_init(&b);
64    mp_init(&c);
65    mp_init(&d);
66    mp_init(&e);
67    mp_init(&f);
68 
69    srand(time(NULL));
70 
71 #if 0
72    // test montgomery
73    printf("Testing montgomery...\n");
74    for (i = 1; i < 10; i++) {
75       printf("Testing digit size: %d\n", i);
76       for (n = 0; n < 1000; n++) {
77          mp_rand(&a, i);
78          a.dp[0] |= 1;
79 
80          // let's see if R is right
81          mp_montgomery_calc_normalization(&b, &a);
82          mp_montgomery_setup(&a, &mp);
83 
84          // now test a random reduction
85          for (ix = 0; ix < 100; ix++) {
86              mp_rand(&c, 1 + abs(rand()) % (2*i));
87              mp_copy(&c, &d);
88              mp_copy(&c, &e);
89 
90              mp_mod(&d, &a, &d);
91              mp_montgomery_reduce(&c, &a, mp);
92              mp_mulmod(&c, &b, &a, &c);
93 
94              if (mp_cmp(&c, &d) != MP_EQ) {
95 printf("d = e mod a, c = e MOD a\n");
96 mp_todecimal(&a, buf); printf("a = %s\n", buf);
97 mp_todecimal(&e, buf); printf("e = %s\n", buf);
98 mp_todecimal(&d, buf); printf("d = %s\n", buf);
99 mp_todecimal(&c, buf); printf("c = %s\n", buf);
100 printf("compare no compare!\n"); exit(EXIT_FAILURE); }
101          }
102       }
103    }
104    printf("done\n");
105 
106    // test mp_get_int
107    printf("Testing: mp_get_int\n");
108    for (i = 0; i < 1000; ++i) {
109       t = ((unsigned long) rand() * rand() + 1) & 0xFFFFFFFF;
110       mp_set_int(&a, t);
111       if (t != mp_get_int(&a)) {
112 	 printf("mp_get_int() bad result!\n");
113 	 return 1;
114       }
115    }
116    mp_set_int(&a, 0);
117    if (mp_get_int(&a) != 0) {
118       printf("mp_get_int() bad result!\n");
119       return 1;
120    }
121    mp_set_int(&a, 0xffffffff);
122    if (mp_get_int(&a) != 0xffffffff) {
123       printf("mp_get_int() bad result!\n");
124       return 1;
125    }
126    // test mp_sqrt
127    printf("Testing: mp_sqrt\n");
128    for (i = 0; i < 1000; ++i) {
129       printf("%6d\r", i);
130       fflush(stdout);
131       n = (rand() & 15) + 1;
132       mp_rand(&a, n);
133       if (mp_sqrt(&a, &b) != MP_OKAY) {
134 	 printf("mp_sqrt() error!\n");
135 	 return 1;
136       }
137       mp_n_root(&a, 2, &a);
138       if (mp_cmp_mag(&b, &a) != MP_EQ) {
139 	 printf("mp_sqrt() bad result!\n");
140 	 return 1;
141       }
142    }
143 
144    printf("\nTesting: mp_is_square\n");
145    for (i = 0; i < 1000; ++i) {
146       printf("%6d\r", i);
147       fflush(stdout);
148 
149       /* test mp_is_square false negatives */
150       n = (rand() & 7) + 1;
151       mp_rand(&a, n);
152       mp_sqr(&a, &a);
153       if (mp_is_square(&a, &n) != MP_OKAY) {
154 	 printf("fn:mp_is_square() error!\n");
155 	 return 1;
156       }
157       if (n == 0) {
158 	 printf("fn:mp_is_square() bad result!\n");
159 	 return 1;
160       }
161 
162       /* test for false positives */
163       mp_add_d(&a, 1, &a);
164       if (mp_is_square(&a, &n) != MP_OKAY) {
165 	 printf("fp:mp_is_square() error!\n");
166 	 return 1;
167       }
168       if (n == 1) {
169 	 printf("fp:mp_is_square() bad result!\n");
170 	 return 1;
171       }
172 
173    }
174    printf("\n\n");
175 
176    /* test for size */
177    for (ix = 10; ix < 128; ix++) {
178       printf("Testing (not safe-prime): %9d bits    \r", ix);
179       fflush(stdout);
180       err =
181 	 mp_prime_random_ex(&a, 8, ix,
182 			    (rand() & 1) ? LTM_PRIME_2MSB_OFF :
183 			    LTM_PRIME_2MSB_ON, myrng, NULL);
184       if (err != MP_OKAY) {
185 	 printf("failed with err code %d\n", err);
186 	 return EXIT_FAILURE;
187       }
188       if (mp_count_bits(&a) != ix) {
189 	 printf("Prime is %d not %d bits!!!\n", mp_count_bits(&a), ix);
190 	 return EXIT_FAILURE;
191       }
192    }
193 
194    for (ix = 16; ix < 128; ix++) {
195       printf("Testing (   safe-prime): %9d bits    \r", ix);
196       fflush(stdout);
197       err =
198 	 mp_prime_random_ex(&a, 8, ix,
199 			    ((rand() & 1) ? LTM_PRIME_2MSB_OFF :
200 			     LTM_PRIME_2MSB_ON) | LTM_PRIME_SAFE, myrng,
201 			    NULL);
202       if (err != MP_OKAY) {
203 	 printf("failed with err code %d\n", err);
204 	 return EXIT_FAILURE;
205       }
206       if (mp_count_bits(&a) != ix) {
207 	 printf("Prime is %d not %d bits!!!\n", mp_count_bits(&a), ix);
208 	 return EXIT_FAILURE;
209       }
210       /* let's see if it's really a safe prime */
211       mp_sub_d(&a, 1, &a);
212       mp_div_2(&a, &a);
213       mp_prime_is_prime(&a, 8, &cnt);
214       if (cnt != MP_YES) {
215 	 printf("sub is not prime!\n");
216 	 return EXIT_FAILURE;
217       }
218    }
219 
220    printf("\n\n");
221 
222    mp_read_radix(&a, "123456", 10);
223    mp_toradix_n(&a, buf, 10, 3);
224    printf("a == %s\n", buf);
225    mp_toradix_n(&a, buf, 10, 4);
226    printf("a == %s\n", buf);
227    mp_toradix_n(&a, buf, 10, 30);
228    printf("a == %s\n", buf);
229 
230 
231 #if 0
232    for (;;) {
233       fgets(buf, sizeof(buf), stdin);
234       mp_read_radix(&a, buf, 10);
235       mp_prime_next_prime(&a, 5, 1);
236       mp_toradix(&a, buf, 10);
237       printf("%s, %lu\n", buf, a.dp[0] & 3);
238    }
239 #endif
240 
241    /* test mp_cnt_lsb */
242    printf("testing mp_cnt_lsb...\n");
243    mp_set(&a, 1);
244    for (ix = 0; ix < 1024; ix++) {
245       if (mp_cnt_lsb(&a) != ix) {
246 	 printf("Failed at %d, %d\n", ix, mp_cnt_lsb(&a));
247 	 return 0;
248       }
249       mp_mul_2(&a, &a);
250    }
251 
252 /* test mp_reduce_2k */
253    printf("Testing mp_reduce_2k...\n");
254    for (cnt = 3; cnt <= 128; ++cnt) {
255       mp_digit tmp;
256 
257       mp_2expt(&a, cnt);
258       mp_sub_d(&a, 2, &a);	/* a = 2**cnt - 2 */
259 
260 
261       printf("\nTesting %4d bits", cnt);
262       printf("(%d)", mp_reduce_is_2k(&a));
263       mp_reduce_2k_setup(&a, &tmp);
264       printf("(%d)", tmp);
265       for (ix = 0; ix < 1000; ix++) {
266 	 if (!(ix & 127)) {
267 	    printf(".");
268 	    fflush(stdout);
269 	 }
270 	 mp_rand(&b, (cnt / DIGIT_BIT + 1) * 2);
271 	 mp_copy(&c, &b);
272 	 mp_mod(&c, &a, &c);
273 	 mp_reduce_2k(&b, &a, 2);
274 	 if (mp_cmp(&c, &b)) {
275 	    printf("FAILED\n");
276 	    exit(0);
277 	 }
278       }
279    }
280 
281 /* test mp_div_3  */
282    printf("Testing mp_div_3...\n");
283    mp_set(&d, 3);
284    for (cnt = 0; cnt < 10000;) {
285       mp_digit r1, r2;
286 
287       if (!(++cnt & 127))
288 	 printf("%9d\r", cnt);
289       mp_rand(&a, abs(rand()) % 128 + 1);
290       mp_div(&a, &d, &b, &e);
291       mp_div_3(&a, &c, &r2);
292 
293       if (mp_cmp(&b, &c) || mp_cmp_d(&e, r2)) {
294 	 printf("\n\nmp_div_3 => Failure\n");
295       }
296    }
297    printf("\n\nPassed div_3 testing\n");
298 
299 /* test the DR reduction */
300    printf("testing mp_dr_reduce...\n");
301    for (cnt = 2; cnt < 32; cnt++) {
302       printf("%d digit modulus\n", cnt);
303       mp_grow(&a, cnt);
304       mp_zero(&a);
305       for (ix = 1; ix < cnt; ix++) {
306 	 a.dp[ix] = MP_MASK;
307       }
308       a.used = cnt;
309       a.dp[0] = 3;
310 
311       mp_rand(&b, cnt - 1);
312       mp_copy(&b, &c);
313 
314       rr = 0;
315       do {
316 	 if (!(rr & 127)) {
317 	    printf("%9lu\r", rr);
318 	    fflush(stdout);
319 	 }
320 	 mp_sqr(&b, &b);
321 	 mp_add_d(&b, 1, &b);
322 	 mp_copy(&b, &c);
323 
324 	 mp_mod(&b, &a, &b);
325 	 mp_dr_reduce(&c, &a, (((mp_digit) 1) << DIGIT_BIT) - a.dp[0]);
326 
327 	 if (mp_cmp(&b, &c) != MP_EQ) {
328 	    printf("Failed on trial %lu\n", rr);
329 	    exit(-1);
330 
331 	 }
332       } while (++rr < 500);
333       printf("Passed DR test for %d digits\n", cnt);
334    }
335 
336 #endif
337 
338 /* test the mp_reduce_2k_l code */
339 #if 0
340 #if 0
341 /* first load P with 2^1024 - 0x2A434 B9FDEC95 D8F9D550 FFFFFFFF FFFFFFFF */
342    mp_2expt(&a, 1024);
343    mp_read_radix(&b, "2A434B9FDEC95D8F9D550FFFFFFFFFFFFFFFF", 16);
344    mp_sub(&a, &b, &a);
345 #elif 1
346 /*  p = 2^2048 - 0x1 00000000 00000000 00000000 00000000 4945DDBF 8EA2A91D 5776399B B83E188F  */
347    mp_2expt(&a, 2048);
348    mp_read_radix(&b,
349 		 "1000000000000000000000000000000004945DDBF8EA2A91D5776399BB83E188F",
350 		 16);
351    mp_sub(&a, &b, &a);
352 #endif
353 
354    mp_todecimal(&a, buf);
355    printf("p==%s\n", buf);
356 /* now mp_reduce_is_2k_l() should return */
357    if (mp_reduce_is_2k_l(&a) != 1) {
358       printf("mp_reduce_is_2k_l() return 0, should be 1\n");
359       return EXIT_FAILURE;
360    }
361    mp_reduce_2k_setup_l(&a, &d);
362    /* now do a million square+1 to see if it varies */
363    mp_rand(&b, 64);
364    mp_mod(&b, &a, &b);
365    mp_copy(&b, &c);
366    printf("testing mp_reduce_2k_l...");
367    fflush(stdout);
368    for (cnt = 0; cnt < (1UL << 20); cnt++) {
369       mp_sqr(&b, &b);
370       mp_add_d(&b, 1, &b);
371       mp_reduce_2k_l(&b, &a, &d);
372       mp_sqr(&c, &c);
373       mp_add_d(&c, 1, &c);
374       mp_mod(&c, &a, &c);
375       if (mp_cmp(&b, &c) != MP_EQ) {
376 	 printf("mp_reduce_2k_l() failed at step %lu\n", cnt);
377 	 mp_tohex(&b, buf);
378 	 printf("b == %s\n", buf);
379 	 mp_tohex(&c, buf);
380 	 printf("c == %s\n", buf);
381 	 return EXIT_FAILURE;
382       }
383    }
384    printf("...Passed\n");
385 #endif
386 
387    div2_n = mul2_n = inv_n = expt_n = lcm_n = gcd_n = add_n =
388       sub_n = mul_n = div_n = sqr_n = mul2d_n = div2d_n = cnt = add_d_n =
389       sub_d_n = 0;
390 
391    /* force KARA and TOOM to enable despite cutoffs */
392    KARATSUBA_SQR_CUTOFF = KARATSUBA_MUL_CUTOFF = 8;
393    TOOM_SQR_CUTOFF = TOOM_MUL_CUTOFF = 16;
394 
395    for (;;) {
396       /* randomly clear and re-init one variable, this has the affect of triming the alloc space */
397       switch (abs(rand()) % 7) {
398       case 0:
399 	 mp_clear(&a);
400 	 mp_init(&a);
401 	 break;
402       case 1:
403 	 mp_clear(&b);
404 	 mp_init(&b);
405 	 break;
406       case 2:
407 	 mp_clear(&c);
408 	 mp_init(&c);
409 	 break;
410       case 3:
411 	 mp_clear(&d);
412 	 mp_init(&d);
413 	 break;
414       case 4:
415 	 mp_clear(&e);
416 	 mp_init(&e);
417 	 break;
418       case 5:
419 	 mp_clear(&f);
420 	 mp_init(&f);
421 	 break;
422       case 6:
423 	 break;			/* don't clear any */
424       }
425 
426 
427       printf
428 	 ("%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu ",
429 	  add_n, sub_n, mul_n, div_n, sqr_n, mul2d_n, div2d_n, gcd_n, lcm_n,
430 	  expt_n, inv_n, div2_n, mul2_n, add_d_n, sub_d_n);
431       fgets(cmd, 4095, stdin);
432       cmd[strlen(cmd) - 1] = 0;
433       printf("%s  ]\r", cmd);
434       fflush(stdout);
435       if (!strcmp(cmd, "mul2d")) {
436 	 ++mul2d_n;
437 	 fgets(buf, 4095, stdin);
438 	 mp_read_radix(&a, buf, 64);
439 	 fgets(buf, 4095, stdin);
440 	 sscanf(buf, "%d", &rr);
441 	 fgets(buf, 4095, stdin);
442 	 mp_read_radix(&b, buf, 64);
443 
444 	 mp_mul_2d(&a, rr, &a);
445 	 a.sign = b.sign;
446 	 if (mp_cmp(&a, &b) != MP_EQ) {
447 	    printf("mul2d failed, rr == %d\n", rr);
448 	    draw(&a);
449 	    draw(&b);
450 	    return 0;
451 	 }
452       } else if (!strcmp(cmd, "div2d")) {
453 	 ++div2d_n;
454 	 fgets(buf, 4095, stdin);
455 	 mp_read_radix(&a, buf, 64);
456 	 fgets(buf, 4095, stdin);
457 	 sscanf(buf, "%d", &rr);
458 	 fgets(buf, 4095, stdin);
459 	 mp_read_radix(&b, buf, 64);
460 
461 	 mp_div_2d(&a, rr, &a, &e);
462 	 a.sign = b.sign;
463 	 if (a.used == b.used && a.used == 0) {
464 	    a.sign = b.sign = MP_ZPOS;
465 	 }
466 	 if (mp_cmp(&a, &b) != MP_EQ) {
467 	    printf("div2d failed, rr == %d\n", rr);
468 	    draw(&a);
469 	    draw(&b);
470 	    return 0;
471 	 }
472       } else if (!strcmp(cmd, "add")) {
473 	 ++add_n;
474 	 fgets(buf, 4095, stdin);
475 	 mp_read_radix(&a, buf, 64);
476 	 fgets(buf, 4095, stdin);
477 	 mp_read_radix(&b, buf, 64);
478 	 fgets(buf, 4095, stdin);
479 	 mp_read_radix(&c, buf, 64);
480 	 mp_copy(&a, &d);
481 	 mp_add(&d, &b, &d);
482 	 if (mp_cmp(&c, &d) != MP_EQ) {
483 	    printf("add %lu failure!\n", add_n);
484 	    draw(&a);
485 	    draw(&b);
486 	    draw(&c);
487 	    draw(&d);
488 	    return 0;
489 	 }
490 
491 	 /* test the sign/unsigned storage functions */
492 
493 	 rr = mp_signed_bin_size(&c);
494 	 mp_to_signed_bin(&c, (unsigned char *) cmd);
495 	 memset(cmd + rr, rand() & 255, sizeof(cmd) - rr);
496 	 mp_read_signed_bin(&d, (unsigned char *) cmd, rr);
497 	 if (mp_cmp(&c, &d) != MP_EQ) {
498 	    printf("mp_signed_bin failure!\n");
499 	    draw(&c);
500 	    draw(&d);
501 	    return 0;
502 	 }
503 
504 
505 	 rr = mp_unsigned_bin_size(&c);
506 	 mp_to_unsigned_bin(&c, (unsigned char *) cmd);
507 	 memset(cmd + rr, rand() & 255, sizeof(cmd) - rr);
508 	 mp_read_unsigned_bin(&d, (unsigned char *) cmd, rr);
509 	 if (mp_cmp_mag(&c, &d) != MP_EQ) {
510 	    printf("mp_unsigned_bin failure!\n");
511 	    draw(&c);
512 	    draw(&d);
513 	    return 0;
514 	 }
515 
516       } else if (!strcmp(cmd, "sub")) {
517 	 ++sub_n;
518 	 fgets(buf, 4095, stdin);
519 	 mp_read_radix(&a, buf, 64);
520 	 fgets(buf, 4095, stdin);
521 	 mp_read_radix(&b, buf, 64);
522 	 fgets(buf, 4095, stdin);
523 	 mp_read_radix(&c, buf, 64);
524 	 mp_copy(&a, &d);
525 	 mp_sub(&d, &b, &d);
526 	 if (mp_cmp(&c, &d) != MP_EQ) {
527 	    printf("sub %lu failure!\n", sub_n);
528 	    draw(&a);
529 	    draw(&b);
530 	    draw(&c);
531 	    draw(&d);
532 	    return 0;
533 	 }
534       } else if (!strcmp(cmd, "mul")) {
535 	 ++mul_n;
536 	 fgets(buf, 4095, stdin);
537 	 mp_read_radix(&a, buf, 64);
538 	 fgets(buf, 4095, stdin);
539 	 mp_read_radix(&b, buf, 64);
540 	 fgets(buf, 4095, stdin);
541 	 mp_read_radix(&c, buf, 64);
542 	 mp_copy(&a, &d);
543 	 mp_mul(&d, &b, &d);
544 	 if (mp_cmp(&c, &d) != MP_EQ) {
545 	    printf("mul %lu failure!\n", mul_n);
546 	    draw(&a);
547 	    draw(&b);
548 	    draw(&c);
549 	    draw(&d);
550 	    return 0;
551 	 }
552       } else if (!strcmp(cmd, "div")) {
553 	 ++div_n;
554 	 fgets(buf, 4095, stdin);
555 	 mp_read_radix(&a, buf, 64);
556 	 fgets(buf, 4095, stdin);
557 	 mp_read_radix(&b, buf, 64);
558 	 fgets(buf, 4095, stdin);
559 	 mp_read_radix(&c, buf, 64);
560 	 fgets(buf, 4095, stdin);
561 	 mp_read_radix(&d, buf, 64);
562 
563 	 mp_div(&a, &b, &e, &f);
564 	 if (mp_cmp(&c, &e) != MP_EQ || mp_cmp(&d, &f) != MP_EQ) {
565 	    printf("div %lu %d, %d, failure!\n", div_n, mp_cmp(&c, &e),
566 		   mp_cmp(&d, &f));
567 	    draw(&a);
568 	    draw(&b);
569 	    draw(&c);
570 	    draw(&d);
571 	    draw(&e);
572 	    draw(&f);
573 	    return 0;
574 	 }
575 
576       } else if (!strcmp(cmd, "sqr")) {
577 	 ++sqr_n;
578 	 fgets(buf, 4095, stdin);
579 	 mp_read_radix(&a, buf, 64);
580 	 fgets(buf, 4095, stdin);
581 	 mp_read_radix(&b, buf, 64);
582 	 mp_copy(&a, &c);
583 	 mp_sqr(&c, &c);
584 	 if (mp_cmp(&b, &c) != MP_EQ) {
585 	    printf("sqr %lu failure!\n", sqr_n);
586 	    draw(&a);
587 	    draw(&b);
588 	    draw(&c);
589 	    return 0;
590 	 }
591       } else if (!strcmp(cmd, "gcd")) {
592 	 ++gcd_n;
593 	 fgets(buf, 4095, stdin);
594 	 mp_read_radix(&a, buf, 64);
595 	 fgets(buf, 4095, stdin);
596 	 mp_read_radix(&b, buf, 64);
597 	 fgets(buf, 4095, stdin);
598 	 mp_read_radix(&c, buf, 64);
599 	 mp_copy(&a, &d);
600 	 mp_gcd(&d, &b, &d);
601 	 d.sign = c.sign;
602 	 if (mp_cmp(&c, &d) != MP_EQ) {
603 	    printf("gcd %lu failure!\n", gcd_n);
604 	    draw(&a);
605 	    draw(&b);
606 	    draw(&c);
607 	    draw(&d);
608 	    return 0;
609 	 }
610       } else if (!strcmp(cmd, "lcm")) {
611 	 ++lcm_n;
612 	 fgets(buf, 4095, stdin);
613 	 mp_read_radix(&a, buf, 64);
614 	 fgets(buf, 4095, stdin);
615 	 mp_read_radix(&b, buf, 64);
616 	 fgets(buf, 4095, stdin);
617 	 mp_read_radix(&c, buf, 64);
618 	 mp_copy(&a, &d);
619 	 mp_lcm(&d, &b, &d);
620 	 d.sign = c.sign;
621 	 if (mp_cmp(&c, &d) != MP_EQ) {
622 	    printf("lcm %lu failure!\n", lcm_n);
623 	    draw(&a);
624 	    draw(&b);
625 	    draw(&c);
626 	    draw(&d);
627 	    return 0;
628 	 }
629       } else if (!strcmp(cmd, "expt")) {
630 	 ++expt_n;
631 	 fgets(buf, 4095, stdin);
632 	 mp_read_radix(&a, buf, 64);
633 	 fgets(buf, 4095, stdin);
634 	 mp_read_radix(&b, buf, 64);
635 	 fgets(buf, 4095, stdin);
636 	 mp_read_radix(&c, buf, 64);
637 	 fgets(buf, 4095, stdin);
638 	 mp_read_radix(&d, buf, 64);
639 	 mp_copy(&a, &e);
640 	 mp_exptmod(&e, &b, &c, &e);
641 	 if (mp_cmp(&d, &e) != MP_EQ) {
642 	    printf("expt %lu failure!\n", expt_n);
643 	    draw(&a);
644 	    draw(&b);
645 	    draw(&c);
646 	    draw(&d);
647 	    draw(&e);
648 	    return 0;
649 	 }
650       } else if (!strcmp(cmd, "invmod")) {
651 	 ++inv_n;
652 	 fgets(buf, 4095, stdin);
653 	 mp_read_radix(&a, buf, 64);
654 	 fgets(buf, 4095, stdin);
655 	 mp_read_radix(&b, buf, 64);
656 	 fgets(buf, 4095, stdin);
657 	 mp_read_radix(&c, buf, 64);
658 	 mp_invmod(&a, &b, &d);
659 	 mp_mulmod(&d, &a, &b, &e);
660 	 if (mp_cmp_d(&e, 1) != MP_EQ) {
661 	    printf("inv [wrong value from MPI?!] failure\n");
662 	    draw(&a);
663 	    draw(&b);
664 	    draw(&c);
665 	    draw(&d);
666 	    mp_gcd(&a, &b, &e);
667 	    draw(&e);
668 	    return 0;
669 	 }
670 
671       } else if (!strcmp(cmd, "div2")) {
672 	 ++div2_n;
673 	 fgets(buf, 4095, stdin);
674 	 mp_read_radix(&a, buf, 64);
675 	 fgets(buf, 4095, stdin);
676 	 mp_read_radix(&b, buf, 64);
677 	 mp_div_2(&a, &c);
678 	 if (mp_cmp(&c, &b) != MP_EQ) {
679 	    printf("div_2 %lu failure\n", div2_n);
680 	    draw(&a);
681 	    draw(&b);
682 	    draw(&c);
683 	    return 0;
684 	 }
685       } else if (!strcmp(cmd, "mul2")) {
686 	 ++mul2_n;
687 	 fgets(buf, 4095, stdin);
688 	 mp_read_radix(&a, buf, 64);
689 	 fgets(buf, 4095, stdin);
690 	 mp_read_radix(&b, buf, 64);
691 	 mp_mul_2(&a, &c);
692 	 if (mp_cmp(&c, &b) != MP_EQ) {
693 	    printf("mul_2 %lu failure\n", mul2_n);
694 	    draw(&a);
695 	    draw(&b);
696 	    draw(&c);
697 	    return 0;
698 	 }
699       } else if (!strcmp(cmd, "add_d")) {
700 	 ++add_d_n;
701 	 fgets(buf, 4095, stdin);
702 	 mp_read_radix(&a, buf, 64);
703 	 fgets(buf, 4095, stdin);
704 	 sscanf(buf, "%d", &ix);
705 	 fgets(buf, 4095, stdin);
706 	 mp_read_radix(&b, buf, 64);
707 	 mp_add_d(&a, ix, &c);
708 	 if (mp_cmp(&b, &c) != MP_EQ) {
709 	    printf("add_d %lu failure\n", add_d_n);
710 	    draw(&a);
711 	    draw(&b);
712 	    draw(&c);
713 	    printf("d == %d\n", ix);
714 	    return 0;
715 	 }
716       } else if (!strcmp(cmd, "sub_d")) {
717 	 ++sub_d_n;
718 	 fgets(buf, 4095, stdin);
719 	 mp_read_radix(&a, buf, 64);
720 	 fgets(buf, 4095, stdin);
721 	 sscanf(buf, "%d", &ix);
722 	 fgets(buf, 4095, stdin);
723 	 mp_read_radix(&b, buf, 64);
724 	 mp_sub_d(&a, ix, &c);
725 	 if (mp_cmp(&b, &c) != MP_EQ) {
726 	    printf("sub_d %lu failure\n", sub_d_n);
727 	    draw(&a);
728 	    draw(&b);
729 	    draw(&c);
730 	    printf("d == %d\n", ix);
731 	    return 0;
732 	 }
733       }
734    }
735    return 0;
736 }
737 
738 /* $Source: /cvs/libtom/libtommath/demo/demo.c,v $ */
739 /* $Revision: 1.3 $ */
740 /* $Date: 2005/06/24 11:32:07 $ */
741