libhebimath.h (11909B)
1 #include <hebimath.h> 2 3 #include <setjmp.h> 4 #include <stddef.h> 5 #include <stdint.h> 6 #include <stdio.h> 7 #include <stdlib.h> 8 9 #define BIGINT_LIBRARY "hebimath" 10 #define HEBIMATH 11 12 typedef hebi_int z_t[1]; 13 14 static z_t _0, _1, _2, _add_sub_a, _add_sub_b, _cmp_b, _pow_a, _pow_m; 15 static z_t _trunc_btest_a, _bits_lsb_a, _str_a, _str_b, _str_m, _bset_a; 16 static z_t _logic_a, _logic_b, _logic_x, _not_a, _not_b, _gcd_a, _gcd_b; 17 static z_t _shift_b, _div_a, _div_b, _divmod; 18 static int error; 19 static jmp_buf jbuf; 20 21 #define FAST_RANDOM 0 22 #define SECURE_RANDOM 0 23 #define DEFAULT_RANDOM 0 24 #define FASTEST_RANDOM 0 25 #define LIBC_RAND_RANDOM 0 26 #define LIBC_RANDOM_RANDOM 0 27 #define LIBC_RAND48_RANDOM 0 28 #define QUASIUNIFORM 0 29 #define UNIFORM 1 30 #define MODUNIFORM 2 31 32 #ifdef ZAHL_UNSAFE 33 # define try(expr) (expr) 34 #else 35 # define try(expr) do if ((error = (expr))) longjmp(jbuf, 1); while (0) 36 #endif 37 38 #define zsave(a, s) zstr(a, s, sizeof(s) - 1) 39 #define zload(a, s) zsets(a, s) 40 41 static inline void 42 zsetup(jmp_buf env) 43 { 44 *jbuf = *env; 45 hebi_init(_0); 46 hebi_init(_1); 47 hebi_init(_2); 48 hebi_init(_add_sub_a); 49 hebi_init(_add_sub_b); 50 hebi_init(_cmp_b); 51 hebi_init(_pow_a); 52 hebi_init(_pow_m); 53 hebi_init(_trunc_btest_a); 54 hebi_init(_bits_lsb_a); 55 hebi_init(_str_a); 56 hebi_init(_str_b); 57 hebi_init(_str_m); 58 hebi_init(_bset_a); 59 hebi_init(_logic_a); 60 hebi_init(_logic_b); 61 hebi_init(_logic_x); 62 hebi_init(_not_a); 63 hebi_init(_not_b); 64 hebi_init(_gcd_a); 65 hebi_init(_gcd_b); 66 hebi_init(_shift_b); 67 hebi_init(_div_a); 68 hebi_init(_div_b); 69 hebi_init(_divmod); 70 try(hebi_set_u(_0, 0)); 71 try(hebi_set_u(_1, 1)); 72 try(hebi_set_u(_2, 2)); 73 } 74 75 static inline void 76 zunsetup(void) 77 { 78 hebi_destroy(_0); 79 hebi_destroy(_1); 80 hebi_destroy(_2); 81 hebi_destroy(_add_sub_a); 82 hebi_destroy(_add_sub_b); 83 hebi_destroy(_cmp_b); 84 hebi_destroy(_pow_a); 85 hebi_destroy(_pow_m); 86 hebi_destroy(_trunc_btest_a); 87 hebi_destroy(_bits_lsb_a); 88 hebi_destroy(_str_a); 89 hebi_destroy(_str_b); 90 hebi_destroy(_str_m); 91 hebi_destroy(_bset_a); 92 hebi_destroy(_logic_a); 93 hebi_destroy(_logic_b); 94 hebi_destroy(_logic_x); 95 hebi_destroy(_not_a); 96 hebi_destroy(_not_b); 97 hebi_destroy(_gcd_a); 98 hebi_destroy(_gcd_b); 99 hebi_destroy(_shift_b); 100 hebi_destroy(_div_a); 101 hebi_destroy(_div_b); 102 hebi_destroy(_divmod); 103 } 104 105 static inline void 106 zperror(const char *str) 107 { 108 const char *serr; 109 switch (error) { 110 case hebi_badrange: serr = "hebi_badrange"; break; 111 case hebi_badvalue: serr = "hebi_badvalue"; break; 112 case hebi_nomem: serr = "hebi_nomem"; break; 113 default: serr = "unknown error"; break; 114 } 115 if (str && *str) 116 fprintf(stderr, "%s: %s\n", str, serr); 117 else 118 fprintf(stderr, "%s\n", serr); 119 } 120 121 static inline void 122 zinit(z_t a) 123 { 124 hebi_init(a); 125 } 126 127 static inline void 128 zfree(z_t a) 129 { 130 hebi_destroy(a); 131 } 132 133 static inline void 134 zset(z_t r, const z_t a) 135 { 136 try(hebi_set(r, a)); 137 } 138 139 static inline void 140 zsetu(z_t r, unsigned long long int a) 141 { 142 try(hebi_set_u(r, a)); 143 } 144 145 static inline void 146 zseti(z_t r, long long int a) 147 { 148 try(hebi_set_i(r, a)); 149 } 150 151 static inline void 152 zneg(z_t r, const z_t a) 153 { 154 try(hebi_neg(r, a)); 155 } 156 157 static inline void 158 zabs(z_t r, const z_t a) 159 { 160 if (hebi_sign(a) < 0) 161 zneg(r, a); 162 else 163 zset(r, a); 164 } 165 166 static inline void 167 zadd(z_t r, const z_t a, const z_t b) 168 { 169 try(hebi_add(r, a, b)); 170 } 171 172 static inline void 173 zsub(z_t r, const z_t a, const z_t b) 174 { 175 try(hebi_sub(r, a, b)); 176 } 177 178 static inline void 179 zadd_unsigned(z_t r, const z_t a, const z_t b) 180 { 181 zabs(_add_sub_a, a); 182 zabs(_add_sub_b, b); 183 zadd(r, _add_sub_a, _add_sub_b); 184 } 185 186 static inline void 187 zsub_unsigned(z_t r, const z_t a, const z_t b) 188 { 189 zabs(_add_sub_a, a); 190 zabs(_add_sub_b, b); 191 zsub(r, _add_sub_a, _add_sub_b); 192 } 193 194 static inline int 195 zeven(const z_t a) 196 { 197 return ~hebi_get_u(a) & 1; 198 } 199 200 static inline int 201 zodd(const z_t a) 202 { 203 return hebi_get_u(a) & 1; 204 } 205 206 static inline int 207 zeven_nonzero(const z_t a) 208 { 209 return zeven(a); 210 } 211 212 static inline int 213 zodd_nonzero(const z_t a) 214 { 215 return zodd(a); 216 } 217 218 static inline void 219 zswap(z_t a, z_t b) 220 { 221 hebi_swap(a, b); 222 } 223 224 static inline int 225 zcmpmag(const z_t a, const z_t b) 226 { 227 return hebi_cmp_mag(a, b); 228 } 229 230 static inline int 231 zcmp(const z_t a, const z_t b) 232 { 233 return hebi_cmp(a, b); 234 } 235 236 static inline int 237 zcmpi(const z_t a, long long int b) 238 { 239 zseti(_cmp_b, b); 240 return zcmp(a, _cmp_b); 241 } 242 243 static inline int 244 zcmpu(const z_t a, long long int b) 245 { 246 zsetu(_cmp_b, b); 247 return zcmp(a, _cmp_b); 248 } 249 250 static inline int 251 zsignum(const z_t a) 252 { 253 return zcmp(a, _0); 254 } 255 256 static inline int 257 zzero(const z_t a) 258 { 259 return !zsignum(a); 260 } 261 262 static inline void 263 zmul(z_t r, const z_t a, const z_t b) 264 { 265 try(hebi_mul(r, a, b)); 266 } 267 268 static inline void 269 zsqr(z_t r, const z_t a) 270 { 271 zmul(r, a, a); 272 } 273 274 static inline void 275 zsets(z_t a, const char *s) 276 { 277 try(hebi_set_str(a, s, 0, 10)); 278 } 279 280 static inline void 281 zdivmod(z_t q, z_t r, const z_t a, const z_t b) 282 { 283 try(hebi_div(q, r, a, b)); 284 } 285 286 static inline void 287 zdiv(z_t q, const z_t a, const z_t b) 288 { 289 zdivmod(q, 0, a, b); 290 } 291 292 static inline void 293 zmod(z_t r, const z_t a, const z_t b) 294 { 295 zdivmod(_divmod, r, a, b); 296 } 297 298 static inline void 299 zmodmul(z_t r, const z_t a, const z_t b, const z_t m) 300 { 301 zmul(r, a, b); 302 zmod(r, r, m); 303 } 304 305 static inline void 306 zmodsqr(z_t r, const z_t a, const z_t m) 307 { 308 zsqr(r, a); 309 zmod(r, r, m); 310 } 311 312 static inline void 313 zpowu(z_t r, const z_t a, unsigned long long int b) 314 { 315 if (!b) { 316 if (zzero(a)) { 317 error = hebi_badvalue; 318 longjmp(jbuf, 1); 319 } 320 zsetu(r, 1); 321 return; 322 } 323 324 zset(_pow_a, a); 325 zsetu(r, 1); 326 for (; b; b >>= 1) { 327 if (b & 1) 328 zmul(r, r, _pow_a); 329 zsqr(_pow_a, _pow_a); 330 } 331 } 332 333 static inline void 334 zpow(z_t r, const z_t a, const z_t b) 335 { 336 zpowu(r, a, hebi_get_u(b)); 337 } 338 339 static inline void 340 zmodpowu(z_t r, const z_t a, unsigned long long int b, const z_t m) 341 { 342 if (!b) { 343 if (zzero(a) || zzero(m)) { 344 error = hebi_badvalue; 345 longjmp(jbuf, 1); 346 } 347 zsetu(r, 1); 348 return; 349 } else if (zzero(m)) { 350 error = hebi_badvalue; 351 longjmp(jbuf, 1); 352 } 353 354 zmod(_pow_a, a, m); 355 zset(_pow_m, m); 356 zsetu(r, 1); 357 for (; b; b >>= 1) { 358 if (b & 1) 359 zmodmul(r, r, _pow_a, _pow_m); 360 zmodsqr(_pow_a, _pow_a, _pow_m); 361 } 362 } 363 364 static inline void 365 zmodpow(z_t r, const z_t a, const z_t b, const z_t m) 366 { 367 zmodpowu(r, a, hebi_get_u(b), m); 368 } 369 370 static inline void 371 zlsh(z_t r, const z_t a, size_t b) 372 { 373 try(hebi_shl(r, a, b)); 374 } 375 376 static inline void 377 zrsh(z_t r, const z_t a, size_t b) 378 { 379 try(hebi_shr(r, a, b)); 380 } 381 382 static inline void 383 ztrunc(z_t r, const z_t a, size_t b) 384 { 385 zrsh(_trunc_btest_a, a, b); 386 zlsh(_trunc_btest_a, _trunc_btest_a, b); 387 zsub(r, a, _trunc_btest_a); 388 } 389 390 static inline int 391 zbtest(z_t a, size_t b) 392 { 393 zlsh(_trunc_btest_a, a, b); 394 return zodd(a); 395 } 396 397 static inline size_t 398 zbits(const z_t a) 399 { 400 hebi_uword x = x; 401 size_t rc = 0; 402 zset(_bits_lsb_a, a); 403 while (zsignum(_bits_lsb_a)) { 404 x = hebi_get_u(_bits_lsb_a); 405 zrsh(_bits_lsb_a, _bits_lsb_a, 8 * sizeof(x)); 406 rc += 8 * sizeof(x); 407 } 408 if (rc) 409 rc -= 8 * sizeof(x); 410 while (x) { 411 x >>= 1; 412 rc += 1; 413 } 414 return rc; 415 } 416 417 static inline size_t 418 zlsb(const z_t a) 419 { 420 hebi_uword x; 421 size_t rc = 0; 422 if (zzero(a)) 423 return SIZE_MAX; 424 zset(_bits_lsb_a, a); 425 while (!(x = hebi_get_u(_bits_lsb_a))) { 426 zrsh(_bits_lsb_a, _bits_lsb_a, 8 * sizeof(x)); 427 rc += 8 * sizeof(x); 428 } 429 while (~x & 1) { 430 x >>= 1; 431 rc += 1; 432 } 433 return rc; 434 } 435 436 static inline void 437 zptest(z_t w, const z_t a, int t) 438 { 439 static int gave_up = 0; 440 if (!gave_up) { 441 gave_up = 1; 442 printf("I'm sorry, primality test has not been implemented.\n"); 443 } 444 (void) w; 445 (void) a; 446 (void) t; 447 } 448 449 static inline size_t 450 zstr_length(const z_t a, unsigned long long int radix) 451 { 452 size_t size_total = 1, size_temp; 453 zset(_str_a, a); 454 while (!zzero(_str_a)) { 455 zsetu(_str_m, radix); 456 zset(_str_b, _str_m); 457 size_temp = 1; 458 while (zcmpmag(_str_m, _str_a) <= 0) { 459 zset(_str_b, _str_m); 460 zsqr(_str_m, _str_m); 461 size_temp <<= 1; 462 } 463 size_temp >>= 1; 464 size_total += size_temp; 465 zdiv(_str_a, _str_a, _str_b); 466 } 467 return size_total + (zsignum(a) < 0); 468 } 469 470 static inline char * 471 zstr(const z_t a, char *s, size_t n) 472 { 473 if (!n) 474 n = zstr_length(a, 10) * 2 + 1; 475 try(hebi_get_str(s, n, a, 10)); 476 return s; 477 } 478 479 static inline void 480 zsplit(z_t high, z_t low, const z_t a, size_t brk) 481 { 482 if (low == a) { 483 zrsh(high, a, brk); 484 ztrunc(low, a, brk); 485 } else { 486 ztrunc(low, a, brk); 487 zrsh(high, a, brk); 488 } 489 } 490 491 static inline void 492 zbset(z_t r, const z_t a, size_t bit, int mode) 493 { 494 zrsh(_bset_a, a, bit); 495 if (mode && zeven(_bset_a)) { 496 zlsh(_bset_a, _1, bit); 497 zadd(r, a, _bset_a); 498 } else if (mode <= 0 && zodd(_bset_a)) { 499 zlsh(_bset_a, _1, bit); 500 zsub(r, a, _bset_a); 501 } else { 502 zset(r, a); 503 } 504 } 505 506 static inline void 507 zrand(z_t r, int dev, int dist, const z_t n) 508 { 509 static int gave_up[] = {0, 0, 0}; 510 if (!gave_up[dist]) { 511 gave_up[dist] = 1; 512 printf("I'm sorry, prng has not been implemented.\n"); 513 } 514 (void) r; 515 (void) dev; 516 (void) n; 517 } 518 519 static inline void 520 zand(z_t r, const z_t a, const z_t b) 521 { 522 int neg = hebi_sign(a) < 0 && hebi_sign(b) < 0; 523 hebi_uword x; 524 size_t i = 0; 525 zset(_logic_a, a); 526 zset(_logic_b, b); 527 zsetu(r, 0); 528 while (zsignum(_logic_a) && zsignum(_logic_b)) { 529 x = hebi_get_u(_logic_a) & hebi_get_u(_logic_b); 530 zsetu(_logic_x, x); 531 zlsh(_logic_x, _logic_x, i++ * 8 * sizeof(x)); 532 zadd(r, r, _logic_x); 533 zrsh(_logic_a, _logic_a, 8 * sizeof(x)); 534 zrsh(_logic_b, _logic_b, 8 * sizeof(x)); 535 } 536 if (neg) 537 zneg(r, r); 538 } 539 540 static inline void 541 zor(z_t r, const z_t a, const z_t b) 542 { 543 int neg = hebi_sign(a) < 0 || hebi_sign(b) < 0; 544 hebi_uword x; 545 size_t i = 0; 546 zset(_logic_a, a); 547 zset(_logic_b, b); 548 zsetu(r, 0); 549 while (zsignum(_logic_a) || zsignum(_logic_b)) { 550 x = hebi_get_u(_logic_a) | hebi_get_u(_logic_b); 551 zsetu(_logic_x, x); 552 zlsh(_logic_x, _logic_x, i++ * 8 * sizeof(x)); 553 zadd(r, r, _logic_x); 554 zrsh(_logic_a, _logic_a, 8 * sizeof(x)); 555 zrsh(_logic_b, _logic_b, 8 * sizeof(x)); 556 } 557 if (neg) 558 zneg(r, r); 559 } 560 561 static inline void 562 zxor(z_t r, const z_t a, const z_t b) 563 { 564 int neg = (hebi_sign(a) < 0) ^ (hebi_sign(b) < 0); 565 hebi_uword x; 566 size_t i = 0; 567 zset(_logic_a, a); 568 zset(_logic_b, b); 569 zsetu(r, 0); 570 while (zsignum(_logic_a) || zsignum(_logic_b)) { 571 x = hebi_get_u(_logic_a) ^ hebi_get_u(_logic_b); 572 zsetu(_logic_x, x); 573 zlsh(_logic_x, _logic_x, i++ * 8 * sizeof(x)); 574 zadd(r, r, _logic_x); 575 zrsh(_logic_a, _logic_a, 8 * sizeof(x)); 576 zrsh(_logic_b, _logic_b, 8 * sizeof(x)); 577 } 578 if (neg) 579 zneg(r, r); 580 } 581 582 static inline void 583 zgcd(z_t r, const z_t a, const z_t b) 584 { 585 size_t shifts, a_lsb, b_lsb; 586 int neg, cmpmag; 587 588 if (zzero(a)) { 589 if (r != b) 590 zset(r, b); 591 return; 592 } 593 if (zzero(b)) { 594 if (r != a) 595 zset(r, a); 596 return; 597 } 598 599 neg = hebi_sign(a) < 0 && hebi_sign(b) < 0; 600 601 a_lsb = zlsb(a); 602 b_lsb = zlsb(b); 603 shifts = a_lsb < b_lsb ? a_lsb : b_lsb; 604 zrsh(_gcd_a, a, a_lsb); 605 zrsh(_gcd_b, b, b_lsb); 606 607 for (;;) { 608 if ((cmpmag = zcmpmag(_gcd_a, _gcd_b)) >= 0) { 609 if (cmpmag == 0) 610 break; 611 zswap(_gcd_a, _gcd_b); 612 } 613 zsub(_gcd_b, _gcd_b, _gcd_a); 614 zrsh(_gcd_b, _gcd_b, zlsb(_gcd_b)); 615 } 616 617 zlsh(r, _gcd_a, shifts); 618 if (neg) 619 zneg(r, r); 620 } 621 622 static inline void 623 znot(z_t r, const z_t a) 624 { 625 size_t bits = zbits(a); 626 zsetu(_not_b, 0); 627 zsetu(_not_a, 1); 628 zlsh(_not_a, _not_a, bits); 629 zadd(_not_b, _not_b, _logic_a); 630 zsub(_not_b, _not_b, _1); 631 zxor(r, a, _not_b); 632 zneg(r, r); 633 } 634 635 /* Prototype declared, but implementation missing, in hebimath */ 636 637 int 638 hebi_shl(hebi_int *r, const hebi_int *a, unsigned int b) 639 { 640 zpowu(_shift_b, _2, b); 641 zmul(r, a, _shift_b); 642 return hebi_success; 643 } 644 645 int 646 hebi_shr(hebi_int *r, const hebi_int *a, unsigned int b) 647 { 648 zpowu(_shift_b, _2, b); 649 zdiv(r, a, _shift_b); 650 return hebi_success; 651 } 652 653 int 654 hebi_div(hebi_int *q, hebi_int *r, const hebi_int *a, const hebi_int *b) 655 { 656 int neg = zsignum(a) < 0; 657 zset(q, _0); 658 zabs(_div_a, a); 659 zabs(_div_b, b); 660 if (zzero(b)) { 661 error = hebi_badvalue; 662 longjmp(jbuf, 1); 663 } 664 while (zcmpmag(_div_a, _div_b) >= 0) { 665 zadd(q, q, _1); 666 zsub(_div_a, _div_a, _div_b); 667 } 668 if (neg) 669 zneg(q, q); 670 if (r) 671 zset(r, _div_a); 672 return hebi_success; 673 }