---
tags: 進階電腦系統理論與實作
---
# 2020q3 Homework3 (quiz3)
contributed by < `RinHizakura` >
>[2020q3 第 3 週測驗題](https://hackmd.io/@sysprog/2020-quiz3?fbclid=IwAR2VKM07QHTneLdK1VBa9XgZ9ca10sZsEKUgVXTXumv0PBhKUGibXLKxEe0)
> [GitHub](https://github.com/RinHizakura/NCKU-System-Software-Quiz/tree/master/quiz3)
## 測驗 `1`
根據 C 語言規格書,對**有號型別的負數**算術右移(Arithmetic right shift) 在 C 語言中是 [unspecified behavior](https://en.wikipedia.org/wiki/Unspecified_behavior)。
> The result of E1 >> E2 is E1 right-shifted E2 bit positions.
> ...
> If E1 has a signed type and a negative value, the resulting value is implementation-defined.
意即其行為結果雖在語言規範中並無規定,但可以根據平台或者編譯器的參考文件得知。而算術右移無非就是 2 進位的最左位元補 0 或者補 1 兩種。由題意描述可知,期待的結果是如果原本數字為負數的情形下,右移後最左邊補 1。
因此,如果 `(m >> n)` 的結果已經如同預期(左邊補 1),則 `(fix ^ (fix >> n));` 可以為 0,也就是 `fix` 可以為 0。由此可以反推 `logical` 是在判斷右移的結果為何,因此 `OP1` = `(b) >>1` :
* 如果 -1 算術右移最左邊補 0,(-1) >> 1 是正數,> 0 為 true,則 logical = 1
* 如果 -1 算術右移最左邊補 1,(-1) >> 1 是負數,> 0 為 false,則 logical = 0。此情況下 `fix` 必為 0,符合我們預想的 `m >> n` 的結果就如預期,不需額外更動
為求圖解精簡,圖例以 8 bits 為例,推廣到 32 bits 則同理。
即使負數右移後左邊是補 0,但在 `m` 為正數的情況下,規格書明確的說 `m >> n` 的結果會等同是最左位元補 0。
> The result of E1 >> E2 is E1 right-shifted E2 bit positions. If E1 has an unsigned type or if E1 has a signed type and a nonnegative value, the value of the result is the integral part of the quotient of E1 / $2^{E2}$.
所以 `fix` 仍可為 0,由此可以反推 `OP2` = `(c) m < 0`。只有在滿足
1. 負數的右移結果為最左位元補 0
2. m 為負數
的情況下,`fixu` 會有真正的內容,且此時 `fixu` 會是 0x11111111,轉到 int 表示的 `fix` 就是 -1,則 `| (fix ^ (fix >> n))` 就會把右移的位元都補成 1,讓結果符合期待。
:::warning
問題:
```cpp
unsigned int fixu = -(logical & (m < 0));
int fix = *(int *) &fixu;
```
和
```cpp
int fix = -(logical & (m < 0));
```
兩個寫法的差異性?
> 前者可避免編譯器進行過度 (aggressive) 最佳化
> :notes: jserv
:::
### 驗證正確性
因為在我的電腦上對負數算術右移是左邊補 1,為了確認程式在左邊補 0 的狀況下同樣正確,所以我就先實作了一個很在我的電腦中可以保證左邊補 0 的 right shift,再透過 `asr_i` 去比較和左邊補 1 的 right shift 結果是否有差異,檢驗 `asr_i` 是否是可以適用於跨平台的實作。
:::warning
:warning: 顯然這個方法很暴力,而且也欠缺對 data model 的考量,不過因為重點只是確保 `asr_i` 的正確性所以就無傷大雅。
:::
```c=
int zero_rs(int a, int n)
{
for (int i = 0; i < n; i++) {
a = (a >> 1) & 0x7FFFFFFF;
}
return a;
}
int asr_i(signed int m, unsigned int n)
{
const int logical = zero_rs(((int) -1), 1) > 0;
unsigned int fixu = -(logical & (m < 0));
int fix = *(int *) &fixu;
return zero_rs(m, n) | (fix ^ zero_rs(fix, n));
}
```
### 通用的 ASR
實作可以針對不同資料寬度的 ASR。
```c=
#define asr_i(m, n) \
_Generic((m), \
int8_t: asr_i8, \
int16_t: asr_i16, \
int32_t: asr_i32, \
int64_t: asr_i64 \
)(m,n)
#define asr(type) \
const type logical = (((type) -1) >> 1) > 0; \
u##type fixu = -(logical & (m < 0)); \
type fix = *(type *) &fixu; \
return (m >> n) | (fix ^ (fix >> n))
int8_t asr_i8(int8_t m, unsigned int n)
{
asr(int8_t);
}
int16_t asr_i16(int16_t m, unsigned int n)
{
asr(int16_t);
}
int32_t asr_i32(int32_t m, unsigned int n)
{
asr(int32_t);
}
int64_t asr_i64(int64_t m, unsigned int n)
{
asr(int64_t);
}
```
根據幾點考量而如此實作:
* 根據 [data model](https://en.wikipedia.org/wiki/64-bit_computing#64-bit_data_models) 的不同,`int`、`long` 等等型別的長度可能有所差異,我認為根據 bits 實際的長度選擇正確寬度的 ASR 是比較好的作法
* 透過 `u##type`([Concatenation](https://gcc.gnu.org/onlinedocs/cpp/Concatenation.html)) 技巧連接字串
* 透過 [_Generic](https://hackmd.io/@sysprog/c-preprocessor#_Generic-C11) 達到類似 C++ template 的效果
:::warning
延伸問題:
乍想之下,無論型別都直接轉型成 64 bits 的版本處理,計算出結果後再轉回原本的型別,似乎正確性上是沒有問題的?這樣想是否完全正確呢?排除正確性的話,這樣方式的實作可能有甚麼缺點?
> 不是每款 C 語言編譯器都能正確處理 64 位元,且 data model 會涉及 ABI 議題
> :notes: jserv
:::
## 測驗 `2`
* `num > 0`: 因為 $4^n$ 必定 >= 1
* `(num & (num - 1))`: 可以判斷非負數的 `num` 是否為 $2^n$(n 為整數),因為如果 `num` 不是 $2^n$ 就必不為 $4^n$
* 在[你所不知道的 C 語言:數值系統篇](https://hackmd.io/@sysprog/c-numerics#%E9%81%8B%E7%94%A8-bit-wise-operator)有提及相關內容,或者可以參考[這個筆記](https://hackmd.io/9tGfpJtLTwyyOzDvlyJS2w?view#x-amp-x---1--0)
* 如果是 $4^n$:
| Column 1 | Column 2 |
| -------- |:--------- |
| $4^0$ | 0b1 |
| 4 | 0b100 |
| $4^2$ | 0b10000 |
| $4^3$ | 0b1000000 |
在滿足前兩個條件的前提下,從右往左數到第 1 個 1 之前的 0 總數(也就是`ctz`)為偶數個,因此 `OPQ` = `(f) & 0x1`
### 不同的實作
要判斷是否是 4 的倍數,可以有許多不同的實作方式。
```c=
bool isPowerOfFour_new(int num)
{
return !(__builtin_popcount(num) ^ 1) && __builtin_ffs(num) & 1;
}
```
一種思路是可以擅用 popcount 和 find first set,因為要判斷是否是 $4^n$ 只要滿足:
1. 2 進位表示中僅有一個 1
2. 且此 1 的位置為奇數(注意是從 1 開始數,所以 `__builtin_ffs(1) == 1`。事實上,`ctz` 計算後其實就可以達到同等 `ffs` 的效果)
或者也可以改寫成: (注意 `&&` 左右的 statement 是不能交換的,因為這裡的小技巧是通過 `&&` 左邊如果為 0 的狀況就不會去執行右邊的特性,避免 num >> 32 這樣的非法操作。)
```c=
bool isPowerOfFour_new(int num)
{
int ffs = __builtin_ffs(num);
return (ffs & 1) && !(num >> ffs);
}
```
### [1009. Complement of Base 10 Integer](https://leetcode.com/problems/complement-of-base-10-integer/)
題解為:
```c=
int bitwiseComplement(int N){
if(N == 0)
return 1;
int mask = (1 << (31-__builtin_clz(N)))- 1;
return ~N & mask;
}
```
根據題意,目標是將足以表示該數字的 2 進位表示做 complement(例如 5 可以表示成 0b101,所以有效部份為 3 個 bit)。因此,我們只需找到這個位元數 `n`,取直接對數字做 complement 的後 `n`,其他則維持為 0,就可以得到正確答案。
我們可以透過 `clz` 來計算,透過 `(1 << (31-__builtin_clz(N)))- 1` 可以得到需要的位元位置為 1,其他為 0 的 mask(舉例來說,5 = 0b000...00101,則 mask 就會是 0b000...00111),對 5 取補數後再和 `mask` 做 and 運算,就可以得到正確的結果。
附上在 leetcode 中執行正確且通過的證明:
![](https://i.imgur.com/9ijqTI8.png)
### [41. First Missing Positive](https://leetcode.com/problems/first-missing-positive/)
題解為:
```c=
#define XORSWAP_UNSAFE(a, b) \
((a) ^= (b), (b) ^= (a), \
(a) ^= (b))
int firstMissingPositive(int* nums, int numsSize){
for(int i=0;i<numsSize;i++){
if(nums[i] > 0 && nums[i] < numsSize){
int pos = nums[i] - 1;
if(pos != i && nums[i] != nums[pos]){
XORSWAP_UNSAFE(nums[i],nums[pos]);
i--;
}
}
}
for(int i = 0;i < numsSize;i++){
if(nums[i] != i + 1)
return i + 1;
}
return numsSize + 1;
}
```
此題的解題思路是:因為目標是找到最小不在陣列中的正整數解,因此我們可以經歷整個陣列,如果找到 `1` 就放在陣列的 `0` 、找到 `2`,就放在陣列的 `1`,依照此規則擺放。則再從頭檢查陣列時,如果陣列 `0` 位置不是 `1`,就表示原始的陣列中沒有 `1`,為正確答案,依照同樣的規則,就知道原始的陣列是不是有 `2`、`3`...,依此類推。如果陣列中的內容都照規則擺放,則最小的正整數解就是 `numSize + 1`。
附上在 leetcode 中執行正確且通過的證明:
![](https://i.imgur.com/ngJIqr6.png)
:::warning
想不到這題用甚麼解法會需要 `clz`
> 最後一個迴圈內的結束條件,可用到 `clz`
> :notes: jserv
:::
## Golomb coding
[Golomb coding](https://en.wikipedia.org/wiki/Golomb_coding) 是一種無損的資料壓縮方法,屬於 [Prefix code](https://en.wikipedia.org/wiki/Prefix_code) 的分組編碼系統,可以對非負整數進行編碼。其演算法流程為:
1. 根據一個定義的參數 M,對於編碼目標 N,計算出 N / M 的商 q 以及餘數 r
2. 最終的編碼會是 `<對 q 的編碼><對 r 的編碼>` 。對於 q 的編碼,使用 [unary coding](https://en.wikipedia.org/wiki/Unary_coding),對於 r 則使用 [Truncated binary encoding](https://en.wikipedia.org/wiki/Truncated_binary_encoding)
3. 因此對 q 的編碼為: q 個 1 加上 1 個 0(或者 q 個 0 加上 1 個 1)
4. 對 r 的編碼,根據情形會有不同的截斷量, 令 b = $\lfloor log_{2}M \rfloor$,則:
* 如果 r < $2^{b+1}$ - M,則用 b 個 bits 表示 r 作為編碼
* 否則,用 b + 1 個 bits 表示 r + $2^{b+1}$ - M 作為編碼
因此,每一個商 q 最多可以對應 M 種結果。詳情參照 [Wikipedia 範例](https://en.wikipedia.org/wiki/Golomb_coding#Example)可以更容易理解。
### Golomb-Rice coding
Golomb-Rice coding 是 Golomb coding 的特例,規定前面提到的參數 M 必須是 $2^N$ (N 為整數),這可以帶來幾個好處:
* 可以透過 bitwise operation 計算餘數(`r = N & (M - 1)`)
* 對於 r 的編碼,一致取 r 二進位表示的低 $log_{2}m$ 位元來表示
### Exponential-Golomb coding
[Exponential-Golomb coding](https://en.wikipedia.org/wiki/Exponential-Golomb_coding)(order-k exp-Golomb code) 是另一種對 Golomb coding 的延伸。其演算法流程為:
1. 對於一個數字 A,去掉最低的 k 個位元後加一,得到數字 B
2. 在數字 B 的前面,加上等同 B 的位元數量減 1 個零
3. 將前面去掉的 k 個位元補回尾端
與前面的方法最大的不同點是,Exponential-Golomb coding 的每一個 group 成員數量不同,會成指數型的成長。
Reference: [Golomb 及指数哥伦布编码原理介绍及实现
](https://www.cnblogs.com/wangguchangqing/p/6297792.html)
### x-compressor
[x-compressor](https://github.com/jserv/x-compressor) 是以 Golomb-Rice coding 為基礎的資料壓縮器。這裡我們先特別關注其中 Golomb coding 相關的實作。
#### `generic_io_t`
```cpp
typedef struct {
uint32_t *ptr; /* pointer to memory */
void *end;
uint32_t b; /* bit buffer */
size_t c; /* bit counter */
} generic_io_t;
```
為了可以更好的理解程式碼,大致敘述 `generic_io_t` 的用途:
* 編碼使用的位元數目不固定,但整數型別有固定的長寬
* 因此解法為:每次 append 新的 bit 時,對於那些不足以一個 32 位元的表示的 bit 就放在 `b` 中,並用 `c` 紀錄多少部份是那些 append 的 bit
* 每當 append 滿 32 bits,就把 `b` 轉過去 `ptr` 中,並清空 b 和 c,以供新的 bit 置放(透過 `flush_buffer` 函式來完成,因為非此章重點就暫不詳細說明)
#### `write_golomb`
在壓縮時,會透過 write_golomb 來填入對應的 bit
```cpp
static void write_golomb(generic_io_t *io, size_t k, uint32_t N)
{
assert(k < 32);
write_unary(io, N >> k);
write_bits(io, N, k);
}
```
* 對於 N / $2^k$ 的商部份,使用 unary code 編碼
* 對於 N % $2^k$ ,也就是 N / $2^k$ 餘數部份,計算相應結果後寫入
#### `write_unary`
```cpp
static void write_unary(generic_io_t *io, uint32_t N)
{
for (; N > 32; N -= 32)
write_zero_bits(io, 32);
write_zero_bits(io, N);
put_nonzero_bit(io);
}
```
* `write_unary` 會根據商的結果寫入 q 個 0 加上一個 1
* `write_zero_bits` 一次最多允許寫 32 個 bits,因此如果要寫的總數超過,要透過 for 迴圈依次進行
* `put_nonzero_bit` 寫入一個 1
#### `write_bits`
```cpp
static void write_bits(generic_io_t *io, uint32_t b, size_t n)
{
assert(n <= 32);
for (int i = 0; i < 2; ++i) {
assert(io->c < 32);
size_t m = MIN(32 - io->c, n);
io->b |= (b & (((uint32_t) 1 << m) - 1)) << io->c;
io->c += m;
if (io->c == 32)
flush_buffer(io);
b >>= m;
n -= m;
if (n == 0)
return;
}
}
```
* 把 b 的後 n 個 bits 寫入,事實上就等同在取 `b % n` 的結果做編碼,但需考慮到餘數編碼可能會橫跨兩個 byte 儲存,需要兩個迴圈去做對應的處理
#### `read_golom`
反之將資料解壓縮
```cpp
static uint32_t read_golom(generic_io_t *io, size_t k)
{
uint32_t N = read_unary(io) << k;
N |= read_bits(io, k);
return N;
}
```
* 根據 `read_unary`,從 leading zero 到第一個 1 的 prefix 可以得出原本數字的商
* 結合 read_bits 得到的餘數部份,則可以解壓縮出原本的數字 `N`
下面來檢驗一下我們對程式的理解是否有錯誤:
* 假設壓縮一個文字檔,其中內容只有 `A` (ASCll = 65) 和一個換行符號 (ASCll = 10)
* 對於 65 / 8 = 8 ... 1 因此對於 `A` 的編碼為 <font color="#00f">001</font><font color="#f00">100000000</font>(注意要從右往左看)
* 對於 10 / 8 = 1 ... 2 因此對於換行符號的編碼為<font color="#00f">010</font><font color="#f00">10</font>
* 因此整合兩個結果會是 <font color="#00f">0 1010</font> <font color="#f00">0011 0000 0000</font>,因此前 4 個 byte 會是 a300,這一致於 hexdump 輸出的結果
```shell
$ hexdump output
0000000 a300 0000 0000 0002
0000008
```
* 最後還會對 EOF(ASCll = 256) 做編碼,這裡就不繼續說明
### 改寫 x-compressor 為 Exponential-Golomb coding
原本的壓縮是使用 Golomb-Rice coding,這裡我們基於原始的實作,根據維基百科上對演算法的敘述,抽換壓縮技巧為 Exponential-Golomb coding。
#### `write_expgolomb`
```cpp
static uint32_t reverse(uint32_t x)
{
x = (((x & 0xaaaaaaaa) >> 1) | ((x & 0x55555555) << 1));
x = (((x & 0xcccccccc) >> 2) | ((x & 0x33333333) << 2));
x = (((x & 0xf0f0f0f0) >> 4) | ((x & 0x0f0f0f0f) << 4));
x = (((x & 0xff00ff00) >> 8) | ((x & 0x00ff00ff) << 8));
return((x >> 16) | (x << 16));
}
static void write_expgolomb(generic_io_t *io, size_t order, uint32_t N)
{
// 1. Encode ⌊x/2k⌋ using order-0 exp-Golomb code
uint32_t B = (N >> order) + 1;
assert(B != 0);
int zero_write = 32 - __builtin_clz(B) - 1;
// since zero_write won't larger then 32, so we can write zero bits directly
write_zero_bits(io, zero_write);
int total_bits = 32 - __builtin_clz(B);
write_bits(io, reverse(B) >> (32 - total_bits), total_bits);
// 2. Encode x mod 2k in binary
uint32_t mod = N & ((1 << order) - 1);
write_bits(io, mod, order);
}
```
將 `x_compress` 中兩處的 `write_golomb` 換成 `write_expgolomb`,其中有幾點需注意:
* order 為多少才可以最好的壓縮我們的資料呢?之後會再根據文獻或者實驗去對參數的選擇做修改。
* 為使解碼正確,先將 `B` 反轉再編碼,如何反轉?編碼的部份如何有效的 bitwise 處理? 目前尚未周延考慮何種才是最佳的作法
下面來檢驗一下程式是否有錯誤(以 order 0 為例):
* 假設壓縮一個文字檔,其中內容只有 `A`(ASCll = 65) 和一個換行符號(ASCll = 10)
* 對於 65 + 1 = 0b1000010 因此對於 A 的編碼為 <font color="#00f">0100001</font><font color="#f00">000000</font>(注意要從右往左看)
* 對於 10 + 1 = 0b1011 因此對於換行符號的編碼為<font color="#00f">1101</font><font color="#f00">000</font>
* 因此整合兩個結果會是 <font color="#00f">1101 000</font><font color="#f00">0 1000 0100 0000</font>,因此前 5 個 byte 會是 d0840,這一致於 hexdump 輸出的結果
```shell
$ hexdump out
0000000 0840 000d 0000 0010
0000008
```
#### `read_expgolomb`
```cpp
static uint32_t read_expgolom(generic_io_t *io, size_t order)
{
uint32_t N = read_unary(io);
uint32_t rb =read_bits(io, N);
uint32_t num = reverse((rb << 1) | 1) >> (32 - N - 1);
rb = read_bits(io, order);
return (((num - 1) << order) | rb ) ;
}
```
解讀的部份則顯而易懂:首先藉由 0 的數量決定後續要讀出的數量,經過反轉和位移處理後則可以得到原本的數字。根據 order 的不同,後續也要做出相對應的處理。
詳細請參考 [RinHizakura/x-compressor](https://github.com/RinHizakura/x-compressor)(目前用 `make check` 測試 order = 0, 1, 2 都沒問題)
:::warning
TODO:
1. 依據資料特性,比較 Exponential-Golomb coding 及 Golomb-Rice coding 對資料壓縮的影響;
2. 參照 `git log`,原本 `x-compressor` 支援多個 layer,用以強化壓縮率
:notes: jserv
:::
## 測驗 `3`
思考一下各步驟在 bitwise 操作上的意義:
* 如果是奇數:表示最後一個位元是 1,減 1 是把最後一個 bit 設為 0
* 如果是偶數:表示最後一個位元是 0,除 2 就是向右位移
把演算法的流程圖畫出來就像這樣:
```graphviz
digraph g{
node [shape="box", fontcolor=black,fontsize=16,style=filled, fillcolor=aquamarine ]
start[shape=ellipse, fontcolor=white,style=filled, fillcolor=tomato,label="start"]
condition[label="last bit == 1?",shape="diamond", fontcolor=white, fillcolor=darkorange]
terminate[label="total bit == 1?",shape="diamond", fontcolor=white, fillcolor=darkorange]
C[label="total bit = ??",fillcolor=yellow]
{
rank = same
A[label="right shift \n(step++, total bit -= 1)"]
B[label="set last bit to 0 \n (step ++)"]
end[shape=ellipse, fontcolor=white,style=filled, fillcolor=tomato]
}
start->C->terminate
condition->B[label=" yes", fontcolor=blue,fontsize=24,color=green]
B->terminate
terminate:w->condition[label=" no", fontcolor=blue,fontsize=24]
condition:w->A[label=" no", fontcolor=blue,fontsize=24,color=red]
terminate->end[label=" yes", fontcolor=blue,fontsize=24]
A->condition
}
```
終止條件是 `total bit == 1`(因為最後一步必然是 1 變成 0,表示最後一個位元不需要做 right shift),所以紅色的 edge 會走剛好整整 `total bit - 1` 次,而走過綠色 edge 的次數則是根據整個 `num` 中有幾個 1,也就是 `popcount` 在計算的數量。至於 `total bit` 實際上是多少呢?就是表達該數字需要的最少位元數量(比如說 8 = 0b1000,所以 total bit = 4),而這個數量可以根據 32 - `clz(num)` 就可以被計算出來。
所以此題的答案是 `__builtin_popcount(num) + 32 - __builtin_clz(num) - 1` = `__builtin_popcount(num) + 31 - __builtin_clz(num)` 。`AAA` = `(a) __builtin_popcount(num)`,`BBB` = `__builtin_clz(num)`
## 測驗 `4`
根據 [Greatest Common Divisor 特性和實作考量](https://hackmd.io/@sysprog/gcd-impl) 中的 [Binary GCD Algorithm](https://iq.opengenus.org/binary-gcd-algorithm/) 來閱讀程式。
* `for (shift = 0; !((u | v) & 1); shift++)` 迴圈在判斷是否 u 和 v 皆為偶數,如果 u 和 v 皆為偶數,則 gcd(u, v) = 2 * gcd(u/2, v/2)
* 經過前面步驟的 u 和 v 的數值可能會改變(假設是 $u'$、$v'$),所以接下來要先解 gcd($u'$, $v'$),最後的結果再位移先前除 2 的次數
* `while (!(u & 1))` 在計算: 如果 $u'$ 為偶數,則 gcd($u'$, $v'$) = gcd($u'$/2, $v'$)
* `do` 迴圈的第一個 iteration,`while (!(v & 1))` 在計算: 如果 $v'$ 為偶數,則 gcd($u'$, $v'$) = gcd($u', v'$/2)
* 然後剩下的其實就是輾轉相除法(Euclidean algorithm),並且會確保進入下一個迴圈的 v' 是比較小的那個值,重複直到**進入這次輾轉相除時的 u 和 v 相同**,因此 `XXX` = `(b) v`(因為在判斷式時 v 已經是前面計算的 t)
* 則這個相同的數字就是最大公因數,也就是 $u'$,再位移最開始除 2 的次數,所以 `YYY` = `(e) u << shift`
### 透過 `__builtin_ctz` 改寫
注意到下面的實驗中會將 quiz3 中題目的作法命名為 `fast gcd`,使用 `__builtin_ctz` 調整後的版本則稱為 `faster gcd` 以方便解釋。
前面的 `for (shift = 0; !((u | v) & 1); shift++)` 迴圈可以使用 `__builtin_ctz` 來改寫為如下:
```cpp
#ifndef min
#define min(x, y) ((x) < (y) ? (x) : (y))
#endif
uint64_t faster_gcd64(uint64_t u, uint64_t v){
if (!u || !v)
return u | v;
int shift = min(__builtin_ctz(u), __builtin_ctz(v));
u >>= shift; v>>= shift;
while (!(u & 1))
u >>= 1;
do {
while (!(v & 1))
v >>= 1;
if (u < v) {
v -= u;
} else {
uint64_t t = u - v;
u = v;
v = t;
}
} while (v);
return u << shift;
}
```
則可以將 for 迴圈中的操作化簡,減少程式分支和對此操作的總指令數量。而為了比較兩者的效能區別而設計對應的實驗:對於相同 ctz 的 a 和 b,進行 50 次實驗並測量時間,注意一些細節:
* 為了確保平等的比較兩個方法,原本的 fast gcd 中的 `/ 2` 都替換成 `<< 1` (雖然經過最佳化後可能沒有差別)
* 這裡我使用了一些技巧來確保 a 和 b 可以有固定數量的 trailing zero,當然時間的比較會基於相同 a、b 來計算
* 此外,這個實驗中 a 和 b 會是同樣數量的 trailing zero,但是也該考慮 trailing zero 的不同排列組合來設計實驗(暫時還沒想到如何視覺化可以比較好的觀察結果)
實驗程式以最佳化參數 `-O2` 進行編譯,初步我嘗試用 3 維的圖形來觀察希望可以更好的視覺化實驗結果,不過似乎不是那麼理想QQ 但還是可以從圖中觀察觀察一些資訊
* 紅色的點(透過 `__builtin_ctz` 改寫的版本) 整體上會有較佳的執行時間
* 隨著數字的 ctz 越大,Euclidean algorithm 的操作次數就有機會變得更少,所以 y 軸是向 ctz of and b = 30 的地方傾斜的
![](https://i.imgur.com/TKTNFBt.png)
重新設計實驗:對兩個方法進行 1000 次的實驗,排除標準差 2 外的 outlier 後繪制 boxplot 來視覺化實驗結果,則可以進一步看出其中差距
![](https://i.imgur.com/wZgxS3Y.png)
可以看到藍色的 fasted_gcd 執行的效率較佳,同時也可以看出前面所說時間會往 ctz = 30 傾斜的狀況。
## 測驗 `5`
`bitmap` 中有 `bitmapsize` 個 64 位元的整數,因此可以用來表示一串 64 * `bitmapsize` 的數字。而 `naive` 程式就是在計算這些數字中哪些位置為 1,儲存在 `out` 中,而 `pos` 則是結果為 1 的總數。
在 `naive` 的解法中,我們會逐一檢查每一個 bit 是否為 1,去紀錄該位置。但事實上,如果我們把問題轉換成: **找出最右邊的 1**(顯然可以透過前面的 `ctz` 做到) + 將最右邊的 1 設成 0,則重複此步驟,我們就可以不必依序檢查每個 bits,得到相同的結果。
如何將最右邊的 1 設成 0 呢?我們知道如果 A & $\bar{A}$ = 0(每個 bit 都相異),$A$ & $(\bar{A} + 1)$ 則會得到和 A 最右邊的 1 相同位置為 1,其他為 0 的 mask `t`,與此 mask 做 XOR 運算就可以消除最右邊的 1 。因此此題答案為:`KKK` = `(e) bitset & -bitset`。