2022q1 Homework2 (quiz2)

contributed by ???

第 2 周測驗題

測驗`1` 兩個無號整數取平均值

說明及補完程式碼

考慮以下對二個無號整數取平均值的程式碼:
```
#include <stdint.h>
uint32_t average(uint32_t a, uint32_t b)
{
    return (a + b) / 2;
}
```
這個直覺的解法會有 overflow 的問題，若我們已知 a, b 數值的大小，可用下方程式避免 overflow:
```
#include <stdint.h>
uint32_t average(uint32_t low, uint32_t high)
{
    return low + (high - low) / 2;
}
```
接著我們可改寫為以下等價的實作:
```
#include <stdint.h>
uint32_t average(uint32_t a, uint32_t b)
{
    return (a >> 1) + (b >> 1) + (EXP1);
}
```
其中 EXP1 為 a, b, 1 (數字) 進行某種特別 bitwise 操作，限定用 OR, AND, XOR 這三個運算子。

兩個整數 a,b 做 / 運算，在 C 語言當中為整數除法，e.g. 5 / 2 = 2.
而對兩數分別作 >> 運算，形同對兩數作 / 2 運算。因此若是遇到奇數，作除法後小數會被捨去。本來 （a+b）/ 2 只捨去一次，現在分別對a,b 作除法，最多可能被捨去兩次，因此 EXP1 是為了補足被捨去超過一次的情況。
在 a/2 + b/2 之下，探討 a,b 的奇偶：

a 的最右邊bit b 的最右邊bit 小數被捨去次數期望補上的值

0 0 0 0

1 0 1 0

0 1 1 0

1 1 2 1

期望補上的值，可以看出是 a 的最右邊bit & b 的最右邊bit ，取得 a 的最右邊bit 為 a & 1 ，因此 EXP1 = (a & 1) & (b & 1) = a & b & 1

參考 laneser
我們再次改寫為以下等價的實作:
```
uint32_t average(uint32_t a, uint32_t b)
{
    return (EXP2) + ((EXP3) >> 1);
}
```
其中 EXP2 和 EXP3 為 a 和 b 進行某種特別 bitwise 操作，限定用 OR, AND, XOR 這三個運算子。
EXP2 和 EXP3 必定包含 a 和 b 字元

參考你所不知道的 C 語言：數值系統篇其中的「算術完全可用數位邏輯」節，a + b 可拆成兩個運算，a & b 算的是進位（carry），a ^ b 算的是相加但不進位，可用以下真值表說明：

a b a&b a^b

0 1 0 1

1 0 0 1

0 0 0 0

1 1 1 0

因此 a + b = ((a & b) << 1) + (a ^ b)，而根據題目，要求的是 (a + b)/2 ，所以將上述式子 / 2：(((a & b) << 1) + (a ^ b)) / 2 = (((a & b) << 1) + (a ^ b)) >> 1 = (a & b) + ((a ^ b) >> 1)，所以 EXP2 = a & b ， EXP3 = a ^ b

a 的最右邊bit	b 的最右邊bit	小數被捨去次數	期望補上的值
0	0	0	0
1	0	1	0
0	1	1	0
1	1	2	1

a	b	a&b	a^b
0	1	0	1
1	0	0	1
0	0	0	0
1	1	1	0

TODO:

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

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

進入測驗 `2` 前，先來看取絕對值與 min / MAX 問題

參考解讀計算機編碼其中的「常數時間實作」篇，提到絕對值的 bit-wise 操作：

a b s (x) = {\begin{cases} x & if x \geq 0 \\ (\sim x) + 1 = \sim (x - 1) & if x < 0 \end{cases}

對應的實作可以為 (注意此實作 x 不能為 INT32_MIN )：

#include <stdint.h>
int32_t abs(int32_t x) {
    int32_t mask = (x >> 31);
    return (x + mask) ^ mask; //if x < 0, return ~(x-1), x 對 -1 作 ^, 形同 x 反轉位元
}

另外也提到給定兩個有號整數，找出其中較小的數值（注意這實作僅在 INT32_MIN <= (a - b) <= INT32_MAX 才成立）：

#include <stdint.h>
int32_t min(int32_t a, int32_t b) {
    int32_t diff = a - b;
    return b + (diff & (diff >> 31));
}

若 a - b = diff < 0，則 a 較小，diff >> 31 為 0xffffffff，與 diff & 運算後仍為 diff，回傳 b + diff = a
若 a - b = diff > 0，則 b 較小，diff >> 31 後為0x00000000，與 diff & 運算後為 0x00000000，回傳 b + 0x00000000 = b

測驗`2` 找 MAX

說明及補完程式碼

改寫〈解讀計算機編碼〉一文的「不需要分支的設計」一節提供的程式碼 min，我們得到以下實作 (max):
```
#include <stdint.h>
uint32_t max(uint32_t a, uint32_t b)
{
    return a ^ ((EXP4) & -(EXP5));
}
```
其中 EXP4 為 a 和 b 進行某種特別 bitwise 操作，限定用 OR, AND, XOR 這三個運算子。
EXP5 為 a 和 b 的比較操作，亦即 logical operator，限定用 >, <, ==, >=, <=, != 這 6 個運算子。
EXP4, EXP5 必定包含 a 和 b 字元。

EXP5 前有個負號，直覺 EXP5 應為 0 或 1，而若要找較大者，必定是比較 a, b 大小，假設 EXP5 = a < b，再探討下述兩種情況：
a > b：應該回傳 a = a ^ 0（a 對 0 xor，形同無作用），所以 0 = (EXP4) & -(EXP5) = (EXP4) & 0 –- 一式。
a < b：應該回傳 b = a ^ (a ^ b )（b 對 a 作兩次 xor 形同沒作用），所以 a ^ b = (EXP4) & -(EXP5) = (EXP4) & -1 –- 二式。
解二式 EXP4 = a ^ b。代入一式 (a ^ b) & 0 也確實等於 0 。
因此答案為：
```
return a ^ ((a ^ b) & -(a < b)); 
```

TODO:

針對 32 位元無號/有號整數，撰寫同樣 branchless 的實作
不論是有號或無號上述的方法都可以使用，只需要將型別改成有號即可：uint32_t -> int32_t

測驗`3` GCD

說明及補完程式碼

考慮以下 64 位元 GCD (greatest common divisor, 最大公因數) 求值函式:

#include <stdint.h>
uint64_t gcd64(uint64_t u, uint64_t v)
{
    if (!u || !v) return u | v;//回傳非 0 的參數
    /*輾轉相除法概念*/
    while (v) {                               
        uint64_t t = v;
        v = u % v;
        u = t;
    }
    return u;
}

註：GCD 演算法
1. If both
  $x$ and
  $y$ are
  $0$ , gcd is
  $0$ .
  $g c d (0, 0) = 0$ ,
  $g c d (0, 0) = 0$ .
2. $g c d (x, 0) = x$ and
  $g c d (0, y) = y$ , because everything divides
  $0$ .
3. 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.
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.

改寫為以下等價實作:

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

第一個 if 為上述註的第一及第二點。
執行 for loop 的條件為 u,v 兩數皆為偶數（u, v 任一數二進制的最低位元為 1，條件便無法成立），若皆為偶數兩數都 /= 2。此為上述註的第三點。
while loop 執行的為上述註的第四點。（
$x$ ,
$y$ 對應 u,v）
do while loop 當中，第一個 while loop 執行的是註的第五點。接下來的 if else 為輾轉相除法的實作，結束條件為當較小的數為 0 的時後，這裡較小的數為 v，因此 COND = v，且當較小數為 0 時，最大公因數即為另一較大的數。但需要注意的是，在執行上面的 for loop 時，有計算多少 2 的冪 = shift ，使得
$2^{s h i f t}$ 為公因數，因此 RET = u << shift。

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

參考你所不知道的 C 語言：數值系統篇其中的「Count Leading Zero」節，__builtin_ctz 算的是一數的二進制表示法中，從最低位開始遇到第一個 1 之前，共有幾個 0。
知道一數的 ctz（count tailing zero）後，便可以知道其中的因數為 2的多少次冪。
以相同上述實作的邏輯改寫後：

































#include <stdint.h>
uint64_t mygcd64(uint64_t u, uint64_t v) 
{
    if (!u || !v) return u | v;

    int ctz_u, ctz_v, shift;
    ctz_u = __builtin_ctz(u);
    ctz_v = __builtin_ctz(v);

    if (ctz_u > ctz_v) {
        shift = ctz_u - (ctz_u - ctz_v);
    } else {
        shift = ctz_v - (ctz_v - ctz_u);
    }
    u = u >> shift;
    v = v >> shift;

    if (ctz_u = __builtin_ctz(u)) u = u >> ctz_u;

    do {
        if (ctz_v = __builtin_ctz(v)) v = v >> ctz_v;

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

    return u << shift;
}

第 10 到第 16 行，做的是上份實作的 for loop 的功能，算出公因數： 2 的多少次冪。
第 18 及 21 行，執行的是上份實作兩個 while loop 的功能，若有偶數，將 2 的因數除掉。

參考Linux效能分析工具：Perf，根據以下 main function ：

int main() 
{
uint64_t u = 64, v = 52;

printf("%ld\n",gcd64(u,v));
//printf("%ld\n",mygcd64(u,v));

return 0;
}

改寫前實作效能如下：

$ perf stat --repeat 5 -e cache-misses,cache-references,instructions,cycles ./quiz

> Performance counter stats for './quiz' (5 runs):

    22,466      cache-misses              #   35.419 % of all cache refs      ( +-  4.84% )
    63,429      cache-references                                              ( +-  1.88% )
   825,880      instructions              #    0.56  insn per cycle           ( +-  0.40% )
 1,474,393      cycles                                                        ( +-  3.79% )

 0.0007654 +- 0.0000326 seconds time elapsed  ( +-  4.26% )

改寫後實作效能如下：

$ perf stat --repeat 5 -e cache-misses,cache-references,instructions,cycles ./mine

> Performance counter stats for './mine' (5 runs):

    23,806      cache-misses              #   36.953 % of all cache refs      ( +-  2.95% )
    64,423      cache-references                                              ( +-  5.30% )
   830,068      instructions              #    0.69  insn per cycle           ( +-  0.50% )
 1,199,807      cycles                                                        ( +-  5.34% )

  0.001786 +- 0.000135 seconds time elapsed  ( +-  7.53% )

兩者數據結果差不多，接著把執行次數提升到 100，並多偵測 branch event：

$ perf stat --repeat 100 -e cache-misses,cache-references,instructions,cycles,branch-misses,branch-instructions ./quiz

> Performance counter stats for './quiz' (100 runs):

    22,362      cache-misses              #   40.160 % of all cache refs      ( +-  1.78% )
    55,682      cache-references                                              ( +-  0.76% )
   830,198      instructions              #    0.62  insn per cycle           ( +-  0.10% )
 1,339,491      cycles                                                        ( +-  1.86% )
     6,600      branch-misses             #    3.91% of all branches          ( +-  0.39% )
   168,835      branch-instructions                                           ( +-  0.10% )

 0.0007404 +- 0.0000235 seconds time elapsed  ( +-  3.17% )

$ perf stat --repeat 100 -e cache-misses,cache-references,instructions,cycles,branch-misses,branch-instructions ./mine

> Performance counter stats for './mine' (100 runs):

    22,054      cache-misses              #   40.183 % of all cache refs      ( +-  1.38% )
    54,885      cache-references                                              ( +-  0.61% )
   830,690      instructions              #    0.65  insn per cycle           ( +-  0.12% )
 1,270,811      cycles                                                        ( +-  1.71% )
     6,630      branch-misses             #    3.92% of all branches          ( +-  0.42% )
   168,972      branch-instructions                                           ( +-  0.11% )

 0.0006431 +- 0.0000144 seconds time elapsed  ( +-  2.23% )

兩份實作 cache miss 約 40.1%，branch miss 約 3.9％，數據結果差異不大，改寫後的程式尚須精進。

TODO:

為何改寫後效能仍尚未改進？
Linux 核心中也內建 GCD (而且還不只一種實作)，例如 lib/math/gcd.c，請解釋其實作手法和探討在核心內的應用場景。

測驗`4` 計算 bitmap 中有多少個1, 及分別在第幾位

說明及補完程式碼

在影像處理中，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;
}

此段程式中，pos 代表整個 bitmap 中有多少個 1。out 紀錄的是每個為 1 的 bit 在 bitmap 的第幾位，這裡可以看出即使 bitmap 為 64 * bitmapsize，在儲存空間上仍是一維的。

考慮 GNU extension 的 __builtin_ctzll 的行為是回傳由低位往高位遇上連續多少個 0 才碰到 1。
用以改寫的程式碼如下：
















#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 的效益就更高。
Hint：需考慮到 -bitwise 的特性

根據此目的，可以先算出 bitset 的 ctz，把 1L 左移 ctz 個，如此可以得出只保留最低位 1 的 bitset。因此 EXP6 = 1L << __builtin_ctzll(bitset)。

但若要利用 -bitwise 的特性，需要換個作法：若把 bitset 取負，與 bitset 差異為除了最低位的 1 以及更低位的 0 相同之外，其他位皆相反。因此若要得到只保留最低位 1 的 bitset，把 bitset & -bitset 即可。EXP6 = bitset & -bitset。

第 12 行做的便是把 bitset 最低位 1 消除，讓下一個迴圈再次去找最低位 1 在哪。

TODO:

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

2022q1 Homework2 (quiz2)

測驗1 兩個無號整數取平均值

說明及補完程式碼

進入測驗 2 前，先來看取絕對值與 min / MAX 問題

測驗2 找 MAX

說明及補完程式碼

測驗3 GCD

說明及補完程式碼

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

測驗4 計算 bitmap 中有多少個1, 及分別在第幾位

說明及補完程式碼

Read more

2022q1 Homework3 (fibdrv)

2022q1 Homework5 (quiz5)

2022q1 Homework1 (lab0)

2022q1 Homework1 (quiz1)

測驗`1` 兩個無號整數取平均值

進入測驗 `2` 前，先來看取絕對值與 min / MAX 問題

測驗`2` 找 MAX

測驗`3` GCD

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

測驗`4` 計算 bitmap 中有多少個1, 及分別在第幾位