# Linux 核心專題: 位元運算 > 執行人: NeedToDebugMyLife ### Reviewed by `jserv` 說好的進度呢？ ### Reviewed by `Denny0097` >圖片顯示錯誤，可能是權限問題。 ### Reviewed by `ginsengAttack` 幾何平均是否可以用2項式定理實作。 ## TODO: 以定點數實作開平方根 ### `說明` #### 定點數 定點數是指用固定整數位數表達分數的格式，考慮一個 32 bits 的數，如果使用 Q16.16 格式，就代表這個數的前 16 位元用於表示整數部分，後 16 位元用於表示小數部分 ![image](https://hackmd.io/_uploads/HJtSfHSBxe.png) 例如: | 十進位 | 二進位 | 定點數(Q16.16) | | ------ |:---------:|:---------------------------------:| | 1.5 | $(1.1)_2$ | 0000000000000001 1000000000000000 | #### 定點數運算 ##### `加法` `減法` 和整數加減法相同 ##### `乘法` 考慮 2 浮點數 $f_1$、$f_2$，分別對應 Q16.16 定點數 $n_1$、$n_2$。 兩定點數可表示為 $$ \begin{align} n_1=f_1\times 2^{16} \\ n_2=f_2\times 2^{16} \end{align} $$ 兩定點數相乘後得到 $$ \begin{align} n_1\times n_2 &= (f_1\times 2^{16})(f_2\times 2^{16}) \\ &= f_1\times f_2\times 2^{32} \\ \end{align} $$ 但將兩浮點數相乘後再表示為定點數卻會得到 $$ \begin{align} n_{f_{1\times 2}} &= \left(f_1\times f_2\right)\times 2^{16} \\ &\neq n_1\times n_2 \\ \end{align} $$ 所以還需要再將 $n_1\times n_2$ 的結果除以 $2^{16}$ <small>(右移 16 位元)</small> 才會得到正確的乘法運算結果，即 $$ n_1\times n_2 = f_1\times f_2\times 2^{-16} $$ ##### `除法` 和 `乘法` 類似 兩定點數相除後得到 $$ \begin{align} n_1\div n_2 &= (f_1\times 2^{16})\div(f_2\times 2^{16}) \\ &= f_1\div f_2 \end{align} $$ 但將兩浮點數相除後再表示為定點數卻會得到 $$ \begin{align} n_{f_{1\div 2}} &= \left(f_1\div f_2\right)\times 2^{16} \\ &\neq n_1\times n_2 \end{align} $$ 所以還需要再將 $n_1\times n_2$ 的結果乘以 $2^{16}$ <small>(左移 16 位元)</small> 才會得到正確的乘法運算結果，即 $$ n_1\div n_2 = f_1\div f_2\times 2^{16} $$ <br> #### 牛頓法求根號 [參考資料](https://hackmd.io/@sysprog/sqrt#%E7%89%9B%E9%A0%93%E6%B3%95%E6%B1%82%E5%B9%B3%E6%96%B9%E6%A0%B9) 牛頓法迭代公式: $$ x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)} $$ 考慮求 $\sqrt {n}:$ 可以將原題目轉換為求 $f(x)=x^2-n$ 的根 (因為 $x=\sqrt n$ 為函式 $f(x)$ 的解)後得到: $$ \begin{align} f(x)&=x^2-n \\ f'(x)&=2x \end{align} $$ 接著將 $f(x)$ 和 $f'(x)$ 代入到迭代公式 得到: $$ \begin{align} 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}\left(x_n-\frac{n}{x_n}\right) \\ &=\frac{1}{2}\left(x_n+\frac{n}{x_n}\right) \end{align} $$ ### 情境 1 > 輸入是 IEEE float (單精度)，輸出是 int 型態 (採用 LP64 data model) > 針對上述情境，皆不可使用 FPU，只能藉由定點數予以計算，務必降低誤差並比較 glibc 的效能表現，提出改進方案 浮點數結構: $n=2^S*M*2^E=M*2^E$ :::info 由於要求根號，所以只考慮 f 為正的情況 (即 S=0) ::: (負數以例外狀況處理) 開平方後得到 $$ \begin{align} \sqrt n&=\sqrt {M*2^E} \\ &=\sqrt M*{2^\frac{E}{2}} \\ \end{align} $$ 如果 $E$ 是偶數，即 $E=2k$，則 $\sqrt n=\sqrt M*{2^k}$ 如果 $E$ 是奇數，即 $E=2k+1$，則 $\sqrt n=\sqrt M*\sqrt 2*{2^k}=\sqrt {2M}*{2^k}$ 可以發現只需要計算出 $\sqrt M$，就可以快速得到 $\sqrt n$ <br> 在計算前，我們還需要先得到浮點數的 $E$ 和 $M$ ![image](https://hackmd.io/_uploads/rJkfMBSSge.png) * `E` : 為第 23 位元到第 30 位元，所以需要右移 23 位，並且保留需要的 8 位。因為是讀出操作，所以需要再減去 127 * `M` : 為第 0 位元到第 22 位元，不需要做任何移位操作，但仍需要保留需要的 23 位 ```c E = ((n >> 23) & 0xFF) - 127; M = n & 0x7FFFFF; ``` 因為浮點數不能直接做位移運算，所以需要使用 union 將浮點數的 32 位元表示讀取到一個 uint32_t 中 ```c typedef union { float f; uint32_t i; } fi_union; fi_union u = f; ``` 使用 `u.i` 求出 $E$ 和 $M$ (需要額外補上整數部分的 `1`) ```c 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 位元為小數部分)。 ```c uint64_t M_fixed = (uint64_t) M << 9 ; ``` 前面提到， 如果 $E$ 是偶數，$\sqrt n=\sqrt M*{2^k}$ 如果 $E$ 是奇數，$\sqrt n=\sqrt {2M}*{2^k}$ 所以計算 $M$ 前還要檢查 $E$ 是否為奇數。 ```c // 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}\left(x_n+\frac{n}{x_n}\right) $$ ```c 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` ```c uint64_t result = sqrt_M_fixed << (E >> 1); ``` 因為回傳的值為 `int`，所以只需要保留前 32 位元 ```c uint64_t result = (sqrt_M_fixed << (E >> 1)) >> 32; ``` 整理後得到 ```c uint64_t sqrt_M_fixed = sqrt_Newton_64(M_fixed); uint64_t result = sqrt_M_fixed >> (32 - (E >> 1)); ``` 回傳結果並強制轉型為 `int` <br><br> ### 情境2 > 輸入是 IEEE float，輸出也是 float 型態 > 針對上述情境，皆不可使用 FPU，只能藉由定點數予以計算，務必降低誤差並比較 glibc 的效能表現，提出改進方案 執行流程和[情境1](#%E6%83%85%E5%A2%83-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$ 個位元。 ![image](https://hackmd.io/_uploads/SkswMSSBex.png) 對齊後去除多餘的整數部分 (`1`) ```c uint32_t M_new = (sqrt_M_fixed >> 9) & 0x7FFFFF; ``` * `E_new` : 前面提到，$\sqrt M$ 是一個落在 $[1.0, 2.0)$ 區間的數，以 1.OOO 來表示，再代回原算式: $$ \begin{align} \sqrt n&=\sqrt M*{2^k} \\ &=1.OOO*{2^k} \end{align} $$ 可以發現這就是一個浮點數表示法，所以浮點數的 exponent 就等於 $k$ (`E >> 1`)。但因為是寫入操作，所以要再加上 127 ```c int32_t E_new = (E >> 1) + 127; ``` 串接 `E_new` 和 `M_new` ![image](https://hackmd.io/_uploads/ryKGcEHrgl.png) ```c fi_union result; result.i = (E_new << 23) | M_new; ``` 回傳 `result.f` <br><br> ### 情境3 > 輸入和輸出都是 double (倍精度) 型態 > 針對上述情境，皆不可使用 FPU，只能藉由定點數予以計算，務必降低誤差並比較 glibc 的效能表現，提出改進方案 <br><br> ## TODO: 計算幾何平均數 ### `說明` 幾何平均數 (Geometric Mean): $$GM=\sqrt[n]{x_1\cdot x_2\cdot x_3\cdot \ldots\cdot x_n}$$ 計算的方法很簡單，但問題也很明顯。 * 如果數值都是小於 1 的小數，乘積就會趨近於 0，超出浮點數可表示的精度 * 如果數值都是很大數，乘積就會很容易超出最大可表示的範圍 要解決溢位的問題，可以將原本幾何平均數的公式取個對數 即: $$ln(GM)=ln\left(\sqrt[n]{x_1\cdot x_2\cdot x_3\cdot \ldots\cdot x_n}\right)=\frac{1}{n}\sum\limits_{i = 1}^nln\left({x_i}\right)$$ 這樣就能將原本的連乘運算轉成連加運算。 但因為要求的值是 $GM$，不是 $ln(GM)$，所以還要再做個指數運算 $$GM=exp(ln(GM))=exp\left(\frac{1}{n}\sum\limits_{i = 1}^nln\left({x_i}\right)\right)$$ 才能得出原本的 $GM$ ### 情境1 > 給定若干 int 型態的資料，在 O(1) 空間複雜度的前提，計算幾何平均，要考慮到 overflow，並探討限制和改進方案 ### 情境2 > 改為 float 型態，其餘同情境 1 > 針對上述情境，皆不可使用 FPU，只能藉由定點數予以計算，務必降低誤差並比較 glibc 的效能表現，提出改進方案 ## TODO: 以定點數實作快速傅利葉轉換 (FFT) > 務必降低誤差並持續提出改進方案 ## TODO: 探討上述在 Linux 核心的應用 > 搭配課程教材，理解開平方根和幾何平均和 FFT 如何應用於 Linux 核心，予以探討其實作手法