ANALYSIS OF A SPLITTING PROCESS ARISING IN PROBABILISTIC

COUNTING AND OTHER RELATED ALGORITHMS ^{∗}

May 10, 1996

Peter Kirschenhofer  Helmut Prodinger  Wojciech Szpankowski^{†}
Dept. Algebra&Discrete Math. Dept. Algebra&Discrete Math. Dept. Computer Science
Technical University of Vienna Technical University of Vienna Purdue University
A-1040 Vienna  A-1040 Vienna  W. Lafayette, IN 47907
Austria  Austria  U.S.A.

Abstract

We present an analytical method of analyzing a class of “splitting algorithms” that include probabilistic counting, selecting the leader, estimating the number of questions necessary to identify distinct objects, searching algorithms based on digital tries, approximate counting, and so forth. In our discussion we concentrate on the analysis of a generalized probabilistic counting algorithm. Our technique belongs to the toolkit of the analytical analysis of algorithms, and it involves solutions of functional equations, analytical poissonization and depoissonization, Mellin transform, etc. In particular, we deal with a functional equation of the form g(z) = βa(z)g(z/2) + b(z) where a(z) and b(z) are given functions, and β < 1 is a constant. With respect to our generalized probabilistic counting algorithm, we obtain asymptotic expansions of the first two moments of an estimate of the cardinality of a set that is computed by the algorithm. We also derive the asymptotic distribution of this estimate, and observe that it actually fluctuates, leading to a conclusion that its limiting distribution does not exist.

∗A preliminary version of parts of this paper was presented in ICALP 92, Vienna [26]. This work was supported by NSF Grants INT-8912631 and CCR-9201078.

†This author was additionally was supported by NSF Grants NSF Grants NCR-9206315 and NCR-9415491, and in part by NATO Collaborative Grant CRG.950060.

1. INTRODUCTION

In several database algorithms a major determinant of efficiency is the cardinality of the underlying set. Therefore, one needs an efficient way of estimating the cardinality of large (multi)sets of data. Clearly, the trivial method of building a list of elements without replications is unacceptable due to the cost of disk access and auxiliary memory. Knowing this fact, Flajolet and Martin [9] proposed an algorithm that is probabilistic in its nature. It works as follows: In order to estimate the cardinality N of a set M (with replications) every element x ∈ M is hashed into a binary string of size m (the choice of m is easy, and m = 5 + log N suffices [4]). The bitwise OR-composition of the modified hashed strings in which only the least significant (rightmost) “1” survives is used to build the so called bitmap and to obtain the estimate R_{N} of the cardinality N of M. More precisely, the position R_{N} of the leftmost zero in the bitmap string (where we start enumerating positions with 0) approximates log_{2} N as shown by Flajolet and Martin in [9] (cf. also [4]). Observe that, equivalently, R_{N} is the length of the longest run of ones at the beginning of the bitmap string.

In fact, probabilistic counting is a special case of a splitting process that also includes such algorithms as selecting the loser [7, 31], estimating the number of questions necessary to identify distinct objects [30], searching algorithms based on digital tries [16, 17, 18, 26, 28, 32, 36], approximate counting [24], conflict resolution algorithms for multiaccess communication [17, 34, 35] and so forth. Using a digital tree representation one can describe the probabilistic counting in terms of this splitting process as follows: Imagine N persons flipping a fair coin, and those who get 1 (a hit) discontinue throwing (move to the right in the digital tree representation) while the rest continue throwing (i.e., move to the left in the digital tree) until they get a hit. The process continues until all remaining persons flip a 0. It should be clear that the number of rounds in the flipping procedure equals one plus the quantity R_{N} from above.

Let us now consider the following more general version of the splitting process, where the coin flipping ends as soon as at most d persons have flipped a 1 in the same round, where d is a given parameter. This time we denote the number of rounds by 1 + R_{N,d}. Observe that d = 0 corresponds to the original situation. Again we will expect R_{N,d} to be an estimate of log_{2} N. We should point out that a similar extension of the primary leader election algorithm was recently considered by Grabner [13].

Returning to the probabilistic counting algorithm: We can interpret the generalized situation as follows: We start from an empty bitmap string, that is, with all positions filled by zeros. We further assume that N different objects (e.g., data, persons, etc.) can randomly insert (hit) a 1 at any position of the bitmap, however, the probability of hitting the j ≥ 0 position is equal to 2^{−j−1}. (In terms of the modified hashed strings discussed above this means that the probability of the occurrence of the pattern like 0^{j}1 ⋯ is equal to 2^{−j−1} since 0 and 1 are equally likely.) Every object can hit only one time. We count the number of hits in any position of the bitmap, but we count the number of hits only up to the value d + 1, where d is a given parameter. In other words, the bitmap is a d + 2-ary string now. Clearly, the bitmap will contain many d + 1’s at the beginning of the string. This is due to the fact that the probability of a hit decreases exponentially fast with the increase of the position in the bitmap. It is easy to see that the length of the longest run of d + 1’s in the front of the bitmap equals R_{N,d} defined above. More precisely, we have:

R_{N,d} = min{k : bitmap(k) < d + 1 and for all 0 ≤ i < k bitmap(i) = d + 1} . (1)

Observe that this time replications in the object set are not allowable since they might influence the result of the bitmap procedure if d ≥ 1.

For example, for d = 3 we may have bitmap = 4444444444434100000. As the estimate R_{N,d} we take the position of the leftmost value smaller than d + 1 (e.g., in the example above we have R_{N,3} = 11).

In the present paper we analyze the probabilistic behavior of the parameter R_{N,d}. In particular, we prove that ER_{N,d} = log_{2} N − 1 − C_{d} + P (log_{2} N) + O(\sqrt{N}) where C_{d} is a constant such that C_{d} = log_{2} d − \frac{1}{2} + O(1/d) as d → ∞, and P (x) is a periodic function of period 1 and small amplitude (cf. Theorem 1(i)). More importantly, we show that the variance of R_{N,d} fulfills (neglecting periodic fluctuations of mean zero) Var R_{N,d} ∼ \frac{1}{\sqrt{πd}} as d → ∞ (cf. Theorem 1(ii). Finally, we derive the asymptotic distribution of R_{N,d} which confirms our believes that the estimate R_{N,d} is well concentrated around its mean (cf. Theorem 2). But, we observe that the limiting distribution of R_{N,d} does not exist due to a tiny oscillation of the asymptotic distribution. A preliminary version of some of these results can be found in [26].

It must be stressed that the novelity of our paper also lies in the area of the mathematical analysis of algorithms, and we think of the splitting process (in particular, the probabilistic counting) as a motivation for a more general study. In fact, the proposed technique of analysis can be applied to several other algorithms. As indicated above, the study of the splitting process can be reduced to a problem on digital trees. Evaluating a parameter of such a tree, we obtain the following functional equation

g(z) = βa(z/2)g(z/2) + b(z) .  (2)

In general, this equation does not have an explicit finite solution unless a(z) is a very particular function. For example, a(z) = e^{z} was discussed in Knuth [28] (cf. also [36]), Flajolet and Saheb [10], and Szpankowski [35], [36] (cf. also [28]) investigated the case a(z) = e^{z} + 1. For these two cases explicit finite solutions of (2) were obtained, but it is rather hopeless to expect the existence of a useful closed form solution for general a(z) and b(z). An asymptotic solution to (2) was suggested in Jacquet and Szpankowski [17] who analyzed the case a(z) = 1 + (1 + z)e^{−z}. More recently, Flajolet and Richmond [11] investigated asymptotically the above equation with a(z) = 2z^{b}/(1 + z)^{b} for some integer b > 0. In the analysis of the algorithms introduced in this paper, we must investigate among others the equation (2) with a(z) = 1 − e_{d}(z)e^{−z} where e_{d}(z) = ∑_{i=0}^{d} z^{i}/i!. Even more complicated functions for a(z) and b(z) are involved in the analysis of the variance and the asymptotic distribution (cf. Sections 2 and 3).

As mentioned above, the case d = 0 was analyzed by Flajolet and Martin [9]. In addition, Greenberg et al. [14] investigated the instance d = 1 in the analysis of a conflict resolution algorithm for multiple channels. It must be stressed that the analysis of Flajolet and Martin [9] is based on the inclusion-exclusion rule, and seems to be very unlikely extendable to d ≥ 1. This is partially evident from the analysis in [14].

The plan for the present paper is as follows: In the next section we present our main results concerning the generalized counting algorithm and the associated splitting process, and their algorithmic consequences. The proofs are delayed till Section 3 where – besides of proving our findings – we also present a general methodology of dealing with functional equations like (2). In particular, we use a depoissonization lemma that allows us to obtain an estimate of a sequence from its exponential generating function. This approach turned to be very useful in many other situations encountered in the analysis of algorithms (cf. [7, 12, 17, 33, 32]).

2. MAIN RESULTS

Let F_{N}(u) = Eu^{R_{N,d}} be the probability generating function of the estimator R_{N,d}. That is, [u^{j}]F_{N}(u) (the jth coefficient at F_{N}(u)) is the probability that either the splitting process terminates after j + 1 rounds or that the result of the bitmap is (d + 1)^{j}c … where c ≤ d, and (d + 1)^{j} denotes j consecutive values “d + 1”. To derive the recurrence for F_{N}(u) we observe: If at the first step k persons flip 1, then either the contribution to the generating function F_{N}(u) is u^{0} ∑_{k≤d} \binom{N}{k}2^{−N} for k ≤ d, or the contribution becomes u ∑_{d+1≤k≤N} \binom{N}{k}2^{−N}F_{N−k}(u) for k > d. In short, we obtain the following recurrence

F_{N}(u) = ∑_{k=0}^{d} \binom{N}{k}2^{−N} + u ∑_{k=0}^{N−d−1} \binom{N}{k}2^{−N}F_{k}(u) ,  (3)

which is valid for all N ≥ 0.

Using (3) we can derive recurrences for the first and second factorial moments of R_{N,d}, namely: ER_{N,d} = F′_{N}(1) and F′′_{N}(1). Observe that Var R_{N,d} = F′′_{N}(1) + F′_{N}(1) − (F′_{N}(1))^{2}. In particular, let L(z) and W (z) be the exponential generating function for F′_{N}(1) and F′′_{N}(1), respectively. That is,

L(z) = ∑_{N≥0} F′_{N}(1)z^{N}/N! ,
W (z) = ∑_{N≥0} F′′_{N}(1)z^{N}/N! .

Define also L̃(z) = e^{−z}L(z) and W̃(z) = e^{−z}W (z) as the Poisson generating functions. These Poisson generating functions represent the first two moments in a variation of the probabilistic counting problem in which the number of people (objects) is not fixed but is a random variable distributed according to a Poisson process with mean z. Then, after some algebra (3) implies

L̃(z) = f_{d}(z/2)L̃(z/2) + f_{d}(z/2) ,  (4)
W̃(z) = f_{d}(z/2)W̃(z/2) + 2L̃(z)f_{d}(z/2)  (5)

where

f_{d}(z) = 1 − e_{d}(z)e^{−z}  and  e_{d}(z) = 1 + \frac{z^{1}}{1!} + ⋯ + \frac{z^{d}}{d!} .  (6)

The above functional equations are solved in the next section (cf. Section 3.1). Their solutions are expressed in terms of the following function that is needed to articulate our main results below:

φ(x) = ∏_{j=0}^{∞} f_{d}(x2^{j}) = ∏_{j=0}^{∞} (1 − e_{d}(x2^{j})e^{−x2^{j}} ) .  (7)

We now can present our first main result which is proved in Section 3.2.

Theorem 1. Consider the generalized splitting algorithms described above. Then the estimator R_{N,d} behaves as follows:

(i) The average value of R_{N,d} becomes asymptotically as N → ∞

ER_{N,d} = log_{2} N − 1 − C_{d} + P_{1}(log_{2} N) + O (N^{−1/2+ε})  (8)

for any ε > 0 with

C_{d} = \frac{1}{L^{2}} ∫_{0}^{∞} e^{−x}e_{d}(x)φ(2x)\frac{log x}{x}dx
= log_{2} d − \frac{1}{2} + O (\frac{1}{d})  as d → ∞  (9)

where L = log 2 and the function φ(x) is defined in (7). Furthermore, with χ_{k} = 2kπi/L,

P_{1}(x) = −\frac{1}{L^{2}} ∑_{k≠0} Φ^{∗′}(−χ_{k})e^{2πikx}  (10)

where Φ^{∗}(s) is the Mellin transform of φ(2x) − φ(x) and Φ^{∗′}(s) is the first derivative of Φ(s) (cf. Section 3 for details). In particular, P_{1}(x) is a periodic function of mean zero and small amplitude for reasonable values of d.

(ii) The variance of R_{N,d} is

Var R_{N,d} = D_{d} − C_{d}^{2} − [P_{1}^{2}]_{0} + P_{2}(log_{2} N) + O (N^{−1/2+ε})  (11)

for any ε > 0, where

D_{d} − C_{d}^{2} − [P_{1}^{2}]_{0} = \frac{1}{\sqrt{πd}} + O (\frac{1}{d}) for  d → ∞ .  (12)

More precisely

D_{d} = \frac{1}{L^{3}} ∫_{0}^{∞} e_{d}(x)e^{−x}φ(2x)\frac{log^{2} x}{x}dx
= −C_{d} − 1/6 + log_{2}^{2} d + O (\frac{1}{d}) as d → ∞ ,  (13)

so that

D_{d} − C_{d}^{2} = \frac{1}{12} + O (\frac{1}{d})  as d → ∞ .  (14)

The quantity [P_{1}^{2}]_{0} is the (integral) mean of the square of P_{1}(x), and for large d behaves as

[P_{1}^{2}]_{0} = \frac{1}{12} − \frac{1}{\sqrt{πd}} + O (\frac{1}{d})  as d → ∞ .  (15)

Finally, P_{2}(x) is a periodic function with period 1 and small amplitude.

Using numerical integration, we computed for some values of d the constant terms (for N → ∞) in ER_{N,d} − log_{2} N − P_{1}(log_{2} N) and Var R_{N,d} − P_{2}(log_{2} N). They are displayed in Table 1. We note that the variance initially rapidly decreases (for d ≤ 2) and then the decrease slows down.

Table 1: Performance of the algorithm for small values of d (neglecting fluctuations of mean zero).

\begin{matrix}d & ER_{N,d} − log_{2} N & Var R_{N,d} \\ 0 & -0.37 & 1.26 \\ 1 & -1.40 & 0.78 \\ 2 & -2.00 & 0.59 \\ 3 & -2.42 & 0.49 \\ 4 & -2.74 & 0.43 \\ 5 & -3.02 & 0.38 \\ 6 & -3.24 & 0.35 \\ 7 & -3.43 & 0.32 \end{matrix}

As an immediate consequence of Theorem 1, we observe that R_{N,d} converges in probability (pr.) to ER_{N,d}. Indeed, by Chebyshev’s inequality we have

Pr\{|R_{N,d}/ER_{N,d} − 1| > ε\} ≤ \frac{Var R_{N,d}}{ε^{2}E^{2}R_{N,d}} = O (\frac{1}{log^{2} N} ) → 0

as N → ∞. This proves the announced convergence, but does not warrant the almost sure (a.s.) convergence since the Borel-Cantelli Lemma cannot be applied. However, observe that R_{N,d} is a nondecreasing function of N, in the sense that on every sample path we have R_{N+1,d} ≥ R_{N,d}. To prove the (a.s.) convergence one should apply the Borel-Cantelli Lemma along a doubly exponentially increasing sequence such as N = 2^{2^{k}} (cf. [21]), thus R_{N,d} → ER_{N,d} (a.s.), too.

Our technique also allows to deal with the asymptotic distribution for R_{N,d} and the limiting generating function. For this, we introduce the Poisson transform of F_{N}(z) as

G̃(z, u) = e^{−z} ∑_{N=0}^{∞} F_{N}(u)\frac{z^{N}}{N!} .  (16)

Observe that G̃(z, u) can be interpreted as the generating function of R_{N,d} when the cardinality of a set is Poisson distributed with mean z.

The Poisson generating function G̃(z, u) satisfies the following functional equation

G̃(z, u) = uf_{d}(z/2)G̃(z/2, u) + (u − 1)(f_{d}(z/2) − 1) .  (17)

Using the technique from Section 3.1, we can solve it for z → ∞, and then by depoissonization (that is, after recovering F_{N}(u) from G̃(z, u)) we prove in Section 3.3 our second main finding.

Theorem 2. Let t be a fixed real number. Then, for large N the asymptotic distribution of R_{N,d} becomes for N → ∞ (d fixed)

Pr{R_{N,d} ≤ log_{2} N + t − 1} = 1 − φ (2^{−t−{log_{2} N+t}}) + O(N^{−1/2+ε}) ,  (18)

for any ε > 0, where φ(x) is defined in (7), and {log_{2} N +t} = log_{2} N +t−⌊log_{2} N +t⌋. Since {log_{2} N + t} is dense in [0, 1] but not uniformly distributed in [0, 1] the limiting distribution of R_{N,d} does not exist. In fact,

lim inf_{N→∞} Pr{R_{N,d} ≤ log_{2} N + t − 1} = 1 − φ (2^{−t−1})  (19)
lim sup_{N→∞} Pr{R_{N,d} ≤ log_{2} N + t − 1} = 1 − φ (2^{−t})  (20)

for real t and all d ≥ 0.

Remark. With a more involved proof, using the complex Mellin transform it can be shown that the remainder term in (18) is O(N^{−1/2}), and the error terms in (8) resp. (11) are O(log N/N) resp. O(log^{2} N/N). □

The nonexistence of the limiting distribution could be expected. Observe that the estimate R_{N,d} is related to an extreme statistics of a set of N discrete random variables (in fact, geometrically distributed). Anderson [3] in 1970 predicted such a behavior for M_{n} = max{X_{1}, …, X_{N}} where X_{i} are i.i.d. discrete random variables. In fact, when X_{i} are geometrically distributed with parameter 1/2, then the asymptotic distribution of M_{n} − log_{2} N oscillates between two double exponential distributions, namely e^{−2^{−t−1}} and e^{−2^{−t}}. The bounds (19) and (20) resemble the double exponential distribution. Indeed, consider d = 0, and take only the first two terms of the infinite product φ(2^{−t}). Then, φ(2^{−t}) ≈ ∑_{j≥0} e^{−2^{−t+j}}. This is illustrated in Figure 1 where we show the upper “envelope” for the asymptotic distribution of R_{N,d}−log_{2} N (i.e., 1−φ(2^{−t})) for d = 0 and d = 3). In passing, we should point out that recently Fill et. al. [7] used the depoissonization technique to obtain a similar oscillation in the height of the incomplete trie. A thorough account on analytical depoissonization can be found in Jacquet and Szpankowski [19]. Our technique might be applied to analyze several other algorithms that are based on the “splitting process” described at the beginning of this section. More precisely, a splitting algorithm divides randomly N objects according to some rule. For example, consider a multiaccess system with a large number of users. The following “conflict resolution” algorithm found to be successful for some multiaccess systems (cf. [14], [17]). Let N be the size of a conflict, that is, the number of users that simultaneously sent packets. To split it, every user flips a fair coin, and those who get 1 “move right” in a digital tree represented this process, while the others “move left”. This splitting process continues until all users sent successfully their packets (i.e., subtrees containing conflicting users are of size one). This algorithm can be analyzed in a unified manner by our approach through the functional equation of type (21).

Here is another example that generalizes PATRICIA tries (cf. [28], [36]). Consider N infinite strings that are built over a V -ary alphabet. Using a trie construction as in Knuth [28] (see also [36]) we can split all strings in such a manner that no two of them will share the common prefix. A trie built in such a manner may contain some internal nodes that are unary. To avoid such a waste of storage, one may compress this trie to obtain a PATRICIA trie. In the PATRICIA all internal nodes have degree greater or equal to two. But a generalization of a PATRICIA trie might be also explored by imposing that all internal nodes have degree of size greater or equal to d ≤ V. What is the search time (i.e., the length of a path to an external node) in such a trie? What about the average height, etc.? The analysis of such a data structure seems to be a nontrivial task, and our approach can be applied to solve some of these problems.

3. ANALYSIS

In this section we present the proof of Theorem 1 (cf. Section 3.2) and Theorem 2 (cf. Section 3.3). To streamline our analysis, in Section 3.1 we discuss a general solution of a functional equation that arises on many occasions in this paper (cf. (2)), and in general in the analysis of algorithms (cf. [4, 7, 10, 11, 14, 22, 24, 25, 28, 31, 32, 33, 34, 35, 36, 37]). In addition, we offer a unified approach to depoissonization: a technique that allows to analyze the Poisson model instead of the more difficult Bernoulli model (i.e., when the number of objects is fixed). Poissonization is a standard probabilistic technique (cf. Aldous [2]), however, depoissonization causes usually some problems. We cope with it in Section 3.1 where we present an analytical approach to the depoissonization (cf. also [7, 17, 18, 19, 32, 33]).

3.1 General Solution and Depoissonization

The functional equations (4) and (5) for the Poisson generating functions of L̃(z) and W̃(z), as well as (17) on G̃(z, u) satisfy the following general functional equation

g̃(z, u) = βa(z/2, u)g̃(z/2, u) + b(z, u) ,  (21)

where |β| ≤ 1, and u is either fixed (cf. (4) and (5)) and in this case we simplify the notation to g̃(z), or u belongs to a compact neighborhood U(u_{0}) of u_{0} (cf. in (17) we have u_{0} = 0). The factor β may depend on u. Iterating (21), we obtain a general solution in the form of (cf. [17, 35])

g̃(z, u) = ∑_{n=0}^{∞} β^{n}b(z2^{−n}, u) ∏_{k=1}^{n} a(z2^{−k}, u)  (22)

if g̃(0, u) = 0. Define

φ(z, u) = ∏_{j=0}^{∞} a(z2^{j}, u) ,  (23)

provided the infinite product in (23) converges. Then, the general solution (22) can be rewritten as

g̃(z, u) = ∑_{n=0}^{∞} β^{n}b(z2^{−n}, u) ∏_{j=0}^{n−1} φ(z2^{−j}, u) .  (24)

or even more conveniently as

g̃(z, u)φ(z, u) = ∑_{n=0}^{∞} β^{n}b(z2^{−n}, u)φ(z2^{−n}, u) .  (25)

The last form is the clue to our asymptotic solution of the functional equation (22), and in particular to solve (4), (5) and (17). This is due to the fact that the sum in (25) falls under the so called harmonic sum that can be easily handled by the Mellin transform technique (cf. [12, 34]): Using the Mellin transform we obtain the asymptotic behavior of the above function for z → ∞, and then we “depoissonize” it to recover the original sequence.

In view of this, the important part of our analysis relies on depoissonization that we discuss next. To assure enough generality, let g_{N}(k) be a double sequence of N ≥ 0 and k ≥ 0 (e.g., g_{N}(k) = Pr{R_{N,d} = k}). Define g_{N}(u) = ∑_{k≥0} u^{k}g_{N}(k), the exponential generating function g(z, u) = ∑_{N≥0} g_{N}(u)z^{N}/N!, and the Poisson generating function g̃(z, u) = e^{−z}g(z, u). When the sequence depends only on N (e.g., g_{N} = ER_{N,d}), we write g̃(z) = e^{−z} ∑_{N≥0} g_{N}z^{N}/N!. On several instances the Poisson generating function g̃(z, u) (or g̃(z)) is easier to handle, and one can obtain the asymptotic behavior of g̃(z, u) for z → ∞. The question is how to recover g_{N}(u) (or g_{N}). This reverse process is called depoissonization, and the lemma below was proved in Rais et. al. [32] (cf. also [16, 18, 33]).

Depoissonization Lemma. Let S_{θ} be a cone S_{θ} = {z : |argz| < θ , 0 < θ < π/2}. As z → ∞ the following two conditions are assumed to hold for all u in a compact set U:

(I) For z ∈ S_{θ}

g̃(z, u) = O(|z|^{ε})  (26)

for some ε > 0; (O) For z ∉ S_{θ}

g̃(z, u)e^{z} = O(e^{α|z|})  (27)

for some α < 1.

Then for large N uniformly in u ∈ U the generating function g_{N}(u) = [z^{N}/N!]g̃(z, u)e^{z} satisfies

g_{N}(u) = g̃(N, u) + O(N^{−1/2+ε}) .  (28)

Remark (i) The error term in (28) depends only on the bounds on g̃(z, u) (cf. (26) and (27)) and does not depend directly on the function g̃(z, u) itself.

(ii) A more careful analysis reveals that the error term in (28) is O(N^{−1+ε}) (cf. [19]). □

Using the Depoissonization Lemma, we can “depoissonize” functional equation (21), and prove the following general result.

Theorem 3. Consider the functional equation (21) for u in a compact set U such that g̃(z, u) is bounded in U. Choose large B > 0 such that for |z| > B and any constants β_{1}, β_{2}, δ, δ_{1} > 0 the following conditions hold: (I) or z ∈ S_{θ}

\frac{|βa(z/2, u)|}{2^{ε}} < 1 − δ , |b(z, u)| < β_{1}δ|z|^{ε} ;  (29)

(O) for z ∉ S_{θ}

|βa(z/2, u)e^{z/2}| < (1 − δ_{1})e^{α|z/2|} , |b(z, u)e^{z}| < β_{2}δ_{1}e^{α|z|} ,  (30)

for α < 1. Then, q

g_{N}(u) = g̃(N, u) + O(N^{−1/2+ε})  (31)

for u ∈ U.

Proof. It suffices to prove that under (29) and (30) the two conditions (26) and (27) of the Depoissonization Lemma hold. The proof is by induction along the so called increasing domains as used in [16] and defined precisely in [18]: Suppose 1 < λ < 2 is a given real number, and then define increasing domains D_{m} for m = 0, 1, … as

D_{m} = {z : B ≤ |z| < Bλ^{m+1}} .

Note that if z ∈ D_{m+1} − D_{m}, then z/2 ∈ D_{m}, for m = 0, 1, …. Thus, we can apply induction with respect to m.

We first show that (26) holds provided inequalities (29) are fulfilled for z ∈ S_{θ}. Define D̂ = D_{m} ∩ S_{θ}. Clearly (26) holds in D̂_{0} as long as a solution of the functional equation is bounded in a compact set. Let us now assume that condition (26) is satisfied in D̂_{m}, and we prove that it also holds in D̂_{m+1}. If z ∈ D̂_{m}, then we are done. If z ∈ D̂_{m+1} − D̂_{m}, then z/2 ∈ D̂_{m}, and we can apply induction hypothesis, that is, |g(z/2, u)| ≤ β_{1}|z|^{ε}/2^{ε} for large z ∈ S_{θ}. Then, from the induction hypothesis, inequalities (29), and equation (21) we have

|g̃(z, u)| ≤ β_{1}(1 − δ)|z|^{ε} + β_{1}δ|z|^{ε} = β_{1}|z|^{ε} .

Thus, (26) holds in the larger domain D_{m+1}, and by induction, in the whole cone S_{θ}.

In a similar manner we prove that the second set of inequalities (30) imply (27) by noting that |e^{z}| = e^{ℜ(z)} ≤ e^{α|z|} where α = cos θ < 1. This time we use domains \overline{D}_{m} = D_{m} ∩ \overline{S}_{θ} where \overline{S}_{θ} is the complimentary set to S_{θ}. The rest follows the same line of arguments as above. This completes the proof of Theorem 3.

3.2 Analysis of Moments

We now prove Theorem 1. We consider two cases: bounded d, and large d (d → ∞).

Case A: Bounded d

We start with the average value ER_{N,d}. Its Poisson generating function L̃(z) satisfied the functional equation (4) which we repeat below

L̃(z) = f_{d}(z/2)L̃(z/2) + f_{d}(z/2)  (32)

where f_{d}(z) = 1 − e_{d}(z)e^{−z} and e_{d}(z) = 1 + z^{1}/1! + ⋯ + z^{d}/d!. This equation falls under the general functional equation (21), thus by (25) we have

L̃(z)φ(z) = ∑_{n=0}^{∞} φ(z2^{−n−1}) ,  (33)

where

φ(z) = ∏_{j=0}^{∞} f_{d}(z2^{j}) = ∏_{j=0}^{∞} (1 − e_{d}(z2^{j})e^{−z2^{j}} ) .  (34)

Observe that, by d’Alembert’s criterion, the final expression converges absolutely for all complex numbers z.

The idea of our analysis is to obtain the asymptotic behavior of L̃(z) as z → ∞, and then to apply depoissonization. To obtain the asymptotics of L̃(z) we make use of the Mellin transform technique. The Mellin approach is by now a standard technique in the analysis of algorithms. The interested reader can find more on Mellin transform in a recent survey [12] (cf. also [34]). For the reader’s convenience, however, we provide below some properties of the transform that we shall use in our analysis.

The following properties are used in this paper: (P1) Definition and Inverse of the Real Mellin Transform

f^{∗}(s) = ∫_{0}^{∞} f(x)x^{s−1}dx , f(x) = \frac{1}{2πi} ∫_{c−i∞}^{c+i∞} f^{∗}(s)x^{−s}ds  (35)

where c belongs to the fundamental strip defined below.

(P2) Fundamental Strip The Mellin transform of f(x) exists in the fundamental strip ℜ(s) ∈< −α, −β > where

f(x) = \{ \begin{matrix}O(x^{α}) & x → 0 \\ O(x^{β}) & x → ∞ . \end{matrix}

(P3) Harmonic sum

f(x) = ∑_{k≥0} λ_{k}g(μ_{k}x) ⇔ f^{∗}(s) = g^{∗}(s) ∑_{k≥0} λ_{k}μ_{k}^{−s}  (36)

provided the above sum converges. In the above, we used the following property

f(x) = g(ax) ⇔ f^{∗}(s) = a^{−s}g^{∗}(s)  (37)

in the fundamental strip of f^{∗}(s).

(P4) Asymptotic Expansion If f^{∗}(s) satisfies certain smallness-conditions towards i∞ (cf. [12]) for −β < ℜ(s) ≤ M, (M > 0) and

f^{∗}(s) = ∑_{k=0}^{K} \frac{d_{k}}{(s − b)^{k+1}} ,  (38)

then as x → ∞

f(x) = ∑_{k=0}^{K} \frac{d_{k}}{k!} x^{−b}(− log x)^{k} + O(x^{−M}) .  (39)

(P5) Mellin Transform in Complex Plane (cf. [5, 6, 8]) If F (z) is analytic in a cone θ_{1} ≤ arg(z) ≤ θ_{2} with θ_{1} < 0 < θ_{2}, then the Mellin transform F^{∗}(s) can be defined by replacing the path of integration [0, ∞[ by any curve starting at z = 0 and going to ∞ inside the cone, and it is identical with the real transform f^{∗}(s) of f(z) = F (z)|_{z∈R} from (35). In particular, if f^{∗}(s) fulfills an asymptotic expansion (38), then (39) holds for F (z) as z → ∞ in the cone.

Now we are ready to solve asymptotically (33) since we recognize it as a harmonic sum. Let Q_{1}(z) = L̃(z)φ(z) and Q_{1}^{∗}(s) be its Mellin transform which due to (P2) exists in the strip ℜ(s) ∈ (−∞, 0). By the harmonic sum formula (P3) we have Q_{1}(s) = 2^{s}/(1−2^{s})φ^{∗}(s) where the Mellin φ^{∗}(s) of φ(x) exists in ℜ(s) ∈ (−∞, 0). There is a problem, however, since the Mellin transform of φ(s) does not exists at s = 0, and we need a more precise estimate of φ^{∗}(s) at s = 0.

Thus, we proceed as follows. Define

Φ(x) = φ(2x) − φ(x) = e_{d}(x)e^{−x}φ(2x) .  (40)

Observe now that the Mellin transform of Φ(x) exists in (−∞, ∞), and by (37) we finally obtain

Q_{1}^{∗}(s) = \frac{2^{s}}{1 − 2^{s}} φ^{∗}(s) = (\frac{2^{s}}{1 − 2^{s}} )^{2} Φ^{∗}(s) ,  (41)

where Φ^{∗}(s) is an entire function, as we noticed above, and defined as

Φ^{∗}(s) = ∫_{0}^{∞} e_{d}(x)e^{−x}φ(2x)x^{s−1}dx .  (42)

To assess asymptotically L̃(z) we use property (P4) of the Mellin transform, that is, we must estimate residues of (41). Note that χ_{k} = 2πik/L (k = 0, ±1, ±2, …) are the solutions of 1 − 2^{−s} = 0, and the main contribution to the asymptotics comes from χ_{0} = 0. Since φ(z) ∼ 1 + O(z^{−M}) for any M > 0, we obtain, for z → ∞

L̃(z) = log_{2} z · \frac{Φ^{∗}(0)}{L} − \frac{Φ^{∗}(0)}{L} − \frac{Φ^{∗′}(0)}{L^{2}} + P_{1}(log_{2} z) + O(z^{−M})  (43)

where

Φ^{∗}(0) = ∫_{0}^{∞} e^{−x}e_{d}(x)φ(2x)\frac{dx}{x} ,  (44)
Φ^{∗′}(0) = ∫_{0}^{∞} e^{−x}e_{d}(x)φ(2x)\frac{log x}{x}dx ,  (45)

and

P_{1}(x) = ∑_{k≠0} (\frac{Φ^{∗}(−χ_{k})(x − 1)}{L} − \frac{Φ^{∗′}(−χ_{k})}{L^{2}}) e^{2kπix} .  (46)

To complete the proof, we need to evaluate Φ^{∗}(0),Φ^{∗}(χ_{k}), Φ^{∗′}(0), and Φ^{∗′}(χ_{k}). The first one is rather easy, and by properties of the Mellin transform we find out that Φ^{∗}(0) = L = log 2. The derivatives are harder to estimate. In order to treat the other constants, we first present another evaluation of Φ^{∗}(0).

Define the function 1(x) as follows

1(x) = \{ \begin{matrix}1 & if x ≥ 1 \\ 0 & if x < 1 . \end{matrix}  (47)

Then, we can write

Φ^{∗}(0) = ∫_{0}^{∞}(φ(2x) − φ(x))\frac{dx}{x}
= ∫_{0}^{∞}(φ(2x) − 1(2x))\frac{dx}{x} + ∫_{0}^{∞}(1(2x) − 1(x))\frac{dx}{x} + ∫_{0}^{∞}(1(x) − φ(x))\frac{dx}{x} (48)
= ∫_{0}^{∞}(1(2x) − 1(x))\frac{dx}{x} = ∫_{1/2}^{1} \frac{dx}{x} = L ,

since after the substitution u = 2x in the second line of the above display, the first and the second integrals cancel. A similar derivation shows that Φ^{∗}(χ_{k}) = 0, thus one proves formula (10) for P_{1}(x) in Theorem 1(i).

To complete the proof of Theorem 1(i), we must establish depoissonization of L̃(z) which follows directly from Theorem 3. For the reader convenience, we repeat the arguments below: Observe that due to (43) condition (28) of the Depoissonization Lemma (growth estimate inside a cone S_{θ}) is automatically satisfied. Assume now z ∉ S_{θ} for some 0 < θ < π/2, and define 1 > α > cos θ := c. We apply induction and the increasing domains as in the proof of Theorem 3. Consider (32), and by induction hypothesis assume that |L̃(z/2)e^{z/2}| ≤ βe^{α|z/2|} for some β > 1. Then:

|L̃(z)e^{z}| ≤ |f_{d}(z/2)e^{z/2}||L̃(z/2)e^{z/2}| + |f_{d}(z/2)e^{z}|
≤ |f_{d}(z/2)|e^{c|z/2|}βe^{α|z/2|} + |f_{d}(z/2)|e^{c|z|}
≤ βe^{α|z/2|}

where the last inequality holds as long as we choose ξ such that for |z| > ξ we have |f_{d}(z/2)| ≤ e^{(α−c)|z/2|}. This completes the depoissonization argument, and after an application of the depoissonization Lemma Theorem 1(i) is finally proved.

The proof of Theorem 1(ii) concerning the variance of R_{N,d} is more intricate but it proceeds along the same lines as the above derivation. In particular, we deal now with the functional equation (5) which we repeat below

W̃(z) = f_{d}(z/2)W̃(z/2) + 2L̃(z)f_{d}(z/2) .  (49)

Using the idea from Section 3.1 (cf. (21) – (25)) we find immediately an explicit solution to (49), namely

W̃(z)φ(z) = 2 ∑_{n=0}^{∞} L̃(z2^{−n−1})φ(z2^{−n−1}) ,

which becomes, after inserting solution (33) for L̃(z)

W̃(z)φ(z) = 2 ∑_{j=1}^{∞} ∑_{n=1}^{∞} φ(z2^{−n−j}) .  (50)

To establish the asymptotics of W̃(z), thus also for F′′_{N}(1), we proceed as before through the Mellin transform and depoissonization. Let Q_{2}(z) = W̃(z)φ(z) and Q_{2}^{∗}(s) be its Mellin transform. Then, from the harmonic sum formula

Q_{2}^{∗}(s) = 2 (\frac{2^{s}}{1 − 2^{s}} )^{2} φ^{∗}(s) = 2 (\frac{2^{s}}{1 − 2^{s}} )^{3} Φ^{∗}(s)

where Φ^{∗}(s) is defined as before. Noting that

\frac{1}{(1 − 2^{−s})^{3}} = \frac{1}{s^{3}L^{3}} + \frac{3}{2s^{2}L^{2}} + \frac{1}{sL} + O(1) ,

and computing residues, we easily obtain the following for z → ∞ and any M > 0

Q_{2}(z) = log_{2}^{2} z − (3 + 2C_{d}) log_{2} z + 2 + 3C_{d} + D_{d} + P_{2}(log_{2} z) + O(z^{−M})

with C_{d} from (8), P_{2}(x) is a periodic function with period 1, and

D_{d} = \frac{1}{L^{3}} Φ^{∗′′}(0) = \frac{1}{L^{3}} ∫_{0}^{∞} e_{d}(x)e^{−x}φ(2x)\frac{log^{2} x}{x}dx .

The rest is easy. We first depoissonize Q_{2}(z) using the same arguments as above for L̃(z). Since φ(z) = 1+O(z^{−M}) for any M > 0, we easily find the second factorial moment of R_{N,d}. Using Depoissonization Lemma, the proof of Theorem 1, apart from the asymptotics of the constant term for d → ∞, is completed.

Case B: Large d

Now, we analyze the case when d → ∞. Considering more carefully Φ^{∗′}(0), we obtain the following (with 1(x) defined in (49))

Φ^{∗′}(0) = −L^{2}/2 + L ∫_{0}^{∞}(1(x) − φ(x))\frac{dx}{x}
= −L^{2}/2 + L ∫_{0}^{1} φ(x)\frac{dx}{x} + L ∫_{1}^{∞}(1 − φ(x))\frac{dx}{x} .  (51)

We now estimate the above two integrals of (51) for large d. The first one is easy. Note that φ(x) ≤ 1 − e_{d}(x)e^{−x} = e^{−x} ∑_{k=1}^{∞} x^{d+k}/(d + k)!. Thus,

0 ≤ ∫_{0}^{1} φ(x)\frac{dx}{x} ≤ ∑_{k=1}^{∞} \frac{1}{(d + k)!(d + k)} = O (\frac{1}{d!}) .

In order to estimate the second integral of (51) – which we call I(d) – we recognize its similarity to the incomplete gamma function (cf. [1]). We shall derive separately a lower and an upper bound on I(d). For the lower bound, we observe that

I(d) = ∫_{1}^{∞}(1 − φ(x))\frac{dx}{x} ≥ ∫_{1}^{∞} e_{d}(x)e^{−x} \frac{dx}{x} = ∫_{1}^{∞} \frac{e^{−t}}{t} dt + ∑_{i=1}^{d} \frac{1}{i!}(Γ(i) − γ(i, 1))
= E_{1}(1) + ∑_{i=1}^{d} \frac{1}{i!}(Γ(i) − γ(i, 1))  (52)

where E_{1}(1) = −γ − ∑_{n≥1}(−1)^{n}/(nn!) is a particular value of the exponential integral E_{1}(z) = ∫_{z}^{∞} t^{−1}e^{−t} dt, and γ(a, z) = ∫_{0}^{z} e^{−t}t^{a−1}dt the incomplete gamma function (cf. [1]).

Now we must estimate the second term of (52). Note that

\frac{γ(i + 1, 1)}{i!} = \frac{1}{e} ∑_{k=i+1}^{∞} \frac{1}{k!} ,

and (cf. Knuth [27] Ex. 1.2.7–1.2.13),

\frac{1}{e} ∑_{k=1}^{∞} \frac{H_{k}}{k!} = e ∑_{n=1}^{∞} \frac{(−1)^{n+1}}{nn!} ,

where H_{k} is the kth harmonic number. Hence

I(d) ≥ H_{d} − γ + \frac{1}{e} ∑_{k≥1} \frac{H_{k}}{k!} − \frac{1}{e} ∑_{i=1}^{d} \frac{1}{i} ∑_{k≥i} \frac{1}{k!}
= H_{d} − γ + \frac{1}{e} ∑_{k>d}(H_{k} − H_{d})
≥ H_{d} − γ = log d + O(d^{−1}),  d → ∞.

An upper bound on I(d) is more intricate. First, we rewrite I(d) as

I(d) = ∫_{1}^{d+1}(1 − φ(x))\frac{dx}{x} + ∫_{d+1}^{∞}(1 − φ(x))\frac{dx}{x} .  (53)

The first integral above denoted as I_{1}(d) fulfills

I_{1}(d) ≤ ∫_{1}^{∞} \frac{dx}{x} = log(d + 1) = log d + O(d^{−1}), d → ∞.

In the following we estimate the second integral of (53) which we denote as I_{2}(d). We start from the observation that for 0 ≤ a_{k} ≤ 1

1 − ∑_{k≥0} a_{k} ≤ ∏_{k≥0}(1 − a_{k}),

which can easily be proved taking logarithms. Therefore

I_{2}(d) ≤ ∑_{k≥0} I_{2,k}(d),

with

I_{2,k}(d) = ∫_{d+1}^{∞} e_{d}(x2^{k})e^{−x2^{k}} \frac{dx}{x} = ∫_{(d+1)2^{k}}^{∞} e_{d}(u)e^{−u} \frac{du}{u} .

Now, for u ≥ d,

e_{d}(u) ≤ (d + 1)\frac{u^{d}}{d!} ,

so that

I_{2,k}(d) ≤ \frac{d + 1}{d!} ∫_{(d+1)2^{k}}^{∞} u^{d−1}e^{−u} \frac{du}{u} .

Since the function f(u) = u^{d}e^{−u} takes its maximum for u = d and is monotonically decreasing for u ≥ d, we have for u ≥ (d + 1)2^{k}

u^{d−1}e^{−u} ≤ u^{−2} ((d + 1)2^{k})^{d+1} e^{−(d+1)2^{k}} .

Hence

I_{2,k}(d) ≤ \frac{d + 1}{d!} ((d + 1)2^{k})^{d+1} e^{−(d+1)2^{k}} ∫_{(d+1)2^{k}}^{∞} \frac{du}{u^{2}}
= \frac{d + 1}{d!} (\frac{d + 1}{e})^{d} 2^{k(d+1)}e^{−(d+1)(2^{k}−1)}
≤ \frac{d + 1}{d!} (\frac{d + 1}{e})^{d} (\frac{2}{e})^{k(d+1)} ,

so that ∑_{k≥0} I_{2,k}(d) is exponentially small in d → ∞. Altogether we have proved that

I(d) = log d + O(d^{−1}),  d → ∞ .  (54)

Putting everything together, we finally obtained the asymptotics for C_{d} described in Theorem 1(i), that is,

C_{d} = \frac{Φ^{∗′}(0)}{L^{2}} = log_{2} d − 1/2 + O(d^{−1}) .  (55)

Now, can turn our attention to establishing Theorem 1(ii) for large d. We recall that to compute the variance of R_{N,d} we need a precise estimate of (L̃(N))^{2}. But the asymptotic formula on (L̃(N))^{2} involves P_{1}^{2}(log_{2} N), where the zeroth term in the Fourier expansion of P_{1}^{2}(x) is not equal to zero (in contrast to P_{1}(x)). To recover this zeroth Fourier coefficient [P_{1}^{2}]_{0}, we need a new representation of P_{1}(x) which we recall below

P_{1}(x) = −\frac{1}{L^{2}} ∑_{k≠0} Φ^{∗′}(−χ_{k})e^{2kπix} .

Proceeding as with Φ^{∗′}(0) above we find

\frac{1}{L^{2}} Φ^{∗′}(χ_{k}) = −\frac{1}{Lχ_{k}} − \frac{1}{L} ∫_{0}^{d+1} φ(t)\frac{dt}{t^{1+χ_{k}}} + \frac{1}{L} ∫_{0}^{d+1}(1 − φ(t))\frac{dt}{t^{1+χ_{k}}} .

Using integration by parts the whole expression turns out to be

\frac{1}{L^{2}} Φ^{∗′}(χ_{k}) = −\frac{1}{Lχ_{k}} ∫_{0}^{∞} φ′(t)\frac{dt}{t^{χ_{k}}} .  (56)

For later use we note that in the last integral φ′(t) can be replaced by f′_{d}(t) = \frac{t^{d}}{d!} e^{−t} with an exponentially small error in d: We have

∫_{0}^{∞} φ′(t)t^{−χ_{k}}dt = ∫_{0}^{∞} f′_{d}(t)t^{−χ_{k}}dt + R(k, d),

with

R(k, d) = ∫_{0}^{∞} f′_{d}(t)(φ(2t) − 1)t^{−χ_{k}}dt + 2 ∫_{0}^{∞} φ_{0}(t)φ′(2t)t^{−χ_{k}}dt ,

so that, integrating by parts,

|R(k, d)| ≤ 2 ∫_{0}^{∞} f′_{d}(t)(φ(2t) − 1)dt .

Splitting up the range of integration according to the regions where either one of the factors of the integrand is small, |R(k, d)| can be estimated similarly to Φ^{∗′}(0), and R(k, d) turns out to be exponentially small in d (uniformly in k). Now,

∫_{0}^{∞} f′_{d}(t)t^{−χ_{k}}dt = \frac{1}{d!} ∫_{0}^{∞} t^{d−χ_{k}}e^{−t}dt = \frac{Γ(d + 1 − χ_{k})}{d!} ,

so that

−\frac{1}{L^{2}} Φ^{∗′}(−χ_{k}) = \frac{1}{Lχ_{k}} (\frac{Γ(d + 1 + χ_{k})}{d!} + R(k, d))

and we finally obtain

P_{1}(x) = \frac{1}{L} ∑_{k≠0} \frac{1}{χ_{k}} (\frac{Γ(d + 1 + χ_{k})}{d!} + R(k, d)) e^{2kπix} .  (57)

Finally, we can return to the evaluation of the zeroth Fourier coefficient [P_{1}^{2}]_{0} of P_{1}^{2}(x). According to the uniform exponential smallness of R(k, d) for d getting large we have

[P_{1}^{2}]_{0} = \frac{1}{L^{2}d!^{2}} ∑_{k≠0} \frac{1}{|χ_{k}|^{2}} |Γ(d + 1 + χ_{k})|^{2} + exp. small terms in d .  (58)

In order to evaluate the series occurring in (58) we observe that

\frac{1}{χ_{k}} \frac{Γ(d + 1 + χ_{k})}{d!} = \binom{d + χ_{k}}{d}Γ(χ_{k}) = ∑_{λ=0}^{d} \binom{λ + χ_{k} − 1}{λ}Γ(χ_{k}) .  (59)

Therefore the main term of P_{1}(x) coincides perfectly with the periodic function P_{1}(x) analyzed in [25] (cf. Theorem 3 of [25]) if d is replaced by d + 1. The mean [P_{1}^{2}]_{0} of the square of the latter fluctuating function can be analyzed by different methods. In [25] this quantity is treated using Hankel contour integrals with the “kernel”

\frac{Γ(j + z)Γ(j − z)}{e^{Lz} − 1} .

An earlier approach, using transformation results due to Ramanujan, may be found in [23]. From [25] we have the asymptotic result

[P_{1}^{2}]_{0} = \frac{1}{12} − \frac{1}{\sqrt{πd}} + O(\frac{1}{d}), d → ∞ .

This completes the proof of part (ii) of Theorem 1.

3.3 Asymptotic Limiting Distributions

In this subsection, we prove Theorem 2 concerning the asymptotic distribution of R_{N,d}. The Poisson generating function G̃(z, u) of R_{N,d} satisfies the functional equation (17). But, the Mellin transform of G̃(z, u) with respect to z does not exist since G̃(z, u) = O(1) for both z → ∞ and z → 0. Therefore, we introduce the new function H̃(z, u) = G̃(z, u) − 1 whose Mellin transform exists in ℜ(s) ∈ (−1, 0). Observe that H̃(z, u) fulfills the following functional equation

H̃(z, u) = uf_{d}(z/2)H̃(z/2, u) + (u − 1)f_{d}(z/2) .  (60)

This equation falls under our general equation (21) from Section 3.1. In particular, using the same arguments as before, we can write the solution for H̃(z, x) as

φ(z)H̃(z, u) = (u − 1) ∑_{n=0}^{∞} u^{n}φ(z2^{−n−1}) ,  (61)

where φ(z) is defined in (34).

We can now easily obtain asymptotic expansion for the probability generating function F_{N}(z). To derive it we apply the Mellin transform technique and then depoissonization since (60) satisfies all conditions of Theorem 3. As a result, we can prove the following

F_{N}(u) = u^{log_{2} N} \frac{Φ^{∗}(− log_{2} u)}{L(u − 1)} + P_{3}(log_{2} N) + O(N^{−1/2+ε}) .  (62)

where Φ^{∗}(s) is the Mellin transform of φ(2x) − φ(x) as discussed in Section 3.2 and P_{3} is a fluctuating function.

However, to derive the asymptotic distribution of R_{N,d}, we proceed in a different manner. Note that from (61) we obtain

\frac{1 − G̃(z, u)}{1 − u} = ∑_{k=0}^{∞} u^{k} \frac{φ(z2^{−k−1})}{φ(z)} .

But, using the definition of G̃(z, u) and comparing coefficients at u^{k} we arrive at

h̃_{k}(z) = [u^{k}]\frac{H̃(z, u)}{1 − u} = e^{−z} ∑_{N=0}^{∞} Pr\{R_{N,d} > k\}\frac{z^{N}}{N!} = \frac{φ(z2^{−k−1})}{φ(z)} = ∏_{j=1}^{k+1} f_{d}(z/2^{j})  (63)

with f_{d}(z) = 1 − e_{d}(z)e^{−z}. In the following we depoissonize (63).

We start with an estimate on h̃_{k}(z) inside the cone, thus we assume z ∈ S_{θ}. Then, ℜz = |z| cos(arg z) ≥ |z| cos θ, and

|f_{d}(z)| = |1 − e_{d}(z)e^{−z}| ≤ 1 + e_{d}(|z|)|e^{−z}| = 1 + e_{d}(|z|)e^{−ℜz}
≤ 1 + e_{d}(|z|)e^{−|z| cos θ} .

If we assume θ ≤ \frac{π}{3} , the last estimate reads

|f_{d}(z)| ≤ 1 + e_{d}(|z|)e^{−|z|/2} for all z ∈ S_{θ} .

Let ε > 0 be arbitrarily small. Then there exists K = K(ε) ≥ 1 such that e_{d}(x)e^{−x/2} < \frac{ε}{2} for all x > K, resp.

|f_{d}(z)| ≤ 1 + \frac{ε}{2}  for z ∈ S_{θ}, |z| > K(ε) .  (64)

Let us now consider the region |z| ≤ K: Since 0 ≤ f_{d}(x) = 1 − e_{d}(x)e^{−x} ≤ f_{d}(K) < 1 for x ∈ [0, K] and the function f_{d}(z) is continuous, there exists an open neighbourhood U(x) of each point x ∈ [0, K] with |f_{d}(z)| ≤ 1 for z ∈ U(x). Due to the compactness of [0, K], it follows that we can choose θ = θ(ε) small enough such that

|f_{d}(z)| ≤ 1  for z ∈ S_{θ(ε)}, |z| ≤ K(ε) .  (65)

By (65) it follows that

|f_{d}(z/2^{j})| ≤ 1  for j > log_{2}(|z|/K(ε)), z ∈ S_{θ(ε)}, |z| ≤ K(ε) .  (66)

To estimate h̃_{k}(z) we split the product from (63) into the regions 1 ≤ j ≤ min(k + 1, log_{2}(|z|/K(ε)) and j > min(k + 1, log_{2}(|z|/K(ε)), and we obtain

|h̃_{k}(z)| ≤ (1 + \frac{ε}{2})^{log_{2}(|z|/K(ε))} = (\frac{|z|}{K(ε)})^{log_{2}(1+ \frac{ε}{2} )}

i.e.,

|h̃_{k}(z)| ≤ |z|^{ε} for z ∈ S_{θ(ε)}, |z| ≥ 1 .

To complete the depoissonization, we need to estimate h̃_{k}(z) outside the cone. Thus, assume z ∉ S_{θ}, so that ℜz = |z| cos(arg z) ≤ |z| cos θ. Therefore, with C_{1} = C_{1}(d)

|e^{z}f_{d}(z)| = |e^{z} − e_{d}(z)| ≤ e^{ℜz} + e_{d}(|z|) ≤ (1 + C_{1})e^{|z| cos θ} .  (67)

Again we need a sharper estimate for |e^{z/2^{j}}f_{d}(z/2^{j})| if j is large: Since e^{z}f_{d}(z) is continuous with e^{0}f_{d}(0) = 0, there exists a constant C_{3} = C_{3}(d) > 0 such that

|e^{z}f_{d}(z)| < 1  for |z| < C_{3}.  (68)

From (68) it follows that

|e^{z/2^{j}}f_{d}(z/2^{j})| < 1  for j > log_{2}(|z|/C_{3}) .  (69)

Now we proceed as follows: Since ∑_{j=1}^{k+1} \frac{1}{2^{j}} = 1 − \frac{1}{2^{k+1}} , we have

|e^{z} ∏_{j=1}^{k+1} f_{d}(z/2^{j})| = |e^{z/2^{k+1}} ∏_{j=1}^{k+1}(e^{z/2^{j}} − e_{d}(z/2^{j}))| ,  (70)

which is estimated by

e^{ℜz/2^{k+1}} ∏_{1≤j≤K}(1 + C_{1}) e^{(|z| cos θ)/2^{j}} ,

where K = min(k + 1, log_{2}(|z|/C_{3})) . Therefore the estimate is less than or equal to

e^{(|z| cos θ)/2^{k+1}}(1 + C_{1})^{K}e^{|z| cos θ (1− \frac{1}{2^{k}})} ≤ (1 + C_{1})^{log_{2}(|z|/C_{3})}e^{|z| cos θ (1+ \frac{1}{2^{k+1}} − \frac{1}{2^{k}})}

≤ (\frac{|z|}{C_{3}} )^{log_{2}(1+C_{1})} e^{|z| cos θ} < e^{α|z|}

for cos θ < α < 1 and z → ∞, z ∉ S_{θ}.

Thus, we have established the conditions of the Depoissonization Lemma, and consequently get that

Pr\{R_{N,d} > k\} = \frac{φ(N2^{−k−1})}{φ(N)} + O (N^{−1/2+ε}) ,

for any ε > 0 where the O–constant is independent of k.

To complete the proof of Theorem 2, it suffices to set k = ⌊log_{2} N + t⌋ − 1, and notice that Pr{R_{N,d} > log_{2} N + t − 1} = Pr{R_{N,d} > ⌊log_{2} N + t⌋ − 1} .

ACKNOWLEDGEMENT

We thank Philippe Jacquet (INRIA, Rocquencourt) for several valuable discussions concerning the depoissonization issues.

References

[1] M. Abramowitz, and I. Stegun (Eds.), Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, John Wiley & Sons, New York 1972.

[2] D. Aldous, Probability Approximations via the Poisson Clumping Heuristic, Springer Verlag, New York 1989.

[3] Anderson, C. (1970). Extreme value theory for a class of discrete distributions with applications to some stochastic processes. J. Appl. Probability, vol. 7, pp. 99–113.

[4] G. Brassard and P. Bratley, Algorithmics. Theory and Practice, Prentice Hall, Englewood Cliffs, 1988.

[5] B. Davies, Integral Transforms and Their Applications, Springer-Verlag, 1978.

[6] G. Doetsch, Handbuch der Laplace Transformation, Vol. 1–3, Birkhäuser Verlag, Basel, 1955.

[7] J. Fill, H. Mahmoud and W. Szpankowski, On the Distribution for the Duration of a Randomized Leader Election Algorithm, Annals of Applied Probability, to appear.

[8] P. Flajolet, P. Jacquet, and W. Szpankowski, Mellin Transform of a Complex Variable: A Simple Extension, Unpublished Manuscript, 1995.

[9] P. Flajolet and G.N. Martin, Probabilistic Counting Algorithms for Data Base Applications, J. Comp. Syst. Sci., 31, 182-209, 1985.

[10] P. Flajolet and N. Saheb, The Complexity of Generating an Exponentially Distributed Variate, J. Algorithms, 7, 463-488, 1986.

[11] P. Flajolet and B. Richmond, Generalized Digital Trees and Their Difference-Differential Equations, Random Structures and Algorithms, 3, 305-320, 1992.

[12] P. Flajolet, X. Gourdon and P. Dumas, Mellin Transforms and Asymptotics: Harmonic Sums, Theoretical Computer Science, 144, 3-58, 1995.

[13] P. Grabner, Searching for losers, Random Structures and Algorithms, 4, 99–110, 1993.

[14] A. Greenberg, P. Flajolet and R. Ladner, Estimating the Multiplicities of Conflicts to Speed Their Resolution in Multiple Access Channels, JACM, 34, 289-325, 1987.

[15] P. Henrici, Applied and Computational Complex Analysis, vol. 2, John Wiley & Sons 1977.

[16] P. Jacquet, M. Régnier, Trie partitioning process: limiting distributions. Lecture Notes in Computer Science, vol. 214, pp. 196–210, Springer, New York 1986.

[17] P. Jacquet and W. Szpankowski, Ultimate Characterizations of the Burst Response of an Interval Searching Algorithm: A Study of a Functional Equation, SIAM J. Computing, 18, 777-791, 1989.

[18] P. Jacquet and W. Szpankowski, Asymptotic behavior of the Lempel-Ziv Parsing Scheme and Digital Search Trees, Theoretical Computer Science, 144, 161-197, 1995.

[19] P. Jacquet and W. Szpankowski, Analytical depoissonization and its applications, preprint.

[20] M. Kendall and A. Stuart, The Advanced Theory of Statistics, 4th Ed. vol. 1, Charles and Griffin and Co. Ltd., London (1977).

[21] J.F.C. Kingman, Subadditive processes, in Ecole d’Eté de Probabilités de Saint-Flour V-1975, Lecture Notes in Mathematics, 539, Springer-Verlag, Berlin (1976).

[22] P. Kirschenhofer and H. Prodinger, On the Analysis of Probabilistic Counting, Lecture Notes in Mathematics, 1452 (Eds. E. Hlawka and R. Tichy), 117-120, Springer Verlag 1990.

[23] P. Kirschenhofer and H. Prodinger, On some applications of formulæ of Ramanujan in the analysis of algorithms, Mathematika, 38, 14-33, 1991.

[24] P. Kirschenhofer and H. Prodinger, Approximate Counting: An Alternative Approach, Informatique Théorique et Applications/Theoretical Informatics and Applications, 25, 43-48, 1991.

[25] P. Kirschenhofer and H. Prodinger, A Result in Order Statistics Related to Probabilistic Counting, Computing, 51, 15-27, 1993.

[26] P. Kirschenhofer, H. Prodinger, and W. Szpankowski, How to count quickly and accurately: A unified analysis of probabilistic counting and other related problems, International Conference on Automata, Languages, and Programming (ICALP’92), Vienna, LNCS, No. 623, 211-222, 1992.

[27] D. Knuth, The Art of Computer Programming: Fundamental Algorithms, vo. 1., Addison-Wesley, Reading 1973.

[28] D.E. Knuth, The Art of Computer Programming. Sorting and Searching, vol. 3., Addison-Wesley, Reading, MA 1973.

[29] R. Morris, Counting Large Numbers of Events in Small Registers, Comm. ACM, 21, 840-842, 1978.

[30] B. Pittel and H. Rubin, How Many Random Questions Are Necessary to Identify n Distinct Objects?, J. Combinatorial Theory, Ser. A, 55, 292-312, 1990.

[31] H. Prodinger, How to Select a Loser, Discrete Math., 120, 149–159, 1993.

[32] B. Rais, P. Jacquet, and W. Szpankowski, Limiting Distribution for the Depth in Patricia Tries, SIAM J. Discrete Mathematics, 6, 197-213, 1993.

[33] M. Régnier and P. Jacquet, Normal Limiting Distribution of the Size of Tries, Proc. Performance’87, Amsterdam: North-Holland, 209-223, 1987.

[34] U. Schmid, The Average CRI-length of a Tree Collision Resolution Algorithm in Presence of Multiplicity-Dependent Capture Effects, Proc. ICALP 92, Vienna, 223-234, 1992.

[35] W. Szpankowski, Solution of a Linear Recurrence Equation Arising in the Analysis of Some Algorithms, SIAM J. Alg. Disc. Methods, 8, 233-250, 1987.

[36] W. Szpankowski, Patricia Tries Again Revisited, JACM, 37, 691-711, 1990.

[37] W. Szpankowski and V. Rego, Yet Another Application of a Binomial Recurrence. Order Statistics, Computing, 43, 401-410, 1990.