diff -r b641c8181188 .hgignore --- a/.hgignore Wed Feb 15 14:46:40 2012 +0100 +++ b/.hgignore Fri Feb 17 11:16:53 2012 +0100 @@ -56,6 +56,22 @@ Makefile ^tests/mpz/logic$ ^tests/mpz/reuse$ +^mini-gmp/tests/t-add$ +^mini-gmp/tests/t-sub$ +^mini-gmp/tests/t-mul$ +^mini-gmp/tests/t-invert$ +^mini-gmp/tests/t-div$ +^mini-gmp/tests/t-div_2exp$ +^mini-gmp/tests/t-double$ +^mini-gmp/tests/t-gcd$ +^mini-gmp/tests/t-lcm$ +^mini-gmp/tests/t-sqrt +^mini-gmp/tests/t-powm$ +^mini-gmp/tests/t-logops$ +^mini-gmp/tests/t-bitops$ +^mini-gmp/tests/t-scan$ +^mini-gmp/tests/t-str$ + ^\.libs .*\.a$ diff -r b641c8181188 Makefile.am --- a/Makefile.am Wed Feb 15 14:46:40 2012 +0100 +++ b/Makefile.am Fri Feb 17 11:16:53 2012 +0100 @@ -309,13 +309,13 @@ install-data-hook: # the .h files are not properly expressed for the various objects that use # them. -EXTRA_DIST += dumbmp.c +EXTRA_DIST += bootstrap.c fac_ui.h: gen-fac_ui$(EXEEXT_FOR_BUILD) ./gen-fac_ui $(GMP_LIMB_BITS) $(GMP_NAIL_BITS) >fac_ui.h || (rm -f fac_ui.h; exit 1) BUILT_SOURCES += fac_ui.h -gen-fac_ui$(EXEEXT_FOR_BUILD): gen-fac_ui$(U_FOR_BUILD).c dumbmp.c +gen-fac_ui$(EXEEXT_FOR_BUILD): gen-fac_ui$(U_FOR_BUILD).c bootstrap.c $(CC_FOR_BUILD) `test -f 'gen-fac_ui$(U_FOR_BUILD).c' || echo '$(srcdir)/'`gen-fac_ui$(U_FOR_BUILD).c -o gen-fac_ui$(EXEEXT_FOR_BUILD) DISTCLEANFILES += gen-fac_ui$(EXEEXT_FOR_BUILD) EXTRA_DIST += gen-fac_ui.c @@ -329,7 +329,7 @@ mpn/fib_table.c: gen-fib$(EXEEXT_FOR_BUI ./gen-fib table $(GMP_LIMB_BITS) $(GMP_NAIL_BITS) >mpn/fib_table.c || (rm -f mpn/fib_table.c; exit 1) BUILT_SOURCES += mpn/fib_table.c -gen-fib$(EXEEXT_FOR_BUILD): gen-fib$(U_FOR_BUILD).c dumbmp.c +gen-fib$(EXEEXT_FOR_BUILD): gen-fib$(U_FOR_BUILD).c bootstrap.c $(CC_FOR_BUILD) `test -f 'gen-fib$(U_FOR_BUILD).c' || echo '$(srcdir)/'`gen-fib$(U_FOR_BUILD).c -o gen-fib$(EXEEXT_FOR_BUILD) DISTCLEANFILES += gen-fib$(EXEEXT_FOR_BUILD) EXTRA_DIST += gen-fib.c @@ -343,7 +343,7 @@ mpn/mp_bases.c: gen-bases$(EXEEXT_FOR_BU ./gen-bases table $(GMP_LIMB_BITS) $(GMP_NAIL_BITS) >mpn/mp_bases.c || (rm -f mpn/mp_bases.c; exit 1) BUILT_SOURCES += mpn/mp_bases.c -gen-bases$(EXEEXT_FOR_BUILD): gen-bases$(U_FOR_BUILD).c dumbmp.c +gen-bases$(EXEEXT_FOR_BUILD): gen-bases$(U_FOR_BUILD).c bootstrap.c $(CC_FOR_BUILD) `test -f 'gen-bases$(U_FOR_BUILD).c' || echo '$(srcdir)/'`gen-bases$(U_FOR_BUILD).c -o gen-bases$(EXEEXT_FOR_BUILD) $(LIBM_FOR_BUILD) DISTCLEANFILES += gen-bases$(EXEEXT_FOR_BUILD) EXTRA_DIST += gen-bases.c @@ -353,7 +353,7 @@ trialdivtab.h: gen-trialdivtab$(EXEEXT_F ./gen-trialdivtab $(GMP_LIMB_BITS) 8000 >trialdivtab.h || (rm -f trialdivtab.h; exit 1) BUILT_SOURCES += trialdivtab.h -gen-trialdivtab$(EXEEXT_FOR_BUILD): gen-trialdivtab$(U_FOR_BUILD).c dumbmp.c +gen-trialdivtab$(EXEEXT_FOR_BUILD): gen-trialdivtab$(U_FOR_BUILD).c bootstrap.c $(CC_FOR_BUILD) `test -f 'gen-trialdivtab$(U_FOR_BUILD).c' || echo '$(srcdir)/'`gen-trialdivtab$(U_FOR_BUILD).c -o gen-trialdivtab$(EXEEXT_FOR_BUILD) $(LIBM_FOR_BUILD) DISTCLEANFILES += gen-trialdivtab$(EXEEXT_FOR_BUILD) EXTRA_DIST += gen-trialdivtab.c @@ -373,12 +373,14 @@ mpn/perfsqr.h: gen-psqr$(EXEEXT_FOR_BUIL ./gen-psqr $(GMP_LIMB_BITS) $(GMP_NAIL_BITS) >mpn/perfsqr.h || (rm -f mpn/perfsqr.h; exit 1) BUILT_SOURCES += mpn/perfsqr.h -gen-psqr$(EXEEXT_FOR_BUILD): gen-psqr$(U_FOR_BUILD).c dumbmp.c +gen-psqr$(EXEEXT_FOR_BUILD): gen-psqr$(U_FOR_BUILD).c bootstrap.c $(CC_FOR_BUILD) `test -f 'gen-psqr$(U_FOR_BUILD).c' || echo '$(srcdir)/'`gen-psqr$(U_FOR_BUILD).c -o gen-psqr$(EXEEXT_FOR_BUILD) $(LIBM_FOR_BUILD) DISTCLEANFILES += gen-psqr$(EXEEXT_FOR_BUILD) EXTRA_DIST += gen-psqr.c - +# Distribute mini-gmp. Test sources copied by dist-hook. +EXTRA_DIST += mini-gmp/README mini-gmp/mini-gmp.c mini-gmp/mini-gmp.h \ + mini-gmp/tests/Makefile mini-gmp/tests/run-tests # Avoid: CVS - cvs directories # *~ - emacs backups @@ -390,6 +392,7 @@ EXTRA_DIST += gen-psqr.c dist-hook: -find $(distdir) \( -name CVS -type d \) -o -name "*~" -o -name ".#*" \ | xargs rm -rf + cp "$(srcdir)"/mini-gmp/tests/*.[ch] "$(distdir)/mini-gmp/tests" # grep -F $(VERSION) $(srcdir)/Makefile.am \ # | grep -q "^# *$(VERSION) *$(LIBGMP_LT_CURRENT):$(LIBGMP_LT_REVISION):$(LIBGMP_LT_AGE) *$(LIBGMPXX_LT_CURRENT):$(LIBGMPXX_LT_REVISION):$(LIBGMPXX_LT_AGE) *$(LIBMP_LT_CURRENT):$(LIBMP_LT_REVISION):$(LIBMP_LT_AGE)" # test -z "`sed -n 's/^# *[0-9]*\.[0-9]*\.[0-9]* *\([0-9]*:[0-9]*:[0-9]*\) *\([0-9]*:[0-9]*:[0-9]*\) *\([0-9]*:[0-9]*:[0-9]*\).*/A\1\nB\2\nC\3/p' $(srcdir)/Makefile.am | grep -v 'A6:3:3\|B3:5:0\|C4:7:1' | sort | uniq -d`" diff -r b641c8181188 bootstrap.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/bootstrap.c Fri Feb 17 11:16:53 2012 +0100 @@ -0,0 +1,135 @@ +/* Functions needed for bootstrapping the gmp build, based on mini-gmp. + +Copyright 2001, 2002, 2004, 2011, 2012 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 3 of the License, or (at your +option) any later version. + +The GNU MP Library is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public +License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ + + +#include "mini-gmp/mini-gmp.c" + +#define MIN(l,o) ((l) < (o) ? (l) : (o)) +#define PTR(x) ((x)->_mp_d) +#define SIZ(x) ((x)->_mp_size) + +#define xmalloc gmp_default_alloc + +int +isprime (unsigned long int t) +{ + unsigned long int q, r, d; + + if (t < 32) + return (0xa08a28acUL >> t) & 1; + if ((t & 1) == 0) + return 0; + + if (t % 3 == 0) + return 0; + if (t % 5 == 0) + return 0; + if (t % 7 == 0) + return 0; + + for (d = 11;;) + { + q = t / d; + r = t - q * d; + if (q < d) + return 1; + if (r == 0) + break; + d += 2; + q = t / d; + r = t - q * d; + if (q < d) + return 1; + if (r == 0) + break; + d += 4; + } + return 0; +} + +int +log2_ceil (int n) +{ + int e; + assert (n >= 1); + for (e = 0; ; e++) + if ((1 << e) >= n) + break; + return e; +} + +/* Set inv to the inverse of d, in the style of invert_limb, ie. for + udiv_qrnnd_preinv. */ +void +mpz_preinv_invert (mpz_t inv, mpz_t d, int numb_bits) +{ + mpz_t t; + int norm; + assert (SIZ(d) > 0); + + norm = numb_bits - mpz_sizeinbase (d, 2); + assert (norm >= 0); + mpz_init_set_ui (t, 1L); + mpz_mul_2exp (t, t, 2*numb_bits - norm); + mpz_tdiv_q (inv, t, d); + mpz_set_ui (t, 1L); + mpz_mul_2exp (t, t, numb_bits); + mpz_sub (inv, inv, t); + + mpz_clear (t); +} + +/* Calculate r satisfying r*d == 1 mod 2^n. */ +void +mpz_invert_2exp (mpz_t r, mpz_t a, unsigned long n) +{ + unsigned long i; + mpz_t inv, prod; + + assert (mpz_odd_p (a)); + + mpz_init_set_ui (inv, 1L); + mpz_init (prod); + + for (i = 1; i < n; i++) + { + mpz_mul (prod, inv, a); + if (mpz_tstbit (prod, i) != 0) + mpz_setbit (inv, i); + } + + mpz_mul (prod, inv, a); + mpz_tdiv_r_2exp (prod, prod, n); + assert (mpz_cmp_ui (prod, 1L) == 0); + + mpz_set (r, inv); + + mpz_clear (inv); + mpz_clear (prod); +} + +/* Calculate inv satisfying r*a == 1 mod 2^n. */ +void +mpz_invert_ui_2exp (mpz_t r, unsigned long a, unsigned long n) +{ + mpz_t az; + mpz_init_set_ui (az, a); + mpz_invert_2exp (r, az, n); + mpz_clear (az); +} diff -r b641c8181188 dumbmp.c --- a/dumbmp.c Wed Feb 15 14:46:40 2012 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,990 +0,0 @@ -/* dumbmp mini GMP compatible library. - -Copyright 2001, 2002, 2004, 2011 Free Software Foundation, Inc. - -This file is part of the GNU MP Library. - -The GNU MP Library is free software; you can redistribute it and/or modify -it under the terms of the GNU Lesser General Public License as published by -the Free Software Foundation; either version 3 of the License, or (at your -option) any later version. - -The GNU MP Library is distributed in the hope that it will be useful, but -WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY -or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public -License for more details. - -You should have received a copy of the GNU Lesser General Public License -along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ - - -/* The code here implements a subset (a very limited subset) of the main GMP - functions. It's designed for use in a few build-time calculations and - will be slow, but highly portable. - - None of the normal GMP configure things are used, nor any of the normal - gmp.h or gmp-impl.h. To use this file in a program just #include - "dumbmp.c". - - mp_limb_t here is an unsigned long, since that's a sensible type - everywhere we know of, with 8*sizeof(unsigned long) giving the number of - bits in the type (that not being true for instance with int or short on - Cray vector systems.) - - Only the low half of each mp_limb_t is used, so as to make carry handling - and limb multiplies easy. GMP_LIMB_BITS is the number of bits used. */ - -#include -#include -#include - - -typedef unsigned long mp_limb_t; - -typedef struct { - int _mp_alloc; - int _mp_size; - mp_limb_t *_mp_d; -} mpz_t[1]; - -#define GMP_LIMB_BITS (sizeof (mp_limb_t) * 8 / 2) - -#define ABS(x) ((x) >= 0 ? (x) : -(x)) -#define MIN(l,o) ((l) < (o) ? (l) : (o)) -#define MAX(h,i) ((h) > (i) ? (h) : (i)) - -#define ALLOC(x) ((x)->_mp_alloc) -#define PTR(x) ((x)->_mp_d) -#define SIZ(x) ((x)->_mp_size) -#define ABSIZ(x) ABS (SIZ (x)) -#define LOMASK ((1L << GMP_LIMB_BITS) - 1) -#define LO(x) ((x) & LOMASK) -#define HI(x) ((x) >> GMP_LIMB_BITS) - -#define ASSERT(cond) \ - do { \ - if (! (cond)) \ - { \ - fprintf (stderr, "Assertion failure\n"); \ - abort (); \ - } \ - } while (0) - - -char * -xmalloc (int n) -{ - char *p; - p = malloc (n); - if (p == NULL) - { - fprintf (stderr, "Out of memory (alloc %d bytes)\n", n); - abort (); - } - return p; -} - -mp_limb_t * -xmalloc_limbs (int n) -{ - return (mp_limb_t *) xmalloc (n * sizeof (mp_limb_t)); -} - -void -mem_copyi (char *dst, char *src, int size) -{ - int i; - for (i = 0; i < size; i++) - dst[i] = src[i]; -} - -static int -isprime (unsigned long int t) -{ - unsigned long int q, r, d; - - if (t < 32) - return (0xa08a28acUL >> t) & 1; - if ((t & 1) == 0) - return 0; - - if (t % 3 == 0) - return 0; - if (t % 5 == 0) - return 0; - if (t % 7 == 0) - return 0; - - for (d = 11;;) - { - q = t / d; - r = t - q * d; - if (q < d) - return 1; - if (r == 0) - break; - d += 2; - q = t / d; - r = t - q * d; - if (q < d) - return 1; - if (r == 0) - break; - d += 4; - } - return 0; -} - -int -log2_ceil (int n) -{ - int e; - ASSERT (n >= 1); - for (e = 0; ; e++) - if ((1 << e) >= n) - break; - return e; -} - -void -mpz_realloc (mpz_t r, int n) -{ - if (n <= ALLOC(r)) - return; - - ALLOC(r) = n; - PTR(r) = (mp_limb_t *) realloc (PTR(r), n * sizeof (mp_limb_t)); - if (PTR(r) == NULL) - { - fprintf (stderr, "Out of memory (realloc to %d)\n", n); - abort (); - } -} - -void -mpn_normalize (mp_limb_t *rp, int *rnp) -{ - int rn = *rnp; - while (rn > 0 && rp[rn-1] == 0) - rn--; - *rnp = rn; -} - -void -mpn_copyi (mp_limb_t *dst, mp_limb_t *src, int n) -{ - int i; - for (i = 0; i < n; i++) - dst[i] = src[i]; -} - -void -mpn_zero (mp_limb_t *rp, int rn) -{ - int i; - for (i = 0; i < rn; i++) - rp[i] = 0; -} - -void -mpz_init (mpz_t r) -{ - ALLOC(r) = 1; - PTR(r) = xmalloc_limbs (ALLOC(r)); - PTR(r)[0] = 0; - SIZ(r) = 0; -} - -void -mpz_clear (mpz_t r) -{ - free (PTR (r)); - ALLOC(r) = -1; - SIZ (r) = 0xbadcafeL; - PTR (r) = (mp_limb_t *) 0xdeadbeefL; -} - -int -mpz_sgn (mpz_t a) -{ - return (SIZ(a) > 0 ? 1 : SIZ(a) == 0 ? 0 : -1); -} - -int -mpz_odd_p (mpz_t a) -{ - if (SIZ(a) == 0) - return 0; - else - return (PTR(a)[0] & 1) != 0; -} - -int -mpz_even_p (mpz_t a) -{ - if (SIZ(a) == 0) - return 1; - else - return (PTR(a)[0] & 1) == 0; -} - -size_t -mpz_sizeinbase (mpz_t a, int base) -{ - int an = ABSIZ (a); - mp_limb_t *ap = PTR (a); - int cnt; - mp_limb_t hi; - - if (base != 2) - abort (); - - if (an == 0) - return 1; - - cnt = 0; - for (hi = ap[an - 1]; hi != 0; hi >>= 1) - cnt += 1; - return (an - 1) * GMP_LIMB_BITS + cnt; -} - -void -mpz_set (mpz_t r, mpz_t a) -{ - mpz_realloc (r, ABSIZ (a)); - SIZ(r) = SIZ(a); - mpn_copyi (PTR(r), PTR(a), ABSIZ (a)); -} - -void -mpz_init_set (mpz_t r, mpz_t a) -{ - mpz_init (r); - mpz_set (r, a); -} - -void -mpz_set_ui (mpz_t r, unsigned long ui) -{ - int rn; - mpz_realloc (r, 2); - PTR(r)[0] = LO(ui); - PTR(r)[1] = HI(ui); - rn = 2; - mpn_normalize (PTR(r), &rn); - SIZ(r) = rn; -} - -void -mpz_init_set_ui (mpz_t r, unsigned long ui) -{ - mpz_init (r); - mpz_set_ui (r, ui); -} - -void -mpz_setbit (mpz_t r, unsigned long bit) -{ - int limb, rn, extend; - mp_limb_t *rp; - - rn = SIZ(r); - if (rn < 0) - abort (); /* only r>=0 */ - - limb = bit / GMP_LIMB_BITS; - bit %= GMP_LIMB_BITS; - - mpz_realloc (r, limb+1); - rp = PTR(r); - extend = (limb+1) - rn; - if (extend > 0) - mpn_zero (rp + rn, extend); - - rp[limb] |= (mp_limb_t) 1 << bit; - SIZ(r) = MAX (rn, limb+1); -} - -int -mpz_tstbit (mpz_t r, unsigned long bit) -{ - int limb; - - if (SIZ(r) < 0) - abort (); /* only r>=0 */ - - limb = bit / GMP_LIMB_BITS; - if (SIZ(r) <= limb) - return 0; - - bit %= GMP_LIMB_BITS; - return (PTR(r)[limb] >> bit) & 1; -} - -int -popc_limb (mp_limb_t a) -{ - int ret = 0; - while (a != 0) - { - ret += (a & 1); - a >>= 1; - } - return ret; -} - -unsigned long -mpz_popcount (mpz_t a) -{ - unsigned long ret; - int i; - - if (SIZ(a) < 0) - abort (); - - ret = 0; - for (i = 0; i < SIZ(a); i++) - ret += popc_limb (PTR(a)[i]); - return ret; -} - -void -mpz_add (mpz_t r, mpz_t a, mpz_t b) -{ - int an = ABSIZ (a), bn = ABSIZ (b), rn; - mp_limb_t *rp, *ap, *bp; - int i; - mp_limb_t t, cy; - - if ((SIZ (a) ^ SIZ (b)) < 0) - abort (); /* really subtraction */ - if (SIZ (a) < 0) - abort (); - - mpz_realloc (r, MAX (an, bn) + 1); - ap = PTR (a); bp = PTR (b); rp = PTR (r); - if (an < bn) - { - mp_limb_t *tp; int tn; - tn = an; an = bn; bn = tn; - tp = ap; ap = bp; bp = tp; - } - - cy = 0; - for (i = 0; i < bn; i++) - { - t = ap[i] + bp[i] + cy; - rp[i] = LO (t); - cy = HI (t); - } - for (i = bn; i < an; i++) - { - t = ap[i] + cy; - rp[i] = LO (t); - cy = HI (t); - } - rp[an] = cy; - rn = an + 1; - - mpn_normalize (rp, &rn); - SIZ (r) = rn; -} - -void -mpz_add_ui (mpz_t r, mpz_t a, unsigned long int ui) -{ - mpz_t b; - - mpz_init (b); - mpz_set_ui (b, ui); - mpz_add (r, a, b); - mpz_clear (b); -} - -void -mpz_sub (mpz_t r, mpz_t a, mpz_t b) -{ - int an = ABSIZ (a), bn = ABSIZ (b), rn; - mp_limb_t *rp, *ap, *bp; - int i; - mp_limb_t t, cy; - - if ((SIZ (a) ^ SIZ (b)) < 0) - abort (); /* really addition */ - if (SIZ (a) < 0) - abort (); - - mpz_realloc (r, MAX (an, bn) + 1); - ap = PTR (a); bp = PTR (b); rp = PTR (r); - if (an < bn) - { - mp_limb_t *tp; int tn; - tn = an; an = bn; bn = tn; - tp = ap; ap = bp; bp = tp; - /* This needs sign change, not done so abort. */ - abort (); - } - - cy = 0; - for (i = 0; i < bn; i++) - { - t = ap[i] - bp[i] - cy; - rp[i] = LO (t); - cy = LO (-HI (t)); - } - for (i = bn; i < an; i++) - { - t = ap[i] - cy; - rp[i] = LO (t); - cy = LO (-HI (t)); - } - rp[an] = cy; - rn = an + 1; - - if (cy != 0) - { - cy = 0; - for (i = 0; i < rn; i++) - { - t = -rp[i] - cy; - rp[i] = LO (t); - cy = LO (-HI (t)); - } - SIZ (r) = -rn; - return; - } - - mpn_normalize (rp, &rn); - SIZ (r) = rn; -} - -void -mpz_sub_ui (mpz_t r, mpz_t a, unsigned long int ui) -{ - mpz_t b; - - mpz_init (b); - mpz_set_ui (b, ui); - mpz_sub (r, a, b); - mpz_clear (b); -} - -void -mpz_mul (mpz_t r, mpz_t a, mpz_t b) -{ - int an = ABSIZ (a), bn = ABSIZ (b), rn; - mp_limb_t *scratch, *tmp, *ap = PTR (a), *bp = PTR (b); - int ai, bi; - mp_limb_t t, cy; - - scratch = xmalloc_limbs (an + bn); - tmp = scratch; - - for (bi = 0; bi < bn; bi++) - tmp[bi] = 0; - - for (ai = 0; ai < an; ai++) - { - tmp = scratch + ai; - cy = 0; - for (bi = 0; bi < bn; bi++) - { - t = ap[ai] * bp[bi] + tmp[bi] + cy; - tmp[bi] = LO (t); - cy = HI (t); - } - tmp[bn] = cy; - } - - rn = an + bn; - mpn_normalize (scratch, &rn); - free (PTR (r)); - PTR (r) = scratch; - SIZ (r) = (SIZ (a) ^ SIZ (b)) >= 0 ? rn : -rn; - ALLOC (r) = an + bn; -} - -void -mpz_mul_ui (mpz_t r, mpz_t a, unsigned long int ui) -{ - mpz_t b; - - mpz_init (b); - mpz_set_ui (b, ui); - mpz_mul (r, a, b); - mpz_clear (b); -} - -void -mpz_mul_2exp (mpz_t r, mpz_t a, unsigned long int bcnt) -{ - mpz_set (r, a); - -#if 0 - { - unsigned long int lcnt; - - lcnt = bcnt / GMP_LIMB_BITS; - if (lcnt > 0) - { - int rn = ABSIZ (r); - mp_limb_t *rp; - int i; - - mpz_realloc (r, rn + lcnt); - rp = PTR (r); - - for (i = rn - 1; i >= 0; i--) - rp[i + lcnt] = rp[i]; - for (i = lcnt - 1; i >= 0; i--) - rp[i] = 0; - - rn += lcnt; - SIZ (r) = SIZ (r) >= 0 ? rn : -rn; - - bcnt %= GMP_LIMB_BITS; - } - } -#endif - - while (bcnt) - { - mpz_add (r, r, r); - bcnt -= 1; - } -} - -void -mpz_ui_pow_ui (mpz_t r, unsigned long b, unsigned long e) -{ - unsigned long i; - mpz_t bz; - - mpz_init (bz); - mpz_set_ui (bz, b); - - mpz_set_ui (r, 1L); - for (i = 0; i < e; i++) - mpz_mul (r, r, bz); - - mpz_clear (bz); -} - -void -mpz_tdiv_q_2exp (mpz_t r, mpz_t a, unsigned long int bcnt) -{ - int as, rn; - int cnt, tnc; - int lcnt; - mp_limb_t high_limb, low_limb; - int i; - - as = SIZ (a); - lcnt = bcnt / GMP_LIMB_BITS; - rn = ABS (as) - lcnt; - if (rn <= 0) - SIZ (r) = 0; - else - { - mp_limb_t *rp, *ap; - - mpz_realloc (r, rn); - - rp = PTR (r); - ap = PTR (a); - - cnt = bcnt % GMP_LIMB_BITS; - if (cnt != 0) - { - ap += lcnt; - tnc = GMP_LIMB_BITS - cnt; - high_limb = *ap++; - low_limb = high_limb >> cnt; - - for (i = rn - 1; i != 0; i--) - { - high_limb = *ap++; - *rp++ = low_limb | LO (high_limb << tnc); - low_limb = high_limb >> cnt; - } - *rp = low_limb; - rn -= low_limb == 0; - } - else - { - ap += lcnt; - mpn_copyi (rp, ap, rn); - } - - SIZ (r) = as >= 0 ? rn : -rn; - } -} - -void -mpz_tdiv_r_2exp (mpz_t r, mpz_t a, unsigned long int bcnt) -{ - int rn, bwhole; - - mpz_set (r, a); - rn = ABSIZ(r); - - bwhole = bcnt / GMP_LIMB_BITS; - bcnt %= GMP_LIMB_BITS; - if (rn > bwhole) - { - rn = bwhole+1; - PTR(r)[rn-1] &= ((mp_limb_t) 1 << bcnt) - 1; - mpn_normalize (PTR(r), &rn); - SIZ(r) = (SIZ(r) >= 0 ? rn : -rn); - } -} - -int -mpz_cmp (mpz_t a, mpz_t b) -{ - mp_limb_t *ap, *bp, al, bl; - int as = SIZ (a), bs = SIZ (b); - int i; - int sign; - - if (as != bs) - return as > bs ? 1 : -1; - - sign = as > 0 ? 1 : -1; - - ap = PTR (a); - bp = PTR (b); - for (i = ABS (as) - 1; i >= 0; i--) - { - al = ap[i]; - bl = bp[i]; - if (al != bl) - return al > bl ? sign : -sign; - } - return 0; -} - -int -mpz_cmp_ui (mpz_t a, unsigned long b) -{ - mpz_t bz; - int ret; - mpz_init_set_ui (bz, b); - ret = mpz_cmp (a, bz); - mpz_clear (bz); - return ret; -} - -void -mpz_tdiv_qr (mpz_t q, mpz_t r, mpz_t a, mpz_t b) -{ - mpz_t tmpr, tmpb; - unsigned long cnt; - - ASSERT (mpz_sgn(a) >= 0); - ASSERT (mpz_sgn(b) > 0); - - mpz_init_set (tmpr, a); - mpz_init_set (tmpb, b); - mpz_set_ui (q, 0L); - - if (mpz_cmp (tmpr, tmpb) >= 0) - { - cnt = mpz_sizeinbase (tmpr, 2) - mpz_sizeinbase (tmpb, 2) + 1; - mpz_mul_2exp (tmpb, tmpb, cnt); - - for ( ; cnt > 0; cnt--) - { - mpz_mul_2exp (q, q, 1); - mpz_tdiv_q_2exp (tmpb, tmpb, 1L); - if (mpz_cmp (tmpr, tmpb) >= 0) - { - mpz_sub (tmpr, tmpr, tmpb); - mpz_add_ui (q, q, 1L); - ASSERT (mpz_cmp (tmpr, tmpb) < 0); - } - } - } - - mpz_set (r, tmpr); - mpz_clear (tmpr); - mpz_clear (tmpb); -} - -void -mpz_tdiv_qr_ui (mpz_t q, mpz_t r, mpz_t a, unsigned long b) -{ - mpz_t bz; - mpz_init_set_ui (bz, b); - mpz_tdiv_qr (q, r, a, bz); - mpz_clear (bz); -} - -void -mpz_tdiv_q (mpz_t q, mpz_t a, mpz_t b) -{ - mpz_t r; - - mpz_init (r); - mpz_tdiv_qr (q, r, a, b); - mpz_clear (r); -} - -void -mpz_tdiv_r (mpz_t r, mpz_t a, mpz_t b) -{ - mpz_t q; - - mpz_init (q); - mpz_tdiv_qr (q, r, a, b); - mpz_clear (q); -} - -void -mpz_tdiv_q_ui (mpz_t q, mpz_t n, unsigned long d) -{ - mpz_t dz; - mpz_init_set_ui (dz, d); - mpz_tdiv_q (q, n, dz); - mpz_clear (dz); -} - -/* Set inv to the inverse of d, in the style of invert_limb, ie. for - udiv_qrnnd_preinv. */ -void -mpz_preinv_invert (mpz_t inv, mpz_t d, int numb_bits) -{ - mpz_t t; - int norm; - ASSERT (SIZ(d) > 0); - - norm = numb_bits - mpz_sizeinbase (d, 2); - ASSERT (norm >= 0); - mpz_init_set_ui (t, 1L); - mpz_mul_2exp (t, t, 2*numb_bits - norm); - mpz_tdiv_q (inv, t, d); - mpz_set_ui (t, 1L); - mpz_mul_2exp (t, t, numb_bits); - mpz_sub (inv, inv, t); - - mpz_clear (t); -} - -/* Remove leading '0' characters from the start of a string, by copying the - remainder down. */ -void -strstrip_leading_zeros (char *s) -{ - char c, *p; - - p = s; - while (*s == '0') - s++; - - do - { - c = *s++; - *p++ = c; - } - while (c != '\0'); -} - -char * -mpz_get_str (char *buf, int base, mpz_t a) -{ - static char tohex[] = "0123456789abcdef"; - - mp_limb_t alimb, *ap; - int an, bn, i, j; - char *bp; - - if (base != 16) - abort (); - if (SIZ (a) < 0) - abort (); - - if (buf == 0) - buf = xmalloc (ABSIZ (a) * (GMP_LIMB_BITS / 4) + 3); - - an = ABSIZ (a); - if (an == 0) - { - buf[0] = '0'; - buf[1] = '\0'; - return buf; - } - - ap = PTR (a); - bn = an * (GMP_LIMB_BITS / 4); - bp = buf + bn; - - for (i = 0; i < an; i++) - { - alimb = ap[i]; - for (j = 0; j < GMP_LIMB_BITS / 4; j++) - { - bp--; - *bp = tohex [alimb & 0xF]; - alimb >>= 4; - } - ASSERT (alimb == 0); - } - ASSERT (bp == buf); - - buf[bn] = '\0'; - - strstrip_leading_zeros (buf); - return buf; -} - -void -mpz_out_str (FILE *file, int base, mpz_t a) -{ - char *str; - - if (file == 0) - file = stdout; - - str = mpz_get_str (0, 16, a); - fputs (str, file); - free (str); -} - -/* Calculate r satisfying r*d == 1 mod 2^n. */ -void -mpz_invert_2exp (mpz_t r, mpz_t a, unsigned long n) -{ - unsigned long i; - mpz_t inv, prod; - - ASSERT (mpz_odd_p (a)); - - mpz_init_set_ui (inv, 1L); - mpz_init (prod); - - for (i = 1; i < n; i++) - { - mpz_mul (prod, inv, a); - if (mpz_tstbit (prod, i) != 0) - mpz_setbit (inv, i); - } - - mpz_mul (prod, inv, a); - mpz_tdiv_r_2exp (prod, prod, n); - ASSERT (mpz_cmp_ui (prod, 1L) == 0); - - mpz_set (r, inv); - - mpz_clear (inv); - mpz_clear (prod); -} - -/* Calculate inv satisfying r*a == 1 mod 2^n. */ -void -mpz_invert_ui_2exp (mpz_t r, unsigned long a, unsigned long n) -{ - mpz_t az; - mpz_init_set_ui (az, a); - mpz_invert_2exp (r, az, n); - mpz_clear (az); -} - -/* x=y^z */ -void -mpz_pow_ui (mpz_t x, mpz_t y, unsigned long z) -{ - mpz_t t; - - mpz_init_set_ui (t, 1); - for (; z != 0; z--) - mpz_mul (t, t, y); - mpz_set (x, t); - mpz_clear (t); -} - -/* x=x+y*z */ -void -mpz_addmul_ui (mpz_t x, mpz_t y, unsigned long z) -{ - mpz_t t; - - mpz_init (t); - mpz_mul_ui (t, y, z); - mpz_add (x, x, t); - mpz_clear (t); -} - -/* x=floor(y^(1/z)) */ -void -mpz_root (mpz_t x, mpz_t y, unsigned long z) -{ - mpz_t t, u; - - if (mpz_sgn (y) < 0) - { - fprintf (stderr, "mpz_root does not accept negative values\n"); - abort (); - } - if (mpz_cmp_ui (y, 1) <= 0) - { - mpz_set (x, y); - return; - } - mpz_init (t); - - /* One-bit initial approximation */ - mpz_init_set_ui (u, 1); - mpz_mul_2exp (u, u, ((mpz_sizeinbase (y, 2) - 1) / z) + 1); - - for (;;) - { - mpz_pow_ui (t, u, z - 1); - mpz_tdiv_q (t, y, t); - mpz_addmul_ui (t, u, z - 1); - mpz_tdiv_q_ui (t, t, z); - if (mpz_cmp (t, u) >= 0) - break; - mpz_set (u, t); - } - - mpz_set (x, u); - mpz_clear (t); - mpz_clear (u); -} - -/* x=floor(y^(1/2)) */ -void -mpz_sqrt (mpz_t x, mpz_t y) -{ - mpz_t t, u; - - if (mpz_sgn (y) < 0) - { - fprintf (stderr, "mpz_sqrt does not accept negative values\n"); - abort (); - } - if (mpz_cmp_ui (y, 1) <= 0) - { - mpz_set (x, y); - return; - } - mpz_init (t); - - /* One-bit initial approximation */ - mpz_init_set_ui (u, 1); - mpz_mul_2exp (u, u, ((mpz_sizeinbase (y, 2) - 1) / 2) + 1); - - for (;;) - { - mpz_tdiv_q (t, y, u); - mpz_add (t, t, u); - mpz_tdiv_q_2exp (t, t, 1); - if (mpz_cmp (t, u) >= 0) - break; - mpz_set (u, t); - } - - mpz_set (x, u); - mpz_clear (t); - mpz_clear (u); -} diff -r b641c8181188 gen-bases.c --- a/gen-bases.c Wed Feb 15 14:46:40 2012 +0100 +++ b/gen-bases.c Fri Feb 17 11:16:53 2012 +0100 @@ -20,7 +20,7 @@ along with the GNU MP Library. If not, #include -#include "dumbmp.c" +#include "bootstrap.c" int chars_per_limb; diff -r b641c8181188 gen-fac_ui.c --- a/gen-fac_ui.c Wed Feb 15 14:46:40 2012 +0100 +++ b/gen-fac_ui.c Fri Feb 17 11:16:53 2012 +0100 @@ -20,7 +20,7 @@ along with the GNU MP Library. If not, #include #include -#include "dumbmp.c" +#include "bootstrap.c" /* sets x=y*(y+2)*(y+4)*....*(y+2*(z-1)) */ diff -r b641c8181188 gen-fib.c --- a/gen-fib.c Wed Feb 15 14:46:40 2012 +0100 +++ b/gen-fib.c Fri Feb 17 11:16:53 2012 +0100 @@ -18,7 +18,7 @@ You should have received a copy of the G along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ #include -#include "dumbmp.c" +#include "bootstrap.c" mpz_t *f; int fnum, fib_limit, luc_limit; @@ -34,7 +34,7 @@ generate (int numb_bits) /* fib(2n) > 2^n, so use 2n as a limit for the table size */ falloc = 2 * numb_bits; - f = (mpz_t *) xmalloc (falloc * sizeof (*f)); + f = xmalloc (falloc * sizeof (*f)); mpz_init_set_ui (f[0], 1L); /* F[-1] */ mpz_init_set_ui (f[1], 0L); /* F[0] */ @@ -43,7 +43,7 @@ generate (int numb_bits) for (i = 2; ; i++) { - ASSERT (i < falloc); + assert (i < falloc); /* F[i] = F[i-1] + F[i-2] */ mpz_init (f[i]); diff -r b641c8181188 gen-psqr.c --- a/gen-psqr.c Wed Feb 15 14:46:40 2012 +0100 +++ b/gen-psqr.c Fri Feb 17 11:16:53 2012 +0100 @@ -20,7 +20,7 @@ along with the GNU MP Library. If not, #include #include -#include "dumbmp.c" +#include "bootstrap.c" /* The aim of this program is to choose either mpn_mod_34lsub1 or mpn_mod_1 @@ -152,9 +152,9 @@ f_cmp_fraction (const void *parg, const accordingly. */ #define COLLAPSE_ELEMENT(array, idx, narray) \ do { \ - mem_copyi ((char *) &(array)[idx], \ - (char *) &(array)[idx+1], \ - ((narray)-((idx)+1)) * sizeof (array[0])); \ + memmove (&(array)[idx], \ + &(array)[idx+1], \ + ((narray)-((idx)+1)) * sizeof (array[0])); \ (narray)--; \ } while (0) @@ -173,7 +173,7 @@ mul_2exp_mod (int n, int p, int m) int neg_mod (int n, int m) { - ASSERT (n >= 0 && n < m); + assert (n >= 0 && n < m); return (n == 0 ? 0 : m-n); } @@ -202,7 +202,7 @@ generate_sq_res_0x100 (int limb_bits) int i, res; nsq_res_0x100 = (0x100 + limb_bits - 1) / limb_bits; - sq_res_0x100 = (mpz_t *) xmalloc (nsq_res_0x100 * sizeof (*sq_res_0x100)); + sq_res_0x100 = xmalloc (nsq_res_0x100 * sizeof (*sq_res_0x100)); for (i = 0; i < nsq_res_0x100; i++) mpz_init_set_ui (sq_res_0x100[i], 0L); @@ -233,9 +233,8 @@ generate_mod (int limb_bits, int nail_bi /* no more than limb_bits many factors in a one limb modulus (and of course in reality nothing like that many) */ factor_alloc = limb_bits; - factor = (struct factor_t *) xmalloc (factor_alloc * sizeof (*factor)); - rawfactor = (struct rawfactor_t *) - xmalloc (factor_alloc * sizeof (*rawfactor)); + factor = xmalloc (factor_alloc * sizeof (*factor)); + rawfactor = xmalloc (factor_alloc * sizeof (*rawfactor)); if (numb_bits % 4 != 0) { @@ -301,7 +300,7 @@ generate_mod (int limb_bits, int nail_bi } while (mpz_sgn (r) == 0); - ASSERT (nrawfactor < factor_alloc); + assert (nrawfactor < factor_alloc); rawfactor[nrawfactor].divisor = i; rawfactor[nrawfactor].multiplicity = multiplicity; nrawfactor++; @@ -341,7 +340,7 @@ generate_mod (int limb_bits, int nail_bi break; mpz_set (pp, new_pp); - ASSERT (nrawfactor < factor_alloc); + assert (nrawfactor < factor_alloc); rawfactor[nrawfactor].divisor = i; rawfactor[nrawfactor].multiplicity = 1; nrawfactor++; @@ -377,7 +376,7 @@ generate_mod (int limb_bits, int nail_bi for (i = 0; i < nrawfactor; i++) { int j; - ASSERT (nfactor < factor_alloc); + assert (nfactor < factor_alloc); factor[nfactor].divisor = 1; for (j = 0; j < rawfactor[i].multiplicity; j++) factor[nfactor].divisor *= rawfactor[i].divisor; diff -r b641c8181188 gen-trialdivtab.c --- a/gen-trialdivtab.c Wed Feb 15 14:46:40 2012 +0100 +++ b/gen-trialdivtab.c Fri Feb 17 11:16:53 2012 +0100 @@ -36,7 +36,7 @@ along with the GNU MP Library. If not, #include #include -#include "dumbmp.c" +#include "bootstrap.c" int sumspills (mpz_t, mpz_t *, int); void mpn_mod_1s_4p_cps (mpz_t [7], mpz_t); diff -r b641c8181188 mini-gmp/README --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mini-gmp/README Fri Feb 17 11:16:53 2012 +0100 @@ -0,0 +1,46 @@ +This is "mini-gmp", a small implementation of a subset of GMP's mpn +and mpz interfaces. + +It is intended for applications which need arithmetic on numbers +larger than a machine word, but which don't need to handle very large +numbers very efficiently. Those applications can include a copy of +mini-gmp to get a GMP-compatible interface with small footprint. One +can also arrange for optional linking with the real GMP library, using +mini-gmp as a fallback when for some reason GMP is not available, or +not desired as a dependency. + +The supported GMP subset is declared in mini-gmp.h. The implemented +functions are fully compatible with the corresponding GMP functions, +as specified in the GMP manual, with a few exceptions: + + mpz_export and mpz_import support only NAILS = 0. + + The REALLOC_FUNC and FREE_FUNC registered with + mp_set_memory_functions does not get the correct size of the + allocated block in the corresponding argument. mini-gmp always + passes zero for these rarely used arguments. + +The implementation is a single file, mini-gmp.c. + +The performance target for mini-gmp is to be at most 10 times slower +than the real GMP library, for numbers of size up to a few hundred +bits. No asymptotically fast algorithms are included in mini-gmp, so +it will be many orders of magnitude slower than GMP for very large +numbers. + +You should never "install" mini-gmp. Applications can either just +#include mini-gmp.c (but then, beware that it defines several macros +and functions outside of the advertised interface). Or compile +mini-gmp.c as a separate compilation unit, and use the declarations in +mini-gmp.h. + +The tests subdirectory contains a testsuite. To use it, you need GMP +and GNU make. Just run make check in the tests directory. If the +hard-coded compiler settings are not right, you have to either edit the +Makefile or pass overriding values on the make command line (e.g., +make CC=cc check). Testing is not (yet) as thorough as for the real +GMP. + +The current version was put together by Niels Möller +, with a fair amount of copy-and-paste from the +GMP sources. diff -r b641c8181188 mini-gmp/mini-gmp.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mini-gmp/mini-gmp.c Fri Feb 17 11:16:53 2012 +0100 @@ -0,0 +1,4041 @@ +/* mini-gmp, a minimalistic implementation of a GNU GMP subset. + + Contributed to the GNU project by Niels Möller + +Copyright 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1999, 2000, 2001, +2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012 Free +Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 3 of the License, or (at your +option) any later version. + +The GNU MP Library is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public +License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ + +/* NOTE: All functions in this file which are not declared in + mini-gmp.h are internal, and are not intended to be compatible + neither with GMP nor with future versions of mini-gmp. */ + +/* Much of the material copied from GMP files, including: gmp-impl.h, + longlong.h, mpn/generic/add_n.c, mpn/generic/addmul_1.c, + mpn/generic/lshift.c, mpn/generic/mul_1.c, + mpn/generic/mul_basecase.c, mpn/generic/rshift.c, + mpn/generic/sbpi1_div_qr.c, mpn/generic/sub_n.c, + mpn/generic/submul_1.c. */ + +#include +#include +#include +#include +#include +#include + +#include "mini-gmp.h" + + +/* Macros */ +#define GMP_LIMB_BITS (sizeof(mp_limb_t) * CHAR_BIT) + +#define GMP_LIMB_MAX (~ (mp_limb_t) 0) +#define GMP_LIMB_HIGHBIT ((mp_limb_t) 1 << (GMP_LIMB_BITS - 1)) + +#define GMP_HLIMB_BIT ((mp_limb_t) 1 << (GMP_LIMB_BITS / 2)) +#define GMP_LLIMB_MASK (GMP_HLIMB_BIT - 1) + +#define GMP_ULONG_BITS (sizeof(unsigned long) * CHAR_BIT) +#define GMP_ULONG_HIGHBIT ((unsigned long) 1 << (GMP_ULONG_BITS - 1)) + +#define GMP_ABS(x) ((x) >= 0 ? (x) : -(x)) + +#define GMP_MIN(a, b) ((a) < (b) ? (a) : (b)) +#define GMP_MAX(a, b) ((a) > (b) ? (a) : (b)) + +#define gmp_assert_nocarry(x) do { \ + mp_limb_t __cy = x; \ + assert (__cy == 0); \ + } while (0) + +#define gmp_clz(count, x) do { \ + mp_limb_t __clz_x = (x); \ + unsigned __clz_c; \ + for (__clz_c = 0; \ + (__clz_x & ((mp_limb_t) 0xff << (GMP_LIMB_BITS - 8))) == 0; \ + __clz_c += 8) \ + __clz_x <<= 8; \ + for (; (__clz_x & GMP_LIMB_HIGHBIT) == 0; __clz_c++) \ + __clz_x <<= 1; \ + (count) = __clz_c; \ + } while (0) + +#define gmp_ctz(count, x) do { \ + mp_limb_t __ctz_x = (x); \ + unsigned __ctz_c = 0; \ + gmp_clz (__ctz_c, __ctz_x & - __ctz_x); \ + (count) = GMP_LIMB_BITS - 1 - __ctz_c; \ + } while (0) + +#define gmp_add_ssaaaa(sh, sl, ah, al, bh, bl) \ + do { \ + mp_limb_t __x; \ + __x = (al) + (bl); \ + (sh) = (ah) + (bh) + (__x < (al)); \ + (sl) = __x; \ + } while (0) + +#define gmp_sub_ddmmss(sh, sl, ah, al, bh, bl) \ + do { \ + mp_limb_t __x; \ + __x = (al) - (bl); \ + (sh) = (ah) - (bh) - ((al) < (bl)); \ + (sl) = __x; \ + } while (0) + +#define gmp_umul_ppmm(w1, w0, u, v) \ + do { \ + mp_limb_t __x0, __x1, __x2, __x3; \ + unsigned __ul, __vl, __uh, __vh; \ + mp_limb_t __u = (u), __v = (v); \ + \ + __ul = __u & GMP_LLIMB_MASK; \ + __uh = __u >> (GMP_LIMB_BITS / 2); \ + __vl = __v & GMP_LLIMB_MASK; \ + __vh = __v >> (GMP_LIMB_BITS / 2); \ + \ + __x0 = (mp_limb_t) __ul * __vl; \ + __x1 = (mp_limb_t) __ul * __vh; \ + __x2 = (mp_limb_t) __uh * __vl; \ + __x3 = (mp_limb_t) __uh * __vh; \ + \ + __x1 += __x0 >> (GMP_LIMB_BITS / 2);/* this can't give carry */ \ + __x1 += __x2; /* but this indeed can */ \ + if (__x1 < __x2) /* did we get it? */ \ + __x3 += GMP_HLIMB_BIT; /* yes, add it in the proper pos. */ \ + \ + (w1) = __x3 + (__x1 >> (GMP_LIMB_BITS / 2)); \ + (w0) = (__x1 << (GMP_LIMB_BITS / 2)) + (__x0 & GMP_LLIMB_MASK); \ + } while (0) + +#define gmp_udiv_qrnnd_preinv(q, r, nh, nl, d, di) \ + do { \ + mp_limb_t _qh, _ql, _r, _mask; \ + gmp_umul_ppmm (_qh, _ql, (nh), (di)); \ + gmp_add_ssaaaa (_qh, _ql, _qh, _ql, (nh) + 1, (nl)); \ + _r = (nl) - _qh * (d); \ + _mask = -(mp_limb_t) (_r > _ql); /* both > and >= are OK */ \ + _qh += _mask; \ + _r += _mask & (d); \ + if (_r >= (d)) \ + { \ + _r -= (d); \ + _qh++; \ + } \ + \ + (r) = _r; \ + (q) = _qh; \ + } while (0) + +#define gmp_udiv_qr_3by2(q, r1, r0, n2, n1, n0, d1, d0, dinv) \ + do { \ + mp_limb_t _q0, _t1, _t0, _mask; \ + gmp_umul_ppmm ((q), _q0, (n2), (dinv)); \ + gmp_add_ssaaaa ((q), _q0, (q), _q0, (n2), (n1)); \ + \ + /* Compute the two most significant limbs of n - q'd */ \ + (r1) = (n1) - (d1) * (q); \ + gmp_sub_ddmmss ((r1), (r0), (r1), (n0), (d1), (d0)); \ + gmp_umul_ppmm (_t1, _t0, (d0), (q)); \ + gmp_sub_ddmmss ((r1), (r0), (r1), (r0), _t1, _t0); \ + (q)++; \ + \ + /* Conditionally adjust q and the remainders */ \ + _mask = - (mp_limb_t) ((r1) >= _q0); \ + (q) += _mask; \ + gmp_add_ssaaaa ((r1), (r0), (r1), (r0), _mask & (d1), _mask & (d0)); \ + if ((r1) >= (d1)) \ + { \ + if ((r1) > (d1) || (r0) >= (d0)) \ + { \ + (q)++; \ + gmp_sub_ddmmss ((r1), (r0), (r1), (r0), (d1), (d0)); \ + } \ + } \ + } while (0) + +/* Swap macros. */ +#define MP_LIMB_T_SWAP(x, y) \ + do { \ + mp_limb_t __mp_limb_t_swap__tmp = (x); \ + (x) = (y); \ + (y) = __mp_limb_t_swap__tmp; \ + } while (0) +#define MP_SIZE_T_SWAP(x, y) \ + do { \ + mp_size_t __mp_size_t_swap__tmp = (x); \ + (x) = (y); \ + (y) = __mp_size_t_swap__tmp; \ + } while (0) +#define MP_BITCNT_T_SWAP(x,y) \ + do { \ + mp_bitcnt_t __mp_bitcnt_t_swap__tmp = (x); \ + (x) = (y); \ + (y) = __mp_bitcnt_t_swap__tmp; \ + } while (0) +#define MP_PTR_SWAP(x, y) \ + do { \ + mp_ptr __mp_ptr_swap__tmp = (x); \ + (x) = (y); \ + (y) = __mp_ptr_swap__tmp; \ + } while (0) +#define MP_SRCPTR_SWAP(x, y) \ + do { \ + mp_srcptr __mp_srcptr_swap__tmp = (x); \ + (x) = (y); \ + (y) = __mp_srcptr_swap__tmp; \ + } while (0) + +#define MPN_PTR_SWAP(xp,xs, yp,ys) \ + do { \ + MP_PTR_SWAP (xp, yp); \ + MP_SIZE_T_SWAP (xs, ys); \ + } while(0) +#define MPN_SRCPTR_SWAP(xp,xs, yp,ys) \ + do { \ + MP_SRCPTR_SWAP (xp, yp); \ + MP_SIZE_T_SWAP (xs, ys); \ + } while(0) + +#define MPZ_PTR_SWAP(x, y) \ + do { \ + mpz_ptr __mpz_ptr_swap__tmp = (x); \ + (x) = (y); \ + (y) = __mpz_ptr_swap__tmp; \ + } while (0) +#define MPZ_SRCPTR_SWAP(x, y) \ + do { \ + mpz_srcptr __mpz_srcptr_swap__tmp = (x); \ + (x) = (y); \ + (y) = __mpz_srcptr_swap__tmp; \ + } while (0) + + +/* Memory allocation and other helper functions. */ +static void +gmp_die (const char *msg) +{ + fprintf (stderr, "%s\n", msg); + abort(); +} + +static void * +gmp_default_alloc (size_t size) +{ + void *p; + + assert (size > 0); + + p = malloc (size); + if (!p) + gmp_die("gmp_default_alloc: Virtual memory exhausted."); + + return p; +} + +static void * +gmp_default_realloc (void *old, size_t old_size, size_t new_size) +{ + mp_ptr p; + + p = realloc (old, new_size); + + if (!p) + gmp_die("gmp_default_realoc: Virtual memory exhausted."); + + return p; +} + +static void +gmp_default_free (void *p, size_t size) +{ + free (p); +} + +static void * (*gmp_allocate_func) (size_t) = gmp_default_alloc; +static void * (*gmp_reallocate_func) (void *, size_t, size_t) = gmp_default_realloc; +static void (*gmp_free_func) (void *, size_t) = gmp_default_free; + +void +mp_get_memory_functions (void *(**alloc_func) (size_t), + void *(**realloc_func) (void *, size_t, size_t), + void (**free_func) (void *, size_t)) +{ + if (alloc_func) + *alloc_func = gmp_allocate_func; + + if (realloc_func) + *realloc_func = gmp_reallocate_func; + + if (free_func) + *free_func = gmp_free_func; +} + +void +mp_set_memory_functions (void *(*alloc_func) (size_t), + void *(*realloc_func) (void *, size_t, size_t), + void (*free_func) (void *, size_t)) +{ + if (!alloc_func) + alloc_func = gmp_default_alloc; + if (!realloc_func) + realloc_func = gmp_default_realloc; + if (!free_func) + free_func = gmp_default_free; + + gmp_allocate_func = alloc_func; + gmp_reallocate_func = realloc_func; + gmp_free_func = free_func; +} + +#define gmp_xalloc(size) ((*gmp_allocate_func)((size))) +#define gmp_free(p) ((*gmp_free_func) ((p), 0)) + +static mp_ptr +gmp_xalloc_limbs (mp_size_t size) +{ + return gmp_xalloc (size * sizeof (mp_limb_t)); +} + +static mp_ptr +gmp_xrealloc_limbs (mp_ptr old, mp_size_t size) +{ + assert (size > 0); + return (*gmp_reallocate_func) (old, 0, size * sizeof (mp_limb_t)); +} + + +/* MPN interface */ + +void +mpn_copyi (mp_ptr d, mp_srcptr s, mp_size_t n) +{ + mp_size_t i; + for (i = 0; i < n; i++) + d[i] = s[i]; +} + +void +mpn_copyd (mp_ptr d, mp_srcptr s, mp_size_t n) +{ + while (n-- > 0) + d[n] = s[n]; +} + +int +mpn_cmp (mp_srcptr ap, mp_srcptr bp, mp_size_t n) +{ + for (; n > 0; n--) + { + if (ap[n-1] < bp[n-1]) + return -1; + else if (ap[n-1] > bp[n-1]) + return 1; + } + return 0; +} + +static int +mpn_cmp4 (mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn) +{ + if (an > bn) + return 1; + else if (an < bn) + return -1; + else + return mpn_cmp (ap, bp, an); +} + +static mp_size_t +mpn_normalized_size (mp_srcptr xp, mp_size_t n) +{ + for (; n > 0 && xp[n-1] == 0; n--) + ; + return n; +} + +#define mpn_zero_p(xp, n) (mpn_normalized_size ((xp), (n)) == 0) + +mp_limb_t +mpn_add_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b) +{ + mp_size_t i; + + assert (n > 0); + + for (i = 0; i < n; i++) + { + mp_limb_t r = ap[i] + b; + /* Carry out */ + b = (r < b); + rp[i] = r; + } + return b; +} + +mp_limb_t +mpn_add_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n) +{ + mp_size_t i; + mp_limb_t cy; + + for (i = 0, cy = 0; i < n; i++) + { + mp_limb_t a, b, r; + a = ap[i]; b = bp[i]; + r = a + cy; + cy = (r < cy); + r += b; + cy += (r < b); + rp[i] = r; + } + return cy; +} + +mp_limb_t +mpn_add (mp_ptr rp, mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn) +{ + mp_limb_t cy; + + assert (an >= bn); + + cy = mpn_add_n (rp, ap, bp, bn); + if (an > bn) + cy = mpn_add_1 (rp + bn, ap + bn, an - bn, cy); + return cy; +} + +mp_limb_t +mpn_sub_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b) +{ + mp_size_t i; + + assert (n > 0); + + for (i = 0; i < n; i++) + { + mp_limb_t a = ap[i]; + /* Carry out */ + mp_limb_t cy = a < b;; + rp[i] = a - b; + b = cy; + } + return b; +} + +mp_limb_t +mpn_sub_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n) +{ + mp_size_t i; + mp_limb_t cy; + + for (i = 0, cy = 0; i < n; i++) + { + mp_limb_t a, b; + a = ap[i]; b = bp[i]; + b += cy; + cy = (b < cy); + cy += (a < b); + rp[i] = a - b; + } + return cy; +} + +mp_limb_t +mpn_sub (mp_ptr rp, mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn) +{ + mp_limb_t cy; + + assert (an >= bn); + + cy = mpn_sub_n (rp, ap, bp, bn); + if (an > bn) + cy = mpn_sub_1 (rp + bn, ap + bn, an - bn, cy); + return cy; +} + +mp_limb_t +mpn_mul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl) +{ + mp_limb_t ul, cl, hpl, lpl; + + assert (n >= 1); + + cl = 0; + do + { + ul = *up++; + gmp_umul_ppmm (hpl, lpl, ul, vl); + + lpl += cl; + cl = (lpl < cl) + hpl; + + *rp++ = lpl; + } + while (--n != 0); + + return cl; +} + +mp_limb_t +mpn_addmul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl) +{ + mp_limb_t ul, cl, hpl, lpl, rl; + + assert (n >= 1); + + cl = 0; + do + { + ul = *up++; + gmp_umul_ppmm (hpl, lpl, ul, vl); + + lpl += cl; + cl = (lpl < cl) + hpl; + + rl = *rp; + lpl = rl + lpl; + cl += lpl < rl; + *rp++ = lpl; + } + while (--n != 0); + + return cl; +} + +mp_limb_t +mpn_submul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl) +{ + mp_limb_t ul, cl, hpl, lpl, rl; + + assert (n >= 1); + + cl = 0; + do + { + ul = *up++; + gmp_umul_ppmm (hpl, lpl, ul, vl); + + lpl += cl; + cl = (lpl < cl) + hpl; + + rl = *rp; + lpl = rl - lpl; + cl += lpl > rl; + *rp++ = lpl; + } + while (--n != 0); + + return cl; +} + +mp_limb_t +mpn_mul (mp_ptr rp, mp_srcptr up, mp_size_t un, mp_srcptr vp, mp_size_t vn) +{ + assert (un >= vn); + assert (vn >= 1); + + /* We first multiply by the low order limb. This result can be + stored, not added, to rp. We also avoid a loop for zeroing this + way. */ + + rp[un] = mpn_mul_1 (rp, up, un, vp[0]); + rp += 1, vp += 1, vn -= 1; + + /* Now accumulate the product of up[] and the next higher limb from + vp[]. */ + + while (vn >= 1) + { + rp[un] = mpn_addmul_1 (rp, up, un, vp[0]); + rp += 1, vp += 1, vn -= 1; + } + return rp[un - 1]; +} + +void +mpn_mul_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n) +{ + mpn_mul (rp, ap, n, bp, n); +} + +void +mpn_sqr (mp_ptr rp, mp_srcptr ap, mp_size_t n) +{ + mpn_mul (rp, ap, n, ap, n); +} + +mp_limb_t +mpn_lshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt) +{ + mp_limb_t high_limb, low_limb; + unsigned int tnc; + mp_size_t i; + mp_limb_t retval; + + assert (n >= 1); + assert (cnt >= 1); + assert (cnt < GMP_LIMB_BITS); + + up += n; + rp += n; + + tnc = GMP_LIMB_BITS - cnt; + low_limb = *--up; + retval = low_limb >> tnc; + high_limb = (low_limb << cnt); + + for (i = n - 1; i != 0; i--) + { + low_limb = *--up; + *--rp = high_limb | (low_limb >> tnc); + high_limb = (low_limb << cnt); + } + *--rp = high_limb; + + return retval; +} + +mp_limb_t +mpn_rshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt) +{ + mp_limb_t high_limb, low_limb; + unsigned int tnc; + mp_size_t i; + mp_limb_t retval; + + assert (n >= 1); + assert (cnt >= 1); + assert (cnt < GMP_LIMB_BITS); + + tnc = GMP_LIMB_BITS - cnt; + high_limb = *up++; + retval = (high_limb << tnc); + low_limb = high_limb >> cnt; + + for (i = n - 1; i != 0; i--) + { + high_limb = *up++; + *rp++ = low_limb | (high_limb << tnc); + low_limb = high_limb >> cnt; + } + *rp = low_limb; + + return retval; +} + + +/* MPN division interface. */ +mp_limb_t +mpn_invert_3by2 (mp_limb_t u1, mp_limb_t u0) +{ + mp_limb_t r, p, m; + unsigned ul, uh; + unsigned ql, qh; + + /* First, do a 2/1 inverse. */ + /* The inverse m is defined as floor( (B^2 - 1 - u1)/u1 ), so that 0 < + * B^2 - (B + m) u1 <= u1 */ + assert (u1 >= GMP_LIMB_HIGHBIT); + + ul = u1 & GMP_LLIMB_MASK; + uh = u1 >> (GMP_LIMB_BITS / 2); + + qh = ~u1 / uh; + r = ((~u1 - (mp_limb_t) qh * uh) << (GMP_LIMB_BITS / 2)) | GMP_LLIMB_MASK; + + p = (mp_limb_t) qh * ul; + /* Adjustment steps taken from udiv_qrnnd_c */ + if (r < p) + { + qh--; + r += u1; + if (r >= u1) /* i.e. we didn't get carry when adding to r */ + if (r < p) + { + qh--; + r += u1; + } + } + r -= p; + + /* Do a 3/2 division (with half limb size) */ + p = (r >> (GMP_LIMB_BITS / 2)) * qh + r; + ql = (p >> (GMP_LIMB_BITS / 2)) + 1; + + /* By the 3/2 method, we don't need the high half limb. */ + r = (r << (GMP_LIMB_BITS / 2)) + GMP_LLIMB_MASK - ql * u1; + + if (r >= (p << (GMP_LIMB_BITS / 2))) + { + ql--; + r += u1; + } + m = ((mp_limb_t) qh << (GMP_LIMB_BITS / 2)) + ql; + if (r >= u1) + { + m++; + r -= u1; + } + + if (u0 > 0) + { + mp_limb_t th, tl; + r = ~r; + r += u0; + if (r < u0) + { + m--; + if (r >= u1) + { + m--; + r -= u1; + } + r -= u1; + } + gmp_umul_ppmm (th, tl, u0, m); + r += th; + if (r < th) + { + m--; + if (r > u1 || (r == u1 && tl > u0)) + m--; + } + } + + return m; +} + +struct gmp_div_inverse +{ + /* Normalization shift count. */ + unsigned shift; + /* Normalized divisor (d0 unused for mpn_div_qr_1) */ + mp_limb_t d1, d0; + /* Inverse, for 2/1 or 3/2. */ + mp_limb_t di; +}; + +static void +mpn_div_qr_1_invert (struct gmp_div_inverse *inv, mp_limb_t d) +{ + unsigned shift; + + assert (d > 0); + gmp_clz (shift, d); + inv->shift = shift; + inv->d1 = d << shift; + inv->di = mpn_invert_limb (inv->d1); +} + +static void +mpn_div_qr_2_invert (struct gmp_div_inverse *inv, + mp_limb_t d1, mp_limb_t d0) +{ + unsigned shift; + + assert (d1 > 0); + gmp_clz (shift, d1); + inv->shift = shift; + if (shift > 0) + { + d1 = (d1 << shift) | (d0 >> (GMP_LIMB_BITS - shift)); + d0 <<= shift; + } + inv->d1 = d1; + inv->d0 = d0; + inv->di = mpn_invert_3by2 (d1, d0); +} + +static void +mpn_div_qr_invert (struct gmp_div_inverse *inv, + mp_srcptr dp, mp_size_t dn) +{ + assert (dn > 0); + + if (dn == 1) + mpn_div_qr_1_invert (inv, dp[0]); + else if (dn == 2) + mpn_div_qr_2_invert (inv, dp[1], dp[0]); + else + { + unsigned shift; + mp_limb_t d1, d0; + + d1 = dp[dn-1]; + d0 = dp[dn-2]; + assert (d1 > 0); + gmp_clz (shift, d1); + inv->shift = shift; + if (shift > 0) + { + d1 = (d1 << shift) | (d0 >> (GMP_LIMB_BITS - shift)); + d0 = (d0 << shift) | (dp[dn-3] >> (GMP_LIMB_BITS - shift)); + } + inv->d1 = d1; + inv->d0 = d0; + inv->di = mpn_invert_3by2 (d1, d0); + } +} + +/* Not matching current public gmp interface, rather corresponding to + the sbpi1_div_* functions. */ +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) +{ + mp_limb_t d, di; + mp_limb_t r; + mp_ptr tp = NULL; + + if (inv->shift > 0) + { + tp = gmp_xalloc_limbs (nn); + r = mpn_lshift (tp, np, nn, inv->shift); + np = tp; + } + else + r = 0; + + d = inv->d1; + di = inv->di; + while (nn-- > 0) + { + mp_limb_t q; + + gmp_udiv_qrnnd_preinv (q, r, r, np[nn], d, di); + if (qp) + qp[nn] = q; + } + if (inv->shift > 0) + gmp_free (tp); + + return r >> inv->shift; +} + +static mp_limb_t +mpn_div_qr_1 (mp_ptr qp, mp_srcptr np, mp_size_t nn, mp_limb_t d) +{ + assert (d > 0); + + /* Special case for powers of two. */ + if (d > 1 && (d & (d-1)) == 0) + { + unsigned shift; + mp_limb_t r = np[0] & (d-1); + gmp_ctz (shift, d); + if (qp) + mpn_rshift (qp, np, nn, shift); + + return r; + } + else + { + struct gmp_div_inverse inv; + mpn_div_qr_1_invert (&inv, d); + return mpn_div_qr_1_preinv (qp, np, nn, &inv); + } +} + +static void +mpn_div_qr_2_preinv (mp_ptr qp, mp_ptr rp, mp_srcptr np, mp_size_t nn, + const struct gmp_div_inverse *inv) +{ + unsigned shift; + mp_size_t i; + mp_limb_t d1, d0, di, r1, r0; + mp_ptr tp; + + assert (nn >= 2); + shift = inv->shift; + d1 = inv->d1; + d0 = inv->d0; + di = inv->di; + + if (shift > 0) + { + tp = gmp_xalloc_limbs (nn); + r1 = mpn_lshift (tp, np, nn, shift); + np = tp; + } + else + r1 = 0; + + r0 = np[nn - 1]; + + for (i = nn - 2; i >= 0; i--) + { + mp_limb_t n0, q; + n0 = np[i]; + gmp_udiv_qr_3by2 (q, r1, r0, r1, r0, n0, d1, d0, di); + + if (qp) + qp[i] = q; + } + + if (shift > 0) + { + assert ((r0 << (GMP_LIMB_BITS - shift)) == 0); + r0 = (r0 >> shift) | (r1 << (GMP_LIMB_BITS - shift)); + r1 >>= shift; + + gmp_free (tp); + } + + rp[1] = r1; + rp[0] = r0; +} + +#if 0 +static void +mpn_div_qr_2 (mp_ptr qp, mp_ptr rp, mp_srcptr np, mp_size_t nn, + mp_limb_t d1, mp_limb_t d0) +{ + struct gmp_div_inverse inv; + assert (nn >= 2); + + mpn_div_qr_2_invert (&inv, d1, d0); + mpn_div_qr_2_preinv (qp, rp, np, nn, &inv); +} +#endif + +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) +{ + mp_size_t i; + + mp_limb_t d1, d0; + mp_limb_t cy, cy1; + mp_limb_t q; + + assert (dn > 2); + assert (nn >= dn); + assert ((dp[dn-1] & GMP_LIMB_HIGHBIT) != 0); + + d1 = dp[dn - 1]; + d0 = dp[dn - 2]; + + /* Iteration variable is the index of the q limb. + * + * We divide + * by + */ + + for (i = nn - dn; i >= 0; i--) + { + mp_limb_t n0 = np[dn-1+i]; + + if (n1 == d1 && n0 == d0) + { + q = GMP_LIMB_MAX; + mpn_submul_1 (np+i, dp, dn, q); + n1 = np[dn-1+i]; /* update n1, last loop's value will now be invalid */ + } + else + { + gmp_udiv_qr_3by2 (q, n1, n0, n1, n0, np[dn-2+i], d1, d0, dinv); + + cy = mpn_submul_1 (np + i, dp, dn-2, q); + + cy1 = n0 < cy; + n0 = n0 - cy; + cy = n1 < cy1; + n1 = n1 - cy1; + np[dn-2+i] = n0; + + if (cy != 0) + { + n1 += d1 + mpn_add_n (np + i, np + i, dp, dn - 1); + q--; + } + } + + if (qp) + qp[i] = q; + } + + np[dn - 1] = n1; +} + +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) +{ + assert (dn > 0); + assert (nn >= dn); + + if (dn == 1) + np[0] = mpn_div_qr_1_preinv (qp, np, nn, inv); + else if (dn == 2) + mpn_div_qr_2_preinv (qp, np, np, nn, inv); + else + { + mp_limb_t nh; + unsigned shift; + + assert (dp[dn-1] & GMP_LIMB_HIGHBIT); + + shift = inv->shift; + if (shift > 0) + nh = mpn_lshift (np, np, nn, shift); + else + nh = 0; + + assert (inv->d1 == dp[dn-1]); + assert (inv->d0 == dp[dn-2]); + + mpn_div_qr_pi1 (qp, np, nn, nh, dp, dn, inv->di); + + if (shift > 0) + gmp_assert_nocarry (mpn_rshift (np, np, dn, shift)); + } +} + +static void +mpn_div_qr (mp_ptr qp, mp_ptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn) +{ + struct gmp_div_inverse inv; + mp_ptr tp = NULL; + + assert (dn > 0); + assert (nn >= dn); + + mpn_div_qr_invert (&inv, dp, dn); + if (dn > 2 && inv.shift > 0) + { + tp = gmp_xalloc_limbs (dn); + gmp_assert_nocarry (mpn_lshift (tp, dp, dn, inv.shift)); + dp = tp; + } + mpn_div_qr_preinv (qp, np, nn, dp, dn, &inv); + if (tp) + gmp_free (tp); +} + + +/* MPN base conversion. */ +static unsigned +mpn_base_power_of_two_p (unsigned b) +{ + switch (b) + { + case 2: return 1; + case 4: return 2; + case 8: return 3; + case 16: return 4; + case 32: return 5; + case 64: return 6; + case 128: return 7; + case 256: return 8; + default: return 0; + } +} + +struct mpn_base_info +{ + /* bb is the largest power of the base which fits in one limb, and + exp is the corresponding exponent. */ + unsigned exp; + mp_limb_t bb; +}; + +static void +mpn_get_base_info (struct mpn_base_info *info, mp_limb_t b) +{ + mp_limb_t m; + mp_limb_t p; + unsigned exp; + + m = GMP_LIMB_MAX / b; + for (exp = 1, p = b; p <= m; exp++) + p *= b; + + info->exp = exp; + info->bb = p; +} + +static mp_bitcnt_t +mpn_limb_size_in_base_2 (mp_limb_t u) +{ + unsigned shift; + + assert (u > 0); + gmp_clz (shift, u); + return GMP_LIMB_BITS - shift; +} + +static size_t +mpn_get_str_bits (unsigned char *sp, unsigned bits, mp_srcptr up, mp_size_t un) +{ + unsigned char mask; + size_t sn, j; + mp_size_t i; + int shift; + + sn = ((un - 1) * GMP_LIMB_BITS + mpn_limb_size_in_base_2 (up[un-1]) + + bits - 1) / bits; + + mask = (1U << bits) - 1; + + for (i = 0, j = sn, shift = 0; j-- > 0;) + { + unsigned char digit = up[i] >> shift; + + shift += bits; + + if (shift >= GMP_LIMB_BITS && ++i < un) + { + shift -= GMP_LIMB_BITS; + digit |= up[i] << (bits - shift); + } + sp[j] = digit & mask; + } + return sn; +} + +/* We generate digits from the least significant end, and reverse at + the end. */ +static size_t +mpn_limb_get_str (unsigned char *sp, mp_limb_t w, + const struct gmp_div_inverse *binv) +{ + mp_size_t i; + for (i = 0; w > 0; i++) + { + mp_limb_t h, l, r; + + h = w >> (GMP_LIMB_BITS - binv->shift); + l = w << binv->shift; + + gmp_udiv_qrnnd_preinv (w, r, h, l, binv->d1, binv->di); + assert ( (r << (GMP_LIMB_BITS - binv->shift)) == 0); + r >>= binv->shift; + + sp[i] = r; + } + return i; +} + +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) +{ + struct gmp_div_inverse binv; + size_t sn; + size_t i; + + mpn_div_qr_1_invert (&binv, base); + + sn = 0; + + if (un > 1) + { + struct gmp_div_inverse bbinv; + mpn_div_qr_1_invert (&bbinv, info->bb); + + do + { + mp_limb_t w; + size_t done; + w = mpn_div_qr_1_preinv (up, up, un, &bbinv); + un -= (up[un-1] == 0); + done = mpn_limb_get_str (sp + sn, w, &binv); + + for (sn += done; done < info->exp; done++) + sp[sn++] = 0; + } + while (un > 1); + } + sn += mpn_limb_get_str (sp + sn, up[0], &binv); + + /* Reverse order */ + for (i = 0; 2*i + 1 < sn; i++) + { + unsigned char t = sp[i]; + sp[i] = sp[sn - i - 1]; + sp[sn - i - 1] = t; + } + + return sn; +} + +size_t +mpn_get_str (unsigned char *sp, int base, mp_ptr up, mp_size_t un) +{ + unsigned bits; + + assert (un > 0); + assert (up[un-1] > 0); + + bits = mpn_base_power_of_two_p (base); + if (bits) + return mpn_get_str_bits (sp, bits, up, un); + else + { + struct mpn_base_info info; + + mpn_get_base_info (&info, base); + return mpn_get_str_other (sp, base, &info, up, un); + } +} + +static mp_size_t +mpn_set_str_bits (mp_ptr rp, const unsigned char *sp, size_t sn, + unsigned bits) +{ + mp_size_t rn; + size_t j; + unsigned shift; + + for (j = sn, rn = 0, shift = 0; j-- > 0; ) + { + if (shift == 0) + { + rp[rn++] = sp[j]; + shift += bits; + } + else + { + rp[rn-1] |= (mp_limb_t) sp[j] << shift; + shift += bits; + if (shift >= GMP_LIMB_BITS) + { + shift -= GMP_LIMB_BITS; + if (shift > 0) + rp[rn++] = (mp_limb_t) sp[j] >> (bits - shift); + } + } + } + rn = mpn_normalized_size (rp, rn); + return rn; +} + +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) +{ + mp_size_t rn; + mp_limb_t w; + unsigned first; + unsigned k; + size_t j; + + first = 1 + (sn - 1) % info->exp; + + j = 0; + w = sp[j++]; + for (k = 1; k < first; k++) + w = w * b + sp[j++]; + + rp[0] = w; + + for (rn = (w > 0); j < sn;) + { + mp_limb_t cy; + + w = sp[j++]; + for (k = 1; k < info->exp; k++) + w = w * b + sp[j++]; + + cy = mpn_mul_1 (rp, rp, rn, info->bb); + cy += mpn_add_1 (rp, rp, rn, w); + if (cy > 0) + rp[rn++] = cy; + } + assert (j == sn); + + return rn; +} + +mp_size_t +mpn_set_str (mp_ptr rp, const unsigned char *sp, size_t sn, int base) +{ + unsigned bits; + + if (sn == 0) + return 0; + + bits = mpn_base_power_of_two_p (base); + if (bits) + return mpn_set_str_bits (rp, sp, sn, bits); + else + { + struct mpn_base_info info; + + mpn_get_base_info (&info, base); + return mpn_set_str_other (rp, sp, sn, base, &info); + } +} + + +/* MPZ interface */ +void +mpz_init (mpz_t r) +{ + r->_mp_alloc = 1; + r->_mp_size = 0; + r->_mp_d = gmp_xalloc_limbs (1); +} + +/* The utility of this function is a bit limited, since many functions + assings the result variable using mpz_swap. */ +void +mpz_init2 (mpz_t r, mp_bitcnt_t bits) +{ + mp_size_t rn; + + bits -= (bits != 0); /* Round down, except if 0 */ + rn = 1 + bits / GMP_LIMB_BITS; + + r->_mp_alloc = rn; + r->_mp_size = 0; + r->_mp_d = gmp_xalloc_limbs (rn); +} + +void +mpz_clear (mpz_t r) +{ + gmp_free (r->_mp_d); +} + +static void * +mpz_realloc (mpz_t r, mp_size_t size) +{ + if (size < 1) + size = 1; + + r->_mp_d = gmp_xrealloc_limbs (r->_mp_d, size); + r->_mp_alloc = size; + + if (GMP_ABS (r->_mp_size) > size) + r->_mp_size = 0; + + return r->_mp_d; +} + +/* Realloc for an mpz_t WHAT if it has less than NEEDED limbs. */ +#define MPZ_REALLOC(z,n) ((n) > (z)->_mp_alloc \ + ? mpz_realloc(z,n) \ + : (z)->_mp_d) + +/* MPZ assignment and basic conversions. */ +void +mpz_set_si (mpz_t r, signed long int x) +{ + if (x > 0) + { + r->_mp_size = 1; + r->_mp_d[0] = x; + } + else if (x < 0) + { + r->_mp_size = -1; + r->_mp_d[0] = -x; + } + else + r->_mp_size = 0; +} + +void +mpz_set_ui (mpz_t r, unsigned long int x) +{ + if (x > 0) + { + r->_mp_size = 1; + r->_mp_d[0] = x; + } + else + r->_mp_size = 0; +} + +void +mpz_set (mpz_t r, const mpz_t x) +{ + /* Allow the NOP r == x */ + if (r != x) + { + mp_size_t n; + mp_ptr rp; + + n = GMP_ABS (x->_mp_size); + rp = MPZ_REALLOC (r, n); + + mpn_copyi (rp, x->_mp_d, n); + r->_mp_size = x->_mp_size; + } +} + +void +mpz_init_set_si (mpz_t r, signed long int x) +{ + mpz_init (r); + mpz_set_si (r, x); +} + +void +mpz_init_set_ui (mpz_t r, unsigned long int x) +{ + mpz_init (r); + mpz_set_ui (r, x); +} + +void +mpz_init_set (mpz_t r, const mpz_t x) +{ + mpz_init (r); + mpz_set (r, x); +} + +int +mpz_fits_slong_p (const mpz_t u) +{ + mp_size_t us = u->_mp_size; + + if (us == 0) + return 1; + else if (us == 1) + return u->_mp_d[0] < (GMP_LIMB_HIGHBIT / 2); + else if (us == -1) + return u->_mp_d[0] <= (GMP_LIMB_HIGHBIT / 2); + else + return 0; +} + +int +mpz_fits_ulong_p (const mpz_t u) +{ + mp_size_t us = u->_mp_size; + + return us == 0 || us == 1; +} + +long int +mpz_get_si (const mpz_t u) +{ + mp_size_t us = u->_mp_size; + + if (us > 0) + return (long) (u->_mp_d[0] & ~GMP_LIMB_HIGHBIT); + else if (us < 0) + return (long) (- u->_mp_d[0] | GMP_LIMB_HIGHBIT); + else + return 0; +} + +unsigned long int +mpz_get_ui (const mpz_t u) +{ + return u->_mp_size == 0 ? 0 : u->_mp_d[0]; +} + +size_t +mpz_size (const mpz_t u) +{ + return GMP_ABS (u->_mp_size); +} + +mp_limb_t +mpz_getlimbn (const mpz_t u, mp_size_t n) +{ + if (n >= 0 && n < GMP_ABS (u->_mp_size)) + return u->_mp_d[n]; + else + return 0; +} + + +/* Conversions and comparison to double. */ +void +mpz_set_d (mpz_t r, double x) +{ + int sign; + mp_ptr rp; + mp_size_t rn, i; + double B; + double Bi; + mp_limb_t f; + + /* x != x is true when x is a NaN, and x == x * 0.5 is true when x is + zero or infinity. */ + if (x == 0.0 || x != x || x == x * 0.5) + { + r->_mp_size = 0; + return; + } + + if (x < 0.0) + { + x = - x; + sign = 1; + } + else + sign = 0; + + if (x < 1.0) + { + r->_mp_size = 0; + return; + } + B = 2.0 * (double) GMP_LIMB_HIGHBIT; + Bi = 1.0 / B; + for (rn = 1; x >= B; rn++) + x *= Bi; + + rp = MPZ_REALLOC (r, rn); + + f = (mp_limb_t) x; + x -= f; + assert (x < 1.0); + rp[rn-1] = f; + for (i = rn-1; i-- > 0; ) + { + x = B * x; + f = (mp_limb_t) x; + x -= f; + assert (x < 1.0); + rp[i] = f; + } + + r->_mp_size = sign ? - rn : rn; +} + +void +mpz_init_set_d (mpz_t r, double x) +{ + mpz_init (r); + mpz_set_d (r, x); +} + +double +mpz_get_d (const mpz_t u) +{ + mp_size_t un; + double x; + double B = 2.0 * (double) GMP_LIMB_HIGHBIT; + + un = GMP_ABS (u->_mp_size); + + if (un == 0) + return 0.0; + + x = u->_mp_d[--un]; + while (un > 0) + x = B*x + u->_mp_d[--un]; + + if (u->_mp_size < 0) + x = -x; + + return x; +} + +int +mpz_cmp_d (const mpz_t x, double d) +{ + mp_size_t xn = x->_mp_size; + int sign; + double B, Bi; + mp_size_t i; + + if (xn == 0) + { + if (d < 0) + return 1; + else if (d > 0) + return -1; + else + return 0; + } + + if (xn < 0) + { + xn = -xn; + d = -d; + sign = -1; + } + else + sign = 1; + + if (d < 1.0) + return sign; + + B = 2.0 * (double) GMP_LIMB_HIGHBIT; + Bi = 1.0 / B; + for (i = 1; i < xn; i++) + { + d *= Bi; + if (d < 1.0) + return sign; + } + if (d >= B) + return -sign; + + for (i = xn; i-- > 0; i++) + { + mp_limb_t f, xl; + + f = (mp_limb_t) d; + xl = x->_mp_d[i]; + if (xl > f) + return sign; + else if (xl < f) + return -sign; + d = B * (d - f); + } + + if (d > 0) + return -sign; + else + return 0; +} + + +/* MPN comparisons and the like. */ +int +mpz_sgn (const mpz_t u) +{ + mp_size_t usize = u->_mp_size; + + if (usize > 0) + return 1; + else if (usize < 0) + return -1; + else + return 0; +} + +int +mpz_cmp_si (const mpz_t u, long v) +{ + mp_size_t usize = u->_mp_size; + + if (usize > 1) + return 1; + else if (usize < -1) + return -1; + else if (usize == 0) + { + if (v > 0) + return -1; + else if (v < 0) + return 1; + } + else + { + mp_limb_t ul = u->_mp_d[0]; + if (usize == 1) + { + if (v <= 0 || (mp_limb_t) v < ul) + return 1; + else if ((mp_limb_t) v > ul) + return -1; + } + else /* usize == -1 */ + { + if (v >= 0 || (mp_limb_t) -v < ul) + return -1; + else if ((mp_limb_t) -v > ul) + return 1; + } + } + return 0; +} + +int +mpz_cmp_ui (const mpz_t u, unsigned long v) +{ + mp_size_t usize = u->_mp_size; + + if (usize > 1) + return 1; + else if (usize < 0) + return -1; + else + { + mp_limb_t ul = (usize > 0) ? u->_mp_d[0] : 0; + if (ul > v) + return 1; + else if (ul < v) + return -1; + } + return 0; +} + +int +mpz_cmp (const mpz_t a, const mpz_t b) +{ + mp_size_t asize = a->_mp_size; + mp_size_t bsize = b->_mp_size; + + if (asize > bsize) + return 1; + else if (asize < bsize) + return -1; + else if (asize > 0) + return mpn_cmp (a->_mp_d, b->_mp_d, asize); + else if (asize < 0) + return -mpn_cmp (a->_mp_d, b->_mp_d, asize); + else + return 0; +} + +int +mpz_cmpabs_ui (const mpz_t u, unsigned long v) +{ + mp_size_t un = GMP_ABS (u->_mp_size); + mp_limb_t ul; + + if (un > 1) + return 1; + + ul = (un == 1) ? u->_mp_d[0] : 0; + + if (ul > v) + return 1; + else if (ul < v) + return -1; + else + return 0; +} + +int +mpz_cmpabs (const mpz_t u, const mpz_t v) +{ + return mpn_cmp4 (u->_mp_d, GMP_ABS (u->_mp_size), + v->_mp_d, GMP_ABS (v->_mp_size)); +} + +void +mpz_abs (mpz_t r, const mpz_t u) +{ + if (r != u) + mpz_set (r, u); + + r->_mp_size = GMP_ABS (r->_mp_size); +} + +void +mpz_neg (mpz_t r, const mpz_t u) +{ + if (r != u) + mpz_set (r, u); + + r->_mp_size = -r->_mp_size; +} + +void +mpz_swap (mpz_t u, mpz_t v) +{ + MP_SIZE_T_SWAP (u->_mp_size, v->_mp_size); + MP_SIZE_T_SWAP (u->_mp_alloc, v->_mp_alloc); + MP_PTR_SWAP (u->_mp_d, v->_mp_d); +} + + +/* MPZ addition and subtraction */ + +/* Adds to the absolute value. Returns new size, but doesn't store it. */ +static mp_size_t +mpz_abs_add_ui (mpz_t r, const mpz_t a, unsigned long b) +{ + mp_size_t an; + mp_ptr rp; + mp_limb_t cy; + + an = GMP_ABS (a->_mp_size); + if (an == 0) + { + r->_mp_d[0] = b; + return b > 0; + } + + rp = MPZ_REALLOC (r, an + 1); + + cy = mpn_add_1 (rp, a->_mp_d, an, b); + rp[an] = cy; + an += (cy > 0); + + return an; +} + +/* Subtract from the absolute value. Returns new size, (or -1 on underflow), + but doesn't store it. */ +static mp_size_t +mpz_abs_sub_ui (mpz_t r, const mpz_t a, unsigned long b) +{ + mp_size_t an = GMP_ABS (a->_mp_size); + mp_ptr rp = MPZ_REALLOC (r, an); + + if (an == 0) + { + rp[0] = b; + return -(b > 0); + } + else if (an == 1 && a->_mp_d[0] < b) + { + rp[0] = b - a->_mp_d[0]; + return -1; + } + else + { + gmp_assert_nocarry (mpn_sub_1 (rp, a->_mp_d, an, b)); + return mpn_normalized_size (rp, an); + } +} + +void +mpz_add_ui (mpz_t r, const mpz_t a, unsigned long b) +{ + if (a->_mp_size >= 0) + r->_mp_size = mpz_abs_add_ui (r, a, b); + else + r->_mp_size = -mpz_abs_sub_ui (r, a, b); +} + +void +mpz_sub_ui (mpz_t r, const mpz_t a, unsigned long b) +{ + if (a->_mp_size < 0) + r->_mp_size = -mpz_abs_add_ui (r, a, b); + else + r->_mp_size = mpz_abs_sub_ui (r, a, b); +} + +void +mpz_ui_sub (mpz_t r, unsigned long a, const mpz_t b) +{ + if (b->_mp_size < 0) + r->_mp_size = mpz_abs_add_ui (r, b, a); + else + r->_mp_size = -mpz_abs_sub_ui (r, b, a); +} + +static mp_size_t +mpz_abs_add (mpz_t r, const mpz_t a, const mpz_t b) +{ + mp_size_t an = GMP_ABS (a->_mp_size); + mp_size_t bn = GMP_ABS (b->_mp_size); + mp_ptr ap = a->_mp_d; + mp_ptr bp = b->_mp_d; + mp_ptr rp; + mp_limb_t cy; + + if (an < bn) + MPN_PTR_SWAP (ap, an, bp, bn); + + rp = MPZ_REALLOC (r, an + 1); + + cy = mpn_add (rp, ap, an, bp, bn); + rp[an] = cy; + an += (cy > 0); + + return an; +} + +static mp_size_t +mpz_abs_sub (mpz_t r, const mpz_t a, const mpz_t b) +{ + mp_size_t an = GMP_ABS (a->_mp_size); + mp_size_t bn = GMP_ABS (b->_mp_size); + int cmp; + mp_ptr rp; + + cmp = mpn_cmp4 (a->_mp_d, an, b->_mp_d, bn); + if (cmp > 0) + { + rp = MPZ_REALLOC (r, an); + gmp_assert_nocarry (mpn_sub (rp, a->_mp_d, an, b->_mp_d, bn)); + return mpn_normalized_size (rp, an); + } + else if (cmp < 0) + { + rp = MPZ_REALLOC (r, bn); + gmp_assert_nocarry (mpn_sub (rp, b->_mp_d, bn, a->_mp_d, an)); + return -mpn_normalized_size (rp, bn); + } + else + return 0; +} + +void +mpz_add (mpz_t r, const mpz_t a, const mpz_t b) +{ + mp_size_t rn; + + if ( (a->_mp_size ^ b->_mp_size) >= 0) + rn = mpz_abs_add (r, a, b); + else + rn = mpz_abs_sub (r, a, b); + + r->_mp_size = a->_mp_size >= 0 ? rn : - rn; +} + +void +mpz_sub (mpz_t r, const mpz_t a, const mpz_t b) +{ + mp_size_t rn; + + if ( (a->_mp_size ^ b->_mp_size) >= 0) + rn = mpz_abs_sub (r, a, b); + else + rn = mpz_abs_add (r, a, b); + + r->_mp_size = a->_mp_size >= 0 ? rn : - rn; +} + + +/* MPZ multiplication */ +void +mpz_mul_si (mpz_t r, const mpz_t u, long int v) +{ + int sign; + mp_size_t un; + mpz_t t; + mp_ptr tp; + mp_limb_t cy; + + if (v == 0) + { + r->_mp_size = 0; + return; + } + sign = (u->_mp_size ^ v) < 0; + + un = GMP_ABS (u->_mp_size); + mpz_init (t); + + tp = MPZ_REALLOC (t, un + 1); + cy = mpn_mul_1 (tp, u->_mp_d, un, GMP_ABS (v)); + tp[un] = cy; + un += cy > 0; + + t->_mp_size = sign ? - un : un; + mpz_swap (r, t); + mpz_clear (t); +} + +void +mpz_mul_ui (mpz_t r, const mpz_t u, unsigned long int v) +{ + mp_size_t un; + mpz_t t; + mp_ptr tp; + mp_limb_t cy; + + un = GMP_ABS (u->_mp_size); + + if (un == 0 || v == 0) + { + r->_mp_size = 0; + return; + } + + mpz_init (t); + + tp = MPZ_REALLOC (t, un + 1); + cy = mpn_mul_1 (tp, u->_mp_d, un, v); + tp[un] = cy; + + t->_mp_size = un + (cy > 0); + + mpz_swap (r, t); + mpz_clear (t); +} + +void +mpz_mul (mpz_t r, const mpz_t u, const mpz_t v) +{ + int sign; + mp_size_t un, vn, rn; + mpz_t t; + mp_ptr tp; + + un = GMP_ABS (u->_mp_size); + vn = GMP_ABS (v->_mp_size); + + if (un == 0 || vn == 0) + { + r->_mp_size = 0; + return; + } + + sign = (u->_mp_size ^ v->_mp_size) < 0; + + mpz_init (t); + + tp = MPZ_REALLOC (t, un + vn); + if (un >= vn) + mpn_mul (tp, u->_mp_d, un, v->_mp_d, vn); + else + mpn_mul (tp, v->_mp_d, vn, u->_mp_d, un); + + rn = un + vn; + rn -= tp[rn-1] == 0; + + t->_mp_size = sign ? - rn : rn; + mpz_swap (r, t); + mpz_clear (t); +} + +void +mpz_mul_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bits) +{ + mp_size_t un, rn; + mp_size_t limbs; + unsigned shift; + mp_ptr rp; + + un = GMP_ABS (u->_mp_size); + if (un == 0) + { + r->_mp_size = 0; + return; + } + + limbs = bits / GMP_LIMB_BITS; + shift = bits % GMP_LIMB_BITS; + + rn = un + limbs + (shift > 0); + rp = MPZ_REALLOC (r, rn); + if (shift > 0) + { + mp_limb_t cy = mpn_lshift (rp + limbs, u->_mp_d, un, shift); + rp[rn-1] = cy; + rn -= (cy == 0); + } + else + mpn_copyd (rp + limbs, u->_mp_d, un); + + while (limbs > 0) + rp[--limbs] = 0; + + r->_mp_size = (u->_mp_size < 0) ? - rn : rn; +} + + +/* MPZ division */ +enum mpz_div_round_mode { DIV_FLOOR, DIV_CEIL, DIV_TRUNC }; + +/* Allows q or r to be zero. Returns 1 iff remainder is non-zero. */ +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) +{ + mp_size_t ns, ds, nn, dn, qs; + ns = n->_mp_size; + ds = d->_mp_size; + + if (ds == 0) + gmp_die("mpz_div_qr: Divide by zero."); + + if (ns == 0) + { + if (q) + q->_mp_size = 0; + if (r) + r->_mp_size = 0; + return 0; + } + + nn = GMP_ABS (ns); + dn = GMP_ABS (ds); + + qs = ds ^ ns; + + if (nn < dn) + { + if (mode == DIV_CEIL && qs >= 0) + { + /* q = 1, r = n - d */ + if (r) + mpz_sub (r, n, d); + if (q) + mpz_set_ui (q, 1); + } + else if (mode == DIV_FLOOR && qs < 0) + { + /* q = -1, r = n + d */ + if (r) + mpz_add (r, n, d); + if (q) + mpz_set_si (q, -1); + } + else + { + /* q = 0, r = d */ + if (r) + mpz_set (r, n); + if (q) + q->_mp_size = 0; + } + return 1; + } + else + { + mp_ptr np, qp; + mp_size_t qn, rn; + mpz_t tq, tr; + + mpz_init (tr); + mpz_set (tr, n); + np = tr->_mp_d; + + qn = nn - dn + 1; + + if (q) + { + mpz_init (tq); + qp = MPZ_REALLOC (tq, qn); + } + else + qp = NULL; + + mpn_div_qr (qp, np, nn, d->_mp_d, dn); + + if (qp) + { + qn -= (qp[qn-1] == 0); + + tq->_mp_size = qs < 0 ? -qn : qn; + } + rn = mpn_normalized_size (np, dn); + tr->_mp_size = ns < 0 ? - rn : rn; + + if (mode == DIV_FLOOR && qs < 0 && rn != 0) + { + if (q) + mpz_sub_ui (tq, tq, 1); + if (r) + mpz_add (tr, tr, d); + } + else if (mode == DIV_CEIL && qs >= 0 && rn != 0) + { + if (q) + mpz_add_ui (tq, tq, 1); + if (r) + mpz_sub (tr, tr, d); + } + + if (q) + { + mpz_swap (tq, q); + mpz_clear (tq); + } + if (r) + mpz_swap (tr, r); + + mpz_clear (tr); + + return rn != 0; + } +} + +void +mpz_cdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d) +{ + mpz_div_qr (q, r, n, d, DIV_CEIL); +} + +void +mpz_fdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d) +{ + mpz_div_qr (q, r, n, d, DIV_FLOOR); +} + +void +mpz_tdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d) +{ + mpz_div_qr (q, r, n, d, DIV_TRUNC); +} + +void +mpz_cdiv_q (mpz_t q, const mpz_t n, const mpz_t d) +{ + mpz_div_qr (q, NULL, n, d, DIV_CEIL); +} + +void +mpz_fdiv_q (mpz_t q, const mpz_t n, const mpz_t d) +{ + mpz_div_qr (q, NULL, n, d, DIV_FLOOR); +} + +void +mpz_tdiv_q (mpz_t q, const mpz_t n, const mpz_t d) +{ + mpz_div_qr (q, NULL, n, d, DIV_TRUNC); +} + +void +mpz_cdiv_r (mpz_t r, const mpz_t n, const mpz_t d) +{ + mpz_div_qr (NULL, r, n, d, DIV_CEIL); +} + +void +mpz_fdiv_r (mpz_t r, const mpz_t n, const mpz_t d) +{ + mpz_div_qr (NULL, r, n, d, DIV_FLOOR); +} + +void +mpz_tdiv_r (mpz_t r, const mpz_t n, const mpz_t d) +{ + mpz_div_qr (NULL, r, n, d, DIV_TRUNC); +} + +static void +mpz_div_q_2exp (mpz_t q, const mpz_t u, mp_bitcnt_t bit_index, + enum mpz_div_round_mode mode) +{ + mp_size_t un, qn; + mp_size_t limb_cnt; + mp_ptr qp; + mp_limb_t adjust; + + un = u->_mp_size; + if (un == 0) + { + q->_mp_size = 0; + return; + } + limb_cnt = bit_index / GMP_LIMB_BITS; + qn = GMP_ABS (un) - limb_cnt; + bit_index %= GMP_LIMB_BITS; + + adjust = 0; + if ( (mode == DIV_FLOOR && un < 0) || (mode == DIV_CEIL && un > 0)) + /* Note: Below, the final indexing at limb_cnt is valid because at + that point we have qn > 0. */ + adjust = (!mpn_zero_p (u->_mp_d, limb_cnt) + || qn <= 0 + || (u->_mp_d[limb_cnt] + & (((mp_limb_t) 1 << bit_index) - 1))); + else + adjust = 0; + + if (qn <= 0) + qn = 0; + + else + { + qp = MPZ_REALLOC (q, qn); + + if (bit_index != 0) + { + mpn_rshift (qp, u->_mp_d + limb_cnt, qn, bit_index); + qn -= qp[qn - 1] == 0; + } + else + { + mpn_copyi (qp, u->_mp_d + limb_cnt, qn); + } + } + + q->_mp_size = qn; + + mpz_add_ui (q, q, adjust); + if (un < 0) + mpz_neg (q, q); +} + +static void +mpz_div_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bit_index, + enum mpz_div_round_mode mode) +{ + mp_size_t us, un, rn; + mp_ptr rp; + mp_limb_t mask; + + us = u->_mp_size; + if (us == 0 || bit_index == 0) + { + r->_mp_size = 0; + return; + } + rn = (bit_index + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS; + assert (rn > 0); + + rp = MPZ_REALLOC (r, rn); + un = GMP_ABS (us); + + mask = GMP_LIMB_MAX >> (rn * GMP_LIMB_BITS - bit_index); + + if (rn > un) + { + /* Quotient (with truncation) is zero, and remainder is + non-zero */ + if ( (mode == DIV_FLOOR && us < 0) || (mode == DIV_CEIL && us > 0)) + { + /* Have to negate and sign extend. */ + mp_size_t i; + mp_limb_t cy; + + for (cy = 1, i = 0; i < un; i++) + { + mp_limb_t s = ~u->_mp_d[i] + cy; + cy = s < cy; + rp[i] = s; + } + assert (cy == 0); + for (; i < rn - 1; i++) + rp[i] = GMP_LIMB_MAX; + + rp[rn-1] = mask; + us = -us; + } + else + { + /* Just copy */ + if (r != u) + mpn_copyi (rp, u->_mp_d, un); + + rn = un; + } + } + else + { + if (r != u) + mpn_copyi (rp, u->_mp_d, rn - 1); + + rp[rn-1] = u->_mp_d[rn-1] & mask; + + if ( (mode == DIV_FLOOR && us < 0) || (mode == DIV_CEIL && us > 0)) + { + /* If r != 0, compute 2^{bit_count} - r. */ + mp_size_t i; + + for (i = 0; i < rn && rp[i] == 0; i++) + ; + if (i < rn) + { + /* r > 0, need to flip sign. */ + rp[i] = ~rp[i] + 1; + for (i++; i < rn; i++) + rp[i] = ~rp[i]; + + rp[rn-1] &= mask; + + /* us is not used for anything else, so we can modify it + here to indicate flipped sign. */ + us = -us; + } + } + } + rn = mpn_normalized_size (rp, rn); + r->_mp_size = us < 0 ? -rn : rn; +} + +void +mpz_cdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt) +{ + mpz_div_q_2exp (r, u, cnt, DIV_CEIL); +} + +void +mpz_fdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt) +{ + mpz_div_q_2exp (r, u, cnt, DIV_FLOOR); +} + +void +mpz_tdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt) +{ + mpz_div_q_2exp (r, u, cnt, DIV_TRUNC); +} + +void +mpz_cdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt) +{ + mpz_div_r_2exp (r, u, cnt, DIV_CEIL); +} + +void +mpz_fdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt) +{ + mpz_div_r_2exp (r, u, cnt, DIV_FLOOR); +} + +void +mpz_tdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt) +{ + mpz_div_r_2exp (r, u, cnt, DIV_TRUNC); +} + +void +mpz_divexact (mpz_t q, const mpz_t n, const mpz_t d) +{ + gmp_assert_nocarry (mpz_div_qr (q, NULL, n, d, DIV_TRUNC)); +} + +int +mpz_divisible_p (const mpz_t n, const mpz_t d) +{ + return mpz_div_qr (NULL, NULL, n, d, DIV_TRUNC) == 0; +} + +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) +{ + mp_size_t ns, qn; + mp_ptr qp; + mp_limb_t rl; + mp_size_t rs; + + ns = n->_mp_size; + if (ns == 0) + { + if (q) + q->_mp_size = 0; + if (r) + r->_mp_size = 0; + return 0; + } + + qn = GMP_ABS (ns); + if (q) + qp = MPZ_REALLOC (q, qn); + else + qp = NULL; + + rl = mpn_div_qr_1 (qp, n->_mp_d, qn, d); + assert (rl < d); + + rs = rl > 0; + rs = (ns < 0) ? -rs : rs; + + if (rl > 0 && ( (mode == DIV_FLOOR && ns < 0) + || (mode == DIV_CEIL && ns >= 0))) + { + if (q) + gmp_assert_nocarry (mpn_add_1 (qp, qp, qn, 1)); + rl = d - rl; + rs = -rs; + } + + if (r) + { + r->_mp_d[0] = rl; + r->_mp_size = rs; + } + if (q) + { + qn -= (qp[qn-1] == 0); + assert (qn == 0 || qp[qn-1] > 0); + + q->_mp_size = (ns < 0) ? - qn : qn; + } + + return rl; +} + +unsigned long +mpz_cdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (q, r, n, d, DIV_CEIL); +} + +unsigned long +mpz_fdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (q, r, n, d, DIV_FLOOR); +} + +unsigned long +mpz_tdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (q, r, n, d, DIV_TRUNC); +} + +unsigned long +mpz_cdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (q, NULL, n, d, DIV_CEIL); +} + +unsigned long +mpz_fdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (q, NULL, n, d, DIV_FLOOR); +} + +unsigned long +mpz_tdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (q, NULL, n, d, DIV_TRUNC); +} + +unsigned long +mpz_cdiv_ui (const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (NULL, NULL, n, d, DIV_CEIL); +} + +unsigned long +mpz_fdiv_ui (const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (NULL, NULL, n, d, DIV_FLOOR); +} + +unsigned long +mpz_tdiv_ui (const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (NULL, NULL, n, d, DIV_TRUNC); +} + +void +mpz_divexact_ui (mpz_t q, const mpz_t n, unsigned long d) +{ + gmp_assert_nocarry (mpz_div_qr_ui (q, NULL, n, d, DIV_TRUNC)); +} + +int +mpz_divisible_ui_p (const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (NULL, NULL, n, d, DIV_TRUNC) == 0; +} + + +/* GCD */ +static mp_limb_t +mpn_gcd_11 (mp_limb_t u, mp_limb_t v) +{ + unsigned shift; + + assert ( (u | v) > 0); + + if (u == 0) + return v; + else if (v == 0) + return u; + + gmp_ctz (shift, u | v); + + u >>= shift; + v >>= shift; + + if ( (u & 1) == 0) + MP_LIMB_T_SWAP (u, v); + + while ( (v & 1) == 0) + v >>= 1; + + while (u != v) + { + if (u > v) + { + u -= v; + do + u >>= 1; + while ( (u & 1) == 0); + } + else + { + v -= u; + do + v >>= 1; + while ( (v & 1) == 0); + } + } + return u << shift; +} + +unsigned long +mpz_gcd_ui (mpz_t g, const mpz_t u, unsigned long v) +{ + mp_size_t un; + + if (v == 0) + { + if (g) + mpz_abs (g, u); + } + else + { + un = GMP_ABS (u->_mp_size); + if (un != 0) + v = mpn_gcd_11 (mpn_div_qr_1 (NULL, u->_mp_d, un, v), v); + + if (g) + mpz_set_ui (g, v); + } + + return v; +} + +static mp_bitcnt_t +mpz_make_odd (mpz_t r, const mpz_t u) +{ + mp_size_t un, rn, i; + mp_ptr rp; + unsigned shift; + + un = GMP_ABS (u->_mp_size); + assert (un > 0); + + for (i = 0; u->_mp_d[i] == 0; i++) + ; + + gmp_ctz (shift, u->_mp_d[i]); + + rn = un - i; + rp = MPZ_REALLOC (r, rn); + if (shift > 0) + { + mpn_rshift (rp, u->_mp_d + i, rn, shift); + rn -= (rp[rn-1] == 0); + } + else + mpn_copyi (rp, u->_mp_d + i, rn); + + r->_mp_size = rn; + return i * GMP_LIMB_BITS + shift; +} + +void +mpz_gcd (mpz_t g, const mpz_t u, const mpz_t v) +{ + mpz_t tu, tv; + mp_bitcnt_t uz, vz, gz; + + if (u->_mp_size == 0) + { + mpz_abs (g, v); + return; + } + if (v->_mp_size == 0) + { + mpz_abs (g, u); + return; + } + + mpz_init (tu); + mpz_init (tv); + + uz = mpz_make_odd (tu, u); + vz = mpz_make_odd (tv, v); + gz = GMP_MIN (uz, vz); + + if (tu->_mp_size < tv->_mp_size) + mpz_swap (tu, tv); + + mpz_tdiv_r (tu, tu, tv); + if (tu->_mp_size == 0) + { + mpz_swap (g, tv); + } + else + for (;;) + { + int c; + + mpz_make_odd (tu, tu); + c = mpz_cmp (tu, tv); + if (c == 0) + { + mpz_swap (g, tu); + break; + } + if (c < 0) + mpz_swap (tu, tv); + + if (tv->_mp_size == 1) + { + mp_limb_t vl = tv->_mp_d[0]; + mp_limb_t ul = mpz_tdiv_ui (tu, vl); + mpz_set_ui (g, mpn_gcd_11 (ul, vl)); + break; + } + mpz_sub (tu, tu, tv); + } + mpz_clear (tu); + mpz_clear (tv); + mpz_mul_2exp (g, g, gz); +} + +void +mpz_gcdext (mpz_t g, mpz_t s, mpz_t t, const mpz_t u, const mpz_t v) +{ + mpz_t tu, tv, s0, s1, t0, t1; + mp_bitcnt_t uz, vz, gz; + mp_bitcnt_t power; + + if (u->_mp_size == 0) + { + /* g = 0 u + sgn(v) v */ + mpz_abs (g, v); + if (s) + mpz_set_ui (s, 0); + if (t) + mpz_set_si (t, mpz_sgn (v)); + return; + } + + if (v->_mp_size == 0) + { + /* g = sgn(u) u + 0 v */ + mpz_abs (g, u); + if (s) + mpz_set_si (s, mpz_sgn (u)); + if (t) + mpz_set_ui (t, 0); + return; + } + + mpz_init (tu); + mpz_init (tv); + mpz_init (s0); + mpz_init (s1); + mpz_init (t0); + mpz_init (t1); + + uz = mpz_make_odd (tu, u); + vz = mpz_make_odd (tv, v); + gz = GMP_MIN (uz, vz); + + uz -= gz; + vz -= gz; + + /* Cofactors corresponding to odd gcd. gz handled later. */ + if (tu->_mp_size < tv->_mp_size) + { + mpz_swap (tu, tv); + MPZ_SRCPTR_SWAP (u, v); + MPZ_PTR_SWAP (s, t); + MP_BITCNT_T_SWAP (uz, vz); + } + + /* Maintain + * + * u = t0 tu + t1 tv + * v = s0 tu + s1 tv + * + * where u and v denote the inputs with common factors of two + * eliminated, and det (s0, t0; s1, t1) = 2^p. Then + * + * 2^p tu = s1 u - t1 v + * 2^p tv = -s0 u + t0 v + */ + + /* After initial division, tu = q tv + tu', we have + * + * u = 2^uz (tu' + q tv) + * v = 2^vz tv + * + * or + * + * t0 = 2^uz, t1 = 2^uz q + * s0 = 0, s1 = 2^vz + */ + + mpz_setbit (t0, uz); + mpz_tdiv_qr (t1, tu, tu, tv); + mpz_mul_2exp (t1, t1, uz); + + mpz_setbit (s1, vz); + power = uz + vz; + + if (tu->_mp_size > 0) + { + mp_bitcnt_t shift; + shift = mpz_make_odd (tu, tu); + mpz_mul_2exp (t0, t0, shift); + mpz_mul_2exp (s0, s0, shift); + power += shift; + + for (;;) + { + int c; + c = mpz_cmp (tu, tv); + if (c == 0) + break; + + if (c < 0) + { + /* tv = tv' + tu + * + * u = t0 tu + t1 (tv' + tu) = (t0 + t1) tu + t1 tv' + * v = s0 tu + s1 (tv' + tu) = (s0 + s1) tu + s1 tv' */ + + mpz_sub (tv, tv, tu); + mpz_add (t0, t0, t1); + mpz_add (s0, s0, s1); + + shift = mpz_make_odd (tv, tv); + mpz_mul_2exp (t1, t1, shift); + mpz_mul_2exp (s1, s1, shift); + } + else + { + mpz_sub (tu, tu, tv); + mpz_add (t1, t0, t1); + mpz_add (s1, s0, s1); + + shift = mpz_make_odd (tu, tu); + mpz_mul_2exp (t0, t0, shift); + mpz_mul_2exp (s0, s0, shift); + } + power += shift; + } + } + + /* Now tv = odd part of gcd, and -s0 and t0 are corresponding + cofactors. */ + + mpz_mul_2exp (tv, tv, gz); + mpz_neg (s0, s0); + + /* 2^p g = s0 u + t0 v. Eliminate one factor of two at a time. To + adjust cofactors, we need u / g and v / g */ + + mpz_divexact (s1, v, tv); + mpz_abs (s1, s1); + mpz_divexact (t1, u, tv); + mpz_abs (t1, t1); + + while (power-- > 0) + { + /* s0 u + t0 v = (s0 - v/g) u - (t0 + u/g) v */ + if (mpz_odd_p (s0) || mpz_odd_p (t0)) + { + mpz_sub (s0, s0, s1); + mpz_add (t0, t0, t1); + } + mpz_divexact_ui (s0, s0, 2); + mpz_divexact_ui (t0, t0, 2); + } + + /* Arrange so that |s| < |u| / 2g */ + mpz_add (s1, s0, s1); + if (mpz_cmpabs (s0, s1) > 0) + { + mpz_swap (s0, s1); + mpz_sub (t0, t0, t1); + } + if (u->_mp_size < 0) + mpz_neg (s0, s0); + if (v->_mp_size < 0) + mpz_neg (t0, t0); + + mpz_swap (g, tv); + if (s) + mpz_swap (s, s0); + if (t) + mpz_swap (t, t0); + + mpz_clear (tu); + mpz_clear (tv); + mpz_clear (s0); + mpz_clear (s1); + mpz_clear (t0); + mpz_clear (t1); +} + +void +mpz_lcm (mpz_t r, const mpz_t u, const mpz_t v) +{ + mpz_t g; + + if (u->_mp_size == 0 || v->_mp_size == 0) + { + r->_mp_size = 0; + return; + } + + mpz_init (g); + + mpz_gcd (g, u, v); + mpz_divexact (g, u, g); + mpz_mul (r, g, v); + + mpz_clear (g); + mpz_abs (r, r); +} + +void +mpz_lcm_ui (mpz_t r, const mpz_t u, unsigned long v) +{ + if (v == 0 || u->_mp_size == 0) + { + r->_mp_size = 0; + return; + } + + v /= mpz_gcd_ui (NULL, u, v); + mpz_mul_ui (r, u, v); + + mpz_abs (r, r); +} + +int +mpz_invert (mpz_t r, const mpz_t u, const mpz_t m) +{ + mpz_t g, tr; + int invertible; + + if (u->_mp_size == 0 || mpz_cmpabs_ui (m, 1) <= 0) + return 0; + + mpz_init (g); + mpz_init (tr); + + mpz_gcdext (g, tr, NULL, u, m); + invertible = (mpz_cmp_ui (g, 1) == 0); + + if (invertible) + { + if (tr->_mp_size < 0) + { + if (m->_mp_size >= 0) + mpz_add (tr, tr, m); + else + mpz_sub (tr, tr, m); + } + mpz_swap (r, tr); + } + + mpz_clear (g); + mpz_clear (tr); + return invertible; +} + + +/* Higher level operations (sqrt and pow) */ + +/* Compute s = floor(sqrt(u)) and r = u - s^2. Allows r == NULL */ +void +mpz_sqrtrem (mpz_t s, mpz_t r, const mpz_t u) +{ + mpz_t x, t, dx, q; + + if (u->_mp_size < 0) + gmp_die ("mpz_sqrtrem: Negative argument."); + if (u->_mp_size == 0) + { + mpz_set_ui (s, 0); + mpz_set_ui (r, 0); + return; + } + + mpz_init (x); + mpz_init (t); + mpz_init (dx); + mpz_init (q); + + /* Make x > sqrt(a). This will be invariant through the loop. */ + mpz_setbit (x, (mpz_sizeinbase (u, 2) + 1) / 2); + + for (;;) + { + /* Compute x^2 - u */ + mpz_mul (t, x, x); + mpz_sub (t, t, u); + assert (t->_mp_size > 0); + + mpz_mul_2exp (dx, x, 1); + mpz_tdiv_q (q, t, dx); + if (q->_mp_size == 0) + break; + mpz_sub (x, x, q); + assert (x->_mp_size > 0); + } + /* We have 0 < u - x^2 < 2x, which implies that sqrt(a) < x < 1 + + sqrt(1+a), and x-1 = floor(sqrt(a)). Then r = a - (x-1)^2 = 2x - + 1 - t. */ + mpz_sub_ui (x, x, 1); + assert (x->_mp_size > 0); + + if (r) + { + mpz_sub_ui (dx, dx, 1); + mpz_sub (t, dx, t); + assert (t->_mp_size >= 0); + + mpz_swap (t, r); + } + mpz_swap (s, x); + + mpz_clear (x); + mpz_clear (t); + mpz_clear (dx); + mpz_clear (q); +} + +void +mpz_sqrt (mpz_t s, const mpz_t u) +{ + mpz_sqrtrem (s, NULL, u); +} + +void +mpz_ui_pow_ui (mpz_t r, unsigned long b, unsigned long e) +{ + unsigned long bit; + mpz_set_ui (r, 1); + + for (bit = GMP_ULONG_HIGHBIT; bit > 0; bit >>= 1) + { + mpz_mul (r, r, r); + if (e & bit) + mpz_mul_ui (r, r, b); + } +} + +void +mpz_powm (mpz_t r, const mpz_t b, const mpz_t e, const mpz_t m) +{ + mpz_t tr; + mpz_t base; + mp_size_t en, mn; + mp_srcptr mp; + struct gmp_div_inverse minv; + unsigned shift; + mp_ptr tp = NULL; + + en = GMP_ABS (e->_mp_size); + mn = GMP_ABS (m->_mp_size); + if (mn == 0) + gmp_die ("mpz_powm: Zero modulo."); + + if (en == 0) + { + mpz_set_ui (r, 1); + return; + } + + mp = m->_mp_d; + mpn_div_qr_invert (&minv, mp, mn); + shift = minv.shift; + + if (shift > 0) + { + /* To avoid shifts, we do all our reductions, except the final + one, using a *normalized* m. */ + minv.shift = 0; + + tp = gmp_xalloc_limbs (mn); + gmp_assert_nocarry (mpn_lshift (tp, mp, mn, shift)); + mp = tp; + } + + mpz_init (base); + + if (e->_mp_size < 0) + { + if (!mpz_invert (base, b, m)) + gmp_die ("mpz_powm: Negative exponent and non-invertibe base."); + } + else + { + mp_size_t bn; + mpz_abs (base, b); + + bn = base->_mp_size; + if (bn >= mn) + { + mpn_div_qr_preinv (NULL, base->_mp_d, base->_mp_size, mp, mn, &minv); + bn = mn; + } + + /* We have reduced the absolute value. Now take care of the + sign. Note that we get zero represented non-canonically as + m. */ + if (b->_mp_size < 0) + { + mp_ptr bp = MPZ_REALLOC (base, mn); + gmp_assert_nocarry (mpn_sub (bp, mp, mn, bp, bn)); + bn = mn; + } + base->_mp_size = mpn_normalized_size (base->_mp_d, bn); + } + mpz_init_set_ui (tr, 1); + + while (en-- > 0) + { + mp_limb_t w = e->_mp_d[en]; + mp_limb_t bit; + + for (bit = GMP_LIMB_HIGHBIT; bit > 0; bit >>= 1) + { + mpz_mul (tr, tr, tr); + if (w & bit) + mpz_mul (tr, tr, base); + if (tr->_mp_size > mn) + { + mpn_div_qr_preinv (NULL, tr->_mp_d, tr->_mp_size, mp, mn, &minv); + tr->_mp_size = mpn_normalized_size (tr->_mp_d, mn); + } + } + } + + /* Final reduction */ + if (tr->_mp_size >= mn) + { + minv.shift = shift; + mpn_div_qr_preinv (NULL, tr->_mp_d, tr->_mp_size, mp, mn, &minv); + tr->_mp_size = mpn_normalized_size (tr->_mp_d, mn); + } + if (tp) + gmp_free (tp); + + mpz_swap (r, tr); + mpz_clear (tr); + mpz_clear (base); +} + + +/* Logical operations and bit manipulation. */ + +/* Numbers are treated as if represented in two's complement (and + infinitely sign extended). For a negative values we get the two's + complement from -x = ~x + 1, where ~ is bitwise complementt. + Negation transforms + + xxxx10...0 + + into + + yyyy10...0 + + where yyyy is the bitwise complement of xxxx. So least significant + bits, up to and including the first one bit, are unchanged, and + the more significant bits are all complemented. + + To change a bit from zero to one in a negative number, subtract the + corresponding power of two from the absolute value. This can never + underflow. To change a bit from one to zero, add the corresponding + power of two, and this might overflow. E.g., if x = -001111, the + two's complement is 110001. Clearing the least significant bit, we + get two's complement 110000, and -010000. */ + +int +mpz_tstbit (const mpz_t d, mp_bitcnt_t bit_index) +{ + mp_size_t limb_index; + unsigned shift; + mp_size_t ds; + mp_size_t dn; + mp_limb_t w; + int bit; + + ds = d->_mp_size; + dn = GMP_ABS (ds); + limb_index = bit_index / GMP_LIMB_BITS; + if (limb_index >= dn) + return ds < 0; + + shift = bit_index % GMP_LIMB_BITS; + w = d->_mp_d[limb_index]; + bit = (w >> shift) & 1; + + if (ds < 0) + { + /* d < 0. Check if any of the bits below is set: If so, our bit + must be complemented. */ + if (shift > 0 && (w << (GMP_LIMB_BITS - shift)) > 0) + return bit ^ 1; + while (limb_index-- > 0) + if (d->_mp_d[limb_index] > 0) + return bit ^ 1; + } + return bit; +} + +static void +mpz_abs_add_bit (mpz_t d, mp_bitcnt_t bit_index) +{ + mp_size_t dn, limb_index; + mp_limb_t bit; + mp_ptr dp; + + dn = GMP_ABS (d->_mp_size); + + limb_index = bit_index / GMP_LIMB_BITS; + bit = (mp_limb_t) 1 << (bit_index % GMP_LIMB_BITS); + + if (limb_index >= dn) + { + mp_size_t i; + /* The bit should be set outside of the end of the number. + We have to increase the size of the number. */ + dp = MPZ_REALLOC (d, limb_index + 1); + + dp[limb_index] = bit; + for (i = dn; i < limb_index; i++) + dp[i] = 0; + dn = limb_index + 1; + } + else + { + mp_limb_t cy; + + dp = d->_mp_d; + + cy = mpn_add_1 (dp + limb_index, dp + limb_index, dn - limb_index, bit); + if (cy > 0) + { + dp = MPZ_REALLOC (d, dn + 1); + dp[dn++] = cy; + } + } + + d->_mp_size = (d->_mp_size < 0) ? - dn : dn; +} + +static void +mpz_abs_sub_bit (mpz_t d, mp_bitcnt_t bit_index) +{ + mp_size_t dn, limb_index; + mp_ptr dp; + mp_limb_t bit; + + dn = GMP_ABS (d->_mp_size); + dp = d->_mp_d; + + limb_index = bit_index / GMP_LIMB_BITS; + bit = (mp_limb_t) 1 << (bit_index % GMP_LIMB_BITS); + + assert (limb_index < dn); + + gmp_assert_nocarry (mpn_sub_1 (dp + limb_index, dp + limb_index, + dn - limb_index, bit)); + dn -= (dp[dn-1] == 0); + d->_mp_size = (d->_mp_size < 0) ? - dn : dn; +} + +void +mpz_setbit (mpz_t d, mp_bitcnt_t bit_index) +{ + if (!mpz_tstbit (d, bit_index)) + { + if (d->_mp_size >= 0) + mpz_abs_add_bit (d, bit_index); + else + mpz_abs_sub_bit (d, bit_index); + } +} + +void +mpz_clrbit (mpz_t d, mp_bitcnt_t bit_index) +{ + if (mpz_tstbit (d, bit_index)) + { + if (d->_mp_size >= 0) + mpz_abs_sub_bit (d, bit_index); + else + mpz_abs_add_bit (d, bit_index); + } +} + +void +mpz_combit (mpz_t d, mp_bitcnt_t bit_index) +{ + if (mpz_tstbit (d, bit_index) ^ (d->_mp_size < 0)) + mpz_abs_sub_bit (d, bit_index); + else + mpz_abs_add_bit (d, bit_index); +} + +void +mpz_com (mpz_t r, const mpz_t u) +{ + mpz_neg (r, u); + mpz_sub_ui (r, r, 1); +} + +void +mpz_and (mpz_t r, const mpz_t u, const mpz_t v) +{ + mp_size_t un, vn, rn, i; + mp_ptr up, vp, rp; + + mp_limb_t ux, vx, rx; + mp_limb_t uc, vc, rc; + mp_limb_t ul, vl, rl; + + un = GMP_ABS (u->_mp_size); + vn = GMP_ABS (v->_mp_size); + if (un < vn) + { + MPZ_SRCPTR_SWAP (u, v); + MP_SIZE_T_SWAP (un, vn); + } + if (vn == 0) + { + r->_mp_size = 0; + return; + } + + up = u->_mp_d; + vp = v->_mp_d; + + uc = u->_mp_size < 0; + vc = v->_mp_size < 0; + rc = uc & vc; + + ux = -uc; + vx = -vc; + rx = -rc; + + /* If the smaller input is positive, higher limbs don't matter. */ + rn = vx ? un : vn; + + rp = MPZ_REALLOC (r, rn + rc); + + for (i = 0; i < vn; i++) + { + ul = (up[i] ^ ux) + uc; + uc = ul < uc; + + vl = (vp[i] ^ vx) + vc; + vc = vl < vc; + + rl = ( (ul & vl) ^ rx) + rc; + rc = rl < rc; + rp[i] = rl; + } + assert (vc == 0); + + for (; i < rn; i++) + { + ul = (up[i] ^ ux) + uc; + uc = ul < uc; + + rl = ( (ul & vx) ^ rx) + rc; + rc = rl < rc; + rp[i] = rl; + } + if (rc) + rp[rn++] = rc; + else + rn = mpn_normalized_size (rp, rn); + + r->_mp_size = rx ? -rn : rn; +} + +void +mpz_ior (mpz_t r, const mpz_t u, const mpz_t v) +{ + mp_size_t un, vn, rn, i; + mp_ptr up, vp, rp; + + mp_limb_t ux, vx, rx; + mp_limb_t uc, vc, rc; + mp_limb_t ul, vl, rl; + + un = GMP_ABS (u->_mp_size); + vn = GMP_ABS (v->_mp_size); + if (un < vn) + { + MPZ_SRCPTR_SWAP (u, v); + MP_SIZE_T_SWAP (un, vn); + } + if (vn == 0) + { + mpz_set (r, u); + return; + } + up = u->_mp_d; + vp = v->_mp_d; + + uc = u->_mp_size < 0; + vc = v->_mp_size < 0; + rc = uc | vc; + + ux = -uc; + vx = -vc; + rx = -rc; + + /* If the smaller input is negative, by sign extension higher limbs + don't matter. */ + rn = vx ? vn : un; + + rp = MPZ_REALLOC (r, rn + rc); + for (i = 0; i < vn; i++) + { + ul = (up[i] ^ ux) + uc; + uc = ul < uc; + + vl = (vp[i] ^ vx) + vc; + vc = vl < vc; + + rl = ( (ul | vl) ^ rx) + rc; + rc = rl < rc; + rp[i] = rl; + } + assert (vc == 0); + + for (; i < rn; i++) + { + ul = (up[i] ^ ux) + uc; + uc = ul < uc; + + rl = ( (ul | vx) ^ rx) + rc; + rc = rl < rc; + rp[i] = rl; + } + if (rc) + rp[rn++] = rc; + else + rn = mpn_normalized_size (rp, rn); + + r->_mp_size = rx ? -rn : rn; +} + +void +mpz_xor (mpz_t r, const mpz_t u, const mpz_t v) +{ + mp_size_t un, vn, i; + mp_ptr up, vp, rp; + + mp_limb_t ux, vx, rx; + mp_limb_t uc, vc, rc; + mp_limb_t ul, vl, rl; + + un = GMP_ABS (u->_mp_size); + vn = GMP_ABS (v->_mp_size); + if (un < vn) + { + MPZ_SRCPTR_SWAP (u, v); + MP_SIZE_T_SWAP (un, vn); + } + if (vn == 0) + { + mpz_set (r, u); + return; + } + + up = u->_mp_d; + vp = v->_mp_d; + + uc = u->_mp_size < 0; + vc = v->_mp_size < 0; + rc = uc ^ vc; + + ux = -uc; + vx = -vc; + rx = -rc; + + rp = MPZ_REALLOC (r, un + rc); + for (i = 0; i < vn; i++) + { + ul = (up[i] ^ ux) + uc; + uc = ul < uc; + + vl = (vp[i] ^ vx) + vc; + vc = vl < vc; + + rl = (ul ^ vl ^ rx) + rc; + rc = rl < rc; + rp[i] = rl; + } + assert (vc == 0); + + for (; i < un; i++) + { + ul = (up[i] ^ ux) + uc; + uc = ul < uc; + + rl = (ul ^ ux) + rc; + rc = rl < rc; + rp[i] = rl; + } + if (rc) + rp[un++] = rc; + else + un = mpn_normalized_size (rp, un); + + r->_mp_size = rx ? -un : un; +} + +static unsigned +gmp_popcount_limb (mp_limb_t x) +{ + unsigned c; + + /* Do 16 bits at a time, to avoid limb-sized constants. */ + for (c = 0; x > 0; x >>= 16) + { + unsigned w = ((x >> 1) & 0x5555) + (x & 0x5555); + w = ((w >> 2) & 0x3333) + (w & 0x3333); + w = ((w >> 4) & 0x0f0f) + (w & 0x0f0f); + w = (w >> 8) + (w & 0x00ff); + c += w; + } + return c; +} + +mp_bitcnt_t +mpz_popcount (const mpz_t u) +{ + mp_size_t un, i; + mp_bitcnt_t c; + + un = u->_mp_size; + + if (un < 0) + return ~(mp_bitcnt_t) 0; + + for (c = 0, i = 0; i < un; i++) + c += gmp_popcount_limb (u->_mp_d[i]); + + return c; +} + +mp_bitcnt_t +mpz_hamdist (const mpz_t u, const mpz_t v) +{ + mp_size_t un, vn, i; + mp_limb_t uc, vc, ul, vl, comp; + mp_srcptr up, vp; + mp_bitcnt_t c; + + un = u->_mp_size; + vn = v->_mp_size; + + if ( (un ^ vn) < 0) + return ~(mp_bitcnt_t) 0; + + if (un < 0) + { + assert (vn < 0); + un = -un; + vn = -vn; + uc = vc = 1; + comp = - (mp_limb_t) 1; + } + else + uc = vc = comp = 0; + + up = u->_mp_d; + vp = v->_mp_d; + + if (un < vn) + MPN_SRCPTR_SWAP (up, un, vp, vn); + + for (i = 0, c = 0; i < vn; i++) + { + ul = (up[i] ^ comp) + uc; + uc = ul < uc; + + vl = (vp[i] ^ comp) + vc; + vc = vl < vc; + + c += gmp_popcount_limb (ul ^ vl); + } + assert (vc == 0); + + for (; i < vn; i++) + { + ul = (up[i] ^ comp) + uc; + uc = ul < uc; + + c += gmp_popcount_limb (ul ^ comp); + } + + return c; +} + +mp_bitcnt_t +mpz_scan1 (const mpz_t u, mp_bitcnt_t starting_bit) +{ + mp_ptr up; + mp_size_t us, un, i; + mp_limb_t limb, ux, uc; + unsigned cnt; + + up = u->_mp_d; + us = u->_mp_size; + un = GMP_ABS (us); + i = starting_bit / GMP_LIMB_BITS; + + /* Past the end there's no 1 bits for u>=0, or an immediate 1 bit + for u<0. Notice this test picks up any u==0 too. */ + if (i >= un) + return (us >= 0 ? ~(mp_bitcnt_t) 0 : starting_bit); + + if (us < 0) + { + ux = GMP_LIMB_MAX; + uc = mpn_zero_p (up, i); + } + else + ux = uc = 0; + + limb = (ux ^ up[i]) + uc; + uc = limb < uc; + + /* Mask to 0 all bits before starting_bit, thus ignoring them. */ + limb &= (GMP_LIMB_MAX << (starting_bit % GMP_LIMB_BITS)); + + while (limb == 0) + { + i++; + if (i == un) + { + assert (uc == 0); + /* For the u > 0 case, this can happen only for the first + masked limb. For the u < 0 case, it happens when the + highest limbs of the absolute value are all ones. */ + return (us >= 0 ? ~(mp_bitcnt_t) 0 : un * GMP_LIMB_BITS); + } + limb = (ux ^ up[i]) + uc; + uc = limb < uc; + } + gmp_ctz (cnt, limb); + return (mp_bitcnt_t) i * GMP_LIMB_BITS + cnt; +} + +mp_bitcnt_t +mpz_scan0 (const mpz_t u, mp_bitcnt_t starting_bit) +{ + mp_ptr up; + mp_size_t us, un, i; + mp_limb_t limb, ux, uc; + unsigned cnt; + + up = u->_mp_d; + us = u->_mp_size; + un = GMP_ABS (us); + i = starting_bit / GMP_LIMB_BITS; + + /* When past end, there's an immediate 0 bit for u>=0, or no 0 bits for + u<0. Notice this test picks up all cases of u==0 too. */ + if (i >= un) + return (us >= 0 ? starting_bit : ~(mp_bitcnt_t) 0); + + if (us < 0) + { + ux = GMP_LIMB_MAX; + uc = mpn_zero_p (up, i); + } + else + ux = uc = 0; + + limb = (ux ^ up[i]) + uc; + uc = limb < uc; + + /* Mask to 1 all bits before starting_bit, thus ignoring them. */ + limb |= ((mp_limb_t) 1 << (starting_bit % GMP_LIMB_BITS)) - 1; + + while (limb == GMP_LIMB_MAX) + { + i++; + if (i == un) + { + assert (uc == 0); + return (us >= 0 ? un * GMP_LIMB_BITS : ~(mp_bitcnt_t) 0); + } + limb = (ux ^ up[i]) + uc; + uc = limb < uc; + } + gmp_ctz (cnt, ~limb); + return (mp_bitcnt_t) i * GMP_LIMB_BITS + cnt; +} + + +/* MPZ base conversion. */ + +size_t +mpz_sizeinbase (const mpz_t u, int base) +{ + mp_size_t un; + mp_srcptr up; + mp_ptr tp; + mp_bitcnt_t bits; + struct gmp_div_inverse bi; + size_t ndigits; + + assert (base >= 2); + assert (base <= 36); + + un = GMP_ABS (u->_mp_size); + if (un == 0) + return 1; + + up = u->_mp_d; + + bits = (un - 1) * GMP_LIMB_BITS + mpn_limb_size_in_base_2 (up[un-1]); + switch (base) + { + case 2: + return bits; + case 4: + return (bits + 1) / 2; + case 8: + return (bits + 2) / 3; + case 16: + return (bits + 3) / 4; + case 32: + return (bits + 4) / 5; + /* FIXME: Do something more clever for the common case of base + 10. */ + } + + tp = gmp_xalloc_limbs (un); + mpn_copyi (tp, up, un); + mpn_div_qr_1_invert (&bi, base); + + for (ndigits = 0; un > 0; ndigits++) + { + mpn_div_qr_1_preinv (tp, tp, un, &bi); + un -= (tp[un-1] == 0); + } + gmp_free (tp); + return ndigits; +} + +char * +mpz_get_str (char *sp, int base, const mpz_t u) +{ + unsigned bits; + const char *digits; + mp_size_t un; + size_t i, sn; + + if (base >= 0) + { + digits = "0123456789abcdefghijklmnopqrstuvwxyz"; + } + else + { + base = -base; + digits = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ"; + } + if (base <= 1) + base = 10; + if (base > 36) + return NULL; + + sn = 1 + mpz_sizeinbase (u, base); + if (!sp) + sp = gmp_xalloc (1 + sn); + + un = GMP_ABS (u->_mp_size); + + if (un == 0) + { + sp[0] = '0'; + sp[1] = '\0'; + return sp; + } + + i = 0; + + if (u->_mp_size < 0) + sp[i++] = '-'; + + bits = mpn_base_power_of_two_p (base); + + if (bits) + /* Not modified in this case. */ + sn = i + mpn_get_str_bits ((unsigned char *) sp + i, bits, u->_mp_d, un); + else + { + struct mpn_base_info info; + mp_ptr tp; + + mpn_get_base_info (&info, base); + tp = gmp_xalloc_limbs (un); + mpn_copyi (tp, u->_mp_d, un); + + sn = i + mpn_get_str_other ((unsigned char *) sp + i, base, &info, tp, un); + gmp_free (tp); + } + + for (; i < sn; i++) + sp[i] = digits[(unsigned char) sp[i]]; + + sp[sn] = '\0'; + return sp; +} + +int +mpz_set_str (mpz_t r, const char *sp, int base) +{ + unsigned bits; + mp_size_t rn, alloc; + mp_ptr rp; + size_t sn; + size_t dn; + int sign; + unsigned char *dp; + + assert (base == 0 || (base >= 2 && base <= 36)); + + while (isspace( (unsigned char) *sp)) + sp++; + + if (*sp == '-') + { + sign = 1; + + do + sp++; + while (isspace( (unsigned char) *sp)); + } + else + sign = 0; + + if (base == 0) + { + if (*sp == '0') + { + sp++; + if (*sp == 'x') + { + base = 16; + sp++; + } + else if (*sp == 'b') + { + base = 2; + sp++; + } + else + base = 8; + } + else + base = 10; + } + + sn = strlen (sp); + dp = gmp_xalloc (sn); + + for (dn = 0; *sp; sp++) + { + unsigned digit; + + if (isspace ((unsigned char) *sp)) + continue; + if (*sp >= '0' && *sp <= '9') + digit = *sp - '0'; + else if (*sp >= 'a' && *sp <= 'z') + digit = *sp - 'a' + 10; + else if (*sp >= 'A' && *sp <= 'Z') + digit = *sp - 'A' + 10; + else + { + fail: + gmp_free (dp); + r->_mp_size = 0; + return -1; + } + + if (digit >= base) + goto fail; + + dp[dn++] = digit; + } + + bits = mpn_base_power_of_two_p (base); + + if (bits > 0) + { + alloc = (sn * bits + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS; + rp = MPZ_REALLOC (r, alloc); + rn = mpn_set_str_bits (rp, dp, dn, bits); + } + else + { + struct mpn_base_info info; + mpn_get_base_info (&info, base); + alloc = (sn + info.exp - 1) / info.exp; + rp = MPZ_REALLOC (r, alloc); + rn = mpn_set_str_other (rp, dp, dn, base, &info); + } + assert (rn <= alloc); + gmp_free (dp); + + r->_mp_size = sign ? - rn : rn; + + return 0; +} + +size_t +mpz_out_str (FILE *stream, int base, const mpz_t x) +{ + char *str; + size_t len; + + str = mpz_get_str (NULL, base, x); + len = strlen (str); + len = fwrite (str, 1, len, stream); + gmp_free (str); + return len; +} + + +static int +gmp_detect_endian (void) +{ + static const int i = 1; + const unsigned char *p = (const unsigned char *) &i; + if (*p == 1) + /* Little endian */ + return -1; + else + /* Big endian */ + return 1; +} + +/* Import and export. Does not support nails. */ +void +mpz_import (mpz_t r, size_t count, int order, size_t size, int endian, + size_t nails, const void *src) +{ + const unsigned char *p; + ptrdiff_t word_step; + mp_ptr rp; + mp_size_t rn; + + /* The current (partial) limb. */ + mp_limb_t limb; + /* The number of bytes already copied to this limb (starting from + the low end). */ + size_t bytes; + /* The index where the limb should be stored, when completed. */ + mp_size_t i; + + if (nails != 0) + gmp_die ("mpz_import: Nails not supported."); + + assert (order == 1 || order == -1); + assert (endian >= -1 && endian <= 1); + + if (endian == 0) + endian = gmp_detect_endian (); + + p = (unsigned char *) src; + + /* Process bytes from the least significant end, so point p at the + least significant word. */ + if (order == 1) + { + p += size * (count - 1); + word_step = -(ptrdiff_t) size; + } + else + word_step = size; + + /* And at east significant byte of that word. */ + if (endian == 1) + p += (size - 1); + + rn = (size * count + sizeof(mp_limb_t) - 1) / sizeof(mp_limb_t); + rp = MPZ_REALLOC (r, rn); + + for (limb = 0, bytes = 0, i = 0; count > 0; count--, p += word_step) + { + size_t j; + for (j = 0; j < size; j++, p -= (ptrdiff_t) endian) + { + limb |= *p << (bytes++ * CHAR_BIT); + if (bytes == sizeof(mp_limb_t)) + { + rp[i++] = limb; + bytes = 0; + } + } + } + if (bytes > 0) + rp[i++] = limb; + assert (i == rn); + + r->_mp_size = rn; +} + +void * +mpz_export (void *r, size_t *countp, int order, size_t size, int endian, + size_t nails, const mpz_t u) +{ + unsigned char *p; + ptrdiff_t word_step; + size_t count, k; + mp_size_t un; + + /* The current (partial) limb. */ + mp_limb_t limb; + /* The number of bytes left to to in this limb. */ + size_t bytes; + /* The index where the limb was read. */ + mp_size_t i; + + if (nails != 0) + gmp_die ("mpz_import: Nails not supported."); + + assert (order == 1 || order == -1); + assert (endian >= -1 && endian <= 1); + + un = GMP_ABS (u->_mp_size); + if (!un) + return r; + + count = (un * sizeof (mp_limb_t) + size - 1) / size; + if (!r) + r = gmp_xalloc (count * size); + + if (endian == 0) + endian = gmp_detect_endian (); + + p = (unsigned char *) r; + + /* Process bytes from the least significant end, so point p at the + least significant word. */ + if (order == 1) + { + p += size * (count - 1); + word_step = -(ptrdiff_t) size; + } + else + word_step = size; + + /* And at east significant byte of that word. */ + if (endian == 1) + p += (size - 1); + + for (limb = 0, bytes = 0, i = 0, k = 0; i < un; k++, p += word_step) + { + size_t j; + for (j = 0; j < size; j++, p -= (ptrdiff_t) endian) + { + if (bytes == 0) + { + if (i < un) + limb = u->_mp_d[i++]; + bytes = sizeof (mp_limb_t); + } + *p = limb; + limb >>= CHAR_BIT; + bytes--; + } + } + assert (i == un); + assert (k == count); + + if (countp) + *countp = count; + + return r; +} diff -r b641c8181188 mini-gmp/mini-gmp.h --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mini-gmp/mini-gmp.h Fri Feb 17 11:16:53 2012 +0100 @@ -0,0 +1,238 @@ +/* mini-gmp, a minimalistic implementation of a GNU GMP subset. + +Copyright 2011, 2012 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 3 of the License, or (at your +option) any later version. + +The GNU MP Library is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public +License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ + +/* About mini-gmp: This is a minimal implementation of a subset of the + GMP interface. It is intended for inclusion into applications which + have modest bignums needs, as a fallback when the real GMP library + is not installed. + + This file defines the public interface. */ + +#ifndef __MINI_GMP_H__ +#define __MINI_GMP_H__ + +/* For size_t */ +#include + +#if defined (__cplusplus) +extern "C" { +#endif + +void mp_set_memory_functions (void *(*) (size_t), + void *(*) (void *, size_t, size_t), + void (*) (void *, size_t)); + +void mp_get_memory_functions (void *(**) (size_t), + void *(**) (void *, size_t, size_t), + void (**) (void *, size_t)); + +typedef unsigned long mp_limb_t; +typedef long mp_size_t; +typedef unsigned long mp_bitcnt_t; + +typedef mp_limb_t *mp_ptr; +typedef const mp_limb_t *mp_srcptr; + +typedef struct +{ + int _mp_alloc; /* Number of *limbs* allocated and pointed + to by the _mp_d field. */ + int _mp_size; /* abs(_mp_size) is the number of limbs the + last field points to. If _mp_size is + negative this is a negative number. */ + mp_limb_t *_mp_d; /* Pointer to the limbs. */ +} __mpz_struct; + +typedef __mpz_struct mpz_t[1]; + +typedef __mpz_struct *mpz_ptr; +typedef const __mpz_struct *mpz_srcptr; + +void mpn_copyi (mp_ptr, mp_srcptr, mp_size_t); +void mpn_copyd (mp_ptr, mp_srcptr, mp_size_t); + +int mpn_cmp (mp_srcptr, mp_srcptr, mp_size_t); + +mp_limb_t mpn_add_1 (mp_ptr, mp_srcptr, mp_size_t, mp_limb_t); +mp_limb_t mpn_add_n (mp_ptr, mp_srcptr, mp_srcptr, mp_size_t); +mp_limb_t mpn_add (mp_ptr, mp_srcptr, mp_size_t, mp_srcptr,mp_size_t); + +mp_limb_t mpn_sub_1 (mp_ptr, mp_srcptr, mp_size_t, mp_limb_t); +mp_limb_t mpn_sub_n (mp_ptr, mp_srcptr, mp_srcptr, mp_size_t); +mp_limb_t mpn_sub (mp_ptr, mp_srcptr, mp_size_t, mp_srcptr,mp_size_t); + +mp_limb_t mpn_mul_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b); +mp_limb_t mpn_addmul_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b); +mp_limb_t mpn_submul_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b); + +mp_limb_t mpn_mul (mp_ptr rp, mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn); +void mpn_mul_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n); +void mpn_sqr (mp_ptr rp, mp_srcptr ap, mp_size_t n); + +mp_limb_t mpn_lshift (mp_ptr, mp_srcptr, mp_size_t, unsigned int); +mp_limb_t mpn_rshift (mp_ptr, mp_srcptr, mp_size_t, unsigned int); + +mp_limb_t mpn_invert_3by2 (mp_limb_t, mp_limb_t); +#define mpn_invert_limb(x) mpn_invert_3by2 ((x), 0) + +size_t mpn_get_str (unsigned char *, int, mp_ptr, mp_size_t); +mp_size_t mpn_set_str (mp_ptr, const unsigned char *, size_t, int); + +void mpz_init (mpz_t); +void mpz_init2 (mpz_t, mp_bitcnt_t); +void mpz_clear (mpz_t); + +#define mpz_odd_p(z) (((z)->_mp_size != 0) & (int) (z)->_mp_d[0]) +#define mpz_even_p(z) (! mpz_odd_p (z)) + +int mpz_sgn (const mpz_t); +int mpz_cmp_si (const mpz_t, long); +int mpz_cmp_ui (const mpz_t, unsigned long); +int mpz_cmp (const mpz_t, const mpz_t); +int mpz_cmpabs_ui (const mpz_t, unsigned long); +int mpz_cmpabs (const mpz_t, const mpz_t); +int mpz_cmp_d (const mpz_t, double); + +void mpz_abs (mpz_t, const mpz_t); +void mpz_neg (mpz_t, const mpz_t); +void mpz_swap (mpz_t, mpz_t); + +void mpz_add_ui (mpz_t, const mpz_t, unsigned long); +void mpz_add (mpz_t, const mpz_t, const mpz_t); +void mpz_sub_ui (mpz_t, const mpz_t, unsigned long); +void mpz_ui_sub (mpz_t, unsigned long, const mpz_t); +void mpz_sub (mpz_t, const mpz_t, const mpz_t); + +void mpz_mul_si (mpz_t, const mpz_t, long int); +void mpz_mul_ui (mpz_t, const mpz_t, unsigned long int); +void mpz_mul (mpz_t, const mpz_t, const mpz_t); +void mpz_mul_2exp (mpz_t, const mpz_t, mp_bitcnt_t); + +void mpz_cdiv_qr (mpz_t, mpz_t, const mpz_t, const mpz_t); +void mpz_fdiv_qr (mpz_t, mpz_t, const mpz_t, const mpz_t); +void mpz_tdiv_qr (mpz_t, mpz_t, const mpz_t, const mpz_t); +void mpz_cdiv_q (mpz_t, const mpz_t, const mpz_t); +void mpz_fdiv_q (mpz_t, const mpz_t, const mpz_t); +void mpz_tdiv_q (mpz_t, const mpz_t, const mpz_t); +void mpz_cdiv_r (mpz_t, const mpz_t, const mpz_t); +void mpz_fdiv_r (mpz_t, const mpz_t, const mpz_t); +void mpz_tdiv_r (mpz_t, const mpz_t, const mpz_t); + +void mpz_cdiv_q_2exp (mpz_t, const mpz_t, mp_bitcnt_t); +void mpz_fdiv_q_2exp (mpz_t, const mpz_t, mp_bitcnt_t); +void mpz_tdiv_q_2exp (mpz_t, const mpz_t, mp_bitcnt_t); +void mpz_cdiv_r_2exp (mpz_t, const mpz_t, mp_bitcnt_t); +void mpz_fdiv_r_2exp (mpz_t, const mpz_t, mp_bitcnt_t); +void mpz_tdiv_r_2exp (mpz_t, const mpz_t, mp_bitcnt_t); + +void mpz_divexact (mpz_t q, const mpz_t n, const mpz_t d); + +int mpz_divisible_p (const mpz_t, const mpz_t); + +unsigned long mpz_cdiv_qr_ui (mpz_t, mpz_t, const mpz_t, unsigned long); +unsigned long mpz_fdiv_qr_ui (mpz_t, mpz_t, const mpz_t, unsigned long); +unsigned long mpz_tdiv_qr_ui (mpz_t, mpz_t, const mpz_t, unsigned long); +unsigned long mpz_cdiv_q_ui (mpz_t, const mpz_t, unsigned long); +unsigned long mpz_fdiv_q_ui (mpz_t, const mpz_t, unsigned long); +unsigned long mpz_tdiv_q_ui (mpz_t, const mpz_t, unsigned long); +unsigned long mpz_cdiv_ui (const mpz_t, unsigned long); +unsigned long mpz_fdiv_ui (const mpz_t, unsigned long); +unsigned long mpz_tdiv_ui (const mpz_t, unsigned long); + +void mpz_divexact_ui (mpz_t q, const mpz_t n, unsigned long d); + +int mpz_divisible_ui_p (const mpz_t, unsigned long); + +unsigned long mpz_gcd_ui (mpz_t, const mpz_t, unsigned long); +void mpz_gcd (mpz_t, const mpz_t, const mpz_t); +void mpz_gcdext (mpz_t, mpz_t, mpz_t, const mpz_t, const mpz_t); +void mpz_lcm_ui (mpz_t, const mpz_t, unsigned long); +void mpz_lcm (mpz_t, const mpz_t, const mpz_t); +int mpz_invert (mpz_t, const mpz_t, const mpz_t); + +void mpz_sqrtrem (mpz_t, mpz_t, const mpz_t); +void mpz_sqrt (mpz_t, const mpz_t); + +void mpz_ui_pow_ui (mpz_t, unsigned long, unsigned long); +void mpz_powm (mpz_t, const mpz_t, const mpz_t, const mpz_t); + +int mpz_tstbit (const mpz_t, mp_bitcnt_t); +void mpz_setbit (mpz_t, mp_bitcnt_t); +void mpz_clrbit (mpz_t, mp_bitcnt_t); +void mpz_combit (mpz_t, mp_bitcnt_t); + +void mpz_com (mpz_t, const mpz_t); +void mpz_and (mpz_t, const mpz_t, const mpz_t); +void mpz_ior (mpz_t, const mpz_t, const mpz_t); +void mpz_xor (mpz_t, const mpz_t, const mpz_t); + +mp_bitcnt_t mpz_popcount (const mpz_t); +mp_bitcnt_t mpz_hamdist (const mpz_t, const mpz_t); +mp_bitcnt_t mpz_scan0 (const mpz_t, mp_bitcnt_t); +mp_bitcnt_t mpz_scan1 (const mpz_t, mp_bitcnt_t); + +int mpz_fits_slong_p (const mpz_t); +int mpz_fits_ulong_p (const mpz_t); +long int mpz_get_si (const mpz_t); +unsigned long int mpz_get_ui (const mpz_t); +double mpz_get_d (const mpz_t); +size_t mpz_size (const mpz_t); +mp_limb_t mpz_getlimbn (const mpz_t, mp_size_t); + +void mpz_set_si (mpz_t, signed long int); +void mpz_set_ui (mpz_t, unsigned long int); +void mpz_set (mpz_t, const mpz_t); +void mpz_set_d (mpz_t, double); + +void mpz_init_set_si (mpz_t, signed long int); +void mpz_init_set_ui (mpz_t, unsigned long int); +void mpz_init_set (mpz_t, const mpz_t); +void mpz_init_set_d (mpz_t, double); + +size_t mpz_sizeinbase (const mpz_t, int); +char *mpz_get_str (char *, int, const mpz_t); +int mpz_set_str (mpz_t, const char *, int); + +/* This long list taken from gmp.h. */ +/* For reference, "defined(EOF)" cannot be used here. In g++ 2.95.4, + defines EOF but not FILE. */ +#if defined (FILE) \ + || defined (H_STDIO) \ + || defined (_H_STDIO) /* AIX */ \ + || defined (_STDIO_H) /* glibc, Sun, SCO */ \ + || defined (_STDIO_H_) /* BSD, OSF */ \ + || defined (__STDIO_H) /* Borland */ \ + || defined (__STDIO_H__) /* IRIX */ \ + || defined (_STDIO_INCLUDED) /* HPUX */ \ + || defined (__dj_include_stdio_h_) /* DJGPP */ \ + || defined (_FILE_DEFINED) /* Microsoft */ \ + || defined (__STDIO__) /* Apple MPW MrC */ \ + || defined (_MSL_STDIO_H) /* Metrowerks */ \ + || defined (_STDIO_H_INCLUDED) /* QNX4 */ \ + || defined (_ISO_STDIO_ISO_H) /* Sun C++ */ +size_t mpz_out_str (FILE *, int, const mpz_t); +#endif + +void mpz_import (mpz_t, size_t, int, size_t, int, size_t, const void *); +void *mpz_export (void *, size_t *, int, size_t, int, size_t, const mpz_t); + +#if defined (__cplusplus) +} +#endif +#endif /* __MINI_GMP_H__ */ diff -r b641c8181188 mini-gmp/tests/Makefile --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mini-gmp/tests/Makefile Fri Feb 17 11:16:53 2012 +0100 @@ -0,0 +1,40 @@ +# Note: Requires GNU make + +CC = gcc +OPTFLAGS = -O +CFLAGS = $(OPTFLAGS) -Wall -g -I .. +LDFLAGS = + +LIBS = -lgmp -lm -lmcheck + +CHECK_PROGRAMS = t-add t-sub t-mul t-invert t-div t-div_2exp \ + t-double t-gcd t-lcm \ + t-sqrt t-powm t-logops t-bitops t-scan t-str + +MISC_OBJS = hex-random.o mini-random.o mini-gmp.o + +all: + +clean: + rm -f *.o $(CHECK_PROGRAMS) + +%: %.c +.c: + +# Keep object files +.PRECIOUS: %.o + +%.o: %.c ../mini-gmp.h hex-random.h mini-random.h + $(CC) $(CFLAGS) -c $< -o $@ + +mini-gmp.o: ../mini-gmp.c ../mini-gmp.h + $(CC) $(CFLAGS) -c ../mini-gmp.c -o mini-gmp.o + +%: %.o $(MISC_OBJS) + $(CC) $(LDFLAGS) $^ $(LIBS) -o $@ + +# Missing tests include: +# mpz_cmp_d, mpz_popcount, mpz_hamdist, mpz_ui_pow_ui + +check: $(CHECK_PROGRAMS) + ./run-tests $(CHECK_PROGRAMS) diff -r b641c8181188 mini-gmp/tests/hex-random.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mini-gmp/tests/hex-random.c Fri Feb 17 11:16:53 2012 +0100 @@ -0,0 +1,386 @@ +/* + +Copyright 2011, Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 3 of the License, or (at your +option) any later version. + +The GNU MP Library is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public +License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ + +#include +#include + +#include +#include + +#include "gmp.h" + +#include "hex-random.h" + +static gmp_randstate_t state; + +void +hex_random_init (void) +{ + unsigned long seed; + char *env_seed; + + env_seed = getenv("GMP_CHECK_RANDOMIZE"); + if (env_seed && env_seed[0]) + { + seed = strtoul (env_seed, NULL, 0); + if (seed) + printf ("Re-seeding with GMP_CHECK_RANDOMIZE=%lu\n", seed); + else + { + seed = time(NULL) + getpid(); + printf ("Seed GMP_CHECK_RANDOMIZE=%lu (include this in bug reports)\n", seed); + } + } + else + seed = 4711; + + gmp_randinit_default (state); + gmp_randseed_ui (state, seed); +} + +char * +hex_urandomb (unsigned long bits) +{ + char *res; + mpz_t x; + + mpz_init (x); + mpz_urandomb (x, state, bits); + gmp_asprintf (&res, "%Zx", x); + mpz_clear (x); + return res; +} + +char * +hex_rrandomb (unsigned long bits) +{ + char *res; + mpz_t x; + + mpz_init (x); + mpz_rrandomb (x, state, bits); + gmp_asprintf (&res, "%Zx", x); + mpz_clear (x); + return res; +} + +void +hex_random_op (enum hex_random_op op, unsigned long maxbits, + char **ap, char **bp, char **rp) +{ + mpz_t a, b, r; + unsigned long abits, bbits; + unsigned signs; + + mpz_init (a); + mpz_init (b); + mpz_init (r); + + abits = gmp_urandomb_ui (state, 32) % maxbits; + bbits = gmp_urandomb_ui (state, 32) % maxbits; + + mpz_rrandomb (a, state, abits); + mpz_rrandomb (b, state, bbits); + + signs = gmp_urandomb_ui (state, 3); + if (signs & 1) + mpz_neg (a, a); + if (signs & 2) + mpz_neg (b, b); + + switch (op) + { + default: + abort (); + case OP_ADD: + mpz_add (r, a, b); + break; + case OP_SUB: + mpz_sub (r, a, b); + break; + case OP_MUL: + mpz_mul (r, a, b); + break; + case OP_GCD: + if (signs & 4) + { + /* Produce a large gcd */ + unsigned long gbits = gmp_urandomb_ui (state, 32) % maxbits; + mpz_rrandomb (r, state, gbits); + mpz_mul (a, a, r); + mpz_mul (b, b, r); + } + mpz_gcd (r, a, b); + break; + case OP_LCM: + if (signs & 4) + { + /* Produce a large gcd */ + unsigned long gbits = gmp_urandomb_ui (state, 32) % maxbits; + mpz_rrandomb (r, state, gbits); + mpz_mul (a, a, r); + mpz_mul (b, b, r); + } + mpz_lcm (r, a, b); + break; + case OP_AND: + mpz_and (r, a, b); + break; + case OP_IOR: + mpz_ior (r, a, b); + break; + case OP_XOR: + mpz_xor (r, a, b); + break; + } + + gmp_asprintf (ap, "%Zx", a); + gmp_asprintf (bp, "%Zx", b); + gmp_asprintf (rp, "%Zx", r); + + mpz_clear (a); + mpz_clear (b); + mpz_clear (r); +} + +void +hex_random_op4 (enum hex_random_op op, unsigned long maxbits, + char **ap, char **bp, char **cp, char **dp) +{ + mpz_t a, b, c, d; + unsigned long abits, bbits; + unsigned signs; + + mpz_init (a); + mpz_init (b); + mpz_init (c); + mpz_init (d); + + if (op == OP_POWM) + { + unsigned long cbits; + abits = gmp_urandomb_ui (state, 32) % maxbits; + bbits = 1 + gmp_urandomb_ui (state, 32) % maxbits; + cbits = 2 + gmp_urandomb_ui (state, 32) % maxbits; + + mpz_rrandomb (a, state, abits); + mpz_rrandomb (b, state, bbits); + mpz_rrandomb (c, state, cbits); + + signs = gmp_urandomb_ui (state, 3); + if (signs & 1) + mpz_neg (a, a); + if (signs & 2) + { + mpz_t g; + + /* If we negate the exponent, must make sure that gcd(a, c) = 1 */ + if (mpz_sgn (a) == 0) + mpz_set_ui (a, 1); + else + { + mpz_init (g); + + for (;;) + { + mpz_gcd (g, a, c); + if (mpz_cmp_ui (g, 1) == 0) + break; + mpz_divexact (a, a, g); + } + mpz_clear (g); + } + mpz_neg (b, b); + } + if (signs & 4) + mpz_neg (c, c); + + mpz_powm (d, a, b, c); + } + else + { + unsigned long qbits; + bbits = 1 + gmp_urandomb_ui (state, 32) % maxbits; + qbits = gmp_urandomb_ui (state, 32) % maxbits; + abits = bbits + qbits; + if (abits > 30) + abits -= 30; + else + abits = 0; + + mpz_rrandomb (a, state, abits); + mpz_rrandomb (b, state, bbits); + + signs = gmp_urandomb_ui (state, 2); + if (signs & 1) + mpz_neg (a, a); + if (signs & 2) + mpz_neg (b, b); + + switch (op) + { + default: + abort (); + case OP_CDIV: + mpz_cdiv_qr (c, d, a, b); + break; + case OP_FDIV: + mpz_fdiv_qr (c, d, a, b); + break; + case OP_TDIV: + mpz_tdiv_qr (c, d, a, b); + break; + } + } + gmp_asprintf (ap, "%Zx", a); + gmp_asprintf (bp, "%Zx", b); + gmp_asprintf (cp, "%Zx", c); + gmp_asprintf (dp, "%Zx", d); + + mpz_clear (a); + mpz_clear (b); + mpz_clear (c); + mpz_clear (d); +} + +void +hex_random_bit_op (enum hex_random_op op, unsigned long maxbits, + char **ap, unsigned long *b, char **rp) +{ + mpz_t a, r; + unsigned long abits, bbits; + unsigned signs; + + mpz_init (a); + mpz_init (r); + + abits = gmp_urandomb_ui (state, 32) % maxbits; + bbits = gmp_urandomb_ui (state, 32) % (maxbits + 100); + + mpz_rrandomb (a, state, abits); + + signs = gmp_urandomb_ui (state, 1); + if (signs & 1) + mpz_neg (a, a); + + switch (op) + { + default: + abort (); + + case OP_SETBIT: + mpz_set (r, a); + mpz_setbit (r, bbits); + break; + case OP_CLRBIT: + mpz_set (r, a); + mpz_clrbit (r, bbits); + break; + case OP_COMBIT: + mpz_set (r, a); + mpz_combit (r, bbits); + break; + case OP_CDIV_Q_2: + mpz_cdiv_q_2exp (r, a, bbits); + break; + case OP_CDIV_R_2: + mpz_cdiv_r_2exp (r, a, bbits); + break; + case OP_FDIV_Q_2: + mpz_fdiv_q_2exp (r, a, bbits); + break; + case OP_FDIV_R_2: + mpz_fdiv_r_2exp (r, a, bbits); + break; + case OP_TDIV_Q_2: + mpz_tdiv_q_2exp (r, a, bbits); + break; + case OP_TDIV_R_2: + mpz_tdiv_r_2exp (r, a, bbits); + break; + } + + gmp_asprintf (ap, "%Zx", a); + *b = bbits; + gmp_asprintf (rp, "%Zx", r); + + mpz_clear (a); + mpz_clear (r); +} + +void +hex_random_scan_op (enum hex_random_op op, unsigned long maxbits, + char **ap, unsigned long *b, unsigned long *r) +{ + mpz_t a; + unsigned long abits, bbits; + unsigned signs; + + mpz_init (a); + + abits = gmp_urandomb_ui (state, 32) % maxbits; + bbits = gmp_urandomb_ui (state, 32) % (maxbits + 100); + + mpz_rrandomb (a, state, abits); + + signs = gmp_urandomb_ui (state, 1); + if (signs & 1) + mpz_neg (a, a); + + switch (op) + { + default: + abort (); + + case OP_SCAN0: + *r = mpz_scan0 (a, bbits); + break; + case OP_SCAN1: + *r = mpz_scan1 (a, bbits); + break; + } + gmp_asprintf (ap, "%Zx", a); + *b = bbits; + + mpz_clear (a); +} + +void +hex_random_str_op (unsigned long maxbits, + int base, char **ap, char **rp) +{ + mpz_t a; + unsigned long abits; + unsigned signs; + + mpz_init (a); + + abits = gmp_urandomb_ui (state, 32) % maxbits; + + mpz_rrandomb (a, state, abits); + + signs = gmp_urandomb_ui (state, 2); + if (signs & 1) + mpz_neg (a, a); + + *ap = mpz_get_str (NULL, 16, a); + *rp = mpz_get_str (NULL, base, a); + + mpz_clear (a); +} diff -r b641c8181188 mini-gmp/tests/hex-random.h --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mini-gmp/tests/hex-random.h Fri Feb 17 11:16:53 2012 +0100 @@ -0,0 +1,43 @@ +/* + +Copyright 2011, Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 3 of the License, or (at your +option) any later version. + +The GNU MP Library is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public +License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ + +enum hex_random_op + { + OP_ADD, OP_SUB, OP_MUL, OP_CDIV, OP_FDIV, OP_TDIV, + OP_CDIV_Q_2, OP_CDIV_R_2, + OP_FDIV_Q_2, OP_FDIV_R_2, + OP_TDIV_Q_2, OP_TDIV_R_2, + OP_GCD, OP_LCM, OP_POWM, OP_AND, OP_IOR, OP_XOR, + OP_SETBIT, OP_CLRBIT, OP_COMBIT, + OP_SCAN0, OP_SCAN1, + }; + +void hex_random_init (void); +char *hex_urandomb (unsigned long bits); +char *hex_rrandomb (unsigned long bits); +void hex_random_op (enum hex_random_op op, unsigned long maxbits, + char **ap, char **bp, char **rp); +void hex_random_op4 (enum hex_random_op op, unsigned long maxbits, + char **ap, char **bp, char **rp, char **qp); +void hex_random_bit_op (enum hex_random_op op, unsigned long maxbits, + char **ap, unsigned long *b, char **rp); +void hex_random_scan_op (enum hex_random_op op, unsigned long maxbits, + char **ap, unsigned long *b, unsigned long *r); +void hex_random_str_op (unsigned long maxbits, + int base, char **ap, char **rp); diff -r b641c8181188 mini-gmp/tests/mini-random.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mini-gmp/tests/mini-random.c Fri Feb 17 11:16:53 2012 +0100 @@ -0,0 +1,117 @@ +/* + +Copyright 2011, Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 3 of the License, or (at your +option) any later version. + +The GNU MP Library is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public +License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ + +#include +#include + +#include "mini-random.h" + +static void +set_str (mpz_t r, const char *s) +{ + if (mpz_set_str (r, s, 16) != 0) + { + fprintf (stderr, "mpz_set_str failed on input %s\n", s); + abort (); + } +} + +void +mini_urandomb (mpz_t r, unsigned long bits) +{ + char *s; + s = hex_urandomb (bits); + set_str (r, s); + free (s); +} + +void +mini_rrandomb (mpz_t r, unsigned long bits) +{ + char *s; + s = hex_rrandomb (bits); + set_str (r, s); + free (s); +} + +void +mini_random_op (enum hex_random_op op, unsigned long maxbits, + mpz_t a, mpz_t b, mpz_t r) +{ + char *ap; + char *bp; + char *rp; + + hex_random_op (op, maxbits, &ap, &bp, &rp); + set_str (a, ap); + set_str (b, bp); + set_str (r, rp); + + free (ap); + free (bp); + free (rp); +} + +void +mini_random_op4 (enum hex_random_op op, unsigned long maxbits, + mpz_t a, mpz_t b, mpz_t c, mpz_t d) +{ + char *ap; + char *bp; + char *cp; + char *dp; + + hex_random_op4 (op, maxbits, &ap, &bp, &cp, &dp); + set_str (a, ap); + set_str (b, bp); + set_str (c, cp); + set_str (d, dp); + + free (ap); + free (bp); + free (cp); + free (dp); +} + +void +mini_random_bit_op (enum hex_random_op op, unsigned long maxbits, + mpz_t a, mp_bitcnt_t *b, mpz_t r) +{ + char *ap; + char *rp; + + hex_random_bit_op (op, maxbits, &ap, b, &rp); + set_str (a, ap); + set_str (r, rp); + + free (ap); + free (rp); +} + +void +mini_random_scan_op (enum hex_random_op op, unsigned long maxbits, + mpz_t a, mp_bitcnt_t *b, mp_bitcnt_t *r) +{ + char *ap; + + hex_random_scan_op (op, maxbits, &ap, b, r); + set_str (a, ap); + + free (ap); +} diff -r b641c8181188 mini-gmp/tests/mini-random.h --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mini-gmp/tests/mini-random.h Fri Feb 17 11:16:53 2012 +0100 @@ -0,0 +1,32 @@ +/* + +Copyright 2011, Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 3 of the License, or (at your +option) any later version. + +The GNU MP Library is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public +License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ + +#include "mini-gmp.h" +#include "hex-random.h" + +void mini_urandomb (mpz_t r, unsigned long bits); +void mini_rrandomb (mpz_t r, unsigned long bits); +void mini_random_op (enum hex_random_op op, unsigned long maxbits, + mpz_t a, mpz_t b, mpz_t r); +void mini_random_op4 (enum hex_random_op op, unsigned long maxbits, + mpz_t a, mpz_t b, mpz_t c, mpz_t e); +void mini_random_scan_op (enum hex_random_op op, unsigned long maxbits, + mpz_t a, mp_bitcnt_t *b, mp_bitcnt_t *r); +void mini_random_bit_op (enum hex_random_op op, unsigned long maxbits, + mpz_t a, mp_bitcnt_t *b, mpz_t r); diff -r b641c8181188 mini-gmp/tests/run-tests --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mini-gmp/tests/run-tests Fri Feb 17 11:16:53 2012 +0100 @@ -0,0 +1,108 @@ +#! /bin/sh + +failed=0 +all=0 + +debug='no' +testflags='' + +if [ -z "$srcdir" ] ; then + srcdir=`pwd` +fi + +export srcdir + +# When used in make rules, we sometimes get the filenames VPATH +# expanded, but usually not. +find_program () { + case "$1" in + */*) + echo "$1" + ;; + *) + if [ -x "$1" ] ; then + echo "./$1" + else + echo "$srcdir/$1" + fi + ;; + esac +} + +env_program () { + if [ -x "$1" ] ; then + if "$1"; then : ; else + echo FAIL: $1 + exit 1 + fi + fi +} + +test_program () { + testname=`basename "$1" .exe` + testname=`basename "$testname" -test` + if [ -z "$EMULATOR" ] || head -1 "$1" | grep '^#!' > /dev/null; then + "$1" $testflags + else + "$EMULATOR" "$1" $testflags + fi + case "$?" in + 0) + echo PASS: $testname + all=`expr $all + 1` + ;; + 77) + echo SKIP: $testname + ;; + *) + echo FAIL: $testname + failed=`expr $failed + 1` + all=`expr $all + 1` + ;; + esac +} + +env_program `find_program setup-env` + +while test $# != 0 +do + case "$1" in + --debug) + debug=yes + ;; + -v) + testflags='-v' + ;; + -*) + echo >&2 'Unknown option `'"$1'" + exit 1 + ;; + *) + break + ;; + esac + shift +done + +if [ $# -eq 0 ] ; then + for f in *-test; do test_program "./$f"; done +else + for f in "$@" ; do test_program `find_program "$f"`; done +fi + +if [ $failed -eq 0 ] ; then + banner="All $all tests passed" +else + banner="$failed of $all tests failed" +fi +dashes=`echo "$banner" | sed s/./=/g` +echo "$dashes" +echo "$banner" +echo "$dashes" + +if [ "x$debug" = xno ] ; then + env_program `find_program teardown-env` +fi + +[ "$failed" -eq 0 ] + diff -r b641c8181188 mini-gmp/tests/t-add.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mini-gmp/tests/t-add.c Fri Feb 17 11:16:53 2012 +0100 @@ -0,0 +1,50 @@ +#include +#include + +#include "mini-random.h" + +#define MAXBITS 400 +#define COUNT 10000 + +static void +dump (const char *label, const mpz_t x) +{ + char *buf = mpz_get_str (NULL, 16, x); + fprintf (stderr, "%s: %s\n", label, buf); + free (buf); +} + +int +main (int argc, char **argv) +{ + unsigned i; + mpz_t a, b, res, ref; + + hex_random_init (); + + mpz_init (a); + mpz_init (b); + mpz_init (res); + mpz_init (ref); + + for (i = 0; i < COUNT; i++) + { + mini_random_op (OP_ADD, MAXBITS, a, b, ref); + mpz_add (res, a, b); + if (mpz_cmp (res, ref)) + { + fprintf (stderr, "mpz_add failed:\n"); + dump ("a", a); + dump ("b", b); + dump ("r", res); + dump ("ref", ref); + abort (); + } + } + mpz_clear (a); + mpz_clear (b); + mpz_clear (res); + mpz_clear (ref); + + return 0; +} diff -r b641c8181188 mini-gmp/tests/t-bitops.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mini-gmp/tests/t-bitops.c Fri Feb 17 11:16:53 2012 +0100 @@ -0,0 +1,96 @@ +#include +#include +#include + +#include "mini-random.h" + +#define MAXBITS 400 +#define COUNT 10000 + +static void +dump (const char *label, const mpz_t x) +{ + char *buf = mpz_get_str (NULL, 16, x); + fprintf (stderr, "%s: %s\n", label, buf); + free (buf); +} + +int +main (int argc, char **argv) +{ + unsigned i; + mpz_t a, res, ref; + mp_bitcnt_t b; + + hex_random_init (); + + mpz_init (a); + mpz_init (res); + mpz_init (ref); + + for (i = 0; i < COUNT; i++) + { + mini_random_bit_op (OP_SETBIT, MAXBITS, a, &b, ref); + mpz_set (res, a); + mpz_setbit (res, b); + if (mpz_cmp (res, ref)) + { + fprintf (stderr, "mpz_setbit failed:\n"); + dump ("a", a); + fprintf (stderr, "b: %lu\n", b); + dump ("r", res); + dump ("ref", ref); + abort (); + } + if (!mpz_tstbit (res, b)) + { + fprintf (stderr, "mpz_tstbit failed (after mpz_setbit):\n"); + dump ("res", a); + fprintf (stderr, "b: %lu\n", b); + abort (); + } + mini_random_bit_op (OP_CLRBIT, MAXBITS, a, &b, ref); + mpz_set (res, a); + mpz_clrbit (res, b); + if (mpz_cmp (res, ref)) + { + fprintf (stderr, "mpz_clrbit failed:\n"); + dump ("a", a); + fprintf (stderr, "b: %lu\n", b); + dump ("r", res); + dump ("ref", ref); + abort (); + } + if (mpz_tstbit (res, b)) + { + fprintf (stderr, "mpz_tstbit failed (after mpz_clrbit):\n"); + dump ("res", a); + fprintf (stderr, "b: %lu\n", b); + abort (); + } + mini_random_bit_op (OP_COMBIT, MAXBITS, a, &b, ref); + mpz_set (res, a); + mpz_combit (res, b); + if (mpz_cmp (res, ref)) + { + fprintf (stderr, "mpz_combit failed:\n"); + dump ("a", a); + fprintf (stderr, "b: %lu\n", b); + dump ("r", res); + dump ("ref", ref); + abort (); + } + if (mpz_tstbit (res, b) == mpz_tstbit (a, b)) + { + fprintf (stderr, "mpz_tstbit failed (after mpz_combit):\n"); + dump ("res", a); + fprintf (stderr, "b: %lu\n", b); + abort (); + } + } + mpz_clear (a); + mpz_clear (res); + mpz_clear (ref); + + return 0; +} diff -r b641c8181188 mini-gmp/tests/t-div.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mini-gmp/tests/t-div.c Fri Feb 17 11:16:53 2012 +0100 @@ -0,0 +1,94 @@ +#include +#include +#include + +#include "mini-random.h" + +#define MAXBITS 400 +#define COUNT 10000 + +static void +dump (const char *label, const mpz_t x) +{ + char *buf = mpz_get_str (NULL, 16, x); + fprintf (stderr, "%s: %s\n", label, buf); + free (buf); +} + +typedef void div_func (mpz_t, mpz_t, const mpz_t, const mpz_t); +typedef unsigned long div_ui_func (mpz_t, mpz_t, const mpz_t, unsigned long); + +int +main (int argc, char **argv) +{ + unsigned i; + mpz_t a, b, q, r, rq, rr; + + hex_random_init (); + + mpz_init (a); + mpz_init (b); + mpz_init (r); + mpz_init (q); + mpz_init (rr); + mpz_init (rq); + + for (i = 0; i < COUNT; i++) + { + unsigned j; + for (j = 0; j < 3; j++) + { + static const enum hex_random_op ops[3] = { OP_CDIV, OP_FDIV, OP_TDIV }; + static const char name[3] = { 'c', 'f', 't'}; + static div_func * const div [3] = + { + mpz_cdiv_qr, mpz_fdiv_qr, mpz_tdiv_qr + }; + static div_ui_func *div_ui[3] = + { + mpz_cdiv_qr_ui, mpz_fdiv_qr_ui, mpz_tdiv_qr_ui + }; + + mini_random_op4 (ops[j], MAXBITS, a, b, rq, rr); + div[j] (q, r, a, b); + if (mpz_cmp (r, rr) || mpz_cmp (q, rq)) + { + fprintf (stderr, "mpz_%cdiv_qr failed:\n", name[j]); + dump ("a", a); + dump ("b", b); + dump ("r ", r); + dump ("rref", rr); + dump ("q ", q); + dump ("qref", rq); + abort (); + } + if (mpz_fits_ulong_p (b)) + { + mp_limb_t rl; + rl = div_ui[j] (q, r, a, mpz_get_ui (b)); + + if (rl != mpz_get_ui (rr) + || mpz_cmp (r, rr) || mpz_cmp (q, rq)) + { + fprintf (stderr, "mpz_%cdiv_qr_ui failed:\n", name[j]); + dump ("a", a); + dump ("b", b); + fprintf(stderr, "rl = %lx\n", rl); + dump ("r ", r); + dump ("rref", rr); + dump ("q ", q); + dump ("qref", rq); + abort (); + } + } + } + } + mpz_clear (a); + mpz_clear (b); + mpz_clear (r); + mpz_clear (q); + mpz_clear (rr); + mpz_clear (rq); + + return 0; +} diff -r b641c8181188 mini-gmp/tests/t-div_2exp.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mini-gmp/tests/t-div_2exp.c Fri Feb 17 11:16:53 2012 +0100 @@ -0,0 +1,75 @@ +#include +#include +#include + +#include "mini-random.h" + +#define MAXBITS 400 +#define COUNT 10000 + +static void +dump (const char *label, const mpz_t x) +{ + char *buf = mpz_get_str (NULL, 16, x); + fprintf (stderr, "%s: %s\n", label, buf); + free (buf); +} + +typedef void div_func (mpz_t, const mpz_t, mp_bitcnt_t); + +int +main (int argc, char **argv) +{ + unsigned i; + mpz_t a, res, ref; + mp_bitcnt_t b; + + hex_random_init (); + + mpz_init (a); + mpz_init (res); + mpz_init (ref); + + for (i = 0; i < COUNT; i++) + { + unsigned j; + for (j = 0; j < 6; j++) + { + static const enum hex_random_op ops[6] = + { + OP_CDIV_Q_2, OP_CDIV_R_2, + OP_FDIV_Q_2, OP_FDIV_R_2, + OP_TDIV_Q_2, OP_TDIV_R_2 + }; + static const char *name[6] = + { + "cdiv_q", "cdiv_r", + "fdiv_q", "fdiv_r", + "tdiv_q", "tdiv_r" + }; + static div_func * const div [6] = + { + mpz_cdiv_q_2exp, mpz_cdiv_r_2exp, + mpz_fdiv_q_2exp, mpz_fdiv_r_2exp, + mpz_tdiv_q_2exp, mpz_tdiv_r_2exp + }; + + mini_random_bit_op (ops[j], MAXBITS, a, &b, ref); + div[j] (res, a, b); + if (mpz_cmp (ref, res)) + { + fprintf (stderr, "mpz_%s_2exp failed:\n", name[j]); + dump ("a", a); + fprintf (stderr, "b: %lu\n", b); + dump ("r", res); + dump ("ref", ref); + abort (); + } + } + } + mpz_clear (a); + mpz_clear (res); + mpz_clear (ref); + + return 0; +} diff -r b641c8181188 mini-gmp/tests/t-double.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mini-gmp/tests/t-double.c Fri Feb 17 11:16:53 2012 +0100 @@ -0,0 +1,108 @@ +#include +#include +#include +#include +#include + +#include "mini-random.h" + +#define GMP_LIMB_BITS (sizeof(mp_limb_t) * CHAR_BIT) + +#define COUNT 10000 + +static void +dump (const char *label, const mpz_t x) +{ + char *buf = mpz_get_str (NULL, 16, x); + fprintf (stderr, "%s: %s\n", label, buf); + free (buf); +} + +static const struct +{ + double d; + const char *s; +} values[] = { + { 0.0, "0" }, + { 0.3, "0" }, + { -0.3, "0" }, + { M_PI, "3" }, + { M_PI*1e15, "b29430a256d21" }, + { -M_PI*1e15, "-b29430a256d21" }, + /* 17 * 2^{200} = + 27317946752402834684213355569799764242877450894307478200123392 */ + {0.2731794675240283468421335556979976424288e62, + "1100000000000000000000000000000000000000000000000000" }, + { 0.0, NULL } +}; + +int +main (int argc, char **argv) +{ + unsigned i; + mpz_t x; + + hex_random_init (); + + mpz_init (x); + + for (i = 0; values[i].s; i++) + { + char *s; + mpz_set_d (x, values[i].d); + s = mpz_get_str (NULL, 16, x); + if (strcmp (s, values[i].s) != 0) + { + fprintf (stderr, "mpz_set_d failed:\n" + "d = %.20g\n" + "s = %s\n" + "r = %s\n", + values[i].d, s, values[i].s); + abort (); + } + free(s); + } + + for (i = 0; i < COUNT; i++) + { + double d, f; + unsigned long m; + int e; + + mini_rrandomb (x, GMP_LIMB_BITS); + m = mpz_get_ui (x); + mini_urandomb (x, 8); + e = mpz_get_ui (x) - 100; + + d = ldexp ((double) m, e); + mpz_set_d (x, d); + f = mpz_get_d (x); + if (f != floor (d)) + { + fprintf (stderr, "mpz_set_d/mpz_get_d failed:\n"); + dump ("x", x); + fprintf (stderr, "m = %lx, e = %i\n", m, e); + fprintf (stderr, "d = %.15g\n", d); + fprintf (stderr, "f = %.15g\n", f); + fprintf (stderr, "d - f = %.5g\n", d - f); + abort (); + } + d = - d; + + mpz_set_d (x, d); + f = mpz_get_d (x); + if (f != ceil (d)) + { + fprintf (stderr, "mpz_set_d/mpz_get_d failed:\n"); + dump ("x", x); + fprintf (stderr, "m = %lx, e = %i\n", m, e); + fprintf (stderr, "d = %.15g\n", d); + fprintf (stderr, "c = %.15g\n", f); + fprintf (stderr, "c - d = %.5g\n", f - d); + abort (); + } + } + + mpz_clear (x); + return EXIT_SUCCESS; +} diff -r b641c8181188 mini-gmp/tests/t-gcd.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mini-gmp/tests/t-gcd.c Fri Feb 17 11:16:53 2012 +0100 @@ -0,0 +1,164 @@ +#include +#include +#include + +#include "mini-random.h" + +#define MAXBITS 400 +#define COUNT 10000 + +static void +dump (const char *label, const mpz_t x) +{ + char *buf = mpz_get_str (NULL, 16, x); + fprintf (stderr, "%s: %s\n", label, buf); + free (buf); +} + +/* Called when g is supposed to be gcd(a,b), and g = s a + t b. */ +static int +gcdext_valid_p (const mpz_t a, const mpz_t b, + const mpz_t g, const mpz_t s, const mpz_t t) +{ + mpz_t ta, tb, r; + + /* It's not clear that gcd(0,0) is well defined, but we allow it and + require that gcd(0,0) = 0. */ + if (mpz_sgn (g) < 0) + return 0; + + if (mpz_sgn (a) == 0) + { + /* Must have g == abs (b). Any value for s is in some sense "correct", + but it makes sense to require that s == 0. */ + return mpz_cmpabs (g, b) == 0 && mpz_sgn (s) == 0; + } + else if (mpz_sgn (b) == 0) + { + /* Must have g == abs (a), s == sign (a) */ + return mpz_cmpabs (g, a) == 0 && mpz_cmp_si (s, mpz_sgn (a)) == 0; + } + + if (mpz_sgn (g) <= 0) + return 0; + + mpz_init (ta); + mpz_init (tb); + mpz_init (r); + + mpz_mul (ta, s, a); + mpz_mul (tb, t, b); + mpz_add (ta, ta, tb); + + if (mpz_cmp (ta, g) != 0) + { + fail: + mpz_clear (ta); + mpz_clear (tb); + mpz_clear (r); + return 0; + } + mpz_tdiv_qr (ta, r, a, g); + if (mpz_sgn (r) != 0) + goto fail; + + mpz_tdiv_qr (tb, r, b, g); + if (mpz_sgn (r) != 0) + goto fail; + + /* Require that 2 |s| < |b/g|, or |s| == 1. */ + if (mpz_cmpabs_ui (s, 1) > 0) + { + mpz_mul_2exp (r, s, 1); + if (mpz_cmpabs (r, tb) > 0) + goto fail; + } + + /* Require that 2 |t| < |a/g| or |t| == 1*/ + if (mpz_cmpabs_ui (t, 1) > 0) + { + mpz_mul_2exp (r, t, 1); + if (mpz_cmpabs (r, ta) > 0) + return 0; + } + return 1; +} + +int +main (int argc, char **argv) +{ + unsigned i; + mpz_t a, b, g, s, t; + + hex_random_init (); + + mpz_init (a); + mpz_init (b); + mpz_init (g); + mpz_init (s); + mpz_init (t); + + for (i = 0; i < COUNT; i++) + { + mini_random_op (OP_GCD, MAXBITS, a, b, s); + mpz_gcd (g, a, b); + if (mpz_cmp (g, s)) + { + fprintf (stderr, "mpz_gcd failed:\n"); + dump ("a", a); + dump ("b", b); + dump ("r", g); + dump ("ref", s); + abort (); + } + } + + for (i = 0; i < COUNT; i++) + { + unsigned flags; + mini_urandomb (a, 32); + flags = mpz_get_ui (a); + mini_rrandomb (a, MAXBITS); + mini_rrandomb (b, MAXBITS); + + if (flags % 37 == 0) + mpz_mul (a, a, b); + if (flags % 37 == 1) + mpz_mul (b, a, b); + + if (flags & 1) + mpz_neg (a, a); + if (flags & 2) + mpz_neg (b, b); + + mpz_gcdext (g, s, t, a, b); + if (!gcdext_valid_p (a, b, g, s, t)) + { + fprintf (stderr, "mpz_gcdext failed:\n"); + dump ("a", a); + dump ("b", b); + dump ("g", g); + dump ("s", s); + dump ("t", t); + abort (); + } + + mpz_gcd (s, a, b); + if (mpz_cmp (g, s)) + { + fprintf (stderr, "mpz_gcd failed:\n"); + dump ("a", a); + dump ("b", b); + dump ("r", g); + dump ("ref", s); + abort (); + } + } + mpz_clear (a); + mpz_clear (b); + mpz_clear (g); + mpz_clear (s); + mpz_clear (t); + + return 0; +} diff -r b641c8181188 mini-gmp/tests/t-invert.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mini-gmp/tests/t-invert.c Fri Feb 17 11:16:53 2012 +0100 @@ -0,0 +1,91 @@ +#include +#include +#include + +#include "mini-random.h" + +#define GMP_LIMB_BITS (sizeof(mp_limb_t) * CHAR_BIT) + +#define COUNT 10000 + +static void +dump (const char *label, const mpz_t x) +{ + char *buf = mpz_get_str (NULL, 16, x); + fprintf (stderr, "%s: %s\n", label, buf); + free (buf); +} + +int +main (int argc, char **argv) +{ + unsigned i; + mpz_t u, m, p, t; + + hex_random_init (); + + mpz_init (u); + mpz_init (m); + mpz_init (p); + mpz_init (t); + + for (i = 0; i < COUNT; i++) + { + mini_urandomb (u, GMP_LIMB_BITS); + mpz_setbit (u, GMP_LIMB_BITS -1); + + mpz_set_ui (m, mpn_invert_limb (u[0]._mp_d[0])); + mpz_setbit (m, GMP_LIMB_BITS); + + mpz_mul (p, m, u); + + mpz_set_ui (t, 0); + mpz_setbit (t, 2* GMP_LIMB_BITS); + mpz_sub (t, t, p); + + /* Should have 0 < B^2 - m u <= u */ + if (mpz_sgn (t) <= 0 || mpz_cmp (t, u) > 0) + { + fprintf (stderr, "mpn_invert_limb failed:\n"); + dump ("u", u); + dump ("m", m); + dump ("p", p); + dump ("t", t); + abort (); + } + } + + for (i = 0; i < COUNT; i++) + { + mini_urandomb (u, 2*GMP_LIMB_BITS); + mpz_setbit (u, 2*GMP_LIMB_BITS -1); + + mpz_set_ui (m, mpn_invert_3by2 (u[0]._mp_d[1], u[0]._mp_d[0])); + + mpz_setbit (m, GMP_LIMB_BITS); + + mpz_mul (p, m, u); + + mpz_set_ui (t, 0); + mpz_setbit (t, 3 * GMP_LIMB_BITS); + mpz_sub (t, t, p); + + /* Should have 0 < B^3 - m u <= u */ + if (mpz_sgn (t) <= 0 || mpz_cmp (t, u) > 0) + { + fprintf (stderr, "mpn_invert_3by2 failed:\n"); + dump ("u", u); + dump ("m", m); + dump ("p", p); + dump ("t", t); + abort (); + } + } + + mpz_clear (u); + mpz_clear (m); + mpz_clear (p); + mpz_clear (t); + + return 0; +} diff -r b641c8181188 mini-gmp/tests/t-lcm.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mini-gmp/tests/t-lcm.c Fri Feb 17 11:16:53 2012 +0100 @@ -0,0 +1,52 @@ +#include +#include +#include + +#include "mini-random.h" + +#define MAXBITS 400 +#define COUNT 10000 + +static void +dump (const char *label, const mpz_t x) +{ + char *buf = mpz_get_str (NULL, 16, x); + fprintf (stderr, "%s: %s\n", label, buf); + free (buf); +} + +int +main (int argc, char **argv) +{ + unsigned i; + mpz_t a, b, g, s; + + hex_random_init (); + + mpz_init (a); + mpz_init (b); + mpz_init (g); + mpz_init (s); + + for (i = 0; i < COUNT; i++) + { + mini_random_op (OP_LCM, MAXBITS, a, b, s); + mpz_lcm (g, a, b); + if (mpz_cmp (g, s)) + { + fprintf (stderr, "mpz_lcm failed:\n"); + dump ("a", a); + dump ("b", b); + dump ("r", g); + dump ("ref", s); + abort (); + } + } + + mpz_clear (a); + mpz_clear (b); + mpz_clear (g); + mpz_clear (s); + + return 0; +} diff -r b641c8181188 mini-gmp/tests/t-logops.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mini-gmp/tests/t-logops.c Fri Feb 17 11:16:53 2012 +0100 @@ -0,0 +1,75 @@ +#include +#include +#include + +#include "mini-random.h" + +#define MAXBITS 400 +#define COUNT 10000 + +static void +dump (const char *label, const mpz_t x) +{ + char *buf = mpz_get_str (NULL, 16, x); + fprintf (stderr, "%s: %s\n", label, buf); + free (buf); +} + +int +main (int argc, char **argv) +{ + unsigned i; + mpz_t a, b, res, ref; + + hex_random_init (); + + mpz_init (a); + mpz_init (b); + mpz_init (res); + mpz_init (ref); + + for (i = 0; i < COUNT; i++) + { + mini_random_op (OP_AND, MAXBITS, a, b, ref); + mpz_and (res, a, b); + if (mpz_cmp (res, ref)) + { + fprintf (stderr, "mpz_and failed:\n"); + dump ("a", a); + dump ("b", b); + dump ("r", res); + dump ("ref", ref); + abort (); + } + + mini_random_op (OP_IOR, MAXBITS, a, b, ref); + mpz_ior (res, a, b); + if (mpz_cmp (res, ref)) + { + fprintf (stderr, "mpz_ior failed:\n"); + dump ("a", a); + dump ("b", b); + dump ("r", res); + dump ("ref", ref); + abort (); + } + + mini_random_op (OP_XOR, MAXBITS, a, b, ref); + mpz_xor (res, a, b); + if (mpz_cmp (res, ref)) + { + fprintf (stderr, "mpz_xor failed:\n"); + dump ("a", a); + dump ("b", b); + dump ("r", res); + dump ("ref", ref); + abort (); + } + } + mpz_clear (a); + mpz_clear (b); + mpz_clear (res); + mpz_clear (ref); + + return 0; +} diff -r b641c8181188 mini-gmp/tests/t-mul.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mini-gmp/tests/t-mul.c Fri Feb 17 11:16:53 2012 +0100 @@ -0,0 +1,51 @@ +#include +#include +#include + +#include "mini-random.h" + +#define MAXBITS 400 +#define COUNT 10000 + +static void +dump (const char *label, const mpz_t x) +{ + char *buf = mpz_get_str (NULL, 16, x); + fprintf (stderr, "%s: %s\n", label, buf); + free (buf); +} + +int +main (int argc, char **argv) +{ + unsigned i; + mpz_t a, b, res, ref; + + hex_random_init (); + + mpz_init (a); + mpz_init (b); + mpz_init (res); + mpz_init (ref); + + for (i = 0; i < COUNT; i++) + { + mini_random_op (OP_MUL, MAXBITS, a, b, ref); + mpz_mul (res, a, b); + if (mpz_cmp (res, ref)) + { + fprintf (stderr, "mpz_mul failed:\n"); + dump ("a", a); + dump ("b", b); + dump ("r", res); + dump ("ref", ref); + abort (); + } + } + mpz_clear (a); + mpz_clear (b); + mpz_clear (res); + mpz_clear (ref); + + return 0; +} diff -r b641c8181188 mini-gmp/tests/t-powm.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mini-gmp/tests/t-powm.c Fri Feb 17 11:16:53 2012 +0100 @@ -0,0 +1,54 @@ +#include +#include +#include + +#include "mini-random.h" + +#define MAXBITS 400 +#define COUNT 1000 + +static void +dump (const char *label, const mpz_t x) +{ + char *buf = mpz_get_str (NULL, 16, x); + fprintf (stderr, "%s: %s\n", label, buf); + free (buf); +} + +int +main (int argc, char **argv) +{ + unsigned i; + mpz_t b, e, m, res, ref; + + hex_random_init (); + + mpz_init (b); + mpz_init (e); + mpz_init (m); + mpz_init (res); + mpz_init (ref); + + for (i = 0; i < COUNT; i++) + { + mini_random_op4 (OP_POWM, MAXBITS, b, e, m, ref); + mpz_powm (res, b, e, m); + if (mpz_cmp (res, ref)) + { + fprintf (stderr, "mpz_powm failed:\n"); + dump ("b", b); + dump ("e", e); + dump ("m", m); + dump ("r", res); + dump ("ref", ref); + abort (); + } + } + mpz_clear (b); + mpz_clear (e); + mpz_clear (m); + mpz_clear (res); + mpz_clear (ref); + + return 0; +} diff -r b641c8181188 mini-gmp/tests/t-scan.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mini-gmp/tests/t-scan.c Fri Feb 17 11:16:53 2012 +0100 @@ -0,0 +1,57 @@ +#include +#include +#include + +#include "mini-random.h" + +#define MAXBITS 400 +#define COUNT 10000 + +static void +dump (const char *label, const mpz_t x) +{ + char *buf = mpz_get_str (NULL, 16, x); + fprintf (stderr, "%s: %s\n", label, buf); + free (buf); +} + +int +main (int argc, char **argv) +{ + unsigned i; + mpz_t a; + mp_bitcnt_t b, res, ref; + + hex_random_init (); + + mpz_init (a); + + for (i = 0; i < COUNT; i++) + { + mini_random_scan_op (OP_SCAN0, MAXBITS, a, &b, &ref); + res = mpz_scan0 (a, b); + if (res != ref) + { + fprintf (stderr, "mpz_scan0 failed:\n"); + dump ("a", a); + fprintf (stderr, "b: %lu\n", b); + fprintf (stderr, "r: %lu\n", res); + fprintf (stderr, "ref: %lu\n", ref); + abort (); + } + mini_random_scan_op (OP_SCAN1, MAXBITS, a, &b, &ref); + res = mpz_scan1 (a, b); + if (res != ref) + { + fprintf (stderr, "mpz_scan1 failed:\n"); + dump ("a", a); + fprintf (stderr, "b: %lu\n", b); + fprintf (stderr, "r: %lu\n", res); + fprintf (stderr, "ref: %lu\n", ref); + abort (); + } + } + mpz_clear (a); + + return 0; +} diff -r b641c8181188 mini-gmp/tests/t-sqrt.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mini-gmp/tests/t-sqrt.c Fri Feb 17 11:16:53 2012 +0100 @@ -0,0 +1,75 @@ +#include +#include +#include + +#include "mini-random.h" + +#define MAXBITS 400 +#define COUNT 10000 + +static void +dump (const char *label, const mpz_t x) +{ + char *buf = mpz_get_str (NULL, 16, x); + fprintf (stderr, "%s: %s\n", label, buf); + free (buf); +} + +/* Called when s is supposed to be floor(sqrt(u)), and r = u - s^2 */ +static int +sqrtrem_valid_p (const mpz_t u, const mpz_t s, const mpz_t r) +{ + mpz_t t; + + mpz_init (t); + mpz_mul (t, s, s); + mpz_sub (t, u, t); + if (mpz_sgn (t) < 0 || mpz_cmp (t, r) != 0) + { + mpz_clear (t); + return 0; + } + mpz_add_ui (t, s, 1); + mpz_mul (t, t, t); + if (mpz_cmp (t, u) <= 0) + { + mpz_clear (t); + return 0; + } + + mpz_clear (t); + return 1; +} + +int +main (int argc, char **argv) +{ + unsigned i; + mpz_t u, s, r; + + hex_random_init (); + + mpz_init (u); + mpz_init (s); + mpz_init (r); + + for (i = 0; i < COUNT; i++) + { + mini_rrandomb (u, MAXBITS); + mpz_sqrtrem (s, r, u); + + if (!sqrtrem_valid_p (u, s, r)) + { + fprintf (stderr, "mpz_sqrtrem failed:\n"); + dump ("u", u); + dump ("sqrt", s); + dump ("rem", r); + abort (); + } + } + mpz_clear (u); + mpz_clear (s); + mpz_clear (r); + + return 0; +} diff -r b641c8181188 mini-gmp/tests/t-str.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mini-gmp/tests/t-str.c Fri Feb 17 11:16:53 2012 +0100 @@ -0,0 +1,92 @@ +#include +#include +#include +#include + +#include "mini-random.h" + +#define MAXBITS 400 +#define COUNT 2000 + +static void +dump (const char *label, const mpz_t x) +{ + char *buf = mpz_get_str (NULL, 16, x); + fprintf (stderr, "%s: %s\n", label, buf); + free (buf); +} + +int +main (int argc, char **argv) +{ + unsigned i; + char *ap; + char *bp; + char *rp; + size_t bn, rn, arn; + + mpz_t a, b; + + hex_random_init (); + + mpz_init (a); + mpz_init (b); + + for (i = 0; i < COUNT; i++) + { + int base; + for (base = 2; base <= 36; base++) + { + hex_random_str_op (MAXBITS, base, &ap, &rp); + if (mpz_set_str (a, ap, 16) != 0) + { + fprintf (stderr, "mpz_set_str failed on input %s\n", ap); + abort (); + } + + rn = strlen (rp); + arn = rn - (rp[0] == '-'); + + bn = mpz_sizeinbase (a, base); + if (bn < arn || bn > (arn + 1)) + { + fprintf (stderr, "mpz_sizeinbase failed:\n"); + dump ("a", a); + fprintf (stderr, "r = %s\n", rp); + fprintf (stderr, " base %d, correct size %u, got %u\n", + base, (unsigned) arn, (unsigned)bn); + abort (); + } + bp = mpz_get_str (NULL, base, a); + if (strcmp (bp, rp)) + { + fprintf (stderr, "mpz_get_str failed:\n"); + dump ("a", a); + fprintf (stderr, "b = %s\n", bp); + fprintf (stderr, " base = %d\n", base); + fprintf (stderr, "r = %s\n", rp); + abort (); + } + + mpz_set_str (b, rp, base); + + if (mpz_cmp (a, b)) + { + fprintf (stderr, "mpz_set_str failed:\n"); + fprintf (stderr, "r = %s\n", rp); + fprintf (stderr, " base = %d\n", base); + fprintf (stderr, "r = %s\n", ap); + fprintf (stderr, " base = 16\n"); + dump ("b", b); + dump ("r", a); + abort (); + } + free (ap); + free (bp); + } + } + mpz_clear (a); + mpz_clear (b); + + return 0; +} diff -r b641c8181188 mini-gmp/tests/t-sub.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mini-gmp/tests/t-sub.c Fri Feb 17 11:16:53 2012 +0100 @@ -0,0 +1,50 @@ +#include +#include + +#include "mini-random.h" + +#define MAXBITS 400 +#define COUNT 10000 + +static void +dump (const char *label, const mpz_t x) +{ + char *buf = mpz_get_str (NULL, 16, x); + fprintf (stderr, "%s: %s\n", label, buf); + free (buf); +} + +int +main (int argc, char **argv) +{ + unsigned i; + mpz_t a, b, res, ref; + + hex_random_init (); + + mpz_init (a); + mpz_init (b); + mpz_init (res); + mpz_init (ref); + + for (i = 0; i < COUNT; i++) + { + mini_random_op (OP_SUB, MAXBITS, a, b, ref); + mpz_sub (res, a, b); + if (mpz_cmp (res, ref)) + { + fprintf (stderr, "mpz_sub failed:\n"); + dump ("a", a); + dump ("b", b); + dump ("r", res); + dump ("ref", ref); + abort (); + } + } + mpz_clear (a); + mpz_clear (b); + mpz_clear (res); + mpz_clear (ref); + + return 0; +}