C++ による GNU Scientific Library(GSL) 入門

GNU Scientific Library(GSL) は，種々の数値計算が可能な C 言語用の数値計算ライブラリです．
GSL がカバーしている範囲は多岐に渡るので，数学として理解できる部分を少しずつ増やしていこうという戦術で進めています．
これまでに，初等的な数学関数，一元高次方程式，統計，数値微分，といったトピックスを採りあげ，プログラムを走らせて勉強してみました．

GSL は C 言語用の数値計算ライブラリですが，このページでは，C 言語ではなく，C++ 言語で記述しています．

インフォメーション

ウェブサイトとオンラインドキュメント

Debian パッケージ

このページは，Debian GNU/Linux バージョン 13 "trixie" で作成しました．
関連する主なパッケージを示します．

パッケージ名(バージョン)注
libgsl-dev(2.8)GNU 科学ライブラリ（数値解析用ルーチン集）の開発ファイル
gsl-doc-info(2.7)info 形式のリファレンスマニュアル（未インストール）
gsl-doc-pd(2.7)pdf 形式のリファレンスマニュアル（未インストール）
g++(14.2.0)GNU C++ コンパイラ．このページでは C 言語ではなく C++ 言語で練習しています

目次（ページ内リンク）

オンラインドキュメントで目次としてリストされている大項目と，このページのセクションへのリンクとを対応付けてみました．
オンラインドキュメントの項目は膨大です．未着手の項目が連続している箇所は，'___' で区切って 1 行にまとめました．
これらの項目は，オンラインドキュメントでは，"Introduction" の 最初のセクション "Routines available in GSL" に一覧表としてまとめられています．

Using the Library

オンラインドキュメントの "Introduction" の次のページ "Using the Library" を見ると，最初のセクションは An Example Program となっています．
"An Example Program" では，Bessel 関数の値を計算しています．
いきなり難しそうな材料を扱っていますが，内容は，gsl_sf_bessel_J0(x) という一変数関数の x = 5 のときの戻り値を出力しているだけです．

ソースコード

C 言語のコードを C++ 言語に改変しました．pringf を cout に変えただけですが． 


#include <iostream>
#include <iomanip>
#include <gsl/gsl_sf_bessel.h>  //ヘッダファイルはこれを使う

int main (void) {
  double x = 5.0;
  double y = gsl_sf_bessel_J0(x);  //Bessel 関数の値を計算
  std::cout << "J0(" << x << ") = " << std::scientific << std::setprecision(18) << y << std::endl;
  return 0;
}

ビルドと実行

"An Example Program" の次のパラグラフからは，コンパイルとリンク法が記されています．
C++ 言語用には，ページの後半に Compatibility with C++ というセクションが設けられています．
そこには，次のように記されています．
「GSL の関数は，ヘッダファイルで，C++プログラムにインクルードされると自動的に extern "C " リンケージを持つように定義されています．
これにより，C++から関数を直接呼び出すことができます．
ユーザー定義関数の中で C++ の例外処理を使用するには，CFLAGS のコンパイル・オプション -fexceptions を追加してビルドする必要があります．」

当ページでは例外処理をしていないので，ビルド時のリンクオプションは，-lgsl のみとします（-fexceptions は付けないということ）．

端末で，
$ g++ main.cc -lgsl
としてビルド．a.out という実行ファイルが生成するので，
$ ./a.out
として実行したら， 

J0(5) = -1.775967713143382642e-01

と表示されました．"An Example Program" と同じ値が出力されました．

Mathematical Functions（数学関数）

Mathematical Functionsでは，数学定数，無限大と非数，初等関数，…の順で解説が進みます．
調べていないのですが，標準 C ライブラリに無い関数や定数を中心に開発が進められているようです．補完や代替用ということでしょうか．
もちろん重複している関数や定数もあります．
いくつかのパラグラフから，最初に紹介されている定数や関数を使ってコードを書いてみました． 


#include <iostream>
#include <gsl/gsl_math.h>  //ヘッダファイルはこれを使う

int main (void) {
 //"Mathematical Constants" から．自然対数の底 e の値を出力してみる
 std::cout << "e = " << M_E << std::endl;

 //”Elementary Functions" から．ln(1 + 0.01) を計算してみる
 std::cout << "ln(1 + 0.01) = " << gsl_log1p(0.01) << std::endl;

 //"Small integer powers" から．2 の 2.3 乗を計算してみる
 std::cout << "2^(3/2) = " << gsl_pow_2(2.3) << std::endl;

 return 0;
}

ビルドと実行

$ g++ main.cc -lgsl
$ ./a.out
端末への出力を下に記します． 

e = 2.71828
ln(1 + 0.01) = 0.00995033
2^(3/2) = 5.29

Polynominal（多項式）

Polynominalでは，多項式を評価する（x を与えて値を計算する）関数や，二次方程式，三次方程式の解を求める関数，ソルバーを使って高次方程式を解く方法などが紹介されています．
ここでは，二次方程式とを解く関数と，ソルバーを使って高次方程式を解くコードを紹介します．
後者はドキュメントのサンプルコードを改変したコードです．

ソースコード 1：二次方程式を解く関数

ax2 + bx + c = 0
を解く関数は，
int gsl_poly_solve_quadratic(double a, double b, double c, double *x0, double *x1)
です（"quadratic" は，二次の多項式の意）．
暗算でも解けて検算が楽なサンプルコードを作成してみました． 


#include <iostream>
#include <gsl/gsl_poly.h>  //ヘッダファイル

int main() {
 //解こうとしている方程式は，x^2 + 5x + 6 = 0
 double a = 1.0;  //x^2 の係数
 double b = 5.0;  //x の係数
 double c = 6.0;  //定数
 double x0 = 0.0;  //未知数（解を格納するための変数）．アドレスを渡す
 double x1 = 0.0;  //未知数（解を格納するための変数）．アドレスを渡す

 //未知数 x0 と x1 はアドレスを渡す
 std::cout << "実解の個数：" << gsl_poly_solve_quadratic(a, b, c, &x0, &x1) << std::endl;

 std::cout << "解：" << x0 << ", " << x1 << std::endl;

 return 0;
 }

ビルドと実行

$ g++ main.cc -lgsl
$ ./a.out
としたら，端末に実解の数と解が出力されました．
虚数解が 2 個の二次方程式を供したら，戻り値は 0 で，未知数は 2 個とも 0.0 のままでした． 

実解の個数：2
解：-3, -2

ソースコード 2：ソルバーで高次方程式を解く

下のコードはサンプルコードを少し改変したものです．五次方程式
-1 + 0*x + 0*x2 + 0*x3 + 0*x4 + 1*x5 = 0（すなわち x5 = 1）
を解いています．
この方程式を解くソルバーは，
int gsl_poly_complex_solve(const double *a, size_t n, gsl_poly_complex_workspace *w, gsl_complex_packed_ptr z)
です．引数はそれぞれ，

  1. 係数を格納した配列へのポインタ
  2. 係数の個数（方程式の次数 + 1）
  3. ワークスペース（計算用のメモリ領域）
  4. 複素解の実部と虚部を格納するための配列へのポインタ

ワークスペースを確保するための関数は，
gsl_poly_complex_workspace *gsl_poly_complex_workspace_alloc(size_t n)
です．係数の個数を引数として与えると，ワークスペースへのポインタが返されます． 


#include <iostream>
#include <iomanip>
#include <gsl/gsl_poly.h>

int main (void) {
 double a[6] = { -1, 0, 0, 0, 0, 1 };  //方程式の係数
 double z[10];  //複素解の実部と虚部を格納するための配列

 //ソルバー．解は一般に複素数である．配列 z には係数のみ格納される
 gsl_poly_complex_workspace *w = gsl_poly_complex_workspace_alloc(6);  //ワークスペースを確保
 gsl_poly_complex_solve(a, 6, w, z);  //方程式を解く
 gsl_poly_complex_workspace_free(w);  //ワークスペースを解放


 //解を出力
 std::cout << std::setprecision(18);
 for (int i = 0; i < 5; i++)
  std::cout << "z[" << i << "] = " << z[2*i] << ' ' << z[2*i+1] << std::endl;

 return 0;
}

ビルドと実行

$ g++ main.cc -lgsl
$ ./a.out
端末への出力を下に記します．各行に少数が 2 個並んでいますが，左側が解の実部，右側が虚部です．
ドキュメントと同じ数字が出力されました． 

z[0] = -0.809016994374947673 0.587785252292473359
z[1] = -0.809016994374947673 -0.587785252292473359
z[2] = 0.309016994374947507 0.951056516295152976
z[3] = 0.309016994374947507 -0.951056516295152976
z[4] = 0.999999999999999889 0

Vector and Matrices（ベクトルと行列）

Vectors and Matrices では，GSL でのベクトルと行列を表現するデータ型やこれらを操作する関数が解説されています．線形代数のようなもう少し高度な話題は，別のページで解説されています．

ページ冒頭で，ベクトルと行列は，C 言語の配列を使って表現していることと，ベクトルと行列のメモリ管理をおこなうための型 block（以下「ブロック」）が定義されていると述べられています．
ブロック構造体の本質は，ベクトルや行列の要素を格納するための一次元配列です．ベクトル構造体や行列構造体はブロックへのポインタを 1 個有しています．
ブロック型の理解が，GSL でベクトルと行列を処理する第一歩となるようです．

データ型

いろいろな数値型に対するブロック型（とブロック型を扱うための関数）が定義されています，
基本的には double 型を使うようで，各関数には gsl_block，gsl_vector，gsl_matrix という接頭辞を付けます．他の数値型にはその型を明示する接頭辞を付けます．例えば float 型なら gsl_block_float，int 型なら gsl_block_int といった感じです．
ドキュメントに接頭辞の一覧表があります．

ブロック構造体

gsl_block 構造体（つまり double 型の数値から成るブロック）の宣言が紹介されているので，そのまま引用します． 


type gsl_block

typedef struct
{
 size_t size;
 double *data;
} gsl_block;

ブロック構造体は， 配列のサイズと配列へのポインタから構成されています．
後述しますが，ベクトル構造体と行列構造体は内部にブロックへのポインタ gsl_block *block を有しており，ブロックが有する配列へのポインタ double *data がベクトルや行列の要素となります．

ブロックを扱うための関数

ドキュメントでは，ブロックのメモリアロケーションをおこなうための関数と，ファイルと端末に対する入出力をおこなう関数が紹介されています．
下は，メモリアロケーションをおこなうための関数です．接頭辞から double 型であることが判ります．
関数の名前からは，内部で malloc() と free() を使って *data をアロケートしているように見えます．

サンプルコード 1

ドキュメントにはブロックを生成するサンプルコードが載っています．端末への出力部分を C++ 言語に書き直して引用します．
ブロックを生成して，要素数と配列のアドレスを端末に出力するプログラムです． 


#include <iostream>
#include <gsl/gsl_block.h>

int main (void) {
 gsl_block *b = gsl_block_alloc (100);

 std::cout << "length of block = " << b->size << std::endl;
 std::cout << "block data address = " << b->data << std::endl;

 gsl_block_free(b);
 return 0;
}

ビルドと実行

$ g++ main.cc -lgsl
$ ./a.out

実行結果

length of block = 100
block data address = 0x55561436f2d0

ベクトルの宣言

ベクトルの宣言を引用します．
ブロックへのポインタやメモリ関連のオブジェクトを有する構造体ですが，今のレベルではこれらの扱い方が解りません． 


//ベクトル
typedef struct {
  size_t size;  //ベクトルの要素数
  size_t stride;  //物理メモリにおける要素間のステップサイズ
  double *data;  //ベクトルの先頭要素のメモリ上の位置
  gsl_block *block;  //ベクトルの要素が配置されているメモリブロック
  int owner;  //メモリの解放用フラグ．ブロックが複数のベクトルによって所有されている場合は 0，他のベクトルが所有していない場合は 0
} gsl_vector;

ベクトルの生成と削除

宣言に続いて操作関数が紹介されています．ベクトルでの項目を例に挙げると，

  1. Vector allocation
  2. Accessing vector elements
  3. Initializing vector elements
  4. Reading and writing vectors
  5. Vector views
  6. Copying vectors
  7. Exchanging elements
  8. Vector operations
  9. Finding maximum and minimum elements of vectors
  10. Vector properties

ベクトルの生成と削除

上記 "Vector allocation" にはベクトルの生成と削除を行う関数が 3 個紹介されています．
コンストラクタ引数は要素数だけなので，とりあえずは 0 初期化するかどうか，作成したら消去，に留意すればよさそうです．

サンプルコード 2

下のコードは，ドキュメントに載っているサンプルコードの，端末への出力部分を C++ 言語に書き直したものです．
三次元ベクトル [1.23, 2.23, 3.23] を生成して，要素を端末に出力します．
outo of range エラーが発生するコードですが，例外処理をしていないので三次元ベクトル v のメモリは解放されません． 


#include <iostream>
#include <gsl/gsl_vector.h>

int main(void) {
 int i;
 gsl_vector *v = gsl_vector_alloc(3);  //三次元ベクトルを生成

 for (i = 0; i < 3; i++) {
  gsl_vector_set(v, i, 1.23 + i);
 }

 for(i = 0; i < 100; i++) {  //わざと 100 次元までを出力させようとする．当然，例外が発生する
  std::cout << "v_" << i << " = " << gsl_vector_get (v, i) << std::endl;
 }

 gsl_vector_free (v);
 return 0;
}

ビルドと実行

$ g++ main.cc -lgsl
$ ./a.out

実行結果

予定通りエラーが発生しました． 

v_0 = 1.23
v_1 = 2.23
v_2 = 3.23
gsl: ../gsl/gsl_vector_double.h:182: ERROR: index out of range
v_3 = Default GSL error handler invoked.
中止

行列の宣言

ベクトルの宣言を引用します．
ベクトルと同じく，内部にブロックへのポインタを含んでいます．当然ながら行と列を示すサイズオブジェクトが 2 個あります． 


//行列
typedef struct {
  size_t size1;  //行数
  size_t size2;  //列数
  size_t tda;  //物理次元．メモリ上に配置される行列列数
  double *data;  //行列の要素が格納されているメモリへのポインタ
  gsl_block *block;  //行列ブロックが所有するブロック
  int owner;  //所有権フラグ
} gsl_matrix;

続いて，行列操作関数がいろいろ紹介されています．ベクトルの操作関数と似ており，煩雑になるので省略します．

Statistics（統計）

配列を引数とする統計関数なのですが，留意すべき点を紹介しておきます（Statistics ページの冒頭に書いてあります）．

ヘッダファイル

まず，ヘッダファイルが数値型ごとに用意されている，と言うことです．
サンプルコードのように，
#include <gsl/gsl_statistics.h>
とすれば気にすることはありません．

どのようなヘッダファイルがあるかは，
$ ls /usr/include/gsl/gsl_stat*
とすれば確認できます．
$ lv /usr/include/gsl/gsl_statistics.h
とすると，gsl_statistics.h は，他のヘッダファイルをインクルードしているだけということが解ります．

同じような関数が，仕様を異にして定義されている場合があります．これらは別のヘッダファイルで宣言されます．
例えば，"gsl_statistics_int.h" を見ると，平均値を求める関数は，
double gsl_stats_int_mean (const int data[], const size_t stride, const size_t n);
と宣言されています．この関数は，int 型の配列を引数とし，double 型の平均値を返します．
想像できるように，"gsl_statistics_double.h" というヘッダファイルもあります．

stride

次に留意すべき点は，配列の要素を 2 個おきとか 3 個おきとかの間隔でサンプリングできるという点です．
double 型のデータの平均を求める関数は，
double gsl_stats_mean(const double data[], size_t stride, size_t n)
です．引数は各々，

  1. data[]：データが格納されている配列へのポインタ
  2. stride：読み取るデータの間隔
  3. n：読み取られるデータ数

これらの点を確認しつつ，プログラムを走らせてみます．
上の gsl_stats_mean() という平均を求める関数の他に，分散，最大値と最小値を求める関数も使っています．
これらの関数の引数は，上の gsl_stats_mean() と同じ様式で与えます．
下のコードは，ドキュメントのサンプルコードを改変したものです．

サンプルコード 1


#include <iostream>
#include <gsl/gsl_statistics_double.h>  //ドキュメントのサンプルコードでは，gsl/gsl_statistics.h

int main(void) {
 double data[6] = {1.1, 10.1, 1.2, 10.2, 1.3, 10.3};
 double mean, variance, largest, smallest;

 mean     = gsl_stats_mean(data, 1, 6);
 variance = gsl_stats_variance(data, 1, 6);
 largest  = gsl_stats_max(data, 1, 6);
 smallest = gsl_stats_min(data, 1, 6);

 std::cout << "The sample mean is " << mean << std::endl;
 std::cout << "The estimated variance is " << variance << std::endl;
 std::cout << "The largest value is " << largest << std::endl;
 std::cout << "The smallest value is " << smallest << std::endl;

 return 0;
}

ビルドと実行

$ g++ main.cc -lgsl
$ ./a.out
端末への出力を下に記します．
stride == 1, n= 6 と言うことで，全データの平均値等が得られました． 

The sample mean is 5.7
The estimated variance is 24.308
The largest value is 10.3
The smallest value is 1.1

サンプルコード 2

上のサンプルコード 1 を改変（引数を少し変えただけ）して，配列のデータを一つおきにサンプリングして統計諸量を求めてるコードにしました． 


#include <iostream>
#include <gsl/gsl_statistics_double.h>

int main(void) {
 double data[6] = {1.1, 10.1, 1.2, 10.2, 1.3, 10.3};
 double mean, variance, largest, smallest;

 mean     = gsl_stats_mean(data, 2, 3);
 variance = gsl_stats_variance(data, 2, 3);
 largest  = gsl_stats_max(data, 2, 3);
 smallest = gsl_stats_min(data, 2, 3);

 std::cout << "The sample mean is " << mean << std::endl;
 std::cout << "The estimated variance is " << variance << std::endl;
 std::cout << "The largest value is " << largest << std::endl;
 std::cout << "The smallest value is " << smallest << std::endl;

 return 0;
}

ビルドと実行

$ g++ main.cc -lgsl
$ ./a.out
端末への出力を下に記します．
stride == 2, n= 3 と言うことで，0，2，4 番めの 3 個のデータの平均値等が得られました． 

The sample mean is 1.2
The estimated variance is 0.01
The largest value is 1.3
The smallest value is 1.1

Numerical Differentiation（数値微分）

Numerical Differentiation のページを練習します．
「数値」とついているのは，微小幅（差分）を取って {f(x + Δh) - f(x)} / Δh を計算するからです．ただし，計算のアルゴリズムはもう少し複雑ですが．
サンプルコードでは，y = x3/2 の，x = 2 および x= 0 における微分係数を計算しています．
これに対して，高校数学から使っている，「x のべき乗の値を係数にかけて，べき乗の値から 1 引く」といった感じで微分するのは，解析的な微分などといいます．

差分の取り方には Δh を x のどこに置くかによって，「中心差分」，「前方差分」，「後方差分」があります．
サンプルコードの x = 2 では中心差分を取っています．
サンプルコードの x = 0 では前方差分を取っています．これは，x が負数になると y = x3/2 が意味を成さなくなるからです．
サンプルコードでの差分の大きさは 10-8 です． 


#include <iostream>
#include <iomanip>
#include <gsl/gsl_math.h>
#include <gsl/gsl_deriv.h>

double f(double x, void *params) {
 (void)(params); /* avoid unused parameter warning */
 return pow (x, 1.5);
}

int main(void) {
 gsl_function F;
 double result, abserr;

 F.function = &f;
 F.params = 0;

 std::cout << "f(x) = x^(3/2)" << std::endl;
 std::cout << std::fixed << std::setprecision(8);

 gsl_deriv_central (&F, 2.0, 1e-8, &result, &abserr);
 std::cout << "x = 2.0" << std::endl;
 std::cout << "f'(x) = " << result << " +/- " <<  abserr << std::endl;
 std::cout << "解析的に求めた微分係数 = " << 1.5 * sqrt(2.0) << std::endl;

 gsl_deriv_forward (&F, 0.0, 1e-8, &result, &abserr);
 std::cout << "x = 0.0" << std::endl;
 std::cout << "f'(x) = " << result << " +/- " <<  abserr << std::endl;
 std::cout << "解析的に求めた微分係数 = " << 0.0 << std::endl;

 return 0;
}

ビルドと実行

~$ g++ main.cc -lgsl
~$ ./a.out
としたときの結果を下に示します．有効桁数は異なりますが，ドキュメントの値と一致したようです． 

f(x) = x^(3/2)
x = 2.0
f'(x) = 2.12132031 +/- 0.00000050
解析的に求めた微分係数 = 2.12132034
x = 0.0
f'(x) = 0.00000002 +/- 0.00000003
解析的に求めた微分係数 = 0.00000000

ソルバに渡す関数

上のソースコードでは，y = x3/2 を double f(double x, void *params) という関数で表現しています．
この表現は，One Dimensional Root-Finding というページの "Providing the function to solve" というパラグラフで解説されています．
ソルバーに渡される，解かれるべき関数を表すデータ型です（C 言語は勉強したことないのですが，関数ポインタでいいのかな）．
例として，f(x) = 3x2 + 2x + 1 がどう表現されるかが，ソースコードで示されています．すなわち， 


double my_f(double x, void *p) {
 struct my_f_params *params = (struct my_f_params *)p;
 double a = (params->a);
 double b = (params->b);
 double c = (params->c);

 return  (a * x + b) * x + c;
}

gsl_function F;
struct my_f_params params = {3.0, 2.0, 1.0};

F.function = &my_f;
F.params = &params;