summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authornsz <nsz@port70.net>2012-03-27 22:49:37 +0200
committernsz <nsz@port70.net>2012-03-27 22:49:37 +0200
commitcf682072ce16080ebe199e5055f686d81c8416ce (patch)
tree4a860dae892a57c8c2cf75d82371a915ebf3af3b
parentbbfbc7edaf992abe1d3d09868be07c4c1cc44db7 (diff)
downloadmusl-cf682072ce16080ebe199e5055f686d81c8416ce.tar.gz
musl-cf682072ce16080ebe199e5055f686d81c8416ce.tar.bz2
musl-cf682072ce16080ebe199e5055f686d81c8416ce.tar.xz
musl-cf682072ce16080ebe199e5055f686d81c8416ce.zip
math: fix a regression in powl and do some cleanups
previously a division was accidentally turned into integer div (w = -i/NXT;) instead of long double div (w = -i; w /= NXT;)
-rw-r--r--src/math/powl.c23
1 files changed, 11 insertions, 12 deletions
diff --git a/src/math/powl.c b/src/math/powl.c
index 2ee0b10b..ce6274cf 100644
--- a/src/math/powl.c
+++ b/src/math/powl.c
@@ -165,8 +165,6 @@ static const long double R[] = {
6.9314718055994530931447E-1L,
};
-#define douba(k) A[k]
-#define doubb(k) B[k]
#define MEXP (NXT*16384.0L)
/* The following if denormal numbers are supported, else -MEXP: */
#define MNEXP (-NXT*(16384.0L+64.0L))
@@ -300,15 +298,15 @@ long double powl(long double x, long double y)
/* find significand in antilog table A[] */
i = 1;
- if (x <= douba(17))
+ if (x <= A[17])
i = 17;
- if (x <= douba(i+8))
+ if (x <= A[i+8])
i += 8;
- if (x <= douba(i+4))
+ if (x <= A[i+4])
i += 4;
- if (x <= douba(i+2))
+ if (x <= A[i+2])
i += 2;
- if (x >= douba(1))
+ if (x >= A[1])
i = -1;
i += 1;
@@ -319,9 +317,9 @@ long double powl(long double x, long double y)
*
* log(x/a) = log(1+v), v = x/a - 1 = (x-a)/a
*/
- x -= douba(i);
- x -= doubb(i/2);
- x /= douba(i);
+ x -= A[i];
+ x -= B[i/2];
+ x /= A[i];
/* rational approximation for log(1+v):
*
@@ -340,7 +338,8 @@ long double powl(long double x, long double y)
z += x;
/* Compute exponent term of the base 2 logarithm. */
- w = -i / NXT;
+ w = -i;
+ w /= NXT;
w += e;
/* Now base 2 log of x is w + z. */
@@ -397,7 +396,7 @@ long double powl(long double x, long double y)
i = 1;
i = e/NXT + i;
e = NXT*i - e;
- w = douba(e);
+ w = A[e];
z = w * z; /* 2**-e * ( 1 + (2**Hb-1) ) */
z = z + w;
z = scalbnl(z, i); /* multiply by integer power of 2 */