#include "math.h"
#include "fp_private.h"
#if (__WANT_LONG_DOUBLE_FORMAT__ - 0L == 128L)
#include "DD.h"
static const double kLn2 = 6.9314718055994529e-1; static const double kLn2b = 2.3190468138462996e-17; static const double kLn2c = 5.7077084384162121e-34; static const double kLn2d = -3.5824322106018114e-50; static const double k1ByLn2 = 1.4426950408889634; static const double kBig = 6.755399441055744e+15; static const double r13 = 1.6059043836821613e-10;
static const Double kInfinityu = {{0x7ff00000, 0x0}};
static const double coeff[] = {
6227020800.0, 3113510400.0, 1037836800.0, 259459200.0, 51891840.0, 8648640.0,
1235520.0, 154440.0, 17160.0, 1716.0, 156.0, 13.0, 1.0
};
struct ExpTableEntry {
double x;
double fhead;
double ftail;
};
extern uint32_t ExpTableLD[];
long double expl( long double x )
{
double fpenv;
DblDbl u;
double extra;
u.f = x;
if (((u.i[0] & 0x7fffffffu) | u.i[1]) == 0) return 1.0L;
if ((u.i[0] & 0x7ff00000u) != 0x7ff00000u) { FEGETENVD(fpenv);
FESETENVD(0.0);
u.f = _ExpInnerLD(u.d[0], u.d[1], 0.0, &extra, 0);
if ((u.i[0] & 0x7ff00000) == 0x7ff00000)
u.d[1] = 0.0;
FESETENVD(fpenv);
return u.f;
}
if (u.d[0] != u.d[0]) return x;
if ((u.i[0] & 0xfff00000u) == 0x7ff00000u) return x;
else return 0.0L;
}
long double exp2l( long double x )
{
double fpenv;
DblDbl u;
double head, tail, extra, uu, vv, ww, xx, yy, c, top, bottom;
u.f = x;
if (((u.i[0] & 0x7fffffffu) | u.i[1]) == 0) return 1.0L;
if ((u.i[0] & 0x7ff00000u) != 0x7ff00000u) { FEGETENVD(fpenv);
FESETENVD(0.0);
head = u.d[0];
tail = u.d[1];
uu = __FADD( head * kLn2c, tail * kLn2b );
vv = head * kLn2b;
ww = __FMSUB( head, kLn2b, vv );
xx = tail * kLn2;
yy = __FMSUB( tail, kLn2, xx );
c = ww + yy + uu;
top = head * kLn2;
ww = __FMSUB( head, kLn2, top );
uu = __FADD( vv, xx );
if (fabs(xx) > fabs(vv))
c = xx - uu + vv + c;
else
c = vv - uu + xx + c;
bottom = __FADD( uu, ww );
if (fabs(ww) > fabs(uu))
c = ww - bottom + uu + c;
else
c = uu - bottom + ww + c;
head = __FADD( top, bottom );
ww = (top - head) + bottom;
tail = __FADD( ww, c );
c = (ww - tail) + c;
u.f = _ExpInnerLD(head, tail, c, &extra, 3);
if ((u.i[0] & 0x7ff00000) == 0x7ff00000) {
u.d[1] = 0.0;
}
FESETENVD(fpenv);
return u.f;
}
if (u.d[0] != u.d[0]) return x;
if ((u.i[0] & 0xfff00000u) == 0x7ff00000u) return x;
else return 0.0L;
}
long double expm1l( long double x )
{
double fpenv;
DblDbl u;
double extra;
u.f = x;
if (((u.i[0] & 0x7fffffffu) | u.i[1]) == 0) return 0.0L;
if ((u.i[0] & 0x7ff00000u) != 0x7ff00000u) { FEGETENVD(fpenv);
FESETENVD(0.0);
if ((u.i[0] & 0x7ff00000u) >= 0x2ef00000u) u.f = _ExpInnerLD(u.d[0], u.d[1], 0.0, &extra, 2);
if ((u.i[0] & 0x7ff00000) == 0x7ff00000)
u.d[1] = 0.0;
FESETENVD(fpenv);
return u.f;
}
if (u.d[0] != u.d[0]) return x;
if ((u.i[0] & 0xfff00000u) == 0x7ff00000u) return x;
else return -1.0L;
}
long double _ExpInnerLD(double alpha, double beta, double gamma, double *extra,
int exptype)
{
int notnormal, i;
double tail, factor, redarg, b, blow, redarg1, ca, carry;
double arg, arg2a, c, d, arg2, sum, temp, suml, hi, prod, residual;
double sum2, suml2, t1, t2, small, high, prodlow, bump, resnew;
double close, top, residual2, residual3, bottom;
Double test, accndx, en, rscale, scale, power;
DblDbl result;
struct ExpTableEntry *TblPtr = (struct ExpTableEntry *)ExpTableLD + 24;
test.f = alpha;
scale.f = 0.0;
en.f = alpha*k1ByLn2 + kBig;
factor = en.f - kBig;
redarg = __FNMSUB( kLn2, factor, alpha ); notnormal = 0;
if ((test.i[0] & 0x7fffffffu) > 0x40862000u) { if (alpha > 0) {
if (alpha <= 1026.0 * kLn2) {
en.i[1] -= 4;
notnormal = (1023 + 4) * 1048576;
}
else { if ((exptype == 1) || (exptype == 4))
*extra = 0.0;
result.d[0] = kInfinityu.f;
result.d[1] = 0.0;
return result.f;
}
}
else {
if ((alpha >= -800.0) && (exptype != 2)) {
en.i[1] += 256;
notnormal = (1023 - 256)*1048576;
}
else { if ((exptype == 1) || (exptype == 4))
*extra = 0.0;
result.d[0] = (exptype == 2) ? - 1.0 : 0.0;
result.d[1] = 0.0;
return result.f;
}
}
}
b = kLn2b*factor;
blow = __FMSUB( kLn2b, factor, b ); redarg1 = __FNMSUB( kLn2b, factor, beta ); accndx.f = redarg*64.0 + kBig;
ca = kLn2c*factor;
if (fabs(b) > fabs(beta))
carry = beta - (b + redarg1) - blow;
else
carry = (beta - redarg1) - b - blow;
redarg -= TblPtr[(int32_t)accndx.i[1]].x;
arg = __FADD( redarg, redarg1 );
arg2a = (redarg - arg + redarg1);
c = __FSUB( carry, ca );
d = __FNMSUB( kLn2c, factor, ca ) - (ca - (carry - c)) - kLn2d*factor; scale.i[0] = (en.i[1] + 1023) << 20;
arg2 = arg2a + c + d;
if (exptype >= 3)
arg2 += gamma;
sum = coeff[12];
sum = __FMADD( sum, arg, coeff[11] ); sum = __FMADD( sum, arg, coeff[10] ); sum = __FMADD( sum, arg, coeff[9] ); sum = __FMADD( sum, arg, coeff[8] ); temp = __FMADD( sum, arg, coeff[7] ); suml = __FMADD( sum, arg, coeff[7] - temp ); sum = temp;
for (i = 6; i > 0; i--) {
hi = __FMADD( sum, arg, coeff[i] ); suml = __FMADD( suml, arg, __FMADD( sum, arg, coeff[i] - hi ) ); sum = hi;
}
prod = sum*r13;
residual = __FNMSUB( coeff[0], prod, sum ) + suml; suml = residual*r13;
sum2 = __FMADD( prod, arg, 1.0); suml2 = __FMADD( suml, arg, __FMADD( prod, arg, 1.0 - sum2 ) ); sum = __FMADD( suml2, arg, sum2*arg ); suml = __FMADD( suml2, arg, __FMSUB( sum2, arg, sum ) );
t1 = TblPtr[(int32_t)accndx.i[1]].fhead;
t2 = TblPtr[(int32_t)accndx.i[1]].ftail;
small = t2 + sum*(t2 + arg2*t1);
prod = __FMUL( t1, sum );
high = __FADD( prod, t1 );
prodlow = __FMSUB( t1, sum, prod );
residual = (t1 - high) + prod;
if (exptype == 2) {
rscale.f = 0.0;
rscale.i[0] = 2046*1048576 - scale.i[0];
temp = __FSUB( high, rscale.f );
if (temp > 0)
bump = high - temp - rscale.f;
else
bump = high - (temp + rscale.f);
resnew = __FADD( residual, bump );
if (fabs(residual) > fabs(bump))
small = ((residual - resnew) + bump) + small;
else
small = ((bump - resnew) + residual) + small;
high = temp;
residual = resnew;
}
tail = small + arg2*t1 + prodlow + suml*t1;
close = __FADD( tail, residual );
top = __FADD( high, close );
residual2 = (residual - close) + tail;
residual3 = (high - top) + close;
bottom = __FADD( residual3, residual2 );
if ((exptype == 1) || (exptype == 4)) {
if (exptype == 1) scale.f *= 0.5;
if (fabs(residual2) > fabs(residual3))
*extra = (residual2 - bottom + residual3)*scale.f;
else
*extra = (residual3 - bottom + residual2)*scale.f;
}
if (notnormal != 0) {
power.f = 0.0;
power.i[0] = notnormal;
top *= power.f;
bottom *= power.f;
if ((exptype == 1) || (exptype == 4))
*extra *= power.f;
}
result.d[0] = top*scale.f;
result.d[1] = bottom*scale.f;
return result.f;
}
#endif