# 2020q3 Homework3 (quiz3)
contributed by < `nelsonlai1` >
## 測驗 1
```c=
int asr_i(signed int m, unsigned int n)
{
const int logical = (((int) -1) OP1) > 0;
unsigned int fixu = -(logical & (OP2));
int fix = *(int *) &fixu;
return (m >> n) | (fix ^ (fix >> n));
}
```
這題要實作在任何環境下可運作的 ars(arithmetic right shift)，ars 會在右移後，依照原本首位的 bit 補回值，也就是正數補 0，負數補 1。
一開始的 logical 判斷環境是不是 logical right shift (一律補 0)，所以對 -1 右移 1(`OP1 = >> 1`)，如果是 lrs，會變成 0b011...111 > 0 = true = 1，反之則是 0b111...111 < 0 = false = 0。
fixu 用來判斷被除數 m 是不是負數(`OP2 = m < 0`)，當 logical = 1 以及 m < 0 時，fixu = -(1 & 1) = 0b111...111，其他時候fixu = 0b000...000。fix 則是將 fixu(unsigned int) 的值存成 int。
最後 return 的結果，如果在 logical = 1 以及 m < 0 都成立時，fix ^ (fix >> n) 會是前 n 個 bits 為 1，其餘為 0。所以 | (fix ^ (fix >> n)) 就可以把 lrs 吃掉的 1 補回來。
example : m = 0b10110101, n = 4, system = logical
```graphviz
digraph {
m [shape = plaintext]
intm [label = "1|0|1|1|0|1|0|1" shape = record]
}
```
```graphviz
digraph {
m_shift [label = "m >> n" shape = plaintext]
intm_shift [label = "0|0|0|0|1|0|1|1" shape = record]
}
```
```graphviz
digraph {
fix [shape = plaintext]
intfix [label = "1|1|1|1|1|1|1|1" shape = record]
}
```
```graphviz
digraph {
fix_shift [label = "fix ^ (fix >> n)" shape = plaintext]
intfix_shift [label = "1|1|1|1|0|0|0|0" shape = record]
}
```
```graphviz
digraph {
return [shape = plaintext]
intreturn [label = "1|1|1|1|1|0|1|1" shape = record]
}
```
## 測驗 2
```c=
bool isPowerOfFour(int num)
{
return num > 0 && (num & (num - 1))==0 &&
!(__builtin_ctz(num) OPQ);
}
```
這題要輸入的數字是否為 4 的次方，4 的二進位表示是 100，也就是說 4 的次方一定是 1 後面接著偶數個 0 (每次乘四都是 << 2)。
一開始 `num > 0` 判斷是否為正整數，`(num & (num - 1)) == 0` 判斷是否為 2 的次方
($1{\overbrace{000...000}^{n個}}$ & $0{\overbrace{111...111}^{n個}} = 0$) ，所以最後這項 `!(__builtin_ctz(num) OPQ` 可以知道是在判斷後面的 0 的個數是否為偶數了。而任一偶數 & 1 = 0 (偶數最後一項 bit 為 0)，所以 `OPQ = & 0x1`。
### 延伸問題
1. 改寫上述程式碼，提供等價功能的不同實作，儘量降低 branch 的數量
```c=
bool isPowerOfFour(int num){
return num & 0x55555555 && (num & (num - 1)) == 0;
}
```
& 0x55555555 可以一次判斷是否為正數以及 trailing 0-bits 是否為偶數(在 num 是 2 的次方條件下)
2. 練習 LeetCode 1009. Complement of Base 10 Integer 和 41. First Missing Positive，應善用 clz
Complement of Base 10 Integer :
```c=
int bitwiseComplement(int N){
if (N == 0)
return 1;
int shift = __builtin_clz(N);
return ~((uint32_t)N << shift) >> shift;
}
```
利用右移會無條件捨去的特性，先將輸入值左移到最左邊後做 not，再右移一樣的次數(`__builtin_clz(N)`)，就能得到結果。
First Missing Positive :
```c=
int firstMissingPositive(int* nums, int numsSize)
{
for (int i = 0; i < numsSize; i++) {
if (nums[i] > 0 && nums[i] <= numsSize) {
int swap = nums[i] - 1;
if (nums[i] != nums[swap]) {
int tmp = nums[i];
nums[i] = nums[swap];
nums[swap] = tmp;
i--;
}
}
}
for (int i = 0; i < numsSize; i++) {
if (nums[i] != i + 1)
return i + 1;
}
return numsSize + 1;
}
```
用 bucket sort 的概念，將 1, 2, 3... 換到 nums[0], nums[1], nums[2]...。最後再從 nums[0] 開始找起，第一個不符合的數就是答案，如果都符合的話答案為 numsSize + 1 (代表從 1 到 numsSize 都各出現一次)。
值得一提的是，如果直接寫 nums[nums[i] - 1]，會出現 `ERROR: AddressSanitizer: heap-buffer-overflow` 錯誤，推測是因為數據中有小於 1 的數，即使不會執行到這行，系統也會判斷操作到錯誤的 array 位置 (< 0)。
## 測驗 3
```c=
int numberOfSteps (int num)
{
return num ? AAA + 31 - BBB : 0;
}
```
這題計算一個數經過多少次計算可以變為零，而計算只有分兩種 : >> 1 (當前是偶數)以及 - 1 (當前是奇數)。
int 寬度為 32 bits，因此將最左邊的 bit 移到最右邊要 31 次 >> 1，每有一個 1 都需要做一次 - 1，而起始 1 的位置每往右一個 (leading 0-bits 多一)就能少做一次 >> 1。
所以 `AAA = __builtin_popcount(num)`, `BBB = __builtin_clz(num)`
### 延伸問題
1. 改寫上述程式碼，提供等價功能的不同實作並解說
```c=
int numberOfSteps (int num){
return __builtin_popcount(num) + 31 - __builtin_clz(num | 1);
}
```
因為 `__builtin_clz(0)` 不被允許，所以加上 `| 1` 讓即使在 `num == 0` 時也能返回正確值，同時也不會影響其他情況的結果(因為是最後一位)。
2. 避免用到編譯器的擴充功能，只用 C99 特徵及 bitwise 運算改寫 LeetCode 1342. Number of Steps to Reduce a Number to Zero，實作出 branchless 解法，儘量縮減程式碼行數和分支數量
```c=
int popcount(int x)
{
uint32_t v = (uint32_t) x;
uint32_t n;
n = (v >> 1) & 0x77777777;
v -= n;
n = (n >> 1) & 0x77777777;
v -= n;
n = (n >> 1) & 0x77777777;
v -= n;
v = (v + (v >> 4)) & 0x0F0F0F0F;
v *= 0x01010101;
return v >> 24;
}
int clz(int x)
{
x |= x >> 1;
x |= x >> 2;
x |= x >> 4;
x |= x >> 8;
x |= x >> 16;
return 32 - popcount(x);
}
int numberOfSteps (int num){
return popcount(num) + 31 - clz(num | 1);
}
```
popcount 的寫法參照測驗 3 題目敘述的 popcnt_branchless，而 clz 寫法前面五行可以保證 MSB 後的 bits 皆為 1。最後再用 32(total bits) 減去 popcount 就是 leading 0-bits 數量。
(做法參考 [The Aggregate Magic Algorithms](http://aggregate.org/MAGIC/#Leading%20Zero%20Count))
## 測驗 4
64-bit GCD (greatest common divisor, 最大公因數) 求值函式
```c=
#include <stdint.h>
uint64_t gcd64(uint64_t u, uint64_t v) {
if (!u || !v) return u | v;
while (v) {
uint64_t t = v;
v = u % v;
u = t;
}
return u;
}
```
改寫為以下等價實作
```c=
#include <stdint.h>
uint64_t gcd64(uint64_t u, uint64_t v) {
if (!u || !v) return u | v;
int shift;
for (shift = 0; !((u | v) & 1); shift++) {
u /= 2, v /= 2;
}
while (!(u & 1))
u /= 2;
do {
while (!(v & 1))
v /= 2;
if (u < v) {
v -= u;
} else {
uint64_t t = u - v;
u = v;
v = t;
}
} while (XXX);
return YYY;
}
```
上面的 function 在實作輾轉相除法，一開始先判斷如果 u, v 有任一數不存在，就回傳存在的那個數。之後 while 迴圈裡在做的是將較大的一方 u 對較小的一方 v 取餘數(即使一開始 u 比較小，過完一輪後就會跟 v 互換值)，將 v 改為餘數值(小)，將 u 改為原本的 v (大)，算是比較直觀的寫法。
下面的做法可以搭配 [Binary GCD](https://hackmd.io/@sysprog/gcd-impl#Binary-GCD) 來思考。第一個 for 迴圈的內容就是這句的應用
>binary gcd 的精神就是當兩數為偶數時，必有一 2 因子。
而 shift 用來紀錄有幾個 2 因子，所以最後 return 的結果應該要 << shift ($\times 2^{shift}$) 補回來。
之後分別檢查 u, v 是否偶數則是對應這句
>且一數為奇數另一數為偶數，則無 2 因子。
所以透過不斷除以二，直到兩者都是奇數為止。
最後就是做一般的輾轉相除法，只是這次不使用 % 運算子，所以就是比較兩數大小，然後大減小，變得比較像是 [Wikipedia 的展示動畫](https://zh.wikipedia.org/zh-tw/%E8%BC%BE%E8%BD%89%E7%9B%B8%E9%99%A4%E6%B3%95#/media/File:Euclidean_algorithm_252_105_animation_flipped.gif)，而我們也可以發現 v 都是被減數，所以`XXX = v`。
而 return 結果則是依照 $\text{gcd}(x, y) = \text{even_factor}·\text{gcd}((x\gg \text{even_factor}), (y\gg \text{even_factor}))$，把 $\text {even_factor}(2^{shift})$ 乘回來，所以 `YYY = u << shift`。
### 延伸問題
1. 在 x86_64 上透過 __builtin_ctz 改寫 GCD，分析對效能的提升
```c=
uint64_t gcd64(uint64_t u, uint64_t v) {
if (!u || !v) return u | v;
int ctzu = __builtin_ctz(u);
int ctzv = __builtin_ctz(v);
int shift = ctzv ^ ((ctzu ^ ctzv) & -(ctzu < ctzv));
u = u >> ctzu;
v = v >> ctzv;
do {
while (!(v & 1))
v /= 2;
if (u < v) {
v -= u;
} else {
uint64_t t = u - v;
u = v;
v = t;
}
} while (v);
return u << shift;
}
```
第五行使用 [Bit Twiddling Hacks](http://graphics.stanford.edu/~seander/bithacks.html#IntegerMinOrMax) 的方法得到兩數中較小的數。
利用 perf 檢查兩者效能差距 :
由於單一次呼叫 function 不足以有效看出差距，所以使用 for 迴圈來多次呼叫，且因為主要改善在
u, v 皆為偶數時，故只測 u, v 皆為偶數的情況。
```c=
for (uint64_t i = 0; i <= 0x100000;i += 2)
gcd64(i, 0x100000);
```
原始 :
```
Performance counter stats for './quiz3_4_control':
57.54 msec task-clock # 0.995 CPUs utilized
1 context-switches # 0.017 K/sec
0 cpu-migrations # 0.000 K/sec
46 page-faults # 0.799 K/sec
2,3809,7810 cycles # 4.138 GHz
2,2032,2845 instructions # 0.93 insn per cycle
4311,6621 branches # 749.310 M/sec
332,9927 branch-misses # 7.72% of all branches
0.057811581 seconds time elapsed
0.057830000 seconds user
0.000000000 seconds sys
```
使用 __builtin_ctz 改寫 :
```
Performance counter stats for './quiz3_4_test':
40.82 msec task-clock # 0.996 CPUs utilized
0 context-switches # 0.000 K/sec
0 cpu-migrations # 0.000 K/sec
45 page-faults # 0.001 M/sec
1,6873,4941 cycles # 4.134 GHz
1,4741,2793 instructions # 0.87 insn per cycle
3052,7407 branches # 747.908 M/sec
261,0331 branch-misses # 8.55% of all branches
0.040974330 seconds time elapsed
0.041005000 seconds user
0.000000000 seconds sys
```
可以發現在速度上以及 branches 有大約 30% 的改善。若是 u, v 有很大的 2^N 公因數時，效能差距可以更為顯著。
## 測驗 5
在影像處理中，bit array (也稱 bitset) 廣泛使用，考慮以下程式碼:
```c=
#include <stddef.h>
size_t naive(uint64_t *bitmap, size_t bitmapsize, uint32_t *out)
{
size_t pos = 0;
for (size_t k = 0; k < bitmapsize; ++k) {
uint64_t bitset = bitmap[k];
size_t p = k * 64;
for (int i = 0; i < 64; i++) {
if ((bitset >> i) & 0x1)
out[pos++] = p + i;
}
}
return pos;
}
```
可用 clz 改寫為效率更高且等價的程式碼:
```c=
#include <stddef.h>
size_t improved(uint64_t *bitmap, size_t bitmapsize, uint32_t *out)
{
size_t pos = 0;
uint64_t bitset;
for (size_t k = 0; k < bitmapsize; ++k) {
bitset = bitmap[k];
while (bitset != 0) {
uint64_t t = KKK;
int r = __builtin_ctzll(bitset);
out[pos++] = k * 64 + r;
bitset ^= t;
}
}
return pos;
}
```
一個 bitmap 長得像這樣，共有 64 * bitmapsize 個 bits。
for 迴圈裡做的是找出那些 bits 為 1 (`if ((bitset >> i) & 0x1)`)，pos 紀錄 1 的個數，out 紀錄 1 的位置(`out[pos++] = p + i`)。
| | 0 | 1 | 2 | ... | bitmapsize - 1 |
| --- | --- | --- | --- | --- | -------------- |
| 0 | 0 | 1 | 0 | ... | 1 |
| 1 | 0 | 0 | 1 | ... | 0 |
| 2 | 1 | 1 | 1 | ... | 1 |
| ... | ... | ... | ... | ... | ... |
| 63 | 0 | 1 | 0 | ... | 1 |
下方的 code 改寫了後面尋找 1 的部分，從原本逐個 bit 測試，改成直接計算有幾個 trailing 0-bits，再把最後的 1 改成 0，進入下一輪。
而可以做到的 t 是 `bitset & -bitset`，因為 -bitset = ~bitset + 1 (2's complement)，所以 t 會是 bitset 最後一個 1 的位置為 1，而其餘皆是 0 (獨立 LSB)。而 `bitset ^= t` 就能把最後一個 1 改成 0 了。
### 延伸問題
1. 舉出這樣的程式碼用在哪些真實案例中
2. 設計實驗，檢驗 clz 改寫的程式碼相較原本的實作有多少改進？應考慮到不同的 bitmap density
3. 思考進一步的改進空間
在[這篇文章](https://lemire.me/blog/2018/03/08/iterating-over-set-bits-quickly-simd-edition/)中有提到使用 SIMD 指令改寫的方法，在 bitmap 密集時，有明顯的效能增長。
Sign in with Wallet
Connect another wallet