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の精度ってそもそもどのくらいだったっけ？
ということで、今日は、
- doubleの最大値、最小値
- doubleの精度
- 丸め誤差
について、実際に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 件のコメント:
コメントを投稿