contributed by < vacantron
>
考慮以下對二個無號整數取平均值的程式碼:
當兩數都接近 32 位元無號整數的最大值時,若貿然的將兩數相加可能會導致溢位,需要採取其他的做法。若我們知道兩個數大小關係,則可以改寫為以下所示,其數值並不會大於較大的那個數,而免除了溢位的風險:
而使用 bitwise 操作不需事先得知兩數的關係:
其運作原理是將分別將兩數除以二並相加,但因為使用右移運算子會失去最低位元的資訊,故使用 a & b & 1
處理最低位元相加會進位的狀況
而我們可以再次改寫為:
運用除法的分配律,假設 a = 0x1110, b = 0x0111 ,我們可以分項成 a = 8 + 4 + 2, b = 4 + 2 + 1 ,用 & 運算子提出兩數共同的 4, 2 (除以二後會留下一個),再加上兩數各自剩餘部份的一半 (使用右移與 | 運算子)
未開啟最佳化:
開啟最佳化 (-O1):
開啟編譯器最佳化後,可以發現到為 average()
對函式所配置的記憶體堆疊操作被帶有加數與被加數的暫存器操作所取代,省去了記憶體的讀取、寫入,從而提昇程式的效率
Linux kernel 原始碼中的 average.h 實作了 exponentially weighted moving average (EWMA) , EWMA 是給予先前與目前資料不同的權重來計算平均的一種統計方法,其主要組成如下,其中 為權重、 為新舊數值:
權重 為一個介於 0, 1 間的小數,考量到電腦處理浮點數的開銷較大,我們可以把上述式子分子分母同乘 改寫為:
在註解中有特別提到:Note that this parameter must be a power of two for efficiency.
讓我們可以使用位移來取代除法運算,也使用了 BUILD_BUG_ON_NOT_POWER_OF_2
這個
巨集在編譯時期檢查是否為二的冪
static inline void ewma_##name##_add
函式中也使用了 READ_ONCE
與 WRITE_ONCE
來避免資料不一致的問題 (這兩個巨集也可以避免編譯器最佳化造成並行時沒有確實從記憶體中更新變數數值)
ewma 在 linux kernel 中的無線通訊部份重複出現了許多次,主要為計算無線訊號的強度 Received Signal Strength Indication (RSSI)
不清楚在 average.h #58 中將 _precision 指派到新的 變數中以及參數 _precision 的用意為何
以下試著推導數學式:
假設 precision 為 p 、weight 為 w ,成員 internal 在加入第一個數時使用 val << p
,將數值左移後儲存,其餘的使用 (((internal << w) - internal) + (val << p)) >> w
,這段實作可以改為 ,
在讀取 ewma 結構體的 internal 成員時的 internal >> (p)
會將分子的 項給除掉,precision 對這裡的計算看起來沒有任何的幫助,自行使用程式測試代入不同的值也得到相同的結果
20220420: 上面所說的 " 除以 " 並不完全正確,只有呼叫 read()
時才要進行除法 (位移) ,若只是單純的 add()
則是把上次的結果與新的輸入相加,因為在向右位移時會失去低位元的值,所以使用左移精度 p 來保護低位元的數值,所以平常的 internal
的狀態是一個精度為 p 的 定點數 。若沒有精度存在, add()
的行為會變成相加後取不大於結果的最大整數,沒有小數點後的數值就無法處理由小數進位到整數的狀況。之前使用程式測試時看起來沒有錯是因為只有加入兩三個數而已,項次還不夠多的狀況下還沒觸發小數的進位
改寫 〈解讀計算機編碼〉 一文的「不需要分支的設計」一節提供的程式碼 min
,我們得到以下的 max
實作:
試著解析這個函式:當 a > b 時應該要輸出 a,代表 EXP4 & -EXP5
必須為 0。而 a < b 時要輸出 b,我們可以利用 n ^ n = 0
這項性質推導出 EXP4 & -EXP5 = a ^ b
根據提示 EXP4 為 bitwise 操作、EXP5 為比較操作,得出 EXP4 = a ^ b 且 a < b 時 EXP5 的全部位元都是 1、 a > b 時為 0,得出 EXP5 = a ≦ b, 如果 a = b 時輸出 a 或 b 都是正確的,但求最精簡的表示式就省去等於了
嘗試對 32 位元有號整數 max
的實作:
由兩個有號數的符號關係作為 mask。若符號不同,會回傳正的那個數;若符號相同才會比較數值,因為符號相同的有號數去除符號後正好可以直接比較大小,例如:
20220421: 其實根本不用去掉符號,因為已經確定兩數的符號相同,直接讀取數值當成無號數處理就可以了
有號數數值 | 二補數 | 去除符號 (由大到小) |
---|---|---|
-1 | 1111 | 0111 |
-5 | 1011 | 0011 |
-8 | 1000 | 0000 |
7 | 0111 | 0111 |
5 | 0101 | 0101 |
2 | 0010 | 0010 |
若符號不同,-(~(SIGN(x) ^ SIGN(y))))
為 0,代表不去比較除了符號外的數字,由 x ^ (x ^ y) & -(SIGN(x) > SIGN(y))
用符號判斷大小關係;若符號相同,-(~(SIGN(x) ^ SIGN(y)))
是一個全部位元都是 1 的 mask,搭配 x ^ x = 0
的性質及其他運算求得較大的那個數
上面的實作有點冗長,我們換個方式思考:若要回傳 a 或 b 其中之一,那必定會有 a ^ (a ^ b) & ( X )
的結構,而 X 在 a > b 時要為 0、a < b 時是一個全部位元都是 1 的 mask,依照這個思考方向可以得到下面的改寫:
這個 GCD 實作主要是先去除有 2 的冪倍的因數,縮小數值後可以減少輾轉相除的次數。在 while 迴圈中使用 !(x & 1)
作為 while 的判斷條件確認 x 是否為 2 的倍數,去除因數中 2 的成份,因為 ,減小被除數可以減少輾轉相除的次數。最後將一開始移出的 2 的冪使用向左位移歸還回去
跳出 do - while 迴圈的 COND
要為 false (0),推測為輾轉相除法最後一步兩數相減後為零的情況,所以 COND
= v
。而 RET
要將一開始拿掉的 2 的冪倍放回來,因此 RET
= u << shift
在 x86_64 上透過 __builtin_ctz 改寫 GCD,分析對效能的提升:
將原本的 x /= 2
替換成 x >> 1
、不斷除以二直到該數為奇數的操作替換成 x >>= __builtin_ctzl(x)
,因為如果一個數是 的倍數,它的尾端就會有 個 0
將使用 __builtin_ctzl()
和原始的版本對由 rand()
產生的亂數執行一百萬次,透過 perf 可以得知兩者的相對效率,約提昇了 2.5 倍的速度:
在測試在 linux 核心內建的 lib/math/gcd.c #23 其中之一的實作時,會回傳的錯誤的 GCD 數值,另一種實作和從 git 版本記錄的最初版本的 GCD 函式都能回傳正確值,不知道是測試方式有誤還是實作上的問題
將要計算的數值以空白隔開輸入,可得到各個 GCD 實作的計算結果
20220420: 從 lib/math/gcd.c 挑出一些程式碼進行分析:
a | b
: 若用數學式描述 a、b,假設 其中 a > b > … > g,則可以提出 ,得出 ,b 也由相同的手法得到 其中 m > n > … > z 。若 g > z, 可以提出 z,則 a | b
= ,而後兩項多項式因為都含有 的部份所以皆為奇數,這個操作像前面提到的提出 2 的冪降低輾轉相除次數,並且 a | b
尾端有 z 個 0。
r & -r
: 這裡可以視為一個 mask,只留下最低位元的 1 。結合上面的 a | b
,只留下 ,實際上的意義是兩個數最大的 2 的冪的因數
從第 9 行的程式碼的行為反推:若 b 經過某種向右位移後等於 1,表示比某一元位元高位的位元都是 0,且回傳 r | -r
。 b 應該要是一個只有一個位元為 1 的數,可以表示成 2 的冪,因為沒有其他多項式的部份所以不用擔心產生其他非 2 的冪的因數。因此這個位移操作應該要讓唯一的 1 落到 LSB 上
for 迴圈中會進行前面提到的輾轉相除法
最後,輾轉相除法的中止條件為 ,或是 a == b,像前面提到的把一開始移去的 2 的冪倍向左位移回去,這裡應該要左移 z
綜合以上,我認為函式 __ffs()
應該要替換成 __builtin_ctz
或類似的 ctz 函式。經過測試確認替換函式後即可輸出正確的數值
< eecheng >: 以下是可測試的 GCD 程式碼,__ffs() 的表現看起來是對的
還好沒有導致在 perf 原理和實務 中提到的安裝 perf 時顯示錯誤,但還是試著自行重新編譯安裝一次
原本以為在對應的 kernel archive 中找到檔案並編譯完就可以了,但在執行時發現 symbols 會顯示無法直接閱讀的十六進位的 symbol 而不是函式名稱:
在網路上找到了相關的參考資料及解法,後來也注意到編譯時有顯示缺少一些相依套件而沒有啟用一些 feature (螢幕刷新太快一開始沒有注意到):
安裝了 libw-dev 後:
可以看到有關 dwarf 的 feature 被啟用了,根據官網的說明,它像是 GDB 一樣是在編譯時插入 debugger 程式碼時的中間層,像上面的 symbol 轉換就是其中之一的功能
後續把相依套件安裝完重新編譯後即可使用
bit array 由一連串的 0, 1 組成,實際上可以是開關或是表示有無,例如點陣圖的像素儲存,以下使用 __builtin_ctrl()
改寫程式碼,其中 *out
指向一個儲存 bit 為 1 位置的陣列:
用 64 位元的整數裡面的位元儲存 bitset,其中 t 作用是要將最低位元的 1 清除,可以視為一個 mask,因此答案應為 bitset & -bitset
。我們可以把 -bitset 拆解成 (bitset ^ -1) + 1
,將全部的位元反轉可以保留住前面的其他沒有要清除的位元,以及將尾端的 0 全部換成 1,最後由 + 1
會使尾端一連串的 1 進位,直到 bitset 最低位元的 1 的位置才會停止
改寫後的程式碼不用迭代每一個位元去檢查是否為 1,只需要使用 __builtin_ctzl()
就能知道後面有幾個 0
使用 perf 並且加入 probe 測量函式執行所花費的時間:
將編譯過後的可執行檔交給 perf 加工,讓函式在執行時觸發事件 (event),再由 perf 捕捉匯集成記錄
函式 | bitmap density | 花費時間 (μs) |
---|---|---|
naive | 0% | 1347 |
25% | 1573 | |
50% | 1782 | |
75% | 1629 | |
100% | 1395 | |
improved | 0% | 28 |
25% | 84 | |
50% | 140 | |
75% | 207 | |
100% | 274 |
可以發現當 bitmap density 越小時 ( 1 越稀疏時),修改過後 (使用 __builtin_ctzl()
) 的程式明顯比其他還有效率,隨著 bitmap 的 1 增加,__builtin_ctzl()
的影響力相對的就下降了。原始版本的程式的時間分佈為 0, 1 各半時最長,0 或 1 的密度較高時花費時間較短。推測與 branch-misses 的次數有關