The Automatic Solution of Recurrence Relations

I. Linear Recurrences of Finite Order with Constant Coefficients^{∗}

           Roberto Bagnara (bagnara@cs.unipr.it)
    Alessandro Zaccagnini (alessandro.zaccagnini@unipr.it)
              Tatiana Zolo (zolo@cs.unipr.it)

Department of Mathematics University of Parma Italy

∗This work has been partly supported by the Italian-Spanish Azione Integrata IT 229 “Advanced Development Environments for Logic Programs”, by the University of Parma’s FIL scientific research project (ex 60%) “Pure and Applied Mathematics”, and by the MURST projects “Aggregate- and number-reasoning for computing: from decision algorithms to constraint programming with multisets, sets, and maps”, “Constraint Based Verification of Reactive Systems”, and “Zeta and L Functions and Diophantine Problems in Number Theory.”

Series Foreword

Complexity analysis aims at the derivation of upper and lower bounds to the complexity of algorithms, processes and data structures. The results of such analyses, which may be part or wholly automated, can be used, say, to decide whether mobile agents should be allowed to run in a given context, assist programmers in reasoning about the behavior of programs, guide applications of optimized program transformations, and discover efficiency bugs that are otherwise very difficult to detect.

Recurrence relations play an important role in the field of complexity analysis since complexity measures of, e.g., programs, can usually be very elegantly expressed by means of such relations. Therefore there is significant demand for efficient software systems capable of solving, with a high degree of precision, systems of recurrence relations. Moreover, to be really useful, the solvers need to be fully automatic, obtaining such solutions without human intervention.

Although the mathematical and computing literature describes several techniques and software for solving recurrences, these do not fulfill all of the above requirements; on the one hand, many of the methods assume interaction with a human operator and only deal with a rather restricted range of cases, while, on the other hand, the fully automated tools that are available only provide quite crude approximations of the exact solutions.

The PURRS project (Parma University’s Recurrence Relation Solver, see http://www.cs.unipr.it/purrs/) is working at improving the state of the art in this field. The aim of the project is to create a software library that provides the services needed for efficiently computing the solution or an approximation of the solution of a system of recurrence relations that arise in performing fully automated complexity analysis.

Finding exact solutions and/or tight approximations in closed form in an acceptable timescale for a large class of recurrence relations is a challenging task; such an objective requires solutions to a host of unresolved theoretical and practical problems. The current series of papers is devoted to the presentation of all the mathematics behind the PURRS project. As our aim is to provide both rigor and thoroughness, the series will not only include original work but also accounts of known results, describing any useful extensions or modifications.

Abstract

We describe algorithmic techniques for the efficient solution of a wide class of linear recurrences of finite order with constant coefficients. We give an outline of the underlying theory both from an abstract and a more concrete point of view, in an attempt to convey the general principles as clearly as possible yet providing a well marked path for the implementation. In particular, the presentation is thorough and reasonably self-contained, covering topics such as the automatic solution of polynomial equations and efficient, exact symbolic summation of some special functions.

Contents

1 Introduction 1

2 Preliminaries 2

3 Linear Recurrences of Finite Order with Constant Coefficients 3 3.1 The Basic Solution Technique . . . . . . . . . . . . . . . . . . 4 3.2 Abstract Setting of the Problem . . . . . . . . . . . . . . . . 5 3.3 Special Case: Recurrences of Order 1 . . . . . . . . . . . . . . 9

3.3.1 Polynomials and Exponentials . . . . . . . . . . . . . 9 3.4 Special Case: Recurrences of Order 2 . . . . . . . . . . . . . . 11

3.4.1 “Guessing” the Solution in Special Cases . . . . . . . 12 3.5 Higher-Order Linear Recurrences with Constant Coefficients . 13 3.6 Order Reduction . . . . . . . . . . . . . . . . . . . . . . . . . 13 3.7 Summary of the Method . . . . . . . . . . . . . . . . . . . . . 17

4 Systems of Linear Recurrences 17 4.1 Systems of Recurrences of Order 1 . . . . . . . . . . . . . . . 18 4.2 Systems of Recurrences of Higher Order . . . . . . . . . . . . 21

5 Related Work 24 5.1 Lueker, 1980 . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 5.2 Cohen & Katcoff, 1977 . . . . . . . . . . . . . . . . . . . . . . 26

6 Conclusion 28

A Solving Polynomial Equations 30 A.1 Detecting Multiple Roots . . . . . . . . . . . . . . . . . . . . 30 A.2 Detecting Rational Roots . . . . . . . . . . . . . . . . . . . . 31 A.3 Equations of Special Form . . . . . . . . . . . . . . . . . . . . 32 A.4 Solving the General Polynomial Equation of Degree 3 …. 32 A.5 Solving the General Polynomial Equation of Degree 4 …. 33 A.6 The Discriminant and the Resultant . . . . . . . . . . . . . . 34 A.7 Parametric Equations . . . . . . . . . . . . . . . . . . . . . . 35 A.8 Solving a Set of Linear Equations . . . . . . . . . . . . . . . . 36

B Exact Summation of Some Special Functions 36 B.1 Exponentials . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 B.2 Products of a Polynomial and an Exponential . . . . . . . . . 36 B.3 Polynomials . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 B.4 Some More Exact Formulæ . . . . . . . . . . . . . . . . . . . 40 B.5 How To Avoid the Use of Complex Numbers . . . . . . . . . 41 C Examples 42

C.1 Recurrences of Order 1 . . . . . . . . . . . . . . . . . . . . . . 42 C.2 Recurrences of Order 2 . . . . . . . . . . . . . . . . . . . . . . 44 C.3 Recurrences of Higher Order . . . . . . . . . . . . . . . . . . 48 C.4 Order Reduction . . . . . . . . . . . . . . . . . . . . . . . . . 49

1 Introduction

In this part of the series we deal only with linear recurrences of finite order with constant coefficients, in the cases when it is possible to give an exact closed formula for the solution without resorting to complex methods for summation in closed form.^{1}

Linear recurrences of finite order with constant coefficients (LRFOCCs, for short) arise frequently in complexity analysis. For instance, several common algorithms on lists and other kinds of sequences work by manipulating a few elements on, say, the front of the sequence, and invoke themselves recursively to get the job done on the entire list. In the analysis of functional and logic programs this is most often the case, and linear recurrences with constant coefficients of order 1 or 2 cover the vast majority of the recurrences an automatic complexity analysis tool has to solve. LRFOCCs of higher order are also important since, as we will see, they arise in the solution of systems of LRFOCCs of order 1.

Luckily enough, LRFOCCs are relatively easy to solve efficiently, especially if their order is small. In this paper we give an outline of the theoretical results underlying the automatic solution of LRFOCCs. This is presented both from an abstract and a more concrete point of view, in an attempt to convey the general principles as clearly as possible yet providing a well marked path for the implementation. In fact, the illustrated techniques (or variants of them) are used in the PURRS system (Parma University’s Recurrence Relation Solver, see http://www.cs.unipr.it/purrs/). Care is taken to provide a presentation that is thorough, reasonably self-contained, and can be followed both by mathematicians and computer scientists.

The plan of the paper is as follows: Section 2 introduces preliminary concepts and notations; Section 3 is the main one and presents solution techniques for linear recurrences with constant coefficients of order 1, 2 and higher; Section 4 shows how to solve systems of linear recurrences of order 1 with constant coefficients by reducing them to a set of autonomous recurrences (that is, independent of one another) of higher order; Section 5 discusses related work in this field while Section 6 concludes the main body of the paper. The appendix contains the tools required to turn the solution methods into actual algorithms. Procedures for the solution of polynomial equations of degree up to 4 are described in Appendix A. The exact symbolic summation of an important class of functions (linear combinations of products of polynomials and exponential functions) is treated in detail in Appendix B. Finally, several examples are collected in Appendix C so as to show the solution algorithms at work.

1Such methods are best viewed within the context of linear recurrences with variable coefficients [5]. Approximate summation will be the topic of another paper of this series [4]. We will discuss generalized recurrences like those arising from the analysis of “divide et impera” algorithms in [3].

2 Preliminaries

Throughout the paper we will use the following conventions: we denote by ℕ the set of nonnegative integers, so that 0 ∈ ℕ; we denote by C^{(k)} the k-th iterate of the operator C; the value of an empty sum is 0.

We assume that the reader is familiar with the basic properties of polynomials (roots and their multiplicity, division with quotient and remainder, greatest common divisor), and has a working knowledge of calculus and of linear algebra (eigenvalues, eigenvectors, characteristic equation of a matrix, vector spaces, dimension, kernel of a linear operator).

We let ℂ[n] denote the ring of polynomials in the indeterminate n with coefficients in ℂ. This is the smallest set (with respect to inclusion) that contains the set of complex numbers ℂ, the indeterminate n ∉ ℂ, and is closed with respect to addition and multiplication. The rings ℤ[n], ℚ[n] and ℝ[n] are defined similarly. Let p ∈ ℂ[n] \ {0}: we denote by deg(p) the degree of the polynomial p and by lead(p) its leading coefficient, that is, the coefficient of the monomial n^{deg(p)}. We also let coeff_{j}(p) denote the coefficient of the monomial of degree j of the polynomial p. With this notation, lead(p) = coeff_{deg(p)}(p). A more abstract definition of polynomial is given in Mignotte and Ştefănescu [19].

We are particularly interested in deriving, when possible, mathematical formulas that are in closed form, a concept that we define informally as follows: a mathematical expression denoting a function f: ℕ → ℂ is said to be in elementary closed form if it is, syntactically, a linear combination of a fixed number, r, say, of products of polynomials in ℂ[n] and (complex) exponentials. The number r must be an absolute constant, that is, it must be independent of all variables and parameters of the problem. As an example, the left hand side of the identity

∑_{k=0}^{n} \binom{n}{k} = 2^{n}

is not in closed form, while the right hand side is.

Generally speaking, we focus on effectively computable methods that can constitute the basis of practical, efficient algorithms. As a consequence, very often we refer to syntax instead of semantics. Moreover, in an attempt to keep the discussion more fluid, we do that implicitly. Thus, when we say “p is a polynomial” we actually mean “p is an expression that, syntactically, is a polynomial” or “p is, syntactically, of the form a_{0}λ^{k} + a_{1}λ^{k−1} + ⋯ + a_{k−1}λ + a_{k} where … ” Of course, sticking to pure syntactic equality would be too restricting. Hence, mathematical formulas are evaluated modulo a computable equality theory EQ. While it is reasonable to assume that this theory captures associativity and commutativity of ‘+’ and ‘·’ and other simple properties of elementary operations, EQ is otherwise left unspecified.

We will write a ≐ b when we want to emphasize that the expressions a and b are equal under EQ. It must be stressed that, while ‘≐’ is decidable, ordinary mathematical equality is not. For example, it is well-known that it is undecidable whether an expression involving polynomials and the sine function is equal to zero [7, 21, 22].

We prefer to use the term “recurrence relations,” but in the literature it is also possible to encounter the term difference equations. The latter is mainly used in a mathematical context, in order to give prominence to the tie with the differential equations. In fact, difference equations can be considered the discrete analogue of differential equations, and there exists a collection of mathematical tools, called difference calculus, which is quite similar to differential calculus. Throughout the paper we will use both terms.

3 Linear Recurrences of Finite Order with Constant Coefficients

A linear recurrence of finite order with constant coefficients is a recurrence of the form

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

for n ≥ k, where k is a fixed positive integer that is called order of the recurrence relation (3.1). Here and throughout we implicitly assume that a_{k} ≠ 0, for otherwise the order would be smaller. In Equation (3.1) the coefficients a_{j} are real or complex numbers. The function p is defined on ℕ ∩ [k, +∞), and is the non-homogeneous part of the recurrence (3.1). If p ≐ 0 (that is, p vanishes identically from a syntactical point of view), we say that the recurrence (3.1) is homogeneous. Every recurrence of the form (3.1) is given with a set of k initial values x_{0}, x_{1}, …, x_{k−1}, which we assume to be known.

In the next section we present the general solution technique for this class of recurrences, whereas Section 3.2 presents the problem and its solution in a more abstract setting. The problem is then reconsidered under a more concrete point of view. Sections 3.3 and 3.4 are devoted to the “easy” cases of linear recurrences with constant coefficients of order 1 and 2, respectively; the description of tailor-made, more efficient solution methods for these subclasses is justified by the observation that these cases appear to be the most frequent ones in the context of complexity analysis. Section 3.5 deals with the general case of linear recurrences with constant coefficients of any finite order, while in Section 3.6 we show how to solve special recurrences of higher order in a more efficient way.

3.1 The Basic Solution Technique

We associate to any recurrence like (3.1) the homogeneous recurrence

g_{n} = a_{1}g_{n−1} + a_{2}g_{n−2} + ⋯ + a_{k}g_{n−k}.  (3.2)

We first solve this recurrence by means of the characteristic equation that, by definition, is the polynomial equation

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

We use throughout the notation g for the homogeneous recurrence associated to a general recurrence of finite order, and λ for the complex variable in the associated characteristic equation. In general, the characteristic equation (3.3) has k complex roots λ_{1}, …, λ_{k}. If these roots are all distinct (we also say that they are simple), the general solution of the homogeneous recurrence (3.2) is given by

g_{n} = α_{1}λ_{1}^{n} + ⋯ + α_{k}λ_{k}^{n},  (3.4)

where α_{1}, …, α_{k} are complex numbers. It is quite easy to see that α_{i}λ_{i}^{n} is a solution of (3.2), for i = 1, …, k. This is an immediate consequence of the fact that each λ_{i} is a root of (3.3). Furthermore, it is also easy to see that the sum of any two solutions of (3.2) is again a solution, so that the right hand side of (3.4) gives indeed solutions of the homogeneous recurrence (3.2).

If the characteristic equation (3.3) has multiple (or repeated) roots, we collect them: say they are λ_{1} with multiplicity μ_{1}, …, λ_{r} with multiplicity μ_{r}, where the λ_{j}’s are different complex numbers, and the μ_{j}’s are positive integers. In this case the general solution is

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

Of course, (3.5) contains (3.4) as the special case where r = k and μ_{j} = 1 for j = 1, …, r. The proof of this result can be found in Kelley & Peterson [14, Theorem 3.7]. Note that assuming that a_{k} ≠ 0 is tantamount to saying that λ = 0 is not a solution of (3.3), but it is clear from either (3.4) or (3.5) that this does not really matter.

The results above mean that any recurrence relation of the type (3.2) has a solution of the type (3.4) or (3.5), for suitable values of the coefficients α, which are determined by means of the initial values g_{0}, …, g_{k−1}.

Once the general solution of the homogeneous equation (3.2) has been found, we have to solve the complete equation (3.1): by linearity, it is sufficient to find any solution of the complete equation (we call this a particular solution), sum the general solution of the homogeneous equation and impose that the initial conditions hold. This is the so-called superposition principle (see Kelley & Peterson [14, Theorem 3.3]): we show why it is true by means of an analogy. In order to know a straight line in the usual plane, it is sufficient to know which straight line through the origin it is parallel to (and this is uniquely determined) and any point on the original line, it does not matter which one. In our case we need to know the complete solution of the homogeneous equation (3.2) (which plays the role of the line through the origin), and any solution of the complete, non-homogeneous equation (3.1).

3.2 Abstract Setting of the Problem

We now restate the problem we are interested in in a more abstract manner, and fully characterize the set of non-homogeneous parts of a recurrence of the type (3.1) that possess a solution whose elementary closed form can be computed without transcendental summation methods that will be the topic of [5]. Some remarks on how to actually and efficiently compute these closed forms can be found in Appendix B.

In the language of linear algebra, we are trying to solve the equation Ax = p, where A is a linear operator, so we first compute the kernel of A which we denote by ker(A) (that is, we solve the equation Ax = 0). The theory of linear operators guarantees that the general solution of Ax = p is given by the sum of an element of ker(A) and any solution of the complete equation. Furthermore, the dimension of ker(A) is precisely the same as the order of the recurrence, so that (3.5) gives all solutions of the homogeneous recurrence (3.2): in fact, there are exactly k unknown coefficients, and the functions n^{i} · λ_{j}^{n}, for j = 1, …, r and i = 0, …, μ_{j} − 1, are linearly independent over ℂ provided that the values of the λ_{j}’s are different. This means that the dimension of the vector space generated by these functions is precisely k. More concretely, we are talking about the space X of complex-valued sequences defined over ℕ, which, endowed with the obvious operations, can be considered a vector space over ℂ: if (x_{n})_{n∈ℕ}, (y_{n})_{n∈ℕ} ∈ X, and λ, μ ∈ ℂ, then

λ · (x_{n})_{n∈ℕ} + μ · (y_{n})_{n∈ℕ} =^{def} (λx_{n} + μy_{n})_{n∈ℕ} ∈ X .

For a recurrence of order k as in (3.1), the operator A is defined by

A((x_{n})_{n∈ℕ}) =^{def} (x_{n} − (a_{1}x_{n−1} + a_{2}x_{n−2} + ⋯ + a_{k}x_{n−k}))_{n∈ℕ},

with the convention that, if 0 ≤ n < k, the above definition has to be replaced by 0. It is easily seen that A is linear, that is

A(λ · (x_{n})_{n∈ℕ} + μ · (y_{n})_{n∈ℕ}) = λ · A((x_{n})_{n∈ℕ}) + μ · A((y_{n})_{n∈ℕ}),

for any (x_{n})_{n∈ℕ}, (y_{n})_{n∈ℕ} ∈ X, and λ, μ ∈ ℂ. Thus ker(A) is exactly the set of solutions of (3.2). The fact that A is linear is the property that ensures that the superposition principle holds.

Summing up, the solution of a homogeneous linear recurrence of finite order with constant coefficient is completely straightforward from a theoretical point of view, the main difficulty being in the solution of the characteristic equation (3.3). On the contrary, finding the particular solution of the general equation (3.1) is non-trivial, and an elementary closed form can only be obtained in comparatively few cases: in general we have recourse to summation formulæ and other approximations that will be explained in [4].

The ideas expressed above are now rigorously formalized.

Definition 3.1 (Shift and identity operators.) Let X denote the complex vector space of all sequences defined over the set ℕ of natural numbers. We denote by (x_{n})_{n∈ℕ} its generic element (sometimes dropping the subscript for brevity). We consider the shift operator E: X → X defined by

E((x_{n})_{n∈ℕ}) =^{def} (x_{n+1})_{n∈ℕ}.

We also consider the identity operator I: X → X (which does nothing).

We need linear combinations of iterates of these two operators: we define Δ(α) =^{def} E − αI, and introduce the family

F =^{def} { Δ(α_{1}) ∘ ⋯ ∘ Δ(α_{k}) | k ∈ ℕ \ {0}, α_{1}, …, α_{k} ∈ ℂ }.

We remark that the operators Δ commute. With this notation, we first rewrite (3.1) in the form

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

and then

(E^{(k)} − a_{1}E^{(k−1)} − ⋯ − a_{k−1}E − a_{k}I)(x_{n}) = p(n + k).

If p vanishes identically, the formal analogy with the characteristic equation (3.3) is clear. It is also clear that we would like to write the linear operator on the left as a composition of Δ(α) for suitable complex numbers α. This corresponds to solving the characteristic equation, for once we know that its roots are λ_{1} with multiplicity μ_{1}, …, λ_{r} with multiplicity μ_{r}, we have

E^{(k)} − a_{1}E^{(k−1)} − ⋯ − a_{k−1}E − a_{k}I = Δ(λ_{1})^{(μ_{1})} ∘ ⋯ ∘ Δ(λ_{r})^{(μ_{r})}. (3.6)

Here we exploit the commutativity of these operators. In order to solve homogeneous recurrences we need to understand the kernel of these operators.

Lemma 3.2 For any α ∈ ℂ and any integer d ∈ ℕ we have the set identity

ker(Δ(α)^{(d+1)}) = { (q(n)α^{n})_{n∈ℕ} | q ∈ ℂ[n], q = 0 or deg(q) ≤ d }.

Proof. First, assume that q ∈ ℂ[n] is a polynomial of degree ≤ d, and note that

Δ(α)(q(n)α^{n}) = α(q(n + 1) − q(n))α^{n}.

Since deg(q(n + 1) − q(n)) = d − 1, we see that Δ(α)^{(d)}(q(n)α^{n}) = kα^{n} for some constant k, and then it is easy to see that Δ(α)(kα^{n}) vanishes identically.

In the other direction, we assume that α ≠ 0, for otherwise the result is trivial. Now let (x_{n})_{n∈ℕ} ∈ ker(Δ(α)^{(d+1)}). If d = 0, this means that x_{n} satisfies the recurrence x_{n+1} = αx_{n}, whose solution is x_{n} = α^{n}x_{0}. For d > 0 we write

0 = (Δ(α)^{(d+1)})(x_{n}) = (Δ(α)^{(d)})[(Δ(α))(x_{n})].

By the inductive hypothesis, there is a polynomial p ∈ ℂ[n] of degree ≤ d−1 such that

(Δ(α))(x_{n}) = p(n)α^{n}.

The proof is therefore complete if we show that there exists a polynomial q ∈ ℂ[n] of degree deg(p) + 1 such that

α(q(n + 1) − q(n)) = p(n)

identically. This is readily done by invoking Lemma B.1 on page 37 and observing that^{2}

n_{(d)} = \frac{1}{d + 1}[(n + 1)_{(d+1)} − n_{(d+1)}],

since this ensures that we can write every monomial in the polynomial p in the desired form. □

The last displayed formula is the discrete analogue of the identity x^{d} = (d+1)^{−1} \frac{d}{dx}x^{d+1}. Indeed, as we already mentioned above, there is an analogy between recurrence relations and differential equations: in our case the role of the derivative is played by the forward difference operator Δ: X → X, which, by definition, has the following property:

Δ((x_{n})_{n∈ℕ}) =^{def} (x_{n+1} − x_{n})_{n∈ℕ}.

In other words, Δ = Δ(1). As in the case of differential equations, it can be proved (see Kelley & Peterson [14, §2.2]) that if the non-homogeneous part of the linear recurrence p belongs to the rather special class of functions A defined below, then a particular solution of (3.1) can be found in a suitable subclass of A. The analogy with differential equations, though, can not be pushed too far: in the fairly elementary cases treated in this paper most things have an exact counterpart in the field of differential equations, an important exception being the order reduction method of Section 3.6, but in general (for recurrences with variable coefficients, to mention one instance) results and techniques can both be very different.

2n_{(k)} is the falling factorial function defined by n_{(k)} =^{def} k!\binom{n}{k} = n(n − 1) ⋯ (n − k + 1).

Theorem 3.3 For all elements p of the set A defined by

A =^{def} ⋃_{C∈F} ker(C),

it is possible to find an elementary closed formula for

P (n) =^{def} ∑_{k=0}^{n} p(k).

Proof. By Lemma 3.2 every element p ∈ A has the form p(n) = r_{1}(n)α_{1}^{n} + ⋯ + r_{k}(n)α_{k}^{n}, where r_{i} ∈ ℂ[n] for i = 1, …, k, and the α_{i} are distinct complex numbers. It is sufficient to prove the theorem for an expression of the form p(n) = r(n)α^{n}. Now P (n) satisfies the first order recurrence x_{n} = x_{n−1} + p(n) and p ∈ ker(Δ(α)^{(deg(r)+1)}) by Lemma 3.2. Therefore P ∈ ker(Δ(α)^{(deg(r)+1)}) ∘ Δ(1), so that P ∈ A by Lemma 3.2 again. □

For ease of future reference, we state the following immediate Corollary.

Corollary 3.4 Let p ∈ A have the form p(n) = r(n)α^{n}, where r ∈ ℂ[n] has degree d, and α ∈ ℂ \ {0}, and let λ_{1}, …, λ_{r} be the distinct solutions of the characteristic equation (3.3) associated to the recurrence (3.1), with respective multiplicities μ_{1}, …, μ_{r}. Then any solution of the recurrence (3.1) belongs to

\begin{cases}ker(Δ(λ_{1})^{(μ_{1})} ∘ ⋯ ∘ Δ(λ_{r})^{(μ_{r})} ∘ Δ(α)^{(d+1)}), & if α ∉ \{λ_{1}, …, λ_{r}\}; \\ ker(Δ(λ_{1})^{(μ_{1}+d+1)} ∘ ⋯ ∘ Δ(λ_{r})^{(μ_{r})}), & if α = λ_{1}. \end{cases}

We give some examples in Appendix C.

We just saw that the elements of the set A satisfy some homogeneous recurrence relation of finite order with constant coefficients, and that polynomials, exponentials and their products all belong to A. But, since we allow complex values of α, we see that the functions sin(nθ), cos(nθ) and their products also belong to A, and in general, any linear combination of terms of the form n^{k}ρ^{n} cos(nθ) and n^{k}ρ^{n} sin(nθ), for k ∈ ℕ, ρ, θ ∈ ℝ. Indeed, n^{k}ρ^{n} cos(nθ) and n^{k}ρ^{n} sin(nθ) both belong to ker(Δ(ρe^{iθ})^{(k+1)}), and also to the kernel of the real operator

C = (E^{(2)} − 2(ρ cos θ)E + ρ^{2}I)^{(k+1)}.

See Appendix B for the description of algorithms for finding the particular solution of a non-homogeneous recurrence in the cases just stated, and Appendix A.8 for the computation of the coefficients α in the relation (3.5).

3.3 Special Case: Recurrences of Order 1

The next few subsections are devoted to the special but important case of recurrences of order 1, because we want to emphasize again some of the features of the solution of recurrences without having to deal with irrelevant details.

It is quite easy to prove by induction that the solution of the first-order recurrence

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

is given by the formula

x_{n} = α^{n}x_{0} + ∑_{k=1}^{n} α^{n−k}p(k).  (3.7)

We consider the special case x_{0} = 0, and set

\begin{cases}y_{n} = αy_{n−1} + p(n), & for n ≥ 1; \\ y_{0} = 0, \end{cases}

whence, by (3.7),

y_{n} = ∑_{k=1}^{n} α^{n−k}p(k).  (3.8)

Clearly, we just have to solve the last recurrence since then x_{n} = α^{n}x_{0} + y_{n}. We also remark that, by linearity, if p is the sum of several functions, p_{1}+p_{2}, say, by the superposition principle (or, more simply, by (3.8)) it is enough to solve two distinct recurrences, one with p_{1} in place of p, the other with p_{2}, and then add the results. See Appendix C.1 for some concrete examples.

3.3.1 Polynomials and Exponentials

We insist on an important theoretical point: we often need the closed formula for the sum of the product of a polynomial and an exponential. In some cases, the answer is the product of the same exponential function and another polynomial of the same degree, and in some cases the degree increases: see the abstract discussion in Section 3.2 and in particular Corollary 3.4; for concrete cases, see Examples C.2, C.3, C.4, C.6, C.9, C.10.

Concretely, we are looking for a closed formula for the sum in (3.8) for some α ∈ ℂ, where p is q(n)α^{n} and q is a polynomial. In this case we see that the closed formula has again the shape of the product of a polynomial and α^{n}, but the degree of the polynomial has increased by 1. We see why it has to be so by considering, more generally, the recurrence

x_{n} = αx_{n−1} + q(n)β^{n},

where α and β are fixed, non-zero complex numbers, and q is a polynomial. Using (3.8), we look for a closed formula for the sum

∑_{k=1}^{n} α^{n−k}β^{k}q(k) = α^{n} ∑_{k=1}^{n}(\frac{β}{α})^{k}q(k) = α^{n} ∑_{k=1}^{n} γ^{k}q(k),  (3.9)

say. We claim that the case γ = 1 (that is, α = β) is quite special. In fact, if r is a polynomial of degree d, say, then

γ^{k}r(k) − γ^{k−1}r(k − 1) = γ^{k}(r(k) − γ^{−1}r(k − 1)).

The degree of the polynomial on the right is d if and only if γ ≠ 1, and is d − 1 otherwise. To see why this fact is relevant, assume that we can find a polynomial r such that q(k) = r(k) − γ^{−1}r(k − 1) identically. Now we can compute a closed formula for the sum in (3.9), far right:

\begin{aligned}∑_{k=1}^{n} γ^{k}q(k) & = ∑_{k=1}^{n} γ^{k}(r(k) − γ^{−1}r(k − 1)) = ∑_{k=1}^{n} γ^{k}r(k) − ∑_{k=1}^{n} γ^{k−1}r(k − 1) \\ & = ∑_{k=1}^{n} γ^{k}r(k) − ∑_{k=0}^{n−1} γ^{k}r(k) = γ^{n}r(n) − r(0). \end{aligned}

This implies that the closed formula for the right hand side of (3.9) contains a polynomial of the same degree as q if α ≠ β, and of higher degree if α = β. The facts just stated are analogous to the computation of derivatives and primitives. We say more on this topic in the Examples quoted at the beginning of this section.

Lemma 3.5 Let q ∈ ℂ[x] be any polynomial, and γ ∈ ℂ \ {0}. There exists a polynomial r ∈ ℂ[x] such that

q(x) = r(x) − γ^{−1}r(x − 1)

for all x ∈ ℂ. Furthermore, if γ ≠ 1 then deg(r) = deg(q), while if γ = 1 then deg(r) = deg(q) + 1.

Proof. If q is constant, q(x) = a, say, then we may take r(x) = aγ/(γ − 1) if γ ≠ 1, and r(x) = ax if γ = 1. If d = deg(q) ≥ 1 we proceed by induction: let a_{d} = lead(q), and set s(x) = a_{d}γx^{d}/(γ −1) if γ ≠ 1, and s(x) = a_{d}x^{d+1}/d if γ = 1. Let t(x) = q(x) − (s(x) − γ^{−1}s(x − 1)). A simple check reveals that deg(t) < deg(q), and therefore, by the induction hypothesis, there exists a polynomial r_{1} ∈ ℂ[x] such that t(x) = r_{1}(x) − γ^{−1}r_{1}(x − 1) for all x ∈ ℂ. Hence r(x) = s(x) + r_{1}(x) is the desired polynomial. □

Notice that, if q ∈ ℝ[x], then r ∈ ℝ[x] and, similarly, q ∈ ℚ[x] implies r ∈ ℚ[x].

3.4 Special Case: Recurrences of Order 2

The case of recurrences of order 2 is also interesting because it is now possible to have characteristic equations with multiple roots, and again we can put emphasis on quite important facts in a fairly simple situation.

We now want to solve general recurrences of order 2, like

x_{n} = αx_{n−1} + βx_{n−2} + p(n)  (3.10)

where p is a function defined over the natural numbers. Here we assume that β ≠ 0, for otherwise the recurrence would have order 1. We set

\begin{cases}g_{n} = αg_{n−1} + βg_{n−2}, & for n ≥ 2; \\ g_{0} = 1, \\ g_{1} = α, \end{cases}  (3.11)

for the fundamental solution. This is readily solved by means of the technique described in Section 3. It is not difficult to prove by induction that, for n ≥ 2, the solution of (3.10) is

x_{n} = g_{n−1}x_{1} + βg_{n−2}x_{0} + ∑_{k=2}^{n} g_{n−k}p(k).  (3.12)

Here and below, the solutions, as they stand, are only valid for n ≥ 2. We observe that if λ_{1} ≠ λ_{2}, then g_{n} has the explicit expression

g_{n} = \frac{λ_{1}^{n+1} − λ_{2}^{n+1}}{λ_{1} − λ_{2}} .

If λ_{1} = λ_{2}, that is α^{2} + 4β = 0, the formula is different, though the idea is similar; see Example C.6 on page 45. In fact (as in the case of the Fibonacci numbers, see Example C.5 on page 44) we have

g_{n} = aλ_{1}^{n} + bλ_{2}^{n},

for suitable complex numbers a and b, and for n = 0 and n = 1 we find

\begin{cases}a + b = 1, \\ aλ_{1} + bλ_{2} = α, \end{cases}  whence  \begin{cases}a = λ_{1}(λ_{1} − λ_{2})^{−1}, \\ b = −λ_{2}(λ_{1} − λ_{2})^{−1}. \end{cases}

Thus, by (3.12), for n ≥ 2 we have

x_{n} = g_{n−1}x_{1} + βg_{n−2}x_{0} + ∑_{k=2}^{n} \frac{λ_{1}^{n+1−k} − λ_{2}^{n+1−k}}{λ_{1} − λ_{2}} p(k)  (3.13)

= g_{n−1}x_{1} + βg_{n−2}x_{0} + \frac{λ_{1}^{n+1}}{λ_{1} − λ_{2}} ∑_{k=2}^{n} λ_{1}^{−k} p(k) − \frac{λ_{2}^{n+1}}{λ_{1} − λ_{2}} ∑_{k=2}^{n} λ_{2}^{−k} p(k).

If p is in the set A defined in Section 3.2, the sums above also belong to A, and there is an algorithm to compute the closed formula (see Appendix B). A similar result also holds for linear recurrences of any finite order with constant coefficients: see Section 3.5 for the details.

3.4.1 “Guessing” the Solution in Special Cases

We now show how to solve directly equation (3.10) when p is a polynomial of degree d, say: we first compute the homogeneous part x_{1}g_{n−1} + βx_{0}g_{n−2} and seek a particular solution of the complete equation (a different approach is based on the techniques of Appendix B). Corollary 3.4 implies that if p is a polynomial and λ = 1 is not a solution of the characteristic equation, then the particular solution itself is a polynomial of the same degree as p, whose coefficients we have to find. This leads to the solution of a set of linear equations with d + 1 unknowns. If λ = 1 is a solution of the characteristic equation of multiplicity μ, we look for a polynomial of degree d + μ.

Example 3.6 Solve

\begin{cases}x_{n} = αx_{n−1} + βx_{n−2} + p(n), & for n ≥ 2; \\ x_{0} = a, \\ x_{1} = b. \end{cases}  (3.14)

The homogeneous fundamental equation is

\begin{cases}g_{n} = αg_{n−1} + βg_{n−2}, & for n ≥ 2; \\ g_{0} = 1, \\ g_{1} = α. \end{cases}

For the time being, we assume that λ = 1 is not a solution of the characteristic equation, which amounts to assuming that α + β ≠ 1. By Corollary 3.4 we know that there is a polynomial q such that

x_{n} = (b + q(1))g_{n−1} + β(a + q(0))g_{n−2} − q(n),

where q has to be computed as a function of p, which is the non-homogeneous part of the recurrence. If λ = 1 is not a solution of the characteristic equation, q can be determined by imposing that it has the same degree as p and satisfies the functional equation

p(n) = αq(n − 1) + βq(n − 2) − q(n).  (3.15)

This can be verified by induction. In order to find q it is necessary to solve a set of linear equations with d + 1 unknowns, d being the degree of p.

Example 3.7 Solve (3.14) for α = β = 1 and p(n) = n^{2}.

We have to find real numbers a, b and c such that

n^{2} = (a(n − 1)^{2} + b(n − 1) + c) + (a(n − 2)^{2} + b(n − 2) + c) − (an^{2} + bn + c) holds identically, that is, collecting and simplifying,

(a − 1)n^{2} + (b − 6a)n + (5a − 3b + c) = 0.

By the identity principle for polynomials, this implies that a = 1, b = 6, c = 13. Note (in agreement with Lemma 3.5) that if λ = 1 were a solution of the characteristic equation of multiplicity μ, the right hand side of (3.15) would be a polynomial of degree deg(q) − μ: that is why, if μ > 0, we have to look for a polynomial q of degree deg(p) + μ. We give more examples in Appendix C.2: see, in particular, Examples C.11–C.14.

3.5 Higher-Order Linear Recurrences with Constant Coefficients

The theory for the higher-order cases is analogous, but the computations are obviously heavier. Notice though that it is fairly easy to recognize special recurrences of order k > 2 such as x_{n} = a(n)x_{n−k} + p(n), which can be treated as a (much easier and faster to solve) first order recurrence. The point is that in this case x_{n} does not depend on x_{n−1}, x_{n−2}, …, x_{n−k+1}, so that it is possible to use a simpler method. A similar remark holds for x_{n} = x_{n−2} + x_{n−4}, which should be treated as an order 2 recurrence. For more on this topic, see Section 3.6.

We now give the general solution of (3.1); we introduce the fundamental solution of the associated homogeneous equation, which is

\begin{cases}g_{n} = a_{1}g_{n−1} + a_{2}g_{n−2} + ⋯ + a_{k}g_{n−k}, & for n ≥ k; \\ g_{0} = 1, \\ g_{n} = a_{1}g_{n−1} + a_{2}g_{n−2} + ⋯ + a_{n−1}g_{1} + a_{n}g_{0}, & for 1 ≤ n < k. \end{cases}

The general solution of (3.1) is then given by the formula

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}),  (3.16)

which is valid for n ≥ k and can be proved by induction. This is essentially equation (5) of Cohen & Katcoff [10], and it contains as special cases both (3.7) and (3.12). We remark that the two sums in (3.16) correspond naturally to the non-homogeneous part p and to the initial conditions, respectively. Furthermore, the quantities in brackets on the right depend only on the initial values and on the coefficients of the recurrence.

3.6 Order Reduction

Suppose that we want to solve the recurrence

x_{n} = αx_{n−k} + p(n),

where k ≥ 2 is a fixed integer. As we said above, this is essentially a recurrence of order 1. We remark that the variable n runs over integers that are congruent mod k. Hence, setting r =^{def} n mod k, n =^{def} km + r and y_{m}^{(r)} =^{def} x_{km+r}, we consider the k recurrences

\begin{cases}y_{m}^{(r)} =^{def} αy_{m−1}^{(r)} + p(km + r), & for m ≥ 1; \\ y_{0}^{(r)} =^{def} x_{r}, \end{cases}  (3.17)

which are actually of order 1 and can be solved by means of the techniques described in Section 3.3. It is not really necessary to solve k recurrences (one for each possible value of r) since, in a sense, the recurrences (3.17), which may actually be distinct as Example 3.8 below shows, can be considered as a single, parametric recurrence. Here we use the term “parametric” in a quite narrow sense, since the parameter r is only allowed to lie in the finite set {0, 1, …, k − 1}. Elsewhere in this document, a parameter is usually allowed to take any real (or indeed complex) value, without restrictions.

The advantage over the general method is best seen in the light of a few examples. Similar remarks, of course, apply for recurrences of the form x_{n} = αx_{n−k} + βx_{n−2k} + p(n).

Example 3.8 Solve the recurrence

x_{n} = x_{n−2} + n.

The general method gives the solution

x_{n} = \frac{1}{2}(1 + (−1)^{n})x_{0} + \frac{1}{2}(1 − (−1)^{n})x_{1} + \frac{1}{4}n^{2} + \frac{1}{2}n + \frac{3}{8}((−1)^{n} − 1). (3.18)

On the other hand, the substitution above leads to the recurrence

\begin{cases}y_{m}^{(r)} =^{def} y_{m−1}^{(r)} + 2m + r, & for m ≥ 1; \\ y_{0}^{(r)} =^{def} x_{r}, \end{cases}

which is parametric in the parameter r in the above sense. Its solution is

y_{m}^{(r)} = m^{2} + (1 + r)m + y_{0}^{(r)},

so that, substituting again m = \frac{1}{2}(n − r), we obtain

x_{n} = x_{r} + \frac{1}{4}n^{2} + \frac{1}{2}n − \frac{1}{4}r(r + 2),  (3.19)

where r = n mod 2 ∈ {0, 1}. It is immediate to check that this is the same as (3.18). It is also clear that (3.19) is more readable than (3.18), where there is a complex part depending on the initial conditions.

The following example is perhaps even more striking.

Example 3.9 Solve

x_{n} = x_{n−3}.

The reduction method yields at once x_{n} = x_{n mod 3}. The standard method entails the solution of the characteristic equation λ^{3} −1 = 0, whose roots are λ_{1} = 1, λ_{2} = e^{2πi/3} and λ_{3} = e^{−2πi/3}. For simplicity, we write ω = e^{2πi/3}, so that the roots are 1, ω, \overline{ω} = ω^{2}. We also need to solve a system of three linear equations in order to find suitable constants a_{1}, a_{2} and a_{3} so that

a_{1} + a_{2}ω^{n} + a_{3}ω^{2n} = x_{n} for n = 0, 1, 2.

The solution is therefore

x_{n} = \frac{1}{3}(x_{0} + x_{1} + x_{2}) + \frac{1}{3}(x_{0} + ω^{2}x_{1} + ωx_{2})ω^{n} + \frac{1}{3}(x_{0} + ωx_{1} + ω^{2}x_{2})ω^{2n}

= (1 + ω^{n} + ω^{2n})\frac{x_{0}}{3} + (1 + ω^{n+2} + ω^{2n+1})\frac{x_{1}}{3} + (1 + ω^{n+1} + ω^{2n+2})\frac{x_{2}}{3} .

Though both forms of the solution are equally valid, there is no doubt that the latter is much more complicated than the former, and that it is more difficult to find. It is also plain that the reduction of the order allows one to solve x_{n} = x_{n−k} for fixed integral, positive k, with the same amount of computation, regardless of k, whereas the standard method, even in this almost trivial case, requires computations whose cost increases with k. Other examples that stress the convenience of this method are in Appendix C.4.

The remainder of this section is devoted to the transformation of the solution obtained by means of the order reduction into the general solution that would be obtained using the standard method: in other words, we show in general how to deduce the exact form of the solution (3.18) from (3.19).

There are essentially two problems: expressing the terms x_{r} where r = n mod k, the initial conditions that appear in the “reduced” solution, and expressing r itself. The examples above suggest that the first problem is solved introducing suitable linear combinations of k-th roots of unity (the complex solutions of the polynomial equation z^{k} = 1), as the following Lemma proves. The same result solves also the second problem.

Lemma 3.10 For any positive integer k and any r ∈ {0, …, k−1} there is a linear combination of the k-th roots of unity ω_{1}, …, ω_{k}, with coefficients α_{1}(r), …, α_{k}(r) such that for any n ∈ ℤ we have

α_{1}(r)ω_{1}^{n} + ⋯ + α_{k}(r)ω_{k}^{n} = \begin{cases}1 & if n ≡ r mod k, \\ 0 & otherwise. \end{cases}

Proof. It is sufficient to take α_{j}(r) = ω_{j}^{−r}/k, since ω_{1}^{n} + ⋯ + ω_{k}^{n} = k if k divides n and is 0 otherwise. □

It is now possible to express the part depending on the initial conditions. Let

η_{r,k}(n) =^{def} \frac{1}{k} ∑_{h=1}^{k} ω_{h}^{n−r},

so that η_{r,k}(n) = 1 if n ≡ r mod k, and η_{r,k}(n) = 0 otherwise, by the observation in the proof of Lemma 3.10. Therefore we have

x_{r} = ∑_{s=0}^{k−1} η_{s,k}(n)x_{s}.  (3.20)

It might seem at first sight that the right hand side does not depend on r, but on second thoughts it will be seen that η_{r,k} is periodic with period k, since the ω’s are k-th roots of unity, so that the expression on the right depends only on n mod k = r, and the remark after the definition of the functions η_{r,k} shows that (3.20) is actually correct.

Another important consequence of Lemma 3.10 is the fact that

n mod k = ∑_{s=0}^{k−1} sη_{s,k}(n) = \frac{1}{k} ∑_{s=0}^{k−1} s ∑_{h=1}^{k} ω_{h}^{n−s} = \frac{1}{k} ∑_{h=1}^{k} ω_{h}^{n} ∑_{s=0}^{k−1} sω_{h}^{−s}. (3.21)

This expression is really too complicated to be of any practical use, but a neat trick using some identities of Appendix B.2 will help us in finding an equivalent expression which is much easier to handle, providing a feasible solution to our second problem.

Lemma 3.11 Let k be a positive integer, and let 1 = ω_{1}, ω_{2}, …, ω_{k} be the complex k-th roots of unity. For any integer n ∈ ℤ we have

n mod k = \frac{1}{2}(k − 1) + ∑_{h=2}^{k} \frac{ω_{h}}{1 − ω_{h}} ω_{h}^{n}.

Proof. Let z = ω_{h}^{−1} for h ≠ 1 and recall that ω_{h}^{k} = 1 by definition. Since z ≠ 1, by (B.5) with j = 1 we have

∑_{s=0}^{k−1} sz^{s} = \frac{z}{(1 − z)^{2}} (1 − kz^{k−1} + (k − 1)z^{k}) = \frac{kω_{h}}{1 − ω_{h}} .

Therefore, by (3.21)

n mod k = \frac{1}{k} ∑_{h=1}^{k} ω_{h}^{n} ∑_{s=0}^{k−1} sω_{h}^{−s} = \frac{1}{2}(k − 1) + ∑_{h=2}^{k} \frac{ω_{h}}{1 − ω_{h}} ω_{h}^{n},

as stated. □

3.7 Summary of the Method

Let us summarize the overall solution methodology for linear recurrences of finite order with constant coefficients: The very first step is to recognize whether we can apply the order reduction method, which is always beneficial as it allows to work with polynomials of lower degree. After order reduction, the homogeneous part of the recurrence is solved by the characteristic equation method, which involves finding the roots of a polynomial along with their multiplicities. Now, to obtain the general solution of the complete equation, it is sufficient to find any one of its solutions (which is then called the particular solution), add it to the general solution of the homogeneous equation found with the previous step, and impose that the initial conditions hold. For the recurrences of the first or the second order when the characteristic equation has simple roots, a non-closed formula for the solution is simply given by (3.7) and (3.13), respectively. For higherorder recurrences (including order 2 when the characteristic equation has a multiple root), the general solution is given by (3.16), which requires the symbolic solution of a system of equations. Whatever the order is, finding a closed-form solution requires the computation of symbolic sums that, depending on the inhomogeneous term of the recurrence, can be arbitrarily difficult when not plainly impossible (see Appendix B for the symbolic summation of an important class of functions). In some special cases, i.e., when the non homogeneous term p belongs to the class A defined in Section 3.2, we have shown an alternative way for “guessing” the particular solution. Indeed, when p ∈ A, Lemma 3.2 and Corollary 3.4 characterize the shape of the solution: in other words we know that it is the (sum of several terms of the form) product of a polynomial and an exponential function, where we know the degree of the polynomial and the base of the exponential function. The actual values of the coefficients in the polynomial can be determined by imposing that the solution really satisfies the initial recurrence. A similar strategy is used in Examples 3.7 and C.11–C.14.

In conclusion, the above procedure allows to reduce the problem of solving LRFOCCs to (1) the solution of polynomial equations and (2) the computation of closed forms for symbolic sums. In the appendix we illustrate some useful techniques for dealing with (1) and (2); other techniques exist and may be worth exploring, since any advancement in tackling problems (1) and (2) will extend the class of recurrences that can be solved by this method.

4 Systems of Linear Recurrences

In this section, we show how to solve systems of linear recurrences of finite order with constant coefficients. We first deal with systems of recurrences where each recurrence is of order 1. These have been studied, e.g., by Cohen & Katcoff [10] and by Kelley & Peterson [14]. We then turn our attention to the solution of systems of recurrences of any order and show that they can all be reconducted to systems of order-1 recurrences.

4.1 Systems of Recurrences of Order 1

Sets of k recurrence relations with constant coefficients of order 1 may be tackled in two different, but essentially equivalent, ways, that are quite similar to the solution of sets of linear equations: one can either solve for one of the “variables” (a recurrence in our case) and then plug the result in the remaining equations, iterating if necessary in order to get only one equation, or see the problem from a more abstract point of view, the system being a linear operator from the set of k-tuples of recurrences to itself, so that the general theory developed in Section 3.2 applies. The first approach yields k recurrences that contain only one “variable,” each recurrence being of order k, all with the same characteristic equation of degree k; the second one needs the solution of the characteristic equation of the linear operator referred to above, which has degree k and is exactly the same as the characteristic equation of the linear operator: that is why the two approaches are essentially equivalent. We give examples of both approaches below. We begin with a set of two first-order recurrences with constant coefficients in two variables, and we show how to transform it into two second-order recurrences, each with only one variable. Suppose we want to solve

\begin{cases}x_{n} = ax_{n−1} + by_{n−1} + p(n), \\ y_{n} = cx_{n−1} + dy_{n−1} + q(n). \end{cases}  (4.1)

The second equation with n replaced by n−1 yields cx_{n−2} = y_{n−1} −dy_{n−2} − q(n − 1). We now replace n by n − 1 in the first equation and multiply throughout by c. Plugging these two relations into the second equation above, after some simplifications we obtain

y_{n} = (a + d)y_{n−1} − (ad − bc)y_{n−2} + q(n) − aq(n − 1) + cp(n − 1).

A similar argument gives

x_{n} = (a + d)x_{n−1} − (ad − bc)x_{n−2} + p(n) − dp(n − 1) + bq(n − 1).

Observe that the equations (4.1) imply that x_{1} = ax_{0} + by_{0} + p(1) and that y_{1} = cx_{0} + dy_{0} + q(1). Notice also that the coefficients of the new recurrences are the trace and the opposite of the determinant of the matrix of the coefficients in the system (4.1). This discussion means that we can always transform a system of two linear recurrences of the first order with constant coefficients into a pair of linear recurrences of the second order, with the variables decoupled.

The procedure is quite general, and, as a more difficult example, we show how to reduce a system of 5 equations of order 1 to a single equation of order 5 (this is Example 12 in Table II of the paper by Cohen & Katcoff [10]). Actually, one should repeat the same procedure in order to get 5 equations, one involving a single variable, but it is also possible to find the solution to one and then substitute in the original problem, repeating until all equations have been solved.

Example 4.1 Solve

\begin{cases}a_{n} = a_{n−1} + b_{n−1} + 2^{n}, \\ b_{n} = 5a_{n−1} + c_{n−1}, \\ c_{n} = −5a_{n−1} + d_{n−1} − 3^{n}, \\ d_{n} = −4a_{n−1} + e_{n−1}, \\ e_{n} = 4a_{n−1}. \end{cases}  (4.2)

The last equation implies that e_{n−1} = 4a_{n−2}. We substitute this into (4.2) and, omitting from now on the last equation, we find

\begin{cases}a_{n} = a_{n−1} + b_{n−1} + 2^{n}, \\ b_{n} = 5a_{n−1} + c_{n−1}, \\ c_{n} = −5a_{n−1} + d_{n−1} − 3^{n}, \\ d_{n} = −4a_{n−1} + 4a_{n−2}. \end{cases}

The last equation shows that d_{n−1} = −4a_{n−2} + 4a_{n−3}, and we substitute again:

\begin{cases}a_{n} = a_{n−1} + b_{n−1} + 2^{n}, \\ b_{n} = 5a_{n−1} + c_{n−1}, \\ c_{n} = −5a_{n−1} − 4a_{n−2} + 4a_{n−3} − 3^{n}. \end{cases}

We obtain c_{n−1} = −5a_{n−2} − 4a_{n−3} + 4a_{n−4} − 3^{n−1}, whence

\begin{cases}a_{n} = a_{n−1} + b_{n−1} + 2^{n}, \\ b_{n} = 5a_{n−1} − 5a_{n−2} − 4a_{n−3} + 4a_{n−4} − 3^{n−1}. \end{cases}

Finally b_{n−1} = 5a_{n−2} − 5a_{n−3} − 4a_{n−4} + 4a_{n−5} − 3^{n−2}, so that we are left with

a_{n} = a_{n−1} + 5a_{n−2} − 5a_{n−3} − 4a_{n−4} + 4a_{n−5} − 3^{n−2} + 2^{n}.

The characteristic equation factors as (λ − 1)^{2}(λ + 1)(λ − 2)(λ + 2) = 0, so that this recurrence is easily solved with the technique described in Section 3. Once a closed formula for a_{n} is found, it can be used to find the recurrences satisfied by the other variables, until all the original equations are solved.

This mechanical procedure shows quite clearly that we can eliminate one variable at the cost of increasing the order of the other equations by 1.

Kelley & Peterson [14, §4.1] describe a different method. Given the set of k linear recurrences with constant coefficients

\begin{cases}x_{n}^{(1)} = a_{11}x_{n−1}^{(1)} + ⋯ + a_{1k}x_{n−1}^{(k)} + p_{1}(n), \\ … \\ x_{n}^{(k)} = a_{k1}x_{n−1}^{(1)} + ⋯ + a_{kk}x_{n−1}^{(k)} + p_{k}(n), \end{cases}

we introduce the vector equation

x_{n} = Ax_{n−1} + p(n),  (4.3)

where

x_{n} =^{def} [ \begin{matrix}x_{n}^{(1)} \\ … \\ x_{n}^{(k)} \end{matrix} ] , A =^{def} [ \begin{matrix}a_{11} & & & … & & & a_{1k} \\ . & . & . & . & . & . & . & . & . & . & . & . & . \\ a_{k1} & & & … & & & a_{kk} \end{matrix} ] , p(n) =^{def} [ \begin{matrix}p_{1}(n) \\ … \\ p_{k}(n) \end{matrix} ] . (4.4)

As above, we first study the homogeneous equation x_{n} = Ax_{n−1}, whose solution is, obviously,

x_{n} = A^{n}x_{0}.  (4.5)

We recall that, by definition, an eigenvector of A relative to the eigenvalue λ is a non-zero vector v with the property that Av = λv. Hence it is immediate that x_{n} = λ^{n}v satisfies (4.5) with the initial condition x_{0} = v. If A has k linearly independent eigenvectors (which certainly is the case if it has k distinct eigenvalues) then there is an invertible matrix M such that A = M^{−1}BM, where B is diagonal, with the eigenvalues of A on the main diagonal. Let’s denote by λ_{1}, …, λ_{k} these eigenvalues, so that B_{ij} = λ_{i} if i = j and 0 otherwise. We remark that B^{n} is easily computed since it is diagonal in its own turn, with the n-th powers of the λ’s on the main diagonal. Since A^{n} = M^{−1}B^{n}M, we have that the solution of the homogeneous equation is

x_{n} = M^{−1}B^{n}Mx_{0}.

Actually, there is no need to find the matrix M explicitly: indeed, if A has k linearly independent eigenvectors, they necessarily form a basis of ℝ^{k} (with some care if there are complex eigenvalues): since by definition the eigenvector v_{j} relative to the eigenvalue λ_{j} has the property Av_{j} = λ_{j}v_{j}, and every vector x ∈ ℝ^{k} is a linear combination of eigenvectors, we can write the solution of (4.5) in the form

x_{n} = μ_{1}λ_{1}^{n}v_{1} + ⋯ + μ_{k}λ_{k}^{n}v_{k},

where x_{0} = μ_{1}v_{1} + ⋯ + μ_{k}v_{k}. Anyway, computing this decomposition is equivalent to the inversion of the matrix M above. We remark that the situation is quite similar to the case of a single recurrence of order k: indeed, it can be proved that they are essentially equivalent.

Now we turn our attention to the complete equation (4.3). Formally, the solution of (4.3) is given by the formula (3.7):

x_{n} = A^{n}x_{0} + ∑_{k=1}^{n} A^{n−k}p(k).

The point here is that, since A satisfies its own characteristic equation det(A − λI_{k}) = 0 (where I_{k} is the identity matrix of order k), A^{n} can be written as a linear combination of A^{0} = I_{k}, A, …, A^{k−1}. For instance, in the case of the system (4.2), we have

A = [\begin{matrix}1 & 1 & 0 & 0 & 0 \\ 5 & 0 & 1 & 0 & 0 \\ −5 & 0 & 0 & 1 & 0 \\ −4 & 0 & 0 & 0 & 1 \\ 4 & 0 & 0 & 0 & 0 \end{matrix}]

whose characteristic equation is the same as above, namely

det(A − λI_{k}) = −(λ − 1)^{2}(λ + 1)(λ − 2)(λ + 2) = 0.

This means that we can avoid the explicit substitution, as we did above.

4.2 Systems of Recurrences of Higher Order

In the previous section we only considered systems of equations of order 1, although the first method developed above entails the transformation of the equations into higher-order ones. From a purely theoretical point of view, it is quite interesting to know that it is possible to transform any system of linear recurrences with constant coefficients of any finite order into an equivalent system with a larger number of variables, with the property that every equation has order exactly 1. From a computational point of view, such a transformation is probably inefficient, so that the first method explained above should usually be the preferred one.

The idea is quite simple: it is possible to reduce by 1 a suitable measure of the complexity of an equation of order at least 2 (not necessarily the order; we give precise definitions below), at the cost of introducing exactly one more variable and one more equation of order 1: this is, essentially, the inverse process of the first method above where we reduced the number of variables by increasing the order of some equations. The total complexity of the system, which is the sum of the complexities of the equations, does not change, while the complexity of an equation of order at least 2 decreases. Iterating the process, it is possible to eliminate all equations of order 2 or more, until we are left with a (larger) system with only order 1 equations.

In order to avoid exceedingly cumbersome notations, we stick to the vector notation introduced in (4.4), with a slight generalization.

Notation 4.2 Let k and K be fixed, positive integers. For r = 1, …, K let

A^{(r)} =^{def} [ \begin{matrix}a_{11}^{(r)} & & … & & & a_{1k}^{(r)} \\ . & . & . & . & . & . & . & . & . & . & . & . & . & . \\ a_{k1}^{(r)} & & … & & & a_{kk}^{(r)} \end{matrix} ] .

We are interested in the system

x_{n} = A^{(1)}x_{n−1} + ⋯ + A^{(K)}x_{n−K} + p(n),  (4.6)

for fixed positive integers k (the number of variables) and K, where, in complete analogy with the discussion of linear recurrences of finite order with constant coefficients, we tacitly assume that A^{(K)} is not the zero matrix. Here it does not quite make sense to say that the order of the system is K, and we give a more suitable definition.

Definition 4.3 (Weight.) For each s = 1, …, K let the weight w_{s} of the s-th equation of the system (4.6) be defined by

w_{s} =^{def} ∑_{i=1}^{K} ∑_{\begin{matrix}j=1 \\ a_{sj}^{(i)}≠0 \end{matrix}}^{k} i.

We let the weight of the system (4.6) be w_{1} + ⋯ + w_{K}.

We are weighing each equation not only with the number of variables that actually occur in it (so that a_{sj}^{(i)} ≠ 0), but also taking into account the order of each occurring variable. The weight, as defined here, is not the same thing as the order, so the above is not necessarily coherent with the case k = 1 (a single equation) nor with the case of k equations of order 1. We only need the concept of weight to prove that the strategy outlined above actually works.

We now formally describe a strategy for transforming a system like (4.6) with at least an equation of order 2 or more (if all equations have order 1 there is nothing to do) into an equivalent one with the same weight, one more variable and one more equation, and such that the weight of one equation has decreased by 1. We may assume that the first equation, say, has order q > 1 and weight w_{1}. The strategy is as follows:

• let i be an index such that a_{1i}^{(q)} ≠ 0;

• add the new equation x_{n}^{(k+1)} = x_{n−1}^{(i)} to the system;

• replace the instance of x_{n−q}^{(i)} in the first equation by x_{n−q+1}^{(k+1)}, which is equivalent by definition.

We have plainly added one equation of order 1 and one variable, and also changed the first equation of the system: we have to prove that

• the weight of the new system is the same as the weight of the old one;

• the weight of the first equation has decreased by 1.

The first thing can be checked easily: in fact, the new equation in the first row has weight w_{1} − 1 (we left every term untouched except for the i-th in the q-th matrix, which has been replaced by one of order one unit less), and the equation we have added has weight 1. Therefore the new system has the same weight as the original one. The second item above has been verified incidentally.

Example 4.4 Using the strategy above, rewrite the system

\begin{cases}x_{n} = x_{n−2}, \\ y_{n} = x_{n−1} + y_{n−1}. \end{cases}

Both equations have weight 2, and the system has weight 4. We add a new variable z_{n} satisfying z_{n} = x_{n−1}. The system becomes

\begin{cases}x_{n} = z_{n−1}, \\ y_{n} = x_{n−1} + y_{n−1}, \\ z_{n} = x_{n−1}. \end{cases}

The equations have weight 1, 2, 1 respectively, and each equation now has order 1, while the system has still weight 4. In terms of matrices, for the initial system we have

A^{(1)} = [ \begin{matrix}0 & 0 \\ 1 & 1 \end{matrix} ] , A^{(2)} = [ \begin{matrix}1 & 0 \\ 0 & 0 \end{matrix} ] ,

and for the transformed system we have

B^{(1)} = [ \begin{matrix}0 & 0 & 1 \\ 1 & 1 & 0 \\ 1 & 0 & 0 \end{matrix} ]  so that  [ \begin{matrix}x_{n} \\ y_{n} \\ z_{n} \end{matrix} ] = B^{(1)} [ \begin{matrix}x_{n−1} \\ y_{n−1} \\ z_{n−1} \end{matrix} ] .

Adding the new equation z_{n} = x_{n−1} corresponds to adding a new row to the matrix A^{(1)} (with a zero entry in all columns, except for the first), and a new column (again, with a zero entry in all rows, except for the first).

Example 4.5 Using the strategy above, rewrite the system

\begin{cases}x_{n} = x_{n−2} + y_{n−3}, \\ y_{n} = x_{n−1} + y_{n−2}, \\ z_{n} = y_{n−1}. \end{cases}

Here the weights of the single equations are 5, 3 and 1, for a total weight of 9. Setting t_{n} = y_{n−1}, we obtain the new system

\begin{cases}x_{n} = x_{n−2} + t_{n−2}, \\ y_{n} = x_{n−1} + y_{n−2}, \\ z_{n} = y_{n−1}, \\ t_{n} = y_{n−1}, \end{cases}

where the equations have weights 4, 3, 1, 1 respectively, for a total weight of 9 again. It is plain that we can repeat the same procedure on the rewritten system until all equations have order 1.

5 Related Work

In this section we briefly outline the differences and similarities between our approach (and implementation) and the work by Lueker [18] and Cohen & Katcoff [10]. We use the notation of (3.1), page 3: in particular, we always denote the non-homogeneous part of the recurrence by p(n). Recall that we denote by C^{(k)} the k-th iterate of the operator C.

5.1 Lueker, 1980

Lueker [18] describes some useful techniques for solving a wide range of recurrences, presumably intended for manual computation.

As far as linear, first-order recurrences (both with constant and variable coefficients) are concerned, our approach is essentially the one of Lueker. The main differences arise when we consider non-homogeneous recurrences: we exploit the closed formula (3.7), (3.12) or (3.16) thereby transforming the problem of “guessing” a particular solution of the non-homogeneous recurrence into the problem of computing a finite sum, by means of the tools explained in Appendix B. Lueker, instead, introduces the shift operator E as we did in Section 3.2 and also (without saying it explicitly) the identity operator I.

The analysis is quite similar to our Section 3.2. Indeed, the point here is that for some forms of p(n) it is possible to find a suitable combination C of the operators E and I that annihilates p, that is C(p) ≡ 0, the constant sequence zero. In the customary language of linear algebra, p ∈ ker(C): see our Corollary 3.4. Then we apply the operator to both sides of the recurrence (3.1), and we are left with a new, homogeneous recurrence, which can be solved by means of the characteristic equation. We know that if p ∈ A, then there is such an operator C by Lemma 3.2.

The main drawback, however, is that the order of the new recurrence is higher than that of the original relation, since, for instance, the annihilator of the function p(n) = n·2^{n} +1 is the operator (E−2I)^{(2)}(E−I), so that the order of the recurrence with p as a non-homogeneous part would increase by 3. It is true that it is quite easy to detect the roots of the characteristic polynomial introduced in this way, but one is left in the end with the problem of determining the unknown coefficients in formula (3.4) or (3.5), that is, with a set of linear equations, whose number of variables is precisely the order of the new recurrence relation. Actually, the determination of the operator C that annihilates p is equivalent to the classification of the function p (polynomial, exponential, a product of a polynomial and an exponential, …). Thus, there is no particular advantage in this method, since we show in Appendix B how to find explicitly the particular solution in all of these cases.

Even though no experimental comparison of our technique and that proposed by Lueker was conducted (we are not aware of any available implementation of the latter), we believe our approach is more efficient, especially for the case of recurrences of low order that most frequently arises in automatic complexity analysis.

We now work out an example to illustrate the method described by Lueker. The same recurrence is tackled in Example C.2 on page 42 using the techniques of Section 3.

Example 5.1 Solve the recurrence

x_{n} = 2x_{n−1} + 2^{n−1} + n + 1.

We start by rewriting it as

x_{n+1} − 2x_{n} = 2^{n} + n + 2,

and then in the equivalent form

(E − 2I)(x_{n}) = 2^{n} + n + 2.

We already know that the operator E − 2I annihilates the sequence 2^{n}, so that

(E − 2I)(E − 2I)(x_{n}) = (E − 2I)(2^{n} + n + 2)
= (2^{n+1} + n + 3) − 2(2^{n} + n + 2)
= −n − 1.

We are left with a polynomial of degree 1: we know that the operator E − I reduces the degree of any polynomial by 1, and we need to apply it twice to get

(E − I)^{(2)}(E − 2I)^{(2)}(x_{n}) = (E − I)^{(2)}(−n − 1)
= (E − I)(−(n + 1) − 1 − (−n − 1))
= (E − I)(−1)
= 0.

We have therefore proved that the operator C =^{def} (E − I)^{(2)}(E − 2I)^{(2)} annihilates the sequence (x_{n}). In other words, we just found an order 4 homogeneous recurrence satisfied by (x_{n}), with the degree 4 characteristic equation

(λ − 1)^{2}(λ − 2)^{2} = 0.

Hence the general solution is

x_{n} = (αn + β) · 2^{n} + γn + δ,

where α, β, γ and δ are suitable real numbers that can be found by taking into account the fact that x_{1} = 2x_{0} + 3, x_{2} = 4x_{0} + 11, x_{3} = 8x_{0} + 30. These values can be computed directly from the recurrence we started with. Finally, we solve the corresponding system of equations (in the parameter x_{0}), that is, we solve

\{ \begin{matrix}& & β & & & + & δ & = & x_{0} \\ 2α & + & 2β & + & γ & + & δ & = & 2x_{0} + 3 \\ 8α & + & 4β & + & 2γ & + & δ & = & 4x_{0} + 11 \\ 24α & + & 8β & + & 3γ & + & δ & = & 8x_{0} + 30 \end{matrix}

and find the coefficients

α = \frac{1}{2},  β = x_{0} + 3,  γ = −1,  δ = −3.

5.2 Cohen & Katcoff, 1977

Cohen & Katcoff [10] describe an interactive computer program for solving difference equations; in this respect, our focus is different as we target completely algorithmic solution techniques. They present two methods for solving recurrences, namely, (1) “guessing” the solution of the non-homogeneous part,^{3} and (2) using the generating function. Since they “feel that neither of the preceding methods is suitable for direct use in the automatic solution of difference equations”, they use a modified form of the generating function method for solving (systems of) linear recurrences. In contrast, in this paper we describe several ways for “guessing” the particular solution of a non-homogeneous linear recurrence with constant coefficients by means of a simple formula, at least in the rather common cases where p belongs to the class A defined in Section 3.2: see Appendix B on page 36 for the details. In particular, there are at least 4 different ways to find the closed form for the sum of values of a polynomial, which are described in Section B.3, each probably with a different range of optimal performance. The main reason for our choice is that we want to solve or approximate a wider set of recurrences: when solving divide-et-impera recurrences [3], for example, we often need to compute symbolic sums, like the ones that appear in (3.16), that do not necessarily arise as solutions of linear recurrences. Therefore, we felt the need to tackle such problems directly.

3That is, using the result contained in our Corollary 3.4 that suggests the shape of the solution.

The overall solution algorithm we use in the PURRS system is quite similar to the one by Cohen & Katcoff, in that we first classify the recurrence (order, variable or constant coefficients, system), but then we use rather more direct ways for solving the simpler recurrences (that, as already noted, occur much more frequently in automatic complexity analysis).

Cohen & Katcoff give several examples of the actual computations carried out by their program [10, Table II]. PURRS easily solves their examples 1–7 (note that the last two contain symbolic coefficients) and should be soon able to solve also the systems 11–13. The more difficult recurrences 8–10 have variable coefficients, and belong to the class considered in [5]. Notice that PURRS gives exact answers whenever possible (including exact representation of rational numbers), whereas Cohen & Katcoff only give approximated coefficients rounded to the third decimal digit. We also point out that the solution they give for their example 3 is incorrect.^{4}

In their final remarks, Cohen & Katcoff discuss possible improvements of their work: apart from the problem of factoring polynomials of degree 5 or more which we discuss in Appendix A (item 1 in their list), we share their feelings about the simplification routines (item 4). We incorporated into PURRS many routines that simplify expressions, in particular those that contain irrational numbers that may arise when solving the characteristic equation associated to a linear recurrence. Finally, we devote a future paper of this series [4] to approximate solution, an interesting problem that they suggest in item 8.

4PURRS computes it correctly as a_{n} = \frac{4}{9} · 2^{n} + (−1)^{n}( \frac{4}{3} n − \frac{4}{9} ) + (−1)^{n}(1 − n)a_{0} + (−1)^{n+1}na_{1}. The correct solution with the given initial conditions a_{0} = 1, a_{1} = 2 is therefore a_{n} = \frac{4}{9} · 2^{n} + (−1)^{n}( \frac{5}{9} − \frac{5}{3} n).

6 Conclusion

We have described some algorithmic techniques for the efficient solution of a wide class of (systems of) linear recurrences of finite order with constant coefficients. The basic theory has been exposed both from an abstract and a more concrete point of view, so as to convey the general principles and to provide, at the same time, a rather precise specification for the implementation. In fact, variants of these techniques, along with the associated overall solution methodology and the tools illustrated in the appendix, have been implemented in PURRS, a system for the automatic solution and approximation of a much wider class of recurrences (see http://www.cs.unipr.it/purrs/).

In this paper, we only tackled the problem of “elementary” summation in closed form. In many interesting applications, though, the non-homogeneous part of the LRFOCC does not belong to the class A, and therefore the results described here can not be used. There is a wide class of non-homogeneous parts (the simplest element of this class being the function p(n) = n^{α}, where α is not necessarily integral) that can be summed approximately, with rather sharp upper and lower bounds for the solution, and we deal with this problem in general in [4]. Another class (that contains for instance some, but not all, rational functions) that admits closed-form solution will be introduced in the paper on “transcendental” summation methods [5] that deals with the recent holonomy theory of hypergeometric summation [20], and contains results on recurrences with variable coefficients.

References

[1] L. M. Adelman and A. M. Odlyzko. Irreducibility testing and factorization of polynomials. Mathematics of Computation, 41:699–709, 1983.

[2] T. M. Apostol. Introduction to Analytic Number Theory. Springer-Verlag, Berlin, 1975.

[3] R. Bagnara, A. Zaccagnini, and T. Zolo. The automatic solution of recurrence relations. II. “Divide et Impera” recurrences. Quaderno, Dipartimento di Matematica, Università di Parma, Italy. To appear.

[4] R. Bagnara, A. Zaccagnini, and T. Zolo. The automatic solution of recurrence relations. III. Approximate solutions. Quaderno, Dipartimento di Matematica, Università di Parma, Italy. To appear.

[5] R. Bagnara, A. Zaccagnini, and T. Zolo. The automatic solution of recurrence relations. IV. Linear recurrences of finite order with variable coefficients. Quaderno, Dipartimento di Matematica, Università di Parma, Italy. To appear.

[6] B. L. Burrows and R. F. Talbot. Sums of powers of integers. American Mathematical Monthly, 91:394–403, 1984.

[7] B. F. Caviness. On canonical forms and simplification. Journal of the ACM, 17(2):385–396, 1970.

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

[9] L. Childs. A Concrete Introduction to Higher Algebra. Springer-Verlag, Berlin, 1979.

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

[11] A. W. F. Edwards. A quick route to sums of powers. American Mathematical Monthly, 93:451–455, 1986.

[12] G. H. Hardy and E. M. Wright. An Introduction to the Theory of Numbers. Oxford Science Publications, 1979. Fifth Edition.

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

[14] W. G. Kelley and A. C. Peterson. Difference Equations. An Introduction with Applications. Academic Press, Inc. Harcourt Brace Jovanovich, 1991.

[15] A. Lenstra, H. Lenstra Jr., and L. Lovasz. Factoring polynomials with rational coefficients. Mathematische Annalen, 261:515–534, 1982.

[16] L. S. Levy. Summation of the series 1^{n} + 2^{n} + ⋯ + x^{n} using elementary calculus. American Mathematical Monthly, 77:840–847, 1970. See also [17].

[17] L. S. Levy. Corrigendum. American Mathematical Monthly, 78:987, 1971. See also [16].

[18] G. S. Lueker. Some techniques for solving recurrences. ACM Computing Surveys, 12(4):419–436, 1980.

[19] M. Mignotte and D. Ştefănescu. Polynomials. An Algorithmic Approach. Springer-Verlag, Berlin, 1999.

[20] M. Petkovšek, H. S. Wilf, and D. Zeilberger. A = B. A. K. Peters, 1996.

[21] D. Richardson. Some undecidable problems involving elementary functions of a real variable. Journal of Symbolic Logic, 33(4):514–520, 1968.

[22] P. S. Wang. The undecidability of the existence of zeros of real elementary functions. Journal of the ACM, 21(4):586–589, 1974.

[23] E. Weisstein. The world of mathematics. A web site full of informations

    on Mathematics. See http://mathworld.wolfram.com/.

A Solving Polynomial Equations

We recall that the general formula for the polynomial equation of degree 5 or more does not exist. A possibility is the following: first look for repeated (multiple) roots (see Appendix A.1), for the simple root x = 0 and for rational roots with “small” numerator and denominator (see Appendix A.2). Failing this, one can give control to an equation solver that detects polynomials of a special shape (see Appendix A.3). As a last resort, one can try to solve the general polynomial equation of degree up to 4. Polynomials of degree 1 or 2 are easy to treat, whereas polynomials of degree 3 (see Appendix A.4) and 4 (see Appendix A.5) need special care.

A.1 Detecting Multiple Roots

In order to detect multiple roots of a polynomial equation p(x) = 0 it is possible to use the square-free decomposition, which yields a partially factored form of p, in the case that p has repeated factors. This can be accomplished by computing gcd(p, p′), the greatest common divisor of the polynomial p and its derivative, since the rule for the derivative of a product immediately yields that if p(x) = g(x)^{k}h(x), say, where g and h are polynomials and k ≥ 2 is an integer, then p′(x) = kg′(x)g(x)^{k−1}h(x) + g(x)^{k}h′(x), so that g(x)^{k−1} is a non-trivial factor of gcd(p, p′). We recall that it is extremely important to detect repeated roots (if any) of the characteristic equation (see Section 3) and their multiplicity, and we do so at a very early stage of the computation. We now give precise definitions.

Definition A.1 (Square-free polynomial.) A polynomial p ∈ ℚ[x] is square-free if it is not divisible by the square of any polynomial q ∈ ℚ[x] of degree larger or equal to 1.

Definition A.2 (Square-free decomposition.) Let p ∈ ℤ[x] be a polynomial. We say that the decomposition

p(x) = a ∏_{i=1}^{r} p_{i}(x)^{α_{i}}

is a square-free decomposition of p if the following conditions hold:

• a ∈ ℤ \ {0};

• α_{i} ∈ ℕ \ {0} for i = 1, …, r and 1 ≤ α_{1} < α_{2} < ⋯ < α_{r};

• p_{i} ∈ ℤ[x] and deg(p_{i}) ≥ 1 for i = 1, …, r;

• gcd(p_{i}, p_{j}) = 1 if i ≠ j;

• p_{i} is square-free for i = 1, …, r.

We might express the last two conditions more concisely by saying that the polynomial

∏_{i=1}^{r} p_{i}(x)

is square-free.

A.2 Detecting Rational Roots

It is also possible to detect rational roots (if any) of any polynomial with integer coefficients. This depends on the following Theorem of Gauss (see Hardy & Wright [12, Theorem 45] for a variant).

Theorem A.3 If a_{0}, a_{1}, …, a_{n} are integers without common factors larger than 1, where a_{0} and a_{n} do not vanish, and x = p/q is a rational root of the polynomial equation

a_{n}x^{n} + ⋯ + a_{1}x + a_{0} = 0

(where p and q ≠ 0 are mutually prime integers) then p divides a_{0} and q divides a_{n}.

Looking for all rational roots would involve factoring the integers a_{0} and a_{n}, a notoriously difficult problem: in fact, at the moment, the best factoring algorithms are non-polynomial. So one can be content with detecting rational roots p/q, where |p| and |q| do not exceed some threshold M = m^{2}, with m ∈ ℕ small. The first step is to consider the divisors of a_{0} and a_{n}. Let x = a_{0} or x = a_{n}, |x| < M, be the integer whose divisors we are looking for. With a loop we check if the integer i ∈ {1,…, m − 1} divides x, and in this case we push onto a suitable structure both i and its “conjugate” divisor |x|/i, if they are different, of course. The process stops when i exceeds the square root of |x|. We remark that although most recurrence relations have a characteristic equation whose leading coefficient is 1, we have to take into account the fact that the coefficients may be rational numbers.

Much more ingenuity would be needed if it were desired to have a general routine for finding all rational roots. In this case, it would probably be better to use a general algorithm for splitting polynomials with integer coefficients such as the one described in Lenstra, Lenstra & Lovasz [15], rather than trying to factor completely the integers a_{0} and a_{n}. A simpler method is described in Childs [9, Part II, Chapters 11 & 13] and another in Adelman & Odlyzko [1]. Another feature we plan to add is an irreducibility test for polynomials with integral coefficients over ℚ: see Childs [9, Part II, Chapter 8] for some simple such tests.

A.3 Equations of Special Form

There is one last important remark: the previous computations might destroy any “structure” of the equation: for example, if p(x) = x^{7} − 1, we detect the rational root x = 1, but after this we are left with an equation of degree 6 which we cannot solve (actually p is a cyclotomic polynomial and it is known that it has exactly two irreducible factors over ℚ, one being x−1). In this case, and more generally when q(x) = r(x^{n}) for some r ∈ ℚ[x] and integral n ≥ 2, probably an alternative approach would work better: before giving up, go back to the original equation, and check whether it has any special shape. It is quite easy to recognize polynomials of this shape. In particular, in our example above it is possible to transform the original equation into another (actually y − 1 = 0) in the variable y = x^{7}. After solving this equation, we compute the complex 7-th roots of all the values y. We perform this task by taking just one 7-th root, and then multiplying it by the 7-th roots of unity, which are given by the Euler-de Moivre formula cos(\frac{2kπ}{7} ) + i sin(\frac{2kπ}{7} ) for k = 0,…, 6. It might even be more efficient to check this at once, before embarking on any computations at all^{5}.

A.4 Solving the General Polynomial Equation of Degree 3

Suppose that the polynomial p(x) = ax^{3} + bx^{2} + cx + d has no repeated roots, and that it has neither the root 0 (in other words, d ≠ 0) nor “small” rational roots. The formula for solving the equation of degree 3 needs the computation of the discriminant Δ of the polynomial, and we take different routes according to the sign of Δ. Notice that Δ = 0 if and only if p has repeated roots, which is not our case. It is well known that the roots are all real if and only if Δ ≤ 0, and that there is a pair of complex-conjugate roots otherwise (see Appendix A.6 for the general definition of the discriminant).

In order to solve

x^{3} + a_{1}x^{2} + a_{2}x + a_{3} = 0,  (A.1)

we set

R =^{def} \frac{1}{54}(9a_{1}a_{2} − 27a_{3} − 2a_{1}^{3}), Q =^{def} \frac{1}{9}(3a_{2} − a_{1}^{2}).

5The observation in this paragraph is in general valid, but we observe that, in the context of the computation of roots of characteristic equations in order to solve recurrences, this could become useless when the order reduction method is applied.

The discriminant is Δ =^{def} Q^{3} + R^{2}. The solutions of (A.1) are of the form

−\frac{1}{3}a_{1} + (R + Δ^{1/2})^{1/3} + (R − Δ^{1/2})^{1/3}.

This yields 9 values, 3 for each possible choice of the complex cube root, but only three actually give roots of (A.1). Here and below we denote by x^{1/3} the set of all possible solutions of the equation z^{3} = x, that is all complex cube roots of the complex number x, and the same applies to x^{1/2}. On the other hand, we denote by \sqrt{x}, \sqrt[3]{x} the “arithmetical” square and cube roots, that is roots of positive quantity that do not give rise to any ambiguities.

If Δ < 0, then Q < 0 and a direct computation aided by some trigonometrical identities yields

x_{1} =^{def} −\frac{1}{3}a_{1} + 2\sqrt{−Q} cos(\frac{1}{3}θ),
x_{2} =^{def} −\frac{1}{3}a_{1} + 2\sqrt{−Q} cos(\frac{1}{3}(θ + 2π)),
x_{3} =^{def} −\frac{1}{3}a_{1} + 2\sqrt{−Q} cos(\frac{1}{3}(θ + 4π)),

where θ =^{def} arccos(−R/(Q\sqrt{−Q})), so that all roots are real. If Δ > 0 we set

\begin{cases}S =^{def} \sqrt[3]{R + \sqrt{Δ}} \\ T =^{def} \sqrt[3]{R − \sqrt{Δ}} \end{cases}  and  \begin{cases}t_{1} =^{def} −\frac{1}{3}a_{1} − \frac{1}{2}(S + T) \\ t_{2} =^{def} \frac{1}{2}(S − T)i\sqrt{3} \end{cases}

(here we take the arithmetical cube root so that S and T are real numbers), where i^{2} = −1. Finally, we compute the roots of equation (A.1) by means of

x_{1} =^{def} −\frac{1}{3}a_{1} + \frac{1}{2}(S + T),
x_{2} =^{def} t_{1} + t_{2},
x_{3} =^{def} t_{1} − t_{2}.

Notice that x_{1} is real. Notice also that, in the procedure outlined above, we only actually take arithmetical square and cube roots. For the details, we refer to Weisstein [23], or Childs [9, Part II, Chapter 3].

A.5 Solving the General Polynomial Equation of Degree 4

We briefly sketch an algorithm in order to solve the equation x^{4} + a_{1}x^{3} +

a_{2}x^{2} + a_{3}x + a_{4} = 0, exposed in http://www.1728.com/quartic2.htm. We

set

f =^{def} a_{2} − (3a_{1}^{2}/8),
g =^{def} a_{3} + (a_{1}^{3}/8) − (a_{1}a_{2}/2),
h =^{def} a_{4} − (3a_{1}^{4}/256) + (a_{1}^{2}a_{2}/16) − (a_{1}a_{3}/4).

We first find the solutions y_{1}, y_{2}, y_{3} of the third-degree polynomial equation y^{3} + (f/2)y^{2} + ((f^{2} − 4h)/16) y − g^{2}/64 = 0 by means of the algorithm described in the previous section. We observe that if g = 0 then 0 is a root of the cubic equation and in this case we can obtain the solution more efficiently calling the formula for equation of the second degree on y^{2} + (f/2)y + ((f^{2} − 4h)/16) = 0.

Let p and q be two non-zero roots of the cubic equation; the solutions of the original equation are

x_{1} =^{def} p + q + r − s,
x_{2} =^{def} p − q − r − s,
x_{3} =^{def} −p + q − r − s,
x_{4} =^{def} −p − q + r − s.

A.6 The Discriminant and the Resultant

We adapt some definitions from Childs [9, Part III, Chapter 15]. See also Mignotte and Ştefănescu [19, §1.5.2].

Definition A.4 (Resultant.) Let f, g ∈ ℂ[x] be polynomials of degrees m and n respectively, and roots α_{1}, …, α_{m}, and β_{1}, …, β_{n} respectively (with the convention that multiple roots are repeated according to multiplicities). Let a and b be their leading coefficients. The resultant of f and g is defined by

R(f, g) =^{def} a^{n}b^{m} ∏_{i=1}^{m} ∏_{j=1}^{n}(α_{i} − β_{j}).

A few remarks: it can be proved that if f and g have integer coefficients, then R(f, g) is an integer. It is also clear that R(f, g) = 0 if, and only if, f and g have a common root. There are several ways to compute R(f, g) without knowing the roots of the polynomials f and g. One of the possibilities is to consider an algorithm based on the following properties:

• R(g, f) = (−1)^{deg(f) deg(g)}R(f, g);

• R(f, g) = a^{deg(g)−deg(r)}R(f, r) if g = fq + r for polynomials q and r. Here a is the leading coefficient of the polynomial f;

• R(f, b) = b^{deg(f)} if b is a scalar.

The second property shows that we can use Euclid’s algorithm to compute the resultant. Another possibility is to compute the resultant as the determinant of the Sylvester matrix (see Weisstein [23]): this is useful when we deal with parametric polynomials, since the method based on Euclid’s algorithm needs the computation of division with quotient and remainder, and therefore the partial results may be rational functions of the parameter. We need parametric polynomials in the actual practice when using Gosper’s algorithm (see [5] and Petkovšek, Wilf and Zeilberger [20, Chapter 5]).

Definition A.5 (Discriminant.) Let f ∈ ℂ[x] be a polynomial of degree n. We define the discriminant of f by means of

D(f) =^{def} (−1)^{n(n−1)/2}R(f, f′).

We remark that D(f) = 0 if, and only if, f and f′ have a common factor, and this happens if, and only if, f has a repeated root.

A.7 Parametric Equations

In the discussion in the last few sections we assumed that the coefficients of the polynomial equations under consideration are integers or rational numbers, but it would be very interesting to allow these coefficients to contain one or more unspecified parameters. Solving parametric polynomial equations of degree 2 or more, though, can be very difficult since we may have to take into account many possible cases, as we now explain.

We saw in Section A.6 that we can associate to every polynomial equation the discriminant d, which is an integer (resp. real or complex number) if the coefficients are integers (resp. real or complex). It can be proved that a polynomial equation has a multiple root if and only if d = 0: the discussion in Section 3 shows that we have to treat differently the case when the characteristic equation has only simple roots from the case of multiple roots. Not only this, we also want to handle possible complex roots with some care. When the degree of the characteristic equation (that is, the order of the recurrence) is 1, there can be no multiple roots; when it is 2, the sign of the discriminant can be used to distinguish among the various cases, and when the degree is 3 we can still meet both requirements (see Appendix A.4).

But when the degree exceeds 3, things become rapidly more difficult: if we were to allow parametric coefficients, we should be able to identify all possible situations for the multiplicities of the roots. As an example, if the degree is 5 and we know that d = 0, the possible multiplicities for the roots are:

5;  4, 1;  3, 2;  3, 1, 1;  2, 2, 1;  2, 1, 1, 1.

If the degree is k, the number is exactly p(k) − 1, where p(k) denotes the number of the unrestricted partitions of the integer k. There is a clash of notation with the remainder of this document, but the notation p for the number of partitions is well-established, and will not occur elsewhere. It is known that p grows very rapidly: indeed

p(n) ∼ \frac{e^{K\sqrt{n}}}{4n\sqrt{3}} as n → ∞,

where K = π(2/3)^{1/2} (see Apostol [2, §14.6]). This means that we should give conditions on the parameters for each one of these possible occurrences, and that would be too difficult.

A.8 Solving a Set of Linear Equations

In order to compute the coefficients α in the general solution (3.5) is necessary to solve systems with n linear equations and n unknowns. For this is possible to use the well-known method of inverse matrix, which computes the inverse of a square matrix, where the elements of the matrix are the coefficients of the recurrence relation.

B Exact Summation of Some Special Functions

The aim of this section is the solution of this problem: given a polynomial p with complex coefficients, and a complex number z, compute an elementary closed formula for the expression

∑_{k=0}^{N} p(k)z^{k}.  (B.1)

In the last subsections, we also deal with a more general problem with the elementary trigonometric functions, so that we essentially solve the problem above for any element of the set A defined in Section 3.

B.1 Exponentials

This is by far the easiest case: it all depends on the well-known formula

∑_{k=0}^{N} z^{k} = \frac{1 − z^{N+1}}{1 − z}  (B.2)

which is valid for any complex z ≠ 1 and any integer N ≥ 0.

B.2 Products of a Polynomial and an Exponential

Here we show how to give a closed formula for (B.1) where p is a polynomial and z ≠ 1 is a complex number. It all hinges on the following important result.

Lemma B.1 Given a polynomial p ∈ ℂ[n] of degree d, there exist complex numbers b_{0}, …, b_{d} such that

p(n) = ∑_{k=0}^{d} b_{k}n_{(k)}  (B.3)

where n_{(k)} is the falling factorial function defined by

n_{(k)} =^{def} k!\binom{n}{k} = n(n − 1) ⋯ (n − k + 1).

In the language of linear algebra, the set { n_{(k)} | k ≥ 0 } is a basis of the ring ℂ[n]. If p ∈ ℝ[n] (resp. ℤ[n]), then the coefficients b_{j} belong to ℝ (resp. ℤ). Proof. We give the proof in the form of an algorithm for computing the coefficients b_{k}, and we show that it needs at most d iterations: in fact, assume that

p(n) = a_{d}n^{d} + ⋯ + a_{0},

where a_{d} ≠ 0. Now set p_{d}(n) =^{def} p(n), b_{d} =^{def} a_{d} and define inductively the polynomial p_{j} and the complex number b_{j} for j = d − 1, …, 0 by means of

p_{j}(n) =^{def} p_{j+1}(n) − b_{j+1}n_{(j+1)},
b_{j} =^{def} coeff_{j}(p_{j}),

where coeff_{j}(q) denotes the coefficient of degree j of the polynomial q. This gives the desired proof because by construction the polynomial p_{j} has degree at most j. □

Actually, we use a quite different method, which appears to be faster: we note that the polynomial P (z) = z_{(k)} vanishes for z = 0, 1, …, k − 1. This implies that b_{0} = p(0), and that p(z) − b_{0} has the factor z. Iterating, we find that b_{1} is (p(z) − b_{0})/z evaluated at z = 1, and so on. In practice, we construct a sequence of polynomials q_{k}(z) as follows: we set

\begin{cases}q_{0}(z) = p(z), \\ b_{0} = p_{0}(0), \end{cases}  and  \begin{cases}q_{k+1}(z) = \frac{q_{k}(z) − b_{k}}{z − k} , \\ b_{k+1} = q_{k+1}(k + 1). \end{cases}

Once p has been written in the form (B.3), we see that (B.1) becomes

∑_{k=0}^{N} p(k)z^{k} = ∑_{k=0}^{N} ∑_{j=0}^{d} b_{j}k_{(j)}z^{k} = ∑_{j=0}^{d} b_{j} ∑_{k=0}^{N} k_{(j)}z^{k}.  (B.4)

The only remaining problem is the summation in closed form of single summands in (B.4), far right. The following Lemma B.2, whose proof is a simple verification and is therefore omitted, is the tool that can be used.

Lemma B.2 For z ∈ ℂ with z ≠ 1 we have the identity

∑_{k=0}^{N} k_{(j)}z^{k} = ∑_{k=0}^{N} k(k − 1) ⋯ (k − j + 1)z^{k} = z^{j} \frac{d^{j}}{dz^{j}} ∑_{k=0}^{N} z^{k} = z^{j} \frac{d^{j}}{dz^{j}} \frac{1 − z^{N+1}}{1 − z} .

(B.5)

B.3 Polynomials

If z = 1 the above idea fails: we suggest several different solutions. The first, more straightforward one, depends upon the following result.

Theorem B.3 For any k ∈ ℕ there exists a polynomial p_{k} ∈ ℚ[n] of degree k + 1 such that for any N ∈ ℕ

∑_{n=1}^{N} n^{k} = p_{k}(N).

This is Theorem 1 of Levy [16], and we give a sketch of the proof below. See also the papers by Burrows & Talbot [6], and Edwards [11], for related results, that also bear on the algorithms. The most familiar instances of this fact are the formulæ

∑_{n=1}^{N} n = \frac{1}{2}N(N + 1),

∑_{n=1}^{N} n^{2} = \frac{1}{6}N(N + 1)(2N + 1),  (B.6)

∑_{n=1}^{N} n^{3} = (\frac{1}{2}N(N + 1))^{2} .

Once we know that it exists, the polynomial p_{k} can be determined by imposing k + 2 conditions on its k + 2 unknown coefficients (actually, one can prove that the leading coefficient is (k + 1)^{−1}, the next one is \frac{1}{2}, and all other coefficients of the monomials x^{j} where j has the same parity as k, as well as the constant term vanish), and solving the arising system of linear equations. For more details see Levy [16]. A rather different method follows from Appendix B.2, and it also yields a proof of Theorem B.3. Proof of Theorem B.3. We first prove by induction on N that for any k ∈ ℕ we have

∑_{n=1}^{N} n_{(k)} = \frac{(N + 1) · N · (N − 1) ⋯ (N − k + 1)}{k + 1}.  (B.7)

This is basically a simple property of binomial coefficients: indeed, it is equivalent to

∑_{n=1}^{N} \binom{n}{k} = \binom{N + 1}{k + 1}.

Hence the decomposition of the polynomial p provided by (B.3) yields the desired result, since we can write the polynomial p(n) = n^{k} as a linear combination with integral coefficients of the polynomials n_{(d)} for d = 0, …, k, and compute ∑_{n=1}^{N} n_{(d)} by means of (B.7), which gives a polynomial in N of degree d + 1. □

We now describe two more approaches to this problem: the first one depends on the computation of the Bernoulli numbers (see for example Apostol [2, §12.12] or Hardy & Wright [12, §7.9]), and the second is an iterative method described in Levy [16].

Definition B.4 (Bernoulli numbers.) The Bernoulli numbers, denoted by B_{n} for n ∈ ℕ, are the coefficients in the Taylor series development

\frac{z}{e^{z} − 1} = 1 − \frac{1}{2}z + \frac{B_{1}}{2!} z^{2} − \frac{B_{2}}{4!} z^{4} + \frac{B_{3}}{6!} z^{6} − ⋯ + (−1)^{k−1} \frac{B_{k}}{(2k)!}z^{2k} + ⋯

which is valid in the circle \{ z ∈ ℂ | |z| < 2π \}. In particular, B_{1} = \frac{1}{6}, B_{2} = \frac{1}{30}, B_{3} = \frac{1}{42}.

Theorem B.5 Setting β_{0} =^{def} 1, β_{1} =^{def} −\frac{1}{2}, β_{2k} =^{def} (−1)^{k−1}B_{k}, β_{2k+1} =^{def} 0 for k ∈ ℕ \ {0}, where the B_{k} are the Bernoulli numbers, we have

∑_{n=1}^{N−1} n^{k} = ∑_{r=0}^{k} \frac{1}{k + 1 − r} \binom{k}{r}N^{k+1−r}β_{r}.  (B.8)

This is formula (7.9.1) of Hardy & Wright [12]. The proof is given comparing coefficients of suitable Taylor series expansions. We remark that the left hand side of (B.8) contains N − 1 summands, while the right hand side is a polynomial in N of degree k + 1 and therefore contains at most k+2 summands. The number β_{n} can be computed by means of the following iterative formula: for n ≥ 2

β_{n} = ∑_{k=0}^{n} \binom{n}{k}β_{k}.

This is Apostol [2, Theorem 12.15]. We remark that there is no established notation for the Bernoulli numbers: in particular, Apostol calls B_{n} the quantities that we call β_{n}. We also remark that it is quite important, for our purposes, that the β_{n} are rational numbers.

The iterative method referred to above depends on another result in Levy [16, Theorem 2].

Theorem B.6 For each k ∈ ℕ \ {0} there is a rational number C_{k} such that, in the notation of Theorem B.3, we have

p_{k}(N) = ∑_{n=1}^{N} n^{k} = k ∫_{0}^{N} p_{k−1}(t) dt + C_{k}N.

The value of C_{k} can be computed putting N = 1, and it is known that C_{k} = 0 whenever k > 1 is odd (see Levy [16, Corollary 3]). We remark that the integral can be easily computed using a Computer Algebra System, since we are only dealing with polynomials.

B.4 Some More Exact Formulæ

For completeness’s sake, we insert some more closed formulæ, which can be useful to a general recurrence relation solver. Their theoretical interest is increased by the fact that they can be used to avoid dealing with complex numbers when studying recurrences whose characteristic equation has a pair of complex conjugate roots. The first couple of formulæ is relevant in the theory of Fourier series. For θ ∉ 2πℤ we have

∑_{k=0}^{n} sin(kθ) = \frac{cos((n + 1/2)θ) − cos(θ/2)}{2 sin(θ/2)},

∑_{k=0}^{n} cos(kθ) = \frac{1}{2} + \frac{sin((n + 1/2)θ)}{2 sin(θ/2)}.

The proof by induction is straightforward. Another (trickier) proof can be given by separating the real and imaginary parts of

∑_{k=0}^{n} e^{ikθ} = \frac{1 − e^{i(n+1)θ}}{1 − e^{iθ}} = \frac{e^{i(n+1)θ/2}}{e^{iθ/2}} · \frac{e^{i(n+1)θ/2} − e^{−i(n+1)θ/2}}{e^{iθ/2} − e^{−iθ/2}}

= e^{inθ/2} · \frac{sin((n + 1)θ/2)}{sin(θ/2)},

where the first equality is just a special case of (B.2). These formulæ imply that both sums are, in absolute value,

≤ min(n + 1, | sin(θ/2)|^{−1}) = O(min(n + 1, ∥θ/π∥^{−1}))

where ∥x∥ denotes the distance of the real number x from the nearest integer, so that ∥x∥ =^{def} min({x}, 1 − {x}). We remark that a similar trick enables one to find the closed formula for sums of the type ∑_{k=0}^{n} λ^{k} cos(kθ), where λ is any real number. For brevity, we simply state the results: For θ ∉ 2πℤ and real λ ≠ 0 we have

∑_{k=0}^{n} λ^{k} sin(kθ) = \frac{λ^{n+2} sin(nθ) − λ^{n+1} sin((n + 1)θ) − λ sin θ}{λ^{2} − 2λ cos θ + 1} ,

∑_{k=0}^{n} λ^{k} cos(kθ) = \frac{λ^{n+2} cos(nθ) − λ^{n+1} cos((n + 1)θ) + 1 − λ cos θ}{λ^{2} − 2λ cos θ + 1} .

Combining these identities with the technique developed in Appendix B.2, one can obtain exact formulæ for sums of the type ∑ p(k)λ^{k} cos(kθ), where p is the product of a polynomial and an exponential. We remark that, because of the decomposition (B.3), we need only consider polynomials of the form p(k) = k_{(j)} for some fixed non negative integer j. Next, we note that an identity similar to (B.5) with x replaced by λe^{iθ}, implies that

∑_{k=0}^{n} k_{(j)}λ^{k}e^{ikθ} = λ^{j} \frac{∂^{j}}{∂λ^{j}} ∑_{k=0}^{n} λ^{k}e^{ikθ}.

Separating the real and imaginary parts by means of the above formulæ, we obtain

∑_{k=0}^{n} k_{(j)}λ^{k} sin(kθ) = λ^{j} \frac{∂^{j}}{∂λ^{j}} \frac{λ^{n+2} sin(nθ) − λ^{n+1} sin((n + 1)θ) − λ sin θ}{λ^{2} − 2λ cos θ + 1} ,(B.9)

∑_{k=0}^{n} k_{(j)}λ^{k} cos(kθ) = λ^{j} \frac{∂^{j}}{∂λ^{j}} \frac{λ^{n+2} cos(nθ) − λ^{n+1} cos((n + 1)θ) + 1 − λ cos θ}{λ^{2} − 2λ cos θ + 1} . (B.10)

B.5 How To Avoid the Use of Complex Numbers

Here we describe how to use the results of the previous section to avoid the use of complex numbers altogether in some interesting cases. In particular, we refer to (3.13) when λ_{1} and λ_{2} are a pair of complex conjugate numbers, with ℑ(λ_{1}) > 0. In this case we write ρ =^{def} |λ_{1}| = |λ_{2}| and choose θ ∈ (0, π) such that λ_{1} = ρe^{iθ}; since λ_{2} = \overline{λ_{1}}, we also have λ_{2} = ρe^{−iθ}. We remark that if λ_{1} = α + iβ, say, then ρ = \sqrt{α^{2} + β^{2}} and cos θ = αρ^{−1}, so that both ρ and θ can be computed staying safely within the set of real numbers. Our definitions imply that

∑_{k=0}^{n} \frac{λ_{1}^{n+1−k} − λ_{2}^{n+1−k}}{λ_{1} − λ_{2}} p(k) = ∑_{k=0}^{n} ρ^{n+1−k} \frac{e^{i(n+1−k)θ} − e^{−i(n+1−k)θ}}{ρe^{iθ} − ρe^{−iθ}} p(k)

= ∑_{k=0}^{n} ρ^{n−k} \frac{2i sin((n + 1 − k)θ)}{2i sin θ} p(k)

= \frac{ρ^{n}}{sin θ} ∑_{k=0}^{n} ρ^{−k} sin((n + 1 − k)θ) p(k).

If p is the product of a polynomial and an exponential, the trigonometrical identity for sin(α−β) transforms the above formula into a linear combination of (B.9) and (B.10), for a suitable value of λ.

C Examples

Here we discuss some concrete examples: we assume familiarity with the topics in Appendix B.2. In particular, here we omit details pertaining to the computation of a closed formula for expressions of the form (B.1). We omit the initial conditions, unless we need them explicitly.

C.1 Recurrences of Order 1

Example C.1 Solve

x_{n} = x_{n−1} + p(n).  (C.1)

The procedure outlined in Section 3 yields g_{n} = α_{1} · 1^{n} = α_{1}, since the characteristic equation is λ = 1. We remark that the general solution of (C.1) is indeed

x_{n} = x_{0} + p(1) + p(2) + ⋯ + p(n) = x_{0} + ∑_{k=1}^{n} p(k).

With the notation of (3.4), g_{n} = α_{1} corresponds to the constant x_{0}. The possibility of finding an exact, closed formula for the right hand side depends heavily on the function p. If p belongs to the set A defined in Section 3.2, we can find this closed expression (see Appendix B for details). In most cases, though, even if p is a rather “simple” function such as p(n) = 1/n, p(n) = \sqrt{n} or p(n) = log n, this is not possible and we can only find approximate solutions by means of suitable summation formulæ. We will devote a future paper of this series to the problem of approximate summation: see [4].

Example C.2 Solve the first-order recurrence^{6}

x_{n} = 2x_{n−1} + 2^{n−1} + n + 1.  (C.2)

6This is taken from the complexity analysis of a computer program.

Here we have α =^{def} 2 and p(n) =^{def} 2^{n−1}+n+1, which we split as p_{1}(n) =^{def} 2^{n−1} and p_{2}(n) =^{def} n + 1. We call the resulting recurrences y_{n} and z_{n} respectively. We compute the solution for p_{1} by means of (3.8), which yields

y_{n} = ∑_{k=1}^{n} 2^{n−k} · 2^{k−1} = ∑_{k=1}^{n} 2^{n−1} = n2^{n−1}.  (C.3)

The trick explained in Appendix B.2 shows that it is convenient to write

z_{n} = ∑_{k=1}^{n} 2^{n−k}(k + 1) = 2^{n} ∑_{k=1}^{n} k2^{−k} + 2^{n} ∑_{k=1}^{n} 2^{−k}

and we find (simplifying somewhat and omitting details)

z_{n} = 2^{n} (−n(\frac{1}{2})^{n} − 2(\frac{1}{2})^{n} + 2) + 2^{n} (2(1 − (\frac{1}{2})^{n+1}) − 1)
= (−n − 2 + 2^{n+1}) + (2^{n} − 1) = 3 · 2^{n} − n − 3.

Summing the three solutions, we finally get

x_{n} = x_{0}2^{n} + n2^{n−1} + 3 · 2^{n} − n − 3.

An easy induction shows that this is, indeed, the solution of our problem.

Example C.3 The same procedure solves the simpler recurrence^{7}

x_{n} = 2x_{n−1} + 2^{n}.

We just sketch the argument: the above discussion shows that the solution has the form x_{n} = x_{0}2^{n} +y_{n}, where y_{n} is the solution of the same recurrence with y_{0} = 0. A computation similar to that in (C.3) gives y_{n} = n2^{n}, so that x_{n} = (x_{0} + n)2^{n}, as an induction proves immediately.

Example C.4 Solve

x_{n} = αx_{n−1} + n^{2}.

Corollary 3.4 implies that the case α = 1 is rather special. By (3.8) we need to find a closed formula for

∑_{k=1}^{n} α^{n−k}k^{2} = α^{n} ∑_{k=1}^{n} α^{−k}k^{2}.

If α = 1, a direct application of (B.6) from Appendix B.3 yields the answer

x_{n} = x_{0} + \frac{1}{6}n(n + 1)(2n + 1).

7This is also taken from the complexity analysis of a computer program.

If α ≠ 1, the technique described in Appendix B.2 suggests to write

∑_{k=1}^{n} α^{−k}k^{2} = ∑_{k=1}^{n} k(k − 1)α^{−k} + ∑_{k=1}^{n} kα^{−k},

so that the answer is

∑_{k=1}^{n} α^{−k}k^{2} = [z^{2} \frac{d^{2}}{dz^{2}} \frac{1 − z^{n+1}}{1 − z} ]_{z=1/α} + [z \frac{d}{dz} \frac{1 − z^{n+1}}{1 − z} ]_{z=1/α}

= [\frac{z + z^{2} − (n + 1)^{2}z^{n+1} + (2n^{2} + 2n − 1)z^{n+2} − n^{2}z^{n+3}}{(1 − z)^{3}} ]_{z=1/α}

= \frac{α^{2} + α + α^{−n}(−n^{2} + (2n^{2} + 2n − 1)α − (n + 1)^{2}α^{2})}{(α − 1)^{3}}.

The closed formula follows at once:

x_{n} = x_{0} · α^{n} + \frac{(α^{n+2} + α^{n+1} − n^{2} + (2n^{2} + 2n − 1)α − (n + 1)^{2}α^{2})}{(α − 1)^{3}}.

C.2 Recurrences of Order 2

Example C.5 Another interesting case is the famous sequence of the Fibonacci numbers: here and henceforward we write f_{n} for the solution of the second order recurrence

\begin{cases}f_{n} = f_{n−1} + f_{n−2}, & for n ≥ 2; \\ f_{0} = 0, \\ f_{1} = 1. \end{cases}

The characteristic equation for this recurrence is

λ^{2} = λ + 1

with roots λ_{1} = \frac{1}{2}(1 + \sqrt{5}) and λ_{2} = \frac{1}{2}(1 − \sqrt{5}) = −λ_{1}^{−1}. The general solution of any such recurrence is thus

αλ_{1}^{n} + βλ_{2}^{n},

but, since f_{0} = 0 and f_{1} = 1, we have

\begin{cases}α + β = 0, \\ αλ_{1} + βλ_{2} = 1, \end{cases}  whence  \begin{cases}α = 5^{−1/2} \\ β = −5^{−1/2} \end{cases}

and the n-th Fibonacci number is

f_{n} = \frac{1}{\sqrt{5}} ((\frac{1 + \sqrt{5}}{2})^{n} − (\frac{1 − \sqrt{5}}{2})^{n})

(this is known as Binet’s formula). We observe a rather general phenomenon: here |λ_{1}| > |λ_{2}| and hence the value of f_{n} is quite close to

\frac{1}{\sqrt{5}}(\frac{1 + \sqrt{5}}{2})^{n}.

Indeed, but this is rather special, since |λ_{2}| < 1, f_{n} is the integer nearest the above expression. More generally, we remark that even if we can not find all roots of the characteristic equation, for the asymptotic analysis of the behavior of the solution of the recurrence it is sufficient to find the ones with the largest modulus. We will develop this topic in a future paper dealing with approximate solutions: see [4].

It is convenient to note another general phenomenon: if x_{n} satisfies x_{n} = x_{n−1} + x_{n−2}, then x_{n} = f_{n}x_{1} + f_{n−1}x_{0}, as can be seen by induction. In other words, the solution with x_{0} = 0 and x_{1} = 1 is quite important, and indeed it is closely allied to the fundamental one as defined in (3.11), which is g_{n} = f_{n+1}.

We want to see how things change if the characteristic equation has a double root.

Example C.6 Solve

\begin{cases}x_{n} = αx_{n−1} + βx_{n−2}, & for n ≥ 2; \\ x_{0} = 0, \\ x_{1} = 1, \end{cases}

where α^{2} + 4β = 0, with β ≠ 0.

We now look for a solution of the form

x_{n} = (a + bn)λ^{n},

where λ satisfies λ^{2} = αλ + β = αλ − \frac{1}{4}α^{2}. Since x_{0} = 0 and x_{1} = 1 we have

\begin{cases}a = 0, \\ (a + b)λ = 1, \end{cases}  whence  \begin{cases}a = 0, \\ b = λ^{−1}. \end{cases}

Thus the solution is

x_{n} = nλ^{n−1} = n(\frac{α}{2})^{n−1},

as can be checked easily.

Example C.7 Solve

\begin{cases}x_{n} = 2x_{n−1} − 2x_{n−2}, & for n ≥ 2; \\ x_{0} = 0, \\ x_{1} = 1. \end{cases}

This is a nice example of a recurrence of order 2 whose characteristic equation has complex roots. Its first few values are 0, 1, 2, 2, 0, −4, −8, −8, 0, 16, 32, 32, 0, … The procedure described in Section 3 yields

x_{n} = \frac{1}{2i}((1 + i)^{n} − (1 − i)^{n}).

Since 1 + i = \sqrt{2} · e^{iπ/4}, and 1 − i = \sqrt{2} · e^{−iπ/4}, we can also write

x_{n} = 2^{n/2} \frac{e^{inπ/4} − e^{−inπ/4}}{2i} = 2^{n/2} sin \frac{nπ}{4} .

Both solution are correct, though the latter can probably be considered easier to understand, as it conveys more readily the relevant informations about x_{n} and avoids complex numbers. The oscillations in absolute value are due to the presence of two complex conjugate roots.

Example C.8 This recurrence arises from a simple combinatorial problem.

x_{n} = x_{n−1} + x_{n−2} + 1.  (C.4)

Here g_{n} = f_{n+1} is the (n + 1)-st Fibonacci number.

x_{n} = x_{1}f_{n} + x_{0}f_{n−1} + (f_{0} + f_{1} + f_{2} + ⋯ + f_{n−1}) = x_{1}f_{n} + x_{0}f_{n−1} + f_{n+1} − 1.

The last equality follows from a well-known property of the Fibonacci numbers.

Example C.9 Solve

\begin{cases}x_{n} = \frac{1}{2}x_{n−1} + \frac{1}{2}x_{n−2} + 1, & for n ≥ 2; \\ x_{0} = 0, \\ x_{1} = 0. \end{cases}

In this case the characteristic equation has a simple root λ = 1. Here g(n) = \frac{2}{3}(1 − (\frac{−1}{2})^{n+1}) (the initial conditions are g_{0} = 1 and g(1) = \frac{1}{2}) and the solution is

x_{n} = \frac{2}{3}n − \frac{4}{9} (1 − (\frac{−1}{2})^{n}) .

Example C.10 Solve

\begin{cases}x_{n} = 2x_{n−1} − x_{n−2} + 1, & for n ≥ 2; \\ x_{0} = 0, \\ x_{1} = 0. \end{cases}

In this case the characteristic equation has the double root λ = 1. According to (3.11) we set

\begin{cases}g_{n} = 2g_{n−1} − g_{n−2}, & for n ≥ 2; \\ g_{0} = 1, \\ g_{1} = 2. \end{cases}

Since λ = 1 is a double root, we look for a solution of the type g(n) = an+b for suitable a, b ∈ ℝ. The initial conditions imply a = b = 1 so that g(n) = n + 1. The presence of the root λ = 1 forces the general solution to a polynomial of degree 2 (see Corollary 3.4), and actually

x_{n} = \frac{1}{2}n(n − 1).

In the next few examples we show how to use the “guessing” procedure described in Section 3.4.1 in several cases. We recall that f_{n} denotes the n-th Fibonacci number.

Example C.11 We start with a linear polynomial. Solve the recurrence

\begin{cases}x_{n} = x_{n−1} + x_{n−2} + n, & for n ≥ 2; \\ x_{0} = 0, \\ x_{1} = 0. \end{cases}

Solving (3.15) for q we find q(n) = n + 3 so that we have

x_{n} = q(1)f_{n} + q(0)f_{n−1} − q(n) = 4f_{n} + 3f_{n−1} − n − 3.

Example C.12 Next, a quadratic polynomial. Solve

\begin{cases}x_{n} = x_{n−1} + x_{n−2} + n^{2}, & for n ≥ 2; \\ x_{0} = 0, \\ x_{1} = 0. \end{cases}

Here we find that the solution of (3.15) is q(n) = n^{2} + 6n + 13 and we have

x_{n} = q(1)f_{n} + q(0)f_{n−1} − q(n) = 20f_{n} + 13f_{n−1} − n^{2} − 6n − 13.

Example C.13 Finally, a cubic polynomial

\begin{cases}x_{n} = x_{n−1} + x_{n−2} + n^{3}, & for n ≥ 2; \\ x_{0} = 0, \\ x_{1} = 0. \end{cases}

In this case q(n) = n^{3} + 9n^{2} + 30n + 54 so that

x_{n} = q(1)f_{n} + q(0)f_{n−1} − q(n) = 94f_{n} + 54f_{n−1} − n^{3} − 9n^{2} − 30n − 54.

Example C.14 More generally, solve

x_{n} = αx_{n−1} + βx_{n−2} + 2^{n}.

We look for a function of the type q(n) = γ · 2^{n} that satisfies the functional equation (3.15). It is quite easy to see that we need

2^{n} = αγ · 2^{n−1} + βγ · 2^{n−2} − γ · 2^{n},

and this implies

γ = \frac{4}{2α + β − 4}.

Note that this fails if 2α+β = 4 (that is, if 2 is a solution of the characteristic equation), and in this case we look for a solution of the form q(n) = γn · 2^{n}. Indeed, if 2α + β = 4 we have

2^{n} = αγ(n − 1) · 2^{n−1} + βγ(n − 2) · 2^{n−2} − γn · 2^{n},

so that

γ = −\frac{2}{α + β} .

Again, this fails if we also have α + β = 0, but the two conditions together imply that α = −β = 4 so that 2 is a double root of the characteristic equation, and in this case we look for a solution of the type q(n) = γn^{2} · 2^{n}. We remark that now the recurrence is indeed

x_{n} = 4x_{n−1} − 4x_{n−2} + 2^{n}.

We want to find γ ∈ ℝ so that

2^{n} = 4γ(n − 1)^{2} · 2^{n−1} − 4γ(n − 2)^{2} · 2^{n−2} − γn^{2} · 2^{n}.

Simplifying, we find γ = −\frac{1}{2}, and the general solution is x_{n} = (1 − n)x_{0} · 2^{n} + nx_{1} · 2^{n−1} − n · 2^{n−2} + \frac{1}{2}n^{2}.

C.3 Recurrences of Higher Order

Example C.15 Solve

x_{n} = 2x_{n−1} − x_{n−2} + 2x_{n−3}.

This is more difficult example than the previous ones, since the recurrence has order 3. The characteristic equation has the roots λ_{1} = 2, λ_{2} = i, λ_{3} = −i, and the general solution is therefore

x_{n} = α · 2^{n} + β · i^{n} + γ · (−i)^{n},

for suitable complex numbers α, β, γ. The computation of these coefficients reveals that the solution has the form

x_{n} = \frac{1}{5}(x_{0} + x_{2}) · 2^{n} + \begin{cases}\frac{1}{5}(4x_{0} − x_{2}) & if n ≡ 0 mod 4, \\ −\frac{1}{5}(2x_{0} − 5x_{1} + 2x_{2}) & if n ≡ 1 mod 4, \\ −\frac{1}{5}(4x_{0} − x_{2}) & if n ≡ 2 mod 4, \\ \frac{1}{5}(2x_{0} − 5x_{1} + 2x_{2}) & if n ≡ 3 mod 4, \end{cases} (C.5)

where, clearly, the rightmost expression is periodic with period 4 and is therefore O_{x_{0},x_{1},x_{2}}(1). As above in Example C.7 we can also transform i and −i into e^{iπ/2} and e^{−iπ/2} respectively, so that i^{n} = cos(nπ/2) + i sin(nπ/2) and (−i)^{n} = cos(nπ/2) − i sin(nπ/2). This implies the following alternative formula for x_{n}:

x_{n} = \frac{1}{5}(x_{0} + x_{2}) · 2^{n} + \frac{1}{5}(4x_{0} − x_{2}) cos(\frac{nπ}{2} ) + \frac{1}{5}(5x_{1} − 2x_{0} − 2x_{2}) sin(\frac{nπ}{2} ),

which, of course, is the same as (C.5). We remark that it is possible to avoid any computations with complex numbers (see Appendix B.5), since we consider only real recurrences, and therefore the characteristic equation has pairs of complex conjugate roots, if any. Here we do not have oscillations in absolute value as in Example C.7, since there is only one root of maximal modulus, and it is real and positive: therefore the term containing 2^{n} dominates the solution, unless x_{0} + x_{2} = 0. If this is the case, the solution is x_{n} = x_{0} cos(nπ/2) + x_{1} sin(nπ/2), and is periodic with period 4, and therefore bounded.

C.4 Order Reduction

Example C.16 Solve

x_{n} = 2x_{n−2} + 1.

Applying the traditional method exposed in Section 3.4, the solution is

x_{n} = 2^{\frac{n}{2} −1}(1 + (−1)^{n})x_{0} + 2^{\frac{n}{2} − \frac{3}{2}} (1 − (−1)^{n})x_{1}

+ \sqrt{2}^{n} (\frac{1}{2} + \frac{1}{4}\sqrt{2}) + (−\sqrt{2})^{n} (\frac{1}{2} − \frac{1}{4}\sqrt{2}) − 1.

Applying the order reduction method in Section 3.6, instead, the solution has the following shape

x_{n} = x_{r} 2^{− \frac{1}{2} r+ \frac{1}{2} n} + 2^{− \frac{1}{2} r+ \frac{1}{2} n} − 1,

where r = n mod 2 ∈ {0, 1}, more readable and concise than the previous one.

Example C.17 Solve

x_{n} = 2x_{n−2} − x_{n−4} + n.

The traditional method yields the solution

x_{n} =\frac{1}{4}(2 + 2(−1)^{n} − (−1)^{n}n − n)x_{0} + \frac{1}{4}(3 − 3(−1)^{n} + (−1)^{n}n − n)x_{1}

+ \frac{1}{4}((−1)^{n}n + n)x_{2} + \frac{1}{4}(−1 + (−1)^{n} − (−1)^{n}n + n)x_{3}

+ \frac{5}{8} − \frac{5}{8}(−1)^{n} + \frac{1}{24}n^{3} + \frac{1}{4}n^{2} + \frac{7}{16}(−1)^{n}n − \frac{53}{48}n.

Applying the order reduction method we obtain the equivalent solution

x_{n} = (1 − \frac{1}{2}n + \frac{1}{2}r) x_{r} + (\frac{1}{2}n − \frac{1}{2}r) x_{r+2}

+ \frac{1}{12}r^{3} + \frac{1}{2}r^{2} + \frac{1}{24}n^{3} + \frac{1}{4}n^{2} − \frac{1}{8}r^{2}n − \frac{3}{4}rn + \frac{2}{3}r − \frac{2}{3}n,

where r = n mod 2 ∈ {0, 1}, which is again much simpler.