1 /* dumbmp mini GMP compatible library. 2 3 Copyright 2001, 2002, 2004 Free Software Foundation, Inc. 4 5 This file is part of the GNU MP Library. 6 7 The GNU MP Library is free software; you can redistribute it and/or modify 8 it under the terms of the GNU Lesser General Public License as published by 9 the Free Software Foundation; either version 3 of the License, or (at your 10 option) any later version. 11 12 The GNU MP Library is distributed in the hope that it will be useful, but 13 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 14 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 15 License for more details. 16 17 You should have received a copy of the GNU Lesser General Public License 18 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ 19 20 21 /* The code here implements a subset (a very limited subset) of the main GMP 22 functions. It's designed for use in a few build-time calculations and 23 will be slow, but highly portable. 24 25 None of the normal GMP configure things are used, nor any of the normal 26 gmp.h or gmp-impl.h. To use this file in a program just #include 27 "dumbmp.c". 28 29 ANSI function definitions can be used here, since ansi2knr is run if 30 necessary. But other ANSI-isms like "const" should be avoided. 31 32 mp_limb_t here is an unsigned long, since that's a sensible type 33 everywhere we know of, with 8*sizeof(unsigned long) giving the number of 34 bits in the type (that not being true for instance with int or short on 35 Cray vector systems.) 36 37 Only the low half of each mp_limb_t is used, so as to make carry handling 38 and limb multiplies easy. GMP_LIMB_BITS is the number of bits used. */ 39 40 #include <stdio.h> 41 #include <stdlib.h> 42 #include <string.h> 43 44 45 typedef unsigned long mp_limb_t; 46 47 typedef struct { 48 int _mp_alloc; 49 int _mp_size; 50 mp_limb_t *_mp_d; 51 } mpz_t[1]; 52 53 #define GMP_LIMB_BITS (sizeof (mp_limb_t) * 8 / 2) 54 55 #define ABS(x) ((x) >= 0 ? (x) : -(x)) 56 #define MIN(l,o) ((l) < (o) ? (l) : (o)) 57 #define MAX(h,i) ((h) > (i) ? (h) : (i)) 58 59 #define ALLOC(x) ((x)->_mp_alloc) 60 #define PTR(x) ((x)->_mp_d) 61 #define SIZ(x) ((x)->_mp_size) 62 #define ABSIZ(x) ABS (SIZ (x)) 63 #define LOMASK ((1L << GMP_LIMB_BITS) - 1) 64 #define LO(x) ((x) & LOMASK) 65 #define HI(x) ((x) >> GMP_LIMB_BITS) 66 67 #define ASSERT(cond) \ 68 do { \ 69 if (! (cond)) \ 70 { \ 71 fprintf (stderr, "Assertion failure\n"); \ 72 abort (); \ 73 } \ 74 } while (0) 75 76 77 char * 78 xmalloc (int n) 79 { 80 char *p; 81 p = malloc (n); 82 if (p == NULL) 83 { 84 fprintf (stderr, "Out of memory (alloc %d bytes)\n", n); 85 abort (); 86 } 87 return p; 88 } 89 90 mp_limb_t * 91 xmalloc_limbs (int n) 92 { 93 return (mp_limb_t *) xmalloc (n * sizeof (mp_limb_t)); 94 } 95 96 void 97 mem_copyi (char *dst, char *src, int size) 98 { 99 int i; 100 for (i = 0; i < size; i++) 101 dst[i] = src[i]; 102 } 103 104 int 105 isprime (int n) 106 { 107 int i; 108 if (n < 2) 109 return 0; 110 for (i = 2; i < n; i++) 111 if ((n % i) == 0) 112 return 0; 113 return 1; 114 } 115 116 int 117 log2_ceil (int n) 118 { 119 int e; 120 ASSERT (n >= 1); 121 for (e = 0; ; e++) 122 if ((1 << e) >= n) 123 break; 124 return e; 125 } 126 127 void 128 mpz_realloc (mpz_t r, int n) 129 { 130 if (n <= ALLOC(r)) 131 return; 132 133 ALLOC(r) = n; 134 PTR(r) = (mp_limb_t *) realloc (PTR(r), n * sizeof (mp_limb_t)); 135 if (PTR(r) == NULL) 136 { 137 fprintf (stderr, "Out of memory (realloc to %d)\n", n); 138 abort (); 139 } 140 } 141 142 void 143 mpn_normalize (mp_limb_t *rp, int *rnp) 144 { 145 int rn = *rnp; 146 while (rn > 0 && rp[rn-1] == 0) 147 rn--; 148 *rnp = rn; 149 } 150 151 void 152 mpn_copyi (mp_limb_t *dst, mp_limb_t *src, int n) 153 { 154 int i; 155 for (i = 0; i < n; i++) 156 dst[i] = src[i]; 157 } 158 159 void 160 mpn_zero (mp_limb_t *rp, int rn) 161 { 162 int i; 163 for (i = 0; i < rn; i++) 164 rp[i] = 0; 165 } 166 167 void 168 mpz_init (mpz_t r) 169 { 170 ALLOC(r) = 1; 171 PTR(r) = xmalloc_limbs (ALLOC(r)); 172 PTR(r)[0] = 0; 173 SIZ(r) = 0; 174 } 175 176 void 177 mpz_clear (mpz_t r) 178 { 179 free (PTR (r)); 180 ALLOC(r) = -1; 181 SIZ (r) = 0xbadcafeL; 182 PTR (r) = (mp_limb_t *) 0xdeadbeefL; 183 } 184 185 int 186 mpz_sgn (mpz_t a) 187 { 188 return (SIZ(a) > 0 ? 1 : SIZ(a) == 0 ? 0 : -1); 189 } 190 191 int 192 mpz_odd_p (mpz_t a) 193 { 194 if (SIZ(a) == 0) 195 return 0; 196 else 197 return (PTR(a)[0] & 1) != 0; 198 } 199 200 int 201 mpz_even_p (mpz_t a) 202 { 203 if (SIZ(a) == 0) 204 return 1; 205 else 206 return (PTR(a)[0] & 1) == 0; 207 } 208 209 size_t 210 mpz_sizeinbase (mpz_t a, int base) 211 { 212 int an = ABSIZ (a); 213 mp_limb_t *ap = PTR (a); 214 int cnt; 215 mp_limb_t hi; 216 217 if (base != 2) 218 abort (); 219 220 if (an == 0) 221 return 1; 222 223 cnt = 0; 224 for (hi = ap[an - 1]; hi != 0; hi >>= 1) 225 cnt += 1; 226 return (an - 1) * GMP_LIMB_BITS + cnt; 227 } 228 229 void 230 mpz_set (mpz_t r, mpz_t a) 231 { 232 mpz_realloc (r, ABSIZ (a)); 233 SIZ(r) = SIZ(a); 234 mpn_copyi (PTR(r), PTR(a), ABSIZ (a)); 235 } 236 237 void 238 mpz_init_set (mpz_t r, mpz_t a) 239 { 240 mpz_init (r); 241 mpz_set (r, a); 242 } 243 244 void 245 mpz_set_ui (mpz_t r, unsigned long ui) 246 { 247 int rn; 248 mpz_realloc (r, 2); 249 PTR(r)[0] = LO(ui); 250 PTR(r)[1] = HI(ui); 251 rn = 2; 252 mpn_normalize (PTR(r), &rn); 253 SIZ(r) = rn; 254 } 255 256 void 257 mpz_init_set_ui (mpz_t r, unsigned long ui) 258 { 259 mpz_init (r); 260 mpz_set_ui (r, ui); 261 } 262 263 void 264 mpz_setbit (mpz_t r, unsigned long bit) 265 { 266 int limb, rn, extend; 267 mp_limb_t *rp; 268 269 rn = SIZ(r); 270 if (rn < 0) 271 abort (); /* only r>=0 */ 272 273 limb = bit / GMP_LIMB_BITS; 274 bit %= GMP_LIMB_BITS; 275 276 mpz_realloc (r, limb+1); 277 rp = PTR(r); 278 extend = (limb+1) - rn; 279 if (extend > 0) 280 mpn_zero (rp + rn, extend); 281 282 rp[limb] |= (mp_limb_t) 1 << bit; 283 SIZ(r) = MAX (rn, limb+1); 284 } 285 286 int 287 mpz_tstbit (mpz_t r, unsigned long bit) 288 { 289 int limb; 290 291 if (SIZ(r) < 0) 292 abort (); /* only r>=0 */ 293 294 limb = bit / GMP_LIMB_BITS; 295 if (SIZ(r) <= limb) 296 return 0; 297 298 bit %= GMP_LIMB_BITS; 299 return (PTR(r)[limb] >> bit) & 1; 300 } 301 302 int 303 popc_limb (mp_limb_t a) 304 { 305 int ret = 0; 306 while (a != 0) 307 { 308 ret += (a & 1); 309 a >>= 1; 310 } 311 return ret; 312 } 313 314 unsigned long 315 mpz_popcount (mpz_t a) 316 { 317 unsigned long ret; 318 int i; 319 320 if (SIZ(a) < 0) 321 abort (); 322 323 ret = 0; 324 for (i = 0; i < SIZ(a); i++) 325 ret += popc_limb (PTR(a)[i]); 326 return ret; 327 } 328 329 void 330 mpz_add (mpz_t r, mpz_t a, mpz_t b) 331 { 332 int an = ABSIZ (a), bn = ABSIZ (b), rn; 333 mp_limb_t *rp, *ap, *bp; 334 int i; 335 mp_limb_t t, cy; 336 337 if ((SIZ (a) ^ SIZ (b)) < 0) 338 abort (); /* really subtraction */ 339 if (SIZ (a) < 0) 340 abort (); 341 342 mpz_realloc (r, MAX (an, bn) + 1); 343 ap = PTR (a); bp = PTR (b); rp = PTR (r); 344 if (an < bn) 345 { 346 mp_limb_t *tp; int tn; 347 tn = an; an = bn; bn = tn; 348 tp = ap; ap = bp; bp = tp; 349 } 350 351 cy = 0; 352 for (i = 0; i < bn; i++) 353 { 354 t = ap[i] + bp[i] + cy; 355 rp[i] = LO (t); 356 cy = HI (t); 357 } 358 for (i = bn; i < an; i++) 359 { 360 t = ap[i] + cy; 361 rp[i] = LO (t); 362 cy = HI (t); 363 } 364 rp[an] = cy; 365 rn = an + 1; 366 367 mpn_normalize (rp, &rn); 368 SIZ (r) = rn; 369 } 370 371 void 372 mpz_add_ui (mpz_t r, mpz_t a, unsigned long int ui) 373 { 374 mpz_t b; 375 376 mpz_init (b); 377 mpz_set_ui (b, ui); 378 mpz_add (r, a, b); 379 mpz_clear (b); 380 } 381 382 void 383 mpz_sub (mpz_t r, mpz_t a, mpz_t b) 384 { 385 int an = ABSIZ (a), bn = ABSIZ (b), rn; 386 mp_limb_t *rp, *ap, *bp; 387 int i; 388 mp_limb_t t, cy; 389 390 if ((SIZ (a) ^ SIZ (b)) < 0) 391 abort (); /* really addition */ 392 if (SIZ (a) < 0) 393 abort (); 394 395 mpz_realloc (r, MAX (an, bn) + 1); 396 ap = PTR (a); bp = PTR (b); rp = PTR (r); 397 if (an < bn) 398 { 399 mp_limb_t *tp; int tn; 400 tn = an; an = bn; bn = tn; 401 tp = ap; ap = bp; bp = tp; 402 } 403 404 cy = 0; 405 for (i = 0; i < bn; i++) 406 { 407 t = ap[i] - bp[i] - cy; 408 rp[i] = LO (t); 409 cy = LO (-HI (t)); 410 } 411 for (i = bn; i < an; i++) 412 { 413 t = ap[i] - cy; 414 rp[i] = LO (t); 415 cy = LO (-HI (t)); 416 } 417 rp[an] = cy; 418 rn = an + 1; 419 420 if (cy != 0) 421 { 422 cy = 0; 423 for (i = 0; i < rn; i++) 424 { 425 t = -rp[i] - cy; 426 rp[i] = LO (t); 427 cy = LO (-HI (t)); 428 } 429 SIZ (r) = -rn; 430 return; 431 } 432 433 mpn_normalize (rp, &rn); 434 SIZ (r) = rn; 435 } 436 437 void 438 mpz_sub_ui (mpz_t r, mpz_t a, unsigned long int ui) 439 { 440 mpz_t b; 441 442 mpz_init (b); 443 mpz_set_ui (b, ui); 444 mpz_sub (r, a, b); 445 mpz_clear (b); 446 } 447 448 void 449 mpz_mul (mpz_t r, mpz_t a, mpz_t b) 450 { 451 int an = ABSIZ (a), bn = ABSIZ (b), rn; 452 mp_limb_t *scratch, *tmp, *ap = PTR (a), *bp = PTR (b); 453 int ai, bi; 454 mp_limb_t t, cy; 455 456 scratch = xmalloc_limbs (an + bn); 457 tmp = scratch; 458 459 for (bi = 0; bi < bn; bi++) 460 tmp[bi] = 0; 461 462 for (ai = 0; ai < an; ai++) 463 { 464 tmp = scratch + ai; 465 cy = 0; 466 for (bi = 0; bi < bn; bi++) 467 { 468 t = ap[ai] * bp[bi] + tmp[bi] + cy; 469 tmp[bi] = LO (t); 470 cy = HI (t); 471 } 472 tmp[bn] = cy; 473 } 474 475 rn = an + bn; 476 mpn_normalize (scratch, &rn); 477 free (PTR (r)); 478 PTR (r) = scratch; 479 SIZ (r) = (SIZ (a) ^ SIZ (b)) >= 0 ? rn : -rn; 480 ALLOC (r) = an + bn; 481 } 482 483 void 484 mpz_mul_ui (mpz_t r, mpz_t a, unsigned long int ui) 485 { 486 mpz_t b; 487 488 mpz_init (b); 489 mpz_set_ui (b, ui); 490 mpz_mul (r, a, b); 491 mpz_clear (b); 492 } 493 494 void 495 mpz_mul_2exp (mpz_t r, mpz_t a, unsigned long int bcnt) 496 { 497 mpz_set (r, a); 498 while (bcnt) 499 { 500 mpz_add (r, r, r); 501 bcnt -= 1; 502 } 503 } 504 505 void 506 mpz_ui_pow_ui (mpz_t r, unsigned long b, unsigned long e) 507 { 508 unsigned long i; 509 mpz_t bz; 510 511 mpz_init (bz); 512 mpz_set_ui (bz, b); 513 514 mpz_set_ui (r, 1L); 515 for (i = 0; i < e; i++) 516 mpz_mul (r, r, bz); 517 518 mpz_clear (bz); 519 } 520 521 void 522 mpz_tdiv_q_2exp (mpz_t r, mpz_t a, unsigned long int bcnt) 523 { 524 int as, rn; 525 int cnt, tnc; 526 int lcnt; 527 mp_limb_t high_limb, low_limb; 528 int i; 529 530 as = SIZ (a); 531 lcnt = bcnt / GMP_LIMB_BITS; 532 rn = ABS (as) - lcnt; 533 if (rn <= 0) 534 SIZ (r) = 0; 535 else 536 { 537 mp_limb_t *rp, *ap; 538 539 mpz_realloc (r, rn); 540 541 rp = PTR (r); 542 ap = PTR (a); 543 544 cnt = bcnt % GMP_LIMB_BITS; 545 if (cnt != 0) 546 { 547 ap += lcnt; 548 tnc = GMP_LIMB_BITS - cnt; 549 high_limb = *ap++; 550 low_limb = high_limb >> cnt; 551 552 for (i = rn - 1; i != 0; i--) 553 { 554 high_limb = *ap++; 555 *rp++ = low_limb | LO (high_limb << tnc); 556 low_limb = high_limb >> cnt; 557 } 558 *rp = low_limb; 559 rn -= low_limb == 0; 560 } 561 else 562 { 563 ap += lcnt; 564 mpn_copyi (rp, ap, rn); 565 } 566 567 SIZ (r) = as >= 0 ? rn : -rn; 568 } 569 } 570 571 void 572 mpz_tdiv_r_2exp (mpz_t r, mpz_t a, unsigned long int bcnt) 573 { 574 int rn, bwhole; 575 576 mpz_set (r, a); 577 rn = ABSIZ(r); 578 579 bwhole = bcnt / GMP_LIMB_BITS; 580 bcnt %= GMP_LIMB_BITS; 581 if (rn > bwhole) 582 { 583 rn = bwhole+1; 584 PTR(r)[rn-1] &= ((mp_limb_t) 1 << bcnt) - 1; 585 mpn_normalize (PTR(r), &rn); 586 SIZ(r) = (SIZ(r) >= 0 ? rn : -rn); 587 } 588 } 589 590 int 591 mpz_cmp (mpz_t a, mpz_t b) 592 { 593 mp_limb_t *ap, *bp, al, bl; 594 int as = SIZ (a), bs = SIZ (b); 595 int i; 596 int sign; 597 598 if (as != bs) 599 return as > bs ? 1 : -1; 600 601 sign = as > 0 ? 1 : -1; 602 603 ap = PTR (a); 604 bp = PTR (b); 605 for (i = ABS (as) - 1; i >= 0; i--) 606 { 607 al = ap[i]; 608 bl = bp[i]; 609 if (al != bl) 610 return al > bl ? sign : -sign; 611 } 612 return 0; 613 } 614 615 int 616 mpz_cmp_ui (mpz_t a, unsigned long b) 617 { 618 mpz_t bz; 619 int ret; 620 mpz_init_set_ui (bz, b); 621 ret = mpz_cmp (a, bz); 622 mpz_clear (bz); 623 return ret; 624 } 625 626 void 627 mpz_tdiv_qr (mpz_t q, mpz_t r, mpz_t a, mpz_t b) 628 { 629 mpz_t tmpr, tmpb; 630 unsigned long cnt; 631 632 ASSERT (mpz_sgn(a) >= 0); 633 ASSERT (mpz_sgn(b) > 0); 634 635 mpz_init_set (tmpr, a); 636 mpz_init_set (tmpb, b); 637 mpz_set_ui (q, 0L); 638 639 if (mpz_cmp (tmpr, tmpb) > 0) 640 { 641 cnt = mpz_sizeinbase (tmpr, 2) - mpz_sizeinbase (tmpb, 2) + 1; 642 mpz_mul_2exp (tmpb, tmpb, cnt); 643 644 for ( ; cnt > 0; cnt--) 645 { 646 mpz_mul_2exp (q, q, 1); 647 mpz_tdiv_q_2exp (tmpb, tmpb, 1L); 648 if (mpz_cmp (tmpr, tmpb) >= 0) 649 { 650 mpz_sub (tmpr, tmpr, tmpb); 651 mpz_add_ui (q, q, 1L); 652 ASSERT (mpz_cmp (tmpr, tmpb) < 0); 653 } 654 } 655 } 656 657 mpz_set (r, tmpr); 658 mpz_clear (tmpr); 659 mpz_clear (tmpb); 660 } 661 662 void 663 mpz_tdiv_qr_ui (mpz_t q, mpz_t r, mpz_t a, unsigned long b) 664 { 665 mpz_t bz; 666 mpz_init_set_ui (bz, b); 667 mpz_tdiv_qr (q, r, a, bz); 668 mpz_clear (bz); 669 } 670 671 void 672 mpz_tdiv_q (mpz_t q, mpz_t a, mpz_t b) 673 { 674 mpz_t r; 675 676 mpz_init (r); 677 mpz_tdiv_qr (q, r, a, b); 678 mpz_clear (r); 679 } 680 681 void 682 mpz_tdiv_q_ui (mpz_t q, mpz_t n, unsigned long d) 683 { 684 mpz_t dz; 685 mpz_init_set_ui (dz, d); 686 mpz_tdiv_q (q, n, dz); 687 mpz_clear (dz); 688 } 689 690 /* Set inv to the inverse of d, in the style of invert_limb, ie. for 691 udiv_qrnnd_preinv. */ 692 void 693 mpz_preinv_invert (mpz_t inv, mpz_t d, int numb_bits) 694 { 695 mpz_t t; 696 int norm; 697 ASSERT (SIZ(d) > 0); 698 699 norm = numb_bits - mpz_sizeinbase (d, 2); 700 ASSERT (norm >= 0); 701 mpz_init_set_ui (t, 1L); 702 mpz_mul_2exp (t, t, 2*numb_bits - norm); 703 mpz_tdiv_q (inv, t, d); 704 mpz_set_ui (t, 1L); 705 mpz_mul_2exp (t, t, numb_bits); 706 mpz_sub (inv, inv, t); 707 708 mpz_clear (t); 709 } 710 711 /* Remove leading '0' characters from the start of a string, by copying the 712 remainder down. */ 713 void 714 strstrip_leading_zeros (char *s) 715 { 716 char c, *p; 717 718 p = s; 719 while (*s == '0') 720 s++; 721 722 do 723 { 724 c = *s++; 725 *p++ = c; 726 } 727 while (c != '\0'); 728 } 729 730 char * 731 mpz_get_str (char *buf, int base, mpz_t a) 732 { 733 static char tohex[] = "0123456789abcdef"; 734 735 mp_limb_t alimb, *ap; 736 int an, bn, i, j; 737 char *bp; 738 739 if (base != 16) 740 abort (); 741 if (SIZ (a) < 0) 742 abort (); 743 744 if (buf == 0) 745 buf = xmalloc (ABSIZ (a) * (GMP_LIMB_BITS / 4) + 3); 746 747 an = ABSIZ (a); 748 if (an == 0) 749 { 750 buf[0] = '0'; 751 buf[1] = '\0'; 752 return buf; 753 } 754 755 ap = PTR (a); 756 bn = an * (GMP_LIMB_BITS / 4); 757 bp = buf + bn; 758 759 for (i = 0; i < an; i++) 760 { 761 alimb = ap[i]; 762 for (j = 0; j < GMP_LIMB_BITS / 4; j++) 763 { 764 bp--; 765 *bp = tohex [alimb & 0xF]; 766 alimb >>= 4; 767 } 768 ASSERT (alimb == 0); 769 } 770 ASSERT (bp == buf); 771 772 buf[bn] = '\0'; 773 774 strstrip_leading_zeros (buf); 775 return buf; 776 } 777 778 void 779 mpz_out_str (FILE *file, int base, mpz_t a) 780 { 781 char *str; 782 783 if (file == 0) 784 file = stdout; 785 786 str = mpz_get_str (0, 16, a); 787 fputs (str, file); 788 free (str); 789 } 790 791 /* Calculate r satisfying r*d == 1 mod 2^n. */ 792 void 793 mpz_invert_2exp (mpz_t r, mpz_t a, unsigned long n) 794 { 795 unsigned long i; 796 mpz_t inv, prod; 797 798 ASSERT (mpz_odd_p (a)); 799 800 mpz_init_set_ui (inv, 1L); 801 mpz_init (prod); 802 803 for (i = 1; i < n; i++) 804 { 805 mpz_mul (prod, inv, a); 806 if (mpz_tstbit (prod, i) != 0) 807 mpz_setbit (inv, i); 808 } 809 810 mpz_mul (prod, inv, a); 811 mpz_tdiv_r_2exp (prod, prod, n); 812 ASSERT (mpz_cmp_ui (prod, 1L) == 0); 813 814 mpz_set (r, inv); 815 816 mpz_clear (inv); 817 mpz_clear (prod); 818 } 819 820 /* Calculate inv satisfying r*a == 1 mod 2^n. */ 821 void 822 mpz_invert_ui_2exp (mpz_t r, unsigned long a, unsigned long n) 823 { 824 mpz_t az; 825 mpz_init_set_ui (az, a); 826 mpz_invert_2exp (r, az, n); 827 mpz_clear (az); 828 } 829 830 /* x=y^z */ 831 void 832 mpz_pow_ui (mpz_t x, mpz_t y, unsigned long z) 833 { 834 mpz_t t; 835 836 mpz_init_set_ui (t, 1); 837 for (; z != 0; z--) 838 mpz_mul (t, t, y); 839 mpz_set (x, t); 840 mpz_clear (t); 841 } 842 843 /* x=x+y*z */ 844 void 845 mpz_addmul_ui (mpz_t x, mpz_t y, unsigned long z) 846 { 847 mpz_t t; 848 849 mpz_init (t); 850 mpz_mul_ui (t, y, z); 851 mpz_add (x, x, t); 852 mpz_clear (t); 853 } 854 855 /* x=floor(y^(1/z)) */ 856 void 857 mpz_root (mpz_t x, mpz_t y, unsigned long z) 858 { 859 mpz_t t, u; 860 861 if (mpz_sgn (y) < 0) 862 { 863 fprintf (stderr, "mpz_root does not accept negative values\n"); 864 abort (); 865 } 866 if (mpz_cmp_ui (y, 1) <= 0) 867 { 868 mpz_set (x, y); 869 return; 870 } 871 mpz_init (t); 872 mpz_init_set (u, y); 873 do 874 { 875 mpz_pow_ui (t, u, z - 1); 876 mpz_tdiv_q (t, y, t); 877 mpz_addmul_ui (t, u, z - 1); 878 mpz_tdiv_q_ui (t, t, z); 879 if (mpz_cmp (t, u) >= 0) 880 break; 881 mpz_set (u, t); 882 } 883 while (1); 884 mpz_set (x, u); 885 mpz_clear (t); 886 mpz_clear (u); 887 } 888