igami.c
Go to the documentation of this file.
1 /*
2  * (C) Copyright John Maddock 2006.
3  * Use, modification and distribution are subject to the
4  * Boost Software License, Version 1.0. (See accompanying file
5  * LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt)
6  */
7 #include "mconf.h"
8 
9 static double find_inverse_s(double, double);
10 static double didonato_SN(double, double, unsigned, double);
11 static double find_inverse_gamma(double, double, double);
12 
13 
14 static double find_inverse_s(double p, double q)
15 {
16  /*
17  * Computation of the Incomplete Gamma Function Ratios and their Inverse
18  * ARMIDO R. DIDONATO and ALFRED H. MORRIS, JR.
19  * ACM Transactions on Mathematical Software, Vol. 12, No. 4,
20  * December 1986, Pages 377-393.
21  *
22  * See equation 32.
23  */
24  double s, t;
25  double a[4] = {0.213623493715853, 4.28342155967104,
26  11.6616720288968, 3.31125922108741};
27  double b[5] = {0.3611708101884203e-1, 1.27364489782223,
28  6.40691597760039, 6.61053765625462, 1};
29 
30  if (p < 0.5) {
31  t = sqrt(-2 * log(p));
32  }
33  else {
34  t = sqrt(-2 * log(q));
35  }
36  s = t - polevl(t, a, 3) / polevl(t, b, 4);
37  if(p < 0.5)
38  s = -s;
39  return s;
40 }
41 
42 
43 static double didonato_SN(double a, double x, unsigned N, double tolerance)
44 {
45  /*
46  * Computation of the Incomplete Gamma Function Ratios and their Inverse
47  * ARMIDO R. DIDONATO and ALFRED H. MORRIS, JR.
48  * ACM Transactions on Mathematical Software, Vol. 12, No. 4,
49  * December 1986, Pages 377-393.
50  *
51  * See equation 34.
52  */
53  double sum = 1.0;
54 
55  if (N >= 1) {
56  unsigned i;
57  double partial = x / (a + 1);
58 
59  sum += partial;
60  for(i = 2; i <= N; ++i) {
61  partial *= x / (a + i);
62  sum += partial;
63  if(partial < tolerance) {
64  break;
65  }
66  }
67  }
68  return sum;
69 }
70 
71 
72 static double find_inverse_gamma(double a, double p, double q)
73 {
74  /*
75  * In order to understand what's going on here, you will
76  * need to refer to:
77  *
78  * Computation of the Incomplete Gamma Function Ratios and their Inverse
79  * ARMIDO R. DIDONATO and ALFRED H. MORRIS, JR.
80  * ACM Transactions on Mathematical Software, Vol. 12, No. 4,
81  * December 1986, Pages 377-393.
82  */
83  double result;
84 
85  if (a == 1) {
86  if (q > 0.9) {
87  result = -log1p(-p);
88  }
89  else {
90  result = -log(q);
91  }
92  }
93  else if (a < 1) {
94  double g = Gamma(a);
95  double b = q * g;
96 
97  if ((b > 0.6) || ((b >= 0.45) && (a >= 0.3))) {
98  /* DiDonato & Morris Eq 21:
99  *
100  * There is a slight variation from DiDonato and Morris here:
101  * the first form given here is unstable when p is close to 1,
102  * making it impossible to compute the inverse of Q(a,x) for small
103  * q. Fortunately the second form works perfectly well in this case.
104  */
105  double u;
106  if((b * q > 1e-8) && (q > 1e-5)) {
107  u = pow(p * g * a, 1 / a);
108  }
109  else {
110  u = exp((-q / a) - SCIPY_EULER);
111  }
112  result = u / (1 - (u / (a + 1)));
113  }
114  else if ((a < 0.3) && (b >= 0.35)) {
115  /* DiDonato & Morris Eq 22: */
116  double t = exp(-SCIPY_EULER - b);
117  double u = t * exp(t);
118  result = t * exp(u);
119  }
120  else if ((b > 0.15) || (a >= 0.3)) {
121  /* DiDonato & Morris Eq 23: */
122  double y = -log(b);
123  double u = y - (1 - a) * log(y);
124  result = y - (1 - a) * log(u) - log(1 + (1 - a) / (1 + u));
125  }
126  else if (b > 0.1) {
127  /* DiDonato & Morris Eq 24: */
128  double y = -log(b);
129  double u = y - (1 - a) * log(y);
130  result = y - (1 - a) * log(u)
131  - log((u * u + 2 * (3 - a) * u + (2 - a) * (3 - a))
132  / (u * u + (5 - a) * u + 2));
133  }
134  else {
135  /* DiDonato & Morris Eq 25: */
136  double y = -log(b);
137  double c1 = (a - 1) * log(y);
138  double c1_2 = c1 * c1;
139  double c1_3 = c1_2 * c1;
140  double c1_4 = c1_2 * c1_2;
141  double a_2 = a * a;
142  double a_3 = a_2 * a;
143 
144  double c2 = (a - 1) * (1 + c1);
145  double c3 = (a - 1) * (-(c1_2 / 2)
146  + (a - 2) * c1
147  + (3 * a - 5) / 2);
148  double c4 = (a - 1) * ((c1_3 / 3) - (3 * a - 5) * c1_2 / 2
149  + (a_2 - 6 * a + 7) * c1
150  + (11 * a_2 - 46 * a + 47) / 6);
151  double c5 = (a - 1) * (-(c1_4 / 4)
152  + (11 * a - 17) * c1_3 / 6
153  + (-3 * a_2 + 13 * a -13) * c1_2
154  + (2 * a_3 - 25 * a_2 + 72 * a - 61) * c1 / 2
155  + (25 * a_3 - 195 * a_2 + 477 * a - 379) / 12);
156 
157  double y_2 = y * y;
158  double y_3 = y_2 * y;
159  double y_4 = y_2 * y_2;
160  result = y + c1 + (c2 / y) + (c3 / y_2) + (c4 / y_3) + (c5 / y_4);
161  }
162  }
163  else {
164  /* DiDonato and Morris Eq 31: */
165  double s = find_inverse_s(p, q);
166 
167  double s_2 = s * s;
168  double s_3 = s_2 * s;
169  double s_4 = s_2 * s_2;
170  double s_5 = s_4 * s;
171  double ra = sqrt(a);
172 
173  double w = a + s * ra + (s_2 - 1) / 3;
174  w += (s_3 - 7 * s) / (36 * ra);
175  w -= (3 * s_4 + 7 * s_2 - 16) / (810 * a);
176  w += (9 * s_5 + 256 * s_3 - 433 * s) / (38880 * a * ra);
177 
178  if ((a >= 500) && (fabs(1 - w / a) < 1e-6)) {
179  result = w;
180  }
181  else if (p > 0.5) {
182  if (w < 3 * a) {
183  result = w;
184  }
185  else {
186  double D = fmax(2, a * (a - 1));
187  double lg = lgam(a);
188  double lb = log(q) + lg;
189  if (lb < -D * 2.3) {
190  /* DiDonato and Morris Eq 25: */
191  double y = -lb;
192  double c1 = (a - 1) * log(y);
193  double c1_2 = c1 * c1;
194  double c1_3 = c1_2 * c1;
195  double c1_4 = c1_2 * c1_2;
196  double a_2 = a * a;
197  double a_3 = a_2 * a;
198 
199  double c2 = (a - 1) * (1 + c1);
200  double c3 = (a - 1) * (-(c1_2 / 2)
201  + (a - 2) * c1
202  + (3 * a - 5) / 2);
203  double c4 = (a - 1) * ((c1_3 / 3)
204  - (3 * a - 5) * c1_2 / 2
205  + (a_2 - 6 * a + 7) * c1
206  + (11 * a_2 - 46 * a + 47) / 6);
207  double c5 = (a - 1) * (-(c1_4 / 4)
208  + (11 * a - 17) * c1_3 / 6
209  + (-3 * a_2 + 13 * a -13) * c1_2
210  + (2 * a_3 - 25 * a_2 + 72 * a - 61) * c1 / 2
211  + (25 * a_3 - 195 * a_2 + 477 * a - 379) / 12);
212 
213  double y_2 = y * y;
214  double y_3 = y_2 * y;
215  double y_4 = y_2 * y_2;
216  result = y + c1 + (c2 / y) + (c3 / y_2) + (c4 / y_3) + (c5 / y_4);
217  }
218  else {
219  /* DiDonato and Morris Eq 33: */
220  double u = -lb + (a - 1) * log(w) - log(1 + (1 - a) / (1 + w));
221  result = -lb + (a - 1) * log(u) - log(1 + (1 - a) / (1 + u));
222  }
223  }
224  }
225  else {
226  double z = w;
227  double ap1 = a + 1;
228  double ap2 = a + 2;
229  if (w < 0.15 * ap1) {
230  /* DiDonato and Morris Eq 35: */
231  double v = log(p) + lgam(ap1);
232  z = exp((v + w) / a);
233  s = log1p(z / ap1 * (1 + z / ap2));
234  z = exp((v + z - s) / a);
235  s = log1p(z / ap1 * (1 + z / ap2));
236  z = exp((v + z - s) / a);
237  s = log1p(z / ap1 * (1 + z / ap2 * (1 + z / (a + 3))));
238  z = exp((v + z - s) / a);
239  }
240 
241  if ((z <= 0.01 * ap1) || (z > 0.7 * ap1)) {
242  result = z;
243  }
244  else {
245  /* DiDonato and Morris Eq 36: */
246  double ls = log(didonato_SN(a, z, 100, 1e-4));
247  double v = log(p) + lgam(ap1);
248  z = exp((v + z - ls) / a);
249  result = z * (1 - (a * log(z) - z - v + ls) / (a - z));
250  }
251  }
252  }
253  return result;
254 }
255 
256 
257 double igami(double a, double p)
258 {
259  int i;
260  double x, fac, f_fp, fpp_fp;
261 
262  if (isnan(a) || isnan(p)) {
263  return NAN;
264  }
265  else if ((a < 0) || (p < 0) || (p > 1)) {
266  sf_error("gammaincinv", SF_ERROR_DOMAIN, NULL);
267  }
268  else if (p == 0.0) {
269  return 0.0;
270  }
271  else if (p == 1.0) {
272  return INFINITY;
273  }
274  else if (p > 0.9) {
275  return igamci(a, 1 - p);
276  }
277 
278  x = find_inverse_gamma(a, p, 1 - p);
279  /* Halley's method */
280  for (i = 0; i < 3; i++) {
281  fac = igam_fac(a, x);
282  if (fac == 0.0) {
283  return x;
284  }
285  f_fp = (igam(a, x) - p) * x / fac;
286  /* The ratio of the first and second derivatives simplifies */
287  fpp_fp = -1.0 + (a - 1) / x;
288  if (isinf(fpp_fp)) {
289  /* Resort to Newton's method in the case of overflow */
290  x = x - f_fp;
291  }
292  else {
293  x = x - f_fp / (1.0 - 0.5 * f_fp * fpp_fp);
294  }
295  }
296 
297  return x;
298 }
299 
300 
301 double igamci(double a, double q)
302 {
303  int i;
304  double x, fac, f_fp, fpp_fp;
305 
306  if (isnan(a) || isnan(q)) {
307  return NAN;
308  }
309  else if ((a < 0.0) || (q < 0.0) || (q > 1.0)) {
310  sf_error("gammainccinv", SF_ERROR_DOMAIN, NULL);
311  }
312  else if (q == 0.0) {
313  return INFINITY;
314  }
315  else if (q == 1.0) {
316  return 0.0;
317  }
318  else if (q > 0.9) {
319  return igami(a, 1 - q);
320  }
321 
322  x = find_inverse_gamma(a, 1 - q, q);
323  for (i = 0; i < 3; i++) {
324  fac = igam_fac(a, x);
325  if (fac == 0.0) {
326  return x;
327  }
328  f_fp = (igamc(a, x) - q) * x / (-fac);
329  fpp_fp = -1.0 + (a - 1) / x;
330  if (isinf(fpp_fp)) {
331  x = x - f_fp;
332  }
333  else {
334  x = x - f_fp / (1.0 - 0.5 * f_fp * fpp_fp);
335  }
336  }
337 
338  return x;
339 }
w
RowVector3d w
Definition: Matrix_resize_int.cpp:3
igamc
double igamc(double a, double x)
Definition: igam.c:169
D
MatrixXcd D
Definition: EigenSolver_EigenSolver_MatrixType.cpp:14
s
RealScalar s
Definition: level1_cplx_impl.h:126
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
didonato_SN
static double didonato_SN(double, double, unsigned, double)
Definition: igami.c:43
b
Scalar * b
Definition: benchVecAdd.cpp:17
x
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
Definition: gnuplot_common_settings.hh:12
Gamma
double Gamma(double x)
Definition: gamma.c:160
igam
double igam(double a, double x)
Definition: igam.c:128
log
const EIGEN_DEVICE_FUNC LogReturnType log() const
Definition: ArrayCwiseUnaryOps.h:128
isnan
#define isnan(X)
Definition: main.h:93
exp
const EIGEN_DEVICE_FUNC ExpReturnType exp() const
Definition: ArrayCwiseUnaryOps.h:97
lgam
double lgam(double x)
Definition: gamma.c:275
result
Values result
Definition: OdometryOptimize.cpp:8
SCIPY_EULER
#define SCIPY_EULER
Definition: mconf.h:128
boost::multiprecision::fabs
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:119
c1
static double c1
Definition: airy.c:54
find_inverse_s
static double find_inverse_s(double, double)
Definition: igami.c:14
Eigen::numext::q
EIGEN_DEVICE_FUNC const Scalar & q
Definition: SpecialFunctionsImpl.h:1984
polevl
static double polevl(double x, const double coef[], int N)
Definition: polevl.h:65
igam_fac
double igam_fac(double a, double x)
Definition: igam.c:231
log1p
double log1p(double x)
Definition: unity.c:49
pybind_wrapper_test_script.z
z
Definition: pybind_wrapper_test_script.py:61
g
void g(const string &key, int i)
Definition: testBTree.cpp:41
y
Scalar * y
Definition: level1_cplx_impl.h:124
Eigen::ArrayBase::pow
const Eigen::CwiseBinaryOp< Eigen::internal::scalar_pow_op< typename Derived::Scalar, typename ExponentDerived::Scalar >, const Derived, const ExponentDerived > pow(const Eigen::ArrayBase< Derived > &x, const Eigen::ArrayBase< ExponentDerived > &exponents)
Definition: GlobalFunctions.h:143
a
ArrayXXi a
Definition: Array_initializer_list_23_cxx11.cpp:1
mconf.h
Eigen::bfloat16_impl::fmax
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 fmax(const bfloat16 &a, const bfloat16 &b)
Definition: BFloat16.h:587
find_inverse_gamma
static double find_inverse_gamma(double, double, double)
Definition: igami.c:72
p
float * p
Definition: Tutorial_Map_using.cpp:9
sf_error
void sf_error(const char *func_name, sf_error_t code, const char *fmt,...)
Definition: sf_error.c:41
c2
static double c2
Definition: airy.c:55
v
Array< int, Dynamic, 1 > v
Definition: Array_initializer_list_vector_cxx11.cpp:1
N
#define N
Definition: igam.h:9
NULL
#define NULL
Definition: ccolamd.c:609
align_3::t
Point2 t(10, 10)
SF_ERROR_DOMAIN
@ SF_ERROR_DOMAIN
Definition: sf_error.h:16
ceres::sqrt
Jet< T, N > sqrt(const Jet< T, N > &f)
Definition: jet.h:418
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
igamci
double igamci(double a, double q)
Definition: igami.c:301
igami
double igami(double a, double p)
Definition: igami.c:257
isinf
#define isinf(X)
Definition: main.h:94

gtsam
Author(s):
autogenerated on Thu Jun 13 2024 03:02:35