ANIMA  4.0
animaBesselFunctions.cxx
Go to the documentation of this file.
1 #include <cmath>
2 
3 #include "animaBesselFunctions.h"
4 #include <animaGammaFunctions.h>
5 #include <algorithm>
6 
7 #include <boost/math/special_functions/bessel.hpp>
8 
9 namespace anima
10 {
11 
12 double log_bessel_i(unsigned int N, double x)
13 {
14  if (x < std::sqrt((double)N+1) / 100)
15  {
16  // tgamma(N) = factorial(N-1)
17  if (N == 0)
18  return -std::log(std::tgamma(N+1));
19  else
20  return -std::log(std::tgamma(N+1)) + N * std::log(x / 2.0);
21  }
22 
23  if (x <= std::max(9.23 * N + 15.934, 11.26 * N - 236.21)) // before was 600; now garantees an absolute error of maximum 1e-4 between true function and approximation
24  return std::log(boost::math::cyl_bessel_i(N,x));
25 
26 // // Too big value, switching to approximation
27 // // Works less and less when orders go up but above that barrier we get non valid values
28 // double resVal = x - 0.5 * log(2 * M_PI * x);
29 //
30 // double secTermSerie = - (4.0 * order * order - 1.0) / (8.0 * x);
31 // double thirdTermSerie = (4.0 * order * order - 1.0) * (4.0 * order * order - 9.0) / (2.0 * (8.0 * x) * (8.0 * x));
32 //
33 // resVal += log(1.0 + secTermSerie + thirdTermSerie);
34 
35  // Correct the problem
36  // First, compute approx for I_0
37  double resVal = x - 0.5 * std::log(2 * M_PI * x);
38 
39  double secTermSerie = 1.0 / (8.0 * x);
40  double thirdTermSerie = 9.0 / (2.0 * (8.0 * x) * (8.0 * x));
41 
42  resVal += std::log(1.0 + secTermSerie + thirdTermSerie);
43  // Then compute log(I_L) as log(I_0) + sum_{i=1}^L log(bessel_ratio_i(x,L))
44  for (unsigned int i = 1;i <= N;++i)
45  resVal += std::log(anima::bessel_ratio_i(x, i));
46 
47  return resVal;
48 }
49 
50 double bessel_i_lower_bound(unsigned int N, double x)
51 {
52  if (x < 1e-2)
53  return std::pow(x / 2.0, N) / std::tgamma(N+1);
54 
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)));
57 
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);
59 
60  return res;
61 }
62 
63 double log_bessel_i_lower_bound(unsigned int N, double x)
64 {
65  if (x < 1e-2)
66  return N * std::log(x / 2.0) - std::log(std::tgamma(N+1));
67 
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)));
70 
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;
72 
73  return res;
74 }
75 
76 double bessel_ratio_i(double x, unsigned int N, unsigned int approx_order)
77 {
78  if (N == 0)
79  return 0;
80 
81  // Ajout valeur asymptotique pour x infini
82  if (x > std::max(70.0, 5.97 * N - 45.25)) // this garantees an absolute error of maximum 1e-4 between true function and approximation
83  {
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));
87  }
88 
89  double pk = 1;
90  double rho = 0;
91  double resVal = pk;
92  for (unsigned int i = 1;i <= approx_order;++i)
93  {
94  double ak = anima::ak_support(x,N,i);
95  rho = - ak * (rho + 1.0) / (1.0 + ak * (1.0 + rho));
96  pk = pk * rho;
97  resVal += pk;
98  }
99 
100  return anima::a0r_support(x,N) * resVal;
101 }
102 
103 double bessel_ratio_i_lower_bound(double x, unsigned int N)
104 {
105  double res = x / (N - 0.5 + std::sqrt(x * x + (N + 0.5) * (N + 0.5)));
106 
107  return res;
108 }
109 
110 double bessel_ratio_i_derivative(double x, unsigned int N, unsigned int approx_order)
111 {
112  double tmpVal = bessel_ratio_i(x, N, approx_order);
113  double firstTerm = tmpVal * bessel_ratio_i(x, N + 1, approx_order);
114  double secondTerm = tmpVal * tmpVal;
115  double thirdTerm = tmpVal / x;
116 
117  return firstTerm - secondTerm + thirdTerm;
118 }
119 
120 double bessel_ratio_i_derivative_approx(double x, unsigned int N)
121 {
122  double tmpVal = bessel_ratio_i_lower_bound(x, N);
123  double firstTerm = tmpVal * bessel_ratio_i_lower_bound(x, N + 1);
124  double secondTerm = tmpVal * tmpVal;
125  double thirdTerm = tmpVal / x;
126 
127  return firstTerm - secondTerm + thirdTerm;
128 }
129 
130 double log_bessel_order_derivative_i(double x, unsigned int order, double emc, unsigned int approx_order)
131 {
132  double y = x;
133  if (y > 2795.0)
134  y = 2795.0;
135 
136  double resVal = log(y / 2.0);
137 
138  double num = 0;
139  double denom = 0;
140  for (unsigned int i = 0;i < approx_order;++i)
141  {
142  double tmpVal = pow(y * y / 4.0, (double)i);
143  tmpVal /= std::tgamma(i+1);
144  tmpVal /= std::tgamma(order+i+1);
145  denom += tmpVal;
146  tmpVal *= anima::psi_function(order+i+1, emc);
147  num += tmpVal;
148  }
149 
150  resVal -= num / denom;
151 
152  return resVal;
153 }
154 
155 double a0r_support(double x, unsigned int N)
156 {
157  return x / (x + 2.0 * N);
158 }
159 
160 double ak_support(double x, unsigned int N, unsigned int k)
161 {
162  if (k == 1)
163  return - x * (N+0.5) / (2.0 * (N + x * 0.5) * (N + x + 0.5));
164  else
165  return -x * (N + k - 0.5) / (2.0 * (N + x + (k - 1) * 0.5) * (N + x + k * 0.5));
166 }
167 
168 } // end namespace anima
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) /...