k1.c
Go to the documentation of this file.
1 /* k1.c
2  *
3  * Modified Bessel function, third kind, order one
4  *
5  *
6  *
7  * SYNOPSIS:
8  *
9  * double x, y, k1();
10  *
11  * y = k1( x );
12  *
13  *
14  *
15  * DESCRIPTION:
16  *
17  * Computes the modified Bessel function of the third kind
18  * of order one of the argument.
19  *
20  * The range is partitioned into the two intervals [0,2] and
21  * (2, infinity). Chebyshev polynomial expansions are employed
22  * in each interval.
23  *
24  *
25  *
26  * ACCURACY:
27  *
28  * Relative error:
29  * arithmetic domain # trials peak rms
30  * IEEE 0, 30 30000 1.2e-15 1.6e-16
31  *
32  * ERROR MESSAGES:
33  *
34  * message condition value returned
35  * k1 domain x <= 0 INFINITY
36  *
37  */
38  /* k1e.c
39  *
40  * Modified Bessel function, third kind, order one,
41  * exponentially scaled
42  *
43  *
44  *
45  * SYNOPSIS:
46  *
47  * double x, y, k1e();
48  *
49  * y = k1e( x );
50  *
51  *
52  *
53  * DESCRIPTION:
54  *
55  * Returns exponentially scaled modified Bessel function
56  * of the third kind of order one of the argument:
57  *
58  * k1e(x) = exp(x) * k1(x).
59  *
60  *
61  *
62  * ACCURACY:
63  *
64  * Relative error:
65  * arithmetic domain # trials peak rms
66  * IEEE 0, 30 30000 7.8e-16 1.2e-16
67  * See k1().
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 x(K1(x) - log(x/2) I1(x))
79  * in the interval [0,2].
80  *
81  * lim(x->0){ x(K1(x) - log(x/2) I1(x)) } = 1.
82  */
83 
84 static double A[] = {
85  -7.02386347938628759343E-18,
86  -2.42744985051936593393E-15,
87  -6.66690169419932900609E-13,
88  -1.41148839263352776110E-10,
89  -2.21338763073472585583E-8,
90  -2.43340614156596823496E-6,
91  -1.73028895751305206302E-4,
92  -6.97572385963986435018E-3,
93  -1.22611180822657148235E-1,
94  -3.53155960776544875667E-1,
95  1.52530022733894777053E0
96 };
97 
98 /* Chebyshev coefficients for exp(x) sqrt(x) K1(x)
99  * in the interval [2,infinity].
100  *
101  * lim(x->inf){ exp(x) sqrt(x) K1(x) } = sqrt(pi/2).
102  */
103 static double B[] = {
104  -5.75674448366501715755E-18,
105  1.79405087314755922667E-17,
106  -5.68946255844285935196E-17,
107  1.83809354436663880070E-16,
108  -6.05704724837331885336E-16,
109  2.03870316562433424052E-15,
110  -7.01983709041831346144E-15,
111  2.47715442448130437068E-14,
112  -8.97670518232499435011E-14,
113  3.34841966607842919884E-13,
114  -1.28917396095102890680E-12,
115  5.13963967348173025100E-12,
116  -2.12996783842756842877E-11,
117  9.21831518760500529508E-11,
118  -4.19035475934189648750E-10,
119  2.01504975519703286596E-9,
120  -1.03457624656780970260E-8,
121  5.74108412545004946722E-8,
122  -3.50196060308781257119E-7,
123  2.40648494783721712015E-6,
124  -1.93619797416608296024E-5,
125  1.95215518471351631108E-4,
126  -2.85781685962277938680E-3,
127  1.03923736576817238437E-1,
128  2.72062619048444266945E0
129 };
130 
131 extern double MINLOG;
132 
133 double k1(double x)
134 {
135  double y, z;
136 
137  if (x == 0.0) {
139  return INFINITY;
140  }
141  else if (x < 0.0) {
142  sf_error("k1", SF_ERROR_DOMAIN, NULL);
143  return NAN;
144  }
145  z = 0.5 * x;
146 
147  if (x <= 2.0) {
148  y = x * x - 2.0;
149  y = log(z) * i1(x) + chbevl(y, A, 11) / x;
150  return (y);
151  }
152 
153  return (exp(-x) * chbevl(8.0 / x - 2.0, B, 25) / sqrt(x));
154 }
155 
156 
157 
158 
159 double k1e(double x)
160 {
161  double y;
162 
163  if (x == 0.0) {
165  return INFINITY;
166  }
167  else if (x < 0.0) {
168  sf_error("k1e", SF_ERROR_DOMAIN, NULL);
169  return NAN;
170  }
171 
172  if (x <= 2.0) {
173  y = x * x - 2.0;
174  y = log(0.5 * x) * i1(x) + chbevl(y, A, 11) / x;
175  return (y * exp(x));
176  }
177 
178  return (chbevl(8.0 / x - 2.0, B, 25) / sqrt(x));
179 }
MINLOG
double MINLOG
Definition: const.c:60
k1e
double k1e(double x)
Definition: k1.c:159
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
A
Definition: test_numpy_dtypes.cpp:298
k1
double k1(double x)
Definition: k1.c:133
pybind_wrapper_test_script.z
z
Definition: pybind_wrapper_test_script.py:61
y
Scalar * y
Definition: level1_cplx_impl.h:124
mconf.h
i1
double i1(double x)
Definition: i1.c:150
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 Fri Nov 1 2024 03:33:00