金子邦彦研究室プログラミングC 言語によるアルゴリズムとデータ構造二次方程式の解

二次方程式の解

1 浮動小数点数

少数を含む計算を行うときには、取り扱う少数の精度を 考慮にいれる必要がある.ここではまず、単精度浮動小数点数と倍精度浮動 小数点数の2つについて説明する.

1.1 単精度浮動小数点数

単精度浮動小数点数はc言語ではfloat型として宣言される、符合1ビット、指数 部8ビット、仮数部23ビットの合計32ビットで構成される数値です.
   0  01111111  00000000000000000000000    
 符号  指数8(e)          仮数23(f)
符合、指数、仮数の各部が示す値の組み合わせによって、以下のような数値を表現することができ ます. 単精度浮動小数点数を用いて、円周率3.14159265358979323846...を表現しよう とした場合、表現できる数値は以下のようになる.
  0/1 10000000 10010010000111111011010 = +/- 3.14159250
  0/1 10000000 10010010000111111011011 = +/- 3.14159274
  0/1 10000000 10010010000111111011100 = +/- 3.14159297
つまり、3.1241592までは正しく表現することが出来ますが、それ以上の精度を 期待することはできません。

1.2 倍精度浮動小数点数

倍精度浮動小数点数は、読んで字のごとく単精度浮動小数点数の倍の精度を持ち ます.符合1ビット、指数11ビット、仮数52ビットの合計64ビットで構成され、C 言語ではdouble型として宣言される.
   0  01111111111  000000000000000000000...00    
  符号  指数11(e)          仮数52(f)
単精度浮動小数点数と同様に、符合、指数、仮数の組合わせで以下のような数 値が表現できる. 単精度浮動小数点数の時と同様に円周率を表現しようとした場合、倍精度浮動小数点数を用 いると以下のような数値が表現可能です
  0/1 10000000000 10010010000111111011010..0111 = +/- 3.14159265358979267
  0/1 10000000000 10010010000111111011010..1000 = +/- 3.14159265358979311
  0/1 10000000000 10010010000111111011010..1001 = +/- 3.14159265358979356
つまり、3.141592653589793までは正しく表現することが出来ます.
このように、コンピュータの内部形式では表現できる浮動小数点数の桁数に限 界があり、それ以上の桁を正確に表現することは出来ません。 この、浮動小数点数を表現の可能な有限桁に丸めたために発生する誤差を「丸め誤差」 と呼びます.

2 二次方程式の解

コンピュータによって表現できる浮動小数点数には限界があることを理解したうえ で、プログラムを作成して二次方程式の解を求めてみます.
二次方程式は下のような形式で表現することができる
  Ax2+By+C=0
このとき、二次方程式の解を求めるプログラムの入出力形式は次のようになりま す。 ここではA,B,Cに相当する変数を受け取り、その変数により構成される二次方程 式を解いた結果が出力として表示されるような形式になる. さらに、プログラム内の構成として、倍精度の浮動小数点数を使用する事を定義 します. 二次方程式の解を求める方法としては以下に示すような、おなじみの公式を用いる のが最も一般的です. この公式を実現するプログラムは以下のようになっています
#include<stdio.h>
#include<math.h>

const int BUFFER_SIZE = 100;

main()
{
    double a, b, c, d;
    double x1, x2;
    int n;
    char buf[BUFFER_SIZE], e;

    printf("Please input a b c : ");
    while(1){
    fgets(buf, BUFFER_SIZE, stdin);
    n = sscanf(buf,"%lf %lf %lf %c", &a, &b, &c, &e);
    if(n == 3) {
        break;
    }
    else if(n == -1) {
        continue;
    }
    else {
        printf("Usage:a b c\n");
        return 0;
    }
    }
    d = b * b - 4 * a * c;        /* 判別式の計算 */
    if(d >= 0){                   /* 判別式が正のとき、解を計算 */
    x1 = (-b + sqrt(d)) / (2.0 * a);
    x2 = (-b - sqrt(d)) / (2.0 * a);
    printf("x1 = %27.20e(%27.20e)\n",x1,a*x1*x1+b*x1+c);
    printf("x2 = %27.20e(%27.20e)\n",x2,a*x2*x2+b*x2+c);
    }
    else{                         /* 判別式が負のとき、虚数解と表示 */
    printf("x1 = 虚数解\n");
    printf("x2 = 虚数解\n");
    }

    return 0;
}
これは、プログラムにおける計算中に「桁落ち」と 呼ばれる現象が起っているために起る誤差です.

3 桁落ち

「桁落ち」とは、計算途中に有効数字桁数が減少することです.ここでは、二次 方程式の公式のうち|B|-√(B2-4AC)の部分でこの桁落ちが生じています. A、Cに比べて、Bの値が大きくなるにつれてB2が4ACに比べて極端に大きく なり、その結果、|B|と√(B2-4AC)の値が極めて似通ったものになり、 |B|-√(B2-4AC)の有効数字桁数が極端に減ってしまいます.
A=1、B=-1000000、C=1のとき

 |B|=|-1000000|
    =1000000.000000000000000 (0x412e848000000000)
 √(B2-4AC)=√(1000000000000.0-4.0)
            =999999.999997999984771 (0x412e847fffffbce4)
 2つの差 (|B|-√(B2-4AC))
            =0.000002000015229 (0x3ec0c70000000000)
                                      この部分の情報が欠落
つまりここでは、Bの絶対値がA,Cに比べて大きくなることで桁落ちが生じ、正確な 解が得られなくなってしまっているわけです.

4 プログラムの改良

桁落ちによる精度の悪化を防ぐために、先ほどのプログラムを改良してやる必要 があります.具体的には|B|-√(B2-4AC)で桁落ちが生じているので、こ の式を使わずに計算を行えばOKです.次のように修正します.
    x1 = (fabs(b)+sqrt(d))/(2.0*a); /*fabs()は絶対値を返す関数*/
    x1 = (b<0.0)?x1:-x1; /*bが負ならx1=x1,bが正ならx1=-x1を返す*/
    x2 = c/(a*x1);
修正したプログラムは次のようになる.
#include<stdio.h>
#include<math.h>

const int BUFFER_SIZE = 100;

main()
{
    double a, b, c, d;
    double x1, x2;
    int n;
    char buf[BUFFER_SIZE], e;

    printf("Please input a b c : ");
    while(1){
    fgets(buf, BUFFER_SIZE, stdin);
    n = sscanf(buf,"%lf %lf %lf %c", &a, &b, &c, &e);
    if(n == 3) {
        break;
    }
    else if(n == -1) {
        continue;
    }
    else{
        printf("Usage:a b c\n");
        return 0;
    }
    }
    if(a == 0){
    printf("x = %27.20e\n",-c/b);
    }
    else{
    if(b == 0 && c == 0){
        printf("x1 = %27.20e\n",0.0);
        printf("x2 = %27.20e\n",0.0);
    }
    else{
        d = b * b - 4 * a * c;                   /* 判別式の計算 */
        if(d >= 0){                              /* 判別式が正のとき、解を計算 */
        x1 = (fabs(b) + sqrt(d)) / (2.0 * a);  /* 修正箇所 */
        if ( b >= 0.0 ) {
            x1 = -x1;
        }
        x2 = c / (a * x1);                     /* 修正箇所 */
        printf("x1 = %27.20e(%27.20e)\n",x1,a*x1*x1+b*x1+c);
        printf("x2 = %27.20e(%27.20e)\n",x2,a*x2*x2+b*x2+c);
        }else{                                    /* 判別式が負のとき虚数解と表示 */
        printf("x1 = 虚数解\n");
        printf("x2 = 虚数解\n");
        }
    }
    }

    return 0;
}


   

これをコンパイルし、実行すると、一つめのプログラムで生じていたx2の誤差が無く なっています.つまり、桁落ちが改善され、精度的にも信頼のおけるプロ グラムが完成したということになる. このように、同じ式を表現するプログラムでも、その記述の仕方によって精度に 違いが出ることがあります.実際にプログラムを作成するときは、この精度につ いても十分に考慮するようにしてください。
  1 -2 1
  x1 =  1.00000000000000000000e+00( 0.00000000000000000000e+00)
  x2 =  1.00000000000000000000e+00( 0.00000000000000000000e+00)

  1 -1000 1
  x1 =  9.99998999999000034222e+02( 1.16415321826934814453e-10)
  x2 =  1.00000100000200006431e-03( 0.00000000000000000000e+00)

  1 -1000000 1
  x1 =  9.99999999998999992386e+05( 0.00000000000000000000e+00)
  x2 =  1.00000000000100008890e-06( 0.00000000000000000000e+00)

  1 -1000000000 1
  x1 =  1.00000000000000000000e+09( 1.00000000000000000000e+00)
  x2 =  1.00000000000000006228e-09( 0.00000000000000000000e+00)