近0端多项式逼近#

sin(x)的泰勒展开：

cos(x)的泰勒展开：

当然，因为展开项数限制和数的精度限制，x的范围只在靠近0的时候才有比较高的精度。

float sinf_poly (float x) {
    float r, t, s = x * x;
    r =         2.86567956e-6f; //  0x1.80a000p-19
    r = r * s - 1.98559923e-4f; // -0x1.a0690cp-13
    r = r * s + 8.33338592e-3f; //  0x1.111182p-07
    r = r * s - 1.66666672e-1f; // -0x1.555556p-03
    t = x * s + 0.0f; // ensure -0 is passed through
    r = r * t + x;
    return r;
}

float cosf_poly (float x) {
    float r, s = x * x;
    r =         2.44677067e-5f;
    r = r * s - 1.38877297e-3f;
    r = r * s + 4.16666567e-2f;
    r = r * s - 5.00000000e-1f;
    r = r * s + 1.00000000e+0f;
    return r;
}

因此，需要进行Reduction，这是这一篇的重点。

Reduction#

Cody-Waite Reduction#

因为上面的多项式要在靠近0的时候才会有比较高的精度，因此，就需要把x归约到0附近。理论上来说，我们实际上只需要把x归约到$[-\frac{\pi}{2}, \frac{\pi}{2}]$附近即可。数学上的方法为：

例如，当x = 16时，有如下的计算过程：

注意这里$\lfloor \frac{x}{\pi/2} \rfloor = 10$，这也就意味着我们得到的$x'$与$x$有：

因此，将x'代入到sin的多项式中，再取负号即可。

再举一个例子，当x=18时，我们有$\lfloor \frac{x}{\pi/2} \rfloor = 11$，此时$x'$与$x$有：

因此此时，是把x'代入到cos的多项式中，再取负号。

按照上面的理论，可以有代码：

float j; int q;
j = a * 6.36619747e-1f + 1.2582912e+7f;  // 6.36619747e-1即为2/pi的近似值
j = j - 1.25829120e+7f;
r = a - j * 1.57078552e+00f;
r = r - j * 1.08043314e-05f;
r = r - j * 2.56334407e-12f;
q = (int)j;

不过可惜的是，如果给的x过大的话，上面的代码会失效。其中的一个重点问题是通过单精度浮点数0x3f22f983表示的$2/\pi$精度是不够的，假设这里的x的值是100000000，那么算出来的j的结果是63661976，但是实际上，如果换用double的话，可以算出j的实际值为63661977，这样一来就造成了偏差。因此，对于大数而言，我们必须另外想办法。

Payne-Hanek Reduction#

对于比较大的数，Cody-Waite的算法就会失效，主要的问题在于Cody-Waite算法只使用到固定精度的浮点数，对于超大的数精度就会不够。1983年Payne与Hanek提出把浮点运算转换为大整数运算，这样来解决比较大型的浮点数的归约问题。但相对于Cody-Waite算法来说，Payne-Hanek算法相当不好理解。需要一些数学上的推理才可以表述清楚。

Mathematical representation of Cody-Waite#

对于给定的数$x$我们希望把$x$给归约到$C$中，那么我们使用的公式为：

Mathematical representation of Payne-Hanek#

Payne与Hanek的论文中主要是针对三角函数的归约问题，这一点很重要，如果不是三角函数的话是不能用这个算法的，也就是说，现在给定$x$我们希望计算：

对于IEEE754下的浮点数$x$，我们规定$x$的阶码$e_x$，尾数位$p$，那么$x$的尾数$X$与$x$有如下的数学关系：

or:

例如对于一个float数384，$1.1_2\times 2^8$，十进制表示为384.0，十六进制表示为0x43c00000，它的尾数是0x400000，但是根据IEE754标准，我们知道要填一个1，所以它的尾数就要变成0xc00000，也就是12582912。如果直接运算的话就是：

现在我们把(2)给代入到(1)当中，可以得到：

然后把$C=\pi/2$代入，把(3)式子变形成：

接下来的问题是我们怎么去处理这个$2/\pi \times 2^{e_x - p}$，我们把$2/\pi$的200位二进制写在下面：

0.10100010111110001100000110110111001001110010001000001010
100101001111111000010011101010111110100011111010100110100
110111011100000011011011011000101001010110011001001111000
100001110010000010000011111111001010001011...

我们把这个$2/\pi$，看成如下的形式：

这样一来我们就有如下的数学表达式：

然后我们把(5)式再代入带(4)当中：

我们考察其中连加式的第i项，是$v_i 2^{e_x-p-i-1}$，这里的$v_i$要么是0要么是1，当然我们不去理会它是0的情况，只考虑它是1，那么这一项就变成了$2^{e_x-p-i-1}\pi$。注意回顾一下我们本来的目的，我们的目的是要把一个大数$x$进行归约。但是要针对三角函数。三角函数都是以$2\pi$为周期的周期函数。那么如果$2^{e_x-p-i-1}\pi$是$2\pi$的倍数，自然也就不需要计算了，换句话说，对于$i\leq e_x-p-1$的项而言，我们是可以不必计算的，只需要计算$i\ge e_x-p-1$之后的项即可。

我们先假定$m > 0$，大到这个连加式的前几项$m-i$都是整数，那么把这个式子变成：

这个式子中的第一个连加项必定是$2\pi$的倍数，因此我们直接用一个$2N\pi$来代替就好，仅仅考虑后面的项：

对于中间的项，实际当中我们不可能加到无穷项，假设我们只加到$m+l$项，那么我们做如下的处理：

现在我们来关注一个问题，这个公式的含义：

实际上代表了$\pi/2$从第m位开始，到第$m+l$位为止的二进制整数。这个整数乘以$X$，再乘以$2^{-(m+l)}$，可以得到一个数，让这个$k$取这个数的整数部分，那么这个数的小数部分乘以$\pi/2$就是我们要求的$y$值了。

Example：#

求$x=1.1_2\times 2^{30}$经过Paney-Hanek算法归约后的数。经过上述的论证，我们可以得到以下的一些变量：

因此我们从$\pi/2$的第6个小数位开始取，连续取64位，取的位数必须超过p的值，否则最后算的时候会发现没有小数部分。取得位数越多越精确。这里是1011111000110000011011011100100111001000100000101010010100111111，它的十进制值是13704574379508278591，然后我们可以得到一个$G$

然后计算$G\times X$，得到37392713.36464738，我们只要令k=37392713，那么小数部分0.36464738乘以$\pi/2$，就可以得到0.5727867650793，这也就是我们要求的y值。

不过上面的论述有一点不足，就是上面我们是假定$m>0$的，有的时候$x$也可能没有达到一定程度，因而$m$也可能小于0，我们重新来论述一下：

这里面的$m$小于0，将其用$-n$来替代：

上面的公式表明我们取数的时候要从$2/\pi$的第一个小数点就要开始取了，假定我们往后取$l$项：

这里面的$\sum_{i=1}^l v_i 2^{l-i}$就代表着$2/\pi$从第1个小数位到第$l$个小数位的整数，然后这个整数要除以$2^{-(n+l)}$，乘以$X$后得到的数跟先前类似，得到的小数部分与$\pi/2$相乘即可得到最终的归约值。

Example#

来看先前$x=384$的例子：

让$l=48$，这样$n+l=64$，取出的数为10100010111110001100000110110111001001110010001000001010，十进制表示为179189285594914，然后运算：

注意这里是63。

然后可以得到小数部分为：0.457558794，乘以$\pi/2$得到0.718731673，非常接近我们要计算的值。