PythonからPyROOT経由でC言語ライブラリを使う具体例としてGSL(GNU Scientific Library)を使ってみる

# 概要

PythonからC言語ライブラリであるGSL(GNU Scientific Library)を試します。試すのはGSLのドキュメントに書かれているexamplesです。
**PyROOT**は**CERN ROOT**のC++インタプリタも使えるためC++のライブラリはそのまま使えますが、Cのライブラリを使う場合`extern "C"`を書かなければなりません。しかし、GSLはC++で使われるときは`extern "C"`が有効になるため、ユーザが`extern "C"`を書かなくても大丈夫です。

なお、Macは**CERN ROOT**をインストールするときは**GSL**も要求されるので、今回のために**GSL**を新たにインストールする必要はありませんでした。ただし、**Linux,Windows**は未確認

### 関連記事

https://qiita.com/h6591/items/4ef61d82454e7e148a85

# 実行環境

```terminal
> sw_vers 
ProductName:		macOS
ProductVersion:		15.3.2
BuildVersion:		24D81

> root --version
ROOT Version: 6.34.06
Built for macosxarm64 on Mar 27 2025, 04:39:07
From tags/6-34-06@6-34-06

> root-config --python-version
3.13.2

> python --version
Python 3.13.2
```


# GSL(GNU Scientific Library)

以下の5つのexamplesを試してみます。使用したexamplesのリンクをつけてました。
さらに、実行結果の後にexamplesのオリジナルのCのコードと結果を折りたたみ形式で載せました。

## 5次方程式の解

 https://www.gnu.org/software/gsl/doc/html/poly.html#examples

$$ x^5 = 1$$

```py:gsl1.py
# PyROOTのimport
import ROOT

import array

# GSLヘッダーファイルのPATH追加
ROOT.gInterpreter.AddIncludePath("/opt/homebrew/Cellar/gsl/2.8/include")
# GSLの動的リンクライブラリのパス追加
ROOT.gSystem.AddDynamicPath("/opt/homebrew/Cellar/gsl/2.8/lib")
# GSLの動的リンクライブラリの追加　
# 標準的な拡張子(so,dll,dylib,..)なら拡張子は不要
ROOT.gSystem.Load("libgsl")
ROOT.gSystem.Load("libgslcblas")

# GSLのヘッダーファイルをinclude
ROOT.gInterpreter.ProcessLine("#include <gsl/gsl_poly.h>")

a = array.array('d',[-1]+[0]*4+[1])
z = array.array('d',[0]*10)

w = ROOT.gsl_poly_complex_workspace_alloc(6)

ROOT.gsl_poly_complex_solve (a, len(a), w, z);

ROOT.gsl_poly_complex_workspace_free (w)

for i in range(5):
    idx = i*2
    print("z{} = {:+.18f} {:+.18f}".format(i, z[idx], z[idx+1]))
print()
for i in range(5):
    idx = i*2
    x = complex(z[idx], z[idx+1])
    print("z{} = ".format(i), pow(x, 5))
```
**GSL**の配列の型は**double**になっているので、**array**か**numpy**を使ってデータタイプを**double**にして使います。
コードの最後で得られた解(x)を5乗して１になるか検証できるようにしています。


`gsl1.py`を実行します。

```terminal
> python gsl1.py
z0 = -0.809016994374947451 +0.587785252292473359
z1 = -0.809016994374947451 -0.587785252292473359
z2 = +0.309016994374947618 +0.951056516295154086
z3 = +0.309016994374947618 -0.951056516295154086
z4 = +1.000000000000000222 +0.000000000000000000

z0 =  (1.0000000000000009-7.686127107015008e-16j)
z1 =  (1.0000000000000009+7.686127107015008e-16j)
z2 =  (1.0000000000000029-7.65965843452121e-17j)
z3 =  (1.0000000000000029+7.65965843452121e-17j)
z4 =  (1.000000000000001+0j)
```

多少の誤差はありますが、ほぼ実数の1になっているのが確認できます。**j**は虚数単位です。

<details><summary>ドキュメントのソース&実行結果</summary>

```c
#include <stdio.h>
#include <gsl/gsl_poly.h>

int
main (void)
{
  int i;
  /* coefficients of P(x) =  -1 + x^5  */
  double a[6] = { -1, 0, 0, 0, 0, 1 };
  double z[10];

  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);

  for (i = 0; i < 5; i++)
    {
      printf ("z%d = %+.18f %+.18f\n",
              i, z[2*i], z[2*i+1]);
    }

  return 0;
}
```
```terminal
z0 = -0.809016994374947673 +0.587785252292473359
z1 = -0.809016994374947673 -0.587785252292473359
z2 = +0.309016994374947507 +0.951056516295152976
z3 = +0.309016994374947507 -0.951056516295152976
z4 = +0.999999999999999889 +0.000000000000000000
```


</details>


## myroot.py

PyROOTとGSLで毎回必要なコードを`myroot.py`にまとめて、それを**import**することにします。なお、概要でも述べたように**CERN ROOT**で**GSL**が使われているため、ライブラリパスやライブラリの指定をしなくても大丈夫のようなのでコメントにしました。

```py:myroot.py
import ROOT

ROOT.gInterpreter.AddIncludePath("/opt/homebrew/Cellar/gsl/2.8/include")
# ROOT.gSystem.AddDynamicPath("/opt/homebrew/Cellar/gsl/2.8/lib")
# ROOT.gSystem.Load("libgsl")
# ROOT.gSystem.Load("libgslcblas")

def include(inc_array):
	for inc in inc_array:
		ROOT.gInterpreter.ProcessLine("#include "+inc)

def declare(str):
	ROOT.gInterpreter.Declare(str)
```

上記`gsl1.py`は下記のようなコードに変わります。

```py:gsl1x.py
from myroot import *
import array

include(["<gsl/gsl_poly.h>"])

a = array.array('d',[-1]+[0]*4+[1])
z = array.array('d',[0]*10)

w = ROOT.gsl_poly_complex_workspace_alloc(6)

ROOT.gsl_poly_complex_solve (a, len(a), w, z);

ROOT.gsl_poly_complex_workspace_free (w)

for i in range(5):
    idx = i*2
    print("z{} = {:+.18f} {:+.18f}".format(i, z[idx], z[idx+1]))
print()
for i in range(5):
    idx = i*2
    x = complex(z[idx], z[idx+1])
    print("z{} = ".format(i), pow(x, 5))
```

**GSL**の関数を呼ぶのはC++のほうがPythonよりも書きやすいので、わざわざPythonで書く必要はなく入力と結果がPythonで使えれば良いとの考えで、下記のコードにしてみました。
`declare (ROOT.gInterpreter.Declare関数)`はC++の関数などを定義できます。**CERN ROOT**はC++のインタプリタを内蔵しています。

```py:gsl2.py
from myroot import *
import array

include(["<gsl/gsl_poly.h>"])

declare("""
void gsl2_func(int len, double *a, double *z) {
    auto w = gsl_poly_complex_workspace_alloc (len);
    gsl_poly_complex_solve (a, len, w, z);
    gsl_poly_complex_workspace_free (w);
}
""")

a = array.array('d',[-1]+[0]*4+[1])
z = array.array('d',[0]*10)
ROOT.gsl2_func(len(a), a, z)

for i in range(5):
    idx = i*2
    print("z{} = {:+.18f} {:+.18f}".format(i, z[idx], z[idx+1]))
print()
for i in range(5):
    idx = i*2
    x = complex(z[idx], z[idx+1])
    print("z{} = ".format(i), pow(x, 5))
```


## BLASで行列演算

上記のようにC++でコードを記述したほうが簡単ですが、これ以降の例はできるだけPythonで書けるところはPythonで書きます。

https://www.gnu.org/software/gsl/doc/html/blas.html#examples


行列の演算
$$
AB=C
$$

```gsl3.py
from myroot import *
import array

include(["<gsl/gsl_blas.h>"])

a = array.array('d', [0.11, 0.12, 0.13,
                      0.21, 0.22, 0.23])
b = array.array('d',[1011, 1012,
                     1021, 1022,
                     1031, 1032])
c = array.array('d', [0.00, 0.00,
                      0.00, 0.00])
A = ROOT.gsl_matrix_view_array(a, 2, 3)
B = ROOT.gsl_matrix_view_array(b, 3, 2)
C = ROOT.gsl_matrix_view_array(c, 2, 2)

#  /* Compute C = A B */
ROOT.gsl_blas_dgemm (ROOT.CblasNoTrans, ROOT.CblasNoTrans,
                  1.0, A.matrix, B.matrix,
                  0.0, C.matrix)
print ("[ {:6.2f}, {:6.2f}".format(c[0], c[1]))
print ("  {:6.2f}, {:6.2f} ]".format(c[2], c[3]))
```
`gsl3.py`の実行

```terminal
> python gsl3.py
[ 367.76, 368.12
  674.06, 674.72 ]
```
小数点２桁で表示しているためドキュメントとほぼ同じ結果になりましたが、四捨五入されているところもあります。

<details><summary>ドキュメントのソース&実行結果</summary>

```c
#include <stdio.h>
#include <gsl/gsl_blas.h>

int
main (void)
{
  double a[] = { 0.11, 0.12, 0.13,
                 0.21, 0.22, 0.23 };

  double b[] = { 1011, 1012,
                 1021, 1022,
                 1031, 1032 };

  double c[] = { 0.00, 0.00,
                 0.00, 0.00 };

  gsl_matrix_view A = gsl_matrix_view_array(a, 2, 3);
  gsl_matrix_view B = gsl_matrix_view_array(b, 3, 2);
  gsl_matrix_view C = gsl_matrix_view_array(c, 2, 2);

  /* Compute C = A B */

  gsl_blas_dgemm (CblasNoTrans, CblasNoTrans,
                  1.0, &A.matrix, &B.matrix,
                  0.0, &C.matrix);

  printf ("[ %g, %g\n", c[0], c[1]);
  printf ("  %g, %g ]\n", c[2], c[3]);

  return 0;
}
```
```terminal
[ 367.76, 368.12
  674.06, 674.72 ]
```

</details>

## 連立1次方程式の解

https://www.gnu.org/software/gsl/doc/html/linalg.html#examples

$$ A\boldsymbol x=\boldsymbol b$$

行列$A$が**a_data**、ベクトル$\boldsymbol b$が**b_data**でベクトル$\boldsymbol x$ を求める。

---

以下の`gsl4.py`でC言語を用いて書いているコードについて
- **stdout**がPythonで扱えなかった
    - **stdout**は `#define stdout __stdoutp`のように定義されていたが、`#define`で定義された識別子はPythonで扱えない
    - **ROOT.__stdoutp**としてはPythonからアクセスできたが、**ROOT.gsl_vector_fprintf**に渡そうとすると型が`FILE*`でないというエラーになった

```gsl4.py
from myroot import *
import array
import numpy as np

include(["<gsl/gsl_linalg.h>"])

declare("""
void vector_fprintf(gsl_vector *x) {
    gsl_vector_fprintf (stdout, x, "%g");
}
""")

a = [0.18, 0.60, 0.57, 0.96,
     0.41, 0.24, 0.99, 0.58,
     0.14, 0.30, 0.97, 0.66,
     0.51, 0.13, 0.19, 0.85]
a_data = array.array('d', a)

b_data = array.array('d', [1.0, 2.0, 3.0, 4.0])

m = ROOT.gsl_matrix_view_array (a_data, 4, 4)
b = ROOT.gsl_vector_view_array (b_data, 4)
x = ROOT.gsl_vector_alloc (4)
s = array.array('i',[0])
p = ROOT.gsl_permutation_alloc (4)
ROOT.gsl_linalg_LU_decomp (m.matrix, p, s)
ROOT.gsl_linalg_LU_solve (m.matrix, p, b.vector, x)

# Pythonでは実行できない
#ROOT.gsl_vector_fprintf (stdout, x, "%g")
ROOT.vector_fprintf(x)

# 求めたベクトルを取り出し、Pythonの配列に代入s
X = np.zeros(4, dtype=np.float64)
for i in range(4):
    X[i] = ROOT.gsl_vector_get(x, i)

ROOT.gsl_permutation_free (p)
ROOT.gsl_vector_free (x)

A = np.array(a)
A2 = A.reshape(4,4)  # 4x4の行列にする
X2 = X.reshape(4,1)  # 縦ベクトルにする
print("\n***********\n\n", A2.dot(X2))
```

`int s;`,`&s`のint型のポインタをPythonで直接に表すことはできなかった。そこで、配列の変数はポインタと同じ扱いになると思い`s`を配列とした。

`gsl4.py`を実行
```terminal
> python gsl4.py
-4.05205
-12.6056
1.66091
8.69377

***********

 [[1.]
 [2.]
 [3.]
 [4.]]
```
求めた解を**numpy**の行列演算を使って検証した結果が**b_data**と同じになり解の正しさが確認できた。`s`の配列定義もうまく行ったようである。

<details><summary>ドキュメントのソース&実行結果</summary>

```c
#include <stdio.h>
#include <gsl/gsl_linalg.h>

int
main (void)
{
  double a_data[] = { 0.18, 0.60, 0.57, 0.96,
                      0.41, 0.24, 0.99, 0.58,
                      0.14, 0.30, 0.97, 0.66,
                      0.51, 0.13, 0.19, 0.85 };

  double b_data[] = { 1.0, 2.0, 3.0, 4.0 };

  gsl_matrix_view m
    = gsl_matrix_view_array (a_data, 4, 4);

  gsl_vector_view b
    = gsl_vector_view_array (b_data, 4);

  gsl_vector *x = gsl_vector_alloc (4);

  int s;

  gsl_permutation * p = gsl_permutation_alloc (4);

  gsl_linalg_LU_decomp (&m.matrix, p, &s);

  gsl_linalg_LU_solve (&m.matrix, p, &b.vector, x);

  printf ("x = \n");
  gsl_vector_fprintf (stdout, x, "%g");

  gsl_permutation_free (p);
  gsl_vector_free (x);
  return 0;
}
```
```terminal
x =
-4.05205
-12.6056
1.66091
8.69377
```
</details>

## 物理定数

https://www.gnu.org/software/gsl/doc/html/const.html#examples


物理定数は`#define`で定義されているので、一度C++の変数に代入してからC++の変数をPythonで使うようにする。

```gsl5.py
from myroot import *

include(["<gsl/gsl_const_mksa.h>"])

declare("""
  constexpr double c  = GSL_CONST_MKSA_SPEED_OF_LIGHT;
// 太陽と地球の間の距離 (メートル)が1au
  constexpr double au = GSL_CONST_MKSA_ASTRONOMICAL_UNIT;
  constexpr double minutes = GSL_CONST_MKSA_MINUTE;
""")

c = ROOT.c
au = ROOT.au 
minutes = ROOT.minutes

#  /* distance stored in meters */
r_earth = 1.00 * au
r_mars  = 1.52 * au

# 火星と地球が最も接近したときの光が届く時間(秒)
t_min = (r_mars - r_earth) / c 
# 火星が太陽を挟んで地球の反対側にいるときの光が届く時間(秒)
t_max = (r_mars + r_earth) / c

print("light travel time from Earth to Mars:")
print("minimum = {:.1f} minutes".format(t_min / minutes))
print("maximum = {:.1f} minutes".format(t_max / minutes))
```
`gsl5.py`を実行
```terminal
> python gsl5.py
light travel time from Earth to Mars:
minimum = 4.3 minutes
maximum = 21.0 minutes
```

<details><summary>ドキュメントのソース&実行結果</summary>

```c
#include <stdio.h>
#include <gsl/gsl_const_mksa.h>

int
main (void)
{
  double c  = GSL_CONST_MKSA_SPEED_OF_LIGHT;
  double au = GSL_CONST_MKSA_ASTRONOMICAL_UNIT;
  double minutes = GSL_CONST_MKSA_MINUTE;

  /* distance stored in meters */
  double r_earth = 1.00 * au;
  double r_mars  = 1.52 * au;

  double t_min, t_max;

  t_min = (r_mars - r_earth) / c;
  t_max = (r_mars + r_earth) / c;

  printf ("light travel time from Earth to Mars:\n");
  printf ("minimum = %.1f minutes\n", t_min / minutes);
  printf ("maximum = %.1f minutes\n", t_max / minutes);

  return 0;
}
```
```terminal
light travel time from Earth to Mars:
minimum = 4.3 minutes
maximum = 21.0 minutes
```
</details>


## 数値積分

https://www.gnu.org/software/gsl/doc/html/integration.html#adaptive-integration-example

$$\int_0^1 x^{-1/2} \log(x)~dx=-4$$

数値積分がどの程度の精度で-4に近くなるかを確かめるプログラムです。

念のため、定積分を解いてみる。

$t=\sqrt{x}$とすると$x=t^2,~dx=2t~dt$なので
$$
\begin{align}
\int_0^1 x^{-1/2} \log(x)~dx &= \int_0^1 t^{-1} \log(t^2) 2t~dt \\\\ 
  &=4\int_0^1 \log (t)~dt \\\\
  &=4~[~t\log(t)-t~]_0^1 \\\\
  &=-4
\end{align}
$$

---

関数**f**内の間接参照はPythonでは出来ないのでC++の関数**convert_double**で実装しました。
また、関数**f**が何回呼ばれるか気になったのでカウントを追加しました。

```py:gsl6.py
from myroot import *
import array
import math

include(["<gsl/gsl_integration.h>"])

declare("""
double convert_double(void *p) {
    return *(double *)p;
}          
""")

called_count = 0 
def f(x, params):
    global called_count
    called_count += 1
    alpha = ROOT.convert_double(params)
    return math.log(alpha*x) / math.sqrt(x)

def make_array(a):
    return array.array('d',[a])

w  = ROOT.gsl_integration_workspace_alloc (1000);

result = make_array(0)
error = make_array(0)
expected = make_array(-4.0)
alpha = make_array(1.0)

F = ROOT.gsl_function()
F.function = f
F.params = alpha

ROOT.gsl_integration_qags (F, 0, 1, 0, 1e-7, 1000,
                        w, result, error)
  
print("result          = {:21.18f}".format(result[0]))
print("exact result    = {:21.18f}".format(expected[0]))
print("estimated error = {:21.18f}".format(error[0]))
print("actual error    = {:21.18f}".format(result[0] - expected[0]))
print("intervals       = {}".format(w.size))

ROOT.gsl_integration_workspace_free (w)

print("called_count    = {}".format(called_count))
```

`gsl6.py`を実行

```terminal
> python gsl6.py
result          = -3.999999999999982681
exact result    = -4.000000000000000000
estimated error =  0.000000000000233591
actual error    =  0.000000000000017319
intervals       = 8
called_count    = 315
```
**result**はドキュメントのexamplesとは異なりましたが、`actual error`の絶対値はこちらのほうが小さいです。

<details><summary>ドキュメントのソース&実行結果</summary>

```c
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>

double f (double x, void * params) {
  double alpha = *(double *) params;
  double f = log(alpha*x) / sqrt(x);
  return f;
}

int
main (void)
{
  gsl_integration_workspace * w
    = gsl_integration_workspace_alloc (1000);

  double result, error;
  double expected = -4.0;
  double alpha = 1.0;

  gsl_function F;
  F.function = &f;
  F.params = &alpha;

  gsl_integration_qags (&F, 0, 1, 0, 1e-7, 1000,
                        w, &result, &error);

  printf ("result          = % .18f\n", result);
  printf ("exact result    = % .18f\n", expected);
  printf ("estimated error = % .18f\n", error);
  printf ("actual error    = % .18f\n", result - expected);
  printf ("intervals       = %zu\n", w->size);

  gsl_integration_workspace_free (w);

  return 0;
}
```
```terminal
result          = -4.000000000000085265
exact result    = -4.000000000000000000
estimated error =  0.000000000000135447
actual error    = -0.000000000000085265
intervals       = 8
```

</details>


# 終わりに

一部はC/C++でないと表現出来ないのですが、以外と簡単にPythonで出来たと思います。
これは**GSL**を題材にPythonからPyROOTを使ってC言語のライブラリを使うときのやり方の一部を試しただけです。別のライブラリを使うときの参考にしてください。
今回やってみて一番悩んだのはint型,double型のポインタを渡すことで、配列を使うことに気づくには時間がかかりました。最初はC++でポインタ変数を作っていました。

# Appendix

最後に`数値積分`のC言語のコードの**main**を通常の関数と考えて、Pythonから入出力を**main**に与えて結果を得るようなコードを作って見ます。
**main**は特殊な名前で一般関数としては使えないため**mainx**とします。また、Pythonとのインターフェイスは`struct rslt`にしています。
<details><summary>コード</summary>

以下のCのソースのバイナリは作りません。つまり**gcc**等でコンパイルして動的リンクライブラリは作りません。

```c:gsl7.c
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>

struct rslt {
    double result;
    double error;
    double expected;
};
        

double f (double x, void * params) {
  double alpha = *(double *) params;
  double f = log(alpha*x) / sqrt(x);
  return f;
}
 
int
mainx (struct rslt& arg)
{
  gsl_integration_workspace * w
    = gsl_integration_workspace_alloc (1000);

  double result, error;
  double expected = -4.0;
  double alpha = 1.0;

  gsl_function F;
  F.function = &f;
  F.params = &alpha;

  gsl_integration_qags (&F, 0, 1, 0, 1e-7, 1000,
                        w, &result, &error);

  printf ("result          = % .18f\n", result);
  printf ("exact result    = % .18f\n", expected);
  printf ("estimated error = % .18f\n", error);
  printf ("actual error    = % .18f\n", result - expected);
  printf ("intervals       = %zu\n", w->size);

  gsl_integration_workspace_free (w);

  arg.result = result;
  arg.expected = expected;
  arg.error = error;
  return 0;
}
```

```py:gsl7.py
from myroot import *

include(['"./gsl7.c"'])

result = ROOT.rslt()
ROOT.mainx(result)
print()
print("result   = {:21.18f}".format(result.result))
print("error    = {:21.18f}".format(result.error))
print("expected = {:21.18f}".format(result.expected))
```
`gsl7.py`を実行

```terminal
> python gsl7.py
result          = -3.999999999999982681
exact result    = -4.000000000000000000
estimated error =  0.000000000000233591
actual error    =  0.000000000000017319
intervals       = 8

result   = -3.999999999999982681
error    =  0.000000000000233591
expected = -4.000000000000000000
```


</details>