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