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

[d-bugmail at puremagic.com]

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


bearophile_hugs at eml.cc changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |bearophile_hugs at eml.cc


--- Comment #6 from bearophile_hugs at eml.cc 2010-05-12 03:12:43 PDT ---
For a generic compile-time function fit to be used with D FP types, it can be
better to use reals everywhere:


/** 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)
 *
 */
real ctfe_log(real x_) {
    if (x_ == 0.0) {
        return -real.infinity;
    }
    if (x_ < 0.0) {
        return real.nan;
    }
    if (x_ == real.infinity) {
        return real.infinity;
    }

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

    if (!(x_ == x_)) { // nan
        return real.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)
 */
real ctfe_exp(real x_) {
    if (x_ == 0.0) {
        return 1.0;
    }

    if (x_ >= 710.0) { // includes +infinity
        return real.infinity;
    }
    if (x_ <= -746.0) { // includes -infinity
        return 0.0;
    }
    if (!(x_ == x_)) { // nan
        return real.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;
    }
}

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