www.digitalmars.com         C & C++   DMDScript  

digitalmars.D.bugs - [Issue 3749] New: cannot evaluate ylog2x at compile time

reply d-bugmail puremagic.com writes:
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



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: -------
Jan 28 2010
next sibling parent d-bugmail puremagic.com writes:
http://d.puremagic.com/issues/show_bug.cgi?id=3749




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: -------
Jan 28 2010
prev sibling next sibling parent d-bugmail puremagic.com writes:
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: -------
Jan 28 2010
prev sibling next sibling parent d-bugmail puremagic.com writes:
http://d.puremagic.com/issues/show_bug.cgi?id=3749


David Simcha <dsimcha yahoo.com> changed:

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



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: -------
Jan 28 2010
prev sibling next sibling parent d-bugmail puremagic.com writes:
http://d.puremagic.com/issues/show_bug.cgi?id=3749


Don <clugdbug yahoo.com.au> changed:

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




 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: -------
Jan 28 2010
prev sibling next sibling parent d-bugmail puremagic.com writes:
http://d.puremagic.com/issues/show_bug.cgi?id=3749




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: -------
Jan 28 2010
prev sibling next sibling parent d-bugmail puremagic.com writes:
http://d.puremagic.com/issues/show_bug.cgi?id=3749




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: -------
Jan 31 2010
prev sibling next sibling parent d-bugmail puremagic.com writes:
http://d.puremagic.com/issues/show_bug.cgi?id=3749


bearophile_hugs eml.cc changed:

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



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: -------
May 12 2010
prev sibling next sibling parent d-bugmail puremagic.com writes:
http://d.puremagic.com/issues/show_bug.cgi?id=3749




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: -------
May 12 2010
prev sibling next sibling parent d-bugmail puremagic.com writes:
http://d.puremagic.com/issues/show_bug.cgi?id=3749




---
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: -------
Nov 15 2010
prev sibling next sibling parent d-bugmail puremagic.com writes:
http://d.puremagic.com/issues/show_bug.cgi?id=3749


Walter Bright <bugzilla digitalmars.com> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |bugzilla digitalmars.com
           Severity|minor                       |enhancement



22:31:08 PST ---
I'm going to defer this until there's a reasonably portable implementation of
the specifics of yl2x and yl2xp1. Probably one can be built out of the C
standard library math functions.

Failing that, we'll need an inline asm version for gcc, and a 64 bit inline
version.

-- 
Configure issuemail: http://d.puremagic.com/issues/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
Dec 05 2010
prev sibling parent d-bugmail puremagic.com writes:
http://d.puremagic.com/issues/show_bug.cgi?id=3749


Martin Nowak <code dawg.eu> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |code dawg.eu



What's the state of this?
Having log/exp functions at compile time would be very useful, e.g. to
pregenerate scientific tables.

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