diff options
Diffstat (limited to 'libm/double/kn.c')
| -rw-r--r-- | libm/double/kn.c | 255 |
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); -} |