Totally rework the math library, this time based on the MacOs X

math library (which is itself based on the math lib from FreeBSD).
 -Erik

1 files changed, 0 insertions, 243 deletions

diff --git a/libm/double/kolmogorov.c b/libm/double/kolmogorov.c
deleted file mode 100644
index 0d6fe92bd..000000000
--- a/libm/double/kolmogorov.c
+++ /dev/null
@@ -1,243 +0,0 @@
-
-/* 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