@@ -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
@@ -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)
@@ -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
@@ -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)
@@ -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)
@@ -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
@@ -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)