三角関数のテーブル引きによる高速化

[2002/08/16]

十分注意して作成していますが，ご自身の責任の下で使用して下さい！

はじめに

三角関数はマクローリン展開などで実装されているため低速です．これを角度をキーとしたテーブル引きで高速化しようというお話です．

作成方針

マニピュレータのシミュレーションなどにも使いたいため粗い近似ではない精度が必要となります．この場合にテーブルの分割数を細かく取らなければならなくなります．細かくするほど精度は向上しますが，確保できるメモリ容量との兼ね合いが出てきます．このためテーブル作成の範囲はなるべく制限して [０〜π/２]とし，適宜この範囲へ算出したい角度を変換することにします．また，作成するテーブルは cos() に関するものだけとして sin() も変換して値を引くようにします．

テーブル作成

次のようにして指定された分割数（num）で９０度分のテーブルを作成します． malloc() で動的に領域を確保し，Cos[i]←cos(i*π/2*num) のテーブルを作ります．（掲げたソースではエラー処理は省いています）．

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "sincos_table.h"

int SinCos_num=-1;
double *SinCos_table=NULL;

/* SinCos_init: ９０度分の cos() テーブルを作成する */
void SinCos_init(int num)
{
    int i;

    SinCos_table=malloc(sizeof(double)*num);
    SinCos_num=num;
    for (i=0; i<num; i++) SinCos_table[i]=cos(i*M_PI/2/num);
}

テーブル引き関数 Cos()

作成したテーブルから指定角度に応じて値を引く関数 Cos() を作成します．ここでの問題はテーブルが [0〜π/2] の範囲でしか作成されていないのに，指定角度は負の値も複数回の回転分の値も取り得ることです．これに対する処理としては次のようになります．

θ＜０のとき，cos(θ)=cos(-θ) よりθ←-θとして正の値にする．
θ←θ/π として以下θ＝１で半回転πを表す．
θ←{θ-int(θ/2)*2} で複数回転を除去．ここで int() は整数化の関数．
θ＞１のとき，θ←{２-θ} で [０〜１]の範囲に収める．
θ＞０．５のとき，θ←{１-θ} で [０〜０．５]の範囲に収める．この場合はテーブルの値に負符号を付けて返す．

これを実装したのが次のリストです．

/* Cos() 関数 */
double Cos(double x)
{
    double sign=1;
    int i;

    if (x<0) x=-x;              /*  cos(x)=cos(-x)     */
    x/=M_PI;
    x=(x-(int)(x/2)*2);         /*  cos(x)=cos(x+2*PI) */
    if (x>1) x=2-x;             /*  cos(x)=cos(x-PI)   */
    if (x>.5) {x=1-x; sign=-1;} /* -cos(x)=cos(PI-x)   */
    i=(int)(x*2*SinCos_num);
    return (sign*SinCos_table[i]);
}

テーブル引きのSin()は sin(x)=cos(x-π/2) で求めることにします．

/* Sin() 関数: sin(x)=cos(x-M_PI/2) */
#define Sin(x) (Cos(x-M_PI/2))

以上です．コンパイルにあたってはコンパイラの最適化オプションを付加しないと高速化につながらない場合があるようです．gcc の場合は -O3 以上が必要のようでした．次のようにしてオブジェクトファイル sincos_table.o を作成しリンクして使用します．リンク時には -lm オプションを付けるのを忘れずに．

gcc -O3 -c sincos_table.c

算出精度・時間について

ざっと見てみたところテーブルの分割数と精度との間には分割数 1e+X ⇔ 精度 1e-(X+1) の関係があるようです．つまり 1000 分割すると 0.0001 程度の精度が得られるということです．分割数については必要とする精度とメモリ容量との兼ね合いで調整していくしかないと思います．

また眼目の処理速度については１５億回の関数呼び出しについて手元のコンピュータでは

関数	時間(s)
cos()	30
Cos()	1

という結果となりました．当然，テーブル化の効果は呼び出し回数が増えるほど有効になります．ただし，これは前述の最適化オプションを付けた場合で，単に関数だけを呼び出しているので，コンパイラによって予想外の最適化（実際は関数を呼び出してないとか！？）をされている可能性もあります．テストプログラムを呼び出す場合は次のようにして -DSINCOS_TABLE_TEST を付けてコンパイルして下さい．

gcc -O3 -DSINCOS_TABLE_TEST -o sincos_table sincos_table.c -lm