Pollaczek contour integrals for the fixed-cycle traffic-light queue

The fixed-cycle traffic-light (FCTL) queue is the standard model for intersections with static signaling, where vehicles arrive, form a queue and depart during cycles controlled by a traffic light. Classical analysis of the FCTL queue based on transform methods requires a computationally challenging step of finding the complex-valued roots of some characteristic equation. Building on the recent work of Oblakova et al. (Exact expected delay and distribution for the fixed-cycle traffic-light model and similar systems in explicit form, 2016), we obtain a contour-integral expression, reminiscent of Pollaczek integrals for bulk-service queues, for the probability generating function of the steady-state FCTL queue. We also show that similar contour integrals arise for generalizations of the FCTL queue introduced in Oblakova et al. (2016) that relax some of the classical assumptions. Our results allow to compute the queue-length distribution and all its moments using algorithms that rely on contour integrals and avoid root-finding procedures.


Introduction
The fixed-cycle traffic-light (FCTL) queue is an intensively studied stochastic model in traffic engineering [3,5,6,11,12,17,19]. Vehicles arrive at an intersection controlled by a traffic light and form a queue. The FCTL queue is traditionally modeled in discrete time, and time is divided into slots of unit length. The green and red periods, of length g and r slots, respectively, and thus the cycle length c = g +r , are assumed to be fixed multiples of one slot. Each slot corresponds to the time needed for a delayed vehicle to depart from the queue. Vehicles that arrive during a red period are delayed, as well as those that arrive during a green period and meet a non-empty queue. Those vehicles that arrive to the queue and are delayed join the queue at the end of the slot in which they arrive. Vehicles that arrive during a green period and meet no other vehicles in the queue are treated according to the following assumption: Definition 1 (FCTL assumption) For those cycles in which the queue clears before the green period terminates, all vehicles that arrive during the residual green period pass through the system and experience no delay whatsoever.
So, in the case of a non-empty queue during the green period, the queue length in the next slot is reduced by one compared to the queue length in the previous slot, but increased by the number of arrivals in that slot, whereas the queue remains empty if the queue was already empty at the start of the slot.
The FCTL assumption is quite realistic for straight-going traffic flows, since vehicles that find no queue during a residual green period will proceed without stopping and can keep driving at free-flow speed. For turning flows, however, the FCTL assumption seems too restrictive. Oblakova et al. [13] considered several generalizations of the FCTL queue that alleviate the classical FCTL assumption to account for rightturns, disruptions of the traffic and uncertainty in departure times. We will cover these extensions in Sect. 2.3.
The FCTL assumption (or its relaxed counterparts in [13]) causes some mathematical challenges. Let X k,n denote the queue length at time k in cycle n (time expressed in slots). Then, in cycle n, X 0,n is the queue length at the beginning of the green period, and X g,n the overflow defined as the queue length at the end of the green period (and thus the beginning of the red period). Let A n denote the total number of vehicles that arrive at the intersection in between the two measurements of the overflow X g,n and X g,n+1 . Thus, A n are the arrivals from the end of the time slot at which X g,n is measured onwards in a consecutive red and green period. Further, A n = A d n + A p n , where A d n denotes the number of delayed vehicles and A p n the number of vehicles that pass without delay on behalf of the FCTL assumption. The overflow queue can then be defined as X g,n+1 = max{X g,n + A d n − g, 0}. (1.1) The fact that A d n depends on both X g,n and the exact specification of when the arrivals occur makes (1.1) hard to analyze. To capture that level of detail, let Y k,n denote the number of vehicles that arrive at the intersection during slot k in cycle n. The random variables Y k,n are assumed independent and identically distributed, for all k and n. Notice that all of the above assumptions together mean that the queue lengths at the end of time slots can be modeled as a discrete-time Markov chain. Using analytic techniques suitable for dealing with such Markov chains, Darroch [5] obtained the probability generating function (PGF) of the steady-state overflow queue (the number of vehicles waiting in front of the traffic light at the end of a green period), and the PGF of the steady-state delay was obtained in van Leeuwaarden [17]. Hence, all information about the distribution of the steady-state overflow queue and steady-state delay in the FCTL queue can be obtained from the results in [5,17], including all moments of the steady-state queue length and delay, and the distribution of the output process (the way vehicles leave the intersection).
The FCTL queue belongs to a large class of cyclic queueing models related to vehicle dispatching with uncertain arrivals and bulk service [14,15,18]. A range of transportation and manufacturing systems can be modeled in this way, including batch production systems, bulk movements of goods in a factory, truck shipments and bus transportation. Within this class of models, many different rules can be considered that apply to customer arrivals and vehicle departures within a cycle. Think of vehiclecancelation policies that hold a vehicle until the queue length reaches a specified threshold. The FCTL assumption can also be viewed as a special rule that influences the dynamics within a cycle.
The classical FCTL treatment in [5,17] comes with computational challenges. This is because the PGF for the stationary queue-length distribution contains g−1 boundary terms that need to be found separately. The traditional way of determining these remaining unknowns consists of two steps: finding the g − 1 complex-valued roots and using these roots as input for a system of linear equations whose solution gives the boundary terms. Both steps can present difficulties, but were somehow considered unavoidable in the mathematically rigorous treatment of the FCTL queue [5,11,12,17,19] and of related bulk-service queues [9,14,15]. This perspective changed with the work of Oblakova et al. [13], which presented a contour-integral expression of the mean stationary queue length. The present paper extends that methodology to the transform domain and results in a contour-integral expression for the PGF of the FCTL queue and some of the generalizations considered in [13].
Paper outline We present in Sect. 2 the main result of the paper, an alternative representation for the transform solution in terms of a contour integral. Theorem 1 deals with the classical FCTL queue and Theorem 2 with the generalized FCTL queues introduced in [13]. We also explain in Sect. 3 how these contour integrals lead to algorithms that can compete with existing algorithms based on root finding. The proof of our main result is presented in Sect. 4 and uses several basic notions from complex analysis. We present some conclusions in Sect. 5.

Main results
We briefly review the standard solution method for obtaining an exact transform solution for the steady-state FCTL queue in Sect. 2.1. We then present in Sect. 2.2 the contour-integral representation for the FCTL queue and in Sect. 2.3 for the generalized FCTL queues.

Standard solution
Let Y be the number of arrivals during one slot and define Y (z) = E[z Y ]. Assume P(Y = 0) > 0, Y (1) < 1 and Y (z) to be analytic in a region |z| < R with R > 1 and R maximal. The key quantity in the mathematical analysis of the FCTL queue is the steady-state overflow queue, defined as X g = lim n→∞ X g,n . Clearly, to have stability, and for X g to be well defined, (2.1) Using the kernel method and transform techniques, the PGF of X g , denoted by X g (z) = E[z X g ], can be obtained using a by now classical line of reasoning.
This expression still contains g unknowns q 0 , . . . , q g−1 , representing the probability that the queue empties in slot k, thus P(X k = 0) = q k , which can be found by exploiting the analytic properties of PGFs. With Rouché's theorem, it can be shown that the denominator of (2.2) has g zeros on or within the unit circle |z| = 1. Because a PGF is well defined in |z| ≤ 1, the numerator of X g (z) should vanish at each of the zeros. This gives g equations. One of the zeros equals 1 and leads to a trivial equation. However, the normalization condition X g (1) = 1 provides an additional equation. That summarizes the highest level of general development for FCTL queue analysis: transform techniques yield an expression for X g (z) that in order to be evaluated demands finding g − 1 roots in the complex plane of the function z g = A(z) and solving a set of g linear equations.

Standard FCTL queue
We now turn to the alternative expression for X g (z), where we change the argument of X g from z to w to distinguish between the standard version of the PGF and the contour-integral version. Here is the main result in this paper:

Theorem 1 (Pollaczek integral for FCTL)
There is an 0 > 0 such that, for all ∈ (0, 0 ), with the principal value of the logarithm.
Here, 0 should satisfy the inequality 0 < min{t 0 , R 0 }, where t 0 = sup{t ∈ R + |Y (t)t − Y (t) ≤ 0} and R 0 is the unique root with smallest modulus of A(z) = z g in (1, ∞). Since E[Y ] < g/c < 1, this root always exists, also in the case Y (z) = y 0 + y 1 z. Sketch of the proof The proof of Theorem 1 finds a way to go from representation (2.2) to contour integrals. A significant start in this direction was made by [13], who rewrote (2.2) as

Remark 1 Formula (2.3) for X g (w) is essentially equivalent to
Denote the g roots of z g = A(z) on and within the unit circle by z 0 = 1, z 1 , . . . , z g−1 . Now, here is where the authors in [13] took an eye-opening step: instead of using the g roots in the traditional manner for finding the unknowns q k and completing transform (2.2), use these roots for factorizing the numerator of (2.2). Notice that this cannot be done immediately, because interpreted as a function of z, the numerator is by no means a polynomial of degree g or less. However, by treating the function Y (z)/z as a variable itself, the summation in the numerator is a polynomial of degree g − 1 and can be factorized as using that X g (z) is well defined in the disk |z| ≤ 1, that z 1 , . . . z g−1 are roots of the denominator and therefore also should be roots of the numerator, and that Y (z)/z is injective (see Sect. 4). After normalization using X g (1) = 1, the factorization in (2.6) leads to the representation Our proof then proceeds by interpreting (2.7) as the outcome of Cauchy's residue theorem, the classical tool from complex analysis to evaluate line integrals of analytic functions over closed curves. An important step is to write 8) 123 and to regard (2.8) as the sum of residues at z = z k . To construct an analytic function that, in conjunction with Cauchy's theorem and the closed curve |z| = 1 + , returns (2.8) and has singularities at z 1 , . . . , z g−1 leads us to consider the integrand in (2.4).
Here, the logarithmic function follows from (2.8) and the singularities with appropriate residues are created through . After careful consideration of the analytic properties of the integrand in (2.4), we then show that Cauchy's theorem gives (2.7), from which (2.4) follows. As said, (2.3) is obtained by manipulating (2.4), using partial integration. The formal proof of Theorem 1 presented in Sect. 4 contains several challenging steps and requires, among other things, a proof that the function Y (z)/z is injective in a region that contains the unit disk, and a way to account for the branch cut caused by the logarithm in (2.9) being taken over negative values.
Historical notes Integrals of this sort go a long way back in the history of queueing theory and were first found in the ground-breaking work of Pollaczek on the classical single-server queue (see [1,4,7] for historical accounts). Let us point out the connection to the well-known Pollaczek-type integral for the discrete bulk-service queue [9], a discrete-time queueing model in which customers upon arrival are placed in a queue and, after some stochastic service time, a (possibly stochastic) number of customers is served, all at once. Served customers leave the system immediately, whereas the remaining customers in the queue have to wait at least one more service period. The analysis of the bulk-service queue is easier than that of the FCTL queue. To see this, observe that the analysis of the FCTL queue would be greatly simplified if all vehicles were delayed [16], so that all vehicles arrive while the queue length is at least one and the complicated random variable A d n in (1.1) can be replaced by A n . In that case, (1.1) becomes a standard stochastic recursion driven by i.i.d. random variables and the FCTL queue reduces to the classical bulk-service queue, a special case of the more general single-server queue investigated by Pollaczek. Let X b denote the steady-state queue length in that bulk-service queue, defined as the solution of the stochastic equation (2.10) Pollaczek's result then says that (see [8,9] for a direct derivation)

Remark 2
The bulk-service queue serves as a popular approximation of the FCTL queue [16]. In fact, for Bernoulli arrivals with per time slot one or no arrival (which is case (i) in the proofs of Lemmas 3 and 4), this approximation becomes exact. To see this, substitute Y (z) = 1 − p + pz into the logarithmic function in (2.4), and observe that this gives the logarithmic function in (2.11). Obviously, when Y can take values larger than one, the bulk-service queue is an approximation and yields an upper bound on the overflow queue.

Generalized FCTL queues
Oblakova et al. [13] have introduced generalized FCTL queues and established contour integrals for the first moment of the steady-state queue length. We now show how contour-integral representations for these generalized FCTL queues follow almost directly from the standard FCTL queue. We start from the definition of X (z) in [13], a generalization of the function X g (z) that contains as special cases several extensions of the FCTL queue.
Definition 2 [13] Consider the function X (z) with X (1) = 1 and where B(z) and A(z) are PGFs and ξ(z) is a function satisfying ξ(1) = 0, ξ(z l ) = 0 with z l = 1 the roots of z g − A(z) inside the unit disk. Assume, moreover, that B (1) < 1; A (1) < g; that, for some δ > 0, the functions A(z) and B(z) are analytic within the disk |w| < 1 + δ; and that X (z) is analytic inside the unit disk and continuous up to the unit circle. Also assume that Here is the main result for the function X (z): Theorem 2 (Pollaczek integrals for generalized FCTL queues) Under Definition 2, there exists an 0 > 0 such that, for all ∈ (0, 0 ), (2.13) for all |z| < 1 + , with the principal value of the logarithm.
Proof We shall express (2.12) as a product of the PGF of the standard FCTL queue and some analytic function. Denote the g − 1 roots of z g − A(z) inside the unit circle by z 1 , . . . , z g−1 . We can then rewrite (2.12), using X (1) = 1, as This gives the result.
Let us now discuss the extensions contained in X (z).
(i) The first extension concerns right-turning traffic, in which case the difference in discharge rate between delayed and non-delayed vehicles almost vanishes. This requires modifying the FCTL assumption in order to put an upper bound on the number of vehicles that pass the traffic light without delay. This upper bound is set to one, so that at most one vehicle can depart per green slot. Following [13], it can be shown that this model for right-turning traffic follows by setting , and the contour-integral expression for the PGF thus follows from Theorem 2.
(ii) Another extension of the classical FCTL queue is one that accounts for disruptions of the traffic flow by, for example, pedestrians. To account for these disruptions, one could extend the red period or shorten the green period for the main stream of vehicles [13]. This extension thus requires a FCTL queue with random (but finite) green and red times, for which we choose g = G, with G denoting the maximum green time. Setting B(z) = Y (z), A(z) = r ,g p r ,g Y (z) r +g z G−g with p rg the probability that a cycle consists of g green and r red slots, and ξ(z) = z − Y (z), then shows that also this extension of interrupted flows is contained in Theorem 2.
(iii) The third extension we mention relates to uncertainty in departure times of vehicles, for instance, due to distracted drivers facing a green light and hence causing the driver to not depart from the queue with some probability p [13]. Assuming a geometrically distributed number of slots before a driver leaves the queue models this situation, with (iv) A fourth extension deals with relaxing the independence assumption on the arrival process during the red slots [13]. In this extension, the arrivals during a red time within a cycle may be dependent (the arrivals during green slots still need to be i.i.d.). For this FCTL queue, we should choose where A r (z) denotes the PGF of the arrival process during the whole red period, and The present paper adds to [13] the PGFs in terms of contour integrals, of which the contour integrals for the mean queue length obtained in [13] follow by evaluating the derivative at one. For insights into the differences between the various FCTL queue extensions, we refer to the elaborate numerical study in [13].

Algorithmic methods
We now discuss the computational challenges that come with calculating the steadystate queue-length distribution, using either the contour integrals in Theorems 1 and 2 or the standard expression in terms of roots. The algorithms using contour integrals in this section are based on representation (2.4) (but one could also take (2.3)). Notice that we only need to expand X g (w) at w = 0 and w = 1, so inside the validity range of (2.4).

From PGF to performance measures
The mean stationary overflow queue EX g is given by X g (1) and takes the form This result was recently obtained in [13] using a direct proof that converted the classical expression for EX g in terms of complex-valued roots into the integral expression (3.1). From the PGF X g (z), we can in principle determine all stationary moments. Define The moments E[X k g ] then follow from symbolically differentiating PGF (2.4), and these derivatives can be expressed as for k = 0, 1, 2, . . .. Using this recursive expression, X (k) g (w) can be expressed in terms of f (w) and the first k derivatives of f (w), denoted by f (1) and g ( j) (w, z) := ∂ j ∂w j g(w, z), for j = 1, 2, . . . , k. After substituting w = 1, we can express the first k moments of X g in terms of k contour integrals that only involve the 123 model primitives and the first k moments of Y . Using f (1) = 0, the variance of X g given by Var(X g ) = h 2 (1) + h 1 (1) − (h 1 (1)) 2 takes the form (3.7) To determine the stationary distribution of the overflow queue, we use that First observe that (3.9) Expressions for the other probabilities P(X g = k) follow in a similar way, but require evaluating the resulting function at w = 0 instead of w = 1 and dividing by k!. As a consequence, P(X g = k) can be expressed in terms of f (0), f (1) (0), . . . , f (k) (0), again an expression that involves explicit contour integrals only.

Roots or integrals?
Compared with root finding, contour integrals have advantages and disadvantages. On the one hand, avoiding the implicitly defined roots is nice, because the integrals are explicit expressions in terms of the model primitives g, r and Y (z). On the other hand, the number of terms required to evaluate f ( j) (w) grows exponentially in j. For tail probabilities, this symbolic differentiation becomes computationally cumbersome.
While in the early queueing literature root finding was considered to be prohibitively difficult, with the computational methods available nowadays, it is possible to find the complex-valued roots of z g − A(z) with great accuracy. In Appendix A, we present the root-finding algorithm that we use in this paper, which after extensive testing was found to be accurate and reliable for all choices of A(z). The simple idea behind the algorithm is to approximate A(z) with its Taylor series A n (z) of order n, reducing the problem to finding roots of polynomial equations, and subsequently use Newton's method to find the roots of z g − A(z) with arbitrary precision. We also present some results that show that the roots of the n-th system converge to the roots of z g − A(z) and provide an explicit characterization of the roots for the case when A(z) is the PGF of a Poisson random variable. In that case, the roots can be written in terms of the Lambert W-function.
Extensive tests with both algorithms did not result in any numerical issues, except for two obvious limitations: for tail probabilities the symbolic differentiation within the integrand becomes a bottleneck, and for root finding loss of accuracy is expected when the number of roots g becomes excessively large (although a thousand roots In terms of computation time, contour integration is generally slower than root finding. For moments there is little difference, because both methods lead to explicit expressions. For queue-length probabilities, however, numerical inversion of the PGF is required (see, for example, [2]), and in that case calculating the many contour integrals becomes computationally expensive.
To illustrate the algorithms, we now show some results for the FCTL queue with g = 20 and c = 50. We consider Poisson arrivals with on average λ vehicles arriving per slot, and four scenarios: λ = 0.2 (light traffic), λ = 0.3 (moderate traffic), λ = 0.36 (heavy traffic) and λ = 0.38 (extreme traffic). These arrival rates correspond to a volume/capacity ratio ρ = λc/g ranging from 0.5 to 0.95. The results are calculated with both roots and contour integrals and are on the scale of the displayed figures indistinguishable. Figure 1a shows the mean queue lengths E[X 0 ], . . . , E[X c ] through one cycle. Observe the strong cyclic behavior and the high sensitivity for ρ. Figure 1b shows the queue-length distribution at cycle start, the moment that the traffic signal turns green and queue lengths are expected to peak. Observe the difference between operating at 75% or 90% of maximal capacity: the probability that more than 20 vehicles are waiting is only 0.002 for λ = 0.3 and 0.32 for λ = 0.38. Figure 1c depicts the distribution of the effective green time G, defined in [16] as the number of slots used for departure of delayed vehicles that arrive throughout the whole cycle. We have (3.10) Since only one delayed vehicle departs per slot, this can also be considered to be the distribution of the platoon length consisting of delayed vehicles departing during one cycle. Observe that P(G = g) is practically zero when ρ = 0.5, but as high as 0.71 when ρ = 0.95, which means that in only 29% of the cycles the green time is long enough to let the queue vanish. Finally, we consider the delay distribution of an arbitrary vehicle arriving in the 10th slot, which is during the green period. The stationary delay of a vehicle arriving in slot k, denoted by D [k] , is defined as the number of slots between arrival and departure, not including the slot of arrival. Figure 1d shows D [10] , which can be computed directly from X 9 , i.e., the number of vehicles waiting at the start of the 10-th slot. If X 9 = 0, we have that D [10] = 0; otherwise, the delay can be expressed as a function of the number of vehicles present at the arrival of the tagged vehicle. This function (studied in detail in [17]) should take into account interruptions due to red periods, which explains the fragmented histograms in Fig. 1d.

Proof of the Pollaczek contour-integral representation
As explained briefly in Sect. 2, the proof of Theorem 1 exploits the factorized form (2.7) and investigates in detail the logarithmic function (2.9). We present some useful properties of the function Y (z)/z, visible in both (2.7) and (2.9). We then proceed to use Cauchy's theorem to obtain the contour-integral representation (2.4) for the case that 1 < w < 1 + , and finally manipulate (2.4) to obtain (2.3) on the full range |w| < 1 + .

Auxiliary results
Before we prove Theorem 1, we present some auxiliary results for the function Y (z)/z. In [13,Theorem 1], it was shown that the function Y (z)/z is injective on the disk |z| ≤ 1, so that all Y (z k )/z k = Y (z l )/z l when z k = z l . For our proof, we also need injectivity, but for the larger disk with radius t 0 > 1. More specifically, let where R is the maximum value such that Y (z) is analytic in the region |z| < R.

Fig. 2
The four components, C 1 , C w , L + and L − , of contour C As a consequence of Lemma 5, taking the principal value logarithm in (2.9) when 1 < w < 1+ < t 0 , we obtain a function of z that is analytic in the open disk |z| < t 0 , with branch cut [1, w].

Contour integral for (2.4)
We next consider the function z g − A(z) that has its zeros in |z| ≤ 1 at z = z 0 = 1, z 1 , . . . , z g−1 , while its other zeros have modulus greater than one. Let R 0 be the zero outside |z| ≤ 1 of the smallest modulus; this R 0 is real and larger than one. Take > 0 such that 1 + < min{t 0 , R 0 } and consider the integral Choose δ > 0 such that δ < 1 2 (w − 1) and δ < 1 + − w, while |z k − 1| > δ, k = 1, . . . , g − 1. Now let C be the positively oriented contour consisting of the circles C 1 (δ) and C w (δ) of radii δ around 1 and w, respectively, together with the line segments Fig. 2 for the positioning of the contour C with its four components in the disk |z| < 1 + and relative to the zeros of z g − A(z). Then, by Cauchy's theorem, On the line segments z = t ± i0, 1 + δ ≤ t ≤ w − δ, we use that (4.19) 123 With the principal value choice for ln, we then get (with Therefore, also using that t g − A(t) > 0, 1 < t < 1 + , where we have used that As to the last integral on the last line of (4.21), we use that with nonvanishing numbers wY (w) − Y (w), Y (w) − w and w g − A(w). Therefore, The middle integral on the last line of (4.21) is more delicate, since both Y (z) − z and z g − A(z) vanish at z = 1. For z = 1 + δe iφ , with 0 < φ < 2π and δ ↓ 0, Next, as z → 1, since g − A (1) > 0. Hence, from (4.29) and (4.30) with z = 1 + δe iφ and dz = iδe iφ dφ in the integral over C 1 , where we have also used that 2π 0 (π − φ) dφ = 0. Using (4.22), (4.27) and (4.31) in (4.21), we get Returning then to (4.17)-(4.18), letting δ ↓ 0, we see that (4.33) by (2.7). Here we have also used that the zeros z k are real or come in conjugate pairs so that, for w ∈ (1, 1 + ) by (4.11) in Lemma 5, both X g (w) and the product g−1 k=1 in (4.33) are real and positive, with of Y . Since X g (w) = exp(I (w)) for 1 < w < 1 + , we then get by analyticity of X g that (4.44) holds for all w, |w| ≤ 1 + 1 and any 1 ∈ (0, ). Then a simple rearrangement of the integrand in (4.44) yields Theorem 1.

Conclusions
We have presented novel formal solutions for the FCTL queue in the form of contour integrals. Theorem 1 presents the contour-integral representation for the PGF of the overflow queue. From this PGF, essentially all relevant information about the stationary behavior of the FCTL queue can be obtained, by taking derivatives at one for the moments, derivatives at zero for the distribution, and by using simple recursions to obtain the queue lengths at all moments within the cycle and the stationary delay distribution. A contour-integral expression for the first moment was obtained recently by Oblakova et al. [13], and the present paper can be seen as an extension of that work. The two papers together present an alternative approach for the FCTL queue and its generalizations, using contour integrals instead of factorizations in terms of complex roots that need to be determined numerically.
A possible thread for future research relates to asymptotics. In classical queueing theory, a prominent line of research is related to heavy traffic, an asymptotic regime in which the traffic intensity approaches 100%. Next to more probabilistic methods such as weak convergence techniques and coupling, another way to obtain heavytraffic results is through the asymptotic evaluation of Pollaczek-type integrals; see, for example, [4,10] for single-server queues and [9] for classical bulk-service queues. Now that Pollaczek-type integrals for the FCTL queue are available, it is worthwhile to explore the possibilities for heavy-traffic analysis.
A(z) (which typically is a non-polynomial function) with its Taylor series A n (z) of order n. Solving this truncated equation boils down to root-finding of a polynomial. If the roots of the truncated equation are sufficiently close to the roots of z g − A(z), we can find the latter roots easily from the roots of z g − A n (z) by using a Newton-Raphson-type method.
We now present two propositions, in support of Algorithm 1. The propositions state that under very mild conditions the number of roots of z g − A(z) on or within the unit circle is equal to the number of roots of the truncated equation z g − A(z) on or within the unit circle and that the roots of the truncated equation converge to the roots of z g − A(z) (when n tends to infinity).

Proposition 6
Let D(z) = z g − A(z) and let D n (z) := z g − A n (z), where A n (z) denotes the n-th order Taylor approximation of A(z). Upon assuming that A(z) is a PGF; that A(z) is analytic in the disk |z| < 1 + δ for some δ > 0; and that g < A (1), D n (z) = 0 has as many roots on or within the unit circle as D(z) (i.e., g).

Proof
Rouché's theorem says that if f and g are analytic inside some region K with closed contour ∂ K and if |g(z)| < | f (z)| on ∂ K , then f and f + g have the same number of zeros inside K .
The conditions that A(z) has to be analytic in |z| < 1 + δ and that g < A (1) together imply (1 + γ ) g > A(1 + γ ), (A.1) for some γ ∈ (0, δ); see, for example, [5]. Assume |z| = 1 + γ . Then where the strict inequality follows from (A.1) and the remaining inequalities from the fact that A(z) is a PGF. So we may apply Rouché's theorem on f (z) = z g and g(z) = −A n (z). As z g has g roots on or within the unit circle, D n (z) will have g roots as well, just as D(z). because a i ≥ 0 and |z j | ≤ 1.
From Proposition 7, we see that if we let n tend to infinity, then D n (z j ) tends to 0. This implies that the roots obtained by using D n (z) will be close to the actual roots of D(z) when n is sufficiently high.