C# で二次方程式の解の桁落ちを見る

Tech
Share

二次方程式の解を公式で解くことで、C# で桁落ちを試します。

桁落ち

cancellation of significant digits

桁落ちとは、浮動小数点演算で、計算結果が 0 に極端に近くなる加減算を行ったときに、有効数字の桁数が極端に少なくなる現象。

（中略）

桁落ち自体による問題はコンピュータとは無関係に発生するが、コンピュータ上での桁落ちは、計算途中の値が分からない・結果の桁数が常に一定なので気づきにくいという特徴がある。

非常に有名な二次方程式の解の公式を用います。 ax2+bx+c=0  (a≠0)ax^{2}+bx+c=0\;(a\neq0) に対して、

x=−b±b2−4ac2ax=\frac{-b\pm\sqrt{b^{2}-4ac}}{2a}

が成り立ちます。ここで、 ∣4ac∣≪1\left|4ac\right|\ll1 のとき、 b2−4ac≃∣b∣\sqrt{b^{2}-4ac}\simeq\left|b\right| となり、解のどちらかでは同じくらいの値同士の引き算、つまり桁落ちが発生してしまいます。 これを回避するためには、桁落ちが起きない方の解からもう片方の解を求める処理を行う必要があります。以下の式を用います。

x1⋅x2=ca∴x1=cax2,  x2=cax1x_{1}\cdot x_{2}=\frac{c}{a}\quad\therefore x_{1}=\frac{c}{ax_{2}},\;x_{2}=\frac{c}{ax_{1}}

今回は float だけ見てみます。ただし、平方根の計算も入るため、桁落ちを回避したとしても誤差は出てきてしまいます。ソースコード等はこちらから

まずは計算処理部分から。

12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
static Tuple<float, float> NormalSolver(float a, float b, float c)
{
    float det = (float)Math.Sqrt(b * b - 4 * a * c);
    float x_1 = (-b + det) / 2.0f / a;
    float x_2 = (-b - det) / 2.0f / a;
    return Tuple.Create(x_1, x_2);
}

static Tuple<float, float> CaredSolver(float a, float b, float c)
{
    float det = (float)Math.Sqrt(b * b - 4 * a * c);
    float x_1, x_2;
    if (b > 0)
    {
        x_2 = (-b - det) / 2.0f / a;
        x_1 = c / (a * x_2);
    }
    else
    {
        x_1 = (-b + det) / 2.0f / a;
        x_2 = c / (a * x_1);
    }
    return Tuple.Create(x_1, x_2);
}

通常は NormalSolver のように計算すると思います。 CaredSolver では b の正負を確かめて、引き算による桁落ちを防ぎます。

計算結果は下のプログラムを使って確認しました。

37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
static void CheckCancellation()
{
    float a = 1.0f, b = -1000.0f, c = 1.0f; // x^2 - 1000x + c = 0

    var normal = NormalSolver(a, b, c);
    var cared = CaredSolver(a, b, c);

    Console.WriteLine("f(x) = x^2 - 1000x + c = 0");

    Console.WriteLine($"normal: x_1= {normal.Item1:F8}, x_2= {normal.Item2:F8}");
    Console.WriteLine($"f(x_2) -> {(a * normal.Item2 * normal.Item2 + b * normal.Item2 + c):F8}");

    Console.WriteLine($"cared : x_1= {cared.Item1:F8}, x_2= {cared.Item2:F8}");
    Console.WriteLine($"f(x_2) -> {(a * cared.Item2 * cared.Item2 + b * cared.Item2 + c):F8}");
}

計算結果は次の通りです。

1
2
3
4
5
f(x) = x^2 - 1000x + c = 0
normal: x_1= 999.99900000, x_2= 0.00100708
f(x_2) -> -0.00707901
cared : x_1= 999.99900000, x_2= 0.00100000
f(x_2) -> 0.00000006

対策の有無で、 x_2 の値に違いが出てきています。 f(x_2) を計算してみると、桁落ちの対策をした場合は0との差が小さくなる結果となりました。 x_1 での結果を見てはいけない…。