@@ -0,0 +1,4120 @@
1+/* mini-gmp, a minimalistic implementation of a GNU GMP subset.
2+
3+ Contributed to the GNU project by Niels Möller
4+
5+Copyright 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1999, 2000, 2001,
6+2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013
7+Free Software Foundation, Inc.
8+
9+This file is part of the GNU MP Library.
10+
11+The GNU MP Library is free software; you can redistribute it and/or modify
12+it under the terms of the GNU Lesser General Public License as published by
13+the Free Software Foundation; either version 3 of the License, or (at your
14+option) any later version.
15+
16+The GNU MP Library is distributed in the hope that it will be useful, but
17+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18+or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
19+License for more details.
20+
21+You should have received a copy of the GNU Lesser General Public License
22+along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
23+
24+/* NOTE: All functions in this file which are not declared in
25+ mini-gmp.h are internal, and are not intended to be compatible
26+ neither with GMP nor with future versions of mini-gmp. */
27+
28+/* Much of the material copied from GMP files, including: gmp-impl.h,
29+ longlong.h, mpn/generic/add_n.c, mpn/generic/addmul_1.c,
30+ mpn/generic/lshift.c, mpn/generic/mul_1.c,
31+ mpn/generic/mul_basecase.c, mpn/generic/rshift.c,
32+ mpn/generic/sbpi1_div_qr.c, mpn/generic/sub_n.c,
33+ mpn/generic/submul_1.c. */
34+
35+#include <assert.h>
36+#include <ctype.h>
37+#include <limits.h>
38+#include <stdio.h>
39+#include <stdlib.h>
40+#include <string.h>
41+
42+#include "mini-gmp.h"
43+
44+
45+/* Macros */
46+#define GMP_LIMB_BITS (sizeof(mp_limb_t) * CHAR_BIT)
47+
48+#define GMP_LIMB_MAX (~ (mp_limb_t) 0)
49+#define GMP_LIMB_HIGHBIT ((mp_limb_t) 1 << (GMP_LIMB_BITS - 1))
50+
51+#define GMP_HLIMB_BIT ((mp_limb_t) 1 << (GMP_LIMB_BITS / 2))
52+#define GMP_LLIMB_MASK (GMP_HLIMB_BIT - 1)
53+
54+#define GMP_ULONG_BITS (sizeof(unsigned long) * CHAR_BIT)
55+#define GMP_ULONG_HIGHBIT ((unsigned long) 1 << (GMP_ULONG_BITS - 1))
56+
57+#define GMP_ABS(x) ((x) >= 0 ? (x) : -(x))
58+#define GMP_NEG_CAST(T,x) (-((T)((x) + 1) - 1))
59+
60+#define GMP_MIN(a, b) ((a) < (b) ? (a) : (b))
61+#define GMP_MAX(a, b) ((a) > (b) ? (a) : (b))
62+
63+#define gmp_assert_nocarry(x) do { \
64+ mp_limb_t __cy = x; \
65+ assert (__cy == 0); \
66+ } while (0)
67+
68+#define gmp_clz(count, x) do { \
69+ mp_limb_t __clz_x = (x); \
70+ unsigned __clz_c; \
71+ for (__clz_c = 0; \
72+ (__clz_x & ((mp_limb_t) 0xff << (GMP_LIMB_BITS - 8))) == 0; \
73+ __clz_c += 8) \
74+ __clz_x <<= 8; \
75+ for (; (__clz_x & GMP_LIMB_HIGHBIT) == 0; __clz_c++) \
76+ __clz_x <<= 1; \
77+ (count) = __clz_c; \
78+ } while (0)
79+
80+#define gmp_ctz(count, x) do { \
81+ mp_limb_t __ctz_x = (x); \
82+ unsigned __ctz_c = 0; \
83+ gmp_clz (__ctz_c, __ctz_x & - __ctz_x); \
84+ (count) = GMP_LIMB_BITS - 1 - __ctz_c; \
85+ } while (0)
86+
87+#define gmp_add_ssaaaa(sh, sl, ah, al, bh, bl) \
88+ do { \
89+ mp_limb_t __x; \
90+ __x = (al) + (bl); \
91+ (sh) = (ah) + (bh) + (__x < (al)); \
92+ (sl) = __x; \
93+ } while (0)
94+
95+#define gmp_sub_ddmmss(sh, sl, ah, al, bh, bl) \
96+ do { \
97+ mp_limb_t __x; \
98+ __x = (al) - (bl); \
99+ (sh) = (ah) - (bh) - ((al) < (bl)); \
100+ (sl) = __x; \
101+ } while (0)
102+
103+#define gmp_umul_ppmm(w1, w0, u, v) \
104+ do { \
105+ mp_limb_t __x0, __x1, __x2, __x3; \
106+ unsigned __ul, __vl, __uh, __vh; \
107+ mp_limb_t __u = (u), __v = (v); \
108+ \
109+ __ul = __u & GMP_LLIMB_MASK; \
110+ __uh = __u >> (GMP_LIMB_BITS / 2); \
111+ __vl = __v & GMP_LLIMB_MASK; \
112+ __vh = __v >> (GMP_LIMB_BITS / 2); \
113+ \
114+ __x0 = (mp_limb_t) __ul * __vl; \
115+ __x1 = (mp_limb_t) __ul * __vh; \
116+ __x2 = (mp_limb_t) __uh * __vl; \
117+ __x3 = (mp_limb_t) __uh * __vh; \
118+ \
119+ __x1 += __x0 >> (GMP_LIMB_BITS / 2);/* this can't give carry */ \
120+ __x1 += __x2; /* but this indeed can */ \
121+ if (__x1 < __x2) /* did we get it? */ \
122+ __x3 += GMP_HLIMB_BIT; /* yes, add it in the proper pos. */ \
123+ \
124+ (w1) = __x3 + (__x1 >> (GMP_LIMB_BITS / 2)); \
125+ (w0) = (__x1 << (GMP_LIMB_BITS / 2)) + (__x0 & GMP_LLIMB_MASK); \
126+ } while (0)
127+
128+#define gmp_udiv_qrnnd_preinv(q, r, nh, nl, d, di) \
129+ do { \
130+ mp_limb_t _qh, _ql, _r, _mask; \
131+ gmp_umul_ppmm (_qh, _ql, (nh), (di)); \
132+ gmp_add_ssaaaa (_qh, _ql, _qh, _ql, (nh) + 1, (nl)); \
133+ _r = (nl) - _qh * (d); \
134+ _mask = -(mp_limb_t) (_r > _ql); /* both > and >= are OK */ \
135+ _qh += _mask; \
136+ _r += _mask & (d); \
137+ if (_r >= (d)) \
138+ { \
139+ _r -= (d); \
140+ _qh++; \
141+ } \
142+ \
143+ (r) = _r; \
144+ (q) = _qh; \
145+ } while (0)
146+
147+#define gmp_udiv_qr_3by2(q, r1, r0, n2, n1, n0, d1, d0, dinv) \
148+ do { \
149+ mp_limb_t _q0, _t1, _t0, _mask; \
150+ gmp_umul_ppmm ((q), _q0, (n2), (dinv)); \
151+ gmp_add_ssaaaa ((q), _q0, (q), _q0, (n2), (n1)); \
152+ \
153+ /* Compute the two most significant limbs of n - q'd */ \
154+ (r1) = (n1) - (d1) * (q); \
155+ gmp_sub_ddmmss ((r1), (r0), (r1), (n0), (d1), (d0)); \
156+ gmp_umul_ppmm (_t1, _t0, (d0), (q)); \
157+ gmp_sub_ddmmss ((r1), (r0), (r1), (r0), _t1, _t0); \
158+ (q)++; \
159+ \
160+ /* Conditionally adjust q and the remainders */ \
161+ _mask = - (mp_limb_t) ((r1) >= _q0); \
162+ (q) += _mask; \
163+ gmp_add_ssaaaa ((r1), (r0), (r1), (r0), _mask & (d1), _mask & (d0)); \
164+ if ((r1) >= (d1)) \
165+ { \
166+ if ((r1) > (d1) || (r0) >= (d0)) \
167+ { \
168+ (q)++; \
169+ gmp_sub_ddmmss ((r1), (r0), (r1), (r0), (d1), (d0)); \
170+ } \
171+ } \
172+ } while (0)
173+
174+/* Swap macros. */
175+#define MP_LIMB_T_SWAP(x, y) \
176+ do { \
177+ mp_limb_t __mp_limb_t_swap__tmp = (x); \
178+ (x) = (y); \
179+ (y) = __mp_limb_t_swap__tmp; \
180+ } while (0)
181+#define MP_SIZE_T_SWAP(x, y) \
182+ do { \
183+ mp_size_t __mp_size_t_swap__tmp = (x); \
184+ (x) = (y); \
185+ (y) = __mp_size_t_swap__tmp; \
186+ } while (0)
187+#define MP_BITCNT_T_SWAP(x,y) \
188+ do { \
189+ mp_bitcnt_t __mp_bitcnt_t_swap__tmp = (x); \
190+ (x) = (y); \
191+ (y) = __mp_bitcnt_t_swap__tmp; \
192+ } while (0)
193+#define MP_PTR_SWAP(x, y) \
194+ do { \
195+ mp_ptr __mp_ptr_swap__tmp = (x); \
196+ (x) = (y); \
197+ (y) = __mp_ptr_swap__tmp; \
198+ } while (0)
199+#define MP_SRCPTR_SWAP(x, y) \
200+ do { \
201+ mp_srcptr __mp_srcptr_swap__tmp = (x); \
202+ (x) = (y); \
203+ (y) = __mp_srcptr_swap__tmp; \
204+ } while (0)
205+
206+#define MPN_PTR_SWAP(xp,xs, yp,ys) \
207+ do { \
208+ MP_PTR_SWAP (xp, yp); \
209+ MP_SIZE_T_SWAP (xs, ys); \
210+ } while(0)
211+#define MPN_SRCPTR_SWAP(xp,xs, yp,ys) \
212+ do { \
213+ MP_SRCPTR_SWAP (xp, yp); \
214+ MP_SIZE_T_SWAP (xs, ys); \
215+ } while(0)
216+
217+#define MPZ_PTR_SWAP(x, y) \
218+ do { \
219+ mpz_ptr __mpz_ptr_swap__tmp = (x); \
220+ (x) = (y); \
221+ (y) = __mpz_ptr_swap__tmp; \
222+ } while (0)
223+#define MPZ_SRCPTR_SWAP(x, y) \
224+ do { \
225+ mpz_srcptr __mpz_srcptr_swap__tmp = (x); \
226+ (x) = (y); \
227+ (y) = __mpz_srcptr_swap__tmp; \
228+ } while (0)
229+
230+
231+/* Memory allocation and other helper functions. */
232+static void
233+gmp_die (const char *msg)
234+{
235+ fprintf (stderr, "%s\n", msg);
236+ abort();
237+}
238+
239+static void *
240+gmp_default_alloc (size_t size)
241+{
242+ void *p;
243+
244+ assert (size > 0);
245+
246+ p = malloc (size);
247+ if (!p)
248+ gmp_die("gmp_default_alloc: Virtual memory exhausted.");
249+
250+ return p;
251+}
252+
253+static void *
254+gmp_default_realloc (void *old, size_t old_size, size_t new_size)
255+{
256+ mp_ptr p;
257+
258+ p = realloc (old, new_size);
259+
260+ if (!p)
261+ gmp_die("gmp_default_realoc: Virtual memory exhausted.");
262+
263+ return p;
264+}
265+
266+static void
267+gmp_default_free (void *p, size_t size)
268+{
269+ free (p);
270+}
271+
272+static void * (*gmp_allocate_func) (size_t) = gmp_default_alloc;
273+static void * (*gmp_reallocate_func) (void *, size_t, size_t) = gmp_default_realloc;
274+static void (*gmp_free_func) (void *, size_t) = gmp_default_free;
275+
276+void
277+mp_get_memory_functions (void *(**alloc_func) (size_t),
278+ void *(**realloc_func) (void *, size_t, size_t),
279+ void (**free_func) (void *, size_t))
280+{
281+ if (alloc_func)
282+ *alloc_func = gmp_allocate_func;
283+
284+ if (realloc_func)
285+ *realloc_func = gmp_reallocate_func;
286+
287+ if (free_func)
288+ *free_func = gmp_free_func;
289+}
290+
291+void
292+mp_set_memory_functions (void *(*alloc_func) (size_t),
293+ void *(*realloc_func) (void *, size_t, size_t),
294+ void (*free_func) (void *, size_t))
295+{
296+ if (!alloc_func)
297+ alloc_func = gmp_default_alloc;
298+ if (!realloc_func)
299+ realloc_func = gmp_default_realloc;
300+ if (!free_func)
301+ free_func = gmp_default_free;
302+
303+ gmp_allocate_func = alloc_func;
304+ gmp_reallocate_func = realloc_func;
305+ gmp_free_func = free_func;
306+}
307+
308+#define gmp_xalloc(size) ((*gmp_allocate_func)((size)))
309+#define gmp_free(p) ((*gmp_free_func) ((p), 0))
310+
311+static mp_ptr
312+gmp_xalloc_limbs (mp_size_t size)
313+{
314+ return gmp_xalloc (size * sizeof (mp_limb_t));
315+}
316+
317+static mp_ptr
318+gmp_xrealloc_limbs (mp_ptr old, mp_size_t size)
319+{
320+ assert (size > 0);
321+ return (*gmp_reallocate_func) (old, 0, size * sizeof (mp_limb_t));
322+}
323+
324+
325+/* MPN interface */
326+
327+void
328+mpn_copyi (mp_ptr d, mp_srcptr s, mp_size_t n)
329+{
330+ mp_size_t i;
331+ for (i = 0; i < n; i++)
332+ d[i] = s[i];
333+}
334+
335+void
336+mpn_copyd (mp_ptr d, mp_srcptr s, mp_size_t n)
337+{
338+ while (n-- > 0)
339+ d[n] = s[n];
340+}
341+
342+int
343+mpn_cmp (mp_srcptr ap, mp_srcptr bp, mp_size_t n)
344+{
345+ while (--n >= 0)
346+ {
347+ if (ap[n] != bp[n])
348+ return ap[n] > bp[n] ? 1 : -1;
349+ }
350+ return 0;
351+}
352+
353+static int
354+mpn_cmp4 (mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn)
355+{
356+ if (an != bn)
357+ return an < bn ? -1 : 1;
358+ else
359+ return mpn_cmp (ap, bp, an);
360+}
361+
362+static mp_size_t
363+mpn_normalized_size (mp_srcptr xp, mp_size_t n)
364+{
365+ for (; n > 0 && xp[n-1] == 0; n--)
366+ ;
367+ return n;
368+}
369+
370+#define mpn_zero_p(xp, n) (mpn_normalized_size ((xp), (n)) == 0)
371+
372+mp_limb_t
373+mpn_add_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b)
374+{
375+ mp_size_t i;
376+
377+ assert (n > 0);
378+ i = 0;
379+ do
380+ {
381+ mp_limb_t r = ap[i] + b;
382+ /* Carry out */
383+ b = (r < b);
384+ rp[i] = r;
385+ }
386+ while (++i < n);
387+
388+ return b;
389+}
390+
391+mp_limb_t
392+mpn_add_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n)
393+{
394+ mp_size_t i;
395+ mp_limb_t cy;
396+
397+ for (i = 0, cy = 0; i < n; i++)
398+ {
399+ mp_limb_t a, b, r;
400+ a = ap[i]; b = bp[i];
401+ r = a + cy;
402+ cy = (r < cy);
403+ r += b;
404+ cy += (r < b);
405+ rp[i] = r;
406+ }
407+ return cy;
408+}
409+
410+mp_limb_t
411+mpn_add (mp_ptr rp, mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn)
412+{
413+ mp_limb_t cy;
414+
415+ assert (an >= bn);
416+
417+ cy = mpn_add_n (rp, ap, bp, bn);
418+ if (an > bn)
419+ cy = mpn_add_1 (rp + bn, ap + bn, an - bn, cy);
420+ return cy;
421+}
422+
423+mp_limb_t
424+mpn_sub_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b)
425+{
426+ mp_size_t i;
427+
428+ assert (n > 0);
429+
430+ i = 0;
431+ do
432+ {
433+ mp_limb_t a = ap[i];
434+ /* Carry out */
435+ mp_limb_t cy = a < b;;
436+ rp[i] = a - b;
437+ b = cy;
438+ }
439+ while (++i < n);
440+
441+ return b;
442+}
443+
444+mp_limb_t
445+mpn_sub_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n)
446+{
447+ mp_size_t i;
448+ mp_limb_t cy;
449+
450+ for (i = 0, cy = 0; i < n; i++)
451+ {
452+ mp_limb_t a, b;
453+ a = ap[i]; b = bp[i];
454+ b += cy;
455+ cy = (b < cy);
456+ cy += (a < b);
457+ rp[i] = a - b;
458+ }
459+ return cy;
460+}
461+
462+mp_limb_t
463+mpn_sub (mp_ptr rp, mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn)
464+{
465+ mp_limb_t cy;
466+
467+ assert (an >= bn);
468+
469+ cy = mpn_sub_n (rp, ap, bp, bn);
470+ if (an > bn)
471+ cy = mpn_sub_1 (rp + bn, ap + bn, an - bn, cy);
472+ return cy;
473+}
474+
475+mp_limb_t
476+mpn_mul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl)
477+{
478+ mp_limb_t ul, cl, hpl, lpl;
479+
480+ assert (n >= 1);
481+
482+ cl = 0;
483+ do
484+ {
485+ ul = *up++;
486+ gmp_umul_ppmm (hpl, lpl, ul, vl);
487+
488+ lpl += cl;
489+ cl = (lpl < cl) + hpl;
490+
491+ *rp++ = lpl;
492+ }
493+ while (--n != 0);
494+
495+ return cl;
496+}
497+
498+mp_limb_t
499+mpn_addmul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl)
500+{
501+ mp_limb_t ul, cl, hpl, lpl, rl;
502+
503+ assert (n >= 1);
504+
505+ cl = 0;
506+ do
507+ {
508+ ul = *up++;
509+ gmp_umul_ppmm (hpl, lpl, ul, vl);
510+
511+ lpl += cl;
512+ cl = (lpl < cl) + hpl;
513+
514+ rl = *rp;
515+ lpl = rl + lpl;
516+ cl += lpl < rl;
517+ *rp++ = lpl;
518+ }
519+ while (--n != 0);
520+
521+ return cl;
522+}
523+
524+mp_limb_t
525+mpn_submul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl)
526+{
527+ mp_limb_t ul, cl, hpl, lpl, rl;
528+
529+ assert (n >= 1);
530+
531+ cl = 0;
532+ do
533+ {
534+ ul = *up++;
535+ gmp_umul_ppmm (hpl, lpl, ul, vl);
536+
537+ lpl += cl;
538+ cl = (lpl < cl) + hpl;
539+
540+ rl = *rp;
541+ lpl = rl - lpl;
542+ cl += lpl > rl;
543+ *rp++ = lpl;
544+ }
545+ while (--n != 0);
546+
547+ return cl;
548+}
549+
550+mp_limb_t
551+mpn_mul (mp_ptr rp, mp_srcptr up, mp_size_t un, mp_srcptr vp, mp_size_t vn)
552+{
553+ assert (un >= vn);
554+ assert (vn >= 1);
555+
556+ /* We first multiply by the low order limb. This result can be
557+ stored, not added, to rp. We also avoid a loop for zeroing this
558+ way. */
559+
560+ rp[un] = mpn_mul_1 (rp, up, un, vp[0]);
561+ rp += 1, vp += 1, vn -= 1;
562+
563+ /* Now accumulate the product of up[] and the next higher limb from
564+ vp[]. */
565+
566+ while (vn >= 1)
567+ {
568+ rp[un] = mpn_addmul_1 (rp, up, un, vp[0]);
569+ rp += 1, vp += 1, vn -= 1;
570+ }
571+ return rp[un - 1];
572+}
573+
574+void
575+mpn_mul_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n)
576+{
577+ mpn_mul (rp, ap, n, bp, n);
578+}
579+
580+void
581+mpn_sqr (mp_ptr rp, mp_srcptr ap, mp_size_t n)
582+{
583+ mpn_mul (rp, ap, n, ap, n);
584+}
585+
586+mp_limb_t
587+mpn_lshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt)
588+{
589+ mp_limb_t high_limb, low_limb;
590+ unsigned int tnc;
591+ mp_size_t i;
592+ mp_limb_t retval;
593+
594+ assert (n >= 1);
595+ assert (cnt >= 1);
596+ assert (cnt < GMP_LIMB_BITS);
597+
598+ up += n;
599+ rp += n;
600+
601+ tnc = GMP_LIMB_BITS - cnt;
602+ low_limb = *--up;
603+ retval = low_limb >> tnc;
604+ high_limb = (low_limb << cnt);
605+
606+ for (i = n; --i != 0;)
607+ {
608+ low_limb = *--up;
609+ *--rp = high_limb | (low_limb >> tnc);
610+ high_limb = (low_limb << cnt);
611+ }
612+ *--rp = high_limb;
613+
614+ return retval;
615+}
616+
617+mp_limb_t
618+mpn_rshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt)
619+{
620+ mp_limb_t high_limb, low_limb;
621+ unsigned int tnc;
622+ mp_size_t i;
623+ mp_limb_t retval;
624+
625+ assert (n >= 1);
626+ assert (cnt >= 1);
627+ assert (cnt < GMP_LIMB_BITS);
628+
629+ tnc = GMP_LIMB_BITS - cnt;
630+ high_limb = *up++;
631+ retval = (high_limb << tnc);
632+ low_limb = high_limb >> cnt;
633+
634+ for (i = n; --i != 0;)
635+ {
636+ high_limb = *up++;
637+ *rp++ = low_limb | (high_limb << tnc);
638+ low_limb = high_limb >> cnt;
639+ }
640+ *rp = low_limb;
641+
642+ return retval;
643+}
644+
645+static mp_bitcnt_t
646+mpn_common_scan (mp_limb_t limb, mp_size_t i, mp_srcptr up, mp_size_t un,
647+ mp_limb_t ux)
648+{
649+ unsigned cnt;
650+
651+ assert (ux == 0 || ux == GMP_LIMB_MAX);
652+ assert (0 <= i && i <= un );
653+
654+ while (limb == 0)
655+ {
656+ i++;
657+ if (i == un)
658+ return (ux == 0 ? ~(mp_bitcnt_t) 0 : un * GMP_LIMB_BITS);
659+ limb = ux ^ up[i];
660+ }
661+ gmp_ctz (cnt, limb);
662+ return (mp_bitcnt_t) i * GMP_LIMB_BITS + cnt;
663+}
664+
665+mp_bitcnt_t
666+mpn_scan1 (mp_srcptr ptr, mp_bitcnt_t bit)
667+{
668+ mp_size_t i;
669+ i = bit / GMP_LIMB_BITS;
670+
671+ return mpn_common_scan ( ptr[i] & (GMP_LIMB_MAX << (bit % GMP_LIMB_BITS)),
672+ i, ptr, i, 0);
673+}
674+
675+mp_bitcnt_t
676+mpn_scan0 (mp_srcptr ptr, mp_bitcnt_t bit)
677+{
678+ mp_size_t i;
679+ i = bit / GMP_LIMB_BITS;
680+
681+ return mpn_common_scan (~ptr[i] & (GMP_LIMB_MAX << (bit % GMP_LIMB_BITS)),
682+ i, ptr, i, GMP_LIMB_MAX);
683+}
684+
685+
686+/* MPN division interface. */
687+mp_limb_t
688+mpn_invert_3by2 (mp_limb_t u1, mp_limb_t u0)
689+{
690+ mp_limb_t r, p, m;
691+ unsigned ul, uh;
692+ unsigned ql, qh;
693+
694+ /* First, do a 2/1 inverse. */
695+ /* The inverse m is defined as floor( (B^2 - 1 - u1)/u1 ), so that 0 <
696+ * B^2 - (B + m) u1 <= u1 */
697+ assert (u1 >= GMP_LIMB_HIGHBIT);
698+
699+ ul = u1 & GMP_LLIMB_MASK;
700+ uh = u1 >> (GMP_LIMB_BITS / 2);
701+
702+ qh = ~u1 / uh;
703+ r = ((~u1 - (mp_limb_t) qh * uh) << (GMP_LIMB_BITS / 2)) | GMP_LLIMB_MASK;
704+
705+ p = (mp_limb_t) qh * ul;
706+ /* Adjustment steps taken from udiv_qrnnd_c */
707+ if (r < p)
708+ {
709+ qh--;
710+ r += u1;
711+ if (r >= u1) /* i.e. we didn't get carry when adding to r */
712+ if (r < p)
713+ {
714+ qh--;
715+ r += u1;
716+ }
717+ }
718+ r -= p;
719+
720+ /* Do a 3/2 division (with half limb size) */
721+ p = (r >> (GMP_LIMB_BITS / 2)) * qh + r;
722+ ql = (p >> (GMP_LIMB_BITS / 2)) + 1;
723+
724+ /* By the 3/2 method, we don't need the high half limb. */
725+ r = (r << (GMP_LIMB_BITS / 2)) + GMP_LLIMB_MASK - ql * u1;
726+
727+ if (r >= (p << (GMP_LIMB_BITS / 2)))
728+ {
729+ ql--;
730+ r += u1;
731+ }
732+ m = ((mp_limb_t) qh << (GMP_LIMB_BITS / 2)) + ql;
733+ if (r >= u1)
734+ {
735+ m++;
736+ r -= u1;
737+ }
738+
739+ if (u0 > 0)
740+ {
741+ mp_limb_t th, tl;
742+ r = ~r;
743+ r += u0;
744+ if (r < u0)
745+ {
746+ m--;
747+ if (r >= u1)
748+ {
749+ m--;
750+ r -= u1;
751+ }
752+ r -= u1;
753+ }
754+ gmp_umul_ppmm (th, tl, u0, m);
755+ r += th;
756+ if (r < th)
757+ {
758+ m--;
759+ m -= ((r > u1) | ((r == u1) & (tl > u0)));
760+ }
761+ }
762+
763+ return m;
764+}
765+
766+struct gmp_div_inverse
767+{
768+ /* Normalization shift count. */
769+ unsigned shift;
770+ /* Normalized divisor (d0 unused for mpn_div_qr_1) */
771+ mp_limb_t d1, d0;
772+ /* Inverse, for 2/1 or 3/2. */
773+ mp_limb_t di;
774+};
775+
776+static void
777+mpn_div_qr_1_invert (struct gmp_div_inverse *inv, mp_limb_t d)
778+{
779+ unsigned shift;
780+
781+ assert (d > 0);
782+ gmp_clz (shift, d);
783+ inv->shift = shift;
784+ inv->d1 = d << shift;
785+ inv->di = mpn_invert_limb (inv->d1);
786+}
787+
788+static void
789+mpn_div_qr_2_invert (struct gmp_div_inverse *inv,
790+ mp_limb_t d1, mp_limb_t d0)
791+{
792+ unsigned shift;
793+
794+ assert (d1 > 0);
795+ gmp_clz (shift, d1);
796+ inv->shift = shift;
797+ if (shift > 0)
798+ {
799+ d1 = (d1 << shift) | (d0 >> (GMP_LIMB_BITS - shift));
800+ d0 <<= shift;
801+ }
802+ inv->d1 = d1;
803+ inv->d0 = d0;
804+ inv->di = mpn_invert_3by2 (d1, d0);
805+}
806+
807+static void
808+mpn_div_qr_invert (struct gmp_div_inverse *inv,
809+ mp_srcptr dp, mp_size_t dn)
810+{
811+ assert (dn > 0);
812+
813+ if (dn == 1)
814+ mpn_div_qr_1_invert (inv, dp[0]);
815+ else if (dn == 2)
816+ mpn_div_qr_2_invert (inv, dp[1], dp[0]);
817+ else
818+ {
819+ unsigned shift;
820+ mp_limb_t d1, d0;
821+
822+ d1 = dp[dn-1];
823+ d0 = dp[dn-2];
824+ assert (d1 > 0);
825+ gmp_clz (shift, d1);
826+ inv->shift = shift;
827+ if (shift > 0)
828+ {
829+ d1 = (d1 << shift) | (d0 >> (GMP_LIMB_BITS - shift));
830+ d0 = (d0 << shift) | (dp[dn-3] >> (GMP_LIMB_BITS - shift));
831+ }
832+ inv->d1 = d1;
833+ inv->d0 = d0;
834+ inv->di = mpn_invert_3by2 (d1, d0);
835+ }
836+}
837+
838+/* Not matching current public gmp interface, rather corresponding to
839+ the sbpi1_div_* functions. */
840+static mp_limb_t
841+mpn_div_qr_1_preinv (mp_ptr qp, mp_srcptr np, mp_size_t nn,
842+ const struct gmp_div_inverse *inv)
843+{
844+ mp_limb_t d, di;
845+ mp_limb_t r;
846+ mp_ptr tp = NULL;
847+
848+ if (inv->shift > 0)
849+ {
850+ tp = gmp_xalloc_limbs (nn);
851+ r = mpn_lshift (tp, np, nn, inv->shift);
852+ np = tp;
853+ }
854+ else
855+ r = 0;
856+
857+ d = inv->d1;
858+ di = inv->di;
859+ while (nn-- > 0)
860+ {
861+ mp_limb_t q;
862+
863+ gmp_udiv_qrnnd_preinv (q, r, r, np[nn], d, di);
864+ if (qp)
865+ qp[nn] = q;
866+ }
867+ if (inv->shift > 0)
868+ gmp_free (tp);
869+
870+ return r >> inv->shift;
871+}
872+
873+static mp_limb_t
874+mpn_div_qr_1 (mp_ptr qp, mp_srcptr np, mp_size_t nn, mp_limb_t d)
875+{
876+ assert (d > 0);
877+
878+ /* Special case for powers of two. */
879+ if ((d & (d-1)) == 0)
880+ {
881+ mp_limb_t r = np[0] & (d-1);
882+ if (qp)
883+ {
884+ if (d <= 1)
885+ mpn_copyi (qp, np, nn);
886+ else
887+ {
888+ unsigned shift;
889+ gmp_ctz (shift, d);
890+ mpn_rshift (qp, np, nn, shift);
891+ }
892+ }
893+ return r;
894+ }
895+ else
896+ {
897+ struct gmp_div_inverse inv;
898+ mpn_div_qr_1_invert (&inv, d);
899+ return mpn_div_qr_1_preinv (qp, np, nn, &inv);
900+ }
901+}
902+
903+static void
904+mpn_div_qr_2_preinv (mp_ptr qp, mp_ptr rp, mp_srcptr np, mp_size_t nn,
905+ const struct gmp_div_inverse *inv)
906+{
907+ unsigned shift;
908+ mp_size_t i;
909+ mp_limb_t d1, d0, di, r1, r0;
910+ mp_ptr tp;
911+
912+ assert (nn >= 2);
913+ shift = inv->shift;
914+ d1 = inv->d1;
915+ d0 = inv->d0;
916+ di = inv->di;
917+
918+ if (shift > 0)
919+ {
920+ tp = gmp_xalloc_limbs (nn);
921+ r1 = mpn_lshift (tp, np, nn, shift);
922+ np = tp;
923+ }
924+ else
925+ r1 = 0;
926+
927+ r0 = np[nn - 1];
928+
929+ i = nn - 2;
930+ do
931+ {
932+ mp_limb_t n0, q;
933+ n0 = np[i];
934+ gmp_udiv_qr_3by2 (q, r1, r0, r1, r0, n0, d1, d0, di);
935+
936+ if (qp)
937+ qp[i] = q;
938+ }
939+ while (--i >= 0);
940+
941+ if (shift > 0)
942+ {
943+ assert ((r0 << (GMP_LIMB_BITS - shift)) == 0);
944+ r0 = (r0 >> shift) | (r1 << (GMP_LIMB_BITS - shift));
945+ r1 >>= shift;
946+
947+ gmp_free (tp);
948+ }
949+
950+ rp[1] = r1;
951+ rp[0] = r0;
952+}
953+
954+#if 0
955+static void
956+mpn_div_qr_2 (mp_ptr qp, mp_ptr rp, mp_srcptr np, mp_size_t nn,
957+ mp_limb_t d1, mp_limb_t d0)
958+{
959+ struct gmp_div_inverse inv;
960+ assert (nn >= 2);
961+
962+ mpn_div_qr_2_invert (&inv, d1, d0);
963+ mpn_div_qr_2_preinv (qp, rp, np, nn, &inv);
964+}
965+#endif
966+
967+static void
968+mpn_div_qr_pi1 (mp_ptr qp,
969+ mp_ptr np, mp_size_t nn, mp_limb_t n1,
970+ mp_srcptr dp, mp_size_t dn,
971+ mp_limb_t dinv)
972+{
973+ mp_size_t i;
974+
975+ mp_limb_t d1, d0;
976+ mp_limb_t cy, cy1;
977+ mp_limb_t q;
978+
979+ assert (dn > 2);
980+ assert (nn >= dn);
981+
982+ d1 = dp[dn - 1];
983+ d0 = dp[dn - 2];
984+
985+ assert ((d1 & GMP_LIMB_HIGHBIT) != 0);
986+ /* Iteration variable is the index of the q limb.
987+ *
988+ * We divide <n1, np[dn-1+i], np[dn-2+i], np[dn-3+i],..., np[i]>
989+ * by <d1, d0, dp[dn-3], ..., dp[0] >
990+ */
991+
992+ i = nn - dn;
993+ do
994+ {
995+ mp_limb_t n0 = np[dn-1+i];
996+
997+ if (n1 == d1 && n0 == d0)
998+ {
999+ q = GMP_LIMB_MAX;
1000+ mpn_submul_1 (np+i, dp, dn, q);
1001+ n1 = np[dn-1+i]; /* update n1, last loop's value will now be invalid */
1002+ }
1003+ else
1004+ {
1005+ gmp_udiv_qr_3by2 (q, n1, n0, n1, n0, np[dn-2+i], d1, d0, dinv);
1006+
1007+ cy = mpn_submul_1 (np + i, dp, dn-2, q);
1008+
1009+ cy1 = n0 < cy;
1010+ n0 = n0 - cy;
1011+ cy = n1 < cy1;
1012+ n1 = n1 - cy1;
1013+ np[dn-2+i] = n0;
1014+
1015+ if (cy != 0)
1016+ {
1017+ n1 += d1 + mpn_add_n (np + i, np + i, dp, dn - 1);
1018+ q--;
1019+ }
1020+ }
1021+
1022+ if (qp)
1023+ qp[i] = q;
1024+ }
1025+ while (--i >= 0);
1026+
1027+ np[dn - 1] = n1;
1028+}
1029+
1030+static void
1031+mpn_div_qr_preinv (mp_ptr qp, mp_ptr np, mp_size_t nn,
1032+ mp_srcptr dp, mp_size_t dn,
1033+ const struct gmp_div_inverse *inv)
1034+{
1035+ assert (dn > 0);
1036+ assert (nn >= dn);
1037+
1038+ if (dn == 1)
1039+ np[0] = mpn_div_qr_1_preinv (qp, np, nn, inv);
1040+ else if (dn == 2)
1041+ mpn_div_qr_2_preinv (qp, np, np, nn, inv);
1042+ else
1043+ {
1044+ mp_limb_t nh;
1045+ unsigned shift;
1046+
1047+ assert (inv->d1 == dp[dn-1]);
1048+ assert (inv->d0 == dp[dn-2]);
1049+ assert ((inv->d1 & GMP_LIMB_HIGHBIT) != 0);
1050+
1051+ shift = inv->shift;
1052+ if (shift > 0)
1053+ nh = mpn_lshift (np, np, nn, shift);
1054+ else
1055+ nh = 0;
1056+
1057+ mpn_div_qr_pi1 (qp, np, nn, nh, dp, dn, inv->di);
1058+
1059+ if (shift > 0)
1060+ gmp_assert_nocarry (mpn_rshift (np, np, dn, shift));
1061+ }
1062+}
1063+
1064+static void
1065+mpn_div_qr (mp_ptr qp, mp_ptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn)
1066+{
1067+ struct gmp_div_inverse inv;
1068+ mp_ptr tp = NULL;
1069+
1070+ assert (dn > 0);
1071+ assert (nn >= dn);
1072+
1073+ mpn_div_qr_invert (&inv, dp, dn);
1074+ if (dn > 2 && inv.shift > 0)
1075+ {
1076+ tp = gmp_xalloc_limbs (dn);
1077+ gmp_assert_nocarry (mpn_lshift (tp, dp, dn, inv.shift));
1078+ dp = tp;
1079+ }
1080+ mpn_div_qr_preinv (qp, np, nn, dp, dn, &inv);
1081+ if (tp)
1082+ gmp_free (tp);
1083+}
1084+
1085+
1086+/* MPN base conversion. */
1087+static unsigned
1088+mpn_base_power_of_two_p (unsigned b)
1089+{
1090+ switch (b)
1091+ {
1092+ case 2: return 1;
1093+ case 4: return 2;
1094+ case 8: return 3;
1095+ case 16: return 4;
1096+ case 32: return 5;
1097+ case 64: return 6;
1098+ case 128: return 7;
1099+ case 256: return 8;
1100+ default: return 0;
1101+ }
1102+}
1103+
1104+struct mpn_base_info
1105+{
1106+ /* bb is the largest power of the base which fits in one limb, and
1107+ exp is the corresponding exponent. */
1108+ unsigned exp;
1109+ mp_limb_t bb;
1110+};
1111+
1112+static void
1113+mpn_get_base_info (struct mpn_base_info *info, mp_limb_t b)
1114+{
1115+ mp_limb_t m;
1116+ mp_limb_t p;
1117+ unsigned exp;
1118+
1119+ m = GMP_LIMB_MAX / b;
1120+ for (exp = 1, p = b; p <= m; exp++)
1121+ p *= b;
1122+
1123+ info->exp = exp;
1124+ info->bb = p;
1125+}
1126+
1127+static mp_bitcnt_t
1128+mpn_limb_size_in_base_2 (mp_limb_t u)
1129+{
1130+ unsigned shift;
1131+
1132+ assert (u > 0);
1133+ gmp_clz (shift, u);
1134+ return GMP_LIMB_BITS - shift;
1135+}
1136+
1137+static size_t
1138+mpn_get_str_bits (unsigned char *sp, unsigned bits, mp_srcptr up, mp_size_t un)
1139+{
1140+ unsigned char mask;
1141+ size_t sn, j;
1142+ mp_size_t i;
1143+ int shift;
1144+
1145+ sn = ((un - 1) * GMP_LIMB_BITS + mpn_limb_size_in_base_2 (up[un-1])
1146+ + bits - 1) / bits;
1147+
1148+ mask = (1U << bits) - 1;
1149+
1150+ for (i = 0, j = sn, shift = 0; j-- > 0;)
1151+ {
1152+ unsigned char digit = up[i] >> shift;
1153+
1154+ shift += bits;
1155+
1156+ if (shift >= GMP_LIMB_BITS && ++i < un)
1157+ {
1158+ shift -= GMP_LIMB_BITS;
1159+ digit |= up[i] << (bits - shift);
1160+ }
1161+ sp[j] = digit & mask;
1162+ }
1163+ return sn;
1164+}
1165+
1166+/* We generate digits from the least significant end, and reverse at
1167+ the end. */
1168+static size_t
1169+mpn_limb_get_str (unsigned char *sp, mp_limb_t w,
1170+ const struct gmp_div_inverse *binv)
1171+{
1172+ mp_size_t i;
1173+ for (i = 0; w > 0; i++)
1174+ {
1175+ mp_limb_t h, l, r;
1176+
1177+ h = w >> (GMP_LIMB_BITS - binv->shift);
1178+ l = w << binv->shift;
1179+
1180+ gmp_udiv_qrnnd_preinv (w, r, h, l, binv->d1, binv->di);
1181+ assert ( (r << (GMP_LIMB_BITS - binv->shift)) == 0);
1182+ r >>= binv->shift;
1183+
1184+ sp[i] = r;
1185+ }
1186+ return i;
1187+}
1188+
1189+static size_t
1190+mpn_get_str_other (unsigned char *sp,
1191+ int base, const struct mpn_base_info *info,
1192+ mp_ptr up, mp_size_t un)
1193+{
1194+ struct gmp_div_inverse binv;
1195+ size_t sn;
1196+ size_t i;
1197+
1198+ mpn_div_qr_1_invert (&binv, base);
1199+
1200+ sn = 0;
1201+
1202+ if (un > 1)
1203+ {
1204+ struct gmp_div_inverse bbinv;
1205+ mpn_div_qr_1_invert (&bbinv, info->bb);
1206+
1207+ do
1208+ {
1209+ mp_limb_t w;
1210+ size_t done;
1211+ w = mpn_div_qr_1_preinv (up, up, un, &bbinv);
1212+ un -= (up[un-1] == 0);
1213+ done = mpn_limb_get_str (sp + sn, w, &binv);
1214+
1215+ for (sn += done; done < info->exp; done++)
1216+ sp[sn++] = 0;
1217+ }
1218+ while (un > 1);
1219+ }
1220+ sn += mpn_limb_get_str (sp + sn, up[0], &binv);
1221+
1222+ /* Reverse order */
1223+ for (i = 0; 2*i + 1 < sn; i++)
1224+ {
1225+ unsigned char t = sp[i];
1226+ sp[i] = sp[sn - i - 1];
1227+ sp[sn - i - 1] = t;
1228+ }
1229+
1230+ return sn;
1231+}
1232+
1233+size_t
1234+mpn_get_str (unsigned char *sp, int base, mp_ptr up, mp_size_t un)
1235+{
1236+ unsigned bits;
1237+
1238+ assert (un > 0);
1239+ assert (up[un-1] > 0);
1240+
1241+ bits = mpn_base_power_of_two_p (base);
1242+ if (bits)
1243+ return mpn_get_str_bits (sp, bits, up, un);
1244+ else
1245+ {
1246+ struct mpn_base_info info;
1247+
1248+ mpn_get_base_info (&info, base);
1249+ return mpn_get_str_other (sp, base, &info, up, un);
1250+ }
1251+}
1252+
1253+static mp_size_t
1254+mpn_set_str_bits (mp_ptr rp, const unsigned char *sp, size_t sn,
1255+ unsigned bits)
1256+{
1257+ mp_size_t rn;
1258+ size_t j;
1259+ unsigned shift;
1260+
1261+ for (j = sn, rn = 0, shift = 0; j-- > 0; )
1262+ {
1263+ if (shift == 0)
1264+ {
1265+ rp[rn++] = sp[j];
1266+ shift += bits;
1267+ }
1268+ else
1269+ {
1270+ rp[rn-1] |= (mp_limb_t) sp[j] << shift;
1271+ shift += bits;
1272+ if (shift >= GMP_LIMB_BITS)
1273+ {
1274+ shift -= GMP_LIMB_BITS;
1275+ if (shift > 0)
1276+ rp[rn++] = (mp_limb_t) sp[j] >> (bits - shift);
1277+ }
1278+ }
1279+ }
1280+ rn = mpn_normalized_size (rp, rn);
1281+ return rn;
1282+}
1283+
1284+static mp_size_t
1285+mpn_set_str_other (mp_ptr rp, const unsigned char *sp, size_t sn,
1286+ mp_limb_t b, const struct mpn_base_info *info)
1287+{
1288+ mp_size_t rn;
1289+ mp_limb_t w;
1290+ unsigned k;
1291+ size_t j;
1292+
1293+ k = 1 + (sn - 1) % info->exp;
1294+
1295+ j = 0;
1296+ w = sp[j++];
1297+ for (; --k > 0; )
1298+ w = w * b + sp[j++];
1299+
1300+ rp[0] = w;
1301+
1302+ for (rn = (w > 0); j < sn;)
1303+ {
1304+ mp_limb_t cy;
1305+
1306+ w = sp[j++];
1307+ for (k = 1; k < info->exp; k++)
1308+ w = w * b + sp[j++];
1309+
1310+ cy = mpn_mul_1 (rp, rp, rn, info->bb);
1311+ cy += mpn_add_1 (rp, rp, rn, w);
1312+ if (cy > 0)
1313+ rp[rn++] = cy;
1314+ }
1315+ assert (j == sn);
1316+
1317+ return rn;
1318+}
1319+
1320+mp_size_t
1321+mpn_set_str (mp_ptr rp, const unsigned char *sp, size_t sn, int base)
1322+{
1323+ unsigned bits;
1324+
1325+ if (sn == 0)
1326+ return 0;
1327+
1328+ bits = mpn_base_power_of_two_p (base);
1329+ if (bits)
1330+ return mpn_set_str_bits (rp, sp, sn, bits);
1331+ else
1332+ {
1333+ struct mpn_base_info info;
1334+
1335+ mpn_get_base_info (&info, base);
1336+ return mpn_set_str_other (rp, sp, sn, base, &info);
1337+ }
1338+}
1339+
1340+
1341+/* MPZ interface */
1342+void
1343+mpz_init (mpz_t r)
1344+{
1345+ r->_mp_alloc = 1;
1346+ r->_mp_size = 0;
1347+ r->_mp_d = gmp_xalloc_limbs (1);
1348+}
1349+
1350+/* The utility of this function is a bit limited, since many functions
1351+ assigns the result variable using mpz_swap. */
1352+void
1353+mpz_init2 (mpz_t r, mp_bitcnt_t bits)
1354+{
1355+ mp_size_t rn;
1356+
1357+ bits -= (bits != 0); /* Round down, except if 0 */
1358+ rn = 1 + bits / GMP_LIMB_BITS;
1359+
1360+ r->_mp_alloc = rn;
1361+ r->_mp_size = 0;
1362+ r->_mp_d = gmp_xalloc_limbs (rn);
1363+}
1364+
1365+void
1366+mpz_clear (mpz_t r)
1367+{
1368+ gmp_free (r->_mp_d);
1369+}
1370+
1371+static void *
1372+mpz_realloc (mpz_t r, mp_size_t size)
1373+{
1374+ size = GMP_MAX (size, 1);
1375+
1376+ r->_mp_d = gmp_xrealloc_limbs (r->_mp_d, size);
1377+ r->_mp_alloc = size;
1378+
1379+ if (GMP_ABS (r->_mp_size) > size)
1380+ r->_mp_size = 0;
1381+
1382+ return r->_mp_d;
1383+}
1384+
1385+/* Realloc for an mpz_t WHAT if it has less than NEEDED limbs. */
1386+#define MPZ_REALLOC(z,n) ((n) > (z)->_mp_alloc \
1387+ ? mpz_realloc(z,n) \
1388+ : (z)->_mp_d)
1389+
1390+/* MPZ assignment and basic conversions. */
1391+void
1392+mpz_set_si (mpz_t r, signed long int x)
1393+{
1394+ if (x >= 0)
1395+ mpz_set_ui (r, x);
1396+ else /* (x < 0) */
1397+ {
1398+ r->_mp_size = -1;
1399+ r->_mp_d[0] = GMP_NEG_CAST (unsigned long int, x);
1400+ }
1401+}
1402+
1403+void
1404+mpz_set_ui (mpz_t r, unsigned long int x)
1405+{
1406+ if (x > 0)
1407+ {
1408+ r->_mp_size = 1;
1409+ r->_mp_d[0] = x;
1410+ }
1411+ else
1412+ r->_mp_size = 0;
1413+}
1414+
1415+void
1416+mpz_set (mpz_t r, const mpz_t x)
1417+{
1418+ /* Allow the NOP r == x */
1419+ if (r != x)
1420+ {
1421+ mp_size_t n;
1422+ mp_ptr rp;
1423+
1424+ n = GMP_ABS (x->_mp_size);
1425+ rp = MPZ_REALLOC (r, n);
1426+
1427+ mpn_copyi (rp, x->_mp_d, n);
1428+ r->_mp_size = x->_mp_size;
1429+ }
1430+}
1431+
1432+void
1433+mpz_init_set_si (mpz_t r, signed long int x)
1434+{
1435+ mpz_init (r);
1436+ mpz_set_si (r, x);
1437+}
1438+
1439+void
1440+mpz_init_set_ui (mpz_t r, unsigned long int x)
1441+{
1442+ mpz_init (r);
1443+ mpz_set_ui (r, x);
1444+}
1445+
1446+void
1447+mpz_init_set (mpz_t r, const mpz_t x)
1448+{
1449+ mpz_init (r);
1450+ mpz_set (r, x);
1451+}
1452+
1453+int
1454+mpz_fits_slong_p (const mpz_t u)
1455+{
1456+ mp_size_t us = u->_mp_size;
1457+
1458+ if (us == 0)
1459+ return 1;
1460+ else if (us == 1)
1461+ return u->_mp_d[0] < GMP_LIMB_HIGHBIT;
1462+ else if (us == -1)
1463+ return u->_mp_d[0] <= GMP_LIMB_HIGHBIT;
1464+ else
1465+ return 0;
1466+}
1467+
1468+int
1469+mpz_fits_ulong_p (const mpz_t u)
1470+{
1471+ mp_size_t us = u->_mp_size;
1472+
1473+ return (us == (us > 0));
1474+}
1475+
1476+long int
1477+mpz_get_si (const mpz_t u)
1478+{
1479+ mp_size_t us = u->_mp_size;
1480+
1481+ if (us > 0)
1482+ return (long) (u->_mp_d[0] & ~GMP_LIMB_HIGHBIT);
1483+ else if (us < 0)
1484+ return (long) (- u->_mp_d[0] | GMP_LIMB_HIGHBIT);
1485+ else
1486+ return 0;
1487+}
1488+
1489+unsigned long int
1490+mpz_get_ui (const mpz_t u)
1491+{
1492+ return u->_mp_size == 0 ? 0 : u->_mp_d[0];
1493+}
1494+
1495+size_t
1496+mpz_size (const mpz_t u)
1497+{
1498+ return GMP_ABS (u->_mp_size);
1499+}
1500+
1501+mp_limb_t
1502+mpz_getlimbn (const mpz_t u, mp_size_t n)
1503+{
1504+ if (n >= 0 && n < GMP_ABS (u->_mp_size))
1505+ return u->_mp_d[n];
1506+ else
1507+ return 0;
1508+}
1509+
1510+
1511+/* Conversions and comparison to double. */
1512+void
1513+mpz_set_d (mpz_t r, double x)
1514+{
1515+ int sign;
1516+ mp_ptr rp;
1517+ mp_size_t rn, i;
1518+ double B;
1519+ double Bi;
1520+ mp_limb_t f;
1521+
1522+ /* x != x is true when x is a NaN, and x == x * 0.5 is true when x is
1523+ zero or infinity. */
1524+ if (x != x || x == x * 0.5)
1525+ {
1526+ r->_mp_size = 0;
1527+ return;
1528+ }
1529+
1530+ sign = x < 0.0 ;
1531+ if (sign)
1532+ x = - x;
1533+
1534+ if (x < 1.0)
1535+ {
1536+ r->_mp_size = 0;
1537+ return;
1538+ }
1539+ B = 2.0 * (double) GMP_LIMB_HIGHBIT;
1540+ Bi = 1.0 / B;
1541+ for (rn = 1; x >= B; rn++)
1542+ x *= Bi;
1543+
1544+ rp = MPZ_REALLOC (r, rn);
1545+
1546+ f = (mp_limb_t) x;
1547+ x -= f;
1548+ assert (x < 1.0);
1549+ i = rn-1;
1550+ rp[i] = f;
1551+ while (--i >= 0)
1552+ {
1553+ x = B * x;
1554+ f = (mp_limb_t) x;
1555+ x -= f;
1556+ assert (x < 1.0);
1557+ rp[i] = f;
1558+ }
1559+
1560+ r->_mp_size = sign ? - rn : rn;
1561+}
1562+
1563+void
1564+mpz_init_set_d (mpz_t r, double x)
1565+{
1566+ mpz_init (r);
1567+ mpz_set_d (r, x);
1568+}
1569+
1570+double
1571+mpz_get_d (const mpz_t u)
1572+{
1573+ mp_size_t un;
1574+ double x;
1575+ double B = 2.0 * (double) GMP_LIMB_HIGHBIT;
1576+
1577+ un = GMP_ABS (u->_mp_size);
1578+
1579+ if (un == 0)
1580+ return 0.0;
1581+
1582+ x = u->_mp_d[--un];
1583+ while (un > 0)
1584+ x = B*x + u->_mp_d[--un];
1585+
1586+ if (u->_mp_size < 0)
1587+ x = -x;
1588+
1589+ return x;
1590+}
1591+
1592+int
1593+mpz_cmpabs_d (const mpz_t x, double d)
1594+{
1595+ mp_size_t xn;
1596+ double B, Bi;
1597+ mp_size_t i;
1598+
1599+ xn = x->_mp_size;
1600+ d = GMP_ABS (d);
1601+
1602+ if (xn != 0)
1603+ {
1604+ xn = GMP_ABS (xn);
1605+
1606+ B = 2.0 * (double) GMP_LIMB_HIGHBIT;
1607+ Bi = 1.0 / B;
1608+
1609+ /* Scale d so it can be compared with the top limb. */
1610+ for (i = 1; i < xn; i++)
1611+ d *= Bi;
1612+
1613+ if (d >= B)
1614+ return -1;
1615+
1616+ /* Compare floor(d) to top limb, subtract and cancel when equal. */
1617+ for (i = xn; i-- > 0;)
1618+ {
1619+ mp_limb_t f, xl;
1620+
1621+ f = (mp_limb_t) d;
1622+ xl = x->_mp_d[i];
1623+ if (xl > f)
1624+ return 1;
1625+ else if (xl < f)
1626+ return -1;
1627+ d = B * (d - f);
1628+ }
1629+ }
1630+ return - (d > 0.0);
1631+}
1632+
1633+int
1634+mpz_cmp_d (const mpz_t x, double d)
1635+{
1636+ if (x->_mp_size < 0)
1637+ {
1638+ if (d >= 0.0)
1639+ return -1;
1640+ else
1641+ return -mpz_cmpabs_d (x, d);
1642+ }
1643+ else
1644+ {
1645+ if (d < 0.0)
1646+ return 1;
1647+ else
1648+ return mpz_cmpabs_d (x, d);
1649+ }
1650+}
1651+
1652+
1653+/* MPZ comparisons and the like. */
1654+int
1655+mpz_sgn (const mpz_t u)
1656+{
1657+ mp_size_t usize = u->_mp_size;
1658+
1659+ return (usize > 0) - (usize < 0);
1660+}
1661+
1662+int
1663+mpz_cmp_si (const mpz_t u, long v)
1664+{
1665+ mp_size_t usize = u->_mp_size;
1666+
1667+ if (usize < -1)
1668+ return -1;
1669+ else if (v >= 0)
1670+ return mpz_cmp_ui (u, v);
1671+ else if (usize >= 0)
1672+ return 1;
1673+ else /* usize == -1 */
1674+ {
1675+ mp_limb_t ul = u->_mp_d[0];
1676+ if ((mp_limb_t)GMP_NEG_CAST (unsigned long int, v) < ul)
1677+ return -1;
1678+ else
1679+ return (mp_limb_t)GMP_NEG_CAST (unsigned long int, v) > ul;
1680+ }
1681+}
1682+
1683+int
1684+mpz_cmp_ui (const mpz_t u, unsigned long v)
1685+{
1686+ mp_size_t usize = u->_mp_size;
1687+
1688+ if (usize > 1)
1689+ return 1;
1690+ else if (usize < 0)
1691+ return -1;
1692+ else
1693+ {
1694+ mp_limb_t ul = (usize > 0) ? u->_mp_d[0] : 0;
1695+ return (ul > v) - (ul < v);
1696+ }
1697+}
1698+
1699+int
1700+mpz_cmp (const mpz_t a, const mpz_t b)
1701+{
1702+ mp_size_t asize = a->_mp_size;
1703+ mp_size_t bsize = b->_mp_size;
1704+
1705+ if (asize != bsize)
1706+ return (asize < bsize) ? -1 : 1;
1707+ else if (asize >= 0)
1708+ return mpn_cmp (a->_mp_d, b->_mp_d, asize);
1709+ else
1710+ return mpn_cmp (b->_mp_d, a->_mp_d, -asize);
1711+}
1712+
1713+int
1714+mpz_cmpabs_ui (const mpz_t u, unsigned long v)
1715+{
1716+ mp_size_t un = GMP_ABS (u->_mp_size);
1717+ mp_limb_t ul;
1718+
1719+ if (un > 1)
1720+ return 1;
1721+
1722+ ul = (un == 1) ? u->_mp_d[0] : 0;
1723+
1724+ return (ul > v) - (ul < v);
1725+}
1726+
1727+int
1728+mpz_cmpabs (const mpz_t u, const mpz_t v)
1729+{
1730+ return mpn_cmp4 (u->_mp_d, GMP_ABS (u->_mp_size),
1731+ v->_mp_d, GMP_ABS (v->_mp_size));
1732+}
1733+
1734+void
1735+mpz_abs (mpz_t r, const mpz_t u)
1736+{
1737+ if (r != u)
1738+ mpz_set (r, u);
1739+
1740+ r->_mp_size = GMP_ABS (r->_mp_size);
1741+}
1742+
1743+void
1744+mpz_neg (mpz_t r, const mpz_t u)
1745+{
1746+ if (r != u)
1747+ mpz_set (r, u);
1748+
1749+ r->_mp_size = -r->_mp_size;
1750+}
1751+
1752+void
1753+mpz_swap (mpz_t u, mpz_t v)
1754+{
1755+ MP_SIZE_T_SWAP (u->_mp_size, v->_mp_size);
1756+ MP_SIZE_T_SWAP (u->_mp_alloc, v->_mp_alloc);
1757+ MP_PTR_SWAP (u->_mp_d, v->_mp_d);
1758+}
1759+
1760+
1761+/* MPZ addition and subtraction */
1762+
1763+/* Adds to the absolute value. Returns new size, but doesn't store it. */
1764+static mp_size_t
1765+mpz_abs_add_ui (mpz_t r, const mpz_t a, unsigned long b)
1766+{
1767+ mp_size_t an;
1768+ mp_ptr rp;
1769+ mp_limb_t cy;
1770+
1771+ an = GMP_ABS (a->_mp_size);
1772+ if (an == 0)
1773+ {
1774+ r->_mp_d[0] = b;
1775+ return b > 0;
1776+ }
1777+
1778+ rp = MPZ_REALLOC (r, an + 1);
1779+
1780+ cy = mpn_add_1 (rp, a->_mp_d, an, b);
1781+ rp[an] = cy;
1782+ an += cy;
1783+
1784+ return an;
1785+}
1786+
1787+/* Subtract from the absolute value. Returns new size, (or -1 on underflow),
1788+ but doesn't store it. */
1789+static mp_size_t
1790+mpz_abs_sub_ui (mpz_t r, const mpz_t a, unsigned long b)
1791+{
1792+ mp_size_t an = GMP_ABS (a->_mp_size);
1793+ mp_ptr rp = MPZ_REALLOC (r, an);
1794+
1795+ if (an == 0)
1796+ {
1797+ rp[0] = b;
1798+ return -(b > 0);
1799+ }
1800+ else if (an == 1 && a->_mp_d[0] < b)
1801+ {
1802+ rp[0] = b - a->_mp_d[0];
1803+ return -1;
1804+ }
1805+ else
1806+ {
1807+ gmp_assert_nocarry (mpn_sub_1 (rp, a->_mp_d, an, b));
1808+ return mpn_normalized_size (rp, an);
1809+ }
1810+}
1811+
1812+void
1813+mpz_add_ui (mpz_t r, const mpz_t a, unsigned long b)
1814+{
1815+ if (a->_mp_size >= 0)
1816+ r->_mp_size = mpz_abs_add_ui (r, a, b);
1817+ else
1818+ r->_mp_size = -mpz_abs_sub_ui (r, a, b);
1819+}
1820+
1821+void
1822+mpz_sub_ui (mpz_t r, const mpz_t a, unsigned long b)
1823+{
1824+ if (a->_mp_size < 0)
1825+ r->_mp_size = -mpz_abs_add_ui (r, a, b);
1826+ else
1827+ r->_mp_size = mpz_abs_sub_ui (r, a, b);
1828+}
1829+
1830+void
1831+mpz_ui_sub (mpz_t r, unsigned long a, const mpz_t b)
1832+{
1833+ if (b->_mp_size < 0)
1834+ r->_mp_size = mpz_abs_add_ui (r, b, a);
1835+ else
1836+ r->_mp_size = -mpz_abs_sub_ui (r, b, a);
1837+}
1838+
1839+static mp_size_t
1840+mpz_abs_add (mpz_t r, const mpz_t a, const mpz_t b)
1841+{
1842+ mp_size_t an = GMP_ABS (a->_mp_size);
1843+ mp_size_t bn = GMP_ABS (b->_mp_size);
1844+ mp_ptr rp;
1845+ mp_limb_t cy;
1846+
1847+ if (an < bn)
1848+ {
1849+ MPZ_SRCPTR_SWAP (a, b);
1850+ MP_SIZE_T_SWAP (an, bn);
1851+ }
1852+
1853+ rp = MPZ_REALLOC (r, an + 1);
1854+ cy = mpn_add (rp, a->_mp_d, an, b->_mp_d, bn);
1855+
1856+ rp[an] = cy;
1857+
1858+ return an + cy;
1859+}
1860+
1861+static mp_size_t
1862+mpz_abs_sub (mpz_t r, const mpz_t a, const mpz_t b)
1863+{
1864+ mp_size_t an = GMP_ABS (a->_mp_size);
1865+ mp_size_t bn = GMP_ABS (b->_mp_size);
1866+ int cmp;
1867+ mp_ptr rp;
1868+
1869+ cmp = mpn_cmp4 (a->_mp_d, an, b->_mp_d, bn);
1870+ if (cmp > 0)
1871+ {
1872+ rp = MPZ_REALLOC (r, an);
1873+ gmp_assert_nocarry (mpn_sub (rp, a->_mp_d, an, b->_mp_d, bn));
1874+ return mpn_normalized_size (rp, an);
1875+ }
1876+ else if (cmp < 0)
1877+ {
1878+ rp = MPZ_REALLOC (r, bn);
1879+ gmp_assert_nocarry (mpn_sub (rp, b->_mp_d, bn, a->_mp_d, an));
1880+ return -mpn_normalized_size (rp, bn);
1881+ }
1882+ else
1883+ return 0;
1884+}
1885+
1886+void
1887+mpz_add (mpz_t r, const mpz_t a, const mpz_t b)
1888+{
1889+ mp_size_t rn;
1890+
1891+ if ( (a->_mp_size ^ b->_mp_size) >= 0)
1892+ rn = mpz_abs_add (r, a, b);
1893+ else
1894+ rn = mpz_abs_sub (r, a, b);
1895+
1896+ r->_mp_size = a->_mp_size >= 0 ? rn : - rn;
1897+}
1898+
1899+void
1900+mpz_sub (mpz_t r, const mpz_t a, const mpz_t b)
1901+{
1902+ mp_size_t rn;
1903+
1904+ if ( (a->_mp_size ^ b->_mp_size) >= 0)
1905+ rn = mpz_abs_sub (r, a, b);
1906+ else
1907+ rn = mpz_abs_add (r, a, b);
1908+
1909+ r->_mp_size = a->_mp_size >= 0 ? rn : - rn;
1910+}
1911+
1912+
1913+/* MPZ multiplication */
1914+void
1915+mpz_mul_si (mpz_t r, const mpz_t u, long int v)
1916+{
1917+ if (v < 0)
1918+ {
1919+ mpz_mul_ui (r, u, GMP_NEG_CAST (unsigned long int, v));
1920+ mpz_neg (r, r);
1921+ }
1922+ else
1923+ mpz_mul_ui (r, u, (unsigned long int) v);
1924+}
1925+
1926+void
1927+mpz_mul_ui (mpz_t r, const mpz_t u, unsigned long int v)
1928+{
1929+ mp_size_t un, us;
1930+ mp_ptr tp;
1931+ mp_limb_t cy;
1932+
1933+ us = u->_mp_size;
1934+
1935+ if (us == 0 || v == 0)
1936+ {
1937+ r->_mp_size = 0;
1938+ return;
1939+ }
1940+
1941+ un = GMP_ABS (us);
1942+
1943+ tp = MPZ_REALLOC (r, un + 1);
1944+ cy = mpn_mul_1 (tp, u->_mp_d, un, v);
1945+ tp[un] = cy;
1946+
1947+ un += (cy > 0);
1948+ r->_mp_size = (us < 0) ? - un : un;
1949+}
1950+
1951+void
1952+mpz_mul (mpz_t r, const mpz_t u, const mpz_t v)
1953+{
1954+ int sign;
1955+ mp_size_t un, vn, rn;
1956+ mpz_t t;
1957+ mp_ptr tp;
1958+
1959+ un = u->_mp_size;
1960+ vn = v->_mp_size;
1961+
1962+ if (un == 0 || vn == 0)
1963+ {
1964+ r->_mp_size = 0;
1965+ return;
1966+ }
1967+
1968+ sign = (un ^ vn) < 0;
1969+
1970+ un = GMP_ABS (un);
1971+ vn = GMP_ABS (vn);
1972+
1973+ mpz_init2 (t, (un + vn) * GMP_LIMB_BITS);
1974+
1975+ tp = t->_mp_d;
1976+ if (un >= vn)
1977+ mpn_mul (tp, u->_mp_d, un, v->_mp_d, vn);
1978+ else
1979+ mpn_mul (tp, v->_mp_d, vn, u->_mp_d, un);
1980+
1981+ rn = un + vn;
1982+ rn -= tp[rn-1] == 0;
1983+
1984+ t->_mp_size = sign ? - rn : rn;
1985+ mpz_swap (r, t);
1986+ mpz_clear (t);
1987+}
1988+
1989+void
1990+mpz_mul_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bits)
1991+{
1992+ mp_size_t un, rn;
1993+ mp_size_t limbs;
1994+ unsigned shift;
1995+ mp_ptr rp;
1996+
1997+ un = GMP_ABS (u->_mp_size);
1998+ if (un == 0)
1999+ {
2000+ r->_mp_size = 0;
2001+ return;
2002+ }
2003+
2004+ limbs = bits / GMP_LIMB_BITS;
2005+ shift = bits % GMP_LIMB_BITS;
2006+
2007+ rn = un + limbs + (shift > 0);
2008+ rp = MPZ_REALLOC (r, rn);
2009+ if (shift > 0)
2010+ {
2011+ mp_limb_t cy = mpn_lshift (rp + limbs, u->_mp_d, un, shift);
2012+ rp[rn-1] = cy;
2013+ rn -= (cy == 0);
2014+ }
2015+ else
2016+ mpn_copyd (rp + limbs, u->_mp_d, un);
2017+
2018+ while (limbs > 0)
2019+ rp[--limbs] = 0;
2020+
2021+ r->_mp_size = (u->_mp_size < 0) ? - rn : rn;
2022+}
2023+
2024+
2025+/* MPZ division */
2026+enum mpz_div_round_mode { GMP_DIV_FLOOR, GMP_DIV_CEIL, GMP_DIV_TRUNC };
2027+
2028+/* Allows q or r to be zero. Returns 1 iff remainder is non-zero. */
2029+static int
2030+mpz_div_qr (mpz_t q, mpz_t r,
2031+ const mpz_t n, const mpz_t d, enum mpz_div_round_mode mode)
2032+{
2033+ mp_size_t ns, ds, nn, dn, qs;
2034+ ns = n->_mp_size;
2035+ ds = d->_mp_size;
2036+
2037+ if (ds == 0)
2038+ gmp_die("mpz_div_qr: Divide by zero.");
2039+
2040+ if (ns == 0)
2041+ {
2042+ if (q)
2043+ q->_mp_size = 0;
2044+ if (r)
2045+ r->_mp_size = 0;
2046+ return 0;
2047+ }
2048+
2049+ nn = GMP_ABS (ns);
2050+ dn = GMP_ABS (ds);
2051+
2052+ qs = ds ^ ns;
2053+
2054+ if (nn < dn)
2055+ {
2056+ if (mode == GMP_DIV_CEIL && qs >= 0)
2057+ {
2058+ /* q = 1, r = n - d */
2059+ if (r)
2060+ mpz_sub (r, n, d);
2061+ if (q)
2062+ mpz_set_ui (q, 1);
2063+ }
2064+ else if (mode == GMP_DIV_FLOOR && qs < 0)
2065+ {
2066+ /* q = -1, r = n + d */
2067+ if (r)
2068+ mpz_add (r, n, d);
2069+ if (q)
2070+ mpz_set_si (q, -1);
2071+ }
2072+ else
2073+ {
2074+ /* q = 0, r = d */
2075+ if (r)
2076+ mpz_set (r, n);
2077+ if (q)
2078+ q->_mp_size = 0;
2079+ }
2080+ return 1;
2081+ }
2082+ else
2083+ {
2084+ mp_ptr np, qp;
2085+ mp_size_t qn, rn;
2086+ mpz_t tq, tr;
2087+
2088+ mpz_init (tr);
2089+ mpz_set (tr, n);
2090+ np = tr->_mp_d;
2091+
2092+ qn = nn - dn + 1;
2093+
2094+ if (q)
2095+ {
2096+ mpz_init2 (tq, qn * GMP_LIMB_BITS);
2097+ qp = tq->_mp_d;
2098+ }
2099+ else
2100+ qp = NULL;
2101+
2102+ mpn_div_qr (qp, np, nn, d->_mp_d, dn);
2103+
2104+ if (qp)
2105+ {
2106+ qn -= (qp[qn-1] == 0);
2107+
2108+ tq->_mp_size = qs < 0 ? -qn : qn;
2109+ }
2110+ rn = mpn_normalized_size (np, dn);
2111+ tr->_mp_size = ns < 0 ? - rn : rn;
2112+
2113+ if (mode == GMP_DIV_FLOOR && qs < 0 && rn != 0)
2114+ {
2115+ if (q)
2116+ mpz_sub_ui (tq, tq, 1);
2117+ if (r)
2118+ mpz_add (tr, tr, d);
2119+ }
2120+ else if (mode == GMP_DIV_CEIL && qs >= 0 && rn != 0)
2121+ {
2122+ if (q)
2123+ mpz_add_ui (tq, tq, 1);
2124+ if (r)
2125+ mpz_sub (tr, tr, d);
2126+ }
2127+
2128+ if (q)
2129+ {
2130+ mpz_swap (tq, q);
2131+ mpz_clear (tq);
2132+ }
2133+ if (r)
2134+ mpz_swap (tr, r);
2135+
2136+ mpz_clear (tr);
2137+
2138+ return rn != 0;
2139+ }
2140+}
2141+
2142+void
2143+mpz_cdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d)
2144+{
2145+ mpz_div_qr (q, r, n, d, GMP_DIV_CEIL);
2146+}
2147+
2148+void
2149+mpz_fdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d)
2150+{
2151+ mpz_div_qr (q, r, n, d, GMP_DIV_FLOOR);
2152+}
2153+
2154+void
2155+mpz_tdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d)
2156+{
2157+ mpz_div_qr (q, r, n, d, GMP_DIV_TRUNC);
2158+}
2159+
2160+void
2161+mpz_cdiv_q (mpz_t q, const mpz_t n, const mpz_t d)
2162+{
2163+ mpz_div_qr (q, NULL, n, d, GMP_DIV_CEIL);
2164+}
2165+
2166+void
2167+mpz_fdiv_q (mpz_t q, const mpz_t n, const mpz_t d)
2168+{
2169+ mpz_div_qr (q, NULL, n, d, GMP_DIV_FLOOR);
2170+}
2171+
2172+void
2173+mpz_tdiv_q (mpz_t q, const mpz_t n, const mpz_t d)
2174+{
2175+ mpz_div_qr (q, NULL, n, d, GMP_DIV_TRUNC);
2176+}
2177+
2178+void
2179+mpz_cdiv_r (mpz_t r, const mpz_t n, const mpz_t d)
2180+{
2181+ mpz_div_qr (NULL, r, n, d, GMP_DIV_CEIL);
2182+}
2183+
2184+void
2185+mpz_fdiv_r (mpz_t r, const mpz_t n, const mpz_t d)
2186+{
2187+ mpz_div_qr (NULL, r, n, d, GMP_DIV_FLOOR);
2188+}
2189+
2190+void
2191+mpz_tdiv_r (mpz_t r, const mpz_t n, const mpz_t d)
2192+{
2193+ mpz_div_qr (NULL, r, n, d, GMP_DIV_TRUNC);
2194+}
2195+
2196+void
2197+mpz_mod (mpz_t r, const mpz_t n, const mpz_t d)
2198+{
2199+ mpz_div_qr (NULL, r, n, d, d->_mp_size >= 0 ? GMP_DIV_FLOOR : GMP_DIV_CEIL);
2200+}
2201+
2202+static void
2203+mpz_div_q_2exp (mpz_t q, const mpz_t u, mp_bitcnt_t bit_index,
2204+ enum mpz_div_round_mode mode)
2205+{
2206+ mp_size_t un, qn;
2207+ mp_size_t limb_cnt;
2208+ mp_ptr qp;
2209+ int adjust;
2210+
2211+ un = u->_mp_size;
2212+ if (un == 0)
2213+ {
2214+ q->_mp_size = 0;
2215+ return;
2216+ }
2217+ limb_cnt = bit_index / GMP_LIMB_BITS;
2218+ qn = GMP_ABS (un) - limb_cnt;
2219+ bit_index %= GMP_LIMB_BITS;
2220+
2221+ if (mode == ((un > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* un != 0 here. */
2222+ /* Note: Below, the final indexing at limb_cnt is valid because at
2223+ that point we have qn > 0. */
2224+ adjust = (qn <= 0
2225+ || !mpn_zero_p (u->_mp_d, limb_cnt)
2226+ || (u->_mp_d[limb_cnt]
2227+ & (((mp_limb_t) 1 << bit_index) - 1)));
2228+ else
2229+ adjust = 0;
2230+
2231+ if (qn <= 0)
2232+ qn = 0;
2233+
2234+ else
2235+ {
2236+ qp = MPZ_REALLOC (q, qn);
2237+
2238+ if (bit_index != 0)
2239+ {
2240+ mpn_rshift (qp, u->_mp_d + limb_cnt, qn, bit_index);
2241+ qn -= qp[qn - 1] == 0;
2242+ }
2243+ else
2244+ {
2245+ mpn_copyi (qp, u->_mp_d + limb_cnt, qn);
2246+ }
2247+ }
2248+
2249+ q->_mp_size = qn;
2250+
2251+ if (adjust)
2252+ mpz_add_ui (q, q, 1);
2253+ if (un < 0)
2254+ mpz_neg (q, q);
2255+}
2256+
2257+static void
2258+mpz_div_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bit_index,
2259+ enum mpz_div_round_mode mode)
2260+{
2261+ mp_size_t us, un, rn;
2262+ mp_ptr rp;
2263+ mp_limb_t mask;
2264+
2265+ us = u->_mp_size;
2266+ if (us == 0 || bit_index == 0)
2267+ {
2268+ r->_mp_size = 0;
2269+ return;
2270+ }
2271+ rn = (bit_index + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS;
2272+ assert (rn > 0);
2273+
2274+ rp = MPZ_REALLOC (r, rn);
2275+ un = GMP_ABS (us);
2276+
2277+ mask = GMP_LIMB_MAX >> (rn * GMP_LIMB_BITS - bit_index);
2278+
2279+ if (rn > un)
2280+ {
2281+ /* Quotient (with truncation) is zero, and remainder is
2282+ non-zero */
2283+ if (mode == ((us > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* us != 0 here. */
2284+ {
2285+ /* Have to negate and sign extend. */
2286+ mp_size_t i;
2287+ mp_limb_t cy;
2288+
2289+ for (cy = 1, i = 0; i < un; i++)
2290+ {
2291+ mp_limb_t s = ~u->_mp_d[i] + cy;
2292+ cy = s < cy;
2293+ rp[i] = s;
2294+ }
2295+ assert (cy == 0);
2296+ for (; i < rn - 1; i++)
2297+ rp[i] = GMP_LIMB_MAX;
2298+
2299+ rp[rn-1] = mask;
2300+ us = -us;
2301+ }
2302+ else
2303+ {
2304+ /* Just copy */
2305+ if (r != u)
2306+ mpn_copyi (rp, u->_mp_d, un);
2307+
2308+ rn = un;
2309+ }
2310+ }
2311+ else
2312+ {
2313+ if (r != u)
2314+ mpn_copyi (rp, u->_mp_d, rn - 1);
2315+
2316+ rp[rn-1] = u->_mp_d[rn-1] & mask;
2317+
2318+ if (mode == ((us > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* us != 0 here. */
2319+ {
2320+ /* If r != 0, compute 2^{bit_count} - r. */
2321+ mp_size_t i;
2322+
2323+ for (i = 0; i < rn && rp[i] == 0; i++)
2324+ ;
2325+ if (i < rn)
2326+ {
2327+ /* r > 0, need to flip sign. */
2328+ rp[i] = ~rp[i] + 1;
2329+ while (++i < rn)
2330+ rp[i] = ~rp[i];
2331+
2332+ rp[rn-1] &= mask;
2333+
2334+ /* us is not used for anything else, so we can modify it
2335+ here to indicate flipped sign. */
2336+ us = -us;
2337+ }
2338+ }
2339+ }
2340+ rn = mpn_normalized_size (rp, rn);
2341+ r->_mp_size = us < 0 ? -rn : rn;
2342+}
2343+
2344+void
2345+mpz_cdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2346+{
2347+ mpz_div_q_2exp (r, u, cnt, GMP_DIV_CEIL);
2348+}
2349+
2350+void
2351+mpz_fdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2352+{
2353+ mpz_div_q_2exp (r, u, cnt, GMP_DIV_FLOOR);
2354+}
2355+
2356+void
2357+mpz_tdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2358+{
2359+ mpz_div_q_2exp (r, u, cnt, GMP_DIV_TRUNC);
2360+}
2361+
2362+void
2363+mpz_cdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2364+{
2365+ mpz_div_r_2exp (r, u, cnt, GMP_DIV_CEIL);
2366+}
2367+
2368+void
2369+mpz_fdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2370+{
2371+ mpz_div_r_2exp (r, u, cnt, GMP_DIV_FLOOR);
2372+}
2373+
2374+void
2375+mpz_tdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2376+{
2377+ mpz_div_r_2exp (r, u, cnt, GMP_DIV_TRUNC);
2378+}
2379+
2380+void
2381+mpz_divexact (mpz_t q, const mpz_t n, const mpz_t d)
2382+{
2383+ gmp_assert_nocarry (mpz_div_qr (q, NULL, n, d, GMP_DIV_TRUNC));
2384+}
2385+
2386+int
2387+mpz_divisible_p (const mpz_t n, const mpz_t d)
2388+{
2389+ return mpz_div_qr (NULL, NULL, n, d, GMP_DIV_TRUNC) == 0;
2390+}
2391+
2392+static unsigned long
2393+mpz_div_qr_ui (mpz_t q, mpz_t r,
2394+ const mpz_t n, unsigned long d, enum mpz_div_round_mode mode)
2395+{
2396+ mp_size_t ns, qn;
2397+ mp_ptr qp;
2398+ mp_limb_t rl;
2399+ mp_size_t rs;
2400+
2401+ ns = n->_mp_size;
2402+ if (ns == 0)
2403+ {
2404+ if (q)
2405+ q->_mp_size = 0;
2406+ if (r)
2407+ r->_mp_size = 0;
2408+ return 0;
2409+ }
2410+
2411+ qn = GMP_ABS (ns);
2412+ if (q)
2413+ qp = MPZ_REALLOC (q, qn);
2414+ else
2415+ qp = NULL;
2416+
2417+ rl = mpn_div_qr_1 (qp, n->_mp_d, qn, d);
2418+ assert (rl < d);
2419+
2420+ rs = rl > 0;
2421+ rs = (ns < 0) ? -rs : rs;
2422+
2423+ if (rl > 0 && ( (mode == GMP_DIV_FLOOR && ns < 0)
2424+ || (mode == GMP_DIV_CEIL && ns >= 0)))
2425+ {
2426+ if (q)
2427+ gmp_assert_nocarry (mpn_add_1 (qp, qp, qn, 1));
2428+ rl = d - rl;
2429+ rs = -rs;
2430+ }
2431+
2432+ if (r)
2433+ {
2434+ r->_mp_d[0] = rl;
2435+ r->_mp_size = rs;
2436+ }
2437+ if (q)
2438+ {
2439+ qn -= (qp[qn-1] == 0);
2440+ assert (qn == 0 || qp[qn-1] > 0);
2441+
2442+ q->_mp_size = (ns < 0) ? - qn : qn;
2443+ }
2444+
2445+ return rl;
2446+}
2447+
2448+unsigned long
2449+mpz_cdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d)
2450+{
2451+ return mpz_div_qr_ui (q, r, n, d, GMP_DIV_CEIL);
2452+}
2453+
2454+unsigned long
2455+mpz_fdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d)
2456+{
2457+ return mpz_div_qr_ui (q, r, n, d, GMP_DIV_FLOOR);
2458+}
2459+
2460+unsigned long
2461+mpz_tdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d)
2462+{
2463+ return mpz_div_qr_ui (q, r, n, d, GMP_DIV_TRUNC);
2464+}
2465+
2466+unsigned long
2467+mpz_cdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d)
2468+{
2469+ return mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_CEIL);
2470+}
2471+
2472+unsigned long
2473+mpz_fdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d)
2474+{
2475+ return mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_FLOOR);
2476+}
2477+
2478+unsigned long
2479+mpz_tdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d)
2480+{
2481+ return mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_TRUNC);
2482+}
2483+
2484+unsigned long
2485+mpz_cdiv_r_ui (mpz_t r, const mpz_t n, unsigned long d)
2486+{
2487+ return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_CEIL);
2488+}
2489+unsigned long
2490+mpz_fdiv_r_ui (mpz_t r, const mpz_t n, unsigned long d)
2491+{
2492+ return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_FLOOR);
2493+}
2494+unsigned long
2495+mpz_tdiv_r_ui (mpz_t r, const mpz_t n, unsigned long d)
2496+{
2497+ return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_TRUNC);
2498+}
2499+
2500+unsigned long
2501+mpz_cdiv_ui (const mpz_t n, unsigned long d)
2502+{
2503+ return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_CEIL);
2504+}
2505+
2506+unsigned long
2507+mpz_fdiv_ui (const mpz_t n, unsigned long d)
2508+{
2509+ return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_FLOOR);
2510+}
2511+
2512+unsigned long
2513+mpz_tdiv_ui (const mpz_t n, unsigned long d)
2514+{
2515+ return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_TRUNC);
2516+}
2517+
2518+unsigned long
2519+mpz_mod_ui (mpz_t r, const mpz_t n, unsigned long d)
2520+{
2521+ return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_FLOOR);
2522+}
2523+
2524+void
2525+mpz_divexact_ui (mpz_t q, const mpz_t n, unsigned long d)
2526+{
2527+ gmp_assert_nocarry (mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_TRUNC));
2528+}
2529+
2530+int
2531+mpz_divisible_ui_p (const mpz_t n, unsigned long d)
2532+{
2533+ return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_TRUNC) == 0;
2534+}
2535+
2536+
2537+/* GCD */
2538+static mp_limb_t
2539+mpn_gcd_11 (mp_limb_t u, mp_limb_t v)
2540+{
2541+ unsigned shift;
2542+
2543+ assert ( (u | v) > 0);
2544+
2545+ if (u == 0)
2546+ return v;
2547+ else if (v == 0)
2548+ return u;
2549+
2550+ gmp_ctz (shift, u | v);
2551+
2552+ u >>= shift;
2553+ v >>= shift;
2554+
2555+ if ( (u & 1) == 0)
2556+ MP_LIMB_T_SWAP (u, v);
2557+
2558+ while ( (v & 1) == 0)
2559+ v >>= 1;
2560+
2561+ while (u != v)
2562+ {
2563+ if (u > v)
2564+ {
2565+ u -= v;
2566+ do
2567+ u >>= 1;
2568+ while ( (u & 1) == 0);
2569+ }
2570+ else
2571+ {
2572+ v -= u;
2573+ do
2574+ v >>= 1;
2575+ while ( (v & 1) == 0);
2576+ }
2577+ }
2578+ return u << shift;
2579+}
2580+
2581+unsigned long
2582+mpz_gcd_ui (mpz_t g, const mpz_t u, unsigned long v)
2583+{
2584+ mp_size_t un;
2585+
2586+ if (v == 0)
2587+ {
2588+ if (g)
2589+ mpz_abs (g, u);
2590+ }
2591+ else
2592+ {
2593+ un = GMP_ABS (u->_mp_size);
2594+ if (un != 0)
2595+ v = mpn_gcd_11 (mpn_div_qr_1 (NULL, u->_mp_d, un, v), v);
2596+
2597+ if (g)
2598+ mpz_set_ui (g, v);
2599+ }
2600+
2601+ return v;
2602+}
2603+
2604+static mp_bitcnt_t
2605+mpz_make_odd (mpz_t r)
2606+{
2607+ mp_bitcnt_t shift;
2608+
2609+ assert (r->_mp_size > 0);
2610+ /* Count trailing zeros, equivalent to mpn_scan1, because we know that there is a 1 */
2611+ shift = mpn_common_scan (r->_mp_d[0], 0, r->_mp_d, 0, 0);
2612+ mpz_tdiv_q_2exp (r, r, shift);
2613+
2614+ return shift;
2615+}
2616+
2617+void
2618+mpz_gcd (mpz_t g, const mpz_t u, const mpz_t v)
2619+{
2620+ mpz_t tu, tv;
2621+ mp_bitcnt_t uz, vz, gz;
2622+
2623+ if (u->_mp_size == 0)
2624+ {
2625+ mpz_abs (g, v);
2626+ return;
2627+ }
2628+ if (v->_mp_size == 0)
2629+ {
2630+ mpz_abs (g, u);
2631+ return;
2632+ }
2633+
2634+ mpz_init (tu);
2635+ mpz_init (tv);
2636+
2637+ mpz_abs (tu, u);
2638+ uz = mpz_make_odd (tu);
2639+ mpz_abs (tv, v);
2640+ vz = mpz_make_odd (tv);
2641+ gz = GMP_MIN (uz, vz);
2642+
2643+ if (tu->_mp_size < tv->_mp_size)
2644+ mpz_swap (tu, tv);
2645+
2646+ mpz_tdiv_r (tu, tu, tv);
2647+ if (tu->_mp_size == 0)
2648+ {
2649+ mpz_swap (g, tv);
2650+ }
2651+ else
2652+ for (;;)
2653+ {
2654+ int c;
2655+
2656+ mpz_make_odd (tu);
2657+ c = mpz_cmp (tu, tv);
2658+ if (c == 0)
2659+ {
2660+ mpz_swap (g, tu);
2661+ break;
2662+ }
2663+ if (c < 0)
2664+ mpz_swap (tu, tv);
2665+
2666+ if (tv->_mp_size == 1)
2667+ {
2668+ mp_limb_t vl = tv->_mp_d[0];
2669+ mp_limb_t ul = mpz_tdiv_ui (tu, vl);
2670+ mpz_set_ui (g, mpn_gcd_11 (ul, vl));
2671+ break;
2672+ }
2673+ mpz_sub (tu, tu, tv);
2674+ }
2675+ mpz_clear (tu);
2676+ mpz_clear (tv);
2677+ mpz_mul_2exp (g, g, gz);
2678+}
2679+
2680+void
2681+mpz_gcdext (mpz_t g, mpz_t s, mpz_t t, const mpz_t u, const mpz_t v)
2682+{
2683+ mpz_t tu, tv, s0, s1, t0, t1;
2684+ mp_bitcnt_t uz, vz, gz;
2685+ mp_bitcnt_t power;
2686+
2687+ if (u->_mp_size == 0)
2688+ {
2689+ /* g = 0 u + sgn(v) v */
2690+ signed long sign = mpz_sgn (v);
2691+ mpz_abs (g, v);
2692+ if (s)
2693+ mpz_set_ui (s, 0);
2694+ if (t)
2695+ mpz_set_si (t, sign);
2696+ return;
2697+ }
2698+
2699+ if (v->_mp_size == 0)
2700+ {
2701+ /* g = sgn(u) u + 0 v */
2702+ signed long sign = mpz_sgn (u);
2703+ mpz_abs (g, u);
2704+ if (s)
2705+ mpz_set_si (s, sign);
2706+ if (t)
2707+ mpz_set_ui (t, 0);
2708+ return;
2709+ }
2710+
2711+ mpz_init (tu);
2712+ mpz_init (tv);
2713+ mpz_init (s0);
2714+ mpz_init (s1);
2715+ mpz_init (t0);
2716+ mpz_init (t1);
2717+
2718+ mpz_abs (tu, u);
2719+ uz = mpz_make_odd (tu);
2720+ mpz_abs (tv, v);
2721+ vz = mpz_make_odd (tv);
2722+ gz = GMP_MIN (uz, vz);
2723+
2724+ uz -= gz;
2725+ vz -= gz;
2726+
2727+ /* Cofactors corresponding to odd gcd. gz handled later. */
2728+ if (tu->_mp_size < tv->_mp_size)
2729+ {
2730+ mpz_swap (tu, tv);
2731+ MPZ_SRCPTR_SWAP (u, v);
2732+ MPZ_PTR_SWAP (s, t);
2733+ MP_BITCNT_T_SWAP (uz, vz);
2734+ }
2735+
2736+ /* Maintain
2737+ *
2738+ * u = t0 tu + t1 tv
2739+ * v = s0 tu + s1 tv
2740+ *
2741+ * where u and v denote the inputs with common factors of two
2742+ * eliminated, and det (s0, t0; s1, t1) = 2^p. Then
2743+ *
2744+ * 2^p tu = s1 u - t1 v
2745+ * 2^p tv = -s0 u + t0 v
2746+ */
2747+
2748+ /* After initial division, tu = q tv + tu', we have
2749+ *
2750+ * u = 2^uz (tu' + q tv)
2751+ * v = 2^vz tv
2752+ *
2753+ * or
2754+ *
2755+ * t0 = 2^uz, t1 = 2^uz q
2756+ * s0 = 0, s1 = 2^vz
2757+ */
2758+
2759+ mpz_setbit (t0, uz);
2760+ mpz_tdiv_qr (t1, tu, tu, tv);
2761+ mpz_mul_2exp (t1, t1, uz);
2762+
2763+ mpz_setbit (s1, vz);
2764+ power = uz + vz;
2765+
2766+ if (tu->_mp_size > 0)
2767+ {
2768+ mp_bitcnt_t shift;
2769+ shift = mpz_make_odd (tu);
2770+ mpz_mul_2exp (t0, t0, shift);
2771+ mpz_mul_2exp (s0, s0, shift);
2772+ power += shift;
2773+
2774+ for (;;)
2775+ {
2776+ int c;
2777+ c = mpz_cmp (tu, tv);
2778+ if (c == 0)
2779+ break;
2780+
2781+ if (c < 0)
2782+ {
2783+ /* tv = tv' + tu
2784+ *
2785+ * u = t0 tu + t1 (tv' + tu) = (t0 + t1) tu + t1 tv'
2786+ * v = s0 tu + s1 (tv' + tu) = (s0 + s1) tu + s1 tv' */
2787+
2788+ mpz_sub (tv, tv, tu);
2789+ mpz_add (t0, t0, t1);
2790+ mpz_add (s0, s0, s1);
2791+
2792+ shift = mpz_make_odd (tv);
2793+ mpz_mul_2exp (t1, t1, shift);
2794+ mpz_mul_2exp (s1, s1, shift);
2795+ }
2796+ else
2797+ {
2798+ mpz_sub (tu, tu, tv);
2799+ mpz_add (t1, t0, t1);
2800+ mpz_add (s1, s0, s1);
2801+
2802+ shift = mpz_make_odd (tu);
2803+ mpz_mul_2exp (t0, t0, shift);
2804+ mpz_mul_2exp (s0, s0, shift);
2805+ }
2806+ power += shift;
2807+ }
2808+ }
2809+
2810+ /* Now tv = odd part of gcd, and -s0 and t0 are corresponding
2811+ cofactors. */
2812+
2813+ mpz_mul_2exp (tv, tv, gz);
2814+ mpz_neg (s0, s0);
2815+
2816+ /* 2^p g = s0 u + t0 v. Eliminate one factor of two at a time. To
2817+ adjust cofactors, we need u / g and v / g */
2818+
2819+ mpz_divexact (s1, v, tv);
2820+ mpz_abs (s1, s1);
2821+ mpz_divexact (t1, u, tv);
2822+ mpz_abs (t1, t1);
2823+
2824+ while (power-- > 0)
2825+ {
2826+ /* s0 u + t0 v = (s0 - v/g) u - (t0 + u/g) v */
2827+ if (mpz_odd_p (s0) || mpz_odd_p (t0))
2828+ {
2829+ mpz_sub (s0, s0, s1);
2830+ mpz_add (t0, t0, t1);
2831+ }
2832+ mpz_divexact_ui (s0, s0, 2);
2833+ mpz_divexact_ui (t0, t0, 2);
2834+ }
2835+
2836+ /* Arrange so that |s| < |u| / 2g */
2837+ mpz_add (s1, s0, s1);
2838+ if (mpz_cmpabs (s0, s1) > 0)
2839+ {
2840+ mpz_swap (s0, s1);
2841+ mpz_sub (t0, t0, t1);
2842+ }
2843+ if (u->_mp_size < 0)
2844+ mpz_neg (s0, s0);
2845+ if (v->_mp_size < 0)
2846+ mpz_neg (t0, t0);
2847+
2848+ mpz_swap (g, tv);
2849+ if (s)
2850+ mpz_swap (s, s0);
2851+ if (t)
2852+ mpz_swap (t, t0);
2853+
2854+ mpz_clear (tu);
2855+ mpz_clear (tv);
2856+ mpz_clear (s0);
2857+ mpz_clear (s1);
2858+ mpz_clear (t0);
2859+ mpz_clear (t1);
2860+}
2861+
2862+void
2863+mpz_lcm (mpz_t r, const mpz_t u, const mpz_t v)
2864+{
2865+ mpz_t g;
2866+
2867+ if (u->_mp_size == 0 || v->_mp_size == 0)
2868+ {
2869+ r->_mp_size = 0;
2870+ return;
2871+ }
2872+
2873+ mpz_init (g);
2874+
2875+ mpz_gcd (g, u, v);
2876+ mpz_divexact (g, u, g);
2877+ mpz_mul (r, g, v);
2878+
2879+ mpz_clear (g);
2880+ mpz_abs (r, r);
2881+}
2882+
2883+void
2884+mpz_lcm_ui (mpz_t r, const mpz_t u, unsigned long v)
2885+{
2886+ if (v == 0 || u->_mp_size == 0)
2887+ {
2888+ r->_mp_size = 0;
2889+ return;
2890+ }
2891+
2892+ v /= mpz_gcd_ui (NULL, u, v);
2893+ mpz_mul_ui (r, u, v);
2894+
2895+ mpz_abs (r, r);
2896+}
2897+
2898+int
2899+mpz_invert (mpz_t r, const mpz_t u, const mpz_t m)
2900+{
2901+ mpz_t g, tr;
2902+ int invertible;
2903+
2904+ if (u->_mp_size == 0 || mpz_cmpabs_ui (m, 1) <= 0)
2905+ return 0;
2906+
2907+ mpz_init (g);
2908+ mpz_init (tr);
2909+
2910+ mpz_gcdext (g, tr, NULL, u, m);
2911+ invertible = (mpz_cmp_ui (g, 1) == 0);
2912+
2913+ if (invertible)
2914+ {
2915+ if (tr->_mp_size < 0)
2916+ {
2917+ if (m->_mp_size >= 0)
2918+ mpz_add (tr, tr, m);
2919+ else
2920+ mpz_sub (tr, tr, m);
2921+ }
2922+ mpz_swap (r, tr);
2923+ }
2924+
2925+ mpz_clear (g);
2926+ mpz_clear (tr);
2927+ return invertible;
2928+}
2929+
2930+
2931+/* Higher level operations (sqrt, pow and root) */
2932+
2933+void
2934+mpz_pow_ui (mpz_t r, const mpz_t b, unsigned long e)
2935+{
2936+ unsigned long bit;
2937+ mpz_t tr;
2938+ mpz_init_set_ui (tr, 1);
2939+
2940+ bit = GMP_ULONG_HIGHBIT;
2941+ do
2942+ {
2943+ mpz_mul (tr, tr, tr);
2944+ if (e & bit)
2945+ mpz_mul (tr, tr, b);
2946+ bit >>= 1;
2947+ }
2948+ while (bit > 0);
2949+
2950+ mpz_swap (r, tr);
2951+ mpz_clear (tr);
2952+}
2953+
2954+void
2955+mpz_ui_pow_ui (mpz_t r, unsigned long blimb, unsigned long e)
2956+{
2957+ mpz_t b;
2958+ mpz_init_set_ui (b, blimb);
2959+ mpz_pow_ui (r, b, e);
2960+ mpz_clear (b);
2961+}
2962+
2963+void
2964+mpz_powm (mpz_t r, const mpz_t b, const mpz_t e, const mpz_t m)
2965+{
2966+ mpz_t tr;
2967+ mpz_t base;
2968+ mp_size_t en, mn;
2969+ mp_srcptr mp;
2970+ struct gmp_div_inverse minv;
2971+ unsigned shift;
2972+ mp_ptr tp = NULL;
2973+
2974+ en = GMP_ABS (e->_mp_size);
2975+ mn = GMP_ABS (m->_mp_size);
2976+ if (mn == 0)
2977+ gmp_die ("mpz_powm: Zero modulo.");
2978+
2979+ if (en == 0)
2980+ {
2981+ mpz_set_ui (r, 1);
2982+ return;
2983+ }
2984+
2985+ mp = m->_mp_d;
2986+ mpn_div_qr_invert (&minv, mp, mn);
2987+ shift = minv.shift;
2988+
2989+ if (shift > 0)
2990+ {
2991+ /* To avoid shifts, we do all our reductions, except the final
2992+ one, using a *normalized* m. */
2993+ minv.shift = 0;
2994+
2995+ tp = gmp_xalloc_limbs (mn);
2996+ gmp_assert_nocarry (mpn_lshift (tp, mp, mn, shift));
2997+ mp = tp;
2998+ }
2999+
3000+ mpz_init (base);
3001+
3002+ if (e->_mp_size < 0)
3003+ {
3004+ if (!mpz_invert (base, b, m))
3005+ gmp_die ("mpz_powm: Negative exponent and non-invertible base.");
3006+ }
3007+ else
3008+ {
3009+ mp_size_t bn;
3010+ mpz_abs (base, b);
3011+
3012+ bn = base->_mp_size;
3013+ if (bn >= mn)
3014+ {
3015+ mpn_div_qr_preinv (NULL, base->_mp_d, base->_mp_size, mp, mn, &minv);
3016+ bn = mn;
3017+ }
3018+
3019+ /* We have reduced the absolute value. Now take care of the
3020+ sign. Note that we get zero represented non-canonically as
3021+ m. */
3022+ if (b->_mp_size < 0)
3023+ {
3024+ mp_ptr bp = MPZ_REALLOC (base, mn);
3025+ gmp_assert_nocarry (mpn_sub (bp, mp, mn, bp, bn));
3026+ bn = mn;
3027+ }
3028+ base->_mp_size = mpn_normalized_size (base->_mp_d, bn);
3029+ }
3030+ mpz_init_set_ui (tr, 1);
3031+
3032+ while (en-- > 0)
3033+ {
3034+ mp_limb_t w = e->_mp_d[en];
3035+ mp_limb_t bit;
3036+
3037+ bit = GMP_LIMB_HIGHBIT;
3038+ do
3039+ {
3040+ mpz_mul (tr, tr, tr);
3041+ if (w & bit)
3042+ mpz_mul (tr, tr, base);
3043+ if (tr->_mp_size > mn)
3044+ {
3045+ mpn_div_qr_preinv (NULL, tr->_mp_d, tr->_mp_size, mp, mn, &minv);
3046+ tr->_mp_size = mpn_normalized_size (tr->_mp_d, mn);
3047+ }
3048+ bit >>= 1;
3049+ }
3050+ while (bit > 0);
3051+ }
3052+
3053+ /* Final reduction */
3054+ if (tr->_mp_size >= mn)
3055+ {
3056+ minv.shift = shift;
3057+ mpn_div_qr_preinv (NULL, tr->_mp_d, tr->_mp_size, mp, mn, &minv);
3058+ tr->_mp_size = mpn_normalized_size (tr->_mp_d, mn);
3059+ }
3060+ if (tp)
3061+ gmp_free (tp);
3062+
3063+ mpz_swap (r, tr);
3064+ mpz_clear (tr);
3065+ mpz_clear (base);
3066+}
3067+
3068+void
3069+mpz_powm_ui (mpz_t r, const mpz_t b, unsigned long elimb, const mpz_t m)
3070+{
3071+ mpz_t e;
3072+ mpz_init_set_ui (e, elimb);
3073+ mpz_powm (r, b, e, m);
3074+ mpz_clear (e);
3075+}
3076+
3077+/* x=trunc(y^(1/z)), r=y-x^z */
3078+void
3079+mpz_rootrem (mpz_t x, mpz_t r, const mpz_t y, unsigned long z)
3080+{
3081+ int sgn;
3082+ mpz_t t, u;
3083+
3084+ sgn = y->_mp_size < 0;
3085+ if ((~z & sgn) != 0)
3086+ gmp_die ("mpz_rootrem: Negative argument, with even root.");
3087+ if (z == 0)
3088+ gmp_die ("mpz_rootrem: Zeroth root.");
3089+
3090+ if (mpz_cmpabs_ui (y, 1) <= 0) {
3091+ mpz_set (x, y);
3092+ if (r)
3093+ r->_mp_size = 0;
3094+ return;
3095+ }
3096+
3097+ mpz_init (u);
3098+ {
3099+ mp_bitcnt_t tb;
3100+ tb = mpz_sizeinbase (y, 2) / z + 1;
3101+ mpz_init2 (t, tb);
3102+ mpz_setbit (t, tb);
3103+ }
3104+
3105+ if (z == 2) /* simplify sqrt loop: z-1 == 1 */
3106+ do {
3107+ mpz_swap (u, t); /* u = x */
3108+ mpz_tdiv_q (t, y, u); /* t = y/x */
3109+ mpz_add (t, t, u); /* t = y/x + x */
3110+ mpz_tdiv_q_2exp (t, t, 1); /* x'= (y/x + x)/2 */
3111+ } while (mpz_cmpabs (t, u) < 0); /* |x'| < |x| */
3112+ else /* z != 2 */ {
3113+ mpz_t v;
3114+
3115+ mpz_init (v);
3116+ if (sgn)
3117+ mpz_neg (t, t);
3118+
3119+ do {
3120+ mpz_swap (u, t); /* u = x */
3121+ mpz_pow_ui (t, u, z - 1); /* t = x^(z-1) */
3122+ mpz_tdiv_q (t, y, t); /* t = y/x^(z-1) */
3123+ mpz_mul_ui (v, u, z - 1); /* v = x*(z-1) */
3124+ mpz_add (t, t, v); /* t = y/x^(z-1) + x*(z-1) */
3125+ mpz_tdiv_q_ui (t, t, z); /* x'=(y/x^(z-1) + x*(z-1))/z */
3126+ } while (mpz_cmpabs (t, u) < 0); /* |x'| < |x| */
3127+
3128+ mpz_clear (v);
3129+ }
3130+
3131+ if (r) {
3132+ mpz_pow_ui (t, u, z);
3133+ mpz_sub (r, y, t);
3134+ }
3135+ mpz_swap (x, u);
3136+ mpz_clear (u);
3137+ mpz_clear (t);
3138+}
3139+
3140+int
3141+mpz_root (mpz_t x, const mpz_t y, unsigned long z)
3142+{
3143+ int res;
3144+ mpz_t r;
3145+
3146+ mpz_init (r);
3147+ mpz_rootrem (x, r, y, z);
3148+ res = r->_mp_size == 0;
3149+ mpz_clear (r);
3150+
3151+ return res;
3152+}
3153+
3154+/* Compute s = floor(sqrt(u)) and r = u - s^2. Allows r == NULL */
3155+void
3156+mpz_sqrtrem (mpz_t s, mpz_t r, const mpz_t u)
3157+{
3158+ mpz_rootrem (s, r, u, 2);
3159+}
3160+
3161+void
3162+mpz_sqrt (mpz_t s, const mpz_t u)
3163+{
3164+ mpz_rootrem (s, NULL, u, 2);
3165+}
3166+
3167+
3168+/* Combinatorics */
3169+
3170+void
3171+mpz_fac_ui (mpz_t x, unsigned long n)
3172+{
3173+ mpz_set_ui (x, n + (n == 0));
3174+ for (;n > 2;)
3175+ mpz_mul_ui (x, x, --n);
3176+}
3177+
3178+void
3179+mpz_bin_uiui (mpz_t r, unsigned long n, unsigned long k)
3180+{
3181+ mpz_t t;
3182+
3183+ mpz_set_ui (r, k <= n);
3184+
3185+ if (k > (n >> 1))
3186+ k = (k <= n) ? n - k : 0;
3187+
3188+ mpz_init (t);
3189+ mpz_fac_ui (t, k);
3190+
3191+ for (; k > 0; k--)
3192+ mpz_mul_ui (r, r, n--);
3193+
3194+ mpz_divexact (r, r, t);
3195+ mpz_clear (t);
3196+}
3197+
3198+
3199+/* Logical operations and bit manipulation. */
3200+
3201+/* Numbers are treated as if represented in two's complement (and
3202+ infinitely sign extended). For a negative values we get the two's
3203+ complement from -x = ~x + 1, where ~ is bitwise complement.
3204+ Negation transforms
3205+
3206+ xxxx10...0
3207+
3208+ into
3209+
3210+ yyyy10...0
3211+
3212+ where yyyy is the bitwise complement of xxxx. So least significant
3213+ bits, up to and including the first one bit, are unchanged, and
3214+ the more significant bits are all complemented.
3215+
3216+ To change a bit from zero to one in a negative number, subtract the
3217+ corresponding power of two from the absolute value. This can never
3218+ underflow. To change a bit from one to zero, add the corresponding
3219+ power of two, and this might overflow. E.g., if x = -001111, the
3220+ two's complement is 110001. Clearing the least significant bit, we
3221+ get two's complement 110000, and -010000. */
3222+
3223+int
3224+mpz_tstbit (const mpz_t d, mp_bitcnt_t bit_index)
3225+{
3226+ mp_size_t limb_index;
3227+ unsigned shift;
3228+ mp_size_t ds;
3229+ mp_size_t dn;
3230+ mp_limb_t w;
3231+ int bit;
3232+
3233+ ds = d->_mp_size;
3234+ dn = GMP_ABS (ds);
3235+ limb_index = bit_index / GMP_LIMB_BITS;
3236+ if (limb_index >= dn)
3237+ return ds < 0;
3238+
3239+ shift = bit_index % GMP_LIMB_BITS;
3240+ w = d->_mp_d[limb_index];
3241+ bit = (w >> shift) & 1;
3242+
3243+ if (ds < 0)
3244+ {
3245+ /* d < 0. Check if any of the bits below is set: If so, our bit
3246+ must be complemented. */
3247+ if (shift > 0 && (w << (GMP_LIMB_BITS - shift)) > 0)
3248+ return bit ^ 1;
3249+ while (limb_index-- > 0)
3250+ if (d->_mp_d[limb_index] > 0)
3251+ return bit ^ 1;
3252+ }
3253+ return bit;
3254+}
3255+
3256+static void
3257+mpz_abs_add_bit (mpz_t d, mp_bitcnt_t bit_index)
3258+{
3259+ mp_size_t dn, limb_index;
3260+ mp_limb_t bit;
3261+ mp_ptr dp;
3262+
3263+ dn = GMP_ABS (d->_mp_size);
3264+
3265+ limb_index = bit_index / GMP_LIMB_BITS;
3266+ bit = (mp_limb_t) 1 << (bit_index % GMP_LIMB_BITS);
3267+
3268+ if (limb_index >= dn)
3269+ {
3270+ mp_size_t i;
3271+ /* The bit should be set outside of the end of the number.
3272+ We have to increase the size of the number. */
3273+ dp = MPZ_REALLOC (d, limb_index + 1);
3274+
3275+ dp[limb_index] = bit;
3276+ for (i = dn; i < limb_index; i++)
3277+ dp[i] = 0;
3278+ dn = limb_index + 1;
3279+ }
3280+ else
3281+ {
3282+ mp_limb_t cy;
3283+
3284+ dp = d->_mp_d;
3285+
3286+ cy = mpn_add_1 (dp + limb_index, dp + limb_index, dn - limb_index, bit);
3287+ if (cy > 0)
3288+ {
3289+ dp = MPZ_REALLOC (d, dn + 1);
3290+ dp[dn++] = cy;
3291+ }
3292+ }
3293+
3294+ d->_mp_size = (d->_mp_size < 0) ? - dn : dn;
3295+}
3296+
3297+static void
3298+mpz_abs_sub_bit (mpz_t d, mp_bitcnt_t bit_index)
3299+{
3300+ mp_size_t dn, limb_index;
3301+ mp_ptr dp;
3302+ mp_limb_t bit;
3303+
3304+ dn = GMP_ABS (d->_mp_size);
3305+ dp = d->_mp_d;
3306+
3307+ limb_index = bit_index / GMP_LIMB_BITS;
3308+ bit = (mp_limb_t) 1 << (bit_index % GMP_LIMB_BITS);
3309+
3310+ assert (limb_index < dn);
3311+
3312+ gmp_assert_nocarry (mpn_sub_1 (dp + limb_index, dp + limb_index,
3313+ dn - limb_index, bit));
3314+ dn -= (dp[dn-1] == 0);
3315+ d->_mp_size = (d->_mp_size < 0) ? - dn : dn;
3316+}
3317+
3318+void
3319+mpz_setbit (mpz_t d, mp_bitcnt_t bit_index)
3320+{
3321+ if (!mpz_tstbit (d, bit_index))
3322+ {
3323+ if (d->_mp_size >= 0)
3324+ mpz_abs_add_bit (d, bit_index);
3325+ else
3326+ mpz_abs_sub_bit (d, bit_index);
3327+ }
3328+}
3329+
3330+void
3331+mpz_clrbit (mpz_t d, mp_bitcnt_t bit_index)
3332+{
3333+ if (mpz_tstbit (d, bit_index))
3334+ {
3335+ if (d->_mp_size >= 0)
3336+ mpz_abs_sub_bit (d, bit_index);
3337+ else
3338+ mpz_abs_add_bit (d, bit_index);
3339+ }
3340+}
3341+
3342+void
3343+mpz_combit (mpz_t d, mp_bitcnt_t bit_index)
3344+{
3345+ if (mpz_tstbit (d, bit_index) ^ (d->_mp_size < 0))
3346+ mpz_abs_sub_bit (d, bit_index);
3347+ else
3348+ mpz_abs_add_bit (d, bit_index);
3349+}
3350+
3351+void
3352+mpz_com (mpz_t r, const mpz_t u)
3353+{
3354+ mpz_neg (r, u);
3355+ mpz_sub_ui (r, r, 1);
3356+}
3357+
3358+void
3359+mpz_and (mpz_t r, const mpz_t u, const mpz_t v)
3360+{
3361+ mp_size_t un, vn, rn, i;
3362+ mp_ptr up, vp, rp;
3363+
3364+ mp_limb_t ux, vx, rx;
3365+ mp_limb_t uc, vc, rc;
3366+ mp_limb_t ul, vl, rl;
3367+
3368+ un = GMP_ABS (u->_mp_size);
3369+ vn = GMP_ABS (v->_mp_size);
3370+ if (un < vn)
3371+ {
3372+ MPZ_SRCPTR_SWAP (u, v);
3373+ MP_SIZE_T_SWAP (un, vn);
3374+ }
3375+ if (vn == 0)
3376+ {
3377+ r->_mp_size = 0;
3378+ return;
3379+ }
3380+
3381+ uc = u->_mp_size < 0;
3382+ vc = v->_mp_size < 0;
3383+ rc = uc & vc;
3384+
3385+ ux = -uc;
3386+ vx = -vc;
3387+ rx = -rc;
3388+
3389+ /* If the smaller input is positive, higher limbs don't matter. */
3390+ rn = vx ? un : vn;
3391+
3392+ rp = MPZ_REALLOC (r, rn + rc);
3393+
3394+ up = u->_mp_d;
3395+ vp = v->_mp_d;
3396+
3397+ i = 0;
3398+ do
3399+ {
3400+ ul = (up[i] ^ ux) + uc;
3401+ uc = ul < uc;
3402+
3403+ vl = (vp[i] ^ vx) + vc;
3404+ vc = vl < vc;
3405+
3406+ rl = ( (ul & vl) ^ rx) + rc;
3407+ rc = rl < rc;
3408+ rp[i] = rl;
3409+ }
3410+ while (++i < vn);
3411+ assert (vc == 0);
3412+
3413+ for (; i < rn; i++)
3414+ {
3415+ ul = (up[i] ^ ux) + uc;
3416+ uc = ul < uc;
3417+
3418+ rl = ( (ul & vx) ^ rx) + rc;
3419+ rc = rl < rc;
3420+ rp[i] = rl;
3421+ }
3422+ if (rc)
3423+ rp[rn++] = rc;
3424+ else
3425+ rn = mpn_normalized_size (rp, rn);
3426+
3427+ r->_mp_size = rx ? -rn : rn;
3428+}
3429+
3430+void
3431+mpz_ior (mpz_t r, const mpz_t u, const mpz_t v)
3432+{
3433+ mp_size_t un, vn, rn, i;
3434+ mp_ptr up, vp, rp;
3435+
3436+ mp_limb_t ux, vx, rx;
3437+ mp_limb_t uc, vc, rc;
3438+ mp_limb_t ul, vl, rl;
3439+
3440+ un = GMP_ABS (u->_mp_size);
3441+ vn = GMP_ABS (v->_mp_size);
3442+ if (un < vn)
3443+ {
3444+ MPZ_SRCPTR_SWAP (u, v);
3445+ MP_SIZE_T_SWAP (un, vn);
3446+ }
3447+ if (vn == 0)
3448+ {
3449+ mpz_set (r, u);
3450+ return;
3451+ }
3452+
3453+ uc = u->_mp_size < 0;
3454+ vc = v->_mp_size < 0;
3455+ rc = uc | vc;
3456+
3457+ ux = -uc;
3458+ vx = -vc;
3459+ rx = -rc;
3460+
3461+ /* If the smaller input is negative, by sign extension higher limbs
3462+ don't matter. */
3463+ rn = vx ? vn : un;
3464+
3465+ rp = MPZ_REALLOC (r, rn + rc);
3466+
3467+ up = u->_mp_d;
3468+ vp = v->_mp_d;
3469+
3470+ i = 0;
3471+ do
3472+ {
3473+ ul = (up[i] ^ ux) + uc;
3474+ uc = ul < uc;
3475+
3476+ vl = (vp[i] ^ vx) + vc;
3477+ vc = vl < vc;
3478+
3479+ rl = ( (ul | vl) ^ rx) + rc;
3480+ rc = rl < rc;
3481+ rp[i] = rl;
3482+ }
3483+ while (++i < vn);
3484+ assert (vc == 0);
3485+
3486+ for (; i < rn; i++)
3487+ {
3488+ ul = (up[i] ^ ux) + uc;
3489+ uc = ul < uc;
3490+
3491+ rl = ( (ul | vx) ^ rx) + rc;
3492+ rc = rl < rc;
3493+ rp[i] = rl;
3494+ }
3495+ if (rc)
3496+ rp[rn++] = rc;
3497+ else
3498+ rn = mpn_normalized_size (rp, rn);
3499+
3500+ r->_mp_size = rx ? -rn : rn;
3501+}
3502+
3503+void
3504+mpz_xor (mpz_t r, const mpz_t u, const mpz_t v)
3505+{
3506+ mp_size_t un, vn, i;
3507+ mp_ptr up, vp, rp;
3508+
3509+ mp_limb_t ux, vx, rx;
3510+ mp_limb_t uc, vc, rc;
3511+ mp_limb_t ul, vl, rl;
3512+
3513+ un = GMP_ABS (u->_mp_size);
3514+ vn = GMP_ABS (v->_mp_size);
3515+ if (un < vn)
3516+ {
3517+ MPZ_SRCPTR_SWAP (u, v);
3518+ MP_SIZE_T_SWAP (un, vn);
3519+ }
3520+ if (vn == 0)
3521+ {
3522+ mpz_set (r, u);
3523+ return;
3524+ }
3525+
3526+ uc = u->_mp_size < 0;
3527+ vc = v->_mp_size < 0;
3528+ rc = uc ^ vc;
3529+
3530+ ux = -uc;
3531+ vx = -vc;
3532+ rx = -rc;
3533+
3534+ rp = MPZ_REALLOC (r, un + rc);
3535+
3536+ up = u->_mp_d;
3537+ vp = v->_mp_d;
3538+
3539+ i = 0;
3540+ do
3541+ {
3542+ ul = (up[i] ^ ux) + uc;
3543+ uc = ul < uc;
3544+
3545+ vl = (vp[i] ^ vx) + vc;
3546+ vc = vl < vc;
3547+
3548+ rl = (ul ^ vl ^ rx) + rc;
3549+ rc = rl < rc;
3550+ rp[i] = rl;
3551+ }
3552+ while (++i < vn);
3553+ assert (vc == 0);
3554+
3555+ for (; i < un; i++)
3556+ {
3557+ ul = (up[i] ^ ux) + uc;
3558+ uc = ul < uc;
3559+
3560+ rl = (ul ^ ux) + rc;
3561+ rc = rl < rc;
3562+ rp[i] = rl;
3563+ }
3564+ if (rc)
3565+ rp[un++] = rc;
3566+ else
3567+ un = mpn_normalized_size (rp, un);
3568+
3569+ r->_mp_size = rx ? -un : un;
3570+}
3571+
3572+static unsigned
3573+gmp_popcount_limb (mp_limb_t x)
3574+{
3575+ unsigned c;
3576+
3577+ /* Do 16 bits at a time, to avoid limb-sized constants. */
3578+ for (c = 0; x > 0; x >>= 16)
3579+ {
3580+ unsigned w = ((x >> 1) & 0x5555) + (x & 0x5555);
3581+ w = ((w >> 2) & 0x3333) + (w & 0x3333);
3582+ w = ((w >> 4) & 0x0f0f) + (w & 0x0f0f);
3583+ w = (w >> 8) + (w & 0x00ff);
3584+ c += w;
3585+ }
3586+ return c;
3587+}
3588+
3589+mp_bitcnt_t
3590+mpz_popcount (const mpz_t u)
3591+{
3592+ mp_size_t un, i;
3593+ mp_bitcnt_t c;
3594+
3595+ un = u->_mp_size;
3596+
3597+ if (un < 0)
3598+ return ~(mp_bitcnt_t) 0;
3599+
3600+ for (c = 0, i = 0; i < un; i++)
3601+ c += gmp_popcount_limb (u->_mp_d[i]);
3602+
3603+ return c;
3604+}
3605+
3606+mp_bitcnt_t
3607+mpz_hamdist (const mpz_t u, const mpz_t v)
3608+{
3609+ mp_size_t un, vn, i;
3610+ mp_limb_t uc, vc, ul, vl, comp;
3611+ mp_srcptr up, vp;
3612+ mp_bitcnt_t c;
3613+
3614+ un = u->_mp_size;
3615+ vn = v->_mp_size;
3616+
3617+ if ( (un ^ vn) < 0)
3618+ return ~(mp_bitcnt_t) 0;
3619+
3620+ comp = - (uc = vc = (un < 0));
3621+ if (uc)
3622+ {
3623+ assert (vn < 0);
3624+ un = -un;
3625+ vn = -vn;
3626+ }
3627+
3628+ up = u->_mp_d;
3629+ vp = v->_mp_d;
3630+
3631+ if (un < vn)
3632+ MPN_SRCPTR_SWAP (up, un, vp, vn);
3633+
3634+ for (i = 0, c = 0; i < vn; i++)
3635+ {
3636+ ul = (up[i] ^ comp) + uc;
3637+ uc = ul < uc;
3638+
3639+ vl = (vp[i] ^ comp) + vc;
3640+ vc = vl < vc;
3641+
3642+ c += gmp_popcount_limb (ul ^ vl);
3643+ }
3644+ assert (vc == 0);
3645+
3646+ for (; i < un; i++)
3647+ {
3648+ ul = (up[i] ^ comp) + uc;
3649+ uc = ul < uc;
3650+
3651+ c += gmp_popcount_limb (ul ^ comp);
3652+ }
3653+
3654+ return c;
3655+}
3656+
3657+mp_bitcnt_t
3658+mpz_scan1 (const mpz_t u, mp_bitcnt_t starting_bit)
3659+{
3660+ mp_ptr up;
3661+ mp_size_t us, un, i;
3662+ mp_limb_t limb, ux;
3663+
3664+ us = u->_mp_size;
3665+ un = GMP_ABS (us);
3666+ i = starting_bit / GMP_LIMB_BITS;
3667+
3668+ /* Past the end there's no 1 bits for u>=0, or an immediate 1 bit
3669+ for u<0. Notice this test picks up any u==0 too. */
3670+ if (i >= un)
3671+ return (us >= 0 ? ~(mp_bitcnt_t) 0 : starting_bit);
3672+
3673+ up = u->_mp_d;
3674+ ux = 0;
3675+ limb = up[i];
3676+
3677+ if (starting_bit != 0)
3678+ {
3679+ if (us < 0)
3680+ {
3681+ ux = mpn_zero_p (up, i);
3682+ limb = ~ limb + ux;
3683+ ux = - (mp_limb_t) (limb >= ux);
3684+ }
3685+
3686+ /* Mask to 0 all bits before starting_bit, thus ignoring them. */
3687+ limb &= (GMP_LIMB_MAX << (starting_bit % GMP_LIMB_BITS));
3688+ }
3689+
3690+ return mpn_common_scan (limb, i, up, un, ux);
3691+}
3692+
3693+mp_bitcnt_t
3694+mpz_scan0 (const mpz_t u, mp_bitcnt_t starting_bit)
3695+{
3696+ mp_ptr up;
3697+ mp_size_t us, un, i;
3698+ mp_limb_t limb, ux;
3699+
3700+ us = u->_mp_size;
3701+ ux = - (mp_limb_t) (us >= 0);
3702+ un = GMP_ABS (us);
3703+ i = starting_bit / GMP_LIMB_BITS;
3704+
3705+ /* When past end, there's an immediate 0 bit for u>=0, or no 0 bits for
3706+ u<0. Notice this test picks up all cases of u==0 too. */
3707+ if (i >= un)
3708+ return (ux ? starting_bit : ~(mp_bitcnt_t) 0);
3709+
3710+ up = u->_mp_d;
3711+ limb = up[i] ^ ux;
3712+
3713+ if (ux == 0)
3714+ limb -= mpn_zero_p (up, i); /* limb = ~(~limb + zero_p) */
3715+
3716+ /* Mask all bits before starting_bit, thus ignoring them. */
3717+ limb &= (GMP_LIMB_MAX << (starting_bit % GMP_LIMB_BITS));
3718+
3719+ return mpn_common_scan (limb, i, up, un, ux);
3720+}
3721+
3722+
3723+/* MPZ base conversion. */
3724+
3725+size_t
3726+mpz_sizeinbase (const mpz_t u, int base)
3727+{
3728+ mp_size_t un;
3729+ mp_srcptr up;
3730+ mp_ptr tp;
3731+ mp_bitcnt_t bits;
3732+ struct gmp_div_inverse bi;
3733+ size_t ndigits;
3734+
3735+ assert (base >= 2);
3736+ assert (base <= 36);
3737+
3738+ un = GMP_ABS (u->_mp_size);
3739+ if (un == 0)
3740+ return 1;
3741+
3742+ up = u->_mp_d;
3743+
3744+ bits = (un - 1) * GMP_LIMB_BITS + mpn_limb_size_in_base_2 (up[un-1]);
3745+ switch (base)
3746+ {
3747+ case 2:
3748+ return bits;
3749+ case 4:
3750+ return (bits + 1) / 2;
3751+ case 8:
3752+ return (bits + 2) / 3;
3753+ case 16:
3754+ return (bits + 3) / 4;
3755+ case 32:
3756+ return (bits + 4) / 5;
3757+ /* FIXME: Do something more clever for the common case of base
3758+ 10. */
3759+ }
3760+
3761+ tp = gmp_xalloc_limbs (un);
3762+ mpn_copyi (tp, up, un);
3763+ mpn_div_qr_1_invert (&bi, base);
3764+
3765+ ndigits = 0;
3766+ do
3767+ {
3768+ ndigits++;
3769+ mpn_div_qr_1_preinv (tp, tp, un, &bi);
3770+ un -= (tp[un-1] == 0);
3771+ }
3772+ while (un > 0);
3773+
3774+ gmp_free (tp);
3775+ return ndigits;
3776+}
3777+
3778+char *
3779+mpz_get_str (char *sp, int base, const mpz_t u)
3780+{
3781+ unsigned bits;
3782+ const char *digits;
3783+ mp_size_t un;
3784+ size_t i, sn;
3785+
3786+ if (base >= 0)
3787+ {
3788+ digits = "0123456789abcdefghijklmnopqrstuvwxyz";
3789+ }
3790+ else
3791+ {
3792+ base = -base;
3793+ digits = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
3794+ }
3795+ if (base <= 1)
3796+ base = 10;
3797+ if (base > 36)
3798+ return NULL;
3799+
3800+ sn = 1 + mpz_sizeinbase (u, base);
3801+ if (!sp)
3802+ sp = gmp_xalloc (1 + sn);
3803+
3804+ un = GMP_ABS (u->_mp_size);
3805+
3806+ if (un == 0)
3807+ {
3808+ sp[0] = '0';
3809+ sp[1] = '\0';
3810+ return sp;
3811+ }
3812+
3813+ i = 0;
3814+
3815+ if (u->_mp_size < 0)
3816+ sp[i++] = '-';
3817+
3818+ bits = mpn_base_power_of_two_p (base);
3819+
3820+ if (bits)
3821+ /* Not modified in this case. */
3822+ sn = i + mpn_get_str_bits ((unsigned char *) sp + i, bits, u->_mp_d, un);
3823+ else
3824+ {
3825+ struct mpn_base_info info;
3826+ mp_ptr tp;
3827+
3828+ mpn_get_base_info (&info, base);
3829+ tp = gmp_xalloc_limbs (un);
3830+ mpn_copyi (tp, u->_mp_d, un);
3831+
3832+ sn = i + mpn_get_str_other ((unsigned char *) sp + i, base, &info, tp, un);
3833+ gmp_free (tp);
3834+ }
3835+
3836+ for (; i < sn; i++)
3837+ sp[i] = digits[(unsigned char) sp[i]];
3838+
3839+ sp[sn] = '\0';
3840+ return sp;
3841+}
3842+
3843+int
3844+mpz_set_str (mpz_t r, const char *sp, int base)
3845+{
3846+ unsigned bits;
3847+ mp_size_t rn, alloc;
3848+ mp_ptr rp;
3849+ size_t sn;
3850+ int sign;
3851+ unsigned char *dp;
3852+
3853+ assert (base == 0 || (base >= 2 && base <= 36));
3854+
3855+ while (isspace( (unsigned char) *sp))
3856+ sp++;
3857+
3858+ sign = (*sp == '-');
3859+ sp += sign;
3860+
3861+ if (base == 0)
3862+ {
3863+ if (*sp == '0')
3864+ {
3865+ sp++;
3866+ if (*sp == 'x' || *sp == 'X')
3867+ {
3868+ base = 16;
3869+ sp++;
3870+ }
3871+ else if (*sp == 'b' || *sp == 'B')
3872+ {
3873+ base = 2;
3874+ sp++;
3875+ }
3876+ else
3877+ base = 8;
3878+ }
3879+ else
3880+ base = 10;
3881+ }
3882+
3883+ sn = strlen (sp);
3884+ dp = gmp_xalloc (sn + (sn == 0));
3885+
3886+ for (sn = 0; *sp; sp++)
3887+ {
3888+ unsigned digit;
3889+
3890+ if (isspace ((unsigned char) *sp))
3891+ continue;
3892+ if (*sp >= '0' && *sp <= '9')
3893+ digit = *sp - '0';
3894+ else if (*sp >= 'a' && *sp <= 'z')
3895+ digit = *sp - 'a' + 10;
3896+ else if (*sp >= 'A' && *sp <= 'Z')
3897+ digit = *sp - 'A' + 10;
3898+ else
3899+ digit = base; /* fail */
3900+
3901+ if (digit >= base)
3902+ {
3903+ gmp_free (dp);
3904+ r->_mp_size = 0;
3905+ return -1;
3906+ }
3907+
3908+ dp[sn++] = digit;
3909+ }
3910+
3911+ bits = mpn_base_power_of_two_p (base);
3912+
3913+ if (bits > 0)
3914+ {
3915+ alloc = (sn * bits + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS;
3916+ rp = MPZ_REALLOC (r, alloc);
3917+ rn = mpn_set_str_bits (rp, dp, sn, bits);
3918+ }
3919+ else
3920+ {
3921+ struct mpn_base_info info;
3922+ mpn_get_base_info (&info, base);
3923+ alloc = (sn + info.exp - 1) / info.exp;
3924+ rp = MPZ_REALLOC (r, alloc);
3925+ rn = mpn_set_str_other (rp, dp, sn, base, &info);
3926+ }
3927+ assert (rn <= alloc);
3928+ gmp_free (dp);
3929+
3930+ r->_mp_size = sign ? - rn : rn;
3931+
3932+ return 0;
3933+}
3934+
3935+int
3936+mpz_init_set_str (mpz_t r, const char *sp, int base)
3937+{
3938+ mpz_init (r);
3939+ return mpz_set_str (r, sp, base);
3940+}
3941+
3942+size_t
3943+mpz_out_str (FILE *stream, int base, const mpz_t x)
3944+{
3945+ char *str;
3946+ size_t len;
3947+
3948+ str = mpz_get_str (NULL, base, x);
3949+ len = strlen (str);
3950+ len = fwrite (str, 1, len, stream);
3951+ gmp_free (str);
3952+ return len;
3953+}
3954+
3955+
3956+static int
3957+gmp_detect_endian (void)
3958+{
3959+ static const int i = 2;
3960+ const unsigned char *p = (const unsigned char *) &i;
3961+ return 1 - *p;
3962+}
3963+
3964+/* Import and export. Does not support nails. */
3965+void
3966+mpz_import (mpz_t r, size_t count, int order, size_t size, int endian,
3967+ size_t nails, const void *src)
3968+{
3969+ const unsigned char *p;
3970+ ptrdiff_t word_step;
3971+ mp_ptr rp;
3972+ mp_size_t rn;
3973+
3974+ /* The current (partial) limb. */
3975+ mp_limb_t limb;
3976+ /* The number of bytes already copied to this limb (starting from
3977+ the low end). */
3978+ size_t bytes;
3979+ /* The index where the limb should be stored, when completed. */
3980+ mp_size_t i;
3981+
3982+ if (nails != 0)
3983+ gmp_die ("mpz_import: Nails not supported.");
3984+
3985+ assert (order == 1 || order == -1);
3986+ assert (endian >= -1 && endian <= 1);
3987+
3988+ if (endian == 0)
3989+ endian = gmp_detect_endian ();
3990+
3991+ p = (unsigned char *) src;
3992+
3993+ word_step = (order != endian) ? 2 * size : 0;
3994+
3995+ /* Process bytes from the least significant end, so point p at the
3996+ least significant word. */
3997+ if (order == 1)
3998+ {
3999+ p += size * (count - 1);
4000+ word_step = - word_step;
4001+ }
4002+
4003+ /* And at least significant byte of that word. */
4004+ if (endian == 1)
4005+ p += (size - 1);
4006+
4007+ rn = (size * count + sizeof(mp_limb_t) - 1) / sizeof(mp_limb_t);
4008+ rp = MPZ_REALLOC (r, rn);
4009+
4010+ for (limb = 0, bytes = 0, i = 0; count > 0; count--, p += word_step)
4011+ {
4012+ size_t j;
4013+ for (j = 0; j < size; j++, p -= (ptrdiff_t) endian)
4014+ {
4015+ limb |= (mp_limb_t) *p << (bytes++ * CHAR_BIT);
4016+ if (bytes == sizeof(mp_limb_t))
4017+ {
4018+ rp[i++] = limb;
4019+ bytes = 0;
4020+ limb = 0;
4021+ }
4022+ }
4023+ }
4024+ assert (i + (bytes > 0) == rn);
4025+ if (limb != 0)
4026+ rp[i++] = limb;
4027+ else
4028+ i = mpn_normalized_size (rp, i);
4029+
4030+ r->_mp_size = i;
4031+}
4032+
4033+void *
4034+mpz_export (void *r, size_t *countp, int order, size_t size, int endian,
4035+ size_t nails, const mpz_t u)
4036+{
4037+ size_t count;
4038+ mp_size_t un;
4039+
4040+ if (nails != 0)
4041+ gmp_die ("mpz_import: Nails not supported.");
4042+
4043+ assert (order == 1 || order == -1);
4044+ assert (endian >= -1 && endian <= 1);
4045+ assert (size > 0 || u->_mp_size == 0);
4046+
4047+ un = u->_mp_size;
4048+ count = 0;
4049+ if (un != 0)
4050+ {
4051+ size_t k;
4052+ unsigned char *p;
4053+ ptrdiff_t word_step;
4054+ /* The current (partial) limb. */
4055+ mp_limb_t limb;
4056+ /* The number of bytes left to to in this limb. */
4057+ size_t bytes;
4058+ /* The index where the limb was read. */
4059+ mp_size_t i;
4060+
4061+ un = GMP_ABS (un);
4062+
4063+ /* Count bytes in top limb. */
4064+ limb = u->_mp_d[un-1];
4065+ assert (limb != 0);
4066+
4067+ k = 0;
4068+ do {
4069+ k++; limb >>= CHAR_BIT;
4070+ } while (limb != 0);
4071+
4072+ count = (k + (un-1) * sizeof (mp_limb_t) + size - 1) / size;
4073+
4074+ if (!r)
4075+ r = gmp_xalloc (count * size);
4076+
4077+ if (endian == 0)
4078+ endian = gmp_detect_endian ();
4079+
4080+ p = (unsigned char *) r;
4081+
4082+ word_step = (order != endian) ? 2 * size : 0;
4083+
4084+ /* Process bytes from the least significant end, so point p at the
4085+ least significant word. */
4086+ if (order == 1)
4087+ {
4088+ p += size * (count - 1);
4089+ word_step = - word_step;
4090+ }
4091+
4092+ /* And at least significant byte of that word. */
4093+ if (endian == 1)
4094+ p += (size - 1);
4095+
4096+ for (bytes = 0, i = 0, k = 0; k < count; k++, p += word_step)
4097+ {
4098+ size_t j;
4099+ for (j = 0; j < size; j++, p -= (ptrdiff_t) endian)
4100+ {
4101+ if (bytes == 0)
4102+ {
4103+ if (i < un)
4104+ limb = u->_mp_d[i++];
4105+ bytes = sizeof (mp_limb_t);
4106+ }
4107+ *p = limb;
4108+ limb >>= CHAR_BIT;
4109+ bytes--;
4110+ }
4111+ }
4112+ assert (i == un);
4113+ assert (k == count);
4114+ }
4115+
4116+ if (countp)
4117+ *countp = count;
4118+
4119+ return r;
4120+}