2014/03/24(月)Visual Studioのsqrt
#include <float.h>
#include <math.h>
#include <stdio.h>
int main() {
volatile double r, x;
_controlfp(_PC_53, _MCW_PC);
x = 2.;
_controlfp(_RC_UP, _MCW_RC);
r = sqrt(x);
_controlfp(_RC_NEAR, _MCW_RC);
printf("%.17g\n", r);
_controlfp(_RC_DOWN, _MCW_RC);
r = sqrt(x);
_controlfp(_RC_NEAR, _MCW_RC);
printf("%.17g\n", r);
return 0;
}
こんな感じのソースで、なぜか表示される値が同じになってしまう。自分の環境はwin7 64bitで、Visual Studio 2008。
Visual Studio 2008 コマンドプロンプト (Win64じゃない方) でcl使ってコンパイルした。
64bitコンパイルすれば問題ない。
なんでだろ。FPU使わずに独自コード使ったりしてんのかな。
(2014年3月25日追記)
windows用のmatlabの古い奴 (R2007b) の32bit版で同じような現象を確認した。
>> format hex
>> system_dependent('setround', inf)
>> 1/3
ans =
3fd5555555555556
>> sqrt(2)
ans =
3ff6a09e667f3bcd
>> system_dependent('setround', -inf)
>> 1/3
ans =
3fd5555555555555
>> sqrt(2)
ans =
3ff6a09e667f3bcd
除算の丸めは変わっているけどsqrtは変わってない。windows用のmatlabはVisual C++でmakeされているのだろうか。2014/03/19(水)kv-0.4.2
ddでのDBL_MAX/3はちゃんと計算できるようになりました。
他にも、区間演算での無限大の扱いとか、方向付き丸めのエミュレートあたりの完成度が上がっています。
問題がありましたら遠慮なくご意見お願いします。
Rump先生にこのページ知られちゃったので、ドキュメントの英語化をしたくなってきた。が、手間かかるなあ。
2014/03/19(水)dd
z=x/yを計算し、その誤差を拾うためにz*yをtwoproductで計算するのだが、
xがオーバーフロー寸前なのでこの計算がオーバーフローしてしまう。
直すのめんどいな。
ddは最初から近似計算だから許されるとしても、interval + dd + rddは精度保証を謳っているのだから、
(NaNなら嘘を言っているわけでは無いにしても)これはまずい。
2014/02/03(月)ベキ級数演算の解説
みなさんの原稿を見ていると、ベキ級数演算の記述の苦戦が目立ちました。
確かにまとまって書かれた原稿が少なく、困るのも無理もないと思ったので、急いでまとめた資料を作りました。せっかくなので、ここにアップロードしておきます。
psa.pdf
kvライブラリの常微分方程式、数値積分などはこの演算が基礎となっています。書き殴りで完成度はあまり高くないかと思われますが、参考になれば幸いです。
2014/01/25(土)Error Free Transformation (EFT) is NOT error-free
def twosum(a, b) :
x = a + b
tmp = x - a
y = (a - (x - tmp)) + (b - tmp)
return x, y
def split(a) :
tmp = a * (2**27 + 1)
x = tmp - (tmp - a)
y = a - x
return x, y
def twoproduct(a, b) :
x = a * b
a1, a2 = split(a)
b1, b2 = split(b)
y = a2 * b2 - (((x - a1 * b1) - a2 * b1) - a1 * b2)
return x, y
python風に書くとこんな感じ。x,y = twosum(a, b) は、xがa+bを浮動小数点数で普通に計算した値、yはその誤差分になり、x+y=a+bが厳密に成立するとされています。
x,y = twoproduct(a, b) は、xがa*bを浮動小数点数で普通に計算した値、yはその誤差分になり、x+y=a*bが厳密に成立するとされています。
しかし、オーバーフローやアンダーフローに注意が必要です。
xがオーバーフローしたら誤差も何も無いのでまあよしとしましょう。
アンダーフローは、twosumの場合、IEEE754標準のように非正規化数を用いた漸近アンダーフローを備えた浮動小数点数では原理的に発生しないので問題ありません。
しかし、twoproductの場合は、yが非正規化数の領域に入る場合にはyの有効桁数が失われてしまい、無誤差ではなくなってしまいます。すなわち、|x| < 2^(-969) (=-1074+106-1) の場合は誤差が発生します。
また、twoproductで使われているsplitにも注意が必要で、(2^27+1)を乗じている部分でオーバーフローする可能性があります。ただしこれは実装上の工夫でどうにでも回避可能。
ここまではよく知られている問題点ですが、今までにあまり注意を払われて来なかったと思われる問題点を見つけたのでメモしておきます。
まず、twoproductで
a = 6.929001713869936e+236
b = 2.5944475251952003e+71
のケース。x=a*bはオーバーフローしませんが、a1*b1がオーバーフローしてしまい正しくyが計算されません。2つ目は、twosumで、
a = 3.5630624444874539e+307
b = -1.7976931348623157e+308
のケース。x=a+bはオーバーフローしませんが、中間変数tmpがオーバーフローしてしまい正しくyが計算されません。twosum, twoproductを用いた方向付き丸めのエミュレート法を考えていて、これらを見つけました。
精度保証付き数値計算では一つの例外も許されないので、こういう現象には本当に神経を使います。他にまだ気付いていない問題点が存在しないと言い切れるのでしょうか。