From: simonm Date: Fri, 5 Jun 1998 14:38:00 +0000 (+0000) Subject: [project @ 1998-06-05 14:37:59 by simonm] X-Git-Tag: Approx_2487_patches~607 X-Git-Url: http://git.megacz.com/?a=commitdiff_plain;h=47adae51c9f4e0a34b760a9d15a8d6b567acd175;p=ghc-hetmet.git [project @ 1998-06-05 14:37:59 by simonm] Import GMP 2.0.2 --- diff --git a/ghc/rts/gmp/mpz/abs.c b/ghc/rts/gmp/mpz/abs.c new file mode 100644 index 0000000..080cac6 --- /dev/null +++ b/ghc/rts/gmp/mpz/abs.c @@ -0,0 +1,51 @@ +/* mpz_abs(dst, src) -- Assign the absolute 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_abs (mpz_ptr w, mpz_srcptr u) +#else +mpz_abs (w, u) + mpz_ptr w; + mpz_srcptr u; +#endif +{ + mp_ptr wp, up; + mp_size_t size; + + size = ABS (u->_mp_size); + + if (u != w) + { + 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 = size; +} diff --git a/ghc/rts/gmp/mpz/add.c b/ghc/rts/gmp/mpz/add.c new file mode 100644 index 0000000..10dd970 --- /dev/null +++ b/ghc/rts/gmp/mpz/add.c @@ -0,0 +1,120 @@ +/* mpz_add -- Add 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_add (mpz_ptr w, mpz_srcptr u, mpz_srcptr v) +#else +mpz_add (w, u, v) + mpz_ptr w; + mpz_srcptr u; + mpz_srcptr v; +#endif +#else /* BERKELEY_MP */ +void +#if __STDC__ +madd (mpz_srcptr u, mpz_srcptr v, mpz_ptr w) +#else +madd (u, v, w) + mpz_srcptr u; + mpz_srcptr v; + mpz_ptr w; +#endif +#endif /* BERKELEY_MP */ +{ + mp_srcptr up, vp; + mp_ptr wp; + mp_size_t usize, vsize, wsize; + mp_size_t abs_usize; + mp_size_t abs_vsize; + + usize = u->_mp_size; + vsize = v->_mp_size; + abs_usize = ABS (usize); + abs_vsize = ABS (vsize); + + if (abs_usize < abs_vsize) + { + /* Swap U and V. */ + {const __mpz_struct *t = u; u = v; v = t;} + {mp_size_t t = usize; usize = vsize; vsize = t;} + {mp_size_t t = abs_usize; abs_usize = abs_vsize; abs_vsize = t;} + } + + /* True: ABS_USIZE >= ABS_VSIZE. */ + + /* If not space for w (and possible carry), increase space. */ + wsize = abs_usize + 1; + if (w->_mp_alloc < wsize) + _mpz_realloc (w, wsize); + + /* These must be after realloc (u or v may be the same as w). */ + up = u->_mp_d; + vp = v->_mp_d; + wp = w->_mp_d; + + if ((usize ^ vsize) < 0) + { + /* U and V have different sign. Need to compare them to determine + which operand to subtract from which. */ + + /* This test is right since ABS_USIZE >= ABS_VSIZE. */ + if (abs_usize != abs_vsize) + { + mpn_sub (wp, up, abs_usize, vp, abs_vsize); + wsize = abs_usize; + MPN_NORMALIZE (wp, wsize); + if (usize < 0) + wsize = -wsize; + } + else if (mpn_cmp (up, vp, abs_usize) < 0) + { + mpn_sub_n (wp, vp, up, abs_usize); + wsize = abs_usize; + MPN_NORMALIZE (wp, wsize); + if (usize >= 0) + wsize = -wsize; + } + else + { + mpn_sub_n (wp, up, vp, abs_usize); + wsize = abs_usize; + MPN_NORMALIZE (wp, wsize); + if (usize < 0) + wsize = -wsize; + } + } + else + { + /* U and V have same sign. Add them. */ + mp_limb_t cy_limb = mpn_add (wp, up, abs_usize, vp, abs_vsize); + wp[abs_usize] = cy_limb; + wsize = abs_usize + cy_limb; + if (usize < 0) + wsize = -wsize; + } + + w->_mp_size = wsize; +} diff --git a/ghc/rts/gmp/mpz/add_ui.c b/ghc/rts/gmp/mpz/add_ui.c new file mode 100644 index 0000000..a1e4306 --- /dev/null +++ b/ghc/rts/gmp/mpz/add_ui.c @@ -0,0 +1,84 @@ +/* mpz_add_ui -- Add an mpz_t and an unsigned one-word 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_add_ui (mpz_ptr w, mpz_srcptr u, unsigned long int v) +#else +mpz_add_ui (w, u, v) + mpz_ptr w; + mpz_srcptr u; + unsigned long int v; +#endif +{ + mp_srcptr up; + mp_ptr wp; + mp_size_t usize, wsize; + mp_size_t abs_usize; + + usize = u->_mp_size; + abs_usize = ABS (usize); + + /* If not space for W (and possible carry), increase space. */ + wsize = abs_usize + 1; + if (w->_mp_alloc < wsize) + _mpz_realloc (w, wsize); + + /* These must be after realloc (U may be the same as W). */ + up = u->_mp_d; + wp = w->_mp_d; + + if (abs_usize == 0) + { + wp[0] = v; + w->_mp_size = v != 0; + return; + } + + if (usize >= 0) + { + mp_limb_t cy; + cy = mpn_add_1 (wp, up, abs_usize, v); + wp[abs_usize] = cy; + wsize = abs_usize + cy; + } + else + { + /* The signs are different. Need exact comparison to determine + which operand to subtract from which. */ + if (abs_usize == 1 && up[0] < v) + { + wp[0] = v - up[0]; + wsize = 1; + } + else + { + mpn_sub_1 (wp, up, abs_usize, v); + /* Size can decrease with at most one limb. */ + wsize = -(abs_usize - (wp[abs_usize - 1] == 0)); + } + } + + w->_mp_size = wsize; +} diff --git a/ghc/rts/gmp/mpz/and.c b/ghc/rts/gmp/mpz/and.c new file mode 100644 index 0000000..838d4b1 --- /dev/null +++ b/ghc/rts/gmp/mpz/and.c @@ -0,0 +1,278 @@ +/* mpz_and -- Logical and. + +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_and (mpz_ptr res, mpz_srcptr op1, mpz_srcptr op2) +#else +mpz_and (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) + { + res_size = MIN (op1_size, op2_size); + /* 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; + + /* Handle allocation, now then we know exactly how much space is + needed for the result. */ + if (res->_mp_alloc < res_size) + { + _mpz_realloc (res, res_size); + op1_ptr = op1->_mp_d; + op2_ptr = op2->_mp_d; + res_ptr = res->_mp_d; + } + + /* Second loop computes the real result. */ + for (i = res_size - 1; i >= 0; i--) + res_ptr[i] = op1_ptr[i] & op2_ptr[i]; + + 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; + mp_size_t res_alloc; + + /* 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 */ + + /* It might seem as we could end up with an (invalid) result with + a leading zero-limb here when one of the operands is of the + type 1,,0,,..,,.0. But some analysis shows that we surely + would get carry into the zero-limb in this situation... */ + + op1_size = -op1_size; + op2_size = -op2_size; + + res_alloc = 1 + MAX (op1_size, op2_size); + + opx = (mp_ptr) TMP_ALLOC (op1_size * BYTES_PER_MP_LIMB); + mpn_sub_1 (opx, op1_ptr, op1_size, (mp_limb_t) 1); + op1_ptr = opx; + + opx = (mp_ptr) TMP_ALLOC (op2_size * BYTES_PER_MP_LIMB); + mpn_sub_1 (opx, op2_ptr, op2_size, (mp_limb_t) 1); + op2_ptr = opx; + + if (res->_mp_alloc < res_alloc) + { + _mpz_realloc (res, res_alloc); + res_ptr = res->_mp_d; + /* Don't re-read OP1_PTR and OP2_PTR. They point to + temporary space--never to the space RES->_mp_D used + to point to before reallocation. */ + } + + if (op1_size >= op2_size) + { + MPN_COPY (res_ptr + op2_size, op1_ptr + op2_size, + op1_size - op2_size); + for (i = op2_size - 1; i >= 0; i--) + res_ptr[i] = op1_ptr[i] | op2_ptr[i]; + res_size = op1_size; + } + else + { + MPN_COPY (res_ptr + op1_size, op2_ptr + op1_size, + op2_size - op1_size); + for (i = op1_size - 1; i >= 0; i--) + res_ptr[i] = op1_ptr[i] | op2_ptr[i]; + res_size = op2_size; + } + + cy = mpn_add_1 (res_ptr, res_ptr, res_size, (mp_limb_t) 1); + if (cy) + { + res_ptr[res_size] = cy; + res_size++; + } + + res->_mp_size = -res_size; + TMP_FREE (marker); + return; + } + else + { + /* We should compute -OP1 & OP2. Swap OP1 and OP2 and fall + through to the code that handles OP1 & -OP2. */ + {mpz_srcptr 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;} + } + + } + + { +#if ANDNEW + mp_size_t op2_lim; + mp_size_t count; + + /* OP2 must be negated as with infinite precision. + + Scan from the low end for a non-zero limb. The first non-zero + limb is simply negated (two's complement). Any subsequent + limbs are one's complemented. Of course, we don't need to + handle more limbs than there are limbs in the other, positive + operand as the result for those limbs is going to become zero + anyway. */ + + /* Scan for the least significant. non-zero OP2 limb, and zero the + result meanwhile for those limb positions. (We will surely + find a non-zero limb, so we can write the loop with one + termination condition only.) */ + for (i = 0; op2_ptr[i] == 0; i++) + res_ptr[i] = 0; + op2_lim = i; + + op2_size = -op2_size; + + if (op1_size <= op2_size) + { + /* The ones-extended OP2 is >= than the zero-extended OP1. + RES_SIZE <= OP1_SIZE. Find the exact size. */ + for (i = op1_size - 1; i > op2_lim; i--) + if ((op1_ptr[i] & ~op2_ptr[i]) != 0) + break; + res_size = i + 1; + for (i = res_size - 1; i > op2_lim; i--) + res_ptr[i] = op1_ptr[i] & ~op2_ptr[i]; + res_ptr[op2_lim] = op1_ptr[op2_lim] & -op2_ptr[op2_lim]; + /* Yes, this *can* happen! */ + MPN_NORMALIZE (res_ptr, res_size); + } + else + { + /* The ones-extended OP2 is < than the zero-extended OP1. + RES_SIZE == OP1_SIZE, since OP1 is normalized. */ + res_size = op1_size; + MPN_COPY (res_ptr + op2_size, op1_ptr + op2_size, op1_size - op2_size); + for (i = op2_size - 1; i > op2_lim; i--) + res_ptr[i] = op1_ptr[i] & ~op2_ptr[i]; + res_ptr[op2_lim] = op1_ptr[op2_lim] & -op2_ptr[op2_lim]; + } + + res->_mp_size = res_size; +#else + + /* OP1 is positive and zero-extended, + OP2 is negative and ones-extended. + The result will be positive. + OP1 & -OP2 = OP1 & ~(OP2 - 1). */ + + mp_ptr opx; + + op2_size = -op2_size; + opx = (mp_ptr) TMP_ALLOC (op2_size * BYTES_PER_MP_LIMB); + mpn_sub_1 (opx, op2_ptr, op2_size, (mp_limb_t) 1); + op2_ptr = opx; + + if (op1_size > op2_size) + { + /* The result has the same size as OP1, since OP1 is normalized + and longer than the ones-extended OP2. */ + res_size = op1_size; + + /* Handle allocation, now then we know exactly how much space is + needed for the result. */ + if (res->_mp_alloc < res_size) + { + _mpz_realloc (res, res_size); + res_ptr = res->_mp_d; + op1_ptr = op1->_mp_d; + /* Don't re-read OP2_PTR. It points to temporary space--never + to the space RES->_mp_D used to point to before reallocation. */ + } + + MPN_COPY (res_ptr + op2_size, op1_ptr + op2_size, + res_size - op2_size); + for (i = op2_size - 1; i >= 0; i--) + res_ptr[i] = op1_ptr[i] & ~op2_ptr[i]; + + res->_mp_size = res_size; + } + else + { + /* Find out the exact result size. Ignore the high limbs of OP2, + OP1 is zero-extended and would make the result zero. */ + for (i = op1_size - 1; i >= 0; i--) + if ((op1_ptr[i] & ~op2_ptr[i]) != 0) + break; + res_size = i + 1; + + /* Handle allocation, now then we know exactly how much space is + needed for the result. */ + if (res->_mp_alloc < res_size) + { + _mpz_realloc (res, res_size); + res_ptr = res->_mp_d; + op1_ptr = op1->_mp_d; + /* Don't re-read OP2_PTR. It points to temporary space--never + to the space RES->_mp_D used to point to before reallocation. */ + } + + for (i = res_size - 1; i >= 0; i--) + res_ptr[i] = op1_ptr[i] & ~op2_ptr[i]; + + res->_mp_size = res_size; + } +#endif + } + TMP_FREE (marker); +} diff --git a/ghc/rts/gmp/mpz/array_init.c b/ghc/rts/gmp/mpz/array_init.c new file mode 100644 index 0000000..8b2e85c --- /dev/null +++ b/ghc/rts/gmp/mpz/array_init.c @@ -0,0 +1,48 @@ +/* mpz_array_init (array, array_size, size_per_elem) -- + +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_array_init (mpz_ptr arr, mp_size_t arr_size, mp_size_t nbits) +#else +mpz_array_init (arr, arr_size, nbits) + mpz_ptr arr; + mp_size_t arr_size; + mp_size_t nbits; +#endif +{ + register mp_ptr p; + register size_t i; + mp_size_t nlimbs; + + nlimbs = (nbits + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB; + p = (mp_ptr) (*_mp_allocate_func) (arr_size * nlimbs * BYTES_PER_MP_LIMB); + + for (i = 0; i < arr_size; i++) + { + arr[i]._mp_alloc = nlimbs + 1; /* Yes, lie a little... */ + arr[i]._mp_size = 0; + arr[i]._mp_d = p + i * nlimbs; + } +} diff --git a/ghc/rts/gmp/mpz/cdiv_q.c b/ghc/rts/gmp/mpz/cdiv_q.c new file mode 100644 index 0000000..860a232 --- /dev/null +++ b/ghc/rts/gmp/mpz/cdiv_q.c @@ -0,0 +1,51 @@ +/* mpz_cdiv_q -- Division rounding the quotient towards +infinity. The + remainder gets the opposite sign as the denominator. + +Copyright (C) 1994, 1995, 1996 Free Software Foundation, Inc. + +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_cdiv_q (mpz_ptr quot, mpz_srcptr dividend, mpz_srcptr divisor) +#else +mpz_cdiv_q (quot, dividend, divisor) + mpz_ptr quot; + mpz_srcptr dividend; + mpz_srcptr divisor; +#endif +{ + mp_size_t dividend_size = dividend->_mp_size; + mp_size_t divisor_size = divisor->_mp_size; + mpz_t rem; + TMP_DECL (marker); + + TMP_MARK (marker); + + MPZ_TMP_INIT (rem, 1 + ABS (dividend_size)); + + mpz_tdiv_qr (quot, rem, dividend, divisor); + + if ((divisor_size ^ dividend_size) >= 0 && rem->_mp_size != 0) + mpz_add_ui (quot, quot, 1L); + + TMP_FREE (marker); +} diff --git a/ghc/rts/gmp/mpz/cdiv_q_ui.c b/ghc/rts/gmp/mpz/cdiv_q_ui.c new file mode 100644 index 0000000..7b6cfd7 --- /dev/null +++ b/ghc/rts/gmp/mpz/cdiv_q_ui.c @@ -0,0 +1,64 @@ +/* mpz_cdiv_q_ui -- Division rounding the quotient towards +infinity. The + remainder gets the opposite sign as the denominator. In order to make it + always fit into the return type, the negative of the true remainder is + returned. + +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_cdiv_q_ui (mpz_ptr quot, mpz_srcptr dividend, unsigned long int divisor) +#else +mpz_cdiv_q_ui (quot, dividend, divisor) + mpz_ptr quot; + mpz_srcptr dividend; + unsigned long int divisor; +#endif +{ + mp_size_t dividend_size; + mp_size_t size; + mp_ptr quot_ptr; + mp_limb_t remainder_limb; + + dividend_size = dividend->_mp_size; + size = ABS (dividend_size); + + if (quot->_mp_alloc < size) + _mpz_realloc (quot, size); + + quot_ptr = quot->_mp_d; + + remainder_limb = mpn_divmod_1 (quot_ptr, dividend->_mp_d, size, + (mp_limb_t) divisor); + + if (remainder_limb != 0 && dividend_size >= 0) + { + mpn_add_1 (quot_ptr, quot_ptr, size, (mp_limb_t) 1); + remainder_limb = divisor - remainder_limb; + } + + size -= size != 0 && quot_ptr[size - 1] == 0; + quot->_mp_size = dividend_size >= 0 ? size : -size; + + return remainder_limb; +} diff --git a/ghc/rts/gmp/mpz/cdiv_qr.c b/ghc/rts/gmp/mpz/cdiv_qr.c new file mode 100644 index 0000000..bf7d6da --- /dev/null +++ b/ghc/rts/gmp/mpz/cdiv_qr.c @@ -0,0 +1,62 @@ +/* mpz_cdiv_qr -- Division rounding the quotient towards +infinity. The + remainder gets the opposite sign as the denominator. + +Copyright (C) 1994, 1995, 1996 Free Software Foundation, Inc. + +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_cdiv_qr (mpz_ptr quot, mpz_ptr rem, mpz_srcptr dividend, mpz_srcptr divisor) +#else +mpz_cdiv_qr (quot, rem, dividend, divisor) + mpz_ptr quot; + 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 quotient and + remainder have been preliminary calculated. We have to copy it to + temporary space if it's the same variable as either QUOT or REM. */ + if (quot == divisor || rem == divisor) + { + MPZ_TMP_INIT (temp_divisor, ABS (divisor_size)); + mpz_set (temp_divisor, divisor); + divisor = temp_divisor; + } + + mpz_tdiv_qr (quot, rem, dividend, divisor); + + if ((divisor_size ^ dividend->_mp_size) >= 0 && rem->_mp_size != 0) + { + mpz_add_ui (quot, quot, 1L); + mpz_sub (rem, rem, divisor); + } + + TMP_FREE (marker); +} diff --git a/ghc/rts/gmp/mpz/cdiv_qr_ui.c b/ghc/rts/gmp/mpz/cdiv_qr_ui.c new file mode 100644 index 0000000..a780e77 --- /dev/null +++ b/ghc/rts/gmp/mpz/cdiv_qr_ui.c @@ -0,0 +1,68 @@ +/* mpz_cdiv_qr_ui -- Division rounding the quotient towards +infinity. The + remainder gets the opposite sign as the denominator. In order to make it + always fit into the return type, the negative of the true remainder is + returned. + +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_cdiv_qr_ui (mpz_ptr quot, mpz_ptr rem, mpz_srcptr dividend, unsigned long int divisor) +#else +mpz_cdiv_qr_ui (quot, rem, dividend, divisor) + mpz_ptr quot; + mpz_ptr rem; + mpz_srcptr dividend; + unsigned long int divisor; +#endif +{ + mp_size_t dividend_size; + mp_size_t size; + mp_ptr quot_ptr; + mp_limb_t remainder_limb; + + dividend_size = dividend->_mp_size; + size = ABS (dividend_size); + + if (quot->_mp_alloc < size) + _mpz_realloc (quot, size); + + quot_ptr = quot->_mp_d; + + remainder_limb = mpn_divmod_1 (quot_ptr, dividend->_mp_d, size, + (mp_limb_t) divisor); + + if (remainder_limb != 0 && dividend_size >= 0) + { + mpn_add_1 (quot_ptr, quot_ptr, size, (mp_limb_t) 1); + remainder_limb = divisor - remainder_limb; + } + + size -= size != 0 && quot_ptr[size - 1] == 0; + quot->_mp_size = dividend_size >= 0 ? size : -size; + + rem->_mp_d[0] = remainder_limb; + rem->_mp_size = -(remainder_limb != 0); + + return remainder_limb; +} diff --git a/ghc/rts/gmp/mpz/cdiv_r.c b/ghc/rts/gmp/mpz/cdiv_r.c new file mode 100644 index 0000000..d34d138 --- /dev/null +++ b/ghc/rts/gmp/mpz/cdiv_r.c @@ -0,0 +1,59 @@ +/* mpz_cdiv_r -- Division rounding the quotient towards +infinity. The + remainder gets the opposite 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" + +void +#if __STDC__ +mpz_cdiv_r (mpz_ptr rem, mpz_srcptr dividend, mpz_srcptr divisor) +#else +mpz_cdiv_r (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 ((divisor_size ^ dividend->_mp_size) >= 0 && rem->_mp_size != 0) + mpz_sub (rem, rem, divisor); + + TMP_FREE (marker); +} diff --git a/ghc/rts/gmp/mpz/cdiv_r_ui.c b/ghc/rts/gmp/mpz/cdiv_r_ui.c new file mode 100644 index 0000000..757a3f5 --- /dev/null +++ b/ghc/rts/gmp/mpz/cdiv_r_ui.c @@ -0,0 +1,54 @@ +/* mpz_cdiv_r_ui -- Division rounding the quotient towards +infinity. The + remainder gets the opposite sign as the denominator. In order to make it + always fit into the return type, the negative of the true remainder is + returned. + +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_cdiv_r_ui (mpz_ptr rem, mpz_srcptr dividend, unsigned long int divisor) +#else +mpz_cdiv_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/cdiv_ui.c b/ghc/rts/gmp/mpz/cdiv_ui.c new file mode 100644 index 0000000..df841ed --- /dev/null +++ b/ghc/rts/gmp/mpz/cdiv_ui.c @@ -0,0 +1,50 @@ +/* mpz_cdiv_ui -- Division rounding the quotient towards +infinity. The + remainder gets the opposite sign as the denominator. In order to make it + always fit into the return type, the negative of the true remainder is + returned. + +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_cdiv_ui (mpz_srcptr dividend, unsigned long int divisor) +#else +mpz_cdiv_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/clear.c b/ghc/rts/gmp/mpz/clear.c new file mode 100644 index 0000000..00f3cfd --- /dev/null +++ b/ghc/rts/gmp/mpz/clear.c @@ -0,0 +1,35 @@ +/* mpz_clear -- de-allocate the space occupied by the dynamic digit space of + an 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" + +void +#if __STDC__ +mpz_clear (mpz_ptr m) +#else +mpz_clear (m) + mpz_ptr m; +#endif +{ + (*_mp_free_func) (m->_mp_d, m->_mp_alloc * BYTES_PER_MP_LIMB); +} diff --git a/ghc/rts/gmp/mpz/clrbit.c b/ghc/rts/gmp/mpz/clrbit.c new file mode 100644 index 0000000..59d9565 --- /dev/null +++ b/ghc/rts/gmp/mpz/clrbit.c @@ -0,0 +1,114 @@ +/* mpz_clrbit -- clear a specified bit. + +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_clrbit (mpz_ptr d, unsigned long int bit_index) +#else +mpz_clrbit (d, bit_index) + mpz_ptr d; + unsigned long int bit_index; +#endif +{ + mp_size_t dsize = d->_mp_size; + mp_ptr dp = d->_mp_d; + mp_size_t limb_index; + + limb_index = bit_index / BITS_PER_MP_LIMB; + if (dsize >= 0) + { + if (limb_index < dsize) + { + dp[limb_index] &= ~((mp_limb_t) 1 << (bit_index % BITS_PER_MP_LIMB)); + MPN_NORMALIZE (dp, dsize); + d->_mp_size = dsize; + } + else + ; + } + else + { + mp_size_t zero_bound; + + /* Simulate two's complement arithmetic, i.e. simulate + 1. Set OP = ~(OP - 1) [with infinitely many leading ones]. + 2. clear the bit. + 3. Set OP = ~OP + 1. */ + + dsize = -dsize; + + /* No upper bound on this loop, we're sure there's a non-zero limb + sooner ot later. */ + for (zero_bound = 0; ; zero_bound++) + if (dp[zero_bound] != 0) + break; + + if (limb_index > zero_bound) + { + if (limb_index < dsize) + dp[limb_index] |= (mp_limb_t) 1 << (bit_index % BITS_PER_MP_LIMB); + else + { + /* Ugh. The bit should be cleared outside of the end of the + number. We have to increase the size of the number. */ + if (d->_mp_alloc < limb_index + 1) + { + _mpz_realloc (d, limb_index + 1); + dp = d->_mp_d; + } + MPN_ZERO (dp + dsize, limb_index - dsize); + dp[limb_index] = (mp_limb_t) 1 << (bit_index % BITS_PER_MP_LIMB); + d->_mp_size = -(limb_index + 1); + } + } + else if (limb_index == zero_bound) + { + dp[limb_index] = ((dp[limb_index] - 1) + | ((mp_limb_t) 1 << (bit_index % BITS_PER_MP_LIMB))) + 1; + if (dp[limb_index] == 0) + { + mp_size_t i; + for (i = limb_index + 1; i < dsize; i++) + { + dp[i] += 1; + if (dp[i] != 0) + goto fin; + } + /* We got carry all way out beyond the end of D. Increase + its size (and allocation if necessary). */ + dsize++; + if (d->_mp_alloc < dsize) + { + _mpz_realloc (d, dsize); + dp = d->_mp_d; + } + dp[i] = 1; + d->_mp_size = -dsize; + fin:; + } + } + else + ; + } +} diff --git a/ghc/rts/gmp/mpz/cmp.c b/ghc/rts/gmp/mpz/cmp.c new file mode 100644 index 0000000..37be334 --- /dev/null +++ b/ghc/rts/gmp/mpz/cmp.c @@ -0,0 +1,75 @@ +/* mpz_cmp(u,v) -- Compare U, V. Return postive, zero, or negative + based on if U > V, U == V, or U < V. + +Copyright (C) 1991, 1993, 1994, 1996 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" + +#ifndef BERKELEY_MP +int +#if __STDC__ +mpz_cmp (mpz_srcptr u, mpz_srcptr v) +#else +mpz_cmp (u, v) + mpz_srcptr u; + mpz_srcptr v; +#endif +#else /* BERKELEY_MP */ +int +#if __STDC__ +mcmp (mpz_srcptr u, mpz_srcptr v) +#else +mcmp (u, v) + mpz_srcptr u; + mpz_srcptr v; +#endif +#endif /* BERKELEY_MP */ +{ + mp_size_t usize = u->_mp_size; + mp_size_t vsize = v->_mp_size; + mp_size_t size; + mp_srcptr up, vp; + int cmp; + + if (usize != vsize) + return usize - vsize; + + if (usize == 0) + return 0; + + size = ABS (usize); + + up = u->_mp_d; + vp = v->_mp_d; + + cmp = mpn_cmp (up, vp, size); + + if (cmp == 0) + return 0; + + if ((cmp < 0) == (usize < 0)) + return 1; + else + return -1; +} diff --git a/ghc/rts/gmp/mpz/cmp_si.c b/ghc/rts/gmp/mpz/cmp_si.c new file mode 100644 index 0000000..8063fd3 --- /dev/null +++ b/ghc/rts/gmp/mpz/cmp_si.c @@ -0,0 +1,66 @@ +/* mpz_cmp_si(u,v) -- Compare an integer U with a single-word int V. + Return positive, zero, or negative based on if U > V, U == V, or U < V. + +Copyright (C) 1991, 1993, 1994, 1995 Free Software Foundation, Inc. + +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" + +/* gmp.h defines a macro for mpz_cmp_si. */ +#undef mpz_cmp_si + +int +#if __STDC__ +mpz_cmp_si (mpz_srcptr u, signed long int v_digit) +#else +mpz_cmp_si (u, v_digit) + mpz_srcptr u; + signed long int v_digit; +#endif +{ + mp_size_t usize = u->_mp_size; + mp_size_t vsize; + mp_limb_t u_digit; + + vsize = 0; + if (v_digit > 0) + vsize = 1; + else if (v_digit < 0) + { + vsize = -1; + v_digit = -v_digit; + } + + if (usize != vsize) + return usize - vsize; + + if (usize == 0) + return 0; + + u_digit = u->_mp_d[0]; + + if (u_digit == v_digit) + return 0; + + if (u_digit > v_digit) + return usize; + else + return -usize; +} diff --git a/ghc/rts/gmp/mpz/cmp_ui.c b/ghc/rts/gmp/mpz/cmp_ui.c new file mode 100644 index 0000000..1a50b96 --- /dev/null +++ b/ghc/rts/gmp/mpz/cmp_ui.c @@ -0,0 +1,56 @@ +/* mpz_cmp_ui.c -- Compare a mpz_t a with an mp_limb_t b. Return positive, + zero, or negative based on if a > b, a == b, or a < b. + +Copyright (C) 1991, 1993, 1994, 1995 Free Software Foundation, Inc. + +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" + +/* gmp.h defines a macro for mpz_cmp_ui. */ +#undef mpz_cmp_ui + +int +#if __STDC__ +mpz_cmp_ui (mpz_srcptr u, unsigned long int v_digit) +#else +mpz_cmp_ui (u, v_digit) + mpz_srcptr u; + unsigned long int v_digit; +#endif +{ + mp_size_t usize = u->_mp_size; + + if (usize == 0) + return -(v_digit != 0); + + if (usize == 1) + { + mp_limb_t u_digit; + + u_digit = u->_mp_d[0]; + if (u_digit > v_digit) + return 1; + if (u_digit < v_digit) + return -1; + return 0; + } + + return (usize > 0) ? 1 : -1; +} diff --git a/ghc/rts/gmp/mpz/com.c b/ghc/rts/gmp/mpz/com.c new file mode 100644 index 0000000..559f2b6 --- /dev/null +++ b/ghc/rts/gmp/mpz/com.c @@ -0,0 +1,93 @@ +/* mpz_com(mpz_ptr dst, mpz_ptr src) -- Assign the bit-complemented value of + SRC to DST. + +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_com (mpz_ptr dst, mpz_srcptr src) +#else +mpz_com (dst, src) + mpz_ptr dst; + mpz_srcptr src; +#endif +{ + mp_size_t size = src->_mp_size; + mp_srcptr src_ptr; + mp_ptr dst_ptr; + + if (size >= 0) + { + /* As with infinite precision: one's complement, two's complement. + But this can be simplified using the identity -x = ~x + 1. + So we're going to compute (~~x) + 1 = x + 1! */ + + if (dst->_mp_alloc < size + 1) + _mpz_realloc (dst, size + 1); + + src_ptr = src->_mp_d; + dst_ptr = dst->_mp_d; + + if (size == 0) + { + /* Special case, as mpn_add wants the first arg's size >= the + second arg's size. */ + dst_ptr[0] = 1; + dst->_mp_size = -1; + return; + } + + { + mp_limb_t cy; + + cy = mpn_add_1 (dst_ptr, src_ptr, size, (mp_limb_t) 1); + if (cy) + { + dst_ptr[size] = cy; + size++; + } + } + + /* Store a negative size, to indicate ones-extension. */ + dst->_mp_size = -size; + } + else + { + /* As with infinite precision: two's complement, then one's complement. + But that can be simplified using the identity -x = ~(x - 1). + So we're going to compute ~~(x - 1) = x - 1! */ + size = -size; + + if (dst->_mp_alloc < size) + _mpz_realloc (dst, size); + + src_ptr = src->_mp_d; + dst_ptr = dst->_mp_d; + + mpn_sub_1 (dst_ptr, src_ptr, size, (mp_limb_t) 1); + size -= dst_ptr[size - 1] == 0; + + /* Store a positive size, to indicate zero-extension. */ + dst->_mp_size = size; + } +} diff --git a/ghc/rts/gmp/mpz/divexact.c b/ghc/rts/gmp/mpz/divexact.c new file mode 100644 index 0000000..b897448 --- /dev/null +++ b/ghc/rts/gmp/mpz/divexact.c @@ -0,0 +1,112 @@ +/* mpz_divexact -- finds quotient when known that quot * den == num && den != 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. */ + +/* Ken Weber (kweber@mat.ufrgs.br, kweber@mcs.kent.edu) + + Funding for this work has been partially provided by Conselho Nacional + de Desenvolvimento Cienti'fico e Tecnolo'gico (CNPq) do Brazil, Grant + 301314194-2, and was done while I was a visiting reseacher in the Instituto + de Matema'tica at Universidade Federal do Rio Grande do Sul (UFRGS). + + References: + T. Jebelean, An algorithm for exact division, Journal of Symbolic + Computation, v. 15, 1993, pp. 169-180. */ + +#include "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" + +void +#if __STDC__ +mpz_divexact (mpz_ptr quot, mpz_srcptr num, mpz_srcptr den) +#else +mpz_divexact (quot, num, den) + mpz_ptr quot; + mpz_srcptr num; + mpz_srcptr den; +#endif +{ + mp_ptr qp, tp; + mp_size_t qsize, tsize; + + mp_srcptr np = num->_mp_d; + mp_srcptr dp = den->_mp_d; + mp_size_t nsize = ABS (num->_mp_size); + mp_size_t dsize = ABS (den->_mp_size); + TMP_DECL (marker); + + /* Generate divide-by-zero error if dsize == 0. */ + if (dsize == 0) + { + quot->_mp_size = 1 / dsize; + return; + } + + if (nsize == 0) + { + quot->_mp_size = 0; + return; + } + + qsize = nsize - dsize + 1; + if (quot->_mp_alloc < qsize) + _mpz_realloc (quot, qsize); + qp = quot->_mp_d; + + TMP_MARK (marker); + + /* QUOT <-- NUM/2^r, T <-- DEN/2^r where = r number of twos in DEN. */ + while (dp[0] == 0) + np += 1, nsize -= 1, dp += 1, dsize -= 1; + tsize = MIN (qsize, dsize); + if (dp[0] & 1) + { + if (qp != dp) + MPN_COPY (qp, np, qsize); + if (qp == dp) /* QUOT and DEN overlap. */ + { + tp = (mp_ptr) TMP_ALLOC (sizeof (mp_limb_t) * tsize); + MPN_COPY (tp, dp, tsize); + } + else + tp = (mp_ptr) dp; + } + else + { + unsigned long int r; + tp = (mp_ptr) TMP_ALLOC (sizeof (mp_limb_t) * tsize); + count_trailing_zeros (r, dp[0]); + mpn_rshift (tp, dp, tsize, r); + if (dsize > tsize) + tp[tsize-1] |= dp[tsize] << (BITS_PER_MP_LIMB - r); + mpn_rshift (qp, np, qsize, r); + if (nsize > qsize) + qp[qsize-1] |= np[qsize] << (BITS_PER_MP_LIMB - r); + } + + /* Now QUOT <-- QUOT/T. */ + mpn_bdivmod (qp, qp, qsize, tp, tsize, qsize * BITS_PER_MP_LIMB); + MPN_NORMALIZE (qp, qsize); + + quot->_mp_size = (num->_mp_size < 0) == (den->_mp_size < 0) ? qsize : -qsize; + + TMP_FREE (marker); +} diff --git a/ghc/rts/gmp/mpz/fac_ui.c b/ghc/rts/gmp/mpz/fac_ui.c new file mode 100644 index 0000000..a170060 --- /dev/null +++ b/ghc/rts/gmp/mpz/fac_ui.c @@ -0,0 +1,157 @@ +/* mpz_fac_ui(result, n) -- Set RESULT to N!. + +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. */ + +#ifdef DBG +#include +#endif + +#include "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" + +void +#if __STDC__ +mpz_fac_ui (mpz_ptr result, unsigned long int n) +#else +mpz_fac_ui (result, n) + mpz_ptr result; + unsigned long int n; +#endif +{ +#if SIMPLE_FAC + + /* Be silly. Just multiply the numbers in ascending order. O(n**2). */ + + unsigned long int k; + + mpz_set_ui (result, 1L); + + for (k = 2; k <= n; k++) + mpz_mul_ui (result, result, k); +#else + + /* Be smarter. Multiply groups of numbers in ascending order until the + product doesn't fit in a limb. Multiply these partial product in a + balanced binary tree fashion, to make the operand have as equal sizes + as possible. When the operands have about the same size, mpn_mul + becomes faster. */ + + unsigned long int p, k; + mp_limb_t p1, p0; + + /* Stack of partial products, used to make the computation balanced + (i.e. make the sizes of the multiplication operands equal). The + topmost position of MP_STACK will contain a one-limb partial product, + the second topmost will contain a two-limb partial product, and so + on. MP_STACK[0] will contain a partial product with 2**t limbs. + To compute n! MP_STACK needs to be less than + log(n)**2/log(BITS_PER_MP_LIMB), so 30 is surely enough. */ +#define MP_STACK_SIZE 30 + mpz_t mp_stack[MP_STACK_SIZE]; + + /* TOP is an index into MP_STACK, giving the topmost element. + TOP_LIMIT_SO_FAR is the largets value it has taken so far. */ + int top, top_limit_so_far; + + /* Count of the total number of limbs put on MP_STACK so far. This + variable plays an essential role in making the compututation balanced. + See below. */ + unsigned int tree_cnt; + + top = top_limit_so_far = -1; + tree_cnt = 0; + p = 1; + for (k = 2; k <= n; k++) + { + /* Multiply the partial product in P with K. */ + umul_ppmm (p1, p0, (mp_limb_t) p, (mp_limb_t) k); + + /* Did we get overflow into the high limb, i.e. is the partial + product now more than one limb? */ + if (p1 != 0) + { + tree_cnt++; + + if (tree_cnt % 2 == 0) + { + mp_size_t i; + + /* TREE_CNT is even (i.e. we have generated an even number of + one-limb partial products), which means that we have a + single-limb product on the top of MP_STACK. */ + + mpz_mul_ui (mp_stack[top], mp_stack[top], p); + + /* If TREE_CNT is divisable by 4, 8,..., we have two + similar-sized partial products with 2, 4,... limbs at + the topmost two positions of MP_STACK. Multiply them + to form a new partial product with 4, 8,... limbs. */ + for (i = 4; (tree_cnt & (i - 1)) == 0; i <<= 1) + { + mpz_mul (mp_stack[top - 1], + mp_stack[top], mp_stack[top - 1]); + top--; + } + } + else + { + /* Put the single-limb partial product in P on the stack. + (The next time we get a single-limb product, we will + multiply the two together.) */ + top++; + if (top > top_limit_so_far) + { + if (top > MP_STACK_SIZE) + abort(); + /* The stack is now bigger than ever, initialize the top + element. */ + mpz_init_set_ui (mp_stack[top], p); + top_limit_so_far++; + } + else + mpz_set_ui (mp_stack[top], p); + } + + /* We ignored the last result from umul_ppmm. Put K in P as the + first component of the next single-limb partial product. */ + p = k; + } + else + /* We didn't get overflow in umul_ppmm. Put p0 in P and try + with one more value of K. */ + p = p0; /* bogus if long != mp_limb_t */ + } + + /* We have partial products in mp_stack[0..top], in descending order. + We also have a small partial product in p. + Their product is the final result. */ + if (top < 0) + mpz_set_ui (result, p); + else + mpz_mul_ui (result, mp_stack[top--], p); + while (top >= 0) + mpz_mul (result, result, mp_stack[top--]); + + /* Free the storage allocated for MP_STACK. */ + for (top = top_limit_so_far; top >= 0; top--) + mpz_clear (mp_stack[top]); +#endif +} diff --git a/ghc/rts/gmp/mpz/fdiv_q.c b/ghc/rts/gmp/mpz/fdiv_q.c new file mode 100644 index 0000000..3da943a --- /dev/null +++ b/ghc/rts/gmp/mpz/fdiv_q.c @@ -0,0 +1,51 @@ +/* mpz_fdiv_q -- Division rounding the quotient towards -infinity. + The remainder gets the same sign as the denominator. + +Copyright (C) 1994, 1995, 1996 Free Software Foundation, Inc. + +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_q (mpz_ptr quot, mpz_srcptr dividend, mpz_srcptr divisor) +#else +mpz_fdiv_q (quot, dividend, divisor) + mpz_ptr quot; + mpz_srcptr dividend; + mpz_srcptr divisor; +#endif +{ + mp_size_t dividend_size = dividend->_mp_size; + mp_size_t divisor_size = divisor->_mp_size; + mpz_t rem; + TMP_DECL (marker); + + TMP_MARK (marker); + + MPZ_TMP_INIT (rem, 1 + ABS (dividend_size)); + + mpz_tdiv_qr (quot, rem, dividend, divisor); + + if ((divisor_size ^ dividend_size) < 0 && rem->_mp_size != 0) + mpz_sub_ui (quot, quot, 1L); + + TMP_FREE (marker); +} diff --git a/ghc/rts/gmp/mpz/fdiv_q_2exp.c b/ghc/rts/gmp/mpz/fdiv_q_2exp.c new file mode 100644 index 0000000..3f56baf --- /dev/null +++ b/ghc/rts/gmp/mpz/fdiv_q_2exp.c @@ -0,0 +1,94 @@ +/* mpz_fdiv_q_2exp -- Divide an integer by 2**CNT. Round the quotient + towards -infinity. + +Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc. + +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_q_2exp (mpz_ptr w, mpz_srcptr u, unsigned long int cnt) +#else +mpz_fdiv_q_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 wsize; + mp_size_t abs_usize = ABS (usize); + mp_size_t limb_cnt; + mp_ptr wp; + mp_limb_t round = 0; + + limb_cnt = cnt / BITS_PER_MP_LIMB; + wsize = abs_usize - limb_cnt; + if (wsize <= 0) + { + wp = w->_mp_d; + wsize = 0; + /* Set ROUND since we know we skip some non-zero words in this case. + Well, if U is zero, we don't, but then this will be taken care of + below, since rounding only really takes place for negative U. */ + round = 1; + wp[0] = 1; + w->_mp_size = -(usize < 0); + return; + } + else + { + mp_size_t i; + mp_ptr up; + + /* Make sure there is enough space. We make an extra limb + here to account for possible rounding at the end. */ + if (w->_mp_alloc < wsize + 1) + _mpz_realloc (w, wsize + 1); + + wp = w->_mp_d; + up = u->_mp_d; + + /* Set ROUND if we are about skip some non-zero limbs. */ + for (i = 0; i < limb_cnt && round == 0; i++) + round = up[i]; + + cnt %= BITS_PER_MP_LIMB; + if (cnt != 0) + { + round |= mpn_rshift (wp, up + limb_cnt, wsize, cnt); + wsize -= wp[wsize - 1] == 0; + } + else + { + MPN_COPY_INCR (wp, up + limb_cnt, wsize); + } + } + + if (usize < 0 && round != 0) + { + mp_limb_t cy; + cy = mpn_add_1 (wp, wp, wsize, 1); + wp[wsize] = cy; + wsize += cy; + } + w->_mp_size = usize >= 0 ? wsize : -wsize; +} diff --git a/ghc/rts/gmp/mpz/fdiv_q_ui.c b/ghc/rts/gmp/mpz/fdiv_q_ui.c new file mode 100644 index 0000000..3d6825d --- /dev/null +++ b/ghc/rts/gmp/mpz/fdiv_q_ui.c @@ -0,0 +1,62 @@ +/* mpz_fdiv_q_ui -- Division rounding the quotient towards -infinity. + The remainder gets the same sign as the denominator. + +Copyright (C) 1994, 1995, 1996 Free Software Foundation, Inc. + +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_q_ui (mpz_ptr quot, mpz_srcptr dividend, unsigned long int divisor) +#else +mpz_fdiv_q_ui (quot, dividend, divisor) + mpz_ptr quot; + mpz_srcptr dividend; + unsigned long int divisor; +#endif +{ + mp_size_t dividend_size; + mp_size_t size; + mp_ptr quot_ptr; + mp_limb_t remainder_limb; + + dividend_size = dividend->_mp_size; + size = ABS (dividend_size); + + if (quot->_mp_alloc < size) + _mpz_realloc (quot, size); + + quot_ptr = quot->_mp_d; + + remainder_limb = mpn_divmod_1 (quot_ptr, dividend->_mp_d, size, + (mp_limb_t) divisor); + + if (remainder_limb != 0 && dividend_size < 0) + { + mpn_add_1 (quot_ptr, quot_ptr, size, (mp_limb_t) 1); + remainder_limb = divisor - remainder_limb; + } + + size -= size != 0 && quot_ptr[size - 1] == 0; + quot->_mp_size = dividend_size >= 0 ? size : -size; + + return remainder_limb; +} diff --git a/ghc/rts/gmp/mpz/fdiv_qr.c b/ghc/rts/gmp/mpz/fdiv_qr.c new file mode 100644 index 0000000..2abb16c --- /dev/null +++ b/ghc/rts/gmp/mpz/fdiv_qr.c @@ -0,0 +1,62 @@ +/* mpz_fdiv_qr -- Division rounding the quotient towards -infinity. + The remainder gets the same sign as the denominator. + +Copyright (C) 1994, 1995, 1996 Free Software Foundation, Inc. + +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_qr (mpz_ptr quot, mpz_ptr rem, mpz_srcptr dividend, mpz_srcptr divisor) +#else +mpz_fdiv_qr (quot, rem, dividend, divisor) + mpz_ptr quot; + 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 quotient and + remainder have been preliminary calculated. We have to copy it to + temporary space if it's the same variable as either QUOT or REM. */ + if (quot == divisor || rem == divisor) + { + MPZ_TMP_INIT (temp_divisor, ABS (divisor_size)); + mpz_set (temp_divisor, divisor); + divisor = temp_divisor; + } + + mpz_tdiv_qr (quot, rem, dividend, divisor); + + if ((divisor_size ^ dividend->_mp_size) < 0 && rem->_mp_size != 0) + { + mpz_sub_ui (quot, quot, 1L); + mpz_add (rem, rem, divisor); + } + + TMP_FREE (marker); +} diff --git a/ghc/rts/gmp/mpz/fdiv_qr_ui.c b/ghc/rts/gmp/mpz/fdiv_qr_ui.c new file mode 100644 index 0000000..a22b702 --- /dev/null +++ b/ghc/rts/gmp/mpz/fdiv_qr_ui.c @@ -0,0 +1,66 @@ +/* mpz_fdiv_qr_ui -- Division rounding the quotient towards -infinity. + The remainder gets the same sign as the denominator. + +Copyright (C) 1994, 1995, 1996 Free Software Foundation, Inc. + +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_qr_ui (mpz_ptr quot, mpz_ptr rem, mpz_srcptr dividend, unsigned long int divisor) +#else +mpz_fdiv_qr_ui (quot, rem, dividend, divisor) + mpz_ptr quot; + mpz_ptr rem; + mpz_srcptr dividend; + unsigned long int divisor; +#endif +{ + mp_size_t dividend_size; + mp_size_t size; + mp_ptr quot_ptr; + mp_limb_t remainder_limb; + + dividend_size = dividend->_mp_size; + size = ABS (dividend_size); + + if (quot->_mp_alloc < size) + _mpz_realloc (quot, size); + + quot_ptr = quot->_mp_d; + + remainder_limb = mpn_divmod_1 (quot_ptr, dividend->_mp_d, size, + (mp_limb_t) divisor); + + if (remainder_limb != 0 && dividend_size < 0) + { + mpn_add_1 (quot_ptr, quot_ptr, size, (mp_limb_t) 1); + remainder_limb = divisor - remainder_limb; + } + + size -= size != 0 && quot_ptr[size - 1] == 0; + quot->_mp_size = dividend_size >= 0 ? size : -size; + + rem->_mp_d[0] = remainder_limb; + rem->_mp_size = remainder_limb != 0; + + return remainder_limb; +} diff --git a/ghc/rts/gmp/mpz/fdiv_r.c b/ghc/rts/gmp/mpz/fdiv_r.c new file mode 100644 index 0000000..14e045b --- /dev/null +++ b/ghc/rts/gmp/mpz/fdiv_r.c @@ -0,0 +1,58 @@ +/* mpz_fdiv_r -- 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" + +void +#if __STDC__ +mpz_fdiv_r (mpz_ptr rem, mpz_srcptr dividend, mpz_srcptr divisor) +#else +mpz_fdiv_r (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 ((divisor_size ^ dividend->_mp_size) < 0 && rem->_mp_size != 0) + mpz_add (rem, rem, divisor); + + TMP_FREE (marker); +}