UDN-企业互联网技术人气社区

板块导航

浏览  : 448
回复  : 1

[其它] 数学之美:平方根倒数速算法中的神奇数字 0x5f3759df①

[复制链接]
白青青的头像 楼主
发表于 2016-9-14 10:46:57 | 显示全部楼层 |阅读模式
  「平方根倒数速算法」中的神奇数字 0x5f3759df。Christian 在本文中探讨了该算法诸多有趣的地方,解释了它的原理,对指数在 -1 到 1 的情况进行了拓展,还有一些数学相关的新思路。
1.jpg

  0x5f3759df

  这篇文章讲述一个神奇的常数 0x5f3759df 和一个非常简洁的算法,平方根倒数速算法,也是这个常数的来源。

  来看看平方根倒数速算法:
  1. float FastInvSqrt(float x) {
  2.   float xhalf = 0.5f * x;
  3.   int i = *(int*)&x;         // evil floating point bit level hacking
  4.   i = 0x5f3759df - (i >> 1);  // what the fuck?
  5.   x = *(float*)&i;
  6.   x = x*(1.5f-(xhalf*x*x));
  7.   return x;
  8. }
复制代码

  这段代码可以快速计算
18.jpg

  的高精度近似值。

  现如今这个函数很出名,而它的流行是从出现在 2005 年的 Quake III Arena 的源码中开始的,起初归功于 John Carmack,但后来发现早在 Quake 之前,80 年代中期,Ardent Computer 公司就已使用,比 SGI 和 3dfx 还早,已能追溯到的最初作者是 Greg Walsh。上面具体的代码是 Quake 源码的改编版(Quake 源码是所有讨论的源头)。

  这篇文章探讨了这个算法诸多有趣的地方,解释了它的原理,对指数在 -1 到 1 的情况进行了拓展,还有一些数学相关的新思路。

  (本文确实包含许多数学知识,你可以把这些方程看做注释,无需阅读即可掌握本文要点,但如果你想要完整的阅读本文,验证我的观点的正确性,则需要关注这些方程。)

  为什么?

  为何需要计算倒平方根,甚至于值得实现一个怪异的算法来加速计算?因为它一直是 3D 编程计算中的一部分。在 3D 图形中,你使用平面法线,长度为 1 的三坐标向量,来表示光线和反射。你会使用很多平面法线,计算它们时需要对向量进行标准化。如何标准化一个向量呢?每个坐标分别除以向量的长度,因此,每个坐标需乘上
17.jpg

  计算
16.jpg

   相对来说开销很小,计算平方根并做除法开销很大。进入FastInvSqrt。

  如何做到的?

  这个函数究竟做了什么计算出了结果?它有 4 个主要步骤。首先它把浮点类型输入的二进制表示重定义为一个整形数字的二进制表示。
  1. int i = *(int*)&x;         // evil floating point bit level hack
复制代码

  运用得到的结果进行整形算术计算,得到我们所求值的近似值:
  1. i = 0x5f3759df - (i >> 1);  // what the fuck?
复制代码

  尽管结果不是近似值本身,但如果你把这个整形数字的二进制表示重定义为一个浮点数,那么就会得到近似值。所以代码做了与步骤一相反的变换回到了浮点数:
  1. x = *(float*)&i;
复制代码

  最终进行一次牛顿法迭代来提高精度。
  1. x = x*(1.5f-(xhalf*x*x));
复制代码

  这便得到 x 倒平方根的非常好的近似值。最后一步牛顿法相对比较直观,我便不多花时间讲解。关键步骤是第二步:对原始浮点数转化来的整形数进行算术运算,得到一个有意义的结果,下面重点关注这一步。

  怎么回事?

  这部分解释第二步背后的数学原理(下面推导的第一部分,至常数值的计算,最早由 McEniry 发现的)。

  在进入这一有趣的部分之前,我先快速解释下标准浮点数编码,我只讲解我要用的部分,完整的背景知识请参见维基百科。浮点数由三部分组成:符号,指数和尾数。这是单精度(32 位)浮点数二进制表示:
  1. s e e e e e e e e m m m m m m m m m m m m m m m m m m m m m m m
复制代码

  最高位是符号,接下来的 8 位是指数,底部 23 位是尾数。由于将要计算的平方根只对正数定义,所以从现在起假设符号位是 0 。

  当把一个浮点数看作一堆 0 、1 时,指数和底数便只是普通的正整数,没什么特别的。把它们记为 E 和 M(会经常使用)。另一方面,当把这些 0、1 解释为浮点值时,尾数则是介于 0 到 1 之间的值,全 0 意味着 0,全 1 是非常接近但略微小于 1 的数。并非将指数当做 8 位无符号整形使用,而是减去一个偏量 B,使之成为介于 -127 到 128 之间的有符号整形。把这些值的浮点解释记为 e 和 m。总之,我按照 McEnry 把整形解释相关的值用大写字母表示,把浮点解释相关的值用小写字母表示。

  这两种解释之间的转换很直接:
1.png

  对 32 位浮点数来说, L 是 223,B 是 127。给定 e 和 m,浮点数的值可计算得到:
14.jpg

  相应地整形解释的数是
13.jpg

  现在解释这个算法的所有基础几乎都有了,因此可以开始了,遗漏的部分会在解释过程中引入。给定输入 x,想要计算的是它的倒平方根,或者
12.jpg

  我们先对等式两边取以 2 为底的对数,至于原因很快就会知道:
11.jpg

  既然我们进行计算的值都是浮点数,则可以把 x 和 y 用各自的浮点组成替换:
10.jpg

  对数比较麻烦,幸运的是可以很轻松的摆脱它们,但在这之前,需暂停一下来解释其原理。

  在方程的两边都含有这样的项:
9.jpg

  其中 v 介于 0 到 1 之间,也仅仅当 v 在 0 到 1 之间时,这个函数和一条直线很接近:
8.jpg

  或者方程形式:
7.jpg

  其中 σ 是一个可选常数,尽管不能完美匹配但可以调节 σ 使得二者非常接近。通过它我们可以把上述包含对数的方程近似为一个线性方程,计算起来更轻松:
6.jpg

  到达关键的部分了!这时,用整形解释下的指数和尾数替代浮点解释下的表示进行计算就很方便了:
5.jpg

  通过少许步骤进行移项合并,就会得到很熟悉的形式(细节很乏味,可以跳过):
1.png

  走完最后一步之后,有趣的事情发生了:在这些杂乱的项中,方程两边都有整形解释的值:
3.jpg

  换句话说,y 的整形解释等于某个常数减去 x 的整形解释的一半,或者用 C 代码表示:
  1. i = K - (i >> 1);
复制代码

  K 看起来很熟悉对吧?

  剩下的工作就是找到这个常数。我们已经知道 B 和 L 的值,σ 的值还不知道。记住,σ 是对对数函数近似的调节因子,所以它的取值有些自由。我选取最初应用过的一个值,0.0450465,得到:
2.jpg

  想知道这个值的 16 进制表示?0x5f3759df。(当然理应是它,因为选的特定的 σ 得到的),这个常数不是一个位模式,这点你可能通过它被写成 16 进制形式而推知,它只是一个普通计算的结果化为整形的形式。

  但是正如 Knuth 所说,目前为止我们只是证明这个方法应该有效,还没检测。为了对这个近似的精确程度有直观的认识,用它的图线和精确的倒平方根的图线进行对比:
1.jpg

  这里输入的值介于 0 到 100 之间。很精确对吗?理应如此。这不仅仅是神奇,正如上面介绍的其推导,它的计算也很奇怪,但都是清晰、有意义的操作——浮点和整形之间的位转换

相关帖子

发表于 2016-9-14 14:05:57 | 显示全部楼层
赞一个
使用道具 举报

回复

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

关于我们
联系我们
  • 电话:010-86393388
  • 邮件:udn@yonyou.com
  • 地址:北京市海淀区北清路68号
移动客户端下载
关注我们
  • 微信公众号:yonyouudn
  • 扫描右侧二维码关注我们
  • 专注企业互联网的技术社区
版权所有:用友网络科技股份有限公司82041 京ICP备05007539号-11 京公网网备安1101080209224 Powered by Discuz!
快速回复 返回列表 返回顶部