Wrapper C ++ para algoritmo de localização de raiz GSL com derivada

Portanto, enquanto estou feliz em encontrar muitas respostas no Stack Overflow, decidi que é hora de fazer uma pergunta.
Estou tentando usarum algoritmo de localização de raiz com derivadas. De acordo com o GSL, tenho que definir a função e sua derivada antecipadamente. Mas eu me pergunto se isso pode ser feito mais elegante com um invólucro.

Algum tempo atrás, eu achei um muito útilmodelo (wrapper GSL C ++) que funciona bem para uma função, p. integrar e fazer uso pesado dele.

Agora, estou me perguntando se essa abordagem pode ser estendida para fornecer duas funções para o GSL, a saber, a própria função e sua derivada.

Editar: Solução

template <typename F, typename G>
class gsl_root_deriv : public gsl_function_fdf
{
private:
    const F&    _f;
    const G&    _df;

    static double invoke_f(double x, void* params){
        return static_cast<gsl_root_deriv*>(params) -> _f(x);
    }

    static double invoke_df(double x, void* params){
        return static_cast<gsl_root_deriv*>(params) -> _df(x);
    }

    static void     invoke_fdf(double x, void* params, double* f, double* df){
        (*f)    = static_cast<gsl_root_deriv*>(params)  ->  _f(x);
        (*df)   = static_cast<gsl_root_deriv*>(params)  ->  _df(x);
    }

public:
    gsl_root_deriv(const F& f_init, const G& df_init)
    : _f(f_init), _df(df_init)
    {
        f               = &gsl_root_deriv::invoke_f;
        df          = &gsl_root_deriv::invoke_df;
        fdf         = &gsl_root_deriv::invoke_fdf;
        params  = this;
    }
};

E um exemplo mínimo que se assemelhao exemplo do GSL:

#include <iostream>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_roots.h>
#include <memory>

#include "gsl_root_deriv.h"

int main()
{
auto f = [](double x) -> double { return 0.25 * x*x - 1.0; };
auto df = [](double x) -> double { return 0.5 * x; };

gsl_root_deriv<decltype(f),decltype(df)> Fp(f,df);

int status(0), iter(0), max_iter(100);

const gsl_root_fdfsolver_type* T( gsl_root_fdfsolver_newton);

std::unique_ptr<gsl_root_fdfsolver,void(*)(gsl_root_fdfsolver*)>
    s(gsl_root_fdfsolver_alloc(T),gsl_root_fdfsolver_free);

double x_0(0.0), x(5.0);

gsl_root_fdfsolver_set( s.get(), &Fp, x );

do
{
    iter++;
    std::cout << "Iteration " << iter << std::endl;
    status = gsl_root_fdfsolver_iterate( s.get() );
    x_0 = x;
    x = gsl_root_fdfsolver_root( s.get() );
    status = gsl_root_test_delta( x, x_0, 0.0, 1.0e-3 );
} while( status == GSL_CONTINUE && iter < max_iter );

std::cout << "Converged to root " << x << std::endl;

return 0;
}

Atenciosamente,
- Klaus

questionAnswers(1)

yourAnswerToTheQuestion