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）在進行密集的數學運算時，查表法對高速緩存利用率高；對于其他常用情況，則容易遇到緩存命中失敗的情況。