[2405.03588] Effective Quadratic Error Bounds for Floating-Point Algorithms Computing the Hypotenuse Functionという論文を見つけました。この論文では浮動小数点数演算の誤差の上界を自動で計算するソフトが紹介されており、実際にhypot関数を精度良く計算するいくつかの方法について誤差を解析しています。著者のウェブサイトで、そのソフトが公開されています（The BoundRoundingError package – Bruno Salvy's Web Site）。これで遊んでみました。

Mapleのインストール

Mapleは、カナダの会社が開発している数式処理ソフトです。有料ソフトですが、15日間だけつかえる評価版があるので、それをインストールしました（Windows版を選びました）。インストールが終わってMapleを開くと、Start_ja.mwというのが開いているので、ここから新規ドキュメントを選ぶと、ノートブック的な感じのものが出てきます。

使い方

上記ウェブサイトで公開されているBoundRoundingError.zipをダウンロードして、zipを展開します。そうするとBoundRoundingError.mlaというファイルが現れるので、これを先ほど作ったノートブック上にドラッグ＆ドロップします。すると、「このアーカイブには実行可能なコマンドが含まれています。実行しますか？」と出てくるので、「はい (Y)」を選びます。その後、BoundRoundingError:-version;と打ち込み、エンターを押すと、0.3と出てきます。これでこのライブラリ（？）は取り込めたのでしょう。

つづいて、例題を入力していきます。上記ウェブサイトで紹介されている、以下をコピーペーストし、エンターを押します。

Algo:=[
   Input(x=0..2^16,y=0..2^16,_u=0..1/4),
   sx=RN(x^2),
   sy=RN(y^2),
   sigma=RN(sx+sy),
   rho=RN(sqrt(sigma))];

すると、数式っぽいものが描画されます。そこで、BoundRoundingError(Algo);と打つと、gfunというライブラリのバージョンが古くてダメのようです。画面は図1のような感じです（入力した文字が全部イタリックになっていて、これで合っているのかな？と言う感じになったので、画面を紹介します）。

図1: Maple画面のスクリーンショット。gfunのバージョンが古いとのエラーが出た。

そこで、エラーメッセージに記載されている通り、The gfun package – Bruno Salvy's Web Siteからgfun 4.10をダウンロードします。このライブラリを使用する場合は、ウェブサイトに記載の通り、ライブラリの優先順位が重要です。 Mapleにはデフォルトでgfun 3.10が入っているので、図1のようにlibnameを設定してしまうと、（左が優先なので）gfun 4.10が使われません。そこで、libname := "C:\\Users\\lpha\\Downloads\\BoundRoundingError\\BoundRoundingError\\BoundRoundingError.mla", "C:\\Users\\lpha\\Downloads\\gfun\\gfun\\gfun.mla", "C:\\Program Files\\Maple 2025\\lib"のように入力します。あとは同じです。

a + bの誤差上界

以下のように入力したところ、 $u - \frac45 u^2$ と出ました。

図2: Maple画面のスクリーンショット。RN(a+b)の誤差限界がu - 4/5 u²だと出力された

唐突に $\frac45$ が出てきたのが気になりましたが、どうもこれは $u_{\mathrm{max}}$ に依存するようで、 $\frac1{u_{\mathrm{max}}+1}$ になるようです。無限級数をある次数（ここでは二次）までで打ち切る代わりに残りの項の影響をその次数の係数に押し付ける、Type II PSAと呼ばれる方法でしょうか（参考：ベキ級数演算について柏木雅英）。

ここに出力されているものは、ULPで見た誤差ではなく、相対誤差である点に注意します。そもそもここで解析しているものは単位丸め（unit roundoff）が $u_{\mathrm{max}}=\frac14$ であるものです。つまり仮数部がJビット（暗黙の最上位1のこと。いわゆるケチビット）を除いて1ビットの浮動小数点数を解析しています。仮数部が1ビットなので、仮数部は1.0か1.5しかとれません。最大の誤差が発生するのは、1.0+1.5を計算して2.0が得られる場合（絶対誤差は0.5、相対誤差は1/5 (20%)、ULPで見た誤差は0.25ULP）で、たしかに $u - \frac45 u^2$ で正しそうです。

a+b+c+dの誤差上界

計算手順によって、誤差がどのように発生するかが変わることがあります。それを確かめてみます。

a+b+c+dを普通に逐次に計算してしまうと、以下のように $3u$ です。これは、丸めが3回発生するから、と単純に考えられます（入力は全て非負と仮定しているため、桁落ちが発生しないことに注意します）。

このような計算については、Pairwise summationと呼ばれる、二分木状に足し合わせる方式が、誤差の蓄積が少なくてよいとされています。実際、解析してみると $2u -\frac{24}{25}u^2$ だそうです。入力から出力へのパス上に含まれる丸めは、最大2回なので、なんとなくよさそうです。

図4: Maple画面のスクリーンショット。Pairwise summationについて誤差解析している。

入力に負の数が含まれる場合はどうでしょうか？ $-1\le x\le1, 0\le y\le1, 0\le z\le1, 0\le w\le1$ の場合、少し解析に時間がかかりましたが、 $2u -\frac{24}{25}u^2$ だそうです。本当でしょうか？ $x = -1, y = \varepsilon, z = 1, w = 0$ とかで相対誤差1が発生する気がするのですが、私の考察が間違っているのでしょうか。

$-1\le x\le1, -1\le y\le1, 0\le z\le1, 0\le w\le1$ の場合も、 $2u -\frac{24}{25}u^2$ だそうです。 $-1\le x\le1, -1\le y\le1, -1\le z\le1, 0\le w\le1$ の場合、エラーになりました。致命的な桁落ちが発生して評価不能ということでしょうか。

絶対誤差の評価

$0\le y\le \frac1{32}$ に対して、 $1\le x \le \frac{63}{32}$ の範囲でx + yの計算の絶対誤差は、 $u$ に収まるそうです。一方で、 $1\le x \le 2$ の範囲だと、 $\frac{65}{32}u-Ku^2$ になりました。ここで、 $K$ は以下の表のとおりです。

$u_{\mathrm{max}}$	$K$
$2^{-2}$	$\frac{13}{8}$
$2^{-3}$	$\frac{65}{36}$
$2^{-4}$	$\frac{65}{34}$
$2^{-5}$	$\frac{65}{33}$

$K=\frac{65}{32(u_{\mathrm{max}}+1)}$ でしょうか。どういう計算が行われたのかはよくわかりませんが、計算結果が2以下になることが確定する場合と、2を超えることがあり得る場合とで、丸めの絶対誤差が変わる、という点においては間違っていなさそうです。相対誤差だと、 $1\le x \le \frac{63}{32}$ の範囲で $u$ 、 $1\le x \le 2$ の範囲だと $u - \frac45u^2$ と逆転してしまっているので、何か挙動が変です（タイトな評価になっていないだけで間違ってはいなさそうですが）。

入力の比率

論文の31ページ冒頭のRN(x * alpha, 0)のようにやることで、入力xに対して比率がalphaであるyを作ることができる（第二引数の0で、ここの丸めの誤差を考慮しないことを指定）できるようです。

$1\le x\le2$ に対して、 $4\le\alpha\le5$ の時、x / yの誤差は $\frac58 u$ です。 $5\le\alpha\le6$ の時、x / yの誤差は $\frac34 u$ 、 $6\le\alpha\le7$ の時は $\frac78 u$ 、 $7\le\alpha\le8$ の時は $u$ 、 $4\le\alpha\le8$ の時も $u$ 、しかし $4\le\alpha\le9$ の時は $u-2u^2$ だそうです。指数部が特定できる場合と特定できない場合で誤差評価ルーチンが異なるということでしょうか。

Sterbenz lemma

Sterbenz lemmaとは、「比率が二倍以内の同符号の浮動小数点数の引き算は、無誤差である」というやつです。

$0\le x\le 1, 1\le\alpha\le2$ に対してx - yの誤差が、ちゃんと $0$ であると出力されます。他にも、 $1\le x\le 2, 1\le y\le2$ に対しても、x - yの誤差が、ちゃんと $0$ であると出力されます。

しかし、次の例はダメでした。

$2\le x\le3$
$2\le y\le3$
z = y + 1; res = x - z;

Error, (in BoundRoundingError:-SterbenzApplies) cannot determine if this expression is true or false: -2*_u < -3 and 0 < 4+2*_u

だそうです。やっぱりSterbenz lemmaは特殊扱いしているようですね。左側の条件式はyが最小の時Sterbenz lammaが適用可能かを（一回目の丸め誤差込みで）判定する式、右側の条件式はyが最大の時の式、になっています。境界のギリギリを攻めているからダメと言うわけでもないらしく、 $\frac94\le y\le\frac{11}4$ とかでもダメでした。このくらいの条件式の充足可能性を評価することがそんなに難しいとは思えませんが……。

平方根を含む式

この解析ソフトウェアの本領は、平方根を含む式の解析です。例えば、以下のコードを解析してみます。

z = sqrt(y); res = x + z;

$2\le x \le3, 2\le y\le8$ の場合

誤差は $\left(3-\sqrt2\right)u + \frac{14800-10480\sqrt2}{17}u^2$ だそうです（ $u_{\mathrm{max}}=2^{-4}$ の場合）。

$3-\sqrt2$ は1.5858くらいですが、実際にここで発生する誤差は0.5ULP+0.25ULP=0.75ULP止まり（相対誤差は1.5uが上界）のはずで、ちょっと過大評価ですがそれほど的外れではないでしょう。

$2\le x \le3, 2\le y\le7$ の場合

誤差は $\left(\frac{10}3-\frac{2\sqrt7}3\right)u + \frac{30304+6144\sqrt{14}-21504\sqrt2-26016\sqrt7}{51}u^2$ だそうです（ $u_{\mathrm{max}}=2^{-4}$ の場合）。

$\frac{10}3-\frac{2\sqrt7}3$ は1.5695くらいで、やはり過大評価ですが、過大評価の程度は大きくありません。

$2\le x \le3, 2\le y\le4$ の場合

Error, (in BoundRoundingError:-ratapprox) Digits cannot exceed 38654705646とエラーが出ました。値を小さくしているのにこんなエラーが出るのは謎です。ぴったりすぎて等価性判定が無限ループしたとかでしょうか。

$2\le x \le3, 2\le y\le3$ の場合

誤差は $\left(2-\frac{\sqrt2}2\right)u-\frac{8\sqrt2}{17}u^2$ だそうです（ $u_{\mathrm{max}}=2^{-4}$ の場合）。

$2-\frac{\sqrt2}2$ は1.2929くらいです。 ULPで見た誤差の上界は0.5ULP+0.25ULP=0.75ULPですが、これに近い誤差が発生するのはresが4未満の時であり、相対誤差が1.5u付近になることはありません。相対誤差が最も大きくなるのはresが4付近の時であり、その上界は0.5ULP+0.125ULP=0.625ULPで1.25uです。 1.5uよりちゃんと小さくなっている点は、ちゃんと評価されているのだなぁという感じです。それでも微妙な過大評価が残っています。

$2\le x \le3, 2\le y\le\frac94$ の場合

応答が返ってきませんでした。

無誤差変換

FastTwoSumやTwoProdなどは、そのままではまともに取り扱うことができません。しかしながら、28ページの上の方の数式で出ているように、sigmal=RN(sxh+syh-sigmah,0)のように無誤差を指定する数式で誤差項を強制的に導入することで、解析が可能になります。

例題：Sloppy加算の絶対誤差評価

Sloppy加算とは、以下の手順でdouble-doubleの和を計算する方法です。

DD SloppyAdd( DD x, DD y ) {
    auto [hi, lo] = TwoSum( x.hi, y.hi );
    return FastTwoSum( hi, lo + x.lo + y.lo );
}

これは相対誤差の意味でまずい計算方法であり、最悪の場合相対誤差が1になってしまう（真の結果が非0なのに、0を返してしまう）ことがあることが知られています。

しかしながら、そもそもそんな桁落ちが発生するような計算をしている時点で相対誤差が大きくなるのは当たり前であり、Sloppy加算を責めるのは筋違いに感じます。ここは絶対誤差が重要と思われますが、それを評価した論文は見当たりませんでした。なのでずっとこれを評価したいと考えていたのですが、私の手には負えなくて困っていたところでした。

まず、次のように入力を作りました。

Algo := [
  Input(xh=-1..1, yh=-1..1, alpha=-2^(-53)..2^(-53), beta=-2^(-53)..2^(-53), _u=0..2^(-53)),
  xl = RN(alpha*xh, 0),
  yl = RN(beta*yh, 0),
  hi = RN(xh + yh),
  lo = RN(xh + yh - hi, 0),
  tmp1 = RN(xl + yl),
  tmp2 = RN(lo + tmp1),
  res = RN(xh + yh + xl + yl - hi - tmp2, 0)
]

まず、xlとylをalpha倍やbeta倍で作り出します。正規化されているdouble-doubleは、下位パートの絶対値が上位パートの絶対値の $2^{-53}$ 倍までに収まるので、それを使います。次に、TwoSumの下位パート計算部分も無誤差項導入で作り出します。最後に、評価したい誤差項を、無誤差項導入で作り出します。これで本当にいいんでしょうか……！？

そうすると、以下の答えが返ってきます（verboseレベルを3に設定した場合）。

BoundRoundingError:
            "step-by-step analysis in 287.296 sec."

BoundLinearTerm: bounding the linear part of the absolute error
Optimize: computing upper bound for (-alpha*xh-beta*yh)*_epsilon[tmp1]+(-alpha*xh-beta*yh)*_epsilon[tmp2]
Optimize: bound: 1/2251799813685248
BoundRoundingError:
      "linear part: 1/2251799813685248*_u, in 2.844 sec."

BoundQuadraticTerm: bounding the quadratic part of the absolute error
Error, (in BoundRoundingError:-BoundQuadraticTerm) relative error of 0 undefined

ここにでてくる2251799813685248は $2^{51}$ であり、要するに入力の絶対値が1以下の数同士のSloppy加算の最大絶対誤差が $2^{-104}$ であることを示しているようです。やはり、絶対誤差で見れば、double-doubleとして許容できる程度の精度を持つようです。

解析している数式を眺めてて気づいたんですが、Sloppy加算のアルゴリズム中で、回収されない丸め誤差が生じるのは2回だけなので、そんなに証明は難しくなかったのかもしれません。

正規化されていない場合

Sloppy加算で、入力のdouble-doubleが正規化されていない場合はどうでしょうか？例えば、 $|\alpha|, |\beta| \le 2^{-50}$ に設定してみると以下のようになりました。

BoundRoundingError:
            "step-by-step analysis in 180.609 sec."

BoundLinearTerm: bounding the linear part of the absolute error
Optimize: computing upper bound for (-alpha*xh-beta*yh)*_epsilon[tmp1]+(-alpha*xh-beta*yh)*_epsilon[tmp2]
Optimize: bound: 1/281474976710656
BoundRoundingError:
       "linear part: 1/281474976710656*_u, in 3.125 sec."

BoundQuadraticTerm: bounding the quadratic part of the absolute error
Error, (in BoundRoundingError:-BoundQuadraticTerm) relative error of 0 undefined

ここにでてくる281474976710656は $2^{48}$ であり、入力の正規化されていない度合いが3bit増えたことに応じて出力の誤差も3bit悪化するということのようです。

乗算

Sloppy乗算の絶対誤差

加算をやったので、せっかくなので乗算も見ていきましょう。試してみるのは、Tight and rigorous error bounds for basic building blocks of double-word arithmeticという論文*1のAlgorithm 11です。 $u^2$ 以下の大きさである $x_{lo}y_{lo}$ 項を省略する方式のdouble-double乗算で、最大誤差は $6u^2$ で抑えられると証明されています（実験的には $5u^2$ くらいの誤差しか生じないらしいですが……）。

Algo := [
  Input(xh=-1..1, yh=-1..1, alpha=-2^(-53)..2^(-53), beta=-2^(-53)..2^(-53), _u=0..2^(-53)),
   xl = RN(alpha*xh, 0),
   yl = RN(beta*yh, 0),
   hi = RN(xh*yh),
   tmp1 = RN(xh*yh - hi, 0),
   tmp2 = RN(xh*yl),
   tmp3 = RN(xl*yh + tmp2),
   lo = RN(tmp3 + tmp1),
   err = RN((xh + xl)*(yh + yl) - hi - lo, 0)
]

step-by-step analysisに276.984秒、一次の項の導出に9.969秒、二次の項の導出に63.25秒かかって、以下の結果を得ました。

$\frac5{2^{53}}u + \frac{2^{159}-2^{55}-2^{53}-5}{\left(2^{53}+1\right)^3}u^2$

ちゃんと二次の項まで見ると、たしかに $6u^2$ に収まっていることがわかります。ただ、これは相対誤差ではなく絶対誤差なので、当該論文で示されている結果より弱い結果になっています（相対誤差で見ると $24u^2$ を示したことになっています）。

Sloppy乗算の相対誤差

次のルーチンを解析することで、相対誤差を評価します。

Algo := [
  Input(xh=-1..1, yh=-1..1, alpha=-2^(-53)..2^(-53), beta=-2^(-53)..2^(-53), _u=0..2^(-53)),
  xl = RN(alpha*xh, 0),
  yl = RN(beta*yh, 0),
  hi = RN(xh*yh),
  tmp1 = RN(xh*yh - hi, 0),
  tmp2 = RN(xh*yl),
  tmp3 = RN(xl*yh + tmp2),
  lo = RN(tmp1 + tmp3),
  err = RN(1 - (hi + lo)/((xh + xl)*(yh + yl)), 0)
]

step-by-step analysisに156.203秒、一次の項の導出に1.422秒、二次の項の導出に18.015秒かかって、以下の結果を得ました。

$\frac{5\times2^{53}}{\left(2^{53}-1\right)^2}u+\frac{2^{53}}{\left(2^{53}-1\right)^2}u^2$

ほぼ $6u^2$ ですが、わずかに大きいです。これはおそらく、加算・乗算の丸め誤差を $u$ と評価するか、 $\frac u{1+u}$ と正確に評価するか、という部分の違いだと思います。

入力範囲を分割するとよりタイトな結果が得られることがある（論文24ページ）らしいので、いくつか分割を試してみましたが微妙でした。積が繰り上がるか繰り上がらないが重要な分岐点になると思ったのですが、そうではなかったようです。

$1\le x\le\frac43, 1\le y\le\frac32$ としても、上記と同じ（ $6u^2$ よりわずかに大きい）
$\frac32\le x\le2, \frac32\le y\le2$ とした場合、 $\left(5+\frac89\right)u^2$ よりわずかに大きいくらい。多少厳しくなった

このあたりのルーチンを解析するだけでも5分とかかかってなかなかいろいろ試せなくてもどかしいです。あと、同じノートブックで続けて解析していると妙に遅くなるような気がします。メモ化とかをしていてその照合に時間がかかっているのでしょうか？

xhやyhの範囲を1～2に固定したら、αやβを定義せずに直にxlやylを定義したほうがよりタイトな評価になりそうなので、また今度やってみます。

まとめ

SalvyのBoundRoundingErrorで遊びました。全ての浮動小数点数プログラムを自動的にタイトに誤差評価できる――というところには程遠いですが、遊んでみるといろいろ非自明な結果が得られて、数式処理ソフトウェアの威力を体感できました。また、絶対値が1以下のdouble-double同士のSloppy加算の最大絶対誤差は $2^{-104}$ であるという結果を得ました。

ぽちぽち入力すると何かの応答が返ってくるという意味で、LLMに通じるものがある気がします。 AIの言いなりになるのではなく、アイデアの種を生み出すくらいの使い方にとどめて、自分の頭を使って考えたいですね。

*1:Mioara Maria Joldes, Jean-Michel Muller, Valentina Popescu: "Tight and rigorous error bounds for basic building blocks of double-word arithmetic", ACM Transactions on Mathematical Software, 44(2) pp. 1-27 (2017).