nettle  3.7.3
About: Nettle is a low-level cryptographic library.
  Fossies Dox: nettle-3.7.3.tar.gz  ("unofficial" and yet experimental doxygen-generated source code documentation)  

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