2024q1 Homework4 (quiz3+4)

contributed by <csotaku0926>

第三週

測驗一

• 原理

根據 digit-by-digit calculation 原理，假設欲求平方根為

$N$ ，可將

$N^2$ 化為 2 的冪總和:

$$N^2=(a_0+a_1+...+a_n{)}^2,a_m=2^m(or{2}^m=0),m\in \{0,1,...,n\}$$
已知

$(a+b{)}^2=a^2+2ab+b^2$，則有以下遞迴關係

(單純看式子很難懂，可以用平方矩陣去想比較好理解)

$$\left[\begin{array}{ccc}{a}_0{a}_0& a_0{a}_1& a_0{a}_2\\ a_1{a}_0& a_1{a}_1& a_1{a}_2\\ a_2{a}_0& a_2{a}_1& a_2{a}_2\end{array}\right]$$

$N^2=((a_0+a_1+...+a_{n-1})+a_n{)}^2=a_n^2+2a_n\sum _{k=0}^{n-1}{a}_k+(a_0+a_1+...+a_{n-1}{)}^2$

$=a_n^2+(2P_n+a_{n-1})a_{n-1}+...+(2P_1+a_0)a_0$

令平方根逼近值

$P_m=\sum _{i=m}^n{a}_i$，

$m\in \{0,1,...,n\}$

$P_0=\sum _{i=0}^n{a}_i$ 即是欲求平方根

顯然

$P_{m-1}=P_m+a_{m-1}$
因此對於每個

$m\in \{n,...,1,0\}$，於

$P_{m-1}$ 加上

$a_m$，由於是 2 的冪，因此只有以下兩種可能:

$P_{m-1}=P_m+a_{m-1},(P_m^2<N^2)$,

$P_{m-1}=P_m,(otherwise)$

但我們的程式並不是用這個較直白的方法判別是否加

$a_{m-1}$ ，而是用差值的方式
令

$X_{n+1}=N^2$ (考慮從 n 開始),

$Y_m=(2P_m+a_{m-1})a_{m-1}$
對於每個 m:

$X_m=X_{m+1}-Y_m=N^2-P_m^2$

令我費解的是，程式又進一步將

$Y_m$ 拆解成

$c_m$,

$d_m$
想了很久，原來只是單純將加法的兩項拆開來看，即:

$c_m=P_{m+1}{2}^{m+1}$,

$d_m=(a_m{)}^2$

$Y_m=c_m+d_m$

每次迴圈更新:

$\begin{array}{c}{c}_{m-1}=P_m{2}^m=(P_{m+1}+a_m)2^m=\{\begin{array}{l}{c}_m/2+d_{m}if{a}_m=2^m\\ c_m/2if{a}_m=0\end{array}\end{array}$

$d_{m-1}=d_m/4$

計算到最後

$c_{-1}=P_0=N$,

$d_0=0$

太難理解了，單單根號運算為了追求表現，竟然會變得這麼複雜 (我資質駑鈍)

• 以 ffs 改寫

commit e413059

原先 (31 - __builtin_clz(x)) 的部份是為了算 MSB，現在用 strings.h 提供的 ffs 達成

• Linux kernel 的平方根運算

要理解 linux 核心程式碼並不是容易的事情，光是說明提供參考的 block 目錄就有許多檔案，要如何在其中找到平方根實作？

透過 google 關鍵字搜尋，總算是在 block/blk-wbt.c 這個檔案中找到有關平方根的動作:

- If the minimum latency in the above window exceeds some target, increment
 *   scaling step and scale down queue depth by a factor of 2x. The monitoring
 *   window is then shrunk to 100 / sqrt(scaling step + 1).

由註解看出這是一段有關調整 monitor window size 的程式碼，有點像網路程式設計學到的滑動視窗 (sliding window) ，可以根據 scaling step 逐漸縮小

實際呼叫平方根的函式為 rwb_arm_timer :

static void rwb_arm_timer(struct rq_wb *rwb)
{
	struct rq_depth *rqd = &rwb->rq_depth;

	if (rqd->scale_step > 0) {
		/*
		 * We should speed this up, using some variant of a fast
		 * integer inverse square root calculation. Since we only do
		 * this for every window expiration, it's not a huge deal,
		 * though.
		 */
		rwb->cur_win_nsec = div_u64(rwb->win_nsec << 4,
					int_sqrt((rqd->scale_step + 1) << 8));

而定義 int_sqrt 的函式位於 /lib/math/int_sqrt.c

unsigned long int_sqrt(unsigned long x)
{
	unsigned long b, m, y = 0;

	if (x <= 1)
		return x;

	m = 1UL << (__fls(x) & ~1UL);
	while (m != 0) {
		b = y + m;
		y >>= 1;

		if (x >= b) {
			x -= b;
			y += m;
		}
		m >>= 2;
	}

	return y;
}
EXPORT_SYMBOL(int_sqrt);

除了獲取 MBS 的方式是使用 _fls, 將迴圈改為 while loop 之外，整體演算法是相同的

• 閱讀論文與改進
• 論文1
• 論文2

測驗二

• 程式解讀

根據教師所言，這題目是希望引導學員思考針對 RV32I/RV64I 這類沒有整數乘法和除法指令的指令集，該如何撰寫有效的程式

希望以 bitwise operation 實作除法，但若除數 (這個例子中, 10) 包含非 2 的冪的因數，會有不精確的問題
項目說明 中透過推導證明，若存在

$n\in N$，選定一

$x$ 使

$\frac{n}{x}$ 直到小數第一位都是精準的，即

$a.b\le \frac{n}{x}\le a.b9$
則對於任何小於 n 的非負整數 l ，

$\frac{l}{x}$ 直到小數第一位也是精準的

考慮輸入不大於 19 的狀況下，找出能夠維持小數第一位精確度的除數 x :

$1.9\le \frac{19}{x}\le 1.99$ –>

$9.55\le x\le 10$

由於是 bitwise operation，除數的形式應該要是

$\frac{a}{2^n}$，a 是任意配對的正整數

接下來在測試程式處，有個匪夷所思的地方:

q = ((((n >> 3) + (n >> 1) + n) << 3) + d0 + d1 + d2) >> 7;

這段程式碼想要求得 13 * tmp ，但他的步驟是先求得 13 * tmp / 8 後透過左移變回 13 * tmp，又此時因經過右移導致部份位元遭捨棄，因此又用 d0, d1, d2 補回。
為什麼不直接求 13 * tmp 就好了？ 可以進一步寫成以下形式，節省變數運用:

q = ((n << 3) + (n << 2) + n) >> 7;

可見測試題中的程式碼只是其中一種可行方法，不代表最優解

而包裝程式又換了一種寫法:

uint32_t x = (in | 1) - (in >> 2); /* div = in/10 ==> div = 0.75*in/8 */
uint32_t q = (x >> 4) + x;

這一行初步認為是將除數 x 換成

$0.75\times in$，而 q 的值指派為

$\frac{51\times in}{64}=0.796\times in$

接下來又是一個匪夷所思的部份:

q = (q >> 8) + x;

這行程式碼重複了四次，也就是說，in 不斷加入 in >> 8 數值，若 in 是小於 256 的正整數，這些程式對 in 毫無作用。也許要結合後段程式碼才會知道這樣寫的目的

最後，q 需要再右移一個數，才是最終的商。考慮

$1.9\le \frac{19}{x}\le 1.99$ –>

$9.55\le x\le 10$
需要將

$\frac{64}{51}$ 乘上適當的 2 的冪，使其符合 x 的條件
而

$\frac{64\times 8}{51}=10.03$，最符合上述條件，可知 *div = (q >> 3);

最後一行是在算餘式，已知

$div=q>>3$，

$mod=n-div\times 10$
q & ~0x7 相當於

$div\times 8$
所以: *mod = in - ((q & ~0x7) + (*div << 1));

• instruction cycle 計算

項目參考
一行行看指令對照表有些費時，如果可以做個程式自動計算 cycle 應該很不錯

想要計算這部份，第一步是要確認電腦的處理器版本，我使用的型號為 Intel 12th gen Core

  Model name:   12th Gen Intel(R) Core(TM) i5-12450H

對照圖表，並沒有發現直接對應的微處理器，先以型號最接近的 "Nehalem" 探討

先以 gcc div10.c -o div10 產出執行檔後，以 objdump -d div10 > div10.asm 觀察反組譯的結果 (-O 是 gcc 最佳化選項，預設是 -O0 )

預設最佳化等級取決於 GNU toolchain 編譯時的選項

jserv

首先是以除法和模數寫成的直覺版本:

void div10_naive(uint32_t n, uint32_t *q, uint32_t *r)
{
    *q = n / 10;
    *r = n % 10;
}

實際查看指令表格，有許多欄位。觀看敘述得知 "latency" 以 cpu cycle 為時間單位測量指令延遲，因此要對照這個欄位，以下是出現在反組譯檔案中的指令之對應 cycle 數:

指令	cycle	count
push r	3	1
mov r r/i	1	18
imul	3	2
add/sub	1	3
shr/shl	1	4
pop r	2	1

佔用 CPU 週期總計 36 個

至於只使用位元運算的除法:

void div10(uint32_t n, uint32_t *q, uint32_t *r)
{
    *q = ((n << 3) + (n << 2) + n) >> 7; /* n / (128/13) */
    *r = n - (((*q << 2) + *q) << 1);
}

指令	cycle	count
push r	3	1
mov r r/i	1	18
lea	1	3
add/sub	1	4
shr/shl	1	2
pop r	2	1

佔用 CPU 週期總計 32 個，較先前版本略為減少，發現 imul 被換成 lea 指令

參考 var-x 的 writeup ，可以利用 perf stat <command> 這個命令測量 cpu cycle

• 註: 第一次使用 perf stat 會遇到權限問題，需要修改 /proc/sys/kernel/perf_event_paranoid
• 可以利用 sudo sysctl <parameter>=<value> 快速修改，ref

於同個 c 檔案定義好兩個除法函式後，以 main() 呼叫之。編譯出執行檔後測試 perf stat ./div10 :

int main(void)
{
    uint32_t r, q;
    div10(123, &q, &r);
}

初始版本的除法耗費 977,489 個 cycle，而不用除法模數的版本只花費 880,667 個，這是很大的進步

• 撰寫 modulo 5, modulo 9 (不依賴除法指令)

注意先前的 div10 指令，預期輸入是不比 19 大的正整數，但現在有效範圍必須是 uint32_t 以內

Modulus without Division, a tutorial 一文針對模數處理，提出了不需除法指令的實作

commit fd9cce

考慮 modulo 5，已知

$8\mathrm{mod}5=3$，因此符合以下條件:

$a\mathrm{mod}5=(\frac{3a}{8}+a\mathrm{mod}8)\mathrm{mod}5$
可以用這個式子，加上迴圈判斷數值是否還需要做減去的動作。
這個原則可以運用於任一正整數模數，即:

$$a\mathrm{mod}m=((2^n\mathrm{mod}m)\frac{a}{2^n}+a\mathrm{mod}{2}^n)\mathrm{mod}m$$

其中 m 為任一正整數，n 為非負整數

但是這麼做的時間成本還是太高，其中一種作法是用 mod 15 推算出 mod 5

$$a\mathrm{mod}5=(a\mathrm{mod}15)\mathrm{mod}5$$

因為 mod 15 存在快速用位元運算計算的方式

commit b408910

至於 mod 9，目前想不到較為直接的方式計算
我想到的一種方法為先計算 mod 63，再推算出 mod 9

測驗三

程式碼回傳 log2 數值
顯然，只須回傳最高位數 (MSB) 是第幾位即可
對於 32 位元的無號整數，則可以透過 GNU extension 的 __builtin_clz() 達成

linux 核心 log2 案例

透過 google 搜尋，找到 bootlin 這個網站，裡面記載了許多核心的程式碼
其中可以找到位於核心目錄的 /include/linux/log2.h

針對 arch 架構，允許透過在 asm/bitops.h 撰寫比使用 fls 更有效率的 log2 實作
這也解釋為什麼許多標頭檔都會加上 #ifndef...endif 這類的巨集

#ifndef CONFIG_ARCH_HAS_ILOG2_U32
static __always_inline __attribute__((const))
int __ilog2_u32(u32 n)
{
	return fls(n) - 1;
}
#endif

可以發現具有常數 log2 的實作方式

#define const_ilog2(n)				\
(						\
	__builtin_constant_p(n) ? (		\
		(n) < 2 ? 0 :			\
		(n) & (1ULL << 63) ? 63 :	\
		(n) & (1ULL << 62) ? 62 :	\
		// ...
		(n) & (1ULL <<  2) ?  2 :	\
		1) :				\
	-1)

根據 gnu 官網，__bulitin_constant_p() 可用於檢驗表達式在編譯時間是否為常數

參考 編譯器和最佳化原理篇，在編譯器最佳化的策略上，經常可以看到將迴圈內容展開的步驟，稱作 loop unrolling

也有提供非常數的實作:

#define ilog2(n) \
( \
	__builtin_constant_p(n) ?	\
	((n) < 2 ? 0 :			\
	 63 - __builtin_clzll(n)) :	\
	(sizeof(n) <= 4) ?		\
	__ilog2_u32(n) :		\
	__ilog2_u64(n)			\
 )

原來在 linux 核心追求的程式設計，即使是像 log2 這樣基本的函式，不只有正確性的考量，也需要考量硬體配置的差異 (如 32 位元及 64 位元)

測驗四

根據 EWMA 公式:

$\begin{array}{c}{S}_t=\{\begin{array}{l}{Y}_{0}t=0\\ \alpha Y_t+(1-\alpha )S_{t-1}t>0\end{array}\end{array}$

• $\alpha$ 為資料遞減的程度 (
  $0\le \alpha <1$)
• $Y_t$ 為在時間為 t 時的資料數值
• $S_t$ 為在時間為 t 時的 EWMA

對應到程式碼，有以下架構:

struct ewma {
    unsigned long internal;
    unsigned long factor;
    unsigned long weight;
};

internal 對應到

$S_t$，
由於平均數涉及到小數點，這裡使用 factor 表示定點數，weight 則對應

$\alpha$
於是 ewma_add 對應到上方公式

avg->internal = avg->internal
                    ? (((avg->internal << avg->weight) - avg->internal) +
                       (val << avg->factor)) >> avg->weight
                    : (val << avg->factor);

linux 核心案例

測驗五

TODO

第四週