Klaus Dohmen
email: dohmen@htwm.de

André Pönitz
email: poenitz@htwm.de

Peter Tittmann

A new two-variable generalization of the chromatic polynomial

Keywords: chromatic polynomial, set partition, broken circuit, pathwidth, chromatic symmetric function

Let P(G; x, y) be the number of vertex colorings φ : V → {1, ..., x} of an undirected graph G = (V, E) such that for all edges {u, v} ∈ E the relations φ(u) ≤ y and φ(v) ≤ y imply φ(u) = φ(v). We show that P(G; x, y) is a polynomial in x and y which is closely related to Stanley's chromatic symmetric function, and which simultaneously generalizes the chromatic polynomial, the independence polynomial, and the matching polynomial of G. We establish two general expressions for this new polynomial: one in terms of the broken circuit complex, and one in terms of the lattice of forbidden colorings. Finally, we give explicit expressions for the generalized chromatic polynomial of complete graphs, complete bipartite graphs, paths, and cycles, and show that P(G; x, y) can be evaluated in polynomial time for trees and graphs of restricted pathwidth.

Introduction

All graphs in this paper are assumed to be finite, undirected and without loops or multiple edges. We write G = (V, E) to denote that G is a graph having vertex set V and edge set E.

The well-known chromatic polynomial P(G; y) of a graph G = (V, E) gives the number of vertexcolorings of G with at most y colors such that adjacent vertices receive different colors. For an introduction to chromatic polynomials the reader is referred to READ [START_REF] Read | An introduction to chromatic polynomials[END_REF].

TUTTE [START_REF] Tutte | Graph Theory[END_REF][START_REF] Tutte | A ring in graph theory[END_REF] has generalized the chromatic polynomial to the Tutte polynomial T (G; x, y). We propose a different generalization of the chromatic polynomial by weakening the requirements for proper colorings: Let X = Y ∪ Z, Y ∩ Z = / 0, be the set of available colors with |X | = x and |Y | = y. Then, a generalized proper coloring of G is a map φ : V → X such that for all {u, v} ∈ E with φ(u) ∈ Y and φ(v) ∈ Y the relation φ(u) = φ(v) holds. Consequently, adjacent vertices can be colored alike only if the color is chosen from Z. In order to distinguish between these two sets of colors we call the colors of Y proper and the colors of Z improper. The number of generalized proper colorings of G is denoted by P(G; x, y), and as we will see in Theorem 1 below, this number turns out to be a polynomial in x and y, which we refer to as the generalized chromatic polynomial of G. Despite its simple definition, this new two-variable polynomial shows some interesting properties which are discussed later on in more detail:

• It generalizes the chromatic polynomial, the independence polynomial, and the matching polynomial.

• It is closely related to Stanley's chromatic symmetric function [START_REF] Stanley | Enumerative Combinatorics[END_REF].

• It distinguishes all non-isomorphic trees with up to nine vertices.

• It satisfies both an edge decomposition formula and a vertex decomposition formula.

In subsequent sections, we give explicit expressions for the generalized chromatic polynomial of complete graphs, complete bipartite graphs, paths, and cycles, and show that our generalized chromatic polynomial can be evaluated in polynomial time for trees and graphs of restricted pathwidth.

Basic properties

Let G = (V, E) be a graph, and let n = |V | and m = |E| denote the number of vertices and edges of G, respectively. By considering Y = / 0 we get P(G; x, 0) = x n . The case Z = / 0 shows that P(G; y, y) coincides with the usual chromatic polynomial P(G; y). For any edge e ∈ E let Ge and G/e be the graphs obtained from G by deleting resp. contracting e and then, in the resulting multigraph, replacing each class of parallel edges by a single edge. By the edge decomposition formula for the usual chromatic polynomial we have P(G; y, y) = P(Ge; y, y) -P(G/e; y, y).

(

Since vertices belonging to different components of G can be colored independently, we find that if G 1 , ..., G k are the connected components of G, then

P (G; x, y) = k ∏ i=1 P (G i ; x, y) . (2)
The following theorem shows that P (G; x, y) can be expressed in terms of chromatic polynomials of subgraphs of G. Let G -X be the subgraph obtained from G by removing all vertices of a vertex subset X ⊆ V . For short, we write Gv instead of G -{v}.

Theorem 1 Let G be a graph. Then, P(G; x, y) = ∑ X⊆V (xy) |X| P (G -X; y) .

Proof. Every generalized proper coloring of G can be obtained by first choosing a subset X of V that is colored with colors of Z. There are (xy) |X| different colorings for these vertices. The remaining subgraph has to be colored properly using colors of Y , for which there are P(G -X; y) possibilities. ✷ Corollary 2 Let a i denote the number of independent vertex sets of cardinality ni where n denotes the number of vertices of G. Then the independence polynomial of G, defined by I (G, x) := ∑ n i=0 a i x i , satisfies

I (G, x) = P (G; x + 1, 1) .
In particular, P(G; 2, 1) is the number of independent vertex sets of G.

Proof. By Theorem 1 we have

P(G; x + 1, 1) = ∑ X⊆V x |X| P (G -X; 1) .
Note that the chromatic polynomial P (G -X; 1) equals 1 if V -X is an independent set (a set of isolated vertices) of G; otherwise, P (G -X; 1) = 0. Consequently, the coefficient of x i in P(G; x + 1, 1) counts the independent sets of cardinality ni in G. ✷

Theorem 3

The generalized chromatic polynomial P(G; x, y) satisfies

∂ ∂x P(G; x, y) = ∑ v∈V P(G -v; x, y).
Proof. From Theorem 1 it follows that

∑ v∈V P(G -v; x, y) = ∑ v∈V ∑ W ⊆V \{v} (x -y) |W | P(G -W -v; y) = ∑ X⊆V |X|(x -y) |X|-1 P(G -X; y) = ∂ ∂x ∑ X⊆V (x -y) |X| P(G -X; y) = ∂ ∂x P(G; x, y),
which gives the result. ✷

We close this section by remarking that our generalized chromatic polynomial is not an evaluation of the Tutte polynomial [START_REF] Tutte | Graph Theory[END_REF][START_REF] Tutte | A ring in graph theory[END_REF]: The Tutte polynomial T (G; x, y) equals x n-1 for every tree having n vertices regardless of the structure of the tree. However, there are non-isomorphic trees with different generalized chromatic polynomials. The smallest example for such a pair of trees is presented in Figure 1. The generalized chromatic polynomials of these two trees P 4 and S 3 are respectively given by

P(P 4 ; x, y) = x 4 -3x 2 y + y 2 + 2xy -y, (3)
P(S 3 ; x, y) = x 4 -3x 2 y + 3xy -y. ( 4
)
Fig. 1: Non-isomorphic trees P 4 and S 3 having four vertices and distinct generalized chromatic polynomials.

The partition representation of P(G; x, y)

Let G = (V, E) be a graph and Π(V ) be the partition lattice of the vertex set V . A subset W ⊆ V is connected if the subgraph of G induced by the vertices of W is connected. Let Π G be the set of all partitions π of V such that every block of π is connected. Then Π G forms a geometric sublattice of Π(V ) (see e.g. [START_REF] Stanley | Enumerative Combinatorics[END_REF]). The unique minimal element 0 in Π G is the finest partition consisting only of singletons.

We show that Π G corresponds to the lattice of forbidden colorings of G. Let f (π) be the number of colorings of V with y proper colors and xy improper colors such that 1. the vertices of each block B of π all get the same color, 2. if |B| = 1 then that color can be any color of X while if |B| ≥ 2 then that color must be from Y , and 3. no two blocks which are connected by an edge may be assigned the same color of Y .

Consequently, no proper coloring of G contributes to f (π) for π > 0. We denote by |π| the number of blocks of π. Let k 1 (π) be the number of singletons of π. We obtain

x k 1 (π) y |π|-k 1 (π) = ∑ σ∈Π G σ≥π f (σ).
Our generalized chromatic polynomial P(G; x, y) corresponds to f ( 0). Thus, by Möbius inversion we obtain the following theorem.

Theorem 4

The generalized chromatic polynomial can be expressed as a sum over the lattice Π G :

P(G; x, y) = ∑ π∈Π G µ( 0, π)x k 1 (π) y |π|-k 1 (π) . (5)

Fig. 2: A graph with four vertices

Consider the graph G shown in Figure 2. The lattice Π G of forbidden colorings of G and the values of the Möbius function µ( 0, π) for this lattice is presented in Figure 3. The resulting polynomial is

P(G; x, y) = x 4 -4x 2 y + 4xy + y 2 -2y.
In Eq. ( 5) the sum is extended over all blocks π ∈ Π G . In the worst case this gives B(n) terms, where B(n) denotes the nth Bell number. A combinatorial interpretation of the coefficients is given in the next section. STANLEY [START_REF] Stanley | Enumerative Combinatorics[END_REF] defines the chromatic symmetric function for any graph G = (V, E) by

X G := ∑ κ ∏ i x |κ -1 (i)| i
where the sum is over all proper colorings κ : V → {1, . . . , n} of G. Thus, the coefficient of

x λ 1 1 • • • x λ r r in X G is the number of proper colorings of G such that λ k vertices are colored k for k = 1, . . . , r. Stanley shows that X G = ∑ π∈Π G µ 0, π p λ(π)
where λ (π) denotes the type of π, that is, λ(π) = (λ 1 , . . . , λ r ) where the λ k 's denote the sizes of the blocks of π, and where p denotes the power sum symmetric function, which for λ = (λ 1 , . . . , λ r ) is defined by (cf. [START_REF] Stanley | Enumerative Combinatorics[END_REF])

p λ := |λ| ∏ k=1 ∑ i x λ k i ,
where |λ| denotes the number of parts of λ. In view of this and Theorem 4, the coefficients of our generalized chromatic polynomial and of Stanley's chromatic symmetric function are related via the lattice Π G .

We proceed by establishing a connection between our polynomial and the matching polynomial. Recall that a matching of a graph G = (V, E) is a subset F of E such that no two edges of F share a common vertex. Let m k be the number of matchings of cardinality k in G. The matching polynomial of G is defined by

M (G; x) := ⌊n/2⌋ ∑ k=0 (-1) k m k x n-2k .

Corollary 5

The matching polynomial is related to the generalized chromatic polynomial via z n-2k M (G; z) = x n-2k y k P (G; x, y) , where z n-2k M (G; z) and x n-2k y k P (G; x, y) denote the coefficients of z n-2k in M (G; z) and x n-2k y k in P (G; x, y), respectively.

Proof. By Theorem 4 we have

x n-2k y k P (G; x, y) = ∑ π∈Π G k 1 (π)=n-2k |π|=n-k µ( 0, π), (6)
where the Möbius function of an interval [ 0, π] in the partition lattice is given by (c.f. ROTA [START_REF] Rota | On the Foundations of Combinatorial Theory: I. Theory of Möbius functions[END_REF])

µ( 0, π) = (-1) n-|π| ∏ A∈π (|A| -1)! . ( 7
)
Here A ∈ π indicates that A is a block of π. A partition of {1, ..., n} with exactly n -2k singletons and nk blocks must contain exactly k two-element blocks. If π is such a partition, then Eq. ( 7) yields

µ( 0, π) = (-1) k .
This is also the value of the Möbius function in Π G since the interval 0, π coincides with the corresponding interval of the partition lattice. Thus, Eq. ( 6) becomes

x n-2k y k P (G; x, y) = (-1) k |{π ∈ Π G : k 1 (π) = n -2k, |π| = n -k}| .
Since there is a one-to-one correspondence between the set of partitions on the right-hand side of the formula and the set of matchings of cardinality k, the statement of the corollary immediately follows. ✷ Thus the generalized chromatic polynomial P (G; x, y) incorporates the matching polynomial of G as well as the independence polynomial of G. The matching polynomial and the independence polynomial are also related by the line graph via the identity (cf. [START_REF] Godsil | Algebraic Graph Theory[END_REF])

M (G; x) = (-1) n x n-2m I L (G) ; -x 2 ,
but using P(G; x, y) eliminates the need for using L(G). Thus, our new polynomial may be considered as a generalization of the chromatic polynomial, the independence polynomial, and the matching polynomial.

We proceed by establishing a representation of our new polynomial in terms of falling factorials. A partition π of the vertex set V (G) is called independent if every block of π is an independent set of G. We use λ ⊢ n to denote that λ is a number partition of n. We obtain the following partition representation:

Theorem 6 For each λ ⊢ n = |V (G)|,
let a λ be the number of independent partitions of G of type λ. Then, the generalized chromatic polynomial satisfies

P(G; x, y) = ∑ λ⊢n a λ y |λ|-k 1 (λ) k 1 (λ) ∑ k=0 k 1 (λ) k (x -y) k (y -|λ| + k 1 (λ)) k 1 (λ)-k ,
where y i denotes the falling factorial y (y -1)

• • • (y -i + 1).
Proof. The formula counts all proper generalized colorings of G. To avoid multiple counting, we color all vertices of each block having more than one vertex by one proper color of Y . Different blocks containing two or more vertices are colored differently. This gives y |λ|-k 1 (λ) possibilities to color these blocks. The remaining singletons of the partition may be colored arbitrarily with proper colors or improper ones. When using proper colors the singletons must be colored differently. This gives the second sum on the right-hand side, and we are done. ✷

Corollary 7

Let a i j be the number of independent partitions of G with exactly i singletons and j blocks having two or more vertices. Then, the following equation holds:

P (G; x, y) = n ∑ i=0 n ∑ j=0 a i j i ∑ k=0 i k (x -y) k y i+ j-k .
Corollary 7 suggests the introduction of a simpler two-variable polynomial of the form

Q (G; x, y) = n ∑ i=0 n ∑ j=0 a i j x i y j , (8)
which comprises the same amount of information as our generalized chromatic polynomial. However, the simplicity of defining the polynomial in this way has to be paid for with loss of many nice properties, e.g., the multiplicity with respect to components.

The coefficients in terms of broken circuits

One of the most important results about the usual chromatic polynomial is Whitney's broken circuit theorem [START_REF] Whitney | A logical expansion in mathematics[END_REF], which expresses the coefficients of the chromatic polynomial in terms of broken circuits. Many results can be deduced directly from Whitney's broken circuit theorem; see e.g., LAZEBNIK [START_REF] Lazebnik | Some corollaries of a theorem of Whitney on the chromatic polynomial[END_REF].

Let G = (V, E) be a graph where |V | = n and E is linearly ordered. A broken circuit of G is obtained from the edge set of a cycle of G by removing its maximum edge. The broken circuit complex of G, abbreviated to BC (G), is the set of all non-empty subsets of the edge set not including any broken circuit as a subset. The definition of a broken circuit goes back to WHITNEY [START_REF] Whitney | A logical expansion in mathematics[END_REF], while the broken circuit complex was initiated by WILF [START_REF] Wilf | Which polynomials are chromatic?[END_REF] and further investigated by BRYLAWSKI and OXLEY [START_REF] Brylawski | The broken circuit complex[END_REF][START_REF] Brylawski | The broken-circuit-complex: its structure and factorizations[END_REF]. In fact, BC (G) is an abstract simplicial complex in the sense of topology, whence we refer to each I ∈ BC (G) as a face of BC (G).

Whitney's broken circuit theorem [START_REF] Whitney | A logical expansion in mathematics[END_REF] states that for any y ∈ N,

P(G; y) = n ∑ k=0 (-1) k b k y n-k ,
where b 0 = 1 and b k , k > 0, counts the faces of cardinality k in the broken circuit complex of G. Note that while the definition of the broken circuit complex depends on the ordering of the edges, the same does not apply to the b k 's.

We now generalize Whitney's broken circuit theorem to our new two-variable polynomial P(G; x, y). The proof of this generalization (and thus of Whitney's original result) is facilitated by applying the following inclusion-exclusion variant, which is an immediate consequence of [START_REF] Dohmen | Improved Bonferroni inequalities via union-closed set systems[END_REF]Corollary 3.5].

Proposition 8

Let {A e } e∈E be a finite family of finite sets, where E is endowed with a linear ordering relation. Furthermore, let F be a set of non-empty subsets of E such that for any F ∈ F , i∈F

A i ⊆ j>max F A j . Then, e∈E A e = ∑ I⊆E,I = / 0 I ⊇F(∀F∈F ) (-1) |I|-1 i∈I A i .
For any graph G = (V, E) and any subset I of E we use G[I] to denote the graph having vertex set {v ∈ V | v ∈ e for some e ∈ I} and edge set I.

We are now ready to state the main result of this section, which coincides with Whitney's broken circuit theorem in the diagonal case where x = y. Theorem 9 Let G = (V, E) be a graph with n vertices and m edges and whose edge set is endowed with a linear ordering relation. Furthermore, let x ∈ N and y ∈ {0, . . . , x}. Then,

P(G; x, y) = m ∑ k=0 k ∑ l=0 (-1) k b kl x n-k-l y l , ( 9
)
where b 00 = 1 and b kl , k > 0, counts all faces I of cardinality k in the broken circuit complex of G such that G[I] has exactly l connected components.

Proof. Define F as the set of broken circuits of G and for any edge e of G define

A e := { f : V → {1, . . . , x} | f (u) = f (v) ≤ y}, e = {u, v}.
Then, the requirements of Proposition 8 are satisfied, and thus we obtain

P(G; x, y) = x n - e∈E A e = x n + m ∑ k=1 (-1) k ∑ I∈BC (G) |I|=k i∈I A i . (10
A i = x n-n I y c I = x n-m I -c I y c I . (11)
Now, put ( 11) into [START_REF] Stanley | Enumerative Combinatorics[END_REF] and take into account that c I ≤ m I and m I = |I|. ✷

In view of the results in [START_REF] Dohmen | Improved Bonferroni inequalities via union-closed set systems[END_REF] it can even be shown that for any r ∈ N,

P (G; x, y) ≤ r ∑ k=0 k ∑ l=0 (-1) k b kl x n-k-l y l (r odd), P (G; x, y) ≥ r ∑ k=0 k ∑ l=0 (-1) k b kl x n-k-l y l (r even),
where b kl is defined as in Theorem 9. The proof of this 'Bonferroni-like' result is left to the reader. We further remark that the independence polynomial I(G; x) satisfies

I(G; x) = m ∑ k=0 k ∑ l=0 (-1) k b kl (x + 1) n-k-l ,
which immediately follows from I(G; x) = P(G; x + 1, 1) and Theorem 9, and that (9) can be restated as

P(G; x, y) = x n Z G; 1 x , y x ,
where

Z(G; x, y) := m ∑ k=0 k ∑ l=0 (-1) k b kl x k y l .
Thus, y n Z(G; y -1 , 1) and (x + 1) n Z(G; (x + 1) -1 , (x + 1) -1 ) give the chromatic polynomial and the independence polynomial of G, respectively.

As an example, consider the path P 4 and the star S 3 in Figure 1. Let b kl and b ′ kl denote the coefficients of P(P 4 ; x, y) and P(S 3 ; x, y), respectively. According to the interpretation of the coefficients provided by Theorem 9 we find that

b 00 = 1, b 22 = 1, b ′ 00 = 1, b ′ 22 = 0, b 10 = 0, b 30 = 0, b ′ 10 = 0, b ′ 30 = 0, b 11 = 3, b 31 = 1, b ′ 11 = 3, b ′ 31 = 1, b 20 = 0, b 32 = 0, b ′ 20 = 0, b ′ 32 = 0, b 21 = 2, b 33 = 0, b ′ 21 = 3, b ′ 33 = 0.
Putting these values into (9) we again obtain (3) and ( 4). The corresponding Z-polynomials are

Z(P 4 ; x, y) = 1 -3xy + 2x 2 y + x 2 y 2 -x 3 y , Z(S 3 ; x, y) = 1 -3xy + 3x 2 y -x 3 y .
5 Special graphs

The complete graph

The lattice Π G for the complete graph G = K n coincides with the partition lattice Π(V ) of V , where V denotes the vertex set of G. From (5) we obtain the equation

P(K n ; x, y) = ∑ π∈Π(V ) µ( 0, π)x k 1 (π) y |π|-k 1 (π) , ( 12
)
where µ is the Möbius function of the partition lattice as defined in Eq. [START_REF] Read | An introduction to chromatic polynomials[END_REF]. A closer look at Eq. ( 12) reveals that the terms of the sum do not depend on the partition π but only on its type. In the following, we use a second representation for a number partition λ. Besides λ = (λ 1 , . . . , λ r ) we also write λ = (1

k 1 2 k 2 ...n k n )
where each k i gives the number of parts of size i in λ. Consequently,

n ∑ i=1 k i = |λ| and n ∑ i=1 ik i = n ,
and hence,

P(K n ; x, y) = n! ∑ λ⊢n (-1) n-|λ| ∏ |λ| i=1 λ i ∏ n i=1 k i ! x k 1 y |λ|-k 1 .
A simpler representation of this formula follows immediately from Corollary 7:

P(K n ; x, y) = n ∑ k=0 n k (x -y) k y n-k .

The complete bipartite graph

Let K m,n = (V ∪ W, E) be the complete bipartite graph with m + n = |V | + |W | vertices. We choose k vertices from V and color these vertices with proper colors. Let k j denote the Stirling number of the second kind. For j different colors there are k j y j different colorings. The remaining mk vertices of V may be colored with improper colors. No vertex of W may receive one of the j proper colors that are used for V . Hence, we count (xj) n different colorings for vertices of W . Summing up all possibilities, we obtain

P(K m,n ; x, y) = m ∑ k=0 m k (x -y) m-k k ∑ j=0 k j y j (x -j) n .
As a special case, we obtain for S n = K 1,n :

P(S n ; x, y) = x n (x -y) + y(x -1) n .

The path

Let P n be the path with n vertices. As an example, we first consider the path P 5 as shown in Figure 4. The lattice of forbidden colorings for P 5 is shown in Figure 5. We observe that this lattice is isomorphic to the Boolean lattice B 4 . More generally, the lattice Π P n is isomorphic to the Boolean lattice B n-1 (in fact, Π T is isomorphic to B n-1 for any tree T having n vertices). The proof of this fact is rather easy. A connected partition of the vertex set of the path P n with k blocks is generated by choosing k -1 separators Fig. 5: The lattice of forbidden colorings for P 5 between the elements of the linearly ordered set {v 1 , ..., v n }. A separator corresponds to a removed edge of the path. In the following, we refer to such a partition as a linear partition. Now, by Theorem 4,

P(P n ; x, y) = ∑ i, j (-1) n-i-j l i j x i y j ,
where l i j denotes the number of linear partitions π of v 1 , v 2 , . . . , v n with i singletons blocks and j nonsingleton blocks. All such π can be constructed as follows. Start with a diagram of i + j -1 slashes which create i + j spaces (between, before and after slashes) in which to place dots. Place one dot in i of the spaces and two dots in the remaining j. This can be done in i+ j i ways. Now distribute ni -2 j dots among the spaces that already have two dots. The number of ways of doing this is the number of ways of choosing ni -2 j objects from j objects with repetition which equals

j + (n -i -2 j) -1 n -i -2 j = n -i -j -1 n -i -2 j .
Now replace the n dots from left to right with the labels v 1 , v 2 , . . . , v n . The result is a linear partition with the specified number of singleton and non-singleton blocks and every such partition is obtainable in this way. Note that the total number of choices is

l i j = i + j i n -i -j -1 n -i -2 j
as desired. Hence,

P(P n ; x, y) = ∑ 0<i+2 j≤n (-1) n-i-j i + j i n -i -j -1 n -i -2 j x i y j .
For example, the polynomials of paths of up to five vertices are P(P 1 ; x, y) = x, P(P 2 ; x, y) = x 2y, P(P 3 ; x, y) = x 3 -2xy + y, P(P 4 ; x, y) = x 4 -3x 2 y + y 2 + 2xyy, P(P 5 ; x, y) = x 5 -4x 3 y + 3xy 2 + 3x 2 y -2y 2 -2xy + y .

Section 7 gives a method for obtaining the generalized chromatic polynomial of any tree by using a polynomial-time recursive algorithm.

The cycle

The last special graph to be considered here is the cycle C n . The construction of the lattice of forbidden colorings differs only slightly from the construction of the corresponding lattice for the path. Let V = {v 1 , ..., v n } be the vertex set of C n in correspondence with the order of traversing the cycle. A cyclic partition π of V is one obtained by removing edges from C n and using the resulting connected subsets of vertices as blocks of π. Note that the lattice of cyclic partitions of V is isomorphic to B n with the rank just below the maximum element removed. By Theorem 4,

P(C n ; x, y) = (n -1)(-1) n-1 y + ∑ (i, j) =(0,1) (-1) n-i-j c i j x i y j ,
where c i j denotes the number of cyclic partitions which have i singletons and j non-singletons. These can be obtained as follows. As in Section 5.3, build a dot diagram which can be done in l i j ways where l i j again denotes the number of linear partitions. Now turn the diagram into a cyclic partition by picking one of the n dots and then replacing the dots starting with that dot as v 1 , moving right, and then wrapping around circularly to get the dots to the left of v 1 . Thus we get nl i j labeled diagrams. We claim that every cyclic partition is obtained exactly i + j times. This is because any circular partition can be written with the element v 1 in any of the i + j blocks. So nl i j = (i + j)c i j and hence,

c i j = n i + j i + j i n -i -j -1 n -i -2 j .
Therefore,

P(C n ; x, y) = (-1) n y + n ∑ 0<i+2 j≤n (-1) n-i-j i + j i + j i n -i -j -1 n -i -2 j x i y j .
For example, the generalized chromatic polynomial for cycles of up to five vertices are P(C 3 ; x, y) = x 3 -3xy + 2y, P(C 4 ; x, y) = x 4 -4x 2 y + 4xy + 2y 2 -3y, P(C 5 ; x, y) = x 5 -5x 3 y + 5xy 2 + 5x 2 y -5y 2 -5xy + 4y.

Non-isomorphic graphs

Does the generalized chromatic polynomial distinguish non-isomorphic graphs? Two non-isomorphic graphs having the same generalized chromatic polynomial must have the same chromatic polynomial, the same matching polynomial, and the same independence polynomial. Consequently, the number of vertices, edges, and components must be the same for each such pair. The smallest non-isomorphic graphs having the same generalized chromatic polynomial are depicted in Figure 6. We found by complete enumeration that there is no pair of non-isomorphic trees with up to nine vertices that have the same generalized chromatic polynomial. We found such a pair with 10 vertices which is presented in Figure 7.

Note that the chromatic symmetric function, as defined by STANLEY [START_REF] Stanley | Enumerative Combinatorics[END_REF], does not coincide for these two trees. The question whether the chromatic symmetric function distinguishes non-isomorphic trees is still open. We hope that our polynomial could be of some value to decide this question. The search for pairs of non-isomorphic trees having the same chromatic symmetric function can now be restricted to those pairs of trees that have the same generalized chromatic polynomial. In this way (and using a computer) we found that there are no non-isomorphic trees with up to 15 vertices that have the same chromatic symmetric function.

As with the usual chromatic polynomial, the computation of the generalized chromatic polynomial is an NP-hard (in fact: #P-complete) counting problem. As a consequence, polynomial-time algorithms can only be found for some restricted classes of graphs. In this section, we consider the restricted class of trees.

Let P(G,W ; x, y) be the generalized chromatic polynomial of G with the additional restriction that there is one proper color forbidden for all vertices of some given subset W of V . Let C(G, v) be the set of components of Gv. The following theorem provides an algorithm for computing P (G; x, y) for any tree G.

Theorem 10 Let T be a tree and v ∈ V (T ). Then, for each T i ∈ C(G, v) the set N(v) ∩ V (T i ) consists of only one vertex, and the polynomial P(T ; x, y) can be computed via the following recursion:

P(T ; x, y) = (x -y) ∏ T i ∈C(T,v) P(T i ; x, y) + y ∏ T i ∈C(T,v) P(T i , N(v) ∩V (T i ); x, y) , P(T, {v}; x, y) = (x -y) ∏ T i ∈C(T,v) P(T i ; x, y) + (y -1) ∏ T i ∈C(T,v) P(T i , N(v) ∩V (T i ); x, y) .

If T only consists of v, then

P(T ; x, y) = x, P(T, {v}; x, y) = x -1 .
Proof. Consider the first equation. Each vertex of V can be colored with proper or improper colors. There are exactly xy possibilities to choose an improper color for v in which case the coloring of the neighbors of v is independent of the color of v. Consequently, v can be removed from G. The number of remaining colorings follows from the product formula [START_REF] Brylawski | The broken circuit complex[END_REF]. The second product of the first equation arises from the y different colorings of v with proper colors. Now, the neighbors of v have to be colored differently from v. The neighbor of v in the subtree T i is given by N(v) ∩ V (T i ). Indeed, this set contains exactly one vertex. The number of remaining colorings of the subtree T i are given by P(T i , N(v) ∩V (T i ); x, y). This shows the first formula. In order to prove the second formula we first remark that there are again xy different colorings using only improper colors, which cause no restriction on the components of Tv. However, if we color v with a proper color, then we have only y -1 possibilities. The remaining part is proved in a similar way. ✷ Figure 8 illustrates the computation of the generalized chromatic polynomial for an example tree. The darkened vertices in this figure correspond to the vertices for which one color is forbidden. The decomposition vertex at each step is denoted by v. The resulting polynomial for the depicted tree is:

P(T ; x, y) = (x -y)x 3 (x -y) x 2 + y (x -1) 2 + y(x -1) 3 (x -y) x 2 + (y -1)(x -1) 2 = x 7 -6x 5 y + 9x 4 y + 6x 3 y 2 -11x 3 y -9x 2 y 2 + 10x 2 y + 5xy 2 -5xy -y 2 + y.
There exists a vertex v in every tree with n vertices such that no component of Tv has more than n/2 vertices. This vertex can be found in polynomial time as follows. Start with a randomly chosen decomposition vertex w. If Tw contains a component T i with more the n/2 vertices, then next choose

v y x -y x x - y y - x -y x -y ! ! x x - v v
Fig. 8: Computation of the generalized chromatic polynomial of a tree the neighbor of w belonging to T i . Repeating this procedure, we obtain after at most n/2 steps a maximum component with at most n/2 vertices. By the appropriate combination of components, we can always achieve a decomposition into three subtrees each consisting of at most n/2 vertices. Let f (n) be the time complexity of the recursive algorithm as described in Theorem 10 for the computation of the generalized chromatic polynomial of any input tree of size n = |V |. We find that

f (n) ≤ 3 f n 2 + p(n) , (13)
where p(n) is a polynomial that incorporates the effort of the search for a decomposition vertex. For each polynomial p(n) inequality (13) leads to solutions which are polynomial in n. Consequently, we obtain a polynomial time algorithm for the computation of the generalized chromatic polynomial of a tree:

Corollary 11

The generalized chromatic polynomial of a tree can be computed in polynomial time.

Corollary 12

The generalized chromatic polynomial of a path P n , n > 1, satisfies the recurrence relation

P (P n ; x, y) = (x -y) P (P n-1 ; x, y) + y(y -1) n-2 (x -1) + y (x -y) n-1 ∑ i=2 (y -1) i-2 P(P n-i ; x, y)
with the initial condition P(P 1 ; x, y) = x.

Corollary 12 is an immediate consequence of Theorem 10.

Graphs of restricted pathwidth

Many NP-hard graph problems become simple in graphs of restricted treewidth or pathwidth. A nice introduction to treewidth and related problems is given by BODLAENDER [START_REF] Bodlaender | A partial k-arboretum of graphs with bounded treewidth[END_REF]. Here we present a polynomialtime algorithm for the computation of the generalized chromatic polynomial for graphs of restricted pathwidth. The main ideas may also be extended to cover graphs of restricted treewidth. The problem of determining the pathwidth (treewidth) is NP-hard in general. If the pathwidth is known to be at most k then there is a polynomial-time algorithm for finding a path decomposition (composition order) of width k.

A path decomposition of a graph G = (V, E) is a sequence J = (J 0 , ..., J r ) of subsets of V such that

1. r i=0 J i = V , J 0 = J r = / 0,
2. every edge of G has both its vertices contained in at least one subset of J , 3. for all i, j, k with 0

≤ i < j < k ≤ r the relation J i ∩ J k ⊆ J j holds.
The width of a path decomposition J is the maximal cardinality of a subset of J . The pathwidth of the graph G is the minimum width over all path decompositions of G. A path decomposition J = (J 0 , ..., J r ) is called nice if |(J i+1 J i ) ∪ (J i J i+1 )| = 1 for all i = 0, ..., r -1. A nice path decomposition J may be described in terms of a sequence S (J ) of signed vertices of G: The notation +v indicates the inclusion of the vertex v into the set J i while the notation -v stands for the removal of v from J i . For example, a nice path decomposition of the bridge graph presented in Figure 9 is

J = ( / 0, {1}, {1, 2}, {1, 2, 3}, {2, 3}, {2, 3, 4}, {3, 4}, {4}, / 0) .
which is uniquely determined by the following sequence of signed vertices S (J ) = (+1, +2, +3, -1, +4, -2, -3, -4) .

It follows from the definition of path decomposition that each vertex appears exactly once with positive sign and exactly once with negative sign in S(J ). Let J = (J 0 , ..., J r ) a path decomposition of width k of G. Then J can be transformed in O(n) time into a nice path decomposition J of width at most k (see [START_REF] Scheffler | Die Baumweite von Graphen als ein Maß für die Kompliziertheit algorithmischer Probleme[END_REF] and [START_REF] Thilikos | Constructive linear time algorithms for small cutwidth and carving-width[END_REF]).

A composition order of a graph G = (V, E) is a sequence s = (s 1 , ..., s t ) of signed vertices and of edges of G such that 1. if s i ∈ s represents an edge with end vertices u and v then there are indices j, k, m, p with j, k < i < m, p and

s j = u, s k = v, s m = -u, s p = -v,
2. every edge of G appears exactly once in s,

3. the removal of all edges from s yields a nice path decomposition of G.

A composition order of the graph in Figure 9 is

s = (+1, +2, {1, 2}, +3, {1, 3}, -1, {2, 3}, +4, {2, 4}, -2, {3, 4}, -3, -4) . ( 14
)
We refer to the objects in a composition order as steps or, to be more precise, as vertex activation steps (positively signed vertices), vertex deactivation steps (negatively signed vertices), and edge insertion steps. The composition order gives complete information about the graph structure, that is, no other graph representation is required for the computation of the generalized chromatic polynomial. Corollary 7 shows that the numbers a pq of independent partitions consisting of p singletons and q other blocks give all the information which is necessary to determine the generalized chromatic polynomial of G. The algorithm presented below computes these numbers directly. Some additional notation is needed. Let G i be the graph consisting of all activated vertices and all edges inserted up to and including step i in a composition order s, let V i be the vertex set of G i , and A i be the set of active vertices of G i , i.e., V i minus the set of vertices deactivated up to and including step i. A labeled partition π of A i is a set partition of A i where blocks may have a label that will be denoted by a star ( * ).

For each independent partition σ of V i we define the induced index of σ as a triple I = (π, p, q) consisting of a labeled partition π of A i and two integers p and q where p equals the number of singletons of σ and q is the number of blocks of σ with more than one vertex. The blocks of π correspond to the restriction of σ to A i , and a block of π is labeled iff the corresponding block of σ is not a singleton. We note that different partitions of V i may have the same induced index and refer to these partitions as represented by the index. Let c be the number of partitions of V i represented by the induced index I. A pair (I, c) will be called a state consisting of the index I and the value c, and the set of all states at step i defines the state set Z i .

The main idea behind the algorithm is to compute the state set Z i from the previous state set Z i-1 according to certain rules depending on the step s i of the given composition order.

Obviously, the very first state set Z 0 consists only of a single state (( / 0, 0, 0), 1), since V 0 = / 0 and since there is only one partition of an empty set, and this partition contains neither singleton nor non-singleton blocks. More interestingly, the very last state set Z r contains all the numbers a pq . The set of active vertices will be empty since each vertex has been activated and deactivated once. The partitions in the induced indices are partitions of this active set and hence empty partitions. The values c corresponding to indices I = (π, p, q) gives the number a pq of partitions of the full vertex set (since V 2|V |+|E| = V ) with p singleton and q non-singleton blocks.

The transformation rules for the state sets are as follows: For each step s i , all states of the state set Z i-1 are considered in turn leading to some (possibly empty) set of intermediate successor states. These states represent the same set of partitions as the original state but use the new active vertex set for their description. After the individual states are transformed, there might be sets of such intermediate successor states sharing the same index. These sets are replaced by single states using the common index and a value equal to the sum of the values of states in the set. The states remaining after this clean-up form the next state set Z i .

Vertex activation (+v): After activation step s i activating some vertex v, this vertex v is a singleton in G i , since no edge incident to v could have been introduced before v became active. So v is independent of all vertices in G i-1 and hence, it can extend any independent set of G i-1 . Suppose we need to transform a state z = (I, c) with I = (π, p, q). Let S = {σ 1 , σ 2 , . . .} be the set of independent partitions consisting of p singletons and q non-singleton blocks that induce I. We then proceed as follows:

• First of all, note that the new vertex v can form a singleton block of its own in any of the σ k .

In this case, the number of singletons increases by one, and v forms a new unlabeled block of π. Everything else remains unchanged, leading to a single intermediate successor state of ((π|v, p + 1, q), c), where π|v is the partition obtained by augmenting π by a single block {v}.

• Secondly, since v is independent of all vertices in G i-1 , it is certainly independent of the vertices that constitute any of the q non-singleton blocks of some σ k . Hence, v can be added to any such block, which results in a block of cardinality larger than one. Let X be such a non-singleton block and define Y := X ∩ A i-1 . We have to distinguish two cases: (a) If Y = / 0, then v will be a singleton in (X + {v}) ∩ A i . Since |X ∪ {v} | > 1 it must be marked, and thus the single successor state is ((π|v * , p, q), c). (b) If Y = / 0, then we already have parts of X in the active vertex set, namely the vertices of Y . In this case, the successor state is ((π -Y | (Y ∪ {v}) * , p, q), c).

• Similarly to the above, the new vertex v can be added to all p singleton blocks. This works in a similar way with one exception: the singleton grows to a non-singleton block, whence we have to decrement p and increment q by one in the successor state.

Vertex deactivation (-v): If I = (π, p, q) describes a set of independent partitions of G i-1 , then (πv, p, q) describes the same set for G i where πv is π -{v} in case that v is a singleton resp. π -X|X {v} in case that v ∈ X for a larger block X of π. So for deactivation, every state has a single intermediate successor state, which is obtained by removing the vertex for the partition in the index partition.

Edge insertion {u, v}: Let I, w and π be defined as in the vertex activation step above. If π is a partition containing a block X with {u, v} ⊆ X then this block forms no independent set in G i as G i results from G i-1 by inserting edge {u, v}. Consequently, the state (I, c) has no successor states as there are no independent partitions inducing π in G i .

As an example, we compute the generalized chromatic polynomial of the bridge graph depicted in Figure 9. The resulting transformations using the composition order given by ( 14) are presented in Tab. 1: Computation of the generalized chromatic polynomial of the graph in Figure 9.

1. At the end, two states remain corresponding to the coefficients a 40 = 1 and a 21 = 1. The polynomial in equation ( 8) is thus x 4 + x 2 y, which can easily be transformed via Corollary 7 into P (G; x, y) = x 4 -5x 2 y + 6xy -4y + 2y 2 .

A minor modification of the presented algorithm permits the computation of a function Q * (G; x, y) that is defined analogously to Q (G; x, y) with the exception that the coefficient a pq counts the number of partitions consisting of p singletons and q other blocks such that each block forms a clique.

The computation of Q * (G; x, y) again requires a nice path decomposition of G. However, the edge insertion step is not necessary. A state z = (π, f ) consists of a labeled partition π as index and a bivariate polynomial f as value. This time, the label of a block of π indicates that this block belongs to a clique that contains already deactivated vertices. Such a labeled block can not be extended to a larger clique. The coefficient [x p y q ] f of the polynomial f gives the number of partitions represented by the index π that consist of p singletons and q other cliques. The initial state for the empty graph is ( / 0, 1). We use the following elementary transformations:

Vertex activation (+v): A new vertex may always form a one-vertex clique. Thus, every state (π, f ) induces a new state (π|v, x f ). If v is adjacent to a unlabeled singleton {w} of π then {v, w} forms a two-vertex clique of G producing a new state π -{w} | {v, w} , x -1 y f . If v is adjacent to all vertices of a non-singleton block W that is not labeled then an existing clique may be enlarged by v generating a new state (π -W |W ∪ {v} , f ).

Vertex deactivation (-v):

The vertex v has to be removed from the index (i.e. the labeled partition). If v is removed from a non-singleton block then the remaining set becomes a labeled block. The values remain unchanged. A reduction of the number of states is possible by collecting states with identical index into one state the value of which is the sum of the values of the collection.

Fig. 3 :

 3 Fig. 3: Lattice of forbidden colorings for the graph in Figure 2.

)

  For any subset I ∈ BC (G) let m I , n I and c I denote the number of edges, vertices and connected components of the edge-subgraph G[I], respectively. Since G[I] is cycle-free for any I ∈ BC (G), we conclude that m In I + c I = 0 and hence, i∈I

Fig. 4 :

 4 Fig. 4: The path P 5

Fig. 6 :

 6 Fig. 6: Non-isomorphic graphs having the same generalized chromatic polynomial.

Fig. 7 :

 7 Fig. 7: Non-isomorphic trees having the same generalized chromatic polynomial.

Fig. 9 :

 9 Fig. 9: Bridge graph.

Table Step Partition

 Step

	+1	/ 0 1	p q Count 0 0 1 1 0 1	Step Partition p q Count 24 * |3 2 1 1 2|34 * 2 1 1
	+2	1|2 12 *	2 0 0 1	1 1	{2, 4}	2|3|4 2|3|4 *	4 0 2 1	1 1
	{1, 2} 1|2	2 0	1		2|34 *	2 1	1
	+3	1|2|3 13 * |2 23 * |1	3 0 1 1 1 1	1 1 1	-2	3|4 3|4 * 34 *	4 0 2 1 2 1	1 1 1
	{1, 3} 1|2|3 23 * |1	3 0 1 1	1 1	{3, 4}	3|4 3|4 *	4 0 2 1	1 1
	-1 {2, 3} 2|3 2|3 23 * +4 2|3|4	3 0 1 1 3 0 4 0	1 1 1 1	-3 -4	4 4 * / 0 / 0	4 0 2 1 4 0 2 1	1 1 1 1
		2|3|4 *	2 1	1

Acknowledgement

The authors wish to express their gratitude to the referees for the many helpful suggestions that resulted in an improvement of the paper. In particular, the derivations of the generalized chromatic polynomial of a path (Subsection 5.3) and a cycle (Subsection 5.4) were greatly simplified by one of the referees, whose solution we adopted with minor changes and whom we especially would like to thank for his worthy contribution.

The application of this algorithm to the graph of Figure 9 yields

Let Ḡ be the complement of the graph G. It follows from the definition of Q * that Q (G; x, y) = Q * Ḡ; x, y . This gives the following result:

Theorem [START_REF] Tutte | Graph Theory[END_REF] The generalized chromatic polynomial can be computed in polynomial time for graphs of restricted pathwidth and for graphs whose complement is of restricted pathwidth. In particular, the matching polynomial, the independence polynomial, and the chromatic polynomial can be computed in polynomial time for graphs of restricted pathwidth and for graphs whose complement is of restricted pathwidth.


Charles Knessl

Numerical Studies of the Asymptotic Height Distribution in Binary Search Trees

Keywords:

) for the tails of the distribution. The asymptotic height distribution is shown to be unimodal with highly asymmetric tails.

Introduction

A binary search tree is a fundamental data structure used in searching and sorting. It is defined as follows. There are n elements to be stored in the tree. A root node is created for the first element. Subsequent elements are directed to the left or right subtree according to whether they are less than or greater than the element in the root node. By this construction, the left and right subtrees are also binary search trees by themselves. Many classic sorting algorithms (such as QUICKSORT) can be conveniently represented by binary search trees (BST).

It is well known that the worst search time for this model is O(n), but the average search time is only O(log n). We consider the average case performance and introduce the following probabilistic model.

We take all n! permutations of the n elements to be equally likely and analyze the height H n of a BST constructed from n elements. The height is the longest path in the randomly built tree. Clearly H n cannot exceed n and must exceed log 2 n. In view of the probabilistic assumption H n is a random variable and we set L k n = Prob{H n ≤ k}. The support lies in the range k ≤ n < 2 k . There has been a lot of previous work on computing various aspects of this probability distribution, in the limit n → ∞. [START_REF] Pittel | On growing random binary trees[END_REF] showed that (almost surely) H n log n → A 0 as n → ∞, with A 0 ≤ A = 4.31107 . . . . [START_REF] Devroye | A note on the height of binary search trees[END_REF] established that E[H n ] ∼ A log n as n → ∞. This was refined to E[H n ] = A log n + O(log log n) by [START_REF] Devroye | On the variance of the height of random binary search trees[END_REF].

In the past it had been conjectured that E[H n ] -A log n ∼ -δ A

A-1 log log n with δ = 1 2 , but recent results of [START_REF] Reed | How tall is a tree[END_REF] (see also [START_REF] Knessl | The height of a binary search tree: the limiting distribution perspective[END_REF]) show that the correct value is δ = 3 2 . There has also been some work on estimating the variance var[H n ]. Experimental studies of [START_REF] Robson | The height of binary search trees[END_REF] show 1). This has been established rigorously by [START_REF] Drmota | An analytic approach to the height of binary search trees[END_REF][START_REF] Drmota | The variance of the height of binary search trees[END_REF] and by [START_REF] Reed | How tall is a tree[END_REF].

that E H n -E[H n ] is bounded, suggesting that var[H n ] = O(
Very little seems to be known about the full distribution L k n . In [START_REF] Knessl | The height of a binary search tree: the limiting distribution perspective[END_REF] we used singular perturbation methods to analyze a recurrence relation satisfied by the distribution. Under some assumptions about the forms of various asymptotic expansions, we obtained expressions for L k n and 1-L k n for n → ∞ and various ranges of k. In the range where most of the probability mass accumulates, we showed that L k n can be approximated by the solution of a non-linear integral equation. This was related to a functional equation studied by [START_REF] Drmota | An analytic approach to the height of binary search trees[END_REF]. We established some asymptotic properties of the solution to this integral equation, but could not solve it exactly. Recently, [START_REF] Drmota | An analytic approach to the height of binary search trees II[END_REF] established rigorously that the height distribution function satisfies this integral equation, in the limit n → ∞. It was also shown that the equation has a unique solution that satisfies a certain auxiliary condition (cf. ( 4) and ( 7)).

In this note we supplement the results of [START_REF] Knessl | The height of a binary search tree: the limiting distribution perspective[END_REF] by numerically analyzing the integral equation. We thus obtain the shape of the asymptotic height distribution numerically. We state the problem more precisely in section 2, and the numerical results are discussed in section 3 and in the Figures and Tables therein.

Problem Statement

Let us denote by H n the height of a binary search tree that stores n elements. Its probability distribution

L k n = Prob{H n ≤ k} (1)
satisfies the non-linear recurrence

L k+1 n+1 = 1 n + 1 n ∑ ℓ=0 L k ℓ L k n-ℓ
(2) subject to the initial condition L 0 n = δ(n, 0). Setting

z = k -A log n + 3 2 A A -1 log log n + c (3)
and assuming that L k n ∼ f (z) we derived in [START_REF] Knessl | The height of a binary search tree: the limiting distribution perspective[END_REF] the following non-linear integral equation for f (z):

f (z + 1) = Z 1 0 f (z -A log x) f z -A log(1 -x) dx, -∞ < z < ∞ (4) f (-∞) = 0, f (∞) = 1.
In (3) A = 4.31107 . . . is the unique solution to (x/2) x = e x-1 in the range x > 1. We observe that if f 0 (z) is a solution to (4) then any translation (i.e., f 0 (z + C)) is also a solution. Thus by retaining in the right side of (3) the arbitrary constant c, we can choose a convenient way to normalize the solution to (4) so as to make it unique. Recently, [START_REF] Drmota | An analytic approach to the height of binary search trees II[END_REF] established rigorously that (4) has a unique solution, modulo the translation.

We let f (z) = 1g(z) where clearly g(z) will be small for z → ∞. Then we can approximate, for z large, g(z) ≈ g L (z) where g L satisfies the linearized equation

g L (z + 1) = 2 Z 1 0 g L (z -A log x)dx = 2 Z ∞ 0 g L (z + At)e -t dt.
(5) Now ( 5) admits exponential solutions of the form e -νz provided that ν satisfies the characteristic equation

e -ν = 2 1 + νA . ( 6
)
We can easily show that ν = 1 -1/A is a double root of ( 6) and that this is the only real solution for ν > 0. There exist infinitely many complex solutions to ( 6) and these can be used to construct solutions to ( 4), e.g., by the method of successive iterations. However, the numerical and analytic studies in [START_REF] Knessl | The height of a binary search tree: the limiting distribution perspective[END_REF] show that these lead to solutions that are inappropriate (they typically oscillate and/or become negative). This again follows more rigorously from the work of [START_REF] Drmota | An analytic approach to the height of binary search trees II[END_REF]. Thus we write the general "acceptable" solution to (5) as

g L (z) = exp -1 - 1 A z (αz + β). (7)
Knessl and Szpankowski ( 2002) also showed that if α = 0 then the solution to the non-linear problem (4) becomes negative for -z sufficiently large. Thus we have α > 0 and in view of (3) normalize our solution by setting α = 1. Now we use ( 7) to construct a solution to the non-linear problem, with f (z) = 1g(z), in the form

g(z) = (z + β)e -az + ∞ ∑ m=2 e -maz P m (z) (8) = ∞ ∑ m=1 e -maz P m (z), a = 1 - 1 A ,
where P m (z) is a polynomial of degree m; we write

P m (z) = m ∑ j=0 F(m, j)z j . (9)
Using ( 8) in ( 4) leads to e -ma P m (z + 1) -2

Z 1 0 x maA P m (z -A log x)dx (10) = - m-1 ∑ ℓ=1 Z 1 0 x ℓaA (1 -x) (m-ℓ)aA P ℓ (z -A log x)P m-ℓ (z -A log(1 -x))dx,
for m ≥ 2. Then by using (9) and comparing coefficients of z M we are led to

m ∑ J=M J M F(m, J) e -m(1-1/A) - 2(J -M)!A J-M [1 + m(A -1)] J-M+1 (11) = -[z M ] m-1 ∑ ℓ=1 ℓ ∑ k 1 =0 m-ℓ ∑ k 2 =0 z k 1 +k 2 ℓ ∑ i=k 1 m-ℓ ∑ j=k 2 i k 1 j k 2 ×D(ℓ, m -ℓ, i -k 1 , j -k 2 )F(ℓ, i)F(m -ℓ, j).
Here [z M ] denotes the coefficient of z M (we may replace

[z M ]z k 1 +k 2 by δ(k 1 + k 2 , M)) and D(α, β, γ, δ) = A γ+δ Z 1 0 x α(A-1) (1 -x) β(A-1) (-log x) γ [-log(1 -x)] δ dx. ( 12
)
The above can be expressed in terms of derivatives of the Beta function. We also note that

F(1, 1) = 1 and F(1, 0) = β.
Approximations to g(z), and hence f (z), may be obtained simply by truncating the sum in (8) at some large value m = N. However, this leads to problems for z negative and -z sufficiently large, as discussed in section 3. We also note that given P 1 (z) = z + β, each F(m, j) is a polynomial in β in view of ( 11).

Finally we represent f (z) in the contour integral form

f (z) = 1 2πi Z i∞ -i∞ exp η A e -z/A F (η)dη, (13)
which says that, after an appropriate variable change, F (η) is the (two-sided) Laplace transform of f (z).

Then from (4) it follows that

-F ′ (η) = e -2/A [F (ηe -1/A )] 2 . ( 14
)
This is a functional-differential equation studied by [START_REF] Drmota | An analytic approach to the height of binary search trees[END_REF][START_REF] Drmota | The variance of the height of binary search trees[END_REF], who used the normalization condition F (0) = 1, with which (14) has a unique analytic solution about η = 0, that is in fact an entire function (this is rigorously shown in [START_REF] Drmota | An analytic approach to the height of binary search trees[END_REF][START_REF] Drmota | The variance of the height of binary search trees[END_REF]). We note that our normalization (which took α = 1 in ( 7)) is different from Drmota's; unfortunately it seems that neither can be used to infer the true value of c in (3). An important difference is that while (14) has a unique solution, our problem still has a one-parameter infinity of solutions, with β indexing the family. However we show numerically that only one value of β leads to a solution that can satisfy the condition f (-∞) = 0 (g(-∞) = 1). The other solutions grow very rapidly as z → -∞ and apparently do not have Laplace transforms. Hence they are excluded from ( 14) by the form (13).

In [START_REF] Knessl | The height of a binary search tree: the limiting distribution perspective[END_REF] we also established that as

z → -∞ f (z) satisfies f (z) ∼ 2 2κ π √ A log 2 A log 2 -1 e -ωz/2 exp(-κe -ωz ) (15)
where ω = log 2 A log 2 -1 = .3486294060 . . . and κ is a constant, which is made unique by choosing α = 1. Our derivation of (15) made certain assumptions about the asymptotic form; the numerical studies here provide more justification for this and also estimate the constant κ.

Numerical Results

We define

f N (z) = 1 - N ∑ m=1 e -amz m ∑ j=0 F(m, j)z j , N ≥ 1 (1)
with g N (z) = 1f N (z). These correspond to approximate solutions to (4). The exact solution must also satisfy f (-∞) = 0 and f ′ (z) > 0 for all z. Some of the problems arising in the convergence of f N (z) to f (z) are illustrated by discussing f N for a particular N, and we consider N = 6 in detail below.

A plot of f 6 (z) = f 6 (z; β) for various β shows that typically both f 6 and f ′ 6 have zeros and hence lead to unacceptable approximations to a probability distribution. Our goal is to define a criteria and choose an optimal value of β that somehow minimizes this "unacceptability". Then we shall increase N and obtain a sequence of optimal β that converges to the unique β for which we have f (-∞) = 0 in (4).

In Table 1 we give the zeros of both f 6 and f ′ 6 for various values of β. We note that for general N we have, as z → -∞, f N (z) ∼ -e -aNz F(N, N)z N . We can easily show that F(m, m) > 0 for all m ≥ 1 and thus as z → -∞ f N (z) → +∞ (resp. -∞) for N odd (resp. N even). The data in Table 1 show that f 6 has exactly three zeros if β ∈ (-∞, β ′ ) or β ∈ (β * , ∞), and a single zero if β ∈ (β ′ , β * ). Here β ′ ∈ (-1.5, -1.4) and β * ∈ (-.9112, -.9111). β > β, with β ∈ (-1.4, -1.3). We define an optimal β as follows. For a given β we consider the minimum value of z such that f 6 (z) and f ′ 6 (z) are both positive for all z exceeding this value. More precisely we let z (6) * (β) = max{z : f 6 (z) = 0} subject to the constraints that f 6 > 0 and f ′ 6 > 0 for z > z (6) * (β). Then β opt is defined as the value of β that minimizes z (6) * (β). Note that z (6) * (β) may or may not exist for a particular β. When N = 6 Table 1 shows that it exists for all β exceeding about -.9111. We can define a more general z (N) * (β) by setting it equal to the largest value of z where f ′ N vanishes (if z (N) * (β) fails to exist by the previous definition). In either case f N (z) can be an acceptable approximation to a probability distribution only for z > z (N) * (β). Our computational experience has shown that β opt always corresponds to two roots of f N coalescing into a double root (and thus a root of f ′ N ). When N = 6, β opt ≈ -.9112 and with this value f 6 has a double zero at z min ≈ -2.637. Also, if we plot f 6 (z + 1)f 6 (z), which would be an approximation to Prob{H n = k + 1}, we find that it has a zero at z min ≈ -2.9917. Thus f 6 (z + 1)f 6 (z) is an acceptable approximation to a probability density (or discrete distribution) only for z > z min .

f N (z) z = -1 z = -2 z = -3 z = -4 N =
For arbitrary N we again compute β opt , z min and z min . These are summarized in Table 2 for N ≤ 15.

The data suggest that β opt converges rapidly to the value -.9119760 . . . min | decreasing with N. The sequence of functions h N (z) = f N (z + 1)f N (z) is converging to some unimodal positive function h(z), with a maximum value of h max = .1743 . . ., which is attained at z = z max = .2918 . . . . The convergence near z = z max is very rapid, as illustrated in Table 3. As z becomes negative, the convergence of f N (z) becomes much slower as |z| increases. Also, there is a lot of cancelation in the sum in (1); for N = 14 and 15 we had to increase the precision to 20 digits in order to accurately do the calculation. In Table 4 we illustrate the convergence of f N (z) for (fixed) negative values of z. The value N = 15 is not sufficient to see convergence at z = -5. The minimum value of z that sees f 15 (z) settling to its limit is about z = -4.3. We find that f (-4.1) = .00058157 . . ., f (-4.2) = .00042828 . . . and f (-4.3) = .00031059 . . . . We also see from Table 4 that once f N settles to its limit value, it does so very quickly.

Next we test the asymptotic formula (15), which applies for z → -∞. The difficulties described above preclude us from computing f (z) for large negative values. The constant κ in (15) could not be determined analytically. We can estimate it simply by numerically computing f (z 0 ) for a certain negative z 0 (our results allow only z 0 ≥ -4.3), comparing this to the right side of (15) (with z = z 0 ), and solving the resulting transcendental equation for κ. In Table 5 z Fig. 1: The approximation to f (z) using N = 15, for z ∈ (-5, 8) κ ≈ 2.13. It should be noted that the three values for z < -4 are more sensitive to error. Also, f (z) should be very small due to the doubly-exponential convergence to zero as z → -∞. However, when z = -4, -ωz ≈ 1.4 which is not particularly large. A really accurate calculation of κ (and verification of ( 15)) would probably require that we calculate f (z) accurately for values much more negative than -4.3.

In Figure 1 we plot f 15 (z) for z ∈ [-5, 8] and in Figure 2 we plot f 15 (z + 1)f 15 (z) for the same range. These are our final approximations to f (z) and f (z + 1)f (z). The second figure clearly illustrates the shape of the "density", showing its unimodal structure, the (roughly) exponential right tail and the very thin (roughly double exponential) left tail. These figures used β = -.9119760150. Note that we are approximating the discrete distribution (1) (or L k+1 n -L k n = Prob{H n = k + 1}) by the continuous function h(z). For a given large n we can choose several values of k in (3) to make z = "O(1)" and the corresponding values of L k n should lie close to the curve in Figure 1, for some appropriate c. The values of L k+1 n -L k n should then lie close to the curve in Figure 2 for this value of c. We have no analytic method for estimating the value of c. In [START_REF] Knessl | The height of a binary search tree: the limiting distribution perspective[END_REF] it was shown that if we had an asymptotic approximation to L k n valid on the scale k, n → ∞ with k/ log n fixed and > A, then we could use asymptotic matching to infer the value of c. However, we could not completely analyze this scale, which we refer to as the "near right tail" of the distribution. There 1 -L k n is algebraically small in n (for a fixed k/ log n ∈ (A, ∞)).

To summarize, we have presented an efficient numerical method for calculating the asymptotic height distribution in binary search trees. Our results yield the distribution's shape, but there is still the arbitrary translation arising from c in (3). Our results also suggest that the non-linear integral equation has, up to a translation, a unique solution that can represent a probability distribution. This was recently established rigorously by [START_REF] Drmota | An analytic approach to the height of binary search trees II[END_REF]. If α = 1 in (7) this solution corresponds to β = -.9119760 . . . . There are still many issues regarding the convergence of L k n to f (z) that need further work.

  min are converging to -∞, but much more slowly, with the "gaps" |z

Fig. 2 :

 2 Fig. 2:The approximation to h(z) = f (z + 1)f (z) using N = 15, for z ∈(-5, 8)

  The derivative f ′ 6 has four zeros if β < β and two zeros if

	N	β opt	z min	z min
	6 7 8 9 10 11 12 13 14 15	-.9111950 -.9117765 -.9119242 -.9119624 -.9119724 -.9119750693 -.9119757674 -.9119759527 -.9119760019 -.9119760150	-2.638 -3.052 -3.425 -3.764 -4.074 -4.360 -4.625 -4.873 -5.105 -5.324	-2.991 -3.398 -3.763 -4.097 -4.403 -4.686 -4.950 -5.196 -5.427 -5.645
	N	z max	h max = h(z max )
	4	.3092	.1741
	5	.2922	.1743
	6	.2918	.1743
	8	.2918	.1743
	10	.2918	.1743

Tab. 2: β opt , z min and z min for N ≤ 15.

Tab. 3: Convergence near z = z max .

  we estimate κ using various z. It would appear that

		1
		0.8
		0.6
		f(z)
		0.4
		0.2
	-4	-2	2	4	6	8
	z	f (z)		κ
	-1 -2 -3 -4 -4.1 -4.2 -4.3	.12953 .042753 .0083862 .00077929 .00058157 .00042828 .00031059	2.0495 2.0898 2.1100 2.1219 2.1236 2.1257 2.1287
		Tab. 5: Estimating κ using various z.

† This work was supported by NSF Grants DMS-99-71656 and DMS-02-02815, and NSA Grant MDA 904-03-1-0036. 1365-8050


Wei-Mei Chen

Hsien-Kuei Hwang

Tsung-Hsi Tsai

Efficient maxima-finding algorithms for random planar samples

Keywords: Maximal points, computational geometry, probabilistic analysis of algorithms, Pareto optimality, sieve algorithms, dominance knapsack problems: see [2, 27, 39], 7. longest common subsequence problem: see [3, 37], 8. data mining: see [13, 31], 9. graph drawing: see [7, 40], 10. heuristic search in artificial intelligence: see [1, 22, 49, 57], 11. analysis of linear programming: see [15, 16], 12. transportation: see [18, 20, 21], 13. gene finding: see [32, 38], 14. statistical decision theory: see [55, 63], 15. ecological process models: see [54, 61]

We collect major known algorithms in the literature for finding the maxima of multi-dimensional points and provide a simple classification. Several new algorithms are proposed. In particular, we give a new maxima-finding algorithm with expected complexity n + O( √ n log n) when the input is a sequence of points uniformly chosen at random from general planar regions. We also give a sequential algorithm, very efficient for practical purposes.

Introduction

Finding the maxima of a given point set is a fundamental problem in computational geometry, which also arises in many applications such as Pareto-optimality in bargaining games or multicriteria optimization, analysis of linear programming, analysis of unbounded knapsack problem, statistical decision theory, etc.; see below for a brief list of instances. In the special case when d = 1, it reduces to finding the maximum value of n numbers. We propose in this paper a simple classification of several known algorithms for finding the maxima, together with several new algorithms; among these are two efficient algorithms-one with expected complexity n + O( √ n log n) when the point samples are issued from some planar regions, and another more efficient than existing ones.

A point p = (x 1 , . . . , x d ) in R d is said to dominate another point q = (y 1 , . . . , y d ) if x i > y i for all i = 1, . . . , d. For convenience, we write p ≻ q. The nondominated points in a point set are called maximal points or maxima. Similar terms like Pareto optimality (or Pareto front, Pareto solutions, etc.), admissibility, efficiency, and noninferiority have been used in other fields. Dominance of multidimensional points is a natural and one of the most widely used partial orders in diverse fields such as 1. multicriteria decision analysis: see [START_REF] Greco | Rough approximation by dominance relations[END_REF][START_REF] Pomerol | Multicriterion Decision in Management[END_REF]; 2. Pareto optimality in multicriteria optimization: see [START_REF] Ehrgott | Multicriteria Optimization[END_REF][START_REF] O'neill | The number of outcomes in the Pareto-optimal set of discrete bargaining games[END_REF]; 3. evolutionary computation: see [START_REF] Van Veldhuizen | Multiobjective evolutionary algorithms: analyzing the state-of-the-art[END_REF][START_REF] Zitzler | Proceedings of the First International Conference on Evolutionary Multi-Criterion Optimization[END_REF]; 4. VLSI design: see [START_REF] Kathail | PICO: automatically designing custom computers[END_REF][START_REF] Lillis | Timing optimization for multisource nets: characterization and optimal repeater insertion[END_REF][START_REF] Sun | An area minimizer for floorplans with L-shaped regions[END_REF]; Apart from the preceding classification, we are mainly interested in this paper in algorithms of incremental type, which are in general simple and have good expected performance. We briefly describe in the next section known algorithms for d = 1, and then give a corresponding version for higher dimensions. Our new algorithms for general planar regions are given in Sections 3 and 4. A comparative discussion of the latter algorithm with the list algorithm of Bentley et al. [START_REF] Bentley | Fast linear expected-time algorithms for computing maxima and convex hulls[END_REF] is also given. We then discuss various aspects of different design paradigms in the last section.

2 Maxima-finding algorithms: from d = 1 to d > 1

How to find the maximum of n numbers? We collect below several simple algorithms and state a natural extension to finding the maxima in higher dimensions.

Sequential algorithms. The most straightforward algorithm for finding the maximum value of n numbers has the following incremental, on-line, one-loop pattern:

Algorithm SEQUENTIAL //input = {a 1 , . . . , a n }, a j ∈ R, j = 1, . . . , n// m := a 1 for i := 2 to n do if a i > m then m := a i //output m = max{a 1 , . . . , a n }//
Obviously, the number of comparisons used by this algorithm is n -1. Such a pattern has a natural extension to multidimensional inputs (see [START_REF] Bentley | Fast linear expected-time algorithms for computing maxima and convex hulls[END_REF][START_REF] Kung | On finding the maxima of a set of vectors[END_REF]):

Algorithm SEQUENTIAL(d) //input = {a 1 , . . . , a n }, a j ∈ R d , j = 1, . . . , n, d ≥ 2// M := {a 1 } for i := 2 to n do if a i dominates some points in M then delete these points; M := M ∪ {a i } elseif ∀b ∈ M, a i is incomparable with b then M := M ∪ {a i } //output M =maxima{a 1 , . . . , a n }//
Here we say that two points p and q are incomparable if neither p dominates q nor q dominates p. Kung et al. [START_REF] Kung | On finding the maxima of a set of vectors[END_REF] used balanced search trees to maintain the maxima found so far, while Bentley et al. [START_REF] Bentley | Fast linear expected-time algorithms for computing maxima and convex hulls[END_REF] used a list to implement the steps inside the for-loop, and apply the move-to-front heuristic to the dominating point in the list of maxima when new-inserting points are dominated. We will see later that an implementation using random search trees turns out to be very efficient in practice.

Divide-and-conquer algorithms. Another simple way to find the maximum of n numbers is to use the half-half divide-and-conquer paradigm:

Algorithm DIVIDE-AND-CONQUER //input {a 1 , . . . , a n }, a j ∈ R, j = 1, . . . , n// maximum{a 1 , . . . , a n } (n ≥ 1) if n = 1 then return a n m 1 := maximum{a 1 , . . . , a ⌊n/2⌋ } m 2 := maximum{a ⌊n/2⌋+1 , . . . , a n } if m 1 ≥ m 2 then return m 1 else return m 2 //output max{a 1 , . . . , a n }//
It is easily seen that the number of comparisons used is also n -1. Indeed, the solution to the recurrence f n = f λ(n) + f n-λ(n) + 1 with f 1 = 0 and λ(n) ∈ {1, . . . , n -1} satisfies f n = n -1. Thus the half-half dividing rule can be replaced by other dividing rules without changing the (optimum) worst-case complexity.

This design pattern, although less useful for d = 1, turned out to be very fruitful for higher dimensions; see [START_REF] Bentley | Multidimensional divide-and-conquer[END_REF][START_REF] Bentley | On the average number of maxima in a set of vectors and applications[END_REF][START_REF] Bentley | Divide and conquer for linear expected time[END_REF][START_REF] Devroye | Moment inequalities for random variables in computational geometry[END_REF][START_REF] Flajolet | Exact asymptotics of divide-and-conquer recurrences[END_REF][START_REF] Kung | On finding the maxima of a set of vectors[END_REF][START_REF] Preparata | Computational Geometry: An Introduction[END_REF].

Bucket algorithms. The top-down description of the divide-and-conquer algorithm suggests another algorithm, bucket-select, for finding selected order statistics in a given sequence; see Mahmoud et al. [START_REF] Mahmoud | Analytic variations on bucket selection and sorting[END_REF]. The generalization of this algorithm to finding the maxima in higher dimensions is also straightforward; see Devroye [START_REF] Devroye | A note on the expected time required to construct the outer layer[END_REF][START_REF] Devroye | Lecture Notes on Bucket Algorithms[END_REF]. No description of these bucket algorithms is given here since we are not discussing them in this paper. Note that unlike sequential and divide-and-conquer algorithms (deterministic in nature), the performance of bucket algorithms is more sensitive to the distribution of the input.

Selection algorithms. Yet another algorithm for finding the maximum of n numbers is quickselect when always selecting the largest order statistics:

Algorithm QUICKSELECT //input P := {a 1 , . . . , a n }, a j ∈ R, j = 1, . . . , n// quickselect(P) (|P| ≥ 1) if |P| = 1 then return P x := a "pivot" element from P Q := {p ∈ P : p > x} if |Q| = 0 then return x else quickselect(Q) //output max{a 1 , . . . , a n }//
Although the worst-case of this algorithm is quadratic, its expected number of comparisons (assuming that the input is a sequence of independent and identically distributed random variables with a common continuous distribution) is equal to 2n -2 ∑ 1≤ j≤n 1/ j; see Knuth [START_REF] Knuth | Mathematical analysis of algorithms[END_REF].

This simple pattern can be easily extended to finding the maxima in R 2 . One natural pattern is the following sieve algorithm.

Algorithm SELECTION //input = P := {a 1 , . . . , a n }, a j ∈ R 2 , j = 1, . . . , n// selection(P)
if |P| = 1 then return P choose a maximal point, say x in P as the sieve discard points dominated by x M := {x} for each nonempty quadrant

Q of x do M := M∪ selection(Q) return M //output M = maxima{a 1 , . . . , a n }//
The expected complexity of this algorithm can be quickly seen from considering the special case when the points are drawn uniformly at random, one independently from the other, from the right triangle with corners (0, 0), (0, 1) and (1, 0). A good sieve, suggested by the probabilistic analysis in [START_REF] Bai | Limit theorems for the number of maxima in random samples from planar regions[END_REF], is the point, say x, that maximizes the sum of both coordinates. The algorithm thus proceeds by first locating x, and then continues recursively for points lying in the second and fourth quadrants of x, respectively. The expected number of comparisons µ n satisfies µ 0 := 0 and the recurrence (see [START_REF] Bai | Limit theorems for the number of maxima in random samples from planar regions[END_REF])

µ n = βn + c + (n -1)! Γ(n + 1/2) ∑ 2≤ j<n Γ( j + 1/2) j! µ j (n ≥ 1),
(assuming that βn + c comparisons on average are used to split the problem into two smaller problems) where Γ denotes the Gamma function. The exact solution is given by (see [START_REF] Bai | Limit theorems for the number of maxima in random samples from planar regions[END_REF])

µ n = 3βn + 2β -c -(2β -c) √ π n! Γ(n + 1/2) (n ≥ 1) = 3βn + O( √ n).
Thus the expected complexity of such a selection algorithm is linear. The limit distribution of the cost can also be derived by the approach used in [START_REF] Bai | Limit theorems for the number of maxima in random samples from planar regions[END_REF] if more information on the "splitting toll cost" is available. Note that the probabilistic analysis is in general less transparent for random points from other regions. For multidimensional points, a similar algorithm can be devised but an additional merge as that used in divide-and-conquer algorithms is required.

Another sequential selection algorithm is discussed in Clarkson [START_REF] Clarkson | More output-sensitive geometric algorithms[END_REF] in which maxima are found one after another by the following steps: compare each incomparable element p thus far to the set of maxima (initially empty); drop p if it is dominated; otherwise p remains incomparable and is then used to find a new maximal element by scanning through the remaining incomparable elements, updating p if p is dominated (again, dominated elements are dropped in the process). The procedure is reminiscent of the selection sort. Clarkson's algorithm works for all dimensions, but becomes less useful if maxima are abundant.

Sieve algorithms.

Selection algorithms spend some effort in finding first a good sieve point; when more information on the input is known (in addition to the dimension, like the underlying distribution), one may quickly identify a sieve point with almost no additional cost, and with good sieving power. A naive pattern for d = 1 is as follows.

Algorithm SIEVE //input P := {a 1 , . . . , a n }, a j ∈ R, j = 1, . . . , n// M := / 0 choose a sieve point, say x //x ∈ P// for i := 1 to n do if a i > x then M := M ∪ {a i } if |M| > 0 then
find the maximum value in M else use other maximum-finding algorithm //output max{a 1 , . . . , a n }// If the underlying distribution of the input is Uniform(0, 1), then we can take x = 1 -(log n)/n, so that M is empty with probability proportional to n -1 . If that really happens, then we can use any algorithm like SEQUENTIAL to find the maximum from zero, and the expected cost so introduced is bounded. The algorithm has thus an overall expected cost satisfying n

+ O(log n) if SEQUENTIAL is used when |M| > 0.
Such a sieve-based approach is, although not useful for d = 1, easily rewritten for higher dimensional problem:

Algorithm SIEVE(d) //input P = {a 1 , . . . , a n }, a j ∈ R d , j = 1, . . . , n, d ≥ 2// M := / 0
choose a sieve point, say x //x ∈ P// for i := 1 to n do if a i is not dominated by x then M := M ∪ {a i } if all maxima are included in M then find these maxima else use other maxima-finding algorithm //output maxima{a 1 , . . . , a n }//

The two other algorithms proposed in [START_REF] Bentley | Fast linear expected-time algorithms for computing maxima and convex hulls[END_REF] share the same pattern as above, and as good sieves for random points issuing from [0, 1] d , they take

1. x = (1 -δ n , . . . , 1 -δ n ) ∈ R d , where δ n = n -1/d (log n) 1/d ; 2. x = (ξ [1] n , . . . , ξ [d] n ) ∈ R d , where ξ [i]
n denotes the ⌊nδ n ⌋-th largest element in the i-th dimension. Both choices of p guarantee that all but O(n 1-1/d (log n) 1/d ) points on average are sieved out by p when the points are independently and identically chosen from [0, 1] d . To check the condition (or certificate) "all maxima included in M", one can divide points in M into two sets when filtering the points: points dominating p and points incomparable with p. Then the condition "all maxima included in M" holds if there are points dominating p; otherwise, we can resort to any algorithm with even O(n 2 ) worst-case time bound.

A new sieve algorithm

The efficiency of sieve algorithms relies heavily on the quality of the sieve, which in turn depends on the underlying distribution. By construction, both sieve points for SIEVE(d) given in [START_REF] Bentley | Fast linear expected-time algorithms for computing maxima and convex hulls[END_REF] are based on the assumption that the underlying shape (for the input points) is a hypercube (or, more generally, having an upper-right corner † ; see [START_REF] Bai | Limit theorems for the number of maxima in random samples from planar regions[END_REF]). It is obvious that when the underlying planar shape does not have an upper-right corner such as a circle or a triangle of the form ❅ ❅ , no point inside the region has the power of dominating all but O(n 1-ε ) points.

We propose a new sieve algorithm based on the same design pattern as SIEVE(d) but using a curve as a sieve instead of a point, so that almost all points can be discarded. While the idea is straightforward, finding a good certificate for assuring that all maxima are in the sieved region is non-trivial.

Let X 1 , X 2 , ..., X n be independent and identically distributed random vectors in R 2 . Assume that X 1 is uniformly distributed in a region

D := {(x, y) : 0 ≤ x ≤ 1, 0 ≤ y ≤ f (x)}, (1)
where f (x) is a nonincreasing function on [0, 1] with f (0) > 0.

The reason of considering the shape D is that for general planar regions most maxima lie in the upperright part; see [START_REF] Bai | Limit theorems for the number of maxima in random samples from planar regions[END_REF] for more details. Also the case when f (x) = 1x is a prototype for general regions; see [START_REF] Bai | Berry-Esseen bounds for the number of maxima in planar regions[END_REF]. Our result on D can be easily extended to general planar regions R provided that the boundary of the upper-right part of R has the form of D.

Without loss of generality, we may assume that |D| = 1 and that

f (x) > 0, for 0 ≤ x < 1; = 0, for x ≥ 1. Let S := {(x, y) ∈ D : y + δ < f (x + δ)};
see Figure 3. If δ > 0 is suitably chosen, then there is no maxima in S and the number of points in .5 ). Then we apply the bucket algorithm (see [START_REF] Devroye | Lecture Notes on Bucket Algorithms[END_REF]) for points in D -S and read off the maximal points there in a suitable order; see below for more details.

D -S is O( √ n log n) with probability 1 -o(n -1
To state the algorithm precisely, we need the following lemma.

Lemma 1 Let Z := (x, y) ∈ D : y + δ 2 ≥ f x + δ 2 ,
and α := inf x :

f x + δ 2 ≤ δ 2 . If (0, α) ⊂ [ (x i ,y i )∈Z x i - δ 4 , x i + δ 4 , (2)
where X i = (x i , y i ), then there is no maxima in S.

Proof. Suppose that there is a maximum in S. Then there exists a point (x 0 , y 0 ) such that y 0 + δ = f (x 0 + δ) and there is no point in the first quadrant of (x 0 , y 0 ). We now find a point (x ′ 0 , y 0 ) ∈ ∂Z, i.e.

y 0 + δ 2 = f (x ′ 0 + δ 2 ). Since f (x 0 + δ) > f (x ′ 0 + δ 2 )
and f is nonincreasing,

x ′ 0 -x 0 > δ 2 ;
and none of the points {X 1 , . . . , X n } is in In view of this lemma, the pattern of our algorithm is then as follows.

Z ∩ {(x, y) : x 0 < x < x ′ 0 }. These contradict (2). ✷ S Z (x0, y0) δ δ δ/2 δ/2 (x ′ 0 , y0) f (x) α
Algorithm CURVE-SIEVE // input = {a 1 , . . . , a n }, a j ∈ R 2 , j = 1, . . . , n// M := / 0 [START_REF] Andonov | Unbounded knapsack problem: Dynamic programming revisited[END_REF] holds then find these maxima else use other maxima-finding algorithm //output maxima{a 1 , . . . , a n }// Here a i .x and a i .y denotes respectively the x and y coordinates of a i .

for i := 1 to n do if a i .y + δ > f (a i .x + δ) then M := M ∪ {a i } if
Assume that comparisons of the type a i .y + δ > f (a i .x + δ) takes unit cost. Let C n denote the number of comparisons used by the algorithm CURVE-SIEVE for finding the maxima of {X 1 , . . . , X n }.

Theorem 1 Take δ = 4 log n n .
Then the average number of comparisons used by the CURVE-SIEVE algorithm satisfies

E(C n ) = n + O n log n .
Proof. The number of comparisons used to test if a point lies in S or D -S is n. Assume at the moment that the expected number of comparisons used in checking ( 2) is proportional to the number of points in D -S. Then we have

C n = n + O(the number of points in D -S), if (2) holds; O n 2 , otherwise.
Now we show that the probability of the event that (2) fails is less than o(n -1.5 ). Let m = ⌈4α/δ⌉ + 1 and define the strips

Z i := Z ∩ {(x, y) : δ(i -1)/4 < x < δi/4}, 1 ≤ i ≤ m.
If (2) fails, then there is at least a strip Z j , j ∈ {0, 1, . . . , m}, in which there is no point. Observe that due to our choice of α, the height of Z i is no less than δ/2, and thus |Z i | ≥ δ 2 /8. It follows that the probability that (2) fails is bounded above by

∑ 1≤i≤m (1 -|Z i |) n ≤ m 1 -δ 2 /8 n = o(n -1.5 ).
Observe that

|D -S| = Z 1 0 ( f (x) -max{0, f (x + δ) -δ}) dx ≤ (1 + f (0))δ = O log n n , (3)
where the O-term depends on f (0), which is finite by our assumption. Since the probability of the event that the number of points in

D -S is more than O( √ n log n) is less than o(n -1 ), we have E(C n ) = n + O n log n 1 -o(n -1 ) + O (n) o(n -1 ) + O n 2 o(n -1.5 ) = n + O n log n .
We still need to show that (i) the maxima in D -S can be found in time O( √ n log n), and (ii) the condition (2) can be checked within the same time bound.

The problem (ii) is easy: we divide the unit interval into ⌈2/δ⌉ buckets, find points in each interval with the maximum and minimum x-coordinates, respectively, and then the condition (2) can be easily checked. This is essentially the so-called maximal gap problem; see [25, p. 85].

For problem (i), we use the bucket argument used in Devroye [START_REF] Devroye | Lecture Notes on Bucket Algorithms[END_REF] by dividing the region D -S into small buckets each of which has at most O(1) expected number of points; this is achievable by uniformity of the input points in D and by [START_REF] Apostolico | The longest common subsequence problem revisited[END_REF]. We scan the strips from right to left and for the buckets in each strip we go from top to bottom, checking for the maximal points and each time keeping track of the topmost nonempty bucket in each strip (so that we need only to check buckets with higher values of y-coordinates for the strip left to the current one).

✷ List algorithm vs tree algorithm. While a general analysis of the expected number of comparisons used by our tree algorithm is not easy, there is a special case when the analysis becomes straightforward: it is the case when C tree n (A) ≤ C list n (A), where A is the set of input points, C tree n denotes the number of comparison used by our tree algorithm, and C list n that by the list algorithm of Bentley et al. when the "move-to-front" heuristic is not used. When the input points are taken independently and identically at random from the unit square, Devroye [START_REF] Devroye | A note on the expected time for finding maxima by list algorithms[END_REF] showed that Simulations also suggested that the number of comparisons has large variance (of the order of mean square), so that larger scale simulations are needed in order to obtain more meaningful approximations.

Special cases of f (x). (i) When f (x) = 1 for 0 ≤ x ≤ 1, our sieve equation becomes y > 1 -δ, if x < 1 -δ; = 0, if 1 -δ ≤ x ≤ 1. (ii) if f (x) = 2(1 -x), then α = 1 -3δ/4
E(C list n ) = n + o(n),

Comparative discussions

We briefly discuss the preceding algorithms from several different angles. Maxima-finding vs sorting. Knuth [45] classified major sorting algorithms by the main operations used: insertion, exchanging, selection, merging, and distribution. An alternative classification, similar to that of Knuth's, of maxima-finding algorithms is given in Table 3.

Paradigms

Sorting before divide-and-conquer vs sorting after divide-and-conquer. Sorting is usually used before divide-and-conquer to achieve better worst-case performance. The scenario changes when considering average-case complexity because in most cases the number of maxima is at most O(n 1-ε ) on average, and thus sorting after divide-and-conquer is preferable, especially when the number of maxima is small like log n or √ n.

Probabilistic vs deterministic. Once a good sieve for sieve algorithms is selected, we may completely neglect the case when the "certificate" fails, giving a probabilistic algorithm that finds all maxima with high probability. This adds to the flexibility of sieve algorithms. However, the certificate like (2) may in general not easy to be devised.

Known vs unknown underlying distribution. Different cost is needed when finding the maximum of n numbers with given information. Typical examples include:

1. Maximum-finding without any a priori information on the input requires n -1 comparisons;

2. Zero comparison is needed if one knows that the given input is a sorted sequence;

3. Maximum-finding takes Θ( √ n) expected time if the input is a random walk; see Odlyzko [START_REF] Odlyzko | Search for the maximum of a random walk[END_REF];

4. A probabilistic algorithm (with an error probability O(n -c ) of finding the maximum) with cost satisfying O((log n) 2 ) is given in Ting and Yao [START_REF] Ting | A randomized algorithm for finding maximum with O((log n) 2 ) polynomial tests[END_REF];

5. O(log n) expected time is needed for finding the maximum value in a binary search tree.

Our sieve algorithm CURVE-SIEVE assumes that the input points are randomly taken from some planar region with about O( √ n) maxima on average; are there good sieve algorithms when this a priori information of the input points is not available? Worst-case vs typical case. The two-dimensional maxima problem may be regarded as an easy problem as far as worst-case complexity is concerned because one can simply sort the points according to their x-coordinates and then read off their y-coordinates from right to left, looking for right-to-left maximal elements. This results in an O(n log n) algorithm. Our main interests here lie in maxima-finding algorithms with linear expected cost and of practical value.

There are data structures in the literature capable of maintaining dynamically (insertion and deletion) the current maxima in worst-case time bound log n, but their practical usefulness is unclear; see [START_REF] Kapoor | Dynamic maintenance of maxima of 2-d point sets[END_REF] and the references therein. For other notions like output-sensitive and distribution-sensitive algorithms, see [START_REF] Kirkpatrick | Output-size sensitive algorithms for finding maximal vectors[END_REF][START_REF] Sen | Distribution-sensitive algorithms[END_REF]. Again the main focus of these algorithms is on worst-case instead of average-case as we emphasized in this paper.

Square vs other regions. Most probabilistic analysis of maxima-finding algorithms assume that the input is uniformly distributed in a unit square. The results so obtained only reflect the tip of an iceberg since the number of maxima is very different for other general planar regions. Also efficient randomized algorithms for square may become less efficient for other inputs. See Fukuda et al. [START_REF] Fukuda | Interval finding and its application to data mining[END_REF] for a concrete situation where different distributions arise.

  1365-8050 c 2003 Discrete Mathematics and Theoretical Computer Science (DMTCS), Nancy, France Tab. 1: A classification of known maxima-finding algorithms according to the design pattern for d = 1.

Fig. 1 :

 1 Fig. 1: Dissection of the region D.

  and the sieve equation is y < 2(1x) -3δ for x ≤ α. //input: b ∈ {left, right} and T := a tree to be properly pruned. prune(b, T ) check := T next := check.b while (next = ∅) do if (next.key lies in the third quadrant of check.key) then next := next.b remove check.b and check.b's b subtrees //left := right and right := left// else b := b next := check.b //output: more compact T

  and the same time bound applies to E(C tree

				n	).
	n	list	list+MTF	tree		n	list	list+MTF	tree
	10	1.499	1.307	1.316		10	2.076	1.792	1.655
	10 2 1.808	1.560	1.406		10 2	4.695	3.683	2.442
	10 3 1.573	1.309	1.187		10 3	7.249	4.901	2.567
	10 4 1.370	1.201	1.069		10 4	9.672	5.237	2.543
	10 5 1.261	1.136	1.028		10 5 12.629	7.019	2.512
	10 6 1.198	1.084	1.010		10 6 14.351	9.189	2.513
	10 7 1.137	1.025	1.003		10 7 17.585	10.918	2.517
	10 8 1.095	1.016	1.0002	10 8 21.326	13.866	2.576
	Tab. 2: Average number of comparisons per element (by simulations) when the underlying region is the unit square
	(left) and the right triangle with corners (0, 0), (0, 1) and (1, 0) (right); MTF denotes "move-to-front."

† Roughly, an upper-right corner of a bounded region R is a point on the boundary of R dominating all points in and on R.

Acknowledgements

We thank Ralph Neininger and Michael Fuchs for useful suggestions.

Extension to higher dimensions. The same sieving procedure can be applied to higher dimensions with suitable modifications, but the bucket argument (used to guarantee a linear expected time bound) does not extend easily. One may apply the divide-and-conquer algorithm (see [START_REF] Kung | On finding the maxima of a set of vectors[END_REF][START_REF] Preparata | Computational Geometry: An Introduction[END_REF]) instead of bucket algorithm for points not sieved out, but the error term becomes slightly bigger O(n 1-1/d (log n) c ) for some constant c > 1/d depending on d.

A sequential tree algorithm

Our description of SEQUENTIAL(d) is completely general in that we can implement the steps inside the for-loop by different types of data structures. The original implementation by Bentley et al. [START_REF] Bentley | Fast linear expected-time algorithms for computing maxima and convex hulls[END_REF] used a list (together with the move-to-front heuristic), and it is proved by Devroye [START_REF] Devroye | A note on the expected time for finding maxima by list algorithms[END_REF] that the expected number of vector comparisons used by the list algorithm (with or without the move-to-front rule) satisfies n + o(n) if the input points are taken uniformly at random from [0, 1] d , one independently of the others; see also Clarkson [17], Golin [START_REF] Golin | A provably fast linear-expected-time maxima-finding algorithm[END_REF] for d = 2.

When the underlying distribution changes from the hypercubes to other regions, simulations suggest that the list algorithm may become less efficient, although rigorous analysis is still lacking; see [START_REF] Bentley | Fast linear expected-time algorithms for computing maxima and convex hulls[END_REF][START_REF] Golin | A provably fast linear-expected-time maxima-finding algorithm[END_REF]. We propose in this section an implementation using random search trees for d = 2 (instead of balanced search trees as in [START_REF] Kung | On finding the maxima of a set of vectors[END_REF]). The tree is essentially a quadtree, but without the first and the third quadrants. The algorithm is very efficient in practice and robust to different distributions, but the extension to higher dimensions and the analysis of its expected cost are not clear and left open.

//output: T := the tree containing all maximal vectors of P //input: A point p i and T := the tree containing all maximal vectors of {p