1 #include <lib9.h>
2 #include <mp.h>
3 #include <libsec.h>
4
5 char *sfactors[] =
6 { "3", "5", "17", "257", "641", "65537", "274177", "2424833", "6700417", "45592577",
7 "6487031809", "67280421310721", "1238926361552897", "59649589127497217",
8 "5704689200685129054721", "4659775785220018543264560743076778192897",
9 "7455602825647884208337395736200454918783366342657",
10 "93461639715357977769163558199606896584051237541638188580280321",
11
12 "741640062627530801524787141901937474059940781097519023905821316144415759504705008092818711693940737",
13
14 "130439874405488189727484768796509903946608530841611892186895295776832416251471863574140227977573104895898783928842923844831149032913798729088601617946094119449010595906710130531906171018354491609619193912488538116080712299672322806217820753127014424577"
15 };
16
17 long start;
18
19 void
printmp(mpint * b,char * tag)20 printmp(mpint *b, char *tag)
21 {
22 int n;
23 char *p;
24
25 print("%s (%d) ", tag, b->top);
26 p = mptoa(b, 10, nil, 0);
27 write(1, p, strlen(p));
28 free(p);
29 print("\n");
30 }
31
32 int
timing(void)33 timing(void)
34 {
35 long now, span;
36
37 now = time(0);
38 span = now-start;
39 start = now;
40
41 return span;
42 }
43
44 int expdebug;
45
46
47 void
main(int argc,char ** argv)48 main(int argc, char **argv)
49 {
50 mpint *p, *k, *d, *b, *e, *x, *r;
51 int i;
52
53 start = time(0);
54 fmtinstall('B', mpconv);
55 mpsetminbits(2*Dbits);
56
57 x = mpnew(0);
58 e = mpnew(0);
59 r = mpnew(0);
60 p = mpnew(0);
61
62 // b = 2^32
63 b = mpcopy(mpone);
64 mpleft(b, 32, b);
65
66 // 2^29440
67 p = mpcopy(mpone);
68 mpleft(p, 29440, p);
69 // 2^27392
70 k = mpcopy(mpone);
71 mpleft(k, 27392, k);
72 // k = 2^29440 - 2^27392
73 mpsub(p, k, k);
74 // p = 2^29440 - 2^27392 + 1
75 mpadd(k, mpone, p);
76
77 // if(!probably_prime(p, 18)){
78 // print("not a prime\n");
79 // exits(0);
80 // }
81 // print("probably prime\n");
82
83 mpright(k, 10, k);
84 printmp(k, "k =");
85
86 expdebug = 1;
87 mpexp(b, k, p, x);
88 printmp(x, "x =");
89 print("timing %d\n", timing());
90
91 for(i = 0; i < nelem(sfactors); i++){
92 d = strtomp(sfactors[i], nil, 10, nil);
93 // e = k/d
94 mpdiv(k, d, e, r);
95 printmp(r, "r =");
96
97 // x = b^e mod p
98 mpexp(b, e, p, x);
99 printmp(x, "x =");
100 print("timing %d\n", timing());
101 }
102 }
103