2011年6月20日,星期一

快速近似对数,指数,幂和逆根

受启发 伊恩·斯蒂芬森(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-1ubuntu1libc 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-1ubuntu1libc 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-1ubuntu1libc 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-1ubuntu1libc 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 查找逆根比使用一般幂函数更准确。

买者自负

除了精度较低外,这些例程还不能处理浮点数 特殊价值 像标准例程一样。另外,如果您不适当地向这些例程提供负输入,则会发生奇怪的事情。

22条评论:

  1. 你知道沟日堆栈吗".693147"最终的线索是什么?将LN_2倒入安装在横梁上的管中,然后温度开关打开门。

    回复删除
  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

    回复删除
  3. 丹尼斯,我怀疑是上面介绍的fastexp的双精度变体,但是我还没有'请仔细阅读您提供的参考资料。

    我(此处未介绍)基于此处的fastexp和fastexp函数开发了fastsigmoid和fastssigmoid。有趣的是,由于S型(和tanh)如何使用指数,因此它们倾向于减小近似值的相对误差。鉴于参考文献是由神经网络驱动的,我可以看到他将如何以相对较低的精度摆脱困境。

    对于LDA,我们需要近似lgamma和digamma,我认为附加精度会有所帮助。

    回复删除
  4. 顺便说一句's是[Knuth vol。 1]:log_e + log_10。它'精确到99%。 (不是那样's any use to you!)

    回复删除
  5. 哇那's hard to beat.

    5(log_11-log_100)

    是更好的log_e近似值,但难记。

    回复删除
  6. 保罗,你好

    我对您的快速指数算法fastpow2有疑问。对于负输入,它似乎不太准确,例如对于输入-0.01,您的算法会产生1.505,而真实值应为0.99

    回复删除
  7. @HotLikeAToaster:我只是在我的机器上尝试过,然后得到:

    fastpow2(-0.01)= 0.993092

    回复删除
  8. 非常感谢!我在那里找到了丹尼斯提到的Schraudolph技术,但是通过使用有理多项式对残差进行更复杂的建模,可以明显改善其结果。我当时在想我该怎么做,但是您的这个程序似乎正是这样。这应该是众所周知的,并且随处可见!

    回复删除
  9. I'm试图使fastexp()以双精度工作,但收效甚微。
    I'd如果有人能指出正确的方向,我将感激不尽。
    谢谢!

    回复删除
    回覆
    1. @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'一旦设计出更高的精度,它可能不会比标准库快。

      删除
  10. 大约1.0时,fastlog和fastlog的近似值很差(相对误差>对于fastlog为13,对于fastlog更糟糕)

    回复删除
  11. 基于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;
    }

    回复删除
  12. 我们可以将其用于开源项目吗?

    回复删除
  13. 对于接近零的自变量,指数的幂级数收敛更快。因此,将您的参数除以2的幂以使其接近零,执行截断的幂级数,然后将结果平方适当的次数。我不'不知道这是否会击败股票Exp函数,但是值得一看。

    回复删除
    回覆
    1. It'自从我已经有一段时间了'我们已经考虑过这一点,但是位操作本质上是提取2的幂的倍数,因此,这可能比显式除法和平方更好。

      删除
  14. 感谢您的话题。我的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)
    最终功能日志



    回复删除
  15. FYI: "WARNING: cdn.mathjax.org has been retired. Check //www.mathjax.org/cdn-shutting-down/ for migration tips."

    您可能需要更新mathjax cdn,甚至切换到KaTex以获得更快的性能! ;)

    回复删除
  16. 您能否解释一下那些神奇数字((1<<23),121.2740838f,27.7280233f,4.85522568f,1.49012907f)中的快速近似值?

    我猜(1<<23)它是浮点尾数长度(23位)的偏移。我对吗?

    谢谢!

    回复删除
  17. 当对残差(log_2(1 + z)-z)进行有理函数逼近时,将显示幻数。

    是的,转变是提取Manitissa:表示的指数部分是"easy to log".

    回复删除