1 | /*! |
---|
2 | * \file |
---|
3 | * \brief Implementation of modified Bessel functions of order zero |
---|
4 | * \author Tony Ottosson |
---|
5 | * |
---|
6 | * ------------------------------------------------------------------------- |
---|
7 | * |
---|
8 | * IT++ - C++ library of mathematical, signal processing, speech processing, |
---|
9 | * and communications classes and functions |
---|
10 | * |
---|
11 | * Copyright (C) 1995-2009 (see AUTHORS file for a list of contributors) |
---|
12 | * |
---|
13 | * This program is free software; you can redistribute it and/or modify |
---|
14 | * it under the terms of the GNU General Public License as published by |
---|
15 | * the Free Software Foundation; either version 2 of the License, or |
---|
16 | * (at your option) any later version. |
---|
17 | * |
---|
18 | * This program is distributed in the hope that it will be useful, |
---|
19 | * but WITHOUT ANY WARRANTY; without even the implied warranty of |
---|
20 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
---|
21 | * GNU General Public License for more details. |
---|
22 | * |
---|
23 | * You should have received a copy of the GNU General Public License |
---|
24 | * along with this program; if not, write to the Free Software |
---|
25 | * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA |
---|
26 | * |
---|
27 | * ------------------------------------------------------------------------- |
---|
28 | * |
---|
29 | * This is slightly modified routine from the Cephes library: |
---|
30 | * http://www.netlib.org/cephes/ |
---|
31 | */ |
---|
32 | |
---|
33 | #include <itpp/base/bessel/bessel_internal.h> |
---|
34 | |
---|
35 | |
---|
36 | /* |
---|
37 | * |
---|
38 | * Modified Bessel function of order zero |
---|
39 | * |
---|
40 | * double x, y, i0(); |
---|
41 | * |
---|
42 | * y = i0( x ); |
---|
43 | * |
---|
44 | * |
---|
45 | * DESCRIPTION: |
---|
46 | * |
---|
47 | * Returns modified Bessel function of order zero of the |
---|
48 | * argument. |
---|
49 | * |
---|
50 | * The function is defined as i0(x) = j0( ix ). |
---|
51 | * |
---|
52 | * The range is partitioned into the two intervals [0,8] and |
---|
53 | * (8, infinity). Chebyshev polynomial expansions are employed |
---|
54 | * in each interval. |
---|
55 | * |
---|
56 | * |
---|
57 | * ACCURACY: |
---|
58 | * |
---|
59 | * Relative error: |
---|
60 | * arithmetic domain # trials peak rms |
---|
61 | * IEEE 0,30 30000 5.8e-16 1.4e-16 |
---|
62 | * |
---|
63 | */ |
---|
64 | |
---|
65 | /* |
---|
66 | * Modified Bessel function of order zero, |
---|
67 | * exponentially scaled |
---|
68 | * |
---|
69 | * double x, y, i0e(); |
---|
70 | * |
---|
71 | * y = i0e( x ); |
---|
72 | * |
---|
73 | * DESCRIPTION: |
---|
74 | * |
---|
75 | * Returns exponentially scaled modified Bessel function |
---|
76 | * of order zero of the argument. |
---|
77 | * |
---|
78 | * The function is defined as i0e(x) = exp(-|x|) j0( ix ). |
---|
79 | * |
---|
80 | * ACCURACY: |
---|
81 | * |
---|
82 | * Relative error: |
---|
83 | * arithmetic domain # trials peak rms |
---|
84 | * IEEE 0,30 30000 5.4e-16 1.2e-16 |
---|
85 | * See i0(). |
---|
86 | */ |
---|
87 | |
---|
88 | /* |
---|
89 | Cephes Math Library Release 2.8: June, 2000 |
---|
90 | Copyright 1984, 1987, 2000 by Stephen L. Moshier |
---|
91 | */ |
---|
92 | |
---|
93 | |
---|
94 | /* Chebyshev coefficients for exp(-x) I0(x) |
---|
95 | * in the interval [0,8]. |
---|
96 | * |
---|
97 | * lim(x->0){ exp(-x) I0(x) } = 1. |
---|
98 | */ |
---|
99 | |
---|
100 | static double A[] = { |
---|
101 | -4.41534164647933937950E-18, |
---|
102 | 3.33079451882223809783E-17, |
---|
103 | -2.43127984654795469359E-16, |
---|
104 | 1.71539128555513303061E-15, |
---|
105 | -1.16853328779934516808E-14, |
---|
106 | 7.67618549860493561688E-14, |
---|
107 | -4.85644678311192946090E-13, |
---|
108 | 2.95505266312963983461E-12, |
---|
109 | -1.72682629144155570723E-11, |
---|
110 | 9.67580903537323691224E-11, |
---|
111 | -5.18979560163526290666E-10, |
---|
112 | 2.65982372468238665035E-9, |
---|
113 | -1.30002500998624804212E-8, |
---|
114 | 6.04699502254191894932E-8, |
---|
115 | -2.67079385394061173391E-7, |
---|
116 | 1.11738753912010371815E-6, |
---|
117 | -4.41673835845875056359E-6, |
---|
118 | 1.64484480707288970893E-5, |
---|
119 | -5.75419501008210370398E-5, |
---|
120 | 1.88502885095841655729E-4, |
---|
121 | -5.76375574538582365885E-4, |
---|
122 | 1.63947561694133579842E-3, |
---|
123 | -4.32430999505057594430E-3, |
---|
124 | 1.05464603945949983183E-2, |
---|
125 | -2.37374148058994688156E-2, |
---|
126 | 4.93052842396707084878E-2, |
---|
127 | -9.49010970480476444210E-2, |
---|
128 | 1.71620901522208775349E-1, |
---|
129 | -3.04682672343198398683E-1, |
---|
130 | 6.76795274409476084995E-1 |
---|
131 | }; |
---|
132 | |
---|
133 | |
---|
134 | |
---|
135 | /* Chebyshev coefficients for exp(-x) sqrt(x) I0(x) |
---|
136 | * in the inverted interval [8,infinity]. |
---|
137 | * |
---|
138 | * lim(x->inf){ exp(-x) sqrt(x) I0(x) } = 1/sqrt(2pi). |
---|
139 | */ |
---|
140 | |
---|
141 | static double B[] = { |
---|
142 | -7.23318048787475395456E-18, |
---|
143 | -4.83050448594418207126E-18, |
---|
144 | 4.46562142029675999901E-17, |
---|
145 | 3.46122286769746109310E-17, |
---|
146 | -2.82762398051658348494E-16, |
---|
147 | -3.42548561967721913462E-16, |
---|
148 | 1.77256013305652638360E-15, |
---|
149 | 3.81168066935262242075E-15, |
---|
150 | -9.55484669882830764870E-15, |
---|
151 | -4.15056934728722208663E-14, |
---|
152 | 1.54008621752140982691E-14, |
---|
153 | 3.85277838274214270114E-13, |
---|
154 | 7.18012445138366623367E-13, |
---|
155 | -1.79417853150680611778E-12, |
---|
156 | -1.32158118404477131188E-11, |
---|
157 | -3.14991652796324136454E-11, |
---|
158 | 1.18891471078464383424E-11, |
---|
159 | 4.94060238822496958910E-10, |
---|
160 | 3.39623202570838634515E-9, |
---|
161 | 2.26666899049817806459E-8, |
---|
162 | 2.04891858946906374183E-7, |
---|
163 | 2.89137052083475648297E-6, |
---|
164 | 6.88975834691682398426E-5, |
---|
165 | 3.36911647825569408990E-3, |
---|
166 | 8.04490411014108831608E-1 |
---|
167 | }; |
---|
168 | |
---|
169 | |
---|
170 | double i0(double x) |
---|
171 | { |
---|
172 | double y; |
---|
173 | |
---|
174 | if (x < 0) |
---|
175 | x = -x; |
---|
176 | if (x <= 8.0) { |
---|
177 | y = (x / 2.0) - 2.0; |
---|
178 | return(exp(x) * chbevl(y, A, 30)); |
---|
179 | } |
---|
180 | |
---|
181 | return(exp(x) * chbevl(32.0 / x - 2.0, B, 25) / sqrt(x)); |
---|
182 | |
---|
183 | } |
---|
184 | |
---|
185 | |
---|
186 | double i0e(double x) |
---|
187 | { |
---|
188 | double y; |
---|
189 | |
---|
190 | if (x < 0) |
---|
191 | x = -x; |
---|
192 | if (x <= 8.0) { |
---|
193 | y = (x / 2.0) - 2.0; |
---|
194 | return(chbevl(y, A, 30)); |
---|
195 | } |
---|
196 | |
---|
197 | return(chbevl(32.0 / x - 2.0, B, 25) / sqrt(x)); |
---|
198 | |
---|
199 | } |
---|