4.2 四元数与万向节死锁
上一节的 ZYX 欧拉角把朝向交给三个角度，直观，却在某些姿态下会「塌掉」——这就是万向节死锁。本节先把死锁从一句口号变成三个可断言的数字（Jacobian 行列式、旋转轴夹角、世界朝向差），再给出四元数的表示与乘法，说明它为什么不退化。四元数的工程用途远不止旋转：游戏里整条骨骼层级靠它避免逐骨骼累积的欧拉误差，那部分实践见 骨骼动画里的四元数 ；本节只关心一件事——如何用整数把「死锁」演示出来，并让一个可证伪的断言盯着它。
四元数的表示与乘法
单位四元数 q = (w, x, y, z) 表示绕单位轴 n 转 θ：
q = (cos(θ/2), sin(θ/2) · n)
注意是半角。在二进制角度下 θ/2 就是 (a_bam >> 1) & (ONE - 1)——右移一位，越过 180° 时自动折回，不需要任何特殊处理；cos / sin 直接查上一节的表。旋转一个向量用 v' = q·v·q*（v 取实部为 0 的纯四元数），两次 Hamilton 乘即可；朝向也能转成 4×4 矩阵，与前面几章的管线无缝对接。
Hamilton 乘法的分量式是四元数的核心：
(w₁w₂ - x₁x₂ - y₁y₂ - z₁z₂,
w₁x₂ + x₁w₂ + y₁z₂ - z₁y₂,
w₁y₂ - x₁z₂ + y₁w₂ + z₁x₂,
w₁z₂ + x₁y₂ - y₁x₂ + z₁w₂)
它把「先转 A 再转 B」合成一个四元数，且旋转方向与矩阵右乘一致。
四个分量各自是什么，值得单独说清。实部 w = cos(θ/2) 只记录转了多大角（半角），三个虚部 (x, y, z) = sin(θ/2)·n 记录绕哪根轴——整根向量 (x, y, z) 的方向是转轴，长度是 sin(θ/2)。于是单位约束 w² + x² + y² + z² = 1 不是附加条件，而是「半角的正弦平方加余弦平方恒为 1」的必然结果。反过来读：任何单位四元数都能还原出 θ = 2·acos(w) 与轴 (x, y, z)/sin(θ/2)，表示与朝向之间是可逆的。
整段实现如下（工作区文件 lantern/sim/cg_raster_04_quat.py）：
#!/usr/bin/env python3
# file: lantern/sim/cg_raster_04_quat.py
# tier: S
# verify: /opt/homebrew/bin/python3 sim/cg_raster_04_2_quaternion.py
"""四元数（Q16.16，全整数）：表示、乘法、旋转向量、转矩阵、nlerp、夹角。
约定：q = (w, x, y, z)，单位四元数；绕轴 n 转 θ 的表示为
q = (cos(θ/2), sin(θ/2)·n)。旋转向量 v 用 v' = q·v·q*（v 取 w=0 的纯四元数）。
与 TypeScript 侧 core/test/cg_raster_04_orientation.test.ts 的 quat 实现逐字对齐。
"""
import math
import os
import sys
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import cg_raster_04_trig as trig # noqa: E402
from cglib import ONE, fdiv, fmul, fsqrt # noqa: E402
def q_from_axis_angle(axis, a_bam):
"""绕 axis（单位向量）转 a_bam（二进制角度）的四元数。"""
half = (a_bam >> 1) & (ONE - 1)
c, s = trig.cos_bam(half), trig.sin_bam(half)
return (c, fmul(s, axis[0]), fmul(s, axis[1]), fmul(s, axis[2]))
def q_mul(a, b):
"""Hamilton 积。"""
aw, ax, ay, az = a
bw, bx, by, bz = b
return (
fmul(aw, bw) - fmul(ax, bx) - fmul(ay, by) - fmul(az, bz),
fmul(aw, bx) + fmul(ax, bw) + fmul(ay, bz) - fmul(az, by),
fmul(aw, by) - fmul(ax, bz) + fmul(ay, bw) + fmul(az, bx),
fmul(aw, bz) + fmul(ax, by) - fmul(ay, bx) + fmul(az, bw),
)
def q_conj(q):
return (q[0], -q[1], -q[2], -q[3])
def q_dot(a, b):
return fmul(a[0], b[0]) + fmul(a[1], b[1]) + fmul(a[2], b[2]) + fmul(a[3], b[3])
def q_norm(q):
n = fsqrt(q_dot(q, q))
if n == 0:
return (ONE, 0, 0, 0)
return tuple(fdiv(c, n) for c in q)
def q_rotate(q, v):
"""用单位四元数 q 旋转向量 v。"""
r = q_mul(q_mul(q, (0, v[0], v[1], v[2])), q_conj(q))
return (r[1], r[2], r[3])
def q_to_mat(q):
"""单位四元数 -> 4x4 旋转矩阵（行主序 Q16.16）。"""
w, x, y, z = q
x2, y2, z2 = fmul(x, x), fmul(y, y), fmul(z, z)
return (
ONE - 2 * (y2 + z2), 2 * (fmul(x, y) - fmul(w, z)), 2 * (fmul(x, z) + fmul(w, y)), 0,
2 * (fmul(x, y) + fmul(w, z)), ONE - 2 * (x2 + z2), 2 * (fmul(y, z) - fmul(w, x)), 0,
2 * (fmul(x, z) - fmul(w, y)), 2 * (fmul(y, z) + fmul(w, x)), ONE - 2 * (x2 + y2), 0,
0, 0, 0, ONE,
)
def q_neg(q):
return (-q[0], -q[1], -q[2], -q[3])
def q_shortest(q0, q1):
"""双覆盖：若 dot(q0,q1) < 0，翻转 q1 走最短路径。"""
return q_neg(q1) if q_dot(q0, q1) < 0 else q1
def q_nlerp(q0, q1, t):
"""整数 nlerp：先在单位球面上取最短路径，再线性插值并归一化。t ∈ [0, ONE]。"""
q1 = q_shortest(q0, q1)
lerp = tuple(q0[i] + fmul(q1[i] - q0[i], t) for i in range(4))
return q_norm(lerp)
def q_angle_bam(q0, q1):
"""两个单位四元数所代表朝向之间的夹角（二进制角度，∈ [0, 2·QUARTER]）。
用相对四元数 r = q0* · q1 的向量部：|r_vec| = sin(θ/2)，故 θ = 2·asin(|r_vec|)。
小角度下 |r_vec| 与 θ 近似线性，分辨率远优于用 dot 反算 acos（后者在 θ→0 时陡峭）。
"""
r = q_mul(q_conj(q0), q1)
# 向量部各分量的平方是 Q32.32；其和开方即 Q16.16 的 |vec|（不能用 fmul，小分量会被 >>16 抹掉）
v = math.isqrt(r[1] * r[1] + r[2] * r[2] + r[3] * r[3])
if v >= ONE:
return 2 * trig.QUARTER
return 2 * trig.asin_bam(v)
万向节：这个名字从哪来
「万向节」（gimbal）原本是机械词：陀螺仪里三层同心圆环，外环绕一个轴转、中环绕第二个轴、内环绕第三个轴，串起来让最里层的物体可以指向任意方向。问题出在中间那层环转到 ±90° 时——外环与内环的转轴变成共线，两个环的转动合成同一个效果，整套机构只剩两个自由度。航海陀螺仪因此会在特定姿态「失去一个轴」，这就是 gimbal lock 的原始场景。欧拉角的 ZYX 参数 (α, β, γ) 与这三层环一一对应：β 就是中环。所以死锁不是欧拉角实现的缺陷（换一种数值精度也救不回来），而是「用三个角度串三个轴」这套表示法的固有拓扑性质——它属于表示，不属于代码。
把「死锁」量化为三个数字
其一，自由度重合。 世界角速度对参数速度 (α̇, β̇, γ̇) 的 Jacobian 行列式恰好等于 -cos β。β = 0° 时 |det| = 1，三个参数各自对应一条线性无关的世界旋转轴；β = 90° 时 det = 0——三个轴降到秩 2，α 与 γ 的旋转轴从互相垂直（夹角 90°）变成反向共线（夹角 180°），两个自由度重合。这就是「死锁」的严格含义：不是数值不稳，而是参数到朝向的映射在这一点降秩。
其二，死区。 β = 90° 时 R 只通过 (γ - α) 依赖 α 与 γ。于是让 α 扫过 90°：若 γ 同时 +90°，世界朝向完全不动（差 0°）；若 γ 取 -90°，世界朝向却转了 180°。同样的参数运动量（各 180°），世界结果相差 180°——参数空间里走得一样远，朝向空间里一个没动、一个转了大半圈。
其三，四元数无此病。 朝向到四元数的映射唯一（只差 q 与 -q 的符号），不存在一整条流形塌到一点的退化。从 q0 到 q1 做 nlerp，逐帧世界转角单调不减、端点精确命中，没有任何一个参数方向的运动会「白走」。
死锁的精确值与常见触发
-cos β 这个行列式在 Q16.16 下有精确取值：脚本产出的 jacobian_det_beta0_q16 = -65536 就是 -ONE（即 -1.000），jacobian_det_beta90_q16 = 0 就是 0——不是「接近零」，是位模式全零。β = ±90° 是死锁点，|β| 越接近 90° 越接近退化，所以粗略地说，在 80°~100° 这个区间做欧拉插值都要格外小心。三种常见触发场景：相机俯仰拉到正上方或正下方直视（β = ±90°）、角色抬头看天或低头看脚、以及两个欧拉角速度分量同时打得很大——最后一种在死锁点附近会让朝向突然「翻面」，因为主导项 (γ - α) 的符号变了，世界结果从「不动」跳到「转 180°」。
四元数接进管线
四元数最终要变回矩阵才能喂给顶点管线：q_to_mat 把单位四元数展开成 4×4，形式与上一节的三轴旋转矩阵同类（行主序、最后一行与最后一列是单位）。它比「先算欧拉角再连乘三个矩阵」少两次 mat_mul，也少两次 >> 16 的舍入——旋转一个顶点时，矩阵路径是「三矩阵相乘 + 一次 mat_vec4」，四元数路径是「两次 Hamilton 乘 + 一次 mat_vec4」，中间量更少、误差源更少。本卷后续凡是需要「朝向」的地方（相机、骨骼）都以四元数为内部表示，只在需要出示给用户或落盘时才转成欧拉角。
实跑输出
$ /opt/homebrew/bin/python3 sim/cg_raster_04_2_quaternion.py
β=0° α 与 γ 旋转轴夹角 : 90° Jacobian det = -1.000
β=90° α 与 γ 旋转轴夹角 : 180° Jacobian det = +0.000 <- 降秩
死区 β=90°，α 各扫 90°：γ=+90° -> 世界差 0°；γ=-90° -> 世界差 180°
（参数运动量都是 180°，世界结果相差 180°）
四元数 nlerp：单调=True 端点误差=0° 逐帧步长 672~1302 BAM
断言全部通过 (exit 0)
读数。det 从 -1.000 掉到 +0.000 是本节最硬的证据：它不是「看起来会抖」，而是一个行列式从满秩归零。夹角 90° → 180° 说明两条旋转轴从垂直变成反向共线。死区那一行把「参数动、世界不动」与「同样参数动、世界却转 180°」并排摆出来，一眼可见哪个更像 bug。四元数一侧，nlerp 逐帧步长 672~1302 BAM（约 3.7°~7.2°），单调、端点精确。
P 层：四元数逐位一致与负例
TypeScript 侧用 BigInt 独立重写四元数内核，断言轴角构造、乘法、旋转向量、转矩阵与 nlerp 中点都与 Python 侧逐位相同。
// snippet: 4.2 P 层四元数断言（完整文件 core/test/cg_raster_04_orientation.test.ts）
test("四元数：轴角构造、乘法、旋转向量、转矩阵与 Python 侧逐位一致", () => {
const q1 = qFromAxisAngle(AX_Y, deg2bam(90));
assert.deepEqual(q1, [46341n, 0n, 46341n, 0n]);
const rx = qRotate(q1, AX_X); // 绕 Y 转 90° 把 +X 送到 -Z（容差 2 LSB）
assert.ok(rx[0] === 0n && rx[1] === 0n && rx[2] <= -65536n && rx[2] >= -65538n);
});
$ node --test 'core/test/cg_raster_04_orientation.test.ts'
✔ 四元数：轴角构造、乘法、旋转向量、转矩阵与 Python 侧逐位一致
✔ nlerp：端点精确命中，中点与 Python 侧逐位一致
✔ 双覆盖：q 与 -q 同朝向，最短路径必须翻转符号
ℹ tests 7 pass 7 fail 0
这里值得记一笔：qRotate 把 +X 送到 -Z 时结果是 -65538 而不是 -65536，差 2 LSB——四元数旋转经两次 Hamilton 乘、共 32 次 fmul，每次 >> 16 都向下取整，误差就这么攒出来。所以断言写成区间（-65538 ≤ z ≤ -65536）而不是等值。这正是「派生量用容差、逐字节量用 sha256」这条分界在四元数上的体现。
为什么还在用欧拉角
既然四元数更好，为什么游戏引擎的编辑器里还到处是欧拉角？因为人。美术给一个「稍微向右看一点」的调整，用 yaw += 10° 一句话就说清；用四元数得描述一根轴和一个角，直觉成本高得多。于是工程上的分工是：编辑期用欧拉角（可读、可增量调），运行期转成四元数（可插值、不退化），需要时再转回欧拉角显示。这条「authoring 用欧拉角、runtime 用四元数」的分界，与本卷「S 层脚本用 Python、顶点管线用 BigInt」的分工是同一种思路——表示形式跟着使用场景走，不追求一种表示打天下。另外，欧拉角还是可数的：三个数各有明确物理含义（俯仰、偏航、滚转），调试时打印出来一眼能读懂；四个分量的四元数打印出来只有一堆难以直觉的数。
半角为什么是必须的
一个自然的疑问：为什么用 θ/2 而不是 θ？因为要让四元数乘法同态于旋转的合成。两个旋转先 A 后 B，合成角是 θ_A + θ_B，而半角恰好让余弦的和角公式成立：cos(θ_A/2)·cos(θ_B/2) - sin(θ_A/2)·sin(θ_B/2) = cos((θ_A + θ_B)/2)。若改用整角，两次旋转的合成会得到 2θ 的余弦，乘法的语义立刻崩掉。半角还顺带解释了双覆盖：θ 与 θ + 360° 表示同一旋转，但半角相差 180°，四元数恰好差一个符号——q 与 -q 同朝向，根子就在这里。
欧拉角与四元数：一张对照表
把两种表示放在一起对比，取舍一目了然：
| 维度 | ZYX 欧拉角 | 单位四元数 |
|---|---|---|
| 参数个数 | 3 | 4（含 1 个单位约束） |
| 表示是否唯一 | 否（无穷多组等价） | 是（只差 ± 符号） |
| 是否有退化点 | 有（β = ±90°） | 无 |
| 插值是否均匀 | 否（受死锁影响） | nlerp 略不均匀，slerp 恒定 |
| 可读性 | 高（俯仰 / 偏航 / 滚转） | 低（四个难直觉的数） |
| 旋转一个顶点 | 三矩阵相乘 | 两次 Hamilton 乘 |
本节的演示集中在前四行：det 归零是「有退化点」的直接证据，死区是「表示不唯一」在插值上的后果，四元数 nlerp 的单调性对应「无退化点」。后两行则是工程取舍——可读性交给编辑期，运行时用四元数。
这条链为什么可证伪
本节的每一个结论都对应一个可复算的整数。det = -ONE 与 det = 0 是位模式级别的断言，不是「看起来接近」；死区的 0° 与 180° 是把两个欧拉角三元组分别转成四元数再量夹角得到的；四元数的单调性与端点精确性由逐帧累积转角验证。判据有区分力：把 β 从 90° 挪开一点点，det 立刻离开 0；把最短路径的符号翻转去掉，nlerp 的终点仍对但过程角度会翻倍。能拒绝一个「看起来合理」的改动的断言，pass 才可信——这也是本卷 S 层脚本一律带负例的原因。
小结
万向节死锁不是玄学，是 Jacobian 行列式 -cos β 在 β = ±90° 处归零、两个自由度重合的结果；它的可观测后果是死区——参数走一样远，世界朝向却不动或转 180°。四元数用四个数加一个单位约束表示朝向，映射唯一（只差符号），nlerp 路径单调、端点精确，从根上避开这条退化。
一句话记住两种表示的分工：欧拉角负责让人看懂，四元数负责让机器算对。
下一节把「从 q0 插值到 q1」做完整：nlerp 与 slerp 的速度差异、q 与 -q 的双覆盖，以及为什么必须先取最短路径。
阅读导航：上一节：4.1 旋转矩阵与欧拉角 · 下一节：4.3 朝向的插值 。
继续阅读
探索更多技术文章
浏览归档，发现更多关于系统设计、工具链和工程实践的内容。