四元数计算中小浮点数的Python除法

Question

考虑数学方程：y = (sin(x/2)/x)，它出现在旋转计算中。在 x->0 的限制下，其值为 1/2。在 Python（或任何具有类似 64 位浮点类型的语言）中实现这一点似乎对于除 x = 0.0 之外的所有值都可以正常工作。

然而，当这个数学出现在书中时，通常建议使用不存在被零除问题的截断泰勒级数展开式来实现该函数。这是 SciPy 中的一个示例，其中为此逻辑选择了看似任意的 1e-3 弧度阈值。

基本实现

y = np.sin(x/2)/x

似乎对于较小的 x 值工作得很好（显然除了 x = 0.0 之外）。这里真的需要担心精度问题吗？如果是这样，显示它的测试用例是什么？

检查 0.0 作为特殊情况的代码比使用具有任意阈值的截断泰勒级数更清晰。

Answer 1

基本实现
y = np.sin(x/2)/x
似乎对于较小的 x 值工作得很好（显然除了 x = 0.0 之外）。这里真的需要担心精度问题吗？如果是这样，显示它的测试用例是什么？

是的。您可以在以下代码示例中看到这一点：

>>> angle = 5e-324
>>> np.sin(angle/2)/angle
0.0

更一般地说，这是从 -5e-322 到 5e-322 绘制的函数

np.sin(x/2)/x

。回想一下，sin(x/2)/x 接近 0 时的极限是 1/2，因此该图在 y=0.5 处应该是一条几乎笔直的水平线，中间有一个缺失值。

（缺少 X 轴，因为当您尝试绘制极小的 X 值时，matplotlib 会崩溃。）

数字 5e-324 是一个极端情况，因为 5e-324 / 2 不能表示为浮点数，但它后面的数字也不是很大。直到 x=2.5e-321 左右，计算

np.sin(x/2)/x

才始终精确到小数点后三位。

检查 0.0 作为特殊情况的代码比使用具有任意阈值的截断泰勒级数更清晰。

我不同意 - 使用浮点数学并检查与零完全相等的代码几乎总是代码味道。通常，如果当输入恰好为零时代码被零除，那么它在接近零的值时在数值上不稳定。