diff options
Diffstat (limited to 'mpfr/mul_ui.c')
-rw-r--r-- | mpfr/mul_ui.c | 136 |
1 files changed, 58 insertions, 78 deletions
diff --git a/mpfr/mul_ui.c b/mpfr/mul_ui.c index e92a79155..c997c41d9 100644 --- a/mpfr/mul_ui.c +++ b/mpfr/mul_ui.c @@ -1,20 +1,20 @@ /* mpfr_mul_ui -- multiply a floating-point number by a machine integer -Copyright (C) 1999 Free Software Foundation. +Copyright (C) 1999, 2001 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. */ @@ -26,26 +26,19 @@ MA 02111-1307, USA. */ #include "mpfr.h" #include "mpfr-impl.h" -void -#if __STDC__ -mpfr_mul_ui(mpfr_ptr y, mpfr_srcptr x, unsigned long int u, mp_rnd_t rnd_mode) -#else -mpfr_mul_ui(y, x, u, rnd_mode) - mpfr_ptr y; - mpfr_srcptr x; - unsigned long int u; - mp_rnd_t rnd_mode; -#endif +int +mpfr_mul_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int u, mp_rnd_t rnd_mode) { - mp_limb_t carry, *my, *old_my, *my2; unsigned long c; - unsigned long xsize, ysize, cnt, dif, ex, sh; + mp_limb_t *my, *old_my; + unsigned long xn, yn, cnt, c; + int inexact; TMP_DECL(marker); if (MPFR_IS_NAN(x)) { MPFR_CLEAR_FLAGS(y); MPFR_SET_NAN(y); - return; + return 1; /* a NaN is always inexact */ } if (MPFR_IS_INF(x)) @@ -54,81 +47,68 @@ mpfr_mul_ui(y, x, u, rnd_mode) if (u) { MPFR_SET_INF(y); - if (MPFR_SIGN(y) != MPFR_SIGN(x)) { MPFR_CHANGE_SIGN(y); } - return; + if (MPFR_SIGN(y) != MPFR_SIGN(x)) + MPFR_CHANGE_SIGN(y); + return 0; /* infinity is exact */ + } + else /* 0 * infinity */ + { + MPFR_SET_NAN(y); + return 1; /* NaN is inexact */ } - else { MPFR_SET_NAN(y); return; } } - MPFR_CLEAR_FLAGS(y); + if (MPFR_IS_ZERO(x) || !u) + { + MPFR_CLEAR_FLAGS(y); + MPFR_SET_ZERO(y); + return 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); ex = MPFR_EXP(x); - ysize = (MPFR_PREC(y)-1)/BITS_PER_MP_LIMB + 1; - xsize = (MPFR_PREC(x)-1)/BITS_PER_MP_LIMB + 1; + my = 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_my = my; - if (ysize < xsize) { - my = (mp_ptr) TMP_ALLOC (xsize * BYTES_PER_MP_LIMB); - dif=0; - } - else dif=ysize-xsize; + if (yn < xn + 1) + my = (mp_ptr) TMP_ALLOC ((xn + 1) * BYTES_PER_MP_LIMB); - carry = mpn_mul_1(my+dif, MPFR_MANT(x), xsize, u); - MPN_ZERO(my, dif); + my[xn] = mpn_mul_1 (my, MPFR_MANT(x), xn, u); + + /* x * u is stored in my[xn], ..., my[0] */ - /* WARNING: count_leading_zeros is undefined for carry=0 */ - if (carry) count_leading_zeros(cnt, carry); - else cnt=BITS_PER_MP_LIMB; + /* 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); + count_leading_zeros(cnt, my[xn]); + mpn_lshift (my, my, xn + 1, cnt); - /* Warning: the number of limbs used by x and the lower part - of y may differ */ - sh = (MPFR_PREC(x)+BITS_PER_MP_LIMB-1)/BITS_PER_MP_LIMB - - (MPFR_PREC(y)+cnt-1)/BITS_PER_MP_LIMB; + /* now my[xn], ..., my[0] is msb-normalized too, and has at most + PREC(x) + (BITS_PER_MP_LIMB - cnt) non-zero bits */ - /* Warning: if all significant bits are in the carry, one has to - be careful */ - if (cnt + MPFR_PREC(y) < BITS_PER_MP_LIMB) - { - /* Quick 'n dirty */ - - if (xsize > ysize) { - my2 = (mp_ptr) TMP_ALLOC ((xsize + 1) * BYTES_PER_MP_LIMB); - my2[xsize] = mpn_lshift(my2, my, xsize, cnt) - | (carry << (cnt ? BITS_PER_MP_LIMB - cnt : 0)); - } - else { - my2 = (mp_ptr) TMP_ALLOC ((ysize + 1) * BYTES_PER_MP_LIMB); - my2[ysize] = mpn_lshift(my2, my, ysize, cnt) - | (carry << (cnt ? BITS_PER_MP_LIMB - cnt : 0)); - } - - my = my2; ex += BITS_PER_MP_LIMB - cnt; - carry = 0; cnt = BITS_PER_MP_LIMB; - } + c = mpfr_round_raw (old_my, my, (xn + 1) * BITS_PER_MP_LIMB, + (MPFR_SIGN(x) < 0), MPFR_PREC(y), rnd_mode, &inexact); - c = mpfr_round_raw(my+sh, my, MPFR_PREC(x), (MPFR_SIGN(x)<0), - MPFR_PREC(y)-BITS_PER_MP_LIMB+cnt, rnd_mode); + MPFR_EXP(y) = MPFR_EXP(x) + BITS_PER_MP_LIMB - cnt; - /* If cnt = 1111111111111 and c = 1 we shall get depressed */ - if (c && (carry == (((mp_limb_t)1) << (cnt ? BITS_PER_MP_LIMB - cnt : 0)) - - 1)) - { - cnt--; - mpn_rshift(my, my, ysize, BITS_PER_MP_LIMB - cnt); - my[ysize - 1] |= ((mp_limb_t) 1) << (BITS_PER_MP_LIMB - 1); - } - else + if (c) /* rounded result is 1.0000000000000000... */ { - /* Warning: mpn_rshift is undefined for shift=0 */ - if (cnt!=BITS_PER_MP_LIMB) - mpn_rshift(my, my, ysize, BITS_PER_MP_LIMB - cnt); - my[ysize - 1] |= (carry << cnt); + old_my[yn-1] = MP_LIMB_T_HIGHBIT; + MPFR_EXP(y) ++; } - MPFR_EXP(y) = ex + BITS_PER_MP_LIMB - cnt; - if (ysize < xsize) MPN_COPY(old_my, my, ysize); + /* set sign */ - if (MPFR_SIGN(y) * MPFR_SIGN(x) < 0) MPFR_CHANGE_SIGN(y); + if (MPFR_SIGN(y) * MPFR_SIGN(x) < 0) + MPFR_CHANGE_SIGN(y); + TMP_FREE(marker); + + return inexact; } |