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