面試題精選:求根號2簡單？高級算法你確定不會

時間 2020-07-25

原文原文鏈接

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

這首先要從一個詭異的常數提及—— __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)