1 // Original code derived from Boost and is distributed here
2 // under the Boost license (licenses/boost-license.txt)
3 // Copyright (c) 2006 Xiaogang Zhang
4 // Copyright (c) 2007, 2017 John Maddock
5 // Secondary code copyright by its author and is distributed here
6 // under the BSD-3 license (LICENSE.md)
7
8 #ifndef STAN_MATH_PRIM_FUN_LOG_MODIFIED_BESSEL_FIRST_KIND_HPP
9 #define STAN_MATH_PRIM_FUN_LOG_MODIFIED_BESSEL_FIRST_KIND_HPP
10
11 #include <stan/math/prim/meta.hpp>
12 #include <stan/math/prim/err.hpp>
13 #include <stan/math/prim/fun/constants.hpp>
14 #include <stan/math/prim/fun/inv.hpp>
15 #include <stan/math/prim/fun/is_inf.hpp>
16 #include <stan/math/prim/fun/lgamma.hpp>
17 #include <stan/math/prim/fun/log.hpp>
18 #include <stan/math/prim/fun/log_sum_exp.hpp>
19 #include <stan/math/prim/fun/log1p.hpp>
20 #include <stan/math/prim/fun/log1p_exp.hpp>
21 #include <stan/math/prim/fun/multiply_log.hpp>
22 #include <stan/math/prim/fun/sqrt.hpp>
23 #include <stan/math/prim/fun/square.hpp>
24 #include <stan/math/prim/functor/apply_scalar_binary.hpp>
25 #include <boost/math/tools/rational.hpp>
26 #include <cmath>
27
28 namespace stan {
29 namespace math {
30
31 /* Log of the modified Bessel function of the first kind,
32 * which is better known as the Bessel I function. See
33 * modified_bessel_first_kind.hpp for the function definition.
34 * The derivatives are known to be incorrect for v = 0 because a
35 * simple constant 0 is returned.
36 *
37 * @tparam T1 type of the order (v)
38 * @tparam T2 type of argument (z)
39 * @param v Order, can be a non-integer but must be at least -1
40 * @param z Real non-negative number
41 * @throws std::domain_error if either v or z is NaN, z is
42 * negative, or v is less than -1
43 * @return log of Bessel I function
44 */
45 template <typename T1, typename T2,
46 require_all_stan_scalar_t<T1, T2>* = nullptr>
log_modified_bessel_first_kind(const T1 v,const T2 z)47 inline return_type_t<T1, T2, double> log_modified_bessel_first_kind(
48 const T1 v, const T2 z) {
49 check_not_nan("log_modified_bessel_first_kind", "first argument (order)", v);
50 check_not_nan("log_modified_bessel_first_kind", "second argument (variable)",
51 z);
52 check_nonnegative("log_modified_bessel_first_kind",
53 "second argument (variable)", z);
54 check_greater_or_equal("log_modified_bessel_first_kind",
55 "first argument (order)", v, -1);
56
57 using boost::math::tools::evaluate_polynomial;
58 using std::log;
59 using std::pow;
60 using std::sqrt;
61
62 using T = return_type_t<T1, T2, double>;
63
64 if (z == 0) {
65 if (v == 0) {
66 return 0.0;
67 }
68 if (v > 0) {
69 return NEGATIVE_INFTY;
70 }
71 return INFTY;
72 }
73 if (is_inf(z)) {
74 return z;
75 }
76 if (v == 0) {
77 // WARNING: will not autodiff for v = 0 correctly
78 // modified from Boost's bessel_i0_imp in the double precision case,
79 // which refers to:
80 // Modified Bessel function of the first kind of order zero
81 // we use the approximating forms derived in:
82 // "Rational Approximations for the Modified Bessel Function of the
83 // First Kind -- I0(x) for Computations with Double Precision"
84 // by Pavel Holoborodko, see
85 // http://www.advanpix.com/2015/11/11/rational-approximations-for-the-modified-bessel-function-of-the-first-kind-i0-computations-double-precision
86 // The actual coefficients used are [Boost's] own, and extend
87 // Pavel's work to precisions other than double.
88
89 if (z < 7.75) {
90 // Bessel I0 over[10 ^ -16, 7.75]
91 // Max error in interpolated form : 3.042e-18
92 // Max Error found at double precision = Poly : 5.106609e-16
93 // Cheb : 5.239199e-16
94 static const double P[]
95 = {1.00000000000000000e+00, 2.49999999999999909e-01,
96 2.77777777777782257e-02, 1.73611111111023792e-03,
97 6.94444444453352521e-05, 1.92901234513219920e-06,
98 3.93675991102510739e-08, 6.15118672704439289e-10,
99 7.59407002058973446e-12, 7.59389793369836367e-14,
100 6.27767773636292611e-16, 4.34709704153272287e-18,
101 2.63417742690109154e-20, 1.13943037744822825e-22,
102 9.07926920085624812e-25};
103 return log1p_exp(multiply_log(2, z) - log(4.0)
104 + log(evaluate_polynomial(P, 0.25 * square(z))));
105 }
106 if (z < 500) {
107 // Max error in interpolated form : 1.685e-16
108 // Max Error found at double precision = Poly : 2.575063e-16
109 // Cheb : 2.247615e+00
110 static const double P[]
111 = {3.98942280401425088e-01, 4.98677850604961985e-02,
112 2.80506233928312623e-02, 2.92211225166047873e-02,
113 4.44207299493659561e-02, 1.30970574605856719e-01,
114 -3.35052280231727022e+00, 2.33025711583514727e+02,
115 -1.13366350697172355e+04, 4.24057674317867331e+05,
116 -1.23157028595698731e+07, 2.80231938155267516e+08,
117 -5.01883999713777929e+09, 7.08029243015109113e+10,
118 -7.84261082124811106e+11, 6.76825737854096565e+12,
119 -4.49034849696138065e+13, 2.24155239966958995e+14,
120 -8.13426467865659318e+14, 2.02391097391687777e+15,
121 -3.08675715295370878e+15, 2.17587543863819074e+15};
122 return z + log(evaluate_polynomial(P, inv(z))) - multiply_log(0.5, z);
123 }
124 // Max error in interpolated form : 2.437e-18
125 // Max Error found at double precision = Poly : 1.216719e-16
126 static const double P[] = {3.98942280401432905e-01, 4.98677850491434560e-02,
127 2.80506308916506102e-02, 2.92179096853915176e-02,
128 4.53371208762579442e-02};
129 return z + log(evaluate_polynomial(P, inv(z))) - multiply_log(0.5, z);
130 }
131 if (v == 1) { // WARNING: will not autodiff for v = 1 correctly
132 // modified from Boost's bessel_i1_imp in the double precision case
133 // see credits above in the v == 0 case
134 if (z < 7.75) {
135 // Bessel I0 over[10 ^ -16, 7.75]
136 // Max error in interpolated form: 5.639e-17
137 // Max Error found at double precision = Poly: 1.795559e-16
138
139 static const double P[]
140 = {8.333333333333333803e-02, 6.944444444444341983e-03,
141 3.472222222225921045e-04, 1.157407407354987232e-05,
142 2.755731926254790268e-07, 4.920949692800671435e-09,
143 6.834657311305621830e-11, 7.593969849687574339e-13,
144 6.904822652741917551e-15, 5.220157095351373194e-17,
145 3.410720494727771276e-19, 1.625212890947171108e-21,
146 1.332898928162290861e-23};
147 T a = square(z) * 0.25;
148 T Q[3] = {1, 0.5, evaluate_polynomial(P, a)};
149 return log(z) + log(evaluate_polynomial(Q, a)) - LOG_TWO;
150 }
151 if (z < 500) {
152 // Max error in interpolated form: 1.796e-16
153 // Max Error found at double precision = Poly: 2.898731e-16
154
155 static const double P[]
156 = {3.989422804014406054e-01, -1.496033551613111533e-01,
157 -4.675104253598537322e-02, -4.090895951581637791e-02,
158 -5.719036414430205390e-02, -1.528189554374492735e-01,
159 3.458284470977172076e+00, -2.426181371595021021e+02,
160 1.178785865993440669e+04, -4.404655582443487334e+05,
161 1.277677779341446497e+07, -2.903390398236656519e+08,
162 5.192386898222206474e+09, -7.313784438967834057e+10,
163 8.087824484994859552e+11, -6.967602516005787001e+12,
164 4.614040809616582764e+13, -2.298849639457172489e+14,
165 8.325554073334618015e+14, -2.067285045778906105e+15,
166 3.146401654361325073e+15, -2.213318202179221945e+15};
167 return z + log(evaluate_polynomial(P, inv(z))) - multiply_log(0.5, z);
168 }
169 // Max error in interpolated form: 1.320e-19
170 // Max Error found at double precision = Poly: 7.065357e-17
171 static const double P[]
172 = {3.989422804014314820e-01, -1.496033551467584157e-01,
173 -4.675105322571775911e-02, -4.090421597376992892e-02,
174 -5.843630344778927582e-02};
175 return z + log(evaluate_polynomial(P, inv(z))) - multiply_log(0.5, z);
176 }
177 if (z > 100) {
178 // Boost does something like this in asymptotic_bessel_i_large_x
179 T lim = pow((square(v) + 2.5) / (2 * z), 3) / 24;
180 if (lim < (EPSILON * 10)) {
181 T s = 1;
182 T mu = 4 * square(v);
183 T ex = 8 * z;
184 T num = mu - 1;
185 T denom = ex;
186 s -= num / denom;
187 num *= mu - 9;
188 denom *= ex * 2;
189 s += num / denom;
190 num *= mu - 25;
191 denom *= ex * 3;
192 s -= num / denom;
193 s = z - log(sqrt(z * TWO_PI)) + log(s);
194 return s;
195 }
196 }
197
198 return_type_t<T2> log_half_z = log(0.5 * z);
199 return_type_t<T1> lgam = v > -1 ? lgamma(v + 1.0) : 0;
200 T lcons = (2.0 + v) * log_half_z;
201 T out;
202 if (v > -1) {
203 out = log_sum_exp(v * log_half_z - lgam, lcons - lgamma(v + 2));
204 lgam += log1p(v);
205 } else {
206 out = lcons;
207 }
208
209 double m = 2;
210 double lfac = 0;
211 T old_out;
212 do {
213 lfac += log(m);
214 lgam += log(v + m);
215 lcons += 2 * log_half_z;
216 old_out = out;
217 out = log_sum_exp(out, lcons - lfac - lgam); // underflows eventually
218 m++;
219 } while (out > old_out || out < old_out);
220 return out;
221 }
222
223 /**
224 * Enables the vectorised application of the log_modified_bessel_first_kind
225 * function, when the first and/or second arguments are containers.
226 *
227 * @tparam T1 type of first input
228 * @tparam T2 type of second input
229 * @param a First input
230 * @param b Second input
231 * @return log_modified_bessel_first_kind function applied to the two inputs.
232 */
233 template <typename T1, typename T2, require_any_container_t<T1, T2>* = nullptr>
log_modified_bessel_first_kind(const T1 & a,const T2 & b)234 inline auto log_modified_bessel_first_kind(const T1& a, const T2& b) {
235 return apply_scalar_binary(a, b, [&](const auto& c, const auto& d) {
236 return log_modified_bessel_first_kind(c, d);
237 });
238 }
239
240 } // namespace math
241 } // namespace stan
242
243 #endif
244