HackMD
2022q1 Homework2 (quiz2)
contributed by < KHLee529 >
測試環境
$ gcc --version
gcc (Ubuntu 9.3.0-17ubuntu1~20.04) 9.3.0
$ lscpu
Architecture: x86_64
CPU op-mode(s): 32-bit, 64-bit
Byte Order: Little Endian
Address sizes: 39 bits physical, 48 bits virtual
CPU(s): 8
On-line CPU(s) list: 0-7
Thread(s) per core: 2
Core(s) per socket: 4
Socket(s): 1
NUMA node(s): 1
Vendor ID: GenuineIntel
CPU family: 6
Model: 94
Model name: Intel(R) Core(TM) i7-6770HQ CPU @ 2.60GHz
Stepping: 3
CPU MHz: 2600.000
CPU max MHz: 3500.0000
CPU min MHz: 800.0000
BogoMIPS: 5199.98
Virtualization: VT-x
L1d cache: 128 KiB
L1i cache: 128 KiB
L2 cache: 1 MiB
L3 cache: 6 MiB
NUMA node0 CPU(s): 0-7
測驗 1
- 解釋上述程式碼運作的原理
- 比較上述實作在編譯器最佳化開啟的狀況，對應的組合語言輸出，並嘗試解讀 (可研讀 CS:APP 第 3 章)
- 研讀 Linux 核心原始程式碼 include/linux/average.h，探討其 Exponentially weighted moving average (EWMA) 實作，以及在 Linux 核心的應用
移動平均（Moving average)，又稱滾動平均值、滑動平均，在統計學中是種藉由建立整個資料集合中不同子集的一系列平均數，來分析資料點的計算方法。
移動平均通常與時間序列資料一起使用，以消除短期波動，突出長期趨勢或周期。短期和長期之間的閾值取決於應用，移動平均的參數將相應地設置。例如，它通常用於對財務數據進行技術分析，如股票價格、收益率或交易量。它也用於經濟學中研究國內生產總值、就業或其他宏觀經濟時間序列。
運作原理
這個直覺的解法會有 overflow 的問題，若我們已知 a, b 數值的大小，可用下方程式避免 overflow:
#include <stdint.h>
uint32_t average(uint32_t low, uint32_t high)
{
return low + (high - low) / 2;
}
接著我們可改寫為以下等價的實作:
#include <stdint.h>
uint32_t average(uint32_t a, uint32_t b)
{
return (a >> 1) + (b >> 1) + (a & b & 1);
}
其中，這個實作方式是透過將 a, b 先除二後再相加，並考慮若 a, b 同時為奇數時除二的餘數相加會在平均中再加 1。
因此，透過 a & b & 1 得到是否 a, b 兩數最末位同時為 1。若是，則結果為 1 ，若不是則結果為 0。
我們再次改寫為以下等價的實作:
uint32_t average(uint32_t a, uint32_t b)
{
return (a & b) + ((a ^ b) >> 1);
}
其中，這個實作方式則是需要考慮到二進位加法與位元運算的相關聯。
在二進位的加法當中，對於其中一個位元的運算方式為
由此規則可以發現，在兩數 a、b 相加時，對應任一位，相加後的結果與 a ^ b 相同，而進位的部分與 a & b 相同。
故 a + b 可以表達成 (a ^ b) + (a & b) << 1。
而在取平均時，即為將此一結果再右移一個位元，即為 (a & b) + (a ^ b) >> 1
編譯優化
原始程式碼
// q1.c
#include <stdint.h>
uint32_t average_v1(uint32_t a, uint32_t b) {
return (a + b) / 2;
}
uint32_t average_v2(uint32_t low, uint32_t high) {
return low + (high - low) / 2;
}
uint32_t average_v3(uint32_t a, uint32_t b) {
return (a >> 1) + (b >> 1) + (a & b & 1);
}
uint32_t average_v4(uint32_t a, uint32_t b) {
return (a & b) + ((a ^ b) >> 1);
}
開啟編譯優化前 (-O0)
$ gcc -c q1.c
$ objdump -d q1.o
Disassembly of section .text:
0000000000000000 <average_v1>:
0: f3 0f 1e fa endbr64
4: 55 push %rbp
5: 48 89 e5 mov %rsp,%rbp
8: 89 7d fc mov %edi,-0x4(%rbp)
b: 89 75 f8 mov %esi,-0x8(%rbp)
e: 8b 55 fc mov -0x4(%rbp),%edx
11: 8b 45 f8 mov -0x8(%rbp),%eax
14: 01 d0 add %edx,%eax
16: d1 e8 shr %eax
18: 5d pop %rbp
19: c3 retq
000000000000001a <average_v2>:
1a: f3 0f 1e fa endbr64
1e: 55 push %rbp
1f: 48 89 e5 mov %rsp,%rbp
22: 89 7d fc mov %edi,-0x4(%rbp)
25: 89 75 f8 mov %esi,-0x8(%rbp)
28: 8b 45 f8 mov -0x8(%rbp),%eax
2b: 2b 45 fc sub -0x4(%rbp),%eax
2e: d1 e8 shr %eax
30: 89 c2 mov %eax,%edx
32: 8b 45 fc mov -0x4(%rbp),%eax
35: 01 d0 add %edx,%eax
37: 5d pop %rbp
38: c3 retq
0000000000000039 <average_v3>:
39: f3 0f 1e fa endbr64
3d: 55 push %rbp
3e: 48 89 e5 mov %rsp,%rbp
41: 89 7d fc mov %edi,-0x4(%rbp)
44: 89 75 f8 mov %esi,-0x8(%rbp)
47: 8b 45 fc mov -0x4(%rbp),%eax
4a: d1 e8 shr %eax
4c: 89 c2 mov %eax,%edx
4e: 8b 45 f8 mov -0x8(%rbp),%eax
51: d1 e8 shr %eax
53: 01 c2 add %eax,%edx
55: 8b 45 fc mov -0x4(%rbp),%eax
58: 23 45 f8 and -0x8(%rbp),%eax
5b: 83 e0 01 and $0x1,%eax
5e: 01 d0 add %edx,%eax
60: 5d pop %rbp
61: c3 retq
0000000000000062 <average_v4>:
62: f3 0f 1e fa endbr64
66: 55 push %rbp
67: 48 89 e5 mov %rsp,%rbp
6a: 89 7d fc mov %edi,-0x4(%rbp)
6d: 89 75 f8 mov %esi,-0x8(%rbp)
70: 8b 45 fc mov -0x4(%rbp),%eax
73: 23 45 f8 and -0x8(%rbp),%eax
76: 89 c2 mov %eax,%edx
78: 8b 45 fc mov -0x4(%rbp),%eax
7b: 33 45 f8 xor -0x8(%rbp),%eax
7e: d1 e8 shr %eax
80: 01 d0 add %edx,%eax
82: 5d pop %rbp
83: c3 retq
開啟編譯優化後 (-O1)
$ gcc -O1 -c q1.c
$ objdump -d q1.o
註：編譯優化等級
-O1到-O3皆呈現相同成果，推論在-O2與-O3開啟的優化選項中並無影響這段程式編譯結果的部分。
q1_1.o: file format elf64-x86-64
Disassembly of section .text:
0000000000000000 <average_v1>:
0: f3 0f 1e fa endbr64
4: 8d 04 37 lea (%rdi,%rsi,1),%eax
7: d1 e8 shr %eax
9: c3 retq
000000000000000a <average_v2>:
a: f3 0f 1e fa endbr64
e: 89 f0 mov %esi,%eax
10: 29 f8 sub %edi,%eax
12: d1 e8 shr %eax
14: 01 f8 add %edi,%eax
16: c3 retq
0000000000000017 <average_v3>:
17: f3 0f 1e fa endbr64
1b: 89 f8 mov %edi,%eax
1d: d1 e8 shr %eax
1f: 89 f2 mov %esi,%edx
21: d1 ea shr %edx
23: 01 d0 add %edx,%eax
25: 21 f7 and %esi,%edi
27: 83 e7 01 and $0x1,%edi
2a: 01 f8 add %edi,%eax
2c: c3 retq
000000000000002d <average_v4>:
2d: f3 0f 1e fa endbr64
31: 89 f8 mov %edi,%eax
33: 31 f0 xor %esi,%eax
35: d1 e8 shr %eax
37: 21 f7 and %esi,%edi
39: 01 f8 add %edi,%eax
3b: c3 retq
由反組譯結果可以看出，在未開啟編譯優化時的編譯結果，基本上都是按照原始程式碼的邏輯進行編譯。
意即對於一個 + 運算便會有一個 add 指令，有一個 - 運算便有一個 sub 指另與其對應⋯⋯等等。
而在編譯優化後，可以由 average_v1 的反組譯結果看出，編譯器對於兩個數值相加的運算簡單的以取址的 lea 指令進行，其餘在運算部分並無特別大的差異。
但在另一方面，由反組譯結果可以看到，在編譯優化開啟與否影響到實際指令長度的關鍵差異並非在運算部分的指令，而是未開啟編譯優化時，編譯器會先將所有的引數(arguments) 放入 stack 當中，
再對其進行相關運算。這個變化待閱讀 gcc manual 中有關編譯優化選項的部分後補上。而為何 -O0 的編譯結果會將引數放入 stack 中再進行運算目前沒有頭緒待研究。
EWMA in linux/average.h
在 linux/average.h 當中 EWMA 的實作以巨集宣告函數的方式實作，其中包含三個函數 ewma_name_init、ewma_name_read 與 ewma_name_add。
而這三個函數的意義則從 EWMA 的計算方式出發
由此定義可以理解為何會有 add 相關函數的存在，對於每一個新的值，皆會以 ewma_name_add 函數即可得到加入新值的結果。
而在 Linux 的實作當中可以將各個函數主要分成兩個部分，參數錯誤預防與運算部分。
-
錯誤預防部分
針對本巨集的使用，除了精度與權重的數值必須以編譯時期即可得到數值的常數方式宣告以外，亦訂定了對於精確度的規定以及權重需為
其中的BUILD_BUG_ON(condition)巨集得以在編譯時期針對condition為true的情況進行編譯錯誤的警告。將其展開後可以得到BUILD_BUG_ON_MSG(condition, "BUILD_BUG_ON failed : " #condition)在展開為
compiletime_assert(!(condition), "BUILD_BUG_ON failed : " #condition)最後展開為
do { __attribute__((__noreturn__)) extern void __compiletime_assert_COUNT(void) __attribute__((__error__("BUILD_BUG_ON failed : " #condition))); if (!(!(condition))) __compiletime_assert_COUNT(void); } while (0)其中
COUNT的部分會由編譯器提供的__COUNT__巨集得到唯一數值。
透過此一實作，對於編譯時期即知道真偽結果的condition，透過編譯器優化決定是否跳過第 5 行的函數呼叫。
當condition為真時，會保留第 5 行的函數呼叫，而由於該函數並沒有實際實作，便會引起連結器錯誤，亦或是透過__attribute__((__error__(msg)))的功能，於函數呼叫時即會引發編譯錯誤，使得該檢查於編譯時期即有效。
而在這部分當中，由於針對此一巨集當中，只要其中一個編譯錯誤發生即會停止編譯，而為何需要於三個函數皆導入相同的錯誤預防機制則需往後繼續研究。 -
運算部分
若將巨集的實作當中錯誤預防的部分去除後剩下以下部分。#define DECLARE_EWMA(name, _precision, _weight_rcp) \ struct ewma_##name { \ unsigned long internal; \ }; \ static inline void ewma_##name##_init(struct ewma_##name *e) \ { \ e->internal = 0; \ } \ static inline unsigned long \ ewma_##name##_read(struct ewma_##name *e) \ { \ return e->internal >> (_precision); \ } \ static inline void ewma_##name##_add(struct ewma_##name *e, \ unsigned long val) \ { \ unsigned long internal = READ_ONCE(e->internal); \ unsigned long weight_rcp = ilog2(_weight_rcp); \ unsigned long precision = _precision; \ \ WRITE_ONCE(e->internal, internal ? \ (((internal << weight_rcp) - internal) + \ (val << precision)) >> weight_rcp : \ (val << precision)); \ }透過此一簡化過後的結果可以看出，當需要使用 EWMA 的功能時，需要先呼叫
init函數對資料初始化後透過add加入各資料點後以read讀取其中資料。
而在主要進行 EWMA 功能運算的add當中可以看到其透過大量的運算縮短運算時間。以下為運算主體程式碼。WRITE_ONCE(e->internal, internal ? (((internal << weight_rcp) - internal) + (val << precision)) >> weight_rcp : (val << precision));其中
WRITE_ONCE巨集為確保資料寫入記憶體的並行相關處理。在此先將其簡化為以下賦值表示式以關注其運算方式。e->internal = internal ? (((internal << weight_rcp) - internal) + (val << precision)) >> weight_rcp : (val << precision);由以上程式碼可看出，若當
internal數值為 0 時，即甫初始化結束第一次加入數值時，會將結果設定為(val << precision)
由此可以看出其e->internal的內部儲存結構為由此簡易的浮點數表示方式，使得加權平均的分數運算得以快速的進行。接著僅需要在加入每一個數值時搭配
在此部分由於巨集以將 \alpha 限制為val << precision)，再將其右移(val << precision) >> weight_rcp)。 而這部分對於巨集作者在
internal左移再右移回來而非直接減去右移的結果有些疑惑仍待研究。
由於此部分若internal的整數部份儲存實際使用到的位元數較大，leading zeros 的數量小於weight_rcp時，即會在左移時產生 overflow 產生錯誤。而初步推測原因為，若weight_rcp數值很大時，每一次 EWMA 幾乎即等於其新加入的值，使得移動平均的意義與效果大幅削弱，因此weight_rcp應不會太大，且若每次加入的precision以避免 overflow。此部分待往後研究其使用案例時進行推測。
而在 Linux 核心當中主要使用到 EWMA 為無線網路相關程式碼。
$ grep -r DECLARE_EWMA
net/mac80211/ieee80211_i.h:DECLARE_EWMA(beacon_signal, 4, 4)
net/mac80211/sta_info.h:DECLARE_EWMA(avg_signal, 10, 8)
net/mac80211/sta_info.h:DECLARE_EWMA(mesh_fail_avg, 20, 8)
net/mac80211/sta_info.h:DECLARE_EWMA(mesh_tx_rate_avg, 8, 16)
net/mac80211/sta_info.h:DECLARE_EWMA(signal, 10, 8)
net/batman-adv/types.h:DECLARE_EWMA(throughput, 10, 8)
drivers/net/wireless/ath/ath5k/ath5k.h:DECLARE_EWMA(beacon_rssi, 10, 8)
drivers/net/wireless/realtek/rtw89/core.h:DECLARE_EWMA(tp, 10, 2);
drivers/net/wireless/realtek/rtw89/core.h:DECLARE_EWMA(rssi, 10, 16);
drivers/net/wireless/realtek/rtw89/core.h:DECLARE_EWMA(thermal, 4, 4);
drivers/net/wireless/realtek/rtw88/main.h:DECLARE_EWMA(tp, 10, 2);
drivers/net/wireless/realtek/rtw88/main.h:DECLARE_EWMA(rssi, 10, 16);
drivers/net/wireless/realtek/rtw88/main.h:DECLARE_EWMA(thermal, 10, 4);
drivers/net/wireless/realtek/rtw88/main.h:DECLARE_EWMA(evm, 10, 4);
drivers/net/wireless/realtek/rtw88/main.h:DECLARE_EWMA(snr, 10, 4);
drivers/net/wireless/ralink/rt2x00/rt2x00.h:DECLARE_EWMA(rssi, 10, 8)
drivers/net/wireless/mediatek/mt7601u/mt7601u.h:DECLARE_EWMA(rssi, 10, 4);
drivers/net/wireless/mediatek/mt76/mt76.h:DECLARE_EWMA(signal, 10, 8);
drivers/net/wireless/mediatek/mt76/mt7921/mt7921.h:DECLARE_EWMA(rssi, 10, 8);
drivers/net/wireless/mediatek/mt76/mt76x02_mac.h:DECLARE_EWMA(pktlen, 8, 8);
drivers/net/wireless/intel/iwlwifi/mvm/mvm.h:DECLARE_EWMA(rate, 16, 16)
drivers/net/virtio_net.c:DECLARE_EWMA(pkt_len, 0, 64)
drivers/gpu/drm/drm_self_refresh_helper.c:DECLARE_EWMA(psr_time, 4, 4)
drivers/gpu/drm/i915/gt/intel_context_types.h:DECLARE_EWMA(runtime, 3, 8);
drivers/gpu/drm/i915/gt/intel_engine_types.h:DECLARE_EWMA(_engine_latency, 6, 4)
include/linux/average.h:#define DECLARE_EWMA(name, _precision, _weight_rcp)
測驗 2
- 解釋上述程式碼運作的原理
- 針對 32 位元無號/有號整數，撰寫同樣 branchless 的實作
- Linux 核心也有若干 branchless / branch-free 的實作，例如 lib/sort.c:
請在 Linux 核心原始程式碼中，找出更多類似的實作手法。請善用 git log 檢索。
/* * Logically, we're doing "if (i & lsbit) i -= size;", but since the * branch is unpredictable, it's done with a bit of clever branch-free * code instead. */ __attribute_const__ __always_inline static size_t parent(size_t i, unsigned int lsbit, size_t size) { i -= size; i -= size & -(i & lsbit); return i / 2; }
運作原理
#include <stdint.h>
uint32_t max(uint32_t a, uint32_t b)
{
return a ^ ((a ^ b) & -(a < b));
}
從 Xor 的幾個特性中可以看出端倪，分別為
x ^ 0 = xx ^ x = 0x ^ y ^ y = x
故，若欲回傳 a 時，即((EXP4) & -(EXP5)) 的部分需為 0。
而若要回傳 b 時，即((EXP4) & -(EXP5)) 的部分需為 a ^ b。
透過 EXP5 為一個比較操作且前面加上一個負號的提示可知當 EXP5 為真時，((EXP4) & -(EXP5)) 的結果為 EXP4，而當 EXP5 為偽時，((EXP4) & -(EXP5)) 的結果為 0。整理如下，
EXP5 |
-(EXP5) |
(EXP4) & -(EXP5) |
|---|---|---|
True (1) |
-1 (0xFFFFFFFF) |
EXP4 |
False (0) |
0 (0x00000000) |
0 |
在配合回上方推論，即可得知 EXP4 為 a ^ b，且當 a 較大時 EXP5 需為偽，故推論 EXP5 為 a < b。
有號/無號數對應實作
由於運作原理中使用到的特性對於有號與無號數皆適用，故對於有號數亦可直接套用。得到結果如下
#include <stdint.h>
int32_t max(int32_t a, int32_t b)
{
return a ^ ((a ^ b) & -(a < b));
}
int32_t min(int32_t a, int32_t b)
{
return a ^ ((a ^ b) & -(a > b));
}
Branchless/Branch-free
透過以下指令搜尋在 git log 當中 commit message 包含 "branchless" 或 "branch-free" 相關內容的。
$ git log --oneline --name-only -E --grep "branchless|branch-free"
e7de0023e123 powerpc/ebpf32: Rework 64 bits shifts to avoid tests and branches
arch/powerpc/net/bpf_jit_comp32.c
efdb25efc764 arm64/lib: improve CRC32 performance for deep pipelines
arch/arm64/lib/crc32.S
6c786bcb29dd xfrm: branchless addr4_match() on 64-bit
include/net/xfrm.h
6bb69c9b69c3 KVM: MMU: simplify last_pte_bitmap
arch/x86/include/asm/kvm_host.h
arch/x86/kvm/mmu.c
97ec8c067d32 KVM: Add SMAP support when setting CR4
arch/x86/kvm/cpuid.h
arch/x86/kvm/mmu.c
arch/x86/kvm/mmu.h
arch/x86/kvm/paging_tmpl.h
arch/x86/kvm/x86.c
747f49ba67b8 Replace int2float() with an optimized version.
drivers/gpu/drm/radeon/r600_blit.c
97d64b788114 KVM: MMU: Optimize pte permission checks
arch/x86/include/asm/kvm_host.h
arch/x86/kvm/mmu.c
arch/x86/kvm/mmu.h
arch/x86/kvm/paging_tmpl.h
arch/x86/kvm/x86.c
很不幸的在這當中大部分的程式碼變更相對時間比較久遠導致都已經被修改掉了，因此嘗試以其他關鍵字搜尋相關的 branchless 操作。
而在搜尋關鍵字 "eliminate branch" 時看到如下的結果，幸運的找到了目前還保留的程式碼。
$ git log --oneline --name-only --grep "eliminate branch"
6d874eca6595 lib: bitmap: eliminate branch in __bitmap_shift_left
lib/bitmap.c
9d8a6b2a02c5 lib: bitmap: eliminate branch in __bitmap_shift_right
lib/bitmap.c
diff --git a/lib/bitmap.c b/lib/bitmap.c
index 74bdf3601245..36e380da00c5 100644
--- a/lib/bitmap.c
+++ b/lib/bitmap.c
@@ -169,15 +169,14 @@ void __bitmap_shift_left(unsigned long *dst, const unsigned
long *src,
* word below and make them the bottom rem bits of result.
*/
if (rem && k > 0)
- lower = src[k - 1];
+ lower = src[k - 1] >> (BITS_PER_LONG - rem);
else
lower = 0;
upper = src[k];
if (left && k == lim - 1)
upper &= (1UL << left) - 1;
- dst[k + off] = upper << rem;
- if (rem)
- dst[k + off] |= lower >> (BITS_PER_LONG - rem);
+ upper <<= rem;
+ dst[k + off] = lower | upper;
if (left && k + off == lim - 1)
dst[k + off] &= (1UL << left) - 1;
}
diff --git a/lib/bitmap.c b/lib/bitmap.c
index 45e7d14ebdfd..a7a8bc02892d 100644
--- a/lib/bitmap.c
+++ b/lib/bitmap.c
@@ -129,13 +129,13 @@ void __bitmap_shift_right(unsigned long *dst, const unsigned long *src,
upper = src[off + k + 1];
if (off + k + 1 == lim - 1 && left)
upper &= mask;
+ upper <<= (BITS_PER_LONG - rem);
}
lower = src[off + k];
if (left && off + k == lim - 1)
lower &= mask;
- dst[k] = lower >> rem;
- if (rem)
- dst[k] |= upper << (BITS_PER_LONG - rem);
+ lower >>= rem;
+ dst[k] = lower | upper;
if (left && k == lim - 1)
dst[k] &= mask;
}
透過這個兩個 patch 紀錄可以清楚的看到它是如何的從 if statement 轉換成 branchless 的程式碼。
測驗 3
- 解釋上述程式運作原理;
- 在 x86_64 上透過
__builtin_ctz改寫 GCD，分析對效能的提升; - Linux 核心中也內建 GCD (而且還不只一種實作)，例如 lib/math/gcd.c，請解釋其實作手法和探討在核心內的應用場景。
運作原理
#include <stdint.h>
uint64_t gcd64(uint64_t u, uint64_t v)
{
if (!u || !v) return u | v;
int shift;
for (shift = 0; !((u | v) & 1); shift++) {
u /= 2, v /= 2;
}
while (!(u & 1))
u /= 2;
do {
while (!(v & 1))
v /= 2;
if (u < v) {
v -= u;
} else {
uint64_t t = u - v;
u = v;
v = t;
}
} while (v);
return u << shift;
}
這個實作方式是將公質因數中 2 的部分先行提出，大幅降低數值的大小範圍，使得後面的迭代次數得以降低。
int shift;
for (shift = 0; !((u | v) & 1); shift++) {
u /= 2, v /= 2;
}
首先這個部分即為將兩者公因數中 2 的次方數部分先提出，並以右移的次數(除以 2 的次數)作為記錄方式。
while (!(u & 1))
u /= 2;
再來在這個部分，若 u 的質因數分解中仍有 2 的次方部分，該部分必不為兩者的公因數，故可以將其剔除。
do {
while (!(v & 1))
v /= 2;
if (u < v) {
v -= u;
} else {
uint64_t t = u - v;
u = v;
v = t;
}
} while (v);
最後的部分就如同原本的輾轉相除法，只是由於先前已將質因數中 2 的次方部分提出，所以可以在過程中，再將質因數中 2 的次方部分剔除儘速降低數值大小。
也因此，由於此部分的計算即為輾轉相除法，結束判斷的標準仍為 v 是否等於 0。故 COND 部分為 v。
return u << shift;
如同輾轉相除法的結果，最後 u 即為最大公因數，但由於原先還有 2 的次方部分已經被提出，在回傳時需透過左移位元恢復該部分。故 RET 部分為 u << shift
使用 __builtin_ctz 改寫
將以上實作方式以 __builtin_ctz 相關函數進行改寫後，結果如下
int min(int a, int b) {
return a ^ ((a ^ b) & -(a > b));
}
uint64_t gcd64_use_builtin_ctz(uint64_t u, uint64_t v)
{
if (!u || !v) return u | v;
int shift, u_tz, v_tz;
u_tz = __builtin_ctzll(u);
v_tz = __builtin_ctzll(v);
shift = min(u_tz, v_tz);
u >>= u_tz;
v >>= v_tz;
do {
v >>= __builtin_ctzll(v);
if (u < v) {
v -= u;
} else {
uint64_t t = u - v;
u = v;
v = t;
}
} while (v);
return u << shift;
}
效能測試相關程式碼
#include <stdint.h>
#include "gcd.h"
uint64_t gcd64(uint64_t u, uint64_t v)
{
if (!u || !v) return u | v;
int shift;
for (shift = 0; !((u | v) & 1); shift++) {
u /= 2, v /= 2;
}
while (!(u & 1))
u /= 2;
do {
while (!(v & 1))
v /= 2;
if (u < v) {
v -= u;
} else {
uint64_t t = u - v;
u = v;
v = t;
}
} while (v);
return u << shift;
}
int min(int a, int b) {
return a ^ ((a ^ b) & -(a > b));
}
uint64_t gcd64_use_builtin_ctz(uint64_t u, uint64_t v)
{
if (!u || !v) return u | v;
int shift, u_tz, v_tz;
u_tz = __builtin_ctzll(u);
v_tz = __builtin_ctzll(v);
shift = min(u_tz, v_tz);
u >>= u_tz;
v >>= v_tz;
do {
v >>= __builtin_ctzll(v);
if (u < v) {
v -= u;
} else {
uint64_t t = u - v;
u = v;
v = t;
}
} while (v);
return u << shift;
}
#ifndef _GCD_H_
#define _GCD_H_
typedef uint64_t (gcd_f)(uint64_t, uint64_t);
uint64_t gcd64(uint64_t, uint64_t);
uint64_t gcd64_use_builtin_ctz(uint64_t, uint64_t);
#endif /* _GCD_H_ */
本測試程式碼實作方式取自 Merge Sort 與他的變化
#include <stdio.h>
#include <stdint.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <sys/time.h>
#include "util.h"
#include "gcd.h"
int main() {
srand(time(0));
uint64_t u, v;
const int epochs = 1000;
int interval = epochs / 10;
struct timespec start, end;
gcd_f *gcds[] = {
gcd64,
gcd64_use_builtin_ctz,
};
const int sort_len = sizeof(gcds) / sizeof(gcd_f *);
time_t measures[epochs][sort_len];
for (int i = 0; i < epochs; i++)
{
u = ((uint64_t) rand() << 32) | rand();
v = ((uint64_t) rand() << 32) | rand();
for (size_t j = 0; j < sort_len; j++)
{
clock_gettime(CLOCK_MONOTONIC, &start);
gcds[j](u, v);
clock_gettime(CLOCK_MONOTONIC, &end);
measures[i][j] = diff(start, end);
}
if ((i + 1) % interval == 0)
printf("Benchmark %-3d%%\n", 100 * (i + 1) / epochs);
}
FILE *file = fopen("./benchmark.txt", "w+");
for (size_t i = 0; i < epochs; i++)
{
for (size_t j = 0; j < sort_len; j++)
{
char *token = j == sort_len - 1 ? j == epochs - 1 ? "" : "\n" : "\t";
fprintf(file, "%ld%s", measures[i][j], token);
}
}
return 0;
}
將兩函數進行效能測試結果如下
此結果圖有調整過 y 軸範圍，由於原始數據有數筆花費時間高達數萬 nsec 使得圖形被擠壓至下方，故調整 y 軸範圍以呈現主要現象。
從結果可以看出，使用 __builtin_ctz 系列函數計算結尾的 0 的數量較直接右移至最末位不為零的實作方式快。其原因待研究，目前先以編譯結果來做一個初步的假設。
gcd.o: file format elf64-x86-64
Disassembly of section .text:
0000000000000000 <gcd64>:
0: f3 0f 1e fa endbr64
4: 48 89 f8 mov %rdi,%rax
7: 48 09 f0 or %rsi,%rax
a: 48 85 ff test %rdi,%rdi
d: 74 5a je 69 <gcd64+0x69>
f: 48 85 f6 test %rsi,%rsi
12: 74 55 je 69 <gcd64+0x69>
14: 31 c9 xor %ecx,%ecx
16: a8 01 test $0x1,%al
18: 75 19 jne 33 <gcd64+0x33>
1a: 66 0f 1f 44 00 00 nopw 0x0(%rax,%rax,1)
20: 48 d1 ef shr %rdi
23: 48 d1 ee shr %rsi
26: 83 c1 01 add $0x1,%ecx
29: 48 89 f8 mov %rdi,%rax
2c: 48 09 f0 or %rsi,%rax
2f: a8 01 test $0x1,%al
31: 74 ed je 20 <gcd64+0x20>
33: 40 f6 c7 01 test $0x1,%dil
37: 75 17 jne 50 <gcd64+0x50>
39: 0f 1f 80 00 00 00 00 nopl 0x0(%rax)
40: 48 d1 ef shr %rdi
43: 40 f6 c7 01 test $0x1,%dil
47: 74 f7 je 40 <gcd64+0x40>
49: 0f 1f 80 00 00 00 00 nopl 0x0(%rax)
50: 40 f6 c6 01 test $0x1,%sil
54: 74 1a je 70 <gcd64+0x70>
56: 48 39 f7 cmp %rsi,%rdi
59: 73 1d jae 78 <gcd64+0x78>
5b: 48 29 fe sub %rdi,%rsi
5e: 48 85 f6 test %rsi,%rsi
61: 75 ed jne 50 <gcd64+0x50>
63: 48 89 f8 mov %rdi,%rax
66: 48 d3 e0 shl %cl,%rax
69: c3 retq
6a: 66 0f 1f 44 00 00 nopw 0x0(%rax,%rax,1)
70: 48 d1 ee shr %rsi
73: eb db jmp 50 <gcd64+0x50>
75: 0f 1f 00 nopl (%rax)
78: 48 89 f8 mov %rdi,%rax
7b: 48 89 f7 mov %rsi,%rdi
7e: 48 29 f0 sub %rsi,%rax
81: 48 89 c6 mov %rax,%rsi
84: eb d8 jmp 5e <gcd64+0x5e>
86: 66 2e 0f 1f 84 00 00 nopw %cs:0x0(%rax,%rax,1)
8d: 00 00 00
0000000000000090 <min>:
90: f3 0f 1e fa endbr64
94: 31 c0 xor %eax,%eax
96: 39 f7 cmp %esi,%edi
98: 0f 9f c0 setg %al
9b: 31 fe xor %edi,%esi
9d: f7 d8 neg %eax
9f: 21 f0 and %esi,%eax
a1: 31 f8 xor %edi,%eax
a3: c3 retq
a4: 66 66 2e 0f 1f 84 00 data16 nopw %cs:0x0(%rax,%rax,1)
ab: 00 00 00 00
af: 90 nop
00000000000000b0 <gcd64_use_builtin_ctz>:
b0: f3 0f 1e fa endbr64
b4: 48 89 f8 mov %rdi,%rax
b7: 48 85 ff test %rdi,%rdi
ba: 74 6c je 128 <gcd64_use_builtin_ctz+0x78>
bc: 48 85 f6 test %rsi,%rsi
bf: 74 67 je 128 <gcd64_use_builtin_ctz+0x78>
c1: 31 c9 xor %ecx,%ecx
c3: 31 d2 xor %edx,%edx
c5: f3 48 0f bc cf tzcnt %rdi,%rcx
ca: f3 48 0f bc d6 tzcnt %rsi,%rdx
cf: 31 ff xor %edi,%edi
d1: 39 d1 cmp %edx,%ecx
d3: 41 89 d0 mov %edx,%r8d
d6: 40 0f 9f c7 setg %dil
da: 41 31 c8 xor %ecx,%r8d
dd: 48 d3 e8 shr %cl,%rax
e0: f7 df neg %edi
e2: 44 21 c7 and %r8d,%edi
e5: 31 cf xor %ecx,%edi
e7: 89 d1 mov %edx,%ecx
e9: 48 d3 ee shr %cl,%rsi
ec: eb 0d jmp fb <gcd64_use_builtin_ctz+0x4b>
ee: 66 90 xchg %ax,%ax
f0: 48 89 d6 mov %rdx,%rsi
f3: 48 29 c6 sub %rax,%rsi
f6: 48 85 f6 test %rsi,%rsi
f9: 74 20 je 11b <gcd64_use_builtin_ctz+0x6b>
fb: 31 c9 xor %ecx,%ecx
fd: 48 89 f2 mov %rsi,%rdx
100: f3 48 0f bc ce tzcnt %rsi,%rcx
105: 48 d3 ea shr %cl,%rdx
108: 48 39 c2 cmp %rax,%rdx
10b: 77 e3 ja f0 <gcd64_use_builtin_ctz+0x40>
10d: 48 29 d0 sub %rdx,%rax
110: 48 89 c6 mov %rax,%rsi
113: 48 89 d0 mov %rdx,%rax
116: 48 85 f6 test %rsi,%rsi
119: 75 e0 jne fb <gcd64_use_builtin_ctz+0x4b>
11b: 89 f9 mov %edi,%ecx
11d: 48 d3 e0 shl %cl,%rax
120: c3 retq
121: 0f 1f 80 00 00 00 00 nopl 0x0(%rax)
128: 48 09 f0 or %rsi,%rax
12b: c3 retq
由編譯結果的組合語言可以看出，透過 __builtin_ctz 函數，gcc 得以直接將其編譯成 tzcnt 指令。使得所有計算結尾零數量的部份得以從未知次數的迴圈簡化至接近 constant time 的時間結果。
原先希望透過以搜尋 tzcnt 的 CPI(Cycle per Instruction) 解釋此現象，於網路上搜尋過後發現許多人對於 CPI 的相關資料都是取自Agner Fog. 的 Insturction Table，發現在大多數的處理器當中，tzcnt 的 CPI 大多為 add, sub, shr 等基本運算指令的 2 ~ 4 倍。故推論，若迴圈次數大於2次後，使用迴圈的次數應會較慢一些。
為驗證這部份推論，設計簡易實驗，僅測試 __builtin_ctz 功能與迴圈差異。結果如下
實驗相關程式碼
#include <stdio.h>
#include <stdint.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include "util.h"
int func(unsigned int a){
int shift;
for (shift = 0; !(a & 1); shift++)
a >>= 1;
return shift;
}
int func1(unsigned int a){
return __builtin_ctz(a);
}
typedef int (func_f)(unsigned int);
void test(unsigned int (*a_gen)(void), const char *out_fn) {
unsigned int a;
const int epochs = 1000;
int interval = epochs / 10;
struct timespec start, end;
func_f *funcs[] = {
func,
func1,
};
const int sort_len = sizeof(funcs) / sizeof(func_f *);
time_t measures[epochs][sort_len];
for (int i = 0; i < epochs; i++)
{
a = a_gen();
for (size_t j = 0; j < sort_len; j++)
{
clock_gettime(CLOCK_MONOTONIC, &start);
funcs[j](a);
clock_gettime(CLOCK_MONOTONIC, &end);
measures[i][j] = diff(start, end);
}
if ((i + 1) % interval == 0)
printf("Benchmark %-3d%%\n", 100 * (i + 1) / epochs);
}
FILE *file = fopen(out_fn, "w+");
for (size_t i = 0; i < epochs; i++)
{
for (size_t j = 0; j < sort_len; j++)
{
char *token = j == sort_len - 1 ? j == epochs - 1 ? "" : "\n" : "\t";
fprintf(file, "%ld%s", measures[i][j], token);
}
}
}
unsigned int r1() {
return (unsigned int) 1;
}
unsigned int rmax() {
return (unsigned int) 1 << 31;
}
unsigned int rrand(){
return (unsigned int) rand();
}
int main(void) {
srand(time(0));
printf("RAND_MAX = %d\n", RAND_MAX);
test(r1, "./txt/1.txt");
test(rmax, "./txt/0x80000000.txt");
test(rrand, "./txt/rand.txt");
return 0;
}
由實驗結果可以看出，使用 __builtin_ctz 函數的計算使用時間接近常數，而使用迴圈的計算方式則會根據結尾的零的數量有些許變動，但由於對於一 64 bit 整數而言最多需要計算的迴圈數量就僅有 63 個，所以就使用到這個子函數的程式整體而言使用迴圈仍可算是常數時間的運算。但是用 __builtin_ctz 則可更有效率。
測驗 4
- 解釋上述程式運作原理，並舉出這樣的程式碼用在哪些真實案例中;
- 設計實驗，檢驗 ctz/clz 改寫的程式碼相較原本的實作有多少改進？應考慮到不同的 bitmap density;
- 思考進一步的改進空間;
- 閱讀 Data Structures in the Linux Kernel 並舉出 Linux 核心使用 bitmap 的案例;
運作原理
#include <stddef.h>
size_t improved(uint64_t *bitmap, size_t bitmapsize, uint32_t *out)
{
size_t pos = 0;
uint64_t bitset;
for (size_t k = 0; k < bitmapsize; ++k) {
bitset = bitmap[k];
while (bitset != 0) {
uint64_t t = bitset & (~bitset + 1);
int r = __builtin_ctzll(bitset);
out[pos++] = k * 64 + r;
bitset ^= t;
}
}
return pos;
}
本函數的目的為將在 bitmap 當中位元為 1 的位元編號紀錄於 out 並回傳總數量。
其中在 improved 版本的函數中，k 代表一個 64 位元的 bitset 在整個 bitmap 當中的編號，而 r 為在 bitmap 當中末端零的數量。
而 k * 64 + r 即為從 bitmap 最開頭的位元標記為 0 號開始往後數的編號，若該位元為 1 則將其紀錄於 out 陣列當中。
而透過第 12 行的 bitset ^= t 將最末一個為 1 的位元清零後即可在下一次迭代中取得下一個 1 的位置，便可達到本函數目的。
可以使用到本函數的最基本情況為，對於一個佔用容量龐大但位元當中 0, 1 數量差異較大的 bit array 資料，可將其資料壓縮為各個 1 或 0 的位置，進行儲存，可以達到基本的壓縮容量效果。
測驗 5
- 解釋上述程式碼運作原理，指出其中不足，並予以改進
例如，判斷負號只要寫作 bool isNegative = numerator < 0 ^ denominator < 0;
搭配研讀 The simple math behind decimal-binary conversion algorithms - 在 Linux 核心原始程式碼的 mm/ 目錄 (即記憶體管理) 中，找到特定整數比值 (即 fraction) 和 bitwise 操作的程式碼案例，予以探討和解說其應用場景
運作原理
#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "list.h"
struct rem_node {
int key;
int index;
struct list_head link;
};
static int find(struct list_head *heads, int size, int key)
{
struct rem_node *node;
int hash = key % size;
list_for_each_entry (node, &heads[hash], link) {
if (key == node->key)
return node->index;
}
return -1;
}
char *fractionToDecimal(int numerator, int denominator)
{
int size = 1024;
char *result = malloc(size);
char *p = result;
if (denominator == 0) {
result[0] = '\0';
return result;
}
if (numerator == 0) {
result[0] = '0';
result[1] = '\0';
return result;
}
/* using long long type make sure there has no integer overflow */
long long n = numerator;
long long d = denominator;
/* deal with negtive cases */
if (n < 0)
n = -n;
if (d < 0)
d = -d;
bool sign = (float) numerator / denominator >= 0;
if (!sign)
*p++ = '-';
long long remainder = n % d;
long long division = n / d;
sprintf(p, "%ld", division > 0 ? (long) division : (long) -division);
if (remainder == 0)
return result;
p = result + strlen(result);
*p++ = '.';
/* Using a map to record all of reminders and their position.
* if the reminder appeared before, which means the repeated loop begin,
*/
char *decimal = malloc(size);
memset(decimal, 0, size);
char *q = decimal;
size = 1333;
struct list_head *heads = malloc(size * sizeof(*heads));
for (int i = 0; i < size; i++)
INIT_LIST_HEAD(&heads[i]);
for (int i = 0; remainder; i++) {
int pos = find(heads, size, remainder);
if (pos >= 0) {
while (PPP > 0)
*p++ = *decimal++;
*p++ = '(';
while (*decimal != '\0')
*p++ = *decimal++;
*p++ = ')';
*p = '\0';
return result;
}
struct rem_node *node = malloc(sizeof(*node));
node->key = remainder;
node->index = i;
MMM(&node->link, EEE);
*q++ = (remainder * 10) / d + '0';
remainder = (remainder * 10) % d;
}
strcpy(p, decimal);
return result;
}
作答
PPP=pos--MMM=list_add_tailEEE=&heads[remainder % size]
首先 fractionToDecimal 函數最一開始先進行較簡單的資料處理，當分子或分母為零時回傳值皆是固定所以先行進行處理。接著由於負數的負號必須紀錄於回傳值的最一開頭，故將正負號進行處理，並在往後運算除去正負號影響因素。其中注意由於 -INT_MIN = INT_MIN 的特性，有可能直接進行負號計算會產生 overflow ，故先將其轉為 long long 再進行計算。
而接著便可以藉由整除與否對接下來的結果進行計算，並且先將整數部分紀錄至結果，最後再考慮小數部分結果。
而在小數部分，為處理循環小數部分，必須有方法紀錄循環小數的開頭。為處理這個部分，可以透過循環小數的特性–由相同的餘數開始除–來進行計算。
透過此一特性，基於長除法的計算方式，紀錄每一次除法運算後的餘數以及其相應的小數位置，以便往後尋找。
為此，建立一個雜湊表 (hash map)，以每次的餘數為 hash key (暫時找不到一個合適的中文翻譯，先以這個方式紀錄)，紀錄其相對應的小數位數。對於每一個往後出現的餘數，皆可以先搜尋其是否出現過。以便得到小數的循環部分。
以上程式碼未被註解掉的部分為較複雜的小數部分，其中 decimal 紀錄小數點後的數字。14至17行的迴圈進行雜湊表的初始化。19行的迴圈開始進行長除法計算，由 37, 38 行可以看到對於每一個剩下的餘數皆進行乘以 10 再除一次取餘數的長除法，而 20 ~ 35 行間即是對於餘數的紀錄。
而在 20 ~ 30 行之間為進行當餘數已開始重複時的判斷，由 decimal 開始第 pos 位為循環小數的開始，便使用 22, 23 行進行循環小數前的部分複製，加上題目需求的括號後再複製剩餘的循環小數部分，便可得到回傳的結果。而 31 ~ 35 行間為將餘數紀錄至雜湊表的過程。若除法結果並非一循環小數，則會在最後脫離迴圈，使用最後 41 行的方式複製至回傳結果，變達成此函數的需求。
功能改進
以上程式碼有些許部分有些不足，以列點方式於底下呈現。
- 在函數中動態配置的記憶體除了回傳的字串外皆應釋放，以避免記憶體洩漏問題。應在處理小數部分中兩個可能的回傳處前加上外加的釋放記憶體函數。
- 99 行中的
strcpy函數由於缺乏邊界確認，故被視為很可能會產生 buffer overflow 的不安全函數，應避免使用。應將其更改為逐一字元拷貝迴圈式複製或者使用strncpy將其替換掉。
可能的程式碼為或是while (*decimal) *p++ = *decimal++; *p = '\0';strncpy(p, decimal, 1024 - strlen(result)); - 程式碼中的字串長度的表示稍嫌不清楚，若能將其定義為巨集，便可以更清楚理解該常數意義
可能修改方式為甚至可將前半部分 /* ... */ + #define RESULT_LEN 1024 /* length of result is guaranteed to be not over 1000 */ + #define HASH_LEN 1333 /* ... */ - int size = 1024; + int size = RESULT_LEN; /* ... */ - size = 1333; + size = HASH_LEN; /* ... */size皆取代為RESULT_LEN，後半部分計算小數部分中的size皆取代為HASH_LEN以減少模糊的變數意義。 - 由於處理器對於浮點數運算大多較整數運算需要更多時間，且判斷正負號的部分可以直接透過分子分母的正負號得到，在該部分可簡化為整數運算。
可能修改方式為此方式透過直接判斷分子分母同號與否來決定結果的正負號。- bool sign = (float) numerator / denominator >= 0; - if (!sign) + if ((numerator ^ denominator) >> 31) *p++ = '-'; - 本題當中
MMM部分巨集，可以從其上下文得知需求為將新製作出之rem_node加入雜湊表當中，而由於本實作方式使用 Linux 風格的 doubly-linked-list，對於增加於頭尾兩端花費皆為常數時間。而應將其加在頭或尾端則是直得探討的部分。若是對於大部分的循環小數部分，循環部分的長度會大於未循環的部分，則應將新產生的節點加入雜湊表各鍊的尾端，若小於則應加於尾端。但巨觀來看，雜湊表的性質來說，一個理想的雜湊表應可以達到
測驗 6
- 解釋上述程式碼運作原理
- 在 Linux 核心原始程式碼中找出
__alignof__的使用案例 2 則，並針對其場景進行解說 - 在 Linux 核心源使程式碼找出
ALIGN,ALIGN_DOWN,ALIGN_UP等巨集，探討其實作機制和用途，並舉例探討 (可和上述第二點的案例重複)。思索能否對 Linux 核心提交貢獻，儘量共用相同的巨集
運作原理
/*
* ALIGNOF - get the alignment of a type
* @t: the type to test
*
* This returns a safe alignment for the given type.
*/
#define ALIGNOF(t) \
((char *)(&((struct { char c; t _h; } *)0)->_h) - (char *)0)
本巨集目的在求得資料型態 t 於記憶體中會需要對齊幾個位元。透過 struct 安排的結果來進行計算。
先將巨集大致拆解後可以看出其結構為
((char *)Y - (char *)X)
故為透過兩個指標相差的距離以計算其對齊方式的技巧。
接著在細看 Y 的部分，可以看到其透過一個 struct 排列兩個元素，由此結構可以推論欲使用的策略如下圖所示
對於任意一個資料類型與 char c 放入一個結構中時，記憶體會呈現上圖的排列模式，若 struct 開頭位於 addr 位址時且資料類型 t 需對其特定數值時，在元素 c 與 _h 當中會加入填充用的記憶體空間 padding 以達到對齊效果。透過這個方式，當 addr 被設定為 0 時，元素 _h 的記憶體位址即為其對齊所需的位元數。
因此，在上方 Y 的部分取得 _h 元素在 struct 中的相對位址，再減去其開頭 (也就是 0) 即可得到其對齊位元數。
測驗 7
- 解釋上述程式運作原理並評估 naive.c 和 bitwise.c 效能落差
- 避免 stream I/O 帶來的影響，可將 printf 更換為 sprintf
- 分析 Faster remainders when the divisor is a constant: beating compilers and libdivide 一文的想法 (可參照同一篇網誌下方的評論)，並設計另一種 bitmask，如「可被 3 整除則末位設為 1」「可被 5 整除則倒數第二位設定為 1」，然後再改寫 bitwise.c 程式碼，試圖運用更少的指令來實作出 branchless;
- 研讀 The Fastest FizzBuzz Implementation 及其 Hacker News 討論，提出 throughput (吞吐量) 更高的 Fizzbuzz 實作
- 解析 Linux 核心原始程式碼 kernel/time/timekeeping.c 裡頭涉及到除法運算的機制，探討其快速除法的實作 (注意: 你可能要對照研讀 kernel/time/ 目錄的標頭檔和程式碼)
過程中，你可能會發現可貢獻到 Linux 核心的空間，請充分討論
運作原理
static inline bool is_divisible(uint32_t n, uint64_t M)
{
return n * M <= M - 1;
}
static uint64_t M3 = UINT64_C(0xFFFFFFFFFFFFFFFF) / 3 + 1;
static uint64_t M5 = UINT64_C(0xFFFFFFFFFFFFFFFF) / 5 + 1;
int main(int argc, char **argv)
{
for (size_t i = 1; i <= 100; i++) {
uint8_t div3 = is_divisible(i, M3);
uint8_t div5 = is_divisible(i, M5);
unsigned int length = (2 << div3) << div5;
char fmt[9];
strncpy(fmt, &"FizzBuzz%u"[(9 >> div5) >> (KK3)], length);
fmt[length] = '\0';
printf(fmt, i);
printf("\n");
}
return 0;
}
