この前計算機室にいたら，学生が「long ではなくて float を使うのかな」などと話しておりました。
# float より double を使ってくれ...
コンピュータの内部で使用される数は，2進数で表される。 以下のプログラムを考えてみよう。
#include <stdio.h> #include <math.h> void put_float(char *string, float z); #define N 100 int main(void) { float x, sum; int i; sum = 0.0; x = 1.0/N; /* (1) */ for (i = 0; i < N; i++) { /* (2) */ sum += x; } put_float("sum = ", sum); if (sum == 1.0) { printf("合計は１になりました\n"); } else { printf("合計は１になりませんでした\n"); } return 0; } void put_float(char *string, float z) { printf("%s%.16g\n", string, z); }
floatというのは，4バイトを使って近似的に数値を表すもので，10進数では 有効桁数は6桁程度しかない。sum = 0.9999993443489075 合計は１になりませんでした
float x, sum;をdouble x, sum;にしても根本的な解決にはならない。
このプログラムが正しい答えを出さないのは，“コンピュータ内では実数値は 近似値でしかない”ということが原因である。sum = 1.000000000000001 合計は１になりませんでした
#define N 100の代わりに
#define N 128のように，Nを2のべき乗にすると
のように期待される答えを出す。1/128というのは，10進数では0.0078125であるが， 2進数でも近似値ではなく正確な数値として格納されるからである。ちなみに，整数値 (小数部分を持たない数)は全て正確な数値として格納される。sum = 1 合計は１になりました
大きさの同程度の数値の差を取ると，有効桁は極端に減少する。例えば，
123.456789−123.456788 = 0.000001となる。2数の有効桁は9桁あるのに，
答えの有効桁は1桁になってしまう。
以下に示したプログラムは二次方程式の解を求めるものであるが，誰でも知っている
“解の公式”を用いて(1), (2) のように計算すると不都合が生ずる場合がある。
このとき，a = 1 b = 1.0000000000000001 c = 0.0000000000000001
となり，disc = 0.9999999999999996 −b+disc = −4.440892098500626e−16 −b−disc = −2
−b+discの計算で桁落ちが生ずる。その結果は以下のようになり，
小さい方の解の精度は不十分なものである。
桁落ちを避けるために，(3),(4)のようにして，大きい方の解をまず求め，次いで 解と係数の公式から小さい方の解を求めるようにするのが正しいやり方である。 (3)あるいは(5)の計算では，桁落ちが生じない。解の公式による解 x1 = −1.110223024625157e−16 x2 = −0.9999999999999999
解の公式と解と係数の関係による解 x1 = −1e−16 x2 = −0.9999999999999999
#include <stdio.h> #include <math.h> double get_double(char *string); void put_double(char *string, double z); int main(void) { double a, b, c, disc, x1, x2; a = get_double("a = "); b = get_double("b = "); c = get_double("c = "); disc = b*b-4.0*a*c; if (disc > 0.0) { printf("解の公式による解\n"); x1 = (-b+sqrt(disc))/(2.0*a); /* (1) */ x2 = (-b-sqrt(disc))/(2.0*a); /* (2) */ put_double("x1 = ", x1); put_double("x2 = ", x2); printf("解の公式と解と係数の関係による解\n"); if (b > 0.0) { x2 = (-b-sqrt(disc))/(2.0*a); /* (3) */ x1 = (c/a)/x2; /* (4) */ put_double("x1 = ", x1); put_double("x2 = ", x2); } else { x1 = (-b+sqrt(disc))/(2.0*a); /* (5) */ x2 = (c/a)/x1; /* (6) */ put_double("x1 = ", x1); put_double("x2 = ", x2); } return 0; } else { printf("実数解はありません\n"); return 1; } } double get_double(char *string) { char line[128]; printf(string); gets(line); return atof(line); } void put_double(char *string, double z) { printf("%s%.16g\n", string, z); }
コンピュータの内部では，有限の桁数を用いて近似計算が行われている。
我々が筆算を有効桁3桁で切り捨て計算を行うときを考えてみよう。
12.3+0.456=12.756であるが，結果は12.7となり，0.056は蒸発してしまう。
もし，四捨五入により計算したときは，結果は12.8となり，0.044だけ積み過ぎてしまう。
1 2 . 3
+ 0 . 4 5 6
----------------
1 2 . 7
n
平均値 =ΣXi/n ・・・・・・・・・・・・・・・・（1）
i=1
n
分散 = Σ(Xi−平均値)2/(n−1) ・・・・・・・・（2）
i=1
n
二乗和 = ΣXi2 ・・・・・・・・・・・・・・・・（3）
i=1
分散 = （二乗和−n・平均値2）/（n−1） ・・・・・（4）
平均値と分散の計算結果
呼出し(a) 0から100まで prog1: mean = 49.5 variance = 841.6666870117188
呼出し(b) 10000から10100まで prog1: mean = 10049.5 variance = 858.8181762695312
呼出し(c) 0から100まで prog2: mean = 49.5 variance = 841.6666870117188
呼出し(d) 10000から10100まで prog2: mean = 10049.5 variance = 841.6666870117188
841.6666666666666であり，
単精度での答えは当然ではあるが有効桁数は少ない。
#include <stdio.h> #include <math.h> float get_float(char *string); void put_float(char *string, float z); void prog1(long base); void prog2(long base); #define N 100 int main(void) { prog1(0); /* 呼出し(a) */ prog1(10000L); /* 呼出し(b) */ prog2(0); /* 呼出し(c) */ prog2(10000L); /* 呼出し(d) */ return 0; } void prog1(long base) { float mean, variance; long i; mean = variance = 0; for (i = base; i < base+N; i++) { mean += i; variance += (float)i*i; /* 二乗和を計算する */ } mean / = N; variance = (variance-mean*mean*N)/(N-1); /* 危険な方法 */ put_float("prog1: mean = ", mean); put_float("prog1: variance = ", variance); } void prog2(long base) { float mean, variance; long i; mean = variance = 0; for (i = base; i < base+N; i++) { mean += i; } mean / = N; for (i = base; i < base+N; i++) { /* 安全な方法 */ variance += pow(i-mean, 2.0); } variance / = N-1; put_float("prog2: mean = ", mean); put_float("prog2: variance = ", variance); } void put_float(char *string, float z) { printf("%s%.16g\n", string, z); }
以上からの教訓 ● float は使うな，double を使え ● いい加減に公式を使うな，よいアルゴリズムを使え
直前のページへ戻る E-mail to Shigenobu AOKI