1 #include <u.h>
2 #include <libc.h>
3 #include <mp.h>
4 #include "dat.h"
5
6 int loops = 1;
7
8 long randomreg;
9
10 void
srand(long seed)11 srand(long seed)
12 {
13 randomreg = seed;
14 }
15
16 long
lrand(void)17 lrand(void)
18 {
19 randomreg = randomreg*104381 + 81761;
20 return randomreg;
21 }
22
23 void
prng(uchar * p,int n)24 prng(uchar *p, int n)
25 {
26 while(n-- > 0)
27 *p++ = lrand();
28 }
29
30 void
testconv(char * str)31 testconv(char *str)
32 {
33 mpint *b;
34 char *p;
35
36 b = strtomp(str, nil, 16, nil);
37
38 p = mptoa(b, 10, nil, 0);
39 print("%s = ", p);
40 strtomp(p, nil, 10, b);
41 free(p);
42 print("%B\n", b);
43
44 p = mptoa(b, 16, nil, 0);
45 print("%s = ", p);
46 strtomp(p, nil, 16, b);
47 free(p);
48 print("%B\n", b);
49
50 p = mptoa(b, 32, nil, 0);
51 print("%s = ", p);
52 strtomp(p, nil, 32, b);
53 free(p);
54 print("%B\n", b);
55
56 p = mptoa(b, 64, nil, 0);
57 print("%s = ", p);
58 strtomp(p, nil, 64, b);
59 free(p);
60 print("%B\n", b);
61
62 mpfree(b);
63 }
64
65 void
testshift(char * str)66 testshift(char *str)
67 {
68 mpint *b1, *b2;
69 int i;
70
71 b1 = strtomp(str, nil, 16, nil);
72 b2 = mpnew(0);
73 for(i = 0; i < 64; i++){
74 mpleft(b1, i, b2);
75 print("%2.2d %B\n", i, b2);
76 }
77 for(i = 0; i < 64; i++){
78 mpright(b2, i, b1);
79 print("%2.2d %B\n", i, b1);
80 }
81 mpfree(b1);
82 mpfree(b2);
83 }
84
85 void
testaddsub(char * str)86 testaddsub(char *str)
87 {
88 mpint *b1, *b2;
89 int i;
90
91 b1 = strtomp(str, nil, 16, nil);
92 b2 = mpnew(0);
93 for(i = 0; i < 16; i++){
94 mpadd(b1, b2, b2);
95 print("%2.2d %B\n", i, b2);
96 }
97 for(i = 0; i < 16; i++){
98 mpsub(b2, b1, b2);
99 print("%2.2d %B\n", i, b2);
100 }
101 mpfree(b1);
102 mpfree(b2);
103 }
104
105 void
testvecdigmuladd(char * str,mpdigit d)106 testvecdigmuladd(char *str, mpdigit d)
107 {
108 mpint *b, *b2;
109 int i;
110 vlong now;
111
112 b = strtomp(str, nil, 16, nil);
113 b2 = mpnew(0);
114
115 mpbits(b2, (b->top+1)*Dbits);
116 now = nsec();
117 for(i = 0; i < loops; i++){
118 memset(b2->p, 0, b2->top*Dbytes);
119 mpvecdigmuladd(b->p, b->top, d, b2->p);
120 }
121 if(loops > 1)
122 print("%lld ns for a %d*%d vecdigmul\n", (nsec()-now)/loops, b->top*Dbits, Dbits);
123 mpnorm(b2);
124 print("0 + %B * %ux = %B\n", b, d, b2);
125
126 mpfree(b);
127 mpfree(b2);
128 }
129
130 void
testvecdigmulsub(char * str,mpdigit d)131 testvecdigmulsub(char *str, mpdigit d)
132 {
133 mpint *b, *b2;
134 int i;
135 vlong now;
136
137 b = strtomp(str, nil, 16, nil);
138 b2 = mpnew(0);
139
140 mpbits(b2, (b->top+1)*Dbits);
141 now = nsec();
142 for(i = 0; i < loops; i++){
143 memset(b2->p, 0, b2->top*Dbytes);
144 mpvecdigmulsub(b->p, b->top, d, b2->p);
145 }
146 if(loops > 1)
147 print("%lld ns for a %d*%d vecdigmul\n", (nsec()-now)/loops, b->top*Dbits, Dbits);
148 mpnorm(b2);
149 print("0 - %B * %ux = %B\n", b, d, b2);
150
151 mpfree(b);
152 mpfree(b2);
153 }
154
155 void
testmul(char * str)156 testmul(char *str)
157 {
158 mpint *b, *b1, *b2;
159 vlong now;
160 int i;
161
162 b = strtomp(str, nil, 16, nil);
163 b1 = mpcopy(b);
164 b2 = mpnew(0);
165
166 now = nsec();
167 for(i = 0; i < loops; i++)
168 mpmul(b, b1, b2);
169 if(loops > 1)
170 print("%lld µs for a %d*%d mult\n", (nsec()-now)/(loops*1000),
171 b->top*Dbits, b1->top*Dbits);
172 print("%B * %B = %B\n", b, b1, b2);
173
174 mpfree(b);
175 mpfree(b1);
176 mpfree(b2);
177 }
178
179 void
testmul2(mpint * b,mpint * b1)180 testmul2(mpint *b, mpint *b1)
181 {
182 mpint *b2;
183 vlong now;
184 int i;
185
186 b2 = mpnew(0);
187
188 now = nsec();
189 for(i = 0; i < loops; i++)
190 mpmul(b, b1, b2);
191 if(loops > 1)
192 print("%lld µs for a %d*%d mult\n", (nsec()-now)/(loops*1000), b->top*Dbits, b1->top*Dbits);
193 print("%B * ", b);
194 print("%B = ", b1);
195 print("%B\n", b2);
196
197 mpfree(b2);
198 }
199
200 void
testdigdiv(char * str,mpdigit d)201 testdigdiv(char *str, mpdigit d)
202 {
203 mpint *b;
204 mpdigit q;
205 int i;
206 vlong now;
207
208 b = strtomp(str, nil, 16, nil);
209 now = nsec();
210 for(i = 0; i < loops; i++)
211 mpdigdiv(b->p, d, &q);
212 if(loops > 1)
213 print("%lld ns for a %d / %d div\n", (nsec()-now)/loops, 2*Dbits, Dbits);
214 print("%B / %ux = %ux\n", b, d, q);
215 mpfree(b);
216 }
217
218 void
testdiv(mpint * x,mpint * y)219 testdiv(mpint *x, mpint *y)
220 {
221 mpint *b2, *b3;
222 vlong now;
223 int i;
224
225 b2 = mpnew(0);
226 b3 = mpnew(0);
227 now = nsec();
228 for(i = 0; i < loops; i++)
229 mpdiv(x, y, b2, b3);
230 if(loops > 1)
231 print("%lld µs for a %d/%d div\n", (nsec()-now)/(1000*loops),
232 x->top*Dbits, y->top*Dbits);
233 print("%B / %B = %B %B\n", x, y, b2, b3);
234 mpfree(b2);
235 mpfree(b3);
236 }
237
238 void
testmod(mpint * x,mpint * y)239 testmod(mpint *x, mpint *y)
240 {
241 mpint *r;
242 vlong now;
243 int i;
244
245 r = mpnew(0);
246 now = nsec();
247 for(i = 0; i < loops; i++)
248 mpmod(x, y, r);
249 if(loops > 1)
250 print("%lld µs for a %d/%d mod\n", (nsec()-now)/(1000*loops),
251 x->top*Dbits, y->top*Dbits);
252 print("%B mod %B = %B\n", x, y, r);
253 mpfree(r);
254 }
255
256 void
testinvert(mpint * x,mpint * y)257 testinvert(mpint *x, mpint *y)
258 {
259 mpint *r, *d1, *d2;
260 vlong now;
261 int i;
262
263 r = mpnew(0);
264 d1 = mpnew(0);
265 d2 = mpnew(0);
266 now = nsec();
267 mpextendedgcd(x, y, r, d1, d2);
268 mpdiv(x, r, x, d1);
269 mpdiv(y, r, y, d1);
270 for(i = 0; i < loops; i++)
271 mpinvert(x, y, r);
272 if(loops > 1)
273 print("%lld µs for a %d in %d invert\n", (nsec()-now)/(1000*loops),
274 x->top*Dbits, y->top*Dbits);
275 print("%B**-1 mod %B = %B\n", x, y, r);
276 mpmul(r, x, d1);
277 mpmod(d1, y, d2);
278 print("%B*%B mod %B = %B\n", x, r, y, d2);
279 mpfree(r);
280 mpfree(d1);
281 mpfree(d2);
282 }
283
284 void
testsub1(char * a,char * b)285 testsub1(char *a, char *b)
286 {
287 mpint *b1, *b2, *b3;
288
289 b1 = strtomp(a, nil, 16, nil);
290 b2 = strtomp(b, nil, 16, nil);
291 b3 = mpnew(0);
292 mpsub(b1, b2, b3);
293 print("%B - %B = %B\n", b1, b2, b3);
294 }
295
296 void
testmul1(char * a,char * b)297 testmul1(char *a, char *b)
298 {
299 mpint *b1, *b2, *b3;
300
301 b1 = strtomp(a, nil, 16, nil);
302 b2 = strtomp(b, nil, 16, nil);
303 b3 = mpnew(0);
304 mpmul(b1, b2, b3);
305 print("%B * %B = %B\n", b1, b2, b3);
306 }
307
308 void
testexp(char * base,char * exp,char * mod)309 testexp(char *base, char *exp, char *mod)
310 {
311 mpint *b, *e, *m, *res;
312 int i;
313 uvlong now;
314
315 b = strtomp(base, nil, 16, nil);
316 e = strtomp(exp, nil, 16, nil);
317 res = mpnew(0);
318 if(mod != nil)
319 m = strtomp(mod, nil, 16, nil);
320 else
321 m = nil;
322 now = nsec();
323 for(i = 0; i < loops; i++)
324 mpexp(b, e, m, res);
325 if(loops > 1)
326 print("%ulldµs for a %d to the %d bit exp\n", (nsec()-now)/(loops*1000),
327 b->top*Dbits, e->top*Dbits);
328 if(m != nil)
329 print("%B ^ %B mod %B == %B\n", b, e, m, res);
330 else
331 print("%B ^ %B == %B\n", b, e, res);
332 mpfree(b);
333 mpfree(e);
334 mpfree(res);
335 if(m != nil)
336 mpfree(m);
337 }
338
339 void
testgcd(void)340 testgcd(void)
341 {
342 mpint *a, *b, *d, *x, *y, *t1, *t2;
343 int i;
344 uvlong now, then;
345 uvlong etime;
346
347 d = mpnew(0);
348 x = mpnew(0);
349 y = mpnew(0);
350 t1 = mpnew(0);
351 t2 = mpnew(0);
352
353 etime = 0;
354
355 a = strtomp("4EECAB3E04C4E6BC1F49D438731450396BF272B4D7B08F91C38E88ADCD281699889AFF872E2204C80CCAA8E460797103DE539D5DF8335A9B20C0B44886384F134C517287202FCA914D8A5096446B40CD861C641EF9C2730CB057D7B133F4C2B16BBD3D75FDDBD9151AAF0F9144AAA473AC93CF945DBFE0859FB685D5CBD0A8B3", nil, 16, nil);
356 b = strtomp("C41CFBE4D4846F67A3DF7DE9921A49D3B42DC33728427AB159CEC8CBBDB12B5F0C244F1A734AEB9840804EA3C25036AD1B61AFF3ABBC247CD4B384224567A863A6F020E7EE9795554BCD08ABAD7321AF27E1E92E3DB1C6E7E94FAAE590AE9C48F96D93D178E809401ABE8A534A1EC44359733475A36A70C7B425125062B1142D", nil, 16, nil);
357 mpextendedgcd(a, b, d, x, y);
358 print("gcd %B*%B+%B*%B = %B?\n", a, x, b, y, d);
359 mpfree(a);
360 mpfree(b);
361
362 for(i = 0; i < loops; i++){
363 a = mprand(2048, prng, nil);
364 b = mprand(2048, prng, nil);
365 then = nsec();
366 mpextendedgcd(a, b, d, x, y);
367 now = nsec();
368 etime += now-then;
369 mpmul(a, x, t1);
370 mpmul(b, y, t2);
371 mpadd(t1, t2, t2);
372 if(mpcmp(d, t2) != 0)
373 print("%d gcd %B*%B+%B*%B != %B\n", i, a, x, b, y, d);
374 // else
375 // print("%d euclid %B*%B+%B*%B == %B\n", i, a, x, b, y, d);
376 mpfree(a);
377 mpfree(b);
378 }
379
380 mpfree(x);
381 mpfree(y);
382 mpfree(d);
383 mpfree(t1);
384 mpfree(t2);
385
386 if(loops > 1)
387 print("binary %llud\n", etime);
388 }
389
390 extern int _unnormalizedwarning = 1;
391
392 void
main(int argc,char ** argv)393 main(int argc, char **argv)
394 {
395 mpint *x, *y;
396
397 ARGBEGIN{
398 case 'n':
399 loops = atoi(ARGF());
400 break;
401 }ARGEND;
402
403 fmtinstall('B', mpconv);
404 fmtinstall('Q', mpconv);
405 srand(0);
406 mpsetminbits(2*Dbits);
407 testshift("1111111111111111");
408 testaddsub("fffffffffffffffff");
409 testdigdiv("1234567812345678", 0x76543218);
410 testdigdiv("1ffff", 0xffff);
411 testdigdiv("ffff", 0xffff);
412 testdigdiv("fff", 0xffff);
413 testdigdiv("effffffff", 0xffff);
414 testdigdiv("ffffffff", 0x1);
415 testdigdiv("ffffffff", 0);
416 testdigdiv("200000000", 2);
417 testdigdiv("ffffff00fffffff1", 0xfffffff1);
418 testvecdigmuladd("fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff", 2);
419 testconv("0");
420 testconv("-abc0123456789abcedf");
421 testconv("abc0123456789abcedf");
422 testvecdigmulsub("fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff", 2);
423 testsub1("1FFFFFFFE00000000", "FFFFFFFE00000001");
424 testmul1("ffffffff", "f");
425 testmul("ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff");
426 testmul1("100000000000000000000000000000000000000000000000000000002000000000000000000000000000000000000000000000000000000030000000000000000000000000000000000000000000000000000000400000000000000000000000000000000000000000000000000000004FFFFFFFFFFFFFFFE0000000200000000000000000000000000000003FFFFFFFFFFFFFFFE0000000200000000000000000000000000000002FFFFFFFFFFFFFFFE0000000200000000000000000000000000000001FFFFFFFFFFFFFFFE0000000200000000000000000000000000000000FFFFFFFFFFFFFFFE0000000200000000FFFFFFFE00000001", "FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF");
427 testmul1("1000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001", "1000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001");
428 testmul1("1000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001", "1000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001");
429 testmul1("FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000", "FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000");
430 x = mprand(256, prng, nil);
431 y = mprand(128, prng, nil);
432 testdiv(x, y);
433 x = mprand(2048, prng, nil);
434 y = mprand(1024, prng, nil);
435 testdiv(x, y);
436 // x = mprand(4*1024, prng, nil);
437 // y = mprand(4*1024, prng, nil);
438 // testmul2(x, y);
439 testsub1("677132C9", "-A26559B6");
440 testgcd();
441 x = mprand(512, prng, nil);
442 x->sign = -1;
443 y = mprand(256, prng, nil);
444 testdiv(x, y);
445 testmod(x, y);
446 x->sign = 1;
447 testinvert(y, x);
448 testexp("111111111", "222", "1000000000000000000000");
449 testexp("ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff", "ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff", "100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000");
450 #ifdef asdf
451 #endif adsf
452 print("done\n");
453 exits(0);
454 }
455