# 2024q1 Homework4 (quiz3+4)
contributed by < [`padaray`](https://github.com/padaray) >

## quiz3

### 測驗 `1`
#### 計算平方根
**版本一：**
對輸入的參數 `N` 以 2 為底數取 log，意義是當我們將 `N` 取 log2 時，可以得到 most significant bit (msb)。
對 1 左移 msb 個位元後存入變數 `a` ，意義是找出一個必定大於 `N` 的平方根的數，利用此數來逼近找到 `N` 的平方根。
最後使用 while 迴圈，計算 `a + result` 的平方是否大於 `N`，若小於 `N` 則將此數存入 `result` 變數，對 `a` 右移 1 bit，繼續計算，直到 `a` = 0。

</br>

**版本二：**
版本二的差別是不使用 `math.h` 的 `log2` 函式，直接對參數 `N` 進行右移，右移至參數為 0 時，則代表找到 most significant bit。此方法將程式的可移植性增強，並且不影響原本的函式原型 (prototype) 和精度 (precision)

- int msb = (int) log2(N);
+ int msb = 0;
+ int n = N;
+ while (n > 1) {
+     n >>= 1;
+     msb++;
+ }

</br>

**版本三：**

完整程式碼(無遮蔽處)

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

`__builtin_clz(x)` ：是 gcc 內建處理二進位的其中一個函式巨集，此函式會回傳 counting leading zero，幫助我們找到 most significant bit。

在 for 迴圈中左移 `31 - __builtin_clz(x)` bit 數得到 msb，`& ~1UL` 可以確保最後一個 bit 為 0，理由是變數 `m` 必須是偶數對應到 $d_n = (2^n)^2$。而 m 右移兩個 bit 是對應到 $d_{m-1} = d_m/4$

</br>

### 測驗 `2`
#### 針對正整數在相鄰敘述進行 mod 10 和 div 10 操作，如何減少運算成本？

除法原理：$f=g \times Q + r$
* $f$: 被除數
* $g$: 除數
* $Q$: 商
* $r$: 餘數

程式碼如下：

carry = tmp / 10;  // carry 代表商(Q)
tmp = tmp - carry * 10;  // tmp 代表餘數(r)

如果我們想要使用 bitwise operation 來實作除法，如果被除數包含了除了 2 以外的因數，就會造成結果不精確，所以使用以下方式。
假設除數為 10，我們必須找到 $\dfrac{1}{10} \approx \dfrac{a}{2^N}$，最後發現是 $\approx \dfrac{13}{128}$
程式碼如下：

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

`(((tmp >> 3) + (tmp >> 1) + tmp) << 3)` 使用 bitwise operation 取得一個分子有 13 的分數，因此用 $\dfrac{tmp}{8} + \dfrac{tmp}{2} + tmp = \dfrac{13tmp}{8}$，之後再左移三位消掉分母 8。
`d0`、`d1`、`d2` 是為了將 `tmp` 右移時損失的 bit 數加回，加回後再右移 7 位也就是除以 128 得到商數。
`(((q << 2) + q) << 1)` 相當於 $(q \times 4 + q) \times 2$，也就是最一開始程式碼的 `carry * 10`。

最終程式碼可包裝為以下：

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

</br>

### 測驗 `3`
#### 考慮 ilog2 可計算以 2 為底的對數，且其輸入和輸出皆為整數。

**版本一：**
初始化變數 `log` 為負一，直接對參數 `i` 進行右移，每次右移將 `log++`，右移至參數 `i` 為 0 時

</br>

**版本二：**
完整程式碼(無遮蔽處)

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

為了避免數字過大每次位移一個 bit 太慢，例如：65536 需要位移 16 個 bit，因此加入條件判斷，若大於特定 2 的冪次方，直接移動該次方的位元數。

</br>

**版本三：**

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

</br>

### 測驗 `4`
#### 解釋 EWMA 的實作
Exponentially Weighted Moving Average (EWMA) 是一種統計方法，透過對每筆資料給予不同的權重計算平均，通常越新的資料權重就越高，越舊的資料權重越小，好處是可以適應新的資料變化，又保留一定的歷史資料影響。

數學定義：
$S_t = \begin{cases}
Y_0, & t=0 \\
\alpha Y_t + (1-\alpha)\cdot S_{t-1}, & t > 0
\end{cases}$

$\alpha$ : 舊資料加權降低的程度，越高代表降低越快。介於 0 ~ 1 之間
$Y_t$ : 時間 t 時的資料
$S_t$ : 時間 t 時算出的 EWMA

以下為程式碼實作：

 * 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 << avg->weight) - avg->internal) +
                           (val << avg->factor)) >> avg->weight
                        : (val << avg->factor);
    return avg;
}

`internal` 代表目前計算的 ewma (數學式中的 $S_t$) ， `(avg->internal << avg->weight)` 代表數學式中的 $\alpha = \dfrac{1}{2^{weight}}$ (**但不知道原因為何**)，因此 $1- \alpha = \dfrac{2^{weight} - 1}{2^{weight}}$。
當 `internal` 為 0 時，`internal` 代入 `(val << avg->factor)` 寫成公式為 $2^{factor} \cdot val$。
當 `internal` 不為 0 時，`internal` 代入 `(((avg->internal << avg->weight) - avg->internal) +
(val << avg->factor)) >> avg->weight` 寫成公式為 $\dfrac{2^{weight}-1}{2^{weight}} \cdot internal + \dfrac {2^{factor}}{2^{weight}} \cdot val$，轉換一下可以得到 $(1-\alpha) \cdot internal + {2^{factor}} \cdot \alpha \cdot val$

</br>

### 測驗 `5`
#### 延伸測驗 `3` 使用下方程式碼計算 $\lceil \log_2 x\rceil$

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` 是否大於 0xFFFF，意義是確定 `x` 的 msb 是否位於高位的 16 bit，若為真，將 `x` 右移 16 bit。
第七、八、九行判斷參數 `x` 在這剩下的 16 bit 中，msb 是否位於高位的 8 bit，若為真，將 `x` 右移 8 bit，將右移的位數紀錄在 `r` 中。
第十行到第十四行依此類推。
第四、十五行將 `x` 先 -1 後 + 1 回傳是為了處理 2 的冪次方，若不處理會讓 2 的冪次方多計算一次，例：$2^7$ 經過無條件進位法，會變成 8，因此要先減 1

</br></br>

## quiz4

### 測驗 `1`
#### [popcount(population count)](https://en.wikichip.org/wiki/population_count)