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