Linux 核心快速運算 √x

整數 N 怎麼開平方根？

最簡單的做法，我們可以切一半取平方，然後比大小，夾取範圍。一直做二分逼近，直到我們可以接受的誤差。這樣做使用了大量乘除法，不僅了速度慢，也毫無「品味」。

讓我們來看看 Linux 核心如何實現？

unsigned long int_sqrt(unsigned long x)
{
	unsigned long b, m, y = 0;

	if (x <= 1)
		return x;

	m = 1UL << (__fls(x) & ~1UL);
    while (m != 0) {
        b = y + m;
        y >>= 1;
		if (x >= b) {
			x -= b;
			y += m;
		}
		m >>= 2;
	}

	return y;
}

這品味實在有點難以消化。

讓我們娓娓道來。

數值分析

假設我們知道答案 √N = y，也就是 y² = N。若將 y 表示為 n 項常數的和，則可以寫成 y_n = a₁ + a₂ + … + a_n。也就是說，我們可以透過從 a₁ 猜到 a_n 以求得 N 的平方根。這樣做在每次計算到前 m 項時，都必須得重新計算一次 y_m²，並與 N 做比較，以求得下一項 a_m+1 的數值。

我們試著將前 m 項的平方展開，看看其平方和有沒有規律。

• \(y_{1}^{2} = (0 + a_{1}^{2})\)
• \(y_{2}^{2} = (0 + a_{1}^{2}) + (2( a_{1} ) + a_{2})a_{2}\)
• \(y_{3}^{2} = (0 + a_{1}^{2}) + (2( a_{1} ) + a_{2})a_{2} + (2(a_{1} + a_{2}) + a_{3})a_{3}\)
• \(y_{m}^{2} = \text{…}\)

可以觀察出：\(y_{m}^2 = y_{m-1}^2 + Y_{m}\)，其中

若我們考慮計算到第 m 項的平方和與 N 的誤差 X_m，也就是 X_m = N – y_m²。隨著 m 越接近 n，這項誤差就會越小。並且可依據前面的歸納，得出底下規律：

而我們的目標，就是使得當 m = n 時，X_n 可以逼近 0。

上述式子可以用 n 位元 10 進位數字來替換。此時 y_n = b₁10^n-1 + b₂10^n-1 + … + b_n-110¹ + b_n ，其中 b_i 為 1 ～ 9 自然數 。這樣 Y_m 就會變成：

手算開平方

底下以求 √15200 為例我們一項一項試看看。當 N 等於 15200 時，由於 4 位數的平方就超過，n 取 3 即可。

• a₁ = b₁·100。取 b₁ 使得 (b₁·100)² <= N < ((b₁ + 1)·100)²。可得 b₁ = 1。
  • 此時 y₁² = a₁² = 100² = 10000。與 N 的誤差 X₁ = 5200。
• X₂ = X₁ – Y₂ = 5200 – (20·1 + b₂)·b₂。取 b₂ 使得 X₂ 為最小正數。可得 b₂ = 2。
• X₃ = X₂ – Y₃ = 800 – (20·12 + b₃)·b₃。取 b₃ 使得 X₃ 為最小正數。可得 b₃ = 3。
• 故所求 y_n = 1·100 + 2·10 + 3 = 123。且誤差為 X₃ = 71。

有點難懂？底下用類似長除法的視覺輔助，帶大家來手算開平方看看。

求 √15200。

      1  2  3  <- yn
    √ 1 52 00
      1        <- 1002
X1 ->   52
        44     <- (20·1 + b2)·b2
X2 ->    8 00
         7 29  <- (20·12 + b3)·b3
X2 ->      71

二進位簡化

由於電腦使用二進位的每一位元都是由 0 或者 1。從第 n 位元找到第 m 位元的表示式就可以寫成 y = b_n + b_n-1 + … + b_m+1 + b_m ，其中 b_m 為 2^m 或 0。則前述公式可改寫成：

欲確定 b_m 為 2^m 或 0，一樣可以透過逐項減少與 N 的誤差 \(X_{m} = X_{m+1} – Y_{m}\)。若X_m 仍大於等於 0，則為 2^m 反之為 0 。若令 \(P_{m+1} = \sum_{i=n}^{m+1}b_{i}\) 則此時 \(Y_{m} = 2 P_{m+1} \cdot 2^{m}+ 2^{2m}\) 且我們將 Y_m 的右邊兩項分別令為 c_m、d_m，可得：

分析 c_m、d_m 與下一項（第 m – 1 項）的關係，可得：

留意所求 y = P₀ = c_-1。

將以上數值分析寫成程式。

unsigned long int_sqrt(unsigned long x)
{
    unsigned long d;
    unsigned long c = 0; // c_n
    if (x <= 1)
        return x;
    
    // d_n is the highest bit we need to check
    // such that 2^2n < x
    d = 1UL << (__fls(x) & ~1UL);

    // for each d_m where d_m is d_n ... d_0
    while (d != 0) {
        if (x >= c + d) {      // if X_m = X_(m+1) - Ym ≥ 0, then b_m = 2^m
            x -= c + d;        //   store X_m
            c = (c >> 1) + d;  //   store c_(m-1)
        } else {               // else b_m = 0
            c = c >> 1;        //   store c_(m-1)
        }
        d >>= 2;               // store d_(m-1)
    }
    return c;                  // c_(-1)
}

終於搞定！

回頭對照一開始 Linux 核心的版本，發現還有再進一步簡化，以減少些許重複計算。實在歎為觀止。

後記

也許自己不是資工背景，光是了解數值運算方程就花了我好長時間。好不容易理解了方程，一變成程式碼，腦袋又卡住，無法一一對應前面結論，故撰文重新整理。除對齊所有變數之外，也補充了一些當初自己卡住的細節及說明。

References

• 從 √2 的存在談開平方根的快速運算 – HackMD
• Methods of computing square roots – Wikipedia