summaryrefslogtreecommitdiff
path: root/libm/float/knf.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
committerEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
commit1077fa4d772832f77a677ce7fb7c2d513b959e3f (patch)
tree579bee13fb0b58d2800206366ec2caecbb15f3fc /libm/float/knf.c
parent22358dd7ce7bb49792204b698f01a6f69b9c8e08 (diff)
uClibc now has a math library. muahahahaha!
-Erik
Diffstat (limited to 'libm/float/knf.c')
-rw-r--r--libm/float/knf.c252
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);
+}