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