算是回答一個網友的問題。(真希望回答問題有錢可拿.....)

在沒有標準程式庫的情況下,計算 sin、cos 這些函數值最直覺的方法可能是用級數逼近。但假如我們只想求有限範圍內離散數字的三角函數值,級數逼近就顯得有點複雜了。

舉例來說,假如你在設計火砲射控系統,那麼處理的單位應該會是密位(mil),北約制的密位定義為1圓周 = 6400密位

最簡單的想法是建立一個 6400 個元素的陣列,直接用密位查表,這裡就不示範。但是這個作法相當消耗記憶體,有沒有更節省記憶體的作法呢?

好,現在回溯到國中時期的記憶。在中學時期我們背的三角函數值用五根手指就數得出來,但考試還是會考一些奇怪的角度,還記得是怎麼解題的嗎?這些題目之所以可以解,是因為我們還硬背了很多三角函數公式。假如你有一套0~360間隔20度的表,和一套0~20度間隔1度的表,總共才多少筆資料?想一下用合角公式你可以回答多少角度的 sin、cos 的問題?

這就是查表法精簡的關鍵,運氣好的時候我們可以將問題拆解成若干個小問題各別查表,最後再合成回答案。

回到密位的問題,我們將 6400 分解成 80x80,用兩層 80 個元素的陣列表示。取80的理由很簡單,因為越接近平方根就能用越少筆資料合成同樣多組答案。建表程式如下,我沒有將 sin、cos 分成兩表,而是變成陣列的第二維度,是因為用這個方法總是同時取兩者,放在一起可以增進 locality。

const double PI = 3.1415926535897932;

double TRIG_TABLE1[80][2];
double TRIG_TABLE2[80][2];


void InitTrigTables()
{
    const double t1Unit = (2 * PI) / 80.0;
    const double t2Unit = (2 * PI) / 6400.0;

    for (int i = 0; i < 80; ++i) {
        double r = double(i) * t1Unit;
        TRIG_TABLE1[i][0] = sin(r);
        TRIG_TABLE1[i][1] = cos(r);
    }

    for (int i = 0; i < 80; ++i) {
        double r = double(i) * t2Unit;
        TRIG_TABLE2[i][0] = sin(r);
        TRIG_TABLE2[i][1] = cos(r);
    }
}

密位的 sin 函數:

double SinMil(int mil)
{
    assert(0 <= mil && mil < 6400);

    // sin(a+b) = sin(a)*cos(b) + cos(a)*sin(b)
    double* a = TRIG_TABLE1[mil / 80];
    double* b = TRIG_TABLE2[mil % 80];
    return a[0] * b[1] + a[1] * b[0];
}

cos作法也相同,實用上還會再加一個同時取sin和cos的函數,這裡就不列舉。比起直接建 6400 大小的表,這個方法當然是多了一些計算成本,大約多了兩個浮點乘法、一個加法、還多查一次表。但是省下的空間很可觀,在這裡用了 320 個 double 就能表示 6400 密位的 sin、cos 表,用原始的方法需要 6400 x 2 = 12800 個 double 才能達到同樣的效果(好啦,利用sin和cos的關係應該可以少一點)。

以下使用密位轉換弧度函數來驗證正確性;另外提供弧度轉密位函數,如此一來也可以用上面的程式求取連續值的三角函數值。請注意以下只是示意,沒有特別處理正負範圍和 rounding,還有其中的常數值可以事先算起來加速:

inline double Mil2Rad(int m)
{
    return m * (2.0 * PI) / 6400.0;
}

inline int Rad2Mil(double r)
{
    return (r * 6400.0) / (2.0 * PI);
}
...

InitTrigTables();

cout << SinMil(3251) - sin(Mil2Rad(3251)) << endl;  //--test1
cout << SinMil(1830) - sin(Mil2Rad(1830)) << endl;

cout << SinMil(Rad2Mil(2.5)) - sin(2.5) << endl;    //--test2
cout << SinMil(Rad2Mil(0.8)) - sin(0.8) << endl;

和預期的差不多,test1的誤差低於小數以下15位,已經和內建函數沒太大差別。而test2誤差大約是小數下3~4位,畢竟密位解析度只有 6400 而已。不過想像一下,只要有兩層 1024 的表(別忘了1024可用位元運算代替除餘),解析度將可以達到 20 bit,已經很接近 float 的精度了。對於精度需求不超過七位的應用,這個方法搭配比較「非正規」的 rounding 手法,速度有機會超過內建的三角函數甚至 x86 的 fsin、fcos 指令。


註:密位的原意是千分之一弧度,一圈大約可劃分成 6283.185 密位,但這個數字實在是難記又難算,所以實用上的密位都直接規定為整數值,例如北約取 6400,前蘇聯和對岸用 6000,瑞典之前用過6300。

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