[C語言數值分析] cmath / math.h - sin.cos @ Edison.X. Blog

參考網址

相關公式

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 那段)，有興趣者可自行測試研究。