前两天逛github看到一道很简单的面试题——如何不用库函数快速求出$\sqrt2$的值,精确到小数点后10位! 第一反应这不很简单嘛,大学数据结构课讲二分查找的时候老师还用这个作过示例。但转念一想,能做为大厂的面试题,背后绝对没有那么简单,因而我google了下,结果找到了更巧妙的数学方法,甚至发现了一件奇闻趣事…… 一道简简单单的面试题,不只能考察到候选人的编程能力,还能间接考察到候选人的数学素养,难怪不少大厂都会问这个。。。
回到正题,求$\sqrt2$究竟有多少种解法,咱们由简入难一步步来看下咱们是如何让计算机更快计算sqrt的。html
首先是稍微具有点编程和数据结构基础的人都能想到的二分查找,这里我再也不具体讲解思路,但仍是要编码测试下,主要是测试须要迭代多少次才能到达小数点后10位的精度。java
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。python
数学学得好的同窗确定知道牛顿迭代法,它是求解线性方程近似解的方法,由于有些方程没法快速求出精确解,只能尽量去逼近。git
回到咱们求解$\sqrt2$上,咱们能够构造出多项式方程f(x),使得$f(x)=0$的一个正根是$\sqrt2$,最简单的就是$f(x) = x^2 - 2= 0$,而后咱们就能够运用牛顿迭代求解它的根了。github
$$ 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__,提到这个常数还得提到一个游戏,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
上面说了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的近似求根算法。
下面这段代码能够返回输入值求根后的整数部分,但彻底不知道是什么原理。
_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...