C/C++ 之 GSL 数学运算库使用笔记

Part.I Introduction

本文主要记录一下笔者使用 GSL 过程当中所做的一些笔记。

Chap.I 传送门

一些传送门

• GSL源码（CMakeList 版本-Windows）
• GSL源码（configure 版本-Linux）
• GSL 在线文档
• GSL 文档下载

Chap.II GSL 简介

GSL 全称 GNU Scientific Library，即 GNU 科学计算库。GNU 是什么东西呢？GNU 是一个自由的操作系统，其内容软件完全以 GPL 方式发布。GPL 又是什么东西？不要急~要想理解GNU 操作系统，必须先了解 GUN Project ，GNU计划、GNU项目，随你怎么翻译好了…… 它又被译为 革奴计划 ，是一个自由软件集体协作计划，1983年9月27日由理查德·斯托曼在麻省理工学院公开发起。它的目的是创建一套完全自由的操作系统，称为GNU¹。GNU 全称 GNU’s Not Unix!，递归定义缩写。

GPL 的全称是 GNU General Public License ，被翻译为 GNU通用公共许可协议 ，缩写为 GPL 或 GNU GPL ，说到底，这只不过是一个协议系列，那这一系列协议有啥特点呢？它有两个特点（copyleft）：

• 一方面，它给予了用户充分的自由，允许用户 运行、学习、共享和修改软件 ；
• 而另一方面，它死死限制了用户的一个方面的自由，那就是：GPL的派生作品只能以相同的许可证 发布 。这两个特点结合起来翻译成人话就是“软件随便用，源码也给你，随你怎么copy怎么修改，这都是你的自由，但是！你不能将大家伙的劳动成果变成你一个人私有的！”

就简单将 GNU 理解为一个项目，在这个项目中需要遵守 GPL 协议，即可以运行、学习、共享和修改 GNU 项目中的软件，但是不能据为己有。而 GSL 就是 GNU 中的一个项目，它是一个数学相关的科学计算库。

Part.II GSL 编译与安装

因为笔者使用的是 Windows 系统，在 Windows 下可以使用 『MSYS2 + Visual Studio』 对 GSL 进行安装和编译，也可以使用『CMake + Visual Studio』 对 GSL 进行安装和编译。前者已经有前辈写好了教程，下面主要针对后者的操作进行一下记录。

Chap.I GSL 的编译

这一步的目的是生成 GSL 动态库。

首先安装 cmake，VS等，然后需要查看 cmake 的路径是否添加到了系统变量path中。方法：我的电脑→右键属性→高级系统设置→环境变量→系统变量→path→添加 cmake.exe 的路径进去

下载最新版 gsl-2.7。然后解压在gsl-master目录下新建build文件夹

生成动态库

cmake -G"Visual Studio 16 2019" -DGSL_INSTALL_MULTI_CONFIG=ON -DNO_AMPL_BINDINGS=1 -DCMAKE_INSTALL_PREFIX=%cd%/install -DBUILD_SHARED_LIBS=ON -DMSVC_RUNTIME_DYNAMIC=ON A:/downloads/gsl-master
// 注意最后的路径是你自己的路径哦

编译gsl

cmake --build . --config Release --target install
// 生成 Release 版本的库
cmake --build . --config Debug --target install
// 生成 Debug 版本的库

得到所需文件

最后得到的文件树如下：

gsl
├─bin
│  ├─Debug
│  └─Release
├─include
│  └─gsl
└─lib
    ├─cmake
    │  └─gsl-2.7
    ├─Debug
    ├─pkgconfig
    └─Release

好像 Release 版本和 Debug 版本得到的东西是一样的。下面将介绍如何使用得到的这些库。

Chap.II 基于 VS 可视化操作的使用

新建一个项目，然后依次加入『附加包含目录』、『附加库目录』、『附加依赖项』

// 附加包含目录
A:\aWork\third-party\gsl\include
// 附加库目录
A:\aWork\third-party\gsl\lib\Release
// 附加依赖项
gsl.lib
gslcblas.lib

注意：dll 文件要加到工程目录，或者将其路径加入环境变量。
注意：如果一个项目要依赖其他项目，可以 选中项目→右键『属性』→VC++ 目录→包含目录和库目录把 GSL 的路径添加进去。

测试代码如下：

#include <stdio.h>
#include <gsl/gsl_matrix.h>

int main(void)
{
    int i, j;
    gsl_matrix* m = gsl_matrix_alloc(3, 3);
    for (i = 0; i < 3; i++) {
        for (j = 0; j < 3; j++) {
            gsl_matrix_set(m, i, j, i + j);
        }
    }

    for (i = 0; i < 3; i++) {
        for (j = 0; j < 3; j++) {
            printf("m(%d,%d) = %g\n", i, j, gsl_matrix_get(m, i, j));
        }
    }

    gsl_matrix_free(m);
    return 0;
}

运行结果如下：

Chap.III 基于 CMake 的使用

基于 CMake 的使用其实和上面的差不多，就是“ 加入『附加包含目录』、『附加库目录』、『附加依赖项』”这几步操作变成了 CMake 语句，主要的语句为：

最好用 target_xx ，因为它只会影响当前 target，而另一个会影响所有的 target。

因为文件比较多，不方便都放上来，感兴趣的朋友可以戳我下载（要钱的哦），文件树如下：

.
 ─app_test
│  ├─TEMP_TestGSL
│  └─TEMP_TestPython
├─LibGSL
│  ├─example
│  └─gexport
├─LibPython
│  ├─example
│  │  └─__pycache__
│  └─gexport
└─_doc

Part.III GSL 内置分布

Chap.I 一些常量

GSL_POSINF		// ＋∞
GSL_NEGINF		// -∞
GSL_NAN			// NAN
GSL_POSZERO		// +0
GSL_NEGZERO		// -0
M_PI			// PI

Chap.I gamma 分布

需包含头文件

#include <gsl/gsl_sf.h>

• gsl_sf_gamma(a)，计算得到 $\Gamma (a)$，$a$ 要大于 0
• gsl_sf_gamma_inc(a,x)，计算得到不完全伽马函数 $\Gamma (a,x)$（其积分下限是 $x$）
• gsl_sf_gamma_inc_Q(a,x)，计算得到归一化的不完全伽马函数 $Q(a,x)$ （$\frac{\Gamma (a,x)}{\Gamma (x)}$）
• gsl_sf_gamma_inc_P(a,x)，计算得到互补归一化的不完全伽马函数 $P(a,x)=1-Q(a,x)$

Chap.II t 分布

需包含头文件

#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>

• gsl_ran_tdist_pdf(x,n)，计算得到自由度为 n，变量值为 x 的概率密度
• gsl_ran_tdist_pdf((x - mu) / sigma, n) / sigma，计算得到自由度为 n，变量值为 x，位置参数为 mu，标准差为 sigma 的
• gsl_cdf_tdist_P(x,n)，计算得到自由度为 n 的 t 函数从负无穷到 x 的积分值；或者说得到自由度为 n，变量值为 x 的累积分布函数值
• gsl_cdf_tdist_Q(x,n)，$1-P(x,n)$

关于位置尺度 t 分布
$$f(t)=\frac{\Gamma \left(\frac{n+1}{2}\right)}{\sigma \sqrt{\pi n}\Gamma \left(\frac{n}{2}\right)}{\left(1+\frac{{\left(\frac{t-\mu }{\sigma }\right)}^2}{n}\right)}^{-\frac{n+1}{2}}$$

# python 代码
print(t.pdf(2,3,0,2))		# x, n, mu, sigam
print(t.pdf(1,3)/2)			# x, n

可以用下面一行的代码，使得位置尺度 t 分布与通常意义上的 t 分布进行转换。

Chap.III 其他分布

• Laplace 分布：gsl_ran_laplace_pdf(double x, double a)

GSL 只定义了稳定分布随机数生成，没有定义它的 pdf

double gsl_ran_levy(const gsl_rng * r, const double c, const double alpha);
double gsl_ran_levy_skew(const gsl_rng * r, const double c, const double alpha, const double beta);

Part.IV 随机数的生成

生成随机需要引用头文件#include <gsl/gsl_rng.h>
重要函数

• gsl_rng_uniform：返回一个在[0,1)区间内的随机数
• gsl_rng_uniform_pos：返回一个在(0,1)区间内的随机数
• gsl_ran_exponential：返回一个服从指数分布 $\frac{1}{\mu }{e}^{-\frac{x}{\mu }},x>0$ 的随机数

Chap.I 报错的代码

下面的代码会报错，无法解析的外部符号 gsl_rng_default

#include <iostream>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>


int main()
{
    gsl_rng* r;
    int i, n = 10;
    double low = 1;
    double up = 10;
    const gsl_rng_type* T;

    gsl_rng_env_setup();
    T = gsl_rng_default;
    r = gsl_rng_alloc(T);

    for (i = 0; i < n; i++)
    {
        double u = low + gsl_rng_uniform(r) * (up - low);
        printf("%.5f\n", u);
    }

    gsl_rng_free(r);

    return 0;
}

Chap.II 可以跑通的

#include <iostream>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

using namespace std;

int main()
{
    gsl_rng* rng;
    const gsl_rng_type** r;
    const gsl_rng_type** rngs = gsl_rng_types_setup();

    double a = 1.85361, b = 0.612543, c = 3.78207, d = -0.749357;

    for (r = rngs; *r != 0; r++)
    {
        rng = gsl_rng_alloc(*r);
        double u = gsl_ran_levy_skew(rng, c, a, b);
        cout << u << "  " << rng->type->name << endl;
        gsl_rng_free(rng);
    }
    cout << "---------------------------" << endl;
    rng = gsl_rng_alloc(*rngs);
    for (int i = 0; i < 10; i++)
    {
        double u = gsl_ran_levy_skew(rng, c, a, b);
        cout << u << "  " << rng->type->name << endl;
    }
    gsl_rng_free(rng);
    getchar();
    return 0;
}

Part.V 求积分

求积分需要引用头文件#include <gsl/gsl_integration.h>，下面是一些重要的函数的参数

• gsl_function：输入，GSL 函数，需要构造
• epsabs：输入，绝对误差，一般设为0
• epsrel：输入，相对误差，一般设为1e-7
• limit：输入，子区间的最大数目
• gsl_integration_workspace：输入，积分寄存空间，子间隔及其结果存储在工作区提供的内存中
• result：输出，积分结果
• abserr：输出，绝对误差
• neval：不知道

Chap.I gsl_integration_qng

QNG算法是一种非自适应过程，它使用固定的Gauss-Kronrod-Patterson abscissae在最大87点处对被积函数进行采样。它提供了光滑函数的快速积分。函数原型如下：

int gsl_integration_qng(const gsl_function *f, double a, double b, double epsabs, double epsrel, double *result, double *abserr, size_t *neval)

示例：求积分 $G={\int }_{-1}^1\frac{1}{\sqrt{2\pi }}{e}^{-\frac{x^2}{2}}dx$ （标准正态分布）的值，已知 $G=0.6826$。

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

using namespace std;


double g1(double x, void* params)
{
    double f = 1 / sqrt(2 * M_PI) * exp(-x * x / 2);
    return f;
}

int main()
{
    double a = -1, b = 1;
    gsl_function gf;
    gf.function = g1;

    double r, er;
    size_t n;
    gsl_integration_qng(&gf, a, b, 1e-10, 1e-10, &r, &er, &n);
    cout << r << endl;
    return 0;
}

Chap.II gsl_integration_qag

QAG算法是一种简单的自适应积分算法。将积分区域划分为子区间，在每次迭代中对估计误差最大的子区间进行等分。这可以迅速降低总体误差，因为子区间集中在被积函数的局部困难周围。其函数原型如下：

int gsl_integration_qag(const gsl_function *f, double a, double b, double epsabs, double epsrel, size_t limit, int key, gsl_integration_workspace *workspace, double *result, double *abserr)

Chap.III gsl_integration_qags

当积分区域中存在可积奇点时，自适应程序将在奇点周围集中新的子区间。随着子区间的减小，积分的连续逼近以极限方式收敛。这种接近极限的方法可以使用外推过程来加速。QAGS算法将自适应二分法与Wynn - epsilon算法相结合，加快了对多种可积奇点的积分速度。其函数原型如下：

int gsl_integration_qags(const gsl_function *f, double a, double b, double epsabs, double epsrel, size_t limit, gsl_integration_workspace *workspace, double *result, double *abserr)

示例：求积分 $G={\int }_0^1{x}^{-\frac{1}{2}}log(\alpha x)dx$ 的值，已知 $\alpha =1$时，$G=-4$。

#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);
    getchar();
    return 0;
}

Chap.IV gsl_integration_qagp

其函数原型如下：

int gsl_integration_qagp(const gsl_function *f, double *pts, size_t npts, double epsabs, double epsrel, size_t limit, gsl_integration_workspace *workspace, double *result, double *abserr)

Chap.V gsl_integration_qagi

计算(-∞,+∞)的积分，其函数原型如下：

int gsl_integration_qagi(gsl_function *f, double epsabs, double epsrel, size_t limit, gsl_integration_workspace *workspace, double *result, double *abserr)

Chap.VI gsl_integration_qagiu

计算(a,+∞)的积分， 其函数原型如下：

int gsl_integration_qagiu(gsl_function *f, double a, double epsabs, double epsrel, size_t limit, gsl_integration_workspace *workspace, double *result, double *abserr)

Chap.VII gsl_integration_qagil

计算(-∞,b)的积分，其函数原型如下：

int gsl_integration_qagil(gsl_function *f, double b, double epsabs, double epsrel, size_t limit, gsl_integration_workspace *workspace, double *result, double *abserr)

Reference