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