2022q1 Homework2 (quiz2)

contributed by < tinyynoob >

作業要求

測驗 1

解釋程式碼

這題的目標是對兩數取平均，由於直覺的寫法會有 overflow 的問題，因此這裡試圖用 bitwise operation 來達成。

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

這裡的 >> 1 顯然是除以 2 的意思，那麼這個 EXP1 應該就是來解決兩個奇數分開運算造成的錯誤，例如：

(1 >> 1) + (1 >> 1) = 0 neq 1

要確認兩個都是奇數，我想到的最簡易方法是

EXP1 = a & b & 1

取平均的另外一種等價的寫法

uint32_t average(uint32_t a, uint32_t b)
{
    return (EXP2) + ((EXP3) >> 1);
}

根據你所不知道的 C 語言：數值系統篇裡的提點，我們了解到可以直接在 C 語言上應用數位邏輯設計，嘗試再次推導：

(x + y) / 2
=> (x + y) >> 1
=> ((x ^ y) + cin) >> 1
=> ((x ^ y) + ((x & y) << 1)) >> 1
=> ((x ^ y) >> 1) + (x & y)

C Operator Precedence

注意到不要把這個數學推導過程的左右移想成 C 語言內的左右移，否則會誤以為最後一個步驟不相等。

這裡的加法長得是如此簡單又沒有 ripple 的問題，使我納悶在邏輯設計時為何 carry lookahead adder 要用那麼複雜的方式，或許 shift register 影響此法的硬體可行性？

已解決。
前面推導的過程將 + 誤改成 ^ 導致我想錯了

結論得

EXP2 = a & b
EXP3 = a ^ b

編譯器最佳化下的比較

開啟編譯器最佳化 -O2 ，對於上述兩個函式得到以下結果：






















...
average1:
	endbr64
	movl	%edi, %eax
	movl	%esi, %edx
	andl	%esi, %edi
	shrl	%eax
	shrl	%edx
	andl	$1, %edi
	addl	%edx, %eax
	addl	%edi, %eax
	ret
...
average2:
	endbr64
	movl	%edi, %eax
	andl	%esi, %edi
	xorl	%esi, %eax
	shrl	%eax
	addl	%edi, %eax
	ret
...

兩種 C 語言寫法在 assembly 層級會相差 3 行指令，主要的差異來自於第 5, 8, 10 行。

第一種算法 average1 意圖在 %eax 和 %edx 分別算出 a >> 1 和 b >> 1 ，並在 %edi 算出 a & b & 1 ，最後加總到返回值 %eax 。
可以想見若我們先做 a & b & 1 則可以不必用到 4 個 register ，也有機會減少 1 行指令，變成：

movl	%edi, %eax
andl    %esi, %eax
andl    $1, %eax
shrl    %edi
shrl    %esi
addl    %edi, %eax
addl    %esi, %eax
ret

不過我認為編譯器原先的作法可以降低前後指令的相依關係，帶來更好的平行性。假如以指令的行數來表示，則可以有：

\begin{aligned} 4 \to 7 & \to 10 \\ 5 \to 8 & \to 11 \\ 6 \to 9 & \underset{同 步}{\underset{⏟}{\to}} \end{aligned}

在多處理器的情況下有機會做 3 步就得出結果

第二種算法 average2 雖然生成的組合語言指令較少，但前後的相依性大，在多處理器的情況下可能需要執行更久：

\begin{aligned} 16 \to 18 \to 19 & \to 20 \\ 17 & \underset{同 步}{\underset{⏟}{\to}} \end{aligned}

不過如此短的計算應該還是會在單一處理器上完成，結論而言我認為仍是法二較佳。

不確定我這裡對 multicore 相關的認知是否正確

研讀 kernel 中的 EWMA

linux/average.h

從定義

#define DECLARE_EWMA(name, _precision, _weight_rcp)

可以看出這個 macro 接收了 3 個參數，並且可以從註解中得知

precision 是定點數在小數點以下的精度
weight_rcp 是加權平均的權重，令它為
$w$ 則

$\begin{matrix} (1) & y_{n} = \frac{1}{w} \cdot x_{n} + \frac{w - 1}{w} \cdot y_{n - 1} \end{matrix}$
其中
$x_{n}$ 為原始輸入，
$y_{n}$ 是移動加權後的值

w

被規定須選為 power of 2 ，如此一來除法便可以化為簡單的右移運算。

我們主要需要解釋的是這一段程式碼：

unsigned long weight_rcp = ilog2(_weight_rcp);	
...
WRITE_ONCE(e->internal, internal ?			\
	(((internal << weight_rcp) - internal) +	\
		(val << precision)) >> weight_rcp :	\
	(val << precision));				\

令 weight_rcp =

\log_{2}

_weight_rcp

若 internal 非 0 ，則依據 Equation

(1)

更新 internal 值。若為 0 ，則可以直接省略加號右邊的部份。

WRITE_ONCE() 的功能是在多執行緒的情境下進行正確寫入，把第二個參數的值 assign 給第一個參數。

Linux内核中的READ_ONCE和WRITE_ONCE宏

不知道哪裡可以找到 WRITE_ONCE() 的第一手資料，它不在 compiler.h 裡

我推測這裡是使用定點數運算，不過 >> precision 不是等同於把低位元丟棄嗎，這是如何把 unsigned long 變成定點數的呢？

moving average 的功能主要是讓資料變得更平滑，不因突發的波動而隨之起舞，其實它就是離散時間的低通濾波器。以這題而言，我們可以嘗試由 Z transform 計算 frequency response ：

\begin{aligned} Z {y_{n}} & = Z {\frac{1}{w} \cdot x_{n} + \frac{w - 1}{w} \cdot y_{n - 1}} \\ Y (z) & = \frac{1}{w} X (z) + \frac{w - 1}{w} \cdot z^{- 1} Y (z) \\ H (z) = \frac{Y (z)}{X (z)} & = \frac{z}{1 - w + z w} \end{aligned}

換成頻域得

H (e^{j Ω}) = \frac{e^{j Ω}}{1 - w + w e^{j Ω}} = \frac{1}{1 - (1 - \frac{1}{w}) e^{- j Ω}}

對於頻率響應形式為

\frac{1}{1 - a e^{- j Ω}}

的離散濾波器，有：

filter type {\begin{cases} highpass, & - 1 < a < 0 \\ lowpass, & 0 < a < 1 \\ unstable, & | a | \geq 1 \end{cases}

因為

0 < (1 - \frac{1}{w}) < 1

，我們的 moving average 是一個低通濾波器。

探討 kernel/sched/fair.c 中的 Low-pass filter

studying…

測驗 3

解釋程式碼























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

該程式計算 gcd ，首先處理 u, v 有 0 的情形。

第 6, 7 行基於

\gcd (x, y) = 2 \cdot \gcd (\frac{x}{2}, \frac{y}{2})

，先將答案中 2 的成份 (以質因數分解的角度) 之 power 存到 shift 。

此時必有其一為奇數，接著要利用

\gcd (x, y) = \gcd (\frac{x}{2}, y)

的性質，先嘗試將 u 的大小盡量降低。

最後 11 到 21 行進行 Euclidean algorithm ，不停做相減直到出現 1 ：

由於 v 是每次相減後新出現的值，有可能又是偶數， 12, 13 行先確認能不能再用性質把值降下來。
將大減小放到 v 並將原先較小的數放到 u 。

因此 COND 為 Euclidean algorithm 繼續的條件，即是 v 不等於 1 ，考量在計算過程中 u, v 可能變相同，例如：

\gcd (3, 3) = 3

，因此 v 的中止條件也可能為 0 ，總結來說 COND 為 (v != 1) || (v != 0) 。

答案便是留下來的 u 值再補回早先存的 shift ，因此 RET 為 u << shift 。

拿去測試之後發現有錯，練習使用 gdb 除蟲：

以 u = 3, v = 5 呼叫
break 在 while (上方程式的 21 行)，依序得

`u`	`v`
3	2
1	2
1	0

接著程式卻沒有正確結束

判斷 while 條件有問題，將 COND 更正為 (v != 1) && (v != 0) 後得到理想結果。

透過 `__builtin_ctz` 改寫並比較效能

由於目前程式碼中存在大量判斷偶數與右移的操作，因此可以用 ctz() 來處理，從 gcc doc 找到說明：

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.

improved 版實作程式碼：

#define min(x, y)          \
    ({                     \
        typeof(x) _x = x;  \
        typeof(y) _y = y;  \
        _x < _y ? _x : _y; \
    })

static uint64_t gcd64_imp(uint64_t u, uint64_t v)
{
    if (!u || !v)
        return u | v;
    int shift = min(__builtin_ctz(u), __builtin_ctz(v));

    u >>= __builtin_ctz(u);
    do {
        v >>= __builtin_ctz(v);
        if (u < v) {
            v -= u;
        } else {
            uint64_t t = u - v;
            u = v;
            v = t;
        }
    } while ((v != 1) && (v != 0));
    return u << shift;
}

為了進行實驗，我需要一些隨機的整數輸入 gcd() ， C 常用的 rand() 最大只到 32767 。考量兩種計算 gcd 的方法可能在大數 cases 有更顯著的差距，因此我參考了 stack overflow ，使用 C++ 來產生亂數，程式碼。

我打算分成兩組輸入數值範圍做比較：

0 ~ 65535
65536 ~ 4294967295

這個測試方法好像怪怪的，因為隨機出來的 pair 大部分都會互質，但是我不確定什麼樣的測試方式更好

參考共筆在 user space 量測時間的方式，撰寫用來測量的程式碼。

輸入

$ sudo taskset -c 0 ./3

進行量測

結果

每組 100000 個數據進行測量，得到結果：

	small numbers	big numbers
original	27944393	44117441 (ns)
improved	16084115	23787207 (ns)
the ratio	1.73739	1.85467

small number 組：

big number 組：

顯示出利用 __builtin_ctz() 優化後的版本比原版加快將近一倍，看來 shift operation 的運算量在原版程式中佔了很大的比例；而 big number 組的差距並不比 small number 組來得顯著，與我原先猜想的不同。

解釋 kernel 中的實作

linux/lib/math/gcd.c

測驗 4

解釋程式碼

首先看 naive 版














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

因為看不出這到底在做什麼，直接丟數字進去 run ：

    uint64_t bitmap[2] = {0xFF22, 0x033F};
    uint32_t out[16];
    printf("%d\n", naive(bitmap, 2, out));

得到








(gdb) run
Starting program: /home/tinyynoob/code/linux2022-quiz2/4 
18

Breakpoint 1, main () at 4.c:25
25          return 0;
(gdb) p out
$1 = {1, 5, 8, 9, 10, 11, 12, 13, 14, 15, 64, 65, 66, 67, 68, 69}

與輸入值 {0000 0011 0011 1111, 1111 1111 0010 0010} (左邊為高位) 一起觀察，發現 out 算出了每一個 set bit 的

\log_{2} ()

。而 gdb 第 3 行的 18 則是算出 set bit 的總數。理解程式功能後往下繼續看 improved 版：
















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

由於要使用 ctz() 函式，因此每次計算完低位後就要把它設為 0 ，讓更高位可以繼續使用 ctz() ，以 00001010 舉例來說，首先 out[0] 會被設為 1 ，接著為了讓 out[1] 被設為 3 ，必須在算完 out[0] 之後將 bitset 改為 00001000 ，因此 EXP6 應為 1u << __builtin_ctzll(bitset) ，如果這裡能把第 10 行往前顯然更好。

實驗比較

實驗程式碼 (配合 Makefile 會輕鬆許多)

我用了四種 bitmap density 進行實驗， set bit 的個數分別是 8, 24, 40, 56 ，每一組都跑 100000 個數據。得到如下結果：

8

24

40

56

暫且用 Hamming weight 來稱呼 bitmap 中 set bit 的個數。

可以看到在 weight 8 的情境下， improved 比 naive 快了將近 10 倍，不過當 Hamming weight 變大時，兩條線就越來越靠近，此結果佐證了作業說明中的：

若 bitmap 越鬆散 (即 1 越少)，於是 improved 的效益就更高。

可以見到當 weight = 56 時， improved 雖然仍比 naive 更好，但已經沒有顯著的改善了。

improved() 的效能隨著 Hamming weight 增大而隨之掉落符合我們的預期，值得注意的是在另一方面，我們可以發現 naive() 在 weight 56 的效能好於 weight 40 ，甚至好於 weight 8 ，目前不曉得該如何解釋這個現象

待

進一步改進

測驗 7

探討 faster remainder

https://lemire.me/blog/2019/02/08/faster-remainders-when-the-divisor-is-a-constant-beating-compilers-and-libdivide/

Let

n, d

be constants.
Suppose that we want to do

n \div d

. we may consider

n \div d = n \times \frac{2^{N}}{d} \times \frac{1}{2^{N}}

If we can precompute

2^{N} / d

, then the division can be replaced by multiplication and bit-shift.

由於整數運算到小數點後就捨去，所以我嘗試用一些亂數跑跑看 accuracy 和 precision，整體看起來 precision 並不高，甚至在少部份的情況會出現低 accuracy 的結果：

accuracy and precision

接著我們嘗試分析 is_divisible() 函式

uint64_t c = 1 + UINT64_C(0xffffffffffffffff) / d;

// given precomputed c, checks whether n % d == 0
bool is_divisible(uint32_t n) {
  return n * c <= c - 1; 
}

簡記 UINT64_MAX 為 MAX

RHS = c - 1
    = 1 + MAX / d - 1
    = MAX / d

LHS = n * c
    = n * (1 + MAX / d)
    = n + MAX / d * n
    = n + MAX / d * (n / d * d + n % d)
    = n + (MAX * (n / d)) + (MAX / d) * (n % d)
    = n - (n / d) + (MAX / d) * (n % d)

因此，若

n, d

滿足

$n - n / d \leq MAX / d$

則有：

{\begin{cases} L H S \leq R H S & if d divides n \\ L H S > R H S & if not \end{cases}

假設

n, d

的型態皆為 uint32_t，該條件可被滿足。

補完程式碼

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"[(8 >> div5) >> (KK3)], length);
        fmt[length] = '\0';

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

直接列出全部四種情況：

(div3 && div5)
- start = 0
- length = 8
(div3)
- start = 0
- length = 4
(div5)
- start = 4
- length = 4
Otherwise
- start = 8
- length = 2

於是可以輕易推得：KK1 = div3, KK2 = div5, KK3 = 4 * div3

Jim Huang

2022/04/07 05:04:33

就在 Linux 核心原始程式碼內附的文件: https://www.kernel.org/doc/Documentation/memory-barriers.txt (Edited)

2022/04/07 06:48:56

我們的 moving average 是一個低通濾波器

TODO: 探討 CFS ( kernel/sched/fair.c ) 裡頭的低通濾波器設計 (Edited)

2022/04/07 06:49:38

Update Task's estimated utilization (Edited)

2022q1 Homework2 (quiz2)

測驗 1

解釋程式碼

編譯器最佳化下的比較

研讀 kernel 中的 EWMA

探討 kernel/sched/fair.c 中的 Low-pass filter

測驗 3

解釋程式碼

透過 __builtin_ctz 改寫並比較效能

結果

解釋 kernel 中的實作

測驗 4

解釋程式碼

實驗比較

8

24

40

56

進一步改進

測驗 7

探討 faster remainder

補完程式碼

Read more

讀 CS:APP 簡中版

學習 Linux 流水筆記

2022q1 Homework3 (fibdrv)

fibdrv 2

透過 `__builtin_ctz` 改寫並比較效能