数値計算
誤差
桁落ち

今回取り上げるのは、計算機上で小数を扱う際に重要な問題となってくる『桁落ち』という現象である。本記事は元にした拙記事

ななついろ☆ドロップアウト - あざらしとペンギンの問題

の口調を変えたものである1。本記事の読者が何らかの試験に落ちないことを祈ってタイトルを変えたことは言うまでもない。なお、本記事の内容はほとんど実装言語に依存しないことを断っておく。

結論から言えば、『桁落ち』とは「近接した値同士で引き算をすると結果の不確かさが異常に大きくなる」ということだ。かなり微妙な言い回しだが、ともかく「引き算には気をつけろ」といった意味だと今は思っておいていただきたい。


小数と誤差

計算機上では普通は2進小数が使われるのだが、『桁落ち』および『情報落ち』は特に基数に依存するものではないので、我々にとって馴染み深い10進小数2で説明する。


有効数字と誤差

現代において我々は整数でない数を小数で表し、特に無限小数となる場合は有限桁の小数で打ち切ることを日常的に行っている。例えば、

\begin{align*}

\sqrt{2} &= 1.41421356\dots \\
\pi &= 3.14159265\dots
\end{align*}

のように。あるいは、打ち切る桁より下の位に丸め(四捨五入など)、切り上げ、切り下げのいずれかを行って

\begin{align*}

\sqrt{2} &= 1.4142 \times 10^0 \\
\pi &= 3.1416 \times 10^0
\end{align*}

のように表したりもするだろう。これらの例は6桁目を四捨五入したものだ。前にある小数は『仮数部』(significand)、10 の肩に乗っている数は『指数部』(exponent) という。仮数部は必ず 1 以上 10 未満の範囲に入っていなければならない。この形式は自然科学や工学において不確さを含む値を表記するときと同じである。ただし、この等号は真の等号ではないことに注意しなければならない。

上の例では5桁の数字に意味があるので、『有効数字』5桁の値と言う。

有効数字の桁数は先頭の0でない桁から数える。例えば、

\frac{23}{4567} = 0.0050361287\dots = 5.0361 \times 10^{-3}

のように、元の小数の6桁目を丸めたのでなく、最初に0でない数字が出る桁から数えて6桁目を丸めた値を、有効数字5桁の値として採用する。なお、このような形式の小数は『浮動小数点数』と呼ばれる。計算機における小数の実装のほとんど(C言語の double など)は2進浮動小数点数である。

数に対して有効数字への丸めを行うと、当然のことながら真の値に対して誤差が発生する。これは『丸め誤差』と呼ばれる。例えば、

\pi = 3.1416 \times 10^0

は真の値から $0.0000073464\dots$ だけ大きくなっている。あるいは、これを真の値から $:0.00023384\dots \%$ だけ大きいと表現することもある。前者を『絶対誤差』、後者を『相対誤差』と言う。それぞれの定義は次に回すとして、実数を有限の桁数で表す以上は、誤差は常について回るものだ。

計算を有限の桁数で行うことに起因する誤差は『打ち切り誤差』と呼ばれる。一般に、誤差というと測定誤差などの元々の数値に含まれる誤差も含まれるが、それらは今回扱う問題の範囲外とする。

なお、切り下げ(非負の場合。負の場合は切り上げ)の場合を除いて、有効数字はその桁数の真の値の表示とは一致しないことがある。上の例では5桁目の数字が異なっている。

このように「有効数字xx桁」という表現の意味するところは曖昧である。あくまで、精度を考える上で『不確かさ』を含む位以下を除いた数字を数えているだけであり、その桁数がどういう意味を持つかは、誤差の意味や計算上の丸めモードなどによって異なる。しかし、この記事においては以後も「有効数字xx桁」という表現を使うことにする。桁落ちを説明するにはそれで十分だからである。


絶対誤差と相対誤差

さて、少し前の文で『絶対誤差』と『相対誤差』という言葉が現れた。誤差といえば、例えば、

\begin{align*}

\sqrt{2} &= 1.4142 \times 10^0 \\
\pi &= 3.1416 \times 10^0 \\
\frac{23}{4567} &= 5.0361 \times 10^{-3}
\end{align*}

のそれぞれの右辺の値にはすべて丸め誤差が含まれている。このような誤差は大きく2種類に分けられるということだ。

その2種類とはすなわち『絶対誤差』と『相対誤差』であり、次のように定義される。ここで、$x^*$ は真の値、$\tilde{x}$ は誤差を含む値とする。

絶対誤差:

\epsilon_a = \left| \tilde{x} - x^* \right|

相対誤差:

\epsilon_r = \frac{\epsilon_a}{x^*} = \frac{\left| \tilde{x} - x^* \right|}{x^*}

相対誤差に関する見やすい例では、真の値 100 に対する 100.003、それと真の値 100000 に対する 100003 を誤差の意味で同じと見なす。双方とも言葉では「真の値より 0.003% だけ大きい」と言うことができる。

これらの二つの誤差のうち、浮動小数点数では絶対誤差ではなく相対誤差の方が重要である。というのも、数の大きさが指数部に依存するためだ。前文の例では、有効数字5桁として

\begin{align*}

100.003 &= (1.0000 + 0.00003) \times 10^{2} \\
100003 &= (1.0000 + 0.00003) \times 10^{5}
\end{align*}

であり、両者は仮数部の誤差の意味で全く同じであることがわかる。ここで、1.0000 と 1 が異なる意味を持ち、両者を区別しなければならないことは、学校の理科の授業などで注意されたはずだ。あくまで有効数字の意味の取り方ひとつとして、真の値 $x^*$ と誤差を含む値 $\tilde{x}$ との間で

\begin{align*}

\tilde{x} = 1.0000 &\Leftrightarrow 0.99995 \le \tilde{x} < 1.00005 \\
\tilde{x} = 1 &\Leftrightarrow 0.5 \le \tilde{x} < 1.5
\end{align*}

を意味するものと仮定する。このとき、誤差を含む値の『不確かさ』を有効数字の桁数で表しているため、必然的に区別が必要なのだ。ここでは有効数字5桁としたので、6桁目を四捨五入すると、

\begin{align*}

10.003 &= 1.0000 \times 10^{2} \\
100003 &= 1.0000 \times 10^{5}
\end{align*}

となる。

この通り、浮動小数点数では誤差を相対誤差で見ることにより、仮数部に関するものに限定することができる。逆に、もし計算機が5桁までの値しか扱えないならば、仮数部の 0.00005 の不確かさを除去することができない3

このような計算機の有限性の制限のため生じる誤差の範囲を『マシンイプシロン』と呼ぶことがある。ただし、この言葉は絶対誤差の意味でも使われるため、曖昧な表現であることに注意が必要である。

ところで、今まで『誤差』と『不確かさ』という二つの言葉を使ってきました、この記事では両者を明確に区別して使うことにする。


『誤差』とは本節の最初に定義された『絶対誤差』または『相対誤差』のみを指すものとする。


対して、有効数字から外れるために誤差として扱われる部分の範囲を『不確かさ』と呼ぶことにする。有効数字という言葉が曖昧性を含むように、不確かさという言葉も曖昧性を含むことは注意しておく。先にも断った通り、この曖昧性は本記事において重要ではない。


計算による不確かさの拡大

さて、敢えて『不確かさ』という曖昧な言葉を使うことのいくつかの理由のひとつは、丸め誤差や打ち切り誤差を含む値同士で計算を行ったとき、一般に誤差の範囲が拡大するためだ。例えば、有効数字5桁、6桁目を四捨五入するとして、小数を用いた計算では

\frac{1}{3} + \frac{1}{3} = 3.3333 \times 10^{-1} + 3.3333 \times 10^{-1} = 6.6666 \times 10^{-1}

となるところが、真の値を四捨五入すると

\frac{1}{3} + \frac{1}{3} = \frac{2}{3} = 6.6667 \times 10^{-1}

となることなどがある。第一の計算値については

6.6666\dots \times 10^{-1} < \left(\frac{2}{3} - 0.00005 \right) \times 10^{-1}

であるため、6桁目の四捨五入によって許容される誤差の範囲に収まっていないことが判る。このように、有効数字の桁数はあくまで精度の目安であって、厳密な誤差の範囲を表すものではないことに注意が必要である。

ここまで、誤差を含む値を表すのに考えなく等号を使ってきたが、当然のことながらその等号は正しくない。また、計算による不確かさの拡大があるため、計算順序は無視できない問題である。本題の『桁落ち』とは、最初に書いた通り、計算による不確かさの拡大を意味する。

そこで、計算の推移を $\to$、真の値から小数への変換を $\Rightarrow$ 単に値が近いことを $\approx$、値が右に示される区間に入っていることを $\in$ で表すことにする。上の例では、

\begin{align*}

\frac{1}{3} + \frac{1}{3} &\Rightarrow 3.3333 \times 10^{-1} + 3.3333 \times 10^{-1} \to 6.6666 \times 10^{-1} \\
\frac{1}{3} + \frac{1}{3} \to \frac{2}{3} &\Rightarrow 6.6667 \times 10^{-1} \\
\frac{2}{3} &\approx 6.6667 \times 10^{-1} \\
\frac{2}{3} &\in (6.6667 \pm 0.00005) \times 10^{-1}
\end{align*}

のように表される。


相対誤差の上界

ここでひとつ相対誤差について注意すべきことがある。本来、相対誤差の定義では真の値から見た誤差を測るべきところが、これまでの例では誤差を含む値から見た真の値の存在する範囲を示しているのだ。当然だが両者は同じではない。実を言えば我々は数値計算によって相対誤差それ自体を厳密に知ることはできない。というのも、我々が計算機上で扱う値は誤差を含む値の方だからだ。

とは言っても、相対誤差の取りうる値の上界を見積もることは可能だ。まず、誤差を含む値 $\tilde{x}$ を仮数部と指数部に分けて

\tilde{x} = m \times 10^e

とする。ここで $1 \le |m| < 10$ かつ $e$ は整数である。

このとき、相対誤差の上界は、$\epsilon > 0$ として

x^* \in (m \pm \epsilon) \times 10^e

のとき

\epsilon_r \le \frac{\epsilon \times 10^e}{(|m| - \epsilon) \times 10^e} = \frac{\epsilon}{|m| - \epsilon} < \frac{\epsilon}{1 - \epsilon}

となる。

ここで $|m| = 1$ と置いて上界を取ったことはいささか乱暴に思われるだろう。もし $|m| = 9$ ならば、誤差を9倍以上緩く見積もっていることになる。

この上界の緩さを狭めるには、2進小数を使うことが最良の方法となる。なぜなら、2進小数では $|m| = 1$ に限定されるためだ。その場合、誤差の見積もりは高々2倍以内に収まる。計算機では2進小数が使われるのが普通だが、機械的な実装だけでなく誤差の見積もりの意味でも、2進小数を使う方が有利なのだ。

そのような事情はあるものの、この記事ではあくまでわかりやすさのために、我々にとって馴染み深い10進小数を使い続けることにする。

10進で有効数字 $k \ge 6$ 桁4、$(k + 1)$ 桁目を四捨五入するとした場合、$\epsilon = 5 \times 10^{-(k + 1)}$ なので、相対誤差の上界を乱暴に計算すると

\epsilon_r < \frac{\epsilon}{1 - \epsilon} < 1.000006 \times 5 \times 10^{-k} =  5.00003 \times 10^{-(k + 1)}

となる。


桁落ち

さて、ここまでは準備で、本題はここからだ。

最初に言ったように、『桁落ち』とは「相対的に近接した値同士で引き算をすると相対的な不確かさが異常に大きくなる」ということだ。問題点をより明らかするために、『近接』と『不確かさ』の前に『相対的な』を冠したことに注意されたい。

本節においても、小数は10進有効数字5桁で6桁目を四捨五入するものとする。このように有限個の数字で計算を行うことは、この問題に関して本質的である。


桁落ちの例

まず、次のような大きさの極端に異なる数の足し算を行ってみる。

\sqrt{2} + \frac{23}{4567} \Rightarrow (1.4142 + 0.0050361) \times 10^0 \to 1.41923 \times 10^0 \to 1.4192 \times 10^0

この計算により、誤差そのものは最初の誤差の上限 0.00005 より増えている可能性があるが、この時点では相対的な不確かさの拡大は小さなものに収まっている。実際、

\displaystyle \sqrt{2} + \frac{23}{4567} = 1.4192496911228\dots \in (1.4192 \pm 0.00005) \times 10^0

なので、厳密な相対誤差の意味でも想定内に収まっていることが判る。

問題は、先程求めた値 $1.4192 \times 10^0$ から $\sqrt{2} = 1.4142 \times 10^0$ を引くことによって発生する。注意すべきはこの時点で持っている値が $1.4192 \times 10^0$ と $1.4142 \times 10^0$ 以外の何物でもないことだ。あくまで数値のみで計算を行うため、記号としての $\sqrt{2}$ などは忘却の彼方にあるというわけだ。それでは実際に引いてみましょう。

1.4192 \times 10^0 - 1.4142 \times 10^0 \to 0.0050 \times 10^0 = 5.0 \times 10^{-3}

はい。有効数字2桁になりました。5 $5.0 \times 10^{-3}$ と $5.0000 \times 10^{-3}$ の違いを思い出されたい。真の値を5桁に丸めたものは $5.0361 \times 10^{-3}$ なので、不確かさの拡大は厳密な相対誤差の拡大とも合致している。

改めて全体の計算を書くと、

\begin{align*}

\left( \sqrt{2} + \frac{23}{4567} \right) - \sqrt{2} &\Rightarrow \left( (1.4142 + 0.0050361) - 1.4142 \right) \times 10^0 \\
&\to (1.4192 - 1.4142) \times 10^0 \to 5.0 \times 10^{-3}
\end{align*}

となる。

以上のことをもって『桁落ち』が生じたと言う。


桁落ちはとにかく悪い

桁落ちの何が問題かというと、計算値の相対的な不確かさが他の四則演算に比べて異常に増えることだ。

馴染み深い例は微分を数値的に計算するときに現れる。まず微分係数の定義式が

f'(x) = \lim_{h \to 0} \frac{f(x + h) - f(x)}{h}

であることは言うまでもないが、この式より(悪い)自明な計算方法として $h$ を適当に小さな値に取って式を評価するということが考えられる。ここでも計算は全体を通して有効数字5桁で行われるとする。

以下、$f(x) = \sqrt{x}$ として $x = 2$ ちょうどでの微分係数を数値計算することを考える。真の値は

f'(2) = \frac{1}{2 \sqrt {2}} = 0.35355339059 \in 3.5355339 \times 10^{-1}

である。

$h$ は適当に小さな値として $h = 0.001$ ちょうどを取るとしよう。このときの計算は

\frac{\sqrt{2 + 0.001} - \sqrt{2}}{0.001} \Rightarrow \frac{\sqrt{2.0000 + 0.0010} - \sqrt{2.0000}}{0.0010000} \to (\sqrt{2.0010} - \sqrt{2.0000}) \times 1000.0 \to (1.4146 - 1.4142) \times 1000.0 \to 0.40000 = 4.0000 \times 10^{-1}

となり、大きく異なった値が出ることがわかる。

このように、桁落ちは不確かさが異常に大きくなることにより、どのような意味でも精度の悪い計算結果を生み出してしまうのである。

なお、よく知られた方法として、中心差分を用いた公式

f'(x) = \lim_{h \to 0} \frac{f(x + h) - f(x - h)}{2 h}

を用いて計算すると良い計算値が出せるというものがある。これによれば、同じ有効数字5桁のもとで $f'(2)$ は次のように計算される:

\frac{\sqrt{2 + 0.001} - \sqrt{2 - 0.001}}{2 \times 0.001} \Rightarrow \frac{\sqrt}2.0000 + 0.0010000} - \sqrt{2.0000 - 0.0010000}}{2 \times 0.0010000} \to (\sqrt{2.0010} - \sqrt{1.9990}) \times 500.00 \to (1.4146 - 1.4139) \times 500.00 \to 0.3500 = 3.5000 \times 10^{-1}

先程より少しマシな結果が得られている。結果が2桁まで合っていることが見て取れる。

それにはちゃんとした理由がある。次の図より直感的にはわかってもらえることを期待する。

mid.png

しかし、この設定でこれ以上良い結果を得ることはおそらく不可能だろう。なぜなら、有効数字5桁では $2.000 + 0.0010000 \Rightarrow 2.0010$ となるように、$0.0010000$ のうち2桁より下の情報が失われてしまうためだ。それを次節で説明する。


桁落ちと情報落ち

『桁落ち』と関係する言葉として『情報落ち』というものがある。これは、桁数が有限であるために元々の数が持っていた情報が計算によって失われることを指して言う言葉である。例えば、先に出した例:

\sqrt{2} + \frac{23}{4567} \Rightarrow (1.4142 + 0.0050361) \times 10^0 \to 1.41923 \times 10^0 \to 1.4192 \times 10^0

では後ろの桁の情報 361 が失われていることから『情報落ち』が発生している。これは一般には『積み残し』と呼ばれ、桁違い6に大きさの異なる数を足し引きしたときに発生するものである。

しかし、この段階では有効桁数はほとんど失われていないことに注意が必要である。つまり、『情報落ち』こそ発生しているが、『桁落ち』が発生したとまでは言えない。

『桁落ち』が発生するのは近接した値の引き算を行ったときである。上の結果から $\sqrt{2} \in 1.4142 \times 10^0$ を引いて

1.4192 \times 10^0 - 1.4142 \times 10^0 \to 0.0050 \times 10^0 = 5.0 \times 10^{-3}

としたとき、初めて『桁落ち』が発生する。

『桁落ち』の原因は上位桁が等しい数の有効数字がキャンセルされてしまうことである。すなわち、データとして見れば『情報落ち』の範疇に含まれる。しかし、『桁落ち』は単なる情報落ちに留まらず、有効数字が著しく失われることを特に指して言う言葉である。有効数字がほとんど失われない『情報落ち』も避けるに越したことに違いはないが、『桁落ち』はそれよりずっと重要な問題であり、小さくない犠牲を払ってでも避けられるべきである。


桁落ちの回避

桁落ちを回避する方法は、とにもかくにも、その原因となる「近接した値同士の引き算」を行わないようにすることである。それには、計算する対象についての事前知識があると、有利に進めることができる場合がある。


完全に同じ値の引き算を記号的に除去する

前節のひとつ前例では、完全に同じとわかっている値を記号的に引いてしまってから、小数としての計算を行うことによって、桁落ちを回避できる。

\begin{align*}

\frac{1}{\left( \sqrt{2} + \frac{23}{4567} \right) - \sqrt{2}} &\to \frac{1}{\left( \sqrt{2} - \sqrt{2} \right) + \frac{23}{4567}} \to \frac{4567}{23} \\
&\Rightarrow \frac{4.5670 \times 10^3}{2.3000 \times 10^1} \to 1.9857 \times 10^2
\end{align*}

一見すると馬鹿馬鹿しい例だが、記号のまま処理できる部分を先に記号のまま処理しておくことは、桁落ちを回避する上で重要なことである。あるいは、数学的、物理的な知見からの不変性、対称性が利用できるかもしれない。


とにかく引き算を使わずに済ませる

桁落ちを回避するための最も一般的な方策は、足し算と引き算の選択肢では必ず足し算を選ぶこと、また、足し算と乗除算を組み合わせることで引き算を回避できるならばそれを行うことだ。

最も馴染み深いのは二次方程式の解の公式だろう。中学校で習うように、二次方程式

ax^2 + bx + c = 0\;(a \neq 0)

の二つの解は

x = \{ r, s \} = \frac{-b \pm \sqrt{b^2 - 4ac}}{2a}

という公式によって求められる。ここでは二つの解が実数である場合のみを扱い、$r \ge s$ とする。

以下において、コロン付き等号 $:=$ を用いて定義されている式は計算式であって、それを評価して得られる値とは区別されるものとする。

$b > 0$ のとき、$r$ の直接の計算式

r_1 := \frac{-b + \sqrt{b^2 - 4ac}}{2a}

において、分子の計算

-b + \sqrt{b^2 - 4ac}

は引き算となる。ひとつの場合として $b$ が平方根の中の値に近い場合、すなわち $b^2 \gg 4ac$ の場合に桁落ちが生じる。もうひとつの場合としては、$b^2 \approx 4ac$ のときは平方根の中の引き算で桁落ちが発生する。ここでは前者について考えることにする。

とりあえず、一方の解 $r$ を得るために $r_1$ を直接計算することは可能な限り避けなければならない。

他方、$s$ の計算

s_1 := \displaystyle x = \frac{-b - \sqrt{b^2 - 4ac}}{2a}

においては、分子は負の値同士の足し算となるため、桁落ちが生じないことが判る。

よって、最初に $s_1$ だけを計算し、もう一方を解と係数の関係

rs = \frac{c}{a}

から

r_2 := \frac{c}{a s_1}

の計算によって求めることで、桁落ちを回避することが可能です。

例としてこの場合の問題が生じる二次方程式 A

x^2 + 10 x + 0.001 = 0

について、$r_1$ と $r_2$ の式を32ビット浮動小数点数(10進有効数字7桁程度)の範囲内で計算してみる。真の解の小数値は次の通り:

r = -1.0000100002000050\dots \times 10^{-4}, s = -9.9998999989999799\dots

コードと計算結果:

[https://ideone.com/FFul9l]

\begin{align*}

r_1 &= -1.001358 \times 10^{-4}, s_1 = -9.999900 \times 10^0 \\
r_2 &= -1.000010 \times 10^{-4}
\end{align*}

真の値と比較すると、$r_1$ はかなり精度が悪く、$r_2$ は10進有効数字7桁程度では十分な精度を持った値であることが判る。$s_1$ について後ろの桁の値が真の値と異なっているが、誤差の意味では精度から想定される範囲を逸脱するものではない。


分母、分子の有理化

二次方程式の解の公式のような平方根の引き算を伴う計算では、分母または分子の有理化を行うことによって、桁落ちが回避できる場合がある。

r_1 := \frac{-b + \sqrt{b^2 - 4ac}}{2a}

の分子

-b + \sqrt{b^2 - 4ac}

を有理化するとは、分母と分子の両方に

-b - \sqrt{b^2 - 4ac}

を掛けることである。少し計算すると、分子は

b^2 - (b^2 - 4ac) = 4ac

分母は

-2a \left( b + \sqrt{b^2 - 4ac} \right)

となって、

r_3 := -\frac{2c}{b + \sqrt{b^2 - 4ac}}

という計算式が得られる。

再び前節の二次方程式 A

 x^2 + 10 x + 0.001 = 0

について、$r_1$ と $r_3$ の式を32ビット浮動小数点数(10進有効数字7桁程度)の範囲内で計算してみる。

真の値:$r = -1.0000100002000050\dots \times 10^{-4}$

コードと計算結果:

[https://ideone.com/FFul9l]

\begin{align*}

r_1 &= -1.001358 \times 10^{-4} \\
r_3 &= -1.000010 \times 10^{-4}
\end{align*}

期待通り $r_3$ は十分な精度の値になっていることが判る。

ところで、$r_3$ についてよく見ると、前節と同じ

s_1 := \frac{-b - \sqrt{b^2 - 4ac}}{2a}

を用いて

r_3 = \frac{2c}{2a s_1} = \frac{c}{a s_1}

と書けることが判る。すなわち、$r_3$ は前節の $r_2$ をひとつの式にまとめたものと見なすことができる。計算手順の違い(コンパイラ依存なので特定はできない)はあるものの、桁落ちが生じないように改良された計算では大きな問題にはなっていない。


級数展開を使う

前の二次方程式の解の公式の問題を、級数展開を使って解決することも考えられる。三度、問題の引き算(平方根の外)

-b + \sqrt{b^2 - 4ac}

で桁落ちが発生するときというのは、$b^2 \gg 4ac$ となるときである。少し変形すると、

-b + \sqrt{b^2 - 4ac} = -b + b \sqrt{1 - \frac{4ac}{b^2}} = b \left( \sqrt{1 - \frac{4ac}{b^2}} - 1 \right)

において

\frac{4ac}{b^2} \ll 1

となるときとも言える。このとき、一般化二項定理により得られる級数展開

\sqrt{1 - \frac{4ac}{b^2}} = 1 - \frac{2ac}{b^2} - \frac{2 a^2 c^2}{b^4} \dots

は良い近似になっていると考えられる。式の中に引き算が現れるが、近接した値でないことが保証されているので、桁落ちが問題になる場合ではない。

というわけで、級数を上の3項で打ち切って代入すると、

-b + \sqrt{b^2 - 4ac} = -\frac{2ac}{b} - \frac{2 a^2 c^2}{b^3}

から

r_4 := -\frac{c}{b} \left( 1 + \frac{ac}{b^2} \right)

という計算式が得られます。この式において括弧の中は引き算になる可能性があるが、仮定より近接した値でないことが保証されているので、それによる桁落ちは問題にはならない。

というわけで、三度、問題の二次方程式 A

x^2 + 10 x + 0.001 = 0

について、$r_1$ と $r_4$ の式を32ビット浮動小数点数(10進有効数字7桁程度)の範囲内で計算してみる。

真の値:$r = -1.0000100002000050\dots \times 10^{-4}$

コードと計算結果:

[https://ideone.com/FFul9l]

\begin{align*}

r_1 &= -1.001358 \times 10^{-4} \\
r_4 &= -1.000010 \times 10^{-4}
\end{align*}

期待通り $r_4$ は十分な精度の値になっていることが判る。計算式 $r_4$ の導出には近似を用いているので $r_2, r_3$ とは記号的にも異なっているが、結果としては同程度の精度の値が出ていることは注目に値するだろう。


加速法を使う

上で例として出したように、微分方程式の差分解法のような、微分係数の差分近似を含む計算を行うとき、桁落ちは宿命的な問題となる。この前の記事で微分係数を加速法によって高精度で求める方法を説明してるが、それだけでひとつの記事になってしまうので、現時点ではリンクを貼るだけに留めておく。

アクセルほど微分精度は加速していくよ! - あざらしとペンギンの問題


引き算には気をつけろ

桁落ちを回避する方法はここに書いた以外にも様々あるが、どんな場合でも上手くいくような万能の処方箋はない。とにもかくにも引き算に気をつけるという意識を持つことが第一である。可能ならば引き算を行わないようにすること、回避できない場合は可能な限り回数を減らす、また、ソートを行って順序を変えるなど、工夫して計算を行うことが求められる。


有限性の壁

計算機において有効数字、桁落ち、情報落ちといった概念が出てくる理由は、物質的な制約のために扱えるデータ量が限られることにある。特には浮動小数点演算装置(FPU)のレジスタが(符号、指数を含めて)少ない桁数の分しかないということである。レジスタを大きくすれば桁数を増やすことはできるが、その分だけ数のデータを保持し、また、処理するためにより大きな、あるいは集積度の高い回路が必要となる。回路の素子は原子で構成されていて、宇宙には有限個の原子しかないので、無限の桁数を扱うことは当然ながら不可能だ。それは、数のデータをチップの外、例えばハードディスクやクラウドストレージに入れたとしても同じである。

さらに言えば、データの読み出しはつまるところ測定なので、物理学的な不確かさが本質的に有効数字を限定しているとも考えられる。この不確かさは測定誤差ばかりでなく、量子力学的、統計力学的な不確かさ、つまり、この世界においてどうやっても克服できない不確かさが含まれる。桁数というものが名目上存在しないアナログ計算機でさえ、有効数字の束縛から逃れ得るものではないのである。


まとめ

桁落ちとは、浮動小数点数演算において近接した値の引き算を行うことによって有効数字が失われ、計算値の不確かさが異常に大きくなることであった。桁落ちは数値計算における重要な問題であり、可能な限り回避しなければならないものである。本記事ではいくつかの方法を示した。前の記事に書いた加速法も桁落ちを回避する手段となり得る。その他様々な方法があるが、万能の処方箋はなく、問題を見て工夫すべき場合も少なくない。本記事の教訓を一言で言えば、

** 引き算には気をつけろ **

これは数値計算を行う上で常に意識しなければならないことである。

桁落ち、情報落ちは有限の世界で計算を行う以上は避けられない問題である。一見すると数値計算における細かく面倒な問題のようで、実は計算の本質に関わる問題であるというのが私の思うところである。





  1. 本記事の目的のひとつは機械的な置換を検証するためである。数式に関しては概ね一括変換しても問題なさそうだ。ただ、逐語的な書き換えを行うにあたって、記事本体に少なくない誤りを発見した。これを機械に発見してもらうことはまだ難しいだろう。また、どちらに投稿するのが相応しいか、まだ判断材料に欠けるところである。 



  2. 他の記数法を用いる文化圏では、10進小数が馴染み深いとは言えないかもしれない。 



  3. 多倍長演算によってより多くの桁の小数を扱うことは可能だが、計算の速度は極端に落ちる。 



  4. IEEE 754 の定める 32 ビット浮動小数点数は10進で有効数字7桁程度 



  5. ここでは意図的に口調を変えなかった。 



  6. ここでは読んで字のままの意味