【发布时间】：2010-09-17 19:30:44
【问题描述】：
余切函数有比 return 1.0/tan(x); 更稳定的实现吗？
【问题讨论】：
-
我假设你对 1/tan 的问题是当函数应该返回零时它在 pi/2 处未定义？
-
在接近 pi/2 时进行反转是不可靠的...我想知道是否有更好的方法来做到这一点。
标签： c++ trigonometry
【发布时间】：2010-09-17 19:30:44
【问题描述】：
cot(x) = cos(x)/sin(x) 在接近 π/2 的数值上应该比cot(x) = 1/tan(x) 更稳定。您可以在拥有它的平台上使用sincos 有效地实现它。
另一种可能性是cot(x) = tan(M_PI_2 - x)。这应该比上面的要快（即使sincos可用），但也可能不太准确，因为M_PI_2当然只是超越数π/2的近似值，所以差值M_PI_2 - x会不准确到 double 尾数的整个宽度——事实上，如果你运气不好，它可能只有几个有意义的位。
【讨论】：
-
好的。我想这是我能做的最好的了。我还了解了 sincos，我必须承认我以前从未见过！
TL;DR 编号
根据经验，在寻找不准确的来源时，首先应该关注加法和减法，这可能会导致减法取消的问题。除了增加额外的舍入误差外，乘法和除法通常对准确性无害，但可能会在中间计算中通过溢出和下溢引起问题。
没有机器号 x 可以足够接近 π/2 的倍数而导致 tan(x) 溢出，因此对于任何 IEEE-754 的所有浮点编码，tan(x) 都是明确定义和有限的浮点格式，通过扩展，cot(x) = 1.0 / tan(x) 也是如此。
这很容易通过使用所有数字 float 编码执行详尽测试来证明，因为使用 double 进行详尽测试是不可行的，除非使用当今最大的超级计算机。
使用精确实现 tan() 且最大误差为 ~= 0.5 ulp 的数学库，我们发现计算 cot(x) = 1.0 / tan(x) 产生的最大误差小于 1.5 ulp，其中与tan() 本身是由除法的舍入误差造成的。
使用cot(x) = cos(x) / sin(x) 对所有float 值重复此详尽测试，其中sin() 和cos() 的最大误差为~= 0.5 ulp，我们发现cot() 中的最大误差更小比 2.0 ulps，所以略大。这很容易通过三个误差源来解释，而不是前面公式中的两个。
最后，cot(x) = tan (M_PI_2 - x) 在x 接近 M_PI_2 时会遇到前面提到的减法抵消问题，以及在有限精度浮点运算中，M_PI_2 - x == M_PI_2 当x 足够小的问题震级。这会导致非常大的错误，导致结果中没有有效位。
【讨论】：
如果考虑两个向量（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(θ)
- v 。 w = |v|.|w|.cos(theta)
- cot(theta) = cos(theta) / sin(theta) = (v . w) / |v x w|
【讨论】：