diff options
author | Szabolcs Nagy <nsz@port70.net> | 2020-01-20 20:38:45 +0000 |
---|---|---|
committer | Rich Felker <dalias@aerifal.cx> | 2020-02-21 23:42:12 -0500 |
commit | d20558148d8a2c52229b02668627697e83ca3840 (patch) | |
tree | e4f3373bb8649db621bfeb871ed335e13886b343 /src/math | |
parent | b3797d3b2e10e6fff2a6b04af917e61e95838b08 (diff) | |
download | musl-d20558148d8a2c52229b02668627697e83ca3840.zip musl-d20558148d8a2c52229b02668627697e83ca3840.tar.gz musl-d20558148d8a2c52229b02668627697e83ca3840.tar.bz2 |
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.
Diffstat (limited to 'src/math')
-rw-r--r-- | src/math/__expo2.c | 5 | ||||
-rw-r--r-- | src/math/__expo2f.c | 5 | ||||
-rw-r--r-- | src/math/cosh.c | 2 | ||||
-rw-r--r-- | src/math/coshf.c | 2 | ||||
-rw-r--r-- | src/math/sinh.c | 2 | ||||
-rw-r--r-- | src/math/sinhf.c | 2 |
6 files changed, 10 insertions, 8 deletions
diff --git a/src/math/__expo2.c b/src/math/__expo2.c index 740ac68..248f052 100644 --- 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 5163e41..538eb09 100644 --- 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 100f823..490c15f 100644 --- 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 b09f2ee..e739cff 100644 --- 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 00022c4..a01951a 100644 --- 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 6ad19ea..b9caa79 100644 --- 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; } |