summaryrefslogtreecommitdiff
path: root/libm/double/kn.c
diff options
context:
space:
mode:
Diffstat (limited to 'libm/double/kn.c')
-rw-r--r--libm/double/kn.c255
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);
-}