#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_miser.h>
#include <gsl/gsl_monte_plain.h>
#include <gsl/gsl_monte_vegas.h>
double exact = 1.3932039296856768591842462603255;
double g(double *k, size_t dim, void *params) {
(void)(dim);
(void)(params);
double A = 1.0 / (M_PI * M_PI * M_PI);
return A / (1.0 - cos(k[0]) * cos(k[1]) * cos(k[2]));
}
void display_results(char *title, double result, double error) {
printf("%s ==================\n", title);
printf("result = % .6f\n", result);
printf("sigma = % .6f\n", error);
printf("exact = % .6f\n", exact);
printf("error = % .6f = %.2g sigma\n", result - exact,
fabs(result - exact) / error);
}
int main(void) {
double res, err;
double xl[3] = {0, 0, 0};
double xu[3] = {M_PI, M_PI, M_PI};
const gsl_rng_type *T;
gsl_rng *r;
gsl_monte_function G = {&g, 3, 0};
size_t calls = 500000;
gsl_rng_env_setup();
T = gsl_rng_default;
r = gsl_rng_alloc(T);
{
gsl_monte_plain_state *s = gsl_monte_plain_alloc(3);
gsl_monte_plain_integrate(&G, xl, xu, 3, calls, r, s, &res, &err);
gsl_monte_plain_free(s);
display_results("plain", res, err);
}
{
gsl_monte_miser_state *s = gsl_monte_miser_alloc(3);
gsl_monte_miser_integrate(&G, xl, xu, 3, calls, r, s, &res, &err);
gsl_monte_miser_free(s);
display_results("miser", res, err);
}
{
gsl_monte_vegas_state *s = gsl_monte_vegas_alloc(3);
gsl_monte_vegas_integrate(&G, xl, xu, 3, 10000, r, s, &res, &err);
display_results("vegas warm-up", res, err);
printf("converging...\n");
do {
gsl_monte_vegas_integrate(&G, xl, xu, 3, calls / 5, r, s, &res, &err);
printf(
"result = % .6f sigma = % .6f "
"chisq/dof = %.1f\n",
res, err, gsl_monte_vegas_chisq(s));
} while (fabs(gsl_monte_vegas_chisq(s) - 1.0) > 0.5);
display_results("vegas final", res, err);
gsl_monte_vegas_free(s);
}
gsl_rng_free(r);
return 0;
}