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 11 srand(long seed) 12 { 13 randomreg = seed; 14 } 15 16 long 17 lrand(void) 18 { 19 randomreg = randomreg*104381 + 81761; 20 return randomreg; 21 } 22 23 void 24 prng(uchar *p, int n) 25 { 26 while(n-- > 0) 27 *p++ = lrand(); 28 } 29 30 void 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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