Try   HackMD

2022q1 Homework2 (quiz2)

contributed by < cwl0429 >

測驗 1

觀察下列程式碼

  • a >> 1b >> 1 的功用都是將其除 2,此時會遇到最低位元遺失問題
  • 需要將最低位元補上,方法是 a 和 b 作 AND,接著和 1 作 AND,便可取得遺失的位元,因此 EXP1a & b & 1
#include <stdint.h>
uint32_t average(uint32_t a, uint32_t b)
{
    return (a >> 1) + (b >> 1) + (EXP1);
}

再觀察改寫後的程式碼,發現完全想不到該如何下手,於是嘗試用湊數字的方式,最後得出:

  • EXP2a & b
  • EXP3a ^ b
uint32_t average(uint32_t a, uint32_t b)
{
    return (EXP2) + ((EXP3) >> 1);
}

後來,參照了 laneser 的筆記,得知其實 x + y = (x ^ y) + ((x & y) << 1),而 (x + y) / 2 就等於 ((x ^ y) + ((x & y) << 1)) >> 1 也就是 (x & y) + ((x ^ y) >> 1)
因此得出:

  • EXP2a & b
  • EXP3a ^ b

比較上述在編譯器最佳化開啟的狀況

先放一張來自 x64 文件的圖片,方便檢視暫存器的大小

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

接著嘗試理解下列四種 average 的組合語言,使用 gcc -S 1-1.c -O3 -fno-asynchronous-unwind-tables 生成

  • -O3 開啟編譯器最佳化第四層
  • -fno-asynchronous-unwind-tables 過濾不必要的 .cfi
uint32_t average1(uint32_t a, uint32_t b)
{
    return (a + b) / 2;
}

average1:
	endbr64
	leal	(%rdi,%rsi), %eax
	shrl	%eax
	ret

根據 Introduction to the leal instruction 的描述

Load effective address: leal S,D # D<&S, where D must be a register, and S is a Memory operand.

leal looks like a mov instr, but does not access Memory. Instead, it takes advantage of the addressing circuitry and uses it to do arithmetic (as opposed to generating multiple arithmetic instructions to do arithmetic).

(ex) if edx holds the value of x:
leal (%eax),%ecx # R[%ecx]<&(M[R[%eax]])

# this moves the value stored in %eax to %ecx

The key is that the address of (M[ at address x ]) is x, so this is moving the value stored in %eax to %ecx; there is no memory access in this instruction's execution.

  • leal 會將 (%rdi,%rsi) 的值放至 %eax,其間沒有 memory access 發生
  • 接著利用 shrl%eax 右移一位元,達成除 2 的效果
uint32_t average2(uint32_t low, uint32_t high)
{
    return low + (high - low) / 2;
}

average2:
	endbr64
	movl	%esi, %eax
	subl	%edi, %eax
	shrl	%eax
	addl	%edi, %eax
	ret
  • movlhigh 存入 %eax
  • subl 使 %eaxlow 後,將 %eax 往右移一個 bit
  • addl 將結果相加並存放於 %eax
uint32_t average3(uint32_t a, uint32_t b)
{
    return (a >> 1) + (b >> 1) + (a & b & 1);
}

average3:
	endbr64
	movl	%edi, %eax
	movl	%esi, %edx
	andl	%esi, %edi
	shrl	%eax
	shrl	%edx
	andl	$1, %edi
	addl	%edx, %eax
	addl	%edi, %eax
	ret
  • xymovl 讀至 %eax%edx
  • 接著做對應的位元運算 andlshrl
  • 最後便將三個結果 %edx, %edi%eax 相加,存放在 %eax
uint32_t average4(uint32_t a, uint32_t b)
{
    return (a & b) + ((a ^ b) >> 1);
}

average4:
	endbr64
	movl	%edi, %eax
	andl	%esi, %edi
	xorl	%esi, %eax
	shrl	%eax
	addl	%edi, %eax
	ret
  • a 存入 %eax ,此時 %edi%eax 為 a %esi 為 b
  • 做對應的位元運算,andlxorl
  • 完成後將 %eax 右移
  • 最後將 a & b 結果加至 eax

研讀 Linux 核心原始程式碼 include/linux/average.h,探討其 Exponentially weighted moving average (EWMA) 實作,以及在 Linux 核心的應用

TODO..

測驗 2

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

觀察題目,若 a 較大,則需 return a ^ 0 才能成立,反之則需要 return a ^ (a ^ b) 才能回傳 b

到此根據老師的提示,便可猜出 EXP4 應為 a ^ b

剩下 EXP5 要決定回傳 a ^ 0 還是 a ^ (a ^ b)

  • a ^ 0
    • 代表 a ^ ((a ^ b) & 0x0)
  • a ^ (a ^ b)
    • 代表 a ^ ((a ^ b) & 0xFFFFFFFF)

所以得出 EXP5a < b

針對 32 位元無號/有號整數,撰寫同樣 branchless 的實作

無號數版本能和有號數幾乎一樣,更改回傳型態及參數型態即可

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

Linux 核心也有若干 branchless / branch-free 的實作,請在 Linux 核心原始程式碼中,找出更多類似的實作手法。請善用 git log 檢索。

drivers/gpu/drm/radeon/r600_blit.c

Replace int2float() with an optimized version. We use
__fls() to find the most significant bit. Using that, the
loop can be avoided. A second trick is to use the
behaviour of the rotate instructions to expand
the range of the unsigned int to float
conversion to the full 32 bits in a branchless way.

The routine is now exact up to 2^24. Above that, we truncate which is equivalent to rounding towards zero.

Signed-off-by: Steven Fuerst svfuerst@gmail.com

這份 commit 使用 branchless 的技巧優化了 int2float

主要優化的位置在第 12 行

uint32_t int2float(uint32_t input) { u32 result, i, exponent, fraction; if ((input & 0x3fff) == 0) result = 0; /* 0 is a special case */ else { exponent = 140; /* exponent biased by 127; */ fraction = (input & 0x3fff) << 10; /* cheat and only handle numbers below 2^^15 */ for (i = 0; i < 14; i++) { if (fraction & 0x800000) break; else { fraction = fraction << 1; /* keep shifting left until top bit = 1 */ exponent = exponent - 1; } } result = exponent << 23 | (fraction & 0x7fffff); /* mask off top bit; assumed 1 */ } return result; }

新的 int2float 利用 __fls 找到 msb 的位置,再利用 ror rotate 到正確位置,使我們不需判斷是否低於 2 ^ 24

uint32_t int2float(uint32_t x)
{
	uint32_t msb, exponent, fraction;

	/* Zero is special */
	if (!x) return 0;

	/* Get location of the most significant bit */
	msb = __fls(x);
	/*
	 * Use a rotate instead of a shift because that works both leftwards
	 * and rightwards due to the mod(32) behaviour.  This means we don't
	 * need to check to see if we are above 2^24 or not.
	 */
	fraction = ror32(x, (msb - I2F_FRAC_BITS) & 0x1f) & I2F_MASK;
	exponent = (127 + msb) << I2F_FRAC_BITS;

	return fraction + exponent;
}

測驗 3

GCD 演算法

  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
    brcause everything divides 0.
  3. If x and y are both even,
    gcd(x,y)=2gcd(x/2,y/2)
    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(x/2,y)
    .
  5. On similar lines, if x is odd and y is even, then
    gcd(x,y)=gcd(x,y/2)
    . It is because 2 is not a common divisor.
if (!u || !v) return u | v;

參考老師提供的 GCD 演算法,這段程式碼實作了 1. 及 2. ,也就是 gcd(0, 0), gcd(x, 0) 及 gcd(0, y) 的情況

for (shift = 0; !((u | v) & 1); shift++) {
        u /= 2, v /= 2;
    }

這是 3. 的情況,當二數的最低位元皆為 0 時,便將二數都除 2 並且紀錄除的次數,也就是 shift 的數值

while (!(u & 1))
        u /= 2;

這是 4. 的情況,當 x 為偶數時,gcd(x ,y) = gcd(x / 2, y)

    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 (COND);
return RET;

最後這段則是 5. 的情況及 GCD 的原始算法

這邊需要找出 CONDRET 為何,觀察 COND ,得知其作用為中止 while 的判斷式,那根據 GCD 的原始定義,當兩數相減至其中一數為 0 時則中止,而且得知

  • CONDv
    最後 RET 需要回傳二數中剩餘非 0 的數,也就是 u,但要注意在前面的 shift 代表了二數整除於 2 多少次數,所以需要將 u 左移 shift
  • RETu << shift

在 x86_64 上透過 __builtin_ctz 改寫 GCD,分析對效能的提升

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. 取自 GNU 官方網站

__buildin_ctz 能取得 unsigned int x 從最低位元算起的連續 0-bits 個數,再回頭看 gcd_64 程式碼,得知主要想改善的部份為:

while (!(u & 1))
        u /= 2;

於是將相關操作改成 __buildin_ctz 後,得到以下實作:

uint64_t gcd64_ctz(uint64_t u, uint64_t v)
{
    if (!u || !v) return u | v;

    int shift;

    int u_ctz =__builtin_ctz(u);
    int v_ctz =__builtin_ctz(v);

    shift = u_ctz > v_ctz ? v_ctz : u_ctz;  
    u >> u_ctz;
    v >> v_ctz;
    
    do {
        v = v >> __builtin_ctz(v);
        if (u < v) {
            v -= u;
        } else {
            uint64_t t = u - v;
            u = v;
            v = t;
        }
    } while (v);
    return u << shift;
}

接著利用 perf 進行效能測試

用亂數進行測試, N 的大小為 10000

int main()
{
    srand(time(NULL));

    for(int i = 0; i < N; i++) {
        //gcd64(rand(), rand());
        gcd64_ctz(rand(), rand());
    }

    return 0;
}

使用以下命令進行測試
perf stat --repeat 10 ./1-3

gdb64

0.0055975 +- 0.0000618 seconds time elapsed  ( +-  1.10% )

gdb64_ctz

0.002960 +- 0.000433 seconds time elapsed  ( +- 14.63% )

可看出 ctz 的版本約快了 1.8910 倍

Linux 核心中也內建 GCD (而且還不只一種實作),例如 lib/math/gcd.c,請解釋其實作手法和探討在核心內的應用場景。

研究 ffs 的用途

ffs — find first set bit in word

得知 ffs 能找到從右邊數起的第一個 1 位置,再根據維基百科的解釋, ffs 和之前題目用到的 ctz 的關係能用數學式子表示

在 x 不等於 0 的情況下

ctz(x)=ffs(x)1
舉例來說, x = 0x11000000ctz(x) = 6ffs(x) = 7

但實際在 trace 此函式時,若是根據上述的定義,這個函式有機會進入無窮迴圈,才發現 __ffs 並不等於 ffs

於是重新翻了資料,找到 __ffs 實際定義的地方 asm-generic/bitops/__ffs.h

/**
 * __ffs - find first bit in word.
 * @word: The word to search
 *
 * Undefined if no bit exists, so code should check against 0 first.
 */

因此 __ffsctz 會回傳一樣的結果

另一個需要注意的地方是 r & -r 會保留從右往左數第一個為 1 的位元,其餘位元皆設定為 0

有趣的是,若將 r & -r

log2 ,其實就會是 ctz(r) 的結果
也學到教訓,找資料請盡量找第一手資料

解決上述兩點後,能觀察到這個版本的 GCD 和我們的版本差異為使用了 __ffs 進行加速,演算法相同

unsigned long gcd(unsigned long a, unsigned long b)
{
	unsigned long r = a | b;
        
	if (!a || !b) /*1. If both x and y are 0, gcd is zero $gcd(0, 0) = 0$.*/
		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;
	}
}

接著來探討 lib/math/gcd.c 實際的應用情境

利用 grep -r "#include <linux/gcd.h>" 列出有哪些檔案使用了 lib/math/gcd.c ,發現在 sound, net, drivers 等等都有用到 gcd,這邊我嘗試理解 linux/sound/soc/codecs/arizona.c 內 gcd 的用途

觀察目錄,發現 arizona.c 和音訊解碼有關