#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 kBig = 6.755399441055744e+15; static const double recip15fac = 0.138750138750138750138750e-5; static const double scaleup = 1.3407807929942597e+154; static const double kLog10eHi = 4.3429448190325182e-1; static const double kLog10eMid = 1.0983196502167651e-17; static const double kLog10eLow = 3.7171812331109602e-34; static const double kLog2eHi = 1.4426950408889634e+0; static const double kLog2eMid = 2.0355273740931033e-17; static const double kLog2eLow = -1.0614659956117258e-33;
static const Double kInfinityu = {{0x7ff00000, 0x0}};
static const double coeff[] = {720720.0, -360360.0, 240240.0, -180180.0,
144144.0, -120120.0, 102960.0, -90090.0, 80080.0, -72072.0, 65520.0,
-60060.0, 55440.0, -51480.0, 48048.0, -45045.0};
extern uint32_t LogTableLD[];
struct LogTableEntry {
double x;
double u; double fhead;
double ftail;
};
long double logl( long double x )
{
double fpenv;
DblDbl u;
double extra;
FEGETENVD(fpenv);
FESETENVD(0.0);
u.f = x;
if (((u.i[0] & 0x7fffffffu) | u.i[1]) == 0) { u.d[0] = -1.0/__FABS(u.d[0]);
u.d[1] = 0.0;
FESETENVD(fpenv);
return u.f;
}
if (u.i[0] < 0x7ff00000u) { if ((u.i[0] == 0x3ff00000) && (u.i[1] | u.i[2] | u.i[3]) == 0) {
FESETENVD(fpenv);
return 0.0L; }
u.f = _LogInnerLD(u.d[0], u.d[1], 0.0, &extra, 0);
FESETENVD(fpenv);
return u.f;
}
if (u.d[0] != u.d[0]) {
FESETENVD(fpenv);
return x; }
if ((u.i[0] & 0xfff00000u) == 0x7ff00000u) {
FESETENVD(fpenv);
return x; }
else { u.d[0] = 0.0 * kInfinityu.f + u.d[0];
u.d[1] = 0.0;
FESETENVD(fpenv);
return u.f;
}
}
long double log10l( long double x )
{
double fpenv;
DblDbl u;
double result, resmid, resbot;
double uu, vv, ww, xx, yy, c, top, bottom;
FEGETENVD(fpenv);
FESETENVD(0.0);
u.f = x;
if (((u.i[0] & 0x7fffffffu) | u.i[1]) == 0) { u.d[0] = -1.0/__FABS(u.d[0]);
u.d[1] = 0.0;
FESETENVD(fpenv);
return u.f;
}
if (u.i[0] < 0x7ff00000u) { if ((u.i[0] == 0x3ff00000) && (u.i[1] | u.i[2] | u.i[3]) == 0) {
FESETENVD(fpenv);
return 0.0L; }
u.f = _LogInnerLD(u.d[0], u.d[1], 0.0, &resbot, 1);
result = u.d[0];
resmid = u.d[1];
uu = __FADD( result * kLog10eLow, __FADD( resmid * kLog10eMid, resbot * kLog10eHi ) );
vv = result * kLog10eMid;
ww = __FMSUB( result, kLog10eMid, vv);
xx = resmid * kLog10eHi;
yy = __FMSUB( resmid, kLog10eHi, xx );
c = ww + yy + uu;
top = result * kLog10eHi;
ww = __FMSUB( result, kLog10eHi, 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;
u.d[0] = __FADD( top, bottom );
u.d[1] = top - u.d[0] + bottom + c;
FESETENVD(fpenv);
return u.f;
}
if (u.d[0] != u.d[0]) {
FESETENVD(fpenv);
return x; }
if ((u.i[0] & 0xfff00000u) == 0x7ff00000u) {
FESETENVD(fpenv);
return x; }
else { u.d[0] = 0.0 * kInfinityu.f + u.d[0];
u.d[1] = 0.0;
FESETENVD(fpenv);
return u.f;
}
}
long double log2l( long double x )
{
double fpenv;
DblDbl u;
double result, resmid, resbot;
double uu, vv, ww, xx, yy, c, top, bottom;
FEGETENVD(fpenv);
FESETENVD(0.0);
u.f = x;
if (((u.i[0] & 0x7fffffffu) | u.i[1]) == 0) { u.d[0] = -1.0/__FABS(u.d[0]);
u.d[1] = 0.0;
FESETENVD(fpenv);
return u.f;
}
if (u.i[0] < 0x7ff00000u) { if ((u.i[0] == 0x3ff00000) && (u.i[1] | u.i[2] | u.i[3]) == 0) {
FESETENVD(fpenv);
return 0.0L; }
u.f = _LogInnerLD(u.d[0], u.d[1], 0.0, &resbot, 1);
result = u.d[0];
resmid = u.d[1];
uu = __FADD( result * kLog2eLow, __FADD( resmid * kLog2eMid, resbot * kLog2eHi ) );
vv = result * kLog2eMid;
ww = __FMSUB( result, kLog2eMid, vv );
xx = resmid * kLog2eHi;
yy = __FMSUB( resmid, kLog2eHi, xx );
c = ww + yy + uu;
top = result * kLog2eHi;
ww = __FMSUB( result, kLog2eHi, 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;
u.d[0] = __FADD( top, bottom );
u.d[1] = top - u.d[0] + bottom + c;
FESETENVD(fpenv);
return u.f;
}
if (u.d[0] != u.d[0]) {
FESETENVD(fpenv);
return x; }
if ((u.i[0] & 0xfff00000u) == 0x7ff00000u) {
FESETENVD(fpenv);
return x; }
else { u.d[0] = 0.0 * kInfinityu.f + u.d[0];
u.d[1] = 0.0;
FESETENVD(fpenv);
return u.f;
}
}
long double log1pl( long double x )
{
double fpenv;
DblDbl u;
double high, mid, low, extra, a, b, c, d, e, f;
u.f = x;
if (u.d[0] != u.d[0]) return x;
if (((u.i[0] & 0x7fffffffu) | u.i[1]) == 0) return x;
FEGETENVD(fpenv);
FESETENVD(0.0);
if ((x > -1.0L) && ((u.i[0] & 0x7ff00000) != 0x7ff00000)) { if (1.0 > __FABS(u.d[0])) { high = 1.0;
mid = u.d[0];
low = u.d[1];
}
else if (__FABS(u.d[1]) > 1.0) { high = u.d[0];
mid = u.d[1];
low = 1.0;
}
else { high = u.d[0];
mid = 1.0;
low = u.d[1];
}
a = __FADD( mid, low );
b = (mid - a) + low;
c = __FADD( high, a );
d = (high - c) + a;
e = __FADD( d, b );
f = (d - e) + b;
u.f = _LogInnerLD(c, e, f, &extra, 2);
FESETENVD(fpenv);
return u.f;
}
if ((u.i[0] == 0xbff00000) && (u.i[1] | (u.i[2] & 0x7fffffffu) | u.i[3]) == 0) {
u.d[0] = -INFINITY;
u.d[1] = 0.0;
FESETENVD(fpenv);
return u.f;
}
if ((u.i[0] & 0xfff00000u) == 0x7ff00000u) {
FESETENVD(fpenv);
return x; }
else { u.d[0] = 0.0 * kInfinityu.f + u.d[0];
u.d[1] = 0.0;
FESETENVD(fpenv);
return u.f;
}
}
long double _LogInnerLD(double alpha, double beta, double gamma, double *extra,
int logtype)
{
int i, iexp, ndx;
double head, tail, redarg, diff, diff2, quot, arg, rem;
double reciparg, e, sum, arglow, u, x, y, remlost, arglost, temp, suml;
double argl, argl2, hi, prod, c, enlog2hi, w, enlog2mid, prodl, sum2;
double sum3, enlog2low, small, petit, smallres, summid, resmid, summid1;
double resmid1, summid2, resmid2, top, resmid3, resmid4, bottom, resmid5;
double residual, result, residual2, resmidtemp;
Double test, rscale, scale;
DblDbl out;
struct LogTableEntry *TblPtr = (struct LogTableEntry *)LogTableLD;
test.f = alpha;
tail = beta;
rscale.f = 0.0;
scale.f = kBig;
if (test.i[0] >= 0x07000000u) { iexp = (test.i[0] & 0x7ff00000);
rscale.i[0] = 2046 * 1048576 - iexp;
scale.i[1] = (iexp >> 20);
}
else { test.f *= scaleup;
tail *= scaleup;
if (logtype == 2) gamma *= scaleup;
head = __FADD( test.f, tail );
tail = test.f - head + tail;
iexp = (test.i[0] & 0x7ff00000);
rscale.i[0] = 2046 * 1048576 - iexp;
scale.i[1] = (iexp >> 20);
scale.f -= 512.0;
}
ndx = ((test.i[0] & 0x000fc000) >> 14);
if (rscale.i[0] == 0)
rscale.i[0] = 524288;
head = test.f * rscale.f;
tail = tail * rscale.f;
if (logtype == 2)
gamma = gamma * rscale.f;
scale.f -= kBig + 1023.0;
if (test.i[0] & 0x00080000) {
scale.f += 1.0;
head *= 0.5;
tail *= 0.5;
if (logtype == 2)
gamma *= 0.5;
}
if (ndx == 0)
if ((test.i[0] & 0x000fffffu) < 0x00002000u)
ndx = 63;
redarg = head - TblPtr[ndx].x;
diff = __FADD( redarg, tail );
diff2 = redarg - diff + tail + gamma;
quot = __FMUL( diff, TblPtr[ndx].u );
arg = __FMADD( __FNMSUB( quot, TblPtr[ndx].x, diff ), TblPtr[ndx].u, quot ); remlost = __FNMSUB( arg, TblPtr[ndx].x, diff );
rem = __FADD( remlost, diff2 );
reciparg = 1.0 - arg;
e = arg * arg;
sum = __FMADD( arg, coeff[15], coeff[14] ); arglow = rem * TblPtr[ndx].u;
u = scale.f * kLn2c;
x = scale.f * kLn2b;
y = __FMSUB( scale.f, kLn2b, x );
if (__FABS(remlost) > __FABS(diff2))
remlost = remlost - rem + diff2;
else
remlost = diff2 - rem + remlost;
arglost = __FNMSUB( arglow, TblPtr[ndx].x, rem );
arglost = __FADD( remlost, arglost );
arglost = arglost * TblPtr[ndx].u;
sum = __FMADD( arg, sum, coeff[13] ); reciparg = __FMADD( reciparg, e, reciparg ); sum = __FMADD( arg, sum, coeff[12] ); e = e * e;
sum = __FMADD( arg, sum, coeff[11] ); reciparg = __FMADD( reciparg, e, reciparg ); sum = __FMADD( arg, sum, coeff[10] ); e = e * e;
sum = __FMADD( arg, sum, coeff[ 9] ); reciparg = __FMADD( reciparg, e, reciparg ); temp = __FMADD( arg, sum, coeff[8] ); suml = __FSUB( coeff[8], temp );
suml = __FMADD( arg, sum, suml ); sum = temp;
argl = arglow * reciparg;
argl2 = (__FNMSUB( argl, arg, arglow - argl) + arglost) *reciparg; for (i = 7; i >= 2; i--) {
hi = __FMADD( sum, arg, coeff[i] );
suml = __FMADD( suml, arg, __FMADD( sum, arg, coeff[i] - hi ) ); sum = hi;
}
prod = sum * recip15fac;
c = u + y;
enlog2hi = __FMUL( scale.f, kLn2 );
w = __FMSUB( scale.f , kLn2 , enlog2hi );
enlog2mid = __FADD( x , w );
prodl = (__FNMSUB( prod, coeff[0], sum ) + suml) * recip15fac; sum2 = __FMSUB( prod, arg, 0.5);
suml = __FMADD( prodl, arg, __FMADD( prod, arg, __FMSUB( -1.0, 0.5, sum2 ) ) ); sum3 = __FMADD( sum2, arg, (suml * arg) );
suml = __FMADD( suml, arg, __FMSUB( sum2, arg, sum3 ) ); sum = __FMADD( sum3, arg, (suml * arg) );
suml = __FMADD( suml, arg, __FMSUB( sum3, arg, sum ) ); if (scale.f == 0.0) {
top = __FADD( arg, sum );
resmid = arg - top + sum;
small = __FADD( suml, argl );
smallres = argl - small + suml;
if (TblPtr[ndx].x == 1.0) {
bottom = __FADD( resmid, small );
residual2 = smallres + argl2;
residual = resmid - bottom + small + residual2;
}
else {
hi = __FADD( TblPtr[ndx].fhead, top );
resmid2 = TblPtr[ndx].fhead - hi + top + resmid;
petit = small;
bottom = __FADD( resmid2, petit );
if (__FABS(resmid2) > __FABS(petit))
residual = resmid2-bottom+petit+smallres+argl2+TblPtr[ndx].ftail;
else
residual = petit-bottom+resmid2+smallres+argl2+TblPtr[ndx].ftail;
top = hi;
}
}
else {
if (__FABS(x) > __FABS(w))
enlog2low = x - enlog2mid + w + c;
else
enlog2low = w - enlog2mid + x + c;
small = __FADD( argl, suml );
petit = TblPtr[ndx].ftail + enlog2low + argl2 + small;
smallres = argl - small + suml;
summid = __FADD( sum, enlog2mid );
if (__FABS(sum) > __FABS(enlog2mid))
resmid = sum - summid + enlog2mid;
else
resmid = enlog2mid - summid + sum;
summid1 = __FADD( arg, summid );
resmid1 = arg - summid1 + summid;
summid2 = __FADD( summid1, TblPtr[ndx].fhead );
resmid2 = TblPtr[ndx].fhead - summid2 + summid1 + resmid1;
petit = petit + resmid + smallres;
top = __FADD( enlog2hi, summid2 );
resmid3 = enlog2hi - top + summid2;
resmid4 = __FADD( resmid3, resmid2 );
bottom = __FADD( petit, resmid4 );
resmid5 = resmid3 - resmid4 + resmid2;
residual = resmid4 - bottom + petit + resmid5;
}
result = __FADD( top, bottom );
if (logtype != 1)
resmid = top - result + bottom + residual;
else { resmidtemp = top - result + bottom;
resmid = __FADD( resmidtemp, residual );
*extra = resmidtemp - resmid + residual;
}
out.d[0] = result;
out.d[1] = resmid;
return out.f;
}
#endif