---
tags: linux2023
---

# 2023q1 Homework1 (fibdrv)
contributed by < [Jerejere0808](https://github.com/Jerejere0808) >

新增以下資料結構使 fibonnaci 可以進行大數運算，其中 capacity 是用於分配多餘的記憶體空間，這樣就可以避免多次的 resize (realloc) ， 只有在 capacity 不足時才會需要重新分配記憶體空間。

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。
:::

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 所需多少計算時間

kt = ktime_get();
bn_fib_fdoubling(result, *offset);
kt = ktime_sub(ktime_get(), kt);

參考 [colinyoyo26](https://github.com/colinyoyo26/fibdrv) 的方法用 clock_gettime 在 client 取得 user 所需的時間

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 如下

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;
}

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 實作如下:

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;
}

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 原本沒有

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
還在研究為何換個公式就可以減少執行時間...
:::

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