Home | History | Annotate | Download | only in demo
      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 
     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 
     21 static void draw(mp_int * a)
     22 {
     23    ndraw(a, "");
     24 }
     25 
     26 
     27 unsigned long lfsr = 0xAAAAAAAAUL;
     28 
     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 
     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];
     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