昨天從 GEMM 講起：同一句 C = A × B 有很多算法，其中一條路需要很高的精度，但消費級卡把 FP64 閹得很慘。
然後留下一個比較討厭的問題：
既然 4060 上 FP16 Tensor Core 幾乎跟 INT8 一樣快，為什麼不直接把 FP64 降成 FP16？為什麼還要切成 INT8？
今天會多一點數學。不是要證明論文，是希望看完以後，你可以自己拿一張紙，把一個數字拆開、再拼回去。
AdaptiveGEMM 真正的革新，其實不是又寫了一顆比較勤勞的 matmul。
是硬體潮流把「快」跟「準」拆成兩條路之後，數學只好跟著改：
硬體最會的：INT8 Tensor Core
你要的結果：看起來像 FP64 的 C = A × B
中間那座橋，一邊叫 Ozaki Scheme，一邊叫 中國餘數定理（CRT）。
一個負責切，一個負責拼。
過去很長一段時間，科學計算預設的語言是 FP64。Tesla、A100 這類卡，double 是一等公民。
然後深度學習來了。模型不需要 16 位小數都準，它要的是吞吐。於是硬體開始把算力堆到：
FP16 / BF16 / TF32 / INT8 Tensor Core
消費級 Ada 走得更極端。Day 4 量過：
FP32 : FP64 ≈ 64 : 1
INT8 Tensor Core 如果是 1
FP64 CUDA Core 大概只剩 1/138
所以現在的取捨不是「數學家喜歡哪種精度」，而是：
硬體把精確算力變貴、把低精度算力變便宜。
那數學能不能改口：先在便宜的單元上算，再把答案拼回昂貴的精度？
這就是今天的主題。
這是最直覺的捷徑。FP16 Tensor Core 在 4060 上幾乎跟 INT8 一樣快。那把 double 轉 half，乘完再轉回來，不就結束了？
問題在於 FP16 會丟掉你回不去的東西。
一個正規化的 FP16，有效小數大約 11 bits（10 位 mantissa + 1 位隱藏的 1）。兩個 FP16 相乘，再累加進 FP32，中間每一步都在 rounding。你得到的是「看起來差不多」的浮點數，不是「這個整數唯一對應那個結果」。
對神經網路來說，差不多通常就夠了。
對 GEMM 這種會一路累加 K 次的線性代數來說，rounding 會跟著 K 走，誤差邊界也會跟著走。
更關鍵的是：FP16 沒有提供一個乾淨的逆運算。
double → half → 乘完 → 再轉回 double
中間已經不是同一個數字了。你沒辦法說：
我現在手上這 11 bits
可以唯一拼回原來那 53 bits
因為資訊已經丟了。
所以如果目標是「對 cuBLAS FP64 的 max error 還能待在 1e-9、甚至 1e-15」，FP16 不是一條拼圖，它是一條單向路。
INT8 看起來更糟，只有 8 bits。但它有一個 FP16 沒有的好處：
整數的餘數可以拼回來。
這就是為什麼要繞去 INT8，而不是停在 FP16。
Ozaki Scheme 來自 2012 年的一篇論文：
Ozaki, Ogita, Oishi, Rump.
Error-free transformations of matrix multiplication...
Numerical Algorithms, 2012.
它要解決的事情很樸素：浮點數乘法本身就有 rounding。那能不能先把一個浮點數切成幾片，讓每一片的 bits 少到「兩片相乘是精確的」，再把這些精確的部分乘積加回去？
我們程式裡的切法，本質上就是：
scale = 2^shift
hi = round(v * scale) / scale
lo = v - hi
注意第二行：lo = v - hi。
在浮點運算裡，這個减法是 error-free transformation 的關鍵一步。hi 拿走了 v 比較高位的部分，lo 是剩下來的、比較小的尾巴。兩者加回去，就是原來的 v：
v = hi + lo
用一個具體的數字走一次。假設：
v = 3.1415926535
shift = 8
scale = 256
v * 256 = 803.4473203...
round(...) = 803
hi = 803 / 256 = 3.13671875
lo = 3.1415926535 - 3.13671875
= 0.0048739035
hi 變粗了，但 bits 變少了。lo 很細，但絕對值變小了。
矩陣乘法就變成：
C = A × B
= (A_hi + A_lo)(B_hi + B_lo)
= A_hi B_hi + A_hi B_lo + A_lo B_hi + A_lo B_lo
四個乘積，每一個都比原來那次「完整 FP64 GEMM」更容易被對應到低精度硬體。
我們後面會再乘上對應的 2^{-shift}，把整數結果縮放回浮點。
AdaptiveGEMM 裡，主路徑常用：
high scale sh = 2^17
low scale sl = 2^34
對應四個 pass 的縮放大約是：
hi × hi → 2^{-34}
hi × lo → 2^{-51}
lo × hi → 2^{-51}
lo × lo → 2^{-68}
切愈多片，能保住的 bits 愈多，但乘積數量也愈多。兩片就要 4 次 GEMM，三片就是 9 次。這是 Ozaki 自己的第一個取捨：
片數 s
愈多 → 愈接近 FP64
愈多 → GEMM 次數以 s² 成長
切完還只做了一半。因為就算 hi 已經比較瘦，乘上 2^17 之後，它仍然是個 遠大於 8 bits 的整數。
[-1, 1] 的數字 × 2^17
→ 整數大約落在 ±131072
→ 這已經不是 INT8 了
INT8 Tensor Core 一次只能吃 -128 ~ 127。
所以 Ozaki 負責「切成可以控制誤差的薄片」，但還沒辦法直接餵給 Tensor Core。
這時需要第二座橋。
先忘掉矩陣，只看一個整數。
假設我告訴你：
x ≡ 2 (mod 3) 除以 3 餘 2
x ≡ 3 (mod 5) 除以 5 餘 3
x ≡ 2 (mod 7) 除以 7 餘 2
3、5、7 兩兩互質。中國餘數定理說：在 0 到 3×5×7-1 之間，也就是 0 到 104 之間，剛好只有一個 x 同時滿足這三條。
那個數字是 23。我們一步一步拼出來。
M = 3 × 5 × 7 = 105
x 最後會落在 0 ~ 104，或等一下再折成有號整數。
對每個質數 p_i，先算：
M_i = M / p_i
也就是「把其他質數都乘起來」。
p=3 M_i = 105 / 3 = 35
p=5 M_i = 105 / 5 = 21
p=7 M_i = 105 / 7 = 15
35 可以被 5 跟 7 整除，所以對那兩個模數來說，這塊積木是 0。它只在 mod 3 時還活著。
要找 y_i，使得：
(M_i × y_i) ≡ 1 (mod p_i)
35 ≡ 2 (mod 3) 2 × 2 = 4 ≡ 1 (mod 3) → y = 2
21 ≡ 1 (mod 5) 1 × 1 = 1 → y = 1
15 ≡ 1 (mod 7) 1 × 1 = 1 → y = 1
於是每一塊「重建係數」是：
f_i = M_i × y_i
f_3 = 35 × 2 = 70
f_5 = 21 × 1 = 21
f_7 = 15 × 1 = 15
餘數是 a_i：
a_3 = 2
a_5 = 3
a_7 = 2
x ≡ a_3 f_3 + a_5 f_5 + a_7 f_7 (mod 105)
≡ 2×70 + 3×21 + 2×15
≡ 140 + 63 + 30
≡ 233
≡ 233 - 2×105
≡ 233 - 210
≡ 23
拼回來了。
這就是中國餘數定理在做的事：
一個太大、塞不進 INT8 的整數
↓
拆成好幾份「除以小質數的餘數」
↓
每份都很小，可以當 INT8
↓
算完以後，用 f_i 當積木加回去
↓
mod M，得到原來那個大整數
如果結果可能是負的，最後再對一下中心：
if (x > M/2) x = x - M
這樣 x 就落在大約 [-M/2, M/2]，而不是 [0, M)。
因為餘數對加法和乘法都成立：
(a + b) mod p = (a mod p + b mod p) mod p
(a × b) mod p = (a mod p × b mod p) mod p
所以內積也成立：
(Σ_k A[i,k] B[k,j]) mod p
=
(Σ_k (A[i,k] mod p) (B[k,j] mod p)) mod p
也就是說：
不要在大整數上做一次 GEMM
改成：
對每個質數 p，做一次 INT8 GEMM
得到 C_p = A_p × B_p
然後對每個輸出元素，用 CRT 拼回那個大整數
這就是 AdaptiveGEMM 的骨架。
Ozaki 先把 FP64 變成「bits 受控的整數薄片」。
CRT 再把這些還是太大的整數，變成 7 組塞得進 INT8 的餘數。
Tensor Core 真正去乘的，是餘數，不是 double。
全部都小於 128，這樣餘數才塞得進 int8：
p = 127, 113, 109, 107, 103, 101, 97
它們的乘積是：
M = 168897325606883
≈ 1.69 × 10^14
約 47 bits。這就是這組質數一次能唯一重建的整數範圍。
程式裡預先算好的重建係數 f_i，就是上面那個 M_i × y_i：
p=127 f = 147618922380819
p=113 f = 112099994871825
p=109 f = 134807957135769
p=107 f = 34726552928518
p=103 f = 96747011755399
p=101 f = 130435558389474
p=97 f = 19153304965729
拼的時候就是：
x = ( Σ (C_p[i,j] mod p) × f_p ) mod M
if (x > M/2) x = x - M
C[i,j] += x × 2^{-shift}
7 個質數不一定每次都用滿。我們會先掃每個 tile 的最大值，估計：
threshold ≈ K × max|A_slice| × max|B_slice|
然後看這個 threshold 掉進 M 的哪一段前綴，決定這次用幾個質數 nl = 1..7：
M_1 = 127
M_2 = 127×113 = 14351
M_3 = ... = 1564259
...
M_7 = 168897325606883
數字很小的 tile，3 個質數可能就夠了。K 很大、動態範圍很大的 tile，才需要 7 個。
這就是名字裡 Adaptive 的由來：不是永遠做最重的那條路，而是範圍要多大，就拼多大。
可以把它想成三組旋鈕。每一組都是硬體潮流逼出來的。
ss 愈大
愈接近 FP64
但 GEMM 次數是 s²
兩片 high/low，就是 4 個 pass。再切下去可以更準，也更貴。
消費級卡的時間預算很緊，所以主路徑停在 high/low，而不是切成十片。
nn 愈大
M 愈大，能重建的整數愈大
但 INT8 GEMM 要做 n 次，還要存 n 份餘數矩陣
7 個小於 128 的質數，是「塞進 INT8」跟「M 大約 47 bits」之間的折衷。
少兩個，範圍可能不夠撐 K=4096 的累加；多兩個，precompute 跟記憶體會再漲一截。
FP16
快
但 rounding 是單向的，拼不回唯一的大整數
INT8 + CRT
每次要做比較多顆 GEMM
但餘數是資訊，不是垃圾
所以可以對答案，對到 cuBLAS FP64
這就是昨天那個問題的答案。
INT8 不是因為 8 bits 比較準。
8 bits 很不准。准的是：8 bits 的餘數，配合互質的模數，可以唯一決定一個大整數。
FP16 沒有這件事。
中國餘數定理不是新的。Ozaki Scheme 也不是 2026 年才發明的。
新的是硬體把問題改寫了。
以前
硬體會 FP64
數學就用 FP64
現在（消費級 Ada）
硬體把 FP64 砍成 1/64
卻把 INT8 Tensor Core 堆得很肥
數學只好改成：
切成能控制的整數薄片 ← Ozaki
拆成能塞進 INT8 的餘數 ← CRT
用 Tensor Core 做 GEMM
再拼回看起來像 FP64 的結果
這不是「用比較髒的精度假裝有 FP64」。
這是承認：
你要的精度，跟硬體願意給你的吞吐，已經不在同一條路上。
那就把數字改造成硬體肯算的形狀，算完再變回來。
代價是多做幾次 GEMM、多存幾份餘數、多一次 precompute。
換到的是：在 4060 這種卡上，還能跟 cublasDgemm 對答案。
後面幾天會看到那些代價具體長什麼樣子。今天先把「為什麼這套數學會出現」講完。
可以帶回三句話：
明天把這套數學收成三個可以寫成 kernel 的步驟：建構 CRT、求餘數、把數字拼回來。先分開寫，對答案；再後面才談為什麼要融在一起。
今天如果你願意，用 3、5、7 把 23 再拼一次。拼得回來，後面的 kernel 就只是把這件事做成矩陣。
iThome鐵人賽
看影片追技術