By Long Luo

之前写过一篇介绍 CORDIC 算法 ¹ 的文章，里面提到 CORDIC 算法的 不足 之处，CORDIC 算法的输入角度范围需要在 \([−99.88^{\circ} , 99.88^{\circ}]\) ，那么我们不禁要问，如果输入角度 \(\large {\theta }\) 很大的话，怎么处理呢？

这个问题同样存在于 泰勒展开式(Taylor series) ² 中，比如 \(\large {\sin (x) }\) 和 \(\large {\cos (x) }\) 的泰勒展开式：

\[ \sin(x) = x - \frac {1}{3!}x^3 + \frac {1}{5!}x^5 - \frac {1}{7!} x^7 + \frac {1}{9!} x^9 + o(x^9) \quad \forall x \subset \mathbb{R} \]

\[ \cos(x) = 1 - \frac {1}{2!}x^2 + \frac {1}{4!}x^4 - \frac {1}{6!} x^6 + \frac {1}{8!} x^8 + o(x^8) \quad \forall x \subset \mathbb{R} \]

虽然在整个实数集 \(\large { \mathbb{R}}\) 都成立，但是在实际应用中因为展开项数限制和浮点数的精度限制， \(\large {x}\) 的范围只有在接近 \(\large {0}\) 的时候才有比较高的精度。

但是实际应用中，如果输入 \(\large {x}\) 很大的话，比如 \(\large {2^{32}, 10^{10}, 10^{22} \dots }\) 情况下怎么得到足够精确的值呢？

中学里我们知道三角函数是周期函数，对于比较大的值，我们可以使用下面的公式将值归约到一个比较小的范围内。

\[ x' = x - 2k \pi \quad k \subset \mathbb{Z} \]

这就是我们今天要讲的 参数归约(Argument Reduction) 算法。

从小学计算题开始

参数归约 听起来就很唬人，什么是参数啊，什么归约啊，都是些高大上的名词，听起来云里雾里的！

为了不让大家产生厌倦和畏难心理，我们先从一道小学数学计算题开始：

不借助计算器，计算 \(\large {66600 \times 666000}\) 的值！

对于这道题，大家可能会列出下列算术：

\[ 66600 \times 666000 = 666 \times 666 \times 100000 = 44355600000 \]

但其实呢，我们也可以使用下面的方法：

\[ \begin{aligned} 66600 \times 666000 &= 111^2 \times 4 \times 9 \times 10^5 \\&= 444 \times 999 \times 10^5 \\&= 444 \times (1000 - 1) \times 10^5 \\&= 4443556 \times 10^5 \end{aligned} \]

如果我说上面这 \(\large {2}\) 种方法都用到了参数归约的思想，你可能会感到震惊，什么？这种小学计算题也用到了参数归约算法吗？

什么是参数归约 Argument Reduction ？

上一章计算 \(\large {66600 \times 666000}\) 时，我们将 \(\large {666 \times 666}\) 化简为 \(\large {444 \times (1000 - 1)}\) ，再在结果后面直接加上 \(\large {5}\) 个 \(\large {0}\) ，那么你有没有想过这背后隐含了什么数学思想吗？

下面我们正式进入今天的课题：参数归约(Argument Reduction) 。

为了提高数学函数的计算效率，将初始问题转变或者说缩小到函数更容易计算的域内，这就是参数归约。

已知函数 \(\large {f}\) ，求 \(\large {y = f(x)}\) 的值，可以通过以下 \(\large {3}\) 个步骤进行计算：

1. 将 \(\large {x}\) 转换为缩小的参数 \(\large {x'}\) ；
2. 计算 \(\large {y' = f(x')}\) ；
3. 使用函数恒等式从 \(\large {f(x')}\) 计算出 \(\large {f(x)}\) 。

现在回到上一节的小学数学计算题，我们实际上用到了 \(\large {2}\) 种参数归约：

1. 指数/对数 运算公式

\[ exp(x + y) = \exp(x) \exp(y) \]

\[ \log (xy) = \log (x) + \log (y) \]

2. 相加公式。不过上面小学数学题用的非常简单的分配律和结合律，实际上我们用的更复杂的公式，比如各种三角恒等式：

\[ \sin (x + y) = \sin(x) \cos (y) + \cos (x) \sin (y) \]

\[ \tan (x + y) = \frac {\tan (x) + \tan (y)}{1 - \tan (x) \tan (y)} \]

实际上为了让幂级数更快地收敛，通常我们取 \(\large {x = y}\) 以获得双倍公式，比如 \(\large {e^ {2x} = (e^x)^2}\) ，比如 快速幂算法(Exponentiation by squaring) ³ , 其具体实现可参考这篇文章： Fast Power Algorithm: Binary Exponentiation 。

而计算器中也常用到三倍角公式 \(\large {\sin (3x) = 3 \sin (x) - 4 \sin ^3(x)}\) 去计算三角函数值，具体可参考这个视频： 计算器是如何计算出三角函数和对数的？ 。

可能有同学会问，那二倍角公式 \(\large {\sin (2x) = 2 \sin(x) \cos (x)}\) 就不用了吗？这个谜底待后续章节介绍。

如何对参数进行归约？

这一章我们来讲如何进行参数归约，通常我们区分 \(\large {2}\) 种参数归约：

1. 加法参数归约： \(\large {x' = x - kC}\) ，其中 \(\large {C}\) 是实常数， \(\large {k}\) 是整数。

这种归约可以应用在 \(\large {f(x)}\) 是周期函数的情况，比如三角函数，此时 \(\large {C = 2 \pi}\) ；也可以应用于其他函数，比如小学数学我们知道计算 \(\large { \frac {a}{b}}\) 就是看有多少个 \(\large {b}\) 相加小于等于 \(\large {a}\) ，具体可参考这篇文章：29. Divide Two Integers 。

2. 乘法参数归约：\(\large {x' = \frac{x}{kC}}\)，其中 \(\large {C}\) 是实常数， \(\large {k}\) 是整数。

应用于计算指数函数 \(\large { \exp(x)}\) 时，其中 \(\large {C = 2}\) 。

值得注意的是，对于给定的函数，两种参数归约方式都可能使用。例如，对于 \(\large {\sin (x) }\) ，我们既可以使用三倍角公式 \(\large {\sin (3x) = 3 \sin (x) - 4 \sin^3 (x)}\) 化简，也可以使用加法归约 \(\large {\sin (x + 2 k \pi) = \sin (x)}\) 。

数值分析 Numerical Analysis

通过上面的分析，现在让我们去计算任意输入 \(\large {x}\) 的 \(\large { \sin (x)}\) 、 \(\large {\cos (x)}\) 的值，可以分为下面 \(\large {2}\) 种情况：

1. \(\large {0 < x \leq \frac {\pi}{2}}\) ，使用泰勒展开或者 CORDIC 算法；
2. \(\large {x > \frac {\pi}{2}}\) ，先将 \(\large {x}\) 归约到 \(\large {x' = x + k \frac {\pi}{2}}\) ，再回到第一步计算。

听起来似乎很简单，但事实上远远没有这么容易！

我们的电脑是基于 二进制(Binary) ⁴ 的，本质只是高电平和低电平在电路上切换运行而已。因为 CPU 种的 逻辑运算单元(Arithmetic logic unit) ⁵ 只能做加法和移位操作，因此而诞生了 计算机算术(Computer Arithmetic) ⁶ 这门学科！

数学中有一门学科 数值分析(Numerical Analysis) ⁷ 就是专门研究各种计算的！

虽然三角函数的周棋是 \(\large {2 \pi}\) ，但实际上我们只用归约到 \(\large {[-\frac {\pi}{4},\frac {\pi}{4}]}\) 即可，这里大家可以想想为什么？

之前我以为数值运算对于 \(\large {[-\frac {\pi}{2},\frac {\pi}{2}]}\) 的参数，会使用 CORDIC 算法，但实际上我看了一些数值计算库，发现对于 \(\large {[-\frac {\pi}{2},\frac {\pi}{2}]}\) 还是使用泰勒(Taylor Series) ⁸ 逼近，当然里面用了很多技巧，大家可以看看库函数的具体实现即可解惑（这里先挖个坑，等我彻底看懂了再来这里填坑！）。

那对于 \(\large {x > \frac {\pi}{2}}\) ，如何计算呢？

Cody-Waite 归约算法

我们可以使用下列公式将 \(\large {x}\) 归约到 \(\large {[-\frac {\pi}{4},\frac {\pi}{4}]}\) ：

\[ x' = x - \lfloor \frac {x}{\frac {\pi}{2}} \rfloor \times \frac {\pi}{2} \]

我们可以很容易按照上述思想写出对应的代码，这就是 Cody & Waite 提出的 Cody-Waite 归约算法⁹ 。

但是如果你认为这样就高枕无忧了的话，就太早了！

假如输入 \(\large {x = 1000001}\) 的话，上面的方法就会失效！想一想为什么？

Payne-Hanek 归约算法

上一章提出了一个问题，有效数字 超过 \(\large {15}\) 位的超大数字该如何计算呢？针对这个问题， Payne 与 Hanek ¹⁰ 提出把浮点运算转换为大整数运算，来解决超大数字的浮点数归约问题。

要弄懂 Payne-Hanek 归约算法，需要对数学有比较深的理解，下面一步一步来分析！

对于输入 \(\large {x}\) ：

\[ x = k \cdot (\frac {\pi}{2}) + r \quad k \subset \mathbb{Z}, r \subset [-\frac {\pi}{4},\frac {\pi}{4}] \]

两边同乘 \(\large {\frac{2}{\pi }}\) ，可得：

\[ x \cdot ( \frac{2}{\pi }) = k + r \cdot (\frac{2}{\pi }) \]

因为 \(\left| r \right | \leq \frac {\pi}{4}\) ， \(r \cdot (\frac {2}{\pi }) \leq 0.5\) ，也就是说：

\[ y = x \cdot (\frac {2}{\pi }) \]

即：

\[ k = \left \lfloor y \right \rfloor \]

那么所求浮点数的尾数部分：

\[ f = y − k \]

最终可得到归约之后的结果 \(\large {r}\) ：

\[ r = f \cdot (\frac {\pi }{2}) \]

回到我们的目标，我们需要知道 \(\large {k}\) 的值 和 \(\large {r}\) 的值！

那我们能直接用上述公式计算吗？

我们知道 \(\large {\pi}\) 是超越数，是无法用二进制表示的，在计算机里只能去近似。我们最终要求得的三角函数的误差取决于下面几个方面：

1. 使用多少位数的 \(\large {\pi}\) 近似值；
2. 参数归约时产生的误差；
3. 计算参数归约之后的三角函数时的误差。

对于输入参数 \(\large {x}\) 不是很大的情况，误差主要由参数归约时产生的误差决定，而当输入参数 \(\large {x}\) 很大的情况，参数归约产生的误差就不再是主要因素了！

计算 \(\large {k}\)

由之前的推导，我们知道：

\[ k = \left \lfloor y \right \rfloor = x \cdot (\frac {2}{\pi }) \]

但是由于浮点数的精度限制，我们知道对于 \(\large {x}\) 很大情况，我们不能直接去计算！

由三角函数关系可知，我们实际上只需要计算 \(\large {k \% 4}\) 的值即可，也就是说只需要知道 \(\large {k}\) 的最后 \(\large {2}\) 个 二进制位值即可，这样就可以节省大量运算了！

让我们回到 浮点数标准 ¹¹ ，以 \(\large {32}\) 位单精度浮点数为例，其值可以表示为：

\[ x = (-1)^{b_{31}} \times 2^{(b_{30}b_{29}\dots b_{23})_{2} - 127} \times (1.b_{22}b_{21}\dots b_{0})_{2} \]

即为：

\[ \text {value} = (-1)^{\text {sign}} \times 2^{(E - 127)} \times \left (1 + \sum _{i=1}^{23}b_{23-i} 2^{-i} \right) \]

(原始论文和数值分析具体实现代码太难看懂了，这篇文章写了快 1 个月了！:-( )

小结

这是目前我对参数归约(Argument Reduction) 算法的理解，后续有新的发现、感悟都会更新此文章。

参考文献