| [Top] | [Contents] | [Index] | [ ? ] |
| 1. BNCとは？ | ||
| 2. BNCpackのインストール | BNCのインストール | |
| 3. 機能一覧 | ||
| 4. BNCpackを使ったサンプルプログラム集 | テストプログラム | |
| 5. 今後の予定(あてにしないように) | ||
| 謝辞・参考文献 | ||
| 索引 |
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
BNCpackはIEEE単精度・倍精度浮動小数点数(float, double)，及びGNU MP(GMP)(1)及びGNU MPFR(2)の多倍長浮動小数点数をサポートする基本数値計算ライブラリです。全てANSI C(C++にあらず)で記述されています。
元々は，自分が使って来たCプログラムが貯まってきたので整理しよう，という動機で始めたものです。が，近年のPCの著しい機能向上とコスト低下のおかげで多倍長計算もお手軽な環境下で使い物になってきた，じゃあ，そいつも組み込んでしまえということで，やっつけ仕事で作ってしまった代物です。
以前，私は旧労働省管轄下の特殊法人職員で，社会人向けの有料セミナーの講師をよく勤めていました。プログラミング関係で需要が多かったのは，Internet絡みのものとVisual Basic(3)関係でした。それで少しVisual Basicと遊んでみましたが・・・。使い勝手は良いものの，実行プログラムのスピードが遅く，数値計算用言語としてはあまりお勧めできないことがわかりました。そこでnative Cで書いたライブラリをDLLにして，それを売り物にしようと目論んだこともありましたが，程なく転職してしまいましたので，Visual Basicに義理立てする必要がなくなり，その計画は頓挫しております。もっとも，BNCはGMPを切り離して構築することも可能ですので，IEEE754浮動小数点数演算のみサポートするDLLにするのはそれほど困難ではないと思います(4)。
今回公開するものは，GMPをベースとした多倍長計算サポートというのが一番の目玉(5)です。GMPが主としてUNIX環境をターゲットにしている関係上，本ライブラリもUNIX環境で開発してきました。
具体的なBNCpackの利用方法については機能一覧，もしくはサンプルプログラムの章を参照して下さい。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
最初に述べたように，BNCpackはGMPを利用した多倍長計算が売りのライブラリですが，IEEE754 standard(単精度，倍精度)で動作する関数も用意されています。頭文字を除いて同じ名前の関数が3つ並んでいたら，全く同じアルゴリズムを単精度(f/F)，倍精度(d/D)，多倍長精度(mpf/MPF)で実行するものだと思って間違いないでしょう。
従来多倍長計算はGMPのmpf_t型を利用してきましたが，Version 0.3からはMPFRパッケージも利用することが出来るように改良されました。従って，Version 0.3以降のBNCpackは
という3種類の環境で使用可能です。但し，MPFRパッケージを使う場合，現状ではmpf_t型及びそれを使用する関数群を，全てmpfr_t型とそれを使用する関数群に置き換える(mpf2mpfr.h(6)を利用する)ことで凌いでいますので，基本的にはmpf_t型とmpfr_t型は共存できません(7)。
以下，この3つの環境でBNCpackをインストールする方法を簡単に示します。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
`Makefile'をご自分の環境に合わせて設定した後，
% make |
として，`libnc.a'を作成して下さい。後は`bnc.h'とこの`libbnc.a'を適当なディレクトリにコピーしてお使い下さい。
当然のことながら，このようにして作成された`libbnc.a'は以下で述べる多倍長計算の機能は一切使用できません。
後は，自分のプログラムが`myprog.c'であれば，適切な場所で`bnc.h'をインクルードしておき
% cc -omyprog myprog.c -lbnc -lm |
あるいはダイレクトに
% cc -omyprog myprog.c /libdir/libbnc.a -lm |
としてコンパイル・リンクして使用して下さい。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
GMPのインストール方法等についてはhttp://gmplib.org/をご参照下さい。特段変わった環境でなければ，Version 3.0.1以降については殆どのUNIX環境下でするようです。私の使用しているRPM系のLinuxやPlamo Linux/Slackware Linuxでは
% ./configure; make; su # make install |
でコンパイルとインストールを行うことが出来ています。
GMPがインストールしてあれば，`bnc-x.x.x.tar.gz'をダウンロード後，`Makefile'のコンパイルオプション等を適切に設定した後，
% make |
で，`libbnc.a'が出来上がるはずです。あとは`bnc.h'とこの`libbnc.a'を適当なディレクトリにコピーしてお使い下さい。
自分で作ったプログラムが`myprog.c'であれば，適切な場所で`bnc.h'をインクルードしておき
% cc -omyprog -DUSE_GMP myprog.c -lbnc -lgmp -lm |
あるいはダイレクトに
% cc -omyprog -DUSE_GMP myprog.c /libdir/libbnc.a -lgmp -lm |
としてコンパイル・リンクして下さい。この-DUSE_GMPというオプションがうっとうしいときには，直接プログラム中で
#define USE_GMP #include "bnc.h" |
として下さい。`bnc.h'をインクルードする前に定義しておく必要があります。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
MPFRパッケージ(http://www.mpfr.org/)はGMPをベースに，IEEE754浮動小数点数と互換性を持つように改良された多倍長浮動小数点演算ライブラリです。GMPのmpf_t型との違いは
という点にあります。MPFRパッケージは数値計算用途の機能を拡充したものと言えるでしょう。
mpf_t型と同様，`bnc-x.x.x.tar.gz'をダウンロード後，`Makefile'のコンパイルオプション等を適切に設定した後，
% make |
で，`libbnc.a'が出来上がるはずです。あとは`bnc.h'とこの`libbnc.a'を適当なディレクトリにコピーしてお使い下さい。
MPFRパッケージと共に，自分で作ったプログラム`myprog.c'をコンパイルするのであれば，適切な場所で`bnc.h'をインクルードしておき
% cc -omyprog -DUSE_GMP -DUSE_MPFR myprog.c -lbnc -lmpfr -lgmp -lm |
あるいはダイレクトに
% cc -omyprog -DUSE_GMP -DUSE_MPFR myprog.c /libdir/libbnc.a \ -lmpfr -lgmp -lm |
としてコンパイル・リンクして下さい。この-DUSE_GMP -DUSE_MPFRというオプションがうっとうしいときには，直接プログラム中で
#define USE_GMP #define USE_MPFR #include "bnc.h" |
として下さい。`bnc.h'をインクルードする前に定義しておく必要があります。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
| 3.1 基本データ型定義 | ||
| 3.2 基本関数 | ||
| 3.3 スタック | ||
| 3.4 複素数 | ||
| 3.5 配列 | ||
| 3.6 多項式 | ||
| 3.7 基本線型計算 | ||
| 3.8 LU分解 | ||
| 3.9 Conjugate-Gradient 法 | ||
| 3.10 非線型方程式 | Newton法 | |
| 3.11 常微分方程式の初期値問題(単段法) | ||
| 3.12 DKA法 | ||
| 3.13 疎行列 |
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
構造体には，その構造体ごとの精度を保持する変数unsigned long precを含みます。それ以外の構成はIEEE754浮動小数点数を利用したものと変わりません。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
上から順に，IEEE754単精度浮動小数点数による複素数型，倍精度浮動小数点数による複素数型，多倍長浮動小数点数による複素数型です。
この複素数型は次のように定義されています。
typedef struct{
float re; /* Real part */
float im; /* Imaginary part */
} fcmplx;
|
これらの複素数型を利用した関数群としては，DKA法が存在します。線型計算ではまだ利用できません。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
上から順に，単精度浮動小数点数による多項式型，倍精度係数の多項式型，多倍長係数の多項式型です。
多項式型は次のように定義されています。
typedef struct{
float *coef; /* Coefficients of polynomial */
long deg; /* Degree of polynomial */
} fpoly;
typedef fpoly *FPoly;
|
これら多項式型を利用した関数群は，DKA法に利用されています。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
上から順に，単精度浮動小数点数によるベクトル型，倍精度ベクトル型，多倍長ベクトル型です。
上から順に，単精度浮動小数点数の複素ベクトル型，倍精度複素ベクトル型，多倍長複素ベクトル型です。
ベクトル型は次のように定義されています。原則としては構造体へのポインタである[F,D,MPF,CF,CD,CMPF]Vectorだけを使います。
typedef struct{
float *element; /* Elements of vector */
long dim; /* Dimension of vector */
} fvector;
typedef fvector *FVector;
|
複素ベクトル型を使った関数群はまだ提供されていません。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
上から順に，単精度一般行列型，倍精度一般行列型，多倍長一般行列型です。
上から順に，単精度複素一般行列型，倍精度複素一般行列型，多倍長複素一般行列型です。
一般行列型は次のように定義されています。原則としては構造体へのポインタである[F,D,MPF,CF,CD,CMPF]Matrixだけを使います。
typedef struct {
float *element; /* Elements */
long row_dim, col_dim; /* Rows, Columns */
} fmatrix;
typedef fmatrix *FMatrix;
|
複素一般行列型を扱う関数群はまだ提供されていません。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
上から順に，単精度スタック，倍精度スタック，多倍長スタックです。
スタックは次のように定義されています。原則としては構造体へのポインタである[F,D,MPF]Stackだけを使います。
typedef struct {
float *array; /* stack */
long size; /* Height of stack */
long index; /* pointer to top of stack */
} fstack;
typedef fstack *FStack;
|
スタックを利用した関数群はまだ提供されていません。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
上から順に，単精度，倍精度，多倍長，複素単精度，複素倍精度，複素多倍長配列です。
配列は次のように定義されています。原則としては構造体へのポインタである[F,D,MPF]Arrayだけを使います。
typedef struct {
float *array; /* array */
long int size; /* Length of array */
} farray;
typedef fstack *FArray;
|
配列は，DKA法の関数群で利用されています。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
上から順に，倍精度，多倍長精度の疎行列用データ型です。
疎行列は次のような構造体で定義されています。
typedef struct {
unsigned long prec;
mpf_t *element; // Elements of matrix
long int row_dim, col_dim; // Dimensions of Row and Column
long int **nzero_index; // Indeces of Non-zero elements
long int *nzero_col_dim; // Numbers of non-zero elements in i-th row
long int *nzero_row_dim; // Numbers of non-zero elements in i-th column
long int nzero_total_num; // Total number of non-zero elements
} mpfrsmatrix;
typedef mpfrsmatrix *MPFRSMatrix;
|
疎行列は次のように格納されます。
0 1 2 3 4
A = [a 0 b c 0]0
[0 d 0 0 0]1
[0 e f 0 0]2
[0 0 0 g 0]3
[0 0 h 0 i]4
<--> element = [a b c d e f g h i]
row_dim = 5, col_dim = 5
nzero_index[0] = [0 2 3]
nzero_index[1] = [4]
nzero_index[2] = [1 2]
nzero_index[3] = [3]
nzero_index[4] = [2 4]
nzero_col_dim[0] = 3
nzero_col_dim[1] = 1
nzero_col_dim[2] = 2
nzero_col_dim[3] = 1
nzero_col_dim[4] = 2
nzero_total_num = 9
nzero_row_dim[0] = 1
nzero_row_dim[1] = 2
nzero_row_dim[2] = 3
nzero_row_dim[3] = 2
nzero_row_dim[4] = 1
|
但し，現状では格納段階ではnzero_row_dim[]はセットされません。設定する際にはset_nzero_row_dim関数を呼び出して下さい。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
ropとopを比較して最大値，最小値を返します。
rop^{opを計算します。
BNCpackのmpf_...関数群で使用するデフォルトの精度(仮数部のbit数)をセットします。この関数を呼び出した後のmpf_...関数が影響を受けます。MPFRパッケージを用いる場合はここでデフォルトの丸めモード(RN)もセットします。
MPFRパッケージを読み込んでいる場合のみ，有効となる関数です。 BNCpackのmpf_...関数群で使用するデフォルトの丸めモードをセットします。使用できる丸めモードは次の４種類です。
この関数を呼び出した後のmpf_...関数が影響を受けます。
BNCpackのmpf_...関数群で使用するデフォルトの精度を返します。
円周率と自然対数の底の値を返します。MPFRパッケージを使用する場合は，そこで定義されている初等関数を使用することになります。
床関数です。GMP Ver.3.xよりサポートされましたが，Ver.2.x以前でも有効な関数として作成してあります。
mpf_t型を単精度，倍精度に変換します。
上から順に，正弦関数，余弦関数，指数関数，対数関数，常用対数関数です。Maclaurin展開に基づくルーチンですので，ものすごく速度が遅いです。高速な多倍長初等関数群が必要でしたらmpfrパッケージの使用をご検討下さい。MPFRパッケージを使うことで，これらの関数は全て高速なものに置き換えられます。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
スタックを初期化します。多倍長スタックはset_bnc_default_prec関数で定義された精度で初期化されます。
多倍長スタックを精度指定で初期化します。
スタックを消去します。
値をスタックに積みます。
スタックから値を取り出します。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
複素数を格納する領域を確保します。
prec bitの精度を持つ複素数を格納する領域を確保します。
複素数opの格納されている領域を解放します。
複素数opの実数部を返します。多倍長型はropに実数部が格納されます。
複素数opの虚数部を返します。多倍長型はropに虚数部が格納されます。
実数valを複素数ropの実数部に格納します。
整数valを複素数ropの実数部に格納します。
実数valを複素数ropの虚数部に格納します。
符号なし長整数valを複素数ropの虚数部に格納します。
複素数opを複素数ropに代入します。
複素数ropに0をセットします。
複素数op1と複素数op2の和を複素数ropに格納します。
複素数opに実数valを加え，複素数ropに格納します。
複素数ropに複素数opを加え，結果をropに戻します。
複素数同士の差op1 - op2を計算し，結果をropに格納します。
複素数opと共約な複素数を計算し，結果をropに返します。
複素数op全体の符号を反転し，結果をropに格納します。
複素数op全体の符号を反転し，結果をropに格納します。
複素数opの絶対値を計算して返します。
複素数同士の積op1*op2を計算し，結果をropに格納します。
複素数opと実数valとの積をropに格納します。
複素数opと整数valとの積をropに格納します。
複素数同士の積rop\times \var{opを計算し，結果をropに戻します。
\exp(\var{op\times i)を計算し，複素数ropに格納します。
複素数opの値を標準出力に表示します。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
array_size個のデータが格納できる配列を確保します。
array_size個の，prec bitの精度を持つデータを格納できる配列を確保します。
配列ropの記憶領域を削除します。
配列opのindex番目のデータを返します。
配列ropのindex番目にデータvalを格納します。
配列opの内容を配列ropに代入します。
配列opの内容を標準出力に表示します。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
最大次数degreeの多項式を初期化します。係数が格納される配列の要素数がdegreeで決定されます。
係数がprecbitの精度を持つ最大次数degreeの多項式を初期化します。
多項式polyの記憶領域を解放します。
多項式polyのindex次係数を返します。
多項式polyの次数を返します。
多項式polyのindex次係数に実数valを格納します。
多項式polyのindex次係数に倍精度実数valを格納します。
多項式polyのindex次係数に，文字列で表現されたbase進数の実数valを格納します。
多項式polyの係数を標準出力に表示します。
単精度多項式fpoly, 倍精度多項式dpoly, 多倍長多項式mpfpolyの係数を標準出力に表示します。
2つの多項式poly1とpoly2の和をrpolyに格納します。
2つの多項式poly1とpoly2の差をrpolyに格納します。
多項式polyと実数valとのスカラー積をrpolyに格納します。
多項式polyを多項式rpolyに代入します。
多項式polyの係数を全て0にセットします。
多項式polyの絶対値最大係数の次数を返します。
多項式polyの非零係数の個数を返します。
多項式polyの導関数を計算し，polyに導関数を格納します。
多項式polyの一変数が実数valである時の値を計算して返します。多倍長多項式の場合はevalにその値を格納します。
多項式polyの一変数が実数valである時，その導関数の値を計算して返します。多倍長多項式の場合はevalにその値を格納します。
多項式polyの一変数が複素数valである時の値を計算し，evalにその値を格納します。
多項式polyの一変数が複素数valである時，その導関数の値を計算し，evalに格納します。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
ベクトルを初期化します。多倍長ベクトルはset_bnc_default_prec関数で定義された精度で初期化されます。
多倍長ベクトルを精度指定で初期化します。
ベクトルを解放します。
ベクトルvecのindex番目の要素を取り出します。下のget_...関数を短い名前でマクロにしたものです。
ベクトルvecのindex番目の要素を取り出します。
ベクトルvecのindex番目の要素に値valを代入します。下のset_...関数を短い名前でマクロにしたものです。
ベクトルvecのindex番目の要素に値valを代入します。
多倍長ベクトルの精度を返します。
多倍長ベクトルの要素の中で，最小の精度と最大の精度をそれぞれ返します。
ベクトルを標準出力に書き出します。
行列を初期化します。多倍長行列は，set_bnc_default_prec関数で設定された精度で初期化されます。
多倍長行列を精度指定で初期化します。
行列を解放します。
行列matのrow_index, col_index番目の要素を取り出します。下のget_...関数を短い名前でマクロにしたものです。
行列matのrow_index, col_index番目の要素を取り出します。
行列matのrow_index, col_index番目の要素に値valを代入します。下のset_...関数を短い名前でマクロにしたものです。
行列matのrow_index, col_index番目の要素に値valを代入します。
多倍長行列の精度を返します。
多倍長行列成分の最小精度と最大精度をそれぞれ返します。
行列を標準出力に表示します。
ベクトルの和を計算します。
ベクトルの差を計算します。
rvec + vecを計算し，結果をrvecに代入します。
rvec - vecを計算し，結果をrvecに代入します。
vecのスカラーopt倍を計算します。
rvecのスカラーopt倍を計算し，結果をrvecに代入します。
内積を計算します。
vecのノルムをそれぞれ計算します。
ベクトルvecをベクトルrvecに代入します。
行列の和を計算します。
mat1 - mat2を計算し，rmatに代入します。
mat1 * mat2を計算し，rmatに代入します。
行列matを行列rmatに代入します。
行列rmatに零行列を代入します。
行列rmatを単位行列にします。
行列matとベクトルvecとの積を計算し，結果をベクトルrvecに代入します。
行列matの逆行列を計算します。matは正方行列でなければなりません。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
行列matのLU分解を行います。
LU分解された行列matと定数ベクトルvecを使い，後退代入を行って解ベクトルをrvecに代入していきます。
行列matの部分ピボット選択付きLU分解を行います。ピボット選択後の行の状態は配列iarrayに保存されます。
部分ピボット選択付きでLU分解された行列matと定数ベクトルvecを使い，後退代入を行って解ベクトルをrvecに代入していきます。
行列matの完全ピボット選択付きLU分解を行います。ピボット選択後の行の状態は配列row_iarrayに，列の状態は配列col_iarrayに保存されます。
完全ピボット選択付きでLU分解された行列matと定数ベクトルvecを使い，後退代入を行って解ベクトルをrvecに代入していきます。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
Conjugate-Gradient法を実行します。係数行列をmatに，定数ベクトルをvecに入れておき，残差のノルムを使った停止条件をrepsとaepsに指定します。最大反復回数はmaxtimesで指定します。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
多次元のNewton法を実行します。
1次元のNewton法を実行します。
多次元の簡易Newton法を実行します。
1次元の簡易Newton法を実行します。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
常微分方程式の初期値問題を単段法で解くための関数を呼び出す時には必ずbnc.hに加えてbncode.hを呼び出して下さい。
#include "bnc.h" #include "bncode.h" |
加えて，libbncode.aもリンクして下さい。
常微分方程式の初期値問題は
dy/dx = f(x, y) <---> void f(Vector ret_val, x, Vector y) y(x0) = y0 |
という形で与えられるものとします。陰的解法で必要となるJacobi行列は
df/dx <---> void df(Matrix ret_val, x, Vector y) |
という形で与えられるものとします。中間出力が必要な複数ステップの計算の途中結果はファイルポインタfpで指定された出力先に表示されます。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
陽的Euler法の1ステップ分の計算を行います。
陽的Euler法と固定ステップh (= (x - x0) / div_num)を使って$y(x)$の近似値yを計算します。
陽的Euler法と補外法に基づく局所誤差制御(||y_next - y_old || <= rtol * ||y_old|| + atol)を用いてy(x)の近似値yを計算します。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
陽的Runge-Kutta(RK)法の係数を与えます。coef_typeには
が使用可能です。
陽的RK法の1ステップ分の計算を行います。
陽的RK法と固定ステップh (= (x - x0) / div_num)を使ってy(x)の近似値yを計算します。
陽的RK法と補外法に基づく局所誤差制御(||y_next - y_old || <= rtol * ||y_old|| + atol)を用いてy(x)の近似値yを計算します。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
陰的RK法の係数を与えます。coef_typeには
が使用可能です。
陰的RK法で使用する一時変数を初期化します。
陰的RK法で使用する一時変数を消去します。
陰的RK法の1ステップ分の計算を行います。maxtimesは内部反復(簡易Newton法)の最大反復回数の指定です。
陰的RK法と固定ステップh (= (x - x0) / div_num)を使ってy(x)の近似値yを計算します。
陰的RK法と補外法に基づく局所誤差制御(||y_next - y_old || <= rtol * ||y_old|| + atol)を用いてy(x)の近似値yを計算します。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
調和級数(～_harmonic)とRomberg級数(～_nim)に基づく補外法で任意段数(次数)の近似値を計算します。
void init_mpf_ex_nim(long int stage, long int dimension, unsigned long prec) 補外法の計算に必要な一時変数を宣言します。stageは最大段数の設定です。
補外法の計算に使用した一時変数を消去します
補外法の1ステップ分の計算を行います。局所誤差制御(||y_next - y_old || <= rtol * ||y_old|| + atol)を用います。
局所誤差制御を行う補外法と固定ステップh (= (x - x0) / div_num)を使ってy(x)の近似値yを計算します。。
[テスト用関数] 局所誤差制御を行う補外法を用いてy(x)の近似値yを計算します。。ansfunc関数に真の解を設定すると誤差の計算も行えます。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
実係数多項式polyから，Aberthの初期値決定のための円の中心を求めます。多倍長型の場合はcenterに中心の座標が格納されます。実係数なので中心は必ず実軸上にのります。
実係数多項式polyから，Aberthの初期値を計算するための円の半径を求めます。多倍長型の場合はradに半径が格納されます。
実係数多項式polyから，Aberthの初期値を求め，配列init_arrayに格納します。
poly = 0という実係数代数方程式の近似解をDKA法で計算します。初期値init_arrayから出発し，最大反復回数maxtimes以下で全ての近似解の変化が停止条件式abs_eps + rel_eps answer_array[i]以下になれば，配列answer_arrayに近似解を格納します。返り値は反復回数です。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
使用時には必ずbncsparse.hをインクルードした上でlibbncsparse.aをリンクして下さい。
疎行列を初期化します。
精度指定付きで多倍長疎行列を初期化します。
疎行列matを消去します。
行方向に疎行列matの非零成分を探索し，nzero_row_dim[]に行ごとの非零成分数をセットします。
疎行列matを標準出力に表示します。
密行列org_matを疎行列形式に変換したものを初期化・格納して返します。
ファイル名fnameに格納された疎行列の形式を調べて返します。疎行列初期化に必要な変数が自動的に設定できます。
ファイル名fnameに格納されたリンク情報に基づいて疎行列(グラフ表現)を読み取ります。
疎行列matとベクトルvecとの積を計算し，retに返します。
疎行列matの転置とベクトルvecとの積を計算し，retに返します。
疎行列matの絶対値最大固有値と固有ベクトルをべき乗法を用いて求めます。
ベクトルのスカラー倍を計算します。
ベクトルvecの絶対値最大成分の番号と値を返します。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
| 4.1 基本関数 | ||
| 4.2 複素数 | ||
| 4.3 多項式 | ||
| 4.4 線型計算 | ||
| 4.5 LU分解 | ||
| 4.6 Conjugate-Gradient 法 | ||
| 4.7 Newton 法 | ||
| 4.8 Runge-Kutta法 | 常微分方程式の初期値問題 | |
| 4.9 DKA法 | Durand-Kerner-Aberth法 |
テストプログラムの内容は以下の通りです。
基本関数のテストプログラム
複素数のテストプログラム
多項式のテストプログラム
基本線型計算のテストプログラム
LU分解による連立一次方程式の求解
CG法による連立一次方程式の求解
Newton法による非線型方程式の求解
Runge-Kutta法による常微分方程式の初期値問題の求解
DKA法による代数方程式の求解
コンパイル方法は，前述の通り
%cc test_complex.c -lbnc |
とする。-lbncの代わりに直接/libdir/libbnc.aとライブラリファイルを指定してもよい。
%cc -DUSE_GMP test_complex.c -lbnc -lgmp |
とする。/libdir/libbnc.a /libdir/libgmp.aと直接ライブラリファイルを指定してもよい。
%cc -DUSE_GMP -DUSE_MPFR test_complex.c -lbnc -lmpfr -lgmp |
とする。/libdir/libbnc.a /libdir/libmpfr.a /libdir/libgmp.aと直接ライブラリファイルを指定してもよい。MPFRパッケージを利用する場合は，`mpfr.h'と`mpf2mpfr.h'の両方を使用するので，共に読み込むことの出来るディレクトリ位置に置いておくこと。
とします。
以下の節で，BNCpackの関数群の使い方を紹介します。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
IEEE754の初等関数は出揃っているため，C標準のもの(`math.h'定義のもの)を利用します。従って，ここでは多倍長計算用のサンプルのみ示すことにします。
#include <stdio.h> #include <math.h> #include "bnc.h" |
必要なヘッダファイルを読み込みます。bnc.hは多倍長計算を要求する定義(USE_GMP及びUSE_MPFR)が存在すれば，適宜必要となるヘッダファイル(`gmp.h', `mpfr.h')をその中で読み込んでいます。
main()
{
#ifndef USE_GMP
printf("This program are not available without GMP.\n");
#else
|
このように，多倍長計算のみで実行したい部分は，USE_GMPもしくはUSE_MPFRの定義が在るかどうかをチェックする#ifdefもしくは#ifndef～#endifで括って下さい。
long int i, max_times;
double x_min, x_max, x, h;
mpf_t mp_x_min, mp_x_max, mp_x, mp_h, \
mp_sin, mp_cos, mp_exp, mp_pi, mp_e, mp_ln2, mp_ln, mp_log10;
|
多倍長計算を行うための変数は全てmpf_t型として宣言して下さい。MPFRパッケージ利用の場合は，これがマクロで自動的に全てmpfr_t型に置き換えられます(9)。
set_bnc_default_prec(128); |
デフォルトの精度桁(多倍長浮動小数点数の仮数部桁(2進))を指定しています。この場合は128bit(10進38.5桁)となります。
max_times = 128;
x_min = -30.0;
x_max = 30.0;
mpf_init_set_si(mp_x_min, -30);
mpf_init_set_si(mp_x_max, 30);
h = (x_max - x_min) / max_times;
mpf_init(mp_h);
mpf_init(mp_sin);
mpf_init(mp_cos);
mpf_init(mp_exp);
mpf_init(mp_pi);
mpf_init(mp_e);
mpf_init(mp_ln2);
mpf_init(mp_ln);
mpf_init(mp_log10);
mpf_sub(mp_h, mp_x_max, mp_x_min);
mpf_div_ui(mp_h, mp_h, (unsigned long)max_times);
mpf_init_set(mp_x, mp_x_min);
x = x_min;
|
多倍長計算を使いこなすには当然，GMP及びMPFRで定義されている関数群をそのまま利用する必要が出てきます。詳細については各マニュアルを参照して下さい。
printf(" x sin(x) \
mpf_sin(x)\n");
for(i = 0; i < max_times; i++)
{
printf("%25.17e %25.17e ", x, sin(x));
mpf_sin(mp_sin, mp_x);
// mpf_out_str(stdout, 10, 0, mp_sin);
printf("%25.17e %25.17e", mpf2double(mp_sin), cos(x));
mpf_cos(mp_cos, mp_x);
// mpf_out_str(stdout, 10, 0, mp_cos);
printf("%25.17e %25.17e", mpf2double(mp_cos), exp(x));
mpf_exp(mp_exp, mp_x);
// mpf_out_str(stdout, 10, 0, mp_exp);
printf("%25.17e %25.17e", mpf2double(mp_exp), log(x));
mpf_ln(mp_ln, mp_x);
// mpf_out_str(stdout, 10, 0, mp_ln);
printf("%25.17e %25.17e", mpf2double(mp_ln), \
log(x)/log(10.0));
mpf_log10(mp_log10, mp_x);
// mpf_out_str(stdout, 10, 0, mp_log10);
printf("%25.17e\n", mpf2double(mp_log10));
x += h;
mpf_add(mp_x, mp_x, mp_h);
}
mpf_set_ui(mp_x, 100UL);
printf("sin(%25.17e): ", mpf2double(mp_x));
mpf_sin(mp_sin, mp_x);
mpf_out_str(stdout, 10, 0, mp_sin);
printf(" %25.17e", sin(mpf2double(mp_x)));
printf("\n");
printf("cos(%25.17e): ", mpf2double(mp_x));
mpf_cos(mp_cos, mp_x);
mpf_out_str(stdout, 10, 0, mp_cos);
printf(" %25.17e", cos(mpf2double(mp_x)));
printf("\n");
printf("exp(%25.17e): ", mpf2double(mp_x));
mpf_exp(mp_exp, mp_x);
mpf_out_str(stdout, 10, 0, mp_exp);
printf(" %25.17e", exp(mpf2double(mp_x)));
printf("\n");
printf("ln (%25.17e): ", mpf2double(mp_x));
mpf_ln(mp_ln, mp_x);
mpf_out_str(stdout, 10, 0, mp_ln);
printf(" %25.17e", log(mpf2double(mp_x)));
printf("\n");
printf("log10(%25.17e): ", mpf2double(mp_x));
mpf_log10(mp_log10, mp_x);
mpf_out_str(stdout, 10, 0, mp_log10);
printf(" %25.17e", log(mpf2double(mp_x))/log(10.0));
printf("\n");
printf("PI: ");
mpf_pi(mp_pi);
mpf_out_str(stdout, 10, 0, mp_pi);
printf(" -> %25.17e", mpf2double(mp_pi));
mpf_floor(mp_pi, mp_pi);
printf(" floor(PI), floor(PI*10000):");
mpf_out_str(stdout, 10, 0, mp_pi);
mpf_pi(mp_pi); mpf_mul_ui(mp_pi, mp_pi, 10000UL);
mpf_floor(mp_pi, mp_pi);
printf(", "); mpf_out_str(stdout, 10, 0, mp_pi);
fflush(stdout);
printf("\n");
printf("E : ");
mpf_e(mp_e);
mpf_out_str(stdout, 10, 0, mp_e);
printf(" -> %25.17e", mpf2double(mp_e));
mpf_floor(mp_e, mp_e);
printf(" floor(E):");
mpf_out_str(stdout, 10, 0, mp_e);
printf("\n");
printf("log 2 : ");
mpf_ln_2(mp_ln2);
mpf_out_str(stdout, 10, 0, mp_ln2);
printf(" -> %25.17e", mpf2double(mp_ln2));
mpf_floor(mp_ln2, mp_ln2);
printf(" floor(ln2):");
mpf_out_str(stdout, 10, 0, mp_ln2);
printf("\n");
|
ここで使われているGMPには定義されていない関数mpf_sin, mpf_cos, mpf_exp, mpf_ln, mpf_log10等の関数は，前述の通り，MPFRパッケージが読み込まれると全てそこで定義されているものに置き換えられます。BNCpackでGMP用に作成されたこれらの関数は非常に低速ですので，高速な初等関数が必要な時はMPFRパッケージとの併用をお勧めします。
mpf_clear(mp_h);
mpf_clear(mp_sin);
mpf_clear(mp_cos);
mpf_clear(mp_exp);
mpf_clear(mp_pi);
mpf_clear(mp_e);
mpf_clear(mp_ln2);
mpf_clear(mp_ln);
mpf_clear(mp_log10);
#endif
}
|
最後に，使用した変数を解放します。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
BNCpackで定義した複素数型は独自のもので，C++標準のcomplexクラスとも，GSL(http://sources.redhat.com/gsl/)のそれとも互換性はありません。
とある先生からはC++のcomplexクラスでmpf_tもしくはmpfr_t型が利用できるようにならないか？というご要望を受けています。テンプレートですから，ちょっと苦労すれば出来そうな気もしますが・・・どなたかやりません？
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
DCmplx dca, dcb, dcc; |
IEEE754倍精度複素数DCmplx型の変数を宣言します。単精度の場合はFCmplxと宣言します。
/* init */ dca = init_dcmplx(); dcb = init_dcmplx(); dcc = init_dcmplx(); |
複素数dca, dcb, dccの領域を確保して初期化(0にセット)します。
/* set */
set_real_dcmplx(dca, 1.0);
set_image_dcmplx(dca, 2.0);
set_real_dcmplx(dcb, 3.0);
set_image_dcmplx(dcb, 4.0);
set_real_dcmplx(dcc, (double)rand());
set_image_dcmplx(dcc, (double)rand());
/* print */
printf("dca = "); print_dcmplx(dca);
printf("dcb = "); print_dcmplx(dcb);
printf("dcc = "); print_dcmplx(dcc);
|
dca = 1 + 2i, dcb = 3 + 4i，dccの実数部，虚数部に乱数をセットし，表示します。
/* basic arithmetic */
add_dcmplx(dcc, dca, dcb);
printf("dca + dcb = "); print_dcmplx(dcc);
sub_dcmplx(dcc, dca, dcb);
printf("dca - dcb = "); print_dcmplx(dcc);
add2_dcmplx(dcc, dcb);
printf("dca - dcb + dcb = "); print_dcmplx(dcc);
mul_dcmplx(dcc, dca, dcb);
printf("dca * dcb = "); print_dcmplx(dcc);
div_dcmplx(dcc, dca, dcb);
printf("dca / dcb = "); print_dcmplx(dcc);
mul2_dcmplx(dcc, dcb);
printf("dca / dcb * dcb = "); print_dcmplx(dcc);
printf("~dca = ");
conj_dcmplx(dcc, dca); print_dcmplx(dcc);
printf("|dca| = %25.17e\n", abs_dcmplx(dca));
printf("exp(i*dca) = "); i
exp_dcmplx(dcc, abs_dcmplx(dca)); print_dcmplx(dcc);
|
四則演算，初等関数を実行し，その結果を表示します。
/* clear */ free_dcmplx(dca); free_dcmplx(dcb); free_dcmplx(dcc); |
最後に，確保した複素数型の変数領域を解放します。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
MPFCmplx mpfca, mpfcb, mpfcc; mpf_t tmp; |
多倍長複素数MPFCmplx型の変数を宣言します。
/* init */ set_bnc_default_prec(128); |
デフォルトの精度を128bitとします。
mpfca = init_mpfcmplx(); mpfcb = init_mpfcmplx(); mpfcc = init2_mpfcmplx(256); |
複素数mpfca, mpfcb, mpfccの領域を確保して初期化(0にセット)します。mpfccはデフォルトの精度とは異なり，256bitの精度を持つように初期化されています。
/* set */
set_real_mpfcmplx_ui(mpfca, 1);
set_image_mpfcmplx_ui(mpfca, 2);
set_real_mpfcmplx_ui(mpfcb, 3);
set_image_mpfcmplx_ui(mpfcb, 4);
set_real_mpfcmplx_ui(mpfcc, rand());
set_image_mpfcmplx_ui(mpfcc, rand());
/* print */
printf("mpfca = "); print_mpfcmplx(mpfca);
printf("mpfcb = "); print_mpfcmplx(mpfcb);
printf("mpfcc = "); print_mpfcmplx(mpfcc);
|
dca = 1 + 2i, dcb = 3 + 4i，dccの実数部，虚数部に乱数をセットし，表示します。
/* basic arithmetic */
add_mpfcmplx(mpfcc, mpfca, mpfcb);
printf("mpfca + mpfcb = "); print_mpfcmplx(mpfcc);
sub_mpfcmplx(mpfcc, mpfca, mpfcb);
printf("mpfca - mpfcb = "); print_mpfcmplx(mpfcc);
add2_mpfcmplx(mpfcc, mpfcb);
printf("mpfca - mpfcb + mpfcb = "); print_mpfcmplx(mpfcc);
mul_mpfcmplx(mpfcc, mpfca, mpfcb);
printf("mpfca * mpfcb = "); print_mpfcmplx(mpfcc);
div_mpfcmplx(mpfcc, mpfca, mpfcb);
printf("mpfca / mpfcb = "); print_mpfcmplx(mpfcc);
mul2_mpfcmplx(mpfcc, mpfcb);
printf("mpfca / mpfcb * mpfcb = "); print_mpfcmplx(mpfcc);
printf("~mpfca = ");
conj_mpfcmplx(mpfcc, mpfca); print_mpfcmplx(mpfcc);
mpf_init(tmp); abs_mpfcmplx(tmp, mpfca);
printf("|mpfca| = ");
mpf_out_str(stdout, 10, 0, tmp); printf("\n");
printf("exp(i*mpfca) = ");
iexp_mpfcmplx(mpfcc, tmp); print_mpfcmplx(mpfcc);
|
四則演算，初等関数を実行し，その結果を表示します。
mpf_clear(tmp); /* clear */ free_mpfcmplx(mpfca); free_mpfcmplx(mpfcb); free_mpfcmplx(mpfcc); |
最後に，確保した複素数型の変数領域を解放します。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
#define MAX_POLY_LEN 4096
#define MAX_DEGREE 1024
main()
{
long int i;
DPoly dpa, dpb, dpc;
DCmplx dca, dcret;
|
IEEE754倍精度実係数を持つ多項式型(DPoly)の変数dpa, dpb, dpcと倍精度複素数型(DCmplx)の変数dca, dcretを宣言します。以下，複素数型の説明は省略します。
/* init */ dpa = init_dpoly(MAX_POLY_LEN); dpb = init_dpoly(MAX_POLY_LEN); dpc = init_dpoly(MAX_POLY_LEN); dca = init_dcmplx(); dcret = init_dcmplx(); set_real_dcmplx(dca, 1.0); set_image_dcmplx(dca, 1.0); |
多項式型変数を初期化します。
for(i = 0; i <= MAX_DEGREE; i++)
{
set_dpoly_i(dpa, i, (double)i);
set_dpoly_i(dpb, i, (double)rand());
set_dpoly_i(dpc, i, (double)rand());
|
多項式の係数をセットします。
printf("dpa: O(x^%d)\n", setdegree_dpoly(dpa));print_dpoly(dpa);
printf("dpa (1) = %25.17e\n", eval_dpoly(dpa, 1.0));
printf("dpa'(1) = %25.17e\n", eval_diff_dpoly(dpa, 1.0));
printf("dpa (1+1i) = "); ceval_dpoly(dcret, dpa, dca);
print_dcmplx(dcret);
printf("dpa'(1+i1) = "); ceval_diff_dpoly(dcret, dpa, dca);
print_dcmplx(dcret);
printf("dpa (2) = %25.17e\n", eval_dpoly(dpa, 2.0));
printf("dpa'(2) = %25.17e\n", eval_diff_dpoly(dpa, 2.0));
printf("dpb: \n");print_dpoly(dpb);
printf("dpc: \n");print_dpoly(dpc);
|
多項式及びその導関数を使った計算を行い，結果を出力します。
/* clear */ free_dpoly(dpa); free_dpoly(dpb); free_dpoly(dpc); } |
最後に多項式型変数領域を解放します。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
#define MAX_POLY_LEN 4096
#define MAX_DEGREE 1024
main()
{
long int i;
MPFPoly mpf_pa, mpf_pb, mpf_pc;
mpf_t mpf_x, mpf_ret;
MPFCmplx mpfca, mpf_cret;
|
多倍長実係数を持つ多項式型(MPFPoly)の変数mpf_pa, mpf_pb, mpf_pcと多倍長複素数型(MPFCmplx)の変数mpfca, mpc_cretを宣言します。
/* init */ mpf_pa = init2_mpfpoly(MAX_POLY_LEN, 128); mpf_pb = init2_mpfpoly(MAX_POLY_LEN, 256); mpf_pc = init2_mpfpoly(MAX_POLY_LEN, 1024); mpfca = init2_mpfcmplx(1024); mpf_cret = init2_mpfcmplx(1024); set_real_mpfcmplx_ui(mpfca, 1UL); set_image_mpfcmplx_ui(mpfca, 1UL); |
多項式型変数を初期化します。
for(i = 0; i <= MAX_DEGREE; i++)
{
set_mpfpoly_i_d(mpf_pa, i, (double)i);
set_mpfpoly_i_d(mpf_pb, i, (double)rand());
set_mpfpoly_i_d(mpf_pc, i, (double)rand());
|
多項式の係数をセットします。IEEE754倍精度の実数をセットすると2進数のビットパターンがそのまま転写されますので，精度は倍精度程度で止まってしまうことに留意して下さい(10)。
mpf_init2(mpf_x, 128);
mpf_init2(mpf_ret, 128);
printf("mpf_pa: O(x^%d)\n", setdegree_mpfpoly(mpf_pa));
print_mpfpoly(mpf_pa);
mpf_set_ui(mpf_x, 1UL);
printf("mpf_pa(1) = ");
eval_mpfpoly(mpf_ret, mpf_pa, mpf_x);
mpf_out_str(stdout, 0, 10, mpf_ret);
printf("\n");
printf("mpf_pa'(1) = ");
eval_diff_mpfpoly(mpf_ret, mpf_pa, mpf_x);
mpf_out_str(stdout, 0, 10, mpf_ret);
printf("\n");
mpf_set_ui(mpf_x, 1UL);
printf("mpf_pa(1+1i) = ");
ceval_mpfpoly(mpf_cret, mpf_pa, mpfca);
print_mpfcmplx(mpf_cret);
printf("mpf_pa'(1+1i) = ");
ceval_diff_mpfpoly(mpf_cret, mpf_pa, mpfca);
print_mpfcmplx(mpf_cret);
mpf_set_ui(mpf_x, 2UL);
printf("mpf_pa(2) = ");
eval_mpfpoly(mpf_ret, mpf_pa, mpf_x);
mpf_out_str(stdout, 0, 10, mpf_ret);
printf("\n");
printf("mpf_pa'(2) = ");
eval_diff_mpfpoly(mpf_ret, mpf_pa, mpf_x);
mpf_out_str(stdout, 0, 10, mpf_ret);
printf("\n");
mpf_set_ui(mpf_x, 100000UL);
printf("mpf_pa(100000) = ");
eval_mpfpoly(mpf_ret, mpf_pa, mpf_x);
mpf_out_str(stdout, 0, 10, mpf_ret);
printf("\n");
printf("mpf_pa'(100000) = ");
eval_diff_mpfpoly(mpf_ret, mpf_pa, mpf_x);
mpf_out_str(stdout, 0, 10, mpf_ret);
printf("\n");
printf("mpf_pb: \n");print_mpfpoly(mpf_pb);
printf("mpf_pc: \n");print_mpfpoly(mpf_pc);
|
多項式及びその導関数を使った計算を行い，結果を出力します。
mpf_clear(mpf_x); mpf_clear(mpf_ret); /* clear */ free_mpfpoly(mpf_pa); free_mpfpoly(mpf_pb); free_mpfpoly(mpf_pc); } |
最後に，使用した変数領域を解放します。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
線型計算のサンプルは次のLU分解，CG法の例で代替します。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
DMatrix da; DVector db, dx, dans; long int ret_f, ret_d, ret_mpf; long int row_ch[DIM], col_ch[DIM]; long int i, j; /* initialize */ da = init_dmatrix(DIM, DIM); db = init_dvector(DIM); dx = init_dvector(DIM); dans = init_dvector(DIM); |
IEEE754倍精度実行列型DMatrix，同ベクトル型DVectorの変数を宣言し，初期化します。
/* get problem */ get_dproblem(da, db, dans); print_dmatrix(da); |
テスト問題の係数行列，定数ベクトル，真の解を与え，係数行列を表示します。テスト問題を作る関数get_dproblemは，例えば以下のように記述します。
void get_dproblem(DMatrix a, DVector b, DVector ans)
{
long int i, j, k;
double tmp;
/* Lotkin Matrix */
for(i = 0; i < a->col_dim; i++)
set_dmatrix_ij(a, 0, i, 1.0);
for(i = 1; i < a->row_dim; i++)
{
for(j = 0; j < a->col_dim; j++)
set_dmatrix_ij(a, i, j, 1.0 / (i + j + 1));
}
/* Answer */
for(i = 0; i < ans->dim; i++)
set_dvector_i(ans, i, (double)i);
/* Make constant vector */
mul_dmatrix_dvec(b, a, ans);
}
|
/* run DLUdecomp & SolveDLS */ // ret_d = DLUdecomp(da); // ret_d = DLUdecompP(da, row_ch); ret_d = DLUdecompC(da, row_ch, col_ch); |
係数行列のLU分解を計算します。DLUdecompはpivotingなし，DLUdcompPは部分pivotingのみ実行，DLUdecompCは完全pivotingを実行します。
// ret_d = SolveDLS(dx, da, db);
// ret_d = SolveDLSP(dx, da, db, row_ch);
ret_d = SolveDLSC(dx, da, db, row_ch, col_ch);
/* print */
printf(" i row_ch[i] col_ch[i]\n");
for(i = 0; i < DIM; i++)
printf("%5ld %25.17e %25.17e\n", i, get_dvector_i(dx, i),\
get_dvector_i(dans, i));
|
LU分解された係数行列を用いて後退代入を行い，解を出力します。
/* end */ free_dmatrix(da); free_dvector(db); free_dvector(dx); free_dvector(dans); |
最後に，使用した変数領域を解放します。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
MPFMatrix mpfa; MPFVector mpfb, mpfx, mpfans; mpf_t reps, aeps; long int ret_f, ret_d, ret_mpf; long int row_ch[DIM], col_ch[DIM]; long int i, j; set_bnc_default_prec(256); /* initialize */ mpf_init(reps); mpf_init(aeps); mpfa = init_mpfmatrix(DIM, DIM); // mpfa = init2_mpfmatrix(DIM, DIM, 256); mpfb = init_mpfvector(DIM); // mpfb = init2_mpfvector(DIM, 256); mpfx = init_mpfvector(DIM); // mpfx = init2_mpfvector(DIM, 256); mpfans = init_mpfvector(DIM); |
多倍長実行列型MPFMatrix，同ベクトル型MPFVectorの変数を宣言し，デフォルトの精度を256bitにして初期化します。
/* get problem */ get_mpfproblem(mpfa, mpfb, mpfans); print_mpfmatrix(mpfa); |
テスト問題の係数行列，定数ベクトル，真の解を与え，係数行列を表示します。テスト問題を作る関数get_mpfproblemは，例えば以下のように記述します。
void get_mpfproblem(MPFMatrix a, MPFVector b, MPFVector ans)
{
long int i, j, k;
mpf_t tmp;
mpf_init(tmp);
/* Lotkin Matrix */
for(i = 0; i < a->col_dim; i++)
set_mpfmatrix_ij_d(a, 0, i, 1.0);
for(i = 1; i < a->row_dim; i++)
{
for(j = 0; j < a->col_dim; j++)
{
mpf_set_ui(tmp, 1UL);
mpf_div_ui(tmp, tmp, \
(unsigned long)(i + j + 1));
set_mpfmatrix_ij(a, i, j, tmp);
}
}
/* Answer */
for(i = 0; i < ans->dim; i++)
{
mpf_set_si(tmp, i);
set_mpfvector_i(ans, i, tmp);
}
/* Make constant vector */
mul_mpfmatrix_mpfvec(b, a, ans);
}
|
/* run MPFLUdecomp & SolveMPFLS */ ret_mpf = MPFLUdecomp(mpfa); // ret_mpf = MPFLUdecompP(mpfa, row_ch); // ret_mpf = MPFLUdecompC(mpfa, row_ch, col_ch); |
係数行列のLU分解を計算します。MPFLUdecompはpivotingなし，MPFLUdcompPは部分pivotingのみ実行，MPFLUdecompCは完全pivotingを実行します。
ret_mpf = SolveMPFLS(mpfx, mpfa, mpfb);
// ret_mpf = SolveMPFLSP(mpfx, mpfa, mpfb, row_ch);
// ret_mpf = SolveMPFLSC(mpfx, mpfa, mpfb, row_ch, col_ch);
/* print */
for(i = 0; i < DIM; i++)
{
printf("%5ld ", i);
mpf_out_str(stdout, 10, 0, get_mpfvector_i(mpfx, i));
printf(" ");
mpf_out_str(stdout, 10, 0, get_mpfvector_i(mpfans, i));
printf("\n");
|
LU分解された係数行列を用いて後退代入を行い，解を出力します。
/* end */ mpf_clear(reps); mpf_clear(aeps); free_mpfmatrix(mpfa); free_mpfvector(mpfb); free_mpfvector(mpfx); free_mpfvector(mpfans); |
最後に，使用した変数領域を解放します。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
/* initialize */ da = init_dmatrix(DIM, DIM); db = init_dvector(DIM); dx = init_dvector(DIM); dans = init_dvector(DIM); /* get problem */ get_dproblem(da, db, dans); print_dmatrix(da); |
連立一次方程式の係数行列da，定数ベクトルdb, 近似解を格納するベクトルdx，真の解ベクトルdansを宣言し，値を代入して，係数行列を表示します。
/* run DCG */
itimes_d = DCG(dx, da, db, 1.0e-13, 1.0e-99, DIM * 5);
/* print */
for(i = 0; i < DIM; i++)
printf("%5ld %25.17e %25.17e\n",
i,
get_dvector_i(dx, i),
get_dvector_i(dans, i));
|
CG法(Conjugate-Gradient法)を実行し，数値解を表示します。
/* end */ free_dmatrix(da); free_dvector(db); free_dvector(dx); free_dvector(dans); |
使用した変数領域を解放します。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
set_bnc_default_prec(128); /* initialize */ mpf_init(reps); mpf_init2(reps2, 256); mpf_init2(reps3, 512); mpf_init(aeps); mpf_init2(aeps2, 256); mpf_init2(aeps3, 512); mpfa = init_mpfmatrix(DIM, DIM); mpfa2 = init2_mpfmatrix(DIM, DIM, 256); mpfa3 = init2_mpfmatrix(DIM, DIM, 512); mpfb = init_mpfvector(DIM); mpfb2 = init2_mpfvector(DIM, 256); mpfb3 = init2_mpfvector(DIM, 512); mpfx = init_mpfvector(DIM); mpfx2 = init2_mpfvector(DIM, 256); mpfx3 = init2_mpfvector(DIM, 512); mpfans = init_mpfvector(DIM); mpfans2 = init2_mpfvector(DIM, 256); mpfans3 = init2_mpfvector(DIM, 512); /* get problem */ get_mpfproblem(mpfa, mpfb, mpfans); get_mpfproblem(mpfa2, mpfb2, mpfans2); get_mpfproblem(mpfa3, mpfb3, mpfans3); print_mpfmatrix(mpfa); print_mpfmatrix(mpfa2); print_mpfmatrix(mpfa3); /* run MPFFCG */ mpf_set_d(reps, 1.0e-20); mpf_set_d(reps2, 1.0e-20); mpf_set_d(reps3, 1.0e-20); mpf_set_d(aeps, 1.0e-50); mpf_set_d(aeps2, 1.0e-50); mpf_set_d(aeps3, 1.0e-50); |
連立一次方程式の係数行列mpfa/mpfa2/mpfa3，定数ベクトルmpfb/mpfb2/mpfb3, 近似解を格納するベクトルmpfx/mpfx2/mpfx3，真の解ベクトルmpfans/mpfans2/mpfans3を，それぞれ128bit, 256bit, 512bitの精度を持つように宣言し，値を代入して，係数行列をそれぞれ表示します。
itimes_mpf = MPFCG(mpfx, mpfa, mpfb, reps, aeps, DIM * 5);
itimes_mpf2 = MPFCG(mpfx2, mpfa2, mpfb2, reps2, aeps2, DIM * 5);
itimes_mpf3 = MPFCG(mpfx3, mpfa3, mpfb3, reps3, aeps3, DIM * 5);
/* print */
for(i = 0; i < DIM; i++)
{
printf("%5ld ", i);
mpf_out_str(stdout, 10, 0, get_mpfvector_i(mpfx, i));
printf(" ");
mpf_out_str(stdout, 10, 0, get_mpfvector_i(mpfx2, i));
printf(" ");
mpf_out_str(stdout, 10, 0, get_mpfvector_i(mpfx3, i));
printf(" ");
mpf_out_str(stdout, 10, 0, get_mpfvector_i(mpfans, i));
printf("\n");
|
CG法(Conjugate-Gradient法)を実行し，各精度の数値解を表示します。
/* end */ mpf_clear(reps); mpf_clear(aeps); mpf_clear(reps2); mpf_clear(aeps2); mpf_clear(reps3); mpf_clear(aeps3); free_mpfmatrix(mpfa); free_mpfmatrix(mpfa2); free_mpfmatrix(mpfa3); free_mpfvector(mpfb); free_mpfvector(mpfb2); free_mpfvector(mpfb3); free_mpfvector(mpfx); free_mpfvector(mpfx2); free_mpfvector(mpfx3); free_mpfvector(mpfans); free_mpfvector(mpfans2); free_mpfvector(mpfans3); |
最後に使用した変数領域を解放します。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
double fd(double x)
{
/* x^2 - 2 = 0 */
return x * x - 2.0;
}
double dfd(double x)
{
/* 2*x */
return 2.0 * x;
}
|
Newton法及び簡易Newton法で使用する方程式f(x) = 0の左辺の関数f，及びその導関数dfの定義をしています。
/* double */
times = dnewton_1(&dans, 1.0, fd, dfd, 100, 1.0e-50, 1.0e-10);
printf(" times : %ld\n", times);
printf("dnewton_1: %25.17e\n", dans);
times = dsnewton_1(&dans, 1.0, fd, dfd, 100, 1.0e-50, 1.0e-10);
printf(" times : %ld\n", times);
printf("dsnewton_1: %25.17e\n", dans);
|
一次元のNewton法及び簡易Newton法を実行してその結果を出力します。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
void f(mpf_t ret, mpf_t x)
{
/* x^2 - 2 = 0 */
mpf_mul(ret, x, x);
mpf_sub_ui(ret, ret, 2UL);
}
void df(mpf_t ret, mpf_t x)
{
/* 2*x */
mpf_mul_ui(ret, x, 2UL);
}
|
Newton法及び簡易Newton法で使用する方程式f(x) = 0の左辺の関数f，及びその導関数dfの定義をしています。
/* mpf_t */
set_bnc_default_prec(1024);
mpf_init(ans);
mpf_init(x0);
mpf_init_set_d(aeps, 1.0e-100);
mpf_init_set_d(reps, 1.0e-50);
/* newton_1 */
mpf_set_ui(x0, 1UL);
times = mpf_newton_1(ans, x0, f, df, 100, aeps, reps);
printf(" times : %ld\n", times);
printf("mpf_newton_1: "); mpf_out_str(stdout, 10, 0, ans);
printf("\n");
/* snewton_1 */
times = mpf_snewton_1(ans, x0, f, df, 100, aeps, reps);
printf(" times : %ld\n", times);
printf("mpf_snewton_1: "); mpf_out_str(stdout, 10, 0, ans);
printf("\n");
mpf_sqrt_ui(ans, 2UL);
printf("mpf_sqrt: "); mpf_out_str(stdout, 10, 0, ans); printf("\n");
|
一次元のNewton法及び簡易Newton法を実行してその結果を出力します。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
使用時には必ずbncode.hをインクルードした上でlibbncode.aをリンクして下さい。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
常微分方程式の定義(ddf)と初期値の設定を行います。
void ddf(DVector y, double x0, DVector y0)
{
/* y' = y */
subst_dvector(y, y0);
return;
void initial_dvalue(double *lf_initx, double *lf_maxx, DVector lf_inity)
{
long int i;
/* x0 := 0 */
/* y := [1, ..., 1]^T */
*lf_initx = 0.0;
*lf_maxx = 1.0;
for(i = 0; i < lf_inity->dim; i++)
set_dvector_i(lf_inity, i, 1.0);
|
main関数では次のように関数を呼び出します。
DVector dy, init_dy, lf_dtmp; long int div_num; double dx, init_dx, dstepsize, dabs_tol, drel_tol; dy = init_dvector(DIM); lf_dtmp = init_dvector(DIM); init_dy = init_dvector(DIM); initial_dvalue(&init_dx, &dx, init_dy); dstepsize = 0.5; dabs_tol = 0.0; drel_tol = 1.0e-15; /* 固定ステップ */ div_num = 128; // 128分割 derk_fs(NULL, dx, dy, init_dx, init_dy, div_num, ddf, ERK_BUTCHER1_7_6); /* 局所誤差制御付き */ derk_ex(NULL, dx, dy, init_dx, init_dy, dstepsize, ddf, ERK_BUTCHER1_7_6, drel_tol, dabs_tol); |
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
常微分方程式の定義(ddf)と初期値の設定を行います。
void df(MPFVector y, mpf_t x0, MPFVector y0)
{
/* y' = y */
subst_mpfvector(y, y0);
return;
void initial_value(mpf_t lf_initx, mpf_t lf_maxx, MPFVector lf_inity)
{
long int i;
/* x0 := 0 */
/* y := [1, ..., 1]^T */
mpf_init_set_str(lf_initx, "0.0", 10);
mpf_init_set_str(lf_maxx, "1.0", 10);
for(i = 0; i < lf_inity->dim; i++)
set_mpfvector_i_str(lf_inity, i, "1.0", 10);
|
main関数では次のように関数を呼び出します。
MPFVector y, init_y, lf_tmp; mpf_t x, init_x, stepsize, abs_tol, rel_tol; set_bnc_default_prec(256); /* x := 1 */ mpf_init(init_x); mpf_init(x); y = init_mpfvector(DIM); lf_tmp = init_mpfvector(DIM); init_y = init_mpfvector(DIM); initial_value(init_x, x, init_y); mpf_init_set_d(stepsize, dstepsize); mpf_init_set_d(abs_tol, dabs_tol); mpf_init_set_d(rel_tol, drel_tol); /* 固定ステップ */ div_num = 128; // 128分割 mpf_erk_fs(NULL, x, y, init_x, init_y, div_num, df, ERK_BUTCHER1_7_6); /* 局所誤差制御付き */ mpf_erk_ex(NULL, x, y, init_x, init_y, stepsize, df, ERK_BUTCHER1_7_6, el_tol, abs_tol); |
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
/* init */ dabs_eps = 1.0e-100; drel_eps = 1.0e-14; df = init_dpoly(MAX_LENGTH); cdans = init_cdarray(10); cdinit = init_cdarray(10); |
ここでDKA法の停止則に必要な値をセットします。
/* ff = (x-1)(x-2)...(x-10) */ set_dpoly_i(df, 0, (double)3628800); set_dpoly_i(df, 1, (double)-10628640); set_dpoly_i(df, 2, (double)12753576); set_dpoly_i(df, 3, (double)-8409500); set_dpoly_i(df, 4, (double)3416930); set_dpoly_i(df, 5, (double)-902055); set_dpoly_i(df, 6, (double)157773); set_dpoly_i(df, 7, (double)-18150); set_dpoly_i(df, 8, (double)1320); set_dpoly_i(df, 9, (double)-55); set_dpoly_i(df,10, (double)1); print_dpoly(df); |
多項式の係数を低次項からセットし，表示します。
/* set Aberth's initial value */ ddka_init(cdinit, df); print_cdarray(cdinit); |
Aberthの初期値を計算し，cdinitにセットします。
/* DKA method */ dtimes = ddka(cdans, cdinit, df, 1000, dabs_eps, drel_eps); |
DKA法を実行します。dtimesに反復回数が格納されます。
/* print answer */
printf("Iterative times: %d\n", dtimes);
print_cdarray(cdans);
|
反復回数と近似解(cdans)が表示されます。
/* clear */ free_dpoly(df); free_cdarray(cdans); free_cdarray(cdinit); |
最後に，使用した変数を解放します。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
set_bnc_default_prec(512); /* init */ mpf_init_set_d(mpfabs_eps, 1.0e-300); mpf_init_set_d(mpfrel_eps, 1.0e-25); |
ここで必要となる精度(512bit)と，DKA法の停止則に必要な値をセットします。
mpff = init_mpfpoly(MAX_LENGTH); cmpfans = init_cmpfarray(20); cmpfinit = init_cmpfarray(20); |
多倍長実数係数の多項式(mpff)，多倍長複素数の配列(cmpfans, cmpfinit)を初期化します。
/* ff = (x-1)(x-2)...(x-20) */ set_mpfpoly_i_str(mpff, 0, "2432902008176640000", 10); set_mpfpoly_i_str(mpff, 1, "-8752948036761600000", 10); set_mpfpoly_i_str(mpff, 2, "13803759753640704000", 10); set_mpfpoly_i_str(mpff, 3, "-12870931245150988800", 10); set_mpfpoly_i_str(mpff, 4, "8037811822645051776", 10); set_mpfpoly_i_str(mpff, 5, "-3599979517947607200", 10); set_mpfpoly_i_str(mpff, 6, "1206647803780373360", 10); set_mpfpoly_i_str(mpff, 7, "-311333643161390640", 10); set_mpfpoly_i_str(mpff, 8, "63030812099294896", 10); set_mpfpoly_i_str(mpff, 9, "-10142299865511450", 10); set_mpfpoly_i_str(mpff,10, "1307535010540395", 10); set_mpfpoly_i_str(mpff,11, "-135585182899530", 10); set_mpfpoly_i_str(mpff,12, "11310276995381", 10); set_mpfpoly_i_str(mpff,13, "-756111184500", 10); set_mpfpoly_i_str(mpff,14, "40171771630", 10); set_mpfpoly_i_str(mpff,15, "-1672280820", 10); set_mpfpoly_i_str(mpff,16, "53327946", 10); set_mpfpoly_i_str(mpff,17, "-1256850", 10); set_mpfpoly_i_str(mpff,18, "20615", 10); set_mpfpoly_i_str(mpff,19, "-210", 10); set_mpfpoly_i_str(mpff,20, "1", 10); print_mpfpoly(mpff); |
多項式の係数を低次項からセットし，表示します。
/* set Aberth's initial value */ mpf_dka_init(cmpfinit, mpff); print_cmpfarray(cmpfinit); |
Aberthの初期値を計算し，cmpfinitにセットします。
/* DKA method */ mpftimes = mpf_dka(cmpfans, cmpfinit, mpff, 1000, \ mpfabs_eps, mpfrel_eps); |
DKA法を実行します。mpftimesに反復回数が格納されます。
/* print answer */
printf("Iterative times: %d\n", mpftimes);
print_cmpfarray(cmpfans);
|
反復回数と近似解(cmpfans)が表示されます。
/* clear */ free_mpfpoly(mpff); free_cmpfarray(cmpfans); free_cmpfarray(cmpfinit); |
最後に，使用した変数を解放します。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
今後の予定としては
といったものがありますが，気が向かないと何もしない我が儘な性格なので，あまりあてにしないで下さい。
ご要望を頂いてもすぐにお答えできるかどうか分かりませんが，もし何かご意見などありましたら，マニュアル表紙のメールアドレスまでお願い致します。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
GNU MPなかりせば本ライブラリは存在しませんでした。開発チーム一同に感謝いたします。
また，プリンタがトラブって以来，一貫してFSFを率いつつ，「自由なソフトウェア」(Free Software)の普及活動に邁進し続けているRMSにも感謝いたします。GNU Projectなかりせば，開発環境に使っているLinuxも存在しなかったでしょうから。
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
| Jump to: | A C D E F G I L M N P S T |
|---|
| Jump to: | A C D E F G I L M N P S T |
|---|
| [Top] | [Contents] | [Index] | [ ? ] |
2011年8月現在の最新版はVersion 5.0.2です。URLはhttp://gmplib.org/
2011年8月現在の最新版はVersion 3.0.1です。URLはhttp://www.mpfr.org/
Microsoftの商標です。
今のところ，真面目に取り組むつもりはありません。気が向けば取りかかるかも。・・・でもそーゆー用途にはGSL(GNU Scientific Library, http://sources.redhat.com/gsl/)とかありますのでそちらをお使いになった方がいいでしょう。
今となってはそれほどでもないでしょうけど，自分としては売りのつもり。
MPFRパッケージの マニュアル参照。
実は共用できるようにしてありますが，表では使わないようになっています。
mpf_t型の丸め処理は切り捨てで行われます。
詳細はMPFRパッケージに同封されている`mpf2mpfr.h'を参照のこと。
それで何度失敗したことか・・・。
| [Top] | [Contents] | [Index] | [ ? ] |
| [Top] | [Contents] | [Index] | [ ? ] |
This document was generated by Tomonori Kouya on September, 1 2011 using texi2html 1.76.
The buttons in the navigation panels have the following meaning:
| Button | Name | Go to | From 1.2.3 go to |
|---|---|---|---|
| [ < ] | Back | previous section in reading order | 1.2.2 |
| [ > ] | Forward | next section in reading order | 1.2.4 |
| [ << ] | FastBack | beginning of this chapter or previous chapter | 1 |
| [ Up ] | Up | up section | 1.2 |
| [ >> ] | FastForward | next chapter | 2 |
| [Top] | Top | cover (top) of document | |
| [Contents] | Contents | table of contents | |
| [Index] | Index | index | |
| [ ? ] | About | about (help) |
where the Example assumes that the current position is at Subsubsection One-Two-Three of a document of the following structure:
This document was generated by Tomonori Kouya on September, 1 2011 using texi2html 1.76.