编写接受通用 n 维函数的模板化递归集成方法

问题描述

我正在开发一个大型 C++ 框架来进行一些粒子物理计算，如果有一种方法可以将通用函数与 N 个变量进行 N 维集成，那就太好了。积分发生在 N-cube [0:1]^N 上。 不幸的是，被积函数是一个成员函数，因此我将 std::bind 对象传递给模板化函数（lambdas 也能很好地工作）。当 N > 1 时，积分方法再次绑定固定一个变量的函数，并将新函数传递给 N-1 的方法。

我觉得这可以通过将 N 作为模板参数递归地完成，但这就是我陷入困境的地方：有没有办法“模板化”绑定部分，以便它添加正确的具有函数元数或 N 的占位符数量？然后将新的绑定方法传递给 integral<N-1> 等。我认为问题在于 std::bind 对象没有签名。也许使用 lambdas 可以利用可变参数包扩展，但模板又必须以某种方式解析签名。 有没有办法使这项工作？如果可能，在 c++11 中并且没有 boost 库。

这是我现在所拥有的非常简化的版本，没有递归，而是每个维度的集成方法的不同定义（现在最多 3 个）

#include <iostream>
#include <functional>
#include <cmath>

// integrate between [0:1]
template <typename Functor>
double integral1D(Functor fn,size_t steps = 100) {
        double sum = 0.;
        double step = 1./steps;
        for (size_t s = 0; s < steps; ++s)
                sum += fn(s * step);

        return sum * step;
}

template <typename Functor>
double integral2D(Functor fn,size_t steps = 100) {
        auto subfn = [&](double v) ->double {
                auto subint = std::bind(fn,v,std::placeholders::_1);
                return integral1D(subint,steps);
        };

        return integral1D(subfn,steps);
}

template <typename Functor>
double integral3D(Functor fn,std::placeholders::_1,std::placeholders::_2);
                return integral2D(subint,steps);
}
        
struct A
{
        double gaus1(double x,double range) { // computes jacobian on [-range:range]
                x = range * (2. * x - 1.);
                return 2 * range * std::exp(- x*x);
        }
    
        double gaus2(double x,double y,double range) {
                return gaus1(x,range) * gaus1(y,range);
        }

        double gaus3(double x,double z,double range) {
                return gaus2(x,y,range) * gaus1(z,range);
        }

        double gaus1_integrate(double range) {
                auto func = std::bind(&A::gaus1,this,range);
                return integral1D(func);
        }

        double gaus2_integrate(double range) {
                auto func = std::bind(&A::gaus2,std::placeholders::_2,range);
                return integral2D(func);
        }

        double gaus3_integrate(double range) {
                auto func = std::bind(&A::gaus3,std::placeholders::_3,range);
                return integral3D(func);
        }
};

int main() {
        A a;
        std::cout << "1D integral is " << a.gaus1_integrate(50) << "\n";
        std::cout << "2D integral is " << a.gaus2_integrate(50) << "\n";
        std::cout << "3D integral is " << a.gaus3_integrate(50) << "\n";
        return 0;
}

上面的例子有效并给出了预期的结果。我知道不建议对具有更多维度的更复杂的函数进行 riemann 求和（或等效），但很高兴看到上述类似的方法是否可行。

答案

基于 Guillaume 的有用建议，这里是 c++11 工作示例。需要稍作修改。

#include <iostream>
#include <functional>
#include <cmath>

template<int n,typename Functor>
struct integralNDsubint {
        double v;
        const Functor &fn;

        template<typename... Args>
                double operator()(Args... rest) const {
                        return fn(v,rest...);
                }
};

template <int n,typename Functor,typename std::enable_if<(n > 1),bool>::type = true>
double integralND(Functor fn,size_t steps = 100) {
        auto subfn = [&](double v) -> double {
                auto subint = integralNDsubint<n,Functor> {v,fn};
                // following block can be used in c++14 without
                // need of helper struct integralNDsubint
                //auto subint = [v,&fn](auto... rest) {
                        //return fn(v,rest...);
                //};

                return integralND<n - 1>(subint,steps);
        };

        return integralND<1>(subfn,steps);
}

template <int n,typename std::enable_if<(n == 1),size_t steps = 100) {
        double sum = 0.;
        double step = 1./steps;
        for (size_t s = 0; s < steps; ++s) {
                sum += fn(s * step);
        }

        return sum * step;
}


struct A
{
        double gaus1(double x,double range) {
                x = range * (2. * x - 1.);
                return 2 * range * std::exp(- x*x);
        }

        double gaus2(double x,range);
        }

        double gaus1_integrate(double range) {
                auto func = [=](double x) { return gaus1(x,range); };
                return integralND<1>(func);
        }

        double gaus2_integrate(double range) {
                auto func = [=](double x,double y) { return gaus2(x,range); } ;
                return integralND<2>(func);
        }

        double gaus3_integrate(double range) {
                auto func = [=](double x,double z) { return gaus3(x,z,range); } ;
                return integralND<3>(func);
        }
};

int main() {
        A a;
        std::cout << "1D integral is " << a.gaus1_integrate(50) << "\n";
        std::cout << "2D integral is " << a.gaus2_integrate(50) << "\n";
        std::cout << "3D integral is " << a.gaus3_integrate(50) << "\n";
        return 0;
}

解决方法

您的函数可以泛化为任意数量的维度。

template <int n,typename Functor,std::enable_if_t<(n > 1),int> = 0>
double integralND(Functor fn,size_t steps = 100) {
    auto subfn = [&](double v) -> double {
        auto subint = [v,&fn](auto... rest) {
            return fn(v,rest...);
        };

        return integralND<n - 1>(subint,steps);
    };

    return integralND<1>(subfn,steps);
}

template <int n,std::enable_if_t<(n == 1),size_t steps = 100) {
    double sum = 0.;
    double step = 1./steps;
    for (size_t s = 0; s < steps; ++s) {
        sum += fn(s * step);
    }

    return sum * step;
}

如您所见，它使用通用 lambda。为了与 C++11 兼容，您可以手动编写函数对象：

template<int n,typename F>
struct integralNDsubint {
    double v;
    F const& fn;

    template<typename... Args>
    auto operator()(Args... rest) const -> double {
        return fn(v,rest...);
    }
};

并像这样使用它：

auto subint = integralNDsubint<n,Functor>{v,fn};

您还需要更改 sfinae：

template <int n,typename std::enable_if<(n > 1),int>::type = 0>
double integralND(Functor fn,size_t steps = 100) {
    // ...
}

template <int n,typename std::enable_if<(n == 1),size_t steps = 100) {
    // ...
}