# 2021q1 Homework3 (fibdrv)
contributed by < [`bakudr18`](https://github.com/bakudr18) >
###### tags: `linux2021`
> [fibdrv](https://hackmd.io/@sysprog/linux2021-fibdrv)
## 建立效能分析實驗
### 環境設定
#### 於指定 CPU 核心上執行程式
首先於 `/etc/default/grub` 新增 `isolcpu=7` 以指定 cpu 7 不被用來做 kernel schduler 排程，設定完後執行 `sudo update-grub` 然後重開機，可透過以下命令確認是否設定成功
```shell
$ cat /sys/devices/system/cpu/isolated
7
```
接著以 [taskset](http://man7.org/linux/man-pages/man1/taskset.1.html) 指定 CPU 核心來執行程式：
```shell
$ taskset -c 7 EXECUTABLE
```
#### 排除效能干擾因素
* 抑制 [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`
* 針對 Intel 處理器，關閉 turbo mode
```shell
$ sudo sh -c "echo 1 > /sys/devices/system/cpu/intel_pstate/no_turbo"
```
#### SMP IRQ affinity
> 參考自 [KYG-yaya573142](https://hackmd.io/@KYWeng/rkGdultSU#%E6%8E%92%E9%99%A4%E5%B9%B2%E6%93%BE%E6%95%88%E8%83%BD%E5%88%86%E6%9E%90%E7%9A%84%E5%9B%A0%E7%B4%A0) 的筆記
[SMP IRQ afffinity](https://www.kernel.org/doc/Documentation/IRQ-affinity.txt) 是用來指定 IRQ# 可以使用的 CPU 核心，執行以下 script
```shell
#!/bin/bash
for file in `find /proc/irq -name "smp_affinity"`
do
var=0x`cat ${file}`
var="$(( $var & 0x7f ))"
var=`printf '%.2x' ${var}`
sudo bash -c "echo ${var} > ${file}"
done
sudo bash -c "echo 7f > /proc/irq/default_smp_affinity"
```
其中 irq 0 和 2 會出現 `bash: line 0: echo: write error: Input/output error
` 的錯誤，透過 `cat /proc/irq/2/effective_affinity` 察看可修改的 CPU 發現回傳值為 `00` ，代表此 IRQ 所有 CPU 都不可修改。
可透過 `cat /proc/interrupts` 觀察 CPU 中斷次數是否增加， [x86](https://en.wikipedia.org/wiki/Interrupt_request_(PC_architecture)) PC 架構下 IRQ 0 為 system timer ，IRQ 2 為 cascaded inerrupts 。
### 測量 user space 與 kernel space 時間
#### 在 user space 下
使用 [clock_gettime](https://linux.die.net/man/2/clock_gettime) 下的 `CLOCK_MONOTONIC` 來計時，`CLOCK_MONOTONIC` 是從特定起始時間單調遞增的計時方式。
```c
for (int i = 0; i <= offset; i++) {
lseek(fd, i, SEEK_SET);
clock_gettime(CLOCK_MONOTONIC, &t1);
read(fd, buf, 1);
clock_gettime(CLOCK_MONOTONIC, &t2);
printf("%d %lld\n", i, elapsed(&t1, &t2));
}
```
```c
static inline long long elapsed(struct timespec *t1, struct timespec *t2)
{
return (long long) (t2->tv_sec - t1->tv_sec) * 1e9 +
(long long) (t2->tv_nsec - t1->tv_nsec);
}
```
#### 在 kernel space 下
同作業說明使用 [hrtimer](https://lwn.net/Articles/167897/) 進行量測。
```c
static ktime_t kt;
static long long fib_time_proxy(long long k)
{
kt = ktime_get();
long long result = fib_sequence(k);
kt = ktime_sub(ktime_get(), kt);
return result;
}
```
為了方便開啟 kernel space 計時功能，這裡在 [sysfs](https://www.kernel.org/doc/Documentation/filesystems/sysfs.txt) 底下註冊了一個 `ktime_measure` 變數來開關計時功能。根據 [sysfs](https://www.kernel.org/doc/Documentation/filesystems/sysfs.txt) 文件
> sysfs is a ram-based filesystem initially based on ramfs. It provides
a means to export kernel data structures, their attributes, and the
linkages between them to userspace.
[sysfs](https://www.kernel.org/doc/Documentation/filesystems/sysfs.txt) 提供了一組介面可以讓 userspace 連接到 kernelspace 的 attributes ，因此可以以如下的結構註冊一個 device attribute
```c
struct device_attribute {
struct attribute attr;
ssize_t (*show)(struct device *dev, struct device_attribute *attr,
char *buf);
ssize_t (*store)(struct device *dev, struct device_attribute *attr,
const char *buf, size_t count);
};
```
`show` function pointer 會在讀取操作時執行， `store` 會在寫入操作時執行，相關實作如下
```c
static bool ktime_enable = false;
static long long (*fib_exec)(long long);
static ssize_t ktime_measure_show(struct device *dev,
struct device_attribute *attr,
char *buf)
{
return scnprintf(buf, 2, "%d", ktime_enable);
}
static ssize_t ktime_measure_store(struct device *dev,
struct device_attribute *attr,
const char *buf,
size_t count)
{
ktime_enable = (bool) (buf[0] - '0');
fib_exec = ktime_enable ? &fib_time_proxy : &fib_sequence;
return count;
}
static DEVICE_ATTR(ktime_measure,
S_IWUSR | S_IRUGO,
ktime_measure_show,
ktime_measure_store);
```
透過寫入 `"0"` 或 `"1"` (須注意傳入的是字串）來註冊 `fib_exec` function pointer 來選擇是否執行要量測時間的 fibonacci 函式， 其中 `S_IWUSR` 表示指定的 user 才能寫入，`S_IRUGO` 則是 user, group, other 都可讀取。註冊完後在 `/sys/class/fibonacci/fibonacci` 路徑下可看到 `ktime_measure` ，可透過以下命令讀取和寫入(寫入須為 root 權限)
```shell
$ cat /sys/class/fibonacci/fibonacci/ktime_measure
$ sudo bash -c "echo 1 > /sys/class/fibonacci/fibonacci/ktime_measure"
```
然後透過 `write` 函式取出量測到的時間
```c
static ssize_t fib_write(struct file *file,
const char *buf,
size_t size,
loff_t *offset)
{
return ktime_to_ns(kt);
}
```
因此可透過以下測試程式測量時間
```c
#define KTIME_ENABLE "/sys/class/fibonacci/fibonacci/ktime_measure"
int fd_kt = open(KTIME_ENABLE, O_RDWR);
if (fd_kt < 0) {
perror("Failed to open sysfs");
err = 1;
goto close_fib;
}
write(fd_kt, "1", 2);
for (int i = 0; i <= offset; i++) {
lseek(fd, i, SEEK_SET);
read(fd, buf, 1);
ssize_t kt = write(fd, NULL, 0);
printf("%d %ld\n", i, kt);
}
...
```
:::info
註：在 [sysfs](https://www.kernel.org/doc/Documentation/filesystems/sysfs.txt) 文件中有以下描述，
> static ssize_t show_name(struct device *dev, struct device_attribute *attr,
> char *buf)
> {
> return scnprintf(buf, PAGE_SIZE, "%s\n", dev->name);
> }
>
> static ssize_t store_name(struct device *dev, struct device_attribute *attr,
> const char *buf, size_t count)
> {
> snprintf(dev->name, sizeof(dev->name), "%.*s",
> (int)min(count, sizeof(dev->name) - 1), buf);
> return count;
> }
>
> static DEVICE_ATTR(name, S_IRUGO, show_name, store_name);
>
>
> (Note that the real implementation doesn't allow userspace to set the
> name for a device.)
但實際上 `struct device` 並沒有 `name` 這個 member 可用 (大約在 Linux v2.6 之後就被刪掉了），雖然以下補充有說不應該在 userspace 修改 device `name` ，但顯然這是很過時的範例。
> TODO: 準備提交文件修改到 LKML
> :notes: jserv
:::
### Page faults
經過以上設定可以量測到下圖的結果，可以看到在 userspace 前幾次所花的時間特別高，參考 [KYG-yaya573142](https://hackmd.io/@KYWeng/rkGdultSU#%E9%87%8F%E6%B8%AC%E6%99%82%E9%96%93%E7%9A%84%E6%96%B9%E5%BC%8F) 的筆記中老師的提示，可知道可能原因是發生 page faults 導致。
參考 [Threaded RT-application with memory locking and stack handling example
](https://rt.wiki.kernel.org/index.php/Threaded_RT-application_with_memory_locking_and_stack_handling_example) ，可透過 [mlockall](https://linux.die.net/man/2/mlockall) 來將 process 的 address space 鎖在記憶體中，**儘可能**避免發生 page faults 。
根據 [mlockall](https://linux.die.net/man/2/mlockall) man page
>Real-time processes that are using mlockall() to prevent delays on page faults should reserve enough locked stack pages before entering the time-critical section, so that no page fault can be caused by function calls. This can be achieved by calling a function that allocates a sufficiently large automatic variable (an array) and writes to the memory occupied by this array in order to touch these stack pages. This way, enough pages will be mapped for the stack and can be locked into RAM. The dummy writes ensure that not even copy-on-write page faults can occur in the critical section.
可知可以先透過配置足夠大的記憶體來 touch 到足夠的空間，來使空間被 lock 在 RAM 中作為 stack pages 使用來執行對時間要求較高的程式。
接著討論[連結](https://rt.wiki.kernel.org/index.php/Threaded_RT-application_with_memory_locking_and_stack_handling_example)中的程式碼
```c
static void configure_malloc_behavior(void)
{
/* Now lock all current and future pages
from preventing of being paged */
if (mlockall(MCL_CURRENT | MCL_FUTURE))
perror("mlockall failed:");
/* Turn off malloc trimming.*/
mallopt(M_TRIM_THRESHOLD, -1);
/* Turn off mmap usage. */
mallopt(M_MMAP_MAX, 0);
}
```
[mallopt](https://man7.org/linux/man-pages/man3/mallopt.3.html) 中說明 `M_TRIM_THRESHOLD` 選項是在 `free` 後的 memory size 足夠大時會釋放記憶體資源給系統，設 `-1` 來關閉此功能。`M_MMAP_MAX` 設 `0` 可以避免 **large** allocation requests 使用 `mmap` （看起來暗示著小的記憶體分配還是可能用到 `mmap` )
做完初始設定後，透過 `malloc` 將 address space 載入記憶體中，可逾期此時會發生 page faults ，之後就算 `free(buffer)` 也不會將記憶體 swap out。可以透過 [getrusage](https://man7.org/linux/man-pages/man2/getrusage.2.html) 中的 `ru_majflt` 和 `ru_minflt` 來觀察是否方生 major page faults 與 minor page faults。
```c
static void reserve_process_memory(int size)
{
int i;
char *buffer;
buffer = malloc(size);
/* Touch each page in this piece of memory to get it mapped into RAM */
for (i = 0; i < size; i += sysconf(_SC_PAGESIZE)) {
/* Each write to this buffer will generate a pagefault.
Once the pagefault is handled a page will be locked in
memory and never given back to the system. */
buffer[i] = 0;
}
/* buffer will now be released. As Glibc is configured such that it
never gives back memory to the kernel, the memory allocated above is
locked for this process. All malloc() and new() calls come from
the memory pool reserved and locked above. Issuing free() and
delete() does NOT make this locking undone. So, with this locking
mechanism we can build C++ applications that will never run into
a major/minor pagefault, even with swapping enabled. */
free(buffer);
}
```
但 `mlockall` 不是萬能的，根據 [stackoverflow](https://stackoverflow.com/a/24026447/11425048) 的討論，至少包含以下兩個項目都可能造成 page fault
* [memory compaction ](https://lwn.net/Articles/368869/) 會需要 migrate page frame ，而[migrating mlocked pages](https://www.kernel.org/doc/Documentation/vm/unevictable-lru.txt) 就提到 migrate page 可能會與 `mlock` 競爭造成 `page fault`
* [NUMA scheduling progress](https://lwn.net/Articles/568870/) 提到因 [NUMA balacing](https://documentation.suse.com/sles/15-SP1/html/SLES-all/cha-tuning-numactl.html) 功能，scheduler 會定期去掃描 process's address space 並撤銷 page 在 RAM 的使用權，導致下次 process 存取記憶體時發生 page fault。
而我在作業中實際測試，即便 `mlockall` 之後第一次執行 `clock_gettime` 也確實會發生 minor page fault ，這篇 [stackoverflow](https://stackoverflow.com/a/53366400/11425048) 討論也有提到這件事，但沒明確說明原因，因此目前作法先在測試前先 call `clock_gettime` 一次，而量測結果如下，第一筆的量測時間確實減少了。
#### 探討 `clock_gettime` 造成 page fault 的原因
首先透過 [perf trace](https://man7.org/linux/man-pages/man1/perf-trace.1.html) 追蹤發生 page fault 的程式碼
```shell
$ sudo perf trace -F all --call-graph=dwarf ./utime
11.651 ( 0.000 ms): utime/24527 minfault [__vdso_clock_gettime+0x46] => [vvar]@0x80 (d.)
__vdso_clock_gettime ([vdso])
__GI___clock_gettime (inlined)
main (/home/bakud/linux2021/fibdrv/utime)
__libc_start_main (/usr/lib/x86_64-linux-gnu/libc-2.31.so)
_start (/home/bakud/linux2021/fibdrv/utime)
```
可以看到在執行 `__vdso_clock_gettime` 時發生了 minor page fault ，並且呼叫了 INT 0x80 system call ，接著追蹤 Linux kernel 察看原始碼的實作
```c
/* /arch/x86/entry/vdso/vclock_gettime.c */
extern int __vdso_clock_gettime(clockid_t clock, struct __kernel_timespec *ts);
int __vdso_clock_gettime(clockid_t clock, struct __kernel_timespec *ts)
{
return __cvdso_clock_gettime(clock, ts);
}
int clock_gettime(clockid_t, struct __kernel_timespec *)
__attribute__((weak, alias("__vdso_clock_gettime")));
```
gcc extension 中的 [weak](https://gcc.gnu.org/onlinedocs/gcc-4.3.5/gcc/Function-Attributes.html) 可以使函式在編譯時產生 weak symbol ，讓程式可以只宣告函式不實作，卻不會發生編譯錯誤。可以預期 `__vdso_clock_gettime` 實作會被編譯到一個動態連結函式庫 `.so` 裡，等待執行時期連結。
可以透過 [`ldd`](https://man7.org/linux/man-pages/man1/ldd.1.html) 察看執行檔會連結到的 share libraires
```shell
$ ldd ./utime
linux-vdso.so.1 (0x00007ffcfb36d000)
libc.so.6 => /lib/x86_64-linux-gnu/libc.so.6 (0x00007f2bd5efe000)
/lib64/ld-linux-x86-64.so.2 (0x00007f2bd6109000)
```
```c
/* /lib/vdso/gettimeofday.c */
static __maybe_unused int
__cvdso_clock_gettime(clockid_t clock, struct __kernel_timespec *ts)
{
return __cvdso_clock_gettime_data(__arch_get_vdso_data(), clock, ts);
}
```
剛開始不理解這裡為什麼會發生 page fault ，於是我去閱讀了以下參考資料。
* [Linux 核心設計: 賦予應用程式生命的系統呼叫 [1]](https://hackmd.io/@sysprog/linux-syscall?type=view#Linux-%E6%A0%B8%E5%BF%83%E8%A8%AD%E8%A8%88-%E8%B3%A6%E4%BA%88%E6%87%89%E7%94%A8%E7%A8%8B%E5%BC%8F%E7%94%9F%E5%91%BD%E7%9A%84%E7%B3%BB%E7%B5%B1%E5%91%BC%E5%8F%AB)
* [什麼是 Linux vDSO 與 vsyscall？——發展過程 [2]](https://alittleresearcher.blogspot.com/2017/04/linux-vdso-and-vsyscall-history.html)
首先得先提到 vsyscall ，因一般的 system call 都會涉及到[特權級](https://en.wikipedia.org/wiki/Protection_ring#Privilege_level)轉換，造成很大的 overhead，因此 kernel 作者便想出一種機制來改善 ，vsyscall 是把 kernel mapping 到 user space 的一種機制，這種機制可以讓頻繁執行的 system call 先 mapping 到 user space ，直接存取 user space address 可以使效能獲得大幅提升，然而初期的 `vsyscall` 會把 kernel 程式 mapping 到一個固定的記憶體位置，這會造成很大的資安漏洞。
後來 [ASLR（address space layout randomization）](https://en.wikipedia.org/wiki/Address_space_layout_randomization#Linux) 被引入 Linux kernel 中，而 [vDSO](https://man7.org/linux/man-pages/man7/vdso.7.html) (virtual dynamic shared object) 的實作便參照了 ASLR 的觀念，改善了 `vsyscall` 的缺點。顧名思義， vDSO 會讓 kernel 被 mapping 到隨機的記憶體位置，而這段資料在 kernel space 為可讀可寫，在 user space 僅為可讀。
接著來驗證 ASLR，參考自 [askUbuntu](https://askubuntu.com/a/318476) ，可以透過以下指令察看目前狀態
```shell
$ cat /proc/sys/kernel/randomize_va_space
```
其中 `0` 表示關閉 randomization， `1` 表示僅在S hared libraries, stack, mmap, VDSO, heap 為 randomized， `2` 表示 Full randomization ，可以透過 `echo 2 | sudo tee /proc/sys/kernel/randomize_va_space` 打開，或是在 `/etc/sysctl.d/01-disable-aslr.conf` 中加入 `kernel.randomize_va_space = 2`
開啟後可以輸入以下指令察看，會發現每次印出來的地址都會不一樣。
```shell
$ cat /proc/self/maps
555f68288000-555f6828a000 r--p 00000000 08:13 1311101 /usr/bin/cat
555f6828a000-555f6828f000 r-xp 00002000 08:13 1311101 /usr/bin/cat
555f6828f000-555f68292000 r--p 00007000 08:13 1311101 /usr/bin/cat
555f68292000-555f68293000 r--p 00009000 08:13 1311101 /usr/bin/cat
555f68293000-555f68294000 rw-p 0000a000 08:13 1311101 /usr/bin/cat
555f68294000-555f682b5000 rw-p 00000000 00:00 0 [heap]
7f1834c19000-7f1834c3b000 rw-p 00000000 00:00 0
7f1834c3b000-7f1835a65000 r--p 00000000 08:13 1313911 /usr/lib/locale/locale-archive
7f1835a65000-7f1835a8a000 r--p 00000000 08:13 1315717 /usr/lib/x86_64-linux-gnu/libc-2.31.so
7f1835a8a000-7f1835c02000 r-xp 00025000 08:13 1315717 /usr/lib/x86_64-linux-gnu/libc-2.31.so
7f1835c02000-7f1835c4c000 r--p 0019d000 08:13 1315717 /usr/lib/x86_64-linux-gnu/libc-2.31.so
7f1835c4c000-7f1835c4d000 ---p 001e7000 08:13 1315717 /usr/lib/x86_64-linux-gnu/libc-2.31.so
7f1835c4d000-7f1835c50000 r--p 001e7000 08:13 1315717 /usr/lib/x86_64-linux-gnu/libc-2.31.so
7f1835c50000-7f1835c53000 rw-p 001ea000 08:13 1315717 /usr/lib/x86_64-linux-gnu/libc-2.31.so
7f1835c53000-7f1835c59000 rw-p 00000000 00:00 0
7f1835c6b000-7f1835c6c000 r--p 00000000 08:13 1315501 /usr/lib/x86_64-linux-gnu/ld-2.31.so
7f1835c6c000-7f1835c8f000 r-xp 00001000 08:13 1315501 /usr/lib/x86_64-linux-gnu/ld-2.31.so
7f1835c8f000-7f1835c97000 r--p 00024000 08:13 1315501 /usr/lib/x86_64-linux-gnu/ld-2.31.so
7f1835c98000-7f1835c99000 r--p 0002c000 08:13 1315501 /usr/lib/x86_64-linux-gnu/ld-2.31.so
7f1835c99000-7f1835c9a000 rw-p 0002d000 08:13 1315501 /usr/lib/x86_64-linux-gnu/ld-2.31.so
7f1835c9a000-7f1835c9b000 rw-p 00000000 00:00 0
7ffe14809000-7ffe1482a000 rw-p 00000000 00:00 0 [stack]
7ffe14932000-7ffe14936000 r--p 00000000 00:00 0 [vvar]
7ffe14936000-7ffe14938000 r-xp 00000000 00:00 0 [vdso]
ffffffffff600000-ffffffffff601000 --xp 00000000 00:00 0 [vsyscall]
```
接著繼續看到參考資料 [[2]](https://alittleresearcher.blogspot.com/2017/04/linux-vdso-and-vsyscall-history.html) ，其中提到 vsyscall 有三種核心參數可調整，分別是
> * none：完全取消 vsyscall 頁面，如果有程式存取該頁面，就直接引發記憶體區段錯誤（segmentation fault），最終導致程式關閉。
> * native：保留 vsyscall 頁面，但其內容不再是特製的 vsyscall 函式，而是單純發起正常的系統呼叫。
> * emulate：vsyscall 頁面的內容與 native 相同，但存取權限調整為唯讀且不可執行，當程式存取時，會觸發頁面錯誤（page fault），但核心會對觸發錯誤的位址進行篩選，如果是因為呼叫 vsyscall 頁面而導致，就在核心中模擬對應的 vsyscall 函式，然後正常返回。
可以透過 `cat /usr/src/linux-headers-$(uname -r)/.config | grep VSYSCALL
` 察看是否開啟，我的回傳 `CONFIG_X86_VSYSCALL_EMULATION=y` ，為 emulate 模式，因此當程式首次存取到 `__vdso_data` 時就會引發 page fault !
:::info
註： vsyscall 與 vDSO 相關的程式碼至今仍然變動的很快，例如這個 [patch](https://lkml.org/lkml/2020/12/20/132) 提到 `randomize_va_space` 已經又新增了一種模式，可以印證老師說的 Linux kerenl 是活的，所以要每年都更新教材 :+1:
:::
### 取 95% 信賴區間樣本平均
參考[作業說明](https://hackmd.io/@sysprog/linux2021-fibdrv#%E7%94%A8%E7%B5%B1%E8%A8%88%E6%89%8B%E6%B3%95%E5%8E%BB%E9%99%A4%E6%A5%B5%E7%AB%AF%E5%80%BC)，去除兩個標準差以外的 outliers ，對剩餘樣本取平均 50 次平均計算結果，如下圖表示，程式碼可見 [commit b3fd87](https://github.com/bakudr18/fibdrv/commit/b3fd8753d40cc3e28b05a44bffa506c5f6a3623a) 。
---
## Fibonacci 實作
### Fast doubling
依[作業說明](https://hackmd.io/@sysprog/linux2021-fibdrv#-%E8%B2%BB%E6%B0%8F%E6%95%B8%E5%88%97)，可用 Fast doubling 手法計算得以下公式
$$
\begin{split}
F(2k) &= F(k)[2F(k+1) - F(k)] \\
F(2k+1) &= F(k+1)^2+F(k)^2
\end{split}
$$
參考 [Fast doubling](https://chunminchang.github.io/blog/post/calculating-fibonacci-numbers-by-fast-doubling) 文章說法，再新增一行公式帶入
$$
\begin{split}
F(2k) &= F(k)[2F(k+1) - F(k)] \\
F(2k+1) &= F(k+1)^2+F(k)^2 \\
F(2k+2) &= F(2k+1) + F(2k)
\end{split}
$$
由以上公式可知，只要用 $F(k)$ 與 $F(k+1)$ 就可計算出 $F(2k)$, $F(2k+1)$ 與 $F(2k+2)$ 的結果，因此假設要計算 $F_{10}$ ( $n$ 為偶數) 流程會如下，
$$
\require{AMScd}
\begin{CD}
\left(
\begin{array}{c}
F_0 \\
F_1
\end{array}
\right)
@>{2k+1, 2k+2}>>
\left(
\begin{array}{c}
F_1 \\
F_2
\end{array}
\right)
@>{2k, 2k+1}>>
\left(
\begin{array}{c}
F_2 \\
F_3
\end{array}
\right)
@>{2k+1, 2k+2}>>
\left(
\begin{array}{c}
F_5 \\
F_6
\end{array}
\right)
@>{2k, 2k+1}>>
\left(
\begin{array}{c}
F_{10} \\
F_{11}
\end{array}
\right)
\end{CD}
$$
而若是計算 $F_{11}$ ( $n$ 為奇數) 則為
$$
\require{AMScd}
\begin{CD}
\left(
\begin{array}{c}
F_0 \\
F_1
\end{array}
\right)
@>{2k+1, 2k+2}>>
\left(
\begin{array}{c}
F_1 \\
F_2
\end{array}
\right)
@>{2k, 2k+1}>>
\left(
\begin{array}{c}
F_2 \\
F_3
\end{array}
\right)
@>{2k+1, 2k+2}>>
\left(
\begin{array}{c}
F_5 \\
F_6
\end{array}
\right)
@>{2k+1, 2k+2}>>
\left(
\begin{array}{c}
F_{11} \\
F_{12}
\end{array}
\right)
\end{CD}
$$
可知若要計算的 $n$ 為奇數時，須先計算 $k=(n-1)/2$ 的 $F(k)$ 與 $F(k+1)$ ，若 $n$ 為偶數，則須先計算$k=n/2$ 的 $F(k)$ 與 $F(k+1)$ 。
而由 $n/2$ 到 0 的次數便可知道須計算的次數，由於 C 語言整數除法特性會把小數點後捨去，而除以 2 又相當於 right shift 1，因此可將部份計算流程改以 bitwise operation 操作，而以 iterative 手法實作結果如下
```c
static long long fib_fast_doubling(unsigned int n)
{
// The position of the highest bit of n.
// So we need to loop `h` times to get the answer.
// Example: n = (Dec)50 = (Bin)00110010, then h = 6.
// ^ 6th bit from right side
unsigned int h = 0;
for (unsigned int i = n ; i ; ++h, i >>= 1);
long long a = 0; // F(0) = 0
long long b = 1; // F(1) = 1
// There is only one `1` in the bits of `mask`. The `1`'s position is same
// as the highest bit of n(mask = 2^(h-1) at first), and it will be shifted
// right iteratively to do `AND` operation with `n` to check `n_j` is odd or
// even, where n_j is defined below.
for (unsigned int mask = 1 << (h - 1); mask; mask >>= 1) { // Run h times!
// Let j = h-i (looping from i = 1 to i = h), n_j = floor(n / 2^j) = n
// >> j (n_j = n when j = 0), k = floor(n_j / 2), then a = F(k), b =
// F(k+1) now.
uint64_t c =
a * (2 * b - a); // F(2k) = F(k) * [ 2 * F(k+1) - F(k) ]
uint64_t d = a * a + b * b; // F(2k+1) = F(k)^2 + F(k+1)^2
if (mask & n) { // n_j is odd: k = (n_j-1)/2 => n_j = 2k + 1
a = d; // F(n_j) = F(2k + 1)
b = c + d; // F(n_j + 1) = F(2k + 2) = F(2k) + F(2k + 1)
} else { // n_j is even: k = n_j/2 => n_j = 2k
a = c; // F(n_j) = F(2k)
b = d; // F(n_j + 1) = F(2k + 1)
}
}
return a;
}
```
### Most significant bit
上述提到計算 $n/2$ 到 0 次數，即是在計算 Most significant bit (MSB) 位置，根據 [quiz2](https://hackmd.io/@MR-9Qa9wQLWglSAUyd6UhQ/SkS-Y_lX_#Linux-kernel-%E4%B8%AD%E7%9A%84-power-of-two-%E5%AF%A6%E4%BD%9C) 寫過的共筆，在 Linux kernel 的 [bitops.h](https://github.com/torvalds/linux/blob/v5.8/include/linux/bitops.h) 有提供了 `fls` 相關函式來計算 MSB ，內部是以 CPU 指令實作優化的，因此可將上述程式碼替換為
```diff
- unsigned int h = 0;
- for (unsigned int i = n ; i ; ++h, i >>= 1);
+ unsigned int h = fls(n);
```
接著比較 sequence, fast doubling 與 fast doubling with fls 的結果如下，可以發現兩個 fast doubling 明顯比 sequence 方法快，而加了 `fls` 優化大約快了 8 %
### 加入 clz 比較
[`__builtin_clz`](https://gcc.gnu.org/onlinedocs/gcc/Other-Builtins.html) 是 GCC extension 用於計算 leading zero 的數量，因此新增以下修改
```diff
- unsigned int h = fls(n);
+ unsigned int h = 32 - __builtin_clz(n);
```
實驗結果如下，顯示使用 `__builtin_clz` 並沒有比 `fls` 快。
---
## 大數運算
> 本節以 [KYG-yaya573142](https://hackmd.io/@KYWeng/rkGdultSU#%E5%A6%82%E4%BD%95%E6%B8%9B%E5%B0%91%E5%A4%A7%E6%95%B8%E9%81%8B%E7%AE%97%E7%9A%84%E6%88%90%E6%9C%AC) 的實作為基礎進行延伸並與其比較，量測皆為在 user space 進行，實作程式碼在 [bakudr18/bignum](https://github.com/bakudr18/bignum)。
### `bn_lshift` 改進
首先參考 [bignum](https://github.com/sysprog21/bignum) 以同樣方式實作如下
```c
static bn_digit bn_lshift_mod(const bn_digit *src,
uint64_t size,
unsigned int shift,
bn_digit *dst)
{
if (!size)
return 0;
shift &= BN_DIGIT_BITS - 1;
if (!shift) {
memcpy(dst, src, BN_DIGIT_SIZE * size);
return 0;
}
const unsigned int subp = BN_DIGIT_BITS - shift;
bn_digit carry = 0;
do {
const bn_digit p = *src++;
*dst++ = (p << shift) | carry;
carry = p >> subp;
} while (--size);
return carry;
}
void bn_lshift_original(const bn *src, unsigned int bits, bn *dst)
{
if (unlikely(bits == 0)) {
bn_cpy(dst, src);
return;
}
const unsigned int digits = bits / BN_DIGIT_BITS;
bits %= BN_DIGIT_BITS;
bn_digit cy;
bn_resize(dst, src->size + digits);
cy = bn_lshift_mod(src->num, src->size, bits, dst->num + digits);
bn_digit_zero(dst->num, digits);
if (cy) {
bn_resize(dst, dst->size + 1);
dst->num[dst->size - 1] = cy;
}
}
```
原實作方式是由 lower byte 開始移動到 higher byte ，而在 MSB 需要判斷是否須多配置 1 byte 來儲存 `cy` ，這會造成潛在的 `realloc` 成本。
新的實作方式如下，對需要 `realloc` 的 size 是可以事先計算的，也因此可以從 higher byte 開始移動，避免需要 `realloc` 兩次。另外針對 8 bits 倍數的 left shift 同樣可用 [quiz2](https://hackmd.io/@MR-9Qa9wQLWglSAUyd6UhQ/SkS-Y_lX_#%E6%94%B9%E5%96%84-bitcpy-%E6%95%88%E8%83%BD) 中提到的 `memcpy` 來優化，相比原實作只在 64 bits 倍數使用 `memcpy` ，新實作的 `memcpy` 使用機會更高。
```c
void bn_lshift(const bn *src, unsigned int bits, bn *dst)
{
if (unlikely(bits == 0)) {
bn_cpy(dst, src);
return;
}
const unsigned int digits = bits / BN_DIGIT_BITS;
bits %= BN_DIGIT_BITS;
const unsigned int subp = BN_DIGIT_BITS - bits;
uint64_t orig_size = src->size;
uint64_t new_size = src->size + digits;
if (bits > 0 && src->num[src->size - 1] >> subp != 0)
new_size++;
bn_resize(dst, new_size);
if ((bits & 7) == 0) {
dst->num[new_size - 1] = 0; /* clean uninitialized block */
uint8_t *dst_head = (uint8_t *) (dst->num + digits) + (bits >> 3);
memcpy(dst_head, src->num, BN_DIGIT_SIZE * src->size);
memset(dst->num, 0, BN_DIGIT_SIZE * digits + (bits >> 3));
return;
}
bn_digit *src_tail = &src->num[src->size - 1];
bn_digit *dst_tail = &dst->num[dst->size - 1];
if (src->size + digits + 1 == dst->size) {
*dst_tail-- = *src_tail >> subp;
}
while (--orig_size) {
const bn_digit p = *src_tail--;
*dst_tail-- = (p << bits) | (*src_tail >> subp);
}
*dst_tail = *src_tail << bits;
bn_digit_zero(dst->num, digits);
}
```
針對未對齊 8 bits 的 shift 做效能分析結果如下，可發現在 bignum size 愈大情況下新實作的執行時間改善愈多。
而針對 8 bits 為倍數的 left shift input 做效能分析，執行時間縮短情況明顯更佳。
### `bn_mul_partial` 嘗試優化
`bn_mul_partial` 是實作中用於計算 `c[size] += a[size] * k`, 並回傳進位值 `carry` 的函式，由 [perf-record](https://man7.org/linux/man-pages/man1/perf-record.1.html) 分析 `Fib(300000)` 結果有約 75% 時間都花在這個函式上，原實作方式如下
```c
#define bn_digit_mul(u, v, hi, lo) \
__asm__("mulq %3" : "=a"(lo), "=d"(hi) : "%0"(u), "rm"(v))
bn_digit bn_mul_partial(const bn_digit *a,
uint64_t size,
const bn_digit k,
bn_digit *c)
{
if (k <= 1) {
if (k == 0)
bn_digit_zero(c, size);
else
bn_digit_cpy(a, size, c);
return 0;
}
bn_digit carry = 0;
for (uint64_t i = 0; i < size; i++) {
bn_digit high, low;
bn_digit_mul(a[i], k, high, low);
carry = high + ((low += carry) < carry);
carry += ((c[i] += low) < low);
}
return carry;
}
```
原實作已經使用內嵌組合語言指令 `mulq` 來做乘法優化，此指令為計算 `u * v` 的結果將低位元組存在 `lo` ，並將乘法進位後產生的高位元組存在 `hi`。
原實作技巧已經使用到組合語言指令，比較難在此繼續優化，因此這裡轉念想針對 power of 2 的 `k` 直接做 left shift 做運算，實作方式如下
```c
if (is_power_of_2(k)) {
int bits = ilog2(k);
const int subp = BN_DIGIT_BITS - bits;
for (uint64_t i = 0; i < size; i++) {
bn_digit p = a[i];
c[i] += (p << bits) | carry;
carry = p >> subp;
}
return carry;
}
```
以 `k = 1ULL << 63` 為測資對 `size = 1~100000` 執行效能測試結果如下，反而原作法執行時間更短，可見得 `mulq` 指令的優化效果是非常驚人的。
:::warning
考慮以下程式碼:
```cpp
static inline long long multiply(long long a, long long b)
{
long long res = 0;
/*
* 1. Find the MSB (Most Significant Bit) position of 'b'.
* 2. Accumulate the result with (1 << MSB position).
* 3. Clear the MSB position.
*
* FIXME: Implement the negative number multiplication.
*/
while (b) {
int msb_pos = ilog2(b);
res += (a << msb_pos);
b &= ~(1ULL << msb_pos); /* Clear MSB */
}
return res;
}
```
依據 [Intel® Xeon® Scalable Processor Throughput and Latency](https://software.intel.com/en-us/download/intel-xeon-scalable-processor-throughput-and-latency)，可知 `imul rax rcx` 只需花費 3 個 CPU clock cycles (在 Intel Xeon 處理器，乘法指令只需花費 3-8 CPU clock cycles)。`multiply()` 程式碼 `b` 的位元有效位越多的話，則會花更多的時間運算。所以，使用 clz/fls, bitwise 操作與加法實作整數乘法 (避開乘法操作)，在 Xeon 處理器 (或較新的 Intel 消費型處理器) 上，反而需要花費更多的運算。
:notes: jserv
:::
### `bn_resize` 改進
`bn_resize` 內部使用到 `realloc` 函式來配置新的記憶體空間給加法或乘法等運算須增加的空間，為了減少實際 `realloc` 次數，這裡參考了 [quiz3](https://hackmd.io/@sysprog/linux2021-quiz3#%E6%B8%AC%E9%A9%97-1) 的 `xs_string` 結構，使用到 [bit-field](https://hackmd.io/@sysprog/c-bitfield?type=view) 與 `capacity` 的技巧。使用結構如下
```c
struct bn {
uint64_t size : 57, /* count of num. bytes length = count * BN_DIGIT_SIZE */
capacity : 6, /* actual bytes allocated is (2^capacity) * BN_DIGIT_SIZE */
sign : 1;
bn_digit *num;
};
```
以此結構最大可儲存 $2^{57}-1$ 個 `bn_digit` ，`capacity` 最大值為 63 ，表示最大會配置到 $2^{63}-1$ 個 `bn_digit` 大小的記憶體，但實際只能用到 $2^{57}-1$ ，另外保留 1 bit 紀錄正負號。
依照此結構實作 `bn_resize` 如下
```c
static void bn_resize(bn *src, uint64_t size)
{
int capacity = ilog2(size) + 1;
if (src->capacity >= capacity) {
src->size = size;
return;
}
/* Note that the added memory is not initialized */
bn_digit *tmp = REALLOC(src->num, BN_DIGIT_SIZE * (1ULL << capacity));
src->num = tmp;
src->size = size;
src->capacity = capacity;
}
```
只有在 `capacity` 變大時才會 `realloc` ，可大幅降低實際 `realloc` 的次數，將此實作與 [KYG-yaya573142](https://hackmd.io/@KYWeng/rkGdultSU#%E5%A6%82%E4%BD%95%E6%B8%9B%E5%B0%91%E5%A4%A7%E6%95%B8%E9%81%8B%E7%AE%97%E7%9A%84%E6%88%90%E6%9C%AC) 的實作進行效能比較，在 `-O2` 編譯優化後的執行結果如下，執行時間大幅降低，可見得重複 `realloc` 的時間成本相當顯著。
當然預先配置記憶體的 trade-off 就是會造成記憶體的浪費，使用 [Valgrind massif](https://valgrind.org/docs/manual/ms-manual.html) 紀錄 `F(1000000)` 執行時期所使用的記憶體用量，首先是 [KYG-yaya573142](https://hackmd.io/@KYWeng/rkGdultSU#%E5%A6%82%E4%BD%95%E6%B8%9B%E5%B0%91%E5%A4%A7%E6%95%B8%E9%81%8B%E7%AE%97%E7%9A%84%E6%88%90%E6%9C%AC) 的實作，記憶體使用的峰值約為 462 KB
```shell
KB
462.4^ #
| @#
| @ @#
| @ :: @#
| @ @ :::@#
| :@ ::@ ::::@#
| :::@ :::: @ : ::::@#
| :: : :@ : ::::: @::: ::::::@#
| : :: :@ :: ::::: @:::::::::::@#
| :: :::: :: :@::::@::::::: @:::: ::::::@#
| :: : ::: : :: :@::::@: ::::: @:::: ::::::@#
| @:::::: :::::: : :: :@::::@: ::::: @:::: ::::::@#
| : @: : :: :::::: : :: :@::::@: ::::: @:::: ::::::@#
| : :@: : :: :::::: : :: :@::::@: ::::: @:::: ::::::@#
| ::::::::::@: : :: :::::: : :: :@::::@: ::::: @:::: ::::::@#
| @@ :: : :::: :@: : :: :::::: : :: :@::::@: ::::: @:::: ::::::@#
| @ ::@ ::: : :::: :@: : :: :::::: : :: :@::::@: ::::: @:::: ::::::@#
| @::::@ ::: : :::: :@: : :: :::::: : :: :@::::@: ::::: @:::: ::::::@#
| :@: ::@ ::: : :::: :@: : :: :::::: : :: :@::::@: ::::: @:::: ::::::@#
| :::@: ::@ ::: : :::: :@: : :: :::::: : :: :@::::@: ::::: @:::: ::::::@#
0 +----------------------------------------------------------------------->Mi
0 642.8
```
接著是我的實作，記憶體峰值約來到 611 KB
```shell
KB
611.0^ #
| @#
| :@ @#
| @@ :@:@#
| @ @@ :: ::@:@#
| @ @@:@@ :::::@:::@:@#
| @ @:@ : :@@:@@::::::@:::@:@#
| :::@:@:@: :: :::@@:@@::::::@:::@:@#
| @:: : :@:@:@:::::::::@@:@@::::::@:::@:@#
| :::@: :: ::: :@:@:@:::::::::@@:@@::::::@:::@:@#
| @: :@: : :::: :@:@:@:::::::::@@:@@::::::@:::@:@#
| @@: :@: : :::: :@:@:@:::::::::@@:@@::::::@:::@:@#
| @@ :@@: :@: : :::: :@:@:@:::::::::@@:@@::::::@:::@:@#
| @@ :@ :::@@: :@: : :::: :@:@:@:::::::::@@:@@::::::@:::@:@#
| @ @@ :::@ :::@@: :@: : :::: :@:@:@:::::::::@@:@@::::::@:::@:@#
| ::@:@@ :::@ :::@@: :@: : :::: :@:@:@:::::::::@@:@@::::::@:::@:@#
| @: @:@@ :::@ :::@@: :@: : :::: :@:@:@:::::::::@@:@@::::::@:::@:@#
| ::@: @:@@ :::@ :::@@: :@: : :::: :@:@:@:::::::::@@:@@::::::@:::@:@#
| :::::@: @:@@ :::@ :::@@: :@: : :::: :@:@:@:::::::::@@:@@::::::@:::@:@#
|:::: ::@: @:@@ :::@ :::@@: :@: : :::: :@:@:@:::::::::@@:@@::::::@:::@:@#
0 +----------------------------------------------------------------------->Mi
0 173.1
```
明顯多了約 149 KB 的記憶體用量，也就是多使用了約 24.3% 的記憶體，因此還是得根據使用場合決定適合的 `realloc` 方式。
## Mutex 對 Linux 核心模組的影響
原始程式碼會在 open character device 時會嘗試 lock mutex，在 relese 時釋放，因此可避免有多個 thread 同時使用 character device。在前述章節中提到我在 [sysfs](https://www.kernel.org/doc/Documentation/filesystems/sysfs.txt) 中註冊了一個 `ktime_enable` 的 static variable 用來開啟或關閉 ktime measurement，而 `ktime_enable` 就屬於共用資源。
```c
static bool ktime_enable = false;
static ssize_t fib_write(struct file *file,
const char *buf,
size_t size,
loff_t *offset)
{
return ktime_enable ? ktime_to_ns(kt) : 1;
}
```
考慮以下程式碼，開啟兩個 pthread ，`fib_with_ktime` 在每次執行 `read`, `write` 前都會先開啟 `ktime_enable`， 而 `fib_without_ktime` 則是相反，在關閉 mutex 支援的情況下會有 race condition 發生，`write` 回傳結果很可能會不如預期，如 `fib_with_ktime` 的 thread 量測的 time elapsed 回傳為 1 ns。
```c
#define FIB_DEV "/dev/fibonacci"
#define KTIME_ENABLE "/sys/class/fibonacci/fibonacci/ktime_measure"
void *fib_with_ktime(void *arg)
{
char buf[1];
int offset = 92;
int fd = open(FIB_DEV, O_RDWR);
if (fd < 0) {
perror("Failed to open character device");
return NULL;
}
int fd_kt = open(KTIME_ENABLE, O_RDWR);
if (fd_kt < 0) {
perror("Failed to open sysfs");
close(fd);
return NULL;
}
for (int i = 0; i <= offset; i++) {
write(fd_kt, "1", 2); /* Enable ktime measurement */
lseek(fd, i, SEEK_SET);
read(fd, buf, 1);
ssize_t kt = write(fd, NULL, 0);
printf("With ktime offset: %d, time: %ld ns\n", i, kt);
}
close(fd_kt);
close(fd);
return NULL;
}
void *fib_without_ktime(void *arg)
{
char buf[1];
int offset = 92;
int fd = open(FIB_DEV, O_RDWR);
if (fd < 0) {
perror("Failed to open character device");
return NULL;
}
int fd_kt = open(KTIME_ENABLE, O_RDWR);
if (fd_kt < 0) {
perror("Failed to open sysfs");
close(fd);
return NULL;
}
for (int i = 0; i <= offset; i++) {
write(fd_kt, "0", 2); /* Disable ktime measurement */
lseek(fd, i, SEEK_SET);
read(fd, buf, 1);
ssize_t kt = write(fd, NULL, 0);
printf("Without ktime offset: %d, time: %ld ns\n", i, kt);
}
close(fd_kt);
close(fd);
return NULL;
}
int main(int argc, char **argv)
{
pthread_t enable_ktime_pthread, disable_ktime_pthread;
pthread_create(&enable_ktime_pthread, NULL, &fib_with_ktime, NULL);
pthread_create(&disable_ktime_pthread, NULL, &fib_without_ktime, NULL);
pthread_join(enable_ktime_pthread, NULL);
pthread_join(disable_ktime_pthread, NULL);
return 0;
}
```
在 `With ktime offset: 33` 的結果明顯是錯誤的。
```shell
...
Without ktime offset: 37, time: 1 ns
With ktime offset: 33, time: 1 ns
With ktime offset: 34, time: 118 ns
Without ktime offset: 38, time: 1 ns
With ktime offset: 35, time: 119 ns
...
```
`ktime_enable` 還算是好發現的錯誤，在量測時間所用到的 `kt` 也是共用資源，在多執行緒下就不能保證以下程式碼所計算的 `kt` 為正確值，很有可能相減到其他 thread 計時的 `kt` ，造成量測時間錯誤，且較難以察覺！
```c
static ktime_t kt;
static long long fib_time_proxy(unsigned int k)
{
kt = ktime_get();
long long result = fib_impl(k);
kt = ktime_sub(ktime_get(), kt);
return result;
}
```
Sign in with Wallet
Connect another wallet