Bit Twiddling Hacks 解析 (二)

下一篇：Bits Twiddling Hacks 解析 (三)
上一篇：Bits Twiddling Hacks 解析 (一)

• Bit Twiddling Hacks 解析 (二)
  • 不用位元計算的技巧
    • 用乘法重複一組數字
    • 用乘法壓縮位元
    • 用取餘數計算位數和
    • 小提示
  • Population Count
    • 方法一：清除所有位元
    • 方法二：暴力查表
    • 方法三 (之一)：Divide and Conquer
    • 方法三 (之二)：Divide and Conquer + 乘法
    • 方法四（之一）：加總位數和 (14 位元數值 + 64 位元運算)
    • 方法四（之二）：加總位數和 (24 位元數值 + 64 位元運算)
    • 方法四（之三）：加總位數和 (32 位元數值 + 64 位元運算)
  • Parity Bit
    • 方法一：直覺的方法
    • 方法二：暴力查表
    • 方法三：加總位數和 (8-bit 限定)
    • 方法四：直式乘法
    • 方法五：全部 XOR 起來
  • 反轉所有位元
    • 方法一：暴力查表
    • 方法二：分治法
    • 方法三：「64 位元乘法 + 取餘數」反轉一個 byte
    • 方法四：用 64 位元乘法反轉一個 byte
    • 方法五：不用 64 位元乘法反轉 1 個 byte
  • 取餘數
    • 正數除以 (1 << s) 的餘數
    • 正數除以 (1 << s) - 1 的餘數 (迴圈)
    • 正數除以 (1 << s) - 1 的餘數 (查表)

不用位元計算的技巧

這種技巧主要是建立在算術之上，類似的比如國小計算除以 9 的餘數，只要把每個位數加起來再除以 9 就好，不用真的做一次除法。

下面這邊有一些技巧是建立在類似的觀念上，僅利用加減乘除，沒有用位元運算。在十進位中大部分都有對應的做法，所以可以借助已經熟悉的十進位進行理解。

這種做法可以讓運算子看起來很少，而且有種到處都是 magic number 的感覺，比如 3 個運算計算 14 位元整數的 population count 其中一個解法就有用到。

用乘法重複一組數字

先看十進位的例子：

$$1,001,001,001,001,001\times 987=987,987,987,987,987,987$$

主要就是某個數字，乘上重複的

$100...001$ ，就會把本來的數字重複很多次。不管是寫成科學記號或是用直式乘法來證明都很容易看出來這件事。

接下來看看一個 16 進位的版本：

$$\mathtt{res}\mathtt{=}\mathtt{0}\mathtt{x}\mathtt{200040008001}\mathtt{ULL}\times \mathtt{0}\mathtt{xBAD}$$

這個結果寫成 16 進位的話，會是 0x175a2eb45d68bad ，乍看之下不是很好懂。但如果用直式乘法去思考的話，其實會是：

res = (0xBAD) + 
      (0xBAD << 15) + 
      (0xBAD << 30) + 
      (0xBAD << 45);

原理也是一樣。如果把上面的東西寫成二進位的話，其實是在做：

mask = 0b1000000000000001000000000000001000000000000001;
bad  = 0b101110101101;    // 0xBAD
res  = mask * bad;

因為 64 位元的關係，把完整的數字寫出來實在是很長，所以大部分都會寫 16 進位，然後就會變成一個有點令人費解的樣子。

用乘法壓縮位元

這個用直式乘法非常好證明。舉例來說，現在有一個 64 位元無號數的 2 進位 形式長這樣：

tmp = 0000000H 000000G0 00000F00 0000E000 000D0000 00C00000 0B000000 A0000000

其中 A、B、C、D、E、F、G、H 都是 0 或 1。如果做：

tmp = (tmp << 0)  + 
      (tmp << 8)  + 
      (tmp << 16) + 
      (tmp << 24) + 
      (tmp << 32) + 
      (tmp << 40) + 
      (tmp << 48) + 
      (tmp << 56);

那可以發現計算結果最高位的 8 個位元，會是 (觀察對角線方向)：

0000000H 
000000G0 
00000F00 
0000E000 
000D0000 
00C00000 
0B000000 
A0000000
---------(+)
ABCDEFGH

回到剛剛的計算：那一連串的位移其實可以簡化成：

tmp = tmp * 0x0101010101010101ULL;

所以總結以上，做：

res = (tmp * 0x0101010101010101ULL) >> 56; 

就會有把稀疏的位元壓縮再一起的效果。

用取餘數計算位數和

有一些問題的思路可以從「計算位數和」下手。比如說 population count 就是在算每個位數的和; 而 parity bit 就是算完位數和之後，再看奇偶性 (比如做 x & 1)。

暖身：以 10 為基數的位數和

這件事情的原理，跟國小時判斷一個數是不是 9 的倍數，使用「(位數和) % 9，看看餘數是不是 0」，的原理是很像的。比如說想要計算：

$$9487÷9$$

的餘數，那麼小時候大家都學過，去計算：

$$(9+4+8+7)÷9$$

的餘數就好。這是因為對於任意非負整數

$k$，有：

$${10}^k=\underset{k}{\underbrace{999...9}}+1$$

因此，假定一個數的 10 進位表示法是

$a_k\dots a_2{a}_1{a}_0$，或是說：

$$n=\sum _{i=0}^k{a}_i\times {10}^i$$

那麼可以把所有的

${10}^i$ 拆成

$(\underset{i}{\underbrace{999...9}}+1)$，然後經過一點計算之後：

$$\begin{array}{rlr}n& =\sum _{i=0}^k{a}_i\times (\underset{i}{\underbrace{999\dots 9+1}})& \mathrm{mod}9\\ & =\sum _{i=0}^k{a}_i\times (\underset{i}{\underbrace{999\dots 9}})+(a_i\times 1)& \mathrm{mod}9\\ & =\sum _{i=0}^k{a}_i& \mathrm{mod}9\end{array}$$

但如果位數的總和比 9 小，那麼有 mod 跟沒有 mod 結果是一樣的，那麼除法計算的商根本就是 0。也就是說

$\text{mod}$ 根本可以拿掉，像這樣：

$$\begin{array}{rl}n& =\sum _{i=0}^k{a}_i\mathrm{mod}9\\ & =\sum _{i=0}^k{a}_i(\text{if }\sum _{i=0}^k{a}_i<9)\end{array}$$

比如說，想計算

$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{0}\mathtt{x}\underset{\mathtt{i}}{\underbrace{\mathtt{FF}\mathtt{.}\mathtt{.}\mathtt{.}\mathtt{FF}}}+1$

$$\begin{array}{rlr}n& =\sum _{i=0}^7{b}_i\times {16}^i& \mathrm{mod}15\\ & =\sum _{i=0}^7{b}_i\times (\mathtt{0}\mathtt{x}\underset{\mathtt{i}}{\underbrace{\mathtt{FF}\mathtt{.}\mathtt{.}\mathtt{.}\mathtt{FF}}}+1)& \mathrm{mod}\mathtt{0}\mathtt{xF}\\ & =\sum _{i=0}^7{b}_i& \mathrm{mod}\mathtt{0}\mathtt{xF}\\ & =\sum _{i=0}^7{b}_i& \text{if }\sum _{i=0}^7{b}_i<\mathtt{0}\mathtt{xF}\end{array}$$

所以可以一眼知道位數和就是每個位元加起來，就是

$7$。

小提示

後面有很多利用乘法的招式，大致上都是遵照下面的形式：

1. 在 64 位元中，用乘法重複一個 byte 的內容很多次。
2. 依照特定間隔取出位元，使得該 byte 中每個位元恰好在 mask 後的結果出現一次。
3. 用取餘數計算所有位元的加總(比如 population count 中的 方法四（之一）：加總位數和 (14 位元數值 + 64 位元運算); 或用上面的乘法把所有位元壓縮在一起(比如反轉所有位元中的方法四：用乘法反轉一個-byte)

舉例來說下面三者：

/* 對一個 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。最天真的方法就是用迴圈：

/* @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

for (res = 0; n; n &= n-1) {
    res++;
}

這個是用上面「清除最低位元」的做法去做。在位元很稀疏的時候，有機會比天真的做法快。不過還是免不了使用 branch。

方法二：暴力查表

原文：Counting bits set by lookup table

就是那個傳說中把結果存在月球上的做法啊。

這個做法是建立一個大小是 256 的陣列 table 當作表，枚舉所有 8 位元的整數所有可能數值中，對應的 1 的個數。比如說我想要查 73 這個數字中，有多少個 1，那就去找 table[73]。所以這個有點大的表會像是：

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

不過這樣打起來有點煩瑣，所以可以用巨集的技巧簡化：

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}\mathtt{(}\mathtt{n}\mathtt{(}\mathtt{[}\mathtt{N}\mathtt{:}\mathtt{0}\mathtt{]}\mathtt{)}\{\begin{array}{ll}0+\mathtt{pop}\mathtt{(}\mathtt{n}\mathtt{[}\mathtt{N}\mathtt{-}\mathtt{2}\mathtt{:}\mathtt{0}\mathtt{]}\mathtt{)}& \text{if }\mathtt{n}\mathtt{[}\mathtt{N}\mathtt{:}\mathtt{N}\mathtt{-}\mathtt{1}\mathtt{]}\mathtt{=}\mathtt{0}\\ 1+\mathtt{pop}\mathtt{(}\mathtt{n}\mathtt{[}\mathtt{N}\mathtt{-}\mathtt{2}\mathtt{:}\mathtt{0}\mathtt{]}\mathtt{)}& \text{if }\mathtt{n}\mathtt{[}\mathtt{N}\mathtt{:}\mathtt{N}\mathtt{-}\mathtt{1}\mathtt{]}\mathtt{=}\mathtt{1}\\ 1+\mathtt{pop}\mathtt{(}\mathtt{n}\mathtt{[}\mathtt{N}\mathtt{-}\mathtt{2}\mathtt{:}\mathtt{0}\mathtt{]}\mathtt{)}& \text{if }\mathtt{n}\mathtt{[}\mathtt{N}\mathtt{:}\mathtt{N}\mathtt{-}\mathtt{1}\mathtt{]}\mathtt{=}\mathtt{2}\\ 2+\mathtt{pop}\mathtt{(}\mathtt{n}\mathtt{[}\mathtt{N}\mathtt{-}\mathtt{2}\mathtt{:}\mathtt{0}\mathtt{]}\mathtt{)}& \text{if }\mathtt{n}\mathtt{[}\mathtt{N}\mathtt{:}\mathtt{N}\mathtt{-}\mathtt{1}\mathtt{]}\mathtt{=}\mathtt{3}\end{array}$$

原理就是照最高 2 個位元進行討論。他有可能是 0b00, 0b01, 0b10, 0b11，而 pop count 的結果就是「最高 2 位元的 1 的數目」加上「其餘低位元的 1 的數目」。

對於其他 8 的倍數寬度的數值型別，只要多做幾次就可以了：

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:

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

這個應該是最有名的做法吧？主要的思路是這樣：

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

第一行做法為什麼會對，可以窮舉所有可能的狀況：

0b11 - (0b11 >> 1) = 2
0b10 - (0b10 >> 1) = 1
0b01 - (0b01 >> 1) = 1
0b00 - (0b00 >> 1) = 0

所以就剛好跟原先的第一步是一樣的。接下來的做法大上就是跟原先一樣。

方法三 (之二)：Divide and Conquer + 乘法

原文：Counting bits set, in parallel 的下面

概念大致上跟上面一樣，只是後面用乘法做：

v = v - ((v >> 1) & 0x55555555);
v = (v & 0x33333333) + ((v >> 2) & 0x33333333);
r = ((v + (v >> 4) & 0xF0F0F0F) * 0x1010101) >> 24;

最後面乘法的作用是像下面這樣：

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

這邊使用到另外一個概念：population count = 加總所有位數。所以一個可能的做法是下面這樣：

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. 把取出來的位元加總。

不過在文章中出現的程式碼很短，是這樣：

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 之後會變成：

tmp1 = 0xBAD + 
       0xBAD << 15 + 
       0xBAD << 30 + 
       0xBAD << 45;

反正就是把 0xBAD 重複寫很多次。但關鍵就是：他會以每 15 個位元為週期重複。

2. & 0x111111111111111ULL

如果把每個位元用 0DCBA9876543210 表示每個位元，(比如第 A 代表第 10 個位元、4 代表第 4 位元等等)的話，上面那個過程其實是「讓 n 中每一個位元恰好在結果中出現 1 次」。如果用 0DCBA9876543210 分別代表 0 ~ 15 個二進位中的位數的話，上面的步驟其實是照下面的過程取出對應的位數：

$$\begin{array}{r}\mathtt{0}\mathtt{DC}\underset{\uparrow }{\mathtt{B}}\mathtt{A}\mathtt{98}\underset{\uparrow }{\mathtt{7}}\mathtt{654}\underset{\uparrow }{\mathtt{3}}\mathtt{210}\underset{\uparrow }{\mathtt{0}}\mathtt{DCB}\underset{\uparrow }{\mathtt{A}}\mathtt{987}\underset{\uparrow }{\mathtt{6}}\mathtt{543}\underset{\uparrow }{\mathtt{2}}\mathtt{10}\mathtt{0}\underset{\uparrow }{\mathtt{D}}\mathtt{CBA}\underset{\uparrow }{\mathtt{9}}\mathtt{876}\underset{\uparrow }{\mathtt{5}}\mathtt{432}\underset{\uparrow }{\mathtt{1}}\mathtt{0}\mathtt{0}\mathtt{D}\underset{\uparrow }{\mathtt{C}}\mathtt{BA}\mathtt{9}\underset{\uparrow }{\mathtt{8}}\mathtt{765}\underset{\uparrow }{\mathtt{4}}\mathtt{321}\underset{\uparrow }{\mathtt{0}}\\ \mathtt{000}\underset{\uparrow }{\mathtt{1}}\mathtt{000}\underset{\uparrow }{\mathtt{1}}\mathtt{000}\underset{\uparrow }{\mathtt{1}}\mathtt{000}\underset{\uparrow }{\mathtt{0}}\mathtt{000}\underset{\uparrow }{\mathtt{1}}\mathtt{000}\underset{\uparrow }{\mathtt{1}}\mathtt{000}\underset{\uparrow }{\mathtt{1}}\mathtt{00}\mathtt{0}\underset{\uparrow }{\mathtt{1}}\mathtt{000}\underset{\uparrow }{\mathtt{1}}\mathtt{000}\underset{\uparrow }{\mathtt{1}}\mathtt{000}\underset{\uparrow }{\mathtt{1}}\mathtt{0}\mathtt{00}\underset{\uparrow }{\mathtt{1}}\mathtt{000}\underset{\uparrow }{\mathtt{1}}\mathtt{000}\underset{\uparrow }{\mathtt{1}}\mathtt{000}\underset{\uparrow }{\mathtt{1}}\\ \mathtt{000}\underset{\uparrow }{\mathtt{B}}\mathtt{000}\underset{\uparrow }{\mathtt{7}}\mathtt{000}\underset{\uparrow }{\mathtt{3}}\mathtt{000}\underset{\uparrow }{\mathtt{0}}\mathtt{000}\underset{\uparrow }{\mathtt{A}}\mathtt{000}\underset{\uparrow }{\mathtt{6}}\mathtt{000}\underset{\uparrow }{\mathtt{2}}\mathtt{00}\mathtt{0}\underset{\uparrow }{\mathtt{D}}\mathtt{000}\underset{\uparrow }{\mathtt{9}}\mathtt{000}\underset{\uparrow }{\mathtt{5}}\mathtt{000}\underset{\uparrow }{\mathtt{1}}\mathtt{0}\mathtt{00}\underset{\uparrow }{\mathtt{C}}\mathtt{000}\underset{\uparrow }{\mathtt{8}}\mathtt{000}\underset{\uparrow }{\mathtt{4}}\mathtt{000}\underset{\uparrow }{\mathtt{0}}\end{array}$$

也就是依照：

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{tmp}\mathtt{2}\mathtt{:=}\mathtt{(}\mathtt{n}\mathtt{\ast }\mathtt{0}\mathtt{x}\mathtt{200040008001}\mathtt{ULL}\mathtt{\&}\mathtt{0}\mathtt{x}\mathtt{111111111111111}\mathtt{ULL}\mathtt{)};$$

這個結果會是個位元形如下面這樣的數字：

$$\mathtt{tmp}\mathtt{2}\mathtt{=}\mathtt{0}\mathtt{b}\mathtt{0.}\mathtt{.}\mathtt{.}\underset{\mathtt{4}\times \mathtt{14}\mathtt{=}\mathtt{56}\mathtt{\text{ bits}}}{\underbrace{\mathtt{000}{\mathtt{n}}_{\mathtt{11}}\mathtt{000}{\mathtt{n}}_{\mathtt{7}}\mathtt{000}{\mathtt{n}}_{\mathtt{3}}\dots \mathtt{000}{\mathtt{n}}_{\mathtt{1}}\mathtt{000}{\mathtt{n}}_{\mathtt{12}}\mathtt{000}{\mathtt{n}}_{\mathtt{8}}\mathtt{000}{\mathtt{n}}_{\mathtt{4}}\mathtt{000}{\mathtt{n}}_{\mathtt{0}}}}$$

現在我們有什麼？

1. 所有

   $n$ 的位元都在裡面剛好出現過一次

2. 目前的計算結果是：

   $$\begin{array}{rl}\text{tmp2}=n_0& \cdot {16}^0\\ n_4& \cdot {16}^1& +\\ n_8& \cdot {16}^2& +\\ n_{12}& \cdot {16}^3& +\\ n_1& \cdot {16}^4& +\\ n_4& \cdot {16}^5& +\\ \dots \\ n_7& \cdot {16}^{12}& +\\ n_{11}& \cdot {16}^{13}\end{array}$$

   現在所有 n 的位元都恰好在 tmp2 中出現一次，其他位元都是 0。現在只欠東風：把 tmp 的每個位元加總。

3. % 0xF

因為現在基數是 16，加上位數和最大就是 14 (14 個位元都是 1) < (16 - 1)。用前面計算位數和的結論可以知道：這時做 % (16 - 1) 就可以求出位數和。

方法四（之二）：加總位數和 (24 位元數值 + 64 位元運算)

原文：Counting bits set in 14, 24, or 32-bit words using 64-bit instructions

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

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

最直覺的方法就是做：

while (x) {
    parity = !parity;
    x = x & (x - 1);
}

不過照樣子會希望避開 branch 跟 jump。做法大概就是兩類：

1. 從 popullation count 的做法改來
2. parity bit 就是想辦法計算所有位元 XOR 起來的值。所以往這個方向做也可以。

方法二：暴力查表

原文：Compute parity by lookup table

這個跟 population count 那邊的查表法類似。可以建立一個大小 256 的表，枚舉 0 到 255 中所有值的 parity：

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 的巨集建立方式。所以程式就可以簡化成：

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：

// Variation:
unsigned char * p = (unsigned char *) &v;
parity = ParityTable256[p[0] ^ p[1] ^ p[2] ^ p[3]];

或是先把狀態「壓縮」到 8 bit 裡面，再查表：

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

unsigned char b;  // byte value to compute the parity of
bool parity = 
  (((b * 0x0101010101010101ULL) & 0x8040201008040201ULL) % 0x1FF) & 1;

這個就是從 population count 的方法四（之一）：加總位數和 (14 位元數值 + 64 位元運算)改來的。先做 population count，再看結果的奇偶性。

方法四：直式乘法

原文：Compute parity of word with a multiply

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 位元其實是在做：

v = v + 
    (v << 4) + 
    (v << 8) + 
    (v << 12) + 
    (v << 16) + 
    (v << 20) + 
    (v << 24) + 
    (v << 28)

所以第 28 位元，就會是加總結果的個位數。因此把它取出來，就是 parity bit。

方法五：全部 XOR 起來

原文：Compute parity in parallel

可以直接做 32 次。不過因為 XOR 有結合律跟交換律，所以可以用分治法，一直「對折」來找出 parity bit：

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 去減少次數：

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

問題：給定一個 unsigned int v，反轉 v 的所有位元。即結果 c 的第 i 位元，必須等於 v 的第 31 - i 位元。其中 i 的範圍是 0 ~ 31。

最直覺的做法就是用迴圈下去：

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 的陣列裡面：

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

類似地，這個表也可以用巨集去建立：

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{array}{rl}\mathtt{abcdefgh}=& \mathtt{(}\mathtt{ab}\mathtt{\ast }\mathtt{64}\mathtt{)}\mathtt{+}\\ & \mathtt{(}\mathtt{cd}\mathtt{\ast }\mathtt{16}\mathtt{)}\mathtt{+}\\ & \mathtt{(}\mathtt{ef}\mathtt{\ast }\mathtt{4}\mathtt{)}\mathtt{+}\\ & \mathtt{(}\mathtt{gh}\mathtt{\ast }\mathtt{1}\mathtt{)}\end{array}$$

這個過程可以想像成成要開下面這樣的陣列：

$$\mathtt{int}\mathtt{table}\mathtt{[}\mathtt{4}\mathtt{]}\mathtt{[}\mathtt{4}\mathtt{]}\mathtt{[}\mathtt{4}\mathtt{]}\mathtt{[}\mathtt{4}\mathtt{]}\mathtt{;}$$

其中：

$$\begin{array}{rl}\mathtt{table}\mathtt{[}\mathtt{hg}\mathtt{]}\mathtt{[}\mathtt{fe}\mathtt{]}\mathtt{[}\mathtt{dc}\mathtt{]}\mathtt{[}\mathtt{ba}\mathtt{]}=& \mathtt{(}\mathtt{ab}\mathtt{\ast }\mathtt{64}\mathtt{)}\mathtt{+}\\ & \mathtt{(}\mathtt{cd}\mathtt{\ast }\mathtt{16}\mathtt{)}\mathtt{+}\\ & \mathtt{(}\mathtt{ef}\mathtt{\ast }\mathtt{4}\mathtt{)}\mathtt{+}\\ & \mathtt{(}\mathtt{gh}\mathtt{\ast }\mathtt{1}\mathtt{)}\end{array}$$

另外一方面，手動展開上面的巨集，就會發現：開一個 table[4][4][4][4] 去存取 table[hg][fe][dc][ba]，在 C 語言的 array subscription 下，這件事跟開一個大小 256 的表 table[256] ，然後直接存取 table[hgfedcba] 效果，是一樣的。所以就可以用上面的巨集去做。

會做 8 位元的表之後，就可以做 32 位元的反轉：

c = (BitReverseTable256[v & 0xff] << 24) | 
    (BitReverseTable256[(v >> 8) & 0xff] << 16) | 
    (BitReverseTable256[(v >> 16) & 0xff] << 8) |
    (BitReverseTable256[(v >> 24) & 0xff]);

或是用指標：

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 位元的狀況下，就會看起來像下面這樣：

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

另外，也有迴圈版的：

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 類似的技巧：

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 的餘數。這個可以用前面的計算位數和中的結論。因為b * 0x0202020202ULL & 0x010884422010ULL 有下面這個形式：：

tmp = (00000D0000 << 0)  + 
      (00A000E000 << 10) +
      (000B000F00 << 20) +
      (0000C000G0 << 30) +
      (000000000H << 40)

所以，他可以表示成

$2^{10}$ 為基數的形式。把以

$2^{10}$ 為基數的位數和加總起來，就會發現結果剛好是反轉之後的位元：

res = 00000D0000 + 
      00A000E000 +
      000B000F00 +
      0000C000G0 +
      000000000H
---------------------
      00ABCFEFGH

方法四：用 64 位元乘法反轉一個 byte

這個方法跟上面很像：

b = ((b * 0x80200802ULL) & 0x0884422110ULL) * 0x0101010101ULL >> 32;

第一眼看過去會發現： & 0x0884422110ULL 這步跟前面很像，只是重複 b 的週期不一樣，以及最後是用「乘法」而非取「餘數」。

具體的機制可以把它前幾部展開看看：

                                    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

這件事最前面也提到：壓縮位元可以用乘法來做。把這個過程展開如下：觀察：

res = ((tmp << 0)  + 
       (tmp << 8)  + 
       (tmp << 16) + 
       (tmp << 24) + 
       (tmp << 32)) >> 32

也就是：

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

b = ((b * 0x0802LU & 0x22110LU) | (b * 0x8020LU & 0x88440LU)) * 0x10101LU >> 16; 

跟使用 64 位元的版本比較的話：

b = ((b * 0x80200802ULL) & 0x0884422110ULL) * 0x0101010101ULL >> 32;

會發現這個做法跟上面大致雷同，連 magic number 也一樣，只是拆成兩半做而已。原本使用 64 位元乘法的過程是這樣：

                                    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 位元乘法程式稍微改寫成這樣：

l = b * 0x0802LU & 0x22110LU;
r = b * 0x8020LU & 0x88440LU;
b = (l | r) * 0x10101U >> 16;

那求出 l 與 r 的過程，其實就是 使用 64 位元乘法的版本切成一半而已。像下面這樣：

                                    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) 就會變成：

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) 的餘數

實際上只有一行：

m = n & ((1U << s) - 1); 

這個做法很好懂：取餘數就是取最低的 s 位元，而 (1U << s) - 1 剛好就會生出取最低位元的那個 mask。跟計算機組織算 direct map 是同一個算法。

不過為了避免 integer promotion 各種麻煩還是附上完整的：

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 的餘數」的做法是一樣的：一直加總在那個基數下表示的各個位數，直到小於除數：

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 的大小，調整「柵欄」的寬度：

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]

   ​​​​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)。像這樣：

   ​​​​                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 跟位移量都記起來：

   ​​​​int Q[] = {14, 7, 7, 7};
   ​​​​int R[] = {0x00003fff, 0x0000007f, 0x0000007f, 0x0000007f};
   

   然後用一個迴圈來做：

   ​​​​for (const unsigned int * q = &Q[0], * r = &R[0]; m > d; q++, r++) {
   ​​​​    m = (m >> *q) + (m & *r);
   ​​​​}
   

接著就問：那可不可以從 s = 0 一直到 s = 32 的所有狀況，都做上面的討論，列舉出每一步需要的位移量，以及需要的 mask。然後就可以如法炮製？這個就是裡面程式碼的做法：

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。