手刻一套 C 語言大數運算之物件導向函式庫

contributed by <Shiritai>

開發環境

開發環境

$ gcc --version
gcc (Ubuntu 11.3.0-1ubuntu1~22.04) 11.3.0
$ lscpu
Architecture:            x86_64
  CPU op-mode(s):        32-bit, 64-bit
  Address sizes:         39 bits physical, 48 bits virtual
  Byte Order:            Little Endian
CPU(s):                  12
  On-line CPU(s) list:   0-11
Vendor ID:               GenuineIntel
  Model name:            Intel(R) Core(TM) i7-8750H CPU @ 2.20GHz
    CPU family:          6
    Model:               158
    Thread(s) per core:  2
    Core(s) per socket:  6
    Socket(s):           1
    Stepping:            10
    CPU max MHz:         4100.0000
    CPU min MHz:         800.0000
    BogoMIPS:            4399.99

Arbitrary bit stream，包含經典 bitwise 操作與算數運算。
小至單個位元，大至加減乘法，目標是成為 bit stream 的主人。
另外，本 Repository 目標是可使用於 Linux Kernel 或 user space 中。
透過條件編譯適應不同環境。

實現大數運算有很多方法，比如暴力用一個整數表一位元、字串法，作業文件的實作方法也很有趣。不過為了不被該實作方法影響，我先實作了一個版本，經比較後發現我的作法可能更為激進，且採用諸多課程提及的技巧，為一次非常全面的練習經驗，請容我娓娓道來。

Repository: Shiritai/inf_bits

Check list

對菜鳥如我而言，實作的一路上可說是困難重重，至今尚未完成，以下列表紀錄我走過的路子。

• Repository
  • 建立 repository
  • 引入輔助工具 (clang-format, git commit hook)
  • 整理程式碼並分批 commit (源自於我寫程式的不良習慣，本來應該完成部分功能就 commit 的)
  • 更新遠端至目前版本
• 設計物件導向風格之 inf_bits
  • 成員、方法定義，公開方法包成結構體
  • 切成 big_num.h、big_num.c，前者表給使用者的 API，後者為實作
  • 撰寫 makefile 完成自動化分別編譯
    • 自動化帶參數之 gdb、運行 python 腳本
    • 自動化 valgrind 記憶體檢測
    • 自動化 kernel 內的編譯
  • 利用多型設計有、無號版的 inf_bits
    • 搭配如 Linux 紅黑樹的著色設計決定 signed、unsigned
• 實作與測試
  • 消滅 magic number
  • 物件生命週期相關方法
    • copy_test: construct, finalize, copy and related methods, equals, is_zero, show_hex
  • 輸出方法
  • Zeros 相關方法
  • 單個位元運算
    • single_bit_test: mark, clear, flip, test
  • 位移運算
    • shift_test: all shift methods
  • 加減法與相關運算
    • inc_dec_test with Python script: dump, increment, decrement, shift (in place)
    • sub_test: add, sub
    • dump_fib with Python script: add related, assign
  • 乘法 (還有 bug…)
    • arith_test with Python script: add_in_place, sub, increment, decrement, mul, assign
    • mul_test with Python script: add, shift related, mul
  • 準備測試用工具
  • 調整測試工具為: 針對公開、私有方法驗證
• 永無止境的重構
• 解決在 kernel 編譯的問題
  • 針對 5.19.0 版除了 dump 之外可行
  • 其他版本的解決方法
• 完成 HackMD 筆記
• 效能議題
  • 效能測量與圖表繪製
  • 參考教授、同儕的實作並提出改進方案

設計

大數運算的容器 inf_bits 我以 C 風格之物件導向實現。以下將分為成員與方法來介紹，設計內容都位於 big_num.h，該檔案定義所有欲提供用戶使用的內容。

物件導向是一種態度。jserv
詳細可見你所不知道的 C 語言：物件導向程式設計篇。

物件與成員

由於 64 bits CPU 硬體可以直接支援的資料長度最大為 64 bits，故起出設計時就採用 uint64_t 為資料儲存型別的首選，這與改善方案 4: 善用 64 位元微處理器特性不謀而合。

TODO: 考慮 32/64 位元微處理器架構的特性

jserv

不過相較之下，我希望自己的 big_num 直接模擬一串 bit string，可以以二補數的方式解釋並運算，而無需在意正負號，這點較為激進，但可以因此獲得二補數提供的一些優良性質。

對不同數字需要的位元數不同，故我以 vector 動態陣列資料結構的設計發想，設計 capacity 和 use 兩成員，分別為位元容量和使用位元數。特別注意 use 實際上與 bignum 不同，use 更進一步地表示「從最高位 1 以降，使用多少個位元」。另外為了提升多數運算的效能，我設計成員 low，表最低位之 1 bit。

另外很有趣的是，對於 bit stream 的部分，我使用上課提及之零長度陣列的技巧設計物件。最終結構體定義如下。

struct big_num {
    int capacity;      // bit capacity
    int use;           // bit length using
    int low;           // starting bit which is not zero
    bn_data_t arr[0];  // determines in runtime
};

/**
 * A infinite bit stream container supporting many bitwise and arithmetic
 * operations
 */
typedef struct big_num **inf_bits;

注意用語: a pointer to a pointer 是間接指標，關鍵在於記憶體操作的相依關係，而非「二重」，請查閱辭典。

jserv

特別注意到末行定義 inf_bits 為 struct big_num ** 作為出現在使用者面前的型別。之所以需要間接指標是因為

1. 物件中有零長度陣列，這意味著每次產生新物件時基本上都會呼叫 malloc 等函式分配記憶體，這使本物件至少要是指標型別。雖然也可以在 stack 記憶體上初始化，但這僅限使用者在自己的 scope 上使用，違反物件導向從物件的生成到死亡都由實作方完成的習慣。
2. 設計之 inf_bits 應支援 in place 操作，比如 in place 的左右移、單個位元的設定與加減法等。如果物件的實體只由普通指標管理，由於產生物件後 arr 成員的長度也確定下來，今若需改變長度，只有重新建立物件才行，此時舊的指標會被拋棄，那就需要想辦法讓使用者可以取得新的指標。對此，我們可以提高一層視角，使用指標的指標。第一層指向結構指標，第二層才是真的指向結構。這樣一來就可以讓使用者在不替換變數的情況下修改內部結構的長度。

方法封裝結構

對外公開的方法使用結構體 big_num_op 封裝，並留下一全域實例 bn_methods 以供使用者呼叫方法。

// in file big_num.h
/**
 * inf_bits is a infinite bit stream supporting
 * many bitwise and arithmetic operations such as
 * addition, subtraction and left/right shift
 *
 * Support multiplication in simulative way
 */
struct big_num_ops {
    /* ... */
}

/**
 * Methods of inf_bits object
 * 
 * Initialize bn_methods by calling big_num_init()
 */
extern struct big_num_ops bn_methods;

與你所不知道的 C 語言：物件導向程式設計篇中 stack 的設計不同的是，我將成員與方法分成兩個結構體，是因為將 function pointer 嵌入實例會消耗額外的空間，故獨立出來使用。以下提供用例:

inf_bits bn = bn_methods.construct(0xffffffff); // new inf_bits
bn_methods.show_hex(bn);                        // print
bn_methods.left_shift_in_place(bn, 3);          // left shift
bn_methods.show_hex(bn);                        // print
bn_methods.finalize(bn);                        // delete

看起來滿賞心悅目的?!
Shiritai

避免 camel case，以便融入 Linux 核心的程式碼風格

jserv

我也在猶豫對外的名稱，想說 struct big_num 本身就是 snail case，對外 struct big_num ** 的名稱才改成 camel case。
不過這名字有點沒創意，等想到比較好的就來換 Shiritai
sysprog21/bignum 就用很單調但好理解的 bn，命名風格一致即可。

jserv

Ok, 那就取叫 inf_bits 好了，最符合設計的初衷。對內為 struct big_num，對外為 inf_bitsShiritai

封裝

物件導向一大特色即封裝，規定成員與方法的私有性來避免用戶不正當的使用物件。我不知道 C 語言要使成員私有化應如何完成，不過方法的私有化倒不難。利用 static 關鍵字便可編譯出僅檔案內部知曉函式位置的 object file。使用者透過呼叫我為 big_num.h 留下的函式 big_num_init (實作於 big_num.c) 初始化 bn_methods 的成員函式指標們，將 big_num.c 中想提供給使用者調用的函式位置記在 bn_methods 內，使用者便能正常的呼叫方法。

注意到以下 new 和 delete 方法因為與 C++ 運算子撞名，導致 clang-format 會錯誤的為其前方加上一空白，故改命名為 bn_XXX。

依 jserv 建議，改命名為 construct, finalize，清楚易懂。
Eroiko

// in big_num.h
void big_num_init();

// in big_num.c
struct big_num_ops bn_methods;
void big_num_init()
{
    bn_methods = (struct big_num_ops){
        .construct = construct,
        .finalize = finalize,
        .show_hex = show_hex,
        /* ... */
    }
}

實作

接下來才是重頭戲：盡可能高效的使用奇技淫巧實現眾多方法，並測試其正確性。

目前除錯進行到 inf_bits::mul，對乘數小於 bn_data_t 長度時通過測試，但大於後就會抽搐…細節詳見測試進度。

我認為這應該是自己目前為止學習並使用 bitwise 操作的集大成了。
Shiritai
善用 _Static_assert

jserv

目前我採用巨集展開搭配 gcc predefined macros 增加除錯訊息的可讀性，稍後於測試框架介紹。

然後想請問 _Static_assert 的使用是是否是在實作中而非測試的時候？
Shiritai
參見 第 6 週測驗題 的測驗 4，_Static_assert 可用於實作和測試階段。

jserv

實作前的準備

後續將會使用大量 bitwise 操作，自然會用到很多神奇的位移量、各種常數和 mask。可以預期如果這些 Magic number 充斥整份程式碼，未來調整時將多麽痛苦。為了防止該慘劇，我先定義很可能常用的常數，預期是消滅所有魔法。

#define __BITS_IN_BYTE 8
#define __BITS_IN_ULL (sizeof(bn_data_t) * __BITS_IN_BYTE)

/**
 * Offset of ULL from 0: i.e. (bits of ULL) - 1
 */
#define __ULL_OFFSET (__BITS_IN_ULL - 1)

/**
 * Default 4 words, 2 long long
 */
#define __DEFAULT_ULL_LEN 2
#define __ULL_ONES ((bn_data_t) ~0)
#define __LOG_BITS_IN_ULL (LOG2(__BITS_IN_ULL))

可以注意到當中使用 LOG2 巨集，其引自數值系統的介紹。另外受 Linux 中 list_head API 的啟發，我定義數個巨集函式，讓諸多位元運算的意義一目瞭然，以下舉最重要的兩例。

/**
 * Get index limit (upper bound) of given big_num instance
 * Round up
 */
#define bn_up_index(bits) \
    (((bits) >> __LOG_BITS_IN_ULL) + !!((bits) & __ULL_OFFSET))
/**
 * Get index base (lower bound) of given big_num instance
 * Round down
 */
#define bn_dn_index(bits) ((bits) >> __LOG_BITS_IN_ULL)

bn_up_index 的意義在給定位元數

$bits$，得到可包含該所有位元的 bn_data_t 數。比如有

$65$ 個位元，那便需要

$2$ 個 unsigned long long 才能裝得下。

該巨集的實作為 (((bits) >> __LOG_BITS_IN_ULL) + !!((bits) & __ULL_OFFSET))，分為前後兩段。前者相當於對 bn_data_t 的長度進行整數除法，後者為確認餘數是否非零，非零則透過 !! 轉為 1。

物件生命週期相關方法

下圖所有 in-degree 為零的節點都是公開方法。

很複雜吼，這也許是 pointer-to-pointer 的缺點吧，或者說我不熟悉 pointer-to-pointer 更良好的使用方式。
Shiritai

配置與釋放

__bn_malloc 和 __bn_free 的來歷請見依賴函式庫議題，可簡單視為 malloc 和 free，可以配置、釋放 pointer 或 pointer-to-pointer。__bn_new 和 __finalize 這兩個很 C++ 味的私有方法負責完成 pointer (struct big_num 本體) 的配置、釋放，並調整 capacity 成員。construct 和 finalize 則是簡化參數，並將 pointer 和 pointer-to-pointer 的配置與釋放包在一起，開放給用戶簡潔的方法。以下舉 __bn_new 這個配置 struct big_num 本體的私有方法為例，可見其利用巨集函式增加使用 bitwise 操作的可讀性。

// implementation of private method "__bn_new"
const int sz = sizeof(struct big_num) + bn_ull_as_bytes(num_ull);
struct big_num *ret = __bn_malloc(sz);
memset(ret, 0, sz);
ret->capacity = bn_ull_as_bits(num_ull);
return ret;

複製與所有權共享 (i.e. deep copy and shallow copy)

記憶體管理中，複製物件是一個重要的議題。由於希望對內設計的邏輯盡可能不要重複且功能充足，對外簡單明瞭而無需知道細節。首先是 inner_copy，負責處物件本體的複製，bn_copy 將其包裝間接指標，copy 則隱藏 bn_copy 不必要公開的參數，讓使用者傳入舊的 inf_bits 即可。

// implementation of public method "copy"
return bn_copy(self, bn_up_index((*self)->use));

assign 對外的意義為複製參考，也就是讓新的 inf_bits 共享舊的 inf_bits 的所有權 (並不會轉移)，以費氏數列的計算為例，使用方法如下。

// in test.c::dump_fib, i.e. user perspective
// doing (a, b) = (a + b, a)
inf_bits a_ = bn_methods.copy(a); // a_ = a
bn_methods.add_in_place(a, b);    // a += b
bn_methods.finalize(b);          // deconstruct b
// void assign(struct big_num ***dest, struct big_num **src)
// notice that inf_bits is struct big_num **
bn_methods.assign(&b, a_);        // b = a_

實作上，由於要確保

1. 修改 dest 的值
2. 不修改也不釋放 dest 內部的物件，因為不確定是否有其他 pointer-to-pointer 持有之
3. 不修改 src 及其內部物件的值

可見 dest 需為 pointer-to-pointer 的地址，也就是 ***；src 則保持 pointer-to-pointer 即可。而實作僅一行

// implementation of public method "assign"
*dest = src;

恭喜人生第一次用上 *** 這東西 :)
Shiritai

另外擴容的私有方法 expand 相對簡單，就不介紹了。

輸出方法

提供的方法有二：show_hex 和 dump，皆為公開方法，前者輸出所有成員 (Metadata) 和 16 進制的 bit stream 內容，每

$64$ 位元以 _ 隔開；後者將 16 進制的 bit stream 內容寫入指定檔案。

以下以 show_hex 為例，當中 format string 使用將 i 非零時轉成 1 的技巧來存取字元陣列；同時由於依賴函式庫的議題，印出使用的 __bn_print 將由適當的函式取代。

// implementation of public method "show_hex"
__bn_print(
    "[Metadata] addr: %p cap: %d, use %d, low: %d\n[Data] addr: %p, arr: ",
    &(*self)->capacity, (*self)->capacity, (*self)->use, (*self)->low,
    (*self)->arr);
int cap_ull = bn_up_index((*self)->capacity);
for (int i = cap_ull - 1; i >= 0; --i)
    __bn_print("%016lx%c", (*self)->arr[i], "\n_"[!!i]);

Zeros 相關方法與特殊情況

Zero 相關方法對效能乃至於正確性而言至關重要。一旦平時維護 use、low 兩成員，即可在各式計算過程中略過大量未使用的位元，增加幾乎所有運算的效率。維護所需的成本在不同方法中有所不同。Worse case 複雜度與位元數成正比，__bn_c{l,t}z 私有方法正是在最糟情況下維護該成員們的手段。

上圖中 in-degree 為零者為公開方法。有了 use 和 low 成員，is_zero 的實作便僅需一行。

// implementation of public method "is_zero"
return (*self)->use == (*self)->low;

Counting Zero

實作的根源是 Counting Zero 的實作中我使用的生成技巧，搭配以下兩行完成函式的生成:

clz(bn_data_t, bn_data);  // clz_bn_data
ctz(bn_data_t, bn_data);  // ctz_bn_data

便能使用。亦可使用 __builtin_c{l,t}zll，不過要注意輸入為零時存在的 undefined behavior，我同時採用這兩者，根據條件編譯決定要如何使用，預計用於效能測試。以下舉 __bn_clz 為例，可見 __builtin_c{l,t}z 的版本稍微比較不美觀一點。

以下為版本一，branching 判斷 (*self)->arr[i] 是否為零:

// implementation of private method "__bn_clz"
    int cap_ull = bn_up_index((*self)->capacity);
    int res = 0, i = cap_ull;
    do {
        --i;
#ifdef MY_CZ
        res += clz_bn_data((*self)->arr[i]);
#else
        if ((*self)->arr[i])
            res += __builtin_clzll((*self)->arr[i]);
        else // avoid undefined behavior
            res += __BITS_IN_ULL;
#endif
    } while (i > 0 && !(*self)->arr[i]);
    // automatically refresh use member
    (*self)->use = (*self)->capacity - res;
    return res;

為了避免 branching，我實現以下版本，利用將

$0,1$ 映射為全一、全零的技巧實現。

  #else
+         bn_data_t _is_zero = !(*self)->arr[i];
+         res += -!_is_zero & __builtin_clzll((*self)->arr[i]) + -_is_zero & __BITS_IN_ULL;
-         if ((*self)->arr[i])
-             res += __builtin_clzll((*self)->arr[i]);
-         else // avoid undefined behavior
-             res += __BITS_IN_ULL;
  #endif

經過 __bn_c{l,t}z 方法調整 use 和 low 成員後，bn_c{l,t}z 便能直接取用結果回傳。

// implementation of public method "bn_clz"
/**
 * Happily we've maintained capacity and use,
 * just use them directly :)
 */
return (*self)->capacity - (*self)->use;

另由於我將 branching 改為 branchless，自然會想做實驗確認前後的差異，我以本文最前方的開發環境，以 cc -O2 -std=c99 -S big_num.c 編譯得到 x86 組語，擷取重點如下:

• Version 1 (branch in loop)
  ​​​​.L70:
  ​​​​    movq	8(%rdi,%rdx,8), %rax
  ​​​​    testq	%rax, %rax
  ​​​​    je	.L68
  ​​​​    bsrq	%rax, %rax
  ​​​​    xorq	$63, %rax
  ​​​​    addl	%eax, %ecx
  ​​​​.L69:
  ​​​​    subl	%ecx, %esi
  ​​​​    movl	%esi, 4(%rdi)
  ​​​​    ret
  ​​​​.L68:
  ​​​​    subq	$1, %rdx
  ​​​​    addl	$64, %ecx
  ​​​​    testl	%edx, %edx
  ​​​​    jg	.L70
  ​​​​    jmp	.L69
  

  迴圈的範圍很詭異，竟然是 .L68 和 .L70 之間來回跳，.L69 則是收尾回傳。
• Version 2 (branchless in loop)
  ​​​​.L69:
  ​​​​    movq	16(%r9,%rcx,8), %rax
  ​​​​    testq	%rax, %rax
  ​​​​    sete	%sil
  ​​​​    bsrq	%rax, %rdx
  ​​​​    xorq	$63, %rdx
  ​​​​    movzbl	%sil, %r8d
  ​​​​    negq	%rax
  ​​​​    sbbl	%eax, %eax
  ​​​​    subl	%r8d, %edx
  ​​​​    andl	%edx, %eax
  ​​​​    andl	$64, %eax
  ​​​​    addl	%eax, %edi
  ​​​​    testl	%ecx, %ecx
  ​​​​    setg	%al
  ​​​​    subq	$1, %rcx
  ​​​​    testb	%sil, %al
  ​​​​    jne	.L69
  ​​​​    subl	%edi, %r10d
  ​​​​    movl	%r10d, 4(%r9)
  ​​​​    ret
  

  迴圈的範圍很容易界定，於第 18 行跳回第 1 行就是，其餘地方沒有出現轉跳。另外可見 __builtin_clzll 對應的 x86 bsrq 指令於第 5 行。

上方實驗可見 branchless 版與 branching 版行數相近，不過有趣的是 branching 版的實際上說明我的程式碼依然有改善的空間: 編譯器意識到只有當 (*self)->arr[i] 為零時，迴圈才會繼續，可以發現 branching 版只會呼叫 bsrq 指令一次！

這讓我大開眼界：表示我還能再優化！因為隨機情況下 (*self)->arr[i] 為零的機會比較小，故編譯器為我以 likely non zero 的方式安排組語，實際上這不一定與現實相同。考慮當 inf_bits 進行一定程度的加減乘法後，leading zero 可能會比想像的多，此時若以 unlikely 的邏輯安排有機會增加效率；於此同時，安排為 unlikely 的話程式只會轉跳一次後便結束函式，可以預期最多僅有一次回圈內轉跳。

於是 Version 3 登場，原本的迴圈變為以下:

    while (--i >= 0 && !(*self)->arr[i])
        ++res;
#ifdef MY_CZ
    res = (res << __LOG_BITS_IN_ULL) + clz_bn_data((*self)->arr[i]);
#else
    res = (res << __LOG_BITS_IN_ULL) + __builtin_clzll((*self)->arr[i]);
#endif

對應的迴圈部分組合語言如下:

.L38:
	addl	$1, %eax
.L36:
	subl	$1, %edx
	js	.L37
	movslq	%edx, %rcx
	cmpq	$0, 16(%rsi,%rcx,8)
	je	.L38
.L37:
	sall	$6, %eax
	movslq	%edx, %rdx
	bsrq	16(%rsi,%rdx,8), %rdx
	xorq	$63, %rdx
	addl	%edx, %eax
	subl	%eax, %edi
	movl	%edi, 4(%rsi)
	ret

與原本 Version 1 相比，Version 3 只在 .L38 間來回跳，不像原本在兩個 tag 間來回跳。且回圈內的指令變的更加簡單、轉跳距離也更近，應該是更為理想的程式碼 ;)

Refresh 方法與成員的維護

另外針對 refresh 需要特別說明，該方法有兩種不同的簽名註解。

// in big_num.c (implementation)
/**
 * Refresh "use", "low" members
 * If you're not maintain "use" and "low",
 * please call this method to make sure big_num works correctly
 */
// p.s. the implementation is:
__bn_clz(self);
__bn_ctz(self);

// in big_num.h (user perspective)
/**
 * Refresh members
 * Call this method to make sure big_num works correctly
 * after using single bitwise operations
 */

對內 refresh 方法用途多多，當進行某運算後難以直接維護 use 或 low，便可呼叫來以

$O(\mathrm{\#}ofbn\mathrm{\_}data\mathrm{\_}t)$ 複雜度重新維護之。對外原本理當盡可能不暴露此方法，但礙於單個位元運算中所有操作複雜度都是

$O(1)$，但若要維護該兩成員的話複雜度只能變為

$O(\mathrm{\#}ofbn\mathrm{\_}data\mathrm{\_}t)$，實在划不來。所以權衡一下，將此方法暴露於外，並在所有單個位元運算的簽名註解中加上提示，以 mark 方法為例。

/**
 * Note: Always call "refresh" method to make sure big_num is ready for
 * non-single-bitwise operations after using single bit operations
 *
 * Set a bit at "which" position
 * Does nothing if "which" is out of range
 */
void (*mark)(inf_bits, int);

單個位元運算

mark, clear, flip 和 test 四個方法大同小異，複雜度皆為

$O(1)$。參考你所不知道的 C 語言：bitwise 操作和數值系統閱讀紀錄，後者筆記正好用到前者針對 bitstream 的擴展：

// void Bitmap::Mark(int which)
map[which / BitsInWord] |= 1 << (which % BitsInWord);

為了避免除法跟取餘數，我改善此實作為

// implementation of public method "mark"
(*self)->arr[bn_dn_index(which)] |= 1ULL << (which & __ULL_OFFSET);

當中利用 bitwise 操作進行除法和取餘數，並利用巨集改善可讀性與消滅 Magic Number。

位移運算

i.e. >>, >>=, <<, <<= 與其變體。

圖看似複雜，實則為

• 有無 _in_place (無 _in_place 的方法直接依賴 _in_place 者)
  • 實作為複製一份後執行 _in_place 版
• 有 _long 表位移單位為 bn_data_t，沒有表位移單位為
  $bits$
  • 位移量在
    $[0,64)$ 之間
• _arb 表位移量為任意非負整數

想必讀者已經猜到 _arb 的實作就是有無 _long 的套皮。目前以上全部皆為公開方法，未來考慮只公開有 _arb 的版本，功能最強大且隱藏不必要的細節 (bn_data_t 的長度)。

實作上大量利用 bitwise 操作，舉例如下。

// inside implementation of "left_shift_in_place"
//     after init local variables and handle "expand" case
// special case
bn_data_t carry = n_bn->arr[use_ull - 1] >> (__BITS_IN_ULL - shift);
n_bn->arr[use_ull - 1] <<= shift;
if (use_ull < bn_dn_index(n_bn->capacity))
    n_bn->arr[use_ull] |= carry;  // append carrying bits

for (int i = use_ull - 2; i >= low_ull; --i) {
    bn_data_t carry = n_bn->arr[i] >> (__BITS_IN_ULL - shift);
    n_bn->arr[i] <<= shift;
    n_bn->arr[i + 1] |= carry;  // append carrying bits
}

left_shift_in_place 中 in place 的左移為了避免覆蓋未位移的資料，必須從高位開始，不過需考慮最高索引的特殊情況，故獨立 4~7 行於迴圈之外。

迴圈內先取得要移至高位的 carry，注意到在 shift 範圍為

$[0,64)$ 的條件下，左移 shift 位元至更高位的取值即往右移 __BITS_IN_ULL - shift 位的結果。

對於 use 和 low 成員的維護，left_shift 系列非常容易做到，只要增加對應的數值即可；right_shift 則不然，要考慮位移可能超過使用量而導致最後的值為負數。

我使用以下技巧完成無分支的

$max(num,0)$ 的邏輯。

// inside implementation of "right_shift_in_place"
// bitwise max(use, 0)
n_bn->use -= shift;
n_bn->use &= -!(n_bn->use < 0);

• 當 use 為正數或零，!(use < 0) 為 1，取負數為 0xFFFFFFFF (-1)，拿去 &= 後不改變值。
• 當 use 為負數，!(use < 0) 為 0，取負數為 0 (-0)，拿去 &= 後將數值設為零。

有 _long 的實作相較沒有的簡單很多，不過有一個 bug 因疏忽而抽搐了很久 0.0

// implementation of "right_long_shift_in_place"
 for (int i = low_ull >= n ? low_ull : n; i < use_ull; ++i) {
     n_bn->arr[i - n] = n_bn->arr[i]; // move
+   n_bn->arr[i] = 0ULL;              // clear
 }

加減法與相關運算

i.e. ++, --, +, +=, -，注意沒有 -=，因為我的實作方式無法做到 in place 的減法。

上圖中，bn_minus 是為了隱藏 minus 不必要暴露的參數，功能上相當於 -num；實作減法利用到加法跟取加法反元素 (bn_minus) 的特色。另外由於 minus_in_place 的簽名沒有暴露內部的實作，且有公開的意義，故直接作為公開方法，其餘都是 in-degree 為零的節點才是公開方法。

increment

++ 的實作是 += 的縮影，經依樣畫葫蘆後也能實作出 --。分為三步:

1. 加法
   由於僅需
   $+1$，利用 carry 技巧性地把進位變成迴圈的進行與否與 carry 值掛鉤，並利用邏輯判斷來更新 carry 的值。
   ​​​​const int lim = bn_up_index((*self)->capacity);
   ​​​​_Bool carry = 1;
   ​​​​for (int i = 0; carry && i < lim; ++i) {
   ​​​​    carry = (*self)->arr[i] == __ULL_ONES;
   ​​​​    ++(*self)->arr[i];
   ​​​​}
   

2. 額外進位: 擴展位元的例外
   上述加法沒有考慮位元溢位與否，由於我們的目標是任意位元的 bit stream，溢位時必須擴展。注意需擴展與否一事難以事先得知 (或者說事先得知的效率很低)，而且僅有唯一情況：所有位元為 1，故只需針對此例外單獨處理即可。
   ​​​​if (carry) {
   ​​​​    // overflow, but we're arbitrary bit stream!
   ​​​​    expand(self, lim + 1);
   ​​​​    (*self)->arr[lim] = 1;
   ​​​​}
   

3. 調整 use、low
   注意事先得知 use 和 low 要如何改變，基本上只能等加完才知道 (或者低效的掃描一遍所有位元)，不好直接維護，故直接呼叫 refresh 更實在。
   ​​​​refresh(self);
   

decrement

與 increment 類似，特別的是不用考慮借位不夠借的問題，因為本 bit stream 直接模擬二補數的行為。

// implementation of public method "decrement"
const int lim = bn_up_index((*self)->capacity);
_Bool borrow = 1;
for (int i = 0; borrow && i < lim; ++i) {
    borrow = (*self)->arr[i] == 0;
    --(*self)->arr[i];
}
// it's ok for twos complement that borrow is now true
refresh(self);

加法流程變化

流程與 increment 很相似，但改為提前檢測加法前後是否隱含擴展。

考慮二進制加法，只有當被加數與加數位元數相等時才必然進位。原因說明如下:

1. 說明最高位相等時會進位
   考慮 edge case: 僅最高位為一即可。
   ​​​​    1000...0000
   ​​​​+)  1000...0000
   ​​​​---------------
   ​​​​   1000...00000
   

   最高位相等的條件下，最難進位的兩數相加，即乘二，結果是左移一位。
2. 說明最高位相等最多進一位
   考慮 edge case: 全一即可。
   ​​​​    1111...1111
   ​​​​+)  1111...1111
   ​​​​---------------
   ​​​​   1111...11110
   

   最高位相等的條件下，最可能進位的兩數相加，即乘二，結果是左移一位。

利用 capacity 和 use 成員，是否需要擴展可以快速決定，故將隱含的擴展提至加法運算前。

加法技巧

在加法運算的迴圈中可以看到如下的技巧，當中 mask 為 0x7fffffffffffffff，此技巧是為了擷取「進位」。

// inside implementation of "add_in_place"
_Bool carry = 0;
const bn_data_t mask = ~(1ULL << __ULL_OFFSET);
for (int i = low_b_ull; i < use_b_ull; ++i) {  // add b
    // addition with carry bit detection
    bn_data_t tmp =
        ((*a)->arr[i] & mask) + ((*b)->arr[i] & mask) + carry;
    // highest two bits
    bn_data_t msb2 = ((*a)->arr[i] >> __ULL_OFFSET) +
                     ((*b)->arr[i] >> __ULL_OFFSET) + (tmp >> __ULL_OFFSET);
    (*a)->arr[i] = (tmp & mask) | ((msb2 & 1) << __ULL_OFFSET);
    carry = msb2 >> 1;
}

針對每組 bn_data_t，考慮加法流程變化對加法的簡單分析，得知相同的最高位決定「兩個」二進制數是否進位。雖然加法過程中有另一 carry 的存在，不過 carry 最高為 1，不影響前面的推論。

我使用 bitwise 操作過濾最高位元，tmp 即忽略最高位元的加法，其值不溢位。

因為考慮極端情況:

    0111...1111
    0111...1111
+)            1
---------------
    111...11111

第九行針對最高為進行加法後，以 or 運算將低 63 位和最高位 (msb2 之最低位) 的加法結果合併；新的 carry 便是 msb2 次低位，利用左移方能適應次低位為 0 或 1 兩種情況。

Minus: additive inverse

由於本 bit stream 實現的是二補數的運算，那麼便可以利用對某元素的減法為加上加法反元素的特性完成減法實作，為此衍伸出 minus 系列方法。

實作上與正常的加法器有點像，由低位起將一個 bn_data_t 內容 flip 與加一；另外 use 和 low 兩成員難以事先確認，故在最後需呼叫 refresh。

// implementation of public method "minus_in_place"
const int cap_ull = bn_up_index((*a)->capacity);
_Bool carry = 1;  // default add one
for (int i = 0; i < cap_ull; ++i) {
    // negate and add one
    (*a)->arr[i] = ((*a)->arr[i] ^ __ULL_ONES) + carry;
    // check carry up
    carry = !(*a)->arr[i];
}
refresh(a);

特別的是，由於 minus 預設是提供給減法使用，故要考慮有號長度擴展的議題。對此我實作 bn_minus 完成擴展邏輯。說是如此，也只需要 bn_copy 一下即可，因為 minus_in_place 會照顧到所有位元。

// implementation of private method "bn_minus"
struct big_num **res = bn_copy(a, expand_to);
minus_in_place(res);
return res;

減法

有了 add_in_place 和 minus，減法變得非常簡單：

// implementation of public method "sub"
const int use_a_ull = bn_up_index((*a)->use),
          use_b_ull = bn_up_index((*b)->use),
          expand_to = use_a_ull > use_b_ull ? use_a_ull : use_b_ull;
struct big_num **res = bn_minus(b, expand_to);
add_in_place(res, a);
return res;

TODO: bound 與否的 bug

利用二補數的性質，加減法目前看似天下太平，但其實迴避了一個嚴重的問題: 當前 bit stream 究竟「是不是負數」。
假設 bit stream 的內容為 0xffffffffffffffff，其是表示的是正還是負數無從得知。未來考慮兩種改法:

1. 加入 less_than, greater_than 方法，讓使用者在加減法前確認誰大誰小，以預測加、減法結果

   實際上這不解決問題。雖然利用二補數的性質，負數彼此間的次序跟正數彼此間的次序的比較無二異。問題是正數與負數間的比較要如何正確的完成 0.0

2. 概念上加入新成員 sign，在所有運算中維護此數
   考慮利用 Linux 中紅黑樹的節點著色技巧，由於 capacity 成員最低六個位元實際上沒有用到，可以在那裡紀錄
   $6$ bits 內的資訊。

乘法

乘法是目前最複雜的實作，雖然沒有用上特別精妙的演算法，不過還是用上了資工仔常用的 DP 跟 lazy evaluation。概念上源自最簡單的實作：見乘數某位元為

$1$ 時，將被乘數複製並位移，與回傳結果相累加，worse case 複雜度為

$O(mn)$，

$m$、

$n$ 分別為被乘數和乘數的 bn_data_t 數。

這樣的實作對於長度很大且

$1$ 很多的乘數來說極其低效，同時可預期乘數內不同

$1$ 數量的乘法運行時間大相徑庭，效能難以預期。仔細想會發現，若將乘數切成等長而更細的 bit slice，會發現程式碼存在大量重複的加法，比如下例。假設乘數如下 (方便起見，以一個 byte 為例)：

若兩兩為一組看，變成

如果我們紀錄被乘數乘以各組 bit slice 的結果，便可重複利用乘法、加法的計算結果，將問題切成重複的子問題，那麼就可以掏出 DP 來用。

首先先粗暴的建立表格，大小

測試

其實我覺得測試是最困難的。以前沒有獨自寫過測試程式，所以實作完只能肉眼除錯，或者低效的嘗試一些測試資料，真沒辦法就只能占卜了 0.0
Shiritai

這是 SAT 和 SMT Solving 的範疇，可見 SAT and SMT Solving

Wow, SAT 是演算法課程後面那個只簡單提到的東西…這得花時間好好讀 0.0 Shiritai

測試框架

我撰寫一個簡單的測試框架，以下讓我介紹方法。

Basic

除錯最簡單的即印出資料，show_hex 方法便由此而生，是第一個完成方法。一開始的測試方法非常粗暴:

// in test.c
int main()
{
    inf_bits bn1 = bn_methods.construct(0xffffffff);
    bn_methods.show_hex(bn1);
    bn_inner_methods.left_shift(bn1, 3);
    bn_methods.show_hex(bn1);
    bn_methods.mark(bn1, pos);
    bn_methods.show_hex(bn1);
    printf("%d\n", bn_methods.test(bn1, pos));
    /* ... */
}

這種除錯方法顯然過於低效，如果有專屬 inf_bits 物件的「相等檢查」將能簡化除錯流程，equals 方法因此誕生，該實作盡可能利用 use、low 成員縮減比較範圍。

// main implemantation part of public method "equals"
const int use_ull = bn_up_index((*self)->use),
          low_ull = bn_dn_index((*self)->low);
for (int i = low_ull; i < use_ull; ++i)
    if ((*self)->arr[i] != (*other)->arr[i])
        return 0;

這樣至少就能在不相等時才印出資料，讓除錯訊息不會過於冗長。

以巨集實作友善的 Assert

現在想想，C 語言的巨集和 Preprocessor 真的好方便，能獲得完全不同的寫扣體驗 Shiritai

測試時常用斷言來實現，比如 Rust 語言中執行 cargo test，通常會以有沒有觸發 assert_eq 等來判斷是否運行正確。C 語言有 assert 可用，不過比起直接引發 panic，我希望能印出錯誤訊息並繼續執行其他測試。故我將 equals 利用巨集包裝，其包含對人很友善的錯誤訊息。用法如下：

// in test.c::sub_test
inf_bits sub = bn_methods.sub(added, a);
show_and_abort_if_not_equal(sub, a);

實作如下：

#define show_and_abort_if_not_equal(bn, correct)                         \
    if (!bn_methods.equals(bn, correct)) {                               \
        printf("\n<Assertion failed at line %d in file %s>\n", __LINE__, \
               __FILE__);                                                \
        printf("Incorrect value: ");                                     \
        bn_methods.show_hex(bn);                                         \
        printf("Ground true: ");                                         \
        bn_methods.show_hex(correct);                                    \
        printf("Abort testing.\n\n");                                    \
        return 0;                                                        \
    }

利用 Standard Predefined Macros，印出的訊息包含檔案與行數等編譯時期可確定的訊息，錯誤訊息舉例如下:

<Assertion failed at line 468 in file test.c>
Incorrect value: [Metadata] addr: 0x6093b80 cap: 128, use 65, low: 1
[Data] addr: 0x6093b90, arr: 0000000000000001_4a6d6986c7b0b3f2
Ground true: [Metadata] addr: 0x659b7e0 cap: 128, use 63, low: 1
[Data] addr: 0x659b7f0, arr: 0000000000000000_4a6d6986c7b0b3f2
Abort testing.

以巨集裝飾測試函式

有了前述巨集，撰寫測試函式變的簡單許多，於是測試函式們一一登場。我另外撰寫 test_all 函式來呼叫這些測試函式。

_Bool test_all()
{
    if (!shift_test()) {
        printf("shift_test failed");
        return 0;
    }
    if (!single_bit_test()) {
        printf("single_bit_test failed");
        return 0;
    }
    /* ... */
}

上述程式碼怎麼看怎麼冗余，有沒有更好的使用方法？由於 failed 前的字串即函式名，當然可以選擇為函式包一層函式 (i.e. decorator)，其執行函式本身並印出函式名，這要求該函式需同時傳遞函式與函式名。此時若改用巨集的字串串接可以寫得更加簡潔：

#define run_test(test_name)                           \
    {                                                 \
        printf("INFO: "#test_name"\n");               \
        if (!test_name()){                            \
            printf("ERROR: " #test_name " failed\n"); \
            return 0;                                 \
        }                                             \
        printf("INFO: passed "#test_name" \n");       \
    }

訊息的效果如下：

...
INFO: passed single_bit_test 
INFO: copy_test
INFO: passed copy_test 
INFO: sub_test

<Assertion failed at line 345 in file test.c>
Incorrect value: [Metadata] addr: 0x6520800 cap: 128, use 62, low: 0
[Data] addr: 0x6520810, arr: 0000000000000000_369f8b3100f3f8d9
Ground true: [Metadata] addr: 0x6143f90 cap: 128, use 0, low: 0
[Data] addr: 0x6143fa0, arr: 0000000000000000_0000000000000000
Abort testing.

ERROR: sub_test failed

更進一步的，某些測試函式我想透過 Python 協助驗證，這要求執行的函式多一個檔名參數。利用巨集 __VA_ARGS__ 同樣可以完成要求，且為了避免重複寫近乎相同的巨集，我將前面的巨集拆分為以下奇妙的樣子：

#define __run_test_front(test_name)       \
    {                                     \
        printf("INFO: " #test_name "\n"); \
        if (!

#define __run_test_mid_no_arg(test_name) test_name()
#define __run_test_mid_args(test_name, ...) test_name(__VA_ARGS__)

#define __run_test_end(test_name)                 \
        )                                         \
    {                                             \
        printf("ERROR: " #test_name " failed\n"); \
        return 0;                                 \
    }                                             \
    printf("INFO: passed " #test_name " \n");     \
    }  

如此方能以巨集生成巨集：需參數與不需參數版的 run_test。

#define run_test(test_name)                                      \
    __run_test_front(test_name) __run_test_mid_no_arg(test_name) \
        __run_test_end(test_name)

#define run_test_args(test_name, ...)                                       \
    __run_test_front(test_name) __run_test_mid_args(test_name, __VA_ARGS__) \
        __run_test_end(test_name)

經這一番折騰，功能的測試程式碼變得極其簡單美觀，人類可輕鬆閱讀並擴充：）

_Bool test_all()
{
    run_test(shift_test);
    run_test(single_bit_test);
    ...
    // test that need to dump file
    run_test_args(dump_fib, "output/fib_output.txt");
    run_test_args(inc_dec_test, "output/inc_dec_output.txt");
    ...
    return 1;
}

Makefile 自動化測試

我準備利用三種工具測試程式

物件公開與私有方法測試

測試進度

TODO

依賴函式庫議題

初步解法

可以發現資料的型別為 bn_data_t，原本我將其定義為依賴 stdint.h 的:

#include <stdint.h>
/**
 * Define inner data type
 */
#define bn_data_t uint64_t

但編譯時會產生以下錯誤。

make -C /lib/modules/5.19.0-38-generic/build M=/home/eroiko/repos/linux2023/fibdrv modules
...
In file included from /home/eroiko/repos/linux2023/fibdrv/fibdrv.c:9:
/home/eroiko/repos/linux2023/fibdrv/util.h:69:10: fatal error: stdint.h: No such file or directory
   69 | #include <stdint.h>
      |          ^~~~~~~~~~
compilation terminated.
...
make: *** [Makefile:18: run] Error 2

也就是找不到 stdint.h。這是因為 fibdrv 本身即 kernel module，我們的 include path 中也不存在 C 語言標準函式庫，自然無法編譯。

參照 khttpd，利用 compat 目錄提供的 compatibility layer，解決程式碼在 Linux 核心和 non-freestanding 環境編譯問題。

jserv

參考本文，提到 Linux kernel 其實也定義好諸多型別可供使用，在 asm-generic/int-ll64.h 中。由於我希望此函式庫可同時適用於 kernel 和 user space，最終我利用條件編譯完成在兩種不同環境下的編譯。

#ifdef __KERNEL__

#include <asm-generic/int-ll64.h>
/**
 * Define inner data type
 */
#define bn_data_t u64

#else

#include <stdint.h>
/**
 * Define inner data type
 */
#define bn_data_t uint64_t

#endif

同樣的問題也適用在其他地方，比如配置記憶體、列印資料等原本依賴 C 語言標準函式庫的功能，其在 kernel 內不存在。此時可以利用 lab0-c 中提及的 function hooking 技巧，將 kernel 中對應 C 標準函式庫功能的函式重新包裝，搭配條件編譯便可在 user space 或 kernel 中編譯本函式庫，以下擷取用例。

// in big_num.c
// conditional compilation
#ifdef __KERNEL__

#include <linux/printk.h>
#include <linux/slab.h>

#define __bn_print(fmt, ...) printk(KERN_ALERT fmt, __VA_ARGS__)
#define __bn_malloc(size) kmalloc(size, GFP_KERNEL)
#define __bn_free(ptr) kfree(ptr)

#else

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#define __bn_print(fmt, ...) printf(fmt, __VA_ARGS__)
#define __bn_malloc(size) malloc(size)
#define __bn_free(ptr) free(ptr)

#endif

不過事情沒有我想的簡單，遇到下面的情況目前我無法解決。

我在 big_num.c 中實作 dump 函式，其依賴 fopen, fclose 等 File R/W 標準函式，本討論串提供看似開箱即用的解法，另外此文詳細說明這些 API 的行為，於是我嘗試翻閱原始碼，但…

恩，不是我眼睛業障中吧？！ Shiritai

對，並沒有，但很抱歉那是以前有的 0.0

恭喜你發現 Linux 核心的發展「遺跡」

jserv

可見與核心相關的函式隨版本的更迭可能有所變動，這限制本程式碼於不同 kernel 下編譯的成功與否，故我先在此打住，僅利用條件編譯略過 dump 函式，即 kernel space 中的本函式庫無 inf_bits::dump 的定義與實作。