Annotation of OpenXM_contrib/gmp/mpfr/mul.c, Revision 1.1.1.2
1.1 maekawa 1: /* mpfr_mul -- multiply two floating-point numbers
2:
1.1.1.2 ! ohara 3: Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc.
1.1 maekawa 4:
5: This file is part of the MPFR Library.
6:
7: The MPFR Library is free software; you can redistribute it and/or modify
1.1.1.2 ! ohara 8: it under the terms of the GNU Lesser General Public License as published by
! 9: the Free Software Foundation; either version 2.1 of the License, or (at your
1.1 maekawa 10: option) any later version.
11:
12: The MPFR Library is distributed in the hope that it will be useful, but
13: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1.1.1.2 ! ohara 14: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
1.1 maekawa 15: License for more details.
16:
1.1.1.2 ! ohara 17: You should have received a copy of the GNU Lesser General Public License
1.1 maekawa 18: along with the MPFR Library; see the file COPYING.LIB. If not, write to
19: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
20: MA 02111-1307, USA. */
21:
22: #include "gmp.h"
23: #include "gmp-impl.h"
24: #include "mpfr.h"
1.1.1.2 ! ohara 25: #include "mpfr-impl.h"
1.1 maekawa 26:
1.1.1.2 ! ohara 27: int
! 28: mpfr_mul (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode)
1.1 maekawa 29: {
1.1.1.2 ! ohara 30: int sign_product, cc, inexact, ec, em = 0;
! 31: mp_exp_t bx, cx;
! 32: mp_limb_t *ap, *bp, *cp, *tmp;
! 33: mp_limb_t b1;
! 34: mp_prec_t aq, bq, cq;
! 35: mp_size_t an, bn, cn, tn, k;
! 36: TMP_DECL(marker);
1.1 maekawa 37:
38: /* deal with NaN and zero */
1.1.1.2 ! ohara 39: if (MPFR_IS_NAN(b) || MPFR_IS_NAN(c))
! 40: {
! 41: MPFR_SET_NAN(a);
! 42: MPFR_RET_NAN;
! 43: }
! 44:
! 45: MPFR_CLEAR_NAN(a);
! 46:
! 47: sign_product = MPFR_SIGN(b) * MPFR_SIGN(c);
! 48:
! 49: if (MPFR_IS_INF(b))
! 50: {
! 51: if (MPFR_IS_INF(c) || MPFR_NOTZERO(c))
! 52: {
! 53: if (MPFR_SIGN(a) != sign_product)
! 54: MPFR_CHANGE_SIGN(a);
! 55: MPFR_SET_INF(a);
! 56: MPFR_RET(0); /* exact */
! 57: }
! 58: else
! 59: {
! 60: MPFR_SET_NAN(a);
! 61: MPFR_RET_NAN;
! 62: }
! 63: }
! 64: else if (MPFR_IS_INF(c))
! 65: {
! 66: if (MPFR_NOTZERO(b))
! 67: {
! 68: if (MPFR_SIGN(a) != sign_product)
! 69: MPFR_CHANGE_SIGN(a);
! 70: MPFR_SET_INF(a);
! 71: MPFR_RET(0); /* exact */
! 72: }
! 73: else
! 74: {
! 75: MPFR_SET_NAN(a);
! 76: MPFR_RET_NAN;
! 77: }
! 78: }
! 79:
! 80: MPFR_ASSERTN(MPFR_IS_FP(b) && MPFR_IS_FP(c));
! 81: MPFR_CLEAR_INF(a); /* clear Inf flag */
! 82:
! 83: if (MPFR_IS_ZERO(b) || MPFR_IS_ZERO(c))
! 84: {
! 85: if (MPFR_SIGN(a) != sign_product)
! 86: MPFR_CHANGE_SIGN(a);
! 87: MPFR_SET_ZERO(a);
! 88: MPFR_RET(0); /* 0 * 0 is exact */
! 89: }
! 90:
! 91: bx = MPFR_EXP(b);
! 92: cx = MPFR_EXP(c);
! 93: /* Note: exponent of the result will be bx + cx + ec with ec in {-1,0,1} */
! 94: if (bx >= 0 && cx > 0)
! 95: { /* bx + cx > 0 */
! 96: if (__mpfr_emax < 0 ||
! 97: (mp_exp_unsigned_t) bx + cx > (mp_exp_unsigned_t) __mpfr_emax + 1)
! 98: return mpfr_set_overflow(a, rnd_mode, sign_product);
! 99:
! 100: if ((mp_exp_unsigned_t) bx + cx == (mp_exp_unsigned_t) __mpfr_emax + 1)
! 101: em = 1;
! 102: }
! 103: else if (bx <= 0 && cx < 0)
! 104: { /* bx + cx < 0 */
! 105: if (__mpfr_emin > 0 ||
! 106: (mp_exp_unsigned_t) bx + cx < (mp_exp_unsigned_t) __mpfr_emin - 1)
! 107: return mpfr_set_underflow(a, rnd_mode, sign_product);
! 108:
! 109: if ((mp_exp_unsigned_t) bx + cx == (mp_exp_unsigned_t) __mpfr_emin - 1)
! 110: em = -1;
! 111: }
! 112: else
! 113: { /* bx != 0 and cx doesn't have the same sign */
! 114: if ((bx + cx) - 1 > __mpfr_emax)
! 115: return mpfr_set_overflow(a, rnd_mode, sign_product);
! 116:
! 117: if ((bx + cx) - 1 == __mpfr_emax)
! 118: em = 1;
! 119:
! 120: if ((bx + cx) + 1 < __mpfr_emin)
! 121: return mpfr_set_underflow(a, rnd_mode, sign_product);
! 122:
! 123: if ((bx + cx) + 1 == __mpfr_emin)
! 124: em = -1;
! 125: }
! 126:
! 127: ap = MPFR_MANT(a);
! 128: bp = MPFR_MANT(b);
! 129: cp = MPFR_MANT(c);
! 130:
! 131: aq = MPFR_PREC(a);
! 132: bq = MPFR_PREC(b);
! 133: cq = MPFR_PREC(c);
! 134:
! 135: an = (aq-1)/BITS_PER_MP_LIMB + 1; /* number of significant limbs of a */
! 136: bn = (bq-1)/BITS_PER_MP_LIMB + 1; /* number of significant limbs of b */
! 137: cn = (cq-1)/BITS_PER_MP_LIMB + 1; /* number of significant limbs of c */
1.1 maekawa 138:
1.1.1.2 ! ohara 139: MPFR_ASSERTN((mp_size_unsigned_t) bn + cn <= MP_SIZE_T_MAX);
! 140: k = bn + cn; /* effective nb of limbs used by b*c (= tn or tn+1) below */
! 141:
! 142: MPFR_ASSERTN(bq + cq >= bq); /* no integer overflow */
! 143: tn = (bq + cq - 1) / BITS_PER_MP_LIMB + 1; /* <= k, thus no int overflow */
! 144:
! 145: MPFR_ASSERTN(k <= ((size_t) -1) / BYTES_PER_MP_LIMB);
1.1 maekawa 146: TMP_MARK(marker);
1.1.1.2 ! ohara 147: tmp = (mp_limb_t *) TMP_ALLOC((size_t) k * BYTES_PER_MP_LIMB);
1.1 maekawa 148:
149: /* multiplies two mantissa in temporary allocated space */
1.1.1.2 ! ohara 150: b1 = (bn >= cn) ? mpn_mul (tmp, bp, bn, cp, cn)
! 151: : mpn_mul (tmp, cp, cn, bp, bn);
1.1 maekawa 152:
153: /* now tmp[0]..tmp[k-1] contains the product of both mantissa,
1.1.1.2 ! ohara 154: with tmp[k-1]>=2^(BITS_PER_MP_LIMB-2) */
! 155: b1 >>= BITS_PER_MP_LIMB - 1; /* msb from the product */
! 156:
! 157: tmp += k - tn;
! 158: if (b1 == 0)
! 159: mpn_lshift (tmp, tmp, tn, 1);
! 160: cc = mpfr_round_raw (ap, tmp, bq + cq, sign_product < 0, aq,
! 161: rnd_mode, &inexact);
! 162: if (cc) /* cc = 1 ==> result is a power of two */
! 163: ap[an-1] = GMP_LIMB_HIGHBIT;
! 164:
! 165: TMP_FREE(marker);
! 166:
! 167: ec = b1 - 1 + cc;
! 168:
! 169: if (em == 0)
! 170: {
! 171: mp_exp_t ax = bx + cx;
! 172:
! 173: if (ax == __mpfr_emax && ec > 0)
! 174: return mpfr_set_overflow(a, rnd_mode, sign_product);
! 175:
! 176: if (ax == __mpfr_emin && ec < 0)
! 177: return mpfr_set_underflow(a, rnd_mode, sign_product);
! 178:
! 179: MPFR_EXP(a) = ax + ec;
! 180: }
! 181: else if (em > 0)
! 182: {
! 183: if (ec >= 0)
! 184: return mpfr_set_overflow(a, rnd_mode, sign_product);
! 185:
! 186: MPFR_EXP(a) = __mpfr_emax;
! 187: }
! 188: else
! 189: {
! 190: if (ec <= 0)
! 191: return mpfr_set_underflow(a, rnd_mode, sign_product);
! 192:
! 193: MPFR_EXP(a) = __mpfr_emin;
! 194: }
! 195:
! 196: if (MPFR_SIGN(a) != sign_product)
! 197: MPFR_CHANGE_SIGN(a);
! 198:
! 199: return inexact;
1.1 maekawa 200: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to mul.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpfr