source: trunk/minix/lib/ack/libp/exp.c@ 9

/*
 * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
 * See the copyright notice in the ACK home directory, in the file "Copyright".
 *
 * Author: Ceriel J.H. Jacobs
 */

/* $Header: /cvsup/minix/src/lib/ack/libp/exp.c,v 1.1 2005/10/10 15:27:46 beng Exp $ */
#define __NO_DEFS
#include <math.h>
#include <pc_err.h>
extern  _trp();

#if __STDC__
#include <float.h>
#include <pc_math.h>
#define M_MIN_D DBL_MIN
#define M_MAX_D DBL_MAX
#define M_DMINEXP DBL_MIN_EXP
#endif
#undef HUGE
#define HUGE    1e1000

static double
Ldexp(fl,exp)
        double fl;
        int exp;
{
        extern double _fef();
        int sign = 1;
        int currexp;

        if (fl<0) {
                fl = -fl;
                sign = -1;
        }
        fl = _fef(fl,&currexp);
        exp += currexp;
        if (exp > 0) {
                while (exp>30) {
                        fl *= (double) (1L << 30);
                        exp -= 30;
                }
                fl *= (double) (1L << exp);
        }
        else    {
                while (exp<-30) {
                        fl /= (double) (1L << 30);
                        exp += 30;
                }
                fl /= (double) (1L << -exp);
        }
        return sign * fl;
}

double
_exp(x)
        double x;
{
        /*      Algorithm and coefficients from:
                        "Software manual for the elementary functions"
                        by W.J. Cody and W. Waite, Prentice-Hall, 1980
        */

        static double p[] = {
                0.25000000000000000000e+0,
                0.75753180159422776666e-2,
                0.31555192765684646356e-4
        };

        static double q[] = {
                0.50000000000000000000e+0,
                0.56817302698551221787e-1,
                0.63121894374398503557e-3,
                0.75104028399870046114e-6
        };
        double  xn, g;
        int     n;
        int     negative = x < 0;

        if (x <= M_LN_MIN_D) {
                g = M_MIN_D/4.0;

                if (g != 0.0) {
                        /* unnormalized numbers apparently exist */
                        if (x < (M_LN2 * (M_DMINEXP - 53))) return 0.0;
                }
                else {
                        if (x < M_LN_MIN_D) return 0.0;
                        return M_MIN_D;
                }
        }
        if (x >= M_LN_MAX_D) {
                if (x > M_LN_MAX_D) {
                        _trp(EEXP);
                        return HUGE;
                }
                return M_MAX_D;
        }
        if (negative) x = -x;

        n = x * M_LOG2E + 0.5;  /* 1/ln(2) = log2(e), 0.5 added for rounding */
        xn = n;
        {
                double  x1 = (long) x;
                double  x2 = x - x1;

                g = ((x1-xn*0.693359375)+x2) - xn*(-2.1219444005469058277e-4);
        }
        if (negative) {
                g = -g;
                n = -n;
        }
        xn = g * g;
        x = g * POLYNOM2(xn, p);
        n += 1;
        return (Ldexp(0.5 + x/(POLYNOM3(xn, q) - x), n));
}