🚀 go-pugleaf

        I'm trying to evaluate the following series, and wondering if
it can be reduced further (toward some kind of special function).

sum[k=0 to +inf] of [b^k * e^(-b) * P(k+1,a) / K!]

where k is integer, a and b are positive real, e=2.718281828..., and
P() is the incomplete gamma function.

        The only problem is that both a and b can be (and most of the
time in my applications) as large as 5000. Because of this, I have to
put e^(-b) in each term even it has nothing to do with k, and evaluate
log of each term first and then take its exponential,

exp( k*LN(b) - b + LN(P(k+1,a)) - LN(k!) )

otherwise e^(-5000) is effectively zero, and 5000^k will overflow, on
most computers.

        Please let me know if there is better algorithm to evaluate
this series, or if it has an analytical solution. Thanks for your time
and effort.

Ping