diff options
Diffstat (limited to 'mpfr/fma.c')
| -rw-r--r-- | mpfr/fma.c | 104 |
1 files changed, 10 insertions, 94 deletions
diff --git a/mpfr/fma.c b/mpfr/fma.c index 8abebc743..4ded04629 100644 --- a/mpfr/fma.c +++ b/mpfr/fma.c @@ -1,6 +1,6 @@ /* mpfr_fma -- Floating multiply-add -Copyright (C) 2001, 2002 Free Software Foundation, Inc. +Copyright 2001, 2002 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -18,7 +18,7 @@ License for more details. 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 +the file COPYING. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ @@ -36,9 +36,8 @@ int mpfr_fma (mpfr_ptr s, mpfr_srcptr x, mpfr_srcptr y, mpfr_srcptr z, mp_rnd_t rnd_mode) { - int inexact = 0; - /* Flag calcul exacte */ - int not_exact = 0; + int inexact; + mpfr_t u; /* particular cases */ if (MPFR_IS_NAN(x) || MPFR_IS_NAN(y) || MPFR_IS_NAN(z)) @@ -111,95 +110,12 @@ mpfr_fma (mpfr_ptr s, mpfr_srcptr x, mpfr_srcptr y, mpfr_srcptr z, if (MPFR_IS_ZERO(z)) return mpfr_mul (s, x, y, rnd_mode); - - /* General case */ - /* Detail of the compute */ - - /* u <- x*y */ - /* t <- z+u */ - { - /* Declaration of the intermediary variable */ - mpfr_t t, u; - int d; - int accu = 0; - - /* Declaration of the size variable */ - mp_prec_t Nx = MPFR_PREC(x); /* Precision of input variable */ - mp_prec_t Ny = MPFR_PREC(y); /* Precision of input variable */ - mp_prec_t Nz = MPFR_PREC(z); /* Precision of input variable */ - mp_prec_t Ns = MPFR_PREC(s); /* Precision of output variable */ - mp_prec_t Nt; /* Precision of the intermediary variable */ - long int err; /* Precision of error */ - unsigned int first_pass = 0; /* temporary precision */ - - /* compute the precision of intermediary variable */ - Nt = MAX(MAX(Nx,Ny),Nz); - - /* the optimal number of bits is MPFR_EXP(u)-MPFR_EXP(v)+1 */ - /* but u and v are not yet compute, also we take in account */ - /* just one bit */ - Nt += 1 + _mpfr_ceil_log2(Nt) + 20; - /* initialise the intermediary variables */ - mpfr_init(u); - mpfr_init(t); - - /* First computation of fma */ - do - { - if (accu++ > 2) - { - mpfr_clear(t); - mpfr_clear(u); - - /* General case */ - /* Detail of the compute */ - /* u <- x*y exact */ - /* s <- z+u */ - - /* if we take prec(u) >= prec(x) + prec(y), the product - u <- x*y is always exact */ - mpfr_init2 (u, MPFR_PREC(x) + MPFR_PREC(y)); - mpfr_mul (u, x, y, GMP_RNDN); /* always exact */ - inexact = mpfr_add (s, z, u, rnd_mode); - mpfr_clear(u); - return inexact; - } - - /* reactualisation of the precision */ - mpfr_set_prec(u, Nt); - mpfr_set_prec(t, Nt); - - /* computations */ - not_exact = mpfr_mul (u, x, y, GMP_RNDN); - - not_exact |= mpfr_add (t, z, u, GMP_RNDN); - - /* Nt = Nt + (d+1) + _mpfr_ceil_log2(Nt); */ - d = MPFR_EXP(u) - MPFR_EXP(t); - - /* estimate of the error */ - err = Nt - (d+1); - - /* actualisation of the precision */ - Nt += (1-first_pass) * d + first_pass * 10; - if (Nt < 0) - Nt = 0; - - first_pass = 1; - } - while (not_exact && - ((err < 0) || !mpfr_can_round (t, err, GMP_RNDN, rnd_mode, Ns))); - - inexact = mpfr_set (s, t, rnd_mode); - mpfr_clear(t); - mpfr_clear(u); - } - - if (not_exact == 0 && inexact == 0) - return 0; - - if (not_exact != 0 && inexact == 0) - return 1; + /* if we take prec(u) >= prec(x) + prec(y), the product + u <- x*y is always exact */ + mpfr_init2 (u, MPFR_PREC(x) + MPFR_PREC(y)); + mpfr_mul (u, x, y, GMP_RNDN); /* always exact */ + inexact = mpfr_add (s, z, u, rnd_mode); + mpfr_clear(u); return inexact; } |