提问人：Corwin of Amber 提问时间：2/21/2021 更新时间：10/14/2022 访问量：1755
如何在 Python 3.7 中避免 f（x） = （1-cos（x））/x**2 中小数的灾难性取消？
How do I avoid catastrophic cancellation for small numbers in f(x) = (1-cos(x))/x**2 in Python 3.7?
问：
如何避免因小数字而灾难性取消 Python 3.7 中的 f（x） = （1-cos（x））/x**2？
这就是我到目前为止尝试过的（我知道，关键是一些三角恒等式，它使您能够避免取消，而且我也知道，在使用了 L'Hopital 规则后，f（x） 的极限→0 是 0.5，所以正确的程序输出非常接近 0.5，例如，如果您使用 x = 1.2e-4，这就是您得到的， 但是你会用较小的数字（如 1.2e-8）取消，我需要这样做，这样就不会发生这种情况）。
from math import *
def f(x): #these are all the same function using different identities
a = (1-(sin(x)/tan(x)))/(x**2)
b = (1-(sin(2*x)/(2*sin(x))))/(x**2)
c = (1-((1-((tan(x/2))**2))/(1+(tan(x/2))**2)))/(x**2)
d = (sin(x)**2+cos(x)**2-cos(x))/(x**2)
e = (sin(x)**2+cos(x)**2-(sin(2*x)/(2*sin(x))))/(x**2)
return a, b, c, d, e
print(k(1.2e-8))
#Output: (0.0, 0.7709882115452477, 0.0, 0.0, 0.0) - whereas need 0.5000...
答：
4赞
harold
2/21/2021
#1
喜欢这个：
sin(x)/x * tan(x/2)/x
它一直完成工作到最后，仍然没问题。x = 1e-308
不幸的是，我无法提供太多关于为什么它运行良好的见解。
评论
0赞
Tim
2/21/2021
可能与 vs 的实现有关。请注意，虽然 .同样地。sincosinesin(1e-100)=1e-100cos(1e-100)=1.0tan(1e-100)=1e-100
0赞
Tim
2/21/2021
我猜他们在很小的时候使用了近似值。当接近 0 时，您可以使用接近 1 的 L'Hopital 进行显示sin(x)=xxsin(x)/xx
1赞
Stefan
2/21/2021
你可以用 numpy.sinc（x/numpy.pi） 替换 sin（x）/x - 当 x 接近 0 时，sinc 的极限为 1。Tan（x）/x 也可以用 sinc 表示
0赞
Corwin of Amber
2/22/2021
完善！这太棒了，非常感谢@harold、Tim 和 Stefan。您的意见非常有帮助，值得赞赏。
1赞
soegaard
2/22/2021
使用 Herbie！herbie.uwplse.org/demo/......
0赞
hivert
2/21/2021
#2
问题在于浮点和双精度的有限。您必须使用更精确的算术，例如 mpfr。它可以在 Python 中通过绑定使用，例如
https://pypi.org/project/gmpy2/
下面是一个示例，我通过一个名为 Sagemath 的更高级的环境使用它： 我使用的是 100 位精度：
sage: R = RealField(100) // 100 bits of precision
sage: def f(x): #these are all the same function using different identities
....: a = (1-(sin(x)/tan(x)))/(x**2)
....: b = (1-(sin(2*x)/(2*sin(x))))/(x**2)
....: c = (1-((1-((tan(x/2))**2))/(1+(tan(x/2))**2)))/(x**2)
....: d = (sin(x)**2+cos(x)**2-cos(x))/(x**2)
....: e = (sin(x)**2+cos(x)**2-(sin(2*x)/(2*sin(x))))/(x**2)
....: return a, b, c, d, e
....:
sage: f(R(1.2e-8))
(0.50000000000000264647624775223,
0.49999999999999716827551705076,
0.49999999999999716827551705076,
0.49999999999999716827551705076,
0.49999999999999169007478634929)
评论
0赞
hivert
2/21/2021
礼貌地说反对时有什么问题！
0赞
Corwin of Amber
2/22/2021
非常感谢您的回答！（我没有投反对票，投了赞成票）
0赞
Eric Postpischil
2/28/2021
“使用更精确”并不是浮点型的一般好答案。人们应该了解这些问题。威廉·卡汉（William Kahan）在1980年左右将以下例子归功于让-米歇尔·穆勒（Jean-Michel Muller）。设 f（y， z） = 108−（815−1500/z）/y。给定 x[0] = 4 且 x[1] = 41/4，定义 x[n+1] = f（x[n]， x[n-1]）。计算 x[80]。正确答案在5附近。所有数值实现都是错误的;大多数报告100。请参阅链接的论文，了解原因。OP问题的一个很好的答案是了解误差的来源，并设计一个避免或减少误差的计算方法。
0赞
Eric Postpischil
2/28/2021
需要明确的是，“所有数字实现都是错误的”的含义是，使用更多的位将无济于事。更精确地使用算术不会得到正确的答案。这个例子当然是为此而精心设计的，使用更高的精度通常有助于例行计算。但有时它只是把罐子踢得更远。
0赞
hivert
3/4/2021
@EricPostpischil ： 感谢您的指点。我尝试了精度为 1000 位的计数器示例，得到了 4.999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999所以它比这更复杂。虽然我同意你的观点，但我认为这个问题是一个实现问题，而不是一个数学/数值分析问题，这超出了 SO 的范围。
2赞
Tim
2/21/2021
#3
请改用。 是数据分析和更复杂数学的标准。可以在终端中使用以下命令安装它。numpy.float128numpy
pip install numpy
from numpy import *
def f(x): #these are all the same function using different identities
a = (1-(sin(x)/tan(x)))/(x**2)
b = (1-(sin(2*x)/(2*sin(x))))/(x**2)
c = (1-((1-((tan(x/2))**2))/(1+(tan(x/2))**2)))/(x**2)
d = (sin(x)**2+cos(x)**2-cos(x))/(x**2)
e = (sin(x)**2+cos(x)**2-(sin(2*x)/(2*sin(x))))/(x**2)
return a, b, c, d, e
print(f(float128(1.2e-8)))
这将打印
(0.5003141275115400736, 0.49956120933620291774, 0.49993766842387149567, 0.49993766842387149567, 0.49956120933620291774)
评论
0赞
harold
2/21/2021
float128 在某些平台上等同于 float64
1赞
Tim
2/21/2021
更改为 ，因为我们仍然希望使用 trig 函数的版本。numpy.float128float128numpy
0赞
Tim
2/21/2021
不幸的是，float64 不起作用，因此如果您的平台将 float128 默认为 float64，可能不得不听从@harold的回答
1赞
Stefan
2/21/2021
#4
您可以有条件地返回小 x 的限制，如下所示：
import matplotlib.pyplot as plt
from numpy import *
epsilon=1e-8
def f(x):
if x<epsilon
return 0.5
return (1-sin(x)/tan(x))/x**2
#note: same as (1-cos(x))/x**2
x=arange(0,6,0.01)
y=vectorize(f)
plt.plot(x,y(x))
plt.show()
绘制的曲线看起来很平滑
注意：我更喜欢numpy而不是math。矢量化使得使用数组调用函数成为可能（效率不高，但易于使用）。
2赞
Floating Pontiff
10/14/2022
#5
如何在 Python 3.7 中避免 f（x） = （1-cos（x））/x**2 中小数的灾难性取消？
您可以使用半角度标识：
1 − cos（x） = 2 sin²（x/2）。
（历史上的三角函数会将其识别为正弦函数，或者是哈弗正弦函数的两倍，这很有用，并为此目的制成表格！
如果将 x 计算为 0 附近 x 的 cos（x） 的近似值，则减法 1 中的灾难性抵消可能会极大地放大该近似值中的任何误差——如果四舍五入为 1，则可能会给出一个无意义的答案，因此当它应该接近 1/2 时四舍五入为 0。cos(x)cos(x)cos(x)(1 - cos(x))/x**2
相反，如果你把它改写成
（1 − cos（x））/x² = 2 sin²（x/2）/x²
使用半角恒等式，然后计算接近零的 x，您将在条件良好的域上组合函数。
在二进制浮点中，除以 2 和乘以 2 总是精确的，除非出现溢出或下溢的情况，并且 sin 和平方都处于接近零的良好条件，因此即使在数学库的过程中进行中间舍入和误差，相对误差也相当小。2*(sin(x/2)/x)**2sin
也可以使用另一个半角度标识，
1 − cos（x） = sin（x） tan（x/2），
正如另一个答案所暗示的那样——SIN 和 Tan 都处于接近零的良好条件，因此它也会给出一个很好的近似值。 但是计算 sin 和 tan 的成本比仅仅计算 sin 要高，而且我认为至少对于接近零的输入，误差没有明显的优势。
要在没有条件的情况下处理 x = 0，如果您有一个近似于函数 sinc（x） = sin（x）/x 的过程，并且极限为 sinc（0） = 0，您也可以将其重写为：sinc
（1 − cos（x））/x² = 2 sin²（x/2）/x² = （1/2） sin²（x/2）/（x/2）² = （1/2） sinc²（x/2）。
下一个：如何计算输入的小数点？[复制]
评论
numpyfloat64