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