共计 5096 个字符，预计需要花费 13 分钟才能阅读完成。
前两天逛 github 看到一道很简略的面试题——如何不必库函数疾速求出 $\sqrt2$ 的值，准确到小数点后 10 位！第一反馈这不很简略嘛，大学数据结构课讲二分查找的时候老师还用这个做过示例。但转念一想，能作为大厂的面试题，背地相对没有那么简略，于是我 google 了下，后果找到了更奇妙的数学方法，甚至发现了一件奇闻趣事…… 一道简简单单的面试题，不仅能考查到候选人的编程能力，还能间接考查到候选人的数学素养，难怪很多大厂都会问这个。。。
回到正题，求 $\sqrt2$ 到底有多少种解法，咱们由简入难一步步来看下咱们是如何让计算机更快计算 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记住这个迭代次数 34。
牛顿迭代
数学学得好的同学必定晓得牛顿迭代法，它是求解线性方程近似解的办法，因为有些方程无奈疾速求出准确解，只能尽可能去迫近。
回到咱们求解 $\sqrt2$ 上，咱们能够结构出多项式方程 f(x)，使得 $f(x)=0$ 的一个正根是 $\sqrt2$，最简略的就是 $f(x) = x^2 – 2= 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}
$$
为了更直观些，咱们图例展现下迭代两次的过程。
上图中彩色的曲线是 $f(x)=x^2-2$，咱们最终想要的是它和 x 轴的交点 X，也就是 $\sqrt2$ 的具体值。A 点的坐标是(2,2)，绿色的线是 $f(x)=x^2-2$ 在点 A 处的切线，洋红色线是过 A 点的垂线。咱们选的牛顿迭代初始点就是 B，设 B 点的横坐标是 $x_0$，按求导的定义，AC 点的斜率就是 $f^{\prime}(x)$，AB 的长度就是 $f(x)$，晓得了 AB 的长度，AC 的斜率，BC 的长度也就很好求了，就是 $\frac{f\left(x_{0}\right)}{f^{\prime}\left(x_{0}\right)}$，所以 C 点的横坐标就是 $x_1 = x_0 – \frac{f\left(x_{0}\right)}{f^{\prime}\left(x_{0}\right)}$。
怎么样，迭代一次之后就很迫近咱们真正想要的 X 了，当初咱们在 C 点把同样的事件再做一遍，如下图。
再一次迭代后失去了点 E，点 E 比点 C 更加迫近 X 点，我把上图中点 E 部分放大如下。
能够看到点 E 曾经相当靠近 $\sqrt2$ 了，这只是两次迭代的成果，写个代码来看下到底须要迭代多少次就能达到准确到小数点后 10 位的精度。
具体代码如下，这里我取 x 的初始值为 2.0 因为 $\sqrt2$ 不可能大于 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到这里文章进度条只有一半，你是不是感觉我上面要讲比牛顿迭代更好更快的办法？理论我目前没有找到比牛顿迭代又好又快的算法了，然而我找到一个相干的故事，以及它引出的以就义精度换取速度求 $\frac{1}{\sqrt{x}}$ 的神奇算法，当然它也能够用来求 $\sqrt2$。
神奇的数字 0x5f3759df
这首先要从一个诡异的常数说起—— __0x5f3759df__，提到这个常数还得提到一个游戏，Quake-III Arena (雷神之锤 3)。
这是 90 年代的经典游戏之一。该系列的游戏岂但画面和内容不错，而且即便计算机配置低，也能极其流畅地运行。这款游戏起初在 GPL 协定下开源，github 地址 https://github.com/id-Software/Quake-III-Arena，大家从中发现其中一个可能疾速求解 $\frac{1}{\sqrt{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 中并不反对，所以须要做些改变。
public class CarmackMethod {public static void main(String[] args) {
float x = 2.0f;
float xhalf = 0.5f*x;
// int i = *(int*)&x;
int i= Float.floatToRawIntBits(x);
i = 0x5f3759df - (i >> 1);
// x = *(float*)&i;
x = Float.intBitsToFloat(i);
x = x*(1.5f - xhalf*x*x);
System.out.println(1.0/x);
}
}1.4145671304934657当然这里精度就不行了，只准确到了 0.001。
这里我把下面 3 种算法的精度都升高到 0.001 而后用 JMH 做个不严格的性能测试。这里说不严格是因为我只做 $\sqrt2$ 的测试，而且用的是 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各种编程语言是如何实现 sqrt?
下面说了 3 个解法，那你是不是也很好奇目前各种编程语言库函数里对 sqrt 是如何实现的？我也很好奇，于是咱们帮你们翻了 jdk 源码，发现它基本就不须要本人实现 sqrt，毕竟 sqrt 在计算机领域是有个比拟罕用的计算，所以支流的 CPU 架构都提供了对 sqrt 的硬件反对，只须要一条汇编指令就能够了，在 x86 架构下 sqrt 能够间接用上面这条汇编指令。
sqrtsd %1, %0":"=x"(res) :"xm" (x)在 Risc- v 中能够能够用 fsqrt.s 或fsqrt.d指令，Rics- v 中文手册
硬件能够在一个指令周期内实现一个数的开方，相比那些须要几十甚至成千盈百个 CPU 指令实现的各种算法而言，这速度差别不言而喻。实际上上，古代 CPU 在多外围、流水线、多发射、超标量……等技术的加持下，一般家用 CPU 都能够做到每秒百亿次的浮点运算。
__生存处处有惊喜__，当我关上 python math 模块的源码时，没有发现浮点数的求根(预计也是间接用的 CPU 级指令)，但我发现了一个更骚的对 64 位整数求根的操作，所以这里再补充介绍一个 python 的近似求根算法。
python 中的_approximate_isqrt()
上面这段代码能够返回输出值求根后的整数局部，但齐全不晓得是什么原理。
_approximate_isqrt(uint64_t n)
{uint32_t u = 1U + (n >> 62);
u = (u << 1) + (n >> 59) / u;
u = (u << 3) + (n >> 53) / u;
u = (u << 7) + (n >> 41) / u;
return (u << 15) + (n >> 17) / u;
}源码中有正文，这段 C 代码能够翻译为以下 python 代码，不过我仍旧看不懂。
def isqrt(n):
"""Return the integer part of the square root of the input."""
n = operator.index(n)
if n < 0:
raise ValueError("isqrt() argument must be nonnegative")
if n == 0:
return 0
c = (n.bit_length() - 1) // 2
a = 1
d = 0
for s in reversed(range(c.bit_length())):
# Loop invariant: (a-1)**2 < (n >> 2*(c - d)) < (a+1)**2
e = d
d = c >> s
a = (a << d - e - 1) + (n >> 2*c - e - d + 1) // a
return a - (a*a > n)结语
这篇博客从立题到实现经验了好几天的工夫，期间整顿思路、编码、绘图、查阅材料、批改欠缺总累计耗时近 8h。写作不易，如果文章对你有用欢送素质三连(点赞、珍藏加关注)。
援用链接
[1] 牛顿迭代法: https://baike.baidu.com/item/…
[2] 卡马克疾速平方根倒数算法: http://jcf94.com/2016/01/14/2…
[3] Fast Inverse Square Root: http://read.pudn.com/download…
[4] From 百度百科: https://baike.baidu.com/item/…
[5] 这条汇编指令: https://code.woboq.org/usersp…
[6] Rics- v 中文手册: http://crva.ict.ac.cn/documen… Chinese-v2p1.pdf
[7] python math 模块: https://github.com/python/cpy…
