#include "fenv_private.h"
#include "fp_private.h"
#define TRIG_NAN "33"
struct tableEntry
{
double p;
double f5;
double f4;
double f3;
double f2;
double f1;
double f0;
};
extern const uint32_t tanatantable[];
static const double kPiBy2 = 1.570796326794896619231322; static const double kPiBy2Tail = 6.1232339957367660e-17; static const double k2ByPi = 0.636619772367581382; static const double kRint = 6.755399441055744e15; static const double kRintBig = 2.7021597764222976e16; static const double kPiScale42 = 1.38168706094305449e13; static const double kPiScale53 = 2.829695100811376e16; static const hexdouble infinity = HEXDOUBLE(0x7ff00000, 0x00000000);
static const double ts11 = 8.898406739539066157565e-3,
ts9 = 2.186936821731655951177e-2,
ts7 = 5.396825413618260185395e-2,
ts5 = 0.1333333333332513155016,
ts3 = 0.3333333333333333333333;
static const double tt7 = 0.05396875452473400572649,
tt5 = 0.1333333333304691192925,
tt3 = 0.3333333333333333357614;
static const double k2ToM1023 = 0x1.0p-1023; static const uint32_t indexTable[] =
{
16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26,
27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37,
38, 39, 40, 41, 42, 43, 44, 45, 47, 48, 49,
50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 68, 69, 70, 71, 72,
73, 74, 75, 76, 77, 78, 79, 81, 82, 83, 84,
85, 86, 87, 88, 89, 91, 92, 93, 94, 95, 96,
97, 98, 100, 101, 102, 103, 104, 105, 107, 108, 109,
110, 111, 113, 114, 115, 116, 117, 119, 120, 121, 122,
123, 125, 126, 127, 128, 130, 131, 132, 133, 135, 136,
137, 139, 140, 141, 142, 144, 145, 146, 148, 149, 150,
152, 153, 154, 156, 157, 159, 160, 161, 163, 164, 166,
167, 168, 170, 171, 173, 174, 176, 177, 179, 180, 182,
183, 185, 186, 188, 189, 191, 192, 194, 196, 197, 199,
200, 202, 204, 205, 207, 209, 210, 212, 214, 215, 217,
219, 220, 222, 224, 226, 228, 229, 231, 233, 235, 237,
238, 240, 242, 244, 246, 248, 250, 252, 254, 256, 256
};
static const double kMinNormal = 0x1.0p-1022; static const double kPiDiv4 = 0.7853980064392090732;
double tan ( double x )
{
register double absOfX, intquo ,arg, argtail, aSquared, aThird, aFourth,
temp1, temp2, s, u, v, w, y, result;
int tableIndex;
uint32_t iz;
hexdouble zD, aD, OldEnvironment;
struct tableEntry *tablePointer;
register double FPR_env, FPR_z, FPR_one, FPR_PiDiv4, FPR_256, FPR_piScale;
register double FPR_t, FPR_u, FPR_inf, FPR_pi53, FPR_2divPi, FPR_PiDiv2, FPR_PiDiv2Tail, FPR_kRintBig, FPR_kRint;
register uint32_t GPR_f;
register struct tableEntry *pT;
FEGETENVD( FPR_env );
FPR_z = 0.0; tablePointer = ( struct tableEntry * ) ( tanatantable - ( 16 * 14 ) );
absOfX = __FABS( x );
FPR_PiDiv4 = kPiDiv4; FPR_256 = 256.0;
FPR_kRint = kRint; FPR_one = 1.0;
FPR_PiDiv2 = kPiBy2;
__ENSURE( FPR_PiDiv2, FPR_PiDiv4, FPR_256);
__ENSURE( FPR_z, FPR_kRint, FPR_one);
FESETENVD( FPR_z );
if ( absOfX < FPR_PiDiv4 )
{
FPR_t = __FMADD( FPR_256, absOfX, FPR_kRint );
aD.d = FPR_t;
if (unlikely( absOfX == FPR_z ))
{
FESETENVD( FPR_env ); return x; }
if ( aD.i.lo < 16u ) {
register double FPR_ts3, FPR_ts5, FPR_ts7, FPR_ts9, FPR_ts11, FPR_Min;
aSquared = __FMUL( absOfX, absOfX ); FPR_Min = kMinNormal;
__NOOP;
FPR_ts11 = ts11; FPR_ts9 = ts9;
FPR_ts7 = ts7; FPR_ts5 = ts5;
FPR_ts3 = ts3; aThird = __FMUL( absOfX, aSquared );
temp1 = __FMADD( FPR_ts11, aSquared, FPR_ts9 );
temp1 = __FMADD( temp1, aSquared, FPR_ts7 );
temp1 = __FMADD( temp1, aSquared, FPR_ts5 );
temp1 = __FMADD( temp1, aSquared, FPR_ts3 );
if ( x > FPR_z )
result = __FMADD( temp1, aThird, absOfX );
else
{
aThird = -aThird;
result = __FMSUB( temp1, aThird, absOfX );
}
FESETENVD( FPR_env );
if ( __FABS ( result ) < FPR_Min )
__PROG_UF_INEXACT( FPR_Min );
else
__PROG_INEXACT( FPR_PiDiv2 );
return result;
}
else
{
register double FPR_tt3, FPR_tt5, FPR_tt7, FPR_sqr, FPR_f0;
tableIndex = indexTable[aD.i.lo - 16];
pT = &(tablePointer[tableIndex]);
FPR_f0 = pT->f0;
w = absOfX - FPR_f0; y = pT->p; FPR_tt3 = tt3;
FPR_tt7 = tt7; FPR_tt5 = tt5;
FPR_sqr = __FMUL( w, w );
temp1 = __FMADD( FPR_tt7, FPR_sqr, FPR_tt5 );
temp1 = __FMADD( temp1, FPR_sqr, FPR_tt3 );
temp1 = __FMADD( temp1, __FMUL( FPR_sqr, w ), w );
v = __FMUL(y, temp1 );
aSquared = __FMUL( v, v ); aThird = FPR_one + v;
aFourth = __FMUL( aSquared,aSquared );w = __FMADD( aThird, aSquared, aThird );
aThird = __FMADD( w, aFourth, w ); temp1 = __FMADD( v, y, temp1 );
if ( x > FPR_z ) result = __FMADD( temp1, aThird, y );
else
{
aThird = -aThird;
result = __FMSUB( temp1, aThird, y );
}
FESETENVD( FPR_env ); __PROG_INEXACT( FPR_PiDiv2 );
return result;
}
}
if (unlikely( x != x )) {
FESETENVD( FPR_env ); return ( x );
}
FPR_piScale = kPiScale42; FPR_PiDiv2Tail = kPiBy2Tail;
FPR_2divPi = k2ByPi; FPR_pi53 = kPiScale53;
FPR_inf = infinity.d; FPR_kRintBig = kRintBig;
__ENSURE( FPR_2divPi, FPR_PiDiv2, FPR_piScale );
if (unlikely( absOfX > FPR_piScale ))
{ if ( absOfX == FPR_inf )
{ OldEnvironment.d = FPR_env;
__NOOP;
__NOOP;
__NOOP;
OldEnvironment.i.lo |= SET_INVALID;
FESETENVD_GRP( OldEnvironment.d ); return ( nan ( TRIG_NAN ) ); }
while ( absOfX > FPR_pi53 )
{ intquo = __FMUL( x, FPR_2divPi ); FPR_t = __FNMSUB( intquo, FPR_PiDiv2, x );
x = __FNMSUB( intquo, FPR_PiDiv2Tail, FPR_t );
absOfX = __FABS ( x ) ;
}
FPR_t = __FMADD( x, FPR_2divPi, FPR_kRintBig );
intquo = FPR_t - FPR_kRintBig;
FPR_t = __FNMSUB( intquo, FPR_PiDiv2, x );
x = __FNMSUB( intquo, FPR_PiDiv2Tail, FPR_t );
absOfX = __FABS( x );
}
FPR_t = __FMADD( x, FPR_2divPi, FPR_kRint );
intquo = FPR_t - FPR_kRint;
zD.d = FPR_t;
FPR_t = __FNMSUB( intquo, FPR_PiDiv2, x );
arg = __FNMSUB( intquo, FPR_PiDiv2Tail, FPR_t );
absOfX = __FABS( arg ); FPR_t -= arg;
GPR_f = zD.i.lo;
FPR_u = __FMADD( FPR_256, absOfX, FPR_kRint );
aD.d = FPR_u;
__NOOP;
__NOOP;
argtail = __FNMSUB( intquo, FPR_PiDiv2Tail, FPR_t );
if ( aD.i.lo < 16u ) {
register double FPR_ts3, FPR_ts5, FPR_ts7, FPR_ts9, FPR_ts11;
aSquared = __FMUL( absOfX, absOfX );
u = absOfX;
v = argtail;
FPR_ts11 = ts11; FPR_ts9 = ts9;
FPR_ts7 = ts7; FPR_ts5 = ts5;
FPR_ts3 = ts3; iz = GPR_f & 1u;
temp1 = __FMADD( FPR_ts11, aSquared, FPR_ts9 );
temp1 = __FMADD( temp1, aSquared, FPR_ts7 );
temp1 = __FMADD( temp1, aSquared, FPR_ts5 );
temp1 = __FMADD( temp1, aSquared, FPR_ts3 );
aThird = __FMUL( aSquared, u );
if ( iz == 0 )
{
if ( arg > FPR_z )
{ w = __FMADD( temp1, aThird, v );
result = ( w + u );
}
else
{ w = __FNMSUB( temp1, aThird, v );
result = ( w - u );
}
FESETENVD( FPR_env ); __PROG_INEXACT( FPR_PiDiv2 );
return result;
}
else
{
register double FPR_k2ToM1023;
FPR_k2ToM1023 = k2ToM1023;
if ( arg > FPR_z )
{ s = __FMADD( temp1, aThird, v );
temp2 = s + u;
}
else
{ s = __FMSUB( temp1, aThird, v );
temp2 = s + u;
}
aFourth = ( u - temp2 ) + s;
if ( __FABS ( temp2 ) < k2ToM1023 )
{ if ( arg > FPR_z )
result = ( FPR_one / temp2 );
else
result = ( FPR_one / ( -temp2 ) );
FESETENVD( FPR_env ); __PROG_INEXACT( FPR_PiDiv2 );
return result;
}
u = FPR_one / temp2;
v = __FNMSUB( u, temp2, FPR_one );
temp1 = __FMSUB( aFourth, u, v );
if ( arg > FPR_z )
result = __FMSUB( temp1, u, u );
else
result = __FNMSUB( temp1, u, u );
FESETENVD( FPR_env ); __PROG_INEXACT( FPR_PiDiv2 );
return result;
}
}
tableIndex = indexTable [ aD.i.lo - 16 ];
if ( arg > FPR_z )
{ w = absOfX - tablePointer [tableIndex].f0 + argtail; aSquared = __FMUL( w, w );
v = absOfX - tablePointer [tableIndex].f0 - w + argtail; }
else
{ w = absOfX - tablePointer [tableIndex].f0 - argtail; aSquared = __FMUL( w, w );
v = absOfX - tablePointer [tableIndex].f0 - w - argtail; }
{
register double FPR_tt3, FPR_tt5, FPR_tt7;
FPR_t = v + w; FPR_u = __FMUL( aSquared, w );
y = tablePointer[tableIndex].p; FPR_tt3 = tt3;
FPR_tt7 = tt7; FPR_tt5 = tt5;
temp1 = __FMADD( FPR_tt7, aSquared, FPR_tt5 );
temp1 = __FMADD( temp1, aSquared, FPR_tt3 );
temp1 = __FMADD( temp1, FPR_u, FPR_t );
v = __FMUL( y, temp1 ); }
aSquared = __FMUL( v, v ); iz = GPR_f & 1u; aThird = FPR_one + v;
aFourth = __FMUL( aSquared, aSquared );
w = __FMADD( aThird, aSquared, aThird );
aThird = __FMADD( w, aFourth, w ); s = __FMUL( __FMADD( y, y, FPR_one ), aThird );
if ( iz == 0 )
{ if ( arg > FPR_z )
result = __FMADD( temp1, s, y );
else
{
y = -y;
result = __FNMSUB( temp1, s, y );
}
FESETENVD( FPR_env ); __PROG_INEXACT( FPR_PiDiv2 );
return result;
}
else
{ w = __FMADD( temp1, s, y );
aFourth = __FMADD( temp1, s, ( y - w ) );
u = FPR_one / w;
v = __FNMSUB( u, w, FPR_one);
w = __FMSUB( u, aFourth, v );
if ( arg > FPR_z )
result = __FMSUB( u, w, u );
else
result = __FNMSUB( u, w, u );
FESETENVD( FPR_env ); __PROG_INEXACT( FPR_PiDiv2 );
return result;
}
}