2023q1 Homework3 (fibdrv)

contributed by < fewletter >

實驗環境

$ uname -r
5.15.0-67-generic

$ gcc --version
gcc (Ubuntu 9.4.0-1ubuntu1~20.04.1) 9.4.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):                          6
On-line CPU(s) list:             0-5
Thread(s) per core:              1
Core(s) per socket:              6
Socket(s):                       1
NUMA node(s):                    1
Vendor ID:                       GenuineIntel
CPU family:                      6
Model:                           158
Model name:                      Intel(R) Core(TM) i5-9500 CPU @ 3.00GHz
Stepping:                        10
CPU MHz:                         3000.000
CPU max MHz:                     4400.0000
CPU min MHz:                     800.0000
BogoMIPS:                        6000.00
Virtualization:                  VT-x
L1d cache:                       192 KiB
L1i cache:                       192 KiB
L2 cache:                        1.5 MiB
L3 cache:                        9 MiB
NUMA node0 CPU(s):               0-5

作業要求

研讀上述 Linux 效能分析的提示描述，在自己的實體電腦運作 GNU/Linux，做好必要的設定和準備工作
$\to$ 從中也該理解為何不希望在虛擬機器中進行實驗;
研讀上述費氏數列相關材料 (包含論文)，摘錄關鍵手法，並思考 clz / ctz 一類的指令對 Fibonacci 數運算的幫助。請列出關鍵程式碼並解說
複習 C 語言數值系統和 bitwise operation，思考 Fibonacci 數快速計算演算法的實作中如何減少乘法運算的成本;
學習指出針對大數運算的加速運算和縮減記憶體操作成本的舉措，紀錄你的認知和疑惑
注意到 fibdrv.c 存在著 DEFINE_MUTEX, mutex_trylock, mutex_init, mutex_unlock, mutex_destroy 等字樣，什麼場景中會需要呢？撰寫多執行緒的 userspace 程式來測試，觀察 Linux 核心模組若沒用到 mutex，到底會發生什麼問題。嘗試撰寫使用 POSIX Thread 的程式碼來確認。

開發紀錄

Linux 時間測量和效能分析

參考作業說明：Linux 時間測量和效能分析和 yanjiew 所撰寫的時間測試檔案。

為了要達到時間測量的目的，依照作業說明：Linux 時間測量和效能分析中的步驟實作

抑制 address space layout randomization (ASLR)

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

設定 scaling_governor 為 performance。準備以下 shell script，檔名為 performance.sh:

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

在測試時間時，利用 $ sudo sh performance.sh 執行，可藉由

cat /sys/devices/system/cpu/cpu*/cpufreq/scaling_governor

檢查 cpu 是否已經進入 performance 模式，測試結果：

performance
performance
performance
performance
performance
performance

針對 Intel 處理器，關閉 turbo mode

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

時間測量

為了能夠在測量程式在 kernel 中運行的時間，將沒有用到的函式 fib_write 來作為輸出時間的函式。

首先先了解 ktime_t 此 api 如何運作，從 ktime_t 中我們可以知道其中的 ktime_get() 是用來標記時間，所以只要在程式開始和結束時都使用ktime_get() 就可以將程式執行時間擷取出來，程式碼如下：

static ktime_t kt;

在 fib_read() 中計時

static ssize_t fib_read(struct file *file,
                        char *buf,
                        size_t size,
                        loff_t *offset)
{
    ktime_t start = ktime_get();
    long long result = fib_sequence(*offset);
    kt = ktime_sub(ktime_get(), start);
    return (ssize_t) result;
}

從 fib_write() 中回傳時間

static ssize_t fib_write(struct file *file,
                         const char *buf,
                         size_t size,
                         loff_t *offset)
{
    return (ssize_t) kt;
}

參考 yanjiew 所撰寫的時間測試檔案，並將其修改成可以輸出 data.txt 來紀錄時間，測試結果

data.txt
Time on this sequence:44
Time on this sequence:24
Time on this sequence:32
Time on this sequence:35
Time on this sequence:37
Time on this sequence:39
Time on this sequence:42
Time on this sequence:42
...

在 scripts 中寫一個簡單的程式可以將上述資料畫圖

data = np.loadtxt('data.txt')

if __name__ == '__main__':
    fig = plt.figure()
    plt.title("python scripts plot")
    plt.xlabel("nth fabonacci")
    plt.ylabel("time(ns)")
    plt.plot(data,label = 'fast doubling')
    plt.legend()
    fig.savefig('data.png')

針對原始費氏數列的遞迴程式碼進行時間測量就會得到下面的圖

static long long fib_sequence(long long k)
{
    /* FIXME: C99 variable-length array (VLA) is not allowed in Linux kernel. */
    long long f[k + 2];

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

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

    return f[k];
}

既然可以從輸出 kernel 的執行時間，那接著便是測試 user space 的時間，利用 time.h 中的 clock_gettime 進行時間的測量，從其中文件可以看到 clockid_t 和 timespec 兩個引數型態

int clock_gettime(clockid_t clockid, struct timespec *tp)

往下搜尋可以看到 timespec 這個結構體，記錄著以奈秒(ns)為時間的單位

struct timespec {
	__kernel_time_t  tv_sec; /* seconds */
	long            tv_nsec; /* nanoseconds */
};

要利用此結構體只需要將全部的成員換成 ns 的時間單位就好

clock_gettime(CLOCK_MONOTONIC, &t1);
long long time = t1.tv_sec * 1e9 + t1.tv_nsec;

clockid_t 則是選擇 CLOCK_MONOTONIC，主要是因為只要形成 process 開始和結束時候的時間標記，程式碼實作如下:

clock_gettime(CLOCK_MONOTONIC, &t1);
sz = read(fd, buf, 1);
clock_gettime(CLOCK_MONOTONIC, &t2);
long long ut = (long long)(t2.tv_sec * 1e9 + t2.tv_nsec)
                 - (t1.tv_sec * 1e9 + t1.tv_nsec);

利用統計方法消除極端值

參考作業參考:用統計手法去除極端值和 4ce43c4 來實作 python scripts，實作此 scripts 的目的為降低極端值對圖的影響，如果像上面那張圖只憑計算一次的時間來畫圖，很可能會有極端值出現降低圖的可信度

上圖為原始算法所算的測試時間，在 kernel 中的測試時間隨著費氏數列的項數而增加，但是在 user 中的時間卻變動不大，這部分還在研究中。

修改 Makefile

為了將上述命令不需要在每次測試程式時都重打一次，將上述命令寫進 Makefile 中並且命名為 test

test: all
	sudo sh -c "echo 0 > /proc/sys/kernel/randomize_va_space"
	sudo sh performance.sh
	sudo sh -c "echo 1 > /sys/devices/system/cpu/intel_pstate/no_turbo"
	$(MAKE) unload
	$(MAKE) load
	sudo taskset 0x1 ./test_time
	$(MAKE) unload
	python3 scripts/plot.py

研讀費氏數列相關材料

參考作業說明一和 Fast Fibonacci algorithms ，可以知道費式數列有多種計算方法，第一種為公式解

F (n) = \frac{1}{\sqrt{5}} (\frac{1 + \sqrt{5}}{2})^{n} - \frac{1}{\sqrt{5}} (\frac{1 - \sqrt{5}}{2})^{n}

此方法最大的缺點在於

\sqrt{5}

此數字在電腦中計算時無法取得確切數值，所以當計算

(\frac{1 + \sqrt{5}}{2})^{n}

時，誤差會因為 n 次方而增加。
第二種方法為直觀的加法，利用費氏數列原本的定義來計算

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

此種方法的缺點則在會重複計算某些不必要的數值，只為了形成

F (k - 1)

和

F (k - 2)

，並且每次計算時必須要每項都重新計算一次。
第三種則是 Fast Doubling，為了能夠減少計算量觀察其中規律

[\begin{matrix} F (2) \\ F (1) \end{matrix}] = [\begin{matrix} 1 & 1 \\ 1 & 0 \end{matrix}] [\begin{matrix} F (1) \\ F (0) \end{matrix}]

\to [\begin{matrix} F (n + 1) \\ F (n) \end{matrix}] = {[\begin{matrix} 1 & 1 \\ 1 & 0 \end{matrix}]}^{n} [\begin{matrix} F (1) \\ F (0) \end{matrix}]

\to [\begin{matrix} F (2 n + 1) \\ F (2 n) \end{matrix}] = {[\begin{matrix} 1 & 1 \\ 1 & 0 \end{matrix}]}^{2 n} [\begin{matrix} F (1) \\ F (0) \end{matrix}]

\to [\begin{matrix} F (2 n + 1) \\ F (2 n) \end{matrix}] = [\begin{matrix} F (n + 1) & F (n) \\ F (n) & F (n - 1) \end{matrix}] [\begin{matrix} F (n + 1) & F (n) \\ F (n) & F (n - 1) \end{matrix}] [\begin{matrix} 1 \\ 0 \end{matrix}]

\to [\begin{matrix} F (2 n + 1) \\ F (2 n) \end{matrix}] = [\begin{matrix} F (n + 1)^{2} + F (n)^{2} \\ F (n) F (n + 1) + F (n) F (n - 1) \end{matrix}]

\to [\begin{matrix} F (2 n + 1) \\ F (2 n) \end{matrix}] = [\begin{matrix} F (n + 1)^{2} + F (n)^{2} \\ F (n) [2 F (n + 1) - F (n)] \end{matrix}]

從以上算式可以注意到費氏數列的規則，

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

和

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

，接著仔細觀察所有的

2 n + 1

和

2 n

的費氏數列都可以由

n

，

n + 1

兩種數值所組成，那我們每次只需要迭代目標一半的數值，一直迭代到目標就是我們的答案。

那目前為止從上面的資訊可以推斷程式碼可以是使用迭代或是遞迴，那規則呢 ? 我們何時需要使用

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

或

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

，這可以從目標的二進位值來進行分析，假設我們想求

F (8)

，那

8

的二進位如下:

8_{10} = 1000_{2}

要如何從

0000

變成

1000

，首先要

+ 1

然後左移三位

(0000 + 1) << 3 = 1000

對比到費氏數列中，就可以理解到為什麼在作業說明中，一開始

F (0) F (1)

變成

F (1) F (2)

一定是

+ 1

，因為二進位的最高位要是

1

，才會有大於

0

的值。

	$+ 1$	$1 \times 2$	$2 \times 2$	$4 \times 2$
$F (0)$	$F (1)$	$F (2)$	$F (4)$	$F (8)$
$F (1)$	$F (2)$	$F (3)$	$F (5)$	$F (9)$

以上是項數最低位為

0

的情況，下面舉例項數最低位為

1

的情況，以

F (13)

為例

13_{10} = 1101_{2}

	$+ 1$	$1 \times 2 + 1$	$3 \times 2$	$6 \times 2 + 1$
$F (0)$	$F (1)$	$F (3)$	$F (6)$	$F (13)$
$F (1)$	$F (2)$	$F (4)$	$F (7)$	$F (14)$

從以上情況可以歸納出兩個規則

遇到最低位元為
$0$ 的時候，可以直接利用上面所推導的算式直接算出
$F (2 n)$ 和
$F (2 n + 1)$

$F (2 n) = F (n) [2 F (n + 1) - F (n)] F (2 n + 1) = F (n + 1)^{2} + F (n)^{2}$
遇到最低位元為
$1$ 的時候，除了要計算
$F (2 n)$ 和
$F (2 n + 1)$ 外，還要利用費氏數列的規則去計算
$F (2 n + 2)$

$F (2 n) = F (n) [2 F (n + 1) - F (n)] F (2 n + 1) = F (n + 1)^{2} + F (n)^{2} F (2 n + 2) = F (2 n + 1) + F (2 n)$

改進費氏數列

Fast doubling 程式實作:
以下程式碼將上面的想法用最直接的方式實作一次

for (uint64_t i = 1LL << 63; i; i >>= 1) {
    long long n0 = f[0] * (2 * f[1] - f[0]);
    long long n1 = f[0] * f[0] + f[1] * f[1];
    if (k & i) {
        f[0] = n1;
        f[1] = n0 + n1;
    }
    else {
        f[0] = n0;
        f[1] = n1;
    }
}

從上圖可以看到光是只有將算法從一般的 iterative 換成 fast doubling 運行時間就降低了不少，但是秉持著任何程式碼都有可以進步的空間，所以繼續優化程式碼，上述程式碼看似沒什麼問題，但是我很想把分支拿掉。

去除分支

為了把分支拿掉，思考良久，寫出適合此情況的 bitmask

for (uint64_t i = 1LL << 63; i; i >>= 1) {
    long long n0 = f[0] * (f[1] * 2 - f[0]);
    long long n1 = f[0] * f[0] + f[1] * f[1];

    uint64_t mask = !(k & i)-1;
    f[0] = (n0 & ~mask) + (n1 & mask);
    f[1] = (n0 & mask) + n1;
}

首先先考慮 f[0] 和 f[1] 的組合，可以從有分支的情況下看到

if (k & i) {
    f[0] = n1;
    f[1] = n0 + n1;
}
else {
    f[0] = n0;
    f[1] = n1;
}

代表 f[0] 和 f[1] 是由 n0 和 n1 所組成，而 f[0] 當中 n0 和 n1 是交替出現，在 f[1] 中，n0 是只有 k & i 不等於 0 時出現，所以可以很清楚的看出其中的 bitmask 的關係

f[0] = (n0 & ~mask) + (n1 & mask);
f[1] = (n0 & mask) + n1;

而什麼樣的 bitmask 會依照 &mask 和 &~mask 去保留和棄置此項，答案就是 1111 和 0000，既然目標是利用 bitwise operation 將 k & i 形成 1111 和 0000，所以就開始湊可以形成這樣的算式
假設

k = 7_{10} = 0111_{2}

i = 1000
  k & i = 0000
  !(k & i) - 1 = 0000
  f[0] = n0; f[1] = n1;

i >= 1 = 1100
  k & i = 0100
  !(k & i) - 1 = 1111
  f[0] = n1; f[1] = n0 + n1;

這樣湊出所需要的 bitmask

從上圖看來兩個方法速度上並沒有明顯差別，可能是測試的項數不夠大又或者是有分支的方法速度本身就和 bitwise operation 差不多快，這些都還需要進一步實驗才可得知。

利用 `clz` 改寫

從 gcc 中的 __builtin_ 提供 __builtin_clz() 可以用來計算出現 1 的最高位元在左邊數過來第幾位，利用此函式修改的程式碼如下

for (uint32_t i = 1UL << (31 - __builtin_clz(k)); i; i >>= 1) {
    long long n0 = f[0] * ((f[1] << 1) - f[0]);
    long long n1 = f[0] * f[0] + f[1] * f[1];

    uint64_t mask = !(k & i)-1;
    f[0] = (n0 & ~mask) + (n1 & mask);
    f[1] = (n0 & mask) + n1;

}

測試的結果顯示利用 clz 的程式碼比較慢，並且有隨著項數增加而時間越長的趨勢，所以去找此函式的程式碼了解為何會出現此情況

int clz (int x){
    if (x == 0) return 32;
    int n = 0;
    if ((x & 0xFFFF0000) == 0) {n += 16; x <<= 16;}
    if ((x & 0xFF000000) == 0) {n +=  8; x <<=  8;}
    if ((x & 0xF0000000) == 0) {n +=  4; x <<=  4;}
    if ((x & 0xC0000000) == 0) {n +=  2; x <<=  2;}
    if ((x & 0x80000000) == 0) n +=  1;
    return n;
}

從程式碼中可以看到 __builtin_clz() 在函式內已經使用了 bitwise operation 來計算數值的 leading zeros，所以在函式內移過一次， for 迴圈開始時又移了一次，時間上會比只有在 for 迴圈移一次還多，至於隨著項數增加而時間越長的趨勢則需要去實作大數運算才能確認。

支援大數運算

創建 `bn` 結構體

為了能夠計算超過 long long 位元組所限制的費氏數列，所以利用 bn 結構體來實作費氏數列。

bn 結構體

typedef struct _bn {
    unsigned int *number;
    unsigned int size;
} bn;

從結構體中可以將大數視為

2^{32}

進位的數字，size 則是此大數的大小，單位則為 32 位元，以費氏數列

F (92) = 7540113804746346429

做舉例，此數字超過 32 位元必須要用 long long 才能儲存，以 bn 結構體表示為

number[0] = 7540113804746346429 & 4294967295
number[1] = (7540113804746346429 >> 32) & 4294967295
size = 2

而問題來了，當有一個數值超過了 long long 所能表示的範圍的時候，我們要如何去表示，所以在這裡參考了作業說明：預先計算儲存

F (n)

所需的空間，可以利用公式解先算出此數值的位元數

F (n) b i t = \frac{- 1160964 + 694242 * n}{10^{6}} + 1

int fibn_per_32bit(int k)
{
    return k < 2 ? 1 : ((long) k * 694242 - 1160964) / (1000000) + 1;
}

接著是如何創建一個 bn 結構體，以下是參考自作業說明

bn *bn_alloc(size_t size)
{
    bn *new = kmalloc(sizeof(bn), GFP_KERNEL);
    new->number = kmalloc(sizeof(int) * size, GFP_KERNEL);
    memset(new->number, 0, sizeof(int) * size);
    new->size = size;
    return new;
}

有 bn_alloc 就要有 bn_free

int bn_free(bn *src)
{
    if (src == NULL)
        return -1;
    kfree(src->number);
    kfree(src);
    return 0;
}

`bn` 運算

然後觀察費氏數列需要什麼運算符號，從最原始的算法中可以看到需要加法和互換兩種運算式

bn_add
此程式在做

2^{32}

進位加法，其中的運算與

10

進位加法完全相同，先算低位元的加法，用 unsigned long long int 去儲存結果，若有進位的數值自然會在右移 32 位元之後出現。

void bn_add(bn *a, bn *b, bn *c)
{
    unsigned long long int carry = 0;
    for (int i = 0; i < c->size; i++) {
        unsigned int tmp1 = (i < a->size) ? a->number[i] : 0;
        unsigned int tmp2 = (i < b->size) ? b->number[i] : 0;
        carry += (unsigned long long int) tmp1 + tmp2;
        c->number[i] = carry;
        carry >>= 32;
    }
}

bn_swap
互換兩個 bn 的位置，在迭代時常會用到

void bn_swap(bn *a, bn *b)
{
    bn tmp = *a;
    *a = *b;
    *b = tmp;
}

bn_to_string
透過 ASCII code 轉換，將數字和 '0' 的差以 char 陣列儲存，並且從高位逐一檢查是否有進位，若有進位則處理成沒有進位的數字

char *bn_to_string(bn src)
{
    size_t len = (8 * sizeof(int) * src.size) / 3 + 2;
    char *s = kmalloc(len, GFP_KERNEL);
    char *p = s;

    memset(s, '0', len - 1);
    s[len - 1] = '\0';

    for (int i = src.size - 1; i >= 0; i--) {
        for (unsigned int d = 1U << 31; d; d >>= 1) {
            /* binary -> decimal string */
            int carry = !!(d & src.number[i]);
            for (int j = len - 2; j >= 0; j--) {
                s[j] += s[j] - '0' + carry;  // double it
                carry = (s[j] > '9');
                if (carry)
                    s[j] -= 10;
            }
        }
    }
    // skip leading zero
    while (p[0] == '0' && p[1] != '\0') {
        p++;
    }
    memmove(s, p, strlen(p) + 1);
    return s;
}

bn_multSSA
利用 Schönhage–Strassen Algorithm 來實作乘法，Schönhage–Strassen Algorithm 的主要觀念是講兩個要相乘的數字分成若干的小數字相乘後相加，而在大數中，我們已經有分好的小數字也就是 bn 結構體中的 number，所以應用在大數中的 Schönhage–Strassen Algorithm 會變成以下型式

[]為此大數的第幾個number
                                  a[2]                 a[1]        a[0]
x                                                      b[1]        b[0]
------------------------------------------------------------------------
                             a[2]*b[0]            a[1]*b[0]   a[0]*b[0] 
        a[2]*b[1]            a[1]*b[1]            a[0]*b[1]
------------------------------------------------------------------------
        a[2]*b[1]  a[2]*b[0]+a[1]*b[1]  a[1]*b[0]+a[0]*b[1]   a[0]*b[0]

在大數

a

和大數

b

相乘之中，

(a_{2} b_{1}, a_{2} b_{0} + a_{1} b_{1}, a_{1} b_{0} + a_{0} b_{1}, a_{0} b_{0})

為

a

和

b

的線性卷積，而若大數

c = a b

，那

c

當中的

n u m b e r

就可以對應到

a

和

b

的線性卷積

a[2]*b[1]  a[2]*b[0]+a[1]*b[1]  a[1]*b[0]+a[0]*b[1]   a[0]*b[0]
    |               |                    |                |
    |               |                    |                |    
   c[3]            c[2]                 c[1]             c[0]

從上面

c

與

a b

的對應關係可以看出

c [x] = \sum_{x = 0, y = 0}^{(a \to s i z e) - 1} a [y] b [x - y]

，其對應的程式碼則如下

unsigned long long int carry = 0;
for (int i = 0; i < c->size; i++) {
    int j = a->size;
    while (j >= 0) {
        unsigned int tmp1 = ((j <= i) && j < a->size) ? a->number[j] : 0;
        unsigned int tmp2 = ((i - j) < b->size) ? b->number[i - j] : 0;
        carry += (unsigned long long int) tmp1 * tmp2;
        j--;
    }
    c->number[i] = carry;
    carry >>= 32;
}

而為了不讓大數

c

的大小在 fast doubling 時一直膨脹，所以在乘法開始執行之前，使用下面程式碼來限制

c

的大小

int d = bn_msb(a) + bn_msb(b);
d = DIV_ROUNDUP(d, 32) + !d;  // round up, min size = 1
bn_resize(c, d);

在測試時發現

F (293)

和 The first 500 Fibonacci numbers 當中

F (293)

數字並不符合

293 7654090467756936378415884538348976340768064993978954512095813

F(293) = 7654090467756936378415203973615134463841138244764090975672901

檢查過後發現在進行加法時，

c [2]

的總和超過了 unsigned long long int 的上限，所以將 unsigned long long int 改成 __uint128_t ，由此可計算超過 1000 項費氏數列

-unsigned long long int carry = 0;
+__uint128_t carry = 0;
for (int i = 0; i < c->size; i++) {
    int j = a->size;
    while (j >= 0) {
        unsigned int tmp1 = ((j <= i) && j < a->size) ? a->number[j] : 0;
        unsigned int tmp2 = ((i - j) < b->size) ? b->number[i - j] : 0;
-       carry += (unsigned long long int) tmp1 * tmp2;
+       carry += (__uint128_t) tmp1 * tmp2;
        j--;
    }
    c->number[i] = carry;
    carry >>= 32;
}

`bn` 支援費氏數列運算

使用 bn 重新實作一次原始費氏數列運算

void bn_fib(bn *ret, unsigned int n)
{
    int nsize = fibn_per_32bit(n) / 32 + 1;
    bn_resize(ret, nsize);
    if (n < 2) {
        ret->number[0] = n;
        return;
    }

    bn *f1 = bn_alloc(nsize);
    bn *tmp = bn_alloc(nsize);
    ret->number[0] = 0;
    f1->number[0] = 1;

    for (unsigned int i = 1; i < n + 1; i++) {
        bn_add(ret, f1, tmp);
        bn_swap(f1, ret);
        bn_swap(f1, tmp);
    }

    bn_free(f1);
    bn_free(tmp);
}

修改 fib_read，使得 fib_read 可以回傳字串到 user space 中

bn *fib = bn_alloc(1);
bn_fib(fib, *offset);
char *p = bn_to_string(*fib);
copy_to_user(buf, p, len);
bn_free(fib);

調整測試上限，便可獲得超過原本

F (92)

的數值，測試結果如下

Reading from /dev/fibonacci at offset 493, returned the sequence 4801994309516818408452995190383973646671835537213025106017041525333156522800379742496995285687643305513.
Reading from /dev/fibonacci at offset 494, returned the sequence 7769790006581794836315417026718114982822736329999469724305586980435409187727711750121668968517028834617.
Reading from /dev/fibonacci at offset 495, returned the sequence 12571784316098613244768412217102088629494571867212494830322628505768565710528091492618664254204672140130.
Reading from /dev/fibonacci at offset 496, returned the sequence 20341574322680408081083829243820203612317308197211964554628215486203974898255803242740333222721700974747.
Reading from /dev/fibonacci at offset 497, returned the sequence 32913358638779021325852241460922292241811880064424459384950843991972540608783894735358997476926373114877.
Reading from /dev/fibonacci at offset 498, returned the sequence 53254932961459429406936070704742495854129188261636423939579059478176515507039697978099330699648074089624.
Reading from /dev/fibonacci at offset 499, returned the sequence 86168291600238450732788312165664788095941068326060883324529903470149056115823592713458328176574447204501.
Reading from /dev/fibonacci at offset 500, returned the sequence 139423224561697880139724382870407283950070256587697307264108962948325571622863290691557658876222521294125

甚至可以測到

F (1000)

的費氏數列

Reading from /dev/fibonacci at offset 993, returned the sequence 1497068822810020534399232227533005158498637007856137201747455252741318759630046838167766787649593980126381385456182553867019240362815567640483762190938317523078461741297449202133947902177195835489685370439138.
Reading from /dev/fibonacci at offset 994, returned the sequence 2422308238804407089268820758059737271890428034963283291423237466738923487762205937312441964529423332944990116110066323372235058978516564015277004931960761593992730308301490464065917906637454573375206452747367.
Reading from /dev/fibonacci at offset 995, returned the sequence 3919377061614427623668052985592742430389065042819420493170692719480242247392252775480208752179017313071371501566248877239254299341332131655760767122899079117071192049598939666199865808814650408864891823186505.
Reading from /dev/fibonacci at offset 996, returned the sequence 6341685300418834712936873743652479702279493077782703784593930186219165735154458712792650716708440646016361617676315200611489358319848695671037772054859840711063922357900430130265783715452104982240098275933872.
Reading from /dev/fibonacci at offset 997, returned the sequence 10261062362033262336604926729245222132668558120602124277764622905699407982546711488272859468887457959087733119242564077850743657661180827326798539177758919828135114407499369796465649524266755391104990099120377.
Reading from /dev/fibonacci at offset 998, returned the sequence 16602747662452097049541800472897701834948051198384828062358553091918573717701170201065510185595898605104094736918879278462233015981029522997836311232618760539199036765399799926731433239718860373345088375054249.
Reading from /dev/fibonacci at offset 999, returned the sequence 26863810024485359386146727202142923967616609318986952340123175997617981700247881689338369654483356564191827856161443356312976673642210350324634850410377680367334151172899169723197082763985615764450078474174626.
Reading from /dev/fibonacci at offset 1000, returned the sequence 43466557686937456435688527675040625802564660517371780402481729089536555417949051890403879840079255169295922593080322634775209689623239873322471161642996440906533187938298969649928516003704476137795166849228875.

將上面的算法跑 50 次後算平均可以發現不管是在 user space 還是 kernal to user 所花的時間都大到誇張，表示在程式在處理 copy_to_user 時花了大量的時間，而上圖中 kernel 看似沒有太大的變化，但當我們單獨拉出 kernel 運行時間的圖

就可以發現，kernel 運行時間也是隨著費氏數列的增加而變多

實作大數 fast doubling

void bn_fastdoub_fib(bn *ret, unsigned int n)
{
    int nsize = fibn_per_32bit(n) / 32 + 1;
    bn_resize(ret, nsize);
    if (n < 2) {
        ret->number[0] = n;
        return;
    }

    bn_reset(ret);

    bn *f1 = bn_alloc(nsize);
    f1->number[0] = 1;
    bn *n0 = bn_alloc(nsize);
    bn *n1 = bn_alloc(nsize);
    
    for (uint32_t i = 1UL << (31 - __builtin_clz(n)); i; i >>= 1) {
        bn *tmp1 = bn_alloc(1);
        bn *tmp2 = bn_alloc(1);
        bn *tmp3 = bn_alloc(1);
        /*n0 = f[0] * (f[1] * 2 - f[0]);*/
        bn_add(f1, f1, n0);
        bn_sub(n0, ret, n0);
        bn_multSSA(n0, ret, tmp1);
        /*n1 = f[0] * f[0] + f[1] * f[1];*/
        bn_multSSA(ret, ret, tmp2);
        bn_multSSA(f1, f1, tmp3);
        bn_add(tmp2, tmp3, n1);

        if (n & i) {
            bn_cpy(ret, n1);
            bn_cpy(f1, tmp1);
            bn_add(f1, n1, f1);
        } else {
            bn_cpy(ret, tmp1);
            bn_cpy(f1, n1);
        }
    }
    bn_free(n0);
    bn_free(n1);
    bn_free(f1);
}

從上面兩張圖可以看到相比於原始費氏數列運算速度快上許多，但是在 kernel space 的運行速度和 user space 的運行速度相比還是相差太多，所以開啟 perf 檢查到底是什麼造成 user space 的時間增加

$ sudo perf record -g --call-graph dwarf ./test_time
$ sudo perf report --stdio -g graph,0.5,caller
96.86%     0.00%  test_time  libc-2.31.so       [.] __GI___read (inlined)
vfs_read
|          
 --96.49%--fib_read
           |          
           |--92.59%--bn_to_string
           |          
            --3.46%--bn_fastdoub_fib
                      |          
                      |--1.51%--bn_multSSA
                      |          |          
                      |           --0.51%--bn_resize
                      |          
                       --0.88%--bn_alloc

2023q1 Homework3 (fibdrv)

實驗環境

開發紀錄

Linux 時間測量和效能分析

時間測量

利用統計方法消除極端值

修改 Makefile

研讀費氏數列相關材料

改進費氏數列

去除分支

利用 clz 改寫

支援大數運算

創建 bn 結構體

bn 運算

bn 支援費氏數列運算

Read more

STM32

Linux核心實作筆記

Implement Linux userspace RV32 emulation for RVVM

Assignment3: single-cycle RISC-V CPU

利用 `clz` 改寫

創建 `bn` 結構體

`bn` 運算

`bn` 支援費氏數列運算