version 1.1.1.1, 2000/09/09 14:12:19 |
version 1.1.1.2, 2003/08/25 16:06:06 |
|
|
/* mpfr_div_ui -- divide a floating-point number by a machine integer |
/* 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. |
This file is part of the MPFR Library. |
|
|
The MPFR Library is free software; you can redistribute it and/or modify |
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 |
it under the terms of the GNU Lesser General Public License as published by |
the Free Software Foundation; either version 2 of the License, or (at your |
the Free Software Foundation; either version 2.1 of the License, or (at your |
option) any later version. |
option) any later version. |
|
|
The MPFR Library is distributed in the hope that it will be useful, but |
The MPFR Library is distributed in the hope that it will be useful, but |
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |
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. |
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 |
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, |
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
MA 02111-1307, USA. */ |
MA 02111-1307, USA. */ |
Line 22 MA 02111-1307, USA. */ |
|
Line 22 MA 02111-1307, USA. */ |
|
#include <stdio.h> |
#include <stdio.h> |
#include "gmp.h" |
#include "gmp.h" |
#include "gmp-impl.h" |
#include "gmp-impl.h" |
#include "longlong.h" |
#include "longlong.h" |
#include "mpfr.h" |
#include "mpfr.h" |
|
#include "mpfr-impl.h" |
|
|
/* #define DEBUG */ |
|
|
|
/* returns 0 if result exact, non-zero otherwise */ |
/* returns 0 if result exact, non-zero otherwise */ |
int |
int |
#ifdef __STDC__ |
mpfr_div_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int u, mp_rnd_t rnd_mode) |
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 |
|
{ |
{ |
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); |
TMP_DECL(marker); |
|
|
if (FLAG_NAN(x)) { SET_NAN(y); return 1; } |
if (MPFR_IS_NAN(x)) |
if (u==0) { printf("infinity\n"); return 1; } |
{ |
|
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); |
TMP_MARK(marker); |
xn = (PREC(x)-1)/BITS_PER_MP_LIMB + 1; |
xn = (MPFR_PREC(x) - 1)/BITS_PER_MP_LIMB + 1; |
yn = (PREC(y)-1)/BITS_PER_MP_LIMB + 1; |
yn = (MPFR_PREC(y) - 1)/BITS_PER_MP_LIMB + 1; |
|
|
xp = MANT(x); |
xp = MPFR_MANT(x); |
yp = MANT(y); |
yp = MPFR_MANT(y); |
EXP(y) = EXP(x); |
MPFR_EXP(y) = MPFR_EXP(x); |
if (SIGN(x)!=SIGN(y)) CHANGE_SIGN(y); |
|
|
|
dif = yn+1-xn; |
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 |
|
|
|
/* we need to store yn+1 = xn + dif limbs of the quotient */ |
/* we need to store yn+1 = xn + dif limbs of the quotient */ |
if (ABSSIZE(y)>=yn+1) tmp=yp; |
/* don't use tmp=yp since the mpn_lshift call below requires yp >= tmp+1 */ |
else tmp=TMP_ALLOC((yn+1)*BYTES_PER_MP_LIMB); |
tmp = TMP_ALLOC((yn + 1) * BYTES_PER_MP_LIMB); |
|
|
c = (mp_limb_t) u; |
c = (mp_limb_t) u; |
if (dif>=0) { |
if (dif >= 0) |
/* patch for bug in mpn_divrem_1 */ |
c = mpn_divrem_1 (tmp, dif, xp, xn, c); /* used all the dividend */ |
#if (UDIV_NEEDS_NORMALIZATION==1) |
else /* dif < 0 i.e. xn > yn, don't use the (-dif) low limbs from x */ |
count_leading_zeros(sh, c); |
c = mpn_divrem_1 (tmp, 0, xp - dif, yn + 1, 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 (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 */ |
/* shift left to normalize */ |
count_leading_zeros(sh, tmp[yn]); |
count_leading_zeros(sh, tmp[yn]); |
if (sh) { |
if (sh) |
mpn_lshift(yp, tmp+1, yn, sh); |
{ |
yp[0] += tmp[0] >> (BITS_PER_MP_LIMB-sh); |
mpn_lshift (yp, tmp + 1, yn, sh); |
EXP(y) -= sh; |
yp[0] += tmp[0] >> (BITS_PER_MP_LIMB - sh); |
} |
middle = middle || ((tmp[0] << sh) != 0); |
else MPN_COPY(yp, tmp+1, yn); |
inexact = inexact || ((tmp[0] << sh) != 0); |
#ifdef DEBUG |
MPFR_EXP(y) -= sh; |
printf("y="); mpfr_print_raw(y); putchar('\n'); |
} |
#endif |
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 */ |
/* 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 */ |
*yp ^= d; /* set to zero lowest sh bits */ |
|
|
TMP_FREE(marker); |
TMP_FREE(marker); |
if ((c | d)==0) { |
if ((d == 0) && (inexact == 0)) |
for (i=0; i<-dif && xp[i]==0; i++); |
return 0; /* result is exact */ |
if (i>=-dif) return 0; /* result is exact */ |
|
} |
|
|
|
switch (rnd_mode) { |
switch (rnd_mode) |
case GMP_RNDZ: |
{ |
return 1; /* result is inexact */ |
case GMP_RNDZ: |
case GMP_RNDU: |
MPFR_RET(-MPFR_SIGN(x)); /* result is inexact */ |
if (SIGN(y)>0) mpfr_add_one_ulp(y); |
|
return 1; /* result is inexact */ |
case GMP_RNDU: |
case GMP_RNDD: |
if (MPFR_SIGN(y) > 0) |
if (SIGN(y)<0) mpfr_add_one_ulp(y); |
mpfr_add_one_ulp (y, rnd_mode); |
return 1; /* result is inexact */ |
MPFR_RET(1); /* result is inexact */ |
case GMP_RNDN: |
|
if (d < ((mp_limb_t)1 << (sh-1))) return 1; |
case GMP_RNDD: |
else if (d > ((mp_limb_t)1 << (sh-1))) { |
if (MPFR_SIGN(y) < 0) |
mpfr_add_one_ulp(y); |
mpfr_add_one_ulp (y, rnd_mode); |
} |
MPFR_RET(-1); /* result is inexact */ |
else { /* d = (mp_limb_t)1 << (sh-1) */ |
|
if (c) mpfr_add_one_ulp(y); |
case GMP_RNDN: |
else { |
if (sh && d < (MP_LIMB_T_ONE << (sh - 1))) |
for (i=0; i<-dif && xp[i]==0; i++); |
MPFR_RET(-MPFR_SIGN(x)); |
if (i<-dif) mpfr_add_one_ulp(y); |
else if (sh && d > (MP_LIMB_T_ONE << (sh - 1))) |
else { /* exactly in the middle */ |
{ |
if (*yp & ((mp_limb_t)1 << sh)) mpfr_add_one_ulp(y); |
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; |
MPFR_RET(inexact); /* should never go here */ |
} |
|
return 0; /* to prevent warning from gcc */ |
|
} |
} |
|
|
|
|