Checking and Bounding the Solutions of Some Recurrence Relations^{⋆}

Roberto Bagnara and Alessandro Zaccagnini

Department of Mathematics, University of Parma, Italy

      bagnara@cs.unipr.it, alessandro.zaccagnini@unipr.it

Abstract. Recurrence relations play an important role in the field of complexity analysis since complexity measures can often be elegantly expressed by systems of such relations. This justifies the interest in automatic, precise, efficient and correct systems able to solve or to approximate the solution of systems of recurrence relations. Assume such a system is built. Since closed-form solutions for recurrences of even modest complexity can be so big and complex to be unmanageable, how can confidence on such a system be gained? How can we quickly validate or perhaps disprove its results? And, in those cases where the exact solution is too complex to be of practical use, how can we trade precision for efficiency by approximating them from below and from above? We also concern ourselves with a problem related to the handling of sets of solutions of recurrence relations: how can we confine by means of a lower bound and an upper bound a set of such solutions? We provide some solutions to these problems where we are careful to rely, whenever possible, on fast integer computations and/or conditions that are easy to check in a completely automatic way. The ongoing experimental evaluation of these ideas is giving very promising results, showing order-of-magnitude speedups over the more traditional methods.

1 Introduction

An important aspect of automatic complexity analysis [2, 6–8, 11, 13, 14, 18] is the possibility of automatically solving or approximating the solutions of systems of recurrence relations. In the literature, there are both techniques [4, 10] and software^{1} for solving recurrences, but they only deal with a rather restricted range of cases and sometimes assume the interaction with a human operator. On the other hand, existing, completely automatic systems such as [6, 7] only provide quite rough approximations of the solutions. An attempt at improving the state of the

⋆ The work of R. Bagnara has been partly supported by MURST project “Constraint Based Verification of Reactive Systems.” The work of A. Zaccagnini has been partly supported by MIUR project “Zeta and L Functions and Diophantine Problems in Number Theory.”

1 Such as, for instance, Maxima (http://sourceforge.net/projects/maxima/) and
MuPAD (http://www.mupad.de/).

art in this field is being made in the context of the PURRS project (Parma University’s Recurrence Relation Solver, see http://www.cs.unipr.it/purrs/). This project aims at the creation of a software library providing all the services needed for the efficient solution and approximation of recurrence relations, and is especially targeted at fully automatic complexity analysis. Quickly finding exact solutions and/or tight approximations, in closed form, to recurrence relations is certainly an important and challenging task. In this paper, however, we concentrate on some problems concerning the manipulation of such solutions/approximations, in a way that is totally independent from the techniques with which they have been obtained.

A system for automatic complexity analysis must provide an abstraction of (possibly infinite) sets of sequences of real numbers. A sequence of real numbers expresses the “cost” of one “process” in terms of some “input measure” expressed by means of a natural number. Notice that this definition is quite general: cost may be in terms of elapsed time, number of statements executed, memory used, number of packets exchanged over the network, …; a process may be a piece of software but also a communication protocol; the input measure can be any metric of the input of a program/procedure, or, say, the number of participants to some synchronization protocol. Because of well-known undecidability results, effective complexity analysis can only be based on approximations [5] and this is one of the reasons why we need to deal with sets of sequences. The domain of the analysis is thus given by ℘(ℝ_{∞}^{ℕ}): during the analysis, elements of ℘(ℝ_{∞}^{ℕ}) are generated by imposing a recurrence relation that the sequences must satisfy, by computing approximations of set union (this is required in order to approximate the complexity of conditionals), of addition (in order to approximate sequential composition), and so forth.

Assume that a recurrence relation solver is built and suppose, quite reasonably, that it is sophisticated and complex enough so that proving its overall correctness is beyond the state-of-the-art in program verification. Closed-form solutions for recurrences of modest complexity can require several pages to be printed out so that manually checking them is not an option. As we discovered the hard way, the naive approach of substituting the solution provided by the system into the original recurrence and simplifying the resulting expression is also impractical as it leads to serious efficiency problems. Thus, the first problem we attack in this paper is the following: how can confidence on such a solver be gained? To what extent can we validate its results? This is very important from a software engineering perspective, since it allows automatic regressions testing of the solver^{2} and is also an important tool for the development of critical systems (e.g., hard real-time systems in charge of crucial tasks) where the cost of an undetected bug in the recurrence relation solver could be too high. For both

2 There is an always increasing recognition of the importance of unit testing in software development. The Extreme Programming methodology (http://www. extremeprogramming.org/), for instance, pushes this to the point of requiring that the test is written before the unit.

these applications, it is also important that checking of the candidate solutions can be done quickly.

The problem of checking the correctness of the solution provided by a recurrence relation solver is easily reduced to the following problem: given f: ℕ → ℝ in some restricted, syntactically characterized class of functions, does f vanish identically? That is, do we have ∀n ∈ ℕ : f(n) = 0? The automatic solution of this problem presents two main difficulties: the first is that the computer algebra system may not be able to completely simplify expressions involving one integer variable; the second is that, even when this simplification succeeds, it may require too much computation resources (on the order of several minutes of CPU time even for the verification of modestly sized recurrences). The approach proposed in this paper is much more efficient.

The exact solutions of recurrence relations can also be very big and complex. This, besides the impact on the cost of verification, poses a more fundamental problem: the solution can be so complex so as to be unmanageable and thus useless to the client application. This suggests the need to find upper and lower approximations of the exact solutions by means of suitable, “simpler” functions, and this constitutes the second contribution of this paper (the recognition of the importance of approximating the solutions of recurrence equations dates back to Cohen and Katcoff [4]).

The last problem we tackle in this paper is tightly related to the previous one as it deals with the manipulation of lower and upper bounds in order to define approximations of ℘(ℝ_{∞}^{ℕ}). Formally, for some class of boundary functions 𝔅 ∈ ℘(ℝ_{∞}^{ℕ}), we restrict ourselves to the representation of the subset of ℘(ℝ_{∞}^{ℕ}) defined by

𝔉 =^{def} { F ∈ ℘(ℝ_{∞}^{ℕ}) | ∃l, u ∈ 𝔅 . ∀f ∈ F : l ≤ f ≤ u }.

In order to adopt this approach we need ways to compute the lower and upper bounds required for approximating the union of F_{1}, F_{2} ∈ 𝔉, that is, for computing a (possibly) small F ∈ 𝔉 such that F ⊇ F_{1} ∪ F_{2}. This problem reduces to the following: given b_{1}, b_{2} ∈ 𝔅, approximate, within the chosen class 𝔅 of boundary functions, max{b_{1}, b_{2}} from above and min{b_{1}, b_{2}} from below. (A definition of 𝔅 that suits the needs of the present paper is given in Section 5.)

Our long-term plan is, of course, to solve these three problems for a class of recurrence relations that is as wide as possible. In this paper we concentrate on the case when an exact formula for the solution of the recurrence exists, that is for linear recurrences of finite order with constant coefficients^{3} and also for those recurrences with variable coefficients that possess a hypergeometric solution [12]. For these classes of recurrences (which arise with very high frequency in, e.g., the complexity analysis of high-level languages), this paper presents a unified approach to the solution of the aforementioned problems. The solutions we propose share some aspects, such as the identification of dominant terms.

3 The classical case; see, e.g., [1] for a detailed description of both the theoretical aspects and the implementation issues.

What makes them particularly interesting is the fact that the problem of proving/disproving ∀n ∈ ℕ : f(n) = 0 and the problem of finding g such that ∀n ∈ ℕ : g(n) ≥ max{f_{1}(n), f_{2}(n)} or ∀n ∈ ℕ : g(n) ≤ min{f_{1}(n), f_{2}(n)} are reduced to testing a finite (and usually very small) set of conditions. Moreover, in several cases these conditions are simple comparisons between integers (i.e., they are very fast and exact) or amount to the numerical evaluation of suitable subexpressions of f, f_{1} and f_{2} at a small number of integers. We are not aware of other research works tackling these problems.

Notice that we are mainly interested in exact solutions and in the possibility of extracting from them upper and lower bounds that are valid for all positive integers: this is the reason why we do not deal with asymptotics here, though there are well-known techniques for obtaining them. See for example Wilf [17, Chapter 5], for a general description, or Gourdon and Salvy [9] for a case similar to the one considered here.

The paper is structured as follows: Section 2 recalls the classes of recurrence relations relevant to the present paper; Section 3 presents our proposal to validate candidate solutions to such recurrences; Sections 4 and 5 deal with the problems of computing upper and lower bounds, respectively; Section 6 concludes by presenting our first experimental results and discussing current and future work. The proofs of all the stated results are in Section A.

2 Preliminaries

Consider the linear recurrence with constant coefficients given by

x_{n} = a_{1}x_{n−1} + a_{2}x_{n−2} + ⋯ + a_{k}x_{n−k} + p(n),  (2.1)

for n ≥ k, where k is a fixed positive integer that is called order of the recurrence relation (2.1), and p is a function defined over the set of natural numbers. Here we assume tacitly that the coefficients a_{1}, …, a_{k} are real numbers, that a_{k} ≠ 0, and that the values x_{0}, …, x_{k−1} are given real numbers. Eq. (2.1) is associated to the homogeneous recurrence

\begin{cases}g_{n} = a_{1}g_{n−1} + a_{2}g_{n−2} + ⋯ + a_{k}g_{n−k}, \\ g_{0} = 1, \\ g_{n} = a_{1}g_{n−1} + a_{2}g_{n−2} + ⋯ + a_{n}g_{0}, & 1 ≤ n < k. \end{cases}  (2.2)

This can be solved by means of the characteristic equation

λ^{k} = a_{1}λ^{k−1} + ⋯ + a_{k−1}λ + a_{k}.  (2.3)

Let this polynomial equation have roots λ_{1} with multiplicity μ_{1}, …, λ_{r} with multiplicity μ_{r}, where the λ_{j}’s are distinct complex numbers, and the μ_{j}’s are positive integers. The general solution of (2.2) is

g_{n} = ∑_{j=1}^{r}(α_{j,0} + α_{j,1}n + ⋯ + α_{j,μ_{j}−1}n^{μ_{j}−1})λ_{j}^{n},  (2.4)

and the values of the coefficients α must be determined to fit the initial conditions in (2.2). Now the general solution of (2.1) is given by the formula, valid for n ≥ k,

x_{n} = ∑_{i=k}^{n} g_{n−i}p(i) + ∑_{i=0}^{k−1} g_{n−i}(x_{i} − ∑_{j=1}^{i} a_{j}x_{i−j}).  (2.5)

This is essentially equation (5) of [4]. It is worth stressing that even when it is possible to give a closed form for it, the right hand side of (2.5) can be, and in practice often is, quite complicated. This is partly due to the complex numbers arising from the solution of (2.3): these can be represented by algebraic expressions or in trigonometric form and, in the sequel, we will allow for both alternatives.

The other class of recurrences that we will deal with in this paper is the class of linear recurrences with variable coefficients of order 1, like

x_{n} = a(n)x_{n−1} + p(n),  (2.6)

where a: ℕ → ℝ\{0} is any sequence. The simplest such example is the factorial function where a(n) = n, x_{0} = 1 and p(n) = 0 identically. The general solution of (2.6), for n ≥ 1, is

x_{n} = (x_{0} + ∑_{k=1}^{n} p(k) ∏_{j=1}^{k} a(j)^{−1}) ∏_{k=1}^{n} a(k).  (2.7)

3 Validating Candidate Solutions

Assume an automatic recurrence relation solver produces the closed-form answer X_{n} for the recurrence (2.1) or for (2.6). A natural strategy for possibly validating or rejecting X_{n} as a solution is the following: first check whether X_{n} = x_{n} for n = 0, …, k − 1, i.e., that X_{n} satisfies the initial conditions. Secondly, define f: ℕ → ℝ as

f(n) =^{def} X_{n} − (a_{1}X_{n−1} + ⋯ + a_{k}X_{n−k}) − p(n),

for n ≥ k, or

f(n) =^{def} X_{n} − a(n)X_{n−1} − p(n),

for n ≥ 1, depending on the kind of recurrence we are dealing with, and try to verify or disprove that f vanishes identically.

Of course, it is the second condition that is hard to verify in practice^{4}, since it is possible that the simplification routines provided by the underlying computer

4 The verification of the initial conditions, though usually easier since it is essentially a numerical problem, is not necessarily trivial, especially in the case when the characteristic equation (2.3) has irrational roots, and it seems to be fairly hard in the case of irreducible equations of degree 3.

algebra package fail to simplify completely the expression f. Even in the fairly trivial case of the Fibonacci numbers, checking automatically that Binet’s closed formula actually holds is not really easy:^{5} in fact, it is equivalent to proving that, for all n ∈ ℕ, we have

ϕ^{n} − ϕ^{n−1} − ϕ^{n−2} − (−ϕ)^{−n} + (−ϕ)^{1−n} + (−ϕ)^{2−n} = 0,  (3.1)

where ϕ and −ϕ^{−1} are the two roots of the characteristic equation λ^{2} = λ + 1. It is clear that this really is the case if we collect a factor ϕ^{n} from the first three terms, and (−ϕ)^{−n} from the last three, and then recall that both ϕ and −ϕ^{−1} satisfy λ^{2} = λ + 1. Much more convoluted examples arise in practice, and we thus need an efficient procedure for deciding whether an expression which is not syntactically 0 vanishes identically as a function of the non negative integer n.

The objective is to find a finite (and possibly small) set of conditions that ensure that the expression f vanishes for all n ∈ ℕ (an infinite set of conditions). When possible, this task should be reduced to comparisons between integers or to the numerical evaluation of suitable subexpressions of f at a small number of integers. One possible general approach can be summarized as follows:

1. Deduce from either (2.1) or (2.6) the shape that any solution X_{n} must have, i.e., how it should look like from a purely syntactical point of view.

2. Check that X_{n} (and therefore f) has the appropriate shape.

3. Break down f into suitable subexpressions, and prove that each one vanishes identically, by evaluating it numerically, if necessary, at a small number of integers.

The next few subsections tackle this last problem for various classes of expressions, in increasing order of difficulty. The main idea is to find the dominant term in the expression f because this enables us to prove or disprove that they vanish identically, and also because we will use the same procedure for finding upper and lower bounds.

3.1 Polynomials and Exponentials

We begin with the “easy” case of the linear combination of products of polynomials and exponentials. Let ℂ[n] denote the set of polynomials in the variable n, with coefficients in ℂ; lead(p) denotes the leading coefficient of the polynomial p and deg(p) its degree, while coeff_{k}(p) is the coefficient of the monomial n^{k}. For completeness, we also state the classical result alluded to above.

Theorem 31 If q ∈ ℂ[n] then either q is the zero polynomial or it vanishes for at most deg(q) complex numbers.

5 It is, admittedly, a trivial example but we give it only to make our point clearer: more realistic examples, like the full exact solution of the recurrence x_{n} = x_{n−1} + x_{n−2} + x_{n−3} + n, say, would take far too much space.

We now introduce the class where, in view of (2.4) and (2.5), the solutions lie whenever p, the non-homogeneous part of the recurrence (2.1), is itself in the same class. It contains linear combinations of polynomials and complex exponentials, including polynomials themselves as a special case. The solutions we are concerned with are real-valued, but the definition is more general as there is no simplification in assuming that the elements of the class are also real valued.

Definition 32 (Exponential type functions) Let E denote the set

\{ f: ℕ → ℂ | \begin{matrix}∃r ∈ ℕ . ∀j = 1, …, r : ∃q_{j} ∈ ℂ[n] \\ . ∃α_{j} ∈ ℂ \ \{0\} . f(n) = ∑_{j=1}^{r} q_{j}(n)α_{j}^{n} \end{matrix} \}.

Notice that the complex numbers α_{j} need not be distinct, since this would make things harder in practice: if they are, then the representation in the definition is unique, apart from reordering, that is, in this case the integer r and the polynomials q_{j} are uniquely determined by f.

The following result generalizes Theorem 31 to the class E. It is necessary to assume that the values of the bases of the exponentials in Definition 32 be distinct: this will make the theorem difficult to use, and shows the need for some further manipulation of the expression.

Theorem 33 Let f ∈ E be given as in Definition 32 with the additional proviso that α_{j} ≠ α_{k} if j ≠ k. Then f vanishes for all n ∈ ℕ if, and only if, q_{j} is the constant polynomial 0 for all j = 1, …, r.

Note that the function f(n) = 1^{n} +(−1)^{n} vanishes for arbitrarily large values of n (but not identically for n ∈ ℕ), so that the exact analogue of the situation with polynomials (with respect to Theorem 31) does not hold. This particular f is the solution of the recurrence x(n) = x(n − 2), with x(0) = 2 and x(1) = 0, so this is not a contrived example, but one that may actually appear in practice.

Assuming f is given in such a way as to satisfy the hypothesis, the problem of deciding whether it vanishes reduces, thanks to Theorem 31, to the problem of deciding whether coeff_{k}(q_{j}) = 0 for each j = 1, …, r and each k = 0, …, D, where D is the largest degree of the q_{j}’s. We are thus left with the problem of ensuring that the expression f is of the form given in the statement of Theorem 33, and therefore that the values of α_{j} are distinct. Notice that this does not happen in (3.1). Moreover, a direct approach would require the ability to compare complex numbers given by arbitrary algebraic expressions. This problem can be avoided by rearranging f in a different way: instead of collecting terms with respect to the exponential functions, as in Definition 32, collect with respect to the monomials n^{d}, so that

f(n) = h_{0}(n) + nh_{1}(n) + ⋯ + n^{D}h_{D}(n),  (3.2)

where each h_{j} is just a linear combination of powers of the complex numbers α_{k}, and therefore has the general form

h(n) = β_{1}α_{1}^{n} + ⋯ + β_{r}α_{r}^{n}  (3.3)

for suitable coefficients β_{1}, …, β_{r}. This manipulation can be done quite efficiently and only requires comparisons between non negative integers.

The following lemma shows that it is sufficient to evaluate each h_{j} numerically at r_{j} consecutive integers, where r_{j} is the number of summands in h_{j}. If any of these values is different from 0, then f does not vanish for all integers, and therefore the suggested solution is not correct.

Lemma 34 Let α_{1}, …, α_{r} be distinct non-zero complex numbers, and β_{1}, …, β_{r} ∈ ℂ. If there exists an integer n_{0} such that the expression h defined by (3.3) vanishes for n = n_{0}, …, n_{0} + r − 1, then β_{1} = ⋯ = β_{r} = 0.

For example, in the case of (3.1), we have D = 0 so that we need to evaluate the left-hand side for n = 0, …, 5 to ensure that it actually vanishes for all n. This can be done automatically^{6}, so as to obtain a proof that Binet’s formula actually holds. Of course, if we are able to rewrite the left hand side as ϕ^{n}(1−ϕ^{−1} −ϕ^{−2})+(−ϕ)^{−n}(−1−ϕ+ϕ^{2}), then only two evaluations are needed.

If it could only be proved that h(n) is very small for n = n_{0}, …, n_{0} + r − 1, it is clear from the proof of Lemma 34 in the appendix that all the coefficients β are very small in their own turn.

3.2 Rational and Trigonometric Functions

It is possible to replace “polynomial” by “rational function” in the statement of Theorem 33: this is useful in the case when p, the non-homogeneous part of the recurrence (2.1), is a rational function itself.

Theorem 35 Let f: ℕ \ A → ℂ be defined by

f(n) =^{def} ∑_{j=1}^{r} q_{j}(n)α_{j}^{n},

where r ≥ 1 is an integer, q_{j} ∈ ℂ(n) for j = 1, …, r, α_{j} ∈ ℂ \ {0} for j = 1, …, r, α_{j} ≠ α_{k} if j ≠ k, and A denotes the (finite) set of all zeros of all denominators of the rational functions q_{j}. Then f vanishes for all n ∈ ℕ \ A if, and only if, q_{j} is the constant rational function 0 for all j = 1, …, r.

There is even more in Theorem 33 than meets the eye: In fact, recall the de Moivre formulæ, valid for θ ∈ ℝ:

cos(θz) = \frac{e^{iθz} + e^{−iθz}}{2} ,  sin(θz) = \frac{e^{iθz} − e^{−iθz}}{2i} .  (3.4)

They can be used to transform a function of the form

f(n) =^{def} ∑_{j=1}^{r} q_{j}(n)α_{j}^{n} + ∑_{j=r+1}^{r+s} q_{j}(n)α_{j}^{n} cos(θ_{j}n) + ∑_{j=r+s+1}^{r+s+t} q_{j}(n)α_{j}^{n} sin(ϕ_{j}n)

6 For example, as in our prototype system, using the GiNaC library (http://www.

  ginac.de/).

into a function of the shape of the statement of Definition 32, and vice versa. This can be useful when the characteristic equation (2.3) has complex roots and it is preferred to work only with real functions. Notice that even if the α’s, θ’s and ϕ’s are distinct, the resulting values of the new α’s need not be all different. In other words, the rewritten expression may not satisfy the hypotheses of Theorem 33. The simplest example of this phenomenon is probably the function f(n) = sin(πn), that vanishes for all n ∈ ℤ without being the zero function. Since we only need to know whether f vanishes for all nonnegative integers, we could plug (3.4) into f(n), and then rearrange the resulting expression according to the values of the new α_{j}’s. However, these substitutions do not change the degrees of the polynomials involved and the computations can be arranged as in Section 3.1, without any explicit substitutions, collecting with respect to monomials, and writing f(n) in the form (3.2) where D is the largest degree occurring in any term in f, and h_{d}(n) has the form

h_{d}(n) = ∑_{j=1}^{r} β_{d,j}α_{j}^{n} + ∑_{j=r+1}^{r+s} β_{d,j}α_{j}^{n} cos(θ_{j}n) + ∑_{j=r+s+1}^{r+s+t} β_{d,j}α_{j}^{n} sin(ϕ_{j}n)

for suitable complex numbers β_{d,j}. In other words, each h_{d} is (implicitly) a linear combination of at most r +2s+2t complex exponentials. Lemma 34 implies that g vanishes at all integers if h_{d}(n) vanishes at r + 2s + 2t consecutive integers, for d = 0, …, D.

3.3 Factorials

Formula (2.7) shows that natural candidates as solutions of linear recurrences with variable coefficients like (2.6) include the factorial function, or functions allied to it. In this case the results above can be somewhat generalized, since, for n → +∞, the function n! grows more rapidly than any function in the class E: This is an immediate consequence of Stirling’s formula (see [15, Section 4.42]), which we state in the slightly more general form that is needed in our applications. For any fixed a ∈ ℕ \ {0} and b ∈ ℤ, for n → ∞ we have

log(an + b)! = an log n + a(log a − 1)n + (b + \frac{1}{2}) log n
+ (b + \frac{1}{2}) log a + \frac{1}{2} log(2π) + O(n^{−1}).  (3.5)

Because of the functional equation satisfied by the factorial function (namely (n + 1)! = n!(n + 1)), in some cases it may be necessary to apply the rewrite rule below before it is possible to compare the order of magnitude of expressions containing factorials, exponentials and polynomials. Assume a ∈ ℕ \ {0} and b ∈ ℤ. Then

(an + b)! = \begin{cases}(an)!(an + 1) ⋯ (an + b), & if b ∈ ℕ \ \{0\}; \\ (an)!, & if b = 0; \\ \frac{(an)!}{(an)(an − 1) ⋯ (an + b + 1)}, & if b ∈ ℤ \ ℕ. \end{cases}

Linear Combinations of Factorials Consider a linear combination of expressions like (an + b)!, where a ∈ ℕ \ {0}, and b ∈ ℤ, with coefficients in the form of products of a rational function and an exponential. Once the rewrite rule has been applied on each factorial function in our expression, we clear denominators and are left with an expression of the form

f(n) = ∑_{j=1}^{r} h_{j}(n) (a_{j}n)!,  (3.6)

where h_{j} ∈ E for j = 1, …, r. We may assume, without any loss of generality, that a_{1} = a_{2} = ⋯ = a_{s} > a_{s+1} ≥ a_{s+2} ≥ ⋯ ≥ a_{r}. Once again, only comparisons between integers are required. Stirling’s Formula (3.5) implies that the dominant term in (3.6) is g(n) =^{def} (h_{1}(n) + ⋯ + h_{s}(n))(a_{1}n)!: therefore f(n) is 0 for all n ∈ ℕ only if g(n)/(a_{1}n)! vanishes for all n ∈ ℕ, and we are back to the case studied in Section 3.1, since g(n)/(a_{1}n)! belongs to E. Once we have verified that g vanishes identically, we can consider f − g, and repeat.

Products of Factorials More generally, we consider other combinations of the factorial function, beginning with products of several such functions: once again, formula (2.7) shows that this case is relevant for our problem.

Definition 36 (Factorial type functions) Let F denote the set

\{ f: ℕ \ A → ℂ | \begin{matrix}A is a finite set, \\ ∃r ∈ ℕ . ∀j = 1, …, r \\ : ∃a_{j} ∈ ℕ \ \{0\} . ∃b_{j} ∈ ℤ \\ . f(n) = ∏_{j=1}^{r}(a_{j}n + b_{j})! \end{matrix} \}.

In order to compare rates of growth of f and g in F and to fix notation, let a_{j} ∈ ℕ \ {0} and b_{j} ∈ ℤ for j = 1, …, r, and α_{j} ∈ ℕ \ {0} and β_{j} ∈ ℤ for j = 1, …, s. Let

f(n) =^{def} ∏_{j=1}^{r}(a_{j}n + b_{j})!  and  g(n) =^{def} ∏_{j=1}^{s}(α_{j}n + β_{j})!.

By (3.5) we have that

log ∏_{j=1}^{r}(a_{j}n + b_{j})!  is dominated by  ∑_{j=1}^{r} a_{j}n log n,

so that we can compare rates of growth of f and g by means of a comparison between the integers

∑_{j=1}^{r} a_{j}  and  ∑_{j=1}^{s} α_{j}.  (3.7)

If these integers are equal, we compare the second most important terms in (3.5), that is, we compare

∑_{j=1}^{r} a_{j}(log a_{j} − 1)  and  ∑_{j=1}^{s} α_{j}(log α_{j} − 1),

subject to the quantities in (3.7) being equal. Taking the latter into account and simplifying, this is equivalent to comparing the integers

∏_{j=1}^{r} a_{j}^{a_{j}}  and  ∏_{j=1}^{s} α_{j}^{α_{j}}.  (3.8)

It is actually possible that two different functions f and g ∈ F have the same values in both (3.7) and (3.8), as shown by the example

f(n) = (2n)!^{4}  and  g(n) = (4n)! n!^{4}.

If this happens, we have to look at the third highest term in (3.5), and compare

∑_{j=1}^{r} (b_{j} + \frac{1}{2})  and  ∑_{j=1}^{s} (β_{j} + \frac{1}{2}) ,

and, after multiplication by 2, we need another integer comparison between

r + 2 ∑_{j=1}^{r} b_{j}  and  s + 2 ∑_{j=1}^{s} β_{j}.  (3.9)

The functions f(n) = (2n)!^{7}(2n+1)! and g(n) = (4n)!^{2} n!^{8} have a common value in (3.7)–(3.9), and therefore the ratio f(n)/g(n) tends to a finite, nonzero limit as n → +∞. The other terms in (3.5) only contribute a bounded amount, and they can be used, in general, to compute the value of the limit of f(n)/g(n). Indeed, the limit has the value

{(2π)^{r−s} ∏_{j=1}^{r} a_{j}^{2b_{j}+1} ∏_{j=1}^{s} α_{j}^{−2β_{j}−1}}^{1/2}.  (3.10)

In this case, unless r = s, it is impossible to confine computations within the integers, since π ∉ ℚ.

Products of Factorials, Exponentials and Polynomials The situation where functions belonging to F are freely mixed with polynomials and exponentials is even more complicated: the simplest and most striking example of this fact is possibly given by the functions f(n) = πn(2n)!^{2} and g(n) = n!^{4} · 2^{4n}, whose ratio has the limit 1 as n → +∞. In this case, Stirling’s formula (3.5) is not sufficient to determine the limit of f(n) − g(n) as n → +∞. The question can be settled using the more precise formula where the O(n^{−1}) term is replaced by 1/(12n) + O(n^{−2}).

Let f(n) = q(n)α^{n}g(n), where q ∈ ℂ[n], α ∈ ℂ\{0} and g ∈ F has the shape of Definition 36. If lead(q) > 0, then for sufficiently large n we have that log f(n) is dominated by n log n ∑_{j=1}^{r} a_{j} + n(log α + ∑_{j=1}^{r} a_{j}(log a_{j} − 1)). It is therefore possible to repeat, to some extent, what we did in Section 3.3 and compare the order of magnitude of two such functions, though, in general, we may have to compare quantities involving logarithms of complex numbers.

The following examples show that all situations described earlier in this section may actually arise in practice, as solutions of suitable recurrences: The functions f(n) = n!^{2} and g(n) = (2n)! are the solutions of the recurrences x(n) = n^{2}x(n − 1) with x(0) = 1, and y(n) = 2n(2n − 1)y(n − 1) with y(1) = 1, respectively. They have the common value 2 in (3.7), but then they have the different values 1 and 4 respectively in (3.8). The functions f(n) = (2n)!^{4} and g(n) = (4n)! n!^{4} are the solutions of the recurrences x(n) = 16n^{4}x(n − 1) with x(0) = 1, and y(n) = 4n^{5}(4n − 1)(4n − 2)(4n − 3)y(n − 1) with y(0) = 1 respectively. They have the common value 8 in (3.7), and the common value 256 in (3.8), but then they have the different values 4 and 5 respectively in (3.9). Note that g(n)/f(n) ∼ cn^{1/2} as n → ∞, for a suitable constant c > 0. This is only one of a pattern of similar functions: take f(n) = (kn)!^{k^{k}} and g(n) = (k^{k}n)! n!^{k^{k}(k−1)}, where k ∈ ℕ \ {0}. These functions have the common value k^{k+1} in (3.7), and the common value k^{k^{k+1}} in (3.8), but for k > 1 they have the different values k^{k} and k^{k}(k − 1) + 1 respectively in (3.9). The functions f(n) = (2n)!^{7}(2n + 1)! and g(n) = (4n)!^{2} n!^{8} are the solutions of linear, first order recurrences with variable coefficients similar to the ones in the example above. They agree throughout (3.7)–(3.9). The functions f(n) = πn(2n)!^{2} and g(n) = n!^{4}2^{4n} are the solutions of the recurrences x(n) = 4n^{3}(2n − 1)^{2}x(n − 1)/(n − 1) with x(1) = π, and y(n) = 16n^{4}y(n − 1) with y(0) = 1 respectively. Finally, notice that the identity (n!)! = (n! − 1)! · n! is valid for all n ∈ ℕ.

3.4 Higher Transcendental Functions

We recall that a function f: ℂ → ℂ is entire if it is everywhere differentiable in the complex sense. Let A ⊂ ℂ be a (possibly infinite) set of complex numbers such that the topological closure of A is empty. In other words, let A intersected with any circle be a finite set. For each a ∈ A let m_{a} be a positive integer. A Theorem of Weierstrass (see [15, Section 8.11]) implies that there is an entire function f that has a zero of multiplicity m_{a} at the point a ∈ A, and does not vanish for z ∉ A. This means that it is very difficult to extend the results of the previous section to a larger class of functions, including higher transcendental functions like hypergeometric functions that may arise as solutions of recurrences of the form (2.6), or Euler’s Gamma: for example, the entire function 1/Γ (−z) vanishes if and only if z ∈ ℕ.

4 Upper Bounds and Maxima

Here we deal with the problem of determining upper bounds for the exact solution of a recurrence, or for max{f_{1}(n), f_{2}(n)}, where both f_{1} and f_{2} are exact solutions of some recurrences, or some upper bounds themselves. The problem is more difficult than validation, also because we are looking for rather different answers: we are interested both in tight bounds (and accept the fact that they may be very complicated) and in simple bounds (and accept that they may not be very sharp). In the latter case, we are ready to accept that the bounds may be valid only for all sufficiently large natural numbers. This raises the question of knowing exactly what “sufficiently large” means.

Another feature that makes the problem more difficult than validation is the possibility that, though we deal with real-valued solutions of recurrences like (2.1), the solutions themselves may contain complex numbers arising from the (possibly complex) roots of the characteristic equation (2.3). Terms involving complex roots may be rewritten as explained in Section 3.2 (actually, going the other way around) so that complex numbers disappear, but the problem is that they oscillate in sign.

The ideas explained in Section 3 can be also used to perform these tasks, and we briefly sketch the argument.

4.1 Upper Bounds for Functions in E

Let f ∈ E be written as in (3.2). For each h_{d}, defined as in (3.3), let h_{d}^{+}(n) = ∑_{j} β_{j}α_{j}^{n}, where the sum is restricted to the summands in h_{d} satisfying both β_{j} ∈ ℝ^{+} and α_{j} ∈ ℝ^{+}, and let h_{d}^{−}(n) = ∑_{j} β_{j}α_{j}^{n}, where the sum is restricted to the summands in h_{d} satisfying β_{j} ∈ ℝ^{−} and α_{j} ∈ ℝ^{+}. Observe that the function h_{d}^{+} is positive and that h_{d}^{−} is negative: in other words, they have a definite sign for all n ∈ ℕ. Finally, let h_{d}^{∗} = h_{d} − h_{d}^{+} − h_{d}^{−}: this part, if present, is oscillating. For ease of reference let

f^{+}(n) =^{def} ∑_{j=0}^{D} h_{j}^{+}(n)n^{d}, f^{−}(n) =^{def} ∑_{j=0}^{D} h_{j}^{−}(n)n^{d},  (4.1)

f^{∗}(n) =^{def} ∑_{j=0}^{D} h_{j}^{∗}(n)n^{d}.  (4.2)

An upper bound for f is therefore

f^{+}(n) + f^{−}(n) + |f^{∗}(n)|.  (4.3)

A simpler (but weaker) one is obtained by removing the term f^{−} from the expression above: this is the main reason for introducing the function f^{−}. An upper bound for |f^{∗}(n)| is obtained by replacing every term of the form βα^{n} in any h_{d}^{∗} by |β| · |α|^{n}: the resulting bound for the function in (4.3) belongs to the class E, but does not contain any oscillating term.

4.2 The Maximum of Two Functions in E

If both f_{1} and f_{2} are real-valued elements of E, we can determine an upper bound, U(n), say, for max{f_{1}(n), f_{2}(n)} as follows. Assume that f_{1} and f_{2} have the shape in Definition 32 for suitable values of the integers r_{1} and r_{2}, of the polynomials q_{1,j}, q_{2,j}, and of the complex numbers α_{1,j}, α_{2,j}, for j = 1, …, r_{1} and for j = 1, …, r_{2}, respectively. Assume further that the α’s are distinct, and that they are labeled so that for some s_{1} ≤ r_{1} we have

|α_{1,1}| = ⋯ = |α_{1,s_{1}}| > |α_{1,s_{1}+1}| ≥ ⋯ ≥ |α_{r_{1}}|,

and let d_{1,j} =^{def} deg(q_{1,j}), and similarly for f_{2}. Then f_{1}(n) is dominated by the term

g_{1}(n) =^{def} ∑_{\begin{matrix}j=1 \\ d_{1,j}=D_{1} \end{matrix}}^{s_{1}} β_{1,j}n^{D_{1}}α_{1,j}^{n}  (4.4)

where D_{1} =^{def} max{ d_{1,j} | j = 1, …, s_{1} }, and β_{1,j} = lead(q_{1,j}). In fact

f_{1}(n) = g_{1}(n) + o(n^{D_{1}}|α_{1,1}|^{n})  as n → +∞,  (4.5)

since every term in f_{1}(n) − g_{1}(n) has the form βn^{d}α^{n}, where β ∈ ℂ, and either |α| < |α_{1,1}|, or |α| = |α_{1,1}| and d < D_{1}. Let

g_{1}^{+}(n) =^{def} ∑_{\begin{matrix}j=1 \\ d_{1,j}=D_{1} \end{matrix}}^{s_{1}} b_{1,j}n^{D_{1}}|α_{1,j}|^{n}

where

b_{1,j} =^{def} \begin{cases}β_{1,j} & if α_{1,j} ∈ ℝ^{+}, \\ |β_{1,j}| & if α_{1,j} ∉ ℝ^{+}. \end{cases}

Then g_{1}(n) ≤ g_{1}^{+}(n) for all n ∈ ℕ. For a real number x, let

x^{+} =^{def} max{x, 0}, x^{−} =^{def} min{x, 0}.

Now let U(n) := 0 (where ‘:=’ denotes assignment) and repeat the following procedure.

– If |α_{1,1}| = |α_{2,1}| and D_{1} = D_{2}, let

\begin{cases}U(n) := U(n) + |α_{1,1}|^{n}n^{D_{1}} \\ · max\{(b_{1,1} + ⋯ + b_{1,s_{1}}), (b_{2,1} + ⋯ + b_{2,s_{2}})\}; \\ f_{1}(n) := f_{1}(n) − g_{1}(n); \\ f_{2}(n) := f_{2}(n) − g_{2}(n). \end{cases}

– If |α_{1,1}| < |α_{2,1}| or (|α_{1,1}| = |α_{2,1}| and D_{1} < D_{2}), let

\begin{cases}U(n) := U(n) + |α_{2,1}|^{n}n^{D_{2}}(b_{2,1} + ⋯ + b_{2,s_{2}})^{+}; \\ f_{2}(n) := f_{2}(n) − g_{2}(n). \end{cases}

– If |α_{2,1}| < |α_{1,1}| or (|α_{1,1}| = |α_{2,1}| and D_{2} < D_{1}), let

\begin{cases}U(n) := U(n) + |α_{1,1}|^{n}n^{D_{1}}(b_{1,1} + ⋯ + b_{1,s_{1}})^{+}; \\ f_{1}(n) := f_{1}(n) − g_{1}(n). \end{cases}

Notice that the definition of the coefficients b ensures that the bound obtained is sharper than the one given by replacing every β with its absolute value. In fact, take f(n) = −3^{n}: the latter alternative yields the upper bound f(n) ≤ 3^{n}, whereas our procedure gives f as an upper bound to itself. We need to take positive parts of real numbers to avoid errors in computing the upper bound of, say, f_{1}(n) = 3^{n} and f_{2}(n) = 3^{n} − 2^{n}.

4.3 The Maximum of Two Linear Combinations of Factorials

The discussion in Section 3.3 shows that essentially the same procedure can be used in the case of linear combinations of factorials, since any f of the shape in (3.6) can be broken down into terms of decreasing order of magnitude (as n → +∞), in the form of the product of a function from the set E and (an)!, for some positive integer a. More precisely, if f has the shape in (3.6) and we assume that h_{j} ∈ E for j = 1, …, r, and that a_{1} ≥ a_{2} ≥ ⋯ ≥ a_{s} > a_{s+1} ≥ ⋯ ≥ a_{r}, we use the procedure sketched in the previous section to get an upper bound H_{1}^{+} and a lower bound H_{1}^{−} for the function h_{1}+⋯+h_{s} ∈ E. Notice that in this case there will be no oscillating terms, making our analysis a little simpler. We then use the same ideas iteratively on the function f(n) − (h_{1}(n) + ⋯ + h_{s}(n))(an)!, that has again the shape (3.6) for suitable values of the parameters, adding together the upper bounds (resp. the lower bounds) found in the successive iterations.

4.4 Factorial Type Functions

Here we can use the procedure outlined in Section 3.3 to compare functions belonging to F: deciding which is the larger of f_{1} and f_{2} in F can be done via successive comparisons of the quantities in (3.7)–(3.9). If this is not enough to settle the question, it means that the ratio f_{1}/f_{2} has a non-zero, finite limit as n → +∞, and therefore the constant terms in (3.5) come into play. We then compute the expression in (3.10) and if this is larger than 1, then f_{1} is larger than f_{2} for sufficiently large n; if it is smaller than 1, then f_{2} is larger, and if it is exactly 1, then f_{1} and f_{2} have precisely the same behavior as n → ∞, and then we need a more precise form of Stirling’s Formula, as on a previous occasion in Section 3.3.

5 Lower Bounds and Minima

The strategy for obtaining lower bounds and minima is very similar to the one described in Section 4. Let f ∈ E be written as in (4.1)–(4.2). A lower bound for f is

f^{+}(n) + f^{−}(n) − |f^{∗}(n)|.

For minima, we define g^{−} exactly as g^{+}, but with the coefficients

b_{1,j} =^{def} \begin{cases}β_{1,j} & if α_{1,j} ∈ ℝ^{+}, \\ −|β_{1,j}| & if α_{1,j} ∉ ℝ^{+}. \end{cases}

Then g_{1}(n) ≥ g_{1}^{−}(n) for all n ∈ ℕ. Now let L(n) := 0, and repeat the following procedure.

– If |α_{1,1}| = |α_{2,1}| and D_{1} = D_{2}, let

\begin{cases}L(n) := L(n) + |α_{1,1}|^{n}n^{D_{1}} \\ · min\{(b_{1,1} + ⋯ + b_{1,s_{1}}), (b_{2,1} + ⋯ + b_{2,s_{2}})\}; \\ f_{1}(n) := f_{1}(n) − g_{1}(n); \\ f_{2}(n) := f_{2}(n) − g_{2}(n). \end{cases}

– If |α_{1,1}| < |α_{2,1}| or (|α_{1,1}| = |α_{2,1}| and D_{1} < D_{2}), let

\begin{cases}L(n) := L(n) + |α_{2,1}|^{n}n^{D_{2}}(b_{2,1} + ⋯ + b_{2,s_{2}})^{−}; \\ f_{2}(n) := f_{2}(n) − g_{2}(n). \end{cases}

– If |α_{2,1}| < |α_{1,1}| or (|α_{1,1}| = |α_{2,1}| and D_{2} < D_{1}), let

\begin{cases}L(n) := L(n) + |α_{1,1}|^{n}n^{D_{1}}(b_{1,1} + ⋯ + b_{1,s_{1}})^{−}; \\ f_{1}(n) := f_{1}(n) − g_{1}(n). \end{cases}

The procedures for obtaining the minimum of two linear combinations of factorials, or of two functions in the class F, are explained in Sections 4.3 and 4.4 respectively, since they are quite similar to those needed to compute the maximum.

We are now in position to specify the set 𝔅 of boundary functions that we consider in this paper: it is the set of linear combinations of elements of F with coefficients in

\{ f: ℕ → ℝ | \begin{matrix}∃r ∈ ℕ . ∀j = 1, …, r : ∃q_{j} ∈ ℝ[n] \\ . ∃α_{j} ∈ ℝ^{+} . f(n) = ∑_{j=1}^{r} q_{j}(n)α_{j}^{n} \end{matrix} \},

which is closed with respect to the algorithms we gave for computing lower and upper bounds.

6 Conclusion and Future Work

Recurrence relations play an important role in the field of complexity analysis since complexity measures of programs and procedures can often be elegantly expressed by systems of such relations. Their mechanical solution is, also for this reason, an interesting application of symbolic computation.

In this paper we have studied two problems related to the manipulation of (candidate) solutions of recurrence relations: proving or disproving that a closedform expression actually solves a recurrence relation and computing non-trivial lower and upper bounds to finite sets of such expressions. An important feature of the techniques we propose is that they rely as much as possible on integer computations and, more generally, can be efficiently implemented.

A prototype implementation is being developed in the context of the PURRS project (http://www.cs.unipr.it/purrs/), and the part concerning the correctness test for the solutions is almost finished and in current use. At the time of writing, regression testing of the PURRS solver is done on a set of 1286 (generalized) recurrences collected from various sources. Of those, 612 are linear recurrences of finite order with constant coefficients, 579 of which can be solved by the current version of the system. There are 542 recurrences in this class that can be verified by the method proposed in Section 3. On 2 of them the naive method does not terminate in reasonable time. The average speedup granted by the new method on the remaining 540 recurrences is 16, whereas the maximum speedup measured for the verification of one single recurrence is about 200 (this occurs for the verification of x_{n} = 2x_{n−1} + x_{n−2} + n^{3}). Slowdowns are also possible, but we have only observed them for recurrences where the cost of verification is very small, on the order of milliseconds. On this timescale, our measurements of CPU time are affected by very large relative errors, to the point that a small speedup can be mistaken for a slowdown. This situation arises for 39 recurrences. There are still 37 recurrences for which we are unable to check the solution provided by PURRS for correctness (for 33 of them the verification causes system thrashing using both methods, whereas for the remaining 4 the simplifier is not powerful enough): this is, of course, material for future work.

On the subject of future work, there is often the need, given two sets of sequences F_{1}, F_{2} ∈ 𝔉, to know whether one set majorizes the other, e.g., whether we have, for each f_{1} ∈ F_{1} and each f_{2} ∈ F_{2},

∀n ∈ ℕ : f_{1}(n) ≤ f_{2}(n).  (6.1)

This is required, for instance, in order for an optimizing program transformer to automatically decide whether a candidate transformation resulted into an actual improvement [16]. This problem reduces to the problem of approximating, as precisely as possible, the decision problem of the truth of b_{1} ≤ b_{2} given b_{1}, b_{2} ∈ 𝔅. (Notice the emphasis on approximation: as is customary in the field of static analysis, “yes”, “no” and “don’t know” are all sensible answers.)

Future work includes, in addition, the extension to more general recurrences than (2.1) and (2.6), also in those cases where it is possible to give lower and upper approximations for the solutions, but not a closed form. In this case, checking that approximations are correct involves solving the decision problem (6.1). We expect the work on such inequalities to be generally beneficial. For instance, it can be used to improve the procedure outlined in Section 4.2, which may produce an upper bound that is not the best possible (as the example with f_{1}(n) = n and f_{2}(n) = n^{2} shows).

Acknowledgments We wish to thank Tatiana Zolo and Enea Zaffanella for the help they gave us in the implementation of these ideas. We are also grateful to the reviewers of a previous version of this paper for their careful comments that helped us improve the paper.

References

1. Bagnara, R., Zaccagnini, A., and Zolo, T. The automatic solution of recurrence relations. I. Linear recurrences of finite order with constant coefficients. Quaderno 334, Dipartimento di Matematica, Università di Parma, Italy, 2003.

Available at http://www.cs.unipr.it/Publications/.

2. Blieberger, J., and Lieger, R. Worst-case space and time complexity of recursive procedures. Real-Time Systems 11, 2 (1996), 115–144.

3. Celis, P. Remark: Corrections and errors in John Ivie’s some MACSYMA programs for solving recurrence relations. ACM Transactions on Mathematical Software 10, 4 (1984), 477–478. See [10].

4. Cohen, J., and Katcoff, J. Symbolic solution of finite-difference equations. ACM Transactions on Mathematical Software 3, 3 (1977), 261–271.

5. Cousot, P., and Cousot, R. Abstract interpretation: A unified lattice model for static analysis of programs by construction or approximation of fixpoints. In Proceedings of the Fourth Annual ACM Symposium on Principles of Programming Languages (New York, 1977), ACM Press, pp. 238–252.

6. Debray, S., and Lin, N.-W. Cost analysis of logic programs. ACM Transactions on Programming Languages and Systems 15, 5 (1993), 826–875.

7. Debray, S. K., López-García, P., Hermenegildo, M. V., and Lin, N.-W. Lower bound cost estimation for logic programs. In Logic Programming: Proceedings of the 1997 International Symposium (Port Washington, NY, USA, 1997), J. Małuszyński, Ed., MIT Press Series in Logic Programming, The MIT Press, pp. 291–305.

8. Flajolet, P., Salvy, B., and Zimmermann, P. Automatic average-case analysis of algorithms. Theoretical Computer Science 79, 1 (1991), 37–109.

9. Gourdon, X., and Salvy, B. Effective asymptotics of linear recurrences with rational coefficients. Discrete Mathematics 153, 1-3 (1996), 145–163.

10. Ivie, J. Some MACSYMA programs for solving recurrence relations. ACM Transactions on Mathematical Software 4, 1 (1978), 24–33. See also [3].

11. Le Métayer, D. ACE: An automatic complexity evaluator. ACM Transactions on Programming Languages and Systems 10, 2 (1988), 248–266.

12. Petkovšek, M., Wilf, H. S., and Zeilberger, D. A = B. A. K. Peters, Natick, MA, 1996.

13. Rosendahl, M. Automatic complexity analysis. In Proceedings of the 4th International Conference on Functional Programming Languages and Computer Architecture (Imperial College, London, U.K., 1989), ACM Press, pp. 144–156.

14. Sands, D. Complexity analysis for a lazy higher-order language. In Proceedings of the 3rd European Symposium on Programming (Copenhagen, Denmark, 1990), vol. 432 of Lecture Notes in Computer Science, Springer-Verlag, Berlin, pp. 361– 376.

15. Titchmarsh, E. C. The Theory of Functions. Oxford University Press, London, UK, 1988.

16. Wegbreit, B. Goal-directed program transformation. In Proceedings of the 3rd ACM SIGACT-SIGPLAN Symposium on Principles of Programming Languages (Atlanta, Georgia, 1976), ACM Press, pp. 153–170.

17. Wilf, H. S. Generatingfunctionology. Academic Press, New York, 1994.

18. Zimmermann, P., and Zimmermann, W. Automatic complexity analysis of divide-and-conquer algorithms. In Proceedings of the 6th International Symposium on Computer and Information Sciences (Antalya, Turkey, 1991), M. Baray and B. Ozguc, Eds., Elsevier Science Publishers, pp. 395–404.

A Proofs

Proof of Theorem 33 on page 7. Assume that the α’s are labeled in such a way that for some s ∈ { 1, …, r } we have |α_{1}| = ⋯ = |α_{s}| > |α_{s+1}| ≥ ⋯ ≥ |α_{r}|, and let d_{j} =^{def} deg(q_{j}). For brevity we refer to (4.4) and (4.5) with the suffix 1 removed throughout. As in (4.5), for large n the function f is dominated by g(n), which is defined as in (4.4). If s = 1 then g(n) = lead(q_{1})n^{d_{1}}α_{1}^{n} and therefore f does not vanish for large enough n. If s > 1, then let

g(n) =^{def} ∑_{\begin{matrix}j=1 \\ d_{j}=D \end{matrix}}^{s} β_{j}n^{D}α_{j}^{n}

where D =^{def} max{ d_{j} | j = 1, …, s }, and β_{j} = lead(q_{j}). Let g(n) = n^{D}h(n). Again, we have that (4.5) holds, and this implies that the function f can not vanish for all large n. Indeed, we prove in Lemma A1 below that h(n)/α_{1}^{n} does not have the limit 0 as n → +∞, unless all β’s are 0 (which is against their definition). In other words, there is a positive constant c such that for some arbitrarily large values of n we have |h(n)/α_{1}^{n}| > c, and for these values of n the function f does not vanish. □

Proof of Lemma 34 on page 8. By hypothesis, we have that

(β_{1}α_{1}^{n_{0}})α_{1}^{n} + ⋯ + (β_{r}α_{r}^{n_{0}})α_{r}^{n}

vanishes for n = 0, …, r − 1. Therefore, changing if necessary the value of the coefficient β_{j} by the multiplicative factor α_{j}^{n_{0}} (which does not affect the result), we can assume that n_{0} = 0. But the Vandermonde determinant

det [\begin{matrix}1 & α_{1} & & α_{1}^{2} & … & α_{1}^{r−1} \\ . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . \\ 1 & α_{r} & & α_{r}^{2} & … & α_{r}^{r−1} \end{matrix}] = ∏_{i<j}(α_{j} − α_{i}) ≠ 0

since the α’s are distinct. Therefore the system in the unknowns β_{1}, …, β_{r}

\begin{cases}β_{1} & + & β_{2} & + & … & + & β_{r} & = & 0 \\ . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . \\ β_{1}α_{1}^{r−1} & + & β_{2}α_{2}^{r−1} & + & … & + & β_{r}α_{r}^{r−1} & = & 0 \end{cases}

has only the trivial solution β_{1} = ⋯ = β_{r} = 0. □

Lemma A1 Let ω_{1}, …, ω_{r} be distinct complex numbers such that |ω_{j}| = 1 for all j, and let α_{1}, …, α_{r} be non-zero complex numbers. Then either r = 1 and ω_{1} = 1, or the sequence

h(n) =^{def} α_{1}ω_{1}^{n} + ⋯ + α_{r}ω_{r}^{n}

does not have any limit as n → +∞.

Proof. If r = 1 and h(n) → L, then h(n) = L + o(1) as n → ∞, and |h(n) − h(n + 1)| = |α_{1}ω_{1}^{n} − α_{1}ω_{1}^{n+1}| = |α_{1}(1 − ω_{1})| = o(1): this implies ω_{1} = 1 since α_{1} ≠ 0.

Now assume that r ≥ 2, and that the statement has been proved for r − 1. If r > 1 and h(n) had the limit 0, then consider the sequence g(n) = h(n)ω_{r}^{−n} −α_{r}, which has the shape considered in the statement, with ω′_{j} = ω_{j}ω_{r}^{−1} for j = 1, …, r − 1, and r − 1 in place of r. Since ω′_{j} ≠ 1 for all j, the sequence g(n) has no limit by the induction hypothesis, and therefore h does not have the limit 0. If ω_{r} were 1, then g(n) = h(n) − α_{r} would have the limit L − α_{r}, and this is again impossible by the induction hypothesis.

We may assume from now on that h(n) has the limit L ≠ 0 and that ω_{j} ≠ 1 for all j. We choose ϵ = \frac{1}{3}|L| > 0; for n ≥ n_{0}(ϵ) we have h(n) = L + δ_{n}, where |δ_{n}| < ϵ. Set s_{n,k} = h(n) + ⋯ + h(n + k). On the one hand we have

|s_{n,k}| = |∑_{j=1}^{r} α_{j}ω_{j}^{n} \frac{1 − ω_{j}^{k+1}}{1 − ω_{j}} | ≤ ∑_{j=1}^{r} \frac{2|α_{j}|}{|1 − ω_{r}|},

which is a fixed, finite quantity, independent of n and k. On the other hand, for n ≥ n_{0} we have

|s_{n,k} − (k + 1)L| ≤ (k + 1)ϵ = \frac{1}{3}(k + 1)|L| ⟹ |s_{n,k}| ≥ \frac{2}{3}(k + 1)|L|,

by the triangle inequality, which is absurd for large k. □

Proof of Theorem 35 on page 8. Let q_{j}(z) = a_{j}(z)/b_{j}(z), where a_{j}, b_{j} ∈ ℂ[z] are polynomials, and b_{j} is not the 0 polynomial. Let g(z) =^{def} f(z)b_{1}(z) ⋯ b_{r}(z), so that g has the shape of the statement of Theorem 33. By hypothesis, g(n) vanishes for all large n, and the result follows. □