はじめに
最近(2018年11月)、Togetterまとめ30年前のBASICで「0.01を10000回足したら100.003」になると書いてあったので、今の環境で試してみたら同じ結果だった話がちょっと話題になりました。
例えば、似たような話として、Ruby(x86_64-Linux)で0.01を100回足しても、ぴったり 1 にはならず、少し誤差が出ます ( 次の実行例参照 )。他のプログラミング言語でも大体似たような状況になるでしょう。
$ irb
irb(main):001:0> 100.times.reduce(0){|s,_| s+0.01 }
=> 1.0000000000000007
件のまとめでは「今の時代になって…」という意見も出てたりしますが、浮動小数点が何か知ってる人にはアタリマエの知識であっても、そうでないと単なる小数データとしか見えなくて、下手すると「計算でズレが出るなんてケシカラン!」という印象を抱くことも少なくないようです。
…本来であれば、(大抵の言語に標準で備わっているデータの種類なので)入門書等で知識を身に着けておくものかと思いますが、ちゃんと説明があるか?というとあまり無いようにも思いますので、ざっくり記事にすることにしました。
浮動小数点って何?
固定小数点と浮動小数点
一言で言えば、浮動小数点というのは一般の人が想像するような小数ではなく、計算で誤差が発生することを前提とした数値データのことです。
逆に言えば、浮動小数点を使う時は常に近似計算であることを意識し、結果の精度を考えろということでもあります。
それだけだと、何でそんなイケてないことするの? と疑問に思われるかも知れませんが、それは有限のリソースで計算を行うコンピュータの宿命であって、どこを妥協するかのポイントの違いです。
浮動小数点の逆は固定小数点であり、一般に言う整数データが該当します。これは、必ず1刻みで扱え、整数を扱う限り誤差は出ませんが扱える数値の範囲がそれほど広くありません。
※例えば64bit整数では -9,223,372,036,854,775,808～9,223,372,036,854,775,807 の範囲で、10進数として19桁です。
一方で、浮動小数点は誤差が出る代わりに非常に広い範囲の数を扱うことができます。0にほとんど近いナノレベルの数値から、(最近定義が改定された)アボガドロ定数や宇宙規模の数値でも対応可能です。
※例えば、倍精度と呼ばれる64bit浮動小数点だと、0に最も近いところで $10^{-308}$ 程度、大きいところで$10^{308}$程度、実に10進数で600桁程度の範囲をカバーします。ただし精度は10進で16桁程度であるため、この桁数に収まらない分の計算結果は失われ、誤差になっていきます。
そして、大抵のプログラミング言語では、小数データを浮動小数点で扱うことが一般的なのです。
なぜ小数を浮動小数点で?
そうは言っても、こう思われる方もいるかも知れません。
「いやいや、そんなナノとか宇宙とか要らないから。せいぜい0.01程度扱うのに、いちいち誤差なんて勘弁してくれ」と。
実際問題として、通貨や税率等で出てくる小数しか扱わないならその感覚も尤もです。
しかし、じゃあ、コンピュータとして小数第何位までカバーすれば十分でしょうか?
と、ここで、既に小学校で習っているはずです。例えば円周率は無限小数であることを ( 3.14はあくまで近似値 )。
円周率を持ち出さなくても、ピタゴラスさんが存在を封印しようとした無理数$\sqrt{2}$も無限小数であって、じゃあどこまで桁数を見れば良いのか、これは小数の計算で常に付きまとってきます。
結局、有限で表せない以上どこかで打ち切らざるを得ず、計算精度と誤差を考えようということになるわけです。いわゆる科学技術計算もそうですが、統計、グラフィック処理等々、この問題は避けては通れないのです。
※小学校範囲での「概数」、科学分野での「有効数字」に即した考え方でもあります。
…ということもあり、これは私の推測に過ぎませんが、整数は固定小数点、小数は浮動小数点という棲み分けになったのでしょう。
それでも誤差を嫌うなら
そうは言っても、誤差のことなんて考えたくないという場面もあるでしょう。
それであれば、言語なりライブラリで、誤差の発生しない小数点の取り扱いができるものがある筈です。
例えばレガシーで悪名高いCOBOLでは、packed decimalという形式で、指定した桁数の範囲の小数では誤差のない計算ができるようになっています。流石事務処理用と言われただけはあります。
他の言語では、例えばJavaならBigDecimal、C#ならdecimal、RubyならBigDecimal ( 整数÷整数で済む話ならRationalも ) を使うという選択肢があります。
或いは、例えば消費税率 8% を掛ける位の話なら ( その程度なら浮動小数点でもそうそう誤差は出ませんが )、単純に整数演算に置き換えれば済む話です。(金額)*0.08 ではなくて (金額)*8/100のように。
或いは、1ドル50セントを1.50ドルと表現する代わりに150セントと単位を工夫して、整数の範囲に収まるようにするという方法も考えられるでしょう。
要は、現実世界で小数だからと言ってプログラム上で安易に小数として扱うな、まずは考えろということになるでしょう。
誤差の生まれる仕組み
2進数表現
しかし一方で、「なんで0.01という何でもない数で誤差が出るんだろう?」という疑問があるかも知れません。
次のような、とても大きな数と小さな数を足すことで、表現できる精度に収まらなかったというのは、まだ分かり易いところでしょう。( Ruby(x86_64-Linux)では倍精度浮動小数点として扱うので、精度は10進で16桁程度 )
$ irb
irb(main):001:0> 100_000_000 + 0.000_000_01 # 精度に収まる
=> 100000000.00000001
irb(main):002:0> 1000_000_000 + 0.000_000_001 # 収まらない
=> 1000000000.0
しかし、0.01を100回足して1というのは、別に大きさに開きがあるわけではないし…? と。
この「0.01を100回足しても1にならない」というケースでは、コンピュータが内部的に浮動小数点を2進数として扱っていて、0.01が2進数では無限(循環)小数となることが原因となっています。
※コメント頂いたように10進数ベースの浮動小数点もあり、そちらを使う場合は該当しません。
10進数の小数と、2進数の小数とは、次のような対応があることに注意します。
$$
\begin{eqnarray}
0.5&(10)=&0.1&(2) \\
0.25&(10)=&0.01&(2) \\
~ &~&\vdots&~ \\
0.015625&(10)=&0.000001&(2) \\
0.0078125&(10)=&0.0000001&(2) \\
0.00390625&(10)=&0.00000001&(2) \\
0.001953125&(10)=&0.000000001&(2) \\
~ &~&\vdots&~ \\
\end{eqnarray}
$$
そうすると、0.01は次のような20桁周期での2進循環小数であることが(計算してみると)分かります。
$$
0.01(10)=0.000000101000111101011100001010001111010111\cdots(2)
$$
なので、コンピュータでは桁を途中で打ち切って保存しておくことになります。この打ち切った桁が、計算を重ねるうちに誤差として影響してくるのです。
逆に言えば2進数的にキリのいい(有限小数になる)、分母が2のべき乗になる分数(例えば $\frac 1 {128}=0.0078125(10)=0.0000001(2)$)では誤差が生まれません。
誤差があること前提で取り扱う心づもりは必要ですが、必ず誤差が出るとは限らないのです。
なお、打ち切られると言っても、「丸め」によって端数処理が行われますので、結果が減る(切捨てられる)とは限りません。切り上げであったり、あるいは近い方へ寄せる(四捨五入のような処理)等、ポリシーは処理系によります。
誤差と2進表現のデモ
さて。浮動小数点は、今日ではIEEE754という規格で定められています。
※コメントで頂いた10進の浮動小数点も、2008年版のIEEE754-2008で取り入れられていますが、以下2進ベースの話で進めます。
最もメジャーなのが64bitの倍精度ですが、それ以外にも32bitの単精度、128bitの4倍精度、また昨今流行りの機械学習では、単精度よりも精度を落とした(代わりにハードによっては計算量を稼ぐことができる)16bitの半精度も使われます。
内部の表現については浮動小数点数を理解するが詳しいですが、いずれにせよ、次のような2進数表現を元にしています。
$$
(-1)^s\cdot 1.fff\cdots ff(2)\cdot 2^{e-bias}
$$
ということで、Cのプログラムで浮動小数点の内部表現を元に2進数としてどうなっているかを見て、誤差が生まれていく様子を筆算で見られるものを作成しました。
今回は、(桁数が多いと見辛いので)単精度で。試したWindows10/WSL+Ubuntu16+gcc(x86_64)の環境では、float型が相当します。
※intelの環境だと大体そうなると思いますが
コンパイルと実行の様子は次のようになります。100 を入力として与えることで、$\frac 1 {100}=0.01$の2進数表示を出力した後、筆算の様子を出力します。( 右側には、同等の掛け算を計算した場合の2進数表示も付けます。 )
最後に、足し算の結果も10進小数で出力します。
$ gcc -std=c99 -o float float.c
$ ./float <<< 100
in single precision calculation:
1.0/100 = 0.01(10) = 1.01000111101011100001010(2) * 2^-7
continuous binary column additions
0.000000101000111101011100001010 0.000000101000111101011100001010 (1/100=0.01)
+ 0.000000101000111101011100001010
---------------------------------
0.00000101000111101011100001010 0.00000101000111101011100001010 (2/100=0.02)
+ 0.000000101000111101011100001010
---------------------------------
0.00000111101011100001010001111 0.00000111101011100001010001111 (3/100=0.03)
+ 0.000000101000111101011100001010
---------------------------------
0.0000101000111101011100001010 0.0000101000111101011100001010 (4/100=0.04)
+ 0.000000101000111101011100001010
---------------------------------
0.0000110011001100110011001100 0.0000110011001100110011001101 (5/100=0.05)
…(中略)
0.111110101110000100111101 0.111110101110000101001000 (98/100=0.98)
+ 0.000000101000111101011100001010
---------------------------------
0.111111010111000010011001 0.111111010111000010100100 (99/100=0.99)
+ 0.000000101000111101011100001010
---------------------------------
0.111111111111111111110101 1.00000000000000000000000 (100/100=1)
1.0/100+1.0/100+...+1.0/100(100times) = 0.99999934
単精度の場合、精度は2進数24桁分となります。そのため、足し算を重ねるごとに、保持している合計と足す数との桁のズレが出てきて、打ち切られた分の誤差が蓄積していくことになります。 ( 既に5回目で最後の桁に誤差が出てきているのが分かります )
結果的に、単精度で0.01を100回足した場合は0.99999934と、1を僅かに下回る数値になります。
デモのソース
以下、デモのソース(C99)です。ideoneあたりに貼り付けて実行するのも良いかと思います。
#include <stdio.h>
#include <stdint.h>
#include <string.h>
#include <assert.h>
void test(int n);
int main(void) {
int n=0;
int ret=scanf("%d",&n);
if ( ret!=1||n<=0 )
n=100;
test(n);
}
void binform_s1(char *dst,double d);
int binform_s2(char *dst,double d);
void test(int n) {
const float t=1.0/n;
char buf[256];
binform_s1(buf,t);
printf("in single precision calculation:\n"
" 1.0/%d = %g(10) = %s\n\n",n,t,buf);
float s=0;
int c1=binform_s2(buf,t);
char line[256];
memset(line,'-',c1);
line[c1]='\0';
printf("continuous binary column additions\n\n");
for ( int i=1; i<=n; i++ ) {
s+=t;
if ( i>1 )
printf("+%s\n"
" %s\n",buf,line);
char buf2[256],buf3[256];
int c2=binform_s2(buf2,s);
(void)binform_s2(buf3,(double)i/n);
printf(" %s%*s %s (%d/%d=%g)\n",buf2,c1-c2,"",buf3,i,n,(double)i/n);
}
printf("\n"
" 1.0/%d+1.0/%d+...+1.0/%d(%dtimes) = %.8f\n",n,n,n,n,s);
}
// 以下、補助関数
// type punningによる、整数型(bit列)への変換
static inline int32_t f2i(double d) {
// type punning
union {
float f;
int32_t i;
} u = { .f=d };
return u.i;
}
// 単精度から2進表現文字列 ( 仮数*2^指数 ) を生成
void binform_s1(char *dst,double d) {
int32_t x=f2i(d);
sprintf(dst,"%c1.%23s(2) * 2^%d",
x>>31?'-':' ',"",((x>>23)&255)-127);
for ( int i=22; i>=0; i-- )
dst[25-i]='0'+(x>>i&1);
}
// 単精度から2進表現文字列 ( 1.xx… あるいは 0.xx… を生成 )
int binform_s2(char *dst,double d) {
int32_t x=f2i(d);
int r=((x>>23)&255)-127;
// support only small floats
assert(r<=0);
dst[0]=x>>31?'-':' ';
memset(dst+1,'0',1-r);
dst[2]='.';
dst[r<0?2-r:1]='1';
for ( int i=22; i>=0; i-- )
dst[25-r-i]='0'+((x>>i)&1);
dst[26-r]='\0';
return 26-r;
}
まとめ
最後にまとめです。
- 浮動小数点は誤差の発生を前提とした数値データである
- 代わりに、固定小数点に比べ、非常に広い範囲の数値を扱うことができる
- 大抵のプログラミング言語では、小数は浮動小数点として扱われる
- 浮動小数点による誤差を望まないならば、適切なライブラリ・データ型を使うなり、整数操作で扱えるように考えること
- 浮動小数点は内部的に2進小数で表現されている
- 人間的に(10進数として)キリがいい小数であっても、一部を除き2進数的には無限小数になり、桁が打ち切られるため誤差の原因となる