poch.c
Go to the documentation of this file.
1 /*
2  * Pochhammer symbol (a)_m = gamma(a + m) / gamma(a)
3  */
4 #include "mconf.h"
5 
6 static double is_nonpos_int(double x)
7 {
8  return x <= 0 && x == ceil(x) && fabs(x) < 1e13;
9 }
10 
11 double poch(double a, double m)
12 {
13  double r;
14 
15  r = 1.0;
16 
17  /*
18  * 1. Reduce magnitude of `m` to |m| < 1 by using recurrence relations.
19  *
20  * This may end up in over/underflow, but then the function itself either
21  * diverges or goes to zero. In case the remainder goes to the opposite
22  * direction, we end up returning 0*INF = NAN, which is OK.
23  */
24 
25  /* Recurse down */
26  while (m >= 1.0) {
27  if (a + m == 1) {
28  break;
29  }
30  m -= 1.0;
31  r *= (a + m);
32  if (!isfinite(r) || r == 0) {
33  break;
34  }
35  }
36 
37  /* Recurse up */
38  while (m <= -1.0) {
39  if (a + m == 0) {
40  break;
41  }
42  r /= (a + m);
43  m += 1.0;
44  if (!isfinite(r) || r == 0) {
45  break;
46  }
47  }
48 
49  /*
50  * 2. Evaluate function with reduced `m`
51  *
52  * Now either `m` is not big, or the `r` product has over/underflown.
53  * If so, the function itself does similarly.
54  */
55 
56  if (m == 0) {
57  /* Easy case */
58  return r;
59  }
60  else if (a > 1e4 && fabs(m) <= 1) {
61  /* Avoid loss of precision */
62  return r * pow(a, m) * (
63  1
64  + m*(m-1)/(2*a)
65  + m*(m-1)*(m-2)*(3*m-1)/(24*a*a)
66  + m*m*(m-1)*(m-1)*(m-2)*(m-3)/(48*a*a*a)
67  );
68  }
69 
70  /* Check for infinity */
71  if (is_nonpos_int(a + m) && !is_nonpos_int(a) && a + m != m) {
72  return INFINITY;
73  }
74 
75  /* Check for zero */
76  if (!is_nonpos_int(a + m) && is_nonpos_int(a)) {
77  return 0;
78  }
79 
80  return r * exp(lgam(a + m) - lgam(a)) * gammasgn(a + m) * gammasgn(a);
81 }
gammasgn
double gammasgn(double x)
Definition: gammasgn.c:3
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
exp
const EIGEN_DEVICE_FUNC ExpReturnType exp() const
Definition: ArrayCwiseUnaryOps.h:97
lgam
double lgam(double x)
Definition: gamma.c:275
boost::multiprecision::fabs
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:119
is_nonpos_int
static double is_nonpos_int(double x)
Definition: poch.c:6
isfinite
#define isfinite(X)
Definition: main.h:95
poch
double poch(double a, double m)
Definition: poch.c:11
m
Matrix3f m
Definition: AngleAxis_mimic_euler.cpp:1
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
ceil
const EIGEN_DEVICE_FUNC CeilReturnType ceil() const
Definition: ArrayCwiseUnaryOps.h:495


gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:34:23