summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorSzabolcs Nagy <nsz@port70.net>2020-01-20 20:38:45 +0000
committerRich Felker <dalias@aerifal.cx>2020-02-21 23:42:12 -0500
commitd20558148d8a2c52229b02668627697e83ca3840 (patch)
treee4f3373bb8649db621bfeb871ed335e13886b343
parentb3797d3b2e10e6fff2a6b04af917e61e95838b08 (diff)
downloadmusl-d20558148d8a2c52229b02668627697e83ca3840.tar.gz
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.
-rw-r--r--src/internal/libm.h4
-rw-r--r--src/math/__expo2.c5
-rw-r--r--src/math/__expo2f.c5
-rw-r--r--src/math/cosh.c2
-rw-r--r--src/math/coshf.c2
-rw-r--r--src/math/sinh.c2
-rw-r--r--src/math/sinhf.c2
7 files changed, 12 insertions, 10 deletions
diff --git a/src/internal/libm.h b/src/internal/libm.h
index b5bd26b8..7533f6ba 100644
--- 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 740ac680..248f052b 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 5163e418..538eb09c 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 100f8231..490c15fb 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 b09f2ee5..e739cff9 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 00022c4e..a01951ae 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 6ad19ea2..b9caa793 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;
}