* mpfr: Update to 20020301, except internal_ceil_exp2.c,

	internal_ceil_log2.c, internal_floor_log2.c renamed to i_ceil_exp2.c,
	i_ceil_log2.c, i_floor_log2.c to be unique in DOS 8.3.  And sqrtrem.c
	removed since no longer required.

1 files changed, 160 insertions, 155 deletions

diff --git a/mpfr/fma.c b/mpfr/fma.c
index 0a6c1ce73..8abebc743 100644
--- a/mpfr/fma.c
+++ b/mpfr/fma.c
@@ -1,6 +1,6 @@
 /* mpfr_fma -- Floating multiply-add
 
-Copyright (C) 2001 Free Software Foundation, Inc.
+Copyright (C) 2001, 2002 Free Software Foundation, Inc.
 
 This file is part of the MPFR Library.
 
@@ -22,7 +22,6 @@ the file COPYING.LIB.  If not, write to the Free Software
 Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 MA 02111-1307, USA. */
 
-#include <stdio.h>
 #include "gmp.h"
 #include "gmp-impl.h"
 #include "mpfr.h"
@@ -34,167 +33,173 @@ MA 02111-1307, USA. */
 */
 
 int
-mpfr_fma (mpfr_ptr s, mpfr_srcptr x ,mpfr_srcptr y ,mpfr_srcptr z , mp_rnd_t rnd_mode) 
+mpfr_fma (mpfr_ptr s, mpfr_srcptr x, mpfr_srcptr y, mpfr_srcptr z,
+          mp_rnd_t rnd_mode)
 {
-  int inexact =0;
+  int inexact = 0;
   /* Flag calcul exacte */
-  int not_exact=0;
-
-        /* particular cases */
-        if (MPFR_IS_NAN(x) || MPFR_IS_NAN(y) ||  MPFR_IS_NAN(z))
-	  {
-	    MPFR_SET_NAN(s); 
-	    return 1; 
-	  }
-
-        /* cases Inf*0+z, 0*Inf+z, Inf-Inf */
-        if ((MPFR_IS_INF(x) && MPFR_IS_ZERO(y)) ||
-	    (MPFR_IS_INF(y) && MPFR_IS_ZERO(x)) ||
-	    ((MPFR_IS_INF(x) || MPFR_IS_INF(y)) && MPFR_IS_INF(z) &&
-	     ((MPFR_SIGN(x) * MPFR_SIGN(y)) != MPFR_SIGN(z))))
-	  {
-	    MPFR_SET_NAN(s);
-	    return 1;
-	  }
-
-	MPFR_CLEAR_NAN(s);
-
-        if (MPFR_IS_INF(x) || MPFR_IS_INF(y))
-	  {
-	    if (MPFR_IS_INF(z)) /* case Inf-Inf already checked above */
-	      {
-		MPFR_SET_INF(s);
-		MPFR_SET_SAME_SIGN(s, z);
-		return 0;
-	      }
-	    else /* z is finite */
-	      {
-		MPFR_SET_INF(s);
-		if (MPFR_SIGN(s) != (MPFR_SIGN(x) * MPFR_SIGN(y)))
-		  MPFR_CHANGE_SIGN(s);
-		return 0;
-	      }
-	  }
-
-	/* now x and y are finite */
-	if (MPFR_IS_INF(z))
-	  {
-	    MPFR_SET_INF(s);
-	    MPFR_SET_SAME_SIGN(s, z);
-	    return 0;
-	  }
-
-        MPFR_CLEAR_INF(s);
-
-        if(!MPFR_NOTZERO(x) || !MPFR_NOTZERO(y))
-	  return mpfr_set (s, z, rnd_mode);
-        if (!MPFR_NOTZERO(z))
-	  return mpfr_mul (s, x, y, rnd_mode);
-	
-        /* General case */
-        /* Detail of the compute */
-
-        /* u <- x*y */
-        /* t <- z+u */
+  int not_exact = 0;
+
+  /* particular cases */
+  if (MPFR_IS_NAN(x) || MPFR_IS_NAN(y) || MPFR_IS_NAN(z))
+    {
+      MPFR_SET_NAN(s);
+      MPFR_RET_NAN;
+    }
+
+  if (MPFR_IS_INF(x) || MPFR_IS_INF(y))
+    {
+      /* cases Inf*0+z, 0*Inf+z, Inf-Inf */
+      if ((MPFR_IS_FP(y) && MPFR_IS_ZERO(y)) ||
+          (MPFR_IS_FP(x) && MPFR_IS_ZERO(x)) ||
+          (MPFR_IS_INF(z) && ((MPFR_SIGN(x) * MPFR_SIGN(y)) != MPFR_SIGN(z))))
         {
-          /* Declaration of the intermediary variable */
-          mpfr_t t, u;       
-          int d,fl1,fl2;
-          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 */
-          int 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=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);
-              goto fma_paul;
-            }
-
-            not_exact=0;
-
-            /* reactualisation of the precision */
-            mpfr_set_prec(u,Nt);             
-            mpfr_set_prec(t,Nt);             
-            
-            /* computations */
-            fl1=mpfr_mul(u,x,y,GMP_RNDN);
-            if(fl1) not_exact=1;
-
-            fl2=mpfr_add(t,z,u,GMP_RNDN);
-            if(fl2) not_exact=1;        
-
-            /*Nt=Nt+(d+1)+_mpfr_ceil_log2(Nt); */
-            d = MPFR_EXP(u)-MPFR_EXP(t);
+          MPFR_SET_NAN(s);
+          MPFR_RET_NAN;
+        }
 
-            /* estimation of the error */
-            err=Nt-(d+1);
+      MPFR_CLEAR_NAN(s);
 
-            /* actualisation of the precision */
-            Nt +=(1-first_pass)*d + first_pass*10;
-            if(Nt<0)Nt=0;
+      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);
+          if (MPFR_SIGN(s) != (MPFR_SIGN(x) * MPFR_SIGN(y)))
+            MPFR_CHANGE_SIGN(s);
+          MPFR_RET(0);
+        }
+    }
 
-            first_pass=1;
+  MPFR_CLEAR_NAN(s);
 
-          } while ( (fl1!=fl2) || (err <0) || ((!mpfr_can_round(t,err,GMP_RNDN,rnd_mode,Ns)) && not_exact ));
+  /* now x and y are finite */
+  if (MPFR_IS_INF(z))
+    {
+      MPFR_SET_INF(s);
+      MPFR_SET_SAME_SIGN(s, z);
+      MPFR_RET(0);
+    }
 
-                inexact = mpfr_set (s, t, rnd_mode);
-                mpfr_clear(t);
-                mpfr_clear(u);
+  MPFR_CLEAR_INF(s);
 
-                goto fin;
+  if (MPFR_IS_ZERO(x) || MPFR_IS_ZERO(y))
+    {
+      if (MPFR_IS_ZERO(z))
+        {
+          int sign_p, sign_z;
+          sign_p = MPFR_SIGN(x) * MPFR_SIGN(y);
+          sign_z = MPFR_SIGN(z);
+          if (MPFR_SIGN(s) !=
+              (rnd_mode != GMP_RNDD ?
+               ((sign_p < 0 && sign_z < 0) ? -1 : 1) :
+               ((sign_p > 0 && sign_z > 0) ? 1 : -1)))
+            {
+              MPFR_CHANGE_SIGN(s);
+            }
+          MPFR_SET_ZERO(s);
+          MPFR_RET(0);
         }
-        
-
-
- fma_paul:
-
-
-        /* General case */
-        /* Detail of the compute */
-        /* u <- x*y exact */
-        /* s <- z+u */
-	{
-	  mpfr_t 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;
-
- fin:
-        if (not_exact == 0 && inexact == 0)
-          return 0;
-        
-        if (not_exact != 0 && inexact == 0)
-          return 1;
-        
-          return inexact;
-        
+      else
+        return mpfr_set (s, z, rnd_mode);
+    }
+
+  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;
+
+  return inexact;
 }
-