put back simple algorithm that computes x*y exactly and then

directly calls mpfr_add, to avoid wrong inexact flags


git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@1780 280ebfd0-de03-0410-8827-d642c229c3f4

1 files changed, 8 insertions, 92 deletions

diff --git a/fma.c b/fma.c
index 1f4c8c174..fc542207a 100644
--- a/fma.c
+++ b/fma.c
@@ -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;       
-    mp_exp_t 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 = (mp_exp_t) 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;
 }