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