VC++のpow関数でエラー検出

pow(x,y)やると_clearfp()で_SW_INEXACTが出るよねぇ
あれどういう仕組み？
（Power関数のエラーハンドリングについて調べたり試したりして書いてます。
　実際に使う場合はちゃんとテストしてください。）

#include <math.h>
#include <cfloat>

double x = 2.0;
double y = 3.0;
unsigned int rv = _clearfp();
double ans = pow(x, y);
rv = _statusfp();

if ((rv & _SW_INVALID) > 0) {
   //領域エラー xが負で、yが整数でない有限値
} else if ((rv & _SW_ZERODIVIDE) > 0) {
   //極エラー xがゼロで、y が負
} else if ((rv & _SW_UNDERFLOW) > 0 ||
           (rv & _SW_OVERFLOW ) > 0 || !_finite(ans)) {
   //double範囲外エラー
} else {
   //正常に計算終了（丸め誤差発生「rv == _SW_INEXACT」の場合はエラーにしません）
}

このときrvは_SW_INEXACT（0x0001）が設定されます。これはdouble値の計算によって丸め誤差が発生してるのでansの値は正確ではない（かもしれない）よ、という意味です。2の3乗なんて整数の8に決まってますが倍精度浮動小数点数の制限としてビシィっとした「8」ではないよ、という意味です。これはどうしようもないのでrvが_SW_INEXACT単体のエラーであれば無視してansを計算結果として処理し、対応はアプリの使用者に任せるのが良いです。

エラー値はビットで、_SW_INEXACTと_SW_OVERFLOWが同時に設定されたりすることがあるようです。そのため「==」ではなく「&」演算して値が0でないなら、というような判定をしてます（0でない整数値はtrue扱いなので比較は無用なんですが、うーむ）。

この_clearfp()と_statusfp()は簡単に確認した限り、スレッドセーフだと思います。

オーバーフローに関して！

32ビット版のWindowsで動かすとオーバーフローのビットが設定されません。頑張って調べてますが今のところ理由は不明です。オーバーフロー以外は設定されるのでオーバーフローは _finite(x) でもチェックをかけて32ビット版Windowsでも問題が起きないようにしました。それに伴い、オーバーフローのチェックを最後にしました。

---

まとめ！

この手の計算に関して、判明していることを書き残しておきます。
この手の浮動小数の計算では計算精度とエラー検出をコントロールできます。
_controlfp(newVal, mask)を使います。float.hの定義と睨めっこしながら、maskには変更する対象、newValには変更する値を与えます。ビットなのでOR連結で一気に設定も出来ますが個々にやった方が分かり易いです。

unsigned int rv = 0;
_clearfp(); // ステータス値をクリア
unsigned int def = _controlfp(0, 0); // 現在の値を確保
_controlfp(_MCW_EM, _MCW_EM); // 全部のエラーを検出する
_controlfp(_RC_NEAR, _MCW_RC); // 丸め誤差はNEARとする（他DOWN、UP等）
_controlfp(_PC_64, _MCW_PC ); // 制度は64ビットとする（他53、24ビット）

double ans = pow(x, y);
rv = _clearfp(); // 直前のステータス値を取得し、同時にクリアする

_controlfp(def, _MCW_EM | _MCW_RC | _MCW_PC); // 元の値に戻す

ビット演算なのでOR連結で一気に指定も出来ますが個々にやった方が良いでしょう。例えば、検出するエラーを指定する場合は以下のようにやると良いでしょう（でもハンドルしてないエラーになるので止まるかもしれません）。

_controlfp(_EM_INEXACT | _EM_UNDERFLOW | _EM_OVERFLOW, 
             _MCW_EM); // 検出するエラーを指定する

計算精度についてサンプルコードで確認したところ、デフォルトではWin32ビルドでは53ビットが設定されています。x64ビルドでは_PC_64が設定されているように見えますがこれは_PC_64が0x00なのでそう見えるだけで、x64では_MCW_PCに0以外を与えるのはいけないらしいです。つまりx64での計算精度は53ビット固定です。x64のDebugビルドでは_MCW_PCに0以外を設定しようとするとAssertionエラーが起こり、Releaseビルドでも無効なパラメーターの例外が出るそうです。
MSDN（計算精度の違いで結果が違うことを確認するコードが載ってます。「x64 アーキテクチャでは、浮動小数点の精度の変更はサポートされていません」と書いてあります。）

サンプルコードでは計算後に_clearfp()でステータスを取得しています。直前のステータス値を取得しつつ、クリアできるのでこの方がスマートですが、「あぇ？ステータス取得する前にclearやっちゃってん？おぇ？？」みたいになるかもしれないので_statusfp()で取得→_clearfp()でクリア、の方が読み易いんじゃないかなぁとも思います。

（「_controlfp(_EM_OVERFLOW, _EM_OVERFLOW );」とやってもやっぱり32ビットのWindowsではオーバーフローは検出できないようです。「コンパイラが調子悪い」状態ですTT）

さらに「_EM_DENORMAL」というエラーはなんとなく無視していました。これはつまり「非常に小さい値ですよ」という理解で良いんでしょうか？
「計算できなかった」というエラーではなさそう（？）なので、今のとこアンダーフローじゃないならこのエラーは無視してて良いか…
なんか_controlfpでこのエラーを検出しないで良い、と設定しても検出されます。なんででしょう。