math: fix sinh overflows in non-nearest rounding

The final rounding operation should be done with the correct sign
otherwise huge results may incorrectly get rounded to or away from
infinity in upward or downward rounding modes.

This affected sinh and sinhf which set the sign on the result after
a potentially overflowing mul. There may be other non-nearest rounding
issues, but this was a known long standing issue with large ulp error
(depending on how ulp is defined near infinity).

The fix should have no effect on sinh and sinhf performance but may
have a tiny effect on cosh and coshf.

diff --git a/src/internal/libm.h b/src/internal/libm.h
index b5bd26b8407751e7042f6ad974a79af08226d0f0..7533f6baefe64e28f594c66c3fed7617c83b354c 100644 (file)
--- a/src/internal/libm.h
+++ b/src/internal/libm.h
@@ -236,13 +236,13 @@ hidden int    __rem_pio2(double,double*);
 hidden double __sin(double,double,int);
 hidden double __cos(double,double);
 hidden double __tan(double,double,int);
-hidden double __expo2(double);
+hidden double __expo2(double,double);
 
 hidden int    __rem_pio2f(float,double*);
 hidden float  __sindf(double);
 hidden float  __cosdf(double);
 hidden float  __tandf(double,int);
-hidden float  __expo2f(float);
+hidden float  __expo2f(float,float);
 
 hidden int __rem_pio2l(long double, long double *);
 hidden long double __sinl(long double, long double, int);

diff --git a/src/math/__expo2.c b/src/math/__expo2.c
index 740ac680e857070c950d4e0d52d019ad4c225ab6..248f052b98b65ec64ba19943781b03e7ace02712 100644 (file)
--- a/src/math/__expo2.c
+++ b/src/math/__expo2.c
@@ -5,12 +5,13 @@ static const int k = 2043;
 static const double kln2 = 0x1.62066151add8bp+10;
 
 /* exp(x)/2 for x >= log(DBL_MAX), slightly better than 0.5*exp(x/2)*exp(x/2) */
-double __expo2(double x)
+double __expo2(double x, double sign)
 {
        double scale;
 
        /* note that k is odd and scale*scale overflows */
        INSERT_WORDS(scale, (uint32_t)(0x3ff + k/2) << 20, 0);
        /* exp(x - k ln2) * 2**(k-1) */
-       return exp(x - kln2) * scale * scale;
+       /* in directed rounding correct sign before rounding or overflow is important */
+       return exp(x - kln2) * (sign * scale) * scale;
 }

diff --git a/src/math/__expo2f.c b/src/math/__expo2f.c
index 5163e4180033b1ec45ac6cb5e159c4d607bb10dc..538eb09c0724c2583398054560a6261a6f770441 100644 (file)
--- a/src/math/__expo2f.c
+++ b/src/math/__expo2f.c
@@ -5,12 +5,13 @@ static const int k = 235;
 static const float kln2 = 0x1.45c778p+7f;
 
 /* expf(x)/2 for x >= log(FLT_MAX), slightly better than 0.5f*expf(x/2)*expf(x/2) */
-float __expo2f(float x)
+float __expo2f(float x, float sign)
 {
        float scale;
 
        /* note that k is odd and scale*scale overflows */
        SET_FLOAT_WORD(scale, (uint32_t)(0x7f + k/2) << 23);
        /* exp(x - k ln2) * 2**(k-1) */
-       return expf(x - kln2) * scale * scale;
+       /* in directed rounding correct sign before rounding or overflow is important */
+       return expf(x - kln2) * (sign * scale) * scale;
 }

diff --git a/src/math/cosh.c b/src/math/cosh.c
index 100f8231d8829604b9e293491a9abb1774b4b29a..490c15fb166752d18fa4015401f5a27fd34337b6 100644 (file)
--- a/src/math/cosh.c
+++ b/src/math/cosh.c
@@ -35,6 +35,6 @@ double cosh(double x)
 
        /* |x| > log(DBL_MAX) or nan */
        /* note: the result is stored to handle overflow */
-       t = __expo2(x);
+       t = __expo2(x, 1.0);
        return t;
 }

diff --git a/src/math/coshf.c b/src/math/coshf.c
index b09f2ee5751f4341cfb6cc759055279deda1d7cb..e739cff93b1325f931198d8d75b3313359803752 100644 (file)
--- a/src/math/coshf.c
+++ b/src/math/coshf.c
@@ -28,6 +28,6 @@ float coshf(float x)
        }
 
        /* |x| > log(FLT_MAX) or nan */
-       t = __expo2f(x);
+       t = __expo2f(x, 1.0f);
        return t;
 }

diff --git a/src/math/sinh.c b/src/math/sinh.c
index 00022c4e6ff6d8250ea238207a12199a1e349dac..a01951ae6fc4fdb4954de0b6f24323b230f6ae88 100644 (file)
--- a/src/math/sinh.c
+++ b/src/math/sinh.c
@@ -34,6 +34,6 @@ double sinh(double x)
 
        /* |x| > log(DBL_MAX) or nan */
        /* note: the result is stored to handle overflow */
-       t = 2*h*__expo2(absx);
+       t = __expo2(absx, 2*h);
        return t;
 }

diff --git a/src/math/sinhf.c b/src/math/sinhf.c
index 6ad19ea2b0c8d2d639d49feef32106a80b7f0235..b9caa793c4f927e10dff4d0237c78fe677c54f58 100644 (file)
--- a/src/math/sinhf.c
+++ b/src/math/sinhf.c
@@ -26,6 +26,6 @@ float sinhf(float x)
        }
 
        /* |x| > logf(FLT_MAX) or nan */
-       t = 2*h*__expo2f(absx);
+       t = __expo2f(absx, 2*h);
        return t;
 }