1 /* $NetBSD: demo.c,v 1.1.1.2 2014/04/24 12:45:39 pettai Exp $ */
14 void ndraw(mp_int
* a
, char *name
)
19 mp_toradix(a
, buf
, 10);
23 static void draw(mp_int
* a
)
29 unsigned long lfsr
= 0xAAAAAAAAUL
;
33 if (lfsr
& 0x80000000UL
) {
34 lfsr
= ((lfsr
<< 1) ^ 0x8000001BUL
) & 0xFFFFFFFFUL
;
42 int myrng(unsigned char *dst
, int len
, void *dat
)
46 for (x
= 0; x
< len
; x
++)
47 dst
[x
] = rand() & 0xFF;
53 char cmd
[4096], buf
[4096];
56 mp_int a
, b
, c
, d
, e
, f
;
57 unsigned long expt_n
, add_n
, sub_n
, mul_n
, div_n
, sqr_n
, mul2d_n
, div2d_n
,
58 gcd_n
, lcm_n
, inv_n
, div2_n
, mul2_n
, add_d_n
, sub_d_n
, t
;
60 int i
, n
, err
, cnt
, ix
, old_kara_m
, old_kara_s
;
75 printf("Testing montgomery...\n");
76 for (i
= 1; i
< 10; i
++) {
77 printf("Testing digit size: %d\n", i
);
78 for (n
= 0; n
< 1000; n
++) {
82 // let's see if R is right
83 mp_montgomery_calc_normalization(&b
, &a
);
84 mp_montgomery_setup(&a
, &mp
);
86 // now test a random reduction
87 for (ix
= 0; ix
< 100; ix
++) {
88 mp_rand(&c
, 1 + abs(rand()) % (2*i
));
93 mp_montgomery_reduce(&c
, &a
, mp
);
94 mp_mulmod(&c
, &b
, &a
, &c
);
96 if (mp_cmp(&c
, &d
) != MP_EQ
) {
97 printf("d = e mod a, c = e MOD a\n");
98 mp_todecimal(&a
, buf
); printf("a = %s\n", buf
);
99 mp_todecimal(&e
, buf
); printf("e = %s\n", buf
);
100 mp_todecimal(&d
, buf
); printf("d = %s\n", buf
);
101 mp_todecimal(&c
, buf
); printf("c = %s\n", buf
);
102 printf("compare no compare!\n"); exit(EXIT_FAILURE
); }
109 printf("Testing: mp_get_int\n");
110 for (i
= 0; i
< 1000; ++i
) {
111 t
= ((unsigned long) rand() * rand() + 1) & 0xFFFFFFFF;
113 if (t
!= mp_get_int(&a
)) {
114 printf("mp_get_int() bad result!\n");
119 if (mp_get_int(&a
) != 0) {
120 printf("mp_get_int() bad result!\n");
123 mp_set_int(&a
, 0xffffffff);
124 if (mp_get_int(&a
) != 0xffffffff) {
125 printf("mp_get_int() bad result!\n");
129 printf("Testing: mp_sqrt\n");
130 for (i
= 0; i
< 1000; ++i
) {
133 n
= (rand() & 15) + 1;
135 if (mp_sqrt(&a
, &b
) != MP_OKAY
) {
136 printf("mp_sqrt() error!\n");
139 mp_n_root(&a
, 2, &a
);
140 if (mp_cmp_mag(&b
, &a
) != MP_EQ
) {
141 printf("mp_sqrt() bad result!\n");
146 printf("\nTesting: mp_is_square\n");
147 for (i
= 0; i
< 1000; ++i
) {
151 /* test mp_is_square false negatives */
152 n
= (rand() & 7) + 1;
155 if (mp_is_square(&a
, &n
) != MP_OKAY
) {
156 printf("fn:mp_is_square() error!\n");
160 printf("fn:mp_is_square() bad result!\n");
164 /* test for false positives */
166 if (mp_is_square(&a
, &n
) != MP_OKAY
) {
167 printf("fp:mp_is_square() error!\n");
171 printf("fp:mp_is_square() bad result!\n");
179 for (ix
= 10; ix
< 128; ix
++) {
180 printf("Testing (not safe-prime): %9d bits \r", ix
);
183 mp_prime_random_ex(&a
, 8, ix
,
184 (rand() & 1) ? LTM_PRIME_2MSB_OFF
:
185 LTM_PRIME_2MSB_ON
, myrng
, NULL
);
186 if (err
!= MP_OKAY
) {
187 printf("failed with err code %d\n", err
);
190 if (mp_count_bits(&a
) != ix
) {
191 printf("Prime is %d not %d bits!!!\n", mp_count_bits(&a
), ix
);
196 for (ix
= 16; ix
< 128; ix
++) {
197 printf("Testing ( safe-prime): %9d bits \r", ix
);
200 mp_prime_random_ex(&a
, 8, ix
,
201 ((rand() & 1) ? LTM_PRIME_2MSB_OFF
:
202 LTM_PRIME_2MSB_ON
) | LTM_PRIME_SAFE
, myrng
,
204 if (err
!= MP_OKAY
) {
205 printf("failed with err code %d\n", err
);
208 if (mp_count_bits(&a
) != ix
) {
209 printf("Prime is %d not %d bits!!!\n", mp_count_bits(&a
), ix
);
212 /* let's see if it's really a safe prime */
215 mp_prime_is_prime(&a
, 8, &cnt
);
217 printf("sub is not prime!\n");
224 mp_read_radix(&a
, "123456", 10);
225 mp_toradix_n(&a
, buf
, 10, 3);
226 printf("a == %s\n", buf
);
227 mp_toradix_n(&a
, buf
, 10, 4);
228 printf("a == %s\n", buf
);
229 mp_toradix_n(&a
, buf
, 10, 30);
230 printf("a == %s\n", buf
);
235 fgets(buf
, sizeof(buf
), stdin
);
236 mp_read_radix(&a
, buf
, 10);
237 mp_prime_next_prime(&a
, 5, 1);
238 mp_toradix(&a
, buf
, 10);
239 printf("%s, %lu\n", buf
, a
.dp
[0] & 3);
243 /* test mp_cnt_lsb */
244 printf("testing mp_cnt_lsb...\n");
246 for (ix
= 0; ix
< 1024; ix
++) {
247 if (mp_cnt_lsb(&a
) != ix
) {
248 printf("Failed at %d, %d\n", ix
, mp_cnt_lsb(&a
));
254 /* test mp_reduce_2k */
255 printf("Testing mp_reduce_2k...\n");
256 for (cnt
= 3; cnt
<= 128; ++cnt
) {
260 mp_sub_d(&a
, 2, &a
); /* a = 2**cnt - 2 */
263 printf("\nTesting %4d bits", cnt
);
264 printf("(%d)", mp_reduce_is_2k(&a
));
265 mp_reduce_2k_setup(&a
, &tmp
);
267 for (ix
= 0; ix
< 1000; ix
++) {
272 mp_rand(&b
, (cnt
/ DIGIT_BIT
+ 1) * 2);
275 mp_reduce_2k(&b
, &a
, 2);
276 if (mp_cmp(&c
, &b
)) {
284 printf("Testing mp_div_3...\n");
286 for (cnt
= 0; cnt
< 10000;) {
290 printf("%9d\r", cnt
);
291 mp_rand(&a
, abs(rand()) % 128 + 1);
292 mp_div(&a
, &d
, &b
, &e
);
293 mp_div_3(&a
, &c
, &r2
);
295 if (mp_cmp(&b
, &c
) || mp_cmp_d(&e
, r2
)) {
296 printf("\n\nmp_div_3 => Failure\n");
299 printf("\n\nPassed div_3 testing\n");
301 /* test the DR reduction */
302 printf("testing mp_dr_reduce...\n");
303 for (cnt
= 2; cnt
< 32; cnt
++) {
304 printf("%d digit modulus\n", cnt
);
307 for (ix
= 1; ix
< cnt
; ix
++) {
313 mp_rand(&b
, cnt
- 1);
319 printf("%9lu\r", rr
);
327 mp_dr_reduce(&c
, &a
, (((mp_digit
) 1) << DIGIT_BIT
) - a
.dp
[0]);
329 if (mp_cmp(&b
, &c
) != MP_EQ
) {
330 printf("Failed on trial %lu\n", rr
);
334 } while (++rr
< 500);
335 printf("Passed DR test for %d digits\n", cnt
);
340 /* test the mp_reduce_2k_l code */
343 /* first load P with 2^1024 - 0x2A434 B9FDEC95 D8F9D550 FFFFFFFF FFFFFFFF */
345 mp_read_radix(&b
, "2A434B9FDEC95D8F9D550FFFFFFFFFFFFFFFF", 16);
348 /* p = 2^2048 - 0x1 00000000 00000000 00000000 00000000 4945DDBF 8EA2A91D 5776399B B83E188F */
351 "1000000000000000000000000000000004945DDBF8EA2A91D5776399BB83E188F",
356 mp_todecimal(&a
, buf
);
357 printf("p==%s\n", buf
);
358 /* now mp_reduce_is_2k_l() should return */
359 if (mp_reduce_is_2k_l(&a
) != 1) {
360 printf("mp_reduce_is_2k_l() return 0, should be 1\n");
363 mp_reduce_2k_setup_l(&a
, &d
);
364 /* now do a million square+1 to see if it varies */
368 printf("testing mp_reduce_2k_l...");
370 for (cnt
= 0; cnt
< (1UL << 20); cnt
++) {
373 mp_reduce_2k_l(&b
, &a
, &d
);
377 if (mp_cmp(&b
, &c
) != MP_EQ
) {
378 printf("mp_reduce_2k_l() failed at step %lu\n", cnt
);
380 printf("b == %s\n", buf
);
382 printf("c == %s\n", buf
);
386 printf("...Passed\n");
389 div2_n
= mul2_n
= inv_n
= expt_n
= lcm_n
= gcd_n
= add_n
=
390 sub_n
= mul_n
= div_n
= sqr_n
= mul2d_n
= div2d_n
= cnt
= add_d_n
=
393 /* force KARA and TOOM to enable despite cutoffs */
394 KARATSUBA_SQR_CUTOFF
= KARATSUBA_MUL_CUTOFF
= 8;
395 TOOM_SQR_CUTOFF
= TOOM_MUL_CUTOFF
= 16;
398 /* randomly clear and re-init one variable, this has the affect of triming the alloc space */
399 switch (abs(rand()) % 7) {
425 break; /* don't clear any */
430 ("%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu ",
431 add_n
, sub_n
, mul_n
, div_n
, sqr_n
, mul2d_n
, div2d_n
, gcd_n
, lcm_n
,
432 expt_n
, inv_n
, div2_n
, mul2_n
, add_d_n
, sub_d_n
);
433 fgets(cmd
, 4095, stdin
);
434 cmd
[strlen(cmd
) - 1] = 0;
435 printf("%s ]\r", cmd
);
437 if (!strcmp(cmd
, "mul2d")) {
439 fgets(buf
, 4095, stdin
);
440 mp_read_radix(&a
, buf
, 64);
441 fgets(buf
, 4095, stdin
);
442 sscanf(buf
, "%d", &rr
);
443 fgets(buf
, 4095, stdin
);
444 mp_read_radix(&b
, buf
, 64);
446 mp_mul_2d(&a
, rr
, &a
);
448 if (mp_cmp(&a
, &b
) != MP_EQ
) {
449 printf("mul2d failed, rr == %d\n", rr
);
454 } else if (!strcmp(cmd
, "div2d")) {
456 fgets(buf
, 4095, stdin
);
457 mp_read_radix(&a
, buf
, 64);
458 fgets(buf
, 4095, stdin
);
459 sscanf(buf
, "%d", &rr
);
460 fgets(buf
, 4095, stdin
);
461 mp_read_radix(&b
, buf
, 64);
463 mp_div_2d(&a
, rr
, &a
, &e
);
465 if (a
.used
== b
.used
&& a
.used
== 0) {
466 a
.sign
= b
.sign
= MP_ZPOS
;
468 if (mp_cmp(&a
, &b
) != MP_EQ
) {
469 printf("div2d failed, rr == %d\n", rr
);
474 } else if (!strcmp(cmd
, "add")) {
476 fgets(buf
, 4095, stdin
);
477 mp_read_radix(&a
, buf
, 64);
478 fgets(buf
, 4095, stdin
);
479 mp_read_radix(&b
, buf
, 64);
480 fgets(buf
, 4095, stdin
);
481 mp_read_radix(&c
, buf
, 64);
484 if (mp_cmp(&c
, &d
) != MP_EQ
) {
485 printf("add %lu failure!\n", add_n
);
493 /* test the sign/unsigned storage functions */
495 rr
= mp_signed_bin_size(&c
);
496 mp_to_signed_bin(&c
, (unsigned char *) cmd
);
497 memset(cmd
+ rr
, rand() & 255, sizeof(cmd
) - rr
);
498 mp_read_signed_bin(&d
, (unsigned char *) cmd
, rr
);
499 if (mp_cmp(&c
, &d
) != MP_EQ
) {
500 printf("mp_signed_bin failure!\n");
507 rr
= mp_unsigned_bin_size(&c
);
508 mp_to_unsigned_bin(&c
, (unsigned char *) cmd
);
509 memset(cmd
+ rr
, rand() & 255, sizeof(cmd
) - rr
);
510 mp_read_unsigned_bin(&d
, (unsigned char *) cmd
, rr
);
511 if (mp_cmp_mag(&c
, &d
) != MP_EQ
) {
512 printf("mp_unsigned_bin failure!\n");
518 } else if (!strcmp(cmd
, "sub")) {
520 fgets(buf
, 4095, stdin
);
521 mp_read_radix(&a
, buf
, 64);
522 fgets(buf
, 4095, stdin
);
523 mp_read_radix(&b
, buf
, 64);
524 fgets(buf
, 4095, stdin
);
525 mp_read_radix(&c
, buf
, 64);
528 if (mp_cmp(&c
, &d
) != MP_EQ
) {
529 printf("sub %lu failure!\n", sub_n
);
536 } else if (!strcmp(cmd
, "mul")) {
538 fgets(buf
, 4095, stdin
);
539 mp_read_radix(&a
, buf
, 64);
540 fgets(buf
, 4095, stdin
);
541 mp_read_radix(&b
, buf
, 64);
542 fgets(buf
, 4095, stdin
);
543 mp_read_radix(&c
, buf
, 64);
546 if (mp_cmp(&c
, &d
) != MP_EQ
) {
547 printf("mul %lu failure!\n", mul_n
);
554 } else if (!strcmp(cmd
, "div")) {
556 fgets(buf
, 4095, stdin
);
557 mp_read_radix(&a
, buf
, 64);
558 fgets(buf
, 4095, stdin
);
559 mp_read_radix(&b
, buf
, 64);
560 fgets(buf
, 4095, stdin
);
561 mp_read_radix(&c
, buf
, 64);
562 fgets(buf
, 4095, stdin
);
563 mp_read_radix(&d
, buf
, 64);
565 mp_div(&a
, &b
, &e
, &f
);
566 if (mp_cmp(&c
, &e
) != MP_EQ
|| mp_cmp(&d
, &f
) != MP_EQ
) {
567 printf("div %lu %d, %d, failure!\n", div_n
, mp_cmp(&c
, &e
),
578 } else if (!strcmp(cmd
, "sqr")) {
580 fgets(buf
, 4095, stdin
);
581 mp_read_radix(&a
, buf
, 64);
582 fgets(buf
, 4095, stdin
);
583 mp_read_radix(&b
, buf
, 64);
586 if (mp_cmp(&b
, &c
) != MP_EQ
) {
587 printf("sqr %lu failure!\n", sqr_n
);
593 } else if (!strcmp(cmd
, "gcd")) {
595 fgets(buf
, 4095, stdin
);
596 mp_read_radix(&a
, buf
, 64);
597 fgets(buf
, 4095, stdin
);
598 mp_read_radix(&b
, buf
, 64);
599 fgets(buf
, 4095, stdin
);
600 mp_read_radix(&c
, buf
, 64);
604 if (mp_cmp(&c
, &d
) != MP_EQ
) {
605 printf("gcd %lu failure!\n", gcd_n
);
612 } else if (!strcmp(cmd
, "lcm")) {
614 fgets(buf
, 4095, stdin
);
615 mp_read_radix(&a
, buf
, 64);
616 fgets(buf
, 4095, stdin
);
617 mp_read_radix(&b
, buf
, 64);
618 fgets(buf
, 4095, stdin
);
619 mp_read_radix(&c
, buf
, 64);
623 if (mp_cmp(&c
, &d
) != MP_EQ
) {
624 printf("lcm %lu failure!\n", lcm_n
);
631 } else if (!strcmp(cmd
, "expt")) {
633 fgets(buf
, 4095, stdin
);
634 mp_read_radix(&a
, buf
, 64);
635 fgets(buf
, 4095, stdin
);
636 mp_read_radix(&b
, buf
, 64);
637 fgets(buf
, 4095, stdin
);
638 mp_read_radix(&c
, buf
, 64);
639 fgets(buf
, 4095, stdin
);
640 mp_read_radix(&d
, buf
, 64);
642 mp_exptmod(&e
, &b
, &c
, &e
);
643 if (mp_cmp(&d
, &e
) != MP_EQ
) {
644 printf("expt %lu failure!\n", expt_n
);
652 } else if (!strcmp(cmd
, "invmod")) {
654 fgets(buf
, 4095, stdin
);
655 mp_read_radix(&a
, buf
, 64);
656 fgets(buf
, 4095, stdin
);
657 mp_read_radix(&b
, buf
, 64);
658 fgets(buf
, 4095, stdin
);
659 mp_read_radix(&c
, buf
, 64);
660 mp_invmod(&a
, &b
, &d
);
661 mp_mulmod(&d
, &a
, &b
, &e
);
662 if (mp_cmp_d(&e
, 1) != MP_EQ
) {
663 printf("inv [wrong value from MPI?!] failure\n");
673 } else if (!strcmp(cmd
, "div2")) {
675 fgets(buf
, 4095, stdin
);
676 mp_read_radix(&a
, buf
, 64);
677 fgets(buf
, 4095, stdin
);
678 mp_read_radix(&b
, buf
, 64);
680 if (mp_cmp(&c
, &b
) != MP_EQ
) {
681 printf("div_2 %lu failure\n", div2_n
);
687 } else if (!strcmp(cmd
, "mul2")) {
689 fgets(buf
, 4095, stdin
);
690 mp_read_radix(&a
, buf
, 64);
691 fgets(buf
, 4095, stdin
);
692 mp_read_radix(&b
, buf
, 64);
694 if (mp_cmp(&c
, &b
) != MP_EQ
) {
695 printf("mul_2 %lu failure\n", mul2_n
);
701 } else if (!strcmp(cmd
, "add_d")) {
703 fgets(buf
, 4095, stdin
);
704 mp_read_radix(&a
, buf
, 64);
705 fgets(buf
, 4095, stdin
);
706 sscanf(buf
, "%d", &ix
);
707 fgets(buf
, 4095, stdin
);
708 mp_read_radix(&b
, buf
, 64);
709 mp_add_d(&a
, ix
, &c
);
710 if (mp_cmp(&b
, &c
) != MP_EQ
) {
711 printf("add_d %lu failure\n", add_d_n
);
715 printf("d == %d\n", ix
);
718 } else if (!strcmp(cmd
, "sub_d")) {
720 fgets(buf
, 4095, stdin
);
721 mp_read_radix(&a
, buf
, 64);
722 fgets(buf
, 4095, stdin
);
723 sscanf(buf
, "%d", &ix
);
724 fgets(buf
, 4095, stdin
);
725 mp_read_radix(&b
, buf
, 64);
726 mp_sub_d(&a
, ix
, &c
);
727 if (mp_cmp(&b
, &c
) != MP_EQ
) {
728 printf("sub_d %lu failure\n", sub_d_n
);
732 printf("d == %d\n", ix
);
740 /* Source: /cvs/libtom/libtommath/demo/demo.c,v */
742 /* Date: 2005/06/24 11:32:07 */