2024q1 Homework4 (quiz3+4)

contributed by < yeh-sudo >

第三周測驗題

測驗 1

解釋程式碼運作原理

假設要求的平方根為

$x$ ，並且將

$x^2$ 表示為以下形式：

$$x^2=(000b_0{b}_1{b}_2...b_{n-1}{b}_n{)}^2$$

其中

$000b_0{b}_1{b}_2...b_{n-1}{b}_n$ 為 bit pattern ，

$b_0$ 為最高位的 1 。所以

$x^2$ 可以寫成：

$$x^2=(b_0{2}^n+b_1{2}^{n-1}+b_2{2}^{n-2}+...+b_{n-1}{2}^1+b_n{2}^0{)}^2$$

設

$a_0=b_0{2}^n,a_1=b_1{2}^{n-1},a_2=b_2{2}^{n-2}...a_{n-1}=b_{n-1}{2}^1,a_n=b_n{2}^0$，

$x^2$ 可以表示成：

$$x^2=(a_0+a_1+a_2+...+a_{n-1}+a_n{)}^2$$

根據 Digit-by-digit calculation ：

$$\begin{array}{rl}{x}^2& =(a_n+a_{n-1}+a_{n-2}+...+a_1+a_0{)}^2\\ & =a_n^2+2a_n{a}_{n-1}+a_{n-1}^2+2(a_n+a_{n-1})a_{n-2}+a_{n-2}^2+...+a_1^2+2(\sum _{i=1}^n{a}_i)a_0+a_0^2\\ & =a_n^2+[2a_n+a_{n-1}]a_{n-1}+[2(a_n+a_{n-1})+a_{n-2}]a_{n-2}+...+[2(\sum _{i=1}^n{a}_i)+a_0]a_0\end{array}$$

設

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

$x^2$ 可以表示為：

$$x^2=a_n^2+[2P_n+a_{n-1}]a_{n-1}+[2P_{n-1}+a_{n-2}]a_{n-2}+...+[2P_1+a_0]a_0$$

而

$P_0=a_n+a_{n-1}+a_{n-2}+...+a_1+a_0$ 即為要求出的平方根。所以，從

$m=n$ 一路試到

$m=0$ ，測試

$P_m^2\le x^2$ 是否成立，分為兩種情況：

$$\{\begin{array}{ll}{P}_m=P_{m+1}+a_m& \text{if }{P}_m^2\le x^2\\ P_m=P_{m+1}& \text{otherwise}\end{array}$$

在檢驗時，如果每次都要計算

$P_m$ ，會需要很高的成本，所以可以利用重複利用每一次計算的差值

$X$ ，可以得出：

$$\begin{array}{rl}{X}_n& =x^2-P_n^2\\ & =x^2-a_n^2\end{array}$$

$$\begin{array}{rl}{X}_{n-1}& =x^2-P_{n-1}^2\\ & =x^2-(a_n+a_{n-1}{)}^2\\ & =x^2-a_n^2-a_{n-1}^2-2a_n{a}_{n-1}\\ & =(x^2-a_n^2)-a_{n-1}^2-2a_n{a}_{n-1}\\ & =X_n-a_{n-1}^2-2a_n{a}_{n-1}\\ & =X_n-[2a_n+a_{n-1}]a_{n-1}\\ & =X_n-[2P_n+a_{n-1}]a_{n-1}\end{array}$$

$$\begin{array}{rl}{X}_{n-2}& =x^2-P_{n-2}^2\\ & =x^2-(a_n+a_{n-1}+a_{n-2}{)}^2\\ & =x^2-a_n^2-a_{n-1}^2-a_{n-2}^2-2a_n{a}_{n-1}-2a_n{a}_{n-2}-2a_{n-1}{a}_{n-2}\\ & =X_n-[2P_n+a_{n-1}]a_{n-1}-2a_n{a}_{n-2}-2a_{n-1}{a}_{n-2}-a_{n-2}^2\\ & =X_{n-1}-[2(a_n+a_{n-1})+a_{n-2}]a_{n-2}\\ & =X_{n-1}-[2P_{n-1}+a_{n-2}]a_{n-2}\end{array}$$

$$\begin{array}{r}...\end{array}$$

$$\begin{array}{rl}{X}_m& =X_{m+1}-[2P_{m+1}+a_m]a_m\end{array}$$

設

$Y_m=[2P_{m+1}+a_m]a_m$ ，則可以寫成以下遞迴關係：

$$X_m=X_{m+1}-Y_m$$

將

$Y_m$ 拆成

$c_m$ 與

$d_m$ 兩部份，

$$\begin{array}{rl}{Y}_m& =[2P_{m+1}+a_m]a_m\\ & =2P_{m+1}{a}_m+a_m^2\\ & =2^{m+1}{P}_{m+1}+{(2^m)}^2\\ & =c_m+d_m\end{array}$$

$$\{\begin{array}{ll}{Y}_m=c_m+d_m& \text{if }{a}_m=2^m\\ Y_m=0& \text{otherwise}\end{array}$$

$$\{\begin{array}{l}{c}_m=2^{m+1}{P}_{m+1}\\ d_m={(2^m)}^2\end{array}$$

可以從

$c_m$ 與

$d_m$ 推出下一輪的

$c_{m-1}$ 與

$d_{m-1}$ ，再用

$c_{m-1}$ 與

$d_{m-1}$ 推出

$c_{m-2}$ 與

$d_{m-2}$ ，直到推出

$c_{-1}=P_0$ 即為所求的平方根

$x$ ，而首項

$c_n=0,d_n={(2^n)}^2$ 。

$$\begin{array}{rl}{c}_{m-1}& =2^m{P}_m\\ d_{m-1}& ={(2^{m-1})}^2={(2^m)}^2{2}^{-2}=d_m/4\end{array}$$

此處的

$P_m=P_{m+1}+a_m$ 成立的條件為

$P_m^2\le x^2$ ，若不成立，則

$P_m=P_{m+1}$ ，所以當

$P_m^2\le x^2$ 時，

$a_m=2^m$ ，反之則

$a_m=0$ 。所以最後可以將

$c_m$寫成以下關係：

$$\begin{array}{c}{c}_{m-1}=\{\begin{array}{ll}{c}_m/2+d_m& \text{if }{a}_m=2^m\\ c_m/2& \text{if }{a}_m=0\end{array}\end{array}$$

在程式碼中，

$c_m$ 對應到變數 z ，

$d_m$ 對應到變數 m ， 差值

$X$ 對應到變數 x 。由於是紀錄差值，所以每當差值大於等於

$c_m+d_m$ ，也就是 b = z + m 時，就將差值扣掉 b ，並且，因為差值為正，代表

$P_m^2\le x^2$ ，也就是

$a_m=2^m$ ，所以

$c_{m-1}=c_m/2+d_m$ ，對應到的程式碼操作為 z += m ，持續進行迴圈直到 m 為 0 時，結束迴圈。

int z = 0;
for (int m = 1UL << ((31 - __builtin_clz(x)) & ~1UL); m; m >>= 2) {
    int b = z + m;
    z >>= 1;
    if (x >= b)
        x -= b, z += m;
}

以 ffs / fls 取代 __builtin_clz

程式碼 (31 - __builtin_clz(x)) 的作用是找到最大的 bit 所在的位置，而使用第二週測驗題提供的 __ffs ，可以實作出找出對大 bit 的 fls 函式，實作方法是當 x 不為 0 時，持續尋找 ffs ，找到時就使用測驗題提供的 __clear_bit 將找到的 bit 清為 0 ，當 x 為 0 時，代表最後找到的 ffs 為最大的 bit 。

static inline unsigned long fls(unsigned long x)
{
    int result = 0;
    while (x) {
        unsigned int bit = __ffs(x);
        if (bit > result)
            result = bit;
        __clear_bit(bit, &x);
    }
    return result;
}

另外還有一個方法可以找出 fls ，將 x 的逐位元反轉，再用 31 減掉 ffs 即可得到 fls 。

static inline unsigned long fls(unsigned long x)
{
    unsigned long result = x;
    result = ((result & 0xffff0000) >> 16) | ((result & 0x0000ffff) << 16);
    result = ((result & 0xff00ff00) >> 8) | ((result & 0x00ff00ff) << 8);
    result = ((result & 0xf0f0f0f0) >> 4) | ((result & 0x0f0f0f0f) << 4);
    result = ((result & 0xcccccccc) >> 2) | ((result & 0x33333333) << 2);
    result = ((result & 0xaaaaaaaa) >> 1) | ((result & 0x55555555) << 1);
    return 31 - __ffs(result);
}

提供分支和無分支 (branchless) 的實作

在教材〈解讀計算機編碼〉中，有描述不需要分支的設計，當 INT32_MIN <= (b - a) <= INT32_MAX 條件成立時，可以使用 bitwise 操作改寫分支。

• 改寫 i_sqrt

-        if (x >= b)
-            x -= b, z += m;
+        int tmp = x;
+        x -= (~(tmp - b) >> 31) & b;
+        z += (~(tmp - b) >> 31) & m;

• 改寫 __ffs

 #if BITS_PER_LONG == 64
-    if ((word & 0xffffffff) == 0) {
-        num += 32;
-        word >>= 32;
-    }
+    num += !(word & 0xffffffff) & 32;
+    word >>= 32;
 #endif
-    if ((word & 0xffff) == 0) {
-        num += 16;
-        word >>= 16;
-    }
-    if ((word & 0xff) == 0) {
-        num += 8;
-        word >>= 8;
-    }
-    if ((word & 0xf) == 0) {
-        num += 4;
-        word >>= 4;
-    }
-    if ((word & 0x3) == 0) {
-        num += 2;
-        word >>= 2;
-    }
-    if ((word & 0x1) == 0)
-        num += 1;
+    num += !(word & 0xffff) & 16;
+    word >>= 16;
+    num += !(word & 0xff) & 8;
+    word >>= 8;
+    num += !(word & 0xf) & 4;
+    word >>= 4;
+    num += !(word & 0x3) & 2;
+    word >>= 2;
+    num += !(word & 0x1) & 1;
     return num;

Linux 核心原始程式碼

~/linux/block master ❯ grep sqrt *.c
blk-wbt.c: *   window is then shrunk to 100 / sqrt(scaling step + 1).
blk-wbt.c:                                      int_sqrt((rqd->scale_step + 1) << 8));

在 Linux 核心中， block 這個目錄下的 blk-wbt.c 有使用到整數平方根的操作，使用的函式為 rwb_arm_timer 。

• blk-wbt.c

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));
	} else {

閱讀實作 writeback throttling 這個機制的作者 Jens Axboe 的 commit messages ， throttling 這個機制最主要是讓背景執行的 buffered writeback 不會影響到其他正在執行的程式，作者發現在執行大量的 buffered writeback 時，他想要開啟 Chrome 會打不開。

But for as long as I can remember, heavy buffered writers have not behaved like that. For instance, if I do something like this:

$ dd if=/dev/zero of=foo bs=1M count=10k

on my laptop, and then try and start chrome, it basically won't start
before the buffered writeback is done. Or, for server oriented
workloads, where installation of a big RPM (or similar) adversely
impacts database reads or sync writes. When that happens, I get people
yelling at me.

這個機制有點像是 CoDel networking scheduling algorithm ， blk-wb 監控一段時間中的最小延遲，如果這段時間中，有最小延遲超過設定的目標，就增加 scaling step 與壓縮佇列的深度，接著將監控的時間縮小 100 / sqrt(scaling step + 1) 。其中的整數平方根運算對應的程式碼在 /lib/math/int_sqrt.c ，其實作的演算法也是使用 digit-by-digit calculation 。

• int_sqrt

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

閱讀〈Analysis of the lookup table size for square-rooting〉

閱讀〈Timing Square Root〉