#include "math.h"
#include "fp_private.h"
#if (__WANT_LONG_DOUBLE_FORMAT__ - 0L == 128L)
#include "DD.h"
extern long double rintl(long double);
static const Double kInfinityu = {{0x7ff00000,0x0}};
long double _ExpInnerLD(double, double, double, double *, int);
long double _LogInnerLD(double, double, double, double *, int);
static long double PowerInnerLD(double, double, double, double);
long double powl( long double base, long double exponent )
{
DblDbl u, v;
double basehi, baselo, exphi, explo;
int isExpOddInt;
uint32_t expExp;
double fpenv;
long double ldTemp;
u.f = base;
basehi = u.d[0];
baselo = u.d[1];
v.f = exponent;
exphi = v.d[0];
explo = v.d[1];
if (((v.i[0] & 0x000fffffu) | v.i[1] | (v.i[2] & 0x7fffffffu) | v.i[3]) == 0) {
expExp = v.i[0] & 0xfff00000u;
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; }
FEGETENVD(fpenv);
FESETENVD(0.0);
if ((v.i[0] & 0x7ff00000u) < 0x7ff00000u) { if ((u.i[0] & 0x7ff00000u) < 0x7ff00000u) { if (basehi > 0.0) { ldTemp = PowerInnerLD(basehi, baselo, exphi, explo);
FESETENVD(fpenv);
return ldTemp;
}
if (basehi < 0.0) { if (rintl(exponent) != exponent) { u.d[0] = NAN;
u.d[1] = 0.0;
FESETENVD(fpenv);
return u.f;
}
u.f = PowerInnerLD(-basehi, -baselo, exphi, explo);
if (rintl(0.5*exponent) != 0.5*exponent) u.f = -u.f;
FESETENVD(fpenv);
return u.f;
}
else { isExpOddInt = ((rintl(exponent) == exponent) ?
(rintl(0.5*exponent) != 0.5*exponent) : 0);
if (exphi > 0.0) {
ldTemp = (isExpOddInt) ? base : 0.0;
FESETENVD(fpenv);
return ldTemp;
}
else { u.d[0] = (isExpOddInt) ? 1.0/basehi : 1.0/fabs(basehi);
u.d[1] = 0.0;
FESETENVD(fpenv);
return u.f;
}
}
} else if (basehi != basehi) return base;
else {
if (basehi > 0.0) {
if (exphi > 0.0) return (long double) INFINITY;
else return 0.0L;
}
else {
isExpOddInt = ((rintl(exponent) == exponent) && (rintl(0.5*exponent) != 0.5*exponent));
if (exphi > 0.0) {
ldTemp = (isExpOddInt) ? (long double) -INFINITY : (long double) INFINITY;
FESETENVD(fpenv);
return ldTemp;
}
else { u.d[0] = (isExpOddInt) ? -0.0 : 0.0;
u.d[1] = 0.0;
FESETENVD(fpenv);
return u.f;
}
}
}
}
else if (exphi != exphi) {
if (base == 1.0L)
return base;
else
return base + exponent;
}
else {
if (basehi != basehi)
return base;
ldTemp = fabsl(base);
if ( (exphi > 0.0 && ldTemp > 1.0L) || (exphi < 0.0 && ldTemp < 1.0L) )
return (long double) INFINITY;
else if ( (exphi > 0.0 && ldTemp < 1.0L) || (exphi < 0.0 && ldTemp > 1.0L) )
return 0.0L;
else
return 1.0L;
}
}
static long double PowerInnerLD(double alpha, double beta, double gamma, double delta)
{
DblDbl u;
double extra, ahi, amid, amid2, alow, amid3;
double amid4, ahi1, amid5, amid6, ahi2, amid7, amid8;
u.f = _LogInnerLD(alpha, beta, 0.0, &extra, 1);
ahi = __FMUL( u.d[0], gamma );
amid = __FMSUB( u.d[0], gamma, ahi ); amid2 = __FMUL( u.d[1], gamma );
alow = __FMADD( gamma, extra, __FMSUB( u.d[1], gamma, amid2 ) ); amid3 = __FMUL( u.d[0], delta );
alow = __FMADD( delta, u.d[1], __FMSUB( u.d[0], delta, amid3) + alow ); amid4 = __FADD( amid, amid2 );
ahi1 = __FADD( ahi, amid4 );
amid5 = ahi - ahi1 + amid4;
if (fabs(amid) > fabs(amid2))
alow = amid - amid4 + amid2 + alow;
else
alow = amid2 - amid4 + amid + alow;
amid6 = __FADD( amid3, amid5 );
ahi2 = __FADD( ahi1, amid6 );
if (fabs(amid3) > fabs(amid5))
alow = amid3 - amid6 + amid5 + alow;
else
alow = amid5 - amid6 + amid3 + alow;
amid7 = ahi1 - ahi2 + amid6;
amid8 = __FADD( amid7, alow );
alow = amid7 - amid8 + alow;
if (fabs(ahi) == kInfinityu.f) {
if (ahi > 0.0) {
u.d[0] = kInfinityu.f;
u.d[1] = 0.0;
}
else {
u.d[0] = 0.0;
u.d[1] = 0.0;
}
}
else {
u.f = _ExpInnerLD(ahi2, amid8, alow, &extra, 3);
if ((u.i[0] & 0x7ff00000) == 0x7ff00000)
u.d[1] = 0.0;
}
return u.f;
}
#endif