# 2024q1 Homework4 (quiz3+4)
contributed by < [`kkkkk1109`](https://github.com/kkkkk1109) >
## [2024q1 第3周測驗題](https://hackmd.io/@sysprog/linux2024-quiz3)
### 測驗 `1`
測驗一主要為使用 `bitwise` 的方式來進行平方根的計算，以下為版本一
```clike
#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;
}
```
`msb` 的部分是用來計算 `N` 是由幾個 bit 所組成的，以 `N` = 10 為例
```clike
log2(10) = 3.321928
```
使用 `(int)` 強制轉成整數
```clike
(int)log2(10) = 3
```
接著將 1 左移 3 位，得到 8。
```clike
int a = 1 << msb
```
以 `result + a` 的平方檢查目前是否小於 N，若小於，則代表此數可以構成平方根，並不斷將 a 右移以檢查下個 bit，直到 a 為零。
```clike
int result = 0;
while (a != 0) {
if ((result + a) * (result + a) <= N)
result += a;
a >>= 1;
}
return result;
```
版本二為不需要使用 `log` 的計算，利用不斷右移的方式，計算 bit 個數。
```clike
int msb = 0;
int n = N;
while (n > 1) {
n >>= 1;
msb++;
}
```
版本三是使用 [Digit-by-digit calculation](https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Digit-by-digit_calculation)的方式進行，其說明如下
若要求得 `x` 的平方根，可以將 `x` 表示成以下二進位的形式
$$ x = (000b_0b_1b_2b_3.....b_n)^2 $$
也可以表示成2的冪相加
$$x =(a_n+a_{n-1}+a_{n-2}+....+a_0)^2$$
以 `x = 36` 舉例，
$$ x = (000....110)^2 $$
也可以表示為
$$ x = (2^2+2^1+0\times2^0)^2$$
此方法是將 $(a_n+a_{n-1}+a_{n-2}+....+a_0)^2$ 展開，並從係數來回推$(a_n+a_{n-1}+a_{n-2}+....+a_0)$。
另 $$N = (a_n+a_{n-1}+a_{n-2}+....+a_0)$$
平方為
\begin{split} \\
N^2 =& a_n^2+[2a_n+a_{n-1}]a_{n-1}+[2(a_n+a_{n-1})+a_{n-2}]a_{n-2}+...+[2(\displaystyle\sum_{i=1}^{n}a_i)+a_0]a_0 \\
=& a_n^2+[2P_n+a_{n-1}]a_{n-1}+[2P_{n-1}+a_{n-2}]a_{n-2}+...+[2P_1+a_0]a_0 \\
\end{split}
其中 $P_m= P_{m+1} + a_m$ 為
$$P_m =a_n+a_{n-1}+...+a_m$$
因此，當 $P_0 = (a_n+a_{n-1}+a_{n-2}+....+a_0)$ 時即為所求 $N$
且$$P_m= P_{m+1} + a_m$$
$$a_m = 2^m\;or\;0$$
只要從 $a_n$ 試到 $a_0$，且滿足 $P_m^2<=N^2$，即可求出平方根。
因此可以從 $P_n$ 開始，持續加上 $a_{n-1},a_{n-2}...a_0$，並於每次加上後檢查是否滿足 $P_m^2<=N^2$ 即可。
為了減少運算平方的成本，可以改成以下方式進行
$$X_m = N^2 -P_m^2 = X_{m+1}-Y_m$$
$$Y_m =P_m^2 - P_{m+1}^2$$
以 `36` 舉例
$$ N = (2^2+2^1+0)^2$$
$a_n=4 ,a_n-1 = 2,a_0 = 0$
第一輪:
$$X_2 = N^2 -P_2^2 = 36-16=20$$
$$Y_2 = 20$$
第二輪
$$X_1 = N^2 -P_1^2 = X_2 - Y_1 = 0$$
$$Y_1 = 36-16=20$$
此輪結束時，差值$X_1 = 0$，因此無須進行下一輪。不過在計算 $Y_m$ 時，仍需要計算平方，因此可再將 $Y_m$ 分解
$$Y_m = P_m^2 - P_{m+1}^2 = 2P_{m+1}a_m+a_m^2$$
令 $c_m,d_m$
$$c_m = 2P_{m+1}a_m = P_{m+1}2^{m+1} \;(\because a_m = 2^m)$$
$$d_m = a_m^2 = 2^{2m}$$
可以求得 $c_{m-1},d_{m-1}$
$$c_{m-1} = P_m2^m =(p_{m+1}+2^m)2^m=c_m/2 + d_m $$
$$d_{m-1} = 2^{2m-2}=d_m/4$$
以上情況為$a_m = 2^m$時，若 $a_m=0$，則
$$Y_m = 0, \; c_{m-1} =c_m/2$$
這樣就可以使用移位的方式來迭代，不需要使用到平方
首先，檢查是否有平方根，由於此題輸入輸出為整數，因此只需檢查是否小於等於 1
```clike
if (x <= 1) /* Assume x is always positive */
return x;
```
由於 $d_m$ 在迭代到最後會為 0，因此我們可以以 $d_m$ 來判斷是否迭代結束
以 `m` 為 $d_m$，可利用 `31 - CLZ` 來獲得位數，再進行移位得到 $d_m$。而每次迭代，$d_{m-1} = d_m/4$，因此要右移兩位
```clike
for (int m = 1UL << ((31 - __builtin_clz(x)) & ~1UL); m; m >>= 2) {
...
}
```
以 `z` 為 $c_m$， `b` 為 $Y_m$ ， `x` 為 $X_m$
`b = z + m` 求得 $Y_m$，`z >> 1` 對 $c_{m-1} = c_m/2$ 進行迭代，並以 `(x >= b)`檢查 $a_m$是否為 0。
```clike
int z = 0;
for (int m = 1UL << ((31 - __builtin_clz(x)) & ~1UL); m; m >>= AAAA) {
int b = z + m;
z >>= 1;
if (x >= b)
x -= b, z += m;
}
return z;
```
迭代最後，會發現是平方根為 `z`，解釋如下
$$c_0 = 2P_{1}a_0$$
$$P_0= P_{1} + a_0 , P_{1} = P_0 - a_0$$
$$c_0 = 2(P_0 - a_0)a_0 = 2P_0a_0 -2a_0^2$$
$$P_0 = (c_0+2a_0^2)/2 =c_0/2 + a_0^2$$
這時會發現 $a_0^2$ 就是 $d_m$ ，對應到最後一次迭代會發現
`z`會右移，代表 $c_0/2$，又 `z += m` 代表加上了 $d_m$，即為所求 $P_0$。
使用[第 2 週測驗題-測驗三](https://hackmd.io/@sysprog/linux2024-quiz2#%E6%B8%AC%E9%A9%97-3) 中的`ffs`可以取代 `GNU_extension` 中的 `__builtin_clz`。
我使用了測驗三中的 `ffs` 和 `clear_bit` 來取代利用CLZ來找出最高位。
> [`commit`](https://github.com/kkkkk1109/2024q1-Homework4/commit/6bb08d3814ac399df3c55046f8125e0897e5eaa4)
### 測驗 `2`
此題為進行 mod 10 和 div 10操作時，要如何降低運算成本
由於 $10 = 2\times 5$，使用 bitwise 進行除法時，不能夠完全整除。因此要找到一個除數 `x`，求出來的結果會和以 10 為除數相同。
閱讀《[Hacker's Delight](https://web.archive.org/web/20180517023231/http://www.hackersdelight.org/divcMore.pdf)》，提到可以以找出乘上倒數的方式來進行除法，如
$$x/10=x\times{1/10}$$
以 `3` 做舉例
$$1/3=0.0101 0101 0101 0101 0101 0101 0101 0101.$$
我們可以將被除數透過位移的方式來完成除法
```clike
unsigned divu3(unsigned n) {
unsigned q, r;
q = (n >> 2) + (n >> 4); // q = n*0.0101 (approx).
q = q + (q >> 4); // q = n*0.01010101.
q = q + (q >> 8);
q = q + (q >> 16);
r = n - q*3; // 0 <= r <= 15.
return q + (11*r >> 5);
```
但要注意的是，每次右移都有可能造成誤差，會直接忽略低位的位元，由於上述程式碼右移了 5 次，因此餘數 `r` 會介在 $3*5+(3-1)=17$ 之間，其中 3 為除數。參考同學 [`vax-r`](https://github.com/vax-r) 所做的[實驗](https://hackmd.io/@vax-r/linux2024-homework2#bitwise-operation-%E8%A3%9C%E5%BC%B7)，記錄了最大餘數的變化
以 `hexa` 的形式來看
```shell
number contribute to max reminder 1 ,max is 1
number contribute to max reminder 2 ,max is 2
number contribute to max reminder 3 ,max is 3
number contribute to max reminder 7 ,max is 4
number contribute to max reminder b ,max is 5
number contribute to max reminder f ,max is 6
number contribute to max reminder 1f ,max is 7
number contribute to max reminder 2f ,max is 8
number contribute to max reminder cf ,max is 9
number contribute to max reminder 1cf ,max is 10
number contribute to max reminder 2cf ,max is 11
number contribute to max reminder c9cf ,max is 12
number contribute to max reminder 1c9cf ,max is 13
number contribute to max reminder 2c9cf ,max is 14
number contribute to max reminder c9c6c9cf ,max is 15
```
可以發現確實會隨著數字的增加，最大餘數會持續增加。
以數字 `c9c6c9cf` 來查看每次右移的結果
```clike
q = (i >> 2) + (i >> 4);
q = 3f0e1f0f
q = q + (q >> 4);
q = 42ff00ff
q = q + (q >> 8);
q = 4341ffff
q = q + (q >> 16);
q = 43424340
```
可以發現剛好右移的位數都是 `1`，右移 4 ，0~3 bit 為 1；右移 8 ，0~7 bit 為 1，因此造成了誤差 15，也如同文章所做的實驗
> Thus, using (33), for which q is too low by at most 5, we expect the remainder to be at most 5*3 + 2 = 17. Experiment reveals that it is actually at most 15
回到測驗題中，說明以下程式碼中的餘數是介在 `0` ~ `19`
```
carry = tmp / 10;
tmp = tmp - carry * 10;
```
> 然而，觀察上面的程式碼後可以發現， tmp 不可能會大於 19 ，因此只需要考慮 19~0 的情況即可。
:::info
剛開始不太懂為何是 0 ~ 19，猜測是因為 `carry` 所造成的誤差最大就是 1，而餘數最大值為 9，因此最大可能餘數為 $1\times10+9=19$
:::
再以 `Hacker's Delight`中的 `divu10` 測試，最大餘數也會到達 `13`。
```
number contribute to max reminder 1 ,max is 1
number contribute to max reminder 2 ,max is 2
number contribute to max reminder 3 ,max is 3
number contribute to max reminder 4 ,max is 4
number contribute to max reminder 5 ,max is 5
number contribute to max reminder 6 ,max is 6
number contribute to max reminder 7 ,max is 7
number contribute to max reminder 8 ,max is 8
number contribute to max reminder 9 ,max is 9
number contribute to max reminder a ,max is 10
number contribute to max reminder b ,max is 11
number contribute to max reminder 296a ,max is 12
number contribute to max reminder 296b ,max is 13
```
若我們想將控制精度在小數點後第一位，則可以表示為
$$1.9\le{\frac{19}x}\le1.99, \;9.55\le{\frac{19}x}\le10$$
接著找 `x` 的範圍，可以找到 ${\frac{128}{13}}\approx9.84$
接下來就可以透過以下方式實作
```Clike
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);
```
`d0`、`d1`、`d2` 是用來儲存因為右移而消失的位數，透過右移的方式來進行除以 `13`，再左移來完成乘上 `128`，最後再將 `tmp` - `10q`。
包裝成以下方式
```clike
#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 >> (CCCC=2));
*mod = in - ((q & ~0x7) + (*div << (DDDD=1));
}
```
此方法為改成 $x\times{\frac{0.8}8}$的方式來運算。為了得到0.8，會先近似於0.75，再慢慢逼近到 0.8 最後再除上 8。
將`100000000`作為輸入，並查看每次運算的結果
```clike
in = 100000000
x= 75000001
q= 79687501
q= 79998780
q= 79999996
q= 80000000
q= 80000001
```
由於 `0.8` 在 32位元表示是 `0.11001100110011001100110011001100`，可以先產生前面的`1100`11
以下可以先產生 `0.110011`
```clike
uint32_t x = (in | 1) - (in >> 2); /* div = in/10 ==> div = 0.75*in/8 */
uint32_t q = (x >> 4) + x;
```
再進行 4 次的位移和加法
```clike
x = q; //110011_(6bit)
q = (q >> 8) + x;//11001100110011_(14bit)
q = (q >> 8) + x;//1100110011001100110011_(22bit)
q = (q >> 8) + x;//110011001100110011001100110011_(30bit)
q = (q >> 8) + x;//11001100110011001100110011001100110011_(38bit)
```
最後，因為 `q` 目前是 $0.8x$，右移二來除上 8 得到商 $0.1x$，而餘數則是利用 $x- (8\times0.1x+2\times0.1x)$ 來取得， `q & ~0x7`是為了確保加上 `div << 1` 不會進位。
```clike
*div = (q >> (CCCC=2));
*mod = in - ((q & ~0x7) + (*div << (DDDD=1));
```
### 測驗 `3`
此題為計算以 2 為底的對數，且輸入和輸出均為整數。
**版本一**
逐位檢查是否為 1 ，來計算對數
```clike
int ilog2(int i)
{
int log = -1;
while (i) {
i >>= 1;
log++;
}
return log;
}
```
**版本二**
括號中的數值為$2^{result}$，此方法為一次大範圍的計算對數，比版本一逐位計算，此方法減少了計算成本。
```clike
static size_t ilog2(size_t i)
{
size_t result = 0;
while (i >= (AAAA=65536) {
result += 16;
i >>= 16;
}
while (i >= (BBBB=256){
result += 8;
i >>= 8;
}
while (i >= (CCCC=16)
result += 4;
i >>= 4;
}
while (i >= 2) {
result += 1;
i >>= 1;
}
return result;
}
```
**版本三**
利用 `CLZ` 來求得此數的前導零數量，以 `0x00020100` 舉例， `CLZ` 為 12，可快速透過 $31-12$ 來得到對數。
```clike
int ilog32(uint32_t v)
{
return (31 - __builtin_clz(DDDD=v|1));
}
```
同學 [MathewSu-001](https://github.com/MathewSu-001) 提醒了我，由於 [gcc](https://gcc.gnu.org/onlinedocs/gcc/Other-Builtins.html) 中有說明道，若輸入為 0 ，會造成 undefined result，因此要加上 `v | 1` 來避免此行為發生。
> Returns the number of leading 0-bits in x, starting at the most significant bit position. If x is 0, the result is undefined.
### 測驗 `4`
[Exponential 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.$
* $\alpha$表示歷史資料加權降低的程度，介在 0 ～ 1 之間，越高的$\alpha$會使歷史資料減少的越快
* $Y_t$表示在時間$t$時的資料點
* $S_t$表示在時間$t$時計算出的 EWMA
將其展開來如下
$\begin{split}
S_t &= \alpha Y_t + (1 - \alpha)\cdot S_{t-1} \\
&= \alpha Y_t + (1 - \alpha)(\alpha Y_{t-1} + (1 - \alpha)\cdot S_{t-2}) \\
&= \alpha Y_t + (1 - \alpha)(\alpha Y_{t-1} + (1 - \alpha)(\alpha Y_{t-2} + (1 - \alpha)\cdot S_{t-3})) \\
&= \alpha (Y_t + (1 - \alpha)Y_{t-1} + (1 - \alpha)^2Y_{t-2} + \ldots + (1 - \alpha)^kY_{t-k} + \ldots + Y_0)
\end{split}$
隨著時間 $t$ 越來越長，越久以前的資料權重也會越低。
首先，EWMA的結構由三個參數組成，分別是 `internal`, `factor` 和 `weight`。 `internal` 是最後取的平均值，對應到公式中的 $S_t$， `factor` 為 scaling factor，`weight` 即為加權重量。
```clike
/* Exponentially weighted moving average (EWMA) */
struct ewma {
unsigned long internal;
unsigned long factor;
unsigned long weight;
};
```
我們可以藉由 `ewma_init` 來進行初始化，並檢查其中輸入的 `factor` 和 `weight` 是否為 2 的對數，接下來將他們的對數值存入`factor` 和 `weight`中。
```clike
void ewma_init(struct ewma *avg, unsigned long factor, unsigned long weight)
{
if (!is_power_of_2(weight) || !is_power_of_2(factor))
assert(0 && "weight and factor have to be a power of two!");
avg->weight = ilog2(weight);
avg->factor = ilog2(factor);
avg->internal = 0;
}
```
接著，可以使用 `ewma_add` 來增加樣本。
加入新樣本的公式為
$$S_t = \alpha Y_t+(1-\alpha)\times S_{t-1}$$
$S_{t-1}$ 對應到程式中的 `internal`，$\alpha$ 對應到 `weight`， $Y_t$ 對應到 `val`。整理公式可以變成$$S_t = \alpha (Y_t+ ( \frac{1}{\alpha} -1)\times S_{t-1})$$
即對應到程式中，當 `avg->internal` 存在時的敘述。
```clike
struct ewma *ewma_add(struct ewma *avg, unsigned long val)
{
avg->internal = avg->internal
? (((avg->internal << EEEE = avg->weight) - avg->internal) +
(val << FFFF = avg->factor)) >> avg->weight
: (val << avg->factor);
return avg;
}
```
公式中的 $\alpha$ 介在 0 ~ 1 之間的數，但在此 `weight` 為大於零的 2 的冪，因此公式中的 $\alpha$原本進行乘法和除法，會在程式中相反。如$\frac{1}{\alpha}\times S_{t-1}$ 在程式中卻是進行左移，而括號外的乘上 $\alpha$ 卻變成了右移。
### 測驗 `5`
此測驗是延伸[測驗 `3`](https://hackmd.io/t3xU5PpDRiinwgc37y_kqQ?view#%E6%B8%AC%E9%A9%97-3)，已知輸入必為大於 0 的數值 (即 `x` > 0)，以下程式碼可計算 $\lceil log_2(x) \rceil$，也就是 [ceil](https://man7.org/linux/man-pages/man3/ceil.3.html) 和 [log2](https://man7.org/linux/man-pages/man3/log2.3.html) 的組合並轉型為整數:
以下為程式碼
```clike
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 > (GGG = 1)) + 1;
}
```
首先，由於對數會比位數少一位，如 `32 = 0b100000` 位數是6，但對數是 5，因此要先減 1，以完成後續計算
```clike
x--;
```
在不斷的比較中， `r` 來儲存所有比較的結果， `shift` 儲存當前比對結果，利用 `r |= shift` 來紀錄每次的結果
```clike
r = (x > 0xFFFF) << 4; //r = 10000(16)
x >>= r; // x /= 16
shift = (x > 0xFF) << 3; // shift = 1000(8)
x >>= shift; // x /= 8
r |= shift; // r=1100
shift = (x > 0xF) << 2; // shift = 100(4)
x >>= shift; // x /= 4
r |= shift; // r=1110
shift = (x > 0x3) << 1; // shift = 10(2)
x >>= shift; // x /=2
return (r | shift | x > (GGG = 1)) + 1;
```
在比較到最後，可以發現 `x` 最後一位沒有進行比較，且最後一次的 `shift` 沒有被儲存到 `r` 中， 直接在 `return` 中完成上述操作。
:::info
這邊注意到，最後再return 又再加上了 1 ，而若沒有第一步的 `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);
}
```
兩者測試發現，只有在 `x` 為 2的冪時才會相同，才想到此題還有 `ceil`，才會造成兩者結果不一樣。
## [2024q1 第4周測驗題](https://hackmd.io/@sysprog/linux2024-quiz4)
### 測驗 `1`
[popcount](https://en.wikichip.org/wiki/population_count)可以用來計算二進位的表示式中，有多少位元為1，如 `0b10101010` 的 popcount 為 4。
`popcount_naive` 每次將 `v` 的 LSB 位設為零，直到 `v`為零，使用迴圈計算 1 的個數
```cpp
unsigned popcount_naive(unsigned v)
{
unsigned n = 0;
while (v)
v &= (v - 1), n = -(~n);
return n;
}
```
對於 32 bit 的無號整數，popcount的公式推導如下
$popcount(x)= x - [\frac{x}{2}] - [\frac{x}{4}]-...-[\frac{x}{31}]$
以 `x` = `0b101011`舉例
```
101011
- 10101
_______
10110
- 1010
_______
1100
- 101
_______
0111
- 10
_______
0101
- 1
_______
100(4)
```
計算到最後可以得到 4。
假設 $x = b_{31}...b_3b_2b_1b_0$，以 4 個位元 $x[3:0]$來看，其 popcount為
$$(2^3b_3 + 2^2b_2 + 2^1b_1 + 2^0b_0) -
(2^2b_3 + 2^1b_2 + 2^0b_1) - (2^1b_3 + 2^0b_2) - 2^0b_3$$
整理一下
$$(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$$
最後公式可以整理成
$$popcount(x) = \sum\limits_{n=0}^{31} {}(2^n - \sum\limits_{i=0}^{n-1} 2^{i})b_n = \sum\limits_{n=0}^{31}b_n$$
`popcount_naive` 使用迴圈的方式來進行，而我們可以實作常數時間複雜度的做法 `popcount_branchless`
```cpp
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;
}
```
此方法是將 4 個位元(nibble) 為一個單位計算，因此只要使用$x - \left \lfloor{{\dfrac{x}{2}}}\right \rfloor - \left \lfloor{{\dfrac{x}{4}}}\right \rfloor - \left \lfloor{{\dfrac{x}{8}}}\right \rfloor$計算即可
```cpp
n = (v >> 1) & 0x77777777;
v -= n;
```
首先，相當於做 $v - \left \lfloor{{\dfrac{v}{2}}}\right \rfloor$，右移是為了表示 $\frac v2$，而 `0x77777777`作為遮罩，完成一個 nibble 為單位計算
```
b_31 b_30 b_29 b_28 ... b7 b6 b5 b4 b3 b2 b1 b0 // v
0 b_31 b_30 b_29 ... 0 b7 b6 b4 0 b3 b2 b1 // (v >> 1) & 0x77777777
```
重複上述動作，完成$v - \left \lfloor{{\dfrac{v}{2}}}\right \rfloor - \left \lfloor{{\dfrac{v}{4}}}\right \rfloor - \left \lfloor{{\dfrac{v}{8}}}\right \rfloor$
```cpp
b_31 b_30 b_29 b_28 ... b7 b6 b5 b4 b3 b2 b1 b0 // v
- 0 b_31 b_30 b_29 ... 0 b7 b6 b4 0 b3 b2 b1 // (v >> 1) & 0x77777777
- 0 0 b_31 b_30 ... 0 0 b7 b6 0 0 b3 b2
- 0 0 0 b_31 ... 0 0 0 b7 0 0 0 b3
```
接著，將每組所算出的 popcount 加起來
```cpp
v = (v + (v >> 4)) & 0x0F0F0F0F;
v *= 0x01010101;
```
以 $B_n$ 代表第 n 個 nibble 中的數值
```cpp
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)
```
使用 mask `0x0F0F0F0F`
```cpp
// (v + (v >> 4)) & 0x0F0F0F0F
0 (B7+B6) 0 (B5+B4) 0 (B3+B2) 0 (B1+B0)
```
最後乘上 `0x01010101`，右移 24 就會是最終的答案
```cpp
v * 0x01010101 =
(A + B + C + D) (B + C + D) (C + D) (D)
|<-- 1 byte -->|<-- 1 byte -->|<-- 1 byte -->|<-- 1 byte -->|
```
對於要計算兩數的 [Hamming distance](https://en.wikipedia.org/wiki/Hamming_distance)，將兩數做 `xor` 再取 popcount 即可獲得。
```cpp
int hammingDistance(int x, int y)
{
return __builtin_popcount(x ^ y);
}
```
針對 [LeetCode 477. Total Hamming Distance](https://leetcode.com/problems/total-hamming-distance/)，考慮以下程式碼:
```cpp
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 >> (AAAA = 1);
}
```
因為此方法的 `i` 和 `j` 都是由 `0`~ `numsize-1` ，因此會計算到重覆的 hamming distance，因此最後 `totla` 要右移一位。
此方法的時間複雜度為$O(n^{2})$，我們可以改寫 `j = i` 做內部迴圈，但時間複雜度仍然是$O(n^{2})$
```diff
int totalHammingDistance(int* nums, int numsSize)
{
int total = 0;;
for (int i = 0;i < numsSize;i++)
-- for (int j = 0; j < numsSize;j++)
++ for (int j = i; j < numsSize;j++)
total += __builtin_popcount(nums[i] ^ nums[j]);
return total >> (AAAA = 1);
}
```
換個角度思考，從位元的角度來看
| n 'th bit | 4 | 3 | 2 | 1 | 0 |
| --------- | --- |:---:|:---:|:---:| --- |
| input 7 | 0 | 1 | 1 | 1 | 1 |
| input 5 | 0 | 0 | 1 | 0 | 1 |
| input 10 | 0 | 1 | 0 | 0 | 1 |
| input 17 | 1 | 0 | 0 | 0 | 1 |
* 第 `0` 位: 全 1 ，hamming distance = 0
* 第 `1` 位: 一個 1 ，hamming distance = 1 * `numsize-1`
* 第 `2` 位: 兩個 1 ，hamming distance = 2 * `numsize-2`
我們可以直接在位元內，計算 1 的個數 `n`，hamming distance 為 `n` * `numsize - n`，再逐個位元加起來。時間複雜度縮短到 $O(n)$
```cpp
int totalHammingDistance(int* nums, int numsSize)
{
int total = 0,hamming_dis = 0;
for (int i = 0;i < 32;i++){
hamming_dis = 0;
for (int j = 0; j < numsSize;j++)
hamming_dis += (nums[j] >> (31-i) ) & 1 ;
total += hamming_dis * (numsSize - hamming_dis);
}
return total ;
}
```
### 測驗 `2`
若要不使用除法，就獲得某數除以另一數的餘數，可以使用以下方法，只要除數滿足 $2^K\pm1$。以除數為 3 做舉例， $1 = 1(mod\;3)$，$2 = -1(mod \;3) 可寫成以下形式$$$2^k \equiv \begin{cases}
1 (mod \ \ 3), \ \ k \ \ even\\
-1 (mod \ \ 3), \ \ k \ \ odd\\
\end{cases}$$
而考慮某數 `n` ，以二進位表示式 $b_{n-1}b_{n-2}\cdots b_1 b_0$，$n = \sum_{i=0}^{n-1} b_i\cdot 2^i$，我們可以拆成位元來獲得其餘數，舉例，以 `11` = `0b1011`，餘數為 $-1(2^3)-1(2^1)+1(2^0) = -1$。
因此我們可以利用 `n = popcount(n & 0x55555555) - popcount(n & 0xAAAAAAAA)`，計算偶數位和奇數位的 `popcount`，來求出餘數。可以透過以下方式化簡
$$popcount(x \And \overline{m}) - popcount(x \And m) = popcount(x \bigoplus m) - popcount(m)$$
因此 `n = popcount(n & 0x55555555) - popcount(n & 0xAAAAAAAA)` 可以寫成 `n = popcount(n ^ 0xAAAAAAAA) - 16`，不過為了避免出現負數的情況，因此加上了 3 的倍數 39，再進行計算。
```cpp
int mod3(unsigned n)
{
n = popcount(n ^ 0xAAAAAAAA) + 23;//23 <= n <= 55
n = popcount(n ^ 0x2A) - 3; // -3 <= n <= 2
return n + (((int)n >> 31) & 3);
}
```
第一次的結果，`n` 會落在 23 ~ 55 之間，第二次會落在 -3 ~ 2 之間，最後查看是否為負數，若為負數，則再加上 `3`。
:::info
這邊要注意，在 return 時， 右移運算的 `n` 需加上 (int)，unsigned n 右移的話會是補 0，signed 右移才會補上 1，做 `&` 運算才能得到 `3`。
:::
也可以使用 `lookup table` 來完成
```cpp
int mod3(unsigned n)
{
static char table[33] = {2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1 };
n = popcount(n ^ 0xAAAAAAAA);
return table[n];
}
```
原先`n = popcount(n ^ 0xAAAAAAAA) - 16` 的範圍會介在 -16 ~ 16 之間，而這邊等同於直接加上了 16，會讓範圍介在 0~32之間，因此 lookup table 需要 33 個空間，而可以對應到 $mod(n)=-16+(18)=2$，因此 `table[0]`會從 2 開始，$mod(n)=-15+(15)=0$， `table[1]` 為 0，`table[2]` 為 1 ..... `table[32]` 為 1。
在[Hacker's Delight](https://web.archive.org/web/20180517023231/http://www.hackersdelight.org/divcMore.pdf) 中提到，若不使用 `popcount` ，也可以使用除數轉換的方式來求得餘數
$$n\;mod3 = [\frac{4}{3}n]\;mod\;4$$
證明如下:
令 $n=3k+r$，$k$ 和 $r$ 都是整數，且 $0\le r \le 2$，則$$[\frac{4}{3}(3k+r)]\;mod\;4 = [4k+\frac{4}{3}r]mod\;4=[\frac{4}{3}r]mod\;4$$當 $r$ 為 0、1、2 時都成立。
並且，由於 $4 = 1(mod 3)$，可以直接從二進位的方式來判斷餘數，也就是 `00` ~ `11`，剛好對到十進位 `0` ~ `3`，我們可以以 2 個 bit 為一組來計算；另一個概念是，將原本的數字乘上 $2^k$ ($k$ 為偶數)，所求得的餘數不會變，因為 $2^k=1(mod 3)$，因此我們可以將原先的數字分成許多組來相加計算，並且不會影響除數結果，舉例如下
$N =(10111110)_2=(190)_{10}$，$190 \;mod3=1$，把 $N$ 分成$(10+11+11+10)_2=1010_2=10_{10}$，$10\;mod 3=1$ 來求得餘數，也就是將 `N >> 6` + `(N >> 4)& 0b11` + `(N >> 2) & 0b11` + `N & 0b11`，因為 $10000000_2 = 128_{10}$，$10_2 = 2_{10}$，$128\;mod3 = 10\;mod3=2$，因此可以透過移位相加來計算餘數。
```cpp
int remu3(unsigned n) {
n = (n >> 16) + (n & 0xFFFF); // Max 0x1FFFE.
n = (n >> 8) + (n & 0x00FF); // Max 0x2FD.
n = (n >> 4) + (n & 0x000F); // Max 0x3D.
n = (n >> 2) + (n & 0x0003); // Max 0x11.
n = (n >> 2) + (n & 0x0003); // Max 0x6.
return (0x0924 >> (n << 1)) & 3;
}
```
同理，`3` 個 bit 的話，可以表示 0~7 的餘數，
且 $8^k=1(mod7)$，乘上 8 的倍數就不會改變餘數，看 [tictactoe.c](https://gist.github.com/jserv/9088f6f72a35a1fcaaf1c932399a5624) 中的 `mod7` ，就是使用 $n\;mod7 = [\frac{8}{7}n]\;mod\;8$ 來求得。
```cpp
static inline int mod7(uint32_t x)
{
x = (x >> 15) + (x & UINT32_C(0x7FFF));
/* Take reminder as (mod 8) by mul/shift. Since the multiplier
* was calculated using ceil() instead of floor(), it skips the
* value '7' properly.
* M <- ceil(ldexp(8/7, 29))
*/
return (int) ((x * UINT32_C(0x24924925)) >> 29);
}
```
接下來解釋 [`tictactoe.c`](https://gist.github.com/jserv/9088f6f72a35a1fcaaf1c932399a5624)
首先，`move_mask` 有九步，對應到棋盤上的九宮格。
```clike
move mask
ROW1 ROW2 ROW3 COL1 COL2 COL3 DIA1 DIA2
0100 0000 0000 0100 0000 0000 0100 0000
0010 0000 0000 0000 0100 0000 0000 0000
0001 0000 0000 0000 0000 0100 0000 0100
0100 0000 0010 0000 0000 0000 0000
0010 0000 0000 0010 0000 0010 0010
0001 0000 0000 0000 0010 0000 0000
0100 0001 0000 0000 0000 0001
0010 0000 0001 0000 0000 0000
0001 0000 0000 0001 0001 0000
```
可以看出，分別對應到 8 種獲勝方式，`3 個 row`、 `3 個 column` 和 `2 個 diagonal`。因此，對於是否獲勝，只要看 8 個 `nibble` 是否有達到 `0111`，換句話說，加上 `0x11111111` 檢查各個 nibble 的最高位是否為 1。
```cpp
static inline uint32_t is_win(uint32_t player_board)
{
return (player_board + (BBBB = 0x11111111)) & 0x88888888;
}
```
### 測驗`3`
AVL tree 是一種 [Binary Search Tree](https://en.wikipedia.org/wiki/Binary_search_tree)，藉由計算每個節點的高度，並將節點的子節點的高度相減，以判斷是否平衡；若不平衡，則根據目前的情況，進行對應的 rotation 來達到平衡。
`LL`型
```
A B
/ right rotation / \
B -> C A
/
C
```
`RR`型
```
A B
\ left rotation / \
B -> A C
\
C
```
`LR`型
```
A A C
/ left rotation on B / right rotation / \
B -> C -> B A
\ /
C B
```
`RL`型
```
A A C
\ right rotation on B \ left rotation / \
B -> C -> A B
/ \
C B
```
rbtree 也是一種二元搜索樹，而其遵循以下四個性質
1. 根節點始終是黑色的。
2. 所有末端節點 (leaf) 都是黑色的
3. 如果一個節點是紅色的，那麼它的兩個子節點都是黑色的。
4. 任意一條從根節點到葉子節點的路徑上，黑色節點的數量必須相同，以確保樹的平衡性。
在新增的過程中，若要進行新增或刪除，也是按照前述的 right rotation 和 left rotation 之外，還要調整節點的顏色。
AVL tree 在一般搜尋情況下會比還要快，但在新增和刪減的情況下，rb tree 所需要進行的 rotate 比較少。
對於一個兼具AVL tree 和 rb tree 的 XTree，著重在快速的 insert/delete 的操作和合理的 look up 速度。XTree 的平衡機制與 AVL tree 相似，利用 `hint` 來評估是否需要做 rotate，然而 xt_node 的 `hint` 的計算方式是，其左右節點 `hint` 最大值 +1 並且只有在被呼叫到 `xt_update` 時才會被更新。
接著對 [XTree 的程式碼](https://gist.github.com/jserv/7374b3e031ec3a81441b1b93b008c5d2) 做解釋:
XTree 繼承了 AVL Tree 和 rb tree 的優點，並有著以下四種操作 `xt_rotate_left` 、 `xt_rotate_right` 、 `xt_replace_right` 和 `xt_replace_left`。
`xt_rotate_left` 和 `xt_rotate_right` 主要用於平衡二元樹的部分，因此在**插入**或**刪除**時都會使用；而 `xt_replace_right` 和 `xt_replace_left` 主要用在**刪除**的部分，會將需要被刪除的節點的左右子樹中，找到替代的節點。
而 `xt_update` 則判斷目前 XTree 中是否平衡，旋轉後會持續往親代進行 update。
於 [`xtree.h`](https://gist.github.com/jserv/7374b3e031ec3a81441b1b93b008c5d2#file-xtree-h) 中:
節點 `xt_node` 包含了計算其高度的 `hint`、親代 `parent` 和左右節點 `left` 、`right`。
```cpp
struct xt_node {
short hint;
struct xt_node *parent;
struct xt_node *left, *right;
};
```
整顆二元樹則是包含了根節點 `root` 、 用於比較鍵值的 `cmp` 、 包含鍵值的初始節點 `create_node` 和函式 `destroy_node` 指向被刪除的節點。
```cpp
struct xt_tree {
struct xt_node *root;
cmp_t *cmp;
struct xt_node *(*create_node)(void *key);
void (*destroy_node)(struct xt_node *n);
};
```
接著可以看 [xtree.c](https://gist.github.com/jserv/7374b3e031ec3a81441b1b93b008c5d2#file-xtree-c) 中的各個函式
`xt_create` : 用於創建一棵 XTree
`xt_insert`: 用於新增節點，由 `__xt_find` 由 `key` 值比對該插入到哪個位置。接著將 `create_node` 的鍵值設為 `key`，判斷 `root` 是否存在，成立，則呼叫 `__xt_insert` 將節點 `n` 放置 `p` 的左或右子代，並進行 `xt_update` 來平衡二元樹。
```cpp
int xt_insert(struct xt_tree *tree, void *key)
{
struct xt_node *p = NULL;
enum xt_dir d = NONE;
struct xt_node *n = __xt_find(tree, key, &p, &d);
if (n != NULL)
return -1;
n = tree->create_node(key);
if (xt_root(tree)) {
assert(d != NONE);
__xt_insert(&xt_root(tree), p, n, d);
} else
xt_root(tree) = n;
return 0;
}
```
`xt_update`: 先呼叫 `xt_balance` 判斷目前是否平衡，並進行 `rotation` 平衡後，再呼叫 `xt_max_hint` 來儲存平衡後的 `hint`，再對 `parent` 做平衡。
`xt_remove` 用於刪除節點，先輸入鍵值，並查詢是否存在，接著再呼叫 `__xt_remove` 來進行刪除。
```cpp
int xt_remove(struct xt_tree *tree, void *key)
{
struct xt_node *n = xt_find(tree, key);
if (!n)
return -1;
__xt_remove(&xt_root(tree), n);
tree->destroy_node(n);
return 0;
}
```
在 `__xt_remove` 中，若此刪除節點有右子樹，則會將右子樹中的最小值進行替換
```
7 8
/ \ remove 7 / \
5 10 -----> 5 10
/ \ / \ / \ / \
4 6 8 11 4 6 7 11
```
若沒有右子樹，則找出左子樹的最大值進行替換
```
7 6
/ remove 7 /
5 -----> 5
/ \ / \
4 6 4 7
```
再進行 `xt_update` 來完成平衡。
```cpp
static void __xt_remove(struct xt_node **root, struct xt_node *del)
{ ///// if the node has right child
if (xt_right(del)) {
/////find the least node
struct xt_node *least = xt_first(xt_right(del));
if (del == *root)
*root = least;
xt_replace_right(del, (AAAA = least) );
xt_update(root, (BBBB = xt_right(del)) );
return;
}
/////if the node has left child
if (xt_left(del)) {
////find the biggest node in left child
struct xt_node *most = xt_last(xt_left(del));
if (del == *root)
*root = most;
xt_replace_left(del, (CCCC = most) );
xt_update(root, (DDDD = xt_left(most) );
return;
}
if (del == *root) {
*root = 0;
return;
}
/* empty node */
struct xt_node *parent = xt_parent(del);
if (xt_left(parent) == del)
xt_left(parent) = 0;
else
xt_right(parent) = 0;
xt_update((EEEE = root), (FFFF = parent ));
}
```
Sign in with Wallet
Connect another wallet