#include "math.h"
#include "fp_private.h"
#if (__WANT_LONG_DOUBLE_FORMAT__ - 0L == 128L)
#include "DD.h"
static const unsigned short sqrttable[256] = {
0x6b69, 0x6c68, 0x6e67, 0x6f65, 0x7064, 0x7263, 0x7361, 0x7460,
0x765f, 0x775d, 0x795c, 0x7a5b, 0x7b5a, 0x7d58, 0x7e57, 0x7f56,
0x8155, 0x8254, 0x8352, 0x8551, 0x8650, 0x874f, 0x894e, 0x8a4d,
0x8b4c, 0x8c4b, 0x8e4a, 0x8f48, 0x9047, 0x9246, 0x9345, 0x9444,
0x9543, 0x9742, 0x9841, 0x9940, 0x9a3f, 0x9c3e, 0x9d3d, 0x9e3c,
0x9f3c, 0xa13b, 0xa23a, 0xa339, 0xa438, 0xa637, 0xa736, 0xa835,
0xa934, 0xaa33, 0xac33, 0xad32, 0xae31, 0xaf30, 0xb02f, 0xb12e,
0xb32e, 0xb42d, 0xb52c, 0xb62b, 0xb72a, 0xb92a, 0xba29, 0xbb28,
0xbc27, 0xbd26, 0xbe26, 0xbf25, 0xc124, 0xc223, 0xc323, 0xc422,
0xc521, 0xc621, 0xc720, 0xc81f, 0xca1e, 0xcb1e, 0xcc1d, 0xcd1c,
0xce1c, 0xcf1b, 0xd01a, 0xd11a, 0xd219, 0xd418, 0xd518, 0xd617,
0xd716, 0xd816, 0xd915, 0xda14, 0xdb14, 0xdc13, 0xdd13, 0xde12,
0xdf11, 0xe111, 0xe210, 0xe310, 0xe40f, 0xe50e, 0xe60e, 0xe70d,
0xe80d, 0xe90c, 0xea0b, 0xeb0b, 0xec0a, 0xed0a, 0xee09, 0xef09,
0xf008, 0xf108, 0xf207, 0xf306, 0xf406, 0xf505, 0xf605, 0xf704,
0xf804, 0xf903, 0xfa03, 0xfb02, 0xfc02, 0xfd01, 0xfe01, 0xff00,
0x00ff, 0x01fd, 0x02fb, 0x03f9, 0x04f7, 0x05f5, 0x06f3, 0x07f2,
0x08f0, 0x09ee, 0x0aec, 0x0bea, 0x0ce9, 0x0de7, 0x0ee5, 0x0fe4,
0x10e2, 0x11e0, 0x12df, 0x13dd, 0x14db, 0x15da, 0x16d8, 0x17d7,
0x17d5, 0x18d4, 0x19d2, 0x1ad1, 0x1bcf, 0x1cce, 0x1dcc, 0x1ecb,
0x1fc9, 0x20c8, 0x20c6, 0x21c5, 0x22c4, 0x23c2, 0x24c1, 0x25c0,
0x26be, 0x27bd, 0x27bc, 0x28ba, 0x29b9, 0x2ab8, 0x2bb7, 0x2cb5,
0x2db4, 0x2db3, 0x2eb2, 0x2fb0, 0x30af, 0x31ae, 0x32ad, 0x33ac,
0x33aa, 0x34a9, 0x35a8, 0x36a7, 0x37a6, 0x37a5, 0x38a4, 0x39a3,
0x3aa2, 0x3ba0, 0x3c9f, 0x3c9e, 0x3d9d, 0x3e9c, 0x3f9b, 0x409a,
0x4099, 0x4198, 0x4297, 0x4396, 0x4495, 0x4494, 0x4593, 0x4692,
0x4791, 0x4890, 0x488f, 0x498e, 0x4a8d, 0x4b8c, 0x4b8c, 0x4c8b,
0x4d8a, 0x4e89, 0x4e88, 0x4f87, 0x5086, 0x5185, 0x5284, 0x5283,
0x5383, 0x5482, 0x5581, 0x5580, 0x567f, 0x577e, 0x587e, 0x587d,
0x597c, 0x5a7b, 0x5b7a, 0x5b79, 0x5c79, 0x5d78, 0x5d77, 0x5e76,
0x5f76, 0x6075, 0x6074, 0x6173, 0x6272, 0x6372, 0x6371, 0x6470,
0x656f, 0x656f, 0x666e, 0x676d, 0x686d, 0x686c, 0x696b, 0x6a6a
};
long double __sqrtldextra( double x, double y, double *pextra )
{
DblDbl zDD;
Double xD, guessD, recipD;
double scale, guess, recip, diff, recip2, epsilon, guessnew, gsq, guesslow;
double gsqlow, gg, gmid, gmidlow, recipnew, diff1, diff3, diff3x, diff3a;
double diff3ax, diff4, diff5, glow2, diff6, diff7, firstfix, fixup;
uint32_t ghead, reciphead, ival;
int index;
const Double scaleupD = {{0x5ff00000,0x0}};
const Double scaledownD = {{0x1ff00000,0x0}};
const Double adjustupD = {{0x4ff00000,0x0}};
const Double adjustdownD = {{0x2ff00000,0x0}};
xD.f = x;
scale = 1.0;
if (xD.i[0] < 0x07000000u) { x *= scaleupD.f; y *= scaleupD.f;
xD.f = x;
scale = adjustdownD.f; }
else if (xD.i[0] > 0x7fdfffffu) { x *= scaledownD.f; y *= scaledownD.f; xD.f = x;
scale = adjustupD.f; }
ghead = ((xD.i[0] + 0x3ff00000u) >> 1) & 0x7ff00000u;
reciphead = 0x7fc00000u - ghead;
index = (int)((xD.i[0] >> 13) & 0xffu);
ival = (uint32_t) sqrttable[index];
guessD.i[1] = 0u;
guessD.i[0] = ghead + ((ival & 0xff00u) << 4);
recipD.i[0] = reciphead + ((ival & 0xffu) << 12);
recipD.i[1] = 0u;
guess = guessD.f;
recip = recipD.f;
diff = __FNMSUB( guess, guess, x ); recip2 = recip + recip; guess = __FMADD( diff, recip, guess ); epsilon = __FNMSUB( recip, guess, 0.5 ); diff = __FNMSUB( guess, guess, x ); recip = __FMADD( epsilon, recip2, recip ); guess = __FMADD( recip, diff, guess ); recip2 = recip + recip;
epsilon = __FNMSUB( recip, guess, 0.5 ); diff = __FNMSUB( guess, guess, x ) + y; recip = __FMADD( epsilon, recip2, recip ); guessnew = __FMADD( recip, diff, guess ); recip2 = recip + recip;
gsq = __FMUL( guessnew, guessnew ); guesslow = __FMADD( recip, diff, (guess - guessnew) ); gsqlow = __FMSUB( guessnew, guessnew, gsq ); gg = guessnew + guessnew; gmid = __FMUL( gg, guesslow ); gmidlow = __FMSUB( gg, guesslow, gmid ); epsilon = __FNMSUB( recip, guessnew, 0.5 ); recip2 = recip + recip;
recipnew = __FMADD( epsilon, recip2, recip );
diff1 = x - gsq; diff3 = __FSUB( diff1, gmid ); diff3x = diff1 - diff3 - gmid;
diff3a = __FSUB( diff3, gsqlow ); diff3ax = diff3 - (diff3a + gsqlow); diff4 = diff3a + y;
diff5 = diff4 - gmidlow; glow2 = guesslow*guesslow; diff6 = __FSUB( diff5, glow2 ); diff7 = diff6 + (diff3x + diff3ax); firstfix = recipnew*diff7; fixup = __FADD( guesslow, firstfix ); guess = __FADD( guessnew, fixup ); fixup = guessnew - guess + fixup;
zDD.d[0] = guess*scale; zDD.d[1] = fixup*scale;
*pextra = ((guesslow - fixup) + firstfix)*scale;
return (zDD.f); }
long double sqrtl( long double x )
{
DblDbl xDD;
long double y;
double envD;
double extra;
uint32_t hibits;
xDD.f = x;
hibits = xDD.i[0];
if ((xDD.d[0] != 0.0) && (hibits < 0x7ff00000u)) {
FEGETENVD(envD);
FESETENVD(0.0);
y = __sqrtldextra(xDD.d[0],xDD.d[1],&extra);
FESETENVD(envD); return (y);
}
if ((hibits < 0x80000000u) || (xDD.d[0] == 0.0) || (xDD.d[0] != xDD.d[0]))
return (x);
xDD.i[0] = 0x7ff80000u; xDD.i[1] = 0u;
xDD.i[2] = 0x7ff80000u;
xDD.i[3] = 0u;
return (xDD.f); }
#include "math.h"
long double cbrtl( long double x )
{
DblDbl xDD;
double arg, guess;
long double ldguess, ldguess2, sqldguess;
const long double oneThird = 0.3333333333333333333333333333333L;
xDD.f = x;
arg = xDD.d[0];
if (arg != arg || arg == 0.0 || arg == INFINITY || -arg == INFINITY)
return x;
guess = cbrt(arg);
xDD.d[0] = guess;
xDD.d[1] = 0.0;
ldguess = xDD.f;
sqldguess = ldguess * ldguess;
ldguess2 = ldguess + ldguess;
return (ldguess2 + x/sqldguess)*oneThird;
}
#endif