---
tags: linux2023
---
# 2023q1 Homework1 (fibdrv)
contributed by < [Jerejere0808](https://github.com/Jerejere0808) >
新增以下資料結構使 fibonnaci 可以進行大數運算，其中 capacity 是用於分配多餘的記憶體空間，這樣就可以避免多次的 resize (realloc) ， 只有在 capacity 不足時才會需要重新分配記憶體空間。
```c
typedef struct _bn {
unsigned int *number;
unsigned int size;
int sign;
int capacity;
} bn;
```
將計算 fibonnaci 的方法改成 fast doubling，使時間複雜度可以達到 log(n)。
:::warning
注意就算使用大數運算仍會有運算數字大小的限制。因為 kmalloc 會有分配空間的限制，分配超過某個 size 會導致錯誤的情況，而分配空間的限制跟每台電腦的系統配置有關，如果假設最大的分配空間為 4MB ， 那麼根據 Binet 公式 ， 當 n 為足夠大的正整數時可以推算需要使用的位元數: − 1.160964 + n × 0.694242
當 n = 5761680 ， f(n) 就會需要 4MB 來表示，所以當 n > 5761680 就會產生錯誤。如果是以 struct _bn 裡的 unsigned int 陣列來看，Bytes 數量就必須除以 32 加 + 1 ，也就是 unsigned int *number 最多可以分配 125000 個 unsigned int。
:::
```c
void bn_fib_fdoubling(bn *dest, unsigned int n)
{
bn_resize(dest, 1);
if (n < 2) { // Fib(0) = 0, Fib(1) = 1
dest->number[0] = n;
return;
}
bn *f1 = dest; /* F(k) */
bn *f2 = bn_alloc(1); /* F(k+1) */
f1->number[0] = 0;
f2->number[0] = 1;
bn *k1 = bn_alloc(1);
bn *k2 = bn_alloc(1);
/* walk through the digit of n */
for (unsigned int i = 1U << 31; i; i >>= 1) {
/* F(2k) = F(k) * [ 2 * F(k+1) – F(k) ] */
bn_cpy(k1, f2);
bn_lshift(k1, 1);
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_cpy(k2, f1);
bn_add(k2, f2, k2);
if (n & i) {
bn_cpy(f1, k2);
bn_cpy(f2, k1);
bn_add(f2, k2, f2);
} else {
bn_cpy(f1, k1);
bn_cpy(f2, k2);
}
}
// return f[0]
bn_free(f2);
bn_free(k1);
bn_free(k2);
}
```
可以用以下方法來測量在 kernel 所需多少計算時間
```c
kt = ktime_get();
bn_fib_fdoubling(result, *offset);
kt = ktime_sub(ktime_get(), kt);
```
參考 [colinyoyo26](https://github.com/colinyoyo26/fibdrv) 的方法用 clock_gettime 在 client 取得 user 所需的時間
```c
static inline long long get_nanotime()
{
struct timespec ts;
clock_gettime(CLOCK_REALTIME, &ts);
return ts.tv_sec * 1e9 + ts.tv_nsec;
}
```
兩者相減可以得到 kernel to user 的時間，如下圖:
不過我發現所需的時間比預期 log(n) 還多出許多，fib(300) 就需要 50000 ns 經過檢查，發現我把 bn_to_string ，也就是轉化成字串的部分算進去 kernel time 裡面了 ，把這部分的計算時間去掉之後，時間有明顯變短，如下圖所示，可見把 bn 轉成字串是非常耗時間的。
再持續修改，原本的 bn_fib_fdoubling 有 bn_cpy 也就是需要複製資料，這樣會需要額外成本，因此透過修改 bn_lshift 使 src 可以直接把左移的結果給 dest 就不用先複製之後又再往左位移。
ex: bn_cpy(k1, f2); bn_lshift(k1, 1); 可以改成 bn_lshift2(f2, 1, k1)
bn_cpy 和 舊的 bn_lshift 如下
```c
int bn_cpy(bn *dest, bn *src)
{
if (bn_resize(dest, src->size) < 0)
return -1;
dest->sign = src->sign;
memcpy(dest->number, src->number, src->size * sizeof(bn_data));
return 0;
}
```
```c
void bn_lshift(bn *src, size_t shift)
{
size_t z = bn_clz(src);
shift %= 64; // only handle shift within 32 bits atm
if (!shift)
return;
if (shift > z)
bn_resize(src, src->size + 1);
for (int i = src->size - 1; i > 0; i--)
src->number[i] =
src->number[i] << shift | src->number[i - 1] >> (64 - shift);
src->number[0] <<= shift;
}
```
新的 bn_lshift 和 fast doubling 實作如下:
```c
void bn_lshift(const bn *src, size_t shift, bn *dest)
{
size_t z = bn_clz(src);
shift %= 32; // only handle shift within 32 bits atm
if (!shift)
return;
if (shift > z) {
bn_resize(dest, src->size + 1);
bn_resize(src, src->size + 1);
} else {
bn_resize(dest, src->size);
}
for (int i = src->size - 1; i > 0; i--)
dest->number[i] =
src->number[i] << shift | src->number[i - 1] >> (32 - shift);
dest->number[0] = src->number[0] << shift;
}
```
```c
void bn_fib_fdoubling_nocpy(bn *dest, unsigned int n)
{
bn_resize(dest, 1);
if (n < 2) { // Fib(0) = 0, Fib(1) = 1
dest->number[0] = n;
return;
}
bn *f1 = dest; /* F(k) */
bn *f2 = bn_alloc(1); /* F(k+1) */
f1->number[0] = 0;
f2->number[0] = 1;
bn *k1 = bn_alloc(1);
bn *k2 = bn_alloc(1);
/* walk through the digit of n */
for (unsigned int i = 1U << (31 - __builtin_clz(n)); i; i >>= 1) {
/* F(2k) = F(k) * [ 2 * F(k+1) – F(k) ] */
/* F(2k+1) = F(k)^2 + F(k+1)^2 */
bn_lshift(f2, 1, k1);// k1 = 2 * F(k+1)
bn_sub(k1, f1, k1); // k1 = 2 * F(k+1) – F(k)
bn_mult(k1, f1, k2); // k2 = k1 * f1 = F(2k)
bn_mult(f1, f1, k1); // k1 = F(k)^2
bn_swap(f1, k2); // f1 <-> k2, f1 = F(2k) now
bn_mult(f2, f2, k2); // k2 = F(k+1)^2
bn_add(k1, k2, f2); // f2 = f1^2 + f2^2 = F(2k+1) now
if (n & i) {
bn_swap(f1, f2); // f1 = F(2k+1)
bn_add(f1, f2, f2); // f2 = F(2k+2)
}
}
// return f[0]
bn_free(f2);
bn_free(k1);
bn_free(k2);
}
```
沒有使用 bn_cpy 的執行速度如下圖:
:::info
在作業的說明中 bn_lshift 原本沒有
```c
bn_resize(src, src->size + 1);
```
但是我發現這樣算出來的f(92)的結果會是錯誤的，然後我取到 f(1000) 發現後面某些值也是錯的，所以單獨檢查 f(92)。
上圖為當計算到92的最後一個bit時， bn_lshift2(f2, 1, k1) 計算結果是錯的，也就是 printk編號 0 到 1 的時候 k1 應該要變成 f2 的兩倍(2971215073 * 2 = 5,942,430,146) ， 但卻變為 1647462850，導致後面一系列的錯誤，所以問題出在bn_lshift 。 shift > z 時 src 也要resize， 不然會導致當 shift > z 為真時，bn_resize(dest, src->size + 1) 之後 dest 比 src 多出 32 bit，也就是一個 number 然後 for (int i = src->size - 1; i > 0; i--) i 是從 src->size - 1 開始，也就是 dest 最後面那個 number 沒有被賦予值，最後結果錯誤。
:::
下面是用 fdoubling 計算 fb(1) 到 fb(1000) 並且用 perf 去做效能量測的結果。
```
7353,7380 cycles
1,1866,8345 instructions # 1.61 insn per cycle
9,4991 cache-references
3,3307 cache-misses # 35.063 % of all cache refs
0.043825063 seconds time elapsed
0.011958000 seconds user
0.031888000 seconds sys
```
下面是用 fdoubling 且避免使用bn_ cpy 計算 fb(1) 到 fb(1000) 並且用 perf 去做效能量測的結果。
```
1745,6179 cycles
2678,2969 instructions # 1.53 insn per cycle
6,4014 cache-references
2,4079 cache-misses # 37.615 % of all cache refs
0.020755088 seconds time elapsed
0.000000000 seconds user
0.020713000 seconds sys
```
根據作業提示再進一步調整將原本的fast doubling 公式
$$
\begin{split}
F(2k) &= F(k)[2F(k+1) - F(k)] \\
F(2k+1) &= F(k)^2+F(k+1)^2
\end{split}
$$
利用Q-Matrix
$$
\begin{split}
F(2k-1) &= F(k)^2+F(k-1)^2 \\
F(2k) &= F(k)[2F(k-1) + F(k)]
\end{split}
$$
:::info
還在研究為何換個公式就可以減少執行時間...
:::
```c
void bn_fib_fdoubling_Q_Matrix(bn *dest, unsigned int n)
{
bn_resize(dest, 1);
if (n < 2) { // Fib(0) = 0, Fib(1) = 1
dest->number[0] = n;
return;
}
bn *f1 = bn_alloc(1); // f1 = F(k-1)
bn *f2 = dest; // f2 = F(k) = dest
f1->number[0] = 0;
f2->number[0] = 1;
bn *k1 = bn_alloc(1);
bn *k2 = bn_alloc(1);
for (unsigned int i = 1U << (30 - __builtin_clz(n)); i; i >>= 1) {
/* F(2k-1) = F(k)^2 + F(k-1)^2 */
/* F(2k) = F(k) * [ 2 * F(k-1) + F(k) ] */
bn_lshift2(f1, 1, k1); // k1 = 2 * F(k-1)
bn_add(k1, f2, k1); // k1 = 2 * F(k-1) + F(k)
bn_mult(k1, f2, k2); // k2 = k1 * f2 = F(2k)
bn_mult(f2, f2, k1); // k1 = F(k)^2
bn_swap(f2, k2); // f2 <-> k2, f2 = F(2k) now
bn_mult(f1, f1, k2); // k2 = F(k-1)^2
bn_add(k2, k1, f1); // f1 = k1 + k2 = F(2k-1) now
if (n & i) {
bn_swap(f1, f2); // f1 = F(2k+1)
bn_add(f1, f2, f2); // f2 = F(2k+2)
}
}
bn_free(f1);
bn_free(k1);
bn_free(k2);
}
```
再用perf測量一次 時間有下降一些
```
1649,2759 cycles
2510,1465 instructions # 1.52 insn per cycle
6,2615 cache-references
2,7358 cache-misses # 43.692 % of all cache refs
0.006684335 seconds time elapsed
0.000000000 seconds user
0.006707000 seconds sys
```
下圖為有 bn_cpy 和 沒有 bn_cpy 的fast doubling 執行時間比較圖，可以看到後者的表現好非常多
下圖為沒有 bn_cpy 的 fast doubling 和利用Q-Matrix 且沒有 bn_cpy 的 fast doubling 執行時間比較圖，可以看到後者的時間相對較少。
```
7361,9839 cycles
1,1744,5106 instructions # 1.60 insn per cycle
8,7895 cache-references
4,0258 cache-misses # 45.802 % of all cache refs
0.043562180 seconds time elapsed
0.003964000 seconds user
0.039644000 seconds sys
```
```
1564,6279 cycles
2160,7539 instructions # 1.38 insn per cycle
7,3739 cache-references
3,0741 cache-misses # 41.689 % of all cache refs
0.010380267 seconds time elapsed
0.000000000 seconds user
0.010386000 seconds sys
```
```
1421,4606 cycles
1933,1207 instructions # 1.36 insn per cycle
6,7948 cache-references
2,8622 cache-misses # 42.123 % of all cache refs
0.006762957 seconds time elapsed
0.000114000 seconds user
0.006370000 seconds sys
```
Sign in with Wallet
Connect another wallet