novus log

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

在沒有標準程式庫的情況下，計算 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。