前两天逛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

记住这个迭代次数34。

牛顿迭代

数学学得好的同学肯定知道牛顿迭代法,它是求解线性方程近似解的方法,因为有些方程无法快速求出精确解,只能尽可能去逼近。

回到我们求解 2 \sqrt2 2 上,我们可以构造出多项式方程f(x),使得 f ( x ) = 0 f(x)=0 f(x)=0的一个正根是 2 \sqrt2 2 ,最简单的就是 f ( x ) = x 2 − 2 = 0 f(x) = x^2 - 2= 0 f(x)=x22=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} x1=x0f(x0)f(x0)=x02x0x022

为了更直观些,我们图例展示下迭代两次的过程。
在这里插入图片描述
上图中黑色的曲线是 f ( x ) = x 2 − 2 f(x)=x^2-2 f(x)=x22,我们最终想要的是它和x轴的交点X,也就是 2 \sqrt2 2 的具体值。A点的坐标是(2,2),绿色的线是 f ( x ) = x 2 − 2 f(x)=x^2-2 f(x)=x22在点A处的切线,洋红色线是过A点的垂线。我们选的牛顿迭代初始点就是B,设B点的横坐标是 x 0 x_0 x0,按求导的定义,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(x0)f(x0),所以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)} x1=x0f(x0)f(x0)

怎么样,迭代一次之后就很逼近我们真正想要的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 1的神奇算法,当然它也可以用来求 2 \sqrt2 2

神奇的数字0x5f3759df

这首先要从一个诡异的常数说起—— 0x5f3759df,提到这个常数还得提到一个游戏,Quake-III Arena (雷神之锤3)。

这是90年代的经典游戏之一。该系列的游戏不但画面和内容不错,而且即使计算机配置低,也能极其流畅地运行。这款游戏后来在GPL协议下开源,github地址 https://github.com/id-Software/Quake-III-Arena,大家从中发现其中一个能够快速求解 1 x \frac{1}{\sqrt{x}} x 1的函数,代码如下。

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 RootFrom 百度百科

来看看这个算法的实际运行效果怎么样吧,下面是我修改后的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做个不严格的性能测试。这里说不严格是因为我只做 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

各种编程语言是如何实现sqrt?

上面说了3个解法,那你是不是也很好奇目前各种编程语言库函数里对sqrt是如何实现的? 我也很好奇,于是我们帮你们翻了jdk源码,发现它根本就不需要自己实现sqrt,毕竟sqrt在计算机领域是有个比较常用的计算,所以主流的CPU架构都提供了对sqrt的硬件支持,只需要一条汇编指令就可以了,在x86架构下sqrt可以直接用下面这条汇编指令

  sqrtsd %1, %0" : "=x" (res) : "xm" (x)

在Risc-v中可以可以用fsqrt.sfsqrt.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/%E7%89%9B%E9%A1%BF%E8%BF%AD%E4%BB%A3%E6%B3%95/10887580?fr=aladdin
[2] 卡马克快速平方根倒数算法: http://jcf94.com/2016/01/14/2016-01-14-carmack/
[3] Fast Inverse Square Root: http://read.pudn.com/downloads203/sourcecode/game/955182/3D%20Geometry%20Tuts/FastInverseSqrt.pdf
[4] From 百度百科: https://baike.baidu.com/item/0x5f375a86/10449453?fr=aladdin
[5] 这条汇编指令: https://code.woboq.org/userspace/glibc/sysdeps/x86_64/fpu/e_sqrt.c.html
[6] Rics-v中文手册: http://crva.ict.ac.cn/documents/RISC-V-Reader- Chinese-v2p1.pdf
[7] python math模块: https://github.com/python/cpython/blob/master/Modules/mathmodule.c

Logo

开放原子开发者工作坊旨在鼓励更多人参与开源活动,与志同道合的开发者们相互交流开发经验、分享开发心得、获取前沿技术趋势。工作坊有多种形式的开发者活动,如meetup、训练营等,主打技术交流,干货满满,真诚地邀请各位开发者共同参与!

更多推荐