124 lines
2 KiB
C
124 lines
2 KiB
C
/*
|
|
* (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
|
|
*/
|
|
|
|
/* $Id$ */
|
|
#define __NO_DEFS
|
|
#include <math.h>
|
|
#include "pc.h"
|
|
|
|
#if __STDC__
|
|
#include <float.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(double fl, int exp)
|
|
{
|
|
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(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));
|
|
}
|