Port all the IEEE dependent remaining functions so that it uses a generic way to compute the result if it detects that the double is not in IEEE format ( _GMP_IEEE_FLOATS == 0).

MPFR now should work well on non-IEEE machines.


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

10 files changed, 123 insertions, 56 deletions

diff --git a/exceptions.c b/exceptions.c
index d35e6c10d..93b37d33a 100644
--- a/exceptions.c
+++ b/exceptions.c
@@ -179,7 +179,7 @@ mpfr_inexflag_p (void)
   return __gmpfr_flags & MPFR_FLAGS_INEXACT;
 }
 
-#undef mpfr_set_underflow
+/* #undef mpfr_set_underflow */
 
 /* Note: In the rounding to the nearest mode, mpfr_set_underflow
    always rounds away from 0. In this rounding mode, you must call
@@ -209,7 +209,7 @@ mpfr_set_underflow (mpfr_ptr x, mp_rnd_t rnd_mode, int sign)
   return sign > 0 ? inex : -inex;
 }
 
-#undef mpfr_set_overflow
+/* #undef mpfr_set_overflow */
 
 int
 mpfr_set_overflow (mpfr_ptr x, mp_rnd_t rnd_mode, int sign)
diff --git a/get_d.c b/get_d.c
index adca12929..0dcdf0358 100644
--- a/get_d.c
+++ b/get_d.c
@@ -28,14 +28,6 @@ MA 02111-1307, USA. */
 #include "mpfr.h"
 #include "mpfr-impl.h"
 
-#ifdef XDEBUG
-#undef _GMP_IEEE_FLOATS
-#endif
-
-#ifndef _GMP_IEEE_FLOATS
-#define _GMP_IEEE_FLOATS 0
-#endif
-
 /* "double" NaN and infinities are written as explicit bytes to be sure of
    getting what we want, and to be sure of not depending on libm.
 
diff --git a/mpfr-impl.h b/mpfr-impl.h
index 0a8878499..03a5293ba 100644
--- a/mpfr-impl.h
+++ b/mpfr-impl.h
@@ -264,6 +264,14 @@ long double __gmpfr_longdouble_volatile _MPFR_PROTO ((long double)) ATTRIBUTE_CO
 # endif
 #endif
 
+/* Debug non IEEE floats */
+#ifdef XDEBUG
+# undef _GMP_IEEE_FLOATS
+#endif
+#ifndef _GMP_IEEE_FLOATS
+# define _GMP_IEEE_FLOATS 0
+#endif
+
 /* We want to test this :
  *  (rnd == GMP_RNDU && test) || (rnd == RNDD && !test)
  * This macro does this test faster*/
@@ -421,28 +429,36 @@ void mpfr_inits2 _MPFR_PROTO ((mp_prec_t, mpfr_ptr, ...));
 void mpfr_inits _MPFR_PROTO ((mpfr_ptr, ...));
 void mpfr_clears _MPFR_PROTO ((mpfr_ptr, ...));
 
+
 int mpfr_set_underflow _MPFR_PROTO ((mpfr_ptr, mp_rnd_t, int));
 int mpfr_set_overflow _MPFR_PROTO ((mpfr_ptr, mp_rnd_t, int));
 void mpfr_save_emin_emax _MPFR_PROTO ((void));
 void mpfr_restore_emin_emax _MPFR_PROTO ((void));
+
 int mpfr_add1 _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mp_rnd_t));
 int mpfr_sub1 _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mp_rnd_t));
 int mpfr_can_round_raw _MPFR_PROTO ((mp_limb_t *, mp_size_t, int, mp_exp_t,
 				     mp_rnd_t, mp_rnd_t, mp_prec_t));
+
 double mpfr_get_d3 _MPFR_PROTO ((mpfr_srcptr, mp_exp_t, mp_rnd_t));
 int mpfr_cmp2 _MPFR_PROTO ((mpfr_srcptr, mpfr_srcptr, mp_prec_t *));
-long __gmpfr_ceil_log2 _MPFR_PROTO ((double));
-long __gmpfr_floor_log2 _MPFR_PROTO ((double));
-double __gmpfr_ceil_exp2 _MPFR_PROTO ((double));
+
+long          __gmpfr_ceil_log2 _MPFR_PROTO ((double));
+long          __gmpfr_floor_log2 _MPFR_PROTO ((double));
+double        __gmpfr_ceil_exp2 _MPFR_PROTO ((double));
 unsigned long __gmpfr_isqrt _MPFR_PROTO ((unsigned long));
 unsigned long __gmpfr_cuberoot _MPFR_PROTO ((unsigned long));
+
 int mpfr_exp_2 _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr, mp_rnd_t));
 int mpfr_exp3 _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr, mp_rnd_t));
 int mpfr_powerof2_raw _MPFR_PROTO ((mpfr_srcptr));
+
 void mpfr_setmax _MPFR_PROTO ((mpfr_ptr, mp_exp_t));
 void mpfr_setmin _MPFR_PROTO ((mpfr_ptr, mp_exp_t));
+
 long mpn_exp _MPFR_PROTO ((mp_limb_t *, mp_exp_t *, int, 
 			   mp_exp_t, size_t));
+
 void mpfr_print_binary _MPFR_PROTO ((mpfr_srcptr));
 void mpfr_set_str_binary _MPFR_PROTO ((mpfr_ptr, __gmp_const char *));
 
diff --git a/neg.c b/neg.c
index 185a836c1..a4f6c40cf 100644
--- a/neg.c
+++ b/neg.c
@@ -27,7 +27,7 @@ MA 02111-1307, USA. */
 int
 mpfr_neg (mpfr_ptr a, mpfr_srcptr b, mp_rnd_t rnd_mode)
 {
-  if (a != b)
+  if (MPFR_UNLIKELY(a != b))
     return mpfr_set4 (a, b, rnd_mode, -MPFR_SIGN(b));
   else
     {
diff --git a/powerof2.c b/powerof2.c
index f4bf0d3c6..f549fdfa3 100644
--- a/powerof2.c
+++ b/powerof2.c
@@ -32,7 +32,7 @@ mpfr_powerof2_raw (mpfr_srcptr x)
   mp_limb_t *xp;
   mp_size_t xn;
 
-  MPFR_ASSERTN(MPFR_IS_FP(x));
+  MPFR_ASSERTN(MPFR_IS_PURE_FP(x));
   xp = MPFR_MANT(x);
   xn = (MPFR_PREC(x) - 1) / BITS_PER_MP_LIMB;
   if (NOT_POW2(xp[xn]))
diff --git a/set_d.c b/set_d.c
index 7e5672b14..79cd23af9 100644
--- a/set_d.c
+++ b/set_d.c
@@ -34,16 +34,6 @@ MA 02111-1307, USA. */
 #error "Unsupported value of BITS_PER_MP_LIMB"
 #endif
 
-/* Included from gmp-2.0.2, patched to support denorms */
-
-#ifdef XDEBUG
-#undef _GMP_IEEE_FLOATS
-#endif
-
-#ifndef _GMP_IEEE_FLOATS
-#define _GMP_IEEE_FLOATS 0
-#endif
-
 static int
 __mpfr_extract_double (mp_ptr rp, double d)
      /* e=0 iff BITS_PER_MP_LIMB=32 and rp has only one limb */
diff --git a/set_ld.c b/set_ld.c
index 055be494d..b99d7188d 100644
--- a/set_ld.c
+++ b/set_ld.c
@@ -120,35 +120,35 @@ mpfr_set_ld (mpfr_ptr r, long double d, mp_rnd_t rnd_mode)
           /* since -DBL_MAX <= d <= DBL_MAX, the cast to double should not
              overflow here */
           mpfr_set_d (u, (double) d, GMP_RNDN);
-      /* warning: using MPFR_IS_ZERO will cause a READ_UNINIT_MEM if u=Inf */
-      if (mpfr_cmp_ui (u, 0) == 0 && (d != (long double) 0.0)) /* underflow */
-        {
-          long double div10, div11, div12, div13;
-          div10 = (long double) (double) 5.5626846462680034577255e-309; /* 2^(-2^10) */
-          div11 = div10 * div10; /* 2^(-2^11) */
-          div12 = div11 * div11; /* 2^(-2^12) */
-          div13 = div12 * div12; /* 2^(-2^13) */
-          if (ABS(d) <= div13)
-            {
-              d = d / div13; /* exact */
-              shift_exp -= 8192;
-            }
-          if (ABS(d) <= div12)
-            {
-              d = d / div12; /* exact */
-              shift_exp -= 4096;
-            }
-          if (ABS(d) <= div11)
-            {
-              d = d / div11; /* exact */
-              shift_exp -= 2048;
-            }
-          if (ABS(d) <= div10)
-            {
-              d = d / div10; /* exact */
-              shift_exp -= 1024;
-            }
-        }
+	  /* warning: using MPFR_IS_ZERO will cause a READ_UNINIT_MEM if u=Inf */
+	  if (mpfr_cmp_ui (u, 0) == 0 && (d != (long double) 0.0)) /* underflow */
+	    {
+	      long double div10, div11, div12, div13;
+	      div10 = (long double) (double) 5.5626846462680034577255e-309; /* 2^(-2^10) */
+	      div11 = div10 * div10; /* 2^(-2^11) */
+	      div12 = div11 * div11; /* 2^(-2^12) */
+	      div13 = div12 * div12; /* 2^(-2^13) */
+	      if (ABS(d) <= div13)
+		{
+		  d = d / div13; /* exact */
+		  shift_exp -= 8192;
+		}
+	      if (ABS(d) <= div12)
+		{
+		  d = d / div12; /* exact */
+		  shift_exp -= 4096;
+		}
+	      if (ABS(d) <= div11)
+		{
+		  d = d / div11; /* exact */
+		  shift_exp -= 2048;
+		}
+	      if (ABS(d) <= div10)
+		{
+		  d = d / div10; /* exact */
+		  shift_exp -= 1024;
+		}
+	    }
         }
       mpfr_add (t, t, u, GMP_RNDN); /* exact */
       if (!mpfr_number_p (t))
diff --git a/uceil_exp2.c b/uceil_exp2.c
index 8003a5164..16ad030d7 100644
--- a/uceil_exp2.c
+++ b/uceil_exp2.c
@@ -37,6 +37,25 @@ __gmpfr_ceil_exp2 (double d)
     exp++;
   /* now exp = ceil(d) */
   x.d = 1.0;
+#if _GMP_IEEE_FLOATS
   x.s.exp = exp <= -1022 ? 1 : 1023 + exp;
+#else
+  if (exp >= 0)
+    {
+      while (exp != 0)
+	{
+	  x.d *= 2.0;
+	  exp--;
+	}
+    }
+  else
+    {
+      while (exp != 0)
+	{
+	x.d *= (1.0 / 2.0);
+	exp++;
+	}
+    }
+#endif
   return x.d;
 }
diff --git a/uceil_log2.c b/uceil_log2.c
index f7ba603cd..ab3bcd8ab 100644
--- a/uceil_log2.c
+++ b/uceil_log2.c
@@ -24,11 +24,14 @@ MA 02111-1307, USA. */
 #include "mpfr.h"
 #include "mpfr-impl.h"
 
-/* returns ceil(log(d)/log(2)) */
+/* returns ceil(log(d)/log(2)) if d > 0,
+   -1023 if d = +0,
+   and floor(log(-d)/log(2))+1 if d < 0*/
 long
 __gmpfr_ceil_log2 (double d)
 {
   long exp;
+#if _GMP_IEEE_FLOATS
   union ieee_double_extract x;
 
   x.d = d;
@@ -37,4 +40,26 @@ __gmpfr_ceil_log2 (double d)
   if (x.d != 1.0) /* d: not a power of two? */
     exp++;
   return exp;
+#else
+  double m;
+
+  if (d < 0.0)
+    return __gmpfr_floor_log2(-d)+1;
+  else if (d == 0.0)
+    return -1023;
+  else if (d >= 1.0)
+    {
+      exp = 0;
+      for( m= 1.0 ; m < d ; m *=2.0 )
+	exp++;
+    }
+  else
+    {
+      exp = 1;
+      for( m= 1.0 ; m >= d ; m *= (1.0/2.0) )
+        exp--;      
+    }
+#endif
+  return exp;
 }
+
diff --git a/ufloor_log2.c b/ufloor_log2.c
index 7d8e05494..4b9c4a8a5 100644
--- a/ufloor_log2.c
+++ b/ufloor_log2.c
@@ -24,12 +24,37 @@ MA 02111-1307, USA. */
 #include "mpfr.h"
 #include "mpfr-impl.h"
 
-/* returns floor(log(d)/log(2)) */
+/* returns floor(log(abs(d))/log(2)) */
 long
 __gmpfr_floor_log2 (double d)
 {
+#if _GMP_IEEE_FLOATS
   union ieee_double_extract x;
 
   x.d = d;
   return (long) x.s.exp - 1023;
+#else  
+  long exp;
+  double m;
+
+  /* Get Abs */
+  if (d < 0.0)
+    d = -d;
+
+  if (d == 0.0)
+    return -1023;
+  else if (d >= 1.0)
+    {
+      exp = -1;
+      for( m= 1.0 ; m <= d ; m *=2.0 )
+	exp++;
+    }
+  else
+    {
+      exp = 0;
+      for( m= 1.0 ; m > d ; m *= (1.0/2.0) )
+        exp--;      
+    }
+  return exp;
+#endif
 }