# 重做 quiz3 測驗二的延伸問題，實作資料壓縮器並提出改進 x-compressor 的方案
contributed by < `blueskyson` >
###### tags: `linux2020`
## [quiz3](https://hackmd.io/@sysprog/2020-quiz3) 測驗二的延伸問題
在 [C 語言：數值系統篇](https://hackmd.io/@sysprog/c-numerics) 中，我們介紹過 power of 2 (漢字可寫作「二的冪」)。
[GCC extension](https://gcc.gnu.org/onlinedocs/gcc/Other-Builtins.html) 其中兩個是 ctz 和 clz:
>Built-in Function: `int __builtin_ctz (unsigned int x)`
>- Returns the number of trailing 0-bits in x, starting at the least significant bit position. If x is 0, the result is undefined.
>Built-in Function: `int __builtin_clz (unsigned int x)`
>- Returns the number of leading 0-bits in x, starting at the most significant bit position. If x is 0, the result is undefined.
我們可用 `__builtin_ctz` 來實作 [LeetCode 342. Power of Four](https://leetcode.com/problems/power-of-four/) ，假設 int 為 32 位元寬度。實作程式碼如下:
```cpp=
bool isPowerOfFour(int num)
{
return num > 0 && (num & (num - 1))==0 &&
!(__builtin_ctz(num) OPQ);
}
```
### 1. 解釋上述程式運作原理
首先第一個判斷條件是 `num > 0` 不用多作解釋。
因為 2 的冪次有以下特性:
- 2 = `00010`, (2 - 1) = `00001`, 2 & (2 - 1) = `0`
- 4 = `00100`, (4 - 1) = `00011`, 4 & (4 - 1) = `0`
- 8 = `01000`, (8 - 1) = `00111`, 8 & (8 - 1) = `0`
- 16 = `10000`, (16 - 1) = `01111`, 16 & (16 - 1) = `0`
所以 `(num & (num - 1))==0` 可用來判定 2 的冪次
最後，我們需要從所有 2 的冪次中挑出 4 的冪次。觀察上面 2 的冪次的規律，可以發現只要是 4 的冪次，其 trailing 0-bits 的個數都為偶數，如: 4 的 trailing 0-bits 有 2 個、 16 的 trailing 0-bits 有 4 個。
如果 `num` 的 trailing 0-bits 為偶數，則 `!(__builtin_ctz(num) OPQ)` 回傳 `TRUE` ，很明顯的 ==OPQ== 為 `& 0x1` 。
### 2. 改寫上述程式碼，提供等價功能的不同實作，儘量降低 branch 的數量
首先使用 perf 分析本題解法的 branch 數量，我用以下程式碼做測試
```c=
#include <stdbool.h>
#include <stdlib.h>
bool isPowerOfFour(int num)
{
return num > 0 && (num & (num - 1))==0 && !(__builtin_ctz(num) & 0x1);
}
int main(int argc, char *argv[]) {
int n = atoi(argv[1]);
isPowerOfFour(n);
return 0;
}
```
編譯，並將執行檔命名為 `pow4` 後使用以下腳本計算 branch 數量:
```non
$ cat test.sh
for ((i = 0; i < 32; i++))
do
((num = 1<<i))
perf stat -r 10 -e branches -o "perf/$num" ./pow4 $num
grep branches "perf/$num" >> branch.txt
done
```
執行完這份 shell script ，會得到一份 branch.txt ，紀錄 $2^{0}$, $2^{1}$,..., $2^{31}$ 所用的 branch 數量，可以觀察出原版程式所產生的 branch 數量介於 135000 至 138000 之間:
```
$ tr -d ' ,' < branch.txt
137565branches(+-0.53%)
137251branches(+-0.41%)
135862branches(+-0.53%)
136473branches(+-0.76%)
136823branches(+-0.50%)
136616branches(+-0.55%)
136150branches(+-0.56%)
137170branches(+-0.53%)
136848branches(+-0.34%)
137641branches(+-0.50%)
136450branches(+-0.56%)
135561branches(+-0.38%)
135748branches(+-0.40%)
136889branches(+-0.36%)
137397branches(+-0.45%)
138100branches(+-0.64%)
137186branches(+-0.50%)
135619branches(+-0.51%)
136723branches(+-0.46%)
137738branches(+-0.37%)
135568branches(+-0.61%)
136818branches(+-0.37%)
136434branches(+-0.51%)
135632branches(+-0.45%)
136480branches(+-0.56%)
137798branches(+-0.48%)
136101branches(+-0.35%)
136059branches(+-0.43%)
136315branches(+-0.57%)
136510branches(+-0.64%)
135684branches(+-0.42%)
```
接下來我寫了一個等價功能的實作: 首先檢查 `num` 中為 1 的 bit 的數量是否為 1 ，然後檢查 trailing zero 的數量是否為偶數，如果兩個條件都符合就回傳 `true`。
```c=
#include <stdbool.h>
#include <stdlib.h>
bool isPowerOfFour(int num)
{
return __builtin_popcount(num) == 1 &&
num & 0b01010101010101010101010101010101;
}
int main(int argc, char *argv[]) {
int n = atoi(argv[1]);
isPowerOfFour(n);
return 0;
}
```
一樣編譯成 `pow4` 後透過 shell script 執行，印出檔案，平均 branch 數量介於 135700 至 138000 之間，沒有把 branch 數量壓得更低:
```
136751branches(+-0.50%)
136739branches(+-0.52%)
137065branches(+-0.54%)
138113branches(+-0.45%)
136818branches(+-0.44%)
137006branches(+-0.38%)
136894branches(+-0.51%)
137206branches(+-0.51%)
135915branches(+-0.38%)
137785branches(+-0.53%)
137061branches(+-0.42%)
135992branches(+-0.49%)
136021branches(+-0.39%)
137311branches(+-0.53%)
136040branches(+-0.35%)
137076branches(+-0.47%)
137580branches(+-0.37%)
136060branches(+-0.39%)
136346branches(+-0.44%)
137884branches(+-0.53%)
137547branches(+-0.39%)
136204branches(+-0.40%)
136971branches(+-0.44%)
136752branches(+-0.55%)
137632branches(+-0.52%)
137220branches(+-0.51%)
137152branches(+-0.60%)
136830branches(+-0.47%)
136015branches(+-0.44%)
136820branches(+-0.27%)
135756branches(+-0.45%)
```
## 練習 [LeetCode 1009. Complement of Base 10 Integer](https://leetcode.com/problems/complement-of-base-10-integer/)
### 題目描述
Every non-negative integer `N` has a binary representation. For example, `5` can be represented as `"101"` in binary, `11` as `"1011"` in binary, and so on. Note that except for `N = 0`, there are no leading zeroes in any binary representation.
The complement of a binary representation is the number in binary you get when changing every `1` to a `0` and `0` to a `1`. For example, the complement of `"101"` in binary is `"010"` in binary.
For a given number `N` in base-10, return the complement of it's binary representation as a base-10 integer.
**Example 1:**
**Input:** `5`
**Output:** `2`
**Explanation:** `5` is `"101"` in binary, with complement `"010"` in binary, which is `2` in base-10.
**Example 2:**
**Input:** `7`
**Output:** `0`
**Explanation:** `7` is `"111"` in binary, with complement `"000"` in binary, which is `0` in base-10.
**Example 3:**
**Input:** `10`
**Output:** `5`
**Explanation:** `10` is `"1010"` in binary, with complement `"0101"` in binary, which is `5` in base-10.
### 解法
```c=
int bitwiseComplement(int N){
if (N == 0)
return 1;
unsigned int mask = 0xffffffff << (32 - __builtin_clz(N));
return ~(N | mask);
}
```
### 思路
以 `N = 5` 為例，其二進制表示為 `"101"` ，我想要求得 `5` 的補數，也就是 `"010"`
然而在 32 位元的 `int` 中， `5` 會表示成 `"00...00101"` ，直接對其取補數會得到 `"11...11010"` ，其二補數表示為 `-6` 不是我們想要的結果
所以我要做的事為: 將 N 的 leading zero 全部變成 `1` ，再進行取補數，也就是將 `"00...00101"` 變換為 `"11...11101"` 然後反轉為 `"010"` 即為正確答案
逐步解釋程式碼的意義:
```
Suppose N = 5 then:
1. s = 32 - __builtin_clz(N) = 3
2. mask = 0xffffffff << s = "11111111111111111111111111111000"
3. N | mask = "11111111111111111111111111111101"
4. ~(N | mask) = "00000000000000000000000000000010" #solution
```
### 提交
## 練習 [LeetCode 41. First Missing Positive](https://leetcode.com/problems/first-missing-positive/)
### 題目描述
Given an unsorted integer array `nums`, find the smallest missing positive integer.
**Follow up:** Could you implement an algorithm that runs in `O(n)` time and uses constant extra space?
**Example 1:**
**Input:** `nums = [1,2,0]`
**Output:** `3`
**Example 2:**
**Input:** `nums = [3,4,-1,1]`
**Output:** `2`
**Example 3:**
**Input:** `nums = [7,8,9,11,12]`
**Output:** `1`
### 解法
```c=
void sort(int arr[], int front, int end) {
if(front<end) {
int i, index=front;
for(i=front; i<end; i++) {
if(arr[i]<arr[end]) {
int tmp=arr[i]; //swap
arr[i]=arr[index]; //
arr[index]=tmp; //
index++;
}
}
int tmp=arr[i]; //swap pivot
arr[i]=arr[index]; //
arr[index]=tmp; //
sort(arr,front,index-1);
sort(arr,index+1,end);
}
}
int firstMissingPositive(int* nums, int numsSize){
sort(nums, 0, numsSize - 1);
int i = 0, min = 1;
while(i < numsSize) {
if (nums[i] == min)
min++;
else if (nums[i] > min)
break;
i++;
}
return min;
}
```
### 思路
先對 `nums` 做 quick sort ，然後將 `min` 設為 `1` ，從排序好的 `nums` 中依序檢查由 `1` 開始的連續正整數，每當檢查到一個連續正整數便將 `min` 加 1 ，否則跳出迴圈並回傳 `min`
### 提交
### 另一種解法
使用 bit array ，還沒實作
## 研讀 [2017 年修課學生報告](https://hackmd.io/@3xOSPTI6QMGdj6jgMMe08w/Bk-uxCYxz) ，理解 `clz` 的實作方式，並舉出 [Exponential Golomb coding](https://en.wikipedia.org/wiki/Exponential-Golomb_coding) 的案例說明
>[x-compressor](https://github.com/jserv/x-compressor) 是個以 [Golomb-Rice coding](https://en.wikipedia.org/wiki/Golomb_coding) 為基礎的資料壓縮器，實作中也用到 clz/ctz 指令，可參見 [Selecting the Golomb Parameter in Rice Coding](https://ipnpr.jpl.nasa.gov/progress_report/42-159/159E.pdf)。
### `clz` 效能測試
[2017 年修課學生報告](https://hackmd.io/@3xOSPTI6QMGdj6jgMMe08w/Bk-uxCYxz) 列出了 8 種 `clz` 的實做，分別為:
- binary search (iterative)
- binary search (bit mask)
- binary search (byte shfit)
- binary search (recursive)
- De Bruijn method
- Harley’s algorithm
- Harley’s algorithm inline assembly
- binary search (recursive) inline assembly
因為不知道[2017 年修課學生報告](https://hackmd.io/@3xOSPTI6QMGdj6jgMMe08w/Bk-uxCYxz)效能分析的原始碼，而且也不確定我的 CPU 的行為會不會符合他們效能分析的結果，所以我先用自己的方法測試一次，測試環境:
- CPU: Intel i7-10750H
- OS: Ubuntu 20.04 LTS
我使用 `rdtsc` 取得 clz 運算前後的 time stamp ，再將其相減取得 cycle 數，即下方程式碼的 `end - begin` 。測試時，將參數 `argv[1]` 作為 2 的冪次傳入並將其轉成 `int` 並存入 `shift` ，我針對每個 $2^{shift}$ 作為 `clz` 的 input 測試。
以 binary search (iterative) 為例，其測試程式碼如下:
```c=
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#pragma optimize( "", off )
int clz(uint32_t x) { // l := left, r := right
uint32_t n = 32, c = 16, l, r = x;
do {
l = r >> c;
if (l) { n -= c; r = l; }
c >>= 1;
} while (c);
return (n - r);
}
int main(int argc, char* argv[])
{
int shift = atoi(argv[1]);
/* input */
uint32_t x = 1 << shift;
uint64_t begin, end;
uint32_t begin_lo, begin_hi, end_lo, end_hi;
__asm__ __volatile__ ("rdtsc" : "=a" (begin_lo), "=d" (begin_hi));
/* clz implementation */
clz(x);
__asm__ __volatile__ ("rdtsc" : "=a" (end_lo), "=d" (end_hi));
begin = ((uint64_t)begin_hi << 32) | begin_lo;
end = ((uint64_t)end_hi << 32) | end_lo;
printf("%d %lu\n", shift, end - begin);
return 0;
}
```
將 binary search (iterative) 的測試程式碼命名為 `clz1.c` ，編譯為 `clz1` 並藉由以下 shell script 對 `clz1` 至 `clz8` 進行效能測試，測試方式為對 $2^0$ 至 $2^{31}$ 分別執行 100 次 `clz`:
```non
rm clz*.txt
for (( i = 0; i < 32; i++ ))
do
echo "testing for 1 << $i"
for (( j = 0; j < 100; j++))
do
./clz1 "$i" >> clz1.txt
./clz2 "$i" >> clz2.txt
./clz3 "$i" >> clz3.txt
./clz4 "$i" >> clz4.txt
./clz5 "$i" >> clz5.txt
./clz6 "$i" >> clz6.txt
./clz7 "$i" >> clz7.txt
./clz8 "$i" >> clz8.txt
done
done
```
再來用 gnuplot 畫出 2 的冪次對 cpu cycle 的 scatter ，以下是 gnuplot 指令:
:::spoiler
```non
set size 0.5, 0.5
set term png size 900, 900
set output "figure1.png"
set multiplot layout 2, 2 rowsfirst title "Binary Search CLZ" font ",18"
set xlabel "power of 2" font ",10"
set ylabel "rdtsc cycles" font ",10"
set xtics font ",10"
set ytics font ",10"
#
unset key
set title "iterative" font ",14"
plot [-1:32][:500] 'clz1.txt' with points pt 7 ps 1 lt rgb "red"
#
unset key
set title "bit mask" font ",14"
plot [-1:32][:500] 'clz2.txt' with points pt 7 ps 1 lt rgb "orange"
#
unset key
set title "byte shfit" font ",14"
plot [-1:32][:500] 'clz3.txt' with points pt 7 ps 1 lt rgb "blue"
#
unset key
set title "recursive" font ",14"
plot [-1:32][:500] 'clz4.txt' with points pt 7 ps 1 lt rgb "violet"
unset multiplot
```
:::
上圖 Binary Search 表示的結果大致和 [2017 年修課學生報告](https://hackmd.io/@3xOSPTI6QMGdj6jgMMe08w/Bk-uxCYxz) 測試的結果差不多，效能由優至劣排序依序是：
$$\text{bit mask} \approx \text{byte shift} > \text{iteration} > \text{recursive}$$
接下來再看另外 4 種 `clz` 的時間分佈圖:
:::spoiler
```non
set size 0.5, 0.5
set term png size 900, 900
set output "figure2.png"
set multiplot layout 2, 2 rowsfirst title "Other CLZ" font ",18"
set xlabel "power of 2" font ",10"
set ylabel "rdtsc cycles" font ",10"
set xtics font ",10"
set ytics font ",10"
#
unset key
set title "De Bruijn method" font ",14"
plot [-1:32][:500] 'clz5.txt' with points pt 7 ps 1 lt rgb "red"
#
unset key
set title "Harley’s algorithm" font ",14"
plot [-1:32][:500] 'clz6.txt' with points pt 7 ps 1 lt rgb "orange"
#
unset key
set title "Harley’s algorithm asm" font ",14"
plot [-1:32][:500] 'clz7.txt' with points pt 7 ps 1 lt rgb "blue"
#
unset key
set title "bnary search (recursive) asm" font ",14"
plot [-1:32][:500] 'clz8.txt' with points pt 7 ps 1 lt rgb "violet"
unset multiplot
```
:::
可以發現 De Bruijn method 是 8 種方法中速度最慢的，我推測其原因為函式內使用了一次乘法拖慢了速度，此外，將 Harley’s algorithm 和 binary search (recursive) 用組合語言實做的確有改善效能，不過具體改善多少需要比較他們的平均。
以上實驗是以 2 的冪次作為 input ，會有盲點，未來可能需要加一些隨機性較高的輸入值來觀察這些 `clz` 的實做是否會表現其他特性。
我順手也測試了 `__builtin_clz` 所耗的時間，發現速度比上述實做都要快得多
我將 `__builtin_clz` 包在函式 `clz` 中，並且用 gcc -S 編譯，查看其 x86 的實做:
```c
// c implementation
int clz(int x) {
return __builtin_clz(x);
}
```
```asm
; x86 assembly implementation
clz:
.LFB6:
.cfi_startproc
endbr64
pushq %rbp
.cfi_def_cfa_offset 16
.cfi_offset 6, -16
movq %rsp, %rbp
.cfi_def_cfa_register 6
movl %edi, -4(%rbp)
movl -4(%rbp), %eax ; load input value
bsrl %eax, %eax ; get the position of most significant set bit
xorl $31, %eax ; get complement of eax
popq %rbp
.cfi_def_cfa 7, 8
ret
.cfi_endproc
```
由上面 x86 實做可以看出和 `__builtin_clz` 有關的指令只有 3 行，所以會比 2017 報告中實做的快，至於 `bsrl` 是否有獨立的電路或是有被硬體優化則需要再探討。
最後我畫出 `clz` 消耗的 cycle 的 10% trimed mean 比較圖:
### Exponential Golomb coding
摘錄自維基百科，[Exponential-Golomb coding](https://en.wikipedia.org/wiki/Exponential-Golomb_coding) (或稱 Exp-Golomb code) 是一種 [Universal code](https://en.wikipedia.org/wiki/Universal_code_(data_compression))，也就是 Exp-Golomb code 能夠映射到所有正整數。
假設輸入為 `x` ，其編碼步驟為:
1. 以二進位寫下 `x + 1`
2. 數 `x + 1` 的二進位數字總共有多少位數，並在 `x + 1` 前面補上等同其位數減 1 的數量的 `0`
**Example 1:**
```
x = 5
1. Write down 5 + 1 = 6 as binary "110".
2. "110" has 3 digits, we write 2 "0"s preceding "110".
The Exp-Golomb code of 5 is "00110".
```
**Example 2:**
```
x = 24
1. Write down 24 + 1 = 25 as binary "11001".
2. "11001" has 5 digits, we write 4 "0"s preceding "11001".
The Exp-Golomb code of 24 is "000011001".
```
按照上述規則計算 0 ~ 8 的 Exp-Golomb code:
```
x step 1 step 2
0 ⇒ 1 ⇒ 1
1 ⇒ 10 ⇒ 010
2 ⇒ 11 ⇒ 011
3 ⇒ 100 ⇒ 00100
4 ⇒ 101 ⇒ 00101
5 ⇒ 110 ⇒ 00110
6 ⇒ 111 ⇒ 00111
7 ⇒ 1000 ⇒ 0001000
8 ⇒ 1001 ⇒ 0001001
```
如果要讓 Exp-Golomb code 編碼**有號**整數，只需要讓有號數一一對應無號數即能達到目的:
```
x map step 1 step 2
0 ⇒ 0 ⇒ 1 ⇒ 1
1 ⇒ 1 ⇒ 10 ⇒ 010
−1 ⇒ 2 ⇒ 11 ⇒ 011
2 ⇒ 3 ⇒ 100 ⇒ 00100
−2 ⇒ 4 ⇒ 101 ⇒ 00101
3 ⇒ 5 ⇒ 110 ⇒ 00110
−3 ⇒ 6 ⇒ 111 ⇒ 00111
4 ⇒ 7 ⇒ 1000 ⇒ 0001000
−4 ⇒ 8 ⇒ 1001 ⇒ 0001001
```
我們可以按照順序，每 $2^i$ 個 Exp-Golomb code 分為一組。對於每個編碼，前面有 $i$ 個 `0` 就代表該編碼屬於 $S_i$ ，每個編碼正中央紅色的 `1` 作為分界點，後綴的編碼以二進位表示為第 $S_i$ 組的第 $j$ 個元素，記為 $G_{i,j}$ :
$\ S_0 \quad \quad \quad S_1 \quad \quad \quad \quad \quad \quad \quad \quad \ \ S_2 \\
\{\color{red}1\},\
\{0\color{red}10,\ 0\color{red}11\},\
\{00\color{red}100,\ 00\color{red}101,\ 00\color{red}110,\ 00\color{red}111\},\ ... \\
G_{0,0} \ \ \ G_{1,0} \ \ \ G_{1,1} \quad \quad G_{2,0} \quad \ \ G_{2,1} \quad \ \ G_{2,2} \quad \ G_{2,3}$
### Order-k Exponential Golomb coding
接下來，對 Exp-Golomb code 的編碼方式一般化，前綴的 `0` 的數量一樣代表該編碼在第幾組，但是每一組的元素數量變為 $2^{i+k}$ 個，後綴的編碼也加長 $k$ 位，我們將這種編碼方式稱為 $k$ 階 Exp-Golomb code。以下為範例:
**order-0**
$\ S_0 \quad \quad \quad S_1 \quad \quad \quad \quad \quad \quad \quad \quad \ \ S_2 \\
\{\color{red}1\},\
\{0\color{red}10,\ 0\color{red}11\},\
\{00\color{red}100,\ 00\color{red}101,\ 00\color{red}110,\ 00\color{red}111\},\ ... \\
G_{0,0} \ \ \ G_{1,0} \ \ \ G_{1,1} \quad \quad G_{2,0} \quad \ \ G_{2,1} \quad \ \ G_{2,2} \quad \ G_{2,3}$
**order-1**
$\ \quad S_0 \quad \quad \quad \quad \quad \quad \quad S_1 \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad S_2 \\
\{\color{red}10,\ \color{red}11\},\
\{0\color{red}100,\ 0\color{red}101, \ 0\color{red}110,\ 0\color{red}111\},\
\{00\color{red}1000,\ 00\color{red}1001,\ 00\color{red}1010,\ ... \\
G_{0,0} \ \ G_{0,1} \quad \ G_{1,0} \quad G_{1,1} \quad G_{1,2} \quad G_{1,3} \quad \quad G_{3,0} \quad \quad \ G_{3,1} \quad \quad G_{3,2}$
**order-2**
$\ \quad S_0 \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad S_1 \\
\{\color{red}100,\ \color{red}101,\ \color{red}110,\ \color{red}111\},\
\{0\color{red}1000,\ 0\color{red}1001, \ 0\color{red}1010,\ 0\color{red}1011,\ 0\color{red}1100,\ 0\color{red}1101,\ 0\color{red}1110,\ 0\color{red}1111\},\ ...\\
G_{0,0} \ \ \ G_{0,1} \ \ G_{0,2} \ \ G_{0,3} \quad \ \ \ G_{1,0} \quad \ \ G_{1,1} \quad \ \ G_{1,2} \quad \ \ G_{1,3} \quad \ \ G_{1, 4} \quad \ \ G_{1, 5} \quad \ \ G_{1, 6} \quad \ G_{1, 7}$
產生 $k$ 階 Exp-Golomb code 的步驟為:
1. 以二進位寫下 $x+2^k$
2. 數 $x+2^k$ 的二進位數字總共有多少位數，並在 $x+2^k$ 前面補上等同其位數減 $(k+1)$ 的數量的 `0`
以下是 order-0 到 order-4 的前 13 個編碼:
|x |order-0 |order-1 |order-2 |order-3 |order-4 |
|--|-----------|-----------|-----------|-----------|-----------|
|0 |**1** |**1**0 |**1**00 |**1**000 |**1**0000 |
|1 |0**1**0 |**1**1 |**1**01 |**1**001 |**1**0001 |
|2 |0**1**1 |0**1**00 |**1**10 |**1**010 |**1**0010 |
|3 |00**1**00 |0**1**01 |**1**11 |**1**011 |**1**0011 |
|4 |00**1**01 |0**1**10 |0**1**000 |**1**100 |**1**0100 |
|5 |00**1**10 |0**1**11 |0**1**001 |**1**101 |**1**0101 |
|6 |00**1**11 |00**1**000 |0**1**010 |**1**110 |**1**0110 |
|7 |000**1**000|00**1**001 |0**1**011 |**1**111 |**1**0111 |
|8 |000**1**001|00**1**010 |0**1**100 |0**1**0000 |**1**1000 |
|9 |000**1**010|00**1**011 |0**1**101 |0**1**0001 |**1**1001 |
|10|000**1**011|00**1**100 |0**1**110 |0**1**0010 |**1**1010 |
|11|000**1**100|00**1**101 |0**1**111 |0**1**0011 |**1**1011 |
|12|000**1**101|00**1**110 |00**1**0000|0**1**0100 |**1**1100 |
### 實做 Exponential Golomb coding
**repo:** https://github.com/blueskyson/Exponential-Golomb-coding
**如何使用:**
首先將 repo 中的 `exp-golomb.c` 編譯成 `encode` ，再將 `decode` 連結至 `encode` 即完成編譯:
```
$ gcc exp-golomb.c -o encode
$ ln -s encode decode
```
`encode` 和 `decode` 的參數格式如下，以下將 `sample_text.txt` 轉成 order-4 exp-golomb code 的方法:
```
# Usage
encode [input file] [output file] [order-k]
decode [input file] [output file] [order-k]
# Example
./encode sample_text.txt text.encode 4
./decode text.encode text.decode 4
```
**主程式:**
這個小程式能將資料以 `uint8_t` 為一個單位編碼成 exp-golomb code ，在程式第 7 行，先設定 `MAX_ORDER` ，也就是程式能容許的 order-k 的限度，此程式最大的限度為 order-31 ，因為在轉換的過程是以 1 到 2 個 `uint32_t` 作為中繼的容器，如果 MAX_ORDER 大於等於 32 會造成編碼後的資料跨越 3 個 `uint32_t` ，我的沒有設計相應機制處理。
接下來設定 `BUFFER_SIZE` ，這個 macro 決定暫存陣列的長度。在 `encode` 函式中，暫存陣列為 `uint32_t buffer[BUFFER_SIZE]` ；在 `decode` 函式中，暫存陣列為 `uint8_t buffer[BUFFER_SIZE]`。
當 buffer 的元素達到 `WRITE_SIZE` 時便會將 buffer 的內容寫入 output file 並的清空 buffer 。
接下來我將 input file 以 file descriptor 的形式開啟，之後使用 `mmap` 讀檔、 output file 以 FILE* 的形式開啟，以 `fwrite` 寫入。
```c=
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdint.h>
#include <libgen.h> //basename
#include <fcntl.h> //open
#include <unistd.h> //close
#include <sys/mman.h> //mmap
#include <sys/stat.h>
#define MAX_ORDER 7 // must be smaller than 32
#define BUFFER_SIZE 100
#define WRITE_SIZE (BUFFER_SIZE - 2)
/* function status */
#define FAIL 0
#define SUCCESS 1
int encode(int, FILE*, int);
int decode(int, FILE*, int);
int main(int argc, char *argv[])
{
if (argc < 3) {
puts("Usage:\n"
"encode [input file] [output file] [order-k]\n"
"decode [input file] [output file] [order-k]");
return 0;
}
/* open argv[1] as input file */
int in_fd = open(argv[1], O_RDONLY);
if (in_fd < 0) {
puts("cannot open input file");
return 0;
}
/* open argv[2] as output file */
FILE *out_file = fopen(argv[2], "wb");
if (!out_file) {
puts("cannot open output file");
return 0;
}
/* use argv[3] as order */
long int order = 0;
if (argc >= 4) {
char* endptr;
order = strtol(argv[3], &endptr, 10);
if (endptr == argv[3]) {
puts("order is not a decimal number");
return 0;
}
if (order > MAX_ORDER) {
puts("order larger than max order");
return 0;
}
}
int status;
if (strcmp(basename(argv[0]), "encode") == 0) {
status = encode(in_fd, out_file, order);
} else {
status = decode(in_fd, out_file, order);
}
fclose(out_file);
close(in_fd);
if (status == FAIL) {
remove(argv[2]);
}
return 0;
}
```
**Encode**
以 `uint8_t*` 型態打開 input file 的 `mmap` ，轉換後的 exp-golomb code 會緊密排列在 uint32_t 的陣列中，如下圖示:
這裡所使用的編碼方法跟維基百科提到的一樣:
1. 以二進位寫下 $x+2^k$
2. 數 $x+2^k$ 的二進位數字總共有多少位數，並在 $x+2^k$ 前面補上等同其位數減 $(k+1)$ 的數量的 `0`
當轉換完所有資料後，在檔案尾加上一個 (uint32_t) 1 作為結束檔案的記號。
```c=73
int encode(int in_fd, FILE* out_file, int order)
{
struct stat s;
if (fstat(in_fd, &s) < 0) {
puts("cannot get status of iniput file");
return FAIL;
}
uint8_t *map = (uint8_t*) mmap(0, s.st_size, PROT_READ, MAP_PRIVATE, in_fd, 0);
if (map == MAP_FAILED) {
puts("cannot open mmap");
return FAIL;
}
uint32_t buffer[BUFFER_SIZE];
memset(buffer, 0, sizeof(buffer));
int buffer_index = 0;
int cursor = 32;
uint32_t offset = 1 << order;
/* start to encode */
for (long int i = 0; i < s.st_size; i++) {
uint32_t current = (uint32_t) map[i];
current += offset;
int clz = __builtin_clz(current);
int unary_width = 31 - clz - order;
int binary_width = 32 - clz;
/* write unary */
cursor -= unary_width;
if (cursor <= 0) {
cursor += 32;
buffer_index++;
}
/* write binary */
if (cursor < binary_width) { // truncate in uint32_t
buffer[buffer_index++] |= current >> (binary_width - cursor);
buffer[buffer_index] |= current << (32 - (binary_width - cursor));
cursor = cursor + 32 - binary_width;
} else {
buffer[buffer_index] |= current << (cursor - binary_width);
cursor -= binary_width;
}
/* write buffer */
if (buffer_index >= WRITE_SIZE) {
fwrite(buffer, 4, WRITE_SIZE, out_file);
uint32_t tail = buffer[buffer_index];
memset(buffer, 0, sizeof(buffer));
buffer[0] = tail;
buffer_index = 0;
}
}
/* finalize */
buffer[buffer_index + 1] = (uint32_t) 1; //end signal
fwrite(buffer, 4, buffer_index + 2, out_file);
return SUCCESS;
}
```
**Decode**
這裡只要將 encode 的步驟反著做就好了，以 `uint32_t*` 型態打開 input file 的 `mmap` ，轉換後的資料會轉型成 `uint8_t` 。唯一需要注意的是，讓 `uint32_t` 左移或右移 32 bit 時，會發生 [-Wshift-count-overflow] ，所以在 191 行寫了一個判斷去規避平移 32 bit 。
```c=134
int decode(int in_fd, FILE* out_file, int order)
{
struct stat s;
if (fstat(in_fd, &s) < 0) {
puts("cannot get status of iniput file");
return FAIL;
}
uint32_t *map = (uint32_t*) mmap(0, s.st_size, PROT_READ, MAP_PRIVATE, in_fd, 0);
if (map == MAP_FAILED) {
puts("cannot open mmap");
return FAIL;
}
uint8_t buffer[BUFFER_SIZE];
memset(buffer, 0, sizeof(buffer));
int buffer_index = 0;
int cursor = 32;
uint32_t offset = 1 << order;
int map_index = 0;
uint32_t current = map[0];
while (1) {
/* read unary */
int unary_width = 0;
if (current == 0) { // truncate in unary field
current = map[++map_index];
int clz = __builtin_clz(current);
unary_width = cursor + clz;
cursor = 32 - clz;
} else {
int clz = __builtin_clz(current);
unary_width = clz - (32 - cursor);
cursor = 32 - clz;
}
/* end of file is (uint32_t) 1 ,
* the leading zero of end of file is 31 */
if (unary_width >= 31) {
break;
}
/* read binary */
int binary_width = unary_width + 1 + order;
uint32_t tmp = 0;
if (binary_width > cursor) { // truncate in binary field
tmp = current << (binary_width - cursor);
current = map[++map_index];
binary_width -= cursor;
cursor = 32;
}
tmp |= current >> (cursor - binary_width);
tmp -= offset;
buffer[buffer_index++] = (uint8_t) tmp;
/* be careful for left shift 32 bits */
cursor -= binary_width;
if (cursor == 0) {
current = 0;
} else {
int shift = 32 - cursor;
current = (current << shift) >> shift;
}
if (buffer_index == BUFFER_SIZE) {
fwrite(buffer, 1, BUFFER_SIZE, out_file);
buffer_index = 0;
}
}
/* finalize */
if (buffer_index != 0) {
fwrite(buffer, 1, buffer_index, out_file);
}
return SUCCESS;
}
```
### 使用 valgrind 檢查程式
因為整個程式完全沒使用 malloc ，所以很難發生 memory leak ，使用 `valgrind -q --leak-check=full` 編解碼都沒有被偵測到 memory leak。 下面為 encode 和 decode 的記憶體使用狀況:
### 使用 perf 觀察
**Encode**
```
# perf stat --repeat 5 -e cache-misses,cache-references,instructions,cycles,context-switches,branches,branch-misses ./encode 74-0.txt 74-0.encode
Performance counter stats for './encode 74-0.txt 74-0.encode' (5 runs):
5,4589 cache-misses # 32.056 % of all cache refs ( +- 20.69% ) (45.64%)
17,0295 cache-references ( +- 4.66% ) (97.83%)
2515,6815 instructions # 1.46 insn per cycle ( +- 0.56% )
1718,9863 cycles ( +- 0.80% )
1 context-switches ( +- 40.82% )
267,6514 branches ( +- 3.97% ) (54.36%)
<not counted> branch-misses ( +-100.00% ) (2.17%)
0.006105 +- 0.000246 seconds time elapsed ( +- 4.02% )
```
**Decode**
```
# perf stat --repeat 5 -e cache-misses,cache-references,instructions,cycles,context-switches,branches,branch-misses ./decode 74-0.encode 74-0.decode
Performance counter stats for './decode 74-0.encode 74-0.decode' (5 runs):
7,6270 cache-misses # 32.452 % of all cache refs ( +- 20.52% ) (34.36%)
23,5021 cache-references ( +- 12.40% ) (77.71%)
2955,5866 instructions # 1.32 insn per cycle ( +- 2.98% ) (96.42%)
2245,0086 cycles ( +- 5.05% )
1 context-switches ( +- 25.00% )
298,9814 branches ( +- 2.96% ) (69.22%)
<not counted> branch-misses ( +- 28.11% ) (22.29%)
0.00845 +- 0.00136 seconds time elapsed ( +- 16.08% )
```
為了讓程式執行的夠久，讓 perf 能有效偵測效能，我使用 repo 中的 74-0.txt 進行測試，但是仍然無法測出 branch-mises 不知道是系統權限沒設定好或是要讓程式執行更久才能統計出來。
我對 order-0 的 encode 和 decode 測試 5 次，得到上面兩份數據，很明顯發現我的 cache-misses 非常高，平均高達 32% ，代表我程式寫法對變數的 locality 處理非常差，我認為這是目前值得我去改進的部份。
我初步猜測可能原因是類似 `dict` 作業 copy 和 reference 機制的問題，也就是我常常存取的陣列沒有配置在較近的 heap 中。
### Exponential Golomb coding 效能分析
下圖是 8-bit 數值經過 order 0 到 4 的編碼後所佔據的位元長度，可以看到除了 order-4 之外，其他 exp-golomb code 都在原數值大於 50 之前達到 8 bit 的長度。 ASCII code 中常用的數字和英文字母介於 48 的 `'0'` 到 122 的 `'z'` 之間，幾乎都是會讓 exp-golomb code 超過 8 bit 的區間，若沒有修改 exp-golomb 的機制，根本無法壓縮。
x-compressor 以不斷更新模型，盡量讓出現頻率最高的數字使用位元長度較少的 golomb-rice code 達成壓縮的效果
## 閱讀論文 [Selecting the Golomb Parameter in Rice Coding](https://ipnpr.jpl.nasa.gov/progress_report/42-159/159E.pdf)
這篇論文的目標是針對一個資源，找出 golomb-power-of-2 (GPO2) 最佳的參數 $k$ ，使得此資源的 GPO2 的 bit rate 期望值 (編碼的平均長度) 最小。找出參數 $k$ 的過程必須快、節省計算資源。
### Analysis
**A. Minimizing Code Rate**
首先，我們把**非負整數**設定為目標資源，用 GPO2 編碼非負整數的隨機變數 $δ$ 時，前綴的一進位值需要 $\left\lfloor\dfrac{j}{2^k}\right\rfloor + 1$ 個 bit ，後綴二進值則固定為 $k$ 個 bit ，此時的編碼率 (對所有非負整數編碼的長度) 為:
$$ R_k = \sum_{j=0}^\infty(\left\lfloor\dfrac{j}{2^k}\right\rfloor+1+k) \quad \text{Prob}[δ=j]$$
觀察上面的函數，在固定 $j$ 的情況下 $\left\lfloor\dfrac{j}{2^k}\right\rfloor+1+k$ 為一凸函數，又 $R_k$ 為這些凸函數的總和，所以 $R_k$ 也是凸函數，因此 $R_k$ 的區域最小值亦為全域最小值。任何非負整數 $k^*$ 滿足
$$R_{k^*-1} \leq R_{k^*} \leq R_{k^*+1} \quad \quad (1)$$
即為 GPO2 編碼長度最小的解。為了方便計算，將 $R_{-1}$ 視為 $\infty$ 。
接下來，我們考慮已知 $µ=E[δ]$ 的情況下，如何依據先前的編碼估算參數 $k$ 的值。給定 $µ$ ， GPO2 的最佳參數 $k$ 一定會在區間 $[k^*_{min}(µ), \ k^*_{max}(µ)]$ ，其中
$$\begin{align*}
&k^*_{min}(µ) = \text{max}\left\lbrace 0,\left\lfloor \text{log}_2\left(\dfrac{2}{3}(µ+1)\right) \right\rfloor \right\rbrace \\
&k^*_{max}(µ) = \text{max}\lbrace 0, \lceil \text{log}_2 µ \rceil\rbrace
\end{align*} \quad \quad (2)$$
使得 $R_{k^*}$ 接近最低點 (編碼率最佳) ，下圖顯示了 $k^*_{min}(µ)$ 和 $k^*_{max}(µ)$ 的成長趨勢
再來由 (2) 可以推得 $k^*_{max}(µ)-k^*_{min}(µ) \leq 2 \ \text{for all} \ µ$ ，其中的 $µ$ 是樣本的平均值。以 $µ$ 計算編碼率 $R_k$ 和 $R_{k-1}$ 的差值
$$\begin{align*}
R_k - R_{k-1} &= 1 + \sum_j(\left\lfloor \dfrac{j}{2^k} \right\rfloor - \left\lfloor \dfrac{j}{2^{k-1}} \right\rfloor) \\
&= 1 + \sum_j \left\lfloor \dfrac{1}{2} - \dfrac{j+1}{2^k} \right\rfloor \quad \text{Prob}[δ=j]
\end{align*}$$
其中
$$-\dfrac{1}{2}-\dfrac{j}{2^k} ≤ \left\lfloor \dfrac{1}{2} - \dfrac{j+1}{2^k} \right\rfloor ≤ \dfrac{1}{2}-\dfrac{j+1}{2^k}$$
所以我們得到
$$\dfrac{1}{2} - \dfrac{µ}{2^k} ≤ R_k - R_{k-1} ≤ \dfrac{3}{2} - \dfrac{µ+1}{2^k} \quad \quad (3)$$
我們可以將編碼率表示為
$$\begin{align*}
R_k &= k+1+E\left[ \left\lfloor\dfrac{δ}{2^k}\right\rfloor \right]\\
&= k+1+\dfrac{1}{2^k}\left(E[δ]-E\left[ δ-2^k \left\lfloor\dfrac{δ}{2^k} \right\rfloor \right]\right)\\
&=k+1+\dfrac{1}{2^k}(µ-\bar{r}_k) \quad \quad (4)
\end{align*}$$
其中 $\bar{r}_k$ 的值等於
$$\bar{r}_k=δ-2^k\left\lfloor \dfrac{δ}{2^k} \right\rfloor \quad \quad (5)$$
其值剛好等於 $δ \ \text{mod} \ 2^k$ 。由 (4) 可以進一步推出 $R_k^{\text{up}}$ 和 $R_k^{\text{low}}$ 代表 $R_k$ 的上下界
$$\begin{align*}
&R_k ≤ R_k^{\text{up}}(µ) \triangleq k+1+\dfrac{µ}{2^k}\\
&R_k ≥ R_k^{\text{low}}(µ) \triangleq k+\text{max} \left\lbrace 1, \dfrac{µ+1}{2^k}\right\rbrace
\end{align*}$$
我們由 (4) 推導出的算式進一步推論在滿足 (1) 的 $R_k ≤ R_{k-1}$ 時， $µ$ 的範圍為 (6)
$$\begin{align*}
&k+1+\dfrac{1}{2^k}(µ-\bar{r}_k) ≤ (k-1)+1+\dfrac{1}{2^{k-1}}(µ-\bar{r}_{k-1})\\
\Rightarrow \quad &µ≥2^k+2\bar{r}_{k-1}-\bar{r}_{k} \quad \quad (6)
\end{align*}$$
最後，針對 (6) 我們討論邊界問題: 因為 $\bar{r}_0 = 0$ ，當 $µ < 2^1+2\bar{r}_0-\bar{r}_1$ (即無法滿足 (6)) 時，就將參數 $k$ 設為 0 ；否則，選擇滿足 (6) 的最大的 $k$ 值。如此一來，便可以得到最佳的 $k$ 值，我們可以預先計算每個 $µ$ 在 (6) 中應該要對應哪個 $k$ 值，並做一個表格記錄下來。
**B. Dynamic Range Limit and Uncoded Data**
有兩個因素影響了 GPO2 在現實的應用:
1. 儘管 Golomb code 允許任意大小的輸入，但在實際應用中，樣本通常有最大值，縮小了需要考慮的編碼範圍。
2. 許多 entropy 編碼器可以選擇不要編碼樣本，參數選擇 ( $k$ ) 的過程中應該要包含不予編碼的狀況。
再來，當樣本有 $N$ bit 時， $k ≥ N − 1$ 會使得編碼率大於等於原始樣本，這樣就沒有壓縮的效果，所以我們預先施加約束: $k ≤ N − 2$
### Coding a Geometric Source
**A. Optimum Code Selection**
將 $δ$ 以幾何分布建模如下，評估 $\bar{r}_k$
$$\text{Prob}[δ=j]=(1-α)α^j \quad \quad α ∈ (0,1)$$
## 改進 x-compressor 的方案
### 解析 x-compressor 原始碼
**compress**
首先 input file 藉由 `stdin` 進入程序，在 `x.c` 的 main function 中，宣告 `FILE* istream` 指向 `stdin` 的資料流；`FILE* ostream` 指向 `stdout` 的資料流。
接下來，將資料分成兩個 `layer` ，將 `stdin` 的資料載入 heap 中，並用 `layer[0].data` 指向它，之後 `layer[1].data` 會被分配到一塊夠大的空間儲存轉換後的資料 。透過 `x_init()` 初始化 `table` ，也就是字元分佈的模型，首先按照每個字元的大小順序排好。
然後進入到 `x_compress` 函式，使用 `initiate` 初始化 `io` ， `io` 負責將轉換好的編碼緊密的排在 `layer[1].data` 中， `io.ptr` 指向 `layer[1].data` 中即將被寫入的位址，功能類似我前面的程式的 `buffer_index` 、 `io.c` 指向即將被寫入的 bit ，功能類似我前面的程式的 `cursor`
**todo:**
Sign in with Wallet
Connect another wallet