- if (x == 0.) return 0.;
- // The following line is not thoroughly tested, so it is commented out.
- if (p > 1e3) return .5 * kf_erfc(-M_SQRT1_2 * sqrt(p) * 3. * (pow(x / p, 1./3.) + 1. / (p * 9.) - 1.));
- if (x > 1e8) return 1.;
- if (x <= 1. || x < p) { // series expansion
- c = 1.;
- arg = p * log(x) - x - kf_lgamma(p + 1.);
- ret_val = 1.;
- a = p;
- while (c > 1e-14) {
- a += 1.;
- c = c * x / a;
- ret_val += c;
- }
- arg += log(ret_val);
- ret_val = 0.;
- if (arg >= -88.) ret_val = exp(arg);
- } else { // continued expansion
- arg = p * log(x) - x - kf_lgamma(p);
- a = 1. - p;
- b = a + x + 1.;
- c = 0.;
- pn1 = 1.;
- pn2 = x;
- pn3 = x + 1.;
- pn4 = x * b;
- ret_val = pn3 / pn4;
- while (1) {
- a += 1.; b += 2.; c += 1.;
- an = a * c;
- pn5 = b * pn3 - an * pn1;
- pn6 = b * pn4 - an * pn2;
- if (fabs(pn6) > 0.) {
- rn = pn5 / pn6;
- if (fabs(ret_val - rn) <= fmin(1e-14, rn * 1e-14)) break;
- ret_val = rn;
- }
- pn1 = pn3; pn2 = pn4; pn3 = pn5; pn4 = pn6;
- if (fabs(pn5) >= 1e37)
- pn1 /= 1e37, pn2 /= 1e37, pn3 /= 1e37, pn4 /= 1e37;
- }
- arg += log(ret_val);
- ret_val = 1.;
- if (arg >= -88.) ret_val = 1. - exp(arg);
- }
- return ret_val;
+double kf_gammap(double s, double z)
+{
+ return z <= 1. || z < s? _kf_gammap(s, z) : 1. - _kf_gammaq(s, z);