title: Montecarlo integral optimization
body: <p>I am using the montecarlo method as implemented in the gsl library. I need to compute many repetitions of this integral changing a parameter in the integrand. So I need to make my subroutine fast. It seems that the most time consuming part is the evaluation of the integrand at the random points. How could I make the evaluation faster in my specific case?
  Here is a minimal example:</p>
  
  <pre><code>#include &lt;gsl/gsl_rng.h&gt;
  #include &lt;math.h&gt;
  #include &lt;stdlib.h&gt;
  #include &lt;stdio.h&gt;
  #include &lt;gsl/gsl_math.h&gt;
  #include &lt;gsl/gsl_monte.h&gt;
  #include &lt;gsl/gsl_monte_plain.h&gt;
  #include &lt;gsl/gsl_monte_vegas.h&gt;
  
  double q=0.0;
  double mu=0.001;
  double eta=0.1;
  double kF=1.0;
  double Kcut=10;
  long int Nmax=10000000;
  int Nwu=1000000;
  double w=1;
  
  
  struct my_f_params { double y;};
  
  double 
  g (double *k, size_t dim, void *p)
  {
    double A;
    struct my_f_params * fp = (struct my_f_params *)p;
  
    double PQ=q*q+k[1]*k[1]-2*q*k[1]*cos(k[3])+mu;
    double QK=k[0]*k[0]+k[1]*k[1]-2*k[0]*k[1]* (cos(k[2])*cos(k[3])+cos(k[4])*sin(k[2])*sin(k[3]))+mu;
    double KPQ=q*q+k[0]*k[0]+k[1]*k[1]+2*k[0]*cos(k[2])*(q-k[1]*cos(k[3]))+2*k[1]*       (q*cos(k[3])+k[0]*cos(k[4])*sin(k[2])*sin(k[3]));
    double denFreq=fp-&gt;y-0.5*(k[0]*k[0]+k[1]*k[1]+KPQ);
    double vol=k[0]*k[0]*k[1]*k[1]*sin(k[2])*sin(k[3]);
    if (sqrt(KPQ) &lt; kF) {
      A = vol*denFreq*(1/QK-1/PQ)/(QK*(pow(denFreq,2)+eta*eta));
    }
    else {
      A = 0;
    }
  
    return A; 
  }
  
  
  int
  main (void)
  {
    double res, err;
    double xl[5] = {0, kF, 0, 0, 0};
    double xu[5] = {kF, Kcut, M_PI, M_PI, 2*M_PI};
    const gsl_rng_type *T;
    gsl_rng *r;
    gsl_monte_function G;
    size_t calls = Nmax;
  
    gsl_rng_env_setup ();
    struct my_f_params params;
    T = gsl_rng_default;
    r = gsl_rng_alloc (T);
  
    params.y=w;
    G.f=&amp;g;
    G.dim=5;
    G.params=&amp;params;
  
    {
      gsl_monte_vegas_state *s = gsl_monte_vegas_alloc (5);
      gsl_monte_vegas_integrate (&amp;G, xl, xu, 5, Nwu, r, s,&amp;res, &amp;err);
      do
        {
          gsl_monte_vegas_integrate (&amp;G, xl, xu, 5, calls/5, r, s,&amp;res, &amp;err);
        }
      while (fabs (gsl_monte_vegas_chisq (s) - 1.0) &gt; 0.5);
      gsl_monte_vegas_free (s);
    }
  
    printf ("%.6f   %.6f   %.6f\n", w,res,err);
    gsl_rng_free (r);
  
  
  
    return 0;
  }
  </code></pre>