CMathTools.hpp

Go to the documentation of this file.

#ifndef __CMathTools__
#define __CMathTools__

#include "../Basics/CCountedObject.hpp"
#include <algorithm>
using namespace std;

#ifdef WIN32
    #define IMAN 0
#else
    #ifdef __i386
        #define IMAN 0
    #else
        #define IMAN 1
    #endif
#endif

#ifdef WIN32
    #define exponent_position 1
#else
    #define exponent_position 0
#endif


static float shift23=(1<<23);
static float PowBodge=0.33971f;

//  ===========================================================================

namespace Exponent
{
    namespace MathTools
    {
        class CMathTools
        {
        public:

//  ===========================================================================

            const static float CMATH_HALF_PI_FLOAT;             
            const static float CMATH_PI_FLOAT;                  
            const static float CMATH_2PI_FLOAT;                 
            const static float CMATH_4PI_FLOAT;                 
//  ===========================================================================

            const static double CMATH_HALF_PI_DOUBLE;           
            const static double CMATH_PI_DOUBLE;                
            const static double CMATH_2PI_DOUBLE;               
            const static double CMATH_4PI_DOUBLE;               
//  ===========================================================================

            const static double CMATH_LN2_DOUBLE;               
            const static double CMATH_LN2_INV_DOUBLE;           
//  ===========================================================================

            const static float CMATH_LN2_FLOAT;                 
            const static float CMATH_LN2_INV_FLOAT;             
//  ===========================================================================

            const static double CMATH_ONE_THIRD_DOUBLE;         
            const static double CMATH_TWO_THIRDS_DOUBLE;        
//  ===========================================================================

            const static float CMATH_ONE_THIRD_FLOAT;           
            const static float CMATH_TWO_THIRDS_FLOAT;          
//  ===========================================================================

            static FORCEINLINE double minimum(const double a, const double b)
            {
                return min(a, b);
            }

            static FORCEINLINE float minimum(const float a, const float b)
            {
                return min(a, b);
            }

            static FORCEINLINE long minimum(const long a, const long b)
            {
                return min(a, b);
            }

            static FORCEINLINE double fastMinimum(const double a, const double b)
            {
                return 0.5 * (a + b - fabs(a - b));
            }

            static FORCEINLINE float fastMinimum(const float a, const float b)
            {
                return 0.5f * (a + b - fabsf(a - b));
            }

//  ===========================================================================

            static FORCEINLINE double maximum(const double a, const double b)
            {
                return max(a, b);
            }

            static FORCEINLINE float maximum(const float a, const float b)
            {
                return max(a, b);
            }

            static FORCEINLINE long maximum(const long a, const long b)
            {
                return max(a, b);
            }

            static FORCEINLINE double fastMaximum(const double a, const double b)
            {
                return 0.5 * (a + b + fabs(a - b));
            }

            static FORCEINLINE float fastMaximum(const float a, const float b)
            {
                return 0.5f * (a + b + fabsf(a - b));
            }

//  ===========================================================================

            static FORCEINLINE double fastModulus(const double x)
            {
                return x - floor(x);
            }

            static FORCEINLINE float fastModulus(const float x)
            {
                return x - floorf(x);
            }

//  ===========================================================================

            static FORCEINLINE double fastAtan(const double x)
            {
                return (x / (1.0 + 0.28 * (x * x)));
            }

            static FORCEINLINE float fastAtan(const float x)
            {
                return (x / (1.f + 0.28f * (x * x)));
            }

//  ===========================================================================

            static FORCEINLINE double clamp(const double min, const double max, const double x)
            {
                return 0.5 * (fabs(x - min) + (min + max) - fabs(x - max));
            }

            static FORCEINLINE double clamp(const double x)
            {
                return 0.5 * (fabs(x) + 1.0 - fabs(x - 1.0));
            }

            static FORCEINLINE float clamp(const float min, const float max, const float x)
            {
                return 0.5f * (fabsf(x - min) + (min + max) - fabsf(x - max));
            }

//  ===========================================================================

            static FORCEINLINE double cube(const double x)
            {
                return x * x * x;
            }

            static FORCEINLINE float cube(const float x)
            {
                return x * x * x;
            }

//  ===========================================================================

            static FORCEINLINE double square(const double x)
            {
                return x * x;
            }

            static FORCEINLINE float square(const float x)
            {
                return x * x;
            }

//  ===========================================================================

            static FORCEINLINE long fastDoubleToLong(double x)
            {
            #ifdef WIN32
                x = x + 103079215104;
            #else
                #ifndef __i386
                    x = x + 103079215104.0;
                #else
                    return (long)x;
                #endif
            #endif
                return ((long*)&x)[IMAN] >> 16; 
            }

            static FORCEINLINE long fastDoubleToLongFloor(const double x)
            {
            #ifdef WIN32
                long returnValue;
                __asm
                {
                    fld x
                    fistp returnValue
                }
                return returnValue;
            #else
                return (long)floor(x);
            #endif
            }

//  ===========================================================================
            
            static FORCEINLINE double iZero(const double y)
            {
                double s  = 1.0;
                double ds = 1.0;
                double d  = 0.0;
                do
                {
                    d  = d + 2.0;
                    ds = ds * (y * y) / (d * d);
                    s  = s + ds;
                }while(ds > 1e-7 * s);
                return s;
            }

            static FORCEINLINE double zeroethOrderBessel(const double x)
            {
                double bessel = 1.0;
                double term   = 0.0;
                long i        = 1;

                do
                {
                    // Compute the term
                    term = pow(0.5 * x, (double)i) / factorial((double)i);

                    // Add to bessel
                    bessel += (term * term);

                    // Increment to next
                    i++;
                }while (term > 0.000001 * bessel);

                // Return the bessel value
                return bessel;
            }

            static FORCEINLINE double factorial(const double x)
            {
                // Store some values
                double outputValue = 1.0;   
                double value       = x;

                // Compute factorial
                while (value > 1) 
                {
                    outputValue *= value--;
                }

                // Done :)
                return outputValue;
            }

//  ===========================================================================

//  ===========================================================================

            static FORCEINLINE int fastSign(const float x)
            {
            #ifdef WIN32
                return 1 + (((*(int *) &x) >> 31) << 1);
            #else
                if (x > 0)
                {
                    return 1;
                }
                else if (x < 0)
                {
                    return -1;
                }
                return 0;
            #endif
            }

            static FORCEINLINE int fastSign(const double x)
            {
            #ifdef WIN32
                return fastSign((float)x);
            #else
                if (x > 0)
                {
                    return 1;
                }
                else if (x < 0)
                {
                    return -1;
                }
                return 0;
            #endif
            }

//  ===========================================================================

            static FORCEINLINE double fastLog(const double x)
            {
            #ifdef WIN32
                return (double)fastLog((float)x);
            #else
                return logf(x);
            #endif
            }

            static FORCEINLINE float fastLog(const float x)
            {
            #ifdef WIN32
                return (fastLog2(x) * 0.69314718f);
            #else
                return log(x);
            #endif
            }

//  ===========================================================================

            static FORCEINLINE double fastAbsolute(const double x)
            {
            #ifdef WIN32
                return (double)fastAbsolute((float)x);
            #else
                return fabs(x);
            #endif
            }

            static FORCEINLINE float fastAbsolute(const float x)
            {
            #ifdef WIN32
                int i = ((*(int*)&x)&0x7fffffff);
                return (*(float*)&i);
            #else
                return fabsf(x);
            #endif
            }

//  ===========================================================================


            static FORCEINLINE double fastSquareRoot(const double x)
            {
            #ifdef WIN32
                double y, z, temp;
                unsigned long *tempPtr = ((unsigned long *)&temp)+1;

                temp = x;
                *tempPtr = (0xbfcd4600 - *tempPtr) >> 1;
                y = temp;
                z = x * 0.5;
                y *= 1.5 - y * y * z;
                y *= 1.5 - y * y * z;
                y *= 1.5 - y * y * z;
                y *= 1.5 - y * y * z;
                return y * x;
            #else
                return sqrt(x);
            #endif
            }

            static FORCEINLINE float fastSquareRoot(const float x)
            {
            #ifdef WIN32
                return (float)fastSquareRoot((double)x);
            #else
                return sqrtf(x);
            #endif
            }


            static FORCEINLINE float fastSquareRootVersion2(const float x)
            {
            #ifdef WIN32
                float   result;
                _asm
                {
                    mov eax, x
                    sub eax, 0x3f800000
                    sar eax, 1
                    add eax, 0x3f800000
                    mov result, eax
                }
                return result;
            #else
                return sqrtf(x);
            #endif
            }

//  ===========================================================================

            static FORCEINLINE double fastExponential(const double x)
            {
                return myPow2(x * CMATH_LN2_INV_DOUBLE);
            }


            static FORCEINLINE float myPow2(float i)
            {
                //float PowBodge=0.33971f;
                float x;
                float y=i-floorf(i);
                y=(y-y*y)*PowBodge;

                x=i+127-y;
                x*= shift23; //pow(2,23);
                *(int*)&x=(int)x;
                return x;
            }

            static FORCEINLINE double myPow2(double x)
            {
                // store address of float as long pointer
                long *p = (long*)(&x) + exponent_position;
                // truncated x for integer power of 2
                const long tx = fastDoubleToLongFloor(x);
                // float remainder of power of 2
                x -= static_cast<double>(tx);
                // sqr apporoximation of 2^x in [0, 1]
                //x  = 1.0 + x*(2.0/3.0 + x*1.0/3.0);
                x = 1.0 + x * (0.6 + x * 0.3);
                // add integer power of 2 to exponent
                *p += (tx<<20);
                return x;
            }

            static FORCEINLINE double fastExp2(const double x)
            {
            #ifdef WIN32
                int e;
                double ret;

                if (x >= 0)
                {
                    e = int (x);
                    ret = x - (e - 1);
                    ((*(1 + (int *) &ret)) &= ~(2047 << 20)) += (e + 1023) << 20;
                }
                else
                {
                    e = int (x + 1023);
                    ret = x - (e - 1024);
                    ((*(1 + (int *) &ret)) &= ~(2047 << 20)) += e << 20;
                }
                return (ret);
            #else
                return exp2(x);
            #endif
            }

            static FORCEINLINE float fastExp2(const float x)
            {
            #ifdef WIN32
                return (float)fastExp2((double)x);
            #else
                return exp2f(x);
            #endif
            }

            static FORCEINLINE double fastExp(const double x)
            {
                /*

                y = 2


                exp(x * ln(2)) == pow(2, x)

                I want a fast exponential function....

                pow(2, x) == exp(x * ln(2)
                pow(2, (x / ln2)) == exp(x)

                2^x = e^xln(2)




                
                // exp(x * ln(y)) == pow(y, x)

                */
                // Nameless union
                union
                {
                    double d;
                #ifdef WIN32
                    struct { int j, i; } n;
                #else
                    struct { int i, j; } n;
                #endif
                }eco;

                eco.n.i = (int)(1512775.3951951856938358403823062 * x) + (1072632447);  //(1072693248 - 60801);
                eco.n.j = 0;

                return eco.d;
            }

//  ===========================================================================

            static FORCEINLINE double fastPower4(double x)
            {
            #ifdef WIN32
                return (double)fastPower4((float)x);
            #else
                return x * x * x * x;
            #endif
            }

            static FORCEINLINE float fastPower4(float x)
            {
            #ifdef WIN32
                long *longPointer = (long *)(&x);
                long myLong       = *longPointer;
                myLong -= 0x3F800000l;  // Unbias
                myLong <<= 2;           // ^4
                myLong += 0x3F800000l;  // Rebias
                *longPointer = myLong;
                return x;
            #else
                return x * x * x * x;
            #endif
            }

            static FORCEINLINE double fastPower(const double b, const double e)
            {
                return exp(e * log(b));
            }

//  ===========================================================================

            static FORCEINLINE double fastSin(const double x)
            {
                const double xSquared = x * x;
                double returnValue    = 0.00761;
                returnValue          *= xSquared;
                returnValue          -= 0.16605;
                returnValue          *= xSquared;
                returnValue          += 1.0;
                returnValue          *= x;
                return returnValue;
            }

//  ===========================================================================

            static FORCEINLINE double tanhApproximation(const double x)
            {
                double a = fabs(x);
                a*=(6+a*(3+a));
                return ((x<0)?-1:1)*a/(a+12);
                //a = 6.0 + a * (6.0 + a * (3.0 + a));
                //return ((x < 0) ? -1 : 1) * (a - 6.0) / (a + 6.0);
                //return 1.0 - 2.0 / (fastExponential(2.0 * x) + 1.0);
            }

            static FORCEINLINE double fastTanh(const double x)
            {
                //return 1.0 - 2.0 / (exp(2.0 * x) + 1);
                return 1.0 - 2.0 / (fastExp(2.0 * x) + 1.0);
            }

            static FORCEINLINE float fastTanh(const float x)
            {
                return 1.f - 2.f / (expf(2.f * x) + 1.f);
            }

#ifdef WIN32

            static FORCEINLINE double fastTanh2(const double x)
            {
                double absolute = fabs(x);
                absolute *= (6 + absolute * (3 + absolute));
                return fastSign(x) * absolute / (absolute + 12);
            }

            static FORCEINLINE float fastTanh2(const float x)
            {
                float absolute = fabsf(x);
                absolute *= (6 + absolute * (3 + absolute));
                return fastSign(x) * absolute / (absolute + 12);
            }
#endif

//  ===========================================================================

        protected:

//  ===========================================================================
#ifdef WIN32

            static FORCEINLINE float fastLog2(float x)
            {
                int * const  exp_ptr = reinterpret_cast <int *> (&x);
                int val              = *exp_ptr;
                const int log_2      = ((val >> 23) & 255) - 128;
                val &= ~(255 << 23);
                val += 127 << 23;
                *exp_ptr = val;

                return (x + log_2);
            }
#endif

        };
    }
}
#endif  // End of CMathTools.hpp

---