二次方程式の解

1 浮動小数点数

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

1.1 単精度浮動小数点数

単精度浮動小数点数はc言語ではfloat型として宣言される、符合1ビット、指数部8ビット、仮数部23ビットの合計32ビットで構成される数値です．

   0  01111111  00000000000000000000000    
 符号  指数8(e)          仮数23(f)

符合、指数、仮数の各部が示す値の組み合わせによって、以下のような数値を表現することができます．

指数=0 かつ仮数=0 … 符合付き0(±0)
指数=0 かつ仮数≠0 … 不正規化数
0<指数<255 … 正規化数
指数=255かつ仮数=0 … 符合付き無限大(±∞)
指数=255かつ仮数≠0 … 非数値(0/0の結果など、符合は意味を持たない)

単精度浮動小数点数を用いて、円周率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 かつ仮数=0 … 符合付き0(±0)
指数=0 かつ仮数≠0 … 不正規化数
0<指数<2047 …正規化数
指数=2047かつ仮数=0 … 符合付き無限大(±∞)
指数=2047かつ仮数≠0 … 非数値(0/0の結果など、符合は意味を持たない)

単精度浮動小数点数の時と同様に円周率を表現しようとした場合、倍精度浮動小数点数を用いると以下のような数値が表現可能です

  0/1 10000000000 10010010000111111011010..0111 = +/- 3.14159265358979267
  0/1 10000000000 10010010000111111011010..1000 = +/- 3.14159265358979311
  0/1 10000000000 10010010000111111011010..1001 = +/- 3.14159265358979356

つまり、3.141592653589793までは正しく表現することが出来ます．
このように、コンピュータの内部形式では表現できる浮動小数点数の桁数に限界があり、それ以上の桁を正確に表現することは出来ません。この、浮動小数点数を表現の可能な有限桁に丸めたために発生する誤差を「丸め誤差」と呼びます．

2 二次方程式の解

コンピュータによって表現できる浮動小数点数には限界があることを理解したうえで、プログラムを作成して二次方程式の解を求めてみます．
二次方程式は下のような形式で表現することができる

  Ax²+By+C=0

このとき、二次方程式の解を求めるプログラムの入出力形式は次のようになります。

標準入力 … A,B,C
標準出力 … 方程式の解 x1,x2

ここではA,B,Cに相当する変数を受け取り、その変数により構成される二次方程式を解いた結果が出力として表示されるような形式になる．さらに、プログラム内の構成として、倍精度の浮動小数点数を使用する事を定義します．

Aを記憶する変数 … double型 a
Bを記憶する変数 … double型 b
Cを記憶する変数 … double型 c
判別式を記憶する変数 … double型 d
解を記憶する変数 … double型 x1,x2

二次方程式の解を求める方法としては以下に示すような、おなじみの公式を用いるのが最も一般的です．

x1=(-B+√(B²-4AC)/2A
x2=(-B-√(B²-4AC)/2A

この公式を実現するプログラムは以下のようになっています

#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|-√(B²-4AC)の部分でこの桁落ちが生じています． A、Cに比べて、Bの値が大きくなるにつれてB²が4ACに比べて極端に大きくなり、その結果、|B|と√(B²-4AC)の値が極めて似通ったものになり、 |B|-√(B²-4AC)の有効数字桁数が極端に減ってしまいます．

A=1、B=-1000000、C=1のとき

 |B|=|-1000000|
    =1000000.000000000000000 (0x412e848000000000)
 √(B²-4AC)=√(1000000000000.0-4.0)
            =999999.999997999984771 (0x412e847fffffbce4)
 2つの差 (|B|-√(B²-4AC))
            =0.000002000015229 (0x3ec0c70000000000)
                                      この部分の情報が欠落

つまりここでは、Bの絶対値がA,Cに比べて大きくなることで桁落ちが生じ、正確な解が得られなくなってしまっているわけです．

4 プログラムの改良

桁落ちによる精度の悪化を防ぐために、先ほどのプログラムを改良してやる必要があります．具体的には|B|-√(B²-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)