uClibc now has a math library. muahahahaha!

 -Erik

1 files changed, 243 insertions, 0 deletions

diff --git a/libm/double/kolmogorov.c b/libm/double/kolmogorov.c
new file mode 100644
index 000000000..0d6fe92bd
--- /dev/null
+++ b/libm/double/kolmogorov.c
@@ -0,0 +1,243 @@
+
+/* Re Kolmogorov statistics, here is Birnbaum and Tingey's formula for the
+   distribution of D+, the maximum of all positive deviations between a
+   theoretical distribution function P(x) and an empirical one Sn(x)
+   from n samples.
+
+     +
+    D  =         sup     [P(x) - S (x)]
+     n     -inf < x < inf         n
+
+
+                  [n(1-e)]
+        +            -                    v-1              n-v
+    Pr{D   > e} =    >    C    e (e + v/n)    (1 - e - v/n)
+        n            -   n v
+                    v=0
+
+    [n(1-e)] is the largest integer not exceeding n(1-e).
+    nCv is the number of combinations of n things taken v at a time.  */
+
+
+#include <math.h>
+#ifdef ANSIPROT
+extern double pow ( double, double );
+extern double floor ( double );
+extern double lgam ( double );
+extern double exp ( double );
+extern double sqrt ( double );
+extern double log ( double );
+extern double fabs ( double );
+double smirnov ( int, double );
+double kolmogorov ( double );
+#else
+double pow (), floor (), lgam (), exp (), sqrt (), log (), fabs ();
+double smirnov (), kolmogorov ();
+#endif
+extern double MAXLOG;
+
+/* Exact Smirnov statistic, for one-sided test.  */
+double
+smirnov (n, e)
+     int n;
+     double e;
+{
+  int v, nn;
+  double evn, omevn, p, t, c, lgamnp1;
+
+  if (n <= 0 || e < 0.0 || e > 1.0)
+    return (-1.0);
+  nn = floor ((double) n * (1.0 - e));
+  p = 0.0;
+  if (n < 1013)
+    {
+      c = 1.0;
+      for (v = 0; v <= nn; v++)
+	{
+	  evn = e + ((double) v) / n;
+	  p += c * pow (evn, (double) (v - 1))
+	    * pow (1.0 - evn, (double) (n - v));
+	  /* Next combinatorial term; worst case error = 4e-15.  */
+	  c *= ((double) (n - v)) / (v + 1);
+	}
+    }
+  else
+    {
+      lgamnp1 = lgam ((double) (n + 1));
+      for (v = 0; v <= nn; v++)
+	{
+	  evn = e + ((double) v) / n;
+	  omevn = 1.0 - evn;
+	  if (fabs (omevn) > 0.0)
+	    {
+	      t = lgamnp1
+		- lgam ((double) (v + 1))
+		- lgam ((double) (n - v + 1))
+		+ (v - 1) * log (evn)
+		+ (n - v) * log (omevn);
+	      if (t > -MAXLOG)
+		p += exp (t);
+	    }
+	}
+    }
+  return (p * e);
+}
+
+
+/* Kolmogorov's limiting distribution of two-sided test, returns
+   probability that sqrt(n) * max deviation > y,
+   or that max deviation > y/sqrt(n).
+   The approximation is useful for the tail of the distribution
+   when n is large.  */
+double
+kolmogorov (y)
+     double y;
+{
+  double p, t, r, sign, x;
+
+  x = -2.0 * y * y;
+  sign = 1.0;
+  p = 0.0;
+  r = 1.0;
+  do
+    {
+      t = exp (x * r * r);
+      p += sign * t;
+      if (t == 0.0)
+	break;
+      r += 1.0;
+      sign = -sign;
+    }
+  while ((t / p) > 1.1e-16);
+  return (p + p);
+}
+
+/* Functional inverse of Smirnov distribution
+   finds e such that smirnov(n,e) = p.  */
+double
+smirnovi (n, p)
+     int n;
+     double p;
+{
+  double e, t, dpde;
+
+  if (p <= 0.0 || p > 1.0)
+    {
+      mtherr ("smirnovi", DOMAIN);
+      return 0.0;
+    }
+  /* Start with approximation p = exp(-2 n e^2).  */
+  e = sqrt (-log (p) / (2.0 * n));
+  do
+    {
+      /* Use approximate derivative in Newton iteration. */
+      t = -2.0 * n * e;
+      dpde = 2.0 * t * exp (t * e);
+      if (fabs (dpde) > 0.0)
+	t = (p - smirnov (n, e)) / dpde;
+      else
+	{
+	  mtherr ("smirnovi", UNDERFLOW);
+	  return 0.0;
+	}
+      e = e + t;
+      if (e >= 1.0 || e <= 0.0)
+	{
+	  mtherr ("smirnovi", OVERFLOW);
+	  return 0.0;
+	}
+    }
+  while (fabs (t / e) > 1e-10);
+  return (e);
+}
+
+
+/* Functional inverse of Kolmogorov statistic for two-sided test.
+   Finds y such that kolmogorov(y) = p.
+   If e = smirnovi (n,p), then kolmogi(2 * p) / sqrt(n) should
+   be close to e.  */
+double
+kolmogi (p)
+     double p;
+{
+  double y, t, dpdy;
+
+  if (p <= 0.0 || p > 1.0)
+    {
+      mtherr ("kolmogi", DOMAIN);
+      return 0.0;
+    }
+  /* Start with approximation p = 2 exp(-2 y^2).  */
+  y = sqrt (-0.5 * log (0.5 * p));
+  do
+    {
+      /* Use approximate derivative in Newton iteration. */
+      t = -2.0 * y;
+      dpdy = 4.0 * t * exp (t * y);
+      if (fabs (dpdy) > 0.0)
+	t = (p - kolmogorov (y)) / dpdy;
+      else
+	{
+	  mtherr ("kolmogi", UNDERFLOW);
+	  return 0.0;
+	}
+      y = y + t;
+    }
+  while (fabs (t / y) > 1e-10);
+  return (y);
+}
+
+
+#ifdef SALONE
+/* Type in a number.  */
+void
+getnum (s, px)
+     char *s;
+     double *px;
+{
+  char str[30];
+
+  printf (" %s (%.15e) ? ", s, *px);
+  gets (str);
+  if (str[0] == '\0' || str[0] == '\n')
+    return;
+  sscanf (str, "%lf", px);
+  printf ("%.15e\n", *px);
+}
+
+/* Type in values, get answers.  */
+void
+main ()
+{
+  int n;
+  double e, p, ps, pk, ek, y;
+
+  n = 5;
+  e = 0.0;
+  p = 0.1;
+loop:
+  ps = n;
+  getnum ("n", &ps);
+  n = ps;
+  if (n <= 0)
+    {
+      printf ("? Operator error.\n");
+      goto loop;
+    }
+  /*
+  getnum ("e", &e);
+  ps = smirnov (n, e);
+  y = sqrt ((double) n) * e;
+  printf ("y = %.4e\n", y);
+  pk = kolmogorov (y);
+  printf ("Smirnov = %.15e, Kolmogorov/2 = %.15e\n", ps, pk / 2.0);
+*/
+  getnum ("p", &p);
+  e = smirnovi (n, p);
+  printf ("Smirnov e = %.15e\n", e);
+  y = kolmogi (2.0 * p);
+  ek = y / sqrt ((double) n);
+  printf ("Kolmogorov e = %.15e\n", ek);
+  goto loop;
+}
+#endif