Try   HackMD

2022q1 Homework2 (quiz2)

contributed by < NOVBobLee >

作業要求

測驗 1

對兩個無號 32 位元整數取平均值。

/* method 1 */
#include <stdint.h>
uint32_t average(uint32_t a, uint32_t b)
{
    return (a + b) / 2;
}

方法一會有 overflow 的問題,我們將 a, b 都代入

42949672950xFFFF'FFFF ),結果會得到
2147483647
,是代入值的一半,顯然跟我們期望的平均值不同(期望輸出為
4294967295
),以下用 uint8_t 為例。

uint8_t 二進位 備註
a 1111'1111
b 1111'1111
a + b (1)'1111'1111 溢位
(a + b) / 2 0111'1111

無號整數的溢位會發生在相加超過上限值或者相減變為負數的時候,若我們知道誰的值比較大,那可以使用方法二,避免以上兩者的發生。

l+h2=2ll+h2=l+hl21. hl02. Suppose that  h,luint32_t s.t. α=l+hl2>max(uint32_t)Then, l+hl2=l+h2>max(uint32_t)But, l+h2 has upper bound max(uint32_t) , happens when h=l=max(uint32_t)

/* method 2 */
#include <stdint.h>
uint32_t average(uint32_t low, uint32_t high)
{
    return low + (high - low) / 2;
}

但需要事先知道誰大誰小,而方法三就避免掉減法,也就不需要比大小了。在整數上做除法,得到的是商數,所以還需要對兩者的餘數做處裡,從真值表可以看出該處理 EXP1 應該為 a & b & 1 (用 &1 限定

20 位置的位元)。

a+b2=a2+b2

a & 1 b & 1 (a & 1 + b & 1) / 2
0 0 0
0 1 0
1 0 0
1 1 1
/* method 3 */
#include <stdint.h>
uint32_t average(uint32_t a, uint32_t b)
{
    return (a >> 1) + (b >> 1) + (EXP1);
}

那有沒有可能把方法三再縮短,可以觀察 a + b 的位元計算,將 a 自己拆開成與 b 相同的部份跟不同的部份,同樣對 b 做拆開。

a = a & b + (a ^ b) & a
b = a & b + (a ^ b) & b
    ^^^^^   ^^^^^^^^^^^

然後相加 ab ,可以看到相同的部份就是兩倍,而不同的部份則剛好是互補。

a + b = (a & b) << 1 + (a ^ b) & a + (a ^ b) & b
      = (a & b) << 1 + (a ^ b)

相加可以看到左移一位有溢位的可能,不過我們還會再做一個除以

2 的右移,就避開了溢位的危險,而方法四即使用此運算方式。

(a + b) / 2 = ((a & b) << 1 + (a ^ b)) >> 1 // 示意,左移和右移要先相消
            = (a & b) + (a ^ b) >> 1
/* method 4 */
#include <stdint.h>
uint32_t average(uint32_t a, uint32_t b)
{
    return (EXP2) + ((EXP3) >> 1);
}

延伸問題:

  1. 解釋上述程式碼運作的原理
  2. 比較上述實作在編譯器最佳化開啟的狀況,對應的組合語言輸出,並嘗試解讀 (可研讀 CS:APP 第 3 章)
  3. 研讀 Linux 核心原始程式碼 include/linux/average.h,探討其 Exponentially weighted moving average (EWMA) 實作,以及在 Linux 核心的應用

移動平均(Moving average),又稱滾動平均值、滑動平均,在統計學中是種藉由建立整個資料集合中不同子集的一系列平均數,來分析資料點的計算方法。
移動平均通常與時間序列資料一起使用,以消除短期波動,突出長期趨勢或周期。短期和長期之間的閾值取決於應用,移動平均的參數將相應地設置。例如,它通常用於對財務數據進行技術分析,如股票價格、收益率或交易量。它也用於經濟學中研究國內生產總值、就業或其他宏觀經濟時間序列。

以下為比較有開最佳化和沒開的編譯組合語言結果,使用 AT&T 語法。

方法一

沒有開啟最佳化結果:

; gcc -S averages.c
; uint32_t average1(uint32_t a, uint32_t b) { return (a + b) / 2; }
average1:
.LFB0:
	.cfi_startproc
	endbr64
	pushq	%rbp             ; push frame pointer into stack
	.cfi_def_cfa_offset 16
	.cfi_offset 6, -16
	movq	%rsp, %rbp       ; frame pointer = stack pointer
	.cfi_def_cfa_register 6
	movl	%edi, -4(%rbp)   ; write 1st argument into stack
	movl	%esi, -8(%rbp)
	movl	-4(%rbp), %edx   ; read from stack, assign the value to register
	movl	-8(%rbp), %eax
	addl	%edx, %eax       ; %eax += %edx , i.e. (a + b)
	shrl	%eax             ; %eax >>= 1 , right shift
	popq	%rbp             ; pop frame pointer from stack
	.cfi_def_cfa 7, 8
	ret                     ; return (%eax is result register)
	.cfi_endproc
  1. endbr64 是什麼?
    CET (Intel's Control-Flow Enforcement Technology) , GCC 8 開始加入支援,用來防止 ROP, JOP/COP 攻擊。
    來源 What does the endbr64 instruction actually do?
  2. %rbp, %rsp 是什麼?

    Intel® 64 and IA-32 Architectures Software Developer’s Manual
    If the nesting level is 0, the processor pushes the frame pointer from the BP/EBP/RBP register onto the stack, copies the current stack pointer from the SP/ESP/RSP register into the BP/EBP/RBP register, and loads the SP/ESP/RSP register with the current stack-pointer value minus the value in the size operand. For nesting levels of 1 or greater, the processor pushes additional frame pointers on the stack before adjusting the stack pointer. These additional frame pointers provide the called procedure with access points to other nested frames on the stack.
    Releases the stack frame set up by an earlier ENTER instruction. The LEAVE instruction copies the frame pointer (in the EBP register) into the stack pointer register (ESP), which releases the stack space allocated to the stack frame. The old frame pointer (the frame pointer for the calling procedure that was saved by the ENTER instruction) is then popped from the stack into the EBP register, restoring the calling procedure’s stack frame.

  3. .cfi_* 是什麼?
    CFI (call frame information) directives 是屬於 GNU AS extension , assembler directives 之一,可以管理 call frames ,在有些架構裡需要使用這些 CFI directives 來處理錯誤。
    來源一 What are CFI directives in Gnu Assembler (GAS) used for?
    來源二 CFI directives - Using as

開啟最佳化結果(之後都將 .cfi_*, endbr64 去除),其結果 -O0, -O1, -O2 皆相同:

; gcc -S averages.c -O0
; uint32_t average1(uint32_t a, uint32_t b) { return (a + b) / 2; }
average1:
	leal	(%rdi,%rsi), %eax ; sum of 1st and 2nd argument, i.e. %eax = a + b
	shrl	%eax              ; right shift
	ret

從兩者可以看出指令數的差別,還有就是開了最佳化後,連 stack 都不見了XD 沒開最佳化的時候,他會先把引數存進 stack 裡,再從 stack 讀回到 register 裡面,不知道是不是因為 argument register 不能拿來使用才這樣做?(可以使用,見方法三使用最佳化編譯結果)

這邊有點疑問,在 GCC Optimize Options 是說預設就是 -O0 ,但 -O0 結果跟沒加 flag 的結果不同。

方法二

沒開最佳化,可以看到第 7 到 10 行,剛開始使用 %eax 計算,最後又換到 %edx 上,其實可以直接使用 %edx 即可。

; uint32_t average2(uint32_t low, uint32_t high) { return low + (high - low) / 2; } average2: pushq %rbp movq %rsp, %rbp movl %edi, -4(%rbp) ; write 1st argument into stack movl %esi, -8(%rbp) movl -8(%rbp), %eax ; read from stack, assign the value to register subl -4(%rbp), %eax ; %eax -= -4(%rbp), i.e. %eax = high - low shrl %eax ; %eax >>= 1 movl %eax, %edx movl -4(%rbp), %eax addl %edx, %eax ; %eax += %edx, i.e. %eax += low popq %rbp ret

開最佳化之後( -O0, -O1, -O2 三者結果相同), stack 也不見了。

; uint32_t average2(uint32_t low, uint32_t high) { return low + (high - low) / 2; }
average2:
	movl	%esi, %eax    ; copy 2nd argument to %eax
	subl	%edi, %eax    ; %eax -= %edi, i.e. %eax = high - low
	shrl	%eax          ; %eax >>= 1
	addl	%edi, %eax    ; %eax += %edi, i.e. %eax += low
	ret

方法三

沒開最佳化一樣會使用 stack ,跟開最佳化的結果比較後, stack 還是多餘的動作。

; uint32_t average3(uint32_t a, uint32_t b) { return (a >> 1) + (b >> 1) + (a & b & 1); } average3: pushq %rbp movq %rsp, %rbp movl %edi, -4(%rbp) ; write 1st argument into stack movl %esi, -8(%rbp) movl -4(%rbp), %eax ; %eax = a shrl %eax ; %eax >>= 1 movl %eax, %edx movl -8(%rbp), %eax ; %eax = b shrl %eax ; %eax >>= 1 addl %eax, %edx ; %edx = (a >> 1) + (b >> 1) movl -4(%rbp), %eax ; %eax = a andl -8(%rbp), %eax ; %eax &= b andl $1, %eax ; %eax &= 1 addl %edx, %eax ; %eax += %edx , i.e. (a >> 1) + (b >> 1) + (a & b & 1) popq %rbp ret

開最佳化後( -O0, -O1, -O2 結果一樣),沒有使用 stack ,而 register 的使用上比沒開最佳化的更靈活,像是開最佳化後的第 3 到 7 行對比沒開的第 7 至 12 行,還有開最佳化的第 8 到 10 行對比沒開的第 13 到 16 行。

; uint32_t average3(uint32_t a, uint32_t b) { return (a >> 1) + (b >> 1) + (a & b & 1); } average3: movl %edi, %eax ; %eax = a shrl %eax ; %eax >>= 1 movl %esi, %edx ; %edx = b shrl %edx ; %edx >>= 1 addl %edx, %eax ; %eax += %edx , i.e. (a >> 1) + (b >> 1) andl %esi, %edi ; b &= a ,use argument register andl $1, %edi ; b &= 1 addl %edi, %eax ; add together ret

方法四

沒開最佳化有使用 stack ,在 register 的選擇上會優先使用 %eax ,所以出現第 9 行的交換。

; uint32_t average4(uint32_t a, uint32_t b) { return (a & b) + ((a ^ b) >> 1); } average4: pushq %rbp movq %rsp, %rbp movl %edi, -4(%rbp) ; write argument register into stack movl %esi, -8(%rbp) movl -4(%rbp), %eax ; read from stack, %eax = a andl -8(%rbp), %eax ; %eax &= b , i.e. %eax = a & b movl %eax, %edx movl -4(%rbp), %eax ; %eax = a xorl -8(%rbp), %eax ; %eax ^= b , i.e. %eax = a ^ b shrl %eax ; %eax >>= 1 addl %edx, %eax ; %eax = (a & b) + ((a ^ b) >> 1) popq %rbp ret

這次開啟最佳化 -O2 有些變化,與 -O0, -O1 不一樣,兩個都沒有使用 stack 。

; -O0, -O1
; uint32_t average4(uint32_t a, uint32_t b) { return (a & b) + ((a ^ b) >> 1); }
average4:
	movl	%edi, %eax    ; %eax = a
	xorl	%esi, %eax    ; %eax ^= b
	shrl	%eax          ; %eax >>= 1
	andl	%esi, %edi    ; a &= b
	addl	%edi, %eax    ; add together
	ret

雖然做的事跟 -O0, -O1 一樣,不過 -O2 在 register 上的使用上更靈活,第 4 行將 argument register %edi 的值複製進 %eax 之後就拿來存 a & b 的結果了,換句話說,運用 register 的手段比之前的多。

; -O2 ; uint32_t average4(uint32_t a, uint32_t b) { return (a & b) + ((a ^ b) >> 1); } average4: movl %edi, %eax ; %eax = a andl %esi, %edi ; a &= b xorl %esi, %eax ; %eax ^= b shrl %eax ; %eax >>= 1 addl %edi, %eax ; add together ret

看完之後做比較,方法一二的使用指令數比較少,不過有 overflow 的可能或是需要 a, b 的大小關係。而方法三四以比較 general 的方法,且方法四的使用指令數比較少。

(TODO: 測試)

Linux EWMA 實作與應用

EWMA (Exponentially Weighted Moving Average) 是一種特殊的取平均方式,在一資料集中不同的子集合各取其平均,從名稱中可以看到移動,意指不同的子集,而其平均方式會使用指數權重。給一資料集合序列

Y ,而其 EWMA
S
計算方式如下:

St={Y0if t=0αYt+(1α)St1if t>0

其計算非常簡單,最初的 EWMA

S0 為資料本身
Y0
,因為在之前還沒有計算 EWMA 。而之後每增加一筆新的資料
Yt
會乘上一個權重
α
,再加上之前計算的 EWMA
St1
乘上
1α
,亦即以前的資料也會影響到現在的 EWMA ,可參考下面的圖例(
α=0.2
)。

0123t050100150200EWMA S(t)1. Y(0)=200Data Y(t)
0123t050100150200EWMA S(t)1. Y(0)=2002. Y(1)=150Data Y(t)
0123t050100150200EWMA S(t)1. Y(0)=2002. Y(1)=1503. Y(2)=180Data Y(t)
0123t050100150200EWMA S(t)1. Y(0)=2002. Y(1)=1503. Y(2)=1804. Y(3)=100Data Y(t)

其中

α 又稱 smoothing factor ,取值介於
[0,1]
之間,越低則舊的資料的影響越大,相反地,越高則是新加入的資料影響越大。通常會取比較低的值,使得新加入的資料不會讓現有的變化趨勢波動過大,有 smoothing 的作用。

以下為 Linux 的 EWMA 實作,位置於 include/linux/average.h 。要使用 EWMA 需要先用此 macro 將 EWMA 專用的結構和函式建立起來,不同名稱的 EWMA 會有不同的對應結構名稱與函式。第一個引數 name 為 EWMA 名稱,用來識別不同的 EWMA 。第二個參數 _precision 為 EWMA 的精密度是在小數第幾位(二進位)。第三個參數 _weight_rcp 是權重的倒數,而且要為

2 的次方,直接看數學式
St=1wYt+(11w)St1 , t>0
,使用倒數的原因是權重
1w
範圍可以被限定在
(0,1]
之間了。

/* * Exponentially weighted moving average (EWMA) * * This implements a fixed-precision EWMA algorithm, with both the * precision and fall-off coefficient determined at compile-time * and built into the generated helper funtions. * * The first argument to the macro is the name that will be used * for the struct and helper functions. * * The second argument, the precision, expresses how many bits are * used for the fractional part of the fixed-precision values. * * The third argument, the weight reciprocal, determines how the * new values will be weighed vs. the old state, new values will * get weight 1/weight_rcp and old values 1-1/weight_rcp. Note * that this parameter must be a power of two for efficiency. */ #define DECLARE_EWMA(name, _precision, _weight_rcp) \ struct ewma_##name { \ unsigned long internal; \ }; \ static inline void ewma_##name##_init(struct ewma_##name *e) \ { \ BUILD_BUG_ON(!__builtin_constant_p(_precision)); \ BUILD_BUG_ON(!__builtin_constant_p(_weight_rcp)); \ /* \ * Even if you want to feed it just 0/1 you should have \ * some bits for the non-fractional part... \ */ \ BUILD_BUG_ON((_precision) > 30); \ BUILD_BUG_ON_NOT_POWER_OF_2(_weight_rcp); \ e->internal = 0; \ } \

所有建立的函式名稱和結構名稱都會接上 name 這個名稱,用以區分不同的 EWMA 。最開始會先宣告 struct ewma_name (之後用 name 代表接上的部份),裡面只有一個成員 internal ,紀錄我們的 EWMA 。再來是自己宣告結構變數,然後使用 ewma_name_init 初始化該變數。

在函式中可以看到很多 BUILD_BUG_ON ,那是 Linux 設計在 compile-time 檢查錯誤的 macro 。這邊會檢查 _precision, _weight_rcp 是否為常數表示式, _precision 是否在

30 以內, _weight_rcp 是否為
2
的冪(且不可為零),若沒有符合,在編譯的時候就會跳出錯誤。

這邊我有個問題,為何要三個函式都重複檢查,何不在其中一個函數裡作所有檢查就好?我目前猜測是跳出錯誤訊息的時候,可以明確指出哪些函式會使用到這些含有錯誤的表示式,在除錯時可以比較全盤的掌握錯誤所在。

static inline unsigned long \ ewma_##name##_read(struct ewma_##name *e) \ { \ BUILD_BUG_ON(!__builtin_constant_p(_precision)); \ BUILD_BUG_ON(!__builtin_constant_p(_weight_rcp)); \ BUILD_BUG_ON((_precision) > 30); \ BUILD_BUG_ON_NOT_POWER_OF_2(_weight_rcp); \ 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; \ \ BUILD_BUG_ON(!__builtin_constant_p(_precision)); \ BUILD_BUG_ON(!__builtin_constant_p(_weight_rcp)); \ BUILD_BUG_ON((_precision) > 30); \ BUILD_BUG_ON_NOT_POWER_OF_2(_weight_rcp); \ \ WRITE_ONCE(e->internal, internal ? \ (((internal << weight_rcp) - internal) + \ (val << precision)) >> weight_rcp : \ (val << precision)); \ }

剩下兩個函式一個拿來看現在 EWMA ,另一個拿來加入新資料更新 EWMA 用的。 ewma_name_add 會加入一筆新資料,更新我們的 EWMA ,可以看到 _weight_rcp

2 的冪是為了使用位元的左移右移,會經過 ilog2 轉換成次方數 weight_rcp

S0=Y02p
St=((St12wSt1)+Yt2p)2w , t>0=(112w)St1+12wYt2p , t>0

而精密度的實作是將每一筆資料加入的時候都預先左移,這樣紀錄在 internal 裡的 EWMA 其實是實際值的

2precision 倍,這樣多出來的右側位數就可以拿來計算小數的部份,也因為如此,查看 EWMA 的時候,使用 ewma_name_read 裡會做一個右移,將 internal 還原成實際的 EWMA 。

在 Linux 裡很多網路相關的 driver 有使用 EWMA ,有算晶片溫度、資料吞吐量、訊號強度 (rssi) 等平均,而其他像是網路協議和 GPU 的驅動也有在使用。從搜尋結果來看,這巨集設計也是不錯,可以從 name 看出他的作用為何。

drivers/net/virtio_net.c
54:DECLARE_EWMA(pkt_len, 0, 64)
drivers/net/wireless/ath/ath5k/ath5k.h
1255:DECLARE_EWMA(beacon_rssi, 10, 8)
drivers/net/wireless/realtek/rtw88/main.h
642:DECLARE_EWMA(tp, 10, 2);
738:DECLARE_EWMA(rssi, 10, 16);
1497:DECLARE_EWMA(thermal, 10, 4);
1543:DECLARE_EWMA(evm, 10, 4);
1544:DECLARE_EWMA(snr, 10, 4);
drivers/net/wireless/realtek/rtw89/core.h
860:DECLARE_EWMA(tp, 10, 2);
1832:DECLARE_EWMA(rssi, 10, 16);
2384:DECLARE_EWMA(thermal, 4, 4);
drivers/net/wireless/ralink/rt2x00/rt2x00.h
258:DECLARE_EWMA(rssi, 10, 8)
drivers/net/wireless/mediatek/mt76/mt76x02_mac.h
34:DECLARE_EWMA(pktlen, 8, 8);
drivers/net/wireless/mediatek/mt76/mt76.h
232:DECLARE_EWMA(signal, 10, 8);
drivers/net/wireless/mediatek/mt7601u/mt7601u.h
134:DECLARE_EWMA(rssi, 10, 4);
drivers/net/wireless/intel/iwlwifi/mvm/mvm.h
560:DECLARE_EWMA(rate, 16, 16)
drivers/net/wireless/mediatek/mt76/mt7921/mt7921.h
112:DECLARE_EWMA(rssi, 10, 8);
net/batman-adv/types.h
578:DECLARE_EWMA(throughput, 10, 8)
net/mac80211/sta_info.h
126:DECLARE_EWMA(avg_signal, 10, 8)
373:DECLARE_EWMA(mesh_fail_avg, 20, 8)
374:DECLARE_EWMA(mesh_tx_rate_avg, 8, 16)
433:DECLARE_EWMA(signal, 10, 8)
net/mac80211/ieee80211_i.h
435:DECLARE_EWMA(beacon_signal, 4, 4)
drivers/gpu/drm/drm_self_refresh_helper.c
56:DECLARE_EWMA(psr_time, 4, 4)
drivers/gpu/drm/i915/gt/intel_context_types.h
24:DECLARE_EWMA(runtime, 3, 8);
drivers/gpu/drm/i915/gt/intel_engine_types.h
125:DECLARE_EWMA(_engine_latency, 6, 4)

測驗 2

以無分支方式求出最大值。

#include <stdint.h>
uint32_t max(uint32_t a, uint32_t b)
{
    return a ^ ((EXP4) & -(EXP5));
}

先來看看 XOR 運算的其中一些特性:

  1. (a ^ b) ^ a = (b ^ a) ^ a = b ^ (a ^ a) = b
  2. (a ^ b) ^ b = a ^ (b ^ b) = a
  3. a ^ 0 = a

以上特性在 max 的引數裡兩者選其一可以使用,當要替換 a 變成 b ,那麼 a 後面可以接 a ^ b ,而不想替換的話則是接 0 ,這樣看起來 & 後面接一個 bitmask 可以辦到。所以當 max 為 a 時需要 mask 為 0x0 ,而當 max 為 b 時則需要 mask 為 0xFFFF'FFFF ,而題目有提示 EXP5 要用 relational operator ,那應該填 a > b 。當 a > b 為真,則 (a ^ b) & (-1) 變成 (a ^ b) & 0xFFFF'FFFF ,得到 (a ^ b) ,最後 a ^ (a ^ b) 等於 a 。當 a <= b 同理。其中有用到 integer promotions 。

先看 relational operator 結果的型態,是 int

N1256 (6.5.8-6)
Each of the operators < (less than), > (greater than), <= (less than or equal to), and >= (greater than or equal to) shall yield 1 if the specified relation is true and 0 if it is false. The result has type int.

再來看 integer promotions ,由於 (a ^ b) 結果是 unsigned ,而 -1int ,根據 integer promotions 的規則, -1 會改變型態成 unsigned ,而其值則為 0xFFFF'FFFF ,沒有改變。

N1256 (8.3.1.8-1)
[] Otherwise, the integer promotions are performed on both operands. Then the following rules are applied to the promoted operands:

If both operands have the same type, then no further conversion is needed.
Otherwise, if both operands have signed integer types or both have unsigned integer types, the operand with the type of lesser integer conversion rank is converted to the type of the operand with greater rank.
Otherwise, if the operand that has unsigned integer type has rank greater or equal to the rank of the type of the other operand, then the operand with signed integer type is converted to the type of the operand with unsigned integer type.
Otherwise, if the type of the operand with signed integer type can represent all of the values of the type of the operand with unsigned integer type, then the operand with unsigned integer type is converted to the type of the operand with signed integer type.
Otherwise, both operands are converted to the unsigned integer type corresponding to the type of the operand with signed integer type.

Integer Conversion Rank

N1256 (6.3.1.1-1)
Every integer type has an integer conversion rank defined as follows:
— No two signed integer types shall have the same rank, even if they have the same representation.
— The rank of a signed integer type shall be greater than the rank of any signed integer type with less precision.
— The rank of long long int shall be greater than the rank of long int, which shall be greater than the rank of int, which shall be greater than the rank of short int, which shall be greater than the rank of signed char.
— The rank of any unsigned integer type shall equal the rank of the corresponding signed integer type, if any.
— The rank of any standard integer type shall be greater than the rank of any extended integer type with the same width.
— The rank of char shall equal the rank of signed char and unsigned char.
— The rank of _Bool shall be less than the rank of all other standard integer types.
— The rank of any enumerated type shall equal the rank of the compatible integer type.
— The rank of any extended signed integer type relative to another extended signed integer type with the same precision is implementation-defined, but still subject to the other rules for determining the integer conversion rank.
— For all integer types T1, T2, and T3, if T1 has greater rank than T2 and T2 has greater rank than T3, then T1 has greater rank than T3.

延伸問題:

  1. 解釋上述程式碼運作的原理
  2. 針對 32 位元無號/有號整數,撰寫同樣 branchless 的實作
  3. Linux 核心也有若干 branchless / branch-free 的實作,例如 lib/sort.c:
    /*
    ​​​​ * 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;
    ​​​​}
    
    請在 Linux 核心原始程式碼中,找出更多類似的實作手法。請善用 git log 檢索。

程式碼 uint32_t 也可以換成 int32_t ,其無號有號都適用。若要算 min 則將 bitmask 的生成條件改變即可,即 a < b

測驗 3

使用 Euclidean algorithm 或稱輾轉相除法,求出 64 位元 GCD (greatest common divisor)。

#include <stdint.h>
uint64_t gcd64(uint64_t u, uint64_t v)
{
    if (!u || !v) return u | v;
    while (v) {                               
        uint64_t t = v;
        v = u % v;
        u = t;
    }
    return u;
}

GCD 演算法(又稱 binary GCD alogrithm):

  1. If both
    x
    and
    y
    are
    0
    , gcd is zero
    gcd(0,0)=0
    .
  2. gcd(x,0)=x
    and
    gcd(0,y)=y
    because everything divides
    0
    .
  3. If
    x
    and
    y
    are both even,
    gcd(x,y)=2gcd(x2,y2)
    because
    2
    is a common divisor. Multiplication with
    2
    can be done with bitwise shift operator.
  4. If
    x
    is even and
    y
    is odd,
    gcd(x,y)=gcd(x2,y)
    .
  5. On similar lines, if
    x
    is odd and
    y
    is even, then
    gcd(x,y)=gcd(x,y2)
    . It is because
    2
    is not a common divisor.

以下用可位元操作的改寫。可以看出下面改寫比原算法多了 v/=2, u/=2 的動作,其對應的是 GCD 演算法的第 3 到 5 點,相較於使用比較費力的取模操作,改使用右移與減法更適當。

#include <stdint.h>
uint64_t gcd64(uint64_t u, uint64_t v)
{
    if (!u || !v) return u | v; /* GCD 演算法第 1, 2 點 */
    int shift;
    for (shift = 0; !((u | v) & 1); shift++) {
        u /= 2, v /= 2; /* GCD 演算法第 3 點 */
    }
    while (!(u & 1))
        u /= 2; /* GCD 演算法第 4(5) 點 */
    do {
        while (!(v & 1))
            v /= 2; /* GCD 演算法第 4(5) 點 */
        if (u < v) {
            v -= u;
        } else {
            uint64_t t = u - v;
            u = v;
            v = t;
        }
    } while (COND);
    return RET;
}

在 do-while 的區塊,其本質跟原算法相同,即計算餘數,再拿餘數當除數,換句話說就是輾轉相除法,停止條件跟原算法相同,條件為整除,餘數為零的時候,所以 COND 應填為 v 。而回傳值應該也跟原算法回傳商數一樣,不過我們有先做 GCD 演算法第 3 點,先提出公因數

2shift ,所以要再乘回到我們的商數 u 上,即 RET 應為 u << shift。雖說是輾轉相除,但實作是使用反覆的減法代替除法(意義上跟原本的輾轉相除法的商數、餘數相同)。

延伸問題:

  1. 解釋上述程式運作原理;
  2. 在 x86_64 上透過 __builtin_ctz 改寫 GCD,分析對效能的提升;
  3. Linux 核心中也內建 GCD (而且還不只一種實作),例如 lib/math/gcd.c,請解釋其實作手法和探討在核心內的應用場景。

使用 GCC built-in functions__builtin_ctz 作改寫。

Built-in Function: int __builtin_ctz (unsigned int x)

Returns the number of trailing 0-bits in x, starting at the least significant bit position. If x is 0, the result is undefined.

#include <stdint.h>
uint64_t min(uint64_t a, uint64_t b)
{
    return a ^ ((a ^ b) & -(a < b));
}
uint64_t gcd64(uint64_t u, uint64_t v)
{
    if (!u || !v) return u | v; /* GCD 演算法第 1, 2 點 */
    /* GCD 演算法第 3 點 */
    int shift = min(__builtin_ctz(u), __builtin_ctz(v));
    u >>= shift;
    v >>= shift;
    /* GCD 演算法第 4(5) 點 */
    u >>= __builtin_ctz(u);
    do {
        /* GCD 演算法第 4(5) 點 */
        v >>= __builtin_ctz(v);
        if (u < v) {
            v -= u;
        } else {
            uint64_t t = u - v;
            u = v;
            v = t;
        }
    } while (v);
    return u << shift;
}

在 Linux 內建的 GCD 實作使用大量 bitwise 操作,有分使用 __ffs 的版本跟沒有使用的版本,有使用 __ffs 的版本會比 even/odd algorithm (binary GCD algorithm) 的速度更快。

#if !defined(CONFIG_CPU_NO_EFFICIENT_FFS) /* If __ffs is available, the even/odd algorithm benchmarks slower. */ /** * gcd - calculate and return the greatest common divisor of 2 unsigned longs * @a: first value * @b: second value */ unsigned long gcd(unsigned long a, unsigned long b) { unsigned long r = a | b; if (!a || !b) return r; b >>= __ffs(b); if (b == 1) return r & -r; for (;;) { a >>= __ffs(a); if (a == 1) return r & -r; if (a == b) return a << __ffs(r); if (a < b) swap(a, b); a -= b; } }

第一個是使用 __ffs 的版本, ffs (find first set) 是在找 LSB (least significant bit) 的位數,有分硬體支援函式庫支援,這邊指的是硬體。

方法本質也是輾轉相除法,不過有些改良,第 14 行一樣是 GCD 演算法第 1, 2 點,不過這個 r 就有意思了, r & -r 會得到 LSB ,在我們的算法中,剛開始有一步是算公因數

2shift ,這邊使用 r & -r 代替掉我們的迴圈。

第 16 行是 GCD 演算法的第 4(5) 點(第 3 點因為用 r 來算了,所以主要動機應該是第 4(5) 點),將 b

2 的因數都提出,假如可以被
2
整除(第 17 行),則回傳 r & -r ,這是因為 b 只有
2
這個因數,也就不需要進入輾轉相除的階段。

再來進入迴圈準備輾轉相除,第 22 行跟第 16 行一樣,假如把

2 的因數提光之後互質,則回傳 r & -r ,同時也是輾轉相除法停止的條件之一。下一行 a == b 也是停止條件之一,即商數等於餘數,回傳商數乘以
2
的公因數。

#else /* If normalization is done by loops, the even/odd algorithm is a win. */ unsigned long gcd(unsigned long a, unsigned long b) { unsigned long r = a | b; if (!a || !b) return r; /* Isolate lsbit of r */ r &= -r; while (!(b & r)) b >>= 1; if (b == r) return r; for (;;) { while (!(a & r)) a >>= 1; if (a == r) return r; if (a == b) return a; if (a < b) swap(a, b); a -= b; a >>= 1; if (a & r) a += b; a >>= 1; } } #endif

在沒有使用 ffs 的版本裡,可以看到使用迴圈取代 ffs 。最初第 38, 39 行一樣是 GCD 演算法第 1, 2 點。然後將 r 取只剩 LSB 一個位元,意思是 ra, b 兩者最大的

2 的公因數,即 GCD 演算法第 3 點。

從註解中可以看到 normalization ,在這裡的 norm 是

p-adic norm (
p
-adic absolute value)
,特別指
2
-adic norm 。簡單說
2
-adic order 就是一整數
n
所擁有
2
的因數最高次方數,而
n
2
-adic norm 則是
2
2
-adic order 次方之後取倒數,然後整數的長度和兩整數之間的距離可以由
2
-adic norm 來定義。

回到註釋, normalization 就是以 r 為基準,將 a, b

2-adic norm 變成跟 r 一樣,所以就是將 a, b 右移到三者的 LBS 位置一樣,此時的 a, b, r 三者的
2
-adic norm 相同,最大的
2
的公因數也相同。

Let

p be a prime number.

  1. The
    p
    -adic order for
    Z
    is the function
    νp:ZN
    defined by
    νp(n)={max{kN:pk|n}if n0if n=0
  2. The
    p
    -adic absolute value on
    Q
    is the function
    ||p:QR0
    defined by
    |r|p=pνp(r)

For example,

n=45 gives
ν3(n)=2
and
|n|3=3ν3(n)=19
.

(待續)

測驗 4

以下是 bit array (bitset) 的應用。







%0



expand

0

1

1

0

1

0

0

...

0

1

1

0

0

1

0

1

...

1

...
expand the bitmap



i
i



bitset1

1

1

0

0

1

0

1

...

1



i->bitset1:e1





bitmap_ent

bitmap[0]

bitmap[1]

...

bitmap[bitmapsize - 1]



bitmap_ent:e1->bitset1:e0





bitset0

0

1

1

0

1

0

0

...

0



bitmap_ent:e0->bitset0:e0





bitmap

*bitmap



bitmap->bitmap_ent:e0





bitmap 是一個 64 位元無號整數陣列,而每個整數都代表著一個 bitset ,使用 p + i 來計算 bit 的 index ,是將 bitmap 所有 bits 都攤平,所以 p 會是每加一次就是加 64 ,代表一個整數有 64 位元,然後加上在 bitset 裡的位移 i ,傳到 out 裡面。

#include <stddef.h>
size_t naive(uint64_t *bitmap, size_t bitmapsize, uint32_t *out)
{
    size_t pos = 0;
    for (size_t k = 0; k < bitmapsize; ++k) {
        uint64_t bitset = bitmap[k];
        size_t p = k * 64;
        for (int i = 0; i < 64; i++) {
            if ((bitset >> i) & 0x1)
                out[pos++] = p + i;
        }
    }
    return pos;
}

以下使用 GNU extension 的 __builtin_ctzll 作等價改寫。先觀察程式碼,改變的地方從 p = k * 64 開始,變成 while (bitset != 0) ,然後計算 trailing zero 的數量,將之與 k * 64 取和放進 out ,應該就是對應上面程式碼的 index 計算,然後做 XOR 運算,繼續下一個迴圈。所以推測應該每次迴圈就是計算一個 index ,放進 out 之後就要將那個 bit 給反轉,所以每次 r 都可以更新到左邊一個 bit 的位置,進而到 bitset 變為 0 ,算完所有的 bits ,結束迴圈。

#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 = EXP6;
            int r = __builtin_ctzll(bitset);
            out[pos++] = k * 64 + r;
            bitset ^= t;                           
        }
    }
    return pos;
}

這樣 EXP6 就要填那個準備反轉的 bit ,利用 bitset & -bitset 可以得到 LSB 位置,即該準備反轉的 bit 。

延伸問題:

  1. 解釋上述程式運作原理,並舉出這樣的程式碼用在哪些真實案例中;
  2. 設計實驗,檢驗 ctz/clz 改寫的程式碼相較原本的實作有多少改進?應考慮到不同的 bitmap density;
  3. 思考進一步的改進空間;
  4. 閱讀 Data Structures in the Linux Kernel 並舉出 Linux 核心使用 bitmap 的案例;

新算法可以跳過 0-bit ,只對 1 作計算。根據不同的 bit density ,密度越高對新算法越不利,因每次迴圈的計算量較多。

(檢驗待測)

測驗 5

LeetCode - Fraction to Recurring Decimal ,給一分數的分子與分母,回傳其以字串表示的小數。若有循環小數,則將循環部份用括號括起來。

#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; }

這裡對循環小數的對策是使用 hash table ,將遇到的餘數儲存,每算個新的位數,就與之前儲存的餘數作比對,檢查是否出現計算重複的餘數,即出現循環小數。

剛開始先檢查有沒有 trivial solution ,是否有分子為

0 或分母為
0
的情況。再來做正負號檢測,儲存在 sign 中。在第 35 行做了無謂的檢查,因為 n, p 在第 24, 26 行調整為正數,所以 division 一定是正數,不需再做檢查,即該行應改為 sprintf(p, "%ld", (long)division); 較佳。

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; }

看完前半部份再來到第 54 行的迴圈,第一次迴圈 hash table 還是空的,所以 pos 會是 -1 ,再來將儲存餘數進 hash table ,這邊注意到 keyint ,而 remainderlong long ,所以存進去的餘數有可能存不完整, find 在比對 key (truncated remainder) 與 node->key (truncated remainder) 是否相同時會有問題,應將 key 改為 long long 。之後將位數 i 存進 index ,插入 hash table ,所以 MMM 應該為 list_add ,而位置要跟 find 裡的 hash function 對應,所以 EEE 應該為 &heads[remainder % size] 。最後將該位數的商算出,放進 decimal ,然後更新 remainder ,進入下一個迴圈。

回到剛跳過的第 78 行區塊,在若 find 發現該餘數已經出現過了,則回傳前一次出現的位數 pos ,這時我們知道有循環小數,準備將剛才放在 decimal 的小數外圍用括號括起來,放到我們要回傳的 result 。第 80 行將 decimal 前面部份複製到 result 上,然後複製完後給一個括號,所以推測停止條件裡的 PPP 為 pos-- 。之後將循環部份複製,最後補個下括號和 null terminator 就可以回傳了。

若沒有遇到循環小數,也就是在整除的時候,結束迴圈,將 decimal 複製到 result 回傳結束。

延伸問題:

  1. 解釋上述程式碼運作原理,指出其中不足,並予以改進

    例如,判斷負號只要寫作 bool isNegative = numerator < 0 ^ denominator < 0;
    搭配研讀 The simple math behind decimal-binary conversion algorithms

  2. 在 Linux 核心原始程式碼的 mm/ 目錄 (即記憶體管理) 中,找到特定整數比值 (即 fraction) 和 bitwise 操作的程式碼案例,予以探討和解說其應用場景

測驗 6

由於要從記憶體上快速找到某資料位置,一種加速方式是對資料做對齊,讓現代的 CPU 讀寫時更有效率。而使用 __alignof__ 可以得知某資料的對齊大小,進而讓我們對資料位置的控制更方便。

以下 tint 代入為例,中間的 anonymous structure 第一個成員為 cchar 型態,第二個成員為 _hint 型態,有 4 bytes 。由對齊的特性, char 的地址會在 1 的倍數上,換句話說,任意地址皆可,反而 int 是 4 bytes ,所以地址會在 4 的倍數上。將這個 anonymous structure 放在 0 上, 0 為任意數的倍數,用第一個成員 c 佔據,然後以 char 這個最小單位給 _h 做位移,強制我們要計算對齊的對象 _h 的地址在 1 以後(含)。

這時有趣的事發生了,假如 t 的對齊是 1 的倍數,則此時 _h 就會剛好靠在 char 的旁邊,也就是地址 1 上,但現在我們 t 的對齊為 4 的倍數,所以 _h 不會在地址 1 上,而是出現在下一個 4 的倍數地址上,即地址 4 ,而中間空的部份稱為 padding 。







%0



struct
struct { char c; int _h; }



a

char

____

____

____

int 0

int 1

int 2

int 3



addr0
address = 0



addr0->a:c





p0
padding



p0->a:p0





p0->a:p1





p0->a:p2





addr2
address = 4



addr2->a:i





根據以上的例子,對齊可以由 _h0 的位移量得知,所以 M 應該指 _h 這個成員,得到前者 _h 的地址,而後者 X 對應 anonymous structure 的起始位置 0 當基準,計算位移量。

/*
 * 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)->M) - (char *)X)

這邊有兩個問題,為何要轉換成 (char *) ?還有為何要相減,直接拿 _h 的地址不就好了。第一個問題,因為要做指標的減法,所以需要兩個指標的型態相同。而對於前者的對齊可能很多種,不如直接選 (char *) ,因為任何數都是 1 的倍數,這樣可以避免某些錯誤和編譯器的警告。

第二個問題目前只是猜測,根據 N1256 (6.5.6-9) ,兩者相減得到的型態不會是 (char *) ,而是 ptrdiff_t ,也就是會有轉換型態發生。用數學舉例,在二維空間中,兩點的座標相減得到是向量,而向量跟座標是不一樣的概念。

N1256 (6.5.6-7)
For the purposes of these operators, a pointer to an object that is not an element of an array behaves the same as a pointer to the first element of an array of length one with the type of the object as its element type.

N1256 (6.5.6-9)
When two pointers are subtracted, both shall point to elements of the same array object, or one past the last element of the array object; the result is the difference of the subscripts of the two array elements. The size of the result is implementation-defined, and its type (a signed integer type) is ptrdiff_t defined in the <stddef.h> header. []

N1256 (6.2.5-27)
A pointer to void shall have the same representation and alignment requirements as a pointer to a character type. Similarly, pointers to qualified or unqualified versions of compatible types shall have the same representation and alignment requirements. All pointers to structure types shall have the same representation and alignment requirements as each other . All pointers to union types shall have the same representation and alignment requirements as each other. Pointers to other types need not have the same representation or alignment requirements.

延伸問題:

  1. 解釋上述程式碼運作原理
  2. 在 Linux 核心原始程式碼中找出 __alignof__ 的使用案例 2 則,並針對其場景進行解說
  3. 在 Linux 核心源使程式碼找出 ALIGN, ALIGN_DOWN, ALIGN_UP 等巨集,探討其實作機制和用途,並舉例探討 (可和上述第二點的案例重複)。思索能否對 Linux 核心提交貢獻,儘量共用相同的巨集

測驗 7

FizzBuzz 是一個數學遊戲,剛開始由

1 開始往上數,每當數字可以被
3
整除,這時玩家不該喊數字,而是代替那數字喊 Fizz 。同樣地,數字可以被
5
整除,則是喊 Buzz 。當可以被
3
,
5
同時整除,則是喊 FizzBuzz 兩個同時。

/* naive.c */
#include <stdio.h>
int main() {
    for (unsigned int i = 1; i < 100; i++) {
        if (i % 3 == 0) printf("Fizz");
        if (i % 5 == 0) printf("Buzz");
        //if (i % 15 == 0) printf("FizzBuzz"); /* no need to print it again */
        if ((i % 3) && (i % 5)) printf("%u", i);
        printf("\n");
    }
    return 0;
}

上面為直覺的作法,當數字 i 可以被

3 整除就印出 Fizz ,當數字可以被
5
整除則印出 Buzz ,當可以同時被兩者整除時,則在兩個 if 檢查的時候就會依序印出。最後都不能整除的時候,則是印出該數字 i

string literal: "FizzBuzz%u"
        offset:  0   4   8

觀察上面同時整除的時候,依序印出就像是兩個字接在一起,換個方式想,那可以控制字串的起始位置和長度,來選擇印出的字串為何,就像上面的 "FizzBuzz%u" 那樣,選擇 offset 和字串長度來控制要印出的格式。

/* bitwise.c */
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 << KK1) << KK2;

        char fmt[9];
        strncpy(fmt, &"FizzBuzz%u"[(9 >> div5) >> (KK3)], length);
        fmt[length] = '\0';

        printf(fmt, i);
        printf("\n");
    }
    return 0;
}

先跳過 is_divisible 函式,到後面再來分析,從名字可以看出就是判斷 n 是否能被 bitmask M 對應的整數整除。那開始看主函式,使用我們剛才的想法,需要知道要印出的字串起始位置和長度,有的資訊則是可否被

3,
5
整除,下面用表格表示相對關係。

dev3 dev5 start index length
0 0 8 2
0 1 4 4
1 0 0 4
1 1 0 8

從此表可以看出 length 是每多一個可整除的因數(

3 or
5
)就多增加一倍,所以可以直接 KK1 填 dev3 , KK2 填 dev5 (可以互換)。

而起始位置原為 &"FizzBuzz%u"[(9 >> div5) >> (KK3)] 裡的 9 應改為 8 才對,不然不符合 KK3 的填空要求。從表格觀察,若 div50 ,則起始位置為 8 ,沒有右移,若 div51 ,則起始位置為 4 ,這兩個事件在 (8 >> div5) 就已經寫好了。那麼再來是 dev3 ,若為 0 則保持原狀,若為 1 則歸零,所以 KK3 應該為右移位數為零或者右移位數大於

4 的一個表示式(可以將 8 歸零),所以可以填 dev3 << 2

延伸問題:

  1. 解釋上述程式運作原理並評估 naive.cbitwise.c 效能落差
    • 避免 stream I/O 帶來的影響,可將 printf 更換為 sprintf
  2. 分析 Faster remainders when the divisor is a constant: beating compilers and libdivide 一文的想法 (可參照同一篇網誌下方的評論),並設計另一種 bitmask,如「可被 3 整除則末位設為 1」「可被 5 整除則倒數第二位設定為 1」,然後再改寫 bitwise.c 程式碼,試圖運用更少的指令來實作出 branchless;
  3. 研讀 The Fastest FizzBuzz Implementation 及其 Hacker News 討論,提出 throughput (吞吐量) 更高的 Fizzbuzz 實作
  4. 解析 Linux 核心原始程式碼 kernel/time/timekeeping.c 裡頭涉及到除法運算的機制,探討其快速除法的實作 (注意: 你可能要對照研讀 kernel/time/ 目錄的標頭檔和程式碼)

    過程中,你可能會發現可貢獻到 Linux 核心的空間,請充分討論