---
tags: linux2026
---
# [2026q1](https://wiki.csie.ncku.edu.tw/linux/schedule) 第 4 週測驗題
> [作答表單](https://forms.gle/oAJ3Cn2rLjA5iXtq7)
:::warning
:warning: 回答延伸問題時，務必優先以[課程教材](https://wiki.csie.ncku.edu.tw/linux/schedule)為主要出處，隨後參照 C 語言規格書、Linux 核心原始程式碼及其 git log 和 LKML (Linux Kernel Mailing-List)，和權威素材 (如 GCC 和 glibc 手冊，但絕對不是 CSDN 或者尋常的網路日誌/筆記)
AI 工具是輔助性質，可用來撰寫測試程式碼和收集資訊，但主要的推測、查證、分析，和討論，都該由人類進行。
:::
回應網友在 Ptt BBS 八卦板的提問: [Re: [問卦] 用C語言排出張惠妹的臉很強嗎？](https://www.ptt.cc/bbs/Gossiping/M.1773653862.A.BF3.html)，文中的 [amei.c](https://gist.github.com/jserv/d701cc0d17a0e8f0ddba020cbeaa3221) 藉由內建的影像解壓縮程式，得以輸出天后的容顏為 PNG 檔案。
```
cc -O2 amei.c -lm && ./a.out > amei.png
```
影像資料取自[鏡週刊](https://www.mirrormedia.mg/story/20251104edi024)。
> 故意降低品質，避免侵權。提升畫質的任務，交給學員練習
## 整體架構
有損影像壓縮的主體思路是：先預測 (prediction) 畫素值，再對預測誤差做正交轉換 (orthogonal transform) 以集中能量，接著量化 (quantize) 轉換係數來丟棄人眼不 敏感的高頻細節，最後以熵編碼 (entropy coding) 無失真壓縮剩餘的符號。從資訊理論 (information theory) 的角度看，熵編碼的目標是使每個符號的平均碼長趨近其熵 (entropy)：
$$
H = -\sum_{i} p_i \log_2 p_i
$$
從線性代數的角度看，DCT 是對訊號向量做一組餘弦基底 (cosine basis) 的正交投影 (orthogonal projection)，量化則是將係數向量投射到較粗的格點 (lattice) 上。
解碼器分為兩層：
* 上層：有損 (lossy) 影像解碼，涵蓋四元樹 (quadtree) 分割、預測、DCT 殘差 (residual)
* 下層：無損 (lossless) 熵編碼，涵蓋二元算術編碼 (binary arithmetic coding) 與指數哥倫布編碼 (exponential-Golomb coding)
資料流向與各階段的關鍵運算如下：
| 階段 | 函式 | 關鍵運算 | 算術類型 |
|------|------|---------|---------|
| 位元串流解碼 | `init()` | `(hi << 4) \| lo` 半位元組合併 | 位元操作 |
| 算術解碼 | `decode_bit()` | `range *= 256` 重正規化、區間分割 | 整數乘除 |
| 整數解碼 | `decode_unsigned()` | `value << 1 \| bit` 位元建構 | 位元操作 |
| 係數反量化 | `extract_block()` | `level × q` 均勻中踏量化 | 整數乘法 |
| IDCT | `compute_idct()` | 定點乘加、`>> scale` 縮放 | Q10 定點 |
| 預測補償 | `extract_block()` | `& 7` / `>> 3` 次畫素內插 | 定點內插 |
| 去方塊濾波 | `deblock()` | `>> 2` / `>> 3` 加權混合 | 整數位移 |
| 雙邊濾波 | `bilateral3()` | 梯度門檻、加權平均 | 整數算術 |
| 色彩轉換 | `write_png()` | YCgCo → RGB 加減法 | 整數加減 |
| PNG 輸出 | `write_png()` | CRC32 nibble 查表、Adler-32 | 位元操作 |
整條管線不使用浮點運算，所有數值處理以整數算術和位元操作完成。計算密度最高的階段是 IDCT (每個 $N \times N$ 區塊需要 $2N^3$ 次乘加) 和方向性預測 (每個畫素需要 2 次參考畫素查詢和 1 次內插)。
## 輸入資料與十六進位解碼
### 內嵌資料編碼
`amei.c` 的影像資料內嵌於字串 `amei_hex[]`，採用十六進位 (hexadecimal)
編碼。每 2 個字元對應 1 個位元組 (byte)：
```c
/* amei.c:56-66 init() */
const char *s = amei_hex;
uint8_t *o = amei_data;
while (*s) {
int hi = *s <= '9' ? *s - '0' : *s - 'a' + 10;
int lo = s[1] <= '9' ? s[1] - '0' : s[1] - 'a' + 10;
*o++ = (hi << 4) | lo;
s += 2;
}
```
此處的位元操作 (bitwise operation) 是典型的半位元組合併手法：`hi << 4` 將高半位元組 (high nibble) 左移 4 位元，再以位元或 (bitwise OR) `|` 與低半位元組 (low nibble) 合併，組成一個完整的 8 位元值。解碼後得到 776 位元組的原始位元串流。
### DCT 餘弦表預計算
`init()` 的另一個關鍵工作是建立離散餘弦轉換 (Discrete Cosine Transform, DCT) 所需的餘弦查詢表 (cosine lookup table)。DCT-III (即 IDCT) 的主體需要餘弦值，但此處的查詢表包含 DCT 的標準正規化因子 $\sqrt{2}$：
$$
\text{dct\_coef}[k] = \begin{cases} 1024 & k = 0 \\ \left\lfloor 1024\sqrt{2} \cdot \cos\!\left(\dfrac{\pi k}{64}\right) + 0.5 \right\rfloor & k = 1, \dots, 127 \end{cases}
$$
$k > 0$ 的項使用縮放因子 $1024\sqrt{2} \approx 1448$，而 $k = 0$ (DC 項) 單獨設為 $1024$。這對應 DCT-III 的正規化慣例：DC 基底的振幅為 $1/\sqrt{N}$，AC 基底的振幅為 $\sqrt{2/N}$，兩者差一個 $\sqrt{2}$ 倍。以浮點運算表示，等效的參考實作為：
```c
/* 等效的浮點參考實作 (使用 <math.h>) */
#include <math.h>
for (int k = 0; k < 128; k++)
dct_coef[k] = (int) lrint(1024.0 * sqrt(2.0) * cos(M_PI * k / 64.0));
dct_coef[0] = 1024; /* DC 項不乘 sqrt(2) */
```
`amei.c` 改以整數遞推 (integer recurrence) 產生相同的查詢表，完全消除對 `<math.h>` 的依賴：
```c
/* amei.c:69-77 init() */
int a = 0, b = 74509276;
for (int i = 0; i < 128; i++) {
dct_coef[(i + 96) & 127] = ((a >> 19) + 1) >> 1;
int c = b;
b = (int) ((2144896910LL * b >> 30) - a);
a = c;
}
dct_coef[0] = 1024;
```
這段程式碼利用三角函數的遞推關係。正弦和餘弦共用相同的三項遞推式：
$$
f((k+1)\theta) = 2\cos\theta \cdot f(k\theta) - f((k-1)\theta)
$$
其中 $f$ 可為 $\sin$ 或 $\cos$，由初始條件決定。此處 $\theta = \pi/64$，初始條件 $a_0 = 0$ 對應 $\sin(0) = 0$，而 $a_1 = b$ 對應 $\sin(\pi/64)$，因此遞推產生的是正弦序列
$\sin(k\theta)$。經由索引偏移 96，$\sin(k\pi/64)$ 被存入
$\text{dct\_coef}[(k+96) \bmod 128]$，而 $\cos(m\pi/64) = \sin((m+32)\pi/64)$ (餘弦是正弦的相位偏移)，因此最終查詢表中的值等效於
$$
1024\sqrt{2} \cdot \cos(m\pi/64)
$$
程式碼中的常數對應：
- 初始值 `b = 74509276`：代表 $S \cdot \sin(\pi/64)$，
其中內部定點精度 $S = 1024\sqrt{2} \cdot 2^{20} \approx 2^{30.5}$，
故 $74509276 \approx 2^{30.5} \times \sin(\pi/64) \approx 1518500250 \times 0.04907$
- 乘數 `2144896910`：代表 $2\cos\theta$ 以 Q30 定點數表示，
即 $2144896910 / 2^{30} \approx 1.99759 \approx 2\cos(\pi/64)$
- 運算 `2144896910LL * b >> 30`：計算 $2\cos\theta \cdot a_k$，
64 位元乘法避免溢位，右移 30 位元還原定點格式
- 減法 `... - a`：完成遞推 $a_{k+1} = 2\cos\theta \cdot a_k - a_{k-1}$
- 輸出 `((a >> 19) + 1) >> 1`：將內部精度 $S \approx 2^{30.5}$ 縮放到
輸出格式，右移 19 位元後加 1 再右移 1 位元實現四捨五入，
等效於 $\lfloor a / 2^{20} + 0.5 \rfloor$
索引對映 `(i + 96) & 127` 利用位元及 (bitwise AND) 實作模 128 (mod 128) 的循環索引，相當於 `(i + 96) % 128`，但對 2 的冪 (power of two) 而言，位元及運算比取餘 (modulo) 更有效率。偏移 96 將正弦序列對映為餘弦值：
$\sin(k\pi/64)$ 存入索引 $(k + 96) \bmod 128$，使得查詢表位置 $m$ 存放的值等於 $\sin((m - 96)\pi/64) = \cos((m - 64)\pi/64) = \cos(m\pi/64 - \pi)$，經餘弦週期性化簡後正好是 IDCT 所需的 $\cos(m\pi/64)$ (含符號)。最終 `dct_coef[0] = 1024` 覆蓋遞推產生的 $1448 \approx 1024\sqrt{2}$，
將 DC 項的 $\sqrt{2}$ 因子移除，符合 DCT-III 的正規化慣例。
兩種實作 (浮點和整數遞推) 產生完全相同的 128 項查詢表。整數版本避免浮點硬體依賴，適合 Linux 核心和嵌入式環境。
## 二元算術編碼
算術編碼 (arithmetic coding) 是整個壓縮系統的基石，負責將位元串流無損還原為個別的位元和整數。Shannon 的無雜訊編碼定理 (noiseless coding theorem) 證明，最佳無損編碼的平均碼長下界為訊源熵 $H$。Huffman 編碼受限於每個符號至少 1 位元的整數位長，而算術編碼可以逼近分數碼長 (fractional bit length)，對機率極度偏斜的二元訊源 (binary source) 尤其有效：解碼單一位元消耗的輸入
量約為 $-\log_2 P$ 位元，其中 $P$ 為該位元的預測機率。
### 主體原理
算術編碼將 $[0, 1)$ 區間 (interval) 依據符號的機率 (probability) 遞迴分割。對二元算術編碼而言，區間被分成兩段 (分別對應 0 和 1)，大小正比於各自的預測機率。編碼時，每讀入一個符號就縮小區間；解碼時，追蹤目前的分數 (fraction) 落入哪一段，即可還原對應的位元值。
精確追蹤區間的上下界需要任意精度算術。解碼器的做法是：僅保留區間的大小 `range` 和目前位置 `fraction` 這兩個整數，以有限精度近似區間。當 `range` 縮小到不足以區分符號時，便進行重正規化 (renormalization)，將 `range` 和 `fraction` 同時乘以基數 (radix)，並從輸入讀入新的位元組補充精度。
### 解碼實作
`amei.c` 中的 `decode_bit()` 函式實作算術解碼器 (arithmetic decoder)：
```c
/* amei.c:79-102 decode_bit() */
static int decode_bit(int context)
{
if (range < 256) {
range *= 256;
fraction *= 256;
if (inp < inp_end)
fraction += *inp++;
}
int *counts = models + context * 2;
int observations = *counts + counts[1] + 2,
split = range * (*counts + 1) / observations,
bit = fraction >= split;
if (bit) {
fraction -= split;
range -= split;
} else {
range = split;
}
counts[bit]++;
if (observations > 63) {
*counts /= 2;
counts[1] /= 2;
}
return bit;
}
```
關鍵變數：
- `range`：目前區間大小，初始為 1
- `fraction`：目前分數位置，初始為 0
- `models[166]`：83 組上下文容器 (context bin)，每組 2 個計數器
(0-bit 和 1-bit 的觀測次數)
運作流程分為四個步驟：
1. 重正規化：當 `range < 256` 時，將 `range` 和 `fraction` 各乘以 256 (即左移 8 位元)，並從輸入讀取下一個位元組加到 `fraction`。這裡的 256 即為基數，對應一個位元組的 256 種可能值。乘以 256 等效於 `<< 8`，為精度騰出 8 位元的空間。
2. 機率估計 (probability estimation)：取出上下文 `context` 的計數器 `counts[0]` (0 的次數) 和 `counts[1]` (1 的次數)，以 Laplace 接續規則 (Laplace's rule of succession) 計算分割點 (split point)：
`split = range * (*counts + 1) / observations`。
Laplace 規則在真實計數上各加 1 (反映在 `observations = *counts + counts[1] + 2`)，使 0-bit 的估計機率為 $\hat{p}_0 = (c_0 + 1) / (c_0 + c_1 + 2)$，確保即使某個值從未出現過，其預測機率也不會為零，因為零機率會使該符號永遠無法被解碼。
3. 解碼：比較 `fraction >= split` 決定位元值。若為 1，取上半段 (`fraction -= split; range -= split`)；若為 0，取下半段 (`range = split`)。 這就是區間的遞迴分割，每解碼一個位元，區間就縮小為原來的一部分。
4. 更新與週期性縮放 (periodic scaling)：遞增被解碼位元的計數器。當 `observations > 63` 時，兩個計數器都右移一位 (`/= 2`)，使近期觀測的權重更大。這等效於指數加權移動平均 (exponential weighted moving average, EWMA)：每次縮放都將舊觀測的影響減半，使模型能適應位元串流中統計特性隨位置變化的非穩態 (non-stationary) 特性。
### 上下文模型 (context modeling)
83 個上下文容器 (`models[166]` 即 83 對計數器) 依用途分配如下：
| 上下文 | 用途 | 對應程式碼 |
|--------|------|-----------|
| 0-2 | 四元樹分割決策 (8x8, 16x16, 32x32) | `decode_bit(level - 3)` |
| 3 | DCT 係數正負號 (sign bit) | `decode_bit(3)` |
| 4 | 無號整數 (unsigned integer) 的位元值 | `decode_bit(4)` |
| 5-14 | 亮度 DCT 零遊程 (zero run) 字首長度 | `decode_unsigned(5)` |
| 15-24 | 色度 DCT 零遊程字首長度 | `decode_unsigned(15)` |
| 25-34 | 高頻亮度 DCT 量化等級 | `decode_unsigned(25)` |
| 35-44 | 高頻色度 DCT 量化等級 | `decode_unsigned(35)` |
| 45-54 | 低頻亮度 DCT 量化等級 | `decode_unsigned(45)` |
| 55-64 | 低頻色度 DCT 量化等級 | `decode_unsigned(55)` |
| 61-72 | DCT 係數結束旗標 (end-of-block flag) | `decode_bit(61 + level*2 + is_chroma)` |
| 73-78 | 預測模式 (prediction mode) | `decode_unsigned(73)` |
低頻與高頻的分界以 `coeff < area / 8` 判定，即前 1/8 的係數視為低頻，其餘為高頻。分離上下文讓低頻和高頻係數各自維護獨立的機率模型，因為自然影像中低頻係數通常較大、高頻係數通常較小，混在同一個上下文會降低預測準確度。
## 指數哥倫布編碼
函式 `decode_unsigned()` 在算術解碼之上建構無號整數解碼：
```c
/* amei.c:104-113 decode_unsigned() */
static int decode_unsigned(int context)
{
int length = 0;
while (!decode_bit(context + length))
length++;
int value = 1;
while (length--)
value = value << 1 | decode_bit(4);
return value - 1;
}
```
解碼分兩階段。第一階段計數前導零 (leading zeros)：反覆呼叫 `decode_bit()`，使用上下文 `context + length`，直到解碼出一個 1 為止。前導零的個數 `length` 決定後續要讀取多少資料位元。
第二階段以位元建構 (bit construction) 組合數值：`value << 1 | decode_bit(4)` 是經典的位元操作手法，左移一位 (`<< 1`) 騰出最低位元 (least significant bit, LSB)，再以位元或 `|` 將新解碼的位元填入。所有資料位元統一使用上下文 4，與字首的上下文分開。最後減 1 還原為原始值。
編碼規則：將值加 1，寫成二進位表示，字首 (prefix) 補零的個數等於二進位去掉最高位元後的位數。
| 值 | 二進位 (值+1) | 編碼 |
|----|-------------|------|
| 0 | 1 | 1 |
| 1 | 10 | 010 |
| 2 | 11 | 011 |
| 3 | 100 | 00100 |
| 7 | 1000 | 0001000 |
小數值用較少位元，天然適合影像資料中大量趨近零的係數。以原版 Lena 影像為參考，3501 個整數中有 2420 個 (69%) 為零，僅需 1 個位元；整體平均每個整數消耗約 2.11 個位元。
## 影像標頭 (image header)
`main()` 中首先引導 (bootstrap) 算術解碼器，再解碼四個參數：
```c
/* amei.c:332-350 main() */
range *= 256;
fraction *= 256;
if (inp < inp_end)
fraction += *inp++;
int level = decode_unsigned(5);
width = 1 << level;
int height = width - decode_unsigned(5);
luma = decode_unsigned(5);
chroma = decode_unsigned(5);
```
- `level`：影像寬度的 $\log_2$ 值，`1 << level` 以左移運算 (left shift) 計算
$2^{\text{level}}$，得到寬度 (必為 2 的冪)
- `height`：影像高度 = 寬度 - 差值 (僅支援正方形或橫向影像)
- `luma`：亮度量化步長 (luma quantization step size)
- `chroma`：色度量化步長 (chroma quantization step size)
較大的量化步長代表更高的壓縮率 (compression ratio) 但更低的影像品質。最小的合法量化步長為 1，此時近似於無失真壓縮。
## YCgCo 色彩空間
自然影像的 RGB 三通道之間存在高度相關性 (correlation)：相鄰畫素的 R、G、B 值往往同步變化。從線性代數的觀點，對 RGB 向量做去相關轉換 (decorrelation transform) 等同於將資料投影到主成分 (principal component) 方向，使共變異數矩陣 (covariance matrix) $\text{Cov}(X_i, X_j) \to 0 \;(i \ne j)$，從而讓各通道可以獨立壓縮而不損失效率。
壓縮資料並非直接以 RGB 儲存，而是使用 YCgCo 色彩空間：
- Y：亮度 (luma)，攜帶絕大部分影像細節
- Cg：綠-品紅色度 (green-magenta chrominance)
- Co：橙-藍色度 (orange-blue chrominance)
YCgCo 的轉換矩陣為：
$$
\begin{bmatrix} Y \\ C_g \\ C_o \end{bmatrix} = \begin{bmatrix} \frac{1}{4} & \frac{1}{2} & \frac{1}{4} \\[4pt] -\frac{1}{4} & \frac{1}{2} & -\frac{1}{4} \\[4pt] \frac{1}{2} & 0 & -\frac{1}{2} \end{bmatrix} \begin{bmatrix} R \\ G \\ B \end{bmatrix}
$$
此矩陣近似正交 (approximately orthogonal)，且僅涉及加減和位移，不需乘法。YCgCo 將亮度與色度解相關 (decorrelate)，人眼對亮度遠比色度敏感，因此編碼器可以對色度通道使用較大的量化步長而不明顯損失視覺品質。YCgCo 因只需整數加減法與位移（無須浮點乘法），具備完全無損可逆 (lossless reversible) 的特性，這也是 H.264/HEVC 等現代視訊標準正式採納它的原因。解碼器不做色度次取樣 (chroma subsampling)，而是直接提高色度的量化係數。
三個通道以平面格式 (planar format) 分別儲存在記憶體中，最終從 YCgCo 轉回交錯格式 (interleaved format) 的 RGB 輸出：
```c
/* amei.c:302-308 write_png() 中的色彩轉換 */
int i = y * width + x;
int t = image[0][i] - image[1][i];
U8AD(CLAMP8(t + image[2][i])); /* R = (Y - Cg) + Co */
U8AD(CLAMP8(image[0][i] + image[1][i])); /* G = Y + Cg */
U8AD(CLAMP8(t - image[2][i])); /* B = (Y - Cg) - Co */
```
其中 `image[0]` 為 Y，`image[1]` 為 Cg，`image[2]` 為 Co。
`CLAMP8()` 巨集 (macro) 以三元條件運算子將值夾限 (clamp) 至 [0, 255] 範圍：
`((v) < 0 ? 0 : (v) > 255 ? 255 : (v))`。
## 四元樹分割
四元樹 (quadtree) 是一種遞迴空間分割 (recursive spatial partitioning) 結構：每個內部節點恰好有四個子節點，分別對應矩形區域的四個象限。
「四元」強調空間象限的分割，「四叉」則強調樹的分支因子 (branching factor)。從圖論 (graph theory) 的角度，四元樹是一棵完滿四叉樹 (full 4-ary tree)，每個節點要不有恰好四個子節點，要不就是末端節點；但因適應性分割使各末端節點出現在不同深度，並非完全四叉樹 (complete 4-ary tree)。深度 $d$ 最多可表示 $4^d$ 個末端節點。
相比 JPEG 固定使用 8x8 區塊，四元樹以每個內部節點僅 1 位元的開銷 (分割或不分割) 實現適應性區塊大小 (adaptive block size)：細節豐富的區域遞迴分割為小區塊，平坦區域則保持大區塊。
函式 `extract_block()` 是解碼的主體，以遞迴四元樹方式拆分影像。
### 分割策略
```c
/* amei.c:136-144 extract_block() */
static void extract_block(int left, int top, int level)
{
if (level > 5 || (level > 2 && decode_bit(level - 3))) {
int size = 1 << --level;
for (int corner = 0; corner < 4; corner++)
extract_block(left + corner % 2 * size,
top + corner / 2 * size, level);
return;
}
/* ... 末端節點解碼 ... */
}
```
- 區塊 (block) 大於 32x32 (`level > 5`)：強制分割
- 區塊為 4x4 (`level == 2`)：不再分割
- 中間層級 (8x8, 16x16, 32x32)：解碼一個位元決定是否分割
分割時，`1 << --level` 以左移運算計算子區塊大小。四個子區塊的座標以位元操作推導：`corner % 2` 取最低位元 (LSB) 決定水平偏移，`corner / 2` 取次低位元決定垂直偏移，按左上、右上、左下、右下順序遞迴處理。這產生類似 Z-order 曲線的走訪順序，使空間相鄰的區塊在編碼串流中也相鄰，有利於壓縮效率。
區塊尺寸分布大致為：4x4 約 360 個、8x8 約 122 個、16x16 約 11 個、32x32 為 0 個。Amei 影像使用前述的四元樹分割機制，細節豐富的區域 (如五官) 被切得更細，平坦區域 (如背景) 保持較大區塊。
## 預測機制
預測編碼 (predictive coding) 的數學基礎是條件期望值：給定已知的鄰近畫素，預測目前畫素最小化均方誤差 (MSE) 的最佳估計即為條件期望
$\hat{x} = E[x \mid \text{neighbors}]$。實務上使用線性預測器 (linear predictor) 近似，而殘差 $r = x - \hat{x}$ 的分布比原始畫素值更集中於零附近 (峰度 (kurtosis) 更高)，用更少的位元即可編碼。預測越準確，殘差的熵 $H(r)$ 越低，壓縮效果越好。
到達四元樹的末端節點 (leaf node) 後，解碼器對 Y、Cg、Co 三個通道分別處理。首先解碼一個預測模式 `predictor`，此值在 0-33 之間，三個通道共用。
### DC 預測 (模式 0)
```c
/* amei.c:162-169 extract_block() 內 */
if (!predictor) {
int dc = 0;
for (int pixel = 0; pixel < size; pixel++) {
dc += top ? output[-width + pixel] : 0;
dc += left ? output[pixel * width - 1] : 0;
}
*buffer += left && top ? dc / 2 : dc;
}
```
取目前區塊正上方和正左方已解碼畫素 (pixel) 的平均值，寫入 DCT 的 DC 係數。這是個精簡的技巧：不必逐畫素填入預測值，只需修改 DC 係數，因為 DC 係數在 IDCT 後恰好對應整個區塊的常數偏移 (constant offset)。
### 方向性預測 (directional prediction) (模式 1-33)
HEVC 規範共有 35 個影格內預測模式 (0=Planar, 1=DC, 2-34 為方向性)，此處的 33 個方向性預測模式精確對應其中的方向性子集（剔除 Planar，以模式 0 為 DC，1-33 為方向性）。每個方向性模式以 1/8 畫素為步進定義一個斜率 (slope)，圍繞對角線 (模式 17) 對稱分布。
```c
/* amei.c:174-175 extract_block() 內 */
int is_leftward = predictor < 17,
slope = is_leftward ? 9 - predictor : predictor - 25;
```
- 模式 1-16 (`is_leftward=1`)：主要從左方畫素向右延伸
- 模式 17-33 (`is_leftward=0`)：主要從上方畫素向下延伸
- 模式 9：水平複製左方畫素
- 模式 25：垂直複製上方畫素
- 模式 17：沿對角線方向
對於非整數位置 (sub-pixel position)，解碼器以 1/8 畫素精度在兩個相鄰參考畫素 (reference pixel) 間做線性內插 (linear interpolation)，產生抗鋸齒 (anti-aliasing) 效果：
```c
/* amei.c:180-193 extract_block() 內 */
int ref = y * slope + slope;
weight = ref & 7;
ref = (ref >> 3) + x + nb;
/* ... 參考畫素查詢與邊界處理 ... */
output[is_leftward ? x * width + y : y * width + x] +=
(*buffer * (8 - weight) + buffer[1] * weight + 4) >> 3;
```
`ref & 7` 以位元及取出最低 3 位元作為內插權重 (0-7)，`ref >> 3` 右移 3 位元取得整數部分的參考畫素索引。`>> 3` 等效於除以 8，配合 `(8 - weight)` 和 `weight` 構成以 1/8 為精度的雙線性內插。`+ 4` 在右移前加入四捨五入的偏移量。
## 殘差解碼與 DCT
預測只是近似值。要還原實際影像，還需要加上殘差 (residual)。
殘差以離散餘弦轉換 (Discrete Cosine Transform, DCT) 係數的形式儲存。
DCT 可視為離散傅立葉轉換 (DFT) 的實數版本。對長度 $N$ 的訊號向量
$\mathbf{x}$，DCT-II 定義為：
$$
X[k] = \sum_{n=0}^{N-1} x[n] \cos\!\left(\frac{\pi(2n+1)k}{2N}\right), \quad k = 0, 1, \dots, N-1
$$
其基底向量構成正交矩陣 $\mathbf{C}$，滿足 $\mathbf{C}^T \mathbf{C} = \mathbf{I}$
(單位矩陣)。正交性保證轉換不改變向量的 $L^2$ 範數 (Parseval 定理，$\|\mathbf{X}\|^2 = \|\mathbf{x}\|^2$)，因此量化誤差在頻域和空間域的能量相等。
二維 DCT 因可分離性 (separability) 可拆為先對列再對行 (或反向) 各做一次一維 DCT，將 $O(N^4)$ 降為 $O(N^3)$。自然影像的能量大多集中於低頻 DCT 係數，高頻係數趨近零，這種能量壓縮特性 (energy compaction) 使 DCT 成為有失真壓縮的標準工具。
### 係數解碼
對每個區塊的每個色彩通道 (color channel)：
```c
/* amei.c:149-161 extract_block() 內 */
memset(buffer, 0, (size_t) area * sizeof(int));
for (int coeff = 0; coeff < area; coeff++) {
if (decode_bit(61 + level * 2 + is_chroma))
break; /* 結束旗標 */
coeff += decode_unsigned(5 + is_chroma * 10); /* 零遊程長度 */
int sign = 1 - 2 * decode_bit(3); /* 正負號 */
int lv = (decode_unsigned(
25 + (is_chroma + (coeff < area / 8) * 2) * 10) + 1)
* (channel ? chroma : luma);
if (coeff < area)
buffer[coeff] = sign * lv;
}
```
流程類似 H.263 的 (LAST, RUN, LEVEL) 三元組 (triplet)：
1. 解碼一個位元判斷是否已到最後一個非零係數
2. 解碼零遊程 (run of zeros)，即連續跳過多少個零係數
3. 解碼正負號 (sign bit)：`1 - 2 * decode_bit(3)` 以乘法將 {0, 1} 對映為
{1, -1}，避免分支
4. 解碼量化等級 (quantized level)，乘以量化步長 (亮度用 `luma`，色度用
`chroma`) 還原係數，此為均勻中踏量化 (mid-tread uniform quantization)，
無死區 (dead zone)
上下文選擇中的 `(coeff < area / 8) * 2` 利用布林值轉整數的隱式轉換：比較結果為 0 或 1，乘以 2 後偏移上下文索引，使低頻 (前 1/8) 和高頻係數使用不同的機率模型。
此處的均勻中踏量化 (mid-tread) 將係數 $c$ 對映為
$\hat{c} = \lfloor c/q + 0.5 \rfloor$，其中 $q$ 為量化步長。
反量化 (dequantization) 則是簡單的乘法 $c' = \text{level} \times q$。與 JPEG 使用的逐係數量化矩陣 (quantization matrix) 不同，這個解碼器對同一通道的所有係數使用單一步長，程式碼更精簡，但犧牲對各頻率分量的獨立控制。從率失真理論
(rate-distortion theory) 的角度，對高頻使用較大步長可以在相同位元率 (bitrate) 下降低整體失真 $D$，這裡則透過上下文分離來間接達成類似效果。
係數按列主序 (row-major order) 排列，不同於 JPEG 的鋸齒形掃描 (zigzag scan)。非零係數傾向集中於左上角 (低頻區域)，自然影像的 DCT 係數具有高度稀疏性 (sparsity)。
### 反離散餘弦轉換 (Inverse DCT, IDCT)
函式 `compute_idct()` 實作一維 IDCT，透過可分離性 (separability) 拆為二次呼叫完成二維 IDCT：
```c
/* amei.c:170-171 extract_block() 內 */
compute_idct(buffer + area, 1, buffer, 1, size, size, 10);
compute_idct(output, width, buffer + area, size, 1, size, 10 + level);
```
第一次呼叫對各列 (row) 做一維 IDCT，從 `buffer` 的前半 (頻域係數) 讀入，結果寫到後半。第二次呼叫對各行 (column) 做一維 IDCT，從後半讀入，直接寫回影像平面 `output`。步進參數 (stride) 控制讀寫方向：第一次以 `1` 和 `size` 步進處理列，第二次以 `width` 和 `1` 步進處理行。
IDCT 的主體迴圈使用預計算的餘弦表：
```c
/* amei.c:126-133 compute_idct() */
for (int freq = 0; freq < size; freq++)
sum += coefficients[freq * frequency_stride +
seq * sequence_stride] *
dct_coef[(2 * pt + 1) * freq * DCT_MAX / size % DCT_MAX4];
samples[pt * point_stride + seq * sequence_stride] =
(sum + (1 << (scale - 1))) >> scale;
```
`(2 * pt + 1) * freq * DCT_MAX / size % DCT_MAX4` 以整數算術計算餘弦表的索引，對應數學式 $(2n+1)k \cdot 32/N \bmod 128$，即 DCT-III 主體的
$\cos\!\bigl(\frac{\pi(2n+1)k}{2N}\bigr)$ 在查詢表中的位置。
`% DCT_MAX4` (即 `% 128`) 利用餘弦函數的週期性 ($\cos$ 以 $2\pi$ 為週期) 實作環繞存取。`(1 << (scale - 1))` 在右移 `scale` 位元前加入 $0.5$ 的四捨五入偏移，`>> scale` 執行定點數到整數的縮放，等效於
$$
\lfloor x / 2^{\text{scale}} + 0.5 \rfloor
$$
DCT 的本質是將區塊分解為一組正弦波基底圖樣 (basis pattern) 的加權和 (weighted sum)。左上角的低頻 (low-frequency) 基底對應模糊的漸層，右下角的高頻 (high-frequency) 基底對應細節邊緣。量化 (quantization) 就是在此處引入有失真壓縮：除以量化步長並四捨五入，小係數歸零，減少需要編碼的資料量。
原版使用 DCT-II 變體編碼，解碼器對應實作 DCT-III 進行反轉換。
原版的 IDCT 是整個解碼器中唯一使用浮點運算的地方 (`cos()` 和 `lrint()`)，其餘全部採用定點或整數算術。`amei.c` 以預計算的整數餘弦表取代浮點運算，完全消除對 `<math.h>` 的依賴。
## 後處理 (post-processing)
區塊式壓縮的固有缺陷是方塊效應 (blocking artifact)：各區塊獨立量化導致邊界處出現不連續。從訊號處理的觀點，方塊效應等同於在區塊邊界引入高頻雜訊。後處理的目的是在不過度模糊細節的前提下，以低通濾波 (low-pass filtering) 消除這類人為雜訊。解碼器使用兩階段策略：先在區塊邊界做邊緣感知的去方塊濾波 (edge-aware deblocking)，再對整個影像做選擇性雙邊濾波 (selective bilateral filter)。
### 去方塊濾波
函式 `deblock()` 在區塊邊界 (block boundary) 做平滑處理，消除因獨立編碼各區塊而產生的方塊效應 (blocking artifact)：
```c
/* amei.c:201-219 deblock() */
for (int bs = 4; bs <= 64; bs <<= 1)
for (int dir = 0; dir < 2; dir++) {
/* ... */
int a = p[i - step], b = p[i];
if (abs(a - b) <= 4)
continue;
p[i - step] = (3 * a + b + 2) >> 2;
p[i] = (a + 3 * b + 2) >> 2;
/* ... */
p[i - 2 * step] = (7 * p[i - 2 * step] + a + 4) >> 3;
/* ... */
p[i + step] = (7 * p[i + step] + b + 4) >> 3;
}
```
`bs <<= 1` 以左移逐級倍增區塊大小 (4, 8, 16, 32, 64)，走訪所有四元樹層級的邊界。邊界兩側的畫素以 `>> 2` (除以 4) 實現 3:1 的加權混合，更外側的畫素以 `>> 3` (除以 8) 實現 7:1 的微調。`+ 2` 和 `+ 4` 分別是除以 4 和除以 8 時的四捨五入偏移量。
### 選擇性雙邊濾波 (selective bilateral filter)
函式 `bilateral3()` 對平坦區域做 3x3 加權平均濾波 (weighted average filtering)，降低壓縮雜訊 (compression noise)，
同時保留邊緣 (edge)：
```c
/* amei.c:223-257 bilateral3() */
int idx = y * w + x, c = tmp[idx], grad = 0;
/* 計算中心畫素與四鄰的梯度 (gradient) */
if (x > 0) grad += abs(c - tmp[idx - 1]);
if (x + 1 < w) grad += abs(c - tmp[idx + 1]);
if (y > 0) grad += abs(c - tmp[idx - w]);
if (y + 1 < h) grad += abs(c - tmp[idx + w]);
if (grad > grad_limit)
continue; /* 梯度過大 (邊緣) 則跳過 */
/* 對 3x3 鄰域 (neighborhood)，只混合色差在門檻值 (threshold) 內的畫素 */
int sw = (2 - abs(dx)) * (2 - abs(dy));
```
權重 `sw` 以曼哈頓距離 (Manhattan distance) 計算：中心畫素 4，四鄰 2，對角 1。最終以 `(sum + wt / 2) / wt` 做加權平均，`wt / 2` 提供四捨五入。
亮度通道使用較小的門檻值 (threshold=8, grad_limit=42)，色度通道使用較大的門檻值 (threshold=14, grad_limit=56)，反映色度通道量化更粗糙、雜訊更多的特性。
## PNG 輸出
PNG 格式需要：
1. 8 位元組檔頭簽名 (signature)
2. IHDR 區塊 (chunk)：寬度、高度、位元深度 (bit depth)、色彩類型
3. IDAT 區塊：zlib 壓縮的畫素資料
4. IEND 區塊：結束標記
`amei.c` 使用一組巨集 (`PUT`, `U8C`, `U32C`, `U16LC`, `U8AD`,
`BEGIN`, `END`) 實作精簡的 PNG 寫入器：
```c
/* amei.c:277-283 write_png() 中的巨集 */
#define U8C(u) do { unsigned _v=(u)&255; PUT(_v); \
crc^=_v; crc=(crc>>4)^ct[crc&15]; crc=(crc>>4)^ct[crc&15]; } while(0)
#define U8AD(u) do { unsigned _a=(u)&255; U8C(_a); \
ad_a=(ad_a+_a)%65521; ad_b=(ad_b+ad_a)%65521; } while(0)
```
- `U8C()` 在輸出位元組的同時計算 CRC32。`crc & 15` 以位元及取出最低 4 位元
作為查表索引，`crc >> 4` 右移 4 位元後與查表結果做位元互斥或 (XOR) `^`，
執行兩輪 nibble 查表，等效於逐位元組的 CRC 更新。查詢表 `ct[]` 僅 16 項，
比標準的 256 項表節省大量空間。
- `U8AD()` 在 `U8C()` 基礎上同時計算 Adler-32 總和檢查碼 (checksum)，
以 `% 65521` (最大的小於 $2^{16}$ 的質數) 維護兩個累加器。
- zlib 資料採用不壓縮的 stored block 模式 (`0x78 0x01` 標頭)，每列一個區塊。
採用 stored block 模式讓解碼器完全省去實作 Deflate（LZ77 + Huffman）壓縮
演算法的龐大程式碼開銷，依然能封裝出標準相容的 PNG 檔案。
`U16LC()` 以小端序 (little-endian) 輸出 16 位元長度和其一補數 (one's complement)。
## 資料規模
以 Amei 影像為例：
- `amei_hex[]` 含 1552 個十六進位字元，每 2 字元解碼 1 位元組，
共 776 位元組的有效位元串流
- 解碼後影像：128x128 畫素 RGB
## 定點算術與位元操作
整個解碼器不使用浮點運算，所有數值處理以定點算術和位元操作完成。以下彙整各模組使用的定點格式與位元操作模式 (pattern)，並探討後續最佳化的方向。
### 定點數格式一覽
解碼器在不同模組使用不同精度的定點數：
| 模組 | 定點格式 | 精度 (位元) | 說明 |
|------|---------|------------|------|
| 餘弦表 | Q10 (含 $\sqrt{2}$) | 10 位元小數 | `dct_coef[]`，AC 項含 $\sqrt{2}$ 正規化 |
| 餘弦遞推內部 | Q30.5 | ~30.5 位元小數 | `a`, `b` 變數，$S = 1024\sqrt{2} \cdot 2^{20} \approx 2^{30.5}$ |
| 遞推乘數 | Q30 | 30 位元小數 | `2144896910` $\approx 2\cos(\pi/64) \times 2^{30}$ |
| IDCT 中間值 | Q10+level | 可變 | `scale = 10` (列) 或 `10 + level` (行) |
| 方向性預測 | Q3 | 3 位元小數 | 1/8 畫素精度，`& 7` 取小數、`>> 3` 取整數 |
| 去方塊濾波 | Q2 / Q3 | 2-3 位元 | `>> 2` (3:1 混合)、`>> 3` (7:1 混合) |
Q 格式的命名慣例：Q$n$ 代表 $n$ 位元小數精度，整數值 $v$ 對應實數 $v / 2^n$。定點乘法 $a \times b$ 的結果為 Q$(m+n)$，需右移還原到目標精度。
### 位元操作模式分類
解碼器中的位元操作可歸為以下模式：
1. 位元欄位合併 (bit field merging)：`(hi << 4) | lo`、`value << 1 | bit`。
左移騰出空間，位元或填入新值。這是建構多位元數值的基本手法，消除乘法和加法的需求。
2. 善用 2 的冪除法 (power-of-two division)：`>> scale`、`>> 2`、`>> 3`。
右移 $n$ 位元等效於除以 $2^n$ 並向下取整 (floor division)。對有號整數，C 語言的右移行為在早期標準中屬於依實作定義 (implementation-defined)；此解碼器假設算術右移 (arithmetic right shift)，即符號位元 (sign bit) 填充，這在 x86、ARM、RISC-V 等主流架構上成立。自 C++20 與 C23 標準起，有號整數的二補數 (two's complement) 表示法與算術右移已被正式標準化，使得這類最佳化在現代編譯器上具備完全的可攜性 (portability)。
3. 2 的冪取餘 (power-of-two modulo)：`& 127`、`& 7`。
對 $2^n$ 取餘等效於保留最低 $n$ 位元：$x \bmod 2^n = x \mathbin{\&} (2^n - 1)$。這比一般的取餘指令快一個數量級。
4. 四捨五入偏移 (rounding bias)：`(a >> 19) + 1) >> 1`、`(sum + 4) >> 3`。
右移前加入 $2^{n-1}$ 的偏移量，將向下取整轉為四捨五入：
$\lfloor x / 2^n + 0.5 \rfloor = (x + 2^{n-1}) \gg n$。
5. 無分支符號對映 (branchless sign mapping)：`1 - 2 * bit`。
將 $\{0, 1\}$ 對映為 $\{+1, -1\}$，避免條件分支 (conditional branch)。類似的手法也見於 `(coeff < area / 8) * 2`，利用布林值隱式轉換為 0 或 1 來偏移索引。
### 最佳化方向
目前的實作以可讀性為優先，留有幾個最佳化空間：
1. IDCT 蝴蝶分解 (butterfly decomposition)：目前 `compute_idct()` 使用 $O(N^2)$ 的直接矩陣乘法。對 $N = 4, 8, 16, 32$，可改用蝴蝶形快速演算法
將複雜度降為 $O(N \log N)$。例如 8 點 DCT 僅需 11 次乘法 + 29 次加法 (Loeffler 演算法)，取代直接實作的 64 次乘法 + 56 次加法。Linux 主體的 vicodec 使用 Walsh-Hadamard Transform (WHT)，其蝴蝶結構完全消除乘法 (詳見後文「與 Linux 主體的對照」一節)。
2. 位移取代乘法：算術解碼器中的 `range * (*counts + 1) / observations` 涉及整數除法 (通常比乘法慢 10-40 倍)。若將 `observations` 限制為 2 的冪，可用右移取代除法，但會犧牲機率估計的精度。另一個做法是以查表近似倒數乘法 (reciprocal multiplication)：$x / d \approx x \times (2^k / d) \gg k$。
3. SIMD 向量化：IDCT 的內層迴圈對 $N$ 個頻率分量做乘加累積 (MAC)，適合運用 SIMD 指令 (如 ARM NEON 的 `vmla`、x86 SSE 的 `_mm_madd_epi16`)。Q10 定點數可裝入 16 位元，每個 128 位元暫存器可同時處理 8 個乘加。去方塊濾波和雙邊濾波的畫素迴圈同樣可向量化。
4. 餘弦表壓縮：128 項的查詢表可利用餘弦的對稱性 (偶函數) 和反對稱性
($\cos(\pi - x) = -\cos(x)$) 壓縮為 33 項，以少量位元操作 (符號翻轉、索引對映) 換取快取 (cache) 友善性。HEVC 規格使用類似策略，以 6 位元定點數儲存 32 項餘弦值。
## 與 Linux 主體的對照
Linux 主體本身不包含軟體 DCT/IDCT 實作。硬體 JPEG 編解碼器驅動程式 (如 `drivers/media/platform/renesas/rcar_jpu.c`、
`drivers/media/platform/samsung/s5p-jpeg/`) 將量化表和 Huffman 表寫入硬體暫存器，實際的 DCT/IDCT 在專用晶片上執行。
軟體正交轉換位於 vicodec 測試驅動程式：
### vicodec FWHT：純加減法的正交轉換
`drivers/media/test-drivers/vicodec/codec-fwht.c` 實作 8x8 快速 Walsh-Hadamard 轉換 (Fast Walsh-Hadamard Transform, FWHT)。WHT 是 DCT 的變形，同為正交轉換，但基底函數僅使用 $\{+1, -1\}$ 而非餘弦波，因此完全不需乘法。
正向轉換 `fwht()` 的列處理：
```c
/* drivers/media/test-drivers/vicodec/codec-fwht.c fwht() */
/* 三級蝴蝶：每級只有加法和減法 */
/* 第一級：相鄰配對 */
workspace1[0] = tmp[0] + tmp[1] - add;
workspace1[1] = tmp[0] - tmp[1];
workspace1[2] = tmp[2] + tmp[3] - add;
workspace1[3] = tmp[2] - tmp[3];
/* 第二級：跨步配對 */
workspace2[0] = workspace1[0] + workspace1[2];
workspace2[1] = workspace1[0] - workspace1[2];
/* 第三級：更大跨步 */
out[0] = workspace2[0] + workspace2[4];
out[1] = workspace2[0] - workspace2[4];
```
反向轉換 `ifwht()` 結構完全相同 (WHT 為自反轉換，$H^{-1} = H / N$)，僅在最後以 `>>= 6` 除以 $64 = 8 \times 8$ 做正規化。
### 量化：純位移取代除法
vicodec 的量化表以位移量 (shift count) 表示，避免任何除法指令：
```c
/* drivers/media/test-drivers/vicodec/codec-fwht.c */
static const int quant_table[] = {
2, 2, 2, 2, 2, 2, 2, 2, /* 低頻：>> 2 (除以 4) */
2, 2, 2, 2, 2, 2, 2, 2,
/* ... */
2, 2, 2, 3, 6, 6, 6, 6, /* 高頻：>> 6 (除以 64) */
2, 2, 3, 6, 6, 6, 6, 8, /* 最高頻：>> 8 (除以 256) */
};
/* 量化：*coeff >>= *quant; 反量化：*de_coeff = *coeff << *quant; */
```
高頻位置使用較大的位移量，等效於較粗的量化步長，與 JPEG 使用非均勻量化矩陣的策略相同，但完全以 2 的冪實作，消除除法。
### 三者對照
| 特性 | amei 解碼器 | vicodec FWHT | JPEG (ITU-T.81) |
|------|-----------|-------------|-----------------|
| 轉換 | DCT (餘弦基底) | WHT (Walsh 基底) | DCT (餘弦基底) |
| 區塊大小 | 4x4 至 32x32 (四元樹) | 固定 8x8 | 固定 8x8 |
| 乘法需求 | $2N^2$ 次/區塊 (直接法) | 0 次 (純加減) | $2N^2$ 或快速法 |
| 定點精度 | Q10 (10 位元小數) | 無小數 (純整數) | 視實作而定 (libjpeg 用 Q13) |
| 量化 | 單一步長 × 通道 | 位移表 (2 的冪) | 64 項量化矩陣 |
| 係數掃描 | 列主序 | 鋸齒形 (zigzag) | 鋸齒形 (zigzag) |
| 熵編碼 | 算術 + 指數哥倫布 | RLC + 位元封裝 | Huffman 或算術 |
| 浮點依賴 | 無 | 無 | 視實作而定 |
amei 解碼器在精神上更接近 HEVC 的影格內解碼子集 (可變區塊大小、方向性預測、算術編碼)，而 vicodec 的 FWHT 則走 JPEG 的固定 8x8 路線。兩者共享一個主體設計原則：以整數算術和位元操作消除浮點依賴，使程式碼可在無浮點硬體的環境 (如早期嵌入式處理器或 Linux 主體空間) 執行。
vicodec 的 WHT 將這個原則推到極致：不僅沒有浮點，連整數乘法都省去，蝴蝶分解將 8 點轉換拆為 $3 \times 4 = 12$ 次加減法 (對比直接法的 64 次乘加)。amei 解碼器的 DCT 因使用餘弦基底而無法完全消除乘法，但 Chebyshev 遞推在建表階段以 1 次 64 位元乘法取代 128 次 `cos()` 呼叫，同樣體現以算術技巧換取硬體獨立性的思路。
## 設計思維
這個解碼器的設計展現幾個精巧的取捨：
1. 四元樹取代固定大小區塊：僅需每個節點一個位元的額外開銷 (overhead)，就能適應性配置解析度 (resolution)，平坦區域用大區塊、細節區域用小區塊，比 JPEG 的固定 8x8 更靈活。
2. DC 預測寫入 DCT 係數：避免逐畫素迴圈，將空間域 (spatial domain) 的常數偏移直接對映到頻域 (frequency domain)，程式碼極度精簡。
3. 上下文適應性算術編碼 (context-adaptive arithmetic coding)：83 個上下文容器各自追蹤局部統計，配合週期性縮放機制適應資料的非穩態特性。
4. 可分離 IDCT 與整數餘弦表：一個函式二次呼叫完成二維轉換，以預計算的定點餘弦表完全消除浮點運算依賴。
5. 位元操作貫穿全程：從十六進位解碼的 `<< 4 | lo`、算術解碼的 `*= 256` 重正規化、指數哥倫布的 `<< 1 | bit` 位元建構、IDCT 的 `>> scale` 定點縮放、到方向性預測的 `& 7` 和 `>> 3` 次畫素內插，以位元操作取代除法、取餘、浮點運算等高成本操作，同時精確控制定點數的精度和捨入行為。
整個解碼器完成從熵解碼、係數反量化 (dequantization)、IDCT、幀內預測到後處理濾波的完整影像解碼管線 (pipeline)，可視為 HEVC 解碼器的微型子集。
---
本系列測驗以 `amei.c` 為題材，檢驗數值系統、bitwise operation、線性代數、機率、C 語言規格及程式效率等方面的認知。
## Problem A: 定點數值系統與正交轉換
- [ ] Part 1: Chebyshev 遞推、$\sqrt{2}$ 正規化與相位轉換
`init()` 以三項遞推產生 128 項 DCT 餘弦查詢表：
```c
/* amei.c:69-77 init() */
int a = 0, b = 74509276;
for (int i = 0; i < 128; i++) {
dct_coef[(i + 96) & 127] = ((a >> 19) + 1) >> 1;
int c = b;
b = (int) ((2144896910LL * b >> 30) - a);
a = c;
}
dct_coef[0] = 1024;
```
查詢表包含 DCT-III 的標準正規化因子 $\sqrt{2}$ (參閱 `amei.md` §DCT 餘弦表預計算)：
$$\text{dct\_coef}[k] = \begin{cases} 1024 & k = 0 \\[4pt] \left\lfloor 1024\sqrt{2} \cdot \cos\!\left(\dfrac{\pi k}{64}\right) + 0.5 \right\rfloor & k = 1, \dots, 127 \end{cases}$$
DC 基底的振幅為 $1/\sqrt{N}$，AC 基底的振幅為 $\sqrt{2/N}$，兩者差一個 $\sqrt{2}$ 倍。
* 遞推的初始條件 $a_0 = 0$ 對應 $\sin(0) = 0$，因此遞推產生正弦序列。乘數 `2144896910` 以 Q30 格式表示 $2\cos(\pi/64)$：驗算 $2144896910 / 2^{30} \approx$ __ A01 __ (小數點後五位)。初始值 `b = 74509276` 代表 $S \cdot \sin(\pi/64)$，其中內部定點精度 $S = 1024\sqrt{2} \cdot 2^{20} \approx 2^{\_\_A02\_\_}$ (小數點後一位)。
* 索引映射 `(i + 96) & 127`：遞推將 $\sin(k\pi/64)$ 存入 `dct_coef[(k+96) mod 128]`，使得位置 $m$ 存放 $\sin((m - 96)\pi/64)$。利用 $\cos(x) = \sin(x + \pi/2)$ 且 $\pi/2$ 對應 32 個表項，偏移 96 的效果等於相位推進 $\pi/2$，因為 $96 \equiv$ __ A03 __ $(\text{mod } 128)$。
* 輸出轉換 `((a >> 19) + 1) >> 1` 等效於 $\lfloor a / 2^{20} + 0.5 \rfloor$，將內部精度 $S \approx 2^{30.5}$ 縮放到 Q10 輸出。`dct_coef[0] = 1024` 覆蓋遞推產生的 $\lfloor 1024\sqrt{2} + 0.5 \rfloor \approx 1448$，移除 DC 項的 $\sqrt{2}$ 因子，使 `compute_idct()` 不需對 `freq == 0` 做條件分支。
驗證 `dct_coef[1]`：$1024\sqrt{2} \cdot \cos(\pi/64) \approx 1448.15 \times 0.99879 \approx$ __ A04 __ (取整)。
* 若將遞推常數改為 `2147483648LL` ($= 2^{31}$)，$2\cos\theta = 2$，$\theta = 0$。特徵方程 $r^2 - 2r + 1 = 0$ 有重根 $r = 1$，通解 $f(k) = C_1 + C_2 k$。以 $f(-1) = 0, f(0) = B$ 代入，$f(k) = $ __ A05 __ (以 $B, k$ 表示)。此線性成長序列使查詢表的值單調遞增直到溢位，IDCT 完全失去頻率選擇性。
- [ ] Part 2: IDCT 索引推導與可分離性
```c
/* amei.c:126-133 compute_idct() */
for (int freq = 0; freq < size; freq++)
sum += coefficients[freq * frequency_stride +
seq * sequence_stride] *
dct_coef[(__A06__) * freq * __A07__ / size % __A08__];
samples[pt * point_stride + seq * sequence_stride] =
(sum + (1 << (scale - 1))) >> scale;
```
DCT-III 基底函數 $\cos\!\bigl(\frac{\pi(2n+1)k}{2N}\bigr)$ 對應角度 $\alpha = \frac{(2n+1)k\pi}{2N}$。查詢表以 128 項覆蓋 $2\pi$ 週期 (`DCT_MAX4 = 128`)，每項對應 $\pi/64$ 弧度。離散化：
$$\text{index} = \frac{\alpha}{\pi/64} = (2n+1) \cdot k \cdot \frac{32}{N}$$
其中 $n$ = `pt`，$k$ = `freq`，$N$ = `size`，32 = `DCT_MAX`，以模 128 處理週期性。
* 填空 __ A06 __：將 $n$ = `pt` 代入 $(2n+1)$ 的 C 表示式。
* 填空 __ A07 __：離散化常數，將連續角度轉換為查詢表索引單位 (以巨集名稱表達)。
* 填空 __ A08 __：處理 $2\pi$ 週期性的模數 (以巨集名稱表達)。
* 二維 IDCT 以二次一維呼叫完成，利用可分離性將複雜度從 $O(N^4)$ 降為 $O(N^3)$。對 8×8 區塊，需 $2 \times 8^3 = 1024$ 次乘加。第一次 `scale = 10`，第二次 `scale = 10 + level`。查詢表以 $1024\sqrt{2} \approx 2^{10.5}$ 為縮放因子，二次乘法累積 $2^{21}$ 的定點縮放 (兩個 $\sqrt{2}$ 合成為 2，加上 $2^{10} \times 2^{10} = 2^{20}$)，額外的 `level` 位元右移補償 $N = 2^{\text{level}}$ 的歸一化。Parseval 定理保證量化誤差在頻域和空間域的能量相等。
- [ ] Part 3: 量化與預測
```c
/* amei.c:149-161 extract_block() 量化解碼 */
int sign = 1 - 2 * decode_bit(3);
int lv = (decode_unsigned(
25 + (is_chroma + (coeff < area / 8) * 2) * 10) + 1)
* (channel ? chroma : luma);
```
* `1 - 2 * decode_bit(3)` 將 $\{0, 1\}$ 映射為 $\{+1, -1\}$，無分支產生正負號。
* 上下文索引 `25 + (is_chroma + (coeff < area / 8) * 2) * 10` 對四種頻帶組合產生不同起始上下文，使低頻 (能量集中) 和高頻 (多為零) 各自維護獨立的機率模型。對亮度低頻 (`is_chroma = 0`, `coeff < area / 8`) 計算為 __ A09 __；對色度高頻 (`is_chroma = 1`, `coeff >= area / 8`) 計算為 __ A10 __。
* 影像標頭 `luma = 80`, `chroma = 160`，色度量化步長為亮度的 2 倍，利用人眼對色度敏感度較低的視覺特性。DC 預測修改 `buffer[0]`，`left && top` 時除以 2 取平均。方向性預測以 Q3 定點精度 (1/8 畫素) 線性內插，`weight = 4` 時結果為兩參考畫素的算術平均。
### 延伸問題
* `amei.c` 以 Chebyshev 遞推在軟體產生 128 項餘弦查詢表 (Q30 內部精度，Q10 輸出)。Linux 核心 V4L2 JPEG 硬體編解碼驅動程式採取截然不同的策略：[drivers/media/platform/samsung/s5p-jpeg/jpeg-core.c](https://github.com/torvalds/linux/blob/master/drivers/media/platform/samsung/s5p-jpeg/jpeg-core.c) 的 `s5p_jpeg_set_qtbl()` 將 ITU-T.81 Annex K 的 8×8 量化表逐 byte 寫入硬體暫存器，由硬體內部的定點管線執行 DCT/IDCT；[drivers/media/platform/renesas/rcar_jpu.c](https://github.com/torvalds/linux/blob/master/drivers/media/platform/renesas/rcar_jpu.c) 的 `jpu_set_qtbl()` 同樣載入 64 項預計算表 (`JPU_JPEG_QTBL_SIZE = 0x40`)。以此為出發點完成以下分析：
* 設計 C 程式以 `double` 精度餘弦值為參考，逐項比較 Chebyshev 遞推結果，統計 128 步後的最大偏差與 RMSE，驗證累積誤差是否仍在 Q10 的 $\pm 0.5$ 容許範圍內
* 從線性代數的角度，令 $\Psi$ 為 $N \times N$ DCT-III 矩陣，證明 `dct_coef[0] = 1024` 移除 DC 項的 $\sqrt{2}$ 因子後 $\Psi$ 成為正交矩陣。量化等效於將係數向量投射到格點 (lattice)，說明量化步長矩陣 $Q$ 與 $\Psi$ 的交互作用如何決定各頻率的重建誤差，以及 Parseval 定理為何保證此誤差在頻域和空間域能量相等
* 閱讀 `drivers/media/v4l2-core/v4l2-jpeg.c` 的 `v4l2_jpeg_parse_header()` 原始程式碼，說明 DQT marker 的量化表元素精度 ($Q_k$) 如何受 sample precision 約束：核心以 `pq` 欄位區分 8-bit 與 16-bit 量化表元素，但 16-bit $Q_k$ 僅在 extended sequential DCT 且 sample precision 為 12-bit 時合法；若使用者空間以 8-bit sample 搭配 16-bit $Q_k$，`v4l2_jpeg_parse_header()` 回傳 `-EINVAL`。分析此設計與 `amei.c` 固定 Q10 精度之間的取捨，並以 Parseval 定理說明：若硬體管線以截斷的 8-bit $Q_k$ 處理本應使用 16-bit $Q_k$ 的 12-bit 影像，頻域精度損失如何等量反映為空間域的 MSE 增量
* V4L2 test pattern generator `drivers/media/common/v4l2-tpg/v4l2-tpg-core.c` 以巨集 `COEFF(v, r)` 定義 BT.601 定點色彩轉換係數：`COEFF(0.299, 219)` 展開為 `(int)(0.5 + 0.299 * 219 * 256) = 16763`，隨後以 `>> 16` 還原，等效於 Q8 精度乘以 range 參數。`amei.c` 以 $1024\sqrt{2} \cdot 2^{20}$ 縮放餘弦值，屬 Q10 定點。以此為出發點：
* 計算 BT.601 最小係數 `0.114` 在 Q8 下的量化誤差 $|0.114 - 29/256|$，與 Q10 下的誤差 $|0.114 - \lfloor 0.114 \times 1024 + 0.5\rfloor / 1024|$ 比較
* V4L2 stateless H.264 解碼 API 的 `struct v4l2_ctrl_h264_slice_params` (定義於 `include/uapi/linux/v4l2-controls.h`) 中 `slice_alpha_c0_offset_div2` 為有號 `__s8`。依 ITU-T H.264 §7.4.3，此位元串流語法元素的合法範圍為 $[-6, +6]$，硬體以 `FilterOffsetA = slice_alpha_c0_offset_div2 * 2` 導出實際偏移 (範圍 [-12, +12])。設計一個無分支巨集，以 bitwise operation 先將未經驗證的 `__s8` 輸入 clamp 至 $[-6, +6]$ (提示：`max(-6, min(6, val))` 可用算術右移產生符號遮罩達成)，再以 `* 2` (而非 `<< 1`，因對負數左移為 C99 UB) 輸出 `FilterOffsetA`。說明為何核心以 `__s8` 而非 4-bit bit-field 儲存此欄位 (提示：uAPI 結構跨 user-kernel 邊界時的 ABI 對齊與 IOCTL 安全性)，以及若省略 clamp 直接將使用者空間傳入值寫入硬體暫存器，惡意輸入如何導致去方塊濾波器查表越界
* 說明為何核心在定點格式中偏好 2 的冪縮放因子，並從 `include/linux/log2.h` 的 `roundup_pow_of_two` 和 $\mathbb{Z}/128\mathbb{Z}$ 索引對映 `(i + 96) & 127` 的角度佐證
* V4L2 HEVC stateless API 的 `struct v4l2_ctrl_hevc_sps` 以 `log2_min_luma_coding_block_size_minus3` 和 `log2_diff_max_min_luma_coding_block_size` 兩個 `__u8` 欄位表達四元樹 CTU 分割範圍 (提示：`CtbSizeY = 1 << (field1 + 3 + field2)`，典型值為 64×64)。硬體從位元串流解析各 CTU 的分割決策，核心僅傳遞邊界參數。`amei.c` 的四元樹分割 (`level > 5` 強制分割、`level > 2` 條件分割) 是此機制的簡化版。以此為出發點：
* 從機率的角度，`amei.c` 的 83 組上下文模型中，四元樹分割決策使用上下文 0-2 (各對應 `level = 3, 4, 5`)。HEVC CABAC 使用 154 組上下文，其中 split_cu_flag 依 CTU 深度選擇上下文。對 `amei.c` 這類簡化 JPEG 式編碼器，假設亮度低頻 DCT 係數服從 Laplace 分布 $p(x) = \frac{1}{2b} e^{-|x|/b}$ 是合理的近似；推導 Laplace 參數 $b$ 與量化步長 `luma = 80` 的關係，並計算低頻頻帶的理論 bit rate。進一步說明此假設對 HEVC 為何不足：HEVC 的進階空間預測與多參考時間預測大幅壓縮殘差能量，使殘差分布呈高度尖峰 (leptokurtic)，更適合以 Generalized Gaussian Distribution (GGD，形狀參數 $\beta < 2$) 建模。對比 Laplace ($\beta = 1$) 與 GGD ($\beta \approx 0.5$) 在零值集中度上的差異，說明為何 CABAC 的 154 組上下文能比 `amei.c` 的 83 組更有效地追蹤此類稀疏分布
* 四元樹子區塊以 `corner % 2` / `corner / 2` 走訪，等效於 Z-order (Morton code) 位元交錯 $z = \ldots y_1 x_1 y_0 x_0$。設計 C 函式以 bitwise operation 實作 2D 座標到 Morton code 的轉換。在 `drivers/gpu/drm/` 中，DRM framebuffer modifier `DRM_FORMAT_MOD_*` 支援 tile-based 記憶體佈局 (如 Arm AFBC、Intel Y-tiled)；說明 tile 大小為何是 2 的冪，以及此佈局如何改善 GPU texture fetch 的 cache locality
---
## Problem B: 熵編碼、上下文模型與壓縮效率
- [ ] Part 1: 算術編碼器的適應動態
```c
/* amei.c:79-102 decode_bit() */
static int decode_bit(int context)
{
if (range < 256) {
range *= 256;
fraction *= 256;
if (inp < inp_end)
fraction += *inp++;
}
int *counts = models + context * 2;
int observations = *counts + counts[1] + 2,
split = range * (*counts + 1) / observations,
bit = fraction >= split;
if (bit) {
fraction -= split;
range -= split;
} else {
range = split;
}
counts[bit]++;
if (observations > 63) {
*counts /= 2;
counts[1] /= 2;
}
return bit;
}
```
* 重正規化 `range *= 256` 等效於左移 8 位元，為精度騰出空間。`models[166]` 提供 83 組上下文容器，每組 2 個計數器。
* Laplace rule of succession：`observations = c_0 + c_1 + 2` 等價於 Beta(1,1) 均勻先驗的後驗均值，$\hat{p}_0 = (c_0 + 1) / (c_0 + c_1 + 2)$。初始時 $c_0 = c_1 = 0$ 給出 $\hat{p}_0 = 1/2$，確保任何符號的機率不為零；若改用最大似然估計，未觀測的符號分配到 `range` 為 0 的子區間，永遠無法解碼。
* `observations > 63` 時兩計數器 `/= 2` (右移 1 位元)，等效於 EWMA 機制。假設某上下文已累積 `counts[0] = 50`, `counts[1] = 14` (`observations = 66 > 63`)，縮放後 `counts[0]` = __ B01 __，`counts[1]` = __ B02 __。新的 $\hat{p}_0 = 26/34 \approx 0.76$ (縮放前 $51/66 \approx 0.77$)，差異源自整數截斷。
* 連續半衰的有效記憶窗口：每輪半衰前最多累積 63 個觀測，第 $k$ 輪觀測在 $n$ 輪後權重 $(1/2)^{n-k}$。等效樣本數 $\sum_{j=0}^{\infty} 63 \cdot (1/2)^j = $ __ B03 __。
- [ ] Part 2: 指數哥倫布編碼與壓縮效率
```c
/* amei.c:104-113 decode_unsigned() */
static int decode_unsigned(int context)
{
int length = 0;
while (!decode_bit(context + length))
length++;
int value = 1;
while (length--)
value = value << 1 | decode_bit(4);
return value - 1;
}
```
編碼長度 $2\lfloor\log_2(v+1)\rfloor + 1$ 位元。前綴迴圈中 `length` 遞增使每個前綴長度使用不同上下文 (context, context+1, ...)；資料位元統一使用上下文 4。小數值用極少位元 ($v = 0$ 僅 1 位元)，天然適合 DCT 係數的稀疏分布。
* $v = 15$ 時 $v+1 = 16$ = `10000` (5 位元)，編碼長度 $2 \times 4 + 1 = $ __ B04 __ 位元。
完整上下文分配表 (83 組，參閱 `amei.md` §上下文模型)：
| 上下文索引 | 用途 |
|-----------|------|
| 0-2 | 四元樹分割決策 (`level = 3, 4, 5`) |
| 3 | DCT 係數正負號 |
| 4 | 無號整數資料位元 |
| 5-14 | 亮度零遊程字首 |
| 15-24 | 色度零遊程字首 |
| 25-34 | 高頻亮度量化等級 |
| 35-44 | 高頻色度量化等級 |
| 45-54 | 低頻亮度量化等級 |
| 55-64 | 低頻色度量化等級 |
| 65-72 | 結束旗標 (`61 + level*2 + is_chroma`，`level` = 2–5) |
| 73-82 | 預測模式 |
低頻與高頻以 `coeff < area / 8` 分界 (前 1/8 為低頻)，各自維護獨立機率模型，達成 adaptive bit allocation。
- [ ] Part 3: 四元樹率失真與 Amei 影像壓縮比
```c
/* amei.c:136-144 extract_block() */
if (level > 5 || (level > 2 && decode_bit(level - 3))) {
int size = 1 << --level;
for (int corner = 0; corner < 4; corner++)
extract_block(left + corner % 2 * size,
top + corner / 2 * size, level);
return;
}
```
四元樹分割：`level > 5` 強制分割；`level > 2` 條件分割 (讀取 1 位元)；`level == 2` 為最小末端節點。子區塊以 `corner % 2` / `corner / 2` 分離水平與垂直偏移，走訪順序類似 Z-order 曲線。128×128 影像 (`level = 7`) 最大遞迴深度 5 層；全部分割到 4×4 時有 1024 個末端節點。
* Amei 影像壓縮資料 776 位元組 = 6208 位元，未壓縮 128×128×3 影像 (每畫素 8 位元) 為 393216 位元，壓縮比約 __ B05 __ (取整至個位)。
* 二元訊源 Shannon 熵 $H = -p\log_2 p - (1-p)\log_2(1-p)$。當 $p = 0.9$ 時 $H \approx$ __ B06 __ (小數點後二位)。Huffman 編碼每符號至少 1 位元，算術編碼可逼近此分數碼長。四元樹每個節點僅 1 位元開銷，平坦區域以大區塊編碼、細節區域以小區塊編碼，逼近率失真曲線的更優操作點。
### 延伸問題
* `decode_unsigned()` 的 `value = value << 1 | decode_bit(4)` 以有號 `int` 逐位元累積，`length >= 31` 時觸發 C99 6.5.7 §4 UB (Problem 3 UB-2)。H.264 規格 (ITU-T H.264 §9.1) 的 CAVLC 路徑同樣使用 order-0 Exp-Golomb 編碼語法元素，Linux 核心 `lib/lzo/lzo1x_decompress_safe.c` 的位元建構迴圈則以 bounds check 防禦惡意輸入。以此為出發點：
* 推導 `length` 的最大安全值，設計修正方案 (將 `value` 改為 `unsigned int` 或加入上界檢查)，分析 `unsigned` 方案對 `return value - 1` 在 C99 6.3.1.3 下的隱式轉換語意
* 從機率的角度，假設 DCT 係數服從參數 $p$ 的 geometric 分布 $P(v) = (1-p)^v p$，推導 Shannon 熵 $H = \frac{-(1-p)\log_2(1-p) - p\log_2 p}{p}$，以 $p = 0.7$ 代入，計算指數哥倫布的期望碼長與 $H$ 的差距
* 閱讀 `lib/lzo/lzo1x_decompress_safe.c` 原始程式碼，找出 `NEED_OP()` / `NEED_IP()` 巨集如何在解壓縮迴圈中對輸入長度做邊界檢查，說明此防禦策略能否直接套用於 `decode_unsigned()` 的前綴迴圈
* 算術編碼的計數器半衰 (`/= 2`) 等效於 EWMA 機制 (參閱 `ewma.md`)。HEVC CABAC 使用不同的適應策略：每個上下文狀態 $\sigma \in [0, 63]$ 透過查表 `transIdxMPS[]` / `transIdxLPS[]` 以固定步長向等機率點靠近或遠離，構成 64 態的有限狀態機。以此為出發點：
* 建立數學模型：設 $c_0^{(n)}$ 為 `amei.c` 第 $n$ 次半衰後的 bit-0 計數，推導 $c_0^{(n)}$ 與原始累積計數 $C_0$ 的關係，並求等效 EWMA 衰減因子 $\alpha$
* 與 HEVC CABAC 的有限狀態機比較：在觀測序列 $0, 0, 0, \ldots, 1, 1, 1, \ldots$ (先連續出現 30 個 0 再切換至 1) 的情境下，兩種模型各需多少步才能使 $\hat{p}_1$ 超過 0.5？以此量化兩者的「遺忘速度」差異
* 從 Linux 核心的 PELT 機制出發 ，說明 PELT 的固定半衰期 32 ms 與此算術編碼器的 `observations > 63` 門檻在「以指數衰減折扣歷史」的設計哲學上有何共通之處，以及兩者在重尾負載 / 非平穩訊源下的侷限
* `amei.c` 的 CRC32 以 16 項 nibble 查表實作，核心 `lib/crc32.c` 使用 256 項 byte 查表，DRM/KMS 顯示子系統則依賴硬體 CRC。以此為出發點建立從代數到系統層級的完整分析：
* 從 $\text{GF}(2)[x]$ 的多項式代數出發，證明兩輪 4-bit 查表等效於一輪 8-bit 查表 (提示：模運算的結合律，$R_{8} = ((R \cdot x^4 \bmod G) \cdot x^4 \bmod G) = R \cdot x^8 \bmod G$)，並計算兩者的空間取捨 (64 bytes vs 1024 bytes)
* 閱讀 `drivers/gpu/drm/drm_debugfs_crc.c` 的 `drm_crtc_add_crc_entry()` 和 `drivers/gpu/drm/vkms/vkms_composer.c` 的 CRC 計算路徑，說明 VKMS 為何以 `crc32_le()` (Ethernet 多項式 `0xEDB88320`) 逐列累積而非逐幀一次計算，以及 i915 驅動程式的硬體 CRC 來源 (`plane1`、`pipe` 等) 為何可能使用與 Ethernet 不同的多項式。DRM CRC debugfs 介面 (`/sys/kernel/debug/dri/*/crc/`) 如何抽象化此差異？
* 從率失真理論的角度，高斯 i.i.d. 訊源的 $R(D) = \frac{1}{2}\log_2(\sigma^2 / D)$ bits/sample；以 $\sigma^2 = 2000$、$D = 50$ (MSE) 代入得約 2.66 bits/pixel。Amei 的實際 0.126 bits/pixel 遠低於此值，原因是 DCT + 預測大幅降低殘差方差。設計實驗：解碼 Amei 影像後，分別計算原始畫素方差與 IDCT 前殘差係數方差，以此估算有效 $\sigma_{\text{eff}}^2$ 並重新代入 $R(D)$ 公式驗證
---
## Problem C: 後處理管線與 C 語言合規性
- [ ] Part 1: 去方塊濾波與雙邊濾波
```c
/* amei.c:201-219 deblock() */
for (int bs = 4; bs <= 64; bs <<= 1)
for (int dir = 0; dir < 2; dir++) {
/* ... */
int a = p[i - step], b = p[i];
if (abs(a - b) <= 4)
continue;
p[i - step] = (3 * a + b + 2) >> 2;
p[i] = (a + 3 * b + 2) >> 2;
p[i - 2 * step] = (7 * p[i - 2 * step] + a + 4) >> 3;
p[i + step] = (7 * p[i + step] + b + 4) >> 3;
}
```
`abs(a - b) <= 4` 為平坦區域保護：差異小時視為自然紋理而跳過濾波，僅對差異大的方塊邊界施以平滑。邊界像素 `(3a + b + 2) >> 2` 以 3:1 權重混合 (25%)，次邊界 `(7p + a + 4) >> 3` 以 7:1 權重混合 (12.5%)，漸進衰減將影響限於邊界 ±2 像素。
```c
/* amei.c:229-257 bilateral3() */
int sw = (2 - abs(dx)) * (2 - abs(dy));
```
權重核 `(2 - |dx|)(2 - |dy|)`：中心 4、四鄰 2、對角 1。自適應門檻 `lt = threshold - grad / 6` (下限 `threshold / 2`) 使梯度越大時門檻越嚴格。亮度 `threshold = 8, grad_limit = 42`，色度 `threshold = 14, grad_limit = 56` (色度量化步長 160 為亮度 80 的 2 倍，雜訊更大)。實際分母 `wt` 動態累計，被排除的鄰域不參與平均，實現選擇性平滑而非全域模糊。
- [ ] Part 2: YCgCo 色彩轉換與 PNG 編碼
```c
/* amei.c:302-308 write_png() 色彩轉換 */
int t = image[0][i] - image[1][i];
U8AD(CLAMP8(t + image[2][i])); /* R */
U8AD(CLAMP8(image[0][i] + image[1][i])); /* G */
U8AD(CLAMP8(t - image[2][i])); /* B */
```
YCgCo 逆轉換 (`image[0]` = Y，`image[1]` = Cg，`image[2]` = Co)，僅涉及加減法。令 $t = Y - Cg$，展開可得 $R = Y - Cg + Co$，$G = Y + Cg$，$B = Y - Cg - Co$。相比 YCbCr 的乘法轉換矩陣，YCgCo 僅需加減法，適合定點實作。
PNG 寫入以巨集串流輸出。CRC32 以 16 項 nibble 查表 (每位元組 2 輪) 取代標準 256 項 byte 查表，兩輪 4-bit 查表因 CRC 的線性性代數上等效。Adler-32 使用模數 65521 (小於 $2^{16}$ 的最大質數)。zlib 標頭 `\x78\x01` 指定 deflate 壓縮 (最快等級)，實際 deflate 串流以 stored block 格式逐列輸出，完全不需實作 LZ77/Huffman。
- [ ] Part 3: 未定義行為與系統層級問題
`amei.c` 至少存在三處 C99 undefined behavior：
UB-1：十六進位解碼器。若輸入大寫 `'F'`，`'F' - 'a' + 10 = 70 - 97 + 10 = -17`。隨後 `(hi << 4) | lo` 對負數左移，C99 6.5.7 §4 明確將此列為 undefined behavior。
UB-2：算術編碼器的乘法溢位。`split = range * (*counts + 1) / observations` 中，`range` 重正規化後可接近 $2^{31}$，`*counts + 1` 可達 33，乘積溢位 `int`。修正方式是將 `range` 轉型為 __ C01 __ (C99 標準型別) 再相乘。
UB-3：`CLAMP8` 巨集的多次求值。
```c
#define CLAMP8(v) ((v) < 0 ? 0 : (v) > 255 ? 255 : (v))
```
對 `v` 最多求值 __ C02 __ 次 (當 `0 ≤ v ≤ 255`，三個 `(v)` 皆被求值)。若傳入 `CLAMP8(x++)`，`x` 可能遞增同樣次數，在兩個序列點間修改同一物件超過一次 (C99 6.5 §2)。修正方式是改用 `static inline` 函式。
記憶體配置：影像平面 `img_y`, `img_cg`, `img_co` 各 $2^{20}$ 個 `int` (各 4 MiB，合計 12 MiB)，`bilateral3()` 的 `static int tmp[1 << 20]` 佔 4 MiB，均存於 BSS 區段。總 BSS 空間 __ C03 __ MiB。若 `tmp` 改為非 `static`，4 MiB 配置在 stack 上，Linux 預設 stack 大小 8 MiB 下很可能 stack overflow。
### 延伸問題
* `amei.c` 的去方塊濾波 `(3a + b + 2) >> 2` 與 V4L2 stateless H.264 解碼的去方塊機制有直接對應。閱讀 `include/uapi/linux/v4l2-controls.h` 中 `struct v4l2_ctrl_h264_slice_params` 的 `disable_deblocking_filter_idc` (0=啟用, 1=停用, 2=slice 邊界停用)、`slice_alpha_c0_offset_div2` 和 `slice_beta_offset_div2` 三個欄位，以及 H.264 規格 (ITU-T H.264 §8.7) 的 52 組 alpha/beta 門檻表 (以 QP 為索引)。以此為出發點：
* `amei.c` 以固定門檻 `abs(a - b) <= 4` 決定是否濾波，H.264 以 `alpha[QP + FilterOffsetA]` 動態調整。推導：從浮點數截斷的角度，`(3a + b + 2) >> 2` 中 `+2` 提供四捨五入，若省略此偏移 (改為 `(3a + b) >> 2`)，截斷誤差的期望值從 0 變為多少？推廣至一般形式：定點加權平均 $(w_1 x + w_2 y + r) >> s$ 中，令 $r = 2^{s-1}$ 使期望截斷誤差為零的數學依據為何？
* UB-1 涉及負數左移 (C99 6.5.7 §4)。設計一個無分支的十六進位解碼函式：以 `| 0x20` 統一大小寫 (利用 `'A'` = 0x41 與 `'a'` = 0x61 僅差 bit 5)，說明此技巧對 `'0'`-`'9'` 為何無效 (提示：`'0' | 0x20 = 0x30`)。閱讀核心 `lib/hexdump.c` 的 `hex_to_bin()` 實作，比較其分支策略與你的無分支方案
* `CLAMP8` 巨集的多次求值問題 (UB-3) 與核心 `include/linux/minmax.h` 的 `clamp()` 巨集形成對比。設計測試程式，以 GCC `-O0` 和 `-O2` 分別編譯 `CLAMP8(x++)`，記錄 `x` 的遞增次數差異，從 C99 6.5 §2 (sequence point) 的角度解釋觀察結果
* YCgCo 逆轉換僅需加減法，而 BT.601 YCbCr 需要無理數係數 (如 0.299, 0.587, 0.114) 的定點逼近。Linux 核心 V4L2 以 `include/uapi/linux/videodev2.h` 的 `V4L2_YCBCR_ENC_601` / `V4L2_YCBCR_ENC_709` 列舉多種 YCbCr encoding，test pattern generator `drivers/media/common/v4l2-tpg/v4l2-tpg-core.c` 以 `COEFF(v, r)` 巨集實作定點轉換 (如 `COEFF(0.299, 219) = (int)(0.5 + 0.299 * 219 * 256) = 16763`，Q8 + range 縮放)。以此為出發點完成以下分析：
* 從線性代數的角度，計算 YCgCo 矩陣 $\begin{bmatrix} 1 & -1 & 1 \\ 1 & 1 & 0 \\ 1 & -1 & -1 \end{bmatrix}$ 與 BT.601 RGB→YCbCr 矩陣的條件數 $\kappa_2 = \sigma_{\max}/\sigma_{\min}$ (以奇異值分解)，說明條件數較小者在定點截斷下誤差放大倍率更低
* 閱讀 `v4l2-tpg-core.c` 的 `rgb2ycbcr()` 函式，追蹤 `COEFF` 巨集如何將 BT.601 係數轉為 Q8 定點，以及最終 `>> 16` 還原的完整計算路徑。計算最小係數 `0.114` 在 Q8 下的量化誤差 $|0.114 - 29/256|$ 與在 Q10 下的改善幅度
* Rockchip ISP1 的去雜訊預濾波器 (DPF，定義於 `include/uapi/linux/rkisp1-config.h`) 在結構上等同雙邊濾波：`struct rkisp1_cif_isp_dpf_rb_flt` 的 `spatial_coeff[6]` 對應空間域權重，`struct rkisp1_cif_isp_dpf_nll` 的 `coeff[17]` (支援線性或對數刻度) 將畫素強度對映至雜訊方差。此設計與 `amei.c` 的 `bilateral3()` 直接對應：`(2 - |dx|)(2 - |dy|)` 為空間核，`lt = threshold - grad / 6` 為強度域門檻。從機率的角度，雙邊濾波等效於以 Markov random field 為先驗的 MAP 估計。推導：若雜訊服從 $\mathcal{N}(0, \sigma^2)$，`threshold = 8` 對應約幾個標準差的接受範圍？閱讀 `drivers/media/platform/rockchip/rkisp1/rkisp1-params.c` 的 `rkisp1_dpf_config()`，說明驅動程式如何將使用者空間傳入的 `nll` 表寫入硬體暫存器
* BSS 區段的 16 MiB 靜態配置在 Linux 核心模組語境下受到限制。核心模組的 `.bss` 區段由 `kernel/module/main.c` 的模組載入器透過 `vmalloc` 配置，而使用者空間的 BSS 透過 mmap 延遲配置並支援 demand paging。以此為出發點：
* 閱讀 `kernel/module/main.c` 中模組記憶體配置的相關程式碼 (提示：搜尋 `module_alloc` 和 `MODULE_DATA` 相關定義)，說明 16 MiB BSS 的核心模組在 Arm64 上是否會觸發 `vmalloc` 區間耗盡或單次配置大小限制
* 對比使用者空間 BSS (demand paging，page fault 時才配置實體頁面) 與核心空間 BSS (載入時立即配置) 在記憶體壓力下的行為差異：若系統記憶體不足，使用者空間的 BSS 頁面可被回收 (swap out)，核心模組的 BSS 是否有相同的回收機制？
* `amei.c` 的 `bilateral3()` 使用 `static int tmp[1 << 20]` (4 MiB)。若此函式在核心模組中被多個行程同時呼叫，`static` 變數為所有呼叫者共享，會引發 data race。V4L2 裝置驅動程式 (如 `drivers/media/platform/` 下的各 codec) 如何避免此問題？從 `drivers/media/platform/rockchip/rkisp1/rkisp1-params.c` 的暫存緩衝區配置策略出發，說明核心驅動以 per-context 配置 (而非 `static`) 確保執行緒安全的慣用手法
## 參考題解
:::spoiler
* A01: ==1.99759==
* A02: ==30.5==
* A03: ==-32==
* A04: ==1446==
* A05: ==$B(1 + k)$==
* A06: ==`2*pt+1`== OR `pt*2+1` OR `1+pt*2` OR `1+2*pt`
* A07: ==`DCT_MAX`==
* A08: ==`DCT_MAX4`==
* A09: ==45==
* A10: ==35==
* B01: ==25==
* B02: ==7==
* B03: ==126==
* B04: ==9==
* B05: ==63==
* B06: ==0.47==
* C01: ==`unsigned int`==
* C02: ==3==
* C03: ==16==
:::
Sign in via Facebook
Sign in via X(Twitter)
Sign in via GitHub
Sign in via Dropbox
Sign in with Wallet
Connect another wallet