From e5aea4d602551285165acfdaf31552ecfba20866 Mon Sep 17 00:00:00 2001
From: dannysmith <dannysmith>
Date: Wed, 28 Jul 2004 11:25:26 +0000
Subject: [PATCH] 	* mingwex/math/powl.c (powl): Revert change of
 2004-02-01. 	(__convert_inf_to_maxnum): New.static inline. 	(reducl): Use
 it to protect against Inf - Inf. 	(__fast_ldexpl): New function.  Use in
 lieu of ldexpl.

---
 winsup/mingw/ChangeLog           |  7 ++++++
 winsup/mingw/mingwex/math/powl.c | 49 ++++++++++++++++++++++++++++------------
 2 files changed, 42 insertions(+), 14 deletions(-)

diff --git a/winsup/mingw/ChangeLog b/winsup/mingw/ChangeLog
index dbdeb0cdbd..a7615aa4dc 100644
--- a/winsup/mingw/ChangeLog
+++ b/winsup/mingw/ChangeLog
@@ -1,3 +1,10 @@
+2004-07-28  Danny Smith  <dannysmith@users.sourceforge.net>
+
+	* mingwex/math/powl.c (powl): Revert change of 2004-02-01.
+	(__convert_inf_to_maxnum): New.static inline.
+	(reducl): Use it to protect against Inf - Inf.
+	(__fast_ldexpl): New function.  Use in lieu of ldexpl.
+
 2004-07-27  Danny Smith  <dannysmith@users.sourceforge.net>
 	
 	* mingwex/math/expl.c (expl): Move body of code to new static
diff --git a/winsup/mingw/mingwex/math/powl.c b/winsup/mingw/mingwex/math/powl.c
index 19910bd667..f85e556536 100644
--- a/winsup/mingw/mingwex/math/powl.c
+++ b/winsup/mingw/mingwex/math/powl.c
@@ -411,7 +411,7 @@ long double floorl(), fabsl(), frexpl(), ldexpl();
 long double polevll(), p1evll(), __powil();
 static long double reducl();
 int isnanl(), isfinitel(), signbitl();
-#endif
+#endif  /* __MINGW32__ */
 
 #ifdef INFINITIES
 extern long double INFINITYL;
@@ -428,6 +428,24 @@ extern long double NEGZEROL;
 
 #endif /* __MINGW32__ */
 
+#ifdef __MINGW32__
+
+/* No error checking. We handle Infs and zeros ourselves.  */
+static __inline__ long double
+__fast_ldexpl (long double x, int expn)
+{
+  long double res;
+  __asm__ ("fscale"
+  	    : "=t" (res)
+	    : "0" (x), "u" ((long double) expn));
+  return res;
+}
+
+#define ldexpl  __fast_ldexpl
+
+#endif
+
+
 long double powl( x, y )
 long double x, y;
 {
@@ -690,25 +708,16 @@ Fa = reducl(F);
 Fb = F - Fa;
 
 G = Fa + w * ya;
-if (isinf (G))
-  {
-   /* Bail out: G - reducl(G) will result in NAN
-      that will propagate through rest of calculations */ 
-    _SET_ERRNO (ERANGE);
-    mtherr( fname, OVERFLOW );
-    return( MAXNUML );
-  }
 Ga = reducl(G);
 Gb = G - Ga;
 
 H = Fb + Gb;
 Ha = reducl(H);
-w = ldexpl( Ga+Ha, LNXT );
+w = ldexpl( Ga + Ha, LNXT );
 
 /* Test the power of 2 for overflow */
 if( w > MEXP )
 	{
-/*	printf( "w = %.4Le ", w ); */
 	_SET_ERRNO (ERANGE);
 	mtherr( fname, OVERFLOW );
 	return( MAXNUML );
@@ -716,7 +725,6 @@ if( w > MEXP )
 
 if( w < MNEXP )
 	{
-/*	printf( "w = %.4Le ", w ); */
 	_SET_ERRNO (ERANGE);
 	mtherr( fname, UNDERFLOW );
 	return( 0.0L );
@@ -768,6 +776,15 @@ if( nflg )
 return( z );
 }
 
+static __inline__ long double 
+__convert_inf_to_maxnum(long double x)
+{
+  if (isinf(x))
+    return (x > 0.0L ? MAXNUML : -MAXNUML);
+  else
+    return x;
+}
+
 
 /* Find a multiple of 1/NXT that is within 1/NXT of x. */
 static __inline__ long double reducl(x)
@@ -775,9 +792,13 @@ long double x;
 {
 long double t;
 
-t = ldexpl( x, LNXT );
+/* If the call to ldexpl overflows, set it to MAXNUML.
+   This avoids Inf - Inf = Nan result when calculating the 'small'
+   part of a reduction.  Instead, the small part becomes Inf,
+   causing under/overflow when adding it to the 'large' part.
+   There must be a cleaner way of doing this. */
+t =  __convert_inf_to_maxnum (ldexpl( x, LNXT ));
 t = floorl( t );
 t = ldexpl( t, -LNXT );
 return(t);
 }
-
-- 
2.11.0