use scalbn or *2.0 instead of ldexp, fix fmal

Some code assumed ldexp(x, 1) is faster than 2.0*x,
but ldexp is a wrapper around scalbn which uses
multiplications inside, so this optimization is
wrong.

This commit also fixes fmal which accidentally
used ldexp instead of ldexpl loosing precision.

There are various additional changes from the
work-in-progress const cleanups.

diff --git a/src/math/expl.c b/src/math/expl.c
index 9507fd2e32f1b3779f35675b597ba510d5be4928..b289ffece062140d074c94396d5e8cb62ea2c535 100644 (file)
--- a/src/math/expl.c
+++ b/src/math/expl.c
@@ -102,13 +102,13 @@ long double expl(long double x)
        if (x > MAXLOGL)
                return INFINITY;
        if (x < MINLOGL)
-               return 0.0L;
+               return 0.0;
 
        /* Express e**x = e**g 2**n
         *   = e**g e**(n loge(2))
         *   = e**(g + n loge(2))
         */
-       px = floorl(LOG2EL * x + 0.5L); /* floor() truncates toward -infinity. */
+       px = floorl(LOG2EL * x + 0.5); /* floor() truncates toward -infinity. */
        n = px;
        x -= px * C1;
        x -= px * C2;
@@ -120,8 +120,8 @@ long double expl(long double x)
        xx = x * x;
        px = x * __polevll(xx, P, 2);
        x =  px/(__polevll(xx, Q, 3) - px);
-       x = 1.0L + ldexpl(x, 1);
-       x = ldexpl(x, n);
+       x = 1.0 + 2.0 * x;
+       x = scalbnl(x, n);
        return x;
 }
 #endif

diff --git a/src/math/expm1l.c b/src/math/expm1l.c
index 2f94dfa20071dd115c1d699f7fed1c7dfffa16f3..ca0d141f44f6840d7ecba899c667a4c757762474 100644 (file)
--- a/src/math/expm1l.c
+++ b/src/math/expm1l.c
@@ -97,11 +97,11 @@ long double expm1l(long double x)
                return x;
        /* Minimum value.*/
        if (x < minarg)
-               return -1.0L;
+               return -1.0;
 
        xx = C1 + C2;
        /* Express x = ln 2 (k + remainder), remainder not exceeding 1/2. */
-       px = floorl (0.5 + x / xx);
+       px = floorl(0.5 + x / xx);
        k = px;
        /* remainder times ln 2 */
        x -= px * C1;
@@ -116,7 +116,7 @@ long double expm1l(long double x)
        /* exp(x) = exp(k ln 2) exp(remainder ln 2) = 2^k exp(remainder ln 2).
         We have qx = exp(remainder ln 2) - 1, so
         exp(x) - 1  =  2^k (qx + 1) - 1  =  2^k qx + 2^k - 1.  */
-       px = ldexpl(1.0L, k);
+       px = scalbnl(1.0, k);
        x = px * qx + (px - 1.0);
        return x;
 }

diff --git a/src/math/fma.c b/src/math/fma.c
index 87d450c731ac55d95bd95cbb032c222ec8139f2e..5fb9540601c520dac6ce9a519e8c9abe91aed7d7 100644 (file)
--- a/src/math/fma.c
+++ b/src/math/fma.c
@@ -247,7 +247,7 @@ static inline double add_and_denormalize(double a, double b, int scale)
                        INSERT_WORD64(sum.hi, hibits);
                }
        }
-       return (ldexp(sum.hi, scale));
+       return scalbn(sum.hi, scale);
 }
 
 /*
@@ -364,7 +364,7 @@ double fma(double x, double y, double z)
                }
        }
        if (spread <= DBL_MANT_DIG * 2)
-               zs = ldexp(zs, -spread);
+               zs = scalbn(zs, -spread);
        else
                zs = copysign(DBL_MIN, zs);
 
@@ -390,7 +390,7 @@ double fma(double x, double y, double z)
                 */
                fesetround(oround);
                volatile double vzs = zs; /* XXX gcc CSE bug workaround */
-               return (xy.hi + vzs + ldexp(xy.lo, spread));
+               return xy.hi + vzs + scalbn(xy.lo, spread);
        }
 
        if (oround != FE_TONEAREST) {
@@ -400,13 +400,13 @@ double fma(double x, double y, double z)
                 */
                fesetround(oround);
                adj = r.lo + xy.lo;
-               return (ldexp(r.hi + adj, spread));
+               return scalbn(r.hi + adj, spread);
        }
 
        adj = add_adjusted(r.lo, xy.lo);
        if (spread + ilogb(r.hi) > -1023)
-               return (ldexp(r.hi + adj, spread));
+               return scalbn(r.hi + adj, spread);
        else
-               return (add_and_denormalize(r.hi, adj, spread));
+               return add_and_denormalize(r.hi, adj, spread);
 }
 #endif

diff --git a/src/math/fmal.c b/src/math/fmal.c
index cbaf46ebec926516f3abd6afa7fd76212ff6a4aa..be64f145bbaeee42d5f29b63ce87f31251c15942 100644 (file)
--- a/src/math/fmal.c
+++ b/src/math/fmal.c
@@ -115,7 +115,7 @@ static inline long double add_and_denormalize(long double a, long double b, int
                if (bits_lost != 1 ^ (int)(u.bits.manl & 1))
                        sum.hi = nextafterl(sum.hi, INFINITY * sum.lo);
        }
-       return (ldexp(sum.hi, scale));
+       return scalbnl(sum.hi, scale);
 }
 
 /*
@@ -228,7 +228,7 @@ long double fmal(long double x, long double y, long double z)
                }
        }
        if (spread <= LDBL_MANT_DIG * 2)
-               zs = ldexpl(zs, -spread);
+               zs = scalbnl(zs, -spread);
        else
                zs = copysignl(LDBL_MIN, zs);
 
@@ -254,7 +254,7 @@ long double fmal(long double x, long double y, long double z)
                 */
                fesetround(oround);
                volatile long double vzs = zs; /* XXX gcc CSE bug workaround */
-               return (xy.hi + vzs + ldexpl(xy.lo, spread));
+               return xy.hi + vzs + scalbnl(xy.lo, spread);
        }
 
        if (oround != FE_TONEAREST) {
@@ -264,13 +264,13 @@ long double fmal(long double x, long double y, long double z)
                 */
                fesetround(oround);
                adj = r.lo + xy.lo;
-               return (ldexpl(r.hi + adj, spread));
+               return scalbnl(r.hi + adj, spread);
        }
 
        adj = add_adjusted(r.lo, xy.lo);
        if (spread + ilogbl(r.hi) > -16383)
-               return (ldexpl(r.hi + adj, spread));
+               return scalbnl(r.hi + adj, spread);
        else
-               return (add_and_denormalize(r.hi, adj, spread));
+               return add_and_denormalize(r.hi, adj, spread);
 }
 #endif

diff --git a/src/math/log10l.c b/src/math/log10l.c
index b954cc7773116f8c7ba192d595cde0aa103ff84c..f0eeeafb4a733a916d45d2b485886c35d5843973 100644 (file)
--- a/src/math/log10l.c
+++ b/src/math/log10l.c
@@ -123,9 +123,9 @@ long double log10l(long double x)
 
        if (isnan(x))
                return x;
-       if(x <= 0.0L) {
-               if(x == 0.0L)
-                       return -1.0L / (x - x);
+       if(x <= 0.0) {
+               if(x == 0.0)
+                       return -1.0 / (x - x);
                return (x - x) / (x - x);
        }
        if (x == INFINITY)
@@ -142,12 +142,12 @@ long double log10l(long double x)
        if (e > 2 || e < -2) {
                if (x < SQRTH) {  /* 2(2x-1)/(2x+1) */
                        e -= 1;
-                       z = x - 0.5L;
-                       y = 0.5L * z + 0.5L;
+                       z = x - 0.5;
+                       y = 0.5 * z + 0.5;
                } else {  /*  2 (x-1)/(x+1)   */
-                       z = x - 0.5L;
-                       z -= 0.5L;
-                       y = 0.5L * x  + 0.5L;
+                       z = x - 0.5;
+                       z -= 0.5;
+                       y = 0.5 * x  + 0.5;
                }
                x = z / y;
                z = x*x;
@@ -158,13 +158,13 @@ long double log10l(long double x)
        /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
        if (x < SQRTH) {
                e -= 1;
-               x = ldexpl(x, 1) - 1.0L; /*  2x - 1  */
+               x = 2.0*x - 1.0;
        } else {
-               x = x - 1.0L;
+               x = x - 1.0;
        }
        z = x*x;
        y = x * (z * __polevll(x, P, 6) / __p1evll(x, Q, 7));
-       y = y - ldexpl(z, -1);   /* -0.5x^2 + ... */
+       y = y - 0.5*z;
 
 done:
        /* Multiply log of fraction by log10(e)

diff --git a/src/math/log2l.c b/src/math/log2l.c
index 4339c033084a5af0b3cc5a207fce426015f1f58f..8ebce9c43be4b09a6010762bbc510866cd26a19c 100644 (file)
--- a/src/math/log2l.c
+++ b/src/math/log2l.c
@@ -121,8 +121,8 @@ long double log2l(long double x)
                return x;
        if (x == INFINITY)
                return x;
-       if (x <= 0.0L) {
-               if (x == 0.0L)
+       if (x <= 0.0) {
+               if (x == 0.0)
                        return -INFINITY;
                return NAN;
        }
@@ -139,12 +139,12 @@ long double log2l(long double x)
        if (e > 2 || e < -2) {
                if (x < SQRTH) {  /* 2(2x-1)/(2x+1) */
                        e -= 1;
-                       z = x - 0.5L;
-                       y = 0.5L * z + 0.5L;
+                       z = x - 0.5;
+                       y = 0.5 * z + 0.5;
                } else {  /*  2 (x-1)/(x+1)   */
-                       z = x - 0.5L;
-                       z -= 0.5L;
-                       y = 0.5L * x  + 0.5L;
+                       z = x - 0.5;
+                       z -= 0.5;
+                       y = 0.5 * x + 0.5;
                }
                x = z / y;
                z = x*x;
@@ -155,13 +155,13 @@ long double log2l(long double x)
        /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
        if (x < SQRTH) {
                e -= 1;
-               x = ldexpl(x, 1) - 1.0L; /*  2x - 1  */
+               x = 2.0*x - 1.0;
        } else {
-               x = x - 1.0L;
+               x = x - 1.0;
        }
        z = x*x;
        y = x * (z * __polevll(x, P, 6) / __p1evll(x, Q, 7));
-       y = y - ldexpl(z, -1);   /* -0.5x^2 + ... */
+       y = y - 0.5*z;
 
 done:
        /* Multiply log of fraction by log2(e)

diff --git a/src/math/logl.c b/src/math/logl.c
index ee7ca64a8201a44fdc1bbc04b522947bc9ad4cd5..ffd8103813eade25cc51b70944b93d78a49b8cf8 100644 (file)
--- a/src/math/logl.c
+++ b/src/math/logl.c
@@ -119,8 +119,8 @@ long double logl(long double x)
                return x;
        if (x == INFINITY)
                return x;
-       if (x <= 0.0L) {
-               if (x == 0.0L)
+       if (x <= 0.0) {
+               if (x == 0.0)
                        return -INFINITY;
                return NAN;
        }
@@ -137,12 +137,12 @@ long double logl(long double x)
        if (e > 2 || e < -2) {
                if (x < SQRTH) {  /* 2(2x-1)/(2x+1) */
                        e -= 1;
-                       z = x - 0.5L;
-                       y = 0.5L * z + 0.5L;
+                       z = x - 0.5;
+                       y = 0.5 * z + 0.5;
                } else {  /*  2 (x-1)/(x+1)   */
-                       z = x - 0.5L;
-                       z -= 0.5L;
-                       y = 0.5L * x  + 0.5L;
+                       z = x - 0.5;
+                       z -= 0.5;
+                       y = 0.5 * x  + 0.5;
                }
                x = z / y;
                z = x*x;
@@ -156,14 +156,14 @@ long double logl(long double x)
        /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
        if (x < SQRTH) {
                e -= 1;
-               x = ldexpl(x, 1) - 1.0L; /*  2x - 1  */
+               x = 2.0*x - 1.0;
        } else {
-               x = x - 1.0L;
+               x = x - 1.0;
        }
        z = x*x;
        y = x * (z * __polevll(x, P, 6) / __p1evll(x, Q, 6));
        y = y + e * C2;
-       z = y - ldexpl(z, -1);   /*  y - 0.5 * z  */
+       z = y - 0.5*z;
        /* Note, the sum of above terms does not exceed x/4,
         * so it contributes at most about 1/4 lsb to the error.
         */

diff --git a/src/math/powl.c b/src/math/powl.c
index a6ee275f25eedb4c95f367b48ce893009702f8bc..a1d2f0769eb550ac39f1cb36f93903c55ebf8e10 100644 (file)
--- a/src/math/powl.c
+++ b/src/math/powl.c
@@ -203,44 +203,44 @@ long double powl(long double x, long double y)
        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;
 }