# 2024q1 Homework4 (quiz3+4)
contributed by < `ssheep773` >

## 第 3 週測驗題
[2024q1 第 3 週測驗題](https://hackmd.io/@sysprog/linux2024-quiz3)

### 測驗 `1`
目的: 計算 `N` 的整數平方根

方法一 :
<!-- 首先用計算輸入 `N` 的二進位最高非 0 位元，其中 `msb` 儲存最高位元位置，而 `a` 代表 N 的最大位元的數值。

並在第二個 while 迴圈判斷逼近數 `reesult` 加上 $2^a$ 後的平方是否超過原數字`N` ，透過逐步減少 `a` 來逼近 `N` 的平方根。 -->

通過逐步逼近來計算 `N` 的平方根。它首先找到 `N` 的最大位元位置 `msb`，並將最大位元位置對應的數值儲存在 `a` $(a = 2^{msb})$ 。然後逐步減小 `a` 的值，同時檢查 $(result + a)^2$ 是否小於或等於 `N`。如果是的話，就將 `result` 加上 `a`，以逼近 `N` 的平方根。這樣的操作會持續進行，直到 `a` 的值為 0。

int i_sqrt(int N)
{
    int msb = 0;
    int n = N;
    while (n > 1) {
        n >>= 1;
        msb++;
    }
    int a = 1 << msb;
    int result = 0;
    while (a != 0) {
        if ((result + a) * (result + a) <= N)
            result += a;
        a >>= 1;
    }
    return result;
}

方法二 :
利用 [Digit-by-digit calculation](https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Digit-by-digit_calculation) 提及常見的直式開方法 $(x + y)^2 = x^2 + 2xy + y^2$ 的變形實作開平方根。

若要求 $x$ 的平方根，可以假設 $x = N^2$ ，那麼 $N$ 即為我們要求的平方根數值，接著，將 $N^2$ 拆成

$$
N^2=(a_{n}+a_{n-1}+a_{n-2}+...+{a_0})^2 = [(a_{n}+a_{n-1}+...+{a_1})+({a_0})]^2
$$

經過整理

$$
N^2 = (\sum_{i=0}^na_i)^2 = (\sum_{i=1}^{n} a_i)^2 + 2a_0(\sum_{i=1}^{n}a_i) + a_0^2
$$

這裡假設 $P_m = \sum_{i=m}^n a_i$ ，那麼 $P_0 = a_{n}+a_{n-1}+a_{n-2}+...+{a_0}$ 即為我們要求的 $N$， $P_0$ 可在分解為 $P_0 = (a_{n}+a_{n-1}+...+{a_1})+({a_0}) =P_1 + a_0$
可以得到一般式
$$
P_m = P_{m+1} +a_m
$$
我們的目的是要試出所有 $a_m$ ，$a_m = 2^m or 0$。從 $m = n$ 一路試到 $m = 0$ ，試驗 $P_m^2 \le N^2$ 是否成立來求得 $a_m$

$$
\begin{cases}
P_m = P_{m+1} +2^m \ \ ,if \ \ P_m^2 \le N^2\\
P_m = P_{m+1} \ \ ,if \ \ a_m=0\\
\end{cases}
$$

然而如果每次都試驗 $P_m^2 \le N^2$ 的運算成本會太高，一種想法是透過上一輪$N^2 - P_{m+1}^2$ 的差值推估 $N^2 - P_m$ 相關的證明如下:

因為前面 $P_m$ 的假設 $N$ 可以整理為
$$
N^2 = P_0^2 = a_0^2 + 2a_0P_1 + P_1^2
$$
並可推展出一般式
$$
P_m^2 = a_m^2 + 2a_mP_{m+1} + P_{m+1}^2
$$
經過移項
$$
P_m^2 -P_{m+1}^2 = 2a_mP_{m+1} + a_m^2
$$
假設 $Y_m = P_m^2 -P_{m+1}^2 = 2a_mP_{m+1} + a_m^2$

以及假設 $X_m = N^2 -P_m^2$ ，$X_m\geq0$ 即為每次判斷 $a^m$ 的條件，可以推廣 $X_{m+1} = N^2 -P_{m+1}^2$ 為前一次的差值
$$
Y_m = P_m^2 -P_{m+1}^2 = (N^2 -X_m) - (N^2 -X_{m+1}) = X_{m+1}-X_m
$$
我們就可以得到，從上一輪的差值 $X_{m+1}$ 減去 $Y_m$ 得到 $X_m$
$X_m = N^2 -P_m^2 = X_{m+1} - Y_m$
<!-- $Y_m = P_m^2 -P_{m+1}^2 = 2a_mP_{m+1} + a_m^2$
-->
將原本 $X_m \geq 0$ 的問題轉換為 $(X_{m+1} - Y_m) \geq 0$ 即使用到上一輪的差值來計算。
進一步調整，將 $Y_m$ 拆成 $c_m$ 和 $d_m$
- $c_m = 2P_{m+1}a_{m} = P_{m+1}2^{m+1}$
- $d_m = a_{m}^2 = (2^{m})^2$
- $Y_m =
\begin{cases}
c_m + d_m 　　, if 　a_m = 2^m \\
0　　　　　　, if 　a_m = 0 \\
\end{cases}$

以此推出下一輪
- $c_{m-1} = 2P_m a_{m-1} = P_m 2^m = (P_{m+1} + a_m)2^m = P_{m+1}2^m + a_m2^m =
\begin{cases}
c_m/2 +d_m 　, if 　a_m = 2^m\\
c_m/2　　　　, if 　a_m = 0 \\
\end{cases}$
- $d_{m-1} = (2^{m-1})^2 = \dfrac{d_m}{4}$

這裡說明初始與中止條件，因為是判斷 $(X_{m+1} - Y_m) \geq 0$ ，所以我們需要準備兩者的初始條件，在程式進行時，會由 $a_n$ 往下試
- $X_{n+1} = N$
- $n$ 是最大的位數，使的 $d_n^2 = a_n^2 = (2^n)^2 = 4^n \leq N^2$ 所以用迴圈逼近 $d_n$
- $c_n = 0$ 因為 $c_n = 2P_{n+1}a_{n}$ 而 $a_n$ 是最高位數所以 $P_{n+1} = 0$

因為 $d_m$ 恆大於 $0$ 所以終止條件可以設為 `d != 0` ，又 $c_{-1} = P_02^0 = P_0 = a_n + a_{n-1} +...+ a_0$ 故求出 $c_{-1}$ 則等於求出 $N$ 的平方根。

最後看到程式的部分

int i_sqrt(int x)
{
    if (x <= 1) /* Assume x is always positive */
        return x;

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

首先排除 `x` 為負數及 0 或 1 的情況

if (x <= 1) /* Assume x is always positive */
        return x;

變數 `m` 對應到 $d_m$ ，其中因為 m = $d_m = a_{m}^2 = (2^{m})^2$ ，所以 m 不可能為奇數，因此要 `& ~1UL` ，同時再前一個步驟也先處理 `x = 1` 的情況
因為 $d_{m-1} = \dfrac{d_m}{4}$ 所以每次迴圈 `m` 右移 2 。

for (int m = 1UL << ((31 - __builtin_clz(x)) & ~1UL); m; m >>= 2)

`x` 對應到 $X_m$ 呼應初始條件 $X_{n+1} = N$，變數 `b` 對應到 $Y_m$ ，`z` 對應到 $c_m$ 同時也說明 $Y_m = c_m + d_m$

而計算 $c_{m-1} =
\begin{cases}
c_m/2 +d_m 　, if 　a_m = 2^m\\
c_m/2　　　　, if 　a_m = 0 \\
\end{cases}$ ，不管 $a_m = 2^m$ 是否成立 $c_m$ 都需要除以二，這裡以 `z >>= 1` 實作

且若 則減去 $Y_m$ 得 $X_{m-1}$。

前面提到我們將原本 $X_m \geq 0$ 的問題轉換為 $(X_{m+1} - Y_m) \geq 0$ 同時也等價 $X_{m+1} \geq Y_m$ ，若成立則 $X_{m-1} = X_m - Y_m$

int b = z + m;
z >>= 1;
if (x >= b)
    x -= b, z += m; 

#### ffs / fls 取代 __builtin_clz

`ffs` 是找最低為 1 的位元，所以透過迴圈的方式由低往高位元找，到最高的有效位元的位置為止，並且還需要將其減一才會與 `(31 - __builtin_clz(x)` 相等，若沒有減一，雖對於結果沒有影響，但是會多花費一輪的運算成本。

int i_sqrt(int x)
{
    if (x <= 1) /* Assume x is always positive */
        return x;

+       int tmp = x, 
+       int msb = 0;
+       while (tmp) {
+           int i = ffs(tmp);
+           msb += i;
+           tmp >>= i;
+       }
        
    int z = 0;
+   for (int m = 1UL << ((msb-1) & ~1UL); m; m >>= 2) {
-   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;               
    }
    return z;
}

### 測驗 `2`
針對正整數在相鄰敘述進行 `mod 10` 和 `div 10` 操作，減少運算成本
原始的方法:

static void str_add(char *b, char *a, char *res, size_t size)
{
    int carry = 0;

    for (int i = 0; i < size; i++) {
        int tmp = (b[i] - '0') + (a[i] - '0') + carry;
        carry = tmp / 10;
        tmp = tmp % 10;
        res[i] = tmp + '0';
    }
}

可利用除法原理將 mod 10 和 div 10 合併。根據除法原理：
$f = g \times Q +r$
其中 $f:$ 被除數， $g:$ 除數， $Q:$ 商， $r:$ 餘數
將程式碼改寫為

    carry = tmp / 10;
-   tmp = tmp % 10;
+   tmp = tmp - carry * 10

這邊參考《Hacker's Delight》，採用 bitwise operation 來實作 10 的除法，但是因為 10 有 5 這個因數，無法完全用 2 的冪項來表示，因此結果會是不精確的。

我們的目標是，得到 tmp / 10 且直到小數點後第一位都是精確的。

為此，提出一個猜想，假設我目標的最大數是 $n$ ，而 $l$ 是小於 $n$ 的非負整數，則只要以 $n$ 算出來的除數在 $n$ 除以該除數後在精確度內，則 $l$ 除與該除數仍然會在精確度內，證明如下：
假設 $n=ab$ ( $a$ 是十位數 $b$ 是個位數)，$l=cd$ ( $c$ 是十位數 $d$ 是個位數)。

首先看到 $\frac{n}{x}$，因為我們是要確保小數點後第一位是精確的，所以 $\frac{n}{x}$ 可以表示如下 :

$a.b\leq\frac{n}{x}\leq a.b9\\
\Rightarrow\frac{n}{a.b9}\leq x\leq\frac{n}{a.b}$

透過上面的轉換，我們將問題轉換為確保 $x$ 的範圍，來達到 $\frac{n}{x}$ 在精度以內，下面是對 $x$ 範圍的探討

$x$ 的右不等式 $x \leq \frac{n}{a.b} = \frac{ab}{a.b} = 10$ ，因此一定在精確度內。

$x$ 的左不等式 $\frac{n}{a.b9}\leq x$ 無法直接證明，透過代入 $l$ 來證明，$l$的不等式如下:

$c.d\leq\frac{l}{x}\leq c.d9$ ，透過將 $\frac{n}{a.b9}\leq x$ 代入可表示為

$$
c.d\leq l \times \frac{a.b9}{n} \leq c.d9 \Rightarrow c.d\leq cd \times \frac{a.b9}{ab} \leq c.d9
$$

左不等式：
$$c.d\leq cd\times\frac{a.b9}{ab} = cd\times (\frac{a.b}{ab}+ \frac{0.09}{ab}) = c.d + \frac{0.09cd}{ab}
$$
最後推得 $c.d\leq c.d + \frac{0.09cd}{ab}$ 左式顯然成立。

右不等式：
$$cd \times \frac{a.b9}{ab} \leq c.d9 \Rightarrow c.d + \frac{0.09cd}{ab}\leq c.d9 \Rightarrow c.d + \frac{0.09cd}{ab}\leq c.d + 0.09
$$
其中 $\frac{0.09cd}{ab}\leq 0.09$ ，因為 $ab > cd$ 所以右不等式也成立

於是前述猜想也成立，接下來只需要針對 $n$ 來做討論即可，而因為 `int tmp = (b[i] - '0') + (a[i] - '0') + carry` 所以 `tmp` 最大是 9+9+1=19，$n$ 也即為 19 ，下面針對 `19` 做討論即可

$$
1.9\leq \frac{19}{x}\leq 1.99\\
\Rightarrow 9.55\leq x\leq10
$$

上述不等式可得知除數介於 $9.55$ 至 $10$ 之間皆可程式中達到相同效果，即除數直到小數點後第一位都是精確的，而我們的商會比實際的還要大有正偏差。

我們透過 bitwise operation 找到介於 $9.55$ 至 $10$ 之間的除數，除數的算式必定是 $\frac{2^N}{a}$ ，其中 $N$ 是非負整數， $a$ 是正整數，而計算 $n$ 的商可以寫成 $\frac{an}{2^N}$。

因此只需要透過查看 $2^N$ 在配對適合的 $a$ 即可
這裡 $2^N = 128, a= 13 ,\frac{128}{13} \approx 9.84$ 即為可用的除數解

接著看到實際的 bitwise operation

unsigned d2, d1, d0, q, r;

d0 = q & 0b1;
d1 = q & 0b11;
d2 = q & 0b111;
q = ((((tmp >> 3) + (tmp >> 1) + tmp) << 3) + d0 + d1 + d2) >> 7;
r = tmp - (((q << 2) + q) << 1);

`q` 是要求的商 $\frac{an}{2^N}$ ，在這裡是 $\frac{13tmp}{2^7}$，我們首先要得到 13 ，透過 `(tmp >> 3) + (tmp >> 1) + tmp) << 3` 得到 $\frac{13tmp}{8} = \frac{tmp}{8}+\frac{4tmp}{8}+tmp$ 再 `<< 3` 得到 $13tmp$，最後再除以 128 即得到我們要的商 $\frac{13tmp}{128}$

`(((q << 2) + q) << 1)` 這部分是將 `q * 10` 透過 `(q*4 + q) * 2` 達到

而因為原有的 `tmp` 在右移後，會有部分位元被捨棄，因此要另行處理。

    d0 = q & 0b1;
    d1 = q & 0b11;
    d2 = q & 0b111;

然而我覺得這步驟可以省略，因為後續的操作會在將 `tmp >> 7` 等於是補回來的位元馬上又被捨棄，而我思考若是在產生 13 而右移的位數大於 $2^N$ 的 $N$ 時，則需要處理，而目前分別是 3 跟 7 ，所以可以省略不處理。

在來看到包裝的函式

#include <stdint.h>
void divmod_10(uint32_t in, uint32_t *div, uint32_t *mod)
{
    uint32_t x = (in | 1) - (in >> 2); /* div = in/10 ==> div = 0.75*in/8 */
    uint32_t q = (x >> 4) + x;
    x = q;
    q = (q >> 8) + x;
    q = (q >> 8) + x;
    q = (q >> 8) + x;
    q = (q >> 8) + x;

    *div = (q >> 3);
    *mod = in - ((q & ~0x7) + (*div << 1));   
}

`uint32_t x = (in | 1) - (in >> 2)` 將 $x =\frac{3}{4}in$
再透過 `uint32_t q = (x >> 4) + x` 等於是 $q = \frac{3}{4*2^4}in + \frac{3}{4}in = \frac{102}{128}in$ ，而其中 $\frac{102}{128} \approx 0.797 \approx \frac{8}{10}$

後續再用 `q = (q >> 8) + x` 對 `q` 做逼近增加精度更靠近 $\frac{8}{10}$

`*div = (q >> 3)` 是計算商，而我們前面求的 `q` 是 $\frac{8}{10}in$ ，所以需要再除以 8 也就是 `q >> 3` 。

`*mod = in - ((q & ~0x7) + (*div << 1))` 是計算餘數，其中的 `q & ~0x7` 等於是 `*div << 3` ，那麼就程式碼可以轉換成 `*mod = in - (*(div << 3) + (*div << 1))` 後半部的 `*(div << 3) + (*div << 1)` 等於是 $商\times8 + 商\times2 = 商 \times (8+2) = 商 \times10$ ，根據餘數定理 $餘數 = 被除數(in) - 商(div) \times10 除數(10)$

### 測驗 `3`
ilog2 可計算以 2 為底的對數，且其輸入和輸出皆為整數，這裡注意到因為輸出為整數，所以輸入的最高二進位有效位元位置即為結果。

方法一 :
由最低位元往高位元處找最高二進位有效位元位置

int log = -1;
while (i) {
    i >>= 1;
    log++;
}

方法二 :
透過二元搜尋的方法尋找

static size_t ilog2(size_t i)
{
    size_t result = 0;
    while (i >= 2^16) {
        result += 16;
        i >>= 16;
    }
    while (i >= 2^8) {
        result += 8;
        i >>= 8;
    }
    while (i >= 2^4) {
        result += 4;
        i >>= 4;
    }
    while (i >= 2) {
        result += 1;
        i >>= 1;
    }
    return result;
}

方法三:
利用 GNU extension `__builtin_clz` ，找出最高有效位元前面 0 的數量，因此 `31 - __builtin_clz(v | 1)` 即為最高有效位元的位置，而因為 `__builtin_clz` 輸入若是 0 則無定義，所以使用 `v | 1` 確保輸入不為 0 。

int ilog32(uint32_t v)
{
    return (31 - __builtin_clz(v | 1));
}

在 [linux/log2.h](https://github.com/torvalds/linux/blob/39cd87c4eb2b893354f3b850f916353f2658ae6f/include/linux/log2.h#L30) 中使用 `fls` 找出最高有效位元並且說明輸入 0 是未定義的。

/*
 * non-constant log of base 2 calculators
 * - the arch may override these in asm/bitops.h if they can be implemented
 *   more efficiently than using fls() and fls64()
 * - the arch is not required to handle n==0 if implementing the fallback
 */
#ifndef CONFIG_ARCH_HAS_ILOG2_U32
static __always_inline __attribute__((const))
int __ilog2_u32(u32 n)
{
	return fls(n) - 1;
}
#endif

#ifndef CONFIG_ARCH_HAS_ILOG2_U64
static __always_inline __attribute__((const))
int __ilog2_u64(u64 n)
{
	return fls64(n) - 1;
}
#endif

### 測驗 `4`
[Exponentially Weighted Moving Average](https://en.wikipedia.org/wiki/Moving_average#Exponential_moving_average) (EWMA; 指數加權移動平均) 是種取平均的統計手法，並且使經過時間越久的歷史資料的權重也會越低，以下為 EWMA 的數學定義:

$$
S_t =
\begin{cases}
\ Y_0 \qquad\qquad\qquad\qquad ,t=0\\
\ \alpha Y_t +(1-\alpha) \cdot S_{t-1}
\ \ ,t>0\\
\end{cases}
$$

$\alpha$ 表示歷史資料加權降低的程度，介在 0 ～ 1 之間，越高的 $\alpha$ 會使歷史資料減少的越快
$Y_t$ 表示在時間 $t$ 時的資料點
$S_t$ 表示在時間 $t$ 時計算出的 EWMA

首先看到 ewma 結構，`internal` 儲存我們的 EWMA 值也就是 $S_t$ ，值得注意的是它的型態是 `unsigned long` 是無號整數，所以我們可以知道這裡是使用定點數，並在 `ewma_init()` 註解處中有說明 `factor` 是 Scaling factor ， 其中註解處有提到`internal` 可記錄最大值的計算公式為 `ULONG_MAX / (factor * weight)` ，其中 `factor` 是因為做 scaling 提高數值精度變相會犧牲可容納的對大值，而 `weight` 的部分不太理解，只能夠推測因為程式中有做 `avg->internal << avg->weight` 與 `avg->internal >> avg->weight` 與 `factor` 一樣。

並在 `ewma_read` 可以看到輸出時還原 scaling ，可見 `avg->internal >> avg->factor`

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

接著看到 `ewma_init` 函式，目的是初始化 `ewma`，其中使用 `is_power_of_2` 來確保 `weight` 跟 `factor` 是 2 的冪，因為 `internal` 是定點數若是與 2 的冪做乘除法，可以用位元運算。
`is_power_of_2` 是運用當數值 `n` 是 2 的冪時，`n` 跟 `n-1` 在二進位時，不會有相同的位數會是錯開的，利用這個特點檢測

static inline int is_power_of_2(unsigned long n)
{
    return (n != 0 && ((n & (n - 1)) == 0));
}

void ewma_init(struct ewma *avg, unsigned long factor, unsigned long weight)
{
    if (!is_power_of_2(weight) || !is_power_of_2(factor))
        assert(0 && "weight and factor have to be a power of two!");

    avg->weight = ilog2(weight);
    avg->factor = ilog2(factor);
    avg->internal = 0;
}

$$
S_t =
\begin{cases}
\ Y_0 \qquad\qquad\qquad\qquad ,t=0\\
\ \alpha Y_t +(1-\alpha) \cdot S_{t-1}
\ \ ,t>0\\
\end{cases}
$$

最後看到 `ewma_add` ，這裡可想而知會用到前面提到的 $S_t$ ，先說明 `avg->internal` 等於 0 的情況，將 `avg->internal` 設為輸入的 `val` 對應到 $Y_0$ ，但別忘了 scaling 所以還要 `val << avg->factor`

在來看到大於 0 的情況，第一次看到 `(((avg->internal << avg->weight) - avg->internal) +(val << avg->factor)) >> avg->weight` 時感覺跟 $\alpha Y_t +(1-\alpha) \cdot S_{t-1}$ 公式有關，但是仔細看時發現正負號會有錯，這邊再看過 [ShchangAnderson](https://hackmd.io/@ShchangAnderson/linux2024-homework4#%E6%B8%AC%E9%A9%97%E5%9B%9B) 同學的筆記後才了解，其實 $\alpha = \frac{1}{2^{avg->weight}}$，所以程式是先將
$$
\alpha Y_t +(1-\alpha) \cdot S_{t-1}
= [ \alpha Y_t +(1-\alpha) \cdot S_{t-1} ] \times 2^{avg->weight} \times \frac{1}{2^{avg->weight}}
$$
將 $2^{avg->weight}$ 乘進去原始公式，再簡單地移項方便對證程式碼，可以得到
$$
[ 2^{avg->weight}\cdot S_{t-1} - S_{t-1} + Y_t] \times \frac{1}{2^{avg->weight}}
$$
其中先乘以 $2^{avg->weight}$ 再乘 $\frac{1}{2^{avg->weight}}$ 是為了提高精度，程式即是使用上述的公式做計算，並且在每次輸入 `val` 都會 scaling

struct ewma *ewma_add(struct ewma *avg, unsigned long val)
{
    avg->internal = avg->internal
                        ? (((avg->internal << avg->weight) - avg->internal) +
                           (val << avg->factor)) >> avg->weight
                        : (val << avg->factor);
    return avg;
}

### 測驗 `5`
這個程式碼實做一個函式，用來計算輸入的 32 位元無號整數 $x$ 的 2 次方指數值向上進位的結果。也就是說，對於傳入的參數 $x$ ，回傳最小的整數 $n$，滿足 $x\leq2^n$ 。

int ceil_ilog2(uint32_t x)
{
    uint32_t r, shift;

    x--;
    r = (x > 0xFFFF) << 4;
    x >>= r;

    shift = (x > 0xFF) << 3;
    x >>= shift;
    r |= shift;

    shift = (x > 0xF) << 2;
    x >>= shift;
    r |= shift;

    shift = (x > 0x3) << 1;
    x >>= shift;
    
    return (r | shift | x > 1) + 1;  

在函式的開始將 x 減 1，這樣可以確保當 x 是 2 的冪時的正確性，因為 $x\leq2^n$，當等於時不用進位。

接著看到

shift = (x > 0xFF) << 3;
x >>= shift;
r |= shift;

這裡其實跟測驗 `3` 是相同的概念，將測驗三的變數調整方便比較，其中 `0xFFFF = 2^16` ，並且將數字都以 2 的冪表示。
此方法與測驗三的差異在於紀錄 `r` 的方式是透過 bitwise 操作而非加法運算， `r | shift` 等效於 `r + shift` 。

/*   ilog2   */
while (x >= 2^8) {
    r += 2^3;
    x >>= 2^3;
}

最後 `return (r | shift | x > 1) + 1` 的部分我們分開看，其中 `r | shift` 是將前一部分的 `r |= shift` 合併進來， `x > 1` 是為了要處理前面 `shift = (x > 0x3) << 1` 會忽略 `x= 0x2` 的情況，最後再加上一作為取上界 (ceil) 。

#### 改進程式碼，使其得以處理 x = 0 的狀況，並仍是 branchless

當 `x=0` 時會因為減一變成 `0xFFFFFFFF` ，而這不是我們想要的結果。
我想到在 `likely / unlikely` 中有用到 `!!(x)` 將整數輸入結果控制在 `0` 跟 `1` ， 那麼就可以將 `x--` 變為 `x = x - !!x` ，當 `x > 0` 時減一，而當 `x = 0` 時則不變。

int ceil_ilog2(uint32_t x)
{
    uint32_t r, shift;

-   x--;
+   x = x - !!x
    r = (x > 0xFFFF) << 4;
    x >>= r;

## 第 4 週測驗題
[2024q1 第 4 週測驗題](https://hackmd.io/@sysprog/linux2024-quiz4)

### 測驗 `1`
[population count](https://en.wikichip.org/wiki/population_count) 簡稱 popcount 或叫 sideways sum，計算數值的二進位表示中 1 的個數。

unsigned popcount_naive(unsigned v)
{
    unsigned n = 0;
    while (v)
        v &= (v - 1), n = -(~n);
    return n;
}

`v &= (v - 1)` : reset LSB 是將 最低有效位 (LSB) 設為 0 。這讓我想到這跟前面提到的 `is_power_of_2` 觀念類似，都是判斷 `(n & (n - 1)) == 0` 只是 `is_power_of_2` 只判斷一次，而 `popcount_naive` 是重複判斷來計算 1 的個數。
`n = -(~n)` 等同 `n++`，因為在二補數系統中 $-n = \sim n+1$ ，經過移項 $-(\sim n) = n +1$

#### 針對 [LeetCode 477. Total Hamming Distance](https://leetcode.com/problems/total-hamming-distance/description/) 撰寫出更高效的程式碼。

透過觀察 Total Hamming Distance 的規則

可以得到一般式:
觀察每個 Number 的第 i 個 bit ，若有 `cnt` 個 1 ，其餘皆為 0 ，則 Hamming Distance : `(cnt) * (numsize - cnt)`

我們要走訪所有的 bit ，在每次走訪時計算第 i 個 bit 為 1 的 number 個數(`cnt`)，其中使用 `(nums[j] >> i) & 1` 得知目前的number 第 i 個 bit 是否為 1 ，最後透過公式 `(cnt) * (numsize - cnt)` 計算得到第 i 個 bit 的 Hamming Distance ，最後加總所有位元的 Hamming Distance 即為 nums 中所有整數的 Hamming Distance 。

int totalHammingDistance(int* nums, int numsSize) {
    int total = 0;;
    for (int i = 0;i < 32;i++){
        int cnt = 0;
        for (int j = 0; j < numsSize;j++)
            cnt += ((nums[j] >> i) & 1);
        total += cnt *(numsSize - count);
    }
    return total;
}

### 測驗 `2`
Remainder by Summing digits
如何不使用任何除法就算出某數除以另一個數的餘數呢？若除數符合 $2^k \pm 1$ ，則可運用以下手法來達成這個目的。

以除數為 3 為例，我們可以觀察到 $1 \equiv 1 (mod \ \ 3)$ 且 $2 \equiv -1 (mod \ \ 3)$。 將後者不斷的乘上前者可得

$$2^k \equiv \begin{cases}
1 (mod \ \ 3), \ \ k \ even\\
-1 (mod \ \ 3), \ \ k \ odd \\
\end{cases}$$

因此若 n 的二進位表示為 $b_{n-1}b_{n-2}b_{n-3}...b_1b_0$ ，則 $n = \sum_{i=0}^{n-1} b_i \cdot 2^i$ ，由以上的推論可得 $n \equiv \sum_{i=0}^{n-1} b_i \cdot (-1)^i \ \ (mod \ \ 3)$ 將奇偶數分開 $n \equiv \sum_{i=even} b_i - \sum_{i=odd} b_i \ \ (mod \ \ 3)$

可以利用 population count 來計算位元和，因此可以將上式表示為 `n = popcount(n & 0x55555555) - popcount(n & 0xAAAAAAAA)` ，其中 `& 0x55555555` 是留下 n 的奇數位元， `& 0xAAAAAAAA` 是留下 n 的偶數位元，並利用以下定理可以進一步化簡。

$$
popcount(x \& \overline{m}) - popcount(x \& m) = popcount(x \oplus m) - popcount(m)
$$

於是可以寫為 `n = popcount(n ^ 0xAAAAAAAA) - 16` ，計算結果的範圍會落在 -16 至 16 之間 ，但為了避免 n 變成負數，我們要將它加上一個足夠大的 3 的倍數，文中的例子是加上 39 讓數值介於 23 到 55 ，我們再做一次與上面相同的操作，只是這次的數值較小，相差 32 所以跟 `0x2A = 0b101010` 做運算而 `popcount(0x2A) = 3` 所以是減三，最後區間會介於 `-3 <= n <= 2` ，還需要對負數加 3 ，所以最後 `(n >> 31) & 3` 判斷最高位元判斷正負來決定是否加 3 。

int mod3(unsigned n)
{
    n = popcount(n ^ 0xAAAAAAAA) + 23;  // Now 23 <= n <= 55.
    n = popcount(n ^ 0x2A) - 3;  // Now -3 <= n <= 2
    return n + ((n >> 31) & 3);
}

另一種變形是利用 lookup table ，其中 n 介於 0 到 32 透過查表直接得知。

int mod3(unsigned n)
{
    static char table[33] = {2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1 };
    n = popcount(n ^ 0xAAAAAAAA);
    return table[n];
}

然而在 [HackerDelight](https://web.archive.org/web/20180517023231/http://www.hackersdelight.org/divcMore.pdf) 中提到不使用指令的方法

int remu3(unsigned n) {
 n = (n >> 16) + (n & 0xFFFF); // Max 0x1FFFE.
 n = (n >> 8) + (n & 0x00FF); // Max 0x2FD.
 n = (n >> 4) + (n & 0x000F); // Max 0x3D.
 n = (n >> 2) + (n & 0x0003); // Max 0x11.
 n = (n >> 2) + (n & 0x0003); // Max 0x6.
 return (0x0924 >> (n << 1)) & 3;
}

這段程式碼主要是利用 $4^k \equiv 1 (mod \ \ 3)$ 簡化求 $n \ (mod \ \ 3)$
首先假設 n 前的 16 位元是 $c$ ，後 16 位元是 $d$ ，那麼我們就可以將 n 的值寫成 $n=c\cdot2^{16}+d$ ，再透過 $4^k \equiv 1 (mod \ \ 3)$ 可以將 $n (mod \ \ 3)$表示為

$$
n (mod \ \ 3) =
c\cdot2^{16}+d(mod \ \ 3) =
c+d(mod \ \ 3)
$$

上述的數學式可轉換成 `n = (n >> 16) + (n & 0xFFFF)` ，以此類推重複 5 次將 n 縮小至最大 `0x6`
最後一行類似查表，將餘數值存在 `0x0924` ，透過位元操作取值，其中因為 3 的餘數是 {0,1,2} 可以用 2 個位元即可表示，所以 `0x0924` 是每隔 2 位元存一個餘數值，並透過 `(n << 1)` 與 `& 3` 搭配取得。

0x0924   =  00 00 10 01 00 10 01 00
(n << 1) =     12 10  8  6  4  2  0   
 n       =      6  5  4  3  2  1  0   

`tictactoe.c` 中使用到的 mod7 方法說明
上述提到利用 $4^k \equiv 2^{2k}\equiv 1 (mod \ \ 3)$ 簡化求 $n \ (mod \ \ 3)$ ，同理利用 $8^k \equiv 2^{3k} \equiv 1 (mod \ \ 7)$ 簡化求 $n \ (mod \ \ 7)$ ，這裡 `(x >> 15)` 使用 15 是因為 $2^{3k}$ ，所以是每隔 3 位元有同餘的結果，而我們要將原 32 位元的數縮小，找 16 附近的 3 的倍數 15 。
在這個方法中不像 [HackerDelight](https://web.archive.org/web/20180517023231/http://www.hackersdelight.org/divcMore.pdf) 提到的例子，繼續縮小輸入的值後查表求餘數，而是透過 $n (mod \ 7) \equiv \lfloor\frac{8}{7}n\rfloor (mod \ 8)$ 求 7 的餘數，先將數字乘以 $\lceil\frac{2^{32}}{7}\rceil$ 也就是 `0x2492492` 再右移 29 得到剩餘的三位數即為解。
:::info
不太明白為何 $n (mod \ 7) \equiv \lfloor\frac{8}{7}n\rfloor (mod \ 8)$

:::

static inline int mod7(uint32_t x)
{
    x = (x >> 15) + (x & UINT32_C(0x7FFF));
    /* Take reminder as (mod 8) by mul/shift. Since the multiplier
     * was calculated using ceil() instead of floor(), it skips the
     * value '7' properly.
     *    M <- ceil(ldexp(8/7, 29))
     */
    return (int) ((x * UINT32_C(0x24924925)) >> 29);
}

`tictactoe.c` 是一個井字遊戲，他的特別之處在於儲存棋子的方式，以往可能是用一個矩陣儲存井字中每個位置的棋子狀態，這個程式是用棋子可以連成的線來儲存，這個方法的優點在於可以更快速的判斷贏家，也可以更方便取得連線的資訊，

### 測驗 `3`