落胆がらくた街

多分三日坊主で終わると思うんですけど

UINT32正規化版高速sinは何をしてるのか

suzublog.hatenablog.com
この記事でやってる、0~2PIをUINT32に正規化してやると高速でSinを計算できるという話。
ようやくこいつが何をしてるのか理解できたので忘れない内に書いとく。

元ソース。

float fastsin(UINT32 phase)
{
    const float frf3 = -1.0f / 6.0f;
    const float frf5 = 1.0f / 120.0f;
    const float frf7 = -1.0f / 5040.0f;
    const float frf9 = 1.0f / 362880.0f;
    const float f0pi5 = 1.570796327f;
    float x, x2, asin;
    UINT32 tmp = 0x3f800000 | ( phase >> 7 );
    if (phase & 0x40000000)
        tmp ^= 0x007fffff;
    x = ( *( (float*)&tmp ) - 1.0f ) * f0pi5;
    x2 = x * x;
    asin = ( ( ( ( frf9 * x2 + frf7 ) * x2 + frf5 ) * x2 + frf3 ) * x2 + 1.0f ) * x;
    return ( phase & 0x80000000 ) ? -asin : asin;
}






さて、sinを計算機で求めるってのは、精度が保証されるまでの次数でマクローリン展開を使って近似する事を意味する。

f:id:suzumodoki:20180212181224p:plain
...desmos兄貴は便利だなぁ。

赤線が理想的なsin(x)で、青線は9!まで使ったマクローリン展開*1
マクローリンは項数を増やすほど遠くの値(より大きなx)まで近似できるが、どうせsinはループするので、現実的な項数で十分高精度な近似が可能…という理屈。
このUINT版sinも、結局マクローリンを計算するだけだ。
ここでsinの周期性を考える(下図)と、本当に必要なのはsin(0)からsin(π/2)までである事がわかる。

f:id:suzumodoki:20180212182447j:plain





で、位相θをUINT32にに纏めるってのは、こういう対応付けを意味する。

実数θ UINT32 (2進数表記) UINT32 (16進表記)
0 0 0
π/2 0100_0000_0000_0000_0000_0000_0000_0000 0x40000000
π 1000_0000_0000_0000_0000_0000_0000_0000 0x80000000
3π/4 1100_0000_0000_0000_0000_0000_0000_0000 0xc0000000
0000_0000_0000_0000_0000_0000_0000_0000(見えない33ビット目が1) 0x00000000

ここで重要なのはアタマの2ビットだ。2進数の仕組みを考えれば当たり前の事を言っているだけだが、図で示すとこうなる。

f:id:suzumodoki:20180212183943j:plain
さっきの「π/2までわかればOKの法則」に照らすと、

  • 実質要るのは下30ビット(0~π/2に対応)の値
  • 1番上のビットが立ってたら全体を-1倍
  • 上から2番めのビットが立ってたら軸対称に入れ替え

という事になる。おぉ、なんだかそれっぽくなってきた。
ここで「軸対称に入れ替え」はビット演算においては単なるビットの反転で良い。π/2で対象という事は、π/2を超えた分だけπ/2から引いてやればいいわけだが、 π/2が今0b1000...と言ってるのだから、ここから引くということはπ/2より下のビットを反転する事になる。

さて、肝心のマクローリン展開級数はfloat型で出さなきゃならんので、uint型のphaseをfloat型に直さなきゃいけない。正規化の逆の手順を踏むわけだ。 ここで

float x = (PI/2.0f) * (float)phase / (0x3fffffff);

とかやったら台無しなので、float型変数のビットを直で弄るというカッコイイ技を使う。C++の真骨頂だ。


改めてソースを読み解こう。

float fastsin(UINT32 phase)
{
    const float frf3 = -1.0f / 6.0f;
    const float frf5 = 1.0f / 120.0f;
    const float frf7 = -1.0f / 5040.0f;
    const float frf9 = 1.0f / 362880.0f;
    const float f0pi5 = 1.570796327f;
    float x, x2, asin;
    ...略...

変数定義部分では、後のマクローリン展開で必要になる項を予め定数として算出している。即ち、
 \displaystyle
sin(x)=\sum_{k=1}^{N} (-1)^{k-1}\frac{x^{2k-1}}{(2k-1)!}
のxでない部分だ。f0pi5は…なんだろう?ネーミングの由来が分からん。 こいつはPI/2を意味していて、UINTの下30ビットを逆正規化するときの定数だ。なんでこの名前なんだろう…
微妙に追記:frfはおそらくfloat reciprocal(=逆数) factor(=階上)でf0pi5はfloat 0pi5の事だと思った。そういや英語だと0.5piを0pi5とか書くらしいですね
次に、

UINT32 tmp = 0x3f800000 | ( phase >> 7 );
    if (phase & 0x40000000)
        tmp ^= 0x007fffff;
    ...略...

0x3f80....は2進数で0011_1111_1000....だが、ここでfloat型のビットの意味を考える。

浮動小数点数型と誤差

一番上のビットは符号ビットなので、正を意味する。次の8ビット"011_1111_1"は10進数で127だから、Floatの指数部は2の0乗で1、となる。 頭でやってる

tmp = 0x3f800000 | ( phase >> 7 )

は、floatは仮数部が23ビットなので、UINT32の内どうせ入らない下の7ビットを捨てて仮数部としてぶち込む事を意味する。

次のif文は、よく見ると単純だ。0x40....というのは上から2番目のビットが1という事だから、「π/2の向こう側にいる場合」を判定している。 そしてその場合はビットを反転するのだった。0x007fffffとXORするということは、tmpの下23ビットを反転するという事だが、これは結局floatにおける仮数部を反転しているわけだ。

 x = ( *( (float*)&tmp ) - 1.0f ) * f0pi5;
    x2 = x * x;
    asin = ( ( ( ( frf9 * x2 + frf7 ) * x2 + frf5 ) * x2 + frf3 ) * x2 + 1.0f ) * x;
    return ( phase & 0x80000000 ) ? -asin : asin;
}

そして荒業、「UINT32領域でデザインしたtmpをポインタ経由でfloatとして扱う」作戦。これはあくまで同じビット構成を持つfloat型に変えるという事で、いわゆるキャストとは全く違う*2C/C++ではこういうカッコイイ技がさらっと使えて痺れる。

その後この自作floatから1を引いているのは、floatがそもそも「ケチビット」を付与するよう出来ているからだ(上のリンク参照)。 この頼まなくてもついてくる1を引けば、割り算やキャストを介さなくても「π/2に対するUINT32の位相の割合」が出る。割合だから、π/2(変数名はf0pi5)をかければ、逆正規化出来上がり。 あとは普通にマクローリン項を計算して、最後にphaseの頭、つまり符号ビットを解決してあげれば完成だ。




うーん、分かってみると単純だが、考えた人すごいなぁ。実際めちゃ速いしこれ。ちなみに、最後のphaseの最上位ビット見る処理はifを使わなくても

*( (UINT32*)&asin ) |= ( phase & 0x80000000 );
    return asin;

とかやれば同じことが出来るが別に速くならなかった。コンパイラが賢くやってくれてるのか、ポインタキャストもそんなに速くないのか、言うほど条件分岐が重い処理じゃないのか。