一種簡單，快速，精准的sin/cos函數模擬，及as3實現

最新推荐文章于 2022-02-25 21:57:21 发布

Iam_Rocky

最新推荐文章于 2022-02-25 21:57:21 发布

阅读量2.6k

点赞数

分类专栏： [Math] 文章标签： c 制造

[Math] 专栏收录该内容

9 篇文章 0 订阅

订阅专栏

http://www.devmaster.net/forums/showthread.php?t=5784

http://lab.polygonal.de/2007/07/18/fast-and-accurate-sinecosine-approximation/

在某些情況下我們需要一些更高效的且近似於標准值的 sin 和cos函數。

有時候我們並需要過高的精度，這時 C語言中自帶的三角函數（sinf() 和 cosf() f）計算的精度超出了我們所需要的精度要求，所以其效率很低。我們真正需要的是間於精度和效率的一個折中的方案。眾所周知的取近似值的方法是：泰勒級數（和著名的馬克勞林級數）

代碼是：

x - 1/6 x^3 + 1/120 x^5 - 1/5040 x^7 + ...

我們繪制了下圖:

綠線是標准的sin函數，紅線是4項泰勒級數展開式。這個近似值的效果看起來還不錯，但是如果你仔細觀察後會發現

它在pi/2之前的效果還是很好的，但是超過了pi/2後就開始快速偏離標准sin。它在pi處的值比原來的0多了0.075.用這個方法來模擬波動忽動忽停，看起來很呆板，這個效果當然不是我們想要的。

我們可以繼續增加泰勒級數項的個數來減小誤差，但是這將導致我們的公式非常的冗長。用4項的泰勒級數展開式需要我們進行7次乘法和3次加法來完成。泰勒級數不能同時滿足我對精度和效率的要求。

剛剛近似如果能滿足sin（pi）=0就好了。從上圖我還可以發現另一件事：這個曲線看起來很象拋物線！所以我們來尋找一個盡可能和sin接近的拋物線（公式）。拋物線的范式方程是：A + B x + C x^2.這個公式可以讓們控制三個自由度。顯然我們需要其滿足sine(0) = 0, sine(pi/2) = 1 and sine(pi) = 0. 這樣我們就得到了3個等式。A + B 0 + C 0^2 = 0

A + B pi/2 + C (pi/2)^2 = 1

A + B pi + C pi^2 = 0

解得：A = 0, B = 4/pi, C = -4/pi^2.我們的拋物線誕生啦！

貌似這個的誤差看起來比泰勒級數還要遭。其實不是的！這種方法的最大誤差是0.056.（譯者：而且這種近似值沒有誤差積累）而且這個近似值的繪制出的波動是光滑的，而且只需要3次乘法和一次加法。不過它還不夠完美。下圖是[-pi, pi] 之間的圖像：

顯然我們至少需要它在1個完整的周期內都符合我們要求。但是我們可以看出，我們缺少的另一半是原拋物線的一個映射。它的公式是：4/pi x + 4/pi^2 x^2。所以我們可以直接這樣寫：

Code:

if(x > 0) { y = 4/pi x - 4/pi^2 x^2; } else { y = 4/pi x + 4/pi^2 x^2; }

添加一個條件分支不是一個好的方法。它會讓程序漸漸的變慢。但是觀察一下我們模擬的和標准的圖像是多麼的接近啊！觀察上面兩式子，只是中間的一項正負號不同，我的第一個單純的想法是可以提取x的正負號來消除分支，即使用：x / abs(x)。除法的代價是非常大的，我們來觀察一下現在的公式： 4/pi x - x / abs(x) 4/pi^2 x^2。將除法化簡後我們得到：4/pi x - 4/pi^2 x abs(x).所以只需要多一步運算我們就得到了我們需要的sin逼近值！下圖是結果

接下來我們要考慮cos。有基礎的三角公理可以知道：cos(x) = sin(pi/2 + x).把x多加一個pi/2就可以搞定了？事實上它的某一部分不是我們期望得到的。

我們需要做的就是當x > pi/2時「跳回」。這個可以由減去2 pi來實現。

Code:

x += pi/2;

if(x > pi) // Original x > pi/2 { x -= 2 * pi; // Wrap: cos(x) = cos(x - 2 pi)}

y = sine(x);

又出現了一個分支，我們可以用邏輯「與」來消除它，像是這樣：

x -= (x > pi) & (2 * pi);

Code:

x -= (x > pi) & (2 * pi);

注意這並不是c的源代碼。但是這個應該可以說明它是怎麼樣運行的。當x > pi是false 時，邏輯「與」（&）運算後得到的是0，也就是（x-=0）大小沒有改變，哈哈完美的等價!我會給讀這篇文章的讀者留一些關於這個練習。雖然cos 比sin需要多一些運算，但是相比之下貌似也沒有更好方法可以讓程序更快了。現在我們的最大誤差是0.056 ，四項泰勒級數展開式每一次都會有一點點誤差。再來看看我們sin函數：

現在是不是不能繼續提升精准度了呢？當前的版本已經可以滿足大多度sin函數的應用了。但是對一些要求更高一些的程序現在做的還夠。

仔細觀察圖像，你會注意到我們的近似值總是比真實值大，當然除了0，pi/2 和 pi。所以我們要做的就是在不改變這些點（0，pi/2 和 pi）的情況下，將函數再「按下去」一些。解決方法是利用拋物線的平方。看起來就像這樣：

注意它保持著原來那些關鍵點，不同的是它比真實的sin函數值更低了。所以我們可以用一個加權的平均值來使兩個函數更接近。

Code:

Q (4/pi x - 4/pi^2 x^2) + P (4/pi x - 4/pi^2 x^2)^2

利用Q + P = 1. 你可以靈活的控制絕對誤差或相對誤差。別急我來告訴你取不同的極限結果時Q，P的值。絕對誤差的最佳權值是：Q = 0.775, P = 0.225 ；相對誤差的最佳權值是：Q = 0.782，P = 0.218 。讓我們來看一下結果的圖像。

紅線呢？它幾乎被綠線完全覆蓋了，這足以證明我們的近似十分完美。最大誤差是0.001，50倍的提升！這個公式看起來很長，但是括號裡面的公式最終得到的值是相同的，也就是說括號裡的只需要被計算一次。事實上在原來的基礎上只是增加了額外的2次乘法和2次加法就可以得到現在的結果。

先別高興的太早，我們還要「制造」一個負號出來。我們需要增加一個abs()運算。最終的c代碼是：

Code:

float sine(float x)

{

const float B = 4/pi;

const float C = -4/(pi*pi);

float y = B * x + C * x * abs(x);

#ifdef EXTRA_PRECISION // const float Q = 0.775;

const float P = 0.225;

y = P * (y * abs(y) - y) + y;

// Q * y + P * y * abs(y)

#endif }
所以我們僅僅是需要多加5次乘法和3次加法就可以完成了。如果我們忽略abs()這個仍然是比4項泰勒級數展開式快，更精准！Cos只需要相應的變換一下x就可以了。

(譯者注：後面是匯編程序，不翻譯了)

part2

我選取了最小誤差的情況，用as3運行後發現提升了14倍，而且仍然是非常精准。不過你必須直接使用它，不能把它放到一個函數中，因為每調用一次額外的函數調用會削減執行效率，最終你會得到一個比Math.sin() 和 Math.cos()效率更差的結果。還有這裡會用到的三角定理：

cos(x) = sin(x + pi/2)

cos(x - pi/2) = sin(x)

下載: fastTrig.as.

可以清楚到對比結果，現在你可以用這個替換Math.sin() 和 Math.cos()了

哇哦！！！幾乎是相同的精准度（14倍速度提升）

//always wrap input angle to -PI..PI
if (x< -3.14159265)
    x+= 6.28318531;
else
if (x>  3.14159265)
    x-= 6.28318531;

//compute sine
if (x< 0)
    sin= 1.27323954 * x+ .405284735 * x* x;
else
    sin= 1.27323954 * x- 0.405284735 * x* x;

//compute cosine: sin(x + PI/2) = cos(x)
x+= 1.57079632;
if (x>  3.14159265)
    x-= 6.28318531;

if (x< 0)
    cos= 1.27323954 * x+ 0.405284735 * x* x
else
    cos= 1.27323954 * x- 0.405284735 * x* x;
}

High precision sine/cosine (~8x faster)

//always wrap input angle to -PI..PI
if ( x < - 3 . 14159265 )
    x += 6 . 28318531 ;
else
if ( x > 3 . 14159265 )
    x -= 6 . 28318531 ;

//compute sine
if ( x < 0 )
{
    sin = 1 . 27323954 * x + . 405284735 * x * x ;

    if ( sin < 0 )
        sin = . 225 * ( sin *- sin - sin ) + sin ;
    else
        sin = . 225 * ( sin * sin - sin ) + sin ;
}
else
{
    sin = 1 . 27323954 * x - 0 . 405284735 * x * x ;

    if ( sin < 0 )
        sin = . 225 * ( sin *- sin - sin ) + sin ;
    else
        sin = . 225 * ( sin * sin - sin ) + sin ;
}

//compute cosine: sin(x + PI/2) = cos(x)
x += 1 . 57079632 ;
if ( x > 3 . 14159265 )
    x -= 6 . 28318531 ;

if ( x < 0 )
{
    cos = 1 . 27323954 * x + 0 . 405284735 * x * x ;

    if ( cos < 0 )
        cos = . 225 * ( cos *- cos - cos ) + cos ;
    else
        cos = . 225 * ( cos * cos - cos ) + cos ;
}
else
{
    cos = 1 . 27323954 * x - 0 . 405284735 * x * x ;

    if ( cos < 0 )
        cos = . 225 * ( cos *- cos - cos ) + cos ;
    else
        cos = . 225 * ( cos * cos - cos ) + cos ;
}