#include "math.h"
#include <emmintrin.h>
#include <stdint.h>
#include "xmmLibm_Prefix.h"
#include <float.h>
#if defined( __SSE3__ )
#include <pmmintrin.h>
#endif
static const __m128d log2polynomialCoefficients[6] = {
{ -3.86011461177189878, 11.2166692423937703 },
{-19.7108996153552616, 27.5580934192221441 },
{-28.7390040495347408, 22.2562416071123559 },
{-12.7549704819157692, 5.34170406305677394 },
{-1.59009195684350404, 0.318880468456163145 },
{-0.0386477049377158355, 0.00213962027817084459 }
};
static const double exp2polynomialCoefficients[3] = { 1.00000000011116223,
0.693146680445468386,
0.240559810709772152 };
typedef uint32_t vUInt32 __attribute__ ((__vector_size__ (16)));
static const vUInt32 v1023 = { 1023, 0, 0, 0 };
static const __m128d exponentMask = { __builtin_inf(), __builtin_inf() };
static const __m128d half_sd = { 0.5, 0.0 };
static const __m128d one_sd = { 1.0, 0.0 };
static const __m128d two_pd = { 2.0, 2.0 };
static const __m128i mantissaHiBits = { 0x000FF00000000000ULL, 0ULL };
static const __m128d exp2_table[256] = {
{ 1.000000000000000000000, 0x0.00p0 }, { 1.002711275050202521797, 0x0.01p0 }, { 1.005429901112802726360, 0x0.02p0 }, { 1.008155898118417548304, 0x0.03p0 },
{ 1.010889286051700475255, 0x0.04p0 }, { 1.013630084951489429557, 0x0.05p0 }, { 1.016378314910953095662, 0x0.06p0 }, { 1.019133996077737913666, 0x0.07p0 },
{ 1.021897148654116627142, 0x0.08p0 }, { 1.024667792897135720764, 0x0.09p0 }, { 1.027445949118763746100, 0x0.0ap0 }, { 1.030231637686040979673, 0x0.0bp0 },
{ 1.033024879021228414899, 0x0.0cp0 }, { 1.035825693601957198098, 0x0.0dp0 }, { 1.038634101961378730650, 0x0.0ep0 }, { 1.041450124688316103416, 0x0.0fp0 },
{ 1.044273782427413754803, 0x0.10p0 }, { 1.047105095879289793359, 0x0.11p0 }, { 1.049944085800687210153, 0x0.12p0 }, { 1.052790773004626423415, 0x0.13p0 },
{ 1.055645178360557157049, 0x0.14p0 }, { 1.058507322794512761632, 0x0.15p0 }, { 1.061377227289262092924, 0x0.16p0 }, { 1.064254912884464499001, 0x0.17p0 },
{ 1.067140400676823697168, 0x0.18p0 }, { 1.070033711820241872914, 0x0.19p0 }, { 1.072934867525975555225, 0x0.1ap0 }, { 1.075843889062791047806, 0x0.1bp0 },
{ 1.078760797757119860307, 0x0.1cp0 }, { 1.081685614993215249768, 0x0.1dp0 }, { 1.084618362213309206155, 0x0.1ep0 }, { 1.087559060917769659937, 0x0.1fp0 },
{ 1.090507732665257689675, 0x0.20p0 }, { 1.093464399072885839814, 0x0.21p0 }, { 1.096429081816376882585, 0x0.22p0 }, { 1.099401802630221913759, 0x0.23p0 },
{ 1.102382583307840890896, 0x0.24p0 }, { 1.105371445701741173195, 0x0.25p0 }, { 1.108368411723678725878, 0x0.26p0 }, { 1.111373503344817548211, 0x0.27p0 },
{ 1.114386742595892432206, 0x0.28p0 }, { 1.117408151567369278823, 0x0.29p0 }, { 1.120437752409606746440, 0x0.2ap0 }, { 1.123475567333019897731, 0x0.2bp0 },
{ 1.126521618608241848136, 0x0.2cp0 }, { 1.129575928566288078869, 0x0.2dp0 }, { 1.132638519598719195614, 0x0.2ep0 }, { 1.135709414157805463574, 0x0.2fp0 },
{ 1.138788634756691564576, 0x0.30p0 }, { 1.141876203969561576201, 0x0.31p0 }, { 1.144972144431804172982, 0x0.32p0 }, { 1.148076478840178937801, 0x0.33p0 },
{ 1.151189229952982673311, 0x0.34p0 }, { 1.154310420590215935377, 0x0.35p0 }, { 1.157440073633751120852, 0x0.36p0 }, { 1.160578212027498778980, 0x0.37p0 },
{ 1.163724858777577475522, 0x0.38p0 }, { 1.166880036952481658474, 0x0.39p0 }, { 1.170043769683250189928, 0x0.3ap0 }, { 1.173216080163637320410, 0x0.3bp0 },
{ 1.176396991650281220743, 0x0.3cp0 }, { 1.179586527462875844563, 0x0.3dp0 }, { 1.182784710984341014495, 0x0.3ep0 }, { 1.185991565660993840581, 0x0.3fp0 },
{ 1.189207115002721026897, 0x0.40p0 }, { 1.192431382583151178167, 0x0.41p0 }, { 1.195664392039827328418, 0x0.42p0 }, { 1.198906167074380357818, 0x0.43p0 },
{ 1.202156731452703075647, 0x0.44p0 }, { 1.205416109005123859177, 0x0.45p0 }, { 1.208684323626581624822, 0x0.46p0 }, { 1.211961399276801243374, 0x0.47p0 },
{ 1.215247359980468955243, 0x0.48p0 }, { 1.218542229827408451825, 0x0.49p0 }, { 1.221846032972757400969, 0x0.4ap0 }, { 1.225158793637145526745, 0x0.4bp0 },
{ 1.228480536106870024682, 0x0.4cp0 }, { 1.231811284734075861991, 0x0.4dp0 }, { 1.235151063936933191201, 0x0.4ep0 }, { 1.238499898199816540156, 0x0.4fp0 },
{ 1.241857812073484002013, 0x0.50p0 }, { 1.245224830175257979548, 0x0.51p0 }, { 1.248600977189204819240, 0x0.52p0 }, { 1.251986277866316221719, 0x0.53p0 },
{ 1.255380757024691096291, 0x0.54p0 }, { 1.258784439549716527296, 0x0.55p0 }, { 1.262197350394250738859, 0x0.56p0 }, { 1.265619514578806281691, 0x0.57p0 },
{ 1.269050957191733219886, 0x0.58p0 }, { 1.272491703389402761815, 0x0.59p0 }, { 1.275941778396392001227, 0x0.5ap0 }, { 1.279401207505669324505, 0x0.5bp0 },
{ 1.282870016078778263591, 0x0.5cp0 }, { 1.286348229546025567771, 0x0.5dp0 }, { 1.289835873406665944785, 0x0.5ep0 }, { 1.293332973229089466471, 0x0.5fp0 },
{ 1.296839554651009640551, 0x0.60p0 }, { 1.300355643379650594227, 0x0.61p0 }, { 1.303881265191935812098, 0x0.62p0 }, { 1.307416445934677318164, 0x0.63p0 },
{ 1.310961211524764413738, 0x0.64p0 }, { 1.314515587949354635811, 0x0.65p0 }, { 1.318079601266064049270, 0x0.66p0 }, { 1.321653277603157539133, 0x0.67p0 },
{ 1.325236643159741323217, 0x0.68p0 }, { 1.328829724205954354588, 0x0.69p0 }, { 1.332432547083161500368, 0x0.6ap0 }, { 1.336045138204145832361, 0x0.6bp0 },
{ 1.339667524053302916087, 0x0.6cp0 }, { 1.343299731186835321850, 0x0.6dp0 }, { 1.346941786232945803548, 0x0.6ep0 }, { 1.350593715892034252235, 0x0.6fp0 },
{ 1.354255546936892651289, 0x0.70p0 }, { 1.357927306212901141791, 0x0.71p0 }, { 1.361609020638224754052, 0x0.72p0 }, { 1.365300717204011915484, 0x0.73p0 },
{ 1.369002422974590738036, 0x0.74p0 }, { 1.372714165087668414245, 0x0.75p0 }, { 1.376435970754530169202, 0x0.76p0 }, { 1.380167867260237990479, 0x0.77p0 },
{ 1.383909881963832022578, 0x0.78p0 }, { 1.387662042298529074813, 0x0.79p0 }, { 1.391424375771926236212, 0x0.7ap0 }, { 1.395196909966200271569, 0x0.7bp0 },
{ 1.398979672538311236352, 0x0.7cp0 }, { 1.402772691220204759333, 0x0.7dp0 }, { 1.406575993819015435449, 0x0.7ep0 }, { 1.410389608217270662749, 0x0.7fp0 },
{ 1.414213562373094923430, 0x0.80p0 }, { 1.418047884320415175097, 0x0.81p0 }, { 1.421892602169165575887, 0x0.82p0 }, { 1.425747744105494430045, 0x0.83p0 },
{ 1.429613338391970023267, 0x0.84p0 }, { 1.433489413367788900544, 0x0.85p0 }, { 1.437375997448982367644, 0x0.86p0 }, { 1.441273119128625657126, 0x0.87p0 },
{ 1.445180806977046650275, 0x0.88p0 }, { 1.449099089642035043113, 0x0.89p0 }, { 1.453027995849052622646, 0x0.8ap0 }, { 1.456967554401443765144, 0x0.8bp0 },
{ 1.460917794180647044655, 0x0.8cp0 }, { 1.464878744146405731286, 0x0.8dp0 }, { 1.468850433336981842203, 0x0.8ep0 }, { 1.472832890869367528097, 0x0.8fp0 },
{ 1.476826145939499346227, 0x0.90p0 }, { 1.480830227822471867327, 0x0.91p0 }, { 1.484845165872752614789, 0x0.92p0 }, { 1.488870989524397003834, 0x0.93p0 },
{ 1.492907728291264835008, 0x0.94p0 }, { 1.496955411767235455400, 0x0.95p0 }, { 1.501014069626425584403, 0x0.96p0 }, { 1.505083731623406473332, 0x0.97p0 },
{ 1.509164427593422841412, 0x0.98p0 }, { 1.513256187452609813349, 0x0.99p0 }, { 1.517359041198214741897, 0x0.9ap0 }, { 1.521473018908814589523, 0x0.9bp0 },
{ 1.525598150744538195056, 0x0.9cp0 }, { 1.529734466947286986027, 0x0.9dp0 }, { 1.533881997840955913048, 0x0.9ep0 }, { 1.538040773831656826687, 0x0.9fp0 },
{ 1.542210825407940744114, 0x0.a0p0 }, { 1.546392183141021670068, 0x0.a1p0 }, { 1.550584877684999973724, 0x0.a2p0 }, { 1.554788939777088430105, 0x0.a3p0 },
{ 1.559004400237836929222, 0x0.a4p0 }, { 1.563231289971357629298, 0x0.a5p0 }, { 1.567469639965552774541, 0x0.a6p0 }, { 1.571719481292341402678, 0x0.a7p0 },
{ 1.575980845107886496592, 0x0.a8p0 }, { 1.580253762652824578439, 0x0.a9p0 }, { 1.584538265252493749458, 0x0.aap0 }, { 1.588834384317163950229, 0x0.abp0 },
{ 1.593142151342266776837, 0x0.acp0 }, { 1.597461597908627073394, 0x0.adp0 }, { 1.601792755682693414343, 0x0.aep0 }, { 1.606135656416771029242, 0x0.afp0 },
{ 1.610490331949254283472, 0x0.b0p0 }, { 1.614856814204860491202, 0x0.b1p0 }, { 1.619235135194863728358, 0x0.b2p0 }, { 1.623625327017328867640, 0x0.b3p0 },
{ 1.628027421857347833978, 0x0.b4p0 }, { 1.632441451987274971813, 0x0.b5p0 }, { 1.636867449766964410784, 0x0.b6p0 }, { 1.641305447644006321184, 0x0.b7p0 },
{ 1.645755478153964945776, 0x0.b8p0 }, { 1.650217573920617741834, 0x0.b9p0 }, { 1.654691767656194523184, 0x0.bap0 }, { 1.659178092161616158151, 0x0.bbp0 },
{ 1.663676580326736598181, 0x0.bcp0 }, { 1.668187265130582463968, 0x0.bdp0 }, { 1.672710179641596406341, 0x0.bep0 }, { 1.677245357017878468753, 0x0.bfp0 },
{ 1.681792830507429004072, 0x0.c0p0 }, { 1.686352633448393145699, 0x0.c1p0 }, { 1.690924799269305056626, 0x0.c2p0 }, { 1.695509361489332622597, 0x0.c3p0 },
{ 1.700106353718523477525, 0x0.c4p0 }, { 1.704715809658051250963, 0x0.c5p0 }, { 1.709337763100462925792, 0x0.c6p0 }, { 1.713972247929925973864, 0x0.c7p0 },
{ 1.718619298122477934143, 0x0.c8p0 }, { 1.723278947746273992436, 0x0.c9p0 }, { 1.727951230961837669753, 0x0.cap0 }, { 1.732636182022311066575, 0x0.cbp0 },
{ 1.737333835273706217350, 0x0.ccp0 }, { 1.742044225155156444984, 0x0.cdp0 }, { 1.746767386199168825556, 0x0.cep0 }, { 1.751503353031877985302, 0x0.cfp0 },
{ 1.756252160373299453511, 0x0.d0p0 }, { 1.761013843037583681550, 0x0.d1p0 }, { 1.765788435933272726430, 0x0.d2p0 }, { 1.770575974063554713922, 0x0.d3p0 },
{ 1.775376492526521188253, 0x0.d4p0 }, { 1.780190026515424461806, 0x0.d5p0 }, { 1.785016611318934964814, 0x0.d6p0 }, { 1.789856282321401037549, 0x0.d7p0 },
{ 1.794709075003107168200, 0x0.d8p0 }, { 1.799575024940535117324, 0x0.d9p0 }, { 1.804454167806623932080, 0x0.dap0 }, { 1.809346539371031736820, 0x0.dbp0 },
{ 1.814252175500398633901, 0x0.dcp0 }, { 1.819171112158608494269, 0x0.ddp0 }, { 1.824103385407053190548, 0x0.dep0 }, { 1.829049031404897274200, 0x0.dfp0 },
{ 1.834008086409342430656, 0x0.e0p0 }, { 1.838980586775893710794, 0x0.e1p0 }, { 1.843966568958625984465, 0x0.e2p0 }, { 1.848966069510450838109, 0x0.e3p0 },
{ 1.853979125083385470774, 0x0.e4p0 }, { 1.859005772428820479902, 0x0.e5p0 }, { 1.864046048397788979401, 0x0.e6p0 }, { 1.869099989941238604274, 0x0.e7p0 },
{ 1.874167634110299962558, 0x0.e8p0 }, { 1.879249018056560194267, 0x0.e9p0 }, { 1.884344179032334309909, 0x0.eap0 }, { 1.889453154390938971474, 0x0.ebp0 },
{ 1.894575981586965607306, 0x0.ecp0 }, { 1.899712698176555081275, 0x0.edp0 }, { 1.904863341817674138312, 0x0.eep0 }, { 1.910027950270389629495, 0x0.efp0 },
{ 1.915206561397147178027, 0x0.f0p0 }, { 1.920399213163047402730, 0x0.f1p0 }, { 1.925605943636124806062, 0x0.f2p0 }, { 1.930826790987627106233, 0x0.f3p0 },
{ 1.936061793492294347274, 0x0.f4p0 }, { 1.941310989528640451596, 0x0.f5p0 }, { 1.946574417579233218234, 0x0.f6p0 }, { 1.951852116230978317901, 0x0.f7p0 },
{ 1.957144124175400401455, 0x0.f8p0 }, { 1.962450480208927316994, 0x0.f9p0 }, { 1.967771223233175659217, 0x0.fap0 }, { 1.973106392255234098343, 0x0.fbp0 },
{ 1.978456026387950927869, 0x0.fcp0 }, { 1.983820164850219391894, 0x0.fdp0 }, { 1.989198846967266343100, 0x0.fep0 }, { 1.994592112170940234606, 0x0.ffp0 }
};
static const double exp2_polynomial[3] = { 1.00000000011116223, 0.693146680445468386, 0.240559810709772152 };
static const xFloat huge = { 0x1.0p127f, 0, 0, 0 };
static const xFloat tiny = { 0x1.00001p-127f, 0, 0, 0 };
static const vUInt32 bias = { 1023, 0, 0, 0 };
static const __m128d almost1 = { 1.0 - DBL_EPSILON, 0.0 };
static const __m128d minusZero_sd = { -0.0, 0.0 };
float powf( float x, float y )
{
if( x == 1.0f || y == 0.0f )
return 1.0f;
if( x != x || y != y )
return x + y;
float fabsx = __builtin_fabsf( x );
float fabsy = __builtin_fabsf( y );
uint32_t yFracBits = 0;
uint32_t yOneBit = 0;
int32_t xExp, yExp;
union
{
float f;
uint32_t u;
}uy, ux;
uy.f = fabsy;
ux.f = fabsx;
yExp = (uy.u >> 23) - 127;
xExp = (ux.u >> 23) - 127;
if( fabsy < 1.0f )
yFracBits = 0x007FFFFF;
else if( fabsy < 0x1.0p24f ) {
yOneBit = (0x00800000 >> yExp) & uy.u;
yFracBits = (0x007FFFFF >> yExp) & uy.u;
}
if( x == 0.0f )
{
if( y < 0.0f )
{
x = 1.0f / x;
if( 0 != yFracBits || 0 == yOneBit ) x = __builtin_fabsf( x );
return x;
}
else
{
if( 0 != yFracBits || 0 == yOneBit ) return 0.0f;
return x;
}
}
if( fabsx == __builtin_inff() )
{
if( x == __builtin_inff() )
{
if( y < 0.0f )
return 0.0f;
return x;
}
else
{ if( y < 0.0f )
{
if( 0 != yFracBits || 0 == yOneBit ) return 0.0f;
return -0.0f;
}
else
{ if( 0 != yFracBits || 0 == yOneBit ) return __builtin_inff();
return x;
}
}
}
if( fabsy == __builtin_inff() )
{
if( x == -1.0f)
return 1.0f;
if( y == __builtin_inff() )
{
if( fabsx < 1.0f )
return 0.0f;
return __builtin_inff();
}
if( fabsx < 1.0f )
return __builtin_inff();
return 0.0f;
}
if( x < 0.0f && 0 != yFracBits )
{
SET_INVALID_FLAG();
return __builtin_nanf("37");
}
if( 0 == yFracBits && fabsy <= 0x1.0p16f )
{
int32_t iy = y; int32_t yIsNeg = iy >> 31;
iy = abs( iy );
double dx = x;
double result = iy & 1 ? dx : 1.0;
while( iy >>= 1 )
{
dx *= dx;
if( iy & 1 )
result *= dx;
}
if( yIsNeg )
return (float) (1.0 / result);
return (float) result;
}
if( 0.5f == fabsy )
{
double d = sqrt( (double) x );
if( y < 0.0f )
return (float) ( 1.0 / d );
return (float) d;
}
__m128d ylog2X;
if( fabsy > 52.0f || ( 0.98f < fabsx && fabsx < 1.02f ) )
{
double temp;
asm volatile( "flds (%1)\n flds (%2)\n fyl2x \n fstpl %0" : "=m" (*(&temp)) : "r" (&y), "r" (&fabsx) );
ylog2X = _mm_load_sd( &temp );
}
else
{
__m128d X = FLOAT_2_XDOUBLE( fabsx );
__m128d Y = FLOAT_2_XDOUBLE( y );
__m128d log2X = _mm_cvtepi32_pd( _mm_sub_epi32( _mm_srli_epi64( (__m128i) _mm_and_pd( X, exponentMask), 52 ), v1023 ) ); __m128d mantissa = _mm_or_pd( _mm_andnot_pd( exponentMask, X ), one_sd );
if( 1.0 != XDOUBLE_2_DOUBLE( mantissa ) )
{ __m128d vMantissa = _mm_unpacklo_pd( mantissa, mantissa ); __m128d m1 = _mm_unpacklo_pd( one_sd, mantissa ); __m128d m2 = _mm_mul_pd( vMantissa, vMantissa ); __m128d m4 = _mm_mul_pd( m2, m2 ); __m128d v1 = _mm_mul_pd( m4, log2polynomialCoefficients[5] ); __m128d v = _mm_mul_pd( m4, log2polynomialCoefficients[4] ); m2 = _mm_mul_pd( m2, m1 );
v1 = _mm_add_pd( v1, log2polynomialCoefficients[3] ); v = _mm_add_pd( v , log2polynomialCoefficients[2] ); v1 = _mm_mul_pd( v1, m4 ); v = _mm_mul_pd( v , m4 ); v1 = _mm_add_pd( v1, log2polynomialCoefficients[1] ); v = _mm_add_pd( v , log2polynomialCoefficients[0] ); v1 = _mm_mul_pd( v1, m2 ); v = _mm_mul_pd( v , m1 ); v = _mm_add_pd( v , v1 );
#if defined( __SSE3__ )
v = _mm_hadd_pd( v, _mm_setzero_pd() );
#else
v = _mm_add_sd( v, _mm_unpackhi_pd( v, v ) );
#endif
log2X = _mm_add_sd( log2X, v );
}
ylog2X = _mm_mul_sd( log2X, Y );
}
double q = XDOUBLE_2_DOUBLE( ylog2X );
if( q >= 128.0 )
return XFLOAT_2_FLOAT( REQUIRED_ADD_ss( huge, huge ) );
if( q <= -150.0 )
return XFLOAT_2_FLOAT( REQUIRED_MULTIPLY_ss( tiny, tiny ) );
static const __m128d magic = { 0x1.0p52, 0.0 };
__m128d smagic = _mm_or_pd( magic, _mm_and_pd( ylog2X, minusZero_sd ) );
__m128d vi = _mm_sub_sd( _mm_add_sd( ylog2X, smagic ), smagic ); vi = _mm_sub_sd( vi, _mm_and_pd( _mm_cmpgt_sd( vi, ylog2X ), one_sd ) );
if( _mm_ucomieq_sd( vi, ylog2X ) )
return XDOUBLE_2_FLOAT( (xDouble) _mm_slli_epi64( _mm_add_epi32( _mm_cvtpd_epi32( vi ), bias ), 52 ) );
__m128d vf = _mm_sub_sd( ylog2X, vi );
vf = _mm_min_sd( vf, almost1 );
int n = _mm_cvtsi128_si32( _mm_srli_epi64( _mm_and_si128( (__m128i) _mm_add_sd( vf, one_sd), mantissaHiBits ), 44 ) );
vf = _mm_sub_sdm( vf, (double*)&exp2_table[ n ] + 1 );
xDouble e = _mm_mul_sdm( vf, exp2_polynomial + 2 );
e = _mm_add_sdm( e, exp2_polynomial + 1 );
e = _mm_mul_sd( e, vf );
e = _mm_add_sdm( e, exp2_polynomial );
__m128i ei = _mm_slli_epi64( _mm_cvtpd_epi32( vi ), 52 );
ei = _mm_add_epi64( ei, (__m128i) exp2_table[ n ] );
e = _mm_mul_sd( e, (__m128d) ei );
float result = XDOUBLE_2_FLOAT( e );
if( x < 0.0f && yOneBit )
return -result;
return result;
}