稳定的余切

Question

余切函数的实现是否比返回1.0/tan(x)更稳定;？

Answer 1

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不会精确到a的全宽度double尾数 - 事实上,如果你不幸,它可能只有几个有意义的位.

Answer 2

TL;博士编号

根据经验，在寻找不准确的来源时，应首先关注加法和减法，这可能会导致减法抵消问题。除了添加额外的舍入误差之外，乘法和除法通常对准确性无害，但可能会通过中间计算中的上溢和下溢导致问题。

没有机器号x可以接近 ?/2 的倍数以导致tan(x)溢出，因此tan(x)对于任何 IEEE-754 浮点格式的所有浮点编码都是明确定义和有限的，并且通过扩展，cot(x) = 1.0 / tan(x).

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

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

对所有float值重复这个详尽的测试cot(x) = cos(x) / sin(x)，其中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幅度足够小时的问题。这会导致非常大的错误，导致结果中没有有效位。