expf_logf_powf.c   [plain text]

/*
 *  expf_logf_powf.c
 *  Libm
 *
 *  Written by Ian Ollmann
 *  Polynomials by Ali Sazegari.
 *
 *  Copyright 2007 Apple Inc. All rights reserved.
 *
 */

#include "math.h"
#include <emmintrin.h>
#include <stdint.h>
#include "xmmLibm_Prefix.h"
#include <float.h>

#if defined( __SSE3__ )
    #include <pmmintrin.h>
#endif

                                                         
//11th order minimax polynomial approximation of y = log2(x)
//Valid over the range [1, 2]. Error is 0.160490514e-9.
//Packaged as vector doubles, in order from c0 to c11
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 }
                                                    };


//Second order minimax polynomial approximation of  y = 2**x
//Valid over the range [0.0,0.4e-2]. (We need [0, 2**-8] )
//  error is +- 0.11116223e-9,  We need less than 2**-28 here.
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 };

/*
float log2f( float x )
{
//Edge cases for C99 TC2:

    //Handle NaN
    if( x != x )
        return x;

// log2(+-0)returns -Inf and raises the divide-by-zero floating-point exception. 
// log2(x)returns a NaN and raises the  invalid floating-point exception for x < 0. 
    if( x <= 0.0f )
    {
        if( x == 0.0f )
            return XFLOAT_2_FLOAT( REQUIRED_DIVIDE_ss( FLOAT_2_XFLOAT(-1.0f), _mm_setzero_ps() ));
        else
            return XFLOAT_2_FLOAT( REQUIRED_MULTIPLY_ss( FLOAT_2_XFLOAT(__builtin_inff()), _mm_setzero_ps() ));
    }
    
// log2(+Inf)returns +Inf. 
    if( x == __builtin_inff() )
        return x;

// We are going to try to get the last edge case to just happen in arithmetic
// log2(1)returns +0. 


//Main body

    //Convert to double  
    __m128d X = FLOAT_2_XDOUBLE(x);      //{ 0, x }

    //extract exponent -- denormals can't happen here
    __m128i iexponent = _mm_sub_epi32( _mm_srli_epi64( (__m128i) _mm_and_pd( X, exponentMask), 52 ), v1023 );
    __m128d exponent = _mm_cvtepi32_pd( iexponent );

    //get mantissa of X as value in the range [ 1.0, 2.0 )
    __m128d mantissa = _mm_or_pd( _mm_andnot_pd( exponentMask, X ), one_sd );

    //avoid setting inexact flag for exact powers of two
    if( 1.0 == XDOUBLE_2_DOUBLE( mantissa ) )
        return XDOUBLE_2_FLOAT( exponent );

    __m128d m1, m2, m4, v1, v;

#if 1
    if( 0.98f < x && x < 1.02f )
    {
        float temp;

        //use fyl2x to get y*log2(x)
        asm volatile( "fld1 \n flds (%1)\n fyl2x \n fstps %0" : "=m" (*(&temp)) : "r" (&x) );
            
        return temp;
    }
#else
    //inputs between 0.98 and 1.02 need special care
    if( 0.98f < x && x < 1.02f )
    {
        __m128d XX = _mm_unpacklo_pd( X, X );
        m1 = _mm_unpacklo_pd( one_sd, X );                          // { m, 1 }
        m2 = _mm_mul_pd( XX, XX );                                   //{ m**2, m**2 }
        m4 = _mm_mul_pd( m2, m2 );                                  //{ m**4, m**4 }

        v1 = _mm_mul_pd( m4, log2near1polynomialCoefficients[5] );               // { c11 m**4, c10 m**4 } 
        v  = _mm_mul_pd( m4, log2near1polynomialCoefficients[4] );               // { c9 m**4, c8 m**4 } 
        m2 = _mm_mul_pd( m2, m1 );
        v1 = _mm_add_pd( v1, log2near1polynomialCoefficients[3] );                       // { c7 + c11 m**4, c6 + c10 m**4 }
        v  = _mm_add_pd( v , log2near1polynomialCoefficients[2] );                       // { c5 + c9 m**4, c4 + c8 m**4 }
        v1 = _mm_mul_pd( v1, m4 );                                                  // { c7 m**4 + c11 m**8, c6 m**4 + c10 m**8 }
        v  = _mm_mul_pd( v , m4 );                                                  // { c5 m**4 + c9 m**8, c4 m**4 + c8 m**8 }
        v1 = _mm_add_pd( v1, log2near1polynomialCoefficients[1] );                       // { c3 + c7 m**4 + c11 m**8, c2 + c6 m**4 + c10 m**8 }
        v  = _mm_add_pd( v , log2near1polynomialCoefficients[0] );                       // { c1 + c5 m**4 + c9 m**8, c0 + c4 m**4 + c8 m**8 }
        v1 = _mm_mul_pd( v1, m2 );                                                  // { c3 m**3 + c7 m**7 + c11 m**11, c2 m**2 + c6 m**4 + c10 m**8 }
        v  = _mm_mul_pd( v , m1 );                                                  // { c1 m + c5 m**5 + c9 m**9, c0 + c4 m**4 + c8 m**8 }
        v  = _mm_add_pd( v , v1 );

        //Add the two sub terms horizontally
    #if defined( __SSE3__ )
        v = _mm_hadd_pd(  v, _mm_setzero_pd() );
    #else
        v = _mm_add_sd( v, _mm_unpackhi_pd( v, v ) );
    #endif
    
        return XDOUBLE_2_FLOAT( v );
    }
#endif

    //First make a vector out of the mantissa
    __m128d vMantissa = _mm_unpacklo_pd( mantissa, mantissa );                  //{ mantissa, mantissa }
    m1 = _mm_unpacklo_pd( one_sd, mantissa );                           // { m, 1 }
    m2 = _mm_mul_pd( vMantissa, vMantissa );                            //{ m**2, m**2 }
    m4 = _mm_mul_pd( m2, m2 );                                          //{ m**4, m**4 }
    v1 = _mm_mul_pd( m4, log2polynomialCoefficients[5] );               // { c11 m**4, c10 m**4 } 
    v  = _mm_mul_pd( m4, log2polynomialCoefficients[4] );               // { c9 m**4, c8 m**4 } 
    m2 = _mm_mul_pd( m2, m1 );
    v1 = _mm_add_pd( v1, log2polynomialCoefficients[3] );                       // { c7 + c11 m**4, c6 + c10 m**4 }
    v  = _mm_add_pd( v , log2polynomialCoefficients[2] );                       // { c5 + c9 m**4, c4 + c8 m**4 }
    v1 = _mm_mul_pd( v1, m4 );                                                  // { c7 m**4 + c11 m**8, c6 m**4 + c10 m**8 }
    v  = _mm_mul_pd( v , m4 );                                                  // { c5 m**4 + c9 m**8, c4 m**4 + c8 m**8 }
    v1 = _mm_add_pd( v1, log2polynomialCoefficients[1] );                       // { c3 + c7 m**4 + c11 m**8, c2 + c6 m**4 + c10 m**8 }
    v  = _mm_add_pd( v , log2polynomialCoefficients[0] );                       // { c1 + c5 m**4 + c9 m**8, c0 + c4 m**4 + c8 m**8 }
    v1 = _mm_mul_pd( v1, m2 );                                                  // { c3 m**3 + c7 m**7 + c11 m**11, c2 m**2 + c6 m**4 + c10 m**8 }
    v  = _mm_mul_pd( v , m1 );                                                  // { c1 m + c5 m**5 + c9 m**9, c0 + c4 m**4 + c8 m**8 }
    v  = _mm_add_pd( v , v1 );

    //Add the two sub terms horizontally
#if defined( __SSE3__ )
    v = _mm_hadd_pd(  v, _mm_setzero_pd() );
#else
    v = _mm_add_sd( v, _mm_unpackhi_pd( v, v ) );
#endif


    //add in the exponent
    v = _mm_add_sd( v, exponent );
    
    return XDOUBLE_2_FLOAT( v );
}
*/

static const __m128i mantissaHiBits = { 0x000FF00000000000ULL, 0ULL };

/* 
    exp2_table produced as follows:
    
    #include <emmintrin.h>
    #include <stdio.h>
    #include <stdint.h>
    #include <math.h>


    int main( void )
    {
            int i;

            for( i = 0; i < 256; i++ )
            {
                    union
                    {
                            uint64_t        u;
                            double          d;
                    }u = { 0x3ff0000000000000ULL | ((uint64_t) i << 44 ) };

                    if( 0 == i % 4 )
                            printf( "\n" );

                    printf( " {  %17.21f, 0x0.%2.2xp0 },", exp2( u.d - 1.0 ), i  );

            }
            printf( "\n" );

            return 0;
    }
*/

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 }
                                             };     
                                             
// Minimax polynomial approximation for y = 2**x, valid over the range [ 0, 2**-8 ]. Accurate to within 0.11116223e-9. We need to within 2**(-28) here.                                          
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 };

/*
float exp2f( float x )
{
    //take care of NaN
    if( x != x )
        return x;

    //deal with overflow
    if( x >= 128.0f )
    {
        if( __builtin_inff() == x )
            return x;

        return XFLOAT_2_FLOAT( REQUIRED_ADD_ss( huge, huge ) );
    }

    //deal with underflow
    if( x <= -150.0f )
    {
        if( -__builtin_inff() == x )
            return 0;

        return XFLOAT_2_FLOAT( REQUIRED_MULTIPLY_ss( tiny, tiny ) );
    }



    // exp2(+-0)returns 1. 
    // Which will be taken care of through normal arithmetic

    //Split x into fractional and integer parts
    double xd = x;
    // i = floorf( x )
    double magic = x < 0 ? -0x1.0p52f : 0x1.0p52f;
    double i = (xd + magic) - magic;
    if( i == xd )
        return  XDOUBLE_2_FLOAT( (xDouble) _mm_slli_epi64( _mm_add_epi32( _mm_cvtpd_epi32( DOUBLE_2_XDOUBLE( i ) ), bias ), 52 ) );
    else if( i > xd )
        i -= 1.0f;
    double f = xd - i;

    //We do exp2(x) as:
    //
    //      x = i + f               i = floor( x )     f = x - i
    //
    //      exp2(x) = exp2( i ) * exp2( f )
    //
    //  We further break down f into:
    //
    //      f = n * (2**-8) + ff
    //
    //      exp2(x) = exp2 (i ) * exp2( n * (2**-8) ) * exp2(ff)
    //
    //          exp2( i )           calculated by simple slamming of i + 1023 into exponent of a double (with suitable saturation)
    //          exp2( n*(2**-8) )   calculated by table lookup into exp2_table
    //          exp2( ff )          calculated by minimax polynomial

    //reduce the fraction to the range [ 0, 2**-8 )
    xDouble d = DOUBLE_2_XDOUBLE(f);   
    
    //For very small negative x, we can now have d = 1.0, which causes trouble. 
    //In order to keep things sane we clamp vf to be 1.0 - DBL_EPSILON
    d = _mm_min_sd( d, almost1 );

    //extract the high 8 bits of the mantissa by adding one to d, and using integer operations to read out those bits
    int n = _mm_cvtsi128_si32( _mm_srli_epi64( _mm_and_si128( (__m128i) _mm_add_sd( d, one_sd), mantissaHiBits ), 44 ) );

    //subtract out a representation of the high 8 bits we removed from d
    d = _mm_sub_sdm( d, (double*)&exp2_table[ n ] + 1 );     //reduces d to [ 0, 2**-8 )

    
    //d = exp2(d), accurate to within 0.11116223e-9. We need to within 2**(-28) here 
    xDouble e = _mm_mul_sdm( d, exp2_polynomial + 2 );
    e = _mm_add_sdm( e, exp2_polynomial + 1 );
    e = _mm_mul_sd( e, d );
    e = _mm_add_sdm( e, exp2_polynomial );

    //Calculate exp2(i). Denormals are not possible because we use double precision here.
    __m128i ei = _mm_slli_epi64( _mm_cvtpd_epi32( DOUBLE_2_XDOUBLE( i ) ), 52 );     
    
    ei = _mm_add_epi64( ei, (__m128i) exp2_table[ n ] );    

    //scale by ei
    e = _mm_mul_sd( e, (__m128d) ei );
    
    return XDOUBLE_2_FLOAT( e );
}
*/

static const __m128d minusZero_sd = { -0.0, 0.0 };

float powf( float x, float y )
{
// The pow functions compute x raised to the power y. A domain error occurs if x is finite 
// and negative and y is finite and not an integer value. A range error may occur. A domain 
// error may occur if x is zero and y is zero. A domain error or range error may occur if x 
// is zero and y is less than zero. 

    // pow(+1, y) returns 1 for any y, even a NaN. 
    // pow(x, +-0) returns 1 for any x, even a NaN. 
    if( x == 1.0f || y == 0.0f )
        return 1.0f;

    // deal with NaNs before they cause trouble
    if( x != x || y != y )
        return x + y;

    float fabsx = __builtin_fabsf( x );
    float fabsy = __builtin_fabsf( y );
    
    //Find out if Y is an integer and if so if it is even or odd, without setting the inexact flag
    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 )   // y can only be an odd integer if 1.0f <= |y| < 2**24
    {
        yOneBit = (0x00800000 >> yExp) & uy.u;
        yFracBits = (0x007FFFFF >> yExp) & uy.u;
    }

    if( x == 0.0f  )
    {
        if( y < 0.0f )
        {
            // pow(+-0, y) returns +-Inf and raises the ‘‘divide-by-zero’’ floating-point exception for y an odd integer < 0. 
            // pow(+-0, y) returns +Inf and raises the ‘‘divide-by-zero’’ floating-point exception for y < 0 and not an odd integer. 
            x = 1.0f / x;   //set div/0 flag and x to infinity of right sign

            if( 0 != yFracBits || 0 == yOneBit )    //if y is not an odd integer
                x = __builtin_fabsf( x );

            return x;
        }
        else  
        {   //y > 0. y == 0 case taken care of at the top of the function.

            // pow(+-0, y) returns +-0 for y an odd integer > 0. 
            // pow(+-0, y) returns +0 for y > 0 and not an odd integer. 
            if( 0 != yFracBits || 0 == yOneBit )    //if y is not an odd integer
                return 0.0f;
        
            return x;
        }
    }

    if( fabsx == __builtin_inff() )
    {
        // pow(+Inf, y) returns +0 for y < 0. 
        // pow(+Inf, y) returns +Inf for y > 0. 
        if( x == __builtin_inff() )
        {
            if( y < 0.0f )
                return 0.0f;
        
            return x;
        }
        else
        { //x = -Inf
            if( y < 0.0f )
            {
                // pow(-Inf, y) returns −0 for y an odd integer < 0. 
                // pow(-Inf, y) returns +0 for y < 0 and not an odd integer. 
                if( 0 != yFracBits || 0 == yOneBit )    //if y is not an odd integer
                    return 0.0f;
                    
                return -0.0f;
            }
            else
            { //y > 0. y == 0 case taken care of at the top of the function.
                // pow(−Inf, y) returns -Inf for y an odd integer > 0. 
                // pow(−Inf, y) returns +Inf for y > 0 and not an odd integer. 
                if( 0 != yFracBits || 0 == yOneBit )    //if y is not an odd integer
                    return __builtin_inff();
            
                return x;
            }
        }
    }
    
    if( fabsy == __builtin_inff() )
    {
        // pow(−1, +-Inf) returns 1. 
        if( x == -1.0f)
            return 1.0f;

        if( y == __builtin_inff() )
        {
            // pow(x, +Inf) returns +0 for |x| < 1. 
            if( fabsx < 1.0f )
                return 0.0f;
                
            // pow(x, +inf) returns +Inf for |x| > 1. 
            return __builtin_inff();
        }

        // pow(x, -Inf) returns +Inf for |x| < 1. 
        if( fabsx < 1.0f )
            return __builtin_inff();
        // pow(x, -Inf) returns +0 for |x| > 1. 
        return 0.0f;
    }

    // pow(x, y) returns a NaN and raises the ‘‘invalid’’ floating-point exception for finite x < 0 and finite non-integer y. 
    if( x < 0.0f && 0 != yFracBits )
    {
        SET_INVALID_FLAG();
        return __builtin_nanf("37");
    }

    //if y is an integer, less than 2**16, do iPow
    if( 0 == yFracBits && fabsy <= 0x1.0p16f )
    {
        int32_t iy = y; //should be exact
        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;
        }
    
        //We are using double precision here, so we don't need to worry about range differences between tiny vs huge numbers for negative Y
        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;
    }
    
    //deal with 
        
    //Do pow the hard way

    //Convert x and y to double to do the arithmetic
    __m128d ylog2X;

    //We have two fallover cases for the log2 estimate above. 
    //
    //The first happens for values very near 1. This causes results near 0,
    //which means that our minimax error estimate is way off, because the results
    //are in a much smaller binade than normal, so the minimax error seems that much 
    //larger.
    //
    //The second case is for large y. This is a simple matter of minimax error * y
    //eventually exceeds our error threshold.
    //
    //In each case, we just fall back on the 80-bit log2 instruction, which costs us
    //an extra 40 cycles or so.
    //
    if( fabsy > 52.0f || ( 0.98f < fabsx && fabsx < 1.02f ) )
    {
        double temp;

        //use fyl2x to get y*log2(x)
        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 );

        //Extract the exponent and mantissa. We are essentially doing frexp here, except that the mantissa is [1, 2)
        __m128d log2X = _mm_cvtepi32_pd( _mm_sub_epi32( _mm_srli_epi64( (__m128i) _mm_and_pd( X, exponentMask), 52 ), v1023 ) );     //always exact
        __m128d mantissa = _mm_or_pd( _mm_andnot_pd( exponentMask, X ), one_sd );

        //avoid setting inexact flag for exact powers of two
        if( 1.0 != XDOUBLE_2_DOUBLE( mantissa ) )
        { //take log2( mantissa )
            __m128d vMantissa = _mm_unpacklo_pd( mantissa, mantissa );                  //{ mantissa, mantissa }
            __m128d m1 = _mm_unpacklo_pd( one_sd, mantissa );                           // { m, 1 }
            __m128d m2 = _mm_mul_pd( vMantissa, vMantissa );                            //{ m**2, m**2 }
            __m128d m4 = _mm_mul_pd( m2, m2 );                                          //{ m**4, m**4 }
            __m128d v1 = _mm_mul_pd( m4, log2polynomialCoefficients[5] );               // { c11 m**4, c10 m**4 } 
            __m128d v  = _mm_mul_pd( m4, log2polynomialCoefficients[4] );               // { c9 m**4, c8 m**4 } 
            m2 = _mm_mul_pd( m2, m1 );
            v1 = _mm_add_pd( v1, log2polynomialCoefficients[3] );                       // { c7 + c11 m**4, c6 + c10 m**4 }
            v  = _mm_add_pd( v , log2polynomialCoefficients[2] );                       // { c5 + c9 m**4, c4 + c8 m**4 }
            v1 = _mm_mul_pd( v1, m4 );                                                  // { c7 m**4 + c11 m**8, c6 m**4 + c10 m**8 }
            v  = _mm_mul_pd( v , m4 );                                                  // { c5 m**4 + c9 m**8, c4 m**4 + c8 m**8 }
            v1 = _mm_add_pd( v1, log2polynomialCoefficients[1] );                       // { c3 + c7 m**4 + c11 m**8, c2 + c6 m**4 + c10 m**8 }
            v  = _mm_add_pd( v , log2polynomialCoefficients[0] );                       // { c1 + c5 m**4 + c9 m**8, c0 + c4 m**4 + c8 m**8 }
            v1 = _mm_mul_pd( v1, m2 );                                                  // { c3 m**3 + c7 m**7 + c11 m**11, c2 m**2 + c6 m**4 + c10 m**8 }
            v  = _mm_mul_pd( v , m1 );                                                  // { c1 m + c5 m**5 + c9 m**9, c0 + c4 m**4 + c8 m**8 }
            v  = _mm_add_pd( v , v1 );

            //Add the two sub terms horizontally
        #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 );
    }

    // ylog2X = y * log2( x )
    double q = XDOUBLE_2_DOUBLE( ylog2X );
    
    //deal with overflow
    if( q >= 128.0 )
        return XFLOAT_2_FLOAT( REQUIRED_ADD_ss( huge, huge ) );

    //deal with underflow to zero
    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 ) );    //smagic = copysign( 0x1.0p52, ylog2X )

    // separate log2X into fractional and integer components
    __m128d vi = _mm_sub_sd( _mm_add_sd( ylog2X, smagic ), smagic );            // vi = rint( ylog2X )
    vi = _mm_sub_sd( vi, _mm_and_pd( _mm_cmpgt_sd( vi, ylog2X ), one_sd ) );     // vi = floor( ylog2X )
    
    //Avoid setting inexact if ylog2X is an integer
    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 = fractional part ( 0 <= vf < 1.0 )
    
    //For very small negative ylog2X, we can now have vf = 1.0, which causes trouble. 
    //In order to keep things sane we clamp vf to be 1.0 - EPSILON
    vf = _mm_min_sd( vf, almost1 );

    //reduce the fraction to the range [ 0, 2**-8 )
    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 );     //reduces d to [ 0, 2**-8 )
    
    //d = exp2(d), accurate to within 0.11116223e-9. We need to within 2**(-28) here 
    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 );

    //Calculate exp2(i). Denormals are not possible because we use double precision here.
    __m128i ei = _mm_slli_epi64( _mm_cvtpd_epi32( vi ), 52 );     
    ei = _mm_add_epi64( ei, (__m128i) exp2_table[ n ] );    

    //scale by ei
    e = _mm_mul_sd( e, (__m128d) ei );
    
    float result =  XDOUBLE_2_FLOAT( e );

    if( x < 0.0f && yOneBit )
        return -result;

    return result;
}