aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--math/exp.c23
-rw-r--r--math/exp2.c33
-rw-r--r--math/exp_data.c1
-rw-r--r--math/pow.c25
-rw-r--r--test/testcases/directed/exp.tst1
-rw-r--r--test/testcases/directed/pow.tst1
6 files changed, 48 insertions, 36 deletions
diff --git a/math/exp.c b/math/exp.c
index 9ee17e2..cc72780 100644
--- a/math/exp.c
+++ b/math/exp.c
@@ -40,10 +40,10 @@ specialcase (double_t tmp, uint64_t sbits, uint64_t ki)
if ((ki & 0x80000000) == 0)
{
- /* k > 0, the exponent of scale might have overflowed by <= 85. */
- sbits -= 97ull << 52;
+ /* k > 0, the exponent of scale might have overflowed by <= 460. */
+ sbits -= 1009ull << 52;
scale = asdouble (sbits);
- y = 0x1p97 * (scale + scale * tmp);
+ y = 0x1p1009 * (scale + scale * tmp);
return check_oflow (y);
}
/* k < 0, need special care in the subnormal range. */
@@ -72,9 +72,9 @@ specialcase (double_t tmp, uint64_t sbits, uint64_t ki)
}
static inline uint32_t
-top16 (double x)
+top12 (double x)
{
- return asuint64 (x) >> 48;
+ return asuint64 (x) >> 52;
}
static inline double
@@ -85,18 +85,18 @@ exp_inline (double x, double xtail, int hastail)
/* double_t for better performance on targets with FLT_EVAL_METHOD==2. */
double_t kd, z, r, r2, scale, tail, tmp;
- abstop = top16 (x) & 0x7fff;
- if (unlikely (abstop - top16 (0x1p-54) >= top16 (704.0) - top16 (0x1p-54)))
+ abstop = top12 (x) & 0x7ff;
+ if (unlikely (abstop - top12 (0x1p-54) >= top12 (512.0) - top12 (0x1p-54)))
{
- if (abstop - top16 (0x1p-54) >= 0x80000000)
+ if (abstop - top12 (0x1p-54) >= 0x80000000)
/* Avoid spurious underflow for tiny x. */
/* Note: 0 is common input. */
return WANT_ROUNDING ? 1.0 + x : 1.0;
- if (abstop >= top16 (768.0))
+ if (abstop >= top12 (1024.0))
{
if (asuint64 (x) == asuint64 (-INFINITY))
return 0.0;
- if (abstop >= top16 (INFINITY))
+ if (abstop >= top12 (INFINITY))
return 1.0 + x;
if (asuint64 (x) >> 63)
return __math_uflow (0);
@@ -125,6 +125,7 @@ exp_inline (double x, double xtail, int hastail)
kd -= Shift;
#endif
r = x + kd*NegLn2hiN + kd*NegLn2loN;
+ /* The code assumes 2^-200 < |xtail| < 2^-8/N. */
if (hastail)
r += xtail;
/* 2^(k/N) ~= scale * (1 + tail). */
@@ -148,6 +149,8 @@ exp_inline (double x, double xtail, int hastail)
if (unlikely (abstop == 0))
return specialcase (tmp, sbits, ki);
scale = asdouble (sbits);
+ /* Note: tmp == 0 or |tmp| > 2^-200 and scale > 2^-739, so there
+ is no spurious underflow here even without fma. */
return scale + scale * tmp;
}
diff --git a/math/exp2.c b/math/exp2.c
index b11b1f0..5acabf2 100644
--- a/math/exp2.c
+++ b/math/exp2.c
@@ -38,10 +38,10 @@ specialcase (double_t tmp, uint64_t sbits, uint64_t ki)
if ((ki & 0x80000000) == 0)
{
- /* k > 0, the exponent of scale might have overflowed by <= 85. */
- sbits -= 97ull << 52;
+ /* k > 0, the exponent of scale might have overflowed by 1. */
+ sbits -= 1ull << 52;
scale = asdouble (sbits);
- y = 0x1p97 * (scale + scale * tmp);
+ y = 2 * (scale + scale * tmp);
return check_oflow (y);
}
/* k < 0, need special care in the subnormal range. */
@@ -70,9 +70,9 @@ specialcase (double_t tmp, uint64_t sbits, uint64_t ki)
}
static inline uint32_t
-top16 (double x)
+top12 (double x)
{
- return asuint64 (x) >> 48;
+ return asuint64 (x) >> 52;
}
double
@@ -83,26 +83,27 @@ exp2 (double x)
/* double_t for better performance on targets with FLT_EVAL_METHOD==2. */
double_t kd, r, r2, scale, tail, tmp;
- abstop = top16 (x) & 0x7fff;
- if (unlikely (abstop - top16 (0x1p-54) >= top16 (992.0) - top16 (0x1p-54)))
+ abstop = top12 (x) & 0x7ff;
+ if (unlikely (abstop - top12 (0x1p-54) >= top12 (512.0) - top12 (0x1p-54)))
{
- if (abstop - top16 (0x1p-54) >= 0x80000000)
+ if (abstop - top12 (0x1p-54) >= 0x80000000)
/* Avoid spurious underflow for tiny x. */
/* Note: 0 is common input. */
return WANT_ROUNDING ? 1.0 + x : 1.0;
- if (abstop >= top16 (1088.0))
+ if (abstop >= top12 (1024.0))
{
if (asuint64 (x) == asuint64 (-INFINITY))
return 0.0;
- if (abstop >= top16 (INFINITY))
+ if (abstop >= top12 (INFINITY))
return 1.0 + x;
- if (asuint64 (x) >> 63)
- return __math_uflow (0);
- else
+ if (!(asuint64 (x) >> 63))
return __math_oflow (0);
+ else if (asuint64 (x) >= asuint64 (-1075.0))
+ return __math_uflow (0);
}
- /* Large x is special cased below. */
- abstop = 0;
+ if (2 * asuint64 (x) > 2 * asuint64 (928.0))
+ /* Large x is special cased below. */
+ abstop = 0;
}
/* exp2(x) = 2^(k/N) * 2^r, with 2^r in [2^(-1/2N),2^(1/2N)]. */
@@ -132,6 +133,8 @@ exp2 (double x)
if (unlikely (abstop == 0))
return specialcase (tmp, sbits, ki);
scale = asdouble (sbits);
+ /* Note: tmp == 0 or |tmp| > 2^-65 and scale > 2^-928, so there
+ is no spurious underflow here even without fma. */
return scale + scale * tmp;
}
#if USE_GLIBC_ABI
diff --git a/math/exp_data.c b/math/exp_data.c
index adfca9d..2961a94 100644
--- a/math/exp_data.c
+++ b/math/exp_data.c
@@ -67,6 +67,7 @@ const struct exp_data __exp_data = {
// abs error: 1.555*2^-66
// ulp error: 0.509 (0.511 without fma)
// if |x| < ln2/256+eps
+// abs error if |x| < ln2/256+0x1p-15: 1.09*2^-65
// abs error if |x| < ln2/128: 1.7145*2^-56
0x1.ffffffffffdbdp-2,
0x1.555555555543cp-3,
diff --git a/math/pow.c b/math/pow.c
index ba50153..e118737 100644
--- a/math/pow.c
+++ b/math/pow.c
@@ -36,9 +36,9 @@ ulperr_exp: 0.509 ULP (ULP error of exp)
#define OFF 0x3fe6000000000000
static inline uint32_t
-top16 (double x)
+top12 (double x)
{
- return asuint64 (x) >> 48;
+ return asuint64 (x) >> 52;
}
static inline double_t
@@ -142,10 +142,10 @@ specialcase (double_t tmp, uint64_t sbits, uint64_t ki)
if ((ki & 0x80000000) == 0)
{
- /* k > 0, the exponent of scale might have overflowed by <= 85. */
- sbits -= 97ull << 52;
+ /* k > 0, the exponent of scale might have overflowed by <= 460. */
+ sbits -= 1009ull << 52;
scale = asdouble (sbits);
- y = 0x1p97 * (scale + scale * tmp);
+ y = 0x1p1009 * (scale + scale * tmp);
return check_oflow (y);
}
/* k < 0, need special care in the subnormal range. */
@@ -186,17 +186,17 @@ exp_inline (double x, double xtail, uint32_t sign_bias)
/* double_t for better performance on targets with FLT_EVAL_METHOD==2. */
double_t kd, z, r, r2, scale, tail, tmp;
- abstop = top16 (x) & 0x7fff;
- if (unlikely (abstop - top16 (0x1p-54) >= top16 (704.0) - top16 (0x1p-54)))
+ abstop = top12 (x) & 0x7ff;
+ if (unlikely (abstop - top12 (0x1p-54) >= top12 (512.0) - top12 (0x1p-54)))
{
- if (abstop - top16 (0x1p-54) >= 0x80000000)
+ if (abstop - top12 (0x1p-54) >= 0x80000000)
{
/* Avoid spurious underflow for tiny x. */
/* Note: 0 is common input. */
double_t one = WANT_ROUNDING ? 1.0 + x : 1.0;
return sign_bias ? -one : one;
}
- if (abstop >= top16 (768.0))
+ if (abstop >= top12 (1024.0))
{
/* Note: inf and nan are already handled. */
if (asuint64 (x) >> 63)
@@ -226,6 +226,7 @@ exp_inline (double x, double xtail, uint32_t sign_bias)
kd -= Shift;
#endif
r = x + kd*NegLn2hiN + kd*NegLn2loN;
+ /* The code assumes 2^-200 < |xtail| < 2^-8/N. */
r += xtail;
/* 2^(k/N) ~= scale * (1 + tail). */
idx = 2*(ki % N);
@@ -248,6 +249,8 @@ exp_inline (double x, double xtail, uint32_t sign_bias)
if (unlikely (abstop == 0))
return specialcase (tmp, sbits, ki);
scale = asdouble (sbits);
+ /* Note: tmp == 0 or |tmp| > 2^-200 and scale > 2^-739, so there
+ is no spurious underflow here even without fma. */
return scale + scale * tmp;
}
@@ -282,8 +285,8 @@ pow (double x, double y)
ix = asuint64 (x);
iy = asuint64 (y);
- topx = ix >> 52;
- topy = iy >> 52;
+ topx = top12 (x);
+ topy = top12 (y);
if (unlikely (topx - 0x001 >= 0x7ff - 0x001 || (topy & 0x7ff) - 0x3be >= 0x43e - 0x3be))
{
/* Note: if |y| > 1075 * ln2 * 2^53 ~= 0x1.749p62 then pow(x,y) = inf/0
diff --git a/test/testcases/directed/exp.tst b/test/testcases/directed/exp.tst
index 570669e..b0e1c96 100644
--- a/test/testcases/directed/exp.tst
+++ b/test/testcases/directed/exp.tst
@@ -40,3 +40,4 @@ func=exp op1=40862e42.fefa39ef result=7fefffff.ffffff2a.1b1 errno=0
func=exp op1=40862e42.fefa39f0 result=7ff00000.00000000 errno=ERANGE status=ox
func=exp op1=c0874910.d52d3051 result=00000000.00000001 status=ux
func=exp op1=c0874910.d52d3052 result=00000000.00000000 errno=ERANGE status=ux
+func=exp op1=c085d589.f2fe5107 result=00f00000.000000f1.46b errno=0
diff --git a/test/testcases/directed/pow.tst b/test/testcases/directed/pow.tst
index 52c6991..3aa756c 100644
--- a/test/testcases/directed/pow.tst
+++ b/test/testcases/directed/pow.tst
@@ -448,6 +448,7 @@ func=pow op1=40000000.00000000 op2=c0080000.00000000 result=3fc00000.00000000 er
func=pow op1=40000000.00000000 op2=c0120000.00000000 result=3fa6a09e.667f3bcc.908 errno=0
func=pow op1=40000000.00000000 op2=c0180000.00000000 result=3f900000.00000000 errno=0
func=pow op1=40000000.00000000 op2=c07f3000.00000000 result=20c00000.00000000 errno=0
+func=pow op1=40000000.00000000 op2=c08f3a00.00000000 result=017ae89f.995ad3ad.5e8 errno=0
func=pow op1=40000000.00000000 op2=c090ce00.00000000 result=00000000.00000000.5a8 errno=ERANGE status=u
func=pow op1=40000000.00000000 op2=c3dfffff.ffffffff result=00000000.00000000 errno=ERANGE status=u
func=pow op1=40000000.00000000 op2=c3e00000.00000000 result=00000000.00000000 errno=ERANGE status=u
generated by cgit v1.2.3 (git 2.25.1) at 2024-07-08 14:22:02 +0000