Erlang C for Large Numbers The formula for Erlang C gives the probability of not having your call answered immediately at a call center, so a number between 0 and 1:
$${E_C} = {1 \over {1 + \left( {1 - p} \right){{m!} \over {{u^m}}}\sum\limits_{k = 0}^{m - 1} {{{{u^k}} \over {k!}}} }}$$
p (agent occupancy) is a number between 0 and 1. With m (agents on hand) and u (call intensity, m>u) large (e.g., u = 143, m = 144), 64-bit floating-point arithmetic overflows -- so I'm trying to figure out how to calculate the denominator in log space, i.e.,
$${E_C} = {1 \over {1 + \left( {1 - p} \right)\exp \left( {\ln \left( {{{m!} \over {{u^m}}}\sum\limits_{k = 0}^{m - 1} {{{{u^k}} \over {k!}}} } \right)} \right)}}$$
... and got this far:
$$\displaylines{
  \ln \left( {{{m!} \over {{u^m}}}\sum\limits_{k = 0}^{m - 1} {{{{u^k}} \over {k!}}} } \right) = \ln m! - \ln {u^m} + \ln \left( {\sum\limits_{k = 0}^{m - 1} {{{{u^k}} \over {k!}}} } \right) \cr 
   = \ln m! - m\ln u + \ln \left( {\sum\limits_{k = 0}^{m - 1} {{u^k}} } \right) - \ln \left( {\sum\limits_{k = 0}^{m - 1} {k!} } \right) \cr} $$
The first two terms are no problem (Stirling's Approximation for the first), but am stuck on the last two. I looked at the log identity for a sum
$$\ln \sum\limits_{i = 0}^N {{a_i}}  = \ln {a_0} + \ln \left( {1 + \sum\limits_{i = 1}^N {{{{a_i}} \over {{a_0}}}} } \right)$$
... but it doesn't get me anywhere.
This is an exercise for fun, and I am not a mathematician, but would welcome a suggestion.
In searching, I saw this same question asked at https://forums.adobe.com/thread/841586, but there was a deathly silence. I'm optimistic the person was just asking on the wrong forum.
Thanks for reading.
EDIT: In the alternative, might there be a way to restructure the summation so that the m!/u^m term (a little tiny number) can be calculated incrementally to moderate the gargantuan u^k/k! terms? In that context, I reckon this is more of a numerical analysis or programming problem than math. Would Stack Overflow be a more appropriate forum?
 A: The question author shg has come up with a practical solution and posted a VBA function in the Adobe thread linked above.
This answer describes an alternate solution in terms of built-in functions not requiring iteration. Unfortunately it doesn't seem to work accurately in Excel, but it should work in R or Matlab.
The sum can be expressed in terms of a cumulative distribution function:
\begin{align*}
\frac{m!}{u^m}\sum\limits_{k = 0}^{m - 1} {{{{u^k}} \over {k!}}}
&=\frac{m!}{u^m}e^uF_{\mathrm{Poisson}}(m-1;u)\\
&=\exp(\ln(\Gamma(m+1))-m\ln(u)+u)F_{\mathrm{Poisson}}(m-1;u)
\end{align*}
where $F_{\mathrm{Poisson}}(m-1;u)$ is the probability that a $\mathrm{Poisson}(u)$ variable is at most $m-1$, and $\Gamma(m+1)$ denotes the Gamma function. The $F_{\mathrm{Poisson}}(m-1;u)$ factor is numerically nice because it is not too small, in fact it looks like
$$\tfrac 1 2 < F_{\mathrm{Poisson}}(m-1;u)<1\text{ for $u<m$.}$$
(I don't know a rigorous proof for the lower bound, but by the normal approximation for Poisson variables, for large $u$ the median will be approximately the mean, giving $F_{\mathrm{Poisson}}(m-1;u)\geq F_{\mathrm{Poisson}}(u;u)\to \tfrac 1 2$.)
Excel's POISSON function can be used to evaluate $F_{\mathrm{Poisson}}$, so the whole expression could be expressed as:

EXP(GAMMALN(m+1)-m*LN(u)+u)*POISSON(m-1,u,TRUE)

Unfortunately Excel's POISSON can't be trusted here; my copy of Excel gives POISSON(1e6,1e6,TRUE)=0.863245255, where R gives the more plausible ppois(1e6,lambda=1e6)=0.500266.
A: Since the summation is over $k$, you cannot have $k$ in the result.
$$\sum\limits_{k = 0}^{m - 1} {{{{u^k}} \over {k!}}}=\frac{ \Gamma (m,u)}{\Gamma (m)}\,e^u$$ where appears the incomplete gamma function.
$$\frac {m!}{u^m}\sum\limits_{k = 0}^{m - 1} {{{{u^k}} \over {k!}}}=m \,e^u\, u^{-m} \,\Gamma (m,u)$$
Have a look here for the calculation of the incomplete gamma function using Excel.
