1 files changed, 351 insertions, 0 deletions
diff --git a/gsl-1.9/cdf/gauss.c b/gsl-1.9/cdf/gauss.c
new file mode 100644
index 0000000..235c978
--- /dev/null
+++ b/gsl-1.9/cdf/gauss.c
@@ -0,0 +1,351 @@
+/* cdf/gauss.c
+ *
+ * Copyright (C) 2002, 2004 Jason H. Stover.
+ *
+ * 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307, USA.
+ */
+
+/*
+ * Computes the cumulative distribution function for the Gaussian
+ * distribution using a rational function approximation.  The
+ * computation is for the standard Normal distribution, i.e., mean 0
+ * and standard deviation 1. If you want to compute Pr(X < t) for a
+ * Gaussian random variable X with non-zero mean m and standard
+ * deviation sd not equal to 1, find gsl_cdf_ugaussian ((t-m)/sd).
+ * This approximation is accurate to at least double precision. The
+ * accuracy was verified with a pari-gp script.  The largest error
+ * found was about 1.4E-20. The coefficients were derived by Cody.
+ *
+ * References:
+ *
+ * W.J. Cody. "Rational Chebyshev Approximations for the Error
+ * Function," Mathematics of Computation, v23 n107 1969, 631-637.
+ *
+ * W. Fraser, J.F Hart. "On the Computation of Rational Approximations
+ * to Continuous Functions," Communications of the ACM, v5 1962.
+ *
+ * W.J. Kennedy Jr., J.E. Gentle. "Statistical Computing." Marcel Dekker. 1980.
+ * 
+ *  
+ */
+
+#include <config.h>
+#include <math.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_cdf.h>
+
+#ifndef M_1_SQRT2PI
+#define M_1_SQRT2PI (M_2_SQRTPI * M_SQRT1_2 / 2.0)
+#endif
+
+#define SQRT32 (4.0 * M_SQRT2)
+
+/*
+ * IEEE double precision dependent constants.
+ *
+ * GAUSS_EPSILON: Smallest positive value such that 
+ *                gsl_cdf_gaussian(x) > 0.5.
+ * GAUSS_XUPPER: Largest value x such that gsl_cdf_gaussian(x) < 1.0.
+ * GAUSS_XLOWER: Smallest value x such that gsl_cdf_gaussian(x) > 0.0.
+ */
+
+#define GAUSS_EPSILON  (GSL_DBL_EPSILON / 2)
+#define GAUSS_XUPPER (8.572)
+#define GAUSS_XLOWER (-37.519)
+
+#define GAUSS_SCALE (16.0)
+
+static double
+get_del (double x, double rational)
+{
+  double xsq = 0.0;
+  double del = 0.0;
+  double result = 0.0;
+
+  xsq = floor (x * GAUSS_SCALE) / GAUSS_SCALE;
+  del = (x - xsq) * (x + xsq);
+  del *= 0.5;
+
+  result = exp (-0.5 * xsq * xsq) * exp (-1.0 * del) * rational;
+
+  return result;
+}
+
+/*
+ * Normal cdf for fabs(x) < 0.66291
+ */
+static double
+gauss_small (const double x)
+{
+  unsigned int i;
+  double result = 0.0;
+  double xsq;
+  double xnum;
+  double xden;
+
+  const double a[5] = {
+    2.2352520354606839287,
+    161.02823106855587881,
+    1067.6894854603709582,
+    18154.981253343561249,
+    0.065682337918207449113
+  };
+  const double b[4] = {
+    47.20258190468824187,
+    976.09855173777669322,
+    10260.932208618978205,
+    45507.789335026729956
+  };
+
+  xsq = x * x;
+  xnum = a[4] * xsq;
+  xden = xsq;
+
+  for (i = 0; i < 3; i++)
+    {
+      xnum = (xnum + a[i]) * xsq;
+      xden = (xden + b[i]) * xsq;
+    }
+
+  result = x * (xnum + a[3]) / (xden + b[3]);
+
+  return result;
+}
+
+/*
+ * Normal cdf for 0.66291 < fabs(x) < sqrt(32).
+ */
+static double
+gauss_medium (const double x)
+{
+  unsigned int i;
+  double temp = 0.0;
+  double result = 0.0;
+  double xnum;
+  double xden;
+  double absx;
+
+  const double c[9] = {
+    0.39894151208813466764,
+    8.8831497943883759412,
+    93.506656132177855979,
+    597.27027639480026226,
+    2494.5375852903726711,
+    6848.1904505362823326,
+    11602.651437647350124,
+    9842.7148383839780218,
+    1.0765576773720192317e-8
+  };
+  const double d[8] = {
+    22.266688044328115691,
+    235.38790178262499861,
+    1519.377599407554805,
+    6485.558298266760755,
+    18615.571640885098091,
+    34900.952721145977266,
+    38912.003286093271411,
+    19685.429676859990727
+  };
+
+  absx = fabs (x);
+
+  xnum = c[8] * absx;
+  xden = absx;
+
+  for (i = 0; i < 7; i++)
+    {
+      xnum = (xnum + c[i]) * absx;
+      xden = (xden + d[i]) * absx;
+    }
+
+  temp = (xnum + c[7]) / (xden + d[7]);
+
+  result = get_del (x, temp);
+
+  return result;
+}
+
+/*
+ * Normal cdf for 
+ * {sqrt(32) < x < GAUSS_XUPPER} union { GAUSS_XLOWER < x < -sqrt(32) }.
+ */
+static double
+gauss_large (const double x)
+{
+  int i;
+  double result;
+  double xsq;
+  double temp;
+  double xnum;
+  double xden;
+  double absx;
+
+  const double p[6] = {
+    0.21589853405795699,
+    0.1274011611602473639,
+    0.022235277870649807,
+    0.001421619193227893466,
+    2.9112874951168792e-5,
+    0.02307344176494017303
+  };
+  const double q[5] = {
+    1.28426009614491121,
+    0.468238212480865118,
+    0.0659881378689285515,
+    0.00378239633202758244,
+    7.29751555083966205e-5
+  };
+
+  absx = fabs (x);
+  xsq = 1.0 / (x * x);
+  xnum = p[5] * xsq;
+  xden = xsq;
+
+  for (i = 0; i < 4; i++)
+    {
+      xnum = (xnum + p[i]) * xsq;
+      xden = (xden + q[i]) * xsq;
+    }
+
+  temp = xsq * (xnum + p[4]) / (xden + q[4]);
+  temp = (M_1_SQRT2PI - temp) / absx;
+
+  result = get_del (x, temp);
+
+  return result;
+}
+
+double
+gsl_cdf_ugaussian_P (const double x)
+{
+  double result;
+  double absx = fabs (x);
+
+  if (absx < GAUSS_EPSILON)
+    {
+      result = 0.5;
+      return result;
+    }
+  else if (absx < 0.66291)
+    {
+      result = 0.5 + gauss_small (x);
+      return result;
+    }
+  else if (absx < SQRT32)
+    {
+      result = gauss_medium (x);
+
+      if (x > 0.0)
+        {
+          result = 1.0 - result;
+        }
+
+      return result;
+    }
+  else if (x > GAUSS_XUPPER)
+    {
+      result = 1.0;
+      return result;
+    }
+  else if (x < GAUSS_XLOWER)
+    {
+      result = 0.0;
+      return result;
+    }
+  else
+    {
+      result = gauss_large (x);
+
+      if (x > 0.0)
+        {
+          result = 1.0 - result;
+        }
+    }
+
+  return result;
+}
+
+double
+gsl_cdf_ugaussian_Q (const double x)
+{
+  double result;
+  double absx = fabs (x);
+
+  if (absx < GAUSS_EPSILON)
+    {
+      result = 0.5;
+      return result;
+    }
+  else if (absx < 0.66291)
+    {
+      result = gauss_small (x);
+
+      if (x < 0.0)
+        {
+          result = fabs (result) + 0.5;
+        }
+      else
+        {
+          result = 0.5 - result;
+        }
+
+      return result;
+    }
+  else if (absx < SQRT32)
+    {
+      result = gauss_medium (x);
+
+      if (x < 0.0)
+        {
+          result = 1.0 - result;
+        }
+
+      return result;
+    }
+  else if (x > -(GAUSS_XLOWER))
+    {
+      result = 0.0;
+      return result;
+    }
+  else if (x < -(GAUSS_XUPPER))
+    {
+      result = 1.0;
+      return result;
+    }
+  else
+    {
+      result = gauss_large (x);
+
+      if (x < 0.0)
+        {
+          result = 1.0 - result;
+        }
+
+    }
+
+  return result;
+}
+
+double
+gsl_cdf_gaussian_P (const double x, const double sigma)
+{
+  return gsl_cdf_ugaussian_P (x / sigma);
+}
+
+double
+gsl_cdf_gaussian_Q (const double x, const double sigma)
+{
+  return gsl_cdf_ugaussian_Q (x / sigma);
+}