k0.c
Go to the documentation of this file.
1 /* k0.c
2  *
3  * Modified Bessel function, third kind, order zero
4  *
5  *
6  *
7  * SYNOPSIS:
8  *
9  * double x, y, k0();
10  *
11  * y = k0( x );
12  *
13  *
14  *
15  * DESCRIPTION:
16  *
17  * Returns modified Bessel function of the third kind
18  * of order zero of the argument.
19  *
20  * The range is partitioned into the two intervals [0,8] and
21  * (8, infinity). Chebyshev polynomial expansions are employed
22  * in each interval.
23  *
24  *
25  *
26  * ACCURACY:
27  *
28  * Tested at 2000 random points between 0 and 8. Peak absolute
29  * error (relative when K0 > 1) was 1.46e-14; rms, 4.26e-15.
30  * Relative error:
31  * arithmetic domain # trials peak rms
32  * IEEE 0, 30 30000 1.2e-15 1.6e-16
33  *
34  * ERROR MESSAGES:
35  *
36  * message condition value returned
37  * K0 domain x <= 0 INFINITY
38  *
39  */
40  /* k0e()
41  *
42  * Modified Bessel function, third kind, order zero,
43  * exponentially scaled
44  *
45  *
46  *
47  * SYNOPSIS:
48  *
49  * double x, y, k0e();
50  *
51  * y = k0e( x );
52  *
53  *
54  *
55  * DESCRIPTION:
56  *
57  * Returns exponentially scaled modified Bessel function
58  * of the third kind of order zero of the argument.
59  *
60  *
61  *
62  * ACCURACY:
63  *
64  * Relative error:
65  * arithmetic domain # trials peak rms
66  * IEEE 0, 30 30000 1.4e-15 1.4e-16
67  * See k0().
68  *
69  */
70 
71 /*
72  * Cephes Math Library Release 2.8: June, 2000
73  * Copyright 1984, 1987, 2000 by Stephen L. Moshier
74  */
75 
76 #include "mconf.h"
77 
78 /* Chebyshev coefficients for K0(x) + log(x/2) I0(x)
79  * in the interval [0,2]. The odd order coefficients are all
80  * zero; only the even order coefficients are listed.
81  *
82  * lim(x->0){ K0(x) + log(x/2) I0(x) } = -EUL.
83  */
84 
85 static double A[] = {
86  1.37446543561352307156E-16,
87  4.25981614279661018399E-14,
88  1.03496952576338420167E-11,
89  1.90451637722020886025E-9,
90  2.53479107902614945675E-7,
91  2.28621210311945178607E-5,
92  1.26461541144692592338E-3,
93  3.59799365153615016266E-2,
94  3.44289899924628486886E-1,
95  -5.35327393233902768720E-1
96 };
97 
98 /* Chebyshev coefficients for exp(x) sqrt(x) K0(x)
99  * in the inverted interval [2,infinity].
100  *
101  * lim(x->inf){ exp(x) sqrt(x) K0(x) } = sqrt(pi/2).
102  */
103 static double B[] = {
104  5.30043377268626276149E-18,
105  -1.64758043015242134646E-17,
106  5.21039150503902756861E-17,
107  -1.67823109680541210385E-16,
108  5.51205597852431940784E-16,
109  -1.84859337734377901440E-15,
110  6.34007647740507060557E-15,
111  -2.22751332699166985548E-14,
112  8.03289077536357521100E-14,
113  -2.98009692317273043925E-13,
114  1.14034058820847496303E-12,
115  -4.51459788337394416547E-12,
116  1.85594911495471785253E-11,
117  -7.95748924447710747776E-11,
118  3.57739728140030116597E-10,
119  -1.69753450938905987466E-9,
120  8.57403401741422608519E-9,
121  -4.66048989768794782956E-8,
122  2.76681363944501510342E-7,
123  -1.83175552271911948767E-6,
124  1.39498137188764993662E-5,
125  -1.28495495816278026384E-4,
126  1.56988388573005337491E-3,
127  -3.14481013119645005427E-2,
128  2.44030308206595545468E0
129 };
130 
131 double k0(double x)
132 {
133  double y, z;
134 
135  if (x == 0.0) {
137  return INFINITY;
138  }
139  else if (x < 0.0) {
140  sf_error("k0", SF_ERROR_DOMAIN, NULL);
141  return NAN;
142  }
143 
144  if (x <= 2.0) {
145  y = x * x - 2.0;
146  y = chbevl(y, A, 10) - log(0.5 * x) * i0(x);
147  return (y);
148  }
149  z = 8.0 / x - 2.0;
150  y = exp(-x) * chbevl(z, B, 25) / sqrt(x);
151  return (y);
152 }
153 
154 
155 
156 
157 double k0e(double x)
158 {
159  double y;
160 
161  if (x == 0.0) {
163  return INFINITY;
164  }
165  else if (x < 0.0) {
166  sf_error("k0e", SF_ERROR_DOMAIN, NULL);
167  return NAN;
168  }
169 
170  if (x <= 2.0) {
171  y = x * x - 2.0;
172  y = chbevl(y, A, 10) - log(0.5 * x) * i0(x);
173  return (y * exp(x));
174  }
175 
176  y = chbevl(8.0 / x - 2.0, B, 25) / sqrt(x);
177  return (y);
178 }
k0e
double k0e(double x)
Definition: k0.c:157
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
B
Definition: test_numpy_dtypes.cpp:299
log
const EIGEN_DEVICE_FUNC LogReturnType log() const
Definition: ArrayCwiseUnaryOps.h:128
exp
const EIGEN_DEVICE_FUNC ExpReturnType exp() const
Definition: ArrayCwiseUnaryOps.h:97
chbevl
double chbevl(double x, double array[], int n)
Definition: chbevl.c:63
k0
double k0(double x)
Definition: k0.c:131
A
Definition: test_numpy_dtypes.cpp:298
pybind_wrapper_test_script.z
z
Definition: pybind_wrapper_test_script.py:61
y
Scalar * y
Definition: level1_cplx_impl.h:124
i0
double i0(double x)
Definition: i0.c:149
mconf.h
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
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


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