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