Raspberry Pi QPU における逆数平方根計算と精度補正 #RaspberryPi

まえがき

これは昨日から今日にかけて行っていた活動なんですが，先を越されましたね: http://qiita.com/kaityo256/items/b28f95242a9690d4f896
先方の記事と被っている点が多々ありますがご了承ください．

はじめに

Raspberry Pi QPU における逆数計算と精度補正ではVideoCore IV QPUのSFUを用いて逆数を求め，精度補正した．
SFUには逆数の平方根を求める機能もある．今回はこちらを補正する．

環境

ハードウェアは Raspberry Pi 2 Model B を使用した．

ソフトウェアは社長が書いた py-videocore を使用して書いた．
コードはここに置いておきます: https://github.com/Terminus-IMRC/qpu-sfu-accuracy-correction

無補正の場合

まず無補正の精度の確認．

@qpu
def code_sfu_recipsqrt_0(asm):
    mov(r0, vpm)
    mov(sfu_recipsqrt, r0)
    nop()
    nop()
    mov(vpm, r4)

(1 << 23) (inclusive) 〜 (255 << 23) (exclusive) まで¹の一様乱数を16個生成し，その逆数をSFUで計算させた．

実行例は以下．

input: ['303f0346', '4ce817b2', '70ef3918', '302adb0e', '3db9455d', '283c7108', '7d644504', '6b804de1',
        '65f82de2', '324553ac', '674afaeb', '4e613483', '237f0fe7', '509d743f', '651bdd15', '11bf2e63']
input (float): ['6.9490e-10', '1.2168e+08', '5.9229e+29', '6.2157e-10', '9.0464e-02', '1.0461e-14', '1.8964e+37', '3.1022e+26',
                '1.4650e+23', '1.1486e-08', '9.5855e+23', '9.4458e+08', '1.3827e-17', '2.1133e+10', '4.6003e+22', '3.0163e-28']
output 0:   ['47143000', '38be1f00', '26bb4200', '471cb000', '4054cb00', '4b153000', '20878d00', '297fb100',
             '2c37dc00', '4611cc00', '2b8fc000', '38087900', '4d803d00', '36e6d000', '2ca40b00', '56517800']
output ref: ['47142ee1', '38be1d36', '26bb4264', '471cae3c', '4054c8ff', '4b1530cf', '20878d50', '297fb242',
             '2c37d972', '4611cb07', '2b8fbf9a', '3808789d', '4d803c31', '36e6d12d', '2ca40b00', '56517863']

結果を見ると，SFUの逆数平方根は逆数と同じく，下位8ビットが常に0になっている．その上の7ビット（！）も真値と一致しないものもある．逆数よりも精度が悪いですね．

前回同様，SFUの逆数平方根のVerilog実装と思われるものを見てみると，これも逆数と同じく，fn_interpから返ってきた値をそのまま返している．fn_interpは out_ms << 8 をreturnするため，SFUで求めた逆数平方根の下位8ビットが常に0なのは仕様のようだ．

補正公式

ニュートン法を用いる．すなわち

$$
x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}
$$

において $f(x)=0$ となる $f(x)$ を適当に定めてみる．
以下，$a$ の逆数平方根を考える．

試行1

単純に $f(x) = ax^2 - 1$ とおいてみる．このとき $f'(x) = 2ax$ であり，

$$
\frac{f(x)}{f'(x)} = \frac{x}{2}-\frac{1}{2ax}
$$

となる．

この$f(x)$だと$x$の逆数が出てきてしまい，単純には計算できなくなってしまう（SFUを2回呼び出すのは避けたい）．よって違う手を探す．

試行2

前回の手法を真似て， $f(x) = \frac{1}{ax^2} - 1$ とおいてみる．このとき $f'(x) = -\frac{2}{ax^3}$ であり，

$$
\frac{f(x)}{f'(x)} = -\frac{x(1-ax^2)}{2}
$$

となる．これならいけそう．
ニュートン法の式に代入して

$$
x_{n+1} = x_n + \frac{x_n(1-ax_n^2)}{2} = \frac{x_n}{2}(3-ax_n^2)
$$

綺麗になってくれた$x_{n+1}$君に拍手👏👏👏

補正1回

導出した補正公式を用いて，SFUで求めた逆数平方根を1回補正してみる．

@qpu
def code_sfu_recipsqrt_1(asm):
    mov(r0, vpm)
    mov(sfu_recipsqrt, r0)
    nop()
    nop()
    # 1
    mov(r2, 1.0).fmul(r1, r4, r4)
    fadd(r2, r2, 2.0).fmul(r1, r1, r0)
    fsub(r1, r2, r1).fmul(r2, r4, 0.5)
    fmul(vpm, r1, r2)

命令を綺麗に詰めることができた．美しいですね．add/mul ALUを頻繁に同時に使っており，small immediateを多用してる辺りポイント高いです．QPUプログラミングはこういうところが楽しいですね．
実行例は以下．

input: ['303f0346', '4ce817b2', '70ef3918', '302adb0e', '3db9455d', '283c7108', '7d644504', '6b804de1',
        '65f82de2', '324553ac', '674afaeb', '4e613483', '237f0fe7', '509d743f', '651bdd15', '11bf2e63']
input (float): ['6.9490e-10', '1.2168e+08', '5.9229e+29', '6.2157e-10', '9.0464e-02', '1.0461e-14', '1.8964e+37', '3.1022e+26',
                '1.4650e+23', '1.1486e-08', '9.5855e+23', '9.4458e+08', '1.3827e-17', '2.1133e+10', '4.6003e+22', '3.0163e-28']
output 0:   ['47143000', '38be1f00', '26bb4200', '471cb000', '4054cb00', '4b153000', '20878d00', '297fb100',
             '2c37dc00', '4611cc00', '2b8fc000', '38087900', '4d803d00', '36e6d000', '2ca40b00', '56517800']
output 1:   ['47142ee0', '38be1d36', '26bb4263', '471cae3c', '4054c8ff', '4b1530cf', '20878d4f', '297fb241',
             '2c37d973', '4611cb06', '2b8fbf9a', '3808789c', '4d803c31', '36e6d12d', '2ca40b00', '56517862']
output ref: ['47142ee1', '38be1d36', '26bb4264', '471cae3c', '4054c8ff', '4b1530cf', '20878d50', '297fb242',
             '2c37d972', '4611cb07', '2b8fbf9a', '3808789d', '4d803c31', '36e6d12d', '2ca40b00', '56517863']

かなり真値に近付いたが，最大$\pm 1$の誤差がある．
まあ許容できるが2回補正も行ってみる．

補正2回

1回目の補正で残った $\pm 1$ の誤差を消せないか，2回目の補正で試してみる．

@qpu
def code_sfu_recipsqrt_2(asm):
    mov(r0, vpm)
    mov(sfu_recipsqrt, r0)
    nop()
    nop()
    # 1
    mov(r2, 1.0).fmul(r1, r4, r4)
    fadd(r2, r2, 2.0).fmul(r1, r1, r0)
    fsub(r1, r2, r1).fmul(r2, r4, 0.5)
    fmul(r3, r1, r2)
    # 2
    mov(r2, 1.0).fmul(r1, r3, r3)
    fadd(r2, r2, 2.0).fmul(r1, r1, r0)
    fsub(r1, r2, r1).fmul(r2, r3, 0.5)
    fmul(vpm, r1, r2)

実行例は以下．

input: ['303f0346', '4ce817b2', '70ef3918', '302adb0e', '3db9455d', '283c7108', '7d644504', '6b804de1',
        '65f82de2', '324553ac', '674afaeb', '4e613483', '237f0fe7', '509d743f', '651bdd15', '11bf2e63']
input (float): ['6.9490e-10', '1.2168e+08', '5.9229e+29', '6.2157e-10', '9.0464e-02', '1.0461e-14', '1.8964e+37', '3.1022e+26',
                '1.4650e+23', '1.1486e-08', '9.5855e+23', '9.4458e+08', '1.3827e-17', '2.1133e+10', '4.6003e+22', '3.0163e-28']
output 0:   ['47143000', '38be1f00', '26bb4200', '471cb000', '4054cb00', '4b153000', '20878d00', '297fb100',
             '2c37dc00', '4611cc00', '2b8fc000', '38087900', '4d803d00', '36e6d000', '2ca40b00', '56517800']
output 1:   ['47142ee0', '38be1d36', '26bb4263', '471cae3c', '4054c8ff', '4b1530cf', '20878d4f', '297fb241',
             '2c37d973', '4611cb06', '2b8fbf9a', '3808789c', '4d803c31', '36e6d12d', '2ca40b00', '56517862']
output 2:   ['47142ee1', '38be1d36', '26bb4263', '471cae3c', '4054c8ff', '4b1530cf', '20878d50', '297fb242',
             '2c37d973', '4611cb06', '2b8fbf9a', '3808789c', '4d803c31', '36e6d12d', '2ca40b00', '56517862']
output ref: ['47142ee1', '38be1d36', '26bb4264', '471cae3c', '4054c8ff', '4b1530cf', '20878d50', '297fb242',
             '2c37d972', '4611cb07', '2b8fbf9a', '3808789d', '4d803c31', '36e6d12d', '2ca40b00', '56517863']

結果を見ると，2回目の補正で正しくなるものもあるが，$\pm 1$の誤差が残っているものもある．
また，3回補正しても2回補正と値が変わらなかった．

計算順序の改善

以下の順序で計算すると、$1 - a x_n^2 \simeq 1 - a \cdot (\frac 1 {\sqrt a})^2 = 1 - 1 = 0$ すなわち $x_n$ の誤差部分だけが残り、結果として $x_{n + 1}$ の精度が向上する (加算の回数は1回だけ増えてしまう)。

$$
x_{n+1} = x_n + \frac{x_n(1-ax_n^2)}{2}
$$

補正1回

TBA.

まとめ

SFUの逆数平方根は8〜15ビットの誤差を生じる．
ニュートン法で1，2回近似すると誤差を$\pm 1$以内にできる．

Verilog実装がどのような理論に従ってこのような計算（仮数部の一部のビット列を取り出して2乗してテーブルの値と掛け合わせるなど）を行っているのかがわからないので，詳しい方教えていただけると幸いです．

32ビット浮動小数点数に変換した際に0や非正規化数や無限大やNaNにならないように，指数部を1〜254にした．今回は妥協して正数のみ生成した． ↩