gtsam
3rdparty
cephes
cephes
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:
99
sf_error
(
"kn"
,
SF_ERROR_OVERFLOW
,
NULL
);
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
{
109
sf_error
(
"kn"
,
SF_ERROR_SINGULAR
,
NULL
);
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
) {
205
sf_error
(
"kn"
,
SF_ERROR_UNDERFLOW
,
NULL
);
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 Fri Nov 1 2024 03:33:01