# 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{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$

* `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*{2^k}$
如果 $E$ 是奇數，$\sqrt n=\sqrt {2M}*{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}\left(x_n+\frac{n}{x_n}\right)
$$

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`

<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$ 個位元。
對齊後去除多餘的整數部分 (`1`)

    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

    int32_t E_new = (E >> 1) + 127;

串接 `E_new` 和 `M_new`

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 核心，予以探討其實作手法