# 數值系統 閱讀紀錄

> author: < `Shiritai` >
>
> [原文見此](https://hackmd.io/@sysprog/c-numerics)
>
> 以下對該文核心做摘要 + 補充。

## 平衡的三進制

> [Donald E. Knuth](https://en.wikipedia.org/wiki/Donald_Knuth) 在《[The Art of Computer Programming](https://en.wikipedia.org/wiki/The_Art_of_Computer_Programming)》第 2 卷提到：
> > 也許平衡的三進制表示法是最美麗的數值系統。

透過 `+`, `0`, `-` 表示一個平衡三進制的位元:

* `+` 表 $1\times 3^k$
* `-` 表 $-1\times 3^k$
* `0` 表 $0\times 3^k=0$

此時有一些優美的性質:

> 複習[群論相關的離散數學](https://hackmd.io/t2E724Q_R8W3V3EEW-nu8w?both#%E7%94%A8%E5%BE%97%E5%88%B0%E7%9A%84%E4%BB%A3%E6%95%B8%E7%B3%BB%E7%B5%B1%E8%88%87%E7%BE%A4%E8%AB%96%E6%A6%82%E8%BF%B0)...
* 每個位元與加法擁有以下性質
* Identity: `0`
* Inverse: `+` 與 `-` 互為反元素，即「平衡」三進制的平衡兩字
* 另 Closed, Associative, Communicative 三性質要等具體設計好有限編碼的細節才知道。

可以看出，平衡三進制與生俱來（每個位元都）擁有二進制所沒有的 Inverse 這項性質，後者透過編碼以盡可能獲得此性質，這使平衡三進制的運算變得很簡單，比如反元素 (取負數) 變為 `+`, `-` 倒轉。

## 浮點數不是[阿貝爾群](https://hackmd.io/t2E724Q_R8W3V3EEW-nu8w?both#%E7%94%A8%E5%BE%97%E5%88%B0%E7%9A%84%E4%BB%A3%E6%95%B8%E7%B3%BB%E7%B5%B1%E8%88%87%E7%BE%A4%E8%AB%96%E6%A6%82%E8%BF%B0)

考慮浮點數編碼與加法

| Feature | Yes or no |
|:-------------:|:----------------------------:|
| Closed | Y |
| Associative | N |
| Identity | N (`0`, `-0`) |
| Inverse | Y (考慮 `INF`, `NaN` 的話 N) |
| Communicative | Y |

這使編譯器最佳化的前後可能影響浮點數的運算結果。

## 整數使用不當例子

* 2002 年 FreeBSD [53]

    #define KSIZE 1024
    char kbuf[KSIZE];
    int copy_from_kernel(void *user_dest, int maxlen) {
        int len = KSIZE < maxlen ? KSIZE : maxlen;
        memcpy(user_dest, kbuf, len);
        return len;
    }

為將 kernel 資料搬至 user space 的函式，在 `x86` 下若傳入值 `maxlen` 為負數，則 `len` 為負數，傳入 `memcpy` 第三參數 (被呼叫方解釋為 `size_t`) 後變一巨大的正數，如此可能將 kernel 不應該暴露的資料流出。
* 2002 年 External data representation (XDR) [62]

    void *copy_elements(void *ele_src[], int ele_cnt, int ele_size) {
        void *result = malloc(ele_cnt * ele_size);
        if (result==NULL) return NULL;
        void *next = result;
        for (int i = 0; i < ele_cnt; i++) {
            memcpy(next, ele_src[i], ele_size);
            next += ele_size;
        }
        return result;
    }

雖 `ele_cnt` 和 `ele_size` 都於整數範圍內，不過乘法可能導致溢位，進而影響 malloc 的大小不符合預期而複製 `ele_src` 到不該複製的位置。

## 位元運算的有趣應用

### 二補數中 `x & (x - 1) == 0` 的意義

> 透過去除「最低位 `1`」判斷是否為二的冪。

* 定義 `x`

$$
x := \{X...X\}_{part2}1\{0...0\}_{part1}
$$

* `x - 1` 後各個位元變成

$$
\rightarrow x - 1 = \{X...X\}_{part2}0\{1...1\}_{part1}
$$

其中 $\{\}$ 內兩部分的位元數前後相同。

* 則 `x & (x - 1)` 為 $\{X...X\}_{part2}0\{0...0\}_{part1}$，與 `x` 相比，表去除最低位的 $1$。
* 故 `x & (x - 1) == 0` 為真時，表示 $\{X...X\}_{part2}0\{0...0\}_{part1}$ 為零，也就是高位都為零。換言之，可用來判斷 $x$ 是否為二的冪。

### 位元反轉

unsigned char reverse(unsigned char n) {
    n = (n & 0xF0) >> 4 | (n & 0x0F) << 4;
    n = (n & 0xCC) >> 2 | (n & 0x33) << 2;
    n = (n & 0xAA) >> 1 | (n & 0x55) << 1;
    return n;
}

實作原理十分精妙，讓我來分析其中原委。

要完成位元反轉，一種做法是對 $i \forall i \in [0, 3]$ 位位元，將其位移至 $7-i$ 位，對高 $4$ 位同理。

也就是我們需要建立一個映射，使

$$
i \rightarrow 7 - i, \forall i \in [0, 7]
$$

對於 2 進制，我們可以看得更仔細：

(before:  i) 000 001 010 011 100 101 110 111
(after: 7-i) 111 110 101 100 011 010 001 000

映射結果即原本位置二進制的 bitwise flip。實作上，對於 $i$ 位置位元，為零時向 MSB 位移相應的距離，為一時向 LSB。從組成 $7$ 之二進制數 $4$ 開始，位置 $4$ 為零者向 MSB，位置 $4$ 為一者向 LSB。如果用 bit mask 來區別向 MSB 或向 LSB 的位元們，區別向 MSB 的 mask 應該長得像:

$$
0b00001111 = \tt{0x0F}
$$

同時區別向 LSB 的 mask 應該長得像:
$$
0b11110000 = \tt{0xF0}
$$

其他位置 ($2$, $1$) 同理，如此便解釋了前述玄妙的實作原理。

### ASCII 的編碼細節: 快速大小寫轉換

:::spoiler 用一小程式

void print_128bin(int a) {
    for (int i = 0; i < 8; ++i)
        printf("%d", (a >> (7 - i)) & 1);
}

void print_ascii() {
    for (int j = 0; j < 32; ++j) {
        for (int i = 0; i < 4; ++i) {
            int res = i * 32 + j;
            print_128bin(res);
            printf(" [%c]\t", (char) (isprint(res) ? res : '~'));
        }
        printf("\n");
    }
}

:::

:::spoiler 生成下表，`'~'` 表非列印值。

00000000 [~]    00100000 [ ]    01000000 [@]    01100000 [`]
00000001 [~]    00100001 [!]    01000001 [A]    01100001 [a]
00000010 [~]    00100010 ["]    01000010 [B]    01100010 [b]
00000011 [~]    00100011 [#]    01000011 [C]    01100011 [c]
00000100 [~]    00100100 [$]    01000100 [D]    01100100 [d]
00000101 [~]    00100101 [%]    01000101 [E]    01100101 [e]
00000110 [~]    00100110 [&]    01000110 [F]    01100110 [f]
00000111 [~]    00100111 [']    01000111 [G]    01100111 [g]
00001000 [~]    00101000 [(]    01001000 [H]    01101000 [h]
00001001 [~]    00101001 [)]    01001001 [I]    01101001 [i]
00001010 [~]    00101010 [*]    01001010 [J]    01101010 [j]
00001011 [~]    00101011 [+]    01001011 [K]    01101011 [k]
00001100 [~]    00101100 [,]    01001100 [L]    01101100 [l]
00001101 [~]    00101101 [-]    01001101 [M]    01101101 [m]
00001110 [~]    00101110 [.]    01001110 [N]    01101110 [n]
00001111 [~]    00101111 [/]    01001111 [O]    01101111 [o]
00010000 [~]    00110000 [0]    01010000 [P]    01110000 [p]
00010001 [~]    00110001 [1]    01010001 [Q]    01110001 [q]
00010010 [~]    00110010 [2]    01010010 [R]    01110010 [r]
00010011 [~]    00110011 [3]    01010011 [S]    01110011 [s]
00010100 [~]    00110100 [4]    01010100 [T]    01110100 [t]
00010101 [~]    00110101 [5]    01010101 [U]    01110101 [u]
00010110 [~]    00110110 [6]    01010110 [V]    01110110 [v]
00010111 [~]    00110111 [7]    01010111 [W]    01110111 [w]
00011000 [~]    00111000 [8]    01011000 [X]    01111000 [x]
00011001 [~]    00111001 [9]    01011001 [Y]    01111001 [y]
00011010 [~]    00111010 [:]    01011010 [Z]    01111010 [z]
00011011 [~]    00111011 [;]    01011011 [[]    01111011 [{]
00011100 [~]    00111100 [<]    01011100 [\]    01111100 [|]
00011101 [~]    00111101 [=]    01011101 []]    01111101 [}]
00011110 [~]    00111110 [>]    01011110 [^]    01111110 [~]
00011111 [~]    00111111 [?]    01011111 [_]    01111111 [~]

:::

可觀察到大小寫之差在第五位的有無，故大小寫轉換可以透過

* 添加第五位: 將大寫轉小寫 (原本即小寫者不變)
> `' '`: 純第五位

    ('a' | ' ') // 得到 'a'
    ('A' | ' ') // 得到 'a'

* 刪去第五位: 將小寫轉大寫 (原本即大寫者不變)
> `'_'`: 沒有第五位

    ('a' & '_') // 得到 'A'
    ('A' & '_') // 得到 'A'

* 反轉第五位: 大小寫調換

    ('a' ^ ' ') // 得到 'A'
    ('A' ^ ' ') // 得到 'a'

### 以 XOR 調換數字

> 我第一次接觸到此做法是在[劉教授的用片語學 C 語言程式設計](https://sites.google.com/view/c-programming-2ed/home)

[該書中](https://www.youtube.com/watch?v=3bCrd4soChg&list=PLOvZ8aEg7xDmRtTbeWHGKAkXI8ZkcyAhe&index=132)提到這種 Ninja Code 的存在:

a ^= b ^= a ^= b;
// X= 從右邊結合
a ^= (b ^= (a ^= b));

以上成立是因為 XOR 的交換率，對位元 $a_n, b_n$ 而言，以上操作即

$$
a_n' = a_n \oplus b_n
$$

$$
b_n' = b_n \oplus a_n' = b_n \oplus a_n \oplus b_n \\
= a_n \oplus b_n \oplus b_n = a_n \oplus 0 = a_n
$$

$$
a_n'' = a_n' \oplus b_n' = (a_n \oplus b_n) \oplus (a_n) = b_n \oplus 0 = b_n
$$

故所有位元都得到交換，且全程 in-place，即不另外使用暫存器或記憶體。

> 另外，劉教授提到不應為了炫技而寫出不易閱讀的程式碼，不過我認為若其在實務上有價值，且透過正確地理解原理以及家事適當的註解，亦可活用該程式碼。

### 相加除二: 避免溢位

(x & y) + ((x ^ y) >> 1)

`x & y` 為「進位」，`x ^ y` 為「位元和」，將後者左移一位，前者留在原地，正是將加法結果除以二的表現。

### `0` 的偵測

[原文](https://hackmd.io/@sysprog/c-numerics#%E9%81%8B%E7%94%A8-bit-wise-operator)以 top-down 的方式解釋 [`DETECT`](https://github.com/eblot/newlib/blob/master/newlib/libc/string/strlen.c) 巨集的行為。

#if LONG_MAX == 2147483647L
#define DETECT(X) \
    (((X) - 0x01010101) & ~(X) & 0x80808080)
#else
#if LONG_MAX == 9223372036854775807L
#define DETECT(X) \
    (((X) - 0x0101010101010101) & ~(X) & 0x8080808080808080)
#else
#error long int is not a 32bit or 64bit type.
#endif
#endif

我想以 buttom-up 的方式作補充 。

假設對於 $1\ \rm{byte}$ 而言，我們希望確認其是否為零，可以怎麼做？

還記得前面提到我們有[確認數字是否為二的冪](https://hackmd.io/y_jVeGzDQhetNS6A-XxPmA?both#unsigned-%E4%B8%AD-x-amp-x---1--0-%E7%9A%84%E6%84%8F%E7%BE%A9)的邏輯嗎？當中利用 $x$ 和 $x-1$ 的關係進行位元運算來消除最低位的 `1`。

複習一下這兩數的外型 (1 byte):

* 定義 $x$

$$
x := \{X...X\}_{part2}1\{0...0\}_{part1}
$$

其中**第一部分的長度可為 0~8**。

* $x - 1$ 後各個位元變成

$$
\rightarrow x - 1 = \{X...X\}_{part2}0\{1...1\}_{part1}
$$

注意若 $x$ 為零，則實際上 $x$ 只有第一部分而已，此時 $x-1$ 為全一。若 $x$ 非零，則 $x-1$ 至少有一個位元為 $0$ ，緊跟在第一部分之上。不仿想辦法取出這個有標示性意義的 $0$，假設我們最後希望該零存在時為真，不存在時為假...

* 所有 $x$ 一定有第一部分（長度自 $0$ 至 $8$），設法忽略掉。
* 第二部分不重要，設法忽略掉。想忽略掉相同的東西，常見兩方法
* 取 xor (`x ^ (x - 1)`): 在此副作用是使第一部分和標示性位元全變成 $1$，混淆視聽。
* 取 not 後 and: 正好在消除前面的同時，把 $x$ 的標示性位元留下，有兩種取法
* `(x - 1) & ~x`: $\{0...0\}_{part2}0\{1...1\}_{part1}$
* $x = 0 \rightarrow 11111111$，僅此時最高位為 $1$
* $x \ne 0 \rightarrow 0...01...1$，最高位以降至少有一個 $0$
* `(x - 1) & ~x & 0x80` 為真，表示 $x=0$，符合需求 :)
* `~(x - 1) & x`: $\{0...0\}_{part2}1\{0...0\}_{part1}$
* $x = 0 \rightarrow 00000000$
* $x \ne 0 \rightarrow 0...010...0$，某一位為一，只是不知道是哪位
* `(~(x - 1) & x) == 0` 為真，表示 $x=0$，符合需求 :)

問題來了，這倆哪個比較好？就要看架構了。

根據上方的成果，如果應用在類似 SIMD 的場景，比如將八個字元以 64 bits 一次讀入運算，比如[原文給的例子](http://opensource.apple.com/source/Libc/Libc-583/arm/string/strlen.s)，以及[舊版 glibc 實作 strlen](https://github.com/lattera/glibc/blob/master/string/strlen.c)，希望高效的尋找 ``'\0'`` 的位置，就可以活用此技巧。

> 以上 `glibc` 的程式碼為舊版，[新版](https://github.com/bminor/glibc/blob/master/string/strlen.c)程式碼使用巨集技巧讓程式碼變得更加優雅。

## 位元細緻操作

> [`NachOS`](https://en.wikipedia.org/wiki/Not_Another_Completely_Heuristic_Operating_System) 中的 [`bitmap.c`](https://github.com/Shiritai/NTHU-Operating-System-Eroiko_Blue/blob/master/MP4/NachOS-4.0_MP4/code/lib/bitmap.cc) 應該是不錯的例子。順帶一提，過去實作作業系統作業時，我曾對其部分操作的複雜度進行優化。

### Set a bit

void Bitmap::Mark(int which)
{
    map[which / BitsInWord] |= 1 << (which % BitsInWord);
}

### Clear a bit

void Bitmap::Clear(int which)
{
    map[which / BitsInWord] &= ~(1 << (which % BitsInWord));
}

### Flip a bit

void Bitmap::Flip(int which)
{
    map[which / BitsInWord] ^= 1 << (which % BitsInWord);
}

### Test a bit

bool Bitmap::Test(int which) const
{
    return !!map[which / BitsInWord] & (1 << (which % BitsInWord));
}