imath.c (83161B)
1 /* 2 ** imath.c 3 ** Arbitrary precision integer arithmetic routines. 4 ** See Copyright Notice in klisp.h 5 */ 6 7 /* 8 ** SOURCE NOTE: This is mostly from the IMath library, written by 9 ** M.J. Fromberger. It is adapted to klisp, mainly in the use of the 10 ** klisp allocator and fixing of digit size to 32 bits. 11 ** Imported from version (1.15) updated 01-Feb-2011 at 03:10 PM. 12 */ 13 14 #include "imath.h" 15 16 #if DEBUG 17 #include <stdio.h> 18 #endif 19 20 #include <stdlib.h> 21 #include <string.h> 22 #include <ctype.h> 23 24 #include <assert.h> 25 26 /* Andres Navarro: klisp includes */ 27 #include "kobject.h" 28 #include "kstate.h" 29 #include "kmem.h" 30 #include "kerror.h" 31 32 33 #if DEBUG 34 #define STATIC /* public */ 35 #else 36 #define STATIC static 37 #endif 38 39 /* {{{ Constants */ 40 41 const mp_result MP_OK = 0; /* no error, all is well */ 42 const mp_result MP_FALSE = 0; /* boolean false */ 43 const mp_result MP_TRUE = -1; /* boolean true */ 44 const mp_result MP_MEMORY = -2; /* out of memory */ 45 const mp_result MP_RANGE = -3; /* argument out of range */ 46 const mp_result MP_UNDEF = -4; /* result undefined */ 47 const mp_result MP_TRUNC = -5; /* output truncated */ 48 const mp_result MP_BADARG = -6; /* invalid null argument */ 49 const mp_result MP_MINERR = -6; 50 51 const mp_sign MP_NEG = 1; /* value is strictly negative */ 52 const mp_sign MP_ZPOS = 0; /* value is non-negative */ 53 54 STATIC const char *s_unknown_err = "unknown result code"; 55 STATIC const char *s_error_msg[] = { 56 "error code 0", 57 "boolean true", 58 "out of memory", 59 "argument out of range", 60 "result undefined", 61 "output truncated", 62 "invalid argument", 63 NULL 64 }; 65 66 /* }}} */ 67 68 /* Argument checking macros 69 Use CHECK() where a return value is required; NRCHECK() elsewhere */ 70 #define CHECK(TEST) assert(TEST) 71 #define NRCHECK(TEST) assert(TEST) 72 73 /* {{{ Logarithm table for computing output sizes */ 74 75 /* The ith entry of this table gives the value of log_i(2). 76 77 An integer value n requires ceil(log_i(n)) digits to be represented 78 in base i. Since it is easy to compute lg(n), by counting bits, we 79 can compute log_i(n) = lg(n) * log_i(2). 80 81 The use of this table eliminates a dependency upon linkage against 82 the standard math libraries. 83 */ 84 STATIC const double s_log2[] = { 85 0.000000000, 0.000000000, 1.000000000, 0.630929754, /* 0 1 2 3 */ 86 0.500000000, 0.430676558, 0.386852807, 0.356207187, /* 4 5 6 7 */ 87 0.333333333, 0.315464877, 0.301029996, 0.289064826, /* 8 9 10 11 */ 88 0.278942946, 0.270238154, 0.262649535, 0.255958025, /* 12 13 14 15 */ 89 0.250000000, 0.244650542, 0.239812467, 0.235408913, /* 16 17 18 19 */ 90 0.231378213, 0.227670249, 0.224243824, 0.221064729, /* 20 21 22 23 */ 91 0.218104292, 0.215338279, 0.212746054, 0.210309918, /* 24 25 26 27 */ 92 0.208014598, 0.205846832, 0.203795047, 0.201849087, /* 28 29 30 31 */ 93 0.200000000, 0.198239863, 0.196561632, 0.194959022, /* 32 33 34 35 */ 94 0.193426404, /* 36 */ 95 }; 96 97 /* }}} */ 98 /* {{{ Various macros */ 99 100 /* Return the number of digits needed to represent a static value */ 101 #define MP_VALUE_DIGITS(V) \ 102 ((sizeof(V)+(sizeof(mp_digit)-1))/sizeof(mp_digit)) 103 104 /* Round precision P to nearest word boundary */ 105 #define ROUND_PREC(P) ((mp_size)(2*(((P)+1)/2))) 106 107 /* Set array P of S digits to zero */ 108 #define ZERO(P, S) \ 109 do{ \ 110 mp_size i__ = (S) * sizeof(mp_digit); \ 111 mp_digit *p__ = (P); \ 112 memset(p__, 0, i__); \ 113 } while(0) 114 115 /* Copy S digits from array P to array Q */ 116 #define COPY(P, Q, S) \ 117 do{ \ 118 mp_size i__ = (S) * sizeof(mp_digit); \ 119 mp_digit *p__ = (P), *q__ = (Q); \ 120 memcpy(q__, p__, i__); \ 121 } while(0) 122 123 /* Reverse N elements of type T in array A */ 124 #define REV(T, A, N) \ 125 do{ \ 126 T *u_ = (A), *v_ = u_ + (N) - 1; \ 127 while (u_ < v_) { \ 128 T xch = *u_; \ 129 *u_++ = *v_; \ 130 *v_-- = xch; \ 131 } \ 132 } while(0) 133 134 #define CLAMP(Z) \ 135 do{ \ 136 mp_int z_ = (Z); \ 137 mp_size uz_ = MP_USED(z_); \ 138 mp_digit *dz_ = MP_DIGITS(z_) + uz_ -1; \ 139 while (uz_ > 1 && (*dz_-- == 0)) \ 140 --uz_; \ 141 MP_USED(z_) = uz_; \ 142 } while(0) 143 144 /* Select min/max. Do not provide expressions for which multiple 145 evaluation would be problematic, e.g. x++ */ 146 #define MIN(A, B) ((B)<(A)?(B):(A)) 147 #define MAX(A, B) ((B)>(A)?(B):(A)) 148 149 /* Exchange lvalues A and B of type T, e.g. 150 SWAP(int, x, y) where x and y are variables of type int. */ 151 #define SWAP(T, A, B) \ 152 do{ \ 153 T t_ = (A); \ 154 A = (B); \ 155 B = t_; \ 156 } while(0) 157 158 /* Used to set up and access simple temp stacks within functions. */ 159 #define TEMP(K) (temp + (K)) 160 #define SETUP(E, C) \ 161 do{ \ 162 if ((res = (E)) != MP_OK) \ 163 goto CLEANUP; \ 164 ++(C); \ 165 } while(0) 166 167 /* Compare value to zero. */ 168 #define CMPZ(Z) \ 169 (((Z)->used==1&&(Z)->digits[0]==0)?0:((Z)->sign==MP_NEG)?-1:1) 170 171 /* Multiply X by Y into Z, ignoring signs. Requires that Z have 172 enough storage preallocated to hold the result. */ 173 #define UMUL(K, X, Y, Z) \ 174 do{ \ 175 mp_size ua_ = MP_USED(X), ub_ = MP_USED(Y); \ 176 mp_size o_ = ua_ + ub_; \ 177 ZERO(MP_DIGITS(Z), o_); \ 178 (void) s_kmul((K), MP_DIGITS(X), MP_DIGITS(Y), MP_DIGITS(Z), ua_, ub_); \ 179 MP_USED(Z) = o_; \ 180 CLAMP(Z); \ 181 } while(0) 182 183 /* Square X into Z. Requires that Z have enough storage to hold the 184 result. */ 185 #define USQR(K, X, Z) \ 186 do{ \ 187 mp_size ua_ = MP_USED(X), o_ = ua_ + ua_; \ 188 ZERO(MP_DIGITS(Z), o_); \ 189 (void) s_ksqr((K), MP_DIGITS(X), MP_DIGITS(Z), ua_); \ 190 MP_USED(Z) = o_; \ 191 CLAMP(Z); \ 192 } while(0) 193 194 #define UPPER_HALF(W) ((mp_word)((W) >> MP_DIGIT_BIT)) 195 #define LOWER_HALF(W) ((mp_digit)(W)) 196 #define HIGH_BIT_SET(W) ((W) >> (MP_WORD_BIT - 1)) 197 #define ADD_WILL_OVERFLOW(W, V) ((MP_WORD_MAX - (V)) < (W)) 198 199 /* }}} */ 200 /* {{{ Default configuration settings */ 201 202 /* Default number of digits allocated to a new mp_int */ 203 #if IMATH_TEST 204 mp_size default_precision = MP_DEFAULT_PREC; 205 #else 206 STATIC const mp_size default_precision = MP_DEFAULT_PREC; 207 #endif 208 209 /* Minimum number of digits to invoke recursive multiply */ 210 #if IMATH_TEST 211 mp_size multiply_threshold = MP_MULT_THRESH; 212 #else 213 STATIC const mp_size multiply_threshold = MP_MULT_THRESH; 214 #endif 215 216 /* }}} */ 217 218 /* Allocate a buffer of (at least) num digits, or return 219 NULL if that couldn't be done. */ 220 STATIC mp_digit *s_alloc(klisp_State *K, mp_size num); 221 222 /* Release a buffer of digits allocated by s_alloc(). */ 223 STATIC void s_free(klisp_State *K, void *ptr, mp_size size); 224 225 /* Insure that z has at least min digits allocated, resizing if 226 necessary. Returns true if successful, false if out of memory. */ 227 STATIC int s_pad(klisp_State *K, mp_int z, mp_size min); 228 229 /* Fill in a "fake" mp_int on the stack with a given value */ 230 STATIC void s_fake(mp_int z, mp_small value, mp_digit vbuf[]); 231 232 /* Compare two runs of digits of given length, returns <0, 0, >0 */ 233 STATIC int s_cdig(mp_digit *da, mp_digit *db, mp_size len); 234 235 /* Pack the unsigned digits of v into array t */ 236 STATIC int s_vpack(mp_small v, mp_digit t[]); 237 238 /* Compare magnitudes of a and b, returns <0, 0, >0 */ 239 STATIC int s_ucmp(mp_int a, mp_int b); 240 241 /* Compare magnitudes of a and v, returns <0, 0, >0 */ 242 STATIC int s_vcmp(mp_int a, mp_small v); 243 244 /* Unsigned magnitude addition; assumes dc is big enough. 245 Carry out is returned (no memory allocated). */ 246 STATIC mp_digit s_uadd(mp_digit *da, mp_digit *db, mp_digit *dc, 247 mp_size size_a, mp_size size_b); 248 249 /* Unsigned magnitude subtraction. Assumes dc is big enough. */ 250 STATIC void s_usub(mp_digit *da, mp_digit *db, mp_digit *dc, 251 mp_size size_a, mp_size size_b); 252 253 /* Unsigned recursive multiplication. Assumes dc is big enough. */ 254 STATIC int s_kmul(klisp_State *K, mp_digit *da, mp_digit *db, 255 mp_digit *dc, mp_size size_a, mp_size size_b); 256 257 /* Unsigned magnitude multiplication. Assumes dc is big enough. */ 258 STATIC void s_umul(mp_digit *da, mp_digit *db, mp_digit *dc, 259 mp_size size_a, mp_size size_b); 260 261 /* Unsigned recursive squaring. Assumes dc is big enough. */ 262 /* Andres Navarro: allocs temporaries */ 263 STATIC int s_ksqr(klisp_State *K, mp_digit *da, mp_digit *dc, 264 mp_size size_a); 265 266 /* Unsigned magnitude squaring. Assumes dc is big enough. */ 267 STATIC void s_usqr(mp_digit *da, mp_digit *dc, mp_size size_a); 268 269 /* Single digit addition. Assumes a is big enough. */ 270 STATIC void s_dadd(mp_int a, mp_digit b); 271 272 /* Single digit multiplication. Assumes a is big enough. */ 273 STATIC void s_dmul(mp_int a, mp_digit b); 274 275 /* Single digit multiplication on buffers; assumes dc is big enough. */ 276 STATIC void s_dbmul(mp_digit *da, mp_digit b, mp_digit *dc, 277 mp_size size_a); 278 279 /* Single digit division. Replaces a with the quotient, 280 returns the remainder. */ 281 STATIC mp_digit s_ddiv(mp_int a, mp_digit b); 282 283 /* Quick division by a power of 2, replaces z (no allocation) */ 284 STATIC void s_qdiv(mp_int z, mp_size p2); 285 286 /* Quick remainder by a power of 2, replaces z (no allocation) */ 287 STATIC void s_qmod(mp_int z, mp_size p2); 288 289 /* Quick multiplication by a power of 2, replaces z. 290 Allocates if necessary; returns false in case this fails. */ 291 STATIC int s_qmul(klisp_State *K, mp_int z, mp_size p2); 292 293 /* Quick subtraction from a power of 2, replaces z. 294 Allocates if necessary; returns false in case this fails. */ 295 STATIC int s_qsub(klisp_State *K, mp_int z, mp_size p2); 296 297 /* Return maximum k such that 2^k divides z. */ 298 STATIC int s_dp2k(mp_int z); 299 300 /* Return k >= 0 such that z = 2^k, or -1 if there is no such k. */ 301 STATIC int s_isp2(mp_int z); 302 303 /* Set z to 2^k. May allocate; returns false in case this fails. */ 304 STATIC int s_2expt(klisp_State *K, mp_int z, mp_small k); 305 306 /* Normalize a and b for division, returns normalization constant */ 307 STATIC int s_norm(klisp_State *K, mp_int a, mp_int b); 308 309 /* Compute constant mu for Barrett reduction, given modulus m, result 310 replaces z, m is untouched. */ 311 STATIC mp_result s_brmu(klisp_State *K, mp_int z, mp_int m); 312 313 /* Reduce a modulo m, using Barrett's algorithm. */ 314 STATIC int s_reduce(klisp_State *K, mp_int x, mp_int m, mp_int mu, 315 mp_int q1, mp_int q2); 316 317 /* Modular exponentiation, using Barrett reduction */ 318 STATIC mp_result s_embar(klisp_State *K, mp_int a, mp_int b, mp_int m, 319 mp_int mu, mp_int c); 320 321 /* Unsigned magnitude division. Assumes |a| > |b|. Allocates 322 temporaries; overwrites a with quotient, b with remainder. */ 323 STATIC mp_result s_udiv(klisp_State *K, mp_int a, mp_int b); 324 325 /* Compute the number of digits in radix r required to represent the 326 given value. Does not account for sign flags, terminators, etc. */ 327 STATIC int s_outlen(mp_int z, mp_size r); 328 329 /* Guess how many digits of precision will be needed to represent a 330 radix r value of the specified number of digits. Returns a value 331 guaranteed to be no smaller than the actual number required. */ 332 STATIC mp_size s_inlen(int len, mp_size r); 333 334 /* Convert a character to a digit value in radix r, or 335 -1 if out of range */ 336 STATIC int s_ch2val(char c, int r); 337 338 /* Convert a digit value to a character */ 339 STATIC char s_val2ch(int v, int caps); 340 341 /* Take 2's complement of a buffer in place */ 342 STATIC void s_2comp(unsigned char *buf, int len); 343 344 /* Convert a value to binary, ignoring sign. On input, *limpos is the 345 bound on how many bytes should be written to buf; on output, *limpos 346 is set to the number of bytes actually written. */ 347 STATIC mp_result s_tobin(mp_int z, unsigned char *buf, int *limpos, int pad); 348 349 #if DEBUG 350 /* Dump a representation of the mp_int to standard output */ 351 void s_print(char *tag, mp_int z); 352 void s_print_buf(char *tag, mp_digit *buf, mp_size num); 353 #endif 354 355 /* {{{ mp_int_init(z) */ 356 357 mp_result mp_int_init(mp_int z) 358 { 359 if(z == NULL) 360 return MP_BADARG; 361 362 z->single = 0; 363 z->digits = &(z->single); 364 z->alloc = 1; 365 z->used = 1; 366 z->sign = MP_ZPOS; 367 368 return MP_OK; 369 } 370 371 /* }}} */ 372 373 /* {{{ mp_int_alloc() */ 374 375 mp_int mp_int_alloc(klisp_State *K) 376 { 377 mp_int out = klispM_new(K, mpz_t); 378 379 if(out != NULL) 380 mp_int_init(out); 381 382 return out; 383 } 384 385 /* }}} */ 386 387 /* {{{ mp_int_init_size(z, prec) */ 388 389 mp_result mp_int_init_size(klisp_State *K, mp_int z, mp_size prec) 390 { 391 CHECK(z != NULL); 392 393 if(prec == 0) 394 prec = default_precision; 395 else if(prec == 1) 396 return mp_int_init(z); 397 else 398 prec = (mp_size) ROUND_PREC(prec); 399 400 if((MP_DIGITS(z) = s_alloc(K, prec)) == NULL) 401 return MP_MEMORY; 402 403 z->digits[0] = 0; 404 MP_USED(z) = 1; 405 MP_ALLOC(z) = prec; 406 MP_SIGN(z) = MP_ZPOS; 407 408 return MP_OK; 409 } 410 411 /* }}} */ 412 413 /* {{{ mp_int_init_copy(z, old) */ 414 415 mp_result mp_int_init_copy(klisp_State *K, mp_int z, mp_int old) 416 { 417 mp_result res; 418 mp_size uold; 419 420 CHECK(z != NULL && old != NULL); 421 422 uold = MP_USED(old); 423 if(uold == 1) { 424 mp_int_init(z); 425 } 426 else { 427 mp_size target = MAX(uold, default_precision); 428 429 if((res = mp_int_init_size(K, z, target)) != MP_OK) 430 return res; 431 } 432 433 MP_USED(z) = uold; 434 MP_SIGN(z) = MP_SIGN(old); 435 COPY(MP_DIGITS(old), MP_DIGITS(z), uold); 436 437 return MP_OK; 438 } 439 440 /* }}} */ 441 442 /* {{{ mp_int_init_value(z, value) */ 443 444 mp_result mp_int_init_value(klisp_State *K, mp_int z, mp_small value) 445 { 446 mpz_t vtmp; 447 mp_digit vbuf[MP_VALUE_DIGITS(value)]; 448 449 s_fake(&vtmp, value, vbuf); 450 return mp_int_init_copy(K, z, &vtmp); 451 } 452 453 /* }}} */ 454 455 /* {{{ mp_int_set_value(z, value) */ 456 457 mp_result mp_int_set_value(klisp_State *K, mp_int z, mp_small value) 458 { 459 mpz_t vtmp; 460 mp_digit vbuf[MP_VALUE_DIGITS(value)]; 461 462 s_fake(&vtmp, value, vbuf); 463 return mp_int_copy(K, &vtmp, z); 464 } 465 466 /* }}} */ 467 468 /* {{{ mp_int_clear(z) */ 469 470 void mp_int_clear(klisp_State *K, mp_int z) 471 { 472 if(z == NULL) 473 return; 474 475 if(MP_DIGITS(z) != NULL) { 476 if((void *) MP_DIGITS(z) != (void *) &MP_SINGLE(z)) 477 s_free(K, MP_DIGITS(z), MP_ALLOC(z)); 478 MP_DIGITS(z) = NULL; 479 } 480 } 481 482 /* }}} */ 483 484 /* {{{ mp_int_free(z) */ 485 486 void mp_int_free(klisp_State *K, mp_int z) 487 { 488 NRCHECK(z != NULL); 489 490 mp_int_clear(K, z); 491 klispM_free(K, z); /* note: NOT s_free() */ 492 } 493 494 /* }}} */ 495 496 /* {{{ mp_int_copy(a, c) */ 497 498 mp_result mp_int_copy(klisp_State *K, mp_int a, mp_int c) 499 { 500 CHECK(a != NULL && c != NULL); 501 502 if(a != c) { 503 mp_size ua = MP_USED(a); 504 mp_digit *da, *dc; 505 506 if(!s_pad(K, c, ua)) 507 return MP_MEMORY; 508 509 da = MP_DIGITS(a); dc = MP_DIGITS(c); 510 COPY(da, dc, ua); 511 512 MP_USED(c) = ua; 513 MP_SIGN(c) = MP_SIGN(a); 514 } 515 516 return MP_OK; 517 } 518 519 /* }}} */ 520 521 /* {{{ mp_int_swap(a, c) */ 522 523 void mp_int_swap(mp_int a, mp_int c) 524 { 525 if(a != c) { 526 mpz_t tmp = *a; 527 528 *a = *c; 529 *c = tmp; 530 /* Andres Navarro: bugfix */ 531 /* correct if digits was pointing to single */ 532 if (a->digits == &c->single) 533 a->digits = &a->single; 534 if (c->digits == &a->single) 535 c->digits = &c->single; 536 /* Andres Navarro: /bugfix */ 537 } 538 } 539 540 /* }}} */ 541 542 /* {{{ mp_int_zero(z) */ 543 544 void mp_int_zero(mp_int z) 545 { 546 NRCHECK(z != NULL); 547 548 z->digits[0] = 0; 549 MP_USED(z) = 1; 550 MP_SIGN(z) = MP_ZPOS; 551 } 552 553 /* }}} */ 554 555 /* {{{ mp_int_abs(a, c) */ 556 557 mp_result mp_int_abs(klisp_State *K, mp_int a, mp_int c) 558 { 559 mp_result res; 560 561 CHECK(a != NULL && c != NULL); 562 563 if((res = mp_int_copy(K, a, c)) != MP_OK) 564 return res; 565 566 MP_SIGN(c) = MP_ZPOS; 567 return MP_OK; 568 } 569 570 /* }}} */ 571 572 /* {{{ mp_int_neg(a, c) */ 573 574 mp_result mp_int_neg(klisp_State *K, mp_int a, mp_int c) 575 { 576 mp_result res; 577 578 CHECK(a != NULL && c != NULL); 579 580 if((res = mp_int_copy(K, a, c)) != MP_OK) 581 return res; 582 583 if(CMPZ(c) != 0) 584 MP_SIGN(c) = 1 - MP_SIGN(a); 585 586 return MP_OK; 587 } 588 589 /* }}} */ 590 591 /* {{{ mp_int_add(a, b, c) */ 592 593 mp_result mp_int_add(klisp_State *K, mp_int a, mp_int b, mp_int c) 594 { 595 mp_size ua, ub, uc, max; 596 597 CHECK(a != NULL && b != NULL && c != NULL); 598 599 ua = MP_USED(a); ub = MP_USED(b); uc = MP_USED(c); 600 max = MAX(ua, ub); 601 602 if(MP_SIGN(a) == MP_SIGN(b)) { 603 /* Same sign -- add magnitudes, preserve sign of addends */ 604 mp_digit carry; 605 606 if(!s_pad(K, c, max)) 607 return MP_MEMORY; 608 609 carry = s_uadd(MP_DIGITS(a), MP_DIGITS(b), MP_DIGITS(c), ua, ub); 610 uc = max; 611 612 if(carry) { 613 if(!s_pad(K, c, max + 1)) 614 return MP_MEMORY; 615 616 c->digits[max] = carry; 617 ++uc; 618 } 619 620 MP_USED(c) = uc; 621 MP_SIGN(c) = MP_SIGN(a); 622 623 } 624 else { 625 /* Different signs -- subtract magnitudes, preserve sign of greater */ 626 mp_int x, y; 627 int cmp = s_ucmp(a, b); /* magnitude comparision, sign ignored */ 628 629 /* Set x to max(a, b), y to min(a, b) to simplify later code. 630 A special case yields zero for equal magnitudes. 631 */ 632 if(cmp == 0) { 633 mp_int_zero(c); 634 return MP_OK; 635 } 636 else if(cmp < 0) { 637 x = b; y = a; 638 } 639 else { 640 x = a; y = b; 641 } 642 643 if(!s_pad(K, c, MP_USED(x))) 644 return MP_MEMORY; 645 646 /* Subtract smaller from larger */ 647 s_usub(MP_DIGITS(x), MP_DIGITS(y), MP_DIGITS(c), MP_USED(x), MP_USED(y)); 648 MP_USED(c) = MP_USED(x); 649 CLAMP(c); 650 651 /* Give result the sign of the larger */ 652 MP_SIGN(c) = MP_SIGN(x); 653 } 654 655 return MP_OK; 656 } 657 658 /* }}} */ 659 660 /* {{{ mp_int_add_value(a, value, c) */ 661 662 mp_result mp_int_add_value(klisp_State *K, mp_int a, mp_small value, mp_int c) 663 { 664 mpz_t vtmp; 665 mp_digit vbuf[MP_VALUE_DIGITS(value)]; 666 667 s_fake(&vtmp, value, vbuf); 668 669 return mp_int_add(K, a, &vtmp, c); 670 } 671 672 /* }}} */ 673 674 /* {{{ mp_int_sub(a, b, c) */ 675 676 mp_result mp_int_sub(klisp_State *K, mp_int a, mp_int b, mp_int c) 677 { 678 mp_size ua, ub, uc, max; 679 680 CHECK(a != NULL && b != NULL && c != NULL); 681 682 ua = MP_USED(a); ub = MP_USED(b); uc = MP_USED(c); 683 max = MAX(ua, ub); 684 685 if(MP_SIGN(a) != MP_SIGN(b)) { 686 /* Different signs -- add magnitudes and keep sign of a */ 687 mp_digit carry; 688 689 if(!s_pad(K, c, max)) 690 return MP_MEMORY; 691 692 carry = s_uadd(MP_DIGITS(a), MP_DIGITS(b), MP_DIGITS(c), ua, ub); 693 uc = max; 694 695 if(carry) { 696 if(!s_pad(K, c, max + 1)) 697 return MP_MEMORY; 698 699 c->digits[max] = carry; 700 ++uc; 701 } 702 703 MP_USED(c) = uc; 704 MP_SIGN(c) = MP_SIGN(a); 705 706 } 707 else { 708 /* Same signs -- subtract magnitudes */ 709 mp_int x, y; 710 mp_sign osign; 711 int cmp = s_ucmp(a, b); 712 713 if(!s_pad(K, c, max)) 714 return MP_MEMORY; 715 716 if(cmp >= 0) { 717 x = a; y = b; osign = MP_ZPOS; 718 } 719 else { 720 x = b; y = a; osign = MP_NEG; 721 } 722 723 if(MP_SIGN(a) == MP_NEG && cmp != 0) 724 osign = 1 - osign; 725 726 s_usub(MP_DIGITS(x), MP_DIGITS(y), MP_DIGITS(c), MP_USED(x), MP_USED(y)); 727 MP_USED(c) = MP_USED(x); 728 CLAMP(c); 729 730 MP_SIGN(c) = osign; 731 } 732 733 return MP_OK; 734 } 735 736 /* }}} */ 737 738 /* {{{ mp_int_sub_value(a, value, c) */ 739 740 mp_result mp_int_sub_value(klisp_State *K, mp_int a, mp_small value, mp_int c) 741 { 742 mpz_t vtmp; 743 mp_digit vbuf[MP_VALUE_DIGITS(value)]; 744 745 s_fake(&vtmp, value, vbuf); 746 747 return mp_int_sub(K, a, &vtmp, c); 748 } 749 750 /* }}} */ 751 752 /* {{{ mp_int_mul(a, b, c) */ 753 754 mp_result mp_int_mul(klisp_State *K, mp_int a, mp_int b, mp_int c) 755 { 756 mp_digit *out; 757 mp_size osize, ua, ub, p = 0; 758 mp_sign osign; 759 760 CHECK(a != NULL && b != NULL && c != NULL); 761 762 /* If either input is zero, we can shortcut multiplication */ 763 if(mp_int_compare_zero(a) == 0 || mp_int_compare_zero(b) == 0) { 764 mp_int_zero(c); 765 return MP_OK; 766 } 767 768 /* Output is positive if inputs have same sign, otherwise negative */ 769 osign = (MP_SIGN(a) == MP_SIGN(b)) ? MP_ZPOS : MP_NEG; 770 771 /* If the output is not identical to any of the inputs, we'll write 772 the results directly; otherwise, allocate a temporary space. */ 773 ua = MP_USED(a); ub = MP_USED(b); 774 osize = MAX(ua, ub); 775 osize = 4 * ((osize + 1) / 2); 776 777 if(c == a || c == b) { 778 p = ROUND_PREC(osize); 779 p = MAX(p, default_precision); 780 781 if((out = s_alloc(K, p)) == NULL) 782 return MP_MEMORY; 783 } 784 else { 785 if(!s_pad(K, c, osize)) 786 return MP_MEMORY; 787 788 out = MP_DIGITS(c); 789 } 790 ZERO(out, osize); 791 792 if(!s_kmul(K, MP_DIGITS(a), MP_DIGITS(b), out, ua, ub)) 793 return MP_MEMORY; 794 795 /* If we allocated a new buffer, get rid of whatever memory c was 796 already using, and fix up its fields to reflect that. 797 */ 798 if(out != MP_DIGITS(c)) { 799 if((void *) MP_DIGITS(c) != (void *) &MP_SINGLE(c)) 800 s_free(K, MP_DIGITS(c), MP_ALLOC(c)); 801 MP_DIGITS(c) = out; 802 MP_ALLOC(c) = p; 803 } 804 805 MP_USED(c) = osize; /* might not be true, but we'll fix it ... */ 806 CLAMP(c); /* ... right here */ 807 MP_SIGN(c) = osign; 808 809 return MP_OK; 810 } 811 812 /* }}} */ 813 814 /* {{{ mp_int_mul_value(a, value, c) */ 815 816 mp_result mp_int_mul_value(klisp_State *K, mp_int a, mp_small value, 817 mp_int c) 818 { 819 mpz_t vtmp; 820 mp_digit vbuf[MP_VALUE_DIGITS(value)]; 821 822 s_fake(&vtmp, value, vbuf); 823 824 return mp_int_mul(K, a, &vtmp, c); 825 } 826 827 /* }}} */ 828 829 /* {{{ mp_int_mul_pow2(a, p2, c) */ 830 831 mp_result mp_int_mul_pow2(klisp_State *K, mp_int a, mp_small p2, mp_int c) 832 { 833 mp_result res; 834 CHECK(a != NULL && c != NULL && p2 >= 0); 835 836 if((res = mp_int_copy(K, a, c)) != MP_OK) 837 return res; 838 839 if(s_qmul(K, c, (mp_size) p2)) 840 return MP_OK; 841 else 842 return MP_MEMORY; 843 } 844 845 /* }}} */ 846 847 /* {{{ mp_int_sqr(a, c) */ 848 849 mp_result mp_int_sqr(klisp_State *K, mp_int a, mp_int c) 850 { 851 mp_digit *out; 852 mp_size osize, p = 0; 853 854 CHECK(a != NULL && c != NULL); 855 856 /* Get a temporary buffer big enough to hold the result */ 857 osize = (mp_size) 4 * ((MP_USED(a) + 1) / 2); 858 if(a == c) { 859 p = ROUND_PREC(osize); 860 p = MAX(p, default_precision); 861 862 if((out = s_alloc(K, p)) == NULL) 863 return MP_MEMORY; 864 } 865 else { 866 if(!s_pad(K, c, osize)) 867 return MP_MEMORY; 868 869 out = MP_DIGITS(c); 870 } 871 ZERO(out, osize); 872 873 s_ksqr(K, MP_DIGITS(a), out, MP_USED(a)); 874 875 /* Get rid of whatever memory c was already using, and fix up its 876 fields to reflect the new digit array it's using 877 */ 878 if(out != MP_DIGITS(c)) { 879 if((void *) MP_DIGITS(c) != (void *) &MP_SINGLE(c)) 880 s_free(K, MP_DIGITS(c), MP_ALLOC(c)); 881 MP_DIGITS(c) = out; 882 MP_ALLOC(c) = p; 883 } 884 885 MP_USED(c) = osize; /* might not be true, but we'll fix it ... */ 886 CLAMP(c); /* ... right here */ 887 MP_SIGN(c) = MP_ZPOS; 888 889 return MP_OK; 890 } 891 892 /* }}} */ 893 894 /* {{{ mp_int_div(a, b, q, r) */ 895 896 mp_result mp_int_div(klisp_State *K, mp_int a, mp_int b, mp_int q, mp_int r) 897 { 898 int cmp, last = 0, lg; 899 mp_result res = MP_OK; 900 mpz_t temp[2]; 901 mp_int qout, rout; 902 mp_sign sa = MP_SIGN(a), sb = MP_SIGN(b); 903 904 CHECK(a != NULL && b != NULL && q != r); 905 906 if(CMPZ(b) == 0) 907 return MP_UNDEF; 908 else if((cmp = s_ucmp(a, b)) < 0) { 909 /* If |a| < |b|, no division is required: 910 q = 0, r = a 911 */ 912 if(r && (res = mp_int_copy(K, a, r)) != MP_OK) 913 return res; 914 915 if(q) 916 mp_int_zero(q); 917 918 return MP_OK; 919 } 920 else if(cmp == 0) { 921 /* If |a| = |b|, no division is required: 922 q = 1 or -1, r = 0 923 */ 924 if(r) 925 mp_int_zero(r); 926 927 if(q) { 928 mp_int_zero(q); 929 q->digits[0] = 1; 930 931 if(sa != sb) 932 MP_SIGN(q) = MP_NEG; 933 } 934 935 return MP_OK; 936 } 937 938 /* When |a| > |b|, real division is required. We need someplace to 939 store quotient and remainder, but q and r are allowed to be NULL 940 or to overlap with the inputs. 941 */ 942 if((lg = s_isp2(b)) < 0) { 943 if(q && b != q) { 944 if((res = mp_int_copy(K, a, q)) != MP_OK) 945 goto CLEANUP; 946 else 947 qout = q; 948 } 949 else { 950 qout = TEMP(last); 951 SETUP(mp_int_init_copy(K, TEMP(last), a), last); 952 } 953 954 if(r && a != r) { 955 if((res = mp_int_copy(K, b, r)) != MP_OK) 956 goto CLEANUP; 957 else 958 rout = r; 959 } 960 else { 961 rout = TEMP(last); 962 SETUP(mp_int_init_copy(K, TEMP(last), b), last); 963 } 964 965 if((res = s_udiv(K, qout, rout)) != MP_OK) goto CLEANUP; 966 } 967 else { 968 if(q && (res = mp_int_copy(K, a, q)) != MP_OK) goto CLEANUP; 969 if(r && (res = mp_int_copy(K, a, r)) != MP_OK) goto CLEANUP; 970 971 if(q) s_qdiv(q, (mp_size) lg); qout = q; 972 if(r) s_qmod(r, (mp_size) lg); rout = r; 973 } 974 975 /* Recompute signs for output */ 976 if(rout) { 977 MP_SIGN(rout) = sa; 978 if(CMPZ(rout) == 0) 979 MP_SIGN(rout) = MP_ZPOS; 980 } 981 if(qout) { 982 MP_SIGN(qout) = (sa == sb) ? MP_ZPOS : MP_NEG; 983 if(CMPZ(qout) == 0) 984 MP_SIGN(qout) = MP_ZPOS; 985 } 986 987 if(q && (res = mp_int_copy(K, qout, q)) != MP_OK) goto CLEANUP; 988 if(r && (res = mp_int_copy(K, rout, r)) != MP_OK) goto CLEANUP; 989 990 CLEANUP: 991 while(--last >= 0) 992 mp_int_clear(K, TEMP(last)); 993 994 return res; 995 } 996 997 /* }}} */ 998 999 /* {{{ mp_int_mod(a, m, c) */ 1000 1001 mp_result mp_int_mod(klisp_State *K, mp_int a, mp_int m, mp_int c) 1002 { 1003 mp_result res; 1004 mpz_t tmp; 1005 mp_int out; 1006 1007 if(m == c) { 1008 mp_int_init(&tmp); 1009 out = &tmp; 1010 } 1011 else { 1012 out = c; 1013 } 1014 1015 if((res = mp_int_div(K, a, m, NULL, out)) != MP_OK) 1016 goto CLEANUP; 1017 1018 if(CMPZ(out) < 0) 1019 res = mp_int_add(K, out, m, c); 1020 else 1021 res = mp_int_copy(K, out, c); 1022 1023 CLEANUP: 1024 if(out != c) 1025 mp_int_clear(K, &tmp); 1026 1027 return res; 1028 } 1029 1030 /* }}} */ 1031 1032 /* {{{ mp_int_div_value(a, value, q, r) */ 1033 1034 mp_result mp_int_div_value(klisp_State *K, mp_int a, mp_small value, 1035 mp_int q, mp_small *r) 1036 { 1037 mpz_t vtmp, rtmp; 1038 mp_digit vbuf[MP_VALUE_DIGITS(value)]; 1039 mp_result res; 1040 1041 mp_int_init(&rtmp); 1042 s_fake(&vtmp, value, vbuf); 1043 1044 if((res = mp_int_div(K, a, &vtmp, q, &rtmp)) != MP_OK) 1045 goto CLEANUP; 1046 1047 if(r) 1048 (void) mp_int_to_int(&rtmp, r); /* can't fail */ 1049 1050 CLEANUP: 1051 mp_int_clear(K, &rtmp); 1052 return res; 1053 } 1054 1055 /* }}} */ 1056 1057 /* {{{ mp_int_div_pow2(a, p2, q, r) */ 1058 1059 mp_result mp_int_div_pow2(klisp_State *K, mp_int a, mp_small p2, mp_int q, 1060 mp_int r) 1061 { 1062 mp_result res = MP_OK; 1063 1064 CHECK(a != NULL && p2 >= 0 && q != r); 1065 1066 if(q != NULL && (res = mp_int_copy(K, a, q)) == MP_OK) 1067 s_qdiv(q, (mp_size) p2); 1068 1069 if(res == MP_OK && r != NULL && (res = mp_int_copy(K, a, r)) == MP_OK) 1070 s_qmod(r, (mp_size) p2); 1071 1072 return res; 1073 } 1074 1075 /* }}} */ 1076 1077 /* {{{ mp_int_expt(a, b, c) */ 1078 1079 mp_result mp_int_expt(klisp_State *K, mp_int a, mp_small b, mp_int c) 1080 { 1081 mpz_t t; 1082 mp_result res; 1083 unsigned int v = abs(b); 1084 1085 CHECK(c != NULL); 1086 1087 CHECK(b >= 0 && c != NULL); 1088 1089 if((res = mp_int_init_copy(K, &t, a)) != MP_OK) 1090 return res; 1091 1092 (void) mp_int_set_value(K, c, 1); 1093 while(v != 0) { 1094 if(v & 1) { 1095 if((res = mp_int_mul(K, c, &t, c)) != MP_OK) 1096 goto CLEANUP; 1097 } 1098 1099 v >>= 1; 1100 if(v == 0) break; 1101 1102 if((res = mp_int_sqr(K, &t, &t)) != MP_OK) 1103 goto CLEANUP; 1104 } 1105 1106 CLEANUP: 1107 mp_int_clear(K, &t); 1108 return res; 1109 } 1110 1111 /* }}} */ 1112 1113 /* {{{ mp_int_expt_value(a, b, c) */ 1114 1115 mp_result mp_int_expt_value(klisp_State *K, mp_small a, mp_small b, mp_int c) 1116 { 1117 mpz_t t; 1118 mp_result res; 1119 unsigned int v = abs(b); 1120 1121 CHECK(b >= 0 && c != NULL); 1122 1123 if((res = mp_int_init_value(K, &t, a)) != MP_OK) 1124 return res; 1125 1126 (void) mp_int_set_value(K, c, 1); 1127 while(v != 0) { 1128 if(v & 1) { 1129 if((res = mp_int_mul(K, c, &t, c)) != MP_OK) 1130 goto CLEANUP; 1131 } 1132 1133 v >>= 1; 1134 if(v == 0) break; 1135 1136 if((res = mp_int_sqr(K, &t, &t)) != MP_OK) 1137 goto CLEANUP; 1138 } 1139 1140 CLEANUP: 1141 mp_int_clear(K, &t); 1142 return res; 1143 } 1144 1145 /* }}} */ 1146 1147 /* {{{ mp_int_expt_full(a, b, c) */ 1148 1149 mp_result mp_int_expt_full(klisp_State *K, mp_int a, mp_int b, mp_int c) 1150 { 1151 mpz_t t; 1152 mp_result res; 1153 int ix, jx; 1154 1155 CHECK(a != NULL && b != NULL && c != NULL); 1156 1157 if ((res = mp_int_init_copy(K, &t, a)) != MP_OK) 1158 return res; 1159 1160 (void) mp_int_set_value(K, c, 1); 1161 for (ix = 0; ix < MP_USED(b); ++ix) { 1162 mp_digit d = b->digits[ix]; 1163 1164 for (jx = 0; jx < MP_DIGIT_BIT; ++jx) { 1165 if (d & 1) { 1166 if ((res = mp_int_mul(K, c, &t, c)) != MP_OK) 1167 goto CLEANUP; 1168 } 1169 1170 d >>= 1; 1171 if (d == 0 && ix + 1 == MP_USED(b)) 1172 break; 1173 if ((res = mp_int_sqr(K, &t, &t)) != MP_OK) 1174 goto CLEANUP; 1175 } 1176 } 1177 1178 CLEANUP: 1179 mp_int_clear(K, &t); 1180 return res; 1181 } 1182 1183 /* }}} */ 1184 1185 /* {{{ mp_int_compare(a, b) */ 1186 1187 int mp_int_compare(mp_int a, mp_int b) 1188 { 1189 mp_sign sa; 1190 1191 CHECK(a != NULL && b != NULL); 1192 1193 sa = MP_SIGN(a); 1194 if(sa == MP_SIGN(b)) { 1195 int cmp = s_ucmp(a, b); 1196 1197 /* If they're both zero or positive, the normal comparison 1198 applies; if both negative, the sense is reversed. */ 1199 if(sa == MP_ZPOS) 1200 return cmp; 1201 else 1202 return -cmp; 1203 1204 } 1205 else { 1206 if(sa == MP_ZPOS) 1207 return 1; 1208 else 1209 return -1; 1210 } 1211 } 1212 1213 /* }}} */ 1214 1215 /* {{{ mp_int_compare_unsigned(a, b) */ 1216 1217 int mp_int_compare_unsigned(mp_int a, mp_int b) 1218 { 1219 NRCHECK(a != NULL && b != NULL); 1220 1221 return s_ucmp(a, b); 1222 } 1223 1224 /* }}} */ 1225 1226 /* {{{ mp_int_compare_zero(z) */ 1227 1228 int mp_int_compare_zero(mp_int z) 1229 { 1230 NRCHECK(z != NULL); 1231 1232 if(MP_USED(z) == 1 && z->digits[0] == 0) 1233 return 0; 1234 else if(MP_SIGN(z) == MP_ZPOS) 1235 return 1; 1236 else 1237 return -1; 1238 } 1239 1240 /* }}} */ 1241 1242 /* {{{ mp_int_compare_value(z, value) */ 1243 1244 int mp_int_compare_value(mp_int z, mp_small value) 1245 { 1246 mp_sign vsign = (value < 0) ? MP_NEG : MP_ZPOS; 1247 int cmp; 1248 1249 CHECK(z != NULL); 1250 1251 if(vsign == MP_SIGN(z)) { 1252 cmp = s_vcmp(z, value); 1253 1254 if(vsign == MP_ZPOS) 1255 return cmp; 1256 else 1257 return -cmp; 1258 } 1259 else { 1260 if(value < 0) 1261 return 1; 1262 else 1263 return -1; 1264 } 1265 } 1266 1267 /* }}} */ 1268 1269 /* {{{ mp_int_exptmod(a, b, m, c) */ 1270 1271 mp_result mp_int_exptmod(klisp_State *K, mp_int a, mp_int b, mp_int m, 1272 mp_int c) 1273 { 1274 mp_result res; 1275 mp_size um; 1276 mpz_t temp[3]; 1277 mp_int s; 1278 int last = 0; 1279 1280 CHECK(a != NULL && b != NULL && c != NULL && m != NULL); 1281 1282 /* Zero moduli and negative exponents are not considered. */ 1283 if(CMPZ(m) == 0) 1284 return MP_UNDEF; 1285 if(CMPZ(b) < 0) 1286 return MP_RANGE; 1287 1288 um = MP_USED(m); 1289 SETUP(mp_int_init_size(K, TEMP(0), 2 * um), last); 1290 SETUP(mp_int_init_size(K, TEMP(1), 2 * um), last); 1291 1292 if(c == b || c == m) { 1293 SETUP(mp_int_init_size(K, TEMP(2), 2 * um), last); 1294 s = TEMP(2); 1295 } 1296 else { 1297 s = c; 1298 } 1299 1300 if((res = mp_int_mod(K, a, m, TEMP(0))) != MP_OK) goto CLEANUP; 1301 1302 if((res = s_brmu(K, TEMP(1), m)) != MP_OK) goto CLEANUP; 1303 1304 if((res = s_embar(K, TEMP(0), b, m, TEMP(1), s)) != MP_OK) 1305 goto CLEANUP; 1306 1307 res = mp_int_copy(K, s, c); 1308 1309 CLEANUP: 1310 while(--last >= 0) 1311 mp_int_clear(K, TEMP(last)); 1312 1313 return res; 1314 } 1315 1316 /* }}} */ 1317 1318 /* {{{ mp_int_exptmod_evalue(a, value, m, c) */ 1319 1320 mp_result mp_int_exptmod_evalue(klisp_State *K, mp_int a, mp_small value, 1321 mp_int m, mp_int c) 1322 { 1323 mpz_t vtmp; 1324 mp_digit vbuf[MP_VALUE_DIGITS(value)]; 1325 1326 s_fake(&vtmp, value, vbuf); 1327 1328 return mp_int_exptmod(K, a, &vtmp, m, c); 1329 } 1330 1331 /* }}} */ 1332 1333 /* {{{ mp_int_exptmod_bvalue(v, b, m, c) */ 1334 1335 mp_result mp_int_exptmod_bvalue(klisp_State *K, mp_small value, mp_int b, 1336 mp_int m, mp_int c) 1337 { 1338 mpz_t vtmp; 1339 mp_digit vbuf[MP_VALUE_DIGITS(value)]; 1340 1341 s_fake(&vtmp, value, vbuf); 1342 1343 return mp_int_exptmod(K, &vtmp, b, m, c); 1344 } 1345 1346 /* }}} */ 1347 1348 /* {{{ mp_int_exptmod_known(a, b, m, mu, c) */ 1349 1350 mp_result mp_int_exptmod_known(klisp_State *K, mp_int a, mp_int b, mp_int m, 1351 mp_int mu, mp_int c) 1352 { 1353 mp_result res; 1354 mp_size um; 1355 mpz_t temp[2]; 1356 mp_int s; 1357 int last = 0; 1358 1359 CHECK(a && b && m && c); 1360 1361 /* Zero moduli and negative exponents are not considered. */ 1362 if(CMPZ(m) == 0) 1363 return MP_UNDEF; 1364 if(CMPZ(b) < 0) 1365 return MP_RANGE; 1366 1367 um = MP_USED(m); 1368 SETUP(mp_int_init_size(K, TEMP(0), 2 * um), last); 1369 1370 if(c == b || c == m) { 1371 SETUP(mp_int_init_size(K, TEMP(1), 2 * um), last); 1372 s = TEMP(1); 1373 } 1374 else { 1375 s = c; 1376 } 1377 1378 if((res = mp_int_mod(K, a, m, TEMP(0))) != MP_OK) goto CLEANUP; 1379 1380 if((res = s_embar(K, TEMP(0), b, m, mu, s)) != MP_OK) 1381 goto CLEANUP; 1382 1383 res = mp_int_copy(K, s, c); 1384 1385 CLEANUP: 1386 while(--last >= 0) 1387 mp_int_clear(K, TEMP(last)); 1388 1389 return res; 1390 } 1391 1392 /* }}} */ 1393 1394 /* {{{ mp_int_redux_const(m, c) */ 1395 1396 mp_result mp_int_redux_const(klisp_State *K, mp_int m, mp_int c) 1397 { 1398 CHECK(m != NULL && c != NULL && m != c); 1399 1400 return s_brmu(K, c, m); 1401 } 1402 1403 /* }}} */ 1404 1405 /* {{{ mp_int_invmod(a, m, c) */ 1406 1407 mp_result mp_int_invmod(klisp_State *K, mp_int a, mp_int m, mp_int c) 1408 { 1409 mp_result res; 1410 mp_sign sa; 1411 int last = 0; 1412 mpz_t temp[2]; 1413 1414 CHECK(a != NULL && m != NULL && c != NULL); 1415 1416 if(CMPZ(a) == 0 || CMPZ(m) <= 0) 1417 return MP_RANGE; 1418 1419 sa = MP_SIGN(a); /* need this for the result later */ 1420 1421 for(last = 0; last < 2; ++last) 1422 mp_int_init(TEMP(last)); 1423 1424 if((res = mp_int_egcd(K, a, m, TEMP(0), TEMP(1), NULL)) != MP_OK) 1425 goto CLEANUP; 1426 1427 if(mp_int_compare_value(TEMP(0), 1) != 0) { 1428 res = MP_UNDEF; 1429 goto CLEANUP; 1430 } 1431 1432 /* It is first necessary to constrain the value to the proper range */ 1433 if((res = mp_int_mod(K, TEMP(1), m, TEMP(1))) != MP_OK) 1434 goto CLEANUP; 1435 1436 /* Now, if 'a' was originally negative, the value we have is 1437 actually the magnitude of the negative representative; to get the 1438 positive value we have to subtract from the modulus. Otherwise, 1439 the value is okay as it stands. 1440 */ 1441 if(sa == MP_NEG) 1442 res = mp_int_sub(K, m, TEMP(1), c); 1443 else 1444 res = mp_int_copy(K, TEMP(1), c); 1445 1446 CLEANUP: 1447 while(--last >= 0) 1448 mp_int_clear(K, TEMP(last)); 1449 1450 return res; 1451 } 1452 1453 /* }}} */ 1454 1455 /* {{{ mp_int_gcd(a, b, c) */ 1456 1457 /* Binary GCD algorithm due to Josef Stein, 1961 */ 1458 mp_result mp_int_gcd(klisp_State *K, mp_int a, mp_int b, mp_int c) 1459 { 1460 int ca, cb, k = 0; 1461 mpz_t u, v, t; 1462 mp_result res; 1463 1464 CHECK(a != NULL && b != NULL && c != NULL); 1465 1466 ca = CMPZ(a); 1467 cb = CMPZ(b); 1468 if(ca == 0 && cb == 0) 1469 return MP_UNDEF; 1470 else if(ca == 0) 1471 return mp_int_abs(K, b, c); 1472 else if(cb == 0) 1473 return mp_int_abs(K, a, c); 1474 1475 mp_int_init(&t); 1476 if((res = mp_int_init_copy(K, &u, a)) != MP_OK) 1477 goto U; 1478 if((res = mp_int_init_copy(K, &v, b)) != MP_OK) 1479 goto V; 1480 1481 MP_SIGN(&u) = MP_ZPOS; MP_SIGN(&v) = MP_ZPOS; 1482 1483 { /* Divide out common factors of 2 from u and v */ 1484 int div2_u = s_dp2k(&u), div2_v = s_dp2k(&v); 1485 1486 k = MIN(div2_u, div2_v); 1487 s_qdiv(&u, (mp_size) k); 1488 s_qdiv(&v, (mp_size) k); 1489 } 1490 1491 if(mp_int_is_odd(&u)) { 1492 if((res = mp_int_neg(K, &v, &t)) != MP_OK) 1493 goto CLEANUP; 1494 } 1495 else { 1496 if((res = mp_int_copy(K, &u, &t)) != MP_OK) 1497 goto CLEANUP; 1498 } 1499 1500 for(;;) { 1501 s_qdiv(&t, s_dp2k(&t)); 1502 1503 if(CMPZ(&t) > 0) { 1504 if((res = mp_int_copy(K, &t, &u)) != MP_OK) 1505 goto CLEANUP; 1506 } 1507 else { 1508 if((res = mp_int_neg(K, &t, &v)) != MP_OK) 1509 goto CLEANUP; 1510 } 1511 1512 if((res = mp_int_sub(K, &u, &v, &t)) != MP_OK) 1513 goto CLEANUP; 1514 1515 if(CMPZ(&t) == 0) 1516 break; 1517 } 1518 1519 if((res = mp_int_abs(K, &u, c)) != MP_OK) 1520 goto CLEANUP; 1521 if(!s_qmul(K, c, (mp_size) k)) 1522 res = MP_MEMORY; 1523 1524 CLEANUP: 1525 mp_int_clear(K, &v); 1526 V: mp_int_clear(K, &u); 1527 U: mp_int_clear(K, &t); 1528 1529 return res; 1530 } 1531 1532 /* }}} */ 1533 1534 /* {{{ mp_int_egcd(a, b, c, x, y) */ 1535 1536 /* This is the binary GCD algorithm again, but this time we keep track 1537 of the elementary matrix operations as we go, so we can get values 1538 x and y satisfying c = ax + by. 1539 */ 1540 mp_result mp_int_egcd(klisp_State *K, mp_int a, mp_int b, mp_int c, 1541 mp_int x, mp_int y) 1542 { 1543 int k, last = 0, ca, cb; 1544 mpz_t temp[8]; 1545 mp_result res; 1546 1547 CHECK(a != NULL && b != NULL && c != NULL && 1548 (x != NULL || y != NULL)); 1549 1550 ca = CMPZ(a); 1551 cb = CMPZ(b); 1552 if(ca == 0 && cb == 0) 1553 return MP_UNDEF; 1554 else if(ca == 0) { 1555 if((res = mp_int_abs(K, b, c)) != MP_OK) return res; 1556 mp_int_zero(x); (void) mp_int_set_value(K, y, 1); return MP_OK; 1557 } 1558 else if(cb == 0) { 1559 if((res = mp_int_abs(K, a, c)) != MP_OK) return res; 1560 (void) mp_int_set_value(K, x, 1); mp_int_zero(y); return MP_OK; 1561 } 1562 1563 /* Initialize temporaries: 1564 A:0, B:1, C:2, D:3, u:4, v:5, ou:6, ov:7 */ 1565 for(last = 0; last < 4; ++last) 1566 mp_int_init(TEMP(last)); 1567 TEMP(0)->digits[0] = 1; 1568 TEMP(3)->digits[0] = 1; 1569 1570 SETUP(mp_int_init_copy(K, TEMP(4), a), last); 1571 SETUP(mp_int_init_copy(K, TEMP(5), b), last); 1572 1573 /* We will work with absolute values here */ 1574 MP_SIGN(TEMP(4)) = MP_ZPOS; 1575 MP_SIGN(TEMP(5)) = MP_ZPOS; 1576 1577 { /* Divide out common factors of 2 from u and v */ 1578 int div2_u = s_dp2k(TEMP(4)), div2_v = s_dp2k(TEMP(5)); 1579 1580 k = MIN(div2_u, div2_v); 1581 s_qdiv(TEMP(4), k); 1582 s_qdiv(TEMP(5), k); 1583 } 1584 1585 SETUP(mp_int_init_copy(K, TEMP(6), TEMP(4)), last); 1586 SETUP(mp_int_init_copy(K, TEMP(7), TEMP(5)), last); 1587 1588 for(;;) { 1589 while(mp_int_is_even(TEMP(4))) { 1590 s_qdiv(TEMP(4), 1); 1591 1592 if(mp_int_is_odd(TEMP(0)) || mp_int_is_odd(TEMP(1))) { 1593 if((res = mp_int_add(K, TEMP(0), TEMP(7), TEMP(0))) != MP_OK) 1594 goto CLEANUP; 1595 if((res = mp_int_sub(K, TEMP(1), TEMP(6), TEMP(1))) != MP_OK) 1596 goto CLEANUP; 1597 } 1598 1599 s_qdiv(TEMP(0), 1); 1600 s_qdiv(TEMP(1), 1); 1601 } 1602 1603 while(mp_int_is_even(TEMP(5))) { 1604 s_qdiv(TEMP(5), 1); 1605 1606 if(mp_int_is_odd(TEMP(2)) || mp_int_is_odd(TEMP(3))) { 1607 if((res = mp_int_add(K, TEMP(2), TEMP(7), TEMP(2))) != MP_OK) 1608 goto CLEANUP; 1609 if((res = mp_int_sub(K, TEMP(3), TEMP(6), TEMP(3))) != MP_OK) 1610 goto CLEANUP; 1611 } 1612 1613 s_qdiv(TEMP(2), 1); 1614 s_qdiv(TEMP(3), 1); 1615 } 1616 1617 if(mp_int_compare(TEMP(4), TEMP(5)) >= 0) { 1618 if((res = mp_int_sub(K, TEMP(4), TEMP(5), TEMP(4))) != MP_OK) goto CLEANUP; 1619 if((res = mp_int_sub(K, TEMP(0), TEMP(2), TEMP(0))) != MP_OK) goto CLEANUP; 1620 if((res = mp_int_sub(K, TEMP(1), TEMP(3), TEMP(1))) != MP_OK) goto CLEANUP; 1621 } 1622 else { 1623 if((res = mp_int_sub(K, TEMP(5), TEMP(4), TEMP(5))) != MP_OK) goto CLEANUP; 1624 if((res = mp_int_sub(K, TEMP(2), TEMP(0), TEMP(2))) != MP_OK) goto CLEANUP; 1625 if((res = mp_int_sub(K, TEMP(3), TEMP(1), TEMP(3))) != MP_OK) goto CLEANUP; 1626 } 1627 1628 if(CMPZ(TEMP(4)) == 0) { 1629 if(x && (res = mp_int_copy(K, TEMP(2), x)) != MP_OK) goto CLEANUP; 1630 if(y && (res = mp_int_copy(K, TEMP(3), y)) != MP_OK) goto CLEANUP; 1631 if(c) { 1632 if(!s_qmul(K, TEMP(5), k)) { 1633 res = MP_MEMORY; 1634 goto CLEANUP; 1635 } 1636 1637 res = mp_int_copy(K, TEMP(5), c); 1638 } 1639 1640 break; 1641 } 1642 } 1643 1644 CLEANUP: 1645 while(--last >= 0) 1646 mp_int_clear(K, TEMP(last)); 1647 1648 return res; 1649 } 1650 1651 /* }}} */ 1652 1653 /* {{{ mp_int_lcm(a, b, c) */ 1654 1655 mp_result mp_int_lcm(klisp_State *K, mp_int a, mp_int b, mp_int c) 1656 { 1657 mpz_t lcm; 1658 mp_result res; 1659 1660 CHECK(a != NULL && b != NULL && c != NULL); 1661 1662 /* Since a * b = gcd(a, b) * lcm(a, b), we can compute 1663 lcm(a, b) = (a / gcd(a, b)) * b. 1664 1665 This formulation insures everything works even if the input 1666 variables share space. 1667 */ 1668 if((res = mp_int_init(&lcm)) != MP_OK) 1669 return res; 1670 if((res = mp_int_gcd(K, a, b, &lcm)) != MP_OK) 1671 goto CLEANUP; 1672 if((res = mp_int_div(K, a, &lcm, &lcm, NULL)) != MP_OK) 1673 goto CLEANUP; 1674 if((res = mp_int_mul(K, &lcm, b, &lcm)) != MP_OK) 1675 goto CLEANUP; 1676 1677 res = mp_int_copy(K, &lcm, c); 1678 1679 CLEANUP: 1680 mp_int_clear(K, &lcm); 1681 1682 return res; 1683 } 1684 1685 /* }}} */ 1686 1687 /* {{{ mp_int_divisible_value(a, v) */ 1688 1689 int mp_int_divisible_value(klisp_State *K, mp_int a, mp_small v) 1690 { 1691 mp_small rem = 0; 1692 1693 if(mp_int_div_value(K, a, v, NULL, &rem) != MP_OK) 1694 return 0; 1695 1696 return rem == 0; 1697 } 1698 1699 /* }}} */ 1700 1701 /* {{{ mp_int_is_pow2(z) */ 1702 1703 int mp_int_is_pow2(mp_int z) 1704 { 1705 CHECK(z != NULL); 1706 1707 return s_isp2(z); 1708 } 1709 1710 /* }}} */ 1711 1712 /* {{{ mp_int_root(a, b, c) */ 1713 1714 /* Implementation of Newton's root finding method, based loosely on a 1715 patch contributed by Hal Finkel <half@halssoftware.com> 1716 modified by M. J. Fromberger. 1717 */ 1718 mp_result mp_int_root(klisp_State *K, mp_int a, mp_small b, mp_int c) 1719 { 1720 mp_result res = MP_OK; 1721 mpz_t temp[5]; 1722 int last = 0; 1723 int flips = 0; 1724 1725 CHECK(a != NULL && c != NULL && b > 0); 1726 1727 if(b == 1) { 1728 return mp_int_copy(K, a, c); 1729 } 1730 if(MP_SIGN(a) == MP_NEG) { 1731 if(b % 2 == 0) 1732 return MP_UNDEF; /* root does not exist for negative a with even b */ 1733 else 1734 flips = 1; 1735 } 1736 1737 SETUP(mp_int_init_copy(K, TEMP(last), a), last); 1738 SETUP(mp_int_init_copy(K, TEMP(last), a), last); 1739 SETUP(mp_int_init(TEMP(last)), last); 1740 SETUP(mp_int_init(TEMP(last)), last); 1741 SETUP(mp_int_init(TEMP(last)), last); 1742 1743 (void) mp_int_abs(K, TEMP(0), TEMP(0)); 1744 (void) mp_int_abs(K, TEMP(1), TEMP(1)); 1745 1746 for(;;) { 1747 if((res = mp_int_expt(K, TEMP(1), b, TEMP(2))) != MP_OK) 1748 goto CLEANUP; 1749 1750 if(mp_int_compare_unsigned(TEMP(2), TEMP(0)) <= 0) 1751 break; 1752 1753 if((res = mp_int_sub(K, TEMP(2), TEMP(0), TEMP(2))) != MP_OK) 1754 goto CLEANUP; 1755 if((res = mp_int_expt(K, TEMP(1), b - 1, TEMP(3))) != MP_OK) 1756 goto CLEANUP; 1757 if((res = mp_int_mul_value(K, TEMP(3), b, TEMP(3))) != MP_OK) 1758 goto CLEANUP; 1759 if((res = mp_int_div(K, TEMP(2), TEMP(3), TEMP(4), NULL)) != MP_OK) 1760 goto CLEANUP; 1761 if((res = mp_int_sub(K, TEMP(1), TEMP(4), TEMP(4))) != MP_OK) 1762 goto CLEANUP; 1763 1764 if(mp_int_compare_unsigned(TEMP(1), TEMP(4)) == 0) { 1765 if((res = mp_int_sub_value(K, TEMP(4), 1, TEMP(4))) != MP_OK) 1766 goto CLEANUP; 1767 } 1768 if((res = mp_int_copy(K, TEMP(4), TEMP(1))) != MP_OK) 1769 goto CLEANUP; 1770 } 1771 1772 if((res = mp_int_copy(K, TEMP(1), c)) != MP_OK) 1773 goto CLEANUP; 1774 1775 /* If the original value of a was negative, flip the output sign. */ 1776 if(flips) 1777 (void) mp_int_neg(K, c, c); /* cannot fail */ 1778 1779 CLEANUP: 1780 while(--last >= 0) 1781 mp_int_clear(K, TEMP(last)); 1782 1783 return res; 1784 } 1785 1786 /* }}} */ 1787 1788 /* {{{ mp_int_to_int(z, out) */ 1789 1790 mp_result mp_int_to_int(mp_int z, mp_small *out) 1791 { 1792 mp_usmall uv = 0; 1793 mp_size uz; 1794 mp_digit *dz; 1795 mp_sign sz; 1796 1797 CHECK(z != NULL); 1798 1799 /* Make sure the value is representable as an int */ 1800 sz = MP_SIGN(z); 1801 if((sz == MP_ZPOS && mp_int_compare_value(z, MP_SMALL_MAX) > 0) || 1802 mp_int_compare_value(z, MP_SMALL_MIN) < 0) 1803 return MP_RANGE; 1804 1805 uz = MP_USED(z); 1806 dz = MP_DIGITS(z) + uz - 1; 1807 1808 while(uz > 0) { 1809 uv <<= MP_DIGIT_BIT/2; 1810 uv = (uv << (MP_DIGIT_BIT/2)) | *dz--; 1811 --uz; 1812 } 1813 1814 if(out) 1815 *out = (sz == MP_NEG) ? -(mp_small)uv : (mp_small)uv; 1816 1817 return MP_OK; 1818 } 1819 1820 /* }}} */ 1821 1822 /* {{{ mp_int_to_uint(z, *out) */ 1823 1824 mp_result mp_int_to_uint(mp_int z, mp_usmall *out) 1825 { 1826 mp_usmall uv = 0; 1827 mp_size uz; 1828 mp_digit *dz; 1829 mp_sign sz; 1830 1831 CHECK(z != NULL); 1832 1833 /* Make sure the value is representable as an int */ 1834 sz = MP_SIGN(z); 1835 if(!(sz == MP_ZPOS && mp_int_compare_value(z, UINT_MAX) <= 0)) 1836 return MP_RANGE; 1837 1838 uz = MP_USED(z); 1839 dz = MP_DIGITS(z) + uz - 1; 1840 1841 while(uz > 0) { 1842 uv <<= MP_DIGIT_BIT/2; 1843 uv = (uv << (MP_DIGIT_BIT/2)) | *dz--; 1844 --uz; 1845 } 1846 1847 if(out) 1848 *out = uv; 1849 1850 return MP_OK; 1851 } 1852 1853 /* }}} */ 1854 1855 /* {{{ mp_int_to_string(z, radix, str, limit) */ 1856 1857 mp_result mp_int_to_string(klisp_State *K, mp_int z, mp_size radix, 1858 char *str, int limit) 1859 { 1860 mp_result res; 1861 int cmp = 0; 1862 1863 CHECK(z != NULL && str != NULL && limit >= 2); 1864 1865 if(radix < MP_MIN_RADIX || radix > MP_MAX_RADIX) 1866 return MP_RANGE; 1867 1868 if(CMPZ(z) == 0) { 1869 *str++ = s_val2ch(0, 0); /* changed to lowercase, Andres Navarro */ 1870 } 1871 else { 1872 mpz_t tmp; 1873 char *h, *t; 1874 1875 if((res = mp_int_init_copy(K, &tmp, z)) != MP_OK) 1876 return res; 1877 1878 if(MP_SIGN(z) == MP_NEG) { 1879 *str++ = '-'; 1880 --limit; 1881 } 1882 h = str; 1883 1884 /* Generate digits in reverse order until finished or limit reached */ 1885 for(/* */; limit > 0; --limit) { 1886 mp_digit d; 1887 1888 if((cmp = CMPZ(&tmp)) == 0) 1889 break; 1890 1891 d = s_ddiv(&tmp, (mp_digit)radix); 1892 *str++ = s_val2ch(d, 0); /* changed to lowercase, Andres Navarro */ 1893 } 1894 t = str - 1; 1895 1896 /* Put digits back in correct output order */ 1897 while(h < t) { 1898 char tc = *h; 1899 *h++ = *t; 1900 *t-- = tc; 1901 } 1902 1903 mp_int_clear(K, &tmp); 1904 } 1905 1906 *str = '\0'; 1907 if(cmp == 0) 1908 return MP_OK; 1909 else 1910 return MP_TRUNC; 1911 } 1912 1913 /* }}} */ 1914 1915 /* {{{ mp_int_string_len(z, radix) */ 1916 1917 mp_result mp_int_string_len(mp_int z, mp_size radix) 1918 { 1919 int len; 1920 1921 CHECK(z != NULL); 1922 1923 if(radix < MP_MIN_RADIX || radix > MP_MAX_RADIX) 1924 return MP_RANGE; 1925 1926 len = s_outlen(z, radix) + 1; /* for terminator */ 1927 1928 /* Allow for sign marker on negatives */ 1929 if(MP_SIGN(z) == MP_NEG) 1930 len += 1; 1931 1932 return len; 1933 } 1934 1935 /* }}} */ 1936 1937 /* {{{ mp_int_read_string(z, radix, *str) */ 1938 1939 /* Read zero-terminated string into z */ 1940 mp_result mp_int_read_string(klisp_State *K, mp_int z, mp_size radix, 1941 const char *str) 1942 { 1943 return mp_int_read_cstring(K, z, radix, str, NULL); 1944 } 1945 1946 /* }}} */ 1947 1948 /* {{{ mp_int_read_cstring(z, radix, *str, **end) */ 1949 1950 mp_result mp_int_read_cstring(klisp_State *K, mp_int z, mp_size radix, 1951 const char *str, char **end) 1952 { 1953 int ch; 1954 1955 CHECK(z != NULL && str != NULL); 1956 1957 if(radix < MP_MIN_RADIX || radix > MP_MAX_RADIX) 1958 return MP_RANGE; 1959 1960 /* Skip leading whitespace */ 1961 while(isspace((int)*str)) 1962 ++str; 1963 1964 /* Handle leading sign tag (+/-, positive default) */ 1965 switch(*str) { 1966 case '-': 1967 MP_SIGN(z) = MP_NEG; 1968 ++str; 1969 break; 1970 case '+': 1971 ++str; /* fallthrough */ 1972 default: 1973 MP_SIGN(z) = MP_ZPOS; 1974 break; 1975 } 1976 1977 /* Skip leading zeroes */ 1978 while((ch = s_ch2val(*str, radix)) == 0) 1979 ++str; 1980 1981 /* Make sure there is enough space for the value */ 1982 if(!s_pad(K, z, s_inlen(strlen(str), radix))) 1983 return MP_MEMORY; 1984 1985 MP_USED(z) = 1; z->digits[0] = 0; 1986 1987 while(*str != '\0' && ((ch = s_ch2val(*str, radix)) >= 0)) { 1988 s_dmul(z, (mp_digit)radix); 1989 s_dadd(z, (mp_digit)ch); 1990 ++str; 1991 } 1992 1993 CLAMP(z); 1994 1995 /* Override sign for zero, even if negative specified. */ 1996 if(CMPZ(z) == 0) 1997 MP_SIGN(z) = MP_ZPOS; 1998 1999 if(end != NULL) 2000 *end = (char *)str; 2001 2002 /* Return a truncation error if the string has unprocessed 2003 characters remaining, so the caller can tell if the whole string 2004 was done */ 2005 if(*str != '\0') 2006 return MP_TRUNC; 2007 else 2008 return MP_OK; 2009 } 2010 2011 /* }}} */ 2012 2013 /* {{{ mp_int_count_bits(z) */ 2014 2015 mp_result mp_int_count_bits(mp_int z) 2016 { 2017 mp_size nbits = 0, uz; 2018 mp_digit d; 2019 2020 CHECK(z != NULL); 2021 2022 uz = MP_USED(z); 2023 if(uz == 1 && z->digits[0] == 0) 2024 return 1; 2025 2026 --uz; 2027 nbits = uz * MP_DIGIT_BIT; 2028 d = z->digits[uz]; 2029 2030 while(d != 0) { 2031 d >>= 1; 2032 ++nbits; 2033 } 2034 2035 return nbits; 2036 } 2037 2038 /* }}} */ 2039 2040 /* {{{ mp_int_to_binary(z, buf, limit) */ 2041 2042 mp_result mp_int_to_binary(klisp_State *K, mp_int z, unsigned char *buf, 2043 int limit) 2044 { 2045 static const int PAD_FOR_2C = 1; 2046 2047 mp_result res; 2048 int limpos = limit; 2049 2050 CHECK(z != NULL && buf != NULL); 2051 2052 res = s_tobin(z, buf, &limpos, PAD_FOR_2C); 2053 2054 if(MP_SIGN(z) == MP_NEG) 2055 s_2comp(buf, limpos); 2056 2057 return res; 2058 } 2059 2060 /* }}} */ 2061 2062 /* {{{ mp_int_read_binary(z, buf, len) */ 2063 2064 mp_result mp_int_read_binary(klisp_State *K, mp_int z, unsigned char *buf, 2065 int len) 2066 { 2067 mp_size need, i; 2068 unsigned char *tmp; 2069 mp_digit *dz; 2070 2071 CHECK(z != NULL && buf != NULL && len > 0); 2072 2073 /* Figure out how many digits are needed to represent this value */ 2074 need = ((len * CHAR_BIT) + (MP_DIGIT_BIT - 1)) / MP_DIGIT_BIT; 2075 if(!s_pad(K, z, need)) 2076 return MP_MEMORY; 2077 2078 mp_int_zero(z); 2079 2080 /* If the high-order bit is set, take the 2's complement before 2081 reading the value (it will be restored afterward) */ 2082 if(buf[0] >> (CHAR_BIT - 1)) { 2083 MP_SIGN(z) = MP_NEG; 2084 s_2comp(buf, len); 2085 } 2086 2087 dz = MP_DIGITS(z); 2088 for(tmp = buf, i = len; i > 0; --i, ++tmp) { 2089 s_qmul(K, z, (mp_size) CHAR_BIT); 2090 *dz |= *tmp; 2091 } 2092 2093 /* Restore 2's complement if we took it before */ 2094 if(MP_SIGN(z) == MP_NEG) 2095 s_2comp(buf, len); 2096 2097 return MP_OK; 2098 } 2099 2100 /* }}} */ 2101 2102 /* {{{ mp_int_binary_len(z) */ 2103 2104 mp_result mp_int_binary_len(mp_int z) 2105 { 2106 mp_result res = mp_int_count_bits(z); 2107 int bytes = mp_int_unsigned_len(z); 2108 2109 if(res <= 0) 2110 return res; 2111 2112 bytes = (res + (CHAR_BIT - 1)) / CHAR_BIT; 2113 2114 /* If the highest-order bit falls exactly on a byte boundary, we 2115 need to pad with an extra byte so that the sign will be read 2116 correctly when reading it back in. */ 2117 if(bytes * CHAR_BIT == res) 2118 ++bytes; 2119 2120 return bytes; 2121 } 2122 2123 /* }}} */ 2124 2125 /* {{{ mp_int_to_unsigned(z, buf, limit) */ 2126 2127 mp_result mp_int_to_unsigned(klisp_State *K, mp_int z, unsigned char *buf, 2128 int limit) 2129 { 2130 static const int NO_PADDING = 0; 2131 2132 CHECK(z != NULL && buf != NULL); 2133 2134 return s_tobin(z, buf, &limit, NO_PADDING); 2135 } 2136 2137 /* }}} */ 2138 2139 /* {{{ mp_int_read_unsigned(z, buf, len) */ 2140 2141 mp_result mp_int_read_unsigned(klisp_State *K, mp_int z, unsigned char *buf, int len) 2142 { 2143 mp_size need, i; 2144 unsigned char *tmp; 2145 mp_digit *dz; 2146 2147 CHECK(z != NULL && buf != NULL && len > 0); 2148 2149 /* Figure out how many digits are needed to represent this value */ 2150 need = ((len * CHAR_BIT) + (MP_DIGIT_BIT - 1)) / MP_DIGIT_BIT; 2151 if(!s_pad(K, z, need)) 2152 return MP_MEMORY; 2153 2154 mp_int_zero(z); 2155 2156 dz = MP_DIGITS(z); 2157 for(tmp = buf, i = len; i > 0; --i, ++tmp) { 2158 (void) s_qmul(K, z, CHAR_BIT); 2159 *dz |= *tmp; 2160 } 2161 2162 return MP_OK; 2163 } 2164 2165 /* }}} */ 2166 2167 /* {{{ mp_int_unsigned_len(z) */ 2168 2169 mp_result mp_int_unsigned_len(mp_int z) 2170 { 2171 mp_result res = mp_int_count_bits(z); 2172 int bytes; 2173 2174 if(res <= 0) 2175 return res; 2176 2177 bytes = (res + (CHAR_BIT - 1)) / CHAR_BIT; 2178 2179 return bytes; 2180 } 2181 2182 /* }}} */ 2183 2184 /* {{{ mp_error_string(res) */ 2185 2186 const char *mp_error_string(mp_result res) 2187 { 2188 int ix; 2189 if(res > 0) 2190 return s_unknown_err; 2191 2192 res = -res; 2193 for(ix = 0; ix < res && s_error_msg[ix] != NULL; ++ix) 2194 ; 2195 2196 if(s_error_msg[ix] != NULL) 2197 return s_error_msg[ix]; 2198 else 2199 return s_unknown_err; 2200 } 2201 2202 /* }}} */ 2203 2204 /*------------------------------------------------------------------------*/ 2205 /* Private functions for internal use. These make assumptions. */ 2206 2207 /* {{{ s_alloc(num) */ 2208 2209 STATIC mp_digit *s_alloc(klisp_State *K, mp_size num) 2210 { 2211 mp_digit *out = klispM_malloc(K, num * sizeof(mp_digit)); 2212 2213 assert(out != NULL); /* for debugging */ 2214 #if DEBUG > 1 2215 { 2216 mp_digit v = (mp_digit) 0xdeadbeef; 2217 int ix; 2218 2219 for(ix = 0; ix < num; ++ix) 2220 out[ix] = v; 2221 } 2222 #endif 2223 2224 return out; 2225 } 2226 2227 /* }}} */ 2228 2229 /* {{{ s_realloc(old, osize, nsize) */ 2230 2231 STATIC mp_digit *s_realloc(klisp_State *K, mp_digit *old, mp_size osize, 2232 mp_size nsize) 2233 { 2234 #if DEBUG > 1 2235 mp_digit *new = s_alloc(K, nsize); 2236 int ix; 2237 2238 for(ix = 0; ix < nsize; ++ix) 2239 new[ix] = (mp_digit) 0xdeadbeef; 2240 2241 memcpy(new, old, osize * sizeof(mp_digit)); 2242 #else 2243 mp_digit *new = klispM_realloc_(K, old, osize * sizeof(mp_digit), 2244 nsize * sizeof(mp_digit)); 2245 assert(new != NULL); /* for debugging */ 2246 #endif 2247 return new; 2248 } 2249 2250 /* }}} */ 2251 2252 /* {{{ s_free(ptr) */ 2253 2254 STATIC void s_free(klisp_State *K, void *ptr, mp_size size) 2255 { 2256 klispM_freemem(K, ptr, size * sizeof(mp_digit)); 2257 } 2258 2259 /* }}} */ 2260 2261 /* {{{ s_pad(z, min) */ 2262 2263 STATIC int s_pad(klisp_State *K, mp_int z, mp_size min) 2264 { 2265 if(MP_ALLOC(z) < min) { 2266 mp_size nsize = ROUND_PREC(min); 2267 mp_digit *tmp; 2268 2269 if((void *)z->digits == (void *)&(z->single)) { 2270 if((tmp = s_alloc(K, nsize)) == NULL) 2271 return 0; 2272 2273 COPY(MP_DIGITS(z), tmp, MP_USED(z)); 2274 } else if((tmp = s_realloc(K, MP_DIGITS(z), MP_ALLOC(z), nsize)) == NULL) 2275 return 0; 2276 2277 MP_DIGITS(z) = tmp; 2278 MP_ALLOC(z) = nsize; 2279 } 2280 2281 return 1; 2282 } 2283 2284 /* }}} */ 2285 2286 /* {{{ s_fake(z, value, vbuf) */ 2287 2288 STATIC void s_fake(mp_int z, mp_small value, mp_digit vbuf[]) 2289 { 2290 mp_size uv = (mp_size) s_vpack(value, vbuf); 2291 2292 z->used = uv; 2293 z->alloc = MP_VALUE_DIGITS(value); 2294 z->sign = (value < 0) ? MP_NEG : MP_ZPOS; 2295 z->digits = vbuf; 2296 } 2297 2298 /* }}} */ 2299 2300 /* {{{ s_cdig(da, db, len) */ 2301 2302 STATIC int s_cdig(mp_digit *da, mp_digit *db, mp_size len) 2303 { 2304 mp_digit *dat = da + len - 1, *dbt = db + len - 1; 2305 2306 for(/* */; len != 0; --len, --dat, --dbt) { 2307 if(*dat > *dbt) 2308 return 1; 2309 else if(*dat < *dbt) 2310 return -1; 2311 } 2312 2313 return 0; 2314 } 2315 2316 /* }}} */ 2317 2318 /* {{{ s_vpack(v, t[]) */ 2319 2320 STATIC int s_vpack(mp_small v, mp_digit t[]) 2321 { 2322 mp_usmall uv = (mp_usmall) ((v < 0) ? -v : v); 2323 int ndig = 0; 2324 2325 if(uv == 0) 2326 t[ndig++] = 0; 2327 else { 2328 while(uv != 0) { 2329 t[ndig++] = (mp_digit) uv; 2330 uv >>= MP_DIGIT_BIT/2; 2331 uv >>= MP_DIGIT_BIT/2; 2332 } 2333 } 2334 2335 return ndig; 2336 } 2337 2338 /* }}} */ 2339 2340 /* {{{ s_ucmp(a, b) */ 2341 2342 STATIC int s_ucmp(mp_int a, mp_int b) 2343 { 2344 mp_size ua = MP_USED(a), ub = MP_USED(b); 2345 2346 if(ua > ub) 2347 return 1; 2348 else if(ub > ua) 2349 return -1; 2350 else 2351 return s_cdig(MP_DIGITS(a), MP_DIGITS(b), ua); 2352 } 2353 2354 /* }}} */ 2355 2356 /* {{{ s_vcmp(a, v) */ 2357 2358 STATIC int s_vcmp(mp_int a, mp_small v) 2359 { 2360 mp_digit vdig[MP_VALUE_DIGITS(v)]; 2361 int ndig = 0; 2362 mp_size ua = MP_USED(a); 2363 2364 ndig = s_vpack(v, vdig); 2365 2366 if(ua > ndig) 2367 return 1; 2368 else if(ua < ndig) 2369 return -1; 2370 else 2371 return s_cdig(MP_DIGITS(a), vdig, ndig); 2372 } 2373 2374 /* }}} */ 2375 2376 /* {{{ s_uadd(da, db, dc, size_a, size_b) */ 2377 2378 STATIC mp_digit s_uadd(mp_digit *da, mp_digit *db, mp_digit *dc, 2379 mp_size size_a, mp_size size_b) 2380 { 2381 mp_size pos; 2382 mp_word w = 0; 2383 2384 /* Insure that da is the longer of the two to simplify later code */ 2385 if(size_b > size_a) { 2386 SWAP(mp_digit *, da, db); 2387 SWAP(mp_size, size_a, size_b); 2388 } 2389 2390 /* Add corresponding digits until the shorter number runs out */ 2391 for(pos = 0; pos < size_b; ++pos, ++da, ++db, ++dc) { 2392 w = w + (mp_word) *da + (mp_word) *db; 2393 *dc = LOWER_HALF(w); 2394 w = UPPER_HALF(w); 2395 } 2396 2397 /* Propagate carries as far as necessary */ 2398 for(/* */; pos < size_a; ++pos, ++da, ++dc) { 2399 w = w + *da; 2400 2401 *dc = LOWER_HALF(w); 2402 w = UPPER_HALF(w); 2403 } 2404 2405 /* Return carry out */ 2406 return (mp_digit)w; 2407 } 2408 2409 /* }}} */ 2410 2411 /* {{{ s_usub(da, db, dc, size_a, size_b) */ 2412 2413 STATIC void s_usub(mp_digit *da, mp_digit *db, mp_digit *dc, 2414 mp_size size_a, mp_size size_b) 2415 { 2416 mp_size pos; 2417 mp_word w = 0; 2418 2419 /* We assume that |a| >= |b| so this should definitely hold */ 2420 assert(size_a >= size_b); 2421 2422 /* Subtract corresponding digits and propagate borrow */ 2423 for(pos = 0; pos < size_b; ++pos, ++da, ++db, ++dc) { 2424 w = ((mp_word)MP_DIGIT_MAX + 1 + /* MP_RADIX */ 2425 (mp_word)*da) - w - (mp_word)*db; 2426 2427 *dc = LOWER_HALF(w); 2428 w = (UPPER_HALF(w) == 0); 2429 } 2430 2431 /* Finish the subtraction for remaining upper digits of da */ 2432 for(/* */; pos < size_a; ++pos, ++da, ++dc) { 2433 w = ((mp_word)MP_DIGIT_MAX + 1 + /* MP_RADIX */ 2434 (mp_word)*da) - w; 2435 2436 *dc = LOWER_HALF(w); 2437 w = (UPPER_HALF(w) == 0); 2438 } 2439 2440 /* If there is a borrow out at the end, it violates the precondition */ 2441 assert(w == 0); 2442 } 2443 2444 /* }}} */ 2445 2446 /* {{{ s_kmul(da, db, dc, size_a, size_b) */ 2447 2448 STATIC int s_kmul(klisp_State *K, mp_digit *da, mp_digit *db, 2449 mp_digit *dc, mp_size size_a, mp_size size_b) 2450 { 2451 mp_size bot_size; 2452 2453 /* Make sure b is the smaller of the two input values */ 2454 if(size_b > size_a) { 2455 SWAP(mp_digit *, da, db); 2456 SWAP(mp_size, size_a, size_b); 2457 } 2458 2459 /* Insure that the bottom is the larger half in an odd-length split; 2460 the code below relies on this being true. 2461 */ 2462 bot_size = (size_a + 1) / 2; 2463 2464 /* If the values are big enough to bother with recursion, use the 2465 Karatsuba algorithm to compute the product; otherwise use the 2466 normal multiplication algorithm 2467 */ 2468 if(multiply_threshold && 2469 size_a >= multiply_threshold && 2470 size_b > bot_size) { 2471 2472 mp_digit *t1, *t2, *t3, carry; 2473 2474 mp_digit *a_top = da + bot_size; 2475 mp_digit *b_top = db + bot_size; 2476 2477 mp_size at_size = size_a - bot_size; 2478 mp_size bt_size = size_b - bot_size; 2479 mp_size buf_size = 2 * bot_size; 2480 2481 /* Do a single allocation for all three temporary buffers needed; 2482 each buffer must be big enough to hold the product of two 2483 bottom halves, and one buffer needs space for the completed 2484 product; twice the space is plenty. 2485 */ 2486 if((t1 = s_alloc(K, 4 * buf_size)) == NULL) return 0; 2487 t2 = t1 + buf_size; 2488 t3 = t2 + buf_size; 2489 ZERO(t1, 4 * buf_size); 2490 2491 /* t1 and t2 are initially used as temporaries to compute the inner product 2492 (a1 + a0)(b1 + b0) = a1b1 + a1b0 + a0b1 + a0b0 2493 */ 2494 carry = s_uadd(da, a_top, t1, bot_size, at_size); /* t1 = a1 + a0 */ 2495 t1[bot_size] = carry; 2496 2497 carry = s_uadd(db, b_top, t2, bot_size, bt_size); /* t2 = b1 + b0 */ 2498 t2[bot_size] = carry; 2499 2500 /* t3 = t1 * t2 */ 2501 (void) s_kmul(K, t1, t2, t3, bot_size + 1, bot_size + 1); 2502 2503 /* Now we'll get t1 = a0b0 and t2 = a1b1, and subtract them out so that 2504 we're left with only the pieces we want: t3 = a1b0 + a0b1 2505 */ 2506 ZERO(t1, buf_size); 2507 ZERO(t2, buf_size); 2508 (void) s_kmul(K, da, db, t1, bot_size, bot_size); /* t1 = a0 * b0 */ 2509 (void) s_kmul(K, a_top, b_top, t2, at_size, bt_size); /* t2 = a1 * b1 */ 2510 2511 /* Subtract out t1 and t2 to get the inner product */ 2512 s_usub(t3, t1, t3, buf_size + 2, buf_size); 2513 s_usub(t3, t2, t3, buf_size + 2, buf_size); 2514 2515 /* Assemble the output value */ 2516 COPY(t1, dc, buf_size); 2517 carry = s_uadd(t3, dc + bot_size, dc + bot_size, 2518 buf_size + 1, buf_size); 2519 assert(carry == 0); 2520 2521 carry = s_uadd(t2, dc + 2*bot_size, dc + 2*bot_size, 2522 buf_size, buf_size); 2523 assert(carry == 0); 2524 2525 /* note t2 and t3 are just internal pointers to t1 */ 2526 s_free(K, t1, 4 * buf_size); 2527 } 2528 else { 2529 s_umul(da, db, dc, size_a, size_b); 2530 } 2531 2532 return 1; 2533 } 2534 2535 /* }}} */ 2536 2537 /* {{{ s_umul(da, db, dc, size_a, size_b) */ 2538 2539 STATIC void s_umul(mp_digit *da, mp_digit *db, mp_digit *dc, 2540 mp_size size_a, mp_size size_b) 2541 { 2542 mp_size a, b; 2543 mp_word w; 2544 2545 for(a = 0; a < size_a; ++a, ++dc, ++da) { 2546 mp_digit *dct = dc; 2547 mp_digit *dbt = db; 2548 2549 if(*da == 0) 2550 continue; 2551 2552 w = 0; 2553 for(b = 0; b < size_b; ++b, ++dbt, ++dct) { 2554 w = (mp_word)*da * (mp_word)*dbt + w + (mp_word)*dct; 2555 2556 *dct = LOWER_HALF(w); 2557 w = UPPER_HALF(w); 2558 } 2559 2560 *dct = (mp_digit)w; 2561 } 2562 } 2563 2564 /* }}} */ 2565 2566 /* {{{ s_ksqr(da, dc, size_a) */ 2567 2568 STATIC int s_ksqr(klisp_State *K, mp_digit *da, mp_digit *dc, 2569 mp_size size_a) 2570 { 2571 if(multiply_threshold && size_a > multiply_threshold) { 2572 mp_size bot_size = (size_a + 1) / 2; 2573 mp_digit *a_top = da + bot_size; 2574 mp_digit *t1, *t2, *t3, carry; 2575 mp_size at_size = size_a - bot_size; 2576 mp_size buf_size = 2 * bot_size; 2577 2578 if((t1 = s_alloc(K, 4 * buf_size)) == NULL) return 0; 2579 t2 = t1 + buf_size; 2580 t3 = t2 + buf_size; 2581 ZERO(t1, 4 * buf_size); 2582 2583 (void) s_ksqr(K, da, t1, bot_size); /* t1 = a0 ^ 2 */ 2584 (void) s_ksqr(K, a_top, t2, at_size); /* t2 = a1 ^ 2 */ 2585 2586 (void) s_kmul(K, da, a_top, t3, bot_size, at_size); /* t3 = a0 * a1 */ 2587 2588 /* Quick multiply t3 by 2, shifting left (can't overflow) */ 2589 { 2590 int i, top = bot_size + at_size; 2591 mp_word w, save = 0; 2592 2593 for(i = 0; i < top; ++i) { 2594 w = t3[i]; 2595 w = (w << 1) | save; 2596 t3[i] = LOWER_HALF(w); 2597 save = UPPER_HALF(w); 2598 } 2599 t3[i] = LOWER_HALF(save); 2600 } 2601 2602 /* Assemble the output value */ 2603 COPY(t1, dc, 2 * bot_size); 2604 carry = s_uadd(t3, dc + bot_size, dc + bot_size, 2605 buf_size + 1, buf_size); 2606 assert(carry == 0); 2607 2608 carry = s_uadd(t2, dc + 2*bot_size, dc + 2*bot_size, 2609 buf_size, buf_size); 2610 assert(carry == 0); 2611 2612 /* note that t2 and t2 are internal pointers only */ 2613 s_free(K, t1, 4 * buf_size); 2614 } 2615 else { 2616 s_usqr(da, dc, size_a); 2617 } 2618 2619 return 1; 2620 } 2621 2622 /* }}} */ 2623 2624 /* {{{ s_usqr(da, dc, size_a) */ 2625 2626 STATIC void s_usqr(mp_digit *da, mp_digit *dc, mp_size size_a) 2627 { 2628 mp_size i, j; 2629 mp_word w; 2630 2631 for(i = 0; i < size_a; ++i, dc += 2, ++da) { 2632 mp_digit *dct = dc, *dat = da; 2633 2634 if(*da == 0) 2635 continue; 2636 2637 /* Take care of the first digit, no rollover */ 2638 w = (mp_word)*dat * (mp_word)*dat + (mp_word)*dct; 2639 *dct = LOWER_HALF(w); 2640 w = UPPER_HALF(w); 2641 ++dat; ++dct; 2642 2643 for(j = i + 1; j < size_a; ++j, ++dat, ++dct) { 2644 mp_word t = (mp_word)*da * (mp_word)*dat; 2645 mp_word u = w + (mp_word)*dct, ov = 0; 2646 2647 /* Check if doubling t will overflow a word */ 2648 if(HIGH_BIT_SET(t)) 2649 ov = 1; 2650 2651 w = t + t; 2652 2653 /* Check if adding u to w will overflow a word */ 2654 if(ADD_WILL_OVERFLOW(w, u)) 2655 ov = 1; 2656 2657 w += u; 2658 2659 *dct = LOWER_HALF(w); 2660 w = UPPER_HALF(w); 2661 if(ov) { 2662 w += MP_DIGIT_MAX; /* MP_RADIX */ 2663 ++w; 2664 } 2665 } 2666 2667 w = w + *dct; 2668 *dct = (mp_digit)w; 2669 while((w = UPPER_HALF(w)) != 0) { 2670 ++dct; w = w + *dct; 2671 *dct = LOWER_HALF(w); 2672 } 2673 2674 assert(w == 0); 2675 } 2676 } 2677 2678 /* }}} */ 2679 2680 /* {{{ s_dadd(a, b) */ 2681 2682 STATIC void s_dadd(mp_int a, mp_digit b) 2683 { 2684 mp_word w = 0; 2685 mp_digit *da = MP_DIGITS(a); 2686 mp_size ua = MP_USED(a); 2687 2688 w = (mp_word)*da + b; 2689 *da++ = LOWER_HALF(w); 2690 w = UPPER_HALF(w); 2691 2692 for(ua -= 1; ua > 0; --ua, ++da) { 2693 w = (mp_word)*da + w; 2694 2695 *da = LOWER_HALF(w); 2696 w = UPPER_HALF(w); 2697 } 2698 2699 if(w) { 2700 *da = (mp_digit)w; 2701 MP_USED(a) += 1; 2702 } 2703 } 2704 2705 /* }}} */ 2706 2707 /* {{{ s_dmul(a, b) */ 2708 2709 STATIC void s_dmul(mp_int a, mp_digit b) 2710 { 2711 mp_word w = 0; 2712 mp_digit *da = MP_DIGITS(a); 2713 mp_size ua = MP_USED(a); 2714 2715 while(ua > 0) { 2716 w = (mp_word)*da * b + w; 2717 *da++ = LOWER_HALF(w); 2718 w = UPPER_HALF(w); 2719 --ua; 2720 } 2721 2722 if(w) { 2723 *da = (mp_digit)w; 2724 MP_USED(a) += 1; 2725 } 2726 } 2727 2728 /* }}} */ 2729 2730 /* {{{ s_dbmul(da, b, dc, size_a) */ 2731 2732 STATIC void s_dbmul(mp_digit *da, mp_digit b, mp_digit *dc, mp_size size_a) 2733 { 2734 mp_word w = 0; 2735 2736 while(size_a > 0) { 2737 w = (mp_word)*da++ * (mp_word)b + w; 2738 2739 *dc++ = LOWER_HALF(w); 2740 w = UPPER_HALF(w); 2741 --size_a; 2742 } 2743 2744 if(w) 2745 *dc = LOWER_HALF(w); 2746 } 2747 2748 /* }}} */ 2749 2750 /* {{{ s_ddiv(da, d, dc, size_a) */ 2751 2752 STATIC mp_digit s_ddiv(mp_int a, mp_digit b) 2753 { 2754 mp_word w = 0, qdigit; 2755 mp_size ua = MP_USED(a); 2756 mp_digit *da = MP_DIGITS(a) + ua - 1; 2757 2758 for(/* */; ua > 0; --ua, --da) { 2759 w = (w << MP_DIGIT_BIT) | *da; 2760 2761 if(w >= b) { 2762 qdigit = w / b; 2763 w = w % b; 2764 } 2765 else { 2766 qdigit = 0; 2767 } 2768 2769 *da = (mp_digit)qdigit; 2770 } 2771 2772 CLAMP(a); 2773 return (mp_digit)w; 2774 } 2775 2776 /* }}} */ 2777 2778 /* {{{ s_qdiv(z, p2) */ 2779 2780 STATIC void s_qdiv(mp_int z, mp_size p2) 2781 { 2782 mp_size ndig = p2 / MP_DIGIT_BIT, nbits = p2 % MP_DIGIT_BIT; 2783 mp_size uz = MP_USED(z); 2784 2785 if(ndig) { 2786 mp_size mark; 2787 mp_digit *to, *from; 2788 2789 if(ndig >= uz) { 2790 mp_int_zero(z); 2791 return; 2792 } 2793 2794 to = MP_DIGITS(z); from = to + ndig; 2795 2796 for(mark = ndig; mark < uz; ++mark) 2797 *to++ = *from++; 2798 2799 MP_USED(z) = uz - ndig; 2800 } 2801 2802 if(nbits) { 2803 mp_digit d = 0, *dz, save; 2804 mp_size up = MP_DIGIT_BIT - nbits; 2805 2806 uz = MP_USED(z); 2807 dz = MP_DIGITS(z) + uz - 1; 2808 2809 for(/* */; uz > 0; --uz, --dz) { 2810 save = *dz; 2811 2812 *dz = (*dz >> nbits) | (d << up); 2813 d = save; 2814 } 2815 2816 CLAMP(z); 2817 } 2818 2819 if(MP_USED(z) == 1 && z->digits[0] == 0) 2820 MP_SIGN(z) = MP_ZPOS; 2821 } 2822 2823 /* }}} */ 2824 2825 /* {{{ s_qmod(z, p2) */ 2826 2827 STATIC void s_qmod(mp_int z, mp_size p2) 2828 { 2829 mp_size start = p2 / MP_DIGIT_BIT + 1, rest = p2 % MP_DIGIT_BIT; 2830 mp_size uz = MP_USED(z); 2831 mp_digit mask = (1 << rest) - 1; 2832 2833 if(start <= uz) { 2834 MP_USED(z) = start; 2835 z->digits[start - 1] &= mask; 2836 CLAMP(z); 2837 } 2838 } 2839 2840 /* }}} */ 2841 2842 /* {{{ s_qmul(z, p2) */ 2843 2844 STATIC int s_qmul(klisp_State *K, mp_int z, mp_size p2) 2845 { 2846 mp_size uz, need, rest, extra, i; 2847 mp_digit *from, *to, d; 2848 2849 if(p2 == 0) 2850 return 1; 2851 2852 uz = MP_USED(z); 2853 need = p2 / MP_DIGIT_BIT; rest = p2 % MP_DIGIT_BIT; 2854 2855 /* Figure out if we need an extra digit at the top end; this occurs 2856 if the topmost `rest' bits of the high-order digit of z are not 2857 zero, meaning they will be shifted off the end if not preserved */ 2858 extra = 0; 2859 if(rest != 0) { 2860 mp_digit *dz = MP_DIGITS(z) + uz - 1; 2861 2862 if((*dz >> (MP_DIGIT_BIT - rest)) != 0) 2863 extra = 1; 2864 } 2865 2866 if(!s_pad(K, z, uz + need + extra)) 2867 return 0; 2868 2869 /* If we need to shift by whole digits, do that in one pass, then 2870 to back and shift by partial digits. 2871 */ 2872 if(need > 0) { 2873 from = MP_DIGITS(z) + uz - 1; 2874 to = from + need; 2875 2876 for(i = 0; i < uz; ++i) 2877 *to-- = *from--; 2878 2879 ZERO(MP_DIGITS(z), need); 2880 uz += need; 2881 } 2882 2883 if(rest) { 2884 d = 0; 2885 for(i = need, from = MP_DIGITS(z) + need; i < uz; ++i, ++from) { 2886 mp_digit save = *from; 2887 2888 *from = (*from << rest) | (d >> (MP_DIGIT_BIT - rest)); 2889 d = save; 2890 } 2891 2892 d >>= (MP_DIGIT_BIT - rest); 2893 if(d != 0) { 2894 *from = d; 2895 uz += extra; 2896 } 2897 } 2898 2899 MP_USED(z) = uz; 2900 CLAMP(z); 2901 2902 return 1; 2903 } 2904 2905 /* }}} */ 2906 2907 /* {{{ s_qsub(z, p2) */ 2908 2909 /* Compute z = 2^p2 - |z|; requires that 2^p2 >= |z| 2910 The sign of the result is always zero/positive. 2911 */ 2912 STATIC int s_qsub(klisp_State *K, mp_int z, mp_size p2) 2913 { 2914 mp_digit hi = (1 << (p2 % MP_DIGIT_BIT)), *zp; 2915 mp_size tdig = (p2 / MP_DIGIT_BIT), pos; 2916 mp_word w = 0; 2917 2918 if(!s_pad(K, z, tdig + 1)) 2919 return 0; 2920 2921 for(pos = 0, zp = MP_DIGITS(z); pos < tdig; ++pos, ++zp) { 2922 w = ((mp_word) MP_DIGIT_MAX + 1) - w - (mp_word)*zp; 2923 2924 *zp = LOWER_HALF(w); 2925 w = UPPER_HALF(w) ? 0 : 1; 2926 } 2927 2928 w = ((mp_word) MP_DIGIT_MAX + 1 + hi) - w - (mp_word)*zp; 2929 *zp = LOWER_HALF(w); 2930 2931 assert(UPPER_HALF(w) != 0); /* no borrow out should be possible */ 2932 2933 MP_SIGN(z) = MP_ZPOS; 2934 CLAMP(z); 2935 2936 return 1; 2937 } 2938 2939 /* }}} */ 2940 2941 /* {{{ s_dp2k(z) */ 2942 2943 STATIC int s_dp2k(mp_int z) 2944 { 2945 int k = 0; 2946 mp_digit *dp = MP_DIGITS(z), d; 2947 2948 if(MP_USED(z) == 1 && *dp == 0) 2949 return 1; 2950 2951 while(*dp == 0) { 2952 k += MP_DIGIT_BIT; 2953 ++dp; 2954 } 2955 2956 d = *dp; 2957 while((d & 1) == 0) { 2958 d >>= 1; 2959 ++k; 2960 } 2961 2962 return k; 2963 } 2964 2965 /* }}} */ 2966 2967 /* {{{ s_isp2(z) */ 2968 2969 STATIC int s_isp2(mp_int z) 2970 { 2971 mp_size uz = MP_USED(z), k = 0; 2972 mp_digit *dz = MP_DIGITS(z), d; 2973 2974 while(uz > 1) { 2975 if(*dz++ != 0) 2976 return -1; 2977 k += MP_DIGIT_BIT; 2978 --uz; 2979 } 2980 2981 d = *dz; 2982 while(d > 1) { 2983 if(d & 1) 2984 return -1; 2985 ++k; d >>= 1; 2986 } 2987 2988 return (int) k; 2989 } 2990 2991 /* }}} */ 2992 2993 /* {{{ s_2expt(z, k) */ 2994 2995 STATIC int s_2expt(klisp_State *K, mp_int z, mp_small k) 2996 { 2997 mp_size ndig, rest; 2998 mp_digit *dz; 2999 3000 ndig = (k + MP_DIGIT_BIT) / MP_DIGIT_BIT; 3001 rest = k % MP_DIGIT_BIT; 3002 3003 if(!s_pad(K, z, ndig)) 3004 return 0; 3005 3006 dz = MP_DIGITS(z); 3007 ZERO(dz, ndig); 3008 *(dz + ndig - 1) = (1 << rest); 3009 MP_USED(z) = ndig; 3010 3011 return 1; 3012 } 3013 3014 /* }}} */ 3015 3016 /* {{{ s_norm(a, b) */ 3017 3018 STATIC int s_norm(klisp_State *K, mp_int a, mp_int b) 3019 { 3020 mp_digit d = b->digits[MP_USED(b) - 1]; 3021 int k = 0; 3022 3023 while(d < (mp_digit) (1 << (MP_DIGIT_BIT - 1))) { /* d < (MP_RADIX / 2) */ 3024 d <<= 1; 3025 ++k; 3026 } 3027 3028 /* These multiplications can't fail */ 3029 if(k != 0) { 3030 (void) s_qmul(K, a, (mp_size) k); 3031 (void) s_qmul(K, b, (mp_size) k); 3032 } 3033 3034 return k; 3035 } 3036 3037 /* }}} */ 3038 3039 /* {{{ s_brmu(z, m) */ 3040 3041 STATIC mp_result s_brmu(klisp_State *K, mp_int z, mp_int m) 3042 { 3043 mp_size um = MP_USED(m) * 2; 3044 3045 if(!s_pad(K, z, um)) 3046 return MP_MEMORY; 3047 3048 s_2expt(K, z, MP_DIGIT_BIT * um); 3049 return mp_int_div(K, z, m, z, NULL); 3050 } 3051 3052 /* }}} */ 3053 3054 /* {{{ s_reduce(x, m, mu, q1, q2) */ 3055 3056 STATIC int s_reduce(klisp_State *K, mp_int x, mp_int m, mp_int mu, 3057 mp_int q1, mp_int q2) 3058 { 3059 mp_size um = MP_USED(m), umb_p1, umb_m1; 3060 3061 umb_p1 = (um + 1) * MP_DIGIT_BIT; 3062 umb_m1 = (um - 1) * MP_DIGIT_BIT; 3063 3064 if(mp_int_copy(K, x, q1) != MP_OK) 3065 return 0; 3066 3067 /* Compute q2 = floor((floor(x / b^(k-1)) * mu) / b^(k+1)) */ 3068 s_qdiv(q1, umb_m1); 3069 UMUL(K, q1, mu, q2); 3070 s_qdiv(q2, umb_p1); 3071 3072 /* Set x = x mod b^(k+1) */ 3073 s_qmod(x, umb_p1); 3074 3075 /* Now, q is a guess for the quotient a / m. 3076 Compute x - q * m mod b^(k+1), replacing x. This may be off 3077 by a factor of 2m, but no more than that. 3078 */ 3079 UMUL(K, q2, m, q1); 3080 s_qmod(q1, umb_p1); 3081 (void) mp_int_sub(K, x, q1, x); /* can't fail */ 3082 3083 /* The result may be < 0; if it is, add b^(k+1) to pin it in the 3084 proper range. */ 3085 if((CMPZ(x) < 0) && !s_qsub(K, x, umb_p1)) 3086 return 0; 3087 3088 /* If x > m, we need to back it off until it is in range. 3089 This will be required at most twice. */ 3090 if(mp_int_compare(x, m) >= 0) { 3091 (void) mp_int_sub(K, x, m, x); 3092 if(mp_int_compare(x, m) >= 0) 3093 (void) mp_int_sub(K, x, m, x); 3094 } 3095 3096 /* At this point, x has been properly reduced. */ 3097 return 1; 3098 } 3099 3100 /* }}} */ 3101 3102 /* {{{ s_embar(a, b, m, mu, c) */ 3103 3104 /* Perform modular exponentiation using Barrett's method, where mu is 3105 the reduction constant for m. Assumes a < m, b > 0. */ 3106 STATIC mp_result s_embar(klisp_State *K, mp_int a, mp_int b, mp_int m, 3107 mp_int mu, mp_int c) 3108 { 3109 mp_digit *db, *dbt, umu, d; 3110 mpz_t temp[3]; 3111 mp_result res; 3112 int last = 0; 3113 3114 umu = MP_USED(mu); db = MP_DIGITS(b); dbt = db + MP_USED(b) - 1; 3115 3116 while(last < 3) { 3117 SETUP(mp_int_init_size(K, TEMP(last), 4 * umu), last); 3118 ZERO(MP_DIGITS(TEMP(last - 1)), MP_ALLOC(TEMP(last - 1))); 3119 } 3120 3121 (void) mp_int_set_value(K, c, 1); 3122 3123 /* Take care of low-order digits */ 3124 while(db < dbt) { 3125 int i; 3126 3127 for(d = *db, i = MP_DIGIT_BIT; i > 0; --i, d >>= 1) { 3128 if(d & 1) { 3129 /* The use of a second temporary avoids allocation */ 3130 UMUL(K, c, a, TEMP(0)); 3131 if(!s_reduce(K, TEMP(0), m, mu, TEMP(1), TEMP(2))) { 3132 res = MP_MEMORY; goto CLEANUP; 3133 } 3134 mp_int_copy(K, TEMP(0), c); 3135 } 3136 3137 3138 USQR(K, a, TEMP(0)); 3139 assert(MP_SIGN(TEMP(0)) == MP_ZPOS); 3140 if(!s_reduce(K, TEMP(0), m, mu, TEMP(1), TEMP(2))) { 3141 res = MP_MEMORY; goto CLEANUP; 3142 } 3143 assert(MP_SIGN(TEMP(0)) == MP_ZPOS); 3144 mp_int_copy(K, TEMP(0), a); 3145 3146 3147 } 3148 3149 ++db; 3150 } 3151 3152 /* Take care of highest-order digit */ 3153 d = *dbt; 3154 for(;;) { 3155 if(d & 1) { 3156 UMUL(K, c, a, TEMP(0)); 3157 if(!s_reduce(K, TEMP(0), m, mu, TEMP(1), TEMP(2))) { 3158 res = MP_MEMORY; goto CLEANUP; 3159 } 3160 mp_int_copy(K, TEMP(0), c); 3161 } 3162 3163 d >>= 1; 3164 if(!d) break; 3165 3166 USQR(K, a, TEMP(0)); 3167 if(!s_reduce(K, TEMP(0), m, mu, TEMP(1), TEMP(2))) { 3168 res = MP_MEMORY; goto CLEANUP; 3169 } 3170 (void) mp_int_copy(K, TEMP(0), a); 3171 } 3172 3173 CLEANUP: 3174 while(--last >= 0) 3175 mp_int_clear(K, TEMP(last)); 3176 3177 return res; 3178 } 3179 3180 /* }}} */ 3181 3182 /* {{{ s_udiv(a, b) */ 3183 3184 /* Precondition: a >= b and b > 0 3185 Postcondition: a' = a / b, b' = a % b 3186 */ 3187 STATIC mp_result s_udiv(klisp_State *K, mp_int a, mp_int b) 3188 { 3189 mpz_t q, r, t; 3190 mp_size ua, ub, qpos = 0; 3191 mp_digit *da, btop; 3192 mp_result res = MP_OK; 3193 int k, skip = 0; 3194 3195 /* Force signs to positive */ 3196 MP_SIGN(a) = MP_ZPOS; 3197 MP_SIGN(b) = MP_ZPOS; 3198 3199 /* Normalize, per Knuth */ 3200 k = s_norm(K, a, b); 3201 3202 ua = MP_USED(a); ub = MP_USED(b); btop = b->digits[ub - 1]; 3203 if((res = mp_int_init_size(K, &q, ua)) != MP_OK) return res; 3204 if((res = mp_int_init_size(K, &t, ua + 1)) != MP_OK) goto CLEANUP; 3205 3206 da = MP_DIGITS(a); 3207 r.digits = da + ua - 1; /* The contents of r are shared with a */ 3208 r.used = 1; 3209 r.sign = MP_ZPOS; 3210 r.alloc = MP_ALLOC(a); 3211 ZERO(t.digits, t.alloc); 3212 3213 /* Solve for quotient digits, store in q.digits in reverse order */ 3214 while(r.digits >= da) { 3215 assert(qpos <= q.alloc); 3216 3217 if(s_ucmp(b, &r) > 0) { 3218 r.digits -= 1; 3219 r.used += 1; 3220 3221 if(++skip > 1 && qpos > 0) 3222 q.digits[qpos++] = 0; 3223 3224 CLAMP(&r); 3225 } 3226 else { 3227 mp_word pfx = r.digits[r.used - 1]; 3228 mp_word qdigit; 3229 3230 // Bugfix (was pfx <= btop in imath <= 1.17) Andres Navarro 3231 if(r.used > 1 && pfx < btop) { 3232 pfx <<= MP_DIGIT_BIT / 2; 3233 pfx <<= MP_DIGIT_BIT / 2; 3234 pfx |= r.digits[r.used - 2]; 3235 } 3236 3237 qdigit = pfx / btop; 3238 if(qdigit > MP_DIGIT_MAX) { 3239 qdigit = MP_DIGIT_MAX; 3240 } 3241 3242 s_dbmul(MP_DIGITS(b), (mp_digit) qdigit, t.digits, ub); 3243 t.used = ub + 1; CLAMP(&t); 3244 while(s_ucmp(&t, &r) > 0) { 3245 --qdigit; 3246 (void) mp_int_sub(K, &t, b, &t); /* cannot fail */ 3247 } 3248 3249 s_usub(r.digits, t.digits, r.digits, r.used, t.used); 3250 CLAMP(&r); 3251 3252 q.digits[qpos++] = (mp_digit) qdigit; 3253 ZERO(t.digits, t.used); 3254 skip = 0; 3255 } 3256 } 3257 3258 /* Put quotient digits in the correct order, and discard extra zeroes */ 3259 q.used = qpos; 3260 REV(mp_digit, q.digits, qpos); 3261 CLAMP(&q); 3262 3263 /* Denormalize the remainder */ 3264 CLAMP(a); 3265 if(k != 0) 3266 s_qdiv(a, k); 3267 3268 mp_int_copy(K, a, b); /* ok: 0 <= r < b */ 3269 mp_int_copy(K, &q, a); /* ok: q <= a */ 3270 3271 mp_int_clear(K, &t); 3272 CLEANUP: 3273 mp_int_clear(K, &q); 3274 return res; 3275 } 3276 3277 /* }}} */ 3278 3279 /* {{{ s_outlen(z, r) */ 3280 3281 STATIC int s_outlen(mp_int z, mp_size r) 3282 { 3283 mp_result bits; 3284 double raw; 3285 3286 assert(r >= MP_MIN_RADIX && r <= MP_MAX_RADIX); 3287 3288 bits = mp_int_count_bits(z); 3289 raw = (double)bits * s_log2[r]; 3290 3291 return (int)(raw + 0.999999); 3292 } 3293 3294 /* }}} */ 3295 3296 /* {{{ s_inlen(len, r) */ 3297 3298 STATIC mp_size s_inlen(int len, mp_size r) 3299 { 3300 double raw = (double)len / s_log2[r]; 3301 mp_size bits = (mp_size)(raw + 0.5); 3302 3303 return (mp_size)((bits + (MP_DIGIT_BIT - 1)) / MP_DIGIT_BIT); 3304 } 3305 3306 /* }}} */ 3307 3308 /* {{{ s_ch2val(c, r) */ 3309 3310 STATIC int s_ch2val(char c, int r) 3311 { 3312 int out; 3313 3314 if(isdigit((unsigned char) c)) 3315 out = c - '0'; 3316 else if(r > 10 && isalpha((unsigned char) c)) 3317 out = toupper(c) - 'A' + 10; 3318 else 3319 return -1; 3320 3321 return (out >= r) ? -1 : out; 3322 } 3323 3324 /* }}} */ 3325 3326 /* {{{ s_val2ch(v, caps) */ 3327 3328 STATIC char s_val2ch(int v, int caps) 3329 { 3330 assert(v >= 0); 3331 3332 if(v < 10) 3333 return v + '0'; 3334 else { 3335 char out = (v - 10) + 'a'; 3336 3337 if(caps) 3338 return toupper(out); 3339 else 3340 return out; 3341 } 3342 } 3343 3344 /* }}} */ 3345 3346 /* {{{ s_2comp(buf, len) */ 3347 3348 STATIC void s_2comp(unsigned char *buf, int len) 3349 { 3350 int i; 3351 unsigned short s = 1; 3352 3353 for(i = len - 1; i >= 0; --i) { 3354 unsigned char c = ~buf[i]; 3355 3356 s = c + s; 3357 c = s & UCHAR_MAX; 3358 s >>= CHAR_BIT; 3359 3360 buf[i] = c; 3361 } 3362 3363 /* last carry out is ignored */ 3364 } 3365 3366 /* }}} */ 3367 3368 /* {{{ s_tobin(z, buf, *limpos) */ 3369 3370 STATIC mp_result s_tobin(mp_int z, unsigned char *buf, int *limpos, int pad) 3371 { 3372 mp_size uz; 3373 mp_digit *dz; 3374 int pos = 0, limit = *limpos; 3375 3376 uz = MP_USED(z); dz = MP_DIGITS(z); 3377 while(uz > 0 && pos < limit) { 3378 mp_digit d = *dz++; 3379 int i; 3380 3381 for(i = sizeof(mp_digit); i > 0 && pos < limit; --i) { 3382 buf[pos++] = (unsigned char)d; 3383 d >>= CHAR_BIT; 3384 3385 /* Don't write leading zeroes */ 3386 if(d == 0 && uz == 1) 3387 i = 0; /* exit loop without signaling truncation */ 3388 } 3389 3390 /* Detect truncation (loop exited with pos >= limit) */ 3391 if(i > 0) break; 3392 3393 --uz; 3394 } 3395 3396 if(pad != 0 && (buf[pos - 1] >> (CHAR_BIT - 1))) { 3397 if(pos < limit) 3398 buf[pos++] = 0; 3399 else 3400 uz = 1; 3401 } 3402 3403 /* Digits are in reverse order, fix that */ 3404 REV(unsigned char, buf, pos); 3405 3406 /* Return the number of bytes actually written */ 3407 *limpos = pos; 3408 3409 return (uz == 0) ? MP_OK : MP_TRUNC; 3410 } 3411 3412 /* }}} */ 3413 3414 /* {{{ s_print(tag, z) */ 3415 3416 #if DEBUG 3417 void s_print(char *tag, mp_int z) 3418 { 3419 int i; 3420 3421 fprintf(stderr, "%s: %c ", tag, 3422 (MP_SIGN(z) == MP_NEG) ? '-' : '+'); 3423 3424 for(i = MP_USED(z) - 1; i >= 0; --i) 3425 fprintf(stderr, "%0*X", (int)(MP_DIGIT_BIT / 4), z->digits[i]); 3426 3427 fputc('\n', stderr); 3428 3429 } 3430 3431 void s_print_buf(char *tag, mp_digit *buf, mp_size num) 3432 { 3433 int i; 3434 3435 fprintf(stderr, "%s: ", tag); 3436 3437 for(i = num - 1; i >= 0; --i) 3438 fprintf(stderr, "%0*X", (int)(MP_DIGIT_BIT / 4), buf[i]); 3439 3440 fputc('\n', stderr); 3441 } 3442 #endif 3443 3444 /* }}} */ 3445 3446 /* HERE THERE BE DRAGONS */