# 2023q1 Homework3 (fibdrv)
contributed by < [cin-cout](https://github.com/cin-cout) >
Github: [cin-cout/fibdrv](https://github.com/cin-cout/fibdrv)

## 開發環境

$ 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):                          12
On-line CPU(s) list:             0-11
Thread(s) per core:              2
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) i7-9750H CPU @ 2.60GHz
Stepping:                        10
CPU MHz:                         1201.821
CPU max MHz:                     4500.0000
CPU min MHz:                     800.0000
BogoMIPS:                        5199.98
Virtualization:                  VT-x
L1d cache:                       192 KiB
L1i cache:                       192 KiB
L2 cache:                        1.5 MiB
L3 cache:                        12 MiB
NUMA node0 CPU(s):               0-11
Vulnerability Itlb multihit:     KVM: Mitigation: Split huge pages
Vulnerability L1tf:              Mitigation; PTE Inversion; VMX conditional cach
                                 e flushes, SMT vulnerable
Vulnerability Mds:               Mitigation; Clear CPU buffers; SMT vulnerable
Vulnerability Meltdown:          Mitigation; PTI
Vulnerability Mmio stale data:   Mitigation; Clear CPU buffers; SMT vulnerable
Vulnerability Retbleed:          Mitigation; IBRS
Vulnerability Spec store bypass: Mitigation; Speculative Store Bypass disabled v
                                 ia prctl and seccomp
Vulnerability Spectre v1:        Mitigation; usercopy/swapgs barriers and __user
                                  pointer sanitization
Vulnerability Spectre v2:        Mitigation; IBRS, IBPB conditional, RSB filling
                                 , PBRSB-eIBRS Not affected
Vulnerability Srbds:             Mitigation; Microcode
Vulnerability Tsx async abort:   Not affected
Flags:                           fpu vme de pse tsc msr pae mce cx8 apic sep mtr
                                 r pge mca cmov pat pse36 clflush dts acpi mmx f
                                 xsr sse sse2 ss ht tm pbe syscall nx pdpe1gb rd
                                 tscp lm constant_tsc art arch_perfmon pebs bts 
                                 rep_good nopl xtopology nonstop_tsc cpuid aperf
                                 mperf pni pclmulqdq dtes64 monitor ds_cpl vmx e
                                 st tm2 ssse3 sdbg fma cx16 xtpr pdcm pcid sse4_
                                 1 sse4_2 x2apic movbe popcnt tsc_deadline_timer
                                  aes xsave avx f16c rdrand lahf_lm abm 3dnowpre
                                 fetch cpuid_fault epb invpcid_single pti ssbd i
                                 brs ibpb stibp tpr_shadow vnmi flexpriority ept
                                  vpid ept_ad fsgsbase tsc_adjust bmi1 avx2 smep
                                  bmi2 erms invpcid mpx rdseed adx smap clflusho
                                 pt intel_pt xsaveopt xsavec xgetbv1 xsaves dthe
                                 rm ida arat pln pts hwp hwp_notify hwp_act_wind
                                 ow hwp_epp md_clear flush_l1d arch_capabilities

## 開發過程

### 大數運算
#### bn API
大數資料結構：

struct bn {
    int cur_size;
    int max_size;
    unsigned long long *digits;
};

1. bn_init

static struct bn *bn_init(int k)
{
    if (k <= 0) {
        return NULL;
    }

    struct bn *new = kmalloc(sizeof(struct bn), GFP_KERNEL);
    if (!new) {
        return NULL;
    }

    new->digits = kmalloc(sizeof(long long) * k, GFP_KERNEL);
    if (!new->digits) {
        kfree(new);
        return NULL;
    }

    memset(new->digits, 0ULL, sizeof(long long) * k);
    new->max_size = k;
    new->cur_size = 1;
    

2. bn_release

static void bn_release(struct bn* a){
    kfree(a->digits);
    kfree(a);
}

3. bn_swap

參考 [sysprog21/bignum/bignum.c](https://github.com/sysprog21/bignum/blob/master/bignum.c) 中的作法。

static void bn_swap(struct bn *a, struct bn *b)
{
    if (!a || !b) {
        return;
    }

    struct bn tmp = *a;
    *a = *b;
    *b = tmp;
}

4. bn_add

判斷 carry 的方法參考[作業解說](https://hackmd.io/@sysprog/linux2023-fibdrv/%2F%40sysprog%2Flinux2023-fibdrv-b)。

static void bn_add(struct bn* dst, struct bn* a, struct bn* b){
    /*b->cur_size must bigger than a->cur_size*/
    if(a->cur_size > b->cur_size){
        bn_swap(a, b);
    }
    int carry = 0;

    for(int i=0; i<b->cur_size; i++){
        if(i>=a->cur_size){
            dst->digits[i] = b->digits[i] + carry;
            if(b->digits[i] + carry == 0ULL){carry = 1;}
            else{carry = 0;}
        }
        else{
            dst->digits[i] = a->digits[i] + b->digits[i] + carry;
            if(a->digits[i] + carry > ~b->digits[i]){carry = 1;}
            else{carry = 0;}
        }
    }
    dst->cur_size = b->cur_size;

    if(carry){
        dst->digits[b->cur_size] = carry;
        dst->cur_size++;
    }
}

verse 2
原本 for 迴圈內的 else ，因為 digits 跟 a b 有可能會是同個位置，所以不能先給 digits 值再進行下一項的 carry 的判斷。
如果用 tmp 暫存 a 或 digits 的值，需要 64 bits的空間。
最後決定改為佔存 last_carry 需要空間較小。

static void bn_add(struct bn *dst, struct bn *a, struct bn *b)
{
    if (!dst || !a || !b) {
        return;
    }
    /*b->cur_size must bigger than a->cur_size*/
    if (a->cur_size > b->cur_size) {
        bn_swap(a, b);
    }
    int carry = 0;

    for (int i = 0; i < b->cur_size; i++) {
        if (i >= a->cur_size) {
            dst->digits[i] = b->digits[i] + carry;
            if (b->digits[i] + carry == 0ULL) {
                carry = 1;
            } else {
                carry = 0;
            }
        } else {
            int last_carry = carry;
            if (a->digits[i] + carry > ~b->digits[i]) {
                carry = 1;
            } else {
                carry = 0;
            }
            dst->digits[i] = a->digits[i] + b->digits[i] + last_carry;
        }
    }
    dst->cur_size = b->cur_size;
    if (carry) {
        dst->digits[b->cur_size] = 1ULL;
        dst->cur_size++;
    }
}

5. bn_print

參考 [ray90514](https://github.com/ray90514/fibdrv/blob/master/client.c)。

void print_bn(unsigned long long *digits, long long size)
{
    if (!digits || size <= 0) {
        return;
    }
    unsigned __int128 div = 0;
    for (int i = 0; i < size; i++) {
        div = 0;
        digits[size] = 0;
        for (int j = size - 1; j >= i; j--) {
            div = div << 64 | digits[j];
            digits[j + 1] = div / MAX_DIV;
            div = div % MAX_DIV;
        }
        digits[i] = div;
        size += digits[size] > 0;
    }

    printf("%llu", digits[size - 1]);
    for (int i = size - 2; i >= 0; i--) {
        printf("%019llu", digits[i]);
    }
    printf(".\n");
}

接下來此段為 fast doubling 會使用的 API

6. bn_shift

k 為要 shift 多少次，意即乘多少次2。但 fast doubling 只須 *2 所以通常 k 都等於 1。 OF_CHECK 為 2^63 目的為判斷 long long 最高位是否為 1 ，若為 1 下一項就須額外加 1。

static void bn_shift(struct bn *dst, int k, struct bn *a)
{
    int carry;
    for (int i = 0; i < k; i++) {
        carry = 0;
        for (int j = 0; j < a->cur_size; j++) {
            dst->digits[j] = a->digits[j] << 1;
            if (carry) {
                dst->digits[j]++;
            }
            if (a->digits[j] & OF_CHECK) {
                carry = 1;
            }
            else{
                carry = 0;
            }
        }

        dst->cur_size = a->cur_size;
        if (carry) {
            dst->digits[a->cur_size] = 1ULL;
            dst->cur_size++;
        }
    }
}

7. bn_sub

原本有做一版使用先把 b 轉為二補數，在使用原本 bn_add 直接做相加，但後來考慮到二補數相加必須建立在加數與被加數相同 bit 數的情況，不適用於大數的資料結構，故改使用基本的直式相減。

static void bn_sub(struct bn *dst, struct bn *a, struct bn *b)
{
    if (a->cur_size < b->cur_size) {
        bn_swap(a, b);
    }

    int borrow = 0;
    unsigned long long mll = -1;
    for (int i = 0; i < a->cur_size; i++) {
        if (i >= b->cur_size) {
            dst->digits[i] = a->digits[i];
            if (borrow) {
                if (a->digits[i] == 0) {
                    borrow = 1;
                } else {
                    borrow = 0;
                }
                dst->digits[i]--;
            }
        } else {
            if (a->digits[i] - borrow < b->digits[i] ||
                (borrow == 1 && a->digits[i] == 0)) {
                dst->digits[i] =
                    (mll - b->digits[i] - borrow) + a->digits[i] + 1;
                borrow = 1;
            } else {
                dst->digits[i] = a->digits[i] - b->digits[i] - borrow;
                borrow = 0;
            }
        }
    }

    for (int i = dst->max_size - 1; i >= 0; i--) {
        if (dst->digits[i]) {
            dst->cur_size = i + 1;
            break;
        }
    }
}

8. bn_mul

亦是使用直式乘法的邏輯，限制為因為將值寫入 dst 後，還需繼續使用 a b ，故不能像其他 API 一樣 dst 可以等於 a 或 b，另外 dst 初始值必須為 0。 因為 64 bits * 64 bits 結果最大可能為 128bits ，所以使用 unsigned __int128 res 紀錄乘法結果。做乘法時一定需要先做資料型態的轉換意即 (unsigned __int128)，不然相乘結果會是 long long 型態，則無法求得 carry 。

static void bn_mul(struct bn *dst, struct bn *a, struct bn *b)
{
    /*a or b cannot equal to dst*/
    unsigned __int128 res;
    unsigned long long carry;
    for (int i = 0; i < b->cur_size; i++) {
        res = 0;
        carry = 0ULL;
        for (int j = 0; j < a->cur_size; j++) {
            res = (unsigned __int128)a->digits[j] * (unsigned __int128)b->digits[i] + carry;
            carry = res >> 64;
            if (dst->digits[j + i] > (long long) ~res) {
                carry++;
            }
            dst->digits[j + i] += (long long)res;
        }
        if (carry) {
            dst->digits[a->cur_size + i] = carry;
        }
    }

    for (int i = dst->max_size - 1; i >= 0; i--) {
        if (dst->digits[i]) {
            dst->cur_size = i + 1;
            break;
        }
    }
}

### 實做 fib
#### 迭代
先依照範例的迭代方法，只是使用不同的資料結構，來完成大數運算。

首先要考慮的事在進行迭代時，要知道怎麼定義 bn 的 size。
參考了兩個方法，分別來自於 [KYG-yaya573142](https://github.com/KYG-yaya573142/fibdrv/blob/master/bn_kernel.c) 和 [上課講義](https://hackmd.io/@sysprog/linux2023-fibdrv/%2F%40sysprog%2Flinux2023-fibdrv-b)。

[KYG-yaya573142](https://github.com/KYG-yaya573142/fibdrv/blob/master/bn_kernel.c)：

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

    bn *a = bn_alloc(1);
    bn *b = bn_alloc(1);
    dest->number[0] = 1;

    for (unsigned int i = 1; i < n; i++) {
        bn_cpy(b, dest);
        bn_add(dest, a, dest);
        bn_swap(a, b);
    }
    bn_free(a);
    bn_free(b);
}

程式在 bn_add() 時會引導到 bn_do_add():

static void bn_do_add(const bn *a, const bn *b, bn *c)
{
    // max digits = max(sizeof(a) + sizeof(b)) + 1
    int d = MAX(bn_msb(a), bn_msb(b)) + 1;
    d = DIV_ROUNDUP(d, 32) + !d;
    bn_resize(c, d);  // round up, min size = 1

    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;
    }

    if (!c->number[c->size - 1] && c->size > 1)
        bn_resize(c, c->size - 1);
}

[KYG-yaya573142](https://github.com/KYG-yaya573142/fibdrv/blob/master/bn_kernel.c) 在決定 bn 的 size 是由不斷動態調整（一直 resize bn 的 size ）來決定，所以需要一直 realloc。

[上課講義](https://hackmd.io/@sysprog/linux2023-fibdrv/%2F%40sysprog%2Flinux2023-fibdrv-b)：

#define BUFFSIZE (8 * sizeof(int) * LENGTH / 3 + 2)
#define LOGPHI (0.20898764025)
#define LOGSQRT5 (0.34948500216)

int main()
{
    int offset = 100; /* TODO: try test something bigger than the limit */
    int fd = open(FIB_DEV, O_RDWR);
    if (fd < 0) {
        perror("Failed to open character device");
        exit(1);
    }
    float digits;
    int size;
    for (int i = 0; i <= offset; i++) {
        // calculate how many digits are needed for fib(i)
        // digits needed for fib(n) = n*LOG(Phi) - (LOG √5)
        digits = i * LOGPHI - LOGSQRT5;
        float buf_size = floor(digits / 9);
        size = (int) buf_size;
        char *buf = malloc(sizeof(char) * (BUFFSIZE * size));
        lseek(fd, i, SEEK_SET);
        long long time1 = read(fd, buf, 0);
        long long time2 = read(fd, buf, 1);
        printf("%d %lld %s\n", i, time1, buf); 
        free(buf);
    }

透過Fibonnaci Number 有一個近似值的公式 0.4472135955∗1.61803398875n 利用這個公式我們可以近似出所需的空間。

如此一來一開始就可以分配好空間，不需要不斷 realloc bn 的 size。

因為 linux kernel 為了加速運算，所以不能進行浮點數運算，在 kernel 也無法使用 math.h 函式庫。[上課講義](https://hackmd.io/@sysprog/linux2023-fibdrv/%2F%40sysprog%2Flinux2023-fibdrv-b) 提供了不用浮點數的版本 (−181400625+n×108475312)/10^10 + 1
粗估為 (-181+n*109）/10^4 + 1

最終迭代版本的 fib 為下：
詳細 API 請見上方。

static struct bn *fib_sequence(long long k)
{
    if (k < 0) {
        return NULL;
    }
    int ll_size = (-181 + k * 109) / 1000 + 1;
    struct bn *a = bn_init(ll_size);
    struct bn *b = bn_init(ll_size);
    b->digits[0] = 1ULL;

    for (int i = 2; i <= k; i++) {
        bn_add(a, a, b);
        bn_swap(a, b);
    }

    if (k == 0) {
        bn_swap(a, b);
    }

    bn_release(a);
    return b;
}

#### Result
成功計算到 f(500)。

-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.
-Reading from /dev/fibonacci at offset 500, returned the sequence 139423224561697880139724382870407283950070256587697307264108962948325571622863290691557658876222521294125.
-Reading from /dev/fibonacci at offset 499, returned the sequence 86168291600238450732788312165664788095941068326060883324529903470149056115823592713458328176574447204501.
-Reading from /dev/fibonacci at offset 498, returned the sequence 53254932961459429406936070704742495854129188261636423939579059478176515507039697978099330699648074089624.

#### Fast Doubling
參考 [Calculating Fibonacci Numbers by Fast Doubling](https://chunminchang.github.io/blog/post/calculating-fibonacci-numbers-by-fast-doubling) 寫的非常完整，也有非常完整的程式碼，但需要花功夫對大數運算方面做修改。

邏輯與文章中最後一段的實作相似，使用迭代的方式。需要注意的是 因為數列 input 值是使用 long long 紀錄。所以 mask 須改為 long long 的型態。此外 bn_mul 的限制為因為將值寫入 dst 後，還需繼續使用 a b ，故不能像其他 API 一樣 dst 可以等於 a 或 b，另外 dst 初始值必須為 0。 所以在放入參數時，須考慮限制。

static struct bn *fib_fast_doubling(long long k)
{
    unsigned int h = 0;
    for (long long i = k; i; h++, i >>= 1)
        ;

    int ll_size = (-181 + k * 109) / 1000 + 1;
    struct bn *a = bn_init(ll_size);
    struct bn *b = bn_init(ll_size);
    struct bn *c = bn_init(ll_size);
    struct bn *tmp = bn_init(ll_size);
    struct bn *d = bn_init(ll_size);
    b->digits[0] = 1ULL;
    
    for (long long mask = 1 << (h - 1); mask; mask >>= 1) {
        // c = a * (2 * b - a);188
        // d = a * a + b * b189
        memset(c->digits, 0ULL, sizeof(long long) * ll_size);
        bn_shift(tmp, 1, b);
        bn_sub(tmp, tmp, a);
        bn_mul(c, a, tmp);

        memset(d->digits, 0ULL, sizeof(long long) * ll_size);
        memset(tmp->digits, 0ULL, sizeof(long long) * ll_size);
        bn_mul(d, a, a);
        bn_mul(tmp, b, b);
        bn_add(d, d, tmp);

        if (mask & k) {
            bn_add(b, c, d);
            bn_swap(a, d);
        } else {
            bn_swap(a, c);
            bn_swap(b, d);
        }
    }
    
    bn_release(b);
    bn_release(d);
    bn_release(tmp);
    bn_release(c);
    return a;

}

#### Result
與迭代法相同，亦可成功計算到 f(500)。

-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.
-Reading from /dev/fibonacci at offset 500, returned the sequence 139423224561697880139724382870407283950070256587697307264108962948325571622863290691557658876222521294125.
-Reading from /dev/fibonacci at offset 499, returned the sequence 86168291600238450732788312165664788095941068326060883324529903470149056115823592713458328176574447204501.
-Reading from /dev/fibonacci at offset 498, returned the sequence 53254932961459429406936070704742495854129188261636423939579059478176515507039697978099330699648074089624.

### 效能分析
#### 前置設定

為了避免會影響的效能測量的因素，需要做前置設定，前置設定參考 [上課講義](https://hackmd.io/@sysprog/linux2023-fibdrv/%2F%40sysprog%2Flinux2023-fibdrv-c) 和 [KYWeng](https://hackmd.io/@KYWeng/rkGdultSU#%E9%87%8F%E6%B8%AC%E6%99%82%E9%96%93%E7%9A%84%E6%96%B9%E5%BC%8F)，若後者有指令無法執行，則以上課講義為主。

因為每次開機都需要重新設定一次，仿照 [KYWeng](https://hackmd.io/@KYWeng/rkGdultSU#%E9%87%8F%E6%B8%AC%E6%99%82%E9%96%93%E7%9A%84%E6%96%B9%E5%BC%8F) ，自己寫了一個 [measure.sh](https://github.com/cin-cout/fibdrv/blob/master/scripts/measure.sh)，作用為：前置設定->執行程式並畫圖->恢復預設值。

#!/bin/bash
CPUID=5
ORIG_ASLR=`cat /proc/sys/kernel/randomize_va_space`
ORIG_GOV=`cat /sys/devices/system/cpu/cpu$CPUID/cpufreq/scaling_governor`
ORIG_TURBO=`cat /sys/devices/system/cpu/intel_pstate/no_turbo`
ORIG_SMT=`cat /sys/devices/system/cpu/smt/control`

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

python3 ./scripts/driver.py

sudo sh -c "echo $ORIG_ASLR >  /proc/sys/kernel/randomize_va_space"
sudo sh -c "echo $ORIG_GOV > /sys/devices/system/cpu/cpu$CPUID/cpufreq/scaling_governor"
sudo sh -c "echo $ORIG_TURBO > /sys/devices/system/cpu/intel_pstate/no_turbo"
sudo sh -c "echo $ORIG_SMT > /sys/devices/system/cpu/smt/control"

要特別注意，[KYWeng](https://hackmd.io/@KYWeng/rkGdultSU#%E9%87%8F%E6%B8%AC%E6%99%82%E9%96%93%E7%9A%84%E6%96%B9%E5%BC%8F) 的例子中 CPUID 是使用 7，但是每台電腦 CPU 數量都不同，我的電腦只有 6 個 CPU 所以使用 5 。

若要確認 CPU 數量，可執行：

$ nproc

必須先跑完前置設定執行後的值才是可以設定的量。這是因為在前面的操作中，關閉了超線程和 CPU Turbo Boost，導致可用 CPU 的數量減少。nproc 命令顯示的是可用 CPU 的邏輯數量，而不是實際的物理 CPU 數量。在某些情況下，這兩個數量是不同的。

執行時因為路徑的關係，必須在 ./scripts 外以下方式執行：

$ ./scripts/measure.sh

#### 測量時間效能
參考[上課講義](https://hackmd.io/@sysprog/linux2023-fibdrv/%2F%40sysprog%2Flinux2023-fibdrv-c)的建議，使用原本沒有使用到的 write function 作為時間測量，用 size 決定使用何種 method。

static ssize_t fib_write(struct file *file,
                         const char *buf,
                         size_t size,
                         loff_t *offset)
{
    ktime_t kt;
    struct bn *a;
    switch (size) {
    case 0:
        kt = ktime_get();
        a = fib_sequence(*offset);
        kt = ktime_sub(ktime_get(), kt);
        break;
    case 1:
        kt = ktime_get();
        a = fib_fast_doubling(*offset);
        kt = ktime_sub(ktime_get(), kt);
        break;
    }
    bn_release(a);
    return (ssize_t) ktime_to_ns(kt);
}

client 方面把得到的時間寫檔。

    /*fprint iterative*/
    for (int i = 0; i <= offset; i++) {
        lseek(fd, i, SEEK_SET);
        sz = write(fd, write_buf, 0);
        fprintf(fptr_i, "%d %llu\n", i, sz);
    }

    /*fprint fast doubling*/
    for (int i = 0; i <= offset; i++) {
        lseek(fd, i, SEEK_SET);
        sz = write(fd, write_buf, 1);
        fprintf(fptr_f, "%d %llu\n", i, sz);
    }

也有實做利用 gnuplot 來畫圖，但後來考慮到要做統計以及刪除異常值，直接用 python script 的形式會比起使用 c 來的快速許多，畫圖就直接使用 python 了。

set title "perf"
set xlabel "nth fibonacci"
set ylabel "time (ns)"
set terminal png font " Times_New_Roman,12 "
set output "time.png"
set key left

plot\
"iterative.txt" using 1:2 with linespoints linewidth 2 title "iterative", \
"fast_doubling.txt" using 1:2 with linespoints linewidth 2 title "fast doubling"

#### 統計方法
使用 python script 的方法是基於 [colinyoyo26](https://github.com/colinyoyo26/fibdrv/blob/master/scripts/driver.py) 的想法上，配合自己儲存資料的型態去做更改。作用為將執行100次的結果去除 95% 區間之外的值再做平均。make for_py 是前置的一些編譯。

import subprocess
import numpy as np
import matplotlib.pyplot as plt
runs = 100

def outlier_filter(data, threshold = 2):
    z = np.abs((data - data.mean()) / data.std())
    return data[z < threshold]

def data_processing(data):
    n_f = data[0].shape[1]
    final = np.zeros(n_f)

    for n in range(n_f):
        final[n] = outlier_filter(data[:,0,n]).mean()
    
    return final


if __name__ == "__main__":

    Yi = []
    Yf = []
    subprocess.run('make for_py', shell=True)
    for i in range(runs):
        subprocess.run('sudo taskset -c 5 ./client > /dev/null', shell=True)
        output_i = np.loadtxt('./scripts/iterative.txt', dtype='int').T
        output_f = np.loadtxt('./scripts/fast_doubling.txt', dtype='int').T
        Yi.append(np.delete(output_i, 0, 0))
        Yf.append(np.delete(output_f, 0, 0))

    X = output_i[0]

    Yi = np.array(Yi)
    Yi = data_processing(Yi)

    Yf = np.array(Yf)
    Yf = data_processing(Yf)

    fig, ax = plt.subplots(1, 1, sharey = True)
    ax.set_title('perf', fontsize = 16)
    ax.set_xlabel(r'$n_{th}$ fibonacci', fontsize = 16)
    ax.set_ylabel('time (ns)', fontsize = 16)

    ax.plot(X, Yi, marker = '+', markersize = 3, label = 'iterative')
    ax.plot(X, Yf, marker = '*', markersize = 3, label = 'fast_doubling')
    ax.legend(loc = 'upper left')

    plt.savefig('time_5000.png')
    subprocess.run('sudo rmmod fibdrv || true >/dev/null', shell=True)

#### Result
[fast doubling]((https://chunminchang.github.io/blog/post/calculating-fibonacci-numbers-by-fast-doubling)) 雖然也要迭代，但是次數遠比每項迭代來的少許多。而基本的迭代法會隨著 n 線性成長。

前500項比較，與預期結果相同。

前5000項比較，與預期結果相同。

### 額外參考資料
[Character devices](https://hackmd.io/@combo-tw/ryRp--nQS)
[fseek() ](http://tw.gitbook.net/c_standard_library/c_function_fseek.html)
[費式數列計算機](https://codepen.io/tales00/pen/jKroMJ)