公開日：2023/02/18

24th Season

#24-6　浮動小数の開平

　前回、オリジナル32ビットCPUコアNihonium用として設計した単精度浮動小数演算ユニット(S-FPU)を使う応用として、平方根を計算するプログラムを作ってみようと思います。
　いくつかアルゴリズムを検討した上で、最終的に実用的な方法を採用したいと思います。高速かつ、精度良く平方根を計算することを目標とします。

　まずは、以下の記事で使った、整数開平アルゴリズムから検討して行きます。

　　「#23-1　整数演算で開平してみよう」

二分探索

　以前(「#23-1　整数演算で開平してみよう」)、浮動小数演算を使わずに整数のまま開平🚪する事はできない物か？と、考えて採用した方法が、二分探索でした。
　まずは、この方法を浮動小数計算に応用することを考えてみましょう。

　二分探索は、求める値√xを、L≦√x≦Rの範囲に絞り込んでおいて、M＝(L＋R)÷2と、√xの大小を比較しながら、1探索当たり1ビットを確定して行く方法です。この方法で整数xを開平した時(「#23-1　整数演算で開平してみよう」)は、√xを求めるのにlog₂ (x)回の探索が必要でした。

　今回、浮動小数(下図)で開平するためには、どのくらいの探索回数が必要でしょうか？

　二分探索では、探索1回当たり1ビットの数値が確定して行きます。単精度浮動小数float型の仮数部は24ビットあるので、探索回数は24回必要となります。
　探索1回当たり、以下の式で大小関係を比較する処理が発生します。

　　y＝√x　……　両辺二乗
　　y²＝x　……　左辺のyを移項
　　y＝x÷y

　y²＝xの式で大小比較すると、一般に遅いとされるコンピュータ除算を使わずに判定が可能です。但し、yを二乗することで、整数演算であれば整数の最大値を超えたり(こちらを参照)、浮動小数であっても、二乗した結果オーバフロウすると無限大INFという特殊な値になってしまうので、正常に計算できなくなってしまいます。
　もう一段変形した、y＝x÷yの式で判定するとコンピュータ除算が必要となりますが、オーバフロウすることはないので、正しく計算することができます。

　つまり、24ビットの仮数部を求めるためには24回の探索が必要で、1探索当たり除算1回と比較1回の演算が必要です。
　まだ他の方法を検討していないので定かなことは言えないのですが、感覚的に処理コストが大きく、演算速度が遅い方法なのではないでしょうか？

テイラー展開

　以下の記事で、半精度浮動小数演算ユニットH-FPUを使って数学関数を整備しました。

　　「#16-9　数学関数を整備しよう」

　上の記事中では開平sqrt()関数を作らなかったのですが、以下の式で演算可能であると、ご紹介だけしておきました。

　指数exp()関数も自然対数log()関数も、テイラー展開を使って近似値を求めることができます。「#16-9　数学関数を整備しよう」に於いて、半精度浮動小数half型で実装した時は、以下の項数が必要でした。

指数exp()関数：

half型では、テイラー展開の項数n＝0～4。

自然対数log()関数：

half型では、テイラー展開の項数n＝0～2。

　今回のfloat型では、仮数部のビット数が24ビットに増加しているので、もう少し多くの項数まで計算しないと、必要な精度が得られません。ちゃんとした評価はしていませんが、ザっと見積もったところ、以下の項数が必要だと考えられます。

指数exp()関数：

float型では、テイラー展開の項数n＝0～7。

自然対数log()関数：

float型では、テイラー展開の項数n＝0～6。

　2個の関数合わせて、およそ15項ほどの演算です。テイラー展開の演算は浮動小数の積和演算となって、1項当たり乗算2回と加算1回がコストです。二分探索を使うよりは、速いかな？って感じですね。

　計算に関数2個分の時間を要する(遅い)だけではなく、近似値を重ねることで演算結果に誤差が蓄積する(精度が低下する)ことも不利な点として挙げられます。

下準備

　どちらも今一つなので、更に他の方法を探ってみましょう。

　何にしても、まず浮動小数float型の特徴を使って、求めるべき数値xの範囲を限定しておくことは重要です。テイラー展開にしても、xの絶対値を小さくする工夫をしてから、級数に適用しています(例えばこちら)。

　float型のバイナリ値は以下図のような形式です。

　平方根√xは、実数の範囲ではx≧0だけで定義されるので、符号部Sx＝0(正)と考えて差し支えありません。

　　x＝Fx×2^Ex

　√xを代数計算で、できるところまで展開してみましょう。

√x＝√(Fx×2^Ex)

xの定義域は、0≦x＜∞(INF)

ここで、Exの値で場合分けする。

Exが偶数の場合：

Ex＝2×Exkとおく。ただし、Exkは整数。
√x＝√(Fx×2^2×Exk)＝√Fx×2^Exk
2の指数項が開平できた(平方根の外に出た)。
平方根の中に残った数値をx’と置くと、

1≦x’＝Fx＜2である。
float型の数値xを使うと、xの指数部をEx’＝0x7f(Ex＝0)に偽装すればx’となる。

Exが奇数の場合：

Ex＝2×Exk＋1とおく。ただし、Exkは整数。
√x＝√(Fx×2^2×Exk＋1)＝√(2×Fx)×2^Exk
2の指数項が開平できた(平方根の外に出た)。
平方根の中に残った数値をx’と置くと、

2≦x’＝2×Fx＜4である。
float型の数値xを使うと、xの指数部をEx’＝0x80(Ex＝1)に偽装すればx’となる。

定義域：1≦x’＜4の平方根ｙ’＝√x’を求めると、以下式によって、元のy＝√xを知ることができる。

y＝√x＝√x’×2^Exk＝y’×2^Exk
float型の数値y’を使うと、×2^Exkは、y’の指数部ExにExkを加算することと等価である。

ここで、－63≦Exk≦63である。
1≦y’＝√x’＜2なので、y’の指数部Ex＝0。
∴－63≦Ex＋Exk≦63である。この演算で指数部がオーバフロウやアンダフロウすることはない。

　結果、1≦x’＜4範囲のy’＝√x’を求める問題に帰着できました。定義域も狭まったし、絶対値も小さくなりました。仮に、二分探索やテイラー展開へ適用する場合にも有利です。

Newton法

　ようやく真打ち登場です。
　Newton法とかNewton-Raphson法とか呼ばれている方法で、漸化式を使って、求める値に漸近して行く計算方法です。Newton法で開平する方法として、順当に√xを求める漸化式と、間接的に1/√xを求める漸化式が存在します。

√x：

漸化式：　y<n＋1>＝(y<n>＋x÷y<n>)÷2

n番目の数値yを、y<n>と記載する。nが定数の場合、y<0>の代わりにy0等と記す。

本来は、√x＝y<∞>だが、必要な精度が得られた時点で止めて良い。

ただの反復計算なので、Excel等の表計算でも、簡単にシミュレーションできる。

実際に試してみると、Newton法の漸化式は、ものすごく収束が速い。
float型程度の精度であれば、3～4回反復すると、十分な精度を得られる。

反復1回当たり、浮動小数除算1回と加算1回がコスト。
最後の÷2は、整数型ならシフト演算だが、float型の場合は指数部Ex－1で良い。

前節の下準備によって、1≦x＜4より、1≦y＝√x＜2なので、yの指数部Ex’＝0x7f(Ex＝0)付近であるはず。

－1くらいでは、オーバフロウもアンダフロウもない。

従って、int．y－0x0080_0000の整数演算で、÷2が計算できる。

int．yは、yの小数点以下を切り捨てて整数化する事ではなくて、float型yのバイナリ値を、そのまま32ビット整数値と見なして整数演算する事を示す。

注意）加減算によって指数部がオーバフロウまたはアンダフロウする場合、単純に加減算してはならない。

オーバフロウしたら±INF、アンダフロウしたら±0に置き換える必要がある。

1/√x：

漸化式：　y<n＋1>＝y<n>×(3－x×y<n>×y<n>)÷2
y＝1/√xと置いて、右辺を有理化すると、y＝√x÷xより、√x＝y×xで求められる。

用途次第では1/√xのまま使用する場合もある。
例）ベクタV＝(vx，vy，vz)の正規化

Va＝vx²＋vy²＋vz²
vx＝vx÷√Va＝vx×(1/√Va)
vy＝vy÷√Va＝vy×(1/√Va)
vz＝vz÷√Va＝vz×(1/√Va)
Newton法で1/√Vaを求めておいて、(vx，vy，vz)の各成分に乗ずれば良い。

一般的に低速なコンピュータ除算を使わず、一貫して乗算と加減算だけで決着できることが特長。
反復1回当たり、浮動小数乗算3回と減算1回がコスト。
最後の÷2は、予め丸カッコの中に入れて、3÷2＝1.5(定数)およびx÷2としても良い。

　1/√xの漸化式は、一般的に低速とされるコンピュータ除算を使わずに高速演算できる事が謳い文句です。ただし、反復1回当たり3回の乗算が必要となります。それでも「普通は」、1/√xの漸化式を使った方が速いと思います。

　しかし、こと、Nihonium＋S-FPUに関しては、こちらの表にまとめた様に、乗算fmul命令×3回より、除算fdiv命令×1回の方が速い(※)のです（１＾＾；）。よって、普通に√xの漸化式を使います。

※S-FPUはCPU(Nihonium)の外部にコプロセッサとして接続するユニットのため、CPU⇔S-FPU間のデータ転送がクリティカルパスとならない様、往復2サイクルたっぷり確保して、タイミングルースな設計になっています。このオーバヘッドが、演算サイクル数の小さい乗算にとっては大きく響いて、この様になっています。

漸化式の初期値y0

　漸化式なので、初期値y0の値を決めなくてはなりません。漸化式は以下の通りなので、初期値y0はアルゴリズム上、y0＞0ならば何でも構いません。

　　y<n＋1>＝(y<n>＋x÷y<n>)÷2

　y0＝1だろうが、y0＝234だろうが、反復する度に求める値に漸近して行きます。ただ、適切なy0を選ばないと、収束するまでの反復回数が増加します。高速に計算するためには、なるべく反復回数を少なくしたいので、適切なy0を使う必要があります。

　　最適解は、y0＝√xです。

　と、言っても、そもそも√xを求めたくて、この方法を使うのですから、「鶏と卵」になっているわけです。

　ここで、前の節の下準備が活きてきます。下準備によって1≦x＜4の範囲に追い込んであるので、1≦√x＜2より、例えばy0＝(1＋2)÷2＝1.5(定数)に固定しても、十分高速に収束させることができます。

　もう一歩踏み込んで見ましょう。1≦x＜4の範囲で√xのグラフを描画すると、以下図となります。

　本来、y＝√xのグラフは曲線(上図、青色プロット)なのですが、xの範囲が狭いので、無理やり一次回帰直線(上図、赤点線)で近似してしまいます。

　　y0＝√x≒0.3259×x＋0.7407、ただし1≦x＜4

　Excelを使ったシミュレーションでは、y0＝1.5(定数)固定の場合に比べて、上式でy0を決定した場合、Newton法の反復1回分削減することができました。
　計算量も少ない(乗算1回と加算1回)し、反復の初期値y0としては申し分ありません。

プログラム

　ここまでの節で検討して来た、開平アルゴリズムの全体像を以下にまとめます。float型のバイナリデータを整数と見なして処理するなど、少しトリッキィな部分もありますが、全て本ページで解説してありますので、「ん？」となったところは、上の方に戻って復習してみてください。

　　開平sqrtルーチン：

if　x＝±0：　return　x
if　x＜0 || x＝NaN：　return　NaN
if　x＝INF：　return　INF
Ex’＝(int．x>>23)＆0xff
Ex＝Ex’－127
if　Ex＆0x01＝0：　// 偶数の場合

Exk＝Ex>>1
x’＝int．x＆(～0x7f80_0000)｜0x3f80_0000

if　Ex＆0x01＝1：　// 奇数の場合

Exk＝(Ex－1)>>1
x’＝int．x＆(～0x7f80_0000)｜0x4000_0000

// この時点で、1≦x’＜4、－63≦Exk≦63である。
y’＝√x’をNewton法で求める。

y0＝0.3259×x’＋0.7407
repeat　3：　// n＝0～2

y’＝y<n＋1>＝(y<n>＋x’÷y<n>)÷2
　　＝int．(y<n>＋x’÷y<n>)－0x0080_0000

return　y＝√x＝√x’×2^Exk＝int．y’＋(Exk<<23)

※int．xは、float型xの小数点以下を切り捨てて整数化する事ではなく、float型xのバイナリ値を、そのまま32ビット整数値と見なして整数演算する事を示します。
※int．xの演算は、float型がIEEE754 binary32形式(除く非正規化数、かつビッグエンディアン)であることを前提としています。Nihonium＋S-FPUでは全く問題ありませんが、他の環境で実装する場合には注意してください。

　以下に評価用プログラムをリンクしておきますので、具体的なコーディングについては、以下のプログラムを参照してください。

　　完成プログラム：niho_sfpu_sqrt.asm

　Nihonium用アセンブラについては、こちらをご参照ください。アセンブルには、最新版のOlive+が必要となります。ダウンロードコーナなどからダウンロードしてWindowsにインストールしてください。

実験はシミュレーションで

　前回設計したNihonium＋S-FPUを搭載したNHOC113FS0408等のデザインを以下リンクからダウンロードしてください。

　今回は、作成した開平sqrtルーチンを評価するだけですので、FPGAボードには書き込みません。iverilogを使ったシミュレーションで評価します。

　前節のプログラム(niho_sfpu_sqrt.asm)をアセンブルして、シミュレーションしましょう。

　　実験サンプル：【】内はCygwinシェルウインドウからの実行コマンド
　　【iverilog -f test.lst ; if ( ! $status ) vvp a.out】
　　sqrt
　　：
　　PLOT
　　00000000 00000000
　　3da14285 3e8fabb0
　　3e214285 3ecb2e56
　　3e71e3c7 3ef8d84a
　　：
　　：　(中略)
　　：
　　411c3870 4047fb28
　　411d7af5 4048c92d
　　411ebd7a 4049965f
　　411fffff 404a62c1
　　END
　　Test Pass
　　CPI:　15254/　8320(　4262)=　1833/1000
　　finish:　611720 ns

　ログのPLOT～ENDの間に表示された8桁の16進数×2個は、1個目が横軸xの値、2個目が縦軸y＝√xの値です。16進数はfloat型のバイナリデータなので、こちらの図に従って浮動小数値に変換してグラフを描くと、以下図の様になりました。

　青色プロット(左軸)がNewton法で計算したy＝√x、赤色プロット(右軸)はExcelのsqrt(x)関数を期待値とした場合の誤差です。

　　考察：