#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include "gsl/gsl_math.h"
#include "gsl/gsl_monte.h"
#include "gsl/gsl_monte_vegas.h"

#define check_finite(v) assert(gsl_finite(v) && !gsl_isnan(v))

static const double inf = 100;

double hypot_sq(double x, double y) {
    return gsl_pow_2(x) + gsl_pow_2(y);
}

double integrand(double* coords, size_t ndims, void* params) {
    double k = 0.8;
    double k2 = k*k;
    double lx = coords[0];
    double ly = coords[1];
    double lpx = coords[2];
    double lpy = coords[3];
    double lppx = coords[4];
    double lppy = coords[5];
    double F1 = exp(-hypot_sq(lppx, lppy));
    double F2 = exp(-hypot_sq(k - lpx - lppx, -lpy - lppy));
    double F3 = exp(-hypot_sq(k - lx - lppx, -ly - lppy));
    double result;
    assert(ndims == 6);
    
    result = F1 * F2 * F3;
    check_finite(result);
    return result;
}


int main(const int argc, const char** argv) {
    double res, err;
    double min[] = {-inf, -inf, -inf, -inf, -inf, -inf};
    double max[] = { inf,  inf,  inf,  inf,  inf,  inf};
    const size_t ndim = 6;
    const gsl_rng_type* T;
    gsl_rng *r;

    T = gsl_rng_default;
    r = gsl_rng_alloc(T);
    gsl_rng_set(r, 100);

    gsl_monte_function function_package;
    function_package.f = integrand;
    function_package.dim = ndim;
    function_package.params = NULL;
    gsl_monte_vegas_state *s = gsl_monte_vegas_alloc(ndim);
    gsl_monte_vegas_integrate (&function_package, (double*)min, (double*)max, (size_t)ndim, 50000, r, s, &res, &err);
    fprintf(stderr, "Expected output to demonstrate bug:\n result = -nan; errbound =  0.000000; chisq/dof = nan\n");
    fprintf(stderr, "Actual output:\n result = % .6f; errbound = % .6f; chisq/dof = %.2f\n", res, err, gsl_monte_vegas_chisq(s));

    gsl_monte_vegas_free(s);
    gsl_rng_free(r);
    return 0;
}