近似計算について
SSEのところで現れた逆数や逆数平方根の近似式について少し考えてみます.
No. 001
逆数平方根アルゴリズムは次の通りでした.
y = aの逆数平方根の近似値
z = 3/2 * y - 1/2 * a * y3
これはどの程度の精度があるのでしょうか.
まず誤差をどう表すかを考えます. r = 1 / √a を正しい値とします.
単純に |y - r| を誤差としても構いませんが,
a が 0 に近いときや大きいときでも同じように評価するために
ε = (y - r) / r
を誤差と定義することにします. 例えば y が11bitまで正しければ ε は
2-11 程度になります.
y = r + rε と変形して z に代入します. a = r-2に注意
(r + rε)3 = r3 + 3 r2(rε) + 3 r (rε)2 + (rε)3
となりますが r に比べて rε は十分小さいので3乗の項は無視します.
z = 3/2 * (r + rε) - 1/2 * r-2 * (r3 + 3 r2 rε + 3 r (rε)2)
= 3/2 * (r + rε) - 1/2 * (r + 3rε + 3 rε2)
= r - 3/2 * rε2
ここで求めたzとrとの誤差を考えると
|(z - r) / r| = 3/2 * |ε|2
誤差はε2程度の大きさになることが分かります.
つまり y が11bitまで正しければ z は倍の22bit程度
(実際には途中の演算誤差があるので多少はそれよりも落ちますが)
保証されることがわかりました. SIMD計算で使うには十分であることが分かります.
同時にこの式を繰り返すとどんどん精度が良くなっていくことも分かります.
| y = x(赤) と y = 3/2x - 1/2ax3(緑)のグラフ |
x座標の点(an,0)を上に上っていって緑のグラフと交わった点の y 座標が 3/2an-1/2a(an)3になります. そこから右に移動して赤の線と交わったところを降りてきた点をan+1とします. anよりはan+1の方が1/√aに近くなっていること, これを繰り返すと交点に近づくことが予想できます.
さてmp3エンコーダでは非線形量子化のために0.75乗を計算することが多いです.
無論これを高速に求めたいものです.
現在は
aの近似逆数平方根 → 精度向上で1/√a → 1/√aの近似逆数平方根 → 精度向上でa0.25 →3乗してa0.75
という方式で計算しています.
これでも既にC言語レベルに比べて10倍以上高速なのですが,
精度向上を1回で済ませることが出来ないかを考えてみます.
a に近似逆数平方根を2回, 近似逆数を1回して近似 a-0.25 を求めた後,
1回精度向上をすればいけそうな気がします.
その計算式と精度の検証についてはまた次回.
No. 002
少し回り道をしてどんな近似式になるのか推測してみます.
今までのを振り返ると
| y = x-1/2 → y = 3/2x - 1/2ax3 |
でした. これから一般にaの逆数 c 乗根の近似式は
| y = x-1/c → y = px - qaxr, p, q, rは不明 |
という形だろうと仮定してみます.
p, q, rの条件はなんでしょうか.
まず上のグラフから考えるとこの式とy = xの交点が
(a-1/c, a-1/c) でなければなりません.
x = p x - q a xr を変形すると p - 1 = q a xr - 1.
これにx = a-1/cを代入して p - 1 = q a(1 - (r-1)/c) ---(*)
ちょっと分かりにくいですが, p や q が a
に依存して欲しくないのでa(1 - (r-1)/c)は定数になって欲しいです.
よって1 - (r-1)/c = 0, つまり r = c + 1 となります.
このとき(*)に戻って考えると p - 1 = q.
更に戻って y = (q+1) x - q a xc+1 という形になりました. ---(**)
q を決定しましょう.
この式がもとの近似値より精度よくなって欲しいわけですからその条件を調べます.
真の値を α = a-1/c, 近似値を x = a-1/c + ε とします.
x を(**)に代入すると
(q + 1)(α + ε) - q a(α + ε)c + 1.
ここで(α + ε)c + 1 = αc + 1 + (c + 1)αcε + ...
であることと ε2以降は無視できる大きさであることを考えて`...'の部分は切り捨てます.
すると元の式
≒ (q + 1)α - q a αc + 1 + (q + 1 - q a(c + 1)αc)ε
= (q + 1)a-1/c - q a a(-(c+1)/c) + (q + 1 - q(c+1))ε (がんばるんだ...(^^;)
= α + (q + 1 - q(c + 1))ε
となりました.
もとの値より良くなって欲しいのですから ε の係数 = 0, つまり
よって q = 1/c となりました.
一番最初に戻ると結局
| y = x-1/c → y = (1 + 1/c)x - (1/c)a xc+1 |
と結論できます. c = 1 や 2 のとき逆数, 逆数平方根の近似式になることを
確認してください.
こんな泥臭いやり方でも求められるものですね.
No. 003
本質的な内容は同じですがグラフから攻める方法も考えてみましょう.
α = a-1/c は y = x-c - a の根です.
近似値 an での接線を考えます.
y' = -c x-c-1 ですから, 接線は y = -c an-c-1 (x - an) + an-c - a.
グラフを書いてみましょう.
| y = x-c - a(赤) と x = an での接線(緑)のグラフ |
この緑の線と x 座標との交点 an+1 は α により近くなっていますし,
この作業を繰り返すとどんどん α に近づくのも想像できます.
an+1 は
0 = -c an-c-1 (x - an) + an-c - a
となる x です.
つまり(ちょっとがんばって計算すると)
an+1 = (1 + 1/c) an - (1/c) a anc+1
となります. なんとこれは上で求めた式と同じです!
これがニュートン法と呼ばれる近似計算方法です.
No. 004
では実際に適用してみましょう. 最初の近似値を求めるには
No. 011を利用します.
#define BITLEN 8
#define SIZE (1<<BITLEN)
static float idxTbl[256], powTbl[SIZE];
void init_powM0p25(void)
{
int i;
const float A = -0.25;
for (i = 0; i < 256; i++) {
idxTbl[i] = pow(2, (i - 127) * A);
}
for (i = 0; i < SIZE; i++) {
powTbl[i] = pow(1 + i / (double)SIZE, A);
}
}
/* x^(-0.25)を求める. x >= 0であること */
float powM0p25(float x)
{
int s, t;
float y;
s = t = *(int *)&x;
s >>= 23; /* 指数部 */
t = (t & ((1<<23)-1)) >> (23 - BITLEN); /* 仮数部 */
y = idxTbl[s] * powTbl[t];
return (1.25 - 0.25 * x * (y * y) * (y * y)) * y;
}
PIII-700で powM0p25(x) * x と sqrt(sqrt(x) * x) を比較したところ実行時間は100clk : 170clk 程度でした.
SIMD命令では return 命令の直前までを3命令(約6clock)で行い,
かつテーブルも不要ですからかなり高速化されます.
なお, idxTbl[s] というテーブル引きの代わりに
| s-0.25 = ((s & 3) + (s & ~3))-1/4 = (s & 3)-1/4 * (s >> 2)-1 |
を利用してテーブルサイズを小さくする方法もあります.
例えば
float pow_fast1(float x)
{
int s, t;
const float exponentTbl[] = {
1.00000000 / (1UL<<30),
1.18920712 / (1UL<<30),
1.41421356 / (1UL<<30),
1.68179283 / (1UL<<30)
};
float y;
s = t = *(int *)&x;
s >>= 23; /* 指数部 */
t = (t & ((1<<23)-1)) >> (23 - BITLEN); /* 仮数部 */
s = 247 - s;
y = exponentTbl[s & 3] * (1UL<<(s>>2)) * powTbl[t];
return (1.25 - 0.25 * x * (y * y) * (y * y)) * y;
}
こんな感じ.
ただ x の値によっては (1UL<<(s>>2)) が overflow するので使い方に注意が要ります.
メモリアクセスのコストと若干の論理演算のコストとのどちらが重いかはケースバイケースでしょう.
以下は午後のこ〜だのソース(quantize.nas)の pow075_SSE の一部抜粋です.
0.75乗以外の処理も含めて1サンプルあたり約6clockで実行します.
;------x^(3/4)------
; input xm1 = x
rsqrtps xm4, xm1
movaps xm5, xm1
rsqrtps xm4, xm4
mulps xm5, xm7 ; -x/4
rcpps xm4, xm4 ; y = approx. x^-0.25
movaps xm6, xm4
mulps xm4, xm4 ; y^2
mulps xm1, xm6 ; yx
mulps xm5, xm4
mulps xm5, xm4
xorps xm4, xm4
addps xm5, [Q_f1p25] ; 5/4 - 1/4 * x (y^4)
mulps xm1, xm5
maxps xm0, xm4 ; if (x == 0) then x^(3/4) = 0
maxps xm1, xm4
御意見は光成滋生<herumi@nifty.com>までお願いします