2024q1 Homework4 (quiz3+4)

contributed by < Appmedia06 >

第三周測驗題

測驗一

完整題目

• 實作計算開平方根
• 解釋上述程式碼運作原理並重新整理數學式 (原題目故意略去某些細節)，並嘗試用第 2 週測驗題提到的 ffs / fls 取代 __builtin_clz ，使程式不依賴 GNU extension，且提供分支和無分支 (branchless) 的實作。
• 在 Linux 核心原始程式碼找出對整數進行平方根運算的程式碼，並解說其應用案例，至少該包含 block 目錄的程式碼。

版本一

第一個版本，我們將使用不斷迭代累加的方法來找到一個數的平方根。首先，我們需要確定一個初始值。一個有效的方法是從最高有效位(Most Significant Bit)開始，逐步將其減半，這樣每次迭代都能將範圍縮小一半，類似於二分搜尋的應用。為了找到最高有效位，我們可以使用 log2() 函數。根據計算出的 MSB，我們將變數 a 賦值為在最高有效位位置上為 1 的數。

接下來，我們進行迭代。每次迭代中，result + a 代表當前嘗試的平方根值。如果這個值的平方小於或等於 N，那麼這個值就是合理的，我們將其加到 result 上。每次迭代後，我們將 a 右移一位，相當於將其除以 2。這樣，我們逐步縮小範圍，最終得出平方根的近似值。

#include <math.h>
int i_sqrt(int N)
{
    int msb = (int) log2(N);
    int a = 1 << msb;
    int result = 0;
    while (a != 0) {
        if ((result + a) * (result + a) <= N)
            result += a;
        a >>= 1;
    }
    return result;
}

版本二

和版本一的思維方式是一樣的，只是變成透過迴圈以及位元做來取得 MSB 。這樣做的理由為：

• 可移植性: 不需要使用 math.h 函數庫，減少了對外部庫的依賴。無需依賴特定平臺對 log2() 函式的實現。
• 精度: 直接基於位操作進行計算，不涉及浮點數運算，避免了可能的浮點數精度問題

int i_sqrt(int N)
{
-   int msb = (int) log2(N);
+   int msb = 0;
+   int n = N;
+   while (n > 1) {
+       n >>= 1;
+       msb++;
+   }
    int a = 1 << msb;
    int result = 0;
    while (a != 0) {
        if ((result + a) * (result + a) <= N)
            result += a;
        a >>= 1;
    }
    return result;
}

版本三

版本三使用 Digit-by-digit calculation 來實作開平方根。以下說明:

要對

$x$ 計算其平方根

$N$ ，我們可以寫成

$x=N^2$ 。假設

$N^2=(21{)}^2$ ，寫成二進制就會變成是

$N^2=(10101{)}^2$ 。將其用十進制展開，就會得到

$N^2=(2^4+0+2^2+0+2^0{)}^2$ ，可以注意到這邊的

$a_i=2^i$ or

$0$ 。也就是說我們可以將

$N^2$ 寫成，其中

$n$ 是位元數量:

$$N^2=(a_n+a_{n-1}+a_{n-2}+...+a_0{)}^2,a_m=2^{m}or{a}_m=0$$

可以先把

$N^2$ 整理一下:

$$N^2=(\sum _{i=0}^n{a}_i{)}^2$$

利用

$(x+y{)}^2=x^2+2xy+y^2$ 將上式展開， 令

$N^2=(a_0+\sum _{i=1}^n{a}_i{)}^2$:

$$N^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=a_n+a_{n-1}+...+a_m$ ，使

$P_m$ 為從索引 m 到 n 的

$a_i$ 的累加。這樣就可以再將

$P_m$ 上式代入上式:

$$N^2=P_0^2=a_0^2+2a_0{P}_1+P_1^2$$

這樣我們就得到了

$P_m$ 對應的關係式:

$$P_m^2=a_m+2a_m{P}_{m+1}+P_{m+1}^2$$

嘗試整理方程式，將

$a_m$ 提出來:

$$P_m^2=a_m(a_m+2P_{m+1})+P_{m+1}^2$$

我們令

$Y_m=a_m(a_m+2P_{m+1})$ ，代入後得:

$$P_m^2=Y_m+P_{m+1}^2$$

做到這裡，便可以從

$m=n$ 到

$m=0$ 一路去計算，透過

$P_{m+1}$ 來得到

$P_m$ ，得到

$P_m$ 要幹嘛呢？ 我們說

$P_m$ 是從索引 m 到 n 的

$a_i$ 的累加。因此，只要判斷

$P_m^2\le N^2$ 是否成立，就可以知道

$a_m$ 是否為

$2^m$ 還是為

$0$ 。

然而，如果要每次都計算

$P_m^2$ 成本太高了。因此令

$X_m=N^2-P_m^2$ ，就可以結合上式變成:

$$X_m=N^2-P_m^2=N^2-(Y_m+P_{m+1}^2)$$

再把它改寫為只有

$X$ 和

$Y$ 的樣子:

$$X_m=X_{m+1}-Y_m$$

所以我們將紀錄上一輪

$P_{m+1}$ 來計算這一輪的

$Y_m$ 。

到這邊在複習一下

$Y_m=a_m(a_m+2P_{m+1})=a_m^2+2P_{m+1}{a}_m$ ，而

$a_m=2^m$ 。因此可以將

$Y_m$ 透過一些指數操作拆成

$c_m$ 和

$d_m$ 。

$$c_m=2P_{m+1}{a}_m=2P_{m+1}{2}^m=P_{m+1}{2}^{m+1}$$

$$d_m=a_m^2=(2^m{)}^2=2^{m+1}$$

最終就可以得到:

$$Y_m=\{\begin{array}{ll}{c}_m+d_m& \text{if }{a}_m=2^m\\ 0& \text{if }{a}_m=0\end{array}$$

從上面

$Y_m$ 的式子就可以知道，我們可以從

$c_m$ 和

$d_m$ 推出

$c_{m-1}$ 和

$d_{m-1}$ 就可以來判斷

$2^m$ 的值。

$$c_{m-1}=2P_m{a}_{m-1}=P_m{2}^m=(P_{m+1}+a_m)2^m=P_{m+1}{2}^m+a_m{2}^m=\{\begin{array}{l}{c}_m/2+d_m ,if a_m=2^m\\ c_m/2 ,if a_m=0\end{array}$$

$$d_{m-1}=a_{m-1}^2=(2^{m-1}{)}^2={\displaystyle \frac{(2^m{)}^2}{2^2}}={\displaystyle \frac{d_m}{4}}$$

算到這裡終於可以來找我們的解答啦，也就是

$P_0=\sum _{i=m}^n{a}_i=a_n+a_{n-1}+...+a_0$ 。

要從

$m=n$ 到

$m=0$ 來計算，我們需要先設定初始值:

• $c_n$ 初始化
  • 因為
    $P_{n+1}=0$ ，所以
    $c_n$ 也就為
    $0$ 。程式碼中 z 就是
    $c_n$。

int z = 0;

• $d_n$ 初始化
  • $d_n=a_n^2=(2^n{)}^2=4^n$。程式碼中 m 就是
    $d_n$。
  • __builtin_clz(x) 函式回傳 x 的最高有效位前面連續為 0 的數量(leading zero)，因此 31 - __builtin_clz(x) 就是最高有效位(MSB)的位置。
  • & ~1UL 的作用就是將計算出來的 MSB 變成偶數。因為
    $d_n=4^n$ 是不可能為奇數的。
  • 將計算出的偶數位移量左移 1UL ，從而得到一個二進制，只有最高有效位為 1。

int m = 1UL << ((31 - __builtin_clz(x)) & ~1UL)

• 迭代
  • 程式碼中 b 就是
    $Y_m$。
    $Y_m=c_m+d_m$
  • $c_{m-1}$ 需要
    $c_m/2$ ， z >>= 1 。
  • $d_{m-1}={\displaystyle \frac{d_m}{4}}$ ，因此在迴圈結束前 m >>= 2 。
  • 判斷
    $X_m\ge Y_m$ 是否成立，若成立，則下面兩個方程式的上面條件成立。將
    $X_m$ 減去
    $Y_m$ 並使
    $c_{m-1}=c_m/2+d_m$ 。
    
    $$Y_m=\{\begin{array}{ll}{c}_m+d_m& \text{if }{a}_m=2^m\\ 0& \text{if }{a}_m=0\end{array}$$

$$c_{m-1}=\{\begin{array}{l}{c}_m/2+d_m ,if a_m=2^m\\ c_m/2 ,if a_m=0\end{array}$$

for (int m = 1UL << ((31 - __builtin_clz(x)) & ~1UL); m; m >>= 2) { 
    int b = z + m;
    z >>= 1; // BBBB = 1
    if (x >= b)
        x -= b, z += m;               
}

完整程式碼:

int i_sqrt(int x)
{
    if (x <= 1) /* Assume x is always positive */
        return x;

    int z = 0;
    for (int m = 1UL << ((31 - __builtin_clz(x)) & ~1UL); m; m >>= 2) {
        int b = z + m;
        z >>= 1;
        if (x >= b)
            x -= b, z += m;               
    }
    return z;
}

用 ffs, fls 取代 __builtin_clz

• ffs: 找到 x 中最低有效位為 1 的位置。

The ffs(), ffsl() and ffsll() functions find the first bit set (beginning with the least significant bit) in value and return the index of that bit.

• fls: 找到 x 中最高有效位為 1 的位置。

The fls(), flsl() and flsll() functions find the last bit set in value and return the index of that bit.

從上面可以知道

__builtin_clz(x) 函式回傳 x 的最高有效位前面連續為 0 的數量(leading zero)。
因此 31 - __builtin_clz(x) 就是最高有效位(MSB)的位置

因此透過 fls(x) 函式找到 x 中有效位為 1 的位置。但因為 fls(x) 是從 1 開始算到第 n 位。因此需要將其減一，這樣就可以得到和 31 - __builtin_clz(x) 一樣的數值了。

-for (int m = 1UL << ((31 - __builtin_clz(x)) & ~1UL); m; m >>= 2) {
+for (int m = 1UL << ((fls(x) - 1) & ~1UL); m; m >>= 2) {

Linux 核心的 buffered writeback throttling

用於管理和控制緩存寫回操作速度的一種機制。這種機制確保系統在進行大量數據寫入磁碟
時，不會因為寫入操作過多而導致系統性能下降或資源耗盡。

當應用程式寫入數據到檔案時，數據首先寫入到頁緩存（page cache）中，而不是直接寫入到磁碟。這種方式被稱為 "緩衝寫入"（buffered write）。這樣可以加快寫入速度，因為寫入緩存比寫入磁盤要快得多。如下圖所示。

然而，這些緩存數據最終需要寫回到磁盤上，這個過程稱為 "寫回"（writeback）。如果緩存數據積累過多，寫回操作會佔用大量 I/O 頻寬，影響系統其他操作的性能。

linux/block/blk-wbt.c

Linux 中的 buffered writeback throttling 是一種基於 CoDel（Controlled Delay）的機制，但它針對的是 I/O 調度，而不是網絡封包。因此，它不能通過丟棄封包來控制延遲，而是通過調整 I/O 排隊深度來實現。

系統會定義好一個時間窗口內監控 I/O 操作的延遲，而延遲就是從寫入緩存到實際寫入磁盤的時間。如果時間窗口內的最小延遲超過目標值（target），則放大縮放步驟（scaling step），並將排隊深度（queue depth）縮減一半（2x）。縮放步驟越大，縮小時間窗口的速度越快。監控窗口的大小將縮小為 100 / sqrt(scaling step + 1) 。也就用到本題的開平方根了。

在 linux/block/blk-wbt.c 中的 rwb_arm_timer 函式為根據 I/O 請求的當前狀態來設置定時器的時間窗口。在第 12, 13 行，根據 scale_step 的值來動態調整當前時間窗口（cur_win_nsec）的大小。而將其左移 8 位，也就是乘以 256 是為了在接下來的平方根計算中保持數值的整數形式。

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);
}

測驗二

完整題目

• 將 除以 10 (div 10) 以及 取餘 10 (mod 10) 合併成一個函式並降低運算成本。
• 參照《Hacker's Delight》和 CS:APP 第二章，解釋上述程式碼運作原理，並對照 Instruction tables，以分析最初具備除法操作的程式碼和最終完全不依賴除法和乘法的實作程式碼裡頭，指令序列佔用的 CPU 週期數量。

算法推導

題目的要求是要降低除以 10 (div 10) 以及 取餘 10 (mod 10) 合併成一個函式並降低運算成本。以最簡單的方式就是直接使用 / 和 % 運算子來操作。對應到下面的程式碼。

carry = tmp / 10;
tmp = tmp - carry * 10;

然而，我們是不是可以提出一種方式不要使用到 / 和 % 運算子就可以達到目的。參考《Hacker's Delight》，採用 bitwise operation 來求得商和餘數。

和數位邏輯設計中的減法操作想法類似，如果不想要有除法操作，就換成乘法操作就好啦。因此可以把 (除以

$10$) 改成 (乘以

$\frac{1}{10}$ )。但是

$10$ 不符合

$2^k\pm 1$ 的形式，換句話說，因為

$10$ 有包含

$5$ 這個因數，無法完全用

$2$ 的冪項來表示，因此無法以二進位表示其數值。因此可以取一個非常近似於

$\frac{1}{10}$ 的數值

$\frac{a_n}{2^N}$ 的倒數來表示。此題的

$2^N=128,a_n=13$ ，所以

$\frac{2^N}{a_n}=\frac{128}{13}\simeq 9.84$ 。

有了想法後就要開始思考要如何將被除數乘以

$\frac{13}{128}$ 。

• 要如何乘以
  $13$ 呢?

因為

$13$ 並不是

$2$ 的幂，無法直接用簡單直觀的 bitwise 操作。但仍然能透過多個

$2$ 的幂的位移操作的累加來湊到

$13$。

$$\frac{tmp}{8}+\frac{tmp}{2}+tmp=\frac{13\cdot tmp}{8}$$

對應的程式碼為:

((tmp >> 3) + (tmp >> 1) + tmp)

但可以看到如果直接右移 3 位，會導致這三位的位元直接不見。因此利用 3 個變數分別用 mask 儲存這三位的位元。

d0 = q & 0b1;
d1 = q & 0b11;
d2 = q & 0b111;

最後先左移 3 位(乘以 8) ，並將這三位加上去後就可以得到

$13\cdot tmp$

((((tmp >> 3) + (tmp >> 1) + tmp) << 3) + d0 + d1 + d2)

• 除以
  $128$

因為

$128$ 是

$2$ 的幂 ，因此只要將其右移 7 位就可以了。這樣就把商給求出來了。

q = ((((tmp >> 3) + (tmp >> 1) + tmp) << 3) + d0 + d1 + d2) >> 7;

• 求餘數

根據餘式定理： 被除數

$=$ 除數

$\times$ 商數

$+$ 餘數
因此 餘數

$=$ 被除數

$-$ 除數

$\times$ 商數 ，在本題中：餘數

$=$ 被除數

$-$

$10$

$\times$

$q$

和剛剛一樣的問題，要如何乘以

$10$ 呢？ 要用

$2$ 的幂的倍數去湊，我們可以寫成：

$$10\times q=(q\times 4+q)\times 2$$

對應的程式碼為:

(((q << 2) + q) << 1);

最後被除數 tmp 減去

$10\times q$ 即為餘數。

r = tmp - (((q << 2) + q) << 1);

包裝函式

現在要包裝成函式，但可以看到實作和上面說明不一樣。

uint32_t x = (in | 1) - (in >> 2); /* div = in/10 ==> div = 0.75*in/8 */

• (in | 1) : 將輸入轉成奇數，若原本就是奇數則不影響。
• (in >> 2) : 變成
  $\frac{in}{4}$
• (in | 1) - (in >> 2) : 實際上就在計算
  $\frac{3\cdot in}{4}$ ，如下式說明

$$in-\frac{in}{4}=\frac{3\cdot in}{4}$$

uint32_t q = (x >> 4) + x;

• x =
  $\frac{3\cdot in}{4}$ ， x >> 4 =
  $\frac{3\cdot in}{64}$
• 兩個相加後，就變成

$$\frac{3\cdot in}{64}+\frac{3\cdot in}{4}=\frac{102\cdot in}{128}\simeq \frac{8\cdot in}{10}$$

x = q;
q = (q >> 8) + x;
q = (q >> 8) + x;
q = (q >> 8) + x;
q = (q >> 8) + x;

• 利用 4 次的 (q >> 8) + x 來減少誤差值

*div = (q >> 3);

• 將結果除以 8 ，就得到商，也就是將上式變成:
  
  $$\frac{8\cdot in}{10}\times \frac{1}{8}=\frac{in}{10}$$

*mod = in - ((q & ~0x7) + (*div << 1));   

• q & ~0x7 : 清除 q 的最低3位。
• *div << 1 : 將商左移1位，相當於乘以2。
• in 減去這個值，得到餘數。

完整程式碼:

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));   
}

比較有無除法操作

利用 perf 工具來計算兩個操作的指令數量以及執行的 cycle 數。

• 無使用除法、取餘操作 （本題）

Performance counter stats for './divmod_10' (5 runs):

          993,919      instructions  # 1.25  insn per cycle   ( +-  0.45% )
          795,099      cycles                                                                  ( +-  2.27% )

         0.001584 +- 0.000383 seconds time elapsed  ( +- 24.21% )

• 使用除法、取餘操作

Performance counter stats for './divmod2_10' (5 runs):

          992,694      instructions  # 1.23  insn per cycle   ( +-  0.24% )
          808,583      cycles                                                                  ( +-  2.51% )

         0.001533 +- 0.000343 seconds time elapsed  ( +- 22.36% )

令我感到意外的是，無使用除法、取餘操作的指令數量以及執行的 cycle 數竟然和使用了差不多。原本的猜想是直接使用除法操作會使用大量的 CPU cycle ，實際上我自己設計 CPU 時，實作除法模塊需要的 clock cycle 就比其他指令多很多，甚至要暫停流水線來等除法指令跑完。而本題改用 bitwise 的操作來避免掉它，效能上應該要有提昇。

直到看到了 vax-r 筆記裡面老師的留言，這邊引用一下：

CPU 週期數量「沒有明顯差異」，就說明此舉的效益，當指令集受限時，勢必要調整實作手法，而且看似指令數量增加，但避開較長週期的特定指令，有機會在亂序 (out-of-order) 執行的處理器上，爭取到更高的 IPC。

jserv

也就是說雖然指令數量稍微增加了，但因為每個指令的執行時間長短不一，這樣的操作可以避免掉長週期的特定指令，就有機會在特定架構的指令集下取得更好的 IPC (Instruction per clock) 。此外，像是 RISC-V 指令集中 ， RV32I 是並沒有除法和取餘指令的（在 RV32M 中），因此這樣的操作也可以讓不支援除法指令的處理器也能實現除法。

TODO: 練習撰寫不依賴任何除法指令的 % 9 (modulo 9) 和 % 5 (modulo 5) 程式碼，並證明在 uint32_t 的有效範圍可達到等效。

測驗三

完整題目

• 解釋 ilog2 函式實作
• 在 Linux 核心原始程式碼找出 log2 的相關程式碼 (或類似的形式)，並解說應用案例

版本一

我們要計算

$lo{g}_2^i=N$ ，可以透過每次迭代時將 i 除以 2，只要能除以 2 就把 log 加 1 。要注意一開始 log 的值是設成 -1 ，使得第一輪右移操作後 log 可以為 0 ，從而符合

$lo{g}_2^1=0$ 。

int ilog2(int i)
{
    int log = -1;
    while (i) {
        i >>= 1;
        log++;
    }
    return log;
}

這個版本有個缺點，那就是時間複雜度會隨著不同輸入而改變。int 型態正數的最大值為

$2^{31}-1$ ，若輸入為

$2^{31}-1$ ，也就是最高有效位為第 31 位元時，就需要跑 31 次迴圈。

版本二

為了解決上述的問題，我們希望可以把時間複雜度變成

$O(1)$ 。通過逐步大幅度右移來快速計算整數的二進制對數。它首先將 i 右移 16 位，然後是 8 位，接著是 4 位，最後是 1 位。每次右移後，根據移位的次數來累加結果中的對數值。這種方法比簡單逐位右移的效率更高，特別適用於較大的整數。

static size_t ilog2(size_t i)
{
    size_t result = 0;
    while (i >= 65536) { // 2^16
        result += 16;
        i >>= 16;
    }
    while (i >= 256) { // 2^8
        result += 8;
        i >>= 8;
    }
    while (i >= 16) { // 2^4
        result += 4;
        i >>= 4;
    }
    while (i >= 2) { // 2^2
        result += 1;
        i >>= 1;
    }
    return result;
}

size_t 補充

size_t 實際上就是 long unsigned int 。它由 typedef 定義在 <stddef.h> 中。

#define __SIZE_TYPE__ long unsigned int
...
typedef __SIZE_TYPE__ size_t;

• 那為什麼要使用 size_t 呢？

size_t 的意思是 size type ，是一種計數類型。取值範圍與計算機架構與作業系統相關。 32 位元機器一般是 unsigned int，佔 4 位元組；而 64 位元計算機一般是unsigned long，佔 8 位元組。有就和我的電腦是一致的。

我們常用的 sizeof 運算子，以及一些字串操作都是使用 size_t 。而它就是為了表達執行中的計算機所能容納建立最大物件的位元組數大小的意義。也因為這樣設計這樣的型別就是為了能讓 C 語言的程式能夠達成跨設備的兼容性。

版本三

利用 GCC 內建函式 __builtin_clz 實作。

__builtin_clz(x) 函式回傳 x 的最高有效位前面連續為 0 的數量(leading zero)。

特別注意的是 GCC 文檔 Other Built-in Functions Provided by GCC 中提到

Built-in Function: int __builtin_clz (unsigned int x)
Returns the number of leading 0-bits in x, starting at the most significant bit position. If x is 0, the result is undefined.

若輸入 0 進去會得到 undefined befavior 。因此輸入需要或上 1 ，當輸入為 1 時，可以變成 1 。也順便對應了

$lo{g}_2^1=0$ 。

最後將 31 減掉取得的最高有效位即為對數的值了。

int ilog32(uint32_t v)
{
    return (31 - __builtin_clz(v | 1));
}

Linux 核心的 log2

在 Linux 核心中使用 log2 函數進行計算在多個應用中是常見的，因為二進制對數在計算資源分配、數據結構優化和性能調優方面非常有用。例如確定記憶體分配的級別，或是計算紅黑數的深度等等。

在 include/linux/log2.h 中實作了 log2 。想法和上面的版本三類似，用 fls(n) 找到最高有效位，將其減一，就可以算出對數了。

int __ilog2_u32(u32 n)
{
    return fls(n) - 1;
}
int __ilog2_u64(u64 n)
{
    return fls64(n) - 1;
}

測驗四

完整題目

• 解釋 EWMA (指數加權移動平均) 實作
• 在 Linux 核心原始程式碼找出 EWMA 的相關程式碼 (或類似的形式)，並解說應用案例，至少該涵蓋無線網路裝置驅動程式。

背景知識

指數加權移動平均（Exponentially Weighted Moving Average, EWMA）是一種取平均的統計手法。它可以用來消除數據中的隨機波動，突出數據的趨勢。EWMA 通過對時間序列中的數據賦予不同的權重來計算平均值，較新的數據點被賦予較高的權重，而較舊的數據點被賦予較低的權重。這種權重是指數衰減的，隨著數據點的遠離，權重指數性減小。

數學定義如下:

$$S_t=\{\begin{array}{ll}{Y}_0& t=0\\ \alpha Y_t+(1-\alpha )\cdot S_{t-1}& t>0\end{array}$$

其中:

• $S_t$ 為第
  $t$ 個時間點的 EWMA 值
• $Y_t$ 為第
  $t$ 個時間點的觀測值
• $\alpha$ 為歷史資料加權常數，介於0與1之間
  • $\alpha$ 值越大，EWMA 對最近的數據點反應越靈敏（變得更快），對歷史數據的依賴性越小。
  • $\alpha$ 值越小，EWMA 越平滑（變得更慢），對歷史數據的依賴性越大。

可以將上面式子展開成以下型態:

$$S_t=\alpha (Y_t+(1-\alpha )Y_{t-1}+(1-\alpha {)}^2{Y}_{t-2}+...+(1-\alpha {)}^k{Y}_{t-k}+...+Y_0)$$

從結果可以得知第前

$k$ 筆資料為

$(1-\alpha {)}^k{Y}_{t-k}$ 其加權為

$(1-\alpha {)}^k$ 為指數的模樣，這也是為何該手法被稱為 EWMA。

結構體

• internal:
  $S_t$, 第
  $t$ 個時間點的 EWMA 值
• factor: 用於縮放內部表示的縮放因子
  • 最高平均值公式 ULONG_MAX /（factor * weight）
• weight: 指數加權的權重，用於控制舊值影響衰減的速度

/* Exponentially weighted moving average (EWMA) */
struct ewma {
    unsigned long internal;
    unsigned long factor;
    unsigned long weight;
};

ewma_init

• 檢查 factor 和 weight 是否為 2 的冪。如果不是，程式會觸發 assert 。
  • 為了優化性能，縮放因數必須是 2 的冪。
• 使用測驗三的 ilog2 函數計算 weight 和 factor 的對數值，並將結果儲存到 avg 結構中。
  • 為了表示縮放級別，因此取其對數，在後面的操作皆由左移右移操作。
• 將
  $S_t$ 設置為 0 。

void ewma_init(struct ewma *avg, unsigned long factor, unsigned long weight)
{
    if (!is_power_of_2(weight) || !is_power_of_2(factor))
        assert(0 && "weight and factor have to be a power of two!");

    avg->weight = ilog2(weight);
    avg->factor = ilog2(factor);
    avg->internal = 0;
}

ewma_add

將一個新樣本值 val 添加到 EWMA 中，並更新結構體 ewma 的

$S_t$ (internal)。

• 當 internal 為 0 時，代表這是第一個數值。

(val << avg->factor)

新樣本值 val 左移 avg->factor 位數，這樣做是為了將新樣本值放大，適應內部表示的範圍。

• 非第一個數值時

我們將其細分為幾個部份：

1. 當前的內部值 internal 左移 weight 位數，這樣做是為了準備後續的加權計算。

avg->internal << avg->weight

2. 放大的內部值減去原始內部值，這樣做是為了計算衰減部分。

(avg->internal << avg->weight) - avg->internal

3. 將衰減後的內部值加上放大的新樣本值，這一步將新樣本值整合到加權平均值中。

((avg->internal << avg->weight) - avg->internal) + (val << avg->factor)

4. 將結果右移 weight 位數，這樣做是為了將結果縮小回適當的範圍。

(((avg->internal << avg->weight) - avg->internal) + (val << avg->factor)) >> avg->weight

ewma_read

返回當前的平均值。將 internal 右移 factor 位，獲得真實的平均值。

unsigned long ewma_read(const struct ewma *avg)
{
    return avg->internal >> avg->factor;
}

is_power_of_2

對於一個數值 x 為 2 的幂，其二進位表示中一定只會有一個 1 。則 x 減去 1 後，x 和 x-1 的位元表示將不會有共同都為 1 的位，做 AND 運算結果為0。

bool is_power_of_2(unsigned long x) {
    return (x != 0) && ((x & (x - 1)) == 0);
}

EWMA 實作

設置 factor 為 16 (縮放倍率為 4)，weight 為 4 (指數加權的權重為 2)。輸入一組數列如下:

$$10,20,30,25,15,40,35,30,45,50$$

經過 EWMA 取平均之後，便得到以下圖表：

可以看到原本數據(藍線)是比較震盪的，經過 EWMA 處理後(橘線)，明顯變得平滑許多。

Linux 核心的無線網路裝置驅動程式

DECLARE_EWMA

在 linux/include/linux/average.h 中定義了 DECLARE_EWMA 的巨集。

• 輸入三個參數:
  1. name: 根據輸入的 name 來決定結構體以及後面函式的名稱
  2. _precision: 對應到前面的 factor ，縮放因子
  3. _weight_rcp: 對應到前面的 weight ，指數加權的權重
• 和本題一樣，實作了 init, read, add 三個函式。

ath5k 是基於 MadWifi 和 OpenHAL 的 Linux 版本 Atheros 無線驅動（ FreeBSD 的版本叫做 ar5k ）。 ath5k 主要支援 AR24xx 和 AR5xxx 系列晶片。它是 Linux 無線網路子系統的一部分，提供對多種 Atheros 無線芯片的支持。ath5k 驅動程式的目標是提供穩定、高效的無線網絡支持。

ath5k.h

在 linux/drivers/net/wireless/ath/ath5k/ath5k.h 中就使用到 DECLARE_EWMA 來平滑 RSSI 值是一種常見的方法，旨在減少信號質量測量中的短期波動，使得信號強度讀數更穩定和可靠。

RSSI (Received Signal Strength Indicator)

RSSI 值表示接收到的無線信號強度，因環境影響，如障礙物、移動設備等，RSSI 值會頻繁波動。這種波動可能會影響無線網絡的性能和用戶體驗。因此，需要對 RSSI 值進行平滑處理，以便獲得更穩定的信號強度測量。

因此，對應到程式碼中的 DECLARE_EWMA 巨集，它將 factor 設置為

$10$ ，而 weight 設置為

$8$。

DECLARE_EWMA(beacon_rssi, 10, 8)

測驗五

完整題目

• 解釋 ceil_ilog2 函式實作 (
  $\lceil lo{g}_2(x)\rceil$)
• 改進程式碼，使其得以處理 x = 0 的狀況，並仍是 branchless
• 在 Linux 核心原始程式碼找出 ceil 和 log2 的組合 (類似的形式)，並解說應用案例 (提示: 在 Linux 核心排程器就有)

算法推導

要計算給定無符號32位整數 x 的

$\lceil lo{g}_2(x)\rceil$，即求 x 的二進制位數的對數，並將結果向上取整。如果 x 是 2 的冪，結果即為

$lo{g}_2(x)$；否則，結果是大於 x 的最小的 2 的冪次方的對數。

和 ilog2() 的版本二有些類似，我們一樣透過二的幂來位移數值。可以看到 r 以及 shift 都是先判斷是否大於二的幂，若是的話，則右移 x 其指數位。最後的結果將 r 和 shift 或上，因為這兩個變數一定是二的幂，且若是 x 此時還大於 1 ，則要加 1，為了符合 ceil 特性。最後在加上一開始減掉的 1。

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; // GGGG = 1     
}

改進程式碼

上面程式碼在一開始將 x-- ，這樣可能會導致當 x = 0 時，再去減一會得到 x = 0xFFFFFFFF ，這樣就導致後面計算毫無意義。因此我們改進程式碼，使其得以處理 x = 0 的狀況，並仍是 branchless 。

• 一個簡單的想法就是一開始就判斷其數值是否等於 0

- x--;
+ x = x - (x > 0);

但是我們要的是 branchless 的情況阿，因此仍需調整。

• 利用 C 語言強轉布林代數的技巧 !!x ，來判斷 x 是否為 0 。

- x = x - (x > 0);
+ x = x - (!!x);

commit 347a81d

!!x 的作用：

• 第一次 !x 將 x 轉換為其邏輯非值。
• 第二次 ! 再次將結果轉換為其邏輯非值。

這樣可以將任何非零值轉換為 1，將 0 保持為 0。

第四周測驗題

測驗一

完整題目

• 解釋 population count 函式實作。
• 不依賴 popcount，嘗試以上述歸納的規則，針對 LeetCode 477. Total Hamming Distance 撰寫出更高效的程式碼。

population count 簡稱 popcount 或叫 sideways sum，是計算數值的二進位表示中，有多少位元是 1 。

版本一

最簡單的想法就是逐一將 LSB 清除，同時累加 n 。注意到這邊有兩個 bitwise 的操作:

• v &= (v - 1): 將 LSB 設為 0 ，舉例來說 144 的 LSB 在第 4 位:

1001 0000 // 144
1000 1111 // 144 - 1
1000 0000 // 144 & (144 - 1)

• n = -(~n): 將 n 加一：
  
  $$-n=\sim n+1\Rightarrow -(\sim n)=n+1$$

unsigned popcount_naive(unsigned v)
{
    unsigned n = 0;
    while (v)
        v &= (v - 1), n = -(~n);
    return n;
}

版本二

但和前面一樣，我們希望可以不要使用到分支，也就是用常數時間複雜度來實作。因此考慮以下算法。

對於一個 32 bit 的無號整數，popcount 可以寫成以下數學式:

$$popcount(x)=x-\lfloor \frac{x}{2}\rfloor -\lfloor \frac{x}{4}\rfloor -...-\lfloor \frac{x}{2^{31}}\rfloor$$

這邊有些複雜，我們先把問題縮減。假設只看 x 的最低四位元

$x[3:0]$ ，也就可以把上式寫成:

$$popcount(x)=x-\lfloor \frac{x}{2}\rfloor -\lfloor \frac{x}{4}\rfloor -\lfloor \frac{x}{8}\rfloor$$
這樣就可以計算 4 個位元中 1 的個數了。因此在 popcount_branchless 函式的實作就是以每 4 個位元 (nibble) 為一個單位計算 1 的個數。

1. 計算
   $v-\frac{v}{2}-\frac{v}{4}-\frac{v}{8}$

n = (v >> 1) & 0x77777777;
v -= n;
n = (n >> 1) & 0x77777777;
v -= n;
n = (n >> 1) & 0x77777777;
v -= n;

bitwise 操作如下圖所示:

B31 B30 B29 B28 ... B7 B6 B5 B4 B3 B2 B1 B0 # v
  0 B31 B30 B29 ...  0 B7 B6 B5 B4 B3 B2 B1 # (v >> 1) & 0x77777777
  0   0 B31 B30 ...  0  0 B7 B6 B5 B4 B3 B2 # (n >> 1) & 0x77777777
  0   0   0 B31 ...  0  0  0 B7 B6 B5 B4 B3 # (n >> 1) & 0x77777777

舉個例子，假設有輸入 v 的 十六進制為 32'h20030416 ，換算成十進制就是 32'd537068566，再轉成二進制為 32'b0010,0000,0000,0011,0000,0100,0001,0011。

經過此操作後，數值就變成 32'b0001,0000,0000,0010,0000,0001,0001,0010，有意思的地方就是每 4 個位元表示的就是原本數值每 4 個位元 1 的數量。例如

$v[31:28]=0001$ 有一個 1，而處理的前四位的數值就是

$0001$ 也就是 1。

2. 我們需要將所有 nibble 相加。

假設

$B_n$ 代表第 n 個 nibble ，以 4 個位元為單位

B7 B6 B5 B4 B3 B2 B1 B0  // v
 0 B7 B6 B5 B4 B3 B2 B1  // (v >> 4)

將其加起來

// (v + (v >> 4))
B7 (B7+B6) (B6+B5) (B5+B4) (B4+B3) (B3+B2) (B2+B1) (B1+B0)

使用 0x0F0F0F0F 做 mask 可得

// (v + (v >> 4)) & 0x0F0F0F0F
0 (B7+B6) 0 (B5+B4) 0 (B3+B2) 0 (B1+B0)

最後我們使用乘法來取得

$B_0$ 到

$B_7$ 的加總。先令

$A_0=B_0+B_1$ 以此類推。將 v 乘上 0x01010101:

                         0 A6  0 A4  0 A2  0 A0
                      x  0  1  0  1  0  1  0  1
---------------------------------------------------
                         0 A6  0 A4  0 A2  0 A0
                   0 A6  0 A4  0 A2  0 A0  0
             0 A6  0 A4  0 A2  0 A0  0
       0 A6  0 A4  0 A2  0 A0  0
---------------------------------------------------
                            ↑_______________________A6+A4+A2+A0

可以看到在第 7 個單位，也就是第 24 位元的地方就是

$B_0$ 到

$B_7$ 的加總了。因此最終我們將結果右移 24 位，就是答案啦。

• 完整程式碼

unsigned popcount_branchless(unsigned v)
{
    unsigned n;
    n = (v >> 1) & 0x77777777;
    v -= n;
    n = (n >> 1) & 0x77777777;
    v -= n;
    n = (n >> 1) & 0x77777777;
    v -= n;

    v = (v + (v >> 4)) & 0x0F0F0F0F;
    v *= 0x01010101;                                     

    return v >> 24;
}

Hamming distance

兩個整數間的 Hamming distance 為其二進位的每個位元的差

Example :

數字(十進制）	bit3	bit2	bit1	bit0
1	0	0	0	1
4	0	1	0	0

Hamming distance = 2

參考 LeetCode 477. Total Hamming Distance :

A ^ B　就為二進位中兩數字的每個位元的差，再透過　popcount() 就可以得到答案。

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;
}

最後的返回值是 total 右移一位的值，也就是將 total 除以2。這是因為在上述雙重迴圈中，每一對數字的漢明距離被計算了兩次（nums[i] 和 nums[j] 以及 nums[j] 和 nums[i]），所以需要除以2得到正確的總和。

上面實作的時間複雜度為

$O(n^2)$ ，在 leetcode 上找到一個更快速的解法，其時間複雜度僅需

$O(32n)$ ，以下來解釋：

• total_dist 代表最終的漢明距離
• int 型態是 4 個 byte，也就是 32 個 bit。所以我們從第一個位元開始計算到第 32 個位元
• 走訪整個陣列，去計算每個數字的第 i 個位元是否是 1，若是的話就累加到 bits。
• 在每個位元的漢明距離就是有幾個數字被設為 1 乘上有幾個數字被設為 0 ，舉個例子

數字(十進制）	bit3	bit2	bit1	bit0
1	0	0	0	1
2	0	0	1	0
3	0	0	1	1
4	0	1	0	0
5	0	1	0	1

在 bit0 裡面可以看到

$\{1,3,5\}$ 是 1 ，因此 bits 就等於 3 。

$\{2,4\}$ 是 0，因此(numsSize - bits) 是 2。所以在第 0 個位元的漢明距離就是

$3\times 2=6$。以此類推，直到 32 個位元皆跑過一次。

• 完整程式碼

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; // 計算第 i 個位元是否是 1
        }
        total_dist += bits * (numsSize - bits); // 計算 hamming distance
        bits = 0;
    }
    return total_dist;
}

比較效能

使用 perf 工具來看優化後的指令數量和時脈數量。設定陣列大小為 10000 筆資料。

• 原始程式碼:
  $O(n^2)$

 Performance counter stats for './HammingDistance' (5 runs):

     4,202,435,439      instructions                     #    3.71  insn per cycle              ( +-  0.00% )
     1,131,624,752      cycles                                                                  ( +-  0.08% )

          0.271950 +- 0.000485 seconds time elapsed  ( +-  0.18% )

• 優化後:
  $O(32n)$

 Performance counter stats for './HammingDistance_improve' (5 runs):

         7,146,108      instructions                     #    1.88  insn per cycle              ( +-  0.06% )
         3,791,584      cycles                                                                  ( +-  0.43% )

          0.003398 +- 0.000130 seconds time elapsed  ( +-  3.83% )

可以看出指令和時脈數量都大量減少。

測驗二

完整題目

• 解釋 Remainder by Summing digits
• 將上述手法應用於第三次作業中，追求更有效的賽局判定機制。

本題希望提出不使用除法就可以算出餘數。

模運算和同餘關係

模運算（modulus operation）用來求一個數除以另一個數的餘數。同餘關係（congruence relation）表示兩個數在除以特定的除數後具有相同的餘數。下式代表 a 和 b 除以 m 之後的餘數相同。

$$a\equiv b(modm)$$

算法推導

以除數為 3 為例，

$1\equiv 1(mod3)$ 且

$2\equiv -1(mod3)$，，將後者不斷的乘上前者可得:

$$2^k\equiv \{\begin{array}{l}1(mod3),k\text{is even}\\ -1(mod3),k\text{is odd}\end{array}$$

有一整數 n 用二進位表示為

$b_{n-1}{b}_{n-2}...b_1{b}_0$ ，有上式推導可以將 n 改寫成：

$$n\equiv \sum _{n=0}^{n-1}{b}_i\ast (-1{)}^i\ast (mod3)$$

將每次得到的位元和遞迴的應用上式直到得到一個小於 3 的數即是餘數，若變成負數則要加上 3 的倍數。

而位元和通常又可以利用 population count 這類的函式來得到。我們可以寫成偶數位的 set bit 數量減掉奇數位的 set bit 數量。可以寫成：

• $5$ 的二進制是
  $0101$ ，和 n 做 AND 運算就可以獲得偶數位
• $A$ 的二進制是
  $1010$ ，和 n 做 AND 運算就可以獲得奇數位

n = popcount(n & 0x55555555) - popcount(n & 0xAAAAAAAA);

進一步化簡，其中 x 就是 n，而 m 就是 0xAAAAAAAA:

$$popcount(x\wedge \bar{m})-popcount(x\wedge m)=popcount(x\oplus m)-popcount(m)$$

因此可以將其改寫成:

n = popcount(n ^ 0xAAAAAAAA) - 16;

注意到經過此計算後 n 可能變成負數，所以我們要將其轉正再去取餘數。

int mod3(unsigned n)
{
    int diff = popcount(n ^ 0xAAAAAAAA) - 16;

    int result = diff % 3;
    if (result < 0) {
        result += 3;
    }

    return result;
}

但很明顯這樣太過複雜，要把分支去掉，因此修正為:

-    int result = diff % 3;
-    if (result < 0) {
-        result += 3;
-    }

-    return result;
+    return (diff % 3 + 3) % 3;

Lookup table

另一個方式就是利用查表的方式，只要計算奇數位的 set bit 去查表即可。表中的數字重複

$\{2,0,1\}$ 。

int mod3(unsigned n)
{
    static char table[33] = {2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1 };
    n = popcount(n ^ 0xAAAAAAAA);
    return table[n];
}

圈圈叉叉

模擬 100 萬局隨機圈圈叉叉遊戲，並統計和列印每個位置的勝率、每個玩家的勝率、平手次數和總耗時。

每一次動作，都會從 3x3 大小的棋盤中選一個可放的位置，一場最多就是 9 個動作。而如何知道棋盤上的動作呢？

board 是一個 uint32_t 型態存放著棋盤資訊，每次移動都會做以下操作，board |= move_masks[move]。我們知道 32 位元數值以 16 進制表示會有 8 個位元。而圈圈叉叉中勝利的線剛好就是 8 根(直線 3 根、橫線 3 根、斜線 2 根)。因此每個位元就代表一條線。

static const uint32_t move_masks[9] = {
    0x40040040, 0x20004000, 0x10000404, 
    0x04020000, 0x02002022, 0x01000200, 
    0x00410001, 0x00201000, 0x00100110,
};

透過以上遮罩做或運算，就可以得到下表:

勝利條件	board
左邊直線	0x44470041
中間直線	0x22207022
右邊直線	0x11100714
上面橫線	0x70044444
中間橫線	0x07022222
下面橫線	0x00711111
右上左下斜線	0x12412427
左上右下斜線	0x42142172

可以看出某方的勝利條件就是 board 有

$7$ 。因此檢查勝利的條件就可以寫成先將 board 每個 4 個位元組（16 進制）加一， AND 上 0x88888888 ，若有任一位元組為 8 即勝利。

static inline uint32_t is_win(uint32_t player_board)
{
    return (player_board + 0x11111111) & 0x88888888;
}

mod7

本題的原理是使用

$n(mod7)\equiv \lfloor \frac{8}{7}n\rfloor (mod8)$ 的原理來求取餘 7。可以先將 x 乘上

$\lfloor \frac{32}{7}n\rfloor$ ，也就是 0x24924925 。又知道說

$2^{32}=2^{29}\ast 8$ ，所以最終右移 29 位即為取餘 7 的答案。

static inline int mod7(uint32_t x)
{   
    x = (x >> 15) + (x & UINT32_C(0x7FFF));
    /* Take reminder as (mod 8) by mul/shift. Since the multiplier
     * was calculated using ceil() instead of floor(), it skips the
     * value '7' properly.
     *    M <- ceil(ldexp(8/7, 29))
     */
    return (int) ((x * UINT32_C(0x24924925)) >> 29);
}

測驗三

完整題目

• 解釋 XTree 運作原理
• 指出上述程式碼可改進之處，特別跟 AVL tree 和 red-black tree 相比，並予以實作

  考慮到循序輸入和亂序輸入的狀況

• 設計效能評比程式，探討上述程式碼和 Linux 核心 red-black tree 效能落差

紅黑樹

紅黑樹(Red-Black Tree)是一種自平衡樹。它的新增、移除、搜尋的時間複雜度均維持在

$O(logn)$ 。AVL Tree 在每次插入和刪除都會調整樹高以致平衡。而紅黑樹則是當違反以下規則時才會調整樹的結構:

1. 每個節點一定是紅色或是黑色
2. 根節點 (root node) 一定要是黑色的
3. 葉節點 (leaf node) 一定要是 NULL (黑色)
4. 紅色節點連接到的節點一定要是黑色的
5. 從根節點到葉節點每個路徑的黑色節點數量要一樣

而紅黑樹的平衡是可以分成以下三種狀況，參考資料結構與演算法：Red Black Tree 紅黑樹 part 1:

1. uncle 節點是紅色

2. uncle 節點是黑色， 插入節點是反向

3. uncle 節點是黑色， 插入節點是順向

每個情況的解法詳細可以看我錄製的說明:

XTree

XTree 是個嘗試兼具 AVL tree 和 red-black tree 部分特性的 binary search tree (BST) 實作，著重在快速的 insert/delete 操作且合理的 lookup 速度。

XTree 的平衡機制與 AVL tree 相似，利用 hint 來評估是否需要做 rotate，然而 xt_node 的 hint 的計算方式是，其左右節點 hint 最大值 +1 並且只有在被呼叫到 xt_update 時才會被更新。

• 結構體

struct xt_tree {
    struct xt_node *root;
    cmp_t *cmp;
    struct xt_node *(*create_node)(void *key);
    void (*destroy_node)(struct xt_node *n);
};

• __xt_destroy: 刪除整個 XTree。使用遞迴的方式刪除左右兩邊的分支。
• xt_first: 獲得子樹依中序(inorder)走訪的第一個節點
• xt_last: 獲得子樹依中序(inorder)走訪的最後一個節點
• xt_rotate_left: 對節點做左旋轉，如下圖所示:

• xt_rotate_right: 對節點做右旋轉，如下圖所示:

• xt_update: 是 XTree 自平衡的機制。

首先會利用 xt_balance 函式取得左子樹的 hint 減去右子數 hint 的差值。

1. 若差值 b 小於 -1 ，代表 XTree 右偏，因此做左旋轉。
2. 若差值 b 大於 1 ，代表 XTree 左偏，因此做右旋轉。

最後去檢查此節點的 hint 是否有改變到，以及此節點是否變成葉節點。若有一成立，則再對此節點的親代節點做 xt_update 。

• xt_insert: 插入一個值為 key 的節點。

先利用 __xt_find 函式尋找節點應該在的位置。其中 p 是節點指標、 d 則是指示左右位置。接著利用 p 和 d 進行插入。最後用 xt_update 去更新 XTree 。

• xt_remove: 刪除值為 key 的節點。

先用 xt_find 找到要刪除的節點。接著用 __xt_remove 來刪除節點：

1. 如果要刪除的節點有右子節點，就從要被刪除的節點的右子節點中的子樹的中序走訪的第一個節點來取代要被刪除的節點。

if (xt_right(del)) {
    struct xt_node *least = xt_first(xt_right(del));
    if (del == *root)
        *root = least;

    xt_replace_right(del, least);
    xt_update(root, least);
    return;
}

2. 若第一點不成立，則就從要被刪除的節點的左子節點中的子樹的中序走訪的最後一個節點來取代要被刪除的節點。

if (xt_left(del)) {
    struct xt_node *most = xt_last(xt_left(del));
    if (del == *root)
        *root = most;

    xt_replace_left(del, most);
    xt_update(root, xt_left(most));
    return;
}

3. 檢查節點是否為根節點，若是，則直接將根節點設為 0 。
4. 若上面皆不成立，代表要被刪除的是葉節點，就直接將其刪除後，對其親代節點做 xt_update 。

/* empty node */
struct xt_node *parent = xt_parent(del);

if (xt_left(parent) == del)
    xt_left(parent) = 0;
else
    xt_right(parent) = 0;

xt_update(root, parent);

main 函式裡的 argc 和 argv

常常在主函式中看到有兩個參數 argc 和 argv ，若需要在執行時的命令列將數值輸入進程式中時使用。

以本題為例，主程式寫成以下樣子：

int main(int argc, char *argv[])

編譯後，可以用以下命令執行:

$ ./build/treeint 1000 0

可以看到上面命令是由三個部份組成，./build/treeint 為執行檔名稱， 1000 和 0 是兩個引入參數。因此可以得到 argc 就是 3 ，而 *argv[] 就是存放者三個參數字串的陣列。

對此就可以對參數做操作，本題就是將第一個參數（不包含執行名稱）作為 XTree 的大小，將第二個參數作為隨機種子。

if (!sscanf(argv[1], "%ld", &tree_size))
...
if (!sscanf(argv[2], "%ld", &seed))

實測效能

這邊將每次的插入、尋找、刪除時間都做平均。並去計算 XTree 的高度。

• 循序輸入

$ ./build/treeint 1000 0
Average insertion time of XTree: 195.669000
Average find time of XTree: 116.198000
Average remove time of XTree: 144.050000
The height of XTree: 9

• 亂序輸入

 ./build/treeint 1000 3
Average insertion time of XTree: 225.347000
Average find time of XTree: 129.095000
Average remove time of XTree: 163.543000
The height of XTree: 10