90 lines
2.3 KiB
C
90 lines
2.3 KiB
C
#include <gsl/gsl_monte_plain.h>
|
|
#include <gsl/gsl_monte_miser.h>
|
|
#include <gsl/gsl_monte_vegas.h>
|
|
#include <gsl/gsl_monte.h>
|
|
#include <math.h>
|
|
|
|
|
|
// Exact value of the integral.
|
|
double exact = 1.71828182845905;
|
|
|
|
|
|
// Function which must be passed to the gsl_function.
|
|
double function (double * x, size_t dim, void * params){
|
|
// Avoid warning when no parameters are passed to the function.
|
|
(void)(params);
|
|
return exp(x[0]);
|
|
}
|
|
|
|
|
|
// Function which prints out the results when called.
|
|
void results (char *title, double result, double error)
|
|
{
|
|
printf ("---------------------------\n%s\n", title);
|
|
printf ("result: % .10f\n", result);
|
|
printf ("sigma: % .10f\n", error);
|
|
printf ("exact: % .10f\n", exact);
|
|
}
|
|
|
|
|
|
int main(){
|
|
// Some useful variables.
|
|
double integral, error;
|
|
size_t calls = 50000000;
|
|
size_t dims = 1;
|
|
double lower[1] = {0}; // Integration lower limit
|
|
double upper[1] = {1}; // Integration upper limit
|
|
|
|
// Initialize an RNG.
|
|
gsl_rng_env_setup();
|
|
gsl_rng *r = gsl_rng_alloc(gsl_rng_default);
|
|
|
|
// Define GSL function.
|
|
gsl_monte_function expo = {&function, dims, NULL};
|
|
expo.f = &function;
|
|
expo.dim = 1;
|
|
expo.params = NULL;
|
|
|
|
{
|
|
gsl_monte_vegas_state *s = gsl_monte_vegas_alloc (dims);
|
|
gsl_monte_vegas_integrate (&expo, lower, upper, dims, calls,
|
|
r, s, &integral, &error);
|
|
do { gsl_monte_vegas_integrate (&expo, lower, upper, dims,
|
|
calls/5, r, s, &integral, &error);
|
|
} while (fabs (gsl_monte_vegas_chisq (s) - 1.0) > 0.5);
|
|
|
|
// Free memory.
|
|
gsl_monte_vegas_free(s);
|
|
|
|
// Display results
|
|
results("Vegas", integral, error);
|
|
}
|
|
|
|
{
|
|
gsl_monte_miser_state *s = gsl_monte_miser_alloc (dims);
|
|
gsl_monte_miser_integrate (&expo, lower, upper, dims, calls,
|
|
r, s, &integral, &error);
|
|
// Free memory.
|
|
gsl_monte_miser_free (s);
|
|
|
|
// Display results
|
|
results ("Miser", integral, error);
|
|
}
|
|
|
|
{
|
|
gsl_monte_plain_state *s = gsl_monte_plain_alloc (dims);
|
|
gsl_monte_plain_integrate (&expo, lower, upper, dims, calls,
|
|
r, s, &integral, &error);
|
|
// Display results
|
|
results ("Montecarlo", integral, error);
|
|
|
|
// Free memory.
|
|
gsl_monte_plain_free (s);
|
|
}
|
|
|
|
// Free memory.
|
|
gsl_rng_free(r);
|
|
|
|
return EXIT_SUCCESS;
|
|
}
|