contributed by < AndybnACT
>
linux2020
第一次執行的時候發現費氏數列的計算結果不如預期。
推測是因為整數相加溢位的緣故,打開輸出檔案 out
可以確認這個現象:
理論上 ;然而這邊輸出竟然和前一步一樣。打開 fibdrv.c
之後可以發現程式碼為了避免溢位,透過 #define MAX_LENGTH 92
將最大的可輸出範圍限制在 。為了輸出後面的費氏數列,我們必須實作大數運算。這次作業的實作先嘗試將兩個無號長整數接在一起來表示一個 128 位元的無號大整數。
gcc 有 128 位元整數的擴充特徵,可見 128-bit Integers:
__int128
unsigned __int128
另外可參照 recipes/basic/int128.h 得知 x86_64 對於 128-bit 整數運算相關的 inline assembly。
:notes: jserv
因為 read()
和 write()
回傳值型態 ssize_t
在 LP64 之下只有 64 位元,如果 client.c
還是用一樣的介面讀取資料的話,將無法透過一次 read()
拿到費氏數列的值。所以,我們改變 fibdrv
的介面,用 read()
的第二個參數、在 kernel 中透過 copy_to_user
將計算結果回傳給 user space 程式。此外,為了之後測量時間方便,read()
的回傳值為費氏數列函數的計算時間。
變更程式介面之後,client.c
也要做對應的修改,這邊為了省去進制轉換的麻煩,直接以兩個 %llx
做輸出。然後,感恩讚嘆 python
,讓 verify.py
讀懂這樣的輸出只需要做很簡單的修改。
client.c
:verify.py
:兩個大整數的相加可以先拆成高位與低位,然後再處理低位相加的 overflow 問題。其中,判斷式 ~op2.lower
代表 op2.lower
與該無號整數上界的距離。所以若 op1.lower
大於該距離,則代表相加會導致溢位。
完成之後,就可以來驗證費氏數列的輸出。經過測試,兩個無號長整數組合而成 128 位元的無號整數可以正確地輸出費氏數列直到 。
同理,兩個大整數的相減可以先拆成高位與低位,然後再處理低位相減的 underflow 問題。當被減數小於減數時,處理 underflow 的方法是將輸出的高位再減去 1。
左移的方法是將高位低位都執行邏輯左移(logical shift),然後把低位應該移入高位的位元補上。需要注意的是 C 語言在處理 shift operand 大於等於變數的寬度時,行為未定義。
詳請見 C99 spec
If the value of the right operand is negative or is greater than or equal to the width of the promoted left operand, the behavior is undefined.
考慮以下變更:
要點:
uint8_t
表達 128 位元的範圍;unlikely
巨集,提示編譯器產生對 branch predictor 友善的程式碼;else
,而且程式縮排更精簡:notes: jserv
已經改善,見 commit 0ade862
同理,右移之前,先判斷移動量。若移動量為零則直接將輸入複製到輸出;右移移動量大於變數寬度時,直接將輸出的高位清空,然後低位即等於原本高位右移 shift - 64
的值。
TODO: 解釋 popcountll 的作用和考量點
:notes: jserv
這邊實作乘法的演算法是最簡單的二進位直式乘法;把所有的乘法換成加法和左移的組合。因為已經實作 128 位元的加法和左移,所以直接呼叫就可以。函式前半段用 popcountll
來抓出 1
位元比較少的數當做乘法直式中的乘數,以減少加法運算的次數。但是 gcc(9.2)
在編譯 kernel module、最後 linking 階段找不到這個函式。所以先暫時把 op1
和 op2
直接指定為 lo
和 hi
。
b101011
,則 popcount 就會回傳 4。而其考量點可以先考慮下列兩個二進位直式乘法因乘法交換率的緣故,兩邊的結果是相同的。然而右邊的被乘數含有 1
的位元明顯大於乘數,使得要重複執行許多次加法才能夠獲得和左邊只使用一次加法就能得到的結果。
實作完以上函式之後,就可以用來測試費氏數列的 Fast doubling () 的演算法。
看完〈How to implement bignum arithmetic〉之後,得知透過使用更大的資料型態來儲存較小數值的相加結果,可以用來保存 carry。為了正確的實作這種加法概念,我們需要先將 struct uint128
的順序改成低位先、高位在後。
加法的實作如下:
但是相較於原本的做法,這樣做只省去了一個 branching、卻換來額外兩個 addition。除非 64 位元的機器處理 32 位元的加法有比較快、而且這些運算可以被放在 32 位元的加法當中,可能才會有效能改進。
初步查了一下 x86 各指令的 cycle count,發現 64 位元加法和 32 位元加法的花費是一樣的。
Comba multiplication
在該篇演講中,作者也有提到他實作大數乘法的方式。一樣也是直式乘法,然而他並沒有將乘法拆成加法,因為硬體本身還是有乘法器可以用,只是寬度不夠而已。解決寬度不夠的方法一樣也是用 64 位元的資料型態儲存兩個 32 位元的乘法結果。
相較於加法的情境,這麼運用顯得比較有效率,因為兩個 32 位元的整數相加最多只會變成 33 位元、儲存在 64 位元的空間中浪費了 31 位元;然而相乘的結果有機會變成 64 位元,減少浪費的情形。
src
Comba 演算法其實也是直式乘法,相較於一般的直式乘法是 row-based,一行一行將結果相加到結果中;Comba 採用 column-based,一列一列寫到結果裡面。這麼做在寫入時對 wb-cache 的壓力相對降低許多,因此可能有比較好的效能。
為了比較 Comba 和一般的直式乘法,先實作(運用乘法單元的)直式乘法
mul128_school
不是很好理解的名稱,可改為 mul128_naive
(不是 native
,而是 naive
,後者的意思是「天真幼稚」)。
:notes: jserv
已經改善,見 commit 0ade862
而採用 Comba 演算法的實作如下:
實作完成之後,可以在 user-space 透過 ADDRESS_SANITIZER
檢查有沒有越界存取,編譯時選項需要加上 -fsanitize=address -fno-omit-frame-pointer -fno-common
。
gcc
也有提供 128 位元整數的操作,[unsigned] __int128
,透過實驗可以觀察其底層實作。實驗的方式為寫一段用到 __int128
運算的 C 程式碼,編譯成執行檔之後觀察執行結果的正確性,然後用反組譯工具 objdump -d -M intel
分析實作方式:
__int128
加法看來 gcc
用符合預期的寬度來儲存 __int128
型態的變數、而且加法結果正確。
反組譯的結果可以看到一道特別的指令 adc
,根據 Intel® 64 and IA-32 Architectures
Software Developer’s Manual 顯示
Adds the destination operand (first operand), the source operand (second operand), and the carry (CF) flag and stores the result in the destination operand.
…
The ADC instruction does not distinguish between signed or unsigned operands. Instead, the processor evaluates the result for both data types and sets the OF and CF flags to indicate a carry in the signed or unsigned result, respectively. The SF flag indicates the sign of the signed result.
比較讓人注意到的是,adc
這個指令不用管輸入是無號還是有號整數。CPU 兩個都會做,然後再用不同的 FLAG
來存放有號、無號整數相加後的 carry-bit。回頭看看一般的 add
指令,發現其實也有這樣的敘述。多虧了 2 補數的表示法,讓有號、無號數運算方式是一致的,所以硬體不需要另外再做一次加法。
由此也可以看出 gcc
對 128 位元變數的安排一樣也是先放低位再放高位($rax/$rcx
為低位,低位做完之後 adc
將低位的 carry 加進高位)。
__int128
位元左右移為了測試左右移,我們嘗試對 x
左移 8 bits。實際測試和反組譯的結果如下:
又發現一道有趣的指令 shld
,而 Intel® 64 and IA-32 Architectures
Software Developer’s Manual 描述:
The SHLD instruction is used for multi-precision shifts of 64 bits or more.
The instruction shifts the first operand (destination operand) to the left the number of bits specified by the third operand (count operand). The second operand (source operand) provides bits to shift in from the right (starting with bit 0 of the destination operand).
這告訴我們,shld
會把目標暫存器左移指定的位元數、而且還會幫我們把來源暫存器的高位也一併移入目標暫存器中。需要注意移動位元數不得超過暫存器寬度。
此外,實驗發現在左移超過 64 位元時,gcc
產生的程式碼與我們的實作類似,將低位寫入高位之後,左移 shift - 64
位元,然後再將低位清空。
__int128
減法在查詢指令的時候有看到 sbb
,似乎也是與 carry 有關的減法指令。實際執行測試時,令 x = 0
, y = 1
,結果如下:
而開發手冊是這麼寫的:
SBB—Integer Subtraction with Borrow:
Operation:
DEST ← (DEST – (SRC + CF)
因此,在執行減法的時候,低位先相減,然後用 sbb
減去高位。低位相減產生的 carry 會拿去高位減。
__int128
乘法乘法的測試中,我們令 x
為無號 64 位元整數的最大值,y
則是 0xDEADBEEF
,兩者相乘。結果如下:
經查發現其實 mul
指令在 operand 大小為 64 位元的資料時,會用 rdx:rax
的組合來儲存 128 位元的乘法結果。
這邊的乘法實作方式是將 128 位元的整數切成兩個 64 位元的整數,然後用一般的直式乘法將兩數相乘。理論上兩個 128 位元的整數相乘應該會需要 256 位元的整數來儲存其結果。但是因為輸出(z
)也只有 128 位元,所以這裡沒有做 x.hi * y.hi
、也把x.lo * y.hi
以及 y.lo * x.hi
高於 128 位元的部分捨去。
吸收完 gcc
對 __int128
操作的美妙之後,決定用類似的方式實作真正任意大小的大數運算,姑且依標頭檔 abn.h
命名為 ABN。為了讓測試更方便,我們先在 user space 開發測試。
一樣也是用陣列的方式將多個無號長整數組合起來,表達一個大的無號長整數「大數」。
其中,cnt
代表目前大數佔用了幾個無號長整數的空間;cap
代表目前的變數儲存空間允許置放多少個長整數;而 heap
及 stack
則代表儲存空間的實際位置。存放方式一樣是低位在前高位在後。
顧名思義,為了讓數字在相對小的時候減少記憶體配置器(allocator)的開銷,我們預先在 struct
裡面分配 uint64_t * STACKALLOC
的大小。待該空間已經不夠存放的時候,重新用記憶體配置器於 heap 中分配可儲存的空間。
如果在使用 stack 的時候也把儲存 heap
指標的空間拿來存放大數的話,空間運用會更有效率。但為了初期開發方便,這邊分開用兩個變數儲存。
之後,在實際執行運算時,我們需要一個能判斷大數存放在 heap 還是 stack 的方法。我們的做法是:先將 heap
初始化為 NULL
,重新分配空間之後 heap
都是有值的變數。因此,取數只需要這麼做即可:
而重新分配空間的方式和 lib C realloc
的方式差不多,先配置指定大小的空間,將該空間初始化,然後帶入舊空間的值之後,將舊的空間釋放:
基本的資料配置、存取完成之後,就可以實作大數運算。所有大數運算的邏輯都是這樣:先判斷目標運算元的空間是否足夠,不足則配置空間;然後按照順序執行運算的同時,考慮長整數邊界的進位問題;最後目標寫入完成之後,更新目標運算元大數結構內的 cnt
。
開始進行加法之前,我們先透過 cnt
找出位元數比較少的數當作被加數。這樣一來,當低位都相加完之後,我們就只剩下低位加法產生的 carry,和高位要處理。而加法運算本身,為了直接拿到計算而產生的 carry,我們透過 inline assembly 來完成。
不像 __int128
可以直接在緊鄰的高位運算用 adc
指令將 carry 帶入(除非展開迴圈),我們的 inline assembly 透過 setc
這個指令(set byte if CF == 1
)把每次計算的 carry 帶回 C 程式碼,然後迴圈下一步再帶入上次得到的 carry。因為 carry 最多就是一,整個加法最多只會進位一次,所以可以一次把他們加起來。
減法的邏輯與加法一致,不過在減法中我們不需要找出比較小的數,而是直接將兩個運算原由低位置高位逐一相減。
在減法的過程中,因為不確定減法結果的有效位元為何,所以只要減出來的結果不為零,就更新有效位元 cnt
。
在不失一般性的前提下,這邊以左移為例。我們可以善用 shld
這個有趣的指令,在左移高位的同時,幫忙把低位補進來。但要記得的是左移之前,必須要先把高位暫存出來。不然在高位左移後的下一步,我們會找不到要補進更高位的那些位元。
這段程式碼被以 -O2
以上編譯的時候,會有預期之外的行為,主要原因是:於除非在 :
區段標注清楚, gcc
實際上無從直接得知 inline assembly 裡面對暫存器、記憶體位置等的使用情形。因此,當 inline assembly 裡面用到輸入、輸出以外的儲存空間時(暫存器、記憶體),需要在 clobber list 中通知 gcc
「不能假設」在 inline assembly 前後,這些「空間」的值是一致的。因此,暫時能解決問題,但不是很完善的方式為:在 clobber list 裡面加上 "memory"
。「不是很完善」的原因是,告訴 gcc
整個記憶體在 inline assembly 操作之後都不能假設其一致性的成本似乎過大。
-O2
)乘法沿用 Comba multiplication,不過在發現 mul
指令做 64 位元乘法時可以保留 128 位元的運算結果後,我們用 inline assembly 把這樣的結果拉回 C 程式語言使用。
乘法和加法類似,先找出位元數( cnt
)比較少的數值當作被乘數,寫在直式乘法下半,然後才由結果的低位開始計算,依序往高位移動、寫入。因為每次加回結果的數都是 64 位元,我們可以直接沿用 __add_ll()
來幫忙這項操作、並且兼顧 carry 的問題;同時,記得在乘法開始之前先將結果初始化為零。
在計算費氏數列時,若採用 Fast Doubling 的演算法,可以較快速地()迭代至要求目標 。原因在於這種演算法將等式左右的關係矩陣化,然後找出一個 對應到 的關係。而二進位的數字又很容易用上這種關係,舉例來說, 的二進位表達式為 b1101
,也就是
這樣一來,我們只需要從 的 MSB 開始往右迭代,每次都計算 和 ,若遇到位元是一的話只要再算 ,然後繼續迭代。而 clz
是硬體提供的指令,count leading zeros,可以讓我們快速地找到 裡面,MSB 所在的位置。
fib_sequence_ll
:若計算結果可以保持在 64 位元的無號長整數內的話,一般連加演算法可以直接從定義寫出來:fib_sequence
:若考慮大數,以我們目前的實作可以延長到 128 位元。因為是自定義的型態,呼叫函式做加法運算:fibseq_doubling_ll
:Fast Doubling 快速演算法,若使用內建型態只需要直接使用內建運算子。需要注意的是 __builtin_clzll()
在輸入為零的時候,行為未定義,所以需要預先處理:fibseq_doubling
: Fast Doubling 考慮大數的實作如下:根據作業描述提到的方式設定實驗環境。
isolcpus=0
為了方便切換實驗環境,將上述設定寫在一個 script 裡面:
因為執行時期只有 client
一個 process,所以只有將一個 CPU 的組態設定為 performance。
在 fibdrv.c
裡面,我們透過 write()
介面變更計算費波那契數的函式。目前實作的四種函式一使用的資料型態和演算法分類如下,其中,目前自訂 uint128_t
的算數用前述(add128
, mul128
)的方式實作。
自訂 uint128_t | 內建 u64 | |
---|---|---|
一般演算法(Regular) | fib_sequence() |
fib_sequence_ll() |
快速演算法(Fast Doubling) | fibseq_doubling() |
fibseq_doubling_ll() |
uint128_t
加法中,會用到兩次加法和一次 branching;所以也代表 branch overhead 其實不高。聚焦在下圖,可以看到 Fast doubling 的計算效能確實較一般演算法來得好、也大致呈現 的成長。然而,在上圖中,即便迭代次數明顯較少,使用 mul128
大數乘法實作的 fibseq_doubling()
竟也大致呈 成長、且效能明顯不佳。
add128_auto_carry()
和原先實作的加法 add128()
效能 :add128()
已經使用 loop-unrolling 的最佳化手法加深 pipeline ,還是顯著地慢上許多。透過 objdump -M intel
可以看到在加法的部分其實也是使用 64 位元的加法實作。fibseq_doubling()
,置換裡面的乘法函式來做測試。乘法迴圈一樣也採用 loop unrolling 的最佳化手法:TODO: 評估 sysprog21/bignum 的表現並探討實作手法
:notes: jserv