# 2024q1 Homework4 (quiz3+4)
contributed by < `chloe0919` >
## 第三週測驗題
### 測驗一
測驗中的程式碼透過多項式乘法逐步拆解,例如:
$$(a+b+c)^2 = a^2+b^2+c^2+2(ab+bc+ac)$$
也就是將 $a^2$、$b^2$、$c^2$ 相加後再加上兩倍的兩兩元素相乘
現在逐步將 $N^2$ 展開 :
$$
\begin{split}
N^2 & = (a_n+a_{n-1}+a_{n-2}+...+a_0)^2,a_m=2^mor \quad a_m=0 (對應第 m 個bit是1還是0)
\\
& =a_n^2+a_{n-1}^2+a_{n-2}^2+...+a_0^2
+2[a_n(a_{n-1}+a_{n-2}+...+a_0)+a_{n-1}(a_{n-2}+a_{n-3}+...+a_0)+...+a_1a_0]
\\
& =a_n^2+a_{n-1}^2+a_{n-2}^2+...+a_0^2
+2[a_na_{n-1}+(a_n+a_{n-1})a_{n-2}+(a_n+a_{n-1}+a_{n-2})a_{n-3}+...+(a_n+a_{n-1}...a_2+a_1)a_0]
\\
& =a_0^2+\{a_n^2+a_{n-1}^2+a_{n-2}^2+...+a_1^2
+2[a_na_{n-1}+(a_n+a_{n-1})a_{n-2}+...+(a_n+a_{n-1}...+a_2)a_1]\}+2(a_n+a_{n-1}+...+a_2+a_1)a_0
\\
& =a_0^2+2P_1a_0+P_1^2
\\
\\
P_m& =a_n+a_{n-1}+...+a_m
\end{split}
$$
所以最後可以獲得以下式子,也就是 $P_0$ 即為所求平方根
$$
\begin{split}
N^2 & = a_0^2+2P_1a_0+P_1^2 = (a_0+P_1)^2 = P_0^2
\end{split}
$$
而且 $$ P_m=P_{m+1} + a_m $$
因為目標是要求出 $P_0$,也就是從 $a_0$ 累加到 $a_n$,所以要試出從 $m=n$ 到 $m=0$ 所有的 $a_m$,$a_m = 2^m$ 還是 $a_m=0$, 只要計算出 $N^2-P_m^2$ 即求得解,但因為運算成本過高,所以提供了一個想法是利用以下推導將此輪的差值 ($N^2-P_m^2$) 用上一輪的差值 ($X_{m+1}$) 和 $Y_m$ 表示
$$
\begin{split}
N^2-P_m^2 & = N^2-(P_{m+1}+a_m)^2
\\
& = N^2-(P_{m+1}^2+2P_{m+1}a_m+a_m^2)
\\
& =(N^2-P_{m+1}^2)+(2P_{m+1}a_m+a_m^2)
\\
&=X_{m+1} - (2P_{m+1}a_m+a_m^2)
\\
&= X_{m+1} - Y_m
\end{split}
$$
接下來將 $Y_m$ 拆解成 $c_m$ 和 $d_m$,$c_m = 2P_{m+1}a_m$,$d_m = a_m^2$
$$
Y_m =
\begin{cases}
c_m+d_m \quad &if\quad a_m = 2^m
\\
0 \quad &if\quad a_m = 0
\end{cases}
$$
藉由位元運算可以從 $c_m$,$d_m$ 推出下一輪的 $c_{m-1}$,$d_{m-1}$
(假設 $a_{m-1} = 2^{m-1}$,$c_m = 2(2^m)P_{m+1}$,$d_m = (2^m)^2$)
* $$
c_{m-1} = 2P_ma_{m-1} = 2(2^{m-1})P_m = 2^mP_m = 2^m(P_{m+1}+a_m) =
\begin{cases}
c_m/2+d_m \quad &if \quad a_m=2^m
\\
c_m/2 \quad &if \quad a_m=0
\end{cases}
$$
* $d_{m-1} = a_{m-1}^2 = (2^{m-1})^2 = d_m/4$
綜合上述的推導後,演算法求 $P_0$ 從 $a_n$ 開始算到 $a_0$,所以需要算到 $c_{-1}$,即為所求的 $N$
將程式碼改為用以上推導使用到的變數做表示,因為首項 $d_n = (2^n)^2 = 2^{2n}$ 所以先用 `(31 - __builtin_clz(n)` 算出 MSB,再來 `& ~1UL` 是為了要把將 MSB 扣回最接近的偶數 ($2n$),最後再向左移 $2n$ 個 bit 算出 $2^{2n}$,而且每次迭代都會將 `d >>= 2` 對應到 $d_m/4$ 的數學公式
`y = c + d` 對應到的是執行 $c_m/2+d_m$ 後算出 $y$,而每次迭代也會執行 `c >>= 1` 將 $c$ 除以 2
```c
int i_sqrt(int n)
{
if (n <= 1) /* Assume x is always positive */
return n;
int c = 0;
for (int d = 1UL << ((31 - __builtin_clz(n)) & ~1UL); d; d >>= 2) {
int y = c + d;
c >>= 1;
if (n>= y)
n -= y, c += d;
}
return c;
}
```
舉例 $n = 36$
| n | d | c | y | c(after shift) | n>=y | n | c(after add) |
| --- | --- | --- | --- | -------------- | ---- | --- | ------------ |
| 36 | 16 | 0 | 16 | 0 | True | 20 | 16 |
| 20 | 4 | 16 | 20 | 8 | True | 0 | 12 |
| 0 | 1 | 12 | 13 | 6 | False | - | - |
最後 return 6
### 測驗二
```c
carry = tmp / 10;
tmp = tmp - carry * 10;
```
根據除法原理:$f=g \times Q + r$,其中
* $f$: 被除數 (對應上述程式碼中的 $tmp$)
* $g$: 除數 (對應上述程式碼中的 $10$)
* $Q$: 商 (對應上述程式碼中的 $carry$)
* $r$: 餘數 (對應上述程式碼中的 $tmp-carry*10$)
採用 bitwise operation 來實作除法,但是因為 10 無法直接用 2 的冪項來表示,而且在以下程式碼中因為 `b[i]` 和 `a[i]` 都是介於 0-9 的個位數,`carry` 則為 0 or 1,所以 `tmp` 的範圍介於 `19`~`0`
```c
int tmp = (b[i] - '0') + (a[i] - '0') + carry;
```
因此為了得到精確的 $tmp / 10$,要先找出除數 ($x$) 的範圍,假設需要計算的被除數 ($n$) 為 $n=ab$ ($a$ 為十位數 $b$ 為個位數),且 $l=cd$ ($c$ 為十位數 $d$ 為個位數) 為比 $n$ 小的非負整數,我們可以得到此不等式:
$$
a.b \le \frac{n}{x} \le a.b9 \Rightarrow \frac{n}{a.b9} \le x \le \frac{n}{a.b}
$$
因為 $\frac{n}{a.b} = \frac{ab}{a.b} = 10$,所以一定在精確度內
再來需要證明左邊的式子成立:
$$
\frac{n}{a.b9} \le x \Rightarrow \frac{a.b9}{n} \ge \frac{1}{x}
$$
一樣先假設
$$
c.d \le \frac{l}{x} \le c.d9
$$
因為需要證明的是 $\frac{1}{x}$ 最大只有到 $\frac{a.b9}{n}$,所以假設 $\frac{1}{x} = \frac{a.b9}{n}$ 代入上述式子中也會成立:
$$
c.d \le l \times \frac{a.b9}{n} \le c.d9 \Rightarrow c.d \le cd \times \frac{a.b9}{ab} \le c.d9
$$
接下來分別去證明左右兩邊的假設皆成立:
* $c.d \le cd \times \frac{a.b9}{ab} \Rightarrow 0.1 \le \frac{a.b9}{ab} \Rightarrow \frac{a.b}{ab} \le \frac{a.b9}{ab}$ 成立
* $cd \times \frac{a.b9}{ab} \le c.d9 \Rightarrow \frac{cd(0.1\times ab+0.09)}{ab} \le c.d9 \Rightarrow c.d +\frac{0.09 \times cd}{ab}\le c.d+0.09$ 因為 $ab>cd$ 所以也成立
最後可以證明出剛開始的假設是成立的,而且只需要針對 `19` 討論,代表 $x$ 在此範圍間得到的 `tmp / x` 直到小數點後第一位都是精確的:
$$
a.b \le \frac{n}{x} \le a.b9\Rightarrow 1.9 \le \frac{19}{x} \le 1.99 \Rightarrow 9.55 \le x \le 10
$$
因為10 無法直接用 2 的冪項來表示,所以需要找到接近 $\frac{1}{10}$ 的數字 $\frac{a}{2^N}$ 來表示,此測驗中挑選 $2^N = 128$,$a=13$,計算出來$\frac{128}{13}\approx9.84$ 有落在上述範圍內
接下來透過 bitwise operation 進行計算,首先執行 `tmp >> 3) + (tmp >> 1) + tmp` 計算出 $\frac{tmp}{8}+\frac{tmp}{2}+tmp = \frac{13 \times tmp}{8}$,之後再 `<< 3` 算出 $13 \times tmp$,因為直行右移的動作,所以需要將後面被移掉的補回去 (`d0`、`d1`、`d2`),最後再往右 7 bit 算出 $\frac{13 \times tmp}{128}$,也就是 $tmp \div 10$ 的結果
```c
d0 = q & 0b1;
d1 = q & 0b11;
d2 = q & 0b111;
q = ((((tmp >> 3) + (tmp >> 1) + tmp) << 3) + d0 + d1 + d2) >> 7;
```
而餘數的計算方式為 $tmp-q \times 10$,也就是先 `(q << 2) + q)` 計算出 $q \times 5$ 後再左移 1 bit 得到 $q \times 10$
```c
r = tmp - (((q << 2) + q) << 1);
```
#### 解釋包裝成函式後
一開始 `(in | 1) - (in >> 2)` 是再計算 $\frac{3}{4} \times in$,而且如果 `in` 為偶數需要加一,但仍不清楚為何需要此操作,之後 `q = (x >> 4) + x` 是要計算 $\frac{3}{64} \times in + \frac{3}{4} \times in \approx 0.79\times in$,為了使算出來的 $0.79$ 更接近 $0.8$ 接下來連續做了 4 次逼近,最後得到一個非常接近 $0.8$的值,再將 $\frac{8}{10} \times in$ 除以 8 後得到 $\frac{1}{10} \times in$
此外使用 `q & ~0x7` 保留除了最後三個 bits 以外的位元,也就是相當於執行 `*div << 3`,例如 $q = 00110011$,$*div = 00000110$,而 `q & ~0x7` 為 $00110000$,相當於 `*div << 3` = $00110000$。`(q & ~0x7) + (*div << 1)` 就是在算出 `*div * 10`,`*mod` 會等於 $in - *div * 10$ 的結果
```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));
}
```
### 測驗三
#### `ilog2`
`i` 每次迭代都會往右移 1 bit,直到找到最高位元的位置即結束,最後計算出來的 `log` 也就是 `i` 以 2 為底的對數
```c
int ilog2(int i)
{
int log = -1;
while (i) {
i >>= 1;
log++;
}
return log;
}
```
#### `size_t ilog2`
此操作是從 $x =2^{16}$ 開始一直逼近到 $2^1$,假設有大於 $x$ 的情況就進行累加後再右移
```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;
}
```
#### 運用 GNU extension
`__builtin_clz(v)` 可以計算出 v 的最高位左邊連續 0 的個數
```c
int ilog32(uint32_t v)
{
return (31 - __builtin_clz(v));
}
```
### 測驗四
> 參考 [stevendd543](https://hackmd.io/@stevennnnnn/linux2024-homework4#%E7%AC%AC%E5%9B%9B%E9%A1%8C)
>
**Exponentially Weighted Moving Average**
是一種取平均的統計手法,採取 Exponential 主要是為了使越舊的歷史資料的權重更低,以下為計算公式:
$$
S_t =
\begin{cases}
Y_0 \quad\quad &t=0
\\
\alpha Y_t + (1-\alpha) \cdot S_{t-1} \quad\quad &t>0
\end{cases}
$$
* $\alpha$ 表示歷史資料加權降低的程度,介在 0 ~ 1 之間,越高的 $\alpha$ 會使歷史資料減少的越快
* $Y_t$ 表示在時間 $t$ 時的資料點
* $S_t$ 表示在時間 $t$ 時計算出的 EWMA
用以下函式算出 `n` 是否為 2 的冪,實際操作是去看兩種情況皆滿足時則 return True :
* `n` 不為 0
* `n` 和 `n-1` 做 and 運算時為 0
其中,只有當 `n` 為 2 的冪時,`n-1` 才會是除了最高位為 0 其他皆 1 的 情況,所以 `n & (n - 1) ` 就會是 0,例如 $n=8=1000$,$n=7=0111$,`n & (n - 1)` 為 0
```c
static inline int is_power_of_2(unsigned long n)
{
return (n != 0 && ((n & (n - 1)) == 0));
}
```
**struct ewma**
```c
struct ewma {
unsigned long internal;
unsigned long factor;
unsigned long weight;
};
```
* internal : 表示在時間 $t$ 時計算出的 EWMA,也就是 $S_t$
* factor : scaling factor,用來轉換成定點數時所需的縮放位元
* weight : exponential weight,但這邊儲存的是 $\alpha$ 的倒數,也就是 $\alpha^{\prime}$
**ewma_init**
```c
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;
}
```
因為 `weight` 和 `factor` 都要是 2 的冪,所以要先用 `is_power_of_2` 判斷,接下來要將 `weight` 和 `factor` 都分別使用 `ilog2`,也就是只儲存指數的部分
**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;
}
```
接下來 `ewma_add` 主要的操作是假設 $\alpha^{\prime} \alpha = 1$,並且利用 $\alpha^{\prime}$ 將公式改寫成 :
$$
\begin{split}
\alpha^{\prime}S_t &=
\alpha^{\prime}(\alpha Y_t + (1-\alpha) \cdot S_{t-1}) \quad\quad t>0
\\
&=
\alpha^{\prime}\alpha Y_t + \alpha^{\prime}S_{t-1}-\alpha^{\prime}\alpha S_{t-1} \quad\quad t>0
\\
&=Y_t+\alpha^{\prime}S_{t-1}-S_{t-1}
\end{split}
$$
`avg->internal` 對應到的公式是 $S_{t-1}$ (還沒計算出 $S_{t}$,所以目前取得的是 $S_{t-1}$ ),`avg->weight` 對應到的是 $\alpha^{\prime}$ 的指數部分,這是為了進行 bitwise operation
操作對應如下,特別要注意的是 `val` 需要先進行定點數的轉換,也就是需要左移 `avg->factor` (scaling factor) 個 bit:
| $\alpha^{\prime}S_{t-1}$| $S_{t-1}$ | $Y_t$ |
| -------- | -------- | -------- |
| (avg->internal << avg->weight) | avg->internal | val << avg->factor |
而最後將整個公式 `>> avg->weight`,就能算出 $S_t$
### 測驗五
此測驗和測驗三不同的是這邊已知輸入必為 `uint32_t`,所以不需要和測驗三一樣採用 while 去做判斷,其中 `r = (x > 0xFFFF) << 4` 會先去計算出 `x` 大於 $2^{16}-1$ 時需要移動的 bit 指數的值,若 `x` 大於 $2^{16}-1$ 則 `(x > 0xFFFF)` 會為 1 並且將 1 往右移 4 bits,也就是乘上 $2^4=16$,最後將 `x` 往右 16($r$) bits,其中 `r |= shift` 對應到的是 `result += 8` 的操作
最後 return 會判斷 `x>1` 是因為目標要計算的是 log 向上取整數的值,所以要避免忽略掉 `x==0x3` 和 `x==0x2` 的情況
而這邊會先將 `x--` 再去做判斷,是因為如果不加的話結果會多一
```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;
}
```
## 第四週測驗題
### 測驗一
popcount 用來計算數值的二進位表示中位元是 1 的個數,假設 x = $b_{31}b_{30}b_{29}...b_{1}b_0$,可以得知 popcount 的數學式為 : $$popcount(x) = \sum_{n=0}^{31} b_n$$
假設每次都將 x 往右移 1 bit :
$$\begin{split}
x &\rightarrow 2^{31}b_{31}+2^{30}b_{30}+...+2^{2}b_{2}+2^{1}b_{1}+b_0
\\
x>>1 &\rightarrow 2^{30}b_{31}+2^{29}b_{30}+...+2^{1}b_{2}+2^{0}b_{1}
\\
x>>2 &\rightarrow 2^{29}b_{31}+2^{28}b_{30}+...+2^{0}b_{2}
\\
.
\\
.
\\
.
\\
x>>31 &\rightarrow 2^0b_{31}
\end{split}$$
相減後,可以整理出以下數學式 :
$$\begin{split}
popcount(x)&=x-\lfloor \cfrac{x}{2} \rfloor- \lfloor \cfrac{x}{4} \rfloor - \cdots - \lfloor \cfrac{x}{2^{31}} \rfloor
\\
&= (2^{31}-2^{30}-2^{29}-...-2^0)b_{31}+(2^{30}-2^{29}-2^{28}-...-2^0)b_{30}+ ... +(2^1-2^0)b_1+b_0
\\
&= \sum_{n=0}^{31}(2^n-\sum_{i=0}^{n-1}2^i)b_n
\end{split}$$
又因為 $2^n-\sum_{i=0}^{n-1}2^i$ 為 1,所以可以用 $x-\lfloor \cfrac{x}{2} \rfloor- \lfloor \cfrac{x}{4} \rfloor - \cdots - \lfloor \cfrac{x}{2^{31}} \rfloor$ 去計算 popcount
**popcount_branchless** 原理 :
為了避免檢查 32 次,`popcount_branchless` 的實作會以 4 個位元為單位並去計算 $x-\lfloor \cfrac{x}{2} \rfloor- \lfloor \cfrac{x}{4} \rfloor - \lfloor \cfrac{x}{8} \rfloor$
,最後使用 `v = (v + (v >> 4)) & 0x0F0F0F0F` 將每 4 位元中的 set bit 的個數全部加到 byte 中
**Hamming distance** :
`A ^ B` 會計算出 A 和 B 二進位的每個位元的差,之後再以 `__builtin_popcount(A ^ B)` 就能計算出 `A ^ B` 總共的 set bit 個數
```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;
}
```
針對 [Total Hamming Distance](https://leetcode.com/problems/total-hamming-distance/description/),要計算出 `nums` 中所有的 Hamming Distance,所以要以兩個迴圈去計算出整個陣列的 Hamming Distance,最後 `return total >> 1` 是因為多計算了一次
#### Total Hamming Distance 改進
利用 `mask` 表示每次需要取出的 bit 的遮罩,並且在 for 迴圈中計算出 `nums` 中元素該 bit 所有為 1 的個數,並且 `total+=(one*(numsSize-one))` 是去計算為 1 和為 0 的個數相乘後的結果再累加
分析時間:
上述方法會使用到兩個 for 迴圈對 `nums` 中任兩個元素作計算,還有每次都要計算 `nums[i] ^ nums[j]` 也就是需要 $n$ 的位元長度,因此時間複雜度超過 $O(n*numsSize^2)$
採用改進方法只需要時間複雜度 $O(n*numsSize)$
(其中 $n$ 為 32 bits)
```c
int totalHammingDistance(int* nums, int numsSize) {
int total=0;
__uint32_t mask=1;
while(mask!=0){
int one = 0;
for(int i = 0 ;i<numsSize;i++){
__uint32_t bit = nums[i] & mask;
if(bit) one++;
}
total+=(one*(numsSize-one));
mask <<= 1;
}
return total;
}
```
### 測驗二
**Remainder by Summing digits**
當 n 以二進位進行表示時,可以寫成 $b_{n-1}b_{n-2}...b_1b_0$
並且可以利用以下推導將 n 改寫
$$
2^k \equiv
\begin{cases}
1(mod\;3), &k : even
\\
-1(mod\;3), &k : odd
\end{cases}
$$
改寫後的 n 為以下的表達式,也就是將 n 的**偶數位元相加 - 奇數位元相加** ,對應程式碼也就是 `n = popcount(n & 0x55555555) - popcount(n & 0xAAAAAAAA)`。
$$\begin{split}
n &= b_{n-1} \cdot 2^{n-1}+b_{n-2} \cdot 2^{n-2}+...+b_2 \cdot 2^2 + b_1 \cdot 2^1 + b_0
\\
&\equiv ...-b_3+b_2-b_1+b_0(mod \; 3)
\end{split}
$$
接下來為了將程式碼進一步簡化,用以下的定理進行推導(其中 $x$ 代表 $n$,而 $m$ 為 $0xAAAAAAAA$,$\overline{m}$ 為 $0x55555555$) :
* 首先可以將表達式中的**奇數位元相加**改寫成以下,也就是先將 $x$ 的奇數位元先視為全 1 後,再扣掉 $x$ 實際奇數位元為 0 的總數(以 $k$ 代表),最後再整理式子得到 :
$$\begin{split}
&popcount(x \And \overline{m}) - popcount(x \And m)
\\
&=popcount(x \And \overline{m})-(popcount(m)-k)
\\
&=popcount(x \And \overline{m})+k) - popcount(m)
\end{split}$$
* 可以發現其中 $popcount(x \And \overline{m})+k$ 也就是 **x 的偶數位元為 1 的總數 + x 奇數位元為 0 的總數**,並且觀察 $xor$ 的真值表可以發現當 $B$ 為 0 時 $A \oplus B = A \oplus 0 =A$ 而 $B$ 為 1 時 $A \oplus B = A \oplus 1 =\overline{A}$
| $A$ | $B$ | $A \oplus B$ |
| --- | --- | ------------ |
| 0 | 0 | 0 |
| 0 | 1 | 1 |
| 1 | 0 | 1 |
| 1 | 1 | 0 |
* 所以計算 $x \oplus m= x \oplus 0xAAAAAAAA$ 時,偶數位元都會是 $x$ 的原始偶數位元,奇數位元則是 $x$ 的原始奇數位元 $NOT$ 後的結果,$popcount(x \oplus m)$ 也就是計算 **x 的偶數位元為 1 的總數 + x 奇數位元為 0 的總數**,所以可以得到以下等式:
$$\begin{split}
&(popcount(x \And \overline{m})+k) - popcount(m)
\\
&=popcount(x \oplus m) - popcount(m)
\end{split}$$
最後可以將程式碼改寫成 `n = popcount(n ^ 0xAAAAAAAA) - 16`
##### 程式碼實現
###### 方法一
將上面推導出來的數學式運用後會發現 `n = popcount(n ^ 0xAAAAAAAA) - 16` 計算出來的結果範圍在 $-16 \le n \le 16$ 之間,在此範例中為了避免 $n (mod \; 3)$ 算出來的結果為負數,所以我們要將它加上一個足夠大的 3 的倍數(例如 39),現在 $n$ 的範圍就會落在 $23 \le n \le 55$ 之間,再經過一個 $mod \; 3$ 後(這次只需要 6 bit,因為 55 只需要 6 bit 就能表示了),$n$ 的範圍為 $2 \le n \le 4$,最後 `-3` 是將範圍限定在 $-1 \le n \le 1$。
但因為 $-1 \equiv 2(mod \;3)$,所以最後我們想將範圍限定在 0 到 2 之間,當 $n<0$ 經過 `(n >> 31) & 3` 後會為 3,$n \ge 0$ 則是會得到 0,所以最後就能達成 $0 \le n \le 2$
```c
int mod3(unsigned n)
{
n = popcount(n ^ 0xAAAAAAAA) + 23;
n = popcount(n ^ 0x2A) - 3;
return n + ((n >> 31) & 3);
}
```
###### 方法二 - look up table
算出 `k = popcount(n ^ 0xAAAAAAAA)` 後會得到 $k$ 的範圍 : $0 \le k \le 32$,也就是對應到 table 中的 33 個元素,例如 $k=0$ 時 $n(mod \;3) = k-16 = -16 \equiv 2 (mod \;3)$
```c
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];
}
```
##### tictactoe 程式碼理解
> [tictactoe.c](https://gist.github.com/jserv/9088f6f72a35a1fcaaf1c932399a5624)
利用 `move_masks` 定義棋盤中九個位置各自的連線狀態,視每四個 bit 為一組,分別表達該位置在連線中的順序,共 8 組則表達 8 條連線。
```c
static const uint32_t move_masks[9] = {
0x40040040, 0x20004000, 0x10000404, 0x04020000, 0x02002022,
0x01000200, 0x00410001, 0x00201000, 0x00100110,
};
```
其中,32 bit :
* $b_{31}b_{30}b_{29}b_{28}$ 代表左上到右上的連線
* $b_{27}b_{26}b_{25}b_{24}$ 代表左中到右中的連線
* $b_{23}b_{22}b_{21}b_{20}$ 代表左下到右下的連線
* $b_{19}b_{18}b_{17}b_{16}$ 代表左上到左下的連線
* $b_{15}b_{14}b_{13}b_{12}$ 代表中上到中下的連線
* $b_{11}b_{10}b_{9}b_{8}$ 代表右上到右下的連線
* $b_{7}b_{6}b_{5}b_{4}$ 代表左上到右下的連線
* $b_{3}b_{2}b_{1}b_{0}$ 代表左下到右上的連線
以 `0x40040040` 為例,0x40040040 = 0100 0000 0000 0100 0000 0000 0100 0000 表示和棋盤中**第一個位置**有關的三種連線方式,包括第一條、第四條、第七條,其中為 0100 是因為 a 位於個個連線中的順序一:
第一條 :
| | 1 | 2 | 3 |
|:----- | ----- | ----- |:----- |
| **1** | ==a== | ==b== | ==c== |
| **2** | d | e | f |
| **3** | g | h | i |
第四條 :
| | 1 | 2 | 3 |
|:----- | ----- | ----- |:----- |
| **1** | ==a== | ==b== | ==c== |
| **2** | d | e | f |
| **3** | g | h | i |
第七條 :
| | 1 | 2 | 3 |
|:----- | ----- | ----- |:----- |
| **1** | ==a== | ==b== | ==c== |
| **2** | d | e | f |
| **3** | g | h | i |
之後會以 `board |= move_masks[move];` 取得了選定走法 move 對應的連線狀態,然後更新玩家的棋盤狀態,由此可見當任四個 bit 出現 0111 時,代表任一條連線的三個位置都有棋子,也就是玩家勝利,所以利用 `player_board + 0x11111111` 將任四個 bit 中的 0111 轉成 1000 判斷是否玩家勝利
```c
static inline uint32_t is_win(uint32_t player_board)
{
return (player_board + 0x11111111) & 0x88888888;
}
```