# Bit Twiddling Hacks 解析 (二)
下一篇：[Bits Twiddling Hacks 解析 (三)](https://hackmd.io/@0xff07/WRYYYYYYYYYY)
上一篇：[Bits Twiddling Hacks 解析 (一)](https://hackmd.io/@d4cWS3kPSNiNdRbaqbKmqQ/ORAORAORAORA)
[TOC]
## 不用位元計算的技巧
這種技巧主要是建立在算術之上，類似的比如國小計算除以 9 的餘數，只要把每個位數加起來再除以 9 就好，不用真的做一次除法。
下面這邊有一些技巧是建立在類似的觀念上，僅利用加減乘除，沒有用位元運算。在十進位中大部分都有對應的做法，所以可以借助已經熟悉的十進位進行理解。
這種做法可以讓運算子看起來很少，而且有種到處都是 magic number 的感覺，比如 [3 個運算計算 14 位元整數的 population count](https://graphics.stanford.edu/~seander/bithacks.html#CountBitsSet64) 其中一個解法就有用到。
### 用乘法重複一組數字
先看十進位的例子：
$$
1,001,001,001,001,001 \times 987 = 987,987,987,987,987,987
$$
主要就是某個數字，乘上重複的 $100...001$ ，就會把本來的數字重複很多次。不管是寫成科學記號或是用直式乘法來證明都很容易看出來這件事。
接下來看看一個 16 進位的版本：
$$
\mathtt{res = 0x200040008001ULL} \times \mathtt{0xBAD}
$$
這個結果寫成 16 進位的話，會是 `0x175a2eb45d68bad` ，乍看之下不是很好懂。但如果用直式乘法去思考的話，其實會是：
```c=
res = (0xBAD) +
(0xBAD << 15) +
(0xBAD << 30) +
(0xBAD << 45);
```
原理也是一樣。如果把上面的東西寫成二進位的話，其實是在做：
```c=
mask = 0b1000000000000001000000000000001000000000000001;
bad = 0b101110101101; // 0xBAD
res = mask * bad;
```
因為 64 位元的關係，把完整的數字寫出來實在是很長，所以大部分都會寫 16 進位，然後就會變成一個有點令人費解的樣子。
### 用乘法壓縮位元
這個用直式乘法非常好證明。舉例來說，現在有一個 64 位元無號數的 ==2 進位== 形式長這樣：
```c
tmp = 0000000H 000000G0 00000F00 0000E000 000D0000 00C00000 0B000000 A0000000
```
其中 A、B、C、D、E、F、G、H 都是 `0` 或 `1`。如果做：
```c
tmp = (tmp << 0) +
(tmp << 8) +
(tmp << 16) +
(tmp << 24) +
(tmp << 32) +
(tmp << 40) +
(tmp << 48) +
(tmp << 56);
```
那可以發現計算結果最高位的 8 個位元，會是 (觀察對角線方向)：
```c
0000000H
000000G0
00000F00
0000E000
000D0000
00C00000
0B000000
A0000000
---------(+)
ABCDEFGH
```
回到剛剛的計算：那一連串的位移其實可以簡化成：
```c
tmp = tmp * 0x0101010101010101ULL;
```
所以總結以上，做：
```c
res = (tmp * 0x0101010101010101ULL) >> 56;
```
就會有把稀疏的位元壓縮再一起的效果。
### 用取餘數計算位數和
有一些問題的思路可以從「計算位數和」下手。比如說 population count 就是在算每個位數的和; 而 parity bit 就是算完位數和之後，再看奇偶性 (比如做 `x & 1`)。
==暖身：以 10 為基數的位數和==
這件事情的原理，跟國小時判斷一個數是不是 9 的倍數，使用「(位數和) % 9，看看餘數是不是 0」，的原理是很像的。比如說想要計算：
$$
9487 \div 9
$$
的餘數，那麼小時候大家都學過，去計算：
$$
(9 + 4 + 8 + 7) \div 9
$$
的餘數就好。這是因為對於任意非負整數 $k$，有：
$$
10^k = \underbrace{999...9}_{k} + 1
$$
因此，假定一個數的 10 進位表示法是 $a_k \dots a_2a_1a_0$，或是說：
$$
n = \sum_{i = 0}^{k}a_i\times 10^i
$$
那麼可以把所有的 $10^i$ 拆成 $(\underbrace{999...9}_{i} + 1)$，然後經過一點計算之後：
$$
\begin{align}
n &= \sum_{i = 0}^{k}a_i\times (\underbrace{999\dots 9 + 1}_{i}) & \mod 9\newline
&= \sum_{i = 0}^{k}a_i\times (\underbrace{999\dots 9}_{i}) + (a_i \times 1) & \mod 9 \newline
&= \sum_{i = 0}^{k} a_i & \mod 9
\end{align}$$
但如果==位數的總和比 9 小==，那麼有 `mod` 跟沒有 `mod` 結果是一樣的，那麼除法計算的商根本就是 0。也就是說 $\text{mod}$ 根本可以拿掉，像這樣：
$$
\begin{align}
n &= \sum_{i = 0}^{k}a_i \mod 9\newline
&= \sum_{i = 0}^{k}a_i \quad \left(\text{if }\sum_{i = 0}^{k}a_i < 9\right)
\end{align}
$$
比如說，想計算 $112211$ 的位數和是 $1 + 1 + 2 + 2 + 1 + 1 = 8$，剛好就會是它除以 $9$ 的餘數。
==例子：16 進位的位數和==
現在我們把基數從 $10$ 換成 $16$。比如說，想要計算 `n = 0x11101111` 的位數和。這個數字其實是：
$$
n = \sum_{i = 0}^{7} b_i \times 16^i
$$
其中，裡面任意 $b_i$ 只會是 $0$ 或 $1$。仿照上面的拆法。不過因為現在 base 是 16，所以就把 $16^{i}$ 拆成$\mathtt{0x\underbrace{FF...FF}_{i}} + 1$
$$
\begin{align}
n &= \sum_{i = 0}^{7} b_i \times 16^i &\mod 15\newline
&= \sum_{i = 0}^{7} b_i \times (\mathtt{0x\underbrace{FF...FF}_{i}} + 1) &\mod \mathtt{0xF}\newline
&= \sum_{i = 0}^{7} b_i &\mod \mathtt{0xF} \newline
&= \sum_{i = 0}^{7} b_i & \text{if } \sum_{i = 0}^7{b_i} < \mathtt{0xF}
\end{align}
$$
所以可以一眼知道位數和就是每個位元加起來，就是 $7$。
### 小提示
:::info
後面有很多利用乘法的招式，大致上都是遵照下面的形式：
1. 在 64 位元中，用乘法重複一個 byte 的內容很多次。
2. 依照特定間隔取出位元，使得該 byte 中每個位元恰好在 mask 後的結果出現一次。
3. 用取餘數計算所有位元的加總(比如 population count 中的 [方法四（之一）：加總位數和 (14 位元數值 + 64 位元運算)](https://hackmd.io/@d4cWS3kPSNiNdRbaqbKmqQ/MUDAMUDAMUDA#方法四（之一）：加總位數和-14-位元數值--64-位元運算); 或用上面的乘法把所有位元壓縮在一起(比如反轉所有位元中的[方法四：用乘法反轉一個-byte](https://hackmd.io/@d4cWS3kPSNiNdRbaqbKmqQ/MUDAMUDAMUDA#方法四：用乘法反轉一個-byte))
舉例來說下面三者：
```c=
/* 對一個 byte 做 population count */
(b * 0x200040008001ULL & 0x111111111111111ULL) % 0xf;
/* 計算一個 byte b 的 parity bit */
(((b * 0x0101010101010101ULL) & 0x8040201008040201ULL) % 0x1FF) & 1;
/* 反轉一個 byte b */
((b * 0x80200802ULL) & 0x0884422110ULL) * 0x0101010101ULL >> 32;
```
:::
## Population Count
這個經典的問題是這樣：給定一組 bit mask，要問裡面有多少個 1。最天真的方法就是用迴圈：
```c
/* @n : the integer to be counted */
/* @res: the pop count result */
while(res = 0; n ; n >>= 1) {
res += (n & 1);
}
```
不過這樣明顯會出現很多分支，所以接下來就要思考改良的版本。
### 方法一：清除所有位元
原文：[Counting bits set, Brian Kernighan's way](https://graphics.stanford.edu/~seander/bithacks.html#CountBitsSetKernighan)
```c=
for (res = 0; n; n &= n-1) {
res++;
}
```
這個是用上面「清除最低位元」的做法去做。在位元很稀疏的時候，有機會比天真的做法快。不過還是免不了使用 branch。
### 方法二：暴力查表
原文：[Counting bits set by lookup table](https://graphics.stanford.edu/~seander/bithacks.html#CountBitsSetTable)
~~就是那個傳說中把結果存在月球上的做法啊。~~
這個做法是建立一個大小是 256 的陣列 `table` 當作表，枚舉所有 8 位元的整數所有可能數值中，對應的 1 的個數。比如說我想要查 `73` 這個數字中，有多少個 1，那就去找 `table[73]`。所以這個有點大的表會像是：
```c
unsigned char table[256] = {
0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4,
1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
4, 5, 5, 6, 5, 6, 6, 7, 5, 6, 6, 7, 6, 7, 7, 8
};
res = table[n];
```
不過這樣打起來有點煩瑣，所以可以用巨集的技巧簡化：
```c=
unsigned char table[256] = {
# define B2(n) n, n+1, n+1, n+2
# define B4(n) B2(n), B2(n+1), B2(n+1), B2(n+2)
# define B6(n) B4(n), B4(n+1), B4(n+1), B4(n+2)
# define B8(n) B6(n), B6(n+1), B6(n+1), B6(n+2)
B8(0)
};
res = table[n];
```
這個巨集其實是用下面的遞迴結構：
$$
\mathtt{pop(n([N:0])}\begin{cases}
0 + \mathtt{pop(n[N-2:0])} & \text{if }\quad \mathtt{n[N:N-1] = 0} \newline
1 + \mathtt{pop(n[N -2:0])} & \text{if }\quad \mathtt{n[N:N-1] = 1} \newline
1 + \mathtt{pop(n[N-2:0])} & \text{if }\quad \mathtt{n[N:N-1] = 2} \newline
2 + \mathtt{pop(n[N-2:0])} & \text{if }\quad \mathtt{n[N:N-1] = 3}
\end{cases}
$$
原理就是照最高 2 個位元進行討論。他有可能是 `0b00`, `0b01`, `0b10`, `0b11`，而 pop count 的結果就是「最高 2 位元的 1 的數目」加上「其餘低位元的 1 的數目」。
對於其他 8 的倍數寬度的數值型別，只要多做幾次就可以了：
```c=
unsigned char table[256] = {
# define B2(n) n, n+1, n+1, n+2
# define B4(n) B2(n), B2(n+1), B2(n+1), B2(n+2)
# define B6(n) B4(n), B4(n+1), B4(n+1), B4(n+2)
# define B8(n) B6(n), B6(n+1), B6(n+1), B6(n+2)
B8(0)
};
res = table[n & 0xff] +
table[(n >> 8) & 0xff] +
table[(n >> 16) & 0xff] +
table[(n >> 24) & 0xff];
```
另外一個看起來很酷的變形是用一個 `unsigned char` 指標去指要計算的整數，然後就可以用中括號去存取每個 byte:
```c
unsigned char * p = (unsigned char *) &n;
res = table[p[0]] +
table[p[1]] +
table[p[2]] +
table[p[3]];
```
### 方法三 (之一)：Divide and Conquer
原文：[Counting bits set, in parallel](https://graphics.stanford.edu/~seander/bithacks.html#CountBitsSetParallel)
這個應該是最有名的做法吧？主要的思路是這樣：
```c
v = (v & 0x55555555) + ((v >> 1) & 0x55555555);
v = (v & 0x33333333) + ((v >> 2) & 0x33333333);
v = (v & 0x0F0F0F0F) + ((v >> 4) & 0x0F0F0F0F);
v = (v & 0x00FF00FF) + ((v >> 8) & 0x00FF00FF);
v = (v & 0x0000FFFF) + ((v >> 16) & 0x0000FFFF);
res = v;
```
但裡面給出一個等價，但是運算次數更少的版本：
```c=
c = v - ((v >> 1) & 0x55555555);
c = ((c >> 2) & 0x33333333) + (c & 0x33333333);
c = ((c >> 4) + c ) & 0x0F0F0F0F;
c = ((c >> 8) + c ) & 0x00FF00FF;
c = ((c >> 16) + c) & 0x0000FFFF;
```
第一行做法為什麼會對，可以窮舉所有可能的狀況：
```c
0b11 - (0b11 >> 1) = 2
0b10 - (0b10 >> 1) = 1
0b01 - (0b01 >> 1) = 1
0b00 - (0b00 >> 1) = 0
```
所以就剛好跟原先的第一步是一樣的。接下來的做法大上就是跟原先一樣。
### 方法三 (之二)：Divide and Conquer + 乘法
原文：[Counting bits set, in parallel](https://graphics.stanford.edu/~seander/bithacks.html#CountBitsSetParallel) 的下面
概念大致上跟上面一樣，只是後面用乘法做：
```c=
v = v - ((v >> 1) & 0x55555555);
v = (v & 0x33333333) + ((v >> 2) & 0x33333333);
r = ((v + (v >> 4) & 0xF0F0F0F) * 0x1010101) >> 24;
```
最後面乘法的作用是像下面這樣：
```c
c = v - ((v >> 1) & 0x55555555);
c = ((c >> 2) & 0x33333333) + (c & 0x33333333);
c = ((c >> 4) + c ) & 0x0F0F0F0F;
r = (c + (c << 8) + (c << 16) + (c << 24)) >> 24;
```
### 方法四（之一）：加總位數和 (14 位元數值 + 64 位元運算)
原文：[Counting bits set in 14, 24, or 32-bit words using 64-bit instructions](https://graphics.stanford.edu/~seander/bithacks.html#CountBitsSet64)
這邊使用到另外一個概念：==population count = 加總所有位數==。所以一個可能的做法是下面這樣：
```c=
int res = 0;
unsigned long long tmp = (n) + (n << 15) + (n << 30) + (n << 45);
for (int i = 0; i < 14; i++) {
res += (tmp & 1);
tmp >>= 4;
}
```
大致上就是：
1. 重複這個數字很多次
2. 利用「原數字重複的週期」跟「取位元的週期」互質來做到「恰好每個位元都出現一次」
3. 把取出來的位元加總。
不過在文章中出現的程式碼很短，是這樣：
```c
res = (n * 0x200040008001ULL & 0x111111111111111ULL) % 0xf;
```
如果可以一眼看出來在做什麼的話，那可以跳過。如果看不出來的話，下面接著解釋。這段程式只有 3 個運算，分別就代表上面的 3 個步驟：
1. `* 0x200040008001ULL`：重複原來的數字
2. `& 0x111111111111111ULL`：利用「原數字重複的週期」跟「取出位元的週期」互質來做到「恰好每個位元都出現一次」。這邊利用「原數字重複的週期」是 15; 「取出位元的週期」是 4。
3. `% 0xf`：把取出來的位元加總。
#### 1. n * 0x200040008001ULL
如果觀察一下 `0x200040008001ULL`，會發現他是不斷重複的 `0b000 0000 0000 0001`(14 個 `0`，1 個 `1`)。比如說如果數字是 `0xBAD`，那乘上 `0x200040008001ULL` 之後會變成：
```c
tmp1 = 0xBAD +
0xBAD << 15 +
0xBAD << 30 +
0xBAD << 45;
```
反正就是把 0xBAD 重複寫很多次。但關鍵就是：他會以==每 15 個位元為週期重複==。
#### 2. & 0x111111111111111ULL
如果把每個位元用 `0DCBA9876543210` 表示每個位元，(==比如第 `A` 代表第 10 個位元、`4` 代表第 4 位元等等==)的話，上面那個過程其實是「讓 `n` 中每一個位元恰好在結果中出現 1 次」。如果用 `0DCBA9876543210` 分別代表 0 ~ 15 個二進位中的位數的話，上面的步驟其實是照下面的過程取出對應的位數：
$$
\begin{align}
\mathtt{0DC\underset{\uparrow}{B}A98\underset{\uparrow}{7}654\underset{\uparrow}{3}210\ \underset{\uparrow}{0}DCB\underset{\uparrow}{A}987\underset{\uparrow}{6}543\underset{\uparrow}{2}10\ 0\underset{\uparrow}{D}CBA\underset{\uparrow}{9}876\underset{\uparrow}{5}432\underset{\uparrow}{1}0\ 0D\underset{\uparrow}{C}BA9\underset{\uparrow}{8}765\underset{\uparrow}{4}321\underset{\uparrow}{0}}\newline
\mathtt{000\underset{\uparrow}{1}000\underset{\uparrow}{1}000\underset{\uparrow}{1}000\ \underset{\uparrow}{0}000\underset{\uparrow}{1}000\underset{\uparrow}{1}000\underset{\uparrow}{1}00\ 0\underset{\uparrow}{1}000\underset{\uparrow}{1}000\underset{\uparrow}{1}000\underset{\uparrow}{1}0\ 00\underset{\uparrow}{1}000\underset{\uparrow}{1}000\underset{\uparrow}{1}000\underset{\uparrow}{1}}\newline
\hline
\mathtt{000\underset{\uparrow}{B}000\underset{\uparrow}{7}000\underset{\uparrow}{3}000\ \underset{\uparrow}{0}000\underset{\uparrow}{A}000\underset{\uparrow}{6}000\underset{\uparrow}{2}00\ 0\underset{\uparrow}{D}000\underset{\uparrow}{9}000\underset{\uparrow}{5}000\underset{\uparrow}{1}0\ 00\underset{\uparrow}{C}000\underset{\uparrow}{8}000\underset{\uparrow}{4}000\underset{\uparrow}{0}}
\end{align}
$$
也就是依照：
```
n[0], n[4], n[8], n[12],
n[1], n[5], n[9], n[13],
n[2], n[6], n[10], n[14],
n[3], n[7], n[11]
```
的順序，把第 `i` 位元取出。其中 `n[i]` 是由最低位數來第 `i` 個位元。因為想不到用比較好的表達方法，所以就借用類似 verilog 的語法來表打。
關鍵是：==4 跟 15 互質，所以這樣就保證恰好取到 `n` 裡面的所有位數==。為了方便，暫時用 `tmp2` 來稱呼到目前為止的計算結果：
$$
\mathtt{tmp2 := (n * 0x200040008001ULL\ \&\ 0x111111111111111ULL)};
$$
這個結果會是個位元形如下面這樣的數字：
$$
\mathtt{tmp2 = 0b0...\underbrace{000n_{11}000n_{7}000n_{3} \dots000n_1000n_{12}000n_8000n_4000n_0}_{4 \times 14 = 56 \text{ bits}}}
$$
現在我們有什麼？
1. 所有 $n$ 的位元都在裡面剛好出現過一次
2. 目前的計算結果是：
$$
\begin{align}
\text{tmp2} = n_0 &\cdot 16^0 \newline
n_4 &\cdot 16^1 &+ \newline
n_8 &\cdot 16^2 &+ \newline
n_{12} &\cdot 16^{3} &+ \newline
n_{1} &\cdot 16^{4} &+ \newline
n_{4} & \cdot 16^{5} &+ \newline
\dots \newline
n_{7} &\cdot 16^{12} &+ \newline
n_{11} & \cdot 16^{13}
\end{align}
$$
現在所有 `n` 的位元都恰好在 `tmp2` 中出現一次，其他位元都是 0。現在只欠東風：==把 `tmp` 的每個位元加總。==
#### 3. % 0xF
因為現在基數是 16，加上位數和最大就是 14 (14 個位元都是 1) < (16 - 1)。用前面[計算位數和](https://hackmd.io/@d4cWS3kPSNiNdRbaqbKmqQ/ORAORAORAORA#計算位數和)的結論可以知道：這時做 `% (16 - 1)` 就可以求出位數和。
### 方法四（之二）：加總位數和 (24 位元數值 + 64 位元運算)
原文：[Counting bits set in 14, 24, or 32-bit words using 64-bit instructions](https://graphics.stanford.edu/~seander/bithacks.html#CountBitsSet64)
```c
res = ((n & 0xfff) * 0x1001001001001ULL & 0x84210842108421ULL) % 0x1f;
res += (((n & 0xfff000) >> 12) * 0x1001001001001ULL & 0x84210842108421ULL)% 0x1f;
```
原理跟上面一樣，但他把「12 位元數字的，相同的構造方法，重複做 2 次。」加起來就會得到得到 24 位元的 population count。
處理 12 位元的 population count 時，重複寫原來數字對週期，從 15 位元，變成 12 位元; 以及 mask 的週期，由 4 位元變成 5 位元。其餘概念相同。
### 方法四（之三）：加總位數和 (32 位元數值 + 64 位元運算)
原文：[Counting bits set in 14, 24, or 32-bit words using 64-bit instructions](https://graphics.stanford.edu/~seander/bithacks.html#CountBitsSet64)
```c
res = ((n & 0xfff) * 0x1001001001001ULL & 0x84210842108421ULL) % 0x1f;
res += (((n & 0xfff000) >> 12) * 0x1001001001001ULL & 0x84210842108421ULL) % 0x1f;
res += ((n >> 24) * 0x1001001001001ULL & 0x84210842108421ULL) % 0x1f;
```
原理跟上面一樣，變成做 12 位元的狀況做 3 次，一次是計算 32 位元中，低 12 位的位數和; 一次是計算次高 12 位元的位數和; 最後計算最高 8 位的位數和。不過反正 12 位元可以用，把 8 位元數字當成 12 位元數字照樣可以用。所以看起來就是 12 位元的做法做 3 次。
## Parity Bit
parity bit 的做法很多會跟 population count 很像，因為能做 population count 就可以用計算的結果直接去找出 parity bit。
### 方法一：直覺的方法
原文：[Computing parity the naive way](https://graphics.stanford.edu/~seander/bithacks.html#ParityNaive)
最直覺的方法就是做：
```c=
while (x) {
parity = !parity;
x = x & (x - 1);
}
```
不過照樣子會希望避開 branch 跟 jump。做法大概就是兩類：
1. 從 popullation count 的做法改來
2. parity bit 就是想辦法計算所有位元 `XOR` 起來的值。所以往這個方向做也可以。
### 方法二：暴力查表
原文：[Compute parity by lookup table](https://graphics.stanford.edu/~seander/bithacks.html#ParityLookupTable)
這個跟 population count 那邊的查表法類似。可以建立一個大小 256 的表，枚舉 `0` 到 `255` 中所有值的 parity：
```c
char b;
static const bool P[] = {
0, 1, 1, 0, 1, 0, 0, 1,
1, 0, 0, 1, 0, 1, 1, 0,
1, 0, 0, 1, 0, 1, 1, 0,
0, 1, 1, 0, 1, 0, 0, 1,
1, 0, 0, 1, 0, 1, 1, 0,
0, 1, 1, 0, 1, 0, 0, 1,
0, 1, 1, 0, 1, 0, 0, 1,
1, 0, 0, 1, 0, 1, 1, 0,
1, 0, 0, 1, 0, 1, 1, 0,
0, 1, 1, 0, 1, 0, 0, 1,
0, 1, 1, 0, 1, 0, 0, 1,
1, 0, 0, 1, 0, 1, 1, 0,
0, 1, 1, 0, 1, 0, 0, 1,
1, 0, 0, 1, 0, 1, 1, 0,
1, 0, 0, 1, 0, 1, 1, 0,
0, 1, 1, 0, 1, 0, 0, 1
};
bool parity = P[b];
```
而這個表也有類似前面 population count 的巨集建立方式。所以程式就可以簡化成：
```c=
static const bool ParityTable256[256] =
{
# define P2(n) n, n^1, n^1, n
# define P4(n) P2(n), P2(n^1), P2(n^1), P2(n)
# define P6(n) P4(n), P4(n^1), P4(n^1), P4(n)
P6(0), P6(1), P6(1), P6(0)
};
unsigned char b; // byte value to compute the parity of
bool parity = ParityTable256[b];
```
會做一個 byte 的 parity bit，那就可以用一樣的方法去做一個 byte 的整數倍位元的 parity bit：
```c=
// Variation:
unsigned char * p = (unsigned char *) &v;
parity = ParityTable256[p[0] ^ p[1] ^ p[2] ^ p[3]];
```
或是先把狀態「壓縮」到 8 bit 裡面，再查表：
```c=
unsigned int v;
v ^= v >> 16;
v ^= v >> 8;
bool parity = ParityTable256[v & 0xff];
```
### 方法三：加總位數和 (8-bit 限定)
原文：[Compute parity of a byte using 64-bit multiply and modulus division](https://graphics.stanford.edu/~seander/bithacks.html#ParityWith64Bits)
```c=
unsigned char b; // byte value to compute the parity of
bool parity =
(((b * 0x0101010101010101ULL) & 0x8040201008040201ULL) % 0x1FF) & 1;
```
這個就是從 population count 的[方法四（之一）：加總位數和 (14 位元數值 + 64 位元運算)](https://hackmd.io/@d4cWS3kPSNiNdRbaqbKmqQ/MUDAMUDAMUDA#方法四（之一）：加總位數和-14-位元數值--64-位元運算)改來的。先做 population count，再看結果的奇偶性。
### 方法四：直式乘法
原文：[Compute parity of word with a multiply](https://graphics.stanford.edu/~seander/bithacks.html#ParityMultiply)
```c=
unsigned int v; // 32-bit word
v ^= v >> 1;
v ^= v >> 2;
v = (v & 0x11111111U) * 0x11111111U;
return (v >> 28) & 1;
```
做完 `v^= v>> 1` 跟 `v ^= v >> 2` 之後，每 0, 4, 8, 12 ... 24, 28 位元，都會剛好是 0 ~ 3, 4 ~ 7, 8 ~ 12 ... 24 ~ 27, 28 ~ 31 位元的 parity bit。
接著把這些位元取出來，用類似「直式乘法」的方式去思考：乘上 `0x11111111U` 時，第 28 位元其實是在做：
```c=
v = v +
(v << 4) +
(v << 8) +
(v << 12) +
(v << 16) +
(v << 20) +
(v << 24) +
(v << 28)
```
所以第 28 位元，就會是加總結果的個位數。因此把它取出來，就是 parity bit。
### 方法五：全部 XOR 起來
原文：[Compute parity in parallel](https://graphics.stanford.edu/~seander/bithacks.html#ParityParallel)
可以直接做 32 次。不過因為 XOR 有結合律跟交換律，所以可以用分治法，一直「對折」來找出 parity bit：
```c
unsigned int v; // word value to compute the parity of
v ^= v >> 16;
v ^= v >> 8;
v ^= v >> 4;
v ^= v >> 2;
v ^= v >> 1;
parity = v;
```
但是可以用一些 magic number 去減少次數：
```c=
unsigned int v; // word value to compute the parity of
v ^= v >> 16;
v ^= v >> 8;
v ^= v >> 4;
v &= 0xf;
return (0x6996 >> v) & 1;
```
大做法大致上跟上面差不多，不過只對折到寬度為 4 就不繼續對折了，而是用一個奇怪的數字 `0x9669` 去位移。這個數字其實是 `0110 1001 1001 0110`，就是上面的 `ParityTable256` 前 16 個數字壓縮進 16 位元裡面。這樣一來，只要位移並取出第 1 個位元，就可以做到跟前面查表一樣的效果。
## 反轉所有位元
原文：[Reverse bits the obvious way](https://graphics.stanford.edu/~seander/bithacks.html#BitReverseObvious)
問題：給定一個 `unsigned int v`，反轉 `v` 的所有位元。即結果 `c` 的第 `i` 位元，必須等於 `v` 的第 `31 - i` 位元。其中 `i` 的範圍是 `0 ~ 31`。
最直覺的做法就是用迴圈下去：
```c=
unsigned int v; // input bits to be reversed
unsigned int r = v; // r will be reversed bits of v; first get LSB of v
int s = sizeof(v) * CHAR_BIT - 1; // extra shift needed at end
for (v >>= 1; v; v >>= 1)
{
r <<= 1;
r |= v & 1;
s--;
}
r <<= s; // shift when v's highest bits are zero
```
### 方法一：暴力查表
跟前面兩個做法類似。先討論 8 位元的狀況，然後枚舉 0 ~ 255 中，每一個數值反轉之後的結果，存在一個大小 256 的陣列裡面：
```c=
unsigned char b;
static const unsigned char BitReverseTable256[256] = {
0, 128, 64, 192, 32, 160, 96, 224, 16, 144, 80, 208, 48, 176, 112, 240, 8, 136, 72, 200, 40, 168, 104, 232, 24, 152, 88, 216, 56, 184, 120, 248,
4, 132, 68, 196, 36, 164, 100, 228, 20, 148, 84, 212, 52, 180, 116, 244, 12, 140, 76, 204, 44, 172, 108, 236, 28, 156, 92, 220, 60, 188, 124, 252,
2, 130, 66, 194, 34, 162, 98, 226, 18, 146, 82, 210, 50, 178, 114, 242, 10, 138, 74, 202, 42, 170, 106, 234, 26, 154, 90, 218, 58, 186, 122, 250,
6, 134, 70, 198, 38, 166, 102, 230, 22, 150, 86, 214, 54, 182, 118, 246, 14, 142, 78, 206, 46, 174, 110, 238, 30, 158, 94, 222, 62, 190, 126, 254,
1, 129, 65, 193, 33, 161, 97, 225, 17, 145, 81, 209, 49, 177, 113, 241, 9, 137, 73, 201, 41, 169, 105, 233, 25, 153, 89, 217, 57, 185, 121, 249,
5, 133, 69, 197, 37, 165, 101, 229, 21, 149, 85, 213, 53, 181, 117, 245, 13, 141, 77, 205, 45, 173, 109, 237, 29, 157, 93, 221, 61, 189, 125, 253,
3, 131, 67, 195, 35, 163, 99, 227, 19, 147, 83, 211, 51, 179, 115, 243, 11, 139, 75, 203, 43, 171, 107, 235, 27, 155, 91, 219, 59, 187, 123, 251,
7, 135, 71, 199, 39, 167, 103, 231, 23, 151, 87, 215, 55, 183, 119, 247, 15, 143, 79, 207, 47, 175, 111, 239, 31, 159, 95, 223, 63, 191, 127, 255,
};
res = table[b];
```
類似地，這個表也可以用巨集去建立：
```c=
tatic const unsigned char BitReverseTable256[256] =
{
# define R2(n) n, n + 2*64, n + 1*64, n + 3*64
# define R4(n) R2(n), R2(n + 2*16), R2(n + 1*16), R2(n + 3*16)
# define R6(n) R4(n), R4(n + 2*4 ), R4(n + 1*4 ), R4(n + 3*4 )
R6(0), R6(2), R6(1), R6(3)
};
unsigned int v; // reverse 32-bit value, 8 bits at time
unsigned int c; // c will get v reversed
```
這個巨集運作的原理可以用這種方法思考：假定本來的 byte 的二進位形式是：
$$
\mathtt{hgfedcba}
$$
那麼這個 byte 反轉的結果，可以用下面的方法計算出來：
$$
\begin{align}
\mathtt{abcdefgh} =& \mathtt{(ab * 64)\ + } \newline
& \mathtt{(cd * 16)\ + } \newline
& \mathtt{(ef * 4 )\ + } \newline
& \mathtt{(gh * 1 )}
\end{align}
$$
這個過程可以想像成成要開下面這樣的陣列：
$$
\mathtt {int\ table[4][4][4][4];}
$$
其中：
$$
\begin{align}
\mathtt{table[hg][fe][dc][ba]} =& \mathtt{(ab * 64)\ + } \newline
& \mathtt{(cd * 16)\ + } \newline
& \mathtt{(ef * 4 )\ + } \newline
& \mathtt{(gh * 1 )}
\end{align}
$$
另外一方面，手動展開上面的巨集，就會發現：開一個 `table[4][4][4][4]` 去存取 `table[hg][fe][dc][ba]`，在 C 語言的 array subscription 下，這件事跟開一個大小 256 的表 `table[256]` ，然後直接存取 `table[hgfedcba]` 效果，是一樣的。所以就可以用上面的巨集去做。
會做 8 位元的表之後，就可以做 32 位元的反轉：
```c=
c = (BitReverseTable256[v & 0xff] << 24) |
(BitReverseTable256[(v >> 8) & 0xff] << 16) |
(BitReverseTable256[(v >> 16) & 0xff] << 8) |
(BitReverseTable256[(v >> 24) & 0xff]);
```
或是用指標：
```c=
unsigned char * p = (unsigned char *) &v;
unsigned char * q = (unsigned char *) &c;
q[3] = BitReverseTable256[p[0]];
q[2] = BitReverseTable256[p[1]];
q[1] = BitReverseTable256[p[2]];
q[0] = BitReverseTable256[p[3]];
```
### 方法二：分治法
位元逆轉其實可以用分治法去做：假設現在要反轉的是 `v[2n-1:0]`，而且現在已經有左半部的反轉結果 `left = rev(v[2n-1:n)]` 跟右半部的反轉結果 `right = rev(v[n-1:0])`，那只要把這兩半部對調，即 `{righ, left}` 就是反轉的結果。
舉例來說：假設現在想要反轉的東西是這樣：
```=
* * * * * * *
* * * * * *
* * * * *
* * * *
* * *
* *
*
```
如果已經有上下兩半反轉結果，那把這上下兩半交換，就會得到全部的反轉結果：
```=
* * * * * *
* * * * * * * *
* * * * * * * * * *
* * * * * * * * --exchange each half---> * * * *
* to get reversed seq. * * * * *
* * * * * * * *
* * * * * * * * * *
* * * * * * * * * * * *
```
然後遞迴下去：就可以知道可以用下面的順序找出 32 位元的反轉結果：
1. 找長度為 2 的反轉結果 X 16 組。接著利用上面的推論會找到
2. 長度為 4 的反轉結果 8 組，接著
3. 長度為 8 的反轉結果 4 組，接著
4. 長度為 16 的反轉結果 2 組，接著
5. 長度為 32 的反轉結果 1 組
舉實例來說，如果要反轉的位元是 `hgfe dcba`，過程會像下面這個樣子：
```
[h][g][f][e][d][c][b][a] --長度為 2 的反轉結果 4 組-->
[gh] [ef] [cd] [ab] --長度為 4 的反轉結果 2 組-->
[efgh] [abcd] --長度為 8 的反轉結果 1 組-->
[abcdefgh]
```
所以，在 32 位元的狀況下，就會看起來像下面這樣：
```c=
unsigned int v; // 32-bit word to reverse bit order
v = ((v >> 1) & 0x55555555) | ((v & 0x55555555) << 1);
v = ((v >> 2) & 0x33333333) | ((v & 0x33333333) << 2);
v = ((v >> 4) & 0x0F0F0F0F) | ((v & 0x0F0F0F0F) << 4);
v = ((v >> 8) & 0x00FF00FF) | ((v & 0x00FF00FF) << 8);
v = ( v >> 16 ) | ( v << 16);
```
另外，也有迴圈版的：
```c=
unsigned int s = 32 // bit size; must be power of 2
unsigned int mask = ~0;
while ((s >>= 1) > 0)
{
mask ^= (mask << s);
v = ((v >> s) & mask) | ((v << s) & ~mask);
}
```
這個版本差別是他用 top-down 的方法做 (16 位元 X2、8 位元 X4、4 位元 X8 ...)，並且把 `mask` 在迴圈裡面計算。
### 方法三：「64 位元乘法 + 取餘數」反轉一個 byte
跟前面的 population count 類似的技巧：
```c=
unsigned char b; // reverse this (8-bit) byte
b = (b * 0x0202020202ULL & 0x010884422010ULL) % 1023;
```
大致可以看出來，乘上 `0x0202020202` 就是要讓原本的 8 位元數字在 64 位元中「循環很多次」; 而 `% 1023` 看起來跟前面「計算位數和」中的技巧有點像。但那個 `& 0x010884422010` 這個數字就有點不好懂。不過，如果把整個過程展開來，就可以發現他的功用：
```
HGFEDCBA (b)
00000000 00000010 00000010 00000010 00000010 00000010 (0x0202020202)
------------------------------------------------------------------------- (*)
0000000H GFEDCBAH GFEDCBAH GFEDCBAH GFEDCBAH GFEDCBA0
00000001 00001000 10000100 01000010 00100000 00010000 (0x010884422010)
------------------------------------------------------------------------- (&)
0000000H 0000C000 G0000B00 0F0000A0 00E00000 000D0000
```
接著取 `1023` 的餘數。這個可以用前面的[計算位數和](https://hackmd.io/@d4cWS3kPSNiNdRbaqbKmqQ/MUDAMUDAMUDA#計算位數和)中的結論。因為`b * 0x0202020202ULL & 0x010884422010ULL` 有下面這個形式：：
```c
tmp = (00000D0000 << 0) +
(00A000E000 << 10) +
(000B000F00 << 20) +
(0000C000G0 << 30) +
(000000000H << 40)
```
所以，他可以表示成 $2^{10}$ 為基數的形式。把以 $2^{10}$ 為基數的位數和加總起來，就會發現結果剛好是反轉之後的位元：
```c
res = 00000D0000 +
00A000E000 +
000B000F00 +
0000C000G0 +
000000000H
---------------------
00ABCFEFGH
```
### 方法四：用 64 位元乘法反轉一個 byte
這個方法跟上面很像：
```c=
b = ((b * 0x80200802ULL) & 0x0884422110ULL) * 0x0101010101ULL >> 32;
```
第一眼看過去會發現： `& 0x0884422110ULL` 這步跟前面很像，只是重複 `b` 的週期不一樣，以及最後是用「乘法」而非取「餘數」。
具體的機制可以把它前幾部展開看看：
```c
HGFEDCBA (b)
00000000 10000000 00100000 00001000 00000010 (0x0080200802)
------------------------------------------------------------------------- (*)
0HGFEDCB A00HGFED CBA00HGF EDCBA00H GFEDCBA0
00001000 10000100 01000010 00100001 00010000 (0x0884422110)
------------------------------------------------------------------------- (&)
0000E000 A0000F00 0B0000G0 00C0000H 000D0000 (=:tmp)
```
如果稍微觀察一下，可以發現如果把它視為 256 進位，並且對每一個「位數」(也就是組成這個數的每個 byte 都拿出來)做加總的話，其實就是位元組逆轉後的結果：
```
000D0000
00C0000H
0B0000G0
A0000F00
0000E000
---------(+)
ABCDEFGH
```
這件事最前面也提到：壓縮位元可以用乘法來做。把這個過程展開如下：觀察：
```clike
res = ((tmp << 0) +
(tmp << 8) +
(tmp << 16) +
(tmp << 24) +
(tmp << 32)) >> 32
```
也就是：
```clike
tmp *= 0x010101010101ULL
```
這個過程用直式展開：
```
0000E000 A0000F00 0B0000G0 00C0000H 000D0000
0000E000 A0000F00 0B0000G0 00C0000H 000D0000
0000E000 A0000F00 0B0000G0 00C0000H 000D0000
0000E000 A0000F00 0B0000G0 00C0000H 000D0000
0000E000 A0000F00 0B0000G0 00C0000H 000D0000
```
注意最中間的一行，出現了我們需要的「加總」：
```
000D0000
00C0000H
0B0000G0
A0000F00
0000E000
```
而且可以觀察到：整個過程都不會產生進位。所以最後只要把這個結果取出來就好。而可以觀察到：要取出來就是做 `>> 32`。
本來看起來應該還要用 `& 0xFF` 去取出最低的 8 位元，但因為「無號數超出範圍會自動取模數」的行為，所以會自動留下最後 8 位元的結果。
### 方法五：不用 64 位元乘法反轉 1 個 byte
```c=
b = ((b * 0x0802LU & 0x22110LU) | (b * 0x8020LU & 0x88440LU)) * 0x10101LU >> 16;
```
跟使用 64 位元的版本比較的話：
```c
b = ((b * 0x80200802ULL) & 0x0884422110ULL) * 0x0101010101ULL >> 32;
```
會發現這個做法跟上面大致雷同，連 magic number 也一樣，只是拆成兩半做而已。原本使用 64 位元乘法的過程是這樣：
```c
HGFEDCBA (b)
00000000 10000000 00100000 00001000 00000010 (0x0080200802)
------------------------------------------------------------------------- (*)
0HGFEDCB A00HGFED CBA00HGF EDCBA00H GFEDCBA0
00001000 10000100 01000010 00100001 00010000 (0x0884422110)
------------------------------------------------------------------------- (&)
0000E000 A0000F00 0B0000G0 00C0000H 000D0000 (=:tmp)
```
如果把上面不用 64 位元乘法程式稍微改寫成這樣：
```c
l = b * 0x0802LU & 0x22110LU;
r = b * 0x8020LU & 0x88440LU;
b = (l | r) * 0x10101U >> 16;
```
那求出 `l` 與 `r` 的過程，其實就是 使用 64 位元乘法的版本切成一半而已。像下面這樣：
```c
HGFEDCBA (b)
00000000 10000000 00100000 00001000 00000010 (0x8020 0802)
-------------------------- ----------------------------------------------- (*)
0HGFEDCB A00HGFED CBA00HGF EDCBA00H GFEDCBA0
00001000 10000100 01000010 00100001 00010000 (0x08844 22110)
-------------------------- ----------------------------------------------- (&)
0000E000 A0000F00 0B0000G0(l) 00C0000H 000D0000 (r)
```
而 `(l | r)` 就會變成：
```c
0000E000 A0000F00 0B0000G0 (l)
00C0000H 000D0000 (r)
----------------------------------(OR)
0000E000 A0C00F0H 0B0D00G0 (l|r)
```
接下來做直式乘法的時候，就會變成：
```
0000E000 A0C00F0H 0B0D00G0
0000E000 A0C00F0H 0B0D00G0
0000E000 A0C00F0H 0B0D00G0
----------------------------------------------(+)
ABCDEFGH
```
## 取餘數
### 正數除以 (1 << s) 的餘數
實際上只有一行：
```c=
m = n & ((1U << s) - 1);
```
這個做法很好懂：取餘數就是取最低的 `s` 位元，而 `(1U << s) - 1` 剛好就會生出取最低位元的那個 `mask`。跟計算機組織算 direct map 是同一個算法。
不過為了避免 integer promotion 各種麻煩還是附上完整的：
```c=
const unsigned int n; // numerator
const unsigned int s;
const unsigned int d = 1U << s;
unsigned int m; // m will be n % d
m = n & ((1U << s) - 1);
```
### 正數除以 (1 << s) - 1 的餘數 (迴圈)
這個技巧就跟「十進位算除以 9 的餘數」「十六進位算除以 15 的餘數」的做法是一樣的：一直加總在那個基數下表示的各個位數，直到小於除數：
```c=
unsigned int n; // numerator
const unsigned int s; // s > 0
const unsigned int d = (1 << s) - 1;
unsigned int m; // n % d goes here.
for (m = n; n > d; n = m)
for (m = 0; n; n >>= s)
m += n & d;
// Now m is a value from 0 to d, but since with modulus division
// we want m to be 0 when it is d.
m = m == d ? 0 : m;
```
### 正數除以 (1 << s) - 1 的餘數 (查表)
既然是計算位數的總和，能不能仿照 population count 那邊的方法，用分治法來計算？答案是可以的。只是 magic number 要先建好。這邊的第一步跟 population count 的第一步是類似的，只是依照 `s` 的大小，調整「柵欄」的寬度：
```c=
unsigned int n; // numerator
const unsigned int s; // s > 0
const unsigned int d = (1 << s) - 1;
unsigned int m; // n % d goes here.
static const unsigned int M[] =
{
0x00000000, 0x55555555, 0x33333333, 0xc71c71c7,
0x0f0f0f0f, 0xc1f07c1f, 0x3f03f03f, 0xf01fc07f,
0x00ff00ff, 0x07fc01ff, 0x3ff003ff, 0xffc007ff,
0xff000fff, 0xfc001fff, 0xf0003fff, 0xc0007fff,
0x0000ffff, 0x0001ffff, 0x0003ffff, 0x0007ffff,
0x000fffff, 0x001fffff, 0x003fffff, 0x007fffff,
0x00ffffff, 0x01ffffff, 0x03ffffff, 0x07ffffff,
0x0fffffff, 0x1fffffff, 0x3fffffff, 0x7fffffff
};
m = (n & M[s]) + ((n >> s) & M[s]);
```
用 `s = 7` 舉例。為了方便，用英文字母幫這個整數的有效位數做劃分。劃分成 `EEEE DDDDDDD CCCCCCC BBBBBBB AAAAAAA`。
1. 首先，執行 `m = (n & M[s]) + ((n >> s) & M[s]`
```c
EEEE DDDDDDD CCCCCCC BBBBBBB AAAAAAA (n)
EEEE 0000000 CCCCCCC 0000000 AAAAAAA (n & M[s])
0000 DDDDDDD 0000000 BBBBBBB ((n & s) >> & M[s])
--------------------------------------------------------------------------(+)
EEEE 000000F FFFFFFF 000000G GGGGGGG (m = (n & M[s]) + ((n >> s) & M[s]))
```
2. 接著要加總 `F FFFFFFFF` 及 `G GGGGGGGG`。`(m >> 14) + (m & 0x3fff)`。像這樣：
```c
EEEE 000000F FFFFFFF (m >> 14)
000000G GGGGGGG (m & 0x3fff)
--------------------------------------------------------------------------(+)
EEEE 00000HH HHHHHHH (m = (m >> 14) + (m & 0x3fff))
```
3. 然後進行一個觀察：如果要繼續這樣加總，直到有效位數的寬度比 7 小，那麼只要做 `m = (m >> 7) + (m & 0x7f)` 3 次就好。
4. 可以把 2 ~ 3 用到的 mask 跟位移量都記起來：
```c
int Q[] = {14, 7, 7, 7};
int R[] = {0x00003fff, 0x0000007f, 0x0000007f, 0x0000007f};
```
然後用一個迴圈來做：
```c
for (const unsigned int * q = &Q[0], * r = &R[0]; m > d; q++, r++) {
m = (m >> *q) + (m & *r);
}
```
接著就問：那可不可以從 `s = 0` 一直到 `s = 32` 的所有狀況，都做上面的討論，列舉出每一步需要的位移量，以及需要的 `mask`。然後就可以如法炮製？這個就是裡面程式碼的做法：
```c=
static const unsigned int M[] = {
0x00000000, 0x55555555, 0x33333333, 0xc71c71c7,
0x0f0f0f0f, 0xc1f07c1f, 0x3f03f03f, 0xf01fc07f,
0x00ff00ff, 0x07fc01ff, 0x3ff003ff, 0xffc007ff,
0xff000fff, 0xfc001fff, 0xf0003fff, 0xc0007fff,
0x0000ffff, 0x0001ffff, 0x0003ffff, 0x0007ffff,
0x000fffff, 0x001fffff, 0x003fffff, 0x007fffff,
0x00ffffff, 0x01ffffff, 0x03ffffff, 0x07ffffff,
0x0fffffff, 0x1fffffff, 0x3fffffff, 0x7fffffff
};
static const unsigned int Q[][6] = {
{ 0, 0, 0, 0, 0, 0}, {16, 8, 4, 2, 1, 1}, {16, 8, 4, 2, 2, 2},
{15, 6, 3, 3, 3, 3}, {16, 8, 4, 4, 4, 4}, {15, 5, 5, 5, 5, 5},
{12, 6, 6, 6 , 6, 6}, {14, 7, 7, 7, 7, 7}, {16, 8, 8, 8, 8, 8},
{ 9, 9, 9, 9, 9, 9}, {10, 10, 10, 10, 10, 10}, {11, 11, 11, 11, 11, 11},
{12, 12, 12, 12, 12, 12}, {13, 13, 13, 13, 13, 13}, {14, 14, 14, 14, 14, 14},
{15, 15, 15, 15, 15, 15}, {16, 16, 16, 16, 16, 16}, {17, 17, 17, 17, 17, 17},
{18, 18, 18, 18, 18, 18}, {19, 19, 19, 19, 19, 19}, {20, 20, 20, 20, 20, 20},
{21, 21, 21, 21, 21, 21}, {22, 22, 22, 22, 22, 22}, {23, 23, 23, 23, 23, 23},
{24, 24, 24, 24, 24, 24}, {25, 25, 25, 25, 25, 25}, {26, 26, 26, 26, 26, 26},
{27, 27, 27, 27, 27, 27}, {28, 28, 28, 28, 28, 28}, {29, 29, 29, 29, 29, 29},
{30, 30, 30, 30, 30, 30}, {31, 31, 31, 31, 31, 31}
};
static const unsigned int R[][6] = {
{0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000},
{0x0000ffff, 0x000000ff, 0x0000000f, 0x00000003, 0x00000001, 0x00000001},
{0x0000ffff, 0x000000ff, 0x0000000f, 0x00000003, 0x00000003, 0x00000003},
{0x00007fff, 0x0000003f, 0x00000007, 0x00000007, 0x00000007, 0x00000007},
{0x0000ffff, 0x000000ff, 0x0000000f, 0x0000000f, 0x0000000f, 0x0000000f},
{0x00007fff, 0x0000001f, 0x0000001f, 0x0000001f, 0x0000001f, 0x0000001f},
{0x00000fff, 0x0000003f, 0x0000003f, 0x0000003f, 0x0000003f, 0x0000003f},
{0x00003fff, 0x0000007f, 0x0000007f, 0x0000007f, 0x0000007f, 0x0000007f},
{0x0000ffff, 0x000000ff, 0x000000ff, 0x000000ff, 0x000000ff, 0x000000ff},
{0x000001ff, 0x000001ff, 0x000001ff, 0x000001ff, 0x000001ff, 0x000001ff},
{0x000003ff, 0x000003ff, 0x000003ff, 0x000003ff, 0x000003ff, 0x000003ff},
{0x000007ff, 0x000007ff, 0x000007ff, 0x000007ff, 0x000007ff, 0x000007ff},
{0x00000fff, 0x00000fff, 0x00000fff, 0x00000fff, 0x00000fff, 0x00000fff},
{0x00001fff, 0x00001fff, 0x00001fff, 0x00001fff, 0x00001fff, 0x00001fff},
{0x00003fff, 0x00003fff, 0x00003fff, 0x00003fff, 0x00003fff, 0x00003fff},
{0x00007fff, 0x00007fff, 0x00007fff, 0x00007fff, 0x00007fff, 0x00007fff},
{0x0000ffff, 0x0000ffff, 0x0000ffff, 0x0000ffff, 0x0000ffff, 0x0000ffff},
{0x0001ffff, 0x0001ffff, 0x0001ffff, 0x0001ffff, 0x0001ffff, 0x0001ffff},
{0x0003ffff, 0x0003ffff, 0x0003ffff, 0x0003ffff, 0x0003ffff, 0x0003ffff},
{0x0007ffff, 0x0007ffff, 0x0007ffff, 0x0007ffff, 0x0007ffff, 0x0007ffff},
{0x000fffff, 0x000fffff, 0x000fffff, 0x000fffff, 0x000fffff, 0x000fffff},
{0x001fffff, 0x001fffff, 0x001fffff, 0x001fffff, 0x001fffff, 0x001fffff},
{0x003fffff, 0x003fffff, 0x003fffff, 0x003fffff, 0x003fffff, 0x003fffff},
{0x007fffff, 0x007fffff, 0x007fffff, 0x007fffff, 0x007fffff, 0x007fffff},
{0x00ffffff, 0x00ffffff, 0x00ffffff, 0x00ffffff, 0x00ffffff, 0x00ffffff},
{0x01ffffff, 0x01ffffff, 0x01ffffff, 0x01ffffff, 0x01ffffff, 0x01ffffff},
{0x03ffffff, 0x03ffffff, 0x03ffffff, 0x03ffffff, 0x03ffffff, 0x03ffffff},
{0x07ffffff, 0x07ffffff, 0x07ffffff, 0x07ffffff, 0x07ffffff, 0x07ffffff},
{0x0fffffff, 0x0fffffff, 0x0fffffff, 0x0fffffff, 0x0fffffff, 0x0fffffff},
{0x1fffffff, 0x1fffffff, 0x1fffffff, 0x1fffffff, 0x1fffffff, 0x1fffffff},
{0x3fffffff, 0x3fffffff, 0x3fffffff, 0x3fffffff, 0x3fffffff, 0x3fffffff},
{0x7fffffff, 0x7fffffff, 0x7fffffff, 0x7fffffff, 0x7fffffff, 0x7fffffff}
};
unsigned int n; // numerator
const unsigned int s; // s > 0
const unsigned int d = (1 << s) - 1;
unsigned int m; // n % d goes here.
m = (n & M[s]) + ((n >> s) & M[s]);
for (const unsigned int * q = &Q[s][0], * r = &R[s][0]; m > d; q++, r++) {
m = (m >> *q) + (m & *r);
}
m = m == d ? 0 : m;
```
這個表看起來很恐怖，不過大致上就是依照 `s` 的數值，枚舉上面 2. ~ 3. 步所需要的位移量 (`M`) 跟 mask (`R`)。所以可以發現： `M[i][j]` 的數字是多少，`R[i][j]` 對應的 mask 就是「最低的 `M[i][j]` 位元是 `1`，而這以上的位元都是 `0`」的 mask。
Sign in with Wallet
Connect another wallet