#ifdef __APPLE_CC__
#if __APPLE_CC__ > 930
#include "math.h"
#include "fenv.h"
#include "fp_private.h"
#include "fenv_private.h"
#define POWER_NAN "37"
#define SET_INVALID 0x01000000
static const double ln2 = 6.931471805599452862e-1;
static const double ln2Tail = 2.319046813846299558e-17;
static const double kRint = 6.7553994410557440000e+15; static const double maxExp = 7.0978271289338397000e+02; static const double min2Exp = -7.4513321910194122000e+02; static const double denormal = 2.9387358770557188000e-39; static const double twoTo128 = 3.402823669209384635e+38; static const double twoTo52 = 4503599627370496.0; static const double oneOverLn2 = 1.4426950408889633000e+00; static const hexdouble huge = HEXDOUBLE(0x7ff00000, 0x00000000);
static const double c2 = -0.5000000000000000042725;
static const double c3 = 0.3333333333296328456505;
static const double c4 = -0.2499999999949632195759;
static const double c5 = 0.2000014541571685319141;
static const double c6 = -0.1666681511199968257150;
static const double d2 = -0.5000000000000000000000; static const double d3 = 0.3333333333333360968212; static const double d4 = -0.2500000000000056849516; static const double d5 = 0.1999999996377879912068; static const double d6 = -0.1666666662009609743117; static const double d7 = 0.1428690115662276259345; static const double d8 = -0.1250122079043555957999;
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;
extern unsigned long logTable[];
extern unsigned long expTable[];
struct logTableEntry
{
double X;
double G;
hexdouble F;
};
struct expTableEntry
{
double x;
double f;
};
static double PowerInner ( double, double );
static double _NearbyInt ( double x );
double pow ( double base, double exponent )
{
register long int isExpOddInt;
double absBase, result;
hexdouble u, v, OldEnvironment;
register unsigned long int expExp;
v.d = exponent;
u.d = base;
if ( ( ( v.i.hi & 0x000fffff ) | v.i.lo ) == 0 )
{ expExp = v.i.hi & 0xfff00000ul;
if ( expExp == 0x40000000 )
return base*base; else if ( exponent == 0.0 )
return 1.0; else if ( expExp == 0x3ff00000 )
return base; else if ( expExp == 0xbff00000 )
return 1.0/base; }
if ( ( v.i.hi & 0x7ff00000ul ) < 0x7ff00000ul )
{ if ( ( u.i.hi & 0x7ff00000ul ) < 0x7ff00000ul )
{ if ( base > 0.0 )
return PowerInner( base, exponent ); else if ( base < 0.0 )
{
if ( _NearbyInt ( exponent ) != exponent ) {
fegetenvd( OldEnvironment.d );
OldEnvironment.i.lo |= SET_INVALID;
result = nan ( POWER_NAN );
fesetenvd( OldEnvironment.d );
return result;
}
result = PowerInner( -base, exponent );
if ( _NearbyInt ( 0.5 * exponent ) != 0.5 * exponent ) result = - result;
return ( result );
}
else
{ isExpOddInt = ( ( _NearbyInt( exponent ) == exponent ) ?
( _NearbyInt( 0.5 * exponent ) != 0.5 * exponent ) : 0 );
if ( exponent > 0.0 )
return ( ( isExpOddInt ) ? base : 0.0 );
else return ( ( isExpOddInt ) ? 1.0/base : 1.0/__fabs( base ) );
}
}
else if (base != base)
return base;
else
{ if ( base > 0 )
return ( exponent > 0 ) ? huge.d : +0.0;
else
{
isExpOddInt = ( ( _NearbyInt ( exponent ) == exponent ) ?
( _NearbyInt ( 0.5 * exponent ) != 0.5 * exponent ) : 0 );
return ( exponent > 0 ) ?
( isExpOddInt ? -huge.d : +huge.d ) : ( isExpOddInt ? -0.0 : +0.0 );
}
}
}
else if ( exponent != exponent )
return base + exponent;
else
{ if ( base != base )
return base;
absBase = __fabs( base );
if ( ( exponent > 0 && absBase > 1 ) || ( exponent < 0 && absBase < 1 ) )
return huge.d;
else if ( ( exponent > 0 && absBase < 1 ) || ( exponent < 0 && absBase > 1 ) )
return +0.0;
else return 1.0;
}
}
static double _NearbyInt ( double x )
{
hexdouble xInHex;
double OldEnvironment;
xInHex.d = x;
if ( ( xInHex.i.hi & 0x7ffffffful ) < 0x43300000ul )
{ fegetenvd( OldEnvironment ); fesetenvd( 0.0 );
if ( xInHex.i.hi & 0x80000000ul ) x = ( ( x - twoTo52 ) + twoTo52 );
else
x = ( ( x + twoTo52 ) - twoTo52 );
fesetenvd( OldEnvironment );
return x;
}
else return ( x );
}
static double PowerInner ( double base, double exponent )
{
hexdouble u, v, OldEnvironment, scale, exp;
register long int i, n;
register unsigned long int head;
register double z, zLow, high, low, zSquared, temp1, temp2, temp3, result, tail,
d, x, xLow, y, yLow, power;
struct logTableEntry *logTablePointer = ( struct logTableEntry * ) logTable;
struct expTableEntry *expTablePointer = ( struct expTableEntry * ) expTable + 177;
u.d = base;
head = u.i.hi;
if ( head >= 0x000ffffful ) d = ( u.i.hi >> 20 ) - 1023.0;
else
{ u.d = base * twoTo128; d = ( u.i.hi >> 20 ) - 1151.0;
}
i = u.i.hi & 0x000fffff;
u.i.hi = i | 0x3ff00000;
if ( ( ( u.i.hi & 0x000fffff ) | u.i.lo ) == 0 ) {
z = d * exponent;
if ( z == _NearbyInt( z ) )
{
if (z > 2097)
n = 2098;
else if (z < -2098)
n = -2099;
else
n = z;
return scalbn( 1.0, n );
}
}
if ( u.i.hi & 0x00080000 )
{ n = 1;
u.d *= 0.5;
i = ( i >> 13 ) & 0x3f; }
else
{ n = 0;
i = ( i >> 12 ) + 64; }
fegetenvd( OldEnvironment.d ); fesetenvd( 0.0 );
z = ( u.d - logTablePointer[i].X ) * logTablePointer[i].G;
zLow = ( ( u.d - logTablePointer[i].X ) - z * logTablePointer[i].X ) * logTablePointer[i].G;
zSquared = z * z;
if ( n == 0 )
{ v.d = base;
if ( ( v.i.hi == 0x3ff00000 ) && ( v.i.lo == 0x0 ) )
{ fesetenvd( OldEnvironment.d );
return 1.0;
}
temp1 = d8 * zSquared + d6;
temp2 = d7 * zSquared + d5;
temp1 = temp1 * zSquared + d4;
temp2 = temp2 * zSquared + d3;
temp1 = temp1 * z + temp2;
temp2 = temp1 * z + d2;
temp3 = z + logTablePointer[i].F.d;
low = logTablePointer[i].F.d - temp3 + z + ( zLow - z * zLow );
result = ( temp2 * zSquared + low ) + temp3;
tail = temp3 - result + ( temp2 * zSquared + low );
}
else if ( logTablePointer[i].F.i.hi != 0 )
{ low = ln2Tail + zLow;
high = z + logTablePointer[i].F.d;
zLow = logTablePointer[i].F.d - high + z;
temp1 = ln2 + low;
low = ( ln2 - temp1 ) + low;
temp3 = c6 * zSquared + c4;
temp2 = c5 * zSquared + c3;
temp3 = temp3 * zSquared;
temp2 = ( temp2 * z + temp3 ) + c2;
temp3 = high + temp1;
temp1 = ( temp1 - temp3 ) + high;
result = ( ( ( temp2 * zSquared + low ) + temp1 ) + zLow ) + temp3;
tail = temp3 - result + zLow + temp1 + ( temp2 * zSquared + low );
}
else
{ low = ln2Tail + zLow;
temp3 = ln2 + low;
low = ( ln2 - temp3 ) + low;
temp1 = d8 * zSquared + d6;
temp2 = d7 * zSquared + d5;
temp1 = temp1 * zSquared + d4;
temp2 = temp2 * zSquared + d3;
temp1 = temp1 * z + temp2;
temp2 = temp1 * z + d2;
result = ( ( temp2 * zSquared + low ) + z ) + temp3;
tail = temp3 - result + z + ( temp2 * zSquared + low );
}
temp1 = d * ln2;
temp2 = __fmsub(d, ln2, temp1); z = temp1 + result;
zLow = temp1 - z + result + __fmadd(d, ln2Tail, tail + temp2); temp3 = exponent * zLow;
x = __fmadd(exponent, z, temp3); xLow = __fmsub(exponent, z, x) + temp3;
if ( base == _NearbyInt(base) && exponent > 0 && exponent == _NearbyInt(exponent))
{
}
else
OldEnvironment.i.lo |= FE_INEXACT;
u.d = x;
if ( ( u.i.hi & 0x7ff00000ul ) < 0x40800000ul )
{ scale.d = 0.0;
exp.d = x * oneOverLn2 + kRint;
scale.i.hi = ( exp.i.lo + 1023 ) << 20;
exp.d -= kRint;
y = x - ln2 * exp.d;
yLow = ln2Tail * exp.d;
u.d = 512.0 * y + kRint;
i = u.i.lo;
d = y - expTablePointer[i].x;
z = d - yLow;
zLow = d - z - yLow + xLow;
zSquared = z * z;
temp1 = cc0 * zSquared + cc2;
temp2 = cc1 * zSquared + cc3;
temp1 = temp1 * zSquared + cc4;
temp2 = temp1 + temp2 * z;
temp1 = temp2 * zSquared + z;
result = scale.d * expTablePointer[i].f;
#if 1
temp2 = __prod( result , ( temp1 + ( zLow + zLow * temp1 ) )); #else
temp2 = result * ( temp1 + ( zLow + zLow * temp1 ) );
#endif
result += temp2;
fesetenvd( OldEnvironment.d );
return result;
}
if ( ( x <= maxExp ) && ( x > min2Exp ) )
{
scale.d = 0.0;
exp.d = x * oneOverLn2 + kRint;
if ( x >= 512.0 )
{
power = 2.0;
scale.i.hi = ( exp.i.lo + 1023 - 1 ) << 20;
}
else
{
power = denormal;
scale.i.hi = ( exp.i.lo + 1023 + 128 ) << 20;
}
exp.d -= kRint;
y = x - ln2 * exp.d;
yLow = ln2Tail * exp.d;
u.d = 512.0 * y + kRint;
i = u.i.lo;
d = y - expTablePointer[i].x;
z = d - yLow;
zLow = d - z - yLow + xLow;
zSquared = z * z;
temp1 = cc0 * zSquared + cc2;
temp2 = cc1 * zSquared + cc3;
temp1 = temp1 * zSquared + cc4;
temp2 = temp1 + temp2 * z;
temp1 = temp2 * zSquared + z;
result = scale.d * expTablePointer[i].f;
temp2 = __prod( result , ( temp1 + ( zLow + zLow * temp1 ) ) ); u.d = ( result + temp2 ) * power;
}
else if ( x > maxExp )
u.d = huge.d;
else
u.d = +0.0;
i = u.i.hi & 0x7ff00000;
if ( i == 0x7ff00000 )
OldEnvironment.i.lo |= (FE_OVERFLOW | FE_INEXACT); else if ( i == 0x0 )
OldEnvironment.i.lo |= (FE_UNDERFLOW | FE_INEXACT);
fesetenvd( OldEnvironment.d );
return u.d;
}
#ifdef notdef
float powf ( float x, float y )
{
return (float)pow ( x, y );
}
#endif
#else
#error Version gcc-932 or higher required. Compilation terminated.
#endif
#endif