求平方根倒数的算法

下面这个求1/\sqrt{x}的函数号称比直接调用sqrt库函数快4倍,来自游戏Quake III的源代码。

float InvSqrt (float x){
    float xhalf = 0.5f*x;
    int i = *(int*)&x;
    i = 0x5f3759df - (i>>1);
    y = *(float*)&i;
    y = y*(1.5f - xhalf*y*y);
    return x;
}

我们这里分析一下它的原理(指程序的正确性,而不是解释为何快)。

分析程序之前,我们必须解释一下float数据在计算机里的表示方式。一般而言,一个float数据x共32个bit,和int数据一样。其中前23位为有效数字M_x,后面接着一个8位数据E_x表示指数,最后一位表示符号,由于这里被开方的数总是大于0,所以我们暂不考虑最后一个符号位。此时

x=1.M_x 2^{E_x-127}

如果我们把计算机内的浮点数x看做一个整数I_x,那么

I_x = 2^{23}E_x+M_x

现在开始逐步分析函数。这个函数的主体有四个语句,分别的功能是:

int i = *(int*)&x; 这条语句把x转成i=I_x

i = 0x5f3759df - (i>>1); 这条语句从I_x计算I_{1/\sqrt{x}}

y = *(float*)&i; 这条语句将I_{1/\sqrt{x}}转换为1/\sqrt{x}

y = y*(1.5f - xhalf*y*y); 这时候的y是近似解;此步就是经典的牛顿迭代法。迭代次数越多越准确。

关键是第二步 i = 0x5f3759df - (i>>1); 这条语句从I_x计算I_{1/\sqrt{x}},原理:

y=1/\sqrt{x},用x=(1+m_x)2^{e_x}y=(1+m_y)2^{e_y}带入之后两边取对数,再利用近似表示\log_2(1+z)\sim z+\delta,算一算就得到

I_y = \frac{2}{3}(127-\delta)2^{23}-I_x/2

若取\delta=0.0450465679168701171875\frac{2}{3}(127-\delta)2^{23}就是程序里所用的常量0x5f3759df。至于为何选择这个\delta,则应该是曲线拟合实验的结果。

关于, »
  • Perfect Shuffle的算法 珍爱生命,远离政治。我们继续讨论算法。 2008/04/01补充:此算法有重大缺陷。详情请见留言部分。 一年前,我们讨论过一个算法问题,perfect shuffle,...
  • 有序矩阵的中位数算法 给定$$n\times n$$的实数矩阵,每行和每列都是递增的,求这$$n^2$$个数的中位数。 使用类似Tarjan的线性中位数的方法,每次找每列中位数,然后找中位数...
  • 理论计算机初步:算法和计算模型 下面是wikipedia上算法的定义: 算法是指完成一个任务所需要的具体步骤和方法。也就是说给定初始状态或输入数据,经过计算机程序的有限次运算...
  • 算法百科全书 - Encyclopedia of Algorithms Xie Xie推荐了一本今年出版的一本新书,Encyclopedia of Algorithms,全书1200页,涵盖各类有名问题的算法,概率算法,近似算法,量子算法等。 翻了一下,...
  • 一个算法面试题 & 面试题库 一个面试题,号称是微软的 输入$$a_1, a_2, ..., a_n, b_1, b_2, ..., b_n$$,如何在O(n)的时间,用O(1)的空间,将这个序列顺序改为$$a_1, b_1, ..., a_n, b_n$$。 刚一...
  • 一个简单图的三染色算法问题 注: 这个问题来自China Theory Week 2008的Open Problems Session。 我们知道在数学里证明一个东西的存在性的时候,有时候只证明了“存在性”,而且在证明过程...
  • 魔方里的数学 今天香港中文大学的Prof. Cai给我们上graph algorithm。第一节课上教我们玩魔方,先给每人发了一个。我喜欢这样的教学...
  • 递归算法的复杂度 递归算法的复杂度通常很难衡量,一般都认为是每次递归分支数的递归深度次方。但通常情况下没有这个大,如果我们可以保存每次子递归的结果的话...
  • 图同构问题属于P? 更新:证明的关键一步发现错误,作者更新了论文,结论甚至论文标题都改了(废话),新版本On the graph isomorphism problem。 提交论文到arxiv不需要审阅...
  • 编程的核心是数据结构,而不是算法 Rob Pike, 最伟大的 C 语言大师之一 , 在Notes on C Programming(英文原文)中从另一个稍微不同的角度表述了 Unix 的哲学: 你无法断定程序会在什么地方耗费运...
20条留言 -> 跳到留言表格
  • At 2008.09.19 14:35, wilderwein said:

    国王发布一求婚题,若应求者能在5次内猜出公主的生日,就把公主许配给他。在5次猜测当中,考官只会回答对或者否。
    这个题目你会做嘛

    • At 2008.09.19 15:29, zhiqiang said:

      看来我可以娶到公主呢。

      • At 2008.09.19 17:09, cuckoo said:

        日期可猜,月份就没办法了。你还是取不到。 :bigsmile

        • At 2008.09.21 10:59, 阿诺 2.0 said:

          第1问:日期加月数大于等于21么?

          是:公主生于9号之后。

          不是:公主生于1-9日.

          如果第一问是,

          第二问:公主的生日那月的最后一天的最后一位数等于他的他的生日日期第二位数么?

          是:那公主生于1*日.月份为1 3 5 7 8 10 12 .

          不是:公主生于2* 3*日月份是2 4 6 9 11.

          如果第一问不是,

          第二问:公主的生日那月的最后一天的最后一位数等于他的他的生日日期第数么?

          是:那公主生于1日.月份为1 3 5 7 8 10 12这几个数三问月数就出来了 .

          不是:公主生于2-9日 月份是2 4 6 9 11.同样三问之内可以问出的.

          如果第二问是:那公主生于1*日.月份为1 3 5 7 8 10 12 .

          接着问日期数的个位加月份数大于17么.

          大于 生于个位数6-9日.且月份为12.两问确定日期个位结束。

          如果小于 可以排除12月和个位是5以上和10月个位是8以上的

          那就再问个位加上月份数大于等于8么?

          大于: 1月个位8-9 3月6-9 7月1-2 8月1

          小于: 1月个位1-7 2月1-6 3月1-5 和7月1日

          • At 2008.09.21 17:00, zhiqiang said:

            没看懂 :thinking

            • At 2009.01.19 23:15, wilderwein said:

              土人,不过我也买看懂

          • At 2008.09.22 11:35, cuckoo said:

            第一问的推理就不对呀,比如公主生于1月10日。

        • At 2009.01.19 23:16, wilderwein said:

          大胆土鳖想当驸马爷

        • At 2008.09.19 18:39, stone said:

          2次就行了:)
          1.你将对我的第二个问题回答“是”吗?
          是-->公主的生日是1月1日吗?
          否-->公主的生日不是1月1日吧?

        • At 2008.09.19 17:10, cuckoo said:

          This is and interesting but old problem. :P

          • At 2008.09.22 01:30, e12e said:

            :-?

          • At 2008.09.20 01:10, 阿道夫 said:

            这个是根号的倒数 不是根号的逆运算 根号的逆运算是乘方

            • At 2008.09.23 00:22, bellqin said:

              回去读博士去了???还是不是一个人啊。。。汗

              • At 2008.10.19 10:18, Fdream said:

                你这问题后面跟的一个问题真是很汗……

              • At 2008.09.26 01:43, wilderwein said:

                楼上关于公主生日的兄弟的回答,有才。。。

                • At 2008.09.26 01:43, wilderwein said:

                  为此,奖励西域公主一个

                  • At 2008.09.30 23:42, 严酷的魔王 said:

                    这个求平方根倒数的源码有点火星了……
                    在Matrix67上看到过
                    还有一篇论文的下载地址
                    专门讲这个比较神奇的数0x5f3759df是怎么取到的……

                    • At 2008.10.01 09:39, zhiqiang said:

                      的确如此,我只不过详细逐步介绍了这段代码的作用和原理,其实也就是看的那篇论文,忘了给引用了。

                    • At 2008.10.18 21:13, hacker47 said:

                      还是以算法优化为先吧,
                      喜欢这些东西的(底层优化)的,
                      可以看hacker's delight.

                      • At 2009.04.12 17:48, eric said:

                        从理论上来讲 要五次绝对获得生日信息是不可能的 但是应该有使获得答案的概率最大的问法

                        (Required)
                        (Required, not published)

                          B | I | U | D | 添加链接 | 插入引用 | 插入代码 | 插入表情 | | + | ?
                        guest | 注册 | BBS | 管理 | English | 繁體 | https

                        阅微堂

                        zhiqiang's personal blog
                        Loading...
                        Loading...
                        Loading...