xref: /dragonfly/contrib/gmp/dumbmp.c (revision 6e278935)
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 }
923