Ozaki Scheme I 与中国剩余定理：从位切分到 Residue Number System

这篇是 Ozaki Scheme II：用 CRT 把高精度 GEMM 映射成少量 INT8 GEMM 的前置数学笔记。

建议阅读顺序：Ozaki I → CRT → Ozaki II。

0. 一句话先建立心智模型

Ozaki Scheme I 和 II 都在做同一件事：

不用慢的高精度乘法器直接做 GEMM，而是把高精度数换一种表示，再调用很多次极快、但可以保证“精确”的低精度 GEMM。

区别只在“表示方式”：

• Ozaki I：positional / 位段表示，把 significand 一段一段切开；
• Ozaki II：residue / 模数表示，把整数投影到多个互素模数下。

Scheme I 的表示天然产生 cross terms，所以 GEMM 数量近似二次增长；Scheme II 用中国剩余定理把表示改成多个独立的 residue channels，从数学结构上消掉 cross terms。

---

Part I：Ozaki Scheme I 到底是什么？

1. 从一个数开始：不是“降精度”，而是“把精度拆开保存”

假设一个十进制数

$$x=12.345678.$$

如果低精度硬件一次只能保留约 3 位有效数字，最粗暴的方法是

$$x\approx 12.3,$$

后面的信息直接丢掉。

Ozaki 的思想不是这样，而是写成

$$x=12.3+0.045+0.000678.$$

也就是

$$x=x_1+x_2+x_3.$$

每个 $x_i$ 只携带一段有效位，但是所有 slice 加起来仍然包含原数的高精度信息。

二进制中可以粗略想成：

FP64 significand

[ 最高 b bit ][ 接下来 b bit ][ 再接下来 b bit ] ...
      A1              A2              A3

如果每片保留 $b$ bit，则概念上

$$A\approx A_1+A_2+\cdots+A_k.$$

真正的 Ozaki splitting 比“直接切 raw mantissa”严谨：它会用 2 的幂做 shared-place scaling，把一组元素的有效位对齐到低精度表示可安全处理的位置。现代针对 Tensor Core / IMMU 的实现同样遵循“提取低精度矩阵、移动其位权、执行低精度 GEMM、再把位权移回来并高精度累加”的结构。

---

2. 为什么一定会出现 $A_iB_j$？

假设

$$A=A_1+A_2+A_3,$$ $$B=B_1+B_2+B_3.$$

那么

$$AB=(A_1+A_2+A_3)(B_1+B_2+B_3).$$

展开：

$$\begin{aligned} AB={}&A_1B_1\\ &+(A_1B_2+A_2B_1)\\ &+(A_1B_3+A_2B_2+A_3B_1)\\ &+(A_2B_3+A_3B_2)\\ &+A_3B_3. \end{aligned}$$

这不是 Ozaki 特有的麻烦，而是 positional decomposition 的代数必然结果。

十进制也一样：

$$1234=1200+34,$$ $$5678=5600+78.$$

于是

$$1234\times5678 =(1200+34)(5600+78),$$

必须出现

$$1200\times78,\qquad 34\times5600.$$

所以只算

$$A_1B_1+A_2B_2+A_3B_3$$

是不可能恢复 $AB$ 的。

---

3. 更数学一点：把 significand 看成 radix-$2^b$ 展开

为了看清“为什么是三角形”，把一个缩放后的高精度数写成

$$x\approx\sum_{r=0}^{k-1} d_r 2^{-br},$$

其中每个 digit $d_r$ 只有约 $b$ bit。

对应到矩阵：

$$A\approx\sum_{r=0}^{k-1}2^{-br}A_r,$$ $$B\approx\sum_{s=0}^{k-1}2^{-bs}B_s.$$

于是

$$AB \approx \sum_{r=0}^{k-1}\sum_{s=0}^{k-1} 2^{-b(r+s)}A_rB_s.$$

注意权重只取决于

$$r+s.$$

因此：

• $A_0B_0$：最高位；
• $A_0B_1,A_1B_0$：下一档；
• $A_0B_2,A_1B_1,A_2B_0$：再下一档；
• ……

矩阵乘积项自然排成反对角线：

                 B0        B1        B2        B3

A0             level 0   level 1   level 2   level 3
A1             level 1   level 2   level 3   level 4
A2             level 2   level 3   level 4   level 5
A3             level 3   level 4   level 5   level 6

所以如果只需要前 $k$ 档有效信息，就只算

$$r+s\le k-1.$$

这就是常见的三角区域。

需要的低精度 GEMM 数：

$$1+2+\cdots+k =\frac{k(k+1)}2.$$

例如 $k=8$：

$$\frac{8\cdot9}{2}=36.$$

这正是 Scheme II 想消灭的结构性开销。

---

4. Ozaki I 最重要的问题：低精度 GEMM 为什么能“精确”？

这是 Ozaki 的核心，而不是“切片”本身。

设一个 slice 中的整数 significand 只有 $b$ bit，粗略有

$$|a|,|b|<2^b.$$

单次乘积需要大约

$$2b$$

bit。

矩阵乘法一个输出元素是长度为 $q$ 的 dot product：

$$c_{ij}=\sum_{\ell=1}^{q}a_{i\ell}b_{\ell j}.$$

最坏情况下，多加 $q$ 项还需要约

$$\lceil\log_2q\rceil$$

bit。

因此只要 accumulator 能精确容纳

$$2b+\lceil\log_2q\rceil$$

bit，就可以让整个低精度 dot product 不发生舍入或整数溢出。

可以记成一个非常实用的心智公式：

$$\boxed{ 2b+\lceil\log_2q\rceil\lesssim p_{\mathrm{acc}} }$$

其中 $p_{\mathrm{acc}}$ 是 accumulator 的有效整数位宽/精确位数。

INT8 Tensor Core 的例子

若用 signed INT8：

$$|a|,|b|\le128\approx2^7,$$

所以

$$|ab|\lesssim2^{14}.$$

使用 INT32 accumulator，希望

$$q2^{14}<2^{31}.$$

得到

$$q<2^{17}.$$

在这个数量级以内，INT8 × INT8 → INT32 的 dot product 可以被当成一个 exact integer arithmetic primitive。

这说明：

Tensor Core 的输入位宽低，不等于最终算法一定“低数值精度”。Ozaki 是主动把每一片限制到硬件可以精确计算的范围，再通过多片重建高精度结果。

---

5. shared-place / scaling 到底在干什么？

真实 FP64 元素不只有 significand，还有 exponent。

例如同一行可能同时出现

$$2^{30},\qquad 2^{-20}.$$

如果直接按相同 absolute bit position 切，两个数根本对不齐。

因此 Ozaki 会选择 2 的幂作为尺度，把一组元素移动到统一的 fixed-point 位置。概念上可以写成

$$A_r=D_A^{-1}\,2^{-br}\,\widehat A_r,$$

其中 $D_A$ 是由 2 的幂组成的对角/共享尺度，$\widehat A_r$ 则是适合低精度 GEMM 的小整数/低精度 significand。

类似地：

$$B_s=2^{-bs}\,\widehat B_s\,D_B^{-1}.$$

于是

$$A_rB_s = D_A^{-1} 2^{-b(r+s)} (\widehat A_r\widehat B_s) D_B^{-1}.$$

真正交给高速低精度单元的是

$$\widehat A_r\widehat B_s,$$

算完后再把

$$D_A^{-1},\quad D_B^{-1},\quad 2^{-b(r+s)}$$

对应的位权恢复回来。

为什么特别爱用 $2^e$？因为二进制浮点乘 $2^e$ 本质上是在调整 exponent；在不发生 overflow/underflow 的正常范围内，它不会像一般实数乘法那样额外破坏 significand。

---

6. 所以 Ozaki I 的完整流水线是什么？

现代 GEMM-oriented Ozaki 可以压缩成：

Step 1：splitting

把高精度矩阵拆成

$$A\approx A_1+\cdots+A_k,$$ $$B\approx B_1+\cdots+B_k,$$

并通过 shared-place scaling 保证每个 slice 能安全进入低精度算术。

Step 2：低精度 GEMM

计算重要的

$$A_iB_j,\qquad i+j\le k+1.$$

每个低精度 product 被设计成尽可能 exact。

Step 3：恢复位权

把每个 $A_iB_j$ 按对应的 power-of-two scale 移回原来的数值尺度。

Step 4：高精度累加

$$D=\sum_{i+j\le k+1}A_iB_j.$$

若需要模拟完整 GEMM 接口，再做

$$C\leftarrow\alpha D+\beta C.$$

Ozaki I 的关键优势是：绝大部分计算都变成高吞吐 GEMM；关键缺点是：slice 数增加时，$A_iB_j$ 数量近似二次增长。

---

Part II：中国剩余定理到底说了什么？

7. 先从“余数能不能表示一个整数”开始

选择两两互素的整数

$$m_1,m_2,\dots,m_s,$$

满足

$$\gcd(m_i,m_j)=1\qquad(i\ne j).$$

定义

$$M=\prod_{i=1}^{s}m_i.$$

对于整数 $x$，记录它在每个模数下的余数：

$$r_i=x\bmod m_i.$$

也就是把一个整数表示成向量

$$\boxed{ x\longmapsto(r_1,r_2,\dots,r_s) }$$

这就是 Residue Number System 的核心表示。

中国剩余定理说：

当 $m_i$ 两两互素时，任意一组余数 $(r_1,\dots,r_s)$ 都对应唯一的 $x\pmod M$。

注意关键词是：

$$\boxed{\text{唯一 modulo }M}$$

不是“在所有整数中唯一”。

---

8. 为什么是唯一的？

假设两个整数 $x,y$ 在所有模数下余数都一样：

$$x\equiv y\pmod{m_i},\qquad\forall i.$$

那么

$$m_i\mid(x-y),\qquad\forall i.$$

因为所有 $m_i$ 两两互素，所以它们的乘积也整除 $x-y$：

$$M=m_1m_2\cdots m_s\mid(x-y).$$

因此

$$x\equiv y\pmod M.$$

所以在区间

$$0\le x<M$$

里只有一个可能的整数。

如果处理有符号整数，更常用 symmetric representative：

$$-\frac M2\le x<\frac M2.$$

只要真实答案落在这个区间里，也只有一个可能值。

---

9. CRT 的构造公式是怎么来的？

已知

$$x\equiv r_i\pmod{m_i},\qquad i=1,\dots,s.$$

定义

$$M_i=\frac{M}{m_i}.$$

因为 $m_i$ 与其他模数都互素，所以

$$\gcd(M_i,m_i)=1.$$

于是 $M_i$ 在模 $m_i$ 下存在乘法逆元。记

$$y_i=M_i^{-1}\pmod{m_i},$$

也就是说

$$M_i y_i\equiv1\pmod{m_i}.$$

构造

$$e_i=M_i y_i.$$

这个 $e_i$ 有一个非常漂亮的性质：

对自己的模数

$$e_i\equiv1\pmod{m_i}.$$

对其他模数 $m_j$，$j\ne i$

因为 $M_i$ 包含因子 $m_j$，所以

$$e_i=M_i y_i\equiv0\pmod{m_j}.$$

因此 $e_i$ 就像一个“模数空间里的 one-hot basis”：

$$e_i\equiv \begin{cases} 1\pmod{m_i},\\ 0\pmod{m_j},\quad j\ne i. \end{cases}$$

所以直接令

$$\boxed{ x=\sum_{i=1}^{s}r_iM_i y_i\pmod M }$$

就行。

因为模 $m_j$ 看过去，除了第 $j$ 项，其他项全部变成 0：

$$x \equiv r_jM_jy_j \equiv r_j \pmod{m_j}.$$

这就是 CRT reconstruction 公式的来源。

---

10. 一个真正算一遍的 CRT 例子

取

$$m_1=5,\qquad m_2=7.$$

于是

$$M=35.$$

假设我们知道

$$x\equiv2\pmod5,$$ $$x\equiv5\pmod7.$$

先算：

$$M_1=35/5=7.$$

要找 $7$ 在模 $5$ 下的逆元：

$$7\equiv2\pmod5.$$

而

$$2\times3=6\equiv1\pmod5,$$

所以

$$y_1=3.$$

再算

$$M_2=35/7=5.$$

$5$ 在模 $7$ 下的逆元也是 $3$，因为

$$5\times3=15\equiv1\pmod7.$$

所以

$$\begin{aligned} x &\equiv r_1M_1y_1+r_2M_2y_2\pmod{35}\\ &=2\times7\times3+5\times5\times3\\ &=42+75\\ &=117\\ &\equiv12\pmod{35}. \end{aligned}$$

检查：

$$12\bmod5=2,$$ $$12\bmod7=5.$$

完全一致。

---

Part III：为什么 CRT 特别适合矩阵乘法？

11. 模运算是一个 ring homomorphism

这是从 Ozaki I 跳到 Scheme II 的真正数学核心。

对任意整数 $a,b$：

$$(a+b)\bmod m = ((a\bmod m)+(b\bmod m))\bmod m,$$

以及

$$(ab)\bmod m = ((a\bmod m)(b\bmod m))\bmod m.$$

也就是说，映射

$$\phi_m:\mathbb Z\to\mathbb Z/m\mathbb Z$$

同时保持加法和乘法结构。

抽象一点：

$$\phi_m(a+b)=\phi_m(a)+\phi_m(b),$$ $$\phi_m(ab)=\phi_m(a)\phi_m(b).$$

这就是 ring homomorphism。

而矩阵乘法本质上只有“乘法 + 加法”：

$$(AB)_{ij}=\sum_{k}a_{ik}b_{kj}.$$

因此：

$$\begin{aligned} (AB)_{ij}\bmod m &= \left(\sum_k a_{ik}b_{kj}\right)\bmod m\\ &= \sum_k \left((a_{ik}\bmod m)(b_{kj}\bmod m)\right) \bmod m. \end{aligned}$$

所以整个矩阵层面：

$$\boxed{ AB\bmod m = (A\bmod m)(B\bmod m)\bmod m }$$

这就是 Ozaki II 可以直接“一模数一 GEMM”的数学原因。

---

12. 为什么 CRT 表示没有 $A_iB_j$ cross terms？

现在把两种表示并排看。

Ozaki I：positional representation

$$A=\sum_i2^{e_i}A_i,$$ $$B=\sum_j2^{f_j}B_j.$$

乘起来：

$$AB = \sum_i\sum_j 2^{e_i+f_j}A_iB_j.$$

必然产生笛卡尔积。

Ozaki II：residue representation

$$A\longmapsto (A\bmod m_1,\dots,A\bmod m_s),$$ $$B\longmapsto (B\bmod m_1,\dots,B\bmod m_s).$$

第 $i$ 个 channel 里直接得到：

$$C_i =(A\bmod m_i)(B\bmod m_i)\bmod m_i =(AB)\bmod m_i.$$

每个 channel 已经包含完整乘积 $AB$ 在该模数下的信息。

根本不存在

$$(A\bmod m_i)(B\bmod m_j),\qquad i\ne j$$

这种东西，因为不同 residue channel 属于不同的 ring；CRT reconstruction 只在最后把同一个整数的各个 residue 拼回来。

所以：

$$\boxed{ \text{positional digits：乘法会产生 cross terms} }$$ $$\boxed{ \text{residue channels：乘法逐 channel 封闭} }$$

这就是 Scheme II 把低精度 GEMM 数从近似 $O(k^2)$ 降成 $O(s)$ 的本质，而不仅仅是“用了一个 CRT 技巧”。

---

13. 但 CRT 只恢复 modulo $M$：怎么得到真实乘积？

这是 Scheme II 最容易漏掉的条件。

CRT reconstruct 得到的是

$$X\equiv A'B'\pmod M.$$

如果

$$X=42,$$

那么

$$42,\ 42+M,\ 42-M,\dots$$

都有相同 residue。

所以必须提前知道真实结果范围。

若能保证每个整数矩阵乘积元素满足

$$|(A'B')_{ij}|\le c_{\max},$$

并让

$$\boxed{2c_{\max}<M},$$

那么真实答案一定落在

$$(-M/2,M/2)$$

里，而这个区间中每个 residue class 只出现一次。

于是 symmetric CRT representative 就是真正的整数乘积，而不是仅仅一个 modulo-$M$ 的代表元。

一个保守的范围估计

若

$$|A'_{ik}|<2^{k_A},$$ $$|B'_{kj}|<2^{k_B},$$

内积长度为 $q$，那么

$$|(A'B')_{ij}| \le \sum_{k=1}^{q}|A'_{ik}||B'_{kj}| <q2^{k_A+k_B}.$$

因此可以要求

$$q2^{k_A+k_B}<M/2.$$

也就是

$$\boxed{ k_A+k_B < \log_2M-\log_2q-1 }.$$

这条式子解释了 Scheme II 的精度控制：

$$M=\prod_i m_i.$$

多加一个模数 $m_{s+1}$，就让

$$M\leftarrow Mm_{s+1},$$

可表示整数的 dynamic range 变大，于是可以让 $k_A,k_B$ 更大，保留更多输入 significand bit。

而代价只增加一个 residue GEMM。

---

14. 从代数角度看 CRT，Scheme II 会更漂亮

中国剩余定理其实给出了一个 ring isomorphism：

$$\boxed{ \mathbb Z/M\mathbb Z \cong \prod_{i=1}^{s}\mathbb Z/m_i\mathbb Z }$$

其中

$$M=m_1m_2\cdots m_s,$$

且 $m_i$ 两两互素。

左边是一个“大整数模 $M$”的 ring。

右边是 $s$ 个“小整数 ring”的直积。

因此一个大整数运算：

$$A'B'\pmod M$$

可以完全等价地变成 $s$ 个互不依赖的小整数运算：

$$A'_iB'_i\pmod{m_i},\qquad i=1,\dots,s.$$

这句话就是 Ozaki Scheme II 的数学抽象：

把一个宽整数 GEMM 同构成多个窄整数 GEMM。

GPU 上恰好有极快的 INT8 Tensor Core，所以把 $m_i$ 控制到 INT8 能装下的范围，就可以把每个小 ring 的矩阵乘法交给 Tensor Core。

---

Part IV：Scheme I → II 应该怎么记？

15. 两个算法其实是在解决同一个“位宽不匹配”问题

我们想计算 FP64 GEMM，但硬件最强的是 INT8/FP16 Tensor Core。

Scheme I 的办法

把一个宽 significand 表成多个 positional slices：

$$x=\sum_i d_i2^{e_i}.$$

优点：直观，和浮点有效位天然对应。

缺点：

$$xy = \sum_i\sum_jd_i e_j2^{e_i+f_j},$$

所以乘法产生 cross terms。

Scheme II 的办法

先整数化，再用 residue vector 表示宽整数：

$$x\longleftrightarrow (x\bmod m_1,\dots,x\bmod m_s).$$

优点：模映射保持乘法结构，所以：

$$xy\bmod m_i =(x\bmod m_i)(y\bmod m_i)\bmod m_i.$$

每个 channel 独立，完全没有 cross terms。

---

16. 最终心智图

                    “高精度数装不进最快的低精度 GEMM”
                                  │
                  ┌───────────────┴───────────────┐
                  │                               │
             Ozaki Scheme I                 Ozaki Scheme II
                  │                               │
          positional splitting              residue splitting
                  │                               │
        A = Σ 2^ei Ai                     A -> (A mod m1, ...)
        B = Σ 2^fj Bj                     B -> (B mod m1, ...)
                  │                               │
      AB = Σi Σj 2^(ei+fj) AiBj          AB mod mi = Ai Bi mod mi
                  │                               │
            cross terms                      channels independent
                  │                               │
              O(k²) GEMM                       O(k) GEMM
                                                  │
                                                  ↓
                                                  CRT
                                                  │
                                                  ↓
                                             wide integer

---

17. 面试时可以怎么讲

如果对方问 Ozaki，可以压缩成下面这段：

Ozaki Scheme I 用 error-free/shared-place splitting 把高精度矩阵的 significand 拆成多个低位宽 slice，使每个 slice pair 的低精度 dot product 能在底层 accumulator 中精确执行，再以高精度恢复 scale 并累加。问题是 positional decomposition 的乘法天然产生 $A_iB_j$ cross terms，所以为了逐步恢复更多有效位，需要计算三角形数量的 GEMM，复杂度约 $k(k+1)/2$。Scheme II 改用 CRT residue representation。因为 reduction modulo $m_i$ 是 ring homomorphism，$(AB)\bmod m_i=(A\bmod m_i)(B\bmod m_i)\bmod m_i$，所以每个 modulus 只需要一个独立 GEMM，最后 CRT 重建宽整数；只要总模数 $M$ 大于真实结果动态范围的两倍，就能从 modulo-$M$ 结果唯一恢复有符号真实乘积。

如果这段能自己推出来，Ozaki I 和 II 的核心就基本掌握了。

---

References

• K. Ozaki, T. Ogita, S. Oishi, S. M. Rump. Error-free transformations of matrix multiplication by using fast routines of matrix multiplication and its applications. Numerical Algorithms 59(1), 95–118, 2012. DOI: 10.1007/s11075-011-9478-1.
• K. Ozaki, T. Ogita, S. Oishi, S. M. Rump. Generalization of error-free transformation for matrix multiplication and its application. NOLTA, IEICE 4(1), 2–11, 2013. DOI: 10.1587/nolta.4.2.
• H. Ootomo, K. Ozaki, R. Yokota. DGEMM on integer matrix multiplication unit. International Journal of High Performance Computing Applications 38(4), 297–313, 2024. DOI: 10.1177/10943420241239588.
• Y. Uchino, K. Ozaki, T. Imamura. Performance enhancement of the Ozaki Scheme on integer matrix multiplication unit. International Journal of High Performance Computing Applications 39(3), 2025. DOI: 10.1177/10943420241313064.
• K. Ozaki, Y. Uchino, T. Imamura. Ozaki Scheme II: A GEMM-oriented emulation of floating-point matrix multiplication using an integer modular technique. arXiv:2504.08009; journal version first published online July 28, 2026.