2024q1 Homework5 (assessment)

contributed by < Lccgth >

紀錄閱讀〈因為自動飲料機而延畢的那一年〉的啟發

文章中的疑惑

O2O 模式是甚麼

O2O means "Online To Offline" but also "Offline to Online", indicating the two-way flow between the online and the physical world, especially retail and ecommerce, but also between brand marketing and shopper or point-of-sale marketing efforts to influence purchase decisions.

文章中提到作者和飲料店合作，讓消費者可以在網站上調製屬於自己的飲料，並且到實體店購買，這就是一種透過線上行銷手法帶動實體通路的消費( Online To Offline )模式，

JavaScript 中 Promise 如何解決 Callback Hell 問題

Callback Hell 是指回呼函式被嵌在其他回呼函式中，當深度達到一定程度時會造成程式碼難以理解和維護，而 Promise 是一種處理非同步流程的 API，它的主要特色之一為通過 .then() 處理每個非同步操作的結果，可以將巢狀結構轉化為垂直的鏈式結構。

function addTea() {
    return new Promise((resolve, reject) => {
        console.log("加茶中...");
        
        if (Math.random() > 0.1) { // 判斷是否成功的條件
            console.log("茶已添加");
            resolve();
        } else {
            reject("添加茶失敗");
        }
    });
}

function addSugar() {
    // ...
}

function addIce() {
    // ...
}

detectCup()
  .then(addTea)
  .then(addSugar)
  .then(addIce)
  .then(() => {
    console.log('飲料制作完成！');
  })
  .catch((error) => {
    console.error('發生錯誤:', error);
  });

Singleton Pattern 是甚麼

保證一個類別只會產生一個實例，而且有統一方法存取該實例，廣泛應用在需要控制資源的存取或是共享資源的情形中，以下為參考 單例模式 Singleton 的一種簡單實現方法。

public class Singleton {
    // 一開始就建立物件
    private static Singleton instance  = new Singleton();

    // 把 constructor 設為 private，其他物件就沒辦法直接用 new 來獲取此類別的實例
    private Singleton(){}

    // 提供方法讓其他程式調用這個類別
    public static Singleton getInstance(){
        return instance;
    }
}

紀錄啟發

在閱讀完這一系列的文章後，我身受啟發並在很多地方感同身受，首先作者有了一個目標 「作一台自動飲料機」，乍聽之下聽起來很簡單，但其中的細節卻數不勝數，從最後的成品應具備甚麼樣的功能、目標受眾的族群、到能否從中獲得利益等，這些都是必須在開始動工前先決定好的大方向，而這也跟軟體設計的流程有著類似之處。

軟體開發初期，首要的任務為清楚了解需求，接著是規劃系統的架構，決定要用什麼工具來開發，然後才是開始根據規範來實作程式碼，最終進行測試及維護。我曾參與實驗室的一個關於數位孿生的機器人計畫，團隊成員會在初期先開會決定軟體的部分需要實現哪些功能，以及必須準備哪些相應的硬體設備，隨後的會議則是根據功能的實作方式與成果來進行討論和調整，這樣的經驗讓我深刻體會到從構想到成品的每一步都不可忽視，每一個細節都精心計畫與執行。

此外我也體會到人脈的重要性，一個人能力畢竟有限，在自動飲料機的開發過程中，各方面的工作需要不同專長的人協作完成，例如有人負責網頁的架設、有人負責電路的設計、有人負責機械架構的開發，同時，還有許多朋友和老師給予的鼓勵及建議，每當其中任何一環遇到問題，都有可能影響到整個產品的成效，因此，學會如何根據每個人的專長來分工合作是非常關鍵的。在這樣的合作過程中，我們不僅可以學習到以前不熟悉的知識，還能夠獲得他人的獨特見解，這不僅豐富了我們的經驗，也讓我們在專業技能上得到提升。

期末專題

位元操作

在上這堂課之前，我在實作函式時所需要的計算大多使用了如 *、/、% 等運算子，而未曾考慮運用位元操作來提升效率。因此，我希望能透過此期末專題來深入了解位元操作，期望未來能將這些技巧應用於運算相關的程式碼中，從而提升效率。此外，鑑於 Linux kernel 中多處使用位元操作來實作各種函式，我相信這個專題也將加速我對 Linux 原始碼的理解。

透過 Netfilter 自動過濾廣告

隨著網路技術的發展，網頁廣告變得無所不在，尤其是在觀看影片或瀏覽網站時，不相關的廣告常常彈出，令人感到煩擾。目前市面上有許多瀏覽器擴充功能可以對付這些廣告，例如著名的 AdBlock 和 Fadblock，但這些工具大多僅限於 Chrome 瀏覽器。因此，我想了解如何在更底層，也就是操作系統核心層級過濾網路廣告，以便所有應用程式都能從中受益。

研讀教材

心得

從〈軟體缺失導致的危害〉中我體會到身為資訊系的學生，程式碼的嚴謹性遠比我先前想像的要重要。從文中我了解到，軟體相關的錯誤有時會導致嚴重的人員傷亡，這凸顯了軟體工程的重要性。隨著科技的進步，越來越多的設備和系統將依賴於程式的控制，因此作為開發者我們必須保持高度的警覺。在撰寫程式碼時，我們應該以一種負責任的心態來面對，即這些程式碼將被數千甚至數萬的使用者使用。因此，我們需要仔細地關注每一個細節，確保軟體的安全性和穩定性，而不僅僅是讓程式能夠運行。

疑惑

問題 1

在〈你所不知道的 C 語言: linked list 和非連續記憶體〉中有一段敘述為 node = L1->val < L2->val? &L1: &L2 不能寫成 node = &(L1->val < L2->val? L1: L2)，不懂後者在 C 中不成立的原因?

struct ListNode *mergeTwoLists(struct ListNode *L1, struct ListNode *L2) {
    struct ListNode *head = NULL, **ptr = &head, **node;

    for (node = NULL; L1 && L2; *node = (*node)->next) {
        node = (L1->val < L2->val) ? &L1: &L2;
        *ptr = *node;
        ptr = &(*ptr)->next;
    }
    *ptr = (struct ListNode *)((uintptr_t) L1 | (uintptr_t) L2);
    return head;
}

根據 c99 規格書 6.5.15 的敘述:

If an attempt is made to modify the result of a conditional operator or to access it after the next sequence point,the behavior is undefined.
A conditional expression does not yield an lvalue.

lvalue 指向記憶體的一固定位置，而 rvalue 為暫時的值，沒有固定的記憶體位置，因此由三元運算子產生的結果不能被取址，會導致編譯錯誤或未定義行為。

問題 2

在 〈Linux 核心的紅黑樹〉中提到 AVL 的樹高約為

$1.44\times log(n+2)$，而 rbtree 則為

$2\times log(n+1)$，不懂這兩個數值是如何得出的?

AVL tree

首先先證明在高度為

$h$ 的 AVL tree 中，最少有

$F_h+1$ 個節點，而

$F_h$ 為費波那契數列。

費波那契數列的定義為

$F_0=0$，且

$F_1=1$，而對於所有

$n>=2$，

$F_n=F_{n-1}+F_{n-2}$ 。

根據數學歸納法，可觀察基本情況為當樹高

$h=1$ 時節點至少為

$1$ 個而當樹高

$h=2$ 時節點至少為

$2$ 個，接著假設當

$h<k$ 時的所有

$h$ 皆成立。當

$h=k$ 時所有節點的總數為

$1+F_{k-1}+F_{k-2}$，再根據費波那契數列的特性

$F_k=F_{k-1}+F_{k-2}$，可得證在高度為

$h$ 的 AVL tree 中，最少有

$F_h+1$ 個節點。

接著根據數學公式推論 AVL 的樹高約為

$1.44\times log(n+2)$。

根據費波那契數列的定義

$F_n=\frac{{\phi }^n}{\sqrt{5}}$，當有一棵包含 n 個節點的 AVL tree，可以反推其樹高。

$$n\cong F_h-1$$

$$n\cong \frac{{\phi }^h}{\sqrt{5}}-1$$

$${\phi }^h\cong (n+1)\sqrt{5}$$

$$h\cong lo{g}_{\phi }[(n+1)\sqrt{5}]\cong lo{g}_{\phi }(n+2)\cong 1.44\times lo{g}_2(n+2)$$

$$lo{g}_2(\phi )\cong lo{g}_2(1.618)\cong 0.69$$

$$\frac{1}{0.69}\cong 1.44$$

rbtree

假設 rbtree 中根節點到葉節點經過的黑色節點數量為

$bh(T)$，當樹中沒有紅色節點時根據樹高與節點數的關係

$$n=2^{bh(T)}-1$$
當樹中有紅色節點時

$$n\ge 2^{bh(T)}-1$$
根據紅黑樹不會有連續紅色節點的特性，可得到

$bh(T)$ 至少為樹高的一半

$$bh(T)\ge \frac{h}{2}$$

$$n\ge 2^{bh(T)}-1\ge 2^{\frac{h}{2}}-1$$

$$2^{\frac{h}{2}}\le n+1$$

$$\frac{h}{2}\le log(n+1)$$

$$h\le 2\times log(n+1)$$

問題 3

在 〈Linux 核心的紅黑樹〉中提到 commit b0687c 在 rb_set_parent_color 中以 + 取代 | 可以在 x86 善用 LEA 指令減少總指令數，不懂其中的運作原理。

在閱讀此 commit 的敘述時，對其中的一段敘述感到疑惑

The benefit to x86 is it change the codegen for setting a node to block from mov %r0, %r1; or $RB_BLACK, %r1 to lea RB_BLACK(%r0), %r1 which
saves an instructions.

此處的 block 看起來是筆誤，替換成 black 語意較通順。

為何指令等效

以 lea RB_BLACK(%r0), %r1 為例，此指令的功能為將 RB_BLACK 和 %r0 作相加後將結果存在 %r1 中，而原本的 mov %r0, %r1 是將 %r0 的值複製到 %r1，or $RB_BLACK, %r1 是讓 RB_BLACK 和 %r1 作 | 並將結果存在 %r1 中，而因為 RB_BLACK 的值固定為 1，且 %r0 最後一個位元固定為 0，所以在此情況下 + 和 | 等效。

---

clone 系統呼叫的作用？

https://man7.org/linux/man-pages/man2/clone.2.html

clone 系統呼叫是用來創立一個子行程，比起 fork 提供了更精確的控制，例如控制行程間是否要共享虛擬位址( CLONE_VM )、檔案描述表( CLONE_FILES )、信號處理表( CLONE_SIGHAND )等等，還可以將子行程放在不同的命名空間 ( CLONE_NEWNS )。

	fork	clone	clone3
功能	創立行程	創立行程或執行緒	創立行程或執行緒
共享性	不共享	可選擇是否要共享	可選擇是否要共享
參數	無	多個控制參數	結構化參數(clone_args)
堆疊	自動指派	需指定堆疊位置	需指定堆疊位置
回傳值	子行程為0，父行程為子行程 PID	子行程為0，父行程為子行程 PID	子行程為0，父行程為子行程 PID

---

TODO: 閱讀 https://hackmd.io/@sysprog/ByS9w_lHh (含錄影) 並紀錄你的疑惑

問題 1

在閱讀 __ffs 程式碼時，不懂為什麼 num 要宣告成 int，此函式的回傳型態為 unsigned long，在回傳時需要進行一次轉型，若一開始就將 num 宣告成 unsigned long，是否可省略此轉型成本?

static __always_inline unsigned long __ffs(unsigned long word)
{
	int num = 0;

       // 計算 ffs 過程
    
	return num;
}

此問題在 commit 3eb3c33 被修正了，而此 commit 中關於 bitops 還有另一項貢獻：

bitops: Change function return types from long to int

原本 bitops 相關的程式碼回傳型態包含 unsigned long 和 unsigned int，存在不一致的問題，在這次的 commit 中將回傳型態統一為 unsigned int，雖然這兩種資料型態在不同處理器架構下位元數量可能會不同，例如在 64 位元的系統中，unsigned int 為 32 位元，而 unsigned long 為 64 位元，但是這些更動的函式回傳值已經被限定在 0 ～ 63 之間，所以可以統一將回傳型態更改為 unsigned int。

問題 2

__gen_ffs 的參數中有一項為 sizeof(in_type) * 8，而這也是和 gen_ffs 唯一的參數不同之處，但這裡的計算是仰賴於其中一個參數 in_type， 若是在 __gen_ffs 的巨集展開函式中宣告 bit_len，是否就能省略此參數，提高 __gen_ffs 和 gen__ffs 的一致性?

原本程式碼

#define __gen_ffs(in_type, out_type, start_from, fn_name, bit_len) \
    static __always_inline out_type fn_name(in_type x)             \
    {                                                              \
                                                                   \
        shift = bit_len >> 1;                                      \
        mask = (in_type) ~(_ONES << (bit_len >> 1));               \
                                                                   \
        return bits;                                               \
    }

#define gen_fls(in_type, out_type, start_from, fn_name) \
    __gen_fls(in_type, out_type, start_from, fn_name, (sizeof(in_type) * 8))

嘗試改進

#define __gen_ffs(in_type, out_type, start_from, fn_name)          \
    static __always_inline out_type fn_name(in_type x)             \
    {                                                              \
        const int bit_len = sizeof(in_type) * 8;                   \
        shift = bit_len >> 1;                                      \
        mask = (in_type) ~(_ONES << (bit_len >> 1));               \
                                                                   \
        return bits;                                               \
    }

#define gen_fls(in_type, out_type, start_from, fn_name) \
    __gen_fls(in_type, out_type, start_from, fn_name)

作者在提問清單中回答了相似問題，他說明他的設計想法主要有兩個優點:

1. 呼叫者不需要輸入 bit_len 參數，因為這可以從輸入型別直接推得。
2. bit_len 這項輸入在 __gen_fls 中用上兩次。

嘗試改進的程式碼保留以上兩點，且讓呼叫時更加簡潔，接著需做實驗驗證是否會有效能上的改變。

原本程式碼

	cache-misses	cache-references	instructions	cycles
10k	1,4224	6,2514	1391,9494	1165,3715
20k	2,4392	6,4771	2679,5783	2230,2978
30k	2,0021	7,6756	3972,1235	3295,0284
40k	2,1108	6,7162	5259,1896	4355,7414
50k	2,0070	6,8844	6548,1820	5395,5590

嘗試改進後

	cache-misses	cache-references	instructions	cycles
10k	1,4291	6,2514	1421,1424	1150,7252
20k	2,1252	7,4009	2743,9896	2205,4977
30k	1,9854	7,0670	4061,2386	3243,0903
40k	1,5812	6,7534	5377,3521	4285,8376
50k	1,9591	7,1855	6698,0794	5336,1987

根據實驗結果，雖然效能上沒有顯著的提昇，但改進的程式碼能在不影響效能的情況下讓呼叫時能更加的方便與簡潔。

問題 3

__ffs 程式碼的功能為尋找傳入值的第一個位元，在此情形下有固定的回傳值範圍，而這裡的回傳型態為 unsigned long，如果使用 unsigned char 等占用較小記憶體空間的型態，是否可以節省記憶體空間?

// __ffs.h
/**
 * __ffs - find first bit in word.
 * @word: The word to search
 *
 */
static __always_inline unsigned long __ffs(unsigned long word)

目前想法為是否是因為普遍做 4 byte alignment，所以就算其真正使用的大小只有 1 byte，電腦也會給他 4 byte 的空間?

---

重做第 3 週測驗題

測驗 1

解釋測驗程式碼

設所求的平方根答案為 N，透過 digit-by-digit calculation 將 N 拆成 2 的冪相加。

$$N^2=(000b_0{b}_1...b_{n-1}{b}_n{)}^2$$
其中

$b_0$ 為最高位的 1，可再將 N 表示為

$$N^2=(a_n+a_{n-1}+...a_0{)}^2$$
其中

$a_i$ 為 2 的冪或 0，再依照

$(x+y{)}^2=x^2+2xy+y^2$ 的公式逐層拆解。

$\begin{array}{rl}{N}^2& =({\displaystyle \sum _{i=0}^n{a}_i{)}^2}\\ & =a_n^2+[2a_n+a_{n-1}]a_{n-1}+...[2{\displaystyle \sum _{i=1}^n{a}_i+a_0]a_0}\\ & =a_n^2+[2P_n+a_{n-1}]a_{n-1}+...+[2P_1+a_0]a_0\\ P_m& =a_n+a_{n-1}+...+a_0\\ & =P_{m+1}+a_m\\ P_m^2& =P_{m+1}^2+2a_m{P}_{m+1}+a_m^2\end{array}$

接著設

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

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

$X_m=N^2-P_m^2=X_{m+1}-Y_m$，再把

$Y_m$ 拆解為

$c_m+d_m$。

$$c_m=2P_{m+1}{2}^m，d_m=(2^m{)}^2$$

$$c_{m-1}=P_m{2}^m=(P_{m+1}+a_m)2^m=P_{m+1}{2}^m+a_m{2}^m={\displaystyle \frac{c_m}{2}}+d_m$$
而

$c_{-1}=P_0{2}^0=N$ 就是所求的 N。

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

此函式中的 b 對應到數學式的

$Y_m$，z 對應到

$c_m$，而 m 對應到

$d_m$。

使用 ffs / fls 取代 __builtin_clz

fls 的功能為從 lsb 開始往高位尋找最後一個為 1 的位元的位置，而在版本三中的 31 - __builtin_clz(x) 等價於 fls(x) - 1，所以可使用 fls 替換掉 __builtin_clz 使程式不依賴 GNU extension。

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

在 Linux 核心原始程式碼找出對整數進行平方根運算的程式碼

在 block/blk-wbt.c 中的 rwb_arm_timer 函式中使用到了 int_sqrt 來計算整數的平方根，註解說明此函式為根據 CoDel 的 bufferd writeback throttling 機制的一部分，會根據 rqd->scale_step 的值來調整當前的 window 大小。

如果 rqd->scale_step 的值大於 0，註解說明應該要基於 fast inverse square root 來縮小 window，但因為執行次數不高，所以沒有做，對這部分感到疑惑。

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

其中 rwb->win_nsec 為預設的 window 大小。

RWB_WINDOW_NSEC = 100 * 1000 * 1000ULL;
rwb->win_nsec = RWB_WINDOW_NSEC;

CoDel 是甚麼

在網路路由中，通過設置緩衝區來克服硬體中的緩衝膨脹，當封包從較快速的網路要傳輸到較慢速的網路時，設計了一個緩衝區讓快速網路能將封包先儲存在其中，而慢速網路就能以自己的速度讀取封包，當兩者速度差距太大時，緩衝區內的封包會越來越多，導致延遲增加。

CoDel 會持續測量緩衝區中的封包延遲時間，若高於設定的延遲時間 (5 毫秒)，則認定發生了緩衝膨脹，開始丟棄封包讓延遲能低於設定的延遲時間。

fast inverse square root 是甚麼

是一種快速計算

$x$ 的

$\frac{1}{\sqrt{x}}$ 的演算法，參考〈從 √2 的存在談開平方根的快速運算〉和 Fast Inverse Square Root — A Quake III Algorithm 嘗試理解其原理。

float InvSqrt(float x)
{  
    float xhalf = 0.5f * x;
    int i = *(int *) &x;
    i = 0x5f3759df - (i >> 1);
    x = *(float *) &i;
    x = x * (1.5f - xhalf * x * x);
    return x;
}

因為除法操作非常耗時，所以採用近似值來計算。

解析程式碼

一開始先將 0.5 倍的 x 的值存在 xhalf 變數中供後續的計算使用。

float xhalf = 0.5f * x;

接著將 x 的型態從 float 轉型為 int，如此以來才方便使用位元操作來減少計算成本，但這裡特別的是一般在轉型時的寫法為 i = (int) x，這兩者的不同之處在於:

float x = 3.33;       // x  = 0 10000000 10101010001111010111000

int i1 = (int) x;     // i1 = 00000000 00000000 00000000 00000011
int i2 = *(int *) &x; // i2 = 01000000 01010101 00011110 10111000

第二種轉型首先會從 &x 中取得 float x 的記憶體位址，接著通過 (int *) 會將 x 的記憶體位址從 float 型態改為 int 型態，雖然記憶體位址本身並沒有更改，但對 C 而言此記憶體位址中的值已經變成 int，於是在最後對此記憶體位址取值時就會直接將裡面的 32 位元讀成 int。

這時會注意到程式碼中出現了一個 magic number(0x5f3759df)，接下來會一步步推導此數值出現的原因。在 IEEE 754 單精度的表示法中 32 個位元被分成了三個部分，以浮點數 3.33 為例:

|   0   | 10000000 | 10101010001111010111000 |
| sign  | exponent |        fraction         |

此值在十進位的表示中，為

$$(1+\frac{frac}{2^{23}})\times 2^{exp-127}$$
此時將此值取

$log$

$$lo{g}_2((1+\frac{frac}{2^{23}})\times 2^{exp-127})=lo{g}_2(1+\frac{frac}{2^{23}})+exp-127$$
因為

$\frac{frac}{2^{23}}$ 的值介於 0 到 1 之間，因此

$$lo{g}_2(1+\frac{frac}{2^{23}})\approx \frac{frac}{2^{23}}+\mu (\mu \approx 0.0450466)$$
把此結果帶回上述算式

$$\begin{array}{rl}\frac{frac}{2^{23}}+\mu +exp-127& =\frac{frac}{2^{23}}+exp+\mu -127\\ & =\frac{frac}{2^{23}}+\frac{2^{23}\times exp}{2^{23}}+\mu -127\\ & =\frac{frac+2^{23}\times exp}{2^{23}}+\mu -127\end{array}$$
有趣的是

$frac+2^{23}\times exp$ 轉成二進位表示時:

frac         = 0 00000000    frac(23 bits)
2^{23} * exp = 0 exp(8 bits) 0000....0
------------------------------------------
sum          = 0 exp(8 bits) frac(23 bits) 

所以如果要對浮點數取

$log$，只需將其除以

$2^{23}$ 後再加上常數

$\mu -127$ 即可。

此時設所求的解為

$\mathrm{\Gamma }$ :

$$\begin{array}{rl}\mathrm{\Gamma }& =\frac{1}{\sqrt{x}}\\ log(\mathrm{\Gamma })& =log(\frac{1}{\sqrt{x}})=-\frac{1}{2}log(x)\\ \frac{1}{2^{23}}(fra{c}_{\mathrm{\Gamma }}+2^{23}\times ex{p}_{\mathrm{\Gamma }})+\mu -127& =-\frac{1}{2}(\frac{1}{2^{23}}(fra{c}_x+2^{23}\times ex{p}_x))+\mu -127)\\ fra{c}_{\mathrm{\Gamma }}+2^{23}\times ex{p}_{\mathrm{\Gamma }}& =-2^{22}(\frac{1}{2^{23}}(fra{c}_x+2^{23}\times ex{p}_x))-2^{22}(\mu -127)-2^{23}(\mu -127)\\ & =2^{22}(127-\mu )+2\times 2^{22}(127-\mu )-\frac{1}{2}(fra{c}_x+2^{23}\times ex{p}_x)\\ & =3\times 2^{22}(127-\mu )-\frac{1}{2}(fra{c}_x+2^{23}\times ex{p}_x)\\ & =0x5f3759df-(x>>1)\end{array}$$
此行程式碼就是依照以上推導而得出的，用於計算

$i$ 的

$\frac{1}{\sqrt{i}}$。

i = 0x5f3759df - (i >> 1);

計算完畢後再將型態轉回 float。

x = *(float *) &i;

接著再利用 Newton's method 來改進近似值，Newton's method 的基本公式為:

$$x_{n+1}=x_n-\frac{f(x_n)}{f^{\prime }(x_n)}$$
這裡所要求的是

$\frac{1}{\sqrt{x}}$，設:

$$y=\frac{1}{\sqrt{x}}$$
希望能找到 y 的值令

$$f(y)=\frac{1}{y^2}-x=0$$
其導數為

$$f^{\prime }(y)=-\frac{2}{y^3}$$
將兩者帶入到 Newtons's method 中

$$\begin{array}{rl}{y}_{n+1}& =y_n-\frac{f(y_n)}{f^{\prime }(y_n)}\\ & =y_n-\frac{\frac{1}{y_n^2}-x}{-\frac{2}{y_n^3}}\\ & =y_n-\frac{1-x{y}_n^2}{y_n^2}\times (-\frac{y_n^3}{2})\\ & =y_n+\frac{y_n-x{y}_n^3}{2}\\ & =y_n+\frac{y_n}{2}-\frac{x{y}_n^3}{2}\\ & =\frac{3y_n}{2}-\frac{x{y}_n^3}{2}\\ & =y_n(\frac{3}{2}-\frac{x{y}_n^2}{2})\end{array}$$
對應到程式碼中的

x = x * (1.5f - xhalf * x * x);

利用 fast inverse square root 改進 linux kernel

使用 fast inverse square root 讓程式碼中不依賴開根號與除法，預期將 div_u64 移除。

rwb->cur_win_nsec = div_u64(rwb->win_nsec << 4,
                int_sqrt((rqd->scale_step + 1) << 8));

首先要小幅修改上述介紹的 fast inverse square root，因為上述程式碼推導的是基於 32 位元，而此處使用到的是基於 64 位元，推導過程與上述相同，只須將 IEEE 754 單精度換成雙精度，這裡就不贅述，直接展示修改好的程式碼。

double InvSqrt64(double x)
{
    double xhalf = 0.5 * x;
    int64_t i = *(int64_t *)&x;  
    i = 0x5fe6eb50c7b537a9 - (i >> 1);
    x = *(double *)&i;
    x = x * (1.5 - xhalf * x * x);
    return x;
}

接著利用此函式對原本的程式碼進行修改。

rwb->cur_win_nsec = rwb->win_nsec << 4;
double scale_factor = InvSqrt64((double)((rqd->scale_step + 1) << 8));
rwb->cur_win_nsec = (u64)(rwb->cur_win_nsec * scale_factor);

做實驗測試效能差別

測試用命令

$ sudo perf stat --repeat 100 -e cache-misses,cache-references,instructions,cycles

原本程式碼

	Result 1	Result 2	Result 3	Result 4
cache-misses	1826	1597	2064	5973
cache-references	6,1509	6,2647	6,1889	6,1523
instrcuciotns	99,2204	99,1723	99,1330	99,3471
cycles	89,4694	87,9199	88,3323	99,6138
time (sec)	0.0014	0.0014	0.0014	0.0016

使用 fast inverse square root 改寫

	Result 1	Result 2	Result 3	Result 4
cache-misses	1441	1913	2637	4867
cache-references	6,1086	6,1883	6,1737	6,1851
instrcuciotns	99,5354	99,4980	99,4249	99,4466
cycles	90,0680	88,5366	89,7526	99,0297
time (sec)	0.0013	0.0011	0.0010	0.0012

從實驗結果可以得知，修改後的程式碼執行速度比原本快了 0.0001 ~ 0.0004 秒，相當於提升了 7% ~ 28%。

測驗 2

解釋測驗程式碼

為了減少運算成本，題目中想利用位元運算來代替 / 和 % 兩種運算元，因為除以 10 相等於乘

${\displaystyle \frac{1}{10}}$，所以利用 2 的冪的倒數和來逼近。

$${\displaystyle \frac{1}{128}}+{\displaystyle \frac{1}{32}}+{\displaystyle \frac{1}{16}}={\displaystyle \frac{13}{128}}\approx {\displaystyle \frac{1}{10}}$$
而

$13=8+4+1=2^3+2^2+2^0$，因此用位元運算可得到

${\displaystyle \frac{13}{8}}n$。

q = ((((n >> 3) + (n >> 1) + n) << 3)  // q = 13n

接著將做位元運算時捨棄的值加回後再右移 7 (除以 128)。

d0 = n & 0b1;
d1 = n & 0b11;
d2 = n & 0b111;
q = ((((n >> 3) + (n >> 1) + n) << 3) + d0 + d1 + d2) >> 7;  // q = 13n / 128

而

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

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

而另一種算法先將 q 逼近到 0.8 in。

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

$$x={\displaystyle \frac{3}{4}}in$$

$$q={\displaystyle \frac{3}{64}}in+{\displaystyle \frac{3}{4}}in={\displaystyle \frac{51}{64}}in\approx 0.8in$$
將 0.8 右移 3 位可得到 0.1 ( div )。

*div = (q >> 3);

再計算 mod = in - div * 10。

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

改進測驗程式碼

在計算

$q=13n$ 時是先計算出

$\frac{13n}{8}$ 後再乘 8。

q = ((((n >> 3) + (n >> 1) + n) << 3

可以修改為直接使用位元操作計算出

$13n=8n+4n+n$。

q = (n << 3) + (n << 2) + n 

而原本程式碼會記錄位元運算時捨棄的值(d0、d1、d2)，但再加完這些值後會做一次右移 7，導致這些位元會再次被捨去，於是這個操作可以被省略。

- unsigned d2, d1, d0, q, r;
+ unsigned q,r;
- d0 = n & 0b1;
- d1 = n & 0b11;
- d2 = n & 0b111;
- q = ((((n >> 3) + (n >> 1) + n) << 3) + d0 + d1 + d2) >> 7;
+ q = ((n << 3) + (n << 2) + n) >> 7;
  r = n - (((q << 2) + q) << 1);

測試結果:

q: 0 r: 0
q: 0 r: 1
q: 0 r: 2
q: 0 r: 3
q: 0 r: 4
q: 0 r: 5
q: 0 r: 6
q: 0 r: 7
q: 0 r: 8
q: 0 r: 9
q: 1 r: 0
q: 1 r: 1
q: 1 r: 2
q: 1 r: 3
q: 1 r: 4
q: 1 r: 5
q: 1 r: 6
q: 1 r: 7
q: 1 r: 8
q: 1 r: 9

結果和原本的程式碼一致。

撰寫不依賴除法指令的 %5 和 %9

%5

2 的冪 mod 5 有以下的特性

$$\begin{array}{r}{2}^{1}mod5=2\\ 2^{2}mod5=4\\ 2^{3}mod5=3\\ 2^{4}mod5=1\\ 2^{5}mod5=2\\ ...\end{array}$$
可觀察到會以四個為一組循環(2、4、3、1)，接著因為

$2^{16}、2^8、2^4$ mod 5 都等於 1，所以對於任意整數 a

$$a\times 2^{16}(mod5)=a\times 1(mod5)=a(mod5)$$
於是我們可以將 32 位元的整數 (x) 拆解為兩個部分，分別是前 16 位元 (a) 和後 16 位元 (b)

$$x=a\times 2^{16}+b$$
根據取模的分配律

$$\begin{array}{rl}x(mod5)& =a\ast 2^{16}(mod5)+b(mod5)\\ & =a(mod5)+b(mod5)\end{array}$$
可根據位元運算來逐步壓縮 x 的數值。

x = (x >> 16) + (x & 0xFFFF);
x = (x >> 8) + (x & 0xFF);
x = (x >> 4) + (x & 0xF);

壓縮完的數值會介在 0 ~ 63 之間，可以直接建一個表紀錄每個值 mod 5 的結果即可。

int table[63] = {0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0,
                 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1,
                 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2,
                 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1};
return table[x];

%9

2 的冪 mod 5 有以下的特性

$$\begin{array}{r}{2}^{1}mod9=2\\ 2^{2}mod9=4\\ 2^{3}mod9=8\\ 2^{4}mod9=7\\ 2^{5}mod9=5\\ 2^{6}mod9=1\\ 2^{7}mod9=2\\ ...\end{array}$$
觀察到

$2^{18}、2^{12}、2^6$ mod 9 都等於 1，於是修改壓縮數值的程式碼。

x = (x >> 18) + (x & 0x3FFFF);
x = (x >> 12) + (x & 0xFFF);
x = (x >> 6) + (x & 0x3F);

接著依照和 %5 相同的方式建表，將壓縮完的值進行查表即可。

測驗 3

解釋測驗程式碼

版本二函式依序判斷目前的 i 是否大於或等於 65536 (

$2^{16}$ )、256 (

$2^8$ )、16 (

$2^4$ )、2 (

$2^1$ )，如果成立的話就將 i 進行右移並且增加 result，如此一來可找到最大的 result 使得

$2^{result}<=i$。

while (i >= 65536) {
    result += 16;
    i >>= 16;
}
while (i >= 256) {
    result += 8;
    i >>= 8;
}
while (i >= 16) {
    result += 4;
    i >>= 4;
}
while (i >= 2) {
    result += 1;
    i >>= 1;
}
return result;

版本三函式使用 GCC 內建的函式 __builtin_clz 計算傳入的無號整數從最高位開始連續 0 的數量，再使用 31 減 __builtin_clz 的結果，為了處理輸入為 0 的情況，會將 __builtin_clz 的輸入參數和 1 做 or。

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

參考 linux/log2.h 其中 31 - __builtin_clz(v) 等價於 fls(v) - 1，fls 會回傳從最低位往最高位找到的最後一個 1。

static __always_inline __attribute__((const))
int __ilog2_u32(u32 n)
{
	return fls(n) - 1;
}

Linux 核心中 log2 相關程式碼

在 linux/include/linux/log2.h 中存在與 log2 相關的函式與巨集，接下來會逐一說明各自的功能。

__ilog2_u32(u32 n) 和 __ilog2_u64(u64 n)

對於輸入的 32 或 64 位元無號整數 n，找出其以 2 為底的對數，使用 fls() 和 fls64() 找出最高位的 1 的位置後減 1 並返回。

is_power_of_2(unsigned long n)

檢查輸入的 n 是否為 2 的冪，判斷 n & (n - 1) 的結果，如果結果為 0，代表 n 是 2 的冪。

在 2 進位表示中，如果 n 是 2 的冪，則只有 1 個位元會是 1，而 n - 1 會將原本最高位的 1 變成 0，最高位前的 0 變為 1，以 n = 64 為例:

n           = 100 0000  (64)
n - 1       = 011 1111  (63)
n & (n - 1) = 000 0000  (0)

__roundup_pow_of_two(unsigned long n)

將輸入的 n 向上取到最接近的 2 的冪。

__rounddown_pow_of_two(unsigned long n)

將輸入的 n 向下取到最接近的 2 的冪。

const_ilog2(n)

此巨集會先判斷 n 編譯時是否為常數，再來從最高位依序往最低位檢查每個位元是否為 1，若找到第一個位元為 1 的位置時就回傳該位置。

ilog2(n)

此巨集會根據 n 的大小來決定要使用 __ilog2_u32(u32 n) 或 __ilog2_u64(u64 n) 來取得 n 的以 2 為底的對數。

roundup_pow_of_two(n)

此巨集會先判斷 n 編譯時是否為常數，並根據結果決定要使用 ilog2(n) 或 __roundup_pow_of_two(unsigned long n) 將 n 向上取到最接近的 2 的冪。

rounddown_pow_of_two(n)

此巨集會先判斷 n 編譯時是否為常數，並根據結果決定要使用 ilog2(n) 或 __rounddown_pow_of_two(unsigned long n) 將 n 向下取到最接近的 2 的冪。

__order_base_2(unsigned long n)

回傳 n 向上取到 2 的冪後為 2 的幾次方。

ob2(0) = 0     //   0 -> 1
ob2(1) = 0     //   1 -> 1
ob2(2) = 1     //   2 -> 2
ob2(3) = 2     //   3 -> 4
ob2(4) = 2     //   4 -> 4
ob2(5) = 3     //   5 -> 8

order_base_2(n)

此巨集會先判斷 n 編譯時是否為常數，並根據結果決定要使用 ilog2(n) 或 __order_base_2(unsigned long n) 回傳 n 向上取到 2 的冪後為 2 的幾次方。

__bits_per(unsigned long n)

回傳 n 需要用幾個位元來表示。

ob2(0) = 1     //     0
ob2(1) = 1     //     1
ob2(2) = 2     //    10
ob2(3) = 2     //    11
ob2(4) = 3     //   100
ob2(5) = 3     //   101

bits_per(n)

此巨集會先判斷 n 編譯時是否為常數，並根據結果決定要使用 ilog2(n) 或 __bits_per(unsigned long n) 回傳 n 需要用幾個位元來表示。

測驗 4

EWMA 是一種用於平滑時間序列數據的統計手法，讓越久以前的資料所佔的權重越低。

解釋測驗程式碼

ewma_init

初始化 EWMA 結構，設定 factor 和 weight，為了使用位元運算而讓這兩個參數的值為 2 的冪，factor 用於在有限的數值範圍內維持盡可能高的精確度，當 factor 越大時可以提高精確度，但同時也會降低 EWMA 能表示的最大值，而 weight 用於控制舊值對於新值的影響程度，當 weight 越大時會使平均值對於新值的影響越小。

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

ewma_add

向 EWMA 中添加一個新值，如果 internal 不為 0，則按照已有的累計平均值計算新的加權平均值，如果 internal 為 0，則直接根據 factor 調整新值。

avg->internal = avg->internal
                    ? (((avg->internal << avg->weight) - avg->internal) +
                       (val << avg->factor)) >> avg->weight
                    : (val << avg->factor);

ewma_read

返回 EWMA 當前的加權平均值，通過對 internal 依 factor 反向縮放得到。

return avg->internal >> avg->factor;

測驗 5

解釋測驗程式碼

設定 r 和 shift 紀錄各步驟的位移量和最終的對數結果，接著將傳入的 x 值減 1，這是為了處理當 x 為 2 的冪時的情況，因為此函式回傳的是

$\lceil {\log }_2(x)\rceil$，如果不先減 1，會導致回傳的結果比預期要多 1。

再檢查 x 是否大於 0xFFFF，如果為真說明 x 的對數最小為 16，將此資訊紀錄在 r 中，並將 x 右移 16 位。

r = (x > 0xFFFF) << 4;                               
x >>= r;

再檢查當前 x 是否大於 0xFF，如果為真說明當前 x 的對數最小為 8，將此資訊紀錄在 shift 中，並將 x 右移 8 位，並將 shift 的值合併到 r 中，累計目前為止的對數結果。

shift = (x > 0xFF) << 3;                             
x >>= shift;
r |= shift;

根據相同的邏輯依序判斷當前 x 是否大於 0xF 和 0x3，設定 shift 和通過 r 累計結果，最後將判斷 x 是否大於 0x1 的步驟合併在最後一行中，並統計以上步驟中的對數結果，再將一開始減的 1 補上。

return (r | shift | x > 1) + 1;

改進程式碼使其得以處理 x = 0 的狀況，並仍是 branchless

新增判斷傳入的 x 值是否為 0，若 x 為 0 則不會讓 x -= 1，使 x 不會產生 underflow。

  uint32_t r, shift;

- x--;
+ x -= !!x

  return (r | shift | (x > 1)) + 1;