kashiの日記

検索条件

日付で絞り込み 1月(1) 3月(3) 4月(1) 6月(4) 8月(1) 9月(1) 10月(4) 11月(4) 12月(3)

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

2015/06/24(水)INTLAB version 9をubuntu上のoctaveで使う

INTLAB version 9の重要性の記事で書いたように、INTLAB version 9を使えばmatlabに頼らずにオープンソースな環境で精度保証付き数値計算を行うことができます。今回、ubuntu 14.04にoctaveをインストールし、その上でINTLABを動かしてみたので、その様子を記録しておきます。

INTLABは、http://www.ti3.tu-harburg.de/rump/intlab/で、50ユーロで購入できます。有償なのが残念ですが、PayPal経由で支払いを済ませると一時間だけ有効なダウンロードリンクがメールで送られてきます。ダウンロードしたファイルは完全にソースコードを見ることができます。また、Linux/Mac/Windows用、あるいはmatlab用/octave用といった区別はなく、すべてこのファイルでOKです。

octaveのインストールは、ubuntu 14.04で、

sudo apt-get install octave
sudo apt-get install octave-doc
sudo apt-get install octave-htmldoc
sudo apt-get install liboctave-dev

のように行いました。特に最後のliboctave-devが必要なことに注意してください。また、ubuntu 14.04でBLASを使うの記事のようにATLASやopenblasを導入すれば、octaveから使われるBLASがそれに置き換わります。openblasを入れることをお勧めします。

ダウンロードした INTLAB.zip を任意のdirectoryに展開します。例えば ~/Intlab_V9/ としましょう。まず、octaveを起動し、必要ならば cd Intlab_V9 などとしてINTLAB directoryに移動し、

startintlab

と実行します。これは、そこのdirectoryにある startintlab.m を実行しています。うまく行けば、これだけでintlabが使える状態になる筈です。初回起動時は、裏で

Intlab_V9/setround/octave/setround.mex
Intlab_V9/INTLAB_Data_Intlab_Version_9_x86_64-pc-linux-gnu3_8_1.mat

の2つのファイルが作られました。1つめは丸めの向きを変えるためのファイルで、2つめは数学関数の誤差を調査し記録しておくためのファイルです。1つめのファイルを作成するために、先に入れた liboctave-dev が必要になります。C++コンパイラもインストールされている必要があると思います。2つめのファイルは内部にdirectory情報を含んでいるらしく、intlab一式の場所を移してしまうと動かなくなります。その場合はそのファイルを削除して再度 startintlab を実行すれば大丈夫です。

実際に区間演算をしてみます。T先生からの挑戦状で計算した例題です。

octave:1> format long
octave:2> x=intval(1);
octave:3> y=intval(1);
octave:4> for i=2:10; z=(1+2*y)/(x*y^2); x=y; y=z; i, z, end
i =  2
intval z = 
   3.00000000000000
i =  3
intval z = 
   0.77777777777777
i =  4
intval z = 
   1.40816326530612
i =  5
intval z = 
   2.47448015122873
i =  6
intval z = 
   0.68995395922102
i =  7
intval z = 
   2.0203934747424_
i =  8
intval z = 
   1.7898074714308_
i =  9
intval z = 
   0.7075878600520_
i =  10
intval z = 
   2.69514211796___
octave:5>

のように計算出来ました。最初の intval を affari に変更すればaffine arithmeticで計算させることもできます。

また、行列積を計算させてみると、

octave:1> a=rand(5000);b=rand(5000);
octave:2> tic;a*b;toc
Elapsed time is 3.20024 seconds.
octave:3> tic;intval(a)*b;toc
Elapsed time is 8.66678 seconds.

くらいでした。core i7 4770Kですがvmwareを挟んでいるので生で使うと更に速いでしょう。

なお、octave起動時に毎回startintlabを実行させるのは面倒なので常時実行させるには、

/usr/share/octave/site/m/startup/octaverc
home directoryの.octaverc
current directoryの.octaverc

のいずれかのファイルに、

addpath('/intlabを入れたパス/Intlab_V9')
startintlab

と書いておけば良さそうです。

余談ですが、octave起動時に

octave --force-gui

とすると、

のような感じで起動します。変数一覧とか見られて便利そう。

openblasで丸め変更に問題はないかとか、検証する必要はあるかと思います。

自分はmatlab嫌いなものでINTLABの使い方はよく分かっていませんが、お楽しみ下さい。

コメント（0件）

2015/06/06(土)kv-0.4.21

精度保証

kvライブラリを0.4.21にアップデートしました。

今回は、一応今までもひっそりと機能はあった、端点特異性を持つ関数の数値積分のコードを大幅に書き直し、機能も追加し、また詳細なwebドキュメントを書きました。端点特異性を持つ関数のうち、

除去可能特異点を持つもの。例えば $\int_0^1 \frac{\sin(x)}{x} dx$
ベキ型特異点を持つもので、x^pg(x)の形のもの。例えば $\int_0^1 \sqrt{x}\cos(x)dx$
ベキ型特異点を持つもので、f(x)^pの形のもの。例えば $\int_0^1 \sqrt{\sin(x)}dx$
ベキ型特異点を持つもので、f(x)^pg(x)の形のもの。例えば $\int_0^1 \sqrt{\sin(x)}\cos(x)dx$
log型特異点を持つもので、log(x)g(x)の形のもの。例えば $\int_0^1 \frac{\log(x)}{x+1}dx$
log型特異点を持つもので、log(f(x))の形のもの。例えば $\int_0^1 \log(\sin(x))dx$
log型特異点を持つもので、log(f(x))g(x)の形のもの。例えば $\int_0^1 \log(\sin(x))\cos(x)dx$

などを精度保証付きで計算できるようになりました。

その他、conj(共役複素数)を追加したり、全解探索における反復改良が稀にうまくいかなかったのを直したりしました。

最近少しずつユーザが増えてきて、最新のmacのclangで計算値が正しくないことがありそうとか、win7 32bit環境とVC++2012でddの表示がうまくいかないとか、報告を受けています。実機が用意できなかったりこちらでうまく再現しなかったりしてなかなか直せませんが、極力対応していきたいと考えていますので、こちらの環境(ubuntu 64bit)と異なる環境での動作報告、お待ちしております。

コメント（0件）

2015/04/16(木)kv-0.4.20

精度保証

kvライブラリを0.4.20にアップデートしました。

今回は、主に端点に無限大を持つような区間に関するbug-fixです。区間演算ライブラリは、端点に無限大を持たせることによって半無限区間や全区間を取り扱えるように設計されています。区間演算の実装は、上端下端型と、中心半径型の2種類が考えられますが、本ライブラリでは半無限区間や全区間を扱えるという上端下端型のメリットを重視し、上端下端型を採用しています。

端点に無限大が入ると、無限大-無限大、無限大/無限大、無限大*0などでNaNが発生する危険性があり、これらを注意深く排除する必要があります。今回、乗算において思慮が浅かった部分が発見され、きちんと考えなおした上で修正しました。詳細なアルゴリズムは端点に無限大を持つ区間演算まとめを見て下さい。無限大を扱わない場合は問題ありませんが、無限大を扱った計算において、0.4.19以前のライブラリの計算結果は精度保証されていなかった可能性があります。0.4.20へのアップデートをお勧めします。

また、sinh, coshに-無限大を入れた場合の不具合も修正しました。

コメント（0件）

2015/03/23(月)kv-0.4.19

精度保証

kvライブラリを0.4.19にアップデートしました。

主な変更点は、

mpfrである程度以上仮数部の長い数を扱ったとき、e, ln2, piの定数の精度が不十分なため区間幅が小さくならない問題を修正。
atanhの端点処理をちゃんとやるようにした。
ddとmpfrにceilとldexpを追加。

です。特に最初のが重要な修正です。

区間演算では、テンプレートを使って、interval<double>ならdoubleにふさわしい精度で、interval<dd>ならdouble-doubleにふさわしい精度で数学関数の計算を行う仕組みになっています。この仕組みで困るのが定数の計算で、定数は入力を持たないのでどの精度で計算すべきか判断できません。そこで、(numeric_limits<T>の使い勝手に倣って) constants<T>::pi(), constants<T>::e(), constants<T>::ln2()を用意し、Tの型によって計算精度を変える仕組みを採用しています。で、これらの実装は、例えばpiなら、

template <class T> struct constants< interval<T> > {
    static interval<T> pi() {
        static const interval<T> tmp(
            "3.1415926535897932384626433832795028841971693993751",
            "3.1415926535897932384626433832795028841971693993752"
        );
        return tmp;
    }
}

のように文字列で上端下端を与えそれを上下に丸めることによって得ていました。この定数は10進数50桁で円周率の上界下界になっています。doubleとddしか無ければこれだけ精度があれば十分だったのですが、mpfrで50桁以上の精度を指定した場合、(精度保証は出来ているものの)定数の区間幅が広すぎてまともな精度が出ない事態になっていました。ここの仕組みを変えて、きちんとその型なりの定数の計算方法を書くようにしました。

コメント（0件）

2015/03/12(木)kv-0.4.18

精度保証

前回からほとんど時間が経っていませんが、kv-0.4.18にアップデートしました。

0.4.17で、lgammaにdouble以外の型を入れるとコンパイルできないバグを混ぜてしまっていたので修正しました。いつもならもう少し修正が貯まってから公開するのですが、某所でlgammaが必要だったのですぐに公開しました。ついでに、第1種Bessel関数(整数次のみ)の簡単な実装を追加しました。

コメント（0件）