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