kn.c
Go to the documentation of this file.
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  * IEEE 0,30 90000 1.8e-8 3.0e-10
32  *
33  * Error is high only near the crossover point x = 9.55
34  * between the two expansions used.
35  */
36 
37 
38 /*
39  * Cephes Math Library Release 2.8: June, 2000
40  * Copyright 1984, 1987, 1988, 2000 by Stephen L. Moshier
41  */
42 
43 
44 /*
45  * Algorithm for Kn.
46  * n-1
47  * -n - (n-k-1)! 2 k
48  * K (x) = 0.5 (x/2) > -------- (-x /4)
49  * n - k!
50  * k=0
51  *
52  * inf. 2 k
53  * n n - (x /4)
54  * + (-1) 0.5(x/2) > {p(k+1) + p(n+k+1) - 2log(x/2)} ---------
55  * - k! (n+k)!
56  * k=0
57  *
58  * where p(m) is the psi function: p(1) = -EUL and
59  *
60  * m-1
61  * -
62  * p(m) = -EUL + > 1/k
63  * -
64  * k=1
65  *
66  * For large x,
67  * 2 2 2
68  * u-1 (u-1 )(u-3 )
69  * K (z) = sqrt(pi/2z) exp(-z) { 1 + ------- + ------------ + ...}
70  * v 1 2
71  * 1! (8z) 2! (8z)
72  * asymptotically, where
73  *
74  * 2
75  * u = 4 v .
76  *
77  */
78 
79 #include "mconf.h"
80 #include <float.h>
81 
82 #define EUL 5.772156649015328606065e-1
83 #define MAXFAC 31
84 extern double MACHEP, MAXLOG;
85 
86 double kn(int nn, double x)
87 {
88  double k, kf, nk1f, nkf, zn, t, s, z0, z;
89  double ans, fn, pn, pk, zmn, tlg, tox;
90  int i, n;
91 
92  if (nn < 0)
93  n = -nn;
94  else
95  n = nn;
96 
97  if (n > MAXFAC) {
98  overf:
100  return (INFINITY);
101  }
102 
103  if (x <= 0.0) {
104  if (x < 0.0) {
105  sf_error("kn", SF_ERROR_DOMAIN, NULL);
106  return NAN;
107  }
108  else {
110  return INFINITY;
111  }
112  }
113 
114 
115  if (x > 9.55)
116  goto asymp;
117 
118  ans = 0.0;
119  z0 = 0.25 * x * x;
120  fn = 1.0;
121  pn = 0.0;
122  zmn = 1.0;
123  tox = 2.0 / x;
124 
125  if (n > 0) {
126  /* compute factorial of n and psi(n) */
127  pn = -EUL;
128  k = 1.0;
129  for (i = 1; i < n; i++) {
130  pn += 1.0 / k;
131  k += 1.0;
132  fn *= k;
133  }
134 
135  zmn = tox;
136 
137  if (n == 1) {
138  ans = 1.0 / x;
139  }
140  else {
141  nk1f = fn / n;
142  kf = 1.0;
143  s = nk1f;
144  z = -z0;
145  zn = 1.0;
146  for (i = 1; i < n; i++) {
147  nk1f = nk1f / (n - i);
148  kf = kf * i;
149  zn *= z;
150  t = nk1f * zn / kf;
151  s += t;
152  if ((DBL_MAX - fabs(t)) < fabs(s))
153  goto overf;
154  if ((tox > 1.0) && ((DBL_MAX / tox) < zmn))
155  goto overf;
156  zmn *= tox;
157  }
158  s *= 0.5;
159  t = fabs(s);
160  if ((zmn > 1.0) && ((DBL_MAX / zmn) < t))
161  goto overf;
162  if ((t > 1.0) && ((DBL_MAX / t) < zmn))
163  goto overf;
164  ans = s * zmn;
165  }
166  }
167 
168 
169  tlg = 2.0 * log(0.5 * x);
170  pk = -EUL;
171  if (n == 0) {
172  pn = pk;
173  t = 1.0;
174  }
175  else {
176  pn = pn + 1.0 / n;
177  t = 1.0 / fn;
178  }
179  s = (pk + pn - tlg) * t;
180  k = 1.0;
181  do {
182  t *= z0 / (k * (k + n));
183  pk += 1.0 / k;
184  pn += 1.0 / (k + n);
185  s += (pk + pn - tlg) * t;
186  k += 1.0;
187  }
188  while (fabs(t / s) > MACHEP);
189 
190  s = 0.5 * s / zmn;
191  if (n & 1)
192  s = -s;
193  ans += s;
194 
195  return (ans);
196 
197 
198 
199  /* Asymptotic expansion for Kn(x) */
200  /* Converges to 1.4e-17 for x > 18.4 */
201 
202  asymp:
203 
204  if (x > MAXLOG) {
206  return (0.0);
207  }
208  k = n;
209  pn = 4.0 * k * k;
210  pk = 1.0;
211  z0 = 8.0 * x;
212  fn = 1.0;
213  t = 1.0;
214  s = t;
215  nkf = INFINITY;
216  i = 0;
217  do {
218  z = pn - pk * pk;
219  t = t * z / (fn * z0);
220  nk1f = fabs(t);
221  if ((i >= n) && (nk1f > nkf)) {
222  goto adone;
223  }
224  nkf = nk1f;
225  s += t;
226  fn += 1.0;
227  pk += 2.0;
228  i += 1;
229  }
230  while (fabs(t / s) > MACHEP);
231 
232  adone:
233  ans = exp(-x) * sqrt(M_PI / (2.0 * x)) * s;
234  return (ans);
235 }
kn
double kn(int nn, double x)
Definition: kn.c:86
MACHEP
double MACHEP
Definition: const.c:54
s
RealScalar s
Definition: level1_cplx_impl.h:126
MAXLOG
double MAXLOG
Definition: kn.c:84
fn
static double fn[10]
Definition: fresnl.c:103
x
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
Definition: gnuplot_common_settings.hh:12
SF_ERROR_OVERFLOW
@ SF_ERROR_OVERFLOW
Definition: sf_error.h:12
log
const EIGEN_DEVICE_FUNC LogReturnType log() const
Definition: ArrayCwiseUnaryOps.h:128
exp
const EIGEN_DEVICE_FUNC ExpReturnType exp() const
Definition: ArrayCwiseUnaryOps.h:97
boost::multiprecision::fabs
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:119
SF_ERROR_UNDERFLOW
@ SF_ERROR_UNDERFLOW
Definition: sf_error.h:11
n
int n
Definition: BiCGSTAB_simple.cpp:1
pybind_wrapper_test_script.z
z
Definition: pybind_wrapper_test_script.py:61
EUL
#define EUL
Definition: kn.c:82
mconf.h
MAXFAC
#define MAXFAC
Definition: kn.c:83
sf_error
void sf_error(const char *func_name, sf_error_t code, const char *fmt,...)
Definition: sf_error.c:41
SF_ERROR_SINGULAR
@ SF_ERROR_SINGULAR
Definition: sf_error.h:10
NULL
#define NULL
Definition: ccolamd.c:609
M_PI
#define M_PI
Definition: mconf.h:117
align_3::t
Point2 t(10, 10)
nn
idx_t * nn
Definition: include/metis.h:207
SF_ERROR_DOMAIN
@ SF_ERROR_DOMAIN
Definition: sf_error.h:16
ceres::sqrt
Jet< T, N > sqrt(const Jet< T, N > &f)
Definition: jet.h:418
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9


gtsam
Author(s):
autogenerated on Sat Nov 16 2024 04:02:38