BNCpackはIEEE単精度・倍精度浮動小数点数(float, double)，及びMPFR/GNU MP(GMP)1の多倍長浮動小数点数(mpfr_tまたはmpf_t)をサポートする基本数値計算ライブラリです。全てANSI C(C++にあらず)で記述されています。
元々は，自分が使って来たCプログラムが貯まってきたので整理しよう，という動機で始めたものです。が，近年のPCの著しい機能向上とコスト低下のおかげで多倍長計算もお手軽な環境下で使い物になってきた，じゃあ，そいつも組み込んでしまえということで，やっつけ仕事で作ってしまった代物です。それがいつしか並列計算可能なものに進化してしまいました。
今回公開するものは，GMPをベースとした多倍長計算サポートというのが一番の目玉2です。GMPが主としてUNIX環境をターゲットにしている関係上，本ライブラリもUNIX環境で開発してきました。最後に，その環境を記しておきます。
具体的なBNCpackの利用方法については機能一覧，もしくはサンプルプログラムの章を参照して下さい。
最初に述べたように，BNCpackはGMPを利用した多倍長計算が売りのライブラリですが，IEEE754 standard(単精度，倍精度)で動作する関数も用意されています。頭文字を除いて同じ名前の関数が3つ並んでいたら，全く同じアルゴリズムを単精度(f/F)，倍精度(d/D)，多倍長精度(mpf/MPF)で実行するものだと思って間違いないでしょう。但し，今後は倍精度と多倍長精度を中心に拡充を行っていきます。
従来多倍長計算はGMPのmpf_t型を利用してきましたが，Version 0.3からは，表面上はmpft_t型宣言のまま，MPFRパッケージも利用することが出来るように改良されました。従って，Version 0.3以降のBNCpackは
以下，この3つの環境でBNCpackをインストールする方法を簡単に示します。
bncmake.incをご自分の環境に合わせて設定します。INCLUDESとLIBSの設定だけしておけば大丈夫でしょう。
INCLUDES = -I/usr/local/include -I/usr/include -I../ -I./
LIBS = -L/usr/local/lib:/usr/lib:./ $(LIBBNC) -lm
この後，
% make
として，src/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
としてコンパイル・リンクして使用して下さい。
GMPのインストール方法等についてはhttp://swox.com/gmp/をご参照下さい。特段変わった環境でなければ，Version 3.0.1以降については殆どのUNIX環境下で動作するようです。私の使用しているRPM系のLinuxでは
% ./configure; make; su
# make install
でコンパイルとインストールを行うことが出来ました。また，gfortranは内部でMPFRを利用しているので，GCC-4.x以降のコンパイラがインストールされている環境にはGMPもMPFRも既に入っていることが多いと思われます。
GMPがインストールしてあれば，bnc-x.x.x.tar.gzをダウンロード後，bncmake.incの環境変数等を適切に設定します。GMPPREFIXでGMPがインストールされているディレクトリを指定しています。
GMPPREFIX=/usr/local
INCLUDES = -DUSE_GMP -I$(GMPPREFIX)/include -I/usr/include -I../ -I./
LIBS = -L$(GMPPREFIX)/lib:/usr/lib:./ $(LIBBNC) -lgmp -lm
後は
% make
で，src/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をインクルードする前に定義しておく必要があります。
MPFRパッケージ(http://www.mpfr.org/)はGMPをベースに，IEEE754浮動小数点数と互換性を持つように改良された多倍長浮動小数点演算ライブラリです。GMPのmpf_t型との違いは
MPFRをコンパイルする際にはGMPが必要になります。GMPが/usr/localにインストールされている場合は
% ./configure --with-gmp=/usr/local; make; su
# make install
でコンパイルとインストールを行うことができます。なお，GMP 4.1.xに同梱されているMPFRは古いバージョンですので使わない方がいいでしょう。
あとはmpf_t型と同様，bnc-x.x.x.tar.gzをダウンロード後，Makefileのコンパイルオプション等を適切に設定します。この場合は，GMPPREFIXとMPFRPREFIXでGMPとMPFRのインストールディレクトリを指定する必要があります。
GMPPREFIX=/usr/local
MPFRPREFIX=/usr/local
DEFS = -DUSE_GMP -DUSE_MPFR
INCLUDES = $(DEFS) -I$(GMPPREFIX)/include -I$(MPFRPREFIX)/include \
-I$(MPIPREFIX) -I/usr/local/include -I/usr/include -I../ -I./
LIBSDIR = -L$(GMPPREFIX)/lib -L$(MPFRPREFIX)/lib -L$(MPIPREFIX)/lib \
-L/usr/local/lib -L/usr/lib -L./
LIBS = $(LIBSDIR) $(LIBBNC) -lmpfr -lgmp -lm
その後，
% make
で，src/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をインクルードする前に定義しておく必要があります。
構造体には，その構造体ごとの精度を保持する変数unsigned long precを含みます。それ以外の構成はIEEE754浮動小数点数を利用したものと変わりません。
上から順に，IEEE754単精度浮動小数点数による複素数型，倍精度浮動小数点数による複素数型，多倍長浮動小数点数による複素数型です。
この複素数型は次のように定義されています。
typedef struct{
float re; /* Real part */
float im; /* Imaginary part */
} fcmplx;
これらの複素数型を利用した関数群としては，DKA法が存在します。線型計算ではまだ利用できません。
上から順に，単精度浮動小数点数による多項式型，倍精度係数の多項式型，多倍長係数の多項式型です。
多項式型は次のように定義されています。
typedef struct{
float *coef; /* Coefficients of polynomial */
long deg; /* Degree of polynomial */
} fpoly;
typedef fpoly *FPoly;
これら多項式型を利用した関数群は，DKA法に利用されています。
上から順に，単精度浮動小数点数によるベクトル型，倍精度ベクトル型，多倍長ベクトル型です。
上から順に，単精度浮動小数点数の複素ベクトル型，倍精度複素ベクトル型，多倍長複素ベクトル型です。
ベクトル型は次のように定義されています。原則としては構造体へのポインタである[F,D,MPF,CF,CD,CMPF]Vectorだけを使います。
typedef struct{
float *element; /* Elements of vector */
long dim; /* Dimension of vector */
} fvector;
typedef fvector *FVector;
複素ベクトル型を使った関数群はまだ提供されていません。
上から順に，単精度複素一般行列型，倍精度複素一般行列型，多倍長複素一般行列型です。
一般行列型は次のように定義されています。原則としては構造体へのポインタである[F,D,MPF,CF,CD,CMPF]Matrixだけを使います。
typedef struct {
float *element; /* Elements */
long row_dim, col_dim; /* Rows, Columns */
} fmatrix;
typedef fmatrix *FMatrix;
複素一般行列型を扱う関数群はまだ提供されていません。
スタックは次のように定義されています。原則としては構造体へのポインタである[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;
スタックを利用した関数群はまだ提供されていません。
上から順に，単精度，倍精度，多倍長，複素単精度，複素倍精度，複素多倍長配列です。
配列は次のように定義されています。原則としては構造体へのポインタである[F,D,MPF]Arrayだけを使います。
typedef struct {
float *array; /* array */
long int size; /* Length of array */
} farray;
typedef fstack *FArray;
配列は，DKA法の関数群で利用されています。
UNIX上では
times関数を，Windows上ではQueryPerformanceCounter関数を使って，時間計測を行います。UNIX上では，flag!=0の時はtms構造体の内容を全て標準出力に表示します。get_secv関数はget_sec(0)を呼び出します。
ファイルを介して行列，ベクトル，多項式を読み書きする際に使用する関数群です。データは一行ごとに記述され，“index, 要素データ"の順に並んでいなければなりません。要素データは全てASCIIテキスト形式の10進浮動小数点数で表示されていなければなりません。具体的には次のようになります。
行列データは
行番号i, 列番号j, ij要素
となり，ベクトルデータは
番号i, i要素
となり，多項式データは
次数i, i次係数
となります。
openされたファイルポインタfpから行列あるいはベクトルの要素データを読み取り，matあるいはvecに格納します。ファイル内の要素データは全てASCIIテキストで表現された10進浮動小数点数形式で記述されていなければなりません。
fnameというファイル名を持つファイルをopenし，要素データを読み出します。
matあるいはvecを，openされているファイルfpへ，10進浮動小数点形式のASCIIテキストで書き込みます。
fnameというファイル名のファイルをopenし，matあるいはvecを10進浮動小数点形式のASCIIテキストで書き込みます。
openされたファイルポインタfpから多項式の実係数データを最大次数maxdeg分まで読み取り，polyに格納します。ファイル内の係数データは全てASCIIテキストで表現された10進浮動小数点数形式で記述されていなければなりません。
fnameというファイル名を持つファイルから多項式の係数データを読み取ります。
ropとopを比較して最大値，最小値を返します。
ropのop乗を計算します。
BNCpackのmpf_...関数群で使用するデフォルトの精度(仮数部のbit数)をセットします。この関数を呼び出した後のmpf_...関数が影響を受けます。MPFRパッケージを用いる場合はここでデフォルトの丸めモード(RN)もセットします。
デフォルトの精度を10進桁数で指定し，丸めモードをRNにセットします。実際に指定されるのはropの精度を確保できる最小のbit数です。
MPFRパッケージを読み込んでいる場合のみ，有効となる関数です。 BNCpackのmpf_...関数群で使用するデフォルトの丸めモードをセットします。使用できる丸めモードは次の４種類です。
この関数を呼び出した後のmpf_...関数が影響を受けます。
- GMP_RNDN: Round to Nearest
- GMP_RNDU: Round to +Infinity
- GMP_RNDD: Round to -Infinity
- GMP_RNDZ: Round to Zero
円周率と自然対数の底の値を返します。MPFRパッケージを使用する場合は，そこで定義されている初等関数を使用することになります。
床関数です。GMP Ver.3.xよりサポートされましたが，Ver.2.x以前でも有効な関数として作成してあります。
mpf_t型を単精度，倍精度に変換します。
上から順に，正弦関数，余弦関数，指数関数，対数関数，常用対数関数です。Maclaurin展開に基づくルーチンですので，ものすごく速度が遅いです。高速な多倍長初等関数群が必要でしたらmpfrパッケージの使用をご検討下さい。MPFRパッケージを使うことで，これらの関数は全て高速なものに置き換えられます。また，MPFRで定義されている他の初等関数・特殊関数群もmpf_という接頭語で呼び出せます。
スタックを初期化します。多倍長スタックは
set_bnc_default_prec関数で定義された精度で初期化されます。
スタックを消去します。
値をスタックに積みます。
スタックから値を取り出します。
複素数を格納する領域を確保します。
複素数opの格納されている領域を解放します。
複素数opの実数部を返します。多倍長型はropに実数部が格納されます。
複素数opの虚数部を返します。多倍長型は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 opを計算し，結果をropに戻します。
複素数同士の除算op1 / op2を計算し，結果をropに格納します。
\exp(\varop\times i)を計算し，複素数ropに格納します。
複素数opの値を標準出力に表示します。
array_size個のデータが格納できる配列を確保します。
array_size個の，prec bitの精度を持つデータを格納できる配列を確保します。
配列ropの記憶領域を削除します。
配列opのindex番目のデータを返します。
配列ropのindex番目にデータvalを格納します。
配列opの内容を配列ropに代入します。
配列opの内容を標準出力に表示します。
最大次数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に格納します。
2つの多項式poly1とpoly2の積をrpolyに格納します。
多項式polyを多項式rpolyに代入します。
多項式polyの係数を全て0にセットします。
多項式polyの絶対値最大係数の次数を返します。
多項式polyの非零係数の個数を返します。
多項式polyの導関数を計算し，polyに導関数を格納します。
多項式polyの一変数が実数valである時の値を計算して返します。多倍長多項式の場合はevalにその値を格納します。
多項式polyの一変数が実数valである時，その導関数の値を計算して返します。多倍長多項式の場合はevalにその値を格納します。
多項式polyの一変数が複素数valである時の値を計算し，evalにその値を格納します。
多項式polyの一変数が複素数valである時，その導関数の値を計算し，evalに格納します。
ここに取り上げた関数のうち，IEEE754 double及び多倍長浮動小数点数要素のベクトル・行列を扱うものは，MPIによる並列化がなされています。詳細はmpiディレクトリのソースプログラム，もしくはMPIBNCpackマニュアルをご覧下さい。
ベクトルを初期化します。多倍長ベクトルは
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の転置をrmatに代入します。
行列matの逆行列を計算します。matは正方行列でなければなりません。
行列matのFrobeniusノルム, 無限大ノルム, 1ノルムを計算します。
行列matのLU分解を行います。
LU分解された行列matと定数ベクトルvecを使い，後退代入を行って解ベクトルをrvecに代入していきます。
行列matの部分ピボット選択付きLU分解を行います。ピボット選択後の行の状態は配列iarrayに保存されます。
部分ピボット選択付きでLU分解された行列matと定数ベクトルvecを使い，後退代入を行って解ベクトルをrvecに代入していきます。
行列matの完全ピボット選択付きLU分解を行います。ピボット選択後の行の状態は配列row_iarrayに，列の状態は配列col_iarrayに保存されます。
完全ピボット選択付きでLU分解された行列matと定数ベクトルvecを使い，後退代入を行って解ベクトルをrvecに代入していきます。
残差ベクトルb - A * xを計算します。
以下で紹介する関数はまだテスト段階のものです。実行すると，収束過程が全て標準出力に表示されます。
一般の正方行列を係数行列に持つ連立一次方程式に対し，Jacobi反復法を実行します。係数行列をmatに，定数ベクトルをvecに入れておき，残差のノルムを使った停止条件をrepsとaepsに指定します。最大反復回数はmaxtimesで指定します。
一般の正方行列を係数行列に持つ連立一次方程式に対し，Gauss-Seidel法を実行します。引数の意味は全てJacobi反復法の関数と同じです。
一般の正方行列を係数行列に持つ連立一次方程式に対し，加速係数omegaのSOR法を実行します。その他の引数の意味は全てJacobi反復法の関数と同じです。
ここに取り上げた関数のうち，IEEE754 double及び多倍長浮動小数点数要素の連立一次方程式を扱うものは，MPIによる並列化がなされています。詳細はmpiディレクトリのソースプログラム，もしくはMPIBNCpackマニュアルをご覧下さい。
正定値対称行列を係数行列に持つ連立一次方程式に対し，Conjugate-Gradient法を実行します。係数行列をmatに，定数ベクトルをvecに入れておき，残差のノルムを使った停止条件をrepsとaepsに指定します。最大反復回数はmaxtimesで指定します。正常に終了する際には，返り値は反復回数になります。
一般の正方行列を係数行列に持つ連立一次方程式に対し，BiCG法を実行します。係数行列をmatに，定数ベクトルをvecに入れておき，残差のノルムを使った停止条件をrepsとaepsに指定します。最大反復回数はmaxtimesで指定します。
一般の正方行列を係数行列に持つ連立一次方程式に対し，CGS法を実行します。引数の意味は全てBiCG法の関数と同じです。
一般の正方行列を係数行列に持つ連立一次方程式に対し，BiCGSTAB法を実行します。
一般の正方行列を係数行列に持つ連立一次方程式に対し，GPBiCG法を実行します。
多次元非線型方程式func(x) = 0に対し，初期値をx_initに設定したNewton法を実行します。反復の際に使用するfuncのJacobi行列はjfuncで計算します。収束判定は相対許容度reps, 絶対許容度aepsを用いて行います。収束しない場合は最大反復回数maxtimesに達した時に停止します。
1次元非線型方程式func(x) = 0に対してNewton法を実行します。func関数の導関数はdfuncに指定します。
多次元非線型方程式に対し，簡易Newton法を実行します。
1次元非線型方程式に対し，簡易Newton法を実行します。
ここに取り上げた関数のうち，IEEE754 double及び多倍長浮動小数点数係数の代数方程式を扱うものは，MPIによる並列化がなされています。詳細はmpiディレクトリのソースプログラム，もしくはMPIBNCpackマニュアルをご覧下さい。
実係数多項式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に近似解を格納します。返り値は反復回数です。
テストプログラムの内容は以下の通りです。
コンパイル方法は，前述の通り
%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の関数群の使い方を紹介します。
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型に置き換えられます6。
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
}
最後に，使用した変数を解放します。
BNCpackで定義した複素数型は独自のもので，C++標準のcomplexクラスとも，GSL(http://sources.redhat.com/gsl/)のそれとも互換性はありません。
とある先生からはC++のcomplexクラスでmpf_tもしくはmpfr_t型が利用できるようにならないか？というご要望を受けています。テンプレートですから，ちょっと苦労すれば出来そうな気もしますが・・・どなたかやりません？
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);
最後に，確保した複素数型の変数領域を解放します。
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);
最後に，確保した複素数型の変数領域を解放します。
#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);
}
最後に多項式型変数領域を解放します。
#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進数のビットパターンがそのまま転写されますので，精度は倍精度程度で止まってしまうことに留意して下さい7。
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);
}
最後に，使用した変数領域を解放します。
線型計算のサンプルは次のLU分解，CG法の例で代替します。
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);
最後に，使用した変数領域を解放します。
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);
最後に，使用した変数領域を解放します。
/* 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);
使用した変数領域を解放します。
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);
最後に使用した変数領域を解放します。
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法を実行してその結果を出力します。
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法を実行してその結果を出力します。
/* 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);
最後に，使用した変数を解放します。
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);
最後に，使用した変数を解放します。
今後の予定としては
といったものがありますが，気が向かないと何もしない我が儘な性格なので，あまりあてにしないで下さい。
ご要望を頂いてもすぐにお答えできるかどうか分かりませんが，もし何かご意見などありましたら，マニュアル表紙のメールアドレスまでお願い致します。
GNU MPとMPFRなかりせば本ライブラリは存在しませんでした。開発チーム一同に感謝いたします。
本マニュアルの作成にあたり，角藤氏のW32TeXパッケージのTexinfoを使用致しました。神業かと思うような迅速なupdateを継続してこられたことに，敬意の意と感謝の念を表します。
プリンタがトラブって以来，一貫してFSFを率いつつ，「自由なソフトウェア」(Free Software)の普及活動に邁進し続けているRMSにも感謝いたします。GNU Projectなかりせば，開発環境に使っているLinuxも存在しなかったでしょうから。
Copyright © 2000,2001,2002 Free Software Foundation, Inc.
51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA
Everyone is permitted to copy and distribute verbatim copies
of this license document, but changing it is not allowed.
The purpose of this License is to make a manual, textbook, or other functional and useful document free in the sense of freedom: to assure everyone the effective freedom to copy and redistribute it, with or without modifying it, either commercially or noncommercially. Secondarily, this License preserves for the author and publisher a way to get credit for their work, while not being considered responsible for modifications made by others.
This License is a kind of “copyleft”, which means that derivative works of the document must themselves be free in the same sense. It complements the GNU General Public License, which is a copyleft license designed for free software.
We have designed this License in order to use it for manuals for free software, because free software needs free documentation: a free program should come with manuals providing the same freedoms that the software does. But this License is not limited to software manuals; it can be used for any textual work, regardless of subject matter or whether it is published as a printed book. We recommend this License principally for works whose purpose is instruction or reference.
This License applies to any manual or other work, in any medium, that contains a notice placed by the copyright holder saying it can be distributed under the terms of this License. Such a notice grants a world-wide, royalty-free license, unlimited in duration, to use that work under the conditions stated herein. The “Document”, below, refers to any such manual or work. Any member of the public is a licensee, and is addressed as “you”. You accept the license if you copy, modify or distribute the work in a way requiring permission under copyright law.
A “Modified Version” of the Document means any work containing the Document or a portion of it, either copied verbatim, or with modifications and/or translated into another language.
A “Secondary Section” is a named appendix or a front-matter section of the Document that deals exclusively with the relationship of the publishers or authors of the Document to the Document's overall subject (or to related matters) and contains nothing that could fall directly within that overall subject. (Thus, if the Document is in part a textbook of mathematics, a Secondary Section may not explain any mathematics.) The relationship could be a matter of historical connection with the subject or with related matters, or of legal, commercial, philosophical, ethical or political position regarding them.
The “Invariant Sections” are certain Secondary Sections whose titles are designated, as being those of Invariant Sections, in the notice that says that the Document is released under this License. If a section does not fit the above definition of Secondary then it is not allowed to be designated as Invariant. The Document may contain zero Invariant Sections. If the Document does not identify any Invariant Sections then there are none.
The “Cover Texts” are certain short passages of text that are listed, as Front-Cover Texts or Back-Cover Texts, in the notice that says that the Document is released under this License. A Front-Cover Text may be at most 5 words, and a Back-Cover Text may be at most 25 words.
A “Transparent” copy of the Document means a machine-readable copy, represented in a format whose specification is available to the general public, that is suitable for revising the document straightforwardly with generic text editors or (for images composed of pixels) generic paint programs or (for drawings) some widely available drawing editor, and that is suitable for input to text formatters or for automatic translation to a variety of formats suitable for input to text formatters. A copy made in an otherwise Transparent file format whose markup, or absence of markup, has been arranged to thwart or discourage subsequent modification by readers is not Transparent. An image format is not Transparent if used for any substantial amount of text. A copy that is not “Transparent” is called “Opaque”.
Examples of suitable formats for Transparent copies include plain ascii without markup, Texinfo input format, LaTeX input format, SGML or XML using a publicly available DTD, and standard-conforming simple HTML, PostScript or PDF designed for human modification. Examples of transparent image formats include PNG, XCF and JPG. Opaque formats include proprietary formats that can be read and edited only by proprietary word processors, SGML or XML for which the DTD and/or processing tools are not generally available, and the machine-generated HTML, PostScript or PDF produced by some word processors for output purposes only.
The “Title Page” means, for a printed book, the title page itself, plus such following pages as are needed to hold, legibly, the material this License requires to appear in the title page. For works in formats which do not have any title page as such, “Title Page” means the text near the most prominent appearance of the work's title, preceding the beginning of the body of the text.
A section “Entitled XYZ” means a named subunit of the Document whose title either is precisely XYZ or contains XYZ in parentheses following text that translates XYZ in another language. (Here XYZ stands for a specific section name mentioned below, such as “Acknowledgements”, “Dedications”, “Endorsements”, or “History”.) To “Preserve the Title” of such a section when you modify the Document means that it remains a section “Entitled XYZ” according to this definition.
The Document may include Warranty Disclaimers next to the notice which states that this License applies to the Document. These Warranty Disclaimers are considered to be included by reference in this License, but only as regards disclaiming warranties: any other implication that these Warranty Disclaimers may have is void and has no effect on the meaning of this License.
You may copy and distribute the Document in any medium, either commercially or noncommercially, provided that this License, the copyright notices, and the license notice saying this License applies to the Document are reproduced in all copies, and that you add no other conditions whatsoever to those of this License. You may not use technical measures to obstruct or control the reading or further copying of the copies you make or distribute. However, you may accept compensation in exchange for copies. If you distribute a large enough number of copies you must also follow the conditions in section 3.
You may also lend copies, under the same conditions stated above, and you may publicly display copies.
If you publish printed copies (or copies in media that commonly have printed covers) of the Document, numbering more than 100, and the Document's license notice requires Cover Texts, you must enclose the copies in covers that carry, clearly and legibly, all these Cover Texts: Front-Cover Texts on the front cover, and Back-Cover Texts on the back cover. Both covers must also clearly and legibly identify you as the publisher of these copies. The front cover must present the full title with all words of the title equally prominent and visible. You may add other material on the covers in addition. Copying with changes limited to the covers, as long as they preserve the title of the Document and satisfy these conditions, can be treated as verbatim copying in other respects.
If the required texts for either cover are too voluminous to fit legibly, you should put the first ones listed (as many as fit reasonably) on the actual cover, and continue the rest onto adjacent pages.
If you publish or distribute Opaque copies of the Document numbering more than 100, you must either include a machine-readable Transparent copy along with each Opaque copy, or state in or with each Opaque copy a computer-network location from which the general network-using public has access to download using public-standard network protocols a complete Transparent copy of the Document, free of added material. If you use the latter option, you must take reasonably prudent steps, when you begin distribution of Opaque copies in quantity, to ensure that this Transparent copy will remain thus accessible at the stated location until at least one year after the last time you distribute an Opaque copy (directly or through your agents or retailers) of that edition to the public.
It is requested, but not required, that you contact the authors of the Document well before redistributing any large number of copies, to give them a chance to provide you with an updated version of the Document.
You may copy and distribute a Modified Version of the Document under the conditions of sections 2 and 3 above, provided that you release the Modified Version under precisely this License, with the Modified Version filling the role of the Document, thus licensing distribution and modification of the Modified Version to whoever possesses a copy of it. In addition, you must do these things in the Modified Version:
If the Modified Version includes new front-matter sections or appendices that qualify as Secondary Sections and contain no material copied from the Document, you may at your option designate some or all of these sections as invariant. To do this, add their titles to the list of Invariant Sections in the Modified Version's license notice. These titles must be distinct from any other section titles.
You may add a section Entitled “Endorsements”, provided it contains nothing but endorsements of your Modified Version by various parties—for example, statements of peer review or that the text has been approved by an organization as the authoritative definition of a standard.
You may add a passage of up to five words as a Front-Cover Text, and a passage of up to 25 words as a Back-Cover Text, to the end of the list of Cover Texts in the Modified Version. Only one passage of Front-Cover Text and one of Back-Cover Text may be added by (or through arrangements made by) any one entity. If the Document already includes a cover text for the same cover, previously added by you or by arrangement made by the same entity you are acting on behalf of, you may not add another; but you may replace the old one, on explicit permission from the previous publisher that added the old one.
The author(s) and publisher(s) of the Document do not by this License give permission to use their names for publicity for or to assert or imply endorsement of any Modified Version.
You may combine the Document with other documents released under this License, under the terms defined in section 4 above for modified versions, provided that you include in the combination all of the Invariant Sections of all of the original documents, unmodified, and list them all as Invariant Sections of your combined work in its license notice, and that you preserve all their Warranty Disclaimers.
The combined work need only contain one copy of this License, and multiple identical Invariant Sections may be replaced with a single copy. If there are multiple Invariant Sections with the same name but different contents, make the title of each such section unique by adding at the end of it, in parentheses, the name of the original author or publisher of that section if known, or else a unique number. Make the same adjustment to the section titles in the list of Invariant Sections in the license notice of the combined work.
In the combination, you must combine any sections Entitled “History” in the various original documents, forming one section Entitled “History”; likewise combine any sections Entitled “Acknowledgements”, and any sections Entitled “Dedications”. You must delete all sections Entitled “Endorsements.”
You may make a collection consisting of the Document and other documents released under this License, and replace the individual copies of this License in the various documents with a single copy that is included in the collection, provided that you follow the rules of this License for verbatim copying of each of the documents in all other respects.
You may extract a single document from such a collection, and distribute it individually under this License, provided you insert a copy of this License into the extracted document, and follow this License in all other respects regarding verbatim copying of that document.
A compilation of the Document or its derivatives with other separate and independent documents or works, in or on a volume of a storage or distribution medium, is called an “aggregate” if the copyright resulting from the compilation is not used to limit the legal rights of the compilation's users beyond what the individual works permit. When the Document is included in an aggregate, this License does not apply to the other works in the aggregate which are not themselves derivative works of the Document.
If the Cover Text requirement of section 3 is applicable to these copies of the Document, then if the Document is less than one half of the entire aggregate, the Document's Cover Texts may be placed on covers that bracket the Document within the aggregate, or the electronic equivalent of covers if the Document is in electronic form. Otherwise they must appear on printed covers that bracket the whole aggregate.
Translation is considered a kind of modification, so you may distribute translations of the Document under the terms of section 4. Replacing Invariant Sections with translations requires special permission from their copyright holders, but you may include translations of some or all Invariant Sections in addition to the original versions of these Invariant Sections. You may include a translation of this License, and all the license notices in the Document, and any Warranty Disclaimers, provided that you also include the original English version of this License and the original versions of those notices and disclaimers. In case of a disagreement between the translation and the original version of this License or a notice or disclaimer, the original version will prevail.
If a section in the Document is Entitled “Acknowledgements”, “Dedications”, or “History”, the requirement (section 4) to Preserve its Title (section 1) will typically require changing the actual title.
You may not copy, modify, sublicense, or distribute the Document except as expressly provided for under this License. Any other attempt to copy, modify, sublicense or distribute the Document is void, and will automatically terminate your rights under this License. However, parties who have received copies, or rights, from you under this License will not have their licenses terminated so long as such parties remain in full compliance.
The Free Software Foundation may publish new, revised versions of the GNU Free Documentation License from time to time. Such new versions will be similar in spirit to the present version, but may differ in detail to address new problems or concerns. See http://www.gnu.org/copyleft/.
Each version of the License is given a distinguishing version number. If the Document specifies that a particular numbered version of this License “or any later version” applies to it, you have the option of following the terms and conditions either of that specified version or of any later version that has been published (not as a draft) by the Free Software Foundation. If the Document does not specify a version number of this License, you may choose any version ever published (not as a draft) by the Free Software Foundation.
To use this License in a document you have written, include a copy of the License in the document and put the following copyright and license notices just after the title page:
Copyright (C) year your name.
Permission is granted to copy, distribute and/or modify this document
under the terms of the GNU Free Documentation License, Version 1.2
or any later version published by the Free Software Foundation;
with no Invariant Sections, no Front-Cover Texts, and no Back-Cover
Texts. A copy of the license is included in the section entitled ``GNU
Free Documentation License''.
If you have Invariant Sections, Front-Cover Texts and Back-Cover Texts, replace the “with...Texts.” line with this:
with the Invariant Sections being list their titles, with
the Front-Cover Texts being list, and with the Back-Cover Texts
being list.
If you have Invariant Sections without Cover Texts, or some other combination of the three, merge those two alternatives to suit the situation.
If your document contains nontrivial examples of program code, we recommend releasing these examples in parallel under your choice of free software license, such as the GNU General Public License, to permit their use in free software.
abs_dcmplx: Complex functionsabs_fcmplx: Complex functionsabs_mpfcmplx: Complex functionsadd2_dcmplx: Complex functionsadd2_dvector: Linear functionsadd2_fcmplx: Complex functionsadd2_fvector: Linear functionsadd2_mpfcmplx: Complex functionsadd2_mpfvector: Linear functionsadd_dcmplx: Complex functionsadd_dcmplx_d: Complex functionsadd_dmatrix: Linear functionsadd_dpoly: Polynomial functionsadd_dvector: Linear functionsadd_fcmplx: Complex functionsadd_fcmplx_f: Complex functionsadd_fmatrix: Linear functionsadd_fpoly: Polynomial functionsadd_fvector: Linear functionsadd_mpfcmplx: Complex functionsadd_mpfcmplx_mpf: Complex functionsadd_mpfmatrix: Linear functionsadd_mpfpoly: Polynomial functionsadd_mpfvector: Linear functionsceval_diff_dpoly: Polynomial functionsceval_diff_fpoly: Polynomial functionsceval_diff_mpfpoly: Polynomial functionsceval_dpoly: Polynomial functionsceval_fpoly: Polynomial functionsceval_mpfpoly: Polynomial functionscmul2_dvector: Linear functionscmul2_fvector: Linear functionscmul2_mpfvector: Linear functionscmul_dpoly: Polynomial functionscmul_dvector: Linear functionscmul_fpoly: Polynomial functionscmul_fvector: Linear functionscmul_mpfpoly: Polynomial functionscmul_mpfvector: Linear functionsconj_dcmplx: Complex functionsconj_fcmplx: Complex functionsconj_mpfcmplx: Complex functionsDBiCG: CGDBiCGSTAB: CGDCG: CGDCGS: CGddka: DKAddka_center: DKAddka_init: DKAddka_radius: DKADGPBiCG: CGdgs: Iterationdiff_dpoly: Polynomial functionsdiff_fpoly: Polynomial functionsdiff_mpfpoly: Polynomial functionsdiv_dcmplx: Complex functionsdiv_fcmplx: Complex functionsdiv_mpfcmplx: Complex functionsdjacobi: IterationDLUdecomp: LUDLUdecompC: LUDLUdecompP: LUdmax: Elementary functionsdmin: Elementary functionsdnewton: Newtondnewton_1: Newtondpower: Elementary functionsdsnewton: Newtondsnewton_1: Newtondsor: Iterationeval_diff_dpoly: Polynomial functionseval_diff_fpoly: Polynomial functionseval_diff_mpfpoly: Polynomial functionseval_dpoly: Polynomial functionseval_fpoly: Polynomial functionseval_mpfpoly: Polynomial functionsFCG: CGfdka: DKAfdka_center: DKAfdka_init: DKAfdka_radius: DKAfget_sec: Time measurementfget_secv: Time measurementfgs: Iterationfjacobi: IterationFLUdecomp: LUFLUdecompC: LUFLUdecompP: LUfmax: Elementary functionsfmin: Elementary functionsfnewton: Newtonfnewton_1: Newtonfpower: Elementary functionsfread_dmatrix: IO functionsfread_dmatrix_fname: IO functionsfread_dpolycoef: IO functionsfread_dpolycoef_fname: IO functionsfread_dvector: IO functionsfread_dvector_fname: IO functionsfread_mpfmatrix: IO functionsfread_mpfmatrix_fname: IO functionsfread_mpfpolycoef: IO functionsfread_mpfpolycoef_fname: IO functionsfread_mpfvector: IO functionsfread_mpfvector_fname: IO functionsfree_cdarray: Array functionsfree_cfarray: Array functionsfree_cmpfarray: Array functionsfree_darray: Array functionsfree_dcmplx: Complex functionsfree_dmatrix: Linear functionsfree_dpoly: Polynomial functionsfree_dstack: Stack functionsfree_dvector: Linear functionsfree_farray: Array functionsfree_fcmplx: Complex functionsfree_fmatrix: Linear functionsfree_fpoly: Polynomial functionsfree_fstack: Stack functionsfree_fvector: Linear functionsfree_mpfarray: Array functionsfree_mpfcmplx: Complex functionsfree_mpfmatrix: Linear functionsfree_mpfpoly: Polynomial functionsfree_mpfstack: Stack functionsfree_mpfvector: Linear functionsfsnewton: Newtonfsnewton_1: Newtonfsor: Iterationfwrite_dfvector: IO functionsfwrite_dmatrix: IO functionsfwrite_dmatrix_fname: IO functionsfwrite_dvector_fname: IO functionsfwrite_mpfmatrix: IO functionsfwrite_mpfmatrix_fname: IO functionsfwrite_mpfvector: IO functionsfwrite_mpfvector_fname: IO functionsgdmij: Linear functionsgdvi: Linear functionsget_bnc_default_prec: Elementary functionsget_cdarray_i: Array functionsget_cfarray_i: Array functionsget_cmpfarray_i: Array functionsget_darray_i: Array functionsget_dmatrix_ij: Linear functionsget_dpoly_i: Polynomial functionsget_dvector_i: Linear functionsget_farray_i: Array functionsget_fmatrix_ij: Linear functionsget_fpoly_i: Polynomial functionsget_fvector_i: Linear functionsget_image_dcmplx: Complex functionsget_image_fcmplx: Complex functionsget_image_mpfcmplx: Complex functionsget_mpfarray_i: Array functionsget_mpfmatrix_ij: Linear functionsget_mpfpoly_i: Polynomial functionsget_mpfvector_i: Linear functionsget_real_dcmplx: Complex functionsget_real_fcmplx: Complex functionsget_real_mpfcmplx: Complex functionsget_residual_dvector: Iterationget_residual_fvector: Iterationget_residual_mpfvector: Iterationget_sec: Time measurementget_secv: Time measurementgfmij: Linear functionsgfvi: Linear functionsgmpfmij: Linear functionsgmpfvi: Linear functionsiexp_dcmplx: Complex functionsiexp_fcmplx: Complex functionsiexp_mpfcmplx: Complex functionsinit2_cmpfarray: Array functionsinit2_mpfarray: Array functionsinit2_mpfcmplx: Complex functionsinit2_mpfmatrix: Linear functionsinit2_mpfpoly: Polynomial functionsinit2_mpfstack: Stack functionsinit2_mpfvector: Linear functionsinit_cdarray: Array functionsinit_cfarray: Array functionsinit_cmpfarray: Array functionsinit_darray: Array functionsinit_dcmplx: Complex functionsinit_dmatrix: Linear functionsinit_dpoly: Polynomial functionsinit_dstack: Stack functionsinit_dvector: Linear functionsinit_farray: Array functionsinit_fcmplx: Complex functionsinit_fmatrix: Linear functionsinit_fpoly: Polynomial functionsinit_fstack: Stack functionsinit_fvector: Linear functionsinit_mpfarray: Array functionsinit_mpfcmplx: Complex functionsinit_mpfmatrix: Linear functionsinit_mpfpoly: Polynomial functionsinit_mpfstack: Stack functionsinit_mpfvector: Linear functionsinv_dmatrix: Linear functionsinv_fmatrix: Linear functionsinv_mpfmatrix: Linear functionsip_dvector: Linear functionsip_fvector: Linear functionsip_mpfvector: Linear functionslmax: Elementary functionslmin: Elementary functionslpower: Elementary functionsmax_abscoef_dpoly: Polynomial functionsmax_abscoef_fpoly: Polynomial functionsmax_abscoef_mpfpoly: Polynomial functionsmaxprec_mpfmatrix: Linear functionsmaxprec_mpfvector: Linear functionsminprec_mpfmatrix: Linear functionsminprec_mpfvector: Linear functionsmpf2double: Elementary functionsmpf2float: Elementary functionsmpf_cos: Elementary functionsmpf_dka: DKAmpf_dka_center: DKAmpf_dka_init: DKAmpf_dka_radius: DKAmpf_e: Elementary functionsmpf_exp: Elementary functionsmpf_floor: Elementary functionsmpf_gs: Iterationmpf_jacobi: Iterationmpf_ln: Elementary functionsmpf_log10: Elementary functionsmpf_max: Elementary functionsmpf_min: Elementary functionsmpf_newton: Newtonmpf_newton_1: Newtonmpf_pi: Elementary functionsmpf_power: Elementary functionsmpf_sin: Elementary functionsmpf_snewton: Newtonmpf_snewton_1: Newtonmpf_sor: IterationMPFBiCG: CGMPFBiCGSTAB: CGMPFCG: CGMPFCGS: CGMPFGPBiCG: CGMPFLUdecomp: LUMPFLUdecompC: LUMPFLUdecompP: LUmul2_dcmplx: Complex functionsmul2_fcmplx: Complex functionsmul2_mpfcmplx: Complex functionsmul_dcmplx: Complex functionsmul_dcmplx_d: Complex functionsmul_dmatrix: Linear functionsmul_dmatrix_dvec: Linear functionsmul_dpoly: Polynomial functionsmul_fcmplx: Complex functionsmul_fcmplx_f: Complex functionsmul_fmatrix: Linear functionsmul_fmatrix_fvec: Linear functionsmul_mpfcmplx: Complex functionsmul_mpfcmplx_mpf: Complex functionsmul_mpfcmplx_ui: Complex functionsmul_mpfmatrix: Linear functionsmul_mpfmatrix_mpfvec: Linear functionsmul_mpfpoly: Polynomial functionsnorm1_dmatrix: Linear functionsnorm1_dvector: Linear functionsnorm1_fvector: Linear functionsnorm1_mpfmatrix: Linear functionsnorm1_mpfvector: Linear functionsnorm2_dvector: Linear functionsnorm2_fvector: Linear functionsnorm2_mpfvector: Linear functionsnormf_dmatrix: Linear functionsnormf_mpfmatrix: Linear functionsnormi_dmatrix: Linear functionsnormi_dvector: Linear functionsnormi_fvector: Linear functionsnormi_mpfmatrix: Linear functionsnormi_mpfvector: Linear functionsnum_nonzero_dpoly: Polynomial functionsnum_nonzero_fpoly: Polynomial functionsnum_nonzero_mpfpoly: Polynomial functionspop_dstack: Stack functionspop_fstack: Stack functionspop_mpfstack: Stack functionsprec_mpfmatrix: Linear functionsprec_mpfvector: Linear functionsprint_cdarray: Array functionsprint_cfarray: Array functionsprint_cmpfarray: Array functionsprint_darray: Array functionsprint_dcmplx: Complex functionsprint_dmatrix: Linear functionsprint_dpoly: Polynomial functionsprint_dvector: Linear functionsprint_farray: Array functionsprint_fcmplx: Complex functionsprint_fdmpfpoly: Polynomial functionsprint_fmatrix: Linear functionsprint_fpoly: Polynomial functionsprint_fvector: Linear functionsprint_mpfarray: Array functionsprint_mpfcmplx: Complex functionsprint_mpfmatrix: Linear functionsprint_mpfpoly: Polynomial functionsprint_mpfvector: Linear functionspush_dstack: Stack functionspush_fstack: Stack functionspush_mpfstack: Stack functionssdmij: Linear functionssdvi: Linear functionsset0_dcmplx: Complex functionsset0_dmatrix: Linear functionsset0_dpoly: Polynomial functionsset0_fcmplx: Complex functionsset0_fmatrix: Linear functionsset0_fpoly: Polynomial functionsset0_mpfcmplx: Complex functionsset0_mpfmatrix: Linear functionsset0_mpfpoly: Polynomial functionsset_bnc_default_prec: Elementary functionsset_bnc_default_prec_decmal: Elementary functionsset_bnc_rounding_mode: Elementary functionsset_cdarray_i: Array functionsset_cfarray_i: Array functionsset_cmpfarray_i: Array functionsset_darray_i: Array functionsset_dmatrix_ij: Linear functionsset_dpoly_i: Polynomial functionsset_dvector_i: Linear functionsset_farray_i: Array functionsset_fmatrix_ij: Linear functionsset_fpoly_i: Polynomial functionsset_fvector_i: Linear functionsset_image_dcmplx: Complex functionsset_image_fcmplx: Complex functionsset_image_mpfcmplx: Complex functionsset_image_mpfcmplx_ui: Complex functionsset_mpfarray_i: Array functionsset_mpfmatrix_ij: Linear functionsset_mpfmatrix_ij_d: Linear functionsset_mpfmatrix_ij_str: Linear functionsset_mpfpoly_i: Polynomial functionsset_mpfpoly_i_d: Polynomial functionsset_mpfpoly_i_str: Polynomial functionsset_mpfvector_i: Linear functionsset_mpfvector_i_d: Linear functionsset_mpfvector_i_str: Linear functionsset_real_dcmplx: Complex functionsset_real_fcmplx: Complex functionsset_real_mpfcmplx: Complex functionsset_real_mpfcmplx_ui: Complex functionssetdegree_dpoly: Polynomial functionssetdegree_fpoly: Polynomial functionssetI_dmatrix: Linear functionssetI_fmatrix: Linear functionssetI_mpfmatrix: Linear functionssfmij: Linear functionssfvi: Linear functionssign2_dcmplx: Complex functionssign2_fcmplx: Complex functionssign2_mpfcmplx: Complex functionssign_dcmplx: Complex functionssign_fcmplx: Complex functionssign_mpfcmplx: Complex functionssmpfmij: Linear functionssmpfmijd: Linear functionssmpfmijs: Linear functionssmpfvi: Linear functionssmpfvid: Linear functionssmpfvis: Linear functionsSolveDLS: LUSolveDLSC: LUSolveDLSP: LUSolveFLS: LUSolveFLSC: LUSolveFLSP: LUSolveMPFLS: LUSolveMPFLSC: LUSolveMPFLSP: LUsub2_dvector: Linear functionssub2_fvector: Linear functionssub2_mpfvector: Linear functionssub_dcmplx: Complex functionssub_dmatrix: Linear functionssub_dpoly: Polynomial functionssub_dvector: Linear functionssub_fcmplx: Complex functionssub_fmatrix: Linear functionssub_fpoly: Polynomial functionssub_fvector: Linear functionssub_mpfcmplx: Complex functionssub_mpfmatrix: Linear functionssub_mpfpoly: Polynomial functionssub_mpfvector: Linear functionssubst_cdarray: Array functionssubst_cfarray: Array functionssubst_cmpfarray: Array functionssubst_darray: Array functionssubst_dcmplx: Complex functionssubst_dmatrix: Linear functionssubst_dpoly: Polynomial functionssubst_dvector: Linear functionssubst_farray: Array functionssubst_fcmplx: Complex functionssubst_fmatrix: Linear functionssubst_fpoly: Polynomial functionssubst_fvector: Linear functionssubst_mpfarray: Array functionssubst_mpfcmplx: Complex functionssubst_mpfmatrix: Linear functionssubst_mpfpoly: Polynomial functionssubst_mpfvector: Linear functionstranspose_dmatrix: Linear functionstranspose_fmatrix: Linear functionstranspose_mpfmatrix: Linear functionsCDArray: Type definitionsCDMatrix: Type definitionsCDVector: Type definitionsCFArray: Type definitionsCFMatrix: Type definitionsCFVector: Type definitionsCMPFArray: Type definitionsCMPFMatrix: Type definitionsCMPFVector: Type definitionsDArray: Type definitionsdcmplx: Type definitionsDMatrix: Type definitionsDPoly: Type definitionsDStack: Type definitionsDVector: Type definitionsFArray: Type definitionsfcmplx: Type definitionsFMatrix: Type definitionsFPoly: Type definitionsFStack: Type definitionsFVector: Type definitionsMPFArray: Type definitionsmpfcmplx: Type definitionsMPFMatrix: Type definitionsMPFPoly: Type definitionsMPFStack: Type definitionsMPFVector: Type definitions[1] 2006年3月現在の最新版はMPFR 2.2.0, GMP 4.2です。URLはそれぞれhttp://www.mpfr.org/, http://swox.com/gmp/です.
[2] 今となってはそれほどでもないでしょうけど，自分としては売りのつもり。
[3] MPFRパッケージの マニュアル参照。
[4] 実は共用できるようにしてありますが，表では使わないようになっています。
[5] GMP単独では，mpf_t型の丸め処理は切り捨てで行われます。
[6] 詳細はMPFRパッケージに同封されているmpf2mpfr.hを参照のこと。
[7] それで何度失敗したことか・・・。