うさぎでもわかる計算機システム Part04 桁落ち・情報落ち・丸め誤差・打ち切り誤差について [基本情報対応]

大きい解(理論値：-0.001)を求めるときに誤差が発生していますね。ここで公式の分子を思い出してみましょう。\( b = 1000.001 \), \( c = 1 \) を代入し、 \[ -1000.001 \pm \sqrt - 4> \] となります。さらに\[ \sqrt - 4> \fallingdotseq 1000.001 \] ですね。よって、\[ -1000.001 + \sqrt - 4>\] の計算時にほぼ同じ2つの値の減算を行っているため桁落ちが発生します*1。

\( b \) が負となるバージョンもやってみましょう。\[x^ - 1000.001x + 1 = 0\]の解をもとめます。

2つの解は、理論的には\[x^ - 1000.001x + 1 = (x - 0.001)(x - 1000) = 0\]となり、\( x_ = 0.001 \), \( x_ = 1000 \) となります。

（ここを書き換え）の行の部分を double a = 1, b = -1000.001, c = 1; と書き換えると、以下の実行結果を得られます。

☆解の公式で求めた解☆ 1つ目 … 0.00099999999997635314 2つ目 … 1000.00000000000000000000

小さい解(理論値：0.001)を求めるときに誤差が発生していますね。また公式の分子を思い出してみましょう。\( b = -1000.001 \), \( c = 1 \) を代入し、 \[ 1000.001 \pm \sqrt - 4> \] となります。さらに\[ \sqrt - 4> \fallingdotseq 1000.001 \] ですね。よって、\[ 1000.001 - \sqrt - 4>\] の計算時にほぼ同じ2つの値の減算を行っているため桁落ちが発生します*2。

2次方程式の解の公式適応時の桁落ち条件

2次方程式 \[ a x^ + b x + c =0 \]の解 \( x_ \), \( x_ \), \( (x_ \lt x_ ) \) を解の公式 \[x = \frac \] で求める場合

\( b \gt 0 \) のとき\[ b \fallingdotseq \sqrt -4 ac> \] のとき、大きい解である \( x_ \) で誤差発生。

\( b \lt 0 \) のとき\[ b \fallingdotseq - \sqrt -4 ac> \] のとき、小さい解である \( x_ \) で誤差発生。

桁落ちをなくすためには 2次方程式の解と係数の関係

2次方程式 \[ a x^ + b x + c =0 \]の解 \( x_ \), \( x_ \), \( (x_ \lt x_ ) \) とすると、\[ x_ + x_ = - \frac \ \ \ x_ x_ = \frac\]が成立する。

\( b > 0 \) のとき桁落ちが発生する可能性があるのは、大きい解の \( x_ \) であるので、小さい解 \( x_ \) を解の公式で算出後、大きい解 \( x_ \) を次の式で求めます。\[x_ = \frac\]

\( b \lt 0 \) のとき桁落ちが発生する可能性があるのは、小さい解の \( x_ \) であるので、大きい解 \( x_ \) を解の公式で算出後、小さい解 \( x_ \) を次の式で求めます。\[x_ = \frac \]

#include #include int main() < double a = 1, b = 1000.001, c = 1; // ax^2 + bx + c = 0 の入る (ここを書き換え) double d = pow(b,2) - 4.0 * a * c; // 判別式 double ans_small = [3]-1.0 * b) - sqrt(d / (2.0 * a); double ans_big = [4]-1.0 * b) + sqrt(d / (2.0 * a); printf("☆解の公式で求めた解☆\n"); printf("1つ目 … %27.20f\n",ans_small); printf("2つ目 … %27.20f\n\n",ans_big); double ans2; if(b >= 0) < ans2 = c / (a *ans_small); printf("☆解と係数の関係から求めた解☆\n"); printf("1つ目 … %27.20f\n",ans_small); printf("2つ目 … %27.20f\n",ans2); >else < ans2 = c / (a *ans_big); printf("☆解と係数の関係から求めた解☆\n"); printf("1つ目 … %27.20f\n",ans2); printf("2つ目 … %27.20f\n\n",ans_big); >return 0; >

実行結果： \( b = 1000.001 \gt 0 \) の場合（他の変数はプログラム通り）

☆解の公式で求めた解☆ 1つ目 … -1000.00000000000000000000 2つ目 … -0.00099999999997635314 ☆解と係数の関係から求めた解☆ 1つ目 … -1000.00000000000000000000 2つ目 … -0.00100000000000000000

実行結果： \( b = -1000.001 \lt 0 \) の場合（他の変数はプログラム通り）

☆解の公式で求めた解☆ 1つ目 … 0.00099999999997635314 2つ目 … 1000.00000000000000000000 ☆解と係数の関係から求めた解☆ 1つ目 … 0.00100000000000000000 2つ目 … 1000.00000000000000000000 (2) 具体例2 数値微分

もう少し具体例を出してみましょう。例えば、\( f(x) = \sin x \) の \( x = 1 \) における微分係数 \( f'(1) \) を計算機上で求めてみます。

\( x = x_ \) における微分係数微分の定義式は \[ \lim_ \frac \] ですね。

プログラムをすることによって求めた微分係数と、解析的（実際に導関数を求め、 \( x = 1 \) を代入、今回の場合は \( f'(x) = \cos x \) に \( x = 1 \) を代入）に求めたプログラムの誤差を調べます。

ただしプログラムの場合、解析的に \( h \to 0 \) を考えるのは難しいため、\( h \) を有限の範囲内で限りなく小さくして考えます*3。

#include #include int main() < float x = 1; float d1; double d0 = cos(x); // cos(1) の理論値 float y0 = sin(x),h,y1; for(int i = 0; i %2.8f\n",i,fabs(d1 - d0)); > > h = 2^(- 0) -> 0.47247586 h = 2^(- 1) -> 0.22825423 h = 2^(- 2) -> 0.11024764 h = 2^(- 3) -> 0.05392936 h = 2^(- 4) -> 0.02663901 h = 2^(- 5) -> 0.01323321 h = 2^(- 6) -> 0.00659564 h = 2^(- 7) -> 0.00329211 h = 2^(- 8) -> 0.00163653 h = 2^(- 9) -> 0.00081256 h = 2^(-10) -> 0.00038531 h = 2^(-11) -> 0.00014117 h = 2^(-12) -> 0.00001910 h = 2^(-13) -> 0.00022504 h = 2^(-14) -> 0.00071332 h = 2^(-15) -> 0.00071332 h = 2^(-16) -> 0.00266644 h = 2^(-17) -> 0.00657269 h = 2^(-18) -> 0.00657269 h = 2^(-19) -> 0.02219769 h = 2^(-20) -> 0.02219769

確かに、hが小さいうちは、h を小さくすればするほど誤差は小さくなっていますね。しかし、\( h \gt 12 \) からはhが小さくなればなるほどなぜか誤差が逆に大きくなっています。なぜでしょうか？

(1) 打ち切り誤差

最初のうちは、hの値を 1/2 にするたびに誤差が約 1/2 になっていますね。hをより0に近づければ近づけるほど、「0に限りなく近づけた値」に近づき、打ち切り誤差は減少します。

(2) 桁落ち

ほぼ等しい2つの数の引き算 をすることで桁落ちは発生しますね。

\( h \) を小さくすればするほど \( f(x_+h) \), \( f(x_) \) の差は小さくなっていきますね。

小さくなった状態で \( f(x_+h) - f(x_) \) を計算すると、ほぼ等しい2つの数の減算を行うことで計算結果の有効数字が減少し、誤差の原因となります。

以上2つが、数値微分における誤差の原因となります。数値微分をする場合、結果との誤差が小さくなり続けてる中で最も小さい \( h \) の値を採用するのが数値微分における誤差を小さくするコツといえます。今回の場合だと \( h = 12 \) のときとなります。

また、実際に数値微分をする際には、\( \sin x \) のような解析的な微分結果が自明なものではなく、複雑な関数で解析的に微分ができないようなものに行います。

そのような場合には、\[\lim_ \frac \] など、\( h \) の値を変えて計算を行い、上で行った\[ \lim_ \frac \] との差を比べ、差が小さくなり続けてる中で最も小さい \( h \) の値を採用するのがただしい数値微分の結果を得るコツとも言えます。

２．情報落ち

(1) 概要

情報落ちは、 差が極端にある大きい数と小さい数の加算/減算のときに小さい数の計算が反映されない 現象を表します。

10進数の場合で例を出してみましょう。\[ 3.0453 \times 10^ + 8.1019 \times 10^ \]このような計算をするとき、皆さんは大きい方の指数部にそろえて計算しますよね。しかし、揃えようとすると、\[ 3.0453 \times 10^ + 0.0000 \times 10^ \]となってしまい、\( 8.1019 \times 10^ \) の有効数字がなくなってしまい、消えてしまいますね。

計算機上でも 差が極端にあるもの同士 を計算すると、有効数字がなくなってしまい、小さい数が反映されなくなってしまうのです。

(2) 実例1 (1) 単精度 #include #include // 単精度(float)の場合 int main()< float a = 1; float b,sum; for(int i = 0; i >= -10; i -= 1) < b = pow(10,i); sum = a + b; if(a == sum) < printf("b = 10^(%3d) Yes: %.10f\n",i,sum); >else < printf("b = 10^(%3d) No : %.10f\n",i,sum); >> > b = 10^( 0) No : 2.0000000000 b = 10^( -1) No : 1.1000000238 b = 10^( -2) No : 1.0099999905 b = 10^( -3) No : 1.0010000467 b = 10^( -4) No : 1.0001000166 b = 10^( -5) No : 1.0000100136 b = 10^( -6) No : 1.0000009537 b = 10^( -7) No : 1.0000001192 b = 10^( -8) Yes: 1.0000000000 b = 10^( -9) Yes: 1.0000000000 b = 10^(-10) Yes: 1.0000000000 (2) 倍精度

プログラム（単精度で誤差が発生した i = -8からチェックしています）

#include #include // 倍精度の場合 int main()< double a = 1; double b,sum; for(int i = -8; i >= -20; i -= 1) < b = pow(10,i); sum = a + b; if(a == sum) < printf("b = 10^(%3d) Yes: %.20f\n",i,sum); > else < printf("b = 10^(%3d) No : %.20f\n",i,sum); > > > b = 10^( -8) No : 1.00000000999999990000 b = 10^( -9) No : 1.00000000100000010000 b = 10^(-10) No : 1.00000000010000000000 b = 10^(-11) No : 1.00000000001000000000 b = 10^(-12) No : 1.00000000000100010000 b = 10^(-13) No : 1.00000000000009990000 b = 10^(-14) No : 1.00000000000001000000 b = 10^(-15) No : 1.00000000000000110000 b = 10^(-16) Yes: 1.00000000000000000000 b = 10^(-17) Yes: 1.00000000000000000000 b = 10^(-18) Yes: 1.00000000000000000000 b = 10^(-19) Yes: 1.00000000000000000000 b = 10^(-20) Yes: 1.00000000000000000000 (3) 実例2

今回は \( n = 20 \) までとして、昇順 \( n = 1,2,3 \) と加算した場合と、降順 \( n = 20,19,18 \) と加算した場合とで値がどう変わるかを比較します。

\( n = 20 \) までしか計算を行わないため、打ち切り誤差が昇順の場合と降順の場合で発生しますが、両者とも誤差の量は同じ*6になるはずなので今回は考えないことにします。

昇順（大きい値）から足す場合は、\[\frac + \frac + \frac + \cdots\]と計算していくため、足していくにつれて、和に対して加算する \( 1/3^ \) が あまりにも小さくなる ため、情報落ちが発生してしまいます。

一方降順（小さい値）から足した場合は、小さい数から大きい数と計算していくため、足していった和と加算する \( 1/3^ \) の差が大きくないため、情報落ちは発生しません*7。

なので、浮動小数点を用いて和の計算をする際には、大きい数から順番に足すのではなく、 小さい数から順番に足していく と情報落ちを防ぐことができます。

３．丸め誤差

(1) 概要

まずは10進数の場合を考えてみましょう。たとえば、変数aの値が810.364364 という結果だとします。これを小数第3位で四捨五入して \( a = 810.36 \) と計算するとします。

例えば、\( 100a \) の値を小数第2位まで求めるとします。このとき丸めずに計算した場合は、81036.44 と求めることができます。しかし、途中で四捨五入した \( a = 810.36 \) で計算すると、81036.00 となってしまい、0.44の誤差が発生してしまいます。

(2) 具体例 #include #include #define N 20 int main() < float sum_asc = 0; // floatは単精度 float sum_des = 0; for(int i = 1; i printf("ascend : %.20f\n",sum_asc); // 大きい値から加算 printf("descend: %.20f\n",sum_des); // 小さい値から加算 > sum0.1: 9.99999999999998050000 sum0.5: 50.00000000000000000000

ここで、0.1を2進数に直して見ましょう。すると 0.0001 1001 1001 1001… と無限小数になりますね。なので、とある桁で丸められ（切り捨てられ）、誤差が発生してしまいます。

(3) 打ち切り誤差との違い

上にも出てきましたが、打ち切り誤差は、本来無限に続く計算（例：無限等比級数の計算、極限の計算）を有限回で止めたことで発生する誤差です。

例えば、\[ \sum_^ \frac> = \frac \]は、無限回計算した場合に得られる答えですが、計算機上では無限回足せないため、\[ \sum_^ \frac> \]のような有限回の加算となってしまい、誤差が発生してしまいます。

また、\[\lim_ \frac - 3x - 2> - 2x - 8> \] のような極限を考えるときでも、 \( h \) を2に近づける操作を有限の値で打ち切ってしまう*8必要があるため、打ち切り誤差が発生してしまいます。

ただし、打ち切り誤差は、足す回数 \( n \) （極限の場合はよりその値に近づける）などを増やせば増やすほどより無限回の場合（理論値）に近づくので、打ち切り誤差をへらすことができます。

４．2進数の小数表記の際の誤差まとめ

桁落ち（浮動小数点のみ）

ほぼ等しい2つの数の減算で発生する誤差。（例1：20190621.364364 - 20190621.364363）（例2：- 1000.001 + 1000.000）

情報落ち（浮動小数点のみ）丸め誤差（固定・浮動小数点両方）打切り誤差（固定・浮動小数点両方）

５．練習問題

練習1 練習2 練習3 練習4

ア：\( 0.123 \times 10^ + 0.124 \times 10^ \)イ：\( 0.234 \times 10^ - 0.221 \times 10^ \)ウ：\( 0.556 \times 10^ +0.552 \times 10^ \)エ：\( 0.556 \times 10^ - 0.552 \times 10^ \)

練習5