From: simonm Date: Fri, 5 Jun 1998 14:38:03 +0000 (+0000) Subject: [project @ 1998-06-05 14:38:01 by simonm] X-Git-Tag: Approx_2487_patches~605 X-Git-Url: http://git.megacz.com/?a=commitdiff_plain;h=3d6fbc2cd6881054ebce7991fc72e6454e0fdbe4;p=ghc-hetmet.git [project @ 1998-06-05 14:38:01 by simonm] Import GMP 2.0.2 --- diff --git a/ghc/rts/gmp/mpz/fdiv_r_2exp.c b/ghc/rts/gmp/mpz/fdiv_r_2exp.c new file mode 100644 index 0000000..04190b1 --- /dev/null +++ b/ghc/rts/gmp/mpz/fdiv_r_2exp.c @@ -0,0 +1,88 @@ +/* mpz_fdiv_r_2exp -- Divide a integer by 2**CNT and produce a remainder. + +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 +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 +License for more details. + +You should have received a copy of the GNU Library 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_fdiv_r_2exp (mpz_ptr res, mpz_srcptr in, unsigned long int cnt) +#else +mpz_fdiv_r_2exp (res, in, cnt) + mpz_ptr res; + mpz_srcptr in; + unsigned long int cnt; +#endif +{ + mp_size_t in_size = ABS (in->_mp_size); + mp_size_t res_size; + mp_size_t limb_cnt = cnt / BITS_PER_MP_LIMB; + mp_srcptr in_ptr = in->_mp_d; + + if (in_size > limb_cnt) + { + /* The input operand is (probably) greater than 2**CNT. */ + mp_limb_t x; + + x = in_ptr[limb_cnt] & (((mp_limb_t) 1 << cnt % BITS_PER_MP_LIMB) - 1); + if (x != 0) + { + res_size = limb_cnt + 1; + if (res->_mp_alloc < res_size) + _mpz_realloc (res, res_size); + + res->_mp_d[limb_cnt] = x; + } + else + { + res_size = limb_cnt; + MPN_NORMALIZE (in_ptr, res_size); + + if (res->_mp_alloc < res_size) + _mpz_realloc (res, res_size); + + limb_cnt = res_size; + } + } + 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. */ + res_size = in_size; + if (res->_mp_alloc < res_size) + _mpz_realloc (res, res_size); + + limb_cnt = res_size; + } + + if (res != in) + MPN_COPY (res->_mp_d, in->_mp_d, limb_cnt); + res->_mp_size = res_size; + if (in->_mp_size < 0 && res_size != 0) + { + /* Result should be 2^CNT - RES */ + mpz_t tmp; + MPZ_TMP_INIT (tmp, limb_cnt + 1); + mpz_set_ui (tmp, 1L); + mpz_mul_2exp (tmp, tmp, cnt); + mpz_sub (res, tmp, res); + } +} diff --git a/ghc/rts/gmp/mpz/fdiv_r_ui.c b/ghc/rts/gmp/mpz/fdiv_r_ui.c new file mode 100644 index 0000000..c4c3749 --- /dev/null +++ b/ghc/rts/gmp/mpz/fdiv_r_ui.c @@ -0,0 +1,52 @@ +/* mpz_fdiv_r_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. + +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 +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 +License for more details. + +You should have received a copy of the GNU Library 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_fdiv_r_ui (mpz_ptr rem, mpz_srcptr dividend, unsigned long int divisor) +#else +mpz_fdiv_r_ui (rem, dividend, divisor) + mpz_ptr rem; + mpz_srcptr dividend; + unsigned long int divisor; +#endif +{ + mp_size_t dividend_size; + mp_size_t size; + mp_limb_t remainder_limb; + + dividend_size = dividend->_mp_size; + size = ABS (dividend_size); + + remainder_limb = mpn_mod_1 (dividend->_mp_d, size, (mp_limb_t) divisor); + + if (remainder_limb != 0 && dividend_size < 0) + remainder_limb = divisor - remainder_limb; + + rem->_mp_d[0] = remainder_limb; + rem->_mp_size = remainder_limb != 0; + + return remainder_limb; +} diff --git a/ghc/rts/gmp/mpz/fdiv_ui.c b/ghc/rts/gmp/mpz/fdiv_ui.c new file mode 100644 index 0000000..4d018a2 --- /dev/null +++ b/ghc/rts/gmp/mpz/fdiv_ui.c @@ -0,0 +1,48 @@ +/* mpz_fdiv_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. + +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 +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 +License for more details. + +You should have received a copy of the GNU Library 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_fdiv_ui (mpz_srcptr dividend, unsigned long int divisor) +#else +mpz_fdiv_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; + + dividend_size = dividend->_mp_size; + size = ABS (dividend_size); + + remainder_limb = mpn_mod_1 (dividend->_mp_d, size, (mp_limb_t) divisor); + + if (remainder_limb != 0 && dividend_size < 0) + remainder_limb = divisor - remainder_limb; + + return remainder_limb; +} diff --git a/ghc/rts/gmp/mpz/gcd.c b/ghc/rts/gmp/mpz/gcd.c new file mode 100644 index 0000000..f93030c --- /dev/null +++ b/ghc/rts/gmp/mpz/gcd.c @@ -0,0 +1,178 @@ +/* mpz/gcd.c: Calculate the greatest common divisor of two integers. + +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 +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 +License for more details. + +You should have received a copy of the GNU Library 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 (); + +#ifndef BERKELEY_MP +void +#if __STDC__ +mpz_gcd (mpz_ptr g, mpz_srcptr u, mpz_srcptr v) +#else +mpz_gcd (g, u, v) + mpz_ptr g; + mpz_srcptr u; + mpz_srcptr v; +#endif +#else /* BERKELEY_MP */ +void +#if __STDC__ +gcd (mpz_srcptr u, mpz_srcptr v, mpz_ptr g) +#else +gcd (u, v, g) + mpz_ptr g; + mpz_srcptr u; + mpz_srcptr v; +#endif +#endif /* BERKELEY_MP */ + +{ + unsigned long int g_zero_bits, u_zero_bits, v_zero_bits; + mp_size_t g_zero_limbs, u_zero_limbs, v_zero_limbs; + mp_ptr tp; + mp_ptr up = u->_mp_d; + mp_size_t usize = ABS (u->_mp_size); + mp_ptr vp = v->_mp_d; + mp_size_t vsize = ABS (v->_mp_size); + mp_size_t gsize; + TMP_DECL (marker); + + /* GCD(0, V) == V. */ + if (usize == 0) + { + g->_mp_size = vsize; + if (g == v) + return; + if (g->_mp_alloc < vsize) + _mpz_realloc (g, vsize); + MPN_COPY (g->_mp_d, vp, vsize); + return; + } + + /* GCD(U, 0) == U. */ + if (vsize == 0) + { + g->_mp_size = usize; + if (g == u) + return; + if (g->_mp_alloc < usize) + _mpz_realloc (g, usize); + MPN_COPY (g->_mp_d, up, usize); + return; + } + + if (usize == 1) + { + g->_mp_size = 1; + g->_mp_d[0] = mpn_gcd_1 (vp, vsize, up[0]); + return; + } + + if (vsize == 1) + { + g->_mp_size = 1; + g->_mp_d[0] = mpn_gcd_1 (up, usize, vp[0]); + return; + } + + TMP_MARK (marker); + + /* Eliminate low zero bits from U and V and move to temporary storage. */ + while (*up == 0) + up++; + u_zero_limbs = up - u->_mp_d; + usize -= u_zero_limbs; + count_trailing_zeros (u_zero_bits, *up); + tp = up; + up = (mp_ptr) TMP_ALLOC (usize * BYTES_PER_MP_LIMB); + if (u_zero_bits != 0) + { + mpn_rshift (up, tp, usize, u_zero_bits); + usize -= up[usize - 1] == 0; + } + else + MPN_COPY (up, tp, usize); + + while (*vp == 0) + vp++; + v_zero_limbs = vp - v->_mp_d; + vsize -= v_zero_limbs; + count_trailing_zeros (v_zero_bits, *vp); + tp = vp; + vp = (mp_ptr) TMP_ALLOC (vsize * BYTES_PER_MP_LIMB); + if (v_zero_bits != 0) + { + mpn_rshift (vp, tp, vsize, v_zero_bits); + vsize -= vp[vsize - 1] == 0; + } + else + MPN_COPY (vp, tp, vsize); + + if (u_zero_limbs > v_zero_limbs) + { + g_zero_limbs = v_zero_limbs; + g_zero_bits = v_zero_bits; + } + else if (u_zero_limbs < v_zero_limbs) + { + g_zero_limbs = u_zero_limbs; + g_zero_bits = u_zero_bits; + } + else /* Equal. */ + { + g_zero_limbs = u_zero_limbs; + g_zero_bits = MIN (u_zero_bits, v_zero_bits); + } + + /* Call mpn_gcd. The 1st argument must not have more bits than the 2nd. */ + 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); + + /* Here G <-- V << (g_zero_limbs*BITS_PER_MP_LIMB + g_zero_bits). */ + gsize = vsize + g_zero_limbs; + if (g_zero_bits != 0) + { + mp_limb_t cy_limb; + gsize += (vp[vsize - 1] >> (BITS_PER_MP_LIMB - g_zero_bits)) != 0; + if (g->_mp_alloc < gsize) + _mpz_realloc (g, gsize); + MPN_ZERO (g->_mp_d, g_zero_limbs); + + tp = g->_mp_d + g_zero_limbs; + cy_limb = mpn_lshift (tp, vp, vsize, g_zero_bits); + if (cy_limb != 0) + tp[vsize] = cy_limb; + } + else + { + if (g->_mp_alloc < gsize) + _mpz_realloc (g, gsize); + MPN_ZERO (g->_mp_d, g_zero_limbs); + MPN_COPY (g->_mp_d + g_zero_limbs, vp, vsize); + } + + g->_mp_size = gsize; + TMP_FREE (marker); +} diff --git a/ghc/rts/gmp/mpz/gcd_ui.c b/ghc/rts/gmp/mpz/gcd_ui.c new file mode 100644 index 0000000..388ab05 --- /dev/null +++ b/ghc/rts/gmp/mpz/gcd_ui.c @@ -0,0 +1,64 @@ +/* mpz_gcd_ui -- Calculate the greatest common divisior of two integers. + +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 +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 +License for more details. + +You should have received a copy of the GNU Library 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_gcd_ui (mpz_ptr w, mpz_srcptr u, unsigned long int v) +#else +mpz_gcd_ui (w, u, v) + mpz_ptr w; + mpz_srcptr u; + unsigned long int v; +#endif +{ + mp_size_t size; + mp_limb_t res; + + size = ABS (u->_mp_size); + + if (size == 0) + res = v; + else if (v == 0) + { + if (w != NULL && u != w) + { + if (w->_mp_alloc < size) + _mpz_realloc (w, size); + + MPN_COPY (w->_mp_d, u->_mp_d, size); + } + w->_mp_size = size; + /* We can't return any useful result for gcd(big,0). */ + return size > 1 ? 0 : w->_mp_d[0]; + } + else + res = mpn_gcd_1 (u->_mp_d, size, v); + + if (w != NULL) + { + w->_mp_d[0] = res; + w->_mp_size = 1; + } + return res; +} diff --git a/ghc/rts/gmp/mpz/gcdext.c b/ghc/rts/gmp/mpz/gcdext.c new file mode 100644 index 0000000..adf66b0 --- /dev/null +++ b/ghc/rts/gmp/mpz/gcdext.c @@ -0,0 +1,88 @@ +/* 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. + +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 +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 +License for more details. + +You should have received a copy of the GNU Library 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" + +/* Botch: SLOW! */ + +void +#if __STDC__ +mpz_gcdext (mpz_ptr g, mpz_ptr s, mpz_ptr t, mpz_srcptr a, mpz_srcptr b) +#else +mpz_gcdext (g, s, t, a, b) + mpz_ptr g; + mpz_ptr s; + mpz_ptr t; + mpz_srcptr a; + mpz_srcptr b; +#endif +{ + mpz_t s0, s1, q, r, x, d0, d1; + + 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); + + while (d1->_mp_size != 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); + } + + if (t != NULL) + { + mpz_mul (x, s0, a); + mpz_sub (x, d0, x); + if (b->_mp_size == 0) + t->_mp_size = 0; + else + mpz_tdiv_q (t, x, b); + } + mpz_set (s, s0); + mpz_set (g, d0); + if (g->_mp_size < 0) + { + g->_mp_size = -g->_mp_size; + s->_mp_size = -s->_mp_size; + if (t != NULL) + t->_mp_size = -t->_mp_size; + } + + mpz_clear (s0); + mpz_clear (s1); + mpz_clear (q); + mpz_clear (r); + mpz_clear (x); + mpz_clear (d0); + mpz_clear (d1); +} diff --git a/ghc/rts/gmp/mpz/get_d.c b/ghc/rts/gmp/mpz/get_d.c new file mode 100644 index 0000000..0fd7916 --- /dev/null +++ b/ghc/rts/gmp/mpz/get_d.c @@ -0,0 +1,54 @@ +/* double mpz_get_d (mpz_t src) -- Return the double approximation to SRC. + +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 +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 +License for more details. + +You should have received a copy of the GNU Library 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" + +double +#if __STDC__ +mpz_get_d (mpz_srcptr src) +#else +mpz_get_d (src) + mpz_srcptr src; +#endif +{ + double res; + mp_size_t size, i, n_limbs_to_use; + int negative; + mp_ptr qp; + + size = SIZ(src); + if (size == 0) + return 0.0; + + negative = size < 0; + 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]; + + res = __gmp_scale2 (res, (size - n_limbs_to_use) * BITS_PER_MP_LIMB); + + return negative ? -res : res; +} diff --git a/ghc/rts/gmp/mpz/get_si.c b/ghc/rts/gmp/mpz/get_si.c new file mode 100644 index 0000000..45e0e5a --- /dev/null +++ b/ghc/rts/gmp/mpz/get_si.c @@ -0,0 +1,43 @@ +/* mpz_get_si(integer) -- Return the least significant digit from INTEGER. + +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 +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 +License for more details. + +You should have received a copy of the GNU Library 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" + +signed long int +#if __STDC__ +mpz_get_si (mpz_srcptr op) +#else +mpz_get_si (op) + mpz_srcptr op; +#endif +{ + mp_size_t size = op->_mp_size; + mp_limb_t low_limb = op->_mp_d[0]; + + if (size > 0) + return low_limb % ((mp_limb_t) 1 << (BITS_PER_MP_LIMB - 1)); + else if (size < 0) + /* This convoluted expression is necessary to properly handle 0x80000000 */ + return ~((low_limb - 1) % ((mp_limb_t) 1 << (BITS_PER_MP_LIMB - 1))); + else + return 0; +} diff --git a/ghc/rts/gmp/mpz/get_str.c b/ghc/rts/gmp/mpz/get_str.c new file mode 100644 index 0000000..8ccf3ef --- /dev/null +++ b/ghc/rts/gmp/mpz/get_str.c @@ -0,0 +1,118 @@ +/* mpz_get_str (string, base, mp_src) -- Convert the multiple precision + number MP_SRC to a string STRING of base BASE. If STRING is NULL + allocate space for the result. In any case, return a pointer to the + result. If STRING is not NULL, the caller must ensure enough space is + available to store the result. + +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 +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 +License for more details. + +You should have received a copy of the GNU Library 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" + +char * +#if __STDC__ +mpz_get_str (char *res_str, int base, mpz_srcptr x) +#else +mpz_get_str (res_str, base, x) + char *res_str; + int base; + mpz_srcptr x; +#endif +{ + mp_ptr xp; + mp_size_t x_size = x->_mp_size; + unsigned char *str; + char *return_str; + size_t str_size; + char *num_to_text; + int i; + TMP_DECL (marker); + + TMP_MARK (marker); + if (base >= 0) + { + if (base == 0) + base = 10; + num_to_text = "0123456789abcdefghijklmnopqrstuvwxyz"; + } + else + { + base = -base; + num_to_text = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ"; + } + + /* We allways allocate space for the string. If the caller passed a + NULL pointer for RES_STR, we allocate permanent space and return + a pointer to that to the caller. */ + str_size = ((size_t) (ABS (x_size) * BITS_PER_MP_LIMB + * __mp_bases[base].chars_per_bit_exactly)) + 3; + if (res_str == 0) + { + /* We didn't get a string from the user. Allocate one (and return + a pointer to it). */ + res_str = (char *) (*_mp_allocate_func) (str_size); + /* Make str, the variable used for raw result from mpn_get_str, + point to the same string, but just after a possible minus sign. */ + str = (unsigned char *) res_str + 1; + } + else + { + /* Use TMP_ALLOC to get temporary space, since we need a few extra bytes + that we can't expect to caller to supply us with. */ + str = (unsigned char *) TMP_ALLOC (str_size); + } + + return_str = res_str; + + if (x_size == 0) + { + res_str[0] = '0'; + res_str[1] = 0; + TMP_FREE (marker); + return res_str; + } + if (x_size < 0) + { + *res_str++ = '-'; + x_size = -x_size; + } + + /* Move the number to convert into temporary space, since mpn_get_str + clobbers its argument + needs one extra high limb.... */ + xp = (mp_ptr) TMP_ALLOC ((x_size + 1) * BYTES_PER_MP_LIMB); + MPN_COPY (xp, x->_mp_d, x_size); + + str_size = mpn_get_str (str, base, xp, x_size); + + /* mpn_get_str might make some leading zeros. Skip them. */ + while (*str == 0) + { + str_size--; + str++; + } + + /* Translate result to printable chars and move result to RES_STR. */ + for (i = 0; i < str_size; i++) + res_str[i] = num_to_text[str[i]]; + res_str[str_size] = 0; + + TMP_FREE (marker); + return return_str; +} diff --git a/ghc/rts/gmp/mpz/get_ui.c b/ghc/rts/gmp/mpz/get_ui.c new file mode 100644 index 0000000..4bfb5e1 --- /dev/null +++ b/ghc/rts/gmp/mpz/get_ui.c @@ -0,0 +1,37 @@ +/* mpz_get_ui(integer) -- Return the least significant digit from INTEGER. + +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 +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 +License for more details. + +You should have received a copy of the GNU Library 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_get_ui (mpz_srcptr integer) +#else +mpz_get_ui (integer) + mpz_srcptr integer; +#endif +{ + if (integer->_mp_size == 0) + return 0; + else + return integer->_mp_d[0]; +} diff --git a/ghc/rts/gmp/mpz/getlimbn.c b/ghc/rts/gmp/mpz/getlimbn.c new file mode 100644 index 0000000..c7a234b --- /dev/null +++ b/ghc/rts/gmp/mpz/getlimbn.c @@ -0,0 +1,38 @@ +/* mpz_getlimbn(integer,n) -- Return the N:th limb from INTEGER. + +Copyright (C) 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 +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 +License for more details. + +You should have received a copy of the GNU Library 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" + +mp_limb_t +#if __STDC__ +mpz_getlimbn (mpz_srcptr integer, mp_size_t n) +#else +mpz_getlimbn (integer, n) + mpz_srcptr integer; + mp_size_t n; +#endif +{ + if (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 new file mode 100644 index 0000000..58c9273 --- /dev/null +++ b/ghc/rts/gmp/mpz/hamdist.c @@ -0,0 +1,62 @@ +/* mpz_hamdist(mpz_ptr op1, mpz_ptr op2) -- Compute the hamming distance + between OP1 and OP2. If one of the operands is negative, return ~0. (We + could make the function well-defined when both operands are negative, but + that would probably not be worth the trouble. + +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 +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 +License for more details. + +You should have received a copy of the GNU Library 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_hamdist (mpz_srcptr u, mpz_srcptr v) +#else +mpz_hamdist (u, v) + mpz_srcptr u; + mpz_srcptr v; +#endif +{ + mp_srcptr up, vp; + mp_size_t usize, vsize, size; + unsigned long int count; + + usize = u->_mp_size; + vsize = v->_mp_size; + + if ((usize | vsize) < 0) + return ~ (unsigned long int) 0; + + up = u->_mp_d; + vp = v->_mp_d; + + if (usize > vsize) + { + count = mpn_popcount (up + vsize, usize - vsize); + size = vsize; + } + else + { + count = mpn_popcount (vp + usize, vsize - usize); + size = usize; + } + + return count + mpn_hamdist (up, vp, size); +} diff --git a/ghc/rts/gmp/mpz/init.c b/ghc/rts/gmp/mpz/init.c new file mode 100644 index 0000000..f8d8e20 --- /dev/null +++ b/ghc/rts/gmp/mpz/init.c @@ -0,0 +1,36 @@ +/* mpz_init() -- Make a new multiple precision number with value 0. + +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 +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 +License for more details. + +You should have received a copy of the GNU Library 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_init (mpz_ptr x) +#else +mpz_init (x) + mpz_ptr x; +#endif +{ + x->_mp_alloc = 1; + x->_mp_d = (mp_ptr) (*_mp_allocate_func) (BYTES_PER_MP_LIMB); + x->_mp_size = 0; +} diff --git a/ghc/rts/gmp/mpz/inp_raw.c b/ghc/rts/gmp/mpz/inp_raw.c new file mode 100644 index 0000000..e1cec1d --- /dev/null +++ b/ghc/rts/gmp/mpz/inp_raw.c @@ -0,0 +1,101 @@ +/* mpz_inp_raw -- Input a mpz_t in raw, but endianess, and wordsize + independent format (as output by mpz_out_raw). + +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 +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 +License for more details. + +You should have received a copy of the GNU Library 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 + +#include "gmp.h" +#include "gmp-impl.h" + +size_t +#if __STDC__ +mpz_inp_raw (mpz_ptr x, FILE *stream) +#else +mpz_inp_raw (x, stream) + mpz_ptr x; + FILE *stream; +#endif +{ + int i; + mp_size_t s; + mp_size_t xsize; + mp_ptr xp; + unsigned int c; + mp_limb_t x_limb; + mp_size_t in_bytesize; + int neg_flag; + + if (stream == 0) + stream = stdin; + + /* Read 4-byte size */ + in_bytesize = 0; + for (i = 4 - 1; i >= 0; i--) + { + c = fgetc (stream); + in_bytesize = (in_bytesize << BITS_PER_CHAR) | c; + } + + /* Size is stored as a 32 bit word; sign extend in_bytesize for non-32 bit + machines. */ + if (sizeof (mp_size_t) > 4) + in_bytesize |= (-(in_bytesize < 0)) << 31; + + neg_flag = in_bytesize < 0; + in_bytesize = ABS (in_bytesize); + xsize = (in_bytesize + BYTES_PER_MP_LIMB - 1) / BYTES_PER_MP_LIMB; + + if (xsize == 0) + { + x->_mp_size = 0; + return 4; /* we've read 4 bytes */ + } + + if (x->_mp_alloc < xsize) + _mpz_realloc (x, xsize); + xp = x->_mp_d; + + x_limb = 0; + for (i = (in_bytesize - 1) % BYTES_PER_MP_LIMB; i >= 0; i--) + { + c = fgetc (stream); + x_limb = (x_limb << BITS_PER_CHAR) | c; + } + xp[xsize - 1] = x_limb; + + for (s = xsize - 2; s >= 0; s--) + { + x_limb = 0; + for (i = BYTES_PER_MP_LIMB - 1; i >= 0; i--) + { + c = fgetc (stream); + x_limb = (x_limb << BITS_PER_CHAR) | c; + } + xp[s] = x_limb; + } + + if (c == EOF) + return 0; /* error */ + + MPN_NORMALIZE (xp, xsize); + x->_mp_size = neg_flag ? -xsize : xsize; + return in_bytesize + 4; +} diff --git a/ghc/rts/gmp/mpz/inp_str.c b/ghc/rts/gmp/mpz/inp_str.c new file mode 100644 index 0000000..7159062 --- /dev/null +++ b/ghc/rts/gmp/mpz/inp_str.c @@ -0,0 +1,138 @@ +/* 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. + +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 +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 +License for more details. + +You should have received a copy of the GNU Library 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 +#include +#include "gmp.h" +#include "gmp-impl.h" + +static int +digit_value_in_base (c, base) + int c; + int base; +{ + int digit; + + if (isdigit (c)) + digit = c - '0'; + else if (islower (c)) + digit = c - 'a' + 10; + else if (isupper (c)) + digit = c - 'A' + 10; + else + return -1; + + if (digit < base) + return digit; + return -1; +} + +size_t +#if __STDC__ +mpz_inp_str (mpz_ptr dest, FILE *stream, int base) +#else +mpz_inp_str (dest, stream, base) + mpz_ptr dest; + FILE *stream; + int base; +#endif +{ + char *str; + size_t alloc_size, str_size; + int c; + int negative; + mp_size_t dest_size; + 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. */ + do + { + c = getc (stream); + nread++; + } + while (isspace (c)); + + negative = 0; + if (c == '-') + { + negative = 1; + c = getc (stream); + } + + if (digit_value_in_base (c, base == 0 ? 10 : base) < 0) + return 0; /* error if no digits */ + + /* If BASE is 0, try to find out the base by looking at the initial + characters. */ + if (base == 0) + { + base = 10; + if (c == '0') + { + base = 8; + c = getc (stream); + nread++; + if (c == 'x' || c == 'X') + { + base = 16; + c = getc (stream); + nread++; + } + } + } + + for (;;) + { + int dig; + if (str_size >= alloc_size) + { + size_t old_alloc_size = alloc_size; + alloc_size = alloc_size * 3 / 2; + str = (char *) (*_mp_reallocate_func) (str, old_alloc_size, alloc_size); + } + dig = digit_value_in_base (c, base); + if (dig < 0) + break; + str[str_size++] = dig; + c = getc (stream); + } + + 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); + + dest_size = mpn_set_str (dest->_mp_d, (unsigned char *) str, str_size, base); + dest->_mp_size = negative ? -dest_size : dest_size; + + (*_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 new file mode 100644 index 0000000..ff1d6d9 --- /dev/null +++ b/ghc/rts/gmp/mpz/invert.c @@ -0,0 +1,43 @@ +/* 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. + +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 +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 +License for more details. + +You should have received a copy of the GNU Library 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" + +int +#if __STDC__ +mpz_invert (mpz_ptr inverse, mpz_srcptr x, mpz_srcptr n) +#else +mpz_invert (inverse, x, n) + mpz_ptr inverse; + 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; +} diff --git a/ghc/rts/gmp/mpz/ior.c b/ghc/rts/gmp/mpz/ior.c new file mode 100644 index 0000000..77facfd --- /dev/null +++ b/ghc/rts/gmp/mpz/ior.c @@ -0,0 +1,243 @@ +/* mpz_ior -- Logical inclusive or. + +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 +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 +License for more details. + +You should have received a copy of the GNU Library 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_ior (mpz_ptr res, mpz_srcptr op1, mpz_srcptr op2) +#else +mpz_ior (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; + 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; + } + + 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, so will be the result. + -((-OP1) | (-OP2)) = -(~(OP1 - 1) | ~(OP2 - 1)) = + = ~(~(OP1 - 1) | ~(OP2 - 1)) + 1 = + = ((OP1 - 1) & (OP2 - 1)) + 1 */ + + op1_size = -op1_size; + op2_size = -op2_size; + + res_size = MIN (op1_size, op2_size); + + /* Possible optimization: Decrease mpn_sub precision, + as we won't use the entire res of both. */ + opx = (mp_ptr) TMP_ALLOC (res_size * BYTES_PER_MP_LIMB); + mpn_sub_1 (opx, op1_ptr, res_size, (mp_limb_t) 1); + op1_ptr = opx; + + opx = (mp_ptr) TMP_ALLOC (res_size * BYTES_PER_MP_LIMB); + mpn_sub_1 (opx, op2_ptr, res_size, (mp_limb_t) 1); + op2_ptr = opx; + + if (res->_mp_alloc < res_size) + { + _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 + to point to before reallocation. */ + } + + /* First loop finds the size of the result. */ + for (i = res_size - 1; i >= 0; i--) + if ((op1_ptr[i] & op2_ptr[i]) != 0) + break; + res_size = i + 1; + + if (res_size != 0) + { + /* Second loop computes the real result. */ + for (i = res_size - 1; i >= 0; i--) + res_ptr[i] = op1_ptr[i] & op2_ptr[i]; + + cy = mpn_add_1 (res_ptr, res_ptr, res_size, (mp_limb_t) 1); + if (cy) + { + res_ptr[res_size] = cy; + res_size++; + } + } + else + { + res_ptr[0] = 1; + res_size = 1; + } + + 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 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;} + } + } + + { + mp_ptr opx; + mp_limb_t cy; + mp_size_t res_alloc; + 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; + + res_alloc = 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; + + 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) + { + /* We can just ignore the part of OP1 that stretches above OP2, + because the result limbs are zero there. */ + + /* First loop finds the size of the result. */ + for (i = op2_size - 1; i >= 0; i--) + if ((~op1_ptr[i] & op2_ptr[i]) != 0) + break; + res_size = i + 1; + count = res_size; + } + else + { + res_size = op2_size; + + /* Copy the part of OP2 that stretches above OP1, to RES. */ + MPN_COPY (res_ptr + op1_size, op2_ptr + op1_size, op2_size - op1_size); + count = op1_size; + } + + if (res_size != 0) + { + /* Second loop computes the real result. */ + for (i = count - 1; i >= 0; i--) + res_ptr[i] = ~op1_ptr[i] & op2_ptr[i]; + + cy = mpn_add_1 (res_ptr, res_ptr, res_size, (mp_limb_t) 1); + if (cy) + { + res_ptr[res_size] = cy; + res_size++; + } + } + else + { + res_ptr[0] = 1; + res_size = 1; + } + + res->_mp_size = -res_size; + } + TMP_FREE (marker); +} diff --git a/ghc/rts/gmp/mpz/iset.c b/ghc/rts/gmp/mpz/iset.c new file mode 100644 index 0000000..c8a17dc --- /dev/null +++ b/ghc/rts/gmp/mpz/iset.c @@ -0,0 +1,49 @@ +/* mpz_init_set (src_integer) -- Make a new multiple precision number with + a value copied from SRC_INTEGER. + +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 +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 +License for more details. + +You should have received a copy of the GNU Library 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_init_set (mpz_ptr w, mpz_srcptr u) +#else +mpz_init_set (w, u) + mpz_ptr w; + mpz_srcptr u; +#endif +{ + mp_ptr wp, up; + mp_size_t usize, size; + + usize = u->_mp_size; + size = ABS (usize); + + w->_mp_alloc = MAX (size, 1); + w->_mp_d = (mp_ptr) (*_mp_allocate_func) (w->_mp_alloc * BYTES_PER_MP_LIMB); + + wp = w->_mp_d; + up = u->_mp_d; + + MPN_COPY (wp, up, size); + w->_mp_size = usize; +} diff --git a/ghc/rts/gmp/mpz/iset_d.c b/ghc/rts/gmp/mpz/iset_d.c new file mode 100644 index 0000000..41e5c4f --- /dev/null +++ b/ghc/rts/gmp/mpz/iset_d.c @@ -0,0 +1,39 @@ +/* mpz_init_set_d(integer, val) -- Initialize and assign INTEGER with a double + value VAL. + +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 +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 +License for more details. + +You should have received a copy of the GNU Library 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_init_set_d (mpz_ptr dest, double val) +#else +mpz_init_set_d (dest, val) + mpz_ptr dest; + double val; +#endif +{ + dest->_mp_alloc = 1; + dest->_mp_d = (mp_ptr) (*_mp_allocate_func) (BYTES_PER_MP_LIMB); + dest->_mp_size = 0; + mpz_set_d (dest, val); +} diff --git a/ghc/rts/gmp/mpz/iset_si.c b/ghc/rts/gmp/mpz/iset_si.c new file mode 100644 index 0000000..af51f05 --- /dev/null +++ b/ghc/rts/gmp/mpz/iset_si.c @@ -0,0 +1,49 @@ +/* mpz_init_set_si(val) -- Make a new multiple precision number with + value val. + +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 +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 +License for more details. + +You should have received a copy of the GNU Library 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_init_set_si (mpz_ptr x, signed long int val) +#else +mpz_init_set_si (x, val) + mpz_ptr x; + signed long int val; +#endif +{ + x->_mp_alloc = 1; + x->_mp_d = (mp_ptr) (*_mp_allocate_func) (BYTES_PER_MP_LIMB); + if (val > 0) + { + x->_mp_d[0] = val; + x->_mp_size = 1; + } + else if (val < 0) + { + x->_mp_d[0] = -val; + x->_mp_size = -1; + } + else + x->_mp_size = 0; +} diff --git a/ghc/rts/gmp/mpz/iset_str.c b/ghc/rts/gmp/mpz/iset_str.c new file mode 100644 index 0000000..e04ad5d --- /dev/null +++ b/ghc/rts/gmp/mpz/iset_str.c @@ -0,0 +1,44 @@ +/* mpz_init_set_str(string, base) -- Convert the \0-terminated string + STRING in base BASE to a multiple precision integer. Return a MP_INT + structure representing the integer. Allow white space in the + string. If BASE == 0 determine 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, 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 +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 +License for more details. + +You should have received a copy of the GNU Library 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_init_set_str (mpz_ptr x, const char *str, int base) +#else +mpz_init_set_str (x, str, base) + mpz_ptr x; + const char *str; + int base; +#endif +{ + x->_mp_alloc = 1; + x->_mp_d = (mp_ptr) (*_mp_allocate_func) (BYTES_PER_MP_LIMB); + + 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 new file mode 100644 index 0000000..dc39f59 --- /dev/null +++ b/ghc/rts/gmp/mpz/iset_ui.c @@ -0,0 +1,39 @@ +/* mpz_init_set_ui(val) -- Make a new multiple precision number with + value val. + +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 +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 +License for more details. + +You should have received a copy of the GNU Library 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_init_set_ui (mpz_ptr x, unsigned long int val) +#else +mpz_init_set_ui (x, val) + mpz_ptr x; + unsigned long int val; +#endif +{ + x->_mp_alloc = 1; + x->_mp_d = (mp_ptr) (*_mp_allocate_func) (BYTES_PER_MP_LIMB); + x->_mp_d[0] = val; + x->_mp_size = val != 0; +} diff --git a/ghc/rts/gmp/mpz/jacobi.c b/ghc/rts/gmp/mpz/jacobi.c new file mode 100644 index 0000000..409f622 --- /dev/null +++ b/ghc/rts/gmp/mpz/jacobi.c @@ -0,0 +1,53 @@ +/* mpz_jacobi (op1, op2). + Contributed by Bennet Yee (bsy) at Carnegie-Mellon University + +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 +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 +License for more details. + +You should have received a copy of the GNU Library 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" + +/* Precondition: both p and q are positive */ + +int +#if __STDC__ +mpz_jacobi (mpz_srcptr pi, mpz_srcptr qi) +#else +mpz_jacobi (pi, qi) + mpz_srcptr pi, qi; +#endif +{ +#if GCDCHECK + int retval; + mpz_t gcdval; + + mpz_init (gcdval); + mpz_gcd (gcdval, pi, qi); + if (!mpz_cmp_ui (gcdval, 1L)) + { + /* J(ab,cb) = J(ab,c)J(ab,b) = J(ab,c)J(0,b) = J(ab,c)*0 */ + retval = 0; + } + else + retval = mpz_legendre (pi, qi); + mpz_clear (gcdval); + return retval; +#else + return mpz_legendre (pi, qi); +#endif +} diff --git a/ghc/rts/gmp/mpz/legendre.c b/ghc/rts/gmp/mpz/legendre.c new file mode 100644 index 0000000..4de16a6 --- /dev/null +++ b/ghc/rts/gmp/mpz/legendre.c @@ -0,0 +1,184 @@ +/* mpz_legendre (op1, op2). + Contributed by Bennet Yee (bsy) at Carnegie-Mellon University + +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 +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 +License for more details. + +You should have received a copy of the GNU Library 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" + +#if defined (DEBUG) +#include +#endif + +/* Precondition: both p and q are positive */ + +int +#if __STDC__ +mpz_legendre (mpz_srcptr pi, mpz_srcptr qi) +#else +mpz_legendre (pi, qi) +mpz_srcptr pi, qi; +#endif +{ + mpz_t p, q, qdiv2; +#ifdef Q_MINUS_1 + mpz_t q_minus_1; +#endif + mpz_ptr mtmp; + register mpz_ptr pptr, qptr; + register int retval = 1; + register unsigned long int s; + + pptr = p; + mpz_init_set (pptr, pi); + qptr = q; + mpz_init_set (qptr, qi); + +#ifdef Q_MINUS_1 + mpz_init (q_minus_1); +#endif + mpz_init (qdiv2); + +tail_recurse2: +#ifdef DEBUG + printf ("tail_recurse2: p="); + mpz_out_str (stdout, 10, pptr); + printf ("\nq="); + mpz_out_str (stdout, 10, qptr); + putchar ('\n'); +#endif + s = mpz_scan1 (qptr, 0); + if (s) mpz_tdiv_q_2exp (qptr, qptr, s); /* J(a,2) = 1 */ +#ifdef DEBUG + printf ("2 factor decomposition: p="); + mpz_out_str (stdout, 10, pptr); + printf ("\nq="); + mpz_out_str (stdout, 10, qptr); + putchar ('\n'); +#endif + /* postcondition q odd */ + if (!mpz_cmp_ui (qptr, 1L)) /* J(a,1) = 1 */ + goto done; + mpz_mod (pptr, pptr, qptr); /* J(a,q) = J(b,q) when a == b mod q */ +#ifdef DEBUG + printf ("mod out by q: p="); + mpz_out_str (stdout, 10, pptr); + printf ("\nq="); + mpz_out_str (stdout, 10, qptr); + putchar ('\n'); +#endif + /* quick calculation to get approximate size first */ + /* precondition: p < q */ + if ((mpz_sizeinbase (pptr, 2) + 1 >= mpz_sizeinbase (qptr,2)) + && (mpz_tdiv_q_2exp (qdiv2, qptr, 1L), mpz_cmp (pptr, qdiv2) > 0)) + { + /* p > q/2 */ + mpz_sub (pptr, qptr, pptr); + /* J(-1,q) = (-1)^((q-1)/2), q odd */ + if (mpz_get_ui (qptr) & 2) + retval = -retval; + } + /* p < q/2 */ +#ifdef Q_MINUS_1 + mpz_sub_ui (q_minus_q, qptr, 1L); +#endif +tail_recurse: /* we use tail_recurse only if q has not changed */ +#ifdef DEBUG + printf ("tail_recurse1: p="); + mpz_out_str (stdout, 10, pptr); + printf ("\nq="); + mpz_out_str (stdout, 10, qptr); + putchar ('\n'); +#endif + /* + * J(0,q) = 0 + * this occurs only if gcd(p,q) != 1 which is never true for + * Legendre function. + */ + if (!mpz_cmp_ui (pptr, 0L)) + { + retval = 0; + goto done; + } + + if (!mpz_cmp_ui (pptr, 1L)) + { + /* J(1,q) = 1 */ + /* retval *= 1; */ + goto done; + } +#ifdef Q_MINUS_1 + if (!mpz_cmp (pptr, q_minus_1)) + { + /* J(-1,q) = (-1)^((q-1)/2) */ + if (mpz_get_ui (qptr) & 2) + retval = -retval; + /* else retval *= 1; */ + goto done; + } +#endif + /* + * we do not handle J(xy,q) except for x==2 + * since we do not want to factor + */ + if ((s = mpz_scan1 (pptr, 0)) != 0) + { + /* + * J(2,q) = (-1)^((q^2-1)/8) + * + * Note that q odd guarantees that q^2-1 is divisible by 8: + * Let a: q=2a+1. q^2 = 4a^2+4a+1, (q^2-1)/8 = a(a+1)/2, qed + * + * Now, note that this means that the low two bits of _a_ + * (or the low bits of q shifted over by 1 determines + * the factor). + */ + mpz_tdiv_q_2exp (pptr, pptr, s); + + /* even powers of 2 gives J(2,q)^{2n} = 1 */ + if (s & 1) + { + s = mpz_get_ui (qptr) >> 1; + s = s * (s + 1); + if (s & 2) + retval = -retval; + } + goto tail_recurse; + } + /* + * we know p is odd since we have cast out 2s + * precondition that q is odd guarantees both odd. + * + * quadratic reciprocity + * J(p,q) = (-1)^((p-1)(q-1)/4) * J(q,p) + */ + if ((s = mpz_scan1 (pptr, 1)) <= 2 && (s + mpz_scan1 (qptr, 1)) <= 2) + retval = -retval; + + mtmp = pptr; pptr = qptr; qptr = mtmp; + goto tail_recurse2; +done: + mpz_clear (p); + mpz_clear (q); + mpz_clear (qdiv2); +#ifdef Q_MINUS_1 + mpz_clear (q_minus_1); +#endif + return retval; +} diff --git a/ghc/rts/gmp/mpz/mod.c b/ghc/rts/gmp/mpz/mod.c new file mode 100644 index 0000000..b2b8b39 --- /dev/null +++ b/ghc/rts/gmp/mpz/mod.c @@ -0,0 +1,63 @@ +/* mpz_mod -- The mathematical mod function. + +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 +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 +License for more details. + +You should have received a copy of the GNU Library 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_mod (mpz_ptr rem, mpz_srcptr dividend, mpz_srcptr divisor) +#else +mpz_mod (rem, dividend, divisor) + mpz_ptr rem; + mpz_srcptr dividend; + mpz_srcptr divisor; +#endif +{ + mp_size_t divisor_size = divisor->_mp_size; + mpz_t temp_divisor; /* N.B.: lives until function returns! */ + TMP_DECL (marker); + + TMP_MARK (marker); + + /* We need the original value of the divisor after the remainder has been + preliminary calculated. We have to copy it to temporary space if it's + the same variable as REM. */ + if (rem == divisor) + { + MPZ_TMP_INIT (temp_divisor, ABS (divisor_size)); + mpz_set (temp_divisor, divisor); + divisor = temp_divisor; + } + + mpz_tdiv_r (rem, dividend, divisor); + + if (rem->_mp_size != 0) + { + if (dividend->_mp_size < 0) + if (divisor->_mp_size < 0) + mpz_sub (rem, rem, divisor); + else + mpz_add (rem, rem, divisor); + } + + TMP_FREE (marker); +} diff --git a/ghc/rts/gmp/mpz/mul.c b/ghc/rts/gmp/mpz/mul.c new file mode 100644 index 0000000..47ce8e3 --- /dev/null +++ b/ghc/rts/gmp/mpz/mul.c @@ -0,0 +1,127 @@ +/* mpz_mul -- Multiply two integers. + +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 +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 +License for more details. + +You should have received a copy of the GNU Library 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" + +#ifndef BERKELEY_MP +void +#if __STDC__ +mpz_mul (mpz_ptr w, mpz_srcptr u, mpz_srcptr v) +#else +mpz_mul (w, u, v) + mpz_ptr w; + mpz_srcptr u; + mpz_srcptr v; +#endif +#else /* BERKELEY_MP */ +void +#if __STDC__ +mult (mpz_srcptr u, mpz_srcptr v, mpz_ptr w) +#else +mult (u, v, w) + mpz_srcptr u; + mpz_srcptr v; + mpz_ptr w; +#endif +#endif /* BERKELEY_MP */ +{ + mp_size_t usize = u->_mp_size; + mp_size_t vsize = v->_mp_size; + mp_size_t wsize; + mp_size_t sign_product; + mp_ptr up, vp; + mp_ptr wp; + mp_ptr free_me = NULL; + size_t free_me_size; + mp_limb_t cy_limb; + TMP_DECL (marker); + + TMP_MARK (marker); + sign_product = usize ^ vsize; + usize = ABS (usize); + vsize = ABS (vsize); + + if (usize < vsize) + { + /* Swap U and V. */ + {const __mpz_struct *t = u; u = v; v = t;} + {mp_size_t t = usize; usize = vsize; vsize = t;} + } + + up = u->_mp_d; + vp = v->_mp_d; + wp = w->_mp_d; + + /* Ensure W has space enough to store the result. */ + wsize = usize + vsize; + if (w->_mp_alloc < wsize) + { + if (wp == up || wp == vp) + { + free_me = wp; + free_me_size = w->_mp_alloc; + } + else + (*_mp_free_func) (wp, w->_mp_alloc * BYTES_PER_MP_LIMB); + + w->_mp_alloc = wsize; + wp = (mp_ptr) (*_mp_allocate_func) (wsize * BYTES_PER_MP_LIMB); + w->_mp_d = wp; + } + else + { + /* Make U and V not overlap with W. */ + if (wp == up) + { + /* W and U are identical. Allocate temporary space for U. */ + up = (mp_ptr) TMP_ALLOC (usize * BYTES_PER_MP_LIMB); + /* Is V identical too? Keep it identical with U. */ + if (wp == vp) + vp = up; + /* Copy to the temporary space. */ + MPN_COPY (up, wp, usize); + } + else if (wp == vp) + { + /* W and V are identical. Allocate temporary space for V. */ + vp = (mp_ptr) TMP_ALLOC (vsize * BYTES_PER_MP_LIMB); + /* Copy to the temporary space. */ + MPN_COPY (vp, wp, vsize); + } + } + + if (vsize == 0) + { + wsize = 0; + } + else + { + cy_limb = mpn_mul (wp, up, usize, vp, vsize); + wsize = usize + vsize; + wsize -= cy_limb == 0; + } + + w->_mp_size = sign_product < 0 ? -wsize : wsize; + if (free_me != NULL) + (*_mp_free_func) (free_me, free_me_size * BYTES_PER_MP_LIMB); + TMP_FREE (marker); +} diff --git a/ghc/rts/gmp/mpz/mul_2exp.c b/ghc/rts/gmp/mpz/mul_2exp.c new file mode 100644 index 0000000..4d66a98 --- /dev/null +++ b/ghc/rts/gmp/mpz/mul_2exp.c @@ -0,0 +1,76 @@ +/* mpz_mul_2exp -- Multiply a bignum by 2**CNT + +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 +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 +License for more details. + +You should have received a copy of the GNU Library 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_mul_2exp (mpz_ptr w, mpz_srcptr u, unsigned long int cnt) +#else +mpz_mul_2exp (w, u, cnt) + mpz_ptr w; + mpz_srcptr u; + unsigned long int cnt; +#endif +{ + mp_size_t usize = u->_mp_size; + mp_size_t abs_usize = ABS (usize); + mp_size_t wsize; + mp_size_t limb_cnt; + mp_ptr wp; + mp_limb_t wlimb; + + if (usize == 0) + { + w->_mp_size = 0; + return; + } + + limb_cnt = cnt / BITS_PER_MP_LIMB; + wsize = abs_usize + limb_cnt + 1; + if (w->_mp_alloc < wsize) + _mpz_realloc (w, wsize); + + wp = w->_mp_d; + wsize = abs_usize + limb_cnt; + + cnt %= BITS_PER_MP_LIMB; + if (cnt != 0) + { + wlimb = mpn_lshift (wp + limb_cnt, u->_mp_d, abs_usize, cnt); + if (wlimb != 0) + { + wp[wsize] = wlimb; + wsize++; + } + } + else + { + MPN_COPY_DECR (wp + limb_cnt, u->_mp_d, abs_usize); + } + + /* Zero all whole limbs at low end. Do it here and not before calling + mpn_lshift, not to lose for U == W. */ + MPN_ZERO (wp, limb_cnt); + + w->_mp_size = usize >= 0 ? wsize : -wsize; +} diff --git a/ghc/rts/gmp/mpz/neg.c b/ghc/rts/gmp/mpz/neg.c new file mode 100644 index 0000000..0b48e5c --- /dev/null +++ b/ghc/rts/gmp/mpz/neg.c @@ -0,0 +1,53 @@ +/* mpz_neg(mpz_ptr dst, mpz_ptr src) -- Assign the negated value of SRC to DST. + +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 +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 +License for more details. + +You should have received a copy of the GNU Library 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_neg (mpz_ptr w, mpz_srcptr u) +#else +mpz_neg (w, u) + mpz_ptr w; + mpz_srcptr u; +#endif +{ + mp_ptr wp, up; + mp_size_t usize, size; + + usize = u->_mp_size; + + if (u != w) + { + size = ABS (usize); + + if (w->_mp_alloc < size) + _mpz_realloc (w, size); + + wp = w->_mp_d; + up = u->_mp_d; + + MPN_COPY (wp, up, size); + } + + w->_mp_size = -usize; +} diff --git a/ghc/rts/gmp/mpz/out_raw.c b/ghc/rts/gmp/mpz/out_raw.c new file mode 100644 index 0000000..35d311b --- /dev/null +++ b/ghc/rts/gmp/mpz/out_raw.c @@ -0,0 +1,89 @@ +/* mpz_out_raw -- Output a mpz_t in binary. Use an endianess and word size + independent format. + +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 +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 +License for more details. + +You should have received a copy of the GNU Library 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 + +#include "gmp.h" +#include "gmp-impl.h" + +size_t +#if __STDC__ +mpz_out_raw (FILE *stream, mpz_srcptr x) +#else +mpz_out_raw (stream, x) + FILE *stream; + mpz_srcptr x; +#endif +{ + int i; + mp_size_t s; + mp_size_t xsize = ABS (x->_mp_size); + mp_srcptr xp = x->_mp_d; + mp_size_t out_bytesize; + mp_limb_t hi_limb; + int n_bytes_in_hi_limb; + + if (stream == 0) + stream = stdout; + + if (xsize == 0) + { + for (i = 4 - 1; i >= 0; i--) + fputc (0, stream); + return ferror (stream) ? 0 : 4; + } + + hi_limb = xp[xsize - 1]; + for (i = BYTES_PER_MP_LIMB - 1; i > 0; i--) + { + if ((hi_limb >> i * BITS_PER_CHAR) != 0) + break; + } + n_bytes_in_hi_limb = i + 1; + out_bytesize = BYTES_PER_MP_LIMB * (xsize - 1) + n_bytes_in_hi_limb; + if (x->_mp_size < 0) + out_bytesize = -out_bytesize; + + /* Make the size 4 bytes on all machines, to make the format portable. */ + for (i = 4 - 1; i >= 0; i--) + fputc ((out_bytesize >> (i * BITS_PER_CHAR)) % (1 << BITS_PER_CHAR), + stream); + + /* Output from the most significant limb to the least significant limb, + with each limb also output in decreasing significance order. */ + + /* Output the most significant limb separately, since we will only + output some of its bytes. */ + for (i = n_bytes_in_hi_limb - 1; i >= 0; i--) + fputc ((hi_limb >> (i * BITS_PER_CHAR)) % (1 << BITS_PER_CHAR), stream); + + /* Output the remaining limbs. */ + for (s = xsize - 2; s >= 0; s--) + { + mp_limb_t x_limb; + + x_limb = xp[s]; + for (i = BYTES_PER_MP_LIMB - 1; i >= 0; i--) + fputc ((x_limb >> (i * BITS_PER_CHAR)) % (1 << BITS_PER_CHAR), stream); + } + return ferror (stream) ? 0 : ABS (out_bytesize) + 4; +} diff --git a/ghc/rts/gmp/mpz/out_str.c b/ghc/rts/gmp/mpz/out_str.c new file mode 100644 index 0000000..909f533 --- /dev/null +++ b/ghc/rts/gmp/mpz/out_str.c @@ -0,0 +1,108 @@ +/* mpz_out_str(stream, base, integer) -- Output to STREAM the multi prec. + integer INTEGER in base BASE. + +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 +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 +License for more details. + +You should have received a copy of the GNU Library 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 +#include "gmp.h" +#include "gmp-impl.h" + +size_t +#if __STDC__ +mpz_out_str (FILE *stream, int base, mpz_srcptr x) +#else +mpz_out_str (stream, base, x) + FILE *stream; + int base; + mpz_srcptr x; +#endif +{ + mp_ptr xp; + mp_size_t x_size = x->_mp_size; + unsigned char *str; + size_t str_size; + size_t i; + size_t written; + char *num_to_text; + TMP_DECL (marker); + + if (stream == 0) + stream = stdout; + + if (base >= 0) + { + if (base == 0) + base = 10; + num_to_text = "0123456789abcdefghijklmnopqrstuvwxyz"; + } + else + { + base = -base; + num_to_text = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ"; + } + + if (x_size == 0) + { + fputc ('0', stream); + return ferror (stream) ? 0 : 1; + } + + written = 0; + + if (x_size < 0) + { + fputc ('-', stream); + x_size = -x_size; + written = 1; + } + + TMP_MARK (marker); + str_size = ((size_t) (x_size * BITS_PER_MP_LIMB + * __mp_bases[base].chars_per_bit_exactly)) + 3; + str = (unsigned char *) TMP_ALLOC (str_size); + + /* Move the number to convert into temporary space, since mpn_get_str + clobbers its argument + needs one extra high limb.... */ + xp = (mp_ptr) TMP_ALLOC ((x_size + 1) * BYTES_PER_MP_LIMB); + MPN_COPY (xp, x->_mp_d, x_size); + + str_size = mpn_get_str (str, base, xp, x_size); + + /* mpn_get_str might make some leading zeros. Skip them. */ + while (*str == 0) + { + str_size--; + str++; + } + + /* Translate to printable chars. */ + for (i = 0; i < str_size; i++) + str[i] = num_to_text[str[i]]; + str[str_size] = 0; + + { + size_t fwret; + fwret = fwrite ((char *) str, 1, str_size, stream); + written += fwret; + } + + TMP_FREE (marker); + return ferror (stream) ? 0 : written; +} diff --git a/ghc/rts/gmp/mpz/perfsqr.c b/ghc/rts/gmp/mpz/perfsqr.c new file mode 100644 index 0000000..cdf1b5a --- /dev/null +++ b/ghc/rts/gmp/mpz/perfsqr.c @@ -0,0 +1,41 @@ +/* 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. + +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 +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 +License for more details. + +You should have received a copy of the GNU Library 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_perfect_square_p (mpz_srcptr a) +#else +mpz_perfect_square_p (a) + mpz_srcptr a; +#endif +{ + mp_size_t asize = a->_mp_size; + + /* No negative numbers are perfect squares. */ + if (asize < 0) + return 0; + + 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 new file mode 100644 index 0000000..a979380 --- /dev/null +++ b/ghc/rts/gmp/mpz/popcount.c @@ -0,0 +1,42 @@ +/* mpz_popcount(mpz_ptr op) -- Population count of OP. If the operand is + negative, return ~0 (a novel representation of infinity). + +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 +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 +License for more details. + +You should have received a copy of the GNU Library 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_popcount (mpz_srcptr u) +#else +mpz_popcount (u) + mpz_srcptr u; +#endif +{ + mp_size_t usize; + + usize = u->_mp_size; + + if ((usize) < 0) + return ~ (unsigned long int) 0; + + return mpn_popcount (u->_mp_d, usize); +} diff --git a/ghc/rts/gmp/mpz/pow_ui.c b/ghc/rts/gmp/mpz/pow_ui.c new file mode 100644 index 0000000..d8cf7a6 --- /dev/null +++ b/ghc/rts/gmp/mpz/pow_ui.c @@ -0,0 +1,129 @@ +/* mpz_pow_ui(res, base, exp) -- Set RES to BASE**EXP. + +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 +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 +License for more details. + +You should have received a copy of the GNU Library 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. */ + +#ifdef BERKELEY_MP +#include "mp.h" +#endif +#include "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" + +#ifndef BERKELEY_MP +void +#if __STDC__ +mpz_pow_ui (mpz_ptr r, mpz_srcptr b, unsigned long int e) +#else +mpz_pow_ui (r, b, e) + mpz_ptr r; + mpz_srcptr b; + unsigned long int e; +#endif +#else /* BERKELEY_MP */ +void +#if __STDC__ +rpow (const MINT *b, signed short int e, MINT *r) +#else +rpow (b, e, r) + const MINT *b; + signed short int e; + MINT *r; +#endif +#endif /* BERKELEY_MP */ +{ + mp_ptr rp, bp, tp, xp; + mp_size_t rsize, bsize; + int cnt, i; + mp_limb_t blimb; + TMP_DECL (marker); + + bsize = ABS (b->_mp_size); + + /* 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) + { + r->_mp_d[0] = 1; + r->_mp_size = 1; + return; + } + if (bsize == 0 +#ifdef BERKELEY_MP + || e < 0 +#endif + ) + { + r->_mp_size = 0; + return; + } + + bp = b->_mp_d; + + blimb = bp[bsize - 1]; + if (bsize == 1 && blimb < 0x100) + { + /* 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); + } + else + { + /* Over-estimate space requirements somewhat. */ + count_leading_zeros (cnt, blimb); + rsize = bsize * e - cnt * e / BITS_PER_MP_LIMB + 1; + } + + 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); + + MPN_COPY (rp, bp, bsize); + rsize = bsize; + 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; + + if ((e & ((mp_limb_t) 1 << i)) != 0) + { + rsize = rsize + bsize - (mpn_mul (tp, rp, rsize, bp, bsize) == 0); + xp = tp; tp = rp; rp = xp; + } + } + + /* Now then we know the exact space requirements, reallocate if + necessary. */ + if (r->_mp_alloc < rsize) + _mpz_realloc (r, rsize); + + MPN_COPY (r->_mp_d, rp, rsize); + r->_mp_size = (e & 1) == 0 || b->_mp_size >= 0 ? rsize : -rsize; + TMP_FREE (marker); +} diff --git a/ghc/rts/gmp/mpz/powm.c b/ghc/rts/gmp/mpz/powm.c new file mode 100644 index 0000000..5dcd1b1 --- /dev/null +++ b/ghc/rts/gmp/mpz/powm.c @@ -0,0 +1,276 @@ +/* mpz_powm(res,base,exp,mod) -- Set RES to (base**exp) mod MOD. + +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 +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 +License for more details. + +You should have received a copy of the GNU Library 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" + +#ifndef BERKELEY_MP +void +#if __STDC__ +mpz_powm (mpz_ptr res, mpz_srcptr base, mpz_srcptr exp, mpz_srcptr mod) +#else +mpz_powm (res, base, exp, mod) + mpz_ptr res; + mpz_srcptr base; + mpz_srcptr exp; + mpz_srcptr mod; +#endif +#else /* BERKELEY_MP */ +void +#if __STDC__ +pow (mpz_srcptr base, mpz_srcptr exp, mpz_srcptr mod, mpz_ptr res) +#else +pow (base, exp, mod, res) + mpz_srcptr base; + mpz_srcptr exp; + 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; + + rp = res->_mp_d; + ep = exp->_mp_d; + + if (msize == 0) + msize = 1 / msize; /* provoke a signal */ + + if (esize == 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; + 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. */ + + 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); + + bsize = ABS (base->_mp_size); + if (bsize > msize) + { + /* The base is larger than the module. Reduce it. */ + + /* 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); + } + else + bp = base->_mp_d; + + if (bsize == 0) + { + res->_mp_size = 0; + TMP_FREE (marker); + return; + } + + if (res->_mp_alloc < size) + { + /* 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) + { + free_me = rp; + free_me_size = res->_mp_alloc; + } + 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); + } + if (rp == ep) + { + /* 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); + } + if (rp == mp) + { + /* 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); + } + } + + 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. */ + + for (;;) + { + while (c != 0) + { + mp_ptr tp; + mp_size_t xsize; + + mpn_mul_n (xp, rp, rp, rsize); + xsize = 2 * rsize; + if (xsize > msize) + { + mpn_divmod (xp + msize, xp, xsize, mp, msize); + xsize = msize; + } + + tp = rp; rp = xp; xp = tp; + rsize = xsize; + + if ((mp_limb_signed_t) e < 0) + { + 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; + } + 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 (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); + TMP_FREE (marker); +} diff --git a/ghc/rts/gmp/mpz/powm_ui.c b/ghc/rts/gmp/mpz/powm_ui.c new file mode 100644 index 0000000..596815a --- /dev/null +++ b/ghc/rts/gmp/mpz/powm_ui.c @@ -0,0 +1,234 @@ +/* mpz_powm_ui(res,base,exp,mod) -- Set RES to (base**exp) mod MOD. + +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 +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 +License for more details. + +You should have received a copy of the GNU Library 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_powm_ui (mpz_ptr res, mpz_srcptr base, unsigned long int exp, mpz_srcptr mod) +#else +mpz_powm_ui (res, base, exp, mod) + mpz_ptr res; + mpz_srcptr base; + unsigned long int exp; + mpz_srcptr mod; +#endif +{ + mp_ptr rp, mp, bp; + mp_size_t 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); + + msize = ABS (mod->_mp_size); + size = 2 * msize; + + rp = res->_mp_d; + + if (msize == 0) + msize = 1 / msize; /* provoke a signal */ + + if (exp == 0) + { + rp[0] = 1; + res->_mp_size = (msize == 1 && (mod->_mp_d)[0] == 1) ? 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. */ + + 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); + + bsize = ABS (base->_mp_size); + if (bsize > msize) + { + /* The base is larger than the module. Reduce it. */ + + /* 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); + } + else + bp = base->_mp_d; + + if (bsize == 0) + { + res->_mp_size = 0; + TMP_FREE (marker); + return; + } + + if (res->_mp_alloc < size) + { + /* 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 == mp || rp == bp) + { + free_me = rp; + free_me_size = res->_mp_alloc; + } + 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); + } + if (rp == mp) + { + /* 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); + } + } + + MPN_COPY (rp, bp, bsize); + rsize = bsize; + + { + 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 = (exp & 1) && base->_mp_size < 0; + + e = exp; + 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. */ + + while (c != 0) + { + mp_ptr tp; + mp_size_t xsize; + + mpn_mul_n (xp, rp, rp, rsize); + xsize = 2 * rsize; + if (xsize > msize) + { + mpn_divmod (xp + msize, xp, xsize, mp, msize); + xsize = msize; + } + + tp = rp; rp = xp; xp = tp; + rsize = xsize; + + if ((mp_limb_signed_t) e < 0) + { + 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; + } + e <<= 1; + c--; + } + + /* 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); + } + + res->_mp_size = negative_result == 0 ? rsize : -rsize; + + if (free_me != NULL) + (*_mp_free_func) (free_me, free_me_size * BYTES_PER_MP_LIMB); + TMP_FREE (marker); +} diff --git a/ghc/rts/gmp/mpz/pprime_p.c b/ghc/rts/gmp/mpz/pprime_p.c new file mode 100644 index 0000000..494de14 --- /dev/null +++ b/ghc/rts/gmp/mpz/pprime_p.c @@ -0,0 +1,115 @@ +/* mpz_probab_prime_p -- + 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. + +Copyright (C) 1991, 1993, 1994 Free Software Foundation, Inc. +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 +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 +License for more details. + +You should have received a copy of the GNU Library 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" + +static int +possibly_prime (n, n_minus_1, x, y, q, k) + mpz_srcptr n; + mpz_srcptr n_minus_1; + mpz_ptr x; + mpz_ptr y; + mpz_srcptr q; + unsigned long int k; +{ + unsigned long int i; + + /* find random x s.t. 1 < x < n */ + do + { + mpz_random (x, mpz_size (n)); + mpz_mmod (x, x, n); + } + while (mpz_cmp_ui (x, 1L) <= 0); + + mpz_powm (y, x, q, n); + + if (mpz_cmp_ui (y, 1L) == 0 || mpz_cmp (y, n_minus_1) == 0) + return 1; + + for (i = 1; i < k; i++) + { + mpz_powm_ui (y, y, 2L, n); + if (mpz_cmp (y, n_minus_1) == 0) + return 1; + if (mpz_cmp_ui (y, 1L) == 0) + return 0; + } + return 0; +} + +int +#if __STDC__ +mpz_probab_prime_p (mpz_srcptr m, int reps) +#else +mpz_probab_prime_p (m, reps) + mpz_srcptr m; + int reps; +#endif +{ + mpz_t n, n_minus_1, x, y, q; + int i, is_prime; + unsigned long int k; + + mpz_init (n); + /* Take the absolute value of M, to handle positive and negative primes. */ + mpz_abs (n, m); + + if (mpz_cmp_ui (n, 3L) <= 0) + { + mpz_clear (n); + return mpz_cmp_ui (n, 1L) > 0; + } + + if ((mpz_get_ui (n) & 1) == 0) + { + mpz_clear (n); + return 0; /* even */ + } + + mpz_init (n_minus_1); + mpz_sub_ui (n_minus_1, n, 1L); + mpz_init (x); + mpz_init (y); + + /* 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); + + 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); + return is_prime; +} diff --git a/ghc/rts/gmp/mpz/random.c b/ghc/rts/gmp/mpz/random.c new file mode 100644 index 0000000..ab41eef --- /dev/null +++ b/ghc/rts/gmp/mpz/random.c @@ -0,0 +1,56 @@ +/* mpz_random -- Generate a random mpz_t of specified size. + This function is non-portable and generates poor random numbers. + +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 +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 +License for more details. + +You should have received a copy of the GNU Library 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 "urandom.h" + +void +#if __STDC__ +mpz_random (mpz_ptr x, mp_size_t size) +#else +mpz_random (x, size) + mpz_ptr x; + mp_size_t size; +#endif +{ + mp_size_t i; + mp_limb_t ran; + mp_ptr xp; + mp_size_t abs_size; + + abs_size = ABS (size); + + if (x->_mp_alloc < abs_size) + _mpz_realloc (x, abs_size); + + xp = x->_mp_d; + + for (i = 0; i < abs_size; i++) + { + ran = urandom (); + xp[i] = ran; + } + + MPN_NORMALIZE (xp, abs_size); + x->_mp_size = size < 0 ? -abs_size : abs_size; +} diff --git a/ghc/rts/gmp/mpz/random2.c b/ghc/rts/gmp/mpz/random2.c new file mode 100644 index 0000000..725a8b4 --- /dev/null +++ b/ghc/rts/gmp/mpz/random2.c @@ -0,0 +1,48 @@ +/* mpz_random2 -- Generate a positive random mpz_t of specified size, with + long runs of consecutive ones and zeros in the binary representation. + Meant for testing of other MP routines. + +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 +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 +License for more details. + +You should have received a copy of the GNU Library 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_random2 (mpz_ptr x, mp_size_t size) +#else +mpz_random2 (x, size) + mpz_ptr x; + mp_size_t size; +#endif +{ + mp_size_t abs_size; + + abs_size = ABS (size); + if (abs_size != 0) + { + if (x->_mp_alloc < abs_size) + _mpz_realloc (x, abs_size); + + mpn_random2 (x->_mp_d, abs_size); + } + + x->_mp_size = size; +}