/*
Harrison Ainsworth / HXA7241 : 2004-2011.

Redistribution and use in source and binary forms, with or without modification,
are permitted provided that the following conditions are met:

* Redistributions of source code must retain the above copyright notice, this
  list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright notice, this
  list of conditions and the following disclaimer in the documentation and/or
  other materials provided with the distribution.
* The name of the author may not be used to endorse or promote products derived
  from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE AUTHOR "AS IS" AND ANY EXPRESS OR IMPLIED
WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
OF SUCH DAMAGE.
*/

#include <stdint.h>
#include "math.h"
#include "FastMath.h"

#define FastPowPrecision    15
#if (FastPowPrecision==11)
int FastPowTable[2048];
#endif
#if (FastPowPrecision==12)
int FastPowTable[2048*2];
#endif
#if (FastPowPrecision==13)
int FastPowTable[2048*4];
#endif
#if (FastPowPrecision==14)
int FastPowTable[2048*8];
#endif
#if (FastPowPrecision==15)
int FastPowTable[2048*16];
#endif

__attribute__((__section__(".FastCode")))
float fastlog2 (float x)
{
  union { float f; uint32_t i; } vx = { x };
  union { uint32_t i; float f; } mx = { (vx.i & 0x007FFFFF) | (0x7e << 23) };
  float y = vx.i;
  y *= 1.0 / (1 << 23);

  return
    y - 124.22544637f - 1.498030302f * mx.f - 1.72587999f / (0.3520887068f + mx.f);
}

__attribute__((__section__(".FastCode")))
float fastlog (float x)
{
  return 0.69314718f * fastlog2 (x);
}

static const float _2p23 = 8388608.0f;

static void powFastSetTable
(
   unsigned int* const pTable,
   const unsigned int  precision
)
{
    int   i;
    /* step along table elements and x-axis positions */
    float zeroToOne = 1.0f / ((float)(1 << precision) * 2.0f);
    
    for( i = 0;  i < (1 << precision);  ++i )
    {
        /* make y-axis value for table element */
        //must be volatile, otherwise when optimizing is turned on, running will fail
        volatile float f = ((float)pow( 2.0f, zeroToOne ) - 1.0f) * _2p23;
        pTable[i] = (unsigned int)( f < _2p23 ? f : (_2p23 - 1.0f) );
        zeroToOne += 1.0f / (float)(1 << precision);
    }
}


__attribute__((__section__(".FastCode")))
static float powFastLookup
(
   const float         val,
   const float         ilog2,
   unsigned int* const pTable,
   const unsigned int  precision
)
{
   /* build float bits */
   const int i = (int)( (val * (_2p23 * ilog2)) + (127.0f * _2p23) );

   /* replace mantissa with lookup */
   const int it = (i & 0xFF800000) | pTable[(i & 0x7FFFFF) >> (23 - precision)];

   /* convert bits to float */
   union { int i; float f; } pun;
   return pun.i = it,  pun.f;
}


typedef struct {
   unsigned int  precision_m;
   unsigned int* pTable_m;
} PowFast;



PowFast PowFastInitData;
void powFastCreate(void)
{
    PowFast* pPowFast;
    unsigned int precision;
    float a,b,c;
    
    
    pPowFast =&PowFastInitData;
    precision = FastPowPrecision;

    pPowFast->precision_m = precision;
    pPowFast->pTable_m    = (unsigned int *)FastPowTable;

    powFastSetTable( pPowFast->pTable_m, pPowFast->precision_m );

    a=0.1f; b=0.2f;
    
#if (FastPowPrecision==11)
    c=powFast2(a);  if(c!=1.07166469f) while(1) {};
    c=powFastE(a);  if(c!=1.10518444f) while(1) {};
    c=powFast10(a); if(c!=1.25899744f) while(1) {};
    c=powFast(a,b); if(c!=0.63099187f) while(1) {};
#endif

#if (FastPowPrecision==12)
    c=powFast2(a);  if(c!=1.07175529f) while(1) {};
    c=powFastE(a);  if(c!=1.10509086f) while(1) {};
    c=powFast10(a); if(c!=1.25889099f) while(1) {};
    c=powFast(a,b); if(c!=0.63093847f) while(1) {};
#endif

}

#if 0
__attribute__((__section__(".FastCode")))
float powFast2
(
   float       f
)
{
   return powFastLookup( f, 1.0f, PowFastInitData.pTable_m, PowFastInitData.precision_m );
}

__attribute__((__section__(".FastCode")))
float powFastE
(
   float       f
)
{
   return powFastLookup( f, 1.44269504088896f, PowFastInitData.pTable_m, PowFastInitData.precision_m);
}

__attribute__((__section__(".FastCode")))
float powFast10
(
   float       f
)
{
   return powFastLookup( f, 3.32192809488736f, PowFastInitData.pTable_m, PowFastInitData.precision_m);
}
#endif

__attribute__((__section__(".FastCode")))
float powFast
(
   float       r,
   float       f
)
{
   return powFastLookup( f, (fastlog(r) * 1.44269504088896f), PowFastInitData.pTable_m, PowFastInitData.precision_m );
}



