// Example program
#include <iostream>
#include <string>

#include <tuple>
#include <functional>

template<typename... Fn>
class LinearCombination {
public:
    template<typename... Un>
    LinearCombination(Un&&... fs)
    : functions(std::forward<Un>(fs)...)
    {
        coefs.fill(0.0);
    }
    
    double operator()(double x) const 
    {
        return evaluateImpl(x, std::integral_constant<size_t, sizeof...(Fn) - 1>());
    }
    
    void setCoef(size_t i, double c)
    {
        coefs[i] = c;
    }
    
private:
    template<size_t I>
    double evaluateImpl(double x, std::integral_constant<size_t, I> index) const
    {
        return evaluateOne(x, index) + evaluateImpl<I - 1>(x, std::integral_constant<size_t, I - 1>());
    }
    
    template<size_t I>
    double evaluateImpl(double x, std::integral_constant<size_t, 0> index) const
    {
        return evaluateOne(x, index);
    }

    template<size_t I>
    double evaluateOne(double x, std::integral_constant<size_t, I>) const
    {
        auto coef = coefs[I];
        return coef == 0.0 ? 0.0 : coef * std::get<I>(functions)(x);
    }
    

    std::tuple<Fn...> functions;
    std::array<double, sizeof...(Fn)> coefs;
};

template<typename... Fn>
auto make_linear_combination(Fn&&... fn)
{
    return LinearCombination<Fn...>{std::forward<Fn>(fn)...};
}

double A(double) 
{
    return 1.0;
}

/// Integration 1.0
double integrate3D(double (*f)(double, double, double))
{
    return f(1, 2, 3);
}

struct YHelper {
    static double Y(double x, double y, double z)
    {
        return (f(x+y) - f(x)) * (f(y+z) - f(y)) * A(z+y);
    }
    
    static std::function<double(double)> f;
};

std::function<double(double)> YHelper::f;


/// Integration 2.0
template<typename Integrable>
double integrate3D_2(Integrable&& f)
{
    return f(1, 2, 3);
}

int main()
{
    auto f1 = [](double x) { return x; };
    auto f2 = [](double x) { return 2 *x; };
    
    auto lc = make_linear_combination(std::move(f1), std::move(f2));
    lc.setCoef(0, 1.0);
    lc.setCoef(1, -1.0);
    
    std::cout << lc(2.0) << "\n";

    YHelper::f = std::ref(lc);
    std::cout << integrate3D(&YHelper::Y) << "\n";

    auto Y = [&lc](double x, double y, double z) { return (lc(x+y) - lc(x)) * (lc(y+z) - lc(y)) * A(z+y); };
    std::cout << integrate3D_2(Y) << "\n";
}