2.4 定点数开方与超越函数实现方法

计算机在实现不能通过有限次的四则运算计算的函数时，会用多项式/有理式拟合法、迭代法和查表法。因为定点数的除法性能消耗较大，所以这里我们主要使用多项式拟合。当然，也会探讨其他方法的优劣。

2.4.1 多项式拟合

在区间[a,b] 上可以用n次多项式P n (x)近似求解我们要求解的函数f(x)。多项式P n (x)和要求解的函数f(x)的误差绝对值也是一个函数|Pn(x)-f(x)|。多项式P n (x)有多种选择，在泰勒展开式中截取前n+1项：

就是一个例子。泰勒展开式的误差是

可以看出，x如果离a很远的话，误差也会很大。

本书的主编叶劲峰提出了下面这种方法，误差绝对值函数在区间[a, b]上有最大值：

对于给定的n，我们的目标就是求出Pn的系数，使得误差绝对值函数的最大值最小：

这里使用数学软件Mathematica求拟合多项式的系数。在Mathematica中依次输入如下命令：

    ＜＜ FunctionApproximations` // 表示导入函数拟合工具包
    CoefficientList[
        MiniMaxApproximation[Exp[x], {x, {1, 2}, 4, 0}][[2, 1]], x]

其中，{x, {1, 2}, 4, 0}表示在区间[1,2]上逼近，用分子4次多项式、分母0次多项式的有理式进行逼近。得到的结果为

    {1.17974, 0.380661, 1.32348, -0.350357, 0.184801}

意思是在区间[1,2]上，

    Floor[CoefficientList[
        MiniMaxApproximation[Exp[x], {x, {1, 2}, 4, 0}][[2, 1]], x]＊2^32]

乘以232再取整就是为了把系数表示为定点数形式。具体的其他用法请参考软件的帮助文档。

需要再次强调的是，使用MiniMaxApproximation求解出的多项式都只是在指定区间上接近我们要求解的函数。所以，如果要求的点落在这个区间以外，那么就要利用函数的性质转换成这个区间内的点。

2.4.2 正弦/余弦函数

我们可以利用正弦/余弦函数的周期性、对称性等，把一般的角度x转换成小区间范围内的角度进行计算：

我们可以让k1对8取余：

那么角度的终边就落在[kπ/4, (k+1)π/4]范围内。换言之，我们把平面均分成八块，需要确认角度的终边落在哪一块。接下来的思路就是求解终边和最近的坐标轴所在直线的夹角的正弦值和余弦值，再转换为我们要求解的角度的正弦值和余弦值。注意：这里强调的是和坐标轴所在直线的夹角，而非和坐标轴的夹角，因为和坐标轴的夹角一般指的是和坐标轴正向所成的角度，而和坐标轴所在直线的夹角指的是小于或等于直角的那个角度。这里记角的终边和最近的坐标轴所在直线的夹角为θ。如果k是偶数，那么θ=x 1；如果k是奇数，那么θ=π/4-x 1。我们利用拟合多项式求出θ/2的正弦值和余弦值，再利用二倍角公式得到sinθ, cosθ：

然后根据角的终边所在的位置，我们就可以得到所要求的三角函数和sinθ, cosθ之间的关系。比如k1=5，角的终边落在[5π/4,3π/2]范围内，那么就有

如图2.1所示是定点数与编译器自带的正弦函数的相对误差（这里x的范围非常大，所以图中采用了对数刻度）。

图2.1 定点数与编译器自带的正弦函数的相对误差

2.4.3 指数函数

指数函数增长非常迅速，幸好我们可以把输入的x值拆分成整数部分和小数部分：

其中，是e的整数次方，很容易求解。而{x}属于[0,1]区间，在这个区间范围内指数函数可以方便地用多项式拟合。

如图2.2所示是定点数与编译器自带的指数函数的相对误差，使用这种方法求出Exp的误差对比（这里x的范围比较小，图中使用了均匀刻度）。

图2.2 定点数与编译器自带的指数函数的相对误差

2.4.4 对数函数

因为这里的定点数是基于整数的二进制表示的，所以求解以2为底的对数会比较方便，一般的对数也可以使用换底公式进行转换：

x可以约化为

那么

因为这里的定点数是基于整数的二进制实现的，所以知道最高位的1的位置就可以求出n。余下的部分可以用拟合多项式求解。

2.4.5 开方运算

在3D运算中，有大量的求向量长度的运算，其中会用到开方。开方有很多方法，如下所示。

（1）牛顿迭代法：后面一项可以通过前面一项简单地运算得出。对于定点数运算来说，这里有除法，不是理想的计算方式。

（2）借用浮点数开方运算：已知定点数y，我们要求定点数x使得x2=y。根据之前对定点数的运算和它对应的整数运算的关系，有：

其中，f(x),f(y)分别是定点数x,y 对应的整数，我们可以先把整数2nf(y)转换为浮点数开根号再取整，那么f(x)只能取附近一定范围内的整数。我们选取平方后和2nf(y)最接近的一项，在大部分情况下比较两项即可（数学上说只能是这两项之一，但浮点数开根号运算也不是完全精确的）。

（3）整数开方法：可以根据之前论述的定点数的值和它对应的整数值之间的关系，对其对应的整数开方。整数开方可以参考Jack W. Crenshaw的论文Integer Square Roots [1]，把整数视为一个2的幂次多项式，使用类似于长除法的手段，过程类似于手算开根号（日常的手算开根号方法是基于整数的十进制表示的），时间复杂度直接和整数位数相关。

（4）本章提倡的MinMax多项式拟合法。

2.4.6 开方求倒数

在3D运算中，有大量的向量单位化的运算，其中会用到开方求倒数。以向量V(x, y, z)为例：

将其归一化（单位化）的方法是V 乘以向量x, y, z分量平方和的开方倒数。引擎通常提供一个函数InvSqrt()。还记得卡马克实现的版本吗？先求一个近似值，然后通过牛顿迭代法再计算一次，这样得到一个相对精确的值。这里我们先用4次多项式模拟找到一个较精确的值（最大误差为0.001%），如图2.3所示。

图2.3 InvSqrt：多项式拟合

然后通过牛顿迭代法再计算一次（相对误差能缩小到0.00001%以内），如图2.4所示。

图2.4 InvSqrt：多项式拟合+牛顿迭代一次

2.4.7 为什么不用查表法

需要注意的是，使用查表法时，表中的数值是整数，最好将其放在代码中，或者通过配置加载，但不能通过初始时的浮点数运算再转换得来。只有在没有较好的多项式拟合方法时才会选择使用查表法。查表法的优点是实现简单，缺点有两个：

（1）非常依赖定点数的表示，一旦对定点数的表示格式做出微小调整，表中所有数据都要更新。

（2）在进行密集的数学运算时，查表法对高速缓存利用率高；对于其他常用情况，则容易遇到缓存命中失败的情况。