1 files changed, 151 insertions, 0 deletions
diff --git a/gsl-1.9/monte/plain.c b/gsl-1.9/monte/plain.c
new file mode 100644
index 0000000..cab3ae2
--- /dev/null
+++ b/gsl-1.9/monte/plain.c
@@ -0,0 +1,151 @@
+/* monte/plain.c
+ * 
+ * Copyright (C) 1996, 1997, 1998, 1999, 2000 Michael Booth
+ * 
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or (at
+ * your option) any later version.
+ * 
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+/* Plain Monte-Carlo. */
+
+/* Author: MJB */
+
+#include <config.h>
+#include <math.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_monte_plain.h>
+
+int
+gsl_monte_plain_integrate (const gsl_monte_function * f,
+                           const double xl[], const double xu[],
+                           const size_t dim,
+                           const size_t calls,
+                           gsl_rng * r,
+                           gsl_monte_plain_state * state,
+                           double *result, double *abserr)
+{
+  double vol, m = 0, q = 0;
+  double *x = state->x;
+  size_t n, i;
+
+  if (dim != state->dim)
+    {
+      GSL_ERROR ("number of dimensions must match allocated size", GSL_EINVAL);
+    }
+
+  for (i = 0; i < dim; i++)
+    {
+      if (xu[i] <= xl[i])
+        {
+          GSL_ERROR ("xu must be greater than xl", GSL_EINVAL);
+        }
+
+      if (xu[i] - xl[i] > GSL_DBL_MAX)
+        {
+          GSL_ERROR ("Range of integration is too large, please rescale",
+                     GSL_EINVAL);
+        }
+    }
+
+  /* Compute the volume of the region */
+
+  vol = 1;
+
+  for (i = 0; i < dim; i++)
+    {
+      vol *= xu[i] - xl[i];
+    }
+
+  for (n = 0; n < calls; n++)
+    {
+      /* Choose a random point in the integration region */
+
+      for (i = 0; i < dim; i++)
+        {
+          x[i] = xl[i] + gsl_rng_uniform_pos (r) * (xu[i] - xl[i]);
+        }
+
+      {
+        double fval = GSL_MONTE_FN_EVAL (f, x);
+
+        /* recurrence for mean and variance */
+
+        double d = fval - m;
+        m += d / (n + 1.0);
+        q += d * d * (n / (n + 1.0));
+      }
+    }
+
+  *result = vol * m;
+
+  if (calls < 2)
+    {
+      *abserr = GSL_POSINF;
+    }
+  else
+    {
+      *abserr = vol * sqrt (q / (calls * (calls - 1.0)));
+    }
+
+  return GSL_SUCCESS;
+}
+
+gsl_monte_plain_state *
+gsl_monte_plain_alloc (size_t dim)
+{
+  gsl_monte_plain_state *s =
+    (gsl_monte_plain_state *) malloc (sizeof (gsl_monte_plain_state));
+
+  if (s == 0)
+    {
+      GSL_ERROR_VAL ("failed to allocate space for state struct",
+                     GSL_ENOMEM, 0);
+    }
+
+  s->x = (double *) malloc (dim * sizeof (double));
+
+  if (s->x == 0)
+    {
+      free (s);
+      GSL_ERROR_VAL ("failed to allocate space for working vector",
+                     GSL_ENOMEM, 0);
+    }
+
+  s->dim = dim;
+
+  return s;
+}
+
+/* Set some default values and whatever */
+
+int
+gsl_monte_plain_init (gsl_monte_plain_state * s)
+{
+  size_t i;
+
+  for (i = 0; i < s->dim; i++)
+    {
+      s->x[i] = 0.0;
+    }
+
+  return GSL_SUCCESS;
+}
+
+void
+gsl_monte_plain_free (gsl_monte_plain_state * s)
+{
+  free (s->x);
+  free (s);
+}