#include "math.h"
#include "fenv_private.h"
#include "fp_private.h"
static const double rintFactor = 6.7553994410557440000e+15;
static const double oneOverLn2 = 1.4426950408889633000e+00;
static const double ln2Head = 6.9314718055994530000e-01;
static const double ln2Tail = 2.3190468138462996000e-17;
static const double maxExp = 7.0978271289338397000e+02;
static const double denormal = 2.9387358770557188000e-39;
static const hexdouble infinity = HEXDOUBLE(0x7ff00000, 0x00000000);
static const double cm1 = 8.3333336309523691e-03;
static const double c0 = 4.1666668402777808000e-02;
static const double c1 = 1.6666666666666655000e-01;
static const double c2 = 4.9999999999999955000e-01;
extern const uint32_t expTable[];
struct expTableEntry
{
double x;
double f;
};
static const double kMinNormal = 0x1.0p-1022; static const double kMaxNormal = 1.7976931348623157e308;
#if !defined(BUILDING_FOR_CARBONCORE_LEGACY)
static const double minExp = -7.4513321910194122000e+02;
static const double oneOverDenorm = 3.402823669209384635e+38;
static const double cc4 = 0.5000000000000000000;
static const double cc3 = 0.16666666666663809522820;
static const double cc2 = 0.04166666666666111110939;
static const double cc1 = 0.008333338095239329810170;
static const double cc0 = 0.001388889583333492938381;
double exp ( double x )
{
hexdouble scale, xInHex, yInHex;
register double d, y, yTail, z, zTail, z2, temp1, temp2, power, result;
register int32_t i;
register struct expTableEntry *tablePointer, *pT;
register double FPR_oneOverLn2, FPR_rintFactor, FPR_ln2Head, FPR_ln2Tail, FPR_h, FPR_z, FPR_512, FPR_scale;
register double FPR_env, FPR_x, FPR_f, FPR_cm1, FPR_c0, FPR_c1, FPR_c2;
FPR_x = __FABS( x );
FPR_z = 0.0; FPR_512 = 512.0;
FPR_oneOverLn2 = oneOverLn2; FPR_rintFactor = rintFactor;
FPR_ln2Head = ln2Head; FPR_ln2Tail = ln2Tail;
tablePointer = ( struct expTableEntry * ) expTable + 177;
FEGETENVD ( FPR_env ); __ENSURE( FPR_z, FPR_512, FPR_oneOverLn2 ); __ENSURE( FPR_ln2Head, FPR_ln2Tail, FPR_rintFactor );
FESETENVD ( FPR_z );
FPR_h = __FMADD( x, FPR_oneOverLn2, FPR_rintFactor );
yInHex.d = FPR_h;
FPR_h -= FPR_rintFactor;
if (likely( FPR_x < FPR_512 ))
{
if (likely( x != FPR_z ))
{
scale.i.hi = ( yInHex.i.lo + 1023 ) << 20; scale.i.lo = 0;
y = __FNMSUB( FPR_ln2Head, FPR_h, x ); yTail = __FMUL( FPR_ln2Tail, FPR_h );
xInHex.d = __FMADD( FPR_512, y, FPR_rintFactor );
__NOOP;
__NOOP;
__NOOP;
i = xInHex.i.lo;
pT = &(tablePointer[i]);
FPR_x = pT->x; FPR_f = pT->f;
d = y - FPR_x;
z = d - yTail;
FPR_cm1 = cm1; FPR_c1 = c1;
FPR_c0 = c0; FPR_c2 = c2;
z2 = __FMUL( z, z ); zTail = d - z - yTail;
temp1 = __FMADD( FPR_cm1, z2, FPR_c1 ); temp2 = __FMADD( FPR_c0, z2, FPR_c2 );
temp1 = __FMADD( temp1, z, temp2 );
FPR_scale = scale.d;
temp2 = __FMADD( temp1, z2, z );
temp1 = __FMADD( zTail, temp2, zTail ); result = __FMUL( FPR_scale, FPR_f );
temp1 = temp1 + temp2;
result = __FMADD( result, temp1, result );
FESETENVD ( FPR_env );
__PROG_INEXACT( FPR_cm1 );
return result;
}
else
{
FESETENVD ( FPR_env );
return 1.0;
}
}
if (likely( ( x <= maxExp ) && ( x > minExp ) ))
{
if ( x >= FPR_512 )
{
power = 2.0;
scale.i.hi = ( yInHex.i.lo + 1022 ) << 20; scale.i.lo = 0;
}
else
{
power = denormal;
scale.i.hi = ( yInHex.i.lo + 1023+128 ) << 20; scale.i.lo = 0;
}
y = __FNMSUB( FPR_ln2Head, FPR_h, x ); yTail = __FMUL( FPR_ln2Tail, FPR_h );
FPR_scale = scale.d;
xInHex.d = __FMADD( FPR_512, y, FPR_rintFactor );
__NOOP;
__NOOP;
__NOOP;
i = xInHex.i.lo;
pT = &(tablePointer[i]);
FPR_x = pT->x; FPR_f = pT->f;
d = y - FPR_x; result = __FMUL( FPR_scale, FPR_f );
z = d - yTail;
FPR_cm1 = cm1; FPR_c1 = c1;
FPR_c0 = c0; FPR_c2 = c2;
z2 = __FMUL( z, z ); zTail = d - z - yTail;
temp1 = __FMADD( FPR_cm1, z2, FPR_c1 ); temp2 = __FMADD( FPR_c0, z2, FPR_c2 );
temp1 = __FMADD( temp1, z, temp2 );
temp2 = __FMADD( temp1, z2, z );
temp1 = __FMADD( zTail, temp2, zTail );
temp1 = temp1 + temp2;
result = __FMADD( result, temp1, result );
result = __FMUL( result, power );
FESETENVD ( FPR_env );
__PROG_INEXACT( FPR_oneOverLn2 );
return result;
}
FESETENVD ( FPR_env );
if ( x != x )
return x;
else if ( x == infinity.d )
return infinity.d;
else if ( x == -infinity.d )
return FPR_z;
else if ( x > maxExp )
{
__PROG_OF_INEXACT( kMaxNormal );
return infinity.d;
}
else
{
__PROG_UF_INEXACT( kMinNormal );
return FPR_z;
}
}
static const hexdouble k2M55 = HEXDOUBLE(0x3c800000, 0x00000000);
double expm1 ( double x )
{
hexdouble scale, invScale, xInHex, yInHex;
register double d, y, yTail, z, zTail, z2, temp1, temp2, power, result, f;
register int32_t i;
register struct expTableEntry *tablePointer, *pT;
register double FPR_oneOverLn2, FPR_rintFactor, FPR_ln2Head, FPR_ln2Tail, FPR_h, FPR_z, FPR_512, FPR_scale;
register double FPR_env, FPR_x, FPR_f, FPR_cc0, FPR_cc1, FPR_cc2, FPR_cc3, FPR_cc4, FPR_iscale;
FPR_x = __FABS( x );
FPR_z = 0.0; FPR_512 = 512.0;
FPR_oneOverLn2 = oneOverLn2; FPR_rintFactor = rintFactor;
FPR_ln2Head = ln2Head; FPR_ln2Tail = ln2Tail;
tablePointer = ( struct expTableEntry * ) expTable + 177;
FEGETENVD ( FPR_env ); __ENSURE( FPR_z, FPR_512, FPR_oneOverLn2 ); __ENSURE( FPR_ln2Head, FPR_ln2Tail, FPR_rintFactor );
FESETENVD ( FPR_z );
FPR_h = __FMADD( x, FPR_oneOverLn2, FPR_rintFactor );
yInHex.d = FPR_h;
FPR_h -= FPR_rintFactor;
if (likely( FPR_x < FPR_512 ))
{
if (unlikely( FPR_x < k2M55.d ))
{
FESETENVD ( FPR_env );
if ( x == FPR_z )
{
}
else
{
if ( FPR_x < kMinNormal )
__PROG_UF_INEXACT( kMinNormal );
else
__PROG_INEXACT( FPR_oneOverLn2 );
}
return x;
}
scale.i.hi = ( yInHex.i.lo + 1023 ) << 20; scale.i.lo = 0;
y = __FNMSUB( FPR_ln2Head, FPR_h, x ); yTail = __FMUL( FPR_ln2Tail, FPR_h );
xInHex.d = __FMADD( FPR_512, y, FPR_rintFactor );
invScale.i.hi = 0x7fe00000 - scale.i.hi; invScale.i.lo = 0;
i = xInHex.i.lo;
pT = &(tablePointer[i]);
FPR_x = pT->x; FPR_f = pT->f;
d = y - FPR_x;
FPR_cc4 = cc4;
z = d - yTail;
FPR_cc0 = cc0; FPR_cc2 = cc2;
FPR_cc1 = cc1; FPR_cc3 = cc3;
z2 = __FMUL( z, z ); zTail = d - z - yTail;
temp1 = __FMADD( FPR_cc0, z2, FPR_cc2 ); temp2 = __FMADD( FPR_cc1, z2, FPR_cc3 );
temp1 = __FMADD( temp1, z2, FPR_cc4 );
temp2 = __FMADD( temp2, z, temp1 );
FPR_iscale = invScale.d;
temp1 = __FMADD( temp2, z2, z );
temp2 = __FMADD( zTail, temp1, zTail );
temp2 = temp1 + temp2; d = FPR_f - FPR_iscale;
FPR_scale = scale.d;
temp1 = __FMADD( FPR_f, temp2, d );
result = __FMUL( temp1, FPR_scale );
FESETENVD ( FPR_env );
__PROG_INEXACT( FPR_oneOverLn2 );
return result;
}
if (likely( ( x <= maxExp ) && ( x > minExp ) ))
{
if ( x >= FPR_512 )
{
power = 2.0;
f = 0.5;
scale.i.hi = ( yInHex.i.lo + 1022 ) << 20; scale.i.lo = 0;
}
else
{
power = denormal;
f = oneOverDenorm;
scale.i.hi = ( yInHex.i.lo + 1023+128 ) << 20; scale.i.lo = 0;
if ( scale.i.hi < ( 168<<20 ) )
{
FESETENVD ( FPR_env );
__PROG_INEXACT( FPR_oneOverLn2 );
return -1.0;
}
}
y = __FNMSUB( FPR_ln2Head, FPR_h, x ); yTail = __FMUL( FPR_ln2Tail, FPR_h );
FPR_scale = scale.d;
xInHex.d = __FMADD( FPR_512, y, FPR_rintFactor );
invScale.i.hi = 0x7fe00000 - scale.i.hi; invScale.i.lo = 0;
i = xInHex.i.lo;
pT = &(tablePointer[i]);
FPR_x = pT->x; FPR_f = pT->f;
d = y - FPR_x;
FPR_iscale = invScale.d; FPR_cc4 = cc4;
z = d - yTail;
FPR_cc0 = cc0; FPR_cc2 = cc2;
FPR_cc1 = cc1; FPR_cc3 = cc3;
z2 = __FMUL( z, z ); zTail = d - z - yTail;
temp1 = __FMADD( FPR_cc0, z2, FPR_cc2 ); temp2 = __FMADD( FPR_cc1, z2, FPR_cc3 );
temp1 = __FMADD( temp1, z2, FPR_cc4 ); d = __FNMSUB( f, FPR_iscale, FPR_f );
temp2 = __FMADD( temp2, z, temp1 );
temp1 = __FMADD( temp2, z2, z );
temp2 = __FMADD( zTail, temp1, zTail );
temp2 = temp1 + temp2;
temp1 = __FMADD( FPR_f, temp2, d );
result = __FMUL( temp1, FPR_scale );
result = __FMUL( result, power );
FESETENVD ( FPR_env );
__PROG_INEXACT( FPR_oneOverLn2 );
return result;
}
FESETENVD ( FPR_env );
if ( x != x )
return x;
else if ( x == infinity.d )
return infinity.d;
else if ( x == -infinity.d )
return -1.0;
else if ( x > maxExp )
{
__PROG_OF_INEXACT( kMaxNormal );
return infinity.d;
}
else
{
#if 0
__PROG_UF_INEXACT( kMinNormal );
#else
__PROG_INEXACT( FPR_oneOverLn2 );
#endif
return -1.0;
}
}
#else
static const double maxExp2 = 1024.0;
static const double minNormExp2 = -1022.0;
static const double minExp2 = -1075.0;
double exp2 ( double x )
{
hexdouble OldEnvironment, scale, xInHex, yInHex;
register double d, y, yTail, z, zTail, z2, temp1, temp2, power, result;
register int32_t i;
register struct expTableEntry *tablePointer, *pT;
register double FPR_oneOverLn2, FPR_rintFactor, FPR_ln2Head, FPR_ln2Tail, FPR_h, FPR_z, FPR_512, FPR_scale;
register double FPR_env, FPR_diff, FPR_x, FPR_f, FPR_cm1, FPR_c0, FPR_c1, FPR_c2;
FPR_x = __FABS( x );
FPR_z = 0.0; FPR_512 = 512.0;
FPR_oneOverLn2 = oneOverLn2; FPR_rintFactor = rintFactor;
FPR_ln2Head = ln2Head; FPR_ln2Tail = ln2Tail;
tablePointer = ( struct expTableEntry * ) expTable + 177;
FEGETENVD ( FPR_env ); __ENSURE( FPR_z, FPR_512, FPR_oneOverLn2 ); __ENSURE( FPR_ln2Head, FPR_ln2Tail, FPR_rintFactor );
FESETENVD ( FPR_z );
FPR_h = x + FPR_rintFactor;
yInHex.d = FPR_h;
FPR_h -= FPR_rintFactor;
FPR_diff = x - FPR_h;
if (likely( FPR_x < FPR_512 ))
{
if (likely( FPR_x != FPR_z ))
{
scale.i.hi = ( yInHex.i.lo + 1023 ) << 20; scale.i.lo = 0;
y = __FMUL( FPR_ln2Head, FPR_diff ); yTail = __FMUL( FPR_ln2Tail, FPR_diff );
xInHex.d = __FMADD( FPR_512, y, FPR_rintFactor );
__NOOP;
__NOOP;
__NOOP;
i = xInHex.i.lo;
pT = &(tablePointer[i]);
FPR_x = pT->x; FPR_f = pT->f;
d = y - FPR_x;
z = d - yTail;
FPR_cm1 = cm1; FPR_c1 = c1;
FPR_c0 = c0; FPR_c2 = c2;
z2 = __FMUL( z, z ); zTail = d - z - yTail;
temp1 = __FMADD( FPR_cm1, z2, FPR_c1 ); temp2 = __FMADD( FPR_c0, z2, FPR_c2 );
temp1 = __FMADD( temp1, z, temp2 );
temp2 = __FMADD( temp1, z2, z );
FPR_scale = scale.d;
temp1 = __FMADD( zTail, temp2, zTail );
temp1 = temp1 + temp2; result = __FMUL( FPR_scale, FPR_f );
result = __FMADD( result, temp1, result );
FESETENVD ( FPR_env );
if ( FPR_diff != FPR_z)
__PROG_INEXACT( FPR_oneOverLn2 );
return result;
}
else
{
FESETENVD ( FPR_env );
return 1.0;
}
}
if (likely( ( x < maxExp2 ) && ( x > minExp2 ) ))
{
if ( x >= FPR_512 )
{
power = 2.0;
scale.i.hi = ( yInHex.i.lo + 1022 ) << 20; scale.i.lo = 0;
}
else
{
power = denormal;
scale.i.hi = ( yInHex.i.lo + 1023+128 ) << 20; scale.i.lo = 0;
}
FPR_diff = x - FPR_h;
FPR_scale = scale.d;
y = __FMUL( FPR_ln2Head, FPR_diff ); yTail = __FMUL( FPR_ln2Tail, FPR_diff );
xInHex.d = __FMADD( FPR_512, y, FPR_rintFactor );
__NOOP;
__NOOP;
__NOOP;
i = xInHex.i.lo;
pT = &(tablePointer[i]);
FPR_x = pT->x; FPR_f = pT->f;
d = y - FPR_x; result = __FMUL( FPR_scale, FPR_f );
z = d - yTail;
FPR_cm1 = cm1; FPR_c1 = c1;
FPR_c0 = c0; FPR_c2 = c2;
z2 = __FMUL( z, z ); zTail = d - z - yTail;
temp1 = __FMADD( FPR_cm1, z2, FPR_c1 ); temp2 = __FMADD( FPR_c0, z2, FPR_c2 );
temp1 = __FMADD( temp1, z, temp2 );
temp2 = __FMADD( temp1, z2, z );
temp1 = __FMADD( zTail, temp2, zTail );
temp1 = temp1 + temp2;
result = __FMADD( result, temp1, result );
result = __FMUL( result, power );
FESETENVD ( FPR_env );
if ( x < minNormExp2 )
{
OldEnvironment.d = FPR_env;
__NOOP;
__NOOP;
__NOOP;
OldEnvironment.i.lo |= FE_UNDERFLOW;
if ( FPR_h != x )
OldEnvironment.i.lo |= FE_INEXACT;
FESETENVD_GRP ( OldEnvironment.d );
}
else if ( FPR_h != x )
__PROG_INEXACT( FPR_oneOverLn2 );
return result;
}
FESETENVD ( FPR_env );
if ( x != x )
return x;
else if ( x == infinity.d )
return infinity.d;
else if ( x == -infinity.d )
return FPR_z;
else if ( x > maxExp )
{
__PROG_OF_INEXACT( kMaxNormal );
return infinity.d;
}
else
{
__PROG_UF_INEXACT( kMinNormal );
return FPR_z;
}
}
#endif