diff options
Diffstat (limited to 'ghc/rts/gmp/mpz')
122 files changed, 4664 insertions, 980 deletions
diff --git a/ghc/rts/gmp/mpz/Makefile.am b/ghc/rts/gmp/mpz/Makefile.am new file mode 100644 index 0000000000..cd6fec4e21 --- /dev/null +++ b/ghc/rts/gmp/mpz/Makefile.am @@ -0,0 +1,58 @@ +## Process this file with automake to generate Makefile.in + +# Copyright (C) 1996, 1998, 1999, 2000 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 2.1 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; see the file COPYING.LIB. If not, write to +# the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +# MA 02111-1307, USA. + + +AUTOMAKE_OPTIONS = gnu no-dependencies + +SUBDIRS = tests + +INCLUDES = -I$(top_srcdir) -DOPERATION_$* + +noinst_LTLIBRARIES = libmpz.la +libmpz_la_SOURCES = \ + abs.c add.c add_ui.c addmul_ui.c and.c array_init.c \ + bin_ui.c bin_uiui.c cdiv_q.c \ + cdiv_q_ui.c cdiv_qr.c cdiv_qr_ui.c cdiv_r.c cdiv_r_ui.c cdiv_ui.c \ + clear.c clrbit.c cmp.c cmp_si.c cmp_ui.c cmpabs.c cmpabs_ui.c com.c \ + divexact.c dump.c fac_ui.c fdiv_q.c fdiv_q_2exp.c fdiv_q_ui.c \ + fdiv_qr.c fdiv_qr_ui.c fdiv_r.c fdiv_r_2exp.c fdiv_r_ui.c fdiv_ui.c \ + fib_ui.c fits_sint_p.c fits_slong_p.c fits_sshort_p.c fits_uint_p.c \ + fits_ulong_p.c fits_ushort_p.c gcd.c gcd_ui.c gcdext.c get_d.c get_si.c \ + get_str.c get_ui.c getlimbn.c hamdist.c init.c inp_raw.c inp_str.c \ + invert.c ior.c iset.c iset_d.c iset_si.c iset_str.c iset_ui.c \ + jacobi.c kronsz.c kronuz.c kronzs.c kronzu.c \ + lcm.c legendre.c mod.c mul.c mul_2exp.c neg.c nextprime.c \ + out_raw.c out_str.c perfpow.c perfsqr.c popcount.c pow_ui.c powm.c \ + powm_ui.c pprime_p.c random.c random2.c realloc.c remove.c root.c rrandomb.c \ + scan0.c scan1.c set.c set_d.c set_f.c set_q.c set_si.c set_str.c \ + set_ui.c setbit.c size.c sizeinbase.c sqrt.c sqrtrem.c sub.c \ + sub_ui.c swap.c tdiv_ui.c tdiv_q.c tdiv_q_2exp.c tdiv_q_ui.c tdiv_qr.c \ + tdiv_qr_ui.c tdiv_r.c tdiv_r_2exp.c tdiv_r_ui.c tstbit.c ui_pow_ui.c \ + urandomb.c urandomm.c xor.c + +EXTRA_DIST = mul_siui.c +nodist_libmpz_la_SOURCES = mul_si.c mul_ui.c +CLEANFILES = $(nodist_libmpz_la_SOURCES) + +mul_si.c: $(srcdir)/mul_siui.c + cp $(srcdir)/mul_siui.c mul_si.c +mul_ui.c: $(srcdir)/mul_siui.c + cp $(srcdir)/mul_siui.c mul_ui.c diff --git a/ghc/rts/gmp/mpz/README b/ghc/rts/gmp/mpz/README new file mode 100644 index 0000000000..06b481d770 --- /dev/null +++ b/ghc/rts/gmp/mpz/README @@ -0,0 +1,23 @@ +This directory contains functions for GMP's integer function layer. + +In this version of GMP, integers are represented like in the figure below. +(Please note that the format might change between every version, and that +depending on the internal format in any way is a bad idea.) + + most least +significant significant + limb limb + + _mp_d + / + / + \/ + ____ ____ ____ ____ ____ + |____|____|____|____|____| + + <------- _mp_size -------> + + +The most significant limb will be non-zero. The _mp_size field's sign +reflects the sign of the number. Its absolute value is the count of limbs +in the number. diff --git a/ghc/rts/gmp/mpz/abs.c b/ghc/rts/gmp/mpz/abs.c index 080cac60aa..0b5eab1ce6 100644 --- a/ghc/rts/gmp/mpz/abs.c +++ b/ghc/rts/gmp/mpz/abs.c @@ -5,16 +5,16 @@ Copyright (C) 1991, 1993, 1994, 1995 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/add.c b/ghc/rts/gmp/mpz/add.c index 10dd9704c2..a22c3778fb 100644 --- a/ghc/rts/gmp/mpz/add.c +++ b/ghc/rts/gmp/mpz/add.c @@ -1,26 +1,29 @@ /* mpz_add -- Add two integers. -Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc. +Copyright (C) 1991, 1993, 1994, 1996, 2000 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ #include "gmp.h" #include "gmp-impl.h" +#ifdef BERKELEY_MP +#include "mp.h" +#endif #ifndef BERKELEY_MP void @@ -58,9 +61,9 @@ madd (u, v, w) if (abs_usize < abs_vsize) { /* Swap U and V. */ - {const __mpz_struct *t = u; u = v; v = t;} - {mp_size_t t = usize; usize = vsize; vsize = t;} - {mp_size_t t = abs_usize; abs_usize = abs_vsize; abs_vsize = t;} + MPZ_SRCPTR_SWAP (u, v); + MP_SIZE_T_SWAP (usize, vsize); + MP_SIZE_T_SWAP (abs_usize, abs_vsize); } /* True: ABS_USIZE >= ABS_VSIZE. */ diff --git a/ghc/rts/gmp/mpz/add_ui.c b/ghc/rts/gmp/mpz/add_ui.c index a1e4306319..28dbd71f45 100644 --- a/ghc/rts/gmp/mpz/add_ui.c +++ b/ghc/rts/gmp/mpz/add_ui.c @@ -1,20 +1,20 @@ /* mpz_add_ui -- Add an mpz_t and an unsigned one-word integer. -Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc. +Copyright (C) 1991, 1993, 1994, 1996, 1999 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -59,7 +59,7 @@ mpz_add_ui (w, u, v) if (usize >= 0) { mp_limb_t cy; - cy = mpn_add_1 (wp, up, abs_usize, v); + cy = mpn_add_1 (wp, up, abs_usize, (mp_limb_t) v); wp[abs_usize] = cy; wsize = abs_usize + cy; } @@ -74,7 +74,7 @@ mpz_add_ui (w, u, v) } else { - mpn_sub_1 (wp, up, abs_usize, v); + mpn_sub_1 (wp, up, abs_usize, (mp_limb_t) v); /* Size can decrease with at most one limb. */ wsize = -(abs_usize - (wp[abs_usize - 1] == 0)); } diff --git a/ghc/rts/gmp/mpz/addmul_ui.c b/ghc/rts/gmp/mpz/addmul_ui.c new file mode 100644 index 0000000000..7b38d3624d --- /dev/null +++ b/ghc/rts/gmp/mpz/addmul_ui.c @@ -0,0 +1,214 @@ +/* mpz_addmul_ui(prodsum, multiplier, small_multiplicand) -- + Add MULTIPLICATOR times SMALL_MULTIPLICAND to PRODSUM. + +Copyright (C) 1997, 2000 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 2.1 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; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +static mp_limb_t mpn_neg1 _PROTO ((mp_ptr, mp_size_t)); + +#if 0 +#undef MPN_NORMALIZE +#define MPN_NORMALIZE(DST, NLIMBS) \ + do { \ + while (--(NLIMBS) >= 0 && (DST)[NLIMBS] == 0) \ + ; \ + (NLIMBS)++; \ + } while (0) +#undef MPN_NORMALIZE_NOT_ZERO +#define MPN_NORMALIZE_NOT_ZERO(DST, NLIMBS) \ + do { \ + while ((DST)[--(NLIMBS)] == 0) \ + ; \ + (NLIMBS)++; \ + } while (0) +#endif + +void +#if __STDC__ +mpz_addmul_ui (mpz_ptr rz, mpz_srcptr az, unsigned long int bu) +#else +mpz_addmul_ui (rz, az, bu) + mpz_ptr rz; + mpz_srcptr az; + unsigned long int bu; +#endif +{ + mp_size_t rn, an; + mp_ptr rp, ap; + + an = SIZ (az); + + /* If either multiplier is zero, result is unaffected. */ + if (bu == 0 || an == 0) + return; + + rn = SIZ (rz); + + if (rn == 0) + { + mp_limb_t cy; + + an = ABS (an); + if (ALLOC (rz) <= an) + _mpz_realloc (rz, an + 1); + rp = PTR (rz); + ap = PTR (az); + cy = mpn_mul_1 (rp, ap, an, (mp_limb_t) bu); + rp[an] = cy; + an += cy != 0; + SIZ (rz) = SIZ (az) >= 0 ? an : -an; + return; + } + + if ((an ^ rn) >= 0) + { + /* Sign of operands are the same--really add. */ + an = ABS (an); + rn = ABS (rn); + if (rn > an) + { + mp_limb_t cy; + if (ALLOC (rz) <= rn) + _mpz_realloc (rz, rn + 1); + rp = PTR (rz); + ap = PTR (az); + cy = mpn_addmul_1 (rp, ap, an, (mp_limb_t) bu); + cy = mpn_add_1 (rp + an, rp + an, rn - an, cy); + rp[rn] = cy; + rn += cy != 0; + SIZ (rz) = SIZ (rz) >= 0 ? rn : -rn; + return; + } + else + { + mp_limb_t cy; + if (ALLOC (rz) <= an) + _mpz_realloc (rz, an + 1); + rp = PTR (rz); + ap = PTR (az); + cy = mpn_addmul_1 (rp, ap, rn, (mp_limb_t) bu); + if (an != rn) + { + mp_limb_t cy2; + cy2 = mpn_mul_1 (rp + rn, ap + rn, an - rn, (mp_limb_t) bu); + cy = cy2 + mpn_add_1 (rp + rn, rp + rn, an - rn, cy); + } + rn = an; + rp[rn] = cy; + rn += cy != 0; + SIZ (rz) = SIZ (rz) >= 0 ? rn : -rn; + return; + } + } + else + { + /* Sign of operands are different--actually subtract. */ + an = ABS (an); + rn = ABS (rn); + if (rn > an) + { + mp_limb_t cy; + rp = PTR (rz); + ap = PTR (az); + cy = mpn_submul_1 (rp, ap, an, (mp_limb_t) bu); + cy = mpn_sub_1 (rp + an, rp + an, rn - an, cy); + if (cy != 0) + { + mpn_neg1 (rp, rn); + MPN_NORMALIZE_NOT_ZERO (rp, rn); + } + else + { + MPN_NORMALIZE (rp, rn); + rn = -rn; + } + + SIZ (rz) = SIZ (rz) >= 0 ? -rn : rn; + return; + } + else + { + /* Tricky case. We need to subtract an operand that might be larger + than the minuend. To avoid allocating temporary space, we compute + a*b-r instead of r-a*b and then negate. */ + mp_limb_t cy; + if (ALLOC (rz) <= an) + _mpz_realloc (rz, an + 1); + rp = PTR (rz); + ap = PTR (az); + cy = mpn_submul_1 (rp, ap, rn, (mp_limb_t) bu); + if (an != rn) + { + mp_limb_t cy2; + cy -= mpn_neg1 (rp, rn); + cy2 = mpn_mul_1 (rp + rn, ap + rn, an - rn, (mp_limb_t) bu); + if (cy == ~(mp_limb_t) 0) + cy = cy2 - mpn_sub_1 (rp + rn, rp + rn, an - rn, (mp_limb_t) 1); + else + cy = cy2 + mpn_add_1 (rp + rn, rp + rn, an - rn, cy); + rp[an] = cy; + rn = an + (cy != 0); + rn -= rp[rn - 1] == 0; + } + else if (cy != 0) + { + cy -= mpn_neg1 (rp, rn); + rp[an] = cy; + rn = an + 1; + MPN_NORMALIZE_NOT_ZERO (rp, rn); + } + else + { + rn = an; + MPN_NORMALIZE (rp, rn); + rn = -rn; + } + + SIZ (rz) = SIZ (rz) >= 0 ? -rn : rn; + return; + } + } +} + +static mp_limb_t +#if __STDC__ +mpn_neg1 (mp_ptr rp, mp_size_t rn) +#else +mpn_neg1 (rp, rn) + mp_ptr rp; + mp_size_t rn; +#endif +{ + mp_size_t i; + + while (rn != 0 && rp[0] == 0) + rp++, rn--; + + if (rn != 0) + { + rp[0] = -rp[0]; + for (i = 1; i < rn; i++) + rp[i] = ~rp[i]; + return 1; + } + return 0; +} diff --git a/ghc/rts/gmp/mpz/and.c b/ghc/rts/gmp/mpz/and.c index 838d4b1db1..354e9455bf 100644 --- a/ghc/rts/gmp/mpz/and.c +++ b/ghc/rts/gmp/mpz/and.c @@ -1,20 +1,21 @@ /* mpz_and -- Logical and. -Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc. +Copyright (C) 1991, 1993, 1994, 1996, 1997, 2000 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -116,7 +117,7 @@ mpz_and (res, op1, op2) _mpz_realloc (res, res_alloc); res_ptr = res->_mp_d; /* Don't re-read OP1_PTR and OP2_PTR. They point to - temporary space--never to the space RES->_mp_D used + temporary space--never to the space RES->_mp_d used to point to before reallocation. */ } @@ -152,9 +153,8 @@ mpz_and (res, op1, op2) { /* We should compute -OP1 & OP2. Swap OP1 and OP2 and fall through to the code that handles OP1 & -OP2. */ - {mpz_srcptr t = op1; op1 = op2; op2 = t;} - {mp_srcptr t = op1_ptr; op1_ptr = op2_ptr; op2_ptr = t;} - {mp_size_t t = op1_size; op1_size = op2_size; op2_size = t;} + MPZ_SRCPTR_SWAP (op1, op2); + MPN_SRCPTR_SWAP (op1_ptr,op1_size, op2_ptr,op2_size); } } @@ -173,7 +173,7 @@ mpz_and (res, op1, op2) operand as the result for those limbs is going to become zero anyway. */ - /* Scan for the least significant. non-zero OP2 limb, and zero the + /* Scan for the least significant non-zero OP2 limb, and zero the result meanwhile for those limb positions. (We will surely find a non-zero limb, so we can write the loop with one termination condition only.) */ @@ -237,7 +237,7 @@ mpz_and (res, op1, op2) res_ptr = res->_mp_d; op1_ptr = op1->_mp_d; /* Don't re-read OP2_PTR. It points to temporary space--never - to the space RES->_mp_D used to point to before reallocation. */ + to the space RES->_mp_d used to point to before reallocation. */ } MPN_COPY (res_ptr + op2_size, op1_ptr + op2_size, @@ -264,7 +264,7 @@ mpz_and (res, op1, op2) res_ptr = res->_mp_d; op1_ptr = op1->_mp_d; /* Don't re-read OP2_PTR. It points to temporary space--never - to the space RES->_mp_D used to point to before reallocation. */ + to the space RES->_mp_d used to point to before reallocation. */ } for (i = res_size - 1; i >= 0; i--) diff --git a/ghc/rts/gmp/mpz/array_init.c b/ghc/rts/gmp/mpz/array_init.c index 8b2e85c90b..1c22046986 100644 --- a/ghc/rts/gmp/mpz/array_init.c +++ b/ghc/rts/gmp/mpz/array_init.c @@ -5,16 +5,16 @@ Copyright (C) 1991, 1993, 1994, 1995 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/bin_ui.c b/ghc/rts/gmp/mpz/bin_ui.c new file mode 100644 index 0000000000..a7a6c98218 --- /dev/null +++ b/ghc/rts/gmp/mpz/bin_ui.c @@ -0,0 +1,141 @@ +/* mpz_bin_uiui - compute n over k. + +Copyright (C) 1998, 1999, 2000 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 2.1 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; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" + + +/* This is a poor implementation. Look at bin_uiui.c for improvement ideas. + In fact consider calling mpz_bin_uiui() when the arguments fit, leaving + the code here only for big n. + + The identity bin(n,k) = (-1)^k * bin(-n+k-1,k) can be found in Knuth vol + 1 section 1.2.6 part G. */ + + +/* Enhancement: use mpn_divexact_1 when it exists */ +#define DIVIDE() \ + ASSERT (SIZ(r) > 0); \ + ASSERT_NOCARRY (mpn_divrem_1 (PTR(r), (mp_size_t) 0, \ + PTR(r), SIZ(r), kacc)); \ + SIZ(r) -= (PTR(r)[SIZ(r)-1] == 0); + +void +#if __STDC__ +mpz_bin_ui (mpz_ptr r, mpz_srcptr n, unsigned long int k) +#else +mpz_bin_ui (r, n, k) + mpz_ptr r; + mpz_srcptr n; + unsigned long int k; +#endif +{ + mpz_t ni; + mp_limb_t i; + mpz_t nacc; + mp_limb_t kacc; + mp_size_t negate; + + if (mpz_sgn (n) < 0) + { + /* bin(n,k) = (-1)^k * bin(-n+k-1,k), and set ni = -n+k-1 - k = -n-1 */ + mpz_init (ni); + mpz_neg (ni, n); + mpz_sub_ui (ni, ni, 1L); + negate = (k & 1); /* (-1)^k */ + } + else + { + /* bin(n,k) == 0 if k>n + (no test for this under the n<0 case, since -n+k-1 >= k there) */ + if (mpz_cmp_ui (n, k) < 0) + { + mpz_set_ui (r, 0L); + return; + } + + /* set ni = n-k */ + mpz_init (ni); + mpz_sub_ui (ni, n, k); + negate = 0; + } + + /* Now wanting bin(ni+k,k), with ni positive, and "negate" is the sign (0 + for positive, 1 for negative). */ + mpz_set_ui (r, 1L); + + /* Rewrite bin(n,k) as bin(n,n-k) if that is smaller. In this case it's + whether ni+k-k < k meaning ni<k, and if so change to denominator ni+k-k + = ni, and new ni of ni+k-ni = k. */ + if (mpz_cmp_ui (ni, k) < 0) + { + unsigned long tmp; + tmp = k; + k = mpz_get_ui (ni); + mpz_set_ui (ni, tmp); + } + + kacc = 1; + mpz_init_set_ui (nacc, 1); + + for (i = 1; i <= k; i++) + { + mp_limb_t k1, k0; + +#if 0 + mp_limb_t nacclow; + int c; + + nacclow = PTR(nacc)[0]; + for (c = 0; (((kacc | nacclow) & 1) == 0); c++) + { + kacc >>= 1; + nacclow >>= 1; + } + mpz_div_2exp (nacc, nacc, c); +#endif + + mpz_add_ui (ni, ni, 1); + mpz_mul (nacc, nacc, ni); + umul_ppmm (k1, k0, kacc, i); + if (k1 != 0) + { + /* Accumulator overflow. Perform bignum step. */ + mpz_mul (r, r, nacc); + mpz_set_ui (nacc, 1); + DIVIDE (); + kacc = i; + } + else + { + /* Save new products in accumulators to keep accumulating. */ + kacc = k0; + } + } + + mpz_mul (r, r, nacc); + DIVIDE (); + SIZ(r) = (SIZ(r) ^ -negate) + negate; + + mpz_clear (nacc); + mpz_clear (ni); +} diff --git a/ghc/rts/gmp/mpz/bin_uiui.c b/ghc/rts/gmp/mpz/bin_uiui.c new file mode 100644 index 0000000000..b37541ba54 --- /dev/null +++ b/ghc/rts/gmp/mpz/bin_uiui.c @@ -0,0 +1,120 @@ +/* mpz_bin_uiui - compute n over k. + +Copyright (C) 1998, 1999, 2000 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 2.1 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; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" + + +/* Avoid reallocs by rounding up any new size */ +#define ROUNDUP_MASK 15 + +/* Enhancement: use mpn_divexact_1 when it exists */ +#define MULDIV() \ + MPZ_REALLOC (r, (SIZ(r)+1)|ROUNDUP_MASK); \ + PTR(r)[SIZ(r)] = mpn_mul_1 (PTR(r), PTR(r), SIZ(r), nacc); \ + ASSERT_NOCARRY (mpn_divrem_1 (PTR(r), (mp_size_t) 0, \ + PTR(r), SIZ(r)+1, kacc)); \ + SIZ(r) += (PTR(r)[SIZ(r)] != 0); + +void +#if __STDC__ +mpz_bin_uiui (mpz_ptr r, unsigned long int n, unsigned long int k) +#else +mpz_bin_uiui (r, n, k) + mpz_ptr r; + unsigned long int n; + unsigned long int k; +#endif +{ + unsigned long int i, j; + mp_limb_t nacc, kacc; + unsigned long int cnt; + + /* bin(n,k) = 0 if k>n. */ + if (n < k) + { + mpz_set_ui (r, 0); + return; + } + + /* Rewrite bin(n,k) as bin(n,n-k) if that is smaller. */ + k = MIN (k, n-k); + + /* bin(n,0) = 1 */ + if (k == 0) + { + mpz_set_ui (r, 1); + return; + } + + j = n - k + 1; + mpz_set_ui (r, j); + + /* Initialize accumulators. */ + nacc = 1; + kacc = 1; + + cnt = 0; + for (i = 2; i <= k; i++) + { + mp_limb_t n1, n0, k1, k0; + + j++; +#if 0 + /* Remove common multiples of 2. This will allow us to accumulate + more in nacc and kacc before we need a bignum step. It would make + sense to cancel factors of 3, 5, etc too, but this would be best + handled by sieving out factors. Alternatively, we could perform a + gcd of the accumulators just as they have overflown, and keep + accumulating until the gcd doesn't remove a significant factor. */ + while (((nacc | kacc) & 1) == 0) + { + nacc >>= 1; + kacc >>= 1; + } +#else + cnt = ((nacc | kacc) & 1) ^ 1; + nacc >>= cnt; + kacc >>= cnt; +#endif + /* Accumulate next multiples. */ + umul_ppmm (n1, n0, nacc, j); + umul_ppmm (k1, k0, kacc, i); + if (n1 != 0) + { + /* Accumulator overflow. Perform bignum step. */ + MULDIV (); + nacc = j; + kacc = i; + } + else + { + if (k1 != 0) abort (); + /* Save new products in accumulators to keep accumulating. */ + nacc = n0; + kacc = k0; + } + } + + /* Take care of whatever is left in accumulators. */ + MULDIV (); +} diff --git a/ghc/rts/gmp/mpz/cdiv_q.c b/ghc/rts/gmp/mpz/cdiv_q.c index 860a232f05..b15ba8aaa9 100644 --- a/ghc/rts/gmp/mpz/cdiv_q.c +++ b/ghc/rts/gmp/mpz/cdiv_q.c @@ -1,21 +1,21 @@ /* mpz_cdiv_q -- Division rounding the quotient towards +infinity. The remainder gets the opposite sign as the denominator. -Copyright (C) 1994, 1995, 1996 Free Software Foundation, Inc. +Copyright (C) 1994, 1995, 1996, 2000 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -40,7 +40,7 @@ mpz_cdiv_q (quot, dividend, divisor) TMP_MARK (marker); - MPZ_TMP_INIT (rem, 1 + ABS (dividend_size)); + MPZ_TMP_INIT (rem, ABS (divisor_size)); mpz_tdiv_qr (quot, rem, dividend, divisor); diff --git a/ghc/rts/gmp/mpz/cdiv_q_ui.c b/ghc/rts/gmp/mpz/cdiv_q_ui.c index 7b6cfd7af8..74f3a90b83 100644 --- a/ghc/rts/gmp/mpz/cdiv_q_ui.c +++ b/ghc/rts/gmp/mpz/cdiv_q_ui.c @@ -3,21 +3,21 @@ always fit into the return type, the negative of the true remainder is returned. -Copyright (C) 1994, 1996 Free Software Foundation, Inc. +Copyright (C) 1994, 1996, 1999 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -40,6 +40,9 @@ mpz_cdiv_q_ui (quot, dividend, divisor) mp_ptr quot_ptr; mp_limb_t remainder_limb; + if (divisor == 0) + DIVIDE_BY_ZERO; + dividend_size = dividend->_mp_size; size = ABS (dividend_size); @@ -53,7 +56,7 @@ mpz_cdiv_q_ui (quot, dividend, divisor) if (remainder_limb != 0 && dividend_size >= 0) { - mpn_add_1 (quot_ptr, quot_ptr, size, (mp_limb_t) 1); + mpn_incr_u (quot_ptr, (mp_limb_t) 1); remainder_limb = divisor - remainder_limb; } diff --git a/ghc/rts/gmp/mpz/cdiv_qr.c b/ghc/rts/gmp/mpz/cdiv_qr.c index bf7d6dacb1..29c7c41a4e 100644 --- a/ghc/rts/gmp/mpz/cdiv_qr.c +++ b/ghc/rts/gmp/mpz/cdiv_qr.c @@ -1,21 +1,21 @@ /* mpz_cdiv_qr -- Division rounding the quotient towards +infinity. The remainder gets the opposite sign as the denominator. -Copyright (C) 1994, 1995, 1996 Free Software Foundation, Inc. +Copyright (C) 1994, 1995, 1996, 2000 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -35,6 +35,7 @@ mpz_cdiv_qr (quot, rem, dividend, divisor) #endif { mp_size_t divisor_size = divisor->_mp_size; + mp_size_t xsize; mpz_t temp_divisor; /* N.B.: lives until function returns! */ TMP_DECL (marker); @@ -50,9 +51,10 @@ mpz_cdiv_qr (quot, rem, dividend, divisor) divisor = temp_divisor; } + xsize = dividend->_mp_size ^ divisor_size;; mpz_tdiv_qr (quot, rem, dividend, divisor); - if ((divisor_size ^ dividend->_mp_size) >= 0 && rem->_mp_size != 0) + if (xsize >= 0 && rem->_mp_size != 0) { mpz_add_ui (quot, quot, 1L); mpz_sub (rem, rem, divisor); diff --git a/ghc/rts/gmp/mpz/cdiv_qr_ui.c b/ghc/rts/gmp/mpz/cdiv_qr_ui.c index a780e772ec..a7873c6e20 100644 --- a/ghc/rts/gmp/mpz/cdiv_qr_ui.c +++ b/ghc/rts/gmp/mpz/cdiv_qr_ui.c @@ -3,21 +3,21 @@ always fit into the return type, the negative of the true remainder is returned. -Copyright (C) 1994, 1995, 1996 Free Software Foundation, Inc. +Copyright (C) 1994, 1995, 1996, 1999 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -41,6 +41,9 @@ mpz_cdiv_qr_ui (quot, rem, dividend, divisor) mp_ptr quot_ptr; mp_limb_t remainder_limb; + if (divisor == 0) + DIVIDE_BY_ZERO; + dividend_size = dividend->_mp_size; size = ABS (dividend_size); @@ -50,11 +53,11 @@ mpz_cdiv_qr_ui (quot, rem, dividend, divisor) quot_ptr = quot->_mp_d; remainder_limb = mpn_divmod_1 (quot_ptr, dividend->_mp_d, size, - (mp_limb_t) divisor); + (mp_limb_t) divisor); if (remainder_limb != 0 && dividend_size >= 0) { - mpn_add_1 (quot_ptr, quot_ptr, size, (mp_limb_t) 1); + mpn_incr_u (quot_ptr, (mp_limb_t) 1); remainder_limb = divisor - remainder_limb; } diff --git a/ghc/rts/gmp/mpz/cdiv_r.c b/ghc/rts/gmp/mpz/cdiv_r.c index d34d138ce7..e96ce7e677 100644 --- a/ghc/rts/gmp/mpz/cdiv_r.c +++ b/ghc/rts/gmp/mpz/cdiv_r.c @@ -6,16 +6,16 @@ Copyright (C) 1994, 1995, 1996 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/cdiv_r_ui.c b/ghc/rts/gmp/mpz/cdiv_r_ui.c index 757a3f5f7b..e17e2381c0 100644 --- a/ghc/rts/gmp/mpz/cdiv_r_ui.c +++ b/ghc/rts/gmp/mpz/cdiv_r_ui.c @@ -8,16 +8,16 @@ Copyright (C) 1994, 1995, 1996 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -39,6 +39,9 @@ mpz_cdiv_r_ui (rem, dividend, divisor) mp_size_t size; mp_limb_t remainder_limb; + if (divisor == 0) + DIVIDE_BY_ZERO; + dividend_size = dividend->_mp_size; size = ABS (dividend_size); diff --git a/ghc/rts/gmp/mpz/cdiv_ui.c b/ghc/rts/gmp/mpz/cdiv_ui.c index df841ede88..63547a78c0 100644 --- a/ghc/rts/gmp/mpz/cdiv_ui.c +++ b/ghc/rts/gmp/mpz/cdiv_ui.c @@ -8,16 +8,16 @@ Copyright (C) 1994, 1995, 1996 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/clear.c b/ghc/rts/gmp/mpz/clear.c index 00f3cfdc62..5224553f9e 100644 --- a/ghc/rts/gmp/mpz/clear.c +++ b/ghc/rts/gmp/mpz/clear.c @@ -6,16 +6,16 @@ Copyright (C) 1991, 1993, 1994, 1995 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/clrbit.c b/ghc/rts/gmp/mpz/clrbit.c index 59d9565358..865d84902f 100644 --- a/ghc/rts/gmp/mpz/clrbit.c +++ b/ghc/rts/gmp/mpz/clrbit.c @@ -5,16 +5,16 @@ Copyright (C) 1991, 1993, 1994, 1995 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/cmp.c b/ghc/rts/gmp/mpz/cmp.c index 37be3349a6..60628348e5 100644 --- a/ghc/rts/gmp/mpz/cmp.c +++ b/ghc/rts/gmp/mpz/cmp.c @@ -6,16 +6,16 @@ Copyright (C) 1991, 1993, 1994, 1996 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/cmp_si.c b/ghc/rts/gmp/mpz/cmp_si.c index 8063fd386e..0c2212fbe9 100644 --- a/ghc/rts/gmp/mpz/cmp_si.c +++ b/ghc/rts/gmp/mpz/cmp_si.c @@ -1,21 +1,22 @@ /* mpz_cmp_si(u,v) -- Compare an integer U with a single-word int V. Return positive, zero, or negative based on if U > V, U == V, or U < V. -Copyright (C) 1991, 1993, 1994, 1995 Free Software Foundation, Inc. +Copyright (C) 1991, 1993, 1994, 1995, 1996, 2000 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -23,14 +24,11 @@ MA 02111-1307, USA. */ #include "gmp.h" #include "gmp-impl.h" -/* gmp.h defines a macro for mpz_cmp_si. */ -#undef mpz_cmp_si - int #if __STDC__ -mpz_cmp_si (mpz_srcptr u, signed long int v_digit) +_mpz_cmp_si (mpz_srcptr u, signed long int v_digit) #else -mpz_cmp_si (u, v_digit) +_mpz_cmp_si (u, v_digit) mpz_srcptr u; signed long int v_digit; #endif @@ -56,10 +54,10 @@ mpz_cmp_si (u, v_digit) u_digit = u->_mp_d[0]; - if (u_digit == v_digit) + if (u_digit == (mp_limb_t) (unsigned long) v_digit) return 0; - if (u_digit > v_digit) + if (u_digit > (mp_limb_t) (unsigned long) v_digit) return usize; else return -usize; diff --git a/ghc/rts/gmp/mpz/cmp_ui.c b/ghc/rts/gmp/mpz/cmp_ui.c index 1a50b96e1a..fd84f301c1 100644 --- a/ghc/rts/gmp/mpz/cmp_ui.c +++ b/ghc/rts/gmp/mpz/cmp_ui.c @@ -1,21 +1,21 @@ /* mpz_cmp_ui.c -- Compare a mpz_t a with an mp_limb_t b. Return positive, zero, or negative based on if a > b, a == b, or a < b. -Copyright (C) 1991, 1993, 1994, 1995 Free Software Foundation, Inc. +Copyright (C) 1991, 1993, 1994, 1995, 1996 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -23,14 +23,11 @@ MA 02111-1307, USA. */ #include "gmp.h" #include "gmp-impl.h" -/* gmp.h defines a macro for mpz_cmp_ui. */ -#undef mpz_cmp_ui - int #if __STDC__ -mpz_cmp_ui (mpz_srcptr u, unsigned long int v_digit) +_mpz_cmp_ui (mpz_srcptr u, unsigned long int v_digit) #else -mpz_cmp_ui (u, v_digit) +_mpz_cmp_ui (u, v_digit) mpz_srcptr u; unsigned long int v_digit; #endif diff --git a/ghc/rts/gmp/mpz/cmpabs.c b/ghc/rts/gmp/mpz/cmpabs.c new file mode 100644 index 0000000000..037d7a9145 --- /dev/null +++ b/ghc/rts/gmp/mpz/cmpabs.c @@ -0,0 +1,57 @@ +/* mpz_cmpabs(u,v) -- Compare U, V. Return postive, zero, or negative + based on if U > V, U == V, or U < V. + +Copyright (C) 1991, 1993, 1994, 1996, 1997, 2000 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 2.1 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; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +int +#if __STDC__ +mpz_cmpabs (mpz_srcptr u, mpz_srcptr v) +#else +mpz_cmpabs (u, v) + mpz_srcptr u; + mpz_srcptr v; +#endif +{ + mp_size_t usize = u->_mp_size; + mp_size_t vsize = v->_mp_size; + mp_size_t size; + mp_srcptr up, vp; + int cmp; + + usize = ABS (usize); + vsize = ABS (vsize); + + if (usize != vsize) + return usize - vsize; + + if (usize == 0) + return 0; + + up = u->_mp_d; + vp = v->_mp_d; + + cmp = mpn_cmp (up, vp, usize); + + return cmp; +} diff --git a/ghc/rts/gmp/mpz/cmpabs_ui.c b/ghc/rts/gmp/mpz/cmpabs_ui.c new file mode 100644 index 0000000000..db816b5820 --- /dev/null +++ b/ghc/rts/gmp/mpz/cmpabs_ui.c @@ -0,0 +1,56 @@ +/* mpz_cmpabs_ui.c -- Compare a mpz_t a with an mp_limb_t b. Return positive, + zero, or negative based on if a > b, a == b, or a < b. + +Copyright (C) 1991, 1993, 1994, 1995, 1997, 2000 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 2.1 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; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +int +#if __STDC__ +mpz_cmpabs_ui (mpz_srcptr u, unsigned long int v_digit) +#else +mpz_cmpabs_ui (u, v_digit) + mpz_srcptr u; + unsigned long int v_digit; +#endif +{ + mp_size_t usize = u->_mp_size; + + if (usize == 0) + return -(v_digit != 0); + + usize = ABS (usize); + + if (usize == 1) + { + mp_limb_t u_digit; + + u_digit = u->_mp_d[0]; + if (u_digit > v_digit) + return 1; + if (u_digit < v_digit) + return -1; + return 0; + } + + return 1; +} diff --git a/ghc/rts/gmp/mpz/com.c b/ghc/rts/gmp/mpz/com.c index 559f2b6b2b..18d6427779 100644 --- a/ghc/rts/gmp/mpz/com.c +++ b/ghc/rts/gmp/mpz/com.c @@ -6,16 +6,16 @@ Copyright (C) 1991, 1993, 1994, 1996 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/divexact.c b/ghc/rts/gmp/mpz/divexact.c index b897448046..c2970454fd 100644 --- a/ghc/rts/gmp/mpz/divexact.c +++ b/ghc/rts/gmp/mpz/divexact.c @@ -1,20 +1,21 @@ /* mpz_divexact -- finds quotient when known that quot * den == num && den != 0. -Copyright (C) 1991, 1993, 1994, 1995 Free Software Foundation, Inc. +Copyright (C) 1991, 1993, 1994, 1995, 1996, 1997, 1998, 2000 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -27,8 +28,8 @@ MA 02111-1307, USA. */ de Matema'tica at Universidade Federal do Rio Grande do Sul (UFRGS). References: - T. Jebelean, An algorithm for exact division, Journal of Symbolic - Computation, v. 15, 1993, pp. 169-180. */ + T. Jebelean, An algorithm for exact division, Journal of Symbolic + Computation, v. 15, 1993, pp. 169-180. */ #include "gmp.h" #include "gmp-impl.h" @@ -46,30 +47,42 @@ mpz_divexact (quot, num, den) { mp_ptr qp, tp; mp_size_t qsize, tsize; - - mp_srcptr np = num->_mp_d; - mp_srcptr dp = den->_mp_d; - mp_size_t nsize = ABS (num->_mp_size); - mp_size_t dsize = ABS (den->_mp_size); + mp_srcptr np, dp; + mp_size_t nsize, dsize; TMP_DECL (marker); - /* Generate divide-by-zero error if dsize == 0. */ - if (dsize == 0) - { - quot->_mp_size = 1 / dsize; - return; - } + nsize = ABS (num->_mp_size); + dsize = ABS (den->_mp_size); + + qsize = nsize - dsize + 1; + if (quot->_mp_alloc < qsize) + _mpz_realloc (quot, qsize); + + np = num->_mp_d; + dp = den->_mp_d; + qp = quot->_mp_d; if (nsize == 0) { + if (dsize == 0) + DIVIDE_BY_ZERO; quot->_mp_size = 0; return; } - qsize = nsize - dsize + 1; - if (quot->_mp_alloc < qsize) - _mpz_realloc (quot, qsize); - qp = quot->_mp_d; + if (dsize <= 1) + { + if (dsize == 1) + { + mpn_divmod_1 (qp, np, nsize, dp[0]); + qsize -= qp[qsize - 1] == 0; + quot->_mp_size = (num->_mp_size ^ den->_mp_size) >= 0 ? qsize : -qsize; + return; + } + + /* Generate divide-by-zero error since dsize == 0. */ + DIVIDE_BY_ZERO; + } TMP_MARK (marker); @@ -77,36 +90,36 @@ mpz_divexact (quot, num, den) while (dp[0] == 0) np += 1, nsize -= 1, dp += 1, dsize -= 1; tsize = MIN (qsize, dsize); - if (dp[0] & 1) + if ((dp[0] & 1) != 0) { - if (qp != dp) - MPN_COPY (qp, np, qsize); - if (qp == dp) /* QUOT and DEN overlap. */ + if (quot == den) /* QUOT and DEN overlap. */ { - tp = (mp_ptr) TMP_ALLOC (sizeof (mp_limb_t) * tsize); + tp = (mp_ptr) TMP_ALLOC (tsize * BYTES_PER_MP_LIMB); MPN_COPY (tp, dp, tsize); } else tp = (mp_ptr) dp; + if (qp != np) + MPN_COPY_INCR (qp, np, qsize); } else { - unsigned long int r; - tp = (mp_ptr) TMP_ALLOC (sizeof (mp_limb_t) * tsize); + unsigned int r; + tp = (mp_ptr) TMP_ALLOC (tsize * BYTES_PER_MP_LIMB); count_trailing_zeros (r, dp[0]); mpn_rshift (tp, dp, tsize, r); if (dsize > tsize) - tp[tsize-1] |= dp[tsize] << (BITS_PER_MP_LIMB - r); + tp[tsize - 1] |= dp[tsize] << (BITS_PER_MP_LIMB - r); mpn_rshift (qp, np, qsize, r); if (nsize > qsize) - qp[qsize-1] |= np[qsize] << (BITS_PER_MP_LIMB - r); + qp[qsize - 1] |= np[qsize] << (BITS_PER_MP_LIMB - r); } /* Now QUOT <-- QUOT/T. */ mpn_bdivmod (qp, qp, qsize, tp, tsize, qsize * BITS_PER_MP_LIMB); MPN_NORMALIZE (qp, qsize); - quot->_mp_size = (num->_mp_size < 0) == (den->_mp_size < 0) ? qsize : -qsize; + quot->_mp_size = (num->_mp_size ^ den->_mp_size) >= 0 ? qsize : -qsize; TMP_FREE (marker); } diff --git a/ghc/rts/gmp/mpz/dump.c b/ghc/rts/gmp/mpz/dump.c new file mode 100644 index 0000000000..dc318ac8cf --- /dev/null +++ b/ghc/rts/gmp/mpz/dump.c @@ -0,0 +1,44 @@ +/* mpz_dump - Dump an integer to stdout. + + THIS IS AN INTERNAL FUNCTION WITH A MUTABLE INTERFACE. IT IS NOT SAFE TO + CALL THIS FUNCTION DIRECTLY. IN FACT, IT IS ALMOST GUARANTEED THAT THIS + FUNCTION WILL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE. + + +Copyright (C) 1999, 2000 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 2.1 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; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include <stdio.h> +#include "gmp.h" +#include "gmp-impl.h" + +void +#if __STDC__ +mpz_dump (mpz_srcptr u) +#else +mpz_dump (u) + mpz_srcptr u; +#endif +{ + char *str; + + str = mpz_get_str (0, 10, u); + printf ("%s\n", str); + (*_mp_free_func) (str, 0);/* ??? broken alloc interface, pass what size ??? */ +} diff --git a/ghc/rts/gmp/mpz/fac_ui.c b/ghc/rts/gmp/mpz/fac_ui.c index a170060a54..85f40f271c 100644 --- a/ghc/rts/gmp/mpz/fac_ui.c +++ b/ghc/rts/gmp/mpz/fac_ui.c @@ -5,16 +5,16 @@ Copyright (C) 1991, 1993, 1994, 1995 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/fdiv_q.c b/ghc/rts/gmp/mpz/fdiv_q.c index 3da943af91..9d75ca33d2 100644 --- a/ghc/rts/gmp/mpz/fdiv_q.c +++ b/ghc/rts/gmp/mpz/fdiv_q.c @@ -1,21 +1,21 @@ /* mpz_fdiv_q -- Division rounding the quotient towards -infinity. The remainder gets the same sign as the denominator. -Copyright (C) 1994, 1995, 1996 Free Software Foundation, Inc. +Copyright (C) 1994, 1995, 1996, 2000 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -40,7 +40,7 @@ mpz_fdiv_q (quot, dividend, divisor) TMP_MARK (marker); - MPZ_TMP_INIT (rem, 1 + ABS (dividend_size)); + MPZ_TMP_INIT (rem, ABS (divisor_size)); mpz_tdiv_qr (quot, rem, dividend, divisor); diff --git a/ghc/rts/gmp/mpz/fdiv_q_2exp.c b/ghc/rts/gmp/mpz/fdiv_q_2exp.c index 3f56baff2c..8e02180ecc 100644 --- a/ghc/rts/gmp/mpz/fdiv_q_2exp.c +++ b/ghc/rts/gmp/mpz/fdiv_q_2exp.c @@ -1,21 +1,22 @@ /* mpz_fdiv_q_2exp -- Divide an integer by 2**CNT. Round the quotient towards -infinity. -Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc. +Copyright (C) 1991, 1993, 1994, 1996, 1998, 1999 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -86,9 +87,18 @@ mpz_fdiv_q_2exp (w, u, cnt) if (usize < 0 && round != 0) { mp_limb_t cy; - cy = mpn_add_1 (wp, wp, wsize, 1); - wp[wsize] = cy; - wsize += cy; + if (wsize != 0) + { + cy = mpn_add_1 (wp, wp, wsize, (mp_limb_t) 1); + wp[wsize] = cy; + wsize += cy; + } + else + { + /* We shifted something negative to zero. The result is -1. */ + wp[0] = 1; + wsize = 1; + } } w->_mp_size = usize >= 0 ? wsize : -wsize; } diff --git a/ghc/rts/gmp/mpz/fdiv_q_ui.c b/ghc/rts/gmp/mpz/fdiv_q_ui.c index 3d6825d4fe..55d2498693 100644 --- a/ghc/rts/gmp/mpz/fdiv_q_ui.c +++ b/ghc/rts/gmp/mpz/fdiv_q_ui.c @@ -1,21 +1,21 @@ /* mpz_fdiv_q_ui -- Division rounding the quotient towards -infinity. The remainder gets the same sign as the denominator. -Copyright (C) 1994, 1995, 1996 Free Software Foundation, Inc. +Copyright (C) 1994, 1995, 1996, 1999 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -38,6 +38,9 @@ mpz_fdiv_q_ui (quot, dividend, divisor) mp_ptr quot_ptr; mp_limb_t remainder_limb; + if (divisor == 0) + DIVIDE_BY_ZERO; + dividend_size = dividend->_mp_size; size = ABS (dividend_size); @@ -51,7 +54,7 @@ mpz_fdiv_q_ui (quot, dividend, divisor) if (remainder_limb != 0 && dividend_size < 0) { - mpn_add_1 (quot_ptr, quot_ptr, size, (mp_limb_t) 1); + mpn_incr_u (quot_ptr, (mp_limb_t) 1); remainder_limb = divisor - remainder_limb; } diff --git a/ghc/rts/gmp/mpz/fdiv_qr.c b/ghc/rts/gmp/mpz/fdiv_qr.c index 2abb16ccbd..06ce50607b 100644 --- a/ghc/rts/gmp/mpz/fdiv_qr.c +++ b/ghc/rts/gmp/mpz/fdiv_qr.c @@ -1,21 +1,21 @@ /* mpz_fdiv_qr -- Division rounding the quotient towards -infinity. The remainder gets the same sign as the denominator. -Copyright (C) 1994, 1995, 1996 Free Software Foundation, Inc. +Copyright (C) 1994, 1995, 1996, 2000 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -35,6 +35,7 @@ mpz_fdiv_qr (quot, rem, dividend, divisor) #endif { mp_size_t divisor_size = divisor->_mp_size; + mp_size_t xsize; mpz_t temp_divisor; /* N.B.: lives until function returns! */ TMP_DECL (marker); @@ -50,9 +51,10 @@ mpz_fdiv_qr (quot, rem, dividend, divisor) divisor = temp_divisor; } + xsize = dividend->_mp_size ^ divisor_size;; mpz_tdiv_qr (quot, rem, dividend, divisor); - if ((divisor_size ^ dividend->_mp_size) < 0 && rem->_mp_size != 0) + if (xsize < 0 && rem->_mp_size != 0) { mpz_sub_ui (quot, quot, 1L); mpz_add (rem, rem, divisor); diff --git a/ghc/rts/gmp/mpz/fdiv_qr_ui.c b/ghc/rts/gmp/mpz/fdiv_qr_ui.c index a22b702a8d..600c0dacfc 100644 --- a/ghc/rts/gmp/mpz/fdiv_qr_ui.c +++ b/ghc/rts/gmp/mpz/fdiv_qr_ui.c @@ -1,21 +1,21 @@ /* mpz_fdiv_qr_ui -- Division rounding the quotient towards -infinity. The remainder gets the same sign as the denominator. -Copyright (C) 1994, 1995, 1996 Free Software Foundation, Inc. +Copyright (C) 1994, 1995, 1996, 1999 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -39,6 +39,9 @@ mpz_fdiv_qr_ui (quot, rem, dividend, divisor) mp_ptr quot_ptr; mp_limb_t remainder_limb; + if (divisor == 0) + DIVIDE_BY_ZERO; + dividend_size = dividend->_mp_size; size = ABS (dividend_size); @@ -48,11 +51,11 @@ mpz_fdiv_qr_ui (quot, rem, dividend, divisor) quot_ptr = quot->_mp_d; remainder_limb = mpn_divmod_1 (quot_ptr, dividend->_mp_d, size, - (mp_limb_t) divisor); + (mp_limb_t) divisor); if (remainder_limb != 0 && dividend_size < 0) { - mpn_add_1 (quot_ptr, quot_ptr, size, (mp_limb_t) 1); + mpn_incr_u (quot_ptr, (mp_limb_t) 1); remainder_limb = divisor - remainder_limb; } diff --git a/ghc/rts/gmp/mpz/fdiv_r.c b/ghc/rts/gmp/mpz/fdiv_r.c index 14e045bbe5..a3652838d2 100644 --- a/ghc/rts/gmp/mpz/fdiv_r.c +++ b/ghc/rts/gmp/mpz/fdiv_r.c @@ -6,16 +6,16 @@ Copyright (C) 1994, 1995, 1996 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/fdiv_r_2exp.c b/ghc/rts/gmp/mpz/fdiv_r_2exp.c index 04190b1b0a..d0cd03945f 100644 --- a/ghc/rts/gmp/mpz/fdiv_r_2exp.c +++ b/ghc/rts/gmp/mpz/fdiv_r_2exp.c @@ -1,20 +1,21 @@ /* mpz_fdiv_r_2exp -- Divide a integer by 2**CNT and produce a remainder. -Copyright (C) 1991, 1993, 1994, 1995 Free Software Foundation, Inc. +Copyright (C) 1991, 1993, 1994, 1995, 1998, 1999 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -65,7 +66,8 @@ mpz_fdiv_r_2exp (res, in, cnt) else { /* The input operand is smaller than 2**CNT. We perform a no-op, - apart from that we might need to copy IN to RES. */ + apart from that we might need to copy IN to RES, and may need + to round the result. */ res_size = in_size; if (res->_mp_alloc < res_size) _mpz_realloc (res, res_size); @@ -75,14 +77,77 @@ mpz_fdiv_r_2exp (res, in, cnt) if (res != in) MPN_COPY (res->_mp_d, in->_mp_d, limb_cnt); + in_size = in->_mp_size; res->_mp_size = res_size; - if (in->_mp_size < 0 && res_size != 0) + if (in_size < 0 && res_size != 0) { /* Result should be 2^CNT - RES */ mpz_t tmp; - MPZ_TMP_INIT (tmp, limb_cnt + 1); + MPZ_TMP_INIT (tmp, cnt/BITS_PER_MP_LIMB + 2); mpz_set_ui (tmp, 1L); mpz_mul_2exp (tmp, tmp, cnt); mpz_sub (res, tmp, res); } } + +/* This is an alternative ending of the above function using just low-level + functions. Tested, but perhaps excessive? */ +#if 0 + if (in->_mp_size < 0 && res_size != 0) + { + /* Result should be 2^CNT - RES */ + + mp_ptr rp; + + limb_cnt = cnt / BITS_PER_MP_LIMB; + + if (res->_mp_alloc <= limb_cnt) + _mpz_realloc (res, limb_cnt + 1); + rp = PTR(res); + if (res_size > limb_cnt) + { + mpn_nz_neg (rp, rp, res_size); + rp[limb_cnt] &= ~(~(mp_limb_t) 0 << cnt % BITS_PER_MP_LIMB); + MPN_NORMALIZE_NOT_ZERO (rp, res_size); + } + else + { + mp_size_t i; + mpn_nz_neg (rp, rp, res_size); + for (i = res_size; i < limb_cnt; i++) + rp[i] = ~ (mp_limb_t) 0; + res_size = limb_cnt; + if (cnt % BITS_PER_MP_LIMB != 0) + { + rp[res_size] = ((mp_limb_t) 1 << (cnt % BITS_PER_MP_LIMB)) - 1; + res_size++; + } + else + MPN_NORMALIZE_NOT_ZERO (rp, res_size); + } + } + SIZ(res) = res_size; +} + +static void +mpn_nz_neg (rp, sp, n) + mp_ptr rp, sp; + mp_size_t n; +{ + mp_size_t i; + mp_limb_t x; + + x = sp[0]; + rp[0] = -x; + for (i = 1; x == 0; i++) + { + x = sp[i]; + rp[i] = -x; + } + + for (; i < n; i++) + { + rp[i] = ~sp[i]; + } +} +#endif diff --git a/ghc/rts/gmp/mpz/fdiv_r_ui.c b/ghc/rts/gmp/mpz/fdiv_r_ui.c index c4c37498ce..dd5c743d27 100644 --- a/ghc/rts/gmp/mpz/fdiv_r_ui.c +++ b/ghc/rts/gmp/mpz/fdiv_r_ui.c @@ -6,16 +6,16 @@ Copyright (C) 1994, 1995, 1996 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -37,6 +37,9 @@ mpz_fdiv_r_ui (rem, dividend, divisor) mp_size_t size; mp_limb_t remainder_limb; + if (divisor == 0) + DIVIDE_BY_ZERO; + dividend_size = dividend->_mp_size; size = ABS (dividend_size); diff --git a/ghc/rts/gmp/mpz/fdiv_ui.c b/ghc/rts/gmp/mpz/fdiv_ui.c index 4d018a2291..f937b5f6d0 100644 --- a/ghc/rts/gmp/mpz/fdiv_ui.c +++ b/ghc/rts/gmp/mpz/fdiv_ui.c @@ -6,16 +6,16 @@ Copyright (C) 1994, 1995, 1996 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/fib_ui.c b/ghc/rts/gmp/mpz/fib_ui.c new file mode 100644 index 0000000000..4bebb80d94 --- /dev/null +++ b/ghc/rts/gmp/mpz/fib_ui.c @@ -0,0 +1,165 @@ +/* mpz_fib_ui(result, n) -- Set RESULT to the Nth Fibonacci number. + +Copyright (C) 1998, 1999, 2000 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 2.1 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; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +/* This is fast, but could be made somewhat faster and neater. + The timing is somewhat fluctuating for even/odd sizes because + of the extra hair used to save variables and operations. Here + are a few things one might want to address: + 1. Avoid using 4 intermediate variables in mpz_fib_bigcase. + 2. Call mpn functions directly. Straightforward for these functions. + 3. Merge the three functions into one. + +Said by Kevin: + Consider using the Lucas numbers L[n] as an auxiliary sequence, making + it possible to do the "doubling" operation in mpz_fib_bigcase with two + squares rather than two multiplies. The formulas are a little more + complicated, something like the following (untested). + + F[2n] = ((F[n]+L[n])^2 - 6*F[n]^2 - 4*(-1)^n) / 2 + L[2n] = 5*F[n]^2 + 2*(-1)^n + + F[2n+1] = (F[2n] + L[2n]) / 2 + L[2n+1] = (5*F[2n] + L[2n]) / 2 + + The Lucas number that comes for free here could even be returned. + + Maybe there's formulas with two squares using just F[n], but I don't + know of any. +*/ + +/* Determine the needed storage for Fib(n). */ +#define FIB_SIZE(n) (((mp_size_t) ((n)*0.695)) / BITS_PER_MP_LIMB + 2) + +static void mpz_fib_bigcase _PROTO ((mpz_t, mpz_t, unsigned long int)); +static void mpz_fib_basecase _PROTO ((mpz_t, mpz_t, unsigned long int)); + + +#ifndef FIB_THRESHOLD +#define FIB_THRESHOLD 60 +#endif + +void +#if __STDC__ +mpz_fib_ui (mpz_t r, unsigned long int n) +#else +mpz_fib_ui (r, n) + mpz_t r; + unsigned long int n; +#endif +{ + if (n == 0) + mpz_set_ui (r, 0); + else + { + mpz_t t1; + mpz_init (t1); + if (n < FIB_THRESHOLD) + mpz_fib_basecase (t1, r, n); + else + mpz_fib_bigcase (t1, r, n); + mpz_clear (t1); + } +} + +static void +#if __STDC__ +mpz_fib_basecase (mpz_t t1, mpz_t t2, unsigned long int n) +#else +mpz_fib_basecase (t1, t2, n) + mpz_t t1; + mpz_t t2; + unsigned long int n; +#endif +{ + unsigned long int m, i; + + mpz_set_ui (t1, 0); + mpz_set_ui (t2, 1); + m = n/2; + for (i = 0; i < m; i++) + { + mpz_add (t1, t1, t2); + mpz_add (t2, t1, t2); + } + if ((n & 1) == 0) + { + mpz_sub (t1, t2, t1); + mpz_sub (t2, t2, t1); /* trick: recover t1 value just overwritten */ + } +} + +static void +#if __STDC__ +mpz_fib_bigcase (mpz_t t1, mpz_t t2, unsigned long int n) +#else +mpz_fib_bigcase (t1, t2, n) + mpz_t t1; + mpz_t t2; + unsigned long int n; +#endif +{ + unsigned long int n2; + int ni, i; + mpz_t x1, x2, u1, u2; + + ni = 0; + for (n2 = n; n2 >= FIB_THRESHOLD; n2 /= 2) + ni++; + + mpz_fib_basecase (t1, t2, n2); + + mpz_init (x1); + mpz_init (x2); + mpz_init (u1); + mpz_init (u2); + + for (i = ni - 1; i >= 0; i--) + { + mpz_mul_2exp (x1, t1, 1); + mpz_mul_2exp (x2, t2, 1); + + mpz_add (x1, x1, t2); + mpz_sub (x2, x2, t1); + + mpz_mul (u1, t2, x1); + mpz_mul (u2, t1, x2); + + if (((n >> i) & 1) == 0) + { + mpz_sub (t1, u1, u2); + mpz_set (t2, u1); + } + else + { + mpz_set (t1, u1); + mpz_mul_2exp (t2, u1, 1); + mpz_sub (t2, t2, u2); + } + } + + mpz_clear (x1); + mpz_clear (x2); + mpz_clear (u1); + mpz_clear (u2); +} diff --git a/ghc/rts/gmp/mpz/fits_sint_p.c b/ghc/rts/gmp/mpz/fits_sint_p.c new file mode 100644 index 0000000000..82e32a24d5 --- /dev/null +++ b/ghc/rts/gmp/mpz/fits_sint_p.c @@ -0,0 +1,50 @@ +/* int mpz_fits_X_p (mpz_t src) -- Return whether src fits the C type X. + +Copyright (C) 1997, 2000 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 2.1 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; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +int +#if __STDC__ +mpz_fits_sint_p (mpz_srcptr src) +#else +mpz_fits_sint_p (src) + mpz_srcptr src; +#endif +{ + mp_size_t size; + mp_limb_t mpl; + + mpl = PTR(src)[0]; + size = SIZ(src); + if (size > 0) + { + if (size > 1) + return 0; + return mpl < ~((~(unsigned int) 0) >> 1); + } + else + { + if (size < -1) + return 0; + return mpl <= ~((~(unsigned int) 0) >> 1); + } +} diff --git a/ghc/rts/gmp/mpz/fits_slong_p.c b/ghc/rts/gmp/mpz/fits_slong_p.c new file mode 100644 index 0000000000..e0669b5aaa --- /dev/null +++ b/ghc/rts/gmp/mpz/fits_slong_p.c @@ -0,0 +1,50 @@ +/* int mpz_fits_X_p (mpz_t src) -- Return whether src fits the C type X. + +Copyright (C) 1997, 2000 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 2.1 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; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +int +#if __STDC__ +mpz_fits_slong_p (mpz_srcptr src) +#else +mpz_fits_slong_p (src) + mpz_srcptr src; +#endif +{ + mp_size_t size; + mp_limb_t mpl; + + mpl = PTR(src)[0]; + size = SIZ(src); + if (size > 0) + { + if (size > 1) + return 0; + return mpl < ~((~(unsigned long int) 0) >> 1); + } + else + { + if (size < -1) + return 0; + return mpl <= ~((~(unsigned long int) 0) >> 1); + } +} diff --git a/ghc/rts/gmp/mpz/fits_sshort_p.c b/ghc/rts/gmp/mpz/fits_sshort_p.c new file mode 100644 index 0000000000..5b8e31afae --- /dev/null +++ b/ghc/rts/gmp/mpz/fits_sshort_p.c @@ -0,0 +1,50 @@ +/* int mpz_fits_X_p (mpz_t src) -- Return whether src fits the C type X. + +Copyright (C) 1997, 2000 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 2.1 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; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +int +#if __STDC__ +mpz_fits_sshort_p (mpz_srcptr src) +#else +mpz_fits_sshort_p (src) + mpz_srcptr src; +#endif +{ + mp_size_t size; + mp_limb_t mpl; + + mpl = PTR(src)[0]; + size = SIZ(src); + if (size > 0) + { + if (size > 1) + return 0; + return mpl <= (((unsigned short int) ~(unsigned int) 0) >> 1); + } + else + { + if (size < -1) + return 0; + return mpl <= (((unsigned short int) ~(unsigned int) 0) >> 1) + 1; + } +} diff --git a/ghc/rts/gmp/mpz/fits_uint_p.c b/ghc/rts/gmp/mpz/fits_uint_p.c new file mode 100644 index 0000000000..72f62fa723 --- /dev/null +++ b/ghc/rts/gmp/mpz/fits_uint_p.c @@ -0,0 +1,41 @@ +/* int mpz_fits_X_p (mpz_t src) -- Return whether src fits the C type X. + +Copyright (C) 1997, 2000 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 2.1 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; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +int +#if __STDC__ +mpz_fits_uint_p (mpz_srcptr src) +#else +mpz_fits_uint_p (src) + mpz_srcptr src; +#endif +{ + mp_size_t size; + mp_limb_t mpl; + + mpl = PTR(src)[0]; + size = SIZ(src); + if (size < 0 || size > 1) + return 0; + return mpl <= (~(unsigned int) 0); +} diff --git a/ghc/rts/gmp/mpz/fits_ulong_p.c b/ghc/rts/gmp/mpz/fits_ulong_p.c new file mode 100644 index 0000000000..92eb42e86e --- /dev/null +++ b/ghc/rts/gmp/mpz/fits_ulong_p.c @@ -0,0 +1,41 @@ +/* int mpz_fits_X_p (mpz_t src) -- Return whether src fits the C type X. + +Copyright (C) 1997, 2000 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 2.1 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; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +int +#if __STDC__ +mpz_fits_ulong_p (mpz_srcptr src) +#else +mpz_fits_ulong_p (src) + mpz_srcptr src; +#endif +{ + mp_size_t size; + mp_limb_t mpl; + + mpl = PTR(src)[0]; + size = SIZ(src); + if (size < 0 || size > 1) + return 0; + return mpl <= (~(unsigned long int) 0); +} diff --git a/ghc/rts/gmp/mpz/fits_ushort_p.c b/ghc/rts/gmp/mpz/fits_ushort_p.c new file mode 100644 index 0000000000..bde0edae6e --- /dev/null +++ b/ghc/rts/gmp/mpz/fits_ushort_p.c @@ -0,0 +1,41 @@ +/* int mpz_fits_X_p (mpz_t src) -- Return whether src fits the C type X. + +Copyright (C) 1997, 2000 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 2.1 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; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +int +#if __STDC__ +mpz_fits_ushort_p (mpz_srcptr src) +#else +mpz_fits_ushort_p (src) + mpz_srcptr src; +#endif +{ + mp_size_t size; + mp_limb_t mpl; + + mpl = PTR(src)[0]; + size = SIZ(src); + if (size < 0 || size > 1) + return 0; + return mpl <= ((unsigned short int) ~(unsigned int) 0); +} diff --git a/ghc/rts/gmp/mpz/gcd.c b/ghc/rts/gmp/mpz/gcd.c index f93030c8f4..0d950dd609 100644 --- a/ghc/rts/gmp/mpz/gcd.c +++ b/ghc/rts/gmp/mpz/gcd.c @@ -1,20 +1,20 @@ /* mpz/gcd.c: Calculate the greatest common divisor of two integers. -Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc. +Copyright (C) 1991, 1993, 1994, 1996, 2000 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -22,8 +22,10 @@ MA 02111-1307, USA. */ #include "gmp.h" #include "gmp-impl.h" #include "longlong.h" +#ifdef BERKELEY_MP +#include "mp.h" +#endif -void *_mpz_realloc (); #ifndef BERKELEY_MP void @@ -145,10 +147,10 @@ gcd (u, v, g) g_zero_bits = MIN (u_zero_bits, v_zero_bits); } - /* Call mpn_gcd. The 1st argument must not have more bits than the 2nd. */ + /* Call mpn_gcd. The 2nd argument must not have more bits than the 1st. */ vsize = (usize < vsize || (usize == vsize && up[usize-1] < vp[vsize-1])) - ? mpn_gcd (vp, up, usize, vp, vsize) - : mpn_gcd (vp, vp, vsize, up, usize); + ? mpn_gcd (vp, vp, vsize, up, usize) + : mpn_gcd (vp, up, usize, vp, vsize); /* Here G <-- V << (g_zero_limbs*BITS_PER_MP_LIMB + g_zero_bits). */ gsize = vsize + g_zero_limbs; diff --git a/ghc/rts/gmp/mpz/gcd_ui.c b/ghc/rts/gmp/mpz/gcd_ui.c index 388ab05cc9..f3bec58829 100644 --- a/ghc/rts/gmp/mpz/gcd_ui.c +++ b/ghc/rts/gmp/mpz/gcd_ui.c @@ -1,24 +1,25 @@ /* mpz_gcd_ui -- Calculate the greatest common divisior of two integers. -Copyright (C) 1994, 1996 Free Software Foundation, Inc. +Copyright (C) 1994, 1996, 1999, 2000 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ +#include <stdio.h> /* for NULL */ #include "gmp.h" #include "gmp-impl.h" @@ -53,7 +54,7 @@ mpz_gcd_ui (w, u, v) return size > 1 ? 0 : w->_mp_d[0]; } else - res = mpn_gcd_1 (u->_mp_d, size, v); + res = mpn_gcd_1 (u->_mp_d, size, (mp_limb_t) v); if (w != NULL) { diff --git a/ghc/rts/gmp/mpz/gcdext.c b/ghc/rts/gmp/mpz/gcdext.c index adf66b04af..3ba04c84ff 100644 --- a/ghc/rts/gmp/mpz/gcdext.c +++ b/ghc/rts/gmp/mpz/gcdext.c @@ -1,30 +1,30 @@ /* mpz_gcdext(g, s, t, a, b) -- Set G to gcd(a, b), and S and T such that g = as + bt. -Copyright (C) 1991, 1993, 1994, 1995 Free Software Foundation, Inc. +Copyright (C) 1991, 1993, 1994, 1995, 1996, 1997, 2000 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ +#include <stdio.h> /* for NULL */ #include "gmp.h" #include "gmp-impl.h" -/* Botch: SLOW! */ - void #if __STDC__ mpz_gcdext (mpz_ptr g, mpz_ptr s, mpz_ptr t, mpz_srcptr a, mpz_srcptr b) @@ -37,52 +37,101 @@ mpz_gcdext (g, s, t, a, b) mpz_srcptr b; #endif { - mpz_t s0, s1, q, r, x, d0, d1; + mp_size_t asize, bsize, usize, vsize; + mp_srcptr ap, bp; + mp_ptr up, vp; + mp_size_t gsize, ssize, tmp_ssize; + mp_ptr gp, sp, tmp_gp, tmp_sp; + mpz_srcptr u, v; + mpz_ptr ss, tt; + __mpz_struct stmp, gtmp; + TMP_DECL (marker); + + TMP_MARK (marker); + + /* mpn_gcdext requires that U >= V. Therefore, we often have to swap U and + V. This in turn leads to a lot of complications. The computed cofactor + will be the wrong one, so we have to fix that up at the end. */ + + asize = ABS (SIZ (a)); + bsize = ABS (SIZ (b)); + ap = PTR (a); + bp = PTR (b); + if (asize > bsize || (asize == bsize && mpn_cmp (ap, bp, asize) > 0)) + { + usize = asize; + vsize = bsize; + up = (mp_ptr) TMP_ALLOC ((usize + 1) * BYTES_PER_MP_LIMB); + vp = (mp_ptr) TMP_ALLOC ((vsize + 1) * BYTES_PER_MP_LIMB); + MPN_COPY (up, ap, usize); + MPN_COPY (vp, bp, vsize); + u = a; + v = b; + ss = s; + tt = t; + } + else + { + usize = bsize; + vsize = asize; + up = (mp_ptr) TMP_ALLOC ((usize + 1) * BYTES_PER_MP_LIMB); + vp = (mp_ptr) TMP_ALLOC ((vsize + 1) * BYTES_PER_MP_LIMB); + MPN_COPY (up, bp, usize); + MPN_COPY (vp, ap, vsize); + u = b; + v = a; + ss = t; + tt = s; + } - mpz_init_set_ui (s0, 1L); - mpz_init_set_ui (s1, 0L); - mpz_init (q); - mpz_init (r); - mpz_init (x); - mpz_init_set (d0, a); - mpz_init_set (d1, b); + tmp_gp = (mp_ptr) TMP_ALLOC ((usize + 1) * BYTES_PER_MP_LIMB); + tmp_sp = (mp_ptr) TMP_ALLOC ((usize + 1) * BYTES_PER_MP_LIMB); - while (d1->_mp_size != 0) + if (vsize == 0) { - mpz_tdiv_qr (q, r, d0, d1); - mpz_set (d0, d1); - mpz_set (d1, r); - - mpz_mul (x, s1, q); - mpz_sub (x, s0, x); - mpz_set (s0, s1); - mpz_set (s1, x); + tmp_sp[0] = 1; + tmp_ssize = 1; + MPN_COPY (tmp_gp, up, usize); + gsize = usize; } + else + gsize = mpn_gcdext (tmp_gp, tmp_sp, &tmp_ssize, up, usize, vp, vsize); + ssize = ABS (tmp_ssize); + + PTR (>mp) = tmp_gp; + SIZ (>mp) = gsize; - if (t != NULL) + PTR (&stmp) = tmp_sp; + SIZ (&stmp) = (tmp_ssize ^ SIZ (u)) >= 0 ? ssize : -ssize; + + if (tt != NULL) { - mpz_mul (x, s0, a); - mpz_sub (x, d0, x); - if (b->_mp_size == 0) - t->_mp_size = 0; + if (SIZ (v) == 0) + SIZ (tt) = 0; else - mpz_tdiv_q (t, x, b); + { + mpz_t x; + MPZ_TMP_INIT (x, ssize + usize + 1); + mpz_mul (x, &stmp, u); + mpz_sub (x, >mp, x); + mpz_tdiv_q (tt, x, v); + } } - mpz_set (s, s0); - mpz_set (g, d0); - if (g->_mp_size < 0) + + if (ss != NULL) { - g->_mp_size = -g->_mp_size; - s->_mp_size = -s->_mp_size; - if (t != NULL) - t->_mp_size = -t->_mp_size; + if (ALLOC (ss) < ssize) + _mpz_realloc (ss, ssize); + sp = PTR (ss); + MPN_COPY (sp, tmp_sp, ssize); + SIZ (ss) = SIZ (&stmp); } - mpz_clear (s0); - mpz_clear (s1); - mpz_clear (q); - mpz_clear (r); - mpz_clear (x); - mpz_clear (d0); - mpz_clear (d1); + if (ALLOC (g) < gsize) + _mpz_realloc (g, gsize); + gp = PTR (g); + MPN_COPY (gp, tmp_gp, gsize); + SIZ (g) = gsize; + + TMP_FREE (marker); } diff --git a/ghc/rts/gmp/mpz/get_d.c b/ghc/rts/gmp/mpz/get_d.c index 0fd7916eaf..6a7c5856bb 100644 --- a/ghc/rts/gmp/mpz/get_d.c +++ b/ghc/rts/gmp/mpz/get_d.c @@ -1,26 +1,49 @@ /* double mpz_get_d (mpz_t src) -- Return the double approximation to SRC. -Copyright (C) 1996 Free Software Foundation, Inc. +Copyright (C) 1996, 1997, 2000 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ #include "gmp.h" #include "gmp-impl.h" +#include "longlong.h" + + +static int +#if __STDC__ +mpn_zero_p (mp_ptr p, mp_size_t n) +#else +mpn_zero_p (p, n) + mp_ptr p; + mp_size_t n; +#endif +{ + mp_size_t i; + + for (i = 0; i < n; i++) + { + if (p[i] != 0) + return 0; + } + + return 1; +} + double #if __STDC__ @@ -31,9 +54,11 @@ mpz_get_d (src) #endif { double res; - mp_size_t size, i, n_limbs_to_use; + mp_size_t size; int negative; mp_ptr qp; + mp_limb_t hz, lz; + int cnt; size = SIZ(src); if (size == 0) @@ -43,12 +68,61 @@ mpz_get_d (src) size = ABS (size); qp = PTR(src); - res = qp[size - 1]; - n_limbs_to_use = MIN (LIMBS_PER_DOUBLE, size); - for (i = 2; i <= n_limbs_to_use; i++) - res = res * MP_BASE_AS_DOUBLE + qp[size - i]; + if (size == 1) + { + res = qp[size - 1]; + } + else if (size == 2) + { + res = MP_BASE_AS_DOUBLE * qp[size - 1] + qp[size - 2]; + } + else + { + count_leading_zeros (cnt, qp[size - 1]); - res = __gmp_scale2 (res, (size - n_limbs_to_use) * BITS_PER_MP_LIMB); +#if BITS_PER_MP_LIMB == 32 + if (cnt == 0) + { + hz = qp[size - 1]; + lz = qp[size - 2]; + } + else + { + hz = (qp[size - 1] << cnt) | (qp[size - 2] >> BITS_PER_MP_LIMB - cnt); + lz = (qp[size - 2] << cnt) | (qp[size - 3] >> BITS_PER_MP_LIMB - cnt); + } +#if _GMP_IEEE_FLOATS + /* Take bits from less significant limbs, but only if they may affect + the result. */ + if ((lz & 0x7ff) == 0x400) + { + if (cnt != 0) + lz += ((qp[size - 3] << cnt) != 0 || ! mpn_zero_p (qp, size - 3)); + else + lz += (! mpn_zero_p (qp, size - 2)); + } +#endif + res = MP_BASE_AS_DOUBLE * hz + lz; + res = __gmp_scale2 (res, (size - 2) * BITS_PER_MP_LIMB - cnt); +#endif +#if BITS_PER_MP_LIMB == 64 + if (cnt == 0) + hz = qp[size - 1]; + else + hz = (qp[size - 1] << cnt) | (qp[size - 2] >> BITS_PER_MP_LIMB - cnt); +#if _GMP_IEEE_FLOATS + if ((hz & 0x7ff) == 0x400) + { + if (cnt != 0) + hz += ((qp[size - 2] << cnt) != 0 || ! mpn_zero_p (qp, size - 2)); + else + hz += (! mpn_zero_p (qp, size - 1)); + } +#endif + res = hz; + res = __gmp_scale2 (res, (size - 1) * BITS_PER_MP_LIMB - cnt); +#endif + } return negative ? -res : res; } diff --git a/ghc/rts/gmp/mpz/get_si.c b/ghc/rts/gmp/mpz/get_si.c index 45e0e5aecf..8a5d0e4803 100644 --- a/ghc/rts/gmp/mpz/get_si.c +++ b/ghc/rts/gmp/mpz/get_si.c @@ -5,16 +5,16 @@ Copyright (C) 1991, 1993, 1994, 1995 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/get_str.c b/ghc/rts/gmp/mpz/get_str.c index 8ccf3efcd0..c7278afb52 100644 --- a/ghc/rts/gmp/mpz/get_str.c +++ b/ghc/rts/gmp/mpz/get_str.c @@ -9,16 +9,16 @@ Copyright (C) 1991, 1993, 1994, 1996 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/get_ui.c b/ghc/rts/gmp/mpz/get_ui.c index 4bfb5e14db..a8ec9e01a4 100644 --- a/ghc/rts/gmp/mpz/get_ui.c +++ b/ghc/rts/gmp/mpz/get_ui.c @@ -5,16 +5,16 @@ Copyright (C) 1991, 1993, 1994, 1995 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/getlimbn.c b/ghc/rts/gmp/mpz/getlimbn.c index c7a234b5ab..b772ed05c4 100644 --- a/ghc/rts/gmp/mpz/getlimbn.c +++ b/ghc/rts/gmp/mpz/getlimbn.c @@ -1,20 +1,20 @@ /* mpz_getlimbn(integer,n) -- Return the N:th limb from INTEGER. -Copyright (C) 1993, 1994, 1995 Free Software Foundation, Inc. +Copyright (C) 1993, 1994, 1995, 1996 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -31,7 +31,7 @@ mpz_getlimbn (integer, n) mp_size_t n; #endif { - if (integer->_mp_size <= n || n < 0) + if (ABS (integer->_mp_size) <= n || n < 0) return 0; else return integer->_mp_d[n]; diff --git a/ghc/rts/gmp/mpz/hamdist.c b/ghc/rts/gmp/mpz/hamdist.c index 58c927316a..b039a653d2 100644 --- a/ghc/rts/gmp/mpz/hamdist.c +++ b/ghc/rts/gmp/mpz/hamdist.c @@ -8,16 +8,16 @@ Copyright (C) 1994, 1996 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/init.c b/ghc/rts/gmp/mpz/init.c index f8d8e207e4..2e8e4d2cbd 100644 --- a/ghc/rts/gmp/mpz/init.c +++ b/ghc/rts/gmp/mpz/init.c @@ -5,16 +5,16 @@ Copyright (C) 1991, 1993, 1994, 1995 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/inp_raw.c b/ghc/rts/gmp/mpz/inp_raw.c index e1cec1d743..15e601229d 100644 --- a/ghc/rts/gmp/mpz/inp_raw.c +++ b/ghc/rts/gmp/mpz/inp_raw.c @@ -6,16 +6,16 @@ Copyright (C) 1991, 1993, 1994, 1995 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/inp_str.c b/ghc/rts/gmp/mpz/inp_str.c index 715906237e..7aa5e1fc30 100644 --- a/ghc/rts/gmp/mpz/inp_str.c +++ b/ghc/rts/gmp/mpz/inp_str.c @@ -1,21 +1,21 @@ /* mpz_inp_str(dest_integer, stream, base) -- Input a number in base BASE from stdio stream STREAM and store the result in DEST_INTEGER. -Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc. +Copyright (C) 1991, 1993, 1994, 1996, 1998, 2000 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -26,9 +26,13 @@ MA 02111-1307, USA. */ #include "gmp-impl.h" static int +#if __STDC__ +digit_value_in_base (int c, int base) +#else digit_value_in_base (c, base) int c; int base; +#endif { int digit; @@ -48,10 +52,10 @@ digit_value_in_base (c, base) size_t #if __STDC__ -mpz_inp_str (mpz_ptr dest, FILE *stream, int base) +mpz_inp_str (mpz_ptr x, FILE *stream, int base) #else -mpz_inp_str (dest, stream, base) - mpz_ptr dest; +mpz_inp_str (x, stream, base) + mpz_ptr x; FILE *stream; int base; #endif @@ -60,15 +64,12 @@ mpz_inp_str (dest, stream, base) size_t alloc_size, str_size; int c; int negative; - mp_size_t dest_size; + mp_size_t xsize; size_t nread; if (stream == 0) stream = stdin; - alloc_size = 100; - str = (char *) (*_mp_allocate_func) (alloc_size); - str_size = 0; nread = 0; /* Skip whitespace. */ @@ -84,6 +85,7 @@ mpz_inp_str (dest, stream, base) { negative = 1; c = getc (stream); + nread++; } if (digit_value_in_base (c, base == 0 ? 10 : base) < 0) @@ -105,9 +107,26 @@ mpz_inp_str (dest, stream, base) c = getc (stream); nread++; } + else if (c == 'b' || c == 'B') + { + base = 2; + c = getc (stream); + nread++; + } } } + /* Skip leading zeros. */ + while (c == '0') + { + c = getc (stream); + nread++; + } + + alloc_size = 100; + str = (char *) (*_mp_allocate_func) (alloc_size); + str_size = 0; + for (;;) { int dig; @@ -126,12 +145,22 @@ mpz_inp_str (dest, stream, base) ungetc (c, stream); - dest_size = str_size / __mp_bases[base].chars_per_limb + 1; - if (dest->_mp_alloc < dest_size) - _mpz_realloc (dest, dest_size); + /* Make sure the string is not empty, mpn_set_str would fail. */ + if (str_size == 0) + { + x->_mp_size = 0; + (*_mp_free_func) (str, alloc_size); + return nread; + } + + xsize = (((mp_size_t) (str_size / __mp_bases[base].chars_per_bit_exactly)) + / BITS_PER_MP_LIMB + 2); + if (x->_mp_alloc < xsize) + _mpz_realloc (x, xsize); - dest_size = mpn_set_str (dest->_mp_d, (unsigned char *) str, str_size, base); - dest->_mp_size = negative ? -dest_size : dest_size; + /* Convert the byte array in base BASE to our bignum format. */ + xsize = mpn_set_str (x->_mp_d, (unsigned char *) str, str_size, base); + x->_mp_size = negative ? -xsize : xsize; (*_mp_free_func) (str, alloc_size); return str_size + nread; diff --git a/ghc/rts/gmp/mpz/invert.c b/ghc/rts/gmp/mpz/invert.c index ff1d6d937d..749a0969fc 100644 --- a/ghc/rts/gmp/mpz/invert.c +++ b/ghc/rts/gmp/mpz/invert.c @@ -1,27 +1,28 @@ /* mpz_invert (inv, x, n). Find multiplicative inverse of X in Z(N). If X has an inverse, return non-zero and store inverse in INVERSE, - otherwise, return 0 and put garbage in X. + otherwise, return 0 and put garbage in INVERSE. -Copyright (C) 1996 Free Software Foundation, Inc. +Copyright (C) 1996, 1997, 1998, 1999, 2000 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ #include "gmp.h" +#include "gmp-impl.h" int #if __STDC__ @@ -32,12 +33,45 @@ mpz_invert (inverse, x, n) mpz_srcptr x, n; #endif { - mpz_t gcd; - int rv; - - mpz_init (gcd); - mpz_gcdext (gcd, inverse, (mpz_ptr) 0, x, n); - rv = gcd->_mp_size == 1 && (gcd->_mp_d)[0] == 1; - mpz_clear (gcd); - return rv; + mpz_t gcd, tmp; + mp_size_t xsize, nsize, size; + TMP_DECL (marker); + + xsize = SIZ (x); + nsize = SIZ (n); + xsize = ABS (xsize); + nsize = ABS (nsize); + size = MAX (xsize, nsize) + 1; + + /* No inverse exists if the leftside operand is 0. Likewise, no + inverse exists if the mod operand is 1. */ + if (xsize == 0 || (nsize == 1 && (PTR (n))[0] == 1)) + return 0; + + TMP_MARK (marker); + + MPZ_TMP_INIT (gcd, size); + MPZ_TMP_INIT (tmp, size); + mpz_gcdext (gcd, tmp, (mpz_ptr) 0, x, n); + + /* If no inverse existed, return with an indication of that. */ + if (gcd->_mp_size != 1 || (gcd->_mp_d)[0] != 1) + { + TMP_FREE (marker); + return 0; + } + + /* Make sure we return a positive inverse. */ + if (SIZ (tmp) < 0) + { + if (SIZ (n) < 0) + mpz_sub (inverse, tmp, n); + else + mpz_add (inverse, tmp, n); + } + else + mpz_set (inverse, tmp); + + TMP_FREE (marker); + return 1; } diff --git a/ghc/rts/gmp/mpz/ior.c b/ghc/rts/gmp/mpz/ior.c index 77facfd3d4..0bb5a806dc 100644 --- a/ghc/rts/gmp/mpz/ior.c +++ b/ghc/rts/gmp/mpz/ior.c @@ -1,20 +1,21 @@ /* mpz_ior -- Logical inclusive or. -Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc. +Copyright (C) 1991, 1993, 1994, 1996, 1997, 2000 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -126,7 +127,7 @@ mpz_ior (res, op1, op2) _mpz_realloc (res, res_size); res_ptr = res->_mp_d; /* Don't re-read OP1_PTR and OP2_PTR. They point to - temporary space--never to the space RES->_mp_D used + temporary space--never to the space RES->_mp_d used to point to before reallocation. */ } @@ -163,9 +164,8 @@ mpz_ior (res, op1, op2) { /* We should compute -OP1 | OP2. Swap OP1 and OP2 and fall through to the code that handles OP1 | -OP2. */ - {mpz_srcptr t = op1; op1 = op2; op2 = t;} - {mp_srcptr t = op1_ptr; op1_ptr = op2_ptr; op2_ptr = t;} - {mp_size_t t = op1_size; op1_size = op2_size; op2_size = t;} + MPZ_SRCPTR_SWAP (op1, op2); + MPN_SRCPTR_SWAP (op1_ptr,op1_size, op2_ptr,op2_size); } } @@ -187,6 +187,7 @@ mpz_ior (res, op1, op2) opx = (mp_ptr) TMP_ALLOC (op2_size * BYTES_PER_MP_LIMB); mpn_sub_1 (opx, op2_ptr, op2_size, (mp_limb_t) 1); op2_ptr = opx; + op2_size -= op2_ptr[op2_size - 1] == 0; if (res->_mp_alloc < res_alloc) { @@ -194,7 +195,7 @@ mpz_ior (res, op1, op2) op1_ptr = op1->_mp_d; res_ptr = res->_mp_d; /* Don't re-read OP2_PTR. It points to temporary space--never - to the space RES->_mp_D used to point to before reallocation. */ + to the space RES->_mp_d used to point to before reallocation. */ } if (op1_size >= op2_size) diff --git a/ghc/rts/gmp/mpz/iset.c b/ghc/rts/gmp/mpz/iset.c index c8a17dc67a..114bc2d542 100644 --- a/ghc/rts/gmp/mpz/iset.c +++ b/ghc/rts/gmp/mpz/iset.c @@ -6,16 +6,16 @@ Copyright (C) 1991, 1993, 1994, 1996 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/iset_d.c b/ghc/rts/gmp/mpz/iset_d.c index 41e5c4f963..502a8933e2 100644 --- a/ghc/rts/gmp/mpz/iset_d.c +++ b/ghc/rts/gmp/mpz/iset_d.c @@ -6,16 +6,16 @@ Copyright (C) 1996 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/iset_si.c b/ghc/rts/gmp/mpz/iset_si.c index af51f0578a..842db140ef 100644 --- a/ghc/rts/gmp/mpz/iset_si.c +++ b/ghc/rts/gmp/mpz/iset_si.c @@ -1,21 +1,21 @@ /* mpz_init_set_si(val) -- Make a new multiple precision number with value val. -Copyright (C) 1991, 1993, 1994, 1995 Free Software Foundation, Inc. +Copyright (C) 1991, 1993, 1994, 1995, 2000 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -41,7 +41,7 @@ mpz_init_set_si (x, val) } else if (val < 0) { - x->_mp_d[0] = -val; + x->_mp_d[0] = (unsigned long) -val; x->_mp_size = -1; } else diff --git a/ghc/rts/gmp/mpz/iset_str.c b/ghc/rts/gmp/mpz/iset_str.c index e04ad5da22..dfb8c6b230 100644 --- a/ghc/rts/gmp/mpz/iset_str.c +++ b/ghc/rts/gmp/mpz/iset_str.c @@ -5,21 +5,21 @@ i.e. 0xhh...h means base 16, 0oo...o means base 8, otherwise assume base 10. -Copyright (C) 1991, 1993, 1994, 1995 Free Software Foundation, Inc. +Copyright (C) 1991, 1993, 1994, 1995, 2000 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -40,5 +40,8 @@ mpz_init_set_str (x, str, base) x->_mp_alloc = 1; x->_mp_d = (mp_ptr) (*_mp_allocate_func) (BYTES_PER_MP_LIMB); + /* if str has no digits mpz_set_str leaves x->_mp_size unset */ + x->_mp_size = 0; + return mpz_set_str (x, str, base); } diff --git a/ghc/rts/gmp/mpz/iset_ui.c b/ghc/rts/gmp/mpz/iset_ui.c index dc39f59ff4..759182c556 100644 --- a/ghc/rts/gmp/mpz/iset_ui.c +++ b/ghc/rts/gmp/mpz/iset_ui.c @@ -6,16 +6,16 @@ Copyright (C) 1991, 1993, 1994, 1995 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/jacobi.c b/ghc/rts/gmp/mpz/jacobi.c index 409f622fce..9d49e1d0c6 100644 --- a/ghc/rts/gmp/mpz/jacobi.c +++ b/ghc/rts/gmp/mpz/jacobi.c @@ -6,16 +6,16 @@ Copyright (C) 1991, 1996 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/kronsz.c b/ghc/rts/gmp/mpz/kronsz.c new file mode 100644 index 0000000000..c8c6752224 --- /dev/null +++ b/ghc/rts/gmp/mpz/kronsz.c @@ -0,0 +1,126 @@ +/* mpz_si_kronecker -- Kronecker/Jacobi symbol. */ + +/* +Copyright (C) 1999, 2000 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 2.1 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; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. +*/ + +#include "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" + + +int +#if __STDC__ +mpz_si_kronecker (long a, mpz_srcptr b) +#else +mpz_si_kronecker (a, b) + long a; + mpz_srcptr b; +#endif +{ + int b_abs_size; + mp_srcptr b_ptr; + mp_limb_t b_low; + int twos; + int result_bit1; + + b_abs_size = ABSIZ (b); + if (b_abs_size == 0) + return JACOBI_S0 (a); /* (a/0) */ + + b_ptr = PTR(b); + b_low = b_ptr[0]; + + /* (0/b) = 1 if b=+/-1, 0 otherwise */ + if (a == 0) + return (b_abs_size == 1) & (b_low == 1); + + /* account for the effect of the sign of b, so can then ignore it */ + result_bit1 = JACOBI_BSGN_SZ_BIT1 (a, b); + + if ((b_low & 1) == 0) + { + /* b even */ + + if ((a & 1) == 0) + return 0; /* (a/b)=0 if both a,b even */ + + /* Require MP_BITS_PER_LIMB even, so that (a/2)^MP_BITS_PER_LIMB = 1, + and so that therefore there's no need to account for how many zero + limbs are stripped. */ + ASSERT ((BITS_PER_MP_LIMB & 1) == 0); + + MPN_STRIP_LOW_ZEROS_NOT_ZERO (b_ptr, b_abs_size); + b_low = b_ptr[0]; + + if ((b_low & 1) == 0) + { + /* odd a, even b */ + + mp_limb_t b_shl_bit1; + + count_trailing_zeros (twos, b_low); + + /* b_shl_bit1 is b>>twos, but with only bit 1 guaranteed */ + if (twos == BITS_PER_MP_LIMB-1) + b_shl_bit1 = (b_abs_size == 1) ? 0 : (b_ptr[1] << 1); + else + b_shl_bit1 = (b_low >> twos); + + result_bit1 ^= JACOBI_ASGN_SU_BIT1 (a, b_shl_bit1); + a = ABS(a); + + if (a == 1) + return JACOBI_BIT1_TO_PN (result_bit1); /* (1/b)=1 */ + + /* twos (a/2), reciprocity to (b/a), and (b/a) = (b mod a / b) */ + return mpn_jacobi_base (mpn_mod_1_rshift (b_ptr, b_abs_size, + twos, a), + a, + result_bit1 + ^ JACOBI_TWOS_U_BIT1 (twos, a) + ^ JACOBI_RECIP_UU_BIT1 (a, b_shl_bit1)); + } + } + + /* b odd */ + + result_bit1 ^= JACOBI_ASGN_SU_BIT1 (a, b_low); + a = ABS(a); + + /* (a/1) = 1 for any a */ + if (b_abs_size == 1 && b_low == 1) + return JACOBI_BIT1_TO_PN (result_bit1); + + /* Note a is cast to unsigned because 0x80..00 doesn't fit in a signed. */ + if ((a & 1) == 0) + { + count_trailing_zeros (twos, a); + a = ((unsigned long) a) >> twos; + result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, b_low); + } + + if (a == 1) + return JACOBI_BIT1_TO_PN (result_bit1); /* (1/b)=1 */ + + /* reciprocity to (b/a), and (b/a) == (b mod a / a) */ + return mpn_jacobi_base (mpn_mod_1 (b_ptr, b_abs_size, a), a, + result_bit1 ^ JACOBI_RECIP_UU_BIT1 (a, b_low)); +} diff --git a/ghc/rts/gmp/mpz/kronuz.c b/ghc/rts/gmp/mpz/kronuz.c new file mode 100644 index 0000000000..b877e6f64c --- /dev/null +++ b/ghc/rts/gmp/mpz/kronuz.c @@ -0,0 +1,115 @@ +/* mpz_ui_kronecker -- Kronecker/Jacobi symbol. */ + +/* +Copyright (C) 1999, 2000 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 2.1 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; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. +*/ + +#include "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" + + +int +#if __STDC__ +mpz_ui_kronecker (unsigned long a, mpz_srcptr b) +#else +mpz_ui_kronecker (a, b) + unsigned long a; + mpz_srcptr b; +#endif +{ + int b_abs_size; + mp_srcptr b_ptr; + mp_limb_t b_low; + int twos; + int result_bit1; + + /* (a/0) */ + b_abs_size = ABSIZ (b); + if (b_abs_size == 0) + return JACOBI_U0 (a); + + /* (a/-1)=1 when a>=0, so the sign of b is ignored */ + b_ptr = PTR(b); + b_low = b_ptr[0]; + + /* (0/1)=1; (0/-1)=1; (0/b)=0 for b!=+/-1 + (1/b)=1, for any b */ + if (a <= 1) + return (a == 1) | ((b_abs_size == 1) & (b_low == 1)); + + if (b_low & 1) + { + /* (a/1) = 1 for any a */ + if (b_abs_size == 1 && b_low == 1) + return 1; + + count_trailing_zeros (twos, a); + a >>= twos; + if (a == 1) + return JACOBI_TWOS_U (twos, b_low); /* powers of (2/b) only */ + + /* powers of (2/b); reciprocity to (b/a); (b/a) == (b mod a / a) */ + return mpn_jacobi_base (mpn_mod_1 (b_ptr, b_abs_size, a), + a, + JACOBI_TWOS_U_BIT1 (twos, b_low) + ^ JACOBI_RECIP_UU_BIT1 (b_low, a)); + } + + /* b is even; (a/2)=0 if a is even */ + if ((a & 1) == 0) + return 0; + + /* Require MP_BITS_PER_LIMB even, so (a/2)^MP_BITS_PER_LIMB = 1, and so we + don't have to pay attention to how many trailing zero limbs are + stripped. */ + ASSERT ((BITS_PER_MP_LIMB & 1) == 0); + + MPN_STRIP_LOW_ZEROS_NOT_ZERO (b_ptr, b_abs_size); + b_low = b_ptr[0]; + + if (b_low & 1) + /* reciprocity to (b/a); (b/a) == (b mod a / a) */ + return mpn_jacobi_base (mpn_mod_1 (b_ptr, b_abs_size, a), + a, + JACOBI_RECIP_UU_BIT1 (b_low, a)); + + count_trailing_zeros (twos, b_low); + + /* reciprocity to get (b/a) */ + if (twos == BITS_PER_MP_LIMB-1) + { + if (b_abs_size == 1) + { + /* b==0x800...00, one limb high bit only, so (a/2)^(BPML-1) */ + return JACOBI_TWOS_U (BITS_PER_MP_LIMB-1, a); + } + + /* b_abs_size > 1 */ + result_bit1 = JACOBI_RECIP_UU_BIT1 (a, b_ptr[1] << 1); + } + else + result_bit1 = JACOBI_RECIP_UU_BIT1 (a, b_low >> twos); + + /* powers of (a/2); reciprocity to (b/a); (b/a) == (b mod a / a) */ + return mpn_jacobi_base (mpn_mod_1_rshift (b_ptr, b_abs_size, twos, a), + a, + JACOBI_TWOS_U_BIT1 (twos, a) ^ result_bit1); +} diff --git a/ghc/rts/gmp/mpz/kronzs.c b/ghc/rts/gmp/mpz/kronzs.c new file mode 100644 index 0000000000..edfb465976 --- /dev/null +++ b/ghc/rts/gmp/mpz/kronzs.c @@ -0,0 +1,74 @@ +/* mpz_kronecker_si -- Kronecker/Jacobi symbol. */ + +/* +Copyright (C) 1999, 2000 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 2.1 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; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. +*/ + +#include "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" + + +/* This function is expected to be often used with b odd, so there's a test + for this before invoking count_trailing_zeros(). + + After the absolute value of b is established it's treated as an unsigned + long, because 0x80..00 doesn't fit in a signed long. */ + +int +#if __STDC__ +mpz_kronecker_si (mpz_srcptr a, long b) +#else +mpz_kronecker_si (a, b) + mpz_srcptr a; + long b; +#endif +{ + int result_bit1; + int twos; + + if (b == 0) + return JACOBI_Z0 (a); + + result_bit1 = JACOBI_BSGN_ZS_BIT1(a, b); + b = ABS (b); + + if (b == 1) + return JACOBI_BIT1_TO_PN (result_bit1); /* (a/1) = 1 for any a */ + + if (b & 1) + return mpn_jacobi_base (mpz_fdiv_ui (a, b), b, result_bit1); + + /* result 0 if both a,b even */ + if (mpz_even_p (a)) + return 0; + + /* (a/2)=(2/a) when a odd */ + count_trailing_zeros (twos, b); + result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, PTR(a)[0]); + + b = ((unsigned long) b) >> twos; + if (b == 1) + return JACOBI_BIT1_TO_PN (result_bit1); + else + return mpn_jacobi_base (mpz_fdiv_ui (a, b), b, result_bit1); +} + + diff --git a/ghc/rts/gmp/mpz/kronzu.c b/ghc/rts/gmp/mpz/kronzu.c new file mode 100644 index 0000000000..749be5df07 --- /dev/null +++ b/ghc/rts/gmp/mpz/kronzu.c @@ -0,0 +1,66 @@ +/* mpz_kronecker_ui -- Kronecker/Jacobi symbol. */ + +/* +Copyright (C) 1999, 2000 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 2.1 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; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. +*/ + +#include "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" + + +/* This function is expected to be often used with b an odd prime, so the + code for odd b is nice and short. */ + +int +#if __STDC__ +mpz_kronecker_ui (mpz_srcptr a, unsigned long b) +#else +mpz_kronecker_ui (a, b) + mpz_srcptr a; + unsigned long b; +#endif +{ + int twos; + + if (b & 1) + { + if (b != 1) + return mpn_jacobi_base (mpz_fdiv_ui (a, b), b, 0); + else + return 1; /* (a/1)=1 for any a */ + } + + if (b == 0) + return JACOBI_Z0 (a); + + /* (a/2)=0 if a even */ + if (mpz_even_p (a)) + return 0; + + /* (a/2)=(2/a) when a odd */ + count_trailing_zeros (twos, b); + b >>= twos; + if (b == 1) + return JACOBI_TWOS_U (twos, PTR(a)[0]); + + return mpn_jacobi_base (mpz_fdiv_ui (a, b), b, + JACOBI_TWOS_U_BIT1(twos, PTR(a)[0])); +} diff --git a/ghc/rts/gmp/mpz/lcm.c b/ghc/rts/gmp/mpz/lcm.c new file mode 100644 index 0000000000..ca8c15cac9 --- /dev/null +++ b/ghc/rts/gmp/mpz/lcm.c @@ -0,0 +1,56 @@ +/* mpz/lcm.c: Calculate the least common multiple of two integers. + +Copyright (C) 1996 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 2.1 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; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" + +void *_mpz_realloc (); + +void +#if __STDC__ +mpz_lcm (mpz_ptr r, mpz_srcptr u, mpz_srcptr v) +#else +mpz_lcm (r, u, v) + mpz_ptr r; + mpz_srcptr u; + mpz_srcptr v; +#endif +{ + mpz_t g; + mp_size_t usize, vsize, size; + + usize = ABS (SIZ (u)); + vsize = ABS (SIZ (v)); + + if (usize == 0 || vsize == 0) + { + SIZ (r) = 0; + return; + } + + size = MAX (usize, vsize); + MPZ_TMP_INIT (g, size); + + mpz_gcd (g, u, v); + mpz_divexact (g, u, g); + mpz_mul (r, g, v); +} diff --git a/ghc/rts/gmp/mpz/legendre.c b/ghc/rts/gmp/mpz/legendre.c index 4de16a6314..ab665f70d0 100644 --- a/ghc/rts/gmp/mpz/legendre.c +++ b/ghc/rts/gmp/mpz/legendre.c @@ -6,16 +6,16 @@ Copyright (C) 1992, 1996 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/mod.c b/ghc/rts/gmp/mpz/mod.c index b2b8b39bbe..87033b333b 100644 --- a/ghc/rts/gmp/mpz/mod.c +++ b/ghc/rts/gmp/mpz/mod.c @@ -5,16 +5,16 @@ Copyright (C) 1991, 1993, 1994, 1995, 1996 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/mul.c b/ghc/rts/gmp/mpz/mul.c index 47ce8e3548..7854788e50 100644 --- a/ghc/rts/gmp/mpz/mul.c +++ b/ghc/rts/gmp/mpz/mul.c @@ -1,26 +1,30 @@ /* mpz_mul -- Multiply two integers. -Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc. +Copyright (C) 1991, 1993, 1994, 1996, 2000 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ +#include <stdio.h> /* for NULL */ #include "gmp.h" #include "gmp-impl.h" +#ifdef BERKELEY_MP +#include "mp.h" +#endif #ifndef BERKELEY_MP void diff --git a/ghc/rts/gmp/mpz/mul_2exp.c b/ghc/rts/gmp/mpz/mul_2exp.c index 4d66a9877e..abea5fed2c 100644 --- a/ghc/rts/gmp/mpz/mul_2exp.c +++ b/ghc/rts/gmp/mpz/mul_2exp.c @@ -5,16 +5,16 @@ Copyright (C) 1991, 1993, 1994, 1996 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/mul_siui.c b/ghc/rts/gmp/mpz/mul_siui.c new file mode 100644 index 0000000000..9849cd41b0 --- /dev/null +++ b/ghc/rts/gmp/mpz/mul_siui.c @@ -0,0 +1,81 @@ +/* mpz_mul_ui/si (product, multiplier, small_multiplicand) -- Set PRODUCT to + MULTIPLICATOR times SMALL_MULTIPLICAND. + +Copyright (C) 1991, 1993, 1994, 1996, 2000 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 2.1 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; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + + +#ifdef OPERATION_mul_ui +#define FUNCTION mpz_mul_ui +#define MULTIPLICAND_UNSIGNED unsigned +#define MULTIPLICAND_ABS(x) x +#else +#ifdef OPERATION_mul_si +#define FUNCTION mpz_mul_si +#define MULTIPLICAND_UNSIGNED +#define MULTIPLICAND_ABS(x) ABS(x) +#else +Error, error, unrecognised OPERATION +#endif +#endif + + +void +#if __STDC__ +FUNCTION (mpz_ptr prod, mpz_srcptr mult, + MULTIPLICAND_UNSIGNED long int small_mult) +#else +FUNCTION (prod, mult, small_mult) + mpz_ptr prod; + mpz_srcptr mult; + MULTIPLICAND_UNSIGNED long int small_mult; +#endif +{ + mp_size_t size = mult->_mp_size; + mp_size_t sign_product = size; + mp_limb_t cy; + mp_size_t prod_size; + mp_ptr prod_ptr; + + if (size == 0 || small_mult == 0) + { + prod->_mp_size = 0; + return; + } + size = ABS (size); + + prod_size = size + 1; + if (prod->_mp_alloc < prod_size) + _mpz_realloc (prod, prod_size); + + prod_ptr = prod->_mp_d; + + cy = mpn_mul_1 (prod_ptr, mult->_mp_d, size, + (mp_limb_t) MULTIPLICAND_ABS (small_mult)); + if (cy != 0) + { + prod_ptr[size] = cy; + size++; + } + + prod->_mp_size = ((sign_product < 0) ^ (small_mult < 0)) ? -size : size; +} diff --git a/ghc/rts/gmp/mpz/neg.c b/ghc/rts/gmp/mpz/neg.c index 0b48e5c7b4..566c3a95aa 100644 --- a/ghc/rts/gmp/mpz/neg.c +++ b/ghc/rts/gmp/mpz/neg.c @@ -5,16 +5,16 @@ Copyright (C) 1991, 1993, 1994, 1995 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/nextprime.c b/ghc/rts/gmp/mpz/nextprime.c new file mode 100644 index 0000000000..f024dd1206 --- /dev/null +++ b/ghc/rts/gmp/mpz/nextprime.c @@ -0,0 +1,120 @@ +/* mpz_nextprime(p,t) - compute the next prime > t and store that in p. + +Copyright (C) 1999, 2000 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 2.1 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; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +#if __STDC__ +mpz_nextprime (mpz_ptr p, mpz_srcptr t) +#else +mpz_nextprime (p, t) + mpz_ptr p; + mpz_srcptr t; +#endif +{ + mpz_add_ui (p, t, 1L); + while (! mpz_probab_prime_p (p, 5)) + mpz_add_ui (p, p, 1L); +} + +#if 0 +/* This code is not yet tested. Will be enabled in 3.1. */ + +status unsigned short primes[] = +{ +3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97, +101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181, +191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277, +281,283,293,307,311,313,317,331,337,347,349,353,359,367,373,379,383, +389,397,401,409,419,421,431,433,439,443,449,457,461,463,467,479,487, +491,499,503,509,521,523,541,547,557,563,569,571,577,587,593,599,601, +607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701,709, +719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,827, +829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941,947, +953,967,971,977,983,991,997 +}; + +#define NUMBER_OF_PRIMES 167 + +void +#if __STDC__ +mpz_nextprime (mpz_ptr p, mpz_srcptr n) +#else +mpz_nextprime (p, n) + mpz_ptr p; + mpz_srcptr n; +#endif +{ + mpz_t tmp; + unsigned short *moduli; + unsigned long difference; + int i; + int composite; + + /* First handle tiny numbers */ + if (mpz_cmp_ui (n, 2) < 0) + { + mpz_set_ui (p, 2); + return; + } + mpz_add_ui (p, n, 1); + mpz_setbit (p, 0); + + if (mpz_cmp_ui (p, 7) <= 0) + return; + + prime_limit = NUMBER_OF_PRIMES - 1; + if (mpz_cmp_ui (p, primes[prime_limit]) <= 0) + /* Just use first three entries (3,5,7) of table for small numbers */ + prime_limit = 3; + if (prime_limit) + { + /* Compute residues modulo small odd primes */ + moduli = (unsigned short *) TMP_ALLOC (prime_limit * sizeof moduli[0]); + for (i = 0; i < prime_limit; i++) + moduli[i] = mpz_fdiv_ui (p, primes[i]); + } + for (difference = 0; ; difference += 2) + { + composite = 0; + + /* First check residues */ + for (i = 0; i < prime_limit; i++) + { + int acc, pr; + composite |= (moduli[i] == 0); + acc = moduli[i] + 2; + pr = primes[i]; + moduli[i] = acc >= pr ? acc - pr : acc; + } + if (composite) + continue; + + mpz_add_ui (p, p, difference); + difference = 0; + + /* Miller-Rabin test */ + if (mpz_millerrabin (p, 2)) + break; + } +} +#endif diff --git a/ghc/rts/gmp/mpz/out_raw.c b/ghc/rts/gmp/mpz/out_raw.c index 35d311b438..62709479c5 100644 --- a/ghc/rts/gmp/mpz/out_raw.c +++ b/ghc/rts/gmp/mpz/out_raw.c @@ -6,16 +6,16 @@ Copyright (C) 1995 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/out_str.c b/ghc/rts/gmp/mpz/out_str.c index 909f53346d..bf971b0057 100644 --- a/ghc/rts/gmp/mpz/out_str.c +++ b/ghc/rts/gmp/mpz/out_str.c @@ -6,16 +6,16 @@ Copyright (C) 1991, 1993, 1994, 1996 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/perfpow.c b/ghc/rts/gmp/mpz/perfpow.c new file mode 100644 index 0000000000..e71670a0be --- /dev/null +++ b/ghc/rts/gmp/mpz/perfpow.c @@ -0,0 +1,272 @@ +/* mpz_perfect_power_p(arg) -- Return non-zero if ARG is a perfect power, + zero otherwise. + +Copyright (C) 1998, 1999, 2000 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 2.1 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; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +/* + We are to determine if c is a perfect power, c = a ^ b. + Assume c is divisible by 2^n and that codd = c/2^n is odd. + Assume a is divisible by 2^m and that aodd = a/2^m is odd. + It is always true that m divides n. + + * If n is prime, either 1) a is 2*aodd and b = n + or 2) a = c and b = 1. + So for n prime, we readily have a solution. + * If n is factorable into the non-trivial factors p1,p2,... + Since m divides n, m has a subset of n's factors and b = n / m. + + BUG: Should handle negative numbers, since they can be odd perfect powers. +*/ + +/* This is a naive approach to recognizing perfect powers. + Many things can be improved. In particular, we should use p-adic + arithmetic for computing possible roots. */ + +#include <stdio.h> /* for NULL */ +#include "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" + +static unsigned long int gcd _PROTO ((unsigned long int a, unsigned long int b)); +static int isprime _PROTO ((unsigned long int t)); + +static const unsigned short primes[] = +{ 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, + 59, 61, 67, 71, 73, 79, 83, 89, 97,101,103,107,109,113,127,131, + 137,139,149,151,157,163,167,173,179,181,191,193,197,199,211,223, + 227,229,233,239,241,251,257,263,269,271,277,281,283,293,307,311, + 313,317,331,337,347,349,353,359,367,373,379,383,389,397,401,409, + 419,421,431,433,439,443,449,457,461,463,467,479,487,491,499,503, + 509,521,523,541,547,557,563,569,571,577,587,593,599,601,607,613, + 617,619,631,641,643,647,653,659,661,673,677,683,691,701,709,719, + 727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,827, + 829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941, + 947,953,967,971,977,983,991,997,0 +}; +#define SMALLEST_OMITTED_PRIME 1009 + + +int +#if __STDC__ +mpz_perfect_power_p (mpz_srcptr u) +#else +mpz_perfect_power_p (u) + mpz_srcptr u; +#endif +{ + unsigned long int prime; + unsigned long int n, n2; + int i; + unsigned long int rem; + mpz_t u2, q; + int exact; + mp_size_t uns; + TMP_DECL (marker); + + if (mpz_cmp_ui (u, 1) <= 0) + return 0; + + n2 = mpz_scan1 (u, 0); + if (n2 == 1) + return 0; + + TMP_MARK (marker); + + uns = ABSIZ (u) - n2 / BITS_PER_MP_LIMB; + MPZ_TMP_INIT (q, uns); + MPZ_TMP_INIT (u2, uns); + + mpz_tdiv_q_2exp (u2, u, n2); + + if (isprime (n2)) + goto n2prime; + + for (i = 1; primes[i] != 0; i++) + { + prime = primes[i]; + rem = mpz_tdiv_ui (u2, prime); + if (rem == 0) /* divisable? */ + { + rem = mpz_tdiv_q_ui (q, u2, prime * prime); + if (rem != 0) + { + TMP_FREE (marker); + return 0; + } + mpz_swap (q, u2); + for (n = 2;;) + { + rem = mpz_tdiv_q_ui (q, u2, prime); + if (rem != 0) + break; + mpz_swap (q, u2); + n++; + } + + n2 = gcd (n2, n); + if (n2 == 1) + { + TMP_FREE (marker); + return 0; + } + + /* As soon as n2 becomes a prime number, stop factoring. + Either we have u=x^n2 or u is not a perfect power. */ + if (isprime (n2)) + goto n2prime; + } + } + + if (mpz_cmp_ui (u2, 1) == 0) + { + TMP_FREE (marker); + return 1; + } + + if (n2 == 0) + { + unsigned long int nth; + /* We did not find any factors above. We have to consider all values + of n. */ + for (nth = 2;; nth++) + { + if (! isprime (nth)) + continue; +#if 0 + exact = mpz_padic_root (q, u2, nth, PTH); + if (exact) +#endif + exact = mpz_root (q, u2, nth); + if (exact) + { + TMP_FREE (marker); + return 1; + } + if (mpz_cmp_ui (q, SMALLEST_OMITTED_PRIME) < 0) + { + TMP_FREE (marker); + return 0; + } + } + } + else + { + unsigned long int nth; + /* We found some factors above. We just need to consider values of n + that divides n2. */ + for (nth = 2; nth <= n2; nth++) + { + if (! isprime (nth)) + continue; + if (n2 % nth != 0) + continue; +#if 0 + exact = mpz_padic_root (q, u2, nth, PTH); + if (exact) +#endif + exact = mpz_root (q, u2, nth); + if (exact) + { + TMP_FREE (marker); + return 1; + } + if (mpz_cmp_ui (q, SMALLEST_OMITTED_PRIME) < 0) + { + TMP_FREE (marker); + return 0; + } + } + + TMP_FREE (marker); + return 0; + } + +n2prime: + exact = mpz_root (NULL, u2, n2); + TMP_FREE (marker); + return exact; +} + +static unsigned long int +#if __STDC__ +gcd (unsigned long int a, unsigned long int b) +#else +gcd (a, b) + unsigned long int a, b; +#endif +{ + int an2, bn2, n2; + + if (a == 0) + return b; + if (b == 0) + return a; + + count_trailing_zeros (an2, a); + a >>= an2; + + count_trailing_zeros (bn2, b); + b >>= bn2; + + n2 = MIN (an2, bn2); + + while (a != b) + { + if (a > b) + { + a -= b; + do + a >>= 1; + while ((a & 1) == 0); + } + else /* b > a. */ + { + b -= a; + do + b >>= 1; + while ((b & 1) == 0); + } + } + + return a << n2; +} + +static int +#if __STDC__ +isprime (unsigned long int t) +#else +isprime (t) + unsigned long int t; +#endif +{ + unsigned long int q, r, d; + + if (t < 3 || (t & 1) == 0) + return t == 2; + + for (d = 3, r = 1; r != 0; d += 2) + { + q = t / d; + r = t - q * d; + if (q < d) + return 1; + } + return 0; +} diff --git a/ghc/rts/gmp/mpz/perfsqr.c b/ghc/rts/gmp/mpz/perfsqr.c index cdf1b5a271..92e8d08ea9 100644 --- a/ghc/rts/gmp/mpz/perfsqr.c +++ b/ghc/rts/gmp/mpz/perfsqr.c @@ -1,21 +1,21 @@ /* mpz_perfect_square_p(arg) -- Return non-zero if ARG is a perfect square, zero otherwise. -Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc. +Copyright (C) 1991, 1993, 1994, 1996, 2000 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -37,5 +37,9 @@ mpz_perfect_square_p (a) if (asize < 0) return 0; + /* Zero is a perfect square. */ + if (asize == 0) + return 1; + return mpn_perfect_square_p (a->_mp_d, asize); } diff --git a/ghc/rts/gmp/mpz/popcount.c b/ghc/rts/gmp/mpz/popcount.c index a979380309..3105258e26 100644 --- a/ghc/rts/gmp/mpz/popcount.c +++ b/ghc/rts/gmp/mpz/popcount.c @@ -6,16 +6,16 @@ Copyright (C) 1994, 1996 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/pow_ui.c b/ghc/rts/gmp/mpz/pow_ui.c index d8cf7a6106..96ca114e4d 100644 --- a/ghc/rts/gmp/mpz/pow_ui.c +++ b/ghc/rts/gmp/mpz/pow_ui.c @@ -1,20 +1,20 @@ /* mpz_pow_ui(res, base, exp) -- Set RES to BASE**EXP. -Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc. +Copyright (C) 1991, 1993, 1994, 1996, 1997 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -49,7 +49,7 @@ rpow (b, e, r) #endif /* BERKELEY_MP */ { mp_ptr rp, bp, tp, xp; - mp_size_t rsize, bsize; + mp_size_t ralloc, rsize, bsize; int cnt, i; mp_limb_t blimb; TMP_DECL (marker); @@ -82,13 +82,13 @@ rpow (b, e, r) /* Estimate space requirements accurately. Using the code from the `else' path would over-estimate space requirements wildly. */ float lb = __mp_bases[blimb].chars_per_bit_exactly; - rsize = 2 + ((mp_size_t) (e / lb) / BITS_PER_MP_LIMB); + ralloc = 3 + ((mp_size_t) (e / lb) / BITS_PER_MP_LIMB); } else { /* Over-estimate space requirements somewhat. */ count_leading_zeros (cnt, blimb); - rsize = bsize * e - cnt * e / BITS_PER_MP_LIMB + 1; + ralloc = bsize * e - cnt * e / BITS_PER_MP_LIMB + 2; } TMP_MARK (marker); @@ -97,8 +97,8 @@ rpow (b, e, r) product for mpn_mul. (This scheme is used to fulfill the requirements of mpn_mul; that the product space may not be the same as any of the input operands.) */ - rp = (mp_ptr) TMP_ALLOC (rsize * BYTES_PER_MP_LIMB); - tp = (mp_ptr) TMP_ALLOC (rsize * BYTES_PER_MP_LIMB); + rp = (mp_ptr) TMP_ALLOC (ralloc * BYTES_PER_MP_LIMB); + tp = (mp_ptr) TMP_ALLOC (ralloc * BYTES_PER_MP_LIMB); MPN_COPY (rp, bp, bsize); rsize = bsize; diff --git a/ghc/rts/gmp/mpz/powm.c b/ghc/rts/gmp/mpz/powm.c index 5dcd1b1282..e6af855a71 100644 --- a/ghc/rts/gmp/mpz/powm.c +++ b/ghc/rts/gmp/mpz/powm.c @@ -1,20 +1,21 @@ /* mpz_powm(res,base,exp,mod) -- Set RES to (base**exp) mod MOD. -Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc. +Copyright (C) 1991, 1993, 1994, 1996, 1997, 2000 Free Software Foundation, Inc. +Contributed by Paul Zimmermann. 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -22,255 +23,342 @@ MA 02111-1307, USA. */ #include "gmp.h" #include "gmp-impl.h" #include "longlong.h" +#ifdef BERKELEY_MP +#include "mp.h" +#endif + + +/* set c <- (a*b)/R^n mod m c has to have at least (2n) allocated limbs */ +static void +#if __STDC__ +mpz_redc (mpz_ptr c, mpz_srcptr a, mpz_srcptr b, mpz_srcptr m, mp_limb_t Nprim) +#else +mpz_redc (c, a, b, m, Nprim) + mpz_ptr c; + mpz_srcptr a; + mpz_srcptr b; + mpz_srcptr m; + mp_limb_t Nprim; +#endif +{ + mp_ptr cp, mp = PTR (m); + mp_limb_t cy, cout = 0; + mp_limb_t q; + size_t j, n = ABSIZ (m); + + ASSERT (ALLOC (c) >= 2 * n); + + mpz_mul (c, a, b); + cp = PTR (c); + j = ABSIZ (c); + MPN_ZERO (cp + j, 2 * n - j); + for (j = 0; j < n; j++) + { + q = cp[0] * Nprim; + cy = mpn_addmul_1 (cp, mp, n, q); + cout += mpn_add_1 (cp + n, cp + n, n - j, cy); + cp++; + } + cp -= n; + if (cout) + { + cy = cout - mpn_sub_n (cp, cp + n, mp, n); + while (cy) + cy -= mpn_sub_n (cp, cp, mp, n); + } + else + MPN_COPY (cp, cp + n, n); + MPN_NORMALIZE (cp, n); + SIZ (c) = SIZ (c) < 0 ? -n : n; +} +/* average number of calls to redc for an exponent of n bits + with the sliding window algorithm of base 2^k: the optimal is + obtained for the value of k which minimizes 2^(k-1)+n/(k+1): + + n\k 4 5 6 7 8 + 128 156* 159 171 200 261 + 256 309 307* 316 343 403 + 512 617 607* 610 632 688 + 1024 1231 1204 1195* 1207 1256 + 2048 2461 2399 2366 2360* 2396 + 4096 4918 4787 4707 4665* 4670 +*/ + #ifndef BERKELEY_MP void #if __STDC__ -mpz_powm (mpz_ptr res, mpz_srcptr base, mpz_srcptr exp, mpz_srcptr mod) +mpz_powm (mpz_ptr res, mpz_srcptr base, mpz_srcptr e, mpz_srcptr mod) #else -mpz_powm (res, base, exp, mod) +mpz_powm (res, base, e, mod) mpz_ptr res; mpz_srcptr base; - mpz_srcptr exp; + mpz_srcptr e; mpz_srcptr mod; #endif #else /* BERKELEY_MP */ void #if __STDC__ -pow (mpz_srcptr base, mpz_srcptr exp, mpz_srcptr mod, mpz_ptr res) +pow (mpz_srcptr base, mpz_srcptr e, mpz_srcptr mod, mpz_ptr res) #else -pow (base, exp, mod, res) +pow (base, e, mod, res) mpz_srcptr base; - mpz_srcptr exp; + mpz_srcptr e; mpz_srcptr mod; mpz_ptr res; #endif #endif /* BERKELEY_MP */ { - mp_ptr rp, ep, mp, bp; - mp_size_t esize, msize, bsize, rsize; - mp_size_t size; - int mod_shift_cnt; - int negative_result; - mp_limb_t *free_me = NULL; - size_t free_me_size; - TMP_DECL (marker); - - esize = ABS (exp->_mp_size); - msize = ABS (mod->_mp_size); - size = 2 * msize; + mp_limb_t invm, *ep, c, mask; + mpz_t xx, *g; + mp_size_t n, i, K, j, l, k; + int sh; + int use_redc; + +#ifdef POWM_DEBUG + mpz_t exp; + mpz_init (exp); +#endif - rp = res->_mp_d; - ep = exp->_mp_d; + n = ABSIZ (mod); - if (msize == 0) - msize = 1 / msize; /* provoke a signal */ + if (n == 0) + DIVIDE_BY_ZERO; - if (esize == 0) + if (SIZ (e) == 0) { /* Exponent is zero, result is 1 mod MOD, i.e., 1 or 0 - depending on if MOD equals 1. */ - rp[0] = 1; - res->_mp_size = (msize == 1 && (mod->_mp_d)[0] == 1) ? 0 : 1; + depending on if MOD equals 1. */ + SIZ(res) = (ABSIZ (mod) == 1 && (PTR(mod))[0] == 1) ? 0 : 1; + PTR(res)[0] = 1; return; } - TMP_MARK (marker); - - /* Normalize MOD (i.e. make its most significant bit set) as required by - mpn_divmod. This will make the intermediate values in the calculation - slightly larger, but the correct result is obtained after a final - reduction using the original MOD value. */ + /* Use REDC instead of usual reduction for sizes < POWM_THRESHOLD. + In REDC each modular multiplication costs about 2*n^2 limbs operations, + whereas using usual reduction it costs 3*K(n), where K(n) is the cost of a + multiplication using Karatsuba, and a division is assumed to cost 2*K(n), + for example using Burnikel-Ziegler's algorithm. This gives a theoretical + threshold of a*KARATSUBA_SQR_THRESHOLD, with a=(3/2)^(1/(2-ln(3)/ln(2))) ~ + 2.66. */ + /* For now, also disable REDC when MOD is even, as the inverse can't + handle that. */ + +#ifndef POWM_THRESHOLD +#define POWM_THRESHOLD ((8 * KARATSUBA_SQR_THRESHOLD) / 3) +#endif - mp = (mp_ptr) TMP_ALLOC (msize * BYTES_PER_MP_LIMB); - count_leading_zeros (mod_shift_cnt, mod->_mp_d[msize - 1]); - if (mod_shift_cnt != 0) - mpn_lshift (mp, mod->_mp_d, msize, mod_shift_cnt); - else - MPN_COPY (mp, mod->_mp_d, msize); + use_redc = (n < POWM_THRESHOLD && PTR(mod)[0] % 2 != 0); + if (use_redc) + { + /* invm = -1/m mod 2^BITS_PER_MP_LIMB, must have m odd */ + modlimb_invert (invm, PTR(mod)[0]); + invm = -invm; + } - bsize = ABS (base->_mp_size); - if (bsize > msize) + /* determines optimal value of k */ + l = ABSIZ (e) * BITS_PER_MP_LIMB; /* number of bits of exponent */ + k = 1; + K = 2; + while (2 * l > K * (2 + k * (3 + k))) { - /* The base is larger than the module. Reduce it. */ + k++; + K *= 2; + } - /* Allocate (BSIZE + 1) with space for remainder and quotient. - (The quotient is (bsize - msize + 1) limbs.) */ - bp = (mp_ptr) TMP_ALLOC ((bsize + 1) * BYTES_PER_MP_LIMB); - MPN_COPY (bp, base->_mp_d, bsize); - /* We don't care about the quotient, store it above the remainder, - at BP + MSIZE. */ - mpn_divmod (bp + msize, bp, bsize, mp, msize); - bsize = msize; - /* Canonicalize the base, since we are going to multiply with it - quite a few times. */ - MPN_NORMALIZE (bp, bsize); + g = (mpz_t *) (*_mp_allocate_func) (K / 2 * sizeof (mpz_t)); + /* compute x*R^n where R=2^BITS_PER_MP_LIMB */ + mpz_init (g[0]); + if (use_redc) + { + mpz_mul_2exp (g[0], base, n * BITS_PER_MP_LIMB); + mpz_mod (g[0], g[0], mod); } else - bp = base->_mp_d; + mpz_mod (g[0], base, mod); - if (bsize == 0) + /* compute xx^g for odd g < 2^k */ + mpz_init (xx); + if (use_redc) { - res->_mp_size = 0; - TMP_FREE (marker); - return; + _mpz_realloc (xx, 2 * n); + mpz_redc (xx, g[0], g[0], mod, invm); /* xx = x^2*R^n */ } - - if (res->_mp_alloc < size) + else { - /* We have to allocate more space for RES. If any of the input - parameters are identical to RES, defer deallocation of the old - space. */ - - if (rp == ep || rp == mp || rp == bp) + mpz_mul (xx, g[0], g[0]); + mpz_mod (xx, xx, mod); + } + for (i = 1; i < K / 2; i++) + { + mpz_init (g[i]); + if (use_redc) { - free_me = rp; - free_me_size = res->_mp_alloc; + _mpz_realloc (g[i], 2 * n); + mpz_redc (g[i], g[i - 1], xx, mod, invm); /* g[i] = x^(2i+1)*R^n */ } else - (*_mp_free_func) (rp, res->_mp_alloc * BYTES_PER_MP_LIMB); - - rp = (mp_ptr) (*_mp_allocate_func) (size * BYTES_PER_MP_LIMB); - res->_mp_alloc = size; - res->_mp_d = rp; - } - else - { - /* Make BASE, EXP and MOD not overlap with RES. */ - if (rp == bp) { - /* RES and BASE are identical. Allocate temp. space for BASE. */ - bp = (mp_ptr) TMP_ALLOC (bsize * BYTES_PER_MP_LIMB); - MPN_COPY (bp, rp, bsize); + mpz_mul (g[i], g[i - 1], xx); + mpz_mod (g[i], g[i], mod); } - if (rp == ep) + } + + /* now starts the real stuff */ + mask = (mp_limb_t) ((1<<k) - 1); + ep = PTR (e); + i = ABSIZ (e) - 1; /* current index */ + c = ep[i]; /* current limb */ + count_leading_zeros (sh, c); + sh = BITS_PER_MP_LIMB - sh; /* significant bits in ep[i] */ + sh -= k; /* index of lower bit of ep[i] to take into account */ + if (sh < 0) + { /* k-sh extra bits are needed */ + if (i > 0) { - /* RES and EXP are identical. Allocate temp. space for EXP. */ - ep = (mp_ptr) TMP_ALLOC (esize * BYTES_PER_MP_LIMB); - MPN_COPY (ep, rp, esize); + i--; + c = (c << (-sh)) | (ep[i] >> (BITS_PER_MP_LIMB + sh)); + sh += BITS_PER_MP_LIMB; } - if (rp == mp) + } + else + c = c >> sh; +#ifdef POWM_DEBUG + printf ("-1/m mod 2^%u = %lu\n", BITS_PER_MP_LIMB, invm); + mpz_set_ui (exp, c); +#endif + j=0; + while (c % 2 == 0) + { + j++; + c = (c >> 1); + } + mpz_set (xx, g[c >> 1]); + while (j--) + { + if (use_redc) + mpz_redc (xx, xx, xx, mod, invm); + else { - /* RES and MOD are identical. Allocate temporary space for MOD. */ - mp = (mp_ptr) TMP_ALLOC (msize * BYTES_PER_MP_LIMB); - MPN_COPY (mp, rp, msize); + mpz_mul (xx, xx, xx); + mpz_mod (xx, xx, mod); } } - MPN_COPY (rp, bp, bsize); - rsize = bsize; - - { - mp_size_t i; - mp_ptr xp = (mp_ptr) TMP_ALLOC (2 * (msize + 1) * BYTES_PER_MP_LIMB); - int c; - mp_limb_t e; - mp_limb_t carry_limb; - - negative_result = (ep[0] & 1) && base->_mp_size < 0; - - i = esize - 1; - e = ep[i]; - count_leading_zeros (c, e); - e = (e << c) << 1; /* shift the exp bits to the left, lose msb */ - c = BITS_PER_MP_LIMB - 1 - c; - - /* Main loop. - - Make the result be pointed to alternately by XP and RP. This - helps us avoid block copying, which would otherwise be necessary - with the overlap restrictions of mpn_divmod. With 50% probability - the result after this loop will be in the area originally pointed - by RP (==RES->_mp_d), and with 50% probability in the area originally - pointed to by XP. */ +#ifdef POWM_DEBUG + printf ("x^"); mpz_out_str (0, 10, exp); + printf ("*2^%u mod m = ", n * BITS_PER_MP_LIMB); mpz_out_str (0, 10, xx); + putchar ('\n'); +#endif - for (;;) - { - while (c != 0) - { - mp_ptr tp; - mp_size_t xsize; + while (i > 0 || sh > 0) + { + c = ep[i]; + sh -= k; + l = k; /* number of bits treated */ + if (sh < 0) + { + if (i > 0) + { + i--; + c = (c << (-sh)) | (ep[i] >> (BITS_PER_MP_LIMB + sh)); + sh += BITS_PER_MP_LIMB; + } + else + { + l += sh; /* may be less bits than k here */ + c = c & ((1<<l) - 1); + } + } + else + c = c >> sh; + c = c & mask; - mpn_mul_n (xp, rp, rp, rsize); - xsize = 2 * rsize; - if (xsize > msize) - { - mpn_divmod (xp + msize, xp, xsize, mp, msize); - xsize = msize; - } + /* this while loop implements the sliding window improvement */ + while ((c & (1 << (k - 1))) == 0 && (i > 0 || sh > 0)) + { + if (use_redc) mpz_redc (xx, xx, xx, mod, invm); + else + { + mpz_mul (xx, xx, xx); + mpz_mod (xx, xx, mod); + } + if (sh) + { + sh--; + c = (c<<1) + ((ep[i]>>sh) & 1); + } + else + { + i--; + sh = BITS_PER_MP_LIMB - 1; + c = (c<<1) + (ep[i]>>sh); + } + } - tp = rp; rp = xp; xp = tp; - rsize = xsize; +#ifdef POWM_DEBUG + printf ("l=%u c=%lu\n", l, c); + mpz_mul_2exp (exp, exp, k); + mpz_add_ui (exp, exp, c); +#endif - if ((mp_limb_signed_t) e < 0) + /* now replace xx by xx^(2^k)*x^c */ + if (c != 0) + { + j = 0; + while (c % 2 == 0) + { + j++; + c = c >> 1; + } + /* c0 = c * 2^j, i.e. xx^(2^k)*x^c = (A^(2^(k - j))*c)^(2^j) */ + l -= j; + while (l--) + if (use_redc) mpz_redc (xx, xx, xx, mod, invm); + else { - mpn_mul (xp, rp, rsize, bp, bsize); - xsize = rsize + bsize; - if (xsize > msize) - { - mpn_divmod (xp + msize, xp, xsize, mp, msize); - xsize = msize; - } - - tp = rp; rp = xp; xp = tp; - rsize = xsize; + mpz_mul (xx, xx, xx); + mpz_mod (xx, xx, mod); } - e <<= 1; - c--; - } - - i--; - if (i < 0) - break; - e = ep[i]; - c = BITS_PER_MP_LIMB; - } - - /* We shifted MOD, the modulo reduction argument, left MOD_SHIFT_CNT - steps. Adjust the result by reducing it with the original MOD. - - Also make sure the result is put in RES->_mp_d (where it already - might be, see above). */ - - if (mod_shift_cnt != 0) - { - carry_limb = mpn_lshift (res->_mp_d, rp, rsize, mod_shift_cnt); - rp = res->_mp_d; - if (carry_limb != 0) - { - rp[rsize] = carry_limb; - rsize++; - } - } - else - { - MPN_COPY (res->_mp_d, rp, rsize); - rp = res->_mp_d; - } - - if (rsize >= msize) - { - mpn_divmod (rp + msize, rp, rsize, mp, msize); - rsize = msize; - } - - /* Remove any leading zero words from the result. */ - if (mod_shift_cnt != 0) - mpn_rshift (rp, rp, rsize, mod_shift_cnt); - MPN_NORMALIZE (rp, rsize); - } + if (use_redc) + mpz_redc (xx, xx, g[c >> 1], mod, invm); + else + { + mpz_mul (xx, xx, g[c >> 1]); + mpz_mod (xx, xx, mod); + } + } + else + j = l; /* case c=0 */ + while (j--) + { + if (use_redc) + mpz_redc (xx, xx, xx, mod, invm); + else + { + mpz_mul (xx, xx, xx); + mpz_mod (xx, xx, mod); + } + } +#ifdef POWM_DEBUG + printf ("x^"); mpz_out_str (0, 10, exp); + printf ("*2^%u mod m = ", n * BITS_PER_MP_LIMB); mpz_out_str (0, 10, xx); + putchar ('\n'); +#endif + } - if (negative_result && rsize != 0) + /* now convert back xx to xx/R^n */ + if (use_redc) { - if (mod_shift_cnt != 0) - mpn_rshift (mp, mp, msize, mod_shift_cnt); - mpn_sub (rp, mp, msize, rp, rsize); - rsize = msize; - MPN_NORMALIZE (rp, rsize); + mpz_set_ui (g[0], 1); + mpz_redc (xx, xx, g[0], mod, invm); + if (mpz_cmp (xx, mod) >= 0) + mpz_sub (xx, xx, mod); } - res->_mp_size = rsize; + mpz_set (res, xx); - if (free_me != NULL) - (*_mp_free_func) (free_me, free_me_size * BYTES_PER_MP_LIMB); - TMP_FREE (marker); + mpz_clear (xx); + for (i = 0; i < K / 2; i++) + mpz_clear (g[i]); + (*_mp_free_func) (g, K / 2 * sizeof (mpz_t)); } diff --git a/ghc/rts/gmp/mpz/powm_ui.c b/ghc/rts/gmp/mpz/powm_ui.c index 596815a0ea..00f70bd563 100644 --- a/ghc/rts/gmp/mpz/powm_ui.c +++ b/ghc/rts/gmp/mpz/powm_ui.c @@ -1,24 +1,26 @@ /* mpz_powm_ui(res,base,exp,mod) -- Set RES to (base**exp) mod MOD. -Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc. +Copyright (C) 1991, 1993, 1994, 1996, 1997, 2000 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ +#include <stdio.h> /* for NULL */ #include "gmp.h" #include "gmp-impl.h" #include "longlong.h" @@ -49,12 +51,14 @@ mpz_powm_ui (res, base, exp, mod) rp = res->_mp_d; if (msize == 0) - msize = 1 / msize; /* provoke a signal */ + DIVIDE_BY_ZERO; if (exp == 0) { - rp[0] = 1; + /* Exponent is zero, result is 1 mod MOD, i.e., 1 or 0 + depending on if MOD equals 1. */ res->_mp_size = (msize == 1 && (mod->_mp_d)[0] == 1) ? 0 : 1; + rp[0] = 1; return; } @@ -166,6 +170,7 @@ mpz_powm_ui (res, base, exp, mod) mpn_mul_n (xp, rp, rp, rsize); xsize = 2 * rsize; + xsize -= xp[xsize - 1] == 0; if (xsize > msize) { mpn_divmod (xp + msize, xp, xsize, mp, msize); @@ -179,6 +184,7 @@ mpz_powm_ui (res, base, exp, mod) { mpn_mul (xp, rp, rsize, bp, bsize); xsize = rsize + bsize; + xsize -= xp[xsize - 1] == 0; if (xsize > msize) { mpn_divmod (xp + msize, xp, xsize, mp, msize); @@ -226,7 +232,15 @@ mpz_powm_ui (res, base, exp, mod) MPN_NORMALIZE (rp, rsize); } - res->_mp_size = negative_result == 0 ? rsize : -rsize; + if (negative_result && rsize != 0) + { + if (mod_shift_cnt != 0) + mpn_rshift (mp, mp, msize, mod_shift_cnt); + mpn_sub (rp, mp, msize, rp, rsize); + rsize = msize; + MPN_NORMALIZE (rp, rsize); + } + res->_mp_size = rsize; if (free_me != NULL) (*_mp_free_func) (free_me, free_me_size * BYTES_PER_MP_LIMB); diff --git a/ghc/rts/gmp/mpz/pprime_p.c b/ghc/rts/gmp/mpz/pprime_p.c index 494de14cae..0217d8f5f0 100644 --- a/ghc/rts/gmp/mpz/pprime_p.c +++ b/ghc/rts/gmp/mpz/pprime_p.c @@ -2,114 +2,241 @@ An implementation of the probabilistic primality test found in Knuth's Seminumerical Algorithms book. If the function mpz_probab_prime_p() returns 0 then n is not prime. If it returns 1, then n is 'probably' - prime. The probability of a false positive is (1/4)**reps, where - reps is the number of internal passes of the probabilistic algorithm. - Knuth indicates that 25 passes are reasonable. + prime. If it returns 2, n is surely prime. The probability of a false + positive is (1/4)**reps, where reps is the number of internal passes of the + probabilistic algorithm. Knuth indicates that 25 passes are reasonable. -Copyright (C) 1991, 1993, 1994 Free Software Foundation, Inc. -Contributed by John Amanatides. +Copyright (C) 1991, 1993, 1994, 1996, 1997, 1998, 1999, 2000 Free Software +Foundation, Inc. Miller-Rabin code contributed by John Amanatides. 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ #include "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" -static int -possibly_prime (n, n_minus_1, x, y, q, k) +static int isprime _PROTO ((unsigned long int t)); +static int mpz_millerrabin _PROTO ((mpz_srcptr n, int reps)); + +int +#if __STDC__ +mpz_probab_prime_p (mpz_srcptr n, int reps) +#else +mpz_probab_prime_p (n, reps) mpz_srcptr n; - mpz_srcptr n_minus_1; - mpz_ptr x; - mpz_ptr y; - mpz_srcptr q; - unsigned long int k; + int reps; +#endif { - unsigned long int i; + mp_limb_t r; - /* find random x s.t. 1 < x < n */ - do + /* Handle small and negative n. */ + if (mpz_cmp_ui (n, 1000000L) <= 0) { - mpz_random (x, mpz_size (n)); - mpz_mmod (x, x, n); + int is_prime; + if (mpz_sgn (n) < 0) + { + /* Negative number. Negate and call ourselves. */ + mpz_t n2; + mpz_init (n2); + mpz_neg (n2, n); + is_prime = mpz_probab_prime_p (n2, reps); + mpz_clear (n2); + return is_prime; + } + is_prime = isprime (mpz_get_ui (n)); + return is_prime ? 2 : 0; } - while (mpz_cmp_ui (x, 1L) <= 0); - mpz_powm (y, x, q, n); + /* If n is now even, it is not a prime. */ + if ((mpz_get_ui (n) & 1) == 0) + return 0; + + /* Check if n has small factors. */ + if (UDIV_TIME > (2 * UMUL_TIME + 6)) + r = mpn_preinv_mod_1 (PTR(n), SIZ(n), (mp_limb_t) PP, (mp_limb_t) PP_INVERTED); + else + r = mpn_mod_1 (PTR(n), SIZ(n), (mp_limb_t) PP); + if (r % 3 == 0 || r % 5 == 0 || r % 7 == 0 || r % 11 == 0 || r % 13 == 0 + || r % 17 == 0 || r % 19 == 0 || r % 23 == 0 || r % 29 == 0 +#if BITS_PER_MP_LIMB == 64 + || r % 31 == 0 || r % 37 == 0 || r % 41 == 0 || r % 43 == 0 + || r % 47 == 0 || r % 53 == 0 +#endif + ) + { + return 0; + } - if (mpz_cmp_ui (y, 1L) == 0 || mpz_cmp (y, n_minus_1) == 0) - return 1; + /* Do more dividing. We collect small primes, using umul_ppmm, until we + overflow a single limb. We divide our number by the small primes product, + and look for factors in the remainder. */ + { + unsigned long int ln2; + unsigned long int q; + mp_limb_t p1, p0, p; + unsigned int primes[15]; + int nprimes; + + nprimes = 0; + p = 1; + ln2 = mpz_sizeinbase (n, 2) / 30; ln2 = ln2 * ln2; + for (q = BITS_PER_MP_LIMB == 64 ? 59 : 31; q < ln2; q += 2) + { + if (isprime (q)) + { + umul_ppmm (p1, p0, p, q); + if (p1 != 0) + { + r = mpn_mod_1 (PTR(n), SIZ(n), p); + while (--nprimes >= 0) + if (r % primes[nprimes] == 0) + { + if (mpn_mod_1 (PTR(n), SIZ(n), (mp_limb_t) primes[nprimes]) != 0) + abort (); + return 0; + } + p = q; + nprimes = 0; + } + else + { + p = p0; + } + primes[nprimes++] = q; + } + } + } + + /* Perform a number of Miller-Rabin tests. */ + return mpz_millerrabin (n, reps); +} - for (i = 1; i < k; i++) +static int +#if __STDC__ +isprime (unsigned long int t) +#else +isprime (t) + unsigned long int t; +#endif +{ + unsigned long int q, r, d; + + if (t < 3 || (t & 1) == 0) + return t == 2; + + for (d = 3, r = 1; r != 0; d += 2) { - mpz_powm_ui (y, y, 2L, n); - if (mpz_cmp (y, n_minus_1) == 0) + q = t / d; + r = t - q * d; + if (q < d) return 1; - if (mpz_cmp_ui (y, 1L) == 0) - return 0; } return 0; } -int +static int millerrabin _PROTO ((mpz_srcptr n, mpz_srcptr nm1, + mpz_ptr x, mpz_ptr y, + mpz_srcptr q, unsigned long int k)); + +static int #if __STDC__ -mpz_probab_prime_p (mpz_srcptr m, int reps) +mpz_millerrabin (mpz_srcptr n, int reps) #else -mpz_probab_prime_p (m, reps) - mpz_srcptr m; +mpz_millerrabin (n, reps) + mpz_srcptr n; int reps; #endif { - mpz_t n, n_minus_1, x, y, q; - int i, is_prime; + int r; + mpz_t nm1, x, y, q; unsigned long int k; + gmp_randstate_t rstate; + int is_prime; + TMP_DECL (marker); + TMP_MARK (marker); - mpz_init (n); - /* Take the absolute value of M, to handle positive and negative primes. */ - mpz_abs (n, m); + MPZ_TMP_INIT (nm1, SIZ (n) + 1); + mpz_sub_ui (nm1, n, 1L); - if (mpz_cmp_ui (n, 3L) <= 0) - { - mpz_clear (n); - return mpz_cmp_ui (n, 1L) > 0; - } + MPZ_TMP_INIT (x, SIZ (n)); + MPZ_TMP_INIT (y, 2 * SIZ (n)); /* mpz_powm_ui needs excessive memory!!! */ - if ((mpz_get_ui (n) & 1) == 0) + /* Perform a Fermat test. */ + mpz_set_ui (x, 210L); + mpz_powm (y, x, nm1, n); + if (mpz_cmp_ui (y, 1L) != 0) { - mpz_clear (n); - return 0; /* even */ + return 0; + TMP_FREE (marker); } - mpz_init (n_minus_1); - mpz_sub_ui (n_minus_1, n, 1L); - mpz_init (x); - mpz_init (y); + MPZ_TMP_INIT (q, SIZ (n)); + + /* Find q and k, where q is odd and n = 1 + 2**k * q. */ + k = mpz_scan1 (nm1, 0L); + mpz_tdiv_q_2exp (q, nm1, k); - /* find q and k, s.t. n = 1 + 2**k * q */ - mpz_init_set (q, n_minus_1); - k = mpz_scan1 (q, 0); - mpz_tdiv_q_2exp (q, q, k); + gmp_randinit (rstate, GMP_RAND_ALG_DEFAULT, 32L); is_prime = 1; - for (i = 0; i < reps && is_prime; i++) - is_prime &= possibly_prime (n, n_minus_1, x, y, q, k); - - mpz_clear (n_minus_1); - mpz_clear (n); - mpz_clear (x); - mpz_clear (y); - mpz_clear (q); + for (r = 0; r < reps && is_prime; r++) + { + do + mpz_urandomb (x, rstate, mpz_sizeinbase (n, 2) - 1); + while (mpz_cmp_ui (x, 1L) <= 0); + + is_prime = millerrabin (n, nm1, x, y, q, k); + } + + gmp_randclear (rstate); + + TMP_FREE (marker); return is_prime; } + +static int +#if __STDC__ +millerrabin (mpz_srcptr n, mpz_srcptr nm1, mpz_ptr x, mpz_ptr y, + mpz_srcptr q, unsigned long int k) +#else +millerrabin (n, nm1, x, y, q, k) + mpz_srcptr n; + mpz_srcptr nm1; + mpz_ptr x; + mpz_ptr y; + mpz_srcptr q; + unsigned long int k; +#endif +{ + unsigned long int i; + + mpz_powm (y, x, q, n); + + if (mpz_cmp_ui (y, 1L) == 0 || mpz_cmp (y, nm1) == 0) + return 1; + + for (i = 1; i < k; i++) + { + mpz_powm_ui (y, y, 2L, n); + if (mpz_cmp (y, nm1) == 0) + return 1; + if (mpz_cmp_ui (y, 1L) == 0) + return 0; + } + return 0; +} diff --git a/ghc/rts/gmp/mpz/random.c b/ghc/rts/gmp/mpz/random.c index ab41eef926..60d9113991 100644 --- a/ghc/rts/gmp/mpz/random.c +++ b/ghc/rts/gmp/mpz/random.c @@ -6,16 +6,16 @@ Copyright (C) 1991, 1993, 1994, 1995 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/random2.c b/ghc/rts/gmp/mpz/random2.c index 725a8b4914..a90af115e9 100644 --- a/ghc/rts/gmp/mpz/random2.c +++ b/ghc/rts/gmp/mpz/random2.c @@ -7,16 +7,16 @@ Copyright (C) 1991, 1993, 1994, 1996 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/realloc.c b/ghc/rts/gmp/mpz/realloc.c index 2c2a5da927..0b9e447ec3 100644 --- a/ghc/rts/gmp/mpz/realloc.c +++ b/ghc/rts/gmp/mpz/realloc.c @@ -5,16 +5,16 @@ Copyright (C) 1991, 1993, 1994, 1995 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/remove.c b/ghc/rts/gmp/mpz/remove.c new file mode 100644 index 0000000000..bc6675f972 --- /dev/null +++ b/ghc/rts/gmp/mpz/remove.c @@ -0,0 +1,93 @@ +/* mpz_remove -- divide out a factor and return its multiplicity. + +Copyright (C) 1998, 1999, 2000 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 2.1 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; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +unsigned long int +#if __STDC__ +mpz_remove (mpz_ptr dest, mpz_srcptr src, mpz_srcptr f) +#else +mpz_remove (dest, src, f) + mpz_ptr dest; + mpz_srcptr src; + mpz_srcptr f; +#endif +{ + mpz_t fpow[40]; /* inexhaustible...until year 2020 or so */ + mpz_t x, rem; + unsigned long int pwr; + int p; + + if (mpz_cmp_ui (f, 1) <= 0 || mpz_sgn (src) == 0) + DIVIDE_BY_ZERO; + if (mpz_cmp_ui (f, 2) == 0) + { + unsigned long int s0; + s0 = mpz_scan1 (src, 0); + mpz_div_2exp (dest, src, s0); + return s0; + } + + /* We could perhaps compute mpz_scan1(src,0)/mpz_scan1(f,0). It is an + upper bound of the result we're seeking. We could also shift down the + operands so that they become odd, to make intermediate values smaller. */ + + mpz_init (rem); + mpz_init (x); + + pwr = 0; + mpz_init (fpow[0]); + mpz_set (fpow[0], f); + mpz_set (dest, src); + + /* Divide by f, f^2, ..., f^(2^k) until we get a remainder for f^(2^k). */ + for (p = 0;; p++) + { + mpz_tdiv_qr (x, rem, dest, fpow[p]); + if (SIZ (rem) != 0) + break; + mpz_init (fpow[p + 1]); + mpz_mul (fpow[p + 1], fpow[p], fpow[p]); + mpz_set (dest, x); + } + + pwr = (1 << p) - 1; + + mpz_clear (fpow[p]); + + /* Divide by f^(2^(k-1)), f^(2^(k-2)), ..., f for all divisors that give a + zero remainder. */ + while (--p >= 0) + { + mpz_tdiv_qr (x, rem, dest, fpow[p]); + if (SIZ (rem) == 0) + { + pwr += 1 << p; + mpz_set (dest, x); + } + mpz_clear (fpow[p]); + } + + mpz_clear (x); + mpz_clear (rem); + return pwr; +} diff --git a/ghc/rts/gmp/mpz/root.c b/ghc/rts/gmp/mpz/root.c new file mode 100644 index 0000000000..0920bf22d3 --- /dev/null +++ b/ghc/rts/gmp/mpz/root.c @@ -0,0 +1,183 @@ +/* mpz_root(root, u, nth) -- Set ROOT to floor(U^(1/nth)). + Return an indication if the result is exact. + +Copyright (C) 1999, 2000 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 2.1 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; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +/* Naive implementation of nth root extraction. It would probably be a + better idea to use a division-free Newton iteration. It is insane + to use full precision from iteration 1. The mpz_scan1 trick compensates + to some extent. It would be natural to avoid representing the low zero + bits mpz_scan1 is counting, and at the same time call mpn directly. */ + +#include <stdio.h> /* for NULL */ +#include "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" + +int +#if __STDC__ +mpz_root (mpz_ptr r, mpz_srcptr c, unsigned long int nth) +#else +mpz_root (r, c, nth) + mpz_ptr r; + mpz_srcptr c; + unsigned long int nth; +#endif +{ + mpz_t x, t0, t1, t2; + __mpz_struct ccs, *cc = &ccs; + unsigned long int nbits; + int bit; + int exact; + int i; + unsigned long int lowz; + unsigned long int rl; + + /* even roots of negatives provoke an exception */ + if (mpz_sgn (c) < 0 && (nth & 1) == 0) + SQRT_OF_NEGATIVE; + + /* root extraction interpreted as c^(1/nth) means a zeroth root should + provoke a divide by zero, do this even if c==0 */ + if (nth == 0) + DIVIDE_BY_ZERO; + + if (mpz_sgn (c) == 0) + { + if (r != NULL) + mpz_set_ui (r, 0); + return 1; /* exact result */ + } + + PTR(cc) = PTR(c); + SIZ(cc) = ABSIZ(c); + + nbits = (mpz_sizeinbase (cc, 2) - 1) / nth; + if (nbits == 0) + { + if (r != NULL) + mpz_set_ui (r, 1); + if (mpz_sgn (c) < 0) + { + if (r != NULL) + SIZ(r) = -SIZ(r); + return mpz_cmp_si (c, -1L) == 0; + } + return mpz_cmp_ui (c, 1L) == 0; + } + + mpz_init (x); + mpz_init (t0); + mpz_init (t1); + mpz_init (t2); + + /* Create a one-bit approximation. */ + mpz_set_ui (x, 0); + mpz_setbit (x, nbits); + + /* Make the approximation better, one bit at a time. This odd-looking + termination criteria makes large nth get better initial approximation, + which avoids slow convergence for such values. */ + bit = nbits - 1; + for (i = 1; (nth >> i) != 0; i++) + { + mpz_setbit (x, bit); + mpz_tdiv_q_2exp (t0, x, bit); + mpz_pow_ui (t1, t0, nth); + mpz_mul_2exp (t1, t1, bit * nth); + if (mpz_cmp (cc, t1) < 0) + mpz_clrbit (x, bit); + + bit--; /* check/set next bit */ + if (bit < 0) + { + /* We're done. */ + mpz_pow_ui (t1, x, nth); + goto done; + } + } + mpz_setbit (x, bit); + mpz_set_ui (t2, 0); mpz_setbit (t2, bit); mpz_add (x, x, t2); + +#if DEBUG + /* Check that the starting approximation is >= than the root. */ + mpz_pow_ui (t1, x, nth); + if (mpz_cmp (cc, t1) >= 0) + abort (); +#endif + + mpz_add_ui (x, x, 1); + + /* Main loop */ + do + { + lowz = mpz_scan1 (x, 0); + mpz_tdiv_q_2exp (t0, x, lowz); + mpz_pow_ui (t1, t0, nth - 1); + mpz_mul_2exp (t1, t1, lowz * (nth - 1)); + mpz_tdiv_q (t2, cc, t1); + mpz_sub (t2, x, t2); + rl = mpz_tdiv_q_ui (t2, t2, nth); + mpz_sub (x, x, t2); + } + while (mpz_sgn (t2) != 0); + + /* If we got a non-zero remainder in the last division, we know our root + is too large. */ + mpz_sub_ui (x, x, (mp_limb_t) (rl != 0)); + + /* Adjustment loop. If we spend more care on rounding in the loop above, + we could probably get rid of this, or greatly simplify it. */ + { + int bad = 0; + lowz = mpz_scan1 (x, 0); + mpz_tdiv_q_2exp (t0, x, lowz); + mpz_pow_ui (t1, t0, nth); + mpz_mul_2exp (t1, t1, lowz * nth); + while (mpz_cmp (cc, t1) < 0) + { + bad++; + if (bad > 2) + abort (); /* abort if our root is far off */ + mpz_sub_ui (x, x, 1); + lowz = mpz_scan1 (x, 0); + mpz_tdiv_q_2exp (t0, x, lowz); + mpz_pow_ui (t1, t0, nth); + mpz_mul_2exp (t1, t1, lowz * nth); + } + } + + done: + exact = mpz_cmp (t1, cc) == 0; + + if (r != NULL) + { + mpz_set (r, x); + if (mpz_sgn (c) < 0) + SIZ(r) = -SIZ(r); + } + + mpz_clear (t2); + mpz_clear (t1); + mpz_clear (t0); + mpz_clear (x); + + return exact; +} diff --git a/ghc/rts/gmp/mpz/rrandomb.c b/ghc/rts/gmp/mpz/rrandomb.c new file mode 100644 index 0000000000..7d78243674 --- /dev/null +++ b/ghc/rts/gmp/mpz/rrandomb.c @@ -0,0 +1,117 @@ +/* mpz_rrandomb -- Generate a positive random mpz_t of specified bit size, with + long runs of consecutive ones and zeros in the binary representation. + Meant for testing of other MP routines. + +Copyright (C) 2000 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 2.1 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; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +static void gmp_rrandomb _PROTO ((mp_ptr rp, gmp_randstate_t rstate, unsigned long int nbits)); + +void +#if __STDC__ +mpz_rrandomb (mpz_ptr x, gmp_randstate_t rstate, unsigned long int nbits) +#else +mpz_rrandomb (x, rstate, nbits) + mpz_ptr x; + gmp_randstate_t rstate; + unsigned long int nbits; +#endif +{ + mp_size_t nl = 0; + + if (nbits != 0) + { + mp_ptr xp; + nl = (nbits + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB; + if (x->_mp_alloc < nl) + _mpz_realloc (x, nl); + + xp = PTR(x); + gmp_rrandomb (xp, rstate, nbits); + MPN_NORMALIZE (xp, nl); + } + + SIZ(x) = nl; +} + +#define BITS_PER_CHUNK 4 + +static void +#if __STDC__ +gmp_rrandomb (mp_ptr rp, gmp_randstate_t rstate, unsigned long int nbits) +#else +gmp_rrandomb (rp, rstate, nbits) + mp_ptr rp; + gmp_randstate_t rstate; + unsigned long int nbits; +#endif +{ + int nb; + int bit_pos; + mp_size_t limb_pos; + mp_limb_t ran, ranm; + mp_limb_t acc; + mp_size_t n; + + bit_pos = nbits % BITS_PER_MP_LIMB; + limb_pos = nbits / BITS_PER_MP_LIMB; + if (bit_pos == 0) + { + bit_pos = BITS_PER_MP_LIMB; + limb_pos--; + } + + acc = 0; + while (limb_pos >= 0) + { + _gmp_rand (&ranm, rstate, BITS_PER_CHUNK + 1); + ran = ranm; + nb = (ran >> 1) + 1; + if ((ran & 1) != 0) + { + /* Generate a string of ones. */ + if (nb > bit_pos) + { + rp[limb_pos--] = acc | ((((mp_limb_t) 1) << bit_pos) - 1); + bit_pos += BITS_PER_MP_LIMB; + bit_pos -= nb; + acc = (~(mp_limb_t) 0) << bit_pos; + } + else + { + bit_pos -= nb; + acc |= ((((mp_limb_t) 1) << nb) - 1) << bit_pos; + } + } + else + { + /* Generate a string of zeroes. */ + if (nb > bit_pos) + { + rp[limb_pos--] = acc; + acc = 0; + bit_pos += BITS_PER_MP_LIMB; + } + bit_pos -= nb; + } + } +} diff --git a/ghc/rts/gmp/mpz/scan0.c b/ghc/rts/gmp/mpz/scan0.c index 8e45aa3248..6c59cf8939 100644 --- a/ghc/rts/gmp/mpz/scan0.c +++ b/ghc/rts/gmp/mpz/scan0.c @@ -5,16 +5,16 @@ Copyright (C) 1996 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/scan1.c b/ghc/rts/gmp/mpz/scan1.c index 6ecb0aab38..3b84e3420c 100644 --- a/ghc/rts/gmp/mpz/scan1.c +++ b/ghc/rts/gmp/mpz/scan1.c @@ -5,16 +5,16 @@ Copyright (C) 1996 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/set.c b/ghc/rts/gmp/mpz/set.c index d94ab7485b..06b2eef511 100644 --- a/ghc/rts/gmp/mpz/set.c +++ b/ghc/rts/gmp/mpz/set.c @@ -5,16 +5,16 @@ Copyright (C) 1991, 1993, 1994, 1995 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/set_d.c b/ghc/rts/gmp/mpz/set_d.c index c09b1d97e6..e90ed9bc2f 100644 --- a/ghc/rts/gmp/mpz/set_d.c +++ b/ghc/rts/gmp/mpz/set_d.c @@ -1,20 +1,20 @@ /* mpz_set_d(integer, val) -- Assign INTEGER with a double value VAL. -Copyright (C) 1995 Free Software Foundation, Inc. +Copyright (C) 1995, 1996, 2000 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -32,9 +32,9 @@ mpz_set_d (r, d) #endif { int negative; - mp_size_t size; mp_limb_t tp[3]; mp_ptr rp; + mp_size_t rn; negative = d < 0; d = ABS (d); @@ -49,19 +49,20 @@ mpz_set_d (r, d) return; } - size = __gmp_extract_double (tp, d); + rn = __gmp_extract_double (tp, d); - if (ALLOC(r) < size) - _mpz_realloc (r, size); + if (ALLOC(r) < rn) + _mpz_realloc (r, rn); rp = PTR (r); #if BITS_PER_MP_LIMB == 32 - switch (size) + switch (rn) { default: - MPN_ZERO (rp, size - 3); - rp += size - 3; + MPN_ZERO (rp, rn - 3); + rp += rn - 3; + /* fall through */ case 3: rp[2] = tp[2]; rp[1] = tp[1]; @@ -72,22 +73,24 @@ mpz_set_d (r, d) rp[0] = tp[1]; break; case 1: + /* handled in "small aguments" case above */ abort (); } #else - switch (size) + switch (rn) { default: - MPN_ZERO (rp, size - 2); - rp += size - 2; + MPN_ZERO (rp, rn - 2); + rp += rn - 2; + /* fall through */ case 2: - rp[1] = tp[1]; - rp[0] = tp[0]; + rp[1] = tp[1], rp[0] = tp[0]; break; case 1: + /* handled in "small aguments" case above */ abort (); } #endif - SIZ(r) = negative ? -size : size; + SIZ(r) = negative ? -rn : rn; } diff --git a/ghc/rts/gmp/mpz/set_f.c b/ghc/rts/gmp/mpz/set_f.c index 9547907bb3..2273953dfd 100644 --- a/ghc/rts/gmp/mpz/set_f.c +++ b/ghc/rts/gmp/mpz/set_f.c @@ -5,16 +5,16 @@ Copyright (C) 1996 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/set_q.c b/ghc/rts/gmp/mpz/set_q.c index 61bd5c70d9..72d3222a80 100644 --- a/ghc/rts/gmp/mpz/set_q.c +++ b/ghc/rts/gmp/mpz/set_q.c @@ -6,16 +6,16 @@ Copyright (C) 1996 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/set_si.c b/ghc/rts/gmp/mpz/set_si.c index 82a90a3cc1..9ba2fbaf30 100644 --- a/ghc/rts/gmp/mpz/set_si.c +++ b/ghc/rts/gmp/mpz/set_si.c @@ -1,20 +1,20 @@ /* mpz_set_si(integer, val) -- Assign INTEGER with a small value VAL. -Copyright (C) 1991, 1993, 1994, 1995 Free Software Foundation, Inc. +Copyright (C) 1991, 1993, 1994, 1995, 2000 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -40,7 +40,7 @@ mpz_set_si (dest, val) } else if (val < 0) { - dest->_mp_d[0] = -val; + dest->_mp_d[0] = (unsigned long) -val; dest->_mp_size = -1; } else diff --git a/ghc/rts/gmp/mpz/set_str.c b/ghc/rts/gmp/mpz/set_str.c index d1334b1e00..3ab79c0e89 100644 --- a/ghc/rts/gmp/mpz/set_str.c +++ b/ghc/rts/gmp/mpz/set_str.c @@ -4,34 +4,40 @@ the base in the C standard way, i.e. 0xhh...h means base 16, 0oo...o means base 8, otherwise assume base 10. -Copyright (C) 1991, 1993, 1994, Free Software Foundation, Inc. +Copyright (C) 1991, 1993, 1994, 1996, 1997, 1998, 2000 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ +#include <string.h> #include <ctype.h> #include "gmp.h" #include "gmp-impl.h" #include "longlong.h" static int +#if __STDC__ +digit_value_in_base (int c, int base) +#else digit_value_in_base (c, base) int c; int base; +#endif { int digit; @@ -96,13 +102,30 @@ mpz_set_str (x, str, base) base = 16; c = *str++; } + else if (c == 'b' || c == 'B') + { + base = 2; + c = *str++; + } } } + /* Skip leading zeros. */ + while (c == '0') + c = *str++; + /* Make sure the string does not become empty, mpn_set_str would fail. */ + if (c == 0) + { + x->_mp_size = 0; + return 0; + } + TMP_MARK (marker); str_size = strlen (str - 1); s = begs = (char *) TMP_ALLOC (str_size + 1); + /* Remove spaces from the string and convert the result from ASCII to a + byte array. */ for (i = 0; i < str_size; i++) { if (!isspace (c)) @@ -120,10 +143,12 @@ mpz_set_str (x, str, base) str_size = s - begs; - xsize = str_size / __mp_bases[base].chars_per_limb + 1; + xsize = (((mp_size_t) (str_size / __mp_bases[base].chars_per_bit_exactly)) + / BITS_PER_MP_LIMB + 2); if (x->_mp_alloc < xsize) _mpz_realloc (x, xsize); + /* Convert the byte array in base BASE to our bignum format. */ xsize = mpn_set_str (x->_mp_d, (unsigned char *) begs, str_size, base); x->_mp_size = negative ? -xsize : xsize; diff --git a/ghc/rts/gmp/mpz/set_ui.c b/ghc/rts/gmp/mpz/set_ui.c index 73f6cf54ad..d6097c170a 100644 --- a/ghc/rts/gmp/mpz/set_ui.c +++ b/ghc/rts/gmp/mpz/set_ui.c @@ -5,16 +5,16 @@ Copyright (C) 1991, 1993, 1994, 1995 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/setbit.c b/ghc/rts/gmp/mpz/setbit.c index af59e2cc78..d4249a434e 100644 --- a/ghc/rts/gmp/mpz/setbit.c +++ b/ghc/rts/gmp/mpz/setbit.c @@ -1,20 +1,21 @@ /* mpz_setbit -- set a specified bit. -Copyright (C) 1991, 1993, 1994, 1995 Free Software Foundation, Inc. +Copyright (C) 1991, 1993, 1994, 1995, 1997, 1999 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -63,7 +64,7 @@ mpz_setbit (d, bit_index) /* Simulate two's complement arithmetic, i.e. simulate 1. Set OP = ~(OP - 1) [with infinitely many leading ones]. - 2. set the bit. + 2. Set the bit. 3. Set OP = ~OP + 1. */ dsize = -dsize; @@ -108,6 +109,11 @@ mpz_setbit (d, bit_index) } } else - ; + { + mpn_decr_u (dp + limb_index, + (mp_limb_t) 1 << (bit_index % BITS_PER_MP_LIMB)); + dsize -= dp[dsize - 1] == 0; + d->_mp_size = -dsize; + } } } diff --git a/ghc/rts/gmp/mpz/size.c b/ghc/rts/gmp/mpz/size.c index 0b09fbeb07..6574756783 100644 --- a/ghc/rts/gmp/mpz/size.c +++ b/ghc/rts/gmp/mpz/size.c @@ -6,16 +6,16 @@ Copyright (C) 1991, 1993, 1994, 1995 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/sizeinbase.c b/ghc/rts/gmp/mpz/sizeinbase.c index 51bd5558ac..734f9c4532 100644 --- a/ghc/rts/gmp/mpz/sizeinbase.c +++ b/ghc/rts/gmp/mpz/sizeinbase.c @@ -7,16 +7,16 @@ Copyright (C) 1991, 1993, 1994, 1995 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/sqrt.c b/ghc/rts/gmp/mpz/sqrt.c index 44c554eb5c..fe82fe407a 100644 --- a/ghc/rts/gmp/mpz/sqrt.c +++ b/ghc/rts/gmp/mpz/sqrt.c @@ -1,24 +1,25 @@ /* mpz_sqrt(root, u) -- Set ROOT to floor(sqrt(U)). -Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc. +Copyright (C) 1991, 1993, 1994, 1996, 2000 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ +#include <stdio.h> /* for NULL */ #include "gmp.h" #include "gmp-impl.h" @@ -40,7 +41,7 @@ mpz_sqrt (root, op) TMP_MARK (marker); op_size = op->_mp_size; if (op_size < 0) - op_size = 1 / (op_size > 0); /* Divide by zero for negative OP. */ + SQRT_OF_NEGATIVE; /* The size of the root is accurate after this simple calculation. */ root_size = (op_size + 1) / 2; diff --git a/ghc/rts/gmp/mpz/sqrtrem.c b/ghc/rts/gmp/mpz/sqrtrem.c index 757cc5dd56..99a6453122 100644 --- a/ghc/rts/gmp/mpz/sqrtrem.c +++ b/ghc/rts/gmp/mpz/sqrtrem.c @@ -1,27 +1,31 @@ /* mpz_sqrtrem(root,rem,x) -- Set ROOT to floor(sqrt(X)) and REM to the remainder, i.e. X - ROOT**2. -Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc. +Copyright (C) 1991, 1993, 1994, 1996, 2000 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ +#include <stdio.h> /* for NULL */ #include "gmp.h" #include "gmp-impl.h" +#ifdef BERKELEY_MP +#include "mp.h" +#endif #ifndef BERKELEY_MP void @@ -54,7 +58,7 @@ msqrt (op, root, rem) TMP_MARK (marker); op_size = op->_mp_size; if (op_size < 0) - op_size = 1 / (op_size > 0); /* Divide by zero for negative OP. */ + SQRT_OF_NEGATIVE; if (rem->_mp_alloc < op_size) _mpz_realloc (rem, op_size); diff --git a/ghc/rts/gmp/mpz/sub.c b/ghc/rts/gmp/mpz/sub.c index 56ef1a12a2..f3ae7c23a0 100644 --- a/ghc/rts/gmp/mpz/sub.c +++ b/ghc/rts/gmp/mpz/sub.c @@ -1,26 +1,29 @@ /* mpz_sub -- Subtract two integers. -Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc. +Copyright (C) 1991, 1993, 1994, 1996, 2000 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ #include "gmp.h" #include "gmp-impl.h" +#ifdef BERKELEY_MP +#include "mp.h" +#endif #ifndef BERKELEY_MP void @@ -58,9 +61,9 @@ msub (u, v, w) if (abs_usize < abs_vsize) { /* Swap U and V. */ - {const __mpz_struct *t = u; u = v; v = t;} - {mp_size_t t = usize; usize = vsize; vsize = t;} - {mp_size_t t = abs_usize; abs_usize = abs_vsize; abs_vsize = t;} + MPZ_SRCPTR_SWAP (u, v); + MP_SIZE_T_SWAP (usize, vsize); + MP_SIZE_T_SWAP (abs_usize, abs_vsize); } /* True: ABS_USIZE >= ABS_VSIZE. */ diff --git a/ghc/rts/gmp/mpz/sub_ui.c b/ghc/rts/gmp/mpz/sub_ui.c index 7dea4b6c03..327add8503 100644 --- a/ghc/rts/gmp/mpz/sub_ui.c +++ b/ghc/rts/gmp/mpz/sub_ui.c @@ -1,20 +1,20 @@ /* mpz_sub_ui -- Subtract an unsigned one-word integer from an MP_INT. -Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc. +Copyright (C) 1991, 1993, 1994, 1996, 1999 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -59,7 +59,7 @@ mpz_sub_ui (w, u, v) if (usize < 0) { mp_limb_t cy; - cy = mpn_add_1 (wp, up, abs_usize, v); + cy = mpn_add_1 (wp, up, abs_usize, (mp_limb_t) v); wp[abs_usize] = cy; wsize = -(abs_usize + cy); } @@ -74,7 +74,7 @@ mpz_sub_ui (w, u, v) } else { - mpn_sub_1 (wp, up, abs_usize, v); + mpn_sub_1 (wp, up, abs_usize, (mp_limb_t) v); /* Size can decrease with at most one limb. */ wsize = abs_usize - (wp[abs_usize - 1] == 0); } diff --git a/ghc/rts/gmp/mpz/swap.c b/ghc/rts/gmp/mpz/swap.c new file mode 100644 index 0000000000..0070d6ff24 --- /dev/null +++ b/ghc/rts/gmp/mpz/swap.c @@ -0,0 +1,52 @@ +/* mpz_swap (dest_integer, src_integer) -- Swap U and V. + +Copyright (C) 1997, 1998 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 2.1 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; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +#if __STDC__ +mpz_swap (mpz_ptr u, mpz_ptr v) +#else +mpz_swap (u, v) + mpz_ptr u; + mpz_ptr v; +#endif +{ + mp_ptr up, vp; + mp_size_t usize, vsize; + mp_size_t ualloc, valloc; + + ualloc = u->_mp_alloc; + valloc = v->_mp_alloc; + v->_mp_alloc = ualloc; + u->_mp_alloc = valloc; + + usize = u->_mp_size; + vsize = v->_mp_size; + v->_mp_size = usize; + u->_mp_size = vsize; + + up = u->_mp_d; + vp = v->_mp_d; + v->_mp_d = up; + u->_mp_d = vp; +} diff --git a/ghc/rts/gmp/mpz/tdiv_q.c b/ghc/rts/gmp/mpz/tdiv_q.c index b4d36368b8..21db4ab385 100644 --- a/ghc/rts/gmp/mpz/tdiv_q.c +++ b/ghc/rts/gmp/mpz/tdiv_q.c @@ -1,20 +1,20 @@ /* mpz_tdiv_q -- divide two integers and produce a quotient. -Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc. +Copyright (C) 1991, 1993, 1994, 1996, 2000 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -33,101 +33,59 @@ mpz_tdiv_q (quot, num, den) mpz_srcptr den; #endif { - mp_srcptr np, dp; - mp_ptr qp, rp; - mp_size_t nsize = num->_mp_size; - mp_size_t dsize = den->_mp_size; - mp_size_t qsize, rsize; - mp_size_t sign_quotient = nsize ^ dsize; - unsigned normalization_steps; - mp_limb_t q_limb; + mp_size_t ql; + mp_size_t ns, ds, nl, dl; + mp_ptr np, dp, qp, rp; TMP_DECL (marker); - nsize = ABS (nsize); - dsize = ABS (dsize); + ns = SIZ (num); + ds = SIZ (den); + nl = ABS (ns); + dl = ABS (ds); + ql = nl - dl + 1; - /* Ensure space is enough for quotient. */ + if (dl == 0) + DIVIDE_BY_ZERO; - qsize = nsize - dsize + 1; /* qsize cannot be bigger than this. */ - if (qsize <= 0) + if (ql <= 0) { - quot->_mp_size = 0; + SIZ (quot) = 0; return; } - if (quot->_mp_alloc < qsize) - _mpz_realloc (quot, qsize); - - qp = quot->_mp_d; - np = num->_mp_d; - dp = den->_mp_d; - - /* Optimize division by a single-limb divisor. */ - if (dsize == 1) - { - mpn_divmod_1 (qp, np, nsize, dp[0]); - qsize -= qp[qsize - 1] == 0; - quot->_mp_size = sign_quotient >= 0 ? qsize : -qsize; - return; - } + MPZ_REALLOC (quot, ql); TMP_MARK (marker); + qp = PTR (quot); + rp = (mp_ptr) TMP_ALLOC (dl * BYTES_PER_MP_LIMB); + np = PTR (num); + dp = PTR (den); - rp = (mp_ptr) TMP_ALLOC ((nsize + 1) * BYTES_PER_MP_LIMB); - - count_leading_zeros (normalization_steps, dp[dsize - 1]); + /* FIXME: We should think about how to handle the temporary allocation. + Perhaps mpn_tdiv_qr should handle it, since it anyway often needs to + allocate temp space. */ - /* Normalize the denominator, i.e. make its most significant bit set by - shifting it NORMALIZATION_STEPS bits to the left. Also shift the - numerator the same number of steps (to keep the quotient the same!). */ - if (normalization_steps != 0) + /* Copy denominator to temporary space if it overlaps with the quotient. */ + if (dp == qp) { mp_ptr tp; - mp_limb_t nlimb; - - /* Shift up the denominator setting the most significant bit of - the most significant word. Use temporary storage not to clobber - the original contents of the denominator. */ - tp = (mp_ptr) TMP_ALLOC (dsize * BYTES_PER_MP_LIMB); - mpn_lshift (tp, dp, dsize, normalization_steps); + tp = (mp_ptr) TMP_ALLOC (dl * BYTES_PER_MP_LIMB); + MPN_COPY (tp, dp, dl); dp = tp; - - /* Shift up the numerator, possibly introducing a new most - significant word. Move the shifted numerator in the remainder - meanwhile. */ - nlimb = mpn_lshift (rp, np, nsize, normalization_steps); - if (nlimb != 0) - { - rp[nsize] = nlimb; - rsize = nsize + 1; - } - else - rsize = nsize; } - else + /* Copy numerator to temporary space if it overlaps with the quotient. */ + if (np == qp) { - /* The denominator is already normalized, as required. Copy it to - temporary space if it overlaps with the quotient. */ - if (dp == qp) - { - dp = (mp_ptr) TMP_ALLOC (dsize * BYTES_PER_MP_LIMB); - MPN_COPY ((mp_ptr) dp, qp, dsize); - } - - /* Move the numerator to the remainder. */ - MPN_COPY (rp, np, nsize); - rsize = nsize; + mp_ptr tp; + tp = (mp_ptr) TMP_ALLOC (nl * BYTES_PER_MP_LIMB); + MPN_COPY (tp, np, nl); + np = tp; } - q_limb = mpn_divmod (qp, rp, rsize, dp, dsize); + mpn_tdiv_qr (qp, rp, 0L, np, nl, dp, dl); - qsize = rsize - dsize; - if (q_limb) - { - qp[qsize] = q_limb; - qsize += 1; - } + ql -= qp[ql - 1] == 0; - quot->_mp_size = sign_quotient >= 0 ? qsize : -qsize; + SIZ (quot) = (ns ^ ds) >= 0 ? ql : -ql; TMP_FREE (marker); } diff --git a/ghc/rts/gmp/mpz/tdiv_q_2exp.c b/ghc/rts/gmp/mpz/tdiv_q_2exp.c index e70d810e7f..03d1e01f89 100644 --- a/ghc/rts/gmp/mpz/tdiv_q_2exp.c +++ b/ghc/rts/gmp/mpz/tdiv_q_2exp.c @@ -6,16 +6,16 @@ Copyright (C) 1991, 1993, 1994, 1996 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/tdiv_q_ui.c b/ghc/rts/gmp/mpz/tdiv_q_ui.c index 9048e0aef5..a2e3462b76 100644 --- a/ghc/rts/gmp/mpz/tdiv_q_ui.c +++ b/ghc/rts/gmp/mpz/tdiv_q_ui.c @@ -1,21 +1,21 @@ /* mpz_tdiv_q_ui(quot, dividend, divisor_limb) -- Divide DIVIDEND by DIVISOR_LIMB and store the result in QUOT. -Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc. +Copyright (C) 1991, 1993, 1994, 1996, 1998 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -23,7 +23,7 @@ MA 02111-1307, USA. */ #include "gmp.h" #include "gmp-impl.h" -void +unsigned long int #if __STDC__ mpz_tdiv_q_ui (mpz_ptr quot, mpz_srcptr dividend, unsigned long int divisor) #else @@ -36,16 +36,14 @@ mpz_tdiv_q_ui (quot, dividend, divisor) mp_size_t dividend_size; mp_size_t size; mp_ptr quot_ptr; + mp_limb_t remainder_limb; + + if (divisor == 0) + DIVIDE_BY_ZERO; dividend_size = dividend->_mp_size; size = ABS (dividend_size); - if (size == 0) - { - quot->_mp_size = 0; - return; - } - /* No need for temporary allocation and copying if QUOT == DIVIDEND as the divisor is just one limb, and thus no intermediate remainders need to be stored. */ @@ -55,9 +53,12 @@ mpz_tdiv_q_ui (quot, dividend, divisor) quot_ptr = quot->_mp_d; - mpn_divmod_1 (quot_ptr, dividend->_mp_d, size, (mp_limb_t) divisor); + remainder_limb + = mpn_divmod_1 (quot_ptr, dividend->_mp_d, size, (mp_limb_t) divisor); /* The quotient is SIZE limbs, but the most significant might be zero. */ - size -= quot_ptr[size - 1] == 0; + size -= size != 0 && quot_ptr[size - 1] == 0; quot->_mp_size = dividend_size >= 0 ? size : -size; + + return remainder_limb; } diff --git a/ghc/rts/gmp/mpz/tdiv_qr.c b/ghc/rts/gmp/mpz/tdiv_qr.c index 500e19939a..d66f57d9e5 100644 --- a/ghc/rts/gmp/mpz/tdiv_qr.c +++ b/ghc/rts/gmp/mpz/tdiv_qr.c @@ -1,21 +1,21 @@ /* mpz_tdiv_qr(quot,rem,dividend,divisor) -- Set QUOT to DIVIDEND/DIVISOR, and REM to DIVIDEND mod DIVISOR. -Copyright (C) 1991, 1993, 1994 Free Software Foundation, Inc. +Copyright (C) 1991, 1993, 1994, 2000 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -23,6 +23,12 @@ MA 02111-1307, USA. */ #include "gmp.h" #include "gmp-impl.h" #include "longlong.h" +#ifdef BERKELEY_MP +#include "mp.h" +#endif + + +#ifndef BERKELEY_MP void #if __STDC__ @@ -35,5 +41,90 @@ mpz_tdiv_qr (quot, rem, num, den) mpz_srcptr den; #endif -#define COMPUTE_QUOTIENT -#include "dmincl.c" +#else /* BERKELEY_MP */ + +void +#if __STDC__ +mdiv (mpz_srcptr num, mpz_srcptr den, mpz_ptr quot, mpz_ptr rem) +#else +mdiv (num, den, quot, rem) + mpz_srcptr num; + mpz_srcptr den; + mpz_ptr quot; + mpz_ptr rem; +#endif + +#endif /* BERKELEY_MP */ +{ + mp_size_t ql; + mp_size_t ns, ds, nl, dl; + mp_ptr np, dp, qp, rp; + TMP_DECL (marker); + + ns = SIZ (num); + ds = SIZ (den); + nl = ABS (ns); + dl = ABS (ds); + ql = nl - dl + 1; + + if (dl == 0) + DIVIDE_BY_ZERO; + + MPZ_REALLOC (rem, dl); + + if (ql <= 0) + { + if (num != rem) + { + mp_ptr np, rp; + np = PTR (num); + rp = PTR (rem); + MPN_COPY (rp, np, nl); + SIZ (rem) = SIZ (num); + } + /* This needs to follow the assignment to rem, in case the + numerator and quotient are the same. */ + SIZ (quot) = 0; + return; + } + + MPZ_REALLOC (quot, ql); + + TMP_MARK (marker); + qp = PTR (quot); + rp = PTR (rem); + np = PTR (num); + dp = PTR (den); + + /* FIXME: We should think about how to handle the temporary allocation. + Perhaps mpn_tdiv_qr should handle it, since it anyway often needs to + allocate temp space. */ + + /* Copy denominator to temporary space if it overlaps with the quotient + or remainder. */ + if (dp == rp || dp == qp) + { + mp_ptr tp; + tp = (mp_ptr) TMP_ALLOC (dl * BYTES_PER_MP_LIMB); + MPN_COPY (tp, dp, dl); + dp = tp; + } + /* Copy numerator to temporary space if it overlaps with the quotient or + remainder. */ + if (np == rp || np == qp) + { + mp_ptr tp; + tp = (mp_ptr) TMP_ALLOC (nl * BYTES_PER_MP_LIMB); + MPN_COPY (tp, np, nl); + np = tp; + } + + mpn_tdiv_qr (qp, rp, 0L, np, nl, dp, dl); + + ql -= qp[ql - 1] == 0; + MPN_NORMALIZE (rp, dl); + + SIZ (quot) = (ns ^ ds) >= 0 ? ql : -ql; + SIZ (rem) = ns >= 0 ? dl : -dl; + TMP_FREE (marker); +} diff --git a/ghc/rts/gmp/mpz/tdiv_qr_ui.c b/ghc/rts/gmp/mpz/tdiv_qr_ui.c index cb5041cdc5..10368cd340 100644 --- a/ghc/rts/gmp/mpz/tdiv_qr_ui.c +++ b/ghc/rts/gmp/mpz/tdiv_qr_ui.c @@ -2,21 +2,21 @@ Set QUOT to DIVIDEND / SHORT_DIVISOR and REM to DIVIDEND mod SHORT_DIVISOR. -Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc. +Copyright (C) 1991, 1993, 1994, 1996, 1998 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -24,7 +24,7 @@ MA 02111-1307, USA. */ #include "gmp.h" #include "gmp-impl.h" -void +unsigned long int #if __STDC__ mpz_tdiv_qr_ui (mpz_ptr quot, mpz_ptr rem, mpz_srcptr dividend, unsigned long int divisor) #else @@ -40,16 +40,12 @@ mpz_tdiv_qr_ui (quot, rem, dividend, divisor) mp_ptr quot_ptr; mp_limb_t remainder_limb; + if (divisor == 0) + DIVIDE_BY_ZERO; + dividend_size = dividend->_mp_size; size = ABS (dividend_size); - if (size == 0) - { - quot->_mp_size = 0; - rem->_mp_size = 0; - return; - } - /* No need for temporary allocation and copying if QUOT == DIVIDEND as the divisor is just one limb, and thus no intermediate remainders need to be stored. */ @@ -60,7 +56,7 @@ mpz_tdiv_qr_ui (quot, rem, dividend, divisor) quot_ptr = quot->_mp_d; remainder_limb = mpn_divmod_1 (quot_ptr, dividend->_mp_d, size, - (mp_limb_t) divisor); + (mp_limb_t) divisor); if (remainder_limb == 0) rem->_mp_size = 0; @@ -73,6 +69,8 @@ mpz_tdiv_qr_ui (quot, rem, dividend, divisor) } /* The quotient is SIZE limbs, but the most significant might be zero. */ - size -= quot_ptr[size - 1] == 0; + size -= size != 0 && quot_ptr[size - 1] == 0; quot->_mp_size = dividend_size >= 0 ? size : -size; + + return remainder_limb; } diff --git a/ghc/rts/gmp/mpz/tdiv_r.c b/ghc/rts/gmp/mpz/tdiv_r.c index 813a0d4d79..9eb87dfabf 100644 --- a/ghc/rts/gmp/mpz/tdiv_r.c +++ b/ghc/rts/gmp/mpz/tdiv_r.c @@ -1,20 +1,20 @@ /* mpz_tdiv_r(rem, dividend, divisor) -- Set REM to DIVIDEND mod DIVISOR. -Copyright (C) 1991, 1993, 1994 Free Software Foundation, Inc. +Copyright (C) 1991, 1993, 1994, 2000 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -32,6 +32,67 @@ mpz_tdiv_r (rem, num, den) mpz_srcptr num; mpz_srcptr den; #endif +{ + mp_size_t ql; + mp_size_t ns, ds, nl, dl; + mp_ptr np, dp, qp, rp; + TMP_DECL (marker); -#undef COMPUTE_QUOTIENT -#include "dmincl.c" + ns = SIZ (num); + ds = SIZ (den); + nl = ABS (ns); + dl = ABS (ds); + ql = nl - dl + 1; + + if (dl == 0) + DIVIDE_BY_ZERO; + + MPZ_REALLOC (rem, dl); + + if (ql <= 0) + { + if (num != rem) + { + mp_ptr np, rp; + np = PTR (num); + rp = PTR (rem); + MPN_COPY (rp, np, nl); + SIZ (rem) = SIZ (num); + } + return; + } + + TMP_MARK (marker); + qp = (mp_ptr) TMP_ALLOC (ql * BYTES_PER_MP_LIMB); + rp = PTR (rem); + np = PTR (num); + dp = PTR (den); + + /* FIXME: We should think about how to handle the temporary allocation. + Perhaps mpn_tdiv_qr should handle it, since it anyway often needs to + allocate temp space. */ + + /* Copy denominator to temporary space if it overlaps with the remainder. */ + if (dp == rp) + { + mp_ptr tp; + tp = (mp_ptr) TMP_ALLOC (dl * BYTES_PER_MP_LIMB); + MPN_COPY (tp, dp, dl); + dp = tp; + } + /* Copy numerator to temporary space if it overlaps with the remainder. */ + if (np == rp) + { + mp_ptr tp; + tp = (mp_ptr) TMP_ALLOC (nl * BYTES_PER_MP_LIMB); + MPN_COPY (tp, np, nl); + np = tp; + } + + mpn_tdiv_qr (qp, rp, 0L, np, nl, dp, dl); + + MPN_NORMALIZE (rp, dl); + + SIZ (rem) = ns >= 0 ? dl : -dl; + TMP_FREE (marker); +} diff --git a/ghc/rts/gmp/mpz/tdiv_r_2exp.c b/ghc/rts/gmp/mpz/tdiv_r_2exp.c index 99d617e9f3..91de170f5c 100644 --- a/ghc/rts/gmp/mpz/tdiv_r_2exp.c +++ b/ghc/rts/gmp/mpz/tdiv_r_2exp.c @@ -5,16 +5,16 @@ Copyright (C) 1991, 1993, 1994, 1995 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ diff --git a/ghc/rts/gmp/mpz/tdiv_r_ui.c b/ghc/rts/gmp/mpz/tdiv_r_ui.c index 0428b52c3a..2ea411fda1 100644 --- a/ghc/rts/gmp/mpz/tdiv_r_ui.c +++ b/ghc/rts/gmp/mpz/tdiv_r_ui.c @@ -1,21 +1,21 @@ /* mpz_tdiv_r_ui(rem, dividend, divisor_limb) -- Set REM to DIVDEND mod DIVISOR_LIMB. -Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc. +Copyright (C) 1991, 1993, 1994, 1996, 1998 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -23,7 +23,7 @@ MA 02111-1307, USA. */ #include "gmp.h" #include "gmp-impl.h" -void +unsigned long int #if __STDC__ mpz_tdiv_r_ui (mpz_ptr rem, mpz_srcptr dividend, unsigned long int divisor) #else @@ -37,15 +37,12 @@ mpz_tdiv_r_ui (rem, dividend, divisor) mp_size_t size; mp_limb_t remainder_limb; + if (divisor == 0) + DIVIDE_BY_ZERO; + dividend_size = dividend->_mp_size; size = ABS (dividend_size); - if (size == 0) - { - rem->_mp_size = 0; - return; - } - /* No need for temporary allocation and copying if QUOT == DIVIDEND as the divisor is just one limb, and thus no intermediate remainders need to be stored. */ @@ -61,4 +58,6 @@ mpz_tdiv_r_ui (rem, dividend, divisor) rem->_mp_size = dividend_size >= 0 ? 1 : -1; rem->_mp_d[0] = remainder_limb; } + + return remainder_limb; } diff --git a/ghc/rts/gmp/mpz/tdiv_ui.c b/ghc/rts/gmp/mpz/tdiv_ui.c new file mode 100644 index 0000000000..7a40a6a7f7 --- /dev/null +++ b/ghc/rts/gmp/mpz/tdiv_ui.c @@ -0,0 +1,53 @@ +/* mpz_tdiv_ui(dividend, divisor_limb) + -- Return DIVDEND mod DIVISOR_LIMB. + +Copyright (C) 1991, 1993, 1994, 1996, 1997, 1998 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 2.1 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; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +unsigned long int +#if __STDC__ +mpz_tdiv_ui (mpz_srcptr dividend, unsigned long int divisor) +#else +mpz_tdiv_ui (dividend, divisor) + mpz_srcptr dividend; + unsigned long int divisor; +#endif +{ + mp_size_t dividend_size; + mp_size_t size; + mp_limb_t remainder_limb; + + if (divisor == 0) + DIVIDE_BY_ZERO; + + dividend_size = dividend->_mp_size; + size = ABS (dividend_size); + + /* No need for temporary allocation and copying if QUOT == DIVIDEND as + the divisor is just one limb, and thus no intermediate remainders + need to be stored. */ + + remainder_limb = mpn_mod_1 (dividend->_mp_d, size, (mp_limb_t) divisor); + + return remainder_limb; +} diff --git a/ghc/rts/gmp/mpz/tstbit.c b/ghc/rts/gmp/mpz/tstbit.c new file mode 100644 index 0000000000..b0a8b0b31a --- /dev/null +++ b/ghc/rts/gmp/mpz/tstbit.c @@ -0,0 +1,70 @@ +/* mpz_tstbit -- test a specified bit. Simulate 2's complement representation. + +Copyright (C) 1997 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 2.1 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; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +int +#if __STDC__ +mpz_tstbit (mpz_srcptr d, unsigned long int bit_index) +#else +mpz_tstbit (d, bit_index) + mpz_srcptr d; + unsigned long int bit_index; +#endif +{ + mp_size_t dsize = d->_mp_size; + mp_ptr dp = d->_mp_d; + mp_size_t limb_index; + + limb_index = bit_index / BITS_PER_MP_LIMB; + if (dsize >= 0) + { + if (limb_index < dsize) + return (dp[limb_index] >> (bit_index % BITS_PER_MP_LIMB)) & 1; + else + /* Testing a bit outside of a positive number. */ + return 0; + } + else + { + mp_size_t zero_bound; + + dsize = -dsize; + + /* Locate the least significant non-zero limb. */ + for (zero_bound = 0; dp[zero_bound] == 0; zero_bound++) + ; + + if (limb_index > zero_bound) + { + if (limb_index < dsize) + return (~dp[limb_index] >> (bit_index % BITS_PER_MP_LIMB)) & 1; + else + /* Testing a bit outside of a negative number. */ + return 1; + } + else if (limb_index == zero_bound) + return (-dp[limb_index] >> (bit_index % BITS_PER_MP_LIMB)) & 1; + else + return 0; + } +} diff --git a/ghc/rts/gmp/mpz/ui_pow_ui.c b/ghc/rts/gmp/mpz/ui_pow_ui.c index 19baca14b4..edd2dee625 100644 --- a/ghc/rts/gmp/mpz/ui_pow_ui.c +++ b/ghc/rts/gmp/mpz/ui_pow_ui.c @@ -1,20 +1,21 @@ /* mpz_ui_pow_ui(res, base, exp) -- Set RES to BASE**EXP. -Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc. +Copyright (C) 1991, 1993, 1994, 1996, 1997, 2000 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 Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -23,6 +24,9 @@ MA 02111-1307, USA. */ #include "gmp-impl.h" #include "longlong.h" + +static void mpz_pow2 _PROTO ((mpz_ptr r, mp_limb_t blimb, unsigned long int e, mp_limb_t rl)); + void #if __STDC__ mpz_ui_pow_ui (mpz_ptr r, unsigned long int b, unsigned long int e) @@ -33,79 +37,103 @@ mpz_ui_pow_ui (r, b, e) unsigned long int e; #endif { - mp_ptr rp, tp, xp; - mp_size_t rsize; - int cnt, i; mp_limb_t blimb = b; - TMP_DECL (marker); + mp_limb_t rl; - /* Single out cases that give result == 0 or 1. These tests are here - to simplify the general code below, not to optimize. */ if (e == 0) { + /* For x^0 we return 1, even if x is 0. */ r->_mp_d[0] = 1; r->_mp_size = 1; return; } - if (blimb == 0) - { - r->_mp_size = 0; - return; - } - if (blimb < 0x100) + /* Compute b^e as (b^n)^(e div n) * b^(e mod n), where n is chosen such that + the latter factor is the largest number small enough to fit in a limb. */ + + rl = 1; + while (e != 0 && blimb < ((mp_limb_t) 1 << BITS_PER_MP_LIMB/2)) { - /* Estimate space requirements accurately. Using the code from the - `else' path would over-estimate space requirements wildly. */ - float lb = __mp_bases[blimb].chars_per_bit_exactly; - rsize = 2 + ((mp_size_t) (e / lb) / BITS_PER_MP_LIMB); + if ((e & 1) != 0) + rl = rl * blimb; + blimb = blimb * blimb; + e = e >> 1; } - else + + /* rl is now b^(e mod n). (I.e., the latter factor above.) */ + + if (e == 0) { - /* Over-estimate space requirements somewhat. */ - count_leading_zeros (cnt, blimb); - rsize = e - cnt * e / BITS_PER_MP_LIMB + 1; + r->_mp_d[0] = rl; + r->_mp_size = rl != 0; + return; } + mpz_pow2 (r, blimb, e, rl); +} + +/* Multi-precision part of expontialization code. */ +static void +#if __STDC__ +mpz_pow2 (mpz_ptr r, mp_limb_t blimb, unsigned long int e, mp_limb_t rl) +#else +mpz_pow2 (r, blimb, e, rl) + mpz_ptr r; + mp_limb_t blimb; + unsigned long int e; + mp_limb_t rl; +#endif +{ + mp_ptr rp, tp; + mp_size_t ralloc, rsize; + int cnt, i; + TMP_DECL (marker); + TMP_MARK (marker); - /* The two areas are used to alternatingly hold the input and recieve the - product for mpn_mul. (This scheme is used to fulfill the requirements - of mpn_mul; that the product space may not be the same as any of the - input operands.) */ - rp = (mp_ptr) TMP_ALLOC (rsize * BYTES_PER_MP_LIMB); - tp = (mp_ptr) TMP_ALLOC (rsize * BYTES_PER_MP_LIMB); + /* Over-estimate temporary space requirements somewhat. */ + count_leading_zeros (cnt, blimb); + ralloc = e - cnt * e / BITS_PER_MP_LIMB + 1; + + /* The two areas are used to alternatingly hold the input and receive the + product for mpn_mul. (Needed since mpn_mul_n requires that the product + is distinct from either input operand.) */ + rp = (mp_ptr) TMP_ALLOC (ralloc * BYTES_PER_MP_LIMB); + tp = (mp_ptr) TMP_ALLOC (ralloc * BYTES_PER_MP_LIMB); rp[0] = blimb; rsize = 1; - count_leading_zeros (cnt, e); + count_leading_zeros (cnt, e); for (i = BITS_PER_MP_LIMB - cnt - 2; i >= 0; i--) { mpn_mul_n (tp, rp, rp, rsize); rsize = 2 * rsize; rsize -= tp[rsize - 1] == 0; - xp = tp; tp = rp; rp = xp; + MP_PTR_SWAP (rp, tp); if ((e & ((mp_limb_t) 1 << i)) != 0) { mp_limb_t cy; - cy = mpn_mul_1 (tp, rp, rsize, blimb); - if (cy != 0) - { - tp[rsize] = cy; - rsize++; - } - xp = tp; tp = rp; rp = xp; + cy = mpn_mul_1 (rp, rp, rsize, blimb); + rp[rsize] = cy; + rsize += cy != 0; } } - /* Now then we know the exact space requirements, reallocate if - necessary. */ - if (r->_mp_alloc < rsize) - _mpz_realloc (r, rsize); + /* We will need rsize or rsize+1 limbs for the result. */ + if (r->_mp_alloc <= rsize) + _mpz_realloc (r, rsize + 1); + + /* Multiply the two factors (in rp,rsize and rl) and put the final result + in place. */ + { + mp_limb_t cy; + cy = mpn_mul_1 (r->_mp_d, rp, rsize, rl); + (r->_mp_d)[rsize] = cy; + rsize += cy != 0; + } - MPN_COPY (r->_mp_d, rp, rsize); r->_mp_size = rsize; TMP_FREE (marker); } diff --git a/ghc/rts/gmp/mpz/urandomb.c b/ghc/rts/gmp/mpz/urandomb.c new file mode 100644 index 0000000000..caca086e05 --- /dev/null +++ b/ghc/rts/gmp/mpz/urandomb.c @@ -0,0 +1,49 @@ +/* mpz_urandomb (rop, state, n) -- Generate a uniform pseudorandom + integer in the range 0 to 2^N - 1, inclusive, using STATE as the + random state previously initialized by a call to gmp_randinit(). + +Copyright (C) 1999, 2000 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 2.1 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; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +#if __STDC__ +mpz_urandomb (mpz_t rop, gmp_randstate_t rstate, unsigned long int nbits) +#else +mpz_urandomb (rop, rstate, nbits) + mpz_t rop; + gmp_randstate_t rstate; + unsigned long int nbits; +#endif +{ + mp_ptr rp; + mp_size_t size; + + size = (nbits + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB; + if (ALLOC (rop) < size) + _mpz_realloc (rop, size); + + rp = PTR (rop); + + _gmp_rand (rp, rstate, nbits); + MPN_NORMALIZE (rp, size); + SIZ (rop) = size; +} diff --git a/ghc/rts/gmp/mpz/urandomm.c b/ghc/rts/gmp/mpz/urandomm.c new file mode 100644 index 0000000000..aa57784081 --- /dev/null +++ b/ghc/rts/gmp/mpz/urandomm.c @@ -0,0 +1,73 @@ +/* mpz_urandomm (rop, state, n) -- Generate a uniform pseudorandom + integer in the range 0 to N-1, using STATE as the random state + previously initialized by a call to gmp_randinit(). + +Copyright (C) 2000 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 2.1 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; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" + +void +#if __STDC__ +mpz_urandomm (mpz_t rop, gmp_randstate_t rstate, mpz_t n) +#else +mpz_urandomm (rop, rstate, n) + mpz_t rop; + gmp_randstate_t rstate; + mpz_t n; +#endif +{ + mpz_t t, p, m; + mp_ptr tp; + mp_size_t nbits, size; + int count; + + /* FIXME: Should check for n == 0 and report error */ + + size = SIZ (n); + count_leading_zeros (count, PTR (n)[size - 1]); + nbits = size * BITS_PER_MP_LIMB - count; + + /* Allocate enough for any mpz function called since a realloc of + these will fail. */ + MPZ_TMP_INIT (t, size); + MPZ_TMP_INIT (m, size + 1); + MPZ_TMP_INIT (p, size + 1); + + /* Let m = highest possible random number plus 1. */ + mpz_set_ui (m, 0); + mpz_setbit (m, nbits); + + /* Let p = floor(m / n) * n. */ + mpz_fdiv_q (p, m, n); + mpz_mul (p, p, n); + + tp = PTR (t); + do + { + _gmp_rand (tp, rstate, nbits); + MPN_NORMALIZE (tp, size); /* FIXME: Really necessary? */ + SIZ (t) = size; + } + while (mpz_cmp (t, p) >= 0); + + mpz_mod (rop, t, n); +} diff --git a/ghc/rts/gmp/mpz/xor.c b/ghc/rts/gmp/mpz/xor.c new file mode 100644 index 0000000000..69898d1791 --- /dev/null +++ b/ghc/rts/gmp/mpz/xor.c @@ -0,0 +1,217 @@ +/* mpz_xor -- Logical xor. + +Copyright (C) 1991, 1993, 1994, 1996, 1997, 2000 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 2.1 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; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +#if __STDC__ +mpz_xor (mpz_ptr res, mpz_srcptr op1, mpz_srcptr op2) +#else +mpz_xor (res, op1, op2) + mpz_ptr res; + mpz_srcptr op1; + mpz_srcptr op2; +#endif +{ + mp_srcptr op1_ptr, op2_ptr; + mp_size_t op1_size, op2_size; + mp_ptr res_ptr; + mp_size_t res_size, res_alloc; + mp_size_t i; + TMP_DECL (marker); + + TMP_MARK (marker); + op1_size = op1->_mp_size; + op2_size = op2->_mp_size; + + op1_ptr = op1->_mp_d; + op2_ptr = op2->_mp_d; + res_ptr = res->_mp_d; + + if (op1_size >= 0) + { + if (op2_size >= 0) + { + if (op1_size >= op2_size) + { + if (res->_mp_alloc < op1_size) + { + _mpz_realloc (res, op1_size); + op1_ptr = op1->_mp_d; + op2_ptr = op2->_mp_d; + res_ptr = res->_mp_d; + } + + if (res_ptr != op1_ptr) + MPN_COPY (res_ptr + op2_size, op1_ptr + op2_size, + op1_size - op2_size); + for (i = op2_size - 1; i >= 0; i--) + res_ptr[i] = op1_ptr[i] ^ op2_ptr[i]; + res_size = op1_size; + } + else + { + if (res->_mp_alloc < op2_size) + { + _mpz_realloc (res, op2_size); + op1_ptr = op1->_mp_d; + op2_ptr = op2->_mp_d; + res_ptr = res->_mp_d; + } + + if (res_ptr != op2_ptr) + MPN_COPY (res_ptr + op1_size, op2_ptr + op1_size, + op2_size - op1_size); + for (i = op1_size - 1; i >= 0; i--) + res_ptr[i] = op1_ptr[i] ^ op2_ptr[i]; + res_size = op2_size; + } + + MPN_NORMALIZE (res_ptr, res_size); + res->_mp_size = res_size; + return; + } + else /* op2_size < 0 */ + { + /* Fall through to the code at the end of the function. */ + } + } + else + { + if (op2_size < 0) + { + mp_ptr opx; + mp_limb_t cy; + + /* Both operands are negative, the result will be positive. + (-OP1) ^ (-OP2) = + = ~(OP1 - 1) ^ ~(OP2 - 1) = + = (OP1 - 1) ^ (OP2 - 1) */ + + op1_size = -op1_size; + op2_size = -op2_size; + + /* Possible optimization: Decrease mpn_sub precision, + as we won't use the entire res of both. */ + opx = (mp_ptr) TMP_ALLOC (op1_size * BYTES_PER_MP_LIMB); + mpn_sub_1 (opx, op1_ptr, op1_size, (mp_limb_t) 1); + op1_ptr = opx; + + opx = (mp_ptr) TMP_ALLOC (op2_size * BYTES_PER_MP_LIMB); + mpn_sub_1 (opx, op2_ptr, op2_size, (mp_limb_t) 1); + op2_ptr = opx; + + res_alloc = MAX (op1_size, op2_size); + if (res->_mp_alloc < res_alloc) + { + _mpz_realloc (res, res_alloc); + res_ptr = res->_mp_d; + /* Don't re-read OP1_PTR and OP2_PTR. They point to + temporary space--never to the space RES->_mp_d used + to point to before reallocation. */ + } + + if (op1_size > op2_size) + { + MPN_COPY (res_ptr + op2_size, op1_ptr + op2_size, + op1_size - op2_size); + for (i = op2_size - 1; i >= 0; i--) + res_ptr[i] = op1_ptr[i] ^ op2_ptr[i]; + res_size = op1_size; + } + else + { + MPN_COPY (res_ptr + op1_size, op2_ptr + op1_size, + op2_size - op1_size); + for (i = op1_size - 1; i >= 0; i--) + res_ptr[i] = op1_ptr[i] ^ op2_ptr[i]; + res_size = op2_size; + } + + MPN_NORMALIZE (res_ptr, res_size); + res->_mp_size = res_size; + TMP_FREE (marker); + return; + } + else + { + /* We should compute -OP1 ^ OP2. Swap OP1 and OP2 and fall + through to the code that handles OP1 ^ -OP2. */ + MPZ_SRCPTR_SWAP (op1, op2); + MPN_SRCPTR_SWAP (op1_ptr,op1_size, op2_ptr,op2_size); + } + } + + { + mp_ptr opx; + mp_limb_t cy; + mp_size_t count; + + /* Operand 2 negative, so will be the result. + -(OP1 ^ (-OP2)) = -(OP1 ^ ~(OP2 - 1)) = + = ~(OP1 ^ ~(OP2 - 1)) + 1 = + = (OP1 ^ (OP2 - 1)) + 1 */ + + op2_size = -op2_size; + + opx = (mp_ptr) TMP_ALLOC (op2_size * BYTES_PER_MP_LIMB); + mpn_sub_1 (opx, op2_ptr, op2_size, (mp_limb_t) 1); + op2_ptr = opx; + + res_alloc = MAX (op1_size, op2_size) + 1; + if (res->_mp_alloc < res_alloc) + { + _mpz_realloc (res, res_alloc); + op1_ptr = op1->_mp_d; + res_ptr = res->_mp_d; + /* Don't re-read OP2_PTR. It points to temporary space--never + to the space RES->_mp_d used to point to before reallocation. */ + } + + if (op1_size > op2_size) + { + MPN_COPY (res_ptr + op2_size, op1_ptr + op2_size, op1_size - op2_size); + for (i = op2_size - 1; i >= 0; i--) + res_ptr[i] = op1_ptr[i] ^ op2_ptr[i]; + res_size = op1_size; + } + else + { + MPN_COPY (res_ptr + op1_size, op2_ptr + op1_size, op2_size - op1_size); + for (i = op1_size - 1; i >= 0; i--) + res_ptr[i] = op1_ptr[i] ^ op2_ptr[i]; + res_size = op2_size; + } + + cy = mpn_add_1 (res_ptr, res_ptr, res_size, (mp_limb_t) 1); + if (cy) + { + res_ptr[res_size] = cy; + res_size++; + } + + MPN_NORMALIZE (res_ptr, res_size); + res->_mp_size = -res_size; + TMP_FREE (marker); + } +} |