稳定余切

Question

是否有比return 1.0/tan(X)更稳定的余切函数实现？

Answer 1

与cot(x) = 1/tan(x)相比，cot(x) = cos(x)/sin(x)在数值上应该更接近π/2。您可以在拥有它的平台上使用sincos高效地实现这一点。

另一种可能性是cot(x) = tan(M_PI_2 - x)。这应该比上面的更快(即使sincos可用)，但它也可能不太准确，因为M_PI_2当然只是超越数π/2的近似值，所以差M_PI_2 - x不会精确到double尾数的整个宽度--实际上，如果您不走运，它可能只有几个有意义的位。

Answer 2

TL；

作为一个经验法则，当寻找不准确的来源时，应该首先关注加法和减法，这可能会导致减法取消的问题。除了增加额外的舍入误差外，乘法和除法通常对精度无害，但在中间计算中可能会通过上溢和下溢而导致问题。

任何机器编号x都不能接近π/2的倍数而导致tan(x)溢出，因此tan(x)是定义良好的，对于任何IEEE754浮点格式的所有浮点编码都是有限的，cot(x) = 1.0 / tan(x)也是如此。

这很容易通过对所有数值float编码执行穷举测试来证明，因为使用double的穷举测试是不可行的，除非使用当今存在的最大的超级计算机。

使用具有最大误差为~= 0.5 ulp的tan()的精确实现的数学库，我们发现计算cot(x) = 1.0 / tan(x)产生的最大误差小于1.5ULP，其中与tan()本身相比的额外误差是由除法的舍入误差贡献的。

使用cot(x) = cos(x) / sin(x)对所有float值重复此详尽测试，其中计算sin()和cos()的最大误差为~= 0.5ulp，我们发现cot()中的最大误差小于2.0ULPS，因此略大。这很容易解释为有三个误差源，而不是前一个公式中的两个。

最后，cot(x) = tan (M_PI_2 - x)遇到了前面提到的当x接近M_PI_2时的减法消除问题，以及在有限精度浮点算术中，当x足够小时的M_PI_2 - x == M_PI_2问题。这可能会导致非常大的错误，使我们在结果中没有有效位。

Answer 3

如果考虑两个向量(v和w)之间的角度，还可以按如下方式获得余切线(使用Eigen::Vector3d)：

inline double cot(Eigen::Vector3d v, Eigen::Vector3d w) { 
    return( v.dot(w) / (v.cross(w).norm()) ); 
};

用θ表示v和w之间的角度，上面的函数是正确的，因为：

• |v x w| = |v|.|w|.sin(theta)
• v。W= |v|.|w|.cos(theta)
• cot(theta) =cos(θ)/sin(θ)= (v .w) / |v x w|