Effective Asymptotics of Linear Recurrences with Rational Coefficients

Xavier Gourdon^{∗}

             Xavier.Gourdon@inria.fr

Bruno Salvy^{∗}

                                          Bruno.Salvy@inria.fr

Abstract

We give algorithms to compute the asymptotic expansion of solutions of linear recurrences with rational coefficients and rational initial conditions in polynomial time in the order of the recurrence.

Introduction

We investigate sequences defined by a recurrence of the form

a_{k}u_{n+k} + a_{k−1}u_{n+k−1} + ⋯ + a_{0}u_{n} = 0,  (1)

where the coefficients a_{k} and the initial conditions belong to ℚ. This is probably the most simple type of recurrence one may encounter. Recurrences of this type are ubiquitous in many fields of applications (see [?] for numerous examples and references). Among the approximately 2300 sequences listed in Sloane’s book [?], one can estimate that about 13% are of this type [?]. In the rest of this paper “linear recurrence” always means “linear recurrence with rational coefficients” and we shall refer to u_{n} as a “linear recurrent sequence”.

Surprisingly, some problems related to linear recurrences remain open [?], and specially problems related to effectivity. Our aim in this paper is to describe an algorithm that computes an asymptotic expansion of a sequence obeying (??) in polynomial time in the order k of the recurrence. It is quite simple to find the asymptotic expansion of Fibonacci numbers with traditional tools, but these tools break down for illconditioned recurrences. The algorithm we describe works without any limitation on the value of k or those of the coefficients.

Given a recurrence such as (??), one usually computes its general term as a sum of exponential polynomials of the form ∑_{k=0}^{N} p_{k}n^{k}λ^{n}, where λ is an algebraic number. In Section 1 we shall describe an algorithm computing the coefficients p_{k} without factoring any polynomial. This general term does not solve the problem of asymptotic behaviour. To form a proper asymptotic expansion one has to order the moduli of the algebraic numbers λ occurring in the general terms. The problem which will occupy most of this paper is: How can one perform such an ordering exactly, i.e. we prove that the algorithms we propose work on the whole class of recurrences (??). We shall use techniques from computer algebra to free ourselves from problems of ill-conditioning related to the use of floating-point values. The result is an algorithm which, given a positive integer p and a linear recurrence (??) together with its initial conditions—or equivalently a rational function in ℚ(x) (see below)—outputs the p first exponential polynomials of the asymptotic expansion of the solution u_{n} of (??) as n tends to infinity.

As an example consider the recurrence

14448u_{n} + 15136u_{n+1} − 5761u_{n+2} − 990u_{n+3} + 315u_{n+4} = 0.  (2)

Its characteristic polynomial has four simple roots:

α_{1} = \frac{3}{2} + \sqrt{\frac{53}{12}}, α_{2} = \frac{3}{2} − \sqrt{\frac{53}{12}}, α_{3} = \sqrt{\frac{454}{35}}, α_{4} = −α_{3}.

∗Algorithms Project, INRIA, 78153 Le Chesnay Cedex, France.

To compute the asymptotic behaviour of u_{n}, it is necessary to decide which roots have the largest modulus. In this example the roots are given explicitly and this make it possible (but not trivial) to sort them by modulus. The techniques we develop do not require an explicit expression of the roots.

We describe two essentially different decision procedures to compute the asymptotic expansion. The first approach, purely algebraic, completely avoids factorizations. It is made expensive by the increase of degrees due to resultant computations. Currently this is the most natural computer algebra approach to the problem, and the most easily implemented. However, as soon as p ⩾ 2, its cost becomes potentially exponential in the order of the recurrence.

The second approach is perhaps more natural: it consists in using numerical approximations. In our example, the approximations

α_{1} = 3.601586702 …, α_{2} = −0.601586702 …, α_{3} = 3.601586951 …

make it possible to conclude. This example also shows that these numerical approximations may require a lot of precision. To avoid ill-conditioning, we use guaranteed numerical approximations and we show that this can be done in polynomial time in the order of the recurrence. Numerical approximations have long been banned from computer algebra because of the reluctance inherited from fixed precision routines. However, with the arbitrary precision provided by most computer algebra systems, we feel that it is time for floating point numbers to be rehabilitated in computer algebra.

The first step of the algorithm is to compute a suitable partial fraction decomposition of the generating function of u_{n}. Since factorization of polynomials is known to be polynomial-time but depressingly expensive, we shall avoid factorization and rely instead on a recent decomposition algorithm [?]. This is described in Section 1. In Section 2 and 3, we address the problem of comparing the moduli of the singularities (corresponding to the roots of the characteristic polynomial). As opposed to what happens usually in most algorithms involving algebraic numbers, we have to distinguish between roots of a given polynomial. A first method is described in Section 2, based on an algorithm [?] for comparing real algebraic numbers. At this stage, we can produce the desired asymptotic expansion. Section 3 describes a numerical alternative to the algebraic algorithms of Section 2, where we show how to get exact information from numerical values. We prove that this can be done with a cost that is lower than that of the algebraic method. In Section 4, we study optimizations that can be applied to subparts of our algorithm in practical cases. In particular we show there how rough numerical estimates can be used fruitfully. We conclude in Section 5 with a few examples taken from classical combinatorics.

1 Outline of the algorithm

1.1 Generating function

One can translate (??) into the rational generating function ∑ u_{n}z^{n} with O(k^{2}) rational operations: the generating function of the sequence (??) is

\begin{matrix}∑ & {k \atop i=0} & a_{i} ∑_{j=0}^{i−1} u_{j}z^{k−i+j} \\ ∑ & {k \atop i=0} & a_{k−i}z^{i} \end{matrix} .

The reciprocal conversion is also easy.

¿From the asymptotic point of view, the generating function approach enables us to use tools from complex analysis, like residue computation, which prove very effective. Because of the low cost of the conversion from a linear recurrence to the generating function, from now on we shall be concerned with rational functions only. Thus the input of our algorithm is a function f ∈ ℚ(z) regular at the origin, together with a positive integer p, and its output consists of the first p terms of the asymptotic expansion of [z^{n}]f(z)—the nth Taylor coefficient of f at the origin—as n tends to infinity.

1.2 Exact formula

In this section, we derive an exact formula for [z^{n}]f(z), based on a partial fraction decomposition that does not require factorization. This allows for both an efficient implementation and possible future extensions to rational functions with parameters or non-rational coefficients.

Algorithm 1 (Exact formula) Let f(z) = P (z)/Q(z) ∈ ℚ(z), with P and Q two relatively prime polynomials and deg(P) < deg(Q). To compute [z^{n}]f(z),

1. Compute Q = D_{1}D_{2}^{2} ⋯ D_{n}^{n} the square-free decomposition of Q. (Each D_{i} is a square-free polynomial.)

2. Using the decomposition algorithm [?], compute polynomials P_{i,j} ∈ ℚ[z] such that

f(z) = \frac{P (z)}{Q(z)} = ∑_{i=1}^{n} ∑_{j=1}^{i} ∑_{D_{i}(α)=0} \frac{P_{i,j}(α)}{(z − α)^{j}} ,  (3)

with deg(P_{i,j}) < deg(D_{i}). This requires only gcd computations.

3. For each (i, j) such that gcd(P_{i,j}, D_{i}) ≠ 1, write D_{i} = G_{i} H_{i}, where G_{i} = gcd(P_{i,j}, D_{i}) and rewrite all terms in (??) involving the polynomial D_{i} as

∑_{D_{i}(α)=0} \frac{P_{i,j}(α)}{(z − α)^{j}} = ∑_{G_{i}(α)=0} \frac{A_{i,j}(α)}{(z − α)^{j}} + ∑_{H_{i}(α)=0} \frac{B_{i,j}(α)}{(z − α)^{j}}

where A_{i,j} and B_{i,j} are obtained by Euclidian division of P_{i,j} by G_{i} and H_{i}. Repeat this process until all gcd’s are units. This gives a factorization of each D_{i} in the form D_{i} = D_{i,1} ⋯ D_{i,n_{i}} and the partial fraction decomposition has the form

f(z) = ∑_{i=1}^{n} ∑_{j=1}^{n_{i}} ∑_{k=1}^{i} ∑_{D_{i,j}(α)=0} \frac{P_{i,j,k}(α)}{(z − α)^{k}},  (4)

each P_{i,j,k} being a polynomial with rational coefficients, deg(P_{i,j,k}) < deg(D_{i,j}).

4. From this we get the value of [z^{n}]f(z):

[z^{n}]f(z) = ∑_{i=1}^{n} ∑_{j=1}^{n_{i}} ∑_{k=1}^{i} ∑_{D_{i,j}(α)=0} \frac{(n + 1) ⋯ (n + k − 1)}{(k − 1)!} · \frac{P_{i,j,k}(α)}{α^{n+k}}.  (5)

Step 1 is computed by repetitively differentiating and computing gcd’s. Step 3 guarantees that each P_{i,j,k}(α) in (??) is non zero. Step 4 is a consequence of the usual series expansion of (α − z)^{−k}. It is clear that Algorithm 1 runs in polynomial time. We do not worry about its complexity since the following algorithms are much more expensive. Example Let f(z) be the following input

\frac{z^{2} + 2}{(1 − z)^{2}(1 + z − 2z^{2} − z^{3} − 2z^{5})}.

In this case, the square-free decomposition of the denominator Q of f is given by D_{2} = 1 − z, D_{1} = 1+z −2z^{2} −z^{3} −2z^{5} (note that D_{1} is not irreducible). Step (2) of the algorithm then produces the following decomposition

f(z) = −\frac{1}{(1 − z)^{2}} − \frac{14/3}{1 − z} + ∑_{D_{2}(α)=0} \frac{h(α)}{α − z} ,

where h(α) = (11530α^{4} − 7778α^{3} + 11325α^{2} + 3889α − 8080)/93. Since h and D_{2} are relatively prime, Step 3 does not do anything, and then Step 4 produces the result:

[z^{n}]f(z) = −(n + 1) − \frac{14}{3} + ∑_{D_{2}(α)=0} h(α)α^{−n−1}.

To get an asymptotic expansion from this, we have to sort the moduli of the roots of D_{2} and compare them to 1. This is addressed in the following sections.

1.3 Specification of the algorithm

As already mentioned, the asymptotic expansion is governed by the successive “layers” of singularities of the rational function, sorted by increasing moduli. We shall need several ways to describe these moduli.

Notation For P (z) ∈ ℚ[z], we note ρ_{1}(P) < ρ_{2}(P) < ⋯ < ρ_{k}(P) the distinct moduli of the roots of P in increasing order. When there is no ambiguity, we simply denote these numbers ρ_{m}, 1 ⩽ m ⩽ k. We also note Z_{m}(P) the set of zeroes of P whose modulus is ρ_{m}(P), and Z_{m}^{+}(P) the subset of Z_{m}(P) whose elements have positive imaginary part.

We first state the form of the output on the example of f(z) from our previous section. The first four terms of the asymptotic expansion of [z^{n}]f(z) as given by our algorithm are

[(−1)^{n+1}h(−ρ_{1})] \frac{1}{ρ_{1}^{n+1}} + [h(ρ_{2}) + (−1)^{n+1}h(−ρ_{2})] \frac{1}{ρ_{2}^{n+1}} + [\frac{−(n + 1)}{ρ_{3}} − \frac{14}{3}] \frac{1}{ρ_{3}^{n+1}}

+[ ∑_{β∈Z_{3}^{+}(D_{2})} (h(β) + h(\overline{β})) cos[(n + 1) arg(β)]]\frac{1}{ρ_{4}^{n+1}} + o(\frac{1}{ρ_{4}^{n}} ) ,

together with the following information ρ_{1} = ρ_{1}(D_{2}), ρ_{2} = ρ_{2}(D_{2}), ρ_{3} = 1, ρ_{4} = ρ_{3}(D_{2}) and |Z_{3}^{+}(D_{2})| = 1. If requested, we can also give numerical approximations of the ρ_{i} and β.

All the features of the general case are present in this example. We now state precisely the specification of the algorithm, which we encourage the casual reader to skip. The input of our algorithm consists of f(z) = P (z)/Q(z) ∈ ℚ(z) and an integer p ⩾ 1. The output is the following asymptotic expansion of [z^{n}]f(z):

[z^{n}]f(z) = ∑_{ℓ=1}^{p} \frac{c_{ℓ}(n)}{ρ_{ℓ}^{n}} + o(\frac{1}{ρ_{p}^{n}} ) ,  (6)

where

c_{ℓ}(n) = H_{ℓ,0}(n) + (−1)^{n}H_{ℓ,1}(n) + ∑_{j ⩾ 2} ∑_{β∈Z_{n_{ℓ,j}}^{+} (D_{ℓ,j})} ∑_{k} \frac{H_{ℓ,j,k}(n)}{ρ_{ℓ}^{k}} cos[(n + k) arg(β)]

and explicit values are given for the coefficients of the polynomials H_{ℓ,0} and H_{ℓ,1} (in ℚ(ρ_{ℓ})[z]), D_{ℓ,j} (in ℚ[z]) and H_{ℓ,j,k} (in ℚ(β, \overline{β})[z]), as well as for the number of elements of the Z^{+} involved and a definition of ρ_{ℓ} either explicit or as ρ_{ℓ} = ρ_{ℓ′}(D_{i,j}) for some divisor of Q.

In practice, the program will be able to give numerical approximations of the moduli and the roots implicitly defined. Note that, because of the trigonometric functions involved in the coefficients, the expansion (??) is not of Poincaré type. Instead, we resort to the extended definition of Schmidt [?] (see also [?]), according to which an asymptotic expansion is a sum of the form

f(x) = ∑_{k=1}^{p} a_{k}(n) · f_{k}(n) + r(n)  (7)

where as n tends to infinity f_{k+1}(n) = o [f_{k}(n)], (1 ⩽ k ⩽ p − 1); r(n) = o [f_{p}(n)]; and a_{k}(n) are bounded functions of n that do not tend to zero.

1.4 Main algorithm

Starting from the partial fraction decomposition (??), we need to order the moduli of the roots of the D_{i,j} in (??) and find those roots that are purely real along with their signs. Our algorithm is based on the resolution of the two following computational problems.

Task 1 (Ordering the moduli) Given Q = ∏_{i,j} D_{i,j}^{i} a square-free decomposition of Q ∈ ℚ[z] and p a non-negative integer, compute for each (i, j) and for each k, 1 ⩽ k ⩽ p, the number of roots of D_{i,j} of modulus ρ_{k}(Q).

Task 2 (Real roots and their signs) Given P ∈ ℚ[z] a square-free polynomial and k a non-negative integer, compute the number p ∈ {0, 1} (resp. n ∈ {0, 1}) of positive (resp. negative) real roots of P of modulus ρ_{k}(P).

Most of the rest of this paper is devoted to algorithmic solutions to these tasks. Based on these, our main algorithm is as follows.

Algorithm 2 (Main algorithm) Let f(z) = P (z)/Q(z) ∈ ℚ(z) be a rational function, with deg(P) < deg(Q). Let p be a non-negative integer. To compute the p first terms of the asymptotic expansion of the coefficients of f(z),

1. Compute the partial fraction decomposition (??) by Algorithm ??.

2. Perform Task 1 to compute, for each (i, j) and for each ℓ, 1 ⩽ ℓ ⩽ p , the number m_{i,j,ℓ} of roots of D_{i,j} of modulus ρ_{ℓ} = ρ_{ℓ}(Q).

3. Select those terms in the expansion (??) for which |α| ∈ {ρ_{1}, …, ρ_{p}}, and rewrite (??) in the form

[z^{n}]f(z) = ∑_{ℓ=1}^{p} ∑_{{(i,j) \atop m_{i,j,ℓ}≠0}} ∑_{{D_{i,j} (α)=0 \atop |α|=ρ_{ℓ}}} ∑_{k=1}^{i} \binom{n + k − 1}{n}\frac{P_{i,j,k}(α)}{α^{n+k}} + o(\frac{1}{ρ_{p}^{n}} ) .  (8)

4. Perform Task 2 to compute, for each (i, j) and for each ℓ, 1 ⩽ ℓ ⩽ p, the number p_{i,j,ℓ} ∈ {0, 1} (resp. n_{i,j,ℓ}) of positive (resp. negative) real roots of D_{i,j} of modulus ρ_{ℓ}.

5. Rewrite relation (??) in the following form which is exactly the sought expansion (??):

∑_{ℓ=1}^{p}\{ ∑_{{(i,j) \atop m_{i,j,ℓ}≠0}}[p_{i,j,ℓ} ∑_{k=1}^{i} \binom{n + k − 1}{n}\frac{P_{i,j,k}(ρ_{ℓ})}{ρ_{ℓ}^{k}} + (−1)^{n}n_{i,j,ℓ} ∑_{k=1}^{i} \binom{n + k − 1}{n}\frac{P_{i,j,k}(−ρ_{ℓ})}{(−ρ_{ℓ})^{k}}

+ ∑_{{D_{i,j} (ρ_{ℓ}e^{iθ})=0 \atop sin θ>0}}(∑_{k=1}^{i} \binom{n + k − 1}{n} (P_{i,j,k}(ρ_{ℓ}e^{iθ}) + P_{i,j,k}(ρ_{ℓ}e^{−iθ})) \frac{cos[(n + k)θ]}{ρ_{ℓ}^{k}})]\}\frac{1}{ρ_{ℓ}^{n}} + o(\frac{1}{ρ_{p}^{n}} ) .

2 The algebraic method

To complete our main algorithm, there still remains to exhibit algorithms that perform Tasks 1 and 2. We describe in this section how this can be done purely algebraically. We rely principally on three tools:

(1) a method to order real algebraic numbers due to M. Coste and M.-F. Roy [?], based on Sturm sequences;

(2) a resultant computation that, given two polynomials P and Q produces a polynomial P ⊗ Q whose roots are the pairwise products of the roots of P and Q. In particular the smallest non-negative real root of P ⊗ P is the square of ρ_{1}(P) the smallest modulus of the roots of P;

(3) the Graeffe process: G_{k}(P) has for roots the kth power of the roots of P;

(4) the construction of a polynomial P_{k}(P) whose roots are the products α_{i_{1}} ⋯ α_{i_{k}} for i_{1} < ⋯ < i_{k}, the α_{i}’s being the roots of P.

Using (1) and (2) we can compare the smallest moduli ρ_{1}(P) and ρ_{1}(Q) of the roots of two polynomials P and Q. This will be done in Section ??. Using (2) and (3), we can produce the polynomials P_{k}(P) the modulus of the smallest root of which is |α_{1}| ⋯ |α_{k}|. This in turn enables us to compare any pair of moduli ρ_{i}(P) and ρ_{j}(Q) as will be shown in Section ??.

Note that other methods than Coste-Roy’s algorithm are known to compare real algebraic numbers (see, e.g. [?]). One of the reasons for our choice is that the complexity of Coste-Roy’s algorithm is known [?].

The polynomials mentioned above are computed by the formulas:

P ⊗ Q(y) = Resultant_{z}(P (z), z^{deg(Q)}Q(y/z)) , G_{k}(P)(z^{k}) = ∏_{j=0}^{k−1} P (e^{2ijπ/k}z),

∀k, 1 ⩽ k ⩽ n,  [P_{k}(P)]^{k} = \frac{∏_{i=0}^{⌊(k−1)/2⌋} P_{k−(2i+1)}(P) ⊗ G_{2i+1}(P)}{∏_{i=1}^{⌊k/2⌋} P_{k−2i}(P) ⊗ G_{2i}(P)},  (9)

where by convention we set P_{0}(P)(z) = z − 1. Apart from the last one, these polynomials are well known. That the last polynomial has the roots we expect is not difficult to check. All these polynomials have coefficients in the same field as P and Q.

2.1 Sorting the moduli

Given the polynomial Q, its factors D_{i,j} and an integer p, we need to determine the number of roots of each factor belonging to Z_{k}(Q), 1 ⩽ k ⩽ p. To simplify our description, we first concentrate on the case p = 1, corresponding to the first order estimate of the asymptotic expansion.

2.1.1 First order estimate

In this case (p = 1), our task can be performed in polynomial time in the degree of Q by Algorithm ?? below (which is an extension of an algorithm communicated to us by M.-F. Roy, taking into account multiplicities). We first describe an algorithm to compute the number of roots of smallest modulus of a polynomial.

Algorithm 3 (Number of roots of smallest modulus) Let P ∈ ℚ[z].

1. Compute P_{1}P_{2}^{2} ⋯ P_{n}^{n} the square-free decomposition of P ⊗ P.

2. Using Coste-Roy’s algorithm, find i_{0} such that P_{i_{0}} has the smallest non-negative real root.

3. Then |Z_{1}(P)| = i_{0}.

Proof. By construction, the smallest non-negative real root of P ⊗ P (z) is ρ_{1}^{2}(P). Moreover, its order of multiplicity is the number of roots of P of smallest modulus. Computing square-free decompositions in Step ?? ensures that only one of the polynomials P_{i} has the smallest non negative real root. □

Algorithm 4 (Smallest moduli comparison) Let P and Q ∈ ℚ[X].

1. Compute P_{i_{0}} and Q_{j_{0}} as in Algorithm ??. Their smallest non-negative real roots are ρ_{1}^{2}(P) and ρ_{1}^{2}(Q).

2. Applying Coste-Roy’s algorithm to P_{i_{0}} and Q_{j_{0}}, compare ρ_{1}^{2}(P) and ρ_{1}^{2}(Q).

3. The number of roots of P (resp. Q) of modulus ρ_{1}(P Q) = min(ρ_{1}(P), ρ_{1}(Q)) is given by i_{0} (resp. j_{0}) if ρ_{1}(P Q) is equal to ρ_{1}(P) (resp. ρ_{1}(Q)), and 0 otherwise.

Proof. This algorithm works for the same reason as Algorithm ??. □

Applying this algorithm to the polynomials D_{i,j} and Q gives the result we are after. Task 1 is therefore solved for p = 1. ¿From the complexity estimates in [?], it follows that the complexity of Algorithm ?? is O(n^{20}(n + log |P | + log |Q|)^{2}), where n = max(deg(P), deg(Q)) and |P | denotes the sum of the absolute values of the coefficients of the polynomial P.

2.1.2 Ordering the p smallest moduli

We now want to compute for each k, 1 ⩽ k ⩽ p and each (i, j), the number of roots of D_{i,j} of modulus ρ_{k}(Q). Although all the |α_{i}|^{2} are roots of P ⊗ P, Algorithm ?? does not generalize well because in general P ⊗ P has other non-negative real roots. We first give a generalization of Algorithm ??.

Algorithm 5 (Number of roots of a given modulus) Let P ∈ ℚ[z]. Given an integer q ⩾ 1, and m_{i} = |Z_{i}(P)|, for 1 ⩽ i ⩽ q, such that m_{1} + ⋯ + m_{q} < deg(P), to compute m_{q+1} = |Z_{q+1}(P)|, apply Algorithm ?? to the polynomial P̂ = P_{m_{1}+···+m_{q}+1}(P).

Proof. If P (z) = ∏_{i}(z − α_{i}), then by (??) we have P̂ = ∏_{i_{1}<···<i_{k+1}}(z − α_{i_{1}} ⋯ α_{i_{k+1}}), where k = m_{1} + m_{2} + ⋯ + m_{q}. Since |α_{1}| = ⋯ = |α_{m_{1}}| < |α_{m_{1}+1}| = ⋯ = |α_{m_{1}+m_{2}}| < ⋯ < |α_{k+1}| = ⋯ = |α_{k+m_{q+1}}| < ⋯, the roots of smallest modulus of P̂(z) are α_{k+j} ∏_{i=1}^{k} α_{i}, 1 ⩽ j ⩽ m_{q+1}. □ We can now give the generalization of Algorithm ??:

Algorithm 6 ((q + 1)st smallest moduli comparison) Let P and Q ∈ ℚ[z]. Given an integer q ⩾ 1 and m_{i} (resp. n_{i}) the number of roots of P (resp. Q) of modulus ρ_{i}(P Q) for 1 ⩽ i ⩽ q, such that (m_{1} + ⋯ + m_{q}) + (n_{1} + ⋯ + n_{q}) < deg(P) + deg(Q), to find the number m_{q+1} (resp. n_{q+1}) of roots of P (resp. Q) of modulus ρ_{q+1}(P Q),

– If m_{1} + ⋯ + m_{q} = deg(P) (resp. n_{1} + ⋯ + n_{q} = deg(Q)), then m_{q+1} (resp. n_{q+1}) is 0 and n_{q+1} (resp. m_{q+1}) is given by Algorithm ??.

– Otherwise, these values are obtained by applying Algorithm ?? to P̃ = P_{m_{1}+···+m_{q}+1}(P)⊗P_{n_{1}+···+n_{q}}(Q) and Q̃ = P_{n_{1}+···+n_{q}+1}(Q) ⊗ P_{m_{1}+···+m_{q}}(P).

Proof. The first part is obvious. Denote by α_{i} the roots of P and by β_{j} the roots of Q. Let M_{p} be the number of roots of P whose modulus is the smallest ρ_{i}(P) strictly greater than ρ_{q}(P Q) and define similarly M_{q}. The second part follows from noticing that the polynomial Q̃ has been built so that it has M_{q} roots of smallest modulus, namely β_{k} ∏_{i=1}^{m_{1}+···+m_{q}} α_{i} ∏_{j=1}^{n_{1}+···+n_{q}} β_{j}, n_{1} + ⋯ + n_{q} + 1 ⩽ k ⩽ n_{1} + ⋯ + n_{q} + M_{q}. Writing similarly the M_{p} roots of smallest modulus of P̃, one deduces the result. □

By induction on p ⩾ 1, using Algorithm ??, it is now easy to find for each (i, j) the number of roots of D_{i,j} of modulus ρ_{1}(Q), ρ_{2}(Q), … ρ_{p}(Q) with Q = ∏_{i,j} D_{i,j}^{i}. Task 1 is thus solved.

Because k = m_{1} + ⋯ + m_{q} can take any value between 1 and n, and since the degree of P_{m_{1}+···+m_{q}} is \binom{n}{k}, Algorithm ?? runs in exponential time as soon as p ⩾ 2.

2.2 Finding the real roots

We now attack Task 2: given an integer k, 1 ⩽ k ⩽ p, we want to find for each D_{i,j} the number and the sign of the real roots of D_{i,j} of modulus ρ_{k}(Q). The following algorithm solves this problem.

Algorithm 7 (Real roots and their sign) Let P ∈ ℚ[z] be a square-free polynomial. Given an integer q, the number m_{i} = |Z_{i}(P)| and the number n_{i} ∈ {0, 1} of real negative roots of P of modulus ρ_{i}(P) for 1 ⩽ i ⩽ q, with m_{1} + ⋯ + m_{q} < deg(P); to compute the number p_{q+1} ∈ {0, 1} (resp. n_{q+1}) of real positive (resp. negative) roots of P of modulus ρ_{q+1}(P),

1. Compute the polynomial P̂ = P_{m_{1}+···+m_{q}+1}(P) and m_{q+1} by Algorithm ??.

2. If m_{q+1} is odd, then P has exactly one real root of modulus ρ_{q+1}. To find its sign, compare the smallest positive real roots r of P̂(z) and ρ of P̂(−z) by Coste-Roy’s algorithm. If r > ρ or if P̂(z) has no positive real roots, then p_{q+1} ≡ n_{1} + ⋯ + n_{q} (mod 2) and n_{q+1} ≡ 1 + n_{1} + ⋯ + n_{q} (mod 2). Otherwise, either ρ > r or P̂(−z) has no positive real roots, and then p_{q+1} ≡ 1 + n_{1} + ⋯ + n_{q} (mod 2) and n_{q+1} ≡ n_{1} + ⋯ + n_{q} (mod 2).

3. If m_{q+1} is even, compute R(z^{2}) = gcd(P̂(z), P̂(−z)). If its degree is 0 then p_{q+1} = n_{q+1} = 0, otherwise use Coste-Roy’s algorithm to compare the smallest positive real roots r of R and ρ of P̂ ⊗ P̂. If R has no positive real roots, then p_{q+1} = n_{q+1} = 0. If r = ρ then p_{q+1} = n_{q+1} = 1. Otherwise we must have r > ρ and so p_{q+1} = n_{q+1} = 0.

Proof. First, note that P being square-free, p_{q+1} ∈ {0, 1} and n_{q+1} ∈ {0, 1}. Let P = ∏_{i}(z − α_{i}). The roots of P being either real or coming by pairs of conjugates, the number M = ∏_{i=1}^{m_{1}+···+m_{q}} α_{i} is real and its sign is the sign of (−1)^{n_{1}+···+n_{q}}. The polynomial P̂ has m_{q+1} roots of smallest modulus, namely

M · α_{i},  m_{1} + ⋯ + m_{q} + 1 ⩽ i ⩽ m_{1} + ⋯ + m_{q} + m_{q+1},  (10)

so that if m_{q+1} is odd, then P has exactly one real root of modulus ρ_{q+1}(P) and P̂(z) has only one real root of smallest modulus. Step 2 is now obvious.

When m_{q+1} is even, either p_{q+1} = n_{q+1} = 0 or p_{q+1} = n_{q+1} = 1 for conjugacy reasons. The roots of R(z^{2}) are the roots α of P̂ such that −α is also a root of P̂. Thus if deg(R) = 0 we cannot have p_{q+1} = n_{q+1} = 1 because of (??) . Otherwise, the smallest positive real root of R is the square of the smallest real root β of P̂ such that −β is also a root of P̂. The smallest positive real root of P̂ ⊗ P̂ being the square of the moduli of the roots (??), Step 3 is now clear. □

For the same reasons as Algorithm ??, Algorithm ?? runs in exponential time as soon as p ⩾ 2.

Tasks 1 and 2 have been solved, and we are now able to give the asymptotic expansion of the coefficients of a rational function by Algorithm ??.

3 The numerical method

In the last section, we solved Tasks 1 and 2 using only algebraic tools, which is currently the most natural solution to our problem from the computer algebra point of view. We shall now present an alternative method showing that numerical tools can be used reliably to perform our task in polynomial time. We only need to order the moduli of the roots of a polynomial and find which of them are real. Although there exists algorithms which achieve these tasks (for instance, we could use Graeffe’s method to approximate the moduli of the roots of Q), it is cheaper to find directly all the roots of Q with a sufficiently sharp bound on their errors. Our numerical method will depend on a complex root finding algorithm, that we first describe briefly.

3.1 A root finding algorithm

We want to find the complex roots of a polynomial with rational coefficients with arbitrary precision. Numerous algorithms exist to achieve this task, but only few of them are reliable. Newton’s method does not always converge; Traub and Jenkins’ method [?], usually used for root finding in computer algebra systems, converges theoretically but it turns out that precision control is badly handled in practical implementations. Besides, its complexity is not known to be polynomial. We present here a root finding algorithm for which the precision control has been carefully studied. In the following, |P | = ∑_{i} |a_{i}| denotes the norm of the polynomial P = a_{0} + ⋯ + a_{n}z^{n}. An immediate consequence of a theorem from [?] is the following result.

Proposition 1 (Pan) Let P ∈ ℂ[z] be a monic polynomial, n = deg(P) > 0. All the zeros of P can be computed with absolute error ϵ > 0 using O [n^{2} log n(n log(n) + log(|P |/ϵ))] arithmetic operations.

Unfortunately, the constant term in front of the time bound is very high and therefore the result seems to be only of theoretical importance. For instance, this algorithm relies on FFT techniques, which makes it efficient only for very large degrees.

3.2 Necessary precision

The reason why we can rely on numerical methods to solve our task is that two different roots of a polynomial with integer coefficients cannot be too close. The following result [?] makes this precise.

Proposition 2 (Mahler) Let P (z) = a_{0} + a_{1}z + ⋯ + a_{n}z^{n} = a_{n} ∏_{i=1}^{n}(z − α_{i}) be a polynomial of degree n > 0 with integer coefficients. Then

α_{i} ≠ α_{j} ⟹ |α_{i} − α_{j}| ⩾ \sqrt{3}n^{−(n+2)/2}M(P)^{1−n},

where M(P) = |a_{n}| ∏_{i=1}^{n} max(1, |α_{i}|).

¿From this we deduce the following theorem.

Theorem 1 Let P (z) be a polynomial with integer coefficients, n = deg(P) > 0 and α_{1}, …, α_{n} its roots. Define κ(P) to be the following quantity

κ(P) = \frac{\sqrt{3}}{2} [n(n + 1)/2]^{−[n(n+1)/4+1]} · M(P)^{−n(n^{2}+2n−1)/2},  (11)

then |α_{i}| ≠ |α_{j}| ⟹ ||α_{i}| − |α_{j}|| ⩾ κ(P) and |ℑ(α_{i})| is either 0 or larger than κ(P).

Proof. Let C be the leading coefficient of P. We first prove that the polynomial Q(z) = C^{n+1} ∏_{i ⩽ j}(z−α_{i}α_{j}) has integer coefficients. We can suppose that the polynomial P is primitive, i.e. the gcd of its coefficients is 1. The polynomial P ⊗ P has for roots α_{i}α_{j}, 1 ⩽ i, j ⩽ n, its coefficients are integers, and from classical results on the resultant algorithm, its leading coefficient is C^{2n}. Since the polynomial G_{2}(P)(z) = C^{2} ∏_{i}(z −α_{i}^{2}) has integer coefficients, is primitive and divides P ⊗ P, we deduce that the quotient C^{2(n−1)} ∏_{i≠j}(z − α_{i}α_{j}) = [C^{n−1} ∏_{i<j}(z − α_{i}α_{j})]^{2} has integer coefficients and therefore, so has its square root. Finally Q(z) is the product of this polynomial by G_{2}(P)(z) which implies it has integer coefficients.

Next, let x and y be two distinct roots of Q. Since M(Q) ⩽ M(P)^{n+1}, Mahler’s result applied to Q yields

|x − y| ⩾ γ = \sqrt{3} [n(n + 1)/2]^{−[n(n+1)/4+1]} M(P)^{(n+1)(1−n(n+1)/2)}.  (12)

If α_{i} and α_{j} are roots of P with distinct moduli, then |α_{i}|^{2} and |α_{j}|^{2} are two distinct roots of Q and we have ||α_{i}|^{2} − |α_{j}|^{2}| ⩾ γ, hence ||α_{i}| − |α_{j}|| ⩾ γ/(|α_{i}| + |α_{j}|). As |α_{i}| and |α_{j}| are smaller than M(P), we finally deduce

||α_{i}| − |α_{j}|| ⩾ \frac{γ}{2M(P)} = κ(P).

The last part of the theorem can be derived analogously from the inequality |α_{i}\overline{α_{i}} − α_{i}^{2}| ⩾ γ. □

A sharper lower bound on |ℑ(α_{i})| can be derived by considering only the polynomial P. We do not need this sharper bound since we need to compute the roots with an absolute error κ(P) to sort their moduli.

3.3 Numerical algorithm

Using these results, we now give an algorithm which performs reliably Tasks 1 and 2 by purely numerical methods.

Algorithm 8 (Numerical) Let Q = ∏_{i,j} D_{i,j}^{i} be a square-free decomposition of the polynomial Q. Our aim is to compute, for each (i, j) and for each q the number of roots of D_{i,j} of modulus ρ_{q}(Q) and the number of these that are real along with their signs.

1. For each (i, j), compute the number γ_{i,j} = inf_{(k,ℓ)≠(i,j)} γ(D_{i,j}D_{k,ℓ}) where γ(D_{i,j}D_{k,ℓ}) is defined by

γ(D_{i,j}D_{k,ℓ}) = \frac{\sqrt{3}}{2}[d(d + 1)/2]^{−[d(d+1)/4+1]} · |Q|^{−d(d^{2}+2d−1)/2},

with d = deg(D_{i,j}) + deg(D_{k,ℓ}). (Take D_{k,ℓ} = 1 if D_{i,j} is the only polynomial).

2. Using Pan’s Algorithm [?], compute for each (i, j) the roots of the polynomial D_{i,j} with an absolute error ϵ_{i,j} = γ_{i,j}/4.

3. Let α be a root of D_{i,j}, β a root of D_{k,ℓ}, α̂ and β̂ their approximations found at Step 2. If ||α̂| − |β̂|| < γ_{i,j}/2, then |α| = |β|, else the inequality between |α| and |β| is given by the inequality between |α̂| and |β̂|. This way, all the moduli of the roots of Q are sorted.

4. Let α be a root of some D_{i,j}. If its approximation α̂ satisfies |ℑ(α̂)| > γ_{i,j}/2, then α is not real. Otherwise α is real, and its sign is given by the sign of ℜ(α̂).

Proof. The validity of this algorithm results from Theorem ?? applied to each of the polynomials D_{i,j}D_{k,ℓ} and D_{i,j}, and from the inequalities:

(i, j) ≠ (k, ℓ) ⟹ M(D_{i,j}D_{k,ℓ}) ⩽ M(Q) ⩽ |Q|,
∀(i, j), M(D_{i,j}) ⩽ M(Q) ⩽ |Q|.

The inequality M(Q) ⩽ |Q| is due to Mahler [?]. In Step 4, the fact that the sign of α (when α is real) is the sign of ℜ(α̂) results from the inequality |α| ⩾ γ_{i,j}. (This latter inequality can be derived, for example, from the inequality |α| ⩾ 1/M(D_{i,j}) which is easily proved). □

Proposition 1 Algorithm ?? runs in time O [n^{5} log n log |Q|] where n is the degree of Q.

Proof. Apply Theorem ?? to each of the polynomials D_{i,j} with ϵ = ϵ_{i,j}. □

4 Optimizations

In the last two sections, we presented two methods that achieve Tasks 1 and 2. In practice, these two methods are awfully expensive. We present here another algorithm, which works on most of the rational functions, and which is much quicker. Another advantage of this new algorithm is that we can know whether it works or not. When it does not, then we can revert to one of the previous methods. This method is essentially numerical. We compute approximations of the roots of Q(z) using a root finding algorithm, with a relatively crude absolute error (compared to what it was in the previous section). In most cases though, everything can be deduced from these estimates.

In [?], A. Schönhage gave a root finding algorithm and demonstrated the following result.

Theorem 2 (Schönhage) Let P ∈ ℂ[z] be a monic polynomial, n = deg(P), and ϵ > 0. We can compute n complex numbers v_{1}, ⋯, v_{n} such that

|P − (z − v_{1}) ⋯ (z − v_{n})| < ϵ  (13)

within the time bound of O [(n^{3} log(n) + log(|P |/ϵ)n^{2}) log(n log(|P |/ϵ)) log log(n log(|P |/ϵ))] .

Although this bound seems slightly weaker than the previous one in Proposition ??, this one is in terms of bit complexity. Note that this algorithm does not approach directly the roots with an absolute precision ϵ. But from inequality (??) one can derive absolute error bounds on the roots of P. This algorithm was optimized by Gourdon [?] who implemented it in Maple; the program gives the right result in a reasonable time. We shall rely on this method to approximate roots of polynomials.

Let Q(z) = ∏_{i,j} D_{i,j}^{i} be a square-free decomposition of the polynomial Q. Using Schönhage’s algorithm, we compute for each (i, j) approximations α̂_{1}, …, α̂_{p} of the roots of D_{i,j} such that |D_{i,j}(z)−∏_{k}(z −v_{k})| < ϵ (we can assume D_{i,j} monic), with ϵ = 10^{−n}, where n = deg(Q). We have already seen that from this we can compute for each root α_{k} of D_{i,j} an absolute error bound τ_{k} > 0 such that |α̂_{k} − α_{k}| < τ_{k}. Suppose that the absolute bounds τ_{k} determine which roots are conjugates, which roots are real and what their sign is. To achieve Tasks 1 and 2, it then remains to compare the moduli of the non-conjugate roots. If again, the absolute error bounds τ_{k} make it possible to decide these comparisons, then we have finished. Otherwise, we have a certain number of couples of non-conjugates and distinct roots (α, β) of Q such that, if α̂ and β̂ are the approximations of α and β found and τ and τ′ the absolute error bounds found for these approximations, ||α̂| − |β̂|| < τ + τ′. We call these couples candidates. In this case, we use Algorithm ?? (see below) to test the equality of the moduli of the candidates. If all the candidates have the same modulus (this is often the case), then we have solved Tasks 1 and 2. Else, this algorithm failed and we use one of the previous methods discussed in Sections 2 and 3. The underlying idea is that it is very unlikely that two non-real roots of distinct moduli have the same argument.

Algorithm 9 (Equality of candidates) Let P = ∏_{i=1}^{n}(z − α_{i}) be a square-free polynomial ∈ ℚ[z]. We are given approximations α̂_{i}, absolute error bounds τ_{i} such that |α̂_{i} − α_{i}| < τ_{i}, and for each j, 1 ⩽ j ⩽ n, the number s_{j} of elements of the set Γ_{j} = {i, ||α̂_{i}| − |α̂_{j}|| < τ_{i} + τ_{j}}.

1. Compute the square-free decomposition P ⊗ P = P_{1}P_{2}^{2} ⋯ P_{r}^{r}.

2. By Sturm sequences [?], compute for each k the number m_{k} of non-negative real roots of P_{k}.

3. If (a) m_{1} + 2m_{2} + ⋯ + rm_{r} = n, (b) for all (i, j), either Γ_{i} ∩ Γ_{j} = ∅ or Γ_{i} = Γ_{j}, (c) for all ℓ, ℓm_{ℓ} = |∪_{s_{i}=ℓ}Γ_{i}|, then for all j, all the elements of Γ_{j} have the same modulus.

Proof. Since P ⊗ P = ∏_{i,j}(z − α_{i}α_{j}), the |α_{i}|^{2} are roots of it. If m_{1} + 2m_{2} + ⋯ + rm_{r} = n, then these are its only positive roots. The result is now obvious. □

5 Examples

Denumerants [?, p. 108]: the number of ways to make n francs with coins of 1, 2, 5, and 10 francs has for generating function

f(z) = \frac{1}{(1 − z)(1 − z^{2})(1 − z^{5})(1 − z^{10})}.

The ten singularities have the same modulus, but 1 being a singularity of order 4 is isolated in the decomposition (??) produced by Algorithm ??:

\frac{1/100}{(1 − z)^{4}} + \frac{7/100}{(1 − z)^{3}} + \frac{91/400}{(1 − z)^{2}} + \frac{21/50}{1 − z} + ∑_{β^{4}−β^{3}+β^{2}−β+1=0} \frac{4β^{3} − 3β^{2} + 2β − 6}{100(β − z)}

+ ∑_{α^{5}+2α^{4}+2α^{3}+2α^{2}+2α+1=0} \frac{17α^{4} + 9α^{3} + 17α^{2} + α + 1}{2000(z − α)^{2}} + \frac{α^{4} − 27α^{3} − 14α^{2} − 33α + 3}{500(z − α)}.

¿From this we deduce easily the first terms of the asymptotic expansion of [z^{n}]f(z):

\frac{n^{3}}{600} + \frac{9n^{2}}{200} + (\frac{421}{1200} + \frac{(−1)^{n}}{80} − ∑_{α∈Z_{1}^{+}(z^{4}+z^{3}+z^{2}+z+1)} \frac{α^{3} + \overline{α}^{3} + 2(α + \overline{α}) + 4}{250} cos[(n + 2) arg(α)]) n + O(1).

Sum of powers of Fibonacci numbers Since rational functions (when they are regular at infinity) are closed under Hadamard product, and the sum of a sequence is obtained by multiplying its generating function by 1/(1 − z), many operations that can be applied to a linear recurrent sequence yield another linear recurrent sequence. We consider here ∑_{k=1}^{n} F_{k}^{p}. It is not difficult (tedious, rather) to show that the generating function of F_{n}^{p} has the following expression for fixed p:

5^{\frac{1−p}{2}} ∑_{k=0}^{\frac{p−1}{2}} \frac{\binom{p}{k}F_{p−2k}z}{1 − (−1)^{k}L_{p−2k}z − z^{2}} , if p is odd;

5^{−p/2} [∑_{k=0}^{\frac{p}{2} −1} \binom{p}{k}(−1)^{k} \frac{2 − (−1)^{k}L_{p−2k}z}{1 − (−1)^{k}L_{p−2k}z + z^{2}} + \frac{(−1)^{p/2}\binom{p}{p/2}}{1 − (−1)^{p/2}z} ] , if p is even;

where L_{n} denote the Lucas numbers. From this we construct the generating function of the sum of the tenth powers of the Fibonacci numbers, which we give in compact form to our algorithm:

\frac{z^{9} − 87 z^{8} − 4047 z^{7} + 42186 z^{6} + 205690 z^{5} + 42186 z^{4} − 4047 z^{3} − 87 z^{2} + z}{z^{11}−89z^{10}−4895z^{9}+83215z^{8}+582505z^{7}−1514513z^{6}−1514513z^{5}+582505z^{4}+83215z^{3}−4895z^{2}−89z+1} .

The first stage of the algorithm produces the decomposition f(z) = ∑_{Q(α)=0} P (α)/(α − z), where Q is the denominator of f and P is the following polynomial:

P (z) =  − \frac{1088331670771}{46065550781250000} + \frac{147618897967}{128854687500000} z + \frac{25161090223051}{98535937500000} z^{2} + \frac{473498073791}{4692187500000} z^{3}

− \frac{48654728411}{541406250000} z^{4} + \frac{178616503}{8789062500} z^{5} − \frac{2060862361}{541406250000} z^{6} − \frac{4685959559}{4692187500000} z^{7}

+ \frac{383607377}{7579687500000} z^{8} + \frac{1645213621}{1675110937500000} z^{9} − \frac{496515521}{46065550781250000} z^{10}.

This decomposition implies that all the singularities are simple poles. The next stage of the algorithm is to determine the number of real and complex roots of each modulus for the first moduli of the roots. This is done by a numerical evaluation of the roots with error bound 10^{−4} which shows that all the roots are real, and yields their signs. For instance, the three first terms of the expansion are

[z^{n}]f(z) = \frac{P (ρ_{1})}{ρ_{1}^{n+1}} + (−1)^{n+1} \frac{P (−ρ_{2})}{ρ_{2}^{n+1}} + \frac{P (ρ_{3})}{ρ_{3}^{n+1}} + o (\frac{1}{ρ_{3}^{n}} ) ,

with ρ_{1} ≃ 0.00812, ρ_{2} ≃ 0.0212 and ρ_{3} ≃ 0.0753.

A large problem This combinatorial problem was considered in [?]. Starting with 1, we write down a sequence of words by counting the number of contiguous identical digits in the previous word. Thus the second word is 11 because there is one 1 in “1”. Then we have two 1s, hence the third word is 21, and so on. The first few words are: 1, 11, 21, 1211, 111221, 312211, 13112221, … We then consider the sequence of lengths of these words: 1, 2, 2, 4, 6, 6, 6, 8, …. What happens is that this sequence is rational of degree 72! From the table in [?, pp. 177–178], it is possible to compute this fraction by solving a linear system. The numerator is found to be

P (z) = 1 + z − z^{2} − z^{3} − z^{4} + z^{5} − 5z^{7} + 6z^{9} + 8z^{10} − 10z^{12} − 5z^{13} + z^{14} + 4z^{15} + z^{16} + 4z^{17} + 7z^{18}
+9z^{19} − 4z^{20} − 22z^{21} − 39z^{22} + 4z^{23} + 52z^{24} + 38z^{25} + 17z^{26} − 68z^{27} − 28z^{28} + 22z^{29} − 12z^{30}
+13z^{31} − 37z^{32} + 45z^{33} + 54z^{34} − 12z^{35} − 34z^{36} − 82z^{37} + 17z^{38} + 89z^{39} + 13z^{40} − 34z^{42} − 89z^{43}
+73z^{44} + z^{45} + 26z^{46} + 31z^{47} − 128z^{48} − 14z^{49} + 49z^{50} + 56z^{51} + 74z^{52} − 99z^{53} − 20z^{54} − 43z^{55}
+33z^{56} + 47z^{57} − 41z^{58} + 18z^{59} + 50z^{60} − 10z^{61} − 13z^{62} − 9z^{63} − 17z^{64} + 38z^{65} − 42z^{66} + 37z^{67}
+8z^{68} − 4z^{69} − 29z^{70} − 19z^{71} + 28z^{72} + 30z^{73} − 22z^{74} − 18z^{76} + 12z^{77},

and the denominator is

Q(z) = 1 − z − z^{2} − z^{3} + z^{4} + 3z^{5} − z^{7} − 2z^{8} + 3z^{13} + 3z^{14} − 2z^{15} − 5z^{16} − 8z^{17} + 7z^{18} + z^{19} + 8z^{20}
− 5z^{22} + 8z^{23} − 12z^{24} − 4z^{25} − z^{26} + 18z^{28} − 4z^{29} + 2z^{30} − 13z^{31} − 7z^{32} + 19z^{33} − 14z^{34} + 14z^{35}
− 6z^{36} − 4z^{37} + 13z^{38} − 9z^{39} − 7z^{40} + 4z^{41} − 8z^{42} + 7z^{43} + 5z^{44} + 7z^{45} − 12z^{46} + 17z^{47} − 22z^{48}
+ 8z^{49} − 7z^{50} + 16z^{51} − 6z^{52} − 7z^{53} − 6z^{54} + 3z^{55} + 19z^{56} − 5z^{57} − 5z^{58} − 14z^{59} + 8z^{60} + 2z^{61}
+ 7z^{62} − 5z^{63} + z^{64} − 8z^{65} + 14z^{66} − 11z^{67} + 16z^{68} − 18z^{69} + 9z^{70} − 9z^{71} + 6z^{72}.

One of the nice theorems in [?] states that this denominator is actually independent of the starting string, provided it different from “22”. Thus in the leading term of the asymptotic expansion, only the constant factor depends on the initial string.

Despite the large degree of this denominator, it turns out that the asymptotic expansion is not too difficult to find. For the sequence we consider, the decomposition of P/Q is

\frac{P (z)}{Q(z)} = R(z) + ∑_{Q(α)=0} \frac{F (α)}{z − α},

where R is a polynomial induced by the first terms, and F is a polynomial of degree 71 with 250-digit rational coefficients. This means that all the singularities are simple poles. If one is only interested in the first order estimate, it then remains to determine the number of roots of smallest modulus. As expected since the coefficients of the generating function are positive, one of these roots is a positive real number. Using the program of X. Gourdon based on A. Schönhage’s algorithm [?], we get that the two smallest moduli are approximately 0.767 and 0.861, with error bounds of the order 10^{−40}, which shows that the root of smallest modulus is alone (and therefore real). Thus, [z^{n}]f(z) ∼ F (ρ_{1})ρ_{1}^{−n−1}, ρ_{1} ≃ 0.767119 and F (ρ_{1}) ≃ 1.566. All the 72 moduli belong to the interval (0.767,1.151), showing the need for caution with numerical estimates.

Conclusion

Algorithm ?? should not be implemented blindly. Although its complexity is polynomial, the constant implied in the O() of Proposition ?? is very large. Thus in our last example above, the precision needed to compute the roots would be approximately 522000 digits. Instead, one should use this algorithm as an upper bound in an adaptative program based on a good numerical program such as [?] and Algorithm ??, increasing the precision if necessary.

Note also that we have never used the fact that in combinatorial contexts, the generating functions have only positive coefficients and thus by Pringsheim’s theorem (see [?]), one of their singularities of smallest modulus is real positive, the other ones having arguments commensurable with π. The computation of the first-order estimate could take advantage of this extra information.

This very simple problem of linear recurrences with rational coefficients is not yet completely solved. It would be useful in practice to have some control over the periodicities that may occur in the asymptotic expansions. This problem is exemplified with the following generating function:

\frac{z^{2} + 2z − 2}{(1 − 2z^{2})(1 − z)^{2}} ,

or equivalently u_{n} = 2u_{n−1} +u_{n−2} −4u_{n−3} +2u_{n−4}, u_{0} = u_{1} = 2, u_{2} = 5, u_{3} = 4. The first few terms are 2, 2, 5, 4, 9, 6, 15, 8, 25, 10, …. The first-order asymptotic approximation obtained from this generating function is (1 + cos nπ)2^{n−1} + o(2^{n}). What happens is that although valid for all positive n, this expression reduces to o(2^{n}) when n is odd. Better precision necessitates to look for further terms in the expansion. The ideal algorithm outputs a list of asymptotic expansions depending on arithmetic properties of n. Cancellation in this context is not a trivial problem. For instance, no algorithm is known to determine whether a linear recurrent sequence takes the value 0 for some index. It is known that when such a sequence cancels infinitely often, the indices where it cancels asymptotically form a finite union of arithmetic progressions that can be computed [?], but our problem is different since we are only concerned with indefinite cancellation of the dominant part, which does not satisfy a linear recurrence in general.

Acknowledgement

This work was supported in part by the ESPRIT III Basic Research Action Programme of the E.C. under contract ALCOM II (#7141).

References

[1] J. Berstel and M. Mignotte, Deux propriétés décidables des suites récurrentes linéaires, Bulletin de la Société Mathématique de France 104 (1976) 175–184.

[2] M. Bronstein and B. Salvy, Full partial fraction decomposition of rational functions, in: M. Bronstein, ed., ISSAC’93 (ACM Press, 1993) 157–160. Proceedings ISSAC’93, Kiev, Ukraine, July 1993.

[3] L. Cerlienco, M. Mignotte and F. Piras, Suites récurrentes linéaires. Propriétés algébriques et arithmétiques, L’Enseignement Mathématique XXXIII Fascicule 1-2 (1987) 67–108.

[4] L. Comtet, Advanced Combinatorics, (Reidel, Dordrecht, 1974).

[5] J.H. Conway, The weird and wonderful chemistry of audioactive decay, in: T.M. Cover and B. Gopinath, eds., Open Problems in Communication and Computation (Springer-Verlag, 1987) 173–188.

[6] M. Coste and M.-F. Roy, Thom’s lemma, the coding of real algebraic numbers and the topology of semi-algebraic sets, Journal of Symbolic Computation 5 (1988) 121–129.

[7] J. Dieudonné, Calcul Infinitésimal, (Hermann, Paris, 1968).

[8] X. Gourdon, Algorithmique du théorème fondamental de l’algèbre, Tech. rep. 1852, Institut National de Recherche en Informatique et en Automatique, (February 1993).

[9] M.A. Jenkins and J.F. Traub, A three-stage variable-shift iteration for polynomial zeros and its relation to generalized Rayleigh iteration, Numerische Mathematik 14 (1970) 252–263.

[10] K. Mahler, An application of Jensen’s formula to polynomials, Mathematika 7 (1960) 98–100.

[11] K. Mahler, An inequality for the discriminant of a polynomial, Michigan Mathematical Journal 11 (1964) 257– 262.

[12] V. Pan, Algebraic complexity of computing polynomial zeros, Computers and Mathematics with Applications 14 (1987) 285–304.

[13] S. Plouffe, Approximations de séries génératrices et quelques conjectures, Master’s thesis, Université du Québec à Montréal, Sept. 1992. Also available as Research Report 92-61, Laboratoire Bordelais de Recherche en Informatique, Bordeaux, France.

[14] R. Rioboo, Real algebraic closure of an ordered field. Implementation in Axiom, in: P.S. Wang, ed., Symbolic and Algebraic Computation (ACM Press, 1992) 130–137. Proceedings of ISSAC’92, Berkeley, July 1992.

[15] M.F. Roy and A. Szpirglas, Complexity of the computation with real algebraic numbers, Journal of Symbolic Computation 10 (1990) 39–51..

[16] H. Schmidt, Beiträge zu einer Theorie der allgemeinen asymptotischen Darstellungen, Mathematische Annalen 113 (1936) 629–656.

[17] A. Schönhage, The fundamental theorem of algebra in terms of computational complexity, Tech. rep., Mathematisches Institut der Universität Tübingen, 1982. Preliminary report.

[18] N.J.A. Sloane, A Handbook of Integer Sequences, (Academic Press, 1973).

[19] C. Sturm, Mémoire sur la résolution des équations numériques, Institut de France de Sciences Mathématiques et Physiques 6 (1835) 271–318.

[20] E.C. Titchmarsh, The Theory of Functions, second ed., (Oxford University Press, 1939).