# 2023 q1 Homework3 (fibdriv)
contributed by<[WeiHongWi](https://github.com/WeiHongWi/fibdrv)>
## 開發環境
```shell
gcc --version
gcc (Ubuntu 11.3.0-1ubuntu1~22.04) 11.3.0
Copyright (C) 2021 Free Software Foundation, Inc.
Architecture: x86_64
CPU op-mode(s): 32-bit, 64-bit
Address sizes: 39 bits physical, 48 bits virtual
Byte Order: Little Endian
CPU(s): 12
On-line CPU(s) list: 0-11
Vendor ID: GenuineIntel
Model name: 11th Gen Intel(R) Core(TM) i5-11400H @ 2.70GHz
CPU family: 6
Model: 141
Thread(s) per core: 2
Core(s) per socket: 6
Socket(s): 1
Stepping: 1
CPU max MHz: 4500.0000
CPU min MHz: 800.0000
BogoMIPS: 5376.00
```
- [ ] 研讀上述 ==Linux 效能分析的提示== 描述，在自己的實體電腦運作 GNU/Linux，做好必要的設定和準備工作 $\to$ 從中也該理解為何不希望在虛擬機器中進行實驗;
- [ ] 研讀上述費氏數列相關材料 (包含論文)，摘錄關鍵手法，並思考 [clz / ctz](https://en.wikipedia.org/wiki/Find_first_set) 一類的指令對 Fibonacci 數運算的幫助。請列出關鍵程式碼並解說
- [ ] 複習 C 語言 [數值系統](https://hackmd.io/@sysprog/c-numerics) 和 [bitwise operation](https://hackmd.io/@sysprog/c-bitwise)，思考 Fibonacci 數快速計算演算法的實作中如何減少乘法運算的成本;
- [ ] 學習指出針對大數運算的加速運算和縮減記憶體操作成本的舉措，紀錄你的認知和疑惑
- [ ] 注意到 `fibdrv.c` 存在著 `DEFINE_MUTEX`, `mutex_trylock`, `mutex_init`, `mutex_unlock`, `mutex_destroy` 等字樣，什麼場景中會需要呢？撰寫多執行緒的 userspace 程式來測試，觀察 Linux 核心模組若沒用到 mutex，到底會發生什麼問題。嘗試撰寫使用 [POSIX Thread](https://en.wikipedia.org/wiki/POSIX_Threads) 的程式碼來確認。 $\to$ 搭配閱讀〈[並行和多執行緒程式設計](https://hackmd.io/@sysprog/concurrency)〉
---
## 費氏數列
### Fast Doubling
\begin{gather*}
F(2k) = F(k)[2F(k+1) - F(k)]
\end{gather*}
\begin{gather*}
F(2k+1) = F(k+1)^2 + F(k)^2
\end{gather*}
Fast Doubling 的遞迴
```graphviz
strict digraph G
{
1[label="F(6)"]
2[label="F(3)"]
3[label="F(4)"]
4[label="F(1)", style=filled]
5[label="F(2)", style=filled]
6[label="F(2)", style=filled]
7[label="F(3)"]
8[label="F(1)", style=filled]
9[label="F(2)", style=filled]
1 -> {2, 3}
2 -> {4, 5}
3 -> {6, 7}
7 -> {8, 9}
}
```
可以發現有重複計算到的fibonacci sequence , 若能夠避免重複運算 , 則遞迴能夠表示成下圖:
```graphviz
strict digraph G
{
1[label="F(6)"]
2[label="F(3)"]
3[label="F(4)"]
4[label="F(1)", style=filled]
5[label="F(2)", style=filled]
6[label=" " ,color=white]
7[label=" " ,color=white]
{rank = same; 2;3;}
{rank = same; 4;5;}
1 -> {2, 3}
2 -> {4, 5}
2 -> 3 [style=dashed; arrowhead=vee]
5 -> 3 [style=dashed; arrowhead=vee]
3 -> {6, 7} [color=white]
}
```
教材所提到的硬體加速的 Fast Doubling , 首先利用 **builtin_clz(),算出幾個 leading zero bits** , 省略這些用不到的 bits ,且遞迴從 F(1) 開始算:
(1) 遇到 " 0 " : 則計算 F(2k), F(2k+1)
(2) 遇到 " 1 " : 則先計算 F(2k),F(2k+1) ,再計算 F(2k+2)
### Fast doubling 實作
依據教材給的 pseudo code , 利用 __builtin_clz() 算出 number of leading zero ,然後 current_digit 則用左移來做成一個 mask 用來判斷當下的 bit 是否為 1 , 判斷式用 k & current_digit.
```clike=
static long long fib_sequence(long long k)
{
long long a = 0, b = 1;
int lz_bit = __builtin_clzll(k);
long long current_digit = 1 << (64 - lz_bit - 1);
for (int i = 64 - lz_bit; i >= 1; i--) {
long long t1 = a * (2 * b - a);
long long t2 = pow1(a, 2) + pow1(b, 2);
a = t1;
b = t2;
if (k & current_digit) {
t1 = a + b;
a = b;
b = t1;
}
current_digit >>= 1;
}
return a;
}
```
## 改變測試條件的最大值
### 解讀 Makefile + Verify.py + Expected.txt
```=
check: all
$(MAKE) unload
$(MAKE) load
sudo ./client > out
$(MAKE) unload
@diff -u out scripts/expected.txt && $(call pass)
@scripts/verify.py
```
其中 line 6 將輸出檔 out file 和在 scripts/expected.txt 做比較如果正確了會出現下面的唯一一句：
```
Passed [-]
```
而接後出現的都為上面 Makefile 段落的 line 7 , 也就是 verify.py 所產生出的結果
```
f(93) fail
input: 7540113804746346429
expected: 12200160415121876738
```
除了更改 client.c 的 offset 變成 500 還要去改變 expected.txt 這個檔案,因為原始的**expected.txt檔的內容是錯誤答案**,其中的內容為 Fibonacci Sequence F(n),其中n = 1~100 且 overflow 發生在 $F(93)$ 的情況
所以需要讓 expected.txt 擁有 F(n),且 n = 1,2,...,500的內容且使得之後可以用 make check 對上正確答案且可以對到 $F(500)$
否則Makefile 會出現
:::danger
make: *** [Makefile:40：check] 錯誤 1
:::
> [name=HongWei][time=Tue,Mar 7, 2023 00:36 AM]
### out
發現改變 fibdriv.c 中的內容之後再去 make check 之後內容會隨之改變, 所以確定為輸出檔案
---
## 改寫 fib_sequence()
### Overflow 的原因
```clike=
long long f[k + 2]
```
已知 **type long long** 可以表達的範圍為$- 2^{64}$ ~ $2^{64}-1$,則 $F(92)$ = 7540113804746346429 約等於 $log$(7.5*$10^{18}$) = $log(7.5)$ + $log(10^{18})$ 約等於 63 , 所以當 F(93) 時發生 Overflow
> [name=HongWei][time=Tue,Mar 7, 2023 01:16 AM]
### Big number 運算 - 數字字串
```clike=
typedef struct Bignum {
char element[128];
} bn;
```
* 因為用 long long int 也會導致 overflow , 所以我們需要用 string 來表示 , 每一個 element 都代表一個位數
> [name=HongWei][time=Tue,Mar 7, 2023 10:16 PM]
* 根據教材,實作 string_add_num,可以發現若我們能夠固定 size 的大小次序 ,則不需要去比較 size_a 和 size_b 的大小,舉例來說, f[n-1]的 strlen() 一定大於等於 f[n-2] 的 strlen().
> [name=HongWei][time=Wed,Mar 8, 2023 4:46 PM]
```clike=
static void string_number_add(char *a, char *b, char *out)
{
reverse_str(a);
reverse_str(b);
int i;
int carry = 0, sum = 0;
size_t size_a = strlen(a);
size_t size_b = strlen(b);
for (i = 0; i < size_b; i++) {
sum = (a[i] - '0') + (b[i] - '0') + carry;
out[i] = '0' + sum % 10;
carry = sum / 10;
}
for (i = size_b; i < size_a; i++) {
sum = (a[i] - '0') + carry;
out[i] = '0' + sum % 10;
carry = sum / 10;
}
if (carry)
out[i++] = '0' + carry;
out[i] = '\0';
reverse_str(a);
reverse_str(b);
reverse_str(out);
}
```
* 最後修改的 fib_sequence() , 其中**copy_to_user**()將想傳給使用者模式 (即運作中的 client) 的字串複製到到 fib_read 的 buf 參數後，client 端方可接收到此字串
```clike=
static long long fib_sequence_str(long long k, char *buf)
{
bn *f = kmalloc(sizeof(bn) * (k + 2), GFP_KERNEL);
strncpy(f[0].element, "0", 1);
f[0].element[1] = '\0';
strncpy(f[1].element, "1", 1);
f[1].element[1] = '\0';
for (int i = 2; i <= k; ++i) {
string_number_add(f[i - 1].element, f[i - 2].element, f[i].element);
}
long long ret_size = strlen(f[k].element) + 1;
unsigned long test = __copy_to_user(buf, f[k].element, ret_size);
if (test) {
printk("The copy from kernel to user is fail.");
return -1;
}
return ret_size;
}
```
:::info
依舊無法正確表示 F(93) 以上,還在釐清原因,不知道是否 client.c 需要修改
> [name=HongWei][time=Tue,Mar 7, 2023 5:53 PM]
當我在 client.c 把 ret_size 給 printf() 之後發現長度卡在 20 就停住了
```
Reading from /dev/fibonacci at offset 93, returned the sequence 7540113804746346429.
20
Reading from /dev/fibonacci at offset 92, returned the sequence 7540113804746346429.
20
```
> 已解決
因為在 loff_t fib_device_lseek() 用的是 MAX_LENGTH 所以需要修改它 ,才能讓 file 讀到下一個位置！！！
> [name=HongWei][time=Tue,Mar 7, 2023 7:13 PM]
:::
結果如下：
```
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.
```
### 時間測量與效能分析
#### 利用 ktime 改寫 fb_write 和 fb_read
利用教材的 fib_time_proxy() , 其中 kt 為 ktime_t 的 object
```c=
static long long fib_time_proxy(long long k, char *buf)
{
kt = ktime_get();
long long result = fib_sequence_str(k, buf);
kt = ktime_sub(ktime_get(), kt);
return result;
}
```
並將 fib_read 和 fib_write() 改寫成
```c=
static ssize_t fib_read(struct file *file,
char *buf,
size_t size,
loff_t *offset)
{
return (ssize_t) fib_time_proxy(*offset, buf);
}
```
```c=
static ssize_t fib_write(struct file *file,
const char *buf,
size_t size,
loff_t *offset)
{
return ktime_to_ns(kt);
}
```
但發現 fb_write() 的結果依舊是:
```
-Writing to /dev/fibonacci, returned the sequence 0
-Writing to /dev/fibonacci, returned the sequence 0
-Writing to /dev/fibonacci, returned the sequence 0
-Writing to /dev/fibonacci, returned the sequence 0
-Writing to /dev/fibonacci, returned the sequence 0
-Writing to /dev/fibonacci, returned the sequence 0
-Writing to /dev/fibonacci, returned the sequence 0
-Writing to /dev/fibonacci, returned the sequence 0
.
.
.
```
也就是寫成
```c=
sz = read(fd, buf, 1);
long long ktime = write(fd, write_buf, strlen(write_buf));
printf("Writing to " FIB_DEV ", returned the sequence %lld\n", ktime);
```
:::info
要改寫 client.c 本身 , 且教材本身的 function 為 fib_time_proxy() 計算完後 kt 的值才會有所變化,此時再去 call fib_write() 才是正確的.
:::
#### copy_to_user()
根據教材 , 可以得知 copy_to_user() 將 kernel 端的字串傳回 user 端 fib_read 中的 argument buf 所需要用的函數,所以**分析時間效能的時候,也要將複製資料的時間考慮進去**,我的做法是另外宣告一個 static ktime_t kt1 來紀錄在**fib_sequence_str**()中的copy_to_user() 的時間,如以下：
```c=
kt1 = ktime_get();
unsigned long test = __copy_to_user(buf, f[k].element, ret_size);
kt1 = ktime_sub(ktime_get(), kt1);
```
而後再把 client.c 改寫 , 將 kernel_to_user , kernel time , user time 都寫進 file data.txt.
```c=
long long utime = write(fd, write_buf, strlen(write_buf));
printf("Writing to " FIB_DEV ", returned the sequence %lld\n", utime);
fprintf(data, "%lld %lld %lld\n", utime, kernel_to_user, utime - kernel_to_user);
```
主要是想減少當要複製的資料變大時,所造成的時間差,但我認為我的作法出錯了,目前結果如下,其中 col 2 為上述計算 __copy_to_user() 的時間,我發現值都停在40幾附近:
```
c0 c1 c2 c3
92 6509 44 6465
93 6772 40 6732
94 7148 41 7107
95 7255 41 7214
96 7025 42 6983
97 7108 41 7067
98 8276 40 8236
99 7237 41 7196
100 7387 41 7346
```
:::info
結果是 kt1 紀錄到的是 kernel space 呼叫 function copy_to_user()所花的時間,所以值才會固定.參閱教材時有發現 user time 都高於 kernel time ,但實作時有個**錯誤的認知** **copy_to_user() 在 kernel space 複製資料到 user space ,所以直覺地認為copy_to_user() 花的也是 kernel time.**
:::
> [name=HongWei][time=Mon,Mar 13, 2023 6:17 AM]
#### User,Kernel,Kernel-to-User
首先我們需要可以在 user (也就是 client.c)量測到 ns 單位的 timer.利用作業要求[]()的 clock_gettime() 來實作:
```c=
static inline long long u_ntime()
{
struct timespec ts;
clock_gettime(CLOCK_REALTIME, &ts);
return ts.tv_sec * 1e9 + ts.tv_nsec;
}
```
其中 structure timespec 有兩個成員,分別是 tv_sec 為 second(s),tv_nsec為 nano second(ns).
而 [clock_get_time()](https://linux.die.net/man/2/clock_gettime)中的 CLOCK_REALTIME 是整個系統通用的時鐘,用以量測現實生活的時間,與 clockid CLOCK_MONOTONIC 的差異是, CLOCK_REALTIME 的設定是可以變動的.
```c=
long long start = u_ntime();
sz = read(fd, buf, 1);
long long utime = u_ntime() - start;
long long ktime = write(fd, write_buf, strlen(write_buf));
printf("Writing to " FIB_DEV ", returned the sequence %lld\n", utime);
fprintf(data, "%d %lld %lld %lld\n",i, utime,ktime, utime - ktime );
```
根據教材,將干擾因素都排除後,實際上的差異為
![]()
#### 更改 fib_read()
原本 bn_to_string() 都放在fib_sequence 系列的 function 裡面,但發現這樣的寫法並不美觀,且會增加 fib_sequence 系列 function 的執行時間,原本的 code 如下
```clike=
...
for (unsigned int i = 1; i < k; i++) {
bn_swap(b, dest);
bn_add(a, b, dest);
bn_swap(a, b);
}
bn_free(a);
bn_free(b);
char *str = bn_to_string(dest);
unsigned long sz = strlen(str);
unsigned long long test = __copy_to_user(buf, str, sz)
...
```
也就是將 line 10 到 line 11 都移至 fib_read() 裡 , 這樣也就不需要教材多寫的 fib_proxy_time() 了
```clike=
...
bool doubling = false;
bn *res = bn_alloc(1);
if (doubling) {
kt = ktime_get();
fib_sequence_fastd(*offset, res);
kt = ktime_sub(ktime_get(), kt);
}
else {
kt = ktime_get();
fib_sequence_bn(*offset, res);
kt = ktime_sub(ktime_get(), kt);
}
char *p = bn_to_string(res);
...
```
> [name=HongWei][time=Fri,Mar 17, 2023 12:47 AM]
### 加速 Fibonacci 運算
教材採用的 bn 架構
```clike=
typedef struct _bn {
unsigned int *number;
unsigned int size;
int sign;
} bn;
```
#### bn_to_string()
因為大數沒辦法以數字的形式表示,所以需要轉成字串表達,以下分段解讀 bn_to_string().
```c
size_t len = (8 * sizeof(int) * src.size) / 3 + 2 + src.sign;
```
(8*sizeof(int)*src.size)用來表達這個 big number 共需要幾個 bits 來表示,再來除以3 ,來表示$log_{10}X$ , 也就是以 10 為底 , 需要幾個 number 去做表示(個,十,百,千 ....) ,如果 big number 為負,則在多加1.
```clike=
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;
}
}
}
```
第一個 for loop 用在尋訪在 big number 裡的全部 integer , 第二個 for loop 則是用來尋訪一個 integer 的所有 bits , 並用 s[j] += s[j] - '0' + carry 來將其轉成字元 '0'~'9'
這邊我一開始想不通的點是,當 src.number[i] 超過 255 時就需要用 9 個 bits 以上來表達,則上述式子因為 bits 長度不符合,導致結果有所差異,舉例來說src.number[i] = 256 and binary representation 100000000~2~ , then **carry will seen as 0 in the statement s[j] += s[j] - '0' + carry** because of the length of the object.
目前能想到的解釋是,因為他是 character array , 所以為**連續的記憶體** , 也就是說它依舊還是可以儲存到 , 我認為這就是第三個 for loop 從尾端 s[len-2] 開始存放個位數 , 因為這樣才能利用 allocated memory 來存放超出8個 bits 的值
#### bn_do_add()
```clike=
int d = MAX(bn_msb(a), bn_msb(b)) + 1;
d = DIV_ROUNDUP(d, 32) + !d;
```
其中 DIV_ROUNDUP(d,size) = (d+size-1)/size , 也就是去計算共有幾個 element 在這個 big number 裡 , 但我**仍未釐清 !d 的用途**
```clike=
bn_resize(bn *src, size_t size);
```
看 [fibdrv/bn_kernel.c](https://github.com/KYG-yaya573142/fibdrv/blob/master/bn_kernel.c#L22)後知道 function 在分配空間前,要確認 src 是否為 NULL , size 是否等於 src->size , 和最後確認輸入的 argument size 是否為0 , 而最後的 if statement 用來**分配記憶體大小 size - src->size 給予新的 object**.以下為程式碼
```c=
if (size > src->size)
memset(src->number + src->size, 0, sizeof(int) * (size - src->size));
```
而後教材所給予的程式碼中 bn_mult() 佔據太多時間的問題
```
79.32% 0.00% client [kernel.kallsyms] [k] entry_SYSCALL_64_after_hwframe
|
---entry_SYSCALL_64_after_hwframe
do_syscall_64
|
|--49.30%--__x64_sys_read
| ksys_read
| vfs_read
| fib_read
| bn_mult
| |
| |--27.22%--kfree
| |
```
:::success
之前的作法把 bn_to_string 也放入 fib_sequence 系列的 function 裡頭 ,這樣去算 ktime_t 相關的時間時,無法真正考量到 fib_sequence 本身的運算,所以改善方案1和2的圖只是參考用,之後再去做修改,且會以 $F(1000)$ 做基準,目前的做完改善方案1和2之後的時間圖如下:
:::
> [name=HongWei][time=Fri,Mar 17, 2023 12:47 AM]
#### 改善方案1 : bn_cpy(),bn_mult()的效能改善
和教材的寫法一樣,其中 f1 : $F(k)$ , f2 : $F(k+1)$
```c=
...
for (unsigned int i = 1U << 31; i; i >>= 1) {
bn_cpy(k1, f2);
bn_lshift(k1, 1);
bn_sub(k1, f1, k1);
bn_mult(k1, f1, k1);
bn_mult(f1, f1, f1);
bn_mult(f2, f2, f2);
bn_cpy(k2, f1);
bn_add(k2, f2, k2);
if (k & i) {
bn_cpy(f1, k2);
bn_cpy(f2, k1);
bn_add(f2, k2, f2);
} else {
bn_cpy(f1, k1);
bn_cpy(f2, k2);
}
}
```
line 3 到 line 6 做的是\begin{gather*}
F(2k) = F(k)[2F(k+1) - F(k)]
\end{gather*}
而 line 8 到 line 11 做的是
\begin{gather*}
F(2k+1) = F(k+1)^2 + F(k)^2
\end{gather*}
而教材提到的為用 bn_swap() 來**避免掉 bn_cpy() 中複製資料的成本**, bn_cpy() 如下
```c=
int bn_cpy(bn *dest, bn *src)
{
...
dest->sign = src->sign;
memcpy(dest->number, src->number, src->size * sizeof(int));
...
}
```
首先利用 bn_add() 將 f2 = $F(k+1)$ 加一倍後放置 k1,以避免第一次的 bn_cpy()
經過$F(2k)=F(k)[2F(k+1)-F(k)]$ 的運算後將其存放在 k1 , 因為 f1 還需要去做 $F(2k+1)$ 的部份,之後再利用 bn_swap(k1,f1) 交換
```clike=
for (unsigned int i = 1U << (31-__builtin_clz(k)); i; i >>= 1) {
/* F(2k) = F(k) * [ 2 * F(k+1) – F(k) ] */
bn_add(f2,f2,k1);
bn_sub(k1, f1, k1);
bn_mult(k1, f1, k1);
//F(2k+1) = F(k)^2 + F(k+1)^2
bn_mult(f1, f1, f1);
bn_mult(f2, f2, f2);
bn_add(f1, f2, f2);
bn_swap(k1,f1);
if (k & i) {
bn_swap(f1,f2);
bn_add(f1, f2, f2);
}
}
```
和原本(v0)做比較的效能圖：
但 bn_mult(f1,f1,f1) 和 bn_mult(f2,f2,f2) 會因為 bn_mult() 中的 source 和 destination 相同而導致會多用記憶體去存放 c .
```clike=
if (c == a || c == b) {
tmp = c; // save c
c = bn_alloc(d);
}
```
:::info
根據教材的改法,值目前是錯的,還要在思考一下.
> [name=HongWei][time=Mon,Mar 15, 2023 6:17 PM]
> 已修正
:::
#### 改善方案2 : Q-matrix
\begin{gather*}
F(2k) = F(k)[2F(k-1) + F(k)]
\end{gather*}
\begin{gather*}
F(2k-1) = F(k-1)^2 + F(k)^2
\end{gather*}
```c=
for (unsigned int i = 1U << (30 - __builtin_clz(n)); i; i >>= 1)
```
#### 改善方式 4 : uint64
bn 的 struct 變為:
```c=
...
typedef uint64_t bn_data
...
typedef struct _bn {
bn_data *number;
bn_data size;
bn_data capacity;
int sign;
} bn;
```
此時遇到的第一個問題為
:::info
warning: left shift count >= width of type
:::
其中問題出在 for(bn_data_tmp d = 1U <<64),根據[stackoverflow](https://stackoverflow.com/questions/4201301/warning-left-shift-count-width-of-type),得知在右式中的1也需要變成 unsigned long long 的資料型態
:::info
在修改為 64 bits 的 struct 之後 , 只要跑 make check 無法正常執行的情況 , 但轉為 32 bits 的 struct 就又正常,目前還在思考哪邊出了問題.
:::
> [name=HongWei][time=Sun,Mar 19, 2023 13:47 AM]
#### 改善方式 5 : 改寫 bn_add() 和 bn_mult()
##### 改寫 bn_do_add()
```c=
...
int d = MAX(bn_msb(a), bn_msb(b));
...
for (int i = 0; i < a->size; ++i) {
unsigned int tmp1 = a->number[i];
unsigned int tmp2 = b->number[i];
carry = (tmp2 += carry) < carry;
carry += (c->number[i] = tmp1 + tmp2) < tmp1;
}
if (b->size > a->size) {
for (int i = a->size; i < b->size; ++i) {
unsigned int tmp = b->number[i];
carry = (tmp += carry) < carry;
c->number[i] = tmp;
}
}
```
照著教材改寫 , 可以發現在 line 7 中 b->size > a->size ,寫這樣是因為我在用到 bn_add() 的場合都有特別去注意 b 要比 a 要大,但這樣**要將 function 移用到其他 function 時,還要多去注意這個條件**,容易造成失誤.
##### 改寫 bn_mult()
和教材中的程式碼相同 , 為 32 bits 的版本
```c=
static unsigned int bn_mult_partial(const bn *a,
unsigned int sz,
unsigned int b,
unsigned int *c)
{
if (b == 0) {
return 0;
}
unsigned long long carry = 0;
for (int i = 0; i < a->size; ++i) {
unsigned int high, low;
unsigned long long product = (unsigned long long) a->number[i] * b;
low = product;
high = product >> 32;
carry = high + ((low += carry) < carry);
carry += ((c[i] += low) < low);
}
return carry;
}
```
> [name=HongWei][time=Sun,Mar 19, 2023 7:01 PM]
#### 改善方案 6 : bn_sqrt()
當重新回顧 fib_sequence_fastd() 時,可以發現使用所有的 bn_mult() 中 ,有兩個 bn_mult() 是 source 的平方 , 也就是在實現 \begin{gather*}
F(2k-1) = F(k-1)^2 + F(k)^2
\end{gather*}
所以根據教材
```
a b c
x a b c
---------------
ac bc cc
ab bb bc
aa ab ac
```
可以發現 bc ac ab 只要個別算一次 , 之後在相同位置乘 2 , 就可以節省再次相乘的成本. 實作如下 , 重點擺在 **bn_mult_partial(&ap[i+1], asize+i, ap[i], cp)** , ap[i] 為乘數,而 &ap[i+1] 大於 ap[i] 一個位置的地方開始做,也就是不去計算平方項 $a^2,b^2,...$ .
```c=
...
unsigned int *cp = c+1;
unsigned int *ap = a;
for (int i = 0; i < size-1; ++i) {
cp[size-1- i] = bn_mult_partial(&ap[i+1], asize+i, ap[i], cp);
cp += 2;
}
/* Double it */
for (int i = 2 * size - 1; i > 0; i--)
c[i] = c[i] << 1 | c[i - 1] >> (31);
c[0] <<= 1;
```
之後再把對角線的 $a^2,b^2,c^2$ 加上 , 其中 cp += 2 為對角線的下一個位置
```c=
for (int i = 0; i < asize; i++) {
bn_data high, low;
unsigned long long product = (unsigned long long)a[i] * a[i];
low = product;
high = product >> 32;
high += (low += carry) < carry;
high += (cp[0] += low) < low;
carry = (cp[1] += high) < high;
cp += 2;
}
```
> [name=HongWei][time=Sun,Mar 19, 2023 11:01 PM]
#### 改善方案 7 : Karatsuba Algorithm
:::info
有看以前的同學所做的作品,發現大家都做在 user mode , 這裡留有一個疑問 , 好奇是否在kernel mode 有限制才在 user mode做.
:::
閱讀程式碼[sysprog21/bignum/mul.c](https://github.com/sysprog21/bignum/blob/master/mul.c) , 以下為當割成一半的 number 依舊大於手動設置的 KARATSUBA_MUL_THRESHOLD 時則再切割一次,所以才放入遞迴.
目標,則是去計算 $z0 , z1 , z2$,其中 $z1 = (x1+x0)*(y1+y0) - x0*y0 - x1*y1$
(1)
一開始先把 data , 切割成 4 等塊 , 所以我們需要去計算 **sz 的一半**
```c=
unsigned int even_sz = (sz >> 1) << 1;
unsigned int half_pos = even_sz / 2
```
(2)
將 z0 和 z2 利用 bn_mult_sup() 計算, 利用pointer 讓 z0 , z2 在正確的位置以此省略左移的操作,舉例來說 z2 的位置從 c + even_sz 開始 , 並一直存放到 c + 2*even_sz - 1
```c=
unsigned int *cp = c;
unsigned int *cph = c + even_sz;
// z0 = x0y0 , z2 = x1y1
bn_mult_sup(x0, half_pos, y0, half_pos, cp);
bn_mult_sup(x1, half_pos, y1, half_pos, cph);
```
(3)計算 (x1-x0)*(y0-y1) 過長,可以至 [WeiHongWi](https://github.com/WeiHongWi/fibdrv/blob/master/bn.c)觀看.
(4)
若一半的長度還是大於我們設定的值,我們可以再切一半遞迴丟入 bn_karatsuba()
```c=
// Calculate the value |(x1-x0)|*|(y0-y1)|
tmp = (unsigned int *) (calloc(even_sz, sizeof(unsigned int)));
if (half_pos > KARATSUBA_MUL_THRESHOLD) {
bn_karatsuba(a_tmp, b_tmp, half_pos, tmp);
} else {
bn_mult_sup(a_tmp, b_tmp, half_pos, half_pos, tmp);
}
```
(5)
若 sz 屬於 odd number , 我利用 sz ^ 1 來做判斷 , 且另外需要加上多出來number 和 number[sz-1] 做相乘
```c=
if (sz ^ 1) {
c[even_sz * 2] = bn_mult_partial();
c[even_sz * 2 + 1] = bn_mult_partial();
}
```
效能圖待補 ...
Sign in with Wallet
Connect another wallet