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-2019 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 the GNU Lesser General Public License as published by
11 the Free Software Foundation; either version 3 of the License, or (at your
12 option) any later version.
14 The GNU MP Library is distributed in the hope that it will be useful, but
15 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
16 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
17 License for more details.
19 You should have received a copy of the GNU Lesser General Public License
20 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
22 /* NOTE: All functions in this file which are not declared in
23 mini-gmp.h are internal, and are not intended to be compatible
24 neither with GMP nor with future versions of mini-gmp. */
26 /* Much of the material copied from GMP files, including: gmp-impl.h,
27 longlong.h, mpn/generic/add_n.c, mpn/generic/addmul_1.c,
28 mpn/generic/lshift.c, mpn/generic/mul_1.c,
29 mpn/generic/mul_basecase.c, mpn/generic/rshift.c,
30 mpn/generic/sbpi1_div_qr.c, mpn/generic/sub_n.c,
31 mpn/generic/submul_1.c. */
42 #if !defined(MINI_GMP_DONT_USE_FLOAT_H)
48 #define GMP_LIMB_BITS (sizeof(mp_limb_t) * CHAR_BIT)
50 #define GMP_LIMB_MAX ((mp_limb_t) ~ (mp_limb_t) 0)
51 #define GMP_LIMB_HIGHBIT ((mp_limb_t) 1 << (GMP_LIMB_BITS - 1))
53 #define GMP_HLIMB_BIT ((mp_limb_t) 1 << (GMP_LIMB_BITS / 2))
54 #define GMP_LLIMB_MASK (GMP_HLIMB_BIT - 1)
56 #define GMP_ULONG_BITS (sizeof(unsigned long) * CHAR_BIT)
57 #define GMP_ULONG_HIGHBIT ((unsigned long) 1 << (GMP_ULONG_BITS - 1))
59 #define GMP_ABS(x) ((x) >= 0 ? (x) : -(x))
60 #define GMP_NEG_CAST(T,x) (-((T)((x) + 1) - 1))
62 #define GMP_MIN(a, b) ((a) < (b) ? (a) : (b))
63 #define GMP_MAX(a, b) ((a) > (b) ? (a) : (b))
65 #define GMP_CMP(a,b) (((a) > (b)) - ((a) < (b)))
67 #if defined(DBL_MANT_DIG) && FLT_RADIX == 2
68 #define GMP_DBL_MANT_BITS DBL_MANT_DIG
70 #define GMP_DBL_MANT_BITS (53)
73 /* Return non-zero if xp,xsize and yp,ysize overlap.
74 If xp+xsize<=yp there's no overlap, or if yp+ysize<=xp there's no
75 overlap. If both these are false, there's an overlap. */
76 #define GMP_MPN_OVERLAP_P(xp, xsize, yp, ysize) \
77 ((xp) + (xsize) > (yp) && (yp) + (ysize) > (xp))
79 #define gmp_assert_nocarry(x) do { \
80 mp_limb_t __cy = (x); \
84 #define gmp_clz(count, x) do { \
85 mp_limb_t __clz_x = (x); \
86 unsigned __clz_c = 0; \
87 int LOCAL_SHIFT_BITS = 8; \
88 if (GMP_LIMB_BITS > LOCAL_SHIFT_BITS) \
90 (__clz_x & ((mp_limb_t) 0xff << (GMP_LIMB_BITS - 8))) == 0; \
92 { __clz_x <<= LOCAL_SHIFT_BITS; } \
93 for (; (__clz_x & GMP_LIMB_HIGHBIT) == 0; __clz_c++) \
98 #define gmp_ctz(count, x) do { \
99 mp_limb_t __ctz_x = (x); \
100 unsigned __ctz_c = 0; \
101 gmp_clz (__ctz_c, __ctz_x & - __ctz_x); \
102 (count) = GMP_LIMB_BITS - 1 - __ctz_c; \
105 #define gmp_add_ssaaaa(sh, sl, ah, al, bh, bl) \
109 (sh) = (ah) + (bh) + (__x < (al)); \
113 #define gmp_sub_ddmmss(sh, sl, ah, al, bh, bl) \
117 (sh) = (ah) - (bh) - ((al) < (bl)); \
121 #define gmp_umul_ppmm(w1, w0, u, v) \
123 int LOCAL_GMP_LIMB_BITS = GMP_LIMB_BITS; \
124 if (sizeof(unsigned int) * CHAR_BIT >= 2 * GMP_LIMB_BITS) \
126 unsigned int __ww = (unsigned int) (u) * (v); \
127 w0 = (mp_limb_t) __ww; \
128 w1 = (mp_limb_t) (__ww >> LOCAL_GMP_LIMB_BITS); \
130 else if (GMP_ULONG_BITS >= 2 * GMP_LIMB_BITS) \
132 unsigned long int __ww = (unsigned long int) (u) * (v); \
133 w0 = (mp_limb_t) __ww; \
134 w1 = (mp_limb_t) (__ww >> LOCAL_GMP_LIMB_BITS); \
137 mp_limb_t __x0, __x1, __x2, __x3; \
138 unsigned __ul, __vl, __uh, __vh; \
139 mp_limb_t __u = (u), __v = (v); \
141 __ul = __u & GMP_LLIMB_MASK; \
142 __uh = __u >> (GMP_LIMB_BITS / 2); \
143 __vl = __v & GMP_LLIMB_MASK; \
144 __vh = __v >> (GMP_LIMB_BITS / 2); \
146 __x0 = (mp_limb_t) __ul * __vl; \
147 __x1 = (mp_limb_t) __ul * __vh; \
148 __x2 = (mp_limb_t) __uh * __vl; \
149 __x3 = (mp_limb_t) __uh * __vh; \
151 __x1 += __x0 >> (GMP_LIMB_BITS / 2);/* this can't give carry */ \
152 __x1 += __x2; /* but this indeed can */ \
153 if (__x1 < __x2) /* did we get it? */ \
154 __x3 += GMP_HLIMB_BIT; /* yes, add it in the proper pos. */ \
156 (w1) = __x3 + (__x1 >> (GMP_LIMB_BITS / 2)); \
157 (w0) = (__x1 << (GMP_LIMB_BITS / 2)) + (__x0 & GMP_LLIMB_MASK); \
161 #define gmp_udiv_qrnnd_preinv(q, r, nh, nl, d, di) \
163 mp_limb_t _qh, _ql, _r, _mask; \
164 gmp_umul_ppmm (_qh, _ql, (nh), (di)); \
165 gmp_add_ssaaaa (_qh, _ql, _qh, _ql, (nh) + 1, (nl)); \
166 _r = (nl) - _qh * (d); \
167 _mask = -(mp_limb_t) (_r > _ql); /* both > and >= are OK */ \
180 #define gmp_udiv_qr_3by2(q, r1, r0, n2, n1, n0, d1, d0, dinv) \
182 mp_limb_t _q0, _t1, _t0, _mask; \
183 gmp_umul_ppmm ((q), _q0, (n2), (dinv)); \
184 gmp_add_ssaaaa ((q), _q0, (q), _q0, (n2), (n1)); \
186 /* Compute the two most significant limbs of n - q'd */ \
187 (r1) = (n1) - (d1) * (q); \
188 gmp_sub_ddmmss ((r1), (r0), (r1), (n0), (d1), (d0)); \
189 gmp_umul_ppmm (_t1, _t0, (d0), (q)); \
190 gmp_sub_ddmmss ((r1), (r0), (r1), (r0), _t1, _t0); \
193 /* Conditionally adjust q and the remainders */ \
194 _mask = - (mp_limb_t) ((r1) >= _q0); \
196 gmp_add_ssaaaa ((r1), (r0), (r1), (r0), _mask & (d1), _mask & (d0)); \
199 if ((r1) > (d1) || (r0) >= (d0)) \
202 gmp_sub_ddmmss ((r1), (r0), (r1), (r0), (d1), (d0)); \
208 #define MP_LIMB_T_SWAP(x, y) \
210 mp_limb_t __mp_limb_t_swap__tmp = (x); \
212 (y) = __mp_limb_t_swap__tmp; \
214 #define MP_SIZE_T_SWAP(x, y) \
216 mp_size_t __mp_size_t_swap__tmp = (x); \
218 (y) = __mp_size_t_swap__tmp; \
220 #define MP_BITCNT_T_SWAP(x,y) \
222 mp_bitcnt_t __mp_bitcnt_t_swap__tmp = (x); \
224 (y) = __mp_bitcnt_t_swap__tmp; \
226 #define MP_PTR_SWAP(x, y) \
228 mp_ptr __mp_ptr_swap__tmp = (x); \
230 (y) = __mp_ptr_swap__tmp; \
232 #define MP_SRCPTR_SWAP(x, y) \
234 mp_srcptr __mp_srcptr_swap__tmp = (x); \
236 (y) = __mp_srcptr_swap__tmp; \
239 #define MPN_PTR_SWAP(xp,xs, yp,ys) \
241 MP_PTR_SWAP (xp, yp); \
242 MP_SIZE_T_SWAP (xs, ys); \
244 #define MPN_SRCPTR_SWAP(xp,xs, yp,ys) \
246 MP_SRCPTR_SWAP (xp, yp); \
247 MP_SIZE_T_SWAP (xs, ys); \
250 #define MPZ_PTR_SWAP(x, y) \
252 mpz_ptr __mpz_ptr_swap__tmp = (x); \
254 (y) = __mpz_ptr_swap__tmp; \
256 #define MPZ_SRCPTR_SWAP(x, y) \
258 mpz_srcptr __mpz_srcptr_swap__tmp = (x); \
260 (y) = __mpz_srcptr_swap__tmp; \
263 const int mp_bits_per_limb = GMP_LIMB_BITS;
266 /* Memory allocation and other helper functions. */
268 gmp_die (const char *msg)
270 fprintf (stderr, "%s\n", msg);
275 gmp_default_alloc (size_t size)
283 gmp_die("gmp_default_alloc: Virtual memory exhausted.");
289 gmp_default_realloc (void *old, size_t unused_old_size, size_t new_size)
293 p = realloc (old, new_size);
296 gmp_die("gmp_default_realloc: Virtual memory exhausted.");
302 gmp_default_free (void *p, size_t unused_size)
307 static void * (*gmp_allocate_func) (size_t) = gmp_default_alloc;
308 static void * (*gmp_reallocate_func) (void *, size_t, size_t) = gmp_default_realloc;
309 static void (*gmp_free_func) (void *, size_t) = gmp_default_free;
312 mp_get_memory_functions (void *(**alloc_func) (size_t),
313 void *(**realloc_func) (void *, size_t, size_t),
314 void (**free_func) (void *, size_t))
317 *alloc_func = gmp_allocate_func;
320 *realloc_func = gmp_reallocate_func;
323 *free_func = gmp_free_func;
327 mp_set_memory_functions (void *(*alloc_func) (size_t),
328 void *(*realloc_func) (void *, size_t, size_t),
329 void (*free_func) (void *, size_t))
332 alloc_func = gmp_default_alloc;
334 realloc_func = gmp_default_realloc;
336 free_func = gmp_default_free;
338 gmp_allocate_func = alloc_func;
339 gmp_reallocate_func = realloc_func;
340 gmp_free_func = free_func;
343 #define gmp_xalloc(size) ((*gmp_allocate_func)((size)))
344 #define gmp_free(p) ((*gmp_free_func) ((p), 0))
347 gmp_xalloc_limbs (mp_size_t size)
349 return (mp_ptr) gmp_xalloc (size * sizeof (mp_limb_t));
353 gmp_xrealloc_limbs (mp_ptr old, mp_size_t size)
356 return (mp_ptr) (*gmp_reallocate_func) (old, 0, size * sizeof (mp_limb_t));
363 mpn_copyi (mp_ptr d, mp_srcptr s, mp_size_t n)
366 for (i = 0; i < n; i++)
371 mpn_copyd (mp_ptr d, mp_srcptr s, mp_size_t n)
378 mpn_cmp (mp_srcptr ap, mp_srcptr bp, mp_size_t n)
383 return ap[n] > bp[n] ? 1 : -1;
389 mpn_cmp4 (mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn)
392 return an < bn ? -1 : 1;
394 return mpn_cmp (ap, bp, an);
398 mpn_normalized_size (mp_srcptr xp, mp_size_t n)
400 while (n > 0 && xp[n-1] == 0)
406 mpn_zero_p(mp_srcptr rp, mp_size_t n)
408 return mpn_normalized_size (rp, n) == 0;
412 mpn_zero (mp_ptr rp, mp_size_t n)
419 mpn_add_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b)
427 mp_limb_t r = ap[i] + b;
438 mpn_add_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n)
443 for (i = 0, cy = 0; i < n; i++)
446 a = ap[i]; b = bp[i];
457 mpn_add (mp_ptr rp, mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn)
463 cy = mpn_add_n (rp, ap, bp, bn);
465 cy = mpn_add_1 (rp + bn, ap + bn, an - bn, cy);
470 mpn_sub_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b)
481 mp_limb_t cy = a < b;
491 mpn_sub_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n)
496 for (i = 0, cy = 0; i < n; i++)
499 a = ap[i]; b = bp[i];
509 mpn_sub (mp_ptr rp, mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn)
515 cy = mpn_sub_n (rp, ap, bp, bn);
517 cy = mpn_sub_1 (rp + bn, ap + bn, an - bn, cy);
522 mpn_mul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl)
524 mp_limb_t ul, cl, hpl, lpl;
532 gmp_umul_ppmm (hpl, lpl, ul, vl);
535 cl = (lpl < cl) + hpl;
545 mpn_addmul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl)
547 mp_limb_t ul, cl, hpl, lpl, rl;
555 gmp_umul_ppmm (hpl, lpl, ul, vl);
558 cl = (lpl < cl) + hpl;
571 mpn_submul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl)
573 mp_limb_t ul, cl, hpl, lpl, rl;
581 gmp_umul_ppmm (hpl, lpl, ul, vl);
584 cl = (lpl < cl) + hpl;
597 mpn_mul (mp_ptr rp, mp_srcptr up, mp_size_t un, mp_srcptr vp, mp_size_t vn)
601 assert (!GMP_MPN_OVERLAP_P(rp, un + vn, up, un));
602 assert (!GMP_MPN_OVERLAP_P(rp, un + vn, vp, vn));
604 /* We first multiply by the low order limb. This result can be
605 stored, not added, to rp. We also avoid a loop for zeroing this
608 rp[un] = mpn_mul_1 (rp, up, un, vp[0]);
610 /* Now accumulate the product of up[] and the next higher limb from
616 rp[un] = mpn_addmul_1 (rp, up, un, vp[0]);
622 mpn_mul_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n)
624 mpn_mul (rp, ap, n, bp, n);
628 mpn_sqr (mp_ptr rp, mp_srcptr ap, mp_size_t n)
630 mpn_mul (rp, ap, n, ap, n);
634 mpn_lshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt)
636 mp_limb_t high_limb, low_limb;
642 assert (cnt < GMP_LIMB_BITS);
647 tnc = GMP_LIMB_BITS - cnt;
649 retval = low_limb >> tnc;
650 high_limb = (low_limb << cnt);
655 *--rp = high_limb | (low_limb >> tnc);
656 high_limb = (low_limb << cnt);
664 mpn_rshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt)
666 mp_limb_t high_limb, low_limb;
672 assert (cnt < GMP_LIMB_BITS);
674 tnc = GMP_LIMB_BITS - cnt;
676 retval = (high_limb << tnc);
677 low_limb = high_limb >> cnt;
682 *rp++ = low_limb | (high_limb << tnc);
683 low_limb = high_limb >> cnt;
691 mpn_common_scan (mp_limb_t limb, mp_size_t i, mp_srcptr up, mp_size_t un,
696 assert (ux == 0 || ux == GMP_LIMB_MAX);
697 assert (0 <= i && i <= un );
703 return (ux == 0 ? ~(mp_bitcnt_t) 0 : un * GMP_LIMB_BITS);
707 return (mp_bitcnt_t) i * GMP_LIMB_BITS + cnt;
711 mpn_scan1 (mp_srcptr ptr, mp_bitcnt_t bit)
714 i = bit / GMP_LIMB_BITS;
716 return mpn_common_scan ( ptr[i] & (GMP_LIMB_MAX << (bit % GMP_LIMB_BITS)),
721 mpn_scan0 (mp_srcptr ptr, mp_bitcnt_t bit)
724 i = bit / GMP_LIMB_BITS;
726 return mpn_common_scan (~ptr[i] & (GMP_LIMB_MAX << (bit % GMP_LIMB_BITS)),
727 i, ptr, i, GMP_LIMB_MAX);
731 mpn_com (mp_ptr rp, mp_srcptr up, mp_size_t n)
738 mpn_neg (mp_ptr rp, mp_srcptr up, mp_size_t n)
748 mpn_com (++rp, ++up, --n);
753 /* MPN division interface. */
755 /* The 3/2 inverse is defined as
757 m = floor( (B^3-1) / (B u1 + u0)) - B
760 mpn_invert_3by2 (mp_limb_t u1, mp_limb_t u0)
768 /* For notation, let b denote the half-limb base, so that B = b^2.
769 Split u1 = b uh + ul. */
770 ul = u1 & GMP_LLIMB_MASK;
771 uh = u1 >> (GMP_LIMB_BITS / 2);
773 /* Approximation of the high half of quotient. Differs from the 2/1
774 inverse of the half limb uh, since we have already subtracted
776 qh = (u1 ^ GMP_LIMB_MAX) / uh;
778 /* Adjust to get a half-limb 3/2 inverse, i.e., we want
780 qh' = floor( (b^3 - 1) / u) - b = floor ((b^3 - b u - 1) / u
781 = floor( (b (~u) + b-1) / u),
785 r = b (~u) + b-1 - qh (b uh + ul)
786 = b (~u - qh uh) + b-1 - qh ul
788 Subtraction of qh ul may underflow, which implies adjustments.
789 But by normalization, 2 u >= B > qh ul, so we need to adjust by
793 r = ((~u1 - (mp_limb_t) qh * uh) << (GMP_LIMB_BITS / 2)) | GMP_LLIMB_MASK;
795 p = (mp_limb_t) qh * ul;
796 /* Adjustment steps taken from udiv_qrnnd_c */
801 if (r >= u1) /* i.e. we didn't get carry when adding to r */
810 /* Low half of the quotient is
812 ql = floor ( (b r + b-1) / u1).
814 This is a 3/2 division (on half-limbs), for which qh is a
817 p = (r >> (GMP_LIMB_BITS / 2)) * qh + r;
818 /* Unlike full-limb 3/2, we can add 1 without overflow. For this to
819 work, it is essential that ql is a full mp_limb_t. */
820 ql = (p >> (GMP_LIMB_BITS / 2)) + 1;
822 /* By the 3/2 trick, we don't need the high half limb. */
823 r = (r << (GMP_LIMB_BITS / 2)) + GMP_LLIMB_MASK - ql * u1;
825 if (r >= (GMP_LIMB_MAX & (p << (GMP_LIMB_BITS / 2))))
830 m = ((mp_limb_t) qh << (GMP_LIMB_BITS / 2)) + ql;
838 /* Now m is the 2/1 inverse of u1. If u0 > 0, adjust it to become a
855 gmp_umul_ppmm (th, tl, u0, m);
860 m -= ((r > u1) | ((r == u1) & (tl > u0)));
867 struct gmp_div_inverse
869 /* Normalization shift count. */
871 /* Normalized divisor (d0 unused for mpn_div_qr_1) */
873 /* Inverse, for 2/1 or 3/2. */
878 mpn_div_qr_1_invert (struct gmp_div_inverse *inv, mp_limb_t d)
885 inv->d1 = d << shift;
886 inv->di = mpn_invert_limb (inv->d1);
890 mpn_div_qr_2_invert (struct gmp_div_inverse *inv,
891 mp_limb_t d1, mp_limb_t d0)
900 d1 = (d1 << shift) | (d0 >> (GMP_LIMB_BITS - shift));
905 inv->di = mpn_invert_3by2 (d1, d0);
909 mpn_div_qr_invert (struct gmp_div_inverse *inv,
910 mp_srcptr dp, mp_size_t dn)
915 mpn_div_qr_1_invert (inv, dp[0]);
917 mpn_div_qr_2_invert (inv, dp[1], dp[0]);
930 d1 = (d1 << shift) | (d0 >> (GMP_LIMB_BITS - shift));
931 d0 = (d0 << shift) | (dp[dn-3] >> (GMP_LIMB_BITS - shift));
935 inv->di = mpn_invert_3by2 (d1, d0);
939 /* Not matching current public gmp interface, rather corresponding to
940 the sbpi1_div_* functions. */
942 mpn_div_qr_1_preinv (mp_ptr qp, mp_srcptr np, mp_size_t nn,
943 const struct gmp_div_inverse *inv)
951 /* Shift, reusing qp area if possible. In-place shift if qp == np. */
952 tp = qp ? qp : gmp_xalloc_limbs (nn);
953 r = mpn_lshift (tp, np, nn, inv->shift);
965 gmp_udiv_qrnnd_preinv (q, r, r, np[nn], d, di);
969 if ((inv->shift > 0) && (tp != qp))
972 return r >> inv->shift;
976 mpn_div_qr_2_preinv (mp_ptr qp, mp_ptr np, mp_size_t nn,
977 const struct gmp_div_inverse *inv)
981 mp_limb_t d1, d0, di, r1, r0;
990 r1 = mpn_lshift (np, np, nn, shift);
1001 gmp_udiv_qr_3by2 (q, r1, r0, r1, r0, n0, d1, d0, di);
1010 assert ((r0 & (GMP_LIMB_MAX >> (GMP_LIMB_BITS - shift))) == 0);
1011 r0 = (r0 >> shift) | (r1 << (GMP_LIMB_BITS - shift));
1020 mpn_div_qr_pi1 (mp_ptr qp,
1021 mp_ptr np, mp_size_t nn, mp_limb_t n1,
1022 mp_srcptr dp, mp_size_t dn,
1037 assert ((d1 & GMP_LIMB_HIGHBIT) != 0);
1038 /* Iteration variable is the index of the q limb.
1040 * We divide <n1, np[dn-1+i], np[dn-2+i], np[dn-3+i],..., np[i]>
1041 * by <d1, d0, dp[dn-3], ..., dp[0] >
1047 mp_limb_t n0 = np[dn-1+i];
1049 if (n1 == d1 && n0 == d0)
1052 mpn_submul_1 (np+i, dp, dn, q);
1053 n1 = np[dn-1+i]; /* update n1, last loop's value will now be invalid */
1057 gmp_udiv_qr_3by2 (q, n1, n0, n1, n0, np[dn-2+i], d1, d0, dinv);
1059 cy = mpn_submul_1 (np + i, dp, dn-2, q);
1069 n1 += d1 + mpn_add_n (np + i, np + i, dp, dn - 1);
1083 mpn_div_qr_preinv (mp_ptr qp, mp_ptr np, mp_size_t nn,
1084 mp_srcptr dp, mp_size_t dn,
1085 const struct gmp_div_inverse *inv)
1091 np[0] = mpn_div_qr_1_preinv (qp, np, nn, inv);
1093 mpn_div_qr_2_preinv (qp, np, nn, inv);
1099 assert (inv->d1 == dp[dn-1]);
1100 assert (inv->d0 == dp[dn-2]);
1101 assert ((inv->d1 & GMP_LIMB_HIGHBIT) != 0);
1105 nh = mpn_lshift (np, np, nn, shift);
1109 mpn_div_qr_pi1 (qp, np, nn, nh, dp, dn, inv->di);
1112 gmp_assert_nocarry (mpn_rshift (np, np, dn, shift));
1117 mpn_div_qr (mp_ptr qp, mp_ptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn)
1119 struct gmp_div_inverse inv;
1125 mpn_div_qr_invert (&inv, dp, dn);
1126 if (dn > 2 && inv.shift > 0)
1128 tp = gmp_xalloc_limbs (dn);
1129 gmp_assert_nocarry (mpn_lshift (tp, dp, dn, inv.shift));
1132 mpn_div_qr_preinv (qp, np, nn, dp, dn, &inv);
1138 /* MPN base conversion. */
1140 mpn_base_power_of_two_p (unsigned b)
1156 struct mpn_base_info
1158 /* bb is the largest power of the base which fits in one limb, and
1159 exp is the corresponding exponent. */
1165 mpn_get_base_info (struct mpn_base_info *info, mp_limb_t b)
1171 m = GMP_LIMB_MAX / b;
1172 for (exp = 1, p = b; p <= m; exp++)
1180 mpn_limb_size_in_base_2 (mp_limb_t u)
1186 return GMP_LIMB_BITS - shift;
1190 mpn_get_str_bits (unsigned char *sp, unsigned bits, mp_srcptr up, mp_size_t un)
1197 sn = ((un - 1) * GMP_LIMB_BITS + mpn_limb_size_in_base_2 (up[un-1])
1200 mask = (1U << bits) - 1;
1202 for (i = 0, j = sn, shift = 0; j-- > 0;)
1204 unsigned char digit = up[i] >> shift;
1208 if (shift >= GMP_LIMB_BITS && ++i < un)
1210 shift -= GMP_LIMB_BITS;
1211 digit |= up[i] << (bits - shift);
1213 sp[j] = digit & mask;
1218 /* We generate digits from the least significant end, and reverse at
1221 mpn_limb_get_str (unsigned char *sp, mp_limb_t w,
1222 const struct gmp_div_inverse *binv)
1225 for (i = 0; w > 0; i++)
1229 h = w >> (GMP_LIMB_BITS - binv->shift);
1230 l = w << binv->shift;
1232 gmp_udiv_qrnnd_preinv (w, r, h, l, binv->d1, binv->di);
1233 assert ((r & (GMP_LIMB_MAX >> (GMP_LIMB_BITS - binv->shift))) == 0);
1242 mpn_get_str_other (unsigned char *sp,
1243 int base, const struct mpn_base_info *info,
1244 mp_ptr up, mp_size_t un)
1246 struct gmp_div_inverse binv;
1250 mpn_div_qr_1_invert (&binv, base);
1256 struct gmp_div_inverse bbinv;
1257 mpn_div_qr_1_invert (&bbinv, info->bb);
1263 w = mpn_div_qr_1_preinv (up, up, un, &bbinv);
1264 un -= (up[un-1] == 0);
1265 done = mpn_limb_get_str (sp + sn, w, &binv);
1267 for (sn += done; done < info->exp; done++)
1272 sn += mpn_limb_get_str (sp + sn, up[0], &binv);
1275 for (i = 0; 2*i + 1 < sn; i++)
1277 unsigned char t = sp[i];
1278 sp[i] = sp[sn - i - 1];
1286 mpn_get_str (unsigned char *sp, int base, mp_ptr up, mp_size_t un)
1291 assert (up[un-1] > 0);
1293 bits = mpn_base_power_of_two_p (base);
1295 return mpn_get_str_bits (sp, bits, up, un);
1298 struct mpn_base_info info;
1300 mpn_get_base_info (&info, base);
1301 return mpn_get_str_other (sp, base, &info, up, un);
1306 mpn_set_str_bits (mp_ptr rp, const unsigned char *sp, size_t sn,
1313 for (j = sn, rn = 0, shift = 0; j-- > 0; )
1322 rp[rn-1] |= (mp_limb_t) sp[j] << shift;
1324 if (shift >= GMP_LIMB_BITS)
1326 shift -= GMP_LIMB_BITS;
1328 rp[rn++] = (mp_limb_t) sp[j] >> (bits - shift);
1332 rn = mpn_normalized_size (rp, rn);
1336 /* Result is usually normalized, except for all-zero input, in which
1337 case a single zero limb is written at *RP, and 1 is returned. */
1339 mpn_set_str_other (mp_ptr rp, const unsigned char *sp, size_t sn,
1340 mp_limb_t b, const struct mpn_base_info *info)
1349 k = 1 + (sn - 1) % info->exp;
1354 w = w * b + sp[j++];
1358 for (rn = 1; j < sn;)
1363 for (k = 1; k < info->exp; k++)
1364 w = w * b + sp[j++];
1366 cy = mpn_mul_1 (rp, rp, rn, info->bb);
1367 cy += mpn_add_1 (rp, rp, rn, w);
1377 mpn_set_str (mp_ptr rp, const unsigned char *sp, size_t sn, int base)
1384 bits = mpn_base_power_of_two_p (base);
1386 return mpn_set_str_bits (rp, sp, sn, bits);
1389 struct mpn_base_info info;
1391 mpn_get_base_info (&info, base);
1392 return mpn_set_str_other (rp, sp, sn, base, &info);
1401 static const mp_limb_t dummy_limb = GMP_LIMB_MAX & 0xc1a0;
1405 r->_mp_d = (mp_ptr) &dummy_limb;
1408 /* The utility of this function is a bit limited, since many functions
1409 assigns the result variable using mpz_swap. */
1411 mpz_init2 (mpz_t r, mp_bitcnt_t bits)
1415 bits -= (bits != 0); /* Round down, except if 0 */
1416 rn = 1 + bits / GMP_LIMB_BITS;
1420 r->_mp_d = gmp_xalloc_limbs (rn);
1427 gmp_free (r->_mp_d);
1431 mpz_realloc (mpz_t r, mp_size_t size)
1433 size = GMP_MAX (size, 1);
1436 r->_mp_d = gmp_xrealloc_limbs (r->_mp_d, size);
1438 r->_mp_d = gmp_xalloc_limbs (size);
1439 r->_mp_alloc = size;
1441 if (GMP_ABS (r->_mp_size) > size)
1447 /* Realloc for an mpz_t WHAT if it has less than NEEDED limbs. */
1448 #define MPZ_REALLOC(z,n) ((n) > (z)->_mp_alloc \
1449 ? mpz_realloc(z,n) \
1452 /* MPZ assignment and basic conversions. */
1454 mpz_set_si (mpz_t r, signed long int x)
1459 if (GMP_LIMB_BITS < GMP_ULONG_BITS)
1461 mpz_set_ui (r, GMP_NEG_CAST (unsigned long int, x));
1467 MPZ_REALLOC (r, 1)[0] = GMP_NEG_CAST (unsigned long int, x);
1472 mpz_set_ui (mpz_t r, unsigned long int x)
1477 MPZ_REALLOC (r, 1)[0] = x;
1478 if (GMP_LIMB_BITS < GMP_ULONG_BITS)
1480 int LOCAL_GMP_LIMB_BITS = GMP_LIMB_BITS;
1481 while (x >>= LOCAL_GMP_LIMB_BITS)
1484 MPZ_REALLOC (r, r->_mp_size)[r->_mp_size - 1] = x;
1493 mpz_set (mpz_t r, const mpz_t x)
1495 /* Allow the NOP r == x */
1501 n = GMP_ABS (x->_mp_size);
1502 rp = MPZ_REALLOC (r, n);
1504 mpn_copyi (rp, x->_mp_d, n);
1505 r->_mp_size = x->_mp_size;
1510 mpz_init_set_si (mpz_t r, signed long int x)
1517 mpz_init_set_ui (mpz_t r, unsigned long int x)
1524 mpz_init_set (mpz_t r, const mpz_t x)
1531 mpz_fits_slong_p (const mpz_t u)
1533 return (LONG_MAX + LONG_MIN == 0 || mpz_cmp_ui (u, LONG_MAX) <= 0) &&
1534 mpz_cmpabs_ui (u, GMP_NEG_CAST (unsigned long int, LONG_MIN)) <= 0;
1538 mpn_absfits_ulong_p (mp_srcptr up, mp_size_t un)
1540 int ulongsize = GMP_ULONG_BITS / GMP_LIMB_BITS;
1541 mp_limb_t ulongrem = 0;
1543 if (GMP_ULONG_BITS % GMP_LIMB_BITS != 0)
1544 ulongrem = (mp_limb_t) (ULONG_MAX >> GMP_LIMB_BITS * ulongsize) + 1;
1546 return un <= ulongsize || (up[ulongsize] < ulongrem && un == ulongsize + 1);
1550 mpz_fits_ulong_p (const mpz_t u)
1552 mp_size_t us = u->_mp_size;
1554 return us >= 0 && mpn_absfits_ulong_p (u->_mp_d, us);
1558 mpz_get_si (const mpz_t u)
1560 unsigned long r = mpz_get_ui (u);
1561 unsigned long c = -LONG_MAX - LONG_MIN;
1563 if (u->_mp_size < 0)
1564 /* This expression is necessary to properly handle -LONG_MIN */
1565 return -(long) c - (long) ((r - c) & LONG_MAX);
1567 return (long) (r & LONG_MAX);
1571 mpz_get_ui (const mpz_t u)
1573 if (GMP_LIMB_BITS < GMP_ULONG_BITS)
1575 int LOCAL_GMP_LIMB_BITS = GMP_LIMB_BITS;
1576 unsigned long r = 0;
1577 mp_size_t n = GMP_ABS (u->_mp_size);
1578 n = GMP_MIN (n, 1 + (mp_size_t) (GMP_ULONG_BITS - 1) / GMP_LIMB_BITS);
1580 r = (r << LOCAL_GMP_LIMB_BITS) + u->_mp_d[n];
1584 return u->_mp_size == 0 ? 0 : u->_mp_d[0];
1588 mpz_size (const mpz_t u)
1590 return GMP_ABS (u->_mp_size);
1594 mpz_getlimbn (const mpz_t u, mp_size_t n)
1596 if (n >= 0 && n < GMP_ABS (u->_mp_size))
1603 mpz_realloc2 (mpz_t x, mp_bitcnt_t n)
1605 mpz_realloc (x, 1 + (n - (n != 0)) / GMP_LIMB_BITS);
1609 mpz_limbs_read (mpz_srcptr x)
1615 mpz_limbs_modify (mpz_t x, mp_size_t n)
1618 return MPZ_REALLOC (x, n);
1622 mpz_limbs_write (mpz_t x, mp_size_t n)
1624 return mpz_limbs_modify (x, n);
1628 mpz_limbs_finish (mpz_t x, mp_size_t xs)
1631 xn = mpn_normalized_size (x->_mp_d, GMP_ABS (xs));
1632 x->_mp_size = xs < 0 ? -xn : xn;
1636 mpz_roinit_normal_n (mpz_t x, mp_srcptr xp, mp_size_t xs)
1639 x->_mp_d = (mp_ptr) xp;
1645 mpz_roinit_n (mpz_t x, mp_srcptr xp, mp_size_t xs)
1647 mpz_roinit_normal_n (x, xp, xs);
1648 mpz_limbs_finish (x, xs);
1653 /* Conversions and comparison to double. */
1655 mpz_set_d (mpz_t r, double x)
1664 /* x != x is true when x is a NaN, and x == x * 0.5 is true when x is
1665 zero or infinity. */
1666 if (x != x || x == x * 0.5)
1681 B = 4.0 * (double) (GMP_LIMB_HIGHBIT >> 1);
1683 for (rn = 1; x >= B; rn++)
1686 rp = MPZ_REALLOC (r, rn);
1702 r->_mp_size = sign ? - rn : rn;
1706 mpz_init_set_d (mpz_t r, double x)
1713 mpz_get_d (const mpz_t u)
1719 double B = 4.0 * (double) (GMP_LIMB_HIGHBIT >> 1);
1721 un = GMP_ABS (u->_mp_size);
1728 m = m + GMP_DBL_MANT_BITS - GMP_LIMB_BITS;
1730 l &= GMP_LIMB_MAX << -m;
1732 for (x = l; --un >= 0;)
1739 l &= GMP_LIMB_MAX << -m;
1744 if (u->_mp_size < 0)
1751 mpz_cmpabs_d (const mpz_t x, double d)
1764 B = 4.0 * (double) (GMP_LIMB_HIGHBIT >> 1);
1767 /* Scale d so it can be compared with the top limb. */
1768 for (i = 1; i < xn; i++)
1774 /* Compare floor(d) to top limb, subtract and cancel when equal. */
1775 for (i = xn; i-- > 0;)
1792 mpz_cmp_d (const mpz_t x, double d)
1794 if (x->_mp_size < 0)
1799 return -mpz_cmpabs_d (x, d);
1806 return mpz_cmpabs_d (x, d);
1811 /* MPZ comparisons and the like. */
1813 mpz_sgn (const mpz_t u)
1815 return GMP_CMP (u->_mp_size, 0);
1819 mpz_cmp_si (const mpz_t u, long v)
1821 mp_size_t usize = u->_mp_size;
1824 return mpz_cmp_ui (u, v);
1825 else if (usize >= 0)
1828 return - mpz_cmpabs_ui (u, GMP_NEG_CAST (unsigned long int, v));
1832 mpz_cmp_ui (const mpz_t u, unsigned long v)
1834 mp_size_t usize = u->_mp_size;
1839 return mpz_cmpabs_ui (u, v);
1843 mpz_cmp (const mpz_t a, const mpz_t b)
1845 mp_size_t asize = a->_mp_size;
1846 mp_size_t bsize = b->_mp_size;
1849 return (asize < bsize) ? -1 : 1;
1850 else if (asize >= 0)
1851 return mpn_cmp (a->_mp_d, b->_mp_d, asize);
1853 return mpn_cmp (b->_mp_d, a->_mp_d, -asize);
1857 mpz_cmpabs_ui (const mpz_t u, unsigned long v)
1859 mp_size_t un = GMP_ABS (u->_mp_size);
1861 if (! mpn_absfits_ulong_p (u->_mp_d, un))
1865 unsigned long uu = mpz_get_ui (u);
1866 return GMP_CMP(uu, v);
1871 mpz_cmpabs (const mpz_t u, const mpz_t v)
1873 return mpn_cmp4 (u->_mp_d, GMP_ABS (u->_mp_size),
1874 v->_mp_d, GMP_ABS (v->_mp_size));
1878 mpz_abs (mpz_t r, const mpz_t u)
1881 r->_mp_size = GMP_ABS (r->_mp_size);
1885 mpz_neg (mpz_t r, const mpz_t u)
1888 r->_mp_size = -r->_mp_size;
1892 mpz_swap (mpz_t u, mpz_t v)
1894 MP_SIZE_T_SWAP (u->_mp_size, v->_mp_size);
1895 MP_SIZE_T_SWAP (u->_mp_alloc, v->_mp_alloc);
1896 MP_PTR_SWAP (u->_mp_d, v->_mp_d);
1900 /* MPZ addition and subtraction */
1904 mpz_add_ui (mpz_t r, const mpz_t a, unsigned long b)
1907 mpz_init_set_ui (bb, b);
1913 mpz_sub_ui (mpz_t r, const mpz_t a, unsigned long b)
1915 mpz_ui_sub (r, b, a);
1920 mpz_ui_sub (mpz_t r, unsigned long a, const mpz_t b)
1923 mpz_add_ui (r, r, a);
1927 mpz_abs_add (mpz_t r, const mpz_t a, const mpz_t b)
1929 mp_size_t an = GMP_ABS (a->_mp_size);
1930 mp_size_t bn = GMP_ABS (b->_mp_size);
1936 MPZ_SRCPTR_SWAP (a, b);
1937 MP_SIZE_T_SWAP (an, bn);
1940 rp = MPZ_REALLOC (r, an + 1);
1941 cy = mpn_add (rp, a->_mp_d, an, b->_mp_d, bn);
1949 mpz_abs_sub (mpz_t r, const mpz_t a, const mpz_t b)
1951 mp_size_t an = GMP_ABS (a->_mp_size);
1952 mp_size_t bn = GMP_ABS (b->_mp_size);
1956 cmp = mpn_cmp4 (a->_mp_d, an, b->_mp_d, bn);
1959 rp = MPZ_REALLOC (r, an);
1960 gmp_assert_nocarry (mpn_sub (rp, a->_mp_d, an, b->_mp_d, bn));
1961 return mpn_normalized_size (rp, an);
1965 rp = MPZ_REALLOC (r, bn);
1966 gmp_assert_nocarry (mpn_sub (rp, b->_mp_d, bn, a->_mp_d, an));
1967 return -mpn_normalized_size (rp, bn);
1974 mpz_add (mpz_t r, const mpz_t a, const mpz_t b)
1978 if ( (a->_mp_size ^ b->_mp_size) >= 0)
1979 rn = mpz_abs_add (r, a, b);
1981 rn = mpz_abs_sub (r, a, b);
1983 r->_mp_size = a->_mp_size >= 0 ? rn : - rn;
1987 mpz_sub (mpz_t r, const mpz_t a, const mpz_t b)
1991 if ( (a->_mp_size ^ b->_mp_size) >= 0)
1992 rn = mpz_abs_sub (r, a, b);
1994 rn = mpz_abs_add (r, a, b);
1996 r->_mp_size = a->_mp_size >= 0 ? rn : - rn;
2000 /* MPZ multiplication */
2002 mpz_mul_si (mpz_t r, const mpz_t u, long int v)
2006 mpz_mul_ui (r, u, GMP_NEG_CAST (unsigned long int, v));
2010 mpz_mul_ui (r, u, v);
2014 mpz_mul_ui (mpz_t r, const mpz_t u, unsigned long int v)
2017 mpz_init_set_ui (vv, v);
2024 mpz_mul (mpz_t r, const mpz_t u, const mpz_t v)
2027 mp_size_t un, vn, rn;
2034 if (un == 0 || vn == 0)
2040 sign = (un ^ vn) < 0;
2045 mpz_init2 (t, (un + vn) * GMP_LIMB_BITS);
2049 mpn_mul (tp, u->_mp_d, un, v->_mp_d, vn);
2051 mpn_mul (tp, v->_mp_d, vn, u->_mp_d, un);
2054 rn -= tp[rn-1] == 0;
2056 t->_mp_size = sign ? - rn : rn;
2062 mpz_mul_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bits)
2069 un = GMP_ABS (u->_mp_size);
2076 limbs = bits / GMP_LIMB_BITS;
2077 shift = bits % GMP_LIMB_BITS;
2079 rn = un + limbs + (shift > 0);
2080 rp = MPZ_REALLOC (r, rn);
2083 mp_limb_t cy = mpn_lshift (rp + limbs, u->_mp_d, un, shift);
2088 mpn_copyd (rp + limbs, u->_mp_d, un);
2090 mpn_zero (rp, limbs);
2092 r->_mp_size = (u->_mp_size < 0) ? - rn : rn;
2096 mpz_addmul_ui (mpz_t r, const mpz_t u, unsigned long int v)
2099 mpz_init_set_ui (t, v);
2106 mpz_submul_ui (mpz_t r, const mpz_t u, unsigned long int v)
2109 mpz_init_set_ui (t, v);
2116 mpz_addmul (mpz_t r, const mpz_t u, const mpz_t v)
2126 mpz_submul (mpz_t r, const mpz_t u, const mpz_t v)
2137 enum mpz_div_round_mode { GMP_DIV_FLOOR, GMP_DIV_CEIL, GMP_DIV_TRUNC };
2139 /* Allows q or r to be zero. Returns 1 iff remainder is non-zero. */
2141 mpz_div_qr (mpz_t q, mpz_t r,
2142 const mpz_t n, const mpz_t d, enum mpz_div_round_mode mode)
2144 mp_size_t ns, ds, nn, dn, qs;
2149 gmp_die("mpz_div_qr: Divide by zero.");
2167 if (mode == GMP_DIV_CEIL && qs >= 0)
2169 /* q = 1, r = n - d */
2175 else if (mode == GMP_DIV_FLOOR && qs < 0)
2177 /* q = -1, r = n + d */
2199 mpz_init_set (tr, n);
2206 mpz_init2 (tq, qn * GMP_LIMB_BITS);
2212 mpn_div_qr (qp, np, nn, d->_mp_d, dn);
2216 qn -= (qp[qn-1] == 0);
2218 tq->_mp_size = qs < 0 ? -qn : qn;
2220 rn = mpn_normalized_size (np, dn);
2221 tr->_mp_size = ns < 0 ? - rn : rn;
2223 if (mode == GMP_DIV_FLOOR && qs < 0 && rn != 0)
2226 mpz_sub_ui (tq, tq, 1);
2228 mpz_add (tr, tr, d);
2230 else if (mode == GMP_DIV_CEIL && qs >= 0 && rn != 0)
2233 mpz_add_ui (tq, tq, 1);
2235 mpz_sub (tr, tr, d);
2253 mpz_cdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d)
2255 mpz_div_qr (q, r, n, d, GMP_DIV_CEIL);
2259 mpz_fdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d)
2261 mpz_div_qr (q, r, n, d, GMP_DIV_FLOOR);
2265 mpz_tdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d)
2267 mpz_div_qr (q, r, n, d, GMP_DIV_TRUNC);
2271 mpz_cdiv_q (mpz_t q, const mpz_t n, const mpz_t d)
2273 mpz_div_qr (q, NULL, n, d, GMP_DIV_CEIL);
2277 mpz_fdiv_q (mpz_t q, const mpz_t n, const mpz_t d)
2279 mpz_div_qr (q, NULL, n, d, GMP_DIV_FLOOR);
2283 mpz_tdiv_q (mpz_t q, const mpz_t n, const mpz_t d)
2285 mpz_div_qr (q, NULL, n, d, GMP_DIV_TRUNC);
2289 mpz_cdiv_r (mpz_t r, const mpz_t n, const mpz_t d)
2291 mpz_div_qr (NULL, r, n, d, GMP_DIV_CEIL);
2295 mpz_fdiv_r (mpz_t r, const mpz_t n, const mpz_t d)
2297 mpz_div_qr (NULL, r, n, d, GMP_DIV_FLOOR);
2301 mpz_tdiv_r (mpz_t r, const mpz_t n, const mpz_t d)
2303 mpz_div_qr (NULL, r, n, d, GMP_DIV_TRUNC);
2307 mpz_mod (mpz_t r, const mpz_t n, const mpz_t d)
2309 mpz_div_qr (NULL, r, n, d, d->_mp_size >= 0 ? GMP_DIV_FLOOR : GMP_DIV_CEIL);
2313 mpz_div_q_2exp (mpz_t q, const mpz_t u, mp_bitcnt_t bit_index,
2314 enum mpz_div_round_mode mode)
2327 limb_cnt = bit_index / GMP_LIMB_BITS;
2328 qn = GMP_ABS (un) - limb_cnt;
2329 bit_index %= GMP_LIMB_BITS;
2331 if (mode == ((un > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* un != 0 here. */
2332 /* Note: Below, the final indexing at limb_cnt is valid because at
2333 that point we have qn > 0. */
2335 || !mpn_zero_p (u->_mp_d, limb_cnt)
2336 || (u->_mp_d[limb_cnt]
2337 & (((mp_limb_t) 1 << bit_index) - 1)));
2345 qp = MPZ_REALLOC (q, qn);
2349 mpn_rshift (qp, u->_mp_d + limb_cnt, qn, bit_index);
2350 qn -= qp[qn - 1] == 0;
2354 mpn_copyi (qp, u->_mp_d + limb_cnt, qn);
2361 mpz_add_ui (q, q, 1);
2367 mpz_div_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bit_index,
2368 enum mpz_div_round_mode mode)
2370 mp_size_t us, un, rn;
2375 if (us == 0 || bit_index == 0)
2380 rn = (bit_index + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS;
2383 rp = MPZ_REALLOC (r, rn);
2386 mask = GMP_LIMB_MAX >> (rn * GMP_LIMB_BITS - bit_index);
2390 /* Quotient (with truncation) is zero, and remainder is
2392 if (mode == ((us > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* us != 0 here. */
2394 /* Have to negate and sign extend. */
2397 gmp_assert_nocarry (! mpn_neg (rp, u->_mp_d, un));
2398 for (i = un; i < rn - 1; i++)
2399 rp[i] = GMP_LIMB_MAX;
2408 mpn_copyi (rp, u->_mp_d, un);
2416 mpn_copyi (rp, u->_mp_d, rn - 1);
2418 rp[rn-1] = u->_mp_d[rn-1] & mask;
2420 if (mode == ((us > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* us != 0 here. */
2422 /* If r != 0, compute 2^{bit_count} - r. */
2423 mpn_neg (rp, rp, rn);
2427 /* us is not used for anything else, so we can modify it
2428 here to indicate flipped sign. */
2432 rn = mpn_normalized_size (rp, rn);
2433 r->_mp_size = us < 0 ? -rn : rn;
2437 mpz_cdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2439 mpz_div_q_2exp (r, u, cnt, GMP_DIV_CEIL);
2443 mpz_fdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2445 mpz_div_q_2exp (r, u, cnt, GMP_DIV_FLOOR);
2449 mpz_tdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2451 mpz_div_q_2exp (r, u, cnt, GMP_DIV_TRUNC);
2455 mpz_cdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2457 mpz_div_r_2exp (r, u, cnt, GMP_DIV_CEIL);
2461 mpz_fdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2463 mpz_div_r_2exp (r, u, cnt, GMP_DIV_FLOOR);
2467 mpz_tdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2469 mpz_div_r_2exp (r, u, cnt, GMP_DIV_TRUNC);
2473 mpz_divexact (mpz_t q, const mpz_t n, const mpz_t d)
2475 gmp_assert_nocarry (mpz_div_qr (q, NULL, n, d, GMP_DIV_TRUNC));
2479 mpz_divisible_p (const mpz_t n, const mpz_t d)
2481 return mpz_div_qr (NULL, NULL, n, d, GMP_DIV_TRUNC) == 0;
2485 mpz_congruent_p (const mpz_t a, const mpz_t b, const mpz_t m)
2490 /* a == b (mod 0) iff a == b */
2491 if (mpz_sgn (m) == 0)
2492 return (mpz_cmp (a, b) == 0);
2496 res = mpz_divisible_p (t, m);
2502 static unsigned long
2503 mpz_div_qr_ui (mpz_t q, mpz_t r,
2504 const mpz_t n, unsigned long d, enum mpz_div_round_mode mode)
2510 mpz_init_set_ui (dd, d);
2511 mpz_div_qr (q, rr, n, dd, mode);
2513 ret = mpz_get_ui (rr);
2523 mpz_cdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d)
2525 return mpz_div_qr_ui (q, r, n, d, GMP_DIV_CEIL);
2529 mpz_fdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d)
2531 return mpz_div_qr_ui (q, r, n, d, GMP_DIV_FLOOR);
2535 mpz_tdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d)
2537 return mpz_div_qr_ui (q, r, n, d, GMP_DIV_TRUNC);
2541 mpz_cdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d)
2543 return mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_CEIL);
2547 mpz_fdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d)
2549 return mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_FLOOR);
2553 mpz_tdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d)
2555 return mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_TRUNC);
2559 mpz_cdiv_r_ui (mpz_t r, const mpz_t n, unsigned long d)
2561 return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_CEIL);
2564 mpz_fdiv_r_ui (mpz_t r, const mpz_t n, unsigned long d)
2566 return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_FLOOR);
2569 mpz_tdiv_r_ui (mpz_t r, const mpz_t n, unsigned long d)
2571 return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_TRUNC);
2575 mpz_cdiv_ui (const mpz_t n, unsigned long d)
2577 return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_CEIL);
2581 mpz_fdiv_ui (const mpz_t n, unsigned long d)
2583 return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_FLOOR);
2587 mpz_tdiv_ui (const mpz_t n, unsigned long d)
2589 return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_TRUNC);
2593 mpz_mod_ui (mpz_t r, const mpz_t n, unsigned long d)
2595 return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_FLOOR);
2599 mpz_divexact_ui (mpz_t q, const mpz_t n, unsigned long d)
2601 gmp_assert_nocarry (mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_TRUNC));
2605 mpz_divisible_ui_p (const mpz_t n, unsigned long d)
2607 return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_TRUNC) == 0;
2613 mpn_gcd_11 (mp_limb_t u, mp_limb_t v)
2617 assert ( (u | v) > 0);
2624 gmp_ctz (shift, u | v);
2630 MP_LIMB_T_SWAP (u, v);
2632 while ( (v & 1) == 0)
2642 while ( (u & 1) == 0);
2649 while ( (v & 1) == 0);
2656 mpz_gcd_ui (mpz_t g, const mpz_t u, unsigned long v)
2659 mpz_init_set_ui(t, v);
2673 mpz_make_odd (mpz_t r)
2677 assert (r->_mp_size > 0);
2678 /* Count trailing zeros, equivalent to mpn_scan1, because we know that there is a 1 */
2679 shift = mpn_common_scan (r->_mp_d[0], 0, r->_mp_d, 0, 0);
2680 mpz_tdiv_q_2exp (r, r, shift);
2686 mpz_gcd (mpz_t g, const mpz_t u, const mpz_t v)
2689 mp_bitcnt_t uz, vz, gz;
2691 if (u->_mp_size == 0)
2696 if (v->_mp_size == 0)
2706 uz = mpz_make_odd (tu);
2708 vz = mpz_make_odd (tv);
2709 gz = GMP_MIN (uz, vz);
2711 if (tu->_mp_size < tv->_mp_size)
2714 mpz_tdiv_r (tu, tu, tv);
2715 if (tu->_mp_size == 0)
2725 c = mpz_cmp (tu, tv);
2734 if (tv->_mp_size == 1)
2736 mp_limb_t vl = tv->_mp_d[0];
2737 mp_limb_t ul = mpz_tdiv_ui (tu, vl);
2738 mpz_set_ui (g, mpn_gcd_11 (ul, vl));
2741 mpz_sub (tu, tu, tv);
2745 mpz_mul_2exp (g, g, gz);
2749 mpz_gcdext (mpz_t g, mpz_t s, mpz_t t, const mpz_t u, const mpz_t v)
2751 mpz_t tu, tv, s0, s1, t0, t1;
2752 mp_bitcnt_t uz, vz, gz;
2755 if (u->_mp_size == 0)
2757 /* g = 0 u + sgn(v) v */
2758 signed long sign = mpz_sgn (v);
2763 mpz_set_si (t, sign);
2767 if (v->_mp_size == 0)
2769 /* g = sgn(u) u + 0 v */
2770 signed long sign = mpz_sgn (u);
2773 mpz_set_si (s, sign);
2787 uz = mpz_make_odd (tu);
2789 vz = mpz_make_odd (tv);
2790 gz = GMP_MIN (uz, vz);
2795 /* Cofactors corresponding to odd gcd. gz handled later. */
2796 if (tu->_mp_size < tv->_mp_size)
2799 MPZ_SRCPTR_SWAP (u, v);
2800 MPZ_PTR_SWAP (s, t);
2801 MP_BITCNT_T_SWAP (uz, vz);
2809 * where u and v denote the inputs with common factors of two
2810 * eliminated, and det (s0, t0; s1, t1) = 2^p. Then
2812 * 2^p tu = s1 u - t1 v
2813 * 2^p tv = -s0 u + t0 v
2816 /* After initial division, tu = q tv + tu', we have
2818 * u = 2^uz (tu' + q tv)
2823 * t0 = 2^uz, t1 = 2^uz q
2827 mpz_setbit (t0, uz);
2828 mpz_tdiv_qr (t1, tu, tu, tv);
2829 mpz_mul_2exp (t1, t1, uz);
2831 mpz_setbit (s1, vz);
2834 if (tu->_mp_size > 0)
2837 shift = mpz_make_odd (tu);
2838 mpz_mul_2exp (t0, t0, shift);
2839 mpz_mul_2exp (s0, s0, shift);
2845 c = mpz_cmp (tu, tv);
2853 * u = t0 tu + t1 (tv' + tu) = (t0 + t1) tu + t1 tv'
2854 * v = s0 tu + s1 (tv' + tu) = (s0 + s1) tu + s1 tv' */
2856 mpz_sub (tv, tv, tu);
2857 mpz_add (t0, t0, t1);
2858 mpz_add (s0, s0, s1);
2860 shift = mpz_make_odd (tv);
2861 mpz_mul_2exp (t1, t1, shift);
2862 mpz_mul_2exp (s1, s1, shift);
2866 mpz_sub (tu, tu, tv);
2867 mpz_add (t1, t0, t1);
2868 mpz_add (s1, s0, s1);
2870 shift = mpz_make_odd (tu);
2871 mpz_mul_2exp (t0, t0, shift);
2872 mpz_mul_2exp (s0, s0, shift);
2878 /* Now tv = odd part of gcd, and -s0 and t0 are corresponding
2881 mpz_mul_2exp (tv, tv, gz);
2884 /* 2^p g = s0 u + t0 v. Eliminate one factor of two at a time. To
2885 adjust cofactors, we need u / g and v / g */
2887 mpz_divexact (s1, v, tv);
2889 mpz_divexact (t1, u, tv);
2894 /* s0 u + t0 v = (s0 - v/g) u - (t0 + u/g) v */
2895 if (mpz_odd_p (s0) || mpz_odd_p (t0))
2897 mpz_sub (s0, s0, s1);
2898 mpz_add (t0, t0, t1);
2900 assert (mpz_even_p (t0) && mpz_even_p (s0));
2901 mpz_tdiv_q_2exp (s0, s0, 1);
2902 mpz_tdiv_q_2exp (t0, t0, 1);
2905 /* Arrange so that |s| < |u| / 2g */
2906 mpz_add (s1, s0, s1);
2907 if (mpz_cmpabs (s0, s1) > 0)
2910 mpz_sub (t0, t0, t1);
2912 if (u->_mp_size < 0)
2914 if (v->_mp_size < 0)
2932 mpz_lcm (mpz_t r, const mpz_t u, const mpz_t v)
2936 if (u->_mp_size == 0 || v->_mp_size == 0)
2945 mpz_divexact (g, u, g);
2953 mpz_lcm_ui (mpz_t r, const mpz_t u, unsigned long v)
2955 if (v == 0 || u->_mp_size == 0)
2961 v /= mpz_gcd_ui (NULL, u, v);
2962 mpz_mul_ui (r, u, v);
2968 mpz_invert (mpz_t r, const mpz_t u, const mpz_t m)
2973 if (u->_mp_size == 0 || mpz_cmpabs_ui (m, 1) <= 0)
2979 mpz_gcdext (g, tr, NULL, u, m);
2980 invertible = (mpz_cmp_ui (g, 1) == 0);
2984 if (tr->_mp_size < 0)
2986 if (m->_mp_size >= 0)
2987 mpz_add (tr, tr, m);
2989 mpz_sub (tr, tr, m);
3000 /* Higher level operations (sqrt, pow and root) */
3003 mpz_pow_ui (mpz_t r, const mpz_t b, unsigned long e)
3007 mpz_init_set_ui (tr, 1);
3009 bit = GMP_ULONG_HIGHBIT;
3012 mpz_mul (tr, tr, tr);
3014 mpz_mul (tr, tr, b);
3024 mpz_ui_pow_ui (mpz_t r, unsigned long blimb, unsigned long e)
3028 mpz_init_set_ui (b, blimb);
3029 mpz_pow_ui (r, b, e);
3034 mpz_powm (mpz_t r, const mpz_t b, const mpz_t e, const mpz_t m)
3040 struct gmp_div_inverse minv;
3044 en = GMP_ABS (e->_mp_size);
3045 mn = GMP_ABS (m->_mp_size);
3047 gmp_die ("mpz_powm: Zero modulo.");
3056 mpn_div_qr_invert (&minv, mp, mn);
3061 /* To avoid shifts, we do all our reductions, except the final
3062 one, using a *normalized* m. */
3065 tp = gmp_xalloc_limbs (mn);
3066 gmp_assert_nocarry (mpn_lshift (tp, mp, mn, shift));
3072 if (e->_mp_size < 0)
3074 if (!mpz_invert (base, b, m))
3075 gmp_die ("mpz_powm: Negative exponent and non-invertible base.");
3082 bn = base->_mp_size;
3085 mpn_div_qr_preinv (NULL, base->_mp_d, base->_mp_size, mp, mn, &minv);
3089 /* We have reduced the absolute value. Now take care of the
3090 sign. Note that we get zero represented non-canonically as
3092 if (b->_mp_size < 0)
3094 mp_ptr bp = MPZ_REALLOC (base, mn);
3095 gmp_assert_nocarry (mpn_sub (bp, mp, mn, bp, bn));
3098 base->_mp_size = mpn_normalized_size (base->_mp_d, bn);
3100 mpz_init_set_ui (tr, 1);
3104 mp_limb_t w = e->_mp_d[en];
3107 bit = GMP_LIMB_HIGHBIT;
3110 mpz_mul (tr, tr, tr);
3112 mpz_mul (tr, tr, base);
3113 if (tr->_mp_size > mn)
3115 mpn_div_qr_preinv (NULL, tr->_mp_d, tr->_mp_size, mp, mn, &minv);
3116 tr->_mp_size = mpn_normalized_size (tr->_mp_d, mn);
3123 /* Final reduction */
3124 if (tr->_mp_size >= mn)
3127 mpn_div_qr_preinv (NULL, tr->_mp_d, tr->_mp_size, mp, mn, &minv);
3128 tr->_mp_size = mpn_normalized_size (tr->_mp_d, mn);
3139 mpz_powm_ui (mpz_t r, const mpz_t b, unsigned long elimb, const mpz_t m)
3143 mpz_init_set_ui (e, elimb);
3144 mpz_powm (r, b, e, m);
3148 /* x=trunc(y^(1/z)), r=y-x^z */
3150 mpz_rootrem (mpz_t x, mpz_t r, const mpz_t y, unsigned long z)
3155 sgn = y->_mp_size < 0;
3156 if ((~z & sgn) != 0)
3157 gmp_die ("mpz_rootrem: Negative argument, with even root.");
3159 gmp_die ("mpz_rootrem: Zeroth root.");
3161 if (mpz_cmpabs_ui (y, 1) <= 0) {
3171 mpz_setbit (t, mpz_sizeinbase (y, 2) / z + 1);
3173 if (z == 2) /* simplify sqrt loop: z-1 == 1 */
3175 mpz_swap (u, t); /* u = x */
3176 mpz_tdiv_q (t, y, u); /* t = y/x */
3177 mpz_add (t, t, u); /* t = y/x + x */
3178 mpz_tdiv_q_2exp (t, t, 1); /* x'= (y/x + x)/2 */
3179 } while (mpz_cmpabs (t, u) < 0); /* |x'| < |x| */
3188 mpz_swap (u, t); /* u = x */
3189 mpz_pow_ui (t, u, z - 1); /* t = x^(z-1) */
3190 mpz_tdiv_q (t, y, t); /* t = y/x^(z-1) */
3191 mpz_mul_ui (v, u, z - 1); /* v = x*(z-1) */
3192 mpz_add (t, t, v); /* t = y/x^(z-1) + x*(z-1) */
3193 mpz_tdiv_q_ui (t, t, z); /* x'=(y/x^(z-1) + x*(z-1))/z */
3194 } while (mpz_cmpabs (t, u) < 0); /* |x'| < |x| */
3200 mpz_pow_ui (t, u, z);
3210 mpz_root (mpz_t x, const mpz_t y, unsigned long z)
3216 mpz_rootrem (x, r, y, z);
3217 res = r->_mp_size == 0;
3223 /* Compute s = floor(sqrt(u)) and r = u - s^2. Allows r == NULL */
3225 mpz_sqrtrem (mpz_t s, mpz_t r, const mpz_t u)
3227 mpz_rootrem (s, r, u, 2);
3231 mpz_sqrt (mpz_t s, const mpz_t u)
3233 mpz_rootrem (s, NULL, u, 2);
3237 mpz_perfect_square_p (const mpz_t u)
3239 if (u->_mp_size <= 0)
3240 return (u->_mp_size == 0);
3242 return mpz_root (NULL, u, 2);
3246 mpn_perfect_square_p (mp_srcptr p, mp_size_t n)
3251 assert (p [n-1] != 0);
3252 return mpz_root (NULL, mpz_roinit_normal_n (t, p, n), 2);
3256 mpn_sqrtrem (mp_ptr sp, mp_ptr rp, mp_srcptr p, mp_size_t n)
3262 assert (p [n-1] != 0);
3266 mpz_rootrem (s, r, mpz_roinit_normal_n (u, p, n), 2);
3268 assert (s->_mp_size == (n+1)/2);
3269 mpn_copyd (sp, s->_mp_d, s->_mp_size);
3273 mpn_copyd (rp, r->_mp_d, res);
3281 mpz_mfac_uiui (mpz_t x, unsigned long n, unsigned long m)
3283 mpz_set_ui (x, n + (n == 0));
3284 if (m + 1 < 2) return;
3286 mpz_mul_ui (x, x, n -= m);
3290 mpz_2fac_ui (mpz_t x, unsigned long n)
3292 mpz_mfac_uiui (x, n, 2);
3296 mpz_fac_ui (mpz_t x, unsigned long n)
3298 mpz_mfac_uiui (x, n, 1);
3302 mpz_bin_uiui (mpz_t r, unsigned long n, unsigned long k)
3306 mpz_set_ui (r, k <= n);
3309 k = (k <= n) ? n - k : 0;
3315 mpz_mul_ui (r, r, n--);
3317 mpz_divexact (r, r, t);
3322 /* Primality testing */
3324 /* Computes Kronecker (a/b) with odd b, a!=0 and GCD(a,b) = 1 */
3325 /* Adapted from JACOBI_BASE_METHOD==4 in mpn/generic/jacbase.c */
3327 gmp_jacobi_coprime (mp_limb_t a, mp_limb_t b)
3333 /* assert (mpn_gcd_11 (a, b) == 1); */
3335 /* Below, we represent a and b shifted right so that the least
3336 significant one bit is implicit. */
3345 /* (2/b) = -1 if b = 3 or 5 mod 8 */
3346 bit ^= c & (b ^ (b >> 1));
3364 return bit & 1 ? -1 : 1;
3368 gmp_lucas_step_k_2k (mpz_t V, mpz_t Qk, const mpz_t n)
3370 mpz_mod (Qk, Qk, n);
3371 /* V_{2k} <- V_k ^ 2 - 2Q^k */
3373 mpz_submul_ui (V, Qk, 2);
3374 mpz_tdiv_r (V, V, n);
3375 /* Q^{2k} = (Q^k)^2 */
3376 mpz_mul (Qk, Qk, Qk);
3379 /* Computes V_k, Q^k (mod n) for the Lucas' sequence */
3380 /* with P=1, Q=Q; k = (n>>b0)|1. */
3381 /* Requires an odd n > 4; b0 > 0; -2*Q must not overflow a long */
3382 /* Returns (U_k == 0) and sets V=V_k and Qk=Q^k. */
3384 gmp_lucas_mod (mpz_t V, mpz_t Qk, long Q,
3385 mp_bitcnt_t b0, const mpz_t n)
3392 assert (Q <= - (LONG_MIN / 2));
3393 assert (Q >= - (LONG_MAX / 2));
3394 assert (mpz_cmp_ui (n, 4) > 0);
3395 assert (mpz_odd_p (n));
3397 mpz_init_set_ui (U, 1); /* U1 = 1 */
3398 mpz_set_ui (V, 1); /* V1 = 1 */
3401 for (bs = mpz_sizeinbase (n, 2) - 1; --bs >= b0;)
3403 /* U_{2k} <- U_k * V_k */
3405 /* V_{2k} <- V_k ^ 2 - 2Q^k */
3406 /* Q^{2k} = (Q^k)^2 */
3407 gmp_lucas_step_k_2k (V, Qk, n);
3409 /* A step k->k+1 is performed if the bit in $n$ is 1 */
3410 /* mpz_tstbit(n,bs) or the the bit is 0 in $n$ but */
3411 /* should be 1 in $n+1$ (bs == b0) */
3412 if (b0 == bs || mpz_tstbit (n, bs))
3414 /* Q^{k+1} <- Q^k * Q */
3415 mpz_mul_si (Qk, Qk, Q);
3416 /* U_{k+1} <- (U_k + V_k) / 2 */
3417 mpz_swap (U, V); /* Keep in V the old value of U_k */
3419 /* We have to compute U/2, so we need an even value, */
3420 /* equivalent (mod n) */
3423 mpz_tdiv_q_2exp (U, U, 1);
3424 /* V_{k+1} <-(D*U_k + V_k) / 2 =
3425 U_{k+1} + (D-1)/2*U_k = U_{k+1} - 2Q*U_k */
3426 mpz_mul_si (V, V, -2*Q);
3428 mpz_tdiv_r (V, V, n);
3430 mpz_tdiv_r (U, U, n);
3433 res = U->_mp_size == 0;
3438 /* Performs strong Lucas' test on x, with parameters suggested */
3439 /* for the BPSW test. Qk is only passed to recycle a variable. */
3440 /* Requires GCD (x,6) = 1.*/
3442 gmp_stronglucas (const mpz_t x, mpz_t Qk)
3446 mp_limb_t maxD, D; /* The absolute value is stored. */
3450 /* Test on the absolute value. */
3451 mpz_roinit_normal_n (n, x->_mp_d, GMP_ABS (x->_mp_size));
3453 assert (mpz_odd_p (n));
3454 /* assert (mpz_gcd_ui (NULL, n, 6) == 1); */
3455 if (mpz_root (Qk, n, 2))
3456 return 0; /* A square is composite. */
3458 /* Check Ds up to square root (in case, n is prime)
3459 or avoid overflows */
3460 maxD = (Qk->_mp_size == 1) ? Qk->_mp_d [0] - 1 : GMP_LIMB_MAX;
3463 /* Search a D such that (D/n) = -1 in the sequence 5,-7,9,-11,.. */
3464 /* For those Ds we have (D/n) = (n/|D|) */
3468 return 1 + (D != GMP_LIMB_MAX); /* (1 + ! ~ D) */
3470 tl = mpz_tdiv_ui (n, D);
3474 while (gmp_jacobi_coprime (tl, D) == 1);
3478 /* n-(D/n) = n+1 = d*2^{b0}, with d = (n>>b0) | 1 */
3479 b0 = mpz_scan0 (n, 0);
3481 /* D= P^2 - 4Q; P = 1; Q = (1-D)/4 */
3482 Q = (D & 2) ? (long) (D >> 2) + 1 : -(long) (D >> 2);
3484 if (! gmp_lucas_mod (V, Qk, Q, b0, n)) /* If Ud != 0 */
3485 while (V->_mp_size != 0 && --b0 != 0) /* while Vk != 0 */
3486 /* V <- V ^ 2 - 2Q^k */
3487 /* Q^{2k} = (Q^k)^2 */
3488 gmp_lucas_step_k_2k (V, Qk, n);
3495 gmp_millerrabin (const mpz_t n, const mpz_t nm1, mpz_t y,
3496 const mpz_t q, mp_bitcnt_t k)
3500 /* Caller must initialize y to the base. */
3501 mpz_powm (y, y, q, n);
3503 if (mpz_cmp_ui (y, 1) == 0 || mpz_cmp (y, nm1) == 0)
3508 mpz_powm_ui (y, y, 2, n);
3509 if (mpz_cmp (y, nm1) == 0)
3511 /* y == 1 means that the previous y was a non-trivial square root
3512 of 1 (mod n). y == 0 means that n is a power of the base.
3513 In either case, n is not prime. */
3514 if (mpz_cmp_ui (y, 1) <= 0)
3520 /* This product is 0xc0cfd797, and fits in 32 bits. */
3521 #define GMP_PRIME_PRODUCT \
3522 (3UL*5UL*7UL*11UL*13UL*17UL*19UL*23UL*29UL)
3524 /* Bit (p+1)/2 is set, for each odd prime <= 61 */
3525 #define GMP_PRIME_MASK 0xc96996dcUL
3528 mpz_probab_prime_p (const mpz_t n, int reps)
3537 /* Note that we use the absolute value of n only, for compatibility
3538 with the real GMP. */
3540 return (mpz_cmpabs_ui (n, 2) == 0) ? 2 : 0;
3542 /* Above test excludes n == 0 */
3543 assert (n->_mp_size != 0);
3545 if (mpz_cmpabs_ui (n, 64) < 0)
3546 return (GMP_PRIME_MASK >> (n->_mp_d[0] >> 1)) & 2;
3548 if (mpz_gcd_ui (NULL, n, GMP_PRIME_PRODUCT) != 1)
3551 /* All prime factors are >= 31. */
3552 if (mpz_cmpabs_ui (n, 31*31) < 0)
3558 /* Find q and k, where q is odd and n = 1 + 2**k * q. */
3561 k = mpz_scan1 (nm1, 0);
3562 mpz_tdiv_q_2exp (q, nm1, k);
3565 mpz_init_set_ui (y, 2);
3566 is_prime = gmp_millerrabin (n, nm1, y, q, k) && gmp_stronglucas (n, y);
3567 reps -= 24; /* skip the first 24 repetitions */
3569 /* Use Miller-Rabin, with a deterministic sequence of bases, a[j] =
3570 j^2 + j + 41 using Euler's polynomial. We potentially stop early,
3571 if a[j] >= n - 1. Since n >= 31*31, this can happen only if reps >
3572 30 (a[30] == 971 > 31*31 == 961). */
3574 for (j = 0; is_prime & (j < reps); j++)
3576 mpz_set_ui (y, (unsigned long) j*j+j+41);
3577 if (mpz_cmp (y, nm1) >= 0)
3579 /* Don't try any further bases. This "early" break does not affect
3580 the result for any reasonable reps value (<=5000 was tested) */
3584 is_prime = gmp_millerrabin (n, nm1, y, q, k);
3594 /* Logical operations and bit manipulation. */
3596 /* Numbers are treated as if represented in two's complement (and
3597 infinitely sign extended). For a negative values we get the two's
3598 complement from -x = ~x + 1, where ~ is bitwise complement.
3607 where yyyy is the bitwise complement of xxxx. So least significant
3608 bits, up to and including the first one bit, are unchanged, and
3609 the more significant bits are all complemented.
3611 To change a bit from zero to one in a negative number, subtract the
3612 corresponding power of two from the absolute value. This can never
3613 underflow. To change a bit from one to zero, add the corresponding
3614 power of two, and this might overflow. E.g., if x = -001111, the
3615 two's complement is 110001. Clearing the least significant bit, we
3616 get two's complement 110000, and -010000. */
3619 mpz_tstbit (const mpz_t d, mp_bitcnt_t bit_index)
3621 mp_size_t limb_index;
3630 limb_index = bit_index / GMP_LIMB_BITS;
3631 if (limb_index >= dn)
3634 shift = bit_index % GMP_LIMB_BITS;
3635 w = d->_mp_d[limb_index];
3636 bit = (w >> shift) & 1;
3640 /* d < 0. Check if any of the bits below is set: If so, our bit
3641 must be complemented. */
3642 if (shift > 0 && (mp_limb_t) (w << (GMP_LIMB_BITS - shift)) > 0)
3644 while (--limb_index >= 0)
3645 if (d->_mp_d[limb_index] > 0)
3652 mpz_abs_add_bit (mpz_t d, mp_bitcnt_t bit_index)
3654 mp_size_t dn, limb_index;
3658 dn = GMP_ABS (d->_mp_size);
3660 limb_index = bit_index / GMP_LIMB_BITS;
3661 bit = (mp_limb_t) 1 << (bit_index % GMP_LIMB_BITS);
3663 if (limb_index >= dn)
3666 /* The bit should be set outside of the end of the number.
3667 We have to increase the size of the number. */
3668 dp = MPZ_REALLOC (d, limb_index + 1);
3670 dp[limb_index] = bit;
3671 for (i = dn; i < limb_index; i++)
3673 dn = limb_index + 1;
3681 cy = mpn_add_1 (dp + limb_index, dp + limb_index, dn - limb_index, bit);
3684 dp = MPZ_REALLOC (d, dn + 1);
3689 d->_mp_size = (d->_mp_size < 0) ? - dn : dn;
3693 mpz_abs_sub_bit (mpz_t d, mp_bitcnt_t bit_index)
3695 mp_size_t dn, limb_index;
3699 dn = GMP_ABS (d->_mp_size);
3702 limb_index = bit_index / GMP_LIMB_BITS;
3703 bit = (mp_limb_t) 1 << (bit_index % GMP_LIMB_BITS);
3705 assert (limb_index < dn);
3707 gmp_assert_nocarry (mpn_sub_1 (dp + limb_index, dp + limb_index,
3708 dn - limb_index, bit));
3709 dn = mpn_normalized_size (dp, dn);
3710 d->_mp_size = (d->_mp_size < 0) ? - dn : dn;
3714 mpz_setbit (mpz_t d, mp_bitcnt_t bit_index)
3716 if (!mpz_tstbit (d, bit_index))
3718 if (d->_mp_size >= 0)
3719 mpz_abs_add_bit (d, bit_index);
3721 mpz_abs_sub_bit (d, bit_index);
3726 mpz_clrbit (mpz_t d, mp_bitcnt_t bit_index)
3728 if (mpz_tstbit (d, bit_index))
3730 if (d->_mp_size >= 0)
3731 mpz_abs_sub_bit (d, bit_index);
3733 mpz_abs_add_bit (d, bit_index);
3738 mpz_combit (mpz_t d, mp_bitcnt_t bit_index)
3740 if (mpz_tstbit (d, bit_index) ^ (d->_mp_size < 0))
3741 mpz_abs_sub_bit (d, bit_index);
3743 mpz_abs_add_bit (d, bit_index);
3747 mpz_com (mpz_t r, const mpz_t u)
3749 mpz_add_ui (r, u, 1);
3754 mpz_and (mpz_t r, const mpz_t u, const mpz_t v)
3756 mp_size_t un, vn, rn, i;
3759 mp_limb_t ux, vx, rx;
3760 mp_limb_t uc, vc, rc;
3761 mp_limb_t ul, vl, rl;
3763 un = GMP_ABS (u->_mp_size);
3764 vn = GMP_ABS (v->_mp_size);
3767 MPZ_SRCPTR_SWAP (u, v);
3768 MP_SIZE_T_SWAP (un, vn);
3776 uc = u->_mp_size < 0;
3777 vc = v->_mp_size < 0;
3784 /* If the smaller input is positive, higher limbs don't matter. */
3787 rp = MPZ_REALLOC (r, rn + (mp_size_t) rc);
3795 ul = (up[i] ^ ux) + uc;
3798 vl = (vp[i] ^ vx) + vc;
3801 rl = ( (ul & vl) ^ rx) + rc;
3810 ul = (up[i] ^ ux) + uc;
3813 rl = ( (ul & vx) ^ rx) + rc;
3820 rn = mpn_normalized_size (rp, rn);
3822 r->_mp_size = rx ? -rn : rn;
3826 mpz_ior (mpz_t r, const mpz_t u, const mpz_t v)
3828 mp_size_t un, vn, rn, i;
3831 mp_limb_t ux, vx, rx;
3832 mp_limb_t uc, vc, rc;
3833 mp_limb_t ul, vl, rl;
3835 un = GMP_ABS (u->_mp_size);
3836 vn = GMP_ABS (v->_mp_size);
3839 MPZ_SRCPTR_SWAP (u, v);
3840 MP_SIZE_T_SWAP (un, vn);
3848 uc = u->_mp_size < 0;
3849 vc = v->_mp_size < 0;
3856 /* If the smaller input is negative, by sign extension higher limbs
3860 rp = MPZ_REALLOC (r, rn + (mp_size_t) rc);
3868 ul = (up[i] ^ ux) + uc;
3871 vl = (vp[i] ^ vx) + vc;
3874 rl = ( (ul | vl) ^ rx) + rc;
3883 ul = (up[i] ^ ux) + uc;
3886 rl = ( (ul | vx) ^ rx) + rc;
3893 rn = mpn_normalized_size (rp, rn);
3895 r->_mp_size = rx ? -rn : rn;
3899 mpz_xor (mpz_t r, const mpz_t u, const mpz_t v)
3901 mp_size_t un, vn, i;
3904 mp_limb_t ux, vx, rx;
3905 mp_limb_t uc, vc, rc;
3906 mp_limb_t ul, vl, rl;
3908 un = GMP_ABS (u->_mp_size);
3909 vn = GMP_ABS (v->_mp_size);
3912 MPZ_SRCPTR_SWAP (u, v);
3913 MP_SIZE_T_SWAP (un, vn);
3921 uc = u->_mp_size < 0;
3922 vc = v->_mp_size < 0;
3929 rp = MPZ_REALLOC (r, un + (mp_size_t) rc);
3937 ul = (up[i] ^ ux) + uc;
3940 vl = (vp[i] ^ vx) + vc;
3943 rl = (ul ^ vl ^ rx) + rc;
3952 ul = (up[i] ^ ux) + uc;
3955 rl = (ul ^ ux) + rc;
3962 un = mpn_normalized_size (rp, un);
3964 r->_mp_size = rx ? -un : un;
3968 gmp_popcount_limb (mp_limb_t x)
3972 /* Do 16 bits at a time, to avoid limb-sized constants. */
3973 int LOCAL_SHIFT_BITS = 16;
3976 unsigned w = x - ((x >> 1) & 0x5555);
3977 w = ((w >> 2) & 0x3333) + (w & 0x3333);
3979 w = ((w >> 8) & 0x000f) + (w & 0x000f);
3981 if (GMP_LIMB_BITS > LOCAL_SHIFT_BITS)
3982 x >>= LOCAL_SHIFT_BITS;
3990 mpn_popcount (mp_srcptr p, mp_size_t n)
3995 for (c = 0, i = 0; i < n; i++)
3996 c += gmp_popcount_limb (p[i]);
4002 mpz_popcount (const mpz_t u)
4009 return ~(mp_bitcnt_t) 0;
4011 return mpn_popcount (u->_mp_d, un);
4015 mpz_hamdist (const mpz_t u, const mpz_t v)
4017 mp_size_t un, vn, i;
4018 mp_limb_t uc, vc, ul, vl, comp;
4026 return ~(mp_bitcnt_t) 0;
4028 comp = - (uc = vc = (un < 0));
4040 MPN_SRCPTR_SWAP (up, un, vp, vn);
4042 for (i = 0, c = 0; i < vn; i++)
4044 ul = (up[i] ^ comp) + uc;
4047 vl = (vp[i] ^ comp) + vc;
4050 c += gmp_popcount_limb (ul ^ vl);
4056 ul = (up[i] ^ comp) + uc;
4059 c += gmp_popcount_limb (ul ^ comp);
4066 mpz_scan1 (const mpz_t u, mp_bitcnt_t starting_bit)
4069 mp_size_t us, un, i;
4074 i = starting_bit / GMP_LIMB_BITS;
4076 /* Past the end there's no 1 bits for u>=0, or an immediate 1 bit
4077 for u<0. Notice this test picks up any u==0 too. */
4079 return (us >= 0 ? ~(mp_bitcnt_t) 0 : starting_bit);
4085 if (starting_bit != 0)
4089 ux = mpn_zero_p (up, i);
4091 ux = - (mp_limb_t) (limb >= ux);
4094 /* Mask to 0 all bits before starting_bit, thus ignoring them. */
4095 limb &= GMP_LIMB_MAX << (starting_bit % GMP_LIMB_BITS);
4098 return mpn_common_scan (limb, i, up, un, ux);
4102 mpz_scan0 (const mpz_t u, mp_bitcnt_t starting_bit)
4105 mp_size_t us, un, i;
4109 ux = - (mp_limb_t) (us >= 0);
4111 i = starting_bit / GMP_LIMB_BITS;
4113 /* When past end, there's an immediate 0 bit for u>=0, or no 0 bits for
4114 u<0. Notice this test picks up all cases of u==0 too. */
4116 return (ux ? starting_bit : ~(mp_bitcnt_t) 0);
4122 limb -= mpn_zero_p (up, i); /* limb = ~(~limb + zero_p) */
4124 /* Mask all bits before starting_bit, thus ignoring them. */
4125 limb &= GMP_LIMB_MAX << (starting_bit % GMP_LIMB_BITS);
4127 return mpn_common_scan (limb, i, up, un, ux);
4131 /* MPZ base conversion. */
4134 mpz_sizeinbase (const mpz_t u, int base)
4140 struct gmp_div_inverse bi;
4144 assert (base <= 62);
4146 un = GMP_ABS (u->_mp_size);
4152 bits = (un - 1) * GMP_LIMB_BITS + mpn_limb_size_in_base_2 (up[un-1]);
4158 return (bits + 1) / 2;
4160 return (bits + 2) / 3;
4162 return (bits + 3) / 4;
4164 return (bits + 4) / 5;
4165 /* FIXME: Do something more clever for the common case of base
4169 tp = gmp_xalloc_limbs (un);
4170 mpn_copyi (tp, up, un);
4171 mpn_div_qr_1_invert (&bi, base);
4177 mpn_div_qr_1_preinv (tp, tp, un, &bi);
4178 un -= (tp[un-1] == 0);
4187 mpz_get_str (char *sp, int base, const mpz_t u)
4194 digits = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz";
4198 digits = "0123456789abcdefghijklmnopqrstuvwxyz";
4202 else if (base >= -1)
4211 sn = 1 + mpz_sizeinbase (u, base);
4213 sp = (char *) gmp_xalloc (1 + sn);
4215 un = GMP_ABS (u->_mp_size);
4226 if (u->_mp_size < 0)
4229 bits = mpn_base_power_of_two_p (base);
4232 /* Not modified in this case. */
4233 sn = i + mpn_get_str_bits ((unsigned char *) sp + i, bits, u->_mp_d, un);
4236 struct mpn_base_info info;
4239 mpn_get_base_info (&info, base);
4240 tp = gmp_xalloc_limbs (un);
4241 mpn_copyi (tp, u->_mp_d, un);
4243 sn = i + mpn_get_str_other ((unsigned char *) sp + i, base, &info, tp, un);
4248 sp[i] = digits[(unsigned char) sp[i]];
4255 mpz_set_str (mpz_t r, const char *sp, int base)
4257 unsigned bits, value_of_a;
4258 mp_size_t rn, alloc;
4264 assert (base == 0 || (base >= 2 && base <= 62));
4266 while (isspace( (unsigned char) *sp))
4269 sign = (*sp == '-');
4276 if (sp[1] == 'x' || sp[1] == 'X')
4281 else if (sp[1] == 'b' || sp[1] == 'B')
4298 dp = (unsigned char *) gmp_xalloc (strlen (sp));
4300 value_of_a = (base > 36) ? 36 : 10;
4301 for (dn = 0; *sp; sp++)
4305 if (isspace ((unsigned char) *sp))
4307 else if (*sp >= '0' && *sp <= '9')
4309 else if (*sp >= 'a' && *sp <= 'z')
4310 digit = *sp - 'a' + value_of_a;
4311 else if (*sp >= 'A' && *sp <= 'Z')
4312 digit = *sp - 'A' + 10;
4314 digit = base; /* fail */
4316 if (digit >= (unsigned) base)
4332 bits = mpn_base_power_of_two_p (base);
4336 alloc = (dn * bits + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS;
4337 rp = MPZ_REALLOC (r, alloc);
4338 rn = mpn_set_str_bits (rp, dp, dn, bits);
4342 struct mpn_base_info info;
4343 mpn_get_base_info (&info, base);
4344 alloc = (dn + info.exp - 1) / info.exp;
4345 rp = MPZ_REALLOC (r, alloc);
4346 rn = mpn_set_str_other (rp, dp, dn, base, &info);
4347 /* Normalization, needed for all-zero input. */
4349 rn -= rp[rn-1] == 0;
4351 assert (rn <= alloc);
4354 r->_mp_size = sign ? - rn : rn;
4360 mpz_init_set_str (mpz_t r, const char *sp, int base)
4363 return mpz_set_str (r, sp, base);
4367 mpz_out_str (FILE *stream, int base, const mpz_t x)
4372 str = mpz_get_str (NULL, base, x);
4374 len = fwrite (str, 1, len, stream);
4381 gmp_detect_endian (void)
4383 static const int i = 2;
4384 const unsigned char *p = (const unsigned char *) &i;
4388 /* Import and export. Does not support nails. */
4390 mpz_import (mpz_t r, size_t count, int order, size_t size, int endian,
4391 size_t nails, const void *src)
4393 const unsigned char *p;
4394 ptrdiff_t word_step;
4398 /* The current (partial) limb. */
4400 /* The number of bytes already copied to this limb (starting from
4403 /* The index where the limb should be stored, when completed. */
4407 gmp_die ("mpz_import: Nails not supported.");
4409 assert (order == 1 || order == -1);
4410 assert (endian >= -1 && endian <= 1);
4413 endian = gmp_detect_endian ();
4415 p = (unsigned char *) src;
4417 word_step = (order != endian) ? 2 * size : 0;
4419 /* Process bytes from the least significant end, so point p at the
4420 least significant word. */
4423 p += size * (count - 1);
4424 word_step = - word_step;
4427 /* And at least significant byte of that word. */
4431 rn = (size * count + sizeof(mp_limb_t) - 1) / sizeof(mp_limb_t);
4432 rp = MPZ_REALLOC (r, rn);
4434 for (limb = 0, bytes = 0, i = 0; count > 0; count--, p += word_step)
4437 for (j = 0; j < size; j++, p -= (ptrdiff_t) endian)
4439 limb |= (mp_limb_t) *p << (bytes++ * CHAR_BIT);
4440 if (bytes == sizeof(mp_limb_t))
4448 assert (i + (bytes > 0) == rn);
4452 i = mpn_normalized_size (rp, i);
4458 mpz_export (void *r, size_t *countp, int order, size_t size, int endian,
4459 size_t nails, const mpz_t u)
4465 gmp_die ("mpz_import: Nails not supported.");
4467 assert (order == 1 || order == -1);
4468 assert (endian >= -1 && endian <= 1);
4469 assert (size > 0 || u->_mp_size == 0);
4477 ptrdiff_t word_step;
4478 /* The current (partial) limb. */
4480 /* The number of bytes left to to in this limb. */
4482 /* The index where the limb was read. */
4487 /* Count bytes in top limb. */
4488 limb = u->_mp_d[un-1];
4491 k = (GMP_LIMB_BITS <= CHAR_BIT);
4495 int LOCAL_CHAR_BIT = CHAR_BIT;
4496 k++; limb >>= LOCAL_CHAR_BIT;
4497 } while (limb != 0);
4499 /* else limb = 0; */
4501 count = (k + (un-1) * sizeof (mp_limb_t) + size - 1) / size;
4504 r = gmp_xalloc (count * size);
4507 endian = gmp_detect_endian ();
4509 p = (unsigned char *) r;
4511 word_step = (order != endian) ? 2 * size : 0;
4513 /* Process bytes from the least significant end, so point p at the
4514 least significant word. */
4517 p += size * (count - 1);
4518 word_step = - word_step;
4521 /* And at least significant byte of that word. */
4525 for (bytes = 0, i = 0, k = 0; k < count; k++, p += word_step)
4528 for (j = 0; j < size; ++j, p -= (ptrdiff_t) endian)
4530 if (sizeof (mp_limb_t) == 1)
4539 int LOCAL_CHAR_BIT = CHAR_BIT;
4543 limb = u->_mp_d[i++];
4544 bytes = sizeof (mp_limb_t);
4547 limb >>= LOCAL_CHAR_BIT;
4553 assert (k == count);