1 /* kn.c
2 *
3 * Modified Bessel function, third kind, integer order
4 *
5 *
6 *
7 * SYNOPSIS:
8 *
9 * double x, y, kn();
10 * int n;
11 *
12 * y = kn( n, x );
13 *
14 *
15 *
16 * DESCRIPTION:
17 *
18 * Returns modified Bessel function of the third kind
19 * of order n of the argument.
20 *
21 * The range is partitioned into the two intervals [0,9.55] and
22 * (9.55, infinity). An ascending power series is used in the
23 * low range, and an asymptotic expansion in the high range.
24 *
25 *
26 *
27 * ACCURACY:
28 *
29 * Relative error:
30 * arithmetic domain # trials peak rms
31 * DEC 0,30 3000 1.3e-9 5.8e-11
32 * IEEE 0,30 90000 1.8e-8 3.0e-10
33 *
34 * Error is high only near the crossover point x = 9.55
35 * between the two expansions used.
36 */
37
38
39 /*
40 Cephes Math Library Release 2.1: November, 1988
41 Copyright 1984, 1987, 1988 by Stephen L. Moshier
42 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
43
44 */
45
46
47 /*
48 Algorithm for Kn.
49 n-1
50 -n - (n-k-1)! 2 k
51 K (x) = 0.5 (x/2) > -------- (-x /4)
52 n - k!
53 k=0
54
55 inf. 2 k
56 n n - (x /4)
57 + (-1) 0.5(x/2) > {p(k+1) + p(n+k+1) - 2log(x/2)} ---------
58 - k! (n+k)!
59 k=0
60
61 where p(m) is the psi function: p(1) = -EUL and
62
63 m-1
64 -
65 p(m) = -EUL + > 1/k
66 -
67 k=1
68
69 For large x,
70 2 2 2
71 u-1 (u-1 )(u-3 )
72 K (z) = sqrt(pi/2z) exp(-z) { 1 + ------- + ------------ + ...}
73 v 1 2
74 1! (8z) 2! (8z)
75 asymptotically, where
76
77 2
78 u = 4 v .
79
80 */
81
82 #include "mconf.h"
83 #include "cephes.h"
84
85 #define EUL 5.772156649015328606065e-1
86 #define MAXFAC 31
87
88 extern double MACHEP, MAXNUM, MAXLOG, PI;
89
kn(nn,x)90 double kn( nn, x )
91 int nn;
92 double x;
93 {
94 double k, kf, nk1f, nkf, zn, t, s, z0, z;
95 double ans, fn, pn, pk, zmn, tlg, tox;
96 int i, n;
97
98 if( nn < 0 )
99 n = -nn;
100 else
101 n = nn;
102
103 if( n > MAXFAC )
104 {
105 overf:
106 mtherr( "kn", OVERFLOW );
107 return( MAXNUM );
108 }
109
110 if( x <= 0.0 )
111 {
112 if( x < 0.0 )
113 mtherr( "kn", DOMAIN );
114 else
115 mtherr( "kn", SING );
116 return( MAXNUM );
117 }
118
119
120 if( x > 9.55 )
121 goto asymp;
122
123 ans = 0.0;
124 z0 = 0.25 * x * x;
125 fn = 1.0;
126 pn = 0.0;
127 zmn = 1.0;
128 tox = 2.0/x;
129
130 if( n > 0 )
131 {
132 /* compute factorial of n and psi(n) */
133 pn = -EUL;
134 k = 1.0;
135 for( i=1; i<n; i++ )
136 {
137 pn += 1.0/k;
138 k += 1.0;
139 fn *= k;
140 }
141
142 zmn = tox;
143
144 if( n == 1 )
145 {
146 ans = 1.0/x;
147 }
148 else
149 {
150 nk1f = fn/n;
151 kf = 1.0;
152 s = nk1f;
153 z = -z0;
154 zn = 1.0;
155 for( i=1; i<n; i++ )
156 {
157 nk1f = nk1f/(n-i);
158 kf = kf * i;
159 zn *= z;
160 t = nk1f * zn / kf;
161 s += t;
162 if( (MAXNUM - fabs(t)) < fabs(s) )
163 goto overf;
164 if( (tox > 1.0) && ((MAXNUM/tox) < zmn) )
165 goto overf;
166 zmn *= tox;
167 }
168 s *= 0.5;
169 t = fabs(s);
170 if( (zmn > 1.0) && ((MAXNUM/zmn) < t) )
171 goto overf;
172 if( (t > 1.0) && ((MAXNUM/t) < zmn) )
173 goto overf;
174 ans = s * zmn;
175 }
176 }
177
178
179 tlg = 2.0 * log( 0.5 * x );
180 pk = -EUL;
181 if( n == 0 )
182 {
183 pn = pk;
184 t = 1.0;
185 }
186 else
187 {
188 pn = pn + 1.0/n;
189 t = 1.0/fn;
190 }
191 s = (pk+pn-tlg)*t;
192 k = 1.0;
193 do
194 {
195 t *= z0 / (k * (k+n));
196 pk += 1.0/k;
197 pn += 1.0/(k+n);
198 s += (pk+pn-tlg)*t;
199 k += 1.0;
200 }
201 while( fabs(t/s) > MACHEP );
202
203 s = 0.5 * s / zmn;
204 if( n & 1 )
205 s = -s;
206 ans += s;
207
208 return(ans);
209
210
211
212 /* Asymptotic expansion for Kn(x) */
213 /* Converges to 1.4e-17 for x > 18.4 */
214
215 asymp:
216
217 if( x > MAXLOG )
218 {
219 mtherr( "kn", UNDERFLOW );
220 return(0.0);
221 }
222 k = n;
223 pn = 4.0 * k * k;
224 pk = 1.0;
225 z0 = 8.0 * x;
226 fn = 1.0;
227 t = 1.0;
228 s = t;
229 nkf = MAXNUM;
230 i = 0;
231 do
232 {
233 z = pn - pk * pk;
234 t = t * z /(fn * z0);
235 nk1f = fabs(t);
236 if( (i >= n) && (nk1f > nkf) )
237 {
238 goto adone;
239 }
240 nkf = nk1f;
241 s += t;
242 fn += 1.0;
243 pk += 2.0;
244 i += 1;
245 }
246 while( fabs(t/s) > MACHEP );
247
248 adone:
249 ans = exp(-x) * sqrt( PI/(2.0*x) ) * s;
250 return(ans);
251 }
252