4.1 旋转矩阵与欧拉角
上一章的 MVP 链只用到平移、缩放与透视投影，顶点一旦要「转身」就绕不开旋转。旋转在整数世界里比想象中难：config/ 禁浮点，而 sin / cos 天生是超越函数。本节先用纯整数造一张定点正弦表，再用它搭出绕 x / y / z 三轴的旋转矩阵，并断言两条不变量——sin² + cos² = ONE、旋转矩阵正交。如果你更关心四元数如何驱动蒙皮、在骨骼层级里省下逐骨骼的矩阵乘法，可以顺路看 骨骼动画里的四元数
；本节反过来，把旋转压成一组可断言的整数，每一个 >> 16 的舍入都要经得起复算。
为什么不能直接调 sin
浮点 sin 有两个问题。其一是口径：config/ 明确禁浮点，本卷所有数值走 Q16.16，且必须跨语言逐位可复现——Python 的 int 与 TypeScript 的 BigInt 算同一串整数要得到同一串结果，而 libm 的 sin 在不同平台连最后一位都不保证相同。其二是粒度：我们真正需要的不是「任意实数的 sin」，而是「二进制角度（BAM）的 sin」——把一整圈 360° 定义为 ONE = 65536 个单位，角度本身就是定点数，加减、折象限全用位运算，连弧度换算都省了。
于是问题变成：怎样在 Q16.16 上、只用整数运算，造出一张足够准的正弦表。
二进制角度：一圈就是 ONE
把整圈定义为 ONE 的做法叫二进制角度（BAM，binary angular measurement）。一圈 = 65536 个单位，于是 90° = 16384（ONE >> 2）、180° = 32768、270° = 49152，单位是 2⁻¹⁶ 圈。它带来两个直接好处。第一，角度加减就是整数加减，越界靠 & (ONE - 1) 折回，不需要 while 循环或取模分支。第二，象限折叠退化成位运算：a >> 14 给出象限号 0~3，a & (QUARTER - 1) 给出象限内的余角，两个操作就把任意角折进第一象限，再靠「正弦 ↔ 余弦互换、符号翻转」还原。这也是为什么本卷的 sin_bam / cos_bam 里没有一处 if 是数值比较——全是位运算。
代价是角度的分辨率固定为 360° / 65536 ≈ 0.0055°，对光栅化绰绰有余；真要更细，把 ONE 左移即可，本卷不需要。
定点正弦表：嵌套泰勒 + 查表 + 线性插值
分三步。第一步，只在 [0, 90°] 上直接算 sin，用嵌套泰勒级数（Horner 形式）截断到 u⁸/9!：
sin(u) ≈ u · (1 - u²/6 · (1 - u²/20 · (1 - u²/42 · (1 - u²/72))))
分母 6 = 3!、20 = 5!/3!、42 = 7!/5!、72 = 9!/7!，逐层用 fdiv 得到 Q16.16 的商，全程没有一次浮点，也没有一次对 math.sin 的调用。第二步，在 [0, 90°] 上取 257 个等距样本（N_SEG = 256）存成四分之一波表；sin 的其余三个象限靠位运算折回——sin(90°+r) = cos(r)、sin(180°+r) = -sin(r)、sin(270°+r) = -cos(r)。第三步，落在样本之间的角度用线性插值：表步长是 64 = 2⁶，于是 idx = t >> 6、frac = t & 63，而 frac / 64 的 Q16.16 表示恰好是 frac << 10，连一次除法都不需要。
反三角（acos / asin）用二分查表：cos 在 [0, 90°] 上单调递减，二分找最大的角度使 cos(a) ≥ x 即可。整段实现如下（工作区文件 lantern/sim/cg_raster_04_trig.py）：
#!/usr/bin/env python3
# file: lantern/sim/cg_raster_04_trig.py
# tier: S
# verify: /opt/homebrew/bin/python3 sim/cg_raster_04_1_rotation.py
"""定点三角函数：全整数构造，无浮点、无 math、无第三方库。
`config/` 禁浮点，但旋转需要 sin / cos。本模块的**全整数**方案：
1. 用**嵌套泰勒级数**在 Q16.16 定点上直接构造四分之一波正弦表（不调用 math.sin）；
2. 查表 + **线性插值**得到任意角度的 sin / cos；
3. 角度用**二进制角度**（BAM）：整个圆周 = ONE = 65536 个单位，即一圈 360°；
4. 反向的 acos 用**二分查表**得到，同样全整数。
与 TypeScript 侧 core/test/cg_raster_04_orientation.test.ts 的表构造逐字对齐
（两侧都「先乘后 >> 16」），故表内容逐位相同。
"""
import os
import sys
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
from cglib import ONE, fdiv, fmul # noqa: E402
TURN = ONE # 360° = ONE（二进制角度）
QUARTER = ONE >> 2 # 16384 = 90°
N_SEG = 256 # 四分之一波表段数（257 个条目）
_SHIFT = 8 # N_SEG = 2**_SHIFT；表步长 = QUARTER / N_SEG = 64 = 2^6
HALF_PI = 102944 # π/2 的 Q16.16 整数近似（round(π·65536/2)）
def sin_taylor_q16(u):
"""u ∈ [0, HALF_PI]（Q16.16 弧度），返回 sin(u) 的 Q16.16。
嵌套泰勒（截断到 u⁸/9!）：
sin(u) ≈ u·(1 - u²/6·(1 - u²/20·(1 - u²/42·(1 - u²/72))))
分母 6=3!、20=5!/3!、42=7!/5!、72=9!/7!；u²/72 等用 fdiv 得到 Q16.16 的商。
"""
u2 = fmul(u, u)
t = ONE - fdiv(u2, 72 * ONE)
t = ONE - fmul(fdiv(u2, 42 * ONE), t)
t = ONE - fmul(fdiv(u2, 20 * ONE), t)
t = ONE - fmul(fdiv(u2, 6 * ONE), t)
return fmul(u, t)
def build_sin_table():
"""四分之一波正弦表：257 个 Q16.16 条目，θ_k = k·(π/2)/256。
步长用整数右移取整 (k·HALF_PI) >> _SHIFT，保证 θ_256 恰好 = HALF_PI。
"""
return [sin_taylor_q16((k * HALF_PI) >> _SHIFT) for k in range(N_SEG + 1)]
SIN_TABLE = build_sin_table()
def _interp_quarter(t):
"""t ∈ [0, QUARTER]，返回该象限正弦的 Q16.16（线性插值）。
表步长 64 = 2^6，故 idx = t >> 6、frac = t & 63，frac/64 的 Q16.16 = frac << 10。
"""
idx = t >> 6
frac = t & 0x3F
s0 = SIN_TABLE[idx]
if frac == 0 or idx >= N_SEG:
return s0
s1 = SIN_TABLE[idx + 1]
return s0 + fmul(s1 - s0, frac << 10)
def sin_bam(a):
"""a 为二进制角度，返回 sin(a) 的 Q16.16。象限折叠全用位运算。"""
a &= ONE - 1
q = a >> 14 # 象限 0..3（每象限 16384 = 2^14）
r = a & (QUARTER - 1)
if q == 0:
return _interp_quarter(r)
if q == 1:
return _interp_quarter(QUARTER - r) # sin(90°+r) = cos(r)
if q == 2:
return -_interp_quarter(r) # sin(180°+r) = -sin(r)
return -_interp_quarter(QUARTER - r) # sin(270°+r) = -cos(r)
def cos_bam(a):
"""cos(a) = sin(a + 90°)。"""
return sin_bam((a + QUARTER) & (ONE - 1))
def acos_bam(x):
"""x ∈ [0, ONE] -> acos(x) 的二进制角度 ∈ [0, QUARTER]，二分查表。"""
if x >= ONE:
return 0
if x <= 0:
return QUARTER
lo, hi = 0, QUARTER
while lo < hi: # cos 在 [0, QUARTER] 单调递减，取最大的 a 使 cos(a) >= x
mid = (lo + hi + 1) >> 1
if cos_bam(mid) >= x:
lo = mid
else:
hi = mid - 1
return lo
def asin_bam(x):
"""x ∈ [0, ONE] -> asin(x) 的二进制角度 ∈ [0, QUARTER]，二分查表。
sin 在 [0, QUARTER] 单调递增，且近零处 sin(a) ≈ a，故小角度分辨率约 1 BAM。
"""
if x >= ONE:
return QUARTER
if x <= 0:
return 0
lo, hi = 0, QUARTER
while lo < hi: # 取最大的 a 使 sin(a) <= x
mid = (lo + hi + 1) >> 1
if sin_bam(mid) <= x:
lo = mid
else:
hi = mid - 1
return lo
这里有两个坑值得写下来。其一，表步长决定线性插值的误差：N_SEG 取 64 时 sin²+cos² 的最大偏差实测 14 LSB，超出容差；加密到 256 后降到 6 LSB。其二，表构造里那步右移 (k · HALF_PI) >> _SHIFT 的位数必须与查表时的 >> 6 严格匹配，两者都收敛到常量 _SHIFT = 8；一旦只改一处，就会算出「表最大误差 125631 LSB」这种一眼假的数字——那是插值步长与索引位数不一致导致的错位，不是精度问题。
三轴旋转矩阵
有了 sin_bam / cos_bam，绕坐标轴的旋转矩阵就是课本上的标准形式，只是把 cos / sin 换成定点调用：
# snippet: 4.1 三轴旋转矩阵（完整文件 sim/cg_raster_04_1_rotation.py）
def rot_x(a):
c, s = trig.cos_bam(a), trig.sin_bam(a)
return (ONE, 0, 0, 0, 0, c, -s, 0, 0, s, c, 0, 0, 0, 0, ONE)
def rot_y(a):
c, s = trig.cos_bam(a), trig.sin_bam(a)
return (c, 0, s, 0, 0, ONE, 0, 0, -s, 0, c, 0, 0, 0, 0, ONE)
def rot_z(a):
c, s = trig.cos_bam(a), trig.sin_bam(a)
return (c, -s, 0, 0, s, c, 0, 0, 0, 0, ONE, 0, 0, 0, 0, ONE)
s 与 -s 的符号位置决定了旋转的正方向（右手系、逆时针为正）。矩阵是行主序的 4×4，最后一行与最后一列保持单位，因为绕轴的纯旋转不动原点。
定点舍入：>> 与 // 不是一回事
本节所有降精度都用 >> 16（算术右移），它在负数上等价于「向下取整」；而 Python 的 // 也是向下取整、BigInt 的 / 却是向零截断——三者只在正数上一致。跨语言逐位复现要求两侧选同一种语义，本卷统一规定：降精度一律 >>。定点除法另起一个 fdiv，显式实现向零截断（先取绝对值相除、再补符号），Python 与 TypeScript 各写一份、行为对齐。第 2 章冒烟阶段实测过这个坑：用裸 / 会让两侧在负数上分叉，像素缓冲的 sha256 立刻不同。这里没有新东西，只是提醒——三角函数的表构造同样吃这条约定，sin_taylor_q16 里的每一次 fdiv 都依赖它。
欧拉角顺序：ZYX 与 XYZ 不是同一件事
把三个轴的旋转乘起来就得到欧拉角，但矩阵乘法不交换，顺序不同结果就不同。本卷采用 ZYX：R(α, β, γ) = Rz(α)·Ry(β)·Rx(γ)，作用在右侧顶点上——先转 x、再转 y、最后转 z（最右边的先作用）。换成 XYZ 顺序，同样的三个角度会给出另一个朝向。这个差异不是理论上的：Rz(30°)Ry(45°)Rx(60°) 与 Rx(60°)Ry(45°)Rz(30°) 的矩阵元素最大差达到 54647 LSB（见下方实跑输出），同一根探针向量在两种顺序下会落到相距约 2.1 个世界单位的两点。
实跑输出
$ /opt/homebrew/bin/python3 sim/cg_raster_04_1_rotation.py
表最大误差 : 2 LSB
sin²+cos² 最大偏差 : 6 LSB (在 15°) 容差 8
正交性 R·Rᵀ 最大偏差 : 4 LSB 容差 8
旋转保持长度 : 146542 -> 146542
欧拉顺序最大元素差 : 54647 LSB
欧拉顺序探针落点距离 : 137291 (Q16.16)
断言全部通过 (exit 0)
读数。sin²+cos² 的最大偏差 6 LSB、正交性偏差 4 LSB，都在声明容差 8 LSB 之内——注意这里说的是容差而不是「逐位精确」：旋转矩阵是感知量，允许它在几个 LSB 内浮动，只要不累积成可见形变。这里有一个刻意的区分：像素缓冲那种逐字节比对的目标必须 sha256 完全相同，而旋转矩阵这类派生量只承诺落在声明容差内，两者不能混为一谈。探针向量 (2, 1, 0) 绕 z 转 37° 后长度从 146542 变到 146542，一位不差，说明旋转确实保长。最后两行量化了顺序敏感性。
表误差、可视化与产物
S 层脚本除断言外还落三份产物到 lantern/out/：cg_raster_04_1_rotation.csv 逐角度列出 sin / cos 与 ident_err_lsb（即 sin²+cos²-ONE），cg_raster_04_1_rotation.json 汇总全部断言量，cg_raster_04_1_rotation.svg 用纯 xml.etree 画出正弦、余弦两条曲线（不引第三方绘图库）。把图与 CSV 一起看，可以逐点核对某个角度上的取值——偏差峰值出现在约 15° 处，正是某个插值段的中间，与「线性插值误差在段中央最大」的预期一致，不是随机噪声。可视化在这里不是装饰，而是让「表在哪个角度最不准」可被肉眼定位。
P 层：表的 sha256 与 Python 侧逐位一致
TypeScript 侧用 BigInt 独立重写同一套表构造（同样的「先乘后 >> 16」），断言 257 个条目的 sha256 与 Python 侧相同，并对 sin²+cos² 做容差断言。
// snippet: 4.1 P 层核心断言（完整文件 core/test/cg_raster_04_orientation.test.ts）
const TABLE_SHA256 = "efe8abce7db97b849eec5b5ccb624d846ee35ab9db9a52cb320980cb413912fd";
test("定点正弦表：257 个条目的 sha256 与 Python 侧逐位一致", () => {
assert.equal(createHash("sha256").update(SIN_TABLE.join(",")).digest("hex"), TABLE_SHA256);
assert.equal(SIN_TABLE[0], 0n);
assert.equal(SIN_TABLE[128], 46341n);
assert.equal(SIN_TABLE[256], 65536n);
});
$ node --test 'core/test/cg_raster_04_orientation.test.ts'
✔ 定点正弦表：257 个条目的 sha256 与 Python 侧逐位一致
✔ sin²+cos² = ONE：遍历全部 65536 个二进制角度，容差 8 LSB
✔ 负例：更粗的表（N_SEG=64）指纹不同，说明金标准确实在约束实现
ℹ tests 7 pass 7 fail 0
小结
本节把「旋转需要三角函数」这件事在禁浮点的前提下解决掉：嵌套泰勒在 Q16.16 上造出四分之一波表，查表加线性插值给出任意角的 sin / cos，二分查表给出反三角，全整数、可跨语言逐位复现。旋转矩阵由此直接搭出，两条不变量（sin²+cos²=ONE、R·Rᵀ≈I）在声明容差内成立，旋转保长一位不差。欧拉角顺序敏感，ZYX 与 XYZ 给出不同朝向——这也埋下伏笔：三个角度在某些姿态下会退化，即万向节死锁。
阅读导航：上一节：3.3 投影矩阵与 MVP 合成 · 下一节：4.2 四元数与万向节死锁 。
继续阅读
探索更多技术文章
浏览归档，发现更多关于系统设计、工具链和工程实践的内容。