更新资料:这里的代码示例已过期;小问题已修复,您应该从 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所以这有点令人困惑,但是这里的``positive p''情况意味着要取一个逆$ p $ th的根,而``逆根p''实际上意味着试图让这个例程陷入计算$ y = x ^ p $。因此,这些数字表示(出于我不完全理解的原因) fastinvproot 查找逆根比使用一般幂函数更准确。
你知道沟日堆栈吗".693147"最终的线索是什么?将LN_2倒入安装在横梁上的管中,然后温度开关打开门。
回复删除保罗,您好-请阅读您的博客。
回复删除您对快速近似pow()的讨论使我想起了Schraudolph通过利用IEEE表示实现非常快的exp()的经典技巧。避风港'经过测试,看它在现代CPU上的性能如何's与您的方法比较,但我怀疑这会很有竞争力-'基本上只有3 FLOPS,并具有基于联合的一些巧妙的位调整,例如:
静态联合[
双d;
struct [int j,i; ] n;
生态
#定义EXPA(1048576 / LN(2))
#定义EXPC 60801
#定义EXP(y)(eco.n.i = EXPA *(y)+(1072693248-EXPC),eco.d)
See: http://cnl.salk.edu/~schraudo/pubs/Schraudolph99.pdf
丹尼斯,我怀疑是上面介绍的fastexp的双精度变体,但是我还没有'请仔细阅读您提供的参考资料。
回复删除我(此处未介绍)基于此处的fastexp和fastexp函数开发了fastsigmoid和fastssigmoid。有趣的是,由于S型(和tanh)如何使用指数,因此它们倾向于减小近似值的相对误差。鉴于参考文献是由神经网络驱动的,我可以看到他将如何以相对较低的精度摆脱困境。
对于LDA,我们需要近似lgamma和digamma,我认为附加精度会有所帮助。
顺便说一句's是[Knuth vol。 1]:log_e + log_10。它'精确到99%。 (不是那样's any use to you!)
回复删除哇那 's hard to beat.
回复删除5(log_11-log_100)
是更好的log_e近似值,但难记。
保罗,你好
回复删除我对您的快速指数算法fastpow2有疑问。对于负输入,它似乎不太准确,例如对于输入-0.01,您的算法会产生1.505,而真实值应为0.99
@HotLikeAToaster:我只是在我的机器上尝试过,然后得到:
回复删除fastpow2(-0.01)= 0.993092
非常感谢!我在那里找到了丹尼斯提到的Schraudolph技术,但是通过使用有理多项式对残差进行更复杂的建模,可以明显改善其结果。我当时在想我该怎么做,但是您的这个程序似乎正是这样。这应该是众所周知的,并且随处可见!
回复删除I'm试图使fastexp()以双精度工作,但收效甚微。
回复删除I'd如果有人能指出正确的方向,我将感激不尽。
谢谢!
@DatsMe:
删除It'尚不清楚是否要使用双精度但近似精度相同的精度更高的精度。
前者很容易(如果有点荒谬):垂头丧气地浮起,调用例程,然后将结果转换为两倍。
The latter will require adjusting the constants in the approximation. The mathematica notebook (//code.google.com/p/fastapprox/source/browse/trunk/fastapprox/tests/fastapprox.nb) contains the derivation which could be adapted to the proper mantissa, base, bitmask, etc. for doubles.
It'一旦设计出更高的精度,它可能不会比标准库快。
大约1.0时,fastlog和fastlog的近似值很差(相对误差>对于fastlog为13,对于fastlog更糟糕)
回复删除基于ln的反正切函数的幂级数展开,这也是对浮点数的非常快速的近似(对于double也是合理的,将log_t更改为double)。编译器可以对其向量化,从而大大加快了速度。
回复删除#包括
/ * M_LN2需要* /
#定义__USE_MISC 1
#包括
#包括
/ *预定义常量* /
#define i3 0.33333333333333333
#定义i5 0.2
#define i7 0.14285714285714285
#define i9 0.11111111111111111
#define i11 0.09090909090909091
#define i13 0.07692307692307693
#define i15 0.06666666666666667
/ *使用核心的类型别名以进行更精确的近似* /
typedef float log_t;
/ **
这是一个对数逼近函数。这是通过将尾数和
指数。然后使用标识ln z = 2 * arctan(z-1)/(z + 1)
** /
/ **小数的对数
它计算ln x,其中x<= 1.0
** /
静态内联log_t log1ds(log_t z){
log_t pz =(z-1)/(z + 1);
log_t pz2 = pz * pz;
log_t pz3 = pz2 * pz;
log_t pz5 = pz3 * pz2;
log_t pz7 = pz5 * pz2;
log_t pz9 = pz7 * pz2;
log_t pz11 = pz9 * pz2;
log_t pz13 = pz11 * pz2;
log_t pz15 = pz13 * pz2;
log_t y = 2 *(pz + pz3 * i3 + pz5 * i5 + pz7 * i7 + pz9 * i9 + pz11 * i11 + pz13 * i13 + pz15 * i15);
返回y;
}
#定义通话1
// 0000
#define TILL 1000000000
/ **自然对数** /
log_t logd(log_t x){
int exp;
log_t分数= frexpf((log_t)x,&exp);
log_t res = log1ds(分数);
返回res + exp * M_LN2;
}
我们可以将其用于开源项目吗?
回复删除是的,它是根据新的BSD许可发布的。
删除对于接近零的自变量,指数的幂级数收敛更快。因此,将您的参数除以2的幂以使其接近零,执行截断的幂级数,然后将结果平方适当的次数。我不'不知道这是否会击败股票Exp函数,但是值得一看。
回复删除It'自从我已经有一段时间了'我们已经考虑过这一点,但是位操作本质上是提取2的幂的倍数,因此,这可能比显式除法和平方更好。
删除感谢您的话题。我的C很烂(但是我是55岁的律师)。我使用了“ 2分频”的想法,试图为SIMPOL编写一个相当准确且不太慢的日志函数,'具有内置的日志功能。一旦介于0.5到1之间,我就得出了一些多项式近似值。
回复删除文字是SIMPOL,它是我发现输入错误和逻辑错误的第一种方法。
常数log10_2"0.30102999566398119521373889472449"
函数main()
n号
数字日志
对数= 0
整数e
e = 0
字符串信息,标题
anyvalue提示
讯息="Entry a number."
标题="log test"
提示=""
n = .toval(getuserinput(消息,标题,提示,错误= e),"",10)
log = 杰克_log(n)
结束函数.tostr(log,10)
函数jdk_log(数字n)
数字log,log10,log3
数字评估,转移
整数p
p = 0
评估= n
平移= 0
如果n< 1/2
评估=评估* 100
平移= 2
万一
评估时> 1
评估=评估/ 2
p = p +1
结束片刻
数字x,x2,x3,x4,x5,x6,x7,x8,x9,x10
编号test10,test3
x =评估
x2 = x * x
x3 = x2 * x
x4 = x3 * x
x5 = x4 * x
x6 = x5 * x
x7 = x6 * x
x8 = x7 * x
x9 = x8 * x
x10 = x9 * x
test10 =-(1436/1000)
test10 = test10 +(5541/1000 * x)
test10 = test10-(12204/1000 * x2)
test10 = test10 +(13987/1000 * x3)
test10 = test10 +(0304/1000 * x4)
test10 = test10-(17075/1000 * x5)
test10 = test10 +(4459/1000 * x6)
test10 = test10 +(30204/1000 * x7)
test10 = test10-(42185/1000 * x8)
test10 = test10 +(23255/1000 * x9)
test10 = test10-(4849/1000 * x10)
test3 =-(944/1000)
test3 = test3 +(1814/1000 * x)
test3 = test3-(1241/1000 * x2)
test3 = test3 +(371/1000 * x3)
log10 = p * .toval(log10_2,"",10)+轮(test10,1 / 1000000)-移位
log3 = p * .toval(log10_2,"",10)+轮(test3,1 / 1000000)-移位
log =轮((log10 + log3)/ 2,1 / 10000)
最终功能日志
好吧,我'd从未听说过SIMPOL!
删除FYI: "WARNING: cdn.mathjax.org has been retired. Check //www.mathjax.org/cdn-shutting-down/ for migration tips."
回复删除您可能需要更新mathjax cdn,甚至切换到KaTex以获得更快的性能! ;)
感谢您的注意。
删除您能否解释一下那些神奇数字((1<<23),121.2740838f,27.7280233f,4.85522568f,1.49012907f)中的快速近似值?
回复删除我猜(1<<23)它是浮点尾数长度(23位)的偏移。我对吗?
谢谢!
当对残差(log_2(1 + z)-z)进行有理函数逼近时,将显示幻数。
回复删除是的,转变是提取Manitissa:表示的指数部分是"easy to log".