受启发
伊恩·斯蒂芬森(Ian Stephenson) ,我发现了一个更快更准确的$ \ log_2 $近似值,这反过来又提高了我对$ p $根的反向估计。所产生的技巧通常对于计算能力$ x ^ p $运作良好,这也是计算$ L_p $范数的一部分。还有一个快速的近似对数和指数,它们在广义上相关,因为这些运算在实践中被广泛地用于
表示概率.
更新资料 :这里的代码示例已过期;小问题已修复,您应该从
Google代码存储库 。 (
更新资料 2:Google代码已失效,但是
版本存在于github )
快速近似对数
快速的\ log_2 $近似值意味着对自然对数或编译时已知的任何其他基数的快速近似值,因为差是一个乘法。因此,考虑一个浮点数$ x $和值$ I_x $,它们将$ x $的字节表示形式解释为整数\ [
\ begin {aligned}
x&=(1 + m_x)2 ^ {e_x},\\
\ log_2(x)&= e_x + \ log_2(1 + m_x),\\
I_x&=(e_x + B)L + m_x L,\\
I_x / L-B&= e_x + m_x,\\
\ log_2(x)&= I_x / L-B + \ log_2(1 + m_x)-m_x。
\ end {aligned}
\] 因此,需要一种从$ x $中提取尾数并近似为[0,1)$中的$ z \ g_2(1 + z)-z $的尾数的方法。我采用了有理近似值,这似乎是准确性和速度的最佳结合点,
// WARNING: this code has been updated. Do not use this version.
// Instead, see http://code.google.com/p/fastapprox for the latest version.
inline 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);
}
inline float
fastlog (float x)
{
return 0.69314718f * fastlog2 (x);
}
如果您想要更快,更粗糙,则$ g(z)$的平均值为$ \ frac {-2 + \ log(8)} {\ log(4)} $
// WARNING: this code has been updated. Do not use this version.
// Instead, see http://code.google.com/p/fastapprox for the latest version.
inline float
fasterlog2 (float x)
{
union { float f; uint32_t i; } vx = { x };
float y = vx.i;
y *= 1.0 / (1 << 23);
return y - 126.94269504f;
}
inline float
fasterlog (float x)
{
return 0.69314718f * fasterlog2 (x);
}
时间和准确性
时序测试是通过编译完成的
-O3 -finline功能-fast-math,在运行64位Ubuntu Lucid(所以
gcc 4:4.4.3-1ubuntu1 和
libc 2.11.1-0ubuntu7.6)。我还测量了$ x \ in [1/100,10] $中的平均相对精度。
fastlog2 relative accuracy = 2.09352e-05
fastlog relative accuracy = 2.09348e-05
fasterlog2 relative accuracy = 0.0130367
fasterlog relative accuracy = 0.0130367
fastlog2 million calls per second = 160.141
fastlog million calls per second = 143.552
fasterlog2 million calls per second = 218.345
fasterlog million calls per second = 210.435
log2f million calls per second = 40.8511
快速近似指数
同样,2是最方便的基础。这里$ 2 ^ x $是我要近似的值,$ I_ {2 ^ x} $是与将$ 2 ^ x $的字节模式解释为整数相关的值。 \ [
\ begin {aligned}
2 ^ x&= 2 ^ {x-\ lfloor x \ rfloor} 2 ^ {\ lfloor x \ rfloor} \\
&=(1 + 2 ^ {x-\ lfloor x \ rfloor}-1)2 ^ {\ lfloor x \ rfloor},\\
I_ {2 ^ x}&=(\ lfloor x \ rfloor + B)L +(2 ^ {x-\ lfloor x \ rfloor}-1)L \\
&= L(x + B-1 + 2 ^ {x-\ lfloor x \ rfloor}-x + \ lfloor x \ rfloor)。
\ end {aligned}
\] 校正项的格式为$ g(z)= 2 ^ z-z $ for $ z \ in [0,1)$,并且我又一次从有理函数逼近中得到了最大的收益。这导致
// WARNING: this code has been updated. Do not use this version.
// Instead, see http://code.google.com/p/fastapprox for the latest version.
inline float
fastpow2 (float p)
{
union { float f; uint32_t i; } vp = { p };
int sign = (vp.i >> 31);
int w = p;
float z = p - w + sign;
union { uint32_t i; float f; } v = { (1 << 23) * (p + 121.2740838f + 27.7280233f / (4.84252568f - z) - 1.49012907f * z) };
return v.f;
}
inline float
fastexp (float p)
{
return fastpow2 (1.442695040f * p);
}
如果您想要更快,更粗糙,则$ g(z)$的平均值为$ \ frac {2-\ log(2)} {2 \ log(2)} $
// WARNING: this code has been updated. Do not use this version.
// Instead, see http://code.google.com/p/fastapprox for the latest version.
inline float
fasterpow2 (float p)
{
union { uint32_t i; float f; } v = { (1 << 23) * (p + 126.94269504f) };
return v.f;
}
inline float
fasterexp (float p)
{
return fasterpow2 (1.442695040f * p);
}
时间和准确性
时序测试是通过编译完成的
-O3 -finline功能-fast-math,在运行64位Ubuntu Lucid(所以
gcc 4:4.4.3-1ubuntu1 和
libc 2.11.1-0ubuntu7.6)。我还在[1/20,20] $中测量了$ p的平均相对准确度。当我测量相对准确度时,我尝试了每个$ p $和每个$ p $逆根(即$ -1 / p $)。
fastpow2 relative accuracy (positive p) = 1.58868e-05
fastexp relative accuracy (positive p) = 1.60712e-05
fasterpow2 relative accuracy (positive p) = 0.0152579
fasterexp relative accuracy (positive p) = 0.0152574
fastpow2 relative accuracy (inverse root p) = 1.43517e-05
fastexp relative accuracy (inverse root p) = 1.7255e-05
fasterpow2 relative accuracy (inverse root p) = 0.013501
fasterexp relative accuracy (inverse root p) = 0.0111832
fastpow2 million calls per second = 153.561
fastexp million calls per second = 143.311
fasterpow2 million calls per second = 215.006
fasterexp million calls per second = 214.44
expf million calls per second = 4.16527
快速近似功率
快速近似对数和快速近似指数的组合产生快速近似幂。
// WARNING: this code has been updated. Do not use this version.
// Instead, see http://code.google.com/p/fastapprox for the latest version.
inline float
fastpow (float x,
float p)
{
return fastpow2 (p * fastlog2 (x));
}
时间和准确性
时序测试是通过编译完成的
-O3 -finline功能-fast-math,在运行64位Ubuntu Lucid(所以
gcc 4:4.4.3-1ubuntu1 和
libc 2.11.1-0ubuntu7.6)。我还测量了$(x,p)\ [1/200,5] \ times [1/40,10] $中的平均相对准确度。当我测量相对准确度时,我尝试了每个$ p $和每个$ p $逆根(即$ -1 / p $)。
fastpow relative accuracy (positive p) = 0.000165618
fastpow relative accuracy (inverse root p) = 0.00011997
fastpow million calls per second = 58.8533
powf million calls per second = 8.44577
快速近似逆根
通过反根,我的意思是为$ p \ geq 1 $解决$ y = x ^ {-1 / p} $。这与快速近似功率不一样吗?是的,但是在这种情况下,以下方法相当准确且速度更快。初始近似值基于修改后的
反平方根hack ,\ [
I_y =-\ frac {1} {p} I_x + \ frac {p + 1} {p}(B-\ sigma)L。
\] 使用近似精炼
哈雷法 在$ f(y)= \ log(y)+ \ frac {1} {p} \ log(x)$,\ [
\ begin {aligned}
y_ {n + 1}&= y_n \ frac {1-\ delta} {1 + \ delta} \ approx y_n(1-\ delta)^ 2 \\
\ delta&= \ frac {\ log(2)} {2} \ left(\ log_2(y)+ \ frac {1} {p} \ log_2(x)\ right)
\ end {aligned}
\] $ \ log_2 $是使用上面的快速近似对数计算的。总共变成
// WARNING: this code has been updated. Do not use this version.
// Instead, see http://code.google.com/p/fastapprox for the latest version.
inline float
fasterinvproot (float x,
float p)
{
union { float f; uint32_t i; } v = { x };
unsigned int R = 0x5F375A84U; // R = (3/2) (B - sigma) L
unsigned int pR = ((2 * (p + 1)) / (3 * p)) * R;
unsigned int sub = v.i / p;
v.i = pR - sub;
return v.f;
}
inline float
fastinvproot (float x,
float p)
{
float y = fasterinvproot (x, p);
float log2x = fastlog2 (x);
float log2y = fastlog2 (y);
float err = 1.0 - 0.34657359f * (log2x / p + log2y);
return y * err * err;
}
注意:
fastinvproot 是改进(更快,更准确)的
fast_inverse_p_one 呈现在我的
以前的博客文章.
时间和准确性
时序测试是通过编译完成的
-O3 -finline功能-fast-math,在运行64位Ubuntu Lucid(所以
gcc 4:4.4.3-1ubuntu1 和
libc 2.11.1-0ubuntu7.6)。我还测量了$(x,p)\ [1/200,5] \ times [1/40,10] $中的平均相对准确度。当我测量相对准确度时,我尝试了每个$ p $和每个$ p $逆根(即$ -1 / p $)。
fastinvproot relative accuracy (positive p) = 0.000727901
fastinvproot relative accuracy (inverse root p) = 0.00300208
fastinvproot million calls per second = 82.7019
所以这有点令人困惑,但是这里的``正p''情况意味着要取一个反$ p $的根,而``反根p''的意思实际上是试图让这个例程陷入计算$ y = x ^ p $。因此,这些数字表示(出于我不完全理解的原因)
fastinvproot 查找逆根比使用一般幂函数更准确。
买者自负
除了精度较低外,这些例程还不能处理浮点数
特殊价值 像标准例程一样。另外,如果您不适当地向这些例程提供负输入,则会发生奇怪的事情。