112 lines
1.9 KiB
C
112 lines
1.9 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
|
|
*/
|
|
|
|
/* $Header$ */
|
|
#define __NO_DEFS
|
|
#include <math.h>
|
|
#include <pc_err.h>
|
|
extern _trp();
|
|
|
|
static double
|
|
floor(x)
|
|
double x;
|
|
{
|
|
extern double _fif();
|
|
double val;
|
|
|
|
return _fif(x, 1.0, &val) < 0 ? val - 1.0 : val ;
|
|
/* this also works if _fif always returns a positive
|
|
fractional part
|
|
*/
|
|
}
|
|
|
|
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;
|
|
{
|
|
/* 2**x = (Q(x*x)+x*P(x*x))/(Q(x*x)-x*P(x*x)) for x in [0,0.5] */
|
|
/* Hart & Cheney #1069 */
|
|
|
|
static double p[3] = {
|
|
0.2080384346694663001443843411e+07,
|
|
0.3028697169744036299076048876e+05,
|
|
0.6061485330061080841615584556e+02
|
|
};
|
|
|
|
static double q[4] = {
|
|
0.6002720360238832528230907598e+07,
|
|
0.3277251518082914423057964422e+06,
|
|
0.1749287689093076403844945335e+04,
|
|
0.1000000000000000000000000000e+01
|
|
};
|
|
|
|
int negative = x < 0;
|
|
int ipart, large = 0;
|
|
double xsqr, xPxx, Qxx;
|
|
|
|
if (x < M_LN_MIN_D) {
|
|
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;
|
|
}
|
|
x /= M_LN2;
|
|
ipart = floor(x);
|
|
x -= ipart;
|
|
if (x > 0.5) {
|
|
large = 1;
|
|
x -= 0.5;
|
|
}
|
|
xsqr = x * x;
|
|
xPxx = x * POLYNOM2(xsqr, p);
|
|
Qxx = POLYNOM3(xsqr, q);
|
|
x = (Qxx + xPxx) / (Qxx - xPxx);
|
|
if (large) x *= M_SQRT2;
|
|
x = ldexp(x, ipart);
|
|
if (negative) return 1.0/x;
|
|
return x;
|
|
}
|