三角関数の高速化

〜テーブルを使って〜

浮動小数点演算高速化の必要性

古のコンピュータ、とくにＰＣでは、浮動小数点計算は非常に時間の必要な処理でした。
もちろん、今のＰＣには、浮動小数点計算用のハードウェアが搭載されているのが普通で、浮動小数の計算もずいぶんと速くなりました。

けれども、浮動小数点演算ハードウェアが搭載されていないコンピュータもまだまだあるのです。 私にとって身近な例が、『ＰｏｃｋｅｔＰＣ』です。 このコンピュータのＣＰＵは、ＡＲＭという規格を基準にしてＩｎｔｅｌが作ったＸＳｃａｌｅというものですが、 これには、浮動小数点はおろか、整数の割り算命令もないようなのです。

モバイル用ＣＰＵで、省電力という大儀があるので省かれていると思うのですが、 私のやりたい天文計算では、浮動小数点が遅いのはかなり致命的です。 ＸＳｃａｌｅの動作クロックは４００ＭＨｚくらい、いくら昔のデスクトップＰＣよりもずっと早いとはいえ、 プログラムは早いに越したことはありません。

そこで、天文計算の中で特に使用頻度の高い三角関数の高速化について、調べたことを記録しておきます。 使用言語はＣ＋＋、最近はお手軽Ｗｅｂ関係言語が流行で、すっかり時代遅れ言語ですが、私がこれしか喋れないので。

浮動小数点フォーマット

これは基本中の基本ですが、そんなに精度がいらない場合は、doubleではなくfloatを使います。 というのも、floatは４ｂｙｔｅ、３２ｂｉｔで表せて、データ転送にdoubleの半分の時間でいいということと、 いわゆるintと同じ大きさなのでちょっとしたトリックが使えるためです。

詳細に入る前に、浮動小数がどのようにコンピュータの中で実現されているか確認です。
上の図は、IEEEの32bit浮動小数点フォーマットです。1bitの符号、8bitの指数、23bitの仮数で構成されています。 ＣやＣ＋＋でいうfloat型は、この浮動小数点フォーマットを使っています。
このフォーマットは、どのコンピュータ、ＯＳでも同じになっているはずです。そうでないと規格としての意味がありませんね。

表現範囲などはあちこちに資料がありますので、詳しくは述べませんが、２の±１２７乗の範囲で小数を表すことができます。 仮数部は正規化されていて、22bitの前に暗黙の１があるのも有名ですね。

このフォーマットで表される数ｎは、ｎ＝ｓ×１.ｍ×２＾（ｅ−１２７）となります。（＾は、べき乗を表す）

float → int トリック

ここで、浮動小数点の整数部分だけを取りたいと考えます。というのは、三角関数、特にsin、cosを、テーブルを用いて求めたいからです。

三角関数の高速化にはいくつかの方法があり、精度、速度にそれぞれ利点、欠点があります。 一番単純で高速なのが、０〜３６０度までのsin、cosの値をあらかじめ求めてテーブルにしておき、 必要に応じてそれを参照する方法です。 プログラムとしては、sin、cosの値を配列に持っておき、必要な角度をその配列の範囲にして参照することになります。 配列範囲としてよく使われるのは、２５６や１０２４など、コンピュータにとって都合のいい数字ですね。

この方法では、与えられた角度を配列の参照値とするために、値を正規化し、小数を整数にする必要があります。 ３０度の場合、３０×（２５６÷３６０）＝２１．３３３３３・・という値を求め、この整数部分、２１を使って配列を参照するわけです。

３２ｂｉｔ浮動小数フォーマットは、float型として使われるわけですが、ｂｉｔ並びとして考えると、整数としても見ることができます。 ここに、高速なfloat→int変換が可能となるトリックがあるのです｡

戯言：こういう、あるデータを別の形で見る、というのは、プログラムリテラシーのまさにリテラシーだと思うのですが、 あまりプログラミングの本には出てきませんね。bit列眺めてニヤニヤしているのはやっぱり変ですか？（笑）

ここで使うトリックが、

1. 与えられた角度の値に、浮動小数で１．５×２＾２３を足す
2. 求まった値から、整数で０ｘ４ｂ４０００００を引く

というものです。

なぜこれで整数部分が求まるのかというと、

1. １．５×２＾２３を足すことによって

   元の数の小数部分がアンダーフローして整数部分しか残らない

   符号部分と指数部分が０ｘ４ｂ４と決まった値になる

2. 整数で０ｘ４ｂ４０００００を引くことによって

   浮動小数フォーマットの指数部分が０になる

ためなのです。

マジックナンバー０ｘ４ｂ４０００００は、１．５×２＾２３を、浮動小数点フォーマットで表したｂｉｔ並びです｡ これを整数として引き算するので、指数部分が０となるわけです｡

　たとえば、このトリックで４１．２６の整数部分を取り出してみると、

という操作となり、０ｘ２９＝４１と、正しい値が求まります｡

気をつけないといけないのは、このトリックでは±２＾２２の範囲にある浮動小数の整数部分しか、 得ることができない、ということです。それよりも絶対値が大きい場合は、このトリックは使えません。桁があふれちゃうからですね。

sinテーブル参照版

　では、sinをテーブルを使って求めるプログラムを書いてみます。

// sinテーブルの大きさ
const   int     SINTABLESIZE    = 256;
// バイアス
const float     FTOIBIAS        = 12582912.0f;
// ０〜２π正規化用
const float     PI2_F           = 6.2831853f;   // 3.14 * 2;
const float     TWOPISCALE      = (float)SINTABLESIZE / PI2_F;

// sine table.
float   SinTable[SINTABLESIZE];

// float <=> int 共用体
typedef union
{
    int     i;          // as integer
    float   f;          // as float
} INTORFLOAT;

使用する定数を定義して、float＜＝＞intコンバータ共用体を定義します。

// sine table.初期化
void    s24fSinCosInit( void )
{
    for ( int i = 0; i < SINTABLESIZE; ++i) {
        double  theta = (double)i/TWOPISCALE;
        SinTable[i]   = (float)sin(theta);
    }
}

fsinを呼び出す前に、あらかじめテーブルを初期化します。ランタイムライブラリと互換をとりたいので、与える値は０〜３６０ではなく、 ０〜２πとします｡

float   fsin( float theta )
{
    INTORFLOAT      fi;

    // 与えられた値を０〜２πに正規化して１．５＾２３を足す
    fi.f    = theta * TWOPISCALE + bias.f;
    // 整数化＆２５６の範囲に収めるため、２５５でandを取る
    unsigned int    i = fi.i & (SINTABLESIZE - 1);

    return  SinTable[i];
}

与えられた値に、TWOPISCALE　＝　２５６／６．２８３１８５３をかけて、２πが２５５となるよう正規化します。 その後、１．５×２＾２３を足し２５５でandを取って、整数化と適性範囲化を一度にしてしまいます。
上のfloat＜＝＞int変換では、０ｘ４ｂ４０００００で引き算しましたが、今、求める数は正の数ですから、 指数部が０になればいいので、２５５＝０ｘＦＦでandを取ることにより０にしています。
その値でテーブルを参照すれば、必要なsinの値が求まります。 cosは、−９０度ずれたsinと同じですから、最後のテーブル参照をする時に、配列要素が２５６の場合には６４を足して、 ２５５でマスクしてあげれば求まります。

このsinの求め方については、 Game Programming Gems2に詳しく出ていますので、 詳細が知りたい方はGem2を参照ください。 この記事には他にも、浮動小数の高速な比較や符号判定、絶対値、平方根や任意の関数の高速化技法が出ています。
この本、すごくいい本なんですが、高くて厚くて重いのが難点です。Gem3も出ているようですね。

精度と速度

この方法の精度は、３６０度の範囲を２５６や１０２４で割るわけですから非常に悪く、 配列範囲を２５６とすると、約１．４度位になります。 それでも、半径５０ｄｏｔの円の中に模様を描く、ＰｌａｎＣＥの惑星拡大画面のように、 １度の精度すら不要な場合も結構あり、使う場所に必要な精度さえ気にしていれば、使い道はたくさんあります。

ＰｌａｎＣＥの惑星拡大画面の場合、doubleの組み込みsinでは、１枚描くのに２秒ほどかかっていましたが、 このテーブル法を使うことによって、０.２秒まで高速化されました。実に１０倍！　 １．４ＧＨｚのＰｅｎｔｉｕｍ４でも、３倍ほどの高速化がなされました。

もちろん、精度がガタガタなのを忘れてはいけませんが、ライブラリは適材適所で使うのが一番ですね。

戻ります