diff options
Diffstat (limited to 'libm/float/knf.c')
| -rw-r--r-- | libm/float/knf.c | 252 |
1 files changed, 0 insertions, 252 deletions
diff --git a/libm/float/knf.c b/libm/float/knf.c deleted file mode 100644 index 85e297390..000000000 --- a/libm/float/knf.c +++ /dev/null @@ -1,252 +0,0 @@ -/* 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); -} |