2022q1 Homework3 (fibdrv)

contributed by < jim12312321 >

作業說明

實驗環境

$ gcc --version
gcc (Ubuntu 11.2.0-7ubuntu2) 11.2.0

$ lscpu
Architecture:                    x86_64
CPU op-mode(s):                  32-bit, 64-bit
Byte Order:                      Little Endian
Address sizes:                   39 bits physical, 48 bits virtual
CPU(s):                          20
On-line CPU(s) list:             0-19
Thread(s) per core:              1
Core(s) per socket:              14
Socket(s):                       1
NUMA node(s):                    1
Vendor ID:                       GenuineIntel
CPU family:                      6
Model:                           154
Model name:                      12th Gen Intel(R) Core(TM) i7-12700H
Stepping:                        3
CPU MHz:                         2700.000
CPU max MHz:                     6000.0000
CPU min MHz:                     400.0000
BogoMIPS:                        5376.00
Virtualization:                  VT-x
L1d cache:                       336 KiB
L1i cache:                       224 KiB
L2 cache:                        8.8 MiB
NUMA node0 CPU(s):               0-19

自我檢查清單

研讀上述 Linux 效能分析的提示描述，在自己的實體電腦運作 GNU/Linux，做好必要的設定和準備工作
$\to$ 從中也該理解為何不希望在虛擬機器中進行實驗;
研讀上述費氏數列相關材料 (包含論文)，摘錄關鍵手法，並思考 clz / ctz 一類的指令對 Fibonacci 數運算的幫助。請列出關鍵程式碼並解說
複習 C 語言數值系統和 bitwise operation，思考 Fibonacci 數快速計算演算法的實作中如何減少乘法運算的成本;
研讀 KYG-yaya573142 的報告，指出針對大數運算，有哪些加速運算和縮減記憶體操作成本的舉措？
lsmod 的輸出結果有一欄名為 Used by，這是 "each module's use count and a list of referring modules"，但如何實作出來呢？模組間的相依性和實際使用次數 (reference counting) 在 Linux 核心如何追蹤呢？

搭配閱讀 The Linux driver implementer’s API guide » Driver Basics
注意到 fibdrv.c 存在著 DEFINE_MUTEX, mutex_trylock, mutex_init, mutex_unlock, mutex_destroy 等字樣，什麼場景中會需要呢？撰寫多執行緒的 userspace 程式來測試，觀察 Linux 核心模組若沒用到 mutex，到底會發生什麼問題。嘗試撰寫使用 POSIX Thread 的程式碼來確認。

前置準備工作

排除干擾效能的因素

孤立特定的 cpu
- 在 /etc/default/grub 中用 sudo 權限編輯，並在 GRUB_CMDLINE_LINUX 後加入 isolcpus=0 ，以空出第 0 個 cpu
```
GRUB_CMDLINE_LINUX="isolcpus=0"
```
- 更新 grub (也要使用 sudo 權限)
```
＄sudo update-grub
```
- 重開機
- 大功告成！
```
pid 1's current affinity list: 1-19
```

抑制 address space layout randomization

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

設定 scaling_governor 為 performance

寫腳本，檔名為 performance.sh

for i in /sys/devices/system/cpu/cpu*/cpufreq/scaling_governor
do
echo performance > ${i}
done

執行腳本

$ sudo sh performance.sh

關閉 turbo mode

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

把上述設定全部整合在 performance.sh 中，便於後續重新執行。

初探 fibdrv

將 K04: fibdrv 中提到可以測量時間的方法加進 fibdrv.c
- 需注意不須使用 DEFINE_KTIME 進行初始化，於這個 PATCH 時，DEFINE_KTIME 已被移除。
修改 Makefile 會造成錯誤的指令刪除
```
@diff -u out scripts/expected.txt
```
- 這道指令雖然能讓我們看到輸出檔案跟 scripts/expected.txt 的差別，但是 scripts/expected.txt 中的資料也只是完全沒修改時會印出的內容罷了，留著這到指令會造成make check 時錯誤(我不知道為什麼，不使用 make check 單獨使用這道指令時又沒問題！)且會讓版面很亂，故刪除。

將 ./client 指定在特定的 CPU 運行並紀錄結果

$ sudo taskset -c 0 ./client >./output/output.txt

在 client.c 中加入對 userspace time 的測量
耗費時間散佈圖
Image Not Showing Possible Reasons
- The image file may be corrupted
- The server hosting the image is unavailable
- The image path is incorrect
- The image format is not supported
Learn More →

摘錄 Fibonacci 數關鍵手法，並思考 clz / ctz 一類的指令對 Fibonacci 數運算的幫助

費式數列 - 基本定義

F_{n} = {\begin{array}{r} 0, n = 0 \\ 1, n = 1 \\ F_{n - 2} + F_{n - 1}, n > 1 \end{array}

因此由費式數列的基本定義可以得知，除了前兩項，費式數列的算法為

F_{n} = F_{n - 2} + F_{n - 1}

。

fast Doubling

更詳細的推導在 K04: fibdrv 中有說明，這邊就不重複，僅列出結論。

F (2 k) = F (k) [2 F (k + 1) - F (k)]

F (2 k + 1) = F (k + 1)^{2} + F (k)^{2}

這表示我們可以用

F (k)

和

F (k + 1)

求出

F (2 k)

和

F (2 k + 1)

，舉例來說：

用 fast doubling 求出 F(10):
此時我們可知
$k = 5$ ，表示 F(10) 可透過 F(5) 和 F(6) 求出，而 F(5) 和 F(6) 又分別可透過 F(2),F(3) 和 F(3),F(4) 求出…，以此類推，最後我們可以得知要求出 F(10) ，可透過 F(0),F(1),F(2),F(3),F(4),F(5),F(6) 求出答案。
比起原先定義中從 F(0) 一路算到 F(8),F(9) 才能求出 F(10) 明顯更快。

所以費式數列的定義可進一步變成：

F_{n} = {\begin{array}{r} 0, n = 0 \\ 1, n = 1 \\ F (k) [2 F (k + 1) - F (k)], n i s e v e n a n d k = \frac{n}{2} \\ F (k + 1)^{2} + F (k)^{2}, n i s o d d a n d k = ⌊ \frac{n}{2} ⌋ \end{array}

考量到硬體加速的
$F (n)$

省略
$F (0)$ ，直接從
$F (1)$ 開始;

clz/ffs: 先去除數字 MSB 起算的開頭 0 位元，因為真正的數字不包含 leading 0s，所以計算它也沒用，因此需要 clz 計算有多少 leading 0s
* 遇到 0
$\to$ 進行 fast doubling，也就是求
$F (2 n)$ 和
$F (2 n + 1)$
* 遇到 1
$\to$ 進行 fast doubling，也就是先求
$F (2 n)$ 和
$F (2 n + 1)$ ，再求
$F (2 n + 2)$

關於 clz / ctz 一類的指令如何在運算費式數列中加速，如同上面 K04: fibdrv 中的說明，可以先去除 leading 0s ，直接從 first set bit 開始算起。

耗費時間圖
- n
  $\to$ normal, f
  $\to$ fast doubling
- k
  $\to$ kernel space, u
  $\to$ user space
  Image Not Showing Possible Reasons
  The image file may be corrupted
  The server hosting the image is unavailable
  The image path is incorrect
  The image format is not supported
  Learn More →

從前 92 個的計算結果看來，有無進行 fast doubling 的差異不大。

normal fib seq 程式碼















static long long fib_sequence(long long k)
{
    /* a before b */
    long long a,b;

    a = 0;
    b = 1;

    for (int i = 1; i <= k; i++) {
        a += b;
        _swap(a,b);
    }

    return a;
}

fast doubling fib seq 程式碼

























static long long fast_fib_seq(long long k)
{
    if(k==0)
        return 0;
    /*a before b*/
    long long a,b,tmp1,tmp2;

    a = 0;
    b = 1;
    int idx = __builtin_clzll(k)+1;
    for(;idx<=MAX_LL;idx++){
        tmp1 = a*(2*b-a);
        tmp2 = b*b + a*a;
        a = tmp1;
        b = tmp2;
        if(_get_bit(k,(MAX_LL-idx))){ // current bit == 1
            tmp1 = a+b;
            a = b;
            b = tmp1;
        }
    }

    return a;
}

減少乘法運算的成本

以下是 fast doubling 的虛擬碼 - from K04: fibdrv

Fast_Fib(n)
    a = 0; b = 1;       // m = 0
    for i = (number of binary digit in n) to 1
        t1 = a*(2*b - a);
        t2 = b^2 + a^2;
        a = t1; b = t2; // m *= 2
        if (current binary digit == 1)
            t1 = a + b; // m++
            a = b; b = t1;
    return a;

嘗試改寫平方運算

可以發現其中大部分的乘法運用都是在做平方運算，所以如果能針對平方運算做優化，則理論上能減少整體運算成本。
























#define MAX_INT 32
#define SET_BIT(x,n) (x | 1 << n)
#define CHECK_BIT(x,n) (x >> n & 1)

/* return a*a
 * main idea is using 多項式平方.
 * (a+b+c)^2 = a^2 + b^2 + c^2 + 2ab + 2ac + 2bc.
 */
int square_opt(int a){
    int ret = 0;
    // idx == int index of first set
    int fs_idx = MAX_INT -__builtin_clz(a)-1;
    int cnt = fs_idx;
    do{
        if(CHECK_BIT(a,cnt)){
            ret += SET_BIT(0,cnt<<cnt);
            for(int cur = cnt+1;cur <= fs_idx;cur++){ // handle 2ab
                if(CHECK_BIT(a,cur))
                    ret += SET_BIT(0,cnt+cur+1); // 2ab
            }
        }
    }while(--cnt >= 0);
    return ret;
}

簡單測試 (使用 Online C Compiler 和 <time.h> 進行測量)

using operator "*" :
5656 * 5656 = 31990336
normal square : 100 ns.
using square_opt:
5656 * 5656 = 31990336
opt square : 472 ns.

目前的結果看來是沒有優化成功的，反而成果越糟

Image Not Showing Possible Reasons

The image file may be corrupted
The server hosting the image is unavailable
The image path is incorrect
The image format is not supported

Learn More →

繼續思考有沒有改進空間。

於是我換了一個思路，利用將十進位換成二進位 (e.g.

5 = 2^{2} + 2^{0}

)，然後再利用乘法分配律相加。















/*
 * a * a 
 * = a * (2^n + ... + 2 ^k)
 * = a << n + ... + a << k
 */
int square_opt_v2(int a){
    int ret = 0;
    int cnt = MAX_INT-__builtin_clz(a)-1;
    do{
        if(CHECK_BIT(a,cnt)){
            ret += a << cnt;
        }
    }while(--cnt >= 0);
    return ret;
}

而經過和上面一樣的測試環境和方法，以下是得到的結果。

using operator "*" :
5656 * 5656 = 31990336
cost time :108 ns
using square_opt_v2 :
5656 * 5656 = 31990336
cost time :167 ns

可以看到成果更接近單純使用 * 進行次方運算的速度了，但整體來說還是略慢。

我覺得 v2 的效能損耗主要在第 11 行中的 if(CHECK_BIT(a,cnt)) ，如果能改成 branchless 應該有機會進一步增進效能，因此基於 v2 的 v3 (branchless 的 v2) 如下：










#define MASK(x,n) (-CHECK_BIT(x,n))

int square_opt_v3(int a){
    int ret = 0;
    int cnt = MAX_INT-__builtin_clz(a)-1;
    do{
        ret += (a << cnt) & MASK(a,cnt);
    }while(--cnt >= 0);
    return ret;
}

利用 CHECK_BIT(x,n) 的輸出只有 0 或 1 的特性，製作 MASK 達成 branchless 。

而以下是當前的測試成果：






using operator "*" :
5656 * 5656 = 31990336
cost time :166 ns
using square_opt_v3 :
5656 * 5656 = 31990336
cost time :153 ns

終於有機會比用 * 還快了

Image Not Showing Possible Reasons

The image file may be corrupted
The server hosting the image is unavailable
The image path is incorrect
The image format is not supported

Learn More →

接著需要經過測試才可以知道是不是穩定的表現更好。

* 與 square_opt_v3 差異圖
其中 diff 為調整前減去調整後，意即值大於 0 表示有改進，反之小於 0 表示比原先表現更差。

計算 123456 ^ 2 1000 次並紀錄每次的 diff

可以從目前的測試結果看到， v3 依舊無法比原先的乘法算的更快，也沒有任何測試中能得到比原先表現更好的數據。

嘗試改寫乘法運算

看著程式碼發呆一下後，突然發現其實我在實作的 square_opt 轉個彎後其實就是在實作二進位的乘法！因此我稍微改一下，附上程式跟測試結果圖。

程式碼










/* return a * b */
long long  mul_opt(long long a,long long b){
    long long ret = 0;
    long long cnt = MAX_LL-__builtin_clzll(a)-1;
    do{
        ret += (b << cnt) & MASK(a,cnt);
    }while(--cnt >= 0);
    return ret;
}

* 與 mul_opt 差異圖

計算 123456 * 654321 1000 次並紀錄每次的 diff

還是無法增進效能，甚至也還無法逼近原先的效能。

擴增費式數列的可運算範圍
- 加法

實作成果

可運算最大值
目前到 500-th 都是正確的

sudo taskset -c 0 ./client >./output/output.txt

...
Reading from /dev/fibonacci at offset 499, returned the sequence 86168291600238450732788312165664788095941068326060883324529903470149056115823592713458328176574447204501. cost time in kernel: 1682468 ns,cost time in userspace: 1682894 ns.
Reading from /dev/fibonacci at offset 500, returned the sequence 139423224561697880139724382870407283950070256587697307264108962948325571622863290691557658876222521294125. cost time in kernel: 1691187 ns,cost time in userspace: 1691595 ns.
...

耗費時間散佈圖（log）

發想與細節

實作程式碼： My fibdrv
在看到作業說明後，覺得大數實作（by AdrianHuang）的作法很有意思因此想自己嘗試做做看利用數字字串相加的功能。

結構體



typedef struct fib_node {
    char data[256];
} fib_node;

就是個單純拿來放數字的結構（現在想想好像用 char 的二維陣列就能實作了）。
當初看到大數實作（by AdrianHuang）中的結構體 xs 很複雜，想簡化成簡單明瞭的形式。本結構體一樣能達到擴增到費氏數列第 500 個的功能，因此我不太清楚當初 xs 為什麼要這樣設計。

字串加法

































#include <linux/string.h>
/*
 * add string a and b mathematically.
 * add the char in a and b backward,and record carry if necessary.
 * reverse the answer in the end.
 */
static void string_add(char *a, char *b, char *out)
{
    short tmp = 0, carry = 0;
    long len_a = strlen(a), len_b = strlen(b);
    long len_max = _max(len_a, len_b);
    for (int i = 0; i < len_max; i++) {
        // cppcheck-suppress shiftTooManyBitsSigned
        tmp = ((len_a - i - 1) >> 63) ? 0 : (a[len_a - i - 1] - '0');
        // cppcheck-suppress shiftTooManyBitsSigned
        tmp += ((len_b - i - 1) >> 63) ? 0 : (b[len_b - i - 1] - '0');
        if (carry)
            tmp += carry;
        if (tmp >= 10) {
            tmp -= 10;
            carry = 1;
        } else {
            carry = 0;
        }
        out[i] = tmp + '0';
    }
    if (carry) {
        out[len_max] = carry + '0';
        len_max++;
    }
    out[len_max] = '\0';
    _reverse(out);
}

發想︰
看到大數實作（by AdrianHuang）的實作中要 reverse 很多次且要確保 a 恆大於 b ，我想要把這些步驟省下來，因此直接從 a 和 b 的尾端加回來，最後再將輸出 reverse 即可，且透過 _max 讓 a 不一定要恆大於 b 也能運算，其他邏輯跟大數實作（by AdrianHuang）相異不大。

實作細節︰

將 a 和 b 從後往前逐字元相加
將得出字串反轉得到答案

擴增費式數列的可運算範圍
- fast doubling

實作成果

可運算最大值
目前到 500-th 都是正確的

sudo taskset -c 0 ./client >./output/output.txt

...
Reading from /dev/fibonacci at offset 499, returned the sequence 86168291600238450732788312165664788095941068326060883324529903470149056115823592713458328176574447204501. cost time in kernel: 619861 ns,cost time in userspace: 620154 ns.
  Reading from /dev/fibonacci at offset 500, returned the sequence 139423224561697880139724382870407283950070256587697307264108962948325571622863290691557658876222521294125. cost time in kernel: 609854 ns,cost time in userspace: 610124 ns.
...

耗費時間散佈圖（log）
與用加法實作費氏數列的差異（log）

可以看到利用 fast doubling 後越是後面的費氏數列運算的越快。

發想與細節

實作程式碼： My fast fibdrv

在 fast doubling 實作中，實作邏輯就按照作業說明中提供的虛擬碼。

十進位轉二進位
在 fast doubling 中需要將 k 轉成 2 進位後從高位逐位判斷，在本實作中用 stack 儲存十進位轉二位的位數。











/*
 * convert decimal to binary.
 * store the results in stack.
 */
static void d2b(long long a)
{
    push(a & 1);
    a >>= 1;
    if (a > 0)
        d2b(a);
}

字串減法
































/*
 * out = a - b
 * Make sure that a is always greater than b.
 */
static void string_sub(char *a, char *b, char *out)
{
    // cppcheck-suppress variableScope
    short tmp = 0, borrow = 0;
    for (int i = 0; i < strlen(a); i++) {
        // cppcheck-suppress shiftTooManyBitsSigned
        tmp = ((strlen(a) - i - 1) >> 63) ? 0 : (a[strlen(a) - i - 1] - '0');
        if (borrow < 0)
            tmp -= 1;
        // cppcheck-suppress shiftTooManyBitsSigned
        tmp -= ((strlen(b) - i - 1) >> 63) ? 0 : (b[strlen(b) - i - 1] - '0');
        if (tmp < 0) {
            tmp += 10;
            borrow = -1;
        } else {
            borrow = 0;
        }
        if ((i + 1 == strlen(a)) && tmp == 0)
            break;
        out[i] = tmp + '0';
    }
    if (tmp != 0) {
        out[strlen(a)] = '\0';
    } else {
        out[strlen(a) - 1] = '\0';
    }
    _reverse(out);
}

實作方式基本跟字串加法相同，與加法不同的是在加法中要紀錄進位 carry ，而在減法中要紀錄借位 borrow 。

字串乘法































/*
 * Make sure out have been initialized with 0 before.
 */
static void string_mul(char *a, char *b, char *out)
{
    short tmp = 0, carry = 0;
    for (int ib = 0; ib < strlen(b); ib++) {
        carry = 0;
        for (int ia = 0; ia < strlen(a); ia++) {
            tmp = b[strlen(b) - ib - 1] - '0';
            tmp *= a[strlen(a) - ia - 1] - '0';
            tmp += out[ia + ib] - '0';
            if (carry)
                tmp += carry;
            if (tmp >= 10) {
                carry = tmp / 10;
                tmp %= 10;
            } else {
                carry = 0;
            }
            out[ia + ib] = tmp + '0';
        }
        out[ib + strlen(a)] = carry + '0';
    }
    /*len(a)-1+len(b)-1+1*/
    tmp = strlen(a) + strlen(b) - 1;
    if (carry)
        tmp++;
    out[tmp] = '\0';
    _reverse(out);
}

整體作法就是直式乘法的思維，須注意 out 要先將裡面的值初始為全為 0 的字串，這樣在抓取先前預算過的值才不會出錯。（tmp += out[ia + ib] - '0';）

大數運算的加速策略與運算邏輯

sysprog21/bignum

在老師實作的大數運算中主要由 apm 和 bn 所組成，其中 apm 負責處理高精度運算而 bn 則是利用 apm 中提供的各項基礎運算，包裝成更高階的 API 以利使用。

Arbitrary-Precision arithmetic (apm)

高精度運算是利用數字陣列來模擬大數運算，舉例來說：

uint8_t 的 MAX 是 255 ，那如果要表達 300 怎麼辦？

令一個 uint8_t *digit 的數字陣列來儲存(同時也可以指定此陣列的 size，本例中假設 size = 3)

已知
$300 = 0 \times 256^{2} + 1 \times 256^{1} + 44 \times 256^{0}$

則 300 可以用 digit[2] = 0,digit[1] = 1,digit[0] = 44 表示

視覺化表示就是

$\begin{array}{rrrr} 300 : & 00000000 & 00000001 & 00101100 \\ d i g i t : & d i g i t [2] & d i g i t [1] & d i g i t [0] \\ v a l u e : & 0 & 1 & 44 \end{array}$
最後要輸出給使用者看的時候要再經過格式轉換，也就是 format.c 在做的事。

加法/減法

在加減法中，運算方法其實與正常數值運算雷同，根據指定的 size 從低位元組運算至高位元組並記錄 carry/borrow 。
以加法為例：

size = 1
240 + 50 = 290

$\begin{array}{rrr} 240 + 50 : \\ 1111 & 0000 \\ + & 0011 & 0010 \\ 290 : & 0010 & 0010 \end{array}$

carry = 1

乘法

演算法採用 Karatsuba algorithm ，要運算

U * V

時，首先會先將 U 和 V 視為

U = U_{1} * 2^{N} + U_{0}

V = V_{1} * 2^{N} + V_{0}

，再進行乘法運算。

U * V = (2^{2 N} + 2^{N}) U_{1} V_{1} + 2^{N} (U_{1} - U_{0}) (V_{0} - V_{1}) + (2^{N} + 1) U_{0} V_{0}

而在 APM 中數值的表達方式，本身就符合

U = U_{1} * 2^{N} + U_{0}

的形式。

舉例來說：
假設 4 個 bits 為一組，則

$\begin{array}{rrr} 28 \\ = & 0001 & 1100 \\ = & 1 * 2^{4} & + 12 \end{array}$

在乘法中首先會設定一個閾值 (KARATSUBA_MUL_THRESHOLD)，若

⌊ \frac{s i z e}{2} ⌋

小於閾值則執行一般直式乘法把

U_{1} V_{1}

U_{0} V_{0}

和

(U_{1} - U_{0}) (V_{0} - V_{1})

分別算出來並放在對應的位置上進行加減法運算，即可得到兩數乘積，否則繼續依據閾值拆分。

舉例來說：
假設 4 個 bits 為一組且 KARATSUBA_MUL_THRESHOLD 為 1
要運算
$7 * 28$ ，則

$s i z e = 2$ ,
$7 = 0 * 2^{4} + 7$ ,
$28 = 1 * 2^{4} + 12$

$\to$
$U_{1} V_{1} = 0$ ,
$U_{0} V_{0} = 84$ ,
$(U_{1} - U_{0}) (V_{0} - V_{1}) = - 77$
(為求簡潔，以下皆用 10 進位取代 2 進位，
$0010 \Rightarrow / 2$ )

$\begin{array}{rrrrr} 7 * 28 \\ / 5 & / 4 & 84 \\ + & / 5 & / 4 & (2^{4} + 1) * 84 \\ - & / 4 & / 13 & - 2^{4} * 77 \\ / 0 & / 12 & / 4 & 196 \end{array}$

$\begin{array}{rr} 196 \\ = & 1100 & 0010 \\ = & 12 * 2^{4} & + 4 \end{array}$

上述例子中 28 和 7 分別被拆成兩項，若是進行正常的大數乘法(逐 digit 相乘相加)，則需要

2^{2}

，考慮兩個大數且都被需拆

N

次，則需要

(N + 1)^{2}

次乘法運算。
而使用 Karatsuba algorithm 則會需要

3 * N

次乘法運算。
因此使用 Karatsuba algorithm 對於大數運算能減少許多乘法運算帶來的成本。

bn

在 bn 中進一步組合 APM 中寫好的高精度運算函式成為面向使用者的 API 。
除此之外，也將在 apm 運算中常用的 digit, size 等封裝成 bn 結構體。






typedef struct {
    apm_digit *digits; /* Digits of number. */
    apm_size size;     /* Length of number. */
    apm_size alloc;    /* Size of allocation. */
    unsigned sign : 1; /* Sign bit. */
} bn, bn_t[1];

KYG-yaya573142 的報告

加速運算

改寫 fast donbling 的運算式

將原先老師提供的 fast doubling (簡化中間過程)，以下簡稱為 ver1：

\begin{aligned} [\begin{array}{c} F (2 n + 1) \\ F (2 n) \end{array}] & = {[\begin{array}{c} 1 & 1 \\ 1 & 0 \end{array}]}^{2 n} [\begin{array}{c} F (1) \\ F (0) \end{array}] \\ = [\begin{array}{c} F (n + 1)^{2} + F (n)^{2} \\ F (n) F (n + 1) + F (n - 1) F (n) \end{array}] \end{aligned}

⇓

\begin{aligned} F (2 k) & = F (k) [2 F (k + 1) - F (k)] \\ F (2 k + 1) & = F (k + 1)^{2} + F (k)^{2} \end{aligned}

使用不同的 Q-Matrix 改寫成 (簡化中間過程)，以下簡稱為 ver2：

\begin{aligned} [\begin{array}{c} F (2 n - 1) \\ F (2 n) \end{array}] & = {[\begin{array}{c} 0 & 1 \\ 1 & 1 \end{array}]}^{2 n} [\begin{array}{c} F (0) \\ F (1) \end{array}] \\ = [\begin{array}{c} F (n)^{2} + F (n - 1)^{2} \\ F (n) F (n) + F (n) F (n - 1) \end{array}] \end{aligned}

⇓

\begin{aligned} F (2 k - 1) & = F (k)^{2} + F (k - 1)^{2} \\ F (2 k) & = F (k) [2 F (k - 1) + F (k)] \end{aligned}

在 KYG-yaya573142 的報告聲稱 ver2 可以少掉一次迴圈及避免使用減法，避免使用減法這點很好理解，ver2 的確沒有減法，又因為在 bn_sub 的實現方法為反轉 sign bit 達成減法效果，因此減少減法，的確會達到優化的目的。
而針對少掉一次迴圈的部分我存有一些疑惑，首先為了在同一個基準上比較，首先應該將 ver1 原先的實現方法也加上 __bulitin_clz ，以下為簡易測試函式。























































unsigned int fib_fdoubling_ver1(unsigned int n)
{
    if (n < 2) { //Fib(0) = 0, Fib(1) = 1
        return n;
    }

    unsigned int f1 = 0; //F(k)
    unsigned int f2 = 1; //F(k+1)
    unsigned int k1;
    unsigned int k2;

    for (unsigned int i = 1U << (31 - __builtin_clz(n)); i; i >>= 1) {
        /* F(2k) = F(k) * [ 2 * F(k+1) – F(k) ] */
        k1 = f1 * (2 * f2 - f1); 
        /* F(2k+1) = F(k)^2 + F(k+1)^2 */
        k2 = f1*f1 + f2*f2;
        if (n & i) {
            f1 = k2;
            f2 = k1;
            f2 += k2;
        } else {
            f1 = k1;
            f2 = k2;
        }
    }
    return f1;
}

unsigned int fib_fdoubling_ver2(unsigned int n)
{
    if (n < 2) { //Fib(0) = 0, Fib(1) = 1
        return n;
    }

    unsigned int f1 = 0; // F(k-1)
    unsigned int f2 = 1; // F(k)
    unsigned int k1;
    unsigned int k2;

    for (unsigned int i = 1U << (30 - __builtin_clz(n)); i; i >>= 1) {
        /* F(2k-1) = F(k)^2 + F(k-1)^2 */
        k1 = f2*f2 + f1*f1;
        /* F(2k) = F(k)[2 * F(k-1) + F(k)] */
        k2 = f2 * (2 * f1 + f2);
        if (n & i) {
            f1 = k2;
            f2 = k1;
            f2 += k2;
        } else {
            f1 = k1;
            f2 = k2;
        }
    }
    return f2;
}

其中關鍵的原因在於 ver1 的

F (2 k)

更新方式需要透過

F (k) * v a l u e

，因此必須從 n 的 last set bit 開始迴圈，才可運算正確的結果，若一樣用 1U << (30 - __builtin_clz(n)) ，則會因為一開始

F (k)

為 0 進而無法正確更迭費式數列的值。

若要使用與 ver2 相同的迴圈起點， 1U << (30 - __builtin_clz(n)) (換句話說就是從 n 的 last set bit 的下一個 bit 開始)，則就必須滿足以下任一條件：

由於跳過 last set bit ，因此必須更改 f1 與 f2 的值
推導其他運算 fast doubling 的算式

在 ver2 中利用改變 Q-Matrix 使得運算

F (2 k)

時不會受到

F (k)

為 0 的影響(條件 (2) 的解法)。
但從條件 (1) 著手一樣可以正確處理費式數列的運算，已知由於

F (k)

為 0 會造成運算錯誤，因此只要將其改成 1 後續即可正確運算 (如同以下的 ver1_1)。











+ unsigned int fib_fdoubling_ver1_1(unsigned int n)
- unsigned int fib_fdoubling_ver1(unsigned int n)
 {
    ...
+    unsigned int f1 = 1; 
-    unsigned int f1 = 0;
    ...
+    for (unsigned int i = 1U << (30 - __builtin_clz(n)); i; i >>= 1) {
-    for (unsigned int i = 1U << (31 - __builtin_clz(n)); i; i >>= 1) {
    ...
 }

另外 bn_fib_fdoubling 中本來就有對 n == 1 or 0 的情況做特殊處理，所以更改 f1 的值也不用擔心 n == 0 時會有錯誤的回傳值。









void bn_fib_fdoubling(bn *dest, unsigned int n)
{
    ...
    if (n < 2) {  // Fib(0) = 0, Fib(1) = 1
        dest->number[0] = n;
        return;
    }
    ...
}

最後，我認為在 KYG-yaya573142 的報告造成時間小幅度減少 3% 的原因，主要是來自引入 ver2 的算法後，必須搭配使用 clz 避免運算結果出錯，而引入 clz 便會在 n 值越小的時候，減少越多次迴圈運算，達成整體效能的優化，而非因為減少一次迴圈及取代減法。

改善 bn_do_add 的實作

透過將原先一次將 c->number[0 ... c->size -1] 全部運算的迴圈改成兩個迴圈，變成一個負責 c->number[0 ... bsize -1] ，若 asize > bsize 則再利用一個迴圈將 c->number[bsize ... asize-1] 算完。
不過我認為改善後的實作仍有改進空間，目前想到的有兩點：

$a s i z e \geq b s i z e$

在改善後的實作中須確保

a s i z e \geq b s i z e

才可以後續運算，因此在 bn_do_add 中使用以下指令確保

a s i z e \geq b s i z e

恆成立。


    if (a->size < b->size)  // a->size >= b->size
        SWAP(a, b);

但有沒有可能基於費式數列的特性，本來就可以確保

a s i z e \geq b s i z e

恆成立？

先來看看 bn_add 會使用在什麼地方，除了 bn_sub 中用來實現減法外，其餘都使用在計算費式數列中，讓我們列出 bn_fib_fdoubling 中跟 bn_add 有關的操作。

bn_add(k1, f2, k1); // k1 = 2 * F(k-1) + F(k)

bn_add(k2, k1, f1); // f1 = k1 + k2 = F(2k-1) now

bn_add(f1, f2, f2); // f2 = F(2k+2)

接著將三個會用到 bn_add 的時機分開討論，先講 2 和 3 ，他們很直覺，在第 2 點中此時的 k1 和 k2 分別為

F (k)^{2}

和

F (k - 1)^{2}

，由費式數列的單調遞增特性可知

k 1 \geq k 2

恆成立，第三點中也可以很輕易的判斷執行 bn_add 時

f 1 \geq f 2

恆成立。
而在第 1 點中， k1 和 f2 分別為

2 * F (k - 1)

和

F (k)

，接著可以經過一些簡單的推導得知，

k 1 \geq f 2

恆成立。

\begin{aligned} 2 * F (k - 1) & \geq & F (k) \\ F (k - 1) + F (k - 1) & \geq & F (k - 1) + F (k - 2) \\ F (k - 1) & \geq & F (k - 2) \end{aligned}

在確保所有使用 bn_do_add 的操作都可以確保

a \geq b

恆成立後，便可以改寫 bn_do_add。










 /* |c| = |a| + |b| 
+ * make sure asize always bigger than or equal to bsize
  */
static void bn_do_add(const bn *a, const bn *b, bn *c)
{
     ...
-    if (a->size < b->size)  // a->size >= b->size
-        SWAP(a, b);
...
}

然後改寫 bn_fib_fdoubling 中的 bn_add (原本的寫法就已經符合

a \geq b

因此僅加上註解)。
















void bn_fib_fdoubling(bn *dest, unsigned int n)
{
    ...
+       // k1 size always bigger than f2 size
        bn_add(k1, f2, k1);    // k1 = 2 * F(k-1) + F(k)
       ...
+       // k2 size always bigger than k1 size
        bn_add(k2, k1, f1);    // f1 = k1 + k2 = F(2k-1) now
        if (n & i) {
            ...
+           // f1 size always bigger than f2 size
            bn_add(f1, f2, f2);  // f2 = F(2k+2)
        }
    }
    ...
}

修改第二個迴圈

第二個迴圈主要的功能為將 a 中剩下的數值複製到 c 中並且加上 carry 。但考慮到進位的次數與 asize - bsize 的差值會隨著費式數列越大而越大，因此可以再拆解第二個迴圈，讓第二個迴圈負責處理 carry 而第三個負責處理剩下的 a->number 複製到 c->number 。
























/* |c| = |a| + |b| */
static void bn_do_add(const bn *a, const bn *b, bn *c)
{
    ...
    bn_data carry;
    carry = _add_partial(a->number, b->number, bsize, c->number);
    if (asize != bsize) {  // deal with the remaining part if asize > bsize
+        int i = bsize;
+        for (; i < asize & !!carry; i++) { // !!carry == 0 if carry == 0, == 1 if others
-        for (int i = bsize; i < asize; i++) {
            bn_data tmp1 = a->number[i];
            carry = (tmp1 += carry) < carry;
            c->number[i] = tmp1;
        }
+        for (; i < asize; i++) {
+            c->number[i] = a->number[i];
+        }
    }

    if (carry) {
        c->number[asize] = carry;
        ++(c->size);
    }
}

縮減記憶體操作成本

降低 malloc 與 memcpy 的使用

原先的實作(節錄片段)












void bn_fib_fdoubling(bn *dest, unsigned int n)
{
    ...
        /* F(2k) = F(k) * [ 2 * F(k+1) – F(k) ] */
        bn_cpy(k1, f2);       // k1 = F(k+1)
        bn_lshift(k1, 1);     // k1 = 2 * F(k+1)
        bn_sub(k1, f1, k1);   // k1 = 2 * F(k+1) - F(k)
        bn_mult(k1, f1, k1);  // k1 = F(k) * [ 2 * F(k+1) - F(k) ]
    ...
        bn_cpy(f1,k1); // f1 = k1
    ...
}

可以發現為了將運算結果都儲存在 k1 中，因此每次運算都會遇到以下問題：

資料複製
乘法運算時 c == a or c == b 動態配置暫存記憶體

而改進的做法便是重構算式，將 k2 也拿來儲存計算結果，從而避免問題 (2) ，再利用 bn_swap 省去 memcpy 的運用解決問題 (1) 。












void bn_fib_fdoubling(bn *dest, unsigned int n)
{
    ...
        /* F(2k) = F(k) * [ 2 * F(k+1) – F(k) ] */
        /* F(2k+1) = F(k)^2 + F(k+1)^2 */
        bn_lshift(f2, 1, k1);// k1 = 2 * F(k+1)
        bn_sub(k1, f1, k1);  // k1 = 2 * F(k+1) – F(k)
        bn_mult(k1, f1, k2); // k2 = k1 * f1 = F(2k)
        bn_mult(f1, f1, k1); // k1 = F(k)^2
        bn_swap(f1, k2);     // f1 <-> k2, f1 = F(2k) now
    ...
}

引入 memory pool 的概念

在 bn 中加入 capacity 儲存實際配置的長度，而 size 變成當前使用的大小，進而減少頻繁呼叫 realloc 所產生的效能損耗。
另外在 KYG-yaya573142 的報告中提到

至少能正確計算至第 100000 項

是因為 unsigned int 能表達的數值上界為 4294967295(

2^{32} - 1

，在 LP64 或 LLP64 系統)，因此能表達的費式數列上界為

2^{32 * (2^{32} - 1)}

，考慮可能不需要計算到這麼後面的費式數列值，還可以引入 bit field 的技巧將 size,capacity 和 sign 整合。

typedef struct _bn {
    unsigned int *number;  /* ptr to number */
-    unsigned int size;     /* length of number */
-   unsigned int capacity; /* total allocated length, size <= capacity */
-    int sign;
+    unsigned int size: SIZE;
+    unsigned int capacity: CAPACITY;
+    unsigned int sign: SIGN;
} bn;

更進一步的節省 bn 的空間，不過這樣改寫的話原先很多函式都要一並改寫，整體利弊尚未可知。

因 word size 制宜

根據當前系統不同的 word size 調整 bn 中各項屬性的資料型態，不僅能使 bn 能夠在不同系統運行，同時也確保當使用 64-bit 的 CPU 時能夠增加計算效能與儲存空間。

TODO:
針對 KYG-yaya573142 的報告減少大數運算的成本中的以下章節撰寫報告

改善 bn_mult 的效能
引入 bn_sqr
實作 Karatsuba algorithm

mutex 在 fibdrv 中的應用場景與移除時會發生的事情

參考了 KYG-yaya573142 的報告撰寫一個多執行緒程式來測試。


























void *worker(void *arg){
    int fd = open("/dev/fibonacci",O_RDWR);
    if (fd < 0){
        perror("Failed to open character device");
        goto end;
    }
    int idx = *((int *)arg);
    char buf[100];
    for (int offset = 80; offset <=92; offset++){
        lseek(fd,offset,SEEK_SET);
        long long result = read(fd,buf,0);
        printf("thread %d F(%d): %lld\n",idx,offset,result);
    }
end:
    close(fd);
}

int main(){
    pthread_t test[2];
    int idx[2] = {1,2};
    pthread_create(&test[0],NULL,worker,&idx[0]);
    pthread_create(&test[1],NULL,worker,&idx[1]);
    for(int i = 0;i < 2; i++)
        pthread_join(test[i],NULL);
    return 0;
}

未移除 mutex

在原先的程式中，利用 mutex_trylock 和 mutex_unlock 來避免程式進入 critical section ，從而保證每個執行緒的內容正確。而其中使用的 mutex_trylock 當發現該區段已經被 lock 時，會直接返回，屬於非阻斷式 I/O 。
而這樣的特點也可以經由實驗發現，當反覆嘗試原本的 fibdrv 後，會出現兩種結果，一種是出現錯誤訊息的，另一種是兩個執行緒都成功的完成任務。

$ sudo ./thread
thread 1 F(80): 23416728348467685
Failed to open character device: Device or resource busy
thread 1 F(81): 37889062373143906
thread 1 F(82): 61305790721611591
thread 1 F(83): 99194853094755497
thread 1 F(84): 160500643816367088

$ sudo ./thread
thread 1 F(80): 23416728348467685
thread 1 F(81): 37889062373143906
thread 1 F(82): 61305790721611591
...
thread 2 F(80): 23416728348467685
thread 2 F(81): 37889062373143906
thread 2 F(82): 61305790721611591
...

第一種情況發生的原因是因為兩個執行緒針對 fibdrv 的讀寫順序是相互交雜的，因此當 fibdrv 已經被 lock 時，另一個執行緒嘗試 trylock 失敗，因此印出錯誤訊並結束該執行緒。
而第二種狀況則是運氣好，兩個執行緒的讀寫沒有交雜，因此兩個執行緒都完成各自的任務。

將 trylock 改成 lock

在改成使用 mutex_lock 後，當執行緒程式存取 mutex 時，若該區段已被 lock ，則會持續等待直到 unlock ，屬於阻斷式 I/O 。

$ sudo ./thread
thread 1 F(80): 23416728348467685
thread 1 F(81): 37889062373143906
thread 1 F(82): 61305790721611591
...
thread 2 F(80): 23416728348467685
thread 2 F(81): 37889062373143906
thread 2 F(82): 61305790721611591
...

去除 mutex

去除 mutex 後， thread 將可以依照原先的順序讀取，不用等待其他 thread 釋出程式驅動。

$ sudo ./thread
thread 2 F(80): 23416728348467685
thread 2 F(81): 37889062373143906
thread 2 F(82): 61305790721611591
thread 2 F(83): 99194853094755497
thread 1 F(80): 23416728348467685
thread 1 F(81): 37889062373143906
thread 1 F(82): 61305790721611591
thread 2 F(84): 160500643816367088
thread 2 F(85): 259695496911122585
thread 1 F(83): 99194853094755497
thread 1 F(84): 160500643816367088

可以看到 thread 1 和 2 交互讀取 fibdrv ，不過有趣的是可以發現其中每個費氏數列的值並沒有出錯，這是因為在原先撰寫的測試程式中， file descriptor 並不是共享資源， thread 間沒有共用同一個 fd 自然也不會有衝突。

在把 fd 改為全域變數之後，可以發現 thread 所得到的值出現錯誤，意味著進入了 critical section 。

$ sudo ./thread
thread 1 F(80): 23416728348467685
thread 1 F(81): 37889062373143906
thread 1 F(82): 61305790721611591
...
thread 1 F(89): 1779979416004714189
thread 2 F(80): 1779979416004714189
thread 2 F(81): 37889062373143906
thread 2 F(82): 61305790721611591

不過猜測由於測試程式僅有兩個執行緒，因此在進行 lseek 和 read 的操作時，要出現如 thread 1 lseek 但還沒 read 時， offset 已經被 tread 2 改變造成結果不正確的 race condition 反而機率不大，如果採用像 KYG-yaya573142 的報告中用 sleep 的方法，就可以很明顯的看到效果。

$ sudo ./thread
thread 1 F(80): 23416728348467685
thread 2 F(80): 37889062373143906
thread 1 F(81): 37889062373143906
thread 1 F(82): 61305790721611591
thread 2 F(81): 99194853094755497
thread 1 F(83): 61305790721611591
thread 1 F(84): 160500643816367088
thread 2 F(82): 259695496911122585

或是單純的增加 thread 也可以較容易觀察到錯誤（增加到 20 個 thread）

sudo ./thread
...
thread 18 F(83): 99194853094755497
thread 6 F(84): 61305790721611591
thread 3 F(84): 160500643816367088
thread 19 F(80): 23416728348467685
...

額外議題

xs 的設計與記憶體開銷

在實作大數加法的時候，我直接利用 char* 來儲存數字，企圖"簡化" xs 的架構來達到同樣可以大數運算的效果，但就在多看幾次 SSO (Small/Short String Optimization) 和聽老師上課後才稍微理解為何 xs 要這樣設計。

xs 的架構






























typedef union {
    /* allow strings up to 15 bytes to stay on the stack
     * use the last byte as a null terminator and to store flags
     * much like fbstring:
     * https://github.com/facebook/folly/blob/master/folly/docs/FBString.md
     */
    char data[16];

    struct {
        uint8_t filler[15],
            /* how many free bytes in this stack allocated string
             * same idea as fbstring
             */
            space_left : 4,
            /* if it is on heap, set to 1 */
            is_ptr : 1, is_large_string : 1, flag2 : 1, flag3 : 1;
    };

    /* heap allocated */
    struct {
        char *ptr;
        /* supports strings up to 2^MAX_STR_LEN_BITS - 1 bytes */
        /* MAX_STR_LEN_BITS is 54 in this case,
         * so xs supports strings up to 2^54 -1 bytes */
        size_t size : MAX_STR_LEN_BITS, 
                      /* capacity is always a power of 2 (unsigned)-1 */
                      capacity : 6;
        /* the last 4 bits are important flags */
    };
} xs;

可以看到 xs 在佔用 16 bytes 且小字串的情況下，可以利用 15 bytes 儲存字串資料，利用 space_left 來計算剩餘可用的 bytes ，且 15 bytes 可用 4 個 bits 就可以儲存資料剩餘大小 (

2^{4} - 1 = 15

) 。最後剩下 4 個 bits 的 is_ptr 和 is_large_string 可以分別對應到 FBString 中提到的中字串與大字串，而 flag2 與 flag3 我認為只是把剩下的 bits 定義完，實際上並不會用到。
至於在處理中字串與大字串時則利用 heap allocated 的 struct 來表達。

xs 記憶體開銷

xs 這樣的設計使得在小字串時可以直接用 stack 中的記憶體，避免動態配置新記憶體，空間利用上也更完善，直接在原先配置的空間中紀錄了 size, capaity 和盡可能的儲存字串資料。

2022q1 Homework3 (fibdrv)

實驗環境

自我檢查清單

前置準備工作

初探 fibdrv

摘錄 Fibonacci 數關鍵手法，並思考 clz / ctz 一類的指令對 Fibonacci 數運算的幫助

費式數列 - 基本定義

fast Doubling

考量到硬體加速的 F(n)

減少乘法運算的成本

嘗試改寫平方運算

嘗試改寫乘法運算

擴增費式數列的可運算範圍 - 加法

實作成果

發想與細節

擴增費式數列的可運算範圍 - fast doubling

實作成果

發想與細節

大數運算的加速策略與運算邏輯

sysprog21/bignum

Arbitrary-Precision arithmetic (apm)

bn

KYG-yaya573142 的報告

加速運算

縮減記憶體操作成本

mutex 在 fibdrv 中的應用場景與移除時會發生的事情

額外議題

xs 的設計與記憶體開銷

考量到硬體加速的
$F (n)$

擴增費式數列的可運算範圍
- 加法

擴增費式數列的可運算範圍
- fast doubling