Avoid excess range overflowing results from cosh, sinh, lgamma (bug 18980).

Various i386 libm functions return values with excess range and
precision; Wilco Dijkstra's patches to make isfinite etc. expand
inline cause this pre-existing issue to result in test failures (when
e.g. a result that overflows float but not long double gets counted as
overflowing for some purposes but not others).

This patch addresses those cases arising from functions defined in C,
adding a math_narrow_eval macro that forces values to memory to
eliminate excess precision if FLT_EVAL_METHOD indicates this is
needed, and is a no-op otherwise.  I'll convert existing uses of
volatile and asm for this purpose to use the new macro later, once
i386 has clean test results again (which requires fixes for .S files
as well).

Tested for x86_64 and x86.  Committed.

	[BZ #18980]
	* sysdeps/generic/math_private.h: Include <float.h>.
	(math_narrow_eval): New macro.
	[FLT_EVAL_METHOD != 0] (excess_precision): Likewise.
	* sysdeps/ieee754/dbl-64/e_cosh.c (__ieee754_cosh): Use
	math_narrow_eval on overflowing return value.
	* sysdeps/ieee754/dbl-64/e_lgamma_r.c (__ieee754_lgamma_r):
	Likewise.
	* sysdeps/ieee754/dbl-64/e_sinh.c (__ieee754_sinh): Likewise.
	* sysdeps/ieee754/flt-32/e_coshf.c (__ieee754_coshf): Likewise.
	* sysdeps/ieee754/flt-32/e_lgammaf_r.c (__ieee754_lgammaf_r):
	Likewise.
	* sysdeps/ieee754/flt-32/e_sinhf.c (__ieee754_sinhf): Likewise.

diff --git a/ChangeLog b/ChangeLog
index 9a89e9d9e5652d32a7102924556097ad21746d23..57b07b36ee85ec8abf055629af382109d7dc12dc 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,19 @@
+2015-09-18  Joseph Myers  <joseph@codesourcery.com>
+
+       [BZ #18980]
+       * sysdeps/generic/math_private.h: Include <float.h>.
+       (math_narrow_eval): New macro.
+       [FLT_EVAL_METHOD != 0] (excess_precision): Likewise.
+       * sysdeps/ieee754/dbl-64/e_cosh.c (__ieee754_cosh): Use
+       math_narrow_eval on overflowing return value.
+       * sysdeps/ieee754/dbl-64/e_lgamma_r.c (__ieee754_lgamma_r):
+       Likewise.
+       * sysdeps/ieee754/dbl-64/e_sinh.c (__ieee754_sinh): Likewise.
+       * sysdeps/ieee754/flt-32/e_coshf.c (__ieee754_coshf): Likewise.
+       * sysdeps/ieee754/flt-32/e_lgammaf_r.c (__ieee754_lgammaf_r):
+       Likewise.
+       * sysdeps/ieee754/flt-32/e_sinhf.c (__ieee754_sinhf): Likewise.
+
 2015-09-18  Wilco Dijkstra  <wdijkstr@arm.com>
 
        * include/math.h: Remove __isinf_ns, __isinf_nsf, __isinf_nsl.

diff --git a/sysdeps/generic/math_private.h b/sysdeps/generic/math_private.h
index 6aea8643daa5e118003a602742a61698cecae7bb..35591eed157c12afd6f92faa6a9a9e3e3c4731b1 100644 (file)
--- a/sysdeps/generic/math_private.h
+++ b/sysdeps/generic/math_private.h
@@ -20,6 +20,7 @@
 #include <stdint.h>
 #include <sys/types.h>
 #include <fenv.h>
+#include <float.h>
 #include <get-rounding-mode.h>
 
 /* The original fdlibm code used statements like:
@@ -405,6 +406,29 @@ extern long double __lgamma_productl (long double t, long double x,
 ({ __typeof (x) __x = (x); __asm __volatile__ ("" : : "m" (__x)); })
 #endif
 
+/* math_narrow_eval reduces its floating-point argument to the range
+   and precision of its semantic type.  (The original evaluation may
+   still occur with excess range and precision, so the result may be
+   affected by double rounding.)  */
+#if FLT_EVAL_METHOD == 0
+# define math_narrow_eval(x) (x)
+#else
+# if FLT_EVAL_METHOD == 1
+#  define excess_precision(type) __builtin_types_compatible_p (type, float)
+# else
+#  define excess_precision(type) (__builtin_types_compatible_p (type, float) \
+                                 || __builtin_types_compatible_p (type, \
+                                                                  double))
+# endif
+# define math_narrow_eval(x)                                   \
+  ({                                                           \
+    __typeof (x) math_narrow_eval_tmp = (x);                   \
+    if (excess_precision (__typeof (math_narrow_eval_tmp)))    \
+      __asm__ ("" : "+m" (math_narrow_eval_tmp));              \
+    math_narrow_eval_tmp;                                      \
+   })
+#endif
+
 
 /* The standards only specify one variant of the fenv.h interfaces.
    But at least for some architectures we can be more efficient if we

diff --git a/sysdeps/ieee754/dbl-64/e_cosh.c b/sysdeps/ieee754/dbl-64/e_cosh.c
index af3910dd6ebdb6b64881c07c366c1121046a0c97..52a5d5007da2258ca38118d4ac9df81ac2ddfdb8 100644 (file)
--- a/sysdeps/ieee754/dbl-64/e_cosh.c
+++ b/sysdeps/ieee754/dbl-64/e_cosh.c
@@ -83,6 +83,6 @@ __ieee754_cosh (double x)
     return x * x;
 
   /* |x| > overflowthresold, cosh(x) overflow */
-  return huge * huge;
+  return math_narrow_eval (huge * huge);
 }
 strong_alias (__ieee754_cosh, __cosh_finite)

diff --git a/sysdeps/ieee754/dbl-64/e_lgamma_r.c b/sysdeps/ieee754/dbl-64/e_lgamma_r.c
index ea8a9b42fb70fdeab5d9d7a87261c9907b79ce4c..da158cba3363d0b8d5bb136d70074d17e605d6c4 100644 (file)
--- a/sysdeps/ieee754/dbl-64/e_lgamma_r.c
+++ b/sysdeps/ieee754/dbl-64/e_lgamma_r.c
@@ -296,7 +296,7 @@ __ieee754_lgamma_r(double x, int *signgamp)
            r = (x-half)*(t-one)+w;
        } else
     /* 2**58 <= x <= inf */
-           r =  x*(__ieee754_log(x)-one);
+           r =  math_narrow_eval (x*(__ieee754_log(x)-one));
        /* NADJ is set for negative arguments but not otherwise,
           resulting in warnings that it may be used uninitialized
           although in the cases where it is used it has always been

diff --git a/sysdeps/ieee754/dbl-64/e_sinh.c b/sysdeps/ieee754/dbl-64/e_sinh.c
index c99d28311d39d65e35747998b870fbd9cf056424..291bfad0b331ac6e2a79a92a57912b9fcd4a9016 100644 (file)
--- a/sysdeps/ieee754/dbl-64/e_sinh.c
+++ b/sysdeps/ieee754/dbl-64/e_sinh.c
@@ -89,6 +89,6 @@ __ieee754_sinh (double x)
     }
 
   /* |x| > overflowthresold, sinh(x) overflow */
-  return x * shuge;
+  return math_narrow_eval (x * shuge);
 }
 strong_alias (__ieee754_sinh, __sinh_finite)

diff --git a/sysdeps/ieee754/flt-32/e_coshf.c b/sysdeps/ieee754/flt-32/e_coshf.c
index dedda47c093dfd0f107808607d1fe131da28ca71..7b223758e161e570f838aaa2392ce4a665b32785 100644 (file)
--- a/sysdeps/ieee754/flt-32/e_coshf.c
+++ b/sysdeps/ieee754/flt-32/e_coshf.c
@@ -58,6 +58,6 @@ __ieee754_coshf (float x)
        if(ix>=0x7f800000) return x*x;
 
     /* |x| > overflowthresold, cosh(x) overflow */
-       return huge*huge;
+       return math_narrow_eval (huge*huge);
 }
 strong_alias (__ieee754_coshf, __coshf_finite)

diff --git a/sysdeps/ieee754/flt-32/e_lgammaf_r.c b/sysdeps/ieee754/flt-32/e_lgammaf_r.c
index 424c4e735820818c109a8ca27083f4671bc96be9..45a62c0dabdc092f45da482ba3d40730e1f0e04d 100644 (file)
--- a/sysdeps/ieee754/flt-32/e_lgammaf_r.c
+++ b/sysdeps/ieee754/flt-32/e_lgammaf_r.c
@@ -232,7 +232,7 @@ __ieee754_lgammaf_r(float x, int *signgamp)
            r = (x-half)*(t-one)+w;
        } else
     /* 2**26 <= x <= inf */
-           r =  x*(__ieee754_logf(x)-one);
+           r =  math_narrow_eval (x*(__ieee754_logf(x)-one));
        /* NADJ is set for negative arguments but not otherwise,
           resulting in warnings that it may be used uninitialized
           although in the cases where it is used it has always been

diff --git a/sysdeps/ieee754/flt-32/e_sinhf.c b/sysdeps/ieee754/flt-32/e_sinhf.c
index 17c2219c0ba2154feda10f7221b28bb48839ad8e..a24fa0c4bce6de1b10f8227e1de7ed8f3cd0e2c9 100644 (file)
--- a/sysdeps/ieee754/flt-32/e_sinhf.c
+++ b/sysdeps/ieee754/flt-32/e_sinhf.c
@@ -59,6 +59,6 @@ __ieee754_sinhf(float x)
        }
 
     /* |x| > overflowthresold, sinh(x) overflow */
-       return x*shuge;
+       return math_narrow_eval (x*shuge);
 }
 strong_alias (__ieee754_sinhf, __sinhf_finite)