提交 f4d9bfb3 编写于 作者: S Szabolcs Nagy

math: fix atanh (overflow and underflow issues)

in atanh exception handling was left to the called log functions,
but the argument to those functions could underflow or overflow.

use double_t and float_t to avoid some useless stores on x86
上级 afa2aacc
......@@ -6,16 +6,24 @@ double atanh(double x)
union {double f; uint64_t i;} u = {.f = x};
unsigned e = u.i >> 52 & 0x7ff;
unsigned s = u.i >> 63;
double_t y;
/* |x| */
u.i &= (uint64_t)-1/2;
x = u.f;
y = u.f;
if (e < 0x3ff - 1) {
/* |x| < 0.5, up to 1.7ulp error */
x = 0.5*log1p(2*x + 2*x*x/(1-x));
if (e < 0x3ff - 32) {
/* handle underflow */
if (e == 0)
FORCE_EVAL((float)y);
} else {
/* |x| < 0.5, up to 1.7ulp error */
y = 0.5*log1p(2*y + 2*y*y/(1-y));
}
} else {
x = 0.5*log1p(2*x/(1-x));
/* avoid overflow */
y = 0.5*log1p(2*(y/(1-y)));
}
return s ? -x : x;
return s ? -y : y;
}
......@@ -5,16 +5,24 @@ float atanhf(float x)
{
union {float f; uint32_t i;} u = {.f = x};
unsigned s = u.i >> 31;
float_t y;
/* |x| */
u.i &= 0x7fffffff;
x = u.f;
y = u.f;
if (u.i < 0x3f800000 - (1<<23)) {
/* |x| < 0.5, up to 1.7ulp error */
x = 0.5f*log1pf(2*x + 2*x*x/(1-x));
if (u.i < 0x3f800000 - (32<<23)) {
/* handle underflow */
if (u.i < (1<<23))
FORCE_EVAL((float)(y*y));
} else {
/* |x| < 0.5, up to 1.7ulp error */
y = 0.5f*log1pf(2*y + 2*y*y/(1-y));
}
} else {
x = 0.5f*log1pf(2*x/(1-x));
/* avoid overflow */
y = 0.5f*log1pf(2*(y/(1-y)));
}
return s ? -x : x;
return s ? -y : y;
}
......@@ -17,11 +17,18 @@ long double atanhl(long double x)
u.i.se = e;
x = u.f;
if (e < 0x3fff - 1) {
/* |x| < 0.5, up to 1.7ulp error */
x = 0.5*log1pl(2*x + 2*x*x/(1-x));
if (e < 0x3ff - 1) {
if (e < 0x3ff - LDBL_MANT_DIG/2) {
/* handle underflow */
if (e == 0)
FORCE_EVAL((float)x);
} else {
/* |x| < 0.5, up to 1.7ulp error */
x = 0.5*log1pl(2*x + 2*x*x/(1-x));
}
} else {
x = 0.5*log1pl(2*x/(1-x));
/* avoid overflow */
x = 0.5*log1pl(2*(x/(1-x)));
}
return s ? -x : x;
}
......
Markdown is supported
0% .
You are about to add 0 people to the discussion. Proceed with caution.
先完成此消息的编辑!
想要评论请 注册