本篇是《机器人炼金术》01 期的完整版包含直觉建立、完整推导、核心代码逐句解释、编写逻辑、四种实现方式的实测对比、以及位姿表示在工业界的全部典型应用场景。可独立阅读无需参考任何其他文章。一、问题定义机器人每时每刻都在回答我在哪、朝哪去机械臂末端在哪、朝向哪边无人机在空中什么姿态人形机器人躯干相对地面歪了多少——这些问题全是同一个形式一个刚体相对另一个刚体的位姿position orientation位置 朝向。位置好办三个数 (x, y, z) 是全部。麻烦的是朝向。朝向为什么难因为转不是三个独立数字能老实描述的东西。你可以用三个欧拉角但飞机俯仰 90° 时偏航和滚转突然变成同一个动作万向锁你可以用四元数四个数很优雅但两个四元数相加平均得到的不是中间姿态。本文要解决的具体问题是给机器人系统选一种位姿表示使得它 (a) 能无歧义地描述任意姿态 (b) 能被插值和求导 (c) 能进入优化器做数值更新 (d) 实现代价可控。这四个要求会淘汰掉大部分看起来能用的方案下面逐一展开。二、直觉层9 个数只有 3 个自由度最诚实的出发点是旋转矩阵一个 3×3 的正交矩阵9 个数。它满足两个约束RᵀR I列向量互相垂直且为单位长—— 6 个独立方程det(R) 1不许把空间翻成镜像—— 1 个方程这些约束不独立实际留给旋转的自由度是3。不必背这个减法直接做实验对旋转矩阵的 9 个元素做数值扰动看哪些方向动得了。# 实验代码配套脚本 se3_demo.py 的验证 [1] np.random.seed(0) Rm exp_so3(np.array([0.3, -1.2, 2.0])) # 造一个合法旋转矩阵 eps 1e-7 cols [] for i in range(3): # 对 3 个旋转参数分别扰动 d np.zeros(3); d[i] eps cols.append((Rm exp_so3(d) - Rm exp_so3(-d)) / (2 * eps)) J np.stack([c.reshape(-1) for c in cols], axis1) # 9x3 数值雅可比 print(np.linalg.matrix_rank(J, tol1e-6)) # 输出 3逐句解释第 2 行exp_so3把旋转向量变成矩阵3.2 节详述先造一个合法的旋转矩阵作为实验对象。第 4~6 行对 3 个旋转自由度分别加正负微小扰动用中心差分算矩阵 9 个元素随 3 个参数的变化率。第 7 行把三个方向的变化率拼成 9×3 的雅可比矩阵。第 8 行矩阵的秩是 3——9 维空间里只有 3 个方向能动其余 6 个方向一推就被正交性弹回来。这个冗余不是浪费是矩阵在替你守住这是一个旋转的承诺随手写一个 9 维向量多半根本不是旋转但随手写一个接近正交的矩阵它仍然是近似的旋转。这个性质后面做优化时有大用。左乘还是右乘新手第一大坑同样是再绕 z 轴转 10°有两种完全不同的含义写法名称物理含义R ← Rz(10°)·R左乘在世界坐标系里转机器人绕世界原点公转位置画圆R ← R·Rz(10°)右乘在本体坐标系里转机器人绕自己身上的轴自转位置不动实测验证配套脚本验证 [4]只加旋转扰动时左乘会顺带改变平移实测平移变化 [0.002462, 0.020649, 0]右乘不会[0, 0, 0]。记忆法矩阵从右往左读越靠右的动作发生在更本体的坐标系里。这不是数学趣味——它直接决定状态估计里为什么 IMU 融合几乎总用右扰动扰动只该影响姿态不该把已估计好的位置带偏第五节场景 2 展开。三、进阶层SE(3) 与李代数3.1 把位姿当成一个群齐次变换矩阵把旋转和平移打包成 4×4T 的上 3×3 块是旋转 R、右列是平移 t、底行 [0 0 0 1]记作 T ∈ SE(3)。SE(3) 是一个李群既是群可乘可逆又是光滑流形可求导、可优化。拆开讲群两个位姿复合还是位姿T₁T₂逆还是位姿T⁻¹。流形位姿的集合不是平地是弯曲曲面。在曲面上做线性插值加法更新都会滑出去——除非沿曲面的切空间走。李代数 se(3) 就是 SE(3) 在单位元处的切空间exp 映射负责从切空间走回曲面log 反之。工程上要记三条铁律铁律 1指数映射有闭式解Rodrigues 公式几行代码搞定这就是机器人库从不存欧拉角的原因——存旋转向量或矩阵更新时一步 exp 就回曲面。铁律 2log 在接近 π 时病态。实测正常角度 exp→log 往返误差 ~1e-15接近 π 时精度骤降且有符号翻转π 与 −π 是同一旋转。优化器姿态突然跳变 180°十有八九撞上了这里3.6 节完整复盘。铁律 3参数残差 ≠ 位姿残差。实测验证|φ| 较大时 log(exp(ξ)) 的参数值可能与 ξ 不同但回代后的位姿完全一致。优化时残差必须定义在位姿或李代数差上不是原始参数向量上——这是自研优化代码的高频翻车点。3.2 核心代码逐句解释exp_so3从旋转向量到旋转矩阵Rodrigues 公式def exp_so3(w): # (1) theta np.linalg.norm(w) # (2) if theta 1e-12: # (3) return np.eye(3) hat(w) # (3a) k w / theta # (4) K hat(k) # (5) return (np.eye(3) # (6) np.sin(theta) * K # (6a) (1 - np.cos(theta)) * (K K)) # (6b)(1)输入 w 是三维旋转向量方向旋转轴长度旋转角rad即李代数元素 [φ]× 的向量形式。(2)取模长得旋转角 θ‖φ‖旋转的大小与轴从这里分离。(3)小角度保护exp 级数 I[φ]×½[φ]ײ… 在 θ→0 时只剩一阶项 (3a)。不保护则 (4) 除零且噪声搅乱结果。阈值 1e-12 下近似误差 ~1e-24完全可忽略。(4)单位轴 k w/θ。(5)hat 把向量变反对称矩阵 [k]×让叉积变矩阵乘法[k]×v k×v。全部推导的枢纽。(6)Rodrigues 公式exp(θ[k]×) I sinθ·[k]× (1−cosθ)·[k]ײ。注意这不是级数截断是精确闭式解——反对称矩阵幂满足 K³−K级数每一项都归并进这三项。这也是机器人学库普遍用它而不用通用矩阵指数的原因快一个数量级且没有特征值分解的数值副作用。exp_se3从 twist 到 4×4 位姿def exp_se3(xi): # (1) rho, phi xi[:3], xi[3:] # (2) theta np.linalg.norm(phi) # (3) J (np.eye(3) # (4) (1 - np.cos(theta)) / theta**2 * K # (4a) (theta - np.sin(theta)) / theta**3 * (K K)) # (4b) T[:3, :3] exp_so3(phi) # (5) T[:3, 3] J rho # (6) return T(1)输入 xi 是 twist六维速度前 3 维平移 ρ后 3 维旋转 φ。顺序约定坑GTSAM/本代码平移在前Sophus 旋转在前对接外部库先核对此事错序是高频 bug。(3)旋转部分模长 θ 仍是总旋转角。(4)矩阵 J 是 SE(3) 左雅可比 J_l(θ)见 (4a)(4b)回答绕固定轴转 θ 的过程中平移积累多少。系数来自对 ∫₀¹ exp(s[φ]×)ds 逐项积分——平移不是把 ρ 抄进第四列因为旋转和平移同时发生完整推导见 3.5。(5)旋转块复用 exp_so3。(6)平移块 J_l(θ)ρ。直觉ρ 是随旋转发生的直线运动J 把它扭成螺旋线。SE(3) 元素的几何本质是螺旋运动Chasles 定理任何刚体位移 一次螺旋。3.3 编写逻辑这两个函数面对的约束写这两个函数时的决策过程比函数本身更值得学约束 a零重依赖。目标读者 clone 即跑所以只允许 numpy——这直接排除了直接调 Sophus/GTSAM的省事路线。约束 b教学透明。每行必须能对应到公式出问题逐行可查——这排除了用 scipy 一行解决的路线对比见 3.4。小角度分支不是可选项。第一版没有它θ1e-13 的测试直接除零崩溃。数值代码的边界分支永远不是防御性编程的装饰是正确性的一部分。输出约定SE(3) 用 4×4 矩阵而非 (R,t) 元组因为下游的复合运算 T₁T₂ 和求逆 T⁻¹ 对矩阵是原生的对元组要拆装——数据结构选型跟着最高频操作走。3.4 四种实现方式实测对比同一功能旋转向量 → 旋转矩阵四种写法的实测数据200 个随机旋转向量角度 0~3 radscipy 1.17 / numpy 2.3 / Python 3.11实现最大误差单次耗时*依赖评价① Rodrigues 公式本文5e-14~2.2 µsnumpy精确闭式解教学与轻量场景首选② 泰勒级数截断8 项2.3e3~4.3 µsnumpy大角度直接发散禁用③ scipy.linalg.expm(hat(w))0基准~2.7 µsscipy最通用但对反对称阵是杀鸡用牛刀④ scipy.Rotation.from_rotvec5e-14~2.5 µsscipy生产环境首选库维护、批量接口、C 实现Python 层耗时被解释器开销主导四种差距不大C 环境里差距拉大一个数量级*——Rodrigues 是十几条浮点指令expm 是特征分解级别的重活。1 kHz 控制环每周期几十次变换C 里这是真实成本。每种实现的区别与边界② 为什么崩矩阵指数级数 Σ Kⁿ/n! 只在 ‖K‖ 小时收敛快。反对称阵的范数就是旋转角θπ 时 8 项远不够误差按 θⁿ/n! 爆炸——实测 2.3e3 的误差就是这个来源。所有正规库都用闭式解或帕德逼近没人用裸级数。③ vs ①expm 内部走帕德逼近缩放平方为一般矩阵设计对反对称阵这些 machinery 全是浪费。但遇到非反对称的一般矩阵指数如时变线性化系统的 exp(At)③ 是正确工具。工具没有高下只有场景。④ vs ①性能无明显差异Python 层但④久经测试、双覆盖/万向锁边界都处理过。生产用库教学手写一遍不可省——只会调库的人读不懂 VIO 源码也写不出流形优化。最终选用及理由本文手写①服务于零依赖教学透明的目标从下一期起涉及批量变换时切换到④并重跑对比表。3.5 完整推导exp_se3 里的 J 矩阵从哪来免费渠道通常只给公式这里给出完整推导链让你能自己重推。出发点twist ξ [ρ; φ] 对应矩阵 [ξ]∧上块 [φ]× 和 ρ底行 0。设机器人从原点出发以恒定 body 速度 ξ 运动位姿满足微分方程Ṫ T·[ξ]∧初值 T(0)I。拆成旋转和平移两个方程Ṙ R·[φ]×标准 SO(3) 指数方程解 R exp(θ[k]×)Rodrigues。ṫ R(t)·ρ把 R 代入后对时间积分。用变量代换 s∈[0,1]body 速度恒定时旋转按比例累积t ∫₀¹ exp(s[φ]×)·ρ ds把被积函数按级数展开并逐项积分t Σₙ [φ]×ⁿ/(n1)! · ρ利用反对称矩阵幂循环 K³−K级数按 n mod 3 归并成三类归并需要两个三角恒等式可用 e^(iθ) 实部/虚部验证sinθ/θ Σ (−1)ⁿθ²ⁿ/(2n1)!(1−cosθ)/θ² Σ (−1)ⁿθ²ⁿ/(2n2)!最终归并出J_l(θ) I (1−cosθ)/θ²·[φ]× (θ−sinθ)/θ³·[φ]ײ——正是代码 (4a)(4b) 两行。三个常见错法重点把 [φ]× 当单位轴矩阵用公式里的 K 含模长 θ。若先归一化了轴系数必须换成等价变体。两种写法都对混用必错。θ²、θ³ 直接除θ→0 时爆炸必须走小角度分支J_l → I ½[φ]×。J_l 和 J_r 混用J_r(θ) J_l(−θ)。积分定义速度在 body 系出 J_l世界系定义速度要用 J_r。错用的症状小角度都对大角度位置越走越歪。3.6 log 的边界病态一次真实的实现翻车理论说完看一次真实的数值事故——它揭示的方法论比结论更值钱。第一版实现教科书写法θ arccos((tr(R)−1)/2)轴由 vee(R−Rᵀ) 2sinθ·k 恢复接近 π 时特判走 RI 2kkᵀ 秩 1 分支。写完跑边界压力测试θ∈[π−1e-6, π) 采样 5 万例最差位姿回代误差: 0.0036 rad 预期 ~1e-9误差超预期7 个数量级。抓最差样本看中间量输入 θ π - 1e-6落在普通分支未触发特判 cosθ 真值: -0.9999999999995001 矩阵算出: -0.9999999999994993 ← 差 8e-16正常舍入 log 返回: 模长 3.1390应为 3.14159差 2.6e-3定位误差放大链条arccos 在 x→−1 处导数发散|d arccos/dx| 1/√(1−x²) 1/sinθ。本例 sinθ≈10⁻⁶角度误差 ≈ 8×10⁻¹⁶/10⁻⁶ 8×10⁻¹⁰ rad。轴恢复系数 θ/(2sinθ) 的分子分母都用了病态的 θδθ 传导给 sinθ 一个 cosθ·δθ 的误差相对 sinθ10⁻⁶ 就是8×10⁻⁴ 的相对误差——正是实测量级。教训误差不在 arccos 一步产生而是在病态量的下游每一步被复利。修分支阈值治标不治本阈值放宽到 1e-3θπ−1e-4 时同样的戏码重演只是误差小一点。根治换全程良态的参数化路径Shepperd 四元数方法tr np.trace(Rm) if tr 0: # (S1) s 2.0 * np.sqrt(tr 1.0) # s 4*qw qw 0.25 * s qx (Rm[2,1] - Rm[1,2]) / s qy (Rm[0,2] - Rm[2,0]) / s qz (Rm[1,0] - Rm[0,1]) / s else: # (S2) i int(np.argmax(np.diag(Rm))) # 选最大对角元分支 j, k (i1) % 3, (i2) % 3 s 2.0 * np.sqrt(1.0 Rm[i,i] - Rm[j,j] - Rm[k,k]) # s 4*q_i qw (Rm[k,j] - Rm[j,k]) / s q[i] 0.25 * s q[j] (Rm[j,i] Rm[i,j]) / s q[k] (Rm[k,i] Rm[i,k]) / s qv np.array([qx, qy, qz]) theta 2.0 * np.arctan2(np.linalg.norm(qv), abs(qw)) # (S3) axis qv / np.linalg.norm(qv) * (1.0 if qw 0 else -1.0) # (S4)(S1)tr(R) 4qw²−1trace0 意味着 |qw|1/2走 qw 主导分支分母离零最远。各 qᵢ 由反对称差 R₂₁−R₁₂ 4qw·qx 类公式恢复——分子小但分母更大商良态。(S2)trace≤0 时 |qw|≤1/2改从三个 qᵢ 里选最大的做分母4qᵢ² 1Rᵢᵢ−Rⱼⱼ−Rₖₖ选最大 Rᵢᵢ 保证根号内不出现危险的小数相减。Shepperd 方法的核心永远用最大的四元数分量当除数。(S3)half-angle 公式 θ 2·atan2(‖qv‖, |qw|)。arctan2 全定义域条件数良好且 qv sin(θ/2)·k 直接编码半角正弦×轴不经过任何病态中间量。(S4)四元数双覆盖q 与 −q 是同一旋转。若 qw0这个 q 表示绕 −k 转 θ轴要翻回来否则 θ∈[0,π] 配错误轴向会给出镜像旋转。修复后实测[boundary] 20000 例 θ∈[π-0.05,π)位姿回代 100% 通过 [normal-zone] 20000 例正常区往返最大误差 1.26e-15旧实现 1.38e-11 [sign-flip] θπ-1e-7 时 log 返回 [3.141593, 0, 0] —— ±π 双值位姿等价正常区精度也提升 4 个数量级——良态化不是只修边界是全区间受益。剩下的双值性是数学给的±π 处 k 与 −k 表示同一旋转SO(3) 基本群是 ℤ₂任何实现都消不掉。工程后果三条优化器符号跳变loss 周期性尖刺尖刺间隔旋转累计到 π 的时刻、多旋转平均要先统一半球否则接近 π 的旋转正负抵消均值趋零、slerp 不处理双覆盖会走远弧相机从朝前到朝后动画绕远路转 350° 而不是 10°。四、三种位姿表示的完整对比表示维数插值万向锁适用场景欧拉角3直接线性插值有俯仰 ±90° 退化人机交互、日志可读四元数4slerp/nlerp无运行时状态、IMU旋转矩阵/李代数9/3流形插值无优化、复合变换实测验证 [5]绕倾斜轴从单位阵转到 90°比较取一半时各方法偏离真实半程姿态的角度——矩阵线性混合(SVD投影后) 0.000°欧拉角插值 10.660°四元数 nlerp 0.000°。欧拉角误差按 1−cos(θ/2) 增长多轴耦合时任何涉及平滑转向的功能云台、动画、轨迹插值必须走四元数或流形插值。五、应用场景位姿表示具体解决什么问题以下五个场景按概念→具体问题→现象→根因→解法→验证完整展开。场景 1ROS/tf2 坐标变换——到位但姿态歪 10°具体问题机械臂要抓零件规划器输出的末端位姿转换到基座标系后总是歪。现象位置准确误差 1mm姿态稳定偏 10° 左右。根因链路上两个环节的欧拉角顺序约定不同一个 ZYX 一个 XYZ或左乘右乘用反。顺序错 → 绕中间轴的旋转顺序颠倒 → 特定姿态下误差恰好是轴交换的量级。解法全链路统一为四元数或矩阵传递打印中间位姿时附上 det(R) 和 ‖RᵀR−I‖ 做健康检查不是 1 和 0 说明上游把矩阵当普通数组操作了误差会随复合指数放大。验证用一个已知姿态如绕 z 转 90°走完整条变换链输出应逐位复现。场景 2IMU 姿态融合——为什么用右扰动具体问题IMU 以 200Hz 输出角速度增量李代数侧融合器维护全局姿态李群侧每次更新用 δφ 修正。具体解决的是什么更新量应该只影响姿态。若用左扰动 R←exp(δφ)R绕世界轴转位置平移分量也会被带偏——实测左乘引起平移变化 0.02m 而右乘为 0。现象对比左扰动实现的融合器静止时位置估计漂移比右扰动版本大数倍且漂移方向与旋转相关。解法误差状态卡尔曼滤波ESKF的标准做法R ← R·exp(δφ)δφ 作为误差状态。这是 VIO/SLAM 库的事实标准。场景 3优化器里的流形更新——越优化越远具体问题手写位姿图优化把位姿存成 9 维向量用普通高斯牛顿更新。现象前几次迭代 loss 下降之后开始震荡甚至发散检查发现矩阵不再正交det 偏离 1。根因加法更新把估计值推出了流形——在曲面上做线性外推必然滑出去。解法位姿存矩阵更新量存李代数R ← R·exp(δξ)δξ∈R⁶ 是优化变量。Ceres/g2o/GTSAM 的流形优化全是这个结构。验证每次迭代后断言 ‖RᵀR−I‖ 1e-10。场景 4轨迹插值——云台运动绕远路具体问题相机云台从朝前平滑转到朝后。现象用欧拉角线性插值云台绕远路转 350°或路径抖动。根因欧拉角不是流形上的坐标其线性组合没有几何意义多轴耦合时中间姿态偏离大圆路径实测单轴就有 10.66° 偏差。解法四元数 slerp处理双覆盖q₀·q₁0 时翻转 q₁或李代数插值 exp(log(R₁ᵀR₂)/2)·R₁。验证中间帧用 arccos(tr(Rᵀ·R_target)/2) 检查角度应单调。场景 5手眼标定——AXXB 的结构从哪来具体问题求相机相对机械臂法兰的安装位姿 X。现象标定残差大、每次标定结果差异大。根因AXXB 里 A、B、X 全是 SE(3) 元素朴素地把 12 个数当未知量线性求解忽略了 X 的流形结构噪声敏感。解法在李代数上求解对数线性化或先用四元数解旋转部分再解平移多组数据平均时必须先统一半球见 3.6 双值性。验证标定后用新姿态做重投影残差应 2px。这是番外案例 A 的主题届时完整实操。六、可复现代码清单文件内容运行方式code/se3_demo.py5 组数值验证自由度、exp/log 往返、左右扰动、插值对比 出图python3 se3_demo.py --plotcode/se3_boundary_test.py2 万例边界压力测试 双值性测试 性能基准python3 se3_boundary_test.py全部仅依赖 numpy/scipy实测输出已写入本文各节。七、小结与下一篇位姿表示的选型逻辑交互看欧拉角、运行看四元数、优化看流形三者用 exp/log 打通数值代码的方法论病态 → 实测误差量级 → 抓最差样本 → 定位放大链条 → 换参数化路径 → 全区间回归下一篇02 期《李群与李代数为什么工程师必须懂 so(3) 的指数映射》——扰动模型的雅可比、J_l⁻¹ 闭式解、在 Ceres 残差块里的具体写法附录代码 文 01 付费版配套log_so3 在 θ→π 边界的压力测试 依赖numpy 运行python3 se3_boundary_test.py import numpy as np np.set_printoptions(precision6, suppressTrue) def hat(w): return np.array([[0.0, -w[2], w[1]], [w[2], 0.0, -w[0]], [-w[1], w[0], 0.0]]) def vee(W): return np.array([W[2, 1], W[0, 2], W[1, 0]]) def exp_so3(w): theta np.linalg.norm(w) if theta 1e-12: return np.eye(3) hat(w) k w / theta K hat(k) return np.eye(3) np.sin(theta) * K (1 - np.cos(theta)) * (K K) def log_so3(Rm): Shepperd 四元数路径与主文件一致全程良态。 tr np.trace(Rm) if tr 0: s 2.0 * np.sqrt(tr 1.0) qw 0.25 * s qx (Rm[2, 1] - Rm[1, 2]) / s qy (Rm[0, 2] - Rm[2, 0]) / s qz (Rm[1, 0] - Rm[0, 1]) / s else: i int(np.argmax(np.diag(Rm))) j, k (i 1) % 3, (i 2) % 3 s 2.0 * np.sqrt(1.0 Rm[i, i] - Rm[j, j] - Rm[k, k]) qw (Rm[k, j] - Rm[j, k]) / s q [0.0, 0.0, 0.0] q[i] 0.25 * s q[j] (Rm[j, i] Rm[i, j]) / s q[k] (Rm[k, i] Rm[i, k]) / s qx, qy, qz q qv np.array([qx, qy, qz]) n np.linalg.norm(qv) if n 1e-12: return 2.0 * qv theta 2.0 * np.arctan2(n, abs(qw)) axis qv / n * (1.0 if qw 0 else -1.0) return theta * axis def test_boundary(n20000): θ ∈ [π-0.05, π) 随机采样判据是位姿回代一致容忍双值 rng np.random.default_rng(42) fails [] param_flips 0 for _ in range(n): w rng.normal(size3) w w / np.linalg.norm(w) * rng.uniform(np.pi - 0.05, np.pi) out log_so3(exp_so3(w)) if not np.allclose(exp_so3(out), exp_so3(w), atol1e-9): fails.append((w, out)) # 参数级双值out 与 w 方向相反也算合法位姿相同 if out w 0 and np.allclose(exp_so3(out), exp_so3(w), atol1e-9): param_flips 1 assert not fails, f{len(fails)} 例位姿回代失败 print(f[boundary] {n} 例边界采样位姿回代 100% 通过 f其中 {param_flips} 例参数符号翻转±π 双值合法) def test_sign_flip(): 双值性直接验证 w np.array([np.pi - 1e-7, 0.0, 0.0]) out log_so3(exp_so3(w)) assert abs(abs(out[0]) - np.pi) 1e-4, out assert np.allclose(exp_so3(out), exp_so3(w), atol1e-9) print(f[sign-flip] θπ-1e-7 时 log 返回 {out} —— ±π 双值位姿等价) def test_perturbation_free_zone(): 对照远离边界的正常区参数回代应严格一致 rng np.random.default_rng(7) worst 0.0 for _ in range(20000): w rng.normal(size3) w w / np.linalg.norm(w) * rng.uniform(0, np.pi - 0.1) worst max(worst, np.linalg.norm(log_so3(exp_so3(w)) - w)) print(f[normal-zone] 正常区 20000 例参数往返最大误差 {worst:.2e}应 ~1e-11 量级) def bench(): 边界分支的正常分支性能对照 import timeit ws [] rng np.random.default_rng(1) for _ in range(2000): w rng.normal(size3) ws.append(w / np.linalg.norm(w) * rng.uniform(0, np.pi)) t timeit.timeit(lambda: [log_so3(exp_so3(w)) for w in ws], number50) / (50 * len(ws)) print(f[bench] explog 单次往返平均 {t*1e6:.2f} µs) if __name__ __main__: test_boundary() test_sign_flip() test_perturbation_free_zone() bench()