statistics - Fast incomplete gamma function -

what fast way calculate incomplete gamma function, or @ least "good" approximation of it, in c++?

background

what need calculate

given number of bernoulli trails n, probability p of success, i'm trying calculate probability of obtaining @ k successes, function of k. cumulative binomial distribution f(k,n,p) gives probability.

the need speed

i need calculate few hundred thousand of these cumulative probabilities per second. calculating cumulative binomial distribution straightforward summation computation-intensive large n. using incomplete beta function lot better, still quite computation intensive.

exploitable constraints

i'm hoping following constraints application domain can speeding calculation:

• p < 0.01 (the distribution skew)
• n > 50

poisson approximation

after experimentation in excel, i've learned poisson approximation excellent under above conditions. i.e. b(n,p) @ k identical pois(np) @ k under conditions of interest. means need function of 2 variables, no longer 3.

i understand cumulative poisson distribution can calculated in terms of incomplete gamma function, which, judging source code in cephes library, seems quite lot simpler calculate original incomplete beta function 1 have had calculate without poisson approximation. still isn't simple , iterative numerical calculation. i'm looking fast way calculate incomplete gamma function. i'm wondering whether there isn't closed-form expression can approximate reasonably well.

required precision

20% relative error quite acceptable on integral/probability (considered every k, in both directions).

i've considered using interpolated table poisson cdf directly, evenly-spaced domain-points less-than-ideal , domain have restricted arbitrary rectangle. analytic function quite number of tweaked parameters i'm hoping find ideally.

instead of using gamma function, concocted approximation transforms poisson variables standard normal variable:

float poisson_z(float x, float mu){     static const float twothirds = 2.0f/3.0f;     float w = sqrt((x+0.5f)/mu) - 1.0f;     float coeff = w>=0.0f ? 0.085f : 0.15f;     return (x-mu+twothirds)/sqrtf(mu*(1.0f+w*(0.68f+w*coeff))); } 

there no shortage of approximations standard normal distribution.