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