# 2024q1 Homework4 (quiz3+4)
contributed by < [`jimmylu890303`](https://github.com/jimmylu890303)>
## 第三週測驗題
> [題目](https://hackmd.io/@sysprog/linux2024-quiz3)
### 測驗一
#### 版本一

#include <math.h>
int i_sqrt(int N)
{
    int msb = (int) log2(N);
    int a = 1 << msb;
    int result = 0;
    while (a != 0) {
        if ((result + a) * (result + a) <= N)
            result += a;
        a >>= 1;
    }
    return result;
}

在版本一中，會透過 `log2(N)` 去尋找 N 的最高位元，再從高位逐步去尋找滿足 `(result + a) * (result + a) <= N` 的位元。

以輸入 N = 36 為例， msb = 5 且 a = 100000

| a | result | (result + a) * (result + a) <= N |
| ---------- |:------:|:--------------------------------:|
| 100000(32) | 0 | X |
| 10000(16) | 0 | X |
| 1000(8) | 0 | x |
| 100(4) | 0 | ==O== |
| 10(2) | ==4== | ==O== |
| 1(1) | ==6== | X |

所以 result 的值為 6 ，在此版本中使用到 log2 ，所以需要在編譯時多新增編譯選項 `-lm`

#### 版本二
在版本二中，去改善版本一使其不依賴 log2 並保持原本的函式原型 (prototype) 和精度 (precision)，而改善的作法就是使用 shift 的方式去找出最高位的位元。

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

#### 版本三
在版本三中使用 `__builtin_clz(x)` 的巨集來找出最高位在哪個位元算出 leading zero 個數。

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

| x | m | z | b | z(after shift) | x >= b | x-=b | z+=m |
| --- |:---------:|:---:| --- |:--------------:|:------:| ---- |:----:|
| 36 | 10000(16) | 0 | 16 | 0 | ==O== | 20 | 16 |
| 20 | 100(4) | 16 | 20 | 8 | ==O== | 0 | 12 |
| 0 | 1(1) | 12 | 13 | 6 | X | 0 | 6 |

### 測驗二
針對正整數在相鄰敘述進行 mod 10 和 div 10 操作，如何減少運算成本？
採用 bitwise operation 來實作除法

先考慮除以 10 的情況，除以 10 相當於乘以 $\dfrac{1}{10}$
但由於 $\dfrac{1}{10}$ 無法以二進位的方式去表示，所以只能透過一個近似的值去替代
而 $\dfrac{1}{10} \approx \dfrac{13}{128}$ ，因為 $\dfrac{128}{13} \approx 9.84$
所以 x 除以 10 相當於 $x \times \dfrac{13}{128}$
更進一步拆解上述的式子 $x \times 13 \div 8 \div 16$

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

考慮以下式子，相當於 $\dfrac{tmp}{8} + \dfrac{tmp}{2} + tmp = \dfrac{13}{8} \times tmp$

((tmp >> 3) + (tmp >> 1) + tmp) 

由於往右作 shift 3個位元會造成最右邊三個位元資訊損失，所以在往左 shift 後要將最右邊三個位元加回來。所以以下操作值會等於 $13 \times tmp$

((((tmp >> 3) + (tmp >> 1) + tmp) << 3) + d0 + d1 + d2)

所以最後往右 shift 7個位元相當於除上128，這樣就會得到 $tmp \div 10$ 的結果

q = ((((tmp >> 3) + (tmp >> 1) + tmp) << 3) + d0 + d1 + d2) >> 7;

而 tmp 取 10 的餘數則利用餘數定理 $r = f - g \cdot Q$

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

q 為商數，則以下操作為 $q \times 10$，即可算出餘數 r

(((q << 2) + q) << 1)

以上的概念可以包裝成以下程式碼

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

在以下操作中， $x = \dfrac{3}{4} \times in$

uint32_t x = (in | 1) - (in >> 2);

$q = \dfrac{3}{4} \times in + \dfrac{3}{64} \times in = \dfrac{51}{64} \times in \approx 0.79 \times in$

uint32_t q = (x >> 4) + x;

而這邊的操作是將 0.79 逼近至 0.8

x = q;
q = (q >> 8) + x;
q = (q >> 8) + x;
q = (q >> 8) + x;
q = (q >> 8) + x;

$div = \dfrac{8}{10} \times in \times \dfrac{1}{8} = \dfrac{1}{10} \times in$

*div = (q >> 3);

考慮以下操作
`q >> 3` 為商數
而 `q & ~0x7` 將最低三位清零，相當於商數往左 shift 3位，即為商數的 8 倍
`*div << 1` 為商數的兩倍
`((q & ~0x7) + (*div << 1))` 所以這裡相當於商數乘上 10

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

### 測驗三
#### 版本一
在以下程式碼中，使用右移操作來計算 log 值

int ilog2(int i)
{
    int log = -1;
    while (i) {
        i >>= 1;
        log++;
    }
    return log;
}

以 8 為例
| round | log | i | i >>= 1 | log++ |
| ----- |:---:|:----:|:-------:|:-----:|
| 1 | -1 | 1000 | 100 | 0 |
| 2 | 0 | 100 | 10 | 1 |
| 3 | 1 | 10 | 1 | 2 |
| 4 | 2 | 1 | 0 | 3 |

以 7 為例
| round | log | i | i >>= 1 | log++ |
| ----- |:---:|:----:|:-------:|:-----:|
| 1 | -1 | 0111 | 011 | 0 |
| 2 | 0 | 011 | 01 | 1 |
| 3 | 1 | 01 | 0 | 2 |

由以上觀察可發現 shift 操作次數是由 msb 在第幾個位元而定，所以這個實作也可以藉由 counting leading zero 來完成

#### 版本二

以下程式碼是將版本一改善，假設輸入資料為 65536 ，在版本一中的 while 迴圈次數需要 65536+1 次，而在版本二中只需要 1 次的迴圈，對於輸入資料很大時，可以減少所需的時間。

- while (i) {
-    i >>= 1;
-    log++;
- }
+ 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;
+    }

#### 版本三
使用 GNU extension `__builtin_clz` 來計算最高位前方有幾個 0 ，在透過 31 - 0 個個數即為最高位元的位置，此值即為 log 的值。

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

### 測驗四
EWMA(指數加權移動平均) 是種取平均的統計手法，並且使經過時間越久的歷史資料的權重也會越低，以下為 EWMA 的數學定義:

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

以下為 ewma 的結構， internal 為 $S_{t-1}$ ，factor 為縮放因子， weight 為權重

/* Exponentially weighted moving average (EWMA) */
struct ewma {
    unsigned long internal;
    unsigned long factor;
    unsigned long weight;
};

考慮以下實作，

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

$\alpha = \dfrac{1}{2^{weight}}$
$1- \alpha = \dfrac{2^{weight} - 1}{2^{weight}}$
若 avg->internal 為 0 ，則為 $2^{factor} \cdot val$
若 avg->internal 不為 0 ，則為
$\dfrac{2^{weight}-1}{2^{weight}} \cdot internal - \dfrac {2^{factor}}{2^{weight}} \cdot val$
### 測驗五
以下程式碼計算 $\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;       
}

以下用於判別 msb 在左邊 16 位元或右 16 位元，若在左 16 位元則 `x > 0xFFFF` 條件成立，則將 x shift 16位元。

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

以下用於判別 msb 在這 16 位元中的左 8 位元或右 8 位元，r 用於紀錄 shift 幾個 bits

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

以此類推即可求出 $\lceil \log_2 x\rceil$ 之值

#### 改進程式碼
> [commit - 8d405ff](https://github.com/jimmylu890303/linux2024_lab/commit/8d405ffb572c53874717f777cf3dce9296f7a8ec)

第一個為 x=0 時，做 x-- 會導致 x 值變成 `0xffffffff`，所以作了以下改進
`x -= (!!x)` 會使 x > 0 時才會作減法的動作。
第二個為 x=1 時，應該要回傳 0 的值，但是在原本的程式中會回傳 1 ，所以去檢查 `!(r | shift | x > 0)` 有無成立，若該條件成立代表 x = 0 或者 x = 1 的狀態，再將回傳值清零。

-    x --;
+    x -= (!!x);
-    return (r | shift | x > 1)+1;
+    x = !(r | shift | x > 0) ? 0 : (r | shift | x > 1)+1;
+    return x;       

#### ceil_log2 in Linux
> [linux/tools/testing/selftests/mm/thuge-gen.c](https://github.com/torvalds/linux/blob/master/tools/testing/selftests/mm/thuge-gen.c#L51)

int ilog2(unsigned long v)
{
	int l = 0;
	while ((1UL << l) < v)
		l++;
	return l;
}

## 第四週測驗題
> [題目](https://hackmd.io/@sysprog/linux2024-quiz4)
### 測驗一

#### popcount

unsigned popcount_naive(unsigned v)
{
    unsigned n = 0;
    while (v)
        v &= (v - 1), n = -(~n);
    return n;
}

在最初的 popcount實作中，是使用迴圈將 LSB 的值清 0 ，
這邊使用特別操作的手法為 `v &= (v-1)` ， 此外 `v &= (v-1)` 的操作也可以用於判斷是否為 2 的冪次方。

0001 0100  # 20          ; LSB in bit position 2
0001 0011  # 20 - 1
0001 0000  # 20 & (20 - 1)

`n = -(~n)` 為 n++ 的操作， $- N = \sim N + 1$ 為二的補數操作，則 $- (\sim N) = N + 1$

#### branchless popcount

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

根據以下公式，我們可以計算 popcount(x)
$popcount(x) = x - \lfloor \dfrac{x}{2} \rfloor - \lfloor \dfrac{x}{4} \rfloor - \lfloor \dfrac{x}{8} \rfloor - ... - \lfloor \dfrac{x}{2^{31}} \rfloor = \sum\limits_{n=0}^{31}(2^n - \sum\limits_{i=0}^{n-1}2^i)b_n = \sum\limits_{n=0}^{31}b_n$

以下操作為計算 $x - \lfloor \dfrac{x}{2} \rfloor - \lfloor \dfrac{x}{4} \rfloor - \lfloor \dfrac{x}{8} \rfloor$
`0x77777777` 的設定是因為這邊是以 4 個 bits 為單位 (nibble) ，為了防止另外一組的 nibble 值 shift 過來。

n = (v >> 1) & 0x77777777;
v -= n;
n = (n >> 1) & 0x77777777;
v -= n;
n = (n >> 1) & 0x77777777;
v -= n;

所以操作完後的值會結果如下(只考慮最右邊的 nibble)，這個 nibble 即代表這 4 個 bits 的 popcount 之值 :
$(2^3-2^2-2^1-2^0)b_3 + (2^2-2^1-2^0)b_2 + (2^1-2^0)b_1 + 2^0b_0 = b_3+b_1+b_1+b_0$

 3       2          1             0
b3 (b2-b3) (b1-b2-b3) (b0-b1-b2-b3)

假設 $B_n$ 為 nibble ，所以以下操作可以看出它的結果

v = (v + (v >> 4)) & 0x0F0F0F0F;

B7 B6 B5 B4 B3 B2 B1 B0  // v
 0 B7 B6 B5 B4 B3 B2 B1  // (v >> 4)
// (v + (v >> 4))
B7 (B7+B6) (B6+B5) (B5+B4) (B4+B3) (B3+B2) (B2+B1) (B1+B0)
// (v + (v >> 4)) & 0x0F0F0F0F
0 (B7+B6) 0 (B5+B4) 0 (B3+B2) 0 (B1+B0)

最後在透過 `v *= 0x01010101` 將各個 nibble 相加

A6 = B7+B6 A4 = B5+B4 A2 = B3+B2 A0 = B1+B0

                         0 A6  0 A4  0 A2  0 A0
                      x  0  1  0  1  0  1  0  1
---------------------------------------------------
                         0 A6  0 A4  0 A2  0 A0
                   0 A6  0 A4  0 A2  0 A0  0
             0 A6  0 A4  0 A2  0 A0  0
       0 A6  0 A4  0 A2  0 A0  0
---------------------------------------------------
                            ↑_______________________A6+A4+A2+A0

可以發現 A6+A4+A2+A0 是我們期望的結果，所以透過 v >> 24 取得。

#### 利用 popcount 計算 Total Hamming distance

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

使用 XOR 運算 可以得出兩數字的差異，在利用巨集 `__builtin_popcount` 即可以算出兩者的 Hamming Distance ，而 Total Hamming distance 就是每個數字都跟其他數字相比，所以時間複雜度會被限制在 $N^2$ 。

#### 改善程式碼

int totalHammingDistance(int* nums, int numsSize)
{
    int total = 0;
    for(int i=0 ;i<32;i++){
        int sum=0;
        for(int n=0;n<numsSize;n++){
            int lsb = (nums[n] & 0x1);
            sum += lsb;
            nums[n]>>=1;
        }
        total += sum * (numsSize-sum);
    }
    return total;
}

假設輸入資料為 4, 14, 4 三筆資料，我們從最低位元到最高位元將每個數字的該位元加總，可以發現以下規律 : 該位元的加總 * (數字個數 - 該位元的加總) = 該位元的 Total Hamming distance ，並且時間複雜度為 $N$

 numsize = 3
---------------------------------------------
 4 = 0100
14 = 1110
 4 = 0100
 
 0th bit : 0 + 0 + 0 = 0 // 0 * (numsize - 0) = 0
 1th bit : 0 + 1 + 0 = 1 // 1 * (numsize - 1) = 2
 2th bit : 1 + 1 + 1 = 3 // 3 * (numsize - 3) = 0
 3th bit : 0 + 1 + 0 = 1 // 1 * (numsize - 1) = 2
 
 total = 0 + 2 + 0 + 2 = 4

### 測驗二