Page List

Search on the blog

2010年9月22日水曜日

多倍長演算 --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の内部構造
指数部の例外処理予約値
Posted by kenken
Labels:

0 件のコメント:

コメントを投稿

登録: コメントの投稿 (Atom)