dumbmp.c

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