diff options
Diffstat (limited to 'mpfr/mul_ui.c')
-rw-r--r-- | mpfr/mul_ui.c | 85 |
1 files changed, 44 insertions, 41 deletions
diff --git a/mpfr/mul_ui.c b/mpfr/mul_ui.c index 8c916273d..72f9decef 100644 --- a/mpfr/mul_ui.c +++ b/mpfr/mul_ui.c @@ -1,6 +1,6 @@ /* mpfr_mul_ui -- multiply a floating-point number by a machine integer -Copyright (C) 1999, 2001, 2002 Free Software Foundation, Inc. +Copyright (C) 1999-2002 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -19,7 +19,6 @@ 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. */ -#include <stdio.h> #include "gmp.h" #include "gmp-impl.h" #include "longlong.h" @@ -29,90 +28,94 @@ MA 02111-1307, USA. */ int mpfr_mul_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int u, mp_rnd_t rnd_mode) { - mp_limb_t *my, *old_my; - unsigned long xn, yn, cnt, c; - int inexact; + mp_limb_t *yp, *old_yp; + mp_size_t xn, yn; + int cnt, c, inexact; TMP_DECL(marker); if (MPFR_IS_NAN(x)) { - MPFR_CLEAR_FLAGS(y); MPFR_SET_NAN(y); - return 1; /* a NaN is always inexact */ + MPFR_RET_NAN; } if (MPFR_IS_INF(x)) { - MPFR_CLEAR_FLAGS(y); - if (u) - { - MPFR_SET_INF(y); - if (MPFR_SIGN(y) != MPFR_SIGN(x)) - MPFR_CHANGE_SIGN(y); - return 0; /* infinity is exact */ + if (u != 0) + { + MPFR_CLEAR_FLAGS(y); + MPFR_SET_INF(y); + MPFR_SET_SAME_SIGN(y, x); + MPFR_RET(0); /* infinity is exact */ } else /* 0 * infinity */ { MPFR_SET_NAN(y); - return 1; /* NaN is inexact */ + MPFR_RET_NAN; } } - if (MPFR_IS_ZERO(x) || !u) + MPFR_CLEAR_FLAGS(y); + + if (u == 0 || MPFR_IS_ZERO(x)) { - MPFR_CLEAR_FLAGS(y); MPFR_SET_ZERO(y); - return 0; /* zero is exact */ + MPFR_SET_SAME_SIGN(y, x); + MPFR_RET(0); /* zero is exact */ } - MPFR_CLEAR_FLAGS(y); - if (u == 1) return mpfr_set (y, x, rnd_mode); TMP_MARK(marker); - my = MPFR_MANT(y); + yp = MPFR_MANT(y); yn = (MPFR_PREC(y) - 1) / BITS_PER_MP_LIMB + 1; xn = (MPFR_PREC(x) - 1) / BITS_PER_MP_LIMB + 1; - old_my = my; + old_yp = yp; + MPFR_ASSERTN(xn < MP_SIZE_T_MAX); if (yn < xn + 1) - my = (mp_ptr) TMP_ALLOC ((xn + 1) * BYTES_PER_MP_LIMB); + yp = (mp_ptr) TMP_ALLOC ((size_t) (xn + 1) * BYTES_PER_MP_LIMB); + + yp[xn] = mpn_mul_1 (yp, MPFR_MANT(x), xn, u); - my[xn] = mpn_mul_1 (my, MPFR_MANT(x), xn, u); - - /* x * u is stored in my[xn], ..., my[0] */ + /* x * u is stored in yp[xn], ..., yp[0] */ /* since the case u=1 was treated above, we have u >= 2, thus - my[xn] >= 1 since x was msb-normalized */ - MPFR_ASSERTN(my[xn] != 0); - cnt = 0; - if ((my[xn] & MP_LIMB_T_HIGHBIT) == 0) + yp[xn] >= 1 since x was msb-normalized */ + MPFR_ASSERTN(yp[xn] != 0); + if ((yp[xn] & MP_LIMB_T_HIGHBIT) == 0) { - count_leading_zeros(cnt, my[xn]); - mpn_lshift (my, my, xn + 1, cnt); + count_leading_zeros(cnt, yp[xn]); + mpn_lshift (yp, yp, xn + 1, cnt); + } + else + { + cnt = 0; } - /* now my[xn], ..., my[0] is msb-normalized too, and has at most + /* now yp[xn], ..., yp[0] is msb-normalized too, and has at most PREC(x) + (BITS_PER_MP_LIMB - cnt) non-zero bits */ - c = mpfr_round_raw (old_my, my, (xn + 1) * BITS_PER_MP_LIMB, + c = mpfr_round_raw (old_yp, yp, (mp_prec_t) (xn + 1) * BITS_PER_MP_LIMB, (MPFR_SIGN(x) < 0), MPFR_PREC(y), rnd_mode, &inexact); - MPFR_EXP(y) = MPFR_EXP(x) + BITS_PER_MP_LIMB - cnt; + cnt = BITS_PER_MP_LIMB - cnt; if (c) /* rounded result is 1.0000000000000000... */ { - old_my[yn-1] = MP_LIMB_T_HIGHBIT; - MPFR_EXP(y) ++; + old_yp[yn-1] = MP_LIMB_T_HIGHBIT; + cnt++; } - /* set sign */ - if (MPFR_SIGN(y) * MPFR_SIGN(x) < 0) - MPFR_CHANGE_SIGN(y); - TMP_FREE(marker); + if (__mpfr_emax < MPFR_EMAX_MIN + cnt || MPFR_EXP(x) > __mpfr_emax - cnt) + return mpfr_set_overflow(y, rnd_mode, MPFR_SIGN(x)); + + MPFR_EXP(y) = MPFR_EXP(x) + cnt; + MPFR_SET_SAME_SIGN(y, x); + return inexact; } |