---
title: 2023 年 Linux 核心設計/實作課程作業 —— fibdrv
image: https://i.imgur.com/KXUuZLV.png
description: 引導學員開發 Linux 核心模組，預期知悉核心模組如何掛載進入 Linux 核心、Virtual File System (VFS) 概況、character device driver 撰寫，和準備對應的測試及效能分析工具
tags: linux2023
---
# L04: fibdrv
> 主講人: [jserv](http://wiki.csie.ncku.edu.tw/User/jserv) / 課程討論區: [2023 年系統軟體課程](https://www.facebook.com/groups/system.software2023/)
:mega: 返回「[Linux 核心設計/實作](http://wiki.csie.ncku.edu.tw/linux/schedule)」課程進度表
## [sysprog21/bignum](https://github.com/sysprog21/bignum) 程式碼分析
[sysprog21/bignum](https://github.com/sysprog21/bignum) 在 API 的設計想法，區分
* 高階 `bn` API
* 低階 `apm` API
`apm` 表示 Arbitrary-Precision arithmetic (高精度運算)，利用數字陣列來進行大數運算，`uint8_t` 的最大值是 255，該如何表示 300 呢？
1. 令一個 `uint8_t *digit` 的數字陣列來儲存(同時也可以指定此陣列的 size，本例中假設 size = 3)
2. 已知 $300 = 0\times256^2+1\times256^1+44\times256^0$
3. 則 300 可以用 digit[2] = 0,digit[1] = 1,digit[0] = 44 表示
視覺化表示就是
\begin{array}{rrrr}
300: & 00000000 & 00000001 & 00101100 \\
digit: & digit[2] & digit[1] & digit[0] \\
value: & 0 & 1 & 44 \\
\end{array}
最後要輸出要經過格式轉換，才能成為十進位，可見 [format.c](https://github.com/sysprog21/bignum/blob/master/format.c)。
由於除以 2 可以看作是進行 `binary code` 右移，結合 `clz / ctz` 的指令搭配 `fastdoubling` (先去除以數字 MSB 起算的開頭 0 位元，因為真正的數字不包含 leading 0s)，最後呈現出 [sysprog21/bignum](https://github.com/sysprog21/bignum/blob/master/fibonacci.c) 計算費氏數列的實作
```c
bn *a1 = fib; /* Use output param fib as a1 */
bn_t a0, tmp, a;
bn_init_u32(a0, 0); /* a0 = 0 */
bn_set_u32(a1, 1); /* a1 = 1 */
bn_init(tmp); /* tmp = 0 */
bn_init(a);
/* Start at second-highest bit set. */
for (uint64_t k = ((uint64_t) 1) << (62 - __builtin_clzll(n)); k; k >>= 1) {
/* Both ways use two squares, two adds, one multipy and one shift. */
bn_lshift(a0, 1, a); /* a03 = a0 * 2 */
bn_add(a, a1, a); /* ... + a1 */
bn_sqr(a1, tmp); /* tmp = a1^2 */
bn_sqr(a0, a0); /* a0 = a0 * a0 */
bn_add(a0, tmp, a0); /* ... + a1 * a1 */
bn_mul(a1, a, a1); /* a1 = a1 * a */
if (k & n) {
bn_swap(a1, a0); /* a1 <-> a0 */
bn_add(a0, a1, a1); /* a1 += a0 */
}
}
/* Now a1 (alias of output parameter fib) = F[n] */
```
### Karatsuba 演算法
Karatsuba 的概念是將 $x$、$y$ 以第 $n$ 位數為界，拆成兩半 $x_1$、$x_0$、$y_1$、$y_0$，把這他們視為較小的數相乘，然後再透過左移補回 $x_1$、$y_1$ 損失的位數，以十進位為例：
$x = x_1 * 10^n+ x_0 \\ y=y_1*10^n+y_0$
則 $x*y$ 可以化為：
$\underbrace{x_1y_1}_{z_2}*10^{2n}+\underbrace{(x_1y_0+y_1x_0)}_{z_1}*10^n+\underbrace{x_0y_0}_{z_0}$
上述算法計算 $z_2$、$z_1$、$z_0$ 需要四次乘法，我們還可以透過以下技巧優化成三次乘法：
觀察 $(x_1+x_0)(y_1+y_0)$ 展開的結果
$(x_1+x_0)(y_1+y_0)=\underbrace{x_1y_1}_{z_2}+\underbrace{x_1y_0+x_0y_1}_{z_1}+\underbrace{x_0y_0}_{z_0}$
移項之後，我們就能利用 $(x_1+x_0)(y_1+y_0)$、$z_0$、$z_2$ 來計算 $z_1$
$z_2=x_1y_1$
$z_0 = x_0y_0$
$z_1=(x_1+x_0)(y_1+y_0)-z_0-z_2$
最後計算 $z_2*10^{2n}+z_1*10^n+z_0$ 便能得到 $x$ $y$ 相乘的結果。
接著以 2 個 8 位元數值相乘變成 16 bit 為例，解說 Karatsuba 演算法（假設所採用的處理器只支援 8 位元乘法）：
令 $x=01001001_2=73_{10}, \ y=10000011_2=131_{10}$
$x=x_1*2^4+x_0=0100*2^4+1001$
$y=y_1*2^4+y_0=1000*2^4+0011$
$z_2=x_1y_1=0100*1000=00100000$
$z_0=x_0y_0=1001*0011=00011011$
$\begin{align*}
z_1 &=(x_1+x_0)(y_1+y_0)-z_0-z_2 \\
&=(0100+1001)(1000+0011)-00100000-00011011 \\
&=10001111-00100000-00011011 \\
&=01010100
\end{align*}$
$\begin{align*}
x*y &=z_2*2^8+z_1*2^4+z_0 \\
&=0010010101011011 \\
&=9563_{10}
\end{align*}$
其中 $\times 2^n$ 可以用左移運算代替。
接續上一個例子，當 $x$、$y$ 超過 8 bit 時，可以透過分治法實作 Karatsuba。$x_1$、$x_0$、$y_1$、$y_0$ 的位元數超出處理器的乘法的位數時，就把他們再切為 $x_{11}$、$x_{10}$、$x_{01}$、$x_{00}$、...，再使用 Karatsuba 計算。以下以兩個 16 bit 數值相乘變成 32 bit 來演示（假設所採用的處理器只支援 8 位元乘法）：
由上圖可以看出計算 $z_2$、$z_1$、$z_0$ 時，透過分治法將 $x_1$、$x_0$、$y_1$、$y_0$ 切成更小的數字執行乘法運算。最後再用左移與加法計算 $z_2*2^{16} + z_1*2^8 + z_0$ 即可求得結果。
至此可透過分治法，運用 Karatsuba 演算法計算任意位數的大數。
使用 `Karatsuba algorithm` 來加速乘法與平方運算，乘法與平方由於運算較複雜，占用程式執行時間的比例也最高，因此`sysprog21/bignum` 使用 `Karatsuba algorithm` 來加速乘法與平方運算。對 `sysprog21/bignum` 使用 `perf` 進行分析如下 (只擷取部分)
```shell
$ sudo perf report --stdio -g graph,0.5,caller
# To display the perf.data header info, please use --header/--header-only options.
#
#
# Total Lost Samples: 0
#
# Samples: 5K of event 'cycles'
# Event count (approx.): 7673070456
#
# Children Self Command Shared Object Symbol >
# ........ ........ ......... ................. ..............................>
#
75.85% 2.07% fibonacci fibonacci [.] main
|
|--73.78%--main
| |
| |--26.61%--bn_sqr
| | |
| | |--9.77%--_apm_mul_base
| | | |
| | | |--2.89%--apm_dmul
| | | |
| | | --2.13%--apm_dmul_add
| | |
| | |--3.40%--__GI___libc_realloc (inlined)
| | | |
| | | --2.85%--_int_realloc
| | | |
| | | |--0.93%--_int_free
| | | |
| | | --0.85%--_int_malloc
| | |
| | |--2.27%--__GI___libc_malloc (inlined)
| | | |
| | | --0.52%--tcache_get (inlined)
| | |
| | |--2.04%--apm_sqr
| | |
| | |--1.97%--_int_free
| | |
| | |--1.37%--__memmove_avx_unaligned_erms
| | |
| | --0.52%--__GI___libc_free (inlined)
| |
| |--19.35%--bn_mul
| | |
| | |--4.56%--_apm_mul_base
| | | |
| | | |--1.31%--apm_dmul
| | | |
| | | --0.97%--apm_dmul_add
| | |
| | |--3.90%--__GI___libc_malloc (inlined)
| | | |
| | | |--1.20%--tcache_get (inlined)
| | | |
| | | --0.56%--checked_request2size (inlin>
| | |
| | |--2.56%--_int_free
| | |
| | |--1.80%--apm_mul
| | |
| | |--1.25%--__GI___libc_realloc (inlined)
| | | |
| | | --0.96%--_int_realloc
| | |
| | |--1.12%--__memmove_avx_unaligned_erms
| | |
| | --0.68%--__GI___libc_free (inlined)
| |
| |--16.17%--bn_add
| | |
| | |--3.82%--apm_add_n
| | |
| | |--3.45%--__GI___libc_realloc (inlined)
| | | |
| | | --2.77%--_int_realloc
| | | |
| | | |--0.90%--_int_malloc
| | | |
| | | --0.60%--_int_free
| | |
| | --3.08%--apm_add
| |
| |--3.87%--bn_lshift
| | |
| | |--1.31%--apm_lshift
| | |
| | --0.77%--__memset_avx2_unaligned_erms
| |
| |--1.99%--bn_init
| | |
| | --1.54%--__libc_calloc
| | |
| | --0.71%--_int_malloc
| |
| |--1.83%--_int_free
| |
| |--1.62%--bn_swap
| |
| --1.21%--bn_init_u32
| |
| --1.08%--bn_init
| |
| --0.91%--__libc_calloc
|
|--1.30%--_start
| __libc_start_main
| main
|
--0.77%--0xffffffffffffffff
main
```
可見使用 Karatsuba algorithm 的 [sysprog21/bignum](https://github.com/sysprog21/bignum)，其乘法與平方是程式執行時間佔比最高的運算，分別為 19.35% 與 26.61%，合計近 46% 的時間占比。
### [Schönhage–Strassen Algorithm](https://en.wikipedia.org/wiki/Sch%C3%B6nhage%E2%80%93Strassen_algorithm)
Schönhage–Strassen 的概念是將 $x$、$y$ 以 $n$ 個位數為單位，將大數拆成若干個較小的數 $x_0$、$x_1$...、$y_0$、$y_1$、...，對 $(x_0, x_1, ...)$ 和 $(y_0, y_1,...)$ 線性摺積，執行完線性摺積後將所有元素依據自己的位數左移並相加就可以完成乘法運算。舉例：
令 $x = 26699, y = 188$
將 $x$、$y$ 以 **2** 個位數為界分割會得到以下兩個數列：
$(2, 66, 99),\ (1, 88)$
對數列線性摺積：
```
2 66 99
x 1 88
-----------------------
176 5808 8712
2 66 99
-----------------------
2 242 5907 8712
```
$(2,242,5907,8712)$ 為 $(2, 66, 99)$ 與 $(1, 88)$ 的線性摺積，接下來透過左移和加法將線性摺積的結果相加即可完成大數乘法：
```
8712
5907
242
+ 2
--------------
5019412
```
Sign in with Wallet
Connect another wallet