Search on the blog

2010年9月23日木曜日

多倍長演算 --doubleの精度!?

今日は、ついに多倍長演算デビュー。
PKU4問解いてやった。C++で自力で実装しようと思ってたが結局Javaでやることにした。信頼性の高い既成品があれば、当然そちらを使うべき。(あー、いつの間にかビジネスチックな思考になってしまった。。)

Javaは2種類の多倍長演算に対応する型を装備している。
まずは、BigInteger。その名の通り大きな整数。これはC++でいうlong long intの最大値より大きな値を扱いたい場合に使用する。long long intの最大値は、(2^64) / 2 - 1でだいたい8.0 * 10^18。
話はちょっとそれるが、『2^10 ≒ 10^3』と覚えておくと累乗の計算はすぐできるので便利。

2つ目は、BigDecimal。これはdoubleより精度の高い小数を扱うときに用いる。
BigIntegerの概念は分かりやすいが、BigDecimalは・・・。はて??doubleの精度ってそもそもどのくらいだったっけ?
ということで、今日は、
  1. doubleの最大値、最小値
  2. doubleの精度
  3. 丸め誤差
について、実際にPC上でプログラムを走らせながら検証します。(ちょっと長くなりそうだな~)

その前に、doubleの内部構造についておさらいします。doubleは1ビットの符号部、11ビットの指数部、52ビットの仮数部から成り立っています。
(詳しくは、参考ページをご参照ください。)
簡単に言うと10進数の指数表記(2.999 * 10e+10みたいなやつ)と同様のことを2進数でやってるだけです。

では、準備もできたので、まずは、doubleの最大値、最小値を求めてみましょう。指数部が11ビットなので概算では大体2^1024くらいがdoubleで表すことのできる最大値です。だいたい10^300くらいですねー。最小値は、2^(-1024)なので、10^(-300)程度だと思います。
では実際に計算しましょう。



#include <float.h>

#define EXP_MAX 2046 // 2047 is booked for exception procedure
#define EXP_MIN 1 // 0 is booked for exception procedure
#define EXP_BIAS 1023 // The biased exponent range: [-1022, 1023]
#define SIG_BIT 52 // significand bits

int main() {
// max
double dblmax = pow(2.0, EXP_MAX - EXP_BIAS);
double significand = 1.0;
REP(i, SIG_BIT)
significand += pow(2.0, -(i+1));
dblmax *= significand;
printf("dblmax = %e, DBL_MAX = %e\n", dblmax, DBL_MAX);

// min
double dblmin = pow(2.0, EXP_MIN - EXP_BIAS);
printf("dblmin = %e, DBL_MIN = %e\n", dblmin, DBL_MIN);

return 0;
}


んー、実際は、最大値が1.797693e+308最小値が2.225074e-308でした。ほぼほぼ予想通りです。指数部の値は、0と1024が例外処理のため予約語となっているようです。(参考ページ参照のこと。)やられた。これに気付かず1時間ぐらい悩んでしまった。。。

次に、doubleの精度を考えてみましょう。指数部が数値のスケール、仮数部が数値の精度を表しているので仮数部に注目すればOKです。52ビットなので、およそ10^15程度(有効数字15桁)の精度は持っているのではないでしょうか。(経験的にはそんなに無い気が・・・)実際にPC上で確認してみましょう。



int main() {
double min_significand = pow(2.0, -SIG_BIT);
printf("%.16f\n", min_significand);

double x = 1e15 + 1.0;
double y = 1e16 + 1.0;
printf("%.1lf\n", x);
printf("%.1lf\n", y);
}


んーと、仮数部の最小値は2.0 * 1e-16程度でした。スケールに対して2.0 * e-16程度の大きさの値までしか表現できないという意味なので、有効数字は15桁と考えていいでしょう。
実際に上のソースで検証しましたが、1.0 e+16に1を加算しても無視されてしまいます。
あれ、そうか!これが世間で言われている『情報落ち』ってやつか!!

最後に「丸め誤差」について。言わずと知れた「10進数を2進数に変換するときにでる誤差」のことです。代表的な例に、0.1を10回足しても1にならないというのがあります。
今回は、0.1をdoubleで表したときにどのような表現になるのか、またその値は10進数ではどの程度なのか調べたいと思います。標準関数frexp()を使えば簡単にできるようですが、なかなか面白そうなテーマなので、今回は自力で実装します。
ポイントは、『仮数部は52ビットあるので全件探索ではNG。何らかの”工夫”が必要』ってことです。まーでもこの場合は工夫と言ってもアレしかないですね(笑)



double getSigVal(LL sigbit) {
double ret = 1.0;

REP(i, SIG_BIT) {
if ((sigbit >> SIG_BIT - 1 - i) & 1)
ret += pow(2.0, -(i+1));
}
return ret;
}

double getDoubleVal(int exp, LL sigbit) {
return getSigVal(sigbit) * pow(2.0, exp - EXP_BIAS);
}

int main() {
double x = 0.1;
LL sig, sigl, sigr;
int expo = 1;

while (pow(2.0, expo + 1 - EXP_BIAS) <= x)
expo++;

sigl = 0LL;
sigr = (1LL << 52) - 1LL;
sig = (sigl+ sigr) / 2;

REP(i, 200) {
double val = getDoubleVal(expo, sig);
if (ABS(val - x) < 1e-18)
break;
else if (val > x)
sigr = sig;
else
sigl = sig;
sig = (sigl + sigr)/2;
}


printf("significand = %lld, exponent = %d\n", sig, expo);
printf("0.1 is expressed something like %.18lf on the PC\n", getDoubleVal(expo, sig));
return 0;
}

52ビット(10進数で10^15程度)を線形探索していたのでは、時間がありません。Binary Searchで解きましょう。精度は1e-18程度でいいでしょう。ループは200回くらい回せば探索範囲が1e-60倍くらいに縮小するので200回くらい回せばこの場合は十分です。(Binary Searchはバグを生みやすいので、無限ループはやめて必ず有限回の探索で打ち切るようにしましょう。
参考ページ:
doubleの内部構造
指数部の例外処理予約値

0 件のコメント:

コメントを投稿