二次方程式の解
1 浮動小数点数
少数を含む計算を行うときには、取り扱う少数の精度を 考慮にいれる必要がある.ここではまず、単精度浮動小数点数と倍精度浮動 小数点数の2つについて説明する.1.1 単精度浮動小数点数
単精度浮動小数点数はc言語ではfloat型として宣言される、符合1ビット、指数 部8ビット、仮数部23ビットの合計32ビットで構成される数値です.0 01111111 00000000000000000000000 符号 指数8(e) 仮数23(f)符合、指数、仮数の各部が示す値の組み合わせによって、以下のような数値を表現することができ ます.
- 指数=0 かつ 仮数=0 … 符合付き0(±0)
- 指数=0 かつ 仮数≠0 … 不正規化数
-
1.40129846*10-45から1.1754921*10-38までの数、つまり、物凄
く小さい数。
(-1)符合*2-127*0.fにより求められる - 0<指数<255 … 正規化数
-
1.1749435*10-38以上の有限の数値。
(-1)符合*2指数-127*1.f により求まる。 - 指数=255かつ仮数=0 … 符合付き無限大(±∞)
- 指数=255かつ仮数≠0 … 非数値(0/0の結果など、符合は意味 を持たない)
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 … 不正規化数
-
4.9406564584124654*10-324から2.2250738585072008*10-308ま
での数、つまり、とんでもなく小さい数。
(-1)符合*2-1022*0.fにより求められる - 0<指数<2047 …正規化数
-
2.2250738585072013*10-308以上の有限の数値。
(-1)符合*2指数-1023*1.f により求まる。 - 指数=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 二次方程式の解
コンピュータによって表現できる浮動小数点数には限界があることを理解したうえ で、プログラムを作成して二次方程式の解を求めてみます.二次方程式は下のような形式で表現することができる
Ax2+By+C=0このとき、二次方程式の解を求めるプログラムの入出力形式は次のようになりま す。
- 標準入力 … A,B,C
- 標準出力 … 方程式の解 x1,x2
- Aを記憶する変数 … double型 a
- Bを記憶する変数 … double型 b
- Cを記憶する変数 … double型 c
- 判別式を記憶する変数 … double型 d
- 解を記憶する変数 … double型 x1,x2
- x1=(-B+√(B2-4AC)/2A
- x2=(-B-√(B2-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|-√(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; } |
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)