# Fibonacci Device
contributed by < `xxiex123` >
> [GitHub](https://github.com/zodf0055980/fibdrv)
環境
> linux 5.8.0-59-generic
```
$ 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): 8
On-line CPU(s) list: 0-7
Thread(s) per core: 2
Core(s) per socket: 4
Socket(s): 1
NUMA node(s): 1
Vendor ID: GenuineIntel
CPU family: 6
Model: 94
Model name: Intel(R) Core(TM) i7-6700HQ CPU @ 2.60GHz
Stepping: 3
CPU MHz: 1990.861
CPU max MHz: 3500.0000
CPU min MHz: 800.0000
BogoMIPS: 5199.98
Virtualization: VT-x
L1d cache: 128 KiB
L1i cache: 128 KiB
L2 cache: 1 MiB
L3 cache: 6 MiB
$ lsmem
RANGE SIZE STATE REMOVABLE BLOCK
0x0000000000000000-0x0000000077ffffff 1.9G online yes 0-14
0x0000000100000000-0x0000000387ffffff 10.1G online yes 32-112
Memory block size: 128M
Total online memory: 12G
Total offline memory: 0B
```
## 背景知識補充
#### Linux Kernel modules
> Kernel modules are pieces of code that can be loaded and unloaded into the kernel upon demand. They extend the functionality of the kernel without the need to reboot the system.
Linux 系統之所以會具備 kernel modules 功能主要有以下兩點
* 由於 linux kernel 常常被不同的應用和設備使用，爲了能夠達到兩者之間的相依性，使用者常常想要加一些針對自己應用或設備的程式在 linux kernel。
* 由於 linux 屬於 monolithic kernel ，簡單來說是一個非常大型的程式，因此若要在它裏面加上一些小功能，再去編譯，會耗費非常大的時間和力氣
因此 linux 有一套 kernel module 機制，可以在不對 kernel 進行編譯的情況下，module 能動態的從 linux kernel insert 或 remove。除此之外，kernel module 還具備各種管理上的方便。
所有被掛載到 linux 上的 module 都放在 `/lib/modules/$(uname -r)/kernel/` 裏面。
#### kernel module 的一些特性：
* 是一種 object code
* 不能獨立運行
* 這種 object code 用來實現一種文件系統，一個驅動程序，或其他内核上層的功能。
* module 不是作爲一个 process 來執行的，與其他 static library 的 kernel function 一樣，它在 kernel space 代表當前的 process 執行。
由於 kernel module 實在 kernel space 執行，因此在寫一個 module 時，就要知道 kernel space 與 user space coding 的差別。
* kernel space 不能調用 user space 定義的函數例如 libc，只能使用在 kernel space 定義的那些資源受到限制的函數，例如 printk()
* kernel space 的程式碼可允許不嚴重的錯誤發生，例如 segmentation fault，但過於嚴重的錯誤會對 process 照成致命。
* kernel space 有內核可用內存的限制，例如 stack 空間一般只有 4kbyte（在系統編譯時就決定好），而 user space 則是可以動態改變。
* 由於 kernel 最後會把 module 裏的函數或變量 export 到 symbol table，爲了避免與 kernel 其他變數或內部函數名稱重疊，要確保每個名稱都是 unique ，因此推薦的命名方式是用自己的 module name 作爲名稱的 prefix 。
* 在 kernel space 因爲效能原因，默認 gcc 不允許浮點數運算，另外，必須自己實做除法和 mod 運算。（這點在此專案照成許多不方便）
#### character device driver
Device driver 分成主要有兩大類，是由傳遞資料的方式和資料量來區別。
1. character device driver（主要）
傳遞資料的方式通常是用 byte by byte 的 pipe 或 serial port，也就是一連串的 byte 資料傳送進來。通常傳遞的資料量不多。
2. block device driver（主要）
傳遞資料的方式是一整個 data block 爲一個單位，例如硬盤、ram disk、file 等等，比較大量的資料傳輸。
3. network device driver
處理網絡資訊用的 device driver ，例如從 Ethernet board 到 TCP/IP 的 IRIX 实现。
由 [character device driver](https://linux-kernel-labs.github.io/refs/heads/master/labs/device_drivers.html)
> the character device drivers receive unaltered system calls made by users over device-type files.Consequently, implementation of a character device driver means implementing the system calls specific to files: open, close, read, write, lseek, mmap, etc.
我們對 character device driver 的操作是透過 implement 檔案的系統呼叫如：open，close，write，lseek等等來實現。這種機制是使用 virtual file system 來達到。
#### symbol table
在 linux kernel 裏，每次編譯都會產生一個對應的 system.map (位於/或者/boot、/usr/src/linux/下)，其記錄所有內核的內部函數或變量實際分配的位置（地址值），也就是把內部函數或變量的名稱都對應到其自身的地址值，作用是當內核運行出錯，kernel 可以由地址值找到對應的變量名稱，顯示出來讓開發者容易解讀（因爲人類看不懂一串地址值代表什麼）。
但 system.map 是在 kernel 編譯的時候所產生的，因此我們掛載 module 時，其內部的 symbol 並不會加入到 system.map，而是會加入到 /proc/kallsyms 裏面（也是一種 symbol table）。
兩者的區別在於：
**system.map** ：文件面向内核，对于内核中的没有导出的变量或者函数名，比如 kthread_create_list ，也有其相应的内核地址，该文件一般是 only read、固定大小的，没有动态添加模块中的变量、函数名。
**kallsyms** ：在内核啓動過程中創建，並實時更新，反映的是系統的當前最新情況，其內部也包含内核或者是已加載 module 所 exported 的函數、變量名稱。
所以我們在掛載 module 時，改變的是 /proc/kallsyms 。
## 開始着手
### 前置準備與驗證
#### 追蹤 make check
輸入指令 sudo make check 可得到輸出
```
make -C /lib/modules/5.8.0-55-generic/build M=/home/yangyang/fibdrv modules
make[1]: Entering directory '/usr/src/linux-headers-5.8.0-55-generic'
make[1]: Leaving directory '/usr/src/linux-headers-5.8.0-55-generic'
make unload
make[1]: Entering directory '/home/yangyang/fibdrv'
sudo rmmod fibdrv || true >/dev/null
rmmod: ERROR: Module fibdrv is not currently loaded
make[1]: Leaving directory '/home/yangyang/fibdrv'
make load
make[1]: Entering directory '/home/yangyang/fibdrv'
sudo insmod fibdrv.ko
make[1]: Leaving directory '/home/yangyang/fibdrv'
sudo ./client > out
make unload
make[1]: Entering directory '/home/yangyang/fibdrv'
sudo rmmod fibdrv || true >/dev/null
make[1]: Leaving directory '/home/yangyang/fibdrv'
Passed [-]
f(93) fail
input: 7540113804746346429
expected: 12200160415121876738
```
其中，Makefile 的做法是
1. 先 unload 確保 fibdiv 是卸載的
2. 再 load 掛載 fibdiv
3. 執行 client 並把結果輸出到 out
4. unload 卸載 fibdiv
5. 查看 out 結果是否與 scripts/expected.txt 一致，若成功輸出綠色的 passed [-] （其中 @ 代表不輸出這行指令）
6. 執行 scripts/verify.py
verify.py 是一個用 python 算出 f(x) , where x <= 100 的費氏數列結果，並驗證 out 內的結果是否與正確答案一致。
爲了能夠驗證 x > 100 的結果，我們利用 MAX_FNUM 把 x 調整成可變。
```python
7 MAX_FNUM = 0
8
9 with open('out', 'r') as f:
10 tmp = f.readline()
11 while (tmp):
12 MAX_FNUM = MAX_FNUM + 1
13 result.append(tmp)
14 tmp = f.readline()
15 f.close()
16 for i in range(2, MAX_FNUM):
17 expect.append(expect[i - 1] + expect[i - 2])
```
由於 python 的針對 [long type](https://en.wikibooks.org/wiki/Python_Programming/Numbers) 的處理
> Long: Integer type with unlimited length. In python 2.2 and later, Ints are automatically turned into long ints when they overflow. Dropped since Python 3.0, use int type instead.
>
理論上我們能計算到 MAX_FNUM 非常大，主要是記憶體有多少。這裏算個 f(700)試試
```
f( 698 )= 33410878289444957618607188493530847778173900120569106083403082529157749504523093168915318986912085240482480870354510205414873169790081851662484849
f( 699 )= 54059936666307888585371224524040479564193340847128274990827350063369752406767284486712908163966342091210712498754683466915904358153636317442639426
```
Makefile 第五步，因爲 expected.txt 是一個寫死的檔案，因此若要可以執行更大的費氏數列的 make check ，可以利用 python 寫出對應的 expected.txt 來做驗證，但我覺得因爲 Makefile 第六步已經做到驗證的效果，因此決定把第5步移除。
如此一來，我們便能夠用 make check 驗證我們的 fibdiv。
#### 加入 big num 資料結構
爲了計算 f(93) 或以上的費氏數列，我們需要實做能夠處理大數的資料結構，考慮到有對於處理無限大的 f(x) 有 [GMP](https://gmplib.org/) 的工具幫助，這裏先利用 __int128 來實做出一個能夠處理足夠大的 f(x) 用以進行演算法的效能評估。
把 fibdiv 裏的 long long 改成 __int128
``` c
static __int128 fib_sequence(long long k)
36 {
37 /* FIXME: use clz/ctz and fast algorithms to speed up */
38 __int128 f[k + 2];
39
40 f[0] = 0;
41 f[1] = 1;
42
43 for (int i = 2; i <= k; i++) {
44 f[i] = f[i - 1] + f[i - 2];
45 }
46
47 return f[k];
48 }
49
50 static __int128 fib_time_proxy(long long k)
51 {
52 kt = ktime_get();
53 __int128 result = fib_sequence(k);
54 kt = ktime_sub(ktime_get(),kt);
55
56 return result;
57 }
```
而由[費氏數列](https://hackmd.io/@sysprog/fibonacci-number)的分析，
> 理論分析方面，我們得先知道 F(n) 大約需要多少個 bit 來表示。若 N(n) 是 F(n) 需要的 bit 數，大約可推算出 N(n)=n×0.69424
我們可得知 128 位元可準確計算出約 f(180) 的結果。
另外，注意到這裏的 kt 計算出的時間是不包括 copy_to_user 的時間成本的。
#### 利用 copy_to_user
爲了從 kernel space 取得資料，把 fib_read 計算出的結果用 copy_to_user 把結果復制到 buf （buf 爲__int128 形態，到了 user 端再轉 string）
```c
75 static ssize_t fib_read(struct file *file,
76 char *buf,
77 size_t size,
78 loff_t *offset)
79 {
80 __int128 tmp = fib_time_proxy(*offset);
81 copy_to_user(buf, &tmp, size);
82 return 1;
83 }
```
由於 fib_write 用來輸出 fib 的時間成本，因此直接透過 return 的方式回傳 value 就行。
```c
86 static ssize_t fib_write(struct file *file,
87 const char *buf,
88 size_t size,
89 loff_t *offset)
90 {
91 return (int) (ktime_to_ns(kt));
92 }
```
#### 更改 client.c 內容以方便測試
把 write 寫在 read 之後，以取得上一次 read 的時間消耗，之後輸出到檔案 time.txt 再利用 gnuplot 呈現圖形
```c=49
for (int i = 0; i <= offset; i++) {
__int128 buf;
lseek(fd, i, SEEK_SET);
read(fd, &buf, BUF_SIZE);
time = write(fd, write_buf, strlen(write_buf));
printf("Reading from " FIB_DEV
" at offset %d, returned the sequence ", i);
printBuff(buf);
printf(".\n");
fprintf(f_time, "%d %d\n", i, time);
}
```
由此我們開始進行測試，先執行 make check
``` bash
yangyang@yangyang-GL552VX:~/fibdrv$ sudo make check
make -C /lib/modules/5.8.0-55-generic/build M=/home/yangyang/fibdrv modules
make[1]: Entering directory '/usr/src/linux-headers-5.8.0-55-generic'
CC [M] /home/yangyang/fibdrv/fibdrv.o
/home/yangyang/fibdrv/fibdrv.c: In function ‘fib_sequence’:
/home/yangyang/fibdrv/fibdrv.c:38:5: warning: ISO C90 forbids variable length array ‘f’ [-Wvla]
38 | __int128 f[k + 2];
| ^~~~~~~~
MODPOST /home/yangyang/fibdrv/Module.symvers
LD [M] /home/yangyang/fibdrv/fibdrv.ko
make[1]: Leaving directory '/usr/src/linux-headers-5.8.0-55-generic'
make unload
make[1]: Entering directory '/home/yangyang/fibdrv'
sudo rmmod fibdrv || true >/dev/null
rmmod: ERROR: Module fibdrv is not currently loaded
make[1]: Leaving directory '/home/yangyang/fibdrv'
make load
make[1]: Entering directory '/home/yangyang/fibdrv'
sudo insmod fibdrv.ko
make[1]: Leaving directory '/home/yangyang/fibdrv'
sudo ./client > out
make unload
make[1]: Entering directory '/home/yangyang/fibdrv'
sudo rmmod fibdrv || true >/dev/null
make[1]: Leaving directory '/home/yangyang/fibdrv'
Passed for f( 181 ).
```
得到 Passed for f( 181 ). 確實成功執行到 f(180) 左右。
利用 gnu plot 解析輸出檔 time.txt
``` g
1 reset
2 set ylabel 'time(ns)'
3 set xlabel 'size'
4 set title 'fibdiv performance comparision'
5 set term png enhanced font 'Arial,12'
6 set key top left
7 set output 'runtime_seq.png'
8
9 plot [:][:]'time.txt' using 1:2 with linespoints linewidth 2 title 'sequential'
```
產生圖檔
觀察可的，曲線近乎線性，若只利用 __int128 的確能算出更大的費氏數列，但不能觀察到大數運算所帶來的額外計算時間成本。
因此，考慮實做出一個能計算到 f(1000) 以上的結構，以觀察其大數運算的額外效應。
### 利用 GMP 代替 `__int128`
由 [GMP Initialization Functions](https://gmplib.org/manual/Initializing-Integers) 所述，
> While n defines the initial space, x will grow automatically in the normal way, if necessary, for subsequent values stored. mpz_init2 makes it possible to avoid such reallocations if a maximum size is known in advance.
mpz 形態會自動調整所需要用到的記憶體空間，這種記憶體的重新分配，會造成不必要的效能上的浪費，因爲我們可以從前面的方法 `N(n)=n×0.69424` 得知，給予一個 x ，應該用多少 bit 去存 f(x)。
因此使用 mpz_init2 這個初始化函數,預先把會用到的記憶體大小分配好，就不會消耗到自動增加記憶體空間的時間。
因爲 user 是從 llseek 去設置要計算的費氏數列，所以我們在 f_pos 被成功設置的時候，直接計算好我們所需要的記憶體空間 bit 數量。
```c
114 if (new_pos > MAX_LENGTH)
115 new_pos = MAX_LENGTH; // max case
116 if (new_pos < 0)
117 new_pos = 0; // min case
118 file->f_pos = new_pos; // This is what we'll use now
119
120 long long temp = new_pos * 0.7 + 1;
121 if (temp < 32)
122 temp = 32;//min bit size
123 fib_result_size_bit = temp;
```
fib_result_size_bit 是一個全域變數
```c
34 static ktime_t kt;
35 static long long fib_result_size_bit = 32;
```
如此一來，我們便可利用它來初始化我們的 mpz 形態，並實做出 fib_gmp
```c
38 static char* fib_gmp(long long k)
39 {
40 //initial the mpz_t
41 mpz_t tmp[3];
42 for (int i = 0; i < 3; i++)
43 mpz_init2(tmp[0], fib_result_size_bit);
44 mpz_set_ui(tmp[0],0);
45 mpz_set_ui(tmp[1],1);
46
47 if (k == 0)
48 return tmp[0];
49 else if (k == 1)
50 return tmp[1];
51
52 int target[3] = {2, 0, 1};
53
54 for (int i = 2; i <= k; i++) {
55 int tmp = target[0];
56 target[0] = target[1];
57 target[1] = target[2];
58 target[2] = tmp;
59
60 mpz_add(tmp[target[2]], tmp[target[0]], tmp[target[1]]);
61 }
62
63 char *result = malloc(mpz_sizeinbase(tmp[target[2]], 10) + 2);
64
65 if (result == NULL){
66 printk("malloc fail\n");
67 return NULL;
68 }
69
70 mpz_get_str(result, 10, tmp[target[2]]);
71 mpz_clears(tmp);
72 return result;
73 }
```
這裏回傳的是 char 形態，代表結果已經轉換成 10 進制，在 client 端只要直接輸出就行。
在一切準備好後，進行編譯發現
```bash
yangyang@yangyang-GL552VX:~/fibdrv$ sudo make check
make -C /lib/modules/5.8.0-55-generic/build -I /usr/include/x86_64-linux-gnu M=/home/yangyang/fibdrv modules -lgmp
make[1]: Entering directory '/usr/src/linux-headers-5.8.0-55-generic'
CC [M] /home/yangyang/fibdrv/fibdrv.o
/home/yangyang/fibdrv/fibdrv.c:14:10: fatal error: gmp.h: No such file or directory
14 | #include <gmp.h>
| ^~~~~~~
compilation terminated.
make[2]: *** [scripts/Makefile.build:286: /home/yangyang/fibdrv/fibdrv.o] Error 1
make[1]: *** [Makefile:1785: /home/yangyang/fibdrv] Error 2
make[1]: Leaving directory '/usr/src/linux-headers-5.8.0-55-generic'
make: *** [Makefile:13: all] Error 2
```
出現找不到 gmp.h 這個 header file，之後發現自己傻了，在前面說明 kernel module 的特性已經提過
> kernel space 不能呼叫 user space 定義的函式例如 libc，只能使用在 kernel space 定义的那些资源受到限制的函数，例如 printk()
因此 gmp 是不能在 kernel module 使用的。
### 實做一個大數結構
實做一個大數結構，在不求效能的情況下，基於對大數做除法會非常慢，最難的部分就是把**二進制轉換到十進制的字串**，因爲這需要找到大數的餘數和商數，所以網絡上有許多直接用字串當做 10 進制的數字來做計算，但這種方式所做的算數運算比起用二進制做的算數運算時間與空間成本要來的高，但優點在於能夠直接輸出字串。
這裏經過考慮，我使用二進制的方式實做。
利用 unsigned long long 作爲儲存單位，先算出我們需要用到多少個單位來做運算和儲存。前面我們已經知道 N(n)=n×0.69424 ，可得到我們需要多少 2 進制的位元來儲存該費氏數列的結果。
因此，如下的空間就夠我們使用
```c
unsigned long long x[*offset * 0.7 / 64 + 1];
```
因爲 kernel 內不能浮點數運算，改爲
```c
unsigned long long x[*offset * 7 / 640 + 1];
```
我這裏先以實做出結果並測試爲主要目的，所以我沒把資料結構包裝成 .h 而是直接實做在 fibdiv.c 。
最主要的部分歸於 fib_BN_add ,fib_BN_sub 和 fib_BN_mul
加法
```c
136 static void fib_BN_add (unsigned long long *x,
137 unsigned long long *y,
138 unsigned long long *z,
139 unsigned long long size) //x + y = z
140 {
141 unsigned long long carry = 0;
142 __int128 tmp;
143 for (int i = 0; i < size; i++) {
144 tmp = (__int128) x[i] + y[i] + carry;
145 carry = (tmp > 0xFFFFFFFFFFFFFFFF) ? 1 : 0;
146 z[i] = tmp;
147 }
148 }
```
減法
```c
150 static void fib_BN_sub (unsigned long long *x,
151 unsigned long long *y,
152 unsigned long long *z,
153 unsigned long long size)//x - y = z
154 {
155 int borrow = 0;
156 int next_borrow;
157
158 for (int i = 0; i < size; i++) {
159 next_borrow = (x[i] < y[i]) ? 1 : 0;
160 z[i] = x[i] - y[i] - borrow;
161 borrow = next_borrow;
162 }
163 }
```
乘法，參考[第六周測驗題](https://hackmd.io/@sysprog/linux2021-quiz6#%E6%B8%AC%E9%A9%97-1)
```c
149 static void fib_BN_mul (unsigned long long *x,
150 unsigned long long *y,
151 unsigned long long *z,
152 unsigned long long size)
153 {
154 unsigned long long tmp[size];
155 unsigned long long row[size];
156 fib_BN_init(tmp, size);
157 fib_BN_init(row, size);
158 fib_BN_init(z, size);
159
160 for (int i = 0; i < size; ++i) {
161 fib_BN_init(row, size);
162
163 for (int j = 0; j < size; ++j) {
164 if (i + j < size) {
165 fib_BN_init(tmp, size);
166 __int128 intermediate = x[i] * (__int128) y[j];
167 tmp[0] = intermediate;
168 tmp[1] = intermediate >> 64;
169 fib_BN_leftshift64(tmp, i + j, size);
170 fib_BN_add(tmp, row, row, size);
171 }
172 }
173
174 fib_BN_add(row, z, z, size);
175 }
176 }
```
這三個基本運算實做出來，我們就已經能把費氏數列的結果算出來了，但最後的目的是輸出成字串，但因爲除法的復雜，這裏先利用 shift 可用於除 2的冪次的特性，實做一個可輸出成 16 進制的字串。
```c
263 static char* fib_BN_toStr16(unsigned long long *x)
264 {
265 char *result = kmalloc(buff_size, GFP_KERNEL);
266 int index = 0;
267 if(result == NULL){
268 printk("kmalloc fail\n");
269 exit(1);
270 }
271
272 while (!fib_BN_equalzero(x)) {
273 result[index] = "0123456789abcdef"[x[0] & 0x000000000000000F];
274 index++;
275 fib_BN_shift(x, 0, 4);//right shift 4 bit
276 }
277
278 int count = 0;
279 for(int i = index - 1; i > count; i--) {
280 char tmp = result[i];
281 result[i] = result[count];
282 result[count] = tmp;
283 count++;
284 }
285 result[index] = '\0';
286
287 return result;
288 }
```
這樣我們便可以開始測試，這裏先利用和上面的 sequence 一樣的眼算法。
```c
294 static char* fib_BN_sequence(long long k)
295 {
296 if (k == 0) {
297 char *a = kmalloc(buff_size, GFP_KERNEL);
298 a[0] = '0';
299 a[1] = '\0';
300 return a;
301 } else if (k == 1) {
302 char *a = kmalloc(buff_size, GFP_KERNEL);
303 a[0] = '1';
304 a[1] = '\0';
305 return a;
306 }
307
308 unsigned long long tmp[3][fib_BN_size];
309
310 for(int i = 0; i < 3; i++)
311 fib_BN_init(tmp[i], fib_BN_size);
312 tmp[1][0] = 1;
313
314 int target[3] = {2, 0, 1};
315 for (int i = 2; i <=k; i++) {
316 int s = target[0];
317 target[0] = target[1];
318 target[1] = target[2];
319 target[2] = s;
320
321 fib_BN_add(tmp[target[0]], tmp[target[1]], tmp[target[2]], fib_BN_size);
322 }
323 char *result;
324 result = fib_BN_toStr16(tmp[target[2]]);
325 return result;
326 }
```
一切準備就緒，在[排除干擾效能分析的因素](https://hackmd.io/@sysprog/linux2021-fibdrv#-Linux-%E6%95%88%E8%83%BD%E5%88%86%E6%9E%90%E7%9A%84%E6%8F%90%E7%A4%BA)過後，設定 cpu afinity 指定行程並執行。
```
sudo make all
sudo make unload
sudo make load
sudo taskset 0x1 ./client > o
python3 scripts/verify.py
```
成功輸出 f（1000）
```
Passed for f( 1001 ).
```
（這裏的 verify.py已改成驗證 16 進制的結果）
分析圖形
已經可從圖形觀察出大數運算的效應。
這裏先利用下圖與 [2019 年同學](https://hackmd.io/@zodf0055980/r1BgdLZ1r#%E5%A4%A7%E6%95%B8%E9%81%8B%E7%AE%97%E7%9A%84%E6%99%82%E9%96%93)所做的結果做對比，發現同樣 f(500) ，這邊的 2 進制大數結構會比其所實做的資料結構快，但缺點就在於，輸出的字串是 16 進制的，如果要輸出 10 進制將會增加相應的成本。
### 加入 fast doubling 演算法
從 [Fast Doubling](https://chunminchang.github.io/blog/post/calculating-fibonacci-numbers-by-fast-doubling) 可得知其已經經過許多優化的過程與結果，將其利用在本實驗。其中，對於
```c
300 static char* fib_BN_fd(unsigned int n)//fast doubling
301 {
302 if (n == 0) {
303 char *result = kmalloc(buff_size, GFP_KERNEL);
304 result[0] = '0';
305 result[1] = '\0';
306 return result;
307 }
308 unsigned int h = 0;
309 for (unsigned int i = n ; i ; ++h, i >>= 1);
310
311 uint64_t a[fib_BN_size]; // F(0) = 0
312 uint64_t b[fib_BN_size]; // F(1) = 1
313 fib_BN_init(a, fib_BN_size);
314 fib_BN_init(b, fib_BN_size);
315 b[0] = 1;
316
317 uint64_t c[fib_BN_size];
318 uint64_t d[fib_BN_size];
319 fib_BN_init(c, fib_BN_size);
320 fib_BN_init(d, fib_BN_size);
321
322 uint64_t tmp[2][fib_BN_size];
323 fib_BN_init(tmp[0],fib_BN_size);
324 fib_BN_init(tmp[1],fib_BN_size);
325
326 for (unsigned int mask = 1 << (h - 1) ; mask ; mask >>= 1) { // Run h times!
327
328 fib_BN_mul(a, a, tmp[0], fib_BN_size);
329 fib_BN_mul(b, b, tmp[1], fib_BN_size);
330 fib_BN_add(tmp[0], tmp[1], d, fib_BN_size);
331
332 fib_BN_shift(b, LEFT, 1);
333 fib_BN_sub(b, a, b, fib_BN_size);
334 fib_BN_mul(a, b, c, fib_BN_size);
335
336 if (mask & n) { // n_j is odd: k = (n_j-1)/2 => n_j = 2k + 1
337 fib_BN_assign(a, d, fib_BN_size);
338 fib_BN_add(c, d, b, fib_BN_size);
339 } else { // n_j is even: k = n_j/2 => n_j = 2k
340 fib_BN_assign(a, c, fib_BN_size);
341 fib_BN_assign(b, d, fib_BN_size);
342 }
343 }
344
345 return fib_BN_toStr16(a);
346 }
```
按照 fast doubling 的特性，能只做 log n 次費氏數列運算得到 f(n) 。因此曲線趨勢應該有 log n 的趨向。
分析結果
以次結果來說整體時間成本是比 sequence 的方法來得高許多，按理說不應該更慢才對，而爲什麼會有這種情況呢？推測是乘法的原因
#### 優化 fast doubling
回過來看
```c
149 static void fib_BN_mul (unsigned long long *x,
150 unsigned long long *y,
151 unsigned long long *z,
152 unsigned long long size)
153 {
154 unsigned long long tmp[size];
155 unsigned long long row[size];
156 fib_BN_init(tmp, size);
157 fib_BN_init(row, size);
158 fib_BN_init(z, size);
159
160 for (int i = 0; i < size; ++i) {
161 fib_BN_init(row, size);
162
163 for (int j = 0; j < size; ++j) {
164 if (i + j < size) {
165 fib_BN_init(tmp, size);
166 __int128 intermediate = x[i] * (__int128) y[j];
167 tmp[0] = intermediate;
168 tmp[1] = intermediate >> 64;
169 fib_BN_leftshift64(tmp, i + j, size);
170 fib_BN_add(tmp, row, row, size);
171 }
172 }
173
174 fib_BN_add(row, z, z, size);
175 }
176 }
```
用以下幾點優化：
* 我們發現 row 其實沒必要存在，可直接用 z 去代替從而減少使用加法的次數
* shift64 其實可以不需要用，可直接對位置做 assign。
* 在 x[i] 和 y[i] 等於 0 時，可以不需要計算。
```c
174 static void fib_BN_mul (uint64_t *x,
175 uint64_t *y,
176 uint64_t *z,
177 uint64_t size)
178 {
179 uint64_t tmp[size];
180 fib_BN_init(tmp, size);
181 fib_BN_init(z, size);
182
183 for (int i = 0; i < size; ++i) {
184 if (x[i] == 0)
185 continue;
186 for (int j = 0; j < size; ++j) {
187 if (y[i] == 0)
188 continue;
189 if (i + j < size) {
190 fib_BN_init(tmp, size);
191 __int128 intermediate = x[i] * (__int128) y[j];
192 tmp[i + j] = intermediate;
193 if (i + j + 1 < size)
194 tmp[i + j + 1] = intermediate >> 64;
195 fib_BN_add(tmp, z, z, size);
196 }
197 }
198 }
199 }
```
經過優化後，再次測試
可以看到明顯的速度提升，但比起 sequence ，還是沒有明顯優勢
Sign in with Wallet
Connect another wallet