浮動小数点数について

雑記 小ネタ

はじめに

この記事はCompetitive Programming Advent Calendar12月10日分の記事です。
誤差が原因でWAになる問題が嫌いな人は決して少なくないと思います。「こんなのは本質的なことではない」という人もいるでしょう。誤差を回避することが本質となるような問題がコンテストに出てもきっと人気は出ないでしょう。自分も以前は誤差死の起きうる問題は嫌いでしたが、最近はそれを回避するのも面白さだと感じるになってきました。
しかし面白いとかつまらないとか好きとか嫌いとか以前に非ACは悔しいものです。なので今回は、誤差死などを回避するために地味なネタである浮動小数点数について書こうと思います。

A+B problem

とりあえず、いろんなオンラインジャッジで最初の1問目として見ることの多いA+B problemを解いてみましょう。よくあるのは整数の和を出力させる問題ですが、今回は小数での計算をすることにします。

  • Input
    • 入力は小数で表されるa, b(0
  • Output
    • a+bの値を四捨五入して小数点以下2桁まで出力せよ。
  • Sample Input 1
    • 3.215 4.2537
  • Sample output 1
    • 7.47
  • Sample Input 2
    • 3.0 5.0
  • Sample Output 2
    • 8.00
  • Sample Input 3
    • 0.0002 0.001
  • Sample Output 3
    • 0.00

何の変哲もない問題に見えます。TopCoderやCodeForcesをメインにしている人には、「小数点以下2桁まで出力」というのに馴染みが薄いかもしれませんが、出力の誤差評価をするのは案外面倒な作業なので、オンラインジャッジに収録されている問題はこのような出力形式になっていることが多いです(AOJは例外的に誤差を許容する問題が多い)。さて、問題に特にトラップがあるようには思いません(無理やり0.004999...9999 + 0.0000...001みたいな意地悪な入力は作れますが、入力はあまり長くないものとしておきましょう。本当はこれも甘い？)。この解答プログラムをC++で書くとこんな感じでしょうか。

#include <cstdio>
using namespace std;

int main(){
	double a, b;
	scanf("%lf %lf", &a, &b);
	printf("%.2f\n", a + b);
	return 0;
}

このプログラムでサンプルは全部通ります。でも、わざわざこんなことを書いているのはこのプログラムを撃墜する入力が存在するからです。少し考えてみましょう。
オンラインジャッジで問題を多く解いている人はすぐ気づくでしょうし、TopCoderやCodeForcesの方をメインでやっている人には気づきにくいかもしれません。
では正解発表です。このプログラムに次の入力を与えてみましょう。

0.01 0.005

この答えを素直に計算すると、0.01+0.005=0.015なので、小数第3位を四捨五入して答えは0.02になってほしいはずです。でも出力はこうなります。

0.01

悲しいことに、このプログラムは正解である0.02を出力してくれません。この理由は、「1.0/3.0*3.0が1.0になりません」のようなFAQを知っている人にはすぐに分かると思います。
普段プログラミングコンテストで小数を扱うときはdouble(倍精度浮動小数点数)型で処理することが多いと思います。ここでIEEE 754やら、マシンイプシロンやらの真面目な話をするつもりはありません。知っておくと便利なことをいくつか挙げてみます。

  1. 64ビットである。
  2. 保持しておける値は、実際の値の近似値であることが多い。これが原因で生じる誤差を丸め誤差という。
  3. intの範囲に収まる整数は誤差なしに保持しておける。また、intの範囲で割りきれる割り算はやはり誤差無しでできるらしい(仕様で決まっているのかどうかは知らない)。
  4. 計算速度は 速い < intの加算減算 <= intの乗算 , doubleの加算減算 <= doubleの乗算除算 < intの除算剰余 < 遅い、といった感じになる。
  5. 大きい数に小さい数を足しても値が変わらないことがある(丸め誤差の一例、複数の和をとるときは小さい方から和をとる方が誤差が少ない)。
  6. 0割などの通常でない操作を施すとNaN(Not a Number)という状態になり、しかもそれは伝播する(Runtime Errorは起きない。また、NaNにも複数の種類がある)。
  7. 値が非常に小さく(0に近く)なると非正規化数という計算が非常に遅い状態になる。

さっきのA+B problemの解答が0.01を返したのは2が原因です。double型は0.015(or 0.01, 0.005)という値を正確に保持することができずに、0.014999..という近似値で保持してしまったので、四捨五入した結果0.01という値を出力してしまったのでした。
原因は分かりました。ではこれをどうやって修正しましょうか。JavaやCsharp、VBは固定小数点数を扱うことができますが、それは最後の手段にしたいものです。手頃な手段としては、0.01499...という値がちゃんと0.01500..になるように、適当な小さな数を足してやる方法が考えられます。このような小さな数はEPSという変数名が付けられることが多いようです(数学で微小量を表すことが多いεに由来しています)。プログラムを修正してみます。

#include <cstdio>
using namespace std;

const double EPS = ??;

int main(){
	double a, b;
	scanf("%lf %lf", &a, &b);
	printf("%.2f\n", a + b + EPS);
	return 0;
}

ここでEPSの値を適切に決めるとこのプログラムは正答を返します。EPSが大きいと繰り上げしてはいけないのに繰り上げしてしまうかもしれません。できるだけ小さい方がいいです。じゃあどれくらい小さくすれば良いのでしょうか。
doubleの性質5に関係する話ですが、あまりに小さいと加えても無視されます。大体、元の値の2^(-53)=10^(-15) ~ 10^(-16)倍くらい小さいと無視されてしまう可能性があります。a+bは大体10^2のオーダーであることを考えると、10^(-12)位までなら小さくしてよさそうです。実際は、そんな細かいことを考えなくても10^(-8)とか10^(-9)とかで固定して使って大丈夫なことが多いです。

二分探索とEPS(1)

TopCoderとかで実数の二分探索がでた場合の撃墜ポイントとしてよくあるのが、二分探索の終了条件を

while(high - low > EPS){
        ...
}

としている人です。これもdoubleが大きい数同士の引き算が苦手なことが原因で、答えが大きくなるような入力を与えると終了せずに無限ループしてTLEしてしまいます。EPSがどのくらいのときに答えがどれくらいの大きさになるようにすればTLEを発生させることができるかも簡単に見積もることができます。
これを防ぐ方法は、(lowとhighの初期値にもよるけれど)100回くらいループさせたり、ちゃんと相対誤差も考慮して

while( min(absolute_error(high, low), relative_error(high, low)) > EPS ){
        ...
}

とするのが考えられますが、後者で処理している人を実際に見たことはありません。

幾何とEPS

幾何の問題では平方根や直線の交点などを求める必要がある場合が多く、整数で処理しきれない場合があります。そのようなときもEPSを使って等式や不等式のチェックをしなければなりません。これを避けるのは本質的に難しく、ロバストな計算幾何学という一研究分野が存在するほどです(競技プログラミング以外の部分でもその需要は大きい)。
既に書きましたが元の値の2^(-53)=10^(-15) ~ 10^(-16)倍くらい小さい値は無視されてしまう可能性があります。それを考慮してEPSを使わなければなりません(等号の判定に相対誤差を使ったり)。
例えば、

const double EPS = 1e-9;

...

... if(x*x > y*y + EPS){ ...

のようなコードはよく見かけます。しかし、xやyが10^4のオーダーだったりするとこのEPSはもう意味をなしません。
これを防ぐためには、比較する数が大きいときは相対誤差を使うことが考えられますが、やはり実際にそうしてる人はあまり見たことがありません。

二分探索とEPS(2)

浮動小数点数の比較には常にEPSを使えばいいというものではありません。EPSを使うと失敗してしまうような場合もあります。それは二分探索の判定部分です。そこで安易にEPSを使ってしまうと、出てくる解が真の解より離れてしまうことがあります。許容誤差が10^(-9)だったりするときにEPS=10^(-9)とかでそれをやってしまうと最後の桁が合わないことが起きるので注意です。
例として、正の数xが与えられた時にx^3をわざわざ二分探索を使って求めるプログラムを見てみましょう。

#include <cmath>
#include <cstdio>
using namespace std;

const double EPS = 1e-10;

int main(){

	double x = 2.0;

	double low = 0.0, high = 100.0;
	for(int loop=0; loop<100; loop++){
		double mid = (low + high) / 2.0;
		if( cbrt(mid) < x - EPS){
			low = mid;
		}
		else{
			high = mid;
		}
	}

	printf("%.9f\n", low);
	return 0;
}

これを実行すると7.999999999となって真の解から誤差がでます。この例だけ見たら何をバカなことをやってるんだと思うかもしれませんが、幾何が絡んだときには意外とこれと同類のミスをやってしまいがちです。

三分探索と誤差

以前に記事を書いたのでそれを参考に。読むのが面倒な人向けに一言でまとめると、三分探索は最大値そのものは非常によい精度で求めることができるけど、最大値を与える値はよい精度で求められないという話です。逆関数の微分係数が大きくて摂動がうんぬんかんぬん。

対数と指数

doubleを使うメリットの一つに、intに比べてオーバーフローしにくい、というのはあると思います。しかし完璧に処理しているわけではありません。
(たくさんの数の積)/(たくさんの数の積)のような計算をするときは、先に被除数だけを計算しようとするとオーバーフローしてしまう可能性があります。このようなときは、exp( (たくさんの数の対数の和) - (たくさんの数の対数の和) )で計算することによってオーバーフローを防ぐことができます。実はこの方法も精度の問題はあるのですが、オーバーフローに比べればマシです。

数学関数

標準のライブラリで提供されている数学関数にも誤差は入ってしまいます。例えば、与えられた正の整数xが立方数かどうかO(1)で判定したいときは

int cb(int x){
        return x*x*x;
}

bool is_cb(int x){
        return x>=0 && cb((int)cbrt(x)) == x;
}

と書く方法が考えられます。しかしこれはx<=cbrt(x*x*x)

bool is_cb(int x){
        return x>=0 && cb((int)(cbrt(x) + 0.5)) == x;
}

とかで十分正しく動きます。手元の環境では[0..10^18]に含まれる全ての立方数に対してtrueを返していました。
数学関数がどれくらい正確な値を返すのか一度調べておくとよいかもしれません。例えばガンマ関数を使って階乗を手軽に求めたいとき誤差が入るかどうかとか。

非正規化数

浮動小数点数の性質7で挙げたように、double型のある値xが0

#include <vector>
using namespace std;

struct FoxListeningToMusic {
	vector<double> getProbabilities(vector<int> length, int T) {
		const int N = length.size();
		const double p = 1.0/N;
		vector< vector<double> > dp(T+2, vector<double>(N,0.0));
		dp[0][0] = 1.0;
		for(int i=0; i<=T; i++){
			for(int j=0; j<N; j++){
				for(int k=0; k<N; k++){
				        int nxt = min(T+1, i+length[k]);
					dp[nxt][k] += dp[i][j]*p;
				}
			}
		}
		return dp[T+1];
	}
};

このプログラムの計算量がO(N^2*T)であることはすぐに分かります。Nの最大値は50、Tの最大値は80000です。50*50*80000=200000000なので危険ではありますが、計算自体は単純そうなので間に合うかもしれません。
実際、このプログラムに最悪ケースと思われるlength={1,1,...,1}(長さ50)、T=80000を与えると、手元の環境ではO2オプションの下で約0.98secで答を返してくれました。2secあれば間に合いそうです。
しかしこの入力をちょっといじってlength={1,1,...,1,80000}(長さ50)、T=80000としてみましょう。一見これで実行時間が変わるようには思えませんが、実際には劇的に変わります。同じ環境のもとで答えを返すまで約37.51secかかりました。10倍どころではありません。
細かい説明は置いといて、先のプログラムにちょっとだけ手を加えてみましょう。

#include <vector>
using namespace std;

struct FoxListeningToMusic {
	vector<double> getProbabilities(vector<int> length, int T) {
		const int N = length.size();
		const double p = 1.0/N;
		vector< vector<double> > dp(T+2, vector<double>(N,0.0));
		dp[0][0] = 1.0;
		for(int i=0; i<=T; i++){
			for(int j=0; j<N; j++){
                                //追加部分
				if(dp[i][j] < 1e-100){
					dp[i][j] = 0.0;
				}
				for(int k=0; k<N; k++){
				        int nxt = min(T+1, i+length[k]);
					dp[nxt][k] += dp[i][j]*p;
				}
			}
		}
		return dp[T+1];
	}
};

if文を加えました。これは枝刈りにはなっていないので計算量は変わっていないことが分かると思います。このif文を付け加えた入力にさっきと同じ入力を与えると、今度は約1.02secで答えを返してくれました。一体どういうことなんでしょうか。
ネタばらしをすると、あのTLEを誘発した入力は途中からDPテーブルの中の値が極端に小さくなります。そしてある閾値を越えて小さくなると、非正規化数という状態になって計算が極端に遅い状態になってしまったのでした。下のプログラムでは、値が答えに大きく影響しない程度に小さくなったら強制的に0に落としてしまうことで非正規化数になってしまうことを防いでいます。
なので、確率のDPなどで値が極端に小さくなってしまうときには非正規化数にならないように注意しましょう。ちなみにこの問題にはちゃんとO(N*T)解があるので、それだとこんな対策をせずとも2sec以内に答えを返します。

脱・浮動小数点数

これまで見たように浮動小数点数の世界は狂気に満ちています。なにせ結合法則すら成り立たない世界なので想像しづらい現象が度々起こってしまいます。なので、浮動小数点数を扱わずにいろんなことを済ませてみることを考えます。
有理数の扱いについて考えてみましょう。
まず、有理数クラスを自作してみるのはどうでしょうか。小さい値の範囲で連立一次方程式を解くときとかは役に立つかもしれません。欠点はすぐにオーバーフローしたり、毎回正規化でgcdを計算するのが重かったりすることです。
次に、全てを最初にn倍して計算して最後にnで割る、という方法があります。例えば幾何の問題なんかは、計算途中で整数の範囲から外れるけど入力で与えられる座標は全て整数、ということはよくあります。面積を求める問題で途中に三角形の重心を求める必要があるような問題があったとします。こんなときは最初に全ての座標を3倍しておいて、出てきた答えを9で割ると元に戻ったりします。浮動小数点数が出てくるのは最後に9で割る部分だけで、計算途中では整数しか出てきません。
そもそも整数の計算をしていて途中で非整数が出てくるのは割り算のせいです。なので割り算の式を見たら両辺の分母をはらって掛け算の形にできないかを考えると良いかもしれません。
無理数が出てきたら、これは浮動小数点数を使わざるを得ない場合が多いです。ただし、n乗根くらいしか出ないのであれば、両辺をn乗することで避けられることがあるかもしれません。
ただし、やりすぎはよくありません。浮動小数点数を避けることによって本質的に難しくなる問題というのも存在します(例えば、行列式の計算とか)。オーバーフローの危険も高まります。浮動小数点数を使うときはメリットとデメリットを比較して良い選択をしましょう。

おわりに

これまでに挙げた全てのトピックで自分はWAをもらっています。実際に体験するとやはり印象に残りやすいものです。
個人的には、浮動小数点数の適切な処理というのはメモリ制限や2倍以下の定数倍最適化と同程度に重要なことだと認識しています。
しかしそれは問題を解く側としての立場であって、問題を作る側に立つとコーナーケースを作る面倒さや検証の大変さから目を背けたくなってしまいます。しかも頑張ってコーナーケース作って撃墜しても不人気になるだろうし。誤差はあまり気にしなくていいですよ的な注釈が付いている問題も結構多いです。
ところでCompetitive Programming Advent Calendar、最初に予想していたよりもずっと濃い記事が多くて幸せです。もう全体の4割が終わってしまったのかー、という気分です。普段は読めないような記事がたくさん読めるのは純粋に嬉しいので、またこのような機会があることを願っております。