gcc/mpfr/mul.c - Issue 3050029: [gcc] GCC 4.5.0=>4.5.1

Unified Diff: gcc/mpfr/mul.c

Issue 3050029: [gcc] GCC 4.5.0=>4.5.1 (Closed) Base URL: ssh://git@gitrw.chromium.org:9222/nacl-toolchain.git

Patch Set: Created 10 years, 5 months ago

Use n/p to move between diff chunks; N/P to move between comments. Draft comments are only viewable by you.

Jump to:

View side-by-side diff with in-line comments

Index: gcc/mpfr/mul.c

diff --git a/gcc/mpfr/mul.c b/gcc/mpfr/mul.c

deleted file mode 100644

index ec1f757e5381514c20961d1810f3cfd4fcdc5d16..0000000000000000000000000000000000000000

--- a/gcc/mpfr/mul.c

+++ /dev/null

@@ -1,513 +0,0 @@

-/* mpfr_mul -- multiply two floating-point numbers

-Contributed by the Arenaire and Cacao projects, INRIA.

-This file is part of the GNU MPFR Library.

-The GNU MPFR Library is free software; you can redistribute it and/or modify

-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 GNU 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 Lesser General Public

-License for more details.

-You should have received a copy of the GNU Lesser General Public License

-along with the GNU MPFR Library; see the file COPYING.LIB. If not, write to

-the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,

-MA 02110-1301, USA. */

-#define MPFR_NEED_LONGLONG_H

-#include "mpfr-impl.h"

-/********* BEGINNING CHECK *************/

-/* Check if we have to check the result of mpfr_mul.

- TODO: Find a better (and faster?) check than using old implementation */

-#ifdef WANT_ASSERT

-# if WANT_ASSERT >= 3

-int mpfr_mul2 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode);

-static int

-mpfr_mul3 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode)

- /* Old implementation */

- int sign_product, cc, inexact;

- mp_exp_t ax;

- mp_limb_t *tmp;

- mp_limb_t b1;

- mp_prec_t bq, cq;

- mp_size_t bn, cn, tn, k;

- MPFR_TMP_DECL(marker);

- /* deal with special cases */

- if (MPFR_ARE_SINGULAR(b,c))

- {

- if (MPFR_IS_NAN(b) || MPFR_IS_NAN(c))

- {

- MPFR_SET_NAN(a);

- MPFR_RET_NAN;

- }

- sign_product = MPFR_MULT_SIGN( MPFR_SIGN(b) , MPFR_SIGN(c) );

- if (MPFR_IS_INF(b))

- {

- if (MPFR_IS_INF(c) || MPFR_NOTZERO(c))

- {

- MPFR_SET_SIGN(a,sign_product);

- MPFR_SET_INF(a);

- MPFR_RET(0); /* exact */

- }

- else

- {

- MPFR_SET_NAN(a);

- MPFR_RET_NAN;

- }

- else if (MPFR_IS_INF(c))

- {

- if (MPFR_NOTZERO(b))

- {

- MPFR_SET_SIGN(a, sign_product);

- MPFR_SET_INF(a);

- MPFR_RET(0); /* exact */

- }

- else

- {

- MPFR_SET_NAN(a);

- MPFR_RET_NAN;

- }

- else

- {

- MPFR_ASSERTD(MPFR_IS_ZERO(b) || MPFR_IS_ZERO(c));

- MPFR_SET_SIGN(a, sign_product);

- MPFR_SET_ZERO(a);

- MPFR_RET(0); /* 0 * 0 is exact */

- }

- MPFR_CLEAR_FLAGS(a);

- sign_product = MPFR_MULT_SIGN( MPFR_SIGN(b) , MPFR_SIGN(c) );

- ax = MPFR_GET_EXP (b) + MPFR_GET_EXP (c);

- bq = MPFR_PREC(b);

- cq = MPFR_PREC(c);

- MPFR_ASSERTD(bq+cq > bq); /* PREC_MAX is /2 so no integer overflow */

- bn = (bq+BITS_PER_MP_LIMB-1)/BITS_PER_MP_LIMB; /* number of limbs of b */

- cn = (cq+BITS_PER_MP_LIMB-1)/BITS_PER_MP_LIMB; /* number of limbs of c */

- k = bn + cn; /* effective nb of limbs used by b*c (= tn or tn+1) below */

- tn = (bq + cq + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB;

- /* <= k, thus no int overflow */

- MPFR_ASSERTD(tn <= k);

- /* Check for no size_t overflow*/

- MPFR_ASSERTD((size_t) k <= ((size_t) -1) / BYTES_PER_MP_LIMB);

- MPFR_TMP_MARK(marker);

- tmp = (mp_limb_t *) MPFR_TMP_ALLOC((size_t) k * BYTES_PER_MP_LIMB);

- /* multiplies two mantissa in temporary allocated space */

- b1 = (MPFR_LIKELY(bn >= cn)) ?

- mpn_mul (tmp, MPFR_MANT(b), bn, MPFR_MANT(c), cn)

- : mpn_mul (tmp, MPFR_MANT(c), cn, MPFR_MANT(b), bn);

- /* now tmp[0]..tmp[k-1] contains the product of both mantissa,

- with tmp[k-1]>=2^(BITS_PER_MP_LIMB-2) */

- b1 >>= BITS_PER_MP_LIMB - 1; /* msb from the product */

- /* if the mantissas of b and c are uniformly distributed in ]1/2, 1],

- then their product is in ]1/4, 1/2] with probability 2*ln(2)-1 ~ 0.386

- and in [1/2, 1] with probability 2-2*ln(2) ~ 0.614 */

- tmp += k - tn;

- if (MPFR_UNLIKELY(b1 == 0))

- mpn_lshift (tmp, tmp, tn, 1); /* tn <= k, so no stack corruption */

- cc = mpfr_round_raw (MPFR_MANT (a), tmp, bq + cq,

- MPFR_IS_NEG_SIGN(sign_product),

- MPFR_PREC (a), rnd_mode, &inexact);

- /* cc = 1 ==> result is a power of two */

- if (MPFR_UNLIKELY(cc))

- MPFR_MANT(a)[MPFR_LIMB_SIZE(a)-1] = MPFR_LIMB_HIGHBIT;

- MPFR_TMP_FREE(marker);

- {

- mp_exp_t ax2 = ax + (mp_exp_t) (b1 - 1 + cc);

- if (MPFR_UNLIKELY( ax2 > __gmpfr_emax))

- return mpfr_overflow (a, rnd_mode, sign_product);

- if (MPFR_UNLIKELY( ax2 < __gmpfr_emin))

- {

- /* In the rounding to the nearest mode, if the exponent of the exact

- result (i.e. before rounding, i.e. without taking cc into account)

- is < __gmpfr_emin - 1 or the exact result is a power of 2 (i.e. if

- both arguments are powers of 2), then round to zero. */

- if (rnd_mode == GMP_RNDN &&

- (ax + (mp_exp_t) b1 < __gmpfr_emin ||

- (mpfr_powerof2_raw (b) && mpfr_powerof2_raw (c))))

- rnd_mode = GMP_RNDZ;

- return mpfr_underflow (a, rnd_mode, sign_product);

- }

- MPFR_SET_EXP (a, ax2);

- MPFR_SET_SIGN(a, sign_product);

- }

- MPFR_RET (inexact);

-int

-mpfr_mul (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode)

- mpfr_t ta, tb, tc;

- int inexact1, inexact2;

- mpfr_init2 (ta, MPFR_PREC (a));

- mpfr_init2 (tb, MPFR_PREC (b));

- mpfr_init2 (tc, MPFR_PREC (c));

- MPFR_ASSERTN (mpfr_set (tb, b, GMP_RNDN) == 0);

- MPFR_ASSERTN (mpfr_set (tc, c, GMP_RNDN) == 0);

- inexact2 = mpfr_mul3 (ta, tb, tc, rnd_mode);

- inexact1 = mpfr_mul2 (a, b, c, rnd_mode);

- if (mpfr_cmp (ta, a) || inexact1*inexact2 < 0

- || (inexact1*inexact2 == 0 && (inexact1|inexact2) != 0))

- {

- fprintf (stderr, "mpfr_mul return different values for %s\n"

- "Prec_a = %lu, Prec_b = %lu, Prec_c = %lu\nB = ",

- mpfr_print_rnd_mode (rnd_mode),

- MPFR_PREC (a), MPFR_PREC (b), MPFR_PREC (c));

- mpfr_out_str (stderr, 16, 0, tb, GMP_RNDN);

- fprintf (stderr, "\nC = ");

- mpfr_out_str (stderr, 16, 0, tc, GMP_RNDN);

- fprintf (stderr, "\nOldMul: ");

- mpfr_out_str (stderr, 16, 0, ta, GMP_RNDN);

- fprintf (stderr, "\nNewMul: ");

- mpfr_out_str (stderr, 16, 0, a, GMP_RNDN);

- fprintf (stderr, "\nNewInexact = %d | OldInexact = %d\n",

- inexact1, inexact2);

- MPFR_ASSERTN(0);

- }

- mpfr_clears (ta, tb, tc, (mpfr_ptr) 0);

- return inexact1;

-# define mpfr_mul mpfr_mul2

-# endif

-#endif

-/****** END OF CHECK *******/

-/* Multiply 2 mpfr_t */

-int

-mpfr_mul (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode)

- int sign, inexact;

- mp_exp_t ax, ax2;

- mp_limb_t *tmp;

- mp_limb_t b1;

- mp_prec_t bq, cq;

- mp_size_t bn, cn, tn, k;

- MPFR_TMP_DECL (marker);

- MPFR_LOG_FUNC (("b[%#R]=%R c[%#R]=%R rnd=%d", b, b, c, c, rnd_mode),

- ("a[%#R]=%R inexact=%d", a, a, inexact));

- /* deal with special cases */

- if (MPFR_ARE_SINGULAR (b, c))

- {

- if (MPFR_IS_NAN (b) || MPFR_IS_NAN (c))

- {

- MPFR_SET_NAN (a);

- MPFR_RET_NAN;

- }

- sign = MPFR_MULT_SIGN (MPFR_SIGN (b), MPFR_SIGN (c));

- if (MPFR_IS_INF (b))

- {

- if (!MPFR_IS_ZERO (c))

- {

- MPFR_SET_SIGN (a, sign);

- MPFR_SET_INF (a);

- MPFR_RET (0);

- }

- else

- {

- MPFR_SET_NAN (a);

- MPFR_RET_NAN;

- }

- else if (MPFR_IS_INF (c))

- {

- if (!MPFR_IS_ZERO (b))

- {

- MPFR_SET_SIGN (a, sign);

- MPFR_SET_INF (a);

- MPFR_RET(0);

- }

- else

- {

- MPFR_SET_NAN (a);

- MPFR_RET_NAN;

- }

- else

- {

- MPFR_ASSERTD (MPFR_IS_ZERO(b) || MPFR_IS_ZERO(c));

- MPFR_SET_SIGN (a, sign);

- MPFR_SET_ZERO (a);

- MPFR_RET (0);

- }

- MPFR_CLEAR_FLAGS (a);

- sign = MPFR_MULT_SIGN (MPFR_SIGN (b), MPFR_SIGN (c));

- ax = MPFR_GET_EXP (b) + MPFR_GET_EXP (c);

- /* Note: the exponent of the exact result will be e = bx + cx + ec with

- ec in {-1,0,1} and the following assumes that e is representable. */

- /* FIXME: Useful since we do an exponent check after ?

- * It is useful iff the precision is big, there is an overflow

- * and we are doing further mults...*/

-#ifdef HUGE

- if (MPFR_UNLIKELY (ax > __gmpfr_emax + 1))

- return mpfr_overflow (a, rnd_mode, sign);

- if (MPFR_UNLIKELY (ax < __gmpfr_emin - 2))

- return mpfr_underflow (a, rnd_mode == GMP_RNDN ? GMP_RNDZ : rnd_mode,

- sign);

-#endif

- bq = MPFR_PREC (b);

- cq = MPFR_PREC (c);

- MPFR_ASSERTD (bq+cq > bq); /* PREC_MAX is /2 so no integer overflow */

- bn = (bq+BITS_PER_MP_LIMB-1)/BITS_PER_MP_LIMB; /* number of limbs of b */

- cn = (cq+BITS_PER_MP_LIMB-1)/BITS_PER_MP_LIMB; /* number of limbs of c */

- k = bn + cn; /* effective nb of limbs used by b*c (= tn or tn+1) below */

- tn = (bq + cq + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB;

- MPFR_ASSERTD (tn <= k); /* tn <= k, thus no int overflow */

- /* Check for no size_t overflow*/

- MPFR_ASSERTD ((size_t) k <= ((size_t) -1) / BYTES_PER_MP_LIMB);

- MPFR_TMP_MARK (marker);

- tmp = (mp_limb_t *) MPFR_TMP_ALLOC ((size_t) k * BYTES_PER_MP_LIMB);

- /* multiplies two mantissa in temporary allocated space */

- if (MPFR_UNLIKELY (bn < cn))

- {

- mpfr_srcptr z = b;

- mp_size_t zn = bn;

- b = c;

- bn = cn;

- c = z;

- cn = zn;

- }

- MPFR_ASSERTD (bn >= cn);

- if (MPFR_LIKELY (bn <= 2))

- {

- if (bn == 1)

- {

- /* 1 limb * 1 limb */

- umul_ppmm (tmp[1], tmp[0], MPFR_MANT (b)[0], MPFR_MANT (c)[0]);

- b1 = tmp[1];

- }

- else if (MPFR_UNLIKELY (cn == 1))

- {

- /* 2 limbs * 1 limb */

- mp_limb_t t;

- umul_ppmm (tmp[1], tmp[0], MPFR_MANT (b)[0], MPFR_MANT (c)[0]);

- umul_ppmm (tmp[2], t, MPFR_MANT (b)[1], MPFR_MANT (c)[0]);

- add_ssaaaa (tmp[2], tmp[1], tmp[2], tmp[1], 0, t);

- b1 = tmp[2];

- }

- else

- {

- /* 2 limbs * 2 limbs */

- mp_limb_t t1, t2, t3;

- /* First 2 limbs * 1 limb */

- umul_ppmm (tmp[1], tmp[0], MPFR_MANT (b)[0], MPFR_MANT (c)[0]);

- umul_ppmm (tmp[2], t1, MPFR_MANT (b)[1], MPFR_MANT (c)[0]);

- add_ssaaaa (tmp[2], tmp[1], tmp[2], tmp[1], 0, t1);

- /* Second, the other 2 limbs * 1 limb product */

- umul_ppmm (t1, t2, MPFR_MANT (b)[0], MPFR_MANT (c)[1]);

- umul_ppmm (tmp[3], t3, MPFR_MANT (b)[1], MPFR_MANT (c)[1]);

- add_ssaaaa (tmp[3], t1, tmp[3], t1, 0, t3);

- /* Sum those two partial products */

- add_ssaaaa (tmp[2], tmp[1], tmp[2], tmp[1], t1, t2);

- tmp[3] += (tmp[2] < t1);

- b1 = tmp[3];

- }

- b1 >>= (BITS_PER_MP_LIMB - 1);

- tmp += k - tn;

- if (MPFR_UNLIKELY (b1 == 0))

- mpn_lshift (tmp, tmp, tn, 1); /* tn <= k, so no stack corruption */

- }

- else

- /* Mulders' mulhigh. Disable if squaring, since it is not tuned for

- such a case */

- if (MPFR_UNLIKELY (bn > MPFR_MUL_THRESHOLD && b != c))

- {

- mp_limb_t *bp, *cp;

- mp_size_t n;

- mp_prec_t p;

- /* Fist check if we can reduce the precision of b or c:

- exact values are a nightmare for the short product trick */

- bp = MPFR_MANT (b);

- cp = MPFR_MANT (c);

- MPFR_ASSERTN (MPFR_MUL_THRESHOLD >= 1);

- if (MPFR_UNLIKELY ((bp[0] == 0 && bp[1] == 0) ||

- (cp[0] == 0 && cp[1] == 0)))

- {

- mpfr_t b_tmp, c_tmp;

- MPFR_TMP_FREE (marker);

- /* Check for b */

- while (*bp == 0)

- {

- bp++;

- bn--;

- MPFR_ASSERTD (bn > 0);

- } /* This must end since the MSL is != 0 */

- /* Check for c too */

- while (*cp == 0)

- {

- cp++;

- cn--;

- MPFR_ASSERTD (cn > 0);

- } /* This must end since the MSL is != 0 */

- /* It is not the faster way, but it is safer */

- MPFR_SET_SAME_SIGN (b_tmp, b);

- MPFR_SET_EXP (b_tmp, MPFR_GET_EXP (b));

- MPFR_PREC (b_tmp) = bn * BITS_PER_MP_LIMB;

- MPFR_MANT (b_tmp) = bp;

- MPFR_SET_SAME_SIGN (c_tmp, c);

- MPFR_SET_EXP (c_tmp, MPFR_GET_EXP (c));

- MPFR_PREC (c_tmp) = cn * BITS_PER_MP_LIMB;

- MPFR_MANT (c_tmp) = cp;

- /* Call again mpfr_mul with the fixed arguments */

- return mpfr_mul (a, b_tmp, c_tmp, rnd_mode);

- }

- /* Compute estimated precision of mulhigh.

- We could use `+ (n < cn) + (n < bn)' instead of `+ 2',

- but does it worth it? */

- n = MPFR_LIMB_SIZE (a) + 1;

- n = MIN (n, cn);

- MPFR_ASSERTD (n >= 1 && 2*n <= k && n <= cn && n <= bn);

- p = n * BITS_PER_MP_LIMB - MPFR_INT_CEIL_LOG2 (n + 2);

- bp += bn - n;

- cp += cn - n;

- /* Check if MulHigh can produce a roundable result.

- We may lost 1 bit due to RNDN, 1 due to final shift. */

- if (MPFR_UNLIKELY (MPFR_PREC (a) > p - 5))

- {

- if (MPFR_UNLIKELY (MPFR_PREC (a) > p - 5 + BITS_PER_MP_LIMB

- || bn <= MPFR_MUL_THRESHOLD+1))

- {

- /* MulHigh can't produce a roundable result. */

- MPFR_LOG_MSG (("mpfr_mulhigh can't be used (%lu VS %lu)\n",

- MPFR_PREC (a), p));

- goto full_multiply;

- }

- /* Add one extra limb to mantissa of b and c. */

- if (bn > n)

- bp --;

- else

- {

- bp = (mp_limb_t*) MPFR_TMP_ALLOC ((n+1) * sizeof (mp_limb_t));

- bp[0] = 0;

- MPN_COPY (bp + 1, MPFR_MANT (b) + bn - n, n);

- }

- if (cn > n)

- cp --; /* FIXME: Could this happen? */

- else

- {

- cp = (mp_limb_t*) MPFR_TMP_ALLOC ((n+1) * sizeof (mp_limb_t));

- cp[0] = 0;

- MPN_COPY (cp + 1, MPFR_MANT (c) + cn - n, n);

- }

- /* We will compute with one extra limb */

- n++;

- p = n * BITS_PER_MP_LIMB - MPFR_INT_CEIL_LOG2 (n + 2);

- /* Due to some nasty reasons we can have only 4 bits */

- MPFR_ASSERTD (MPFR_PREC (a) <= p - 4);

- if (MPFR_LIKELY (k < 2*n))

- {

- tmp = (mp_limb_t*) MPFR_TMP_ALLOC (2 * n * sizeof (mp_limb_t));

- tmp += 2*n-k; /* `tmp' still points to an area of `k' limbs */

- }

- MPFR_LOG_MSG (("Use mpfr_mulhigh (%lu VS %lu)\n", MPFR_PREC (a), p));

- /* Compute an approximation of the product of b and c */

- mpfr_mulhigh_n (tmp + k - 2 * n, bp, cp, n);

- /* now tmp[0]..tmp[k-1] contains the product of both mantissa,

- with tmp[k-1]>=2^(BITS_PER_MP_LIMB-2) */

- b1 = tmp[k-1] >> (BITS_PER_MP_LIMB - 1); /* msb from the product */

- /* If the mantissas of b and c are uniformly distributed in (1/2, 1],

- then their product is in (1/4, 1/2] with probability 2*ln(2)-1

- ~ 0.386 and in [1/2, 1] with probability 2-2*ln(2) ~ 0.614 */

- if (MPFR_UNLIKELY (b1 == 0))

- /* Warning: the mpfr_mulhigh_n call above only surely affects

- tmp[k-n-1..k-1], thus we shift only those limbs */

- mpn_lshift (tmp + k - n - 1, tmp + k - n - 1, n + 1, 1);

- tmp += k - tn;

- MPFR_ASSERTD (MPFR_LIMB_MSB (tmp[tn-1]) != 0);

- if (MPFR_UNLIKELY (!mpfr_round_p (tmp, tn, p+b1-1, MPFR_PREC(a)

- + (rnd_mode == GMP_RNDN))))

- {

- tmp -= k - tn; /* tmp may have changed, FIX IT!!!!! */

- goto full_multiply;

- }

- else

- {

- full_multiply:

- MPFR_LOG_MSG (("Use mpn_mul\n", 0));

- b1 = mpn_mul (tmp, MPFR_MANT (b), bn, MPFR_MANT (c), cn);

- /* now tmp[0]..tmp[k-1] contains the product of both mantissa,

- with tmp[k-1]>=2^(BITS_PER_MP_LIMB-2) */

- b1 >>= BITS_PER_MP_LIMB - 1; /* msb from the product */

- /* if the mantissas of b and c are uniformly distributed in (1/2, 1],

- then their product is in (1/4, 1/2] with probability 2*ln(2)-1

- ~ 0.386 and in [1/2, 1] with probability 2-2*ln(2) ~ 0.614 */

- tmp += k - tn;

- if (MPFR_UNLIKELY (b1 == 0))

- mpn_lshift (tmp, tmp, tn, 1); /* tn <= k, so no stack corruption */

- }

- ax2 = ax + (mp_exp_t) (b1 - 1);

- MPFR_RNDRAW (inexact, a, tmp, bq+cq, rnd_mode, sign, ax2++);

- MPFR_TMP_FREE (marker);

- MPFR_EXP (a) = ax2; /* Can't use MPFR_SET_EXP: Expo may be out of range */

- MPFR_SET_SIGN (a, sign);

- if (MPFR_UNLIKELY (ax2 > __gmpfr_emax))

- return mpfr_overflow (a, rnd_mode, sign);

- if (MPFR_UNLIKELY (ax2 < __gmpfr_emin))

- {

- /* In the rounding to the nearest mode, if the exponent of the exact

- result (i.e. before rounding, i.e. without taking cc into account)

- is < __gmpfr_emin - 1 or the exact result is a power of 2 (i.e. if

- both arguments are powers of 2), then round to zero. */

- if (rnd_mode == GMP_RNDN

- && (ax + (mp_exp_t) b1 < __gmpfr_emin

- || (mpfr_powerof2_raw (b) && mpfr_powerof2_raw (c))))

- rnd_mode = GMP_RNDZ;

- return mpfr_underflow (a, rnd_mode, sign);

- }

- MPFR_RET (inexact);

« no previous file with comments | « gcc/mpfr/mpn_exp.c ('k') | gcc/mpfr/mul_2exp.c » ('j') | no next file with comments »