AAX_FastPow.h

Go to the documentation of this file.

/*================================================================================================*/
/*
 *  Copyright 2009-2015, 2023-2024 Avid Technology, Inc.
 *  All rights reserved.
 *  
 *  This file is part of the Avid AAX SDK.
 *  
 *  The AAX SDK is subject to commercial or open-source licensing.
 *  
 *  By using the AAX SDK, you agree to the terms of both the Avid AAX SDK License
 *  Agreement and Avid Privacy Policy.
 *  
 *  AAX SDK License: https://developer.avid.com/aax
 *  Privacy Policy: https://www.avid.com/legal/privacy-policy-statement
 *  
 *  Or: You may also use this code under the terms of the GPL v3 (see
 *  www.gnu.org/licenses).
 *  
 *  THE AAX SDK IS PROVIDED "AS IS" WITHOUT ANY WARRANTY, AND ALL WARRANTIES, WHETHER
 *  EXPRESSED OR IMPLIED, INCLUDING MERCHANTABILITY AND FITNESS FOR PURPOSE, ARE
 *  DISCLAIMED.
 */
 
/*================================================================================================*/
 
 
 
 
#pragma once
 
#ifndef _AAX_FASTPOW_H_
#define _AAX_FASTPOW_H_
 
#include <cmath> 
#if defined(MAC_VERSION) || defined(LINUX_VERSION)
namespace std  {using ::powf;using ::log10f;using ::fabsf;}
#endif
#include "AAX.h"    // For AAX_RESTRICT
 
namespace AAX
{
static const float _2p23 = 8388608.0f;
const unsigned int kPowExtent = 9;  // 9 results in ~-80dB cancellation
const unsigned int kPowTableSize = 1 << kPowExtent;
    
static inline void PowFastSetTable(
    float* const AAX_RESTRICT pTable,
    const unsigned int precision,
    const unsigned int extent,
    const bool isRound)
{
   // step along table elements and x-axis positions
   float zeroToOne = !isRound ?
      0.0f : (1.0f / (static_cast<float>(1 << precision) * 2.0f));
   for( int i = 0;  i < (1 << extent);  ++i )
   {
      // make y-axis value for table element
      pTable[i] = std::powf( 2.0f, zeroToOne );
 
      zeroToOne += 1.0f / static_cast<float>(1 << precision);
   }
}
    
static inline float PowFastLookup(
    const float         val,
    const float         ilog2,
    const float* const AAX_RESTRICT tableH,
    const float* const AAX_RESTRICT tableL)
{
 
    // build float bits
    const int i = static_cast<int>( (val * (_2p23 * ilog2)) + (127.0f * _2p23) );
 
    // replace mantissa with combined lookups
    const float t  = tableH[(i >> 14) & 0x1FF] * tableL[(i >> 5) & 0x1FF];
    const int   it = (i & 0xFF800000) |
      (*reinterpret_cast<const int*>( &t ) & 0x7FFFFF);
 
    // convert bits to float
    return *reinterpret_cast<const float*>( &it );
 
}
 
static inline float Pow2FastLookup(
    const float         val,
    const float* const AAX_RESTRICT tableH,
    const float* const AAX_RESTRICT tableL)
{
 
    // build float bits
    const int i = static_cast<int>((val + 127.0f) * _2p23);
 
    // replace mantissa with combined lookups
    const float t  = tableH[(i >> 14) & 0x1FF] * tableL[(i >> 5) & 0x1FF];
    const int   it = (i & 0xFF800000) |
      (*reinterpret_cast<const int*>( &t ) & 0x007FFFFF);
 
    // convert bits to float
    return *reinterpret_cast<const float*>( &it );
 
}
 
 
static const int kMantissaBitSize = 23;
static const int kExpBias = 127;
static const int kLogMantissaMSBs = 9;
static const int kLogTableSize = (1 << kLogMantissaMSBs) + 1;// 513;    
 
static inline void LogFastSetTable(
    float* const AAX_RESTRICT logTable, 
    int tableSize = kLogTableSize)
{
    const float kInv = 1.0f/float(tableSize - 1); // 0.001953125
    float mantissaVal = 1.0f;   //Start here.
    const float invLog10Base2 = 1.0f / std::log10f(2.0f);
    for (int j = 0; j < tableSize; j++)
    {
        logTable[j] = std::log10f(mantissaVal) * invLog10Base2; //By log equivalency: log2(x) = log10(x)/log10(2)
        mantissaVal += kInv;
        
    }
 
    return;
}
 
template <int kMantMSBs>
static inline void LogFloatToExpMantissa(
    float number, 
    int * const AAX_RESTRICT outExp, 
    int * const AAX_RESTRICT outMant, 
    float* const AAX_RESTRICT fract)
{   
    const int mantissaSize = 1 << (kMantissaBitSize - kMantMSBs);
    const float invMantissaSize = 1.0f/(float)mantissaSize;
    const int intEquiv = *reinterpret_cast<const int*>(&number);
 
    //bool isNegative = (result & 0x80000000) != 0; //Sign bit
    const int exponent      = (intEquiv & 0x7f800000) >> 23; //Exponent bits
    const int mantissa      =  intEquiv & 0x007fffff; //Mantissa bits.
 
    // 0 000|0 000|0 1 010| 0000 0000 0000 0000 0000
    // 0111 1111 1000 0000 0000 0000 0000 0000 = 7f800000 (Exp)                                            
    // 0000 0000 0111 1111 1111 1111 1111 1111 = 007fffff (Mant)
 
    *outExp = exponent - kExpBias;          //Compensate for exponent bias for caller.
    //To do: add halfLSBShifted to mantissa before shifting; *outMant = (mantissa + halfLSBShifted) >> (kMantissaBitSize - mantMSBs)
    *outMant = mantissa >> (kMantissaBitSize - kMantMSBs);      //Shift to provide only mantMSBs for log LUT.
    // Provide a linear interpolation fraction
    *fract = (float)(mantissa & (mantissaSize - 1)) * invMantissaSize;
 
    return;
}     
 
template <int kPrecision>
static inline float LogFastLookup(
    float num, 
    const float* const AAX_RESTRICT logTable)
{
 
    int exponent, mantissa;
    float fract;
 
    AAX::LogFloatToExpMantissa<kPrecision>(num, &exponent, &mantissa, &fract);
 
    const float mantLog1 = logTable[mantissa];
    const float mantLog2 = logTable[mantissa+1];
 
    const float logOfNum = (mantLog1 * (1.0f - fract) + fract * mantLog2) + exponent;
    return logOfNum;
}  
 
 
 
} // namespace AAX
 
#endif // #ifndef _AAX_FASTPOW_H_