By Long Luo
之前写过一篇介绍 CORDIC 算法 1 的文章，里面提到 CORDIC 算法的 不足 之处，CORDIC 算法的输入角度范围需要在 \([−99.88^{\circ} , 99.88^{\circ}]\) ，那么我们不禁要问，如果输入角度 \(\large {\theta }\) 很大的话，怎么处理呢？
这个问题同样存在于 泰勒展开式(Taylor series) 2 中，比如 \(\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}\) 种方法都用到了参数归约的思想，你可能会感到震惊，什么？这种小学计算题也用到了参数归约算法吗？
上一章计算 \(\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}\) 个步骤进行计算：
现在回到上一节的小学数学计算题，我们实际上用到了 \(\large {2}\) 种参数归约：
\[ exp(x + y) = \exp(x) \exp(y) \]
\[ \log (xy) = \log (x) + \log (y) \]
\[ \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) 3 , 其具体实现可参考这篇文章： 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}\) 种参数归约：
这种归约可以应用在 \(\large {f(x)}\) 是周期函数的情况，比如三角函数，此时 \(\large {C = 2 \pi}\) ；也可以应用于其他函数，比如小学数学我们知道计算 \(\large { \frac {a}{b}}\) 就是看有多少个 \(\large {b}\) 相加小于等于 \(\large {a}\) ，具体可参考这篇文章：29. Divide Two Integers 。
应用于计算指数函数 \(\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)}\) 。
通过上面的分析，现在让我们去计算任意输入 \(\large {x}\) 的 \(\large { \sin (x)}\) 、 \(\large {\cos (x)}\) 的值，可以分为下面 \(\large {2}\) 种情况：
听起来似乎很简单，但事实上远远没有这么容易！
我们的电脑是基于 二进制(Binary) 4 的，本质只是高电平和低电平在电路上切换运行而已。因为 CPU 种的 逻辑运算单元(Arithmetic logic unit) 5 只能做加法和移位操作，因此而诞生了 计算机算术(Computer Arithmetic) 6 这门学科！
数学中有一门学科 数值分析(Numerical Analysis) 7 就是专门研究各种计算的！
虽然三角函数的周棋是 \(\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) 8 逼近，当然里面用了很多技巧，大家可以看看库函数的具体实现即可解惑（这里先挖个坑，等我彻底看懂了再来这里填坑！）。
那对于 \(\large {x > \frac {\pi}{2}}\) ，如何计算呢？
我们可以使用下列公式将 \(\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 归约算法9 。
但是如果你认为这样就高枕无忧了的话，就太早了！
假如输入 \(\large {x = 1000001}\) 的话，上面的方法就会失效！想一想为什么？
上一章提出了一个问题，有效数字 超过 \(\large {15}\) 位的超大数字该如何计算呢？针对这个问题， Payne 与 Hanek 10 提出把浮点运算转换为大整数运算，来解决超大数字的浮点数归约问题。
要弄懂 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}\) 是超越数，是无法用二进制表示的，在计算机里只能去近似。我们最终要求得的三角函数的误差取决于下面几个方面：
对于输入参数 \(\large {x}\) 不是很大的情况，误差主要由参数归约时产生的误差决定，而当输入参数 \(\large {x}\) 很大的情况，参数归约产生的误差就不再是主要因素了！
由之前的推导，我们知道：
\[ k = \left \lfloor y \right \rfloor = x \cdot (\frac {2}{\pi }) \]
但是由于浮点数的精度限制，我们知道对于 \(\large {x}\) 很大情况，我们不能直接去计算！
由三角函数关系可知，我们实际上只需要计算 \(\large {k \% 4}\) 的值即可，也就是说只需要知道 \(\large {k}\) 的最后 \(\large {2}\) 个 二进制位值即可，这样就可以节省大量运算了！
让我们回到 浮点数标准 11 ，以 \(\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) 算法的理解，后续有新的发现、感悟都会更新此文章。
计算机算术↩︎
W. Cody and W. Waite, Software Manual for the Elementary Functions, Prentice-Hall, Englewood Cliffs, N.J., 1980.↩︎
M. Payne and R. Hanek, “Radian Reduction for Trigonometric Functions”, Signum, p19-24, Jan 1983.↩︎
IEEE 754↩︎
此内容由惯性聚合(RSS阅读器)自动聚合整理，仅供阅读参考。 原文来自 — 版权归原作者所有。