contrib/gmp/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 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 }