2022q1 Homework2 (quiz2)

contributed by < freshLiver >

tags: `linux2022`

第一題 - 兩正整數之平均值

直覺的寫法

取平均值最直覺的方法當然是相加後除以二：

uint32_t average(uint32_t a, uint32_t b)
{
    return (a + b) / 2;
}

但上面的會有 Integer Overflow 的問題，因此可以用以下方式改寫來避免該問題：

uint32_t average(uint32_t low, uint32_t high)
{
    return low + (high - low) / 2;
}

由於兩數（假設為 low 以及 high ）的平均值就是在兩數間找到另一數 mid，使得 low 以及 high 與 mid 的差值相等（若 high、low 兩數和為奇數的話，mid 會離 low 比較近一點）。

而找中點 mid 的操作可以從 high、low 兩數的差除 2，然後再加到 low 上，這樣就可以避開將 high、low 相加而造成 Integer Overflow 的問題。

問題 1. 需要注意參數的大小順序

$ cat main.c && gcc main.c && ./a.out 200 100
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>

uint32_t average(uint32_t low, uint32_t high)
{
    return low + (high - low) / 2;
}

int main(int argc, char const *argv[])
{
    uint32_t a = atoi(argv[1]);
    uint32_t b = atoi(argv[2]);
    printf("Average is : %u", average(a, b));
    return 0;
}
Average is : 2147483798

如上面的測試結果，由於 low 跟 high 是無號整數，所以在 (high - low) 的部份如果 low 的值比 high 大的話會造成結果與預期不符。

問題 2. 相較於加法以及乘法，除法需要的 cycle 數較多

不用除法的解法 - 1 EXP1

uint32_t average(uint32_t a, uint32_t b)
{
    // EXP1
    return (a >> 1) + (b >> 1) + (a & b & 1);
}

這種寫法與 (a + b) / 2 比較相近，是先將 a、b 兩正整數右移 1 位元來達成除 2 的效果。但直接右移會損失掉最後一位元，而最後一位元的和卻可能會造成進位，因此需要另外計算是否進位。

       0 1 0 0 1 1 1 d (a)
    &) 1 1 1 1 1 1 0 c (b)
    ------------------
       0 1 0 0 1 1 0 ? (sum)

而最低位元的進位會發生在 c 與 d 皆為 1 的時候，因此可以用 AND 檢查兩位元是否皆為 1，但是需要檢查的位元只有最低位元，因此可以利用 1 除了最低位元以外都是 0 的特性，與 a & b 的結果進行 AND 捨去其他位元的資訊：

       0 1 0 0 1 1 0 ? (sum)
    &) 0 0 0 0 0 0 0 1 
    ------------------
       0 0 0 0 0 0 0 ?

綜合以上可以知道 EXP1 是要處理進位的部份，可以透過 a & b & 1 完成

不用除法的解法 - 2 EXP2, EXP3

上面的方法用了多個 +、>> 以及 &，但還能再簡化：

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

與前一個方法不同，這個寫法不只是把最後一位元的進位另外計算，而是直接使用 AND 以及 XOR 找出需要進位以及不須進位的位元，並透過位移同時完成除二以及進位的效果。

        1 1 0 0 1 1 0 1 (a, 205)
    ^ ) 1 1 0 1 0 0 1 1 (b, 211)
    -------------------
        0 0 0 1 1 1 1 0 (XOR result, 不須進位位元的和)

透過對 a、b 進行 XOR 操作能夠找出 0 與 1 各有一個的位元，而這些位元就是相加後不須進位的位元，像是上方的第 4, 3, 2, 1 個位元。

        1 1 0 0 1 1 0 1 (a)
    & ) 1 1 0 1 0 0 1 1 (b)
    -------------------
        1 1 0 0 0 0 0 1 (AND result, 會引起進位的位元)

而如同上一個寫法，透過 a 與 b 進行 AND 則可以找出會造成進位的位元，像是上方的第 7, 6, 0 個位元。

而有了這兩個資訊的話，只要將 AND 的結果左移就能夠知道 a b 兩數相加後的值了：

     左移消失
        v
        1 1 0 0 0 0 0 1 0  (AND result << 1)
    + )   0 0 0 1 1 1 1 0  (XOR result)
    ---------------------
          1 0 1 0 0 0 0 0  (160, 錯誤的相加結果)

但很明顯的這會造成最高位元的進位消失，因此在這個解法中反而是將 XOR 的結果右移，然後再加上 AND 的結果：

                     右移消失
                        v
        1 1 0 0 0 0 0 1   (AND result)
    + ) 0 0 0 0 1 1 1 1 0 (XOR result >> 1)
    -------------------
        1 1 0 1 0 0 0 0   (208)

這樣雖然會造成 XOR 結果的最低位元消失，但其實最後一位元本來就只需要是否進位這個資訊而已，因此這個作法不僅能夠正確的保留資訊，還能完成兩數取平均（除 2）的要求。

比較以上方法在「編譯器最佳化」開啟時的差異

以下是以上三個方式使用 GCC 開啟 -O3 最佳化時的組合語言輸出 (程式碼以及原始輸出)：

使用除法

uint32_t average_div(uint32_t low, uint32_t high)
{
    return low + (high - low) / 2;
}

average_div:
        mov     eax, esi
        sub     eax, edi
        shr     eax
        add     eax, edi
        ret

將存有 low 的暫存器 esi 的值給 eax 暫存器
將 eax 與存有 high 的暫存器 edi 相加並存到 eax
將 eax 右移一位來達成 eax / 2
將 eax 與存有 low 的暫存器 esi 相加得到最終值

Bitwise

uint32_t average_bw1(uint32_t a, uint32_t b)
{
    return (a >> 1) + (b >> 1) + (a & b & 1);
}

average_bw1:
        mov     eax, edi
        mov     edx, esi
        and     edi, esi
        shr     eax
        shr     edx
        and     edi, 1
        add     eax, edx
        add     eax, edi
        ret

分別把存有 a 跟 b 的暫存器 edi、esi 存進 eax、edx
先對 edi、esi 進行 AND 運算，運算結果會存在 edi
分別對 eax、edx 進行 shr 右移來達到 / 2 的效果
將步驟二的結果與 1 進行 AND 運算，得到最低位元是否須進位
將步驟三的結果相加並存在 eax，得到除二後的相加結果
將儲存步驟四和五結果的 eax、edi 相加可得到結果

這個寫法若是使用 -O1 參數進行最佳化，會得到下方的組合語言：

average_bw1:
        mov     eax, edi
        shr     eax
        mov     edx, esi
        shr     edx
        add     eax, edx
        and     edi, esi
        and     edi, 1
        add     eax, edi
        ret

可以發現和 -O3 最佳化輸出的組合語言相比，執行的指令都沒變，只差在執行順序不同， -O1 輸出的組合語言是依照程式碼由左至右執行，而 -O3 則是將相同的指令放在一起，

uint32_t average_bw2(uint32_t a, uint32_t b)
{
    return (a & b) + ((a ^ b) >> 1);
}

average_bw2:
        mov     eax, edi
        and     edi, esi
        xor     eax, esi
        shr     eax
        add     eax, edi
        ret

將存有 a 的 edi 複製到 eax 並與存有 b 的 esi 進行 XOR 運算來取得不須進位的位元
將運算結果透過 shr 右移來除二
將存有 a 與 b 的 edi esi 進行 AND 運算找出需要進位的位元，並存在 edi 中
將需要進位的位元 edi 與第二步的結果 eax 相加獲得結果

若是將這種寫法用 -O1 進行最佳化，也會像前一種寫法一樣，出現指令順序被重排的情況：

average_bw2:
        mov     eax, edi
        xor     eax, esi
        shr     eax
        and     edi, esi
        add     eax, edi
        ret

TODO : 在自己電腦編譯後產出的組合語言內容不同

	.file	"main.c"
	.text
	.p2align 4
	.globl	average_div
	.type	average_div, @function
average_div:
.LFB0:
	.cfi_startproc
	endbr64
	movl	%esi, %eax
	subl	%edi, %eax
	shrl	%eax
	addl	%edi, %eax
	ret
	.cfi_endproc
.LFE0:
	.size	average_div, .-average_div
	.p2align 4
	.globl	average_bw1
	.type	average_bw1, @function
average_bw1:
.LFB1:
	.cfi_startproc
	endbr64
	movl	%edi, %eax
	movl	%esi, %edx
	andl	%esi, %edi
	shrl	%eax
	shrl	%edx
	andl	$1, %edi
	addl	%edx, %eax
	addl	%edi, %eax
	ret
	.cfi_endproc
.LFE1:
	.size	average_bw1, .-average_bw1
	.p2align 4
	.globl	average_bw2
	.type	average_bw2, @function
average_bw2:
.LFB2:
	.cfi_startproc
	endbr64
	movl	%edi, %eax
	andl	%esi, %edi
	xorl	%esi, %eax
	shrl	%eax
	addl	%edi, %eax
	ret
	.cfi_endproc
.LFE2:
	.size	average_bw2, .-average_bw2
	.ident	"GCC: (Ubuntu 10.3.0-1ubuntu1) 10.3.0"
	.section	.note.GNU-stack,"",@progbits
	.section	.note.gnu.property,"a"
	.align 8
	.long	 1f - 0f
	.long	 4f - 1f
	.long	 5
0:
	.string	 "GNU"
1:
	.align 8
	.long	 0xc0000002
	.long	 3f - 2f
2:
	.long	 0x3
3:
	.align 8
4:

比較效能

探討 Linux 核心中的 EWMA 實作原理及應用

第二題 - branchless `max(a,b)`

uint32_t max(uint32_t a, uint32_t b)
{
    return a ^ ((EXP4) & -(EXP5));
}

從題目可知 EXP5 是 a、b 進行大小比較，因此 -(EXP5) 的值不是 0 就是 -1。而一數值與 0 進行 AND 運算的結果必為 0、與 -1 進行 AND 運算的結果則會與原數值相同。

到這邊為止，可以知道這個函式的結果必為 a ^ (EXP4) 或 a ^ 0，所以可以依照 a、b 的大小關係推測 EXP4, EXP5 的最簡表示應分別為：

$a > b$ : -EXP5 或 EXP4 為 0 才能使結果為 a
$a = b$ : -EXP5 或 EXP4 為 0 才能使結果為 a
$a < b$ : EXP4 應為 a ^ b 且 -EXP5 不可為 0 才能使結果為 b

由於 EXP4 應為 a ^ b，所以 -EXP5 在

a \geq b

時會是 0，而

a < b

時則為 -1，因此 EXP5 應為 a < b 或 a <= b。

實作 32 位元有號整數版

這個函式應該也適用於有號數，不太確定要求是什麼

找出 Linux Kernel 中的 branchless 設計

在 include/linux/bitops.h 有中包含了幾個 rotate 的操作函式，像是 rol32、ror64 以及其他名稱相近的函式：

static inline __u8 rol8(__u8 word, unsigned int shift)
{
	return (word << (shift & 7)) | (word >> ((-shift) & 7));
}

以 rol8 為例，這個函式可以把 word 最高位的 shift 個位元移動到最低 shift 個位元，就像是 word 的最高位與最低位相連一樣，所以若是 word 為 0b11110000 的話，當 shift 為 5 時，結果就會是 0b00011110。

實作方式是將 word 左移 shift 位元後，再與 word 右移 8 - shift 位元的結果進行 OR 組合。

比較特別的是右移 8 - shift 位元的部份是透過二補數的特性來進行：

    -shift = ~shift + 1             (二補數 = 一補數 + 1)
    ~shift = 2^n - 1 - shift
    -shift = 2^n - 1 - (shift - 1)  (第二式帶入第一式)

因此若 shift 不大於 8，(shift - 1) 也不會大於 7，則 (shift - 1) 只會影響 2^n - 1 的最低三位元，所以 ((-shift) & 7) 相當於 7 - (shift - 1) 也就是 8 - shift。

在程式碼中可以看到 (shift & 7) 與 ((-shift) & 7)，我想是因為 shift 的最大值超過 word 的總位元數，因此需要與 7 進行 AND 來避免 shift 超過 7 的情況，但不太明白為什麼要選擇 32 位元的 unsigned int 而不是用 8 位元的 __u8 來儲存 shift 的值。

第三題 - 不使用 modulo 的 GCD

求兩整數的最大公因數 (GDB) 可以透過實作輾轉相除法來完成。












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

上面的程式碼使用了迭代方式實作輾轉相除法，但 12 行的 % 本質上是除法，所以可以嘗試將 % 改用其他位元操作來加速運算：
























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

不使用 % 的輾轉相除法實作概念與使用 % 的寫法相似，都是兩數輪流對另一數取餘數，但不使用 % 的實作方式是透過連續的減法來取餘數，可以從第 12 行的 do...while 區塊看出：

當 u < v 時，就不斷將 v 減去 u，直到 u > v 時就代表原本的 v % u 完成
而當 v % u 完成時，會有兩種狀況：
1. v % u == 0：此時代表已經找到最大公因數了，所以迴圈應該要結束，因此 COND 應為 v
2. v % u != 0：此時因取餘數的操作都是以 v -= u 進行，所以第 17 行的 else 部份負責將 v 與 u 對調，而且不只是單純對調，而是將新的 v 換成 u - v 來減少一輪 v -= u 的操作。

但單純將 % 改用減法取代不一定能夠減少計算時間，所以上面的實作還利用了數次除二的計算來減少運算：

u 與 v 可以分別表示成
$a \times 2^{n}$ 以及
$b \times 2^{m}$ ，因此第 7 行的迴圈移去 u、v 最低位元的
$m i n (n, m)$ 個，並用 shift 紀錄。
在找到
$m i n (n, m)$ 後，u 與 v 就不會有 2 的公因數，因此若是 u 或 v 是二的倍數的話，可以直接移除尾端連續的 0：
- 第 10 行的 while 會先將 u 尾端連續的 0 移除，讓
  $u = a$
- 第 13 行的 while 在第一次進入時會將 v 尾端連續的 0 移除，讓
  $v = b$ 。而之後由於 u、v 可能會互換，所以第二次後進入此迴圈的意義則是將被取餘數的 v 尾端連續的 0 移除（u 尾端的 0 在互換前就被移除掉了）。

最後由於一開始就先找到 shift 並位移掉了，所以需要將 v % u == 0 時的 u 左移 shift 位元才會是正確的最大公因數，因此 RET 應為 u << shift。

使用 `__builtin_ctz` 改寫並分析對效能的提升

#define min(a, b) ((a) ^ (((a) ^ (b)) & -((a) > (b))));

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

    int utz = __builtin_ctzll(u);
    int vtz = __builtin_ctzll(v);
    u >>= utz;

    do {
        v >>= __builtin_ctzll(v);

        if (u < v) {
            v -= u;
        } else {
            uint64_t t = u - v;
            u = v;
            v = t;
        }
    } while (v);
    return u << min(utz, vtz);
}

由於 u、v 的型別是 uint64_t，而 __builtin_ctz 的參數型別是 unsigned int，所以要改用 __builtin_ctzl 或 __builtin_ctzll。

效能測試

int main(int argc, char const *argv[])
{
    uint64_t res;
    for (uint64_t a = S; a < E; ++a)
        for (uint64_t b = S; b < E; ++b)
            res = FUNC(a, b);
    return 0;
}

嘗試使用迴圈對 [S,E) 區間的數兩兩呼叫 GCD 函式進行測試，並使用 perf 對結果進行簡單的分析：

使用 `while` 檢查 ctz（無最佳化）

$ make perf FUNC=gcd64_loop START=10000 END=20000
gcc  -D FUNC=gcd64_loop -D S=10000 -D E=20000 main.c -o main.out
perf stat -r 10 ./main.out

 Performance counter stats for './main.out' (10 runs):

         13,696.83 msec task-clock                #    1.000 CPUs utilized            ( +-  0.11% )
               103      context-switches          #    0.007 K/sec                    ( +- 25.64% )
                13      cpu-migrations            #    0.001 K/sec                    ( +- 11.00% )
                48      page-faults               #    0.004 K/sec                    ( +-  0.88% )
    57,818,688,302      cycles                    #    4.221 GHz                      ( +-  0.09% )
    31,840,178,938      instructions              #    0.55  insn per cycle           ( +-  0.01% )
     6,351,903,123      branches                  #  463.750 M/sec                    ( +-  0.01% )
     1,350,948,860      branch-misses             #   21.27% of all branches          ( +-  0.01% )

           13.7005 +- 0.0145 seconds time elapsed  ( +-  0.11% )

使用 `while` 檢查 ctz（開啟 `-O3` 最佳化）

$ make perf CFLAGS=-O3 FUNC=gcd64_loop START=10000 END=20000
gcc -O3 -D FUNC=gcd64_loop -D S=10000 -D E=20000 main.c -o main.out
perf stat -r 10 ./main.out

 Performance counter stats for './main.out' (10 runs):

          8,124.18 msec task-clock                #    1.000 CPUs utilized            ( +-  0.10% )
                64      context-switches          #    0.008 K/sec                    ( +- 27.21% )
                 9      cpu-migrations            #    0.001 K/sec                    ( +- 12.48% )
                48      page-faults               #    0.006 K/sec                    ( +-  1.07% )
    34,258,837,853      cycles                    #    4.217 GHz                      ( +-  0.06% )
    18,307,940,212      instructions              #    0.53  insn per cycle           ( +-  0.02% )
     7,639,222,748      branches                  #  940.307 M/sec                    ( +-  0.01% )
     1,372,710,376      branch-misses             #   17.97% of all branches          ( +-  0.01% )

           8.12629 +- 0.00756 seconds time elapsed  ( +-  0.09% )

使用 `__builtin_ctzll` 改寫（無最佳化）

$ make perf FUNC=gcd64_ctz START=10000 END=20000
gcc  -D FUNC=gcd64_ctz -D S=10000 -D E=20000 main.c -o main.out
perf stat -r 10 ./main.out

 Performance counter stats for './main.out' (10 runs):

          7,176.32 msec task-clock                #    1.000 CPUs utilized            ( +-  0.14% )
                54      context-switches          #    0.008 K/sec                    ( +- 13.40% )
                 7      cpu-migrations            #    0.001 K/sec                    ( +- 14.32% )
                48      page-faults               #    0.007 K/sec                    ( +-  0.69% )
    30,265,247,336      cycles                    #    4.217 GHz                      ( +-  0.14% )
    18,926,047,559      instructions              #    0.63  insn per cycle           ( +-  0.00% )
     2,943,010,768      branches                  #  410.100 M/sec                    ( +-  0.00% )
       432,757,813      branch-misses             #   14.70% of all branches          ( +-  0.01% )

            7.1798 +- 0.0105 seconds time elapsed  ( +-  0.15% )

使用 `__builtin_ctzll` 改寫（開啟 `-O3` 最佳化）

$ make perf CFLAGS=-O3 FUNC=gcd64_ctz START=10000 END=20000
gcc -O3 -D FUNC=gcd64_ctz -D S=10000 -D E=20000 main.c -o main.out
perf stat -r 10 ./main.out

 Performance counter stats for './main.out' (10 runs):

          3,422.79 msec task-clock                #    1.000 CPUs utilized            ( +-  0.09% )
                31      context-switches          #    0.009 K/sec                    ( +- 25.36% )
                 4      cpu-migrations            #    0.001 K/sec                    ( +- 24.85% )
                48      page-faults               #    0.014 K/sec                    ( +-  0.82% )
    14,442,111,363      cycles                    #    4.219 GHz                      ( +-  0.08% )
    10,451,533,640      instructions              #    0.72  insn per cycle           ( +-  0.00% )
     2,429,728,556      branches                  #  709.867 M/sec                    ( +-  0.00% )
       455,213,440      branch-misses             #   18.74% of all branches          ( +-  0.01% )

           3.42427 +- 0.00298 seconds time elapsed  ( +-  0.09% )

上面是在開啟與不開啟編譯器最佳化的情況下，分別對兩種實作方式進行 10 次 1 億次（

[10000, 20000)

區間內，兩兩取最大公因數）呼叫後的平均結果，可以看到使用 __builtin_ctzll 改寫後：

項目	無最佳化	`-O3`	附註
cycle	$- 47.65 %$	$- 57.84 %$
instruction	$- 40.56 %$	$- 42.91 %$
IPC	$+ 14.54 %$	$+ 35.85 %$
branch	$- 53.67 %$	$- 68.19 %$	`while` 版開啟最佳化後 branch 數增加
branch-miss	$- 67.97 %$	$- 66.84 %$	兩種實作方式開啟最佳化後 branch-miss 皆增加
執行時間	$- 47.60 %$	$- 57.88 %$

探討 Linux 核心中 GCD 實作手法

在 lib/math/gcd.c 使用了兩種方式實作 GCD 函式，這兩個函式大致上與這題不使用 % 實作的概念相似，都是用迴圈進行連續減法來完成取餘數的操作，但一個使用了 __ffs、另一個則是使用位元操作尋找最低位的 1（Find First Set）：

#if !defined(CONFIG_CPU_NO_EFFICIENT_FFS)

/* If __ffs is available, the even/odd algorithm benchmarks slower. */

// 使用 __ffs 的實作

#else

// 使用位元操作的實作
    
#endif

而註解則說明了使用兩種方式實作 GCD 的原因是 __ffs 會比第二種方式快，但不是所有硬體都有支援 __ffs 的指令，因此這邊用 CONFIG_CPU_NO_EFFICIENT_FFS 來檢查要使用哪種實作方式。

使用 `__ffs` 實作























unsigned long gcd(unsigned long a, unsigned long b)
{
	unsigned long r = a | b;

	if (!a || !b)
		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;
	}
}

與上面使用 __builtin_ctz 改寫的部份相似，這邊的實作的概念也是先移除 a、b 尾端連續的 0（

a = x \times 2^{n}, b = y \times 2^{m}

），最後（第 17 行）再透過左移補回去。但不同的是這個實作方式是先將 a、b 進行 OR 操作，所以 r 尾端連續 0 的數量就等同於

m i n (n, m)

，不需要額外進行像前面實作中的 min(utz,vtz)。

另一個不同點則是當

x

與

y

的最大公因數為 0 時的處理方式，在使用 __builtin_ctz 的實作中可以直接用統一的形式 u << min(utz, ctz) 回傳答案，但這邊的實作則是用 r & -r 回傳

2^{m i n (n, m)}

。

`r & -r` 的技巧

假設有一數 r 可以用 0b10010100 表示的話，r & -r 的意義相當於將最低位的 1 以外的位元都設為 0，即 r & -r 的結果為 0b00000100。

而這個操作的原理是使用二補數的特性：

    10010100        (r)
    01101011        (~r, r 的一補數) 
    01101100        (-r, r 的二補數 = ~r + 1)

由於 r 的二補數相當於 r 的一補數再加 1，而當 + 1 這個會造成進位並導致：

最低位的 0 被設為 1 (假設為第 k 個位元)
最低位的 0 右側的 1 (第 0 ~ k-1 位元) 則被設為 0

因此可以將 r 與 -r 拆成三個部份來看：

第 k+1 以上的位元：-r 這部份相當於 ~r，所以與 r AND 後皆為 0
第 k 個位元：r 與 -r 的第 k 個位元都是 1，AND 後也是 1
第 0 ~ k-1 個位元：r 的這部份與 -r 都是 0，AND 後也是 0

因此 r & -r 就相當於將第 k 個位元（最低位的 1）以外的位元設為 0，即相當於 1 << ctz(r)。

使用位元操作實作
































unsigned long gcd(unsigned long a, unsigned long b)
{
	unsigned long r = a | b;

	if (!a || !b)
		return r;

	/* Isolate lsbit of r */
	r &= -r;

	while (!(b & r))
		b >>= 1;
	if (b == r)
		return r;

	for (;;) {
		while (!(a & r))
			a >>= 1;
		if (a == r)
			return r;
		if (a == b)
			return a;

		if (a < b)
			swap(a, b);
		a -= b;
		a >>= 1;
		if (a & r)
			a += b;
		a >>= 1;
	}
}

另一個實作方式也有用到前一個實作方式中的 r & -r 的技巧來取得 1 << ctz(r)，但是並不會將 a（

x \times 2^{n}

）、b （

y \times 2^{m}

）尾端的

m i n (n, m)

個 0 移除，反而是只留下尾端的

m i n (n, m)

個 0，第 11 及 17 行的 while 迴圈即是用來將多餘的 0 移除。

因此當

G C D (x, y) = 1

時，回傳 r 就代表

2^{m i n (n, m)}

，而

G C D (x, y) \neq 1

時則可以直接回傳 a（或 b），而不須進行額外的左移。

另一個特別之處則是在第 27 ~ 30 行的部份。在這個實作中只有第 12 行以及第 25 行兩處會修改到 b 的值：

第 12 行的部份會移除尾端多餘的 0，也就是保證 b >> ctz(r) 必為奇數
第 25 行的 swap 則因為 a 在第 18 行時就會移除尾端多餘的 0，因此即使 a、b 互換 b >> ctz(r) 也必為奇數

因此 b >> ctz(r) 必為奇數，而相似的，a 在經過第 18 行與第 25 行後也會保證 a >> ctz(r) 為奇數，因此第 26 行的 a -= b 後 a 尾端連續 0 的數量必大於 ctz(r)，因此第 27 行可以右移移除一個多餘的 0。

但在移除一個 0 後，a >> ctz(r) 則不能保證必為奇數或是偶數，因此第 28 及 29 行則在 a >> ctz(r) 為奇數時用另一個奇數 b 補回，讓 a >> ctz(r) 變成偶數，因此第 30 行可以保證 a >> ctz(r) 必為偶數、可以進行右移。

第 27 ~ 30 行如果只是要移除 a 尾端多餘的 0 的話大可交給下一輪的第 18 行處理，不太明白這四行的用途，但透過 blame 找到這個函式的相關 commit 訊息，似乎是有進行實驗，待讀。

Linux 核心中 GCD 的應用場景

似乎與 FLL 相關，待補。

第四題 - Bitmap













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

這個函式可以用來找到由 bitmapsize 個 64 位元無號整數構成的陣列 bitmap 中所有的 1，並將這些 1 「相較於 bitmap 開頭的位置」依序（從低位址開始）寫到 out 陣列中（第 k 個 1 的位置會寫到 out[k] 中，即第 9 行的作用）。

而為了依序找到 bitmap 中的 1，這邊是用迴圈重複 64 次右移以及 AND 1 檢查是否為 1，若 AND 結果為 1 就紀錄該位元的位置（第 p 個數的第 i 個位元）。

由於每個數都要進行 64 次右移來檢查，所以當該數 1 的位元的數量少時會很沒效率，因此可以改用 ctz 直接找到 1 的位置：

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 找到最低位的 1 的位置，接著再透過前一題提到的 r & -r 的技巧（EXP6）將最低位的 1 以外的位元設成 0，再與原數 XOR 就能將最低位的 1 設為 0 而其他位元保持不便，因此下一個迭代就能再用 ctz 找到新的最低位的 1。

實際使用案例

檢驗以 ctz/clz 改寫後的效能改進

改進空間

Linux 核心使用 bitmap 的案例

Linux Kernel 的 drivers/md/ 目錄中包含許多 Multiple Device driver 的操作程式碼，而其中有不少操作使用了 bitmap 進行處理。

第五題 - LeetCode 166. Fraction to Recurring Decimal

說明

在這個實作中先依序檢查了分母以及分子是否為 0，接著再紀錄商的正負號並將分母以及分子轉成正數（在 32 位元有號整數中，-2147483648 沒有對應的正數，所以改用 64 位元的有號整數來避免轉成正數時出現 Integer Overflow）。






























char *fractionToDecimal(int numerator, int denominator)
{
    int size = 1024;
    char *result = malloc(size);
    char *p = result;

    if (denominator == 0) {
        result[0] = '\0';
        return result;
    }

    if (numerator == 0) {
        result[0] = '0';
        result[1] = '\0';
        return result;
    }

    // 避免 -(-2147483648) Overflow
    long long n = numerator;
    long long d = denominator;

    // 轉成正數
    if (n < 0)
        n = -n;
    if (d < 0)
        d = -d;

    bool sign = (float) numerator / denominator >= 0;
    if (!sign)
        *p++ = '-';

前置處理完成後就進入計算答案的部份，這邊將除法的結果分成四個部份來看：

整數部份的商
小數點
不循環小數部份的商
循環小數部份的商

其中除了「整數部份的商」之外都可能會不存在，因此在第 35 行的部份若發現能夠整除就直接回傳整數部份的商；而若是無法整除則須補上小數點。









    long long remainder = n % d;
    long long division = n / d;

    sprintf(p, "%ld", division > 0 ? (long) division : (long) -division);
    if (remainder == 0)
        return result;

    p = result + strlen(result);
    *p++ = '.';

接著是小數部份的商，但由於小數部份可能會出現循環小數，而循環部份要寫在 () 內，因此需要另外使用 decimal 紀錄小數部份的商，並使用了 Linux Kernel 的 list API 來實作雜湊表，共有 1333 個欄位，並以餘數為 key、key % size 為雜湊函式紀錄此餘數對應的商在 decimal 中的第幾個字元：


































    char *decimal = malloc(size);
    memset(decimal, 0, size);
    char *q = decimal;

    size = 1333;
    struct list_head *heads = malloc(size * sizeof(*heads));
    for (int i = 0; i < size; i++)
        INIT_LIST_HEAD(&heads[i]);

    for (int i = 0; remainder; i++) {
        int pos = find(heads, size, remainder);
        if (pos >= 0) {
            while (PPP > 0)
                *p++ = *decimal++;
            *p++ = '(';
            while (*decimal != '\0')
                *p++ = *decimal++;
            *p++ = ')';
            *p = '\0';
            return result;
        }
        struct rem_node *node = malloc(sizeof(*node));
        node->key = remainder;
        node->index = i;

        MMM(&node->link, EEE);

        *q++ = (remainder * 10) / d + '0';
        remainder = (remainder * 10) % d;
    }

    strcpy(p, decimal);
    return result;
}

餘數不為 0 時，則須從雜湊表中檢查這個餘數是否重複出現（循環）：
- 沒出現過的餘數：將商加入到 decimal 尾端，並將這個餘數以及對應的商在 decimal 的位置（假設為 pos）加到雜湊表中（第 61 ~ 65 行），因此 MMM 及 EEE 應分別為 list_add 以及 &heads[remainder % size]。
- 重複出現的餘數（假設目前位置為 i），可分成兩個部份：
  - 不循環小數：透過第 52 行的迴圈依序將非循環部份 decimal[0:k] 加入到儲存答案的陣列 p，因此 PPP 應為 pos-- 才能使迴圈結束時的 decimal 指標指向第 pos 個字元。
  - 循環小數：依序加入 (、循環部份 decimal[k:i]、) 加入到儲存答案的字元陣列 p，最後再直接回傳答案。
餘數為 0 時，結束迴圈並將整個 decimal 複製到小數點後（無循環小數）。

原題目中有提到：

It is guaranteed that the length of the answer string is less than 104 for all the given inputs.

但不知道為什麼這個實作中 result 這個儲存答案的字元陣列的大小為 1024。

也不明白為什麼雜湊表的大小為 1333。

尋找並實作可改善的部份

雜湊表欄位數、答案長度可用巨集定義
分子分母可用 bitwise 進行絕對值運算後再轉型成 32 位元無號整數
- 但餘數必須為 64 位元無號整數，否則除以大數時可能會 Integer Overflow
無循環時 pos 會是 -1，可用一補數檢查是否為 0
雜湊表大小可以直接宣告成陣列
strcpy
- 應使用比較安全的 strncpy
- decimal 大小應考慮整數部份以及小數點
  - decimal 大小應小於 result 大小以免在 strcpy 時超出範圍
decimal 以及 heads 應該要在程式結束前釋放記憶體
- 也因此在重新實作的第 96 行的 if 內不可以直接對 fraction 操作





































































































































#include "list.h"

#define MAX_LEN 10000
#define TABLE_SIZE (1 << 10)

#define mod(dvd, dvs) ((dvd) % (dvs))
#define div(dvd, dvs) ((dvd) / (dvs))
#define HASH(key) mod(key, TABLE_SIZE)

struct rem_node {
    int key;
    int index;
    struct list_head link;
};

static inline int find(struct list_head *heads, int key)
{
    struct rem_node *entry;
    list_for_each_entry (entry, &heads[HASH(key)], link)
        if (entry->key == key)
            return entry->index;
    return -1;
}

static inline void add(struct list_head *heads, int key, int pos)
{
    // new node
    struct rem_node *node = malloc(sizeof(*node));
    node->key = key;
    node->index = pos;

    // add node
    list_add(&node->link, &heads[HASH(key)]);
}

static inline void delete_table(struct list_head *heads)
{
    for (int i = 0; i < TABLE_SIZE; ++i) {
        struct rem_node *entry, *next;
        list_for_each_entry_safe (entry, next, &heads[i], link) {
            list_del(&entry->link);
            free(entry);
        }
    }
}

char *fractionToDecimal(int numerator, int denominator)
{
    // init result string
    char *result = calloc(1, MAX_LEN + 1), *ptr = result;

    // 分母不應為 0
    if (!denominator)
        return result;

    // 分子也不應為 0（必須先檢查，不然可能會被加負號）
    if (!numerator) {
        *ptr++ = '0';
        return result;
    }


    // 檢查答案正負以及分子分母取絕對值
    int n_mask = (numerator >> 31), d_mask = (denominator >> 31);
    uint32_t n = (uint32_t)(n_mask ^ numerator) - n_mask;
    uint32_t d = (uint32_t)(d_mask ^ denominator) - d_mask;

    if (n_mask ^ d_mask)
        *ptr++ = '-';


    // 整數部份
    uint64_t rem = mod(n, d), quo = div(n, d);
    sprintf(ptr, "%lu", quo);

    // 是否整除
    if (!rem)
        return result;

    // 無法整除需要加小數點
    ptr = &result[strlen(result)];
    *ptr++ = '.';


    // 小數部份需要另外存以及雜湊表來處理循環部份
    int frac_len = MAX_LEN - (ptr - result) + 1;
    char *fraction = calloc(1, frac_len), *fptr = fraction;

    struct list_head heads[TABLE_SIZE];
    for (int i = 0; i < TABLE_SIZE; ++i)
        INIT_LIST_HEAD(&heads[i]);

    for (int i = 0; rem; ++i) {
        // 檢查是否循環（找得到重複餘數）
        int pos = find(heads, rem);
        if (~pos) {
            // 不循環部份
            fptr = fraction;
            while (pos-- > 0)
                *ptr++ = *fptr++;

            // 循環部份
            *ptr++ = '(';
            while (*fptr)
                *ptr++ = *fptr++;
            *ptr++ = ')';

            goto release;
        }

        // 無循環，新增到雜湊表
        add(heads, rem, i);

        // 新增商到 fraction 並重新檢查餘數
        *fptr++ = div(rem * 10, d) + '0';
        rem = mod(rem * 10, d);
    }
    strncpy(ptr, fraction, frac_len);

release:
    delete_table(heads);
    free(fraction);
    return result;
}

int main(int argc, char const *argv[])
{
    char *res = fractionToDecimal(atoi(argv[1]), atoi(argv[2]));
    printf("%s", res);
    free(res);
    return 0;
}

LeetCode 的 runtime 測量有點神奇，相同的程式碼第一次提交跑 9ms，第二次提交卻 0ms。

尋找並解說 Linux Kernel 的記憶體管理中的使用案例

第六題 - 實作 `alignof`

/*
 * ALIGNOF - get the alignment of a type
 * @t: the type to test
 *
 * This returns a safe alignment for the given type.
 */
#define ALIGNOF(t) \
    ((char *)(&((struct { char c; t _h; } *)0)->M) - (char *)X)

巨集實作說明

在這個巨集中，為了找到型別 t 對齊時需要的位元組大小，使用了與 offsetof 相似的技巧，並將 0 看作是一個有兩個成員的結構體的開頭地址：

第一個成員的型別需要是僅有一位元組大小的型別，才能確保型別 t 對齊需要的大小不會小於第一個成員對齊需要的大小
第二個成員的型別則是 t，若此型別對齊需要的大小比第一個成員對齊需要的大小還大的話，這個結構體實際佔用的空間就會是型別 t 對齊所需的大小的兩倍

接著由於這個結構體物件的起始地址是 0，所以第二個成員的地址（即 _h 相對於此結構體開頭的 offset）就是這個結構體佔用空間的一半，也就是型別 t 對齊所需的大小，因此 M 應為 _h

TODO : struct padding, null pointer dereference

到目前為止，感覺上 _h 的 offset 感覺就是 __align__ 要找的數值了，為什麼還需要轉型成 char * 再減去另一個指標？

首先可以先從規格書 n1570 中的 6.5.6 Additive operators 找到關於減法運算的情境：

For subtraction, one of the following shall hold:
— both operands have arithmetic type;
— both operands are pointers to qualified or unqualified versions of compatible complete object types; or
— the left operand is a pointer to a complete object type and the right operand has integer type.

可以看到減法運算要求：

兩個運算元都是算術型別
兩個運算元為 compatible qualified 或 unqualified 型別的指標
左邊的運算元是 complete object type 的指標、右邊的運算元則是整數型別

而在這個巨集的實作是兩個運算元都是指標，所以應該屬於第二種情境。為了確定它是屬於第二種情境，可以先從規格書中的 6.2.5 Types 看到關於 qualified 及 unqualified 型別的說明：

Any type so far mentioned is an unqualified type . Each unqualified type has several qualified versions of its type,
$^{47)}$ corresponding to the combinations of one, two, or all three of the const, volatile, and restrict qualifiers. The qualified or unqualified versions of a type are distinct types that belong to the same type category and have the same representation and alignment requirements.
$^{48)}$ A derived type is not qualified by the qualifiers (if any) of the type from which it is derived.

這段說明了在這（第 26 項）以前提到的型別都是屬於 unqualified type，而在這之前提到了 char, integer types, floating types 等 basic types 以及 function, pointer 等 derived types，可以看到 char 是屬於 unqualified type。

接著則是 compatible type 的說明，可以在規格書的 6.2.7 Compatible type and composite type 看到相關說明：

Two types have compatible type if their types are the same. (後略)

在第一項的一開始就說明了，若兩個型別相同的話，他們就是 compatible type，因此可以確定這個巨集的減法運算確實是屬於第二種情境（兩個運算元是 compatible unqualified type 的指標）。

接著再回到 6.5.6 Additive operators 找到關於兩指標相減的說明：

For the purposes of these operators, a pointer to an object that is not an element of an array behaves the same as a pointer to the first element of an array of length one with the type of the object as its element type.

首先，在這邊的說明可以知道一個指向非陣列元素的指標等同於指向一個大小為 1 的陣列的首位元素，接著再看兩個指向相同陣列中的元素的指標減法運算的說明：

When two pointers are subtracted, both shall point to elements of the same array object, or one past the last element of the array object; the result is the difference of the subscripts of the two array elements. The size of the result is implementation-defined, and its type (a signed integer type) is ptrdiff_t defined in the <stddef.h> header. If the result is not representable in an object of that type, the behavior is undefined. In other words, if the expressions P and Q point to, respectively, the i-th and j-th elements of an array object, the expression (P)-(Q) has the value i−j provided the value fits in an

從這邊可以看到當兩個指向相同陣列中的元素的指標進行減法運算時，結果會是一個有號整數，而其數值則等同於兩指標指向的 index 的差值。

在看完指標間減法運算的特性後，推測在找到 _h 的 offset 後還另外進行減法運算的用意是為了將 &((struct { char c; t _h; } *)0)->_h 這個結構體中 _h 地址透過有定義的方式轉型成一個數值，因此 X 應為 0 才不會影響到正確地址。

為什麼是轉型成 char * 進行減法運算？

轉形成 long * 再進行減法運算會出現不符合預期的結果？

第 9 點有說到兩指標應該要指向「相同陣列物件」，但實作中卻只是將 0 與 _h 的地址轉型成 char * 而已

Linux 核心原始程式碼中找出並解說 `alignof` 的 2 個使用情境

https://github.com/torvalds/linux/blob/master/arch/s390/mm/vmem.c

找出 ALIGN, ALIGN_DOWN, ALIGN_UP 等巨集，探討其實作機制和用途

// include/uapi/linux/const.h
#define __ALIGN_KERNEL_MASK(x, mask) (((x) + (mask)) & ~(mask))
#define __ALIGN_KERNEL(x, a) __ALIGN_KERNEL_MASK(x, (typeof(x)) (a) -1)

`ALIGN` - include/linux/align.h

#define ALIGN(x, a)		__ALIGN_KERNEL((x), (a))

`ALIGN_DOWN` - include/linux/align.h

#define ALIGN_DOWN(x, a)	__ALIGN_KERNEL((x) - ((a) - 1), (a))

`ALIGN_UP` - tools/testing/selftests/net/tcp_mmap.c

#define ALIGN_UP(x, align_to) (((x) + ((align_to) -1)) & ~((align_to) -1))

TODO

第七題 - FizzBuzz

int main() {
    for (unsigned int i = 1; i < 100; i++) {
        if (i % 3 == 0) printf("Fizz");
        if (i % 5 == 0) printf("Buzz");
        if (i % 15 == 0) printf("FizzBuzz");
        if ((i % 3) && (i % 5)) printf("%u", i);
        printf("\n");
    }
    return 0;
}

上面的程式碼是連續的 if 而不是 if...else...，所以在 15 的倍數時前三個 if 都會成立，並印出 FizzBuzzFizzBuzz，不確定是題目刻意設計的還是寫錯。
























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

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

從迴圈最後看回去會比較好理解。

首先從最後的 strncpy 可以看出是要透過 "FizzBuzz%u" 這個字串來輸出，而根據題目可知：

div5 && div3 時應該要印出 FizzBuzz
div3 == true 時應該要印出 Fizz
div5 == true 時應該要印出 Buzz
其他情況則要印出數字

所以在上述四種情況時應該分別從 FizzBuzz%u 字串的特定位置開始印出特定長度的字串內容：

從第 0 個字元開始印出 8 個字元
從第 0 個字元開始印出 4 個字元
從第 4 個字元開始印出 4 個字元
從第 8 個字元開始印出 2 個字元

第 17 行應該是 8 >> div5，不然非 3, 5 倍數時只會印出 u。

接著就可以透過這個條件知道 KK1, KK2, KK3 應該填入什麼：

KK3 應為 div3 << 2 才能讓 div3 == true 時從 &"FizzBuzz"[0] 開始印出。

KK1, KK2 要讓 length 在：

div3 == true 時為 4 (2 << 1)
div5 == true 時為 4 (2 << 1)
div3 && div5 時為 8 (2 << 2) == (2 << 1) << 1
!div3 && !div5 時為 2 (2 << 0)

因此 KK1, KK2 應分別為 div3 以及 div5。

評估 `naive.c` 和 `bitwise.c` 效能落差

TODO

設計另一種 bitmask 試圖運用更少的指令來實作出 branchless

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

https://github.com/lemire/fastmod

throughput (吞吐量) 更高的 Fizzbuzz 實作

解析 Linux 核心原始程式碼 kernel/time/timekeeping.c 裡頭涉及到除法運算的機制

相似題 - LeetCode 412. FizzBuzz

char **fizzBuzz(int n, int *returnSize)
{
    *returnSize = n;
    char **res = (char **) malloc(n * sizeof(char *));

    for (int i = 0, num = 1; i < n; ++i, ++num) {
        bool div3 = is_divisible(num, M3), 
             div5 = is_divisible(num, M5);

        int length = (2 << (div3 + div5));

        char fmt[9] = {0};
        strncpy(fmt, &"FizzBuzz%d"[(8 >> div5) >> (div3 << 2)], length);

        res[i] = calloc(10, sizeof(char));
        sprintf(res[i], fmt, num);
    }
    return res;
}

大致上跟 quiz 的解法差不多，但是由於題目要求要用 malloc 回傳字串陣列，所以用了 sprintf 輸出格式化後的結果到回傳陣列。