Avoid overflows from long double functions using __kernel_standard.

4 files changed, 37 insertions, 4 deletions

diff --git a/sysdeps/generic/math_private.h b/sysdeps/generic/math_private.h
index 813ad93..e217224 100644
--- a/sysdeps/generic/math_private.h
+++ b/sysdeps/generic/math_private.h
@@ -217,6 +217,7 @@ extern double __ieee754_scalb (double,double);
 /* fdlibm kernel function */
 extern double __kernel_standard (double,double,int);
 extern float __kernel_standard_f (float,float,int);
+extern long double __kernel_standard_l (long double,long double,int);
 extern double __kernel_sin (double,double,int);
 extern double __kernel_cos (double,double);
 extern double __kernel_tan (double,double,int);
diff --git a/sysdeps/ieee754/k_standard.c b/sysdeps/ieee754/k_standard.c
index 5d84543..c3326d9 100644
--- a/sysdeps/ieee754/k_standard.c
+++ b/sysdeps/ieee754/k_standard.c
@@ -16,6 +16,7 @@ static char rcsid[] = "$NetBSD: k_standard.c,v 1.6 1995/05/10 20:46:35 jtc Exp $
 
 #include <math.h>
 #include <math_private.h>
+#include <float.h>
 #include <errno.h>
 
 #include <assert.h>
@@ -998,3 +999,34 @@ __kernel_standard_f(float x, float y, int type)
 {
 	return __kernel_standard(x, y, type);
 }
+
+long double
+__kernel_standard_l (long double x, long double y, int type)
+{
+  double dx, dy;
+  if (isfinite (x))
+    {
+      long double ax = fabsl (x);
+      if (ax > DBL_MAX)
+	dx = __copysignl (DBL_MAX, x);
+      else if (ax > 0 && ax < DBL_MIN)
+	dx = __copysignl (DBL_MIN, x);
+      else
+	dx = x;
+    }
+  else
+    dx = x;
+  if (isfinite (y))
+    {
+      long double ay = fabsl (y);
+      if (ay > DBL_MAX)
+	dy = __copysignl (DBL_MAX, y);
+      else if (ay > 0 && ay < DBL_MIN)
+	dy = __copysignl (DBL_MIN, y);
+      else
+	dy = y;
+    }
+  else
+    dy = y;
+  return __kernel_standard (dx, dy, type);
+}
diff --git a/sysdeps/ieee754/ldbl-128/w_expl.c b/sysdeps/ieee754/ldbl-128/w_expl.c
index f4deda8..10193be 100644
--- a/sysdeps/ieee754/ldbl-128/w_expl.c
+++ b/sysdeps/ieee754/ldbl-128/w_expl.c
@@ -39,9 +39,9 @@ long double __expl(long double x)	/* wrapper exp */
 	if(_LIB_VERSION == _IEEE_) return z;
 	if(__finitel(x)) {
 	    if(x>o_threshold)
-	        return __kernel_standard(x,x,206); /* exp overflow */
+	        return __kernel_standard_l(x,x,206); /* exp overflow */
 	    else if(x<u_threshold)
-	        return __kernel_standard(x,x,207); /* exp underflow */
+	        return __kernel_standard_l(x,x,207); /* exp underflow */
 	}
 	return z;
 #endif
diff --git a/sysdeps/ieee754/ldbl-96/w_expl.c b/sysdeps/ieee754/ldbl-96/w_expl.c
index d61c0a3..55c6846 100644
--- a/sysdeps/ieee754/ldbl-96/w_expl.c
+++ b/sysdeps/ieee754/ldbl-96/w_expl.c
@@ -33,12 +33,12 @@ __expl (long double x)
   if (__builtin_expect (isgreater (x, o_threshold), 0))
     {
       if (_LIB_VERSION != _IEEE_)
-	return __kernel_standard (x, x, 206);
+	return __kernel_standard_l (x, x, 206);
     }
   else if (__builtin_expect (isless (x, u_threshold), 0))
     {
       if (_LIB_VERSION != _IEEE_)
-	return __kernel_standard (x, x, 207);
+	return __kernel_standard_l (x, x, 207);
     }
 
   return __ieee754_expl (x);