uClibc now has a math library. muahahahaha!

 -Erik

1 files changed, 252 insertions, 0 deletions

diff --git a/libm/float/knf.c b/libm/float/knf.c
new file mode 100644
index 000000000..85e297390
--- /dev/null
+++ b/libm/float/knf.c
@@ -0,0 +1,252 @@
+/*							knf.c
+ *
+ *	Modified Bessel function, third kind, integer order
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, y, knf();
+ * int n;
+ *
+ * y = knf( 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:
+ *
+ *          Absolute error, relative when function > 1:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      0,30        10000       2.0e-4      3.8e-6
+ *
+ *  Error is high only near the crossover point x = 9.55
+ * between the two expansions used.
+ */
+
+
+/*
+Cephes Math Library Release 2.2:  June, 1992
+Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier
+Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+
+*/
+
+
+/*
+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
+extern float MACHEPF, MAXNUMF, MAXLOGF, PIF;
+
+#define fabsf(x) ( (x) < 0 ? -(x) : (x) )
+
+float expf(float), logf(float), sqrtf(float);
+
+float knf( int nnn, float xx )
+{
+float x, k, kf, nk1f, nkf, zn, t, s, z0, z;
+float ans, fn, pn, pk, zmn, tlg, tox;
+int i, n, nn;
+
+nn = nnn;
+x = xx;
+if( nn < 0 )
+	n = -nn;
+else
+	n = nn;
+
+if( n > MAXFAC )
+	{
+overf:
+	mtherr( "knf", OVERFLOW );
+	return( MAXNUMF );
+	}
+
+if( x <= 0.0 )
+	{
+	if( x < 0.0 )
+		mtherr( "knf", DOMAIN );
+	else
+		mtherr( "knf", SING );
+	return( MAXNUMF );
+	}
+
+
+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( (MAXNUMF - fabsf(t)) < fabsf(s) )
+				goto overf;
+			if( (tox > 1.0) && ((MAXNUMF/tox) < zmn) )
+				goto overf;
+			zmn *= tox;
+			}
+		s *= 0.5;
+		t = fabsf(s);
+		if( (zmn > 1.0) && ((MAXNUMF/zmn) < t) )
+			goto overf;
+		if( (t > 1.0) && ((MAXNUMF/t) < zmn) )
+			goto overf;
+		ans = s * zmn;
+		}
+	}
+
+
+tlg = 2.0 * logf( 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( fabsf(t/s) > MACHEPF );
+
+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 > MAXLOGF )
+	{
+	mtherr( "knf", 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 = MAXNUMF;
+i = 0;
+do
+	{
+	z = pn - pk * pk;
+	t = t * z /(fn * z0);
+	nk1f = fabsf(t);
+	if( (i >= n) && (nk1f > nkf) )
+		{
+		goto adone;
+		}
+	nkf = nk1f;
+	s += t;
+	fn += 1.0;
+	pk += 2.0;
+	i += 1;
+	}
+while( fabsf(t/s) > MACHEPF );
+
+adone:
+ans = expf(-x) * sqrtf( PIF/(2.0*x) ) * s;
+return(ans);
+}