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

Question

我正在开发一个大型 C++ 框架来进行一些粒子物理计算，如果有一种方法可以对具有 N 个变量的通用函数进行 N 维积分，那就太好了。积分发生在 N 立方体 [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, size_t steps = 100) {
        auto subfn = [&](double v) ->double {
                auto subint = std::bind(fn, v, std::placeholders::_1, std::placeholders::_2);
                return integral2D(subint, steps);
        };

        return integral1D(subfn, 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 y, double z, double range) {
                return gaus2(x, y, range) * gaus1(z, range);
        }

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

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

        double gaus3_integrate(double range) {
                auto func = std::bind(&A::gaus3, this, std::placeholders::_1,
                                                       std::placeholders::_2,
                                                       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;
}

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

回答

根据 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 Functor, typename std::enable_if<(n == 1), bool>::type = true>
double integralND(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;
}


struct A
{
        double gaus1(double x, double 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 y, double z, double range) {
                return gaus2(x, y, range) * gaus1(z, 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, y, range); } ;
                return integralND<2>(func);
        }

        double gaus3_integrate(double range) {
                auto func = [=](double x, double y, double z) { return gaus3(x, y, 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;
}

Answer 1

您的函数可以推广到任意数量的维度。

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, typename Functor, std::enable_if_t<(n == 1), int> = 0>
double integralND(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;
}

如您所见，它使用了通用的 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 Functor, typename std::enable_if<(n > 1), int>::type = 0>
double integralND(Functor fn, size_t steps = 100) {
    // ...
}

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

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

Write templated recursive integration method that accepts generic n-dimensional functions

c++

lambda

templates

stdbind

c++11