參考網址

 

1. wiki - Trigonometric_functions

2. wiki - 三角恆等式

3. glibc - sin_S.c 

4. netlib - k_sin.c

5. polygonal lab. - Fast and accurate sine / cosine approximation ( source )

6. csdn - 如何快速且高精度地求sin函數?

7. csdn - 哪位給個sin/cos的優化版?(C/C++)

8. csdn - 誰能提供一下sin(x) cos(x)的數值計算公式?


 

相關公式

 

1. sin(x) = x - x^3/3! + x^5/5! - x^7/7! + ...

2. cos(x) = 1-x^2/2!+x^4/4! - x^6/6!

3. sin(-x) = -sin(x) , cos(-x) = cos(x)

4. sin(PI - x) = sinx , cos(PI-x) = -cos(x)

5. sin(2x) = 2sinx cosx

6. cos(2x) = cosx^2 - sinx^2 = 2*cosx^2 - 1 = 1 - 2*sinx^2

6. sin(3x) = 3sinx - 4sinx^3 , cos(3x) = 4cosx^3 - 3cosx

7. sin(x)≅ x , cos(x)≅ 1-x,  where x≅ 0

8. sin(nx) = 2cosx * sin[(n-1)x] - sin[(n-2)x]


環境與定義

 

硬體環境:AMD Athlon(tm) II X2 245 Processor 2.91GHz / 2.75G RAM

作業系統:Microsoft Windows XP SP3

編譯器:Visual C++ 2010

編譯參數 (摘自 sin 專案):

/Zi /nologo /W3 /WX- /O2 /Oi /Oy- /GL /D "_MBCS" /Gm- /EHsc /GS /Gy /fp:precise /Zc:wchar_t /Zc:forScope /Fp"Release\sin.pch" /Fa"Release\" /Fo"Release\" /Fd"Release\vc100.pdb" /Gd /analyze- /errorReport:queue

效能比定義:若 sin(x) 執行 T1 秒,mysin(x) 執行 T2 秒,mysin(x) 效能為 sin(x) 之 (T1 * 100/ T2) %。即若 T1 = 5,T2 = 20,則 mysin 效能為 sin 之 500/20 = 25 %。

誤差測試:以 -2000*PI~ +2000*PI step PI*5E-4,分別在精度為 {1E-3,1E-6,1E-9,1E-12,DBL_EPSILON} 情況下,測試、保留與實際 sin(x) 之最大誤差絕對值。

效能測試:以 QueryPerformanceFrequency 高精度計時器計時 (準度約至1E-6 sec),測試 2000*PI~ +2000*PI step PI*5E-4 ,分別在精度為 {1E-3,1E-6,1E-9,1E-12,DBL_EPSILON} 情況下 所需時間,換算效能比。

 

學生作業版

 

double factor(int x)
{
    int i;
    double rst=1.0;
    for(i=2; i<=x; ++i)
        rst*=i;
    return rst;
}

double sinx(double x, int times)
{
    int i;
    double rst=0.0;
    double pwr;
    for(i=0; i<times; ++i){
        pwr = pow(x, 2*i+1);
        if(i%2==1) pwr=-pwr;
        rst = rst + pwr / factor(2*i+1);
    }
    return rst;
}

 

分析:完全不想分析。讓人很疑惑,為何授課老師會接受這種解法。速度慢確實不是探討重點,但準確性很差,另也會出現絕對值大於 1 之情況;且當 x 絕對值太大時,容易出現 I.N.D 。

優點:無腦。

缺點:速度慢、精準度不夠、可接受值域小。

 

基本款

 

一般在開放性論壇裡,所談論之 sin 都已先假定一定範圍再行效能、準確度評估。如先將 sin 調整到 [0, PI] 間,或先將 sin 調整到 [0, 2PI] 間,接下來才會 present 其 fast sin。在實際面上我不確定是不是有些情況必保證 sita 角都會在這麼理想的範圍,也因不確定,故本文所有正式之 function,均會先將 sita 調到 [-2PI, 2PI] 間進行,效能上會慢 7% 左右。

 

做三部份改善。

Step 1 : 先將輸入之 x 做 mod 2PI,使其範圍落在 [-2PI, 2PI] 或 [0, 2PI] 間。

Step 2 : 把 pow、 factor 拿掉,全都在 loop 裡執行。

Step 3 : 改為最小迭代誤差做收斂條件。

Step 1 要補充, 若 sita = x 為負,fmod 傳回值視 compiler 不同,可能為正,也可能為負,但算三角函式而言差異不大。

 

#define PI 3.1415926535897932
#define PI2 6.2831853071795864
double sin_basic(double x, double eps)
{
    double xn, nx2;
    double pwr=6.0, rst, rst2;
    double cnt=4.0;
    rst2=x=fmod(x, 2*PI);
    nx2=-x*x, xn=nx2*x;
    do{
        rst  = rst2;
        rst2 += xn/pwr;
        xn*=nx2;
        pwr*=cnt*(cnt+1.0);
        cnt+=2.0;
    }while(rst2-rst>eps || rst-rst2>eps);
    return rst2;
}

 

準確度測試

max error=7.210e-005(for eps = 1.000000e-003)
max error=4.678e-008(for eps = 1.000000e-006)
max error=3.416e-011(for eps = 1.000000e-009)
max error=2.517e-013(for eps = 1.000000e-012)
max error=2.652e-013(for eps = 2.220446e-016)

效能測試

eps=1.00e-003: 33.76 %
eps=1.00e-006: 28.14 %
eps=1.00e-009: 24.68 %
eps=1.00e-012: 22.24 %
eps=2.22e-016: 20.03 %

優點:

(1) 誤差式收斂
(2) 先行 fmod 2PI 增快收斂速度
(3) power / factor 放在 loop 裡,避開不必要開銷。

 

注意到雖 eps 最小給的是 1E-3,但實際測出來卻在 1E-4 左右,這個值若是在一般螢幕不大的座標轉換,應是夠用的,但在速度上差 math.h 還是差太多。既然 sin 都寫好了,那就順便寫 cos 吧。 

 

double cos_basic(double x, double eps)
{
    double nx2, xn, pwr=2.0;
    double rst, rst2=1.0;
    int    cnt=3;
    x = fmod(x, PI2);
    nx2 = -x*x, xn=nx2;
    do{
        rst = rst2;
        rst2 += xn/pwr;
        xn*=nx2;
        pwr*=cnt*(cnt+1);
        cnt+=2;
    }while(rst2-rst>eps || rst-rst2>eps);
    return rst2;
}

 

細節改善

 

由於細節可改善部份非常多,部份可議性也很大,慢慢敘述。

(1) 算階乘部份:原本公式是除以 x!,故直接改成 1.0/x!,再乘上去。

(2) 將除法改為乘法:在速度上,乘法比除法快,故盡可能先定義出新之常數,將除法以乘法取代。

(3) fmod 部份:整個三角函式寫得好不好,fmod 是重要關鍵之一。fmod 可議性我認為是較大的。一般在除法中,商數範圍落於 int 範圍中的話,可用公式去解決,不會有問題;若商數結果落於 int 範圍外則會出包,這部份視個人予以針著。這份程式碼將取代 fmod 部份。

(4) 調整至 [-PI, PI]:上述是將角度調到 [-2PI, +2PI],而 sin 又有個性質,sin(-x) = -sin(x),之所有縮小其範圍,原因為收斂速度較快。

 

事實上關於 fmod 部份,有另一種取代方案,

x=fmod(x, PI2);  (1)

改成

double inte = floor(x * INV_PI2);  (2)
x = x - inte * PI2;

文中說的是

int inte = (int)(x * INV_PI2);   (3)
x = x  - inte * PI2;

速度上, (1) < (2) < (3) ,且光是這段速度便差很多!只是 (3) 「並非完全安全」而已。

 

以下所有程式,均以此程式為基礎架構,若對溢位問題感到不安,又不想用 fmod ,則可以 floor 方式替代。

 

double sin_basic(double x, double eps)
{
    double xn, nx2;
    double pwr=0.1666666666666667, rst, rst2; // pwr = 1.0/6.0
    double cnt=4.0;    
    int inte  = (int)(x*PI_INV);

    rst2 = x-PI*inte;
    nx2=-rst2*rst2, xn=nx2*rst2;
    do{
        rst  = rst2;
        rst2+=xn*pwr;
        xn*=nx2;
        pwr/=(cnt*(cnt+1.0));
        cnt+=2.0;
    }while(rst2-rst>eps || rst-rst2>eps);
    return inte&1 ? -rst2 : rst2;
}

 

準確度測試

max error=4.260e-005(for eps = 1.000000e-003)
max error=2.359e-008(for eps = 1.000000e-006)
max error=1.680e-011(for eps = 1.000000e-009)
max error=1.788e-013(for eps = 1.000000e-012)
max error=1.745e-013(for eps = 2.220446e-016)

效能測試

eps=1.00e-003: 65.70 %
eps=1.00e-006: 51.64 %
eps=1.00e-009: 42.54 %
eps=1.00e-012: 35.39 %
eps=2.22e-016: 30.62 %

優缺點

(1) 誤差比剛剛再小一點。

(2) 使用不非常正確方式去掉 fmod,若數值範圍超出 int 表示範圍,結果為誤。

(3) 效能提昇明顯。

 

由於程式碼長,整份函式並不適合用 inline,同時用 __force_line 、__attribute((always_inline)) 效果並不會有明顯差異,

甚至噁心一點的程式碼

 

#define _sin_basic(y, x, eps){ \
    double xn, nx2; \
    double pwr=1.0/6.0, rst, rst2; \
    double cnt=4.0; \
    int inte =(int)(x*PI_INV); \
    rst2 = x-PI*(int)(inte); \
    nx2=-rst2*rst2, xn=nx2*rst2; \
    do{ \
        rst = rst2; \
        rst2+=xn*pwr; \
        xn*=nx2; \
        pwr/=(cnt*(cnt+1.0)); \
        cnt+=2.0; \
    }while(rst2-rst>eps || rst-rst2>eps); \
    y2 = inte&1 ?-rst2 : rst2; \
}

 

這也不會有明顯的改善。

 

 

 倍角公式

 

承如上所言,sin 收斂有個性質,當 x 愈大時,收斂愈慢,故若可把角度先縮得愈小便愈有利。

所以我們希望能夠找出可以套用 sin(nx) 或  sin(x/n)的公式,

到時只要縮小 n 倍,便可加快收斂速度,減少迭代次數。

 

幸運的是,已有人給了一份 sin(nx) 公式  sin(nx) = 2cosx * sin[(n-1)x] - sin[(n-2)x]

但不運的是,這份公式若實際撰之,遭遇了二個問題:額外呼叫 cosx 、遞回關係,速度可能較慢。

這公式筆者認為還不太算實用,再另尋出路。

 

一般而言,不少人都說「二倍角很好用」,但對 sin 我覺得很不好用。看一下 sin  之二倍角公式

sin(2x) = 2sinx cosx

重點是卡在 cos,要算 sin 時還要再呼叫 cos,而在呼叫其他 function 時認為效率或多或少會拖慢,

( 整個 cos 應是不可能讓 compiler 用 inline,除非有類似 __force_inline 這麼暴力的方法,但不建議),

但幸運的是,我們也知道 3 倍角公式

sin(3x) = 3sinx - 4sinx^3  = sinx * (3 - 4 sinx)

依此為準則,程式碼約如下述。

 

double sin3(double x, double eps)
{
    double xn, nx2;
    double pwr=0.1666666666666667, rst, rst2;
    double cnt=4.0;
    int inte  = (int)(x*PI_INV);
    rst2 = (x-inte*PI)*0.3333333333333333;
    nx2=-rst2*rst2, xn=nx2*rst2;
    do{
        rst  = rst2;
        rst2 += xn*pwr;
        xn*=nx2;
        pwr/=(cnt*(cnt+1.0));
        cnt+=2.0;
    }while(rst2-rst>eps || rst-rst2>eps);
    return inte&1 ? -rst2*(3.0 - 4*rst2*rst2) : \
           rst2*(3.0 - 4*rst2*rst2);
}

 

準確度測試

max error=2.479e-005(for eps = 1.000000e-003)
max error=3.095e-008(for eps = 1.000000e-006)
max error=9.288e-012(for eps = 1.000000e-009)
max error=1.738e-013(for eps = 1.000000e-012)
max error=1.738e-013(for eps = 2.220446e-016)

效能測試

eps=1.00e-003: 86.65 %
eps=1.00e-006: 69.03 %
eps=1.00e-009: 58.77 %
eps=1.00e-012: 51.29 %
eps=2.22e-016: 44.40 %

 

上面只用到三倍角,接著還有其他想法:把 n=10, 20, ..., 50 展開來如何?這想法是不錯,

不過我認為展開除了需要分析能力外,一些誤差也會開始變大,但有個看起來還算不錯的公式可以試試

令 sin(x) = y

sin(5x) = 16*y^5 - 20*y^3 + 5*y = y*(5+y^2*(-20+16*y^2))

這次再事先展開出第二項,程式碼如下

 

double sin5(double x, double eps)
{
    double xn, nx2, y2;
    double pwr=8.3333333333333333E-3, rst, rst2;
    double cnt=6.0;
    int inte  = (int)(x*PI_INV);
    x = (x-inte*PI)*0.2;
    nx2=-x*x, xn=nx2*nx2*x;
    rst2 = x*(1.0+nx2*0.1666666666666667);
    do{
        rst  = rst2;
        rst2 += xn*pwr;
        xn*=nx2;
        pwr/=(cnt*(cnt+1.0));
        cnt+=2.0;
    }while(rst2-rst>eps || rst-rst2>eps);
    y2 = rst2*rst2;
    return inte&1 ? -rst2*(5.0+y2*(-20.0+16.0*y2)) : \
           rst2*(5.0+y2*(-20.0+16.0*y2));
}

 

準確度測試

max error=4.715e-005(for eps = 1.000000e-003)
max error=1.203e-008(for eps = 1.000000e-006)
max error=4.272e-012(for eps = 1.000000e-009)
max error=3.514e-013(for eps = 1.000000e-012)
max error=3.467e-013(for eps = 2.220446e-016)
y1 = -0.000059, ts=0.486501

效能測試

y1=-0.000106, eps=1.00e-003: 105.96 %
y1=-0.000059, eps=1.00e-006: 79.88 %
y1=-0.000059, eps=1.00e-009: 68.71 %
y1=-0.000059, eps=1.00e-012: 60.71 %
y1=-0.000059, eps=2.22e-016: 52.50 %

 

其實上面的在 1E-3 時效能本來是 85%,但已事先多展開一項,使得效能往上拉,多展開一項就變比較好,這結果未必是一定的,也可能是開了 O2 關係。sin 寫完, cos 呢?關於 cos 有幾個問題要提出來

(1) 關於 6 倍角:事實上 cos 6 倍角公式很漂亮 ,令 cos(x) = y,cos(6x) = 16y^6-16y^4+2y^2+1,但六倍角公式誤差實在是太大了,即使用 math.h 裡之 cos 換算驗證,誤差仍相當大,故目前最大只做到5 倍角。

(2) 對稱性問題:上述之 sin 最後能順利轉換,全鑑於一特性:一、二象限為正,三、四象限為負,故在除 PI 後,只要確定商是正或負便可決定其正負號。但 cos 並非如此,它為一、四象限為正,二、三象限為負。此一特性要得知正負號,必須進行 mod PI/4,且後續步驟較鎖碎繁複,故 cos 並不適合,調整角度至 [-PI, PI],只能調整到 [-2PI,2PI],這不論是速度或誤差面都顯得較差。

鑑於上面兩點敘述和剛剛 sin 多展開一次的結果,這次對 cos 再多展開一次,而調整限於 [-2PI,2PI] 裡。

 

double cos5(double x, double eps)
{
    double nx2, xn, pwr=1.3888888888888889E-3;
    double rst, rst2, y2;
    double cnt=7.0;
    const int inte=(int)(x*PI2_INV);

    x = (x-inte * PI2) * 0.2;
    nx2 = -x*x, xn=nx2*nx2*nx2;
    rst2 = 1.0 + 0.5*nx2*( 1.0 + nx2*0.08333333333333333);
    do{
        rst = rst2;
        rst2 += xn*pwr;
        xn*=nx2;
        pwr/=(cnt*(cnt+1.0));
        cnt+=2;
    }while(rst2-rst>eps || rst-rst2>eps);
    y2=rst2*rst2;
    return rst2*(5.0+y2*(-20.0+16.0*y2));
}

 

準確度測試

max error=9.762e-005(for eps = 1.000000e-003)
max error=3.007e-008(for eps = 1.000000e-006)
max error=3.009e-011(for eps = 1.000000e-009)
max error=3.543e-013(for eps = 1.000000e-012)
max error=3.485e-013(for eps = 2.220446e-016)

效能測試


eps=1.00e-003: 91.84 %
eps=1.00e-006: 76.26 %
eps=1.00e-009: 62.52 %
eps=1.00e-012: 53.50 %
eps=2.22e-016: 45.41 %

由於只調整至 [-2PI, 2PI],改善沒像 sine 來得大,cosine 多展開一項,從原本 82 % 拉到 91 % 。至於再多展開一項結果會不會更好?我想還是有限的。

對 sin  而言,有更快的方法嗎?若在可接受較大誤差 (1E-3~1E-5) 時,我想是有的,最簡單的方式便是將 global search 找到的最大誤差發生之 x value 紀錄下來,同時再以程式看進行了幾次迭代,接下來再把 do-while 去掉,以 Horner Rule 方式展開。但這部份效果也未必比原本的還好,因假設一般平均收斂次數為 3 次,但最差情況為 7 次收斂,每次都固定算 7 次,平均就多了 4 次的計算時間。這部份還是必須以實際資料分佈範圍與視實際情形而定,來得較佳。

 

切割至第一象限 ?

 

目前,筆者將角度全都換到 [-PI, PI] ,乃是根據二性質

 

sin ( 2kpi + x) = sin (x)

sin ( -x) = - sin(x)

 

而當x 趨近於 0 時,頂多發生 +-0 問題,而有種作法是要切到第一象限去,

這點是我納悶、疑問的地方之一。目前看過某些文章,將 sita 角度全都轉換到第一象限去,

也就是將 sita 角取 PI/2 之模數,這樣反倒出現幾個問題:

若 sita 角剛好為 0.5PI 或 1.5 PI 時將發生以下現象

sin(0.5 PI) = 1  , sin(0.5PI % 0.5PI) = sin(0) = 0

得到的結果誤差為 1,要處理這種現象反而不好處理

處理不當, branch 使用過多,反而因此拖垮整體速度,故筆者較不建議直接轉到第一象限去,

除非非常肯定,該三角函數值必為正值再這麼做,此時效果將更佳。

 

建表法趨近

 

這問題之前我曾實作、回答過別人,後來發現有 pdf 放在網路上提供下載。

原問題是將 0~360 度 (使用的是度度量,非弳度量) 切成 1000 份後,

要用查表方式估算 sin 值。

 

事後想想,在沒有提供三角函式 library 環境、不要求精度情況下,

事先先用建表法也不失一對策,至於要用度度量還是弳度量,視當時需求或 coder 習慣為佳。

 

建表法一開始便必須先調用其他精度高的 sin function,初始化建好表後,

日後要查大多便直接以線性內插法方式做查表。以 0~2PI 方式做查表,

欲建立 SIZE = 1024 之表,程式碼大既如下

 

#define SIZE 1024
static double sin_tablePI2[SIZE+1];
static double unit;
void sin_init()
{
    int i;
    unit = PI2/SIZE;
    for(i=1; i<=SIZE; ++i)
        sin_tablePI2[i] = sin(i*PI2/(SIZE+1.0));
}
double sinL(double x)
{
    unsigned i;

    x = fmod(x, PI2);
    if(x<0.0) x+=PI2; // 確保為正 
    i = x/unit;       // 算區間

    // 一次插值
    return sin_tablePI2[i]+ \
        (x-i*unit)*(sin_tablePI2[i+1]-sin_tablePI2[i])/unit;
    // return sin_tablePI2[i]; // 不做內插
}

 

初始化是程式使用前呼叫,以長期調用而言,所佔時間幾乎是省略沒影響。之後便一勞永逸,都用查表方式,上述在修正值之方式為使用一階線性內插法,若去掉效率並沒明顯改善,反倒是精度變更差。效能約只到 50% ,精度約 6.5E-3。較好的做法是,傳回前先判斷 x 值與算出來之 i 值誤差多少,若夠小 (如 1E-3) 則直接傳回 sin_tablePI2[i],反之才進行內差;同樣的概念,也該先檢查 x 與 i+1 之誤差。

但看得出來,這份建表法並沒把效率全都做出來,在前半段裡,可以直接仿 sin5 手法,接下來才做查表動作。如此一來,一樣大小的表能提供之精準度提高將近 10 倍,但速度上卻不會再增加 (頂多就增加 fmod 那段),有興趣者可自行測試研究。

edisonx 發表在 痞客邦 PIXNET 留言(0) 人氣()