diff options
Diffstat (limited to 'libm/float/knf.c')
| -rw-r--r-- | libm/float/knf.c | 252 |
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); +} |