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-2015 Free Software Foundation, Inc.
7 This file is part of the GNU MP Library.
9 The GNU MP Library is free software; you can redistribute it and/or modify
10 it under the terms of either:
12 * the GNU Lesser General Public License as published by the Free
13 Software Foundation; either version 3 of the License, or (at your
14 option) any later version.
18 * the GNU General Public License as published by the Free Software
19 Foundation; either version 2 of the License, or (at your option) any
22 or both in parallel, as here.
24 The GNU MP Library is distributed in the hope that it will be useful, but
25 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
26 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
29 You should have received copies of the GNU General Public License and the
30 GNU Lesser General Public License along with the GNU MP Library. If not,
31 see https://www.gnu.org/licenses/. */
33 /* NOTE: All functions in this file which are not declared in
34 mini-gmp.h are internal, and are not intended to be compatible
35 neither with GMP nor with future versions of mini-gmp. */
37 /* Much of the material copied from GMP files, including: gmp-impl.h,
38 longlong.h, mpn/generic/add_n.c, mpn/generic/addmul_1.c,
39 mpn/generic/lshift.c, mpn/generic/mul_1.c,
40 mpn/generic/mul_basecase.c, mpn/generic/rshift.c,
41 mpn/generic/sbpi1_div_qr.c, mpn/generic/sub_n.c,
42 mpn/generic/submul_1.c. */
55 #define GMP_LIMB_BITS (sizeof(mp_limb_t) * CHAR_BIT)
57 #define GMP_LIMB_MAX (~ (mp_limb_t) 0)
58 #define GMP_LIMB_HIGHBIT ((mp_limb_t) 1 << (GMP_LIMB_BITS - 1))
60 #define GMP_HLIMB_BIT ((mp_limb_t) 1 << (GMP_LIMB_BITS / 2))
61 #define GMP_LLIMB_MASK (GMP_HLIMB_BIT - 1)
63 #define GMP_ULONG_BITS (sizeof(unsigned long) * CHAR_BIT)
64 #define GMP_ULONG_HIGHBIT ((unsigned long) 1 << (GMP_ULONG_BITS - 1))
66 #define GMP_ABS(x) ((x) >= 0 ? (x) : -(x))
67 #define GMP_NEG_CAST(T,x) (-((T)((x) + 1) - 1))
69 #define GMP_MIN(a, b) ((a) < (b) ? (a) : (b))
70 #define GMP_MAX(a, b) ((a) > (b) ? (a) : (b))
72 #define gmp_assert_nocarry(x) do { \
77 #define gmp_clz(count, x) do { \
78 mp_limb_t __clz_x = (x); \
81 (__clz_x & ((mp_limb_t) 0xff << (GMP_LIMB_BITS - 8))) == 0; \
84 for (; (__clz_x & GMP_LIMB_HIGHBIT) == 0; __clz_c++) \
89 #define gmp_ctz(count, x) do { \
90 mp_limb_t __ctz_x = (x); \
91 unsigned __ctz_c = 0; \
92 gmp_clz (__ctz_c, __ctz_x & - __ctz_x); \
93 (count) = GMP_LIMB_BITS - 1 - __ctz_c; \
96 #define gmp_add_ssaaaa(sh, sl, ah, al, bh, bl) \
100 (sh) = (ah) + (bh) + (__x < (al)); \
104 #define gmp_sub_ddmmss(sh, sl, ah, al, bh, bl) \
108 (sh) = (ah) - (bh) - ((al) < (bl)); \
112 #define gmp_umul_ppmm(w1, w0, u, v) \
114 mp_limb_t __x0, __x1, __x2, __x3; \
115 unsigned __ul, __vl, __uh, __vh; \
116 mp_limb_t __u = (u), __v = (v); \
118 __ul = __u & GMP_LLIMB_MASK; \
119 __uh = __u >> (GMP_LIMB_BITS / 2); \
120 __vl = __v & GMP_LLIMB_MASK; \
121 __vh = __v >> (GMP_LIMB_BITS / 2); \
123 __x0 = (mp_limb_t) __ul * __vl; \
124 __x1 = (mp_limb_t) __ul * __vh; \
125 __x2 = (mp_limb_t) __uh * __vl; \
126 __x3 = (mp_limb_t) __uh * __vh; \
128 __x1 += __x0 >> (GMP_LIMB_BITS / 2);/* this can't give carry */ \
129 __x1 += __x2; /* but this indeed can */ \
130 if (__x1 < __x2) /* did we get it? */ \
131 __x3 += GMP_HLIMB_BIT; /* yes, add it in the proper pos. */ \
133 (w1) = __x3 + (__x1 >> (GMP_LIMB_BITS / 2)); \
134 (w0) = (__x1 << (GMP_LIMB_BITS / 2)) + (__x0 & GMP_LLIMB_MASK); \
137 #define gmp_udiv_qrnnd_preinv(q, r, nh, nl, d, di) \
139 mp_limb_t _qh, _ql, _r, _mask; \
140 gmp_umul_ppmm (_qh, _ql, (nh), (di)); \
141 gmp_add_ssaaaa (_qh, _ql, _qh, _ql, (nh) + 1, (nl)); \
142 _r = (nl) - _qh * (d); \
143 _mask = -(mp_limb_t) (_r > _ql); /* both > and >= are OK */ \
156 #define gmp_udiv_qr_3by2(q, r1, r0, n2, n1, n0, d1, d0, dinv) \
158 mp_limb_t _q0, _t1, _t0, _mask; \
159 gmp_umul_ppmm ((q), _q0, (n2), (dinv)); \
160 gmp_add_ssaaaa ((q), _q0, (q), _q0, (n2), (n1)); \
162 /* Compute the two most significant limbs of n - q'd */ \
163 (r1) = (n1) - (d1) * (q); \
164 gmp_sub_ddmmss ((r1), (r0), (r1), (n0), (d1), (d0)); \
165 gmp_umul_ppmm (_t1, _t0, (d0), (q)); \
166 gmp_sub_ddmmss ((r1), (r0), (r1), (r0), _t1, _t0); \
169 /* Conditionally adjust q and the remainders */ \
170 _mask = - (mp_limb_t) ((r1) >= _q0); \
172 gmp_add_ssaaaa ((r1), (r0), (r1), (r0), _mask & (d1), _mask & (d0)); \
175 if ((r1) > (d1) || (r0) >= (d0)) \
178 gmp_sub_ddmmss ((r1), (r0), (r1), (r0), (d1), (d0)); \
184 #define MP_LIMB_T_SWAP(x, y) \
186 mp_limb_t __mp_limb_t_swap__tmp = (x); \
188 (y) = __mp_limb_t_swap__tmp; \
190 #define MP_SIZE_T_SWAP(x, y) \
192 mp_size_t __mp_size_t_swap__tmp = (x); \
194 (y) = __mp_size_t_swap__tmp; \
196 #define MP_BITCNT_T_SWAP(x,y) \
198 mp_bitcnt_t __mp_bitcnt_t_swap__tmp = (x); \
200 (y) = __mp_bitcnt_t_swap__tmp; \
202 #define MP_PTR_SWAP(x, y) \
204 mp_ptr __mp_ptr_swap__tmp = (x); \
206 (y) = __mp_ptr_swap__tmp; \
208 #define MP_SRCPTR_SWAP(x, y) \
210 mp_srcptr __mp_srcptr_swap__tmp = (x); \
212 (y) = __mp_srcptr_swap__tmp; \
215 #define MPN_PTR_SWAP(xp,xs, yp,ys) \
217 MP_PTR_SWAP (xp, yp); \
218 MP_SIZE_T_SWAP (xs, ys); \
220 #define MPN_SRCPTR_SWAP(xp,xs, yp,ys) \
222 MP_SRCPTR_SWAP (xp, yp); \
223 MP_SIZE_T_SWAP (xs, ys); \
226 #define MPZ_PTR_SWAP(x, y) \
228 mpz_ptr __mpz_ptr_swap__tmp = (x); \
230 (y) = __mpz_ptr_swap__tmp; \
232 #define MPZ_SRCPTR_SWAP(x, y) \
234 mpz_srcptr __mpz_srcptr_swap__tmp = (x); \
236 (y) = __mpz_srcptr_swap__tmp; \
239 const int mp_bits_per_limb
= GMP_LIMB_BITS
;
242 /* Memory allocation and other helper functions. */
244 gmp_die (const char *msg
)
246 fprintf (stderr
, "%s\n", msg
);
251 gmp_default_alloc (size_t size
)
259 gmp_die("gmp_default_alloc: Virtual memory exhausted.");
265 gmp_default_realloc (void *old
, size_t old_size
, size_t new_size
)
269 p
= realloc (old
, new_size
);
272 gmp_die("gmp_default_realloc: Virtual memory exhausted.");
278 gmp_default_free (void *p
, size_t size
)
283 static void * (*gmp_allocate_func
) (size_t) = gmp_default_alloc
;
284 static void * (*gmp_reallocate_func
) (void *, size_t, size_t) = gmp_default_realloc
;
285 static void (*gmp_free_func
) (void *, size_t) = gmp_default_free
;
288 mp_get_memory_functions (void *(**alloc_func
) (size_t),
289 void *(**realloc_func
) (void *, size_t, size_t),
290 void (**free_func
) (void *, size_t))
293 *alloc_func
= gmp_allocate_func
;
296 *realloc_func
= gmp_reallocate_func
;
299 *free_func
= gmp_free_func
;
303 mp_set_memory_functions (void *(*alloc_func
) (size_t),
304 void *(*realloc_func
) (void *, size_t, size_t),
305 void (*free_func
) (void *, size_t))
308 alloc_func
= gmp_default_alloc
;
310 realloc_func
= gmp_default_realloc
;
312 free_func
= gmp_default_free
;
314 gmp_allocate_func
= alloc_func
;
315 gmp_reallocate_func
= realloc_func
;
316 gmp_free_func
= free_func
;
319 #define gmp_xalloc(size) ((*gmp_allocate_func)((size)))
320 #define gmp_free(p) ((*gmp_free_func) ((p), 0))
323 gmp_xalloc_limbs (mp_size_t size
)
325 return (mp_ptr
) gmp_xalloc (size
* sizeof (mp_limb_t
));
329 gmp_xrealloc_limbs (mp_ptr old
, mp_size_t size
)
332 return (mp_ptr
) (*gmp_reallocate_func
) (old
, 0, size
* sizeof (mp_limb_t
));
339 mpn_copyi (mp_ptr d
, mp_srcptr s
, mp_size_t n
)
342 for (i
= 0; i
< n
; i
++)
347 mpn_copyd (mp_ptr d
, mp_srcptr s
, mp_size_t n
)
354 mpn_cmp (mp_srcptr ap
, mp_srcptr bp
, mp_size_t n
)
359 return ap
[n
] > bp
[n
] ? 1 : -1;
365 mpn_cmp4 (mp_srcptr ap
, mp_size_t an
, mp_srcptr bp
, mp_size_t bn
)
368 return an
< bn
? -1 : 1;
370 return mpn_cmp (ap
, bp
, an
);
374 mpn_normalized_size (mp_srcptr xp
, mp_size_t n
)
376 while (n
> 0 && xp
[n
-1] == 0)
382 mpn_zero_p(mp_srcptr rp
, mp_size_t n
)
384 return mpn_normalized_size (rp
, n
) == 0;
388 mpn_zero (mp_ptr rp
, mp_size_t n
)
395 mpn_add_1 (mp_ptr rp
, mp_srcptr ap
, mp_size_t n
, mp_limb_t b
)
403 mp_limb_t r
= ap
[i
] + b
;
414 mpn_add_n (mp_ptr rp
, mp_srcptr ap
, mp_srcptr bp
, mp_size_t n
)
419 for (i
= 0, cy
= 0; i
< n
; i
++)
422 a
= ap
[i
]; b
= bp
[i
];
433 mpn_add (mp_ptr rp
, mp_srcptr ap
, mp_size_t an
, mp_srcptr bp
, mp_size_t bn
)
439 cy
= mpn_add_n (rp
, ap
, bp
, bn
);
441 cy
= mpn_add_1 (rp
+ bn
, ap
+ bn
, an
- bn
, cy
);
446 mpn_sub_1 (mp_ptr rp
, mp_srcptr ap
, mp_size_t n
, mp_limb_t b
)
457 mp_limb_t cy
= a
< b
;;
467 mpn_sub_n (mp_ptr rp
, mp_srcptr ap
, mp_srcptr bp
, mp_size_t n
)
472 for (i
= 0, cy
= 0; i
< n
; i
++)
475 a
= ap
[i
]; b
= bp
[i
];
485 mpn_sub (mp_ptr rp
, mp_srcptr ap
, mp_size_t an
, mp_srcptr bp
, mp_size_t bn
)
491 cy
= mpn_sub_n (rp
, ap
, bp
, bn
);
493 cy
= mpn_sub_1 (rp
+ bn
, ap
+ bn
, an
- bn
, cy
);
498 mpn_mul_1 (mp_ptr rp
, mp_srcptr up
, mp_size_t n
, mp_limb_t vl
)
500 mp_limb_t ul
, cl
, hpl
, lpl
;
508 gmp_umul_ppmm (hpl
, lpl
, ul
, vl
);
511 cl
= (lpl
< cl
) + hpl
;
521 mpn_addmul_1 (mp_ptr rp
, mp_srcptr up
, mp_size_t n
, mp_limb_t vl
)
523 mp_limb_t ul
, cl
, hpl
, lpl
, rl
;
531 gmp_umul_ppmm (hpl
, lpl
, ul
, vl
);
534 cl
= (lpl
< cl
) + hpl
;
547 mpn_submul_1 (mp_ptr rp
, mp_srcptr up
, mp_size_t n
, mp_limb_t vl
)
549 mp_limb_t ul
, cl
, hpl
, lpl
, rl
;
557 gmp_umul_ppmm (hpl
, lpl
, ul
, vl
);
560 cl
= (lpl
< cl
) + hpl
;
573 mpn_mul (mp_ptr rp
, mp_srcptr up
, mp_size_t un
, mp_srcptr vp
, mp_size_t vn
)
578 /* We first multiply by the low order limb. This result can be
579 stored, not added, to rp. We also avoid a loop for zeroing this
582 rp
[un
] = mpn_mul_1 (rp
, up
, un
, vp
[0]);
584 /* Now accumulate the product of up[] and the next higher limb from
590 rp
[un
] = mpn_addmul_1 (rp
, up
, un
, vp
[0]);
596 mpn_mul_n (mp_ptr rp
, mp_srcptr ap
, mp_srcptr bp
, mp_size_t n
)
598 mpn_mul (rp
, ap
, n
, bp
, n
);
602 mpn_sqr (mp_ptr rp
, mp_srcptr ap
, mp_size_t n
)
604 mpn_mul (rp
, ap
, n
, ap
, n
);
608 mpn_lshift (mp_ptr rp
, mp_srcptr up
, mp_size_t n
, unsigned int cnt
)
610 mp_limb_t high_limb
, low_limb
;
616 assert (cnt
< GMP_LIMB_BITS
);
621 tnc
= GMP_LIMB_BITS
- cnt
;
623 retval
= low_limb
>> tnc
;
624 high_limb
= (low_limb
<< cnt
);
629 *--rp
= high_limb
| (low_limb
>> tnc
);
630 high_limb
= (low_limb
<< cnt
);
638 mpn_rshift (mp_ptr rp
, mp_srcptr up
, mp_size_t n
, unsigned int cnt
)
640 mp_limb_t high_limb
, low_limb
;
646 assert (cnt
< GMP_LIMB_BITS
);
648 tnc
= GMP_LIMB_BITS
- cnt
;
650 retval
= (high_limb
<< tnc
);
651 low_limb
= high_limb
>> cnt
;
656 *rp
++ = low_limb
| (high_limb
<< tnc
);
657 low_limb
= high_limb
>> cnt
;
665 mpn_common_scan (mp_limb_t limb
, mp_size_t i
, mp_srcptr up
, mp_size_t un
,
670 assert (ux
== 0 || ux
== GMP_LIMB_MAX
);
671 assert (0 <= i
&& i
<= un
);
677 return (ux
== 0 ? ~(mp_bitcnt_t
) 0 : un
* GMP_LIMB_BITS
);
681 return (mp_bitcnt_t
) i
* GMP_LIMB_BITS
+ cnt
;
685 mpn_scan1 (mp_srcptr ptr
, mp_bitcnt_t bit
)
688 i
= bit
/ GMP_LIMB_BITS
;
690 return mpn_common_scan ( ptr
[i
] & (GMP_LIMB_MAX
<< (bit
% GMP_LIMB_BITS
)),
695 mpn_scan0 (mp_srcptr ptr
, mp_bitcnt_t bit
)
698 i
= bit
/ GMP_LIMB_BITS
;
700 return mpn_common_scan (~ptr
[i
] & (GMP_LIMB_MAX
<< (bit
% GMP_LIMB_BITS
)),
701 i
, ptr
, i
, GMP_LIMB_MAX
);
705 /* MPN division interface. */
707 mpn_invert_3by2 (mp_limb_t u1
, mp_limb_t u0
)
713 /* First, do a 2/1 inverse. */
714 /* The inverse m is defined as floor( (B^2 - 1 - u1)/u1 ), so that 0 <
715 * B^2 - (B + m) u1 <= u1 */
716 assert (u1
>= GMP_LIMB_HIGHBIT
);
718 ul
= u1
& GMP_LLIMB_MASK
;
719 uh
= u1
>> (GMP_LIMB_BITS
/ 2);
722 r
= ((~u1
- (mp_limb_t
) qh
* uh
) << (GMP_LIMB_BITS
/ 2)) | GMP_LLIMB_MASK
;
724 p
= (mp_limb_t
) qh
* ul
;
725 /* Adjustment steps taken from udiv_qrnnd_c */
730 if (r
>= u1
) /* i.e. we didn't get carry when adding to r */
739 /* Do a 3/2 division (with half limb size) */
740 p
= (r
>> (GMP_LIMB_BITS
/ 2)) * qh
+ r
;
741 ql
= (p
>> (GMP_LIMB_BITS
/ 2)) + 1;
743 /* By the 3/2 method, we don't need the high half limb. */
744 r
= (r
<< (GMP_LIMB_BITS
/ 2)) + GMP_LLIMB_MASK
- ql
* u1
;
746 if (r
>= (p
<< (GMP_LIMB_BITS
/ 2)))
751 m
= ((mp_limb_t
) qh
<< (GMP_LIMB_BITS
/ 2)) + ql
;
773 gmp_umul_ppmm (th
, tl
, u0
, m
);
778 m
-= ((r
> u1
) | ((r
== u1
) & (tl
> u0
)));
785 struct gmp_div_inverse
787 /* Normalization shift count. */
789 /* Normalized divisor (d0 unused for mpn_div_qr_1) */
791 /* Inverse, for 2/1 or 3/2. */
796 mpn_div_qr_1_invert (struct gmp_div_inverse
*inv
, mp_limb_t d
)
803 inv
->d1
= d
<< shift
;
804 inv
->di
= mpn_invert_limb (inv
->d1
);
808 mpn_div_qr_2_invert (struct gmp_div_inverse
*inv
,
809 mp_limb_t d1
, mp_limb_t d0
)
818 d1
= (d1
<< shift
) | (d0
>> (GMP_LIMB_BITS
- shift
));
823 inv
->di
= mpn_invert_3by2 (d1
, d0
);
827 mpn_div_qr_invert (struct gmp_div_inverse
*inv
,
828 mp_srcptr dp
, mp_size_t dn
)
833 mpn_div_qr_1_invert (inv
, dp
[0]);
835 mpn_div_qr_2_invert (inv
, dp
[1], dp
[0]);
848 d1
= (d1
<< shift
) | (d0
>> (GMP_LIMB_BITS
- shift
));
849 d0
= (d0
<< shift
) | (dp
[dn
-3] >> (GMP_LIMB_BITS
- shift
));
853 inv
->di
= mpn_invert_3by2 (d1
, d0
);
857 /* Not matching current public gmp interface, rather corresponding to
858 the sbpi1_div_* functions. */
860 mpn_div_qr_1_preinv (mp_ptr qp
, mp_srcptr np
, mp_size_t nn
,
861 const struct gmp_div_inverse
*inv
)
869 tp
= gmp_xalloc_limbs (nn
);
870 r
= mpn_lshift (tp
, np
, nn
, inv
->shift
);
882 gmp_udiv_qrnnd_preinv (q
, r
, r
, np
[nn
], d
, di
);
889 return r
>> inv
->shift
;
893 mpn_div_qr_1 (mp_ptr qp
, mp_srcptr np
, mp_size_t nn
, mp_limb_t d
)
897 /* Special case for powers of two. */
898 if ((d
& (d
-1)) == 0)
900 mp_limb_t r
= np
[0] & (d
-1);
904 mpn_copyi (qp
, np
, nn
);
909 mpn_rshift (qp
, np
, nn
, shift
);
916 struct gmp_div_inverse inv
;
917 mpn_div_qr_1_invert (&inv
, d
);
918 return mpn_div_qr_1_preinv (qp
, np
, nn
, &inv
);
923 mpn_div_qr_2_preinv (mp_ptr qp
, mp_ptr rp
, mp_srcptr np
, mp_size_t nn
,
924 const struct gmp_div_inverse
*inv
)
928 mp_limb_t d1
, d0
, di
, r1
, r0
;
939 tp
= gmp_xalloc_limbs (nn
);
940 r1
= mpn_lshift (tp
, np
, nn
, shift
);
953 gmp_udiv_qr_3by2 (q
, r1
, r0
, r1
, r0
, n0
, d1
, d0
, di
);
962 assert ((r0
<< (GMP_LIMB_BITS
- shift
)) == 0);
963 r0
= (r0
>> shift
) | (r1
<< (GMP_LIMB_BITS
- shift
));
975 mpn_div_qr_2 (mp_ptr qp
, mp_ptr rp
, mp_srcptr np
, mp_size_t nn
,
976 mp_limb_t d1
, mp_limb_t d0
)
978 struct gmp_div_inverse inv
;
981 mpn_div_qr_2_invert (&inv
, d1
, d0
);
982 mpn_div_qr_2_preinv (qp
, rp
, np
, nn
, &inv
);
987 mpn_div_qr_pi1 (mp_ptr qp
,
988 mp_ptr np
, mp_size_t nn
, mp_limb_t n1
,
989 mp_srcptr dp
, mp_size_t dn
,
1004 assert ((d1
& GMP_LIMB_HIGHBIT
) != 0);
1005 /* Iteration variable is the index of the q limb.
1007 * We divide <n1, np[dn-1+i], np[dn-2+i], np[dn-3+i],..., np[i]>
1008 * by <d1, d0, dp[dn-3], ..., dp[0] >
1014 mp_limb_t n0
= np
[dn
-1+i
];
1016 if (n1
== d1
&& n0
== d0
)
1019 mpn_submul_1 (np
+i
, dp
, dn
, q
);
1020 n1
= np
[dn
-1+i
]; /* update n1, last loop's value will now be invalid */
1024 gmp_udiv_qr_3by2 (q
, n1
, n0
, n1
, n0
, np
[dn
-2+i
], d1
, d0
, dinv
);
1026 cy
= mpn_submul_1 (np
+ i
, dp
, dn
-2, q
);
1036 n1
+= d1
+ mpn_add_n (np
+ i
, np
+ i
, dp
, dn
- 1);
1050 mpn_div_qr_preinv (mp_ptr qp
, mp_ptr np
, mp_size_t nn
,
1051 mp_srcptr dp
, mp_size_t dn
,
1052 const struct gmp_div_inverse
*inv
)
1058 np
[0] = mpn_div_qr_1_preinv (qp
, np
, nn
, inv
);
1060 mpn_div_qr_2_preinv (qp
, np
, np
, nn
, inv
);
1066 assert (inv
->d1
== dp
[dn
-1]);
1067 assert (inv
->d0
== dp
[dn
-2]);
1068 assert ((inv
->d1
& GMP_LIMB_HIGHBIT
) != 0);
1072 nh
= mpn_lshift (np
, np
, nn
, shift
);
1076 mpn_div_qr_pi1 (qp
, np
, nn
, nh
, dp
, dn
, inv
->di
);
1079 gmp_assert_nocarry (mpn_rshift (np
, np
, dn
, shift
));
1084 mpn_div_qr (mp_ptr qp
, mp_ptr np
, mp_size_t nn
, mp_srcptr dp
, mp_size_t dn
)
1086 struct gmp_div_inverse inv
;
1092 mpn_div_qr_invert (&inv
, dp
, dn
);
1093 if (dn
> 2 && inv
.shift
> 0)
1095 tp
= gmp_xalloc_limbs (dn
);
1096 gmp_assert_nocarry (mpn_lshift (tp
, dp
, dn
, inv
.shift
));
1099 mpn_div_qr_preinv (qp
, np
, nn
, dp
, dn
, &inv
);
1105 /* MPN base conversion. */
1107 mpn_base_power_of_two_p (unsigned b
)
1123 struct mpn_base_info
1125 /* bb is the largest power of the base which fits in one limb, and
1126 exp is the corresponding exponent. */
1132 mpn_get_base_info (struct mpn_base_info
*info
, mp_limb_t b
)
1138 m
= GMP_LIMB_MAX
/ b
;
1139 for (exp
= 1, p
= b
; p
<= m
; exp
++)
1147 mpn_limb_size_in_base_2 (mp_limb_t u
)
1153 return GMP_LIMB_BITS
- shift
;
1157 mpn_get_str_bits (unsigned char *sp
, unsigned bits
, mp_srcptr up
, mp_size_t un
)
1164 sn
= ((un
- 1) * GMP_LIMB_BITS
+ mpn_limb_size_in_base_2 (up
[un
-1])
1167 mask
= (1U << bits
) - 1;
1169 for (i
= 0, j
= sn
, shift
= 0; j
-- > 0;)
1171 unsigned char digit
= up
[i
] >> shift
;
1175 if (shift
>= GMP_LIMB_BITS
&& ++i
< un
)
1177 shift
-= GMP_LIMB_BITS
;
1178 digit
|= up
[i
] << (bits
- shift
);
1180 sp
[j
] = digit
& mask
;
1185 /* We generate digits from the least significant end, and reverse at
1188 mpn_limb_get_str (unsigned char *sp
, mp_limb_t w
,
1189 const struct gmp_div_inverse
*binv
)
1192 for (i
= 0; w
> 0; i
++)
1196 h
= w
>> (GMP_LIMB_BITS
- binv
->shift
);
1197 l
= w
<< binv
->shift
;
1199 gmp_udiv_qrnnd_preinv (w
, r
, h
, l
, binv
->d1
, binv
->di
);
1200 assert ( (r
<< (GMP_LIMB_BITS
- binv
->shift
)) == 0);
1209 mpn_get_str_other (unsigned char *sp
,
1210 int base
, const struct mpn_base_info
*info
,
1211 mp_ptr up
, mp_size_t un
)
1213 struct gmp_div_inverse binv
;
1217 mpn_div_qr_1_invert (&binv
, base
);
1223 struct gmp_div_inverse bbinv
;
1224 mpn_div_qr_1_invert (&bbinv
, info
->bb
);
1230 w
= mpn_div_qr_1_preinv (up
, up
, un
, &bbinv
);
1231 un
-= (up
[un
-1] == 0);
1232 done
= mpn_limb_get_str (sp
+ sn
, w
, &binv
);
1234 for (sn
+= done
; done
< info
->exp
; done
++)
1239 sn
+= mpn_limb_get_str (sp
+ sn
, up
[0], &binv
);
1242 for (i
= 0; 2*i
+ 1 < sn
; i
++)
1244 unsigned char t
= sp
[i
];
1245 sp
[i
] = sp
[sn
- i
- 1];
1253 mpn_get_str (unsigned char *sp
, int base
, mp_ptr up
, mp_size_t un
)
1258 assert (up
[un
-1] > 0);
1260 bits
= mpn_base_power_of_two_p (base
);
1262 return mpn_get_str_bits (sp
, bits
, up
, un
);
1265 struct mpn_base_info info
;
1267 mpn_get_base_info (&info
, base
);
1268 return mpn_get_str_other (sp
, base
, &info
, up
, un
);
1273 mpn_set_str_bits (mp_ptr rp
, const unsigned char *sp
, size_t sn
,
1280 for (j
= sn
, rn
= 0, shift
= 0; j
-- > 0; )
1289 rp
[rn
-1] |= (mp_limb_t
) sp
[j
] << shift
;
1291 if (shift
>= GMP_LIMB_BITS
)
1293 shift
-= GMP_LIMB_BITS
;
1295 rp
[rn
++] = (mp_limb_t
) sp
[j
] >> (bits
- shift
);
1299 rn
= mpn_normalized_size (rp
, rn
);
1304 mpn_set_str_other (mp_ptr rp
, const unsigned char *sp
, size_t sn
,
1305 mp_limb_t b
, const struct mpn_base_info
*info
)
1312 k
= 1 + (sn
- 1) % info
->exp
;
1317 w
= w
* b
+ sp
[j
++];
1321 for (rn
= (w
> 0); j
< sn
;)
1326 for (k
= 1; k
< info
->exp
; k
++)
1327 w
= w
* b
+ sp
[j
++];
1329 cy
= mpn_mul_1 (rp
, rp
, rn
, info
->bb
);
1330 cy
+= mpn_add_1 (rp
, rp
, rn
, w
);
1340 mpn_set_str (mp_ptr rp
, const unsigned char *sp
, size_t sn
, int base
)
1347 bits
= mpn_base_power_of_two_p (base
);
1349 return mpn_set_str_bits (rp
, sp
, sn
, bits
);
1352 struct mpn_base_info info
;
1354 mpn_get_base_info (&info
, base
);
1355 return mpn_set_str_other (rp
, sp
, sn
, base
, &info
);
1366 r
->_mp_d
= gmp_xalloc_limbs (1);
1369 /* The utility of this function is a bit limited, since many functions
1370 assigns the result variable using mpz_swap. */
1372 mpz_init2 (mpz_t r
, mp_bitcnt_t bits
)
1376 bits
-= (bits
!= 0); /* Round down, except if 0 */
1377 rn
= 1 + bits
/ GMP_LIMB_BITS
;
1381 r
->_mp_d
= gmp_xalloc_limbs (rn
);
1387 gmp_free (r
->_mp_d
);
1391 mpz_realloc (mpz_t r
, mp_size_t size
)
1393 size
= GMP_MAX (size
, 1);
1395 r
->_mp_d
= gmp_xrealloc_limbs (r
->_mp_d
, size
);
1396 r
->_mp_alloc
= size
;
1398 if (GMP_ABS (r
->_mp_size
) > size
)
1404 /* Realloc for an mpz_t WHAT if it has less than NEEDED limbs. */
1405 #define MPZ_REALLOC(z,n) ((n) > (z)->_mp_alloc \
1406 ? mpz_realloc(z,n) \
1409 /* MPZ assignment and basic conversions. */
1411 mpz_set_si (mpz_t r
, signed long int x
)
1418 r
->_mp_d
[0] = GMP_NEG_CAST (unsigned long int, x
);
1423 mpz_set_ui (mpz_t r
, unsigned long int x
)
1435 mpz_set (mpz_t r
, const mpz_t x
)
1437 /* Allow the NOP r == x */
1443 n
= GMP_ABS (x
->_mp_size
);
1444 rp
= MPZ_REALLOC (r
, n
);
1446 mpn_copyi (rp
, x
->_mp_d
, n
);
1447 r
->_mp_size
= x
->_mp_size
;
1452 mpz_init_set_si (mpz_t r
, signed long int x
)
1459 mpz_init_set_ui (mpz_t r
, unsigned long int x
)
1466 mpz_init_set (mpz_t r
, const mpz_t x
)
1473 mpz_fits_slong_p (const mpz_t u
)
1475 mp_size_t us
= u
->_mp_size
;
1478 return u
->_mp_d
[0] < GMP_LIMB_HIGHBIT
;
1480 return u
->_mp_d
[0] <= GMP_LIMB_HIGHBIT
;
1486 mpz_fits_ulong_p (const mpz_t u
)
1488 mp_size_t us
= u
->_mp_size
;
1490 return (us
== (us
> 0));
1494 mpz_get_si (const mpz_t u
)
1496 mp_size_t us
= u
->_mp_size
;
1499 return (long) (u
->_mp_d
[0] & ~GMP_LIMB_HIGHBIT
);
1501 return (long) (- u
->_mp_d
[0] | GMP_LIMB_HIGHBIT
);
1507 mpz_get_ui (const mpz_t u
)
1509 return u
->_mp_size
== 0 ? 0 : u
->_mp_d
[0];
1513 mpz_size (const mpz_t u
)
1515 return GMP_ABS (u
->_mp_size
);
1519 mpz_getlimbn (const mpz_t u
, mp_size_t n
)
1521 if (n
>= 0 && n
< GMP_ABS (u
->_mp_size
))
1528 mpz_realloc2 (mpz_t x
, mp_bitcnt_t n
)
1530 mpz_realloc (x
, 1 + (n
- (n
!= 0)) / GMP_LIMB_BITS
);
1534 mpz_limbs_read (mpz_srcptr x
)
1540 mpz_limbs_modify (mpz_t x
, mp_size_t n
)
1543 return MPZ_REALLOC (x
, n
);
1547 mpz_limbs_write (mpz_t x
, mp_size_t n
)
1549 return mpz_limbs_modify (x
, n
);
1553 mpz_limbs_finish (mpz_t x
, mp_size_t xs
)
1556 xn
= mpn_normalized_size (x
->_mp_d
, GMP_ABS (xs
));
1557 x
->_mp_size
= xs
< 0 ? -xn
: xn
;
1561 mpz_roinit_n (mpz_t x
, mp_srcptr xp
, mp_size_t xs
)
1564 x
->_mp_d
= (mp_ptr
) xp
;
1565 mpz_limbs_finish (x
, xs
);
1570 /* Conversions and comparison to double. */
1572 mpz_set_d (mpz_t r
, double x
)
1581 /* x != x is true when x is a NaN, and x == x * 0.5 is true when x is
1582 zero or infinity. */
1583 if (x
!= x
|| x
== x
* 0.5)
1598 B
= 2.0 * (double) GMP_LIMB_HIGHBIT
;
1600 for (rn
= 1; x
>= B
; rn
++)
1603 rp
= MPZ_REALLOC (r
, rn
);
1619 r
->_mp_size
= sign
? - rn
: rn
;
1623 mpz_init_set_d (mpz_t r
, double x
)
1630 mpz_get_d (const mpz_t u
)
1634 double B
= 2.0 * (double) GMP_LIMB_HIGHBIT
;
1636 un
= GMP_ABS (u
->_mp_size
);
1643 x
= B
*x
+ u
->_mp_d
[--un
];
1645 if (u
->_mp_size
< 0)
1652 mpz_cmpabs_d (const mpz_t x
, double d
)
1665 B
= 2.0 * (double) GMP_LIMB_HIGHBIT
;
1668 /* Scale d so it can be compared with the top limb. */
1669 for (i
= 1; i
< xn
; i
++)
1675 /* Compare floor(d) to top limb, subtract and cancel when equal. */
1676 for (i
= xn
; i
-- > 0;)
1693 mpz_cmp_d (const mpz_t x
, double d
)
1695 if (x
->_mp_size
< 0)
1700 return -mpz_cmpabs_d (x
, d
);
1707 return mpz_cmpabs_d (x
, d
);
1712 /* MPZ comparisons and the like. */
1714 mpz_sgn (const mpz_t u
)
1716 mp_size_t usize
= u
->_mp_size
;
1718 return (usize
> 0) - (usize
< 0);
1722 mpz_cmp_si (const mpz_t u
, long v
)
1724 mp_size_t usize
= u
->_mp_size
;
1729 return mpz_cmp_ui (u
, v
);
1730 else if (usize
>= 0)
1732 else /* usize == -1 */
1734 mp_limb_t ul
= u
->_mp_d
[0];
1735 if ((mp_limb_t
)GMP_NEG_CAST (unsigned long int, v
) < ul
)
1738 return (mp_limb_t
)GMP_NEG_CAST (unsigned long int, v
) > ul
;
1743 mpz_cmp_ui (const mpz_t u
, unsigned long v
)
1745 mp_size_t usize
= u
->_mp_size
;
1753 mp_limb_t ul
= (usize
> 0) ? u
->_mp_d
[0] : 0;
1754 return (ul
> v
) - (ul
< v
);
1759 mpz_cmp (const mpz_t a
, const mpz_t b
)
1761 mp_size_t asize
= a
->_mp_size
;
1762 mp_size_t bsize
= b
->_mp_size
;
1765 return (asize
< bsize
) ? -1 : 1;
1766 else if (asize
>= 0)
1767 return mpn_cmp (a
->_mp_d
, b
->_mp_d
, asize
);
1769 return mpn_cmp (b
->_mp_d
, a
->_mp_d
, -asize
);
1773 mpz_cmpabs_ui (const mpz_t u
, unsigned long v
)
1775 mp_size_t un
= GMP_ABS (u
->_mp_size
);
1781 ul
= (un
== 1) ? u
->_mp_d
[0] : 0;
1783 return (ul
> v
) - (ul
< v
);
1787 mpz_cmpabs (const mpz_t u
, const mpz_t v
)
1789 return mpn_cmp4 (u
->_mp_d
, GMP_ABS (u
->_mp_size
),
1790 v
->_mp_d
, GMP_ABS (v
->_mp_size
));
1794 mpz_abs (mpz_t r
, const mpz_t u
)
1797 r
->_mp_size
= GMP_ABS (r
->_mp_size
);
1801 mpz_neg (mpz_t r
, const mpz_t u
)
1804 r
->_mp_size
= -r
->_mp_size
;
1808 mpz_swap (mpz_t u
, mpz_t v
)
1810 MP_SIZE_T_SWAP (u
->_mp_size
, v
->_mp_size
);
1811 MP_SIZE_T_SWAP (u
->_mp_alloc
, v
->_mp_alloc
);
1812 MP_PTR_SWAP (u
->_mp_d
, v
->_mp_d
);
1816 /* MPZ addition and subtraction */
1818 /* Adds to the absolute value. Returns new size, but doesn't store it. */
1820 mpz_abs_add_ui (mpz_t r
, const mpz_t a
, unsigned long b
)
1826 an
= GMP_ABS (a
->_mp_size
);
1833 rp
= MPZ_REALLOC (r
, an
+ 1);
1835 cy
= mpn_add_1 (rp
, a
->_mp_d
, an
, b
);
1842 /* Subtract from the absolute value. Returns new size, (or -1 on underflow),
1843 but doesn't store it. */
1845 mpz_abs_sub_ui (mpz_t r
, const mpz_t a
, unsigned long b
)
1847 mp_size_t an
= GMP_ABS (a
->_mp_size
);
1848 mp_ptr rp
= MPZ_REALLOC (r
, an
);
1855 else if (an
== 1 && a
->_mp_d
[0] < b
)
1857 rp
[0] = b
- a
->_mp_d
[0];
1862 gmp_assert_nocarry (mpn_sub_1 (rp
, a
->_mp_d
, an
, b
));
1863 return mpn_normalized_size (rp
, an
);
1868 mpz_add_ui (mpz_t r
, const mpz_t a
, unsigned long b
)
1870 if (a
->_mp_size
>= 0)
1871 r
->_mp_size
= mpz_abs_add_ui (r
, a
, b
);
1873 r
->_mp_size
= -mpz_abs_sub_ui (r
, a
, b
);
1877 mpz_sub_ui (mpz_t r
, const mpz_t a
, unsigned long b
)
1879 if (a
->_mp_size
< 0)
1880 r
->_mp_size
= -mpz_abs_add_ui (r
, a
, b
);
1882 r
->_mp_size
= mpz_abs_sub_ui (r
, a
, b
);
1886 mpz_ui_sub (mpz_t r
, unsigned long a
, const mpz_t b
)
1888 if (b
->_mp_size
< 0)
1889 r
->_mp_size
= mpz_abs_add_ui (r
, b
, a
);
1891 r
->_mp_size
= -mpz_abs_sub_ui (r
, b
, a
);
1895 mpz_abs_add (mpz_t r
, const mpz_t a
, const mpz_t b
)
1897 mp_size_t an
= GMP_ABS (a
->_mp_size
);
1898 mp_size_t bn
= GMP_ABS (b
->_mp_size
);
1904 MPZ_SRCPTR_SWAP (a
, b
);
1905 MP_SIZE_T_SWAP (an
, bn
);
1908 rp
= MPZ_REALLOC (r
, an
+ 1);
1909 cy
= mpn_add (rp
, a
->_mp_d
, an
, b
->_mp_d
, bn
);
1917 mpz_abs_sub (mpz_t r
, const mpz_t a
, const mpz_t b
)
1919 mp_size_t an
= GMP_ABS (a
->_mp_size
);
1920 mp_size_t bn
= GMP_ABS (b
->_mp_size
);
1924 cmp
= mpn_cmp4 (a
->_mp_d
, an
, b
->_mp_d
, bn
);
1927 rp
= MPZ_REALLOC (r
, an
);
1928 gmp_assert_nocarry (mpn_sub (rp
, a
->_mp_d
, an
, b
->_mp_d
, bn
));
1929 return mpn_normalized_size (rp
, an
);
1933 rp
= MPZ_REALLOC (r
, bn
);
1934 gmp_assert_nocarry (mpn_sub (rp
, b
->_mp_d
, bn
, a
->_mp_d
, an
));
1935 return -mpn_normalized_size (rp
, bn
);
1942 mpz_add (mpz_t r
, const mpz_t a
, const mpz_t b
)
1946 if ( (a
->_mp_size
^ b
->_mp_size
) >= 0)
1947 rn
= mpz_abs_add (r
, a
, b
);
1949 rn
= mpz_abs_sub (r
, a
, b
);
1951 r
->_mp_size
= a
->_mp_size
>= 0 ? rn
: - rn
;
1955 mpz_sub (mpz_t r
, const mpz_t a
, const mpz_t b
)
1959 if ( (a
->_mp_size
^ b
->_mp_size
) >= 0)
1960 rn
= mpz_abs_sub (r
, a
, b
);
1962 rn
= mpz_abs_add (r
, a
, b
);
1964 r
->_mp_size
= a
->_mp_size
>= 0 ? rn
: - rn
;
1968 /* MPZ multiplication */
1970 mpz_mul_si (mpz_t r
, const mpz_t u
, long int v
)
1974 mpz_mul_ui (r
, u
, GMP_NEG_CAST (unsigned long int, v
));
1978 mpz_mul_ui (r
, u
, (unsigned long int) v
);
1982 mpz_mul_ui (mpz_t r
, const mpz_t u
, unsigned long int v
)
1990 if (us
== 0 || v
== 0)
1998 tp
= MPZ_REALLOC (r
, un
+ 1);
1999 cy
= mpn_mul_1 (tp
, u
->_mp_d
, un
, v
);
2003 r
->_mp_size
= (us
< 0) ? - un
: un
;
2007 mpz_mul (mpz_t r
, const mpz_t u
, const mpz_t v
)
2010 mp_size_t un
, vn
, rn
;
2017 if (un
== 0 || vn
== 0)
2023 sign
= (un
^ vn
) < 0;
2028 mpz_init2 (t
, (un
+ vn
) * GMP_LIMB_BITS
);
2032 mpn_mul (tp
, u
->_mp_d
, un
, v
->_mp_d
, vn
);
2034 mpn_mul (tp
, v
->_mp_d
, vn
, u
->_mp_d
, un
);
2037 rn
-= tp
[rn
-1] == 0;
2039 t
->_mp_size
= sign
? - rn
: rn
;
2045 mpz_mul_2exp (mpz_t r
, const mpz_t u
, mp_bitcnt_t bits
)
2052 un
= GMP_ABS (u
->_mp_size
);
2059 limbs
= bits
/ GMP_LIMB_BITS
;
2060 shift
= bits
% GMP_LIMB_BITS
;
2062 rn
= un
+ limbs
+ (shift
> 0);
2063 rp
= MPZ_REALLOC (r
, rn
);
2066 mp_limb_t cy
= mpn_lshift (rp
+ limbs
, u
->_mp_d
, un
, shift
);
2071 mpn_copyd (rp
+ limbs
, u
->_mp_d
, un
);
2073 mpn_zero (rp
, limbs
);
2075 r
->_mp_size
= (u
->_mp_size
< 0) ? - rn
: rn
;
2079 mpz_addmul_ui (mpz_t r
, const mpz_t u
, unsigned long int v
)
2083 mpz_mul_ui (t
, u
, v
);
2089 mpz_submul_ui (mpz_t r
, const mpz_t u
, unsigned long int v
)
2093 mpz_mul_ui (t
, u
, v
);
2099 mpz_addmul (mpz_t r
, const mpz_t u
, const mpz_t v
)
2109 mpz_submul (mpz_t r
, const mpz_t u
, const mpz_t v
)
2120 enum mpz_div_round_mode
{ GMP_DIV_FLOOR
, GMP_DIV_CEIL
, GMP_DIV_TRUNC
};
2122 /* Allows q or r to be zero. Returns 1 iff remainder is non-zero. */
2124 mpz_div_qr (mpz_t q
, mpz_t r
,
2125 const mpz_t n
, const mpz_t d
, enum mpz_div_round_mode mode
)
2127 mp_size_t ns
, ds
, nn
, dn
, qs
;
2132 gmp_die("mpz_div_qr: Divide by zero.");
2150 if (mode
== GMP_DIV_CEIL
&& qs
>= 0)
2152 /* q = 1, r = n - d */
2158 else if (mode
== GMP_DIV_FLOOR
&& qs
< 0)
2160 /* q = -1, r = n + d */
2182 mpz_init_set (tr
, n
);
2189 mpz_init2 (tq
, qn
* GMP_LIMB_BITS
);
2195 mpn_div_qr (qp
, np
, nn
, d
->_mp_d
, dn
);
2199 qn
-= (qp
[qn
-1] == 0);
2201 tq
->_mp_size
= qs
< 0 ? -qn
: qn
;
2203 rn
= mpn_normalized_size (np
, dn
);
2204 tr
->_mp_size
= ns
< 0 ? - rn
: rn
;
2206 if (mode
== GMP_DIV_FLOOR
&& qs
< 0 && rn
!= 0)
2209 mpz_sub_ui (tq
, tq
, 1);
2211 mpz_add (tr
, tr
, d
);
2213 else if (mode
== GMP_DIV_CEIL
&& qs
>= 0 && rn
!= 0)
2216 mpz_add_ui (tq
, tq
, 1);
2218 mpz_sub (tr
, tr
, d
);
2236 mpz_cdiv_qr (mpz_t q
, mpz_t r
, const mpz_t n
, const mpz_t d
)
2238 mpz_div_qr (q
, r
, n
, d
, GMP_DIV_CEIL
);
2242 mpz_fdiv_qr (mpz_t q
, mpz_t r
, const mpz_t n
, const mpz_t d
)
2244 mpz_div_qr (q
, r
, n
, d
, GMP_DIV_FLOOR
);
2248 mpz_tdiv_qr (mpz_t q
, mpz_t r
, const mpz_t n
, const mpz_t d
)
2250 mpz_div_qr (q
, r
, n
, d
, GMP_DIV_TRUNC
);
2254 mpz_cdiv_q (mpz_t q
, const mpz_t n
, const mpz_t d
)
2256 mpz_div_qr (q
, NULL
, n
, d
, GMP_DIV_CEIL
);
2260 mpz_fdiv_q (mpz_t q
, const mpz_t n
, const mpz_t d
)
2262 mpz_div_qr (q
, NULL
, n
, d
, GMP_DIV_FLOOR
);
2266 mpz_tdiv_q (mpz_t q
, const mpz_t n
, const mpz_t d
)
2268 mpz_div_qr (q
, NULL
, n
, d
, GMP_DIV_TRUNC
);
2272 mpz_cdiv_r (mpz_t r
, const mpz_t n
, const mpz_t d
)
2274 mpz_div_qr (NULL
, r
, n
, d
, GMP_DIV_CEIL
);
2278 mpz_fdiv_r (mpz_t r
, const mpz_t n
, const mpz_t d
)
2280 mpz_div_qr (NULL
, r
, n
, d
, GMP_DIV_FLOOR
);
2284 mpz_tdiv_r (mpz_t r
, const mpz_t n
, const mpz_t d
)
2286 mpz_div_qr (NULL
, r
, n
, d
, GMP_DIV_TRUNC
);
2290 mpz_mod (mpz_t r
, const mpz_t n
, const mpz_t d
)
2292 mpz_div_qr (NULL
, r
, n
, d
, d
->_mp_size
>= 0 ? GMP_DIV_FLOOR
: GMP_DIV_CEIL
);
2296 mpz_div_q_2exp (mpz_t q
, const mpz_t u
, mp_bitcnt_t bit_index
,
2297 enum mpz_div_round_mode mode
)
2310 limb_cnt
= bit_index
/ GMP_LIMB_BITS
;
2311 qn
= GMP_ABS (un
) - limb_cnt
;
2312 bit_index
%= GMP_LIMB_BITS
;
2314 if (mode
== ((un
> 0) ? GMP_DIV_CEIL
: GMP_DIV_FLOOR
)) /* un != 0 here. */
2315 /* Note: Below, the final indexing at limb_cnt is valid because at
2316 that point we have qn > 0. */
2318 || !mpn_zero_p (u
->_mp_d
, limb_cnt
)
2319 || (u
->_mp_d
[limb_cnt
]
2320 & (((mp_limb_t
) 1 << bit_index
) - 1)));
2329 qp
= MPZ_REALLOC (q
, qn
);
2333 mpn_rshift (qp
, u
->_mp_d
+ limb_cnt
, qn
, bit_index
);
2334 qn
-= qp
[qn
- 1] == 0;
2338 mpn_copyi (qp
, u
->_mp_d
+ limb_cnt
, qn
);
2345 mpz_add_ui (q
, q
, 1);
2351 mpz_div_r_2exp (mpz_t r
, const mpz_t u
, mp_bitcnt_t bit_index
,
2352 enum mpz_div_round_mode mode
)
2354 mp_size_t us
, un
, rn
;
2359 if (us
== 0 || bit_index
== 0)
2364 rn
= (bit_index
+ GMP_LIMB_BITS
- 1) / GMP_LIMB_BITS
;
2367 rp
= MPZ_REALLOC (r
, rn
);
2370 mask
= GMP_LIMB_MAX
>> (rn
* GMP_LIMB_BITS
- bit_index
);
2374 /* Quotient (with truncation) is zero, and remainder is
2376 if (mode
== ((us
> 0) ? GMP_DIV_CEIL
: GMP_DIV_FLOOR
)) /* us != 0 here. */
2378 /* Have to negate and sign extend. */
2382 for (cy
= 1, i
= 0; i
< un
; i
++)
2384 mp_limb_t s
= ~u
->_mp_d
[i
] + cy
;
2389 for (; i
< rn
- 1; i
++)
2390 rp
[i
] = GMP_LIMB_MAX
;
2399 mpn_copyi (rp
, u
->_mp_d
, un
);
2407 mpn_copyi (rp
, u
->_mp_d
, rn
- 1);
2409 rp
[rn
-1] = u
->_mp_d
[rn
-1] & mask
;
2411 if (mode
== ((us
> 0) ? GMP_DIV_CEIL
: GMP_DIV_FLOOR
)) /* us != 0 here. */
2413 /* If r != 0, compute 2^{bit_count} - r. */
2416 for (i
= 0; i
< rn
&& rp
[i
] == 0; i
++)
2420 /* r > 0, need to flip sign. */
2427 /* us is not used for anything else, so we can modify it
2428 here to indicate flipped sign. */
2433 rn
= mpn_normalized_size (rp
, rn
);
2434 r
->_mp_size
= us
< 0 ? -rn
: rn
;
2438 mpz_cdiv_q_2exp (mpz_t r
, const mpz_t u
, mp_bitcnt_t cnt
)
2440 mpz_div_q_2exp (r
, u
, cnt
, GMP_DIV_CEIL
);
2444 mpz_fdiv_q_2exp (mpz_t r
, const mpz_t u
, mp_bitcnt_t cnt
)
2446 mpz_div_q_2exp (r
, u
, cnt
, GMP_DIV_FLOOR
);
2450 mpz_tdiv_q_2exp (mpz_t r
, const mpz_t u
, mp_bitcnt_t cnt
)
2452 mpz_div_q_2exp (r
, u
, cnt
, GMP_DIV_TRUNC
);
2456 mpz_cdiv_r_2exp (mpz_t r
, const mpz_t u
, mp_bitcnt_t cnt
)
2458 mpz_div_r_2exp (r
, u
, cnt
, GMP_DIV_CEIL
);
2462 mpz_fdiv_r_2exp (mpz_t r
, const mpz_t u
, mp_bitcnt_t cnt
)
2464 mpz_div_r_2exp (r
, u
, cnt
, GMP_DIV_FLOOR
);
2468 mpz_tdiv_r_2exp (mpz_t r
, const mpz_t u
, mp_bitcnt_t cnt
)
2470 mpz_div_r_2exp (r
, u
, cnt
, GMP_DIV_TRUNC
);
2474 mpz_divexact (mpz_t q
, const mpz_t n
, const mpz_t d
)
2476 gmp_assert_nocarry (mpz_div_qr (q
, NULL
, n
, d
, GMP_DIV_TRUNC
));
2480 mpz_divisible_p (const mpz_t n
, const mpz_t d
)
2482 return mpz_div_qr (NULL
, NULL
, n
, d
, GMP_DIV_TRUNC
) == 0;
2486 mpz_congruent_p (const mpz_t a
, const mpz_t b
, const mpz_t m
)
2491 /* a == b (mod 0) iff a == b */
2492 if (mpz_sgn (m
) == 0)
2493 return (mpz_cmp (a
, b
) == 0);
2497 res
= mpz_divisible_p (t
, m
);
2503 static unsigned long
2504 mpz_div_qr_ui (mpz_t q
, mpz_t r
,
2505 const mpz_t n
, unsigned long d
, enum mpz_div_round_mode mode
)
2524 qp
= MPZ_REALLOC (q
, qn
);
2528 rl
= mpn_div_qr_1 (qp
, n
->_mp_d
, qn
, d
);
2532 rs
= (ns
< 0) ? -rs
: rs
;
2534 if (rl
> 0 && ( (mode
== GMP_DIV_FLOOR
&& ns
< 0)
2535 || (mode
== GMP_DIV_CEIL
&& ns
>= 0)))
2538 gmp_assert_nocarry (mpn_add_1 (qp
, qp
, qn
, 1));
2550 qn
-= (qp
[qn
-1] == 0);
2551 assert (qn
== 0 || qp
[qn
-1] > 0);
2553 q
->_mp_size
= (ns
< 0) ? - qn
: qn
;
2560 mpz_cdiv_qr_ui (mpz_t q
, mpz_t r
, const mpz_t n
, unsigned long d
)
2562 return mpz_div_qr_ui (q
, r
, n
, d
, GMP_DIV_CEIL
);
2566 mpz_fdiv_qr_ui (mpz_t q
, mpz_t r
, const mpz_t n
, unsigned long d
)
2568 return mpz_div_qr_ui (q
, r
, n
, d
, GMP_DIV_FLOOR
);
2572 mpz_tdiv_qr_ui (mpz_t q
, mpz_t r
, const mpz_t n
, unsigned long d
)
2574 return mpz_div_qr_ui (q
, r
, n
, d
, GMP_DIV_TRUNC
);
2578 mpz_cdiv_q_ui (mpz_t q
, const mpz_t n
, unsigned long d
)
2580 return mpz_div_qr_ui (q
, NULL
, n
, d
, GMP_DIV_CEIL
);
2584 mpz_fdiv_q_ui (mpz_t q
, const mpz_t n
, unsigned long d
)
2586 return mpz_div_qr_ui (q
, NULL
, n
, d
, GMP_DIV_FLOOR
);
2590 mpz_tdiv_q_ui (mpz_t q
, const mpz_t n
, unsigned long d
)
2592 return mpz_div_qr_ui (q
, NULL
, n
, d
, GMP_DIV_TRUNC
);
2596 mpz_cdiv_r_ui (mpz_t r
, const mpz_t n
, unsigned long d
)
2598 return mpz_div_qr_ui (NULL
, r
, n
, d
, GMP_DIV_CEIL
);
2601 mpz_fdiv_r_ui (mpz_t r
, const mpz_t n
, unsigned long d
)
2603 return mpz_div_qr_ui (NULL
, r
, n
, d
, GMP_DIV_FLOOR
);
2606 mpz_tdiv_r_ui (mpz_t r
, const mpz_t n
, unsigned long d
)
2608 return mpz_div_qr_ui (NULL
, r
, n
, d
, GMP_DIV_TRUNC
);
2612 mpz_cdiv_ui (const mpz_t n
, unsigned long d
)
2614 return mpz_div_qr_ui (NULL
, NULL
, n
, d
, GMP_DIV_CEIL
);
2618 mpz_fdiv_ui (const mpz_t n
, unsigned long d
)
2620 return mpz_div_qr_ui (NULL
, NULL
, n
, d
, GMP_DIV_FLOOR
);
2624 mpz_tdiv_ui (const mpz_t n
, unsigned long d
)
2626 return mpz_div_qr_ui (NULL
, NULL
, n
, d
, GMP_DIV_TRUNC
);
2630 mpz_mod_ui (mpz_t r
, const mpz_t n
, unsigned long d
)
2632 return mpz_div_qr_ui (NULL
, r
, n
, d
, GMP_DIV_FLOOR
);
2636 mpz_divexact_ui (mpz_t q
, const mpz_t n
, unsigned long d
)
2638 gmp_assert_nocarry (mpz_div_qr_ui (q
, NULL
, n
, d
, GMP_DIV_TRUNC
));
2642 mpz_divisible_ui_p (const mpz_t n
, unsigned long d
)
2644 return mpz_div_qr_ui (NULL
, NULL
, n
, d
, GMP_DIV_TRUNC
) == 0;
2650 mpn_gcd_11 (mp_limb_t u
, mp_limb_t v
)
2654 assert ( (u
| v
) > 0);
2661 gmp_ctz (shift
, u
| v
);
2667 MP_LIMB_T_SWAP (u
, v
);
2669 while ( (v
& 1) == 0)
2679 while ( (u
& 1) == 0);
2686 while ( (v
& 1) == 0);
2693 mpz_gcd_ui (mpz_t g
, const mpz_t u
, unsigned long v
)
2704 un
= GMP_ABS (u
->_mp_size
);
2706 v
= mpn_gcd_11 (mpn_div_qr_1 (NULL
, u
->_mp_d
, un
, v
), v
);
2716 mpz_make_odd (mpz_t r
)
2720 assert (r
->_mp_size
> 0);
2721 /* Count trailing zeros, equivalent to mpn_scan1, because we know that there is a 1 */
2722 shift
= mpn_common_scan (r
->_mp_d
[0], 0, r
->_mp_d
, 0, 0);
2723 mpz_tdiv_q_2exp (r
, r
, shift
);
2729 mpz_gcd (mpz_t g
, const mpz_t u
, const mpz_t v
)
2732 mp_bitcnt_t uz
, vz
, gz
;
2734 if (u
->_mp_size
== 0)
2739 if (v
->_mp_size
== 0)
2749 uz
= mpz_make_odd (tu
);
2751 vz
= mpz_make_odd (tv
);
2752 gz
= GMP_MIN (uz
, vz
);
2754 if (tu
->_mp_size
< tv
->_mp_size
)
2757 mpz_tdiv_r (tu
, tu
, tv
);
2758 if (tu
->_mp_size
== 0)
2768 c
= mpz_cmp (tu
, tv
);
2777 if (tv
->_mp_size
== 1)
2779 mp_limb_t vl
= tv
->_mp_d
[0];
2780 mp_limb_t ul
= mpz_tdiv_ui (tu
, vl
);
2781 mpz_set_ui (g
, mpn_gcd_11 (ul
, vl
));
2784 mpz_sub (tu
, tu
, tv
);
2788 mpz_mul_2exp (g
, g
, gz
);
2792 mpz_gcdext (mpz_t g
, mpz_t s
, mpz_t t
, const mpz_t u
, const mpz_t v
)
2794 mpz_t tu
, tv
, s0
, s1
, t0
, t1
;
2795 mp_bitcnt_t uz
, vz
, gz
;
2798 if (u
->_mp_size
== 0)
2800 /* g = 0 u + sgn(v) v */
2801 signed long sign
= mpz_sgn (v
);
2806 mpz_set_si (t
, sign
);
2810 if (v
->_mp_size
== 0)
2812 /* g = sgn(u) u + 0 v */
2813 signed long sign
= mpz_sgn (u
);
2816 mpz_set_si (s
, sign
);
2830 uz
= mpz_make_odd (tu
);
2832 vz
= mpz_make_odd (tv
);
2833 gz
= GMP_MIN (uz
, vz
);
2838 /* Cofactors corresponding to odd gcd. gz handled later. */
2839 if (tu
->_mp_size
< tv
->_mp_size
)
2842 MPZ_SRCPTR_SWAP (u
, v
);
2843 MPZ_PTR_SWAP (s
, t
);
2844 MP_BITCNT_T_SWAP (uz
, vz
);
2852 * where u and v denote the inputs with common factors of two
2853 * eliminated, and det (s0, t0; s1, t1) = 2^p. Then
2855 * 2^p tu = s1 u - t1 v
2856 * 2^p tv = -s0 u + t0 v
2859 /* After initial division, tu = q tv + tu', we have
2861 * u = 2^uz (tu' + q tv)
2866 * t0 = 2^uz, t1 = 2^uz q
2870 mpz_setbit (t0
, uz
);
2871 mpz_tdiv_qr (t1
, tu
, tu
, tv
);
2872 mpz_mul_2exp (t1
, t1
, uz
);
2874 mpz_setbit (s1
, vz
);
2877 if (tu
->_mp_size
> 0)
2880 shift
= mpz_make_odd (tu
);
2881 mpz_mul_2exp (t0
, t0
, shift
);
2882 mpz_mul_2exp (s0
, s0
, shift
);
2888 c
= mpz_cmp (tu
, tv
);
2896 * u = t0 tu + t1 (tv' + tu) = (t0 + t1) tu + t1 tv'
2897 * v = s0 tu + s1 (tv' + tu) = (s0 + s1) tu + s1 tv' */
2899 mpz_sub (tv
, tv
, tu
);
2900 mpz_add (t0
, t0
, t1
);
2901 mpz_add (s0
, s0
, s1
);
2903 shift
= mpz_make_odd (tv
);
2904 mpz_mul_2exp (t1
, t1
, shift
);
2905 mpz_mul_2exp (s1
, s1
, shift
);
2909 mpz_sub (tu
, tu
, tv
);
2910 mpz_add (t1
, t0
, t1
);
2911 mpz_add (s1
, s0
, s1
);
2913 shift
= mpz_make_odd (tu
);
2914 mpz_mul_2exp (t0
, t0
, shift
);
2915 mpz_mul_2exp (s0
, s0
, shift
);
2921 /* Now tv = odd part of gcd, and -s0 and t0 are corresponding
2924 mpz_mul_2exp (tv
, tv
, gz
);
2927 /* 2^p g = s0 u + t0 v. Eliminate one factor of two at a time. To
2928 adjust cofactors, we need u / g and v / g */
2930 mpz_divexact (s1
, v
, tv
);
2932 mpz_divexact (t1
, u
, tv
);
2937 /* s0 u + t0 v = (s0 - v/g) u - (t0 + u/g) v */
2938 if (mpz_odd_p (s0
) || mpz_odd_p (t0
))
2940 mpz_sub (s0
, s0
, s1
);
2941 mpz_add (t0
, t0
, t1
);
2943 mpz_divexact_ui (s0
, s0
, 2);
2944 mpz_divexact_ui (t0
, t0
, 2);
2947 /* Arrange so that |s| < |u| / 2g */
2948 mpz_add (s1
, s0
, s1
);
2949 if (mpz_cmpabs (s0
, s1
) > 0)
2952 mpz_sub (t0
, t0
, t1
);
2954 if (u
->_mp_size
< 0)
2956 if (v
->_mp_size
< 0)
2974 mpz_lcm (mpz_t r
, const mpz_t u
, const mpz_t v
)
2978 if (u
->_mp_size
== 0 || v
->_mp_size
== 0)
2987 mpz_divexact (g
, u
, g
);
2995 mpz_lcm_ui (mpz_t r
, const mpz_t u
, unsigned long v
)
2997 if (v
== 0 || u
->_mp_size
== 0)
3003 v
/= mpz_gcd_ui (NULL
, u
, v
);
3004 mpz_mul_ui (r
, u
, v
);
3010 mpz_invert (mpz_t r
, const mpz_t u
, const mpz_t m
)
3015 if (u
->_mp_size
== 0 || mpz_cmpabs_ui (m
, 1) <= 0)
3021 mpz_gcdext (g
, tr
, NULL
, u
, m
);
3022 invertible
= (mpz_cmp_ui (g
, 1) == 0);
3026 if (tr
->_mp_size
< 0)
3028 if (m
->_mp_size
>= 0)
3029 mpz_add (tr
, tr
, m
);
3031 mpz_sub (tr
, tr
, m
);
3042 /* Higher level operations (sqrt, pow and root) */
3045 mpz_pow_ui (mpz_t r
, const mpz_t b
, unsigned long e
)
3049 mpz_init_set_ui (tr
, 1);
3051 bit
= GMP_ULONG_HIGHBIT
;
3054 mpz_mul (tr
, tr
, tr
);
3056 mpz_mul (tr
, tr
, b
);
3066 mpz_ui_pow_ui (mpz_t r
, unsigned long blimb
, unsigned long e
)
3069 mpz_pow_ui (r
, mpz_roinit_n (b
, &blimb
, 1), e
);
3073 mpz_powm (mpz_t r
, const mpz_t b
, const mpz_t e
, const mpz_t m
)
3079 struct gmp_div_inverse minv
;
3083 en
= GMP_ABS (e
->_mp_size
);
3084 mn
= GMP_ABS (m
->_mp_size
);
3086 gmp_die ("mpz_powm: Zero modulo.");
3095 mpn_div_qr_invert (&minv
, mp
, mn
);
3100 /* To avoid shifts, we do all our reductions, except the final
3101 one, using a *normalized* m. */
3104 tp
= gmp_xalloc_limbs (mn
);
3105 gmp_assert_nocarry (mpn_lshift (tp
, mp
, mn
, shift
));
3111 if (e
->_mp_size
< 0)
3113 if (!mpz_invert (base
, b
, m
))
3114 gmp_die ("mpz_powm: Negative exponent and non-invertible base.");
3121 bn
= base
->_mp_size
;
3124 mpn_div_qr_preinv (NULL
, base
->_mp_d
, base
->_mp_size
, mp
, mn
, &minv
);
3128 /* We have reduced the absolute value. Now take care of the
3129 sign. Note that we get zero represented non-canonically as
3131 if (b
->_mp_size
< 0)
3133 mp_ptr bp
= MPZ_REALLOC (base
, mn
);
3134 gmp_assert_nocarry (mpn_sub (bp
, mp
, mn
, bp
, bn
));
3137 base
->_mp_size
= mpn_normalized_size (base
->_mp_d
, bn
);
3139 mpz_init_set_ui (tr
, 1);
3143 mp_limb_t w
= e
->_mp_d
[en
];
3146 bit
= GMP_LIMB_HIGHBIT
;
3149 mpz_mul (tr
, tr
, tr
);
3151 mpz_mul (tr
, tr
, base
);
3152 if (tr
->_mp_size
> mn
)
3154 mpn_div_qr_preinv (NULL
, tr
->_mp_d
, tr
->_mp_size
, mp
, mn
, &minv
);
3155 tr
->_mp_size
= mpn_normalized_size (tr
->_mp_d
, mn
);
3162 /* Final reduction */
3163 if (tr
->_mp_size
>= mn
)
3166 mpn_div_qr_preinv (NULL
, tr
->_mp_d
, tr
->_mp_size
, mp
, mn
, &minv
);
3167 tr
->_mp_size
= mpn_normalized_size (tr
->_mp_d
, mn
);
3178 mpz_powm_ui (mpz_t r
, const mpz_t b
, unsigned long elimb
, const mpz_t m
)
3181 mpz_powm (r
, b
, mpz_roinit_n (e
, &elimb
, 1), m
);
3184 /* x=trunc(y^(1/z)), r=y-x^z */
3186 mpz_rootrem (mpz_t x
, mpz_t r
, const mpz_t y
, unsigned long z
)
3191 sgn
= y
->_mp_size
< 0;
3192 if ((~z
& sgn
) != 0)
3193 gmp_die ("mpz_rootrem: Negative argument, with even root.");
3195 gmp_die ("mpz_rootrem: Zeroth root.");
3197 if (mpz_cmpabs_ui (y
, 1) <= 0) {
3208 tb
= mpz_sizeinbase (y
, 2) / z
+ 1;
3209 mpz_init2 (t
, tb
+ 1);
3213 if (z
== 2) /* simplify sqrt loop: z-1 == 1 */
3215 mpz_swap (u
, t
); /* u = x */
3216 mpz_tdiv_q (t
, y
, u
); /* t = y/x */
3217 mpz_add (t
, t
, u
); /* t = y/x + x */
3218 mpz_tdiv_q_2exp (t
, t
, 1); /* x'= (y/x + x)/2 */
3219 } while (mpz_cmpabs (t
, u
) < 0); /* |x'| < |x| */
3228 mpz_swap (u
, t
); /* u = x */
3229 mpz_pow_ui (t
, u
, z
- 1); /* t = x^(z-1) */
3230 mpz_tdiv_q (t
, y
, t
); /* t = y/x^(z-1) */
3231 mpz_mul_ui (v
, u
, z
- 1); /* v = x*(z-1) */
3232 mpz_add (t
, t
, v
); /* t = y/x^(z-1) + x*(z-1) */
3233 mpz_tdiv_q_ui (t
, t
, z
); /* x'=(y/x^(z-1) + x*(z-1))/z */
3234 } while (mpz_cmpabs (t
, u
) < 0); /* |x'| < |x| */
3240 mpz_pow_ui (t
, u
, z
);
3250 mpz_root (mpz_t x
, const mpz_t y
, unsigned long z
)
3256 mpz_rootrem (x
, r
, y
, z
);
3257 res
= r
->_mp_size
== 0;
3263 /* Compute s = floor(sqrt(u)) and r = u - s^2. Allows r == NULL */
3265 mpz_sqrtrem (mpz_t s
, mpz_t r
, const mpz_t u
)
3267 mpz_rootrem (s
, r
, u
, 2);
3271 mpz_sqrt (mpz_t s
, const mpz_t u
)
3273 mpz_rootrem (s
, NULL
, u
, 2);
3277 mpz_perfect_square_p (const mpz_t u
)
3279 if (u
->_mp_size
<= 0)
3280 return (u
->_mp_size
== 0);
3282 return mpz_root (NULL
, u
, 2);
3286 mpn_perfect_square_p (mp_srcptr p
, mp_size_t n
)
3291 assert (p
[n
-1] != 0);
3292 return mpz_root (NULL
, mpz_roinit_n (t
, p
, n
), 2);
3296 mpn_sqrtrem (mp_ptr sp
, mp_ptr rp
, mp_srcptr p
, mp_size_t n
)
3302 assert (p
[n
-1] != 0);
3306 mpz_rootrem (s
, r
, mpz_roinit_n (u
, p
, n
), 2);
3308 assert (s
->_mp_size
== (n
+1)/2);
3309 mpn_copyd (sp
, s
->_mp_d
, s
->_mp_size
);
3313 mpn_copyd (rp
, r
->_mp_d
, res
);
3321 mpz_fac_ui (mpz_t x
, unsigned long n
)
3323 mpz_set_ui (x
, n
+ (n
== 0));
3325 mpz_mul_ui (x
, x
, --n
);
3329 mpz_bin_uiui (mpz_t r
, unsigned long n
, unsigned long k
)
3333 mpz_set_ui (r
, k
<= n
);
3336 k
= (k
<= n
) ? n
- k
: 0;
3342 mpz_mul_ui (r
, r
, n
--);
3344 mpz_divexact (r
, r
, t
);
3349 /* Primality testing */
3351 gmp_millerrabin (const mpz_t n
, const mpz_t nm1
, mpz_t y
,
3352 const mpz_t q
, mp_bitcnt_t k
)
3356 /* Caller must initialize y to the base. */
3357 mpz_powm (y
, y
, q
, n
);
3359 if (mpz_cmp_ui (y
, 1) == 0 || mpz_cmp (y
, nm1
) == 0)
3364 mpz_powm_ui (y
, y
, 2, n
);
3365 if (mpz_cmp (y
, nm1
) == 0)
3367 /* y == 1 means that the previous y was a non-trivial square root
3368 of 1 (mod n). y == 0 means that n is a power of the base.
3369 In either case, n is not prime. */
3370 if (mpz_cmp_ui (y
, 1) <= 0)
3376 /* This product is 0xc0cfd797, and fits in 32 bits. */
3377 #define GMP_PRIME_PRODUCT \
3378 (3UL*5UL*7UL*11UL*13UL*17UL*19UL*23UL*29UL)
3380 /* Bit (p+1)/2 is set, for each odd prime <= 61 */
3381 #define GMP_PRIME_MASK 0xc96996dcUL
3384 mpz_probab_prime_p (const mpz_t n
, int reps
)
3393 /* Note that we use the absolute value of n only, for compatibility
3394 with the real GMP. */
3396 return (mpz_cmpabs_ui (n
, 2) == 0) ? 2 : 0;
3398 /* Above test excludes n == 0 */
3399 assert (n
->_mp_size
!= 0);
3401 if (mpz_cmpabs_ui (n
, 64) < 0)
3402 return (GMP_PRIME_MASK
>> (n
->_mp_d
[0] >> 1)) & 2;
3404 if (mpz_gcd_ui (NULL
, n
, GMP_PRIME_PRODUCT
) != 1)
3407 /* All prime factors are >= 31. */
3408 if (mpz_cmpabs_ui (n
, 31*31) < 0)
3411 /* Use Miller-Rabin, with a deterministic sequence of bases, a[j] =
3412 j^2 + j + 41 using Euler's polynomial. We potentially stop early,
3413 if a[j] >= n - 1. Since n >= 31*31, this can happen only if reps >
3414 30 (a[30] == 971 > 31*31 == 961). */
3420 /* Find q and k, where q is odd and n = 1 + 2**k * q. */
3421 nm1
->_mp_size
= mpz_abs_sub_ui (nm1
, n
, 1);
3422 k
= mpz_scan1 (nm1
, 0);
3423 mpz_tdiv_q_2exp (q
, nm1
, k
);
3425 for (j
= 0, is_prime
= 1; is_prime
& (j
< reps
); j
++)
3427 mpz_set_ui (y
, (unsigned long) j
*j
+j
+41);
3428 if (mpz_cmp (y
, nm1
) >= 0)
3430 /* Don't try any further bases. This "early" break does not affect
3431 the result for any reasonable reps value (<=5000 was tested) */
3435 is_prime
= gmp_millerrabin (n
, nm1
, y
, q
, k
);
3445 /* Logical operations and bit manipulation. */
3447 /* Numbers are treated as if represented in two's complement (and
3448 infinitely sign extended). For a negative values we get the two's
3449 complement from -x = ~x + 1, where ~ is bitwise complement.
3458 where yyyy is the bitwise complement of xxxx. So least significant
3459 bits, up to and including the first one bit, are unchanged, and
3460 the more significant bits are all complemented.
3462 To change a bit from zero to one in a negative number, subtract the
3463 corresponding power of two from the absolute value. This can never
3464 underflow. To change a bit from one to zero, add the corresponding
3465 power of two, and this might overflow. E.g., if x = -001111, the
3466 two's complement is 110001. Clearing the least significant bit, we
3467 get two's complement 110000, and -010000. */
3470 mpz_tstbit (const mpz_t d
, mp_bitcnt_t bit_index
)
3472 mp_size_t limb_index
;
3481 limb_index
= bit_index
/ GMP_LIMB_BITS
;
3482 if (limb_index
>= dn
)
3485 shift
= bit_index
% GMP_LIMB_BITS
;
3486 w
= d
->_mp_d
[limb_index
];
3487 bit
= (w
>> shift
) & 1;
3491 /* d < 0. Check if any of the bits below is set: If so, our bit
3492 must be complemented. */
3493 if (shift
> 0 && (w
<< (GMP_LIMB_BITS
- shift
)) > 0)
3495 while (--limb_index
>= 0)
3496 if (d
->_mp_d
[limb_index
] > 0)
3503 mpz_abs_add_bit (mpz_t d
, mp_bitcnt_t bit_index
)
3505 mp_size_t dn
, limb_index
;
3509 dn
= GMP_ABS (d
->_mp_size
);
3511 limb_index
= bit_index
/ GMP_LIMB_BITS
;
3512 bit
= (mp_limb_t
) 1 << (bit_index
% GMP_LIMB_BITS
);
3514 if (limb_index
>= dn
)
3517 /* The bit should be set outside of the end of the number.
3518 We have to increase the size of the number. */
3519 dp
= MPZ_REALLOC (d
, limb_index
+ 1);
3521 dp
[limb_index
] = bit
;
3522 for (i
= dn
; i
< limb_index
; i
++)
3524 dn
= limb_index
+ 1;
3532 cy
= mpn_add_1 (dp
+ limb_index
, dp
+ limb_index
, dn
- limb_index
, bit
);
3535 dp
= MPZ_REALLOC (d
, dn
+ 1);
3540 d
->_mp_size
= (d
->_mp_size
< 0) ? - dn
: dn
;
3544 mpz_abs_sub_bit (mpz_t d
, mp_bitcnt_t bit_index
)
3546 mp_size_t dn
, limb_index
;
3550 dn
= GMP_ABS (d
->_mp_size
);
3553 limb_index
= bit_index
/ GMP_LIMB_BITS
;
3554 bit
= (mp_limb_t
) 1 << (bit_index
% GMP_LIMB_BITS
);
3556 assert (limb_index
< dn
);
3558 gmp_assert_nocarry (mpn_sub_1 (dp
+ limb_index
, dp
+ limb_index
,
3559 dn
- limb_index
, bit
));
3560 dn
= mpn_normalized_size (dp
, dn
);
3561 d
->_mp_size
= (d
->_mp_size
< 0) ? - dn
: dn
;
3565 mpz_setbit (mpz_t d
, mp_bitcnt_t bit_index
)
3567 if (!mpz_tstbit (d
, bit_index
))
3569 if (d
->_mp_size
>= 0)
3570 mpz_abs_add_bit (d
, bit_index
);
3572 mpz_abs_sub_bit (d
, bit_index
);
3577 mpz_clrbit (mpz_t d
, mp_bitcnt_t bit_index
)
3579 if (mpz_tstbit (d
, bit_index
))
3581 if (d
->_mp_size
>= 0)
3582 mpz_abs_sub_bit (d
, bit_index
);
3584 mpz_abs_add_bit (d
, bit_index
);
3589 mpz_combit (mpz_t d
, mp_bitcnt_t bit_index
)
3591 if (mpz_tstbit (d
, bit_index
) ^ (d
->_mp_size
< 0))
3592 mpz_abs_sub_bit (d
, bit_index
);
3594 mpz_abs_add_bit (d
, bit_index
);
3598 mpz_com (mpz_t r
, const mpz_t u
)
3601 mpz_sub_ui (r
, r
, 1);
3605 mpz_and (mpz_t r
, const mpz_t u
, const mpz_t v
)
3607 mp_size_t un
, vn
, rn
, i
;
3610 mp_limb_t ux
, vx
, rx
;
3611 mp_limb_t uc
, vc
, rc
;
3612 mp_limb_t ul
, vl
, rl
;
3614 un
= GMP_ABS (u
->_mp_size
);
3615 vn
= GMP_ABS (v
->_mp_size
);
3618 MPZ_SRCPTR_SWAP (u
, v
);
3619 MP_SIZE_T_SWAP (un
, vn
);
3627 uc
= u
->_mp_size
< 0;
3628 vc
= v
->_mp_size
< 0;
3635 /* If the smaller input is positive, higher limbs don't matter. */
3638 rp
= MPZ_REALLOC (r
, rn
+ rc
);
3646 ul
= (up
[i
] ^ ux
) + uc
;
3649 vl
= (vp
[i
] ^ vx
) + vc
;
3652 rl
= ( (ul
& vl
) ^ rx
) + rc
;
3661 ul
= (up
[i
] ^ ux
) + uc
;
3664 rl
= ( (ul
& vx
) ^ rx
) + rc
;
3671 rn
= mpn_normalized_size (rp
, rn
);
3673 r
->_mp_size
= rx
? -rn
: rn
;
3677 mpz_ior (mpz_t r
, const mpz_t u
, const mpz_t v
)
3679 mp_size_t un
, vn
, rn
, i
;
3682 mp_limb_t ux
, vx
, rx
;
3683 mp_limb_t uc
, vc
, rc
;
3684 mp_limb_t ul
, vl
, rl
;
3686 un
= GMP_ABS (u
->_mp_size
);
3687 vn
= GMP_ABS (v
->_mp_size
);
3690 MPZ_SRCPTR_SWAP (u
, v
);
3691 MP_SIZE_T_SWAP (un
, vn
);
3699 uc
= u
->_mp_size
< 0;
3700 vc
= v
->_mp_size
< 0;
3707 /* If the smaller input is negative, by sign extension higher limbs
3711 rp
= MPZ_REALLOC (r
, rn
+ rc
);
3719 ul
= (up
[i
] ^ ux
) + uc
;
3722 vl
= (vp
[i
] ^ vx
) + vc
;
3725 rl
= ( (ul
| vl
) ^ rx
) + rc
;
3734 ul
= (up
[i
] ^ ux
) + uc
;
3737 rl
= ( (ul
| vx
) ^ rx
) + rc
;
3744 rn
= mpn_normalized_size (rp
, rn
);
3746 r
->_mp_size
= rx
? -rn
: rn
;
3750 mpz_xor (mpz_t r
, const mpz_t u
, const mpz_t v
)
3752 mp_size_t un
, vn
, i
;
3755 mp_limb_t ux
, vx
, rx
;
3756 mp_limb_t uc
, vc
, rc
;
3757 mp_limb_t ul
, vl
, rl
;
3759 un
= GMP_ABS (u
->_mp_size
);
3760 vn
= GMP_ABS (v
->_mp_size
);
3763 MPZ_SRCPTR_SWAP (u
, v
);
3764 MP_SIZE_T_SWAP (un
, vn
);
3772 uc
= u
->_mp_size
< 0;
3773 vc
= v
->_mp_size
< 0;
3780 rp
= MPZ_REALLOC (r
, un
+ rc
);
3788 ul
= (up
[i
] ^ ux
) + uc
;
3791 vl
= (vp
[i
] ^ vx
) + vc
;
3794 rl
= (ul
^ vl
^ rx
) + rc
;
3803 ul
= (up
[i
] ^ ux
) + uc
;
3806 rl
= (ul
^ ux
) + rc
;
3813 un
= mpn_normalized_size (rp
, un
);
3815 r
->_mp_size
= rx
? -un
: un
;
3819 gmp_popcount_limb (mp_limb_t x
)
3823 /* Do 16 bits at a time, to avoid limb-sized constants. */
3824 for (c
= 0; x
> 0; x
>>= 16)
3826 unsigned w
= ((x
>> 1) & 0x5555) + (x
& 0x5555);
3827 w
= ((w
>> 2) & 0x3333) + (w
& 0x3333);
3828 w
= ((w
>> 4) & 0x0f0f) + (w
& 0x0f0f);
3829 w
= (w
>> 8) + (w
& 0x00ff);
3836 mpn_popcount (mp_srcptr p
, mp_size_t n
)
3841 for (c
= 0, i
= 0; i
< n
; i
++)
3842 c
+= gmp_popcount_limb (p
[i
]);
3848 mpz_popcount (const mpz_t u
)
3855 return ~(mp_bitcnt_t
) 0;
3857 return mpn_popcount (u
->_mp_d
, un
);
3861 mpz_hamdist (const mpz_t u
, const mpz_t v
)
3863 mp_size_t un
, vn
, i
;
3864 mp_limb_t uc
, vc
, ul
, vl
, comp
;
3872 return ~(mp_bitcnt_t
) 0;
3874 comp
= - (uc
= vc
= (un
< 0));
3886 MPN_SRCPTR_SWAP (up
, un
, vp
, vn
);
3888 for (i
= 0, c
= 0; i
< vn
; i
++)
3890 ul
= (up
[i
] ^ comp
) + uc
;
3893 vl
= (vp
[i
] ^ comp
) + vc
;
3896 c
+= gmp_popcount_limb (ul
^ vl
);
3902 ul
= (up
[i
] ^ comp
) + uc
;
3905 c
+= gmp_popcount_limb (ul
^ comp
);
3912 mpz_scan1 (const mpz_t u
, mp_bitcnt_t starting_bit
)
3915 mp_size_t us
, un
, i
;
3920 i
= starting_bit
/ GMP_LIMB_BITS
;
3922 /* Past the end there's no 1 bits for u>=0, or an immediate 1 bit
3923 for u<0. Notice this test picks up any u==0 too. */
3925 return (us
>= 0 ? ~(mp_bitcnt_t
) 0 : starting_bit
);
3931 if (starting_bit
!= 0)
3935 ux
= mpn_zero_p (up
, i
);
3937 ux
= - (mp_limb_t
) (limb
>= ux
);
3940 /* Mask to 0 all bits before starting_bit, thus ignoring them. */
3941 limb
&= (GMP_LIMB_MAX
<< (starting_bit
% GMP_LIMB_BITS
));
3944 return mpn_common_scan (limb
, i
, up
, un
, ux
);
3948 mpz_scan0 (const mpz_t u
, mp_bitcnt_t starting_bit
)
3951 mp_size_t us
, un
, i
;
3955 ux
= - (mp_limb_t
) (us
>= 0);
3957 i
= starting_bit
/ GMP_LIMB_BITS
;
3959 /* When past end, there's an immediate 0 bit for u>=0, or no 0 bits for
3960 u<0. Notice this test picks up all cases of u==0 too. */
3962 return (ux
? starting_bit
: ~(mp_bitcnt_t
) 0);
3968 limb
-= mpn_zero_p (up
, i
); /* limb = ~(~limb + zero_p) */
3970 /* Mask all bits before starting_bit, thus ignoring them. */
3971 limb
&= (GMP_LIMB_MAX
<< (starting_bit
% GMP_LIMB_BITS
));
3973 return mpn_common_scan (limb
, i
, up
, un
, ux
);
3977 /* MPZ base conversion. */
3980 mpz_sizeinbase (const mpz_t u
, int base
)
3986 struct gmp_div_inverse bi
;
3990 assert (base
<= 36);
3992 un
= GMP_ABS (u
->_mp_size
);
3998 bits
= (un
- 1) * GMP_LIMB_BITS
+ mpn_limb_size_in_base_2 (up
[un
-1]);
4004 return (bits
+ 1) / 2;
4006 return (bits
+ 2) / 3;
4008 return (bits
+ 3) / 4;
4010 return (bits
+ 4) / 5;
4011 /* FIXME: Do something more clever for the common case of base
4015 tp
= gmp_xalloc_limbs (un
);
4016 mpn_copyi (tp
, up
, un
);
4017 mpn_div_qr_1_invert (&bi
, base
);
4023 mpn_div_qr_1_preinv (tp
, tp
, un
, &bi
);
4024 un
-= (tp
[un
-1] == 0);
4033 mpz_get_str (char *sp
, int base
, const mpz_t u
)
4042 digits
= "0123456789abcdefghijklmnopqrstuvwxyz";
4047 digits
= "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
4054 sn
= 1 + mpz_sizeinbase (u
, base
);
4056 sp
= (char *) gmp_xalloc (1 + sn
);
4058 un
= GMP_ABS (u
->_mp_size
);
4069 if (u
->_mp_size
< 0)
4072 bits
= mpn_base_power_of_two_p (base
);
4075 /* Not modified in this case. */
4076 sn
= i
+ mpn_get_str_bits ((unsigned char *) sp
+ i
, bits
, u
->_mp_d
, un
);
4079 struct mpn_base_info info
;
4082 mpn_get_base_info (&info
, base
);
4083 tp
= gmp_xalloc_limbs (un
);
4084 mpn_copyi (tp
, u
->_mp_d
, un
);
4086 sn
= i
+ mpn_get_str_other ((unsigned char *) sp
+ i
, base
, &info
, tp
, un
);
4091 sp
[i
] = digits
[(unsigned char) sp
[i
]];
4098 mpz_set_str (mpz_t r
, const char *sp
, int base
)
4101 mp_size_t rn
, alloc
;
4107 assert (base
== 0 || (base
>= 2 && base
<= 36));
4109 while (isspace( (unsigned char) *sp
))
4112 sign
= (*sp
== '-');
4120 if (*sp
== 'x' || *sp
== 'X')
4125 else if (*sp
== 'b' || *sp
== 'B')
4138 dp
= (unsigned char *) gmp_xalloc (sn
+ (sn
== 0));
4140 for (sn
= 0; *sp
; sp
++)
4144 if (isspace ((unsigned char) *sp
))
4146 if (*sp
>= '0' && *sp
<= '9')
4148 else if (*sp
>= 'a' && *sp
<= 'z')
4149 digit
= *sp
- 'a' + 10;
4150 else if (*sp
>= 'A' && *sp
<= 'Z')
4151 digit
= *sp
- 'A' + 10;
4153 digit
= base
; /* fail */
4165 bits
= mpn_base_power_of_two_p (base
);
4169 alloc
= (sn
* bits
+ GMP_LIMB_BITS
- 1) / GMP_LIMB_BITS
;
4170 rp
= MPZ_REALLOC (r
, alloc
);
4171 rn
= mpn_set_str_bits (rp
, dp
, sn
, bits
);
4175 struct mpn_base_info info
;
4176 mpn_get_base_info (&info
, base
);
4177 alloc
= (sn
+ info
.exp
- 1) / info
.exp
;
4178 rp
= MPZ_REALLOC (r
, alloc
);
4179 rn
= mpn_set_str_other (rp
, dp
, sn
, base
, &info
);
4181 assert (rn
<= alloc
);
4184 r
->_mp_size
= sign
? - rn
: rn
;
4190 mpz_init_set_str (mpz_t r
, const char *sp
, int base
)
4193 return mpz_set_str (r
, sp
, base
);
4197 mpz_out_str (FILE *stream
, int base
, const mpz_t x
)
4202 str
= mpz_get_str (NULL
, base
, x
);
4204 len
= fwrite (str
, 1, len
, stream
);
4211 gmp_detect_endian (void)
4213 static const int i
= 2;
4214 const unsigned char *p
= (const unsigned char *) &i
;
4218 /* Import and export. Does not support nails. */
4220 mpz_import (mpz_t r
, size_t count
, int order
, size_t size
, int endian
,
4221 size_t nails
, const void *src
)
4223 const unsigned char *p
;
4224 ptrdiff_t word_step
;
4228 /* The current (partial) limb. */
4230 /* The number of bytes already copied to this limb (starting from
4233 /* The index where the limb should be stored, when completed. */
4237 gmp_die ("mpz_import: Nails not supported.");
4239 assert (order
== 1 || order
== -1);
4240 assert (endian
>= -1 && endian
<= 1);
4243 endian
= gmp_detect_endian ();
4245 p
= (unsigned char *) src
;
4247 word_step
= (order
!= endian
) ? 2 * size
: 0;
4249 /* Process bytes from the least significant end, so point p at the
4250 least significant word. */
4253 p
+= size
* (count
- 1);
4254 word_step
= - word_step
;
4257 /* And at least significant byte of that word. */
4261 rn
= (size
* count
+ sizeof(mp_limb_t
) - 1) / sizeof(mp_limb_t
);
4262 rp
= MPZ_REALLOC (r
, rn
);
4264 for (limb
= 0, bytes
= 0, i
= 0; count
> 0; count
--, p
+= word_step
)
4267 for (j
= 0; j
< size
; j
++, p
-= (ptrdiff_t) endian
)
4269 limb
|= (mp_limb_t
) *p
<< (bytes
++ * CHAR_BIT
);
4270 if (bytes
== sizeof(mp_limb_t
))
4278 assert (i
+ (bytes
> 0) == rn
);
4282 i
= mpn_normalized_size (rp
, i
);
4288 mpz_export (void *r
, size_t *countp
, int order
, size_t size
, int endian
,
4289 size_t nails
, const mpz_t u
)
4295 gmp_die ("mpz_import: Nails not supported.");
4297 assert (order
== 1 || order
== -1);
4298 assert (endian
>= -1 && endian
<= 1);
4299 assert (size
> 0 || u
->_mp_size
== 0);
4307 ptrdiff_t word_step
;
4308 /* The current (partial) limb. */
4310 /* The number of bytes left to to in this limb. */
4312 /* The index where the limb was read. */
4317 /* Count bytes in top limb. */
4318 limb
= u
->_mp_d
[un
-1];
4323 k
++; limb
>>= CHAR_BIT
;
4324 } while (limb
!= 0);
4326 count
= (k
+ (un
-1) * sizeof (mp_limb_t
) + size
- 1) / size
;
4329 r
= gmp_xalloc (count
* size
);
4332 endian
= gmp_detect_endian ();
4334 p
= (unsigned char *) r
;
4336 word_step
= (order
!= endian
) ? 2 * size
: 0;
4338 /* Process bytes from the least significant end, so point p at the
4339 least significant word. */
4342 p
+= size
* (count
- 1);
4343 word_step
= - word_step
;
4346 /* And at least significant byte of that word. */
4350 for (bytes
= 0, i
= 0, k
= 0; k
< count
; k
++, p
+= word_step
)
4353 for (j
= 0; j
< size
; j
++, p
-= (ptrdiff_t) endian
)
4358 limb
= u
->_mp_d
[i
++];
4359 bytes
= sizeof (mp_limb_t
);
4367 assert (k
== count
);