Stable Numerical Evaluation of Erlang C

Why Erlang C Needs a Numerical Method

For an M/M/c queue, Erlang C is the stationary probability that an arriving customer must wait because all \(c\) servers are busy. Its formula is compact, but a direct computation can be unreliable when \(c\) is large: terms such as \(a^n/n!\) may overflow or underflow even when their ratios, and the final probability, are ordinary-sized. The remedy is to compute relative weights recursively, scale them before summing, and treat the infinite tail as a geometric series.

Assume \(c\) is a positive integer, the arrival rate is \(\lambda\ge0\), and the service rate of each server is \(\mu>0\). Write \(a=\lambda/\mu\) for offered traffic. When \(\lambda>0\), stability requires \(a<c\); equivalently, the utilization \(\rho=a/c\) is less than \(1\). If \(\lambda=0\), no customer arrives and the delay probability is zero.

The queue-length process \((X_t)_{t\ge0}\) is a continuous-time Markov chain on \(S=\{0,1,2,\ldots\}\), on a probability space \((\Omega,\mathcal F,\mathbb P)\). In state \(n\), arrivals occur at rate \(\lambda\) and departures at rate \(\min(n,c)\mu\). Thus its generator has \(q_{n,n+1}=\lambda\), \(q_{n,n-1}=\min(n,c)\mu\) for \(n\ge1\), and \(q_{n,n}=-(\lambda+\min(n,c)\mu)\); the diagonal makes every generator row sum to zero. The stationary state probabilities are proportional to the weights \(w_n=a^n/n!\) for \(0\le n<c\) and \(w_n=a^n/(c!c^{n-c})\) for \(n\ge c\). These are the stationary weights for the stable birth-death queue; only their ratios are needed below.

Scaling and Recursion

Evaluating powers and factorials separately is unnecessary. For the finite set of weights \(w_0,\ldots,w_{c-1}\), choose an index \(m=\min(\lfloor a\rfloor,c-1)\) and scale by \(w_m\). Define \(x_n=w_n/w_m\). Since \(w_{n+1}/w_n=a/(n+1)\), the relative weights can be generated from \(x_m=1\) by stepping in either direction:

Set \(S_c=\sum_{n=0}^{c-1}x_n\), and extend the recursion once to obtain \(x_c=x_{c-1}a/c\). Because every weight has been divided by the same positive number \(w_m\), the scaled finite sum is \(S_c\), and the scaled tail weight is \(x_c/(1-\rho)\). Therefore the probability can be computed without constructing any unscaled \(w_n\).

For numerical work, it is often preferable to avoid forming the potentially large tail weight \(x_c/(1-\rho)\). Dividing numerator and denominator by that tail weight gives an equivalent expression:

When \(a\) and \(c\) are available, compute \(1-\rho\) as \((c-a)/c\), rather than subtracting a rounded value of \(\rho\) from \(1\). The recursion uses \(O(c)\) multiplications and additions. In floating-point arithmetic, a very small \(x_c\) or an extreme ratio can still exceed the representable range; a log-scale calculation is then safer. The scaled recursion greatly reduces avoidable overflow, but it cannot make every possible input representable.

Proof. For \(0\le n<c\), define \(w_n=a^n/n!\), where \(a>0\). Direct division gives \(w_{n+1}/w_n=a/(n+1)\) and \(w_{n-1}/w_n=n/a\) whenever the indicated indices exist. Starting at \(x_m=1=w_m/w_m\), each upward recursion therefore produces \(x_{n+1}=w_{n+1}/w_m\), and each downward recursion produces \(x_{n-1}=w_{n-1}/w_m\). Induction in each direction proves \(x_n=w_n/w_m\) for all \(0\le n<c\), and the extra upward step proves the same identity for \(x_c\). For \(n\ge c\), the stationary weights have successive ratio \(a/c=\rho\), so their total is \(w_c\sum_{r=0}^{\infty}\rho^r=w_c/(1-\rho)\); this series converges because \(0<\rho<1\). Dividing both finite and tail weights by \(w_m>0\) gives \(S_c\) and \(x_c/(1-\rho)\). Substitution into the defining probability proves the scaled formula. Dividing its numerator and denominator by \(x_c/(1-\rho)>0\) gives \(1/[1+S_c(1-\rho)/x_c]\), as claimed. \(\square\)

This argument also explains the role of the anchor index. The ratio \(w_{n+1}/w_n=a/(n+1)\) is at least \(1\) up to the mode and at most \(1\) after it. The chosen \(m\) is a mode of the finite sequence \(w_0,\ldots,w_{c-1}\), so \(0<x_n\le1\) throughout that sequence. Anchoring at an extreme index instead can create very large relative weights, particularly for low offered traffic.

Worked Calculations and Checks

Common Numerical Pitfalls

Common mistake: counting only the state \(X_t=c\) as a delay. In the three-server example, using \(w_3=4/3\) as the entire tail gives the incorrect calculation \((4/3)/(5+4/3)=4/19\). A delay occurs in every state \(n\ge3\), not just state \(3\). The correct tail is \(w_3(1+\rho+\rho^2+\cdots)=(4/3)/(1-2/3)=4\), giving \(4/(5+4)=4/9\).

Other checks help distinguish arithmetic errors from a misapplied formula:

• Confirm \(c\ge1\), \(\mu>0\), and \(a<c\) before using a stationary Erlang C probability. If \(\lambda=0\), handle that case separately and return zero.
• Generate each relative term from its neighbor; do not independently evaluate large powers and factorials if they can overflow.
• Use the geometric tail \(w_c/(1-\rho)\), not just its first term \(w_c\).
• Check that the finite-state probabilities and the aggregated tail probability sum to \(1\), and compare expected busy servers with \(a\) when performing a stationary-distribution check.
• Near saturation, retain \(1-\rho=(c-a)/c\) as directly as possible. If the final ratios leave the floating-point range, move the calculation to logarithms rather than trusting an underflowed term.

These checks are especially valuable at large server counts. The recursion changes the calculation from a sum of potentially unmanageable factorial terms into a sequence of local ratios. Scaling preserves the probability exactly in real arithmetic, while normalization, utilization, and limiting cases provide independent tests of the numerical result.