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

           Summary: cannot evaluate ylog2x at compile time
           Product: D
           Version: 2.041
          Platform: Other
        OS/Version: Linux
            Status: NEW
          Severity: normal
          Priority: P2
         Component: DMD
        AssignedTo: nobody@puremagic.com
        ReportedBy: baryluk@smp.if.uj.edu.pl


--- Comment #0 from Witold Baryluk <baryluk@smp.if.uj.edu.pl> 2010-01-28 12:19:56 PST ---
import std.stdio;
import std.math;

double iter(double x) {
    static immutable a = log(4.0);
    return x*a;
}

void main() {
    writefln("%s", iter(5.0));
}


/usr/include/d/dmd2-posix/phobos/import/std/math.d(1415): Error: cannot
evaluate yl2x(x,0xb.17217f7d1cf79acp-4L) at compile time
aaaa.d(6): Error: cannot evaluate log(4L) at compile time
aaaa.d(6): Error: cannot evaluate log(4L) at compile time


This also means that currently DMD compiler will not perform constant folding on a.

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

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



--- Comment #1 from Witold Baryluk <baryluk@smp.if.uj.edu.pl> 2010-01-28 12:29:09 PST ---
Same problem is with exp function.


/usr/include/d/dmd2-posix/phobos/import/std/math.d(895): Error: cannot evaluate
exp2(0xb.8aa3b295c17f0bcp-3L * x) at compile time
/usr/include/d/dmd2-posix/phobos/import/std/math.d(903): Error: cannot evaluate
exp(cast(real)x) at compile time
aaaa.d(6): Error: cannot evaluate exp(4F) at compile time
aaaa.d(6): Error: cannot evaluate exp(4F) at compile time

I know it can solved by using CTFE functions, but for example sqrt or sin, cos are working correctly. I really don't want to put by hand obscure numerical constants, and then don't know from where they came :)


Releated to bug1475.

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

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


Witold Baryluk <baryluk@smp.if.uj.edu.pl> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
           Severity|normal                      |minor


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

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


David Simcha <dsimcha@yahoo.com> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |dsimcha@yahoo.com


--- Comment #2 from David Simcha <dsimcha@yahoo.com> 2010-01-28 12:46:12 PST ---
A patch has just been submitted fairly recently for if(__ctfe).  Once that's released, people (myself included) will start submitting patches to make compile-time versions of most of std.math.  Most of the runtime implementations use assembly language somewhere for efficiency, or call the C standard lib implementations.  We'll have to write more naive compile-time versions of these functions.

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

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


Don <clugdbug@yahoo.com.au> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |clugdbug@yahoo.com.au


--- Comment #3 from Don <clugdbug@yahoo.com.au> 2010-01-28 13:01:46 PST ---
(In reply to comment #2)
> A patch has just been submitted fairly recently for if(__ctfe).  Once that's released, people (myself included) will start submitting patches to make compile-time versions of most of std.math.  Most of the runtime implementations use assembly language somewhere for efficiency, or call the C standard lib implementations.  We'll have to write more naive compile-time versions of these functions.

if(__ctfe) is in the next release. yl2x() probably should have built-in
support, though. Since it's an intrinsic, if (__ctfe) won't work for it.

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

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



--- Comment #4 from Witold Baryluk <baryluk@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: -------

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



--- Comment #5 from Witold Baryluk <baryluk@smp.if.uj.edu.pl> 2010-01-31 08:03:28 PST ---
Small but important error in unittest (c was not multiplied correctly. Also -300..300 range was somehow too big, compilations was very long.

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 ) {
                ;
            } else {
                r = r + 1;
            }
        }
        c = c * 2.0;
    }
    return r;
}

unittest {
    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],
        -30, 30);
    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: -------

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


bearophile_hugs@eml.cc changed:

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


--- Comment #6 from bearophile_hugs@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: -------

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



--- Comment #7 from bearophile_hugs@eml.cc 2010-05-12 04:10:38 PDT ---
See bug 4177 for a blocker of this.

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

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



--- Comment #8 from simon <s.d.hammett@googlemail.com> 2010-11-15 05:27:01 PST ---
Created an attachment (id=811)
PATCH against rev 755 add support for yl2x, yl2xp1

Implements yl2x, yl2xp1 as builtins for DMC, VisualStudio. Needs implementation for GCC

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