参数归约算法(Argument Reduction):如何在浮点数环境下计算超大数字的三角函数值?

By Long Luo

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

这个问题同样存在于 泰勒展开式(Taylor series) [2] 中,比如 sin(x)\large {\sin (x) }cos(x)\large {\cos (x) } 的泰勒展开式:

sin(x)=x13!x3+15!x517!x7+19!x9+o(x9)xR\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)=112!x2+14!x416!x6+18!x8+o(x8)xR\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}

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

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

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

x=x2kπkZx' = x - 2k \pi \quad k \subset \mathbb{Z}

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

从小学计算题开始

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

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

不借助计算器,计算 66600×666000\large {66600 \times 666000} 的值!

对于这道题,大家可能会列出下列算术:

66600×666000=666×666×100000=4435560000066600 \times 666000 = 666 \times 666 \times 100000 = 44355600000

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

66600×666000=1112×4×9×105=444×999×105=444×(10001)×105=4443556×105\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}

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

什么是参数归约 Argument Reduction ?

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

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

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

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

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

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

  1. 指数/对数 运算公式

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

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

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

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

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

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

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

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

如何对参数进行归约?

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

  1. 加法参数归约: x=xkC\large {x' = x - kC} ,其中 C\large {C} 是实常数, k\large {k} 是整数。

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

  1. 乘法参数归约:x=xkC\large {x' = \frac{x}{kC}},其中 C\large {C} 是实常数, k\large {k} 是整数。

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

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

数值分析 Numerical Analysis

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

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

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

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

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

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

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

那对于 x>π2\large {x > \frac {\pi}{2}} ,如何计算呢?

Cody-Waite 归约算法

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

x=xxπ2×π2x' = x - \lfloor \frac {x}{\frac {\pi}{2}} \rfloor \times \frac {\pi}{2}

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

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

假如输入 x=1000001\large {x = 1000001} 的话,上面的方法就会失效!想一想为什么?

假如你知道浮点数的话,你就知道为什么了!

按照 IEEE 754 浮点数标准 [9] 制定的 浮点数运算法则 [10], float 类型的单精度浮点数 [11] 的尾数部分有 23\large {23} 位二进制数,如下图所示:

IEEE_754_Single_Floating_Point_Format

在十进制下,大致相当于 log102236.9\large {\log_{10}{2^{23}} \approx 6.9} ,有效数字大约有 7\large {7} 位。

所以当 x=1000001\large {x = 1000001} 时,我们应该使用 double 类型的双精度浮点数 [12] ,这样才能保证结果有足够的精度

双精度浮点数的尾数部分有 52\large {52} 位,如下图所示:

IEEE_754_Double_Floating_Point_Format

在十进制中大致相当于 log1025215.95\large {\log_{10}{2^{52}} \approx 15.95} ,也就是说当 x\large {x} 有效数字是 [7,15]\large {[7, 15]} 时,我们应该使用 double 类型的双精度浮点数可以保证精度!

但这仍然有个问题,那就是 x\large {x} 有效数字 超过 15\large {15} 位,应该怎么办?

Payne-Hanek 归约算法

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

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

对于输入 x\large {x}

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

两边同乘 2π\large {\frac{2}{\pi }} ,可得:

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

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

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

即:

k=yk = \left \lfloor y \right \rfloor

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

f=ykf = y − k

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

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

回到我们的目标,我们需要知道 k\large {k} 的值 和 r\large {r} 的值!

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

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

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

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

计算 k\large {k}

由之前的推导,我们知道:

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

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

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

让我们回到 浮点数标准 [9:1] ,以 32\large {32} 位单精度浮点数为例,其值可以表示为:

x=(1)b31×2(b30b29b23)2127×(1.b22b21b0)2x = (-1)^{b_{31}} \times 2^{(b_{30}b_{29}\dots b_{23})_{2} - 127} \times (1.b_{22}b_{21}\dots b_{0})_{2}

即为:

value=(1)sign×2(E127)×(1+i=123b23i2i)\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) 算法的理解,后续有新的发现、感悟都会更新此文章。

参考文献


  1. CORDIC算法:一种高效计算三角函数值的方法 ↩︎

  2. 泰勒展开式(Taylor series) ↩︎ ↩︎

  3. Exponentiation by squaring ↩︎

  4. Binary number ↩︎

  5. Arithmetic logic unit ↩︎

  6. 计算机算术 ↩︎

  7. Numerical Analysis ↩︎

  8. W. Cody and W. Waite, Software Manual for the Elementary Functions, Prentice-Hall, Englewood Cliffs, N.J., 1980. ↩︎

  9. IEEE 754 ↩︎ ↩︎

  10. Floating-point arithmetic ↩︎

  11. Single-precision floating-point format ↩︎

  12. Double-precision floating-point format ↩︎

  13. M. Payne and R. Hanek, “Radian Reduction for Trigonometric Functions”, Signum, p19-24, Jan 1983. ↩︎