#include <cmath>
#include <gsl/gsl_integration.h>

template<typename F>
class Integrator
{
private:
    gsl_integration_workspace *w;
    gsl_function gF;
public:
    bool use_log;
    int limit, key;
    double epsabs, epsrel, result, error;
    void setLog(bool t_use_log=true) {
        use_log = t_use_log;
        if (use_log) {
            gF.function = &functor_gsl_wrapper_log;
        } else {
            gF.function = &functor_gsl_wrapper;
        }
    }
    static double functor_gsl_wrapper_log(double logx, void* functor)
    {
        F* f = (F *) functor;
        double x = std::exp(logx);
        return (*f)(x) * x;
    }
    static double functor_gsl_wrapper(double x, void* functor)
    {
        F* f = (F *) functor;
        return (*f)(x);
    }
    Integrator(bool t_log=false, double t_epsabs=0.0, double t_epsrel=1e-4,
               int t_limit=100000, int t_key=4)
        : use_log {t_log}, limit {t_limit}, key {t_key},
          epsabs {t_epsabs}, epsrel {t_epsrel}, result {0}, error {0} {
        w = gsl_integration_workspace_alloc(t_limit);
        setLog(use_log);
    }
    double operator()(const F& f, double t_lo, double t_up) {
        gF.params = (void *) &f;

        double lo = use_log ? std::log(t_lo) : t_lo;
        double up = use_log ? std::log(t_up) : t_up;
        gsl_integration_qag(&gF, lo, up,
                            epsabs, epsrel, limit, key, w, &result, &error);
        return result;
    }
    ~Integrator() {
        gsl_integration_workspace_free(w);
    }
};
