以下が流れてきて気になったので考えてみました。

係数が（その文脈で考えている）浮動小数点数として表せる二次方程式が与えられたとき、その実数解の correctly-rounded な値を得る方法はあるか？
— えびちゃん🍑🍝🦃 (@rsk0315_h4x) 2025年8月25日

解きたい方程式は $g(x) = ax^2 + bx+c = 0$ であるとします。

私の結論は、「単精度浮動小数点数係数二次方程式の実数解の correctly-rounded な値を得るためには、倍精度浮動小数点数の積和演算（fma）があれば十分で、倍倍精度演算や多倍長演算は不要」です。ただ、細部まで証明したわけではないので、穴があるかもしれません。以下、まず何が難しいかを説明し、その後手順の概要を説明し、最後に手順の詳細を説明します。

難しいところ

係数は単精度浮動小数点数ですが、a*x*x+b*x+cとか(-b+sqrt(b*b-4*a*c))/(2*a)とか計算してしまうと丸め誤差が何回も生じる点が厄介です。二次方程式を解く手法は大きく分けて、二分法*1、ニュートン法、解の公式、などがありますが、いずれも丸め誤差の影響を受けた場合に正しく答えが求まるかは明らかではありません（一般には正しく求まりません）。

二分法では、 $ax^2 + bx+c$ の符号を一回でも間違えてしまうと、全くおかしなところに収束してしまいます。 $ax^2 + bx+c$ の計算は桁落ちになることがあり、桁落ちすれば最終結果の符号が変わってしまうようなことは全然あり得る（例えばRump(るんぷ)の例題とか……）ので、単純な方法ではうまくいきそうにありません。

ニュートン法では、 $x_{i+1} \leftarrow x_i + \delta$ のような更新が出てきますが、ここが丸めに弱いです。 $\delta$ が $x_i$ の0.5ULPよりも小さくて積み残しになってしまって無限ループになってしまうミスをよく見かけます*2。

解の公式は言わずもがな、丸め誤差の影響を受けます。桁落ちとその回避策の例示として、二次方程式の解の公式の愚直適用と分母の有理化（1/xの二次方程式とみなして解く）が紹介されることがよくあります。なので分母を有理化すればいいのでは？と思うかもしれませんが、そもそも分母の有理化をしたところで複数回の計算の丸め誤差が積み重なる点は変わっていないので、“大きな”誤差は生じなくなるでしょうが、correctly-rounded な値を得ることは困難です。

大まかな方針

手順

以下の手順によります。ただし、指数部範囲が十分広くとれるとの理想化を行いました。これは、オーバーフローやアンダーフローが起こるような入力が来ないものと仮定したという意味に解釈して問題ありません。また、一般性を失わずに $a > 0$ としました。

判別式Dを計算し、負なら実数解なしである。NaNを返すなど、適宜仕様通りの動作をすればよい。
判別式Dがちょうど0であれば、重解を一つもち、その真の値は $x_T = -\frac b {2a}$ である。 $x_T$ の最近接浮動小数点数 $f_T$ を返せばよい。
判別式Dが正であれば、実数解を二つ持つ。
1. $g(f_T) \gt 0$ であれば、二つある実数解のうち少なくとも一つについて、最近接の浮動小数点数が $f_T$ であるからこれを返せばよい。
2. $g(f_T) = 0$ であれば、 $f_T$ が実数解なので、これを返せばよい。
3. $g(f_T) \lt 0$ であれば、 $f(\mathrm{LARGE}) > 0$ と組み合わせて二分法を実行し、実数解を挟む1ULP差の二つの浮動小数点数 $f_L, f_R$ を手に入れることができる。 $g(\frac{f_L+f_R}2)$ の符号に応じて $f_L$ か $f_R$ を返せばよい。

最後の部分のアイデアは、以下によります。

x の次の浮動小数点数を x+e として、例えば f(x) < 0 < f(x+e) であれば、 f(x+(e/2)) の正負を調べれば分かりそうです。
2 解がとても近い場合などコーナーケースがあるのでそのへんはうまく処理していただいて・・
— きり (@kiri8128) 2025年8月25日

その他の分岐（特に3-1と3-2）は、コーナーケース対策です。つまり、二分法を実行するためには少なくとも一つ、 $g(f) \lt 0$ を満たす浮動小数点数 $f$ を見つける必要がありますが、これができないケースです。 $\forall f \in \mathbb{F}. f \ne f_T \rightarrow g(f_T) \lt g(f)$ が成り立つので（記事の最後に証明を添付しました）、 $f_T$ だけ確認すれば十分なのが便利です。

ところで、元の問題に対する正しい実装は、実数解が二つあるとき、それぞれに対する最近接浮動小数点数を出力するというものだと思います。一方でこの手順は、実数解が二つあるときは、好きな方の実数解一つに対する最近接浮動小数点数を出力するものになっています。もう片方を出力する方法については、アイデアの段階ですが、以下のような感じでしょうか。

3-1と3-2のケースは、を挟む1ULP差の二つの浮動小数点数（片方はと一致、もう片方をとすれば）を使って3-3と同様にすれば良さそう
- $g(f_T)\lt0$ ではないが、 $f_L$ と $f_R$ の間に $g(x)\lt0$ なる $x$ があることは分かっている（ $x_T$ がそれ）ので、二分法と同じようにできるはず
3-3のケースは、 $f(-\mathrm{LARGE}) \gt 0$ を使えばもう片方の解を手に入れることができそう

示すべき性質

以下の事実は証明しておく必要があるでしょうか。

判別式Dが正で $g(f_T) \gt 0$ の時に、二つある実数解のうち少なくとも一つについて、最近接の浮動小数点数が $f_T$ である

これを証明します。

まず、判別式Dが正なので、 $g(x_T) \lt 0$ が成り立ちます。さらに $g(f_T) \gt 0$ なので、この間に実数解 $x_0$ が存在します。ここで、 $|f_T - x_0| \lt |f_T - x_T| \lt |f_S - x_T| \lt |f_S - x_0|$ なので、 $x_0$ は $f_S$ よりも $f_T$ に近いことになります。 $x_0$ を挟む1ULP差の二つの浮動小数点数は $f_S$ と $f_T$ で、それ以外の浮動小数点数が最近接になることはあり得ないので、 $f_T$ が最近接ということになります。

数式で書くと無意味に冗長ですが、以下の図のようになっているというのが直観的な証明です。

f_T        x_0     x_T                  f_S
_|__________v_______v____________________|_
 |                                       |

示すべき手順

先ほどの手順で、実行可能性が明らかでない部分は以下の通りです。

判別式の符号を正しく判定できるのか。
$x_T$ の最近接浮動小数点数 $f_T$ を求めることができるのか。
浮動小数点数 $f$ に対して、 $g(f)$ の符号を求める（-か0か+かを判定する）ことができるのか。
浮動小数点数とは限らない $m=\frac{f_L+f_R}2$ に対して、 $g(m)$ の符号を求めることができるのか。

このうち、「 $x_T$ の最近接浮動小数点数 $f_T$ を求めることができるのか」については、-bと2*aは浮動小数点数で表せることから、単に浮動小数点数除算を行うことで実現できることがわかります。残り3つについては、実現方法を示していきます。

判別式の符号を正しく判定できるのか。

$b^2$ と $4ac$ はどちらも48ビット整数×二冪の形をしています。よって、(double)b * bや4.0 * a * cのように計算すれば、無誤差です。つまり、doubleの乗算とdoubleの比較を使えば、判別式の符号を正しく判定することは容易です。

なお、この部分はfloatのTwoProdと辞書順比較でもできますね。

浮動小数点数 $f$ に対して、 $g(f)$ の符号を求める（-か0か+かを判定する）ことができるのか。

浮動小数点数とは限らない $m=\frac{f_L+f_R}2$ に対して、 $g(m)$ の符号を求めることができるのか。

結局、仮数部が25ビットの $f$ に対して $g(f)$ の符号を求められれば十分です。倍精度浮動小数点数の積和演算（fma）を用いて $g(f)$ を評価するにあたって、なるべく丸めの回数を少なくしたいです。そのような方法としては、次の三通りが考えられます。ここで、 $f^2$ は50ビット整数×二冪なので、doubleで計算すれば誤差を生じないことに注意します。

fma( fma( a, f, b ), f, c )……手順A
fma( b, x, fma( a, (double)f*f, c ) )……手順B
fma( a, (double)f*f, fma( b, f, c ) )……手順C

この三つの計算手順はどれもfmaを二回含んでおり、一般にはそのどちらもで誤差が生じ得ます。したがって、この手順で計算した結果得られた符号は、真の符号と一致しない可能性があります。実際に問題となるのは、内側のfmaで誤差が生じた後、外側のfmaが桁落ち（絶対値が近い数同士の引き算）になるパターンです。なぜなら、以下の二つが成り立つからです。

主張1 内側のfmaで誤差が生じなければ、符号を簡単に判定できる。
理由の説明 一般に丸め誤差は符号を変えないから、特に外側のfmaで生じる丸め誤差は、それが生じたとしても符号判定において問題とならない。
※アンダーフローが発生しないと仮定しているためこのように結論付けられている点に注意。一般には有限の数が0に丸められることがあり得る。

主張2 外側のfmaで桁落ちにならない（足し算か、絶対値が大きく異なる数の引き算）ならば、符号を簡単に判定できる。
理由の説明 いずれの場合も、計算結果の符号が変わるほど大きな丸め誤差が生じることはあり得ない。
※そのように「桁落ち」を定義したというだけの話。具体的に安全かは都度確認します。

なので、以下では計算手順で「内側のfmaで誤差が生じた後、外側のfmaが桁落ち」が発生しないことを示します。

場合分けの準備

以下、一般性を失わずに $a, f \in [1, 2)$ とします。というのも、まず $f$ について、 $af^2 + bf + c = (a/4)(2f)^2 + (b/2)(2f) + c$ であり、 $a, b$ が浮動小数点数であれば $a/4, b/2$ も浮動小数点数であることを活用すれば、 $f$ の指数部を自由に変更することができます。また、 $a$ についても、全体を二冪でスケーリングすれば、 $a$ の指数部は自由に変更することができます。さらに、 $f$ の符号反転は $b$ の符号反転で、 $a$ の符号反転は $b$ と $c$ の符号反転で、それぞれできます。以上より、わかりやすく $a, f \in [1, 2)$ に取ることができます。

さて、 $af^2$ 、 $bx$ 、 $c$ の符号は、計算せずともわかります。また、 $bf+c$ 、 $af^2+c$ 、 $af+b$ の符号も、計算してみればわかります（丸め誤差は符号判定に影響ないことに注意します）。このことを利用して、以下のように分類してみます。

$b=0$ なら、手順Bを使えば、外側のfmaは値を変えないので判定が可能です。
$c=0$ なら、手順Aを使えば、外側のfmaは値を変えないので判定が可能です。
$a \gt 0, b \gt 0, c \gt 0$ なら、明らかに $g(f) = af^2 + bf + c \gt 0$ です。
の場合
- $bf+c \ge 0$ なら、明らかに $g(f) = af^2 + bf + c \gt 0$ です。
- $af^2+c \ge 0$ なら、明らかに $g(f) = af^2 + bf + c \gt 0$ です。
- 残る $b \gt 0, c \lt 0, bf+c \lt 0, af^2+c \lt 0$ は難しいパターンです（場合1）。
の場合
- $bf+c \ge 0$ なら、明らかに $g(f) = af^2 + bf + c \gt 0$ です。
- $af+b \ge 0$ なら、明らかに $g(f) = af^2 + bf + c \gt 0$ です。
- 残る $b \lt 0, c \gt 0, bf+c \lt 0, af+b \lt 0$ は難しいパターンです（場合2）。
の場合
- $af^2+c\le0$ なら、明らかに $g(f) = af^2 + bf + c \lt 0$ です。
- $af+b\le0$ なら、明らかに $g(f) = af^2 + bf + c \lt 0$ です。
- 残る $b \lt 0, c \lt 0, af^2+c \gt 0, af+b \gt 0$ は難しいパターンです（場合3）。

以下、この三種類の難しいパターンに対して、 $g(f) = af^2 + bf + c$ の符号を判定する手順を説明します。この中では場合1が最も難しいです。そこで、まずは場合2と場合3から示し、最後に場合1を示します。

場合2： $b \lt 0, c \gt 0, bf+c \lt 0, af+b \lt 0$

まずfma(a,f,b)を計算してみます。これが無誤差なら、手順Aで計算すれば符号判定ができます。よって、以下、fma(a,f,b)が誤差を生じる場合を考えます。次にfma(b,f,c)を計算してみます。これが無誤差なら、手順Cで計算すれば符号判定ができます。よって、以下、fma(b,f,c)も誤差を生じる場合を考えます。

さて、fma(a,f,b)が誤差を生じるというのはどういう時でしょうか。 $af$ は49ビット整数×二冪なので、 $b$ がこれを5ビット以上はみ出しているということです。正規化した都合上、 $af$ は1以上です。 $af+b \lt 0$ なので、 $b \lt -1$ であることが従い、上にはみ出していることが確定します。具体的には、 $b \le -32$ となります。 $f \ge 1$ なので $bf \le-32$ です。

また、fma(b,f,c)が誤差を生じるということから、 $bf$ と $c$ は絶対値が少なくとも2倍以上違うということがわかります。 $bf$ と $c$ のうち絶対値の大きい方（ $L$ とする）の絶対値は少なくとも32あり、絶対値が小さい方（ $S$ とする）の絶対値との差は少なくとも16あります。 $af^2$ は高々8なので、 $L$ の符号は、 $S + af^2$ では覆せません。よって、 $L$ の符号が $g(f) = af^2 + bf + c$ の符号です。

以上により、少なくとも倍倍精度演算や多倍長演算に頼らなくても、符号を判定できることがわかりました。

ここまでの手順を確認すると、

$b \gt -32$ ならば手順Aで計算すれば、一回目のfmaが無誤差なので正しく符号判定できる
$b \le -32$ ならば手順Cで計算すれば、一回目のfmaは誤差を生じるかもしれないけれど符号判定には圧倒的に余裕があるので正しく符号判定できる

のような単純な分岐でいけそうな気がします。

場合3： $b \lt 0, c \lt 0, af^2+c \gt 0, af+b \gt 0$

まずfma(a,f,b)を計算してみます。これが無誤差なら、手順Aで計算すれば符号判定ができます。よって、以下、fma(a,f,b)が誤差を生じる場合を考えます。次にfma(a,(double)f*f,c)を計算してみます。これが無誤差なら、手順Bで計算すれば符号判定ができます。よって、以下、fma(a,(double)f*f,c)が誤差を生じる場合を考えます。

まず、fma(a,f,b)が誤差を生じているため、やはり49ビット整数×2冪である $af$ に対して $b$ が5ビット以上はみ出しています。正規化した都合上 $af$ は4未満である一方、 $af+b \gt 0$ なので、 $b \gt -4$ であることが従い、下にはみ出していることが確定します。具体的には、 $b \gt -2^{-28}$ となります。 $1 \le f \lt 2$ なので、 $bf \gt -2^{-27}$ となります。

また、fma(a,(double)f*f,c)が誤差を生じているのはどういう時でしょうか。正規化した都合上 $af^2$ は $2^{-71}$ の倍数です。この時、

$c$ が $2^{-71}$ の倍数にならないなら、 $|c|\lt2^{-48}$ であり、 $af^2+c \gt 1 - 2^{-48} \ge 2^{-18}$
$c$ が $2^{-71}$ の倍数なら $af^2+c$ も $2^{-71}$ の倍数で、誤差を生じている（53ビットに入らなかった）ことから少なくとも $|af^2+c| \ge 2^{-18}$

となります。つまり、 $af^2+c$ は $bf$ と比べてはるかに絶対値が大きく、 $g(f) = af^2 + bf + c$ の符号は $af^2+c$ の符号と一致することがわかりました。

ここまでの手順を確認すると、

$b \le -2^{-28}$ ならば手順Aで計算すれば、一回目のfmaが無誤差なので正しく符号判定できる
$b \gt -2^{-28}$ ならば手順Bで計算すれば、一回目のfmaは誤差を生じるかもしれないけれど符号判定には圧倒的に余裕があるので正しく符号判定できる

のような単純な分岐でいけそうな気がします。

場合1： $b \gt 0, c \lt 0, bf+c \lt 0, af^2+c \lt 0$

ラスボスです。といっても、場合2と場合3を倒していれば、難しくはありません（実際には、最初に場合1に取り掛かったので大混乱しました）。

fma(a,f,b)・fma(a,(double)f*f,c)・fma(b,f,c)を全て計算してみます。いずれかが無誤差であれば、対応する手順で計算すればよいので、いずれもが誤差を生じる場合を考えます。

まず、fma(a,f,b)が誤差を生じているため、49ビット整数×2冪である $af$ に対して $b$ を足すと5ビット以上はみ出します。これは、 $af+b \ge 32$ であるか、または $b\lt 2^{-28}$ であることを意味しています。 $af \lt 4$ なので、前者の場合は $b \gt 28$ です。

ここで、 $b \gt 28$ である場合を考えます。 $bf \gt 28$ です。 fma(b,f,c)が誤差を生じているため、 $bf$ と $c$ は絶対値が少なくとも2倍以上違うということがわかり、場合1の時と同様の議論により、 $bf$ と $c$ のうち絶対値が大きい方の符号が $g(f)$ の符号となることがわかります。

逆に、 $b \lt 2^{-28}$ である場合を考えます。 $bf \lt 2^{-27}$ です。 fma(a,(double)f*f,c)が誤差を生じているため、場合3の時の議論と同様に $|af^2+c| \ge 2^{-18}$ であることがわかり、さらに同様に $af^2+c$ の符号が $g(f)$ の符号となることがわかります。

どちらの場合でも、 $g(f)$ の符号を決定することができました。

ここまでの手順を確認すると、

$2^{-28}\le b \le 28$ ならば手順Aで計算すれば、一回目のfmaが無誤差なので正しく符号判定できる
$b \gt 28$ ならば手順Cで計算すれば、一回目のfmaは誤差を生じるかもしれないけれど符号判定には圧倒的に余裕があるので正しく符号判定できる
$b \lt 2^{-28}$ ならば手順Bで計算すれば、一回目のfmaは誤差を生じるかもしれないけれど符号判定には圧倒的に余裕があるので正しく符号判定できる

のような単純な分岐でいけそうな気がします。

まとめ

単精度浮動小数点数係数二次方程式の実数解の correctly-rounded な値を得る、（倍倍精度演算や多倍長演算が不要な）倍精度浮動小数点数の積和演算（fma）を使った手順の概略とそれが正しく動作する理由の簡易な説明を書き下しました。

今回は、単精度浮動小数点数（仮数部がp=24ビット）係数二次方程式を解くための計算に倍精度浮動小数点数（仮数部が53ビット）を用いたうえで、適切に分岐することで、丸め誤差が問題を引き起こすケースを回避しました。しかしながら、53ビットが「十分」であることは示せた（と考えている）ものの、何ビット必要なのかはよくわかりませんでした。少なくとも、(double)f*fみたいなところがあるので50ビット（2(p+1)ビット）無いと破綻します（なので、floatを二つ並べたdouble-float（仮数部が最小で49ビット）で実行することはできなさそうです）。これにさらに追加で3ビットあれば十分というのが今回の結果ですが、もう少し少ないビット数でいけるのかもしれません。

また、今回は指数部が無制限に取れるという仮定を置いて議論を進めました。実際には、単精度浮動小数点数の指数部（±127くらい）の3倍程度あればよいはずで、倍精度浮動小数点数の指数部（±1023くらい）に十分収まるはずです。しかし、面倒だったのでそのような方向性の証明は行いませんでした。

付録： $\forall f \in \mathbb{F}. f \ne f_T \rightarrow g(f_T) \lt g(f)$ の証明

まずは次の補題を示します。

補題1 浮動小数点数同士の除算結果を丸める際、（指数部範囲が十分広ければ）最近接の方向は一意に定まる。
証明背理法による。つまり、「最近接の方向が一意に定まらない場合がある」と仮定して、矛盾を導く。仮数部がpビットだとする。被除数をa、除数をb、商をcとする（ $\frac a b = c$ ）。最近接の方向が一意に定まらないというのは、商が二つの浮動小数点数のちょうど真ん中にあるということである。したがって、必要であれば適当に二冪でスケーリングして、a,b,cのいずれも整数であると仮定してよい。この時、商cは(p+1)ビット奇数×二冪の形をしている。商と除数の積cbを考えると、この数は(p+1)ビット奇数の倍数である。一方で、被除数aはpビット整数×二冪の形であり、(p+1)ビット奇数の倍数たりえない。 $a = cb$ なので、これは矛盾である。よって仮定「最近接の方向が一意に定まらない場合がある」は誤りであり、最近接の方向は一意に定まる。

さて、 $\forall f \in \mathbb{F}. f \ne f_T \rightarrow g(f_T) \lt g(f)$ の証明に戻ります。一般性を失わずに $a > 0$ とします。

実数に関する全称命題 $\forall x \in \mathbb{R}. x \ne x_T \rightarrow g(x_T) \lt g(x)$ は、高校数学で簡単に証明できたことを思い出します（ $g(x) = (x-x_T)^2+g(x_T)$ と書けて、 $x \ne x_T$ ならば第一項が正になる、ということを使って証明できます）。示したい命題は、これの浮動小数点数版ですが、似た感じで示せます。

除算の最近接の方向は一意に定まるので、 $\forall f\in\mathbb{F}. f \ne f_T \rightarrow |f_T-x_T| \lt |f - x_T|$ が成り立ちます……(1)。 $|f_T-x_T| \lt |f - x_T|$ ならば $(f_T-x_T)^2 \lt (f - x_T)^2$ ……(2)です。 $g(f) = (f-x_T)^2+g(x_T)$ と書くと、 $g(f) - g(f_T) = (f - x_T)^2 - (f_T-x_T)^2$ であり、右辺は(1)(2)により $f \ne f_T$ ならば正です。したがって、 $f \ne f_T$ の時、 $g(f) - g(f_T) \gt 0$ 、つまり $g(f_T) \lt g(f)$ が示せました。