まえがき
以前投稿した記事で丸め誤差を簡単に紹介しました。
本記事では，高精度計算の準備としてエラーフリー変換を紹介します．
エラーフリー変換とは
浮動小数点数を用いた計算では多くの場合誤差が生じることを以前の記事で理解していただけたと思います．誤差は生じてしまうけれど，その誤差情報も計算結果として保持しておこう！というのがエラーフリー変換の考え方です．つまり，以下の数式のように2つの浮動小数点数 a, b の演算結果を近似解 x と誤差項 y に変換します．
a \oplus b = x + y \\
これをエラーフリー変換といいます．
準備
以降で和と積に関するエラーフリー変換と高精度計算方法について紹介していきますが，
その前に紹介中で用いる記号や関数・言葉について紹介しておきます．
また，簡単のため演算中でオーバーフローやアンダーフローは起きないと仮定します．
- $\mathbb{F} \qquad \ $：浮動小数点数の集合
- $\mathrm{u} \qquad \ $：単位相対丸め（float型なら $2^{-24}$, double型なら $2^{-53}$）
-
FMA$\ \ \ ：$ 融合積和演算（Fused-Multiply-Add） $xy+z$ を 1回の丸めで計算できる機能. - flop $\quad \ ：$ 浮動小数点演算回数（Floating-point Operation）
和に関するエラーフリー変換
和に関するエラーフリー変換として Dekker12 に紹介された FastTwoSum というアルゴリズムがあります．
FastTwoSum
$a$, $b \in \mathbb{F}$ について $|a|\geq |b|$ のとき $\rm{FastTwoSum}$ を実行すれば,
a + b = x + y \
が成り立つ.
def FastTwoSum(a, b):
x = a + b
y = b - (x - a)
return [x, y]
FastTwoSum は加算する浮動小数点数の大小関係に制限があるものの 3 flop（通常の和の3倍）で計算できます．また，加算する浮動小数点数の大小関係に仮定を置かないバージョンとして Knuth3 の TwoSum というアルゴリズムもあります．
TwoSum
$a$, $b \in \mathbb{F}$ について$\rm{TwoSum}$ を実行すれば,
a + b = x + y \
が成り立つ.（FastTwoSum の結果と同じ）
def TwoSum(a, b):
x = a + b
tmp = x - a
y = (a - (x - tmp)) + (b - tmp)
return [x, y]
TwoSum では加算する浮動小数点数の制限はなくなりましたが 6 flop（通常の和の6倍）必要になります．しかし，多くの場合で加算する浮動小数点数の大小関係は非自明なので TwoSum を使用するケースが多いです．
積に関するエラーフリー変換
積に関するエラーフリー変換として, またまた Dekker1 により紹介された TwoProduct というアルゴリズムがあります．このアルゴリズムは FastTwoSum/TwoSum に比べると少し複雑ですが以下のようになります．
TwoProduct
$a$, $b \in \mathbb{F}$ について$\rm{TwoProduct}$ を実行すれば,
a * b = x + y
が成り立つ.
def split(a):
# s は浮動小数点数の (仮数部 bit 数 t // 2 ) + 1 で計算できる定数
# 倍精度なら t = 53, 単精度なら t = 24
t = 53
s = (t // 2) + 1
tmp = a * (s + 1)
tmp = s * a
ah = tmp - (tmp - a)
al = a - ah
return [ah, al]
def TwoProduct(a, b):
x = a * b
[ah, al] = split(a)
[bh, bl] = split(b)
y =(al * bl) - ((x - (ah * bh) - al * bh) - ah * bl)
return [x, y]
TwoProduct で使用する split というアルゴリズムもエラーフリー変換で，簡単に説明すると $a \in \mathbb{F}$ を $a = a_h + a_l$ のように浮動小数点数の上位ビット/下位ビットに誤差なく分解するイメージです．
split 中にはハイパーパラメータ $t$ が登場しますが，これは扱う浮動小数点数の仮数部のビット長を意味します．Python では浮動小数点数の精度が倍精度なので $t=53$ としています．
また，TwoProduct では split で 8 flop，その他で 9 flop なので合計 25 flop（通常の積の25倍）必要です．
積のエラーフリー変換には割と演算コストがかかってしまうのですが，FMAを用いることで計算コストを抑えられます．
TwoProductFMA
$a$, $b \in \mathbb{F}$ について$\rm{TwoProductFMA}$ を実行すれば,
a * b = x + y
が成り立つ.（TwoProduct の結果と同じ）
# pyfma はデフォルトでインストールされていないので
# 事前に「pip3 install pyfma」しておく必要があります．
from pyfma import fma
def TwoProductFMA(a, b):
x = a * b
y = fma(a, b, -x)
return [x, y]
FMAを用いることで split も必要なくなり，非常にスッキリとしたアルゴリズムになりました．この方法ではわずか 2 flop（通常の積の2倍）で積のエラーフリー変換ができます．そのため，FMA がハードウェアサポートされている環境4であればこちらのアルゴリズムで高速に計算できます．（サポートされていなくても動作はしますが遅いです．）
実験してみた
これまでに紹介したエラーフリー変換アルゴリズムを使って実際に動作確認してみようと思います．
まず TwoSum への入力を $(a,\ b) = (1,\ 3 \mathbb{u})$ として実行します．
この演算で誤差が生じなければ結果は $1 + 3 \mathbb{u}\ $ですが，実際は丸めにより $1 + 4\mathbb{u}$ となります．つまり，TwoSum により誤差 $y = -\mathbb{u} \ $となれば TwoSum が正しく動作していると言えます．
u = 2**-53
[x, y] = TwoSum(1.0, 3*u)
if (x, y) == (1.0 + 4*u, -u):
print("OK")
OK
正しく動作しました．
続いて TwoProduct を入力 $(a,\ b) = (1+2\mathbb{u},\ 1+2\mathbb{u})\ $として実行します．
この演算で誤差が生じなければ結果は $1 + 4 \mathbb{u} + 4 \mathbb{u}^2\ $ですが，こちらも実際は丸めにより $1 + 4\mathbb{u}\ $となります．つまり，TwoProduct により誤差 $y = 4\mathbb{u}^2 \ $となれば TwoProduct が正しく動作していると言えます．
u = 2**-53
[x, y] = TwoProduct(1+2*u, 1+2*u)
if (x, y) == (1.0 + 4*u, 4*u**2):
print("OK")
OK
TwoProduct も正しく動作しました．
GIthub にエラーフリー変換のソースを置いてありますので，実際に手元で動作させてみたい方は利用してください.（そのうち Python 以外の言語でも実装しようと思います．）
次回これらのエラーフリー変換を用いた高精度計算について紹介しようと思います．
参考
-
TJ Dekker, floating-point technique for extending the available precision, Numer. Math. 18, 224-242, 1971 ↩ ↩2
-
明に示したのは Dekker ですが, 1965年に Kahan が総和に関する論文の中で FastTwoSum をそれとなく使っています． ↩
-
D. E. Kuth, The Art of Computer Programming, Volume 2: Seminunierical Algorithms, volume 2, Addison Wesley, Reading, Massachusetts, 1981 ↩
-
Linux 系環境の方は
lscpu | grep fmaでfmaがハードウェアサポートされているか確認できます． ↩