1 /* mini-gmp, a minimalistic implementation of a GNU GMP subset.
3 Contributed to the GNU project by Niels Möller
5 Copyright 1991-1997, 1999-2017 Free Software Foundation, Inc.
7 This file is part of the GNU MP Library.
9 The GNU MP Library is free software; you can redistribute it and/or modify
10 it under the terms of either:
12 * the GNU Lesser General Public License as published by the Free
13 Software Foundation; either version 3 of the License, or (at your
14 option) any later version.
18 * the GNU General Public License as published by the Free Software
19 Foundation; either version 2 of the License, or (at your option) any
22 or both in parallel, as here.
24 The GNU MP Library is distributed in the hope that it will be useful, but
25 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
26 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
29 You should have received copies of the GNU General Public License and the
30 GNU Lesser General Public License along with the GNU MP Library. If not,
31 see https://www.gnu.org/licenses/. */
33 /* NOTE: All functions in this file which are not declared in
34 mini-gmp.h are internal, and are not intended to be compatible
35 neither with GMP nor with future versions of mini-gmp. */
37 /* Much of the material copied from GMP files, including: gmp-impl.h,
38 longlong.h, mpn/generic/add_n.c, mpn/generic/addmul_1.c,
39 mpn/generic/lshift.c, mpn/generic/mul_1.c,
40 mpn/generic/mul_basecase.c, mpn/generic/rshift.c,
41 mpn/generic/sbpi1_div_qr.c, mpn/generic/sub_n.c,
42 mpn/generic/submul_1.c. */
55 #define GMP_LIMB_BITS (sizeof(mp_limb_t) * CHAR_BIT)
57 #define GMP_LIMB_MAX (~ (mp_limb_t) 0)
58 #define GMP_LIMB_HIGHBIT ((mp_limb_t) 1 << (GMP_LIMB_BITS - 1))
60 #define GMP_HLIMB_BIT ((mp_limb_t) 1 << (GMP_LIMB_BITS / 2))
61 #define GMP_LLIMB_MASK (GMP_HLIMB_BIT - 1)
63 #define GMP_ULONG_BITS (sizeof(unsigned long) * CHAR_BIT)
64 #define GMP_ULONG_HIGHBIT ((unsigned long) 1 << (GMP_ULONG_BITS - 1))
66 #define GMP_ABS(x) ((x) >= 0 ? (x) : -(x))
67 #define GMP_NEG_CAST(T,x) (-((T)((x) + 1) - 1))
69 #define GMP_MIN(a, b) ((a) < (b) ? (a) : (b))
70 #define GMP_MAX(a, b) ((a) > (b) ? (a) : (b))
72 #define GMP_CMP(a,b) (((a) > (b)) - ((a) < (b)))
74 /* Return non-zero if xp,xsize and yp,ysize overlap.
75 If xp+xsize<=yp there's no overlap, or if yp+ysize<=xp there's no
76 overlap. If both these are false, there's an overlap. */
77 #define GMP_MPN_OVERLAP_P(xp, xsize, yp, ysize) \
78 ((xp) + (xsize) > (yp) && (yp) + (ysize) > (xp))
80 #define gmp_assert_nocarry(x) do { \
81 mp_limb_t __cy = (x); \
85 #define gmp_clz(count, x) do { \
86 mp_limb_t __clz_x = (x); \
89 (__clz_x & ((mp_limb_t) 0xff << (GMP_LIMB_BITS - 8))) == 0; \
92 for (; (__clz_x & GMP_LIMB_HIGHBIT) == 0; __clz_c++) \
97 #define gmp_ctz(count, x) do { \
98 mp_limb_t __ctz_x = (x); \
99 unsigned __ctz_c = 0; \
100 gmp_clz (__ctz_c, __ctz_x & - __ctz_x); \
101 (count) = GMP_LIMB_BITS - 1 - __ctz_c; \
104 #define gmp_add_ssaaaa(sh, sl, ah, al, bh, bl) \
108 (sh) = (ah) + (bh) + (__x < (al)); \
112 #define gmp_sub_ddmmss(sh, sl, ah, al, bh, bl) \
116 (sh) = (ah) - (bh) - ((al) < (bl)); \
120 #define gmp_umul_ppmm(w1, w0, u, v) \
122 mp_limb_t __x0, __x1, __x2, __x3; \
123 unsigned __ul, __vl, __uh, __vh; \
124 mp_limb_t __u = (u), __v = (v); \
126 __ul = __u & GMP_LLIMB_MASK; \
127 __uh = __u >> (GMP_LIMB_BITS / 2); \
128 __vl = __v & GMP_LLIMB_MASK; \
129 __vh = __v >> (GMP_LIMB_BITS / 2); \
131 __x0 = (mp_limb_t) __ul * __vl; \
132 __x1 = (mp_limb_t) __ul * __vh; \
133 __x2 = (mp_limb_t) __uh * __vl; \
134 __x3 = (mp_limb_t) __uh * __vh; \
136 __x1 += __x0 >> (GMP_LIMB_BITS / 2);/* this can't give carry */ \
137 __x1 += __x2; /* but this indeed can */ \
138 if (__x1 < __x2) /* did we get it? */ \
139 __x3 += GMP_HLIMB_BIT; /* yes, add it in the proper pos. */ \
141 (w1) = __x3 + (__x1 >> (GMP_LIMB_BITS / 2)); \
142 (w0) = (__x1 << (GMP_LIMB_BITS / 2)) + (__x0 & GMP_LLIMB_MASK); \
145 #define gmp_udiv_qrnnd_preinv(q, r, nh, nl, d, di) \
147 mp_limb_t _qh, _ql, _r, _mask; \
148 gmp_umul_ppmm (_qh, _ql, (nh), (di)); \
149 gmp_add_ssaaaa (_qh, _ql, _qh, _ql, (nh) + 1, (nl)); \
150 _r = (nl) - _qh * (d); \
151 _mask = -(mp_limb_t) (_r > _ql); /* both > and >= are OK */ \
164 #define gmp_udiv_qr_3by2(q, r1, r0, n2, n1, n0, d1, d0, dinv) \
166 mp_limb_t _q0, _t1, _t0, _mask; \
167 gmp_umul_ppmm ((q), _q0, (n2), (dinv)); \
168 gmp_add_ssaaaa ((q), _q0, (q), _q0, (n2), (n1)); \
170 /* Compute the two most significant limbs of n - q'd */ \
171 (r1) = (n1) - (d1) * (q); \
172 gmp_sub_ddmmss ((r1), (r0), (r1), (n0), (d1), (d0)); \
173 gmp_umul_ppmm (_t1, _t0, (d0), (q)); \
174 gmp_sub_ddmmss ((r1), (r0), (r1), (r0), _t1, _t0); \
177 /* Conditionally adjust q and the remainders */ \
178 _mask = - (mp_limb_t) ((r1) >= _q0); \
180 gmp_add_ssaaaa ((r1), (r0), (r1), (r0), _mask & (d1), _mask & (d0)); \
183 if ((r1) > (d1) || (r0) >= (d0)) \
186 gmp_sub_ddmmss ((r1), (r0), (r1), (r0), (d1), (d0)); \
192 #define MP_LIMB_T_SWAP(x, y) \
194 mp_limb_t __mp_limb_t_swap__tmp = (x); \
196 (y) = __mp_limb_t_swap__tmp; \
198 #define MP_SIZE_T_SWAP(x, y) \
200 mp_size_t __mp_size_t_swap__tmp = (x); \
202 (y) = __mp_size_t_swap__tmp; \
204 #define MP_BITCNT_T_SWAP(x,y) \
206 mp_bitcnt_t __mp_bitcnt_t_swap__tmp = (x); \
208 (y) = __mp_bitcnt_t_swap__tmp; \
210 #define MP_PTR_SWAP(x, y) \
212 mp_ptr __mp_ptr_swap__tmp = (x); \
214 (y) = __mp_ptr_swap__tmp; \
216 #define MP_SRCPTR_SWAP(x, y) \
218 mp_srcptr __mp_srcptr_swap__tmp = (x); \
220 (y) = __mp_srcptr_swap__tmp; \
223 #define MPN_PTR_SWAP(xp,xs, yp,ys) \
225 MP_PTR_SWAP (xp, yp); \
226 MP_SIZE_T_SWAP (xs, ys); \
228 #define MPN_SRCPTR_SWAP(xp,xs, yp,ys) \
230 MP_SRCPTR_SWAP (xp, yp); \
231 MP_SIZE_T_SWAP (xs, ys); \
234 #define MPZ_PTR_SWAP(x, y) \
236 mpz_ptr __mpz_ptr_swap__tmp = (x); \
238 (y) = __mpz_ptr_swap__tmp; \
240 #define MPZ_SRCPTR_SWAP(x, y) \
242 mpz_srcptr __mpz_srcptr_swap__tmp = (x); \
244 (y) = __mpz_srcptr_swap__tmp; \
247 const int mp_bits_per_limb = GMP_LIMB_BITS;
250 /* Memory allocation and other helper functions. */
252 gmp_die (const char *msg)
254 fprintf (stderr, "%s\n", msg);
259 gmp_default_alloc (size_t size)
267 gmp_die("gmp_default_alloc: Virtual memory exhausted.");
273 gmp_default_realloc (void *old, size_t old_size, size_t new_size)
277 p = realloc (old, new_size);
280 gmp_die("gmp_default_realloc: Virtual memory exhausted.");
286 gmp_default_free (void *p, size_t size)
291 static void * (*gmp_allocate_func) (size_t) = gmp_default_alloc;
292 static void * (*gmp_reallocate_func) (void *, size_t, size_t) = gmp_default_realloc;
293 static void (*gmp_free_func) (void *, size_t) = gmp_default_free;
296 mp_get_memory_functions (void *(**alloc_func) (size_t),
297 void *(**realloc_func) (void *, size_t, size_t),
298 void (**free_func) (void *, size_t))
301 *alloc_func = gmp_allocate_func;
304 *realloc_func = gmp_reallocate_func;
307 *free_func = gmp_free_func;
311 mp_set_memory_functions (void *(*alloc_func) (size_t),
312 void *(*realloc_func) (void *, size_t, size_t),
313 void (*free_func) (void *, size_t))
316 alloc_func = gmp_default_alloc;
318 realloc_func = gmp_default_realloc;
320 free_func = gmp_default_free;
322 gmp_allocate_func = alloc_func;
323 gmp_reallocate_func = realloc_func;
324 gmp_free_func = free_func;
327 #define gmp_xalloc(size) ((*gmp_allocate_func)((size)))
328 #define gmp_free(p) ((*gmp_free_func) ((p), 0))
331 gmp_xalloc_limbs (mp_size_t size)
333 return (mp_ptr) gmp_xalloc (size * sizeof (mp_limb_t));
337 gmp_xrealloc_limbs (mp_ptr old, mp_size_t size)
340 return (mp_ptr) (*gmp_reallocate_func) (old, 0, size * sizeof (mp_limb_t));
347 mpn_copyi (mp_ptr d, mp_srcptr s, mp_size_t n)
350 for (i = 0; i < n; i++)
355 mpn_copyd (mp_ptr d, mp_srcptr s, mp_size_t n)
362 mpn_cmp (mp_srcptr ap, mp_srcptr bp, mp_size_t n)
367 return ap[n] > bp[n] ? 1 : -1;
373 mpn_cmp4 (mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn)
376 return an < bn ? -1 : 1;
378 return mpn_cmp (ap, bp, an);
382 mpn_normalized_size (mp_srcptr xp, mp_size_t n)
384 while (n > 0 && xp[n-1] == 0)
390 mpn_zero_p(mp_srcptr rp, mp_size_t n)
392 return mpn_normalized_size (rp, n) == 0;
396 mpn_zero (mp_ptr rp, mp_size_t n)
403 mpn_add_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b)
411 mp_limb_t r = ap[i] + b;
422 mpn_add_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n)
427 for (i = 0, cy = 0; i < n; i++)
430 a = ap[i]; b = bp[i];
441 mpn_add (mp_ptr rp, mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn)
447 cy = mpn_add_n (rp, ap, bp, bn);
449 cy = mpn_add_1 (rp + bn, ap + bn, an - bn, cy);
454 mpn_sub_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b)
465 mp_limb_t cy = a < b;
475 mpn_sub_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n)
480 for (i = 0, cy = 0; i < n; i++)
483 a = ap[i]; b = bp[i];
493 mpn_sub (mp_ptr rp, mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn)
499 cy = mpn_sub_n (rp, ap, bp, bn);
501 cy = mpn_sub_1 (rp + bn, ap + bn, an - bn, cy);
506 mpn_mul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl)
508 mp_limb_t ul, cl, hpl, lpl;
516 gmp_umul_ppmm (hpl, lpl, ul, vl);
519 cl = (lpl < cl) + hpl;
529 mpn_addmul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl)
531 mp_limb_t ul, cl, hpl, lpl, rl;
539 gmp_umul_ppmm (hpl, lpl, ul, vl);
542 cl = (lpl < cl) + hpl;
555 mpn_submul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl)
557 mp_limb_t ul, cl, hpl, lpl, rl;
565 gmp_umul_ppmm (hpl, lpl, ul, vl);
568 cl = (lpl < cl) + hpl;
581 mpn_mul (mp_ptr rp, mp_srcptr up, mp_size_t un, mp_srcptr vp, mp_size_t vn)
585 assert (!GMP_MPN_OVERLAP_P(rp, un + vn, up, un));
586 assert (!GMP_MPN_OVERLAP_P(rp, un + vn, vp, vn));
588 /* We first multiply by the low order limb. This result can be
589 stored, not added, to rp. We also avoid a loop for zeroing this
592 rp[un] = mpn_mul_1 (rp, up, un, vp[0]);
594 /* Now accumulate the product of up[] and the next higher limb from
600 rp[un] = mpn_addmul_1 (rp, up, un, vp[0]);
606 mpn_mul_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n)
608 mpn_mul (rp, ap, n, bp, n);
612 mpn_sqr (mp_ptr rp, mp_srcptr ap, mp_size_t n)
614 mpn_mul (rp, ap, n, ap, n);
618 mpn_lshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt)
620 mp_limb_t high_limb, low_limb;
626 assert (cnt < GMP_LIMB_BITS);
631 tnc = GMP_LIMB_BITS - cnt;
633 retval = low_limb >> tnc;
634 high_limb = (low_limb << cnt);
639 *--rp = high_limb | (low_limb >> tnc);
640 high_limb = (low_limb << cnt);
648 mpn_rshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt)
650 mp_limb_t high_limb, low_limb;
656 assert (cnt < GMP_LIMB_BITS);
658 tnc = GMP_LIMB_BITS - cnt;
660 retval = (high_limb << tnc);
661 low_limb = high_limb >> cnt;
666 *rp++ = low_limb | (high_limb << tnc);
667 low_limb = high_limb >> cnt;
675 mpn_common_scan (mp_limb_t limb, mp_size_t i, mp_srcptr up, mp_size_t un,
680 assert (ux == 0 || ux == GMP_LIMB_MAX);
681 assert (0 <= i && i <= un );
687 return (ux == 0 ? ~(mp_bitcnt_t) 0 : un * GMP_LIMB_BITS);
691 return (mp_bitcnt_t) i * GMP_LIMB_BITS + cnt;
695 mpn_scan1 (mp_srcptr ptr, mp_bitcnt_t bit)
698 i = bit / GMP_LIMB_BITS;
700 return mpn_common_scan ( ptr[i] & (GMP_LIMB_MAX << (bit % GMP_LIMB_BITS)),
705 mpn_scan0 (mp_srcptr ptr, mp_bitcnt_t bit)
708 i = bit / GMP_LIMB_BITS;
710 return mpn_common_scan (~ptr[i] & (GMP_LIMB_MAX << (bit % GMP_LIMB_BITS)),
711 i, ptr, i, GMP_LIMB_MAX);
715 mpn_com (mp_ptr rp, mp_srcptr up, mp_size_t n)
722 mpn_neg (mp_ptr rp, mp_srcptr up, mp_size_t n)
732 mpn_com (++rp, ++up, --n);
737 /* MPN division interface. */
739 /* The 3/2 inverse is defined as
741 m = floor( (B^3-1) / (B u1 + u0)) - B
744 mpn_invert_3by2 (mp_limb_t u1, mp_limb_t u0)
746 mp_limb_t r, p, m, ql;
749 assert (u1 >= GMP_LIMB_HIGHBIT);
751 /* For notation, let b denote the half-limb base, so that B = b^2.
752 Split u1 = b uh + ul. */
753 ul = u1 & GMP_LLIMB_MASK;
754 uh = u1 >> (GMP_LIMB_BITS / 2);
756 /* Approximation of the high half of quotient. Differs from the 2/1
757 inverse of the half limb uh, since we have already subtracted
761 /* Adjust to get a half-limb 3/2 inverse, i.e., we want
763 qh' = floor( (b^3 - 1) / u) - b = floor ((b^3 - b u - 1) / u
764 = floor( (b (~u) + b-1) / u),
768 r = b (~u) + b-1 - qh (b uh + ul)
769 = b (~u - qh uh) + b-1 - qh ul
771 Subtraction of qh ul may underflow, which implies adjustments.
772 But by normalization, 2 u >= B > qh ul, so we need to adjust by
776 r = ((~u1 - (mp_limb_t) qh * uh) << (GMP_LIMB_BITS / 2)) | GMP_LLIMB_MASK;
778 p = (mp_limb_t) qh * ul;
779 /* Adjustment steps taken from udiv_qrnnd_c */
784 if (r >= u1) /* i.e. we didn't get carry when adding to r */
793 /* Low half of the quotient is
795 ql = floor ( (b r + b-1) / u1).
797 This is a 3/2 division (on half-limbs), for which qh is a
800 p = (r >> (GMP_LIMB_BITS / 2)) * qh + r;
801 /* Unlike full-limb 3/2, we can add 1 without overflow. For this to
802 work, it is essential that ql is a full mp_limb_t. */
803 ql = (p >> (GMP_LIMB_BITS / 2)) + 1;
805 /* By the 3/2 trick, we don't need the high half limb. */
806 r = (r << (GMP_LIMB_BITS / 2)) + GMP_LLIMB_MASK - ql * u1;
808 if (r >= (p << (GMP_LIMB_BITS / 2)))
813 m = ((mp_limb_t) qh << (GMP_LIMB_BITS / 2)) + ql;
820 /* Now m is the 2/1 invers of u1. If u0 > 0, adjust it to become a
837 gmp_umul_ppmm (th, tl, u0, m);
842 m -= ((r > u1) | ((r == u1) & (tl > u0)));
849 struct gmp_div_inverse
851 /* Normalization shift count. */
853 /* Normalized divisor (d0 unused for mpn_div_qr_1) */
855 /* Inverse, for 2/1 or 3/2. */
860 mpn_div_qr_1_invert (struct gmp_div_inverse *inv, mp_limb_t d)
867 inv->d1 = d << shift;
868 inv->di = mpn_invert_limb (inv->d1);
872 mpn_div_qr_2_invert (struct gmp_div_inverse *inv,
873 mp_limb_t d1, mp_limb_t d0)
882 d1 = (d1 << shift) | (d0 >> (GMP_LIMB_BITS - shift));
887 inv->di = mpn_invert_3by2 (d1, d0);
891 mpn_div_qr_invert (struct gmp_div_inverse *inv,
892 mp_srcptr dp, mp_size_t dn)
897 mpn_div_qr_1_invert (inv, dp[0]);
899 mpn_div_qr_2_invert (inv, dp[1], dp[0]);
912 d1 = (d1 << shift) | (d0 >> (GMP_LIMB_BITS - shift));
913 d0 = (d0 << shift) | (dp[dn-3] >> (GMP_LIMB_BITS - shift));
917 inv->di = mpn_invert_3by2 (d1, d0);
921 /* Not matching current public gmp interface, rather corresponding to
922 the sbpi1_div_* functions. */
924 mpn_div_qr_1_preinv (mp_ptr qp, mp_srcptr np, mp_size_t nn,
925 const struct gmp_div_inverse *inv)
933 tp = gmp_xalloc_limbs (nn);
934 r = mpn_lshift (tp, np, nn, inv->shift);
946 gmp_udiv_qrnnd_preinv (q, r, r, np[nn], d, di);
953 return r >> inv->shift;
957 mpn_div_qr_1 (mp_ptr qp, mp_srcptr np, mp_size_t nn, mp_limb_t d)
961 /* Special case for powers of two. */
962 if ((d & (d-1)) == 0)
964 mp_limb_t r = np[0] & (d-1);
968 mpn_copyi (qp, np, nn);
973 mpn_rshift (qp, np, nn, shift);
980 struct gmp_div_inverse inv;
981 mpn_div_qr_1_invert (&inv, d);
982 return mpn_div_qr_1_preinv (qp, np, nn, &inv);
987 mpn_div_qr_2_preinv (mp_ptr qp, mp_ptr rp, mp_srcptr np, mp_size_t nn,
988 const struct gmp_div_inverse *inv)
992 mp_limb_t d1, d0, di, r1, r0;
1003 tp = gmp_xalloc_limbs (nn);
1004 r1 = mpn_lshift (tp, np, nn, shift);
1017 gmp_udiv_qr_3by2 (q, r1, r0, r1, r0, n0, d1, d0, di);
1026 assert ((r0 << (GMP_LIMB_BITS - shift)) == 0);
1027 r0 = (r0 >> shift) | (r1 << (GMP_LIMB_BITS - shift));
1039 mpn_div_qr_2 (mp_ptr qp, mp_ptr rp, mp_srcptr np, mp_size_t nn,
1040 mp_limb_t d1, mp_limb_t d0)
1042 struct gmp_div_inverse inv;
1045 mpn_div_qr_2_invert (&inv, d1, d0);
1046 mpn_div_qr_2_preinv (qp, rp, np, nn, &inv);
1051 mpn_div_qr_pi1 (mp_ptr qp,
1052 mp_ptr np, mp_size_t nn, mp_limb_t n1,
1053 mp_srcptr dp, mp_size_t dn,
1068 assert ((d1 & GMP_LIMB_HIGHBIT) != 0);
1069 /* Iteration variable is the index of the q limb.
1071 * We divide <n1, np[dn-1+i], np[dn-2+i], np[dn-3+i],..., np[i]>
1072 * by <d1, d0, dp[dn-3], ..., dp[0] >
1078 mp_limb_t n0 = np[dn-1+i];
1080 if (n1 == d1 && n0 == d0)
1083 mpn_submul_1 (np+i, dp, dn, q);
1084 n1 = np[dn-1+i]; /* update n1, last loop's value will now be invalid */
1088 gmp_udiv_qr_3by2 (q, n1, n0, n1, n0, np[dn-2+i], d1, d0, dinv);
1090 cy = mpn_submul_1 (np + i, dp, dn-2, q);
1100 n1 += d1 + mpn_add_n (np + i, np + i, dp, dn - 1);
1114 mpn_div_qr_preinv (mp_ptr qp, mp_ptr np, mp_size_t nn,
1115 mp_srcptr dp, mp_size_t dn,
1116 const struct gmp_div_inverse *inv)
1122 np[0] = mpn_div_qr_1_preinv (qp, np, nn, inv);
1124 mpn_div_qr_2_preinv (qp, np, np, nn, inv);
1130 assert (inv->d1 == dp[dn-1]);
1131 assert (inv->d0 == dp[dn-2]);
1132 assert ((inv->d1 & GMP_LIMB_HIGHBIT) != 0);
1136 nh = mpn_lshift (np, np, nn, shift);
1140 mpn_div_qr_pi1 (qp, np, nn, nh, dp, dn, inv->di);
1143 gmp_assert_nocarry (mpn_rshift (np, np, dn, shift));
1148 mpn_div_qr (mp_ptr qp, mp_ptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn)
1150 struct gmp_div_inverse inv;
1156 mpn_div_qr_invert (&inv, dp, dn);
1157 if (dn > 2 && inv.shift > 0)
1159 tp = gmp_xalloc_limbs (dn);
1160 gmp_assert_nocarry (mpn_lshift (tp, dp, dn, inv.shift));
1163 mpn_div_qr_preinv (qp, np, nn, dp, dn, &inv);
1169 /* MPN base conversion. */
1171 mpn_base_power_of_two_p (unsigned b)
1187 struct mpn_base_info
1189 /* bb is the largest power of the base which fits in one limb, and
1190 exp is the corresponding exponent. */
1196 mpn_get_base_info (struct mpn_base_info *info, mp_limb_t b)
1202 m = GMP_LIMB_MAX / b;
1203 for (exp = 1, p = b; p <= m; exp++)
1211 mpn_limb_size_in_base_2 (mp_limb_t u)
1217 return GMP_LIMB_BITS - shift;
1221 mpn_get_str_bits (unsigned char *sp, unsigned bits, mp_srcptr up, mp_size_t un)
1228 sn = ((un - 1) * GMP_LIMB_BITS + mpn_limb_size_in_base_2 (up[un-1])
1231 mask = (1U << bits) - 1;
1233 for (i = 0, j = sn, shift = 0; j-- > 0;)
1235 unsigned char digit = up[i] >> shift;
1239 if (shift >= GMP_LIMB_BITS && ++i < un)
1241 shift -= GMP_LIMB_BITS;
1242 digit |= up[i] << (bits - shift);
1244 sp[j] = digit & mask;
1249 /* We generate digits from the least significant end, and reverse at
1252 mpn_limb_get_str (unsigned char *sp, mp_limb_t w,
1253 const struct gmp_div_inverse *binv)
1256 for (i = 0; w > 0; i++)
1260 h = w >> (GMP_LIMB_BITS - binv->shift);
1261 l = w << binv->shift;
1263 gmp_udiv_qrnnd_preinv (w, r, h, l, binv->d1, binv->di);
1264 assert ( (r << (GMP_LIMB_BITS - binv->shift)) == 0);
1273 mpn_get_str_other (unsigned char *sp,
1274 int base, const struct mpn_base_info *info,
1275 mp_ptr up, mp_size_t un)
1277 struct gmp_div_inverse binv;
1281 mpn_div_qr_1_invert (&binv, base);
1287 struct gmp_div_inverse bbinv;
1288 mpn_div_qr_1_invert (&bbinv, info->bb);
1294 w = mpn_div_qr_1_preinv (up, up, un, &bbinv);
1295 un -= (up[un-1] == 0);
1296 done = mpn_limb_get_str (sp + sn, w, &binv);
1298 for (sn += done; done < info->exp; done++)
1303 sn += mpn_limb_get_str (sp + sn, up[0], &binv);
1306 for (i = 0; 2*i + 1 < sn; i++)
1308 unsigned char t = sp[i];
1309 sp[i] = sp[sn - i - 1];
1317 mpn_get_str (unsigned char *sp, int base, mp_ptr up, mp_size_t un)
1322 assert (up[un-1] > 0);
1324 bits = mpn_base_power_of_two_p (base);
1326 return mpn_get_str_bits (sp, bits, up, un);
1329 struct mpn_base_info info;
1331 mpn_get_base_info (&info, base);
1332 return mpn_get_str_other (sp, base, &info, up, un);
1337 mpn_set_str_bits (mp_ptr rp, const unsigned char *sp, size_t sn,
1344 for (j = sn, rn = 0, shift = 0; j-- > 0; )
1353 rp[rn-1] |= (mp_limb_t) sp[j] << shift;
1355 if (shift >= GMP_LIMB_BITS)
1357 shift -= GMP_LIMB_BITS;
1359 rp[rn++] = (mp_limb_t) sp[j] >> (bits - shift);
1363 rn = mpn_normalized_size (rp, rn);
1367 /* Result is usually normalized, except for all-zero input, in which
1368 case a single zero limb is written at *RP, and 1 is returned. */
1370 mpn_set_str_other (mp_ptr rp, const unsigned char *sp, size_t sn,
1371 mp_limb_t b, const struct mpn_base_info *info)
1380 k = 1 + (sn - 1) % info->exp;
1385 w = w * b + sp[j++];
1389 for (rn = 1; j < sn;)
1394 for (k = 1; k < info->exp; k++)
1395 w = w * b + sp[j++];
1397 cy = mpn_mul_1 (rp, rp, rn, info->bb);
1398 cy += mpn_add_1 (rp, rp, rn, w);
1408 mpn_set_str (mp_ptr rp, const unsigned char *sp, size_t sn, int base)
1415 bits = mpn_base_power_of_two_p (base);
1417 return mpn_set_str_bits (rp, sp, sn, bits);
1420 struct mpn_base_info info;
1422 mpn_get_base_info (&info, base);
1423 return mpn_set_str_other (rp, sp, sn, base, &info);
1432 static const mp_limb_t dummy_limb = 0xc1a0;
1436 r->_mp_d = (mp_ptr) &dummy_limb;
1439 /* The utility of this function is a bit limited, since many functions
1440 assigns the result variable using mpz_swap. */
1442 mpz_init2 (mpz_t r, mp_bitcnt_t bits)
1446 bits -= (bits != 0); /* Round down, except if 0 */
1447 rn = 1 + bits / GMP_LIMB_BITS;
1451 r->_mp_d = gmp_xalloc_limbs (rn);
1458 gmp_free (r->_mp_d);
1462 mpz_realloc (mpz_t r, mp_size_t size)
1464 size = GMP_MAX (size, 1);
1467 r->_mp_d = gmp_xrealloc_limbs (r->_mp_d, size);
1469 r->_mp_d = gmp_xalloc_limbs (size);
1470 r->_mp_alloc = size;
1472 if (GMP_ABS (r->_mp_size) > size)
1478 /* Realloc for an mpz_t WHAT if it has less than NEEDED limbs. */
1479 #define MPZ_REALLOC(z,n) ((n) > (z)->_mp_alloc \
1480 ? mpz_realloc(z,n) \
1483 /* MPZ assignment and basic conversions. */
1485 mpz_set_si (mpz_t r, signed long int x)
1492 MPZ_REALLOC (r, 1)[0] = GMP_NEG_CAST (unsigned long int, x);
1497 mpz_set_ui (mpz_t r, unsigned long int x)
1502 MPZ_REALLOC (r, 1)[0] = x;
1509 mpz_set (mpz_t r, const mpz_t x)
1511 /* Allow the NOP r == x */
1517 n = GMP_ABS (x->_mp_size);
1518 rp = MPZ_REALLOC (r, n);
1520 mpn_copyi (rp, x->_mp_d, n);
1521 r->_mp_size = x->_mp_size;
1526 mpz_init_set_si (mpz_t r, signed long int x)
1533 mpz_init_set_ui (mpz_t r, unsigned long int x)
1540 mpz_init_set (mpz_t r, const mpz_t x)
1547 mpz_fits_slong_p (const mpz_t u)
1549 mp_size_t us = u->_mp_size;
1552 return u->_mp_d[0] < GMP_LIMB_HIGHBIT;
1554 return u->_mp_d[0] <= GMP_LIMB_HIGHBIT;
1560 mpz_fits_ulong_p (const mpz_t u)
1562 mp_size_t us = u->_mp_size;
1564 return (us == (us > 0));
1568 mpz_get_si (const mpz_t u)
1570 if (u->_mp_size < 0)
1571 /* This expression is necessary to properly handle 0x80000000 */
1572 return -1 - (long) ((u->_mp_d[0] - 1) & ~GMP_LIMB_HIGHBIT);
1574 return (long) (mpz_get_ui (u) & ~GMP_LIMB_HIGHBIT);
1578 mpz_get_ui (const mpz_t u)
1580 return u->_mp_size == 0 ? 0 : u->_mp_d[0];
1584 mpz_size (const mpz_t u)
1586 return GMP_ABS (u->_mp_size);
1590 mpz_getlimbn (const mpz_t u, mp_size_t n)
1592 if (n >= 0 && n < GMP_ABS (u->_mp_size))
1599 mpz_realloc2 (mpz_t x, mp_bitcnt_t n)
1601 mpz_realloc (x, 1 + (n - (n != 0)) / GMP_LIMB_BITS);
1605 mpz_limbs_read (mpz_srcptr x)
1611 mpz_limbs_modify (mpz_t x, mp_size_t n)
1614 return MPZ_REALLOC (x, n);
1618 mpz_limbs_write (mpz_t x, mp_size_t n)
1620 return mpz_limbs_modify (x, n);
1624 mpz_limbs_finish (mpz_t x, mp_size_t xs)
1627 xn = mpn_normalized_size (x->_mp_d, GMP_ABS (xs));
1628 x->_mp_size = xs < 0 ? -xn : xn;
1632 mpz_roinit_n (mpz_t x, mp_srcptr xp, mp_size_t xs)
1635 x->_mp_d = (mp_ptr) xp;
1636 mpz_limbs_finish (x, xs);
1641 /* Conversions and comparison to double. */
1643 mpz_set_d (mpz_t r, double x)
1652 /* x != x is true when x is a NaN, and x == x * 0.5 is true when x is
1653 zero or infinity. */
1654 if (x != x || x == x * 0.5)
1669 B = 2.0 * (double) GMP_LIMB_HIGHBIT;
1671 for (rn = 1; x >= B; rn++)
1674 rp = MPZ_REALLOC (r, rn);
1690 r->_mp_size = sign ? - rn : rn;
1694 mpz_init_set_d (mpz_t r, double x)
1701 mpz_get_d (const mpz_t u)
1705 double B = 2.0 * (double) GMP_LIMB_HIGHBIT;
1707 un = GMP_ABS (u->_mp_size);
1714 x = B*x + u->_mp_d[--un];
1716 if (u->_mp_size < 0)
1723 mpz_cmpabs_d (const mpz_t x, double d)
1736 B = 2.0 * (double) GMP_LIMB_HIGHBIT;
1739 /* Scale d so it can be compared with the top limb. */
1740 for (i = 1; i < xn; i++)
1746 /* Compare floor(d) to top limb, subtract and cancel when equal. */
1747 for (i = xn; i-- > 0;)
1764 mpz_cmp_d (const mpz_t x, double d)
1766 if (x->_mp_size < 0)
1771 return -mpz_cmpabs_d (x, d);
1778 return mpz_cmpabs_d (x, d);
1783 /* MPZ comparisons and the like. */
1785 mpz_sgn (const mpz_t u)
1787 return GMP_CMP (u->_mp_size, 0);
1791 mpz_cmp_si (const mpz_t u, long v)
1793 mp_size_t usize = u->_mp_size;
1798 return mpz_cmp_ui (u, v);
1799 else if (usize >= 0)
1801 else /* usize == -1 */
1802 return GMP_CMP (GMP_NEG_CAST (mp_limb_t, v), u->_mp_d[0]);
1806 mpz_cmp_ui (const mpz_t u, unsigned long v)
1808 mp_size_t usize = u->_mp_size;
1815 return GMP_CMP (mpz_get_ui (u), v);
1819 mpz_cmp (const mpz_t a, const mpz_t b)
1821 mp_size_t asize = a->_mp_size;
1822 mp_size_t bsize = b->_mp_size;
1825 return (asize < bsize) ? -1 : 1;
1826 else if (asize >= 0)
1827 return mpn_cmp (a->_mp_d, b->_mp_d, asize);
1829 return mpn_cmp (b->_mp_d, a->_mp_d, -asize);
1833 mpz_cmpabs_ui (const mpz_t u, unsigned long v)
1835 if (GMP_ABS (u->_mp_size) > 1)
1838 return GMP_CMP (mpz_get_ui (u), v);
1842 mpz_cmpabs (const mpz_t u, const mpz_t v)
1844 return mpn_cmp4 (u->_mp_d, GMP_ABS (u->_mp_size),
1845 v->_mp_d, GMP_ABS (v->_mp_size));
1849 mpz_abs (mpz_t r, const mpz_t u)
1852 r->_mp_size = GMP_ABS (r->_mp_size);
1856 mpz_neg (mpz_t r, const mpz_t u)
1859 r->_mp_size = -r->_mp_size;
1863 mpz_swap (mpz_t u, mpz_t v)
1865 MP_SIZE_T_SWAP (u->_mp_size, v->_mp_size);
1866 MP_SIZE_T_SWAP (u->_mp_alloc, v->_mp_alloc);
1867 MP_PTR_SWAP (u->_mp_d, v->_mp_d);
1871 /* MPZ addition and subtraction */
1873 /* Adds to the absolute value. Returns new size, but doesn't store it. */
1875 mpz_abs_add_ui (mpz_t r, const mpz_t a, unsigned long b)
1881 an = GMP_ABS (a->_mp_size);
1884 MPZ_REALLOC (r, 1)[0] = b;
1888 rp = MPZ_REALLOC (r, an + 1);
1890 cy = mpn_add_1 (rp, a->_mp_d, an, b);
1897 /* Subtract from the absolute value. Returns new size, (or -1 on underflow),
1898 but doesn't store it. */
1900 mpz_abs_sub_ui (mpz_t r, const mpz_t a, unsigned long b)
1902 mp_size_t an = GMP_ABS (a->_mp_size);
1907 MPZ_REALLOC (r, 1)[0] = b;
1910 rp = MPZ_REALLOC (r, an);
1911 if (an == 1 && a->_mp_d[0] < b)
1913 rp[0] = b - a->_mp_d[0];
1918 gmp_assert_nocarry (mpn_sub_1 (rp, a->_mp_d, an, b));
1919 return mpn_normalized_size (rp, an);
1924 mpz_add_ui (mpz_t r, const mpz_t a, unsigned long b)
1926 if (a->_mp_size >= 0)
1927 r->_mp_size = mpz_abs_add_ui (r, a, b);
1929 r->_mp_size = -mpz_abs_sub_ui (r, a, b);
1933 mpz_sub_ui (mpz_t r, const mpz_t a, unsigned long b)
1935 if (a->_mp_size < 0)
1936 r->_mp_size = -mpz_abs_add_ui (r, a, b);
1938 r->_mp_size = mpz_abs_sub_ui (r, a, b);
1942 mpz_ui_sub (mpz_t r, unsigned long a, const mpz_t b)
1944 if (b->_mp_size < 0)
1945 r->_mp_size = mpz_abs_add_ui (r, b, a);
1947 r->_mp_size = -mpz_abs_sub_ui (r, b, a);
1951 mpz_abs_add (mpz_t r, const mpz_t a, const mpz_t b)
1953 mp_size_t an = GMP_ABS (a->_mp_size);
1954 mp_size_t bn = GMP_ABS (b->_mp_size);
1960 MPZ_SRCPTR_SWAP (a, b);
1961 MP_SIZE_T_SWAP (an, bn);
1964 rp = MPZ_REALLOC (r, an + 1);
1965 cy = mpn_add (rp, a->_mp_d, an, b->_mp_d, bn);
1973 mpz_abs_sub (mpz_t r, const mpz_t a, const mpz_t b)
1975 mp_size_t an = GMP_ABS (a->_mp_size);
1976 mp_size_t bn = GMP_ABS (b->_mp_size);
1980 cmp = mpn_cmp4 (a->_mp_d, an, b->_mp_d, bn);
1983 rp = MPZ_REALLOC (r, an);
1984 gmp_assert_nocarry (mpn_sub (rp, a->_mp_d, an, b->_mp_d, bn));
1985 return mpn_normalized_size (rp, an);
1989 rp = MPZ_REALLOC (r, bn);
1990 gmp_assert_nocarry (mpn_sub (rp, b->_mp_d, bn, a->_mp_d, an));
1991 return -mpn_normalized_size (rp, bn);
1998 mpz_add (mpz_t r, const mpz_t a, const mpz_t b)
2002 if ( (a->_mp_size ^ b->_mp_size) >= 0)
2003 rn = mpz_abs_add (r, a, b);
2005 rn = mpz_abs_sub (r, a, b);
2007 r->_mp_size = a->_mp_size >= 0 ? rn : - rn;
2011 mpz_sub (mpz_t r, const mpz_t a, const mpz_t b)
2015 if ( (a->_mp_size ^ b->_mp_size) >= 0)
2016 rn = mpz_abs_sub (r, a, b);
2018 rn = mpz_abs_add (r, a, b);
2020 r->_mp_size = a->_mp_size >= 0 ? rn : - rn;
2024 /* MPZ multiplication */
2026 mpz_mul_si (mpz_t r, const mpz_t u, long int v)
2030 mpz_mul_ui (r, u, GMP_NEG_CAST (unsigned long int, v));
2034 mpz_mul_ui (r, u, (unsigned long int) v);
2038 mpz_mul_ui (mpz_t r, const mpz_t u, unsigned long int v)
2046 if (us == 0 || v == 0)
2054 tp = MPZ_REALLOC (r, un + 1);
2055 cy = mpn_mul_1 (tp, u->_mp_d, un, v);
2059 r->_mp_size = (us < 0) ? - un : un;
2063 mpz_mul (mpz_t r, const mpz_t u, const mpz_t v)
2066 mp_size_t un, vn, rn;
2073 if (un == 0 || vn == 0)
2079 sign = (un ^ vn) < 0;
2084 mpz_init2 (t, (un + vn) * GMP_LIMB_BITS);
2088 mpn_mul (tp, u->_mp_d, un, v->_mp_d, vn);
2090 mpn_mul (tp, v->_mp_d, vn, u->_mp_d, un);
2093 rn -= tp[rn-1] == 0;
2095 t->_mp_size = sign ? - rn : rn;
2101 mpz_mul_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bits)
2108 un = GMP_ABS (u->_mp_size);
2115 limbs = bits / GMP_LIMB_BITS;
2116 shift = bits % GMP_LIMB_BITS;
2118 rn = un + limbs + (shift > 0);
2119 rp = MPZ_REALLOC (r, rn);
2122 mp_limb_t cy = mpn_lshift (rp + limbs, u->_mp_d, un, shift);
2127 mpn_copyd (rp + limbs, u->_mp_d, un);
2129 mpn_zero (rp, limbs);
2131 r->_mp_size = (u->_mp_size < 0) ? - rn : rn;
2135 mpz_addmul_ui (mpz_t r, const mpz_t u, unsigned long int v)
2139 mpz_mul_ui (t, u, v);
2145 mpz_submul_ui (mpz_t r, const mpz_t u, unsigned long int v)
2149 mpz_mul_ui (t, u, v);
2155 mpz_addmul (mpz_t r, const mpz_t u, const mpz_t v)
2165 mpz_submul (mpz_t r, const mpz_t u, const mpz_t v)
2176 enum mpz_div_round_mode { GMP_DIV_FLOOR, GMP_DIV_CEIL, GMP_DIV_TRUNC };
2178 /* Allows q or r to be zero. Returns 1 iff remainder is non-zero. */
2180 mpz_div_qr (mpz_t q, mpz_t r,
2181 const mpz_t n, const mpz_t d, enum mpz_div_round_mode mode)
2183 mp_size_t ns, ds, nn, dn, qs;
2188 gmp_die("mpz_div_qr: Divide by zero.");
2206 if (mode == GMP_DIV_CEIL && qs >= 0)
2208 /* q = 1, r = n - d */
2214 else if (mode == GMP_DIV_FLOOR && qs < 0)
2216 /* q = -1, r = n + d */
2238 mpz_init_set (tr, n);
2245 mpz_init2 (tq, qn * GMP_LIMB_BITS);
2251 mpn_div_qr (qp, np, nn, d->_mp_d, dn);
2255 qn -= (qp[qn-1] == 0);
2257 tq->_mp_size = qs < 0 ? -qn : qn;
2259 rn = mpn_normalized_size (np, dn);
2260 tr->_mp_size = ns < 0 ? - rn : rn;
2262 if (mode == GMP_DIV_FLOOR && qs < 0 && rn != 0)
2265 mpz_sub_ui (tq, tq, 1);
2267 mpz_add (tr, tr, d);
2269 else if (mode == GMP_DIV_CEIL && qs >= 0 && rn != 0)
2272 mpz_add_ui (tq, tq, 1);
2274 mpz_sub (tr, tr, d);
2292 mpz_cdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d)
2294 mpz_div_qr (q, r, n, d, GMP_DIV_CEIL);
2298 mpz_fdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d)
2300 mpz_div_qr (q, r, n, d, GMP_DIV_FLOOR);
2304 mpz_tdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d)
2306 mpz_div_qr (q, r, n, d, GMP_DIV_TRUNC);
2310 mpz_cdiv_q (mpz_t q, const mpz_t n, const mpz_t d)
2312 mpz_div_qr (q, NULL, n, d, GMP_DIV_CEIL);
2316 mpz_fdiv_q (mpz_t q, const mpz_t n, const mpz_t d)
2318 mpz_div_qr (q, NULL, n, d, GMP_DIV_FLOOR);
2322 mpz_tdiv_q (mpz_t q, const mpz_t n, const mpz_t d)
2324 mpz_div_qr (q, NULL, n, d, GMP_DIV_TRUNC);
2328 mpz_cdiv_r (mpz_t r, const mpz_t n, const mpz_t d)
2330 mpz_div_qr (NULL, r, n, d, GMP_DIV_CEIL);
2334 mpz_fdiv_r (mpz_t r, const mpz_t n, const mpz_t d)
2336 mpz_div_qr (NULL, r, n, d, GMP_DIV_FLOOR);
2340 mpz_tdiv_r (mpz_t r, const mpz_t n, const mpz_t d)
2342 mpz_div_qr (NULL, r, n, d, GMP_DIV_TRUNC);
2346 mpz_mod (mpz_t r, const mpz_t n, const mpz_t d)
2348 mpz_div_qr (NULL, r, n, d, d->_mp_size >= 0 ? GMP_DIV_FLOOR : GMP_DIV_CEIL);
2352 mpz_div_q_2exp (mpz_t q, const mpz_t u, mp_bitcnt_t bit_index,
2353 enum mpz_div_round_mode mode)
2366 limb_cnt = bit_index / GMP_LIMB_BITS;
2367 qn = GMP_ABS (un) - limb_cnt;
2368 bit_index %= GMP_LIMB_BITS;
2370 if (mode == ((un > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* un != 0 here. */
2371 /* Note: Below, the final indexing at limb_cnt is valid because at
2372 that point we have qn > 0. */
2374 || !mpn_zero_p (u->_mp_d, limb_cnt)
2375 || (u->_mp_d[limb_cnt]
2376 & (((mp_limb_t) 1 << bit_index) - 1)));
2384 qp = MPZ_REALLOC (q, qn);
2388 mpn_rshift (qp, u->_mp_d + limb_cnt, qn, bit_index);
2389 qn -= qp[qn - 1] == 0;
2393 mpn_copyi (qp, u->_mp_d + limb_cnt, qn);
2400 mpz_add_ui (q, q, 1);
2406 mpz_div_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bit_index,
2407 enum mpz_div_round_mode mode)
2409 mp_size_t us, un, rn;
2414 if (us == 0 || bit_index == 0)
2419 rn = (bit_index + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS;
2422 rp = MPZ_REALLOC (r, rn);
2425 mask = GMP_LIMB_MAX >> (rn * GMP_LIMB_BITS - bit_index);
2429 /* Quotient (with truncation) is zero, and remainder is
2431 if (mode == ((us > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* us != 0 here. */
2433 /* Have to negate and sign extend. */
2436 gmp_assert_nocarry (! mpn_neg (rp, u->_mp_d, un));
2437 for (i = un; i < rn - 1; i++)
2438 rp[i] = GMP_LIMB_MAX;
2447 mpn_copyi (rp, u->_mp_d, un);
2455 mpn_copyi (rp, u->_mp_d, rn - 1);
2457 rp[rn-1] = u->_mp_d[rn-1] & mask;
2459 if (mode == ((us > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* us != 0 here. */
2461 /* If r != 0, compute 2^{bit_count} - r. */
2462 mpn_neg (rp, rp, rn);
2466 /* us is not used for anything else, so we can modify it
2467 here to indicate flipped sign. */
2471 rn = mpn_normalized_size (rp, rn);
2472 r->_mp_size = us < 0 ? -rn : rn;
2476 mpz_cdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2478 mpz_div_q_2exp (r, u, cnt, GMP_DIV_CEIL);
2482 mpz_fdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2484 mpz_div_q_2exp (r, u, cnt, GMP_DIV_FLOOR);
2488 mpz_tdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2490 mpz_div_q_2exp (r, u, cnt, GMP_DIV_TRUNC);
2494 mpz_cdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2496 mpz_div_r_2exp (r, u, cnt, GMP_DIV_CEIL);
2500 mpz_fdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2502 mpz_div_r_2exp (r, u, cnt, GMP_DIV_FLOOR);
2506 mpz_tdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2508 mpz_div_r_2exp (r, u, cnt, GMP_DIV_TRUNC);
2512 mpz_divexact (mpz_t q, const mpz_t n, const mpz_t d)
2514 gmp_assert_nocarry (mpz_div_qr (q, NULL, n, d, GMP_DIV_TRUNC));
2518 mpz_divisible_p (const mpz_t n, const mpz_t d)
2520 return mpz_div_qr (NULL, NULL, n, d, GMP_DIV_TRUNC) == 0;
2524 mpz_congruent_p (const mpz_t a, const mpz_t b, const mpz_t m)
2529 /* a == b (mod 0) iff a == b */
2530 if (mpz_sgn (m) == 0)
2531 return (mpz_cmp (a, b) == 0);
2535 res = mpz_divisible_p (t, m);
2541 static unsigned long
2542 mpz_div_qr_ui (mpz_t q, mpz_t r,
2543 const mpz_t n, unsigned long d, enum mpz_div_round_mode mode)
2562 qp = MPZ_REALLOC (q, qn);
2566 rl = mpn_div_qr_1 (qp, n->_mp_d, qn, d);
2570 rs = (ns < 0) ? -rs : rs;
2572 if (rl > 0 && ( (mode == GMP_DIV_FLOOR && ns < 0)
2573 || (mode == GMP_DIV_CEIL && ns >= 0)))
2576 gmp_assert_nocarry (mpn_add_1 (qp, qp, qn, 1));
2583 MPZ_REALLOC (r, 1)[0] = rl;
2588 qn -= (qp[qn-1] == 0);
2589 assert (qn == 0 || qp[qn-1] > 0);
2591 q->_mp_size = (ns < 0) ? - qn : qn;
2598 mpz_cdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d)
2600 return mpz_div_qr_ui (q, r, n, d, GMP_DIV_CEIL);
2604 mpz_fdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d)
2606 return mpz_div_qr_ui (q, r, n, d, GMP_DIV_FLOOR);
2610 mpz_tdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d)
2612 return mpz_div_qr_ui (q, r, n, d, GMP_DIV_TRUNC);
2616 mpz_cdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d)
2618 return mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_CEIL);
2622 mpz_fdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d)
2624 return mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_FLOOR);
2628 mpz_tdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d)
2630 return mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_TRUNC);
2634 mpz_cdiv_r_ui (mpz_t r, const mpz_t n, unsigned long d)
2636 return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_CEIL);
2639 mpz_fdiv_r_ui (mpz_t r, const mpz_t n, unsigned long d)
2641 return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_FLOOR);
2644 mpz_tdiv_r_ui (mpz_t r, const mpz_t n, unsigned long d)
2646 return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_TRUNC);
2650 mpz_cdiv_ui (const mpz_t n, unsigned long d)
2652 return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_CEIL);
2656 mpz_fdiv_ui (const mpz_t n, unsigned long d)
2658 return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_FLOOR);
2662 mpz_tdiv_ui (const mpz_t n, unsigned long d)
2664 return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_TRUNC);
2668 mpz_mod_ui (mpz_t r, const mpz_t n, unsigned long d)
2670 return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_FLOOR);
2674 mpz_divexact_ui (mpz_t q, const mpz_t n, unsigned long d)
2676 gmp_assert_nocarry (mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_TRUNC));
2680 mpz_divisible_ui_p (const mpz_t n, unsigned long d)
2682 return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_TRUNC) == 0;
2688 mpn_gcd_11 (mp_limb_t u, mp_limb_t v)
2692 assert ( (u | v) > 0);
2699 gmp_ctz (shift, u | v);
2705 MP_LIMB_T_SWAP (u, v);
2707 while ( (v & 1) == 0)
2717 while ( (u & 1) == 0);
2724 while ( (v & 1) == 0);
2731 mpz_gcd_ui (mpz_t g, const mpz_t u, unsigned long v)
2742 un = GMP_ABS (u->_mp_size);
2744 v = mpn_gcd_11 (mpn_div_qr_1 (NULL, u->_mp_d, un, v), v);
2754 mpz_make_odd (mpz_t r)
2758 assert (r->_mp_size > 0);
2759 /* Count trailing zeros, equivalent to mpn_scan1, because we know that there is a 1 */
2760 shift = mpn_common_scan (r->_mp_d[0], 0, r->_mp_d, 0, 0);
2761 mpz_tdiv_q_2exp (r, r, shift);
2767 mpz_gcd (mpz_t g, const mpz_t u, const mpz_t v)
2770 mp_bitcnt_t uz, vz, gz;
2772 if (u->_mp_size == 0)
2777 if (v->_mp_size == 0)
2787 uz = mpz_make_odd (tu);
2789 vz = mpz_make_odd (tv);
2790 gz = GMP_MIN (uz, vz);
2792 if (tu->_mp_size < tv->_mp_size)
2795 mpz_tdiv_r (tu, tu, tv);
2796 if (tu->_mp_size == 0)
2806 c = mpz_cmp (tu, tv);
2815 if (tv->_mp_size == 1)
2817 mp_limb_t vl = tv->_mp_d[0];
2818 mp_limb_t ul = mpz_tdiv_ui (tu, vl);
2819 mpz_set_ui (g, mpn_gcd_11 (ul, vl));
2822 mpz_sub (tu, tu, tv);
2826 mpz_mul_2exp (g, g, gz);
2830 mpz_gcdext (mpz_t g, mpz_t s, mpz_t t, const mpz_t u, const mpz_t v)
2832 mpz_t tu, tv, s0, s1, t0, t1;
2833 mp_bitcnt_t uz, vz, gz;
2836 if (u->_mp_size == 0)
2838 /* g = 0 u + sgn(v) v */
2839 signed long sign = mpz_sgn (v);
2844 mpz_set_si (t, sign);
2848 if (v->_mp_size == 0)
2850 /* g = sgn(u) u + 0 v */
2851 signed long sign = mpz_sgn (u);
2854 mpz_set_si (s, sign);
2868 uz = mpz_make_odd (tu);
2870 vz = mpz_make_odd (tv);
2871 gz = GMP_MIN (uz, vz);
2876 /* Cofactors corresponding to odd gcd. gz handled later. */
2877 if (tu->_mp_size < tv->_mp_size)
2880 MPZ_SRCPTR_SWAP (u, v);
2881 MPZ_PTR_SWAP (s, t);
2882 MP_BITCNT_T_SWAP (uz, vz);
2890 * where u and v denote the inputs with common factors of two
2891 * eliminated, and det (s0, t0; s1, t1) = 2^p. Then
2893 * 2^p tu = s1 u - t1 v
2894 * 2^p tv = -s0 u + t0 v
2897 /* After initial division, tu = q tv + tu', we have
2899 * u = 2^uz (tu' + q tv)
2904 * t0 = 2^uz, t1 = 2^uz q
2908 mpz_setbit (t0, uz);
2909 mpz_tdiv_qr (t1, tu, tu, tv);
2910 mpz_mul_2exp (t1, t1, uz);
2912 mpz_setbit (s1, vz);
2915 if (tu->_mp_size > 0)
2918 shift = mpz_make_odd (tu);
2919 mpz_mul_2exp (t0, t0, shift);
2920 mpz_mul_2exp (s0, s0, shift);
2926 c = mpz_cmp (tu, tv);
2934 * u = t0 tu + t1 (tv' + tu) = (t0 + t1) tu + t1 tv'
2935 * v = s0 tu + s1 (tv' + tu) = (s0 + s1) tu + s1 tv' */
2937 mpz_sub (tv, tv, tu);
2938 mpz_add (t0, t0, t1);
2939 mpz_add (s0, s0, s1);
2941 shift = mpz_make_odd (tv);
2942 mpz_mul_2exp (t1, t1, shift);
2943 mpz_mul_2exp (s1, s1, shift);
2947 mpz_sub (tu, tu, tv);
2948 mpz_add (t1, t0, t1);
2949 mpz_add (s1, s0, s1);
2951 shift = mpz_make_odd (tu);
2952 mpz_mul_2exp (t0, t0, shift);
2953 mpz_mul_2exp (s0, s0, shift);
2959 /* Now tv = odd part of gcd, and -s0 and t0 are corresponding
2962 mpz_mul_2exp (tv, tv, gz);
2965 /* 2^p g = s0 u + t0 v. Eliminate one factor of two at a time. To
2966 adjust cofactors, we need u / g and v / g */
2968 mpz_divexact (s1, v, tv);
2970 mpz_divexact (t1, u, tv);
2975 /* s0 u + t0 v = (s0 - v/g) u - (t0 + u/g) v */
2976 if (mpz_odd_p (s0) || mpz_odd_p (t0))
2978 mpz_sub (s0, s0, s1);
2979 mpz_add (t0, t0, t1);
2981 mpz_divexact_ui (s0, s0, 2);
2982 mpz_divexact_ui (t0, t0, 2);
2985 /* Arrange so that |s| < |u| / 2g */
2986 mpz_add (s1, s0, s1);
2987 if (mpz_cmpabs (s0, s1) > 0)
2990 mpz_sub (t0, t0, t1);
2992 if (u->_mp_size < 0)
2994 if (v->_mp_size < 0)
3012 mpz_lcm (mpz_t r, const mpz_t u, const mpz_t v)
3016 if (u->_mp_size == 0 || v->_mp_size == 0)
3025 mpz_divexact (g, u, g);
3033 mpz_lcm_ui (mpz_t r, const mpz_t u, unsigned long v)
3035 if (v == 0 || u->_mp_size == 0)
3041 v /= mpz_gcd_ui (NULL, u, v);
3042 mpz_mul_ui (r, u, v);
3048 mpz_invert (mpz_t r, const mpz_t u, const mpz_t m)
3053 if (u->_mp_size == 0 || mpz_cmpabs_ui (m, 1) <= 0)
3059 mpz_gcdext (g, tr, NULL, u, m);
3060 invertible = (mpz_cmp_ui (g, 1) == 0);
3064 if (tr->_mp_size < 0)
3066 if (m->_mp_size >= 0)
3067 mpz_add (tr, tr, m);
3069 mpz_sub (tr, tr, m);
3080 /* Higher level operations (sqrt, pow and root) */
3083 mpz_pow_ui (mpz_t r, const mpz_t b, unsigned long e)
3087 mpz_init_set_ui (tr, 1);
3089 bit = GMP_ULONG_HIGHBIT;
3092 mpz_mul (tr, tr, tr);
3094 mpz_mul (tr, tr, b);
3104 mpz_ui_pow_ui (mpz_t r, unsigned long blimb, unsigned long e)
3107 mpz_pow_ui (r, mpz_roinit_n (b, &blimb, 1), e);
3111 mpz_powm (mpz_t r, const mpz_t b, const mpz_t e, const mpz_t m)
3117 struct gmp_div_inverse minv;
3121 en = GMP_ABS (e->_mp_size);
3122 mn = GMP_ABS (m->_mp_size);
3124 gmp_die ("mpz_powm: Zero modulo.");
3133 mpn_div_qr_invert (&minv, mp, mn);
3138 /* To avoid shifts, we do all our reductions, except the final
3139 one, using a *normalized* m. */
3142 tp = gmp_xalloc_limbs (mn);
3143 gmp_assert_nocarry (mpn_lshift (tp, mp, mn, shift));
3149 if (e->_mp_size < 0)
3151 if (!mpz_invert (base, b, m))
3152 gmp_die ("mpz_powm: Negative exponent and non-invertible base.");
3159 bn = base->_mp_size;
3162 mpn_div_qr_preinv (NULL, base->_mp_d, base->_mp_size, mp, mn, &minv);
3166 /* We have reduced the absolute value. Now take care of the
3167 sign. Note that we get zero represented non-canonically as
3169 if (b->_mp_size < 0)
3171 mp_ptr bp = MPZ_REALLOC (base, mn);
3172 gmp_assert_nocarry (mpn_sub (bp, mp, mn, bp, bn));
3175 base->_mp_size = mpn_normalized_size (base->_mp_d, bn);
3177 mpz_init_set_ui (tr, 1);
3181 mp_limb_t w = e->_mp_d[en];
3184 bit = GMP_LIMB_HIGHBIT;
3187 mpz_mul (tr, tr, tr);
3189 mpz_mul (tr, tr, base);
3190 if (tr->_mp_size > mn)
3192 mpn_div_qr_preinv (NULL, tr->_mp_d, tr->_mp_size, mp, mn, &minv);
3193 tr->_mp_size = mpn_normalized_size (tr->_mp_d, mn);
3200 /* Final reduction */
3201 if (tr->_mp_size >= mn)
3204 mpn_div_qr_preinv (NULL, tr->_mp_d, tr->_mp_size, mp, mn, &minv);
3205 tr->_mp_size = mpn_normalized_size (tr->_mp_d, mn);
3216 mpz_powm_ui (mpz_t r, const mpz_t b, unsigned long elimb, const mpz_t m)
3219 mpz_powm (r, b, mpz_roinit_n (e, &elimb, 1), m);
3222 /* x=trunc(y^(1/z)), r=y-x^z */
3224 mpz_rootrem (mpz_t x, mpz_t r, const mpz_t y, unsigned long z)
3229 sgn = y->_mp_size < 0;
3230 if ((~z & sgn) != 0)
3231 gmp_die ("mpz_rootrem: Negative argument, with even root.");
3233 gmp_die ("mpz_rootrem: Zeroth root.");
3235 if (mpz_cmpabs_ui (y, 1) <= 0) {
3245 mpz_setbit (t, mpz_sizeinbase (y, 2) / z + 1);
3247 if (z == 2) /* simplify sqrt loop: z-1 == 1 */
3249 mpz_swap (u, t); /* u = x */
3250 mpz_tdiv_q (t, y, u); /* t = y/x */
3251 mpz_add (t, t, u); /* t = y/x + x */
3252 mpz_tdiv_q_2exp (t, t, 1); /* x'= (y/x + x)/2 */
3253 } while (mpz_cmpabs (t, u) < 0); /* |x'| < |x| */
3262 mpz_swap (u, t); /* u = x */
3263 mpz_pow_ui (t, u, z - 1); /* t = x^(z-1) */
3264 mpz_tdiv_q (t, y, t); /* t = y/x^(z-1) */
3265 mpz_mul_ui (v, u, z - 1); /* v = x*(z-1) */
3266 mpz_add (t, t, v); /* t = y/x^(z-1) + x*(z-1) */
3267 mpz_tdiv_q_ui (t, t, z); /* x'=(y/x^(z-1) + x*(z-1))/z */
3268 } while (mpz_cmpabs (t, u) < 0); /* |x'| < |x| */
3274 mpz_pow_ui (t, u, z);
3284 mpz_root (mpz_t x, const mpz_t y, unsigned long z)
3290 mpz_rootrem (x, r, y, z);
3291 res = r->_mp_size == 0;
3297 /* Compute s = floor(sqrt(u)) and r = u - s^2. Allows r == NULL */
3299 mpz_sqrtrem (mpz_t s, mpz_t r, const mpz_t u)
3301 mpz_rootrem (s, r, u, 2);
3305 mpz_sqrt (mpz_t s, const mpz_t u)
3307 mpz_rootrem (s, NULL, u, 2);
3311 mpz_perfect_square_p (const mpz_t u)
3313 if (u->_mp_size <= 0)
3314 return (u->_mp_size == 0);
3316 return mpz_root (NULL, u, 2);
3320 mpn_perfect_square_p (mp_srcptr p, mp_size_t n)
3325 assert (p [n-1] != 0);
3326 return mpz_root (NULL, mpz_roinit_n (t, p, n), 2);
3330 mpn_sqrtrem (mp_ptr sp, mp_ptr rp, mp_srcptr p, mp_size_t n)
3336 assert (p [n-1] != 0);
3340 mpz_rootrem (s, r, mpz_roinit_n (u, p, n), 2);
3342 assert (s->_mp_size == (n+1)/2);
3343 mpn_copyd (sp, s->_mp_d, s->_mp_size);
3347 mpn_copyd (rp, r->_mp_d, res);
3355 mpz_fac_ui (mpz_t x, unsigned long n)
3357 mpz_set_ui (x, n + (n == 0));
3359 mpz_mul_ui (x, x, --n);
3363 mpz_bin_uiui (mpz_t r, unsigned long n, unsigned long k)
3367 mpz_set_ui (r, k <= n);
3370 k = (k <= n) ? n - k : 0;
3376 mpz_mul_ui (r, r, n--);
3378 mpz_divexact (r, r, t);
3383 /* Primality testing */
3385 gmp_millerrabin (const mpz_t n, const mpz_t nm1, mpz_t y,
3386 const mpz_t q, mp_bitcnt_t k)
3390 /* Caller must initialize y to the base. */
3391 mpz_powm (y, y, q, n);
3393 if (mpz_cmp_ui (y, 1) == 0 || mpz_cmp (y, nm1) == 0)
3398 mpz_powm_ui (y, y, 2, n);
3399 if (mpz_cmp (y, nm1) == 0)
3401 /* y == 1 means that the previous y was a non-trivial square root
3402 of 1 (mod n). y == 0 means that n is a power of the base.
3403 In either case, n is not prime. */
3404 if (mpz_cmp_ui (y, 1) <= 0)
3410 /* This product is 0xc0cfd797, and fits in 32 bits. */
3411 #define GMP_PRIME_PRODUCT \
3412 (3UL*5UL*7UL*11UL*13UL*17UL*19UL*23UL*29UL)
3414 /* Bit (p+1)/2 is set, for each odd prime <= 61 */
3415 #define GMP_PRIME_MASK 0xc96996dcUL
3418 mpz_probab_prime_p (const mpz_t n, int reps)
3427 /* Note that we use the absolute value of n only, for compatibility
3428 with the real GMP. */
3430 return (mpz_cmpabs_ui (n, 2) == 0) ? 2 : 0;
3432 /* Above test excludes n == 0 */
3433 assert (n->_mp_size != 0);
3435 if (mpz_cmpabs_ui (n, 64) < 0)
3436 return (GMP_PRIME_MASK >> (n->_mp_d[0] >> 1)) & 2;
3438 if (mpz_gcd_ui (NULL, n, GMP_PRIME_PRODUCT) != 1)
3441 /* All prime factors are >= 31. */
3442 if (mpz_cmpabs_ui (n, 31*31) < 0)
3445 /* Use Miller-Rabin, with a deterministic sequence of bases, a[j] =
3446 j^2 + j + 41 using Euler's polynomial. We potentially stop early,
3447 if a[j] >= n - 1. Since n >= 31*31, this can happen only if reps >
3448 30 (a[30] == 971 > 31*31 == 961). */
3454 /* Find q and k, where q is odd and n = 1 + 2**k * q. */
3455 nm1->_mp_size = mpz_abs_sub_ui (nm1, n, 1);
3456 k = mpz_scan1 (nm1, 0);
3457 mpz_tdiv_q_2exp (q, nm1, k);
3459 for (j = 0, is_prime = 1; is_prime & (j < reps); j++)
3461 mpz_set_ui (y, (unsigned long) j*j+j+41);
3462 if (mpz_cmp (y, nm1) >= 0)
3464 /* Don't try any further bases. This "early" break does not affect
3465 the result for any reasonable reps value (<=5000 was tested) */
3469 is_prime = gmp_millerrabin (n, nm1, y, q, k);
3479 /* Logical operations and bit manipulation. */
3481 /* Numbers are treated as if represented in two's complement (and
3482 infinitely sign extended). For a negative values we get the two's
3483 complement from -x = ~x + 1, where ~ is bitwise complement.
3492 where yyyy is the bitwise complement of xxxx. So least significant
3493 bits, up to and including the first one bit, are unchanged, and
3494 the more significant bits are all complemented.
3496 To change a bit from zero to one in a negative number, subtract the
3497 corresponding power of two from the absolute value. This can never
3498 underflow. To change a bit from one to zero, add the corresponding
3499 power of two, and this might overflow. E.g., if x = -001111, the
3500 two's complement is 110001. Clearing the least significant bit, we
3501 get two's complement 110000, and -010000. */
3504 mpz_tstbit (const mpz_t d, mp_bitcnt_t bit_index)
3506 mp_size_t limb_index;
3515 limb_index = bit_index / GMP_LIMB_BITS;
3516 if (limb_index >= dn)
3519 shift = bit_index % GMP_LIMB_BITS;
3520 w = d->_mp_d[limb_index];
3521 bit = (w >> shift) & 1;
3525 /* d < 0. Check if any of the bits below is set: If so, our bit
3526 must be complemented. */
3527 if (shift > 0 && (w << (GMP_LIMB_BITS - shift)) > 0)
3529 while (--limb_index >= 0)
3530 if (d->_mp_d[limb_index] > 0)
3537 mpz_abs_add_bit (mpz_t d, mp_bitcnt_t bit_index)
3539 mp_size_t dn, limb_index;
3543 dn = GMP_ABS (d->_mp_size);
3545 limb_index = bit_index / GMP_LIMB_BITS;
3546 bit = (mp_limb_t) 1 << (bit_index % GMP_LIMB_BITS);
3548 if (limb_index >= dn)
3551 /* The bit should be set outside of the end of the number.
3552 We have to increase the size of the number. */
3553 dp = MPZ_REALLOC (d, limb_index + 1);
3555 dp[limb_index] = bit;
3556 for (i = dn; i < limb_index; i++)
3558 dn = limb_index + 1;
3566 cy = mpn_add_1 (dp + limb_index, dp + limb_index, dn - limb_index, bit);
3569 dp = MPZ_REALLOC (d, dn + 1);
3574 d->_mp_size = (d->_mp_size < 0) ? - dn : dn;
3578 mpz_abs_sub_bit (mpz_t d, mp_bitcnt_t bit_index)
3580 mp_size_t dn, limb_index;
3584 dn = GMP_ABS (d->_mp_size);
3587 limb_index = bit_index / GMP_LIMB_BITS;
3588 bit = (mp_limb_t) 1 << (bit_index % GMP_LIMB_BITS);
3590 assert (limb_index < dn);
3592 gmp_assert_nocarry (mpn_sub_1 (dp + limb_index, dp + limb_index,
3593 dn - limb_index, bit));
3594 dn = mpn_normalized_size (dp, dn);
3595 d->_mp_size = (d->_mp_size < 0) ? - dn : dn;
3599 mpz_setbit (mpz_t d, mp_bitcnt_t bit_index)
3601 if (!mpz_tstbit (d, bit_index))
3603 if (d->_mp_size >= 0)
3604 mpz_abs_add_bit (d, bit_index);
3606 mpz_abs_sub_bit (d, bit_index);
3611 mpz_clrbit (mpz_t d, mp_bitcnt_t bit_index)
3613 if (mpz_tstbit (d, bit_index))
3615 if (d->_mp_size >= 0)
3616 mpz_abs_sub_bit (d, bit_index);
3618 mpz_abs_add_bit (d, bit_index);
3623 mpz_combit (mpz_t d, mp_bitcnt_t bit_index)
3625 if (mpz_tstbit (d, bit_index) ^ (d->_mp_size < 0))
3626 mpz_abs_sub_bit (d, bit_index);
3628 mpz_abs_add_bit (d, bit_index);
3632 mpz_com (mpz_t r, const mpz_t u)
3635 mpz_sub_ui (r, r, 1);
3639 mpz_and (mpz_t r, const mpz_t u, const mpz_t v)
3641 mp_size_t un, vn, rn, i;
3644 mp_limb_t ux, vx, rx;
3645 mp_limb_t uc, vc, rc;
3646 mp_limb_t ul, vl, rl;
3648 un = GMP_ABS (u->_mp_size);
3649 vn = GMP_ABS (v->_mp_size);
3652 MPZ_SRCPTR_SWAP (u, v);
3653 MP_SIZE_T_SWAP (un, vn);
3661 uc = u->_mp_size < 0;
3662 vc = v->_mp_size < 0;
3669 /* If the smaller input is positive, higher limbs don't matter. */
3672 rp = MPZ_REALLOC (r, rn + (mp_size_t) rc);
3680 ul = (up[i] ^ ux) + uc;
3683 vl = (vp[i] ^ vx) + vc;
3686 rl = ( (ul & vl) ^ rx) + rc;
3695 ul = (up[i] ^ ux) + uc;
3698 rl = ( (ul & vx) ^ rx) + rc;
3705 rn = mpn_normalized_size (rp, rn);
3707 r->_mp_size = rx ? -rn : rn;
3711 mpz_ior (mpz_t r, const mpz_t u, const mpz_t v)
3713 mp_size_t un, vn, rn, i;
3716 mp_limb_t ux, vx, rx;
3717 mp_limb_t uc, vc, rc;
3718 mp_limb_t ul, vl, rl;
3720 un = GMP_ABS (u->_mp_size);
3721 vn = GMP_ABS (v->_mp_size);
3724 MPZ_SRCPTR_SWAP (u, v);
3725 MP_SIZE_T_SWAP (un, vn);
3733 uc = u->_mp_size < 0;
3734 vc = v->_mp_size < 0;
3741 /* If the smaller input is negative, by sign extension higher limbs
3745 rp = MPZ_REALLOC (r, rn + (mp_size_t) rc);
3753 ul = (up[i] ^ ux) + uc;
3756 vl = (vp[i] ^ vx) + vc;
3759 rl = ( (ul | vl) ^ rx) + rc;
3768 ul = (up[i] ^ ux) + uc;
3771 rl = ( (ul | vx) ^ rx) + rc;
3778 rn = mpn_normalized_size (rp, rn);
3780 r->_mp_size = rx ? -rn : rn;
3784 mpz_xor (mpz_t r, const mpz_t u, const mpz_t v)
3786 mp_size_t un, vn, i;
3789 mp_limb_t ux, vx, rx;
3790 mp_limb_t uc, vc, rc;
3791 mp_limb_t ul, vl, rl;
3793 un = GMP_ABS (u->_mp_size);
3794 vn = GMP_ABS (v->_mp_size);
3797 MPZ_SRCPTR_SWAP (u, v);
3798 MP_SIZE_T_SWAP (un, vn);
3806 uc = u->_mp_size < 0;
3807 vc = v->_mp_size < 0;
3814 rp = MPZ_REALLOC (r, un + (mp_size_t) rc);
3822 ul = (up[i] ^ ux) + uc;
3825 vl = (vp[i] ^ vx) + vc;
3828 rl = (ul ^ vl ^ rx) + rc;
3837 ul = (up[i] ^ ux) + uc;
3840 rl = (ul ^ ux) + rc;
3847 un = mpn_normalized_size (rp, un);
3849 r->_mp_size = rx ? -un : un;
3853 gmp_popcount_limb (mp_limb_t x)
3857 /* Do 16 bits at a time, to avoid limb-sized constants. */
3858 for (c = 0; x > 0; x >>= 16)
3860 unsigned w = ((x >> 1) & 0x5555) + (x & 0x5555);
3861 w = ((w >> 2) & 0x3333) + (w & 0x3333);
3862 w = ((w >> 4) & 0x0f0f) + (w & 0x0f0f);
3863 w = (w >> 8) + (w & 0x00ff);
3870 mpn_popcount (mp_srcptr p, mp_size_t n)
3875 for (c = 0, i = 0; i < n; i++)
3876 c += gmp_popcount_limb (p[i]);
3882 mpz_popcount (const mpz_t u)
3889 return ~(mp_bitcnt_t) 0;
3891 return mpn_popcount (u->_mp_d, un);
3895 mpz_hamdist (const mpz_t u, const mpz_t v)
3897 mp_size_t un, vn, i;
3898 mp_limb_t uc, vc, ul, vl, comp;
3906 return ~(mp_bitcnt_t) 0;
3908 comp = - (uc = vc = (un < 0));
3920 MPN_SRCPTR_SWAP (up, un, vp, vn);
3922 for (i = 0, c = 0; i < vn; i++)
3924 ul = (up[i] ^ comp) + uc;
3927 vl = (vp[i] ^ comp) + vc;
3930 c += gmp_popcount_limb (ul ^ vl);
3936 ul = (up[i] ^ comp) + uc;
3939 c += gmp_popcount_limb (ul ^ comp);
3946 mpz_scan1 (const mpz_t u, mp_bitcnt_t starting_bit)
3949 mp_size_t us, un, i;
3954 i = starting_bit / GMP_LIMB_BITS;
3956 /* Past the end there's no 1 bits for u>=0, or an immediate 1 bit
3957 for u<0. Notice this test picks up any u==0 too. */
3959 return (us >= 0 ? ~(mp_bitcnt_t) 0 : starting_bit);
3965 if (starting_bit != 0)
3969 ux = mpn_zero_p (up, i);
3971 ux = - (mp_limb_t) (limb >= ux);
3974 /* Mask to 0 all bits before starting_bit, thus ignoring them. */
3975 limb &= (GMP_LIMB_MAX << (starting_bit % GMP_LIMB_BITS));
3978 return mpn_common_scan (limb, i, up, un, ux);
3982 mpz_scan0 (const mpz_t u, mp_bitcnt_t starting_bit)
3985 mp_size_t us, un, i;
3989 ux = - (mp_limb_t) (us >= 0);
3991 i = starting_bit / GMP_LIMB_BITS;
3993 /* When past end, there's an immediate 0 bit for u>=0, or no 0 bits for
3994 u<0. Notice this test picks up all cases of u==0 too. */
3996 return (ux ? starting_bit : ~(mp_bitcnt_t) 0);
4002 limb -= mpn_zero_p (up, i); /* limb = ~(~limb + zero_p) */
4004 /* Mask all bits before starting_bit, thus ignoring them. */
4005 limb &= (GMP_LIMB_MAX << (starting_bit % GMP_LIMB_BITS));
4007 return mpn_common_scan (limb, i, up, un, ux);
4011 /* MPZ base conversion. */
4014 mpz_sizeinbase (const mpz_t u, int base)
4020 struct gmp_div_inverse bi;
4024 assert (base <= 62);
4026 un = GMP_ABS (u->_mp_size);
4032 bits = (un - 1) * GMP_LIMB_BITS + mpn_limb_size_in_base_2 (up[un-1]);
4038 return (bits + 1) / 2;
4040 return (bits + 2) / 3;
4042 return (bits + 3) / 4;
4044 return (bits + 4) / 5;
4045 /* FIXME: Do something more clever for the common case of base
4049 tp = gmp_xalloc_limbs (un);
4050 mpn_copyi (tp, up, un);
4051 mpn_div_qr_1_invert (&bi, base);
4057 mpn_div_qr_1_preinv (tp, tp, un, &bi);
4058 un -= (tp[un-1] == 0);
4067 mpz_get_str (char *sp, int base, const mpz_t u)
4074 digits = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz";
4078 digits = "0123456789abcdefghijklmnopqrstuvwxyz";
4082 else if (base >= -1)
4091 sn = 1 + mpz_sizeinbase (u, base);
4093 sp = (char *) gmp_xalloc (1 + sn);
4095 un = GMP_ABS (u->_mp_size);
4106 if (u->_mp_size < 0)
4109 bits = mpn_base_power_of_two_p (base);
4112 /* Not modified in this case. */
4113 sn = i + mpn_get_str_bits ((unsigned char *) sp + i, bits, u->_mp_d, un);
4116 struct mpn_base_info info;
4119 mpn_get_base_info (&info, base);
4120 tp = gmp_xalloc_limbs (un);
4121 mpn_copyi (tp, u->_mp_d, un);
4123 sn = i + mpn_get_str_other ((unsigned char *) sp + i, base, &info, tp, un);
4128 sp[i] = digits[(unsigned char) sp[i]];
4135 mpz_set_str (mpz_t r, const char *sp, int base)
4137 unsigned bits, value_of_a;
4138 mp_size_t rn, alloc;
4144 assert (base == 0 || (base >= 2 && base <= 62));
4146 while (isspace( (unsigned char) *sp))
4149 sign = (*sp == '-');
4156 if (sp[1] == 'x' || sp[1] == 'X')
4161 else if (sp[1] == 'b' || sp[1] == 'B')
4178 dp = (unsigned char *) gmp_xalloc (strlen (sp));
4180 value_of_a = (base > 36) ? 36 : 10;
4181 for (dn = 0; *sp; sp++)
4185 if (isspace ((unsigned char) *sp))
4187 else if (*sp >= '0' && *sp <= '9')
4189 else if (*sp >= 'a' && *sp <= 'z')
4190 digit = *sp - 'a' + value_of_a;
4191 else if (*sp >= 'A' && *sp <= 'Z')
4192 digit = *sp - 'A' + 10;
4194 digit = base; /* fail */
4196 if (digit >= (unsigned) base)
4212 bits = mpn_base_power_of_two_p (base);
4216 alloc = (dn * bits + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS;
4217 rp = MPZ_REALLOC (r, alloc);
4218 rn = mpn_set_str_bits (rp, dp, dn, bits);
4222 struct mpn_base_info info;
4223 mpn_get_base_info (&info, base);
4224 alloc = (dn + info.exp - 1) / info.exp;
4225 rp = MPZ_REALLOC (r, alloc);
4226 rn = mpn_set_str_other (rp, dp, dn, base, &info);
4227 /* Normalization, needed for all-zero input. */
4229 rn -= rp[rn-1] == 0;
4231 assert (rn <= alloc);
4234 r->_mp_size = sign ? - rn : rn;
4240 mpz_init_set_str (mpz_t r, const char *sp, int base)
4243 return mpz_set_str (r, sp, base);
4247 mpz_out_str (FILE *stream, int base, const mpz_t x)
4252 str = mpz_get_str (NULL, base, x);
4254 len = fwrite (str, 1, len, stream);
4261 gmp_detect_endian (void)
4263 static const int i = 2;
4264 const unsigned char *p = (const unsigned char *) &i;
4268 /* Import and export. Does not support nails. */
4270 mpz_import (mpz_t r, size_t count, int order, size_t size, int endian,
4271 size_t nails, const void *src)
4273 const unsigned char *p;
4274 ptrdiff_t word_step;
4278 /* The current (partial) limb. */
4280 /* The number of bytes already copied to this limb (starting from
4283 /* The index where the limb should be stored, when completed. */
4287 gmp_die ("mpz_import: Nails not supported.");
4289 assert (order == 1 || order == -1);
4290 assert (endian >= -1 && endian <= 1);
4293 endian = gmp_detect_endian ();
4295 p = (unsigned char *) src;
4297 word_step = (order != endian) ? 2 * size : 0;
4299 /* Process bytes from the least significant end, so point p at the
4300 least significant word. */
4303 p += size * (count - 1);
4304 word_step = - word_step;
4307 /* And at least significant byte of that word. */
4311 rn = (size * count + sizeof(mp_limb_t) - 1) / sizeof(mp_limb_t);
4312 rp = MPZ_REALLOC (r, rn);
4314 for (limb = 0, bytes = 0, i = 0; count > 0; count--, p += word_step)
4317 for (j = 0; j < size; j++, p -= (ptrdiff_t) endian)
4319 limb |= (mp_limb_t) *p << (bytes++ * CHAR_BIT);
4320 if (bytes == sizeof(mp_limb_t))
4328 assert (i + (bytes > 0) == rn);
4332 i = mpn_normalized_size (rp, i);
4338 mpz_export (void *r, size_t *countp, int order, size_t size, int endian,
4339 size_t nails, const mpz_t u)
4345 gmp_die ("mpz_import: Nails not supported.");
4347 assert (order == 1 || order == -1);
4348 assert (endian >= -1 && endian <= 1);
4349 assert (size > 0 || u->_mp_size == 0);
4357 ptrdiff_t word_step;
4358 /* The current (partial) limb. */
4360 /* The number of bytes left to to in this limb. */
4362 /* The index where the limb was read. */
4367 /* Count bytes in top limb. */
4368 limb = u->_mp_d[un-1];
4373 k++; limb >>= CHAR_BIT;
4374 } while (limb != 0);
4376 count = (k + (un-1) * sizeof (mp_limb_t) + size - 1) / size;
4379 r = gmp_xalloc (count * size);
4382 endian = gmp_detect_endian ();
4384 p = (unsigned char *) r;
4386 word_step = (order != endian) ? 2 * size : 0;
4388 /* Process bytes from the least significant end, so point p at the
4389 least significant word. */
4392 p += size * (count - 1);
4393 word_step = - word_step;
4396 /* And at least significant byte of that word. */
4400 for (bytes = 0, i = 0, k = 0; k < count; k++, p += word_step)
4403 for (j = 0; j < size; j++, p -= (ptrdiff_t) endian)
4408 limb = u->_mp_d[i++];
4409 bytes = sizeof (mp_limb_t);
4417 assert (k == count);