diff options
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/fma.c')
-rw-r--r-- | Build/source/libs/mpfr/mpfr-src/src/fma.c | 297 |
1 files changed, 0 insertions, 297 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/fma.c b/Build/source/libs/mpfr/mpfr-src/src/fma.c deleted file mode 100644 index 44a50553b2b..00000000000 --- a/Build/source/libs/mpfr/mpfr-src/src/fma.c +++ /dev/null @@ -1,297 +0,0 @@ -/* mpfr_fma -- Floating multiply-add - -Copyright 2001-2002, 2004, 2006-2020 Free Software Foundation, Inc. -Contributed by the AriC and Caramba 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 3 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.LESSER. If not, see -https://www.gnu.org/licenses/ or 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" - -/* The fused-multiply-add (fma) of x, y and z is defined by: - fma(x,y,z)= x*y + z -*/ - -/* this function deals with all cases where inputs are singular, i.e., - either NaN, Inf or zero */ -static int -mpfr_fma_singular (mpfr_ptr s, mpfr_srcptr x, mpfr_srcptr y, mpfr_srcptr z, - mpfr_rnd_t rnd_mode) -{ - if (MPFR_IS_NAN(x) || MPFR_IS_NAN(y) || MPFR_IS_NAN(z)) - { - MPFR_SET_NAN(s); - MPFR_RET_NAN; - } - /* now neither x, y or z is NaN */ - else if (MPFR_IS_INF(x) || MPFR_IS_INF(y)) - { - /* cases Inf*0+z, 0*Inf+z, Inf-Inf */ - if ((MPFR_IS_ZERO(y)) || - (MPFR_IS_ZERO(x)) || - (MPFR_IS_INF(z) && - ((MPFR_MULT_SIGN(MPFR_SIGN(x), MPFR_SIGN(y))) != MPFR_SIGN(z)))) - { - MPFR_SET_NAN(s); - MPFR_RET_NAN; - } - else if (MPFR_IS_INF(z)) /* case Inf-Inf already checked above */ - { - MPFR_SET_INF(s); - MPFR_SET_SAME_SIGN(s, z); - MPFR_RET(0); - } - else /* z is finite */ - { - MPFR_SET_INF(s); - MPFR_SET_SIGN(s, MPFR_MULT_SIGN(MPFR_SIGN(x) , MPFR_SIGN(y))); - MPFR_RET(0); - } - } - /* now x and y are finite */ - else if (MPFR_IS_INF(z)) - { - MPFR_SET_INF(s); - MPFR_SET_SAME_SIGN(s, z); - MPFR_RET(0); - } - else if (MPFR_IS_ZERO(x) || MPFR_IS_ZERO(y)) - { - if (MPFR_IS_ZERO(z)) - { - int sign_p; - sign_p = MPFR_MULT_SIGN( MPFR_SIGN(x) , MPFR_SIGN(y) ); - MPFR_SET_SIGN(s, (rnd_mode != MPFR_RNDD ? - (MPFR_IS_NEG_SIGN(sign_p) && MPFR_IS_NEG(z) ? - MPFR_SIGN_NEG : MPFR_SIGN_POS) : - (MPFR_IS_POS_SIGN(sign_p) && MPFR_IS_POS(z) ? - MPFR_SIGN_POS : MPFR_SIGN_NEG))); - MPFR_SET_ZERO(s); - MPFR_RET(0); - } - else - return mpfr_set (s, z, rnd_mode); - } - else /* necessarily z is zero here */ - { - MPFR_ASSERTD(MPFR_IS_ZERO(z)); - return mpfr_mul (s, x, y, rnd_mode); - } -} - -/* s <- x*y + z */ -int -mpfr_fma (mpfr_ptr s, mpfr_srcptr x, mpfr_srcptr y, mpfr_srcptr z, - mpfr_rnd_t rnd_mode) -{ - int inexact; - mpfr_t u; - mp_size_t n; - mpfr_exp_t e; - mpfr_prec_t precx, precy; - MPFR_SAVE_EXPO_DECL (expo); - MPFR_GROUP_DECL(group); - - MPFR_LOG_FUNC - (("x[%Pu]=%.*Rg y[%Pu]=%.*Rg z[%Pu]=%.*Rg rnd=%d", - mpfr_get_prec (x), mpfr_log_prec, x, - mpfr_get_prec (y), mpfr_log_prec, y, - mpfr_get_prec (z), mpfr_log_prec, z, rnd_mode), - ("s[%Pu]=%.*Rg inexact=%d", - mpfr_get_prec (s), mpfr_log_prec, s, inexact)); - - /* particular cases */ - if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x) || MPFR_IS_SINGULAR(y) || - MPFR_IS_SINGULAR(z) )) - return mpfr_fma_singular (s, x, y, z, rnd_mode); - - e = MPFR_GET_EXP (x) + MPFR_GET_EXP (y); - - precx = MPFR_PREC (x); - precy = MPFR_PREC (y); - - /* First deal with special case where prec(x) = prec(y) and x*y does - not overflow nor underflow. Do it only for small sizes since for large - sizes x*y is faster using Mulders' algorithm (as a rule of thumb, - we assume mpn_mul_n is faster up to 4*MPFR_MUL_THRESHOLD). - Since |EXP(x)|, |EXP(y)| < 2^(k-2) on a k-bit computer, - |EXP(x)+EXP(y)| < 2^(k-1), thus cannot overflow nor underflow. */ - if (precx == precy && e <= __gmpfr_emax && e > __gmpfr_emin) - { - if (precx < GMP_NUMB_BITS && - MPFR_PREC(z) == precx && - MPFR_PREC(s) == precx) - { - mp_limb_t umant[2], zmant[2]; - mpfr_t zz; - int inex; - - umul_ppmm (umant[1], umant[0], MPFR_MANT(x)[0], MPFR_MANT(y)[0]); - MPFR_PREC(u) = MPFR_PREC(zz) = 2 * precx; - MPFR_MANT(u) = umant; - MPFR_MANT(zz) = zmant; - MPFR_SIGN(u) = MPFR_MULT_SIGN( MPFR_SIGN(x) , MPFR_SIGN(y) ); - MPFR_SIGN(zz) = MPFR_SIGN(z); - MPFR_EXP(zz) = MPFR_EXP(z); - if (MPFR_PREC(zz) <= GMP_NUMB_BITS) /* zz fits in one limb */ - { - if ((umant[1] & MPFR_LIMB_HIGHBIT) == 0) - { - umant[0] = umant[1] << 1; - MPFR_EXP(u) = e - 1; - } - else - { - umant[0] = umant[1]; - MPFR_EXP(u) = e; - } - zmant[0] = MPFR_MANT(z)[0]; - } - else - { - zmant[1] = MPFR_MANT(z)[0]; - zmant[0] = MPFR_LIMB_ZERO; - if ((umant[1] & MPFR_LIMB_HIGHBIT) == 0) - { - umant[1] = (umant[1] << 1) | - (umant[0] >> (GMP_NUMB_BITS - 1)); - umant[0] = umant[0] << 1; - MPFR_EXP(u) = e - 1; - } - else - MPFR_EXP(u) = e; - } - inex = mpfr_add (u, u, zz, rnd_mode); - /* mpfr_set_1_2 requires PREC(u) = 2*PREC(s), - thus we need PREC(s) = PREC(x) = PREC(y) = PREC(z) */ - return mpfr_set_1_2 (s, u, rnd_mode, inex); - } - else if ((n = MPFR_LIMB_SIZE(x)) <= 4 * MPFR_MUL_THRESHOLD) - { - mpfr_limb_ptr up; - mp_size_t un = n + n; - MPFR_TMP_DECL(marker); - - MPFR_TMP_MARK(marker); - MPFR_TMP_INIT (up, u, un * GMP_NUMB_BITS, un); - up = MPFR_MANT(u); - /* multiply x*y exactly into u */ - mpn_mul_n (up, MPFR_MANT(x), MPFR_MANT(y), n); - if (MPFR_LIMB_MSB (up[un - 1]) == 0) - { - mpn_lshift (up, up, un, 1); - MPFR_EXP(u) = e - 1; - } - else - MPFR_EXP(u) = e; - MPFR_SIGN(u) = MPFR_MULT_SIGN( MPFR_SIGN(x) , MPFR_SIGN(y) ); - /* The above code does not generate any exception. - The exceptions will come only from mpfr_add. */ - inexact = mpfr_add (s, u, z, rnd_mode); - MPFR_TMP_FREE(marker); - return inexact; - } - } - - /* If we take prec(u) >= prec(x) + prec(y), the product u <- x*y - is exact, except in case of overflow or underflow. */ - MPFR_ASSERTN (precx + precy <= MPFR_PREC_MAX); - MPFR_GROUP_INIT_1 (group, precx + precy, u); - MPFR_SAVE_EXPO_MARK (expo); - - if (MPFR_UNLIKELY (mpfr_mul (u, x, y, MPFR_RNDN))) - { - /* overflow or underflow - this case is regarded as rare, thus - does not need to be very efficient (even if some tests below - could have been done earlier). - It is an overflow iff u is an infinity (since MPFR_RNDN was used). - Alternatively, we could test the overflow flag, but in this case, - mpfr_clear_flags would have been necessary. */ - - if (MPFR_IS_INF (u)) /* overflow */ - { - int sign_u = MPFR_SIGN (u); - - MPFR_LOG_MSG (("Overflow on x*y\n", 0)); - MPFR_GROUP_CLEAR (group); /* we no longer need u */ - - /* Let's eliminate the obvious case where x*y and z have the - same sign. No possible cancellation -> real overflow. - Also, we know that |z| < 2^emax. If E(x) + E(y) >= emax+3, - then |x*y| >= 2^(emax+1), and |x*y + z| > 2^emax. This case - is also an overflow. */ - if (sign_u == MPFR_SIGN (z) || e >= __gmpfr_emax + 3) - { - MPFR_SAVE_EXPO_FREE (expo); - return mpfr_overflow (s, rnd_mode, sign_u); - } - } - else /* underflow: one has |x*y| < 2^(emin-1). */ - { - MPFR_LOG_MSG (("Underflow on x*y\n", 0)); - - /* Easy cases: when 2^(emin-1) <= 1/2 * min(ulp(z),ulp(s)), - one can replace x*y by sign(x*y) * 2^(emin-1). Note that - this is even true in case of equality for MPFR_RNDN thanks - to the even-rounding rule. - The + 1 on MPFR_PREC (s) is necessary because the exponent - of the result can be EXP(z) - 1. */ - if (MPFR_GET_EXP (z) - __gmpfr_emin >= - MAX (MPFR_PREC (z), MPFR_PREC (s) + 1)) - { - MPFR_PREC (u) = MPFR_PREC_MIN; - mpfr_setmin (u, __gmpfr_emin); - MPFR_SET_SIGN (u, MPFR_MULT_SIGN (MPFR_SIGN (x), - MPFR_SIGN (y))); - mpfr_clear_flags (); - goto add; - } - - MPFR_GROUP_CLEAR (group); /* we no longer need u */ - } - - /* Let's use UBF to resolve the overflow/underflow issues. */ - { - mpfr_ubf_t uu; - mp_size_t un; - mpfr_limb_ptr up; - MPFR_TMP_DECL(marker); - - MPFR_LOG_MSG (("Use UBF\n", 0)); - - MPFR_TMP_MARK (marker); - un = MPFR_LIMB_SIZE (x) + MPFR_LIMB_SIZE (y); - MPFR_TMP_INIT (up, uu, (mpfr_prec_t) un * GMP_NUMB_BITS, un); - mpfr_ubf_mul_exact (uu, x, y); - mpfr_clear_flags (); - inexact = mpfr_add (s, (mpfr_srcptr) uu, z, rnd_mode); - MPFR_UBF_CLEAR_EXP (uu); - MPFR_TMP_FREE (marker); - } - } - else - { - add: - inexact = mpfr_add (s, u, z, rnd_mode); - MPFR_GROUP_CLEAR (group); - } - - MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags); - MPFR_SAVE_EXPO_FREE (expo); - return mpfr_check_range (s, inexact, rnd_mode); -} |