經常使用的平方根算法詳解與實現

本文從屬於筆者的數據結構與算法系列文章。

SquareRoot

平方根計算一直是計算系統的經常使用算法，本文列舉出幾張簡單易懂的平方根算法講解與實現。其中Java版本的代碼參考這裏

Reference

• 計算平方根的算法

• Wiki-Methods of computing square roots

Babylonian:巴比倫算法/牛頓法

巴比倫算法可能算是最先的用於計算$sqrt{S}$的算法之一，由於其能夠用牛頓法導出，所以在不少地方也被成爲牛頓法。其核心思想在於爲了計算x的平方根，能夠從某個任意的猜想值g開始計算。在真實的運算中，咱們每每將g直接設置爲x，不過也能夠選擇其餘任何的正數值。那麼其計算的迭代過程爲:

1.若是猜想值g已經足夠接近於正確的平方根，算法結束，函數將g做爲結果返回。

2.若是猜想值g不夠精確，那麼使用g和x/g的平均值做爲新的猜想值。由於這兩個值中的一個小於確切的平方根，另外一個則大於確切的平方根，選擇平均值有助於你獲得一個更接近於正確答案的值。

3.把新的猜想值賦予給變量g，重複第一步的判斷。

綜上所述，其計算公式能夠表述爲:

其實現的參考代碼地址爲SquareRoots:

public double Babylonian() {

    double g = this.value;

    while (isApproximate(g)) {
        g = (g + this.value / g) / 2;
    }

    return g;

}

基於泰勒公式的級數逼近

微積分中的泰勒級數能夠表示爲:

在這個公式中，符號a表示某個常量，記號f'、f''和f'''表示函數f的一階、二階和三階導數，以此類推，這個公式稱爲泰勒公式，基於這個公式，咱們平方根公式的展開式爲:

根據該公式咱們能夠在必定精度內逼近真實值，不過這個公式仍然存在一個問題，便是公式的收斂問題。在泰勒級數展開中，平方根函數的公式當且僅當參數值位於一個有效範圍內時纔有效，在該範圍內計算趨於收斂。該範圍便是收斂半徑，當咱們對平方根函數用a=1進行計算時，泰勒級數公式但願x處於範圍:$0<x<2$之間。若是x在收斂半徑以外，則展開式中的項會愈來愈大，泰勒級數離答案也就愈來愈遠。爲了解決該問題，咱們能夠考慮當待開平方數大於4時以4去除它，最後將獲得的數乘以相同次數的2便可。其實現的參考代碼地址爲SquareRoots:

public double TSqrt() {

    //設置修正係數
    double correction = 1;

    //由於要對原值進行縮小,所以設置臨時值
    double tempValue = value;

    while (tempValue >= 2) {
        tempValue = tempValue / 4;
        correction *= 2;
    }

    return this.TSqrtIteration(tempValue) * correction;
}

private double TSqrtIteration(double value) {

    double sum = 0, coffe = 1, factorial = 1, xpower = 1, term = 1;

    int i = 0;

    while (Math.abs(term) > 0.000001) {

        sum += term;

        coffe *= (0.5 - i);

        factorial *= (i + 1);

        xpower *= (value - 1);

        term = coffe * xpower / factorial;

        i++;
    }

    return sum;

}

平方根倒數速算法

首先接收一個32位帶符浮點數，而後將之做爲一個32位整數看待，以將其向右進行一次邏輯移位的方式將之取半，並用十六進制「魔術數字」0x5f3759df減之，如此便可得對輸入的浮點數的平方根倒數的首次近似值；然後從新將其做爲浮點數，以牛頓法反覆迭代，以求出更精確的近似值，直至求出符合精確度要求的近似值。在計算浮點數的平方根倒數的同一精度的近似值時，此算法比直接使用浮點數除法要快四倍。此算法最先被認爲是由約翰·卡馬克所發明，但後來的調查顯示，該算法在這以前就於計算機圖形學的硬件與軟件領域有所應用，如SGI和3dfx就曾在產品中應用此算法。而就如今所知，此算法最先由Gary Tarolli在SGI Indigo的開發中使用。雖然說隨後的相關研究也提出了一些可能的來源，但至今爲止仍未能確切知曉此常數的起源。

其實現的參考代碼地址爲SquareRoots:

public double FastInverseSquareRoot() {

    double tempValue = value;

    double xhalf = 0.5d * tempValue;

    long i = Double.doubleToLongBits(tempValue);

    i = 0x5fe6ec85e7de30daL - (i >> 1);

    tempValue = Double.longBitsToDouble(i);

    tempValue = tempValue * (1.5d - xhalf * tempValue * tempValue);

    tempValue = this.value * tempValue;

    return tempValue;

}

Comparsion:比較

筆者創建了一個專門的單元測試類來比較上述算法的準確度與性能，代碼參考SquareRootsTest，首先在準確度與穩定性測試方面，這幾種算法都能達到較好地穩定性，其中平方根倒數速算法相對而言是較好。

@Test
public void testBabylonian() {

    for (int i = 0; i < 10000; i++) {
        Assert.assertEquals(2.166795861438391, squareRoots.Babylonian(), 0.000001);

    }

}

@Test
public void testTSqrt() {

    for (int i = 0; i < 10000; i++) {

        Assert.assertEquals(2.166795861438391, squareRoots.TSqrt(), 0.000001);

    }
}

@Test
public void testFastInverseSquareRoot() {

    for (int i = 0; i < 10000; i++) {

        Assert.assertEquals(2.1667948388864198, squareRoots.FastInverseSquareRoot(), 0.000001);

    }
}

而在性能測試方面，級數逼近的性能最差，巴比倫算法次之，平方根倒數速算法最好:

@Test
public void benchMark() {

    //巴比倫算法計時器
    long babylonianTimer = 0;

    //級數逼近算法計時器
    long tSqrtTimer = 0;

    //平方根倒數速算法計時器
    long fastInverseSquareRootTimer = 0;

    //隨機數生成器
    Random r = new Random();

    for (int i = 0; i < 100000; i++) {

        double value = r.nextDouble() * 1000;

        SquareRoots squareRoots = new SquareRoots(value);

        long start, stop;

        start = System.currentTimeMillis();

        squareRoots.Babylonian();

        babylonianTimer += (System.currentTimeMillis() - start);

        start = System.currentTimeMillis();

        squareRoots.TSqrt();

        tSqrtTimer += (System.currentTimeMillis() - start);

        start = System.currentTimeMillis();

        squareRoots.FastInverseSquareRoot();

        fastInverseSquareRootTimer += (System.currentTimeMillis() - start);

    }


    System.out.println("巴比倫算法:" + babylonianTimer);

    System.out.println("級數逼近算法:" + tSqrtTimer);

    System.out.println("平方根倒數速算法:" + fastInverseSquareRootTimer);


}


/**
結果爲:
巴比倫算法:17
級數逼近算法:34
平方根倒數速算法:7
**/