2022q1 Homework2 (quiz2)

contributed by yyyyuwen

測驗題

測驗1

題目要求：求兩整數之平均值，並將有可能 overflow 的程式碼改成避免 overflow 的形式。

以下為有可能造成 overflow 的程式碼：






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

版本1





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

注意：確保大小數的順序，不然也有可能會導致輸出結果不同的問題。

當 (high - low) < 0 時，值會是
$⌈ (h i g h - l o w) ⌉$
當 (high - low) > 0 時，值會是
$⌊ (h i g h - l o w) ⌋$

版本2






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

其中 EXP1 為 a, b, 1 (數字) 進行某種特別 bitwise 操作，限定用 OR, AND, XOR 這三個運算子。

a >> 1 和 b >> 1 的意思是將數值除以2。而 a 和 b 都是 uint32_t 的型態，若是奇數無法整除小數會直接被省去。因此若遇上這種狀況時，做完 (a >> 1) + (b >> 1) 之後還要再加上 1。
a & 1 ：檢查 a 的 Least Significant Bit 是否為 1 （即檢查是否為奇數），若是則回傳 1
。
EXP1 ：根據以上推論，需要檢查 a 和 b 是否為奇數，填入 a & b & 1 。

舉例：

a = 5 = 0101, b = 7 = 0111
取平均為 6 -> 0110
a >> 1 = 0010
b >> 1 = 0011
0010 | 0011 = 0011 != 0111

版本3




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

其中 EXP2 和 EXP3 為 a 和 b 進行某種特別 bitwise 操作，限定用 OR, AND, XOR 這三個運算子。
參考加法運算

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 →

a + b 可以變成是 a XOR b + (a AND b) << 1

a XOR b ： a + b 但不進位
(a AND b) << 1 ：檢查 a, b Most Significant Bit是不是同為 1，如果是則進位。

SUM ：取得未進位時的SUM

a	b	a ^ b	a + b
0	0	0	$00_{b}$
0	1	1	$01_{b}$
1	0	1	$01_{b}$
1	1	0	$10_{b}$ (需要進位)

CARRY ：取得是否進位的資訊

a	b	a & b	a + b
0	0	0	$00_{b}$ -> 不進位
0	1	0	$01_{b}$ -> 不進位
1	0	0	$01_{b}$ -> 不進位
1	1	1	$10_{b}$ -> 進位

(a + b) / 2
式子分解：
= ( a + b ) >> 1
= (( a ^ b ) + ( a & b ) << 1) >> 1
= a & b + ( a ^ b ) >> 1

EXP2 ： a & b
EXP3 ： a ^ b

再做實驗時遇到一個問題：當用 uint32_t的時候做 (a + b)/2會 overflow 沒問題，但當如果使用 uint8_t 的時候就不會有 overflow的問題。目前還沒找到原因。

延伸問題:

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

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

測驗2

改寫〈解讀計算機編碼〉一文的「不需要分支的設計」一節提供的程式碼 min，並實作 (max):

以下是該文章提出的 min 程式碼實作




int32_t min(int32_t a, int32_t b) {
    int32_t diff = (a - b);
    return b + (diff & (diff >> 31));
}

注意：這實作僅在 INT32_MIN <= (a - b) <= INT32_MAX 成立時，才會發揮作用。

而我們要實作 max 的方法





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

XOR 的幾個重要功用

x ^ 0 = x
x ^ x = 0
x ^ y ^ x = y

以下為了方便，將 ((EXP4) & -(EXP5)) 設成 x
觀察題目，有兩種狀況：

a >= b
此時應該會 return a ，則算式應為 a ^ 0 = a -> x = 0
a < b
這個狀況則需要 return b ，算式應為 a ^ a ^ b = b -> x = a ^ b

依上述狀況可以判斷

EXP4 可以填入 a ^ b
EXP5 作為判斷 0 或 1
- -(EXP5) 是 0 (
  $00000000000000000000000000000000_{b}$ ) -> a > b
- -(EXP5) 是 -(1) (
  $11111111111111111111111111111111_{b}$ ) -> a < b
- 因此 EXP5 填入 a < b





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

延伸問題:

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 檢索。

測驗3

原本程式











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

改寫成























#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 演算法

If both x and y are 0, gcd is zero
$g c d (0, 0) = 0$ .
$g c d (x, 0) = x$ and
$g c d (0, y) = y$ because everything divides 0.

if (!u || !v)
    return u | v;

If x and y are both even,
$g c d (x, y) = 2 * g c d (\frac{x}{2}, \frac{y}{2})$ because
$2$ is a common divisor. Multiplication with
$2$ can be done with bitwise shift operator.

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

(u | v) & 1 判斷是否同為偶數，shift 紀錄總共做了幾次 right shift 的運算 ( <<1 )。
4. If x is even and y is odd,

g c d (x, y) = g c d (\frac{x}{2}, y)

.
5. On similar lines, if x is odd and y is even, then

g c d (x, y) = g c d (x, \frac{y}{2})

. It is because

2

is not a common divisor.

確保 u 為奇數：

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

輾轉相除法：

if (u < v) {
    v -= u;
} else {
    uint64_t t = u - v;
    u = v;
    v = t;
}

利用 while 去遞迴的找餘數，使用連續減法來找出相除的餘數。
而根據輾轉相除法的定理，中止條件則是除到餘數為0。
在上述程式中
u = v 是商數的部份，而 v = t 是餘數的部份，所以 COND 指的就是餘數的部份也就是當 v 等於0的時候變中止迴圈；因此 COND = v 。
而 u 則是最後的最大公因數，但記得要乘回一開始 shift 的部份，真正的最大公司數應為

2^{s h i f t} * g c d (u, v)

。
因此 RET = u << shift 。

延伸問題:

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

測驗4

題目：在影像處理中，bit array (也稱 bitset) 廣泛使用，考慮以下程式碼:














#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 extenstion 的 __builtin_ctzll的行為是回傳由低位往高位遇上連續多少個 0 才碰到 1。

範例: 當 a = 16
16 這個十進位數值的二進位表示法為 00000000 00000000 00000000 00010000
從低位元 (即右側) 往高位元，我們可發現 0 → 0 → 0 → 0 → 1 ，於是 ctz 就為 4，表示最低位元往高位元有 4 個 0

用以改寫的程式碼如下:
















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

其中第 9 行的作用是找出目前最低位元的 1，並紀錄到 t 變數。舉例來說，若目前 bitset 為 000101_b，那 t 就會變成 000001_b，接著就可以透過 XOR 把該位的 1 清掉，其他保留 (此為 XOR 的特性)。

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

請補完程式碼。書寫規範:

bitwise 運算子和運算元之間用一個半形空白區隔，如: bitset | 0x1
考慮到 -bitwise 的特性
不要包含任何小括號，即 ( 和 )
以最精簡的形式撰寫程式碼

解釋

可以從題目得知該行的目的是為了找出最低位元的1，可以透過二補數的原理來達成目的。即 bitset & -bitset。
舉例

bitset = 6 = 0110
-bitset = 1010
則
bitset & -bitset = 0010 //剩下的即是最低位元的1

關於其他程式碼的解釋：
k 代表的是64位 bitset 在整個 bitmap 中的編號
r 則是 bitmap 中末端為零的數量

延伸問題:

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