7 #include <boost/math/special_functions/bessel.hpp> 14 if (x < std::sqrt((
double)N+1) / 100)
18 return -std::log(std::tgamma(N+1));
20 return -std::log(std::tgamma(N+1)) + N * std::log(x / 2.0);
23 if (x <= std::max(9.23 * N + 15.934, 11.26 * N - 236.21))
24 return std::log(boost::math::cyl_bessel_i(N,x));
37 double resVal = x - 0.5 * std::log(2 * M_PI * x);
39 double secTermSerie = 1.0 / (8.0 * x);
40 double thirdTermSerie = 9.0 / (2.0 * (8.0 * x) * (8.0 * x));
42 resVal += std::log(1.0 + secTermSerie + thirdTermSerie);
44 for (
unsigned int i = 1;i <= N;++i)
53 return std::pow(x / 2.0, N) / std::tgamma(N+1);
55 double C = x / (N + 0.5 + std::sqrt(x * x + (N + 1.5) * (N + 1.5)));
56 double D = x / (N + 1.5 + std::sqrt(x * x + (N + 1.5) * (N + 1.5)));
58 double res = std::sqrt(2/x) * std::pow(N + 1,N + 0.5) / std::tgamma(N+1) * std::exp(x*D) * std::pow(C, N+0.5);
66 return N * std::log(x / 2.0) - std::log(std::tgamma(N+1));
68 double C = x / (N + 0.5 + std::sqrt(x * x + (N + 1.5) * (N + 1.5)));
69 double D = x / (N + 1.5 + std::sqrt(x * x + (N + 1.5) * (N + 1.5)));
71 double res = 0.5 * std::log(2/x) + (N + 0.5) * std::log((N + 1) * C) - std::log(std::tgamma(N+1)) + x * D;
82 if (x > std::max(70.0, 5.97 * N - 45.25))
84 double m1 = 4.0 * N * N;
85 double m0 = 4.0 * (N - 1.0) * (N - 1.0);
86 return 1.0 - (m1 - m0) / (8.0 * x) + ((m0 - 1.0) * (m0 + 7.0) - 2.0 * (m1 - 1.0) * (m0 - 1.0) + (m1 - 1.0) * (m1 - 9)) / (2.0 * (8.0 * x) * (8.0 * x));
92 for (
unsigned int i = 1;i <= approx_order;++i)
95 rho = - ak * (rho + 1.0) / (1.0 + ak * (1.0 + rho));
105 double res = x / (N - 0.5 + std::sqrt(x * x + (N + 0.5) * (N + 0.5)));
113 double firstTerm = tmpVal *
bessel_ratio_i(x, N + 1, approx_order);
114 double secondTerm = tmpVal * tmpVal;
115 double thirdTerm = tmpVal / x;
117 return firstTerm - secondTerm + thirdTerm;
124 double secondTerm = tmpVal * tmpVal;
125 double thirdTerm = tmpVal / x;
127 return firstTerm - secondTerm + thirdTerm;
136 double resVal = log(y / 2.0);
140 for (
unsigned int i = 0;i < approx_order;++i)
142 double tmpVal = pow(y * y / 4.0, (
double)i);
143 tmpVal /= std::tgamma(i+1);
144 tmpVal /= std::tgamma(order+i+1);
150 resVal -= num / denom;
157 return x / (x + 2.0 * N);
163 return - x * (N+0.5) / (2.0 * (N + x * 0.5) * (N + x + 0.5));
165 return -x * (N + k - 0.5) / (2.0 * (N + x + (k - 1) * 0.5) * (N + x + k * 0.5));
double ak_support(double x, unsigned int N, unsigned int k)
Support function for besserl_ratio_i.
double bessel_ratio_i_derivative_approx(double x, unsigned int N)
Computes fast and accurate approximation of the derivative of the ratio of modified Bessel functions ...
double log_bessel_i_lower_bound(unsigned int N, double x)
Computes a lower bound of the log of modified Bessel function of the first kind: I_{N} (N >= 0) ...
double psi_function(unsigned int n, double emc)
double bessel_ratio_i(double x, unsigned int N, unsigned int approx_order)
Computes the ratio of modified Bessel functions of the first kind: I_{N} / I_{N-1} (N >= 1) ...
double bessel_ratio_i_lower_bound(double x, unsigned int N)
Computes a lower bound of the ratio of modified Bessel functions of the first kind: I_{N} / I_{N-1} (...
double log_bessel_i(unsigned int N, double x)
Computes the log of modified Bessel function of the first kind: I_{N} (N >= 0)
double a0r_support(double x, unsigned int N)
Support function for besserl_ratio_i.
double log_bessel_order_derivative_i(double x, unsigned int order, double emc, unsigned int approx_order)
Computes the derivative of the log of modified Bessel function of the first kind w.r.t. its order (emc is the Euler-Mascheroni constant)
double bessel_i_lower_bound(unsigned int N, double x)
Computes a lower bound of the modified Bessel function of the first kind: I_{N} (N >= 0) ...
double bessel_ratio_i_derivative(double x, unsigned int N, unsigned int approx_order)
Computes the derivative of the ratio of modified Bessel functions of the first kind: d/dx( I_{N}(x) /...