Path: csiph.com!newsfeed.hal-mli.net!feeder3.hal-mli.net!newsfeed.hal-mli.net!feeder1.hal-mli.net!border3.nntp.dca.giganews.com!Xl.tags.giganews.com!border1.nntp.dca.giganews.com!nntp.giganews.com!local2.nntp.dca.giganews.com!nntp.earthlink.com!news.earthlink.com.POSTED!not-for-mail
NNTP-Posting-Date: Mon, 16 Apr 2012 17:12:32 -0500
Message-ID: <4F8C9952.2141@mindspring.com>
Date: Mon, 16 Apr 2012 18:12:34 -0400
From: pete <pfiland@mindspring.com>
Reply-To: pfiland@mindspring.com
Organization: PF
X-Mailer: Mozilla 3.04Gold (WinNT; I)
MIME-Version: 1.0
Newsgroups: sci.math.num-analysis,alt.lang.asm,comp.programming
Subject: Re: Exp() function reloaded
References: <jmdbuh$fql$1@dont-email.me> <4F8A5BF2.17C5@mindspring.com> <4f8bc299$0$1377$4fafbaef@reader2.news.tin.it> <4f8bea40$0$1381$4fafbaef@reader1.news.tin.it>
Content-Type: text/plain; charset=us-ascii
Content-Transfer-Encoding: 7bit
Lines: 130
X-Usenet-Provider: http://www.giganews.com
NNTP-Posting-Host: 4.154.221.186
X-Trace: sv3-FtV4XxG/+0tu78uxnayKsPMKXzs1ay6mKw0IXa3tjO9+ciS7sWTJxtF6/IdbDyAzCsm7fxEZ6IyOgME!NhVz0U8mXSojiV3nonpJZT4TpuBgGTBNSCVe0tO/HmEIojf4CHfSDD+dP9o7VHPb3b5+Pm/aG2rf!MDbmEX8qQQU=
X-Abuse-and-DMCA-Info: Please be sure to forward a copy of ALL headers
X-Abuse-and-DMCA-Info: Otherwise we will be unable to process your complaint properly
X-Postfilter: 1.3.40
X-Original-Bytes: 4296
Xref: csiph.com comp.programming:1479

io_x wrote:
> 
> "io_x" <a@b.c.invalid> ha scritto nel messaggio
> news:4f8bc299$0$1377$4fafbaef@reader2.news.tin.it...
> >
> > "pete" <pfiland@mindspring.com> ha scritto nel messaggio
> > news:4F8A5BF2.17C5@mindspring.com...
> 
> fs_expl(1000000) here would go to seg fault...
> it seems sloppy

It *is* sloppy.
I posted it only to try to make a point about the algorithm
concerning the case when there were negative arguments 
in the function call,
and I'm not sure whether or not the point applies 
to the way that hopcode is implementing the exp function.
I wrote it to be portable 
to all conforming freestanding C implementations,
as an academic exercise in C programming.

The real function is here:

http://www.mindspring.com/~pfilandr/C/fs_math/fs_math.c

#include <float.h>

long double fs_expl(long double x);
long double fs_logl(long double x);

long double fs_expl(long double x)
{
    long unsigned n, square; 
    long double b, e;
    static long double x_max, x_min, epsilon;
    static int initialized;

    if (!initialized) {
        initialized = 1;
        x_max = fs_logl(LDBL_MAX);
        x_min = fs_logl(LDBL_MIN);
        epsilon = LDBL_EPSILON / 4;
    }
    if (x_max >= x && x >= x_min) {
        for (square = 0; x > 1; x /= 2) {
            ++square;
        }
        while (-1 > x) {
            ++square;
            x /= 2;
        }
        e = b = n = 1;
        do {
            b /= n++;
            b *= x;
            e += b;
            b /= n++;
            b *= x;
            e += b;
        } while (b > epsilon);
        while (square-- != 0) {
            e *= e;
        }
    } else {
        e = x > 0 ? LDBL_MAX : 0;
    }
    return e;
}

long double fs_logl(long double x)
{
    long int n;
    long double a, b, c, epsilon;
    static long double A, B, C;
    static int initialized;
    
    if (x > 0 && LDBL_MAX >= x) {
        if (!initialized) {
            initialized = 1;
            B = 1.5;
            do {
                A = B;
                B = 1 / A + A / 2;
            } while (A > B);
            B /= 2;             
            C = fs_logl(A);
        }
        for (n = 0; x > A; x /= 2) {
            ++n;
        }
        while (B > x) {
            --n;
            x *= 2;
        }
        a = (x - 1) / (x + 1);
        x = C * n + a;
        c = a * a;
        n = 1;
        epsilon = LDBL_EPSILON * x;
        if (0 > a) {
            if (epsilon > 0) {
                epsilon = -epsilon;
            }
            do {
                n += 2;
                a *= c;
                b = a / n;
                x += b;
            } while (epsilon > b);
        } else {
            if (0 > epsilon) {
                epsilon = -epsilon;
            }
            do {
                n += 2;
                a *= c;
                b = a / n;
                x += b;
            } while (b > epsilon);
        }
        x *= 2;
    } else {
        x = -LDBL_MAX;
    }
    return x;
}


-- 
pete