Linux Final Project

contributed by <56han>

第 7 週測驗題 - 測驗 1

將 2 個 32 位元寬度的整數合併

Question: 為甚麼 x 和 y 都要加 0L，再進行位移？
Answer: 因為 x 和 y 的位元寬是 32 位元 ( uint32_t ) ，左移 32 位會超出它的位元寬，會出現 warning: left shift count >= width of type。long 是 64 位元，左移 32 位在這種情況下是合法的，不會引發警告。

(-1L >> 32) 計算出的值是0000 0000 0000 0000 0000 0000 0000 0000 1111 1111 1111 1111 1111 1111 1111 1111。

static inline uint64_t bit_compound(uint32_t x, uint32_t y)
{
    return ((0L + x) << 32) | ((y + 0L) & (-1L >> 32));
}

---

SIMD within a register (SWAR)

延伸閱讀: SIMD and SWAR Techniques
SIMD 單指令流多資料流（Single Instruction Multiple Data，SIMD）：是一種採用一個控制器來控制多個處理器，同時對一組資料中的每一個分別執行相同的操作，從而實現空間上的並列性的技術。
SWAR (SIMD within a register)：在一個暫存器上進行單指令多資料流操作。它是一種跨 CPU 暫存器各部分應用 SIMD 平行處理的處理模型，在 SIMD（單指令多資料）並行處理模型中，這些 byte-entities（即小於一個字節的資料單位）經常被用來在一個 CPU 寄存器的不同部分上進行並行計算。

---

SWAR 技術在加法運算中的差異和優勢

一般加法

這樣的加法操作會直接對整個 64 位的數字進行運算，包括進位。

uint64_t z = x + y;

SWAR 加法

SWAR 技術則是通過在寄存器內部同時處理多個較小的資料單元來提高計算效率。例如，對於 64 位寄存器，我們可以將其視為包含了 8 個 8 位的 byte，並對每個 byte 進行獨立的加法運算。使用 SWAR 技術，我們可以在一個指令中並行地對所有 byte 進行加法運算。具體的 SWAR 技術加法公式如下：

z = ((x &~H) + (y &~H)) ^ ((x ^ y) & H)

其中，x &~H 和 y &~H 去掉了每個 byte 的最高位，以避免在進行無進位加法時受到影響，再透過 (x ^ y) & H 找出每個 byte 中的進位情況。
如此一來，在一次操作中就能對所有的 byte 進行加法計算，而不需要逐個處理進位，這大大提高了運算效率。

優勢

1. 高效並行計算：SWAR 技術利用並行計算的特性，使得對每個 byte 的加法運算可以同時進行，大大提高了計算速度。
2. 簡化進位處理：傳統的逐字節加法需要逐個處理進位，而 SWAR 技術通過位操作巧妙地合併了無進位加法結果和進位結果，簡化了進位處理的過程。
3. 減少迭代次數：傳統方法需要多次迭代來處理每個 byte 的加法，而 SWAR 技術在一個指令中完成所有 byte 的加法運算，減少了迭代次數。

---

SWAR Arithmetic

For bytewise (rankwise) math inside a 64-bit register, H is 0x8080808080808080 and L is 0x0101010101010101.

z = x + y

z = ((x &~H) + (y &~H)) ^ ((x ^ y) & H)

詳細步驟：

1. 計算非進位和：
   x &~H 和 y &~H 去除了每個 byte 的最高位（sign bit），留下其餘 7 位來進行無進位加法。

(x &~H) + (y &~H)

2. 處理進位：
   x ^ y 會找出 x 和 y 在每個 byte 中不相同的位元（即可能產生進位的位置）。然後與 H 進行 AND 操作來確保只保留最高位，透過這個操作，我們只保留了那些可能進位的位置。：

(x ^ y) & H

3. 最終結果合併：
   最後一步是用 XOR 把去掉最高位進行加法的結果和進位進行合併。XOR 用來處理進位的效果：如果沒有進位，原結果保持不變；如果有進位，則原結果的該位會翻轉。

((x &~H) + (y &~H)) ^ ((x ^ y) & H)

加法例題 x + y = ?

H = 1000 0000
x = 0000 0001
y = 1111 1110

1. 計算非進位和：
   ~H = 0111 1111
x &~H = 0000 0001
y &~H = 0111 1110
(x &~H) + (y &~H) = 0111 1111

2. 處理進位：
   x^y = 1111 1111
(x^y) & H = 1000 0000

3. 最終結果合併：
   ((x &~H) + (y &~H)) ^ ((x ^ y) & H) = 1111 1111

---

解釋上述程式碼運作原理

本程式預期可以在一段 string 中找到一個字元，並印出該字元後的 string。

int main()
{
    const char str[] = "https://wiki.csie.ncku.edu.tw";
    const char ch = '.';
    
    char *ret = memchr_opt(str, ch, strlen(str));
    printf("String after |%c| is - |%s|\n", ch, ret);
    return 0;
}

預期執行結果：

String after |.| is - |.csie.ncku.edu.tw|

1. 檢查記憶體對齊
   檢查指標 X 是否與 long 的邊界對齊。如果 X 沒有對齊，則返回非零值。

#define UNALIGNED(X) ((uintptr_t) X & (sizeof(long) - 1)) // uintptr_t: unsigned long int pointer

為甚麼要檢查記憶體是否對齊？

1. 提高記憶體存取效率：大多數現代處理器在處理對齊的記憶體存取時效率更高，因為對齊的地址允許處理器在一次操作中存取整個記憶體單元。
2. 減少記憶體存取錯誤：一些處理器在遇到未對齊的記憶體存取時會引發異常或陷阱，這會導致程式崩潰。
3. 優化低層記憶體操作：在進行底層記憶體操作（如 memcpy、memset 等）時，處理對齊問題可以使這些操作更加高效。

2. 定義每次迭代載入的位元組數
   定義每次迭代載入的 long 型別的大小（通常為 4 或 8 Byte）。

#define LBLOCKSIZE (sizeof(long))

3. 定義過小的長度
   檢查長度是否小於一個 long 型別的大小。

#define TOO_SMALL(LEN) ((LEN) < LBLOCKSIZE)

4. 檢測空位元組
   根據不同的 long 型別大小，定義檢測空字節的巨集。

#if LONG_MAX == 2147483647L
#define DETECT_NULL(X) (((X) - (0x01010101)) & ~(X) & (0x80808080)) // 32-bit 版本
#else
#if LONG_MAX == 9223372036854775807L
#define DETECT_NULL(X) \
    (((X) - (0x0101010101010101)) & ~(X) & (0x8080808080808080)) // 64-bit 版本
#else
#error long int is not a 32bit or 64bit type.
#endif
#endif

參考這篇文章。strlen 每次處理 1 個字元，因此需要 n 次比較 O(n)，DETECT(X) 在 32-bit 下一次處理 4 個字元，在 64-bit 下一次處理 8 個字元，因此迴圈次數減少為 n/4（32-bit）或 n/8（64-bit）。理論上 DETECT(X) 版本的時間複雜度仍是 O(n)，但它每次能處理多個字元，因此速度大約是 strlen 的 4 到 8 倍。

5. 檢測特定位元組
   檢測 X 是否包含指定的字節，利用 XOR 操作，目的是將 X 中所有與 mask 相同的位元組變為 0。

#define DETECT_CHAR(X, mask) DETECT_NULL((X) ^ (mask)) //AAAA

6. 處理未對齊部分
   在字節未對齊的情況下逐個字節檢查，直到對齊或長度為零。

const unsigned char *src = (const unsigned char *) str;
unsigned char d = c;

while (UNALIGNED(src)) {
    if (!len--)
        return NULL;
    if (*src == d)
        return (void *) src;
    src++;
}

7. 處理長度對齊的部分
   如果長度足夠大且已對齊，則按 long 型別大小逐塊處理：建立 mask，將 d 字節複製到 long 的每個位元組中。

if (!TOO_SMALL(len)) {
    unsigned long *asrc = (unsigned long *) src;
    unsigned long mask = d << 8 | d;
    mask |= mask << 16;
    for (unsigned int i = 32; i < LBLOCKSIZE * 8; i <<= 1)
        mask |= mask << i; //BBBB

第三行的 mask = 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 0010 1110 0010 1110
第四行的 mask = 0000 0000 0000 0000 0000 0000 0000 0000 0010 1110 0010 1110 0010 1110 0010 1110
第六行的 mask = 0010 1110 0010 1110 0010 1110 0010 1110 0010 1110 0010 1110 0010 1110 0010 1110

8. 快速搜尋
   檢測是否包含目標字節。如果找到，則退出循環。

    while (len >= LBLOCKSIZE) {
        if (DETECT_CHAR(*asrc, mask)) //CCCC
            break;
        asrc++;
        len -= LBLOCKSIZE;
    }

9. 進一步尋找剩餘字節中符合的字元

    src = (unsigned char *) asrc;
}

while (len--) {
    if (*src == d)
        return (void *) src;
    src++;
}

比較 Linux 核心原本的實作和 memchr_opt，設計實驗來觀察隨著字串長度和特定 pattern 變化的效能影響

實驗設計

1. 字串長度變化：測試不同長度的字串，例如 1KB、10KB、100KB、1MB、10MB 等。
2. 特定模式變化：在字串中不同位置放置目標字元，例如開頭、中間、結尾，以及不存在的情況。
3. 多次運行：為了減少隨機誤差，每個測試情況下運行多次，取平均值。
4. 記錄時間：使用高解析度計時器（clock_gettime）記錄每次函數運行的時間。

解釋

measure_memchr 和 measure_memchr_opt：分別測量 memchr 和 memchr_opt 的執行時間。
run_experiment：測試不同長度的字串，並在不同位置測試字母的存在。使用多次運行來獲取平均執行時間。其中 void * memset ( void * str , int c , size_t n ) 函數將指定的值 c 複製到 str 所指向的記憶體區域的前 n 個位元組中，這可以用於將記憶體區塊清除或設定為特定值。

void run_experiment(size_t length) {
    char *str = malloc(length);
    if (!str) {
        perror("malloc");
        exit(EXIT_FAILURE);
    }

    // Initialize the string with some pattern, e.g., 'a'
    memset(str, 'a', length);

    // Positions to test: start, middle, end, not present
    size_t positions[] = {0, length / 2, length - 1, length + 1};

    for (int i = 0; i < 4; i++) {
        size_t pos = positions[i];
        char c = (pos < length) ? 'b' : 'z';
        if (pos < length) str[pos] = c;

        double memchr_time = 0, memchr_opt_time = 0;
        int iterations = 10;

        for (int j = 0; j < iterations; j++) {
            memchr_time += measure_memchr(str, c, length);
            memchr_opt_time += measure_memchr_opt(str, c, length);
        }

        memchr_time /= iterations;
        memchr_opt_time /= iterations;

        printf("Length: %zu, Position: %zu, memchr: %f, memchr_opt: %f\n",
               length, pos, memchr_time, memchr_opt_time);

        if (pos < length) str[pos] = 'a'; // Reset character
    }

    free(str);
}

main：呼叫 run_experiment，測試不同的字串長度。

int main() {
    size_t lengths[] = {1 * KB, 10 * KB, 100 * KB, 1 * MB, 10 * MB};

    for (int i = 0; i < 5; i++) {
        run_experiment(lengths[i]);
    }

    return 0;
}

分析結果

輸出：

Length: 1024, Position: 0, memchr: 0.000000, memchr_opt: 0.000000
Length: 1024, Position: 512, memchr: 0.000000, memchr_opt: 0.000000
Length: 1024, Position: 1023, memchr: 0.000000, memchr_opt: 0.000000
Length: 1024, Position: 1025, memchr: 0.000000, memchr_opt: 0.000000

Length: 10240, Position: 0, memchr: 0.000000, memchr_opt: 0.000000
Length: 10240, Position: 5120, memchr: 0.000000, memchr_opt: 0.000001
Length: 10240, Position: 10239, memchr: 0.000000, memchr_opt: 0.000002
Length: 10240, Position: 10241, memchr: 0.000000, memchr_opt: 0.000002

Length: 102400, Position: 0, memchr: 0.000000, memchr_opt: 0.000000
Length: 102400, Position: 51200, memchr: 0.000000, memchr_opt: 0.000016
Length: 102400, Position: 102399, memchr: 0.000000, memchr_opt: 0.000027
Length: 102400, Position: 102401, memchr: 0.000000, memchr_opt: 0.000049

Length: 1048576, Position: 0, memchr: 0.000000, memchr_opt: 0.000000
Length: 1048576, Position: 524288, memchr: 0.000000, memchr_opt: 0.000124
Length: 1048576, Position: 1048575, memchr: 0.000000, memchr_opt: 0.000277
Length: 1048576, Position: 1048577, memchr: 0.000000, memchr_opt: 0.000242

Length: 10485760, Position: 0, memchr: 0.000000, memchr_opt: 0.000000
Length: 10485760, Position: 5242880, memchr: 0.000000, memchr_opt: 0.001668
Length: 10485760, Position: 10485759, memchr: 0.000000, memchr_opt: 0.003597
Length: 10485760, Position: 10485761, memchr: 0.000000, memchr_opt: 0.002807

從結果中可以看到，大部分情況下，memchr 的運行時間比 memchr_opt 短，這可能是由於以下幾個原因：

1. 函數呼叫開銷：
   memchr_opt 在開始時進行了額外的對齊檢查和預處理步驟，這增加了開銷。
2. 對齊檢查：
   memchr_opt 進行了對齊檢查，這會導致在字串長度較小或未對齊時，額外的開銷可能超過了任何潛在的優化收益。當字串較短時，對齊檢查和其他操作可能反而降低了效能。
3. 分支預測失敗：
   memchr_opt 中的條件分支較多，例如對齊檢查和字串長度檢查，可能會導致分支預測失敗，這會影響效能。
4. 記憶體存取模式：
   當處理較大資料塊時，memchr_opt 利用了 SIMD 技術進行批量處理，這可能在更大資料塊下顯示出優勢，但在較小資料塊下，這些優化的開銷可能反而導致性能下降。

研讀 2022 年學生報告，在 Linux 核心原始程式碼找出 x86_64 或 arm64 對應的最佳化實作，跟上述程式碼比較，並嘗試舉出其中的策略和分析

Linux 核心原始程式碼找出 x86 對應的最佳化實作：arch/x86/lib/string_32.c

void *memchr(const void *cs, int c, size_t count)
{
    int d0;
    void *res;
    if (!count)
        return NULL;
    asm volatile("repne\n\t"
        "scasb\n\t"
        "je 1f\n\t"
        "movl $1,%0\n"
        "1:\tdecl %0"
        : "=D" (res), "=&c" (d0)
        : "a" (c), "0" (cs), "1" (count)
        : "memory");
    return res;
}

x86 最佳化實作策略

1. 利用 x86 的指令集特性：利用 x86 指令集中的 repne scasb 指令進行快速字節掃描。
2. Inline assembly 優化：透過使用 inline assembly，程式碼量減少且更加高效。

memchr_opt 實作策略

1. 記憶體對齊檢查：檢查記憶體是否對齊，對齊後進行後續處理。
2. 多字節處理：使用多字節（例如 unsigned long）進行處理，以加快速度。
3. 比對模式設置：設置比對掩碼，使用 DETECT_CHAR 巨集進行高效比對。
4. 回退到單字節處理：在多字節處理失敗時，回退到單字節處理。

整體分析

x86 最佳化實作：利用了 x86 指令集的特性，通過 inline assembly 直接操作 CPU 指令來提高效能，適合 x86 平台，簡單高效。
memchr_opt 實作：使用通用的 C 語言進行優化，具有跨平台的可移植性，但在特定平台上的性能可能不如專門針對該平台進行優化的實作。

第 7 週測驗題 - 測驗 2

若干處理器架構無法有效地處理除法，且我們想避免除法造成的例外狀況 (如 Division by zero)，於是我們考慮以下實作。

解釋上述程式碼運作原理

1. likely 巨集
   在許多原始碼 linux 核心都可以看到這個巨集。使用 !! 將 x 轉換為布林值，並使用 GCC 的內建函數 __builtin_expect，告訴編譯器 x 的值很可能為真，以便進行更好的分支預測和優化。

#define likely(x) __builtin_expect(!!(x), 1)

為什麼要使用 __builtin_expect 函數？它究竟有什麼作用呢？

• 現在處理器都是管線的，系統可以提前取多條指令進行並行處理，但遇到跳轉時，則需要重新取指令，跳轉指令打亂了 CPU 管線。因此，跳轉次數少的程式擁有較高的執行效率。
• 在 C 語言程式設計時，會不可避免地使用 if-else 分支語句，if-else 編譯後，一個分支的編譯程式碼緊接著前面的程式碼，而另一個分支的編譯程式碼需要使用 JMP 指令才能存取。很明顯透過 JMP 存取需要更多的時間，在複雜的程式中，有很多的 if-else 句型，又或者是一個有 if-else 句型的函數庫函數，每秒被呼叫幾萬次，通常程式設計師在分支預測方面做得很糟，編譯器又不能精準的預測每一個分支，這時 JMP 產生的時間浪費就會很大。
• 因此，引入 __builtin_expect 函數來增加條件分支預測的準確性，CPU 會提前裝載後面的指令，遇到條件轉移指令時會提前預測並裝載某個分支的指令。編譯器會產生對應的代碼來優化 CPU 執行效率。

2. 檢查一個數是否為 2 的冪
   如果 n 是 2 的冪，則它的二進位表示中只有一個位為 1，n & (n - 1) 應該等於 0。

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

3. 計算 32 / 64 位元整數 x 的二的對數
   使用內建函式 __builtin_clz 計算前導零的數量，再根據位元數進行轉換。

static inline int __ilog2_u32(uint32_t x)
{
    return x ? sizeof(x) * 8 - __builtin_clz(x) - 1 : 0;
}

static inline int __ilog2_u64(uint64_t x)
{
    uint32_t h = x >> 32;
    if (h)
        return sizeof(h) * 8 - __builtin_clz(h) + 31;
    return x ? sizeof(x) * 8 - __builtin_clz(x) - 1 : 0;
}

__builtin_clz 這個函數作用是傳回輸入數二進位表示：從最高位元開始 ( 左起 ) 連續的 0 的個數；如果傳入 0 則行為未定義。

int __builtin_clz ( unsigned  int x) 
int __builtin_clzl ( unsigned  long ) 
int __builtin_clzll ( unsigned  long  long ) 

三個不同的函數分別用於 unsigned int 、 unsigned long 以及 unsigned long long。

來源：Tom's develop Blog

4. 計算任意整數 n 的二的對數
   其中 __builtin_constant_p 是編譯器 GCC 內建函數，判斷 n 在 compiler time 是否為常數，是的話就 return 1 並使用 (n) < 2 ? 0 : 63 - __builtin_clzll(n) ，否則 return 0 並使用 (sizeof(n) <= 4) ? __ilog2_u32(n) : __ilog2_u64(n) ，但根據 GCC 提供的文件，這不代表 n 不是常數，只是 GCC 無法證明它是在活躍的優化選項集限制內的常數。

The function returns the integer 1 if the argument is known to be a compile-time constant and 0 if it is not known to be a compile-time constant. A return of 0 does not indicate that the value is not a constant, but merely that GCC cannot prove it is a constant within the constraints of the active set of optimization options.

#define ilog2(n)                                                       \
    (__builtin_constant_p(n) ? ((n) < 2 ? 0 : 63 - __builtin_clzll(n)) \
     : (sizeof(n) <= 4)      ? __ilog2_u32(n)                          \
                             : __ilog2_u64(n))

在活躍的優化選項集限制內的常數 ( a constant within the constraints of the active set of optimization options ) 是甚麼？

• 優化選項的影響
  不同的優化選項會影響編譯器如何判斷一個值是否可以被視為常數。例如，某些優化選項可能允許編譯器進行更積極的 constant folding，這是一種在編譯期間計算表達式結果，並用該結果替換表達式的技術。

5. 實現 64 位整數除以 32 位整數的操作

This approach multiplies by the reciprocal of b: to divide n by b, we multiply n by ( p / b ) and then divide by p.

${\displaystyle \frac{n}{b}}$ = n

$\times$

${\displaystyle \frac{p}{b}}$

$\times$

${\displaystyle \frac{1}{p}}$

因為乘法運算比除法運算在 CPU 上執行更快，它需要更少的 CPU 週期和更少的指令數量。透過事先計算除數的倒數，使得除法轉換為乘法，從而提高運算效率。

The efficiency comes from the compiler's ability to optimize much of this code during compile time through constant propagation, leaving behind only a few multiplication operations. This is why we use this large macro, as a static inline function does not consistently achieve the same optimization.

編譯器在編譯期間能夠透過 constant propagation 來優化這段代碼，最終只留下少量的乘法操作，因此提高運算效率。constant propagation：一種編譯器優化技術，編譯器會分析程式碼中的常數，並將這些常數值直接替換到程式碼中的變數，從而減少運行時的計算量。

為甚麼巨集可以達到更好的 constant propagation 效果？而 inline function 不行？

主要原因是巨集在預處理階段進行文本替換，使得編譯器能夠更早、更靈活地對展開的程式碼進行優化。而 inline function 仍然保留了函數的特性和約束，優化效果可能受到限制。因此，在需要充分利用 constant propagation 進行優化的情況下，巨集通常比 inline function 更有效。

---

計算 ___p

___p ： 接近 ___b 的 2 的冪。

uint32_t ___p = 1 << ilog2(___b);

---

計算 ___m，使其接近 ((p << 64) + b - 1) / b 的值

第一步

___m = (~0ULL / ___b) * ___p

• ~0ULL 是
  $2^{64}$ -1
• (
  $2^{64}$ -1) ÷ b 是對
  $2^{64}$ ÷ b 的近似
• 乘以 p 獲得 (
  $2^{64}$-1) × p ÷ b

第二步

___m += (((~0ULL % ___b + 1) * ___p) + ___b - 1) / ___b;

這是一個修正項，確保結果精確：

• (~0ULL % b + 1) 是 (
  $2^{64}$ -1) 除以 b 的餘數加 1
• 乘以 p 後除以 b 是為了補上第一步近似造成的誤差
• (b - 1) 是為了向上取整

兩步合起來，正好等價於 ((p << 64) + b - 1) / b，這個表達式計算了

$\lceil$ (

$2^{64}$ × p) ÷ b

$\rceil$。

---

___x 是能被 ___b 整除的最大數減 1

/* one less than the dividend with highest result */                
___x = ~0ULL / ___b * ___b - 1;  

目的：在進行大數除法時，可以用來檢查和處理極限情況。

---

測試 m 的精確度

• 計算採用的是 res = (___m * ___x) / (2^64 * ___p) 的方式。
• 由於需要 128 位乘法，代碼將 64 位乘法拆分成 4 個 32 位乘法並正確處理進位。
• ___res 之後做為判斷 bias 之標準，如果 ___res = ___x / ___b ，則 m 足夠精確；否則需要 bias 修正。

/* test our ___m with res = m * x / (p << 64) */
___res = ((___m & 0xffffffff) * (___x & 0xffffffff)) >> 32;
___t = ___res += (___m & 0xffffffff) * (___x >> 32);
___res += (___x & 0xffffffff) * (___m >> 32);
___t = (___res < ___t) ? (1ULL << 32) : 0;
___res = (___res >> 32) + ___t;
___res += (___m >> 32) * (___x >> 32);
___res /= ___p;

---

判斷是否需要 bias - 特殊情況的處理

if (~0ULL % (___b / (___b & -___b)) == 0) {
    ___n /= (___b & -___b);
    ___m = ~0ULL / (___b / (___b & -___b));
    ___p = 1;
    ___bias = true;
}

這是檢查除數 ___b 是否符合特殊模式：檢查

$2^{64}$ -1 是否能被 b 的奇數部分整除。

• (___b & -___b) 提取 ___b 的最低位 1 (即 2 的冪次部分)
• ___b / (___b & -___b) 獲得 ___b 的奇數部分
• 檢查 ~0ULL % 奇數部分 == 0，判斷是否可以簡化

若符合，代碼執行特殊處理：

發現了一個數學特性：除數 b 可以分解為 b = 2^k × q，且 q 是一個能夠整除 2^64-1 的奇數。

• 先除掉 2 的冪次部分：___n /= (___b & -___b)

對於 2^k 部分，直接使用位移操作：n /= (b & -b) 等同於 n >>= k

• 僅針對奇數部分計算反元素：___m = ~0ULL / (奇數部分)

• 當奇數 q 能整除
  $2^{64}$ -1 時，有性質：
  $2^{64}$ % q = 1
• 這意味著我們可以使用簡單的乘法逆元：m = ~0ULL / q

• 設置 ___p = 1 和 ___bias = true

---

判斷是否需要 bias - 偏差修正

else if (___res != ___x / ___b)，說明表示 m 不夠精確，需要 bias。也就是說，透過乘法和位移運算來近似除法的過程中，出現了誤差（由於有限精度和位截斷引起）。

It is necessary to introduce a bias to counteract errors caused by bit truncation.

在數字運算中，特別是當我們處理浮點數或大數進行除法、乘法等運算時，計算結果可能超出變數類型所能表示的範圍。為了解決這個問題，我們通常會截斷多餘的位元，這個過程稱為位元截斷（bit truncation）。在這段程式碼中，位元截斷可能會在進行乘法和除法運算時發生。當結果超過了變數所能表示的位數時，多餘的位數會被丟棄，這可能導致數值誤差。為了減少或抵消這些誤差，需要引入一個偏差（Bias）。這個偏差是一個額外的數值，用來調整計算結果，使其更接近真實值。

Without this bias, an extra bit would be required to accurately represent the variable 'm', which would exceed the capacity of a 64-bit variable. To manage this, the solution involves setting 'm' equal to 'p' divided by 'b', and 'n' divided by 'b' as equal to (n * m + m) / p. This approach is taken because avoiding the bias is not feasible due to the limitations imposed by bit truncation errors and the 64-bit variable size constraint.

位元容量限制：如果沒有引入偏差，我們需要額外的位元來表示計算過程中的變數 m。這是因為在計算過程中，可能會出現結果超過 64 位元變數所能表示的範圍。
解決方案：重新計算 m，使用標準公式：

m = (~0ULL / b) * p + ((~0ULL % b + 1) * p) / b

Question: 目前不理解為甚麼在變數 m 、 p 、 b 沒有改變的情況下，這裡需要再計算一次變數 m ？ 使m = (p << 64) / b 。
原本計算過一次 (m = ((p << 64) + b - 1) / b)，差異是 (b - 1)/b。
原本：

___m = (~0ULL / ___b) * ___p;                                       
___m += (((~0ULL % ___b + 1) * ___p) + ___b - 1) / ___b; 

現在：

___m = (~0ULL / ___b) * ___p;                                   
___m += ((~0ULL % ___b + 1) * ___p) / ___b; 

Answer: 帶入數值分析 (b = 10, p = 16)
初始 m 值計算：

• ~0ULL / 10 = (
  $2^{64}$ -1) / 10 = 1844674407370955161
• ~0ULL % 10 = (
  $2^{64}$ -1) % 10 = 5
• (5 + 1) * 16 = 96
• (96 + 10 - 1) / 10 = 105 / 10 = 10.5 → 10
• m = 1844674407370955161 * 16 + 10 = 29514790517935282586

重新計算 m：

• (~0ULL % 10 + 1) * 16 / 10 = 96 / 10 = 9.6 → 9
• m = 1844674407370955161 * 16 + 9 = 29514790517935282585

差值確實是 1，相當於 (b-1)/b 的整數部分。

為什麼需要這種差異？
原因在於精度與偏差（bias）的平衡：

1. 無偏差時：m 需要更精確，包含 (b-1)/b 這個微小調整。
2. 有偏差時：計算變為 (n*m + m)/2^64/p，公式本身已包含補償，所以 m 不需要額外的 (b-1)/b 調整。

數學原理：

• 無偏差時，我們希望 (n*m)/
  $2^{64}$ /p ≈ n/b
• 有偏差時，我們希望 (n*m + m)/
  $2^{64}$ /p ≈ n/b

---

判斷是否需要 bias - 一般情況的處理

如果 m 足夠精確，代碼進入優化階段，試圖簡化 m 和 p 的值（移除 m 和 p 的共有因子）。

在數學上，如果 m/p = k 是一個常數，那麼

$\frac{\frac{m}{c}}{\frac{p}{c}}$ = k 仍然等於相同的常數。這裡代碼試圖找到 m 和 p 的最大公約數（或近似值），然後同時除以它，保持比值不變但降低數值大小。

uint32_t ___bits = -(___m & -___m);
___bits |= ___m >> 32;
___bits = (~___bits) << 1;

if (!___bits) {
    ___p /= (___m & -___m);
    ___m /= (___m & -___m);
} else {
    ___p >>= ilog2(___bits);
    ___m >>= ilog2(___bits);
}

1. ___bits = -(___m & -___m)：
   • (___m & -___m) 獲取 m 的最低位 1
   • 取負數後，得到一個 mask，其中所有高於這個位的位都是 1
2. ___bits |= ___m >> 32：
   • 加入 m 的高 32 位，確保考慮 m 的所有重要位元
3. ___bits = (~___bits) << 1：
   • 取反並左移，創建一個標記位元 mask
4. 根據 ___bits 的值決定如何縮小 ___m 和 ___p：
   • 如果 ___bits == 0：表示必須設置 bit 31，此時除以 m 的最低位 1
   ​​​​___p >>= ilog2(___bits);
   ​​​​___m >>= ilog2(___bits);
   

   • 其他情況：根據 ___bits 的最高位 1 來決定移位量，同時對 m 和 p 右移
   ​​​​___p /= (___m & -___m);
   ​​​​___m /= (___m & -___m);
   

---

兩個條件的組合

Now we have a combination of 2 conditions:
Condition 1 : whether or not we need to apply a bias.
Condition 2 : whether or not there might be an overflow in the cross product determined by (___m & ((1 << 63) | (1 << 31))).
Select the best way to do (m_bias + m * n) / (1 << 64).
From now on there will be actual runtime code generated.

• 條件 1：是否需要偏差
  • __xprod_64 函式中判斷是否需要偏差來抵消因位元截斷造成的誤差。
• 條件 2：交叉乘積是否會發生溢位
  • 透過呼叫 __xprod_64 函式檢查 (___m & ((1 << 63) | (1 << 31))) 來確定是否會發生溢位。
  • 1 << 63 和 1 << 31 分別是第 63 位和第 31 位。
  • 如果 ___m 的第 63 位和第 31 位被設置，可能會導致乘法運算中的溢位問題。

根據這兩個條件，選擇最佳的計算方法來進行 (m_bias + m * n) / (1 << 64) 的運算。這裡 m_bias 是偏差，m 是乘數，n 是被除數的一部分。
__xprod_64 會根據是否應用偏差來決定如何進行運算，以避免溢位並保證結果的準確性。
___p 用於將結果進一步調整到正確的範圍內。

___res = __xprod_64(___m, ___n, ___bias);
___res /= ___p;

---

#define __div64_const32(n, ___b)                                            \
    ({                                                                      \
        uint64_t ___res, ___x, ___t, ___m, ___n = (n);                      \
        bool ___bias;                                                       \
                                                                            \
        /* determine MSB of b */                                            \
        uint32_t ___p = 1 << ilog2(___b);                                   \
                                                                            \
        /* compute m = ((p << 64) + b - 1) / b */                           \
        ___m = (~0ULL / ___b) * ___p;                                       \
        ___m += (((~0ULL % ___b + 1) * ___p) + ___b - 1) / ___b;            \
                                                                            \
        /* one less than the dividend with highest result */                \
        ___x = ~0ULL / ___b * ___b - 1;                                     \
                                                                            \
        /* test our ___m with res = m * x / (p << 64) */                    \
        ___res = ((___m & 0xffffffff) * (___x & 0xffffffff)) >> 32;         \
        ___t = ___res += (___m & 0xffffffff) * (___x >> 32);                \
        ___res += (___x & 0xffffffff) * (___m >> 32);                       \
        ___t = (___res < ___t) ? (1ULL << 32) : 0;                          \
        ___res = (___res >> 32) + ___t;                                     \
        ___res += (___m >> 32) * (___x >> 32);                              \
        ___res /= ___p;                                                     \
                                                                            \
        if (~0ULL % (___b / (___b & -___b)) == 0) {                         \
            /* special case, can be simplified to ... */                    \
            ___n /= (___b & -___b);                                         \
            ___m = ~0ULL / (___b / (___b & -___b));                         \
            ___p = 1;                                                       \
            ___bias = true;                                                 \
        } else if (___res != ___x / ___b) {                                 \ 
            ___bias = true;                                                 \
            /* Compute m = (p << 64) / b */                                 \
            ___m = (~0ULL / ___b) * ___p;                                   \
            ___m += ((~0ULL % ___b + 1) * ___p) / ___b;                     \ //AAAA
        } else {                                                            \
            uint32_t ___bits = -(___m & -___m);                             \
            ___bits |= ___m >> 32;                                          \
            ___bits = (~___bits) << 1;                                      \
            if (!___bits) {                                                 \
                ___p /= (___m & -___m);                                     \
                ___m /= (___m & -___m);                                     \
            } else {                                                        \
                ___p >>= ilog2(___bits);                                    \ //BBBB
                ___m >>= ilog2(___bits);                                    \ //CCCC
            }                                                               \
            /* No bias needed. */                                           \
            ___bias = false;                                                \
        }                                                                   \
                                                                            \ 
        ___res = __xprod_64(___m, ___n, ___bias);                           \
                                                                            \
        ___res /= ___p;                                                     \
    })

---

6. __xprod_64 計算乘積
   執行的是一個 128 位精度的乘法，然後取高 64 位的結果（將 64×64 位乘法分解為四個 32×32 位乘法），考慮是否需要偏差來處理乘法運算中的溢位問題。

計算乘積的低位

retval = ((bias ? m : 0) + m * n) >> 64 根據是否應用偏差進行條件判斷：
Case 1：無偏差。
Case 2：有偏差但無溢出風險。
Case 3：有偏差且有溢出風險。

if (!bias) {
    res = ((uint64_t) m_lo * n_lo) >> 32;
} else if (!(m & ((1ULL << 63) | (1ULL << 31)))) {
    res = (m + (uint64_t) m_lo * n_lo) >> 32; // 加上 m 值（偏差修正），然後取高 32 位
} else {
    res = m + (uint64_t) m_lo * n_lo;
    res_lo = res >> 32;
    res_hi = (res_lo < m_hi);
    res = res_lo | ((uint64_t) res_hi << 32);
}

---

計算乘積的高位部分

計算 m_lo * n_hi 和 m_hi * n_lo，並將結果右移 32 位。
Case 1：無溢出風險。
Case 2：有溢出風險。
加上 m_hi * n_hi 的乘積，最終結果存儲在 res 中並返回。

if (!(m & ((1ULL << 63) | (1ULL << 31)))) { // 檢查 m 的第 63 位和第 31 位是否為 0
    res += (uint64_t) m_lo * n_hi;
    res += (uint64_t) m_hi * n_lo;
    res >>= 32;
} else {
    res += (uint64_t) m_lo * n_hi;
    uint32_t tmp = res >> 32;
    res += (uint64_t) m_hi * n_lo;
    res_lo = res >> 32;
    res_hi = (res_lo < tmp); // 檢測加法是否產生進位
    res = res_lo | ((uint64_t) res_hi << 32);
}
res += (uint64_t) m_hi * n_hi;

---

7. 二進制長除法
   進行 64 位元數字除以 32 位元整數的除法運算，返回餘數。
   透過初步減少高位元和減法迭代操作，避免使用硬體除法指令，非常適合在不支持硬體除法的情況下使用。

快速處理高位元

res = 0;
if (high >= base) {
    high /= base;
    res = (uint64_t) high << 32;
    rem -= (uint64_t) (high * base) << 32;
}

帶入數值分析：計算 12345678901234 ÷ 100：

1. 提取高位：高 32 位 = 2874452
2. 計算高位的商：2874452 ÷ 100 = 28744
3. 設置初始結果：res = 28744 << 32
4. 更新餘數：rem = 原值 - (28744 * 100) << 32

這個初步優化可能大幅減少計算時間。

標準二進制長除法算法

static inline uint32_t __div64_32(uint64_t *n, uint32_t base)
{
    uint64_t rem = *n;
    uint64_t b = base;
    uint64_t res, d = 1;
    uint32_t high = rem >> 32;

    // ...

    /* 將除數 b 左移到盡可能大而又不超過餘數 */
    while ((int64_t) b > 0 && b < rem) {
        b = b + b; // 除數乘 2
        d = d + d; // 對應位的貢獻量也乘 2
    }
    
    // ...
}

帶入數值分析：計算 1234 ÷ 10：
初始值：rem = 1234, b = 10, d = 1, res = 0
第一階段（左移）：

• 迭代 1 : b = 20, d = 2
• 迭代 2 : b = 40, d = 4
• 迭代 3 : b = 80, d = 8
• 迭代 4 : b = 160, d = 16
• 迭代 5 : b = 320, d = 32
• 迭代 6 : b = 640, d = 64
• 迭代 7 : b = 1280, d = 128（停止，因為1280 > 1234）
• 回退: b = 640, d = 64（最大不超過餘數的值）

    
    //...

    /* 從最大的除數開始，逐步右移並嘗試減法 */
    do {
        if (rem >= b) {
            rem -= b; // 如果餘數大於等於除數，則減去
            res += d; // 將對應位的貢獻量加到結果
        }
        b >>= 1; // 除數右移（除2）
        d >>= 1;  // 對應位的貢獻量也右移
    } while (d);

    *n = res;
    return rem;
}

第二階段（右移與減法）：

• 檢查: 1234 >= 640? 是，rem = 594, res = 64
  b = 320, d = 32
• 檢查: 594 >= 320? 是，rem = 274, res = 96
  b = 160, d = 16
• 檢查: 274 >= 160? 是，rem = 114, res = 112
  b = 80, d = 8
• 檢查: 114 >= 80? 是，rem = 34, res = 120
  b = 40, d = 4
• 檢查: 34 < 40? 是（不減），res 不變
  b = 20, d = 2
• 檢查: 34 >= 20? 是，rem = 14, res = 122
  b = 10, d = 1
• 檢查: 14 >= 10? 是，rem = 4, res = 123
  b = 5, d = 0（停止，因為d = 0）

結果：商 = 123，餘數 = 4
這正是 1234 ÷ 10 = 123 餘 4。

---

8. 64 位元除法
   進行 64 位元數字除法，根據除數是否為 2 的冪、是否為常數等情況選擇不同的計算方法，分成以下四個 case，並返回餘數 __rem。
   把除數是否為常數拿出來討論的原因是，如果除數不是常數，編譯器在編譯時無法進行優化，因此必須在運行時執行更通用的除法操作。
   case 1：如果 __base 是常數且是 2 的冪。
   case 2：如果 __base 是常數且不為 0。
   case 3：如果 n 小於 2³²。
   case 4：__base 不是編譯時常數，並且 n 是一個大於 2³² 的大數。

補充：__builtin_constant_p(__base) 會是 false 的情況

1. 變數是一個在運行時返回值的函數

uint32_t base = get_base_value(); // Assume get_base_value() is a function that returns a value at runtime
uint64_t n = 5;
uint32_t rem = div64(n, base);

2. 變數是透過參數傳遞

void my_function(uint32_t base) {
    uint64_t n = 5;
    uint32_t rem = div64(n, base);
}

3. 變數是透過變數賦予

uint32_t base;
scanf("%u", &base); // Get value from standard input
uint64_t n = 5;
uint32_t rem = div64(n, base); 

#define div64(n, base)                                               \
    ({                                                               \
        uint32_t __base = (base);                                    \
        uint32_t __rem;                                              \
        if (__builtin_constant_p(__base) && is_power_of_2(__base)) { \
            __rem = (n) & (__base - 1);                              \
            (n) >>= ilog2(__base);                                   \ //DDDD
        } else if (__builtin_constant_p(__base) && __base != 0) {    \
            uint32_t __res_lo, __n_lo = (n);                         \
            (n) = __div64_const32(n, __base);                        \
            /* the remainder can be computed with 32-bit regs */     \
            __res_lo = (n);                                          \
            __rem = __n_lo - __res_lo * __base;                      \
        } else if (likely(((n) >> 32) == 0)) {                       \
            __rem = (uint32_t) (n) % __base;                         \
            (n) = (uint32_t) (n) / __base;                           \
        } else {                                                     \
            __rem = __div64_32(&(n), __base);                        \
        }                                                            \
        __rem;                                                       \
    })

9. 使用 div64 對 res 進行除法運算，並輸出商數。
   如果想得到餘數，則使用一變數接住巨集回傳的值，即 uint64_t r = div64(res, 3); ，則 r 為餘數。

argc 的值
argc 是命令行參數的數量。當執行程式時，命令行的每一個參數都會被計算在內，包括程式的名稱。
例如，如果在命令行中輸入 ./program arg1 arg2，則 argc 的值為 3（包括程式名和兩個參數）。因此，argc 的值取決於如何執行程式，它不是每次都固定的。
如果以 ./program 執行程式（即沒有額外的參數），argc 是 1。

int main(int argc, char *argv[])
{
    uint64_t res;
    res = argc + 5;
    div64(res, 3);
    printf("%lu\n", res);
}

在 Linux 核心原始程式碼找到相關的程式碼並設計實驗 (使用 Linux 核心模組) 探討程式碼的效益和限制

Linux 核心原始程式碼：include/asm-generic/div64.h

實驗目的

探討 do_div() 巨集在不同除數情境下的表現。

實驗步驟

步驟1：撰寫測試模組
撰寫一個 Linux 核心模組，實現以下功能：

• 定義多個除法測試情境，包括常數除數和非常數除數。
• 已知 BITS_PER_LONG == 32，所以可以分成以下四個情境討論。
• 記錄每個測試情境下 do_div() 的執行時間和結果。

#include <linux/init.h>
#include <linux/module.h>
#include <linux/kernel.h>
#include <linux/timekeeping.h>
#include <linux/types.h>
#include <linux/delay.h>
#include <asm-generic/div64.h>

MODULE_LICENSE("GPL");
MODULE_AUTHOR("Yi-Han Wan");
MODULE_DESCRIPTION("A test module for do_div macro.");
MODULE_VERSION("0.1");

static uint64_t test_division(uint64_t n, uint32_t base) {
    uint64_t start, end;
    uint32_t remainder;

    start = ktime_get_ns();
    remainder = do_div(n, base);
    end = ktime_get_ns();

    printk(KERN_INFO "n: %llu, base: %u, remainder: %u, quotient: %llu, time: %llu ns\n", n, base, remainder, n, end - start);

    return end - start;
}

static int __init div_test_init(void) {
    printk(KERN_INFO "Loading do_div test module...\n");

    // Test case 1: __base is a constant and power of 2
    test_division(123456789012345ULL, 8); // 8 is a power of 2

    // Test case 2: __base is a constant and not power of 2
    test_division(123456789012345ULL, 10); // 10 is not a power of 2

    // Test case 3: n's high 32 bits are 0 and __base is not a constant
    uint32_t dynamic_base = jiffies % 100 + 1; // dynamic_base is not a compile-time constant
    test_division(12345ULL, dynamic_base); // n fits in 32 bits

    // Test case 4: General case where n's high 32 bits are non-zero and __base is not a constant
    dynamic_base = jiffies % 10000 + 1; // dynamic_base is not a compile-time constant
    test_division(98765432109876ULL, dynamic_base); // General case with high bits non-zero

    return 0;
}

static void __exit div_test_exit(void) {
    printk(KERN_INFO "Unloading do_div test module...\n");
}

module_init(div_test_init);
module_exit(div_test_exit);

Jiffies 是什麼？
Jiffies 為 Linux 核心變數( 32 位元變數，unsigned long)，它被用來紀錄系統自開機以來，已經過多少的 tick。每發生一次 timer interrupt，Jiffies 變數會被加一。值得注意的是，Jiffies 於系統開機時，並非初始化成零，而是被設為 -300 HZ (arch/i386/kernel/time.c)，即代表系統於開機五分鐘後，jiffies 便會溢位。那溢位怎麼辦？事實上，Linux 核心定義幾個巨集(timer_after, time_after_eq, time_before 與 time_before_eq)，即便是溢位，也能藉由這幾個巨集正確地取得 jiffies 的內容。

步驟2：編譯並載入模組
編寫 Makefile 以編譯模組：

obj-m += div_test.o

all:
    make -C /lib/modules/$(shell uname -r)/build M=$(PWD) modules

clean:
    make -C /lib/modules/$(shell uname -r)/build M=$(PWD) clean

編譯模組：

make

載入模組：

sudo insmod div_test.ko

查看核心訊息以獲取測試結果：

sudo dmesg | grep "do_div test"

效能評估

• 比較不同情境下 do_div() 的執行時間

[98608.114361] Loading do_div test module...
[98608.114364] n: 15432098626543, base: 8, remainder: 1, quotient: 15432098626543, time: 18 ns
[98608.114366] n: 10000549940, base: 12345, remainder: 3045, quotient: 10000549940, time: 18 ns
[98608.114367] n: 1234, base: 10, remainder: 5, quotient: 1234, time: 21 ns
[98608.114369] n: 29838499126, base: 3310, remainder: 2816, quotient: 29838499126, time: 33 ns

計算常數除數的平均時間與非常數除數的平均時間，然後計算效率提升（Performance Gain）：

$EfficiencyGain={\displaystyle \frac{AvgTime(Non-Constant)-AvgTime(Constant)}{AvgTime(Non-Constant)}}\ast 100$

計算結果顯示，使用常數除數時的效率提升約為 33.3%。

• 分析常數除數和非常數除數對效能的影響

當除數是常數：對於 2 的次方的除數，編譯器可以使用右移操作來代替除法操作（例如 n >> 3 來代替 n / 8），這比實際的除法運算要快得多。case 2 中，雖然除數不是 2 的次方，但由於除數是常數，編譯器可以進行一些數學優化，使得除法運算更高效。

當除數不是常數：
無法使用位元運算或乘法反元素等編譯時期優化技術，必須在運行時進行實際的除法運算，導致運行時間較長。當 n 的值較大且 base 不是常數時，除法運算需要進行完整的 64 位整數除法，這是一個相對複雜且耗時的運算過程，因此運行時間最長。

實驗結果分析

• 在不同 case 下的效益和限制

  case 1：常數除數且為 2 的次方
  效能：右移運算比標準除法運算快很多，因為位運算在硬體上直接支持且執行速度快。
  限制：僅適用於除數為 2 的次方的情況，其他除數不適用。

  case 2：常數除數但不是 2 的次方
  效能：編譯期可以進行一些優化，例如使用乘法和位移來模擬除法，從而提升運行速度。
  限制：僅在除數為常數的情況下有效，變量除數無法應用這種優化。

  case 3：除數是動態變數，且 n 的高 32 位為 0
  效能：直接使用 32 位的除法運算，效率相對較高。
  限制：僅適用於 n 可以用 32 位表示的情況，對於更大的數值無法應用。

  case 4：除數和被除數都是變量且超過 32 位
  效能：需要進行 64 位的除法運算，並且在運行時動態計算，速度較慢。
  限制：無法進行編譯時期優化，完全依賴於運行時的計算，效能較低。