---
tags: linux2022
---
# sysprog21/bignum 的研究和改進
contributed by <`curlyw819` >
## 實驗環境
```
$ gcc --version
gcc (Ubuntu 9.4.0-1ubuntu1~20.04.1) 9.4.0
$ lscpu
架構: x86_64
CPU 作業模式: 32-bit, 64-bit
Byte Order: Little Endian
Address sizes: 39 bits physical, 48 bits virtual
CPU(s): 12
On-line CPU(s) list: 0-11
每核心執行緒數: 2
每通訊端核心數: 6
Socket(s): 1
NUMA 節點: 1
供應商識別號: GenuineIntel
CPU 家族: 6
型號: 158
Model name: Intel(R) Core(TM) i7-8700 CPU @ 3.20GHz
製程: 10
CPU MHz: 3200.000
CPU max MHz: 4600.0000
CPU min MHz: 800.0000
虛擬: 6399.96
Virtualization: VT-x
L1d cache: 192 KiB
L1i cache: 192 KiB
L2 cache: 1.5 MiB
L3 cache: 12 MiB
NUMA node0 CPU(s): 0-11
```
## 作業要求
對比 [`sysprog21/bignum`](https://github.com/sysprog21/bignum) 和自己撰寫的大數運算[2022q1 Homework3(fibdrv)](https://hackmd.io/T9OVVGNcRRS0o0lmfz2ZCA)，同樣在求 `Fibonacci` 數，分析運算速度、記憶體開銷，和潛在的錯誤，並解釋何以有這些落差。研讀 `sysprog21/bignum`，落實原始程式碼中標注的待做事項。
## 開發進度
由於我的 `fibnacci` 運算沒有使用到 `fastdoubling` 與 `clz` 加速，原本在前提不相同的情況下是無法與使用了 `fastdoubling` 與 `clz` 加速過的 `sysprog21/bignum` 去做比對。因此這部分的研究重點會放在重新理解 `fastdoubling` 是如何加速，並還是與自己的程式去做速度上的比對。
* 從計算費式數列的比較圖中可以發現，在線的趨勢上就有明顯的落差，自己的程式隨著計算量的增加，時間耗費也明顯增長。
### `sysprog21/bignum` 程式分析
* 在四則運算的部分，`sysprog21/bignum` 採用高階 `(bn)` 與低階 `(amp)` 的 API 包裝而成
* 在費式數列的計算上，原本使用 `iterative` 的方式計算費式數列 ，迴圈迭代次數與欲求費式數成正比，時間複雜度爲 **O(n)**。使用 `fast doubling` ，由於是 bit 的二進制運算，至多只要迭代 64(或32，依設定有所不同)次，實際上去除 MSB 為 0 的 bit 不用做迭代，時間複雜度為 **O(log n)**。 `sysprog21/bignum` 的 `fast doubling` 使用了不同的 `Q-matrix`(參考自[KYG-yaya573142 的報告](https://hackmd.io/@KYWeng/rkGdultSU#%E5%8E%9F%E6%9C%AC%E7%9A%84%E9%81%8B%E7%AE%97%E6%88%90%E6%9C%AC))，推導如下
$$
\begin{split}
\begin{bmatrix}
F(2n-1) \\
F(2n)
\end{bmatrix} &=
\begin{bmatrix}
0 & 1 \\
1 & 1
\end{bmatrix}^{2n}
\begin{bmatrix}
F(0) \\
F(1)
\end{bmatrix}\\ \\ &=
\begin{bmatrix}
F(n-1) & F(n) \\
F(n) & F(n+1)
\end{bmatrix}
\begin{bmatrix}
F(n-1) & F(n) \\
F(n) & F(n+1)
\end{bmatrix}
\begin{bmatrix}
1 \\
0
\end{bmatrix}\\ \\ &=
\begin{bmatrix}
F(n)^2 + F(n-1)^2\\
F(n)F(n) + F(n)F(n-1)
\end{bmatrix}
\end{split}
$$
整理後可得
$$
\begin{split}
F(2k-1) &= F(k)^2+F(k-1)^2 \\
F(2k) &= F(k)[2F(k-1) + F(k)]
\end{split}
$$
使用上述公式改寫 `fdoubling` 函式，優點是可以少掉一次迴圈的計算及避免使用減法
* 引入 `bn_sqr` (參考自[KYG-yaya573142 的報告](https://hackmd.io/@KYWeng/rkGdultSU#%E5%8E%9F%E6%9C%AC%E7%9A%84%E9%81%8B%E7%AE%97%E6%88%90%E6%9C%AC))
```
a b c
x a b c
-------------------
ac bc cc
ab bb bc
aa ab ac
```
考慮上述 $abc^2$ 的計算過程，會發現數值 $ab$ 、 $ac$ 與 $bc$ 各會重複一次，因此可先計算對角線其中一邊的數值，將數值的總和直接乘二，最終再加上對角線上的 $aa$ 、 $bb$ 與 $cc$。藉由這種方式，平方運算的成本可由本來的 $n^2$ 次乘法降為 $(n^2 - n)/2$ 次乘法
* `clz / ctz`
由於除 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 algorithm` 來加速乘法與平方運算
計算費式數列時，乘法與平方由於運算較複雜，占用程式執行時間的比例也最高，因此`sysprog21/bignum` 使用 `Karatsuba algorithm` 來加速乘法與平方運算。對 `sysprog21/bignum` 使用 `perf` 進行分析如下(只擷取部分重點資訊)
```
$ 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
```
可以看到即使`sysprog21/bignum` 使用 `Karatsuba algorithm` 來加速乘法與平方運算，乘法與平方仍是程式執行時間佔比最高的運算，分別為 19.35% 與 26.61%，合計近 46% 的時間占比。
`Karatsuba algorithm` 原理的部分可以參見 [blueskyson](https://hackmd.io/@blueskyson/linux2022-fibdrv) 同學講解 `Karatsuba algorithm` 的部分，說明得清楚又好理解
### 自己的([curlyw819](https://hackmd.io/T9OVVGNcRRS0o0lmfz2ZCA)) 程式分析
#### 運算速度
首先，自己的程式在費式數列的運算次數上，使用傳統的 `iterative`，時間複雜度為 O(n)，而 `sysprog21/bignum` 為 `fast doubling` 的 O(log n)，撇除費式數列的計算方式不同，在 `bignum` 的實作上也有差異。
`sysprog21/bignum` 採用高階 `(bn)` 與低階 `(amp)` 的 API 包裝而成，其中 `amp` 使用二進制進行 `bitwise` 運算，為整體程式提供許多速度上的優化
apm.c
```c
/* Set u[size] = u[size] + 1, and return the borrow. */
static apm_digit apm_inc(apm_digit *u, apm_size size)
{
if (size == 0)
return 1;
for (; size--; ++u) {
if (++*u)
return 0;
}
return 1;
}
/* Set u[size] = u[size] + v and return the carry. */
apm_digit apm_daddi(apm_digit *u, apm_size size, apm_digit v)
{
if (v == 0 || size == 0)
return v;
return ((u[0] += v) < v) ? apm_inc(&u[1], size - 1) : 0;
}
/* Set w[size] = u[size] + v[size] and return the carry. */
apm_digit apm_add_n(const apm_digit *u,
const apm_digit *v,
apm_size size,
apm_digit *w)
{
ASSERT(u != NULL);
ASSERT(v != NULL);
ASSERT(w != NULL);
apm_digit cy = 0;
while (size--) {
apm_digit ud = *u++;
const apm_digit vd = *v++;
cy = (ud += cy) < cy;
cy += (*w = ud + vd) < vd;
++w;
}
return cy;
}
/* Set w[max(usize, vsize)] = u[usize] + v[vsize] and return the carry. */
apm_digit apm_add(const apm_digit *u,
apm_size usize,
const apm_digit *v,
apm_size vsize,
apm_digit *w)
{
ASSERT(u != NULL);
ASSERT(usize > 0);
ASSERT(u[usize - 1]);
ASSERT(v != NULL);
ASSERT(vsize > 0);
ASSERT(v[vsize - 1]);
if (usize < vsize) {
apm_digit cy = apm_add_n(u, v, usize, w);
if (v != w)
apm_copy(v + usize, vsize - usize, w + usize);
return cy ? apm_inc(w + usize, vsize - usize) : 0;
} else if (usize > vsize) {
apm_digit cy = apm_add_n(u, v, vsize, w);
if (u != w)
apm_copy(u + vsize, usize - vsize, w + vsize);
return cy ? apm_inc(w + vsize, usize - vsize) : 0;
}
/* usize == vsize */
return apm_add_n(u, v, usize, w);
}
/* Set u[usize] = u[usize] + v[vsize] and return the carry. */
apm_digit apm_addi(apm_digit *u,
apm_size usize,
const apm_digit *v,
apm_size vsize)
{
ASSERT(u != NULL);
ASSERT(v != NULL);
ASSERT(usize >= vsize);
apm_digit cy = apm_addi_n(u, v, vsize);
return cy ? apm_inc(u + vsize, usize - vsize) : 0;
}
```
`apm_inc` 會回傳一個 `carry` ，而 `carry` 的計算方式為 : 會有一個由 `LSB` 走往 `MSB` 的指標，直到走到相加不為 0 的位數，則回傳進位
`apm_addi` 去除不必要的資訊
`apm_add` 兩個變數的 usize、vsize 相減，就可以在計算結果把高位元多出來的部分直接移到結果上
`apm_add_n` 若 usize == vsize 則使用 apm_add_n
#### 空間使用
我的程式在計算大位數的數字時就需要開很大的陣列，如下
```c
void bn_add(bignum *x, bignum *y, bignum *dest)
{
if (x->length < y->length)
return bn_add(y, x, dest);
int carry = 0;
for (int i = 0; i < y->length; i++) {
dest->decimal[i] = y->decimal[i] + x->decimal[i] + carry - OFFSET;
if (dest->decimal[i] >= 10 + OFFSET) { // OFFSET = 48
carry = 1;
dest->decimal[i] -= 10;
} else
carry = 0;
}
for (int i = y->length; i < x->length; i++) {
dest->decimal[i] = x->decimal[i] + carry;
if (dest->decimal[i] >= 10 + OFFSET) {
carry = 1;
dest->decimal[i] -= 10;
} else
carry = 0;
}
dest->length = x->length;
dest->decimal[dest->length] = carry + OFFSET;
if (carry)
dest->length++;
dest->decimal[dest->length] = '\0';
}
```
且每個數字都需要用到 2 個 byte 來儲存，若 fib(500) 有 105 個數字，就需要用到 210 byte，記憶體開銷非常大，且有潛在的風險，當fib往上加，就可能會有溢位的問題
而 `sysprog21/bignum` 固定就是 64 bit(or 32 bit)，並且在要輸出結果時才一口氣轉換成 10 進制，節省了不少時間上的開銷
### 落實原始程式碼中標注的待做事項
由於個人能力不足，在經過一陣子嘗試後仍沒有進度，因此無法完成
Sign in with Wallet
Connect another wallet