2022q1 Homework2 (quiz2)

contributed by < Kevin-Shih >

測驗 1 對二個無號整數取平均

EXP1:

考慮以下對二個無號整數取平均值的程式碼:

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

考慮到整數溢位的因素，以下為第一種可行的實做:

#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 >> 1), (b >> 1) 先將 a、b 分別除以 2 再相加，避免了整數溢位的可能性，然而當 a、b 都是奇數時最終的結果會因為兩者在除 2 時各少了 0.5 (右移時均丟失最後 1 bit) 而比正確的平均值少 1 ，因此需要透過加上 EXP1 來補償。

有了這些瞭解就可以得出 EXP1 只有在 a, b 都是奇數，也就是 a, b 最低位均為 1 時為 1 ，其餘為 0 。因此 EXP1 = a & b & 1，其中 & 1 用來取最後一個 bit 的結果。

EXP2, EXP3:

再次改寫為以下等價的實作:

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

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

看到 return 右半的 (EXP3) >> 1 可以把 EXP3 除 2，可以看出 EXP3 主要是做 a + b 的動作，又因題目限用 OR, AND, XOR 這三個運算子，其中只有 XOR 能做到半加器的功能，因此 EXP3 = a ^ b。 接著因半加器缺少進位的動作，所以應該透過加上 EXP2 來補償，因進位只發生在 a, b 均為 1 的 bit 上，故應對 a, b 做 AND 取得需進位的 bit，且最終求的結果是

$(a+b)/2$ 因此進位留在原本的 bit 就好不須左移 1 位 (與除 2 時的右移 1 位抵銷)，依此推論得到 EXP2 = a & b。

比較上述實作在編譯器最佳化開啟時的差異

gcc 版本：

$ gcc --version
gcc (Ubuntu 9.4.0-1ubuntu1~20.04) 9.4.0

gcc -S 所得到的 assembly 預設是 AT&T syntax

• (a + b) / 2
  以 gcc -S -O3 q1_1.c 得到第一種實作經編譯器最佳化後的組合語言輸出如下:

  ​​​​// rdi: @a, rsi: @b
  ​​​​average:
  ​​​​.LFB0:
  ​​​​    .cfi_startproc
  ​​​​    endbr64
  ​​​​    leal	(%rdi,%rsi), %eax // eax = rdi + rsi
  ​​​​    shrl	%eax // eax >>= 1
  ​​​​    ret
  ​​​​    .cfi_endproc
  

  第一種實作得到的結果，很直覺的將 rdi, rsi (a, b) 中的值相加放入 eax 後，將 eax 右移 1 位並返回該結果。 從 rdi, rsi 來看是以 64-bit 來做加法再取 lower 32 bits 存入 eax。

• low + (high - low) / 2
  以 gcc -S -O3 q1_2.c 得到第二種實作經編譯器最佳化後的組合語言輸出如下:

  ​​​​// edi: @low, esi: @high
  ​​​​average:
  ​​​​.LFB0:
  ​    .cfi_startproc
  ​    endbr64
  ​    movl	%esi, %eax    // eax = esi
  ​​​​    subl	%edi, %eax    // eax -= edi
  ​​​​    shrl	%eax          // eax >>= 1
  ​​​​    addl	%edi, %eax    // eax += edi
  ​​​​    ret
  ​​​​    .cfi_endproc
  

  相較於第一種實做，全程採用 32-bit 方式存取暫存器，不像第一種以 64-bit 來做加法再取 lower 32 bits 可能導致整數溢位的問題。

• (a >> 1) + (b >> 1) + (a & b & 1)
  以 gcc -S -O3 q1_3.c 得到第三種實作經編譯器最佳化後的組合語言輸出如下:

  ​​​​// edi: @a, esi: @b
  ​​​​average:
  ​​​​.LFB0:
  ​​​​    .cfi_startproc
  ​​​​    endbr64
  ​​​​    movl	%edi, %eax // eax = edi
  ​​​​    movl	%esi, %edx // edx = esi
  ​​​​    andl	%esi, %edi // edi &= esi (a = a & b)
  ​​​​    shrl	%eax       // eax >>= 1
  ​​​​    shrl	%edx       // edx >>= 1
  ​​​​    andl	$1, %edi   // edi &=1
  ​​​​    addl	%edx, %eax // eax += edx
  ​​​​    addl	%edi, %eax // eax += edi
  ​​​​    ret
  ​​​​    .cfi_endproc
  

  基本上是以 eax, edx, edi 分別來做 (a >> 1), (b >> 1), (a & b & 1) ，但有個疑問是為何 #7 要將 esi 移到 edx ， #12 直接 addl %esi, %eax 應會得到相同的結果。

  TODO
  釐清上述疑問

  雖改用 Shift, AND 操作但最佳化後的組合語言輸出反而比前一種多了 4 行，似乎更耗時了(多了 4 個指令)，有些出乎意料。

• (a & b) + ((a ^ b) >> 1)
  以 gcc -S -O3 q1_4.c 得到第四種實作經編譯器最佳化後的組合語言輸出如下:

  ​​​​// edi: @a, esi: @b
  ​​​​average:
  ​​​​.LFB0:
  ​​​​    .cfi_startproc
  ​​​​    endbr64
  ​​​​    movl	%edi, %eax // eax = edi
  ​​​​    andl	%esi, %edi // edi &= esi
  ​​​​    xorl	%esi, %eax // eax ^= esi
  ​​​​    shrl	%eax       // eax >>= 1
  ​​​​    addl	%edi, %eax // eax += edi
  ​​​​    ret
  ​​​​    .cfi_endproc
  

  以 eax, edi分別來做 a ^ b, (a ^ b) >> 1 ，最後在相加至 eax 並返回該值。

Linux 核心原始程式碼 include/linux/average.h 的 EWMA 實作及 Linux 核心中的應用

EWMA 即 Exponentially weighted moving average ，計算移動平均時較舊的觀測量(數值)權值以指數形式遞減，但永遠不會到 0 ，其函數定義如下：

$S_t=\{\begin{array}{ll}{Y}_0,& t=0\\ \alpha Y_t+(1-\alpha )\cdot S_{t-1},& t>0\end{array}$

• $Y_t$ is the value at a time period
  $t$.
• $S_t$ is the value of the EWMA at any time period
  $t$. – Wiki: EWMA

接著看到 include/linux/average.h 中定義的 macro
#define DECLARE_EWMA(name, _precision, _weight_rcp)， 其三個參數的用途如下

• name: 定義一個對應的 struct ，名為 ewma_##name## ，以及對應的名稱的 init, read, add 函式
• _precision： 計算時 fraction part 要使用多少 bit ，不能大於 30
• _weight_rcp: 相當於公式中的
  $1/\alpha$ ，為了效率考量須為 2 的冪次(可以用 shift 搞定)

DECLARE_EWMA 所定義的 struct

struct ewma_##name {
	unsigned long internal;
};

其實只有一個 unsigned long 變數 internal 用於儲存 EWMA 的結果。

接著看到第一個函式 ewma_##name##_init

static inline void ewma_##name##_init(struct ewma_##name *e)
{
	BUILD_BUG_ON(!__builtin_constant_p(_precision));
	BUILD_BUG_ON(!__builtin_constant_p(_weight_rcp));
    /*
	* Even if you want to feed it just 0/1 you should have
	* some bits for the non-fractional part...
	*/
	BUILD_BUG_ON((_precision) > 30);
	BUILD_BUG_ON_NOT_POWER_OF_2(_weight_rcp);
	e->internal = 0;
}

前 4 行的 BUILD_BUG_ON 均是當輸入的參數不符合要求時可以在 compile 階段就報錯，以保證執行階段時的效率與正確性。 最後就是將 internal 初始化為 0 ，計算 EWMA 前應先呼叫此函式以確保 internal 初始化為 0 。

__builtin_constant_p(x) 用於檢驗是 x 否在 compile time 就已知為常數，若否則中斷 compile 並報錯。 – gcc onlinedocs

BUILD_BUG_ON(condition) break compile if a condition is true. – include/linux/build_bug.h

第二個函式 ewma_##name##_read

static inline unsigned long ewma_##name##_read(struct ewma_##name *e)
{
	BUILD_BUG_ON(!__builtin_constant_p(_precision));
	BUILD_BUG_ON(!__builtin_constant_p(_weight_rcp));
	BUILD_BUG_ON((_precision) > 30);
	BUILD_BUG_ON_NOT_POWER_OF_2(_weight_rcp);
	return e->internal >> (_precision);
}

再度出現同樣 4 行的 BUILD_BUG_ON 有點怪，似乎寫一遍即可。

TODO
find out why

讀出 e->internal 並右移 _precision 位，已得到整數位的 EWMA (因 internal 保留後 _precision bit 作為運算時的 fraction part )。

第三個函式 ewma_##name##_add

static inline void ewma_##name##_add(struct ewma_##name *e,
                                     unsigned long val)
{
	unsigned long internal = READ_ONCE(e->internal);
	unsigned long weight_rcp = ilog2(_weight_rcp);	
	unsigned long precision = _precision;

    BUILD_BUG_ON(!__builtin_constant_p(_precision));
	BUILD_BUG_ON(!__builtin_constant_p(_weight_rcp));
	BUILD_BUG_ON((_precision) > 30);
	BUILD_BUG_ON_NOT_POWER_OF_2(_weight_rcp);

	WRITE_ONCE(e->internal, internal ?
		(((internal << weight_rcp) - internal) +
			(val << precision)) >> weight_rcp :
		(val << precision));
}

計算 EWMA 的主體，用於計算新增一筆資料 val 後的新 EWMA 並存入 e->internal 中。

• #4, #13 使用了 READ_ONCE, WRITE_ONCE 應該是為了保護在多個 thread 的讀寫操做產生 data race 。
• #5 以 ilog2(_weight_rcp) 求得稍後要位移 weight_rcp 個 bit 。
• #6 先將 _precision 的值存入 precision 主要是避免 macro 展開後可能出現的 double evaluation 的情形。

#13 當 internal 為 0 時， e->internal = val << precision 符合 EWMA 定義中

$t=0$ 的情形。 而當 internal 不為 0 時，其敘述相當於

$internal=(internal(2^{weight\mathrm{\_}rcp}-1)+val)/2^{weight\mathrm{\_}rcp}$

$2^{weight\mathrm{\_}rcp}$ = _weight_rcp, 1/_weight_rcp =

$\alpha$

符合 EWMA 定義中

$t>0$ 的情形。

EWMA 在 Linux 核心中的應用
直接查找 DECLARE_EWMA 關鍵字在 kernel 中出現的地方，之前修「通訊網路」時在計算 RTT 時出現過 EWMA ，預期 Linux 核心中也會有類似的應用。

Linux 核心中的應用:

• net/mac80211/ieee80211_i.h, line 435
• net/mac80211/sta_info.h, line 126, 373, 374, 433
• 還有不少也是在 drivers/net/wireless 底下的應用…

在 sta_info.h 中 EWMA 被用於計算 TX bit rate(傳輸的 bit rate), failed MSDUs(link layer 的 payload) 的百分比等用途，。

測驗 2 不需要分支的 max

EXP4, EXP5:

考慮以下「不需分支的 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 個運算子。

題目為了取最大值，因此當 a 較大的時候回傳值應等價於 a ^ 0 = a。
而當 b 較大時回傳值應等價於 a ^ (a ^ b) = 0 ^ b = b

因此得出 EXP4 應為 a ^ b，且 -(EXP5) 在 a > b 時為 0，在 a < b 時為 -1 即 0xFFFFFFFF，其中 0xFFFFFFFF 是透過 - 對 1 取 2 的補數來的，因此 EXP5 在 a > b 時應為 0，在 a < b 時應為 1，故 EXP5 = a < b。

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

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

將輸入的 parameters 改為 int32_t 型態就能以相同的方式的到正確的 max 結果了。

在 Linux 核心原始程式碼中，找出更多類似的實作手法

arch/x86/kvm/mmu.h 中：

/*
 * If CPL < 3, SMAP prevention are disabled if EFLAGS.AC = 1.
 *
 * If CPL = 3, SMAP applies to all supervisor-mode data accesses
 * (these are implicit supervisor accesses) regardless of the value
 * of EFLAGS.AC.
 *
 * This computes (cpl < 3) && (rflags & X86_EFLAGS_AC), leaving
 * the result in X86_EFLAGS_AC. We then insert it in place of
 * the PFERR_RSVD_MASK bit; this bit will always be zero in pfec,
 * but it will be one in index if SMAP checks are being overridden.
 * It is important to keep this branchless.
 */
unsigned long smap = (cpl - 3) & (rflags & X86_EFLAGS_AC);
int index = (pfec >> 1) +
	    (smap >> (X86_EFLAGS_AC_BIT - PFERR_RSVD_BIT + 1));

git log 中還能找到一些其他的實做，但有些似乎已經消失了，在目前 github 上的 repo 找不到了

• Replace int2float() with an optimized version.
• Optimize pte permission checks

測驗 3 64 位元 GCD 求值函式

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

#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 演算法，當 u, v 任一者為 0 時， gcd 即為非零者的值，當均為 0 時，gcd 為 0 。 這三種情形透過 #4 的 if 敘述來處理。

接著 #6 的 for 迴圈會在 u, v 均為偶數時，將其同除以 2 ，並將除的次數紀錄在 shift 以便最後乘回 gcd 中(左移 shift 次)。 #9 的 while 迴圈則是當 u 還是偶數時，除 2 直到 u 變為奇數。

接著 #11 開始的 do while 就是輾轉相除法的本體了，當 u 小於 v 就將 v 減去 u ，大於的話就將 u - v 存到 v ， v 存到 u ，直到其中一者為 0 為止，然而由於這個實做中，相減後的結果都會存到 v ，因此中止條件 COND 只須判斷 v 是否為 0 即可，故 COND = v 。

最終的回傳值 RET 則是將先前輾轉相除法最後留下的 u 值左移 shift 位(先前同除以 2 的次數)得到答案，因此 RET = u << shift 。

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

透過 __builtin_ctz 改寫的 gcd64_ctz

uint64_t gcd64_ctz(uint64_t u, uint64_t v)
{
    if (!u || !v) return u | v;
    int shift = __builtin_ctz(u | 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);
    return u << shift;
}

將原先同除以 2 的部份改為紀錄 u, v 中最低位的 1 作為 shift 。

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

取代為 u 直接右移到無法被 2 整除。 v 的部份則留到迴圈內處理。

其餘部份與原先的 gcd64 相同。

效能測試方式

先以 random 方式準備測試資料

void prepare_data(uint64_t *data1, uint64_t *data2, size_t size){
    srand(time(NULL));
    for (size_t i = 0; i < size; i++) {
        data1[i] = rand();
        data2[i] = rand();
    }
}

兩者分別測試

prepare_data(data1, data2, size);
clock_t start = clock();
for (size_t i = 0; i < size; i++)
    gcd64(data1[i], data2[i]);
clock_t finish = clock();

clock_t start2 = clock();
for (size_t i = 0; i < size; i++)
    gcd64_ctz(data1[i], data2[i]);
clock_t finish2 = clock();

double time1 = (double)(finish - start) / CLOCKS_PER_SEC;
double time2 = (double)(finish2 - start2) / CLOCKS_PER_SEC;

在本機連續執行 10,000,000 次、 3 round 的結果如下
(5 round 取中間 3 round)

	gcd64	gcd64_ctz
1	2.125378 secs	0.936948 secs
2	2.126757 secs	0.936914 secs
3	2.121405 secs	0.935645 secs
avg.	2.124513 secs	0.936502 secs

採用 ctz 的實做 3 round 的平均快了約 2.268562 倍，顯示在測試環境中 ctz 顯然比除法運算快上不少。

請解釋 Linux 核心中內建 GCD (不只一種實作)，的實作手法並探討在核心內的應用場景。

lib/math/gcd.c 中的 GCD 實作手法

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

當 ffs 可用時會採用以上實做，當 a, b 中一者為 0 或均為 0 時，gcd 為 a | b (三種情形分別對應 b, a, 0)。

若 b 是 2 的冪次，則回傳 r & -r。 r & -r 的作法是透過 2 補數特性找出 r 中最低位的 1 並使其他 bit 為 0 (就像測驗 4 中 bitset & -bitset)，而由於 r = a | b 因此相當於找出 min(a & -a, b & -b) ，也就是兩者的最大公因數(因 b 是 2 的冪次，兩者 gcd 是小於等於 b 的 2 的冪次)。

若 b 不是 2 的冪次， b >>= __ffs(b) 也會如 gcd 算法中所述，將 b 除到剩下奇數的狀態。

for 迴圈內有 4 步：

• 將 a 除到剩下奇數的狀態。
• 若 a 是 2 的冪次，若是則 gcd 為 r & -r (同 b 的情形，代表沒有 2 的冪次以外的公因數)。
• 若剩下的 a == b 代表剩下的 a 也是兩者的公因數，將 a 左移先前同除的 2 的次數(ffs(r) = min(ffs(a), ffs(b)))即為 gcd 。 (再減一次就是輾轉相除法，其中一個減到變 0 的情形)
• 最後當上述情形都不成立時，執行輾轉相除法中相減的步驟，若 a < b 就先交換 a, b，確保 a 是被減的那個。

核心內的應用場景

• block/blk-settings.c
• drivers/media/tuners/mt2063.c
• sound/core/pcm_timer.c
• …

其中在 mt2063.c 中一個用於檢測 LO spur 是否發生的函數 IsSpurInBand() 中出現，用於計算某種 scale factor 。

/*
	 ** For each edge (d, c & f), calculate a scale, based on the gcd
	 ** of f_LO1, f_LO2 and the edge value.  Use the larger of this
	 ** gcd-based scale factor or f_Scale.
	 */
	lo_gcd = gcd(f_LO1, f_LO2);
	gd_Scale = max((u32) gcd(lo_gcd, d), f_Scale);
	hgds = gd_Scale / 2;
	gc_Scale = max((u32) gcd(lo_gcd, c), f_Scale);
	hgcs = gc_Scale / 2;
	gf_Scale = max((u32) gcd(lo_gcd, f), f_Scale);
	hgfs = gf_Scale / 2;

另外在 git log --grep 則發現不少 commit 是捨棄 local 的 gcd implementation 而改用 lib/math/gcd.c 的實作，例如剛剛提到的 sound/core/pcm_timer.c 在 commit a96053 時捨棄。

測驗 4 bitset

在影像處理中，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;
}
以 __builtin_ctzll 改寫的程式碼如下:
```c= 
#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];   // a row of bitmap per iteration
        while (bitset != 0) { // till this row has no more `set bit`
            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 的特性)。

題目提到 EXP6 應保留 bitset 中最低位的 1 而其餘均為 0 ，又提示要 -bitset 的特性，以 bitset 為

${010100}_b$ 為例， -bitset 為

${101100}_b$ ，兩者只有原先最低位的 1 的 bit 均為 1 ，因此善用 2 的補數特性， EXP6 = bitset & -bitset。

-bitset 特性:

	binary
bitset	${010100}_b$
~bitset	${101011}_b$
-bitset	${101100}_b$

-bitset = ~bitset + 1 因此只有最低位的 1 與 bitset 均為 1 。

設計實驗，檢驗 ctz/clz 改寫的程式碼相較原本的實作有多少改進？應考慮到不同的 bitmap density

準備測試資料的函式

void prepare_data(uint64_t *bitmap, size_t size, int density){
    srand(time(NULL));
    // control bitmap density with 0 ~ 4 level 0%, 25%, 50%, 75%, 100%
    int count = 16 * size * density;

    // 100% density
    if (density == 4){
        for (size_t i = 0; i < size; i++)
            bitmap[i] = 0xFFFFFFFFFFFFFFFFul;
        return;
    }

    // init bitmap
    for (size_t i = 0; i < size; i++)
        bitmap[i] = 0;
    // 0% density
    if (density == 0)
        return;

    // select `count` (pseudo) random `0` bits set them to `1`
    while (count) {
        int i = rand64() % (size * 64);
        if (bitmap[i / 64] & 1ul << (i % 64))
            continue;
        else{
            bitmap[i / 64] |= 1ul << (i % 64);
            count--;
        }
    }
}

測試方式

size_t bitmapsize = strtoul(argv[1], NULL, 10);
uint64_t *bitmap = (uint64_t *)malloc(bitmapsize * sizeof(uint64_t));
uint32_t *out_naive = (uint32_t *)malloc(bitmapsize * 64 * sizeof(uint32_t));
uint32_t *out_improved = (uint32_t *)malloc(bitmapsize * 64 * sizeof(uint32_t));

for(size_t i = 0; i < 5; i++){  // 5 density level 0%, 25%, 50%, 75%, 100%
    prepare_data(bitmap, bitmapsize, i);
    for (size_t j = 0; j < bitmapsize * 64; j++)
        out_naive[j] = out_improved[j] = 0;

    clock_t start = clock();
    naive(bitmap, bitmapsize, out_naive);
    clock_t finish = clock();

    clock_t start2 = clock();
    improved(bitmap, bitmapsize, out_improved);
    clock_t finish2 = clock();

    double time1 = (double)(finish - start) / CLOCKS_PER_SEC;
    double time2 = (double)(finish2 - start2) / CLOCKS_PER_SEC;
    printf("density: %3ld%%, naive: %f secs, improved: %f secs\n", i * 25, time1, time2);
}
free(bitmap);
free(out_naive);
free(out_improved);

測試 5 種不同的 bitmap density 的情形下 naive 與 ctz 改寫的版本的執行時間。

在我的環境中執行 bitmapsize = 1,000,000 的結果:

density	naive (secs)	improved (secs)
0%	0.129309	0.002487
25%	0.276713	0.069022
50%	0.390655	0.121114
75%	0.273711	0.174631
100%	0.144110	0.229330

在絕大多數時間 ctz 改寫的 improved 版本確實比 naive 快上不少，但有趣的是在 100% density (all set) 的情形中 naive 反而比 improved 快上一些。

兩個版本的主要差異：
naive

for (int i = 0; i < 64; i++) {
    if ((bitset >> i) & 0x1)
        out[pos++] = p + i;
}

improved

while (bitset != 0) { // till this row has no more `set bit`
    uint64_t t = EXP6;
    int r = __builtin_ctzll(bitset);
    out[pos++] = k * 64 + r;
    bitset ^= t;
}

在 100% density 兩者的內層迴圈都須 iterate 64 次且 naive 的 if 條件恆為真，因此 prediction miss 幾乎不會出現，因此迴圈內單純對 out[pos++] 賦值的 naive 反而比較快。

另外， naive 版本在 50% density 時的執行時間最久，應是因為 prediction miss 較高， improved 版本則是隨著 density 上升執行時間隨之遞增 (ctz 所能省下內層迴圈 iteration 次數越來越少)。

思考進一步的改進空間

size_t improved_v2(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) {
            uint64_t t = bitset & -bitset;
            out[pos++] = k << 6 | __builtin_ctzll(bitset);
            bitset ^= t;
        }
    }
    return pos;
}

#9 中避免使用 *, + 運算改用 <<, | ，並將原本 __builtin_ctzll(bitset) 不再存入 r 而直接移至 #9 使用讓程式更簡潔。

測試程式

int main(int argc,char **argv){
    if (argc != 2)
        return 1;

    size_t bitmapsize = strtoul(argv[1], NULL, 10);
    uint64_t *bitmap = (uint64_t *)malloc(bitmapsize * sizeof(uint64_t));
    uint32_t *out_naive = (uint32_t *)malloc(bitmapsize * 64 * sizeof(uint32_t));
    uint32_t *out_improved = (uint32_t *)malloc(bitmapsize * 64 * sizeof(uint32_t));
    uint32_t *out_improved_v2 = (uint32_t *)malloc(bitmapsize * 64 * sizeof(uint32_t));

    for(size_t i = 0; i < 5; i++){ // 5 density level 0%, 25%, 50%, 75%, 100%
        prepare_data(bitmap, bitmapsize, i);
        for (size_t j = 0; j < bitmapsize * 64; j++)
            out_naive[j] = out_improved[j] = out_improved_v2[j] = 0;

        clock_t start = clock();
        naive(bitmap, bitmapsize, out_naive);
        clock_t finish = clock();

        clock_t start2 = clock();
        improved(bitmap, bitmapsize, out_improved);
        clock_t finish2 = clock();

        clock_t start3 = clock();
        improved_v2(bitmap, bitmapsize, out_improved_v2);
        clock_t finish3 = clock();

        // test output correctness
        int equal = 1;
        for (int i = 0; i < bitmapsize * 64; i++){
            if (out_naive[i] != out_improved[i]) {
                printf("improved failed\n");
                equal = 0;
                break;
            }
            if (out_naive[i] != out_improved_v2[i]) {
                printf("improved_v2 failed\n");
                equal = 0;
                break;
            }
        }
        double time1 = (double)(finish - start) / CLOCKS_PER_SEC;
        double time2 = (double)(finish2 - start2) / CLOCKS_PER_SEC;
        double time3 = (double)(finish3 - start3) / CLOCKS_PER_SEC;
        printf("density: %3ld%%, naive: %f secs, improved: %f secs, improved_v2: %f secs\n", i * 25, time1, time2, time3);
    }
    free(bitmap);
    free(out_naive);
    free(out_improved);
    free(out_improved_v2);
    return 0;
}

(與前面測試 ctz 版本的方式相同，但加入檢查 output 正確性的程式，避免 improved_v2 改出 bug )

在我的環境中執行 bitmapsize = 1,000,000 的結果如下

density	naive (secs)	improved (secs)	improved_v2 (secs)
0%	0.135698	0.002625	0.002510
25%	0.291539	0.071203	0.067933
50%	0.395848	0.121448	0.117784
75%	0.274069	0.173835	0.169411
100%	0.143882	0.213872	0.206336

聚焦在 improved, improved_v2 的折線圖

從折線圖可以看到避免使用 *, + 而改用 bitwise operator <<, | 可以節省一小部份的執行時間，而且在 density 較高時比較明顯。

閱讀 Data Structures in the Linux Kernel 並舉出 Linux 核心使用 bitmap 的案例

• drivers/net/ipa/ipa_endpoint.h : commit c1aaa0
• net/mptcp/pm_netlink.c : commit 86e39e
• sound/usb/clock.c #260
• …

大多數是用來儲存一種狀態如 flags, visted 等等，相較使用 bool 佔用 1 byte 儲存 1 個狀態 (true / false)， bitmap 1 byte 最多儲存 8 個狀態，若是需要節省空間， bitmap 是很好的選擇。

include/linux/types.h 中 bitmap 是以 unsigned long array 方式宣告的，因此若要存 65 個狀態就需要 128 bits 似乎也是會浪費不少空間 (但同時也保留了不少未來擴充更多狀態的空間)