contributed by < blueskyson
>
在 fibdrv.c 的 fib_sequence
以 long long
來處理費氏數列的值,運算到 時會 overflow,因此必須以能夠儲存更大數值的資料結構取代 long long
。
struct BigN
取代 long long
首先在 client.c 以及 fibdrv.c 定義結構 struct BigN
並給予別名為 u128
:
我一開始用 u128
表示正整數 ,但是在我無法基數轉換。
參考 hankluo6 的作法,讓 lower 最多只儲存 ,因此 u128
表示的值為 。如此一來數值由 lower
進位到 upper
的規則會類似於 進位,可以表示的數值範圍為 到 ,規律如下:
decimal | upper | lower |
---|---|---|
0 | 0 | 0 |
1 | 0 | 1 |
2 | 0 | 2 |
0 | 99…99 | |
1 | 0 | |
1 | 1 | |
1 | 99…99 | |
2 | 0 |
如此就能夠很輕易的利用 printf
印出字串如下:
套用 u128
改寫的後的 fib_sequence
如下:
接下來再改寫 fib_read
,原本的 fib_read
函式內部直接回傳 的值,亦即透過 read(fd, buf, 1)
的回傳值來告訴 client.c 的值,其型態為 ssize_t
。追溯 ssize_t
的定義:
觀察上述定義以及參閱 Integer-Types 得知 long
為 32 或 64 位元,我測試了自己的電腦得知我的 long
為 64 位元,亦即 ssize_t
為 64 位元。
顯然 ssize_t
傳遞 以後的數值會 overflow,應該改用 copy_to_user
API,將 的資料,複製到 read(fd, buf, 1)
的 buf
中,在 client.c 做基數轉換。改寫後的 fib_read
如下:
接下來改寫 client.c:
測試時我發現程式永遠都顯示 的結果,在 fib_read
中使用 printk("%lld\n", *offset);
然後用 dmesg
觀察紀錄,發現 *offset
會被限制不能超過 92,然後我在 fibdrv.c 中發現巨集:
它在 fib_device_lseek
中限制 lseek
的 offset
。將此巨集改為 100 後即可正常執行到 。執行 python scripts/verify.py
成功通過驗證。
此 BigN 實作最高可以計算 。
截至目前為止的實作在 53f5ddf。
struct BigN
使其能夠計算 改寫 struct BigN
的資料結構,取名為 ubig
,把 upper
, lower
擴增為 cell[size]
如下:
上方的 BIGNSIZE
巨集用來控制 cell
陣列的長度,當 BIGNSIZE
為 12,ubig
即為一個擁有 768 位元的 cell
用以儲存費氏數列的值,可以表示的數值範圍為 到 。cell[0]
儲存最低位、cell[11]
儲存最高位。
接著改寫 fib_sequence
:
為了不讓 client.c 過於複雜,我在 fibdrv.c 就將 BigN 轉換為字串,然後直接用 copy_to_user
回傳 的字串,如此一來不管 的資料結構長的如何,client.c 都能使用同一套程式碼輸出,同時也利於分析轉換字串的效能。
以下函式在 fibdrv.c 中,負責將 a
轉為字串複製到 buf
中:
然後改寫 fib_read
:
最後改寫 client.c 使其直接印出 buf
字串:
此實作最大可以計算 。如果我們可以根據 ,預先計算 的十進位至少有多少位數,進而動態調整 BIGNSIZE
,就可以計算任意 。
截至目前為止的實作在 d0dd8a4。
struct BigN
實作我前面提到 struct BigN
中每個 unsigned long long
達到 就會強迫進位,但是接下來實作 Fast Doubling 時,直接對其左移運算會算出錯誤的值,因為 upper
與 lower
並不是差距 ,所以把 lower
的位元直接左移到 upper
完全是錯誤的。
經過思考後,我回復一開始的想法,在相加以及左移時,使用 u128
表示正整數 ,計算完 後要基數轉換時,再讓 lower
依據 進位、化為字串。此作法在進行加法運算時省去了判斷 lower
是否大於 並進位,所以比原版實作更有效率。
decimal | upper | lower |
---|---|---|
0 | 0 | 0 |
1 | 0 | 1 |
2 | 0 | 2 |
- 1 | 0 | 0xff…ff |
1 | 0 | |
1 | 1 | |
1 | 0xff…ff | |
2 | 0 |
以下為改進後的加法 (ubig_add
) 的實作:
觀摩 arthur-chang 的作法,在轉換為字串的過程,先將字元陣列初始化為 "0"
,每讀取一個 bit 就讓自己加自己(也就是乘以二),如果當前是 set bit 就加 1,然後進到下一個迭代,計算過程中都是透過字元陣列儲存十進位的資訊。以下舉 1010
為例
step | num | computation | result |
---|---|---|---|
0 | 1010 | "0" |
"0" |
1 | 1010 | ("0" + "0") + "1" | "1" |
2 | 1010 | ("1" + "1") + "0" | "2" |
3 | 1010 | ("2" + "2") + "1" | "5" |
4 | 1010 | ("5" + "5") + "0" | "10" |
由上述演算法就可以用簡單的機制將 ubig
基數轉換,arthur-chang 的實作是 uint32_t
的版本,我實作 unsigned long long
的版本:
接下來再套用前述的 function 實作 fibdrv:
到目前為止的程式碼 b2749a6。
首先實作 struct BigN
的減法(ubig_sub
)、左移(ubig_lshift
)、直式乘法(ubig_mul
):
接下來參考作業說明提供的演算法:
因此可得:
pseudo code:
舉例: 求解 : 1010 = 10102
i | start | 4 | 3 | 2 | 1 | result |
---|---|---|---|---|---|---|
n | - | 1010 | 1010 | 1010 | 1010 | - |
F(m) | F(0) | F(0*2+1) | F(1*2) | F(2*2+1) | F(5*2) | F(10) |
a | 0 | 1 | 1 | 5 | 55 | 55 |
b | 1 | 1 | 2 | 8 | 89 | - |
用 fast doubling 實作 fibdrv:
到目前為止的程式碼 473ff4a。
首先修改 grub 設定檔,強制 isolate cpu 11:
編輯 GRUB_CMDLINE_LINUX=""
:
更新 grub
重新開機後,用 taskset
檢查 cpu 11 是否已經被 isolate:
我發現有少數 process 仍然有全部核心的 affinity,例如 rcu_gp
。
將 client
程式固定在 cpu 11 執行:
首先利用 ktime_t
與 ktime_get()
來處理計算費氏數列耗費時間,最後透過 fib_read
的回傳值將耗費的時間傳給 user space,具體實作如下:
然後新增一個 perf-test.c
用來印出 user space、kernel space、kernel to user 的執行時間。此程式會重複計算 至 100 次,並取 10% trimmed mean。
接下來輸出測試的結果到 fast_doubling.txt
:
輸出的結果如下:
gnuplot script:
至此可以畫出 Fast Doubling 的 10% trimmed mean 執行時間:
接下來用同樣的手法畫出 Adding 的 10% trimmed mean 執行時間:
比較 Fast Doubling 和 Adding 的 user space time。
可以發現 Fast Doubling 的效能較 Adding 好。
到目前為止的變更 48956a1。
動態配置記憶體需要注意以下議題:
我將大數加法、乘法、初始化等函式統一定義為:
並將不同版本的實作寫在 ubig_1.h
、ubig_2.h
、…,在 fibdrv.c
藉由切換 #include "ubig_xxx.h"
來切換為不同實作。
根據 wikipedia,費波納西數列的前後項比值會趨近於黃金比例,也就是 。
並且, 也就是說當數列的前後項比值接近 1.61804 時,在十進位表示時,每增加四到五項必定會進一位。
列舉費波納西數列,觀察在第幾項之後 會小於 :
13 | 233 | 1.6180555556 |
14 | 377 | 1.6180257511 |
15 | 610 | 1.6180371353 |
16 | 987 | 1.6180327869 |
17 | 1597 | 1.6180344478 |
從列舉看到當 時, 會收斂到小於 1.61804,所以當 時可以藉由 估算第 的十進位位數。此式可被化簡為 。
假設 ubig
使用了 個 unsigned long long
,計算 得到 ubig
可以表示的最大位數,將常數計算完後化簡為 。
只要使 即可根據 的值動態調整 的值。以下便是估計所需 unsigned long long
數量的函式(因為 kernel module 不允許使用浮點數,所以將預先計算的常數同乘以 100000 後,用整數運算):
接下來改寫 ubig
結構,使其可以根據需求調整 unsigned long long
陣列長度。
利用 kmalloc
和 kfree
分配記憶體給 ubig
:
在開始計算 之前,先分配需要動態配置的物件,如果配置失敗就馬上返回。所以 fast doubling 的實作變為:
在第 11 到 19 行先 kmalloc
配置所有 ubig
,在第 20 行確認所需的 ubig
是否全都配置完畢,若否就將所有 ubig
釋放然後回傳 NULL
。接下來就只會使用這些配置好的 ubig
運算,不再呼叫 kmalloc
:
使用 kmalloc
動態配置記憶體後,照理來說最大可以計算 kmalloc
配置的長度的極限(128KB),即 2048 個 unsigned long long
。可以表達的最大的十進位位數為 。
十進位位數小於等於 39457 的費氏數列項次為 ,因此目前的實作最大可以計算 。
到目前為止的實作在 0ee4f89 的 ubig_2.h
。
copy_to_user
回傳 ubig
的原始資料,在資料回傳至 user space 後,再呼叫 ubig_to_string
基數轉換,藉此減少 kernel to user 複製資料的成本。為了達成上述兩項要求,我首先在 fib_sequence
函式增加一個參數 user_size
來接收 user buffer 的大小,在透過 estimate_size
計算完需要多少 unsigned long long
來儲存數值後,再檢查 buffer 長度夠不夠儲存這麼長的 unsigned long long
陣列:
計算完 之後,就讓 copy_to_user
把 ubig
的 unsigned long long
陣列複製進 user buffer,然後回傳 unsigned long long
陣列長度(每 8 byte 一個單位):
到目前為止的實作在 8e0021f。
近日面試被問到假設 fibdrv 在 userspace 實現,會與 kernel module 有什麼差異,參照 Differences Between Kernel Modules and User Programs,其中提到 Kernel modules have higher execution privilege,故 kernel module 本身被分配到 cpu cycle 執行的機會就比 userspace application 高,所以比較快執行完畢。
以下是我將 fast doubling 做成 user space 的 function 來與 fibdrv 進行執行時間的比較(只計算 ,不包含基數轉換):
可以從圖中看出 kernel module 確實比 userspace application 快。
接下來測試計算費氏數列加上基數轉換所耗費的時間,發現基數轉換的效率非常差,是接下還需要改善的部份。
TODO: 引入 Schönhage–Strassen algorithm 來加速乘法運算,注意 threshold 的設定
bignum
jserv
Schönhage–Strassen 的概念是將 、 以 個位數為單位,將大數拆成若干個較小的數 、…、、、…,對 和 線性卷積,執行完線性卷積後將所有元素依據自己的位數左移並相加就可以完成乘法運算。舉例:
令
將 、 以 2 個位數為界分割會得到以下兩個數列:
對數列線性卷積:
為 與 的線性卷積,接下來透過左移和加法將線性卷積的結果相加即可完成大數乘法:
在 ubig
的乘法實作中,我打算以 32 bit 為界切割數值,並透過 64 bit 的乘法來計算線性卷積。為了方便上述運算,我將 ubig
儲存數值的資料型態由 unsigned long long
改為 unsigned int
。
我稍微調整乘法計算方式,按照下圖 col1
到 col5
順序計算卷積、處理進位,然後直接存到 dest
中:
完整實作在 b21a89a 的 ubig_schonhange_strassen.h
。
Karatsuba 的概念是將 、 以第 位數為界,拆成兩半 、、、,把這他們視為較小的數相乘,然後再透過左移補回 、 損失的位數,以十進位為例:
則 可以化為:
上述算法計算 、、 需要四次乘法,我們還可以透過以下技巧優化成三次乘法:
觀察 展開的結果
移項之後,我們就能利用 、、 來計算
最後計算 便能得到 相乘的結果。
以上便是 Karatsuba 演算法。
接下來我以兩個 8 bit 數值相乘變成 16 bit 來演示 Karatsuba(假 設所採用的處理器只支援 8 bit 乘法):
令
其中 可以用左移運算代替。
接續上一個例子,當 、 超過 8 bit 時,可以透過分治法實作 Karatsuba。、、、 的位元數超出處理器的乘法的位數時,就把他們再切為 、、、、…,再使用 Karatsuba 計算。以下以兩個 16 bit 數值相乘變成 32 bit 來演示(假設所採用的處理器只支援 8 bit 乘法):
由上圖可以看出計算 、、 時,透過分治法將 、、、 切成更小的數字執行乘法運算。最後再用左移與加法計算 即可求得結果。
至此我可以透過分治法使用 Karatsuba 計算任意位數的大數。
我目前的實作並未考慮 ubig
可能需要 resize,所以需要大幅修改 ubig
的賦值、加法、左移等基本操作,才好實作 Karatsuba。
另一種解決辦法是將所有 unsigned int
陣列大小都統一為 128 KiB,這樣就不用在分治的時候處理不同長度的 unsigned int
陣列,但是計算較小的數會浪費記憶體。
此外,計算 、、 需要暫存的 ubig
物件,所以此實作是 not-in-place。以下為基於遞迴的 Karatsuba 的 ubig_mul
實作:
完整實作在 9673a70 的 ubig_karatsuba.h
。
下圖是計算 、、、…、 的執行時間(不包含基數轉換)。用加法計算 需要約 200 毫秒,用長乘法則需要更長時間,因不明原因,在測量長乘法計算 以後的時間失準。
單獨看 Karatsuba 與 Schönhage–Strassen,前者計算 時需要消耗 45 毫秒;後者則需要消耗 1.8 毫秒左右。
接下來看計算較小數字時各演算法消耗的時間:
從圖中可以發現,計算 到 時用 Adding 最快、計算 到 用 Karatsuba 與 Schönhage–Strassen 差不多快、計算之後的數字則是 Schönhage–Strassen 最快。
針對 計算三種演算法的時間複雜度:
unsigned int
陣列的長度次迭代,也就是 次,所以複雜度為 。