はじめに

Raspberry Pi には VideoCore IV という GPU があり，その中に QPU (Quad Processing Unit) というプロセッサがある．
QPU には SFU (Special Function Unit) という，32ビット浮動小数点数の逆数，平方根の逆数，exp2，log2 を求めることのできるユニットがある．この逆数がどのくらいの精度を持つのか確認し，精度補正してみる．

本記事は http://qiita.com/kaityo256/items/df33b64040794473db92 を大いに参考にしました．

環境

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

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

SFUの使い方

SFUは以下の流れで使用できる (VideoCore IV 3D Architecture Reference Guide p.23参照):

SFUのレジスタ (sfu_recip, sfu_recipsqrt, sfu_exp, sfu_log) に 32ビット浮動小数点数を 書き込む
2命令待つ
r4から計算結果を読み出す

無補正の場合

まず無補正の精度の確認．
(コードはコア部のみ示す．詳しい実装は実装コード参照)

def code_sfu_recip_0(asm):
    mov(r0, vpm)
    mov(sfu_recip, r0)
    nop()
    nop()
    mov(vpm, r4)

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

実行例は以下．

input: ['71ae7e4b', '6ade3b69', '77c35720', '1129289c', '391c8cd0', '55b1d5a5', '47384718', '5b09a410',
        '311bc955', '2bd53d28', '3c492ac2', '61f48c2a', '5335ee52', '1b43fb3d', '4bd201d0', '4d3bf5bc']
input as float: ['1.7281e+30', '1.3433e+26', '7.9239e+33', '1.3344e-28', '1.4930e-04', '2.4441e+13', '4.7175e+04', '3.8742e+16',
                 '2.2670e-09', '1.5152e-12', '1.2278e-02', '5.6389e+20', '7.8139e+11', '1.6211e-22', '2.7526e+07', '1.9709e+08']
output (0):   ['0d3bcb00', '14137400', '0727c000', '6dc1b800', '45d15100', '29384400', '37b1d200', '23ee1200',
               '4dd25700', '5319ac00', '42a2e500', '1d05ff00', '2bb41d00', '63a73400', '331c0a00', '31ae5800']
output (ref): ['0d3bca17', '14137310', '0727bf8e', '6dc1b630', '45d15031', '293842db', '37b1d18e', '23ee119f',
               '4dd256d7', '5319ab04', '42a2e3b8', '1d05fe92', '2bb41cc0', '63a73315', '331c0868', '31ae55c6']

結果を見ると，SFUの逆数は常に下位8ビットが0になっている．その上の2ビット程度も正しい値と一致してないことがあり，精度は悪い．

調べてみると²，SFUのVerilog実装をコードに落としたっぽいものが見つかった: https://github.com/broadcomCM/android_hardware_broadcom_brcm_usrlib/blob/ics/dag/vmcsx/middleware/khronos/glsl/2708/glsl_fpu_4.c#L1143
LUT (lookup table) 実装．
詳細は理解できていないが，fn_recip から呼ばれる fn_interp を見ると，関数の最後で out_ms << 8 をreturnしている．下位8ビットが0なのは仕様のようだ．

補正公式

http://qiita.com/kaityo256/items/df33b64040794473db92 をそのまま使用する．
すなわち，$a$の逆数の近似値を$x_n$，補正後の値を$x_{n+1}$とすると

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

である．

補正1回

とりあえずSFUで求めた逆数を1回補正してみる．

def code_sfu_recip_1(asm):
    mov(r0, vpm)
    mov(sfu_recip, r0)
    nop()
    nop()
    fmul(r0, r0, r4)
    fsub(r0, 2.0, r0)
    fmul(vpm, r4, r0)

実行例は以下．

input: ['71ae7e4b', '6ade3b69', '77c35720', '1129289c', '391c8cd0', '55b1d5a5', '47384718', '5b09a410',
        '311bc955', '2bd53d28', '3c492ac2', '61f48c2a', '5335ee52', '1b43fb3d', '4bd201d0', '4d3bf5bc']
output (0):   ['0d3bcb00', '14137400', '0727c000', '6dc1b800', '45d15100', '29384400', '37b1d200', '23ee1200',
               '4dd25700', '5319ac00', '42a2e500', '1d05ff00', '2bb41d00', '63a73400', '331c0a00', '31ae5800']
output (1):   ['0d3bca18', '14137310', '0727bf8d', '6dc1b630', '45d15031', '293842db', '37b1d18f', '23ee119f',
               '4dd256d6', '5319ab05', '42a2e3b7', '1d05fe92', '2bb41cc0', '63a73316', '331c0868', '31ae55c6']
output (ref): ['0d3bca17', '14137310', '0727bf8e', '6dc1b630', '45d15031', '293842db', '37b1d18e', '23ee119f',
               '4dd256d7', '5319ab04', '42a2e3b8', '1d05fe92', '2bb41cc0', '63a73315', '331c0868', '31ae55c6']

結果を見ると，1回補正によってかなり精度が上がっているが，値が最大$\pm 1$ズレている．
http://qiita.com/kaityo256/items/df33b64040794473db92 だと1回の補正でフル精度出てるので，ちょっと悲しい．

補正2回

1回だと$\pm 1$の誤差が出たので，2回補正してみる．

def code_sfu_recip_2(asm):
    mov(r0, vpm)
    mov(sfu_recip, r0)
    nop()
    nop()
    # 1
    fmul(r1, r0, r4)
    fsub(r1, 2.0, r1)
    fmul(r1, r4, r1)
    # 2
    fmul(r2, r0, r1)
    fsub(r2, 2.0, r2)
    fmul(vpm, r1, r2)

実行例は以下．

input: ['71ae7e4b', '6ade3b69', '77c35720', '1129289c', '391c8cd0', '55b1d5a5', '47384718', '5b09a410',
        '311bc955', '2bd53d28', '3c492ac2', '61f48c2a', '5335ee52', '1b43fb3d', '4bd201d0', '4d3bf5bc']
output (0):   ['0d3bcb00', '14137400', '0727c000', '6dc1b800', '45d15100', '29384400', '37b1d200', '23ee1200',
               '4dd25700', '5319ac00', '42a2e500', '1d05ff00', '2bb41d00', '63a73400', '331c0a00', '31ae5800']
output (1):   ['0d3bca18', '14137310', '0727bf8d', '6dc1b630', '45d15031', '293842db', '37b1d18f', '23ee119f',
               '4dd256d6', '5319ab05', '42a2e3b7', '1d05fe92', '2bb41cc0', '63a73316', '331c0868', '31ae55c6']
output (2):   ['0d3bca18', '14137310', '0727bf8e', '6dc1b630', '45d15031', '293842db', '37b1d18f', '23ee119f',
               '4dd256d6', '5319ab05', '42a2e3b8', '1d05fe92', '2bb41cc0', '63a73316', '331c0868', '31ae55c6']
output (ref): ['0d3bca17', '14137310', '0727bf8e', '6dc1b630', '45d15031', '293842db', '37b1d18e', '23ee119f',
               '4dd256d7', '5319ab04', '42a2e3b8', '1d05fe92', '2bb41cc0', '63a73315', '331c0868', '31ae55c6']

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

計算順序の改善

補正の式は以下と変形できる。

$$
x_{n+1} = x_n (2 - ax_n) = x_n + x_n (1 - a x_n)
$$

このとき、$1 - a x_n \simeq 1 - a \cdot \frac 1 a = 1 - 1 = 0$ すなわち $x_n$ の誤差部分のみが残るため、結果として先程の式よりも $x_{n+1}$ の精度が向上する。
この式を実装し、補正をしてみた。

補正1回

補正2回

まとめ

SFUの逆数は8〜10ビットの誤差を生じる．
ニュートン法で1，2回近似すると誤差が2ビット以下になる．

今後はVerilogのLUT実装を解読し，また，近似手法について知見を深めたいです．

32ビット浮動小数点数に変換した際に0や非正規化数や無限大やNaNにならないように，指数部を1〜254にした．今回は妥協して正数のみ生成した． ↩
というよりも，今回の精度補正の試みはこのコードを見つけたことから始まった． ↩