=================================================================== RCS file: /home/cvs/OpenXM_contrib/gmp/mpfr/Attic/div_ui.c,v retrieving revision 1.1.1.1 retrieving revision 1.1.1.2 diff -u -p -r1.1.1.1 -r1.1.1.2 --- OpenXM_contrib/gmp/mpfr/Attic/div_ui.c 2000/09/09 14:12:19 1.1.1.1 +++ OpenXM_contrib/gmp/mpfr/Attic/div_ui.c 2003/08/25 16:06:06 1.1.1.2 @@ -1,20 +1,20 @@ /* mpfr_div_ui -- divide a floating-point number by a machine integer -Copyright (C) 1999 PolKA project, Inria Lorraine and Loria +Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc. This file is part of the MPFR Library. The MPFR Library is free software; you can redistribute it and/or modify -it under the terms of the GNU Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 of the License, or (at your option) any later version. The MPFR Library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY -or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. -You should have received a copy of the GNU Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the MPFR Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -22,113 +22,163 @@ MA 02111-1307, USA. */ #include #include "gmp.h" #include "gmp-impl.h" -#include "longlong.h" +#include "longlong.h" #include "mpfr.h" +#include "mpfr-impl.h" -/* #define DEBUG */ - /* returns 0 if result exact, non-zero otherwise */ int -#ifdef __STDC__ -mpfr_div_ui(mpfr_ptr y, mpfr_srcptr x, unsigned long u, unsigned char rnd_mode) -#else -mpfr_div_ui(y, x, u, rnd_mode) - mpfr_ptr y; - mpfr_srcptr x; - unsigned long u; - unsigned char rnd_mode; -#endif +mpfr_div_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int u, mp_rnd_t rnd_mode) { - int xn, yn, dif, sh, i; mp_limb_t *xp, *yp, *tmp, c, d; + long int xn, yn, dif, sh, i; + mp_limb_t *xp, *yp, *tmp, c, d; + int inexact, middle = 1; TMP_DECL(marker); - if (FLAG_NAN(x)) { SET_NAN(y); return 1; } - if (u==0) { printf("infinity\n"); return 1; } + if (MPFR_IS_NAN(x)) + { + MPFR_SET_NAN(y); + MPFR_RET_NAN; + } + MPFR_CLEAR_NAN(y); /* clear NaN flag */ + + if (MPFR_IS_INF(x)) + { + MPFR_SET_INF(y); + MPFR_SET_SAME_SIGN(y, x); + MPFR_RET(0); + } + + if (u == 0) + { + if (MPFR_IS_ZERO(x)) /* 0/0 is NaN */ + { + MPFR_SET_NAN(y); + MPFR_RET_NAN; + } + else /* x/0 is Inf */ + { + MPFR_SET_INF(y); + MPFR_SET_SAME_SIGN(y, x); + MPFR_RET(0); + } + } + + MPFR_CLEAR_INF(y); + MPFR_SET_SAME_SIGN(y, x); + + if (MPFR_IS_ZERO(x)) + { + MPFR_SET_ZERO(y); + MPFR_RET(0); + } + TMP_MARK(marker); - xn = (PREC(x)-1)/BITS_PER_MP_LIMB + 1; - yn = (PREC(y)-1)/BITS_PER_MP_LIMB + 1; + xn = (MPFR_PREC(x) - 1)/BITS_PER_MP_LIMB + 1; + yn = (MPFR_PREC(y) - 1)/BITS_PER_MP_LIMB + 1; - xp = MANT(x); - yp = MANT(y); - EXP(y) = EXP(x); - if (SIGN(x)!=SIGN(y)) CHANGE_SIGN(y); + xp = MPFR_MANT(x); + yp = MPFR_MANT(y); + MPFR_EXP(y) = MPFR_EXP(x); - dif = yn+1-xn; -#ifdef DEBUG - printf("dif=%d u=%lu xn=%d\n",dif,u,xn); - printf("x="); mpfr_print_raw(x); putchar('\n'); -#endif + dif = yn + 1 - xn; /* we need to store yn+1 = xn + dif limbs of the quotient */ - if (ABSSIZE(y)>=yn+1) tmp=yp; - else tmp=TMP_ALLOC((yn+1)*BYTES_PER_MP_LIMB); + /* don't use tmp=yp since the mpn_lshift call below requires yp >= tmp+1 */ + tmp = TMP_ALLOC((yn + 1) * BYTES_PER_MP_LIMB); c = (mp_limb_t) u; - if (dif>=0) { - /* patch for bug in mpn_divrem_1 */ -#if (UDIV_NEEDS_NORMALIZATION==1) - count_leading_zeros(sh, c); - c <<= sh; - EXP(y) += sh; -#endif - c = mpn_divrem_1(tmp, dif, xp, xn, c); - } - else /* dif < 0 i.e. xn > yn */ - c = mpn_divrem_1(tmp, 0, xp-dif, yn, c); + if (dif >= 0) + c = mpn_divrem_1 (tmp, dif, xp, xn, c); /* used all the dividend */ + else /* dif < 0 i.e. xn > yn, don't use the (-dif) low limbs from x */ + c = mpn_divrem_1 (tmp, 0, xp - dif, yn + 1, c); - if (tmp[yn]==0) { tmp--; sh=0; EXP(y) -= mp_bits_per_limb; } + inexact = (c != 0); + if (rnd_mode == GMP_RNDN) + { + if (2 * c < (mp_limb_t) u) + middle = -1; + else if (2 * c > (mp_limb_t) u) + middle = 1; + else + middle = 0; /* exactly in the middle */ + } + for (i=0; ((inexact == 0) || (middle == 0)) && (i < -dif); i++) + if (xp[i]) + inexact = middle = 1; /* larger than middle */ + + if (tmp[yn] == 0) /* high limb is zero */ + { + tmp--; + sh = 0; + MPFR_EXP(y) -= BITS_PER_MP_LIMB; + } + + /* now we have yn limbs starting from tmp[1], with tmp[yn]<>0 */ + /* shift left to normalize */ count_leading_zeros(sh, tmp[yn]); - if (sh) { - mpn_lshift(yp, tmp+1, yn, sh); - yp[0] += tmp[0] >> (BITS_PER_MP_LIMB-sh); - EXP(y) -= sh; - } - else MPN_COPY(yp, tmp+1, yn); -#ifdef DEBUG -printf("y="); mpfr_print_raw(y); putchar('\n'); -#endif + if (sh) + { + mpn_lshift (yp, tmp + 1, yn, sh); + yp[0] += tmp[0] >> (BITS_PER_MP_LIMB - sh); + middle = middle || ((tmp[0] << sh) != 0); + inexact = inexact || ((tmp[0] << sh) != 0); + MPFR_EXP(y) -= sh; + } + else + MPN_COPY(yp, tmp + 1, yn); - sh = yn*BITS_PER_MP_LIMB - PREC(y); + sh = yn * BITS_PER_MP_LIMB - MPFR_PREC(y); /* it remains sh bits in less significant limb of y */ - d = *yp & (((mp_limb_t)1 << sh) - 1); + d = *yp & ((MP_LIMB_T_ONE << sh) - MP_LIMB_T_ONE); *yp ^= d; /* set to zero lowest sh bits */ TMP_FREE(marker); - if ((c | d)==0) { - for (i=0; i<-dif && xp[i]==0; i++); - if (i>=-dif) return 0; /* result is exact */ - } + if ((d == 0) && (inexact == 0)) + return 0; /* result is exact */ - switch (rnd_mode) { - case GMP_RNDZ: - return 1; /* result is inexact */ - case GMP_RNDU: - if (SIGN(y)>0) mpfr_add_one_ulp(y); - return 1; /* result is inexact */ - case GMP_RNDD: - if (SIGN(y)<0) mpfr_add_one_ulp(y); - return 1; /* result is inexact */ - case GMP_RNDN: - if (d < ((mp_limb_t)1 << (sh-1))) return 1; - else if (d > ((mp_limb_t)1 << (sh-1))) { - mpfr_add_one_ulp(y); - } - else { /* d = (mp_limb_t)1 << (sh-1) */ - if (c) mpfr_add_one_ulp(y); - else { - for (i=0; i<-dif && xp[i]==0; i++); - if (i<-dif) mpfr_add_one_ulp(y); - else { /* exactly in the middle */ - if (*yp & ((mp_limb_t)1 << sh)) mpfr_add_one_ulp(y); + switch (rnd_mode) + { + case GMP_RNDZ: + MPFR_RET(-MPFR_SIGN(x)); /* result is inexact */ + + case GMP_RNDU: + if (MPFR_SIGN(y) > 0) + mpfr_add_one_ulp (y, rnd_mode); + MPFR_RET(1); /* result is inexact */ + + case GMP_RNDD: + if (MPFR_SIGN(y) < 0) + mpfr_add_one_ulp (y, rnd_mode); + MPFR_RET(-1); /* result is inexact */ + + case GMP_RNDN: + if (sh && d < (MP_LIMB_T_ONE << (sh - 1))) + MPFR_RET(-MPFR_SIGN(x)); + else if (sh && d > (MP_LIMB_T_ONE << (sh - 1))) + { + mpfr_add_one_ulp (y, rnd_mode); + MPFR_RET(MPFR_SIGN(x)); } + else /* sh = 0 or d = 1 << (sh-1) */ + { + /* we are in the middle if: + (a) sh > 0 and inexact == 0 + (b) sh=0 and middle=1 + */ + if ((sh && inexact) || (!sh && (middle > 0)) || (*yp & (MP_LIMB_T_ONE << sh))) + { + mpfr_add_one_ulp (y, rnd_mode); + MPFR_RET(MPFR_SIGN(x)); + } + else + MPFR_RET(-MPFR_SIGN(x)); } } - return 1; - } - return 0; /* to prevent warning from gcc */ + MPFR_RET(inexact); /* should never go here */ }