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