C++使用GSL蒙特卡洛积分计算多维函数的参数类型报错问题

1. 函数形参数组退化
   C/C++语法规则中，数组作为函数形参声明时会自动退化为对应元素类型的指针，你声明的void *p[]作为g函数的形参，实际等效于void **p，和GSL要求的gsl_monte_function.f的第三个参数类型void *完全不匹配，这是编译报错的直接原因。
2. 额外参数传递错误
   gsl_monte_function的params成员本身是void *类型，不能直接用初始化列表存放多个double参数，需要先将所有要传入的参数打包到自定义结构体中，再传递结构体地址。
3. 积分上下限参数错误
   GSL的gsl_monte_plain_integrate要求第二、第三个参数为double *类型，分别指向存放各维度积分下限、上限的数组，你当前传入单个double值也会引发后续编译错误。

• 步骤1：自定义参数结构体，封装所有需要传入被积函数的参数
• 步骤2：修改g函数的签名匹配GSL的类型要求，内部将void *类型的入参转换为自定义结构体指针读取参数
• 步骤3：修改integrate_test中的参数传递逻辑，将打包好的参数结构体地址赋值给G.params，同时处理4维积分的上下限数组

#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_plain.h>
#include <gsl/gsl_rng.h>
#include <cmath>

// 自定义参数结构体，封装所有需要传入被积函数的参数
struct IntegralParams {
    double xA0;
    double xB0;
    double yA0;
    double yB0;
    double sigma;
    double t1;
};

double brownian(const double &x,const double &x0,const double & sigma,const double & t) {
    double a1 = (1/sqrt(2.0 * M_PI * sigma * t));
    double b1 = exp(-((x - x0) * (x - x0))/(2.0 * sigma * t));
    double res = a1 * b1;
    return res;
}

// 修改g函数签名匹配GSL要求
double g(double *k, size_t dim, void *p){
    const IntegralParams* params = static_cast<const IntegralParams*>(p);
    double temp_pbxA = brownian(k[0], params->xA0, params->sigma, params->t1);
    double temp_pbxB = brownian(k[1], params->xB0, params->sigma, params->t1);
    double temp_pbyA = brownian(k[2], params->yA0, params->sigma, params->t1);
    double temp_pbyB = brownian(k[3], params->yB0, params->sigma, params->t1);
    return (temp_pbxB * temp_pbyB) * (temp_pbxA * temp_pbyA);
}

double integrate_test(const double xl, const double xu, const double& xA0, const double& xB0, const double& yA0, const double& yB0, const double& t1, const double& sigma){
    double res, err;
    const gsl_rng_type *T;
    gsl_rng *r;

    gsl_monte_function G;
    G.f = &g;
    G.dim = 4;
    // 打包参数到结构体
    IntegralParams params = {xA0, xB0, yA0, yB0, sigma, t1};
    G.params = &params;

    size_t calls = 500000;

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

    // 构造4维积分的上下限数组
    double xl_arr[4] = {xl, xl, xl, xl};
    double xu_arr[4] = {xu, xu, xu, xu};

    gsl_monte_plain_state *s = gsl_monte_plain_alloc (4);
    gsl_monte_plain_integrate (&G, xl_arr, xu_arr, 4, calls, r, s, &res, &err);
    gsl_monte_plain_free (s);
    gsl_rng_free(r);

    return res;
}

编译时需要链接GSL相关依赖库，示例编译命令：
g++ 你的源代码文件名.cpp -o integrate_test -lgsl -lgslcblas -lm

内容的提问来源于stack exchange，提问作者CafféSospeso