2024q1 Homework4 (quiz3+4)
===
contributed by < `wu81177` >
# [第 3 週測驗題](https://hackmd.io/@sysprog/linux2024-quiz3)
## 測驗一
版本一當中，可以理解為使用了 binary search 來找平方根，程式先找到 `msb` 作為初始值，在 `msb` 和兩倍 `msb` 之間不斷二分來勘根，而因為整數型會截斷小數部分，因此 `a` 就會是小於等於 `n` 二的冪次最大值，之後若 `result` 加 `a` 後平方超過 `n` 就把 `a` 除以二查找，而由於 `a` 最小單位是 1 ，對於非平方數只能找到平方根的整數部分。
版本二和一的差別在於 `msb` 起始值的計算方式，版本一是接使用 `math.h` 的 `log2()` ，而版本二是將 n 值不斷除以二直到小於一 ( bitwise 的操作會使 n 最終變 0)，過程中用 `msb` 記數，一樣也能算出 $lg N$ 整數部分。
而版本三的部分，根據 [Digit-by-digit calculation](https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Decimal_(base_10)) ，所求 $N$ 可以寫成 $N^2 = (a_n+a_{n-1}+a_{n-2}+...+a_0)^2$ ，其中 $a_n$ 為二的 $n$ 次冪或0，可整理成:
$$
\begin{split} \\
N^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(\displaystyle\sum_{i=1}^{n}a_i)+a_0]a_0 \\
=& 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_m =& a_n+a_{n-1}+...+a_m \\
\end{split}
$$
可以發現 $P_0$ 肯定是所求，同時也存在某些 $m$ 滿足 $P_m = P_0$ ，因為某些 $a_m$ 是 0 ，令：
- $X_m = N^2 - P_m^2 = X_{m+1} - Y_m$
- $Y_m = P_m^2 - P_{m+1}^2 = 2P_{m+1}a_m+a_m^2$
-
上一輪底數為 $m+1$ 都是已知，同時將 $Y_m$ 拆成 $c_m, d_m$
$$
c_m = P_{m+1}2^{m+1} ,\\
d_m = (2^m)^2 \\
,Y_m=\left.
\begin{cases}
c_m+d_m & \text{if } a_m=2^m \\
0 & \text{if } a_m=0
\end{cases}
\right.
$$
就可以藉由位元運算來計算 $c_m, d_m$，迭代求出 $c_{-1}$ 即為所求 $N$ ，了解原理後，版本三 AAAA 應為 2 ，因為 $d_m = 4^m$， BBBB 應為 1 ，因為 $c_m = P_{m+1}2^{m+1}$
:::spoiler 版本三
```c
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;
}
```
:::
## 測驗二
為了使用 bitwise 取代 `/` 和 `%` ，由於 0.1 寫成二進位會是無窮小數，需要找到一個二進位值近似 0.1 ，目標是在 0~19 範圍內的 `temp` 乘上它到小數點後第一位以前都正確，而測驗說明中證明了只需要檢查範圍內 `temp` 最大值是否達標就好，因為近似值乘上 `temp` 誤差會隨著 `temp` 等比變化，若是最大值在誤差範圍內，絕對值較小值想必也會在範圍內
::: spoiler 證明
假設目標的最大數是 `n` ， `l` 則是比 `n` 還要小的非負整數。假設 $n = ab$ ($a$ 是十位數 b 是個位數)，且 $l = cd$ ($c$ 是十位數，$d$ 是個位數)，則只要以 $n$ 算出來的除數在 $n$ 除以該除數後在精確度內，則 $l$ 除與該除數仍然會在精確度內，證明如下：
假設 $x$ 是近似值，
$$
a.b\leq\frac{n}{x}\leq a.b9\\
\Rightarrow\frac{n}{a.b9}\leq x\leq\frac{n}{a.b}
$$
1. 可見 $x$ 的右邊是 $10$，因此一定在精確度內。
2. $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\\
$$
* $cd\times\frac{a.b9}{ab}$ 的左邊顯然成立
* $cd\times\frac{a.b9}{ab}$ 的右邊：
$$
c.d + \frac{0.09cd}{ab}\leq c.d + 0.09
$$
因為 $ab > cd$ 因此上述式子依然成立。
:::
總之在這個例子中只需要檢查 `temp` 為 19 的情況就好。
0.1 近似值必須滿足一些條件，以有理數 $a/b$ 表示，其中 $b$ 為二的冪，就能用 bitwise 來運算，運算方式如下:
將 $a$ 轉為二進位 $a_n...a_2 a_1 a_0$ ，$b$ 寫為 $2^N$，可以將 $tmp*(a/b)$ 寫為以下通式:
$$
\frac{a_ntmp}{2^{N-n}}+...+\frac{a_2tmp}{2^{N-2}}+\frac{a_1tmp}{2^{N-1}}+\frac{a_0tmp}{2^{N}}
$$
分母 2 的指數就是要右移的位數，以 $13/8$ 為例，$N = 3$，$a = (1101)_2$ ，可以寫為， d 是為了保留右移被捨棄的位元 :
```
d0 = q & 0b1;
d1 = q & 0b11;
d2 = q & 0b111;
((((tmp >> 3) + (tmp >> 1) + tmp) << 3) + d0 + d1 + d2)
```
從上述通式能看出若要找 0.1 近似值， $a$ 取 $b/10$ 四捨五入至整數位誤差最小，但為了滿足小數點後第一位以前正確， $a$ 應該取 $b/10$ 無條件進位值，這樣誤差才會是正的，而 $N$ 越大誤差也就越小，如圖，橫軸為 N 縱軸為 0.1近似值*19的誤差:
可以看到當 $N>=7$ 就能滿足要求，然而 a 實際上可為小數。
接下來看題目實作的函式:
```clike#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 >> CCCC);
*mod = in - ((q & ~0x7) + (*div << DDDD));
}
```
` uint32_t x = (in | 1) - (in >> 2)` 是將 `in` 減去 1/4 存為 `x` ，但不懂為何要 `in | 1`，而第二行 q 約為 $in(\frac{3}{4}*\frac{1}{16}+\frac{3}{4})$ 等於 $0.796875 in$ ，之後的5行讓 q 超級接近 0.8 ，這時我才發現它想讓近似值的分子為 0.8，我在前面的討論沒有意識到分子可為小數。
因為分子是 0.8 ，分母應該要為 8 ，所以依照上面的規則 CCCC 應為 3 。
`q & ~0x7` 的效果是將後面三位設成0，相當於 $0.8 in /8*8 = 8*div$ ，因此後面那項應該要是 $2*div$ 加起來才會是 $10*div$， DDDD 應為 1 ，我的疑惑是不知為何 `(q & ~0x7) + (*div << 1)` 兩項不統一為其中一種寫法就好，不知道這樣有什麼好處。
## 測驗三
版本一的做法是逐位右移，相當於是一直除以二並且計數，直到等於0就能求出二為底的對數
版本二的作法也一樣，差別是對於很大的數檢查如果大於某數就一次右移很多位，這樣計算的次數比起版本一會少很多，所以 AAAA 應為 $2^{16} = 65536$，BBBB 應為 $2^8 = 256$ ， CCCC 為 $2^4 = 16$
```c
static size_t ilog2(size_t i)
{
size_t result = 0;
while (i >= AAAA) {
result += 16;
i >>= 16;
}
while (i >= BBBB) {
result += 8;
i >>= 8;
}
while (i >= CCCC) {
result += 4;
i >>= 4;
}
while (i >= 2) {
result += 1;
i >>= 1;
}
return result;
}
```
版本三使用了 __builtin_clz 來改寫 ，他是 GCC 中的 extention，功能是會計算括號內數左邊0的個數
根據 [GCC 手冊](https://gcc.gnu.org/onlinedocs/gcc/Other-Builtins.html)
> Built-in Function: int __builtin_clz (unsigned int x)
Returns the number of leading 0-bits in x, starting at the most significant bit position. If x is 0, the result is undefined.
所以括號中應該要填 `v` ，但為了避免 `v` 為 0 時 `__builtin_clz` 沒定義的問題， DDDD 應該要填 `v|1` 這樣最小值就是 1 不是 0
## 測驗四
[Exponentially Weighted Moving Average](https://en.wikipedia.org/wiki/Moving_average#Exponential_moving_average) EWMA 是一種 moving average ，越新的樣本權重越大，隨著時間增加權重等比下降，但定義上永遠不會變0，所以和 simple moving average 不同，EWMA 是一種 IIR filter，然而得到新的點不用經過太多計算，因為可以用遞迴關係表示:
$$
S_t =
\begin{cases}
Y_0 \quad\;\quad\;\quad\;\quad\;\quad\;\quad\;\text{, }t=0\\
\alpha Y_t + (1-\alpha) S_{t-1} \quad \text{, }t>0\\
\end{cases}\
$$
其中 $0<\alpha \leq 1$ ，如同其他 ME ，它是個低通濾波器，可以觀察到若 $\alpha$ 值越大，過去的樣本影響越小，因此留下的高頻成分越多，若 $\alpha$ 為最大值 1 ，EWME 將不發揮作用
題目實作程式碼如下：
:::spoiler
```c
#include <assert.h>
/* Exponentially weighted moving average (EWMA) */
struct ewma {
unsigned long internal;
unsigned long factor;
unsigned long weight;
};
/* This section contains universal functions designed for computing the EWMA.
* It maintains a structure housing the EWMA parameters alongside a scaled
* internal representation of the average value, aimed at minimizing rounding
* errors. Initialization of the scaling factor and the exponential weight (also
* known as the decay rate) is achieved through the ewma_init(). Direct access
* to the structure is discouraged; instead, interaction should occur
* exclusively via the designated helper functions.
*/
/**
* ewma_init() - Initialize EWMA parameters
* @avg: Average structure
* @factor: Scaling factor for the internal representation of the value. The
* highest achievable average is determined by the formula
* ULONG_MAX / (factor * weight). To optimize performance, the scaling
* factor must be a power of 2.
* @weight: Exponential weight, also known as the decay rate, dictates the rate
* at which older values' influence diminishes. For efficiency, the weight
* must be set to a power of 2.
*
* Initialize the EWMA parameters for a given struct ewma @avg.
*/
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;
}
/**
* ewma_add() - Exponentially weighted moving average (EWMA)
* @avg: Average structure
* @val: Current value
*
* Add a sample to the average.
*/
struct ewma *ewma_add(struct ewma *avg, unsigned long val)
{
avg->internal = avg->internal
? (((avg->internal << EEEE) - avg->internal) +
(val << FFFF)) >> avg->weight
: (val << avg->factor);
return avg;
}
/**
* ewma_read() - Get average value
* @avg: Average structure
*
* Returns the average value held in @avg.
*/
unsigned long ewma_read(const struct ewma *avg)
{
return avg->internal >> avg->factor;
}
```
:::
一個 `ewma` 結構代表一個 EWMA 濾波器， `internal` 表示要輸出的樣本值， `factor` 是一個放大因子用來放大輸入樣本降低誤差， `weight` 則和 $\alpha$ 有關。
函式部分有 `ewma_init` ，用來初始化一個 EWMA 結構，它先檢查 `weight` 和 `factor` 是不是二的冪，然後將它們取 $lg$ ，方便後面用平移運算取代乘法，還有將 `internal` 設為 0 。
再來是 `ewma_add` 函式:
```clike=
struct ewma *ewma_add(struct ewma *avg, unsigned long val)
{
avg->internal = avg->internal
? (((avg->internal << EEEE) - avg->internal) +
(val << FFFF)) >> avg->weight
: (val << avg->factor);
return avg;
}
```
函式首先判斷 `internal` 是不是非 0，也就是在判斷是不是初始狀態，如果是就將 `val` 放大取 $lg$ 前的 `factor` 倍; 如果不是初始狀態則要進行前面提到的遞迴運算 $\alpha Y_t + (1-\alpha) S_{t-1}$ ，`(avg->internal << EEEE) - avg->internal` 可以對應到 $(1-\alpha) S_{t-1}$ 的部分，整段可以寫成(忽略 factor) $((Y_t*weight-Y_t)+Y_{t-1})/weight$ ，因此 EEEE 應為 avg->weight ，FFFF 應為 avg->factor。
從這邊也就能看的出來 `weight` 和 $\alpha$ 應為倒數關係。
## 測驗五
```clike=
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 > GGG) + 1;
}
```
函式一開始的 `x--` 是為了破壞 `x` 剛好是二的冪的情況，函式中會達到計算 $lg$ 並且捨棄小數的效果，最後 return 的地方會再加上一達成向上取整的功能，如果沒有 `x--` 二的冪的輸出就會多一。
再來的原理和測驗版本二類似，只是寫法更緊湊，需要好好思考; 由於 `x` 是 32 bits 的無號整數型，因此可以不使用迴圈逐步將檢查範圍覆蓋到八個 byte 就好， 依序檢查 x 是否大於 $2^{2^4}$, $2^{2^3}$, $2^{2^2}$,如果是就除掉，用 or operator 將結果依序更新到 `r` ，最後在 return 之前已經檢查了 $2^4 + 2^3 + 2^2 + 2 = 30$ 個 bits ，還剩兩個要檢查，於是我代入函式傳入時 x=1~4 推定 GGG 為 1 。
# [第 4 週測驗題](https://hackmd.io/@sysprog/linux2024-quiz4)
## 測驗一
[Popcount](https://en.wikichip.org/wiki/population_count) 用來計算二進位數值中有幾個 1 ，實作版本一如下:
```clike=
unsigned popcount_naive(unsigned v)
{
unsigned n = 0;
while (v)
v &= (v - 1), n = -(~n);
return n;
}
```
`v &= (v - 1)` 的效果是將最低位的 1 變成 0 ， `n = -(~n)` 的作用相當於 `n++` ，因為在二補數系統中負號是先取反再加一，所以 `－(~n)` = `~(~n)+1` = `n+1`， `n` 會在迴圈中計數直到 `v` 為 0 終止，但這樣程式不是常數時間，可能有資安疑慮，因此可以改進成版本二:
```clike=
unsigned popcount_branchless(unsigned v)
{
unsigned n;
n = (v >> 1) & 0x77777777;
v -= n;
n = (n >> 1) & 0x77777777;
v -= n;
n = (n >> 1) & 0x77777777;
v -= n;
v = (v + (v >> 4)) & 0x0F0F0F0F;
v *= 0x01010101;
return v >> 24;
}
```
這個方式是將每 4 個 bits 平行計算， `(v >> 1) & 0x77777777` 作用是先除以二，且為了保持每 4bits 的獨立性，將每四位第一位改成 0 。
考慮到計算 n bits popcount 公式如下:
$$
nbitspopcount(x) = x - \lfloor \frac{x}{2^1} \rfloor- \lfloor \frac{x}{2^2} \rfloor-...- \lfloor \frac{x}{2^n} \rfloor
$$
所以 4 bits 就是 $x - \lfloor \frac{x}{2^1} \rfloor- \lfloor \frac{x}{2^2} \rfloor- \lfloor \frac{x}{2^3} \rfloor$，第 4~9 行用遞迴的方式達成此式。
但這樣還不是最終結果，還需要把 `v` 每 4bits 分開再加起來，第 11 行會讓第奇數個 nibble 為原本的每兩個相加，第偶數個變 0 ，因為原本每個 nibble 最大是 4 所以不會溢位，然後再乘以 0x01010101 再右移 24 bits 寫成直式就能看出會是所有 nibble 相加的效果。
#### Hamming distance
兩數的 Hamming distance 可以理解為將一個數變成另一個數需要翻轉的位元個數，實作上可以先做 bitwise XOR 再用 popcount 來數:
```clike=
int hammingDistance(int x, int y)
{
return __builtin_popcount(x ^ y);
}
```
其中`__builtin_popcount` 是 GCC 內建函式。
考慮 [LeetCode 477. Total Hamming Distance](https://leetcode.com/problems/total-hamming-distance/) 的解法:
```clike=
int totalHammingDistance(int* nums, int numsSize)
{
int total = 0;;
for (int i = 0;i < numsSize;i++)
for (int j = 0; j < numsSize;j++)
total += __builtin_popcount(nums[i] ^ nums[j]);
return total >> AAAA;
}
```
題目要求加總陣列每兩個數的 Hamming distance ，for loop 算的每一組會重複配對一次，所以結果要除以二， AAAA 為 2 。
## 測驗二
若 $a = b(mod \quad m)$ , $c = d(mod \quad m)$ 則 $ac = bd(mod \quad m)$ , $a+c = b+d(mod \quad m)$ 。
若除數為 3 ，則套用上面第二個性質可以推導出
$$
2^k =
\begin{cases}
1(mod3) \quad\;\text{,k is even }\\
-1 \quad \text{, k is odd }\\
\end{cases}\
$$
考慮正整數 $n$ 以二進位表示 : $b_{n-1}b_{n-2}...b_{0}$， $n = b_{n-1}2^{n-1}+b_{n-2}2^{n-2}+...+b_02^0$ ，套用上面性質得到 $n = \Sigma ^{n-1}_{i=0}b_i(-1)^i(mod\quad3)$ ，求和的部分相當於奇數位和減偶數位和: `popcount(n & 0x55555555) - popcount(n & 0xAAAAAAAA)` ，$mod 3$ 的部分則是重複操作前面直到 n = 0~2 ，類似的方法也可以推廣到任何除數為 $2^k\pm1$ 的時候。
考慮 `tictactoe.c` 中 mod 7 的部分:
```clike=
static inline int mod7(uint32_t x)
{
x = (x >> CCCC) + (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);
}
```
不同於之前的討論，這裡是用 [Hacker's Delight](https://web.archive.org/web/20180517023231/http://www.hackersdelight.org/divcMore.pdf) 中 14 頁的方法，可以避免使用 popcount ，對於每個 $2^k\pm1$ 的除數，都可以從某些地方拆分成前後兩項再加起來而不影響求模結果，在這題由於後面是 & 0x7FFF 是 15 個 1，可以推測 CCCC 是 15。
```clike=
/* Determine if the tic-tac-toe board is in a winning state. */
static inline uint32_t is_win(uint32_t player_board)
{
return (player_board + BBBB) & 0x88888888;
}
```
如果想知道 BBBB 的答案就要先知道如何描述棋盤，它使用一個 32 bits 無號整數來描述棋盤狀態，經過比對 `move_masks` 可以發現，每四個 bits 描述一種連線，總共有八條剛好 32 bits，第一個 bit 為 0，第 2~4 分別代表那條線上三個位置有沒有值，所以 BBBB 是 0x11111111，因為這樣只要存在一條連線那條線的四個 bit 就會是 0111 + 0001 然後進位，接著 & 1000 看有沒有進位就完成檢查。
## 測驗三
AVL tree 為了極致的降低搜尋速度，限制任何左右子樹高度差不能超過一，因此常常需要透過旋轉 (LL, LR, RR, RL)來平衡子樹，過程會消耗許多資源，因此 [XTree](https://gist.github.com/jserv/7374b3e031ec3a81441b1b93b008c5d2) 被設計出來，它比較不會消耗資源去平衡，但樹的分布較不平衡，平衡能力和搜尋時間介於 red-black tree 和 AVL tree 之間。
接著來看移除節點的函式:
```clike
static void __xt_remove(struct xt_node **root, struct xt_node *del)
{
if (xt_right(del)) {
struct xt_node *least = xt_first(xt_right(del));
if (del == *root)
*root = least;
xt_replace_right(del, AAAA);
xt_update(root, BBBB);
return;
}
if (xt_left(del)) {
struct xt_node *most = xt_last(xt_left(del));
if (del == *root)
*root = most;
xt_replace_left(del, CCCC);
xt_update(root, DDDD);
return;
}
if (del == *root) {
*root = 0;
return;
}
/* empty node */
struct xt_node *parent = xt_parent(del);
if (xt_left(parent) == del)
xt_left(parent) = 0;
else
xt_right(parent) = 0;
xt_update(EEEE, FFFF);
}
```
根據註解，如果要刪除的節點有右子節點，則待刪除的節點會先替換為右子樹中遇到的第一個節點，因此 AAAA 為 `least` ，它在前面被設為 `xt_first` ，接著會對新接上節點的右子節點做 `xt_update` 來更新 `hint` ，並且決定要不要平衡，因此 BBBB 是 `xt_right(least)` ，同理 CCCC 為 `most` ， DDDD 為 `xt_left(most)`。
如果要刪除的節點是 leaf 的情況可以直接刪除，但就要檢查它的親代，確保平衡狀況，所以 EEEE 是 `root` ， FFFF 是 `parent` 。
Sign in with Wallet
Connect another wallet