GSL库的学习笔记

本文档整理翻译自GNU Scientific Library (GSL)的官方参考文件gsl-ref.pdf，对应的版本为released 2.7.
本文档License: (GNU Free Documentation License Version 1.3)

本文档的目的: 记录gsl库中的一些函数和定义，作为一个快速的笔记、查阅用的资料。简单的翻译。建议配合官方的gsl-ref.pdf使用。本文档难免有错误处，作者不对正确性作保证。

---

Ch 2 库的使用

2.10: 与C++的兼容性

当库函数被C++程序调用时，头文件会自动添加extern "C"的linkage, 故可以直接调用。

2.11：数组的aliasing问题

库默认所有传入的数组array、向量vector、矩阵matrices均没有aliasing问题、相互使用共同内存块的数据。否则， 某些函数表现将是undefined. 若该参数在函数型中声明为const则可以安全使用。

2.12：线程安全

库函数只要没有使用static的变量，一般就是线程安全thread-safe的。内存与对象关联而非与函数关联。对于有使用workspace的函数，需要为每一个线程独立定义；对于有使用table对象作为只读内存时, 可以安全地在不同线程中同步使用；Table作为参数时总是声明为const，因此（可能）可以多线程使用。static global变量应该先初始化且后面不再更改， 然后再多线程。

2.14 license

GNU GPL v>=3

Ch 3 错误的处理

3.1 错误error的报告方式

其一， 函数返回int型的status表征运行状态，0表示无错误，其他则表示对应错误代码errno的错误。(调用方通过返回status来作错误处理, 发布版本常用)

int status = gsl_function(...); 
if (status) { /* an error occured here，非0，即错误发生了*/
    /* now status equals the error code.此时status的值就是错误码*/
}

其二，错误处理器函数gsl_error(...)（Debug版本常用) 。某函数发生错误后、该函数返回之前被调用。默认的错误处理器函数的行为是: 输出一些错误信息（文件、错误行、错误描述等）然后退出程序abort。打印的内容例如:

gsl: file.c:67: ERROR: invalid argument supplied by user
Default GSL error handler invoked
Aborted

3.2 错误码

定义在gsl_errno.h中，库中使用<=1024的错误区间，>1024留给用户。-2~32翻看源码后发现是定义在enum中。 有GSL_EDOM GSL_ERANGE GSL_ENOMEM GSL_EINVAL等， 分别对应定义域错误、范围错误、内存不可得、非法的参数输入等。 把错误码转为字符串描述的函数: const char *gsl_strerror(const int gsl_errno), 使用例子:

printf("error: %s\n", gsl_strerror(status)); 

3.2 错误处理器

GSL所有错误处理器的类型都应该是gsl_error_handler_t。源代码中这样定义:

typedef void gsl_error_handler_t (const char * reason, const char * file,
                                  int line, int gsl_errno);

因此其实是函数类型，接受四个参数, reason为错误原因字符串，file为文件名，line行号，gsl_errno为错误代码。 file和line在编译阶段由__FILE__和__LINE__宏确定、用户不需传入。定义一个handler:

void myhandler(const char* reason, const char* file, int line, int gsl_errno){...} 

要求使用自定义的错误处理函数handler，并返回旧的handler,可以使用这个函数:

gsl_error_handler_t *gsl_set_error_handler(gsl_error_handler_t * new_handler) 

它接受一个新的handler函数指针，返回旧的处理函数的指针。 例子:

old_handler = gsl_set_error_handler(&myhandler); /* 保存旧的handler, 使用新的handler*/ 
... /*一些代码， 使用myhandler作为错误处理函数*/ 
gsl_set_error_handler(old_handler); /*恢复到旧的错误处理函数*/

当传入函数指针为NULL时，返回默认的handler的指针。

3.2 常用

GSL_ERROR(reason, gsl_errno)：宏，调用异常处理函数，报告异常。例如GSL_ERROR(“residual exceeds tolerance”, GSL_ETOL);

GSL_ERROR_VAL(reason, gsl_errno, value)：宏， 处理GSL_ERROR做的，还在当前函数return value;

gsl_error_handler_t *gsl_set_error_handler_off()关闭该功能，并返回旧的处理器函数指针。

3.n：一点细节

gsl_set_error_handler：GSL错误处理器函数的指针使用static修饰，这样全局只有一个错误处理器函数来接受和处理错误。 多线程程序里面需要注意。

Ch4 数学函数

4.1 常数

M_E	e	M_LOG2E	l o g 2 ( e ) log_2(e) log2​(e)
M_PI	π \pi π	M_LOG10E	l o g 10 ( e ) log_{10}(e) log10​(e)
M_PI_2	π / 2 \pi/2 π/2	M_SQRT2	2 \sqrt{2} 2 ​
M_PI_4	π / 4 \pi/4 π/4	M_SQRT1_2	1 / 2 1/\sqrt{2} 1/2 ​
M_SQRTPI	π \sqrt{\pi} π ​	M_SQRT3	3 \sqrt{3} 3 ​
M_1_PI	1 / π 1/\pi 1/π	M_LN10	ln ⁡ 10 \ln{10} ln10
M_2_PI	2 / π 2/\pi 2/π	M_LN2	ln ⁡ 2 \ln{2} ln2
M_EULER	Euler’s const., γ \gamma γ	M_LNPI	ln ⁡ π \ln{\pi} lnπ

4.2 Nan, Inf

GSL_POSINF: + ∞ +\infty +∞, +1.0/0.0
GSL_NEGINF: − ∞ -\infty −∞, -1.0/0.0
GSL_NAN: NaN,
int gsl_isnan(const doule x): 若nan则返回1， 否则返回0
int gsl_isinf(const double x)：若正无穷则返回+1，负无穷-1，有限则返回0
int gsl_finite(const double x)： 若x为实数返回1， 若nan或者无限则返回1

4.3 基础函数

（以下默认arguments: const double, 默认return: double）
gsl_log1p(x)： ln(1+x) x较小
gsl_exp1m(x): exp(x)-1 x较小
gsl_hypot(x, y): ( x 2 + y 2 ) \sqrt{(x^2 + y^2)} (x2+y2) ​
gsl_hypot3(x, y, z): ( x 2 + y 2 + z 2 ) \sqrt{(x^2 + y^2 + z^2)} (x2+y2+z2) ​
gsl_acosh(x) gsl_asinh(x) gsl_atanh(x) 对应的双曲函数
gsl_ldexp(double x, int e): x ∗ 2 e x*2^e x∗2e, LDEXP
gsl_frexp(double x, int *e): 把x分为f和e， x = f ∗ 2 e x = f*2^e x=f∗2e, 其中 0.5 ≤ f < 1 0.5 \le f<1 0.5≤f<1

4.4 小参数的幂

不检查overflow!
double gsl_pow_int(double x, int n) ： x n x^n xn，
gsl_pow_2(x), 3, 4, 5, …, 9 : 计算幂次，高效率。

4.5-4.7 一些帮助式的函数

GSL_SIGN(x)：宏， >0则为1，否则-1
GSL_IS_ODD(n)：宏， n为奇数则返回1， 否则0. n必须是整数类型
GSL_IS_EVEN(n)：宏， n为偶数则返回1， 否则0. n必须是整数类型
GSL_MAX(a, b)：宏， 返回最大的， 用>定义。 GSL_MIN类似

extern inline double GSL_MAX_DBL(double a, double b)： inline可用时，比宏更安全的方式，因为有类型检查。 同样还有GSL_MIN_DBL, GSL_MAX_INT GSL_MIN_INT, GSL_MAX_LDBL GSL_MIN_LDBL等的定义。

4.8 浮点数的比较

int gsl_fcmp(double x, double y, double epsilon)：x与y的差距在 2 δ , δ   e p s i l o n ∗ 2 k 2\delta, \delta~epsilon*2^k 2δ,δ epsilon∗2k，则返回1，否则为0. 注意，是相对精度的比较，相对误差大概小于epsilon/x就认为相等。

Ch5 复数

在头文件gsl_complex.h中定义了复数gsl_complex结构类型，一些函数和运算定义在gsl_complex_math.h中。

5.1 复数的表示

若>=C11, 且<complex.h>在gsl_complex.h之前导入则，typedef double complex gsl_complex
否则使用gsl_complex.h中的定义:

typedef struct {double dat[2];} gsl_complex; 

实部dat[0]，虚部dat[1].

gsl_complex x = 2 + 5 * I; 
gsl_complex y = x + (3 - 4 * I); 

5.2-5.3 宏与赋值

GSL_REAL(z) GSL_IMAG(z)： 宏，返回实步或虚部内存空间地址(左值)。 gsl_complex x; GSL_REAL(x) = 4; GSL_IMAG(x) = 2; GSL_REAL(y) = GSL_REAL(x); GSL_IMAG(y) = GSL_IMAG(x);

GSL_SET_COMPLEX(zp, x, y)：宏，对zp指向的复数赋值

gsl_complex gsl_complex_rect(double x, double y): z = x + y ∗ i z = x + y*i z=x+y∗i
gsl_complex gsl_complex_polar(double r, double theta): z = r ∗ e i θ z = r * e^{i\theta} z=r∗eiθ

5.4 复数的属性

double gsl_complex_arg(gsl_complex z)： a r g ( z ) , − π < a r g ( z ) < = π arg(z), -\pi < arg(z) <= \pi arg(z),−π<arg(z)<=π
double gsl_complex_abs(gsl_complex z): ∣ z ∣ |z| ∣z∣
double gsl_complex_abs2(gsl_complex z): ∣ z ∣ 2 |z|^2 ∣z∣2
double gsl_complex_logabs(gsl_complex z): l o g ( ∣ z ∣ ) log(|z|) log(∣z∣), 在 ∣ z ∣ ≃ 1 |z|\simeq1 ∣z∣≃1 附近计算有优化。而log(gsl_complex_abs(z))没有。

5.5 复数的运算

gsl_complex gsl_complex_add(gsl_complex a, gsl_complex b): z = a + b
以下默认return gsl_complex:
gsl_complex_sub(a, b): z = a-b
gsl_complex_mul(a, b): z = ab
gsl_complex_div(a, b): z = a/b,
gsl_complex_add_real(a, double x): z = a+x
gsl_complex_sub_real(a, double x): z = a-x
gsl_complex_mul_real(a, double x): z = ax
gsl_complex_div_real(a, double x): z = a/x
gsl_complex_sub_imag…
gsl_complex_conjudate(a): z = a ∗ = x − i y z = a^*= x-iy z=a∗=x−iy
gsl_complex_inverse(z): 1 / z = ( x − i y ) / ( x 2 + y 2 ) 1/z = (x-iy)/(x^2 + y^2) 1/z=(x−iy)/(x2+y2)
gsl_complex_negative(a): z = − a z = -a z=−a

5.6 函数

gsl_complex gsl_complex_sqrt(gsl_complex z): \sqrt{z}, cut: negative real axis, returned: real part >=0.
gsl_complex_sqrt_real(double x): z = x z = \sqrt{x} z=x ​, x可以小于0.
gsl_complex_pow(z, a)： z a = e x p ( l o g ( z ) ∗ a ) z^a = exp(log(z) * a) za=exp(log(z)∗a)
gsl_complex_pow_real(z, double x): z x z^x zx
gsl_complex_exp(z): e z e^z ez
gsl_complex_log(z): log(z), cut: negative real axis,
gsl_complex_log10(z): l o g 10 ( z ) log_{10}(z) log10​(z)
gsl_complex_log_b(z, gsl_complex b): l o g ( z ) / l o g ( b ) log(z)/log(b) log(z)/log(b)

5.7 复数的三角函数, 反三角函数，双曲函数, 反双曲函数

gsl_complex gsl_complex_sin(gsl_complex z): return s i n ( z ) = ( e x p ( i z ) − e x p ( − i z ) ) / ( 2 i ) sin(z) = (exp(iz)-exp(-iz))/(2i) sin(z)=(exp(iz)−exp(−iz))/(2i) 其余类似
cos, tan, sec, csc, cot
gsl_complex_arcsin(z): arcsin(z), cuts: z.real>1 and z.real<-1
gsl_complex_arccos(z): arccos(z), cuts: z.real>1 and z.real<-1
…
gsl_complex_sinh(z): s i n h ( z ) = ( e x p ( z ) − e x p ( − z ) ) / 2 sinh(z) = (exp(z)-exp(-z))/2 sinh(z)=(exp(z)−exp(−z))/2 , 类似地cosh, tanh, sech, csch, coth
…

Ch6 多项式

在头文件gsl_poly.h中。

6.1 计算多项式

P(x) = c[0] + c[1]x + c[2]x^2 +… + c[len-1] x^{len - 1}

使用Horner’s 方法求P(x)， 稳定。 若定义了HAVE_INLINE则下面几个函数使用inline版本

double gsl_poly_eval(const double c[], const int len, const double x): 返回P(x)
gsl_complex gsl_poly_complex_eval(const double c[], const int len, const gsl_complex z)： 实系数、复数z, 计算P(z).
gsl_complex gsl_complex_poly_complex_eval(const gsl_complex c[], const int len, const gsl_complex z)：复系数，复数z，计算P(z)

int gsl_poly_eval_derivs(const double c[], const size_t lenc, const double x, double res[], const size_t lenres)： 计算P(x)和各个阶次的导数，存到一个长度为lenres, 的res数组。依次存d^k/dxk, k=0, 1, …

6.3 二次方程的根

int gsl_poly_solve_quadratic(double a, double b, double c, double *x0, double *x1)： 判别式大于等于0，ax^2 + bx + c = 0的两个根存于x0, x1. 若a=0, 一个根存于x0. 若判别式小于0，没存。
int gsl_poly_complex_solve_quadratic(double x, double b, double c, gsl_complex *z0, gsl_complex *z1)

类似的， cubic： x 3 + a x 2 + b x + c = 0 x^3 + ax^2 + bx + c = 0 x3+ax2+bx+c=0的解。
gsl_poly_solve_cubic(a, b, c, double *z0, *z1, *z2)
gsl_poly_complex_solve_cubic(a, b, c, cpx *z0, *z1, *z2)

6.5 一般的多项式方程的求解

gsl_poly_complex_workspace类型: 寻根时的所需参数

gsl_poly_complex_workspace *gsl_poly_complex_workspace_alloc(size_t n): 申请内存，针对n系数多项式，与gsl_poly_complex_solve()函数适配。

void gsl_complex_workspace_free(… *w)： 空间释放

int gsl_poly_complex_solve(const double *a, size_t n, gsl_poly_complex_workspace w, gsl_complex_packed_ptr z)：计算根: P ( x ) = a 0 + a 1 x + a 2 ∗ x 2 + . . . + a n − 1 ∗ x n − 1 P(x) = a_0 + a_1x + a_2*x^2 + ... + a_{n-1} * x^{n-1} P(x)=a0​+a1​x+a2​∗x2+...+an−1​∗xn−1 ，方法是balanced-QR reduction of the companion matrix. a_(n-1) != 0, roots 由长度为2(n-1)的packed complex array 返回（double, [2i]为i-th根的实部，[2i+1]为虚部）。 returns: 所有根都找到GSL_SUCCESS, QR分解不收敛GSL_EFAILED.

6.5 例子

求解 P ( x ) = x 5 − 1 P(x) = x^5 -1 P(x)=x5−1 的根

#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]; /* size: 2 * (6 - 1) */
    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])
    }
}

Ch7 特殊函数

总gsl_sf.h： gsl_sf_airy.h, gsl_sf_bessel.h …

7.1 用法

一是返回结果型 double y = gsl_sf_bessel_J0(x); /* J 0 ( x ) J_0(x) J0​(x)*/

一是给结果指针指向的地方写入结果, 省掉一次数据拷贝。(但是很奇怪为嘛额外定义了一个gsl_sf_result来传递结果，未来与armadillo或eigen矩阵类整合的话还挺麻烦，，感觉不如直接double*传)

gsl_sf_result result;  
int status = gsl_sf_bessel_J0_e(x, &result); /* with _e */ 

而gsl_sf_result类型: struct {double val; double err; }

7.3 精度控制: modes

type gsl_model_t:
GSL_PREC_DOUBLE： 2E-16精度
GSL_PREC_SINGLE: 1E-7精度
GSL_PREC_APPRX: 5E-4精度

7.4 -7.34 具体函数略， 用时查阅

7.35 例子

#include <stdio.h>
#include <gsl/gsl_sf_bessel.h>
int
main (void)
{
double x = 5.0;
double expected = -0.17759677131433830434739701;
double y = gsl_sf_bessel_J0 (x);
printf ("J0(5.0) = %.18f\n", y);
printf ("exact = %.18f\n", expected);
return 0;
}
/*
 * J0(5.0) = -0.177596771314338264
 * exact = -0.177596771314338292
 */

Ch8 向量和矩阵

8.1 数据类型

Prefix	Type
gsl_block	double
gsl_block_float	float
gsl_block_long_double	long double
gsl_block_int, gsl_block_uint	int, unsigned int
gsl_block_long, gsl_block_ulong	long, unsigned long
gsl_block_short, gsl_block_ushort	short, unsigned short
…_char, …_uchar	char, unsigned char
gsl_block_complex	complex double
gsl_block_complex_float	complex float
gsl_block_complex_long_double	complex long double

8.2 blocks

type gsl_block: 所有的内存都是由gsl_block的结构申请的(定义在gsl_block.h)。向量和矩阵由block&slicing组成。

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

block的申请与释放:
gsl_block *gsl_block_alloc(size_t n): 申请空间，未初始化。
gsl_block *gsl_block_calloc(size_t n)：申请空间，初始化为0.
void gsl_block_free(gsl_block *b)： 释放b

#include <stdio.h>
#include <gsl/gsl_block.h>
int
main (void)
{
gsl_block * b = gsl_block_alloc (100);
printf ("length of block = %zu\n", b->size);
printf ("block data address = %p\n", b->data);
gsl_block_free (b);
return 0;
}

block也可由FILE* stream流等读写。

8.3 向量Vectors

type gsl_vector: 即typedef struct {size_t size; size_t stride; double * data; gsl_block * block; int owner;} gsl_vector; 其中size取0~size-1, stride物理内存上一个元素与下一个元素的内存距离，data: 指向内存中向量第一个元素的地址， block为块地址。owner可选，若vector拥有该block, 且owner设为1， 则该vector释放时将释放block；若指向了一个由另一个vector拥有的block, 则这个vector中owner为0，不释放block. (owner为1的vector负责释放block资源)。

gsl_vector * gsl_vector_alloc(size_t n): 创建vector，owner=1, 新创建一个block并申请大小为n的资源。
gsl_vector *gsl_vector_calloc(size_t n)： 附带初始化为0
void gsl_vector_free(gsl_vector *v): 判断owner的值并选择是否释放。

GSL_RANGE_CHECK_OFF: 宏，定义后，由原来的gsl_vector_get(v, i)替换成v->data[i*v->stride]，无边界检查。无性能损失了。
GSL_C99_INLINE: 若使用C99编译器，extern inline --> inline
int gsl_check_range：

double gsl_vector_get(const gsl_vector* v, const size_t i): 返回i-th 值，若越界则调用错误处理器。（使用inline版本，若定义了HAVE_INLINE）。
void gsl_vector_set(gsl_vector* v, size_t i)

double * gsl_vector_ptr(gsl_vector * v, size_t i)
double * gsl_vector_const_ptr(const gsl_vector * v, size_t i) ：返回指向i-th元素的指针，有边界检查(越界则处理并返回null)

初始化:

void gsl_vector_set_all(gsl_vector *v, double x)： 都设置为x
void gsl_vector_set_zero(gsl_vector *v)： 0
int gsl_vector_set_basis(gsl_vector *v, size_t i): 生成一个基矢vector， 非i-th均为0， i-th为1.

读写IO:

int gsl_vector_fwrite(FILE *stream, const gsl_vector *v): binary format写入。
int gsl_vector_fread(FILE *stream, gsl_vector *v): binary format读取。
int gsl_vector_fprintf(FILE *stream, const gsl_vector *v, const char *format): 一行一行输出，格式format为%g, %e, %f等等
int gsl_vector_fscanf(FILE *stream, gsl_vector *v): 读取

视图vector view:

type gsl_vector_view

type gsl_vector_const_view ：向量视图为临时性的对象，在栈stack上储存，用来操作向量元素的子集。视图下有vector来访问向量资源。最佳访问方法是: 始终使用指针访问vector以避免制造新的vector临时变量。&view.vector即为gsl_vector* 或const gsl_vector*类型的指针。

gsl_vector_view gsl_vector_subvector(gsl_vector *v, size_t offset, size_t n)

gsl_vector_const_view gsl_vector_const_subvector(const gsl_vector *v, size_t offset, size_t n) ：返回向量v的一个视图，偏移offset，总长度n. 视图指向向量的i-th元素为旧向量的(offset + i)-th元素， i从0到n-1。n太大则越界部分指针为null. 视图中的数据内存上与v相同。向量v的生命周期不可以小于视图。

gsl_vector_view gsl_vector_subvector_with_stride(gsl_vector* v, size_t offset, size_t stride, size_t n)：

gsl_vector_const_view gsl_vector_const_subvector_with_stride(const gsl_vector* v, size_t offset, size_t stride, size_t n)： v[offset:stride:…]步长stride，共n个元素。i-th --> (offset + i*stride)-th . 例如:

gsl_vector_view v_even = gsl_vector_subvector_with_stride(v, 0, 2, n/2);
gsl_vector_set_zero( &v_even.vector ); 

函数调用中&v_even.vector与gsl_vector* 类型一致。

gsl_vector_view gsl_vector_complex_real(gsl_vector_complex *v)

gsl_vector_const_view gsl_vector_complex_const_real(const gsl_vector_complex *v): v的实部的视图， 同理_imag…

gsl_vector_view gsl_vector_view_array(double *base, size_t n)

gsl_vector_const_view gsl_vector_const_view_array(const double *base, size_t n): C数组的向量视图，i-th（base[i]）–> v(i)

gsl_vector_view gsl_vector_view_array_with_stride(double *base, size_t stride, size_t n)

gsl_vector_const_view gsl_vector_const_view_array_with_stride(const double *base, size_t stride, size_t n): C数组的视图，带步长。

复制向量

向量的加减乘运算: 可以借助本库对BLAS的支持(见后文)。一下非BLAS函数:

int gsl_vector_memcpy(gsl_vector *dest, const gsl_vector *src): 由src复制到dest, 注意两者必须由相同的长度！

int gsl_vector_complex_conj_memcpy(gsl_vector_complex* dest, const gsl_vector_complex * src): 将src的共轭复制到dest中。两者必须有相同的长度。

int gsl_vector_swap(gsl_vector *v, gsl_vector *w)： 通过复制，交换v和w的元素，两者必须有相同的长度。

int gsl_vector_swap_elements(gsl_vector* v, size_t i, size_t j): 原位置交换i-th j-th的元素

int gsl_vector_reverse(gsl_vector *v)： 反向排列

向量的操作operations

（arguments: gsl_vector* -> vec * ; const gsl_vector* -> cvec *）

int gsl_vector_add(vec* a, cvec *b): a = a + b, 遍历元素，a原址操作, （a b必须等长度）

int gsl_vector_sub(vec *a, cvec *b): a = a - b, 遍历元素，a原址操作 （a b必须等长度）

int gsl_vector_mul(vec *a, cvec *b): a = a * b, 遍历元素，a原址操作 （a b必须等长度）

int gsl_vector_div(vec *a, cvec *b): a = a / b, 遍历元素，a原址操作 （a b必须等长度）

int gsl_vector_scale(vec *a, const double x): a = a * x, 遍历元素，a原址操作

int gsl_vector_add_constatn(vec* a, const double x): a = a + x, 遍历元素，a原址操作

double gsl_vector_sum(cvec *a): a所有元素求和

int gsl_vector_axpby(const double alpha, cvec* x, const double beta, vec *y)： y = alpha * x + beta * y, 遍历元素，y原址操作

向量最值

double gsl_vector_max(cvec * v): v的最大值。 gsl_vector_min(cvec * v)： v的最小值

void gsl_vector_minmax(cvec* v, double *min_out, double *max_out): 最大和最小

size_t gsl_vector_max_index(cvec* v)：最大值下标。 同样gsl_vector_min_index

void gsl_vector_minmax_index(cvec* v, size_t * imin, size_t * imax)： 最大最小的下标。

向量判断

int gsl_vector_isnull(cvec* v), _ispos, _isneg, _isnoneg等：若所有元素都满足该条件则返回1.

int gsl_vector_equal(cvec* u, cvec* v)：若所有元素equal则返回1.

向量举例

初始化、元素访问:

#include <stdio.h>
#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++) /* OUT OF RANGE ERROR */
    {
        printf ("v_%d = %g\n", i, gsl_vector_get (v, i));
    }
    gsl_vector_free (v);
    return 0;
}

向量写入文件:

#include <stdio.h>
#include <gsl/gsl_vector.h>
int
main (void)
{
    int i;
    gsl_vector * v = gsl_vector_alloc (100);
    for (i = 0; i < 100; i++)
    {
        gsl_vector_set (v, i, 1.23 + i);
    }
    
    {
        FILE * f = fopen ("test.dat", "w");
        gsl_vector_fprintf (f, v, "%.5g");
        fclose (f);
    }
gsl_vector_free (v);
return 0;
}

8.4 矩阵gsl_matrix

gsl_matrix.h 中

type gsl_matrix ：该结构体含有6个组分，2个矩阵维度大小，内存间距，矩阵存储位置指针data，gsl_block指针，owner.

注意， gsl_matrix以行优先(row_major)来储存。rows: 0~size1 - 1。 cols: 0~size2 - 1。tda: 矩阵中一行的内存大小！

typedef struct {size_t size1; size_t size2; size_t tda; 
               double * data; 
               gsl_block * block;
               int owner;
               } gsl_matrix; 

矩阵申请与释放

gsl_matrix * gsl_matrix_alloc(size_t n1, size_t n2)： n1行n2列（n1-by-n2）的矩阵。

gsl_matrix * gsl_matrix_calloc(size_t n1, size_t n2)：n1-by-n2, 附带全部初始化为0

void gsl_matrix_free(gsl_matrix * m)： 释放(若owner==1)

矩阵数据访问

GSL_RANGE_CHECK_OFF： 若定义，忽略边界检查. gsl_matrix_get(m, i, j), gsl_matrix_set(m, i, j, x)会使用m->data[i * m->tda + j].

double gsl_matrix_get(const gsl_matrix *m, const size_t i, const size_t j)：(i, j)个元素。若HAVE_INLINE则使用inline版本。

void gsl_matrix_set(gsl_matrix *m, const size_t i, const size_t j, double x)：(i, j)元素设为x

double *gsl_matrix_ptr(gsl_matrix *m, size_t i, size_t j)

const double *gsl_matrix_const_ptr(const gsl_matrix *m, size_t i, size_t j) 元素的指针

初始化

void gsl_matrix_set_all(gsl_matrix *m, double x): 均设为x

void gsl_matrix_set_zero(gsl_matrix *m)： 0

void gsl_matrix_set_identity(gsl_matrix *m): i=j的设为1， 其他为0.

矩阵读写

可以从文件流中读写

int gsl_matrix_fwrite(FILE *stream, const gsl_matrix *m)

int gsl_matrix_fread(FILE *stream, gsl_matrix *m)

int gsl_matrix_fprintf(FILE *stream, const gsl_matrix *m, const char *format)

int gsl_matrix_fscanf(FILE *stream, gsl_matrix *m)

矩阵视图

type gsl_matrix_view

type gsl_matrix_const_view：注意用&view.matrix表示视图下的矩阵对象。视图不拥有矩阵内存。

gsl_matrix_view gsl_matrix_submatrix(gsl_matrix *m, size_t k1, size_t k2, size_t n1, size_t n2)

gsl_matrix_const_view gsl_matrix_const_submatrix(const gsl_matrix m, size_t k1, size_t k2, size_t n1, size_t
n2): 视图中子矩阵n1行n2列，左上角的元素是m的(k1, k2)元素。 子矩阵的(i,j) --> m->data[(k1m->tda + k2) + i*m->tda + j]

(以下函数基本都有const版本， 我们省略它。)

gsl_matrix_view gsl_matrix_view_array(double *base, size_t n1, size_t n2)：一维C数组base的矩阵视图， n1-by-n2。 (i,j ) ==> base[i * n2 + j]

gsl_matrix_view gsl_matrix_view_array_with_tda(double *base, size_t n1, size_t n2, size_t tda)： n1-by-n2, (i,j ) ==> base[i * tda + j]

gsl_matrix_view gsl_matrix_view_vector(gsl_vector *v, size_t n1, size_t n2)： 向量v的矩阵视图。n1-by-n2。(i, j) ==> v->data[i * n2 + j]

gsl_matrix_view gsl_matrix_view_vector_with_tda(gsl_vector *v, size_t n1, size_t n2, size_t tda)： 自定义tda的版本

矩阵的行视图，列视图

gsl_vector_view gsl_matrix_row(gsl_matrix *m, size_t i): 向量视图， m的第i行。

gsl_vector_view gsl_matrix_column(gsl_matrix *m, size_t j): 向量视图， 矩阵的第j列。

gsl_vector_view gsl_matrix_subrow(gsl_matrix *m, size_t i, size_t offset, size_t n): 向量视图， 矩阵的第i行，从offset处开始n个元素

gsl_vector_const_view gsl_matrix_const_subcolumn(const gsl_matrix *m, size_t j, size_t offset, size_t n): 向量视图， 矩阵的第j列，从offset处开始n个元素

gsl_vector_view gsl_matrix_diagonal(gsl_matrix *m): i=j的元素组成的视图， m不一定方阵

gsl_vector_view gsl_matrix_subdiagonal(gsl_matrix *m, size_t k): 第k个下对角线， m不一定方阵

gsl_vector_view gsl_matrix_superdiagonal(gsl_matrix *m, size_t k): 第k个上对角线， m不一定方阵

矩阵的复制

int gsl_matrix_memcpy(gsl_matrix *dest, const gsl_matrix *src): 将src复制到dest

int gsl_matrix_swap(gsl_matrix *m1, gsl_matrix *m2)： 通过复制，交换m1 m2元素

矩阵的行 列复制

int gsl_matrix_get_row(gsl_vector *v, const gsl_matrix *m, size_t i)：将m的第i行，复制到v中。两者长度应该一致

int gsl_matrix_get_col(gsl_vector *v, const gsl_matrix *m, size_t j)：将m的第j列， 复制到v中。两者长度应该一致。

int gsl_matrix_set_row(gsl_matrix *m, size_t i, const gsl_vector *v)：用v设置m中第i行的值。

int gsl_matrix_set_col(gsl_matrix *m, size_t j, const gsl_vector *v)：用v设置m中第j列的值。

交换矩阵的行或列

int gsl_matrix_swap_rows(gsl_matrix *m, size_t i, size_t j): m的i和j行互换

int gsl_matrix_swap_columns(gsl_matrix *m, size_t i, size_t j)：m的i和j列互换

int gsl_matrix_swap_rowcol(gsl_matrix *m, size_t i, size_t j): 第i行与第j列的数据互换，必须为方阵。

int gsl_matrix_transpose_memcpy(gsl_matrix *dest, const gsl_matrix *src): src矩阵转置，结果保存到dest中。

int gsl_matrix_transpose(gsl_matrix *m)：原位置转置。 必须为方阵！

int gsl_matrix_complex_conjtrans_memcpy(gsl_matrix_complex *dest, const gsl_matrix_complex *src): 共轭转置。

矩阵操作

int gsl_matrix_add(gsl_matrix *a, const gsl_matrix *b): a = a + b, 遍历元素，a原址操作。 同样还有_sub, _mul_elements, _div_elements, _add_constant

int gsl_matrix_scale(gsl_matrix *a, const double x)：个元素乘以x

int gsl_matrix_scale_columns(gsl_matrix *A, const gsl_vector *x)：A的第j列乘以x[j], 遍历所有列。

int gsl_matrix_scale_rows(gsl_matrix *A, const gsl_vector *x)：A的第i行乘以x[i], 遍历所有行。

矩阵最值

double gsl_matrix_max(const gsl_matrix *m): 最大。 类似的gsl_matrix_min最小， gsl_matrix_max_index最大值的下标，gsl_matrix_min_index最小值的下标，

void gsl_matrix_minmax_index(const gsl_matrix *m, size_t *imin, size_t *jmin, size_t *imax, size_t *jmax): 最大与最小值的下标。

矩阵属性

int gsl_matrix_isnull(const gsl_matrix *m)：类似的，_ispos, _isneg, isnoneg, 所有值均满足是时返回1.

int gsl_matrix_equal(const gsl_matrix *a, const gsl_matrix *b): 各个元素是否都相等。

double gsl_matrix_norm1(const gsl_matrix *A)：矩阵A的1-norm，||A||_1 = max _{1<=j<=n} {sum {i=1, m} {|A_ij|} }即最大的"列绝对值和"。

Ch9 排列permutation

单词permutation在这里怎么翻译好呢QAQ. 它的意思就是，n个整数的序列(0， 1， 2， …， n-1)，其permulation后得到新的序列中，每个元素仍然只出现一个，元素出现位置任意，大概n!种可能。就类似"排列组合"一词中的排列。

在gsl_permutation.h

9.1 排列的结构体

type gsl_permutation: 含有两个元素的结构体，大小size和permutation数组的指针data

typedef struct {size_t size; size_t * data;} gsl_permutation; 

9.2 排列结构的内存申请

gsl_permutation *gsl_permutation_alloc(size_t n)： 申请大小为n的排列，未初始化

gsl_permutation *gsl_permutation_calloc(size_t n)： 申请大小为n的排列，初始化为恒等排列。 (1, 2, 3) 经过恒等排列后仍为其自身。

void gsl_permutation_init(gsl_permutation *p): 初始化为恒等排列。

void gsl_permutation_free(gsl_permutation *p): 释放p

int gsl_permutation_memcpy(gsl_permutation *dest, const gsl_permutation *src)：把src复制到dest

9.3 访问元素

size_t gsl_permutation_get(const gsl_permutation *p, const size_t i)：返回p的第i个元素，若i超出0~n-1的范围，则调用错误处理器并返回0。

int gsl_permutation_swap(gsl_permutation *p, const size_t i, const size_t j): 交换p的第i和第j个元素

9.4 排列属性

size_t gsl_permutation_size(const gsl_permutation *p): 返回p的size

size_t *gsl_permutation_data(const gsl_permutation *p)：返回指向p中元素数组的指针。

int gsl_permutation_valid(const gsl_permutation *p): 检查p的元素: 0~n-1的整数是否都只出现一次。

9.5 排列相关函数

void gsl_permutation_reverse(gsl_permutation *p): 反序

int gsl_permutation_inverse(gsl_permutation *inv, const gsl_permutation *p): 将p的逆， 赋值到inv中。

int gsl_permutation_next(gsl_permutation *p): 排列p的前一步的排列。若不存在后一步，则返回GSL_FAILURE，p不变。

int gsl_permutation_prev(gsl_permutation *p): 退到排列p的前一步的 排列,（排列的大小: 按照lexicographic字典式顺序）。我的理解是(0, 1, 3, 2) --> (0, 1, 2, 3) 这样？ 若不存在前一步，则返回GSL_FAILURE，p不变。

9.6 应用排列

在gsl_permute.h和gsl_permute_vector.h中。

int gsl_permute(const size_t *p, double *data, size_t stride, size_t n): 将p应用到data, data的stride为stride

int gsl_permute_inverse(const size_t *p, double *data, size_t stride, size_t n)：应用p的逆到数组data中。

int gsl_permute_vector(const gsl_permutation *p, gsl_vector *v): 将排列p应用到向量v, v为行向量。相当于out = vP，P的第j列是单位阵的第p->data[j]列。

int gsl_permute_vector_inverse(const gsl_permutation *p, gsl_vector *v): 相当于out = vP^{T}, P的定义同上

int gsl_permute_matrix(const gsl_permutation *p, gsl_matrix *A): 相当于Aout = AP, P的定义同上

int gsl_permutation_mul(gsl_permutation *p, const gsl_permutation *pa, const gsl_permutation *pb): 将pa和pb操作合并为新的排列p, p = pa * pb, p相当于先应用pb, 然后应用pa.

9.7 排列的读写

int gsl_permutation_fwrite(FILE *stream, const gsl_permutation *p)

int gsl_permutation_fread(FILE *stream, gsl_permutation *p)

int gsl_permutation_fprintf(FILE *stream, const gsl_permutation *p, const char *format)

int gsl_permutation_fscanf(FILE *stream, gsl_permutation *p)

9.8 排列的环表示 略

9.9 例子

#include <stdio.h>
#include <gsl/gsl_permutation.h>
int
main (void)
{
    gsl_permutation * p = gsl_permutation_alloc (3);
    gsl_permutation_init (p);
    do{
        gsl_permutation_fprintf (stdout, p, " %u");
        printf ("\n");
    }
    while (gsl_permutation_next(p) == GSL_SUCCESS);
    gsl_permutation_free (p);
    return 0;
}
/*
0 1 2
0 2 1
1 0 2
1 2 0
2 0 1
2 1 0
*/

ch10 组合

组合: n个元素， 选择k个， k个中不考虑顺序 （ C n k C_n^k Cnk​种可能） gsl_combination.h

10.1 组合结构体

type gsl_combination：含有三个组分的结构体， n, k, data数组指针

typedef struct {size_t n; size_t k; size_t *data;} gsl_combination; 

10.2 组合的申请与释放

gsl_combination *gsl_combination_alloc(size_t n, size_t k): 申请且未初始化

gsl_combination *gsl_combination_calloc(size_t n, size_t k): 申请，且初始化为(字典顺序)第一个组合。

void gsl_combination_init_first(gsl_combination *c): 初始化为(字典顺序)第一个组合。 (0, 1, 2, …, k-1)

void gsl_combination_init_last(gsl_combination *c)

void gsl_combination_free(gsl_combination *c): 释放

int gsl_combination_memcpy(gsl_combination *dest, const gsl_combination *src) : 复制

10.3 访问元素

size_t gsl_combination_get(const gsl_combination *c, const size_t i): 获取第i个元素

10.4 组合的性质

size_t gsl_combination_n(const gsl_combination *c): 返回n

size_t gsl_combination_k(const gsl_combination *c)： 返回k

size_t *gsl_combination_data(const gsl_combination *c)： 返回指向data[0]的指针

int gsl_combination_valid(gsl_combination *c): 检查是否合法

10.5 组合相关的函数

int gsl_combination_next(gsl_combination *c): 原地取下一个的组合，字典顺序。 没有则返回GSL_FAILURE, c不变。

int gsl_combination_prev(gsl_combination *c): 上一个。 没有则返回GSL_FAILURE, c不变。

10.6 组合的读写

int gsl_combination_fwrite(FILE *stream, const gsl_combination *c)

int gsl_combination_fread(FILE *stream, gsl_combination *c)

int gsl_combination_fprintf(FILE *stream, const gsl_combination *c, const char *format): format：推荐"%zu\n" %z对应size_t

int gsl_combination_fscanf(FILE *stream, gsl_combination *c)

10.7 组合的例子

#include <stdio.h>
#include <gsl/gsl_combination.h>
int
main (void)
{
    gsl_combination * c;
    size_t i;
    printf ("All subsets of {0,1,2,3} by size:\n") ;
    for (i = 0; i <= 4; i++)
    {
        c = gsl_combination_calloc (4, i);
        do
        {
        printf ("{");
        gsl_combination_fprintf (stdout, c, " %u");
        printf (" }\n");
        }
        while (gsl_combination_next (c) == GSL_SUCCESS);
        gsl_combination_free (c);
    }
    return 0;
}
/*
All subsets of {0,1,2,3} by size:
{ }
{ 0 }
{ 1 }
{ 2 }
{ 3 }
{ 0 1 }
{ 0 2 }
{ 0 3 }
{ 1 2 }
{ 1 3 }
{ 2 3 }
{ 0 1 2 }
{ 0 1 3 }
{ 0 2 3 }
{ 1 2 3 }
{ 0 1 2 3 }
*/

Ch11 multiset

类似组合，只不过每个元素可以出现一次或多次。 gsl_multiset.h，type gsl_multiset， 其他与组合基本一致。把gsl_combination替换为gsl_multiset即可。

Ch12 排序

使用heapsort堆排序方法，O(N log(N) )。 注意不能保持相等元素的相对顺序–>不稳定排序。

12.1 排序对象

gsl_heapsort.h中。如果有标准的qsort()函数，建议使用该函数而不是这里的。

void gsl_heapsort(void *array, size_t count, size_t size, gsl_comparison_fn_t compare): 对array的count个元素升序排列，每个元素大小size, 比较函数为compare, 类型定义如下

type gsl_comparison_fn_t ：函数指针类型. 函数接受两个const void* 指针a, b, 若*a > *b 返回1， 等返回0， 小于则返回-1. 比如:

int (*gsl_comparison_fn_t) (const void *a, const void *b)

int
compare_doubles (const double * a, const double * b)
{
    if (*a > *b) return 1;
    else if (*a < *b) return -1;
    else return 0;
}
... 
gsl_heapsort (array, count, sizeof(double), compare_doubles);

int gsl_heapsort_index(size_t *p, const void *array, size_t count, size_t size, gsl_comparison_fn_t compare)： p储存该位置应该存储的元素的index.

12.2 向量排序

gsl_sort.h gsl_sort_vector.h: gsl_sort_float.h gsl_sort_vector_float.h ......

void gsl_sort(double *data, const size_t stride, size_t n): 排列data

void gsl_sort2(double *data1, const size_t stride1, double *data2, const size_t stride2, size_t n)：数值顺序排列data1, 同时按照相同的安排来排列data2.

void gsl_sort_vector(gsl_vector *v)： 原位排列v

void gsl_sort_vector2(gsl_vector *v1, gsl_vector *v2): 原位排列v1, 同时v2按照v1的安排原位置排列。

void gsl_sort_index(size_t *p, const double *data, size_t stride, size_t n): 将拍好data的下标安排存到p中。

int gsl_sort_vector_index(gsl_permutation *p, const gsl_vector *v): 将排v的安排存到p, 这样p的第一个元素对应于v的最大位置，…

12.3 前k个最大或最小

O(k*N) , 适应于k<~log(N)的情况。 当N=10000，k= 5000时，不如直接排序。

int gsl_sort_smallest(double *dest, size_t k, const double *src, size_t stride, size_t n): k个最大复制到dest中。

int gsl_sort_largest(double *dest, size_t k, const double *src, size_t stride, size_t n)

int gsl_sort_vector_smallest(double *dest, size_t k, const gsl_vector *v)

int gsl_sort_vector_largest(double *dest, size_t k, const gsl_vector *v)

int gsl_sort_smallest_index(size_t *p, size_t k, const double *src, size_t stride, size_t n)

int gsl_sort_largest_index(size_t *p, size_t k, const double *src, size_t stride, size_t n)

int gsl_sort_vector_smallest_index(size_t *p, size_t k, const gsl_vector *v)

int gsl_sort_vector_largest_index(size_t *p, size_t k, const gsl_vector *v)

Ch13 BLAS支持

本库支持: 1. 低级层，直接与C语言BLAS标准相联系 (CBLAS); 2. 高级层，针对GSL的向量和矩阵的操作。不包括稀疏矩阵相关操作。

gsl_cblas.h：gsl_cblas层的接口， BLAS Technical Forum’s standard for the C interface to legacy BLAS implementations，用户也可以选择CBLAS的其他实现。

BLAS的所有低级的函数表，见Ch41.

13.0 BLAS简单规范

三层结构: Level1 向量的运算， 如y = a*x + y. Level 2 矩阵-向量的运算， 如y = a * A * x + b * y. Level 3 矩阵-矩阵运算

运算的名字: DOT点积 x^T*y; AXPY 向量和， ax+y; MV 矩阵matrix-向量vector的乘积 Ax; SV 矩阵-向量求解 inv(A)x; MM 矩阵矩阵积； SM 矩阵矩阵求解inv(A)B

矩阵的类型: GE general, 普通。GB general band。SY 对称阵。SB 对称band。SP对称packed。HE 厄密。HB厄密band。HP厄密packed。TR三角。TB三角band。TP三角packed。

每个操作都有四种准确度: S 单精度。D 双精度。C 单精度复数型。Z 双精度复数型。

比如 SGEMM --> 单精度、普通矩阵、矩阵乘积运算。

13.1 GSL BLAS 接口

在gsl_blas.h中

Level 1

int gsl_blas_sdsdot(float alpha, const gsl_vector_float *x, const gsl_vector_float *y, float *result)：求和 alpha + xT * y，结果存到result, 单精度

int gsl_blas_sdot(const gsl_vector_float *x, const gsl_vector_float *y, float *result)

int gsl_blas_dsdot(const gsl_vector_float *x, const gsl_vector_float *y, double *result)

int gsl_blas_ddot(const gsl_vector *x, const gsl_vector *y, double *result): 计算xT * y, 存到result中。

int gsl_blas_cdotu(const gsl_vector_complex_float *x, const gsl_vector_complex_float *y, gsl_complex_float
*dotu)

int gsl_blas_zdotu(const gsl_vector_complex *x, const gsl_vector_complex *y, gsl_complex *dotu): 计算xT * y, 存到dotu中。

int gsl_blas_cdotc(const gsl_vector_complex_float *x, const gsl_vector_complex_float *y, gsl_complex_float
*dotc)

int gsl_blas_zdotc(const gsl_vector_complex *x, const gsl_vector_complex *y, gsl_complex *dotc): 计算xH * y, 存到dotc。

float gsl_blas_snrm2(const gsl_vector_float *x)

double gsl_blas_dnrm2(const gsl_vector *x): 计算Euclidean norm ∣ ∣ x ∣ ∣ 2 = ∑ i x i 2 ||x||_2 = \sqrt{\sum_i x_i^2} ∣∣x∣∣2​=∑i​xi2​ ​。类似向量的长度

float gsl_blas_scnrm2(const gsl_vector_complex_float *x)

double gsl_blas_dznrm2(const gsl_vector_complex *x): 计算复数向量的欧氏距离， ∣ ∣ x ∣ ∣ 2 = ∑ i ( x i , r e a l 2 + x i , i m a g 2 ) {||x||}_2 = \sqrt{\sum_i (x_{i,real}^2 + x_{i,imag}^2)} ∣∣x∣∣2​=∑i​(xi,real2​+xi,imag2​) ​

float gsl_blas_sasum(const gsl_vector_float *x)

double gsl_blas_dasum(const gsl_vector *x)： 绝对值的和

float gsl_blas_scasum(const gsl_vector_complex_float *x)

double gsl_blas_dzasum(const gsl_vector_complex *x): 复数的1-norm, ∑ i ∣ x i r e a l ∣ + ∣ x i i m a g ∣ \sum_i |x_{i_real}| + |x_{i_imag}| ∑i​∣xir​eal​∣+∣xii​mag​∣ 对所有i求和

CBLAS_INDEX_t gsl_blas_isamax(const gsl_vector_float *x)：

CBLAS_INDEX_t gsl_blas_idamax(const gsl_vector *x)

CBLAS_INDEX_t gsl_blas_icamax(const gsl_vector_complex_float *x)

CBLAS_INDEX_t gsl_blas_izamax(const gsl_vector_complex *x): x的最大值的下标， 由元素的绝对值 (复数: |re|+|im|)排大小。

int gsl_blas_sswap(gsl_vector_float *x, gsl_vector_float *y)

… 使用的时候查询， 针对gsl向量和矩阵类的初等矩阵运算 基本都有。

Ch14 线性代数

gsl_linalg.h

14.1 LU分解

14.2 QR分解

14.4 LQ分解

14.5 QL分解

14.6 完全的正交分解

14.7 SVD奇异值分解

14.8 Cholesky分解

14.11 LDLT分解

14.12 实矩阵的Tridiagnoal分解

…

Ch15 本征值系统

Ch16 FFT 快速傅里叶变换

16.1 定义

x j = ∑ k = 0 n − 1 z k ∗ exp ⁡ ( − 2 π ∗ I ∗ j ∗ k / n ) xj = \sum_{k=0}^{n-1} z_k * \exp(-2\pi * I* j * k / n) xj=∑k=0n−1​zk​∗exp(−2π∗I∗j∗k/n) -> x = FFT{z}

z j = ( 1 / n ) ∗ ∑ k = 0 n − 1 x k ∗ exp ⁡ ( + 2 π ∗ I ∗ j ∗ k / n ) zj = (1 / n) * \sum_{k=0}^{n-1} x_k * \exp(+2\pi * I * j * k / n) zj=(1/n)∗∑k=0n−1​xk​∗exp(+2π∗I∗j∗k/n) --> x = IFFT{z}

z j b a c k w o r d s = ∑ k = 0 n − 1 x k ∗ exp ⁡ ( + 2 p i ∗ I ∗ j ∗ k / n ) zj^{backwords} = \sum_{k=0}^{n-1} x_k * \exp(+2pi * I * j * k / n) zjbackwords=∑k=0n−1​xk​∗exp(+2pi∗I∗j∗k/n) 省去了1/n因子， 因此不可逆但更快。

16.2 复数数据的FFT

数据的输入与输出: 使用浮点数的packed arrays，实部虚部依次排列:

double x[3*2]; 
gsl_complex_packed_array data = x; 
/* x[0] x[1]为第一个复数的实部和虚部， data一共3个复数 */

stride: z[stride*i] 代替z[i], 这样即使有间隔的数据，也可以inplace地作FFT。

对于gsl_vector_complex *v 数据，FFT中参数这样使用:

gsl_complex_packed_array_data = v->data; 
size_t stride = v->stride; 
size_t n = v->size; 

另外注意FFT结果与实际频率的对应关系。

16.3 Radix-2 FFT，复数数据

数据长度是2^n形式。 gsl_fft_complex.h中

int gsl_fft_complex_radix2_forward(gsl_complex_packed_array data, size_t stride, size_t n)

int gsl_fft_complex_radix2_transform(gsl_complex_packed_array data, size_t stride, size_t n,
gsl_fft_direction sign)

int gsl_fft_complex_radix2_backward(gsl_complex_packed_array data, size_t stride, size_t n)

int gsl_fft_complex_radix2_inverse(gsl_complex_packed_array data, size_t stride, size_t n): forward, backward, inverse FFT， 长度n， stride为stride, data是packed complex array，使用的是in-place radix-2 decimation-in-time算法。符号sign: forward(-1), backward(+1)

int gsl_fft_complex_radix2_dif_forward(gsl_complex_packed_array data, size_t stride, size_t n)
int gsl_fft_complex_radix2_dif_transform(gsl_complex_packed_array data, size_t stride, size_t n,
gsl_fft_direction sign)
int gsl_fft_complex_radix2_dif_backward(gsl_complex_packed_array data, size_t stride, size_t n)
int gsl_fft_complex_radix2_dif_inverse(gsl_complex_packed_array data, size_t stride, size_t n)：使用的是in-place radix-2 decimation-in-frequency算法的版本

16.4 Mixed-radix FFT， 复数数据

16.5-16.6 实数数据的FFT

Ch17 数值积分

在gsl_integration.h中

介绍了多种数值积分的方案， 针对的类型仍然是gsl的基础类型，需要时查阅。

Ch18 随机数发生器

在gsl_rng.h中

提供统一的方法，改变环境变量可以切换不同的随机数源，可以多线程，提供额外的将出均匀分布到高斯、泊松等分布的函数。

18.2 随机数发生器接口

type gsl_rng_type

type gsl_rng： 使用这两个特殊的结构体。gsl_rng_type保存了每种发生器的一些静态static数据。gsl_rng描述了发生器的一个实例，由某种gsl_rng_type创建而来。

18.3 随机数发生器初始化

gsl_rng *gsl_rng_alloc(const gsl_rng_type *T): 返回指向一个新创建的T类型随机数发生器的实例的指针。如 Tauworthe发生器: gsl_rng *r = gsl_rng_alloc(gsl_rng_taus); 将会使用默认的种子gsl_rng_default_seed初始化，默认为0. 可以直接被更改，或者通过GSL_RNG_SEED更改。

void gsl_rng_set(const gsl_rng *r, unsigned long int s): 从此处开始，用s作为r的种子。若s>=1，则不同s得到的随机数流不同。

void gsl_rng_free(gsl_rng *r): 释放

18.4 从随机数发生器获取样本点

均匀分布。 HAVE_INLINE

unsigned long int gsl_rng_get(const gsl_rng *r): 从r返回随机的整数。 [min, max], min和max可以由gsl_rng_max() gsl_rng_min()函数得到。

double gsl_rng_uniform(const gsl_rng *r): 返回[0, 1)区间的double型随机数，均匀分布。

double gsl_rng_uniform_pos(const gsl_rng *r): (0, 1)区间

unsigned long int gsl_rng_uniform_int(const gsl_rng *r, unsigned long int n): [0, n-1]整数，随机数

18.5 随机数发生器辅助函数

const char *gsl_rng_name(const gsl_rng *r): 名字

unsigned long int gsl_rng_max(const gsl_rng *r)

unsigned long int gsl_rng_min(const gsl_rng *r): gsl_rng_get()最大可返回值

void *gsl_rng_state(const gsl_rng *r)

size_t gsl_rng_size(const gsl_rng *r): 返回指向r的state与大小的指针

const gsl_rng_type **gsl_rng_types_setup(void): 指向一个所有可用的随机数发生器类型的数组 的指针，结束于null pointer. 例如

const gsl_rng_type **t, **t0;
t0 = gsl_rng_types_setup ();
printf ("Available generators:\n");
for (t = t0; *t != 0; t++)
{
printf ("%s\n", (*t)->name);
}

18.6 随机数发生器环境变量

GSL_RNG_TYPE：表明默认的发生器，是发生器的name.

GSL_RNG_SEED: 默认的种子

gsl_rng_type *gsl_rng_default: global, 默认的发生器，可以从GSL_RNG_TYPE 初始化，gsl_rng_setup() 。 extern const gsl_rng_type *gsl_rng_default

unsigned long int gsl_rng_default_seed: global, 默认0, 可以由GSL_RNG_SEED初始化， gsl_rng_seutp()。extern unsigned long int gsl_rng_default_seed

const gsl_rng_type *gsl_rng_env_setup(void)： 读取环境变量GSL_RNG_TYPE和GSL_RNG_SEED，设置相关的gsl_rng_default和gsl_rng_default_seed。若无这俩环境变量，使用gsl_rng_my19937, 0.

例如

#include <stdio.h>
#include <gsl/gsl_rng.h>
gsl_rng * r; /* global generator */
int
main (void)
{
    const gsl_rng_type * T;
    gsl_rng_env_setup();
    T = gsl_rng_default;
    r = gsl_rng_alloc (T);
    printf ("generator type: %s\n", gsl_rng_name (r));
    printf ("seed = %lu\n", gsl_rng_default_seed);
    printf ("first value = %lu\n", gsl_rng_get (r));
    gsl_rng_free (r);
    return 0;
}
/*
generator type: mt19937
seed = 0
first value = 4293858116
*/

18.7-18.11 几种随机数发生器的算法与一些函数

18.12 性能

最快的: taus, gfsr4, mt19937 最随机的: RANLUX

1754 k ints/sec, 870 k doubles/sec, taus
1613 k ints/sec, 855 k doubles/sec, gfsr4
1370 k ints/sec, 769 k doubles/sec, mt19937
565 k ints/sec, 571 k doubles/sec, ranlxs0
400 k ints/sec, 405 k doubles/sec, ranlxs1
490 k ints/sec, 389 k doubles/sec, mrg
407 k ints/sec, 297 k doubles/sec, ranlux
243 k ints/sec, 254 k doubles/sec, ranlxd1
251 k ints/sec, 253 k doubles/sec, ranlxs2
238 k ints/sec, 215 k doubles/sec, cmrg
247 k ints/sec, 198 k doubles/sec, ranlux389
141 k ints/sec, 140 k doubles/sec, ranlxd2

18.3 例子

[0.0, 1.0)的均匀随机数

#include <stdio.h>
#include <gsl/gsl_rng.h>
int
main (void)
{
    const gsl_rng_type * T;
    gsl_rng * r;
    int i, n = 10;
    gsl_rng_env_setup();
    T = gsl_rng_default;
    r = gsl_rng_alloc (T);
    for (i = 0; i < n; i++)
    {
        double u = gsl_rng_uniform (r);
        printf ("%.5f\n", u);
    }
    gsl_rng_free (r);
    return 0;
}

Ch19 准随机序列

高维度上的均匀分布(或一些特殊的、略带随机的)随机数。low-discrepancy sequences。gsl_qrng.h。

Ch20 随机数分布

介绍了多个随机分布函数型， 提*生该分布下随机数的方法，提供概率计算所需要的积分。

20.40 随机重排shuffle和采样

void gsl_ran_shuffle(const gsl_rng *r, void *base, size_t n, size_t size): 随机生成n!种排列中的一个，排列base。

int gsl_ran_choose(const gsl_rng *r, void *dest, size_t k, void *src, size_t n, size_t size): 从src[0~n-1]中随机选择对象，大小均为size, 对象只会在dest中出现一次。dest中的顺序与src中一致。gsl_ran_shuffle(r, dest, n, size)则附带乱序

void gsl_ran_sample(const gsl_rng *r, void *dest, size_t k, void *src, size_t n, size_t size): 与choose相比，可能出现一次或多次。

例如:

#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
int
main (void)
{
    const gsl_rng_type * T;
    gsl_rng * r;
    int i, n = 10;
    double mu = 3.0;
    /* create a generator chosen by the
    environment variable GSL_RNG_TYPE */
    gsl_rng_env_setup();
    T = gsl_rng_default;
    r = gsl_rng_alloc (T);
    /* print n random variates chosen from
    the poisson distribution with mean
    parameter mu */
    for (i = 0; i < n; i++)
    {
        unsigned int k = gsl_ran_poisson (r, mu);
        printf (" %u", k);
    }
    printf ("\n");
    gsl_rng_free (r);
    return 0;
}

Ch21 统计

gsl_statistics_double.h gsl_statistics_int.h …

double gsl_stats_mean(const double data[], size_t stride, size_t n): mean

double gsl_stats_variance(const double data[], size_t stride, size_t n)： 方差, 注意是1/(N-1)的归一化因子。

double gsl_stats_variance_m(const double data[], size_t stride, size_t n, double mean): mean给出，求方差

double gsl_stats_sd(const double data[], size_t stride, size_t n)

double gsl_stats_sd_m(const double data[], size_t stride, size_t n, double mean)：标准差， 方差开根号

21.3 高阶矩moments

…

Ch22 实时统计

可以随着数据的积累来统计。

Ch23 移动窗口统计

Ch24 数字滤波器

Ch25 histogram

gsl_histogram.h gsl_histogram2d.h

25.1-25.2 histogram结构体

type gsl_histogram： size_t n, bins的数目； double * range，bins的范围，长度n+1；double * bin，每个bin下的counts数。 bin[i] --> range[i] <= x < range[i+1]

gsl_histogram *gsl_histogram_alloc(size_t n)：获取

int gsl_histogram_set_ranges(gsl_histogram *h, const double range[], size_t size)

int gsl_histogram_set_ranges_uniform(gsl_histogram *h, double xmin, double xmax): range在[xmin, xmax)范围均匀分n个

void gsl_histogram_free(gsl_histogram *h)

int gsl_histogram_memcpy(gsl_histogram *dest, const gsl_histogram *src): 拷贝

gsl_histogram *gsl_histogram_clone(const gsl_histogram *src): 克隆

25.4 元素更新与访问

int gsl_histogram_increment(gsl_histogram *h, double x): 把x所属的bin计数加一。返回0成功，返回GSL_EDOM若不属于整个范围。

int gsl_histogram_accumulate(gsl_histogram *h, double x, double weight): 计数增量带有权值。

double gsl_histogram_get(const gsl_histogram *h, size_t i): 第i个bin的内容

int gsl_histogram_get_range(const gsl_histogram *h, size_t i, double *lower, double *upper): 返回第i个bin的上下限

double gsl_histogram_max(const gsl_histogram *h)

double gsl_histogram_min(const gsl_histogram *h)

size_t gsl_histogram_bins(const gsl_histogram *h)

void gsl_histogram_reset(gsl_histogram *h): counts --> 0

int gsl_histogram_find(const gsl_histogram *h, double x, size_t *i)： 找到x所属于的bing下标i. 失败则返回GSL_EDOM，成功GSL_SUCCESS

25.6 histogram统计

double gsl_histogram_max_val(const gsl_histogram *h)

size_t gsl_histogram_max_bin(const gsl_histogram *h)

double gsl_histogram_min_val(const gsl_histogram *h)

size_t gsl_histogram_min_bin(const gsl_histogram *h)

double gsl_histogram_mean(const gsl_histogram *h)

double gsl_histogram_sigma(const gsl_histogram *h)

double gsl_histogram_sum(const gsl_histogram *h)

25.7 histogram操作

int gsl_histogram_equal_bins_p(const gsl_histogram *h1, const gsl_histogram *h2): 是否相同

int gsl_histogram_add(gsl_histogram *h1, const gsl_histogram *h2): h1 h2的range完全一样，然后h1 = h1 + h2

int gsl_histogram_sub(gsl_histogram *h1, const gsl_histogram *h2)

int gsl_histogram_mul(gsl_histogram *h1, const gsl_histogram *h2)

int gsl_histogram_div(gsl_histogram *h1, const gsl_histogram *h2)

int gsl_histogram_scale(gsl_histogram *h, double scale)： h1 = h1 * scale

int gsl_histogram_shift(gsl_histogram *h, double offset): h1 = h1 + offset

25.10：histogram的概率分布结构体

…

25.12 - 25.22 二维histogram

…

Ch26 N-tuples N-元组

gsl_ntuple.h

26.1 ntuple结构

type gsl_ntuple：ntuple数据存储的文件file，指向当前ntuple数据行的指针，用户定义结构的size

typedef struct {FILE * file; void * ntuple_data; size_t size;} gsl_ntuple; 

26.2 ntuple创建

gsl_ntuple *gsl_ntuple_create(char *filename, void *ntuple_data, size_t size): 创建只写的ntuple文件filename, 针对于大小为size的ntuples，返回这个新建的ntuple结构的指针。已存在==> 截断到0并重写。内存地址void *ntuple_data必须提供，才能从内存IO到文件。

26.3 打开已有

gsl_ntuple *gsl_ntuple_open(char *filename, void *ntuple_data, size_t size)

26.4 写

int gsl_ntuple_write(gsl_ntuple *ntuple): 将当前ntuple->ntuple_data、大小为ntuple->size存到相应文件。

int gsl_ntuple_bookdata(gsl_ntuple *ntuple): 同上

26.5 读

int gsl_ntuple_read(gsl_ntuple *ntuple): 从文件，到ntuple->data

26.6 关闭

int gsl_ntuple_close(gsl_ntuple *ntuple): 关闭且释放临时内存

Ch27 Monte Carlo蒙特卡洛积分

Ch28 simulated annealing模拟退火

Ch29 ODE常微分方程

【TODO】

Ch30 插值

【TODO】

Ch31 数值差分

【TODO】三个差分格式

Ch32 Chebyshev Approximations 函数的Chebyshev多项式逼近

Ch33 级数加速收敛

Levin u-transform方法。

Ch34 小波变换

Ch35 离散Hankel变换

Ch36 一维寻根

【TODO】

Ch37 一维最小化

【TODO】

Ch38 多维寻根

【TODO】

Ch39 多维最小化

【TODO】

Ch40 线性最小二乘法拟合

【TODO】

Ch41 非线性最小二乘法拟合

【TODO】

Ch42 B-Splines

Ch43 稀疏矩阵

Ch44 稀疏BLAS支持

Ch45 稀疏矩阵线性代数

CH46 物理常数

standardMKSAsystem (meters, kilograms,seconds, amperes) gsl_const_mksa.h

CGSM system (centimeters, grams, seconds, gauss), gsl_const_cgsm.h

Dimensionless constants, gsl_const_num.h

46.1 基本常数

GSL_CONST_MKSA_SPEED_OF_LIGHT： c

GSL_CONST_MKSA_VACUUM_PERMEABILITY：\mu_0

GSL_CONST_MKSA_VACUUM_PERMITTIVITY：\varepsilon_0

GSL_CONST_MKSA_PLANCKS_CONSTANT_H

GSL_CONST_MKSA_PLANCKS_CONSTANT_HBAR

GSL_CONST_NUM_AVOGADRO

GSL_CONST_MKSA_FARADAY

GSL_CONST_MKSA_BOLTZMANN

GSL_CONST_MKSA_MOLAR_GAS: R_0

GSL_CONST_MKSA_STANDARD_GAS_VOLUME: V_0

GSL_CONST_MKSA_STEFAN_BOLTZMANN_CONSTANT

GSL_CONST_MKSA_GAUSS: 1Gauss = ? magnetic field

46.2 天文与天文物理

GSL_CONST_MKSA_ASTRONOMICAL_UNIT： 1au=?m

GSL_CONST_MKSA_GRAVITATIONAL_CONSTANT: G

GSL_CONST_MKSA_LIGHT_YEAR: ly in m

GSL_CONST_MKSA_PARSEC

GSL_CONST_MKSA_GRAV_ACCEL: std g,

GSL_CONST_MKSA_SOLAR_MASS

46.3 原子和核物理

GSL_CONST_MKSA_ELECTRON_CHARGE： e

GSL_CONST_MKSA_ELECTRON_VOLT: eV

GSL_CONST_MKSA_UNIFIED_ATOMIC_MASS: atomic mass, 1amu

GSL_CONST_MKSA_MASS_ELECTRON

GSL_CONST_MKSA_MASS_MUON

GSL_CONST_MKSA_MASS_PROTON

GSL_CONST_MKSA_MASS_NEUTRON

GSL_CONST_NUM_FINE_STRUCTURE

GSL_CONST_MKSA_RYDBERG

GSL_CONST_MKSA_BOHR_RADIUS

GSL_CONST_MKSA_ANGSTROM

…

GSL_CONST_MKSA_BOHR_MAGNETON

GSL_CONST_MKSA_NUCLEAR_MAGNETON

GSL_CONST_MKSA_ELECTRON_MAGNETIC_MOMENT

GSL_CONST_MKSA_THOMSON_CROSS_SECTION

GSL_CONST_MKSA_DEBYE

46.4 time

GSL_CONST_MKSA_MINUTE

GSL_CONST_MKSA_HOUR

GSL_CONST_MKSA_DAY

GSL_CONST_MKSA_WEEK

46.6 单位

GSL_CONST_MKSA_INCH

GSL_CONST_MKSA_FOOT

GSL_CONST_MKSA_YARD

GSL_CONST_MKSA_MILE

GSL_CONST_MKSA_MIL

…

GSL_CONST_MKSA_POUND_MASS

GSL_CONST_MKSA_OUNCE_MASS

GSL_CONST_MKSA_TON

GSL_CONST_MKSA_METRIC_TON

GSL_CONST_MKSA_CALORIE

GSL_CONST_MKSA_BAR

GSL_CONST_MKSA_STD_ATMOSPHERE

GSL_CONST_MKSA_TORR

GSL_CONST_MKSA_METER_OF_MERCURY

GSL_CONST_MKSA_PSI

…

Ch47 IEEE浮点数运算

Ch48 debug数值程序

Ch49 Contributors to GSL

Ch50 自动配置宏 autoconf macros

Ch51 GSL CBLAS library

Ch52 GNU GSL >=v3

Ch53 GNU Free Documentation License

---