#include <boost/math/tools/roots.hpp>
#include <boost/math/distributions/gamma.hpp>

typedef boost::math::policies::policy<boost::math::policies::digits10<16> > acc_policy;
typedef boost::math::gamma_distribution<double, acc_policy> gamma_dist;
typedef std::pair<double, double> result_t;

int main() {
    const long unsigned int u = 55;
    const long unsigned int n = 2373;
    const double nu_a = 6456.947055550609;
    const double alpha = 27.21005923114458;
    const double nu_m = ((double)u/n)*nu_a;

    const gamma_dist g_m(nu_m, 1.0/alpha);
    const gamma_dist g_a(nu_a, 1.0/alpha);
    auto zero_func = [&] (const double& x) { 
        return boost::math::pdf(g_m, x) - boost::math::pdf(g_a, x);
    };
    const double lower = boost::math::mean(g_m);
    const double upper = 237.3;
    boost::uintmax_t max_iter = 100;
    boost::math::tools::eps_tolerance<double> tol(16); // 10 bit precition on result
    result_t res = boost::math::tools::toms748_solve(zero_func, lower, upper, tol, max_iter);
    //result_t res = boost::math::tools::bisect(zero_func, lower, upper, tol, max_iter); // TODO: This works as expected
    return 0;
}