# 2024q1 Homework4 (quiz3+4)
contributed by < [MathewSu-001](https://github.com/MathewSu-001) >
## 第一周測驗1
### 程式碼理解
若要求得平方根 `N` ，以 2的冪相加的表示式為
$$N^2 =(a_n+a_{n-1}+a_{n-2}+....+a_0)^2$$
根據 [測驗一](https://hackmd.io/@sysprog/linux2024-quiz3) 以及 [Digit-by-digit calculation](https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Digit-by-digit_calculation) 所提及的定義，$P_m=\displaystyle\sum_{i=0}^{n-m}a_{n-i}$ 是迄今為止找到的近似平方根，則 $P_0=a_n+a_{n-1}+...+a_0$ 即為所求平方根 $N$。
所以我們就是要去計算出所有 $a_{n-i}$ 的總和，從 $a_n$ 開始加到 $a_0$ ，而求得 $a_m$ 只要試驗 $P_m^2 \leq N^2$ 是否成立
\begin{cases}
P_m = P_{m+1} + 2^m, & \text{if $P_m^2 \leq N^2$}\\
P_m = P_{m+1}, & \text{otherwise}
\end{cases}
考量到運算成本問題，一種想法是從上一輪的差值 $X_{m+1}$ 減去 $Y_{m}$ 得到 $X_{m}$，
- $X_m = N^2 -P_m^2 = (N^2 -P_{m+1}^2)-(P_{m}^2 -P_{m+1}^2)=X_{m+1}-Y_m$
- $Y_m=P_m^2 -P_{m+1}^2=(P_{m+1}+a_m)^2-P_{m+1}^2=2P_{m+1}a_m+a_m^2$
當 $X_{n}=0$ 時,找到精確的平方根；如果不是，則 $P_m$ 給出平方根的合適近似值，其中 $X_{n}$ 是近似誤差。
進一步調整，將 $Y_m$ 拆成 $c_m, d_m$
- $c_m = 2P_{m+1}a_m = P_{m+1}2^{m+1}(\text{if } a_m=2^m)$
- $d_m = a_m^2=(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_{m-1}, d_{m-1}$
- $c_{m-1}=P_m2^m=\begin{cases}
c_m/2+d_m & \text{if }a_m=2^m \\
c_m/2 & \text{if }a_m=0
\end{cases}$
- $d_{m-1}=a_{m-1}^2=(2^{m-1})^2=\dfrac{d_m}{4}$
所以可以知道當 $c_{-1}$ 時，得到值 $P_0$ 就是平方根 N 所求。
結合上述方法撰寫演算法來尋求 $P_0=a_n+a_{n-1}+...+a_0$ ，初始化條件
- $X_{n+1}=N^2$
- $c_n=2P_{n+1}a_n$ 且 $a_n$ 已是已知最大，所以 $P_{n+1} = 0 \to c_n=0$
- $d_n=a_n^2=(2^n)^2=2^{2n}$，所以可以設 $d_n=m$ 得到以下程式碼:
```c
int m = 1UL << ((31 - __builtin_clz(x)) & ~1UL)
```
上述程式碼可以得到 x 最近的偶數符合 $2^{2n}=4^n$。
```c
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;
```
在以上程式碼中，$z=c_n ,\ m=d_n,\ b=Y_n,\ x=X_n$。所以依照上述的演算法，不斷的迴圈直到 $m=d_n=0$ 的時候，那所求得的 z 即為 $c_{-1}=P_0$。
### 取代 `__builtin_clz`
使用 linux 提供的 [__fls](https://github.com/torvalds/linux/blob/master/tools/include/asm-generic/bitops/__fls.h) 來進行修改。因為他回傳的 return 值是找到 x 最高位數，所以引用方法如下:
```diff
int z = 0;
-for (int m = 1UL << ( (31 - __builtin_clz(x)) & ~1UL); m; m >>= 2) {
+for (int m = 1UL << ( __fls(x) & ~1UL); m; m >>= 2) {
int b = z + m;
z >>= 1;
if (x >= b)
x -= b, z += m;
}
```
> [58557c9](https://github.com/MathewSu-001/linux_hw4/commit/58557c9b1d4df0dcb7a31bf914c3fe75af9e0170)
## 第一周測驗2
### 程式碼解析
## 第一周測驗3
### 程式碼解析
從版本一開始看起，
```c
int ilog2(int i)
{
int log = -1;
while (i) {
i >>= 1;
log++;
}
return log;
}
```
可以知道找尋 $log_2$ 的做法就是取得最高位元的次方是多少。
所以在版本二中，
```c
static size_t ilog2(size_t i)
{
size_t result = 0;
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;
}
```
首先就是先看 i 的值有沒有超過 $2^{16}$(32 位元的一半)，假如有的話將 result 加上 16，並且數值右移 16 位，然後再依次砍半到 $2^8、2^4、2^1$ 看數值有無超過，這相當於找 fls。
所以版本三中，
```c
int ilog32(size_t v)
{
return (31 - __builtin_clz(v | 1));
}
```
因為 `__builtin_clz`的功用就是回傳 v 中從最高有效位開始的前導零位數，且如果 v 為 0，則結果未定義。所以加上 `| 1` 就是定義如果 v 為零時，回傳 0。
從測驗一中，我學到 fls 的用途就是找到 v 最高位數，與 `31 - __builtin_clz(v | 1)` 作用相同，因此我有進行嘗試並且證實最後結果相同。
> [a1de600](https://github.com/MathewSu-001/linux_hw4/commit/a1de6005a63199a29e3ab48d15d378846c67287f)
## 第一周測驗4
### 程式碼解析
透過[Exponentially Weighted Moving Average](https://en.wikipedia.org/wiki/Moving_average#Exponential_moving_average)(EWMA; 指數加權移動平均)的數學定義:
$$S_t = \left\{
\begin{array}{l}
Y_0& t = 0 \\
\alpha Y_t + (1 - \alpha)\cdot S_{t-1}& t > 0 \\
\end{array}
\right.$$
可以知道程式碼的 `ewma_add` 的計算方式
```c
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;
}
```
在變數上，`val` 為新的一筆資料；`avg` 中的 `internal` 為最後取的平均值，對應到公式中的 $S_t$ ; `weight` 為權重值，對應到公式中的 $\alpha$；最後 `factor` 與 $Y_t$ 的計算有關係，為縮放因子。
當 `avg->internal` 為零時，
$$S_t = Y_0,\ \ t = 0$$
當 `avg->internal` 不為零時，
$$S_t = \alpha Y_t + (1 - \alpha)\cdot S_{t-1}=[Y_t + (\frac{1}{\alpha} -1)\cdot S_{t-1}]\cdot\alpha$$
這邊需要注意的是，$\alpha$ 介在 0 ～ 1 之間，但這邊的 `weight` 為大於零的 2 的冪，所以在作法上，原本的除法要變成乘法( bitewise operation 做右移)，乘法變成除法( bitewise operation 做左移)。
## 第一周測驗5
### 程式碼解析
在一開始我先將程式做改寫，想要知道 `x--` 的用途:
```diff
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;
+ return (r | shift | x > 1) ;
}
```
就發現如果做出這樣的改動的話，程式碼功能就會跟測驗三的 `ilog2` 相同。因此先做 `x--` ，然後 return 值時再加一的用途，就是利用 `ilog2` 會以 2 的冪做為一個區段，原本 $2^k\sim2^{k+1}-1$ return 的值為 k ，但改動後會變成 $2^k+1\sim2^{k+1}$ return 的值為 k+1，達到 [ceil](https://man7.org/linux/man-pages/man3/ceil.3.html) 效果。
### 處理 `x = 0` 的狀況
當 `x = 0` 時，做 `x--` 會取得 `x = 0xFFFFFFFF`，而導致最後回傳的值是不正確的，所以在做法上參考了 [vax-r](https://hackmd.io/@vax-r/linux2024-homework4#%E6%B8%AC%E9%A9%97%E4%BA%94) 同學的做法，只有在 x > 0 才進行減法。
```diff
int ceil_ilog2(uint32_t x)
{
uint32_t r, shift;
- x--;
+ x -= !!x;
```
> [20f56ef](https://github.com/MathewSu-001/linux_hw4/commit/20f56ef20110a67e15f6ee8a0ec2141436f8cab0)
## 第二周測驗1
### 程式碼解析
Sign in with Wallet
Connect another wallet