前两天逛github看到一道很简单的面试题——如何不用库函数快速求出 2 \sqrt2
2
的值，精确到小数点后10位！ 第一反应这不很简单嘛，大学数据结构课讲二分查找的时候老师还用这个做过示例。但转念一想，能作为大厂的面试题，背后绝对没有那么简单，于是我google了下，结果找到了更巧妙的数学方法，甚至发现了一件奇闻趣事…… 一道简简单单的面试题，不仅能考察到候选人的编程能力，还能间接考察到候选人的数学素养，难怪很多大厂都会问这个。。。
回到正题，求 2 \sqrt2
2
究竟有多少种解法，我们由简入难一步步来看下我们是如何让计算机更快计算sqrt的。
二分查找
首先是稍微具备点编程和数据结构基础的人都能想到的二分查找，这里我不再具体讲解思路，但还是要编码测试下，主要是测试需要迭代多少次才能到达小数点后10位的精度。
public class BSearch { static double sqrt2 = 1.4142135624; static double delta = 0.0000000001; public static void main(String[] args) { double l = 1.0; double r = 2.0; int cnt = 0; while (l < r) { double mid = (l + r) / 2; if (Math.abs(l - sqrt2) < delta) { break; } if (mid * mid > 2.0) { r = mid; } else { l = mid; } cnt++; } System.out.println("经过" + cnt + "次迭代后得" + l); } }
经过34次迭代后得1.4142135623260401
1
记住这个迭代次数34。
牛顿迭代
数学学得好的同学肯定知道牛顿迭代法，它是求解线性方程近似解的方法，因为有些方程无法快速求出精确解，只能尽可能去逼近。
回到我们求解2 \sqrt2
2
上，我们可以构造出多项式方程f(x)，使得f ( x ) = 0 f(x)=0f(x)=0的一个正根是2 \sqrt2
2
，最简单的就是f ( x ) = x 2 − 2 = 0 f(x) = x^2 - 2= 0f(x)=x
2
−2=0，然后我们就可以运用牛顿迭代求解它的根了。
x 1 = x 0 − f ( x 0 ) f ′ ( x 0 ) = x 0 − x 0 2 − 2 2 x 0 x_{1}=x_{0}-\frac{f\left(x_{0}\right)}{f^{\prime}\left(x_{0}\right)} = x_0 - \frac{x_0^2-2}{2x_0}
x
1
=x
0
−
f
′
(x
0
)
f(x
0
)
=x
0
−
2x
0
x
0
2
−2
为了更直观些，我们图例展示下迭代两次的过程。
上图中黑色的曲线是f ( x ) = x 2 − 2 f(x)=x^2-2f(x)=x
2
−2，我们最终想要的是它和x轴的交点X，也就是2 \sqrt2
2
的具体值。A点的坐标是(2,2)，绿色的线是f ( x ) = x 2 − 2 f(x)=x^2-2f(x)=x
2
−2在点A处的切线，洋红色线是过A点的垂线。我们选的牛顿迭代初始点就是B，设B点的横坐标是x 0 x_0x
0
，按求导的定义，AC点的斜率就是f ′ ( x ) f^{\prime}(x)f
′
(x)，AB的长度就是f ( x ) f(x)f(x)，知道了AB的长度，AC的斜率，BC的长度也就很好求了，就是f ( x 0 ) f ′ ( x 0 ) \frac{f\left(x_{0}\right)}{f^{\prime}\left(x_{0}\right)}
f
′
(x
0
)
f(x
0
)
，所以C点的横坐标就是x 1 = x 0 − f ( x 0 ) f ′ ( x 0 ) x_1 = x_0 - \frac{f\left(x_{0}\right)}{f^{\prime}\left(x_{0}\right)}x
1
=x
0
−
f
′
(x
0
)
f(x
0
)
。
怎么样，迭代一次之后就很逼近我们真正想要的X了，现在我们在C点把同样的事情再做一遍，如下图。
再一次迭代后得到了点E，点E比点C更加逼近X点，我把上图中点E局部放大如下。
可以看到点E已经相当接近2 \sqrt2
2
了，这只是两次迭代的效果，写个代码来看下到底需要迭代多少次就能达到精确到小数点后10位的精度。
具体代码如下，这里我取x的初始值为2.0 因为2 \sqrt2
2
不可能大于2，我们知道这点就可以取个近似值，减少迭代次数。
public class NewtonMethod { static double sqrt2 = 1.4142135624; static double delta = 0.0000000001; public static void main(String[] args) { double x = 2.0; int cnt = 0; while (Math.abs(x - sqrt2) > delta){ x = x - (x*x - 2)/(2*x); cnt++; } System.out.println("经过" + cnt + "次迭代后得" + x); } }
经过4次迭代后得1.4142135623746899
只需要4次迭代就能达到二分34次迭代的效果了，确实明显快多了。这里再补充一点，实际上x的初始值可以取任意正数，但是会影响到性能，我尝试取1亿，最终需要30次迭代，不过还是比二分快。
为了客观对比下牛顿迭代和二分的性能差异，这里我还是用JMH压测下，结果如下，压测结果仅供参考。
Benchmark Mode Cnt Score Error Units Test.NewtonMethod thrpt 2 96292165.813 ops/s Test.bSearch thrpt 2 11856462.059 ops/s
到这里文章进度条只有一半，你是不是觉得我下面要讲比牛顿迭代更好更快的方法？实际我目前没有找到比牛顿迭代又好又快的算法了，但是我找到一个相关的故事，以及它引出的以牺牲精度换取速度求1 x \frac{1}{\sqrt{x}}
x
的神奇算法，当然它也可以用来求2 \sqrt2
神奇的数字0x5f3759df
这首先要从一个诡异的常数说起—— 0x5f3759df，提到这个常数还得提到一个游戏，Quake-III Arena (雷神之锤3)。
这是90年代的经典游戏之一。该系列的游戏不但画面和内容不错，而且即使计算机配置低，也能极其流畅地运行。这款游戏后来在GPL协议下开源，github地址 https://github.com/id-Software/Quake-III-Arena，大家从中发现其中一个能够快速求解1 x \frac{1}{\sqrt{x}}
x
的函数，代码如下。
float InvSqrt(float x) { float xhalf = 0.5f*x; int i = *(int*)&x; i = 0x5f3759df - (i >> 1); x = *(float*)&i; x = x*(1.5f - xhalf*x*x); return x; }
看完代码的我
其实这个算法就是牛顿迭代单次的近似解法，具体证明请看卡马克快速平方根倒数算法，它能以几十倍的速度优势秒杀其他算法，要知道几十年前的CPU速度可远不及现在的，速度就是绝对优势。这段代码看起来神奇，它的来源也很离奇。
开始大家都以为这个算法是游戏的开发者Carmack发现的，但后来调查发现，该算法在这之前就在计算机图形学的硬件与软件领域中有所应用，如SGI和3dfx就曾在产品中应用此算法，所以至今都无人知晓这个算法是谁发明的。
但传奇并没有结束。Lomont计算出结果以后非常满意，于是拿自己计算出的起始值和卡马克的神秘数字做比赛，看看谁的数字能够更快更精确求得平方根。结果是卡马克赢了… 谁也不知道卡马克是怎么找到这个数字的。最后Lomont怒了，采用暴力方法一个数字一个数字试过来，终于找到一个比卡马克数字要好上那么一丁点的数字，虽然实际上这两个数字所产生的结果非常近似，这个暴力得出的数字是0x5f375a86。Lomont为此写下一篇论文 Fast Inverse Square Root。 From 百度百科
来看看这个算法的实际运行效果怎么样吧，下面是我修改后的java代码，上面的C代码中有些操作java中并不支持，所以需要做些改动。
1.4145671304934657
1
当然这里精度就不行了，只精确到了0.001。
这里我把上面3种算法的精度都降低到0.001然后用JMH做个不严格的性能测试。这里说不严格是因为我只做2 \sqrt2
2
的测试，而且用的是java实现的，而且像是CarmackMethod的实现，可能因为java和c的运行机制的不同，性能会受很大影响，下面这个结果 仅供娱乐，看看就好。
Benchmark Mode Cnt Score Error Units
Test.CarmackMethod thrpt 2 111286742.330 ops/s
Test.bSearch thrpt 2 58705907.393 ops/s
Test.NewtonMethod thrpt 2 110232331.339 ops/s
