1 /* mini-gmp, a minimalistic implementation of a GNU GMP subset. 2 3 Contributed to the GNU project by Niels Möller 4 5 Copyright 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1999, 2000, 2001, 6 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 7 Free Software Foundation, Inc. 8 9 This file is part of the GNU MP Library. 10 11 The GNU MP Library is free software; you can redistribute it and/or modify 12 it under the terms of the GNU Lesser General Public License as published by 13 the Free Software Foundation; either version 3 of the License, or (at your 14 option) any later version. 15 16 The GNU MP Library is distributed in the hope that it will be useful, but 17 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 18 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 19 License for more details. 20 21 You should have received a copy of the GNU Lesser General Public License 22 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ 23 24 /* NOTE: All functions in this file which are not declared in 25 mini-gmp.h are internal, and are not intended to be compatible 26 neither with GMP nor with future versions of mini-gmp. */ 27 28 /* Much of the material copied from GMP files, including: gmp-impl.h, 29 longlong.h, mpn/generic/add_n.c, mpn/generic/addmul_1.c, 30 mpn/generic/lshift.c, mpn/generic/mul_1.c, 31 mpn/generic/mul_basecase.c, mpn/generic/rshift.c, 32 mpn/generic/sbpi1_div_qr.c, mpn/generic/sub_n.c, 33 mpn/generic/submul_1.c. */ 34 35 #include <assert.h> 36 #include <ctype.h> 37 #include <limits.h> 38 #include <stdio.h> 39 #include <stdlib.h> 40 #include <string.h> 41 42 #include "mini-gmp.h" 43 44 45 /* Macros */ 46 #define GMP_LIMB_BITS (sizeof(mp_limb_t) * CHAR_BIT) 47 48 #define GMP_LIMB_MAX (~ (mp_limb_t) 0) 49 #define GMP_LIMB_HIGHBIT ((mp_limb_t) 1 << (GMP_LIMB_BITS - 1)) 50 51 #define GMP_HLIMB_BIT ((mp_limb_t) 1 << (GMP_LIMB_BITS / 2)) 52 #define GMP_LLIMB_MASK (GMP_HLIMB_BIT - 1) 53 54 #define GMP_ULONG_BITS (sizeof(unsigned long) * CHAR_BIT) 55 #define GMP_ULONG_HIGHBIT ((unsigned long) 1 << (GMP_ULONG_BITS - 1)) 56 57 #define GMP_ABS(x) ((x) >= 0 ? (x) : -(x)) 58 #define GMP_NEG_CAST(T,x) (-((T)((x) + 1) - 1)) 59 60 #define GMP_MIN(a, b) ((a) < (b) ? (a) : (b)) 61 #define GMP_MAX(a, b) ((a) > (b) ? (a) : (b)) 62 63 #define gmp_assert_nocarry(x) do { \ 64 mp_limb_t __cy = x; \ 65 assert (__cy == 0); \ 66 } while (0) 67 68 #define gmp_clz(count, x) do { \ 69 mp_limb_t __clz_x = (x); \ 70 unsigned __clz_c; \ 71 for (__clz_c = 0; \ 72 (__clz_x & ((mp_limb_t) 0xff << (GMP_LIMB_BITS - 8))) == 0; \ 73 __clz_c += 8) \ 74 __clz_x <<= 8; \ 75 for (; (__clz_x & GMP_LIMB_HIGHBIT) == 0; __clz_c++) \ 76 __clz_x <<= 1; \ 77 (count) = __clz_c; \ 78 } while (0) 79 80 #define gmp_ctz(count, x) do { \ 81 mp_limb_t __ctz_x = (x); \ 82 unsigned __ctz_c = 0; \ 83 gmp_clz (__ctz_c, __ctz_x & - __ctz_x); \ 84 (count) = GMP_LIMB_BITS - 1 - __ctz_c; \ 85 } while (0) 86 87 #define gmp_add_ssaaaa(sh, sl, ah, al, bh, bl) \ 88 do { \ 89 mp_limb_t __x; \ 90 __x = (al) + (bl); \ 91 (sh) = (ah) + (bh) + (__x < (al)); \ 92 (sl) = __x; \ 93 } while (0) 94 95 #define gmp_sub_ddmmss(sh, sl, ah, al, bh, bl) \ 96 do { \ 97 mp_limb_t __x; \ 98 __x = (al) - (bl); \ 99 (sh) = (ah) - (bh) - ((al) < (bl)); \ 100 (sl) = __x; \ 101 } while (0) 102 103 #define gmp_umul_ppmm(w1, w0, u, v) \ 104 do { \ 105 mp_limb_t __x0, __x1, __x2, __x3; \ 106 unsigned __ul, __vl, __uh, __vh; \ 107 mp_limb_t __u = (u), __v = (v); \ 108 \ 109 __ul = __u & GMP_LLIMB_MASK; \ 110 __uh = __u >> (GMP_LIMB_BITS / 2); \ 111 __vl = __v & GMP_LLIMB_MASK; \ 112 __vh = __v >> (GMP_LIMB_BITS / 2); \ 113 \ 114 __x0 = (mp_limb_t) __ul * __vl; \ 115 __x1 = (mp_limb_t) __ul * __vh; \ 116 __x2 = (mp_limb_t) __uh * __vl; \ 117 __x3 = (mp_limb_t) __uh * __vh; \ 118 \ 119 __x1 += __x0 >> (GMP_LIMB_BITS / 2);/* this can't give carry */ \ 120 __x1 += __x2; /* but this indeed can */ \ 121 if (__x1 < __x2) /* did we get it? */ \ 122 __x3 += GMP_HLIMB_BIT; /* yes, add it in the proper pos. */ \ 123 \ 124 (w1) = __x3 + (__x1 >> (GMP_LIMB_BITS / 2)); \ 125 (w0) = (__x1 << (GMP_LIMB_BITS / 2)) + (__x0 & GMP_LLIMB_MASK); \ 126 } while (0) 127 128 #define gmp_udiv_qrnnd_preinv(q, r, nh, nl, d, di) \ 129 do { \ 130 mp_limb_t _qh, _ql, _r, _mask; \ 131 gmp_umul_ppmm (_qh, _ql, (nh), (di)); \ 132 gmp_add_ssaaaa (_qh, _ql, _qh, _ql, (nh) + 1, (nl)); \ 133 _r = (nl) - _qh * (d); \ 134 _mask = -(mp_limb_t) (_r > _ql); /* both > and >= are OK */ \ 135 _qh += _mask; \ 136 _r += _mask & (d); \ 137 if (_r >= (d)) \ 138 { \ 139 _r -= (d); \ 140 _qh++; \ 141 } \ 142 \ 143 (r) = _r; \ 144 (q) = _qh; \ 145 } while (0) 146 147 #define gmp_udiv_qr_3by2(q, r1, r0, n2, n1, n0, d1, d0, dinv) \ 148 do { \ 149 mp_limb_t _q0, _t1, _t0, _mask; \ 150 gmp_umul_ppmm ((q), _q0, (n2), (dinv)); \ 151 gmp_add_ssaaaa ((q), _q0, (q), _q0, (n2), (n1)); \ 152 \ 153 /* Compute the two most significant limbs of n - q'd */ \ 154 (r1) = (n1) - (d1) * (q); \ 155 gmp_sub_ddmmss ((r1), (r0), (r1), (n0), (d1), (d0)); \ 156 gmp_umul_ppmm (_t1, _t0, (d0), (q)); \ 157 gmp_sub_ddmmss ((r1), (r0), (r1), (r0), _t1, _t0); \ 158 (q)++; \ 159 \ 160 /* Conditionally adjust q and the remainders */ \ 161 _mask = - (mp_limb_t) ((r1) >= _q0); \ 162 (q) += _mask; \ 163 gmp_add_ssaaaa ((r1), (r0), (r1), (r0), _mask & (d1), _mask & (d0)); \ 164 if ((r1) >= (d1)) \ 165 { \ 166 if ((r1) > (d1) || (r0) >= (d0)) \ 167 { \ 168 (q)++; \ 169 gmp_sub_ddmmss ((r1), (r0), (r1), (r0), (d1), (d0)); \ 170 } \ 171 } \ 172 } while (0) 173 174 /* Swap macros. */ 175 #define MP_LIMB_T_SWAP(x, y) \ 176 do { \ 177 mp_limb_t __mp_limb_t_swap__tmp = (x); \ 178 (x) = (y); \ 179 (y) = __mp_limb_t_swap__tmp; \ 180 } while (0) 181 #define MP_SIZE_T_SWAP(x, y) \ 182 do { \ 183 mp_size_t __mp_size_t_swap__tmp = (x); \ 184 (x) = (y); \ 185 (y) = __mp_size_t_swap__tmp; \ 186 } while (0) 187 #define MP_BITCNT_T_SWAP(x,y) \ 188 do { \ 189 mp_bitcnt_t __mp_bitcnt_t_swap__tmp = (x); \ 190 (x) = (y); \ 191 (y) = __mp_bitcnt_t_swap__tmp; \ 192 } while (0) 193 #define MP_PTR_SWAP(x, y) \ 194 do { \ 195 mp_ptr __mp_ptr_swap__tmp = (x); \ 196 (x) = (y); \ 197 (y) = __mp_ptr_swap__tmp; \ 198 } while (0) 199 #define MP_SRCPTR_SWAP(x, y) \ 200 do { \ 201 mp_srcptr __mp_srcptr_swap__tmp = (x); \ 202 (x) = (y); \ 203 (y) = __mp_srcptr_swap__tmp; \ 204 } while (0) 205 206 #define MPN_PTR_SWAP(xp,xs, yp,ys) \ 207 do { \ 208 MP_PTR_SWAP (xp, yp); \ 209 MP_SIZE_T_SWAP (xs, ys); \ 210 } while(0) 211 #define MPN_SRCPTR_SWAP(xp,xs, yp,ys) \ 212 do { \ 213 MP_SRCPTR_SWAP (xp, yp); \ 214 MP_SIZE_T_SWAP (xs, ys); \ 215 } while(0) 216 217 #define MPZ_PTR_SWAP(x, y) \ 218 do { \ 219 mpz_ptr __mpz_ptr_swap__tmp = (x); \ 220 (x) = (y); \ 221 (y) = __mpz_ptr_swap__tmp; \ 222 } while (0) 223 #define MPZ_SRCPTR_SWAP(x, y) \ 224 do { \ 225 mpz_srcptr __mpz_srcptr_swap__tmp = (x); \ 226 (x) = (y); \ 227 (y) = __mpz_srcptr_swap__tmp; \ 228 } while (0) 229 230 231 /* Memory allocation and other helper functions. */ 232 static void 233 gmp_die (const char *msg) 234 { 235 fprintf (stderr, "%s\n", msg); 236 abort(); 237 } 238 239 static void * 240 gmp_default_alloc (size_t size) 241 { 242 void *p; 243 244 assert (size > 0); 245 246 p = malloc (size); 247 if (!p) 248 gmp_die("gmp_default_alloc: Virtual memory exhausted."); 249 250 return p; 251 } 252 253 static void * 254 gmp_default_realloc (void *old, size_t old_size, size_t new_size) 255 { 256 mp_ptr p; 257 258 p = realloc (old, new_size); 259 260 if (!p) 261 gmp_die("gmp_default_realoc: Virtual memory exhausted."); 262 263 return p; 264 } 265 266 static void 267 gmp_default_free (void *p, size_t size) 268 { 269 free (p); 270 } 271 272 static void * (*gmp_allocate_func) (size_t) = gmp_default_alloc; 273 static void * (*gmp_reallocate_func) (void *, size_t, size_t) = gmp_default_realloc; 274 static void (*gmp_free_func) (void *, size_t) = gmp_default_free; 275 276 void 277 mp_get_memory_functions (void *(**alloc_func) (size_t), 278 void *(**realloc_func) (void *, size_t, size_t), 279 void (**free_func) (void *, size_t)) 280 { 281 if (alloc_func) 282 *alloc_func = gmp_allocate_func; 283 284 if (realloc_func) 285 *realloc_func = gmp_reallocate_func; 286 287 if (free_func) 288 *free_func = gmp_free_func; 289 } 290 291 void 292 mp_set_memory_functions (void *(*alloc_func) (size_t), 293 void *(*realloc_func) (void *, size_t, size_t), 294 void (*free_func) (void *, size_t)) 295 { 296 if (!alloc_func) 297 alloc_func = gmp_default_alloc; 298 if (!realloc_func) 299 realloc_func = gmp_default_realloc; 300 if (!free_func) 301 free_func = gmp_default_free; 302 303 gmp_allocate_func = alloc_func; 304 gmp_reallocate_func = realloc_func; 305 gmp_free_func = free_func; 306 } 307 308 #define gmp_xalloc(size) ((*gmp_allocate_func)((size))) 309 #define gmp_free(p) ((*gmp_free_func) ((p), 0)) 310 311 static mp_ptr 312 gmp_xalloc_limbs (mp_size_t size) 313 { 314 return gmp_xalloc (size * sizeof (mp_limb_t)); 315 } 316 317 static mp_ptr 318 gmp_xrealloc_limbs (mp_ptr old, mp_size_t size) 319 { 320 assert (size > 0); 321 return (*gmp_reallocate_func) (old, 0, size * sizeof (mp_limb_t)); 322 } 323 324 325 /* MPN interface */ 326 327 void 328 mpn_copyi (mp_ptr d, mp_srcptr s, mp_size_t n) 329 { 330 mp_size_t i; 331 for (i = 0; i < n; i++) 332 d[i] = s[i]; 333 } 334 335 void 336 mpn_copyd (mp_ptr d, mp_srcptr s, mp_size_t n) 337 { 338 while (n-- > 0) 339 d[n] = s[n]; 340 } 341 342 int 343 mpn_cmp (mp_srcptr ap, mp_srcptr bp, mp_size_t n) 344 { 345 for (; n > 0; n--) 346 { 347 if (ap[n-1] < bp[n-1]) 348 return -1; 349 else if (ap[n-1] > bp[n-1]) 350 return 1; 351 } 352 return 0; 353 } 354 355 static int 356 mpn_cmp4 (mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn) 357 { 358 if (an > bn) 359 return 1; 360 else if (an < bn) 361 return -1; 362 else 363 return mpn_cmp (ap, bp, an); 364 } 365 366 static mp_size_t 367 mpn_normalized_size (mp_srcptr xp, mp_size_t n) 368 { 369 for (; n > 0 && xp[n-1] == 0; n--) 370 ; 371 return n; 372 } 373 374 #define mpn_zero_p(xp, n) (mpn_normalized_size ((xp), (n)) == 0) 375 376 mp_limb_t 377 mpn_add_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b) 378 { 379 mp_size_t i; 380 381 assert (n > 0); 382 383 for (i = 0; i < n; i++) 384 { 385 mp_limb_t r = ap[i] + b; 386 /* Carry out */ 387 b = (r < b); 388 rp[i] = r; 389 } 390 return b; 391 } 392 393 mp_limb_t 394 mpn_add_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n) 395 { 396 mp_size_t i; 397 mp_limb_t cy; 398 399 for (i = 0, cy = 0; i < n; i++) 400 { 401 mp_limb_t a, b, r; 402 a = ap[i]; b = bp[i]; 403 r = a + cy; 404 cy = (r < cy); 405 r += b; 406 cy += (r < b); 407 rp[i] = r; 408 } 409 return cy; 410 } 411 412 mp_limb_t 413 mpn_add (mp_ptr rp, mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn) 414 { 415 mp_limb_t cy; 416 417 assert (an >= bn); 418 419 cy = mpn_add_n (rp, ap, bp, bn); 420 if (an > bn) 421 cy = mpn_add_1 (rp + bn, ap + bn, an - bn, cy); 422 return cy; 423 } 424 425 mp_limb_t 426 mpn_sub_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b) 427 { 428 mp_size_t i; 429 430 assert (n > 0); 431 432 for (i = 0; i < n; i++) 433 { 434 mp_limb_t a = ap[i]; 435 /* Carry out */ 436 mp_limb_t cy = a < b;; 437 rp[i] = a - b; 438 b = cy; 439 } 440 return b; 441 } 442 443 mp_limb_t 444 mpn_sub_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n) 445 { 446 mp_size_t i; 447 mp_limb_t cy; 448 449 for (i = 0, cy = 0; i < n; i++) 450 { 451 mp_limb_t a, b; 452 a = ap[i]; b = bp[i]; 453 b += cy; 454 cy = (b < cy); 455 cy += (a < b); 456 rp[i] = a - b; 457 } 458 return cy; 459 } 460 461 mp_limb_t 462 mpn_sub (mp_ptr rp, mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn) 463 { 464 mp_limb_t cy; 465 466 assert (an >= bn); 467 468 cy = mpn_sub_n (rp, ap, bp, bn); 469 if (an > bn) 470 cy = mpn_sub_1 (rp + bn, ap + bn, an - bn, cy); 471 return cy; 472 } 473 474 mp_limb_t 475 mpn_mul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl) 476 { 477 mp_limb_t ul, cl, hpl, lpl; 478 479 assert (n >= 1); 480 481 cl = 0; 482 do 483 { 484 ul = *up++; 485 gmp_umul_ppmm (hpl, lpl, ul, vl); 486 487 lpl += cl; 488 cl = (lpl < cl) + hpl; 489 490 *rp++ = lpl; 491 } 492 while (--n != 0); 493 494 return cl; 495 } 496 497 mp_limb_t 498 mpn_addmul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl) 499 { 500 mp_limb_t ul, cl, hpl, lpl, rl; 501 502 assert (n >= 1); 503 504 cl = 0; 505 do 506 { 507 ul = *up++; 508 gmp_umul_ppmm (hpl, lpl, ul, vl); 509 510 lpl += cl; 511 cl = (lpl < cl) + hpl; 512 513 rl = *rp; 514 lpl = rl + lpl; 515 cl += lpl < rl; 516 *rp++ = lpl; 517 } 518 while (--n != 0); 519 520 return cl; 521 } 522 523 mp_limb_t 524 mpn_submul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl) 525 { 526 mp_limb_t ul, cl, hpl, lpl, rl; 527 528 assert (n >= 1); 529 530 cl = 0; 531 do 532 { 533 ul = *up++; 534 gmp_umul_ppmm (hpl, lpl, ul, vl); 535 536 lpl += cl; 537 cl = (lpl < cl) + hpl; 538 539 rl = *rp; 540 lpl = rl - lpl; 541 cl += lpl > rl; 542 *rp++ = lpl; 543 } 544 while (--n != 0); 545 546 return cl; 547 } 548 549 mp_limb_t 550 mpn_mul (mp_ptr rp, mp_srcptr up, mp_size_t un, mp_srcptr vp, mp_size_t vn) 551 { 552 assert (un >= vn); 553 assert (vn >= 1); 554 555 /* We first multiply by the low order limb. This result can be 556 stored, not added, to rp. We also avoid a loop for zeroing this 557 way. */ 558 559 rp[un] = mpn_mul_1 (rp, up, un, vp[0]); 560 rp += 1, vp += 1, vn -= 1; 561 562 /* Now accumulate the product of up[] and the next higher limb from 563 vp[]. */ 564 565 while (vn >= 1) 566 { 567 rp[un] = mpn_addmul_1 (rp, up, un, vp[0]); 568 rp += 1, vp += 1, vn -= 1; 569 } 570 return rp[un - 1]; 571 } 572 573 void 574 mpn_mul_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n) 575 { 576 mpn_mul (rp, ap, n, bp, n); 577 } 578 579 void 580 mpn_sqr (mp_ptr rp, mp_srcptr ap, mp_size_t n) 581 { 582 mpn_mul (rp, ap, n, ap, n); 583 } 584 585 mp_limb_t 586 mpn_lshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt) 587 { 588 mp_limb_t high_limb, low_limb; 589 unsigned int tnc; 590 mp_size_t i; 591 mp_limb_t retval; 592 593 assert (n >= 1); 594 assert (cnt >= 1); 595 assert (cnt < GMP_LIMB_BITS); 596 597 up += n; 598 rp += n; 599 600 tnc = GMP_LIMB_BITS - cnt; 601 low_limb = *--up; 602 retval = low_limb >> tnc; 603 high_limb = (low_limb << cnt); 604 605 for (i = n - 1; i != 0; i--) 606 { 607 low_limb = *--up; 608 *--rp = high_limb | (low_limb >> tnc); 609 high_limb = (low_limb << cnt); 610 } 611 *--rp = high_limb; 612 613 return retval; 614 } 615 616 mp_limb_t 617 mpn_rshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt) 618 { 619 mp_limb_t high_limb, low_limb; 620 unsigned int tnc; 621 mp_size_t i; 622 mp_limb_t retval; 623 624 assert (n >= 1); 625 assert (cnt >= 1); 626 assert (cnt < GMP_LIMB_BITS); 627 628 tnc = GMP_LIMB_BITS - cnt; 629 high_limb = *up++; 630 retval = (high_limb << tnc); 631 low_limb = high_limb >> cnt; 632 633 for (i = n - 1; i != 0; i--) 634 { 635 high_limb = *up++; 636 *rp++ = low_limb | (high_limb << tnc); 637 low_limb = high_limb >> cnt; 638 } 639 *rp = low_limb; 640 641 return retval; 642 } 643 644 645 /* MPN division interface. */ 646 mp_limb_t 647 mpn_invert_3by2 (mp_limb_t u1, mp_limb_t u0) 648 { 649 mp_limb_t r, p, m; 650 unsigned ul, uh; 651 unsigned ql, qh; 652 653 /* First, do a 2/1 inverse. */ 654 /* The inverse m is defined as floor( (B^2 - 1 - u1)/u1 ), so that 0 < 655 * B^2 - (B + m) u1 <= u1 */ 656 assert (u1 >= GMP_LIMB_HIGHBIT); 657 658 ul = u1 & GMP_LLIMB_MASK; 659 uh = u1 >> (GMP_LIMB_BITS / 2); 660 661 qh = ~u1 / uh; 662 r = ((~u1 - (mp_limb_t) qh * uh) << (GMP_LIMB_BITS / 2)) | GMP_LLIMB_MASK; 663 664 p = (mp_limb_t) qh * ul; 665 /* Adjustment steps taken from udiv_qrnnd_c */ 666 if (r < p) 667 { 668 qh--; 669 r += u1; 670 if (r >= u1) /* i.e. we didn't get carry when adding to r */ 671 if (r < p) 672 { 673 qh--; 674 r += u1; 675 } 676 } 677 r -= p; 678 679 /* Do a 3/2 division (with half limb size) */ 680 p = (r >> (GMP_LIMB_BITS / 2)) * qh + r; 681 ql = (p >> (GMP_LIMB_BITS / 2)) + 1; 682 683 /* By the 3/2 method, we don't need the high half limb. */ 684 r = (r << (GMP_LIMB_BITS / 2)) + GMP_LLIMB_MASK - ql * u1; 685 686 if (r >= (p << (GMP_LIMB_BITS / 2))) 687 { 688 ql--; 689 r += u1; 690 } 691 m = ((mp_limb_t) qh << (GMP_LIMB_BITS / 2)) + ql; 692 if (r >= u1) 693 { 694 m++; 695 r -= u1; 696 } 697 698 if (u0 > 0) 699 { 700 mp_limb_t th, tl; 701 r = ~r; 702 r += u0; 703 if (r < u0) 704 { 705 m--; 706 if (r >= u1) 707 { 708 m--; 709 r -= u1; 710 } 711 r -= u1; 712 } 713 gmp_umul_ppmm (th, tl, u0, m); 714 r += th; 715 if (r < th) 716 { 717 m--; 718 if (r > u1 || (r == u1 && tl > u0)) 719 m--; 720 } 721 } 722 723 return m; 724 } 725 726 struct gmp_div_inverse 727 { 728 /* Normalization shift count. */ 729 unsigned shift; 730 /* Normalized divisor (d0 unused for mpn_div_qr_1) */ 731 mp_limb_t d1, d0; 732 /* Inverse, for 2/1 or 3/2. */ 733 mp_limb_t di; 734 }; 735 736 static void 737 mpn_div_qr_1_invert (struct gmp_div_inverse *inv, mp_limb_t d) 738 { 739 unsigned shift; 740 741 assert (d > 0); 742 gmp_clz (shift, d); 743 inv->shift = shift; 744 inv->d1 = d << shift; 745 inv->di = mpn_invert_limb (inv->d1); 746 } 747 748 static void 749 mpn_div_qr_2_invert (struct gmp_div_inverse *inv, 750 mp_limb_t d1, mp_limb_t d0) 751 { 752 unsigned shift; 753 754 assert (d1 > 0); 755 gmp_clz (shift, d1); 756 inv->shift = shift; 757 if (shift > 0) 758 { 759 d1 = (d1 << shift) | (d0 >> (GMP_LIMB_BITS - shift)); 760 d0 <<= shift; 761 } 762 inv->d1 = d1; 763 inv->d0 = d0; 764 inv->di = mpn_invert_3by2 (d1, d0); 765 } 766 767 static void 768 mpn_div_qr_invert (struct gmp_div_inverse *inv, 769 mp_srcptr dp, mp_size_t dn) 770 { 771 assert (dn > 0); 772 773 if (dn == 1) 774 mpn_div_qr_1_invert (inv, dp[0]); 775 else if (dn == 2) 776 mpn_div_qr_2_invert (inv, dp[1], dp[0]); 777 else 778 { 779 unsigned shift; 780 mp_limb_t d1, d0; 781 782 d1 = dp[dn-1]; 783 d0 = dp[dn-2]; 784 assert (d1 > 0); 785 gmp_clz (shift, d1); 786 inv->shift = shift; 787 if (shift > 0) 788 { 789 d1 = (d1 << shift) | (d0 >> (GMP_LIMB_BITS - shift)); 790 d0 = (d0 << shift) | (dp[dn-3] >> (GMP_LIMB_BITS - shift)); 791 } 792 inv->d1 = d1; 793 inv->d0 = d0; 794 inv->di = mpn_invert_3by2 (d1, d0); 795 } 796 } 797 798 /* Not matching current public gmp interface, rather corresponding to 799 the sbpi1_div_* functions. */ 800 static mp_limb_t 801 mpn_div_qr_1_preinv (mp_ptr qp, mp_srcptr np, mp_size_t nn, 802 const struct gmp_div_inverse *inv) 803 { 804 mp_limb_t d, di; 805 mp_limb_t r; 806 mp_ptr tp = NULL; 807 808 if (inv->shift > 0) 809 { 810 tp = gmp_xalloc_limbs (nn); 811 r = mpn_lshift (tp, np, nn, inv->shift); 812 np = tp; 813 } 814 else 815 r = 0; 816 817 d = inv->d1; 818 di = inv->di; 819 while (nn-- > 0) 820 { 821 mp_limb_t q; 822 823 gmp_udiv_qrnnd_preinv (q, r, r, np[nn], d, di); 824 if (qp) 825 qp[nn] = q; 826 } 827 if (inv->shift > 0) 828 gmp_free (tp); 829 830 return r >> inv->shift; 831 } 832 833 static mp_limb_t 834 mpn_div_qr_1 (mp_ptr qp, mp_srcptr np, mp_size_t nn, mp_limb_t d) 835 { 836 assert (d > 0); 837 838 /* Special case for powers of two. */ 839 if (d > 1 && (d & (d-1)) == 0) 840 { 841 unsigned shift; 842 mp_limb_t r = np[0] & (d-1); 843 gmp_ctz (shift, d); 844 if (qp) 845 mpn_rshift (qp, np, nn, shift); 846 847 return r; 848 } 849 else 850 { 851 struct gmp_div_inverse inv; 852 mpn_div_qr_1_invert (&inv, d); 853 return mpn_div_qr_1_preinv (qp, np, nn, &inv); 854 } 855 } 856 857 static void 858 mpn_div_qr_2_preinv (mp_ptr qp, mp_ptr rp, mp_srcptr np, mp_size_t nn, 859 const struct gmp_div_inverse *inv) 860 { 861 unsigned shift; 862 mp_size_t i; 863 mp_limb_t d1, d0, di, r1, r0; 864 mp_ptr tp; 865 866 assert (nn >= 2); 867 shift = inv->shift; 868 d1 = inv->d1; 869 d0 = inv->d0; 870 di = inv->di; 871 872 if (shift > 0) 873 { 874 tp = gmp_xalloc_limbs (nn); 875 r1 = mpn_lshift (tp, np, nn, shift); 876 np = tp; 877 } 878 else 879 r1 = 0; 880 881 r0 = np[nn - 1]; 882 883 for (i = nn - 2; i >= 0; i--) 884 { 885 mp_limb_t n0, q; 886 n0 = np[i]; 887 gmp_udiv_qr_3by2 (q, r1, r0, r1, r0, n0, d1, d0, di); 888 889 if (qp) 890 qp[i] = q; 891 } 892 893 if (shift > 0) 894 { 895 assert ((r0 << (GMP_LIMB_BITS - shift)) == 0); 896 r0 = (r0 >> shift) | (r1 << (GMP_LIMB_BITS - shift)); 897 r1 >>= shift; 898 899 gmp_free (tp); 900 } 901 902 rp[1] = r1; 903 rp[0] = r0; 904 } 905 906 #if 0 907 static void 908 mpn_div_qr_2 (mp_ptr qp, mp_ptr rp, mp_srcptr np, mp_size_t nn, 909 mp_limb_t d1, mp_limb_t d0) 910 { 911 struct gmp_div_inverse inv; 912 assert (nn >= 2); 913 914 mpn_div_qr_2_invert (&inv, d1, d0); 915 mpn_div_qr_2_preinv (qp, rp, np, nn, &inv); 916 } 917 #endif 918 919 static void 920 mpn_div_qr_pi1 (mp_ptr qp, 921 mp_ptr np, mp_size_t nn, mp_limb_t n1, 922 mp_srcptr dp, mp_size_t dn, 923 mp_limb_t dinv) 924 { 925 mp_size_t i; 926 927 mp_limb_t d1, d0; 928 mp_limb_t cy, cy1; 929 mp_limb_t q; 930 931 assert (dn > 2); 932 assert (nn >= dn); 933 934 d1 = dp[dn - 1]; 935 d0 = dp[dn - 2]; 936 937 assert ((d1 & GMP_LIMB_HIGHBIT) != 0); 938 /* Iteration variable is the index of the q limb. 939 * 940 * We divide <n1, np[dn-1+i], np[dn-2+i], np[dn-3+i],..., np[i]> 941 * by <d1, d0, dp[dn-3], ..., dp[0] > 942 */ 943 944 for (i = nn - dn; i >= 0; i--) 945 { 946 mp_limb_t n0 = np[dn-1+i]; 947 948 if (n1 == d1 && n0 == d0) 949 { 950 q = GMP_LIMB_MAX; 951 mpn_submul_1 (np+i, dp, dn, q); 952 n1 = np[dn-1+i]; /* update n1, last loop's value will now be invalid */ 953 } 954 else 955 { 956 gmp_udiv_qr_3by2 (q, n1, n0, n1, n0, np[dn-2+i], d1, d0, dinv); 957 958 cy = mpn_submul_1 (np + i, dp, dn-2, q); 959 960 cy1 = n0 < cy; 961 n0 = n0 - cy; 962 cy = n1 < cy1; 963 n1 = n1 - cy1; 964 np[dn-2+i] = n0; 965 966 if (cy != 0) 967 { 968 n1 += d1 + mpn_add_n (np + i, np + i, dp, dn - 1); 969 q--; 970 } 971 } 972 973 if (qp) 974 qp[i] = q; 975 } 976 977 np[dn - 1] = n1; 978 } 979 980 static void 981 mpn_div_qr_preinv (mp_ptr qp, mp_ptr np, mp_size_t nn, 982 mp_srcptr dp, mp_size_t dn, 983 const struct gmp_div_inverse *inv) 984 { 985 assert (dn > 0); 986 assert (nn >= dn); 987 988 if (dn == 1) 989 np[0] = mpn_div_qr_1_preinv (qp, np, nn, inv); 990 else if (dn == 2) 991 mpn_div_qr_2_preinv (qp, np, np, nn, inv); 992 else 993 { 994 mp_limb_t nh; 995 unsigned shift; 996 997 assert (inv->d1 == dp[dn-1]); 998 assert (inv->d0 == dp[dn-2]); 999 assert ((inv->d1 & GMP_LIMB_HIGHBIT) != 0); 1000 1001 shift = inv->shift; 1002 if (shift > 0) 1003 nh = mpn_lshift (np, np, nn, shift); 1004 else 1005 nh = 0; 1006 1007 mpn_div_qr_pi1 (qp, np, nn, nh, dp, dn, inv->di); 1008 1009 if (shift > 0) 1010 gmp_assert_nocarry (mpn_rshift (np, np, dn, shift)); 1011 } 1012 } 1013 1014 static void 1015 mpn_div_qr (mp_ptr qp, mp_ptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn) 1016 { 1017 struct gmp_div_inverse inv; 1018 mp_ptr tp = NULL; 1019 1020 assert (dn > 0); 1021 assert (nn >= dn); 1022 1023 mpn_div_qr_invert (&inv, dp, dn); 1024 if (dn > 2 && inv.shift > 0) 1025 { 1026 tp = gmp_xalloc_limbs (dn); 1027 gmp_assert_nocarry (mpn_lshift (tp, dp, dn, inv.shift)); 1028 dp = tp; 1029 } 1030 mpn_div_qr_preinv (qp, np, nn, dp, dn, &inv); 1031 if (tp) 1032 gmp_free (tp); 1033 } 1034 1035 1036 /* MPN base conversion. */ 1037 static unsigned 1038 mpn_base_power_of_two_p (unsigned b) 1039 { 1040 switch (b) 1041 { 1042 case 2: return 1; 1043 case 4: return 2; 1044 case 8: return 3; 1045 case 16: return 4; 1046 case 32: return 5; 1047 case 64: return 6; 1048 case 128: return 7; 1049 case 256: return 8; 1050 default: return 0; 1051 } 1052 } 1053 1054 struct mpn_base_info 1055 { 1056 /* bb is the largest power of the base which fits in one limb, and 1057 exp is the corresponding exponent. */ 1058 unsigned exp; 1059 mp_limb_t bb; 1060 }; 1061 1062 static void 1063 mpn_get_base_info (struct mpn_base_info *info, mp_limb_t b) 1064 { 1065 mp_limb_t m; 1066 mp_limb_t p; 1067 unsigned exp; 1068 1069 m = GMP_LIMB_MAX / b; 1070 for (exp = 1, p = b; p <= m; exp++) 1071 p *= b; 1072 1073 info->exp = exp; 1074 info->bb = p; 1075 } 1076 1077 static mp_bitcnt_t 1078 mpn_limb_size_in_base_2 (mp_limb_t u) 1079 { 1080 unsigned shift; 1081 1082 assert (u > 0); 1083 gmp_clz (shift, u); 1084 return GMP_LIMB_BITS - shift; 1085 } 1086 1087 static size_t 1088 mpn_get_str_bits (unsigned char *sp, unsigned bits, mp_srcptr up, mp_size_t un) 1089 { 1090 unsigned char mask; 1091 size_t sn, j; 1092 mp_size_t i; 1093 int shift; 1094 1095 sn = ((un - 1) * GMP_LIMB_BITS + mpn_limb_size_in_base_2 (up[un-1]) 1096 + bits - 1) / bits; 1097 1098 mask = (1U << bits) - 1; 1099 1100 for (i = 0, j = sn, shift = 0; j-- > 0;) 1101 { 1102 unsigned char digit = up[i] >> shift; 1103 1104 shift += bits; 1105 1106 if (shift >= GMP_LIMB_BITS && ++i < un) 1107 { 1108 shift -= GMP_LIMB_BITS; 1109 digit |= up[i] << (bits - shift); 1110 } 1111 sp[j] = digit & mask; 1112 } 1113 return sn; 1114 } 1115 1116 /* We generate digits from the least significant end, and reverse at 1117 the end. */ 1118 static size_t 1119 mpn_limb_get_str (unsigned char *sp, mp_limb_t w, 1120 const struct gmp_div_inverse *binv) 1121 { 1122 mp_size_t i; 1123 for (i = 0; w > 0; i++) 1124 { 1125 mp_limb_t h, l, r; 1126 1127 h = w >> (GMP_LIMB_BITS - binv->shift); 1128 l = w << binv->shift; 1129 1130 gmp_udiv_qrnnd_preinv (w, r, h, l, binv->d1, binv->di); 1131 assert ( (r << (GMP_LIMB_BITS - binv->shift)) == 0); 1132 r >>= binv->shift; 1133 1134 sp[i] = r; 1135 } 1136 return i; 1137 } 1138 1139 static size_t 1140 mpn_get_str_other (unsigned char *sp, 1141 int base, const struct mpn_base_info *info, 1142 mp_ptr up, mp_size_t un) 1143 { 1144 struct gmp_div_inverse binv; 1145 size_t sn; 1146 size_t i; 1147 1148 mpn_div_qr_1_invert (&binv, base); 1149 1150 sn = 0; 1151 1152 if (un > 1) 1153 { 1154 struct gmp_div_inverse bbinv; 1155 mpn_div_qr_1_invert (&bbinv, info->bb); 1156 1157 do 1158 { 1159 mp_limb_t w; 1160 size_t done; 1161 w = mpn_div_qr_1_preinv (up, up, un, &bbinv); 1162 un -= (up[un-1] == 0); 1163 done = mpn_limb_get_str (sp + sn, w, &binv); 1164 1165 for (sn += done; done < info->exp; done++) 1166 sp[sn++] = 0; 1167 } 1168 while (un > 1); 1169 } 1170 sn += mpn_limb_get_str (sp + sn, up[0], &binv); 1171 1172 /* Reverse order */ 1173 for (i = 0; 2*i + 1 < sn; i++) 1174 { 1175 unsigned char t = sp[i]; 1176 sp[i] = sp[sn - i - 1]; 1177 sp[sn - i - 1] = t; 1178 } 1179 1180 return sn; 1181 } 1182 1183 size_t 1184 mpn_get_str (unsigned char *sp, int base, mp_ptr up, mp_size_t un) 1185 { 1186 unsigned bits; 1187 1188 assert (un > 0); 1189 assert (up[un-1] > 0); 1190 1191 bits = mpn_base_power_of_two_p (base); 1192 if (bits) 1193 return mpn_get_str_bits (sp, bits, up, un); 1194 else 1195 { 1196 struct mpn_base_info info; 1197 1198 mpn_get_base_info (&info, base); 1199 return mpn_get_str_other (sp, base, &info, up, un); 1200 } 1201 } 1202 1203 static mp_size_t 1204 mpn_set_str_bits (mp_ptr rp, const unsigned char *sp, size_t sn, 1205 unsigned bits) 1206 { 1207 mp_size_t rn; 1208 size_t j; 1209 unsigned shift; 1210 1211 for (j = sn, rn = 0, shift = 0; j-- > 0; ) 1212 { 1213 if (shift == 0) 1214 { 1215 rp[rn++] = sp[j]; 1216 shift += bits; 1217 } 1218 else 1219 { 1220 rp[rn-1] |= (mp_limb_t) sp[j] << shift; 1221 shift += bits; 1222 if (shift >= GMP_LIMB_BITS) 1223 { 1224 shift -= GMP_LIMB_BITS; 1225 if (shift > 0) 1226 rp[rn++] = (mp_limb_t) sp[j] >> (bits - shift); 1227 } 1228 } 1229 } 1230 rn = mpn_normalized_size (rp, rn); 1231 return rn; 1232 } 1233 1234 static mp_size_t 1235 mpn_set_str_other (mp_ptr rp, const unsigned char *sp, size_t sn, 1236 mp_limb_t b, const struct mpn_base_info *info) 1237 { 1238 mp_size_t rn; 1239 mp_limb_t w; 1240 unsigned first; 1241 unsigned k; 1242 size_t j; 1243 1244 first = 1 + (sn - 1) % info->exp; 1245 1246 j = 0; 1247 w = sp[j++]; 1248 for (k = 1; k < first; k++) 1249 w = w * b + sp[j++]; 1250 1251 rp[0] = w; 1252 1253 for (rn = (w > 0); j < sn;) 1254 { 1255 mp_limb_t cy; 1256 1257 w = sp[j++]; 1258 for (k = 1; k < info->exp; k++) 1259 w = w * b + sp[j++]; 1260 1261 cy = mpn_mul_1 (rp, rp, rn, info->bb); 1262 cy += mpn_add_1 (rp, rp, rn, w); 1263 if (cy > 0) 1264 rp[rn++] = cy; 1265 } 1266 assert (j == sn); 1267 1268 return rn; 1269 } 1270 1271 mp_size_t 1272 mpn_set_str (mp_ptr rp, const unsigned char *sp, size_t sn, int base) 1273 { 1274 unsigned bits; 1275 1276 if (sn == 0) 1277 return 0; 1278 1279 bits = mpn_base_power_of_two_p (base); 1280 if (bits) 1281 return mpn_set_str_bits (rp, sp, sn, bits); 1282 else 1283 { 1284 struct mpn_base_info info; 1285 1286 mpn_get_base_info (&info, base); 1287 return mpn_set_str_other (rp, sp, sn, base, &info); 1288 } 1289 } 1290 1291 1292 /* MPZ interface */ 1293 void 1294 mpz_init (mpz_t r) 1295 { 1296 r->_mp_alloc = 1; 1297 r->_mp_size = 0; 1298 r->_mp_d = gmp_xalloc_limbs (1); 1299 } 1300 1301 /* The utility of this function is a bit limited, since many functions 1302 assings the result variable using mpz_swap. */ 1303 void 1304 mpz_init2 (mpz_t r, mp_bitcnt_t bits) 1305 { 1306 mp_size_t rn; 1307 1308 bits -= (bits != 0); /* Round down, except if 0 */ 1309 rn = 1 + bits / GMP_LIMB_BITS; 1310 1311 r->_mp_alloc = rn; 1312 r->_mp_size = 0; 1313 r->_mp_d = gmp_xalloc_limbs (rn); 1314 } 1315 1316 void 1317 mpz_clear (mpz_t r) 1318 { 1319 gmp_free (r->_mp_d); 1320 } 1321 1322 static void * 1323 mpz_realloc (mpz_t r, mp_size_t size) 1324 { 1325 size = GMP_MAX (size, 1); 1326 1327 r->_mp_d = gmp_xrealloc_limbs (r->_mp_d, size); 1328 r->_mp_alloc = size; 1329 1330 if (GMP_ABS (r->_mp_size) > size) 1331 r->_mp_size = 0; 1332 1333 return r->_mp_d; 1334 } 1335 1336 /* Realloc for an mpz_t WHAT if it has less than NEEDED limbs. */ 1337 #define MPZ_REALLOC(z,n) ((n) > (z)->_mp_alloc \ 1338 ? mpz_realloc(z,n) \ 1339 : (z)->_mp_d) 1340 1341 /* MPZ assignment and basic conversions. */ 1342 void 1343 mpz_set_si (mpz_t r, signed long int x) 1344 { 1345 if (x >= 0) 1346 mpz_set_ui (r, x); 1347 else /* (x < 0) */ 1348 { 1349 r->_mp_size = -1; 1350 r->_mp_d[0] = GMP_NEG_CAST (unsigned long int, x); 1351 } 1352 } 1353 1354 void 1355 mpz_set_ui (mpz_t r, unsigned long int x) 1356 { 1357 if (x > 0) 1358 { 1359 r->_mp_size = 1; 1360 r->_mp_d[0] = x; 1361 } 1362 else 1363 r->_mp_size = 0; 1364 } 1365 1366 void 1367 mpz_set (mpz_t r, const mpz_t x) 1368 { 1369 /* Allow the NOP r == x */ 1370 if (r != x) 1371 { 1372 mp_size_t n; 1373 mp_ptr rp; 1374 1375 n = GMP_ABS (x->_mp_size); 1376 rp = MPZ_REALLOC (r, n); 1377 1378 mpn_copyi (rp, x->_mp_d, n); 1379 r->_mp_size = x->_mp_size; 1380 } 1381 } 1382 1383 void 1384 mpz_init_set_si (mpz_t r, signed long int x) 1385 { 1386 mpz_init (r); 1387 mpz_set_si (r, x); 1388 } 1389 1390 void 1391 mpz_init_set_ui (mpz_t r, unsigned long int x) 1392 { 1393 mpz_init (r); 1394 mpz_set_ui (r, x); 1395 } 1396 1397 void 1398 mpz_init_set (mpz_t r, const mpz_t x) 1399 { 1400 mpz_init (r); 1401 mpz_set (r, x); 1402 } 1403 1404 int 1405 mpz_fits_slong_p (const mpz_t u) 1406 { 1407 mp_size_t us = u->_mp_size; 1408 1409 if (us == 0) 1410 return 1; 1411 else if (us == 1) 1412 return u->_mp_d[0] < GMP_LIMB_HIGHBIT; 1413 else if (us == -1) 1414 return u->_mp_d[0] <= GMP_LIMB_HIGHBIT; 1415 else 1416 return 0; 1417 } 1418 1419 int 1420 mpz_fits_ulong_p (const mpz_t u) 1421 { 1422 mp_size_t us = u->_mp_size; 1423 1424 return us == 0 || us == 1; 1425 } 1426 1427 long int 1428 mpz_get_si (const mpz_t u) 1429 { 1430 mp_size_t us = u->_mp_size; 1431 1432 if (us > 0) 1433 return (long) (u->_mp_d[0] & ~GMP_LIMB_HIGHBIT); 1434 else if (us < 0) 1435 return (long) (- u->_mp_d[0] | GMP_LIMB_HIGHBIT); 1436 else 1437 return 0; 1438 } 1439 1440 unsigned long int 1441 mpz_get_ui (const mpz_t u) 1442 { 1443 return u->_mp_size == 0 ? 0 : u->_mp_d[0]; 1444 } 1445 1446 size_t 1447 mpz_size (const mpz_t u) 1448 { 1449 return GMP_ABS (u->_mp_size); 1450 } 1451 1452 mp_limb_t 1453 mpz_getlimbn (const mpz_t u, mp_size_t n) 1454 { 1455 if (n >= 0 && n < GMP_ABS (u->_mp_size)) 1456 return u->_mp_d[n]; 1457 else 1458 return 0; 1459 } 1460 1461 1462 /* Conversions and comparison to double. */ 1463 void 1464 mpz_set_d (mpz_t r, double x) 1465 { 1466 int sign; 1467 mp_ptr rp; 1468 mp_size_t rn, i; 1469 double B; 1470 double Bi; 1471 mp_limb_t f; 1472 1473 /* x != x is true when x is a NaN, and x == x * 0.5 is true when x is 1474 zero or infinity. */ 1475 if (x == 0.0 || x != x || x == x * 0.5) 1476 { 1477 r->_mp_size = 0; 1478 return; 1479 } 1480 1481 if (x < 0.0) 1482 { 1483 x = - x; 1484 sign = 1; 1485 } 1486 else 1487 sign = 0; 1488 1489 if (x < 1.0) 1490 { 1491 r->_mp_size = 0; 1492 return; 1493 } 1494 B = 2.0 * (double) GMP_LIMB_HIGHBIT; 1495 Bi = 1.0 / B; 1496 for (rn = 1; x >= B; rn++) 1497 x *= Bi; 1498 1499 rp = MPZ_REALLOC (r, rn); 1500 1501 f = (mp_limb_t) x; 1502 x -= f; 1503 assert (x < 1.0); 1504 rp[rn-1] = f; 1505 for (i = rn-1; i-- > 0; ) 1506 { 1507 x = B * x; 1508 f = (mp_limb_t) x; 1509 x -= f; 1510 assert (x < 1.0); 1511 rp[i] = f; 1512 } 1513 1514 r->_mp_size = sign ? - rn : rn; 1515 } 1516 1517 void 1518 mpz_init_set_d (mpz_t r, double x) 1519 { 1520 mpz_init (r); 1521 mpz_set_d (r, x); 1522 } 1523 1524 double 1525 mpz_get_d (const mpz_t u) 1526 { 1527 mp_size_t un; 1528 double x; 1529 double B = 2.0 * (double) GMP_LIMB_HIGHBIT; 1530 1531 un = GMP_ABS (u->_mp_size); 1532 1533 if (un == 0) 1534 return 0.0; 1535 1536 x = u->_mp_d[--un]; 1537 while (un > 0) 1538 x = B*x + u->_mp_d[--un]; 1539 1540 if (u->_mp_size < 0) 1541 x = -x; 1542 1543 return x; 1544 } 1545 1546 int 1547 mpz_cmpabs_d (const mpz_t x, double d) 1548 { 1549 mp_size_t xn; 1550 double B, Bi; 1551 mp_size_t i; 1552 1553 xn = x->_mp_size; 1554 d = GMP_ABS (d); 1555 1556 if (xn != 0) 1557 { 1558 xn = GMP_ABS (xn); 1559 1560 B = 2.0 * (double) GMP_LIMB_HIGHBIT; 1561 Bi = 1.0 / B; 1562 1563 /* Scale d so it can be compared with the top limb. */ 1564 for (i = 1; i < xn; i++) 1565 d *= Bi; 1566 1567 if (d >= B) 1568 return -1; 1569 1570 /* Compare floor(d) to top limb, subtract and cancel when equal. */ 1571 for (i = xn; i-- > 0;) 1572 { 1573 mp_limb_t f, xl; 1574 1575 f = (mp_limb_t) d; 1576 xl = x->_mp_d[i]; 1577 if (xl > f) 1578 return 1; 1579 else if (xl < f) 1580 return -1; 1581 d = B * (d - f); 1582 } 1583 } 1584 return - (d > 0.0); 1585 } 1586 1587 int 1588 mpz_cmp_d (const mpz_t x, double d) 1589 { 1590 if (x->_mp_size < 0) 1591 { 1592 if (d >= 0.0) 1593 return -1; 1594 else 1595 return -mpz_cmpabs_d (x, d); 1596 } 1597 else 1598 { 1599 if (d < 0.0) 1600 return 1; 1601 else 1602 return mpz_cmpabs_d (x, d); 1603 } 1604 } 1605 1606 1607 /* MPZ comparisons and the like. */ 1608 int 1609 mpz_sgn (const mpz_t u) 1610 { 1611 mp_size_t usize = u->_mp_size; 1612 1613 if (usize > 0) 1614 return 1; 1615 else if (usize < 0) 1616 return -1; 1617 else 1618 return 0; 1619 } 1620 1621 int 1622 mpz_cmp_si (const mpz_t u, long v) 1623 { 1624 mp_size_t usize = u->_mp_size; 1625 1626 if (usize < -1) 1627 return -1; 1628 else if (v >= 0) 1629 return mpz_cmp_ui (u, v); 1630 else if (usize >= 0) 1631 return 1; 1632 else /* usize == -1 */ 1633 { 1634 mp_limb_t ul = u->_mp_d[0]; 1635 if ((mp_limb_t)GMP_NEG_CAST (unsigned long int, v) < ul) 1636 return -1; 1637 else if ( (mp_limb_t)GMP_NEG_CAST (unsigned long int, v) > ul) 1638 return 1; 1639 } 1640 return 0; 1641 } 1642 1643 int 1644 mpz_cmp_ui (const mpz_t u, unsigned long v) 1645 { 1646 mp_size_t usize = u->_mp_size; 1647 1648 if (usize > 1) 1649 return 1; 1650 else if (usize < 0) 1651 return -1; 1652 else 1653 { 1654 mp_limb_t ul = (usize > 0) ? u->_mp_d[0] : 0; 1655 if (ul > v) 1656 return 1; 1657 else if (ul < v) 1658 return -1; 1659 } 1660 return 0; 1661 } 1662 1663 int 1664 mpz_cmp (const mpz_t a, const mpz_t b) 1665 { 1666 mp_size_t asize = a->_mp_size; 1667 mp_size_t bsize = b->_mp_size; 1668 1669 if (asize > bsize) 1670 return 1; 1671 else if (asize < bsize) 1672 return -1; 1673 else if (asize > 0) 1674 return mpn_cmp (a->_mp_d, b->_mp_d, asize); 1675 else if (asize < 0) 1676 return -mpn_cmp (a->_mp_d, b->_mp_d, -asize); 1677 else 1678 return 0; 1679 } 1680 1681 int 1682 mpz_cmpabs_ui (const mpz_t u, unsigned long v) 1683 { 1684 mp_size_t un = GMP_ABS (u->_mp_size); 1685 mp_limb_t ul; 1686 1687 if (un > 1) 1688 return 1; 1689 1690 ul = (un == 1) ? u->_mp_d[0] : 0; 1691 1692 if (ul > v) 1693 return 1; 1694 else if (ul < v) 1695 return -1; 1696 else 1697 return 0; 1698 } 1699 1700 int 1701 mpz_cmpabs (const mpz_t u, const mpz_t v) 1702 { 1703 return mpn_cmp4 (u->_mp_d, GMP_ABS (u->_mp_size), 1704 v->_mp_d, GMP_ABS (v->_mp_size)); 1705 } 1706 1707 void 1708 mpz_abs (mpz_t r, const mpz_t u) 1709 { 1710 if (r != u) 1711 mpz_set (r, u); 1712 1713 r->_mp_size = GMP_ABS (r->_mp_size); 1714 } 1715 1716 void 1717 mpz_neg (mpz_t r, const mpz_t u) 1718 { 1719 if (r != u) 1720 mpz_set (r, u); 1721 1722 r->_mp_size = -r->_mp_size; 1723 } 1724 1725 void 1726 mpz_swap (mpz_t u, mpz_t v) 1727 { 1728 MP_SIZE_T_SWAP (u->_mp_size, v->_mp_size); 1729 MP_SIZE_T_SWAP (u->_mp_alloc, v->_mp_alloc); 1730 MP_PTR_SWAP (u->_mp_d, v->_mp_d); 1731 } 1732 1733 1734 /* MPZ addition and subtraction */ 1735 1736 /* Adds to the absolute value. Returns new size, but doesn't store it. */ 1737 static mp_size_t 1738 mpz_abs_add_ui (mpz_t r, const mpz_t a, unsigned long b) 1739 { 1740 mp_size_t an; 1741 mp_ptr rp; 1742 mp_limb_t cy; 1743 1744 an = GMP_ABS (a->_mp_size); 1745 if (an == 0) 1746 { 1747 r->_mp_d[0] = b; 1748 return b > 0; 1749 } 1750 1751 rp = MPZ_REALLOC (r, an + 1); 1752 1753 cy = mpn_add_1 (rp, a->_mp_d, an, b); 1754 rp[an] = cy; 1755 an += (cy > 0); 1756 1757 return an; 1758 } 1759 1760 /* Subtract from the absolute value. Returns new size, (or -1 on underflow), 1761 but doesn't store it. */ 1762 static mp_size_t 1763 mpz_abs_sub_ui (mpz_t r, const mpz_t a, unsigned long b) 1764 { 1765 mp_size_t an = GMP_ABS (a->_mp_size); 1766 mp_ptr rp = MPZ_REALLOC (r, an); 1767 1768 if (an == 0) 1769 { 1770 rp[0] = b; 1771 return -(b > 0); 1772 } 1773 else if (an == 1 && a->_mp_d[0] < b) 1774 { 1775 rp[0] = b - a->_mp_d[0]; 1776 return -1; 1777 } 1778 else 1779 { 1780 gmp_assert_nocarry (mpn_sub_1 (rp, a->_mp_d, an, b)); 1781 return mpn_normalized_size (rp, an); 1782 } 1783 } 1784 1785 void 1786 mpz_add_ui (mpz_t r, const mpz_t a, unsigned long b) 1787 { 1788 if (a->_mp_size >= 0) 1789 r->_mp_size = mpz_abs_add_ui (r, a, b); 1790 else 1791 r->_mp_size = -mpz_abs_sub_ui (r, a, b); 1792 } 1793 1794 void 1795 mpz_sub_ui (mpz_t r, const mpz_t a, unsigned long b) 1796 { 1797 if (a->_mp_size < 0) 1798 r->_mp_size = -mpz_abs_add_ui (r, a, b); 1799 else 1800 r->_mp_size = mpz_abs_sub_ui (r, a, b); 1801 } 1802 1803 void 1804 mpz_ui_sub (mpz_t r, unsigned long a, const mpz_t b) 1805 { 1806 if (b->_mp_size < 0) 1807 r->_mp_size = mpz_abs_add_ui (r, b, a); 1808 else 1809 r->_mp_size = -mpz_abs_sub_ui (r, b, a); 1810 } 1811 1812 static mp_size_t 1813 mpz_abs_add (mpz_t r, const mpz_t a, const mpz_t b) 1814 { 1815 mp_size_t an = GMP_ABS (a->_mp_size); 1816 mp_size_t bn = GMP_ABS (b->_mp_size); 1817 mp_size_t rn; 1818 mp_ptr rp; 1819 mp_limb_t cy; 1820 1821 rn = GMP_MAX (an, bn); 1822 rp = MPZ_REALLOC (r, rn + 1); 1823 if (an >= bn) 1824 cy = mpn_add (rp, a->_mp_d, an, b->_mp_d, bn); 1825 else 1826 cy = mpn_add (rp, b->_mp_d, bn, a->_mp_d, an); 1827 1828 rp[rn] = cy; 1829 1830 return rn + (cy > 0); 1831 } 1832 1833 static mp_size_t 1834 mpz_abs_sub (mpz_t r, const mpz_t a, const mpz_t b) 1835 { 1836 mp_size_t an = GMP_ABS (a->_mp_size); 1837 mp_size_t bn = GMP_ABS (b->_mp_size); 1838 int cmp; 1839 mp_ptr rp; 1840 1841 cmp = mpn_cmp4 (a->_mp_d, an, b->_mp_d, bn); 1842 if (cmp > 0) 1843 { 1844 rp = MPZ_REALLOC (r, an); 1845 gmp_assert_nocarry (mpn_sub (rp, a->_mp_d, an, b->_mp_d, bn)); 1846 return mpn_normalized_size (rp, an); 1847 } 1848 else if (cmp < 0) 1849 { 1850 rp = MPZ_REALLOC (r, bn); 1851 gmp_assert_nocarry (mpn_sub (rp, b->_mp_d, bn, a->_mp_d, an)); 1852 return -mpn_normalized_size (rp, bn); 1853 } 1854 else 1855 return 0; 1856 } 1857 1858 void 1859 mpz_add (mpz_t r, const mpz_t a, const mpz_t b) 1860 { 1861 mp_size_t rn; 1862 1863 if ( (a->_mp_size ^ b->_mp_size) >= 0) 1864 rn = mpz_abs_add (r, a, b); 1865 else 1866 rn = mpz_abs_sub (r, a, b); 1867 1868 r->_mp_size = a->_mp_size >= 0 ? rn : - rn; 1869 } 1870 1871 void 1872 mpz_sub (mpz_t r, const mpz_t a, const mpz_t b) 1873 { 1874 mp_size_t rn; 1875 1876 if ( (a->_mp_size ^ b->_mp_size) >= 0) 1877 rn = mpz_abs_sub (r, a, b); 1878 else 1879 rn = mpz_abs_add (r, a, b); 1880 1881 r->_mp_size = a->_mp_size >= 0 ? rn : - rn; 1882 } 1883 1884 1885 /* MPZ multiplication */ 1886 void 1887 mpz_mul_si (mpz_t r, const mpz_t u, long int v) 1888 { 1889 if (v < 0) 1890 { 1891 mpz_mul_ui (r, u, GMP_NEG_CAST (unsigned long int, v)); 1892 mpz_neg (r, r); 1893 } 1894 else 1895 mpz_mul_ui (r, u, (unsigned long int) v); 1896 } 1897 1898 void 1899 mpz_mul_ui (mpz_t r, const mpz_t u, unsigned long int v) 1900 { 1901 mp_size_t un; 1902 mpz_t t; 1903 mp_ptr tp; 1904 mp_limb_t cy; 1905 1906 un = GMP_ABS (u->_mp_size); 1907 1908 if (un == 0 || v == 0) 1909 { 1910 r->_mp_size = 0; 1911 return; 1912 } 1913 1914 mpz_init2 (t, (un + 1) * GMP_LIMB_BITS); 1915 1916 tp = t->_mp_d; 1917 cy = mpn_mul_1 (tp, u->_mp_d, un, v); 1918 tp[un] = cy; 1919 1920 t->_mp_size = un + (cy > 0); 1921 if (u->_mp_size < 0) 1922 t->_mp_size = - t->_mp_size; 1923 1924 mpz_swap (r, t); 1925 mpz_clear (t); 1926 } 1927 1928 void 1929 mpz_mul (mpz_t r, const mpz_t u, const mpz_t v) 1930 { 1931 int sign; 1932 mp_size_t un, vn, rn; 1933 mpz_t t; 1934 mp_ptr tp; 1935 1936 un = GMP_ABS (u->_mp_size); 1937 vn = GMP_ABS (v->_mp_size); 1938 1939 if (un == 0 || vn == 0) 1940 { 1941 r->_mp_size = 0; 1942 return; 1943 } 1944 1945 sign = (u->_mp_size ^ v->_mp_size) < 0; 1946 1947 mpz_init2 (t, (un + vn) * GMP_LIMB_BITS); 1948 1949 tp = t->_mp_d; 1950 if (un >= vn) 1951 mpn_mul (tp, u->_mp_d, un, v->_mp_d, vn); 1952 else 1953 mpn_mul (tp, v->_mp_d, vn, u->_mp_d, un); 1954 1955 rn = un + vn; 1956 rn -= tp[rn-1] == 0; 1957 1958 t->_mp_size = sign ? - rn : rn; 1959 mpz_swap (r, t); 1960 mpz_clear (t); 1961 } 1962 1963 void 1964 mpz_mul_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bits) 1965 { 1966 mp_size_t un, rn; 1967 mp_size_t limbs; 1968 unsigned shift; 1969 mp_ptr rp; 1970 1971 un = GMP_ABS (u->_mp_size); 1972 if (un == 0) 1973 { 1974 r->_mp_size = 0; 1975 return; 1976 } 1977 1978 limbs = bits / GMP_LIMB_BITS; 1979 shift = bits % GMP_LIMB_BITS; 1980 1981 rn = un + limbs + (shift > 0); 1982 rp = MPZ_REALLOC (r, rn); 1983 if (shift > 0) 1984 { 1985 mp_limb_t cy = mpn_lshift (rp + limbs, u->_mp_d, un, shift); 1986 rp[rn-1] = cy; 1987 rn -= (cy == 0); 1988 } 1989 else 1990 mpn_copyd (rp + limbs, u->_mp_d, un); 1991 1992 while (limbs > 0) 1993 rp[--limbs] = 0; 1994 1995 r->_mp_size = (u->_mp_size < 0) ? - rn : rn; 1996 } 1997 1998 1999 /* MPZ division */ 2000 enum mpz_div_round_mode { GMP_DIV_FLOOR, GMP_DIV_CEIL, GMP_DIV_TRUNC }; 2001 2002 /* Allows q or r to be zero. Returns 1 iff remainder is non-zero. */ 2003 static int 2004 mpz_div_qr (mpz_t q, mpz_t r, 2005 const mpz_t n, const mpz_t d, enum mpz_div_round_mode mode) 2006 { 2007 mp_size_t ns, ds, nn, dn, qs; 2008 ns = n->_mp_size; 2009 ds = d->_mp_size; 2010 2011 if (ds == 0) 2012 gmp_die("mpz_div_qr: Divide by zero."); 2013 2014 if (ns == 0) 2015 { 2016 if (q) 2017 q->_mp_size = 0; 2018 if (r) 2019 r->_mp_size = 0; 2020 return 0; 2021 } 2022 2023 nn = GMP_ABS (ns); 2024 dn = GMP_ABS (ds); 2025 2026 qs = ds ^ ns; 2027 2028 if (nn < dn) 2029 { 2030 if (mode == GMP_DIV_CEIL && qs >= 0) 2031 { 2032 /* q = 1, r = n - d */ 2033 if (r) 2034 mpz_sub (r, n, d); 2035 if (q) 2036 mpz_set_ui (q, 1); 2037 } 2038 else if (mode == GMP_DIV_FLOOR && qs < 0) 2039 { 2040 /* q = -1, r = n + d */ 2041 if (r) 2042 mpz_add (r, n, d); 2043 if (q) 2044 mpz_set_si (q, -1); 2045 } 2046 else 2047 { 2048 /* q = 0, r = d */ 2049 if (r) 2050 mpz_set (r, n); 2051 if (q) 2052 q->_mp_size = 0; 2053 } 2054 return 1; 2055 } 2056 else 2057 { 2058 mp_ptr np, qp; 2059 mp_size_t qn, rn; 2060 mpz_t tq, tr; 2061 2062 mpz_init (tr); 2063 mpz_set (tr, n); 2064 np = tr->_mp_d; 2065 2066 qn = nn - dn + 1; 2067 2068 if (q) 2069 { 2070 mpz_init2 (tq, qn * GMP_LIMB_BITS); 2071 qp = tq->_mp_d; 2072 } 2073 else 2074 qp = NULL; 2075 2076 mpn_div_qr (qp, np, nn, d->_mp_d, dn); 2077 2078 if (qp) 2079 { 2080 qn -= (qp[qn-1] == 0); 2081 2082 tq->_mp_size = qs < 0 ? -qn : qn; 2083 } 2084 rn = mpn_normalized_size (np, dn); 2085 tr->_mp_size = ns < 0 ? - rn : rn; 2086 2087 if (mode == GMP_DIV_FLOOR && qs < 0 && rn != 0) 2088 { 2089 if (q) 2090 mpz_sub_ui (tq, tq, 1); 2091 if (r) 2092 mpz_add (tr, tr, d); 2093 } 2094 else if (mode == GMP_DIV_CEIL && qs >= 0 && rn != 0) 2095 { 2096 if (q) 2097 mpz_add_ui (tq, tq, 1); 2098 if (r) 2099 mpz_sub (tr, tr, d); 2100 } 2101 2102 if (q) 2103 { 2104 mpz_swap (tq, q); 2105 mpz_clear (tq); 2106 } 2107 if (r) 2108 mpz_swap (tr, r); 2109 2110 mpz_clear (tr); 2111 2112 return rn != 0; 2113 } 2114 } 2115 2116 void 2117 mpz_cdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d) 2118 { 2119 mpz_div_qr (q, r, n, d, GMP_DIV_CEIL); 2120 } 2121 2122 void 2123 mpz_fdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d) 2124 { 2125 mpz_div_qr (q, r, n, d, GMP_DIV_FLOOR); 2126 } 2127 2128 void 2129 mpz_tdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d) 2130 { 2131 mpz_div_qr (q, r, n, d, GMP_DIV_TRUNC); 2132 } 2133 2134 void 2135 mpz_cdiv_q (mpz_t q, const mpz_t n, const mpz_t d) 2136 { 2137 mpz_div_qr (q, NULL, n, d, GMP_DIV_CEIL); 2138 } 2139 2140 void 2141 mpz_fdiv_q (mpz_t q, const mpz_t n, const mpz_t d) 2142 { 2143 mpz_div_qr (q, NULL, n, d, GMP_DIV_FLOOR); 2144 } 2145 2146 void 2147 mpz_tdiv_q (mpz_t q, const mpz_t n, const mpz_t d) 2148 { 2149 mpz_div_qr (q, NULL, n, d, GMP_DIV_TRUNC); 2150 } 2151 2152 void 2153 mpz_cdiv_r (mpz_t r, const mpz_t n, const mpz_t d) 2154 { 2155 mpz_div_qr (NULL, r, n, d, GMP_DIV_CEIL); 2156 } 2157 2158 void 2159 mpz_fdiv_r (mpz_t r, const mpz_t n, const mpz_t d) 2160 { 2161 mpz_div_qr (NULL, r, n, d, GMP_DIV_FLOOR); 2162 } 2163 2164 void 2165 mpz_tdiv_r (mpz_t r, const mpz_t n, const mpz_t d) 2166 { 2167 mpz_div_qr (NULL, r, n, d, GMP_DIV_TRUNC); 2168 } 2169 2170 void 2171 mpz_mod (mpz_t r, const mpz_t n, const mpz_t d) 2172 { 2173 if (d->_mp_size >= 0) 2174 mpz_div_qr (NULL, r, n, d, GMP_DIV_FLOOR); 2175 else 2176 mpz_div_qr (NULL, r, n, d, GMP_DIV_CEIL); 2177 } 2178 2179 static void 2180 mpz_div_q_2exp (mpz_t q, const mpz_t u, mp_bitcnt_t bit_index, 2181 enum mpz_div_round_mode mode) 2182 { 2183 mp_size_t un, qn; 2184 mp_size_t limb_cnt; 2185 mp_ptr qp; 2186 mp_limb_t adjust; 2187 2188 un = u->_mp_size; 2189 if (un == 0) 2190 { 2191 q->_mp_size = 0; 2192 return; 2193 } 2194 limb_cnt = bit_index / GMP_LIMB_BITS; 2195 qn = GMP_ABS (un) - limb_cnt; 2196 bit_index %= GMP_LIMB_BITS; 2197 2198 if (mode == ((un > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* un != 0 here. */ 2199 /* Note: Below, the final indexing at limb_cnt is valid because at 2200 that point we have qn > 0. */ 2201 adjust = (qn <= 0 2202 || !mpn_zero_p (u->_mp_d, limb_cnt) 2203 || (u->_mp_d[limb_cnt] 2204 & (((mp_limb_t) 1 << bit_index) - 1))); 2205 else 2206 adjust = 0; 2207 2208 if (qn <= 0) 2209 qn = 0; 2210 2211 else 2212 { 2213 qp = MPZ_REALLOC (q, qn); 2214 2215 if (bit_index != 0) 2216 { 2217 mpn_rshift (qp, u->_mp_d + limb_cnt, qn, bit_index); 2218 qn -= qp[qn - 1] == 0; 2219 } 2220 else 2221 { 2222 mpn_copyi (qp, u->_mp_d + limb_cnt, qn); 2223 } 2224 } 2225 2226 q->_mp_size = qn; 2227 2228 mpz_add_ui (q, q, adjust); 2229 if (un < 0) 2230 mpz_neg (q, q); 2231 } 2232 2233 static void 2234 mpz_div_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bit_index, 2235 enum mpz_div_round_mode mode) 2236 { 2237 mp_size_t us, un, rn; 2238 mp_ptr rp; 2239 mp_limb_t mask; 2240 2241 us = u->_mp_size; 2242 if (us == 0 || bit_index == 0) 2243 { 2244 r->_mp_size = 0; 2245 return; 2246 } 2247 rn = (bit_index + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS; 2248 assert (rn > 0); 2249 2250 rp = MPZ_REALLOC (r, rn); 2251 un = GMP_ABS (us); 2252 2253 mask = GMP_LIMB_MAX >> (rn * GMP_LIMB_BITS - bit_index); 2254 2255 if (rn > un) 2256 { 2257 /* Quotient (with truncation) is zero, and remainder is 2258 non-zero */ 2259 if (mode == ((us > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* us != 0 here. */ 2260 { 2261 /* Have to negate and sign extend. */ 2262 mp_size_t i; 2263 mp_limb_t cy; 2264 2265 for (cy = 1, i = 0; i < un; i++) 2266 { 2267 mp_limb_t s = ~u->_mp_d[i] + cy; 2268 cy = s < cy; 2269 rp[i] = s; 2270 } 2271 assert (cy == 0); 2272 for (; i < rn - 1; i++) 2273 rp[i] = GMP_LIMB_MAX; 2274 2275 rp[rn-1] = mask; 2276 us = -us; 2277 } 2278 else 2279 { 2280 /* Just copy */ 2281 if (r != u) 2282 mpn_copyi (rp, u->_mp_d, un); 2283 2284 rn = un; 2285 } 2286 } 2287 else 2288 { 2289 if (r != u) 2290 mpn_copyi (rp, u->_mp_d, rn - 1); 2291 2292 rp[rn-1] = u->_mp_d[rn-1] & mask; 2293 2294 if (mode == ((us > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* us != 0 here. */ 2295 { 2296 /* If r != 0, compute 2^{bit_count} - r. */ 2297 mp_size_t i; 2298 2299 for (i = 0; i < rn && rp[i] == 0; i++) 2300 ; 2301 if (i < rn) 2302 { 2303 /* r > 0, need to flip sign. */ 2304 rp[i] = ~rp[i] + 1; 2305 for (i++; i < rn; i++) 2306 rp[i] = ~rp[i]; 2307 2308 rp[rn-1] &= mask; 2309 2310 /* us is not used for anything else, so we can modify it 2311 here to indicate flipped sign. */ 2312 us = -us; 2313 } 2314 } 2315 } 2316 rn = mpn_normalized_size (rp, rn); 2317 r->_mp_size = us < 0 ? -rn : rn; 2318 } 2319 2320 void 2321 mpz_cdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt) 2322 { 2323 mpz_div_q_2exp (r, u, cnt, GMP_DIV_CEIL); 2324 } 2325 2326 void 2327 mpz_fdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt) 2328 { 2329 mpz_div_q_2exp (r, u, cnt, GMP_DIV_FLOOR); 2330 } 2331 2332 void 2333 mpz_tdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt) 2334 { 2335 mpz_div_q_2exp (r, u, cnt, GMP_DIV_TRUNC); 2336 } 2337 2338 void 2339 mpz_cdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt) 2340 { 2341 mpz_div_r_2exp (r, u, cnt, GMP_DIV_CEIL); 2342 } 2343 2344 void 2345 mpz_fdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt) 2346 { 2347 mpz_div_r_2exp (r, u, cnt, GMP_DIV_FLOOR); 2348 } 2349 2350 void 2351 mpz_tdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt) 2352 { 2353 mpz_div_r_2exp (r, u, cnt, GMP_DIV_TRUNC); 2354 } 2355 2356 void 2357 mpz_divexact (mpz_t q, const mpz_t n, const mpz_t d) 2358 { 2359 gmp_assert_nocarry (mpz_div_qr (q, NULL, n, d, GMP_DIV_TRUNC)); 2360 } 2361 2362 int 2363 mpz_divisible_p (const mpz_t n, const mpz_t d) 2364 { 2365 return mpz_div_qr (NULL, NULL, n, d, GMP_DIV_TRUNC) == 0; 2366 } 2367 2368 static unsigned long 2369 mpz_div_qr_ui (mpz_t q, mpz_t r, 2370 const mpz_t n, unsigned long d, enum mpz_div_round_mode mode) 2371 { 2372 mp_size_t ns, qn; 2373 mp_ptr qp; 2374 mp_limb_t rl; 2375 mp_size_t rs; 2376 2377 ns = n->_mp_size; 2378 if (ns == 0) 2379 { 2380 if (q) 2381 q->_mp_size = 0; 2382 if (r) 2383 r->_mp_size = 0; 2384 return 0; 2385 } 2386 2387 qn = GMP_ABS (ns); 2388 if (q) 2389 qp = MPZ_REALLOC (q, qn); 2390 else 2391 qp = NULL; 2392 2393 rl = mpn_div_qr_1 (qp, n->_mp_d, qn, d); 2394 assert (rl < d); 2395 2396 rs = rl > 0; 2397 rs = (ns < 0) ? -rs : rs; 2398 2399 if (rl > 0 && ( (mode == GMP_DIV_FLOOR && ns < 0) 2400 || (mode == GMP_DIV_CEIL && ns >= 0))) 2401 { 2402 if (q) 2403 gmp_assert_nocarry (mpn_add_1 (qp, qp, qn, 1)); 2404 rl = d - rl; 2405 rs = -rs; 2406 } 2407 2408 if (r) 2409 { 2410 r->_mp_d[0] = rl; 2411 r->_mp_size = rs; 2412 } 2413 if (q) 2414 { 2415 qn -= (qp[qn-1] == 0); 2416 assert (qn == 0 || qp[qn-1] > 0); 2417 2418 q->_mp_size = (ns < 0) ? - qn : qn; 2419 } 2420 2421 return rl; 2422 } 2423 2424 unsigned long 2425 mpz_cdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d) 2426 { 2427 return mpz_div_qr_ui (q, r, n, d, GMP_DIV_CEIL); 2428 } 2429 2430 unsigned long 2431 mpz_fdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d) 2432 { 2433 return mpz_div_qr_ui (q, r, n, d, GMP_DIV_FLOOR); 2434 } 2435 2436 unsigned long 2437 mpz_tdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d) 2438 { 2439 return mpz_div_qr_ui (q, r, n, d, GMP_DIV_TRUNC); 2440 } 2441 2442 unsigned long 2443 mpz_cdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d) 2444 { 2445 return mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_CEIL); 2446 } 2447 2448 unsigned long 2449 mpz_fdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d) 2450 { 2451 return mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_FLOOR); 2452 } 2453 2454 unsigned long 2455 mpz_tdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d) 2456 { 2457 return mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_TRUNC); 2458 } 2459 2460 unsigned long 2461 mpz_cdiv_r_ui (mpz_t r, const mpz_t n, unsigned long d) 2462 { 2463 return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_CEIL); 2464 } 2465 unsigned long 2466 mpz_fdiv_r_ui (mpz_t r, const mpz_t n, unsigned long d) 2467 { 2468 return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_FLOOR); 2469 } 2470 unsigned long 2471 mpz_tdiv_r_ui (mpz_t r, const mpz_t n, unsigned long d) 2472 { 2473 return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_TRUNC); 2474 } 2475 2476 unsigned long 2477 mpz_cdiv_ui (const mpz_t n, unsigned long d) 2478 { 2479 return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_CEIL); 2480 } 2481 2482 unsigned long 2483 mpz_fdiv_ui (const mpz_t n, unsigned long d) 2484 { 2485 return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_FLOOR); 2486 } 2487 2488 unsigned long 2489 mpz_tdiv_ui (const mpz_t n, unsigned long d) 2490 { 2491 return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_TRUNC); 2492 } 2493 2494 unsigned long 2495 mpz_mod_ui (mpz_t r, const mpz_t n, unsigned long d) 2496 { 2497 return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_FLOOR); 2498 } 2499 2500 void 2501 mpz_divexact_ui (mpz_t q, const mpz_t n, unsigned long d) 2502 { 2503 gmp_assert_nocarry (mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_TRUNC)); 2504 } 2505 2506 int 2507 mpz_divisible_ui_p (const mpz_t n, unsigned long d) 2508 { 2509 return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_TRUNC) == 0; 2510 } 2511 2512 2513 /* GCD */ 2514 static mp_limb_t 2515 mpn_gcd_11 (mp_limb_t u, mp_limb_t v) 2516 { 2517 unsigned shift; 2518 2519 assert ( (u | v) > 0); 2520 2521 if (u == 0) 2522 return v; 2523 else if (v == 0) 2524 return u; 2525 2526 gmp_ctz (shift, u | v); 2527 2528 u >>= shift; 2529 v >>= shift; 2530 2531 if ( (u & 1) == 0) 2532 MP_LIMB_T_SWAP (u, v); 2533 2534 while ( (v & 1) == 0) 2535 v >>= 1; 2536 2537 while (u != v) 2538 { 2539 if (u > v) 2540 { 2541 u -= v; 2542 do 2543 u >>= 1; 2544 while ( (u & 1) == 0); 2545 } 2546 else 2547 { 2548 v -= u; 2549 do 2550 v >>= 1; 2551 while ( (v & 1) == 0); 2552 } 2553 } 2554 return u << shift; 2555 } 2556 2557 unsigned long 2558 mpz_gcd_ui (mpz_t g, const mpz_t u, unsigned long v) 2559 { 2560 mp_size_t un; 2561 2562 if (v == 0) 2563 { 2564 if (g) 2565 mpz_abs (g, u); 2566 } 2567 else 2568 { 2569 un = GMP_ABS (u->_mp_size); 2570 if (un != 0) 2571 v = mpn_gcd_11 (mpn_div_qr_1 (NULL, u->_mp_d, un, v), v); 2572 2573 if (g) 2574 mpz_set_ui (g, v); 2575 } 2576 2577 return v; 2578 } 2579 2580 static mp_bitcnt_t 2581 mpz_make_odd (mpz_t r, const mpz_t u) 2582 { 2583 mp_size_t un, rn, i; 2584 mp_ptr rp; 2585 unsigned shift; 2586 2587 un = GMP_ABS (u->_mp_size); 2588 assert (un > 0); 2589 2590 for (i = 0; u->_mp_d[i] == 0; i++) 2591 ; 2592 2593 gmp_ctz (shift, u->_mp_d[i]); 2594 2595 rn = un - i; 2596 rp = MPZ_REALLOC (r, rn); 2597 if (shift > 0) 2598 { 2599 mpn_rshift (rp, u->_mp_d + i, rn, shift); 2600 rn -= (rp[rn-1] == 0); 2601 } 2602 else 2603 mpn_copyi (rp, u->_mp_d + i, rn); 2604 2605 r->_mp_size = rn; 2606 return i * GMP_LIMB_BITS + shift; 2607 } 2608 2609 void 2610 mpz_gcd (mpz_t g, const mpz_t u, const mpz_t v) 2611 { 2612 mpz_t tu, tv; 2613 mp_bitcnt_t uz, vz, gz; 2614 2615 if (u->_mp_size == 0) 2616 { 2617 mpz_abs (g, v); 2618 return; 2619 } 2620 if (v->_mp_size == 0) 2621 { 2622 mpz_abs (g, u); 2623 return; 2624 } 2625 2626 mpz_init (tu); 2627 mpz_init (tv); 2628 2629 uz = mpz_make_odd (tu, u); 2630 vz = mpz_make_odd (tv, v); 2631 gz = GMP_MIN (uz, vz); 2632 2633 if (tu->_mp_size < tv->_mp_size) 2634 mpz_swap (tu, tv); 2635 2636 mpz_tdiv_r (tu, tu, tv); 2637 if (tu->_mp_size == 0) 2638 { 2639 mpz_swap (g, tv); 2640 } 2641 else 2642 for (;;) 2643 { 2644 int c; 2645 2646 mpz_make_odd (tu, tu); 2647 c = mpz_cmp (tu, tv); 2648 if (c == 0) 2649 { 2650 mpz_swap (g, tu); 2651 break; 2652 } 2653 if (c < 0) 2654 mpz_swap (tu, tv); 2655 2656 if (tv->_mp_size == 1) 2657 { 2658 mp_limb_t vl = tv->_mp_d[0]; 2659 mp_limb_t ul = mpz_tdiv_ui (tu, vl); 2660 mpz_set_ui (g, mpn_gcd_11 (ul, vl)); 2661 break; 2662 } 2663 mpz_sub (tu, tu, tv); 2664 } 2665 mpz_clear (tu); 2666 mpz_clear (tv); 2667 mpz_mul_2exp (g, g, gz); 2668 } 2669 2670 void 2671 mpz_gcdext (mpz_t g, mpz_t s, mpz_t t, const mpz_t u, const mpz_t v) 2672 { 2673 mpz_t tu, tv, s0, s1, t0, t1; 2674 mp_bitcnt_t uz, vz, gz; 2675 mp_bitcnt_t power; 2676 2677 if (u->_mp_size == 0) 2678 { 2679 /* g = 0 u + sgn(v) v */ 2680 signed long sign = mpz_sgn (v); 2681 mpz_abs (g, v); 2682 if (s) 2683 mpz_set_ui (s, 0); 2684 if (t) 2685 mpz_set_si (t, sign); 2686 return; 2687 } 2688 2689 if (v->_mp_size == 0) 2690 { 2691 /* g = sgn(u) u + 0 v */ 2692 signed long sign = mpz_sgn (u); 2693 mpz_abs (g, u); 2694 if (s) 2695 mpz_set_si (s, sign); 2696 if (t) 2697 mpz_set_ui (t, 0); 2698 return; 2699 } 2700 2701 mpz_init (tu); 2702 mpz_init (tv); 2703 mpz_init (s0); 2704 mpz_init (s1); 2705 mpz_init (t0); 2706 mpz_init (t1); 2707 2708 uz = mpz_make_odd (tu, u); 2709 vz = mpz_make_odd (tv, v); 2710 gz = GMP_MIN (uz, vz); 2711 2712 uz -= gz; 2713 vz -= gz; 2714 2715 /* Cofactors corresponding to odd gcd. gz handled later. */ 2716 if (tu->_mp_size < tv->_mp_size) 2717 { 2718 mpz_swap (tu, tv); 2719 MPZ_SRCPTR_SWAP (u, v); 2720 MPZ_PTR_SWAP (s, t); 2721 MP_BITCNT_T_SWAP (uz, vz); 2722 } 2723 2724 /* Maintain 2725 * 2726 * u = t0 tu + t1 tv 2727 * v = s0 tu + s1 tv 2728 * 2729 * where u and v denote the inputs with common factors of two 2730 * eliminated, and det (s0, t0; s1, t1) = 2^p. Then 2731 * 2732 * 2^p tu = s1 u - t1 v 2733 * 2^p tv = -s0 u + t0 v 2734 */ 2735 2736 /* After initial division, tu = q tv + tu', we have 2737 * 2738 * u = 2^uz (tu' + q tv) 2739 * v = 2^vz tv 2740 * 2741 * or 2742 * 2743 * t0 = 2^uz, t1 = 2^uz q 2744 * s0 = 0, s1 = 2^vz 2745 */ 2746 2747 mpz_setbit (t0, uz); 2748 mpz_tdiv_qr (t1, tu, tu, tv); 2749 mpz_mul_2exp (t1, t1, uz); 2750 2751 mpz_setbit (s1, vz); 2752 power = uz + vz; 2753 2754 if (tu->_mp_size > 0) 2755 { 2756 mp_bitcnt_t shift; 2757 shift = mpz_make_odd (tu, tu); 2758 mpz_mul_2exp (t0, t0, shift); 2759 mpz_mul_2exp (s0, s0, shift); 2760 power += shift; 2761 2762 for (;;) 2763 { 2764 int c; 2765 c = mpz_cmp (tu, tv); 2766 if (c == 0) 2767 break; 2768 2769 if (c < 0) 2770 { 2771 /* tv = tv' + tu 2772 * 2773 * u = t0 tu + t1 (tv' + tu) = (t0 + t1) tu + t1 tv' 2774 * v = s0 tu + s1 (tv' + tu) = (s0 + s1) tu + s1 tv' */ 2775 2776 mpz_sub (tv, tv, tu); 2777 mpz_add (t0, t0, t1); 2778 mpz_add (s0, s0, s1); 2779 2780 shift = mpz_make_odd (tv, tv); 2781 mpz_mul_2exp (t1, t1, shift); 2782 mpz_mul_2exp (s1, s1, shift); 2783 } 2784 else 2785 { 2786 mpz_sub (tu, tu, tv); 2787 mpz_add (t1, t0, t1); 2788 mpz_add (s1, s0, s1); 2789 2790 shift = mpz_make_odd (tu, tu); 2791 mpz_mul_2exp (t0, t0, shift); 2792 mpz_mul_2exp (s0, s0, shift); 2793 } 2794 power += shift; 2795 } 2796 } 2797 2798 /* Now tv = odd part of gcd, and -s0 and t0 are corresponding 2799 cofactors. */ 2800 2801 mpz_mul_2exp (tv, tv, gz); 2802 mpz_neg (s0, s0); 2803 2804 /* 2^p g = s0 u + t0 v. Eliminate one factor of two at a time. To 2805 adjust cofactors, we need u / g and v / g */ 2806 2807 mpz_divexact (s1, v, tv); 2808 mpz_abs (s1, s1); 2809 mpz_divexact (t1, u, tv); 2810 mpz_abs (t1, t1); 2811 2812 while (power-- > 0) 2813 { 2814 /* s0 u + t0 v = (s0 - v/g) u - (t0 + u/g) v */ 2815 if (mpz_odd_p (s0) || mpz_odd_p (t0)) 2816 { 2817 mpz_sub (s0, s0, s1); 2818 mpz_add (t0, t0, t1); 2819 } 2820 mpz_divexact_ui (s0, s0, 2); 2821 mpz_divexact_ui (t0, t0, 2); 2822 } 2823 2824 /* Arrange so that |s| < |u| / 2g */ 2825 mpz_add (s1, s0, s1); 2826 if (mpz_cmpabs (s0, s1) > 0) 2827 { 2828 mpz_swap (s0, s1); 2829 mpz_sub (t0, t0, t1); 2830 } 2831 if (u->_mp_size < 0) 2832 mpz_neg (s0, s0); 2833 if (v->_mp_size < 0) 2834 mpz_neg (t0, t0); 2835 2836 mpz_swap (g, tv); 2837 if (s) 2838 mpz_swap (s, s0); 2839 if (t) 2840 mpz_swap (t, t0); 2841 2842 mpz_clear (tu); 2843 mpz_clear (tv); 2844 mpz_clear (s0); 2845 mpz_clear (s1); 2846 mpz_clear (t0); 2847 mpz_clear (t1); 2848 } 2849 2850 void 2851 mpz_lcm (mpz_t r, const mpz_t u, const mpz_t v) 2852 { 2853 mpz_t g; 2854 2855 if (u->_mp_size == 0 || v->_mp_size == 0) 2856 { 2857 r->_mp_size = 0; 2858 return; 2859 } 2860 2861 mpz_init (g); 2862 2863 mpz_gcd (g, u, v); 2864 mpz_divexact (g, u, g); 2865 mpz_mul (r, g, v); 2866 2867 mpz_clear (g); 2868 mpz_abs (r, r); 2869 } 2870 2871 void 2872 mpz_lcm_ui (mpz_t r, const mpz_t u, unsigned long v) 2873 { 2874 if (v == 0 || u->_mp_size == 0) 2875 { 2876 r->_mp_size = 0; 2877 return; 2878 } 2879 2880 v /= mpz_gcd_ui (NULL, u, v); 2881 mpz_mul_ui (r, u, v); 2882 2883 mpz_abs (r, r); 2884 } 2885 2886 int 2887 mpz_invert (mpz_t r, const mpz_t u, const mpz_t m) 2888 { 2889 mpz_t g, tr; 2890 int invertible; 2891 2892 if (u->_mp_size == 0 || mpz_cmpabs_ui (m, 1) <= 0) 2893 return 0; 2894 2895 mpz_init (g); 2896 mpz_init (tr); 2897 2898 mpz_gcdext (g, tr, NULL, u, m); 2899 invertible = (mpz_cmp_ui (g, 1) == 0); 2900 2901 if (invertible) 2902 { 2903 if (tr->_mp_size < 0) 2904 { 2905 if (m->_mp_size >= 0) 2906 mpz_add (tr, tr, m); 2907 else 2908 mpz_sub (tr, tr, m); 2909 } 2910 mpz_swap (r, tr); 2911 } 2912 2913 mpz_clear (g); 2914 mpz_clear (tr); 2915 return invertible; 2916 } 2917 2918 2919 /* Higher level operations (sqrt, pow and root) */ 2920 2921 void 2922 mpz_pow_ui (mpz_t r, const mpz_t b, unsigned long e) 2923 { 2924 unsigned long bit; 2925 mpz_t tr; 2926 mpz_init_set_ui (tr, 1); 2927 2928 for (bit = GMP_ULONG_HIGHBIT; bit > 0; bit >>= 1) 2929 { 2930 mpz_mul (tr, tr, tr); 2931 if (e & bit) 2932 mpz_mul (tr, tr, b); 2933 } 2934 mpz_swap (r, tr); 2935 mpz_clear (tr); 2936 } 2937 2938 void 2939 mpz_ui_pow_ui (mpz_t r, unsigned long blimb, unsigned long e) 2940 { 2941 mpz_t b; 2942 mpz_init_set_ui (b, blimb); 2943 mpz_pow_ui (r, b, e); 2944 mpz_clear (b); 2945 } 2946 2947 void 2948 mpz_powm (mpz_t r, const mpz_t b, const mpz_t e, const mpz_t m) 2949 { 2950 mpz_t tr; 2951 mpz_t base; 2952 mp_size_t en, mn; 2953 mp_srcptr mp; 2954 struct gmp_div_inverse minv; 2955 unsigned shift; 2956 mp_ptr tp = NULL; 2957 2958 en = GMP_ABS (e->_mp_size); 2959 mn = GMP_ABS (m->_mp_size); 2960 if (mn == 0) 2961 gmp_die ("mpz_powm: Zero modulo."); 2962 2963 if (en == 0) 2964 { 2965 mpz_set_ui (r, 1); 2966 return; 2967 } 2968 2969 mp = m->_mp_d; 2970 mpn_div_qr_invert (&minv, mp, mn); 2971 shift = minv.shift; 2972 2973 if (shift > 0) 2974 { 2975 /* To avoid shifts, we do all our reductions, except the final 2976 one, using a *normalized* m. */ 2977 minv.shift = 0; 2978 2979 tp = gmp_xalloc_limbs (mn); 2980 gmp_assert_nocarry (mpn_lshift (tp, mp, mn, shift)); 2981 mp = tp; 2982 } 2983 2984 mpz_init (base); 2985 2986 if (e->_mp_size < 0) 2987 { 2988 if (!mpz_invert (base, b, m)) 2989 gmp_die ("mpz_powm: Negative exponent and non-invertibe base."); 2990 } 2991 else 2992 { 2993 mp_size_t bn; 2994 mpz_abs (base, b); 2995 2996 bn = base->_mp_size; 2997 if (bn >= mn) 2998 { 2999 mpn_div_qr_preinv (NULL, base->_mp_d, base->_mp_size, mp, mn, &minv); 3000 bn = mn; 3001 } 3002 3003 /* We have reduced the absolute value. Now take care of the 3004 sign. Note that we get zero represented non-canonically as 3005 m. */ 3006 if (b->_mp_size < 0) 3007 { 3008 mp_ptr bp = MPZ_REALLOC (base, mn); 3009 gmp_assert_nocarry (mpn_sub (bp, mp, mn, bp, bn)); 3010 bn = mn; 3011 } 3012 base->_mp_size = mpn_normalized_size (base->_mp_d, bn); 3013 } 3014 mpz_init_set_ui (tr, 1); 3015 3016 while (en-- > 0) 3017 { 3018 mp_limb_t w = e->_mp_d[en]; 3019 mp_limb_t bit; 3020 3021 for (bit = GMP_LIMB_HIGHBIT; bit > 0; bit >>= 1) 3022 { 3023 mpz_mul (tr, tr, tr); 3024 if (w & bit) 3025 mpz_mul (tr, tr, base); 3026 if (tr->_mp_size > mn) 3027 { 3028 mpn_div_qr_preinv (NULL, tr->_mp_d, tr->_mp_size, mp, mn, &minv); 3029 tr->_mp_size = mpn_normalized_size (tr->_mp_d, mn); 3030 } 3031 } 3032 } 3033 3034 /* Final reduction */ 3035 if (tr->_mp_size >= mn) 3036 { 3037 minv.shift = shift; 3038 mpn_div_qr_preinv (NULL, tr->_mp_d, tr->_mp_size, mp, mn, &minv); 3039 tr->_mp_size = mpn_normalized_size (tr->_mp_d, mn); 3040 } 3041 if (tp) 3042 gmp_free (tp); 3043 3044 mpz_swap (r, tr); 3045 mpz_clear (tr); 3046 mpz_clear (base); 3047 } 3048 3049 void 3050 mpz_powm_ui (mpz_t r, const mpz_t b, unsigned long elimb, const mpz_t m) 3051 { 3052 mpz_t e; 3053 mpz_init_set_ui (e, elimb); 3054 mpz_powm (r, b, e, m); 3055 mpz_clear (e); 3056 } 3057 3058 /* x=trunc(y^(1/z)), r=y-x^z */ 3059 void 3060 mpz_rootrem (mpz_t x, mpz_t r, const mpz_t y, unsigned long z) 3061 { 3062 int sgn; 3063 mpz_t t, u; 3064 3065 sgn = y->_mp_size < 0; 3066 if (sgn && (z & 1) == 0) 3067 gmp_die ("mpz_rootrem: Negative argument, with even root."); 3068 if (z == 0) 3069 gmp_die ("mpz_rootrem: Zeroth root."); 3070 3071 if (mpz_cmpabs_ui (y, 1) <= 0) { 3072 mpz_set (x, y); 3073 if (r) 3074 r->_mp_size = 0; 3075 return; 3076 } 3077 3078 mpz_init (t); 3079 mpz_init (u); 3080 mpz_setbit (t, mpz_sizeinbase (y, 2) / z + 1); 3081 3082 if (z == 2) /* simplify sqrt loop: z-1 == 1 */ 3083 do { 3084 mpz_swap (u, t); /* u = x */ 3085 mpz_tdiv_q (t, y, u); /* t = y/x */ 3086 mpz_add (t, t, u); /* t = y/x + x */ 3087 mpz_tdiv_q_2exp (t, t, 1); /* x'= (y/x + x)/2 */ 3088 } while (mpz_cmpabs (t, u) < 0); /* |x'| < |x| */ 3089 else /* z != 2 */ { 3090 mpz_t v; 3091 3092 mpz_init (v); 3093 if (sgn) 3094 mpz_neg (t, t); 3095 3096 do { 3097 mpz_swap (u, t); /* u = x */ 3098 mpz_pow_ui (t, u, z - 1); /* t = x^(z-1) */ 3099 mpz_tdiv_q (t, y, t); /* t = y/x^(z-1) */ 3100 mpz_mul_ui (v, u, z - 1); /* v = x*(z-1) */ 3101 mpz_add (t, t, v); /* t = y/x^(z-1) + x*(z-1) */ 3102 mpz_tdiv_q_ui (t, t, z); /* x'=(y/x^(z-1) + x*(z-1))/z */ 3103 } while (mpz_cmpabs (t, u) < 0); /* |x'| < |x| */ 3104 3105 mpz_clear (v); 3106 } 3107 3108 if (r) { 3109 mpz_pow_ui (t, u, z); 3110 mpz_sub (r, y, t); 3111 } 3112 mpz_swap (x, u); 3113 mpz_clear (u); 3114 mpz_clear (t); 3115 } 3116 3117 int 3118 mpz_root (mpz_t x, const mpz_t y, unsigned long z) 3119 { 3120 int res; 3121 mpz_t r; 3122 3123 mpz_init (r); 3124 mpz_rootrem (x, r, y, z); 3125 res = r->_mp_size == 0; 3126 mpz_clear (r); 3127 3128 return res; 3129 } 3130 3131 /* Compute s = floor(sqrt(u)) and r = u - s^2. Allows r == NULL */ 3132 void 3133 mpz_sqrtrem (mpz_t s, mpz_t r, const mpz_t u) 3134 { 3135 mpz_rootrem (s, r, u, 2); 3136 } 3137 3138 void 3139 mpz_sqrt (mpz_t s, const mpz_t u) 3140 { 3141 mpz_rootrem (s, NULL, u, 2); 3142 } 3143 3144 3145 /* Combinatorics */ 3146 3147 void 3148 mpz_fac_ui (mpz_t x, unsigned long n) 3149 { 3150 if (n < 2) { 3151 mpz_set_ui (x, 1); 3152 return; 3153 } 3154 mpz_set_ui (x, n); 3155 for (;--n > 1;) 3156 mpz_mul_ui (x, x, n); 3157 } 3158 3159 void 3160 mpz_bin_uiui (mpz_t r, unsigned long n, unsigned long k) 3161 { 3162 mpz_t t; 3163 3164 if (k > n) { 3165 r->_mp_size = 0; 3166 return; 3167 } 3168 mpz_fac_ui (r, n); 3169 mpz_init (t); 3170 mpz_fac_ui (t, k); 3171 mpz_divexact (r, r, t); 3172 mpz_fac_ui (t, n - k); 3173 mpz_divexact (r, r, t); 3174 mpz_clear (t); 3175 } 3176 3177 3178 /* Logical operations and bit manipulation. */ 3179 3180 /* Numbers are treated as if represented in two's complement (and 3181 infinitely sign extended). For a negative values we get the two's 3182 complement from -x = ~x + 1, where ~ is bitwise complementt. 3183 Negation transforms 3184 3185 xxxx10...0 3186 3187 into 3188 3189 yyyy10...0 3190 3191 where yyyy is the bitwise complement of xxxx. So least significant 3192 bits, up to and including the first one bit, are unchanged, and 3193 the more significant bits are all complemented. 3194 3195 To change a bit from zero to one in a negative number, subtract the 3196 corresponding power of two from the absolute value. This can never 3197 underflow. To change a bit from one to zero, add the corresponding 3198 power of two, and this might overflow. E.g., if x = -001111, the 3199 two's complement is 110001. Clearing the least significant bit, we 3200 get two's complement 110000, and -010000. */ 3201 3202 int 3203 mpz_tstbit (const mpz_t d, mp_bitcnt_t bit_index) 3204 { 3205 mp_size_t limb_index; 3206 unsigned shift; 3207 mp_size_t ds; 3208 mp_size_t dn; 3209 mp_limb_t w; 3210 int bit; 3211 3212 ds = d->_mp_size; 3213 dn = GMP_ABS (ds); 3214 limb_index = bit_index / GMP_LIMB_BITS; 3215 if (limb_index >= dn) 3216 return ds < 0; 3217 3218 shift = bit_index % GMP_LIMB_BITS; 3219 w = d->_mp_d[limb_index]; 3220 bit = (w >> shift) & 1; 3221 3222 if (ds < 0) 3223 { 3224 /* d < 0. Check if any of the bits below is set: If so, our bit 3225 must be complemented. */ 3226 if (shift > 0 && (w << (GMP_LIMB_BITS - shift)) > 0) 3227 return bit ^ 1; 3228 while (limb_index-- > 0) 3229 if (d->_mp_d[limb_index] > 0) 3230 return bit ^ 1; 3231 } 3232 return bit; 3233 } 3234 3235 static void 3236 mpz_abs_add_bit (mpz_t d, mp_bitcnt_t bit_index) 3237 { 3238 mp_size_t dn, limb_index; 3239 mp_limb_t bit; 3240 mp_ptr dp; 3241 3242 dn = GMP_ABS (d->_mp_size); 3243 3244 limb_index = bit_index / GMP_LIMB_BITS; 3245 bit = (mp_limb_t) 1 << (bit_index % GMP_LIMB_BITS); 3246 3247 if (limb_index >= dn) 3248 { 3249 mp_size_t i; 3250 /* The bit should be set outside of the end of the number. 3251 We have to increase the size of the number. */ 3252 dp = MPZ_REALLOC (d, limb_index + 1); 3253 3254 dp[limb_index] = bit; 3255 for (i = dn; i < limb_index; i++) 3256 dp[i] = 0; 3257 dn = limb_index + 1; 3258 } 3259 else 3260 { 3261 mp_limb_t cy; 3262 3263 dp = d->_mp_d; 3264 3265 cy = mpn_add_1 (dp + limb_index, dp + limb_index, dn - limb_index, bit); 3266 if (cy > 0) 3267 { 3268 dp = MPZ_REALLOC (d, dn + 1); 3269 dp[dn++] = cy; 3270 } 3271 } 3272 3273 d->_mp_size = (d->_mp_size < 0) ? - dn : dn; 3274 } 3275 3276 static void 3277 mpz_abs_sub_bit (mpz_t d, mp_bitcnt_t bit_index) 3278 { 3279 mp_size_t dn, limb_index; 3280 mp_ptr dp; 3281 mp_limb_t bit; 3282 3283 dn = GMP_ABS (d->_mp_size); 3284 dp = d->_mp_d; 3285 3286 limb_index = bit_index / GMP_LIMB_BITS; 3287 bit = (mp_limb_t) 1 << (bit_index % GMP_LIMB_BITS); 3288 3289 assert (limb_index < dn); 3290 3291 gmp_assert_nocarry (mpn_sub_1 (dp + limb_index, dp + limb_index, 3292 dn - limb_index, bit)); 3293 dn -= (dp[dn-1] == 0); 3294 d->_mp_size = (d->_mp_size < 0) ? - dn : dn; 3295 } 3296 3297 void 3298 mpz_setbit (mpz_t d, mp_bitcnt_t bit_index) 3299 { 3300 if (!mpz_tstbit (d, bit_index)) 3301 { 3302 if (d->_mp_size >= 0) 3303 mpz_abs_add_bit (d, bit_index); 3304 else 3305 mpz_abs_sub_bit (d, bit_index); 3306 } 3307 } 3308 3309 void 3310 mpz_clrbit (mpz_t d, mp_bitcnt_t bit_index) 3311 { 3312 if (mpz_tstbit (d, bit_index)) 3313 { 3314 if (d->_mp_size >= 0) 3315 mpz_abs_sub_bit (d, bit_index); 3316 else 3317 mpz_abs_add_bit (d, bit_index); 3318 } 3319 } 3320 3321 void 3322 mpz_combit (mpz_t d, mp_bitcnt_t bit_index) 3323 { 3324 if (mpz_tstbit (d, bit_index) ^ (d->_mp_size < 0)) 3325 mpz_abs_sub_bit (d, bit_index); 3326 else 3327 mpz_abs_add_bit (d, bit_index); 3328 } 3329 3330 void 3331 mpz_com (mpz_t r, const mpz_t u) 3332 { 3333 mpz_neg (r, u); 3334 mpz_sub_ui (r, r, 1); 3335 } 3336 3337 void 3338 mpz_and (mpz_t r, const mpz_t u, const mpz_t v) 3339 { 3340 mp_size_t un, vn, rn, i; 3341 mp_ptr up, vp, rp; 3342 3343 mp_limb_t ux, vx, rx; 3344 mp_limb_t uc, vc, rc; 3345 mp_limb_t ul, vl, rl; 3346 3347 un = GMP_ABS (u->_mp_size); 3348 vn = GMP_ABS (v->_mp_size); 3349 if (un < vn) 3350 { 3351 MPZ_SRCPTR_SWAP (u, v); 3352 MP_SIZE_T_SWAP (un, vn); 3353 } 3354 if (vn == 0) 3355 { 3356 r->_mp_size = 0; 3357 return; 3358 } 3359 3360 uc = u->_mp_size < 0; 3361 vc = v->_mp_size < 0; 3362 rc = uc & vc; 3363 3364 ux = -uc; 3365 vx = -vc; 3366 rx = -rc; 3367 3368 /* If the smaller input is positive, higher limbs don't matter. */ 3369 rn = vx ? un : vn; 3370 3371 rp = MPZ_REALLOC (r, rn + rc); 3372 3373 up = u->_mp_d; 3374 vp = v->_mp_d; 3375 3376 for (i = 0; i < vn; i++) 3377 { 3378 ul = (up[i] ^ ux) + uc; 3379 uc = ul < uc; 3380 3381 vl = (vp[i] ^ vx) + vc; 3382 vc = vl < vc; 3383 3384 rl = ( (ul & vl) ^ rx) + rc; 3385 rc = rl < rc; 3386 rp[i] = rl; 3387 } 3388 assert (vc == 0); 3389 3390 for (; i < rn; i++) 3391 { 3392 ul = (up[i] ^ ux) + uc; 3393 uc = ul < uc; 3394 3395 rl = ( (ul & vx) ^ rx) + rc; 3396 rc = rl < rc; 3397 rp[i] = rl; 3398 } 3399 if (rc) 3400 rp[rn++] = rc; 3401 else 3402 rn = mpn_normalized_size (rp, rn); 3403 3404 r->_mp_size = rx ? -rn : rn; 3405 } 3406 3407 void 3408 mpz_ior (mpz_t r, const mpz_t u, const mpz_t v) 3409 { 3410 mp_size_t un, vn, rn, i; 3411 mp_ptr up, vp, rp; 3412 3413 mp_limb_t ux, vx, rx; 3414 mp_limb_t uc, vc, rc; 3415 mp_limb_t ul, vl, rl; 3416 3417 un = GMP_ABS (u->_mp_size); 3418 vn = GMP_ABS (v->_mp_size); 3419 if (un < vn) 3420 { 3421 MPZ_SRCPTR_SWAP (u, v); 3422 MP_SIZE_T_SWAP (un, vn); 3423 } 3424 if (vn == 0) 3425 { 3426 mpz_set (r, u); 3427 return; 3428 } 3429 3430 uc = u->_mp_size < 0; 3431 vc = v->_mp_size < 0; 3432 rc = uc | vc; 3433 3434 ux = -uc; 3435 vx = -vc; 3436 rx = -rc; 3437 3438 /* If the smaller input is negative, by sign extension higher limbs 3439 don't matter. */ 3440 rn = vx ? vn : un; 3441 3442 rp = MPZ_REALLOC (r, rn + rc); 3443 3444 up = u->_mp_d; 3445 vp = v->_mp_d; 3446 3447 for (i = 0; i < vn; i++) 3448 { 3449 ul = (up[i] ^ ux) + uc; 3450 uc = ul < uc; 3451 3452 vl = (vp[i] ^ vx) + vc; 3453 vc = vl < vc; 3454 3455 rl = ( (ul | vl) ^ rx) + rc; 3456 rc = rl < rc; 3457 rp[i] = rl; 3458 } 3459 assert (vc == 0); 3460 3461 for (; i < rn; i++) 3462 { 3463 ul = (up[i] ^ ux) + uc; 3464 uc = ul < uc; 3465 3466 rl = ( (ul | vx) ^ rx) + rc; 3467 rc = rl < rc; 3468 rp[i] = rl; 3469 } 3470 if (rc) 3471 rp[rn++] = rc; 3472 else 3473 rn = mpn_normalized_size (rp, rn); 3474 3475 r->_mp_size = rx ? -rn : rn; 3476 } 3477 3478 void 3479 mpz_xor (mpz_t r, const mpz_t u, const mpz_t v) 3480 { 3481 mp_size_t un, vn, i; 3482 mp_ptr up, vp, rp; 3483 3484 mp_limb_t ux, vx, rx; 3485 mp_limb_t uc, vc, rc; 3486 mp_limb_t ul, vl, rl; 3487 3488 un = GMP_ABS (u->_mp_size); 3489 vn = GMP_ABS (v->_mp_size); 3490 if (un < vn) 3491 { 3492 MPZ_SRCPTR_SWAP (u, v); 3493 MP_SIZE_T_SWAP (un, vn); 3494 } 3495 if (vn == 0) 3496 { 3497 mpz_set (r, u); 3498 return; 3499 } 3500 3501 uc = u->_mp_size < 0; 3502 vc = v->_mp_size < 0; 3503 rc = uc ^ vc; 3504 3505 ux = -uc; 3506 vx = -vc; 3507 rx = -rc; 3508 3509 rp = MPZ_REALLOC (r, un + rc); 3510 3511 up = u->_mp_d; 3512 vp = v->_mp_d; 3513 3514 for (i = 0; i < vn; i++) 3515 { 3516 ul = (up[i] ^ ux) + uc; 3517 uc = ul < uc; 3518 3519 vl = (vp[i] ^ vx) + vc; 3520 vc = vl < vc; 3521 3522 rl = (ul ^ vl ^ rx) + rc; 3523 rc = rl < rc; 3524 rp[i] = rl; 3525 } 3526 assert (vc == 0); 3527 3528 for (; i < un; i++) 3529 { 3530 ul = (up[i] ^ ux) + uc; 3531 uc = ul < uc; 3532 3533 rl = (ul ^ ux) + rc; 3534 rc = rl < rc; 3535 rp[i] = rl; 3536 } 3537 if (rc) 3538 rp[un++] = rc; 3539 else 3540 un = mpn_normalized_size (rp, un); 3541 3542 r->_mp_size = rx ? -un : un; 3543 } 3544 3545 static unsigned 3546 gmp_popcount_limb (mp_limb_t x) 3547 { 3548 unsigned c; 3549 3550 /* Do 16 bits at a time, to avoid limb-sized constants. */ 3551 for (c = 0; x > 0; x >>= 16) 3552 { 3553 unsigned w = ((x >> 1) & 0x5555) + (x & 0x5555); 3554 w = ((w >> 2) & 0x3333) + (w & 0x3333); 3555 w = ((w >> 4) & 0x0f0f) + (w & 0x0f0f); 3556 w = (w >> 8) + (w & 0x00ff); 3557 c += w; 3558 } 3559 return c; 3560 } 3561 3562 mp_bitcnt_t 3563 mpz_popcount (const mpz_t u) 3564 { 3565 mp_size_t un, i; 3566 mp_bitcnt_t c; 3567 3568 un = u->_mp_size; 3569 3570 if (un < 0) 3571 return ~(mp_bitcnt_t) 0; 3572 3573 for (c = 0, i = 0; i < un; i++) 3574 c += gmp_popcount_limb (u->_mp_d[i]); 3575 3576 return c; 3577 } 3578 3579 mp_bitcnt_t 3580 mpz_hamdist (const mpz_t u, const mpz_t v) 3581 { 3582 mp_size_t un, vn, i; 3583 mp_limb_t uc, vc, ul, vl, comp; 3584 mp_srcptr up, vp; 3585 mp_bitcnt_t c; 3586 3587 un = u->_mp_size; 3588 vn = v->_mp_size; 3589 3590 if ( (un ^ vn) < 0) 3591 return ~(mp_bitcnt_t) 0; 3592 3593 if (un < 0) 3594 { 3595 assert (vn < 0); 3596 un = -un; 3597 vn = -vn; 3598 uc = vc = 1; 3599 comp = - (mp_limb_t) 1; 3600 } 3601 else 3602 uc = vc = comp = 0; 3603 3604 up = u->_mp_d; 3605 vp = v->_mp_d; 3606 3607 if (un < vn) 3608 MPN_SRCPTR_SWAP (up, un, vp, vn); 3609 3610 for (i = 0, c = 0; i < vn; i++) 3611 { 3612 ul = (up[i] ^ comp) + uc; 3613 uc = ul < uc; 3614 3615 vl = (vp[i] ^ comp) + vc; 3616 vc = vl < vc; 3617 3618 c += gmp_popcount_limb (ul ^ vl); 3619 } 3620 assert (vc == 0); 3621 3622 for (; i < un; i++) 3623 { 3624 ul = (up[i] ^ comp) + uc; 3625 uc = ul < uc; 3626 3627 c += gmp_popcount_limb (ul ^ comp); 3628 } 3629 3630 return c; 3631 } 3632 3633 mp_bitcnt_t 3634 mpz_scan1 (const mpz_t u, mp_bitcnt_t starting_bit) 3635 { 3636 mp_ptr up; 3637 mp_size_t us, un, i; 3638 mp_limb_t limb, ux, uc; 3639 unsigned cnt; 3640 3641 up = u->_mp_d; 3642 us = u->_mp_size; 3643 un = GMP_ABS (us); 3644 i = starting_bit / GMP_LIMB_BITS; 3645 3646 /* Past the end there's no 1 bits for u>=0, or an immediate 1 bit 3647 for u<0. Notice this test picks up any u==0 too. */ 3648 if (i >= un) 3649 return (us >= 0 ? ~(mp_bitcnt_t) 0 : starting_bit); 3650 3651 if (us < 0) 3652 { 3653 ux = GMP_LIMB_MAX; 3654 uc = mpn_zero_p (up, i); 3655 } 3656 else 3657 ux = uc = 0; 3658 3659 limb = (ux ^ up[i]) + uc; 3660 uc = limb < uc; 3661 3662 /* Mask to 0 all bits before starting_bit, thus ignoring them. */ 3663 limb &= (GMP_LIMB_MAX << (starting_bit % GMP_LIMB_BITS)); 3664 3665 while (limb == 0) 3666 { 3667 i++; 3668 if (i == un) 3669 { 3670 assert (uc == 0); 3671 /* For the u > 0 case, this can happen only for the first 3672 masked limb. For the u < 0 case, it happens when the 3673 highest limbs of the absolute value are all ones. */ 3674 return (us >= 0 ? ~(mp_bitcnt_t) 0 : un * GMP_LIMB_BITS); 3675 } 3676 limb = (ux ^ up[i]) + uc; 3677 uc = limb < uc; 3678 } 3679 gmp_ctz (cnt, limb); 3680 return (mp_bitcnt_t) i * GMP_LIMB_BITS + cnt; 3681 } 3682 3683 mp_bitcnt_t 3684 mpz_scan0 (const mpz_t u, mp_bitcnt_t starting_bit) 3685 { 3686 mp_ptr up; 3687 mp_size_t us, un, i; 3688 mp_limb_t limb, ux, uc; 3689 unsigned cnt; 3690 3691 up = u->_mp_d; 3692 us = u->_mp_size; 3693 un = GMP_ABS (us); 3694 i = starting_bit / GMP_LIMB_BITS; 3695 3696 /* When past end, there's an immediate 0 bit for u>=0, or no 0 bits for 3697 u<0. Notice this test picks up all cases of u==0 too. */ 3698 if (i >= un) 3699 return (us >= 0 ? starting_bit : ~(mp_bitcnt_t) 0); 3700 3701 if (us < 0) 3702 { 3703 ux = GMP_LIMB_MAX; 3704 uc = mpn_zero_p (up, i); 3705 } 3706 else 3707 ux = uc = 0; 3708 3709 limb = (ux ^ up[i]) + uc; 3710 uc = limb < uc; 3711 3712 /* Mask to 1 all bits before starting_bit, thus ignoring them. */ 3713 limb |= ((mp_limb_t) 1 << (starting_bit % GMP_LIMB_BITS)) - 1; 3714 3715 while (limb == GMP_LIMB_MAX) 3716 { 3717 i++; 3718 if (i == un) 3719 { 3720 assert (uc == 0); 3721 return (us >= 0 ? un * GMP_LIMB_BITS : ~(mp_bitcnt_t) 0); 3722 } 3723 limb = (ux ^ up[i]) + uc; 3724 uc = limb < uc; 3725 } 3726 gmp_ctz (cnt, ~limb); 3727 return (mp_bitcnt_t) i * GMP_LIMB_BITS + cnt; 3728 } 3729 3730 3731 /* MPZ base conversion. */ 3732 3733 size_t 3734 mpz_sizeinbase (const mpz_t u, int base) 3735 { 3736 mp_size_t un; 3737 mp_srcptr up; 3738 mp_ptr tp; 3739 mp_bitcnt_t bits; 3740 struct gmp_div_inverse bi; 3741 size_t ndigits; 3742 3743 assert (base >= 2); 3744 assert (base <= 36); 3745 3746 un = GMP_ABS (u->_mp_size); 3747 if (un == 0) 3748 return 1; 3749 3750 up = u->_mp_d; 3751 3752 bits = (un - 1) * GMP_LIMB_BITS + mpn_limb_size_in_base_2 (up[un-1]); 3753 switch (base) 3754 { 3755 case 2: 3756 return bits; 3757 case 4: 3758 return (bits + 1) / 2; 3759 case 8: 3760 return (bits + 2) / 3; 3761 case 16: 3762 return (bits + 3) / 4; 3763 case 32: 3764 return (bits + 4) / 5; 3765 /* FIXME: Do something more clever for the common case of base 3766 10. */ 3767 } 3768 3769 tp = gmp_xalloc_limbs (un); 3770 mpn_copyi (tp, up, un); 3771 mpn_div_qr_1_invert (&bi, base); 3772 3773 for (ndigits = 0; un > 0; ndigits++) 3774 { 3775 mpn_div_qr_1_preinv (tp, tp, un, &bi); 3776 un -= (tp[un-1] == 0); 3777 } 3778 gmp_free (tp); 3779 return ndigits; 3780 } 3781 3782 char * 3783 mpz_get_str (char *sp, int base, const mpz_t u) 3784 { 3785 unsigned bits; 3786 const char *digits; 3787 mp_size_t un; 3788 size_t i, sn; 3789 3790 if (base >= 0) 3791 { 3792 digits = "0123456789abcdefghijklmnopqrstuvwxyz"; 3793 } 3794 else 3795 { 3796 base = -base; 3797 digits = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ"; 3798 } 3799 if (base <= 1) 3800 base = 10; 3801 if (base > 36) 3802 return NULL; 3803 3804 sn = 1 + mpz_sizeinbase (u, base); 3805 if (!sp) 3806 sp = gmp_xalloc (1 + sn); 3807 3808 un = GMP_ABS (u->_mp_size); 3809 3810 if (un == 0) 3811 { 3812 sp[0] = '0'; 3813 sp[1] = '\0'; 3814 return sp; 3815 } 3816 3817 i = 0; 3818 3819 if (u->_mp_size < 0) 3820 sp[i++] = '-'; 3821 3822 bits = mpn_base_power_of_two_p (base); 3823 3824 if (bits) 3825 /* Not modified in this case. */ 3826 sn = i + mpn_get_str_bits ((unsigned char *) sp + i, bits, u->_mp_d, un); 3827 else 3828 { 3829 struct mpn_base_info info; 3830 mp_ptr tp; 3831 3832 mpn_get_base_info (&info, base); 3833 tp = gmp_xalloc_limbs (un); 3834 mpn_copyi (tp, u->_mp_d, un); 3835 3836 sn = i + mpn_get_str_other ((unsigned char *) sp + i, base, &info, tp, un); 3837 gmp_free (tp); 3838 } 3839 3840 for (; i < sn; i++) 3841 sp[i] = digits[(unsigned char) sp[i]]; 3842 3843 sp[sn] = '\0'; 3844 return sp; 3845 } 3846 3847 int 3848 mpz_set_str (mpz_t r, const char *sp, int base) 3849 { 3850 unsigned bits; 3851 mp_size_t rn, alloc; 3852 mp_ptr rp; 3853 size_t sn; 3854 size_t dn; 3855 int sign; 3856 unsigned char *dp; 3857 3858 assert (base == 0 || (base >= 2 && base <= 36)); 3859 3860 while (isspace( (unsigned char) *sp)) 3861 sp++; 3862 3863 if (*sp == '-') 3864 { 3865 sign = 1; 3866 sp++; 3867 } 3868 else 3869 sign = 0; 3870 3871 if (base == 0) 3872 { 3873 if (*sp == '0') 3874 { 3875 sp++; 3876 if (*sp == 'x' || *sp == 'X') 3877 { 3878 base = 16; 3879 sp++; 3880 } 3881 else if (*sp == 'b' || *sp == 'B') 3882 { 3883 base = 2; 3884 sp++; 3885 } 3886 else 3887 base = 8; 3888 } 3889 else 3890 base = 10; 3891 } 3892 3893 sn = strlen (sp); 3894 dp = gmp_xalloc (sn + (sn == 0)); 3895 3896 for (dn = 0; *sp; sp++) 3897 { 3898 unsigned digit; 3899 3900 if (isspace ((unsigned char) *sp)) 3901 continue; 3902 if (*sp >= '0' && *sp <= '9') 3903 digit = *sp - '0'; 3904 else if (*sp >= 'a' && *sp <= 'z') 3905 digit = *sp - 'a' + 10; 3906 else if (*sp >= 'A' && *sp <= 'Z') 3907 digit = *sp - 'A' + 10; 3908 else 3909 digit = base; /* fail */ 3910 3911 if (digit >= base) 3912 { 3913 gmp_free (dp); 3914 r->_mp_size = 0; 3915 return -1; 3916 } 3917 3918 dp[dn++] = digit; 3919 } 3920 3921 bits = mpn_base_power_of_two_p (base); 3922 3923 if (bits > 0) 3924 { 3925 alloc = (sn * bits + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS; 3926 rp = MPZ_REALLOC (r, alloc); 3927 rn = mpn_set_str_bits (rp, dp, dn, bits); 3928 } 3929 else 3930 { 3931 struct mpn_base_info info; 3932 mpn_get_base_info (&info, base); 3933 alloc = (sn + info.exp - 1) / info.exp; 3934 rp = MPZ_REALLOC (r, alloc); 3935 rn = mpn_set_str_other (rp, dp, dn, base, &info); 3936 } 3937 assert (rn <= alloc); 3938 gmp_free (dp); 3939 3940 r->_mp_size = sign ? - rn : rn; 3941 3942 return 0; 3943 } 3944 3945 int 3946 mpz_init_set_str (mpz_t r, const char *sp, int base) 3947 { 3948 mpz_init (r); 3949 return mpz_set_str (r, sp, base); 3950 } 3951 3952 size_t 3953 mpz_out_str (FILE *stream, int base, const mpz_t x) 3954 { 3955 char *str; 3956 size_t len; 3957 3958 str = mpz_get_str (NULL, base, x); 3959 len = strlen (str); 3960 len = fwrite (str, 1, len, stream); 3961 gmp_free (str); 3962 return len; 3963 } 3964 3965 3966 static int 3967 gmp_detect_endian (void) 3968 { 3969 static const int i = 1; 3970 const unsigned char *p = (const unsigned char *) &i; 3971 if (*p == 1) 3972 /* Little endian */ 3973 return -1; 3974 else 3975 /* Big endian */ 3976 return 1; 3977 } 3978 3979 /* Import and export. Does not support nails. */ 3980 void 3981 mpz_import (mpz_t r, size_t count, int order, size_t size, int endian, 3982 size_t nails, const void *src) 3983 { 3984 const unsigned char *p; 3985 ptrdiff_t word_step; 3986 mp_ptr rp; 3987 mp_size_t rn; 3988 3989 /* The current (partial) limb. */ 3990 mp_limb_t limb; 3991 /* The number of bytes already copied to this limb (starting from 3992 the low end). */ 3993 size_t bytes; 3994 /* The index where the limb should be stored, when completed. */ 3995 mp_size_t i; 3996 3997 if (nails != 0) 3998 gmp_die ("mpz_import: Nails not supported."); 3999 4000 assert (order == 1 || order == -1); 4001 assert (endian >= -1 && endian <= 1); 4002 4003 if (endian == 0) 4004 endian = gmp_detect_endian (); 4005 4006 p = (unsigned char *) src; 4007 4008 word_step = (order != endian) ? 2 * size : 0; 4009 4010 /* Process bytes from the least significant end, so point p at the 4011 least significant word. */ 4012 if (order == 1) 4013 { 4014 p += size * (count - 1); 4015 word_step = - word_step; 4016 } 4017 4018 /* And at least significant byte of that word. */ 4019 if (endian == 1) 4020 p += (size - 1); 4021 4022 rn = (size * count + sizeof(mp_limb_t) - 1) / sizeof(mp_limb_t); 4023 rp = MPZ_REALLOC (r, rn); 4024 4025 for (limb = 0, bytes = 0, i = 0; count > 0; count--, p += word_step) 4026 { 4027 size_t j; 4028 for (j = 0; j < size; j++, p -= (ptrdiff_t) endian) 4029 { 4030 limb |= (mp_limb_t) *p << (bytes++ * CHAR_BIT); 4031 if (bytes == sizeof(mp_limb_t)) 4032 { 4033 rp[i++] = limb; 4034 bytes = 0; 4035 limb = 0; 4036 } 4037 } 4038 } 4039 if (bytes > 0) 4040 rp[i++] = limb; 4041 assert (i == rn); 4042 4043 r->_mp_size = mpn_normalized_size (rp, i); 4044 } 4045 4046 void * 4047 mpz_export (void *r, size_t *countp, int order, size_t size, int endian, 4048 size_t nails, const mpz_t u) 4049 { 4050 unsigned char *p; 4051 ptrdiff_t word_step; 4052 size_t count, k; 4053 mp_size_t un; 4054 4055 /* The current (partial) limb. */ 4056 mp_limb_t limb; 4057 /* The number of bytes left to to in this limb. */ 4058 size_t bytes; 4059 /* The index where the limb was read. */ 4060 mp_size_t i; 4061 4062 if (nails != 0) 4063 gmp_die ("mpz_import: Nails not supported."); 4064 4065 assert (order == 1 || order == -1); 4066 assert (endian >= -1 && endian <= 1); 4067 assert (size > 0 || u->_mp_size == 0); 4068 4069 un = GMP_ABS (u->_mp_size); 4070 if (un == 0) 4071 { 4072 if (countp) 4073 *countp = 0; 4074 return r; 4075 } 4076 4077 /* Count bytes in top limb. */ 4078 for (limb = u->_mp_d[un-1], k = 0; limb > 0; k++, limb >>= CHAR_BIT) 4079 ; 4080 4081 assert (k > 0); 4082 4083 count = (k + (un-1) * sizeof (mp_limb_t) + size - 1) / size; 4084 4085 if (!r) 4086 r = gmp_xalloc (count * size); 4087 4088 if (endian == 0) 4089 endian = gmp_detect_endian (); 4090 4091 p = (unsigned char *) r; 4092 4093 word_step = (order != endian) ? 2 * size : 0; 4094 4095 /* Process bytes from the least significant end, so point p at the 4096 least significant word. */ 4097 if (order == 1) 4098 { 4099 p += size * (count - 1); 4100 word_step = - word_step; 4101 } 4102 4103 /* And at least significant byte of that word. */ 4104 if (endian == 1) 4105 p += (size - 1); 4106 4107 for (bytes = 0, i = 0, k = 0; k < count; k++, p += word_step) 4108 { 4109 size_t j; 4110 for (j = 0; j < size; j++, p -= (ptrdiff_t) endian) 4111 { 4112 if (bytes == 0) 4113 { 4114 if (i < un) 4115 limb = u->_mp_d[i++]; 4116 bytes = sizeof (mp_limb_t); 4117 } 4118 *p = limb; 4119 limb >>= CHAR_BIT; 4120 bytes--; 4121 } 4122 } 4123 assert (i == un); 4124 assert (k == count); 4125 4126 if (countp) 4127 *countp = count; 4128 4129 return r; 4130 } 4131