contribute by < Shiritai
>
$ 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: 39 bits physical, 48 bits virtual
Byte Order: Little Endian
CPU(s): 12
On-line CPU(s) list: 0-11
Vendor ID: GenuineIntel
Model name: Intel(R) Core(TM) i7-8750H CPU @ 2.20GHz
CPU family: 6
Model: 158
Thread(s) per core: 2
Core(s) per socket: 6
Socket(s): 1
Stepping: 10
CPU max MHz: 4100.0000
CPU min MHz: 800.0000
BogoMIPS: 4399.99
如果純粹將 reposiroty clone 好後使用 VSCode 開啟比如 fibdrv.c
,相信會看到貼心的 IntelliSense 的警告。
解決辦法整理如下。
開始寫作業前都應按流程安裝好 linux kernel header,確認其安裝的路徑,可移動至 /usr/src
下查看。
$ ls /usr/src | grep linux-headers
linux-headers-5.19.0-32-generic
linux-headers-5.19.0-35-generic
由以下指令確認當前 Linux 核心版本,得知第二個是我們感興趣的 kernel header。
uname -r
5.19.0-35-generic
由前車之鑑 1 和 2 得知 IntelliSense 需要調整 C/C++ Extension 設定檔: c_cpp_properties.json
。
準備 Linux 設定檔模板
由 @DokiDokiPB 的建議,對此份作業來說,cStandard
的前綴應該改成 gnu,以協助 intelliSense 解析一些結構。對此我也更新腳本法的實作,詳見 commit df77974。
{
"configurations": [
{
"name": "Linux",
"browse": {
"limitSymbolsToIncludedHeaders": true,
"databaseFilename": ""
},
"intelliSenseMode": "linux-gcc-x64",
"compilerPath": "/usr/bin/gcc",
- "cStandard": "c99",
+ "cStandard": "gnu99",
"cppStandard": "gnu++17"
}
],
"version": 4
}
加入 includePath
設定
...
"name": "Linux",
"includePath": [
"${workspaceFolder}/**",
"/usr/include",
"/usr/local/include",
"/usr/src/LINUX_KERNEL_HEADER/include",
"/usr/src/LINUX_KERNEL_HEADER/arch/x86/include",
"/usr/src/LINUX_KERNEL_HEADER/arch/x86/include/generated",
"/usr/lib/gcc/x86_64-linux-gnu/GCC_VERSION/include"
],
"browse": {
...
要調整兩點。
LINUX_KERNEL_HEADER
換成自己的 kernel header 路徑名。GCC_VERSION
換成自己的 GCC 版本。基本上 includePath
中前五項都是開發 Kernel Module 的常客,當然第七項編譯器本身也是。第六項本次實作會用到,不過未來開發其他模組的話可能會引入其他的 header。
調整 defines
使 MODULE_*
等巨集被正確的識別
defines
是為了解決條件編譯的問題,當中加入的字串會使 IntelliSense 在 indexing 時看見需指定條件編譯後才會編譯到的程式碼。
...
"compilerPath": "/usr/bin/gcc",
"defines": [
"__GNUC__",
"__KERNEL__",
"MODULE"
],
"cStandard": "c99",
...
到此便完成設定,IntelliSense 復活了 :)
Reference
樹狀遞迴版
\[ fib :: \mathbb{R} \rightarrow \mathbb{R} \\ fib(0) = 0 \\ fib(1) = 1 \\ fib(n+2) = fib(n+1) + fib(n) \]
線性遞迴版
\[ fib2 :: \mathbb{R} \rightarrow (\mathbb{R}, \mathbb{R}) \\ fib2(0) = (0, 1) \\ fib2(n+1) = (y, x + y) \\ where\ (x, y) = fib2(n) \]
Binet’s Formula
\[ \dfrac{1}{\sqrt{5}} \Big( \big( \dfrac{1 + \sqrt{5}}{2} \big)^n - \big( \dfrac{1 - \sqrt{5}}{2} \big)^n \Big) \]
當中指數與根號運算造成的計算誤差和複雜度都不容小覷。費氏數列分析一文使用 GNU 的 GMP
, MPFFR
這倆透過盡情使用記憶體資源計算遞迴版 (下述分解法)和本公式解的運算時間,得到前者大勝後者的結果。
最早可追溯至 Nikolai. N. Vorobev 提出的方法。
由費氏數列的遞迴定義,也許你會猜測是否能尋找費氏數列中各項的關聯,分解法便是對此疑問的解答。令費氏數第 \(n\) 項為 \(f_n\),則
\[ f_{n+k} = f_{n-1} \cdot f_k + f_n \cdot f_{k+1} \]
上式對 \(n \ge 1, k \ge 0 \ \land n, k \in \mathbb{R}\) 成立。
證明
如果根據費氏數列定義拆解此式,可以得到如回音般的結果。
\[ f_{n+k} = f_{n-1} \cdot f_k + f_n \cdot f_{k+1} \\ = f_{n-1} \cdot f_k + f_{n} \cdot \big( f_{k} + f_{k-1} \big) \\ = \big( f_{n-1} + f_{n} \big) \cdot f_{k} + f_n + f_{k-1} \\ = f_{n+1} \cdot f_k + f_n \cdot f_{k-1} \]
這樣的結果與原本分解式的展開只差在 \(n\) 和 \(k\) 的互換。並不能給我們額外的結論,也許我們應該另尋其道。通常遇到這種難以直接證明的,我們可以用終極絕招數學歸納法:
Base case: 帶入即正確。
好像沒有寫的必要 :)
\[ f_{1} = f_{1+0} = f_{0} \cdot f_{0} + f_{1} \cdot f_{1} = f{1} = 1 \]
Inductive case
假設分解式對 \(t \le n\) 成立,則當 \(t = n + 1\) 時,利用費氏數列的定義和 \(t \le n\) 的情況,可得
\[ f_{(n+1)+k} \stackrel{\rm{def\ of\ fib}}{=} f_{n+k} + f_{(n-1)+k} \\ \stackrel{\rm{eq\ when\ }t \le n}{=} f_{n-1} \cdot f_k + f_n \cdot f_{k+1} + f_{n-2} \cdot f_k + f_{n-1} \cdot f_{k+1} \\ = \big( f_{n-1} + f_{n-2} \big) \cdot f_k + \big( f_{n} + f_{n-1} \big) \cdot f_{k+1} \\ \stackrel{\rm{def\ of\ fib}}{=} f_{n} \cdot f_k + f_{n+1} \cdot f_{k+1} \]
由數學歸納法得證。
以分解法為基礎,今若要求 \(f_{2n}\),有
\[ f_{2n} = f_{n+n} = f_{n-1} \cdot f_{n} + f_{n} \cdot f_{n+1} \\ = f_{n-1} \cdot f_{n} + f_{n} \cdot (f_{n} + f_{n-1}) \\ = f_{n-1} \cdot 2 f_{n} + f_{n} \cdot f_{n} = f_n \cdot (2f_{n-1} + f_{n}) \]
可以發現 \(f_{2n}\) 需要 \(f_{n-1}, f_{n}\)。那麼 \(f_{2n-1}\) 呢?
\[ f_{2n-1} = f_{n+(n-1)} = f_{n-1} \cdot f_{n-1} + f_{n} \cdot f_{n} = f_{n-1}^2 + f_n^2 \]
也需要 \(f_{n-1}, f_{n}\),真好。
不過其實有其他做法,同樣考慮 \(f_{2n}\),有
\[ f_{2n} \stackrel{\rm{prev\ result}}{=} f_n \cdot (2f_{n-1} + f_{n}) \\ = f_n \cdot (2f_{n+1} - f_{n}) \]
由於我們透過合併改關注 \(n+1\) 項,故改成考慮 \(2n+1\) 項 \(f_{2n+1}\),有
\[ f_{2n+1} = f_{(n+1)+n} = f_{n} \cdot f_{n} + f_{n+1} \cdot f_{n+1} = f_{n}^2 + f_{n+1}^2 \]
後者的結論是 \(f_{2n}\) 和 \(f_{2n+1}\) 都需要 \(f_n\) 和 \(f_{n+1}\),是原文中的方法,兩者的結果都很漂亮,畢竟是等價的東西。
以前者 (\(f_{2n}\), \(f_{2n-1}\) 需 \(f_{n}, f_{n-1}\)) 為例,取得二的冪就不是問題,假設我們想知道 \(f_{16}\),可以藉由以下順序求得:
digraph twos_power {
rankdir=LR;
node[shape=record];
edge[arrowsize=0.5]
12[label="<b> f2|<a> f1"]; 34[label="<b> f4|<a> f3"];
78[label="<b> f8|<a> f7"]; 16[label="f16"];
12:a -> 34:a -> 78:a -> 16;
12:b -> 34:b -> 78:b -> 16;
12:a -> 34:b -> 78:a;
12:b -> 34:a -> 78:b;
}
繼續以前者為例,那如果是任意數呢?假設要計算 \(f_{21} = f_{2n-1}\),我們需要 \(f_n = f_{11}\) 和 \(f_{n-1} = f_{10}\),為此我們還需要 \(f_6, f_5\),同樣我們又需要 \(f_3, f_2\), 最後回到 \(f_2, f_1\) 便皆大歡喜。看來無論哪個正整數,都可以使用此演算法。下圖圖像化這段流程,與上圖不同在於稍作一點簡化。
digraph flow_21 {
rankdir=LR;
node[shape=square];
1[label="f2 \n f1"];
2[label="f3 \n f2"];
3[label="f6 \n f5"];
4[label="f11 \n f10"];
5[label="f22 \n f21"];
1 -> 2 -> 3 -> 4 -> 5;
}
另外前述「前者」與「後者」的差異可透過上圖描述。前者是以數字大者 (比如 \([22, 21]\) 中的 \(22\)) 為主視角,後者則是以數字小者 (\([22, 21]\) 中的 \(21\)) 為主視角,兩者實際上會畫出一樣的關係圖。
由於實作上以後者視角比較好描述,今後將採用該視角,與作業說明同步。
以 \(f_{21}\) 為例,當初我們以 top-down 的角度分析,實作時勢必需要遞迴,如果不想要遞迴呢?我們要如何知道接下來是單純從 \((5, 6)\) 爬兩倍到 \((10, 11)\),還是額外加一,從 \((2, 3)\) 到 \((5, 6)\)?
digraph flow_21 {
rankdir=LR;
node[shape=square];
subgraph 1 {
3[label="f6 \n f5"];
4[label="f11 \n f10"];
3 -> 4;
}
subgraph 2 {
1[label="f3 \n f2"];
2[label="f6 \n f5"];
1 -> 2;
}
}
把之前的流程圖加上是否「加一」的標記:令以 \(1\) 表加一,以 \(0\) 表不加,同時將 \(f_0, f_1\) 這個 trivial 項補回來…
digraph flow_21 {
rankdir=LR;
node[shape=square];
0[label="f1 \n f0"];
1[label="f2 \n f1"];
2[label="f3 \n f2"];
3[label="f6 \n f5"];
4[label="f11 \n f10"];
5[label="f22 \n f21"];
0 -> 1 [label="1"];
1 -> 2 [label="0"];
2 -> 3 [label="1"];
3 -> 4 [label="0"];
4 -> 5 [label="1"];
}
敏銳而有看作業說明的你想必已經觀察到了,標記的 \(10101_{2}\) 不正是二進制版的 \(21_{10}\) 嗎?這顯然不是巧合,因為前述問題實際上等價於:
給定兩種操作:乘二、加一。
如何從 \(0\) 開始用最少計算次數抵達某正整數 \(n\)?
上述是在任意進位制的角度下問的問題,換成熟悉二進制與位元運算的你所熟悉的語言就是:
給定兩種操作:左移一位、加一。
如何從 \(0\) 開始用最少計算次數產生某二進制數 \(n\)?
這一問等於白問,因為答案即藏在 \(n\) 的二進制編碼 XD
不考慮大數運算的情況下,針對 long long
,為了實現 Fast doubling,首先須取得目標費波那契數的二進制編碼,接著由高至低走訪這些編碼完成乘二和加一的邏輯。
取得一個數字精確的二進制編碼,有以下兩個針對 long long
的 gcc 內建函式:
__builtin_clzll
: 其在 intel x86_64 架構下可能使用硬體指令 bit search reverse 找尋首個 1
bit,接著以 63 減去而取得領導的 0
bit 數量。__builtin_ctzll
同理在 intel x86_64 可能使用 bit search forward,最終取得末尾的 0
bit 數量。前者的使用可以加速略過領導之 0
bits 的過程,後者則由於 fast doubling 的「加一與否」牽扯到 branching,而且很難直接避免而轉換成 branchless,不過若所有位元皆為零就不需顧慮 branch 操作,故可以透過後者完成末尾 0
bits 的 branchless 乘二操作,同時降低 branching 時的比較開銷。
具體實作解釋如下。
0
bits 數量,並調整 \(f(n)\) 的 \(n\),將領導的 1
bit 上調至最高位元。
注意 c{l,t}z_long_long
為自己實作的函式,以巨集控制編譯時要採用何者。
long long fn = 0ll, // f(n)
fnp1 = 1ll; // f(n+1) n plus 1
#ifdef MY_CZ
const int lz = clz_long_long(n);
const int tz = ctz_long_long(n);
#else
const int lz = __builtin_clzll(n);
const int tz = __builtin_ctzll(n);
#endif /* MY_CZ */
n <<= lz; // so now there is no leading zeros
1
bit 的存在,即 \(n\) 是否為零。
while (n) { // exists "1" bit
long long f2n = fn * ((fnp1 << 1) - fn),
f2np1 = fn * fn + fnp1 * fnp1;
if (n & 0x8000000000000000ll) {
fn = f2np1, fnp1 = f2np1 + f2n;
} else {
fn = f2n, fnp1 = f2np1;
}
n <<= 1;
}
0
bits 所代表之「乘二」邏輯。
// speed up for tailing zeros: branchless
for (int i = 0; i < tz; ++i) {
long long f2n = fn * ((fnp1 << 1) - fn),
f2np1 = fn * fn + fnp1 * fnp1;
fn = f2n, fnp1 = f2np1;
}
return fn;
迴圈內 if…else… 的部份可用
-!!(target & (1 << i))
作為 mask 的技巧簡化成:
static inline uint64_t fast_doubling_iter(uint64_t target) { if (target <= 2) return !!target; // find first 1 uint8_t count = 63 - __builtin_clzll(target); uint64_t fib_n0 = 1, fib_n1 = 1; for (uint64_t i = count, fib_2n0, fib_2n1, mask; i-- > 0;) { fib_2n0 = fib_n0 * ((fib_n1 << 1) - fib_n0); fib_2n1 = fib_n0 * fib_n0 + fib_n1 * fib_n1; mask = -!!(target & (1UL << i)); fib_n0 = (fib_2n0 & ~mask) + (fib_2n1 & mask); fib_n1 = (fib_2n0 & mask) + fib_2n1; } return fib_n0; }
(target & (1 << i))
的意思為檢查某位元是否為 1
,前面加上 !!
表取得真或假 (0 ro 1),加負號表取得 0xFFFFFFFFFFFFFFFF
。
int absWithoutBranching(int number) {
// 對 number > 0 來說,mask = 0x00000000 (0)
// 對負數來說, mask = 0xffffffff (-1)
const int mask = number >> 31;
// 二補數的「負數取正」為「減一後位元反轉」
// 對負數來說
// 加上 mask 相當於 -1,相當於「減一」
// 對 -1 = 0xffffffff 取 XOR 相當於「位元反轉」得到負數取正
// 對正數來說
// 加上 mask = 0 等於沒做事
// 對 0 取 XOR 相當於沒做事,保持原數
return (mask + number) ^ mask;
}
最後根據此 mask 過濾是否採用 fib_2n0
, fib_2n1
與否,十分精妙。
最後我採納此技巧,將原本的 branching 調整成 branchless,不過還是保留第三階段的 branchless,以達到效能的最大提升。以下展示原 branching 調整的部分。
while (n) { // exists "1" bit
long long f2n = fn * ((fnp1 << 1) - fn), f2np1 = fn * fn + fnp1 * fnp1;
+ long long is_one = -!!(n & 0x8000000000000000ll);
+ fn = (f2np1 & is_one) + (f2n & ~is_one);
+ fnp1 = f2np1 + (f2n & is_one);
- if (n & 0x8000000000000000ll) {
- fn = f2np1, fnp1 = f2np1 + f2n;
- } else {
- fn = f2n, fnp1 = f2np1;
- }
n <<= 1;
}
雖然 __builtin_clzll
和 __builtin_ctzll
好用且高效,但本方法的使用受限於 gcc 與特定型別 (支援 unsigned int
, unsigned long
和 unsigned long long
);且參考 GCC 官方文件,當輸入為 0
時屬於 undefined behavior,可見略有不足之處。於是我受 bitwise 操作教學文件的啟發,設計巨集來生成 clz
和 ctz
函式,支援所有長度在 64 以內的型別。
原本教學文件中的程式碼如下。
int clz (int x){
if (x == 0) return 32;
int n = 0;
if ((x & 0xFFFF0000) == 0) {n += 16; x <<= 16;}
if ((x & 0xFF000000) == 0) {n += 8; x <<= 8;}
if ((x & 0xF0000000) == 0) {n += 4; x <<= 4;}
if ((x & 0xC0000000) == 0) {n += 2; x <<= 2;}
if ((x & 0x80000000) == 0) n += 1;
return n;
}
要改成 long long
很簡單,簡單調整一下即可,如下。
int clz(long long x){
if (x == 0) return 64;
int n = 0;
+ if ((x & 0xFFFFFFFF00000000) == 0) {n += 32; x <<= 32;}
if ((x & 0xFFFF000000000000) == 0) {n += 16; x <<= 16;}
if ((x & 0xFF00000000000000) == 0) {n += 8; x <<= 8;}
if ((x & 0xF000000000000000) == 0) {n += 4; x <<= 4;}
if ((x & 0xC000000000000000) == 0) {n += 2; x <<= 2;}
if ((x & 0x8000000000000000) == 0) n += 1;
return n;
}
計算時間一樣是確定的 (deterministic),不過實在有失優雅。主要有三點:
== 0
邏輯可以簡化第一點即採用 !(x & __)
。第二點的做法可能根據考量點而有所不同,比如此文提供一個很有趣的想法,考量最後一步的效能,他利用 bitwise 操作將分支與加法合併:
- if ((x & 0x8000000000000000) == 0) n += 1;
+ n += (x & 1) ^ 1;
若考慮消去 bad smell 的話,採用迴圈可能是更好的做法。雖然效能上可能會稍微降低,不過程式碼會變的簡潔易讀,且容易擴展。搭配前置處理器技巧,想擴展至更多型別便不是問題。
首先將函式調整成迴圈版。
int clz(long long x){
if (x == 0)
return 64;
int n = 0;
long long mask = 0xFFFFFFFF00000000;
int shift = 32;
while (shift) {
if (!(x & mask)) {
n += shift;
x <<= shift;
}
shift >>= 1;
mask <<= shift;
}
return n;
}
可以發現,透過改成迴圈版,我們將與 long long
型別有關的變因獨立在迴圈之前,那是不是調整這些變因就可以適用於其他型別呢?沒錯,並且搭配前置處理器,程式碼的前段變成:
#define __CZ_MAX_ONES 0xFFFFFFFFFFFFFFFF
#define __clz(type, fn_postfix, bit_len, ones) \
static int clz_##fn_postfix(type x) \
{ \
if (!x) \
return bit_len; \
int bits = 0; \
int shift = bit_len >> 1; \
type mask = (type) (ones << shift); \
/* ... */
/**
* Macro to generate "count leading zeros"
* for type which length is less than int64_t
* Named: clz_"fn_postfix"
*/
#define clz(type, fn_postfix) \
__clz(type, fn_postfix, (sizeof(type) * 8), __CZ_MAX_ONES)
如此一來,我們可以輕鬆產生長度在 64 bits 內的所有 clz
,產生方法如下:
#include <stdint.h>
clz(long long, long_long); // clz_long_long
ctz(long long, long_long); // clz_long_long
clz(int, int); // clz_int
clz(short, short); // clz_short
clz(char, char); // clz_char
clz(int8_t, int8); // clz_int8
clz(int16_t, int16); // clz_int16
clz(int32_t, int32); // clz_int32
clz(int64_t, int64); // clz_int64
clz(uint8_t, uint8); // clz_uint8
clz(uint16_t, uint16); // clz_uint16
clz(uint32_t, uint32); // clz_uint32
clz(uint64_t, uint64); // clz_uint64
上百行的函式就這麼樸實無華的產生了。ctz
同理,只要進行如下的調整…
注意到第一個差異使用 ~
而非反方向位移是因為右移與型別 (type
) 的符號有關,對於有號數,全一再怎麼 (算數) 右移都是全一。
#define __ctz(type, fn_postfix, bit_len, ones) \
static int ctz_##fn_postfix(type x) \
{ \
if (!x) \
return bit_len; \
int bits = 0; \
int shift = bit_len >> 1; \
- type mask = (type) ones << shift; \
+ type mask = (type) ~(ones << shift); \
while (shift) { \
if (!(x & mask)) { \
bits += shift; \
- x <<= shift; \
+ x >>= shift; \
} \
shift >>= 1; \
- mask <<= shift; \
+ mask >>= shift; \
} \
return bits; \
}
如此,我們又可以使出四兩撥千斤之術生成 ctz_XXX
函式們,一勞永逸 :)
測試的部分,可利用與 __builtin_c{l,t}z
來比較輸出,在此省略細節。不過我也是因此注意到 __builtin_c{l,t}z
的 undefined behavior。
恩,我不知道自己多久以後才會完善測試與筆記,還是直接放連結吧。
至於 Repo 中怎麼沒什麼東西…不好意思等我一…下 QwQ
Shiritai
Repo 的部分…
循著作業說明操作的紀錄。
ls -l /dev/fibonacci
crw------- 1 root root 510, 0 三 8 23:51 /dev/fibonacci
可以看見 510, 0
兩數,分別為 device major/minor number。
看完作業說明影片後發現自己讀+實作/實驗太少了,先回去讀…