快速三角函數算法的偏差控制與評估（sin cos）

　　工程應用常涉及三角函數的快速計算。設計者每每須要下降運算精度以提升程序的運行速度。經常使用的快速三角函數算法主要包括CORDIC、泰勒展開式逼近、查表等等。然而，網上的文章大多隻介紹如何實現相應的算法，而忽視了定量分析算法精度並以此指導設計的過程。此外，算法不一樣，偏差分析的數學方法也不盡相同。在多種算法之間比較時，須要耗費一些時間來創建模型。若是創建一個通用的偏差分析框架，框架不侷限於現存的幾類算法，而且實現自動化的篩選，則有可能提升工做效率。

　　本文試圖從開發者的角度，介紹一種通用的數值型偏差分析框架，並以泰勒展開算法爲例分析了其展開式項數與精度之間的定量關係，最終試圖經過這些工做爲快速三角函數算法的設計提供必定參考。

 

1、快速三角函數算法的原理分析

1. 預處理

　　對三角函數求值應充分考慮函數自己的性質。首先是週期性：正餘弦函數的週期爲2π，自變量的全部取值都可變換到單個週期[0, 2π]內，而使因變量保持不變；其次是對稱性：例如正弦函數[0, π]上的圖像與[π, 2π]上的圖像關於X軸對稱。取其中一支[0, π]，還能夠發現[0, π/2]與[π/2, π]上的圖像關於平行Y軸的直線x = π/2對稱。同理餘弦函數也知足相似的性質。

　　綜合上述性質可知，對於某個三角函數，若已知其在[0, π/2]的函數關係及必要的對稱關係，則可經過恆等變換表示該三角函數在整個定義域上的取值。

　　設待求三角函數爲f(x)，f(x)在[0, π/2]的函數關係爲g(x)。上述求值方法可表示爲

　　其中，函數T(t)表示將t變換到單個週期內；函數W(t,y)表示根據t取值，對三角函數的因變量y進行正負變換（即關於X軸對稱的象限變換）；函數Q(t)表示對t進行關於直線t = π/2的對稱變換。

　　以f(x)=sin x爲例說明，W(y,x)、Q(x)、T(x)三個變換函數的定義以下：

　　下圖具體地描述了各變換函數的做用。

 

 

2. g(x)

　　算法的關鍵是經過g(x)來擬合f(x)在[0, π/2]的值。

　　g(x)的具體形式能夠任意選擇，例如選擇CORDIC、泰勒展開、查表等等。g(x)是決定算法的精度與速度的核心因素。

　　以f(x) = sinx爲例，取g(x)爲泰勒4階展開的結果：

　　也能夠取g(x)爲泰勒6階展開的結果：

 

二. 偏差評估模型

 　　因爲預處理機制的存在，咱們只需討論g(x)在[0,π/2]上對f(x)的擬合狀況。爲了便於數值分析，先將自變量x離散化：把X軸上的大區間[0,π/2]均分爲 N=π/2ΔX 個區間，取每一個區間的左端點表明整個區間，這樣咱們便在X軸上得到了 N 個點，每一個點座標爲xi。記每一個點處誤差爲Di = g(xi) - f(xi) ，則標準誤差S表示爲：

　　引入另外一個量Dmax = max{ |Di| }表示峯值誤差。

　　這樣，標準誤差S描述了g(x)的總體偏差，而峯值誤差Dmax則描述了局部偏差，綜合這兩個量便可對算法的偏差進行評估。同時，上述方法與g(x)的具體形式無關，即不管採用何種擬合算法，該方法都適用。

　　

　　仍以f(x) = sinx爲例。令g(x) 爲sinx的泰勒展開式。偏差分析以下。

　　N = 10；泰勒展開階數 = 4階

 　　（注：紫線爲x=π/2。只對區間[0,π/2]進行了偏差統計，下同）

 

　　N = 10；泰勒展開階數 = 6階

 

　　隨着展開式階數增長，標準誤差S和峯值誤差Dmax都呈近似指數減少，擬合精度增長。

 

 

三. 由偏差導出設計

　　到此咱們已經創建了評估偏差的模型，但工做還遠未完成。偏差分析的最終目標是肯定設計，以前的工做只實現了從設計求解偏差；接下來的問題是如何已知偏差，求解設計。

　　事實上，對於g(x)的某一形式而言，若偏差模型的正向解法已經肯定，那麼經過枚舉搜索的方式就可反向求解。算法的輸入爲最大容許標準誤差S或最大容許峯值誤差Dmax；算法對於每一個給定的S或Dmax，經過n次迭代求出當前標準誤差或峯值誤差，當某次迭代得出的精度知足要求，那麼此時n即爲須要的階數；輸出n並結束。

　　利用上述算法繪製出「S的數量級--Taylor階數」圖。給定最大容許標準誤差，由圖可得最小展開階數。例如要求精度達到1e-3，至少應展開9階。

　　實現繪圖的Matlab代碼以下：

 1 n = 14;
 2 expS = logspace(-n,-1, n);
 3 
 4 ords = zeros(size(expS));
 5 for i=1:size(expS, 2)
 6     ords(i) = get_order(expS(i));
 7 end
 8 
 9 expSN = -n:-1;
10 plot(expSN,ords,'bo-.', 'LineWidth', 1);
11 
12 grid on;
13 xlabel('S <= (10^x)');
14 ylabel('Order of taylor');
15 set(gca,'xtick',[-n:1:-1])
16 
17 function v=taylor_n(x, n)
18     syms sx;
19     tayl = taylor(sin(sx), sx, 'Order', n);
20     v =eval(subs(tayl, sx, x));
21 end
22 
23 function v=get_order(expS)
24     n = 1;
25     S = 2;
26     while S > expS
27         t=0:pi/20:pi/2;
28         delta = sin(t)-taylor_n(t, n);
29         S = sqrt(sum(delta.^2) / size(t,2));
30         Dmax = max(abs(delta));
31         n = n + 1;
32     end
33     v = n;
34 end

 

　　對於任意經過若干次迭代來逼近目標函數的算法，上述評估框架仍然適用，只需簡單修改taylor_n函數的內容便可實現。

　　事實上，對於全部精度可變的算法，總存在一個可變因子n，當n增大時，算法精度增長，所以上述代碼可修改成通用的框架：

 1 n = 14; % 指定標準誤差S或峯值誤差Dmax的數量級範圍[10^-1, 10^-n]
 2 expS = logspace(-n,-1, n);
 3 
 4 ords = zeros(size(expS));
 5 for i=1:size(expS, 2)
 6     ords(i) = get_order(expS(i));
 7 end
 8 
 9 expSN = -n:-1;
10 plot(expSN,ords,'bo-.', 'LineWidth', 1);
11 
12 grid on;
13 xlabel('S <= (10^x)');
14 ylabel('Order of taylor');
15 set(gca,'xtick',[-n:1:-1])
16 
17 %%
18 function v=fit_n(x, n)
19     % 在此處添加逼近算法的實現
20 end
21 %%
22 
23 function v=get_order(expS)
24     n = 1;
25     S = 2; Dmax = 2;
26     while S > expS % 約束條件：可改成 Dmax > expDmax 來約束峯值誤差
27         t=0:pi/20:pi/2;
28         delta = sin(t)-fit_n(t, n);
29         S = sqrt(sum(delta.^2) / size(t,2));
30         Dmax = max(abs(delta));
31         n = n + 1;
32     end
33     v = n;
34 end

 

　　進一步地，若將多種不一樣的g(x)算法結合起來，綜合考慮各類算法的資源佔用、運行時間等因素，並編寫程序完成篩選（例如優先選擇精度高的算法，當精度相等時再優先選擇效率高的算法），則設計過程就能夠自動化。

　　相比於依賴設計者經驗來篩選算法，分別對候選算法進行偏差及資源評估，再經過對比作出最終決策的方法而言，本文所述的方法在必定程度上提升了工做效率。同時，因爲偏差評估模型不依賴於具體的算法實現，任意可行的算法均可以添加到框架中來，從而豐富框架的內容。