# 2020q3 Homework (quiz2)
contributed by < `hankluo6` >
### `測驗 1`
```cpp
#include <stdbool.h>
#include <stddef.h>
#include <stdint.h>
#include <string.h>
bool is_ascii(const char str[], size_t size)
{
if (size == 0)
return false;
int i = 0;
while ((i + 8) <= size) {
uint64_t payload;
memcpy(&payload, str + i, 8);
if (payload & MMM)
return false;
i += 8;
}
while (i < size) {
if (str[i] & 0x80)
return false;
i++;
}
return true;
}
```
#### `MMM`
- [x] `(d)` `0x8080808080808080`
#### 說明
`(i + 8) <= size` 用來檢查現在字元數量 `i` 有沒有超過要檢測的字串長度 `size`，並將要檢測的記憶體位置 `str + i` 後面 8 個 bytes 的值複製至 `payload` 中，下方回傳 `false` 可得知 `payload` 用來檢驗是否為 128 內的 ASCII 碼。
從 `uint64_t` 可知 `payload` 為 8 bytes 位元之變數，每個 bytes 代表一個字元 `char`，而每個字元要檢查其值是否在 128 以內，以 16 進為表示就是為 0x00 ~ 0x79。如果該字元與 0x80 做 AND 運算不為 0 的話，表示該字元之 ASCII 碼為 128 以上。
可將上述字元操作擴充成 8 個字元同時檢測，故答案為 `(d)`。
#### 延伸問題
##### 1. 解釋上述程式的運作原理，特別是為何用到 memcpy 呢？提示: 與 data alignment 有關，詳閱 [C 語言：記憶體管理、對齊及硬體特性](https://hackmd.io/@sysprog/c-memory)
* 使用 `memcpy` 可以讓我們一次操作 8 bytes 的資料，如不使用 `memcpy` 而是直接將字串元素與 `(d)` 做 AND 運算，如 `*(str + i) & 0x8080808080808080`，會被解釋為該字元與 `(d)` 的 AND 運算。
* 使用強制轉型可解決上述問題，將 `payload` 改寫成`uint64 *payload = (uint64_t *) (str + i);` 並做 AND 運算，也許你一樣會得到正確的結果，但其實這種結果並非每次都正確。
* 此種記憶體轉型有可能會造成 [Unaligned data access](https://www.kernel.org/doc/Documentation/unaligned-memory-access.txt)，因編譯器不能保證轉型完後，新的 data type 符合它自身的 alignment 需求，在不同的 Arm 處理器下，unaligned data access 造成的行為[不盡相同]((https://heyrick.eu/armwiki/Unaligned_data_access))，但在 intel 處理器下通常能透過多次抓取資料來處理 unaligned data，所以上述轉型才有可能會正確。
* 可以透過編譯指令來提示是否有 unalignmeent data access 的情況發生
1. gcc 8 以上的版本提供編譯選項 `-Wcast-align=strict` 來檢驗是否有 unaligned memory access 的情況發生
```
gcc -Wcast-align=strict align.c
```
3. 或是透過交叉編譯器 gcc-arm-linux-gnueabi 將編譯目標設定在 Arm 處理器上：
```
arm-linux-gnueabi-gcc -Wcast-align align.c
```
上述兩種方法在改寫為 `uint64 *payload = (uint64_t *) (str + i);` 的程式上均會發出下列警告
```
align.c:14:29: warning: cast increases required alignment of target type [-Wcast-align]
uint64_t* payload = (uint64_t *)(str+i);
```
所以此種轉型方法並不通用。
* 回到 `memcpy` 上，使用 `memcpy` 為甚麼不會有 unaligned data access 的問題？參考 glibc 裡 [memcpy.c](https://github.com/lattera/glibc/blob/master/string/memcpy.c) 的實作， `memcpy` 內部會先將 unalignment 的部分使用 byte copy 來對齊，才使用 word copy 來複製接下來的資料，而當剩餘資料長度小於 word 時，再使用 byte copy 複製剩餘資料。
##### 2. 將上述技巧應用於「給予一個已知長度的字串，檢測裡頭是否包含有效的英文大小寫字母」
詳細可參考測驗 5，概念為當字元在 `A-Z` 或 `a-z` 時，將第該字元的第 8 個 bit 設為 1，用該 bit 判斷是否為英文大小寫字母。
```c
bool is_alpha(char *s, size_t n)
{
size_t k = n / 8;
for (size_t i = 0; i < k; i++, s += 8) {
uint64_t *chunk = (uint64_t *) s;
uint64_t A = *chunk + PACKED_BYTE(128 - 'A');
uint64_t Z = *chunk + PACKED_BYTE(127 - 'Z');
uint64_t a = *chunk + PACKED_BYTE(128 - 'a');
uint64_t z = *chunk + PACKED_BYTE(127 - 'z');
if (!(((A ^ Z) & PACKED_BYTE(0x80) | (a ^ z) & PACKED_BYTE(0x80)) & PACKED_BYTE(0x80)))
return false;
}
k = n % 8;
while (k--) {
char c = *(s++);
if (!(('A' <= c && c <= 'Z') || ('a' <= c && c <= 'z')))
return false;
}
return true;
}
```
##### 3. 承 (2)，考慮常見的英文標點符號，判斷輸入字串是否為有效的字元集，避免逐一字元比對
先來定義何謂 **"常見的英文標點符號"**，這邊我使用 [7esl.com](https://7esl.com/punctuation-marks/) 提供的 14 個常用的標點符號。
接著把標點符號與對應的 ASCII Code 列出來觀察
| Marks | Decimal | Hexadecimal |
| ------- | ------- | ------------------- |
| (space) | 32 | 0010 0000 |
| ! | 33 | 0010 0001 |
| " | 34 | 0010 0010 |
| ' | 39 | 0010 0111 |
| () | 40,41 | 0010 1000,0010 1001 |
| , | 44 | 0010 1100 |
| - | 45 | 0010 1101 |
| . | 46 | 0010 1110 |
| / | 47 | 0010 1111 |
| : | 58 | 0011 1010 |
| ; | 59 | 0011 1011 |
| ? | 63 | 0011 1111 |
| \ | 92 | 0101 1010 |
| [] | 91,93 | 0101 1011,0101 1101 |
為了計算方便，我將 `\`也加入至符號中。
:::info
…(Ellipsis) 不在 Ascii Code 裡，故不考慮
:::
接著來看一個 SSE 常用到的 type：`__m128i`，其定義為
```c
typedef union __declspec(intrin_type) __declspec(align(16)) __m128i {
__int8 m128i_i8[16];
__int16 m128i_i16[8];
__int32 m128i_i32[4];
__int64 m128i_i64[2];
unsigned __int8 m128i_u8[16];
unsigned __int16 m128i_u16[8];
unsigned __int32 m128i_u32[4];
unsigned __int64 m128i_u64[2];
} __m128i;
```
可以看到為一個裝有 128 bits 的 union，因為是 union，所以內部記憶體空間是共享的，總共佔有 128 bits。
接著來看使用到的指令，以下指令皆可從 [Intel Intrinsics Guide](https://software.intel.com/sites/landingpage/IntrinsicsGuide/) 中找到。
1. `_mm_setr_epi8` 將參數裡的字元以相反的順序放入到目標變數(`__m128i`)中
e.g. `_mm_setr_epi8(0x10, 0x11, 0x12, ..., 0x24, 0x25);` 會回傳一個 type 為 `__m128i` 且內容為 `0x2524232221...121110` 的數字
2. `_mm_cmpestrm` 會依照所選擇的 mode 來決定如何比對兩個不同的字串，如果符合就在該字元的位置存放 `0xFF`，而 mode 選擇以下：
* `_SIDD_UBYTE_OPS` 表示字元為 8 bits
* `_SIDD_CMP_RANGES` 表示比較 [range](https://doc.rust-lang.org/core/arch/x86_64/constant._SIDD_CMP_RANGES.html)，檢測字元是否在給定的範圍內
* `_SIDD_UNIT_MASK` 回傳 bytes/word 的 mask
3. `_mm_test_all_ones` 如果數字內每個 bit 都為 1 則回傳 1，否則回傳 0
主程式的邏輯很簡單，首先使用 `_mm_setr_epi8` 將合理字元的範圍 pack 成 16 bytes 的變數。接著，每次對 16 bytes 的字串進行 `_mm_cmpestrm`，如果符合定義的範圍(`valid_bits`)，則該函式會回傳 16 bytes 的 1，並透過 `_mm_test_all_ones` 來檢查。
如果剩餘字串少於 16 bytes，則需要另外判斷。我透過 `check` 函式將 16 bytes 的值拆分成兩個 8 bytes，並進行 bitwise 的檢查，當出現不為 1 時表示出現不在範圍內的字元。
:::warning
這邊原本是想直接用函式來檢查未滿 16 bytes 的值是否全為 1，但翻閱 Guide 找不到相關的函式，也不能直接對 `__m128i` 進行 bit 操作，故只能寫成冗長的程式碼，如果有其他想法歡迎提供。
:::
```c
bool check(__m128i var, size_t last_size)
{
uint64_t val[2];
memcpy(val, &var, sizeof(val));
last_size *= 8; //bytes to bits
if (last_size > 64) {
if (val[0] != 0xFFFFFFFFFFFFFFFF)
return false;
last_size -= 64;
while (last_size) {
if (!((var[1] >> (last_size - 1)) & 1))
return false;
--last_size;
}
} else {
while (last_size) {
if (!((var[0] >> (last_size - 1)) & 1))
return false;
--last_size;
}
}
return true;
}
bool is_valid(const char str[], size_t size)
{
__m128i valid_bits = _mm_setr_epi8('a', 'z', 'A', ']', '?', '?', ',', '/', ' ', '"', ':', ';', '\'', ')', 0, 0);
const char mode = _SIDD_UBYTE_OPS | _SIDD_CMP_RANGES | _SIDD_UNIT_MASK;
unsigned int i = 0;
while (i != size) {
size_t last_size = (i + 16) <= size ? 16 : (size - i);
__uint128_t payload;
memcpy(&payload, str + i, last_size);
__m128i *ptr = (__m128i*)&payload;
const __m128i mask = _mm_cmpestrm(valid_bits, 14, *ptr, last_size, mode);
if (last_size != 16 && !check(mask, last_size)) {
return false;
} else if (!_mm_test_all_ones(mask))
return false;
i += last_size;
}
return true;
}
```
:::info
踩坑紀錄：
使用 `_mm_cmpestrm` 時輸入將 `const char mode = _SIDD_UBYTE_OPS | _SIDD_CMP_RANGES | _SIDD_UNIT_MASK;` 放入第五個參數，卻一直跳出 `the fifth argument must be an 8-bit immediate`。因為 gcc 預設不會將在編譯時期決定 function 內的 const 的值，所以才一直報錯。
解決方法為使用任一編譯器優化選項 `-O1~O3` 來讓編譯器在編譯時期就決定 `mode` 的值。
參考 stack overflow：
1. [_mm256_slli_si256: error “last argument must be an 8-bit intermediate](https://stackoverflow.com/questions/31315491/mm256-slli-si256-error-last-argument-must-be-an-8-bit-intermediate)
2. [Is variable defined by const int determined at compilation time?](https://stackoverflow.com/questions/27048025/is-variable-defined-by-const-int-determined-at-compilation-time)
:::
### `測驗 2`
```c
uint8_t hexchar2val(uint8_t in)
{
const uint8_t letter = in & MASK;
const uint8_t shift = (letter >> AAA) | (letter >> BBB);
return (in + shift) & 0xf;
}
```
#### `MASK`
- [x] `(c)` `0x40`
#### `AAA`
- [x] `(a)` `3`
#### `BBB`
- [x] `(b)` `6`
#### 說明
先將各字元的編碼寫下來
| char | 7 | 6 | 5 | 4 | 3 | 2 | 1 | 0 | hex |
| ----- | --- |:---:| --- | --- | --- | --- | --- |:---:| --- |
| `'0'` | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0x30 |
| `'9'` | 0 | 0 | 1 | 1 | 1 | 0 | 0 | 1 | 0x39 |
| `'a'` | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 1 | 0x61 |
| `'f'` | 0 | 1 | 0 | 0 | 0 | 1 | 1 | 0 | 0x66 |
| `'A'` | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 1 | 0x41 |
| `'F'` | 0 | 1 | 1 | 0 | 0 | 1 | 1 | 0 |0x46
可以看到數字 `0` ~ `9` 與英文字 `a` ~ `f` 或 `A` ~ `F` 的差別在第 6 個 bit 是否為 1，而 `a` ~ `f` 與 `A` ~ `F` 的差別則是在第 5 個 bit 是否為 1。
從變數名稱 `letter` 可以猜測應為判斷是否為字母，判斷字母的 bit 為第 6 個 bit，所以 `MASK` 應為 0x40；再來因為 `return` 的值會與 0xf 做 AND，可以只觀察後 4 個 bits，數字的後 4 個 bits 剛好就是對應到的該數字的值，字母的後 4 個 bits 不管是大寫還是小寫皆相同，所以可以一起做判斷。
字母的後 4 個 bits 為該字母所代表的數字之個位數字，亦即：`'a'` 表示 `10`，則後 4 個 bits 為 `1`，`'d'` 表示 `13`，則後 4 個 bits 為 3。所以我們只需將後 4 個 bits 加上數字 10 就行了，用第 6 個 bit 的位移來組合出數字 10，故答案為 `(a)` 與 `(b)`。
#### 延伸問題
##### 1. 將 hexchar2val 擴充為允許 "0xFF", "0xCAFEBABE", "0x8BADF00D" 等字串輸入，且輸出對應的十進位數值 Hexspeak
```c
/* input must start with 0x and the string length need to less than 32 */
uint64_t hex2val(char *s) {
s += 2;
uint64_t val = 0;
while (*s != '\0') {
if (*s == 'x' || *s == 'X')
break;
val <<= 4;
val |= hexchar2val((uint8_t)(*(s++)));
}
return val;
}
```
程式邏輯很簡單，開頭 `s += 2` 用來移除 hex 中 `0x` 的部分，之後從字串頭開始將每個字元轉換為對應的數值。向左位移 4 bits 是因為每個字元在 16 進制表示中只佔用 4 個 bits。所以每次 `*s` 往低位元移動時，轉換的值皆要往高位元移動 4 bits 來維持數字的正確性，而 OR 運算則是將新轉換的字元之值加入到現有的 `val` 當中。
### `測驗 3`
```c
const uint32_t D = 3;
#define M ((uint64_t)(UINT64_C(XXX) / (D) + 1))
/* compute (n mod d) given precomputed M */
uint32_t quickmod(uint32_t n)
{ uint64_t quotient = ((__uint128_t) M * n) >> 64;
return n - quotient * D;
}
```
#### `XXX`
- [x] `(f)` `0xFFFFFFFFFFFFFFFF`
#### 說明
取餘數的公式可以寫成 $n\ mod\ d = n - \lfloor{\frac{n}{d}}\rfloor \times d$，而 $\lfloor\frac{n}{d}\rfloor$ 根據測驗 3 的說明可以改寫為 $\lfloor n \times \frac{\frac{2^N}{d}}{2^N}\rfloor$。對照程式碼，右移 64 bytes 表示除以 $2^{64}$，所以 $N$ 應為 64。
接著我們嘗試將 $\lfloor{\frac{n}{d}}\rfloor$ 改寫為 $\lfloor \frac{n}{2^N} \times \lceil\frac{2^N}{d}\rceil\rfloor$，可以從後面推導回來：
$$
\begin{equation}
\begin{split}
\lfloor \frac{n}{2^N} \times \lceil\frac{2^N}{d}\rceil\rfloor &= \lfloor \frac{n}{2^N} \times \frac{2^N + r}{d}\rfloor \\
&= \lfloor \frac{n \times 2^N}{d \times 2^N} + \frac{n \times r}{d \times 2^N}\rfloor \\
&= \lfloor \frac{n}{d} + \frac{n \times r}{d \times 2^N}\rfloor \\
&= \lfloor q + \frac{k}{d} + \frac{n \times r}{d \times 2^N}\rfloor \\
&= q + \lfloor \frac{k}{d} + \frac{n \times r}{d \times 2^N} \rfloor\\
&= \lfloor\frac{n}{d}\rfloor \quad \text{if} \quad \frac{k}{d} + \frac{n \times r}{d \times 2^N} < 1
\end{split}
\end{equation}
$$
其中：
* $r = -2^N\ mod\ d$，用來將 $\frac{2^N}{d}$ 無條件進位到整數。
* $k = n\ mod\ d$，用來將 $\frac{n}{d}$ 無條件捨去到整數，並用 $q$ 表示。
* $q$ 為一整數，所以可以從 floor function 中提出。
再來就能進行完整的推導：
$$
\begin{equation}
\begin{split}
\lfloor\frac{n}{d}\rfloor &= \lfloor n \times \frac{\frac{2^N}{d}}{2^N}\rfloor \\
&= \lfloor \frac{n}{2^N} \times \lceil\frac{2^N}{d}\rceil\rfloor \quad \text{if} \quad \frac{k}{d} + \frac{n \times r}{d \times 2^N} < 1 \\
&= \lfloor \frac{n}{2^N} \times (\lfloor\frac{2^N - 1}{d}\rfloor+1)\rfloor \\
\end{split}
\end{equation}
$$
上述 ceil function 的轉換可參考 [wikipedia](https://en.wikipedia.org/wiki/Floor_and_ceiling_functions#Quotients) 的 Quotients 的部分。
與程式碼比對，$M$ 為 $\lfloor\frac{2^N - 1}{d}\rfloor + 1$，所以 `XXX` 為 $2^N - 1$，答案為 `(f)`，特別注意 `quotient` type 為 `uint64_t`，與 `UINT64_C`，其目的皆是在實現公式中 floor function 的部分。
我們對再來深入探討條件成立的情形，因 $n < 2^{32},N = 64$，所以：
$$
\begin{equation}
\begin{split}
\frac{k}{d}+\frac{n \times r}{d \times 2^N} &< \frac{k}{d} + \frac{r}{d} \times \frac{2^{32}}{2^{64}} \\
&=\frac{k}{d} + \frac{r}{d} \times \frac{1}{2^{32}} \\
&<\frac{d-1}{d} + \frac{d-1}{d} \times \frac{1}{2^{32}} \\
&= 1
\end{split}
\end{equation}
$$
上述式子要成立，在 $N = 64$ 且 $n = 2^{32}$ 的情況下，$d$ 必須小於 $2^{32} + 1$。
再來探討 $n$ 的範圍，從上面敘述已經得知 $d < 2^{32} + 1$，所以將情況對調，當 $d = 2^{32} + 1$ 時，$n$ 就必須小於 $2^{32}$，隨著 $d$ 減少，$n$ 能容許的範圍就越大。當 $d = 3$ (此例) 時，$n$ 在 $2^{63}$ 以下都能滿足條件。
從以上兩點可知此程式在 32 bits 的範圍下能夠正常運作，我們能將 $d$ 寫成更通用的形式，$d < \frac{2^N}{n} + 1$。當我們想要在更高位元的操作，只需要調整 $N$ 與 $n$ 來滿足需要的位元數。
#### 延伸問題
##### 1. 由 Facebook 公司所維護的 [jemalloc](https://github.com/jemalloc/jemalloc) 是個高效能的記憶體配置器 (memory allocator，即 malloc 和 free 函式對應的實作)，特別在多核多執行緒的環境有優異表現，在其原始程式碼 [include/jemalloc/internal/div.h](https://github.com/jemalloc/jemalloc/blob/dev/include/jemalloc/internal/div.h) 和 [src/div.c](https://github.com/jemalloc/jemalloc/blob/dev/src/div.c) 也運用類似的技巧，請閱讀原始程式碼，並抽離快速除法實作為獨立的程式碼，說明運作原理，也該撰寫測試程式，比較透過處理器的整數除法指令及本技巧的效能差距;
將 `jemalloc` 的除法操作抽取出來
```c
typedef struct div_info_s div_info_t;
struct div_info_s {
uint32_t magic;
};
void div_init(div_info_t *div_info, size_t divisor);
static inline size_t div_compute(div_info_t *div_info, size_t n) {
assert(n <= (uint32_t) - 1);
size_t i = ((uint64_t)n * (uint64_t)div_info->magic) >> 32;
return i;
}
void div_init(div_info_t *div_info, size_t d) {
assert(d != 0);
assert(d != 1);
uint64_t two_to_k = ((uint64_t)1 << 32);
uint32_t magic = (uint32_t)(two_to_k / d);
if (two_to_k % d != 0) {
magic++;
}
div_info->magic = magic;
}
```
其內容與本測驗中的實作非常相似，只有兩個地方不同：
1. `jemalloc` 使用的 $M$ (或是稱為 **Magic Number**) 為 $2^{32}$，而此測驗中用的為 $2^{64}$。
2. 在上面推導中出現的 $\lceil\frac{2^N}{d}\rceil$，因 c 語言沒有有效率的 ceil function，在 `jemalloc` 中使用餘數判斷是否為整除，而此測驗中則是轉換為 floor function 來實作。
:::danger
`jemalloc` 使用 $2^{32}$ 精度太低，註解只有寫道 We bound n < 2^32, and don't support dividing by one.，好像沒有關於 divisor 的限制。根據我在上面的推導，使用 $2^{32}$ 只能保證在 16 bits 下答案正確
:::
這裡給出其中一個例子:$\lfloor\frac{9510000}{16655}\rfloor=570$，但 `jemalloc` 會給出 $571$。
接著量測兩種除法的效率，分別是考慮計算 **Magic Number** 的時間與不考慮的時間。
分別測試 3 種除法之效率：
1. 內建除法運算子
2. 此測驗題中使用 floor function (轉型) 的除法，或稱為 [lkk ](https://arxiv.org/pdf/1902.01961.pdf)
3. `jemalloc` 使用的 `div`
為了讓結果更具有代表性，將本測驗題中的 $N$ 改為 $32$。我將從 $1\sim2^{20}$ 中以 $1000$ 的間隔抽取一次 $n$，並對該數字測試 1000 次來盡量縮小標準差。注意 $d$ 的選擇必須符合公式 $d < \frac{2^N}{n} + 1$。並且為了減少 function call 造成的時間差異 (雖然有 inline，但編譯器不一定會使用)，我將 `jemalloc` 的程式碼展開成 one line 的形式。並將 Magic Number 的精度改為 $2^{64}$
首先是不考慮 Magic Number 的除法效率，可以發現 `jemalloc` 與 `lkk` 在計算除法的部分是一樣的，所以只考慮 `jemalloc` 與除法運算子的效率就好。
測試程式：
```c
#include <stdint.h>
#include <stdbool.h>
#include <stdio.h>
#include <assert.h>
#include <time.h>
#include <stdlib.h>
#define MIN 3
#define MAX 0x80000000
int main() {
srand(time(NULL));
div_info_t div;
struct timespec tt1, tt2;
for(uint32_t i = MIN; i < MAX; i += 100000) {
size_t d = 0;
size_t ans[100001][2];
size_t times = 0;
time_t time1, time2;
while ((d = rand() % i) < 2)
;
jemalloc_init(&div, d);
clock_gettime(CLOCK_MONOTONIC, &tt1);
while (times++ <= 100000)
ans[times][0] = (i * (__uint128_t)div.magic) >> 64;
clock_gettime(CLOCK_MONOTONIC, &tt2);
time1 = (tt2.tv_nsec - tt1.tv_nsec);
times = 0;
clock_gettime(CLOCK_MONOTONIC, &tt1);
while (times++ <= 100000)
ans[times][1] = i / d;
clock_gettime(CLOCK_MONOTONIC, &tt2);
time2 = (tt2.tv_nsec - tt1.tv_nsec);
assert(ans[100000][0] == ans[100000][1]);
printf("test %d / %ld, time on jemalloc = %ldns, time on basic divide = %ldns\n", i, d, time1, time2);
}
}
```
我將 divisor 設為 i 以內，原因是我怕 divisor 大於 dividend 時編譯器會利用比大小來回傳(實際上會不會尚未研究)。每個測資皆跑 100000 次來減少標準差。另外我也將 `jemalloc` 的精度調整為 $2^{64}$。
可以看到 `-O0` 的效能差距相當大，~~而當編譯器加速後，兩者除法之間也有些許的差異。不管怎樣，在當 divisor 不用重新計算時，`jemalloc` 的效率皆比內建除法還好。~~
:::info
原先有測試 `O1~O3` 的效能，但發現編譯器加速的部分過於複雜，擔心測試不夠完整，故先移除。
:::
再來是考慮當 divisor 要重新計算時，計算 **Magic Number** 的除法效率
餘數判斷是非常耗時的操作，在 `jemalloc` 中的 `two_to_k % d != 0` 這個判斷式應該可以再簡化成判斷 `d` 是否為 2 的倍數，改為 `d & (d - 1) == 0` 來看看時間是否比較快。
所以這部分有 4 種待測除法
1. 內建除法運算子
2. `lkk` 使用的 `div`
3. `jemalloc` 使用的 `div`
4. 更改 `jemalloc` 後的 `div`
可以看到如果每次 divisor 都不同，計算 Magic Number 的成本比直接使用一般除法還高，~~而當編譯器加速得越多，不同方法之間的速度差也會越來越少。~~
這邊也順便證明了我去掉模運算的版本速度比原先 `jemalloc` 還要快上許多，但還是比不上 `lkk` 中使用轉型實現 floor function 的方法。
所以只要在固定的 divisor 上跑多次除法，先決定好 Magic Number 的值後，之後的除法效率就能夠超越內建的除法運算，但當每次 divisor 都會變動時，把計算 Magic Number 的成本加進來後，效率就可能比內建的除法運算差！
### `測驗 4`
```c
bool divisible(uint32_t n)
{
return n * M <= YYY;
}
```
#### `YYY`
- [x] `(c)` `M - 1`
#### 說明
證明參考 [sammer1107](https://hackmd.io/@sammer1107/ryHs1sESP#%E6%B8%AC%E9%A9%97-4) 同學的推導
* 當 $n \ {\mid} \ d$ 時，$n \times M <= M - 1$
* 當 $n \ {\nmid} \ d$ 時，$n \times M > M - 1$
分成兩個 case，$n \ {\mid} \ d$ 與 $n \ {\nmid} \ d$，算式內數值需對 $2^{64}$ 取餘數，因 type 為 `uint64_t`。
1. $n \ {\mid} \ d$，$n = 3q$：
$$
\begin{align}
n \times M &= 3q \times \lceil\frac{2^{64}}{3}\rceil \\
&= 3q \times \frac{2^{64} + 2}{3} \\
&= q2^{64} + 2q \quad (\text{mod} \ \ 2^{64})\\
&= 2q
\end{align}
$$
再來思考一下 $q$ 的上界：
$$
\begin{align}
n = 3q < 2^{32} \\
q < \frac{2^{32}}{3}
\end{align}
$$
帶回原式完成證明：
$$
\begin{align}
n \times M &= 2q \\
&< \frac{2^{33}}{3} \\
&< M - 1 \\
&= \frac{2^{64}+2}{3} - 1\\
&= \frac{2^{64} - 1}{3}
\end{align}
$$
2. $n \ {\mid} \ d$，$n = 3q + r，0<r<3$
我們將 $r$ 取 1，因 $3q + 1 > M - 1 \implies 3q + 2 > M - 1$：
$$
\begin{align}
n \times M &= (3q+r) \times \lceil\frac{2^{64}}{3}\rceil \\
&= (3q + 1) \times \frac{2^{64} + 2}{3} \\
&= q2^{64} + 2q + \frac{2^{64} + 2}{3} \quad (\text{mod} \ \ 2^{64})\\
&= 2q + \frac{2^{64} + 2}{3} \\
&\geq \frac{2^{64}+2}{3} \\
&> M - 1 \\
&= \frac{2^{64}-1}{3}
\end{align}
$$
$n \times M$ 是有可能等於 $M$ 的，如果 rhs 為 $M$ 而不是 $M - 1$ 的話，演算法有可能會有誤。如 $n = 1$ 時等號相等，`divisible` 會 return 1 而不是 0。
### `測驗 5`
```c
#include <ctype.h>
#include <stddef.h>
#include <stdint.h>
#define PACKED_BYTE(b) (((uint64_t)(b) & (0xff)) * 0x0101010101010101u)
/* vectorized implementation for in-place strlower */
void strlower_vector(char *s, size_t n)
{
size_t k = n / 8;
for (size_t i = 0; i < k; i++, s += 8) {
uint64_t *chunk = (uint64_t *) s;
if ((*chunk & PACKED_BYTE(VV1)) == 0) { /* is ASCII? */
uint64_t A = *chunk + PACKED_BYTE(128 - 'A' + VV2);
uint64_t Z = *chunk + PACKED_BYTE(128 - 'Z' + VV3);
uint64_t mask = ((A ^ Z) & PACKED_BYTE(VV4)) >> 2;
*chunk ^= mask;
} else
strlower(s, 8);
}
k = n % 8;
if (k)
strlower(s, k);
}
```
#### `VV1`
- [x] `(e)` `0x80`
#### `VV2`
- [x] `(b)` `0`
#### `VV3`
- [x] `(a)` `(-1)`
#### `VV4`
- [x] `(e)` `0x80`
#### 說明
先了解 `PACKED_BYTE` 的作用，從前面的 `& 0xff` 可知前面這部分只有前 1 個 byte 有值，再來後面的乘法可看成 8 個 bytes，每個 bytes 皆為 1，所以 `PACKED_BYTE` 可以知道為將 1 個 byte 的資料擴充成 8 個 bytes。
接著來看 `VV1`，從第一個 `if` 後方的註解 `is ASCII?` 可知該判斷式用來判斷是否為 ASCII code，下方 `else` 裡使用 `strlower`，應是用來處理 extend ASCII 的部分，所以要判斷是否為 ASCII code 的話 `VV1` 需選擇 0x80 (原因可參考測驗1)。
來看 `mask` 的用途，從 `*chunk ^= mask` 可以知道 `mask` 要將 `*chunk` 內 `A-Z` 的 ASCII code 轉為 `a-z`。根據[...]()，對每個 `A-Z` 字元的第 6 個 bit 做 XOR 運算就能轉成 `a-z`。故 `mask` 為在 `*chunk` 內每個為 `A-Z` 字元的部分的第 6 個 bit 為 1，其餘為 0。
由上述可知，`((A ^ Z) & PACKED_BYTE(VV4)) >> 2` 應該會在 `A-Z` 的字元中第 6 個 bit 產生 1，所以 `PACKED_BYTE(VV4)` 在第 8 個 bit 上產生 1，故 `VV4` 為 0x80。而 `A ^ Z` 就是用來檢查是否為 `A-Z` 間的字元，並將結果設置在第 8 個 bit 上。
要讓 XOR 為 True，必須一邊為 True，一邊為 False，但仔細觀察 `A` 與 `Z` 的程式碼，會發現在第 8 個 bit 上，`Z` 為 True 時，`A` 必定為 True。因 `Z` 比 `A` 還要多減了 `'Z' - 'A'` 的 ASCII code。
所以我們目的是讓 `A` 的第 8 個 bit 為 1，`Z` 的第 8 個 bit 為 0。則要讓 `128 - 'A' + VV2 + c` 內當 `c` 的 ASCII code 為 `A` 以上時，第 8 個 bit 為 1，`VV2` 需為 0，而在 `c` 的 ASCII code 大於 `Z` 時，`128 - 'Z' + VV3 + c` 的第 8 個 bit 要為 0，故 `VV3` 為 -1。
#### 延伸問題
##### 1. 將前述測試程式 main 函式的 char str[] 更換為 char *str，會發生什麼事？請參閱 C 語言規格書，說明 string literal 的行為
改成 `char *str` 會造成 Segmentation fault，試著用 pointer to char 指向已經建立 String Literals 上，`char *p_str = str`。程式卻能正常執行。
所以推測問題不是 pointer 的問題，而是建立變數時期的問題，從規格書裡找答案。
在規格書 [6.4.5 String literals](http://c0x.coding-guidelines.com/6.4.5.html) 中提到
> a byte or code of value zero is appended to each multibyte character sequence that results from a string literal or literals. The multibyte character sequence is then used to initialize an array of static storage duration and length just sufficient to contain the sequence.
>
> If the program attempts to modify such an array, the behavior is undefined.
在 [6.7.8 Initialization](http://c0x.coding-guidelines.com/6.7.8.html) 中也有寫道
>**`char *p = "abc";`**
defines p with type “pointer to char” and initializes it to point to an object with type “array of char” with length 4 whose elements are initialized with a character string literal. If an attempt is made to use p to modify the contents of the array, the behavior is undefined.
所以 `p` 為 "pointer to char" 的 type，指向一個 "array of char" 這個 object。在 [Frequently Asked Questions](http://c-faq.com/decl/strlitinit.html) 中也有提到，除了初始化時的變數 type 為 array of char 以外，string literal 會建立一個 unnamed, static array of characters 然後存放在 read-only 的記憶體位置，再將該位置回傳給 pointer to char 這個變數。
統整一下：
1. `char p[] = "hello";` 會被轉譯成 `char p[] = {'h', 'e', 'l', 'l', 'o', '\0'}`，所以操作該變數就跟一般的 array 操作一樣。
2. `char *p = "hello";` 此行為會在 read-only 的記憶體上建立 `hello` 這個字串，並將 `p` 指向該記憶體位置。
一個在 [stack overflow](https://stackoverflow.com/questions/1704407/what-is-the-difference-between-char-s-and-char-s) 上的例子我覺得能幫助理解：
```c
char *s0 = "hello world";
char s1[] = "hello world";
0x01 0x02 0x03 0x04
0x00008000: 'h' 'e' 'l' 'l'
0x00008004: 'o' ' ' 'w' 'o'
0x00008008: 'r' 'l' 'd' 0x00
...
s0: 0x00010000: 0x00 0x00 0x80 0x00
s1: 0x00010004: 'h' 'e' 'l' 'l'
0x00010008: 'o' ' ' 'w' 'o'
0x0001000C: 'r' 'l' 'd' 0x00
```
當使用 `char *s0 = "hello world";`時，
編譯器實際上做了以下事情
```c
static char __unnamed[] = "hello world";
char *c = __unnamed;
```
至於這個存放字串常量的 read-only 空間在哪裡呢？是放在一個名為 rodata 的地方，這邊 ro 指的就是 read-only，在程式內宣告的 `const` 變數，`#define` 的常量皆放置在此空間內，此位置與 [Data segment](https://en.wikipedia.org/wiki/Data_segment) 為兩種不同的空間。
可以寫個小程式來測試：
```c
#include <stdio.h>
char *s = "hello";
int main() {
char *str1 = "hello";
char *str2 = "world";
char str3[] = "!!";
return 0;
}
```
輸出 assembly code 來看編譯器是如何處理這些字串的
```
gcc -S str.c
```
可以看到以下程式碼：
```c
3 .globl s
4 .section .rodata
...
11 s:
12 .quad .LC0
13 .section .rodata
```
可以看到編譯器將字串常量放入 .rodata 裡面，且在 main 裡面的 `str1` 也會指到字串 `"hello"` 所在的地方，意思是 `s == str1`，而使用陣列方式宣告的 `char str3[]` 則是將字串內的值直接存入變數中
```c
35 movw $8481, -11(%rbp)
```
可以印出四個變數所指向的位置來檢查：
```c
s = 0x556e92e397a4
str1 = 0x556e92e397a4
str2 = 0x556e92e397aa
str3 = 0x7ffd2a9f0366
```
與預測的結果符合，`s` 與 `str1` 都指向在 rodata 中的 `"hello"` 字串位置。也可以注意到 `s, str1, str2` 三個位置都為 0x55 開頭，表示放在同樣的區塊，也就是 rodata 區塊，而 `str3` 與其他三著不同，是放置在 stack 區塊。
也可以透過 objdump 來查看 rodata 內的內容
```
objdump -s -j .rodata a.out
-------------------------------------
pointer: file format elf64-x86-64
Contents of section .rodata:
0730 01000200 68656c6c 6f00776f 726c6400 ....hello.world.
```
清楚的看到 `"hello"` 與 `"world"` 都被放到 rodata 中。
### `測驗 6`
```c
int singleNumber(int *nums, int numsSize)
{
int lower = 0, higher = 0;
for (int i = 0; i < numsSize; i++) {
lower ^= nums[i];
lower &= KKK;
higher ^= nums[i];
higher &= JJJ;
}
return lower;
}
```
#### `KKK`
- [x] `(c)` `~higher`
#### `JJJ`
- [x] `(d)` `~lower
`
#### 說明
要用 XOR 來比對出現一次與出現三次的差別，只用一個變數是不夠的，因一次 XOR 與三次 XOR 會造成相同結果。所以我們使用兩個變數 `lower` 與 `higher` 用來記錄出現一次與出現兩次。
我們先只看一個數 `nums[i]`，把 `lower` 與 `higher` 看成兩個 set，用來記錄使用過的 bits。
`lower` 有兩種情形，分別是 `lower` 存有 `nums[i]` 的 bits 與不存有 `nums[i]` 的 bits。經過 `lower ^= nums[i]` 後，`lower` 進到下個狀態，存有 `nums[i]` 的 bits 會被移除，沒有存 `nums[i]` 的話，則會將 `nums[i]` 加入 `lower` 中。
`higher` 同理，有兩種情形，存有 `nums[i]` 的 bits 與不存有 `nums[i]` 的 bits。故可以推測當 `nums[i]` 出現一次時，存到 `lower` 中，出現第二次時，從 `lower` 中移除，並加入至 `higher` 中，第三次出現則從 `higher` 中移除。
這樣答案就很清楚了，可以對照每行程式碼的行為：
1. `lower ^= nums[i]` 將 `nums[i]` 從 `lower` 中加入或移除
2. `lower &= ~higher` 如果 `nums[i]` 已經在 `higher` 當中(第二次出現)，則 `lower` 不需要再加入 `nums[i]`，將其移除
3. `higher ^= nums[i]` 將 `nums[i]` 從 `higher` 中加入或移除
4. `higher &= ~lower` 如果 `nums[i]` 已經在 `lower` 當中(第一次出現)，則 `higher` 不需要再加入 `nums[i]`，將其移除
畫成 state diagram 可以更好理解，可以想像當 `nums[i]` 不同時，出現 3 次的值還是會從 set 中移除，而出現 1 次的值永遠會存在 `lower` 當中，所以最後回傳的 `lower` 就是答案。
```graphviz
digraph G
{
rankdir=LR;
node node1 = [fixedsize=true, width = 2, height = 0.8, label = "lower = {}
higher = {}"] n1;
node node2 = [fixedsize=true, width = 2, height = 0.8, label = "lower = {nums[i]}
higher = {}"] n2;
node node3 = [fixedsize=true, width = 2, height = 0.8, label = "lower = {}
higher = {nums[i]}"] n3;
n1->n2 [label = "first times"];
n2->n3 [label = "second times"];
n3->n1 [label = "third times"];
}
```
#### 延伸問題
##### 1. 請撰寫不同上述的程式碼，並確保在 LeetCode [137. Single Number II](https://leetcode.com/problems/single-number-ii/) 執行結果正確且不超時;
我使用 64 bits 來儲存 32 bits 內每個 bit 出現的情況。
首先，我將每個數字擴充成 64 bits，並且中間間隔 1 個 0 用來存放出現 2 次以上的 bit。
類似下方表格這種轉換：
| number | 4 bits | 8 bits |
| -------- | -------- | -------- |
| 5 | <font color="blue">0101</font> | <font color="#f00">0</font><font color="blue">0</font><font color="#f00">0</font><font color="blue">1</font><font color="#f00">0</font><font color="blue">0</font><font color="#f00">0</font><font color="blue">1</font> |
| 14 | <font color="blue">1110</font> | <font color="#f00">0</font><font color="blue">1</font><font color="#f00">0</font><font color="blue">1</font><font color="#f00">0</font><font color="blue">1</font><font color="#f00">0</font><font color="blue">0</font> |
| 16 | <font color="blue">1111</font> | <font color="#f00">0</font><font color="blue">1</font><font color="#f00">0</font><font color="blue">1</font><font color="#f00">0</font><font color="blue">1</font><font color="#f00">0</font><font color="blue">1</font> |
這樣每個數字原本 32 bits 都擁有兩個 bit 可以儲存每個 bit 出現的次數，接著就是當出現 3 次，也就是該 bit 在擴充的位元裡出現 11 的時候，要將這兩個 bit 清空。
我先對現在累加的數字 `bit_set` 進行 `mask`，`mask` 的值為 `0x101010...1010`，用來檢查該 bit 是否累加超過 2 次以上。如果該結果為 1，再來右移 1 個 bit 檢查前一個 bit 是否也為 1，假如這兩個 bit 都為 1，就從 bit_set 內部清空。
最後再把結果的 64 bits 還原成 32 bits 並回傳。
```c
int singleNumber(int *nums, int numsSize)
{
uint64_t mask = 0xAAAAAAAAAAAAAAAA;
uint64_t bit_set = 0;
int ans = 0;
for (int i = 0; i < numsSize; ++i) {
uint64_t num_bits = 0;
for (int j = 0; j < 32; ++j) {
uint64_t offset = nums[i] & ((uint64_t)1 << j);
num_bits |= (offset << j);
}
bit_set += num_bits;
uint64_t t = (bit_set & mask) >> 1;
t = ((bit_set & t) << 1) | t;
bit_set &= ~t;
}
for (int i = 0; i < 32; ++i) {
ans |= (bit_set & ((uint64_t)1 << (i << 1))) >> i;
}
return ans;
}
```
##### 2. 延伸題意，將重複 3 次改為改為重複 4 次或 5 次，並找出可容易推廣到多次的程式碼;
詳細說明可參考[討論區](https://leetcode.com/problems/single-number-ii/discuss/43295/Detailed-explanation-and-generalization-of-the-bitwise-operation-method-for-single-numbers)，概念上與此題的解法相似，持續相加出現的 bit，直到出現的 bit 到達上限 `k`，將之移除。而 `k` 為每個數字出現的上限，可以知道，要記錄 `k` 個數字的 bit，必須至少要有 $\lceil\log_{2}k\rceil$ 個變數才行。
這邊 `k` 表示同樣數字出現的次數，以下為 `k` 任意數時迴圈內的程式碼：
```c
/* k == 2 */
x2 ^= x1 & i;
x1 ^= i;
mask = ~(x1 & x2);
x2 &= mask;
x1 &= mask;
/* k == 3 */
x3 ^= x2 & x1 & i;
x2 ^= x1 & i;
x1 ^= i;
mask = ~(x1 & ~x2 & x3);
x3 &= mask;
x2 &= mask;
x1 &= mask;
/* k == n */
xn ^= xn-1 & xn-2 ... & x1 & i;
xn-1 ^= xn-2 & xn-3 ... & x1 & i;
...
x1 ^= i;
mask = ~(x1 & x2 ... & xn);
xn &= mask;
xn-1 &= mask;
...
x1 &= mask;
```
要注意上面 `mask` 的值只是參考，實際上要隨著 `k` 值來決定 True 或 False，用來處理當 `xn, xn-1, ..., x1` 在某個 bit 上的值等於 `k` 時的情形，將這些 bit 清空。
假設表格內 `x5, x4, ..., x1` 都為 1 bit 的值，mask 的計算如下：
| k | $x_{5}$ | $x_{4}$ | $x_{3}$ | $x_{2}$ | $x_{1}$ | mask |
| --- | ------- | ------- | ------- | ------- | ------- | -------------------- |
| 1 | False | False | False | False | True | ~(x1) |
| 4 | False | False | True | False | False | ~(x3) |
| 7 | False | False | True | True | True | ~(x3 & x2 & x1) |
| 14 | False | True | True | True | False | ~(x4 & x3 & x2) |
| 29 | True | True | True | False | True | ~(x5 & x4 & x3 & x1) |
擴展到 32 bits 的 integer，當所有變數內的某個位置的 bit 與 `k` 的位元相符合時，表示已經達到上限了，將之從所有變數內刪除(NOT + AND)。
接著就是將這些步驟擴充成能隨著 `k` 來改變，底下程式碼邏輯與上面完全相同，只是用到了迴圈來解決變數數量不同的問題。
```c
int singleNumber(int* nums, int len, int k) {
uint32_t num_variable = ceil(log2(k)); //c99
uint32_t mask;
int* x = calloc(num_variable, sizeof(int));
for (int i = 0; i < len; ++i) {
mask = 0xFFFFFFFF;
for (int v = 0; v < num_variable; ++v) {
int tmp = nums[i];
for (int j = v + 1; j < num_variable; ++j) {
tmp &= *(x + j);
}
*(x + v) ^= tmp;
}
for (int v = 0; v < num_variable; ++v) {
mask &= ((1u << v) & k) > 0 ? *(x + v) : ~*(x + v);
}
mask = ~mask;
for (int v = 0; v < num_variable; ++v) {
*(x + v) &= mask;
}
}
int ans = 0;
for (int v = 0; v < num_variable; ++v)
ans |= *(x + v);
free(x);
return ans;
}
```
:::info
記得使用到 `math.h` 編譯時要輸入 `-lm` 連結至函式庫
:::
###### tags: `linux2020`
Sign in with Wallet
Connect another wallet