float加算時の誤差を調べてみる

Tweet

Check

物理学科で学んでいる1号から「0.1ずつたして4.0で止まるループのはずなのに、4.1までやっちゃうんだけどどうして？」という相談を受けた。
Cで書くとこんな感じ。実際はいろいろ計算が入っているが、問題はそこじゃないので削ってある。

float x=0;
float step=0.1;
while(x<4.0) {
  x += step;
  printf("x=%5.3f\n", x);
}

x=3.9でwhile(x<4.0)の条件をパスして、x=4.000を表示した後、x=4.0でwhile()の条件をパスできなくて終わるはずなんだが、パスしてx=4.100を表示してしまう。

floatだし、誤差の積み上げでx=4.000に見えるけど、実際はx=3.9999xxxなんでしょってのはまぁ推測できる。
実際、printf()の書式指定を%5.3fから%8.6fにすると、x=4.000はx=3.999998となり、4.0より小さいことがわかる。

こういうことが起きるというのを知って、こういう誤差で悩まされないためにどう対処するのかってのを問いている課題なんじゃなかろうか。

とりあえず、
・floatからdoubleにする
・0.1刻みならxを小数点第2位で四捨五入する
・xを整数で扱ってループを回し、計算処理にはxを実数値に変換して使う。(例:1刻みで計算にはx/10.0を使う)
の3つを回避策として提示した。

しかし、
・0.2刻みだと症状が出ない
・0.1はfloatで扱うと0.1よりちょっと大きいのに加算していくと4.0になるはずが3.999998とか意味不明
などと言われて説明に困ってしまった。

まぁ、浮動小数点数の仮数部の桁数の制限があるから桁落ちとかあるんだよとはなんとなくわかるが、ちゃんと説明はできない。

仕方ないので、0.1ずつ足していったときに実際どういうふうに扱われているのかがわかるように表示するプログラムを書いてみた。

出力をざっとみて、2進数の足し算をしてみたりしてなんとなく、桁落ちした部分の切り上げ、切り捨てが発生して0.1を足したいのに0.1より大きかったり、小さかったりしているのがわかる。
また、浮動小数点数の指数部が増えるとき(桁がずれるとき)に差分も変化している。

0.1より実際は大きい数値を加算しているのに40回加算して4.0より小さいのは、差分が0.1より小さい回が多くてそうなったようだ。
実際41回目からは0.1より差分が大きくなっているので、そのうち加算結果は回数/10を超えるかもしれない。

0.2刻みだと違う回数でこの問題が発生するはずだ。

ちなみに、0.25刻みだとこの問題は発生しない。2進数で表すと0.01ちょうどで誤差なく2進数に変換できるから。

#include <stdio.h>
#include <stdlib.h>

//floatの値を加算していくとどうなるのかを検証する
/*
 #  TIME(x)                          xの符号,指数部,仮数部               xの差分
 0, 0.0000000000000000000000000000 /   0.0
 1, 0.1000000014901161193847656250 / + [-004] b)1.10011001100110011001 / 0.1000000014901161193847656250
 2, 0.2000000029802322387695312500 / + [-003] b)1.10011001100110011001 / 0.1000000014901161193847656250
 3, 0.3000000119209289550781250000 / + [-002] b)1.00110011001100110011 / 0.1000000089406967163085937500
*/

/* 参考資料
 * 【C言語】浮動小数点数における「数値⇔内部データ（符号部・指数部・仮数部）」の変換  https://daeudaeu.com/c_float/
 * 浮動小数点数の内部表現を取得してみよう  https://qiita.com/nia_tn1012/items/d26f0fc993895a09b30b
 * 浮動小数点数内部表現シミュレーター  https://tools.m-bsys.com/calculators/ieee754.php
 */


char BitStr[33];  //32bit +終端文字  //floatの内部表現を二進数(文字列)化したものを格納
char FloatInternal[256];             //floatの内部表現を読みやすい形式に変換したものを格納

//指数部と仮数部が0なら0.0なのでtrue(1)を返す。それ以外はfalse(0)。
int isZeroFloat(char *str) {
	int isZero = 0;
	for(int cnt=1; cnt<32; cnt++) { //str[0]は符号なのでスキップ if('0' == str[cnt]) { isZero = 1; } else { isZero = 0; //'0'以外を見つけたら0.0ではないので終了 break; } } return isZero; } //floatの各ビットを'0','1'で表して文字列配列に格納する void Float2BitStr(float pReal) { union { float f; int i; } fandi; int cnt; fandi.f = pReal; for(cnt=31;cnt>=0;cnt--){ //上の桁から処理する。floatは32ビット。 
		BitStr[31-cnt] = ((fandi.i >> cnt) & 1)? '1':'0';
		//printf("%2d - %1d\n", 31-cnt, (fandi.i >> cnt) & 1);
	}
	BitStr[32] = '\0';  //終端させる
	//printf("%16.14f = %s\n", pReal, BitStr);

	
	return;
}


//floatを二進数化(文字列)化したものを読みやすく変換する
void Str2FloatInternal(char *str) {
	if(isZeroFloat(str)) { // 0.0の場合は、'  0.0'を表示して終わり。
		FloatInternal[0] = ' ';
		FloatInternal[1] = ' ';
		FloatInternal[2] = '0';
		FloatInternal[3] = '.';
		FloatInternal[4] = '0';
		FloatInternal[5] = '\0';
		return;
	}

	FloatInternal[0] = ('1'==str[0]) ? '-':'+'; //正または零は符号表示なし

	int base=1;
	int nbits=8;
	int tmp=0;
	for(int cnt=0; cnt<nbits; cnt++){ //指数部を取り出して整数化(10進数化)
		// printf("%2d - %c %1d\n", cnt, str[base+cnt], nbits-cnt-1);
		if('1' == str[base+cnt]) {
			tmp += 1<<(nbits-cnt-1);
		}
	}
	tmp -= 127;
	FloatInternal[1] = ' ';
	//指数部(10進数)を文字列化
	FloatInternal[2] = '[';
	FloatInternal[3] = (tmp<0) ? '-':' ';  //正または零は符号表示なし
	FloatInternal[4] = '0' + (abs(tmp)/100);
	FloatInternal[5] = '0' + (abs(tmp)/10 % 10);
	FloatInternal[6] = '0' + (abs(tmp) % 10);
	FloatInternal[7] = ']';

	FloatInternal[8] = ' ';

	//仮数部を文字列化(2進数)
	FloatInternal[9] = 'b';
	FloatInternal[10] = ')';
	FloatInternal[11] = '1';
	FloatInternal[12] = '.';

	for(int cnt=9, idx=13; cnt<32; cnt++, idx++) {
		FloatInternal[idx] = str[cnt];
	}

	FloatInternal[33] = '\0'; //終端
		
	return;
}


void main() {
	float TIME_INIT=0.0;
	float TIME_LAST=4.0;
	float STEP = 0.1;
	
	//Float2BitStr(STEP);
	//Str2FloatInternal(BitStr);
	//printf("%20.18f = %s / %s\n", STEP, BitStr, FloatInternal);

	int i = 0;
	float x = TIME_INIT;
	float xpre = x;

/*
	Float2BitStr(x);
	Str2FloatInternal(BitStr);
	printf("%2d, %32.30f / %s\n", i, x, FloatInternal);
*/

	while(x<TIME_LAST) {
		Float2BitStr(x);
		Str2FloatInternal(BitStr);
		printf("%2d, %32.30f / %-33s / %32.30f\n", i, x, FloatInternal, x-xpre);
		xpre = x;
		x += STEP;
		i++;
	}

}