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

I am working on a large C++ framework to do some particle physics computation and it would be nice to have a method to do N-dimensional integration of a generic function with N variables. The integration happens on the N-cube [0:1]^N. The integrand is unfortunately a member function and so I am passing an std::bind object to a templated function (lambdas work well too). When N > 1, the integration method binds again the function fixing one variable and passes the new function to the method for N-1.

I have the feeling that this can be done recursively by having N as a template parameter, but this is where I get stuck: is there a way to "templetize" the binding part such that it adds the right amount of placeholders with either the arity of the function or N? Then pass the new bound method to integral<N-1>, etc. I believe the problem lies in that an std::bind object has no signature. Maybe with lambdas one could exploit variadic pack expansions, but again the template has to resolve the signature somehow. Is there a way to make this work? If possible in c++11 and without boost libraries.

This is a very simplified version of what I have now without recursions, but instead different definitions of the integration method for each dimension (up to 3 now)

#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;
}

The above example works and gives expected results. I know it is not recommended to do Riemann sums (or equivalent) for more complicated functions with more dimensions, but it would be nice to see if something like described above can work.