数値的に正しいJIS丸めの実装(1)
JIS丸めとは最近接偶数への丸めのことでありまして、データがランダムの場合、四捨五入と違いバイアスが生じないとされる（世間的にそう言われている）丸め方法であります。*1
普通に書くとつぎよのうなコードになると思われ：
/// JIS丸めを行います。 double mylib::round2(double x) { double fn = floor(x); // 整数部 double fr = x - fn; // 端数 if (fr < 0.5) { return fn; } else if (fr > 0.5) { return fn + 1.0; } else { double r = fmod(fn, 2.0); assert(r == 0.0 || r == 1.0); return (r > 0) ? fn + 1.0 : fn; } }
しかしこれでは人間（特に他人）の作ったfloor()とか、人間（特に他人）の作ったfmod()とか、あるいはx - fnあたりで何か精度の劣化が起きやしないかとかetc.という疑念が渦巻き渦巻き、数値的な問題がほんとうに無いのかかなりｶﾅﾘカナーリ気になる。
そこで今日わ、IEEE 754形式の倍精度浮動小数点数*2を対象として、JIS丸めを整数演算に展開して書いてみて、（漏れが）気分的にすっきり（して年越し）するとともに限界を明らかにしようとぞ思ふ、
IEEE 754形式倍精度浮動小数点数
ｳｨｷでｸﾞｸﾞれば出てくるが、IEEE 754形式倍精度浮動小数点数とは符号ビット1ビット、指数部11ビット、仮数部52ビットで数を表したものであり、それぞれbNeg, exp, fracと表すことにしたとき、NaN、±∞、非正規化数という例外ケースを除き、次の数を約表す：
*3 * 2^((指数部 - 1023) - 52)
という絶対値を約表すというしくみ。（ただし正規化数の場合。また符合の存在は無視した。）
正規化数の代表例を分解して確かめてみると、約こんなコード↓↓↓で、
union DoubleOrLL { double fx; __int64 llx; }; /// IEEE 754形式倍精度浮動小数点数を分解します。 void mylib::DecodeIEEE754Double(double x, bool& bNeg, int& exp, __int64& frac) { DoubleOrLL buf; buf.fx = x; __int64 llx = buf.llx; bNeg = ((llx & 0x8000000000000000LL) != 0); exp = ((llx >> 52) & 0x7ff); frac = (llx & 0x000fffffffffffffLL); } /// IEEE 754形式倍精度浮動小数点数を合成します。 double mylib::GenIEEE754Double(bool bNeg, int exp, __int64 frac) { DoubleOrLL buf; buf.llx = 0LL; if (bNeg) { buf.llx |= 0x8000000000000000LL; } buf.llx |= (((__int64)exp << 52) & (0x7ffLL << 52)); buf.llx |= (frac & 0x000fffffffffffffLL); return buf.fx; } double pow2(int e) { double x = 1.0; if (e > 0) { for (int i = 1; i <= e; ++i) { x *= 2.0; } } else { e = -e; for (int i = 1; i <= e; ++i) { x /= 2.0; } } return x; } void DumpIEEE754Double(const char* eqn, double x) { using namespace std; bool bNeg; int exp; __int64 frac; mylib::DecodeIEEE754Double(x, bNeg, exp, frac); cerr.setf(ios::scientific, ios::floatfield); cerr << boolalpha << setprecision(DBL_DIG + 1) << setfill(' ') << setw(22) << eqn << " = " << setw(26) << x << "; bNeg=" << setw(5) << bNeg << ", exp=" << setw(4) << exp << " (1023+(" << setw(5) << (exp - 1023) << ")), frac=0x" << hex << setfill('0') << setw(13) << frac << dec << std::endl; } #define DUMP_DBL(eqn) DumpIEEE754Double(#eqn, (eqn)) /// いろんな数を分解表示してみる。 void experiment01a() { using namespace std; using namespace mylib; #undef min #undef max DUMP_DBL(DBL_EPSILON); DUMP_DBL(DBL_MAX); DUMP_DBL(DBL_MIN); DUMP_DBL(1); DUMP_DBL(0); DUMP_DBL(-0); DUMP_DBL(-1); DUMP_DBL(numeric_limits<double>::epsilon()); // こう書いてもDBL_EPSILONと同じ DUMP_DBL(numeric_limits<double>::max()); // こう書いてもDBL_MAXと同じ DUMP_DBL(numeric_limits<double>::min()); // こう書いてもDBL_MINと同じ cerr << "#" << endl; DUMP_DBL((1LL << 52) + 0); DUMP_DBL((1LL << 52) + 1); // 上の値と識別される DUMP_DBL((1LL << 53) + 0); DUMP_DBL((1LL << 53) + 1); // 上の値と識別されない cerr << "#" << endl; DUMP_DBL(-((1LL << 52) + 0)); DUMP_DBL(-((1LL << 52) + 1)); // 上の値と識別される DUMP_DBL(-((1LL << 53) + 0)); DUMP_DBL(-((1LL << 53) + 1)); // 上の値と識別されない cerr << "#" << endl; DUMP_DBL(((1LL << 53) - 1) * pow2(1023 - 52)); // DBL_MAXと一致 DUMP_DBL(((1LL << 53) - 0) * pow2(1023 - 52)); // +INF cerr << "#" << endl; DUMP_DBL(-((1LL << 53) - 1) * pow2(1023 - 52)); // -DBL_MAXと一致 DUMP_DBL(-((1LL << 53) - 0) * pow2(1023 - 52)); // -INF cerr << "#" << endl; __int64 all1 = (1LL << 52) - 1; DUMP_DBL(GenIEEE754Double(false, 2046, all1 - 0)); // DBL_MAXと一致 DUMP_DBL(GenIEEE754Double(false, 2046, all1 - 1)); // DBL_MAXの前の数 DUMP_DBL(GenIEEE754Double(false, 1, 1)); // DBL_MINの次の数 DUMP_DBL(GenIEEE754Double(false, 1, 0)); // DBL_MINと一致 cerr << "#" << endl; DUMP_DBL(GenIEEE754Double(true, 1, 0)); // -DBL_MINと一致 DUMP_DBL(GenIEEE754Double(true, 1, 1)); // -DBL_MINの前の数 DUMP_DBL(GenIEEE754Double(true, 2046, all1 - 1)); // -DBL_MAXの次の数 DUMP_DBL(GenIEEE754Double(true, 2046, all1 - 0)); // -DBL_MAXと一致 cerr << "#" << endl; }
約こんな実行結果を得る↓↓↓
DBL_EPSILON = 2.2204460492503131e-016; bNeg=false, exp= 971 (1023+( -52)), frac=0x0000000000000
DBL_MAX = 1.7976931348623157e+308; bNeg=false, exp=2046 (1023+( 1023)), frac=0xfffffffffffff
DBL_MIN = 2.2250738585072014e-308; bNeg=false, exp= 1 (1023+(-1022)), frac=0x0000000000000
1 = 1.0000000000000000e+000; bNeg=false, exp=1023 (1023+( 0)), frac=0x0000000000000
0 = 0.0000000000000000e+000; bNeg=false, exp= 0 (1023+(-1023)), frac=0x0000000000000
-0 = 0.0000000000000000e+000; bNeg=false, exp= 0 (1023+(-1023)), frac=0x0000000000000
-1 = -1.0000000000000000e+000; bNeg= true, exp=1023 (1023+( 0)), frac=0x0000000000000
numeric_limits<double>::epsilon() = 2.2204460492503131e-016; bNeg=false, exp= 971 (1023+( -52)), frac=0x0000000000000
numeric_limits<double>::max() = 1.7976931348623157e+308; bNeg=false, exp=2046 (1023+( 1023)), frac=0xfffffffffffff
numeric_limits<double>::min() = 2.2250738585072014e-308; bNeg=false, exp= 1 (1023+(-1022)), frac=0x0000000000000
#
(1LL << 52) + 0 = 4.5035996273704960e+015; bNeg=false, exp=1075 (1023+( 52)), frac=0x0000000000000
(1LL << 52) + 1 = 4.5035996273704970e+015; bNeg=false, exp=1075 (1023+( 52)), frac=0x0000000000001
(1LL << 53) + 0 = 9.0071992547409920e+015; bNeg=false, exp=1076 (1023+( 53)), frac=0x0000000000000
(1LL << 53) + 1 = 9.0071992547409920e+015; bNeg=false, exp=1076 (1023+( 53)), frac=0x0000000000000
#
-((1LL << 52) + 0) = -4.5035996273704960e+015; bNeg= true, exp=1075 (1023+( 52)), frac=0x0000000000000
-((1LL << 52) + 1) = -4.5035996273704970e+015; bNeg= true, exp=1075 (1023+( 52)), frac=0x0000000000001
-((1LL << 53) + 0) = -9.0071992547409920e+015; bNeg= true, exp=1076 (1023+( 53)), frac=0x0000000000000
-((1LL << 53) + 1) = -9.0071992547409920e+015; bNeg= true, exp=1076 (1023+( 53)), frac=0x0000000000000
#
((1LL << 53) - 1) * pow2(1023 - 52) = 1.7976931348623157e+308; bNeg=false, exp=2046 (1023+( 1023)), frac=0xfffffffffffff
((1LL << 53) - 0) * pow2(1023 - 52) = 1.#INF000000000000e+000; bNeg=false, exp=2047 (1023+( 1024)), frac=0x0000000000000
#
-((1LL << 53) - 1) * pow2(1023 - 52) = -1.7976931348623157e+308; bNeg= true, exp=2046 (1023+( 1023)), frac=0xfffffffffffff
-((1LL << 53) - 0) * pow2(1023 - 52) = -1.#INF000000000000e+000; bNeg= true, exp=2047 (1023+( 1024)), frac=0x0000000000000
#
GenIEEE754Double(false, 2046, all1 - 0) = 1.7976931348623157e+308; bNeg=false, exp=2046 (1023+( 1023)), frac=0xfffffffffffff
GenIEEE754Double(false, 2046, all1 - 1) = 1.7976931348623155e+308; bNeg=false, exp=2046 (1023+( 1023)), frac=0xffffffffffffe
GenIEEE754Double(false, 1, 1) = 2.2250738585072019e-308; bNeg=false, exp= 1 (1023+(-1022)), frac=0x0000000000001
GenIEEE754Double(false, 1, 0) = 2.2250738585072014e-308; bNeg=false, exp= 1 (1023+(-1022)), frac=0x0000000000000
#
GenIEEE754Double(true, 1, 0) = -2.2250738585072014e-308; bNeg= true, exp= 1 (1023+(-1022)), frac=0x0000000000000
GenIEEE754Double(true, 1, 1) = -2.2250738585072019e-308; bNeg= true, exp= 1 (1023+(-1022)), frac=0x0000000000001
GenIEEE754Double(true, 2046, all1 - 1) = -1.7976931348623155e+308; bNeg= true, exp=2046 (1023+( 1023)), frac=0xffffffffffffe
GenIEEE754Double(true, 2046, all1 - 0) = -1.7976931348623157e+308; bNeg= true, exp=2046 (1023+( 1023)), frac=0xfffffffffffff
#(double)(1LL << 52) と (double)((1LL << 52) + 1) が識別されているのに対し、
(double)(1LL << 53) と (double)((1LL << 53) + 1) が識別されていないのに注意されたい
特殊な数
NaN (Not a Number)
std::sqrt(-1.0)等、実数で無い数がNaNになるっぽい
bNegはtrueまたはfalse、指数部は2047、仮数部は非0。*4
NaNは、NaN自身も含め、どの浮動小数点数とのいかなる比較演算も偽となる
±無限大
1/0とか、1/-0とか、無限大になるとき現れる
bNegはtrueまたはfalse、指数部は2047、仮数部は0。
非正規化数
DBL_MINの次の浮動小数点数 - DBL_MIN等、0に近く、かつ互いにきわめて近接した数の差等で生じる。
実際、
DBL_MINの次の浮動小数点数 - DBL_MIN
= (1.0000000000 0000000000 0000000000 0000000000 0000000000 01)_bin * 2^(-1022)
- (1.0000000000 0000000000 0000000000 0000000000 0000000000 00)_bin * 2^(-1022)
= (0.0000000000 0000000000 0000000000 0000000000 0000000000 01)_bin * 2^(-1022)
なんだけど、最後の数値を「仮数部の53 bit目を1とみなす」正規化数の約束事に従って表そうとすると、
(1.0000000000 0000000000 0000000000 0000000000 0000000000 00)_bin * 2^(-1022 - 52)
となって、指数部の値として 1023 + (-1022 - 52) = -51 を要し、アンダーフローしてしまう*5
興味深いのは、この現象はただ端に近接した数同士の差で生じるとは限らず、あくまで差の両辺の数が「0に近い」ことが必要だということ。
例えば、
2^(-970)の次の数 - 2^(-970)
= (1.0000000000 0000000000 0000000000 0000000000 0000000000 01)_bin * 2^(-970)
- (1.0000000000 0000000000 0000000000 0000000000 0000000000 00)_bin * 2^(-970)
= (0.0000000000 0000000000 0000000000 0000000000 0000000000 01)_bin * 2^(-970)
= (1.0000000000 0000000000 0000000000 0000000000 0000000000 00)_bin * 2^(-970 - 52)
となって、指数部の値は 1023 + (-970 - 52) = 1 で正規化数としてぎりぎりオッケー。で、もちろん、2^(-970)より大きい数同士の引き算の結果なら、DBL_MINを最小分解能とする正規化数の体系で余裕を持って記述できる。
0近傍のみが特殊であって、0近傍の数同士の差で生じる指数部のアンダーフローを避けるための0への丸めを「ゼロへのフラッシュ」と言うっぽい
前置きが長くなったが、これを避けるために非正規化数というものがある。
DBL_MINの次の浮動小数点数 - DBL_MIN
= (0.0000000000 0000000000 0000000000 0000000000 0000000000 01)_bin * 2^(-1022)
は、無理に正規化数の約束ごとに従わせずに、そのまんま固定小数点的にビット表現すればいいんでねえが、というアイデアだと考えれば多分間違っていない、
いや知らんけど
bNegはtrueまたはfalse、指数部は0、仮数部は非0で、固定小数点的に解釈される。
で、実例
DUMP_DBL()マクロおよびpow2()、GenIEEE754Double()関数は上掲のものを使うとして、コードは約こんな感じ↓↓↓のコードで
/// 特殊な数を分解して確認する。 void experiment01b() { using namespace std; using namespace mylib; // NaNの例 DUMP_DBL(sqrt(-1.0)); DUMP_DBL(-sqrt(-1.0)); cerr << "#" << endl; // 無限大 DUMP_DBL(((1LL << 53) - 1) * pow2(1023 - 52)); // DBL_MAXと一致 DUMP_DBL(((1LL << 53) - 0) * pow2(1023 - 52)); // +INF cerr << "#" << endl; // 無限小 DUMP_DBL(-((1LL << 53) - 1) * pow2(1023 - 52)); // -DBL_MAXと一致 DUMP_DBL(-((1LL << 53) - 0) * pow2(1023 - 52)); // -INF cerr << "#" << endl; // 非正規化数の例 double x1 = GenIEEE754Double(false, 1, 1); // DBL_MINの次の数 double y1 = GenIEEE754Double(false, 1, 0); // DBL_MIN DUMP_DBL(x1 - y1); // これは非正規化数 double x = GenIEEE754Double(false, 53, 1); // 2^(-970)の次の数 double y = GenIEEE754Double(false, 53, 0); // 2^(-970) DUMP_DBL(x - y); // これは非正規化数にならない cerr << "#" << endl; // 0の例 DUMP_DBL(0); DUMP_DBL(-0); // なぜかこう書いても+0と解釈されるが… DUMP_DBL(GenIEEE754Double(false, 0, 0)); // +0 DUMP_DBL(GenIEEE754Double(true, 0, 0)); // -0 }
こんな感じの出力↓↓↓を得る
sqrt(-1.0) = -1.#IND000000000000e+000; bNeg= true, exp=2047 (1023+( 1024)), frac=0x8000000000000 -sqrt(-1.0) = 1.#QNAN00000000000e+000; bNeg=false, exp=2047 (1023+( 1024)), frac=0x8000000000000 # ((1LL << 53) - 1) * pow2(1023 - 52) = 1.7976931348623157e+308; bNeg=false, exp=2046 (1023+( 1023)), frac=0xfffffffffffff ((1LL << 53) - 0) * pow2(1023 - 52) = 1.#INF000000000000e+000; bNeg=false, exp=2047 (1023+( 1024)), frac=0x0000000000000 # -((1LL << 53) - 1) * pow2(1023 - 52) = -1.7976931348623157e+308; bNeg= true, exp=2046 (1023+( 1023)), frac=0xfffffffffffff -((1LL << 53) - 0) * pow2(1023 - 52) = -1.#INF000000000000e+000; bNeg= true, exp=2047 (1023+( 1024)), frac=0x0000000000000 # x1 - y1 = 4.9406564584124654e-324; bNeg=false, exp= 0 (1023+(-1023)), frac=0x0000000000001 x - y = 2.2250738585072014e-308; bNeg=false, exp= 1 (1023+(-1022)), frac=0x0000000000000 # 0 = 0.0000000000000000e+000; bNeg=false, exp= 0 (1023+(-1023)), frac=0x0000000000000 -0 = 0.0000000000000000e+000; bNeg=false, exp= 0 (1023+(-1023)), frac=0x0000000000000 GenIEEE754Double(false, 0, 0) = 0.0000000000000000e+000; bNeg=false, exp= 0 (1023+(-1023)), frac=0x0000000000000 GenIEEE754Double(true, 0, 0) = -0.0000000000000000e+000; bNeg= true, exp= 0 (1023+(-1023)), frac=0x0000000000000 (Ctrl-Z)
*1:本当は眉唾。この「ランダムの場合」という前置きのが曲者で、現実的な数値演算の仮定（有限桁数）の下では、正の数と負の数が同じ頻度で現れないと微妙にバイアスを生じる。つまり一様分布であっても定義域が正の数ばっかりとか、負の数ばっかりだとバイアスが生じる。後で述べるつもり。
*2:という言い方で良いのか？
*3:bNeg) ? -1 : 1) * ((1LL << 52) | frac) * 2^((exp - 1023) - 52) ＜約表す＞であって＜表す＞でないのは、fracの桁数が有限であることによる。例えば演算途中でfracの桁数で表せない端数が生じたとき、例えばfracに入れるべき数を常に符号なしとみなしてLSB未満の切捨てを行うシステムであれば、(bNeg, exp, frac)で特定される1つのIEEE 754形式倍精度浮動小数点数xは、正確には次の範囲の数のどれかを表しﾁｮﾙと見るのが正しいい： ((1LL << 52) | frac) * 2^((exp - 1023) - 52) ≦ x ＜ (((1LL << 52) | frac) + 1) * 2^((exp - 1023) - 52) (bNegがtrueのとき) -(((1LL << 52) | frac) + 1) * 2^((exp - 1023) - 52) ＜ x ≦ -((1LL << 52) | frac) * 2^((exp - 1023) - 52) (bNegがtrueのとき) 以下、ことわりのない限り、浮動小数点数といえばIEEE 754形式倍精度浮動小数点数を指すとする 申し遅れたが、(1LL << 52)をORする演算が毎回出てくるのは、正規化数において仮数部のMSBの次のビットが常に1である（とみなす）という約束があるため。すわなち、仮数部は52 bitしかないが、MSBの次の53 bit目は正規化数においては常に1と解釈される。 で、指数部には、仮数部の仮想的な53 bit目が2のn乗の桁を表すとして、n+1023が入る。よって、仮数部を52 bitの整数とみなすなら、浮動小数点数全体では、上掲の式と同じことだが ( (1LL << 52) | (仮数部
*4:ウチの環境では、std::sqrt(-1.0)だと仮数部が0x8000000000000のNaNが返される。
*5:正規化数の指数部は常に正になるように+1023の下駄を履かせている。負になってはイケナイ
*6:移植性は考えず、IEEE 754形式な浮動小数点数の環境に限定するなら、浮動小数点型の変数をZeroMemory()とかでゼロクリアすることは不定値避けとしては全く正しいと言える。いや知らんけど
