Hatena::ブログ(Diary)

風と宇宙とプログラム このページをアンテナに追加 RSSフィード

2009-11-08

JavaScriptの大きな数と小さな数の仕組みを理解する 〜 IEEE754入門 〜

JavaScriptでの数値はIEEE754で規定されている倍精度型doubleです。符号部が1ビット、仮数部が52ビット、指数部が11ビットの64ビットで表現される浮動小数です。この辺りは、計算機の初歩の初歩で、当たり前すぎて普段は気にすることはないと思いますが、その境界値や特殊系について調べるといろいろ面白いことがわかります。ここでは、JavaScriptを例にしていますが、内容は一般的なもので、IEEE754の浮動小数入門的な話です。

では問題です。

整数として正確に表現できる最大の値はいくつか?

正確に表現できるというのは、n + 1 が確かに n + 1になることとします。n が非常に大きいときには、n + 1 は桁落ちが発生するので n のままです。考える前に実際にやってみましょう。探す n は 0 から 1e+100 の間にあるのは明らかなので、2分法で探索してみます。

function search_maxint(lo, hi) {
  while (hi - lo > 1) {
    var m = lo + (hi - lo) / 2;
    if (m == lo || m == hi) break;
    if (m + 1 == m) {
      hi = m;
    } else {
      lo = m;
    }
  }
  return lo;
}
var v = serach_maxint(0, 1e+100);
alert("v=" + v);

0から1e+100という大きな範囲での探索ですが、2分法なので一瞬で求まります。vの値は、

v=10295115178936058

になりました。果たして目的の数値でしょうか?

var v1 = 10295115178936058;
var v2 = v1 + 1;
alert("v2=" + v2);

不思議なことに、

v2=10295115178936060

となりました。1を足したのに2増えています。また、v2から1を引いてもv2のままです。理由は後で説明します。

ということは、上のプログラムには誤りがあり、n + 1した結果から1を引いたら元に戻るという条件が必要なようです。

function search_maxint(lo, hi) {
  while (hi - lo > 1) {
    var m = lo + (hi - lo) / 2;
    if (m == lo || m == hi) break;
    var m1 = m + 1;
    if ((m != m1) && (m1 - 1 == m)) {
      lo = m;
    } else {
      hi = m;
    }
  }
  return lo;
}
var v = serach_maxint(0, 1e+100);
alert("v=" + v);

今度は

v=9007199254740991

となりました。

実際に確認してみしょう。

var v1 = 9007199254740991;
var v2 = v1 + 1;             // v2は 9007199254740992
var v3 = v2 + 1;             // v3も 9007199254740992
var v4 = v2 - 1;             // v4は 9007199254740991

確かに求める数のようです。この数値をみて、16進数をすぐにイメージできる人は少ないと思いますが、16進数にすると

0x1FFFFFFFFFFFFF (1F_FFFF_FFFF_FFFF)

です。これに1を足した 0x20000000000000 が整数として正確に表現できる最大の数ということになります。

では、次はどうでしょう?

0より大きい最小の数値はいくつか?

これを求めるには、0と求める数nの間に0とnとは異なる数が存在しない、ということですから、0より十分に大きな数、たとえば1、を半分、半分にどんどんしていき、0かそれ自身になるまで続ければよいことになります。(本当にそれでよいのかちょっと自信はありませんが)

function search_minvalue(hi) {
  while (hi > 0) {
    var m = hi / 2;
    if (m == 0 || m == hi) {
      break;
    } else {
      hi = m;
    }
  }
  return hi;
}
var v = search_minvalue(1);
alert("v=" + v);

実行すると、

v=5e-324

になりました。

ちょっと確認してみましょう。

alert("v=" + v);                // v=5e-324
alert("v=" + 1/v);              // v=Infinity
alert("v=" + 3e-324);           // v=5e-324
alert("v=" + 2e-324);           // v=0
alert("v=" + 7e-324);           // v=5e-324

3e-324も7e-324もどちらも5e-234になります。ミクロな数値は飛び飛びの値しか取れません。なんか量子力学みたいです。

この5e-324という数値はJavaScriptでは Number.MIN_VALUE 定数として定義されています。

では、さらに、次はどうでしょう?

Infinityではない最大の数値はいくつか?

これまでは、2分法で範囲の半分の値を利用していましたが、Inifinityがあるとその半分もInfnityなのでこれまでの方法は使えません。Infinityでない数値から、じわじわとInfinityに近づいていく方法をとります。qを1に近い1より小さい数として、スタートとなる数を(1+q^0)倍、さらに(1+q^0)倍していき、Infinityになった1つ前の数値にもどり、そこからさらに今度は(1+q^1)倍、(1+q^1)倍していき同様にInifinityになるまで繰り返し、次は(1+q^2)倍、(1+q^2)倍していくという方法をとります。つまり、倍数を1+q^n ずつ下げていきます。

ここでは、qを0.9とします。

function search_maxvalue(lo) {
  var r = 1;
  var q = 0.9;
  do {
    var p = 1 + r;
    var m = lo;
    do {
      var hi = m * p;
      if (isFinite(hi)) {
        m = hi;
      } else {
        lo = m;
        break;
      }
    } while (p > 1);
    r *= q;
  } while (p > 1);
  return lo;
}
var v = search_maxvalue(1);
alert("v=" + v);

実行すると

v=1.7976931348623157e+308

になりました。思ったより速く収束します。

本当にこの数値がInfinityではない最も大きな数でしょうか。確認してみます。

alert("v=" + 1.7976931348623157e+308);          // v=1.7976931348623157e+308
alert("v=" + 1.79769313486231571e+308);         // v=1.7976931348623157e+308
alert("v=" + 1.7976931348623158e+308);          // v=1.7976931348623157e+308
alert("v=" + 1.7976931348623156e+308);          // v=1.7976931348623155e+308
alert("v=" + 1.7976931348623159e+308);          // v=Infinity

確かに、なんとなく正しそうです。

実は、この値はJavaScriptでは Number.MAX_VALUE で定義されている定数と一致します。上の計算では実際には桁落ちや丸め誤差が発生するので、qの値や初期値によっては若干異なる結果となる場合があります。

IEEE754の数値表現

長々と当たり前のことを書いてしまいましたが、これらのことはJavaScriptが採用している数値表現であるIEEE754から当然の帰結となります。

IEEE754の数値のビット配列は以下のようになります。

 1     11                           52
+-+-----------+----------------------------------------------------+
|S|   Exp     |                  Fraction                          |
+-+-----------+----------------------------------------------------+
63 62       52 51                                                 0

Sは符号で、負のとき1ビットがセットされます。Expは指数部、Fractionは仮数部です。ここで注意が必要なのは、指数部には負の数が表現できますが、符号ビットを点けてしまうと64ビット全体を単純な数値とみなしたとき、大小比較が面倒になるため、指数部は符号ビットはつけないで、1023をプラスしてシフトさせた値がセットされます。

例として、123.456 はどのように表現されるかみてみましょう。

そのためには2進数で表現する必要がありますが、123は2進数では 1111011 となることは誰でも計算できると思います。では、0.456 を2進数で表現するとどうなるでしょうか。分らないという人も意外に少なくありません。

整数部分を2進数にするには、2で次々に割っていき、余りを逆に並べればよいので、手計算でやると下のようになります。

    123
    ---
     61  … 1
     30  … 1
     15  … 0
      7  … 1
      3  … 1
      1  … 1

123(2) = 1111011

一方、小数部分は次々に2倍していき、整数部分を取り出せばよいことになります。

      0.456
      -----
    0  .912
    1  .824
    1  .648
    1  .296
    0  .592
    1  .184
    0  .368
    0  .736
    1  .472
    0  .944
    1  .888
    1  .776
    1  .552
    1  .104
    0  .208
    0  .416
    ......
0.456(2) = 0.0111010010111100....

2進数は任意の数値を

b_n¥times2^n + b_{n-1}¥times2^{n-1} + ... + b_0¥times2^0 + b_{-1}¥times2^{-1} + b_{-2}¥times2^{-2} +...

で表現することだから当たり前といえば当たり前です。

もちろんこんな手計算をしなくても、JavaScriptを使っているなら、以下のように簡単に求まります。

var v = (123.456).toString(2);
alert("v=" + v);
// v=1111011.0111010010111100011010100111111011111001110111

これを正規化した指数表現で表すと以下のようになります。

1.1110110111010010111100011010100111111011111001110111×2^6

ゼロ以外は整数部分は必ず1になるわけですから、52ビットの仮数部にはこの1は省略して小数部分の52ビットがセットされます。また、指数部は1023プラスされているので1029を2進数にした 10000000101 がセットされます。

つまり、123.456 をIEEE754で表現すると下のようになります。

 S     Exp                    Fraction 
+-+-----------+----------------------------------------------------+
|0|10000000101|1110110111010010111100011010100111111011111001110111|
+-+-----------+----------------------------------------------------+
63 62       52 51                                                 0

本当でしょうか?JavaScript自身で確認したいところですが、JavaScriptはビット演算が32ビットに限定されるため使えません。C言語で確認してみましょう。C言語は仕様的にはIEEE754に準拠しているわけではありませんが、事実上対応しています。

#include <stdio.h>
char* getbitstr(double val)
{
  static char buf[67];
  uint64_t x = *(uint64_t*)&val;
  uint64_t flg = 1LL << 63; 
  int j = 0;
  for (int i = 0; i < 64; i++) {
    buf[j++] = (x & flg) ? '1' : '0';
    if (i == 0 || i == 11) {
      buf[j++] = '_';
    }
    flg >>= 1;
  }
  buf[j] = 0;
  return buf;
}
int main() {
  printf("123.456=%s\n", getbitstr(123.456));
}

実行すると、

123.456=0_10000000101_1110110111010010111100011010100111111011111001110111

と出力されました。ぴったり合っています。


では次に、JavaScriptの最大の数であるNumber.MAX_VALUEと、最小の数であるNumber.MIN_VALUEがどのように表現されているかみてみましょう。

上のCの関数getbitstr()で実際に実行してみます。結果を書くと下のようになります。

MIN_VALUE=0_00000000000_0000000000000000000000000000000000000000000000000001
MAX_VALUE=0_11111111110_1111111111111111111111111111111111111111111111111111

MIN_VALUEの方は仮数の最後のビットが1で、それ以外は全部ゼロになります。確かに最小の数のようです。でもちょっと待ってください。仮数部は整数部分の1が省略されたものですから、上のMIN_VALUEの値は

1.0000000000000000000000000000000000000000000000000001 × 2^-1023

という値になる筈です。しかし、これは 1.1125e-308 付近の値であり、MIN_VALUE の5e-324とは異なります。どういうことでしょう?実は、指数部が全部ゼロのときは、非正規数と呼ばれる数の一種であり、整数部分の1が省略されたいう意味は持ちません。なので、

0.0000000000000000000000000000000000000000000000000001 × 2^-1023
= 1.0 × 2^(-1023-51)
= 5e-324

となるわけです。

一方、MAX_VALUEの方は仮数は全部1ビットが立っていますが、指数部の最後のビットはゼロです。実は、IEEE754では指数部のビットが全部1のときは特殊な数を表現します。指数部が全部1で仮数部が全部ゼロの数はInfinityを意味します。さらに、仮数部にゼロでないものがある場合にはNaNを意味します。getbitstr()関数確認すると確かにそのようになっています。

int main() {
  double inf = 1.0/0.0;
  double nan = 0.0/0.0;
  printf("Inf=%s\n", getbitstr(inf));
  printf("NaN=%s\n", getbitstr(nan));
  printf("NaN=%s\n", getbitstr(-nan));
}
// 実行結果
// Inf=0_11111111111_0000000000000000000000000000000000000000000000000000
// NaN=1_11111111111_1000000000000000000000000000000000000000000000000000
// NaN=0_11111111111_1000000000000000000000000000000000000000000000000000

NaNにも符号がありますが、意味は持ちません。また仮数部はゼロでなければよいので、任意の値をセットすることができ、プロセッサの診断などに使われることがあるようです。

ところで、Number.MAX_VALUEの値は正確にはいくつになるでしょうか。

MAX_VALUE=0_11111111110_1111111111111111111111111111111111111111111111111111
 = 1.1111111111111111111111111111111111111111111111111111 × 2^1023
 = 11111111111111111111111111111111111111111111111111111 × 2^971
 = 17976931348623157081452742373170435679807056752584499659891747680315
   72607800285387605895586327668781715404589535143824642343213268894641
   82768467546703537516986049910576551282076245490090389328944075868508
   45513394230458323690322294816580855933212334827479782620414472316873
   8177180919299881250404026184124858368

確かに、1.7976931348623157e+308 と一致します。

最後に、正確に表現できる最大の整数と、1を足したのに何故2増えたのかを検証してみます。

最大の整数は小数部分がないということですから、小数部分を全部繰り上げるように指数部が設定されていればよいので、下のようになります。

1.1111111111111111111111111111111111111111111111111111 × 2^52
 = 9007199254740991

前に計算で求めた値と一致しました。

次に、1を足したのに2つ増えたというのは、小数計算の丸め誤差と類似のものです。

問題の数 10295115178936058 をIEEE754で表現すると

0_10000110100_0010010010011010110100100101100101001100001101111101

であり、指数部は 53 (1076-1023)です。つまり、仮数部の有効桁より1つ大きいため、10295115178936058は正確には表現できていません。これより1つ大きい10295115178936059は桁落ちしてしまい丸められます。なので、この数以外にも指数部が 53 で仮数部の最下位ビットが1の数は1を足すと2増えることになります。

まとめ

たかが浮動小数点数の数値表現にもこのように奥の深い話があります。最近の新入社員の話を聞くと、初めて覚えた言語はRubyで、C言語は知らないという人も中にはいるようです。C言語を使うときには、intやdoubleなど明確に型を区別し、桁落ちや丸め誤差などは意識すると思いますが、JavaScriptのような型なし言語で、しかも、整数から浮動小数に必要に応じて自動的に変わってしまう言語では思わぬところに落とし穴があります。JavaScriptでこのようなことを気にする必要があるプログラムを書く機会はあまりないと思いますが、IEEE754はJavaScriptが言語仕様として採用している規格であり、知っておいて損はないでしょう。

実際のJavaScriptの実装では、全ての数をIEEE754のdouble型浮動小数で計算しているわけではありません。別の機会にブラウザ毎にこのあたりを確認してみた結果をレポートできればと思います。

bushimichibushimichi 2009/11/09 10:25 実際に利用するシーンがあまりなかったせいもあって
あいまいな認識で適当に流していた箇所だったので
非常に勉強になりました。ありがとうございます。