# 2024q1 Homework 4 (quiz3+4)
**contributed by <`MiohitoKiri5474`>**
## 第三週
### 測驗一
根據 [digit-by-digit calculation](https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Digit-by-digit_calculation) 原理,假設欲求平方根為 $N$,那我們可以把 $N^2$ 化為 2 的冪總和:
$N^2 = (a_0+a_1+ \dots + a_n)^2, a_m = 2^m(or \ 2^m = 0 ), m \in \{0, 1, \dots, n \}$
已知 $( a + b ) ^ 2 = a^2 + 2ab + b^2$,則有以下遞迴關係:
$$N ^ 2 = ((a_0 + a_1 + \dots + a_{n - 1}) + a_n)^2 = {a_n}^2 +2a_n \cdot \sum^{n - 1}_{k = 0}a_k+(a_0 + a_1 + \dots + a_{n - 1})^2$$
$$N ^ 2= {a_n}^2 +(2P_n + a_{n - 1} + \dots + (2P_1 + a_0)a_0$$
同時令平方根逼近值 $P_m = P_m + a_{m - 1}$,因此對於每個 $m \in \{0, 1, \dots, n \}$,於 $P_{m - 1}$ 加上 $a_m$,由於是 $2$ 的冪,因此只會有以下兩種可能:
$$P_{m - 1} = \begin{cases} P_m + a_{m - 1}, & \text{if} \ {P_m}^2 < N^2 \\ P_m , & \text{otherwise} \end{cases}$$
但我們的程式碼不會用這種比較直白的方式判斷是否加上 $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$,令人有些費解,再經過一陣子思考後,發現只是單純將加法的兩項拆開來看,也就是說:
$$\begin{cases} Y_m = c_m + d_m \\ c_m = P_{m+1} \cdot 2^{m + 1} \\ d_m = (a_m)^2 \end{cases}$$
而每次迴圈更新則是:
$c_{m - 1} = P_m \cdot 2^m = (P_{m + 1} + a_m) \cdot 2^m = \begin{cases} \frac{c_m}{2} + d_m, if \ a_m = 2^m \\ \frac{c_m}{2}, if \ a_m = 0 \end{cases}$
$d_{m - 1} = \frac{d_m}{4}$
這個計算會不斷進行下去,直到 $c_{-1} = P_0 = N, d_0 = 0$ 的時候停止。
### 測驗二
為了避免使用 / 和 % 來計算除以 $10$ 的商和餘,因此我們要找到 $\frac{1}{10}$ 的二進位表示法(除以 $10$ 可以看作乘以 $\frac{1}{10}$、同時得出除以 $10$ 的結果之後就能得到餘數),但 $10$ 不符合 $2^k \pm 1$ 的格式,因此無法以二進位或 digit recurssion 來表示,只能以一個非常接近 $\frac{1}{10}$ 的數字 $\frac{a}{2^n}$ 來表示 $10$ 的倒數。
題目這邊挑選 $2^N = 128, a = 13, \frac{128}{13} \approx 9.84,因此我們可以把 `x / 10` 的計算轉換成 `x * 13 / 128`,再進一步拆解成 `x * 13 / 8 / 16`。
其中 `x * 13 / 8` 可以透過 `(x >> 3) + (x >> 1) + x` 得到,為了將其還原成 `x * 13` 會再進行一次左移,不過這樣會遺失最低三位元的值,因此在操作前先將最低的三為元的值先存起來、操作完成後再加回去,得到 `x * 13` 之後再左移七位(除以 $128$)。
> 這邊其實也可以直接再 `(x >> 3) + (x >> 1) + x` 後直接右移四,計算結果是一樣的。
> 反正丟失的三位最後也是會被右移七吃掉,在這個情況下可以不用考慮儲存起來。
```c
unsigned d0, d1, d2, q;
d0 = x & 0b1;
d1 = x & 0b11;
d3 = x & 0b111;
q = ((((x >> 3 ) + (x >> 1) + x) << 3) + d0 + d1 + d2) >> 7;
q = ((x >> 3) + (x >> 1) + x) >> 4;
```
而餘數再使用餘數定理而得:
```c
unsigned r = x - (((q << 2) + q) << 1);
```
完整的函數可以寫成這樣:
```c
void calc_div_mod_10 ( uint32_t in, uint32_t *div, uint32_t *mod ) {
unsigned d0, d1, d2;
d0 = in & 0b1;
d1 = in & 0b11;
d3 = in & 0b111;
*div = ((((in >> 3 ) + (in >> 1) + in) << 3) + d0 + d1 + d2) >> 7;
*mod = in - (((*div << 2) + *div) << 1);
}
```
最後再來看一下包裝的結果:
```c
#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));
}
```
看似跟原先我們計算出來的結果很不相同,好像沒辦法很順利的轉換;那麼我們先來分析一下。
我們檢視前面的計算,發現在函數第一行將 $\frac{3}{4} \cdot in$ 計算出來了,後面再加上一個 $\frac{3}{64} \cdot in$ 也就是 $\frac{51}{64} \cdot in$,而 $\frac{51}{64} \approx 0.797$。
後面再加上 `q >> 8` 也就是 $(\frac{51}{64} \cdot in) \cdot \frac{257}{256} = \frac{13107}{16384} \cdot in \approx 0.7999$,經過額外三次的操作後 $q = \frac{219902325555}{274877906944} \approx 0.8 \cdot in$。
也就是說,到最後要真正計算出 `div` 之前,前面已經將 $in \cdot \frac{13}{16}$ 的部分計算完了,所以最後只要再將最後的 $8$ 除掉就好了。
而計算 `mod` 時,使用 `q & ~0x7` 用意是將最後三個位元的值清除,作用和 `*div << 3` (`*div * 8`)相同,從這邊可以推敲出 `*div << DDDD` 就是要將剩下的值補上。
### 測驗三
原先的作法有點類似二分搜,從高位開始檢查,如果輸入大於 $2^16$ 會將輸入的數字左移 16 位、並將回傳值加上 16;再來檢查是否大於 $2^8$ 依此類推。
另外我們可以發現 `ilog2(n)` 就是在計算 n 從 msb 數來第一個 set bit 的位置,因此就能使用 `__builtin_clz(n)` 先計算 `n` 的 leading zeros,接著用 31 減掉它也就是結果。
### 測驗四
EMWA 公式為:
$$S_t = \begin{cases} \alpha \cdot Y_t + ( 1 - \alpha ) \cdot S_{t - 1}, & t > 0 \\ Y_0, & t = 0 \end{cases}$$
從 `ewma_init` 函數中可以發現,一開始初始化的時候會將 `weight`, `factor` 先對數化,所以平移等量的次數可以視為乘除。
而 `ewma_add` 則是對應到 EWMA 的公式,其中比較特別的是每次資料點會先乘上(左移) `factor` 倍才帶入公式中。
`EEEE` 為 `avg->weight`、`FFFF` 為 `avg->factor`。
### 測驗五
經過觀察發現,測驗三的每個 while 迴圈其實只會執行一次,也就是說我們可以改成 if 的形式實作。
> 舉例來說,如果 $x \ge 2^8$ 這一項會執行兩次,代表 $x \ge 2^8\cdot 2^8 = 2^{16}$,這種情況應該會在前一項 $x \ge 2^{16}$ 時就被篩掉了。
而改為 if 的形式,則代表說我們可以改寫成不會造成分支的寫法,也就是原先題目中給的樣式,並將位移量保留在變數 `r` 中。
## 第四週
### 測驗一
總之,要計算出所有組合的 Hamming Distance 和。
```c
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 >> 1;
}
```
觀察一下會發現,所有組合都會被重複算到(舉例來說,`nums[1], nums[2]` 和 `nums[2], nums[1]` 都會被分別計算一次),因此最後要將答案除以二,也就是右移一。
另外我們可以把 `j` 的起始值設成 `i+1`,這樣就能夠避免重複計算的問題了。
#### 避免使用 popcount 的改進
即便我們從根本上減少一半的計算數量,但還是沒有將複雜度下降(一樣為 $O(N^2)$,因此需要尋找其他作法。
我們換個角度思考,如果我們改成計算出 **某位數為一的數字**,同時也能算出 **某位數不為一的數字**(用全部減掉該位數為一的數量),將這兩個數字相乘,就會是該位數的 Hamming Distance,並將所有位數的 Hamming Distance 加總在一起即可。
所以程式碼可以改寫成這樣:
```c
int HammingDistantLinear(int* nums, int numsSize)
{
int res = 0;
for (int i = 0, cnt = 0, it = 1; i < 32; i++, cnt = 0, it <<= 1) {
for (int j = 0; j < numsSize; j++)
if (nums[j] & it)
cnt++;
res += (numsSize - cnt) * cnt;
}
return res;
}
```
### 測驗二
tictactoe 的棋局一共會有九種可能的選擇,因此我們可以將棋盤以二進位的形式紀錄(對於該玩家而言,哪一位置有他的記號),並且將整個棋盤一共八種路線(三條橫、三條竪、兩條對角線)分別紀錄下來,這樣只需要 24 位元就能夠很詳盡的紀錄完所有資訊,不需要額外做判斷(是否能連成線等等)。
但實際上實作過後才會發現,這樣還是不便於紀錄,那麼我們將其中一條連線的狀態取出觀察,會發現三個位元會是以 $(111)_2$ 的方式呈現,那麼我們便可以為每一條路線新增一位,這樣將這條路線的值 +1 的時候就會成為 $(1000)_2$ 的狀態,可以更方便的做運算;且即便新增一位,總位元數也還在 `uint32_t` 的範圍內,同時每四個位元為一組可以更輕易的轉換成 hex 編碼。
接著往下看,在 `mod7` 函數中,其中 $(24924925)_{16} \approx \frac{2^{32}}{7}$,從 [Hacker's Delight](https://doc.lagout.org/security/Hackers%20Delight.pdf) 中知道可以利用 $8 \equiv 1 \pmod{7}$ 此數學關係是,可以推導出 $n \pmod{7} = popcount(n)$。
例如書中的其中一段程式碼:
```c
int remu7(unsigned n) {
static char table [75] = {0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4};
n = (n >> 15) + (n & 0x7FFF);
n = (n >> 9) + (n & 0x001FF);
n = (n >> 6) + (n & 0x0003F);
return table[n];
}
```
可以發現是利用 $n \pmod{7} = popcount(n)$ 不斷的進行 digit sum,每一步都在重複求 modulo 7。
回來觀察題目中的程式碼,可以發現該程式碼中另外運用了 $n \pmod{7} \equiv \lfloor \frac{8}{7} n \rfloor \pmod{8}$,而我們可以利用 $\lfloor \frac{2^{32}}{7} \rfloor$ 再右移 29 位已取得 $\lfloor \frac{8}{7} n \rfloor \pmod{8}$ 的值,雖然取 floor 會有誤差(例如題目中的程式碼以 `0x24924925` 作為 $\lceil \frac{2^{32}}{7} \rceil$ 的值),不過乘上這個數字後許多高位位元會 overflow 被遺棄掉;同時按照操作將最後的數字右移 29 位元後即是想要的答案。
### 測驗三
按照題目要求補完程式碼,推測是要將目前的節點刪除,並且如果右子樹並非為空則從右子樹中找出適當的節點;如果右子樹為空則從左子樹中尋找;如果兩者皆為空則視此子樹為空。
而按照 AVL 的性質,該節點的右子樹中的所有節點的值、皆要比該節點的值還大,因此從右子樹中需要找出最小的值。
同樣按照 AVL 的性質,該節點中的左子樹中的所有節點的值、皆要比該節點的值還要小,因此從左子樹中需要找出最大的值。
並且將當前節點以上述節點替代時,需要重新平衡該子樹。
```c
if (xt_right(del)) {
struct xt_node *least = xt_first(xt_right(del));
if (del == *root) *root = least;
xt_replace_right(del, least);
xt_update(root, xt_right(least));
return;
}
if (xt_left(del)) {
struct xt_node *most = xt_last(xt_left(del));
if (del == *root) *root = most;
xt_replace_left(del, most);
xt_update(root, xt_left(most));
return;
}
```
```c
/* 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(root, parent);
```