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, 255 deletions

diff --git a/libm/double/kn.c b/libm/double/kn.c
deleted file mode 100644
index 72a1c1a53..000000000
--- a/libm/double/kn.c
+++ /dev/null
@@ -1,255 +0,0 @@
-/*							kn.c
- *
- *	Modified Bessel function, third kind, integer order
- *
- *
- *
- * SYNOPSIS:
- *
- * double x, y, kn();
- * int n;
- *
- * y = kn( n, x );
- *
- *
- *
- * DESCRIPTION:
- *
- * Returns modified Bessel function of the third kind
- * of order n of the argument.
- *
- * The range is partitioned into the two intervals [0,9.55] and
- * (9.55, infinity).  An ascending power series is used in the
- * low range, and an asymptotic expansion in the high range.
- *
- *
- *
- * ACCURACY:
- *
- *                      Relative error:
- * arithmetic   domain     # trials      peak         rms
- *    DEC       0,30         3000       1.3e-9      5.8e-11
- *    IEEE      0,30        90000       1.8e-8      3.0e-10
- *
- *  Error is high only near the crossover point x = 9.55
- * between the two expansions used.
- */
-
-
-/*
-Cephes Math Library Release 2.8:  June, 2000
-Copyright 1984, 1987, 1988, 2000 by Stephen L. Moshier
-*/
-
-
-/*
-Algorithm for Kn.
-                       n-1 
-                   -n   -  (n-k-1)!    2   k
-K (x)  =  0.5 (x/2)     >  -------- (-x /4)
- n                      -     k!
-                       k=0
-
-                    inf.                                   2   k
-       n         n   -                                   (x /4)
- + (-1)  0.5(x/2)    >  {p(k+1) + p(n+k+1) - 2log(x/2)} ---------
-                     -                                  k! (n+k)!
-                    k=0
-
-where  p(m) is the psi function: p(1) = -EUL and
-
-                      m-1
-                       -
-      p(m)  =  -EUL +  >  1/k
-                       -
-                      k=1
-
-For large x,
-                                         2        2     2
-                                      u-1     (u-1 )(u-3 )
-K (z)  =  sqrt(pi/2z) exp(-z) { 1 + ------- + ------------ + ...}
- v                                        1            2
-                                    1! (8z)     2! (8z)
-asymptotically, where
-
-           2
-    u = 4 v .
-
-*/
-
-#include <math.h>
-
-#define EUL 5.772156649015328606065e-1
-#define MAXFAC 31
-#ifdef ANSIPROT
-extern double fabs ( double );
-extern double exp ( double );
-extern double log ( double );
-extern double sqrt ( double );
-#else
-double fabs(), exp(), log(), sqrt();
-#endif
-extern double MACHEP, MAXNUM, MAXLOG, PI;
-
-double kn( nn, x )
-int nn;
-double x;
-{
-double k, kf, nk1f, nkf, zn, t, s, z0, z;
-double ans, fn, pn, pk, zmn, tlg, tox;
-int i, n;
-
-if( nn < 0 )
-	n = -nn;
-else
-	n = nn;
-
-if( n > MAXFAC )
-	{
-overf:
-	mtherr( "kn", OVERFLOW );
-	return( MAXNUM );
-	}
-
-if( x <= 0.0 )
-	{
-	if( x < 0.0 )
-		mtherr( "kn", DOMAIN );
-	else
-		mtherr( "kn", SING );
-	return( MAXNUM );
-	}
-
-
-if( x > 9.55 )
-	goto asymp;
-
-ans = 0.0;
-z0 = 0.25 * x * x;
-fn = 1.0;
-pn = 0.0;
-zmn = 1.0;
-tox = 2.0/x;
-
-if( n > 0 )
-	{
-	/* compute factorial of n and psi(n) */
-	pn = -EUL;
-	k = 1.0;
-	for( i=1; i<n; i++ )
-		{
-		pn += 1.0/k;
-		k += 1.0;
-		fn *= k;
-		}
-
-	zmn = tox;
-
-	if( n == 1 )
-		{
-		ans = 1.0/x;
-		}
-	else
-		{
-		nk1f = fn/n;
-		kf = 1.0;
-		s = nk1f;
-		z = -z0;
-		zn = 1.0;
-		for( i=1; i<n; i++ )
-			{
-			nk1f = nk1f/(n-i);
-			kf = kf * i;
-			zn *= z;
-			t = nk1f * zn / kf;
-			s += t;   
-			if( (MAXNUM - fabs(t)) < fabs(s) )
-				goto overf;
-			if( (tox > 1.0) && ((MAXNUM/tox) < zmn) )
-				goto overf;
-			zmn *= tox;
-			}
-		s *= 0.5;
-		t = fabs(s);
-		if( (zmn > 1.0) && ((MAXNUM/zmn) < t) )
-			goto overf;
-		if( (t > 1.0) && ((MAXNUM/t) < zmn) )
-			goto overf;
-		ans = s * zmn;
-		}
-	}
-
-
-tlg = 2.0 * log( 0.5 * x );
-pk = -EUL;
-if( n == 0 )
-	{
-	pn = pk;
-	t = 1.0;
-	}
-else
-	{
-	pn = pn + 1.0/n;
-	t = 1.0/fn;
-	}
-s = (pk+pn-tlg)*t;
-k = 1.0;
-do
-	{
-	t *= z0 / (k * (k+n));
-	pk += 1.0/k;
-	pn += 1.0/(k+n);
-	s += (pk+pn-tlg)*t;
-	k += 1.0;
-	}
-while( fabs(t/s) > MACHEP );
-
-s = 0.5 * s / zmn;
-if( n & 1 )
-	s = -s;
-ans += s;
-
-return(ans);
-
-
-
-/* Asymptotic expansion for Kn(x) */
-/* Converges to 1.4e-17 for x > 18.4 */
-
-asymp:
-
-if( x > MAXLOG )
-	{
-	mtherr( "kn", UNDERFLOW );
-	return(0.0);
-	}
-k = n;
-pn = 4.0 * k * k;
-pk = 1.0;
-z0 = 8.0 * x;
-fn = 1.0;
-t = 1.0;
-s = t;
-nkf = MAXNUM;
-i = 0;
-do
-	{
-	z = pn - pk * pk;
-	t = t * z /(fn * z0);
-	nk1f = fabs(t);
-	if( (i >= n) && (nk1f > nkf) )
-		{
-		goto adone;
-		}
-	nkf = nk1f;
-	s += t;
-	fn += 1.0;
-	pk += 2.0;
-	i += 1;
-	}
-while( fabs(t/s) > MACHEP );
-
-adone:
-ans = exp(-x) * sqrt( PI/(2.0*x) ) * s;
-return(ans);
-}