Linux 核心專題: 位元運算

執行人: NeedToDebugMyLife

Reviewed by jserv

說好的進度呢？

Reviewed by Denny0097

圖片顯示錯誤，可能是權限問題。

Reviewed by ginsengAttack

幾何平均是否可以用2項式定理實作。

TODO: 以定點數實作開平方根

說明

定點數

定點數是指用固定整數位數表達分數的格式，考慮一個 32 bits 的數，如果使用 Q16.16 格式，就代表這個數的前 16 位元用於表示整數部分，後 16 位元用於表示小數部分

例如:

十進位	二進位	定點數(Q16.16)
1.5	$(1.1{)}_2$	0000000000000001 1000000000000000

定點數運算

加法 減法

和整數加減法相同

乘法

考慮 2 浮點數

$f_1$、

$f_2$，分別對應 Q16.16 定點數

$n_1$、

$n_2$。
兩定點數可表示為

$$\begin{array}{r}{n}_1=f_1\times 2^{16}\\ n_2=f_2\times 2^{16}\end{array}$$
兩定點數相乘後得到

$$\begin{array}{rl}{n}_1\times n_2& =(f_1\times 2^{16})(f_2\times 2^{16})\\ & =f_1\times f_2\times 2^{32}\end{array}$$
但將兩浮點數相乘後再表示為定點數卻會得到

$$\begin{array}{rl}{n}_{f_{1\times 2}}& =(f_1\times f_2)\times 2^{16}\\ & \ne n_1\times n_2\end{array}$$
所以還需要再將

$n_1\times n_2$ 的結果除以

$2^{16}$ (右移 16 位元) 才會得到正確的乘法運算結果，即

$$n_1\times n_2=f_1\times f_2\times 2^{-16}$$

除法

和 乘法 類似
兩定點數相除後得到

$$\begin{array}{rl}{n}_1÷n_2& =(f_1\times 2^{16})÷(f_2\times 2^{16})\\ & =f_1÷f_2\end{array}$$
但將兩浮點數相除後再表示為定點數卻會得到

$$\begin{array}{rl}{n}_{f_{1÷2}}& =(f_1÷f_2)\times 2^{16}\\ & \ne n_1\times n_2\end{array}$$
所以還需要再將

$n_1\times n_2$ 的結果乘以

$2^{16}$ (左移 16 位元) 才會得到正確的乘法運算結果，即

$$n_1÷n_2=f_1÷f_2\times 2^{16}$$

牛頓法求根號

參考資料

牛頓法迭代公式:

$$x_{n+1}=x_n-\frac{f(x_n)}{f^{\prime }(x_n)}$$

考慮求

$\sqrt{n}:$

可以將原題目轉換為求

$f(x)=x^2-n$ 的根 (因為

$x=\sqrt{n}$ 為函式

$f(x)$ 的解)後得到:

$$\begin{array}{rl}f(x)& =x^2-n\\ f^{\prime }(x)& =2x\end{array}$$

接著將

$f(x)$ 和

$f^{\prime }(x)$ 代入到迭代公式
得到:

$$\begin{array}{rl}{x}_{n+1}& =x_n-\frac{x_n^2-n}{2x_n}\\ & =x_n-\frac{1}{2}\left(\frac{x_n^2-n}{x_n}\right)\\ & =x_n-\frac{1}{2}(x_n-\frac{n}{x_n})\\ & =\frac{1}{2}(x_n+\frac{n}{x_n})\end{array}$$

情境 1

輸入是 IEEE float (單精度)，輸出是 int 型態 (採用 LP64 data model)
針對上述情境，皆不可使用 FPU，只能藉由定點數予以計算，務必降低誤差並比較 glibc 的效能表現，提出改進方案

浮點數結構:

$n=2^S\ast M\ast 2^E=M\ast 2^E$

由於要求根號，所以只考慮 f 為正的情況 (即 S=0)

(負數以例外狀況處理)

開平方後得到

$$\begin{array}{rl}\sqrt{n}& =\sqrt{M\ast 2^E}\\ & =\sqrt{M}\ast 2^{\frac{E}{2}}\end{array}$$

如果

$E$ 是偶數，即

$E=2k$，則

$\sqrt{n}=\sqrt{M}\ast 2^k$
如果

$E$ 是奇數，即

$E=2k+1$，則

$\sqrt{n}=\sqrt{M}\ast \sqrt{2}\ast 2^k=\sqrt{2M}\ast 2^k$

可以發現只需要計算出

$\sqrt{M}$，就可以快速得到

$\sqrt{n}$

在計算前，我們還需要先得到浮點數的

$E$ 和

$M$

• E : 為第 23 位元到第 30 位元，所以需要右移 23 位，並且保留需要的 8 位。因為是讀出操作，所以需要再減去 127
• M : 為第 0 位元到第 22 位元，不需要做任何移位操作，但仍需要保留需要的 23 位

E = ((n >> 23) & 0xFF) - 127;
M = n & 0x7FFFFF;

因為浮點數不能直接做位移運算，所以需要使用 union 將浮點數的 32 位元表示讀取到一個 uint32_t 中

typedef union {
    float f;
    uint32_t i;
} fi_union;

fi_union u = f;

使用 u.i 求出

$E$ 和

$M$ (需要額外補上整數部分的 1)

int32_t E = ((u.i >> 23) & 0xFF) - 127;
uint32_t M = u.i & 0x007FFFFF | (1 << 23);

因為只能使用定點數做運算，所以計算

$\sqrt{M}$ 前還需要將

$M$ 轉換成定點數表示。但因為

$M$ 為 Q1.23，所以使用 Q32.32 格式來表示會比較適當 (即 32 位元為整數部分，後 32 位元為小數部分)。

uint64_t M_fixed = (uint64_t) M << 9 ;

前面提到，
如果

$E$ 是偶數，

$\sqrt{n}=\sqrt{M}\ast 2^k$
如果

$E$ 是奇數，

$\sqrt{n}=\sqrt{2M}\ast 2^k$

所以計算

$M$ 前還要檢查

$E$ 是否為奇數。

// multiply M by 2 if E is odd
if (E & 1)
    M_fixed <<= 1;

sqrt_M_fixed = sqrt_Newton_64(M_fixed);

使用牛頓法計算

$\sqrt{M}$

$$x_{n+1}=\frac{1}{2}(x_n+\frac{n}{x_n})$$

uint64_t sqrt_Newton_64(uint64_t fixed) {
  
    // using 1.5 for x0
    uint64_t x = (3ULL << 31);

    // iteration = 6
    for (int i = 0; i < 6; ++i) 
    {
        if (x == 0) 
            return 0;
        
        unsigned __int128 dividend = (unsigned __int128)fixed << 32;
        uint64_t div = dividend / x;
        x = (x + div) >> 1;
    }
    return x;
}

最後使用

$\sqrt{M}$ 和

$k=\frac{E}{2}$ 計算出

$\sqrt{n}$，並強制轉型為 int

uint64_t result = sqrt_M_fixed << (E >> 1);

因為回傳的值為 int，所以只需要保留前 32 位元

uint64_t result = (sqrt_M_fixed << (E >> 1)) >> 32;

整理後得到

uint64_t sqrt_M_fixed = sqrt_Newton_64(M_fixed);
uint64_t result = sqrt_M_fixed >> (32 - (E >> 1));

回傳結果並強制轉型為 int

情境2

輸入是 IEEE float，輸出也是 float 型態
針對上述情境，皆不可使用 FPU，只能藉由定點數予以計算，務必降低誤差並比較 glibc 的效能表現，提出改進方案

執行流程和情境1大致相同，差別在於最後要將 result 存為 float 而不是 int。

在得到計算結果 (sqrt_M_fixed) 後，需要從結果重新提取新的

$E$ 和

$M$

• M_new : 因為原本提取出的 M 是一個落在

  $[1.0,2.0)$ 區間的數，所以開根號後的值 (sqrt_M_fixed) 也必定落在
  $[1.0,2.0)$ 區間，滿足符點數 mantissa 的規範，因此可以直接將開根號結果的小數部分做為浮點數的 mantissa (需要對齊)，對齊部分為
  $32-23=9$ 個位元。
  
  
  對齊後去除多餘的整數部分 (1)
  ​​​​uint32_t M_new = (sqrt_M_fixed >> 9) & 0x7FFFFF;
  

• E_new : 前面提到，

  $\sqrt{M}$ 是一個落在
  $[1.0,2.0)$ 區間的數，以 1.OOO 來表示，再代回原算式:
  
  $$\begin{array}{rl}\sqrt{n}& =\sqrt{M}\ast 2^k\\ & =1.OOO\ast 2^k\end{array}$$
  可以發現這就是一個浮點數表示法，所以浮點數的 exponent 就等於
  $k$ (E >> 1)。但因為是寫入操作，所以要再加上 127
  ​​​​int32_t E_new = (E >> 1) + 127;
  

串接 E_new 和 M_new

fi_union result;
result.i = (E_new << 23) | M_new;

回傳 result.f

情境3

輸入和輸出都是 double (倍精度) 型態
針對上述情境，皆不可使用 FPU，只能藉由定點數予以計算，務必降低誤差並比較 glibc 的效能表現，提出改進方案

TODO: 計算幾何平均數

說明

幾何平均數 (Geometric Mean):

$$GM=\sqrt[n]{x_1\cdot x_2\cdot x_3\cdot \dots \cdot x_n}$$

計算的方法很簡單，但問題也很明顯。

• 如果數值都是小於 1 的小數，乘積就會趨近於 0，超出浮點數可表示的精度
• 如果數值都是很大數，乘積就會很容易超出最大可表示的範圍

要解決溢位的問題，可以將原本幾何平均數的公式取個對數
即:

$$ln(GM)=ln\left(\sqrt[n]{x_1\cdot x_2\cdot x_3\cdot \dots \cdot x_n}\right)=\frac{1}{n}\sum _{i=1}^{n}ln\left(x_i\right)$$
這樣就能將原本的連乘運算轉成連加運算。
但因為要求的值是

$GM$，不是

$ln(GM)$，所以還要再做個指數運算

$$GM=exp(ln(GM))=exp\left(\frac{1}{n}\sum _{i=1}^{n}ln\left(x_i\right)\right)$$
才能得出原本的

$GM$

情境1

給定若干 int 型態的資料，在 O(1) 空間複雜度的前提，計算幾何平均，要考慮到 overflow，並探討限制和改進方案

情境2

改為 float 型態，其餘同情境 1
針對上述情境，皆不可使用 FPU，只能藉由定點數予以計算，務必降低誤差並比較 glibc 的效能表現，提出改進方案

TODO: 以定點數實作快速傅利葉轉換 (FFT)

務必降低誤差並持續提出改進方案

TODO: 探討上述在 Linux 核心的應用

搭配課程教材，理解開平方根和幾何平均和 FFT 如何應用於 Linux 核心，予以探討其實作手法