1 /* mini-gmp, a minimalistic implementation of a GNU GMP subset.
3 Contributed to the GNU project by Niels Möller
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.
9 This file is part of the GNU MP Library.
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.
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.
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/. */
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. */
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. */
46 #define GMP_LIMB_BITS (sizeof(mp_limb_t) * CHAR_BIT)
48 #define GMP_LIMB_MAX (~ (mp_limb_t) 0)
49 #define GMP_LIMB_HIGHBIT ((mp_limb_t) 1 << (GMP_LIMB_BITS - 1))
51 #define GMP_HLIMB_BIT ((mp_limb_t) 1 << (GMP_LIMB_BITS / 2))
52 #define GMP_LLIMB_MASK (GMP_HLIMB_BIT - 1)
54 #define GMP_ULONG_BITS (sizeof(unsigned long) * CHAR_BIT)
55 #define GMP_ULONG_HIGHBIT ((unsigned long) 1 << (GMP_ULONG_BITS - 1))
57 #define GMP_ABS(x) ((x) >= 0 ? (x) : -(x))
58 #define GMP_NEG_CAST(T,x) (-((T)((x) + 1) - 1))
60 #define GMP_MIN(a, b) ((a) < (b) ? (a) : (b))
61 #define GMP_MAX(a, b) ((a) > (b) ? (a) : (b))
63 #define gmp_assert_nocarry(x) do { \
68 #define gmp_clz(count, x) do { \
69 mp_limb_t __clz_x = (x); \
72 (__clz_x & ((mp_limb_t) 0xff << (GMP_LIMB_BITS - 8))) == 0; \
75 for (; (__clz_x & GMP_LIMB_HIGHBIT) == 0; __clz_c++) \
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; \
87 #define gmp_add_ssaaaa(sh, sl, ah, al, bh, bl) \
91 (sh) = (ah) + (bh) + (__x < (al)); \
95 #define gmp_sub_ddmmss(sh, sl, ah, al, bh, bl) \
99 (sh) = (ah) - (bh) - ((al) < (bl)); \
103 #define gmp_umul_ppmm(w1, w0, u, v) \
105 mp_limb_t __x0, __x1, __x2, __x3; \
106 unsigned __ul, __vl, __uh, __vh; \
107 mp_limb_t __u = (u), __v = (v); \
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); \
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; \
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. */ \
124 (w1) = __x3 + (__x1 >> (GMP_LIMB_BITS / 2)); \
125 (w0) = (__x1 << (GMP_LIMB_BITS / 2)) + (__x0 & GMP_LLIMB_MASK); \
128 #define gmp_udiv_qrnnd_preinv(q, r, nh, nl, d, di) \
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 */ \
147 #define gmp_udiv_qr_3by2(q, r1, r0, n2, n1, n0, d1, d0, dinv) \
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)); \
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); \
160 /* Conditionally adjust q and the remainders */ \
161 _mask = - (mp_limb_t) ((r1) >= _q0); \
163 gmp_add_ssaaaa ((r1), (r0), (r1), (r0), _mask & (d1), _mask & (d0)); \
166 if ((r1) > (d1) || (r0) >= (d0)) \
169 gmp_sub_ddmmss ((r1), (r0), (r1), (r0), (d1), (d0)); \
175 #define MP_LIMB_T_SWAP(x, y) \
177 mp_limb_t __mp_limb_t_swap__tmp = (x); \
179 (y) = __mp_limb_t_swap__tmp; \
181 #define MP_SIZE_T_SWAP(x, y) \
183 mp_size_t __mp_size_t_swap__tmp = (x); \
185 (y) = __mp_size_t_swap__tmp; \
187 #define MP_BITCNT_T_SWAP(x,y) \
189 mp_bitcnt_t __mp_bitcnt_t_swap__tmp = (x); \
191 (y) = __mp_bitcnt_t_swap__tmp; \
193 #define MP_PTR_SWAP(x, y) \
195 mp_ptr __mp_ptr_swap__tmp = (x); \
197 (y) = __mp_ptr_swap__tmp; \
199 #define MP_SRCPTR_SWAP(x, y) \
201 mp_srcptr __mp_srcptr_swap__tmp = (x); \
203 (y) = __mp_srcptr_swap__tmp; \
206 #define MPN_PTR_SWAP(xp,xs, yp,ys) \
208 MP_PTR_SWAP (xp, yp); \
209 MP_SIZE_T_SWAP (xs, ys); \
211 #define MPN_SRCPTR_SWAP(xp,xs, yp,ys) \
213 MP_SRCPTR_SWAP (xp, yp); \
214 MP_SIZE_T_SWAP (xs, ys); \
217 #define MPZ_PTR_SWAP(x, y) \
219 mpz_ptr __mpz_ptr_swap__tmp = (x); \
221 (y) = __mpz_ptr_swap__tmp; \
223 #define MPZ_SRCPTR_SWAP(x, y) \
225 mpz_srcptr __mpz_srcptr_swap__tmp = (x); \
227 (y) = __mpz_srcptr_swap__tmp; \
231 /* Memory allocation and other helper functions. */
233 gmp_die (const char *msg)
235 fprintf (stderr, "%s\n", msg);
240 gmp_default_alloc (size_t size)
248 gmp_die("gmp_default_alloc: Virtual memory exhausted.");
254 gmp_default_realloc (void *old, size_t old_size, size_t new_size)
258 p = realloc (old, new_size);
261 gmp_die("gmp_default_realoc: Virtual memory exhausted.");
267 gmp_default_free (void *p, size_t size)
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;
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))
282 *alloc_func = gmp_allocate_func;
285 *realloc_func = gmp_reallocate_func;
288 *free_func = gmp_free_func;
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))
297 alloc_func = gmp_default_alloc;
299 realloc_func = gmp_default_realloc;
301 free_func = gmp_default_free;
303 gmp_allocate_func = alloc_func;
304 gmp_reallocate_func = realloc_func;
305 gmp_free_func = free_func;
308 #define gmp_xalloc(size) ((*gmp_allocate_func)((size)))
309 #define gmp_free(p) ((*gmp_free_func) ((p), 0))
312 gmp_xalloc_limbs (mp_size_t size)
314 return gmp_xalloc (size * sizeof (mp_limb_t));
318 gmp_xrealloc_limbs (mp_ptr old, mp_size_t size)
321 return (*gmp_reallocate_func) (old, 0, size * sizeof (mp_limb_t));
328 mpn_copyi (mp_ptr d, mp_srcptr s, mp_size_t n)
331 for (i = 0; i < n; i++)
336 mpn_copyd (mp_ptr d, mp_srcptr s, mp_size_t n)
343 mpn_cmp (mp_srcptr ap, mp_srcptr bp, mp_size_t n)
347 if (ap[n-1] < bp[n-1])
349 else if (ap[n-1] > bp[n-1])
356 mpn_cmp4 (mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn)
363 return mpn_cmp (ap, bp, an);
367 mpn_normalized_size (mp_srcptr xp, mp_size_t n)
369 for (; n > 0 && xp[n-1] == 0; n--)
374 #define mpn_zero_p(xp, n) (mpn_normalized_size ((xp), (n)) == 0)
377 mpn_add_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b)
383 for (i = 0; i < n; i++)
385 mp_limb_t r = ap[i] + b;
394 mpn_add_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n)
399 for (i = 0, cy = 0; i < n; i++)
402 a = ap[i]; b = bp[i];
413 mpn_add (mp_ptr rp, mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn)
419 cy = mpn_add_n (rp, ap, bp, bn);
421 cy = mpn_add_1 (rp + bn, ap + bn, an - bn, cy);
426 mpn_sub_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b)
432 for (i = 0; i < n; i++)
436 mp_limb_t cy = a < b;;
444 mpn_sub_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n)
449 for (i = 0, cy = 0; i < n; i++)
452 a = ap[i]; b = bp[i];
462 mpn_sub (mp_ptr rp, mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn)
468 cy = mpn_sub_n (rp, ap, bp, bn);
470 cy = mpn_sub_1 (rp + bn, ap + bn, an - bn, cy);
475 mpn_mul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl)
477 mp_limb_t ul, cl, hpl, lpl;
485 gmp_umul_ppmm (hpl, lpl, ul, vl);
488 cl = (lpl < cl) + hpl;
498 mpn_addmul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl)
500 mp_limb_t ul, cl, hpl, lpl, rl;
508 gmp_umul_ppmm (hpl, lpl, ul, vl);
511 cl = (lpl < cl) + hpl;
524 mpn_submul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl)
526 mp_limb_t ul, cl, hpl, lpl, rl;
534 gmp_umul_ppmm (hpl, lpl, ul, vl);
537 cl = (lpl < cl) + hpl;
550 mpn_mul (mp_ptr rp, mp_srcptr up, mp_size_t un, mp_srcptr vp, mp_size_t vn)
555 /* We first multiply by the low order limb. This result can be
556 stored, not added, to rp. We also avoid a loop for zeroing this
559 rp[un] = mpn_mul_1 (rp, up, un, vp[0]);
560 rp += 1, vp += 1, vn -= 1;
562 /* Now accumulate the product of up[] and the next higher limb from
567 rp[un] = mpn_addmul_1 (rp, up, un, vp[0]);
568 rp += 1, vp += 1, vn -= 1;
574 mpn_mul_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n)
576 mpn_mul (rp, ap, n, bp, n);
580 mpn_sqr (mp_ptr rp, mp_srcptr ap, mp_size_t n)
582 mpn_mul (rp, ap, n, ap, n);
586 mpn_lshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt)
588 mp_limb_t high_limb, low_limb;
595 assert (cnt < GMP_LIMB_BITS);
600 tnc = GMP_LIMB_BITS - cnt;
602 retval = low_limb >> tnc;
603 high_limb = (low_limb << cnt);
605 for (i = n - 1; i != 0; i--)
608 *--rp = high_limb | (low_limb >> tnc);
609 high_limb = (low_limb << cnt);
617 mpn_rshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt)
619 mp_limb_t high_limb, low_limb;
626 assert (cnt < GMP_LIMB_BITS);
628 tnc = GMP_LIMB_BITS - cnt;
630 retval = (high_limb << tnc);
631 low_limb = high_limb >> cnt;
633 for (i = n - 1; i != 0; i--)
636 *rp++ = low_limb | (high_limb << tnc);
637 low_limb = high_limb >> cnt;
645 /* MPN division interface. */
647 mpn_invert_3by2 (mp_limb_t u1, mp_limb_t u0)
653 /* First, do a 2/1 inverse. */
654 /* The inverse m is defined as floor( (B^2 - 1 - u1)/u1 ), so that 0 <
655 * B^2 - (B + m) u1 <= u1 */
656 assert (u1 >= GMP_LIMB_HIGHBIT);
658 ul = u1 & GMP_LLIMB_MASK;
659 uh = u1 >> (GMP_LIMB_BITS / 2);
662 r = ((~u1 - (mp_limb_t) qh * uh) << (GMP_LIMB_BITS / 2)) | GMP_LLIMB_MASK;
664 p = (mp_limb_t) qh * ul;
665 /* Adjustment steps taken from udiv_qrnnd_c */
670 if (r >= u1) /* i.e. we didn't get carry when adding to r */
679 /* Do a 3/2 division (with half limb size) */
680 p = (r >> (GMP_LIMB_BITS / 2)) * qh + r;
681 ql = (p >> (GMP_LIMB_BITS / 2)) + 1;
683 /* By the 3/2 method, we don't need the high half limb. */
684 r = (r << (GMP_LIMB_BITS / 2)) + GMP_LLIMB_MASK - ql * u1;
686 if (r >= (p << (GMP_LIMB_BITS / 2)))
691 m = ((mp_limb_t) qh << (GMP_LIMB_BITS / 2)) + ql;
713 gmp_umul_ppmm (th, tl, u0, m);
718 if (r > u1 || (r == u1 && tl > u0))
726 struct gmp_div_inverse
728 /* Normalization shift count. */
730 /* Normalized divisor (d0 unused for mpn_div_qr_1) */
732 /* Inverse, for 2/1 or 3/2. */
737 mpn_div_qr_1_invert (struct gmp_div_inverse *inv, mp_limb_t d)
744 inv->d1 = d << shift;
745 inv->di = mpn_invert_limb (inv->d1);
749 mpn_div_qr_2_invert (struct gmp_div_inverse *inv,
750 mp_limb_t d1, mp_limb_t d0)
759 d1 = (d1 << shift) | (d0 >> (GMP_LIMB_BITS - shift));
764 inv->di = mpn_invert_3by2 (d1, d0);
768 mpn_div_qr_invert (struct gmp_div_inverse *inv,
769 mp_srcptr dp, mp_size_t dn)
774 mpn_div_qr_1_invert (inv, dp[0]);
776 mpn_div_qr_2_invert (inv, dp[1], dp[0]);
789 d1 = (d1 << shift) | (d0 >> (GMP_LIMB_BITS - shift));
790 d0 = (d0 << shift) | (dp[dn-3] >> (GMP_LIMB_BITS - shift));
794 inv->di = mpn_invert_3by2 (d1, d0);
798 /* Not matching current public gmp interface, rather corresponding to
799 the sbpi1_div_* functions. */
801 mpn_div_qr_1_preinv (mp_ptr qp, mp_srcptr np, mp_size_t nn,
802 const struct gmp_div_inverse *inv)
810 tp = gmp_xalloc_limbs (nn);
811 r = mpn_lshift (tp, np, nn, inv->shift);
823 gmp_udiv_qrnnd_preinv (q, r, r, np[nn], d, di);
830 return r >> inv->shift;
834 mpn_div_qr_1 (mp_ptr qp, mp_srcptr np, mp_size_t nn, mp_limb_t d)
838 /* Special case for powers of two. */
839 if (d > 1 && (d & (d-1)) == 0)
842 mp_limb_t r = np[0] & (d-1);
845 mpn_rshift (qp, np, nn, shift);
851 struct gmp_div_inverse inv;
852 mpn_div_qr_1_invert (&inv, d);
853 return mpn_div_qr_1_preinv (qp, np, nn, &inv);
858 mpn_div_qr_2_preinv (mp_ptr qp, mp_ptr rp, mp_srcptr np, mp_size_t nn,
859 const struct gmp_div_inverse *inv)
863 mp_limb_t d1, d0, di, r1, r0;
874 tp = gmp_xalloc_limbs (nn);
875 r1 = mpn_lshift (tp, np, nn, shift);
883 for (i = nn - 2; i >= 0; i--)
887 gmp_udiv_qr_3by2 (q, r1, r0, r1, r0, n0, d1, d0, di);
895 assert ((r0 << (GMP_LIMB_BITS - shift)) == 0);
896 r0 = (r0 >> shift) | (r1 << (GMP_LIMB_BITS - shift));
908 mpn_div_qr_2 (mp_ptr qp, mp_ptr rp, mp_srcptr np, mp_size_t nn,
909 mp_limb_t d1, mp_limb_t d0)
911 struct gmp_div_inverse inv;
914 mpn_div_qr_2_invert (&inv, d1, d0);
915 mpn_div_qr_2_preinv (qp, rp, np, nn, &inv);
920 mpn_div_qr_pi1 (mp_ptr qp,
921 mp_ptr np, mp_size_t nn, mp_limb_t n1,
922 mp_srcptr dp, mp_size_t dn,
937 assert ((d1 & GMP_LIMB_HIGHBIT) != 0);
938 /* Iteration variable is the index of the q limb.
940 * We divide <n1, np[dn-1+i], np[dn-2+i], np[dn-3+i],..., np[i]>
941 * by <d1, d0, dp[dn-3], ..., dp[0] >
944 for (i = nn - dn; i >= 0; i--)
946 mp_limb_t n0 = np[dn-1+i];
948 if (n1 == d1 && n0 == d0)
951 mpn_submul_1 (np+i, dp, dn, q);
952 n1 = np[dn-1+i]; /* update n1, last loop's value will now be invalid */
956 gmp_udiv_qr_3by2 (q, n1, n0, n1, n0, np[dn-2+i], d1, d0, dinv);
958 cy = mpn_submul_1 (np + i, dp, dn-2, q);
968 n1 += d1 + mpn_add_n (np + i, np + i, dp, dn - 1);
981 mpn_div_qr_preinv (mp_ptr qp, mp_ptr np, mp_size_t nn,
982 mp_srcptr dp, mp_size_t dn,
983 const struct gmp_div_inverse *inv)
989 np[0] = mpn_div_qr_1_preinv (qp, np, nn, inv);
991 mpn_div_qr_2_preinv (qp, np, np, nn, inv);
997 assert (inv->d1 == dp[dn-1]);
998 assert (inv->d0 == dp[dn-2]);
999 assert ((inv->d1 & GMP_LIMB_HIGHBIT) != 0);
1003 nh = mpn_lshift (np, np, nn, shift);
1007 mpn_div_qr_pi1 (qp, np, nn, nh, dp, dn, inv->di);
1010 gmp_assert_nocarry (mpn_rshift (np, np, dn, shift));
1015 mpn_div_qr (mp_ptr qp, mp_ptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn)
1017 struct gmp_div_inverse inv;
1023 mpn_div_qr_invert (&inv, dp, dn);
1024 if (dn > 2 && inv.shift > 0)
1026 tp = gmp_xalloc_limbs (dn);
1027 gmp_assert_nocarry (mpn_lshift (tp, dp, dn, inv.shift));
1030 mpn_div_qr_preinv (qp, np, nn, dp, dn, &inv);
1036 /* MPN base conversion. */
1038 mpn_base_power_of_two_p (unsigned b)
1054 struct mpn_base_info
1056 /* bb is the largest power of the base which fits in one limb, and
1057 exp is the corresponding exponent. */
1063 mpn_get_base_info (struct mpn_base_info *info, mp_limb_t b)
1069 m = GMP_LIMB_MAX / b;
1070 for (exp = 1, p = b; p <= m; exp++)
1078 mpn_limb_size_in_base_2 (mp_limb_t u)
1084 return GMP_LIMB_BITS - shift;
1088 mpn_get_str_bits (unsigned char *sp, unsigned bits, mp_srcptr up, mp_size_t un)
1095 sn = ((un - 1) * GMP_LIMB_BITS + mpn_limb_size_in_base_2 (up[un-1])
1098 mask = (1U << bits) - 1;
1100 for (i = 0, j = sn, shift = 0; j-- > 0;)
1102 unsigned char digit = up[i] >> shift;
1106 if (shift >= GMP_LIMB_BITS && ++i < un)
1108 shift -= GMP_LIMB_BITS;
1109 digit |= up[i] << (bits - shift);
1111 sp[j] = digit & mask;
1116 /* We generate digits from the least significant end, and reverse at
1119 mpn_limb_get_str (unsigned char *sp, mp_limb_t w,
1120 const struct gmp_div_inverse *binv)
1123 for (i = 0; w > 0; i++)
1127 h = w >> (GMP_LIMB_BITS - binv->shift);
1128 l = w << binv->shift;
1130 gmp_udiv_qrnnd_preinv (w, r, h, l, binv->d1, binv->di);
1131 assert ( (r << (GMP_LIMB_BITS - binv->shift)) == 0);
1140 mpn_get_str_other (unsigned char *sp,
1141 int base, const struct mpn_base_info *info,
1142 mp_ptr up, mp_size_t un)
1144 struct gmp_div_inverse binv;
1148 mpn_div_qr_1_invert (&binv, base);
1154 struct gmp_div_inverse bbinv;
1155 mpn_div_qr_1_invert (&bbinv, info->bb);
1161 w = mpn_div_qr_1_preinv (up, up, un, &bbinv);
1162 un -= (up[un-1] == 0);
1163 done = mpn_limb_get_str (sp + sn, w, &binv);
1165 for (sn += done; done < info->exp; done++)
1170 sn += mpn_limb_get_str (sp + sn, up[0], &binv);
1173 for (i = 0; 2*i + 1 < sn; i++)
1175 unsigned char t = sp[i];
1176 sp[i] = sp[sn - i - 1];
1184 mpn_get_str (unsigned char *sp, int base, mp_ptr up, mp_size_t un)
1189 assert (up[un-1] > 0);
1191 bits = mpn_base_power_of_two_p (base);
1193 return mpn_get_str_bits (sp, bits, up, un);
1196 struct mpn_base_info info;
1198 mpn_get_base_info (&info, base);
1199 return mpn_get_str_other (sp, base, &info, up, un);
1204 mpn_set_str_bits (mp_ptr rp, const unsigned char *sp, size_t sn,
1211 for (j = sn, rn = 0, shift = 0; j-- > 0; )
1220 rp[rn-1] |= (mp_limb_t) sp[j] << shift;
1222 if (shift >= GMP_LIMB_BITS)
1224 shift -= GMP_LIMB_BITS;
1226 rp[rn++] = (mp_limb_t) sp[j] >> (bits - shift);
1230 rn = mpn_normalized_size (rp, rn);
1235 mpn_set_str_other (mp_ptr rp, const unsigned char *sp, size_t sn,
1236 mp_limb_t b, const struct mpn_base_info *info)
1244 first = 1 + (sn - 1) % info->exp;
1248 for (k = 1; k < first; k++)
1249 w = w * b + sp[j++];
1253 for (rn = (w > 0); j < sn;)
1258 for (k = 1; k < info->exp; k++)
1259 w = w * b + sp[j++];
1261 cy = mpn_mul_1 (rp, rp, rn, info->bb);
1262 cy += mpn_add_1 (rp, rp, rn, w);
1272 mpn_set_str (mp_ptr rp, const unsigned char *sp, size_t sn, int base)
1279 bits = mpn_base_power_of_two_p (base);
1281 return mpn_set_str_bits (rp, sp, sn, bits);
1284 struct mpn_base_info info;
1286 mpn_get_base_info (&info, base);
1287 return mpn_set_str_other (rp, sp, sn, base, &info);
1298 r->_mp_d = gmp_xalloc_limbs (1);
1301 /* The utility of this function is a bit limited, since many functions
1302 assings the result variable using mpz_swap. */
1304 mpz_init2 (mpz_t r, mp_bitcnt_t bits)
1308 bits -= (bits != 0); /* Round down, except if 0 */
1309 rn = 1 + bits / GMP_LIMB_BITS;
1313 r->_mp_d = gmp_xalloc_limbs (rn);
1319 gmp_free (r->_mp_d);
1323 mpz_realloc (mpz_t r, mp_size_t size)
1325 size = GMP_MAX (size, 1);
1327 r->_mp_d = gmp_xrealloc_limbs (r->_mp_d, size);
1328 r->_mp_alloc = size;
1330 if (GMP_ABS (r->_mp_size) > size)
1336 /* Realloc for an mpz_t WHAT if it has less than NEEDED limbs. */
1337 #define MPZ_REALLOC(z,n) ((n) > (z)->_mp_alloc \
1338 ? mpz_realloc(z,n) \
1341 /* MPZ assignment and basic conversions. */
1343 mpz_set_si (mpz_t r, signed long int x)
1350 r->_mp_d[0] = GMP_NEG_CAST (unsigned long int, x);
1355 mpz_set_ui (mpz_t r, unsigned long int x)
1367 mpz_set (mpz_t r, const mpz_t x)
1369 /* Allow the NOP r == x */
1375 n = GMP_ABS (x->_mp_size);
1376 rp = MPZ_REALLOC (r, n);
1378 mpn_copyi (rp, x->_mp_d, n);
1379 r->_mp_size = x->_mp_size;
1384 mpz_init_set_si (mpz_t r, signed long int x)
1391 mpz_init_set_ui (mpz_t r, unsigned long int x)
1398 mpz_init_set (mpz_t r, const mpz_t x)
1405 mpz_fits_slong_p (const mpz_t u)
1407 mp_size_t us = u->_mp_size;
1412 return u->_mp_d[0] < GMP_LIMB_HIGHBIT;
1414 return u->_mp_d[0] <= GMP_LIMB_HIGHBIT;
1420 mpz_fits_ulong_p (const mpz_t u)
1422 mp_size_t us = u->_mp_size;
1424 return us == 0 || us == 1;
1428 mpz_get_si (const mpz_t u)
1430 mp_size_t us = u->_mp_size;
1433 return (long) (u->_mp_d[0] & ~GMP_LIMB_HIGHBIT);
1435 return (long) (- u->_mp_d[0] | GMP_LIMB_HIGHBIT);
1441 mpz_get_ui (const mpz_t u)
1443 return u->_mp_size == 0 ? 0 : u->_mp_d[0];
1447 mpz_size (const mpz_t u)
1449 return GMP_ABS (u->_mp_size);
1453 mpz_getlimbn (const mpz_t u, mp_size_t n)
1455 if (n >= 0 && n < GMP_ABS (u->_mp_size))
1462 /* Conversions and comparison to double. */
1464 mpz_set_d (mpz_t r, double x)
1473 /* x != x is true when x is a NaN, and x == x * 0.5 is true when x is
1474 zero or infinity. */
1475 if (x == 0.0 || x != x || x == x * 0.5)
1494 B = 2.0 * (double) GMP_LIMB_HIGHBIT;
1496 for (rn = 1; x >= B; rn++)
1499 rp = MPZ_REALLOC (r, rn);
1505 for (i = rn-1; i-- > 0; )
1514 r->_mp_size = sign ? - rn : rn;
1518 mpz_init_set_d (mpz_t r, double x)
1525 mpz_get_d (const mpz_t u)
1529 double B = 2.0 * (double) GMP_LIMB_HIGHBIT;
1531 un = GMP_ABS (u->_mp_size);
1538 x = B*x + u->_mp_d[--un];
1540 if (u->_mp_size < 0)
1547 mpz_cmpabs_d (const mpz_t x, double d)
1560 B = 2.0 * (double) GMP_LIMB_HIGHBIT;
1563 /* Scale d so it can be compared with the top limb. */
1564 for (i = 1; i < xn; i++)
1570 /* Compare floor(d) to top limb, subtract and cancel when equal. */
1571 for (i = xn; i-- > 0;)
1588 mpz_cmp_d (const mpz_t x, double d)
1590 if (x->_mp_size < 0)
1595 return -mpz_cmpabs_d (x, d);
1602 return mpz_cmpabs_d (x, d);
1607 /* MPZ comparisons and the like. */
1609 mpz_sgn (const mpz_t u)
1611 mp_size_t usize = u->_mp_size;
1622 mpz_cmp_si (const mpz_t u, long v)
1624 mp_size_t usize = u->_mp_size;
1629 return mpz_cmp_ui (u, v);
1630 else if (usize >= 0)
1632 else /* usize == -1 */
1634 mp_limb_t ul = u->_mp_d[0];
1635 if ((mp_limb_t)GMP_NEG_CAST (unsigned long int, v) < ul)
1637 else if ( (mp_limb_t)GMP_NEG_CAST (unsigned long int, v) > ul)
1644 mpz_cmp_ui (const mpz_t u, unsigned long v)
1646 mp_size_t usize = u->_mp_size;
1654 mp_limb_t ul = (usize > 0) ? u->_mp_d[0] : 0;
1664 mpz_cmp (const mpz_t a, const mpz_t b)
1666 mp_size_t asize = a->_mp_size;
1667 mp_size_t bsize = b->_mp_size;
1671 else if (asize < bsize)
1674 return mpn_cmp (a->_mp_d, b->_mp_d, asize);
1676 return -mpn_cmp (a->_mp_d, b->_mp_d, -asize);
1682 mpz_cmpabs_ui (const mpz_t u, unsigned long v)
1684 mp_size_t un = GMP_ABS (u->_mp_size);
1690 ul = (un == 1) ? u->_mp_d[0] : 0;
1701 mpz_cmpabs (const mpz_t u, const mpz_t v)
1703 return mpn_cmp4 (u->_mp_d, GMP_ABS (u->_mp_size),
1704 v->_mp_d, GMP_ABS (v->_mp_size));
1708 mpz_abs (mpz_t r, const mpz_t u)
1713 r->_mp_size = GMP_ABS (r->_mp_size);
1717 mpz_neg (mpz_t r, const mpz_t u)
1722 r->_mp_size = -r->_mp_size;
1726 mpz_swap (mpz_t u, mpz_t v)
1728 MP_SIZE_T_SWAP (u->_mp_size, v->_mp_size);
1729 MP_SIZE_T_SWAP (u->_mp_alloc, v->_mp_alloc);
1730 MP_PTR_SWAP (u->_mp_d, v->_mp_d);
1734 /* MPZ addition and subtraction */
1736 /* Adds to the absolute value. Returns new size, but doesn't store it. */
1738 mpz_abs_add_ui (mpz_t r, const mpz_t a, unsigned long b)
1744 an = GMP_ABS (a->_mp_size);
1751 rp = MPZ_REALLOC (r, an + 1);
1753 cy = mpn_add_1 (rp, a->_mp_d, an, b);
1760 /* Subtract from the absolute value. Returns new size, (or -1 on underflow),
1761 but doesn't store it. */
1763 mpz_abs_sub_ui (mpz_t r, const mpz_t a, unsigned long b)
1765 mp_size_t an = GMP_ABS (a->_mp_size);
1766 mp_ptr rp = MPZ_REALLOC (r, an);
1773 else if (an == 1 && a->_mp_d[0] < b)
1775 rp[0] = b - a->_mp_d[0];
1780 gmp_assert_nocarry (mpn_sub_1 (rp, a->_mp_d, an, b));
1781 return mpn_normalized_size (rp, an);
1786 mpz_add_ui (mpz_t r, const mpz_t a, unsigned long b)
1788 if (a->_mp_size >= 0)
1789 r->_mp_size = mpz_abs_add_ui (r, a, b);
1791 r->_mp_size = -mpz_abs_sub_ui (r, a, b);
1795 mpz_sub_ui (mpz_t r, const mpz_t a, unsigned long b)
1797 if (a->_mp_size < 0)
1798 r->_mp_size = -mpz_abs_add_ui (r, a, b);
1800 r->_mp_size = mpz_abs_sub_ui (r, a, b);
1804 mpz_ui_sub (mpz_t r, unsigned long a, const mpz_t b)
1806 if (b->_mp_size < 0)
1807 r->_mp_size = mpz_abs_add_ui (r, b, a);
1809 r->_mp_size = -mpz_abs_sub_ui (r, b, a);
1813 mpz_abs_add (mpz_t r, const mpz_t a, const mpz_t b)
1815 mp_size_t an = GMP_ABS (a->_mp_size);
1816 mp_size_t bn = GMP_ABS (b->_mp_size);
1821 rn = GMP_MAX (an, bn);
1822 rp = MPZ_REALLOC (r, rn + 1);
1824 cy = mpn_add (rp, a->_mp_d, an, b->_mp_d, bn);
1826 cy = mpn_add (rp, b->_mp_d, bn, a->_mp_d, an);
1830 return rn + (cy > 0);
1834 mpz_abs_sub (mpz_t r, const mpz_t a, const mpz_t b)
1836 mp_size_t an = GMP_ABS (a->_mp_size);
1837 mp_size_t bn = GMP_ABS (b->_mp_size);
1841 cmp = mpn_cmp4 (a->_mp_d, an, b->_mp_d, bn);
1844 rp = MPZ_REALLOC (r, an);
1845 gmp_assert_nocarry (mpn_sub (rp, a->_mp_d, an, b->_mp_d, bn));
1846 return mpn_normalized_size (rp, an);
1850 rp = MPZ_REALLOC (r, bn);
1851 gmp_assert_nocarry (mpn_sub (rp, b->_mp_d, bn, a->_mp_d, an));
1852 return -mpn_normalized_size (rp, bn);
1859 mpz_add (mpz_t r, const mpz_t a, const mpz_t b)
1863 if ( (a->_mp_size ^ b->_mp_size) >= 0)
1864 rn = mpz_abs_add (r, a, b);
1866 rn = mpz_abs_sub (r, a, b);
1868 r->_mp_size = a->_mp_size >= 0 ? rn : - rn;
1872 mpz_sub (mpz_t r, const mpz_t a, const mpz_t b)
1876 if ( (a->_mp_size ^ b->_mp_size) >= 0)
1877 rn = mpz_abs_sub (r, a, b);
1879 rn = mpz_abs_add (r, a, b);
1881 r->_mp_size = a->_mp_size >= 0 ? rn : - rn;
1885 /* MPZ multiplication */
1887 mpz_mul_si (mpz_t r, const mpz_t u, long int v)
1891 mpz_mul_ui (r, u, GMP_NEG_CAST (unsigned long int, v));
1895 mpz_mul_ui (r, u, (unsigned long int) v);
1899 mpz_mul_ui (mpz_t r, const mpz_t u, unsigned long int v)
1906 un = GMP_ABS (u->_mp_size);
1908 if (un == 0 || v == 0)
1914 mpz_init2 (t, (un + 1) * GMP_LIMB_BITS);
1917 cy = mpn_mul_1 (tp, u->_mp_d, un, v);
1920 t->_mp_size = un + (cy > 0);
1921 if (u->_mp_size < 0)
1922 t->_mp_size = - t->_mp_size;
1929 mpz_mul (mpz_t r, const mpz_t u, const mpz_t v)
1932 mp_size_t un, vn, rn;
1936 un = GMP_ABS (u->_mp_size);
1937 vn = GMP_ABS (v->_mp_size);
1939 if (un == 0 || vn == 0)
1945 sign = (u->_mp_size ^ v->_mp_size) < 0;
1947 mpz_init2 (t, (un + vn) * GMP_LIMB_BITS);
1951 mpn_mul (tp, u->_mp_d, un, v->_mp_d, vn);
1953 mpn_mul (tp, v->_mp_d, vn, u->_mp_d, un);
1956 rn -= tp[rn-1] == 0;
1958 t->_mp_size = sign ? - rn : rn;
1964 mpz_mul_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bits)
1971 un = GMP_ABS (u->_mp_size);
1978 limbs = bits / GMP_LIMB_BITS;
1979 shift = bits % GMP_LIMB_BITS;
1981 rn = un + limbs + (shift > 0);
1982 rp = MPZ_REALLOC (r, rn);
1985 mp_limb_t cy = mpn_lshift (rp + limbs, u->_mp_d, un, shift);
1990 mpn_copyd (rp + limbs, u->_mp_d, un);
1995 r->_mp_size = (u->_mp_size < 0) ? - rn : rn;
2000 enum mpz_div_round_mode { GMP_DIV_FLOOR, GMP_DIV_CEIL, GMP_DIV_TRUNC };
2002 /* Allows q or r to be zero. Returns 1 iff remainder is non-zero. */
2004 mpz_div_qr (mpz_t q, mpz_t r,
2005 const mpz_t n, const mpz_t d, enum mpz_div_round_mode mode)
2007 mp_size_t ns, ds, nn, dn, qs;
2012 gmp_die("mpz_div_qr: Divide by zero.");
2030 if (mode == GMP_DIV_CEIL && qs >= 0)
2032 /* q = 1, r = n - d */
2038 else if (mode == GMP_DIV_FLOOR && qs < 0)
2040 /* q = -1, r = n + d */
2070 mpz_init2 (tq, qn * GMP_LIMB_BITS);
2076 mpn_div_qr (qp, np, nn, d->_mp_d, dn);
2080 qn -= (qp[qn-1] == 0);
2082 tq->_mp_size = qs < 0 ? -qn : qn;
2084 rn = mpn_normalized_size (np, dn);
2085 tr->_mp_size = ns < 0 ? - rn : rn;
2087 if (mode == GMP_DIV_FLOOR && qs < 0 && rn != 0)
2090 mpz_sub_ui (tq, tq, 1);
2092 mpz_add (tr, tr, d);
2094 else if (mode == GMP_DIV_CEIL && qs >= 0 && rn != 0)
2097 mpz_add_ui (tq, tq, 1);
2099 mpz_sub (tr, tr, d);
2117 mpz_cdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d)
2119 mpz_div_qr (q, r, n, d, GMP_DIV_CEIL);
2123 mpz_fdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d)
2125 mpz_div_qr (q, r, n, d, GMP_DIV_FLOOR);
2129 mpz_tdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d)
2131 mpz_div_qr (q, r, n, d, GMP_DIV_TRUNC);
2135 mpz_cdiv_q (mpz_t q, const mpz_t n, const mpz_t d)
2137 mpz_div_qr (q, NULL, n, d, GMP_DIV_CEIL);
2141 mpz_fdiv_q (mpz_t q, const mpz_t n, const mpz_t d)
2143 mpz_div_qr (q, NULL, n, d, GMP_DIV_FLOOR);
2147 mpz_tdiv_q (mpz_t q, const mpz_t n, const mpz_t d)
2149 mpz_div_qr (q, NULL, n, d, GMP_DIV_TRUNC);
2153 mpz_cdiv_r (mpz_t r, const mpz_t n, const mpz_t d)
2155 mpz_div_qr (NULL, r, n, d, GMP_DIV_CEIL);
2159 mpz_fdiv_r (mpz_t r, const mpz_t n, const mpz_t d)
2161 mpz_div_qr (NULL, r, n, d, GMP_DIV_FLOOR);
2165 mpz_tdiv_r (mpz_t r, const mpz_t n, const mpz_t d)
2167 mpz_div_qr (NULL, r, n, d, GMP_DIV_TRUNC);
2171 mpz_mod (mpz_t r, const mpz_t n, const mpz_t d)
2173 if (d->_mp_size >= 0)
2174 mpz_div_qr (NULL, r, n, d, GMP_DIV_FLOOR);
2176 mpz_div_qr (NULL, r, n, d, GMP_DIV_CEIL);
2180 mpz_div_q_2exp (mpz_t q, const mpz_t u, mp_bitcnt_t bit_index,
2181 enum mpz_div_round_mode mode)
2194 limb_cnt = bit_index / GMP_LIMB_BITS;
2195 qn = GMP_ABS (un) - limb_cnt;
2196 bit_index %= GMP_LIMB_BITS;
2198 if (mode == ((un > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* un != 0 here. */
2199 /* Note: Below, the final indexing at limb_cnt is valid because at
2200 that point we have qn > 0. */
2202 || !mpn_zero_p (u->_mp_d, limb_cnt)
2203 || (u->_mp_d[limb_cnt]
2204 & (((mp_limb_t) 1 << bit_index) - 1)));
2213 qp = MPZ_REALLOC (q, qn);
2217 mpn_rshift (qp, u->_mp_d + limb_cnt, qn, bit_index);
2218 qn -= qp[qn - 1] == 0;
2222 mpn_copyi (qp, u->_mp_d + limb_cnt, qn);
2228 mpz_add_ui (q, q, adjust);
2234 mpz_div_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bit_index,
2235 enum mpz_div_round_mode mode)
2237 mp_size_t us, un, rn;
2242 if (us == 0 || bit_index == 0)
2247 rn = (bit_index + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS;
2250 rp = MPZ_REALLOC (r, rn);
2253 mask = GMP_LIMB_MAX >> (rn * GMP_LIMB_BITS - bit_index);
2257 /* Quotient (with truncation) is zero, and remainder is
2259 if (mode == ((us > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* us != 0 here. */
2261 /* Have to negate and sign extend. */
2265 for (cy = 1, i = 0; i < un; i++)
2267 mp_limb_t s = ~u->_mp_d[i] + cy;
2272 for (; i < rn - 1; i++)
2273 rp[i] = GMP_LIMB_MAX;
2282 mpn_copyi (rp, u->_mp_d, un);
2290 mpn_copyi (rp, u->_mp_d, rn - 1);
2292 rp[rn-1] = u->_mp_d[rn-1] & mask;
2294 if (mode == ((us > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* us != 0 here. */
2296 /* If r != 0, compute 2^{bit_count} - r. */
2299 for (i = 0; i < rn && rp[i] == 0; i++)
2303 /* r > 0, need to flip sign. */
2305 for (i++; i < rn; i++)
2310 /* us is not used for anything else, so we can modify it
2311 here to indicate flipped sign. */
2316 rn = mpn_normalized_size (rp, rn);
2317 r->_mp_size = us < 0 ? -rn : rn;
2321 mpz_cdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2323 mpz_div_q_2exp (r, u, cnt, GMP_DIV_CEIL);
2327 mpz_fdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2329 mpz_div_q_2exp (r, u, cnt, GMP_DIV_FLOOR);
2333 mpz_tdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2335 mpz_div_q_2exp (r, u, cnt, GMP_DIV_TRUNC);
2339 mpz_cdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2341 mpz_div_r_2exp (r, u, cnt, GMP_DIV_CEIL);
2345 mpz_fdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2347 mpz_div_r_2exp (r, u, cnt, GMP_DIV_FLOOR);
2351 mpz_tdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2353 mpz_div_r_2exp (r, u, cnt, GMP_DIV_TRUNC);
2357 mpz_divexact (mpz_t q, const mpz_t n, const mpz_t d)
2359 gmp_assert_nocarry (mpz_div_qr (q, NULL, n, d, GMP_DIV_TRUNC));
2363 mpz_divisible_p (const mpz_t n, const mpz_t d)
2365 return mpz_div_qr (NULL, NULL, n, d, GMP_DIV_TRUNC) == 0;
2368 static unsigned long
2369 mpz_div_qr_ui (mpz_t q, mpz_t r,
2370 const mpz_t n, unsigned long d, enum mpz_div_round_mode mode)
2389 qp = MPZ_REALLOC (q, qn);
2393 rl = mpn_div_qr_1 (qp, n->_mp_d, qn, d);
2397 rs = (ns < 0) ? -rs : rs;
2399 if (rl > 0 && ( (mode == GMP_DIV_FLOOR && ns < 0)
2400 || (mode == GMP_DIV_CEIL && ns >= 0)))
2403 gmp_assert_nocarry (mpn_add_1 (qp, qp, qn, 1));
2415 qn -= (qp[qn-1] == 0);
2416 assert (qn == 0 || qp[qn-1] > 0);
2418 q->_mp_size = (ns < 0) ? - qn : qn;
2425 mpz_cdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d)
2427 return mpz_div_qr_ui (q, r, n, d, GMP_DIV_CEIL);
2431 mpz_fdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d)
2433 return mpz_div_qr_ui (q, r, n, d, GMP_DIV_FLOOR);
2437 mpz_tdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d)
2439 return mpz_div_qr_ui (q, r, n, d, GMP_DIV_TRUNC);
2443 mpz_cdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d)
2445 return mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_CEIL);
2449 mpz_fdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d)
2451 return mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_FLOOR);
2455 mpz_tdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d)
2457 return mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_TRUNC);
2461 mpz_cdiv_r_ui (mpz_t r, const mpz_t n, unsigned long d)
2463 return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_CEIL);
2466 mpz_fdiv_r_ui (mpz_t r, const mpz_t n, unsigned long d)
2468 return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_FLOOR);
2471 mpz_tdiv_r_ui (mpz_t r, const mpz_t n, unsigned long d)
2473 return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_TRUNC);
2477 mpz_cdiv_ui (const mpz_t n, unsigned long d)
2479 return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_CEIL);
2483 mpz_fdiv_ui (const mpz_t n, unsigned long d)
2485 return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_FLOOR);
2489 mpz_tdiv_ui (const mpz_t n, unsigned long d)
2491 return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_TRUNC);
2495 mpz_mod_ui (mpz_t r, const mpz_t n, unsigned long d)
2497 return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_FLOOR);
2501 mpz_divexact_ui (mpz_t q, const mpz_t n, unsigned long d)
2503 gmp_assert_nocarry (mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_TRUNC));
2507 mpz_divisible_ui_p (const mpz_t n, unsigned long d)
2509 return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_TRUNC) == 0;
2515 mpn_gcd_11 (mp_limb_t u, mp_limb_t v)
2519 assert ( (u | v) > 0);
2526 gmp_ctz (shift, u | v);
2532 MP_LIMB_T_SWAP (u, v);
2534 while ( (v & 1) == 0)
2544 while ( (u & 1) == 0);
2551 while ( (v & 1) == 0);
2558 mpz_gcd_ui (mpz_t g, const mpz_t u, unsigned long v)
2569 un = GMP_ABS (u->_mp_size);
2571 v = mpn_gcd_11 (mpn_div_qr_1 (NULL, u->_mp_d, un, v), v);
2581 mpz_make_odd (mpz_t r, const mpz_t u)
2583 mp_size_t un, rn, i;
2587 un = GMP_ABS (u->_mp_size);
2590 for (i = 0; u->_mp_d[i] == 0; i++)
2593 gmp_ctz (shift, u->_mp_d[i]);
2596 rp = MPZ_REALLOC (r, rn);
2599 mpn_rshift (rp, u->_mp_d + i, rn, shift);
2600 rn -= (rp[rn-1] == 0);
2603 mpn_copyi (rp, u->_mp_d + i, rn);
2606 return i * GMP_LIMB_BITS + shift;
2610 mpz_gcd (mpz_t g, const mpz_t u, const mpz_t v)
2613 mp_bitcnt_t uz, vz, gz;
2615 if (u->_mp_size == 0)
2620 if (v->_mp_size == 0)
2629 uz = mpz_make_odd (tu, u);
2630 vz = mpz_make_odd (tv, v);
2631 gz = GMP_MIN (uz, vz);
2633 if (tu->_mp_size < tv->_mp_size)
2636 mpz_tdiv_r (tu, tu, tv);
2637 if (tu->_mp_size == 0)
2646 mpz_make_odd (tu, tu);
2647 c = mpz_cmp (tu, tv);
2656 if (tv->_mp_size == 1)
2658 mp_limb_t vl = tv->_mp_d[0];
2659 mp_limb_t ul = mpz_tdiv_ui (tu, vl);
2660 mpz_set_ui (g, mpn_gcd_11 (ul, vl));
2663 mpz_sub (tu, tu, tv);
2667 mpz_mul_2exp (g, g, gz);
2671 mpz_gcdext (mpz_t g, mpz_t s, mpz_t t, const mpz_t u, const mpz_t v)
2673 mpz_t tu, tv, s0, s1, t0, t1;
2674 mp_bitcnt_t uz, vz, gz;
2677 if (u->_mp_size == 0)
2679 /* g = 0 u + sgn(v) v */
2680 signed long sign = mpz_sgn (v);
2685 mpz_set_si (t, sign);
2689 if (v->_mp_size == 0)
2691 /* g = sgn(u) u + 0 v */
2692 signed long sign = mpz_sgn (u);
2695 mpz_set_si (s, sign);
2708 uz = mpz_make_odd (tu, u);
2709 vz = mpz_make_odd (tv, v);
2710 gz = GMP_MIN (uz, vz);
2715 /* Cofactors corresponding to odd gcd. gz handled later. */
2716 if (tu->_mp_size < tv->_mp_size)
2719 MPZ_SRCPTR_SWAP (u, v);
2720 MPZ_PTR_SWAP (s, t);
2721 MP_BITCNT_T_SWAP (uz, vz);
2729 * where u and v denote the inputs with common factors of two
2730 * eliminated, and det (s0, t0; s1, t1) = 2^p. Then
2732 * 2^p tu = s1 u - t1 v
2733 * 2^p tv = -s0 u + t0 v
2736 /* After initial division, tu = q tv + tu', we have
2738 * u = 2^uz (tu' + q tv)
2743 * t0 = 2^uz, t1 = 2^uz q
2747 mpz_setbit (t0, uz);
2748 mpz_tdiv_qr (t1, tu, tu, tv);
2749 mpz_mul_2exp (t1, t1, uz);
2751 mpz_setbit (s1, vz);
2754 if (tu->_mp_size > 0)
2757 shift = mpz_make_odd (tu, tu);
2758 mpz_mul_2exp (t0, t0, shift);
2759 mpz_mul_2exp (s0, s0, shift);
2765 c = mpz_cmp (tu, tv);
2773 * u = t0 tu + t1 (tv' + tu) = (t0 + t1) tu + t1 tv'
2774 * v = s0 tu + s1 (tv' + tu) = (s0 + s1) tu + s1 tv' */
2776 mpz_sub (tv, tv, tu);
2777 mpz_add (t0, t0, t1);
2778 mpz_add (s0, s0, s1);
2780 shift = mpz_make_odd (tv, tv);
2781 mpz_mul_2exp (t1, t1, shift);
2782 mpz_mul_2exp (s1, s1, shift);
2786 mpz_sub (tu, tu, tv);
2787 mpz_add (t1, t0, t1);
2788 mpz_add (s1, s0, s1);
2790 shift = mpz_make_odd (tu, tu);
2791 mpz_mul_2exp (t0, t0, shift);
2792 mpz_mul_2exp (s0, s0, shift);
2798 /* Now tv = odd part of gcd, and -s0 and t0 are corresponding
2801 mpz_mul_2exp (tv, tv, gz);
2804 /* 2^p g = s0 u + t0 v. Eliminate one factor of two at a time. To
2805 adjust cofactors, we need u / g and v / g */
2807 mpz_divexact (s1, v, tv);
2809 mpz_divexact (t1, u, tv);
2814 /* s0 u + t0 v = (s0 - v/g) u - (t0 + u/g) v */
2815 if (mpz_odd_p (s0) || mpz_odd_p (t0))
2817 mpz_sub (s0, s0, s1);
2818 mpz_add (t0, t0, t1);
2820 mpz_divexact_ui (s0, s0, 2);
2821 mpz_divexact_ui (t0, t0, 2);
2824 /* Arrange so that |s| < |u| / 2g */
2825 mpz_add (s1, s0, s1);
2826 if (mpz_cmpabs (s0, s1) > 0)
2829 mpz_sub (t0, t0, t1);
2831 if (u->_mp_size < 0)
2833 if (v->_mp_size < 0)
2851 mpz_lcm (mpz_t r, const mpz_t u, const mpz_t v)
2855 if (u->_mp_size == 0 || v->_mp_size == 0)
2864 mpz_divexact (g, u, g);
2872 mpz_lcm_ui (mpz_t r, const mpz_t u, unsigned long v)
2874 if (v == 0 || u->_mp_size == 0)
2880 v /= mpz_gcd_ui (NULL, u, v);
2881 mpz_mul_ui (r, u, v);
2887 mpz_invert (mpz_t r, const mpz_t u, const mpz_t m)
2892 if (u->_mp_size == 0 || mpz_cmpabs_ui (m, 1) <= 0)
2898 mpz_gcdext (g, tr, NULL, u, m);
2899 invertible = (mpz_cmp_ui (g, 1) == 0);
2903 if (tr->_mp_size < 0)
2905 if (m->_mp_size >= 0)
2906 mpz_add (tr, tr, m);
2908 mpz_sub (tr, tr, m);
2919 /* Higher level operations (sqrt, pow and root) */
2922 mpz_pow_ui (mpz_t r, const mpz_t b, unsigned long e)
2926 mpz_init_set_ui (tr, 1);
2928 for (bit = GMP_ULONG_HIGHBIT; bit > 0; bit >>= 1)
2930 mpz_mul (tr, tr, tr);
2932 mpz_mul (tr, tr, b);
2939 mpz_ui_pow_ui (mpz_t r, unsigned long blimb, unsigned long e)
2942 mpz_init_set_ui (b, blimb);
2943 mpz_pow_ui (r, b, e);
2948 mpz_powm (mpz_t r, const mpz_t b, const mpz_t e, const mpz_t m)
2954 struct gmp_div_inverse minv;
2958 en = GMP_ABS (e->_mp_size);
2959 mn = GMP_ABS (m->_mp_size);
2961 gmp_die ("mpz_powm: Zero modulo.");
2970 mpn_div_qr_invert (&minv, mp, mn);
2975 /* To avoid shifts, we do all our reductions, except the final
2976 one, using a *normalized* m. */
2979 tp = gmp_xalloc_limbs (mn);
2980 gmp_assert_nocarry (mpn_lshift (tp, mp, mn, shift));
2986 if (e->_mp_size < 0)
2988 if (!mpz_invert (base, b, m))
2989 gmp_die ("mpz_powm: Negative exponent and non-invertibe base.");
2996 bn = base->_mp_size;
2999 mpn_div_qr_preinv (NULL, base->_mp_d, base->_mp_size, mp, mn, &minv);
3003 /* We have reduced the absolute value. Now take care of the
3004 sign. Note that we get zero represented non-canonically as
3006 if (b->_mp_size < 0)
3008 mp_ptr bp = MPZ_REALLOC (base, mn);
3009 gmp_assert_nocarry (mpn_sub (bp, mp, mn, bp, bn));
3012 base->_mp_size = mpn_normalized_size (base->_mp_d, bn);
3014 mpz_init_set_ui (tr, 1);
3018 mp_limb_t w = e->_mp_d[en];
3021 for (bit = GMP_LIMB_HIGHBIT; bit > 0; bit >>= 1)
3023 mpz_mul (tr, tr, tr);
3025 mpz_mul (tr, tr, base);
3026 if (tr->_mp_size > mn)
3028 mpn_div_qr_preinv (NULL, tr->_mp_d, tr->_mp_size, mp, mn, &minv);
3029 tr->_mp_size = mpn_normalized_size (tr->_mp_d, mn);
3034 /* Final reduction */
3035 if (tr->_mp_size >= mn)
3038 mpn_div_qr_preinv (NULL, tr->_mp_d, tr->_mp_size, mp, mn, &minv);
3039 tr->_mp_size = mpn_normalized_size (tr->_mp_d, mn);
3050 mpz_powm_ui (mpz_t r, const mpz_t b, unsigned long elimb, const mpz_t m)
3053 mpz_init_set_ui (e, elimb);
3054 mpz_powm (r, b, e, m);
3058 /* x=trunc(y^(1/z)), r=y-x^z */
3060 mpz_rootrem (mpz_t x, mpz_t r, const mpz_t y, unsigned long z)
3065 sgn = y->_mp_size < 0;
3066 if (sgn && (z & 1) == 0)
3067 gmp_die ("mpz_rootrem: Negative argument, with even root.");
3069 gmp_die ("mpz_rootrem: Zeroth root.");
3071 if (mpz_cmpabs_ui (y, 1) <= 0) {
3080 mpz_setbit (t, mpz_sizeinbase (y, 2) / z + 1);
3082 if (z == 2) /* simplify sqrt loop: z-1 == 1 */
3084 mpz_swap (u, t); /* u = x */
3085 mpz_tdiv_q (t, y, u); /* t = y/x */
3086 mpz_add (t, t, u); /* t = y/x + x */
3087 mpz_tdiv_q_2exp (t, t, 1); /* x'= (y/x + x)/2 */
3088 } while (mpz_cmpabs (t, u) < 0); /* |x'| < |x| */
3097 mpz_swap (u, t); /* u = x */
3098 mpz_pow_ui (t, u, z - 1); /* t = x^(z-1) */
3099 mpz_tdiv_q (t, y, t); /* t = y/x^(z-1) */
3100 mpz_mul_ui (v, u, z - 1); /* v = x*(z-1) */
3101 mpz_add (t, t, v); /* t = y/x^(z-1) + x*(z-1) */
3102 mpz_tdiv_q_ui (t, t, z); /* x'=(y/x^(z-1) + x*(z-1))/z */
3103 } while (mpz_cmpabs (t, u) < 0); /* |x'| < |x| */
3109 mpz_pow_ui (t, u, z);
3118 mpz_root (mpz_t x, const mpz_t y, unsigned long z)
3124 mpz_rootrem (x, r, y, z);
3125 res = r->_mp_size == 0;
3131 /* Compute s = floor(sqrt(u)) and r = u - s^2. Allows r == NULL */
3133 mpz_sqrtrem (mpz_t s, mpz_t r, const mpz_t u)
3135 mpz_rootrem (s, r, u, 2);
3139 mpz_sqrt (mpz_t s, const mpz_t u)
3141 mpz_rootrem (s, NULL, u, 2);
3148 mpz_fac_ui (mpz_t x, unsigned long n)
3156 mpz_mul_ui (x, x, n);
3160 mpz_bin_uiui (mpz_t r, unsigned long n, unsigned long k)
3171 mpz_divexact (r, r, t);
3172 mpz_fac_ui (t, n - k);
3173 mpz_divexact (r, r, t);
3178 /* Logical operations and bit manipulation. */
3180 /* Numbers are treated as if represented in two's complement (and
3181 infinitely sign extended). For a negative values we get the two's
3182 complement from -x = ~x + 1, where ~ is bitwise complementt.
3191 where yyyy is the bitwise complement of xxxx. So least significant
3192 bits, up to and including the first one bit, are unchanged, and
3193 the more significant bits are all complemented.
3195 To change a bit from zero to one in a negative number, subtract the
3196 corresponding power of two from the absolute value. This can never
3197 underflow. To change a bit from one to zero, add the corresponding
3198 power of two, and this might overflow. E.g., if x = -001111, the
3199 two's complement is 110001. Clearing the least significant bit, we
3200 get two's complement 110000, and -010000. */
3203 mpz_tstbit (const mpz_t d, mp_bitcnt_t bit_index)
3205 mp_size_t limb_index;
3214 limb_index = bit_index / GMP_LIMB_BITS;
3215 if (limb_index >= dn)
3218 shift = bit_index % GMP_LIMB_BITS;
3219 w = d->_mp_d[limb_index];
3220 bit = (w >> shift) & 1;
3224 /* d < 0. Check if any of the bits below is set: If so, our bit
3225 must be complemented. */
3226 if (shift > 0 && (w << (GMP_LIMB_BITS - shift)) > 0)
3228 while (limb_index-- > 0)
3229 if (d->_mp_d[limb_index] > 0)
3236 mpz_abs_add_bit (mpz_t d, mp_bitcnt_t bit_index)
3238 mp_size_t dn, limb_index;
3242 dn = GMP_ABS (d->_mp_size);
3244 limb_index = bit_index / GMP_LIMB_BITS;
3245 bit = (mp_limb_t) 1 << (bit_index % GMP_LIMB_BITS);
3247 if (limb_index >= dn)
3250 /* The bit should be set outside of the end of the number.
3251 We have to increase the size of the number. */
3252 dp = MPZ_REALLOC (d, limb_index + 1);
3254 dp[limb_index] = bit;
3255 for (i = dn; i < limb_index; i++)
3257 dn = limb_index + 1;
3265 cy = mpn_add_1 (dp + limb_index, dp + limb_index, dn - limb_index, bit);
3268 dp = MPZ_REALLOC (d, dn + 1);
3273 d->_mp_size = (d->_mp_size < 0) ? - dn : dn;
3277 mpz_abs_sub_bit (mpz_t d, mp_bitcnt_t bit_index)
3279 mp_size_t dn, limb_index;
3283 dn = GMP_ABS (d->_mp_size);
3286 limb_index = bit_index / GMP_LIMB_BITS;
3287 bit = (mp_limb_t) 1 << (bit_index % GMP_LIMB_BITS);
3289 assert (limb_index < dn);
3291 gmp_assert_nocarry (mpn_sub_1 (dp + limb_index, dp + limb_index,
3292 dn - limb_index, bit));
3293 dn -= (dp[dn-1] == 0);
3294 d->_mp_size = (d->_mp_size < 0) ? - dn : dn;
3298 mpz_setbit (mpz_t d, mp_bitcnt_t bit_index)
3300 if (!mpz_tstbit (d, bit_index))
3302 if (d->_mp_size >= 0)
3303 mpz_abs_add_bit (d, bit_index);
3305 mpz_abs_sub_bit (d, bit_index);
3310 mpz_clrbit (mpz_t d, mp_bitcnt_t bit_index)
3312 if (mpz_tstbit (d, bit_index))
3314 if (d->_mp_size >= 0)
3315 mpz_abs_sub_bit (d, bit_index);
3317 mpz_abs_add_bit (d, bit_index);
3322 mpz_combit (mpz_t d, mp_bitcnt_t bit_index)
3324 if (mpz_tstbit (d, bit_index) ^ (d->_mp_size < 0))
3325 mpz_abs_sub_bit (d, bit_index);
3327 mpz_abs_add_bit (d, bit_index);
3331 mpz_com (mpz_t r, const mpz_t u)
3334 mpz_sub_ui (r, r, 1);
3338 mpz_and (mpz_t r, const mpz_t u, const mpz_t v)
3340 mp_size_t un, vn, rn, i;
3343 mp_limb_t ux, vx, rx;
3344 mp_limb_t uc, vc, rc;
3345 mp_limb_t ul, vl, rl;
3347 un = GMP_ABS (u->_mp_size);
3348 vn = GMP_ABS (v->_mp_size);
3351 MPZ_SRCPTR_SWAP (u, v);
3352 MP_SIZE_T_SWAP (un, vn);
3360 uc = u->_mp_size < 0;
3361 vc = v->_mp_size < 0;
3368 /* If the smaller input is positive, higher limbs don't matter. */
3371 rp = MPZ_REALLOC (r, rn + rc);
3376 for (i = 0; i < vn; i++)
3378 ul = (up[i] ^ ux) + uc;
3381 vl = (vp[i] ^ vx) + vc;
3384 rl = ( (ul & vl) ^ rx) + rc;
3392 ul = (up[i] ^ ux) + uc;
3395 rl = ( (ul & vx) ^ rx) + rc;
3402 rn = mpn_normalized_size (rp, rn);
3404 r->_mp_size = rx ? -rn : rn;
3408 mpz_ior (mpz_t r, const mpz_t u, const mpz_t v)
3410 mp_size_t un, vn, rn, i;
3413 mp_limb_t ux, vx, rx;
3414 mp_limb_t uc, vc, rc;
3415 mp_limb_t ul, vl, rl;
3417 un = GMP_ABS (u->_mp_size);
3418 vn = GMP_ABS (v->_mp_size);
3421 MPZ_SRCPTR_SWAP (u, v);
3422 MP_SIZE_T_SWAP (un, vn);
3430 uc = u->_mp_size < 0;
3431 vc = v->_mp_size < 0;
3438 /* If the smaller input is negative, by sign extension higher limbs
3442 rp = MPZ_REALLOC (r, rn + rc);
3447 for (i = 0; i < vn; i++)
3449 ul = (up[i] ^ ux) + uc;
3452 vl = (vp[i] ^ vx) + vc;
3455 rl = ( (ul | vl) ^ rx) + rc;
3463 ul = (up[i] ^ ux) + uc;
3466 rl = ( (ul | vx) ^ rx) + rc;
3473 rn = mpn_normalized_size (rp, rn);
3475 r->_mp_size = rx ? -rn : rn;
3479 mpz_xor (mpz_t r, const mpz_t u, const mpz_t v)
3481 mp_size_t un, vn, i;
3484 mp_limb_t ux, vx, rx;
3485 mp_limb_t uc, vc, rc;
3486 mp_limb_t ul, vl, rl;
3488 un = GMP_ABS (u->_mp_size);
3489 vn = GMP_ABS (v->_mp_size);
3492 MPZ_SRCPTR_SWAP (u, v);
3493 MP_SIZE_T_SWAP (un, vn);
3501 uc = u->_mp_size < 0;
3502 vc = v->_mp_size < 0;
3509 rp = MPZ_REALLOC (r, un + rc);
3514 for (i = 0; i < vn; i++)
3516 ul = (up[i] ^ ux) + uc;
3519 vl = (vp[i] ^ vx) + vc;
3522 rl = (ul ^ vl ^ rx) + rc;
3530 ul = (up[i] ^ ux) + uc;
3533 rl = (ul ^ ux) + rc;
3540 un = mpn_normalized_size (rp, un);
3542 r->_mp_size = rx ? -un : un;
3546 gmp_popcount_limb (mp_limb_t x)
3550 /* Do 16 bits at a time, to avoid limb-sized constants. */
3551 for (c = 0; x > 0; x >>= 16)
3553 unsigned w = ((x >> 1) & 0x5555) + (x & 0x5555);
3554 w = ((w >> 2) & 0x3333) + (w & 0x3333);
3555 w = ((w >> 4) & 0x0f0f) + (w & 0x0f0f);
3556 w = (w >> 8) + (w & 0x00ff);
3563 mpz_popcount (const mpz_t u)
3571 return ~(mp_bitcnt_t) 0;
3573 for (c = 0, i = 0; i < un; i++)
3574 c += gmp_popcount_limb (u->_mp_d[i]);
3580 mpz_hamdist (const mpz_t u, const mpz_t v)
3582 mp_size_t un, vn, i;
3583 mp_limb_t uc, vc, ul, vl, comp;
3591 return ~(mp_bitcnt_t) 0;
3599 comp = - (mp_limb_t) 1;
3608 MPN_SRCPTR_SWAP (up, un, vp, vn);
3610 for (i = 0, c = 0; i < vn; i++)
3612 ul = (up[i] ^ comp) + uc;
3615 vl = (vp[i] ^ comp) + vc;
3618 c += gmp_popcount_limb (ul ^ vl);
3624 ul = (up[i] ^ comp) + uc;
3627 c += gmp_popcount_limb (ul ^ comp);
3634 mpz_scan1 (const mpz_t u, mp_bitcnt_t starting_bit)
3637 mp_size_t us, un, i;
3638 mp_limb_t limb, ux, uc;
3644 i = starting_bit / GMP_LIMB_BITS;
3646 /* Past the end there's no 1 bits for u>=0, or an immediate 1 bit
3647 for u<0. Notice this test picks up any u==0 too. */
3649 return (us >= 0 ? ~(mp_bitcnt_t) 0 : starting_bit);
3654 uc = mpn_zero_p (up, i);
3659 limb = (ux ^ up[i]) + uc;
3662 /* Mask to 0 all bits before starting_bit, thus ignoring them. */
3663 limb &= (GMP_LIMB_MAX << (starting_bit % GMP_LIMB_BITS));
3671 /* For the u > 0 case, this can happen only for the first
3672 masked limb. For the u < 0 case, it happens when the
3673 highest limbs of the absolute value are all ones. */
3674 return (us >= 0 ? ~(mp_bitcnt_t) 0 : un * GMP_LIMB_BITS);
3676 limb = (ux ^ up[i]) + uc;
3679 gmp_ctz (cnt, limb);
3680 return (mp_bitcnt_t) i * GMP_LIMB_BITS + cnt;
3684 mpz_scan0 (const mpz_t u, mp_bitcnt_t starting_bit)
3687 mp_size_t us, un, i;
3688 mp_limb_t limb, ux, uc;
3694 i = starting_bit / GMP_LIMB_BITS;
3696 /* When past end, there's an immediate 0 bit for u>=0, or no 0 bits for
3697 u<0. Notice this test picks up all cases of u==0 too. */
3699 return (us >= 0 ? starting_bit : ~(mp_bitcnt_t) 0);
3704 uc = mpn_zero_p (up, i);
3709 limb = (ux ^ up[i]) + uc;
3712 /* Mask to 1 all bits before starting_bit, thus ignoring them. */
3713 limb |= ((mp_limb_t) 1 << (starting_bit % GMP_LIMB_BITS)) - 1;
3715 while (limb == GMP_LIMB_MAX)
3721 return (us >= 0 ? un * GMP_LIMB_BITS : ~(mp_bitcnt_t) 0);
3723 limb = (ux ^ up[i]) + uc;
3726 gmp_ctz (cnt, ~limb);
3727 return (mp_bitcnt_t) i * GMP_LIMB_BITS + cnt;
3731 /* MPZ base conversion. */
3734 mpz_sizeinbase (const mpz_t u, int base)
3740 struct gmp_div_inverse bi;
3744 assert (base <= 36);
3746 un = GMP_ABS (u->_mp_size);
3752 bits = (un - 1) * GMP_LIMB_BITS + mpn_limb_size_in_base_2 (up[un-1]);
3758 return (bits + 1) / 2;
3760 return (bits + 2) / 3;
3762 return (bits + 3) / 4;
3764 return (bits + 4) / 5;
3765 /* FIXME: Do something more clever for the common case of base
3769 tp = gmp_xalloc_limbs (un);
3770 mpn_copyi (tp, up, un);
3771 mpn_div_qr_1_invert (&bi, base);
3773 for (ndigits = 0; un > 0; ndigits++)
3775 mpn_div_qr_1_preinv (tp, tp, un, &bi);
3776 un -= (tp[un-1] == 0);
3783 mpz_get_str (char *sp, int base, const mpz_t u)
3792 digits = "0123456789abcdefghijklmnopqrstuvwxyz";
3797 digits = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
3804 sn = 1 + mpz_sizeinbase (u, base);
3806 sp = gmp_xalloc (1 + sn);
3808 un = GMP_ABS (u->_mp_size);
3819 if (u->_mp_size < 0)
3822 bits = mpn_base_power_of_two_p (base);
3825 /* Not modified in this case. */
3826 sn = i + mpn_get_str_bits ((unsigned char *) sp + i, bits, u->_mp_d, un);
3829 struct mpn_base_info info;
3832 mpn_get_base_info (&info, base);
3833 tp = gmp_xalloc_limbs (un);
3834 mpn_copyi (tp, u->_mp_d, un);
3836 sn = i + mpn_get_str_other ((unsigned char *) sp + i, base, &info, tp, un);
3841 sp[i] = digits[(unsigned char) sp[i]];
3848 mpz_set_str (mpz_t r, const char *sp, int base)
3851 mp_size_t rn, alloc;
3858 assert (base == 0 || (base >= 2 && base <= 36));
3860 while (isspace( (unsigned char) *sp))
3876 if (*sp == 'x' || *sp == 'X')
3881 else if (*sp == 'b' || *sp == 'B')
3894 dp = gmp_xalloc (sn + (sn == 0));
3896 for (dn = 0; *sp; sp++)
3900 if (isspace ((unsigned char) *sp))
3902 if (*sp >= '0' && *sp <= '9')
3904 else if (*sp >= 'a' && *sp <= 'z')
3905 digit = *sp - 'a' + 10;
3906 else if (*sp >= 'A' && *sp <= 'Z')
3907 digit = *sp - 'A' + 10;
3909 digit = base; /* fail */
3921 bits = mpn_base_power_of_two_p (base);
3925 alloc = (sn * bits + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS;
3926 rp = MPZ_REALLOC (r, alloc);
3927 rn = mpn_set_str_bits (rp, dp, dn, bits);
3931 struct mpn_base_info info;
3932 mpn_get_base_info (&info, base);
3933 alloc = (sn + info.exp - 1) / info.exp;
3934 rp = MPZ_REALLOC (r, alloc);
3935 rn = mpn_set_str_other (rp, dp, dn, base, &info);
3937 assert (rn <= alloc);
3940 r->_mp_size = sign ? - rn : rn;
3946 mpz_init_set_str (mpz_t r, const char *sp, int base)
3949 return mpz_set_str (r, sp, base);
3953 mpz_out_str (FILE *stream, int base, const mpz_t x)
3958 str = mpz_get_str (NULL, base, x);
3960 len = fwrite (str, 1, len, stream);
3967 gmp_detect_endian (void)
3969 static const int i = 1;
3970 const unsigned char *p = (const unsigned char *) &i;
3979 /* Import and export. Does not support nails. */
3981 mpz_import (mpz_t r, size_t count, int order, size_t size, int endian,
3982 size_t nails, const void *src)
3984 const unsigned char *p;
3985 ptrdiff_t word_step;
3989 /* The current (partial) limb. */
3991 /* The number of bytes already copied to this limb (starting from
3994 /* The index where the limb should be stored, when completed. */
3998 gmp_die ("mpz_import: Nails not supported.");
4000 assert (order == 1 || order == -1);
4001 assert (endian >= -1 && endian <= 1);
4004 endian = gmp_detect_endian ();
4006 p = (unsigned char *) src;
4008 word_step = (order != endian) ? 2 * size : 0;
4010 /* Process bytes from the least significant end, so point p at the
4011 least significant word. */
4014 p += size * (count - 1);
4015 word_step = - word_step;
4018 /* And at least significant byte of that word. */
4022 rn = (size * count + sizeof(mp_limb_t) - 1) / sizeof(mp_limb_t);
4023 rp = MPZ_REALLOC (r, rn);
4025 for (limb = 0, bytes = 0, i = 0; count > 0; count--, p += word_step)
4028 for (j = 0; j < size; j++, p -= (ptrdiff_t) endian)
4030 limb |= (mp_limb_t) *p << (bytes++ * CHAR_BIT);
4031 if (bytes == sizeof(mp_limb_t))
4043 r->_mp_size = mpn_normalized_size (rp, i);
4047 mpz_export (void *r, size_t *countp, int order, size_t size, int endian,
4048 size_t nails, const mpz_t u)
4051 ptrdiff_t word_step;
4055 /* The current (partial) limb. */
4057 /* The number of bytes left to to in this limb. */
4059 /* The index where the limb was read. */
4063 gmp_die ("mpz_import: Nails not supported.");
4065 assert (order == 1 || order == -1);
4066 assert (endian >= -1 && endian <= 1);
4067 assert (size > 0 || u->_mp_size == 0);
4069 un = GMP_ABS (u->_mp_size);
4077 /* Count bytes in top limb. */
4078 for (limb = u->_mp_d[un-1], k = 0; limb > 0; k++, limb >>= CHAR_BIT)
4083 count = (k + (un-1) * sizeof (mp_limb_t) + size - 1) / size;
4086 r = gmp_xalloc (count * size);
4089 endian = gmp_detect_endian ();
4091 p = (unsigned char *) r;
4093 word_step = (order != endian) ? 2 * size : 0;
4095 /* Process bytes from the least significant end, so point p at the
4096 least significant word. */
4099 p += size * (count - 1);
4100 word_step = - word_step;
4103 /* And at least significant byte of that word. */
4107 for (bytes = 0, i = 0, k = 0; k < count; k++, p += word_step)
4110 for (j = 0; j < size; j++, p -= (ptrdiff_t) endian)
4115 limb = u->_mp_d[i++];
4116 bytes = sizeof (mp_limb_t);
4124 assert (k == count);