Effective Asymptotics of Combinatorial Systems
Carine Pivoteau, Bruno Salvy
TL;DR
This paper extends analytic combinatorics methods to general systems of combinatorial equations, enabling the derivation of asymptotic expansions for a broader class of generating functions under certain conjectural assumptions.
Contribution
It generalizes existing results to all combinatorial systems, providing an almost complete algorithmic framework for asymptotic analysis of their generating functions.
Findings
Applicable to algebraic-logarithmic singularities
Works under Schanuel's conjecture for certain systems
Extends known polynomial system results to general systems
Abstract
Analytic combinatorics studies asymptotic properties of families of combinatorial objects using complex analysis on their generating functions. In their reference book on the subject, Flajolet and Sedgewick describe a general approach that allows one to derive precise asymptotic expansions starting from systems of combinatorial equations. In the situation where the combinatorial system involves only cartesian products and disjoint unions, the generating functions satisfy polynomial systems with positivity constraints for which many results and algorithms are known. We extend these results to the general situation. This produces an almost complete algorithmic chain going from combinatorial systems to asymptotic expansions. Thus, it is possible to compute asymptotic expansions of all generating functions produced by the symbolic method of Flajolet and Sedgewick when they have…
| Bell, Burris, Yeats, well conditioned [2] | |||
|---|---|---|---|
| Strongly connected graph | has no zero-coordinates | ||
| Drmota [9, 10] | |||
| Strongly connected graph | Zero-free system and | ||
| Bell, Burris, Yeats, elementary [3] | |||
| Pivoteau, Salvy, Soria, well founded at 0 [38] | nilpotent | ||
| Zero-free system | |||
| Pivoteau, Salvy, Soria, well founded [38] | Let ; well founded at 0, plus a finiteness condition on its solution | ||
| Banderier, Drmota, analytically well defined [1] | has spectral radius | entire, no polynomial coordinate in | |
| This article, well founded | well defined | nilpotent | |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsAdvanced Combinatorial Mathematics · Polynomial and algebraic computation · Limits and Structures in Graph Theory
Effective Asymptotics of Combinatorial Systems
Carine Pivoteau
Bruno Salvy
Abstract
Analytic combinatorics studies asymptotic properties of families of combinatorial objects using complex analysis on their generating functions. In their reference book on the subject, Flajolet and Sedgewick describe a general approach that allows one to derive precise asymptotic expansions starting from systems of combinatorial equations. In the situation where the combinatorial system involves only cartesian products and disjoint unions, the generating functions satisfy polynomial systems with positivity constraints for which many results and algorithms are known. We extend these results to the general situation. This produces an almost complete algorithmic chain going from combinatorial systems to asymptotic expansions. Thus, it is possible to compute asymptotic expansions of all generating functions produced by the symbolic method of Flajolet and Sedgewick when they have algebraic-logarithmic singularities (which can be decided), under the assumption that Schanuel’s conjecture from number theory holds. That conjecture is not needed for systems that do not involve the constructions of sets and cycles.
This article is dedicated to the memory of Philippe Flajolet, who shared with us his vision of the interaction between analytic combinatorics and computer algebra.
Contents
-
2 Well-Founded Combinatorial Systems and their Characterizations
-
2.2 Constructive Systems in the Symbolic Method of Flajolet & Sedgewick
-
2.3 Characterization of Well-founded Systems by Leading Terms
-
2.5 Characterization of Well-founded Systems by Jacobian Matrices
-
3.4 Radius of Convergence of Constructible Generating Functions
-
5 Radius of Convergence of Constructible Generating Functions using Newton Iteration
-
7 Dominant Singularities of Constructible Generating Functions
1 Introduction
In their book Analytic Combinatorics [19], Flajolet and Sedgewick describe a general approach that starts from a combinatorial description, translates this description into equations satisfied by generating functions, views these generating functions as analytic functions and exploits their singular behavior to deduce asymptotic properties of the combinatorial objects when their size becomes large. The aim of this article is to develop computational tools that automate large parts of this approach. We first outline the main steps and the computational issues that they raise.
Symbolic Method and Well-Founded Systems
The first step of the approach is called the ‘Symbolic Method’.111“The symbolic method is the ‘combinatorics’ in analytic combinatorics” [19, p. 92]. It consists of a precise definition of the type of combinatorial descriptions that can be considered, in the form of a system of recursive equations relating combinatorial families of objects. The general form of such a system is
[TABLE]
that we denote concisely by . Precise statements on what the can be are given in Section 2. For instance, permutations could be described by corresponding to their cycle decomposition, binary trees by the recursive equation and general (Cayley) trees by . In order to automate this approach, a very first step is to determine when such a system is meaningful, excluding systems like that would not have a useful combinatorial meaning. We define a system to be well founded when the iteration started with converges (in a combinatorial sense, see Definition 2.2 and the following paragraphs). A simple algorithm (Algorithm 1) tests this condition. This algorithm is more efficient than our earlier one in [38]; it can be seen as an evolution of Zimmermann’s algorithm [48, Algorithm B, p. 33], where the computation of valuations and the detection of cycles in the graph of dependence of the system are combined in a single iteration.
This combinatorial setting is presented in Section 2, using a small amount of vocabulary from species theory that makes it easier to discuss the well-founded character of combinatorial systems. The combinatorial objects that can be defined with the constructions considered by Flajolet and Sedgewick are called constructible species, to distinguish them from the general case and the corresponding systems are called constructive. The algorithms presented in this article are generally stated for constructible species and their defining systems, while the mathematical results that prove their correctness are often given for arbitrary systems.
Generating Function Systems
The next step of the approach is to introduce generating functions. In this article, we focus on the case of labeled structures and exponential generating functions. This means that the object of interest is the power series
[TABLE]
where is a vector whose th coordinate is the number of objects of type with distinct atoms labeled from 1 to . For constructible species, the translation from to a system of generating function equations is straightforward thanks to an explicit dictionary going from combinatorial constructions to elementary functions [19, Chap. II], see Table 1.
One can also define well-founded systems directly at the level of generating function equations. This allows us to compare our definition with those studied by Banderier and Drmota [1, 9, 10] and Bell, Burris and Yeats [2, 3], that have extra constraints meaningful in the contexts of their analyses. This is presented Section 3.1, which can be used as a starting point by readers who want to focus on the analytic aspects.
Radius of Convergence
Next comes the analytic part of analytic combinatorics. The first step is to determine the radius of convergence of the generating functions, viewed as analytic functions, as this governs the exponential growth of their coefficients. For a given inside the domain of convergence of , the system has as a solution, i.e., the value of the generating function at . It may also have other solutions with nonnegative coordinates inside the domain of convergence of . However, for larger than the radius of convergence of , we prove that no such solution exists (Proposition 3.11). This is a basis for the computability of the radius of convergence by dichotomy, provided one can detect the existence of nonnegative solutions to systems (Theorem 3.14). By the Tarski-Seidenberg theorem, this existence is decidable when the system is a polynomial system. Combinatorially, this corresponds to a context-free grammar, or a combinatorial system that does not involve Cyc or Set. In the more general situation where Cyc and Set are allowed, a result of Macintyre and Wilkie [32] shows that the decision can also be performed, and therefore the radius of convergence computed, provided a conjecture in transcendental number theory called Schanuel’s conjecture holds. As this may well be stronger than needed, we specify more precisely an oracle whose decisions are sufficient for the computation of the radius of convergence, and to which all the decisions in our algorithms reduce.
Newton’s Iteration
Looking more closely at the structure of the systems of generating function equations allows for more precise results, presented in Section 4. Notably, these systems exhibit strong positivity properties, very similar to those considered by Etessami, Stewart, Yannakakis [15, 14, 42] for polynomial systems. This has an important consequence for Newton’s iteration: for where is the radius of convergence of , Newton’s iteration started at converges to [38]. Theorem 4.1 is a converse of this result: if Newton’s iteration started at converges to a point with nonnegative coordinates inside the domain of convergence of , then and . In practice, one still needs to determine when to stop the iteration: it can well be that and that the current iterate is close to but it could also happen that and the ultimate divergence of the iteration has not been detected yet. For this purpose, a Newton-Kantorovich-type theorem bounds the distance of the last iterate to the solution (Theorem 4.4); positivity of the system allows this bound to be obtained from a single evaluation, instead of requiring an upper bound on a disk. Theorem 4.7 bounds the speed of convergence in the spirit of a similar result by Esparza, Kiefer, Luttenberger [13].
Newton’s iteration can also be used to find the radius of convergence. However, the corresponding system does not have a strong positivity property anymore. Unconditional convergence starting from the origin does not hold in that case. Still, Theorems 4.8 and 4.11 show quadratic convergence in a neighborhood of the solution. A difficulty is that the characteristic system having the radius of convergence as a solution (Definition 3.7) may also have other solutions with positive coordinates. Bell, Burris and Yeats [2] have shown that the solution corresponding to the radius of convergence is the unique one where the largest eigenvalue of the Jacobian matrix is equal to 1. This result extends to our more general class of irreducible systems (Theorem 3.8).
Nonreal Dominant Singularities
While the exponential growth of the coefficients of a generating function is governed by its radius of convergence , subexponential terms in their asymptotic behavior are governed by the local behavior of the generating function at all the singularities on the circle , called dominant singularities. By a classical result of Pringsheim [19, Thm. IV.6], as generating functions have nonnegative coefficients, the point itself is a singularity. There might be other ones on the circle . A simple example (see Section 6) is
[TABLE]
a power series with positive integers for coefficients, whose singularities lie at , with . Generating functions coming from combinatorial systems behave more nicely than this example: the arguments of their dominant singularities are rational multiples of (Theorems 6.3 and 6.4) and moreover, these rational numbers are induced by periodicity properties that can be computed from the combinatorial equations. Results of that type are already known for -rational power series [6]; for a class of nonrecursive combinatorial generating functions [18, Prop. 4.6]; for polynomial combinatorial systems [4, Lemma 26]. This is the object of Section 6 and the resulting algorithm for constructible generating functions is given in Section 7.
Singular Behavior
Once the dominant singularities have been located, the next step (Section 8) is to compute the expansions of the generating functions in the neighborhood of these singularities. In the case of positive systems of polynomial or entire functions, Banderier and Drmota have classified the possible singular behaviors at finite singularities [1]. In particular, they proved that the exponents that can occur in the expansions are rational numbers whose denominator is a power of 2.
For combinatorial systems that may involve Cyc and Set, the possibilities are much more diverse. For instance, the generating function of the species of fragmented permutations [19, Ex. VIII.7] has dominant singularity at 1 and behaves like as . The species (that could be called ‘bicycle park’) also has a generating function with dominant singularity at 1 and behaves like as . Theorem 8.5 extends the classification of Banderier and Drmota to the general situation. There is a gap property that was already known in the nonrecursive case [18, Prop. 4.8]: roughly speaking, as , either the generating function grows at least as fast as in this last example, or it has an algebraic-logarithmic behavior, i.e., either it has an infinite limit and is equivalent to
[TABLE]
or it behaves like
[TABLE]
otherwise. In this last situation, as in the classification of Banderier and Drmota, the exponent is a power of 2. In the first one, less regular exponents occur. Theorem 8.5 shows that the generating function of the species222We write for the species .
[TABLE]
has a dominant singularity at , and as , its behavior is
[TABLE]
Another difference with the case of systems of polynomial or entire functions is that the exponents at nonreal dominant singularities can be nonreal as well (Theorem 8.6).
Despite these seemingly intricate exponents, the proofs of Theorems 8.5 and 8.6 give rise to algorithms that compute the singular expansions of constructible generating functions, provided they are of algebraic-logarithmic type. No new equality or inequality involving these constants need to be decided; all the delicate decisions that might possibly require an oracle have been made when computing the radius of convergence.
Asymptotic Expansion of the Coefficients
Expansions of algebraic-logarithmic type (1) are precisely those to which the transfer theorems of Flajolet and Odlyzko [17] apply. The end result of this article (Theorem 9.2) is that asymptotic expansions at arbitrary precision can be computed for the coefficients of all constructible generating functions that have algebraic-logarithmic dominant singularities. The only hypothesis is that Schanuel’s conjecture holds, or that the system contains neither Set nor Cyc. This article is concluded by Section 10, where further effectivity questions in this area are discussed.
2 Well-Founded Combinatorial Systems and their Characterizations
The symbolic method of Flajolet and Sedgewick [19] relies on systems of combinatorial equations, called specifications, between what they call classes of combinatorial structures. Algorithms to determine that such a specification does define a class of structures are known [48, 38]. This section presents two simpler characterizations (Theorem 2.8 with Algorithm 1, and Theorem 2.15). The language of the symbolic method is first recalled in Sections 2.1 and 2.2. It is sufficient to give a first characterization in Section 2.3. Then, the combinatorial derivative of species theory is recalled in Section 2.4 in order to present a second characterization in Section 2.5 with its proof that extends Joyal’s implicit species theorem [29]. This second characterization is the one that translates more easily at the level of analytic properties; it provides us with a useful tool for many proofs. Readers familiar with analytic combinatorics but not with species theory should be able to read this section with an intuitive understanding of species. For more information, the book of Bergeron, Labelle and Leroux [5] is our reference. A short introduction very close to our needs is given in our previous work [38, §1.1].
2.1 Well-founded Combinatorial Systems
The definition of specification used by Flajolet and Sedgewick [19, Def. I.7] uses the words ‘class’, ‘term’, ‘construction’. In the language of species theory [29, 5], the situation becomes simpler as these words all correspond to (possibly multisort) species. The combinatorial class described by the specification is then the species solution of a system of equations on species.
Definition 2.1**.**
A combinatorial system for an –tuple of species is a vector of equations,
[TABLE]
where each denotes a multisort species. (See Table 1 for simple examples of possible .)
Notation**.**
Boldfaced characters are used for vectors, matrices, or tuples of species; for example, a multisort species is written , where stands for the tuple ; a vector of multisort species is consistently written .
When several sorts of atoms are used, they will usually be denoted . Thus, the most general form of combinatorial system with several sorts of atoms considerered here is written .
Gantmacher’s notation is used to denote the -tuple . For instance, the species can also be written if its dimensions need to be made explicit.
The following general definition of well-founded combinatorial systems captures the families of systems considered in earlier literature. For a version of this definition directly at the level of generating functions, see Section 3.1.
Definition 2.2**.**
The combinatorial system is well founded when the recurrence
[TABLE]
is well defined for all and defines a sequence that is convergent. The limit of this sequence is called the combinatorial solution of the system.
In terms of combinatorial structures, ‘well defined’ means that the composition is possible, i.e., for each and each size , there are finitely many combinatorial structures of size in . Convergence means that there is a limit class such that for all , is isomorphic to up to the structures of size , for large enough . Alternatively, the corresponding sequence of generating functions converges as a sequence of power series for total degree.
Equivalently, in the context of species, ‘well defined’ means that for each coordinate , , either the species is polynomial in or . The class being isomorphic to up to the structures of size translates into the species having contact of order with , see [29, 5, 38].
Zero-free Systems and Comparison with [38]
It is natural to focus on combinatorial systems such that the limit of the iteration in Definition 2.2 has no zero coordinate.
Definition 2.3**.**
A combinatorial system is called zero-free if there exists such that for each coordinate , , , with defined by Eq. 3.
Our definitions of well-founded systems in [38] included that condition. Here, we remove this constraint from our definitions, which leads to a slightly larger class of systems that satisfy the conditions of Theorem 2.8 and Theorem 2.15 (e.g., the equation does).
All the results of our earlier work that are used in the present article, for this new definition, have been checked to hold with these extended conditions. Moreover, these coordinates can be detected algorithmically (for instance, using Algorithm 1 below) and removed from the system. This will be used as hypothesis when needed or when it lightens the presentation, knowing that any well-founded system can be split into two sub-systems, one being zero-free and the other one defining only zero-coordinates (see Corollary 2.10).
2.2 Constructive Systems in the Symbolic Method of Flajolet & Sedgewick
In order to distinguish the combinatorial systems used in the symbolic method from the general situation, we call specifications those that are considered by Flajolet and Sedgewick [19, Def. I.8].
Definition 2.4**.**
A specification is a combinatorial system of the form where the multisort species are obtained by combinations of (see Table 1).
In species language, is the species of singletons; is the characteristic species of the empty set; denotes the sum (or disjoint union); the product correspond to the cartesian product; sets, sequences and cycles are easily defined [5, 38]. This is summarized in Table 1, that also recalls how to compose the ‘construction’ species with the species 0 (that corresponds to the empty class) and the species 1 (that corresponds to the class containing only ) and how to obtain the derivatives discussed in the following section. For simplicity, we omit here a discussion on cardinality constraints that can be added to the sequences, sets and cycles. They can also be handled, see the discussion in the conclusion.
Definition 2.5**.**
A well-founded zero-free specification is called a constructive system. The species solution are called constructible.
The name of these systems comes from the fact that they are exactly the systems that define constructible combinatorial classes in the work of Flajolet and Sedgewick [19, Def. I.7].
In our algorithms, it is convenient to make use of systems of a special form similar to the Chomsky normal form for context-free grammars, see also [48, p. 21], [20].
Definition 2.6**.**
A specification is in normal form when its equations are of one of the forms:
[TABLE]
where each and is one of the constructions in (see Table 1).
Any specification can easily be transformed into a system in normal form. Thus all the algorithms presented in this article can take as input a specification in normal form. Most of them also require the zero-free constraint that is equally easy to satisfy.
Example 1. The following specification will be a running example.333A Maple session dedicated to this example can be found in the arXiv along with this article. It defines unordered forests of trees with red, green and blue nodes and various other constraints. Its normal form is given on the right. For this example, it is possible to obtain closed forms for the generating functions. While these are not used by our algorithms, they allow to check the results easily.
[TABLE]
[TABLE]
2.3 Characterization of Well-founded Systems by Leading Terms
A first characterization of the systems that are well founded can be given in terms of leading terms of the iterates. In Section 2.5, a second characterization, based on the Jacobian matrix of , is given. It will be used in the proofs of Section 3.
Definition 2.7**.**
Let be a species. The valuation of is the size of smallest -structure (the valuation of the species 0 is infinite). The leading term of is the species restricted to -structures of size (and the leading term of the species 0 is itself). The leading terms of a vector of species is the vector of the leading terms of . It is denoted .
Theorem 2.8** (Characterization by Leading Terms).**
Let be a vector of species such that the iteration (3) is well defined for . Then the system is well founded if and only if the leading terms of and are equal. In this situation, the combinatorial solution of the system satisfies .
The proof is given in Appendix A.
Thus in at most iterations, one can detect whether a system is well founded or not. This characterization at the level of species yields Algorithm 1 for specifications (where there is only one sort of atom ). This algorithm is close in spirit with early work on these questions [48, 18] for zero-free specifications and results in a much simpler algorithm than our previous one [38]. The algorithm given by Zimmermann [48] (Algorithm B, p. 33) and summarized in [18] first computes the valuations (this corresponds to computing the second coordinate in Algorithm 1). Next, a cycle detection is performed on the dependency graph of the system (for a restricted notion of dependency). As explained in the next section, this corresponds to checking the nilpotence of the Jacobian matrix of at different iterates, which will appear in our second characterization. Our new characterization by leading terms allows us to bound the number of iterations in the computation of valuations and results in a more straightforward algorithm. Also, it does not require removing the zero coordinates beforehand and can be used to detect them (Corollary 2.10 below).
Theorem 2.9**.**
Given as input a specification in normal form, Algorithm 1 returns FAIL if and only if the system is not well founded. Otherwise, it returns a vector of pairs that, for each coordinate, indicate the minimal size of the structures defined by that coordinate (its valuation) and the number of structures of that size.
Proof.
At any step , either Algorithm 1 computes the vector of pairs (count, valuation) corresponding to the leading terms of or the system is not well founded. This is proved by induction using Lemma A.1: if is well defined, ; the rest follows from basic properties of the constructions of Definition 2.4. It is straightforward for most of the constructions:
[TABLE]
an slightly more involved for products:
[TABLE]
and sums:
[TABLE]
where . If the algorithm does not fail, a vector corresponding to the leading terms of , for , is returned. By Theorem 2.8, these are the leading terms of the species solution of the system. ∎
Corollary 2.10**.**
Given a well-founded specification , Algorithm 1 detects its zero-coordinates as those for which it returns . If is such that are the zero-coordinates, then the system defines the same solutions as the nonzero ones in the initial system.
Proof.
The zero coordinates are those that have an infinite valuation. ∎
Example 2. On the system (5) of Definition 2.6, the first iterations of the main loop in Algorithm 1 produce the following vectors, with each pair (coefficient, exponent) written as :
[TABLE]
The last two vectors are equal and therefore, the system of 16 equations is found to be well founded in 6 iterations.
Example 3. If the last two equations of the system defined by Eq. 5 are replaced with another species of trees, , , the system is not well founded anymore. Indeed, the second iteration produces
[TABLE]
which means that (in blue), which will break the condition in line 1 for equation at the next iteration. However, if this equation is also rewritten as , , the condition in the loop now holds. But in this case, the last two vectors after iterations are not equal:
[TABLE]
showing that the new system is not well founded.
Remark 2.11*.*
The number of iterations in Algorithm 1 is bounded by the number of equations in the input system. Using a system in normal form allows for a simple presentation of the algorithm but it can increase significantly the number of iterations, in comparison with the original system (as between the ones of Eq. 4 and Eq. 5). It is always possible to work with a system that is not in normal form (as in [48]), by making explicit how to compute the leading terms of any defined as in Definition 2.4.
2.4 Combinatorial Derivative
The derivative of species theory is crucial in this work. In terms of classes of structures, the derivative of a combinatorial class (with respect to ) is another combinatorial class containing all the possible structures obtained by erasing a in a structure of and replacing it by a distinguished . Then, the product of a derivative with a class can be interpreted as replacing a distinguished in the derivative by a structure in . While Flajolet and Sedgewick do not define a derivative, their closely related pointing construction can be stated in terms of derivatives: .
Partial derivatives are obtained when considering constructions or terms with several sorts of atoms. For example, erasing an atom of the sort in all possible ways in the class taken from Eq. 2 produces the class . Structures of are structures of where a has been replaced by a structure of the derivative of with respect to the sort .
Example 4. Let . The class can be interpreted as the set of sequences of atoms of sort or where one of them as been erased and its position remembered. Similarly, the class is viewed as the set of sequences of or where one is distinguished, and the class as the set of sequences of or and only one
Matrices and vectors of combinatorial classes are defined as usual and the classical notion of Jacobian matrix extends to combinatorial systems. The product of a matrix by a matrix or a vector is also given by the usual rules (with disjoint unions and cartesian products as sums and products), which allows to define nilpotent combinatorial matrices as usual. This notion from species theory is the key element to decide whether a system is well founded.
Example 5. The Jacobian matrix for the combinatorial system of Eq. 4 is
[TABLE]
Definition 2.12**.**
The dependency graph [19, p. 33] of a combinatorial system is the directed graph whose vertices are the s and such that there is an edge from to when appears explicitly on the right-hand side of the equation for . Equivalently, the dependency matrix of this graph is given by the location of the nonzero entries of the Jacobian matrix .
Example 6. For the system of Eq. 4, the dependency graph of is
\mathcal{F}$$\mathcal{T}_{r}$$\mathcal{T}_{b}$$\mathcal{T}_{g}$$\mathcal{G}$$\mathcal{R}$$\mathcal{B}
2.5 Characterization of Well-founded Systems by Jacobian Matrices
The following characterization and its proof use the language of species theory, drawing in part from our previous work [38]. The starting point is the following theorem by Joyal [29].
Theorem 2.13** (Joyal’s implicit species theorem).**
If and is nilpotent, then the system is well founded and has a solution with , that is unique up to isomorphism.
A partial converse of this result (for zero-free systems) holds, showing that its condition is essentially tight.
Proposition 2.14**.**
[38, Thm. 3.6]** If is well founded with and the limit of Iteration (3) has no zero coordinate, then is nilpotent.
These results generalize to the case when is not necessarily . The analogue of Joyal’s theorem is the following.
Theorem 2.15** (Implicit Species Theorem).**
Let be a vector of species. If
letting , for each , is well defined; 2. 2.
the matrix is nilpotent,
then the system is well founded and has a solution with , which is unique up to isomorphism.
Proof.
The sequence satisfies by induction: this is true for and for larger , the inclusion is preserved by composition with . Next, the identity
[TABLE]
shows that any -structure that is not a -structure is an -assembly of -structures, where at least one lies in the difference . Iterating and using the inclusion for all proves
[TABLE]
Since is nilpotent of order , the right-hand side is and .
Let us now consider the shifted system with ; the subtraction is well defined since
[TABLE]
This shifted system satisfies the conditions of Joyal’s Implicit Species Theorem 2.13, that are and nilpotent. Therefore, it admits a solution such that , unique up to isomorphism. Now, is a solution of such that . Any other solution is such that is isomorphic to , which concludes the proof. ∎
Example 7. The system defines binary trees with leaves that are nonempty sequences. Theorem 2.15 shows that this is well founded since the first condition naturally holds and defines and the Jacobian matrix is nilpotent at and .
On the contrary, the system uses possibly empty sequences for the leaves, which can create infinitely many trees of size 0. The system is not well founded, since is now and the Jacobian matrix is not nilpotent at and .
Remark 2.16*.*
It is not sufficient that to ensure that the system is well founded. For instance, the system defined by has , . However, this system produces infinitely many -structures of size . The second condition of the theorem does not hold: the Jacobian matrix has a block decomposition with 1 in the top-left corner. (This is also detected by Algorithm 1 since and ). This example shows also that in the second point of Theorem 2.15, one cannot replace by in the Jacobian matrix.
Actually, the conditions of Theorem 2.15 are not far from optimal, as shown by the following.
Proposition 2.17**.**
If the system is zero-free and well founded, then the conditions of Theorem 2.15 hold.
Proof.
By induction, for all , . Part 1. of Theorem 2.15 follows from , which is well defined.
Next, Theorem 5.5 of [38] shows that, under the assumptions of this proposition, the companion system is well founded; it has a solution such that the following properties hold: ; the limit is ; is polynomial in . By [38, Prop. 4.5], this implies that is nilpotent. This implies that is nilpotent too. Finally, by [38, Thm. 5.7], . ∎
Example 8. Consider again the combinatorial system of Eq. 4. is the vector whose coordinates indexed by are . Here is the Jacobian matrix . It is nilpotent, showing that the system is well founded.
[TABLE]
Example 9. If the last two equations of the system of Eq. 5 are replaced with another species of trees, , , the system is no longer well founded. Indeed, breaks the first condition of Theorem 2.13 (in the equation ). However, if the equation is also rewritten as , , the condition in the loop now holds. But in this case, the resulting Jacobian matrix is not nilpotent.
Example 10. As regards 0-coordinates, it is straightforward to check that the conditions of Theorem 2.15 hold for the system which is well founded and defines 0. Still, the condition ‘zero-free’ in Proposition 2.17 is necessary: the system is also well founded as the sequence defined by Eq. 3 converges to , even though the second condition of Theorem 2.15 does not hold.
2.6 Decomposition into Strongly Connected Components
Some characteristics of a well-founded system can be found on its dependency graph:
Definition 2.18**.**
A well-founded system is recursive when its dependency graph has a cycle, or equivalently when the matrix is not nilpotent.
Definition 2.19**.**
A combinatorial system is irreducible when it is well founded, zero-free, recursive and its dependency graph is strongly connected (every vertex can be reached from any of the other ones).
Note that these definitions would also apply to combinatorial systems that are not well founded but this situation is not considered in the rest of this article.
Example 11. The combinatorial system of Eq. 4 is recursive : the corresponding Jacobian matrix is not nilpotent and, indeed, the corresponding dependency graph has many cycles. It is not irreducible since four strongly connected components can be identified : , and are irreducible while is nonrecursive.
The analysis of well-founded systems is made simpler by focusing on irreducible systems first, and reducing the general case to the irreducible and the reducible cases separately, through the classical decomposition of a directed graph into a directed acyclic graph (DAG) of strongly connected components. Since the components form a DAG, without loss of generality, the species can be relabeled in such a way that the successive parts compose a block-triangular system of the form
[TABLE]
These blocks will be analyzed from the first one to the last, so that any species is analyzed before it occurs in the definitions of other species.
Definition 2.20**.**
By analogy with graphs, the blocks in the system (6) are called components. The ones consisting of a single equation with are the nonrecursive components. The other ones are called the irreducible components of the system. Components of constructive systems are called constructive components. The definition of normal form also extends to any component by including in the set from where can be taken in Definition 2.6.
Example 12. The dependency graph of the system given in Definition 2.12 decomposes into four strongly connected components identified in the left of the following picture. The corresponding DAG is given on the right. Both the subsystems and are irreducible (and thus, irreducible components). Note that the system is an irreducible component defining while the subsystem is not irreducible. The system defining is a strongly connected component in the graph, but being nonrecursive, it does not form an irreducible component.
\mathcal{F}$$\mathcal{T}_{r}$$\mathcal{T}_{b}$$\mathcal{T}_{g}$$\mathcal{G}$$\mathcal{R}$$\mathcal{B}
\mathcal{F}$$\mathcal{T}_{r}\quad\mathcal{T}_{b}$$\mathcal{T}_{g}$$\mathcal{G}$$\mathcal{R}$$\mathcal{B}
3 Well-Founded Systems of Generating Function Equations
A combinatorial system translates into a system over generating functions. In the case of constructive systems, this is direct using the dictionary presented in Table 1. A few examples help understand the type of systems considered in this section.
Example 13. The system of combinatorial equations from Definition 2.6 translates into the system
[TABLE]
where the equations here are sorted by connected components of the system, following the graph in Definition 2.12.
Example 14. Functional graphs [19, II.5.2] decompose as sets of connected components, each consisting of a cycle of trees. This leads to the system of generating function equations
[TABLE]
More generally, we give a definition of well-founded analytic systems that reflects the combinatorial conditions of Theorem 2.15 and extends similar definitions that were considered in the literature [9, 10, 4, 2, 38]. This section then develops the mathematical tools that are used in the development and proof of the algorithms operating on these well-founded systems in the rest of this article.
A difficulty is that for any given , such a system of generating function equations may have many solutions and the only solution of interest is the one that corresponds to the value of the generating function, which is the one we denote by . Even restricting to solutions with nonnegative coordinates at a real may still result in several possible choices. Still, Proposition 3.11 shows that when such solutions exist, then is at most the radius of convergence of the generating function. This is the basis for a dichotomy that reduces the computability of the radius of convergence to the existence of nonnegative real solutions of systems of equations in Section 3.4.
The computation of the radius of convergence — indeed, much of the analysis of these systems — reduces to the special case of irreducible systems. When these are polynomial systems, their analysis is summarized by Flajolet and Sedgewick [19, Ch. VII]; a more general situation is handled by Drmota [9, 10]. A key insight of Bell, Burris and Yeats [2] is that one can distinguish which of the nonnegative real solutions of an irreducible system corresponds to the value of the generating function by considering the dominant eigenvalue of the Jacobian matrix of with respect to . Their hypotheses are slightly more restrictive than those used here, but most of their proofs extend to our more general situation (Theorem 3.8).
Notation** (for power series).**
This section focuses on generating functions in one variable , which explains why the boldfaced is not used. In order to insist on the values taken by the generating functions at real points, sometimes is used instead of . The ring of formal power series with real coefficients is denoted and similarly when the coefficients are complex. Their field of fractions are denoted and and their elements are called formal Laurent series. The subrings and denote the convergent power series, while the subfields and denote power series that are convergent in a punctured neighborhood of 0. When is a positive integer, changing into gives the corresponding rings and fields of Puiseux series. The elements of and converge in a slit neighborhood of 0. Finally, for , changing into gives the rings and fields used to express expansions in the neighborhood of .
3.1 Well-founded Systems
3.1.1 Definition
Systems of generating function equations that are well founded are defined in analogy with Theorem 2.15.
Definition 3.1**.**
Let be analytic in a neighborhood of . The system is called well founded when
each coordinate of has nonnegative real Taylor coefficients at ; 2. 2.
letting , for , belongs to the domain of convergence of the Taylor series of at , allowing one to define ; 3. 3.
the matrix is nilpotent at .
In the common case when , condition 2 is automatically satisfied and .
Definition 3.2**.**
The well-founded system is called recursive when is not nilpotent and irreducible when moreover it is well-founded and zero-free and its dependency graph (defined as in Definition 2.12) is strongly connected. This is equivalent to the irreducibility of the matrix in the sense used in the Perron-Frobenius theory.
Definition 3.1 is analogous to Definition 2.19 in the combinatorial situation. Similarly, nonrecursive and irreducible components are defined in the same way as in Definition 2.20. Also, the system is called constructive when it is the translation of a constructive system in the sense of Definition 2.5 by the rules of Table 1.
3.1.2 Solutions
This definition of well-founded systems is motivated by the following analytic analogue of Theorem 2.15.
Proposition 3.3**.**
If is a well-founded system, then it admits a unique solution analytic in a neighborhood of [math] with given in Definition 3.1. This solution has nonnegative Taylor coefficients at [math].
Definition 3.4**.**
The solution is called the generating function solution of the system. When the system is constructive, it is called a constructible generating function.
Proof.
If , then and the existence of an analytic solution with is granted by the implicit function theorem at (see, e.g., [7, Chap. IV]). Nonnegativity of the Taylor coefficients of comes from the fixed-point iteration over power series, starting with .
In general, the implicit function theorem at shows that the system admits an analytic solution with and again nonnegative Taylor coefficients. If , the function is then a solution of with nonnegative Taylor coefficients and .
It remains to show that . The proof follows the steps of that of Theorem 2.15. All inequalities between vectors mean that they hold for all coordinates.
The sequence satisfies . For , this follows from the nonnegativity of the Taylor coefficients of and for higher by the fact that this nonnegativity makes each coordinate increasing. Next, by the lemma below,
[TABLE]
Iterating and using gives
[TABLE]
Lemma 3.5**.**
Let be as in Definition 3.1 and let and belong to the domain of convergence of with . Then
[TABLE]
Proof.
Since has nonnegative Taylor coefficients at , then so do all its derivatives, which implies in particular that is convex on the intersection of its domain of convergence with . Both inequalities are consequences of that convexity. ∎
3.1.3 Comparison with Previous Definitions
Drmota [9, 10] requires and thus , but also that and the solution has no zero-coordinates. Bell, Burris, Yeats [2, 3] have two definitions: elementary systems of power series equations whose dependency graph is not necessarily strongly connected and where it is only required that and ; and well-conditioned systems with Drmota’s conditions and where moreover should not have zero-coordinates, but with no condition on the solution. This is summarized in Table 2. Despite these variations, some of their results can be used in our context, by the following reductions.
Zero coordinates at and zero-free systems.
As stated earlier, these coordinates can be detected algorithmically (using Algorithm 1) and removed from the system.
Nonzero nilpotent Jacobian matrix.
The case when is nilpotent but not reduces to the case when it is [3, Section 4.4]: writing leads to a new well-founded system for the same ,
[TABLE]
where that property holds. Note that if is the constant , then and then , which shows that this property is preserved by the reduction.
Nonzero value at .
Similarly, if , a new system can be formed for :
[TABLE]
Recall that is obtained as .
This system satisfies , which follows from . It is well founded: the Taylor coefficients of are nonnegative, which can be seen by expansion of ; is nilpotent.
Nonzero value at .
This case is a generalization of the previous one. It contains in particular combinatorial systems that put weights on leaves rather than internal nodes, such as binary trees defined by .
A new system can be formed for defined by
[TABLE]
As before, is obtained as and from there as the solution of the linear system
[TABLE]
Then, satisfies
[TABLE]
By Taylor expansion, . Moreover, this new system is well founded: as before, the Taylor coefficients of are nonnegative; the Jacobian matrix is
[TABLE]
whose value at is , which is nilpotent.
3.1.4 Irreducible Systems
The following lemma ensures that results on irreducible systems apply to irreducible components of well-founded systems.
Lemma 3.6**.**
Let be a well-founded system that is not irreducible and splits into two sub-systems
[TABLE]
where the second one is an irreducible component. Let be the solution of the first component given by Proposition 3.3. Then the system
[TABLE]
is an irreducible system.
Proof.
By definition, is an irreducible system. It is also well founded: (1) has nonnegative Taylor coefficients at by composition of with the generating functions ; (2) recall from Definition 3.1,that and for , . The corresponding iteration for at is defined by and for , . By positivity, for , . The properties required by Definition 3.1 can now be checked. Indeed, as belongs to the domain of convergence of the Taylor series of at , so does for at , for ; (3) the bottom-right block of size in is precisely , while the bottom-left block is zero. Since is nilpotent, so is its bottom-right block. Thus
[TABLE]
3.2 Characteristic System and Dominant Eigenvalues
We consider a well-founded system and denote by the value of the generating function at 0, which is given by Proposition 3.3. The implicit function theorem implies that the solution is analytic on the positive real axis as long as belongs to the domain of convergence of and is not a solution of the characteristic system in the following sense.
Definition 3.7**.**
The characteristic system of the system is the extended system
[TABLE]
The characteristic system may possess several solutions with positive coordinates. An example due to Burris is
[TABLE]
whose characteristic system has both and as solutions [19, Ex.VII.28] and only the first one, with larger value of , is of interest. Thus, even if an algorithm converges to a solution of the characteristic system with positive coordinates, more work is required in order to prove that this solution corresponds to a dominant singularity of the generating function solutions.
An idea due to Bell, Burris and Yeats [2] solves this problem for irreducible systems. We build on their study to design our algorithms. For a well-founded system, at any point inside the domain of convergence of , the Jacobian matrix has nonnegative entries. By the Perron-Frobenius theorem, this implies that its spectral radius is also its largest real eigenvalue, that is called its dominant eigenvalue. In particular, the dominant eigenvalue gives a criterion to check that a solution of the system is the value of the generating function at for irreducible systems. This will be a building block for our algorithm for the general case in the next section.
Theorem 3.8**.**
Let be an irreducible system. Let belong to the domain of convergence of or its boundary and satisfy . Let be the radius of convergence of the generating function solution . Then, and
if is linear in the variable then and ; 2. 2.
in all cases
- (a)
* and ;* 2. (b)
.
Notation**.**
For , we write for .
An application of Part 2 of this theorem is to filter out solutions of the characteristic system (8) that do not correspond to a dominant singularity. This was the motivation for the work of Bell, Burris and Yeats, who use slightly stronger hypotheses. An adaptation of their proof to Theorem 3.8 is given in Appendix B.
Corollary 3.9**.**
Let be an irreducible system. Let be the radius of convergence of the generating function . Assume that is finite. Then either is on the boundary of the domain of convergence of or .
Proof.
When is finite, . By continuity, is either in the domain of convergence of or on its boundary. Item 2a of the theorem implies that . If moreover lies inside the domain of convergence of , then by the implicit function theorem, the singularity occurs at a solution of the characteristic system (8), where . ∎
Corollary 3.10**.**
Theorem 3.8* and Corollary 3.9 hold more generally if the system is an irreducible component of a well-founded system.*
Proof.
By Lemma 3.6. ∎
3.3 Recognizing Points Inside the Domain of Convergence
Using the decomposition of general systems into irreducible components, it is possible to ensure that a point is inside the domain of convergence of the generating function solution of a well-founded system that is not necessarily irreducible, by the following.
Proposition 3.11**.**
Let be a zero-free well-founded system, the domain of convergence of and the radius of convergence of the generating function solution . If there exists such that , then , exists and .
Note however that, except in the linear case, one cannot conclude from the hypotheses that . For instance, the equation for binary trees has the positive solution at . This is not equal to , which is 0. Since this system is irreducible, by Theorem 3.8, detecting that is achieved by checking that . In this example, and .
In general, the inequality is similar to the classical Kleene fixed point theorem. The main use of this proposition however is that it ensures that the system does not have nonnegative solutions beyond the radius of convergence of the generating function solution .
As the proof of Proposition 3.11 uses several of the lemmas proved in Appendix B, it is postponed to Appendix C.
3.4 Radius of Convergence of Constructible Generating Functions
The results on general well-founded systems are now applied to the constructive systems of Flajolet and Sedgewick (Definition 2.5).
3.4.1 Finiteness
It is easy to detect when the radius of convergence is infinite.
Proposition 3.12**.**
The radius of convergence of the constructible generating function solution to is infinite if and only if the system is not recursive and does not contain Seq or Cyc.
Proof.
Let be the generating function solution. When the system is recursive and irreducible, is not nilpotent while its value at is. This implies that is nonzero at for . Since
[TABLE]
it follows by evaluation at that is positive. Then is a nonconstant power series with nonnegative coefficients. If its radius of convergence were infinite, the series would tend to infinity as . Then, so would the dominant eigenvalue , which is bounded by 1, a contradiction.
When the system is recursive and not irreducible, it has an irreducible component to which the previous reasoning applies.
Finally, if the system is not recursive, then is the limit of the iteration of power series with . By Lemma 3.5,
[TABLE]
Since the system is nonrecursive, the matrix is nilpotent and the right-hand side is , which implies that , the th iterate of at . Thus each coordinate of is an explicit combination of (the last two correspond to Seq and Cyc) with additions, multiplications, and the constraints that and are applied to nonzero power series (since the system is zero free) of positive valuation. By induction such a function is entire if and only if neither nor occur in its expression. (See also [18, Prop. 4.2].) ∎
3.4.2 Decidability
Proposition 3.11 gives a basis for an algorithmic characterization of the radius of convergence of generating functions. When the system has a polynomial , then its domain of convergence is all of and the question of the existence of a solution in is decidable by the Tarski-Seidenberg theorem. For the constructive systems of Definition 2.5 that also involve the function , testing that a point belongs to the domain of convergence is also given by a polynomial inequality and thus also decidable.
In general, constructive systems have generating function equations that also involve logarithms and exponentials (see Table 1). Provided the system is written in normal form (Definition 2.6), the hypothesis is expressed by inequalities , for all and for the that are arguments to the functions , . These inequalities, together with the existence of a solution of the system form a formula in the theory of the real numbers with the exponential. This theory was proved decidable by Macintyre and Wilkie [32], under the assumption that the following classical conjecture in number theory holds.
Conjecture 3.13** (Schanuel).**
If are real numbers that are linearly independent over , then the field has transcendence degree at least over .
A direct consequence of these strong results from mathematical logic is the following.
Theorem 3.14**.**
If the radius of convergence of a constructible generating function is finite, then it is a computable number, provided either that its defining system does not contain Cyc or Set, or that Schanuel’s conjecture holds.
Proof.
Recall that by Definition 2.5, constructive systems are zero-free. A simple process sufficient for the computability is a search by dichotomy, using Proposition 3.11 that reduces the decision to the existence of nonnegative solutions of the system in the domain of convergence of . ∎
3.4.3 Constant Oracle
The existence problems that must be decided in the use of Proposition 3.11 are of a specific form and might well be decidable in a more direct way than by the result of Macintyre and Wilkie. For this reason, we state our algorithms in terms of a more restricted oracle.
Definition 3.15**.**
The Constant Oracle is an oracle that takes as input functions with s of one of the forms () and the inequalities for all , and for the coordinates involved in . It is able to decide:
given and the inequalities, whether the system has a solution satisfying the inequalities; 2. 2.
if so, whether it has a solution with ; 3. 3.
given two such systems and with solutions that satisfy their inequalities, which of holds.
A weaker oracle was used in [18], where only nonrecursive equations were considered. The following is a reformulation of Theorem 3.14.
Proposition 3.16**.**
If the radius of convergence of a constructible generating function is finite, then it can be computed with the help of the Constant Oracle.
Again, in the case of context-free systems, i.e., when the combinatorial system contains neither Cyc nor Set, all the numbers involved are real algebraic numbers and the Constant Oracle is computable. Efficient algorithms and implementations are available, such as Safey El Din’s RagLIB package.444Available at https://www-polsys.lip6.fr/~safey/RAGLib/.
In practice, in the more general situation, one can use numerical approximations and assume that if two constants agree for their first few hundred digits, then they coincide. This will very likely be satisfactory unless the system has been designed for the purpose of countering this method.
4 Newton Iteration for Well-Founded Systems
The classical Newton iteration for solving equations or systems of equations numerically is particularly powerful when applied to combinatorial systems, whose positivity allows for more control over the convergence than in a general situation. For inside the disk of convergence of the generating functions, when started at , Newton’s iteration always converges to the values of the generating functions at [38]. A converse holds: if Newton’s iteration started at converges, then , where is the radius of convergence (Section 4.1). Positive systems with polynomial entries have been considered in detail by Etessami, Yannakakis, Esparza, Kiefer, Luttenberger [15, 13, 14, 42]. Their algorithmic approach is similar to ours and faces the same difficulties: if the iteration is stopped prematurely, it could be that while the divergence of the iteration has not manifested itself yet. For this situation, we give a Newton-Kantorovich-type a posteriori bound on the output of the iteration that can be used to certify the convergence (Section 4.2). The positivity of combinatorial systems allows to derive such a bound from only one evaluation, instead of the required bound on a whole disk in the general situation. An analysis of the speed of convergence in the spirit of [13] is also possible (Section 4.3).
With these results, one can compute radii of convergence by dichotomy, but this is not as fast as a direct use of Newton’s iteration. For irreducible systems, the radius of convergence is given as a solution of a characteristic system (Definition 3.7). Unfortunately, this is not a positive system and the unconditional convergence of Newton’s iteration does not hold for these systems. Still, there is a neighborhood of the solution where the convergence is quadratic (Section 4.4). The same property holds for systems expressing that one of the coordinates is 1, that occur when looking for singularities of Seq and Cyc (Section 4.5).
4.1 Convergence and Divergence of Newton’s Iteration
The following strengthens our previous result [38, Thm. 9.13] when is real. It shows that for the values computed by Newton’s iteration, one has more precise information than is provided in general by Proposition 3.11.
Theorem 4.1**.**
Let be a well-founded system and let be the radius of convergence of the generating function of the solution. For , let denote the sequence defined by and
[TABLE]
whenever belongs to the domain of convergence of , and undefined otherwise. Then,
- —
If , then belongs to the domain of convergence of for all and is a strictly increasing sequence converging to ;
- —
conversely, if converges to a vector such that belongs to the domain of convergence of or its boundary, then and ;
- —
if , then is an increasing sequence for until a value such that does not belong to the domain of convergence of or such that .
Proof.
The first part is [38, Thm. 9.13], without the constraint of being zero-free on the system (which can be removed using Theorem 2.15). The second is a consequence of Proposition 3.11 and of the first part.
For the last part, we first show that as long as belongs to the domain of convergence of and , the sequence increases. Indeed, by positivity of , it follows that and by induction showing that the last vector in Eq. 10 has nonnegative coordinates. Next, since , the inverse matrix in the right-hand side of Eq. 10 has nonnegative entries, showing that the sequence is increasing.
Finally, assume by contradiction that , that the sequence is increasing, that for all , belongs to the domain of convergence of and that . If the increasing sequence were bounded, it would converge, but this is impossible by the second part of the theorem. Thus, one of the coordinates of tends to , implying that and depend on this coordinate. The monotonicity of makes this incompatible with the fact that for all . ∎
This theorem gives a practical criterion for detecting that a value is larger than the radius of convergence: the vector computed by the Newton iteration will ultimately have a decreasing coordinate or be outside the domain of convergence of or be such that . An explicit variant of the method deduced from this theorem is given in 2.
Remark 4.2*.*
This method is very similar to the approach taken by Etessami and Yannakakis [15] for positive polynomial systems, which it extends to the case of well-founded systems. Moreover, for systems of polynomial equations coming from recursive Markov chains, they use a ‘decomposed Newton’s method’, where the decomposition of the system into strongly connected components is used to solve one system per component by Newton iteration along the dag of the decomposition. The above theorem shows that, in the context of well-founded combinatorial systems, this is not necessary for convergence. Still, our implementation uses the decomposition, as this leads to smaller Jacobian matrices.
Remark 4.3*.*
A practical implementation of Newton’s iteration (10) uses floating-point numbers. This introduces further issues due to rounding. While the convergence of Newton’s iteration is known to be robust, inequalities like increasing coordinates may fail to hold close to the radius of convergence unless some care is taken. Thus, it is important to use directed rounding towards 0 or (these rounding modes are provided by the IEEE-754 standard) during the evaluation. In addition, it is necessary that the implementation of the evaluation of preserves the inequality and the positivity of . In that case, the monotonicity of the rounding operation implies that one can stop the iteration when a fixed point is reached on floating-point values. We do not consider these issues further in this article.
4.2 A Posteriori Bounds
Theorem 4.1 gives information on the limit of the sequence of Newton iterations. In practice, when the iteration is stopped after a finite number of steps as in 2, the theorem does not help determine whether or not, except when the algorithm returns ‘’.555This is the same problem as in [15, Cor. 7.5 and Footnote 6]. Still, the positivity of our systems is a strong property, that makes it easier to apply the Newton-Kantorovich theorem: in a more general situation, that theorem requires control of the system in a whole neighborhood of the iterate, while here it is only necessary to perform one evaluation at one point in order to certify that .
Theorem 4.4**.**
With the hypotheses and notation of Theorem 4.1, let and , where . If belongs to the domain of convergence of , and
[TABLE]
then and for all ,
[TABLE]
In other words, when the conditions hold, the current iterate is at distance less than twice the length of the Newton step from and moreover, from that point on, the convergence is bounded by that of a quadratically convergent sequence.
Proof.
The proof reduces to a variant of the Newton-Kantorovich theorem [8, §3.2]. In order to simplify notations, we write
[TABLE]
and observe that and .
By positivity of and all its derivatives for nonnegative arguments inside the domain of convergence of , for any such that , for each in ,
[TABLE]
For any two vectors in ,
[TABLE]
where is the Hessian matrix of the th coordinate of at . As is a matrix with nonnegative entries, by triangular inequality, the product on the right-hand side is a vector whose coordinates are bilinear forms in with coefficients that are upper bounded by those obtained by replacing by in its expression. The maximal value of the coordinates of this vector over all of norm 1 is reached at and thus the norm in Eq. 11 is an upper bound on the norm of the operator for .
In this situation, the hypotheses of the Newton-Kantorovich theorem [8, Thm. 88] hold. It follows that:
- –
there exists a unique with such that ;
- –
Newton’s iteration converges to it at the speed given in the last inequality;
- –
is invertible.
It remains to see that . Since is inside the domain of convergence of , so is . Also, since , it follows that . Then follows from the second part of Theorem 4.1. That is a consequence of the invertibility of . ∎
Example 15. Let be the combinatorial system defined by Eq. 4 in Definition 2.6 and its radius of convergence. For , after 13 steps of the iteration, the last step has length longer than . Still, the vector is positive, its entries are larger than the step and the inequality (11) becomes , which is sufficient to conclude that . Conversely, for , Newton’s iteration (10) at produces at least one decreasing coordinate, showing that by Theorem 4.1. Thus and this can be refined by dichotomy. A direct Newton iteration is discussed in Section 4.4.
4.3 Speed of Convergence
Even when and the convergence of Newton’s iteration becomes eventually quadratic as in Theorem 4.4, there is no clear control a priori of the number of steps needed before that quadratic regime takes place. Parts of the analysis of Esparza et al. [13] are now extended to our context where the equations are not necessarily polynomial. The basic result allowing for the generalization is the following.
Lemma 4.5**.**
Let be a well-founded system, and be such that is inside the domain of convergence of , then
[TABLE]
Proof.
Since has nonnegative Taylor coefficients at 0, so do all its derivatives and therefore so does its Taylor expansion at any point with nonnegative real coordinates inside its domain of convergence.
The case of dimension 1 is obtained by Taylor expansion. We write for . By Taylor expansion,
[TABLE]
Reorganizing terms concludes the case of dimension 1.
In higher dimension, it is sufficient to apply this reasoning to each . ∎
This lemma gives control over the speed of convergence of Newton iteration as follows.
Lemma 4.6**.**
Let , and be as in Theorem 4.1 and let be an eigenvector of with positive entries for the eigenvalue . If is such that then .
Proof.
The existence of is given by the Perron-Frobenius theorem. The result is obtained by combining the inequalities as follows
[TABLE]
In the third line, the inequality follows from Lemma 4.5 and the positivity of the inverse of the matrix, which is a consequence of . The next line comes from the induction hypothesis. The following one uses the fact that is an eigenvector or for . ∎
We can finally conclude on the speed of convergence of Newton’s iteration inside the domain of convergence of the generating function.
Theorem 4.7**.**
Let be a well-founded system and let be the radius of convergence of the generating function of the solution. Let be defined by the Newton iteration from Eq. 10.
Let be such that belongs to the domain of convergence of and let be an eigenvector of with positive entries for the eigenvalue . Then there exists such that for all ,
[TABLE]
Proof.
The second inequality comes from Theorem 4.1. Since both and have nonnegative entries, there exists such that . This value gives the result for and the conclusion then follows by induction from Lemma 4.6. ∎
It follows that a linear number of iterations is sufficient to obtain a linear number of digits of the values, however close is to and even for if belongs to the domain of convergence of . A more refined analysis along the lines of that of Esparza et al. [13] could be of interest.
4.4 Radius of Convergence of Irreducible Systems
The characteristic system of Section 3.2 has a set of solutions that contains , where is the radius of convergence of the generating function solution. The computation of plays an important role in random generation by Boltzmann sampling, where it is used in the singular Boltzmann samplers [12, § 7]. That system is not a positive system in the sense of Definition 3.1 and therefore Theorem 4.1 cannot be applied to guarantee convergence of Newton’s iteration. Still, denoting by the matrix , the Jacobian matrix of the characteristic system can be observed to have the block structure
[TABLE]
which can be used to perform a Newton iteration in the irreducible case, thanks to the following.
Theorem 4.8**.**
Let be an irreducible system, let be its generating function solution and its radius of convergence. If is inside the domain of convergence of , then there exists a neighborhood of such that Newton’s iteration for the characteristic system, i.e.,
[TABLE]
converges quadratically to when started in that neighborhood.
Note that if Newton’s iteration for the characteristic system converges to a real positive value , then computing can be used to certify that , by Theorem 3.8.
Proof.
It is classical that it is sufficient to prove that the Jacobian matrix in Eq. 12 is invertible at [8]. By Schur’s formula, its determinant is
[TABLE]
where is the adjugate matrix of .
At for , : it is the value of the characteristic polynomial at , where the last inequality follows from Theorem 3.8.
Since is an irreducible matrix, for each there is a such that the entry of is positive [21, XIII.1.2]. For , as is the sum of nonnegative integer powers of , it has positive entries. Thus has positive entries. This property persists for , where 1 is the dominant eigenvalue of (see [21, XIII.2.2]).
Finally, for and , using Jacobi’s formula, one has
[TABLE]
by positivity of the entries of and monotonicity of the entries of .
Thus both terms of Eq. 13 are nonpositive. When , the first one becomes 0. The second one remains negative: is nonzero, therefore its product with has positive entries and at least one of those of is nonzero by irreducibility, and then it is negative. ∎
Example 16. Consider again Burris’ example from Eq. 9. The characteristic system in that case is
[TABLE]
If one takes for starting point and the corresponding values given by 2, the first few iterates of for the Newton iteration solving the characteristic system are (the value of is always equal to that of )
[TABLE]
The fixed point is thus reached in 11 iterations. That this solution indeed corresponds to the radius of convergence is checked by verifying that at this point (by Theorem 3.8). Note that some of the intermediate values computed during the iteration do not correspond to values of the generating functions (for instance, for small negative , they are negative), and that the iteration uses values of outside the domain of convergence of the generating functions. With the starting point , another solution of the system is reached, namely , which has negative coordinates. Empirically, the radius of convergence is found with any initial point such that .
Theorem 4.8 gives an efficient iteration, but as usual with this type of result, the choice of initial point remains an issue. For constructive systems, experiments suggest to start with a point with and computed by Newton iteration. As is unknown, this seems like a circular argument, but thanks to Theorem 4.1 and Theorem 4.4, one can use dichotomy to find such a point, even when the domain of application of Theorem 4.8 is quite small, as in the following.
Example 17. The grammar
[TABLE]
translates into the system of generating function equations formed of and the irreducible system
[TABLE]
The difficulty here is that the singularity of the logarithm in the equation for is at , which is quite close to the radius . Our implementation first uses a dichotomy to find with 2 digits of precision, leading to the starting point , but there the iteration fails by leaving the domain of convergence of the logarithm in one iteration. Next, a new dichotomy with 4 digits of precision gives a starting point 0.3089 that fails again. Finally with 8 digits, the point is satisfactory and 20 digits of are found in 5 iterations. Correctness of the result is obtained by checking the value of at that point.
Remark 4.9*.*
In practice, the computation of the derivatives of can be performed numerically with the adjugate , even close to the point where , through a method based on the singular value decomposition of [43].
By Lemma 3.6, Theorem 4.8 can also be used for irreducible systems constructed from components of larger systems. With the notation of the lemma, the iteration is performed on the system . At each step of the iteration, the values of the variables are computed at the current value of and then used to evaluate and its derivatives. In practice however, it is easier and empirically more efficient to perform the iteration on the whole system. The convergence is still quadratic, by the following result.
Theorem 4.10**.**
Let be a well-founded system that splits into two sub-systems satisfying the hypotheses of Lemma 3.6 and assume that the radius of convergence of its generating function solution is smaller than that of and is such that is inside the domain of convergence of . Then the conclusion of Theorem 4.8 holds.
Proof.
The proof follows the steps of that of Theorem 4.8. The matrix has a block structure
[TABLE]
where is an irreducible matrix corresponding to the second sub-system and is invertible at for . The block structure of induces a block decomposition of :
[TABLE]
In the factorization , both factors are positive: the first one by invertibility for and the fact that is 1 at , since is nilpotent at 0; the second one for by irreducibility of . This shows that for . Moreover, as in the proof of Theorem 4.8, irreducibility implies that has positive entries, even for . In the case of , only nonnegativity of the entries hold in general. As a consequence, .
The derivatives also decompose by block: for and ,
[TABLE]
It follows form the previous discussion on signs that both terms are .
Thus again, both terms of Eq. 13 are nonpositive. When , the determinant of in the first one tends to 0. In the second term, all factors have constant sign or are 0. Focusing on the last block of variables, as and at least one of the derivatives for is nonzero, it shows that this factor remains negative there, proving the invertibility of . ∎
4.5 Newton Iteration for Other Singularities
Systems of generating function equations coming from constructive systems may contain equations of the form or , where is a coordinate of a vector defined by an earlier strongly connected component . If it has been detected that at its radius of convergence , then the radius of convergence of is the unique positive real solution of . Following the same reasoning as for Theorem 4.8 shows that this can be solved by Newton iteration as well.
Theorem 4.11**.**
Let be an irreducible system, let be its generating function solution and its radius of convergence. Let be a coordinate of such that . Then there exists such that and a neighborhood of such Newton’s iteration for the system converges quadratically to when started in this neighborhood.
Note again that if Newton’s iteration for the system converges to a real positive value , then computing can be used to certify that , by Theorem 3.8.
Proof.
The proof is a simpler version of the previous one. Without loss of generality, assume that is the first coordinate of . The Jacobian matrix of the system has the block structure
[TABLE]
where . By Schur’s formula, the determinant of 14 is
[TABLE]
As in the proof of Theorem 4.8, at for , has positive entries, while has nonnegative entries. This proves that Eq. 15 is positive and therefore that the Jacobian matrix is invertible. ∎
The choice of initial point for constructive systems is dealt with as above. Again, this can be used on a system that is not irreducible but is defined by the construction of Lemma 3.6.
5 Radius of Convergence of Constructible Generating Functions using Newton Iteration
The previous section gives properties of Newton’s iteration for general well-founded systems. These properties are now exploited in the case of constructive systems (Definition 2.5), which restricts the possible analytic functions that are involved.
Theorems 3.14 and 3.16 show that the radius of convergence of a constructible generating function is computable with the Constant Oracle of Definition 3.15. The algorithm underlying the proof is a dichotomy relying strongly on the oracle and therefore not very efficient. Using Newton’s iterations from the previous section leads to an algorithm that still relies on the oracle as those of Section 3.4, but it uses numerical approximations to make many of its decisions. Only equalities are delegated to the oracle. For that reason, the algorithms presented in this section are labeled as ‘Numerical Algorithms’.
The data-structure they use to represent a radius of convergence is made of a system of equations in the radius and possibly other variables, together with an interval for each of its variables, so that the system has only one solution when each variable lies in the corresponding interval. As these intervals isolate one of the possibly many solutions of the system, they are called isolating intervals. Together with the algorithms of the previous section, this data-structure allows one to compute the solutions to arbitrary precision, limiting to equality questions the tasks devoted to the oracle.
Our algorithms are designed for constructive systems in normal form (see Definition 2.6). The main algorithms use subroutines that take as input a component from the decomposition given by Eq. 6. Each of these components can be written as a system in normal form where is the (possibly empty) set of its predecessors in the decomposition. The sub-system defining is also a constructive system in normal form.
For such a system, testing that a point belongs to the domain of convergence of the generating function of is easy: since the system is constructive and in normal form, the domain of convergence is defined by the inequalities for , plus the inequalities for any such that or is the right-hand side of an equation of .
The main result of this section is the following.
Theorem 5.1**.**
Given a Constant Oracle as in Definition 3.15, Algorithm 5 decides whether the generating function solution has a finite radius of convergence or not; if it is finite, it returns a system of equations defining it and a tuple of isolating intervals (one for each variable of the system) where the solution is unique.
We first illustrate the working of the algorithm on a few examples, then deal separately with the basic cases of nonrecursive equations and of irreducible systems, before giving the algorithm that combines them, using the decomposition of systems into their strongly connected components.
5.1 Examples
Example 18. Consider the system (7) satisfied by the generating functions from Definition 2.6. The algorithm works in a bottom-up way on the strongly connected components.
First Component .
The generating function of is defined by the single nonlinear recursive equation
[TABLE]
According to Theorem 3.8, its radius of convergence is either or a solution of the characteristic system . This system has a unique solution such that and are real positive, for which . Thus, by Theorem 3.8, the radius of convergence for this component is
[TABLE]
Since one of the equations is linear, Newton’s iteration from Theorem 4.8 converges in 2 steps to this value. Our algorithm returns the system and a pair of intervals around where it has a unique solution .
Second Component .
The generating functions of and are defined by the irreducible system
[TABLE]
Again, one can use Newton’s iteration on the characteristic system. Starting with the value of at , one gets . As the largest eigenvalue of the Jacobian matrix at is smaller than , by Theorem 3.8 again, the radius of convergence satisfies
[TABLE]
Besides the characteristic system, our algorithm returns this interval and the corresponding intervals where the characteristic system has a unique solution.
Third Component .
The corresponding system of generating function equations is the union of the previous two and
[TABLE]
Again, a consequence of Theorem 3.8 is that the radius of convergence of this system is either on the boundary of the domain of convergence of the system, that is the minimum of and , or the equation is satisfied. The numerical Newton iteration for this system shows that . Thus there exists a solution with real positive values and for the corresponding characteristic system , with the determinant of the Jacobian matrix of the system coming from ,
[TABLE]
Starting with the value of at and using the iteration of Theorem 4.10 gives
[TABLE]
together with
[TABLE]
which gives a domain where the characteristic system has a unique solution. In this purely algebraic setting, one can also use elimination to produce an irreducible polynomial of degree 21 that has as its unique root in the interval above, namely
[TABLE]
Last Component .
Since the generating function of is , the radius of convergence is and it arises from the irreducible component defining . Thus it is returned with the system from that component.
Example 19. We change the system for the component in the previous example by adding cardinality constraints on the nodes of the trees, so that the system for its generating functions becomes
[TABLE]
The first two components and are identical. The computations on the third component are done in the same way as above and yield that and, as before, it is the radius of convergence of the whole system, since . But let us change again the system, to consider ordered trees, setting . We now have to look for a singularity coming from whose generating function is , since, in this case, . Using the numerical Newton iteration to compute a value of such that following Theorem 4.11 and starting from the value of at , we find that the radius of convergence is the unique solution of the system with
[TABLE]
and
.
Example 20. We change again the system for the component so that the generating function equations become
[TABLE]
The Jacobian matrix for this system is
[TABLE]
In this case, the numerical Newton iteration for converges at , where the largest eigenvalue of is , which implies that , from the component . To conclude this last example, recall that , thus the radius of convergence of the whole system is itself and it arises from the irreducible component defining .
5.2 Nonrecursive Equations
The simplest situation for the computation of the radius of convergence is the case of a nonrecursive equation. Then Algorithm 3 is a simple distinction of cases, based on properties of the functions , and . This algorithm is a refinement of ‘Algorithm Radius’ in the old work of Flajolet, Salvy, Zimmermann [18, p. 77]. The main difference is that now, the equation arises inside a larger system, thus information concerning the generating functions on which it depends (denoted in the algorithm) and that have been analyzed before, forms part of the input of the algorithm.
Deciding at line 3 is obtained by comparing 1 with the solution of a system defining . Such a system is obtained by starting from the system defining and adding the equations for the coordinates of that are not in already. The uniqueness of the solution of in the given intervals implies the uniqueness of the solution of the extended system, whose coordinates contain and the values of the , and thus in particular . Comparing the latter value to 1 is achieved by checking whether equality holds, using the Constant Oracle, and refining the interval numerically if they are distinct.
In the pseudo-code, the function NewtonIteration takes as input a combinatorial system and a value of . It implements the iteration of Theorem 4.1. At each iteration, it checks that: the coordinates are increasing; their values belong to the domain of convergence of the corresponding ; the corresponding value of is at most 1. It exits and returns fail if one of these conditions is not satisfied.
5.3 Irreducible Components
Let be an irreducible component input. Algorithm 4 puts together some of the results of the previous sections.
If the radius of convergence of the solution is smaller than the radius of the generating function , then by Corollaries 3.9 and 3.10 either is infinite at (then the system is linear, by Lemma B.2) or it is finite and is a solution of the characteristic system with .
If the system is linear, it is of the form . The generating function solution is . It is singular if and only if either is singular or at the solution of . The latter case is precisely the system computed in Algorithm 4. The hypothesis that is convergent at amounts to and being convergent there, which implies that . The extra condition ensures that by Corollary 3.10.
If the system is not linear, the system computed in Algorithms 4 and 4 is the characteristic system of Eq. 8. If , it has a solution where is convergent and . Then Corollary 3.10 shows that is the radius of convergence for this component.
5.4 General Case
Finally, the general case is described in Algorithm 5. It first decomposes the system into its strongly connected components and then deals with each of them in order, using either of the previous algorithms, depending on the nature of the component.
The isolating intervals that are returned can be refined to any precision using standard numerical solvers. Note that the cases with infinite radius of convergence are straightforward by Proposition 3.12. Ultimately, they are detected by Algorithm 3.
6 Dominant Singularities of Irreducible Systems
The algorithms of the previous section give access to the radius of convergence of the generating function solutions to a combinatorial system.
If is a generating function, the knowledge of its radius of convergence gives the exponential order of the growth of the sequence of coefficients in the sense that
[TABLE]
Since the coefficients are nonnegative, by Pringsheim’s theorem [44] (see also [45, 25]), itself is a singularity of . When it is the only singularity on the circle of convergence, then the estimate above can be greatly refined into an asymptotic expansion by a further analysis of as , as recalled in Section 9. An oscillatory asymptotic behavior occurs when several singularities lie on the circle of convergence. The singularities whose modulus is the radius of convergence are called dominant singularities. When there are finitely many of them, as will be the case in our combinatorial applications, again a local analysis can be performed to capture this oscillatory behavior. We now turn to the determination of these dominant singularities and their arguments.
Sequences of combinatorial origin often have dominant singularities with strongly restricted arguments. For instance, Berstel showed than an -rational function (i.e., one obtained from positive real numbers and a variable by additions, multiplications and applications of ) has all its dominant singularities with arguments that are rational multiples of [6, 40].
Note that positivity by itself is not sufficient for this property to hold. A typical example is the following, related to Chebyshev polynomials [6].
Example 21. For all integers and such that , the rational function
[TABLE]
has positive integer coefficients. It has three dominant singularities, at and , with that is not a rational multiple of in general.
This property — arguments of dominant singularities are rational multiples of — holds for solutions of well-founded systems comprising only one equation, either when they are not recursive and use only the basic species of Table 1 [18, Prop. 4.4], or when they are recursive under an extra assumption [4, Lemma 26]. Moreover, the generating function solution in those cases can be written in the form and the dominant singularities have complex arguments that are integer multiples of ; the integers can be computed directly from the combinatorial definition.
For irreducible systems with more than one equation, that same decomposition property actually holds in the nonlinear case as well. This has been proved by Bell, Burris and Yeats in their study of supports of power series [3, Thm. 4.6].
Definition 6.1**.**
The support of a nonzero power series is the set of indices of its nonzero coefficients. The valuation of is . Its period is . By extension the period of a power series reduced to one monomial is 0.
Regarding generating functions, the definition of the valuation corresponds to the one given for species in Section 2.3.
Theorem 6.2**.**
[3, Thm.4.6, Cor.4.7]** Let be an irreducible system and let be its generating function solution. All coordinates have infinite support and are periodic with the same period . If moreover is not linear, then they are of the form with a power series with at most finitely many coefficients equal to 0.
Proof.
Bell, Burris and Yeats prove this result when and . By the reductions of Section 3.1.3, all well-founded systems can be transformed into systems with those properties, but these reductions do not preserve valuations and periods. Still, the arguments used in their proofs apply. We sketch the proof here and refer to their article [3] for details.
For all , let be the valuation of so that with . Note that a Taylor expansion of at shows that for all . The system
[TABLE]
is well founded: that is a power series follows from the valuations above; the coefficients of are nonnegative since they are formed by positive combinations of coefficients of ; the are well defined because itself is well founded; is nilpotent.
By Lemma B.1, we can assume that for all . Then, the Taylor expansion of contains a nonzero monomial with . Since for all , implies that and , so that the period of , which is that of , divides that of , which is that of . Since this holds for all , they are equal. Applying this argument with , one can further assume that there is an index such that , since otherwise the Jacobian matrix would not be nilpotent. For such a , it follows that , which shows that this support is infinite and then so are those of the other by the inclusion above.
Let be the common period of the . When is not linear, a similar reasoning shows that for some and . The conclusion then follows from [3, Lemma 2.13]. ∎
As shown by Section 6, the theorem above is not sufficient to conclude that the dominant singularities have arguments that are rational multiples of . A result of that type is announced for nonlinear irreducible systems by Bell, Burris, Yeats [3, Prop. 7.2(i)], with a one-line proof from which we could not reconstruct their intent and that seems to be missing a hypothesis. A proof based on Wielandt’s theorem is given below.
Theorem 6.3**.**
In the conditions of Theorem 6.2, let be the period of the coordinates of the generating function solution . If the system is nonlinear, then the dominant singularities are of the form , .
Theorem 6.4**.**
Let be an irreducible linear well-founded system. Let , and be the radii of convergence of and let , and be their sets of dominant singularities.
If the largest real eigenvalue of is smaller than 1, then let and . Otherwise, there exists a unique such that the largest real eigenvalue of is 1. Let be the gcd of the periods of the entries of if and the gcd of the periods of the entries of the diagonal of otherwise and let
[TABLE]
Then, and
[TABLE]
The rest of this section is devoted to a proof of these theorems.
6.1 Periods and Dominant Singularities
Nonnegativity of the coefficients of combinatorial generating functions has a strong impact on their growth. The first one is a consequence of the triangular inequality.
Lemma 6.5**.**
[19, “Daffodil Lemma”]** Let be a power series with nonnegative coefficients , two of which at least being nonzero. Let be such that converges at . If is such that , then there exist relatively prime integers with such that and has period .
The next statement is purely at the level of formal power series. See [3] for similar results.
Lemma 6.6**.**
Let (resp. ) be nonzero power series in one (resp. ) variables with nonnegative coefficients and satisfy for . If has period , then the periods of are all multiples of .
Proof.
Let be the valuations of the power series and their periods. Three basic cases are multiplication by a power of , which follows from the definition of period, sum and product. Without loss of generality assume . Then and
[TABLE]
If has period , the gcd of the elements of this set is . In particular, , as gcd of the elements of the first set in the right-hand side, is a multiple of . Then so is since and finally so is .
Similarly, the product has valuation and
[TABLE]
Again, since both sets contain 0, if this is a subset of , then each of them must be, implying that the gcds of their elements are both multiples of .
By induction on the number of variables, the lemma then holds for each monomial in the variables and then for any sum of these monomials and thus for . The property applies to all because . ∎
6.2 Dominant Eigenvalues
The study of dominant singularities that are not real is in turn related to eigenvalues of the Jacobian matrix outside the positive real axis. The notation is now extended to denote the maximal modulus of an eigenvalue of at . A basic tool is Wielandt’s theorem [21, §XIII.2.3, Lemma 2].
Lemma 6.7** (Wielandt’s Theorem).**
If is a nonnegative, irreducible square matrix and is a matrix with complex entries of the same dimension such that for all , then any eigenvalue of has modulus bounded by that of the dominant eigenvalue of . If for some one has , then can be written with and a diagonal matrix with entries of modulus 1.
This leads to the following result for irreducible systems.
Lemma 6.8**.**
Let be an irreducible system. Let be the radius of convergence of the generating function solution and let be the Jacobian matrix of evaluated at . If is a singularity of with , then . If , each entry of has modulus equal to that of the corresponding entry in .
Proof.
First, for any nonnegative real and real , if converges at , then it also converges at by absolute convergence.
The Jacobian matrix also has a Taylor expansion at 0 with nonnegative coefficients. Thus, again by absolute convergence, for positive real and real , all the entries of have modulus bounded by the corresponding entry in . This implies that the spectral radius of is at most that of by Wielandt’s theorem. By continuity of and , this holds when and : , where the last inequality is provided by Corollary B.6.
If , then the previous inequality implies that and each entry of has modulus equal to the corresponding entry in , again by Wielandt’s theorem. ∎
6.3 Proof of Theorem 6.3
By Lemma B.2, all coordinates are finite and then so are the .
If for all these inequalities are equalities, then the “daffodil lemma” 6.5 concludes. Otherwise, one of these inequalities is strict and then by absolute convergence, lies inside the domain of convergence of . Then being a singularity implies that by the implicit function theorem. Lemma 6.8 then gives that each entry of has modulus equal to that of the corresponding entry in . Since is zero-free (by Definition 3.2) and not linear, one of the entries of depends on a coordinate that has infinite support by Theorem 6.2. Then Lemma 6.5 applies and for a period of that entry. Next, Lemma 6.6 implies that that coordinate itself has period a multiple of and then so do all the coordinates by Theorem 6.2.
6.4 Proof of Theorem 6.4
By the last point of the Definition 3.1 of well-founded systems, the matrix is invertible in a neighborhood of . It follows that the generating function solution is
[TABLE]
and it is sufficient to consider the singularities of both terms in this product and keep those that have minimum modulus. The second term is , with singularities . If the largest eigenvalue of is smaller than 1, then is invertible on and its dominant singularities lie in . Otherwise, since the system is well founded, the largest eigenvalue of is 0 at and larger or equal to 1 at , thus by monotonicity and continuity there is a unique where the largest eigenvalue is 1.
If is a dominant singularity of that is not a singularity of , then has 1 for dominant eigenvalue. Then, Lemma 6.8 applies to the system where in a situation where . Thus and . This shows that the entries of have modulus equal to the corresponding entries in . Let be the gcd of the periods of the entries of . If , then and one of the entries of has a nonzero period and Lemma 6.5 then implies that for some integer .
If then by definition of the periods, the entry of is a monomial for all , so that . Since has 1 for dominant eigenvalue, Wielandt’s theorem implies that there exists a matrix such that . If , extracting entries of this matrix identity yields
[TABLE]
A circuit in the dependency graph of corresponds to a sequence of edges such that each edge starts from the endpoint of the previous one, the last edge ends at the starting point of the first one and each of the coefficients is nonzero. The system being irreducible, it contains at least one such circuit. Summing the identity above along a circuit, the right-hand side telescopes, which shows that is an integer multiple of . Taking the gcd over all circuits concludes the proof. The sum over all circuits starting at vertex of the products of over the edges of the circuit is exactly the th entry in the diagonal of (by the transfer-matrix method [41, §4.7]).
The statement of the theorem is obtained by taking the union of all these sets of potential singularities on the circle .
7 Dominant Singularities of Constructible Generating Functions
Recall that constructive systems (Definition 2.5) are well-founded systems defined using the basic constructions . This specific structure, coupled with the results of the previous section, lets us design algorithms that determine the argument of the dominant singularities other than the radius of convergence.
In the case of nonlinear irreducible systems, Theorem 6.3 shows that the arguments of dominant singularities are of the form where is the period of the generating functions. We first show how to compute this period. The general case of systems that are not necessarily irreducible then follows from the block-decomposition of systems into irreducible components and single nonrecursive equations, with some care for linear cases. This leads to the following analogue of Theorem 5.1 for dominant singularities.
Theorem 7.1**.**
Given a Constant Oracle as in Section 3.4.3, Algorithm 9 computes a finite superset of the dominant singularities of the generating function solution of a constructive system in normal form.
The proof of this theorem splits into several cases depending on the type of the system, detailed in the following subsections.
7.1 Valuations and Periods
The computation of the periods is intertwined with the computation of valuations that was done in Algorithm 1 on page 1.
Denoting by the period of the power series , the following properties are simple consequences of the definitions [3, Prop. 2.9]:
[TABLE]
Thus, following the same template as for the valuations, Algorithm 6 computes the periods. Termination is clear: as the system is constructive, it is zero-free and thus the valuations are all nonnegative integers. Then the iteration taking gcds can only reduce coordinates a finite number of times. This algorithm is slightly unnecessarily technical by taking as input a component and a sub-system, but it is convenient to recover only the periods of the component in other algorithms below.
Note that both algorithms apply to systems defining vectors as well as matrices. This is used in Algorithm 8.
Example 22. Consider the system
[TABLE]
In order to simplify the presentation, we do not rewrite this system in normal form. The valuations computed by Algorithm 1 are . From there, the periods are given as solutions to the system
[TABLE]
Fixed point iteration gives Thus both series have period 2.
7.2 Nonrecursive Equations
In the case of a nonrecursive equation in normal form, the species is one of
[TABLE]
by Definition 2.6. In the first two cases, there are no singularities. In the next two cases (union and cartesian product), the set of singularities is a subset of the union of the sets of singularities of the .666It can be a strict subset. Consider for instance . The sets of singularities of the generating functions of and are and , but the only singularity of the generating function is . In the case of Set, the singularities are those of . Finally, for the last two cases, dominant singularities occur when the generating function reaches 1 on the circle of radius . By Lemma 6.5, the argument of the dominant singularities is of the form with a divisor of the period of , which is computed in Algorithms 7 to 7 of Algorithm 7.
7.3 Irreducible Systems
In the case of irreducible nonlinear systems the computation of periods above is sufficient to locate a finite set containing the dominant singularities of the generating functions by Theorem 6.3.
In the linear case, Theorem 6.4 leads to considering another set of periods.
Example 23. The simple variation
[TABLE]
of the previous example can be treated in the same way. The valuations are again . The system of equations for the periods is
[TABLE]
from where an iteration shows that the periods are . This is correct, but the generating function solution is
[TABLE]
with denominators indicating dominant singularities at for .
Indeed, the algorithm rewrites the system as , with
[TABLE]
and the periods of are those of the solutions of the system . Theorem 2.9 shows that the valuations of satisfy
[TABLE]
and thus by iteration
[TABLE]
Next, the matrix of periods satisfies
[TABLE]
Algorithm 6 starts from the 0 matrix and computes the fixed point
[TABLE]
in two iterations. The gcd of the periods on the diagonal is 5 and the location of possible dominant singularities follows.
Note that in this example, the periods and the valuations (and therefore the leading terms) are computed for a matrix. This is done by flattening the entries of the matrix into a vector in order to use Algorithms 6 and 1 with appropriate input.
Proposition 7.2**.**
Given a Constant Oracle as in Section 3.4.3, Algorithm 8 computes a finite superset of the dominant singularities of the generating function solution of a strongly connected component of a constructive system in normal form.
Proof.
When the system is irreducible and nonlinear, Theorems 6.3 and 6.2 show that the singularities are deduced from the common period of the coordinates of , computed in Algorithm 8.
When the system is irreducible and linear, the location of the dominant singularities is given by Theorem 6.4, implemented in Algorithms 8 to 8. The decision on Algorithm 8 is possible with the oracles at hand, since it was already done in the computation of the radius of convergence (Algorithm 4 of Algorithm 5). ∎
7.4 General Case
As in the case of Algorithm 5, the computation of dominant singularities proceeds by finding them for each of the strongly connected components of the system in order, restricting to those whose generating function have the same radius of convergence as the full system. The result of each strongly connected component is used as input in the next one. This is done by Algorithm 9, which relies on the Algorithms 7 and 8 depending on the nature of the component. The correctness of this approach follows from Proposition 7.2.
Example 24. The generating function has periodicity 1 but 4 dominant singularities. It is the generating function of the last component of the following system in normal form
[TABLE]
The algorithm starts by finding that the radius of convergence of the generating function is the positive root of and that it is the radius of convergence of and . Then NR-DS is called with each of the two equations defining and . For the first one, the valuation of is found to be 4 so that is returned. The second one is a union, resulting in the output .
8 Singular Behavior
While the radius of convergence of a power series conveys information on the exponential growth of its coefficients, singularity analysis obtains a more precise asymptotic behavior by first performing a more detailed analysis of the function in a neighborhood of its dominant singularities.
In their book, Flajolet and Sedgewick ask for an automatic way to do this [19, p. 493]:
It would at least be desirable to determine directly, from a positive (but reducible) system, the type of singular behavior of the solution, but the systematic research involved in such a programme is yet to be carried out.
An answer was provided by Banderier and Drmota [1] for the case of positive systems of polynomial or entire functions under mild assumptions that prevent essential dominant singularities [1]. The class of constructive systems is both smaller (not all positive entire functions occur) and larger (cycles and sequences are allowed). In this section, their results are extended to constructive systems. As exponentials and logarithms can now occur in the systems, there is a larger variety of possible singular behaviors that we classify, focusing on algebraic-logarithmic singularities.
The basis for the analysis of Banderier and Drmota and for ours is the classical Drmota-Lalley-Woods theorem [9];[30];[46];[19, Thm. VII.6], stating that the dominant singularities of ‘nice’ irreducible systems behave in a very predictable way: only square-root-type singularities occur. Section 8.1 gives a version of this result in the case of nonlinear irreducible components of well-founded systems, with few extra hypotheses. This is then used to classify the possible behaviors of all systems at their dominant singularity on the positive real axis. In particular, we recover a gap phenomenon that is known in the case of nonrecursive systems [18, Prop. 4.8]: the behavior is either of the algebraic-logarithmic type (see Definition 8.3) or at least as large as for some . Using this information on the dominant real positive singularity, the possibilities at other dominant singularities can then be classified. There, nonreal exponents can appear, in contrast with the case of entire functions considered by Banderier and Drmota.
8.1 Drmota-Lalley-Woods Theorem
We give a version of the Drmota-Lalley-Woods theorem suited to our context with an irreducible component inside a wider system, as an intermediate step to results for reducible systems. This is the source of the occurrence of Puiseux expansions in the system itself (coming from the bottom-up exploration of the DAG) rather than only in its solution.
Definition 8.1**.**
If is a power series
[TABLE]
then by extension of Definition 6.1, the rational number is called the valuation of and denoted . Similarly, we call valuation at and write for the valuation of .
Proposition 8.2**.**
Let be a nonlinear irreducible component of a constructive system, with depending polynomially on . Let be a dominant singularity of the generating function solution and . If , set ; otherwise, assume that all the generating functions such that have convergent Puiseux expansions in , and that the coefficient of in is nonzero with coordinates in .
If , then all coordinates of have an expansion in . Otherwise, and each coordinate of has an expansion in satisfying
[TABLE]
with all coordinates of and in . If and are respectively a left and a right eigenvector of for the eigenvalue 1 with positive coordinates then
[TABLE]
Note that the condition that depends polynomially on is automatically satisfied for a constructive system in normal form (see Definition 2.6).
The proof is given in Appendix D.
8.2 Algebraic-Logarithmic Singularities
The ‘transfer theorems’ of Flajolet and Odlyzko [17] give precise information on the asymptotic behavior of the coefficients of a generating function given its behavior in the neighborhood of its dominant singularities (see Section 9). Of particular importance is the class of functions with algebraic-logarithmic expansions in the following sense.
Definition 8.3**.**
A function analytic in a slit neighborhood of is said to have an algebraic-logarithmic singularity at if it satisfies
[TABLE]
with and in , .
The results of Flajolet and Odlyzko also hold for exponents that are not natural integers and for iterated logarithms, but these will be proved not to appear in expansions of constructible generating functions (Theorems 8.5 and 8.6).
Constructible generating functions with finite radius of convergence satisfy a gap property: either they have algebraic-logarithmic expansions at their dominant singularities or they have superpolynomial growth there, in the following sense.
Definition 8.4**.**
The function is said to have superpolynomial growth at if
[TABLE]
This gap property is known for the special case of nonrecursive equations [18, Prop. 4.8]; Theorems 8.5 and 8.6 show that it holds for systems, while giving more information on the full local expansion.
8.2.1 Real Positive Dominant Singularity
The positivity of the coefficients of generating functions at the origin induces strong constraints on their behavior at their dominant real positive singularity.
Theorem 8.5**.**
Let be a constructive system and let be the real dominant singularity of the generating function . Let and . As , each coordinate of is either of superpolynomial growth at or can be expressed as a finite sum of expansions of one of the types
[TABLE]
In both cases, is of the form with .
Moreover, when the system does not use the Cyc construction, then does not appear in the expansions ( and ), the exponents that occur in the second type of expansion are rational with denominator and the corresponding sum belongs to .
The second part of the theorem recovers the result of Banderier and Drmota [1, Thm. 5.3] for constructible generating functions. Examples where is 4 or 8 are given by supertrees and supersupertrees [19, p. 412–414].
Example 25. This is a simplified variant of [18, Ex. 4.9]. It gives an example of an expansion of the second type, with an exponent that is not algebraic. The species
[TABLE]
has for generating function
[TABLE]
At the dominant singularity , its behavior is
[TABLE]
The exponent is transcendental by Lindemann’s theorem.
Theorem 8.5 is a classification theorem. Its main point is that no other behavior is possible.
Notation**.**
The set of functions with superpolynomial growth at is denoted . Similarly, denote the sets of functions whose behavior at is a finite sum of expansions of the first type in the theorem, with the maximal value in these expansions. (Functions that are analytic at belong to with .) Finally, denotes the set of functions with an algebraic-logarithmic singularity at with whose behavior at is a finite sum of expansions of the first and second type, with the minimal value and the maximal one for this value of .
Also, we write and .
Proof of Theorem 8.5.
Without loss of generality, we consider the system to be in normal form: the coordinates of the original system form a subset of those of the normalized system.
The proof is by induction on the strongly connected components of the decomposition of the system. The induction hypothesis is that and moreover, if a nonconstant component has an expansion in , then the coefficient of in that expansion is negative (and in particular nonzero).
Irreducible component. In the nonlinear case, the result is a consequence of Proposition 8.2.
Otherwise, the component is of the form , with having entries with expansions in , which is a field. Therefore the coordinates of the solution also belong to . By Lemma B.2, either they all tend to infinity or they all have a finite limit. In the first case, they belong to with rational and with denominator . The second case is when . The expansion is given by
[TABLE]
with , . By irreducibility, the entries of are all positive; by the induction hypothesis at least one of the entries of and has a nonzero coefficient of , which is moreover negative, and thus so do all the entries in the product.
Nonrecursive component. The system is .
If , then the expansions are and in , that satisfy the induction hypothesis.
Sums and products. Let be analytic functions in a slit neighborhood of . Then the cases when is a sum or a product follow from
[TABLE]
In the last case, if both belong to with a negative coefficient of , then so does their product and sum.
Sets. If with or , then . If , then by composition with a convergent Taylor expansion. The coefficient of being negative is preserved by the exponential.
The remaining case for is when (induced by a Cyc in a previous component). Then we have
[TABLE]
with and , and a possibly nonzero with . The exponential of is and thus belongs to .
Sequences and cycles. When with , the generating function belongs to and its limit at is at most 1. If , then by composition with a Taylor series, and the coefficient of is negative. Otherwise, has an expansion
[TABLE]
and the leading coefficient of is 1. Then, for ,
[TABLE]
and has so that this is an expansion in . In the case of one obtains
[TABLE]
which is the sum of the expansion and the expansion , concluding the proof. ∎
8.2.2 Other Dominant Singularities
In view of the application of transfer theorems to the local expansions in the next section, we focus here on generating functions whose behavior along the positive real axis is algebraic-logarithmic. This behavior persists at the other dominant singularities, i.e., they do not have superpolynomial growth there. Still cancellations and complex coefficients can occur in the expansions and in the exponents at nonreal singularities, leading to more diversity in the possible local behaviors.
Theorem 8.6**.**
Let be a constructive system and let be a dominant singularity of the generating function with . Let and . Let be a coordinate of that has an algebraic-logarithmic singularity at . Then as , its behavior can be expressed as a finite sum of expansions of one of the types
[TABLE]
In both cases, is of the form with .
When the system does not use the Cyc construction, then does not appear in the expansions, the exponents that occur in the second type of expansion are rational with denominator and the corresponding sum belongs to .
The main differences with the real positive case of Theorem 8.5 is that in the first type of expansions, the leading term may have to a positive exponent and in the second type, can be a nonreal complex number and the coefficient may be 0.
Example 26. Consider the species . Its generating function has a complex exponent at its singularity :
[TABLE]
This imaginary exponent is the basis for a term in in the asymptotic behavior given in Proposition 9.1 below.
Notation**.**
The set of functions whose behavior at is a finite sum of expansions of the first type is denoted , with the maximal value reached in these expansions. Those whose behavior is a finite sum of expansions of the first and second type is denoted with and as above. Again, and .
Proof of Theorem 8.6.
The proof is similar to the previous one, by an induction on the strongly connected components of the decomposition of the system. The new ingredient is that the behavior of at is dominated by its behavior at , given by Theorem 8.5. In particular if is analytic at , it is also analytic at . The induction hypothesis is that .
Irreducible component. The reasoning is the same as in the real case. If the component is linear, the results follows from the fact that is a field. If the component is not linear, then if , the result is a consequence of Proposition 8.2. Otherwise, and is a singularity of . By Theorem 6.2, all coordinates are of the form with such that for some . Thus, their behavior as is given by their behavior as ; Theorem 8.5 then implies that .
Nonrecursive component. The system is .
If , then the expansions are and in .
Sums and products. The cases when is a sum or product follow from
[TABLE]
Sets. If and , then . Following the different cases in the proof of Theorem 8.6, this implies that belongs to and therefore so does by composition with an analytic function.
The remaining case for is when . Then . By the induction hypothesis, this implies . If then by composition with also . Otherwise,
[TABLE]
with and , and a possibly nonzero with . The exponential of is and thus belongs to . This is where exponents that are not rational with denominator a power of 2 can occur when Cyc is used.
Sequences and Cycles. If with and is not analytic at while it belongs to , then it has limit 1 at and at . By Theorems 6.3 and 6.2, with and with a power series at 0, which is then analytic for . Since at , it follows that , which implies that . Thus the behavior of at is given by its behavior at and by Theorem 8.5,
[TABLE]
with , so that and . ∎
8.3 Algorithm
Theorems 8.5 and 8.6 become effective in Algorithms 11, 12 and 13, whose properties are summarized in the following.
Theorem 8.7**.**
Let be a constructive system in normal form (Definition 2.20) and assume that lies inside or on the boundary of the domain of convergence of the generating function . Algorithm 13 computes truncations of the expansions of all coordinates that are analytic or have an algorithmic-logarithmic singularity at (Definition 8.3). It returns Exp for the other coordinates.
Proof.
The algorithm first decomposes the system into its irreducible components. It then takes each component in turn and deals with systems where may depend on auxiliary variables coming from earlier components. The proof of correctness is by induction and the induction hypothesis is that when is not empty, all its coordinates have expansions at given by Theorems 8.5 and 8.6, for which truncations have been computed recursively.
The nonrecursive cases are dealt with by Algorithm 12, the recursive ones by Algorithm 11. In both cases, the algorithm follows the steps of the proofs of Propositions 8.2, 8.5 and 8.6. ∎
Note that this result does not require any help from an oracle: all difficult decisions have been made during the computation of the dominant singularities.
8.4 Example
For the last time, we consider the system of Eq. 4 given in Definition 2.12. Again, we proceed along the strongly connected components. Recall from Section 5.1 that the radius of convergence is and it arises from the irreducible component defining .
First and Second Components .
Both components match the condition at ligne 11 in Algorithm 11. As is smaller than the radius of convergence of , they have a Taylor expansion there.
In this example, the equations defining the generating functions are polynomials of small degrees, so that these expansions can also be obtained from closed-form expressions:
[TABLE]
In general however, one would only have access to numerical approximations of the values of with arbitrary precision obtained by Newton’s iteration from their defining systems in Section 5.1. From there, the expansions are obtained using either undetermined coefficients or Newton iteration over power series, giving
[TABLE]
where . We display only the first 6 digits of each coefficients and expansions to arbitrary order are easily computed.
Third Component .
This is the singular nonlinear case at line 11 in Algorithm 11, from where the dominant singularity occurs. By Proposition 8.2, the expansion is of the form
[TABLE]
and the values of have been computed by Newton iteration in Section 5.1.
Again, since the equations are polynomial and of low degrees, the system can be solved symbolically in that case, leading to
[TABLE]
The computation performed by the algorithm does not exploit the possible existence of such closed forms. Instead, the expansions are computed in two stages. First, is replaced by where will be used to perform expansions in . Replacing by their expansions (now in ) and by in the last two equations of the system Eq. 16 produces a system with unknowns and variables whose Jacobian matrix is invertible at (see the proof of Proposition 8.2). It can thus be solved in power series in with coefficients that are polynomials in , e.g., by Newton iteration over power series, giving
[TABLE]
where again we display only the first 6 digits of the coefficients and expansions to arbitrary order are easily computed.
Next, these expansions are substituted into the first equation of Eq. 16. Expanding in gives
[TABLE]
The coefficients and are numerically small and are known to be exactly 0. They can thus be discarded, the resulting expansion can be divided by and can then be solved for , e.g., by Newton iteration for power series, once the correct root of the leading term (the negative one) is selected. This gives
[TABLE]
whence
[TABLE]
which can be injected into Eq. 23, giving
[TABLE]
Last Component .
The generating function of is . The expansion is obtain by composition of that of with that of , giving
[TABLE]
9 Asymptotic Behavior
Once the dominant singularities of the generating functions have been computed, together with the local behavior of the generating functions at these dominant singularities, the asymptotic behavior of their coefficients is determined by the method of singularity analysis, described for instance in Chapter 6 of the book by Flajolet and Sedgewick [19]. A version of these results that is sufficient for this article is the following.
Proposition 9.1**.**
(Direct corollary of [19, Thm.VI.5]) If is analytic in , has finitely many singularities on , admits an analytic continuation in a disc with slit along the segments and behaves locally as
[TABLE]
then the asymptotic behavior of is given by
[TABLE]
where
[TABLE]
Example 27. Let be the number of forests with nodes from our running Definition 2.6. Its asymptotic behavior can now be deduced directly from the expansion of Section 8. Indeed, from the expansion of a generating function at its dominant singularity of the form
[TABLE]
with , Proposition 9.1 shows that the th coefficient behaves asymptotically like
[TABLE]
Expanding by Stirling’s formula shows that this is asymptotically
[TABLE]
With the values of Eq. 27, this leads to
[TABLE]
with root of the polynomial from Eq. 18. The coefficients, and , can be computed to arbitrary precision efficiently.
Already with as small as 10, one has the exact value , while the first term of the asymptotic expansion above at is approximately , giving a relative error of . Using one more term gives and a relative error of . Using all three gives and a relative error of . These errors decrease to for .
Example 28. Theorem 8.6 leads to the generating function
[TABLE]
that has 4 dominant singularities at . Its local behavior at these singularities is
[TABLE]
and the conjugate of the last one as . These expansions contribute respectively to the asymptotic behavior of the coefficient of in for
[TABLE]
and the conjugate of the last one. Adding these contributions, expanding them and beautifying the result with the value of [36, Eq. 5.4.3] gives the asymptotic behavior
[TABLE]
with . Full asymptotic expansions can be derived the same way, preferably with the help of a computer algebra system.
That Proposition 9.1 applies in our context follows from the fact that if is a zero-free well-founded system, then the generating function solution is analytic in a disk of positive radius (by Proposition 3.3) and moreover, each of its irreducible components has solutions with finitely many dominant singularities by Theorems 6.3 and 6.4. In summary, the following has been proved.
Theorem 9.2**.**
Let be a constructible generating function. For all coordinates of that have an algebraic-logarithmic dominant singularity, an asymptotic expansion in the form of Eq. 28 can be computed given a Constant Oracle as in Section 3.4.3.
Again, note that when the constructive system contains neither Cyc nor Set, then the oracle is computable, while this computability is dependent on the truth of Schanuel’s conjecture otherwise.
10 Conclusion
The automation of analytic combinatorics was started many years ago [18]. An implementation was developed in Maple and many examples were treated automatically, including applications to the analysis of algorithms [16], which was an important motivation for the development of analytic combinatorics. The main limitation of that implementation was that it could only deal with a restricted set of systems of equations. Either the systems were completely nonrecursive, or the recursive equations were sufficiently simple to be solvable in closed form. The main contribution of the present article is to show how these limitations can be removed. We conclude by reviewing other parts of analytic combinatorics where more automation work is still needed.
Cardinality constraints
For simplicity, we have focused on the species in their versions without cardinality constraints. A large part of the theory naturally extends to indexed versions with a restriction on the cardinality of the sequence, set or cycle. The general case is recovered with for and and with for Cyc. Of course, for effectivity reasons, the set cannot be arbitrary. For the uses in [19], is either a finite set, or the union of a finite set and a finite number of arithmetic progressions. This is easily accommodated by the algorithms presented here. For instance, when is one of , Algorithm 1 would set to with and one of depending on the case and provided ; it would use ( if and is finite; Fail otherwise. Similarly, Algorithm 3 returns for and arbitrary and for and with finite . It uses its other branch in the other cases. Similar modifications can be made to Algorithms 6 and 7. The adjustments to 12 are equally straightforward, but slightly more technical.
Superpolynomial growth
We have focused on algebraic-logarithmic singularities and have not discussed what should be done with generating functions that have a superpolynomial behavior. For generating functions whose growth is at least exponential at their dominant singularity or that are entire and non-polynomial, one can use the saddle-point method [19, Ch. VIII]. This can be automated to some extent, thanks to Hayman’s class of functions, to which membership is easy to check and for which a first order asymptotic expansion is easily computed [27]. A full asymptotic expansion can be obtained for a more restricted class studied by Harris and Schoenfeld [26, 35]. Wyman [47] gives a very nice exposition that is more general, but more difficult to apply in practice. Already in the simple setting of Hayman’s class, from the computer algebra point of view, a difficulty lies in the handling of suitable asymptotic scales, see for instance [39, Sec. 3].
There are also generating functions whose growth lies between the algebraic-logarithmic ones and the exponential ones. In view of Theorem 8.5 and its proof, these are functions that behave at their dominant singularity (normalized here to be 1) like
[TABLE]
or more complicated variants with a polynomial in in the exponential and extra logarithmic factors at the end. We are not aware of results on the asymptotic behavior of their coefficients, which would complete a classification of the possible asymptotic behaviors of the coefficients of the generating functions of constructible species.
Ordinary generating functions
Ordinary generating functions are used for unlabeled enumeration. The dictionary for the translation of a combinatorial specification into ordinary generating functions is very similar to that of Table 1, except for and that translate into
[TABLE]
where is Euler’s totient function [19, Thm. I.1]. Flajolet and Sedgewick also have another construction PowerSet that is not a species, see the discussion in [38, §10.2].
This translation shows that differences between ordinary and exponential generating functions only matter if Set or Cyc are used. Newton’s iteration can be used to obtain numerical values of the generating functions [38], but the general analysis of the location of singularities is more involved. In many cases, there are infinitely many singularities on the circle of convergence, preventing singularity analysis to apply. However, this is not the case when the radius of convergence is less than 1. This can sometimes be detected by comparison with the radius of convergence of a related species. A typical example is the class of Pólya trees whose generating function satisfies
[TABLE]
Its radius of convergence is clearly smaller than , the radius of convergence of the generating function of the Cayley trees, solution to . Then one can compute the radius efficiently (see [34, §15.1],[19, pp. 475–477],[38, §9.4.2]) and ultimately, a full asymptotic expansion of the coefficients. In this setting, we do not have computability results like Theorem 3.14 but it would be interesting to investigate the automation of some of the techniques, along the lines of [22].
Multivariate generating functions
The combinatorial part of analytic combinatorics (the symbolic method) can easily deal with several sorts of atoms, which then translate into equations over multivariate generating functions.
As a typical example, in order to study the expected number of trees in the forests from Definition 2.6, one would replace the first equation by for a new sort of atom . This equation translates into the generating function equation whose coefficient of is the number of forests with nodes using trees , divided by . Differentiating with respect to and setting shows that the expected number of trees is given by . With the expansions computed in Section 8.4, it is easy to compute the behavior of as and from there obtain the asymptotic behavior of the numerator, whence, after division by Eq. 29, the asymptotic expansion of the expected number of trees in a forest drawn uniformly at random among all forests with nodes:
[TABLE]
Again, these coefficients are easily computed to arbitrary precision and more terms of the expansion are easily found.
Many other parameters can be studied in a similar way, see [19, Ch. III,IX]. Analytic combinatorics also allows to predict limiting distributions in many cases. Not much has been done in the automation of these aspects of the theory. A notable exception is the extension of Hayman’s class to the bivariate setting [11].
A direction that has received a lot of attention in recent years is the field of ‘Analytic Combinatorics in Several Variables’. There, the aim is to extract asymptotic information on the coefficient sequence of multivariate generating functions along a diagonal. Already the case of rational generating functions can become involved. In some cases, the approach can be automated [33, 24, 23].777https://github.com/ACSVMath/sage_acsv The recent new edition of the reference book on that topic [37] will be a source of ideas and questions for computer algebraists with an interest in this area.
Implementation
Maple worksheets illustrating some of the examples are available with this article. Still, a complete implementation of the algorithms presented here should be attempted. For a long time, many of the necessary routines were available in the Algolib888https://algo.inria.fr/libraries/ Maple library of the Algorithms Project at Inria, but these have not been maintained and do not run well on modern versions of Maple. The Maple package NewtonGF999https://perso.ens-lyon.fr/bruno.salvy/software/the-newtongf-package/ provides the tools required for the combinatorial part and for generating function expansions [38]. It could serve as a basis to develop the analytic part. We are not aware of implementations of the asymptotic part of analytic combinatorics in other computer algebra systems.
Appendix A Proof of Theorem 2.8 : Characterization by Leading Terms
See 2.8 Leading terms can be computed inductively, thanks to the following.
Lemma A.1**.**
For any species such that the composition is well defined, the leading terms satisfy the equality .
Proof.
This is a direct consequence of the definition of composition of species
[TABLE]
The smallest -structures are therefore obtained from the smallest structures. This applies to each of the coordinates of . ∎
Lemma A.2**.**
Let be a vector of species such that the iteration (3) is well defined for and let be defined as in Theorem 2.15, then for .
Proof.
The proof is by induction. First, . Next, if the property holds for , then
[TABLE]
where the second equality uses Lemma A.1 and the third one uses the fact that for any species , . ∎
Lemma A.3**.**
The composition is well defined if and only if the composition with the leading term, , is well defined.
Proof.
As , the composition is well defined when is.
Conversely, the composition with is well defined if and only if the composition with is. As , if is well defined, then so is . ∎
Lemma A.4**.**
Let and be such that the iteration starting with is well defined for and . Then the iteration is well defined for all and the value of is constant for . If is this value, then the minimal such that is at most the number of coordinates of different from those of .
Proof.
As the iteration is well defined up to , by Lemma A.3, the composition is well defined for all if the leading terms of is . By induction using Lemma A.1, leading terms of can only be built from leading therms of , thus for any . So the iteration is well defined for all .
Without loss of generality, let be minimal such that and let be the set of indices where . By definition of , for , the coordinates for do not depend on . This implies that the vector of coefficients is the same as that of the system in variables less. This new system has , as coefficients change at each iteration, and fewer variables. This construction can be repeated less times than the number of coordinates of different from those of , hence the result. ∎
The next three properties follow directly from the definitions.
Lemma A.5**.**
For any species and ,
; 2. 2.
if then ; 3. 3.
if and then .
Lemma A.6**.**
If the iteration (3) is well defined for and , then let and . Up to renumbering the coordinates, let be such that and the coordinates for . Write and consider the system . Then,
the iteration starting with is well defined for ; 2. 2.
it satisfies ; 3. 3.
the Jacobian matrix is nilpotent if and only if is nilpotent.
Proof.
Note that when , the first two points are irrelevant and the third one just states that the matrix is unconditionally nilpotent.
-
From Lemma A.4, the composition is well defined for all and, from , it follows by induction that , thus is well defined for all .
-
Let be the sequence defined by with initial value . By induction, one has and from the convergence of to , it follows that also converges to .
Next, we prove that . Inclusion is preserved by composition with , thus and all the coordinates of are non empty. And by definition of and the third point of Lemma A.5, ). On the other hand,
[TABLE]
Both equalities come from Lemma A.5 and the inclusion follows from .
Then, by induction, . Indeed, the right-hand side is the previous point. As for the left-hand side, it is true for by definition of and if the property holds for , then it also holds for , using Lemma A.1.
It follows that converges to . Lemma A.4 then proves that the convergence requires at most iterations.
- By Lemma A.2, , with defined as in Theorem 2.15 for . Let be the maximal index where and let . From , it follows that for all such that . Without loss of generality, let these coordinates be the last ones ( by definition of and the hypothesis on ) and consider the new system formed with . Since does not depend on the last coordinates, the Jacobian matrix has a top-left block that is the Jacobian matrix of at , while its next columns are filled with 0. Thus these two matrices are either both nilpotent or both nonnilpotent. Repeating this argument with the new system until is reached or no matrix is left (when ) proves the result. ∎
The proof of Theorem 2.8 can now be detailed.
Proof.
Assume that the system is well founded. Then the iteration (3) is well defined for all and converges to a species . The sequence converges in a finite number of steps: the sequence of valuations is a decreasing sequence of -tuples of nonnegative integers or infinity; once the valuations have converged, the sequence of leading terms is increasing by inclusion and contained in the polynomial species . Thus there exists such that . By Lemma A.4 one can take .
Conversely, assume that the leading terms of and are equal. Then Lemma A.4 shows that the iteration is well defined for all with . This shows that if the system is well founded, then its solution satisfies . All that remains is to prove it is well founded. From now on, without loss of generality, we assume that none of the coordinates of is 0. (Otherwise, the corresponding coordinates of are 0 for all and can be removed from the system.)
If , then by Taylor’s formula and induction, for any ,
[TABLE]
If the matrix is not nilpotent, there exists at least one nonzero structure in , and its size is [math]. Thus, there are indices and structures , , …, , such that all the are nonzero structures. There can be at most such indices that are distinct, thus for some . As none of the coordinates of is 0, one can choose a structure in , and find infinitely many structures of size of the form with belonging to the sequence . This is in contradiction with the convergence of the sequence , showing that the matrix is nilpotent and therefore that the system is well founded when .
If then the hypotheses of Lemma A.6 hold with . When , this lemma defines a new system with unknowns, . By Lemma A.6, the iteration defined by the new system is well defined, it satisfies and this is at . Thus the previous part of the proof applies, showing that this system is well founded. Extending to , Lemma A.6 then shows that , from the original system has a nilpotent Jacobian matrix at . By Lemma A.2, this is equal to , with the sequence defined in Theorem 2.15. This theorem then shows that the system is well founded, which concludes the proof. ∎
Appendix B Proof of Theorem 3.8: Dominant Eigenvalue
See 3.8
The positivity of follows from Proposition 3.3. The rest of the proof relies on the ideas described by Bell, Burris and Yeats in the analysis of their well-conditioned systems [2] (see Table 2).
B.1 Behavior at the Singularity
A basic tool in several proofs for irreducible systems is the following.
Lemma B.1**.**
Let be an irreducible system. Then its solution satisfies a well-founded system where all entries of the Jacobian matrix are nonzero and positive.
Proof.
Let be the dimension of the system. Consider the sequence of systems with defined by and . All these systems have the same domain of convergence and the same solution and so does the system , with
[TABLE]
All the entries of the Jacobian matrix of are nonzero and positive and depend on , since the effect of this sum is to add up all paths of length at most in the dependency graph of the original system, which is irreducible. ∎
The following lemma implies that all coordinates of the solution have a similar behavior at the dominant singularity . It is similar to [2, Prop. 3 (iv)] and parts of [4].
Lemma B.2**.**
Let be an irreducible system and be the radius of convergence of its generating function solution . Either all coordinates of tend to infinity as , or all of them have a finite limit. The former case occurs only for that is linear in .
Proof.
Let be the dimension of the system. Observe first that an irreducible system is well-founded and zero-free, thus since is not nilpotent, while its evaluation at with is.
Lemma B.1 implies that each of the coordinates of depends positively on the other ones and thus that if one of them tends to infinity, all the other ones do.
In the nonlinear case, we follow the argument of [4, Prop. 66]. Since is not linear in , there is an index such that the th equation in Eq. 30 contains a monomial of the form with at least one . Positivity of all the coefficients implies that is larger than that monomial. Dividing by on both sides and letting tend to thus implies that has a finite limit and therefore all of them do by the previous part of the proof. ∎
Note that the last statement of the lemma is not an equivalence. It is possible for a linear system to have generating function solutions with a finite limit at their singularity.
Example 29. Consider the equation with the generating function of the Catalan numbers. Thus, is the generating function of sequences of binary trees. This equation satisfies the hypothesis of the lemma: it is well-founded since ; the system containing only one equation, it is irreducible. This equation is linear in , but at the singularity , the generating function has finite value .
B.2 Ordering Solutions
The following is a variant of [2, Lemma 13].
Lemma B.3**.**
Let be an irreducible system and let be the radius of convergence of its generating function solution . Let . Let belong to the domain of convergence of or its boundary and satisfy .
If for some index , , then ; 2. 2.
if , then ; 3. 3.
if then ; 4. 4.
if then .
where inequalities between vectors mean that they hold for all coordinates.
Proof.
The idea of the proof is taken from the proof of Lemma 13 by Bell, Burris and Yeats [2].
We first prove that . Let be the th iterate of as in the proof of Lemma B.1. By Proposition 3.3, . Up to relabeling the coordinates, assume that for and otherwise. Since is nilpotent at (by definition of a zero-free well-founded system) and for all , it follows by induction that . Then by positivity of for , this implies that is [math] for those values of , i.e., at , is constant wrt . Thus,
[TABLE]
a contradiction, which shows that and thus .
Next, up to relabeling the coordinates, assume for and for (). If , there is nothing to prove in Part 1 of the lemma. Otherwise, by continuity, there exist , all bounded by , such that for . Without loss of generality, they can be ordered as . Since the ’s are increasing, for and the same inequality holds for the other ones by definition of . Let be as in Eq. 30, which satisfies for all . If , then
[TABLE]
a contradiction, which implies that and proves Item 1 of the Lemma.
Next, when , for all and even if ,
[TABLE]
proving Item 2.
The third part of the lemma is the case when . Then
[TABLE]
forces the inequality to be an equality and thus, by irreducibility, monotonicity and dependency in , all the arguments must be equal: for all and . The fourth part of the lemma comes from the dependency in in this inequality. ∎
Part 1 of Theorem 3.8 requires slightly more work.
Corollary B.4**.**
Assume that is an irreducible linear system. Let belong to the domain of convergence of or its boundary and satisfy . Then, and .
Proof.
If tends to infinity as , then all its coordinates do by Lemma B.2 and the conclusion follows from Item 1 and Item 3 of the previous lemma.
If is finite, the point must belong to the boundary of the domain of convergence of , which implies that . Assume by contradiction that tends to 0 as , then from it follows that the second determinant tends to infinity and therefore at least one entry of the inverse matrix tends to infinity, which shows that has an infinite limit as , a contradiction. Thus in the case when has a finite limit as , the matrix is invertible for . Subtracting from gives from which the invertibility of implies . ∎
B.3 Dominant Eigenvalue
We now turn to the rest of the theorem, starting with part 2a () for .
Lemma B.5**.**
Let be an irreducible system and let be the radius of convergence of its generating function solution . If , then .
This is [2, Lemma 10 (b)] for their well-conditioned systems. The proof adapts to our setting.
Proof.
Since lies inside the domain of convergence of , and its partial derivatives are well defined and have nonnegative values at . Differentiation of gives
[TABLE]
By the Perron-Frobenius theorem, there exists a left eigenvector of with positive (and not merely nonnegative) entries for the eigenvalue . Multiplying on the left by this vector leads to an equation over real values
[TABLE]
Since , its value at a positive is nonzero, implying that the right-hand side is positive. This in turn implies that and finally . ∎
Corollary B.6**.**
Let be the radius of convergence of the generating function , solution of the irreducible system . Let . Then .
Proof.
This follows from the previous lemma and the continuity of both and . ∎
Part 2a () of the theorem for then follows from the implicit function theorem (Section 3.2). The other direction of Part 2a is an adaptation of a similar result of Bell, Burris, Yeats [2, Prop. 19].
Lemma B.7**.**
Let be an irreducible system and let be the radius of convergence of its generating function solution . Let belong to the domain of convergence of or its boundary and satisfy . If then and .
Proof.
First, if then by Lemma B.3(1), . Since belongs to the domain of convergence of , this implies that so does . Then by monotonicity of this would imply , a contradiction.
Thus and then, by Lemma B.3(2), . If , Bell, Burris and Yeats cleverly construct a nonnegative matrix with spectral radius such that , while , which implies , a contradiction. The construction is as follows: for each , let , so that and . By the mean value theorem, for some . This rewrites as
[TABLE]
where is the row vector . The matrix is chosen as the matrix with rows . By design, it satisfies and , whence the result. ∎
Finally, if , then this lemma implies that and , while Lemma B.5 shows that , so that , proving Part 2b of Theorem 3.8 and concluding its proof.
Appendix C Proof of Proposition 3.11: Nonnegative Solutions
See 3.11
Proof.
The proof is by induction on the decomposition of the system into irreducible components. If the system consists of a nonrecursive equation , then by hypothesis the point is in the domain of convergence of , which is that of . The only possible value for is then .
If the system is irreducible and linear, then the conclusion comes from part 1 of Theorem 3.8.
If the system is irreducible and nonlinear, let denote , which has finite coordinates, by Lemma B.2. If for some , , then by Lemma B.3(1). Otherwise, . If is not on the boundary of the domain of convergence of , then the implicit function theorem implies that , while part 2(a) of Theorem 3.8 implies that . Thus either is on the boundary of the domain of convergence of or . The first case cannot occur with : is beyond the boundary of the domain of convergence. Finally, it is not possible that and . This is obtained by contradiction, following the proof of [2, Prop. 11]: by Taylor expansion and positivity
[TABLE]
and the inequality is strict for at least one coordinate, by nonlinearity. By the Perron-Frobenius theory, there exists a left eigenvector with positive coordinates of for the dominant eigenvalue . Multiplying by on the left gives
[TABLE]
Since the previous inequality is strict for at least one coordinate, it follows that , a contradiction. Thus for irreducible systems and Lemma B.3(2) shows that .
Finally, if the system is not irreducible, it splits into two sub-systems
[TABLE]
where the second is an irreducible component. Since belongs to the domain of convergence of or its boundary, belongs to the domain of convergence of and thus by induction, is at most the radius of convergence of and moreover .
By Lemma 3.6, is an irreducible system. Since belongs to the domain of convergence of or its boundary, so does and therefore belongs to the domain of convergence of or its boundary. The conditions of the first part of the proof are met and therefore and , which concludes the proof. ∎
Appendix D Proof of Proposition 8.2: Drmota-Lalley-Woods Theorem
See 8.2
Proof.
The proof is an adaptation of Drmota’s [9].
Singularity different from
When , by Lemmas 3.6, 6.3 and 6.2, the singularity is of the form , with and a common period of the coordinates of the solution. Moreover, each coordinate of the generating function is of the form with a function that is analytic in . As , so that the expansion of as is obtained by multiply its expansion as found in the next part of the proof by the expansion of .
Change of variables
By Lemma B.2, all have a finite limit at , which is a real positive number since the system is irreducible and therefore zero-free. All the such that also have a finite limit, as they appear polynomially in . The value of is defined as the solution of . Let and and define similarly, so that has coordinates that are convergent power series in . The functions and being increasing for , the leading coefficients in the expansions of and are positive. The function is solution to the equation
[TABLE]
with
[TABLE]
By design, both sides of Eq. 31 are 0 at and
[TABLE]
Regular case
If , the matrix is invertible at . Then, so is at and therefore Eq. 31 satisfies the hypotheses of the implicit function theorem over convergent power series. It follows that has convergent coordinates in , so that those of lie in , which proves the theorem in that case.
Halving of the exponent in the singular case
If , the matrix is not invertible at . Let be the dimension of the system. The simplest case is when . Then we have
[TABLE]
the last one by nonlinearity and positivity of the power series. Then the Weierstrass preparation theorem [28, Cor. 6.1.2] implies that the solution of the equation has a convergent Puiseux expansion in , giving the theorem in that case.
The case when reduces to the case when by an argument due to Drmota [9, §5],[10, p. 64-65], see also Mac Millan [31] for an ancestor of this idea. Write for and similarly for . A consequence of the Perron-Frobenius theory is that the first minors of are nonzero [21, XIII §2.5]; in particular, is invertible at . Then the implicit function theorem can be applied to solve the system for unique that are convergent power series in and . Moreover, differentiation gives
[TABLE]
Injecting these solutions into the first equation gives a single equation in the unknown ,
[TABLE]
with , a convergent power series in with nonnegative coefficients and . The situation is reduced to the case when , and it is now sufficient to show that and . Differentiation gives
[TABLE]
Since the matrix is singular, it follows that : the Schur complement of the invertible block in this singular matrix is
[TABLE]
Positivity of all the power series involved implies that . The conclusion for is that the solutions of this equation have convergent Puiseux expansions in . Therefore, so do the corresponding by composition. Thus all coordinates belong to .
Appearance of the square-root.
The Taylor expansion of with respect to , at , evaluated at , gives
[TABLE]
By Lemma B.1, we can assume that all entries of are positive. Thus the leading coefficients of each of the coordinates of the first two terms in the right-hand side are nonnegative and the valuation of each coordinate of the right-hand side comes from these two terms. Therefore the valuation of each coordinate is at most that of . Using Lemma B.1 if necessary, we can assume that all the entries of are nonzero, which implies that, on the right-hand side, the minimal valuation among the appears for each coordinate, and therefore also on the left, showing that all these valuations are identical. Let be their common value.
The next step shows that . Let be a left eigenvector of for the eigenvalue 1. It can be taken with positive coordinates by the Perron-Frobenius theorem, so that no cancellation occurs when multiplying Eq. 33 on the left by , giving
[TABLE]
By positivity, , so that the right-hand side has valuation . The linear term of the left-hand side with respect to comes from
[TABLE]
There, (with from the hypothesis) has at least one nonzero coordinate, that is positive, and the corresponding is nonzero and positive too, so that is nonzero. Thus the valuation of all coordinates of is .
Coefficients .
Let finally be such that . Extracting the terms with exponent 1/2 in Eq. 33 gives
[TABLE]
i.e., is a right eigenvector of for the eigenvalue 1. By the Perron-Frobenius theorem, there is a one-dimensional choice of such vectors. Let be a right eigenvector with positive coordinates, and let be such that . (The sign of is dictated by the fact that is increasing for .) The constant factor is obtained by considering the terms with exponent 1 in Eq. 34, giving the quadratic equation
[TABLE]
which concludes the proof of the proposition, in view of Eq. 35. ∎
Acknowledgment
This work has been supported in part by the French project NuSCAP (ANR-20-CE48-0014).
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] Cyril Banderier and Michael Drmota. Formulae and asymptotics for coefficients of algebraic functions. Combinatorics, Probability and Computing , 24:1–53, 1 2015.
- 2[2] Jason Bell, Stanley Burris, and Karen Yeats. Characteristic points of recursive systems. The Electronic Journal of Combinatorics , 17(R 121):1–35, 2010.
- 3[3] Jason P. Bell, Stanley N. Burris, and Karen Yeats. Spectra and systems of equations. In Model theoretic methods in finite combinatorics , volume 558 of Contemp. Math. , pages 43–96. Amer. Math. Soc., Providence, RI, 2011.
- 4[4] Jason P. Bell, Stanley N. Burris, and Karen A. Yeats. Counting rooted trees: the universal law t ( n ) ∼ C ρ − n n − 3 / 2 t(n)\sim C\rho^{-n}n^{-3/2} . Electronic Journal of Combinatorics , 13(1):Research Paper 63, 2006. 64 pp.
- 5[5] F. Bergeron, G. Labelle, and P. Leroux. Combinatorial species and tree-like structures . Number 67 in Encyclopedia of Mathematics and its Applications. Cambridge University Press, Cambridge, 1998.
- 6[6] J. Berstel. Sur les pôles et le quotient de Hadamard de séries N-rationnelles. Comptes-Rendus de l’Académie des Sciences , 272(Série A):1079–1081, 1971.
- 7[7] Henri Cartan. Elementary theory of analytic functions of one or several complex variables . Dover Publications, New York, 1995. Reprint of the 1973 edition, translated from the French version of 1961.
- 8[8] Jean-Pierre Dedieu. Points Fixes, Zéros et la Méthode de Newton . Number 54 in Mathématiques & Applications. Springer, 2006.