- c = lambda / 2.;
- if (lambda < 54.)
- {
- /* Algorithm AS 226. */
- double x0, a0, beta, temp, gx, q, ax, sumq, sum;
- double err_max = 2 * DBL_EPSILON;
- double err_bound;
- int iter_max = 100;
- int iter;
-
- x0 = floor (c - 5.0 * sqrt (c));
- if (x0 < 0.)
- x0 = 0.;
- a0 = a + x0;
- beta = (gsl_sf_lngamma (a0)
- + gsl_sf_lngamma (b)
- - gsl_sf_lngamma (a0 + b));
- temp = gsl_sf_beta_inc (a0, b, x);
- gx = exp (a0 * log (x) + b * log (1. - x) - beta - log (a0));
- if (a0 >= a)
- q = exp (-c + x0 * log (c)) - gsl_sf_lngamma (x0 + 1.);
- else
- q = exp (-c);
- ax = q * temp;
- sumq = 1. - q;
- sum = ax;
-
- iter = 0;
- do
- {
- iter++;
- temp -= gx;
- gx = x * (a + b + iter - 1.) * gx / (a + iter);
- q *= c / iter;
- sumq -= q;
- ax = temp * q;
- sum += ax;
-
- err_bound = (temp - gx) * sumq;
- }
- while (iter < iter_max && err_bound > err_max);
-
- return sum;
- }
- else
- {
- /* Algorithm AS 310. */
- double m, m_sqrt;
- int iter, iter_lower, iter_upper, iter1, iter2, j;
- double t, q, r, psum, beta, s1, gx, fx, temp, ftemp, t0, s0, sum, s;
- double err_bound;
- double err_max = 2 * DBL_EPSILON;
-
- iter = 0;
-
- m = floor (c + .5);
- m_sqrt = sqrt (m);
- iter_lower = m - 5. * m_sqrt;
- iter_upper = m + 5. * m_sqrt;
-
- t = -c + m * log (c) - gsl_sf_lngamma (m + 1.);
- q = exp (t);
- r = q;
- psum = q;
- beta = (gsl_sf_lngamma (a + m)
- + gsl_sf_lngamma (b)
- - gsl_sf_lngamma (a + m + b));
- s1 = (a + m) * log (x) + b * log (1. - x) - log (a + m) - beta;
- fx = gx = exp (s1);
- ftemp = temp = gsl_sf_beta_inc (a + m, b, x);
- iter++;
- sum = q * temp;
- iter1 = m;
-
- while (iter1 >= iter_lower && q >= err_max)
- {
- q = q * iter1 / c;
- iter++;
- gx = (a + iter1) / (x * (a + b + iter1 - 1.)) * gx;
- iter1--;
- temp += gx;
- psum += q;
- sum += q * temp;
- }
-
- t0 = (gsl_sf_lngamma (a + b)
- - gsl_sf_lngamma (a + 1.)
- - gsl_sf_lngamma (b));
- s0 = a * log (x) + b * log (1. - x);
-
- s = 0.;
- for (j = 0; j < iter1; j++)
- {
- double t1;
- s += exp (t0 + s0 + j * log (x));
- t1 = log (a + b + j) - log (a + 1. + j) + t0;
- t0 = t1;
- }
-
- err_bound = (1. - gsl_sf_gamma_inc_P (iter1, c)) * (temp + s);
- q = r;
- temp = ftemp;
- gx = fx;
- iter2 = m;
- for (;;)
- {
- double ebd = err_bound + (1. - psum) * temp;
- if (ebd < err_max || iter >= iter_upper)
- break;
-
- iter2++;
- iter++;
- q = q * c / iter2;
- psum += q;
- temp -= gx;
- gx = x * (a + b + iter2 - 1.) / (a + iter2) * gx;
- sum += q * temp;
- }
-
- return sum;
- }