数学公式快速计算方法

#pragma once


#define fmaxf(a,b)         (((a) > (b)) ? (a) : (b))                   
typedef unsigned int        uint32_t;






//==========================================================================
//                           对数部分
//==========================================================================


///对数（底为2）函数
static inline float fastlog2 (float x)
{
union { float f; uint32_t i; } vx = { x };
union { uint32_t i; float f; } mx = { (vx.i & 0x007FFFFF) | 0x3f000000 };
float y = (float)vx.i;
y *= 1.1920928955078125e-7f;


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


///对数（底为2）函数
static inline float fasterlog2 (float x)
{
union { float f; uint32_t i; } vx = { x };
float y = (float)vx.i;
y *= 1.1920928955078125e-7f;
return y - 126.94269504f;
}


///对数函数 
static inline float fastlog (float x)
{
return 0.69314718f * fastlog2 (x);
}


///对数函数 
static inline float fasterlog (float x)
{
union { float f; uint32_t i; } 


vx = { x };


float y = (float)vx.i;
y *= 8.2629582881927490e-8f;


return y - 87.989971088f;
}


//==========================================================================
//                           指数部分
//==========================================================================


static inline float fastpow2 (float p)
{
float offset = (p < 0) ? 1.0f : 0.0f;
float clipp  = (p < -126) ? -126.0f : p;
int w   = (int)clipp;
float z = clipp - w + offset;


union 
{
uint32_t i; 
float f; 
} 


v = { (unsigned int) ( (1 << 23) * (clipp + 121.2740575f + 27.7280233f
/ (4.84252568f - z) - 1.49012907f * z) ) };


return v.f;
}


static inline float fasterpow2 (float p)
{
float clipp = (p < -126) ? -126.0f : p;


union 
{ 
uint32_t i;
float f; 
}


v = { (unsigned int) ( (1 << 23) * (clipp + 126.94269504f) ) };


return v.f;
}




///exp指数
static inline float fastexp (float p)
{
return fastpow2 (1.442695040f * p);
}


///exp指数
static inline float fasterexp (float p)
{
return fasterpow2 (1.442695040f * p);
}


static inline float fastpow (float x, float p)
{
return fastpow2 (p * fastlog2 (x));
}


static inline float fasterpow (float x,float p)
{
return fasterpow2 (p * fasterlog2 (x));
}




//==========================================================================
//                           三角函数部分(-pi, pi)
//==========================================================================


///超越正弦函数
static inline float fastsinh (float p)
{
return 0.5f * (fastexp (p) - fastexp (-p));
}


///超越正弦函数
static inline float fastersinh (float p)
{
return 0.5f * (fasterexp (p) - fasterexp (-p));
}


///超越余弦函数 
static inline float fastcosh (float p)
{
return 0.5f * (fastexp (p) + fastexp (-p));
}


///超越余弦函数 
static inline float fastercosh (float p)
{
return 0.5f * (fasterexp (p) + fasterexp (-p));
}


///超越正切函数
static inline float fasttanh (float p)
{
return -1.0f + 2.0f / (1.0f + fastexp (-2.0f * p));
}


///超越正切函数
static inline float fastertanh (float p)
{
return -1.0f + 2.0f / (1.0f + fasterexp (-2.0f * p));
}


///正弦函数 
static inline float fastsin (float x)
{
static const float fouroverpi = 1.2732395447351627f;
static const float fouroverpisq = 0.40528473456935109f;
static const float q = 0.78444488374548933f;
union { float f; uint32_t i; } p = { 0.20363937680730309f };
union { float f; uint32_t i; } r = { 0.015124940802184233f };
union { float f; uint32_t i; } s = { -0.0032225901625579573f };


union { float f; uint32_t i; } vx = { x };
uint32_t sign = vx.i & 0x80000000;
vx.i = vx.i & 0x7FFFFFFF;


float qpprox = fouroverpi * x - fouroverpisq * x * vx.f;
float qpproxsq = qpprox * qpprox;


p.i |= sign;
r.i |= sign;
s.i ^= sign;


return q * qpprox + qpproxsq * (p.f + qpproxsq * (r.f + qpproxsq * s.f));
}


///正弦函数 
static inline float fastersin (float x)
{
static const float fouroverpi = 1.2732395447351627f;
static const float fouroverpisq = 0.40528473456935109f;
static const float q = 0.77633023248007499f;
union { float f; uint32_t i; } p = { 0.22308510060189463f };


union { float f; uint32_t i; } vx = { x };
uint32_t sign = vx.i & 0x80000000;
vx.i &= 0x7FFFFFFF;


float qpprox = fouroverpi * x - fouroverpisq * x * vx.f;


p.i |= sign;


return qpprox * (q + p.f * qpprox);
}


///正弦函数 
static inline float fastsinfull (float x)
{
static const float twopi = 6.2831853071795865f;
static const float invtwopi = 0.15915494309189534f;


int k = (int)(x * invtwopi);
float half = (x < 0) ? -0.5f : 0.5f;
return fastsin ((half + k) * twopi - x);
}


///正弦函数 
static inline float fastersinfull (float x)
{
static const float twopi = 6.2831853071795865f;
static const float invtwopi = 0.15915494309189534f;


int k = (int)(x * invtwopi);
float half = (x < 0) ? -0.5f : 0.5f;
return fastersin ((half + k) * twopi - x);
}


///余弦函数
static inline float fastcos (float x)
{
static const float halfpi = 1.5707963267948966f;
static const float halfpiminustwopi = -4.7123889803846899f;
float offset = (x > halfpi) ? halfpiminustwopi : halfpi;
return fastsin (x + offset);
}


///余弦函数
static inline float fastercos (float x)
{
static const float twooverpi = 0.63661977236758134f;
static const float p = 0.54641335845679634f;


union { float f; uint32_t i; } vx = { x };
vx.i &= 0x7FFFFFFF;


float qpprox = 1.0f - twooverpi * vx.f;


return qpprox + p * qpprox * (1.0f - qpprox * qpprox);
}


///余弦函数
static inline float fastcosfull (float x)
{
static const float halfpi = 1.5707963267948966f;
return fastsinfull (x + halfpi);
}


///余弦函数
static inline float fastercosfull (float x)
{
static const float halfpi = 1.5707963267948966f;
return fastersinfull (x + halfpi);
}


///正切函数
static inline float fasttan (float x)
{
static const float halfpi = 1.5707963267948966f;
return fastsin (x) / fastsin (x + halfpi);
}


///正切函数
static inline float fastertan (float x)
{
return fastersin (x) / fastercos (x);
}


///正切函数
static inline float fasttanfull (float x)
{
static const float twopi = 6.2831853071795865f;
static const float invtwopi = 0.15915494309189534f;


int k = (int)(x * invtwopi);
float half = (x < 0) ? -0.5f : 0.5f;
float xnew = x - (half + k) * twopi;


return fastsin (xnew) / fastcos (xnew);
}


///正切函数
static inline float fastertanfull (float x)
{
static const float twopi = 6.2831853071795865f;
static const float invtwopi = 0.15915494309189534f;


int k = (int)(x * invtwopi);
float half = (x < 0) ? -0.5f : 0.5f;
float xnew = x - (half + k) * twopi;


return fastersin (xnew) / fastercos (xnew);
}




//==========================================================================
//                     sigmoid函数部分（1/(1 + exp(x)）
//==========================================================================


static inline float fastsigmoid (float x)
{
return 1.0f / (1.0f + fastexp (-x));
}


static inline float fastersigmoid (float x)
{
return 1.0f / (1.0f + fasterexp (-x));
}


//==========================================================================
//                           朗伯W函数部分
//==========================================================================


static inline float fastlambertw (float x)
{
static const float threshold = 2.26445f;


float c = (x < threshold) ? 1.546865557f : 1.0f;
float d = (x < threshold) ? 2.250366841f : 0.0f;
float a = (x < threshold) ? -0.737769969f : 0.0f;


float logterm = fastlog (c * x + d);
float loglogterm = fastlog (logterm);


float minusw = -a - logterm + loglogterm - loglogterm / logterm;
float expminusw = fastexp (minusw);
float xexpminusw = x * expminusw;
float pexpminusw = xexpminusw - minusw;


return (2.0f * xexpminusw - minusw * (4.0f * xexpminusw - minusw * pexpminusw)) /
(2.0f + pexpminusw * (2.0f - minusw));
}


static inline float fasterlambertw (float x)
{
static const float threshold = 2.26445f;


float c = (x < threshold) ? 1.546865557f : 1.0f;
float d = (x < threshold) ? 2.250366841f : 0.0f;
float a = (x < threshold) ? -0.737769969f : 0.0f;


float logterm = fasterlog (c * x + d);
float loglogterm = fasterlog (logterm);


float w = a + logterm - loglogterm + loglogterm / logterm;
float expw = fasterexp (-w);


return (w * w + expw * x) / (1.0f + w);
}


///朗伯W函数z=W(z).exp(W(z))
static inline float fastlambertwexpx (float x)
{
static const float k = 1.1765631309f;
static const float a = 0.94537622168f;


float logarg = fmaxf (x, k);
float powarg = (x < k) ? a * (x - k) : 0;


float logterm = fastlog (logarg);
float powterm = fasterpow2 (powarg);  // don't need accuracy here


float w = powterm * (logarg - logterm + logterm / logarg);
float logw = fastlog (w);
float p = x - logw;


return w * (2.0f + p + w * (3.0f + 2.0f * p)) /(2.0f - p + w * (5.0f + 2.0f * w));
}


///朗伯W函数z=W(z).exp(W(z))
static inline float fasterlambertwexpx (float x)
{
static const float k = 1.1765631309f;
static const float a = 0.94537622168f;


float logarg = fmaxf (x, k);
float powarg = (x < k) ? a * (x - k) : 0;


float logterm = fasterlog (logarg);
float powterm = fasterpow2 (powarg);


float w = powterm * (logarg - logterm + logterm / logarg);
float logw = fasterlog (w);


return w * (1.0f + x - logw) / (1.0f + w);
}






//==========================================================================
//                           对数Γ函数（伽玛函数）部分
//==========================================================================


///Log Gamma函数
static inline float fastlgamma (float x)
{
float logterm = fastlog (x * (1.0f + x) * (2.0f + x));
float xp3 = 3.0f + x;


return - 2.081061466f - x + 0.0833333f / xp3 - logterm + (2.5f + x) * fastlog (xp3);
}


///Log Gamma函数
static inline float fasterlgamma (float x)
{
return - 0.0810614667f 
- x
- fasterlog (x)
+ (0.5f + x) * fasterlog (1.0f + x);
}


///双伽马函数
static inline float fastdigamma (float x)
{
float twopx = 2.0f + x;
float logterm = fastlog (twopx);


return (-48.0f + x * (-157.0f + x * (-127.0f - 30.0f * x))) / 
(12.0f * x * (1.0f + x) * twopx * twopx) + logterm;
}


///双伽马函数
static inline float fasterdigamma (float x)
{
float onepx = 1.0f + x;


return -1.0f / x - 1.0f / (2 * onepx) + fasterlog (onepx);
}