2020q1 Homework2 (fibdrv)

contributed by < AndybnACT >

tags: linux2020

大數運算

第一次執行的時候發現費氏數列的計算結果不如預期。

Passed [-]
f(93) fail
input: 7540113804746346429
expected: 12200160415121876738

推測是因為整數相加溢位的緣故，打開輸出檔案 out 可以確認這個現象：

...
Reading from /dev/fibonacci at offset 91, returned the sequence 4660046610375530309.
Reading from /dev/fibonacci at offset 92, returned the sequence 7540113804746346429.
Reading from /dev/fibonacci at offset 93, returned the sequence 7540113804746346429.
...

理論上

$F(93)=F(91)+F(92)$ ；然而這邊輸出竟然和前一步一樣。打開 fibdrv.c 之後可以發現程式碼為了避免溢位，透過 #define MAX_LENGTH 92 將最大的可輸出範圍限制在

$k=92$。為了輸出後面的費氏數列，我們必須實作大數運算。這次作業的實作先嘗試將兩個無號長整數接在一起來表示一個 128 位元的無號大整數。

typedef struct __attribute__((packed)) uint128 {
    unsigned long long upper;
    unsigned long long lower;
} uint128_t;

gcc 有 128 位元整數的擴充特徵，可見 128-bit Integers:

• 有號數: __int128
• 無號數: unsigned __int128

另外可參照 recipes/basic/int128.h 得知 x86_64 對於 128-bit 整數運算相關的 inline assembly。

jserv

因為 read() 和 write() 回傳值型態 ssize_t 在 LP64 之下只有 64 位元，如果 client.c 還是用一樣的介面讀取資料的話，將無法透過一次 read() 拿到費氏數列的值。所以，我們改變 fibdrv 的介面，用 read() 的第二個參數、在 kernel 中透過 copy_to_user 將計算結果回傳給 user space 程式。此外，為了之後測量時間方便，read() 的回傳值為費氏數列函數的計算時間。

/* calculate the fibonacci number at given offset */
static ssize_t fib_read(struct file *file,
                        char *buf,
                        size_t size,
                        loff_t *offset)
{
    u64 t1, t2;
    t1 = ktime_get_ns();
    uint128_t result = fib_func(*offset);
    t2 = ktime_get_ns();
    copy_to_user(buf, &result, sizeof(uint128_t));
    return t2 - t1;
}

變更程式介面之後，client.c 也要做對應的修改，這邊為了省去進制轉換的麻煩，直接以兩個 %llx 做輸出。然後，感恩讚嘆 python ，讓 verify.py 讀懂這樣的輸出只需要做很簡單的修改。

• client.c：

​    for (int i = 0; i <= offset; i++) {
​       lseek(fd, i, SEEK_SET);
​       sz = read(fd, fib, 1);
​       printf("Reading from " FIB_DEV
​              " at offset %d, returned the sequence "
​              "0x%llx%016llx.\n",
​              i, fib[0], fib[1]);
​   }

• verify.py：

​f0 = int(result_split[-1][9].split('.')[0], 16)

add128

兩個大整數的相加可以先拆成高位與低位，然後再處理低位相加的 overflow 問題。其中，判斷式 ~op2.lower 代表 op2.lower 與該無號整數上界的距離。所以若 op1.lower 大於該距離，則代表相加會導致溢位。

static inline void add128(uint128_t *out, uint128_t op1, uint128_t op2)
{
    out->upper = op1.upper + op2.upper;
    if (op1.lower > ~op2.lower) {
        out->upper++;
    }
    out->lower = op1.lower + op2.lower;
}

完成之後，就可以來驗證費氏數列的輸出。經過測試，兩個無號長整數組合而成 128 位元的無號整數可以正確地輸出費氏數列直到

$k=187$。

sub128

同理，兩個大整數的相減可以先拆成高位與低位，然後再處理低位相減的 underflow 問題。當被減數小於減數時，處理 underflow 的方法是將輸出的高位再減去 1。

static inline void sub128(uint128_t *out, uint128_t op1, uint128_t op2)
{
    out->upper = op1.upper - op2.upper;
    if (op1.lower < op2.lower) {
        out->upper--;
    }
    out->lower = op1.lower - op2.lower;
}

lsft128

左移的方法是將高位低位都執行邏輯左移（logical shift），然後把低位應該移入高位的位元補上。需要注意的是 C 語言在處理 shift operand 大於等於變數的寬度時，行為未定義。

詳請見 C99 spec

If the value of the right operand is negative or is greater than or equal to the width of the promoted left operand, the behavior is undefined.

static inline void lsft128(uint128_t *out, uint128_t op, unsigned char shift)
{
    if (unlikely(shift == 0)) {
        *out = op;
        return;
    }
    if (shift >= 64) {
        out->upper = op.lower << (shift - 64);
        out->lower = 0ull;
        return;
    }
    out->upper = op.upper << shift;
    out->lower = op.lower << shift;
    out->upper |= op.lower >> (64 - shift);
}

考慮以下變更:

@@ -1,9 +1,11 @@
-static inline void lsft128(uint128_t *out, uint128_t op, unsigned char shift)
+#define unlikely(x) __builtin_expect((x), 0)
+static inline void lsft128(uint128_t *out, uint128_t op, uint8_t shift)
 {
-    if (shift == 0) {
+    if (unlikely(shift == 0)) {
         *out = op;
         return;
-    } else if (shift >= 64) {
+    }	
+    if (shift >= 64) {
         out->upper = op.lower << (shift - 64);
         out->lower = 0ull;
         return;

要點:

1. 使用 uint8_t 表達 128 位元的範圍;
2. 引入 unlikely 巨集，提示編譯器產生對 branch predictor 友善的程式碼;
3. 略過多餘的 else，而且程式縮排更精簡

jserv

已經改善，見 commit 0ade862

rsft128

同理，右移之前，先判斷移動量。若移動量為零則直接將輸入複製到輸出；右移移動量大於變數寬度時，直接將輸出的高位清空，然後低位即等於原本高位右移 shift - 64 的值。

static inline void rsft128(uint128_t *out, uint128_t op, unsigned char shift)
{
    if (unlikely(shift == 0)) {
        *out = op;
        return;
    }
    if (shift >= 64) {
        out->upper = 0ull;
        out->lower = op.upper >> (shift - 64);
        return;
    }
    out->upper = op.upper >> shift;
    // signed bit will not extend for unsigned type
    out->lower = op.lower >> shift;
    out->lower |= op.upper << (64 - shift);
}

mul128

TODO: 解釋 popcountll 的作用和考量點

jserv

這邊實作乘法的演算法是最簡單的二進位直式乘法；把所有的乘法換成加法和左移的組合。因為已經實作 128 位元的加法和左移，所以直接呼叫就可以。函式前半段用 popcountll 來抓出 1 位元比較少的數當做乘法直式中的乘數，以減少加法運算的次數。但是 gcc(9.2) 在編譯 kernel module、最後 linking 階段找不到這個函式。所以先暫時把 op1 和 op2 直接指定為 lo 和 hi。

• popcountll 的作用及考量點：popcount 的作用是算出一個輸入變數中，有幾個位元是一。例如：假設輸入為 b101011 ，則 popcount 就會回傳 4。而其考量點可以先考慮下列兩個二進位直式乘法

                b 11111011                        b 00100100
             x  b 00100100                    x   b 11111011
            ----------------                  ----------------
                 11111011                           00100100
             11111011                              00100100
                                                 00100100
                                                00100100
                                                ...

因乘法交換率的緣故，兩邊的結果是相同的。然而右邊的被乘數含有 1 的位元明顯大於乘數，使得要重複執行許多次加法才能夠獲得和左邊只使用一次加法就能得到的結果。

static inline void mul128(uint128_t *out, uint128_t op1, uint128_t op2)
{
    int op1cnt =
        __builtin_popcountll(op1.lower) + __builtin_popcountll(op1.upper);
    int op2cnt =
        __builtin_popcountll(op2.lower) + __builtin_popcountll(op2.upper);
    uint128_t lo, hi;
    if (op1cnt > op2cnt) {
        lo = op2;
        hi = op1;
    } else {
        lo = op1;
        hi = op2;
    }
    out->lower = 0ull;
    out->upper = 0ull;
    while (lo.lower) {
        uint128_t prod;
        unsigned char ctz = __builtin_ctzll(lo.lower);
        lsft128(&prod, hi, ctz);
        add128(out, prod, *out);
        lo.lower &= ~(1ull << ctz);
        // printf("lower: %llx, ctz=%d\n", lo.lower, ctz);
    }
    while (lo.upper) {
        uint128_t prod;
        unsigned char ctz = __builtin_ctzll(lo.upper);
        lsft128(&prod, hi, ctz + 64);
        add128(out, prod, *out);
        lo.upper &= ~(1ull << ctz);
        // printf("upper: %llx, ctz=%d\n", lo.upper, ctz);
    }
}

實作完以上函式之後，就可以用來測試費氏數列的 Fast doubling (

$O(logn)$) 的演算法。

加法的效能改進

看完〈How to implement bignum arithmetic〉之後，得知透過使用更大的資料型態來儲存較小數值的相加結果，可以用來保存 carry。為了正確的實作這種加法概念，我們需要先將 struct uint128 的順序改成低位先、高位在後。

typedef struct __attribute__((packed)) uint128 {
    unsigned long long lower;
    unsigned long long upper;
} uint128_t;

加法的實作如下：

static inline void add128_auto_carry(uint128_t *out, uint128_t op1, uint128_t op2) {
    unsigned int *a = (unsigned int *) &op1;
    unsigned int *b = (unsigned int *) &op2;
    unsigned int *c = (unsigned int *) out;
    unsigned long long l = 0;
    
#pragma GCC unroll 4
    for (size_t i = 0; i < 4; i++) {
        l += (unsigned long long) a[i] + b[i];
        c[i] = (unsigned int) l;
        l >>= 32;
    }
}

但是相較於原本的做法，這樣做只省去了一個 branching、卻換來額外兩個 addition。除非 64 位元的機器處理 32 位元的加法有比較快、而且這些運算可以被放在 32 位元的加法當中，可能才會有效能改進。

初步查了一下 x86 各指令的 cycle count，發現 64 位元加法和 32 位元加法的花費是一樣的。

乘法的效能改進（Comba）

Comba multiplication
在該篇演講中，作者也有提到他實作大數乘法的方式。一樣也是直式乘法，然而他並沒有將乘法拆成加法，因為硬體本身還是有乘法器可以用，只是寬度不夠而已。解決寬度不夠的方法一樣也是用 64 位元的資料型態儲存兩個 32 位元的乘法結果。

相較於加法的情境，這麼運用顯得比較有效率，因為兩個 32 位元的整數相加最多只會變成 33 位元、儲存在 64 位元的空間中浪費了 31 位元；然而相乘的結果有機會變成 64 位元，減少浪費的情形。

src
Comba 演算法其實也是直式乘法，相較於一般的直式乘法是 row-based，一行一行將結果相加到結果中；Comba 採用 column-based，一列一列寫到結果裡面。這麼做在寫入時對 wb-cache 的壓力相對降低許多，因此可能有比較好的效能。

• 實作 comba square

為了比較 Comba 和一般的直式乘法，先實作（運用乘法單元的）直式乘法

static inline void mul128_naive(uint128_t *out, uint128_t lo, uint128_t hi)
{
    unsigned int *src1 = (unsigned int *) &lo;
    unsigned int *src2 = (unsigned int *) &hi;
    unsigned long long l = 0;
    unsigned int acc[8] = {0, };
    int shift = 0;

#pragma GCC unroll 4
    for (size_t i = 0; i < 4; i++) {
#pragma GCC unroll 4
        for (size_t j = 0; j < 4; j++) {
            shift = i + j;
            if (shift >= 4) {
                continue;
            }
            l = (unsigned long long) src1[i] * src2[j];
            *((unsigned long long*) (acc + shift)) += l;
        }
    }
    memcpy(out, acc, sizeof(uint128_t));
}

mul128_school 不是很好理解的名稱，可改為 mul128_naive (不是 native，而是 naive，後者的意思是「天真幼稚」)。

jserv

已經改善，見 commit 0ade862

而採用 Comba 演算法的實作如下：

static inline void mul128_comba(uint128_t *out, uint128_t lo, uint128_t hi)
{
    unsigned int *src1 = (unsigned int *) &lo;
    unsigned int *src2 = (unsigned int *) &hi;
    unsigned long long l = 0;
    unsigned int acc[8] = {0, };

#pragma GCC unroll 4
    for (int sft = 0; sft < 4; sft++) {
        for (int i = 0, j = sft; i <= sft; i++, j--) {
            l = (unsigned long long) src1[i] * src2[j];
            *((unsigned long long*) (acc + sft)) += l;
        }
    }
#pragma GCC unroll 4
    for (int sft = 4; sft < 8 - 1; sft++) {
        for (int i = sft - 3, j = 3; i < 4 ; i++) {
            l = (unsigned long long) src1[i] * src2[j];
            *((unsigned long long*) (acc + sft)) += l;
        }
    }
    memcpy(out, acc, sizeof(uint128_t));
}

實作完成之後，可以在 user-space 透過 ADDRESS_SANITIZER 檢查有沒有越界存取，編譯時選項需要加上 -fsanitize=address -fno-omit-frame-pointer -fno-common。

GCC 的實作

gcc 也有提供 128 位元整數的操作，[unsigned] __int128，透過實驗可以觀察其底層實作。實驗的方式為寫一段用到 __int128 運算的 C 程式碼，編譯成執行檔之後觀察執行結果的正確性，然後用反組譯工具 objdump -d -M intel 分析實作方式：

__int128 加法

int main () {
    unsigned __int128 x = 0xFFFFFFFFFFFFFFFFull;
    unsigned __int128 y = 0xDEADBEAF;
    unsigned __int128 z = 0x0;
    z = x + y;
    printf("size of __int128 = %zu\n", sizeof(x));
    printbig("z=", z);
    return 0;
}

看來 gcc 用符合預期的寬度來儲存 __int128 型態的變數、而且加法結果正確。

$ gcc __int128.c 
$ ./a.out
size of __int128 = 16
z= 0x000000000000000100000000deadbeae
$ objdump -d -M intel a.out
...
119a:	48 8b 4d c0    mov    rcx,QWORD PTR [rbp-0x40]  # ((char*)&y)
119e:	48 8b 5d c8    mov    rbx,QWORD PTR [rbp-0x38]  # ((char*)&y) + 8
11a2:	48 8b 45 d0    mov    rax,QWORD PTR [rbp-0x30]  # ((char*)&x)
11a6:	48 8b 55 d8    mov    rdx,QWORD PTR [rbp-0x28]  # ((char*)&x) + 8
11aa:	48 01 c8       add    rax,rcx
11ad:	48 11 da       adc    rdx,rbx
...

反組譯的結果可以看到一道特別的指令 adc ，根據 Intel® 64 and IA-32 Architectures
Software Developer’s Manual 顯示

Adds the destination operand (first operand), the source operand (second operand), and the carry (CF) flag and stores the result in the destination operand.
…
The ADC instruction does not distinguish between signed or unsigned operands. Instead, the processor evaluates the result for both data types and sets the OF and CF flags to indicate a carry in the signed or unsigned result, respectively. The SF flag indicates the sign of the signed result.

比較讓人注意到的是，adc 這個指令不用管輸入是無號還是有號整數。CPU 兩個都會做，然後再用不同的 FLAG 來存放有號、無號整數相加後的 carry-bit。回頭看看一般的 add 指令，發現其實也有這樣的敘述。多虧了 2 補數的表示法，讓有號、無號數運算方式是一致的，所以硬體不需要另外再做一次加法。

              b1111 1111                      b0111 1111
            + b0000 0001                    + b0000 0001
            --------------                -----------------
(unsigned)    b0000 0000 -> carry             b1000 0000 -> expected res
  (signed)    b0000 0000 -> expected res      b1000 0000 -> overflow
              ==> (CF/OF)=(1/0)               ==> (CF/OF)=(0/1)

由此也可以看出 gcc 對 128 位元變數的安排一樣也是先放低位再放高位（$rax/$rcx為低位，低位做完之後 adc 將低位的 carry 加進高位）。

__int128 位元左右移

為了測試左右移，我們嘗試對 x 左移 8 bits。實際測試和反組譯的結果如下：

$ ./a.out
size of __int128 = 16
z = x << 8 0x00000000000000ffffffffffffffff00
$ objdump -d -M intel  a.out
...
1196:	48 8b 45 d0      mov    rax,QWORD PTR [rbp-0x30] # ((char*)&x)
119a:	48 8b 55 d8      mov    rdx,QWORD PTR [rbp-0x28] # ((char*)&x) + 8
119e:	48 0f a4 c2 08   shld   rdx,rax,0x8
11a3:	48 c1 e0 08      shl    rax,0x8
...

又發現一道有趣的指令 shld，而 Intel® 64 and IA-32 Architectures
Software Developer’s Manual 描述：

The SHLD instruction is used for multi-precision shifts of 64 bits or more.
The instruction shifts the first operand (destination operand) to the left the number of bits specified by the third operand (count operand). The second operand (source operand) provides bits to shift in from the right (starting with bit 0 of the destination operand).

這告訴我們，shld 會把目標暫存器左移指定的位元數、而且還會幫我們把來源暫存器的高位也一併移入目標暫存器中。需要注意移動位元數不得超過暫存器寬度。
此外，實驗發現在左移超過 64 位元時，gcc 產生的程式碼與我們的實作類似，將低位寫入高位之後，左移 shift - 64 位元，然後再將低位清空。

__int128 減法

在查詢指令的時候有看到 sbb，似乎也是與 carry 有關的減法指令。實際執行測試時，令 x = 0, y = 1，結果如下：

$ ./a.out
size of __int128 = 16
z= 0xffffffffffffffffffffffffffffffff
$ objdump -d -M intel ./a.out
...
1160:	48 c7 45 d0 00 00 00 00	mov    QWORD PTR [rbp-0x30],0x0 # ((char*)&x)
1168:	48 c7 45 d8 00 00 00 00	mov    QWORD PTR [rbp-0x28],0x0 # ((char*)&x) + 8
1170:	48 c7 45 e0 01 00 00 00	mov    QWORD PTR [rbp-0x20],0x1 # ((char*)&y)
1178:	48 c7 45 e8 00 00 00 00	mov    QWORD PTR [rbp-0x18],0x0 # ((char*)&y) + 8
...
1190:	48 8b 45 d0    mov    rax,QWORD PTR [rbp-0x30]
1194:	48 8b 55 d8    mov    rdx,QWORD PTR [rbp-0x28]
1198:	48 2b 45 e0    sub    rax,QWORD PTR [rbp-0x20]
119c:	48 1b 55 e8    sbb    rdx,QWORD PTR [rbp-0x18]

而開發手冊是這麼寫的：

SBB—Integer Subtraction with Borrow:
Operation:
DEST ← (DEST – (SRC + CF)

因此，在執行減法的時候，低位先相減，然後用 sbb 減去高位。低位相減產生的 carry 會拿去高位減。

__int128 乘法

乘法的測試中，我們令 x 為無號 64 位元整數的最大值，y 則是 0xDEADBEEF，兩者相乘。結果如下：

$ ./a.out
size of __int128 = 16
z= 0x00000000deadbeaeffffffff21524151
$ objdump -d -M intel a.out
...
1199:	48 8b 45 d8     mov    rax,QWORD PTR [rbp-0x28] # x.hi
119d:	48 0f af 45 e0  imul   rax,QWORD PTR [rbp-0x20] # rax = (x.hi * y.lo)[0:63]
11a2:	48 89 c2        mov    rdx,rax                  # rdx = rax
11a5:	48 8b 45 e8     mov    rax,QWORD PTR [rbp-0x18] # rax = y.hi
11a9:	48 0f af 45 d0  imul   rax,QWORD PTR [rbp-0x30] # rax = (y.hi * x.lo)[0:63]
11ae:	48 8d 0c 02     lea    rcx,[rdx+rax*1]          # rcx = rax + rdx
11b2:	48 8b 45 e0     mov    rax,QWORD PTR [rbp-0x20] # rax = y.lo
11b6:	48 f7 65 d0     mul    QWORD PTR [rbp-0x30]     # rdx:rax = (y.lo * x.lo)[0:127]
11ba:	48 01 d1        add    rcx,rdx
11bd:	48 89 ca        mov    rdx,rcx

經查發現其實 mul 指令在 operand 大小為 64 位元的資料時，會用 rdx:rax 的組合來儲存 128 位元的乘法結果。
這邊的乘法實作方式是將 128 位元的整數切成兩個 64 位元的整數，然後用一般的直式乘法將兩數相乘。理論上兩個 128 位元的整數相乘應該會需要 256 位元的整數來儲存其結果。但是因為輸出（z）也只有 128 位元，所以這裡沒有做 x.hi * y.hi、也把x.lo * y.hi 以及 y.lo * x.hi 高於 128 位元的部分捨去。

對更大數的支援（Arbitrary Big Number）

• How to implement bignum arithmetic
• YouTube recording

吸收完 gcc 對 __int128 操作的美妙之後，決定用類似的方式實作真正任意大小的大數運算，姑且依標頭檔 abn.h 命名為 ABN。為了讓測試更方便，我們先在 user space 開發測試。

ABN 資料結構

一樣也是用陣列的方式將多個無號長整數組合起來，表達一個大的無號長整數「大數」。
其中，cnt 代表目前大數佔用了幾個無號長整數的空間；cap 代表目前的變數儲存空間允許置放多少個長整數；而 heap 及 stack 則代表儲存空間的實際位置。存放方式一樣是低位在前高位在後。

#define STACKALLOC 5
typedef struct bn {
    int cnt;
    int cap;
    uint64_t *heap;
    uint64_t stack[STACKALLOC];
} bn;

顧名思義，為了讓數字在相對小的時候減少記憶體配置器（allocator）的開銷，我們預先在 struct 裡面分配 uint64_t * STACKALLOC 的大小。待該空間已經不夠存放的時候，重新用記憶體配置器於 heap 中分配可儲存的空間。
如果在使用 stack 的時候也把儲存 heap 指標的空間拿來存放大數的話，空間運用會更有效率。但為了初期開發方便，這邊分開用兩個變數儲存。
之後，在實際執行運算時，我們需要一個能判斷大數存放在 heap 還是 stack 的方法。我們的做法是：先將 heap 初始化為 NULL，重新分配空間之後 heap 都是有值的變數。因此，取數只需要這麼做即可：

uint64_t *bn_getnum(bn *b)
{
    return b->heap ? b->heap : b->stack;
}

而重新分配空間的方式和 lib C realloc 的方式差不多，先配置指定大小的空間，將該空間初始化，然後帶入舊空間的值之後，將舊的空間釋放：

bn *bn_realloc(bn *orig, int newcap)
{
    dprintf(3, "%p bn_realloc %d => %d\n", orig, orig->cap, newcap);
    if (newcap > STACKALLOC) {
        uint64_t *src = bn_getnum(orig);
        uint64_t *heap = (uint64_t *) malloc(newcap * sizeof(uint64_t));
        if (!heap)
            return NULL;

        memset(heap, 0, newcap * sizeof(uint64_t));
        memcpy(heap, src, orig->cap * sizeof(uint64_t));
        bn_free(orig);

        orig->cap = newcap;
        orig->heap = heap;
    }
    return orig;
}

基本的資料配置、存取完成之後，就可以實作大數運算。所有大數運算的邏輯都是這樣：先判斷目標運算元的空間是否足夠，不足則配置空間；然後按照順序執行運算的同時，考慮長整數邊界的進位問題；最後目標寫入完成之後，更新目標運算元大數結構內的 cnt。

ABN 加法內部實作

開始進行加法之前，我們先透過 cnt 找出位元數比較少的數當作被加數。這樣一來，當低位都相加完之後，我們就只剩下低位加法產生的 carry，和高位要處理。而加法運算本身，為了直接拿到計算而產生的 carry，我們透過 inline assembly 來完成。

static inline uint8_t __add_ll(uint64_t *dst,
                               uint64_t src1,
                               uint64_t src2,
                               uint8_t c)
{
    uint8_t rc;
    uint64_t res;
    // Since carry is either 0 or 1, we only need at most 65 (but not 66) bits
    // to hold the result from src1 + src2 + carry. Even if both of them have
    // all bits set.
    // FIXME: should assert rc != 0 || rc != 1
    // FIXME: ugly asm
    __asm__(
        "xor %%rsi, %%rsi\n\t"
        "movb %4, %%sil\n\t"
        "movq %2, %1\n\t"
        "addq %3, %1\n\t"
        "setc %0\n\t"
        "addq %%rsi, %1\n\t"
        "setc %%cl\n\t"
        "orb %%cl, %0\n\t"
        : "=&r"(rc), "=&r"(res)
        : "r"(src1), "r"(src2), "r"(c)
        : "rsi", "cl");
    *dst = res;
    return rc;
}

不像 __int128 可以直接在緊鄰的高位運算用 adc 指令將 carry 帶入（除非展開迴圈），我們的 inline assembly 透過 setc 這個指令（set byte if CF == 1）把每次計算的 carry 帶回 C 程式碼，然後迴圈下一步再帶入上次得到的 carry。因為 carry 最多就是一，整個加法最多只會進位一次，所以可以一次把他們加起來。

    for (i = 0; i < mincap; i++) {
        // d[i] = s1[i] + s2[i] + c;// This CODE is not SAFE
        // if (s1[i] > ~s2[i]) {    // Consider s2 = u64_max, s1 = 0, c = 1
        //     c = 1;
        // } else {
        //     c = 0;
        // }
        c = __add_ll(d + i, s1[i], s2[i], c);
        dprintf(10, "d[%d] = 0x%lx, carry = %d\n", i, d[i], c);
    }
...

ABN 減法內部實作

減法的邏輯與加法一致，不過在減法中我們不需要找出比較小的數，而是直接將兩個運算原由低位置高位逐一相減。

static inline uint8_t __sub_ll(uint64_t *dst,
                               uint64_t src1,
                               uint64_t src2,
                               uint8_t c)
{
    uint8_t rc;
    uint64_t res;
    // do src1 - carry - src1
    __asm__(
        "subq %4, %2\n\t"
        "setc %0\n\t"
        "subq %3, %2\n\t"
        "setc %%sil\n\t"
        "orb %%sil, %0\n\t"
        "movq %2, %1\n\t"
        : "=&r"(rc), "=r"(res)
        : "r"(src1), "r"(src2), "r"((uint64_t) c)
        : "rsi");
    *dst = res;
    return rc;
}

在減法的過程中，因為不確定減法結果的有效位元為何，所以只要減出來的結果不為零，就更新有效位元 cnt。

...
    d = bn_getnum(dst);
    s1 = bn_getnum(src1);
    s2 = bn_getnum(src2);
    for (size_t i = 0; i < src1->cnt; i++) {
        borrow = __sub_ll(d + i, s1[i], s2[i], borrow);
        if (d[i])
            dst->cnt = i + 1;
    }

ABN 左右移內部實作

在不失一般性的前提下，這邊以左移為例。我們可以善用 shld 這個有趣的指令，在左移高位的同時，幫忙把低位補進來。但要記得的是左移之前，必須要先把高位暫存出來。不然在高位左移後的下一步，我們會找不到要補進更高位的那些位元。

    num = bn_getnum(dst);
    for (int i = 0; i <= dst->cnt; i++) {
        uint64_t tmp;
        uint64_t tmp2;
        // dprintf(10, "%llx %llx\n", tmp, num[i]);
        // should clobber memory here, or write to the
        // location elsewhere in C
        __asm__(
            "movq (%2), %0\n\t"
            "shldq %%cl, %1, (%2)"
            : "=&r"(tmp2)
            : "r"(tmp), "r"(num + i), "c"(amount)
            :);
        tmp = tmp2;
        // dprintf(10, "%llx %llx\n", tmp, num[i]);
    }

這段程式碼被以 -O2 以上編譯的時候，會有預期之外的行為，主要原因是：於除非在 : 區段標注清楚， gcc 實際上無從直接得知 inline assembly 裡面對暫存器、記憶體位置等的使用情形。因此，當 inline assembly 裡面用到輸入、輸出以外的儲存空間時（暫存器、記憶體），需要在 clobber list 中通知 gcc「不能假設」在 inline assembly 前後，這些「空間」的值是一致的。因此，暫時能解決問題，但不是很完善的方式為：在 clobber list 裡面加上 "memory"。「不是很完善」的原因是，告訴 gcc 整個記憶體在 inline assembly 操作之後都不能假設其一致性的成本似乎過大。

• 測試 clobber list 加上 memory vs 將結果寫回 C，再由 C 寫回記憶體的開銷 （-O2）

ABN 乘法內部實作

乘法沿用 Comba multiplication，不過在發現 mul 指令做 64 位元乘法時可以保留 128 位元的運算結果後，我們用 inline assembly 把這樣的結果拉回 C 程式語言使用。

static inline struct prod __mul_ll(uint64_t src1, uint64_t src2)
{
    struct prod ret;
    ret.lo = src1;
    // what a neat inline asm is it?
    // https://github.com/chenshuo/recipes/blob/master/basic/int128.h#L26
    __asm__("mulq %2\n\t" : "=d"(ret.hi), "+a"(ret.lo) : "r"(src2) :);
    return ret;
}

乘法和加法類似，先找出位元數（ cnt ）比較少的數值當作被乘數，寫在直式乘法下半，然後才由結果的低位開始計算，依序往高位移動、寫入。因為每次加回結果的數都是 64 位元，我們可以直接沿用 __add_ll() 來幫忙這項操作、並且兼顧 carry 的問題；同時，記得在乘法開始之前先將結果初始化為零。

/*                          var   | loop index |   boundary   |   width
 *                        ===============================================
 *     A3, A2, A1, A0 --> mpcand  |      j     |   *wcand_p   |  larger
 *   x B3, B2, B1, B0 --> mplier  |      i     |   *wlier_p   |  smaller
 *  ---------------------------------------------------------------------
 *                          d     |    index   |        dst->cnt
 */
...
    while (index < width - 1) {
        if (index < *wlier_p) {
            i = index;
            j = 0;
        } else {
            i = *wlier_p - 1;
            j = index - i;
        }

        while (i >= 0 && j < *wcand_p) {
            dprintf(10, "Doing A[%d] x B[%d]\n", j, i);
            if (mpcand[j] == 0 || mplier[i] == 0) {
                dprintf(10, "zero element --> skipping\n");
                goto skip;
            }
            carry = 0;
            prod = __mul_ll(mpcand[j], mplier[i]);
            carry = __add_ll(d + index, d[index], prod.lo, carry);
            carry = __add_ll(d + index + 1, d[index + 1], prod.hi, carry);
            if (carry) {
                d[index + 2]++;
            }
            if (d[index])
                dst->cnt = index + 1;
        skip:
            i--;
            j++;
        }
        dprintf(10, "\n");
        index++;
    }

ABN 用來計算費氏數列

#include <stdint.h>
#include <stdio.h>

#include "abn.h"

bn bn_fib_doubling(uint64_t k)
{
    bn a, b;
    bn t1, t2, bb, tmp, bsquare, asquare;
    int clz;
    int f_i = 0;

    bn_init(&a);
    bn_init(&b);
    bn_init(&t1);
    bn_init(&t2);
    bn_init(&bb);
    bn_init(&tmp);
    bn_init(&bsquare);
    bn_init(&asquare);

    bn_set(&a, 0);
    bn_set(&b, 1);
    if (k <= 1) {
        return k == 0 ? a : b;
    }

    clz = __builtin_clzll(k);
    for (int i = (64 - clz); i > 0; i--) {
        bn_assign(&bb, &b);
        __bn_shld(&bb, 1);
        bn_sub(&tmp, &bb, &a);

        bn_mul_comba(&t1, &a, &tmp);
        bn_mul_comba(&asquare, &a, &a);
        bn_mul_comba(&bsquare, &b, &b);
        bn_add(&t2, &asquare, &bsquare);

        bn_assign(&a, &t1);
        bn_assign(&b, &t2);
        if (k & (1ull << (i - 1))) {  // current bit == 1
            bn_add(&t1, &a, &b);
            bn_assign(&a, &b);
            bn_assign(&b, &t1);
            f_i = (f_i * 2) + 1;
        } else {
            f_i = f_i * 2;
        }
    }

    bn_free(&b);
    bn_free(&t1);
    bn_free(&t2);
    bn_free(&bb);
    bn_free(&tmp);
    bn_free(&bsquare);
    bn_free(&asquare);
    return a;
}

int main(int argc, char const *argv[])
{
    bn f30000;
    f30000 = bn_fib_doubling(30000);
    bn_print(&f30000);
    bn_free(&f30000);
    return 0;
}

sysprog21/bignum 的大數實作

費氏數列

clz/ctz 的運用

在計算費氏數列時，若採用 Fast Doubling 的演算法，可以較快速地（

$O(logn)$）迭代至要求目標

$k$。原因在於這種演算法將等式左右的關係矩陣化，然後找出一個

$k,k+1$ 對應到

$2k,2k+1$ 的關係。而二進位的數字又很容易用上這種關係，舉例來說，

$k=13$ 的二進位表達式為 b1101，也就是

b1101 = (((1*2)+1)*2+0)*2+1
           ^    ^    ^    ^.... bit0 (LSB)
           |    |    |......... bit1
           |    |.............. bit2
           |................... bit3 (MSB)

這樣一來，我們只需要從

$k$ 的 MSB 開始往右迭代，每次都計算

$F(2k)$ 和

$F(2k+1)$，若遇到位元是一的話只要再算

$F(2k+2)$，然後繼續迭代。而 clz 是硬體提供的指令，count leading zeros，可以讓我們快速地找到

$k$ 裡面，MSB 所在的位置。

演算法實作

• fib_sequence_ll：若計算結果可以保持在 64 位元的無號長整數內的話，一般連加演算法可以直接從定義寫出來：

    f[0] = 0;
    f[1] = 1;

    for (int i = 2; i <= k; i++) {
        f[i] = f[i - 1] + f[i - 2];
    }

• fib_sequence：若考慮大數，以我們目前的實作可以延長到 128 位元。因為是自定義的型態，呼叫函式做加法運算：

    f[0] = (uint128_t){.upper = 0, .lower = 0};
    f[1] = (uint128_t){.upper = 0, .lower = 1};

    for (int i = 2; i <= k; i++) {
        add128(f + i, f[i - 1], f[i - 2]);
    }

• fibseq_doubling_ll：Fast Doubling 快速演算法，若使用內建型態只需要直接使用內建運算子。需要注意的是 __builtin_clzll() 在輸入為零的時候，行為未定義，所以需要預先處理：

    int clz = __builtin_clzll(k);
    for (int i = (64 - clz); i > 0; i--) {
        u64 t1, t2;
        t1 = a * ((b << 1) - a);
        t2 = a * a + b * b;

        a = t1;
        b = t2;
        if (k & (1ull << (i - 1))) {  // current bit == 1
            t1 = a + b;
            a = b;
            b = t1;
        }
    }
    ret.lower = a;

• fibseq_doubling： Fast Doubling 考慮大數的實作如下：

    int clz = __builtin_clzll(k);
    for (int i = (64 - clz); i > 0; i--) {
        uint128_t t1, t2, tmp;
        lsft128(&tmp, b, 1);   // tmp = 2b
        sub128(&tmp, tmp, a);  // tmp = tmp -a
        mul128(&t1, a, tmp);   // t1 = a*tmp
        mul128(&a, a, a);   // a = a^2
        mul128(&b, b, b);   // b = b^2
        add128(&t2, a, b);  // t2 = a^2 + b^2
        a = t1;
        b = t2;
        if (k & (1ull << (i - 1))) {  // current bit == 1
            add128(&t1, a, b);
            a = b;
            b = t1;
        }
    }

效能分析

實驗環境

根據作業描述提到的方式設定實驗環境。

• Operating System:

$ cat /proc/version
Linux version 5.5.7-100.fc30.x86_64 ([e-mail]) (gcc version 9.2.1 20190827 (Red Hat 9.2.1-1) (GCC)) #1 SMP Fri Feb 28 17:32:51 UTC 2020

• Toolchain:

$ gcc --version
gcc (GCC) 9.2.1 20190827 (Red Hat 9.2.1-1)

• Additional Linux Boot parameters: isolcpus=0
• CPU: Intel® Core™ i5-2500 CPU @ 3.30GHz
• CPU Turbo Mode: Off

$ sudo sh -c "echo 1 > /sys/devices/system/cpu/intel_pstate/no_turbo"

• CPU Scaling Governor: Performance

#!/bin/bash
for i in /sys/devices/system/cpu/cpu*/cpufreq/scaling_governor
do
    echo performance > $i
done

• ASLR: Off

$ sudo sh -c "echo 0 > /proc/sys/kernel/randomize_va_space"

• Command Line:

$ make load
$ sudo taskset 0x1 ./client 
$ make unload

為了方便切換實驗環境，將上述設定寫在一個 script 裡面：

#!/bin/bash
make unload
make load

CPUID=0
ORIG_SCL=`cat /sys/devices/system/cpu/cpu$CPUID/cpufreq/scaling_governor`
ORIG_NTURBO=`cat /sys/devices/system/cpu/intel_pstate/no_turbo`
ORIG_ASLR=`cat /proc/sys/kernel/randomize_va_space`

sudo sh -c "echo 1 > /sys/devices/system/cpu/intel_pstate/no_turbo"
sudo sh -c "echo performance > /sys/devices/system/cpu/cpu$CPUID/cpufreq/scaling_governor"
sudo sh -c "echo 0 > /proc/sys/kernel/randomize_va_space"

rm fib*.dat
sudo taskset 0x1 ./client

make unload

sudo sh -c "echo $ORIG_NTURBO > /sys/devices/system/cpu/intel_pstate/no_turbo"
sudo sh -c "echo $ORIG_SCL > /sys/devices/system/cpu/cpu$CPUID/cpufreq/scaling_governor"
sudo sh -c "echo $ORIG_ASLR >  /proc/sys/kernel/randomize_va_space"

因為執行時期只有 client 一個 process，所以只有將一個 CPU 的組態設定為 performance。

計算效能

在 fibdrv.c 裡面，我們透過 write() 介面變更計算費波那契數的函式。目前實作的四種函式一使用的資料型態和演算法分類如下，其中，目前自訂 uint128_t 的算數用前述（add128, mul128）的方式實作。

	自訂 uint128_t	內建 u64
一般演算法（Regular）	fib_sequence()	fib_sequence_ll()
快速演算法（Fast Doubling）	fibseq_doubling()	fibseq_doubling_ll()

• 比較上下圖，可以看出自己實作大數加法的執行時間大約是內建加法的兩倍。可以理解原因是 uint128_t 加法中，會用到兩次加法和一次 branching；所以也代表 branch overhead 其實不高。
  因此，從這邊也可以推測出使用小變數、自動進位的方式可能帶來的進步有限。此外，未來如果可以有效地記錄最大位元的位置，則可以省去不必要的加法，帶來較顯著的效能成長。

聚焦在下圖，可以看到 Fast doubling 的計算效能確實較一般演算法來得好、也大致呈現

$O(logn)$ 的成長。然而，在上圖中，即便迭代次數明顯較少，使用 mul128 大數乘法實作的 fibseq_doubling() 竟也大致呈

$O(n)$ 成長、且效能明顯不佳。

• 大數加法中，我們用一般演算法來比較該文作者實作的自動進位加法 add128_auto_carry() 和原先實作的加法 add128() 效能 ：
  
  可以發現，即便 add128() 已經使用 loop-unrolling 的最佳化手法加深 pipeline ，還是顯著地慢上許多。透過 objdump -M intel 可以看到在加法的部分其實也是使用 64 位元的加法實作。

# I'm not pretty sure if this is the actual code segment though. 
# And I only found 3 right shift, but not 4 here
 1c1:   44 8b 20                mov    r12d,DWORD PTR [rax]
 1c4:   8b 48 04                mov    ecx,DWORD PTR [rax+0x4]
 1c7:   8b 50 08                mov    edx,DWORD PTR [rax+0x8]
 1ca:   44 8b 40 0c             mov    r8d,DWORD PTR [rax+0xc]
 1ce:   41 8b 03                mov    eax,DWORD PTR [r11]
 1d1:   45 8b 6b 04             mov    r13d,DWORD PTR [r11+0x4]
 1d5:   41 8b 5b 08             mov    ebx,DWORD PTR [r11+0x8]
 1d9:   45 8b 5b 0c             mov    r11d,DWORD PTR [r11+0xc]
 1dd:   45 8d 34 04             lea    r14d,[r12+rax*1]
 1e1:   4c 01 e0                add    rax,r12
 1e4:   48 c1 e8 20             shr    rax,0x20 #<--------- shift right by 32b
 1e8:   4d 01 d8                add    r8,r11
 1eb:   45 89 34 31             mov    DWORD PTR [r9+rsi*1],r14d
 1ef:   49 89 c4                mov    r12,rax
 1f2:   89 c8                   mov    eax,ecx
 1f4:   4c 01 e8                add    rax,r13
 1f7:   4c 01 e0                add    rax,r12
 1fa:   41 89 44 31 04          mov    DWORD PTR [r9+rsi*1+0x4],eax
 1ff:   48 c1 e8 20             shr    rax,0x20 #<--------- shift right by 32b
 203:   48 89 c1                mov    rcx,rax
 206:   89 d0                   mov    eax,edx
 208:   48 01 d8                add    rax,rbx
 20b:   48 01 c8                add    rax,rcx
 20e:   41 89 44 31 08          mov    DWORD PTR [r9+rsi*1+0x8],eax
 213:   48 c1 e8 20             shr    rax,0x20 #<--------- shift right by 32b

• 大數乘法的部分，我們用 fibseq_doubling()，置換裡面的乘法函式來做測試。乘法迴圈一樣也採用 loop unrolling 的最佳化手法：
  
  可以看得出來，運用硬體乘法單元的費氏數列計算時間成長曲線大致成
  $O(logn)$；而且，寫入時對 cache 較友善的 comba 乘法計算時間幾乎都是傳統演算法的一半。

TODO: 評估 sysprog21/bignum 的表現並探討實作手法

jserv

• sysprog21/bignum 的效能表現

資料傳遞

問題與討論