5月28日追记 其实现在最新版本的Boost库数值积分已经比较完善了，处理常见的绝大多数积分都没有问题，因此C++的数值积分也可以考虑使用Boost，它的接口更加友好、使用也更方便，不需要像这篇文章一样再进行封装了。
C/C++中的数值积分在不考虑引用库Licence类型情况下，可以使用GNU Scientific Library (GSL)这个库来进行处理。如果只是需要计算一两种固定种类的积分，那么直接按照文档中的例子来写没什么问题。但是如果需要计算的积分种类很多的话，比如一个积分中还嵌套着其他好几个积分时，按照文档的写法写出来代码会变得非常丑陋、不可读以及难以维护。因为每算一个积分都需要构造一个函数double f (double x, void * params)，这会造成两个问题：一是函数定义和使用不在同一个位置，由于这个函数是专门为GSL积分所准备的，因此通常只用一次就不会再用到了，这样读代码和维护的时候非常麻烦；二是传参的时候需要自己将一个或多个参数装配成void * params的形式传递进函数，之后再在函数中将自定义参数从void * params中释放出来。
为了解决这两个问题，我决定利用C++的新特性对GSL的数值积分进行一下简单的封装。第一个问题非常适合用lambda表达式来解决，这样定义和使用都在同一个位置，更利于理解和维护。第二个问题其实是可变参数的问题，最适合我们这种情况的是可变模板参数，这样我们既不用手工装配也不用手工释放参数，lambda表达式传进来参数就可以直接使用。
所以我们的模板类可能是这个样子的：
template <typename... T> class Quadgk {
private:
double (*func)(double x, T... params);
double (*gslFunc)(double x, void *params);
public:
Quadgk(double (*_func)(double x, T... params));
double operator()(double a, double b, T... params);
};
不过这样我们遇到了一个难题，由于我们最终是调用GSL进行数值积分运算，因此传给GSL的是double (*gslFunc)(double x, void *params)这种形式的函数，而我们lambda表达式的形式是double (*func)(double x, T... params)那这就意味着要在gslFunc中调用func。但是怎么把T...通过void *传递进去，又该怎么调用这个含有可变模板参数的函数呢？
最终我们在这篇博客『Modern C++: Variadic template parameters and tuples』中找到了答案，至于原理是啥我是看不懂了，反正直接抄过来能用就对了……所以最终封装好的类长这样：
/*
* Copyright (C) 2020 Touuki
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 3 of the License, or (at
* your option) any later version.
*
* This program is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301,
* USA.
*/
#include <gsl/gsl_integration.h>
#include <tuple>
#include <utility>
template <typename... T> class Quadgk {
private:
gsl_integration_workspace *w;
double epsabs;
double epsrel;
double (*func)(double x, T... params);
double (*gslFunc)(double x, void *params);
public:
Quadgk(double (*_func)(double x, T... params), size_t _limit = 1000,
double _epsabs = 0., double _epsrel = 1e-8)
: w(gsl_integration_workspace_alloc(_limit)), epsabs(_epsabs),
epsrel(_epsrel), func(_func),
gslFunc([](double x, void *params) -> double {
void **p = reinterpret_cast<void **>(params);
Quadgk *self = reinterpret_cast<Quadgk *>(*p);
std::tuple<T...> tuple =
*reinterpret_cast<std::tuple<T...> *>(*(p + 1));
return self->callFunc(x, tuple, std::index_sequence_for<T...>());
}){};
~Quadgk() { gsl_integration_workspace_free(w); };
void setEpsabs(double _epsabs) { epsabs = _epsabs; };
void setEpsrel(double _epsrel) { epsrel = _epsrel; };
template <std::size_t... Is>
double callFunc(double x, const std::tuple<T...> &tuple,
std::index_sequence<Is...>) const {
return func(x, std::get<Is>(tuple)...);
};
double operator()(double a, double b, T... params) {
std::tuple<T...> tuple(params...);
void *p[] = {this, &tuple};
gsl_function f = {gslFunc, &p};
double result, abserr;
if (isinf(a)) {
if (isinf(b)) {
gsl_integration_qagi(&f, epsabs, epsrel, w->limit, w, &result,
&abserr);
} else {
gsl_integration_qagil(&f, b, epsabs, epsrel, w->limit, w,
&result, &abserr);
}
} else {
if (isinf(b)) {
gsl_integration_qagiu(&f, a, epsabs, epsrel, w->limit, w,
&result, &abserr);
} else {
gsl_integration_qags(&f, a, b, epsabs, epsrel, w->limit, w,
&result, &abserr);
}
}
return result;
};
};
template <typename... T> class QuadSine {
private:
gsl_integration_qawo_table *wf;
gsl_integration_workspace *workspace;
gsl_integration_workspace *cycle_workspace;
double epsabs;
double epsrel;
double (*func)(double x, T... params);
double (*gslFunc)(double x, void *params);
public:
QuadSine(double (*_func)(double x, T... params),
gsl_integration_qawo_enum _sine, size_t _n = 3,
size_t _limit = 1000, double _epsabs = 1e-8, double _epsrel = 1e-8)
: wf(gsl_integration_qawo_table_alloc(1., 1., _sine, _n)),
workspace(gsl_integration_workspace_alloc(_limit)),
cycle_workspace(gsl_integration_workspace_alloc(_limit)),
epsabs(_epsabs), epsrel(_epsrel), func(_func),
gslFunc([](double x, void *params) -> double {
void **p = reinterpret_cast<void **>(params);
QuadSine *self = reinterpret_cast<QuadSine *>(*p);
std::tuple<T...> tuple =
*reinterpret_cast<std::tuple<T...> *>(*(p + 1));
return self->callFunc(x, tuple, std::index_sequence_for<T...>());
}){};
~QuadSine() {
gsl_integration_qawo_table_free(wf);
gsl_integration_workspace_free(workspace);
gsl_integration_workspace_free(cycle_workspace);
};
void setEpsabs(double _epsabs) { epsabs = _epsabs; };
void setEpsrel(double _epsrel) { epsrel = _epsrel; };
template <std::size_t... Is>
double callFunc(double x, const std::tuple<T...> &tuple,
std::index_sequence<Is...>) const {
return func(x, std::get<Is>(tuple)...);
};
double operator()(double a, double b, double omega, T... params) {
std::tuple<T...> tuple(params...);
void *p[] = {this, &tuple};
gsl_function f = {gslFunc, &p};
double result, abserr;
if (isinf(a)) {
throw std::invalid_argument("a can't be INFINITY.");
} else {
if (isinf(b)) {
if (wf->omega != omega) {
gsl_integration_qawo_table_set(wf, omega, 1., wf->sine);
}
gsl_integration_qawf(&f, a, epsabs, workspace->limit, workspace,
cycle_workspace, wf, &result, &abserr);
} else {
double L = b - a;
if (wf->omega != omega) {
gsl_integration_qawo_table_set(wf, omega, L, wf->sine);
} else if (wf->L != L) {
gsl_integration_qawo_table_set_length(wf, L);
}
gsl_integration_qawo(&f, a, epsabs, epsrel, workspace->limit,
workspace, wf, &result, &abserr);
}
}
return result;
};
};
这样使用起来就非常简单方便了，这里给一个嵌套积分的例子：
#include <gsl/gsl_specfunc.h>
#include <iostream>
int main(int argc, char const *argv[]) {
Quadgk<double> quad([](double chi, double eta) {
static Quadgk<double> intK(
[](double x, double nu) { return gsl_sf_bessel_Knu(nu, x); });
double u = 2. * chi / (eta - 2. * chi);
double u_p = 2. * u / (3. * eta);
return (1. / (1. + u) + 1. + u) * gsl_sf_bessel_Knu(2. / 3., u_p) -
intK(u_p, INFINITY, 1. / 3.);
});
double eta = 10;
std::cout << quad(0., eta / 2., eta) << std::endl;
return 0;
};