2015/10/11(日)数学関数の精度を調べる

精度保証
ふと思い立って、C言語の(というかlibmの)数学関数の精度を調べてみました。

現代の一般的なプログラミング環境では、浮動小数点演算はIEEE 754 Std.に基づいて行われるのが一般的です。そしてIEEE 754 Std.においては、演算に完全精度 (フォーマットで表現可能な範囲内で最善の精度、0.5ulp) を要求しているのは
  • 加算
  • 減算
  • 乗算
  • 除算
  • 平方根
の5種類の演算のみで、それ以外の演算の精度については何も規定されていません。つまり、どんな値が計算されても文句は言えません。

従って、C言語などから使うexp, sinなどの数学関数の返す値は信用することは出来ません。kvライブラリの区間演算でも数学関数は全て自力で実装しています。区間演算ライブラリの一つである PROFIL/BIAS では数学関数をC言語の組み込み関数で計算した後succ/predで区間を広げて計算結果としていますが、これは誤りであり、精度保証付き数値計算を行うためのツールとしては役に立ちません。

Intelの64bit環境では、数学関数は全てソフトウェアで計算されるのが一般的です。従って、計算精度はlibmの実装に依存することになります。たとえ完全精度を出す義務が無かったとしても、長い間のアルゴリズムの改良により完全精度が出るようになっているかも知れません。そこで、
  • CPU: Intel Core M 5Y70
  • OS: Ubuntu 14.04 LTS 64bit
という環境で、数学関数が実際にどのくらいの精度が出ているのか調べてみました。比較対象は、kvライブラリのinterval<dd>の計算結果をinterval<double>に直したもので、その区間の外側の値が計算されたら、誤差が1ulpより大きかったことになります。下に示したような入力範囲で100万個の乱数を発生させて調べました。結果を示します。
関数名下限上限結果
exp-700700OK
log1e-3081e308OK
sin-1e101e10OK
cos-1e101e10OK
tan-1e101e10OK
asin-11OK
acos-11OK
atan-1010OK
sinh-700700NG
cosh-700700OK
tanh-700700NG
asinh-10001000OK
acosh11000OK
atanh-11NG
このように、sinh, tanh, atanhで精度不足が見つかりました。言うまでも無いことですが、OKとなった他の関数について、このテストで精度不足が見つからなかったからと言って精度に問題が無いことの証明にはなりません。また、このテストでは1ulp以上の誤差を検出しますが、一般に完全精度と言えば0.5ulp以下の精度のことを指します。

sinh, tanh, atanhについて、精度不足となった例を示しておきます。

sinh(0.74251904152333736419677734375)=0.8126541801799147535234624228905886...
と計算されたが、真の値はdoubleの区間
[0.8126541801799148645457648854062426,0.8126541801799149755680673479218968]
の間にある。

tanh(0.48066521994769573211669921875)=0.4467762045341044374602290645270841...
と計算されたが、真の値はdoubleの区間
[0.4467762045341044929713802957849111,0.4467762045341045484825315270427382]
の間にある。

atanh(0.2377787786535918712615966796875)=0.2424184572356452571639806592429522...
と計算されたが、真の値はdoubleの区間
[0.2424184572356452016528294279851252,0.2424184572356452294084050436140388]
の間にある。

もう記録は残っていませんが、10年以上前に調べたときはもっと精度に問題があった記憶があるので、libmの品質は上がっているのかもしれません。案外精度良かったんだな、という印象です。
OK キャンセル 確認 その他