文章目錄
- 二分查找
- 牛頓疊代
- 神奇的數字0x5f3759df
- 各種程式設計語言是如何實作sqrt?
- python中的_approximate_isqrt()
- 結語
- 引用連結
前兩天逛github看到一道很簡單的面試題——如何不用庫函數快速求出
的值，精确到小數點後10位！ 第一反應這不很簡單嘛，大學資料結構課講二分查找的時候老師還用這個做過示例。但轉念一想，能作為大廠的面試題，背後絕對沒有那麼簡單，于是我google了下，結果找到了更巧妙的數學方法，甚至發現了一件奇聞趣事…… 一道簡簡單單的面試題，不僅能考察到候選人的程式設計能力，還能間接考察到候選人的數學素養，難怪很多大廠都會問這個。。。
回到正題，求
究竟有多少種解法，我們由簡入難一步步來看下我們是如何讓計算機更快計算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。
牛頓疊代
數學學得好的同學肯定知道牛頓疊代法，它是求解線性方程近似解的方法，因為有些方程無法快速求出精确解，隻能盡可能去逼近。
回到我們求解上，我們可以構造出多項式方程f(x)，使得的一個正根是，最簡單的就是，然後我們就可以運用牛頓疊代求解它的根了。
為了更直覺些，我們圖例展示下疊代兩次的過程。
上圖中黑色的曲線是，我們最終想要的是它和x軸的交點X，也就是的具體值。A點的坐标是(2,2)，綠色的線是在點A處的切線，洋紅色線是過A點的垂線。我們選的牛頓疊代初始點就是B，設B點的橫坐标是，按求導的定義，AC點的斜率就是，AB的長度就是，知道了AB的長度，AC的斜率，BC的長度也就很好求了，就是，是以C點的橫坐标就是。怎麼樣，疊代一次之後就很逼近我們真正想要的X了，現在我們在C點把同樣的事情再做一遍，如下圖。
再一次疊代後得到了點E，點E比點C更加逼近X點，我把上圖中點E局部放大如下。
可以看到點E已經相當接近了，這隻是兩次疊代的效果，寫個代碼來看下到底需要疊代多少次就能達到精确到小數點後10位的精度。
具體代碼如下，這裡我取x的初始值為2.0 因為不可能大于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 到這裡文章進度條隻有一半，你是不是覺得我下面要講比牛頓疊代更好更快的方法？實際我目前沒有找到比牛頓疊代又好又快的算法了，但是我找到一個相關的故事，以及它引出的以犧牲精度換取速度求的神奇算法，當然它也可以用來求。
神奇的數字0x5f3759df
這首先要從一個詭異的常數說起—— 0x5f3759df，提到這個常數還得提到一個遊戲，Quake-III Arena (雷神之錘3)。
這是90年代的經典遊戲之一。該系列的遊戲不但畫面和内容不錯，而且即使計算機配置低，也能極其流暢地運作。這款遊戲後來在GPL協定下開源，github位址 https://github.com/id-Software/Quake-III-Arena，大家從中發現其中一個能夠快速求解的函數，代碼如下。
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做個不嚴格的性能測試。這裡說不嚴格是因為我隻做的測試，而且用的是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)