rgamma.c
Go to the documentation of this file.
1 /* rgamma.c
2  *
3  * Reciprocal Gamma function
4  *
5  *
6  *
7  * SYNOPSIS:
8  *
9  * double x, y, rgamma();
10  *
11  * y = rgamma( x );
12  *
13  *
14  *
15  * DESCRIPTION:
16  *
17  * Returns one divided by the Gamma function of the argument.
18  *
19  * The function is approximated by a Chebyshev expansion in
20  * the interval [0,1]. Range reduction is by recurrence
21  * for arguments between -34.034 and +34.84425627277176174.
22  * 0 is returned for positive arguments outside this
23  * range. For arguments less than -34.034 the cosecant
24  * reflection formula is applied; lograrithms are employed
25  * to avoid unnecessary overflow.
26  *
27  * The reciprocal Gamma function has no singularities,
28  * but overflow and underflow may occur for large arguments.
29  * These conditions return either INFINITY or 0 with
30  * appropriate sign.
31  *
32  * ACCURACY:
33  *
34  * Relative error:
35  * arithmetic domain # trials peak rms
36  * IEEE -30,+30 30000 1.1e-15 2.0e-16
37  * For arguments less than -34.034 the peak error is on the
38  * order of 5e-15 (DEC), excepting overflow or underflow.
39  */
40 
41 /*
42  * Cephes Math Library Release 2.0: April, 1987
43  * Copyright 1985, 1987 by Stephen L. Moshier
44  * Direct inquiries to 30 Frost Street, Cambridge, MA 02140
45  */
46 
47 #include "mconf.h"
48 
49 /* Chebyshev coefficients for reciprocal Gamma function
50  * in interval 0 to 1. Function is 1/(x Gamma(x)) - 1
51  */
52 
53 static double R[] = {
54  3.13173458231230000000E-17,
55  -6.70718606477908000000E-16,
56  2.20039078172259550000E-15,
57  2.47691630348254132600E-13,
58  -6.60074100411295197440E-12,
59  5.13850186324226978840E-11,
60  1.08965386454418662084E-9,
61  -3.33964630686836942556E-8,
62  2.68975996440595483619E-7,
63  2.96001177518801696639E-6,
64  -8.04814124978471142852E-5,
65  4.16609138709688864714E-4,
66  5.06579864028608725080E-3,
67  -6.41925436109158228810E-2,
68  -4.98558728684003594785E-3,
69  1.27546015610523951063E-1
70 };
71 
72 static char name[] = "rgamma";
73 
74 extern double MAXLOG;
75 
76 
77 double rgamma(double x)
78 {
79  double w, y, z;
80  int sign;
81 
82  if (x > 34.84425627277176174) {
83  return exp(-lgam(x));
84  }
85  if (x < -34.034) {
86  w = -x;
87  z = sinpi(w);
88  if (z == 0.0) {
89  return 0.0;
90  }
91  if (z < 0.0) {
92  sign = 1;
93  z = -z;
94  }
95  else {
96  sign = -1;
97  }
98 
99  y = log(w * z) - log(M_PI) + lgam(w);
100  if (y < -MAXLOG) {
102  return (sign * 0.0);
103  }
104  if (y > MAXLOG) {
106  return (sign * INFINITY);
107  }
108  return (sign * exp(y));
109  }
110  z = 1.0;
111  w = x;
112 
113  while (w > 1.0) { /* Downward recurrence */
114  w -= 1.0;
115  z *= w;
116  }
117  while (w < 0.0) { /* Upward recurrence */
118  z /= w;
119  w += 1.0;
120  }
121  if (w == 0.0) /* Nonpositive integer */
122  return (0.0);
123  if (w == 1.0) /* Other integer */
124  return (1.0 / z);
125 
126  y = w * (1.0 + chbevl(4.0 * w - 2.0, R, 16)) / z;
127  return (y);
128 }
w
RowVector3d w
Definition: Matrix_resize_int.cpp:3
name
Annotation for function names.
Definition: attr.h:51
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
sign
const EIGEN_DEVICE_FUNC SignReturnType sign() const
Definition: ArrayCwiseUnaryOps.h:219
SF_ERROR_OVERFLOW
@ SF_ERROR_OVERFLOW
Definition: sf_error.h:12
log
const EIGEN_DEVICE_FUNC LogReturnType log() const
Definition: ArrayCwiseUnaryOps.h:128
exp
const EIGEN_DEVICE_FUNC ExpReturnType exp() const
Definition: ArrayCwiseUnaryOps.h:97
lgam
double lgam(double x)
Definition: gamma.c:275
chbevl
double chbevl(double x, double array[], int n)
Definition: chbevl.c:63
R
static double R[]
Definition: rgamma.c:53
SF_ERROR_UNDERFLOW
@ SF_ERROR_UNDERFLOW
Definition: sf_error.h:11
rgamma
double rgamma(double x)
Definition: rgamma.c:77
pybind_wrapper_test_script.z
z
Definition: pybind_wrapper_test_script.py:61
MAXLOG
double MAXLOG
Definition: const.c:57
y
Scalar * y
Definition: level1_cplx_impl.h:124
mconf.h
sf_error
void sf_error(const char *func_name, sf_error_t code, const char *fmt,...)
Definition: sf_error.c:41
NULL
#define NULL
Definition: ccolamd.c:609
M_PI
#define M_PI
Definition: mconf.h:117
sinpi
double sinpi(double x)
Definition: sinpi.c:11


gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:35:00