kashiの日記

検索条件

全4件 (1/1ページ) 記事本文を表示

2014/03/30(日)Cardanoの公式とFerrariの公式

3次方程式の解法であるCardanoの公式、4次方程式の解法であるFerrariの公式、
数学の歴史で大変有名だけど、中身をよく吟味したことは無かった。
ちょっとしたきっかけで、これらをプログラムする機会があったのだが、
案外ややこしくて苦労した。Cardanoでは、複素数の3乗根を取る部分で3種類の解のうち
適切な1つを選ばないといけない部分があるのだが、こういうところがきちんと
書いてない文献が多かった。Ferrariでも平方根で同じような箇所がある。

最初はmathematica先生にSolve[a4 x^4 + a3 x^3 + a2 x^2 + a1 x + a0 = 0, x]とか投げて
出てきた式を適当に簡約化してこんなんでいいか、なんてやってたが、結果的にそれじゃ
使い物にならなかった。

CardanoやFerrariでは反復によらず直接解を計算するので、全計算過程を区間演算化すれば
容易に精度保証付きの解が得られる。と思って書いてみたが、ギリギリの計算が要求されて
ライブラリの問題点がいろいろ洗い出される結果に。

atan2([0, 1e-10], [-1, 1])のような場合に結果の区間幅が大きくなってしまう。
expに無限大に近いが無限大ではない数を入れたときの結果がおかしい (引数をfloorで丸めたものがintの範囲を超える場合)
複素数の絶対値や平方根で(特に原点を含むような引数で)エラーになることがある。
ddの平方根に0を入れるとnanになる(newton法の失敗)
exp(interval<dd>(-1e308)*2)のような計算で無限ループに陥る

などなど。
複素数とddは実戦投入の経験が少なかったので、ちょうどいい訓練になった。
複素数の0に対してexp(log(0))がちゃんと0(あるいは0を含む微小区間)を返さなければ
うまく動かないようなシチュエーションも出てきて、intlabすらここはきちんと動かなかった。

最終的にはなんとか動くようになったのだが、CardanoはまだマシだがFerrariは素晴らしく
精度が出ない。doubleで平均的には7,8桁くらいか。
計算時間もかかるし、これらの方法が実際の数値計算法としてあまり使われない理由がよく分かった。
何でもやってみるものだなあ。

おかげでライブラリも結構アップデートしたので、ファイルの整理が出来たらversion upしよう。

コメント（0件）

2014/03/24(月)Visual Studioのsqrt

精度保証

Visual Studioの32bitで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されているのだろうか。

コメント（2件）

2014/03/19(水)kv-0.4.2

精度保証

というわけで、変な挙動を放置するのも気持ち悪いので、0.4.2にアップデートしました。
ddでのDBL_MAX/3はちゃんと計算できるようになりました。
他にも、区間演算での無限大の扱いとか、方向付き丸めのエミュレートあたりの完成度が上がっています。

問題がありましたら遠慮なくご意見お願いします。

Rump先生にこのページ知られちゃったので、ドキュメントの英語化をしたくなってきた。が、手間かかるなあ。

コメント（0件）

2014/03/19(水)dd

精度保証

ddで、DBL_MAX/3がNaNになる現象を発見。dd計算はまだ細かいところが甘いな。
z=x/yを計算し、その誤差を拾うためにz*yをtwoproductで計算するのだが、
xがオーバーフロー寸前なのでこの計算がオーバーフローしてしまう。
直すのめんどいな。

ddは最初から近似計算だから許されるとしても、interval + dd + rddは精度保証を謳っているのだから、
(NaNなら嘘を言っているわけでは無いにしても)これはまずい。

コメント（0件）