owned this note
owned this note
Published
Linked with GitHub
# 2024q1 Homework4 (quiz3+4)
contributed by < `vax-r` >
## 第三週測驗
### 測驗一
測驗中程式碼利用 [digit-by-digit calculation](https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Digit-by-digit_calculation) 做為基礎來運算 `x` 的平方根。
基本原理很單純，透過 $(x + y)^2 = x^2 + 2xy + y^2$ 並隨著項次層層拆解，以下假設 $x = N^2$ ， $N$ 即是我們想求的答案，先透過二進位表示法將 $N$ 表示為
$$
N = \sum_{i=0}^n a_i
$$
其中 $a_i = 2^i$ 或 0 ，我們要做的就是從 $a_n$ 到 $a_0$ 逐項找出 $a_i$ 是 $2^i$ 還是 0 。首先預設初始值
$\begin{gather*}
a_n = 2^n, \ \ (2^n)^2 \le N^2
\end{gather*}$
接著我們逐步思考。
$$
N^2 = (\sum_{i=0}^n a_i)^2 = a_n^2 + 2a_n(\sum_{i=0}^{n-1}a_i) + (\sum_{i=0}^{n-1} a_i)^2 = a_0^2 + 2a_0(\sum_{i=1}^{n}a_i) + (\sum_{i=1}^{n} a_i)^2
$$
在這裡我們可以假設 $P_m = \sum_{i=m}^n a_i$ ，因此上式可以表示為
$$
N^2 = P_0^2 = a_0^2 + 2a_0P_1 + P_1^2
$$
不斷推展的話可以得到
$$
P_m^2 = a_m^2 + 2a_mP_{m+1} + P_{m+1}^2 = P_{m+1}^2 + (2P_{m+1} + a_m)a_m
$$
此處我們再假設 $Y_m = P_m^2 - P_{m+1}^2 = (2P_{m+1} + a_m)a_m$ ，又可以把上式表示為
$$
P_m^2 = P_{m+1}^2 + Y_m
$$
到這裡其實就可以利用每次利用 $P_{m+1}$ 推得 $P_m$ 然後計算 $P_m^2 \le N^2$ 是否成立，來判斷 $a_m = 2^m$ 與否，但每次都要進行 $P_m^2$ 成本太昂貴，同時我們又發現 $N^2 - P_m^2$ 的值若表示為 $X_m$ ，則可以發現 $X_m = N^2 - P_m^2$ 且 $X_{m+1} = N^2 - P_{m+1}^2$ ，因此 $X_m - X_{m+1} = -P_m^2 + P_{m+1}^2 = -Y_m$ ，所以每次想要計算的 $N^2 - P_m^2$ 也可以推得以下遞迴式
$$
X_m = X_{m+1} - Y_m
$$
透過 $X_m < 0$ 成立與否來判斷 $Y_m$ 同時也得知 $a_m$ 的值，最再把$Y_m$ 拆解為 $c_m = 2P_{m+1}a_m, d_m = (a_m)^2$
$\begin{gather*}
Y_m = \begin{cases}
c_m + d_m \ \ ,if \ \ a_m=2^m\\
0 \ \ ,if \ \ a_m=0\\
\end{cases}
\end{gather*}$
在每次迴圈更新 $c_m, d_m$ 來計算 $Y_m$ 進而推得 $a_m$
$\begin{gather*}
c_{m-1} = P_m2^m = (P_{m+1} + a_m)2^m = P_{m+1}2^m + a_m2^m \\ = \begin{cases}
\cfrac{c_m}{2} + d_m , \ \ if \ \ a_m=2^m \\
\cfrac{c_m}{2}\\
\end{cases}
d_{m-1} = \cfrac{d_m}{4}
\end{gather*}$
計算到最後得到的 $c_{-1}=P_02^0=N$ 即是答案。
函式中 `z` 即可對應到 $c_m$ ，而 `m` 即可對應 $d_m$, `b` 即是 $Y_m$ ， `x` 則是代表前一輪的差值也就是 $X_{m+1}$。
#### 使用 `ffs()`, `fls()` 進行改寫
- [ ] `ffs()`
從 lsb 開始找第一個 set bit 的位置並回傳 index ， index 由 1 開始算。利用 `ffs()` 和 shift 找出 msb 的位置並減一。
```c
int i_sqrt_ffs(int x)
{
if (x <= 1)
return x;
int tmp = x, msb = 0;
while (tmp) {
int i = ffs(tmp);
msb += i;
tmp >>= i;
}
msb--;
int z = 0;
for (int m = 1UL << (msb & ~1UL); m; m >>= 2) {
int b = z + m;
z >>= 1;
if (x >= b) {
x -= b;
z += m;
}
}
return z;
}
```
> [commit c08e25b](https://github.com/vax-r/bitwise_lab/commit/c08e25be785e6110b505ae637d351a6e824ec8d0)
#### Linux 核心 `block/blk-wbt.c`
此檔案的程式碼是在處理 bottle writeback throttling ，透過一個 time window 來監控延遲時間，如果在該 time window 中的最小延遲時間超過給定 target ，則增加 scaling step 並把 queue depth 除上某個二的倍數，且 monitoring window 會被縮到 `100 / sqrt(scaling step + 1)` 。
對應到程式碼當中即是
`rwb->cur_win_nsec = div_u64(rwb->win_nsec << 4, int_sqrt((rqd->scale_step + 1) << 8));`
但為何 100 是 `rwb->win_nsec << 4` ？ 可以理解 shift 是為了盡量用到 64 位元 ，分母 shift 8 位元是為了開平方根後就剩 shift 4 。
```c
static void rwb_arm_timer(struct rq_wb *rwb)
{
struct rq_depth *rqd = &rwb->rq_depth;
if (rqd->scale_step > 0) {
/*
* We should speed this up, using some variant of a fast
* integer inverse square root calculation. Since we only do
* this for every window expiration, it's not a huge deal,
* though.
*/
rwb->cur_win_nsec = div_u64(rwb->win_nsec << 4,
int_sqrt((rqd->scale_step + 1) << 8));
} else {
/*
* For step < 0, we don't want to increase/decrease the
* window size.
*/
rwb->cur_win_nsec = rwb->win_nsec;
}
blk_stat_activate_nsecs(rwb->cb, rwb->cur_win_nsec);
}
```
#### 更快的整數開平方根運算
:::warning
注意書名號 (即 `〈` 和 `〉`) 的使用，不是「小於」和「大於」符號。
:::
除了上述題目使用的 digit recurrence 以外還可用查表方式來計算開平方根，閱讀〈[Analysis of the lookup table size for square-rooting](https://web.ece.ucsb.edu/~parhami/pubs_folder/parh99-asilo-table-sq-root.pdf)〉以分析查表方式的大小最小需要多少。 查表方式的開平方根運算主要牽涉到利用原本的 radicand $z$ 做為 index 先在表格當中找到 $\sqrt{z}$ 的估算值，假設是 `table[z]` ，接著利用 `table[z]` 和 `talbe[z+1]` 進行線性插值法或二分搜尋。所以一開始的初始估算值非常重要，造成初始估算值的誤差兩大來源分別是
1. 只考慮 radicand $z$ 的部分位數
2. 讀取的值大小並非完整的 word
首先思考如何準確決定 $q = \sqrt{z}$ 的前 $h$ 個位元，我們希望 lookup table 找到的值介於 $(q - r^{-h}, q + r^{-h})$ ，論文中證明了如果希望找到 $q = \sqrt{z}$ 以 $r$ 為底的前 $h$ 位元，則我們需要 $z$ 以 $r$ 為底的 $2h$ 個位元，也就是至少需要讀取 $r^{2h}$ 個 words 。如果是希望獲得 $h$ digits of convergence 則至少需要 $z$ 的 $h$ 個位元，以二進位來說需要 $h \lceil lg(r) \rceil - 1$ 位元。
〈[Timing Square Root](https://web.archive.org/web/20210208132927/http://assemblyrequired.crashworks.org/timing-square-root/)〉則是量測不同開根號實作之間的效能差異，共比較了六種方法
1. 編譯器原本的 `sqrt()` 實作 (編譯為 x87 FSQRT opcode)
2. [SSE](https://web.archive.org/web/20210208013641/https://en.wikipedia.org/wiki/Streaming_SIMD_Extensions) 的 opcode `sqrtss`
3. [magic number approximation technique](https://web.archive.org/web/20210208132927/http://www.beyond3d.com/content/articles/8/)
4. 利用 SSE opcode `rsqrtss` 取得 $\cfrac{1}{\sqrt{x}}$ 的估算值並乘上 $x$
5. 利用方法 4. 並加上 [Newton-Raphson iteration](https://web.archive.org/web/20210330153224/https://en.wikipedia.org/wiki/Newtons_method) 來增加精度。
6. 方法 5. 的第 13 行展開使得每次迴圈都能處理四個浮點數。
實驗結果是方法四 SSE `rsqrtss` 乘上 `x` 的運算時間最短，但平均誤差是 0.0094% ，如果利用方法 5 或 6 會讓運算時間變長，但平均誤差卻可以降到 0% ，並且編譯器內建的 `sqrt()` 是所有方法當中最慢的。
作者提到他沒想到最快的方法居然是先計算開平方根的倒數再乘上原本的數字，利用 SSE `rsqrtss` 乘上 `x` 後再利用牛頓法增加精度，能夠擊敗原先 SSE 的 opcode `sqrtss` 。
閱讀完以上素材我有兩個實作的想法
1. 先對 `x` 進行 `int_sqrt()` 得出 $y$ ，再用查表方式取得 $\sqrt{\cfrac{x}{y^2}} = z$ ，兩者相乘 $y\cdot z$ 即可逼近 $\sqrt{x}$ 。
2. 探討如何在 C 語言使用 SSE `rsqrtss` 取得 $\cfrac{1}{\sqrt{x}}$ ，乘上 $x$ 後再使用牛頓法逼近精度。
實作 fast inverse square root 演算法需要用到一個 magic number `0x5f3759df` ，對於 IEEE 754 規格的浮點數才能直接使用，轉換成二進位會變成 `01011111001101110101100111011111` ，如果採用定點數則 fractional bits 只能用 1 bit ，我認為不是把 magic number 直接轉成二進位，應探討為何 IEEE 754 規格可以推出這個浮點數，我的定點數是否能導出類似的數。
利用 $lg(\sqrt{x}) = \cfrac{1}{2} lg(x)$ 的概念，搭配 [Newton's method](https://en.wikipedia.org/wiki/Newton%27s_method) 進行逼近，實作如下
```c
int ilog2(unsigned n)
{
return (31 - __builtin_clz(n | 1));
}
unsigned i_fastsqrt(unsigned x)
{
int y = (ilog2(x) >> 1);
unsigned z = 1 << y;
z = (z + (x / z)) >> 1;
z = (z + (x / z)) >> 1;
z = (z + (x / z)) >> 1;
return z;
}
```
從 1 到 1000000 和題目的 `i_sqrt()` 整數開平方根做比較，可以發現非常接近，做圖結果如下
利用以下腳本計算兩種不同方法的平均執行時間
```c
#define bench(statement) \
({ \
struct timespec _tt1, _tt2; \
clock_gettime(CLOCK_MONOTONIC, &_tt1); \
statement; \
clock_gettime(CLOCK_MONOTONIC, &_tt2); \
long long time = (long long) (_tt2.tv_sec * 1e9 + _tt2.tv_nsec) - \
(long long) (_tt1.tv_sec * 1e9 + _tt1.tv_nsec); \
time; \
})
int main(void) {
long long isqrt_time = 0;
long long fastsqrt_time = 0;
for (int i = 1; i < 1000000; i++) {
isqrt_time += bench(i_sqrt(i));
fastsqrt_time += bench(i_fastsqrt(i));
}
printf("i_sqrt avg time : %lf\n", (double) isqrt_time / 1000000);
printf("fast sqrt avg time : %lf\n", (double) fastsqrt_time / 1000000);
return 0;
}
```
結果如下
```shell
$ ./out
i_sqrt avg time : 55.248122
fast sqrt avg time : 51.479055
```
如果利用 `$ sudo perf stat --repeat 100 ./out` 分別量測兩種不同實作的 cpu cycle 與 cpu clock utilized 分別可以得到以下結果
* `i_sqrt()`
```
Performance counter stats for './out' (100 runs):
79.11 msec task-clock # 0.996 CPUs utilized ( +- 0.26% )
0 context-switches # 0.000 /sec
0 cpu-migrations # 0.000 /sec
52 page-faults # 657.317 /sec ( +- 0.13% )
2,7698,5996 cycles # 3.501 GHz ( +- 0.09% ) (83.56%)
1,5255,7505 stalled-cycles-frontend # 55.08% frontend cycles idle ( +- 0.09% ) (83.88%)
5742,1358 stalled-cycles-backend # 20.73% backend cycles idle ( +- 0.20% ) (67.79%)
3,4771,1663 指令 # 1.26 insn per cycle
# 0.44 stalled cycles per insn ( +- 0.02% ) (83.91%)
5555,5190 branches # 702.257 M/sec ( +- 0.03% ) (83.88%)
6,2454 branch-misses # 0.11% of all branches ( +- 16.45% ) (80.90%)
0.079415 +- 0.000209 seconds time elapsed ( +- 0.26% )
```
* `i_fastsqrt()`
```
Performance counter stats for './out' (100 runs):
74.05 msec task-clock # 0.996 CPUs utilized ( +- 0.31% )
0 context-switches # 0.000 /sec
0 cpu-migrations # 0.000 /sec
51 page-faults # 688.735 /sec ( +- 0.13% )
2,5955,4261 cycles # 3.505 GHz ( +- 0.06% ) (83.70%)
1,5273,2557 stalled-cycles-frontend # 58.84% frontend cycles idle ( +- 0.06% ) (83.69%)
9279,8175 stalled-cycles-backend # 35.75% backend cycles idle ( +- 0.12% ) (67.39%)
2,5141,9842 指令 # 0.97 insn per cycle
# 0.61 stalled cycles per insn ( +- 0.03% ) (83.69%)
3513,1310 branches # 474.434 M/sec ( +- 0.02% ) (83.69%)
2264 branch-misses # 0.01% of all branches ( +- 6.87% ) (81.53%)
0.074353 +- 0.000227 seconds time elapsed ( +- 0.31% )
```
觀察結果可以發現， `i_fastsqrt()` 用了 259554261 個 cycles ， cpu clock time 是 74.05 msec ，而 `i_sqrt()` 則用了 276985996 個 cycles ， cpu clock time 是 79.11 msec ，是 `fast_sqrt()` 更快，比較指令數目的話， `i_fastsqrt()` 只使用了 251419842 道指令 ，明顯少於 `i_sqrt()` 使用的 347711663 道指令 ，少了幾乎 100000000 道指令。不管從平均執行時間，或者是 cpu cycles 數量甚至到指令數量，我實作的 `i_fastsqrt()` 都少於 `i_sqrt()` 的實作，證明我的實作更快速，而且我的實作還是 branchless 的實作，可以觀察使用到的 branches 數量跟 branch-misses 都是遠小於 `i_sqrt()` 。
:::info
我希望嘗試貢獻這段程式碼實作到 Linux 核心，但不知道該如何將實驗結果表示在 commit message 當中比較妥當。
:::
:::info
做圖到 100000000 的結果，依然幾乎貼合
:::
直接與 Linux 核心的 `int_sqrt` 比較，測試程式碼如下
```c
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
unsigned long int_sqrt(unsigned long x)
{
unsigned long b, m, y = 0;
if (x <= 1)
return x;
m = 1UL << ((31 - __builtin_clz(x)) & ~1UL);
while (m != 0) {
b = y + m;
y >>= 1;
if (x >= b) {
x -= b;
y += m;
}
m >>= 2;
}
return y;
}
int main() {
for (int i = 1; i < 1000000; i++)
int_sqrt(i);
return 0;
}
```
利用 `perf stat` 測試，重複次數為 100
```shell
$ sudo perf stat --repeat 100 ./out
Performance counter stats for './out' (100 runs):
35.37 msec task-clock # 0.992 CPUs utilized ( +- 0.62% )
0 context-switches # 0.000 /sec
0 cpu-migrations # 0.000 /sec
49 page-faults # 1.385 K/sec ( +- 0.14% )
1,2181,3348 cycles # 3.444 GHz ( +- 0.08% ) (76.56%)
6997,1269 stalled-cycles-frontend # 57.44% frontend cycles idle ( +- 0.09% ) (76.69%)
559,1291 stalled-cycles-backend # 4.59% backend cycles idle ( +- 0.20% ) (76.56%)
1,6095,3665 指令 # 1.32 insn per cycle
# 0.43 stalled cycles per insn ( +- 0.06% ) (88.28%)
2551,2990 branches # 721.380 M/sec ( +- 0.04% ) (88.28%)
1,0616 branch-misses # 0.04% of all branches ( +- 2.25% ) (81.91%)
0.035667 +- 0.000220 seconds time elapsed ( +- 0.62% )
```
然後測試我寫的版本
```c
unsigned long i_fastsqrt(unsigned long x)
{
unsigned long y = ((31 - __builtin_clzll(x | 1)) >> 1);
unsigned long z = ((1 << y) + (x >> y)) >> 1;
z = (z + (x / z)) >> 1;
z = (z + (x / z)) >> 1;
return z;
}
int main() {
for (int i = 1; i < 1000000; i++)
i_fastsqrt(i);
return 0;
}
```
再用 `perf stat` 重複測試 100 次
```shell
$ sudo perf stat --repeat 100 ./out
Performance counter stats for './out' (100 runs):
33.96 msec task-clock # 0.991 CPUs utilized ( +- 0.67% )
0 context-switches # 0.000 /sec
0 cpu-migrations # 0.000 /sec
49 page-faults # 1.443 K/sec ( +- 0.14% )
1,1645,7487 cycles # 3.430 GHz ( +- 0.10% ) (78.43%)
7323,1953 stalled-cycles-frontend # 62.88% frontend cycles idle ( +- 0.10% ) (78.40%)
2464,6153 stalled-cycles-backend # 21.16% backend cycles idle ( +- 0.14% ) (69.80%)
5621,0086 指令 # 0.48 insn per cycle
# 1.30 stalled cycles per insn ( +- 0.05% ) (89.21%)
321,0409 branches # 94.543 M/sec ( +- 0.03% ) (89.20%)
2407 branch-misses # 0.07% of all branches ( +- 7.27% ) (84.16%)
0.034267 +- 0.000229 seconds time elapsed ( +- 0.67% )
```
可以發現不管是 task-clock 、 cycles 、指令、 branches 與 branch-misses 都是我撰寫的版本優於 `int_sqrt()` 。
改寫 `x / z` 的操作我參閱了數種方法、不同的逼近方式等尚未找到能取代 division 操作的演算法，由於使用牛頓法至少必須逼近三次，每次 `z` 未必會是 2 的指數次方因此無法直接用 shift 替代。
提交 patch 後發現我的實作 precision 上也有問題，原本的 `int_sqrt()` 和我實作的 `i_fastsqrt()` 利用 gnuplot 做圖看似幾乎貼合，但實際上如果逐個答案和 `sqrt()` 的整數部分進行誤差比較，我的實作誤差比 `int_sqrt()` 大很多，並且我的實作有用到除法這對於某些沒有除法指令的架構是不可行的，應該重新思考一個方法。
### 測驗二
不運用 `/` 和 `%` 來求得除以 10 的商與餘數的方法，參考我在第二週做的筆記 [Integer Division by Constant](https://hackmd.io/@vax-r/linux2024-homework2#bitwise-operation-%E8%A3%9C%E5%BC%B7) 和 [Hacker's Delight](https://web.archive.org/web/20180517023231/http://www.hackersdelight.org/divcMore.pdf) ，除以 10 相同於乘以 $\cfrac{1}{10}$ ，因此我們要找到 $\cfrac{1}{10}$ 的二進位表示法。 由於 10 不符合 $2^k \pm 1$ 的形式因此無法二進位或 digit recurssion 來表示，我們只能盡量以一個非常接近 $\cfrac{1}{10}$ 的數字 $\cfrac{a}{2^N}$ 來表示 10 的倒數，題目挑選 $2^N = 128, a = 13, \cfrac{128}{13} \approx 9.84$ 來使用，因此 $\cfrac{1}{10} \approx \cfrac{13}{128}$ ，計算 `x/10` 相當於 `x * 1/10` 且近似於 `x * 13 / 128` ，可以進一步拆解成 `x * 13 / 8 / 16`。
`x * 13 / 8` 則可以透過 `(x >> 3) + (x >> 1) + x` 得到，而進行 right shift 的時候會遺失最低 3 個位元，所以預先保存最低 3 個位元最後才加回去，得到 `x * 13` 後再除以 128 ，如下
```c
unsigned d0 = x & 0b1;
unsigned d1 = x & 0b11;
unsigned d2 = x & 0b111;
unsigned q = ((x >> 3) + (x >> 1) + x) << 3;
q = q + d0 + d1 + d2; /* restore last 3 bits */
q >>= 7; /* divide by 128 */
```
最後餘數再由餘數定理 $r = f - g\cdot Q$ 而得。
```c
unsigned r = x - (((q << 2) + q) << 1);
```
完整函式如下
```c
void divmod_10(uint32_t in, uint32_t *div, uint32_t *mod)
{
unsigned d0 = x & 0b1;
unsigned d1 = x & 0b11;
unsigned d2 = x & 0b111;
*div = ((((in >> 3) + (in >> 1) + in) << 3) + d0 + d1 + d2) >> 7;
*mod = in - (((*div << 2) + *div) << 1);
}
```
最後包裝成題目的程式碼如下
```c
void divmod_10(uint32_t in, uint32_t *div, uint32_t *mod)
{
uint32_t x = (in | 1) - (in >> 2); /* div = in/10 ==> div = 0.75*in/8 */
uint32_t q = (x >> 4) + x;
x = q;
q = (q >> 8) + x;
q = (q >> 8) + x;
q = (q >> 8) + x;
q = (q >> 8) + x;
*div = (q >> 3);
*mod = in - ((q & ~0x7) + (*div << 1));
}
```
至於是怎麼包裝的呢？原本的程式碼實作跟最後題目給出的函式實作看似毫無關聯，事實上關聯性確實很有限，一開始的 `x = (in | 1) - (in >> 2)` 是在計算 `x * 3/4` ，並且若 `x` 為偶數則加一（不明白原因，尚需探討），接著 `q = (x >> 4) + x` 相同於計算 $\cfrac{3}{4} \cdot in + \cfrac{3}{64} \cdot in = \cfrac{102}{128} \cdot in$ 。 $\cfrac{102}{128} \approx 0.79$ 接近 0.8 也就是 $\cfrac{8}{10}$ ，為了縮小誤差連續進行四次 `q = (q >> 8) + x` 來進行逼近。
最後 `q >> 3` 相當於 $\cfrac{8}{10} \cdot in / 8 = \cfrac{in}{10}$ ，而 `q & ~0x7` 是在保留 `q` 除了最後三個 bits 以外的位元值，用意和 `*div << 3` 相同，因此 `(q & ~0x7) + (*div << 1)` 就是在計算 `*div * 10` 。
#### 比較依賴除法與否的指令集
透過 `lscpu` 命令確認運行電腦的處理器版本，我所使用的是 Intel i7 在 [Instruction tables](https://www.agner.org/optimize/instruction_tables.pdf) 當中對應到第 206 頁的 Intel Nehalem 。 首先實作使用 `/, %` 的實作
```c
static void divmod10(int in, int *div, int *mod)
{
*div = in / 10;
*mod = tmp % 10;
}
```
透過以下命令編譯並透過 `objdump` 觀察對應 assembly instruction (僅列出 `divmod10()` 函式對應指令)
```shell
$ gcc -O0 -o divmod divmod.c -lm -lc
$ objdump -Slz divmod
...
0000000000001149 <divmod10>:
divmod10():
1149: f3 0f 1e fa endbr64
114d: 55 push %rbp
114e: 48 89 e5 mov %rsp,%rbp
1151: 89 7d fc mov %edi,-0x4(%rbp)
1154: 48 89 75 f0 mov %rsi,-0x10(%rbp)
1158: 48 89 55 e8 mov %rdx,-0x18(%rbp)
115c: 8b 45 fc mov -0x4(%rbp),%eax
115f: 48 63 d0 movslq %eax,%rdx
1162: 48 69 d2 67 66 66 66 imul $0x66666667,%rdx,%rdx
1169: 48 c1 ea 20 shr $0x20,%rdx
116d: c1 fa 02 sar $0x2,%edx
1170: c1 f8 1f sar $0x1f,%eax
1173: 29 c2 sub %eax,%edx
1175: 48 8b 45 f0 mov -0x10(%rbp),%rax
1179: 89 10 mov %edx,(%rax)
117b: 8b 4d fc mov -0x4(%rbp),%ecx
117e: 48 63 c1 movslq %ecx,%rax
1181: 48 69 c0 67 66 66 66 imul $0x66666667,%rax,%rax
1188: 48 c1 e8 20 shr $0x20,%rax
118c: c1 f8 02 sar $0x2,%eax
118f: 89 ce mov %ecx,%esi
1191: c1 fe 1f sar $0x1f,%esi
1194: 29 f0 sub %esi,%eax
1196: 89 c2 mov %eax,%edx
1198: 89 d0 mov %edx,%eax
119a: c1 e0 02 shl $0x2,%eax
119d: 01 d0 add %edx,%eax
119f: 01 c0 add %eax,%eax
11a1: 29 c1 sub %eax,%ecx
11a3: 89 ca mov %ecx,%edx
11a5: 48 8b 45 e8 mov -0x18(%rbp),%rax
11a9: 89 10 mov %edx,(%rax)
11ab: 90 nop
11ac: 5d pop %rbp
11ad: c3 ret
...
```
總共有 14 個 `mov` (register to register) 共用了 14 個 cpu cycles ， `push` (register) 一個用了 3 個 cpu cycles ，兩個 `imul` 用了 6 個 cpu cycles ，兩個 `shr` 使用 2 個 cpu cycles ， 4 個 `sar` 用了 4 個 cycles ， 三個 `sub` 共用 3 個 cycles ， 一個 `shl` 用了 1 個 cycles ，總共用了 33 個 cycles 。
改為題目的實作並以相同方法測試
```shell
...
0000000000001149 <divmod10>:
divmod10():
1149: f3 0f 1e fa endbr64
114d: 55 push %rbp
114e: 48 89 e5 mov %rsp,%rbp
1151: 89 7d ec mov %edi,-0x14(%rbp)
1154: 48 89 75 e0 mov %rsi,-0x20(%rbp)
1158: 48 89 55 d8 mov %rdx,-0x28(%rbp)
115c: 8b 45 ec mov -0x14(%rbp),%eax
115f: 83 c8 01 or $0x1,%eax
1162: 89 c2 mov %eax,%edx
1164: 8b 45 ec mov -0x14(%rbp),%eax
1167: c1 e8 02 shr $0x2,%eax
116a: 89 c1 mov %eax,%ecx
116c: 89 d0 mov %edx,%eax
116e: 29 c8 sub %ecx,%eax
1170: 89 45 f8 mov %eax,-0x8(%rbp)
1173: 8b 45 f8 mov -0x8(%rbp),%eax
1176: c1 e8 04 shr $0x4,%eax
1179: 89 c2 mov %eax,%edx
117b: 8b 45 f8 mov -0x8(%rbp),%eax
117e: 01 d0 add %edx,%eax
1180: 89 45 fc mov %eax,-0x4(%rbp)
1183: 8b 45 fc mov -0x4(%rbp),%eax
1186: 89 45 f8 mov %eax,-0x8(%rbp)
1189: 8b 45 fc mov -0x4(%rbp),%eax
118c: c1 e8 08 shr $0x8,%eax
118f: 89 c2 mov %eax,%edx
1191: 8b 45 f8 mov -0x8(%rbp),%eax
1194: 01 d0 add %edx,%eax
1196: 89 45 fc mov %eax,-0x4(%rbp)
1199: 8b 45 fc mov -0x4(%rbp),%eax
119c: c1 e8 08 shr $0x8,%eax
119f: 89 c2 mov %eax,%edx
11a1: 8b 45 f8 mov -0x8(%rbp),%eax
11a4: 01 d0 add %edx,%eax
11a6: 89 45 fc mov %eax,-0x4(%rbp)
11a9: 8b 45 fc mov -0x4(%rbp),%eax
11ac: c1 e8 08 shr $0x8,%eax
11af: 89 c2 mov %eax,%edx
11b1: 8b 45 f8 mov -0x8(%rbp),%eax
11b4: 01 d0 add %edx,%eax
11b6: 89 45 fc mov %eax,-0x4(%rbp)
11b9: 8b 45 fc mov -0x4(%rbp),%eax
11bc: c1 e8 08 shr $0x8,%eax
11bf: 89 c2 mov %eax,%edx
11c1: 8b 45 f8 mov -0x8(%rbp),%eax
11c4: 01 d0 add %edx,%eax
11c6: 89 45 fc mov %eax,-0x4(%rbp)
11c9: 8b 45 fc mov -0x4(%rbp),%eax
11cc: c1 e8 03 shr $0x3,%eax
11cf: 89 c2 mov %eax,%edx
11d1: 48 8b 45 e0 mov -0x20(%rbp),%rax
11d5: 89 10 mov %edx,(%rax)
11d7: 8b 45 fc mov -0x4(%rbp),%eax
11da: 83 e0 f8 and $0xfffffff8,%eax
11dd: 89 c2 mov %eax,%edx
11df: 48 8b 45 e0 mov -0x20(%rbp),%rax
11e3: 8b 00 mov (%rax),%eax
11e5: 01 c0 add %eax,%eax
11e7: 8d 0c 02 lea (%rdx,%rax,1),%ecx
11ea: 8b 45 ec mov -0x14(%rbp),%eax
11ed: 29 c8 sub %ecx,%eax
11ef: 89 c2 mov %eax,%edx
11f1: 48 8b 45 d8 mov -0x28(%rbp),%rax
11f5: 89 10 mov %edx,(%rax)
11f7: 90 nop
11f8: 5d pop %rbp
11f9: c3 ret
...
```
共使用了 44 個 `mov` 與其他若干指令，佔用的 cpu cycles 明顯多於原本使用 `/` 和 `%` 的版本。
用 `perf stat` 分析，使用 `/, %` 的實作共用了 122,3235 cpu cycles ，不使用的版本用了 122,3024 個 cpu cycles ，沒有明顯差異。
> CPU 週期數量「沒有明顯差異」，就說明此舉的效益，當指令集受限時，勢必要調整實作手法，而且看似指令數量增加，但避開較長週期的特定指令，有機會在亂序 (out-of-order) 執行的處理器上，爭取到更高的 IPC。 :notes: jserv
此處讓我很好奇，如果使用使用 `/, %` 的實作反而更好或沒有明顯差異，為何我們要採取 shifting 的策略而避免使用 `/, %` ？編譯器是否對 `/, %` 等指令做了優化的行為？又或者硬體架構上除法指令集本來就設計為使用 shift ，而我撰寫的程式反而讓它多了幾道指令？尚需了解編譯器優化行為還有檢討我的實驗方法。
> 考慮到 RISC-V 一類的處理器，其主體的 ISA (如 RV32I 和 RV64) 甚至沒有乘法指令，上述的調整就是必然。 :notes: jserv
:::info
* i : immediate data
* r : register
* mm : 64 bit mmx register
* xmm : 128 bit xmm register
* mm/x : mmx or xmm register
* sr : segment register
* m : memory
* m32 : 32-bit memory operand
:::
#### modulo 9 and modulo 5 without `/, %`
> [commit cfb7b8e](https://github.com/vax-r/bitwise_lab/commit/cfb7b8e8461e13a09bf75ce95727bef30903010a)
參照自 [Hacker's Delight](https://web.archive.org/web/20180517023231/http://www.hackersdelight.org/divcMore.pdf) 。
### 測驗三
#### 原理解釋
原本的方法是利用二分搜尋，輸入是 32 位元長的資料型態最大可能到 $2^{32} - 1$ ，先判斷 `i` 是否大於 $2^{16}$ ，是的話 $lg(i) \ge 16$ ，接著左移 16 位元在判斷是否大於 $2^{8}$ ，實際上是判斷 `i` 是否大於 $2^{24}$ ，是的話則把答案再加 8 bits 。
可以發現 `ilog2(n)` 就是在計算 `n` 從 msb 數來第一個 set bit 的位置 `__builtin_clz(n)` 先計算 `n` 的 leading zeroes 接著用 31 減掉它也就是從 msb 數來第一個 set bit 的位置。
#### Linux 核心應用案例
在 [include/linux/log2.h](https://elixir.bootlin.com/linux/latest/source/include/linux/log2.h) 可以找到數個和題目程式碼類似的運算，包含 `__ilog2_u32(u32 n), const_ilog2(n), ilog2(n)` 等巨集與函式。
### 測驗四
[Moving average](https://en.wikipedia.org/wiki/Moving_average#Simple_moving_average) 是對一系列資料，有數種變形。最簡單的便是 [simple moving average](https://en.wikipedia.org/wiki/Moving_average#Simple_moving_average) ，把指定時間區間內的資料加總並取平均如下
$$
SMA_k = \cfrac{1}{k}(\sum_{i=n-k+1}^{n} p_i)
$$
如果希望減少時間較久以前的資料，可以利用 [Exponential smoothing](https://en.wikipedia.org/wiki/Exponential_smoothing) 也就是題目中提到的方法，他利用 exponential window function 來對時間序列的資料做 smoothing 。
基本形式如下，假設時間 t 的資料是 $x_t$ ，則經過 exponential smoothing 後的輸出是 $s_t$，初始資料 $s_0 = x_0$，公式如下
$$
s_t = \alpha x_t + (1 - \alpha)s_{t-1}, \ \ t>0
$$
其中 $\alpha$ 稱為 smoothing factor 。
對應到題目的程式碼， `avg` 是用來存取平均數的結構體， `internal` 是目前的 ewma ，值得注意的是 `internal` 是 `unsigned long` 的資料型態，但平均數應該會是浮點數，因此我們可以確認此程式碼利用定點數運算， `factor` 正是這裡定點數設計的 scaling factor ，我們可以從讀取也就是 `ewma_read` 來確認，每次要讀取 `avg->internal` 也就是目前資料的 ewma 值時都會先右移 `avg->factor` 個位元，因此讀取到的是 ewma 的整數部分。
另外寫入時使用 `ewma_add` ，若目前 `avg` 結構當中還沒有資料，則將新加入的 `val` 左移 `avg->factor` 個位元再存入 `avg->internal` 中，和 $s_0 = x_0$ 相同。特別的是 `avg->internal` 不為 0 時，更新 `avg->internal` 的方式居然是 `(((avg->internal << avg->weight) - avg->internal) + (val << avg->factor)) >> avg->weight` ，如果把 `avg->weight` 理解為 $\alpha$ ，不應該是 `avg->internal - avg->weight * avg->internal` 這樣才符合 $(1-\alpha)x_{t-1}$ 嗎？
原來此結構中的 $\alpha$ 其實是 $\cfrac{1}{2^{avg->weight}}$ ， `ewma_add` 當中的運算是先將 $s_t = \alpha x_t + (1 - \alpha) x_{t-1}$ 化為 $2^{avg->weight} \cdot s_t = x_t + (2^{avg->weight} - 1) s_{t-1}$ ，算完之後再除上 $2^{avg->weight}$ 也就是程式當中最後右移 `avg->weight` 個位元的原因。
#### Linux 核心相關應用案例
在 [include/linux/average.h](https://elixir.bootlin.com/linux/latest/source/include/linux/average.h) 定義了 `DECLARE_EWMA(name, _precision, _weight_rcp)` 此巨集，會依據給定的 `name` 給開發者一個 `struct ewma_name` 結構體和 `init`, `read`, `add` 三個函式。此處利用定點數運算 fractional bits 有 `_precision` 個， 而 ewma 中的 $\alpha$ 則由 `1/_weight_rcp` 表示， `_weight_rcp` 必須是二的指數次方。此巨集當中的函式實作方式和我們題目當中非常類似，不過此巨集的結構體只儲存 ewma 值， smoothing factor 會在每次呼叫 `ewma_name_add` 時利用 `ilog2()` 重新計算，我們可以評估是否在 `ewma_name_init` 處就把 `weight_rcp` 當成一個成員儲存起來，如此一來就省去每次呼叫 `ilog2()` 的成本。
而 [drivers/net/wireless/ralink/rt2x00/rt2x00.h](https://elixir.bootlin.com/linux/latest/source/drivers/net/wireless/ralink/rt2x00/rt2x00.h) 當中則是透過 `DECLARE_EWMA(rssi, 10, 8)` 定義一個結構體 `struct ewma_rssi`，定點數的 fractional bits 是 10，smoothing factor 則是 $\cfrac{1}{8}$ ，在 `struct link_ant` 結構體中的成員 `struct ewma_rssi rssi_ant` 就利用 ewma 來算 rssi 的 walking average。
[drivers/net/wireless/ralink/rt2x00/rt2x00.c](https://elixir.bootlin.com/linux/latest/source/drivers/net/wireless/ralink/rt2x00/rt2x00link.c) 當中定義的函式 `rt2x00link_get_avg_rssi` 就是 `ewma_rssi_read` 的包裝。
rt2x00 核心模組當中多處皆運用到 ewma 的運算來更新無線網路的 RSSI。
> [Received signal strength indicator](https://en.wikipedia.org/wiki/Received_signal_strength_indicator) (RSSI)
### 測驗五
```c
int ceil_ilog2(uint32_t x)
{
uint32_t r, shift;
x--;
r = (x > 0xFFFF) << 4;
x >>= r;
shift = (x > 0xFF) << 3;
x >>= shift;
r |= shift;
shift = (x > 0xF) << 2;
x >>= shift;
r |= shift;
shift = (x > 0x3) << 1;
x >>= shift;
return (r | shift | x > 1) + 1;
}
```
`r` 用來記錄 `x` 從 msb 數來第一個 set bit 由右數來的位置，和 `ilog2()` 記錄一樣的數值，但 `ilog2` 是計算 floor of log2 ，同時因為 `r` 和 `shift` 一定是 2 的指數因此加法可以利用 bitwise or `|` 進行 ， `r | shift` 等效於 `r + shift` ，最後還要判斷 `x > 1` 是因為判斷 `x > 0x3` 時若 `x == 0x2, x == 0x3` 則 shift 依舊會是 0 ，會少算 1 ，因此才要多判斷 `x > 1` 。
一開始進行 `x--` 是為了處理 `x` 剛好為 2 的指數次方時的狀況，如果 `x` 剛好為 2 的指數，沒有在一開始減一的話最後輸出會因為最後的 `+1` 而比正確答案多一。
#### 改進程式碼
由於一開始 `x--` ，使得 `x=0` 輸入函式的話經過 `x--` 後會變成 `0xFFFFFFFF` ，最簡單的改進便是在將 `x--` 變為 `x -= !!x` ，也就是在 `x > 0` 才進行減法。
```diff
int ceil_ilog2(uint32_t x)
{
uint32_t r, shift;
- x--;
+ x -= !!x;
```
如此一來便能在 branchless 情況下處理 `x = 0` 的狀況。
> [commit ae309e9](https://github.com/vax-r/bitwise_lab/commit/ae309e920deef7e5156f5cd0167f8445644d8e73)
#### Linux 核心應用案例
TODO
---
## 第四週測驗
### 測驗一
```c
int totalHammingDistance(int* nums, int numsSize)
{
int total = 0;;
for (int i = 0;i < numsSize;i++)
for (int j = 0; j < numsSize;j++)
total += __builtin_popcount(nums[i] ^ nums[j]);
return total >> 1;
}
```
#### 原理解釋
閱讀 [Population Count](https://www.chessprogramming.org/Population_Count) 和 [Hacker's Delight](https://doc.lagout.org/security/Hackers%20Delight.pdf) 。
首先假設以二進位表示的話數字 $x=b_{31}\cdot b_2b_1b_0$ ，以數學式來定義 population count 如下
$$
popcount(x) = \sum_{n=0}^{31} b_n = \sum_{n=0}^{31} (2^n - \sum_{i=0}^{n-1} 2^i) b_n
$$
重新排列後可以推得以下公式
$$
popcount(x) = x - \lfloor \cfrac{x}{2} \rfloor - \lfloor \cfrac{x}{4} \rfloor - \cdots - \lfloor \cfrac{x}{2^{31}} \rfloor
$$
題目的程式實作利用 $x - \lfloor \cfrac{x}{2} \rfloor - \lfloor \cfrac{x}{4} \rfloor - \lfloor \cfrac{x}{8} \rfloor$ 每四個位元為單位進行計算。
`n = (v >> 1) & 0x77777777` 即是在計算 $\lfloor \cfrac{v}{2} \rfloor$
因此程式碼前六行即是在計算 $x - \lfloor \cfrac{x}{2} \rfloor - \lfloor \cfrac{x}{4} \rfloor - \lfloor \cfrac{x}{8} \rfloor$
此時 `v` 當中每四個位元儲存的數字大小即是原本 `v` 當中每四個位元的 population count ，目前的 `v` 可以表達為 $v = \sum_{i=0}^{7} B_i * 16^{i}$ ， $B_i$ 原本 `v` 當中每四個 bits 的 population count 。 `v = (v + (v >> 4)) & 0x0F0F0F0F` 即是將 $B_i, B_{i-1}$ 相加所以可得
```
0 (B7+B6) 0 (B5+B4) 0 (B3+B2) 0 (B1+B0)
```
最後 `v *= 0x01010101` 即是在第四個區塊做 $\sum_{i=0}^{7} B_i$ ，然後結果會為在第四個區塊中所以右移 24 位元 `v >> 24` 。
最後是 [Hamming Distance](https://en.wikipedia.org/wiki/Hamming_distance) 的定義與運用 population count 計算 Hamming Distance 的方法。
在針對 [Leetcode 477. Total Hamming Distance](https://leetcode.com/problems/total-hamming-distance/description/) 的程式碼當中，因為給定程式碼的做法會將每兩個數字的 hamming distance 重複算兩遍，因此加總的最後要除以 2 ，所以向右位移 1 位元即可。
#### 程式碼改進
如果直接把題目程式碼丟進 leetcode ，結果會是 TLE ，要將迴圈順序稍微改變如下
```diff
int totalHammingDistance(int* nums, int numsSize) {
int total = 0;;
for (int i = 0;i < numsSize;i++)
- for (int j = 0; j < numsSize;j++)
+ for (int j = i+1; j < numsSize;j++)
total += popcount_branchless(nums[i] ^ nums[j]);
return total;
}
```
我透過每兩位元 (nibble) 相加可得該二位元的 population count 原理重複計算，程式如下
```c
unsigned popcount_v2(unsigned v)
{
v = v - ((v >> 1) & 0x55555555);
v = (v & 0x33333333) + ((v >> 2) & 0x33333333);
v = (v + (v >> 4)) & 0x0f0f0f0f;
v = (v + (v >> 8)) & 0x00ff00ff;
v = (v + (v >> 16)) & 0x0000ffff;
return v;
}
```
透過 `perf stat --repeat 100` 重複量測 100 次，和題目原先提供的 `popcount_branchless` 相比， task-clock 從 487.32 msec 降至 468.83 msec ， 耗費的 cycles 數量也從 18,2015,8606 降至 17,5005,2307 。
:::info
在 solution 區域看到別人的解法，速度快過 100% submissions ，但原理我尚未理解
```c
int totalHammingDistance(int* nums, int numsSize){
int total_dist = 0;
int bits = 0;
for (int i = 0; i < 32; i++) {
for (int j = 0; j < numsSize; j++) {
bits += (nums[j] >> i) &1;
}
total_dist += bits*(numsSize-bits);
bits = 0;
}
return total_dist;
}
```
:::
### 測驗二
#### 原理解釋
模擬一百萬次棋局，棋盤大小是 3x3 所以最多有 9 種可能的選擇，相較於把下棋看成在九宮格中選擇一格並標記，可以把棋局看成是 8 種可能的路線（三條橫線、三條縱線、兩條對角線），分別對應到 `uint32_t` 以 16 進位表示的不同位元，第一條橫線為 msb ，第二條橫線為 msb - 1 ，依此類推。若滿足其中一個路線大小為 7 則獲勝，因為化為二進位來看 0x7 即是 0111 ，加 1 後會變成 8 ，因此判斷獲勝與否的函式即是
```c
static inline uint32_t is_win(uint32_t player_board)
{
return (player_board + 0x11111111) & 0x88888888;
}
```
每一個下棋的選擇都會影響多條路線的狀態，把選擇第 0 格到第 8 格分別會影響的狀態紀錄在 `move_masks` 當中，例如選擇第 0 格也就是棋盤最左上方的格子，對應到更新的路線會是第一條橫線、第一條縱線、第一條對角線，而第 0 格同時是這些路線最左邊的格子，因此會是 0100 ，對應到該三條路線狀態為 `0x400400400` ，其他元素也是相同概念。
每一局棋局當中一開始的 `available moves` 都是九個，透過 `xorshift()` 產生隨機亂數後再除以 `n_moves` 以隨機找一個選項下棋，其中找餘數的函式是利用 `fastmod()` 。函式當中若除以 2 的指數例如 2、4 或 8 ，可以直接和該數減一進行 bitwise AND ，除以三和七則利用函式 `mod3(), mod7()` 來處理。
觀察 `mod7()` 函式，其中 `0x24924925` 是 $\lceil \cfrac{2^{32}}{7} \rceil$ ，從 [Hacker's Delight](https://doc.lagout.org/security/Hackers%20Delight.pdf) 我們知道利用 $8^k \equiv 1 (mod \ \ 7)$ 此數學關係式，可以推導出 $n \ \ (mod \ \ 7) = popcount(n)$ ，例如書中有給一段使用 lookup table 的程式碼
```c
int remu7(unsigned n) {
static char table [75] = {0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4};
n = (n >> 15) + (n & 0x7FFF);
n = (n >> 9) + (n & 0x001FF);
n = (n >> 6) + (n & 0x0003F);
return table[n];
}
```
原理是利用 $n \ \ (mod \ \ 7) = popcount(n)$ ，不斷的進行 digit sum ，其實每一步都在重複求 modulo 7 ，而我們題目中的程式另外運用了 $n \ \ (mod \ \ 7) \equiv \lfloor \cfrac{8}{7}n \rfloor \ \ (mod \ \ 8)$ 這個關係式， 為了求得 $\lfloor \cfrac{8}{7}n \rfloor \ \ (mod \ \ 8)$ 我們可以計算 $\lfloor \cfrac{2^{32}}{7} \rfloor$ 再右移 29 位元 ，由於取 floor 會有誤差 ，題目的程式碼取 ceil 也就是 $\lceil \cfrac{2^{32}}{7} \rceil$ 結果就如同程式碼中的 `0x24924925` ，乘上這個數後許多 high order bits 會被遺棄，最後右移 29 位元得到的數字即是想要的答案。
### 測驗三
* `treeint.c` 是一支測試程式，當中利用 `struct treeint_ops` 此結構體建立 function hook 並呼叫 `treeint_xt` 當中的函式。
* `treeint_xt.[ch]` 包含了五個開放給使用者使用的 API ，它們都是對應 xtree 操作的包裝，會透過呼叫 xtree 的函式來進行對應操作。
Xtree，AVL tree 和 red-black tree 皆為 self-balanced binary tree ，它使用一個變數 `hint` 來判斷是否有子樹不平衡，特別注意當我們建立一個 `struct xt_tree` 後，該結構實際包含兩個 pointer to function 成員，讓我們自行實作用來建立節點和刪除節點的函式。
* `xt_destroy()`
遞迴呼叫 `__xt_destroy()` 來將整顆樹刪除
* `xt_first(n)`
節點 `n` 代表的子樹中 preorder traversal 的第一個節點
* `xt_last(n)`
節點 `n` 代表的子樹中 preorder traversal 的最後一個節點
* `xt_rotate_left(n)`
對節點 `n` 執行以下操作
```graphviz
digraph structs {
node[shape=circle];
struct0 [label= "p"];
struct1 [label= "n"];
struct2 [label= "l"];
struct3 [label= "r"];
struct4 [label= "ll"];
struct5 [label= "lr"];
node[shape=rarrow];
struct6 [label= ""];
node[shape=circle];
struct7 [label= "p"];
struct11 [label= "ll"];
struct8 [label= "n"];
struct9 [label= "l"];
struct12 [label= "lr"];
struct10 [label= "r"];
struct0 -> struct1;
struct1 -> struct2;
struct1 -> struct3;
struct2 -> struct4;
struct2 -> struct5;
struct7 -> struct9;
struct8 -> struct12;
struct9 -> struct8;
struct9 -> struct11;
struct8 -> struct10;
}
```
* `xt_rotate_right(n)`
對節點 `n` 執行以下操作
```graphviz
digraph structs {
node[shape=circle];
struct0 [label= "p"];
struct1 [label= "n"];
struct2 [label= "l"];
struct3 [label= "r"];
struct4 [label= "rl"];
struct5 [label= "rr"];
node[shape=rarrow];
struct6 [label= ""];
node[shape=circle];
struct7 [label= "p"];
struct11 [label= "n"];
struct8 [label= "rr"];
struct9 [label= "r"];
struct12 [label= "l"];
struct10 [label= "rl"];
struct0 -> struct1;
struct1 -> struct2;
struct1 -> struct3;
struct3 -> struct4;
struct3 -> struct5;
struct7 -> struct9;
struct11 -> struct12;
struct9 -> struct8;
struct9 -> struct11;
struct11 -> struct10;
}
```
* `xt_update()`
透過 `xt_balance()` 計算節點 `n` 左右子節點的 `hint` 差值，若小於負一則進行 `xt_rotate_right` ，若大於一則進行 `xt_rotate_left` ，最後檢查 `n` 的 `hint` 值是否改變還有 `n` 是否變成 leaf node ，若任一者成立則繼續對 `n` 的親代節點進行 `xt_update()`
* `xt_insert()`
嘗試插入一個新節點並且鍵值為 `key` ，首先利用 `__xt_find()` 尋找當前的樹中是否已經存在鍵值為 `key` 的節點，若存在則不進行插入，若不存在， `__xt_find()` 的過程也會順便將節點應該擺放的位置給找好 (由 `p` 和 `d` 共同紀錄) ，之後利用 `__xt_insert()` 進行節點的插入，最後透過 `xt_update()` 檢查樹的平衡狀態並更新。
* `xt_delete()`
先透過 `xt_find()` 找出要刪除的節點，注意到 `xt_find()` 並非呼叫 `__xt_find()` 而是呼叫 `__xt_find2()` ，過程當中省去了紀錄親代節點和尋找方向的紀錄。之後利用 `__xt_remove()` 來進行移除操作，如果要移除的節點 `del` 具有右子節點，則使用在右子節點代表的子樹中 inorder traversal 第一個找到的節點來替代 `del` 的位置，若無則找左子節點代表的子樹中 inorder traversal 最後一個找到的節點來替代（此處原本程式碼的註解有誤，在第 327 行應該是 "the last node found in the left subtree." 而非 "the first node found in the left subtree."），如果 `del` 是 leaf node 則直接移除並對親代節點進行 `xt_update()` 。
:::warning
修正註解的 commit
> [commit e1e9c68](https://github.com/vax-r/bitwise_lab/commit/e1e9c68e83d68aae947a80220a3e050fd9c4510f)
:::
所有 BST 不管是 AVL tree, Red-Black tree 還是一般的 binary search tree ，都可以定義出一個遞迴結構也就是說當 `root` 代表的樹是 BST 時， `root->left, root->right` 也都會是 BST ，同時可能還有數項定義需要遵守， Xtree 是否也有同樣的遞迴結構？有的話我應該要歸納出完整定義。再者，進行插入、刪除等操作的時候，都應該保持幾項 loop invariant ， Xtree 操作的 loop invariant 是什麼？
#### 實作改進
了解 Xtree 各項操作與定義的同時，我發現幾個可以改進的地方。
進行 `xt_insert()` 後會嘗試進行 `xt_update()` ，由於 binary search tree 的特性，新插入的節點在 `xt_update()` 之前一定是 leaf node ，也就是呼叫 `__xt_insert()` 一定會導向 `n->hint == 0` 接著對其親代節點進行 `xt_update()` ，我們可以直接對親代節點做 `xt_update()` 以省去多餘函式呼叫與檢查時間 ( leaf node 檢查 balance factor 一定平衡) 。
```diff
xt_parent(n) = p;
- xt_update(root, n);
+ xt_update(root, p);
```
再來發現一個問題，每當插入的節點 `n` 使得某個子樹結構如下時，會進行旋轉但旋轉後的高度差其實並沒有改變。
```graphviz
digraph structs {
node[shape=circle];
struct1 [label="a"];
struct2 [label="b"];
struct3 [label="c"];
struct4 [label="d"];
struct5 [label="e"];
struct6 [label="n"];
struct1->struct2;
struct1->struct3;
struct2->struct4;
struct2->struct5;
struct5->struct6;
node[shape=rarrow];
struct7 [label=""];
node[shape=circle];
struct8 [label="b"];
struct9 [label="d"];
struct10 [label="a"];
struct11 [label="e"];
struct12 [label="c"];
struct13 [label="n"];
struct8->struct9;
struct8->struct10;
struct10->struct11;
struct10->struct12;
struct11->struct13;
}
```
此時只要用以下方式新增節點 `f` 和 `g` 就能不斷重構上述的子樹結構
```graphviz
digraph structs {
node[shape=circle];
struct8 [label="b"];
struct9 [label="d"];
struct10 [label="a"];
struct11 [label="e"];
struct12 [label="c"];
struct14 [label="f"];
struct13 [label="n"];
struct15 [label="g"];
struct8->struct9;
struct8->struct10;
struct10->struct11;
struct10->struct12;
struct11->struct13;
struct11->struct14;
struct13->struct15;
}
```
在 `g` 插入的之後會進行旋轉，旋轉後的 `xt_update` 會停在節點 `e` 於是這個 Xtree 結構會一直傾斜下去。舉個例子，如果插入節點的順序是 [1000, 500, 300, 0, 400, 350, 450, 425, 475] ，最後產生的 Xtree 會長的像下圖
```graphviz
digraph structs {
node[shape=circle];
struct0 [label= "300"];
struct1 [label= "0"];
struct2 [label= "400"];
struct3 [label= "350"];
struct4 [label= "500"];
struct5 [label= "450"];
struct6 [label= "1000"];
struct7 [label= "425"];
struct8 [label= "475"];
struct0->struct1;
struct0->struct2;
struct2->struct3;
struct2->struct4;
struct4->struct5;
struct4->struct6;
struct5->struct7;
struct5->struct8;
}
```
`n->hint` 代表的是該節點到最遠 leaf node 的邊數， Xtree 希望保持的是左右子節點的 `hint` 差距不超過 1 ，藉此保證左右子數高度的平衡，但在上述例子當中，由於進行旋轉之後 500 的 `hint` 依然保持 2 ，所以不會繼續往上更新親代節點的 `hint` 也就無從得知旋轉後的樹依然是不平衡。事實上解決方法也不該繼續往上更新，如此一來會造成無窮次數的旋轉因為怎麼左右轉都是不平衡的，原因出在造成不平衡的來源是 `n, g` 節點的子樹，但旋轉後 `n, g` 被放入另一邊子樹下面，所以整棵樹依舊不平衡，必須先檢查造成不平衡的原因何在。
解決方法像 AVL tree 一樣將不平衡的類型分類成 RR, RL, LR, LL 四種，如果是 RL 或 LR 則要進行兩次旋轉，例如 RL 要先轉成 RR 然後再轉一次。此處只需要在發現不平衡時再子節點得兩個子節點的 `hint` 值何者較大即可判斷。將程式碼修改為以下
```diff
if (b < -1) {
/* leaning to the right */
+ if (xt_balance(xt_right(n)) > 0)
+ xt_rotate_left(xt_right(n));
if (n == *root)
*root = xt_right(n);
xt_rotate_right(n);
}
else if (b > 1) {
/* leaning to the left */
+ if (xt_balance(xt_left(n)) < 0)
+ xt_rotate_right(xt_left(n));
if (n == *root)
*root = xt_left(n);
xt_rotate_left(n);
}
```
如此一來即可判斷傾斜狀況是何種，例如若 `b < -1` 代表右傾，則進一步檢查 `xt_right(n)` 是左子樹較高還是右子樹，若是左子樹則先對 `xt_right(n)` 進行左旋使得右子樹較高然後對 `n` 進行右旋，左傾也是同樣的邏輯。
進行以上變更後，比較變更前後插入 [1000, 500, 300, 0, 400, 350, 450, 425, 475, 460] 建立的 Xtree 搜尋節點 460 的平均時間（各執行十次取平均）。
變更前 : 99.4
變更後 : 69.2
可以看到改動前後對於距離 root 最遠節點的搜尋時間大幅下降。
> [commit 5fbe5e4](https://github.com/vax-r/bitwise_lab/commit/5fbe5e48805eba330abe8b884c4a414b20513b0d)
不過比較循序輸入與亂序輸入時， `root` 節點的 `hint` 是否正確記錄樹高時，可以發現循序輸入的 `hint` 可以正確反映樹高，亂序卻沒辦法。此處代表實作尚有需要改進的地方，又或者這不是錯誤，而是我該明確定義 Xtree 的平衡性質。
```shell
$ ./build/treeint 100 0
root hint : 6
xtree height 6
$ ./build/treeint 100 3
root hint : 5
xtree height 7
```
#### Compare with AVL tree and Red-Black Tree
和 AVL tree 或 Red-Black Tree 進行比較之前，我們需要先搞清楚這些二元搜尋樹跟一般的搜尋樹有何不同，二元搜尋樹重要的是搜尋時間，如何最小化搜索時間會是關鍵，我們知道按照 General BST 的建構方式，有可能建構出一棵非常傾斜的樹而使搜尋時間跟在一個一維陣列當中搜索一樣，最好的情況則是建立一棵 complete binary tree ，在 complete binary tree 當中第 $i$ 個節點會位在第 $\lfloor lg(i) \rfloor$ 層，也就是從 root 到第 $i$ 個節點共有 $\lfloor lg(i) \rfloor$ 條邊，總共會有 $\sum_{i=1}^{n} \lfloor lg(i) \rfloor = O(n\cdot lg(n))$ 條邊。
如果給定的元素集合為 $a_1, a_2, a_3, \cdots, a_n$ 且 $a_1 < a_2 < \cdots < a_n$ ，則對於此集合來說可以不同的插入順序可以建構出許多不同的二元搜尋樹，所謂的 **Optimal Binary Search Tree** 定義如下，如果搜索元素 $a_i$ 的機率為 $p_i$ ，則搜索某顆 BST 的代價會是
$$
\sum_{i=1}^{n} p_i \cdot level(a_i)
$$
同時我們需要考慮到搜索失敗的可能，假設每個 leaf node 都另外延伸出兩個子節點當作 failure node 也就是搜索失敗會碰到的節點，則我們會有 $n + 1$ 個 failure node ，假設搜索到 $f_i$ 的機率是 $q_i$ 則搜索失敗的機率是
$$
\sum_{i=0}^{n} q_i \cdot (level(f_i) - 1)
$$
**OBST** 的定義是能使得 $\sum_{i=1}^{n} p_i \cdot level(a_i) + \sum_{i=0}^{n} q_i \cdot (level(f_i) - 1)$ 最小化的二元搜尋樹，同時滿足 $\sum_{i=1}^{n} p_i + \sum_{i=0}^{n} q_i = 1$ 。
:::info
證明 external node 數量是 internal node + 1
假設二元樹當中總共有 $n$ 個節點， external node 有 $e$ 個， internal node 有 $i$ 個，則 $n = i + e$ 且因為每個 internal node 都有兩條邊所有邊數 $E = 2i$ 且 $E = V - 1$ 所以 $2i = V - 1$ 得到 $V = 2i + 1 = i + e$ 所以 $i + 1 = e$ 得證。
:::
**Balanced BST** 則是保證此 BST 的 worst case search time complexity 與 average case search time complexity 都保持在 $O(log(n))$
**AVL tree** 的定義是樹 $T$ 若為 height-balanced 若且唯若
1. $T_L, T_R$ 皆為 height-balanced
2. $| h_L - h_R | \le 1$
保證左右子樹的高度差距不超過 1 ，這點和 Xtree 相同，不同的是 AVL tree 會在每個節點當中維護一個成員 `bf` 稱為 balanced factor ，用來記錄當前左右子樹的高度差， Xtree 紀錄的 `hint` 理論上是節點距離最遠 leaf node 的邊數，只要能保證紀錄的值正確（思考如何證明，數學歸納法或許可行），函式 `xt_balance` 就能計算出如同 AVL tree 當中的 `bf` 值。 AVL tree 保證每次 insertion 所做的旋轉次數頂多 2 次， Xtree 能否證明類似性質？
實作 AVL tree 的 insertion 並比較循序與亂序輸入情況下， AVL tree 和 Xtree 的樹高差異。 AVL tree 實作使用 [AVL Tree Program in C](https://www.javatpoint.com/avl-tree-program-in-c) 當中的程式，不過他的實作在節點數量多於 1000 時會因為過多遞迴呼叫造成 segmentation fault ，應改為非遞迴實作。以下測試 100 個節點循序插入與亂序插入
```shell
$ ./build/treeint 100 0
xtree height 6
avl tree height 6
$ ./build/treeint 100 5
xtree height 8
avl tree height 6
```
可以觀察到 Xtree 的實作在亂序輸入的時候樹高會達到 8 ， AVL tree 樹高僅有 6 符合平衡樹樹高 $\lfloor lg(n) \rfloor$，但 Xtree 不符合，雖然沒有嚴謹符合 $\lfloor lg(n) \rfloor$ 的樹高，但也比 Red-Black tree 的最糟情況 $2 * \lfloor lg(n) \rfloor$ 還少。
**Red-black tree**
紅黑樹遵守五個基本定義
1. 所有節點非黑即紅
2. root 一定是黑色
3. 所有 leaf node 一定是紅色
4. 紅色節點的子節點一定是黑色
5. 對於所有節點來說，所有該節點子樹的 leaf node 所接觸到的黑色節點數量一樣多 ( black-height 固定)
參考 [Red-Black Trees](https://www.usna.edu/Users/cs/crabbe/SI321/current/red-black/red-black.html) 內容所述，若把紅黑樹當中所有紅色節點抽離，我們可以得到一棵所有 leaf node 都在同一層的樹，間接證明了紅黑樹的 black-height 是 $O(log(n_b))$ ，最糟的情況也就是紅黑節點相間的話，紅黑樹的高度會是 $2 \cdot log(n_b)$ 也就是兩倍的 black-height 。紅黑樹的 insertion 、 deletion 進行的 fixup （變色或旋轉）有可能持續往 root 進行，此部分和 Xtree 有可能不斷往親代節點做修正類似，但紅黑樹在修正部分也有證明是為了修正哪幾項性質，完成了就不會再往上修正，例如 insertion 時還保證不是違反性質 2 就是 4 ，所有操作都是為了重新讓子樹維持這兩個性質，並且保證最多旋轉兩次。
Xtree 能否推導證明修正後能保照 Xtree 性質被維護？
> 除了上述性質外，針對 Linux 核心的場景，要特別考慮到 [Augmenting Red Black Trees](https://medium.com/swlh/augmenting-red-black-trees-d9b4cd7635f8)，你常聽到的 AR/VR 其全銜的 Augmenting，同一個字。
> :notes: jserv
閱讀 [Augmenting Red Black Trees](https://medium.com/swlh/augmenting-red-black-trees-d9b4cd7635f8) ， Augmenting Data Structure 是在闡述如何擴增原本的資料結構，例如增加成員或函式，使得該資料結構能夠在不改變原本操作的時間複雜度情況下，解決原本解決不了的問題。
以紅黑樹為例，如果要在紅黑樹中找到第 k 大的元素，原本的紅黑樹必須使用 in-order traversal 找到第 k 個存取的節點， worst case time complexity 會是 $O(n)$ ，如果在節點當中新增一個成員 Rank ，用來記錄該節點所形成的子樹的節點數量。尋找第 k 大的元素方法可更改為，若 `k == left.size + 1` 則代表該節點即是第 k 大的元素，如果小於該數值則往左子樹尋找第 k 大元素，否則往右子樹尋找第 `k - left.size + 1` 大的元素。
#### Compare with Linux kernel red-black tree
* [Red-black Trees (rbtree) in Linux](https://www.kernel.org/doc/Documentation/core-api/rbtree.rst)
* [/include/linux/rbtree.h](https://elixir.bootlin.com/linux/v6.8/source/include/linux/rbtree.h)
* [/lib/rbtree.c](https://elixir.bootlin.com/linux/v6.8/source/lib/rbtree.c)
* [/include/linux/rbtree_augmented.h](https://elixir.bootlin.com/linux/v6.8/source/include/linux/rbtree_augmented.h)
Linux 核心 紅黑樹實作可以參考以上四份文件，雖然我覺得程式碼十分混亂例如 rbtree_augmented.h 重複定義了數個在 rbtree.h 就定義的巨集。嘗試引入此紅黑樹實作並和 Xtree 進行比較。
首先定義自己的結構體並把 `struct rb_node` 嵌入其中
```c
struct rbtree_node {
struct rb_node node;
int value;
};
struct rbtree_head {
struct rb_root root;
};
```
利用 [/include/linux/rbtree.h](https://elixir.bootlin.com/linux/v6.8/source/include/linux/rbtree.h) 開放的函式 `rb_find_add()` 與 `rb_find()` 分別實作 insertion 和 find 函式。
> [commit ade5177](https://github.com/vax-r/bitwise_lab/commit/ade5177e968d3229ecce3a3d13839a9a992a89ae)
較為困難的地方在於實作 remove 函式，需要利用到 `rb_erase()` ，此函式事實上包含了 `__rb_erase_augmented()` ，可以注意到在使用 `__rb_insert(), __rb_erase_augmented(), __rb_erase_color()` 利用到了 `dummy_callbacks` ，這是因為在一般紅黑樹 (非 augmented red-black tree) 當中不需要自行實作 `propagate, copy, rotate` 等函式。
比較 Linux 核心 紅黑樹在亂序輸入與循序輸入的插入、搜尋與移除時間，對比 Xtree 如以下
```shell
$ ./build/treeint 1000 0
Red-Black Tree
Average insertion time : 122.658000
Average find time : 46.012000
Average remove time : 43.772000
XTree
Average insertion time : 247.713000
Average find time : 124.979000
Average remove time : 144.251000
$ ./build/treeint 1000 10
Red-Black Tree
Average insertion time : 132.087000
Average find time : 56.848000
Average remove time : 52.963000
XTree
Average insertion time : 242.538000
Average find time : 134.666000
Average remove time : 176.729000
```
可以發現不管是循序輸入還是亂序輸入， Linux 核心版本的紅黑樹平均執行不同操作的時間都遠小於 Xtree 。這裡讓我很好奇， Xtree 的實作應該是比較貼近 AVL tree 的，在插入時間上比紅黑樹長我認為不足為奇，但在搜尋與移除的操作上都是紅黑樹表現更為優秀，這令我思考是 Xtree 維護樹高的操作不正確亦或者是 Linux 核心的紅黑樹設計上有何特別之處？ （此處討論的都是沒有 augmented 的紅黑樹）。
> [commit f1fd1de](https://github.com/vax-r/bitwise_lab/commit/f1fd1de97f0f90c89e69ba0b0b444a39d682d91e)
Google Drive
Gist
Clipboard
Sign in with Wallet
