[Issue 3749] cannot evaluate yl2x (log) and exp functions at compile time

[d-bugmail]

http://d.puremagic.com/issues/show_bug.cgi?id=3749



--- Comment #4 from Witold Baryluk <bary...@smp.if.uj.edu.pl> 2010-01-28 
15:55:35 PST ---
So I release this as public domain. I written this code as workaround to lack
of log and exp. They looks to be accurate to 16 digital digits on whole real
line.
Results can be slightly different than values returned by std.math.{log,exp}
unfortunetly. Anybody want to review this code or know better methods?

/** Calculate natural logarithm of x.
 *
 *
 * Performs reduction of large values to (0..3) inverval, using log(x 3^n) =
log(x) + n*log(3)
 *
 * Then uses truncated Taylor series in variable y=(x-1)/(x+1) for x > 0.
 *
 * For values x < 1, calculate -log(1/x)
 *
 */
double ctfe_log(double x_) {
    if (x_ == 0.0) {
        return -double.infinity;
    }
    if (x_ < 0.0) {
        return double.nan;
    }
    if (x_ == double.infinity) {
        return double.infinity;
    }

    if (x_ == 1.0) {
        return 0.0;
    }

    if (!(x_ == x_)) { // nan
        return double.nan;
    }

    real x = x_;

    if (x > 1.0) {
        uint m = 0;
        // reduce to (1 .. 3) interval
        while (x > 3.0) {
            x = x / 3.0;
            m = m + 1;
        }
        real y = (x-1.0)/(x+1.0);

        real y2 = y*y;

        /* Evaluate Horner's scheme on polynomial
         * log(x) = log((1+y) / (1-y)) = 2 y (1 + y^2/3 + y^4/5  + y^6/7 + ...
y^70/71)
         */
        real temp = 0.0;

        for (int i = 71; i >= 3; i -= 2) {
            temp += 1.0/cast(real)i;
            temp *= y2;
        }
        temp += 1.0;

        y = 2.0*y*temp;
        if (m) {
            return y + m*ctfe_log(3.0);
        } else {
            return y;
        }
    } else {
        return -ctfe_log(1.0/x);
    }
}

/** Compute exponential function of x.
 *
 * Uses truncated Taylor series expeansion of exp function.
 *
 * Performs reduction for |x| > 2, of the form exp(x 2^m) = exp(x)^(2^m)
 */
double ctfe_exp(double x_) {
    if (x_ == 0.0) {
        return 1.0;
    }

    if (x_ >= 710.0) { // includes +infinity
        return double.infinity;
    }
    if (x_ <= -746.0) { // includes -infinity
        return 0.0;
    }
    if (!(x_ == x_)) { // nan
        return double.nan;
    }

    real x = x_;

    int m = 0;

    // reduce to (-2 .. 2) interval
    if (x > 0.0) {
        while (x > 2.0) {
            x = x / 2.0;
            m = m + 1;
        }
    } else {
        while (x < -2.0) {
            x = x / 2.0;
            m = m - 1;
        }
    }


    real temp = 1.0;
    real term = 1.0;

    for (int i = 1; i <= 25; i++) {
        term *= x/cast(real)i;
        temp += term;
    }

    if (m) {
        real exp2 = ctfe_exp(2.0);
        if (m > 0) {
            for (int i = 0; i < m; i++) {
                temp = temp*temp;
            }
        } else {
            for (int i = 0; i < -m; i++) {
                temp = temp*temp;
            }
        }
        return temp;
    } else {
        return temp;
    }
}

int tests(double[] xs, int min, int max) {
    assert(min <= max);
    int r = 0;
    double c = 1.0;
    if (min < 0) {
        for (int i = 0; i < -min; i++) {
            c = c / 2.0;
        }
    }
    if (min > 0) {
        for (int i = 0; i < -min; i++) {
            c = c * 2.0;
        }
    }
    for (int i = min; i <= max; i++) {
        foreach (x0; xs) {
            auto x = c * x0;
            if ( (ctfe_exp(ctfe_log(x)) - x) / x < 1.0e-16 ) {
                // ok
            } else {
                r = r + 1;
            }
        }
    }
    return r;
}

enum c = tests(
    [0.1, 0.1001, 0.11, 0.2, 0.24, 0.3, 0.341, 0.387123, 0.4,
    0.5, 0.55, 0.6, 0.7, 0.732, 0.8, 0.88, 0.9, 0.98, 0.9991, 0.999992],
    -300, 300);
static assert(c == 0);
static assert(ctfe_exp(0.0) == 1.0);
static assert(ctfe_exp(-1000.0) == 0.0);
static assert(ctfe_exp(1000.0) == double.infinity);
static assert(ctfe_log(0.0) == -double.infinity);

-- 
Configure issuemail: http://d.puremagic.com/issues/userprefs.cgi?tab=email
------- You are receiving this mail because: -------