如何错误地计算余弦值

文章背景与核心概要

微积分教材通常会暗示计算机是通过幂级数来计算三角函数的，然而对于远离级数中心的值来说，这是一个常见的误解。例如，尝试使用标准的余弦幂级数来计算 \(\cos(200)\) 会导致巨大的数值不稳定性、上溢问题以及灾难性的相消误差。在实际工程中，精确计算三角函数需要依赖更高级的技术（如精确的区间约简），而不是单纯地进行朴素的级数展开。

本文由知名数学与软件专家 John D. Cook 撰写，深入剖析了初学者在数值计算中常犯的直觉错误。通过具体的 Python 代码示例，文章展示了为什么盲目套用幂级数会产生高达 \(10^{67}\) 的荒谬结果，并简要介绍了工业级数值计算中处理大数和小数运算的正确方向。

幂级数的幻觉

Calculus professors with no experience in numerical computing will tell students that computers calculate trig functions with power series. They don’t. I worked on the implementation of trig functions in hardware, and I can assure you we didn’t just use power series.

没有数值计算经验的微积分教授会告诉学生，计算机是用幂级数来计算三角函数的。事实并非如此。我曾参与过硬件中三角函数的实现，我可以向你保证，我们绝不仅仅使用幂级数。

Power series are an excellent way to calculate functions near the center of the series, such as computing sine for small angles. But the further you get from the center, the less useful power series are.

幂级数是计算靠近级数中心函数值的极佳方法，例如计算小角度的正弦值。但你离中心越远，幂级数的实用性就越差。

一次失败的尝试：计算 \(\cos(200)\)

Let’s suppose you want to calculate \(\cos(200)\) using the power series for cosine. The \(n\)-th term of that series is:

\[\frac{(-1)^n x^{2n}}{(2n)!}\]

This is an alternating series, and so the error in truncating the series after \(n\) terms is bounded by the size of the \(n+1\) term, if you’ve gone far enough out in the series that the terms are monotonically decreasing in absolute value.

假设你想用余弦的幂级数来计算 \(\cos(200)\)。该级数的第 \(n\) 项为：

\[\frac{(-1)^n x^{2n}}{(2n)!}\]

这是一个交错级数，因此如果在第 \(n\) 项后截断级数，其误差受限于第 \(n+1\) 项的大小——前提是你在级数中取得足够远，使得各项的绝对值呈单调递减趋势。

To calculate \(\cos(200)\) to machine precision, i.e., with an error of less than \(2^{-52}\), we’d need to sum the series up to \(n\) where:

\[\left| \frac{200^{2n+2}}{(2n + 2)!} \right| < 2^{-52}\]

Actually, that will ensure that the absolute error is small enough, but not that the relative error is small enough; if the value of \(\cos(200)\) is small, we’d need more terms. Let’s ignore that and assume we’re only concerned with absolute error.

为了以机器精度（即误差小于 \(2^{-52}\)）计算 \(\cos(200)\)，我们需要将级数求和到满足以下条件的 \(n\)：

\[\left| \frac{200^{2n+2}}{(2n + 2)!} \right| < 2^{-52}\]

实际上，这只能确保绝对误差足够小，而不能确保相对误差足够小；如果 \(\cos(200)\) 的值本身很小，我们需要更多的项。让我们暂时忽略这一点，假设我们只关心绝对误差。

It turns out we’d need 287 terms. That’s a lot of terms. But you might say, "That’s fine. I’m not in a hurry, and it’s just more work for the computer, not for me." OK, so let’s try:
from math import *

s = 0
for n in range(288):
    s += (-1)**n * 200**(2*n) / factorial(2*n)
print(s)

事实证明，我们需要 287 项。这是相当多的项。但你可能会说：“没关系。我不急，这只是计算机的工作，又不是我的。”好吧，让我们来试一下：

from math import *

s = 0
for n in range(288):
    s += (-1)**n * 200**(2*n) / factorial(2*n)
print(s)

This prints -3.6840358571084123e+67. You may suspect the answer is incorrect since values of cosine are on the order of \(1\), not on the order of \(10^{67}\). Something went spectacularly bad. On closer inspection, it’s remarkable the code didn’t crash.

这段代码打印出了 -3.6840358571084123e+67。你可能会怀疑这个答案是错误的，因为余弦值应该在 \(1\) 的量级，而不是 \(10^{67}\) 的量级。有些地方出了严重的问题。仔细检查后，代码竟然没有崩溃，这本身就很令人惊奇了。

为什么会失败：上溢与灾难性相消

If you changed 200 to 200.0 above, the code would crash. Calculating 200.0**(2*n) overflows when \(n = 67\). But when we calculate 200**(2*n), the result is an integer. And we’re dividing by factorial(2*n), which is also an integer.

如果你把上面代码中的 200 改为 200.0，代码就会崩溃。当 \(n = 67\) 时，计算 200.0**(2*n) 会发生上溢。但当我们计算 200**(2*n) 时，结果是一个整数。并且我们在除以 factorial(2*n)，这也是一个整数。

Both of these integers become too large to fit in a float, but their ratio has a maximum value of around \(10^{80}\), smaller than the maximum float, which is on the order of \(10^{308}\).

这两个整数都变得太大而无法放入浮点数（float）中，但它们的比值最大值约为 \(10^{80}\)，小于浮点数的最大上限（在 \(10^{308}\) 的量级）。

When we don’t overflow, we have a different problem: catastrophic cancellation. You can’t calculate a number between \(-1\) and \(1\) as an alternating sum of numbers as large as \(10^{80}\). You’d need more than \(80 + 16 = 96\) decimal places of precision to compute the sum accurately, and floating point only gives you between \(15\) and \(16\) decimal places of precision.

即使我们没有遇到上溢，也会面临另一个问题：灾难性相消（Catastrophic Cancellation）。你无法通过将高达 \(10^{80}\) 的数字进行交错求和，来计算出一个介于 \(-1\) 和 \(1\) 之间的数。为了准确计算这个总和，你需要超过 \(80 + 16 = 96\) 位十进制数的精度，而标准浮点数只能提供 \(15\) 到 \(16\) 位的十进制精度。

正确的方法：区间约简

So how would you calculate \(\cos(200)\)? The first step would be to use some sort of range reduction on \(200\). You could reduce \(200 \pmod{2\pi}\) to get a smaller number to work with:
>>> from math import cos, pi
>>> x = 200 % (2*pi)
>>> x
5.221255477432827

那么，究竟该如何计算 \(\cos(200)\) 呢？第一步是对 \(200\) 进行某种形式的区间约简（Range Reduction）。你可以对 \(200 \pmod{2\pi}\) 进行约简，以得到一个更小的数来进行计算：

>>> from math import cos, pi
>>> x = 200 % (2*pi)
>>> x
5.221255477432827

Using a power series to compute the cosine of \(5.221255477432827\) is feasible, but not optimal. There’s also another problem: the naive range reduction above loses some precision.
>>> cos(x)
0.48718767500701254
>>> cos(x) - cos(200)
6.661338147750939e-15

使用幂级数来计算 \(5.221255477432827\) 的余弦值是可行的，但并非最佳选择。这里还有另一个问题：上面这种朴素的区间约简丢失了一些精度。

>>> cos(x)
0.48718767500701254
>>> cos(x) - cos(200)
6.661338147750939e-15

The error is small, but it’s still an order of magnitude larger than machine precision. You can’t simply reduce \(n \pmod{2\pi}\) with ordinary float division because the integer part of \(n / 2\pi\) pushes some digits of precision off the right end. I intend to write about how range reduction works in future posts.

误差虽然很小，但仍然比机器精度大了一个数量级。你不能简单地用普通的浮点数除法来计算 \(n \pmod{2\pi}\)，因为 \(n / 2\pi\) 的整数部分会把右侧的一些精度位数挤掉。我打算在以后的文章中专门探讨区间约简的工作原理。

Source: How not to calculate cosine by John D. Cook.