概要

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は未確認

関連記事

実行環境

> 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次方程式の解

$$ x^5 = 1$$

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を実行します。

> 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は虚数単位です。

ドキュメントのソース&実行結果

#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;
}

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

myroot.py

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

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は下記のようなコードに変わります。

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++のインタプリタを内蔵しています。

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で書きます。

行列の演算
$$
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の実行

> python gsl3.py
[ 367.76, 368.12
  674.06, 674.72 ]

小数点２桁で表示しているためドキュメントとほぼ同じ結果になりましたが、四捨五入されているところもあります。

ドキュメントのソース&実行結果

#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;
}

[ 367.76, 368.12
  674.06, 674.72 ]

連立1次方程式の解

$$ 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を実行

> python gsl4.py
-4.05205
-12.6056
1.66091
8.69377

***********

 [[1.]
 [2.]
 [3.]
 [4.]]

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

ドキュメントのソース&実行結果

#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;
}

x =
-4.05205
-12.6056
1.66091
8.69377

物理定数

物理定数は#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を実行

> python gsl5.py
light travel time from Earth to Mars:
minimum = 4.3 minutes
maximum = 21.0 minutes

ドキュメントのソース&実行結果

#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;
}

light travel time from Earth to Mars:
minimum = 4.3 minutes
maximum = 21.0 minutes

数値積分

$$\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が何回呼ばれるか気になったのでカウントを追加しました。

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を実行

> 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の絶対値はこちらのほうが小さいです。

ドキュメントのソース&実行結果

#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;
}

result          = -4.000000000000085265
exact result    = -4.000000000000000000
estimated error =  0.000000000000135447
actual error    = -0.000000000000085265
intervals       = 8

終わりに

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

Appendix

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

コード

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

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

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を実行

> 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