#include "math.h"
#include "fp_private.h"
#if (__WANT_LONG_DOUBLE_FORMAT__ - 0L == 128L)
#include "DD.h"
void _d3div(double a, double b, double c, double d, double e, double f,
double *high, double *mid, double *low);
void _d3mul(double a, double b, double c, double d, double e, double f,
double *high, double *mid, double *low);
void _d3add(double a, double b, double c, double d, double e, double f,
double *high, double *mid, double *low);
long double _LogInnerLD(double, double, double, double *, int);
long double _ExpInnerLD(double alpha, double beta, double gamma, double *extra,
int exptype);
long double rintl(long double x);
long double sinl(long double x);
typedef double DD[2];
static const uint32_t bernData[] __attribute__ ((aligned(8))) = {
0x3FB55555, 0x55555555, 0x3C555555, 0x55555555,
0xBF66C16C, 0x16C16C17, 0x3BFF49F4, 0x9F49F49F,
0x3F4A01A0, 0x1A01A01A, 0x3B8A01A0, 0x1A01A01A,
0xBF438138, 0x13813814, 0x3BEFB1FB, 0x1FB1FB20,
0x3F4B951E, 0x2B18FF23, 0x3BE5C3A9, 0xCE01B952,
0xBF5F6AB0, 0xD9993C7D, 0x3BFF8255, 0x3C999B0E,
0x3F7A41A4, 0x1A41A41A, 0x3C106906, 0x90690690,
0xBF9E4286, 0xCB0F5398, 0x3C21EFCD, 0xAB896745,
0x3FC6FE96, 0x381E0680, 0xBC279E24, 0x05A71F88,
0xBFF64767, 0x01181F3A, 0x3C724246, 0x319DA678,
0x402ACE44, 0x322CE006, 0xBCC62C2B, 0x1BBCDD32
};
static const uint32_t zetaData[] __attribute__ ((aligned(8))) = {
0x3FD4A34C, 0xC4A60FA6, 0x3C71873D, 0x8912200C,
0xBFB13E00, 0x1A557607, 0x3C5FB68B, 0xE2F8821F,
0x3F951322, 0xAC7D8483, 0x3C3AFC89, 0x088CB729,
0xBF7E404F, 0xC218F5F2, 0x3C1E4A62, 0x7CF1EB34,
0x3F67ADD6, 0xEADB6C30, 0xBBF5B782, 0x8C7FD7F4,
0xBF538AC5, 0xC2BF8E08, 0x3BE8A4C1, 0xCFD9CEC8,
0x3F40B36A, 0xF86396E9, 0xBBE0698D, 0x6C892967,
0xBF2D3FD4, 0xC76D2FC8, 0x3BBC7C55, 0xCFCCBB83,
0x3F1A127B, 0x0F17D65A, 0x3BA9D309, 0xAA700268,
0xBF078DE5, 0xBD7C81EF, 0x3B7A2054, 0x1CDE47A6,
0x3EF580DC, 0xEE66EB02, 0x3B826057, 0x4B258F72,
0xBEE3CBC9, 0x63CE2243, 0x3B8EA56E, 0x6C7D5329,
0x3ED2597A, 0x39F34AAC, 0xBB7BF911, 0x462A7D81,
0xBEC11B2E, 0xB7679541, 0xBB4C76B0, 0xE65AC63A,
0x3EB0064C, 0xDEB22F0F, 0x3B4D0156, 0xAFFDBC11,
0xBE9E2600, 0xD93CFD2F, 0x3B3130AC, 0x39E5C106,
0x3E8C76BB, 0xB3F07A4D, 0x3B2D9A2B, 0x77769B52,
0xBE7AF5A6, 0xCBBF8A97, 0xBB195F22, 0x7E96D83E,
0x3E699B93, 0xC2070B0F, 0x3B003271, 0x64736428,
0xBE5862C7, 0x34DF3EAC, 0xBAFB3280, 0x2BEC0DA0,
0x3E47469D, 0xACCFADCD, 0xBAE369D3, 0x88CEBAA9,
0xBE36434A, 0x8447AEAD, 0xBA8AF72E, 0xDF876FCD,
0x3E2555A8, 0x77FFD2C3, 0xBAC87506, 0x5F26A43B,
0xBE147B16, 0x79258D0E, 0xBAB04F36, 0xE0E854E4,
0x3E03B15D, 0x2B2FC10C, 0xBA9D79F6, 0xFEEEB28B,
0xBDF2F69A, 0x9FABE3E0, 0x3A9A162A, 0xB374C789,
0x3DE24932, 0xA337434C, 0x3A806082, 0x9C24508F,
0xBDD1A7C2, 0x6EC2523C, 0xBA74F4EB, 0xDB4A04B5,
0x3DC11116, 0xE693ED98, 0xBA6C7034, 0xD49E7FC7,
0xBDB08424, 0xCBC543D8, 0xBA440EF8, 0x20DBC9EA,
0x3DA00002, 0x6E3F644F, 0x3A43546A, 0x6054C889,
0xBD8F07C5, 0x14FC9F0A, 0xB9F75B6B, 0xE545AC09,
0x3D7E1E20, 0x26AAFCD8, 0xBA162A85, 0x86538620,
0xBD6D41D5, 0x6E5EE2E2, 0x39F43894, 0xD27CED5E,
0x3D5C71C7, 0xF6F10E37, 0xB9F01074, 0x764D33F2,
0xBD4BACF9, 0xA27BC89B, 0x39D4A5A2, 0x15E0508E,
0x3D3AF287, 0x18A10D6E, 0x39C40D7F, 0x1B842CC3,
0xBD2A41A4, 0x5603E5B6, 0x39C62BE9, 0xCF212D90,
0x3D199999, 0xC0716EE9, 0xB9B39E10, 0xF90435CB,
0xBD08F9C1, 0xA8DF9D78, 0x3989DA56, 0xD4471922,
0x3CF86186, 0x28D28905, 0xB989D7D4, 0xEE5A8893,
0xBCE7D05F, 0x4C31C560, 0xB9871BBA, 0x0B7CC339,
0x3CD745D1, 0x7B56BA4A, 0x3979D38B, 0xC00D6F02,
0xBCC6C16C, 0x1B4D6456, 0xB96AED17, 0x2E5C90F2,
0x3CB642C8, 0x5C023D9D, 0xB95DE052, 0x190D755D,
0xBCA5C988, 0x2D825E9D, 0x3949723F, 0x1BF240B7,
0x3C955555, 0x5698A866, 0x392CF5C8, 0x64970970,
0xBC84E5E0, 0xA8022BC9, 0x39128B9D, 0xC88F5BE2,
0x3C747AE1, 0x4838081F, 0xB91DF461, 0x306465C0,
0xBC641414, 0x146E3E31, 0xB8FE4773, 0xEA1309D6,
0x3C53B13B, 0x13EC2F3A, 0x38C41C5B, 0x07A32CA9,
0xBC43521C, 0xFB520859, 0xB8E225B1, 0x0AA3CE51
};
static const double zero = 0.0;
static const Double pi = {{0x400921FB, 0x54442D18}};
static const Double pib = {{0x3CA1A626, 0x33145C07}};
static const Double pic = {{0xB92F1976, 0xB7ED8FBC}};
static const Double infinity = {{0x7ff00000, 0x00000000}};
static const Double scaleup = {{0x4ff00000, 0x00000000}};
static const Double tiny = {{0x0E000000, 0x00000000}};
static const Double sclog = {{0x40662E42, 0xFEFA39EF}};
static const Double sclog2 = {{0x3CFABC9E, 0x3B39803F}};
static const Double sclog3 = {{0x3987B57A, 0x079A1934}};
static const Double bump = {{0x3FED67F1, 0xC864BEB5}}; static const Double bump2 = {{0xBC865B5A, 0x1B7FF5DF}};
static const Double bump3 = {{0xB91B7F70, 0xC13DC168}};
static const Double ec = {{0x3FDB0EE6, 0x072093CE}}; static const Double ec2 = {{0x3C56CB90, 0x701FBFAB}};
static const Double ec3 = {{0x38F34A95, 0xE3133C51}};
static const Double piln = {{0x3FF250D0, 0x48E7A1BD}}; static const Double pilnb = {{0x3C67ABF2, 0xAD8D5088}};
static const Double pilnc = {{0xB8E6CCF4, 0x32447B52}};
static const Double zeta32 = {{0xB914C685, 0x28DDC956}};
static long double GammaCommon( double dhead, double dtail, int igamma )
{
int signgam, ireflect, increase, i;
double int1, int2, arg, arg2, arg3, argt, arg2t, anew;
double factor, factor1, factor2, factor3, f1, f2, f3;
double prod1, prod2, pextra, pextra1, pextra2, sum1, sum2, sextra;
double sum21, sum22, sextra2, sum31, sum32, sextra3;
double diff1, diff2, dextra, sarg1, sarg2, saextra;
double denom1, denom2, sum, suml, sumt, recip1, recip2, recipsq1, recipsq2;
double extra, extra2, extra3, extrax, extray, extralg, gmextra;
double temp, eexp, diff, z, z2, y, y2, y3, asize, c;
DblDbl gamdd, tempdd, intdd, xhintdd, sinargdd;
DD *bern = (DD *)(bernData - 4);
DD *zeta = (DD *)(zetaData - 8);
signgam = 1;
if (dhead <= 0) { tempdd.d[0] = dhead;
tempdd.d[1] = dtail;
intdd.f = rintl(tempdd.f); int1 = intdd.d[0];
int2 = intdd.d[1];
if (isinf(dhead) || (dhead - int1 + dtail - int2) == 0) {
gamdd.d[0] = INFINITY; gamdd.d[1] = zero;
return gamdd.f;
}
if (dhead > -20.0) {
ireflect = 0; }
else {
ireflect = 1; dhead = -dhead; dtail = -dtail;
}
}
else ireflect = 0;
if (isinf(dhead) || (dhead + dtail) != (dhead + dtail)) {
gamdd.d[0] = dhead + dtail;
gamdd.d[1] = zero;
return gamdd.f;
}
arg = dhead; arg2 = dtail;
arg3 = 0.0;
factor1 = 1.0;
factor2 = 0.0;
factor3 = 0.0;
if (arg > 18.0) { while (arg < 31.0) {
_d3mul( factor1, factor2, factor3, arg, arg2, arg3,
&factor1, &factor2, &factor3 );
arg = arg + 1.0;
} while (arg < 35.0) {
_d3mul( factor1, factor2, factor3, arg, arg2, arg3,
&factor1, &factor2, &factor3 );
argt = __FADD( arg, 1.0 );
arg2t = __FADD( (arg - argt + 1.0), arg2 );
arg3 = ((arg - argt + 1.0) - arg2t + arg2) + arg3;
arg = argt;
arg2 = arg2t;
}
tempdd.f = _LogInnerLD(arg, arg2, 0.0, &f3, 1); f3 = f3 + arg3/arg;
f1 = tempdd.d[0];
f2 = tempdd.d[1];
_d3mul(f1, f2, f3, (arg - 0.5), arg2, arg3,
&prod1, &prod2, &c); _d3add(prod1, prod2, c, -arg, -arg2, -arg3,
&sum1, &sum2, &sextra); _d3add(sum1, sum2, sextra, bump.f, bump2.f, bump3.f,
&sum21, &sum22, &sextra2);
_d3div(1.0, 0.0, 0.0, arg, arg2, arg3, &recip1, &recip2, &extra);
_d3mul(recip1, recip2, extra, recip1, recip2, extra,
&recipsq1, &recipsq2, &extra2);
sum = bern[11][0];
suml = bern[11][1];
for (i = 10; i >= 1; --i) {
temp = __FMADD( sum, recipsq1, bern[i][0] );
suml = __FMADD( sum, recipsq1, bern[i][0] - temp ) +
__FMADD( suml, recipsq1, __FMADD( sum, recipsq2, bern[i][1] ) );
sum = temp;
}
_d3mul(sum, suml, 0.0, recip1, recip2, extra, &prod1, &prod2, &extra3);
_d3add(sum21, sum22, sextra2, prod1, prod2, extra3, &sum31, &sum32, &sextra3);
if ((igamma == 1) && (ireflect != 1)) {
tempdd.f = _ExpInnerLD(sum31, sum32, sextra3, &eexp, 4);
if (factor1 == 1.0) {
gamdd.f = tempdd.f;
gmextra = eexp;
}
else {
_d3div(tempdd.d[0], tempdd.d[1], eexp, factor1, factor2, factor3,
&(gamdd.d[0]), &(gamdd.d[1]), &gmextra);
}
}
else { factor = factor1;
if (factor == 1.0) {
gamdd.d[0] = sum31;
gamdd.d[1] = sum32;
gmextra = sextra3;
}
else {
tempdd.f = _LogInnerLD(factor1, factor2, 0.0, &extrax, 1);
extray = -(extrax + factor3/factor);
_d3add(sum31, sum32, sextra3, -tempdd.d[0], -tempdd.d[1], extray,
&(gamdd.d[0]), &(gamdd.d[1]), &gmextra);
}
}
}
else { arg = dhead;
arg2 = dtail; arg3 = 0.0;
increase = 0; if (arg < 1.5) { increase = -1; factor1 = arg;
factor2 = arg2;
factor3 = 0.0; while (arg < 1.5) {
_d3add(arg, arg2, arg3, 1.0, 0.0, 0.0, &arg, &arg2, &arg3);
if (arg < 1.5) {
_d3mul(factor1, factor2, factor3, arg, arg2, arg3,
&factor1, &factor2, &factor3);
}
}
if (factor1 < 0.0) {
signgam = -1; }
}
else if (arg > 2.5) {
increase = +1; factor1 = 1.0;
factor2 = 0.0;
factor3 = 0.0;
while (arg > 2.5) { arg = arg - 1.0; anew = __FADD( arg, arg2 ); arg2 = arg - anew + arg2; arg = anew;
_d3mul(factor1, factor2, factor3, arg, arg2, 0.0,
&factor1, &factor2, &factor3);
}
arg3 = 0.0;
}
diff = __FSUB( arg, 2.0 ); z = __FADD( (arg - (diff + 2.0)), arg2 ); z2 = (arg - (diff + 2.0)) - z + arg2 + arg3;
y = __FADD( diff, z );
y2 = (diff - y + z) + z2;
y3 = (diff - y + z) - y2 + z2;
sum = zeta[53][0];
suml = zeta[53][1];
for (i = 52; i >= 40; --i) {
sum = __FMADD( sum, y, zeta[i][0] ); }
sumt = __FMADD( sum, y, zeta[39][0] ); suml = __FMADD( sum, y, zeta[39][0] - sumt) + __FMADD( sum, y2, zeta[39][1] ); sum = sumt;
for (i = 38; i >= 3; --i) {
temp = __FMADD( sum, y, zeta[i][0] ); suml = __FMADD( sum, y, zeta[i][0] - temp ) + zeta[i][1] + __FMADD( sum, y2, suml*y );
sum = temp;
}
_d3mul(sum, suml, 0.0, y, y2, y3, &prod1, &prod2, &pextra2);
_d3add(prod1, prod2, pextra2, zeta[2][0], zeta[2][1], zeta32.f,
&sum1, &sum2, &sextra);
_d3mul(sum1, sum2, sextra, y, y2, y3,
&prod1, &prod2, &pextra1); _d3add(prod1, prod2, pextra1, ec.f, ec2.f, ec3.f,
&sum1, &sum2, &sextra); _d3mul(sum1, sum2, sextra, y, y2, y3,
&prod1, &prod2, &pextra);
if (igamma == 1) { tempdd.f = _ExpInnerLD(prod1, prod2, pextra, &eexp, 4);
if (increase == 1) {
_d3mul(tempdd.d[0], tempdd.d[1], eexp, factor1, factor2, factor3,
&(gamdd.d[0]), &(gamdd.d[1]), &gmextra);
}
else if (increase == -1) {
_d3div(tempdd.d[0], tempdd.d[1], eexp, factor1, factor2, factor3,
&(gamdd.d[0]), &(gamdd.d[1]), &gmextra);
}
else {
gamdd.f = tempdd.f;
gmextra = eexp;
}
}
else { if (increase == 0) {
gamdd.d[0] = prod1;
gamdd.d[1] = prod2;
gmextra = pextra;
}
else {
if (signgam < 0) {
factor1 = -factor1;
factor2 = -factor2;
factor3 = -factor3;
}
factor = factor1;
tempdd.f = _LogInnerLD(factor1, factor2, 0.0, &extra, 1);
extra = extra + factor3/factor;
if (increase == -1) { tempdd.d[0] = -tempdd.d[0];
tempdd.d[1] = -tempdd.d[1];
extra = -extra;
}
_d3add(prod1, prod2, pextra, tempdd.d[0], tempdd.d[1], extra,
&(gamdd.d[0]), &(gamdd.d[1]), &gmextra);
}
}
}
if (gamdd.d[0] != gamdd.d[0]) {
gamdd.d[0] = infinity.f;
gamdd.d[1] = zero;
gmextra = 0.0;
}
if (ireflect == 1) { arg = -dhead; arg2 = -dtail;
_d3add(arg, arg2, 0.0, -int1, -int2, 0.0, &diff1, &diff2, &dextra);
asize = fabs(diff1);
if (asize < tiny.f) { diff1 *= scaleup.f; diff2 *= scaleup.f; dextra *= scaleup.f; }
_d3mul(diff1, diff2, dextra, pi.f, pib.f, pic.f, &sarg1, &sarg2, &saextra);
xhintdd.f = 2.0*rintl(0.5*intdd.f);
tempdd.d[0] = sarg1;
tempdd.d[1] = sarg2;
sinargdd.f = sinl(tempdd.f);
if ((xhintdd.d[0] - intdd.d[0] + xhintdd.d[1] - intdd.d[1]) != 0.0) {
sinargdd.f = -sinargdd.f;
}
if (sinargdd.d[0] < 0.0) {
sinargdd.f = -sinargdd.f; signgam = -signgam; }
if (fabs(gamdd.d[0]) == infinity.f){ if(igamma == 1) {
gamdd.d[0] = zero;
gamdd.d[1] = zero;
}
else {
gamdd.d[0] = -infinity.f; gamdd.d[1] = zero; }
return gamdd.f;
}
_d3mul(dhead, dtail, 0.0, sinargdd.d[0], sinargdd.d[1], 0.0,
&prod1, &prod2, &pextra);
tempdd.f = _LogInnerLD(prod1, prod2, 0.0, &extralg, 1); extralg = extralg + pextra/prod1;
_d3add(tempdd.d[0], tempdd.d[1], extralg, gamdd.d[0], gamdd.d[1], gmextra,
&denom1, &denom2, &dextra);
_d3add(piln.f, pilnb.f, pilnc.f, -denom1, -denom2, -dextra,
&(gamdd.d[0]), &(gamdd.d[1]), &gmextra);
if (asize < tiny.f) {
_d3add(gamdd.d[0], gamdd.d[1], gmextra, sclog.f, sclog2.f, sclog3.f,
&(gamdd.d[0]), &(gamdd.d[1]), &gmextra);
}
if (igamma == 1) { tempdd.f = _ExpInnerLD(gamdd.d[0], gamdd.d[1], gmextra, &eexp, 4);
if (signgam == 1) {
gamdd.f = tempdd.f;
}
else {
gamdd.f = -tempdd.f;
}
}
}
return gamdd.f;
}
long double tgammal( long double x )
{
double head, fenv;
DblDbl t;
if (x == 0.0L)
return copysignl( INFINITY, x );
FEGETENVD(fenv);
FESETENVD(0.0);
t.f = x;
t.f = GammaCommon(t.d[0], t.d[1], 1);
head = __FADD( t.d[0], t.d[1] );
t.d[1] = (t.d[0] - head) + t.d[1]; t.d[0] = head;
FESETENVD(fenv);
return t.f;
}
long double lgammal( long double x )
{
double head, fenv;
DblDbl t;
FEGETENVD(fenv);
FESETENVD(0.0);
t.f = x;
t.f = GammaCommon(t.d[0], t.d[1], 0);
head = __FADD( t.d[0], t.d[1] );
t.d[1] = (t.d[0] - head) + t.d[1]; t.d[0] = head;
FESETENVD(fenv);
return t.f;
}
#endif