HackMD
2022q1 Homework3 (fibdrv)
contributed by < NOVBobLee >
環境設置
- Linux 核心版本（Ubuntu）：
- 編譯器版本：
- CPU 架構和規格：
- 確認使用者身份不是
root：
- 確認 UEFI secure boot 已關閉：
- 確認
linux-headers套件已安裝：
- 檢查專案，確認可以執行（出現綠色
Passed [-]）：
獨立實驗環境，排除干擾
參照這次 K04: fibdrv 實驗說明和 KYG 學長的共筆，將干擾實驗的因素排除。
- 保留一顆 CPU 核心不被排其他行程
- 使用
taskset指定行程執行在保留的 CPU 上
- 關閉 ASLR (Address Space Layout Randomization)
關於 sudo sh -c
- 關閉 Intel 處理器的 Turbo mode
- 使用 CPU 的 performace mode
- 修改 SMP IRQ Affinity ，避免保留的 CPU 去處理 IRQ
有些 SMP IRQ Affinity 是無法從 User Space 更改的，比方說 IRQ0 - timer interrupt ，不過從 /proc/interrupts 觀察，剛好保留的 CPU 沒處理過這個 IRQ ，發生次數也不多。
除了 IRQ0 外，另外還有 IRQ2 - cascaded interrupt 的也無法更改。不過 IRQ2 尚未發生過，在更改過 default_smp_affinity 後，若 IRQ2 發生，也不會使用保留的 CPU 。
相關連結：
Is it possible to change which core timer interrupts happen on?
[patch 38/55] genirq: Introduce effective affinity mask
The irq_domain interrupt number mapping library
本共筆測試結果皆參照 用統計手法去除極端值 取樣 1000 次，取 95% 信賴區間之後再計算其平均，可參考程式碼 expt01_userkernel.c 。
測試效果
在原環境下：
在原環境下，每次測試變化大，不易觀察趨勢。
排除以上干擾後：
排除以上干擾後，每次測試的變化變小，趨勢比較容易觀察，穩定度真的上升很多，不過還是會有突波出現。
概觀 fibdrv 模組
- 編譯
fibdrv模組：
fibdrv模組資訊：
- 裝載
fibdrv模組後，可以看到他以 character deivce 的形式出現：
模組主程式 fibdrv.c ：
由 fib_sequence 計算費氏數列。
由於此模組於 VFS 中有產生 character device file ，所以可以藉由檔案操作介面進行操作，有定義的方法如下：
使用 read 讀取費氏數列的計算結果。
計算費氏數列到第幾位是由 offset 控制，可由 lseek 更改 offset 。
藉由 mutex 控制使用此模組的行程數量，從 fib_open 可以看到一次只允許一個行程使用。
write 目前沒有實際功能，之後會拿來改寫成讀取計算時間用。
加入計時器
這裡在 kernel space 使用 hrtimer 計時，而 user space 使用的是 clock_gettime ，最後可以用兩者的差計算出 system call 所佔用的時間。
- 改寫
write，加入計時器，測量在 kernel space 裡計算費氏數列的時間。
- 使用
clock_gettime計算在 user space 裡，計算費氏數列的總耗費時間。
下圖為測試結果，只有計算到第 92 位費氏數，使用的方法為方法一（費氏數列定義），可以看到計算時間正比於費氏數列的位數，符合預期的 \(O(n)\) 時間複雜度。
中間的藍色線為 kernel space 與 user space 的差，幾乎為水平線，為固定值。再做進一步測試，固定計算第零位費氏數，只改變輸入的 size 大小。在實作上是沒有使用到 copy_from_user ，所以應該不會因 size 大小而變動，藍色線應該還會是水平線。
結果符合推測，這結果在之後修改 write 時，引數 size 可以拿來傳遞其他參數。不過可以看到 kernel space 與 user space 之間的差距非常大，差不多是 500 納秒，與在 kernel space 計算所需的時間相比，成本是很高的。其原因是使用 write 這個 system call ，在 kernel space 與 user space 之間做轉換，以下為 perf 的執行結果，程式是執行 write 計算第 92 位費氏數 2,000,000 次，所以總佔時間為 95.84% ，這邊我們只看 write 裡面的成份。
費氏數列實作
方法一（費氏數列定義）
費氏數列的定義為 \(F(n+1)=F(n)+F(n-1)\) ，初始值不只一種，這裡使用的是 \(F(0)=0,\, F(1)=1\) 。
最初在模組裡所用的方法即為費氏數列的定義，先配置一個陣列，大小可容納我們指定的費氏數列所有數量，依費氏數列定義開始計算從小到大的費氏數，最終的計算即為所求。
fibdrv.c 在編譯時會跳出 -Wvla 的警告 warning: ISO C90 forbids variable length array ‘f’ [-Wvla] ，而程式碼裡也有寫 FIXME ，指出 VLA 是不允許在 Linux kernel 中使用的，想到可直接替換的方法包括改用 kmalloc 等動態記憶體配置或是只儲存必要的兩個數。
下面程式碼使用 kmalloc_array 動態配置記憶體，跟原本的程式相比，差異只有 VLA 的部份換成使用動態記憶體配置，計算的部份沒有改變。
這個程式碼只使用必要的兩個數，每次做讀取和寫入的位置剛好與奇偶相關，所以每次陣列位置都需要做判斷奇偶的運算，不過由於陣列只有兩個數，足夠放在 cache 裡，甚至是 register 裡，減少讀取記憶體的時間。
下圖為 kernel space 裡的各方法所需計算時間，可以看出長度為 2 的固定長度陣列（以下用圖中代號 fixed-la 稱呼）所需時間最少，再來是 VLA ，最後是差距很大的 kmalloc ，猜測 fixed-la 比較快是因為不用反覆去記憶體讀寫，而剩下兩個算法相同，都是用連續記憶體，那差距只在動態配置記憶體呼叫上。
相關連結：
The Linux Kernel Is Now VLA-Free: A Win For Security, Less Overhead & Better For Clang
Linux 核心設計: 記憶體管理
方法二（ Exact solution method ）
由費氏數列的定義 \(F_{n+1}=F_n+F_{n-1}\) 與 \(F_n=F_n\) ，可以得出一個矩陣關係式：
\(\begin{pmatrix} F_{n+1} \\ F_n \end{pmatrix} =\begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} F_n \\ F_{n-1} \end{pmatrix}\)
而這個關係式可以進一步推得
\(\begin{pmatrix} F_{n+k+1} \\ F_{n+k} \end{pmatrix} =\begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}^k \begin{pmatrix} F_{n+1} \\ F_{n} \end{pmatrix} \text{ 和 } \begin{pmatrix} F_{n+1} \\ F_n \end{pmatrix} =\begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}^n \begin{pmatrix} F_1 \\ F_0 \end{pmatrix}\)
這樣可以看出若要算出 \(F_n\) ，需要計算出中間矩陣的 \(n\) 次方，由於該矩陣是對稱矩陣，這裡可以使用對角化（ \(\phi_+\) 為黃金比例， \(\phi_-\) 為其另一個特徵值）。
\(\begin{split} \begin{pmatrix} F_{n+1} \\ F_n \end{pmatrix} &= \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}^n \begin{pmatrix} 1 \\ 0 \end{pmatrix} \\ &= \dfrac{1}{\phi_+ - \phi_-} \begin{pmatrix} \phi_+ & \phi_-\\ 1 & 1 \end{pmatrix} \begin{pmatrix} \phi_+ & \\ & \phi_- \end{pmatrix}^n \begin{pmatrix} 1 & -\phi_- \\ -1 & \phi_+ \end{pmatrix} \begin{pmatrix} 1 \\ 0 \end{pmatrix} \\ &= \dfrac{1}{\sqrt{5}} \begin{pmatrix} \phi_+ & \phi_-\\ 1 & 1 \end{pmatrix} \begin{pmatrix} \phi_+^n & \\ & \phi_-^n \end{pmatrix} \begin{pmatrix} 1 \\ -1 \end{pmatrix} \\ &= \dfrac{1}{\sqrt{5}} \begin{pmatrix} \phi_+ & \phi_-\\ 1 & 1 \end{pmatrix} \begin{pmatrix} \phi_+^n \\ -\phi_-^n \end{pmatrix} \\ &= \dfrac{1}{\sqrt{5}} \begin{pmatrix} \phi_+^{n+1} - \phi_-^{n+1} \\ \phi_+^n - \phi_-^n \end{pmatrix} \end{split}\)
所以最後得到 \(F_n=\dfrac{1}{\sqrt{5}}(\phi_+^n - \phi_-^n)\) ，只需要計算冪即可，比較有趣的是 \(\phi_-\) 是個負數，針對他的次方 \(n\) 的奇偶性，會對前面的 \(\phi_+^n\) 一下往上修正，一下往下修正。不過可惜的是 \(\sqrt{5}\) 是無理數（ \(\phi_+\) 和 \(\phi_-\) 都有 \(\sqrt{5}\) ），在資料型態和記憶體大小的限制下，計算一定會產生誤差，而且 \(n\) 越大越明顯。
這裡使用 double 型態，直到 \(F_{71}\) 都沒有問題，然而從下一個結果開始，誤差產生，然後擴大，可能的對應方式是增加計算精準度，針對不同位數的費氏數改變精準度。
| \(F_{72}\) | \(F_{92}\) | |
|---|---|---|
| Fibonacci number | 498454011879264 | 7540113804746346429 |
| Exact-solution method | 498454011879265 | 7540113804746369024 |
| Error | +1 | +22595 |
相關連結：
The Golden Ratio
Double precision - decimal places
IEEE Arithmetic
Why aren't the FPU registers saved and recovered in a “context switch”?
Why am I able to perform floating point operations inside a Linux kernel module?
既然 \(\sqrt{5}\) 是問題根源，是否可以繞過他作計算，觀察 \(\phi_+ = \dfrac{1+\sqrt{5}}{2}\) 還有 \(\phi_- = \dfrac{1-\sqrt{5}}{2}\) ，每個費氏數都是整數，所以計算後的分子應該可以跟前面的分母 \(\sqrt{5}\) 相消，由此可知 \(F_n\) 即 \(\phi_+^n\) 的 \(\sqrt{5}\) 前面係數。
\(\begin{split} F_1 &= \dfrac{1}{\sqrt{5}}(\dfrac{1+\sqrt{5}}{2} - \dfrac{1-\sqrt{5}}{2}) \\ &= \dfrac{1}{\sqrt{5}}\cdot\dfrac{2\sqrt{5}}{2} \\ &= 1 \end{split}\)
從上面推導可以知道，實際上可以不用計算 \(\sqrt{5}\) 。又觀察 \(\dfrac{1+\sqrt{5}}{2}\) 的冪可以發現，最終都可以變成 \(\dfrac{a+b\sqrt{5}}{2}\) 的形式，所以可以得到：
\(\begin{split} (\dfrac{1+\sqrt{5}}{2})^{k-1}(\dfrac{1+\sqrt{5}}{2}) &= \dfrac{a+b\sqrt{5}}{2}\cdot\dfrac{1+\sqrt{5}}{2} \\ &= \dfrac{(a+5b)+(a+b)\sqrt{5}}{4} \\ &= \dfrac{(\dfrac{a+5b}{2})+(\dfrac{a+b}{2})\sqrt{5}}{2} \end{split}\)
以這個方式重複計算到 \(n\) 次方，然後最終 \(\dfrac{a+b}{2}\) 即為解。
因為在運算中用到 (a + 5 * b) ，結果比之前的方法更早發生 overflow ，在計算 \(F_{91}\) 時發生溢位（前面的方法溢位在 \(F_{93}\) ），那麼就試著避免先相加或相乘的情況，利用在 quiz2 裡的測驗一學到的技巧。
\(\big(a,\, b\big) \rightarrow \big(\dfrac{a+5b}{2},\, \dfrac{a+b}{2}\big) = \big(\dfrac{a+b}{2}+2b,\, \dfrac{a+b}{2}\big)\)
這樣溢位延後到 \(F_{92}\) 發生，而在計算時間上，第一改寫結果的計算時間比 fixed-la 高，比 vla 少；而第二次的改寫結果降到跟 fixed-la 一樣，猜測第一個改寫結果所需時間較高，可能跟乘法運算有關係。
方法三（ Fast doubling method ）
以上的時間複雜度都是 \(O(n)\) ，而接下來的第三個方法 Fast doubling method 可以降到 \(O(log\,n)\) 。
由 \(\begin{pmatrix} F_{n+k+1} \\ F_{n+k} \end{pmatrix} = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}^{k} \begin{pmatrix} F_{n+1} \\ F_{n} \end{pmatrix}\) 可以推導出：
\(\begin{split} \begin{pmatrix} F_{2n+1} \\ F_{2n} \end{pmatrix} &= \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}^{2n} \begin{pmatrix} 1 \\ 0 \end{pmatrix} \\ &=\begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}^{n-1} \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}^{n-1} \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} 1 \\ 0 \end{pmatrix} \\ &=\begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}^{n-1} \begin{pmatrix} F_2 & F_1 \\ F_1 & F_0 \end{pmatrix} \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}^{n-1} \begin{pmatrix} F_2 & F_1 \\ F_1 & F_0 \end{pmatrix} \begin{pmatrix} 1 \\ 0 \end{pmatrix} \\ &=\begin{pmatrix} F_{n+1} & F_n \\ F_n & F_{n-1} \end{pmatrix} \begin{pmatrix} F_{n+1} & F_n \\ F_n & F_{n-1} \end{pmatrix} \begin{pmatrix} 1 \\ 0 \end{pmatrix} \\ &=\begin{pmatrix} F_{n+1} & F_n \\ F_n & F_{n-1} \end{pmatrix} \begin{pmatrix} F_{n+1} \\ F_n \end{pmatrix} \\ &= \begin{pmatrix} F_{n+1}^2 + F_n^2 \\ F_n(F_{n+1} + F_{n-1}) \end{pmatrix} \\ &=\begin{pmatrix} F_{n+1}^2 + F_n^2 \\ F_n(2F_{n+1} - F_{n}) \end{pmatrix} \end{split}\)
從推導出的關係式可以看出位數 \(n\) 可以成兩倍跳，但是單只用此關係式只能算出位數 \(n\) 為 \(2\) 的冪左右的費氏數，還需要作一些應變。現在觀察一個數的二進位表示，其組成跟 \(2\) 的冪有關，這裡用 \(13\) 舉例。
| 最高位數 | 十進位表示 | 二進位表示 | 目標的二進位餘下部份 | 動作需求 |
|---|---|---|---|---|
| 目標 | 13 | 1101 | ||
| 0 | 0 | 0 | 1101 | |
| \(2^0\) | 1 | 1 | 101 | \(\times2+1\) |
| \(2^1\) | 3 | 11 | 01 | \(\times2+1\) |
| \(2^2\) | 6 | 110 | 1 | \(\times2\) |
| \(2^3\) | 13 | 1101 | \(\times2+1\) |
由上表看到 \(13\) 的二進位組成方式，剛開始由 \(0\) 開始，若下一個位數是 1 ，則需要 \(\times2+1\) ，所以 \(0\) 才會變成 \(1\) ；反之，下一個位數為 0 ，則只需要 \(\times 2\) ，所以 \(3\) 會變為 \(6\) ，原理跟短除法一樣。
\(\begin{matrix} F_0 \\ F_1 \end{matrix} \xrightarrow[]{\times2+1} \begin{matrix} F_1 \\ F_2 \end{matrix} \xrightarrow[]{\times2+1} \begin{matrix} F_3 \\ F_4 \end{matrix} \xrightarrow[]{\times2} \begin{matrix} F_6 \\ F_7 \end{matrix} \xrightarrow[]{\times2+1} \begin{matrix} F_{13} \\ F_{14} \end{matrix}\)
可以看出每一步都會用到 \(2\) 倍的關係式，而加 1 的部份可以用定義達成，由以上的說明可以寫出一個 bottom-up 的 fast doubling method ：
從測試結果可以看到 fast doubling method 的時間變化疑似水平線，跟之前的方法相比，在計算較大位數的費氏數時，都是 fast doubling method 以更少的時間完成計算。
fast doubling method 迴圈有二，第一個要找最高位數的 bit ，此法是從 long long 的第 62 位開始往右尋找；第二的迴圈則是 fast doubling method 的本體，迭代次數是看總共由幾個位數組成，舉例 13 是由 4 個位數組成，所以會迭代四次。
另外， mask 目前不需要使用到 unsigned long long 的長度，使用 unsigned int 長度即可應付到 \(F_{4294967295}\) （如果記憶體大到足夠使用的話）。
從上面的實驗結果可以看到，深藍色線從 \(F_3\) 開始會進入 fast doubling method 計算，包括第一個迴圈和第二個迴圈。而在第二個迴圈的迭代次數是看組成位數， 3 的組成位數為 2 所以會迭代兩次，然後到最後計算的 \(F_{92}\) ， 92 的組成位數是 6 ，可以明顯看出 fast doubling method 的結果大部分都是在跑第一個迴圈，所以現在針對第一個迴圈去做最佳化。
上面是測試最高位數 bit 的起始點，括號中為起始位數，除了 62 使用 unsigned long long ，其餘都使用 unsigned int 長度。現在是算到 \(F_{92}\) ，所以就測到起始位數 6 ，而這個改變至少可以使計算時間減少一半。
而找尋最高位數 bit 的方法還有其他方式，接下來測試 clz/ffs (count leading zeros/find first set) ，這裡使用 GCC builting-function __builting_clz 、 Linux 核心使用的 fls （方向跟 wiki 使用的相反，方向是從 LSB 往 MSB ，回傳的是位置 (1-based) ），其實作方式有三種（見相關連結），可以用 objdump -dS fibdrv.ko 看是使用哪一種，這裡使用的是 x86 硬體的 bsr 指令；同樣從反組譯的輸出看， __builtin_clz 也是用硬體的 bsr 指令。
Linux 核心中的 __ffs
Linux kernel 裡 (x86) 的 __ffs 雖然名稱跟 wiki 的一樣，不過他的 first 指的是從 LSB 開始往 MSB 的方向，換句話說是在做 ctz (count tailing zeros) ，而他的實作有三種，一種是使用 GCC 的 builtin-function __builtin_ctzl ，一種是使用 x86 指令 bsf ，而最後一種使用二元搜尋法。
在第一次測試中，直接呼叫 __ffs ，結果導致 soft lockup ，最終是用 Magic SysRq key 強制關機，跟核心有關的東西還是小心為上。
相關連結：
Softlockup detector and hardlockup detector (aka nmi_watchdog)
__ffs 使用 bsf 位置： arch/x86/include/asm/bitiops.h
__ffs 使用二元搜尋法位置： include/asm-generic/bitops/__ffs.h
__ffs 使用 __builtin_ctzl 位置： include/asm-generic/bitops/builtin-__ffs.h
從測試可以看到 __builtin_clz 和 fls 跟起始位數 6 的差不多快，剛開始在 \(F_{20}\) 前是 __builtin_clz 與 fls 較快，然後到 \(F_{31}\) 的過程中三者不相上下，之後變成起始位數 6 較快。原因是使用迴圈的次數減少，當在找 \(F_{32}\) 最高位數時只需要右移一次，甚至在 \(F_{64}\) 之後不需要右移。但是我們不會每次都知道該使用什麼起始位數，在計算多種位數的費氏數時，使用 __builtin_clz 或 fls 會更佳。
相關連結：
Calculating Fibonacci Numbers by Fast Doubling
fls 使用 bsr 的位置： arch/x86/include/asm/bitops.h
fls 使用二元搜尋法的位置： include/asm-generic/bitops/fls.h
fls 使用 __builtin_clz 的位置： include/asm-generic/bitops/builtin-fls.h
在 \(F_{92}\) 之前就用上面這個測試當一個小結尾，下一個目標是計算 \(F_{93}\) 以上的費氏數。
計算 \(F_{93}\) 之後的費氏數列 - 大數計算
在計算 \(F_{93}\) 之前要先把 fibdrv.c 的限制提高，這是由於之前 long long 的最高上限只放得下 \(F_{92}\) ，所以在 lseek 的時候會檢查 offset ，如果輸入的參數 offset 大於 92 就會被限制成 92 。
在之前的測試都只計算到 \(F_{92}\) ，是因為計算 \(F_{93}\) 會發生溢位，所以就不放在測試範圍內了。而溢位的原因是型態的限制， Linux 是使用 LP64 ，所以 long long 跟 long 位元是一樣長的，都是 64 bits ，從下表可以看到 long long 跟 unsigned long long 可以存下費氏數的極限，分別是 \(F_{92}\) 跟 \(F_{93}\) ，再多就存不下了，之後就必須靠 GCC 擴充或是大數計算了。
| \(F_{92}\) | \(F_{93}\) | \(F_{94}\) | |
|---|---|---|---|
| Fibonacci | 7540113804746346429 | 12200160415121876738 | 19740274219868223167 |
LONG_MAX |
9223372036854775807 | 9223372036854775807 | 9223372036854775807 |
ULONG_MAX |
18446744073709551615 | 18446744073709551615 | 18446744073709551615 |
相關連結：
The first 300 Fibonacci numbers, factored
The first 500 Fibonacci numbers
參考 KYG-yaya573142 學長的報告，這裡大數計算使用類似的策略，用 uint32_t 儲存資料，在使用加法、乘法、左移運算等會溢位的運算時，轉換成 uint64_t 來計算，最後在依溢位的情況（第 32 個至第 63 個位元）決定是否改變 uint32_t 的數量和儲存方式。
結構 fbn
在 KYG-yaya573142 的實作中，很多方法都是適用於 general case 的，比方說可以計算負數，所以在結構裡會有 sign 這個成員，但在費氏數列中，以我們列舉的方法裡是不會遇到負數的（除了 exact solution method ），所以說有些部份我們可以針對費氏數的計算作特化。因此，這裡使用的結構為
其儲存方式可以用下圖表示，這裡左側為以十六進位表示的目標數值，總共需要 6 bytes 來儲存，可見只用單個 uint32_t 或 uint64_t 都不夠，所以使用多個 uint32_t 來儲存，一個 uint32_t 可以儲存 4 bytes ，即圖右邊每個格子儲存 8 個十六進位的數字。而陣列方向則是順著 x86 的 little endian ，每次新的元素是產生在右邊。
其數學式可以由 \(\text{Big Number} = a_0 + a_1 2^{32} + a_2 2^{2\cdot32}+ a_3 2^{3\cdot32} + \dots + a_n 2^{n\cdot32}\text{ , }\forall\,a_k\in[\,0,\,2^{32} - 1]\) 來表示， \(a_0\) 即陣列第零個元素，以此類推。
加法 fbn_add
而實際計算方式可以從加法開始，而加法也是費氏數列的主要運算之一，先介紹一些用到的函式：
fbn_swap是交換兩個fbn指標
fbn_resize是調整fbn裡的num陣列長度
fbn_swap_content是交換兩fbn所持有的內容，目的是取代比較費工的fbn_copy。在許多大數運算上會指定回傳的fbn，在運算過程中計算的結果可能會在fbn的暫時變數（運算完即被釋放）中，此時使用fbn_copy比較浪費，因為不需要兩個內容都相同，只需要將正確結果換進指定回傳的fbn內即可。
再來解釋加法 fbn_add ，剛開始挑選讓 a 為擁有陣列比較長的 fbn ，以利之後的程式碼簡潔，然後調整回傳結果的 c 的陣列長度，保持跟 a 是一樣長的。
第二步開始進行加法，可以對照下面的數學表示，第一階段是 a, b 都擁有相對應的元素，所以在 bcabinet 的加法裡兩者都會參與，而計算上因為要保留溢位的部份，所以提昇成 u64 的型態來相加，加完以後 bcabinet 可以分為高位元和低位元兩個部份，低位元的部份可以直接存進 c 裡，而 bcabinet 的高位元部份代表現在陣列元素 u32 裡相加溢位的部份，換句話說是進位，影響的部份在下一個陣列元素，所以將之右移 32 位，準備跟下一次迴圈裡的陣列元素相加。
\(A = a_0 + a_1 2^{32} + a_2 2^{2\cdot32} + \dots \\ B = b_0 + b_1 2^{32} + b_2 2^{2\cdot32} + \dots \\ \text{If } a_0 + b_0 > (2^{32}-1)\text{ , then } a_0 + b_0 = (2^{32}-1) + K\text{ , where }K\in[\,1,\,2^{32}-1] \\ \text{i.e. }a_0 + b_0 = 2^{32} + (K-1) \\ \text{So, there's a carry bit from }2^{32}. \\ (a_1 + b_1 + 1)\cdot2^{32} = \dots\)
加法的第二階段跟第一階段相同，不過 b 已經沒有相對應的元素，所以 bcabinet 只需要跟 a 做相加。在最後 a 陣列對應的元素都相加完畢，如果 bcabinet 還有餘數，代表說還需要進位，但 c 的陣列已經到盡頭了，所以還要再加長 c 的陣列長度，將剩餘的 bcabinet 存入，這樣加法就完成了。
在加法的過程中，若 a 跟 c 為同一個指標也是沒有問題的，換句話說可以拿來做 addition assignment 運算（ a += b ），這是由於寫進 a 的某陣列元素後，代表該位數的相加已經結束，不會再回去使用該元素。順代一提， bcabinet 是來自 Biosafety Cabinet 給的印象，由於操作上會有生化危險，所以必須在一個櫃子裡操作。現在生活中充斥著新冠肺炎病毒 Covid-19 ，這裡變數名稱就使用跟病毒有關的 P4 實驗室器具。
大數計算版本的費氏數列 (1) fbn_fib_defi
完成加法後就可以寫出基本的費氏數列計算（使用定義的迭代方式）：
十進位字串輸出 fbn_print
此函式引入 KYG-yaya573142 的實作，原理跟上面 fast doubling method 說明 13 的二進位組成一樣。
改寫 read ，測試大數計算版本的費氏數列計算結果：
與 The first 500 Fibonacci numbers 對照 \(F_{92}\) 之後的結果，確認無誤。
左移運算 fbn_lshift
再來是準備寫 fast doubling method 的大數計算版本，這需要實作左移運算、減法和乘法，那麼先從左移運算開始。
左移運算有實作兩個版本，一個是移動 32 位元以內（超過會取 modulus ）的版本，另一個是通用版本，不過在 fast doubling method 中只會用到左移一位的運算，所以這邊就不放上通用版本了。
減法 fbn_sub
減法運算跟加法 fbn_add 原理一樣，由於沒有使用 sign 這個成員，所以不能使用轉換正負號的方法。最後會將 num 的高位數元素為零的部份截掉，如果最終運算結果是零，則 fbn 的長度還是需要維持為 1 ，不過這在計算費氏數列時不會遇到。
乘法 fbn_mul
目前使用比較直觀的 long multiplication ，跟 fbn_lshift 一樣，因為有使用 fls ，所以要避免遇到 num 高位數元素為零的情況。
過程中 bcabinet 使用 u64 型態解決溢位問題，當每個值處於 u32 最大極限的狀態下， u64 依然可以容納其結果，請見下面數學關係式：
\(\text{Suppose that }a = b = c = bcabinet = 2^{32}-1 \\ \begin{split} \text{Then, } a\cdot b + c + bcabinet &= (2^{64} - 2\cdot 2^{32} + 1) + 2\cdot(2^{32} - 1) \\ &= 2^{64} - 1 \leq (2^{64} - 1) \end{split}\)
大數計算版本的費氏數列 (2) fbn_fib_fastdoubling
有了以上的運算，現在可以實作 fast doubling method 的大數計算版本。
效能測試
下圖為測試兩種方法 fbn_fib_defi （左）和 fbn_fib_fastdoubling （右），分別在 user space （綠線）與 kernel space （紫線）中所花費的時間，而靠在水平線上的是 copy_to_user 加上 kernel/user 空間轉換的時間（藍線）。
兩種方法相比， fast doubling method 稍快（ kernel 是包括 fbn_fib_fastdoubling 和 fbn_print ），但看不出 \(O(log\,n)\) 的趨勢，反而跟左圖趨勢相似，推測主要時間都被 fbn_print 給主導，以下測試將 fbn_print 移到藍線上，從測試結果可以看出推測是正確的。
從以上測試中得知， fbn_print 影響過大。然後在計算 \(F_{1000}\) 以下的費氏數時，大部分時間可以由 kernel space 主導，在之後的測試中，將以測試 kernel space 時間為主。下圖為三個函式在 kernel space 中的執行時間。
從以上測試中，計畫之後的效能改善分為兩個部份：
fbn_print：此函式目的為輸出fbn結構持有數字所代表的十進位表示，影響總體花費時間最大。- 以
fbn_fib_fastdoubling為主的大數運算：在初步實作大數運算的過程中，有想到不少可以改善的空間，之後將以此函式相關部份最佳化。
改善效能
Radix conversion by performing short division
在 fbn_print 中主要計算在裡面的雙迴圈，分別對 bits 和 digits （十進位）走迴圈，其目的是 radix conversion ，從二進位轉換成十進位，原理是二進位的組成方式（見 13 的二進位組成方式）。而另外還有其他種轉換方式，比方說使用短除法，藉由不斷從除法得到的餘數，剛好會是轉換目標十進位表示的每一位數字（順序由低位到高位），所以首先需要一個有效率的大數除法。
舉個例子，將 0x80 除以 10 。一個十六進位數字可以存 4 個位元，所以上限為 0xf ，那麼其最高可以容納一個 10 (dec) 在裡面，如此我們可以以一個十六進位數字為單位進行除法，如下面數學所表示：
\(\begin{split} 0\text{x}80 &= 8\cdot 2^4 + 0\cdot 2^0 \\ &= (0\cdot 10 + 8)\cdot 2^4 + 0\cdot 2^0 \\ &= 0\cdot 10\cdot 2^4 + (8\cdot 2^4 + 0)\cdot 2^0 \\ &= 0\cdot 10\cdot 2^4 + (12\cdot 10 + 8)\cdot 2^0 \\ &= 0\cdot 10\cdot 2^4 + 12\cdot 10\cdot 2^0 + 8\cdot 2^0 \\ &= (0\cdot 2^4 + 12\cdot 2^0)\cdot 10 + 8 \\ &= 0\text{xc}\cdot 10 + 8\ \end{split}\)
剛開始先對高位數的 \(8\) 進行除法，得到商 \(0\) 餘數 \(8\) ，但由於是高位數，所以該餘數的實際值還要乘上他的底 \(2^4\) ，即實際值為 \(8\cdot 2^4\) ，也就是說還有 10 在裡面，所以接下來將他加進比較低的位數（此時兩者基底相同，皆為 \(2^0\) ），變成 \((8\cdot 2^4+0)\) ， 繼續進行除法，得到商 \(12\) 餘數 \(8\) ，此時沒有更低的位數，餘數也比 10 小，所以除法已經完成。如此一來，最終可以算出商 0xc （依舊可以由十六進位表示）與餘數 8 。
此除法可以套用在大數計算上，不過除數有個限制，就是必須可以被 u32 給容納（ num 元素使用 u32 型態）。另外，大數除法需要走過陣列每一個元素，將會耗費不少時間，所以能做越少次越好，參考 sysprog21/bignum 的作法，可以將除數拉高到 \(10^9\) ，這樣總體流程變為
- 大數除法除以 \(10^9\)
- 得到的餘數再做短除法，輸出 9 個十進位數字
- 回到 1.
而高位加低位的除法剛好在 x86-64 裡有指令 div 可以使用，其中 l 是代表 doub(l)eword (4 bytes) ，每個變數皆為 32 bits ， d 是除數， q 為商， r 為餘數，比較特別的是被除數是 64 bits ，由 high, low 組成，所做的事情以下面數學表示：
\(high\cdot2^{32}+low = q\cdot d + r\text{ , where }0\leq q\leq(2^{32}-1)\text{ and }0\leq r<d\)
從上面數學式可以看出這個指令有 integer overflow 的可能，比方說 high = 1, d = 1 的情況，此時 q 會存不下超過 32 bits 的商，不過在這使用的演算法中不會遇到除法溢位，由於 high 是高位數的餘數，上限為 \(10^9-1\) ，此時除法不會遇到溢位問題，請見下面數學：
\(\text{Suppose that }high = 10^9-1,\,low = 2^{32}-1,\,d = 10^9 \\ \begin{split} q &= \lfloor\dfrac{high\cdot2^{32}+low}{d}\rfloor \\ &= \lfloor\dfrac{(10^9-1)\cdot2^{32}+(2^{32}-1)}{10^9}\rfloor \\ &= \lfloor\dfrac{2^{32}\cdot10^9-1}{10^9}\rfloor \\ &= \lfloor\dfrac{(2^{32}-1)\cdot10^9+(10^9-1)}{10^9}\rfloor \\ &= 2^{32}-1 \leq (2^{32}-1) \end{split}\)
以下為步驟 1. 的程式碼，由於想減少呼叫 krealloc ，所以沒有使用 fbn_resize ， obj 不會因為出現零元素就調整長度，這樣的話，需要另外使用 nonzero_len 紀錄現在實際值所需要的長度：
而在 Linux 核心中當然也有 radix conversion 的需要，比方說 printf ，方法也是使用短除法取餘數，不過並不會使用除法和模，取而代之使用定點數，舉個例子， r 除以 10 的商可以轉換成：
\(\begin{split} \dfrac{r}{10} &= \dfrac{r}{10}\cdot\dfrac{2^{15}}{2^{15}} \\ &= (r\cdot\dfrac{2^{15}}{10})\cdot\dfrac{1}{2^{15}} \\ &\approx (r\cdot0\text{xccd})\cdot\dfrac{1}{2^{15}} \end{split}\)
可以看到他將除法轉換成乘法，由於是取近似值，所以還是會有些誤差存在，不過他保證在 \(r<16389\) 的情形下輸出是正確的。這邊將 Linux 核心中的 put_dec 改寫成適合這裡的版本。
程序大概是這樣，先除以 \(10^4\) ，然後除以 \(10\) 四次，分別從餘數得到 4 個 digits （基底為 \(10^0\) ），然後一樣步驟再算出 4 個 digits （基底為 \(10^4\) ），再來因為 put_dec 是輸入除以 \(10^9\) 之後的餘數，所以現在只剩一個 digit （基底為 \(10^8\) ），直接轉換成字元即可。
再來將兩個函式組合起來，先是用 fbn_divten9 做 \(10^9\) 的大數除法，然後使用 put_dec 從餘數算出 9 個 digits ，輸出字串，重複這兩個步驟直到整個的短除法結束。
測試結果真的非常令人滿意，原因是 v0 是疊加二進位的 bit ，每次計算簡單，一個加法和一個 if 裡的減法，但總共計算次數要 bits 數量乘上 digits （十進位）數量（準確說是 kmalloc 的字串長度），這在費氏數越大時，其計算量非常可觀。而 v1 是使用短除法直接計算 digit ，每次計算雖然是比較複雜，除法、乘法和加法（ fbn_divten9 接 put_dec ），不過次數較少， fbn 長度乘上 digits 數量，總共計算量是少上許多。現在使用 v1 執行下面的實驗（ \(F_0\) 算到 \(F_{1000}\) 每次 1000 次），大概 5 秒內可以結束，這在之後開發與測試的效率將可以大幅提昇。
在 kernel space 中三個函式執行時間比較，現在 fbn_printv1 時間比 fbn_fib_fastdoubling 少，不過從趨勢看，之後 fbn_printv1 將會超過 fbn_fib_fastdoubling ，然後 copy_to_user 則是相較之下比較微小的漸增趨勢。
相關連結：
大數除法： sysprog21/bignum
使用定點數做 radix conversion ： put_dec
Reciprocal Multiplication, a tutorial
Ratio of Bits to Decimal Digits
Fast doubling method without subtraction (v0 -> v1)
原本的 fast doubling method 是 \(\begin{pmatrix}F_{n+1} \\ F_n\end{pmatrix}\rightarrow\begin{pmatrix}F_{2n+1} \\F_{2n}\end{pmatrix}\) ，這其中還有幾個地方還能改進，第一個：由於是從 \(\begin{pmatrix}F_1 \\ F_0\end{pmatrix}\) 開始，這在原本的演算法中，在第一次迴圈乘以二的部份會是做白工（ \(\begin{pmatrix}F_1 \\ F_0\end{pmatrix}\xrightarrow[\times2]{}\begin{pmatrix}F_1 \\F_0\end{pmatrix}\) ），只有加一的部份是有效的（ \(\begin{pmatrix}F_1 \\ F_0\end{pmatrix}\xrightarrow[+1]{}\begin{pmatrix}F_2 \\F_1\end{pmatrix}\) ）。第二個：有使用減法，在大數運算的減法中，會將沒有必要的零元素給去除，持續修正 fbn 的陣列長度，這在演算法中會增加呼叫 krealloc 的機會。第三個：若目標是 \(F_k\) ，則每次最終都會多計算一個附加的 \(F_{k+1}\) 。
針對這幾點，這裡使用稍微不同的數學等式，一樣利用 \(\begin{pmatrix} F_{n+2} \\ F_{n+1} \end{pmatrix} = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} F_{n+1} \\ F_{n} \end{pmatrix}\) 推導，不同在於這次是使用 \(2n-1\) 次方：
\(\begin{split} \begin{pmatrix} F_{2n} \\ F_{2n-1} \end{pmatrix} &= \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}^{2n-1} \begin{pmatrix} F_1 \\ F_0 \end{pmatrix} \\ &=\begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}^{n-1} \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}^{n-1} \begin{pmatrix} F_1 \\ F_0 \end{pmatrix} \\ &=\begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}^{n-1} \begin{pmatrix} F_2 & F_1 \\ F_1 & F_0 \end{pmatrix} \begin{pmatrix} F_n \\ F_{n-1} \end{pmatrix} \\ &=\begin{pmatrix} F_{n+1} & F_n \\ F_n & F_{n-1} \end{pmatrix} \begin{pmatrix} F_n \\ F_{n-1} \end{pmatrix} \\ &=\begin{pmatrix} F_n(F_{n+1}+F_{n-1}) \\ F_n^2 + F_{n-1}^2 \end{pmatrix} \\ &= \begin{pmatrix} F_n(2F_{n-1}+F_n) \\ F_n^2 + F_{n-1}^2 \end{pmatrix} \end{split}\)
最終變成 \(\begin{pmatrix}F_n \\ F_{n-1}\end{pmatrix}\rightarrow\begin{pmatrix}F_{2n} \\F_{2n-1}\end{pmatrix}\) ，來比較一下兩種等式的不同之處，以目標 \(F_{13}\) 為例：
原：\(\begin{matrix} F_1 \\ F_0 \end{matrix} \xrightarrow[\times2+1]{} \begin{matrix} F_2 \\ F_1 \end{matrix} \xrightarrow[\times2+1]{} \begin{matrix} F_4 \\ F_3 \end{matrix} \xrightarrow[\times2]{} \begin{matrix} F_7 \\ F_6 \end{matrix} \xrightarrow[\times2+1]{} \begin{matrix} F_{14} \\ F_{13} \end{matrix}\)
新：\(\begin{matrix} F_1 \\ F_0 \end{matrix} \xrightarrow[]{\times2+1} \begin{matrix} F_3 \\ F_2 \end{matrix} \xrightarrow[]{\times2} \begin{matrix} F_6 \\ F_5 \end{matrix} \xrightarrow[]{\times2+1} \begin{matrix} F_{13} \\ F_{12} \end{matrix}\)
可以看出使用新的等式，總體少了一次計算，是第一次迴圈的部份，之後的計算則是跟原本的相同。第二，新的等式最終會計算到 \(F_{13}\) ，不會計算出多餘的 \(F_{14}\) 。第三，減法的部份變成使用加法。
兩種方法進行比較測試，平均減少 260 ns 左右，計算 \(F_{1000}\) 時間從 4112.2986 (v0) 減至 3853.0461 ns (v1) ，減少幅度約為 6% 。
Reduce unnecessary computations (v1 -> v2)
DIV_ROUNDUP32：結合ROUNDUP32巨集與除以 32 的計算，減少不必要的計算，影響的函數有fbn_lshift32和fbn_mul。
- 修改
fbn_lshift32：原函式可以左移 32 位元（含）以下，改為 31 位元（含）以下。
測試結果如下，由於計算費氏數越到後面，使用 fbn_lshift 與 fbn_mul 的次數就越多，所以時間的差距會拉大。計算 \(F_{1000}\) 時間從 3853.0461 (v1) 減少至 3653.0772 ns (v2 - divroundup32) ，減少幅度約 5% 。
Lazy krealloc (memory pool)
類似 Lazy allocation ，只有在需要的時候才會呼叫配置記憶體的函式，目的是減少呼叫 krealloc 的次數。這裡將配置陣列元素的單位增加至 4 個 u32 元素，只有當配置的陣列元素不足時才會呼叫 krealloc ，如此配置的陣列元素只增不減，可以減少呼叫 krealloc 的需求。
首先要修改 fbn 結構，新的成員 cap 是配置的總陣列長度；而新的 len 成員表示有效元素長度，指的是在該長度內元素儲存的值是有效的，因為將配置的單位加大成 4 ，有些元素是沒有意義的零，而 len 就是區分元素是否有效的一個成員，不過要注意若 fbn 要表示零，是使用 len 長度為零來表示：
而 lazy krealloc 的主角 fbn_resize ，將配置的單位加大至 4 個元素，且只在配置的陣列長度不足時才會呼叫 krealloc 。
一樣是測試 fast doubling method ，比較的對象是上一次的 v2 ，計算的費氏數越大，效益也越大。在計算 \(F_{1000}\) 時，時間從 3653.0772 (v2) 減少至 3603.3485 ns (v3) ，幅度約為 1% 。（若有突起，下降比紫線快，綠線總體相較比較平滑， why? ）
另一個方案是先計算最終費氏數需要多少陣列元素，只在最初的時候配置所有陣列元素，所以途中完全不需要使用 krealloc ，而如何計算多少陣列元素，可以參考此式 binary digits = \(\lfloor log_2(F_n)\rfloor + 1\) ，利用總二進位位數換算成 u32 的元素個數，但有兩個問題，一是需要浮點數計算，而 \(n\) 越大，誤差也越大，二是需要先計算出 \(F_n\) ，亦即我們的最終計算目標。
使用 perf 分析
使用 perf 觀察 fbn_fib_fastdoubling 的執行效能，將 fib_read 改成以下樣子，重複計算 fbn_fib_fastdoubling 200000 次，讓該函式變成時間主導項。由於有 lazy realloc 的關係，所以每次迴圈必須重新初始化 fib 這個變數，這裡可以用 caller-callee 關係知道 fbn_alloc 和 fbn_free 是哪一部份的程式，進一步排除該因素。
由 perf 輸出可以知道 fbn_fib_fastdoubling 每部份計算的時間佔比，其中 fbn_mul 佔最多，佔了 84.77% / 95.41% ，其他依序為 fbn_add, fbn_alloc, fbn_free, fbn_lshift31 ，佔比都相對低很多，很清楚可以知道，對 fbn_mul 最佳化效果會是最明顯的。
Karatsuba multiplication
\(U = U_0 + U_1 2^N \\ V = V_0 + V_1 2^N \\ \text{Let }z_0 = U_0V_0,\,z_1 = (U_0 + U_1)(V_0 + V_1),\,z_2 = U_1V_1 \\ \text{Then, }UV = z_0 + (z_1 - z_0 - z_2)2^N + z_1 2^{2N} \\ \text{Let }z_1' = (U_0 - U_1)(V_1 - V_0) = (U_1 - U_0)(V_0 - V_1) \\ \text{Then, }UV = z_0 + (z_1' + z_0 + z_2)2^N + z_1 2^{2N}\)
square instead of mul
(save almost half of multiplications)
