sys/src/libmp/test.c

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