# 2024q1 Homework4 (quiz3+4)
contributed by < [`jimmylu890303`](https://github.com/jimmylu890303)>
## 第三週測驗題
> [題目](https://hackmd.io/@sysprog/linux2024-quiz3)
### 測驗一
#### 版本一
```c
#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 的方式去找出最高位的位元。
```diff
- int msb = (int) log2(N);
+ int msb = 0;
+ int n = N;
+ while (n > 1) {
+ n >>= 1;
+ msb++;
+ }
```
#### 版本三
在版本三中使用 `__builtin_clz(x)` 的巨集來找出最高位在哪個位元算出 leading zero 個數。
```c
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$
```c
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$
```c
((tmp >> 3) + (tmp >> 1) + tmp)
```
由於往右作 shift 3個位元會造成最右邊三個位元資訊損失，所以在往左 shift 後要將最右邊三個位元加回來。所以以下操作值會等於 $13 \times tmp$
```c
((((tmp >> 3) + (tmp >> 1) + tmp) << 3) + d0 + d1 + d2)
```
所以最後往右 shift 7個位元相當於除上128，這樣就會得到 $tmp \div 10$ 的結果
```c
q = ((((tmp >> 3) + (tmp >> 1) + tmp) << 3) + d0 + d1 + d2) >> 7;
```
而 tmp 取 10 的餘數則利用餘數定理 $r = f - g \cdot Q$
```c
r = tmp - (((q << 2) + q) << 1);
```
q 為商數，則以下操作為 $q \times 10$，即可算出餘數 r
```c
(((q << 2) + q) << 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));
}
```
在以下操作中， $x = \dfrac{3}{4} \times in$
```c
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$
```c
uint32_t q = (x >> 4) + x;
```
而這邊的操作是將 0.79 逼近至 0.8
```c
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$
```c
*div = (q >> 3);
```
考慮以下操作
`q >> 3` 為商數
而 `q & ~0x7` 將最低三位清零，相當於商數往左 shift 3位，即為商數的 8 倍
`*div << 1` 為商數的兩倍
`((q & ~0x7) + (*div << 1))` 所以這裡相當於商數乘上 10
```c
*mod = in - ((q & ~0x7) + (*div << 1));
```
### 測驗三
#### 版本一
在以下程式碼中，使用右移操作來計算 log 值
```c
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 次的迴圈，對於輸入資料很大時，可以減少所需的時間。
```diff
- 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 的值。
```c
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 為權重
```c
/* Exponentially weighted moving average (EWMA) */
struct ewma {
unsigned long internal;
unsigned long factor;
unsigned long weight;
};
```
考慮以下實作，
```c
/**
* 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$
```c
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位元。
```c
r = (x > 0xFFFF) << 4;
x >>= r;
```
以下用於判別 msb 在這 16 位元中的左 8 位元或右 8 位元，r 用於紀錄 shift 幾個 bits
```c
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 的狀態，再將回傳值清零。
```diff
- 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)
```c
int ilog2(unsigned long v)
{
int l = 0;
while ((1UL << l) < v)
l++;
return l;
}
```
## 第四週測驗題
> [題目](https://hackmd.io/@sysprog/linux2024-quiz4)
### 測驗一
#### popcount
```c
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
```c
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 過來。
```c
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 ，所以以下操作可以看出它的結果
```c
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
```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;
}
```
使用 XOR 運算 可以得出兩數字的差異，在利用巨集 `__builtin_popcount` 即可以算出兩者的 Hamming Distance ，而 Total Hamming distance 就是每個數字都跟其他數字相比，所以時間複雜度會被限制在 $N^2$ 。
#### 改善程式碼
```c
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
```
### 測驗二
Sign in with Wallet
Connect another wallet