#include "math.h"
#include "fp_private.h"
#if (__WANT_LONG_DOUBLE_FORMAT__ - 0L == 128L)
#include "DD.h"
long double __sqrtldextra( double x, double y, double *pextra);
static const double kPiBy2 = 1.570796326794896619231322; static const double kPiBy2Tail1 = 6.123233995736766e-17; static const double kPiBy2Tail2 = -1.4973849048591698e-33;
static const double kPiBy4 = 7.85398163397448279e-01; static const double kPiBy4Tail1 = 3.061616997868383018e-17; static const double kPiBy4Tail2 = -7.486924524295849165e-34;
static const double k3PiBy4 = 2.356194490192344837; static const double k3PiBy4Tail1 = 9.184850993605148438e-17; static const double k3PiBy4Tail2 = 3.916898464750400322e-33;
static const double kPi = 3.141592653589793116; static const double kPiTail1 = 1.224646799147353207e-16; static const double kPiTail2 = -2.994769809718339666e-33;
static const double kSqrt3By2 = 0.866025403684439;
static const int kNterms = 23;
static const Float kNaNu = {0x7fc00000};
static const uint32_t asinCoeff[] = {
0x3FF00000, 0x00000000, 0x00000000, 0x00000000,
0x3FC55555, 0x55555555, 0x3C655555, 0x55555555,
0x3FB33333, 0x33333333, 0x3C499999, 0x9999999A,
0x3FA6DB6D, 0xB6DB6DB7, 0xBC324924, 0x92492492,
0x3F9F1C71, 0xC71C71C7, 0x3C1C71C7, 0x1C71C71C,
0x3F96E8BA, 0x2E8BA2E9, 0xBC31745D, 0x1745D174,
0x3F91C4EC, 0x4EC4EC4F, 0xBC2D89D8, 0x9D89D89E,
0x3F8C9999, 0x9999999A, 0xBC299999, 0x9999999A,
0x3F87A878, 0x78787878, 0x3C2E1E1E, 0x1E1E1E1E,
0x3F83FDE5, 0x0D79435E, 0x3C2435E5, 0x0D79435E,
0x3F812EF3, 0xCF3CF3CF, 0x3C1E79E7, 0x9E79E79E,
0x3F7DF3BD, 0x37A6F4DF, 0xBC190B21, 0x642C8591,
0x3F7A6863, 0xD70A3D71, 0xBC170A3D, 0x70A3D70A,
0x3F7782DD, 0xA12F684C, 0xBC02F684, 0xBDA12F68,
0x3F751BA3, 0x08D3DCB1, 0xBC1CB08D, 0x3DCB08D4,
0x3F731683, 0xBDEF7BDF, 0xBBE08421, 0x08421084,
0x3F715EE9, 0xD45D1746, 0xBC0745D1, 0x745D1746,
0x3F6FCAF8, 0xFB6DB6DB, 0x3C0B6DB6, 0xDB6DB6DB,
0x3F6D3D2A, 0x8E0DD67D, 0xBC0D67C8, 0xA60DD67D,
0x3F6B026F, 0x57B13B14, 0xBC03B13B, 0x13B13B14,
0x3F690CB7, 0x7F60C7CE, 0x3BD8F9C1, 0x8F9C18FA,
0x3F6750DE, 0x64D7D05F, 0x3C005F41, 0x7D05F418,
0x3F65C5F5, 0x6EFAAAAB, 0xBC055555, 0x55555555,
0x3F6464C0, 0x950F7D47, 0xBBF882B9, 0x31057262,
0x3F632755, 0x86C5F2F0, 0x3C04E5E0, 0xA72F0539,
0x3F6208D3, 0x570AE5A6, 0xBC069696, 0x96969697,
0x3F61052B, 0xC5FA960A, 0xBC05BC60, 0x9A90E7D9,
0x3F6018F9, 0x63C229BF, 0xBBF4F209, 0x4F2094F2,
0x3F5E82BE, 0x60D9127E, 0xBBEF7047, 0xDC11F704,
0x3F5CF7DE, 0xA5B6E830, 0xBBF15B1E, 0x5F75270D,
0x3F5B8D2E, 0x5667CE6C, 0x3BD0C971, 0x4FBCDA3B,
0x3F5A3F1E, 0xF82137EE, 0x3BEC71C7, 0x1C71C71C,
0x3F590A9F, 0x747DB95D, 0xBBE6A56A, 0x56A56A57,
0x3F57ED07, 0x9ED4C037, 0xBBD03D22, 0x6357E16F,
0x3F56E407, 0x90442038, 0x3BCA6F4D, 0xE9BD37A7,
0x3F55ED9A, 0x0BD901B6, 0xBBF040E6, 0xC2B4481D,
0x3F5507F9, 0x4C2470BD, 0x3BE56070, 0x381C0E07,
0x3F543195, 0xBF54E5D7, 0xBBFB9E4B, 0x17E4B17E,
0x3F53690E, 0x51F04536, 0x3BF9D974, 0x5D1745D1,
0x3F52AD29, 0xFCD49D54, 0x3BF82FBF, 0x309B8B57,
0x3F51FCD2, 0x5AE4A26E, 0x3BECA730, 0xFCD6E9E0,
0x3F51570F, 0x16ECE10A, 0x3BFA9E4B, 0xF3A9A378,
0x3F50BB02, 0x0BC1EAA0, 0x3BF8C4B0, 0x78787878,
0x3F5027E3, 0xF7FCD8BF, 0x3BF67F63, 0x2C234F73,
0x3F4F3A03, 0x591C2915, 0x3BE22B2B, 0xE05C0B81,
0x3F4E3373, 0x43F5C1F5, 0x3BEEF13D, 0x589D89D9,
0x3F4D3AF3, 0xC78CE2E4, 0x3BB8EEC6, 0x4A5294A5,
0x3F4C4F7C, 0x88F08B5E, 0x3BD4D455, 0x26BCA1AF,
0x3F4B701D, 0x9E1F038E, 0xBBE2BB23, 0xE09FAB8C,
0x3F4A9BFC, 0xD93A26FD, 0x3BC437EC, 0xD445D174,
0x3F49D253, 0x6C99619A, 0xBB9FD25B, 0x093F5DC8
};
struct asinCoeffEntry {
double Head;
double Tail;
} asinCoeffEntry;
long double asinl( long double x )
{
double arg, argl, argsq, argsq2;
double temp,temp2, temp3;
double sum, sumt, suml;
double xHead, xTail;
double sq, sq2, sq3;
double res, reslow, resmid, restiny, reshi, resbot, resinf;
double h, h2, g2, g3, g4, p, q, t;
double prod, prodlow;
double fenv;
DblDbl ddx;
int i;
uint32_t hibits;
struct asinCoeffEntry *asinPtr = (struct asinCoeffEntry *) asinCoeff;
FEGETENVD(fenv);
FESETENVD(0.0);
ddx.f = x;
xHead = ddx.d[0];
xTail = ddx.d[1];
hibits = ddx.i[0] & 0x7FFFFFFFu;
if (hibits > 0x3ff00000u) { if (xHead != xHead) { FESETENVD(fenv);
return x;
}
else { ddx.d[0] = kNaNu.f; ddx.d[1] = ddx.d[0];
FESETENVD(fenv);
return (ddx.f);
}
}
if ((xHead > 0.0) && (xHead -1.0 + xTail > 0.0)) { ddx.d[0] = kNaNu.f; ddx.d[1] = ddx.d[0];
FESETENVD(fenv);
return (ddx.f);
}
if ((xHead < 0.0) && (xHead +1.0 + xTail < 0.0)) { ddx.d[0] = kNaNu.f; ddx.d[1] = ddx.d[0];
FESETENVD(fenv);
return (ddx.f);
}
if (hibits < 0x3c800000u) { FESETENVD(fenv);
return x;
}
if (hibits <= 0x3fe00000u) {
temp = __FMUL( (2.0*xHead), xTail );
argsq = __FMADD( xHead, xHead, temp );
argsq2 = __FMSUB( xHead, xHead, argsq ) + temp + __FMSUB( (2.0*xHead), xTail, temp );
sum = asinPtr[50].Head; sum = __FMADD( sum, argsq, asinPtr[49].Head ); sum = __FMADD( sum, argsq, asinPtr[48].Head ); sum = __FMADD( sum, argsq, asinPtr[47].Head ); sum = __FMADD( sum, argsq, asinPtr[46].Head );
for (i = 45; i > kNterms+1; i--) sum = __FMADD( sum, argsq, asinPtr[i].Head );
sumt = __FMADD( sum, argsq, asinPtr[kNterms+1].Head ); suml = __FMADD( sum, argsq, asinPtr[kNterms+1].Head - sumt ); sum = sumt;
for (i = 1; i <= kNterms; i++) {
temp = __FMADD( sum, argsq, asinPtr[kNterms-i+1].Head );
suml = __FMADD( sum, argsq, asinPtr[kNterms-i+1].Head - temp ) +
__FMADD( suml, argsq, __FMADD( sum, argsq2, asinPtr[kNterms-i+1].Tail ) );
sum = temp;
}
prodlow = __FMADD( suml, argsq, sum*argsq2 ); prod = __FMADD( sum, argsq, prodlow ); prodlow = __FMSUB( sum, argsq, prod ) + prodlow; temp2 = __FMADD( prodlow, xHead, prod*xTail ); temp = __FMADD( prod, xHead, temp2 ); temp2 = __FMSUB( prod, xHead, temp ) + temp2;
res = __FADD( xHead, temp ); reslow = (xHead - res + temp); resmid = __FADD( xTail, temp2 );
restiny = xTail - resmid + temp2;
p = __FADD( reslow, resmid ); q = reslow - p + resmid; reshi = __FADD( res, p );
resbot = (res - reshi + p) + (q + restiny);
ddx.d[0] = reshi;
ddx.d[1] = resbot;
FESETENVD(fenv);
return (ddx.f);
}
else if (fabs(xHead) < kSqrt3By2) {
h = __FMUL( xHead, xHead ); h2 = __FMSUB( xHead, xHead, h ); g2 = __FMUL( (2.0*xHead), xTail );
g3 = __FMADD( xTail, xTail, __FMSUB( (2.0*xHead), xTail, g2 ) ); t = __FADD( h2 , g2 ); sq = __FADD( h, t );
if (fabs(h2) > fabs(g2)) g4 = h2 - t + g2 + g3; else g4 = g2 - t + h2 + g3; sq2 = __FADD( (h - sq + t), g4 ); sq3 = (h - sq + t) - sq2 + g4; temp = __FNMSUB( 2.0, sq, 1.0 ); temp2 = __FNMSUB( 2.0, sq2, __FNMSUB( 2.0, sq, 1.0 - temp ) ); arg = __FADD( temp, temp2 );
argl = __FNMSUB( 2.0, sq3, temp - arg + temp2 ); temp = 2.0*arg*argl;
argsq = __FMADD( arg, arg, temp ); argsq2 = __FMSUB( arg, arg, argsq ) + temp; sum = asinPtr[50].Head;
for (i = 49; i > kNterms+1; i--) sum = __FMADD( sum, argsq, asinPtr[i].Head );
sumt = __FMADD( sum, argsq, asinPtr[kNterms+1].Head ); suml = __FMADD( sum, argsq, asinPtr[kNterms+1].Head - sumt ); sum = sumt;
for (i = 1; i <= kNterms; i++) { temp = __FMADD( sum, argsq, asinPtr[kNterms-i+1].Head );
suml = __FMADD( sum, argsq, asinPtr[kNterms-i+1].Head - temp ) +
__FMADD( suml, argsq, __FMADD( sum, argsq2, asinPtr[kNterms-i+1].Tail ) );
sum = temp;
}
prodlow = __FMADD( suml, argsq, sum*argsq2 ); prod = __FMADD( sum, argsq, prodlow ); prodlow = __FMSUB( sum, argsq, prod ) + prodlow; temp2 = __FMADD( prodlow, arg, prod*argl ); temp = __FMADD( prod, arg, temp2 ); temp2 = __FMSUB( prod, arg, temp ) + temp2;
res = __FADD( arg, temp ); reslow = (arg - res + temp); resmid = __FADD( argl, temp2 );
restiny = argl - resmid + temp2;
p = __FADD( reslow, resmid ); q = reslow - p + resmid; reshi = __FADD( res, p );
resbot = __FADD( (res - reshi + p), (q + restiny) ); resinf = (res - reshi + p) - resbot + (q + restiny);
temp = kPiBy4;
temp2 = kPiBy4Tail1;
temp3 = kPiBy4Tail2;
res = __FNMSUB( 0.5, reshi, temp ); reslow = __FNMSUB( 0.5, reshi, temp - res ) ; resmid = __FNMSUB( 0.5, resbot, temp2 ); if (fabs(temp2) > fabs(0.5*resbot))
restiny = __FNMSUB( 0.5, resbot, temp2 - resmid ) + __FNMSUB( 0.5, resinf, temp3 );
else
restiny = temp2 - __FMADD( 0.5, resbot, resmid ) + __FNMSUB( 0.5, resinf, temp3 );
p = __FADD( reslow, resmid );
q = reslow - p + resmid;
reshi = __FADD( res, p );
resbot = (res - reshi + p) + (q + restiny);
if (xHead > 0.0) {
ddx.d[0] = reshi;
ddx.d[1] = resbot;
}
else {
ddx.d[0] = -reshi;
ddx.d[1] = -resbot;
}
FESETENVD(fenv);
return (ddx.f);
}
else {
if ((xHead - 1.0 + xTail) == 0.0) { ddx.d[0] = kPiBy2;
ddx.d[1] = kPiBy2Tail1;
FESETENVD(fenv);
return ddx.f;
}
if ((xHead + 1.0 - xTail) == 0.0) { ddx.d[0] = - kPiBy2;
ddx.d[1] = - kPiBy2Tail1;
FESETENVD(fenv);
return ddx.f;
}
if (xHead < 0.0) { temp = -xHead; temp2 = -xTail; }
else { temp = xHead; temp2 = xTail;
}
h = __FNMSUB( 0.5, temp, 0.5 ); argsq = __FNMSUB( 0.5, temp2, h ); argsq2 = __FNMSUB( 0.5, temp2, h - argsq );
ddx.f = __sqrtldextra(argsq, argsq2 , &temp);
sum = asinPtr[25].Head;
for (i = 24; i > kNterms/2+1; i--) sum = __FMADD( sum, argsq, asinPtr[i].Head );
sumt = __FMADD( sum, argsq, asinPtr[kNterms/2+1].Head ); suml = __FMADD( sum, argsq, asinPtr[kNterms/2+1].Head - sumt ); sum = sumt;
for (i = 1; i <= kNterms/2; i++) { temp = __FMADD( sum, argsq, asinPtr[kNterms/2-i+1].Head );
suml = __FMADD( sum, argsq, asinPtr[kNterms/2-i+1].Head - temp ) +
__FMADD( suml, argsq, __FMADD( sum, argsq2, asinPtr[kNterms/2-i+1].Tail ) );
sum = temp;
}
arg = ddx.d[0];
argl = ddx.d[1];
prodlow = __FMADD( suml, argsq, sum*argsq2 ); prod = __FMADD( sum, argsq, prodlow ); prodlow = __FMSUB( sum, argsq, prod ) + prodlow; temp2 = __FMADD( prodlow, arg, prod*argl ); temp = __FMADD( prod, arg, temp2 ); temp2 = __FMSUB( prod, arg, temp ) + temp2;
res = __FADD( arg, temp ); reslow = (arg - res + temp); resmid = __FADD( argl, temp2 );
restiny = argl - resmid + temp2;
p = __FADD( reslow, resmid ); q = reslow - p + resmid; reshi = __FADD( res, p );
resbot = __FADD( (res - reshi + p), (q + restiny) ); resinf = (res - reshi + p) - resbot + (q + restiny);
res = __FNMSUB( 2.0, reshi, kPiBy2 ); reslow = __FNMSUB( 2.0, reshi, kPiBy2 - res ); resmid = __FNMSUB( 2.0, resbot, kPiBy2Tail1 );
if (kPiBy2Tail1 > fabs(2.0*resbot))
restiny = __FNMSUB( 2.0, resbot, kPiBy2Tail1 - resmid ) + __FNMSUB( 2.0, resinf, kPiBy2Tail2 );
else
restiny = kPiBy2Tail1 - __FMADD( 2.0, resbot, resmid ) + __FNMSUB( 2.0, resinf, kPiBy2Tail2 );
p = __FADD( reslow, resmid );
q = reslow - p+resmid;
reshi = __FADD( res, p );
resbot = (res - reshi + p) + (q + restiny);
if (xHead < 0) {
ddx.d[0] = -reshi;
ddx.d[1] = -resbot;
}
else {
ddx.d[0] = reshi;
ddx.d[1] = resbot;
}
FESETENVD(fenv);
return (ddx.f);
}
}
long double acosl( long double x )
{
double arg, argl, argsq, argsq2;
double temp,temp2, temp3;
double sum, sumt, suml;
double xHead, xTail;
double sq, sq2, sq3;
double sqrtextra;
double res, reslow, resmid, restiny, reshi, resbot, resinf;
double h, h2, g2, g3, g4, p, q, t;
double prod, prodlow;
double fenv;
DblDbl ddx;
int i;
uint32_t hibits;
struct asinCoeffEntry *asinPtr = (struct asinCoeffEntry *) asinCoeff;
FEGETENVD(fenv);
FESETENVD(0.0);
ddx.f = x;
xHead = ddx.d[0];
xTail = ddx.d[1];
hibits = ddx.i[0] & 0x7FFFFFFFu;
if (hibits > 0x3ff00000u) { if (xHead != xHead) { FESETENVD(fenv);
return x;
}
else { ddx.d[0] = kNaNu.f; ddx.d[1] = ddx.d[0];
FESETENVD(fenv);
return (ddx.f);
}
}
if ((xHead > 0.0) && (xHead -1.0 + xTail > 0.0)) { ddx.d[0] = kNaNu.f; ddx.d[1] = ddx.d[0];
FESETENVD(fenv);
return (ddx.f);
}
if ((xHead < 0.0) && (xHead +1.0 + xTail < 0.0)) { ddx.d[0] = kNaNu.f; ddx.d[1] = ddx.d[0];
FESETENVD(fenv);
return (ddx.f);
}
if (hibits < 0x3c800000u) { res = kPiBy2;
reslow = -xHead;
resmid = __FSUB( kPiBy2Tail1, xTail );
restiny = kPiBy2Tail1 - resmid - xTail + kPiBy2Tail2;
p = __FADD( reslow, resmid );
reshi = __FADD( res, p );
if (fabs (reslow) > fabs (resmid))
q = reslow - p + resmid;
else
q = resmid - p + reslow;
resbot = (res-reshi+p) + (q+restiny);
ddx.d[0] = __FADD( reshi, resbot );
ddx.d[1] = (reshi - ddx.d[0]) + resbot;
FESETENVD(fenv);
return (ddx.f);
}
if (hibits <= 0x3fe00000u) { temp = __FMUL( (2.0*xHead), xTail ); argsq = __FMADD( xHead, xHead, temp );
argsq2 = __FMSUB( xHead, xHead, argsq ) + temp + __FMSUB( (2.0*xHead), xTail, temp );
sum = asinPtr[50].Head; sum = __FMADD( sum, argsq, asinPtr[49].Head ); sum = __FMADD( sum, argsq, asinPtr[48].Head ); sum = __FMADD( sum, argsq, asinPtr[47].Head ); sum = __FMADD( sum, argsq, asinPtr[46].Head );
for (i = 45; i >= kNterms+2; i--) sum = __FMADD( sum, argsq, asinPtr[i].Head );
sumt = __FMADD( sum, argsq, asinPtr[kNterms+1].Head ); suml = __FMADD( sum, argsq, asinPtr[kNterms+1].Head - sumt ); sum = sumt;
for (i = 1; i <= kNterms; i++) {
temp = __FMADD( sum, argsq, asinPtr[kNterms-i+1].Head );
suml = __FMADD( sum, argsq, asinPtr[kNterms-i+1].Head - temp ) +
__FMADD( suml, argsq, __FMADD( sum, argsq2, asinPtr[kNterms-i+1].Tail ) );
sum = temp;
}
prodlow = __FMADD( suml, argsq, sum*argsq2 ); prod = __FMADD( sum, argsq, prodlow ); prodlow = __FMSUB( sum, argsq, prod ) + prodlow; temp2 = __FMADD( prodlow, xHead, prod*xTail ); temp = __FMADD( prod, xHead, temp2 ); temp2 = __FMSUB( prod, xHead, temp ) + temp2;
res = __FADD( xHead, temp ); reslow = (xHead - res + temp);
resmid = __FADD( xTail, temp2 );
restiny = xTail - resmid + temp2;
p = __FADD( reslow, resmid ); q = reslow - p + resmid; reshi = __FADD( res, p );
resbot = __FADD( (res - reshi + p), (q + restiny) );
resinf = (res - reshi + p) - resbot + (q + restiny);
res = __FSUB( kPiBy2, reshi ); reslow = kPiBy2 - res - reshi;
resmid = __FSUB( kPiBy2Tail1, resbot );
if (kPiBy2Tail1 > fabs (resbot))
restiny = kPiBy2Tail1 - resmid - resbot + (kPiBy2Tail2 - resinf);
else
restiny = kPiBy2Tail1 - (resmid + resbot) + (kPiBy2Tail2 - resinf);
p = __FADD( reslow, resmid );
q = reslow - p + resmid;
reshi = __FADD( res, p );
resbot = (res - reshi + p) + (q + restiny);
ddx.d[0] = __FADD( reshi, resbot );
ddx.d[1] = (reshi - ddx.d[0]) + resbot;
FESETENVD(fenv);
return (ddx.f);
}
if (fabs (xHead) < kSqrt3By2) { h = __FMUL( xHead, xHead ); h2 = __FMSUB( xHead, xHead, h ); g2 = __FMUL( (2.0*xHead), xTail );
g3 = __FMADD( xTail, xTail, __FMSUB( (2.0*xHead), xTail, g2 ) ); t = __FADD( h2 , g2 ); sq = __FADD( h, t );
if (fabs (h2) > fabs (g2)) g4 = h2 - t + g2 + g3; else g4 = g2 - t + h2 + g3; sq2 = __FADD( (h - sq + t), g4 ); sq3 = (h - sq + t) - sq2 + g4; temp = __FMSUB( 2.0, sq, 1.0 ); temp2 = __FMADD( 2.0, sq2, __FMSUB( 2.0, sq, 1.0 + temp ) ); arg = __FADD( temp, temp2 );
argl = __FMADD( 2.0, sq3, temp - arg + temp2 ); temp = 2.0*arg*argl;
argsq = __FMADD( arg, arg, temp ); argsq2 = __FMSUB( arg, arg, argsq ) + temp; sum = asinPtr[50].Head;
for (i = 49;i >= kNterms+2;i--) sum = __FMADD( sum, argsq, asinPtr[i].Head );
sumt = __FMADD( sum, argsq, asinPtr[kNterms+1].Head ); suml = __FMADD( sum, argsq, asinPtr[kNterms+1].Head - sumt ); sum = sumt;
for (i = 1;i <= kNterms;i++) { temp = __FMADD( sum, argsq, asinPtr[kNterms-i+1].Head );
suml = __FMADD( sum, argsq, asinPtr[kNterms-i+1].Head - temp ) +
__FMADD( suml, argsq, __FMADD( sum, argsq2, asinPtr[kNterms-i+1].Tail ) );
sum = temp;
}
prodlow = __FMADD( suml, argsq, sum*argsq2 ); prod = __FMADD( sum, argsq, prodlow ); prodlow = __FMSUB( sum, argsq, prod ) + prodlow; temp2 = __FMADD( prodlow, arg, prod*argl ); temp = __FMADD( prod, arg, temp2 ); temp2 = __FMSUB( prod, arg, temp ) + temp2;
res = __FADD( arg, temp ); reslow = (arg - res + temp); resmid = __FADD( argl, temp2 );
restiny = argl - resmid + temp2;
p = __FADD( reslow, resmid ); q = reslow - p + resmid; reshi = __FADD( res, p );
resbot = __FADD( (res - reshi + p), (q + restiny) ); resinf = (res - reshi + p) - resbot + (q + restiny);
if (xHead>0.0) {
temp = kPiBy4; temp2 = kPiBy4Tail1; temp3 = kPiBy4Tail2;
}
else {
temp = k3PiBy4; temp2 = k3PiBy4Tail1; temp3 = k3PiBy4Tail2; reshi = -reshi;
resbot = -resbot;
resinf = -resinf;
}
res = __FNMSUB( 0.5, reshi, temp ); reslow = __FNMSUB( 0.5, reshi, temp - res ) ; resmid = __FNMSUB( 0.5, resbot, temp2 );
if (fabs (temp2) > fabs (0.5*resbot)) {
restiny = __FNMSUB( 0.5, resbot, temp2 - resmid ) + __FNMSUB( 0.5, resinf, temp3 );
}
else {
restiny = temp2 - __FMADD( 0.5, resbot, resmid ) + __FNMSUB( 0.5, resinf, temp3 );
}
p = __FADD( reslow, resmid );
q = reslow - p + resmid;
reshi = __FADD( res, p );
resbot = (res - reshi + p) + (q + restiny);
ddx.d[0] = __FADD( reshi, resbot );
ddx.d[1] = (reshi - ddx.d[0]) + resbot;
FESETENVD(fenv);
return (ddx.f);
}
else {
if ((xHead - 1.0 + xTail) == 0.0) { FESETENVD(fenv);
return (long double) 0.0;
}
if ((xHead + 1.0 - xTail) == 0.0) { ddx.d[0] = kPi;
ddx.d[1] = kPiTail1;
FESETENVD(fenv);
return ddx.f;
}
if (xHead < 0.0) { temp = -xHead; temp2 = -xTail; }
else { temp = xHead; temp2 = xTail; }
h = __FNMSUB( 0.5, temp, 0.5 ); argsq = __FNMSUB( 0.5, temp2, h ); argsq2 = __FNMSUB( 0.5, temp2, h - argsq );
ddx.f = __sqrtldextra(argsq , argsq2, &sqrtextra);
sum = asinPtr[25].Head;
for (i = 24; i >= kNterms/2+2; i--)
sum = __FMADD( sum, argsq, asinPtr[i].Head );
sumt = __FMADD( sum, argsq, asinPtr[kNterms/2+1].Head ); suml = __FMADD( sum, argsq, asinPtr[kNterms/2+1].Head - sumt ); sum = sumt;
for (i = 1; i <= kNterms/2; i++) { temp = __FMADD( sum, argsq, asinPtr[kNterms/2-i+1].Head );
suml = __FMADD( sum, argsq, asinPtr[kNterms/2-i+1].Head - temp ) +
__FMADD( suml, argsq, __FMADD( sum, argsq2, asinPtr[kNterms/2-i+1].Tail ) );
sum = temp;
}
arg = ddx.d[0];
argl = ddx.d[1];
prodlow = __FMADD( suml, argsq, sum*argsq2 ); prod = __FMADD( sum, argsq, prodlow ); prodlow = __FMSUB( sum, argsq, prod ) + prodlow; temp2 = __FMADD( prodlow, arg, prod*argl ); temp = __FMADD( prod, arg, temp2 ); temp2 = __FMSUB( prod, arg, temp ) + temp2;
res = __FADD( arg, temp ); reslow = (arg - res + temp); resmid = __FADD( argl, temp2 );
restiny = argl - resmid + temp2 + sqrtextra;
p = __FADD( reslow, resmid ); q = reslow - p + resmid; reshi = __FADD( res, p );
resbot = __FADD( (res - reshi + p), (q + restiny) );
if (xHead > 0) { ddx.d[0] = __FMADD( 2.0, reshi, 2.0*resbot ); ddx.d[1] = __FMADD( 2.0, resbot, __FMSUB( 2.0, reshi, ddx.d[0] ) ); FESETENVD(fenv);
return (ddx.f);
}
resinf = (res - reshi +p ) - resbot + (q + restiny); res = __FNMSUB( 2.0, reshi, kPi ); reslow = __FNMSUB( 2.0, reshi, kPi - res ); resmid = __FNMSUB( 2.0, resbot, kPiTail1 );
if (kPiTail1 > fabs(2.0*resbot))
restiny = __FNMSUB( 2.0, resbot, kPiTail1 - resmid ) + __FNMSUB( 2.0, resinf, kPiTail2 );
else
restiny = kPiTail1 - __FMADD( 2.0, resbot, resmid ) + __FNMSUB( 2.0, resinf, kPiTail2 );
p = __FADD( reslow, resmid );
q = reslow - p + resmid;
reshi = __FADD( res, p );
resbot = (res - reshi + p) + (q + restiny);
ddx.d[0] = __FADD( reshi, resbot );
ddx.d[1] = (reshi - ddx.d[0]) + resbot;
FESETENVD(fenv);
return (ddx.f);
}
}
#endif