Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place.

1
#include <stdlib.h>
2
#include <gsl/gsl_math.h>
3
#include <gsl/gsl_monte.h>
4
#include <gsl/gsl_monte_miser.h>
5
#include <gsl/gsl_monte_plain.h>
6
#include <gsl/gsl_monte_vegas.h>
7

8
/* Computation of the integral,
9

10
      I = int (dx dy dz)/(2pi)^3  1/(1-cos(x)cos(y)cos(z))
11

12
   over (-pi,-pi,-pi) to (+pi, +pi, +pi).  The exact answer
13
   is Gamma(1/4)^4/(4 pi^3).  This example is taken from
14
   C.Itzykson, J.M.Drouffe, "Statistical Field Theory -
15
   Volume 1", Section 1.1, p21, which cites the original
16
   paper M.L.Glasser, I.J.Zucker, Proc.Natl.Acad.Sci.USA 74
17
   1800 (1977) */
18

19
/* For simplicity we compute the integral over the region
20
   (0,0,0) -> (pi,pi,pi) and multiply by 8 */
21

22
double exact = 1.3932039296856768591842462603255;
23

24
double g(double *k, size_t dim, void *params) {
25
  (void)(dim); /* avoid unused parameter warnings */
26
  (void)(params);
27
  double A = 1.0 / (M_PI * M_PI * M_PI);
28
  return A / (1.0 - cos(k[0]) * cos(k[1]) * cos(k[2]));
29
}
30

31
void display_results(char *title, double result, double error) {
32
  printf("%s ==================\n", title);
33
  printf("result = % .6f\n", result);
34
  printf("sigma  = % .6f\n", error);
35
  printf("exact  = % .6f\n", exact);
36
  printf("error  = % .6f = %.2g sigma\n", result - exact,
37
         fabs(result - exact) / error);
38
}
39

40
int main(void) {
41
  double res, err;
42

43
  double xl[3] = {0, 0, 0};
44
  double xu[3] = {M_PI, M_PI, M_PI};
45

46
  const gsl_rng_type *T;
47
  gsl_rng *r;
48

49
  gsl_monte_function G = {&g, 3, 0};
50

51
  size_t calls = 500000;
52

53
  gsl_rng_env_setup();
54

55
  T = gsl_rng_default;
56
  r = gsl_rng_alloc(T);
57

58
  {
59
    gsl_monte_plain_state *s = gsl_monte_plain_alloc(3);
60
    gsl_monte_plain_integrate(&G, xl, xu, 3, calls, r, s, &res, &err);
61
    gsl_monte_plain_free(s);
62

63
    display_results("plain", res, err);
64
  }
65

66
  {
67
    gsl_monte_miser_state *s = gsl_monte_miser_alloc(3);
68
    gsl_monte_miser_integrate(&G, xl, xu, 3, calls, r, s, &res, &err);
69
    gsl_monte_miser_free(s);
70

71
    display_results("miser", res, err);
72
  }
73

74
  {
75
    gsl_monte_vegas_state *s = gsl_monte_vegas_alloc(3);
76
    gsl_monte_vegas_integrate(&G, xl, xu, 3, 10000, r, s, &res, &err);
77
    display_results("vegas warm-up", res, err);
78

79
    printf("converging...\n");
80

81
    do {
82
      gsl_monte_vegas_integrate(&G, xl, xu, 3, calls / 5, r, s, &res, &err);
83
      printf(
84
          "result = % .6f sigma = % .6f "
85
          "chisq/dof = %.1f\n",
86
          res, err, gsl_monte_vegas_chisq(s));
87
    } while (fabs(gsl_monte_vegas_chisq(s) - 1.0) > 0.5);
88

89
    display_results("vegas final", res, err);
90

91
    gsl_monte_vegas_free(s);
92
  }
93

94
  gsl_rng_free(r);
95

96
  return 0;
97
}
98

99