---
tags: linux2023
---
# 2023q1 Homework3 (fibdrv)
contributed by < `YSRossi` >
## 開發環境
```shell
$ gcc --version
gcc (Ubuntu 11.3.0-1ubuntu1~22.04) 11.3.0
$ lscpu
Architecture: x86_64
CPU op-mode(s): 32-bit, 64-bit
Address sizes: 48 bits physical, 48 bits virtual
Byte Order: Little Endian
CPU(s): 8
On-line CPU(s) list: 0-7
Vendor ID: AuthenticAMD
Model name: AMD FX-8320E Eight-Core Processor
CPU family: 21
Model: 2
Thread(s) per core: 2
Core(s) per socket: 4
Socket(s): 1
Stepping: 0
Frequency boost: enabled
CPU max MHz: 3200.0000
CPU min MHz: 1400.0000
BogoMIPS: 6429.64
```
## 開發紀錄
### 改成 Fast Doubling 版本（成功計算到 F~92~）
參考作業說明所提供 Fast Doubling 的虛擬碼，搭配 bitwise 實作。
- `n` 的型態為 `long long`，使用 `sizeof(long long) << 3` 計算 `long long` 所佔有的位元數，減去 `n` 的 `leading 0-bits` 數量，得到 `n` 實際上有用到的位元數。其中 `n` 的 `leading 0-bits` 數量，使用 `__builtin_clzll(n)` 獲得。
- 在迴圈中，目前的位元以 `n & (1 << (i - 1))` 表示，將 1 左移 `i - 1` 把對應目前位元的位置設成 1 ，與 `n` 做 `and` 判斷目前位元是否為 1 。
```
Fast_Fib(n)
a = 0; b = 1; // m = 0
for i = (number of binary digit in n) to 1
t1 = a*(2*b - a);
t2 = b^2 + a^2;
a = t1; b = t2; // m *= 2
if (current binary digit == 1)
t1 = a + b; // m++
a = b; b = t1;
return a;
```
- 實驗比較 Fast Doubling 版本與原本方法在 kernel space 所花費的時間（還未排除干擾效能分析的因素）
:::warning
TODO : 嘗試解決在 GRUB 的設定檔中加入 ```isolcpus=7``` 重新開機後仍無法起作用的問題。
```
GRUB_CMDLINE_LINUX_DEFAULT="quiet splash isolcpus=7"
```
> 已解決，使用 ```sudo update-grub``` 更新 GRUB 設定檔
:::
- 實驗比較 Fast Doubling 與原始方法在 kernel space 所花費的時間(已排除干擾效能分析的因素)
```c=
static long long fib_sequence_fast_doubling(long long k)
{
long long a = 0, b = 1;
for (int i = 63; i >= 0; i--) {
long long t1 = a * ((2 * b) - a);
long long t2 = b * b + a * a;
a = t1;
b = t2;
if ((k >> i) & 1) {
t1 = a + b;
a = b;
b = t1;
}
}
return a;
}
```
- 若將上方程式碼的第 5 行改成下方程式碼，計算 fib(0) 在 user space 會花費大量時間。
```c
for (int i = 63 - __builtin_clzll(k); i >= 0; i--) {
```
`write` 回傳 Fast Doubling 版本的 kernel space time
`read` 回傳 Fast Doubling 版本計算的 Fibonacci 值
```c
/* client.c */
for (int i = 0; i <= offset; i++) {
struct timespec start, end;
lseek(fd, i, SEEK_SET);
clock_gettime(CLOCK_MONOTONIC, &start);
sz = read(fd, buf, 1);
clock_gettime(CLOCK_MONOTONIC, &end);
printf("Reading from " FIB_DEV " at offset %d, returned the sequence ",
i);
printf("%lld.\n", sz);
kernel_time = write(fd, write_buf, 0);
user_time = (long long)(end.tv_sec * 1e9 + end.tv_nsec)
- (start.tv_sec * 1e9 + start.tv_nsec);
}
```
- 在計算 fib(0) 會有異常的原因是 [__builtin_clzll](https://gcc.gnu.org/onlinedocs/gcc/Other-Builtins.html#Other-Builtins) 的參數不可以為0，因此在計算 fib(0) 時，直接回傳 0 即可解決問題。
> Returns the number of leading 0-bits in x, starting at the most significant bit position. If x is 0, the result is undefined.
### 排除干擾效能分析的因素
* 限定 CPU 給特定的程式使用
使用 ```isolcpus``` Linux 核心起始參數，讓特定的 CPU 核在開機時就被保留下來，設定完成後使用 ```sudo update-grub``` 更新 grub 設定檔
```shell
GRUB_CMDLINE_LINUX_DEFAULT="quiet splash isolcpus=7"
```
* 設定 [NO_HZ](https://www.kernel.org/doc/Documentation/timers/NO_HZ.txt) 減少影響程式執行的中斷，設定完成後使用 ```sudo update-grub``` 更新 grub 設定檔
```shell
GRUB_CMDLINE_LINUX_DEFAULT="quiet splash isolcpus=7 nohz_full=7"
```
* 抑制 [address space layout randomization](https://en.wikipedia.org/wiki/Address_space_layout_randomization) (ASLR)
```shell
$ sudo sh -c "echo 0 > /proc/sys/kernel/randomize_va_space"
```
* 設定 scaling_governor 為 performance。準備以下 shell script，檔名為 `performance.sh`:
```shell
for i in /sys/devices/system/cpu/cpu*/cpufreq/scaling_governor
do
echo performance > ${i}
done
```
之後再用 `$ sudo sh performance.sh` 執行
* 關閉 [Processor boosting](https://www.kernel.org/doc/Documentation/cpu-freq/boost.txt)
```shell
sudo bash -c "echo 0 > /sys/devices/system/cpu/cpufreq/boost"
```
* 設定 [SMP IRQ Affinity](https://www.kernel.org/doc/html/next/core-api/irq/irq-affinity.html) ，減少目標 cpu 處理 interrupt request ，使用 cpu7 來做此實驗，所以 mask 為 0x7f
```shell
sudo bash -c "echo 7f > /proc/irq/default_smp_affinity"
```
設定完成後，雖然還是有些浮動，但可以測量到相對穩定的實驗結果
### 支援大數運算（成功計算到 F~186~）
#### 大數資料結構
`uppper` 紀錄大數較高的 `64` 位元，`lower` 紀錄大數較低的 `64` 位元
```c
struct BigN {
unsigned long long upper, lower;
};
```
依照前面實作， Fast-doubling 需要加法、減法、乘法、左移一位運算，其中大數乘法運算還會使用右移一位運算。
##### `bn_add` 加法
兩數 `upper` 與 `lower` 分別相加，判斷 `lower` 相加後是否溢位，溢位的話將 `upper` 加一。
```c
struct BigN bn_add(struct BigN a, struct BigN b)
{
struct BigN result;
result.lower = a.lower + b.lower;
result.upper = a.upper + b.upper;
if (a.lower > ~b.lower)
result.upper++;
return result;
}
```
##### `bn_sub` 減法
兩數 `upper` 與 `lower` 分別相減，若 `a.lower < b.lower`，代表需要跟 `upper` 借位，所以 `upper` 減一。
```c
struct BigN bn_sub(struct BigN a, struct BigN b)
{
struct BigN result;
result.upper = a.upper - b.upper;
if (a.lower < b.lower)
result.upper--;
result.lower = a.lower - b.lower;
return result;
}
```
##### `bn_lshift` 左移一位 / `bn_rshift` 右移一位
實作左移一位，將兩數 `upper` 與 `lower` 分別左移一位，若 `lower` 在位移前的最高有效位元為 `1`，則 `upper` 的最低有效位元設為 `1`。實作右移一位與左移相似，兩數 `upper` 與 `lower` 分別右移一位，若 `upper` 在位移前的最低有效位元為 `1`，則 `lower` 的最高有效位元設為 `1`。
```c
struct BigN bn_lshift(struct BigN a)
{
a.upper <<= 1;
if (a.lower & MSB_MASK)
a.upper |= 1;
a.lower <<= 1;
return a;
}
struct BigN bn_rshift(struct BigN a)
{
a.lower >>= 1;
if (a.upper & 1)
a.lower |= MSB_MASK;
a.upper >>= 1;
return a;
}
```
##### `bn_mul` 乘法
每次將 `a` 左移一位，並將 `b` 右移一位，找出 `b` 中值為 `1` 的位元，若找到為 `1` 的位元，則將目前 `a` 的值累加到 `result` ，達到類似乘法直式的運算方式。
例如 : a = 0b1010, b = 0b101
a * b = 0b1010 * 0b101 = 0b1010 + 0b101000 = 0b110010
```c
struct BigN bn_mul(struct BigN a, struct BigN b)
{
struct BigN result = bn_init(0, 0);
while (b.upper || b.lower) {
if (b.lower & 1)
result = bn_add(result, a);
a = bn_lshift(a);
b = bn_rshift(b);
}
return result;
}
```
##### `bn_to_str` 將大數轉成字串
`input` 為要轉換成字串的大數， `num` 為最靠近且小於等於 `input` 的 `10` 的冪，將 `input` 不停減掉 `num` 直到 `input < num` ，減去的次數即目前 `input` 最高位的十進位數。
例如 : `input = 3021`, 最靠近且小於等於 `input` 的 `10` 的冪為 `num = 1000` ， `input` 需要減掉 `num` `3` 次才會滿足 `input < num` 的條件，因此得到最高位的十進位數為 `3` ，後面的位數以此類推。
```c
int bn_to_str(char *buf, char *str)
{
struct BigN input =
bn_init(*(unsigned long long *) buf, *((unsigned long long *) buf + 1));
struct BigN num = bn_init(0, 1);
struct BigN ten = bn_init(0, 10);
/* maximum value of BigN is 340,282,366,920,938,463,463,374,607,431,768,211,455 (39 digits)*/
int pos = 39; /* highest possible digit index*/
int dig_idx = 0;
int dig_num;
while (pos--) {
for (dig_num = 0, num = bn_init(0, 1); dig_num < pos; dig_num++)
num = bn_mul(num, ten);
if (dig_idx || !bn_less(input, num) || !pos) {
for (dig_num = 0; !bn_less(input, num); dig_num++)
input = bn_sub(input, num);
str[dig_idx++] = dig_num + '0';
}
}
str[dig_idx] = 0;
return dig_idx;
}
```
#### 實驗結果
### 支援大數運算(可動態改變大數長度)
#### 大數資料結構
`digits` 紀錄大數的數值，`size` 紀錄大數使用幾個 uint64_t，`alloc` 配置幾個 uint64_t 空間。
```c
typedef u64 bn_digit;
typedef u64 bn_size;
struct {
bn_digit *digits;
bn_size size;
bn_size alloc;
}
```
##### 重新配置空間
若 `size` 超過 `alloc` 時，代表配置的空間不夠，需要 realloc，`BN_MIN_ALLOC(n, s)` 將大數 `n` 重新配置成 `s` 個 uint64_t 的大小，並確保配置的大小為 4 的倍數個 uint64_t。
`BN_SIZE(n, s)` 的功能類似，會額外將大數 `n` 的 `size` 設為 `s`。
```c
#define BN_MIN_ALLOC(n, s) \
do { \
bn *const __n = (n); \
const bn_size __s = (s); \
if (__n->alloc < __s) { \
__n->digits = \
bn_resize(__n->digits, __n->alloc = ((__s + 3) & ~3U)); \
} \
} while (0)
#define BN_SIZE(n, s) \
do { \
bn *const __n = (n); \
__n->size = (s); \
if (__n->alloc < __n->size) { \
__n->digits = \
bn_resize(__n->digits, __n->alloc = ((__n->size + 3) & ~3U)); \
} \
} while (0)
```
##### `bn_add` 加法
將大數 `a` 與大數 `b` 對應的位數相加，若有進位會累加到下一位。為了確保 a = c 與 b = c 時能正確運算不改到 `a` 或 `b` 的值，使用 `tmp` 紀錄結果，最後再將 `tmp` 複製到 `c`。
```c
/* a + b = c */
void bn_add(bn *a, bn *b, bn *c)
{
/* a + 0 or 0 + b */
if (a->size == 0) {
if (b->size == 0)
c->size = 0;
else
bn_set(c, b);
return;
} else if (b->size == 0) {
bn_set(c, a);
return;
}
/* a + b = 2 * a */
if (a == b) {
bn_lshift(a, c);
return;
}
bn_size min_size, max_size;
bn_digit *tmp;
bn *max_bn;
if (a->size <= b->size) {
min_size = a->size;
max_size = b->size;
max_bn = b;
} else {
min_size = b->size;
max_size = a->size;
max_bn = a;
}
tmp = kmalloc((max_size + 1) * BN_DIGITS_SIZE, GFP_KERNEL);
for (int i = 0; i < max_size + 1; i++)
tmp[i] = 0;
bn_digit cy = 0;
for (int i = 0; i < min_size; i++) {
bn_digit u = a->digits[i];
bn_digit v = b->digits[i];
cy = (u += cy) < cy;
tmp[i] = u + v;
cy += tmp[i] < v;
}
if (cy) {
for (int i = min_size; i < max_size; i++) {
bn_digit u = max_bn->digits[i];
cy = (u += cy) < cy;
tmp[i] = u;
}
if (cy)
tmp[max_size++] = cy;
} else
bn_copy(max_bn->digits + min_size, max_size - min_size, tmp + min_size);
BN_SIZE(c, max_size);
bn_copy(tmp, max_size, c->digits);
kfree(tmp);
}
```
##### `bn_lshift` 左移一位
從低位到高位尋訪大數的 `digits`，每次左移 1 位，藉由將 `digits` 右移 63 位獲取該 `digits` 的最高位，並與下一個 `digits` 做 OR 運算，將下一個 `digits` 的最低位設成目前 `digits` 的最高位。
```c
/* c = a << 1 */
void bn_lshift(bn *a, bn *c)
{
if (a->size == 0) {
bn_set(c, a);
return;
}
if(a != c)
BN_SIZE(c, a->size);
const unsigned int subp = 63;
bn_digit cy = 0;
bn_size size = a->size;
for (int i = 0; i < size; i++) {
bn_digit p = *(a->digits + i);
*(c->digits + i) = (p << 1) | cy;
cy = p >> subp;
}
if (cy) {
BN_SIZE(c, c->size + 1);
c->digits[c->size - 1] = cy;
}
}
```
#### `bn_mul` 乘法
大數 `a` 的每個位數與大數 `b` 的每個位數相乘，將結果紀錄在 `tmp` ，二個 64-bit 整數相乘會需要 128-bit ，所以使用 GCC Extension `__uint128_t`，透過右移 64 位，將較高的 64-bit 存在 `tmp` 目前位數的下一位，最後再將 `tmp` 複製到 `c`。
```c
void bn_mul(bn *a, bn *b, bn *c)
{
if (a->size == 0 || b->size == 0)
c->size = 0;
bn_size csize = a->size + b->size;
bn_digit *tmp = kmalloc(csize * BN_DIGITS_SIZE);
for (int i = 0; i < csize; i++)
tmp[i] = 0;
for (int i = 0; i < a->size; i++) {
for (int j = 0; j < b->size; j++) {
__int128 product = (__int128)a->digits[i] * (__int128)b->digits[j];
bn_size low = product;
bn_size high = product >> BN_DIGITS_BITS;
tmp[i + j] += low;
uint32_t cy = tmp[i + j] < low;
tmp[i + j + 1] += high + cy;
if (cy)
cy = tmp[i + j + 1] <= high;
else
cy = tmp[i + j + 1] < high;
for (int k = i + j + 2; cy; k++) {
tmp[k] += cy;
cy = tmp[k] < cy;
}
}
}
while (csize > 0 && tmp[csize - 1] == 0)
csize--;
BN_SIZE(c, csize);
bn_copy(tmp, csize, c->digits);
kfree(tmp);
}
```
##### `bn_to_str` 將大數轉成字串
基本想法為對大數除以 10 取餘數，反覆對其商除以 10 取餘數，將取得的餘數逐一串接起來，最後所得到的就是該數的十進位表示法。
```
123 / 10 = 12...3 <-
12 / 10 = 1...2 <-
1 / 10 = 0...1 <-
```
- 方法一(對應到圖中的 Not divided by max_radix)
每次迴圈將大數除以 10 取餘數，這個餘數就是該數的十進位最低位的數值，保留商到下一次迴圈再除以 10，直到商等於 0 時停止迴圈。
方法一對大數一次只除以 10，需要跑很多次迴圈才可以讓商等於 0，因此有了方法二來加速。
- 方法二(對應到圖中 Divided by max_radix)
`max_radix` 為 10^19^，相當於 uint64_t 可表示值的範圍中，最大的 10 的冪的值，`max_power` 為 19，紀錄 10 的冪次。
每次迴圈藉由將大數除以 `max_radix`，獲得一個 uint64_t 可表示的餘數，再用方法一來取得最終的餘數，重複上述操作直到商等於 0 時停止迴圈。
```c
static char *bn_to_str(const bn_digit *u,
bn_size size,
char *out)
{
while (size && !u[size -1]) /* 取得 u 的有效 size(byte) */
--size;
/* u = 0 ~ 9 */
if (size == 0 || (size == 1 && u[0] < 10)) {
if (!out)
out = malloc(2);
out[0] = size ? radix_chars[u[0]] : '0';
out[1] = '\0';
return out;
}
const bn_digit max_radix = radix_table.max_radix;
const unsigned int max_power = radix_table.max_power;
if (!out)
out = malloc(string_size(size) + 1);
char *outp = out;
bn_digit *tmp = BN_TMP_COPY(u, size);
bn_size tsize = size;
do {
bn_digit r = bn_ddivi(tmp, tsize, max_radix);
tsize -= (tmp[tsize - 1] == 0);
unsigned int i = 0;
do {
bn_digit rq = r / 10;
bn_digit rr = r % 10;
*outp++ = radix_chars[rr];
r = rq;
if (tsize == 0 && r == 0)
break;
} while (++i < max_power);
} while (tsize != 0);
free(tmp);
char *f = outp - 1;
while (*f == '0')
--f;
f[1] = '\0';
for (char *s = out; s < f; ++s, --f)
SWAP(*s, *f);
return out;
}
```
#### 降低 realloc 呼叫次數
在 [sysprog21/bignum](https://github.com/sysprog21/bignum) 中，初始化大數所 malloc 的大小皆相同，在計算 fib(k) 時，k 越大，所需的空間越多，realloc 的次數也越多。
透過預估 fib(k) 所需的位數，在初始化時就分配足夠的空間，降低 realloc 次數。
根據[作業說明](https://hackmd.io/@sysprog/linux2023-fibdrv/%2F%40sysprog%2Flinux2023-fibdrv-b)，當 $n$ 為足夠大的正整數時，
$$
F(n)\approx0.4472135955\times1.61803398875^n
$$
將近似值取 2 為底的對數，可得到需要使用的位元數
$$
log_2(0.4472135955\times1.61803398875^n)\approx-1.160964 + n\times0.694242
$$
除以 64 可得到需要使用的 `uint64` 數量
$$
(-1.160964 + n\times0.694242)\div64 = -0.0181400625 + n\times0.01084753125
$$
$$\lfloor\frac{-181400625 + n\times108475313}{10^{10}}\rfloor+1
$$
```c
int bn_fib_digits_n(int n)
{
if (n < 2)
return 1;
int digits = ((long)n * 108475313 - 181400625) / 10000000000 + 1;
return digits;
}
```
由下圖能觀察到隨著數字增加，預先分配足夠空間的效果越明顯。
### 關於 [sysprog21/bignum](https://github.com/sysprog21/bignum) 的問題
:::info
在 [sysprog21/bignum/bignum.c](https://github.com/sysprog21/bignum/blob/master/bignum.c) 中，macro 使用 do while(0) 將多行內容框起來，我想原因是減少語法錯誤，像是下方在 if 中呼叫 macro，且沒有用大括號括住的情形依舊能正確執行，若 macro 沒有 do while(0) 則會出錯。
在 macro 中不用 if(1) { ... } 而是使用 do while(0) 將內容框起來的原因是單純 do while(0) 後可以加分號，因此使用該 macro 時可以較直觀的在後面加分號，還是有其他考量？
```c
#define BN_MIN_ALLOC(n, s) \
do { \
... \
} while (0)
if (1)
BN_MIN_ALLOC(n, s);
```
:::
:::info
在 [sysprog21/bignum/format.c](https://github.com/sysprog21/bignum/blob/master/format.c) 中，為什麼當 radix 不是 2 的冪時，結果需要額外加 2，而不是加 1？
```c
static size_t apm_string_size(apm_size size, unsigned int radix)
{
ASSERT(radix >= 2);
ASSERT(radix <= 36);
if ((radix & (radix - 1)) == 0) {
unsigned int lg = apm_digit_lsb_shift(radix);
return ((size * APM_DIGIT_BITS + lg - 1) / lg) + 1;
}
/* round up to second next largest integer */
return (size_t)(radix_sizes[radix] * (size * APM_DIGIT_SIZE)) + 2;
}
```
:::
Sign in with Wallet
Connect another wallet