A patch for libm-ieee754

[hjl@lucon.org (H.J. Lu)]

Hi,

It turns out those libm-ieee754 bugs are not in egcs. Here is the patch
for glibc 2.1 to fix a few libm-ieee754 bugs. Ulrich, could you please
take a look?

When you use fpu/cpu to do rounding, you have to mark the variable
volatile. Otherwise, the compiler may do some thing you don't want.

Thanks.


-- 
H.J. Lu (hjl@gnu.org)
----
Sat Feb 21 18:47:44 1998  H.J. Lu  (hjl@gnu.org)

	* sysdeps/libm-ieee754/e_exp.c (__ieee754_exp): Changed type
	of TWO43, TWO52 from float to double. Make 'n' and 't' volatile.
	Use __isinf.

	* sysdeps/libm-ieee754/e_expf.c (__ieee754_expf): Make 'n' and
	't' volatile. Use __isinff.

	* sysdeps/libm-ieee754/s_exp2.c (__ieee754_exp2): Changed type
	of TWO43 from float to double. Make "rx" volatile. Use __isinf.

	* sysdeps/libm-ieee754/s_exp2f.c (__ieee754_exp2f): Make "rx"
	volatile. Use __isinff.

Index: sysdeps/libm-ieee754/e_exp.c
===================================================================
RCS file: /home/work/cvs/gnu/glibc/sysdeps/libm-ieee754/e_exp.c,v
retrieving revision 1.1.1.2
diff -u -p -r1.1.1.2 e_exp.c
--- e_exp.c	1998/02/14 03:01:26	1.1.1.2
+++ e_exp.c	1998/02/22 02:43:20
@@ -71,14 +71,13 @@ static const volatile double TWOM1000 = 
 double
 __ieee754_exp (double x)
 {
-  static const uint32_t a_minf = 0xff800000;
   static const double himark = 709.7827128933840868;
   static const double lomark = -745.1332191019412221;
   /* Check for usual case.  */
   if (isless (x, himark) && isgreater (x, lomark))
     {
-      static const float TWO43 = 8796093022208.0;
-      static const float TWO52 = 4503599627370496.0;
+      static const double TWO43 = 8796093022208.0;
+      static const double TWO52 = 4503599627370496.0;
       /* 1/ln(2).  */
       static const double M_1_LN2 = 1.442695040888963387;
       /* ln(2), part 1 */
@@ -87,7 +86,8 @@ __ieee754_exp (double x)
       static const double M_LN2_1 = 5.497923018708371155e-14;
 
       int tval, unsafe, n_i;
-      double x22, n, t, dely, result;
+      double x22, dely, result;
+      volatile double n, t;
       union ieee754_double ex2_u, scale_u;
       fenv_t oldenv;
 
@@ -166,7 +166,7 @@ __ieee754_exp (double x)
   /* Exceptional cases:  */
   else if (isless (x, himark))
     {
-      if (x == *(const float *) &a_minf)
+      if (__isinf (x))
 	/* e^-inf == 0, with no error.  */
 	return 0;
       else
Index: sysdeps/libm-ieee754/e_expf.c
===================================================================
RCS file: /home/work/cvs/gnu/glibc/sysdeps/libm-ieee754/e_expf.c,v
retrieving revision 1.1.1.2
diff -u -p -r1.1.1.2 e_expf.c
--- e_expf.c	1998/02/14 03:01:27	1.1.1.2
+++ e_expf.c	1998/02/22 02:43:11
@@ -66,7 +66,6 @@ static const volatile float TWO127 = 1.7
 float
 __ieee754_expf (float x)
 {
-  static const uint32_t a_minf = 0xff800000;
   static const float himark = 88.72283935546875;
   static const float lomark = -103.972084045410;
   /* Check for usual case.  */
@@ -82,8 +81,10 @@ __ieee754_expf (float x)
       static const double M_LN2 = .6931471805599452862;
 
       int tval;
-      double x22, t, result, dx;
-      float n, delta;
+      double x22, result, dx;
+      volatile double t;
+      float delta;
+      volatile float n;
       union ieee754_double ex2_u;
       fenv_t oldenv;
 
@@ -144,7 +145,7 @@ __ieee754_expf (float x)
   /* Exceptional cases:  */
   else if (isless (x, himark))
     {
-      if (x == *(const float *) &a_minf)
+      if (__isinff (x))
 	/* e^-inf == 0, with no error.  */
 	return 0;
       else
Index: sysdeps/libm-ieee754/s_exp2.c
===================================================================
RCS file: /home/work/cvs/gnu/glibc/sysdeps/libm-ieee754/s_exp2.c,v
retrieving revision 1.1.1.3
diff -u -p -r1.1.1.3 s_exp2.c
--- s_exp2.c	1998/02/14 03:01:32	1.1.1.3
+++ s_exp2.c	1998/02/22 02:43:24
@@ -42,16 +42,16 @@ static const volatile double TWOM1000 = 
 double
 __ieee754_exp2 (double x)
 {
-  static const uint32_t a_minf = 0xff800000;
   static const double himark = (double) DBL_MAX_EXP;
   static const double lomark = (double) (DBL_MIN_EXP - DBL_MANT_DIG - 1) - 1.0;
 
   /* Check for usual case.  */
   if (isless (x, himark) && isgreater (x, lomark))
     {
-      static const float TWO43 = 8796093022208.0;
+      static const double TWO43 = 8796093022208.0;
       int tval, unsafe;
-      double rx, x22, result;
+      double x22, result;
+      volatile double rx;
       union ieee754_double ex2_u, scale_u;
       fenv_t oldenv;
 
@@ -125,7 +125,7 @@ __ieee754_exp2 (double x)
   /* Exceptional cases:  */
   else if (isless (x, himark))
     {
-      if (x == *(const float *) &a_minf)
+      if (__isinf (x))
 	/* e^-inf == 0, with no error.  */
 	return 0;
       else
Index: sysdeps/libm-ieee754/s_exp2f.c
===================================================================
RCS file: /home/work/cvs/gnu/glibc/sysdeps/libm-ieee754/s_exp2f.c,v
retrieving revision 1.1.1.4
diff -u -p -r1.1.1.4 s_exp2f.c
--- s_exp2f.c	1998/02/14 03:01:33	1.1.1.4
+++ s_exp2f.c	1998/02/22 02:43:16
@@ -43,7 +43,6 @@ static const volatile float TWO127 = 1.7
 float
 __ieee754_exp2f (float x)
 {
-  static const uint32_t a_minf = 0xff800000;
   static const float himark = (float) FLT_MAX_EXP;
   static const float lomark = (float) (FLT_MIN_EXP - FLT_MANT_DIG - 1) - 1.0;
 
@@ -52,7 +51,8 @@ __ieee754_exp2f (float x)
     {
       static const float TWO15 = 32768.0;
       int tval, unsafe;
-      float rx, x22, result;
+      float x22, result;
+      volatile float rx;
       union ieee754_float ex2_u, scale_u;
       fenv_t oldenv;
 
@@ -123,7 +123,7 @@ __ieee754_exp2f (float x)
   /* Exceptional cases:  */
   else if (isless (x, himark))
     {
-      if (x == *(const float *) &a_minf)
+      if (__isinff (x))
 	/* e^-inf == 0, with no error.  */
 	return 0;
       else