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

## 作業表單問題
### 測驗 `4`
在 `x(1,3)` 使用情境會得到 12321 字串輸出，請補完 COND 敘述，使該程式運作符合預期。用最精簡的形式書寫程式碼，不要提交空白字元。

int x(int i, int k) { 
    return (i < k && putchar(i + '0') && COND || putchar(i + '0'));
}

根據現有的程式敘述，當 `i < k` 時會將 `i + '0'` 字元加入到回傳的字串中 (在字串的頭部加一次，在字串的尾部也加一次)，而當 `i >= k` 時則只會將 `i + '0'` 字元從頭部加入一次。

我根據 12321 此輸出字串推斷 COND 應該要呼叫遞迴直到 `i == k`，因此一開始我將 COND 設為 `x(i + 1, k)`，期望輸出和題目要求相等，但實際執行結果為 123，當我仔細檢查時才發現我遺漏了 `&&` 和 `||` 的 [short-circuit evaluation](https://en.wikipedia.org/wiki/Short-circuit_evaluation) 特性，當判斷式為 `a && b` 且 `a` 為假時，就不會執行 `b` 條件，而當判斷式為 `a || b` 且 `a` 為真時，就不會執行 `b` 條件。

當我使用 `x(i + 1, k)` 替換 COND 敘述時，會產生 `||` 前條件都為真的情形，因此不會執行 `||` 後的 `putchar()`，例如當 `x(2, 3)` 時 `i < k` 為真，`putchar(2 + '0')` 為真，`x(3, 3)` 也為真，導致不會執行後方的 `putchar(2 + '0')`，所以應該將 COND 替換為 `!x(i + 1, k)`，才能得到預期的輸出。

## 第 3 週測驗題
### 測驗 `1`
若要求的平方根為 N，透過 [digit-by-digit calculation](https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Digit-by-digit_calculation) 將 `N` 拆成 2 的冪相加。
$$
N^2 = (000b_0b_1...b_{n-1}b_n)^2
$$
其中 $b_0$ 為最高位的 1，可再將 N 表示為
$$
N^2 = (a_n+a_{n-1}+...a_0)^2
$$
其中 $a_i$ 為 2 的冪或 0，再依照 $(x+y)^2=x^2+2xy+y^2$ 的公式逐層拆解。
$\begin{split}
N^2& = (\displaystyle\sum_{i=0}^{n}a_i)^2 \\
&= a_n^2 + [2a_n+a_{n-1}]a_{n-1} +...[2\displaystyle\sum_{i=1}^{n}a_i+a_0]a_0 \\
&= a_n^2 + [2P_n+a_{n-1}]a_{n-1} +...+ [2P_1+a_0]a_0 \\
P_m &= a_n + a_{n1}+...+a_0 \\
&= P_{m+1} + a_m \\
P_m^2 &= P_{m+1}^2 + 2a_mP_{m+1} + a_m^2
\end{split}$

接著設 $Y_m=P_m^2-P_{m+1}^2$，則 $P_m^2=P_{m+1}^2+Y_m$，而 $X_m = N^2-P_m^2=X_{m+1}-Y_m$，再把 $Y_m$ 拆解為 $c_m+d_m$。
$$
c_m = 2P_{m+1}2^m，d_m = (2^m)^2
$$
$$
c_{m-1} = P_m2^m = (P_{m+1} + a_m)2^m = P_{m+1}2^m + a_m2^m = \dfrac{c_m} 2 + d_m
$$
而 $c_{-1} = P_02^0 = N$ 就是所求的 N。

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

此函式中的 `b` 對應到數學式的 $Y_m$，`z` 對應到 $c_m$，而 `m` 對應到 $d_m$。
### 測驗 `2`
為了減少運算成本，題目中想利用位元運算來代替 `/` 和 `%` 兩種運算元，因為除以 10 相等於乘 $\dfrac{1}{10}$，所以利用 2 的冪的倒數和來逼近。
$$
\dfrac{1}{128} + \dfrac{1}{32} + \dfrac{1}{16} = \dfrac{13}{128}\approx \dfrac{1}{10}
$$
而 $13 = 8 + 4 + 1 = 2^3 + 2^2 + 2^0$，因此用位元運算可得到 $\dfrac{13}{8}n$。

q = ((((n >> 3) + (n >> 1) + n) << 3)

接著將做位元運算時捨棄的值加回後再右移 7 (除以 128)。

d0 = q & 0b1;
d1 = q & 0b11;
d2 = q & 0b111;
q = ((((n >> 3) + (n >> 1) + n) << 3) + d0 + d1 + d2) >> 7;

而 $(4q + q) \times 2 = 10q$。

r = n - (((q << 2) + q) << 1);

而另一種算法先將 `q` 逼近到 0.8 in。

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

$$
x = \dfrac{3}{4}in
$$
$$
q= \dfrac{3}{64}in + \dfrac{3}{4}in = \dfrac{51}{64}in\approx0.8in
$$
將 0.8 右移 3 位可得到 0.1 ( `div` )。

*div = (q >> 3);

再計算 `mod = in - div * 10`。

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

### 測驗 `3`
#### 程式碼運作原理
版本二函式依序判斷目前的 `i` 是否大於或等於 65536 ( $2^{16}$ )、256 ( $2^{8}$ )、16 ( $2^{4}$ )、2 ( $2^{1}$ )，如果成立的話就將 `i` 進行右移並且增加 `result`，如此一來可找到最大的 `result` 使得 $2^{result} <= i$。

while (i >= 65536) {
    result += 16;
    i >>= 16;
}
while (i >= 256) {
    result += 8;
    i >>= 8;
}
while (i >= 16) {
    result += 4;
    i >>= 4;
}
while (i >= 2) {
    result += 1;
    i >>= 1;
}
return result;

版本三函式使用 GCC 內建的函式 `__builtin_clz` 計算傳入的無號整數從最高位開始連續 0 的數量，再使用 31 減 `__builtin_clz` 的結果，為了處理輸入為 0 的情況，會將 `__builtin_clz` 的輸入參數和 1 做 or。

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

參考 [linux/log2.h](https://github.com/torvalds/linux/blob/master/include/linux/log2.h) 其中 `31 - __builtin_clz(v)` 等價於 `fls(v) - 1`，`fls` 會回傳從最低位往最高位找到的最後一個 1。

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

### 測驗 `4`
[EWMA](https://en.wikipedia.org/wiki/Moving_average#Exponential_moving_average) 是一種用於平滑時間序列數據的統計手法，讓越久以前的資料所佔的權重越低。
#### 程式碼運作原理
##### `ewma_init`
初始化 EWMA 結構，設定 `factor` 和 `weight`，為了使用位元運算而讓這兩個參數的值為 2 的冪，`factor` 用於在有限的數值範圍內維持盡可能高的精確度，當 `factor` 越大時可以提高精確度，但同時也會降低 EWMA 能表示的最大值，而 `weight` 用於控制舊值對於新值的影響程度，當 `weight` 越大時會使平均值對於新值的影響越小。

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

##### `ewma_add`
向 EWMA 中添加一個新值，如果 `internal` 不為 0，則按照已有的累計平均值計算新的加權平均值，如果 `internal` 為 0，則直接根據 `factor` 調整新值。

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

##### `ewma_read`
返回 EWMA 當前的加權平均值，通過對 `internal` 依 `factor` 反向縮放得到。

return avg->internal >> avg->factor;

### 測驗 `5`
#### 程式碼運作原理
設定 `r` 和 `shift` 紀錄各步驟的位移量和最終的對數結果，接著將傳入的 `x` 值減 1，這是為了處理當 `x` 為 2 的冪時的情況，因為此函式回傳的是 $\lceil\log_2(x)\rceil$，如果不先減 1，會導致回傳的結果比預期要多 1。

再檢查 `x` 是否大於 `0xFFFF`，如果為真說明 `x` 的對數最小為 16，將此資訊紀錄在 `r` 中，並將 `x` 右移 16 位。

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

再檢查當前 `x` 是否大於 `0xFF`，如果為真說明當前 `x` 的對數最小為 8，將此資訊紀錄在 `shift` 中，並將 `x` 右移 8 位，並將 `shift` 的值合併到 `r` 中，累計目前為止的對數結果。

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

根據相同的邏輯依序判斷當前 `x` 是否大於 `0xF` 和 `0x3`，設定 `shift` 和通過 `r` 累計結果，最後將判斷 `x` 是否大於 `0x1` 的步驟合併在最後一行中，並統計以上步驟中的對數結果，再將一開始減的 1 補上。

return (r | shift | x > 1) + 1;

#### 改進程式碼使其得以處理 `x = 0` 的狀況，並仍是 branchless
我新增了一個變數紀錄目前傳入的 `x` 值是否為 0，若 `x` 為 0 則先讓 `x` 值加 1 來抵銷 `x--`，使 `x` 不會產生 underflow，並且在最後不會再將回傳值加 1。

  uint32_t r, shift;
+ int is_zero = !x; 
+ x += is_zero; 
  x--;
    
- return (r | shift | (x > 1)) + 1;
+ return (r | shift | (x > 1)) + !is_zero; 

## 第 4 週測驗題
### 測驗 `1`
[population count](https://en.wikichip.org/wiki/population_count) 是用來計算一個二進位表示的數值中有幾個位元為 1。
#### 程式碼運作原理
`popcount_branchless` 此函式將 32 位元的 `n` 拆分成 4 個位元一組來計算為 1 的位元數。
在觀察程式碼時，一開始不理解為什麼將 `v` 右移一位( 除以 2 )後要再 `& 0x77777777`，後來發現這是為了讓每一組右移完的結果不會被其他組的結果所干擾，所以把每一組中最高位的位元設為 0。

1110 | 1011 | 0111 | 0101 //  v
0111 | 0101 | 1011 | 1010 //  v >> 1
0111 | 0101 | 0011 | 0010 // (v >> 1) & 0x77777777

以上範例可看出 `& 0x77777777` 後的結果才是以 4 個位元為一組右移一位的結果。

[LeetCode 477. Total Hamming Distance](https://leetcode.com/problems/total-hamming-distance/) 參考函式使用 `__builtin_popcount` 來計算一整數陣列中所有兩兩一組的 Hamming Distance 總和，此方法利用巢狀迴圈來計算兩兩一組的 Hamming Distance，但因為會重複計算一次，所以最後在回傳時要將結果除以 2。

for (int i = 0;i < numsSize;i++)
    for (int j = 0; j < numsSize;j++)
        total += __builtin_popcount(nums[i] ^ nums[j]); 
return total >> 1;

#### 針對 LeetCode 477. Total Hamming Distance 撰寫出更高效的程式碼
測驗上範例的時間複雜度為 $O(n^2)$，可將陣列中的每一個數值用二進位觀察，並計算每一個位元有幾個 1，假設第 0 個位元總共有 $x_0$ 個 1，此位元的 Hamming Distance 為 $x_0 \times (numsSize - x_0)$，依序從第 0 個位元計算到第 31 個位元的總和即為此題的解。

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