---
title: '2020q3 Homework3 (quiz3)'
tags: linux2020
---
# 2020q3 Homework3 (quiz3)
contributed by < `ChongMingWei` >
## 開發環境
```shell
$ uname -a
Linux cmw-System-Product-Name 5.4.0-47-generic #51~18.04.1-Ubuntu SMP Sat Sep 5 14:35:50 UTC 2020 x86_64 x86_64 x86_64 GNU/Linux
$ gcc --version
gcc (Ubuntu 7.5.0-3ubuntu1~18.04) 7.5.0
Copyright (C) 2017 Free Software Foundation, Inc.
This is free software; see the source for copying conditions. There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
```
---
## 測驗 1
In some environment, we don't have arithmetic right shift. That is, negative number after `>>` operation will be positive.
So, we have to set bits which become `0` after the `>>` operation.
```c=
int asr_i(signed int m, unsigned int n)
{
const int logical = (((int) -1) >> 1) > 0;//OP1
unsigned int fixu = -(logical & (m<0));//OP2
int fix = *(int *) &fixu;
return (m >> n) | (fix ^ (fix >> n));
}
```
Check if the compiler supports arithmetic right shift.
If yes, `logical = 0`. Otherwise, `logical = 1`
```c=3
const int logical = (((int) -1) ) > 0;//OP1
```
`fixu` will be assigned `0xffffffff` if `logical = 1` and `m<0` (The environment doesn't support arithmetic right shift.)
(`fixu = 4294967295`)
```c=4
unsigned int fixu = -(logical & (m<0));//OP2
```
`fix` will also be assigned `0xffffffff`.
(With the datatype `int`, `fix = -1`)
```c=5
int fix = *(int *) &fixu;
```
`m >> n` will get n `0` bits. (From left to right)
`fix >> n` will also get n `0` bits. (From left to right)
`fix ^ (fix >> n)` will reverse all bits. And we will get n `1` bits(From left to right) and 32 - n `0` bits.
```c=6
return (m >> n) | (fix ^ (fix >> n));
```
:::success
2.練習實作其他資料寬度的 ASR，可參照[ C 語言：前置處理器應用篇 ](https://hackmd.io/@sysprog/c-preprocessor)撰寫通用的巨集以強化程式碼的共用程度;
:::
將原本 `int asr_i` 裡面的程式碼都改寫成 macro 並且能給予不同資料型態。
```c=
#define get_logical(datatype) (((datatype) -1) >> 1) > 0
#define get_fixu(logical, m) -(logical & (m < 0))
#define get_fix(datatype, fixu) (*(datatype *) &(fixu))
#define get_result(fix, m, n) (m >> n) | (fix ^ (fix>>n))
```
## 測驗 2
```cpp
bool isPowerOfFour(int num)
{
return num > 0 && (num & (num - 1))==0 &&
!(__builtin_ctz(num) & 0x1); //OPQ
}
```
### Analysis on return value
- 考慮是否為正數
```cpp
num > 0
```
- 考慮是否有多個 set bits (應當只有一個)，也就是確認是否為` power of 2`。
> 但此種寫法會把0也算進去，應替換成 `num && !(num & (num - 1))`
```cpp
num & (num - 1))==0
```
- 考慮 set bit 是否在正確的位置上(應當要在 bit(2n) 的位置上):
- $4^n = 2^{2n}$，$2n$ 即為 `__builtin_ctz(num)` 會返回的數字，可以發現，只要**返回的數字是大於等於0的偶數**，就會是我們要的答案。
- 因此只要**返回的數字是大於等於2的偶數**，就會是我們要的答案，用 `&0x1` 來做檢查即可。
```cpp
__builtin_ctz(num) & 0x1
```
>測驗時沒有看到 `!` 且只想到要把 bit0 轉成1 (但是前面的 bits 也可能會有值!!)，所以挑了 (g) 選項
:::success
2.改寫上述程式碼，提供等價功能的不同實作，儘量降低 branch 的數量;
:::
- `&& __builtin_popcount(num)==1`: 確認只有一個 set bit。
- `(num&0x55555555)`: 確認 set bit 在正確的位置上。
```c=
bool isPowerOfFour_m(int num)
{
return num > 0 && __builtin_popcount(num)==1 && (num&0x55555555);
}
```
:::success
3.練習[ LeetCode 1009. Complement of Base 10 Integer ](https://leetcode.com/problems/complement-of-base-10-integer/)和[ 41. First Missing Positive](https://leetcode.com/problems/first-missing-positive/)，應善用 clz;
:::
### 1009. Complement of Base 10 Integer
- 想法: 將原本的 input 做 bit inverse ，但使用的是 int 的資料型態，所以做完 inverse 之後，要將和原本 input 無關的 bits 都設回0。
- 如果是0就返回1
- 如果是其他數字，計算 mask: `~((0xfffffffe) << (31-__builtin_clz(N)))` ，
假設現在 N = 5，那麼 mask 會得到`0x000000007` ，
使用此 mask 和 `~N` 做 `AND` 操作，可以將和 input 無關的 bits 設回0。
```c=
int bitwiseComplement(int N){
return N==0?1:~N&~((0xfffffffe) << (31-__builtin_clz(N)));
}
```
### 41. First Missing Positive
> 參考 [LEETCODE 41. First Missing Positive 解题思路分析](https://leetcode.jp/leetcode-41-first-missing-positive-%E8%A7%A3%E9%A2%98%E6%80%9D%E8%B7%AF%E5%88%86%E6%9E%90/)
- 想法: 假設 ^1^ `nums` 裡面的數字是1- `numsSize` ，那麼要回傳 `numsSize+1`，否則^2^回傳的數字一定會落在[1, `numsSize`]的區間內。
- 第一個 for loop: 所以使用一個長度為 `numsSize` 的 array ，若 `index+1` 的數字在 `nums` 中有出現過，那就將 array[index] 設為1。
- 第二個 for loop: 檢查 array ，第一個出現不為1的 index ，將其+1後就是我們要求的數字，如果都有1的話，那就是前面提到的第一種情況，回傳 `numsSize+1`
```c=
int firstMissingPositive(int* nums, int numsSize){
int *tmp = malloc(sizeof(int)*numsSize);
for (int i = 0;i < numsSize; ++i){
if(nums[i]>0&&nums[i]<=numsSize)
tmp[nums[i]-1] = 1;
}
for (int i = 0;i < numsSize; ++i){
if(tmp[i]!=1)
return i+1;
}
return numsSize+1;
}
```
:::success
4.研讀[ 2017 年修課學生報告](https://hackmd.io/@3xOSPTI6QMGdj6jgMMe08w/Bk-uxCYxz)，理解 clz 的實作方式，並舉出[ Exponential Golomb coding ](https://en.wikipedia.org/wiki/Exponential-Golomb_coding)的案例說明，需要有對應的 C 原始程式碼和測試報告;
x-compressor 是個以 Golomb-Rice coding 為基礎的資料壓縮器，實作中也用到 clz/ctz 指令，可參見[ Selecting the Golomb Parameter in Rice Coding](https://ipnpr.jpl.nasa.gov/progress_report/42-159/159E.pdf)。
> [x-compressor](https://github.com/xbarin02/x-compressor)
:::
## 測驗 3
- 如果是0，直接返回數值0
- 如果非0，那麼就要計算減1的次數和除以2的次數，其中前者和後者分別為:
- 1的總個數 `__builtin_popcount(num)`，以及
- MSB 之後的 bits 數量 `31 - __builtin_clz(num)`
```cpp
int numberOfSteps (int num)
{
return num ? __builtin_popcount(num) + 31 - __builtin_clz(num) : 0; //AAA BBB
}
```
:::success
2.改寫上述程式碼，提供等價功能的不同實作並解說;
提示:[ Bit Twiddling Hacks ](http://graphics.stanford.edu/~seander/bithacks.html)中提及許多[ bitwise operation ](https://hackmd.io/@sysprog/ByzoiggIb)的使用，如 bit inverse、abs 等等，是極好的參考材料。
:::
> 參考 [2017q3 Homework4 (改善 clz)](https://hackmd.io/@3xOSPTI6QMGdj6jgMMe08w/Bk-uxCYx)
```c=
int numberOfSteps (int num)
{
/* Compute population count */
uint32_t popcount = num;
uint32_t tmp;
tmp = (popcount>>1)&0x77777777;
popcount -= tmp;
tmp = (tmp>>1)&0x77777777;
popcount -= tmp;
tmp = (tmp>>1)&0x77777777;
popcount -= tmp;
popcount += (popcount >> 4);
popcount &= 0x0f0f0f0f;
popcount *= 0x01010101;
popcount >>= 24;
/* Compute clz */
uint32_t n = 32, c = 16, l, r = num;
do {
l = r >> c;
if (l) {
n -= c;
r = l;
}
c >>= 1;
} while (c);
return num ? popcount + 31 - (n - r) : 0; //AAA BBB
}
```
:::success
3.避免用到編譯器的擴充功能，只用 C99 特徵及 bitwise 運算改寫[ LeetCode 1342. Number of Steps to Reduce a Number to Zero](https://leetcode.com/problems/number-of-steps-to-reduce-a-number-to-zero/)，實作出 branchless 解法，儘量縮減程式碼行數和分支數量;
:::
```c=
int numberOfSteps (int num){
int numofminus=0;//要做-1的次數
int numofdiv=0;
while(num){
int isodd = num&0x1;//確認目前 bit0 是否為1
int iseven = (num&0x1)^0x1;//確認目前 bit0 是否為0
numofminus += isodd;
numofdiv += iseven;
num -= (isodd-iseven)>0;//是奇數做-1
num >>= (iseven-isodd)>0;//是偶數做/2
}
return numofminus+numofdiv;
}
```
## 測驗 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_modified(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 (v); //XXX
return u<<shift; //YYY
}
```
- 檢查 `u` 和 `v` 是否都是偶數，直到一方不為偶數為止，可以取得最後要**向左 shift 回來的次數**
> E.g.
> A = 2^n^\*a, B = 2^n+m^\*b, a and b are odd
> GCD(A, B) = GCD(a, 2^m^\*b) * 2^n^
> n即為此處計算的 **shift**
```cpp=5
for (shift = 0; !((u | v) & 1); shift++) {
u /= 2, v /= 2;
}
```
- 如果 `u` 還是偶數，那就把他除以2直到他變奇數為止，因為當一數為奇數且另一數為偶數時，偶數中不管有幾個因數2，都不會成為兩數最大公因數的因數
> E.g.
> A = 2^n^\*a, B = b, a and b are odd
> GCD(A, B) = GCD(2^n^\*a, b) = GCD(a, b)
> n即為此處計算的 **shift**
```cpp=8
while (!(u & 1))
u /= 2;
```
- Line 11-12: 和上述 `u` 相同的操作，確認 `v` 為奇數
- Line 13-14: 使用了性質 `GCD(A, B) = GCD(A, B-A) (when B > A)`
> 證明參考:
> https://blog.csdn.net/CS171_chengl/article/details/104967588
> Proof 1:
> **If GCD(A, B) = 1, then GCD(A, B - A) = 1**
> Because A and B are coprime, B = q * A + r (r not equal to 0)
> So, B - A = (q - 1) * A + r
> Finally, GCD(A, B - A) = 1.
> Proof 2:
> **GCD(A, B) = GCD(A, B-A)**
> Suppose GCD(A, B) = d
> A = A' * d, B = B' * d, B-A = (B' - A') * d
> GCD(A, B-A) = GCD(A' * d, (B' - A') * d) = d * GCD(A', (B' - A'))
> Because **GCD(A', (B' - A')) = 1** (*by Proof 1*), d * GCD(A', (B' - A')) = d
> So GCD(A, B) = d = GCD(A, B-A)
- Line 15-18: 使用上述相同性質，只是我們需要 `v=0` 來做為迴圈跳出條件，所以要交換位置
```cpp=10
do {
while (!(v & 1))
v /= 2;
if (u < v) {
v -= u;
} else {
uint64_t t = u - v;
u = v;
v = t;
}
} while (v);
```
- `u` 最後要再左移 `shift`
```cpp=21
return u<<shift;
```
:::success
在 x86_64 上透過 `__builtin_ctz` 改寫 GCD，分析對效能的提升;
:::
```c=
uint64_t gcd64_improved(uint64_t u, uint64_t v) {
if (!u || !v) return u | v;
int shift;
shift = __builtin_ctz(u) > __builtin_ctz(v)?__builtin_ctz(v):__builtin_ctz(u);
u >>= shift;
v >>= shift;
u >>= __builtin_ctz(u);
do {
v >>= __builtin_ctz(v);
if (u < v) {
v -= u;
} else {
uint64_t t = u - v;
u = v;
v = t;
}
} while (v);
return u<<shift;
}
```
### 效能分析
實驗設計三種演算法: `gcd64`, `gcd64_modified`, `gcd64_improved` ，分別對相同的兩數做計算，兩數使用亂數產生，計算取 1000 次亂數後得到的平均執行時間如下:
| `gcd64` | `gcd64_modified` | `gcd64_improved` |
|:-------:|:----------------:|:----------------:|
| 1022 ns | 1727 ns | 915 ns |
## 測驗 5
訊號 or 影像處理會用到 (平行化?)
在影像處理中，[bit array](https://en.wikipedia.org/wiki/Bit_array) (也稱 bitset) 廣泛使用，考慮以下程式碼:
>此程式碼是用來計算，`bitmap` 中，set bit 的總個數
```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 改寫為效率更高且等價的程式碼:
```cpp=
#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 = bitset & -bitset;//KKK
int r = __builtin_ctzll(bitset);
out[pos++] = k * 64 + r;
bitset ^= t;
}
}
return pos;
}
```
原本的for 迴圈，需要遍歷64個 bit，修改為以下的while 迴圈之後，看`bitmap` 內set bit 有幾個，就跑幾次迴圈就好。
在 Line 9 的地方，我們使用變數`t` 來存放 LSB，經過 `bitset & -bitset` 的操作之後，只會留下 LSB ，其他的 bit 都會被設為 `0`，如此一來，就可以使用 `XOR` 的操作來將當前的 LSB 消除。 (Line 12)
```c=8
while (bitset != 0) {
uint64_t t = bitset & -bitset;//KKK
int r = __builtin_ctzll(bitset);
out[pos++] = k * 64 + r;
bitset ^= t;
}
```
:::success
延伸問題:
2.設計實驗，檢驗 clz 改寫的程式碼相較原本的實作有多少改進？應考慮到不同的 [bitmap density](http://www.vki.com/2013/Support/docs/vgltools-3.html)
:::
### 實驗設計
- 根據不同的 bitmap density ，計算分別需要多少時間:
這邊的密度是以32-bit 為一組單位，而我們使用的數字資料型態為 `uint64_t` ，測試的時候會重複兩次，換算成16進位數字如下:
| Bitmap density | Hex |
| -------- | -------- |
| `BITMAP_DOT` | `0xc0c0c0c0c0c0c0c0` |
| `BITMAP_DASH` | `0xf0f0f0f0f0f0f0f0` |
| `BITMAP_DOTDASH` | `0xf0c0f0c0f0c0f0c0` |
| `BITMAP_LDASH` | `0xfcfcfcfcfcfcfcfc` |
| `BITMAP_DOTLDASH` | `0xfc30fc30fc30fc30` |
| `BITMAP_DOTDOT` | `0xcccccccccccccccc` |
| `BITMAP_DOTDOTLDASH`| `0xfcccfcccfcccfccc` |
| `BITMAP_LLDASH` | `0xfff0fff0fff0fff0` |
- 執行時間結果如下:
| Bitmap density | Original (unit:ns) | Modified (unit:ns)|Speed up|
| -------- | -------- | -------- |--------|
| `BITMAP_DOT` | 1233 | 565 |2.18x|
| `BITMAP_DASH` | 1877 | 865 |2.17x|
| `BITMAP_DOTDASH` | 1651 | 736 |2.24x|
| `BITMAP_LDASH` | 1555 | 1066|1.16x|
| `BITMAP_DOTLDASH` | 1785 | 776 |2.3x |
| `BITMAP_DOTDOT` | 1083 | 556 |1.95x|
| `BITMAP_DOTDOTLDASH`| 932 | 767 |1.22x|
| `BITMAP_LLDASH` | 1240 | 880 |1.41x|
:::success
延伸問題:
3.思考進一步的改進空間;
:::
### 平行化設計
> 參考 [簡易 Pthreads 平行化範例與效能分析](https://tigercosmos.xyz/post/2020/07/simple-pthread-usage/)
以下程式碼，對 bitset 的部分做平行化。
```c=
//TODO 考慮bitset 數量不能被thread 數量整除的情況
#include <stddef.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <pthread.h>
#define NUMTHRDS 1
#define MAGNIFICATION 300
//先考慮整除的情況，每個thread分配到的都是整數
uint64_t *a;
uint32_t *b;
typedef struct
{
int thread_id;
int start;
int *pos;
int bitmapsize;
} Arg; // 傳入 thread 的參數型別
pthread_t callThd[NUMTHRDS]; // 宣告建立 pthread
pthread_mutex_t mutexsum; // pthread 互斥鎖
// 每個 thread 要做的任務
void *people_count(void *arg)
{
Arg *data = (Arg *)arg;
int thread_id = data->thread_id;
int start = data->start;
int bitmapsize = data->bitmapsize;
int *pos = data->pos;
uint64_t bitset;
for (size_t k = start; k < bitmapsize; ++k) {
bitset = a[k];
while (bitset != 0) {
uint64_t t = bitset & -bitset;//KKK
int r = __builtin_ctzll(bitset);
pthread_mutex_lock(&mutexsum);
b[*pos] = k * 64 + r;
*pos=-(~*pos);//*pos = *pos + 1
pthread_mutex_unlock(&mutexsum);
bitset ^= t;
}
}
pthread_exit((void *)0);
}
int main()
{
a = malloc(sizeof(uint64_t)*MAGNIFICATION);
for(int i =0;i<MAGNIFICATION;++i)
a[i] = 2*i+5;
b = malloc(sizeof(uint32_t)*2048);
// 初始化互斥鎖
pthread_mutex_init(&mutexsum, NULL);
// 設定 pthread 性質是要能 join
pthread_attr_t attr;
pthread_attr_init(&attr);
pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
// 每個 thread 都可以存取的 PI
// 因為不同 thread 都要能存取，故用指標
int *pos = malloc(sizeof(int));
*pos = 0;
int part = MAGNIFICATION / NUMTHRDS;
Arg arg[NUMTHRDS]; // 每個 thread 傳入的參數
for (int i = 0; i < NUMTHRDS; i++){
// 設定傳入參數
arg[i].thread_id = i;
arg[i].start = part * i;
arg[i].pos = pos; // PI 的指標，所有 thread 共用
arg[i].bitmapsize = part * i + part;
// 建立一個 thread，執行 people_count 任務，傳入 arg[i] 指標參數
int temp;
if((temp=pthread_create(&callThd[i], &attr, people_count, (void *)&arg[i]))!= 0) {
printf("can't create thread: %d\n",temp);
return 1;
}
}
// 回收性質設定
pthread_attr_destroy(&attr);
void *status;
for (int i = 0; i < NUMTHRDS; i++) {
// 等待每一個 thread 執行完畢
pthread_join(callThd[i], &status);
}
// 所有 thread 執行完畢，印出結果
printf("result = %d \n", *pos);
// 回收互斥鎖
pthread_mutex_destroy(&mutexsum);
// 離開
pthread_exit(NULL);
return 0;
}
```
> 參考 [Lock-free 程式設計議題](https://hackmd.io/@sysprog/lock-free-prog?fbclid=IwAR1qQIDI5er1cbHQfgGw13B1ozIxDHTahxWF3aHC-bPu4kaxpfv360Q-zaQ)
:::warning
不要忽視 mutex lock 的成本，你嘗試平行化之前，應該適度切割 (partitioning) 數值範圍。
:notes: jserv
:::
Sign in with Wallet
Connect another wallet