use alternate formula for acos asm to avoid loss of precision

[musl] / src / math / powl.c

diff --git a/src/math/powl.c b/src/math/powl.c
index 690f294..a1d2f07 100644 (file)
--- a/src/math/powl.c
+++ b/src/math/powl.c
@@ -84,13 +84,13 @@ long double powl(long double x, long double y)
 /* log(1+x) =  x - .5x^2 + x^3 *  P(z)/Q(z)
  * on the domain  2^(-1/32) - 1  <=  x  <=  2^(1/32) - 1
  */
-static long double P[] = {
+static const long double P[] = {
  8.3319510773868690346226E-4L,
  4.9000050881978028599627E-1L,
  1.7500123722550302671919E0L,
  1.4000100839971580279335E0L,
 };
-static long double Q[] = {
+static const long double Q[] = {
 /* 1.0000000000000000000000E0L,*/
  5.2500282295834889175431E0L,
  8.4000598057587009834666E0L,
@@ -99,7 +99,7 @@ static long double Q[] = {
 /* A[i] = 2^(-i/32), rounded to IEEE long double precision.
  * If i is even, A[i] + B[i/2] gives additional accuracy.
  */
-static long double A[33] = {
+static const long double A[33] = {
  1.0000000000000000000000E0L,
  9.7857206208770013448287E-1L,
  9.5760328069857364691013E-1L,
@@ -134,7 +134,7 @@ static long double A[33] = {
  5.1094857432705833910408E-1L,
  5.0000000000000000000000E-1L,
 };
-static long double B[17] = {
+static const long double B[17] = {
  0.0000000000000000000000E0L,
  2.6176170809902549338711E-20L,
 -1.0126791927256478897086E-20L,
@@ -157,7 +157,7 @@ static long double B[17] = {
 /* 2^x = 1 + x P(x),
  * on the interval -1/32 <= x <= 0
  */
-static long double R[] = {
+static const long double R[] = {
  1.5089970579127659901157E-5L,
  1.5402715328927013076125E-4L,
  1.3333556028915671091390E-3L,
@@ -188,11 +188,9 @@ static long double R[] = {
 static const long double MAXLOGL = 1.1356523406294143949492E4L;
 static const long double MINLOGL = -1.13994985314888605586758E4L;
 static const long double LOGE2L = 6.9314718055994530941723E-1L;
-static volatile long double z;
-static long double w, W, Wa, Wb, ya, yb, u;
 static const long double huge = 0x1p10000L;
 /* XXX Prevent gcc from erroneously constant folding this. */
-static volatile long double twom10000 = 0x1p-10000L;
+static const volatile long double twom10000 = 0x1p-10000L;
 
 static long double reducl(long double);
 static long double powil(long double, int);
@@ -202,45 +200,47 @@ long double powl(long double x, long double y)
        /* double F, Fa, Fb, G, Ga, Gb, H, Ha, Hb */
        int i, nflg, iyflg, yoddint;
        long e;
+       volatile long double z=0;
+       long double w=0, W=0, Wa=0, Wb=0, ya=0, yb=0, u=0;
 
-       if (y == 0.0L)
-               return 1.0L;
+       if (y == 0.0)
+               return 1.0;
        if (isnan(x))
                return x;
        if (isnan(y))
                return y;
-       if (y == 1.0L)
+       if (y == 1.0)
                return x;
 
        // FIXME: this is wrong, see pow special cases in c99 F.9.4.4
-       if (!isfinite(y) && (x == -1.0L || x == 1.0L) )
+       if (!isfinite(y) && (x == -1.0 || x == 1.0) )
                return y - y;   /* +-1**inf is NaN */
-       if (x == 1.0L)
-               return 1.0L;
+       if (x == 1.0)
+               return 1.0;
        if (y >= LDBL_MAX) {
-               if (x > 1.0L)
+               if (x > 1.0)
                        return INFINITY;
-               if (x > 0.0L && x < 1.0L)
-                       return 0.0L;
-               if (x < -1.0L)
+               if (x > 0.0 && x < 1.0)
+                       return 0.0;
+               if (x < -1.0)
                        return INFINITY;
-               if (x > -1.0L && x < 0.0L)
-                       return 0.0L;
+               if (x > -1.0 && x < 0.0)
+                       return 0.0;
        }
        if (y <= -LDBL_MAX) {
-               if (x > 1.0L)
-                       return 0.0L;
-               if (x > 0.0L && x < 1.0L)
+               if (x > 1.0)
+                       return 0.0;
+               if (x > 0.0 && x < 1.0)
                        return INFINITY;
-               if (x < -1.0L)
-                       return 0.0L;
-               if (x > -1.0L && x < 0.0L)
+               if (x < -1.0)
+                       return 0.0;
+               if (x > -1.0 && x < 0.0)
                        return INFINITY;
        }
        if (x >= LDBL_MAX) {
-               if (y > 0.0L)
+               if (y > 0.0)
                        return INFINITY;
-               return 0.0L;
+               return 0.0;
        }
 
        w = floorl(y);
@@ -253,29 +253,29 @@ long double powl(long double x, long double y)
        yoddint = 0;
        if (iyflg) {
                ya = fabsl(y);
-               ya = floorl(0.5L * ya);
-               yb = 0.5L * fabsl(w);
+               ya = floorl(0.5 * ya);
+               yb = 0.5 * fabsl(w);
                if( ya != yb )
                        yoddint = 1;
        }
 
        if (x <= -LDBL_MAX) {
-               if (y > 0.0L) {
+               if (y > 0.0) {
                        if (yoddint)
                                return -INFINITY;
                        return INFINITY;
                }
-               if (y < 0.0L) {
+               if (y < 0.0) {
                        if (yoddint)
-                               return -0.0L;
+                               return -0.0;
                        return 0.0;
                }
        }
 
 
        nflg = 0;       /* flag = 1 if x<0 raised to integer power */
-       if (x <= 0.0L) {
-               if (x == 0.0L) {
+       if (x <= 0.0) {
+               if (x == 0.0) {
                        if (y < 0.0) {
                                if (signbit(x) && yoddint)
                                        return -INFINITY;
@@ -283,12 +283,12 @@ long double powl(long double x, long double y)
                        }
                        if (y > 0.0) {
                                if (signbit(x) && yoddint)
-                                       return -0.0L;
+                                       return -0.0;
                                return 0.0;
                        }
-                       if (y == 0.0L)
-                               return 1.0L;  /*   0**0   */
-                       return 0.0L;  /*   0**y   */
+                       if (y == 0.0)
+                               return 1.0;  /*   0**0   */
+                       return 0.0;  /*   0**y   */
                }
                if (iyflg == 0)
                        return (x - x) / (x - x); /* (x<0)**(non-int) is NaN */
@@ -343,7 +343,7 @@ long double powl(long double x, long double y)
         */
        z = x*x;
        w = x * (z * __polevll(x, P, 3) / __p1evll(x, Q, 3));
-       w = w - ldexpl(z, -1);  /*  w - 0.5 * z  */
+       w = w - 0.5*z;
 
        /* Convert to base 2 logarithm:
         * multiply by log2(e) = 1 + LOG2EA
@@ -355,7 +355,8 @@ long double powl(long double x, long double y)
 
        /* Compute exponent term of the base 2 logarithm. */
        w = -i;
-       w = ldexpl(w, -LNXT); /* divide by NXT */
+       // TODO: use w * 0x1p-5;
+       w = scalbnl(w, -LNXT); /* divide by NXT */
        w += e;
        /* Now base 2 log of x is w + z. */
 
@@ -380,7 +381,7 @@ long double powl(long double x, long double y)
 
        H = Fb + Gb;
        Ha = reducl(H);
-       w = ldexpl( Ga+Ha, LNXT );
+       w = scalbnl( Ga+Ha, LNXT );
 
        /* Test the power of 2 for overflow */
        if (w > MEXP)
@@ -391,9 +392,9 @@ long double powl(long double x, long double y)
        e = w;
        Hb = H - Ha;
 
-       if (Hb > 0.0L) {
+       if (Hb > 0.0) {
                e += 1;
-               Hb -= 1.0L/NXT;  /*0.0625L;*/
+               Hb -= 1.0/NXT;  /*0.0625L;*/
        }
 
        /* Now the product y * log2(x)  =  Hb + e/NXT.
@@ -415,16 +416,16 @@ long double powl(long double x, long double y)
        w = douba(e);
        z = w * z;  /*  2**-e * ( 1 + (2**Hb-1) )  */
        z = z + w;
-       z = ldexpl(z, i);  /* multiply by integer power of 2 */
+       z = scalbnl(z, i);  /* multiply by integer power of 2 */
 
        if (nflg) {
                /* For negative x,
                 * find out if the integer exponent
                 * is odd or even.
                 */
-               w = ldexpl(y, -1);
+               w = 0.5*y;
                w = floorl(w);
-               w = ldexpl(w, 1);
+               w = 2.0*w;
                if (w != y)
                        z = -z;  /* odd exponent */
        }
@@ -438,9 +439,9 @@ static long double reducl(long double x)
 {
        long double t;
 
-       t = ldexpl(x, LNXT);
+       t = scalbnl(x, LNXT);
        t = floorl(t);
-       t = ldexpl(t, -LNXT);
+       t = scalbnl(t, -LNXT);
        return t;
 }
 
@@ -483,18 +484,18 @@ static long double powil(long double x, int nn)
        long double s;
        int n, e, sign, asign, lx;
 
-       if (x == 0.0L) {
+       if (x == 0.0) {
                if (nn == 0)
-                       return 1.0L;
+                       return 1.0;
                else if (nn < 0)
                        return LDBL_MAX;
-               return 0.0L;
+               return 0.0;
        }
 
        if (nn == 0)
-               return 1.0L;
+               return 1.0;
 
-       if (x < 0.0L) {
+       if (x < 0.0) {
                asign = -1;
                x = -x;
        } else
@@ -516,7 +517,7 @@ static long double powil(long double x, int nn)
        e = (lx - 1)*n;
        if ((e == 0) || (e > 64) || (e < -64)) {
                s = (s - 7.0710678118654752e-1L) / (s +  7.0710678118654752e-1L);
-               s = (2.9142135623730950L * s - 0.5L + lx) * nn * LOGE2L;
+               s = (2.9142135623730950L * s - 0.5 + lx) * nn * LOGE2L;
        } else {
                s = LOGE2L * e;
        }
@@ -530,8 +531,8 @@ static long double powil(long double x, int nn)
         * since roundoff error in 1.0/x will be amplified.
         * The precise demarcation should be the gradual underflow threshold.
         */
-       if (s < -MAXLOGL+2.0L) {
-               x = 1.0L/x;
+       if (s < -MAXLOGL+2.0) {
+               x = 1.0/x;
                sign = -sign;
        }
 
@@ -539,7 +540,7 @@ static long double powil(long double x, int nn)
        if (n & 1)
                y = x;
        else {
-               y = 1.0L;
+               y = 1.0;
                asign = 0;
        }
 
@@ -555,7 +556,7 @@ static long double powil(long double x, int nn)
        if (asign)
                y = -y;  /* odd power of negative number */
        if (sign < 0)
-               y = 1.0L/y;
+               y = 1.0/y;
        return y;
 }