小清水 (@curekoshimizu) です。

仕事が忙しくなかなか更新する余力がありませんでしたが一年半振りにブログを更新します！

エラーフリー変換とは？

浮動小数点数の計算というのは、このブログで何度も紹介していますが、丸め誤差を伴います。

簡単な $x+y$ という計算であっても $x + y$ を行ったあとに浮動小数点に変換されるために丸め誤差が発生する可能性があるのです。その他にも

$xy$
$xy + z$
$xy + zw$
$\sum_{n=1}^N a_n$
$\sum_{n=1}^N a_n b_n$

などなど各種計算も同様に計算過程で丸め誤差をともないます。

つまり

(求めたい計算式) = (愚直に浮動小数点体系で計算して得た結果) + (丸めによる誤差)

例えば

$x = 0.300000011920928955078125 = 1.20000004768371581677722\times 2^{-2}$ $y = 0.20000000298023223876953125 = 1.6000000238418585033165\times 2^{-3}$

において $x + y$ を float32 で計算すると $0.5 = 1.0 \times 2^{-1}$

となり、誤差 $0.00000001490116119384765625$ が生じます。

この誤差項も含めて結果を得たいというケースがあるのです。

これを行うのがエラーフリー変換 (Error Free Transformation) と呼ばれる計算です。

こんなこと本当にしたいのでしょうか？

一つの例としてみなさんがお世話になっている glibc でも使われているのだというのを見てみましょう。

glibc でエラーフリー変換が使われている例

たとえば glibc ではソフトウェアFMA を計算するために次のような処理をしています。

https://github.com/lattera/glibc/blob/master/sysdeps/ieee754/dbl-64/s_fma.c の LN188 - LN205 あたり。

  /* Multiplication m1 + m2 = x * y using Dekker's algorithm.  */
  #define C ((1 << (DBL_MANT_DIG + 1) / 2) + 1)
  double x1 = x * C;
  double y1 = y * C;
  double m1 = x * y;
  x1 = (x - x1) + x1;
  y1 = (y - y1) + y1;
  double x2 = x - x1;
  double y2 = y - y1;
  double m2 = (((x1 * y1 - m1) + x1 * y2) + x2 * y1) + x2 * y2;

  /* Addition a1 + a2 = z + m1 using Knuth's algorithm.  */
  double a1 = z + m1;
  double t1 = a1 - z;
  double t2 = a1 - t1;
  t1 = m1 - t1;
  t2 = z - t2;
  double a2 = t1 + t2;

github.com

ここで登場するコメントにも書かれている Dekker の乗算アルゴリズムやKnuth の加算アルゴリズム

はエラーフリー変換の一つです。

  /* Multiplication m1 + m2 = x * y using Dekker's algorithm.  */
  #define C ((1 << (DBL_MANT_DIG + 1) / 2) + 1)
  double x1 = x * C;
  double y1 = y * C;
  double m1 = x * y;
  x1 = (x - x1) + x1;
  y1 = (y - y1) + y1;
  double x2 = x - x1;
  double y2 = y - y1;
  double m2 = (((x1 * y1 - m1) + x1 * y2) + x2 * y1) + x2 * y2;

例えばこの箇所は x*y を計算する m1 = x * y という計算によってどの程度の誤差が生じたのかを知るために m2 を厳密に計算しています。

これにより

m1 + m2 = x * y

となることを数学的に保証するというものです。

こちらの加算も同様で

  /* Addition a1 + a2 = z + m1 using Knuth's algorithm.  */
  double a1 = z + m1;
  double t1 = a1 - z;
  double t2 = a1 - t1;
  t1 = m1 - t1;
  t2 = z - t2;
  double a2 = t1 + t2;

a1 = z1 + m1 で求めたい加算結果を得て、それによって生じた誤差を a2 で手に入れています。

乗算に関するエラーフリー変換と FMA 命令

上の Dekker の乗算アルゴリズムは FMA 命令がないアーキテクチャのためのアルゴリズムとして知られています。

もし FMA命令 (xy+zを計算する命令) を持っている場合には

  /* Multiplication m1 + m2 = x * y using FMA instruction.  */
  double m1 = x * y;
  double m2 = x * y - m1;

これでエラーフリー変換できてしまいます。そのため、現代のアーキテクチャで、この Dekker の乗算アルゴリズムを実際に使う機会は少ないかもしれません。

FMA についてはこちらで詳しく記述しました：

math-koshimizu.hatenablog.jp

加算に関するエラーフリー変換

一方で加算はそうはいきません。FMA命令では解決できないのです。

歴史 -- FastTwoSum(Fast2Sum)・TwoSum アルゴリズム

1971 年に Dekker に紹介された FastTwoSum アルゴリズムというアルゴリズムがあります。(T.J. Dekker: A floating-point technique for extending the available precision)

これが加算に関するエラーフリー変換として最も古いアルゴリズムです。とはいえ、

1965 年に Kahan が総和演算に関する Compensated summation アルゴリズムに関する論文の中でこのアルゴリズムを暗に利用しているのも事実ですので、最古といっていいか悩ましいところはあります…。(https://convexoptimization.com/TOOLS/Kahan.pdf)

この Dekker のアリゴリズムを FastTwoSum と名付けたのは Shewchuk の 1997年の論文です：

f:id:curekoshimizu:20190421191324p:plain — Shewchuk : Adaptive precision floating=point arithmetic and fast robust geometric predicates

ただし、この FastTwoSum は上記通り、条件がついており、引数a, bについて大小比較をする必要があります。

このため、現在のアーキテクチャの高速化技法等を考えると、条件分岐が発生するこのアルゴリズムよりも、

KunuthとMøllerによる条件分岐の発生しない加算エラーフリー変換である、 TwoSum アルゴリズムのほうが一般的に利用されています。

上の glibc で書かれていた

  /* Addition a1 + a2 = z + m1 using Knuth's algorithm.  */
  double a1 = z + m1;
  double t1 = a1 - z;
  double t2 = a1 - t1;
  t1 = m1 - t1;
  t2 = z - t2;
  double a2 = t1 + t2;

この式はこのアルゴリズムです。

命名は先ほど登場した shewchuk の論文です。

論文中に FastTwoSum(a, b), Fast2Sum(a, b), TwoSum(a, b), 2Sum(a, b) などと登場したらこれらのアルゴリズムのことなんだなと思うと良いと思います。

FastTwoSum アルゴリズムの紹介と証明

今回のブログではこの古典的な加算エラーフリー変換アルゴリズムである Dekker の FastTwoSum アルゴリズムについて紹介と証明を与えてみようと思います。

FastTwoSum アルゴリズム

次の計算により $|a| \geq |b|$ の場合に $a+b = s + t = \mathrm{RN}(a+b) + b$ なる $s, t$ が得られます。ただしここでは2進数浮動小数点環境とする：

Dekker FastTwoSum

$s = \mathrm{RN}(a + b)$
$z = \mathrm{RN}(s - a)$
$t = \mathrm{RN}(b - z)$

ただし $\mathrm{RN}$ は最近接偶数方向丸めを表します。この最近接偶数方向丸めについては次の記事を参照ください：

math-koshimizu.hatenablog.jp

FastTwoSum アルゴリズム証明

ここで

$a = M_a 2^{e_a}$
$b = M_b 2^{e_b}$
$s = M_s 2^{e_s} = \mathrm{RN}(a+b)$

とおく、ただし $M_a, M_b, M_s, e_a, e_b, e_s \in \mathbb{Z}$ かつ $|M_a|, |M_b|, |M_c|$ < $2^p$ , $e_{min} \leq e_a, e_b, e_s \leq e_{max}$ とする。

ここで、次のブログの Prop. (整数を用いた $\beta$ 進 $p$ 桁の浮動小数点数の表示) を使用した。この証明ではこの性質を頻繁につかっている。

math-koshimizu.hatenablog.jp

また、

$|a + b| \leq |a| + |b| \leq 2|a| = 2 |M_a| 2^{e_a} = |M_a| 2^{e_a+1}$

ここから、 $M_s$ と $e_s$ のとり方の自由度はあるものの、 $e_s \leq e_a + 1$ と制限した上で $M_s$ を定めても問題ないことがわかる。

ここで次の3つにわけて考える

Case1. $e_s > e_a$ (すなわち $e_s = e_a + 1$ )
Case2. $e_s \leq e_a$ かつ $e_s > e_b$
Case3. $e_s \leq e_a$ かつ $e_s \leq e_b$

このとき、

$s-a$ が浮動小数点として正確に表現できることを示し、 $s-a = (整数)\times 2^{e_b-p+1}$ と表せることを示す。

Case1. $e_s = e_a + 1$ のとき

このとき

$a + b = M_a {2}^{e_a} + M_b {2}^{e_b} = \left( \frac{M_a}{2} + \frac{M_b}{2^{e_a - e_b+1}} \right) \times 2^{e_a+1} = \left( \frac{M_a}{2} + \frac{M_b}{2^{e_a - e_b+1}} \right) \times 2^{e_s}$

となり、

$\left| \frac{M_a}{2} + \frac{M_b}{2^{e_a - e_b+1}} \right| \leq \frac{|M_a|}{2} + \frac{|M_b|}{2} \leq | M_a | < 2^{p}$

と、 $\left( \frac{M_a}{2} + \frac{M_b}{2^{e_a - e_b+1}} \right)$ は $2^p$ で抑えられる実数として表現される。

$a +b$ 最近接偶数方向丸めが $s$ であるから、先の項を整数化することになる。よって

$\left| M_s - \left( \frac{M_a}{2} + \frac{M_b}{2^{e_a - e_b+1}} \right)\right| \leq \frac{1}{2}$

と評価でき、

これより次の式が得られる：

$| 2M_s - M_a | \leq 2 \left| M_s - \left( \frac{M_a}{2} + \frac{M_b}{2^{e_a - e_b+1}} \right)\right| + 2\left| \frac{M_b}{2^{e_a - e_b+1}} \right| \leq 1 + | M_b | \leq 2^p$

ゆえに

$s - a = M_s 2^{e_a+1} - M_a 2^{e_a} = (2 M_s - M_a ) 2^{e_a}$

から、この $2 M_s - M_a$ は整数であって、先の評価から $2^p$ 未満の場合、 $s-a$ は浮動小数点表示される。もちろん $2^{p+1}$ に一致する場合も、指数部を $2^{e_a+1}$ に修正しても、整数部は丸めによる誤差は生じないので、浮動小数点表示されていることがわかる。

そしてこの表示から $e_b \leq e_a$ より $s-a = (整数)\times 2^{e_b}$ と表されることもわかる。

Case2. $e_s \leq e_a$ かつ $e_s > e_b$ のとき

ここで

$a + b = M_a 2^{e_a} + M_b 2^{e_b} = \left( 2^{e_a - e_s} M_a + \frac{M_b}{ 2^{e_s-e_b}} \right) 2^{e_s}$

最近接遇数方向の定義より

$\left|M_s- \left( 2^{e_a - e_s} M_a + \frac{M_b}{ 2^{e_s-e_b}} \right) \right | \leq \frac{1}{2}$

ここから、

$|s-a| = | M_s 2^{e_s} - M_a2^{e_a}| \leq \left| M_s - \left( 2^{e_a - e_s} M_a + \frac{M_b}{ 2^{e_s-e_b}} \right) \right| 2^{e_s} + \frac{|M_b|}{2^{e_s-e_b}} 2^{e_s} \leq \left( \frac{1}{2} + \frac{|M_b|}{2^{e_s-e_b}} \right) 2^{e_s}$

が得られる。 $s-a = \left( M_s 2^{e_s} - M_a 2^{e_a-e_s} \right) 2^{e_s} =: K 2^{e_s}$ と整数部を $K$ と表すと