Effective Bounds for P-Recursive Sequences - Semantic Scholar

EFFECTIVE BOUNDS FOR P-RECURSIVE SEQUENCES MARC MEZZAROBBA AND BRUNO SALVY We describe an algorithm that takes as input a complex sequence given by a linear recurrence relation with polynomial coecients along with initial values, and outputs a simple explicit upper bound (vn ) such that |un | ≤ vn for all n. Generically, the bound is tight, in the sense that its asymptotic behaviour matches that of un . We discuss applications to the evaluation of power series with guaranteed precision.

arXiv:0904.2452v1 [cs.SC] 16 Apr 2009

Abstract.

(un )

1. Introduction A sequence

u ∈ CN

is

polynomially recursive, or P-recursive (over Q) if it satises

a non-trivial linear recurrence relation

p[s] (n) un+s + · · · + p[1] (n) un+1 + p[0] (n) un = 0

(1)

with polynomial coecients power series)

u

is

p[k] ∈ Q[n].

dierentially nite,

Likewise, an analytic function (or a formal

or

D-nite,

if it is solution to a non-trivial

linear dierential equation (2)

p[r] (z) u(r) (z) + · · · + p[1] (z) u0 (z) + p[0] (z) u(z) = 0,

p[k] ∈ Q[z].

The coecients of a D-nite power series form a P-recursive sequence, and conversely, the generating series of a P-recursive sequence is D-nite. Numerous sequences arising in combinatorics are P-recursive, while many elementary and special functions are D-nite. Starting with the works of Stanley (1980), Lipshitz (1989) and Zeilberger (1990), D-niteness relations have gradually been recognized as good

data structures

for

symbolic computation with these analytic objects. This means that many operations of interest may be performed on the implicit representation of sequences and functions provided by an equation such as (1), (2) along with suciently many initial values (see Salvy and Zimmermann, 1994; Stanley, 1999). In recent years, signicant research eorts have been aimed at developing and improving algorithms operating on this data structure. In this article, we describe an algorithm for computing upper bounds on Precursive sequences of complex numbers. Specically, we prove the following theorem (whose vocabulary is made more precise in the sequel).

Theorem 1. Given as input a reversible recurrence relation of the form (1) with rational coecients along with initial values dening a solution (un ) ∈ Q[i]N , Algo¯ ∗ (the set of positive algebraic numbers) rithm 5 computes A ∈ R+ , κ ∈ Q, α ∈ Q + Key words and phrases. Algorithm, bounds, Cauchy-Kovalevskaya majorant, certied evaluation, holonomic functions. This research has been supported in part by the joint Inria-Microsoft research laboratory. 1

2

MARC MEZZAROBBA AND BRUNO SALVY

and φ such that |un | ≤ A n!κ αn φ(n);

∀n ∈ N,

(3)

with φ(n) = eo(n) . Moreover, for generic initial values, κ and α are tight. Asymptotic expansions of P-recursive sequences are a well-studied subject (see, e.g., Odlyzko, 1995; Flajolet and Sedgewick, 2009) and their computation has been largely automated (Wimp and Zeilberger, 1985; Flajolet et al., 1991; Zeilberger, 2008). While an asymptotic estimate gives a precise indication on the behaviour of the sequence for large values of its index, it cannot in general be used to get an estimate for a specic value. Our result lets one obtain explicit bounds valid for any term, while the tightness of the bound with respect to the asymptotic behaviour implies that the bound is not straying too far away from the actual value. These bounds may be useful both inside rigorous numerical algorithms for problems such as D-nite function evaluation or numerical integration, or as standalone results to be reported to the user of a computer algebra system. The problem of accuracy control in several settings covering the evaluation of D-nite functions has been considered by many authors (see Ÿ5.2 for more on this); our main contribution from this viewpoint is to give bounds that are asymptotically tight.

Example.

To get a sense of the kind of bounds we can compute, consider the

following examples.

For readability, the constants appearing in the polynomial

parts of the bounds are replaced by low-precision approximations. (a) Suppose we want to bound

∞

Z

tn e−t

In =

2

−1/t

dt

0 as a function of

n ∈ N.

From the recurrence relation

I0 , I1 , I2 ≤ 1/5,

and the initial conditions

1/2 −n/2

In ≤ n!

2

2In+3 = (n + 2)In+1 + In

Algorithm 5 nds that

  n + 19 · (0.26 n + 0.76) . 19

In ∼ n!1/2 2−n/2−3/4 (π/n)3/4 as n → ∞, so that with the notations −1/2 Theorem 1, κ = 1/2, α = 2 are indeed recovered by our algorithm.

In fact, of

(This example and the following one are adapted from Wimp and Zeilberger (1985, Examples 2.1 and 2.3), who illustrate the computation of asymptotic expansions by the Birkho-Trjitzinsky method.) (b) The number

tn

of involutions of

{1, . . . , n}

satises the recurrence relation

t(n + 2) = (n + 1)t(n) + t(n + 1), tn ∼ (8π)−1/4 n!1/2 e

and

√

n−1/4 −1/4

n

as

t(0) = t(1) = 1,

n → ∞

(see Knuth, 1997, Ÿ5.1.4).

Assume that we wish to bound the probability that a permutation chosen uniformly at random is an involution: the same algorithm leads to

√ 1 8 t(n) ≤ (5n + 13) n!−1/2 [z n ] exp = O(n1/4 n!−1/2 e4 2n ). n! 1000 1−z Compare (Flajolet and Sedgewick, 2009, Example VIII.5). addition to the parameters type

eO(

√

n)

α

and

κ

Notice that, in

of Theorem 1, the subexponential growth

is preserved. However, our algorithm is not designed to preserve

the constant in this

O(·)

term.

EFFECTIVE BOUNDS FOR P-RECURSIVE SEQUENCES

3

Newton polygon, change of variable (Ÿ2)

Recurrence

Dierential

Normalized

Recurrence

recurrence

equation

Poincaré type

0 regular point, nite radius of convergence

yn =

Pn−1 j=0

···

Majorant series (Ÿ3)

Coecient extraction, saddle point method (Ÿ4)

Formula

|un | ≤ n!p/q αn φ(n)

Formula for

Majorant

Majorant

normalized rec.

equation

recurrence

|˜ un | ≤ αn φ(n)

order 1, simple solution

Bounds on tails... Figure 1. Outline of our bound computation method. Solid ar-

rows represent computation steps; dashed arrows indicate proof steps without counterpart in the algorithm.

(c) One of the fastest ways to compute high-precision approximations of

π

resorts

to the following formula due to Chudnovsky and Chudnovsky (1988, p. 389):

∞ X

tk =

k=0

6403203/2 12π

where

tk =

(−1)k (6k)!(13591409 + 545140134k) . (3k)!(k!)3 6403203k

Using the method of Ÿ4.2 on the obvious rst order recurrence relation satised by

(tk ),

our algorithm leads to

∞ X tk ≤ 106 (2.3n3 + 13.6 n2 + 25n + 13.6)αn k=n

where

α=

gives about

1 −14 . We see that each term of the series 151931373056000 ' 0.66 · 10 14 more correct decimal digits of π , and we can easily deduce a

suitable truncation order to compute

π

to any given precision.

(d) Similarly, from the dierential equation

z Si000 (z) + 2 Si00 (z) + z Si0 (z) = 0,

Si(0) = 0, Si0 (0) = 1

the result of our algorithm shows that the Sine integral special function may be approximated with absolute error less than

10−100

on the disk

|z| ≤ 1

by

truncating its Taylor series at origin to the order 73. Outline. Our approach is summarized in Figure 1. Consider a solution

(un ) of Equa-

tion (1). Classical methods involving Newton polygons and characteristic equations allow to extract from the recurrence relation some information on the asymptotic behaviours that

(un )

may assume. We use these methods to factor out the main

asymptotic behaviour, thus reducing the computation of a bound on

|un |

to that

of a bound on a sequence of subexponential growth. This sequence is solution to a normalized recurrence computed in that step. Using the correspondence between

4

MARC MEZZAROBBA AND BRUNO SALVY

P-recursive sequences and D-nite functions, we encode this sequence by a dierential equation satised by its generating function (Ÿ2). Then we adapt the method of Cauchy-Kovalevskaya majorant series to bound this generating function.

The

key point here, in view of the requirement of asymptotic tightness, is to nd a majorant whose disk of convergence extends to the nearest singularity of the equation, thus avoiding the loss of an exponential factor usually associated with the majorant series method (Ÿ3). We show how to deduce several kinds of explicit bounds on and

n n un z

P

un

from the asymptotic behaviour and the majorant series (Ÿ4). Finally,

we introduce our implementation of the algorithms of this article and we briey discuss their use in the context of high-precision numerical evaluation (Ÿ5). Terminology and Notations. We let

Q[n]hSi

be the algebra of recurrence operators

with polynomial coecients, viewed as noncommutative polynomials over the shift operator

S : CN → CN , (un )n∈N 7→ (un+1 )n∈N .

consider are indexed by the nonnegative integers. Similarly, ferentiation of formal power series, and

Q[z]h∂i

Q[n]

in

Note that the sequences we

∂

stands for the dif-

for the algebra of linear dierential

operators with polynomial coecients, written with

∂

on the right. Noncommuta-

tive monomials are written and represented in memory with the coecient on the

S or ∂ on the right. u ∈ C[[z]], we denote by un (or

left and the power of the main variable For any formal power series the coecient of

z

n

in

u.

sometimes by

[z n ]u)

Following van der Hoeven (2003), we also write

u;n =

∞ X

uk z k ,

un; =

n−1 X

uk z k .

k=0

k=n

To avoid ambiguity, most other indexed names are written using bracketed superscripts, like

p[0]

in Equation (1). We use the notations of Graham et al. (1989) for

n−1 n k=0 (x + k) and x = k=0 (x − k). In the statement of algorithms, we employ expressions such as set x ≥ v , to

the rising and falling factorials, namely mean compute an approximation of requirement) and assign it to

v

xn =

Qn−1

Q

by excess (without any precise accuracy

x.

2. Factorial and Exponential behaviour In this section, we collect classical results on the asymptotics of P-recursive sequences. These will both allow us to make precise statements about the tightness of the bounds we compute and serve as a guide to organise the computation in order to meet these requirements.

Moreover, we state eective versions of some

parts of the results, that constitute the rst steps of our algorithm. 2.1.

The Perron-Kreuser theorem.

A linear recurrence relation

p[s] (n)un+s + · · · + p[1] (n)un+1 + p[0] (n)un = 0, P [k] k or the corresponding operator p S , is called nonsingular when p[s] (n) 6= 0 for [0] all n ∈ N. It is called reversible when p (n) 6= 0 for all n ∈ N. [k] Assume that the coecients p (n), k = 0, . . . , s of (4) are sequences such that p[k] (n) ∼n→∞ ck ndk for some ck ∈ C, dk ∈ Z (for instance, they are rational κ functions of n). If (un ) is a solution of (4) with un+1 /un ∼n→∞ λn then for the recurrence equation to hold asymptotically, the maximum value of dk + kκ for k = 0, . . . , s must be reached at least twice, so that the corresponding terms can

(4)

EFFECTIVE BOUNDS FOR P-RECURSIVE SEQUENCES

slope

n

−κ = −1/2

n

5

u ˜n = ψn un (n + q)p ψn+q = ψn S

S

un+3 +un+2 +nun+1

(n+6)(n+4)˜ un+6 +2(n+4)(n+1)˜ un+4 −(n2 −n−5)˜ un+2 −(n+1)˜ un =0

+(n+1)un =0

Figure 2. Newton polygons of recurrence operators, before and

after normalization.

cancel. This means that

polygon

−κ

must be among the slopes of the edges of the

Newton

of the equation.

The Newton polygon of (4) is the upper convex hull of the points

k = 0, . . . , s (see Figure 2). If e is an edge of the polygon, we denote slope. If (t, dt ) is the leftmost point of e, then the algebraic equation X (5) χe (λ) = ck λk−t = 0

(k, dk ) ∈ R2 , by −κ(e) its

(k,dk )∈e is called the

characteristic equation

of

teristic equations sum up to the order

e. Observe that the s of the recurrence.

degrees of the charac-

Theorem 2

(Poincaré, Perron, Kreuser). For each edge e of the Newton polygon of (4), let λe,1 , λe,2 , . . . be the solutions of the characteristic equation χe , counted with multiplicities. (a) If for each e, the moduli |λe,1 | , |λe,2 | , . . . are pairwise distinct, then any solution (un ) that is not ultimately 0 satises un+1 /un ∼n→∞ λe,i nκ(e) for some e and i. (b) If moreover (4) is reversible, then it admits a basis of solutions (u[e,i] )e,1≤i≤deg χe such that

[e,i]

un+1

(6)

[e,i]

un

∼n→∞ λe,i nκ(e) .

(c) If there exists e and i 6= j such that |λe,i | = |λe,j |, results analogous to (a) and (b) hold with the weaker condition u[e,i] 1/n n lim sup κ(e) = |λe,i | . n→∞ n!

(7)

Denition 1

.

(Normalized Recurrences)

(i.e., if after dividing (4) by

[s]

p

the recurrence is said to be of

If all the edges have nonnegative slope

, each coecient tends to a nite limit as

Poincaré type.

In that case, we call it

n → ∞),

normalized

if

the Newton polygon has a horizontal edge. Thus a normalized operator is one whose fastest growing solution has purely exponential (as opposed to factorial) growth. Item (a) above is known as Poincaré's theorem (Poincaré, 1885); Items (b) and (c) are Perron's theorem (Perron, 1909a,b, 1921) in the case of recurrence relations of Poincaré type, and the Perron-Kreuser theorem (Perron, 1910; Kreuser, 1914) in

6

MARC MEZZAROBBA AND BRUNO SALVY

the general case. In addition to the original works, we refer to Meschkowski (1959) and Guelfond (1963) for accessible proofs of Poincaré's and Perron's theorems. Various further extensions and renements of these results are available, see, e.g., Schäfke (1965), Kooman and Tijdeman (1990), Pituk (1997), Buslaev and Buslaeva (2005), and the references therein. In other words, the Perron-Kreuser theorem states that (4) admits a basis of solutions of the form given by Theorem 2 in some neighborhood of innity. The assumption that (4) is reversible ensures that any solution near innity extends to a solution dened on the whole set of nonnegative integers. 2.2.

Dominant Singularities.

P

If

is a polynomial, we denote by

ord(ζ, P )

the

dominant roots of P those of highest multiplicity among its nonzero roots of smallest modulus. We denote by δ(P ) and multiplicity of

ordδ (P )

ζ

P.

as a root of

We call

their modulus and multiplicity, respectively. By convention, the dominant

root of a monomial is

∞.

dominant poles of a rational function the dominant dominant singularities of a dierential operator with

We call

roots of its denominator; and

polynomial coecients the dominant roots of its leading coecient. Besides standard symbolic manipulation routines, we assume that we have at our disposal a few operations on real algebraic numbers represented using the notation

δ(P ), namely δ(P ) = δ(Q)

P, Q ∈ Q[z], whether δ(P ) < δ(Q), δ(P ) > δ(Q) and a procedure to compute arbitrarily good apdefault of δ(P ). The numerical evaluation may for instance be

a function that decides, given or

proximations by

implemented based on Graee's method (see, e.g., Schönhage, 1982, Ÿ14). The comparison can be based on a symbolic-numeric approach as in (Gourdon and Salvy, 1996). More generally, most steps of Algorithms 3 and 4 involving no precise accuracy requirement may be implemented using interval arithmetic or oating-point arithmetic with careful rounding instead of symbolically.

Remark.

Although we work over

Q

all along this paper for clarity, we expect that

most results adapt without diculty to any suciently eective subeld of

C.

However, the way to perform the basic operations we assume available in this section (as well as the details of some algorithms, especially Algorithm 3 below) may dier. 2.3.

Generic Growth of the Solutions.

Let

R ∈ Q[n]hSi

be a nonsingular

R ·u = 0 (u0 , . . . , us−1 ) ∈ Cs . Accordingly, we say that an assertion is true for a generic solution of R · u = 0, or for generic initial values, if it is true for any solution u such that (u0 , . . . , us−1 ) ∈ Cs \ V where V is s a proper linear subspace of C . reversible operator of order

s.

Then any solution of the recurrence relation

is uniquely determined by its initial values

Theorem 2 implies that the factorial and exponential asymptotic behaviour of the fastest growing solutions is determined by the dominant singularities of

R.

We use Algorithm 1 to extract this asymptotic behaviour, which is in fact that of a generic solution of

R · u = 0,

as stated by Proposition 3 below.

Proposition 3 (Factorial and Exponential Growth). Write R as

Ps

k=0

Q[n]hSi and assume b b = 6 0 for some k ∈ {0, . . . , s − 1}. Algorithm (κ, Pα ) = Asympt(R ) such that for any solution (un ) of R · u = 0, u 1/n 1 n ≤α where α = , (8) lim sup κ n! δ(P n→∞ α) [k] [s]

b[k] (n)S k ∈

1 computes

EFFECTIVE BOUNDS FOR P-RECURSIVE SEQUENCES

Algorithm 1: 1

7

Factorial and exponential behaviour

Ps

[k] k k=0 b (n)S ∈ Q[n]hSi) [k] [s] deg b −deg b κ ← maxs−1 . k=0 s−k Ps [s−`] ` Pα ← `=0 bd+`κ z where d = deg b[s] return (κ, Pα )

function Asympt(

2 3 4

with equality in the generic case. Proof. The inequality follows from Theorem 2 since −κ is the slope of the rightmost edge

e

of the Newton polygon of

R

and

Pα

χe . It V = ker R ⊂ CN .

is the reciprocal polynomial of

remains to show that equality holds for generic initial values. Let

u[1] ∈ V such that u[1] 1/n n = α. lim sup κ n→∞ n!

Also by Theorem 2, there exists

u[1] , . . . , u[s]

P V . Let u = k λ[k] u[k] ∈ V . By κ 1/n construction of κ and α, we have the inequality lim sup |un /n! | ≤ α. Up to [1] extraction of a subsequence we can assume (i) that un does not vanish, (ii) that [1] 1/n → β as |un /n!κ |1/n → α and (iii) that there exists β ≤ α such that |un /n!κ | n → ∞. Then [2] [s] 1/n u u β n n [1] → , λ + λ[2] [1] + · · · + λ[s] [1] α un un This can be extended to a basis

so that

β=α

of

unless

[s]

[2]

λ[2] un + · · · + λ[s] un which does not happen for generic

[1] un [k]

λ

→ −λ[1] , 

.

Accordingly tighter results hold if the assumptions of Theorem 2(b) are fullled. 2.4.

Generating Function and Associated P Dierential Equation.

Consider s [k] k [0] [s] b S ∈ Q[n]hSi (with b , b 6= k=0 d 0). Using the θ = z dz , it is classical that the generating series u(z) Pr [k] k of u ∈ ker R cancels the associated dierential operator D = k=0 a θ ∈ Q[z]hθi 1 computed by DieqToRec (Algorithm 2) . Normalizing D , this rewrites again a nonsingular recurrence operator

Euler derivative

R=

a[r−1] r−1 a[1] a[0]  θ + · · · + θ + · u = 0. a[r] a[r] a[r] A point z0 ∈ C is a regular point of (9) if any solution u has polynomial growth O(1) u(z) = 1/ |z − z0 | as z → z0 in a sector with vertex at z0 . Regular points

(9)

encompass



θr +

ordinary points, where the equation is nonsingular and thus has analytic regular singular points. Fuchs' criterion (see,

solutions by Cauchy's theorem, and

e.g., Ince, 1956, Ÿ15.3) states that 0 is a regular point if and only if for all

k,

the

1 Actually, the classical translation of recurrence operators to dierential operators uses g = 1. The multiplication by g in our version comes from our choice to use sequences indexed by N rather than Z.

8

MARC MEZZAROBBA AND BRUNO SALVY

Algorithm 2: 1 2 3

Ps

b[k] S k ∈ Q[n]hSi) Q s g ← Π/ gcd(b[s] , Π) where Π = k=1 (n + k) Ps Ps nj S k ← g R . thus R = k=0 ckj S k (n − k)j k=0 ckj P Ps P r s−k (θ − k)j as D = k=0 a[k] θk expand k=0 j ckj z return D

function DieqToRec(R

4 5 6 7

Recurrence to normalized dierential equation

=

k=0

∈ Q[n]hSi, κ ∈ Q) (p, q) = (0, 1) if κ = 0) Pqs ˆ[k] k ˆ compute the symmetric product R = k=0 b (n)S of R p q (n + q) S − 1 . see, e.g., Stanley (1999, Ÿ6.4) ˆ) return DieqToRec(R

function Normalize(R

p/q ← κ

8

9

(in irreducible form, with

and

a[k] /a[r] of (9) is analytic at 0, while z0 6= 0 is a regular point if and only [r] if each a /a has a pole of order at most r − k in z0 . (This criterion still holds if [k] [r] the a /a are replaced by meromorphic functions.)

coecient

[k]

Lemma 4. If R is normalized (Denition 1), then the origin is a regular point of D, and the reciprocal polynomial Rec a[r] of the leading term of D is the characteristic equation of the horizontal edge of the Newton polygon of R. Proof.

Using the notations of the function DieqToRec() in Algorithm 2, let

d[k] =

[k]

deg b for all k , and m = deg g . Thus r = maxsk=0 d[k] + m. The leading term j −k −k j of θ z as an operator in θ with Laurent polynomial coecients is z θ , hence Ps [r] s−k [s] a (z) = k=0 ckr z . The condition that R is normalized translates into d = s−1 [k] [k] [s] [r] maxk=0 d , that is, d = d = r − m for some k < s. It follows that a (0) = csr 6= 0, and if e

hence

0

is a regular point by Fuchs' criterion. Finally, if

is the edge of its Newton polygon such that

κ(e) = 0,

R

is normalized

then the general

expression

χe (λ) = λ−t

X

ak,d[k] λk

[k]

d +kκ(e) =d[s] +sκ(e) (where

t

is such that

χe (0) 6= 0)

In the general case where

simplies to

χ0 (λ) = λ−t

P

d[k] =r

ak,r λk .



R needs not be normalized, we normalize it by a change

of unknown sequence preserving P-recursiveness before computing the associated dierential equation. This is described in the next proposition. Figure 2 gives an example of normalization of recurrence operators and of its action on their Newton polygons.

Proposition 5. Let

R ∈ Q[n]hSi be nonsingular, reversible, with nonzero constant coecient with respect to S . Let (p/q, Pα ) = Asympt(R) as computed by Algorithm 1, and assume that δ(Pα ) < ∞. Algorithm 2 computes P a normalized n dierential operator D = Normalize(R, p/q) that cancels u˜(z) = ∞ n=0 ψn un z for all sequences ψ and u such that (n + q)p ψn+q = ψn and R · u = 0.

EFFECTIVE BOUNDS FOR P-RECURSIVE SEQUENCES

9

The origin is a regular point of D, and the modulus of the dominant singularities of D equals δ(Pα ). Proof. Let α = 1/δ(Pα ). Let (u[1] , . . . , u[s] ) be a basis of ker R having the asymp-

lim supn→∞ |u[k] /n!p/q |1/n ≤ α for all k . Let (ψ , . . . , ψ ) be the basis of solutions to (n + q)p ψn+q = ψn corresponding [i] to the initial values ψj = δij for 0 ≤ i, j < q , where δij is the Kronecker symbol. ˆ such that for N large enough, the sequences (ψn[j] u[k] Algorithm 2 constructs R n )n≥N 1/n ˆ u)n = 0 for n ≥ N }. For all j and k , lim sup|ψn[j] u[k] form a basis of {ˆ u | (R·ˆ | ≤ α. n P [j,k] [j] [k] ˆ Assume that u ˆ= λ ψ u is solution to R · u ˆ = 0 in some neighborhood of totic behaviours given by (7). In particular

[0]

[q−1]

j,k

lim sup|un |1/n ≤ α (indeed, if  > 0, then |un | ≤ (α + )n for n large [j] [k] 1/n enough). On the other hand lim sup|ψn un | = α for at least one (j, k). Hence, ˆ is normalized and the largest modulus of a root of by Theorem 2, the operator R innity. Then

the characteristic equation associated to the horizontal edge of its Newton polygon is

α.



Applying Lemma 4 concludes the proof.

In the sequel, we will choose as normalizing sequence the solution to p p q) ψn+q = ψn given by ψn = q − q n Γ(n/q + 1)−p .

(n +

3. Subexponential Behaviour: Majorant Series Computation The results of the previous section allow us to compute the generic factorial and exponential asymptotic behaviour of solutions of a linear recurrence relation with polynomial coecients. We now turn to the computation of a bound for the remaining subexponential factor of a particular solution. 3.1.

Majorant Series and the Cauchy-Kovalevskaya Method.

The main tool

we use is a variant of the Cauchy-Kovalevskaya majorant series method, which usually serves to establish the convergence of formal series solutions to dierential and partial dierential equations, but may also be applied to obtain explicit bounds on the tails of these solutions (see also Ÿ5.2 for more on this).

Denition 2

series

of

(Majorant series)

u ∈ C[[z]],

.

and we write

A formal power series

u E v,

if

In particular, the disk of convergence of inside the disk of convergence of

v,

|un | ≤ vn v

v ∈ R+ [[z]] n ∈ N.

is a

majorant

for all

u, and if z lies |un; (z)| ≤ vn; (|z|) for all n ≥ 0.

is contained in that of

we have that

Other immediate properties of majorant series are summarized in the following lemma.

Lemma 6. Assume that

u, u[1] , u[2] ∈ C[[z]], v, v [1] , v [2] ∈ R+ [[z]] E v [2] . Then

u E v , u E v and u du dv E ; u[1] + u[2] E v [1] + v [2] ; dz dz [1]

[1]

[2]

u[1] u[2] E v [1] v [2] ;

are such that

u[2] ◦ u[1] E v [2] ◦ v [1]

where in the last inequality it is assumed that u[1] (0) = v[1] (0) = 0. In the neighborhood of an ordinary point, majorant series for the coecients of a dierential equation like (2) give rise to similar majorants for the solutions. Indeed, if

(

u(r) = a[r−1] u(r−1) + · · · + a[0] u v (r) = b[r−1] v (r−1) + · · · + b[0] v

|u(0)| ≤ v(0), . . . , |u(r−1) (0)| ≤ v (r−1) (0)

10

MARC MEZZAROBBA AND BRUNO SALVY

a[k] , b[k] are analytic functions at 0 such that a[k] E b[k] for all k , then by [k] induction u E v . This result does not hold if one of the a has a pole at 0; however, the method may be adapted to the case where 0 is a regular singular point of the where

dierential equation. We give one way to do this in Ÿ3.3; for a more complete introduction to the usual Cauchy-Kovalevskaya method in the ODE setting covering the regular singular case, see Mezzino and Pinsky (1998), and for a more general statement along these lines, see van der Hoeven (2003, Proposition 3.7).

In any

case, the rst step for obtaining majorant series for the solutions of a dierential equation using the Cauchy-Kovalevskaya method is to compute majorants for its coecients, which in the case we are interested in are rational functions.

Bounds for Rational Functions.

3.2.

rn z n , D(0) 6= 0.

P

The sequence

Consider a rational function r(z) = N (z)/D(z) = (rn ) satises a linear recurrence relation with con-

stant coecients, whose characteristic polynomial is the reciprocal polynomial of

D.

This recurrence can be solved by partial fraction decomposition of

r,

yielding

the explicit expression

ord(ζ,D) (10)

X

X

D(ζ)=0

d=1

rn =

h[ζ,d] ·(n+1)d−1 ·ζ −n ,

n ≥ max(0, deg N −deg D +1),

|rn | ≤ M δ(D)−n nordδ D . Pα and a positive −m integer m. It returns a bound of the form r(z) E M (1−αz) , where α = 1/δ(Pα ). In particular, when Pα = D and m = ordδ (D) this bound is tight. To compute a suitable M , we start with the right-hand side of (10) divided by 1 bn = [z n ] = (n + 1)m−1 · αn . (1 − αz)m

with

h[ζ,d] ∈ Q(ζ).

We are now aiming at a bound of the form

In view of later needs, Algorithm 3 takes as input a polynomial

By applying the triangle inequality, we get a sum

t(n)

of terms of the form

(n + 1)d−1 n λ (n + 1)m−1 where 0 ≤ c, 0 < λ ≤ 1, and m < d only if λ < 1. Such a term is decreasing for n ≥ 1 if d ≤ m and for n ≥ (d − m)/ log(1/λ) otherwise. We compute an index N0 starting from which the inequality |rn /bn | ≤ t(n) is guaranteed to hold and t(n) to be decreasing; then we adjust M from the explicit values of the rst N0 coecients c

and bounds on the tails. For this last part, consider the squarefree decomposition

[ζ,d]

[ζ,d]

D = D1 D22 · · · Dkk . If ζ = hi,d (ζ) · ζ −d for some

Di , then each h may in fact be written h hi,d ∈ Q[ζ] depending only on Di and d. Moreover, in this expression, −1 be bounded by δ(Di ) . Hence we have k i−1 d−1 X X rn −n X −d (n + 1) −n = α ζ hi,d (ζ)ζ bn (n + 1)m−1 i=1 Di (ζ)=0 d=0 ! k X i−1 X X hi,d (ζ) (n + 1)d−1
−n ≤ α δ(Di ) . d m−1 ζ (n + 1) i=1 d=0 D (ζ)=0

is a root of

polynomial

−1

|ζ|

(11)

may

i

We may take for approximation.

t(n)

the right-hand side of (11), or even a suitable numerical

To deal with the sum in parentheses, we may bound

ζ −d hi,d (ζ)

EFFECTIVE BOUNDS FOR P-RECURSIVE SEQUENCES

Algorithm 3: 1 2 3 4

Tight majorant series for rational functions

= N/D ∈ Q(z), Pα ∈ Q[z], m ∈ N∗ ) let A + B/D = r with A, B ∈ Q[z]. 2 k compute the squarefree factorization D = D1 D2 · · · Dk of D compute the coecients hi,d ∈ Q[ζ] of the partial fraction decomposition Pk P Pi hi,d (ζ) B(z) . See, e.g., Bronstein (2005, Ÿ2.7) i=1 Di (ζ)=0 d=1 (ζ−z)d D(z) = for i = 1, . . . , k do for d = 1, . . . , i do P −d set ci,d ≥ Di (ζ)=0 hi,d (ζ)ζ
i−m k set N0 ≥ max 1, 1 + deg A, maxi=m+1 log(δ(Di )/δ(Pα )) Pk Pi−1 n (n+1)d−1 let t(n) = i=1 d=0 ci,d (n+1)m−1 (δ(Pα )/δ(Di )) PN0 −1 compute the truncated series r;N0 (z) = rn z n n=0


N0 −1 n set h(N0 ) ≥ maxn=0 |rn | / n+m−1 m−1 δ(Pα )
return an approximation by excess of max h(N0 ), t(N0 )

function BoundRatpoly(r

5 6 7 8 9 10 11 12

term-by-term, replacing once again on the sign of of

11

Di

`.

ζ`

by

δ(Di )`

or

δ(ζ deg Di Pi (1/ζ))−`

depending

We may also simply compute low-precision enclosures of the roots

and then use interval arithmetic.

The complete procedure is summarized in Algorithm 3.

We have thus proved

the following.

Proposition 7. Given r = N/D ∈ Q(z) (in irreducible form), Pα ∈ Q[z], and m ∈ N∗ ,

such that 0 < δ(Pα ) ≤ δ(D) and δ(Pα ) = δ(D) only if m ≥ ordδ D, Algorithm 3 computes M = BoundRatpoly(r, Pα , m) ∈ Q+ satisfying r(z) E M (1−z/δ(Pα ))−m . To improve

M , we may loop over lines 10 and 11 of Algorithm N0 or t(N0 ) − h(N0 ) reaches some specied value.

3, doubling

N0

each time, until

3.3.

Bounds for D-nite Functions.

We now apply the Cauchy-Kovalevskaya

method proper to deduce a majorant series for of

(un )

u(z) from the asymptotic behaviour

obtained in Ÿ1 and majorant series for the coecients of an associated

dierential equation. The majorant series we obtain is simpler than

u(z)

in the

sense that is always satises a dierential equation of order 1. By Fuchs' criterion, we may isolate the constant term of each coecient of (9), giving (12)

Q(θ) · u = z(˜ a[r−1] θr−1 + · · · + a ˜[1] θ + a ˜[0] ) · u,

Q ∈ Q[X] is a monic polynomial of degree r and the a ˜[k] are rational functions of z . Let mk ∈ N be the maximum multiplicity of a point of the circle |z| = δ(Pα ) as a pole of a ˜[k] and let T = maxr−1 k=0 (mk − r + k). We emphasize that, although Algorithm 4 takes Pα on input, the whole point of the method is that δ(Pα ) may indeed equal the modulus of the dominant singularities of D . In that case, the integer T is sometimes called the Malgrange irregularity of these singularities (see Malgrange, 1974), and by Fuchs' criterion again, T = 0 if and only if the dominant where

12

MARC MEZZAROBBA AND BRUNO SALVY

singularities are all regular. Using Algorithm 3, we compute bounds of the form

[k]

a ˜

(13)

Mk E (1 − αz)r−k+T

i.e.,

for the coecients of the equation, with

  n+r−k+T −1 n [k] an ≤ Mk α ˜ r−k+T −1 α = 1/δ(Pα )

as usual (lines 67 of Algo-

rithm 4). Extracting the coecient of

zn

in (12), we get

Q(n) un =

(14)

n−1 r−1 XX

[k]

a ˜n−1−j j k uj .

j=0 k=0 Since

Q

is monic, let

N1

be such that

Q(n) > 0

for

n ≥ N1 ;

then by (13), for such

n, Q(n) |un | ≤

(15)

n−1 r−1 XX

 Mk

j=0 k=0

 n−1−j + r−k+T − 1 n−1−j k α j |uj | . r−k+T − 1

Lemma 8 (Reduction from order r to order 1). Let M = maxr−1 k=0 Mk / r−1 X k=0

Proof. 

 Mk

r−1 k




; then

   n−1−j+T n−1−j + r−k+T −1 k . j ≤ M nr−1 T r−k+T −1

j ≤ n − 1 and k ≤ r − 1, we have −1   n−1−j+T n−1−j + r−k+T −1 (n − j + T )r−1−k ≤ (n − j)r−1−k ; = T r−k+T −1 (T + 1)r−1−k For

thus



n−1−j+T T

−1 X r−1

 Mk

k=0

 r−1 n−1−j + r−k+T −1 k X j ≤ Mk j k (n − j)r−1−k r−k+T −1 k=0

≤ M nr−1 , 

establishing the lemma.

N2 ≥ N1 be such that M nr ≤ αKQ(n) for n ≥ N2 . (In the common case where Q(n) = nr , we may take K = M/α.) Suppose that some sequence (vn ) satises vn ≥ |un | for 0 ≤ n ≤ N2 and   n−1 n−1−j + T 1X K αn−j vj (16) vn = n j=0 T With

for all

M

as in Lemma 8, choose

n ≥ 1.

K > M/α.

Let

Using (15) and Lemma 8, we get by induction

|un | ≤ vn

for all

n ∈ N.

Now (16) translates into

v 0 (z) =

(17)

αK v(z), (1 − αz)T +1

which admits the simple solutions (18) below. Finally, we adjust the integration constant

n < N2

(lines 1112).

parameter

u;n

A

so as to ensure that

|un | ≤ vn

for

If no specic solution of (9) is given (i.e., if we drop the

of Algorithm 4) we still obtain a result valid up to some multiplica-

tive constant by simply ignoring this last part. The result of this computation is summarized in the following.

EFFECTIVE BOUNDS FOR P-RECURSIVE SEQUENCES

Algorithm 4: 1 2 3 4

Majorant series for normalized D-nite functions

function BoundNormalDieq( for

5 6 7

13

Pr

k=0

k = 0, . . . , r − 1 do c[k] ← (a[k] /a[r] )z=0 (or fail a ˜[k] ← (a[k] /a[r] − c[k] )/z

a[k] θk ∈ Q[z]hθi, Pα ∈ Q[z], u;· )

with error  0 should be a regular point)

T ← max{0; ordδ (den a ˜[k] ) − r + k | 0 ≤ k < r − 1 and δ(den a ˜[k] ) = δ(Pα )}. for k = 0, . . . , r − 1 do Mk ← BoundRatpoly(˜ a[k] , Pα , T + r − k) . thus a ˜[k] E Mk (1 − αz)−T −r+k
r−1 M ← maxr−1 k=0 Mk / k ∗ compute K ∈ N such that K ≥ 2M δ(Pα ) starting with N2 = 1, double N2 until Pr−1 [k] k N2 > (1 − M δ(Pα )/K) N2r k=0 c compute u;N2 +1 and v;N2 +1 where v is given by (18) 2 A ← maxN n=0 |un | /vn return (T, K, A)

8 9 10 11 12 13

Proposition 9. Let D ∈ Q[z]hθi, and let u;n be a function that computes truncated series expansions of a specic u ∈ ker D up to any order n. Let Pα ∈ Q[z]. Assume that 0 is a regular point of D and that the dominant singularities of D are nite and of modulus at least δ(Pα ). Then BoundNormalDieq(D, Pα , u;· ) (Algorithm 4) returns T ∈ N, K ∈ N∗ , A ∈ Q+ such that A (1 − αz)K u(z) E v(z) =  A exp K/T (1 − αz)T   

(18)

In addition to its modulus

α,

if T = 0 otherwise.

Algorithm 4 actually preserves the irregularity

T

of the dominant singularity of the dierential equation, which is connected to the subexponential growth of the coecient sequence.

Remark.

Sometimes all we need is a simple majorant series satisfying the tightness

property of Theorem 1 for the solutions of a dierential equation of the form (2) at an

ordinary

point. Instead of the results of this section, we may then apply the

plain Cauchy-Kovalevskaya method outlined in Ÿ3.1 using a majorant equation of the form

 r−1   α r−k X M r−1 r−k N v = v (k) . (1 − αz)N k 1 − αz k=0
N the majorant series v(z) = exp M/(1 − αz) . If additionally (r)

This gives

the

dominant singularity is regular, we may instead use the Euler equation

v (r) =

r−1 X k=0

yielding

Mk v (k) , (1 − αz)r−k

v(z) = A/(1 − αz)λ where αr λr − Mr−1 αr−1 λr−1 − · · · − M0 = 0. In both M , resp. Mk may be determined using Algorithm 3.

cases suitable parameters

14

MARC MEZZAROBBA AND BRUNO SALVY

Algorithm 5: 1 2 3

Bounds for general P-recursive sequences

Ps

b[k] (n)S k ∈ Q[n]hSi, [u0 , . . . , us−1 ] ∈ Q[i]s ) R ← R · S −m where m = min{k | p[k] 6= 0} (κ, Pα ) ← Asympt(R) . Normalize and encode the subexponential part by a dierential equation D ← Normalize(R, κ) . Bound the solutions of the dierential equation dene a function u ˜;· that unrolls the recurrence relation R · u = 0 starting from u0 , . . . , us−1 to compute Pn u ˜;n = k=0 q −pk/q Γ(k/q + 1)−p uk z k (where p/q = κ) for any n ∈ N (T, K, A) ← BoundNormalDieq(D, Pα , u ˜;· ) return (κ, T, Pα , K, A)

function BoundRec(R

4 5

6 7

=

k=0

4. Explicit Bounds 4.1.

P-Recursive Sequences.

quence

vn

At this point, we are able to bound

given by its generating series

v(z) = L0p,q v˜(z),

where

v˜

un

by a se-

is an explicit

series satisfying a dierential equation of the rst order, and we have denoted

L0p,q v(z) =

∞ X vn n z . ψ n=0 n

(Note that series whose coecients satisfy

recurrence relations

of the rst order,

that is, hypergeometric series, cannot serve as asymptotically tight bounds for normalized D-nite functions because the range of asymptotic beahaviours that their coecient sequences assume is not wide enough: their subexponential asymptotic growth is always polynomial.)

Proposition 10. Given as input a nonsingular reversible recurrence operator R ∈ Q[n]hSi along with initial values u0 , . . . , us−1 ∈ Q[i] dening a solution (un ) ∈ Q[i]N of R · u = 0, the function BoundRec (Algorithm 5) computes p/q ∈ Q, Pα ∈ Q[z], T ∈ N and K, A ∈ R+ such that p n p + 1 v˜n (19) ∀n ∈ N, |un | ≤ vn = q q n Γ q where v˜n is dened as in

(18)

. Additionally, for generic (u0 , . . . , us−1 ), 1/n un lim sup = 1. vn n→∞

Proof.

This follows from combining the statements of Propositions 3, 5 and 9. p

Recall that we have chosen the operator

R

ψn = q − q n Γ(n/q + 1)−p .

After Line 2 of Algorithm 5,

satises the hypotheses of Proposition 5.

computed on Line 4 cancels

u ˜(z) =

P∞

n n=0 ψn un z ,

Hence the operator

and the function

u ˜;·

D

dened on

the next line does indeed compute truncations of this series. By Proposition 9 it follows that

u ˜ E v˜ and,

multiplying the coecients by

ψn−1 ,

that

u E v.

generic initial values,

1/n 1/n un un = lim sup κ n+o(1) O(1) =1 lim sup v n! α n n→∞ n→∞ n

Finally, for

EFFECTIVE BOUNDS FOR P-RECURSIVE SEQUENCES

15



by Proposition 3.

Although this representation (19) is satisfactory for many applications, more

vn

explicit expressions for the coecients

T = 0,

are sometimes desirable. If

it is

readily seen that

v˜n = A α

(20) For

T > 0, the general coecient v˜n

n



 n+K −1 . K −1

still admits a rather complicated closed-form

expression in terms of the general hypergeometric function

F

(see Graham et al.,

1989, Ÿ5.5): one may check that

   k ∞ X K 1 Tk + n − 1 n = A αn TFT v˜n = Aα k! n T k=0

However,

v˜n

n+T T T +1 T

n+T +1 T T +2 T

··· ···



n+2T −1 K T 2T T T

may in turn be bounded by much simpler expressions without losing

the asymptotic tightness (in the sense of Theorem 1) using a simple version of the saddle point method (see, e.g., Flajolet and Sedgewick, 2009, Ÿ4.3).

Since

v˜n ≤ v˜(t)/tn . For xed n, the rightT +1 hand side is minimal for the unique tn ∈ (0; 1) such that Kαtn = n(1 − αtn ) . 1/(T +1) Asymptotically, tn satises 1 − αtn ∼ (K/n) as n → ∞. This approximation

v˜ ∈ R+ [[z]],

for any

t ∈ (0; 1/α),

we have

suits our purposes well: indeed, we set

for

n > K,

and choose



1 α

rn =

(21)

rn ∈ (0; 1/α)

1−

1   K  T +1

n

arbitrarily otherwise. For

T > 0,

we obtain

(22)

    T  1 −n  K  T +1 v˜(rn ) K n T +1 n v˜n ≤ =α 1− exp = αn exp O(nT /(T +1) ), n rn n T K and similarly

v˜n ≤ αn nK K −K (1 − K/n)−n = αn nO(1)

(23) if

T = 0. vn itself, (22) and (23) extend to bounds of the form (3), that make un = n!κ αn eo(n) apparent, by means of the following κ between ψn and n! .

Going back to

the asymptotic behaviour relation

Lemma 11. For q ∈ N \ {0} and n ≥ 3q/2, 1 = Γ(n/q + 1)p q p/q n ≤ ψn

Proof.

Since

Γ(x) q−1 Y

is increasing for

(

(2π)p/q (n/q + 1)p n!p/q , n−p/q n!p/q ,

p>0 p < 0.

x ≥ 3/2,

Γ(n/q + k/q) ≤ Γ(n/q + 1)q ≤

k=0

q−1 Y

Γ(n/q + k/q + 1).

k=0

By Gauÿ' multiplication theorem (see Abramowitz and Stegun, 1972, Formula 6.1.20)

Γ(qz) = (2π)

(1−q)/2 qz−1/2

q

q−1 Y k=0

 k Γ z+ q

(z ∈ C),

! .

16

MARC MEZZAROBBA AND BRUNO SALVY

this implies that

(2π)(q−1)/2 q n Γ(n/q + 1) (2π)(q−1)/2 (n + 1)q−1 ≤ ≤ Γ(n + 1) nq −1/2 q q−1/2 p/q

and the result follows by raising either inequality to the power of the sign of

depending on

p.



This concludes the proof of Theorem 1.

Remark.

If we are content with computing a numerical bound for one coecient

(or one tail, see next section) of a D-nite power seriesthat is, a bound for xed

n,

as opposed to a formula giving a bound as a function of

nthen majorant series

with the same radius of convergence as the coecients of the equation (and thus the method of Ÿ3.3) are not strictly necessary for the bound to become ultimately tight as

0

n approaches innity.

Consider for instance Equation (1) in the case where

ν > α with the notations of Ÿ3.3. Van der p[k] /p[r] E M (ν)/(1 − νz) for k = 0, . . . , r − 1,

is an ordinary point, and assume

Hoeven (2003, Ÿ3.5) proves that if then

C

u(z) E where

C

ν.

does not depend on

(1 −

Also assume that the majorizing procedure for ra-

M (ν)

tional functions used to compute

O nd (α/ν)n




νz)d(M (ν)+1)/νe is tight enough to ensure that

d = maxr−1 k=0 mk ).

(as is Algorithm 3, with

what reminiscent of the saddle-point method, we then choose, say,

(1 + 1/n1/(2d) )α,

M (ν) =

In a manner some-

ν = νn =

hence getting

|un | ≤ vn = αn+n This suggests that it is sensible to take

1−1/(2d)

.

ν = (1 + 1/nΘ(1/d) )α

in the algorithms of

van der Hoeven (2001, 2003). 4.2.

Tails of Power Series.

In Examples 1(c) and (d) on page 3, the sequence

for which we compute an upper bound is the tail whose coecients

un

tn = un; (1)

of a convergent series

are given by a linear recurrence relation of the form (1). In

such a case, the sequence tn is also P-recursive, but its initial values are unknown if we have in mind the evaluation of the sum of the series, these initial values are precisely what we are after.

u(z) E v(z), |un; (1)| ≤ vn; (1). To

However, if

majorant series (Ÿ3) ensure that

the general properties of avoid repeated majorant

computations when working with D-nite power series, notably in the context of numerical analytic continuation (see Ÿ5.2), we actually consider the slightly more general problem of bounding the tails

z

such that

[r]

|z| < δ(p ),

where

[r]

p

(j)

un; (z) of the j -th derivative of u at any point is the leading term of a dierential equation

with polynomial coecients annihilating

u(z).

We assume once again that we have computed

v(z) = L0p,q v˜(z)

(with

p ≤ 0,

κ = p/q

and

v˜

such that

so that the radius of convergence of

v

u(z) E

is positive)

using the algorithms of Ÿ2 and Ÿ3. The formalism of majorant series proves handy here, as we have

(j)

(j)

|un; (z)| ≤ vn; (|z|)

by Lemma 6. Notice that if

may very well lie beyond the disk of convergence of

Proposition 12. Let

p < 0,

the point

z

v˜.

be the modulus of the dominant singularities of the dierential operator dening u. If κ = 0, we assume |z| < α−1 , otherwise z my be arbitrary. As in the previous section, set rn = (|z| + α−1 )/2 when κ = 0 and α−1

EFFECTIVE BOUNDS FOR P-RECURSIVE SEQUENCES

n ≤ K/(1 − α |z|)T +1 and as in n ≥ q/2, (j) (24) un; (z) ≤

Equation

(21)

17

otherwise. Then for all n, j ∈ N with

v˜(j) (rn ) q

−p qn

 |z| n  |z|  . h rn Γ( nq + 1)−p rn

where  q−1 X  p p   exp − x−q/p q q u xu q h(x) = u=0  1   1−x

(25)

if κ < 0 (i.e. p < 0, q > 0) if κ = 0.

The bound is generically tight up to subexponential factors. Proof.

In the case

Using the relation the

n),

κ = 0, the v˜n = ψn vn

choice of

rn

|z| < rn < α−1 for all n. bound v ˜k ≤ v˜(rn )/rnk (notice

ensures that

and the saddle-point

we obtain

(j) (j) (|z|) ≤ un; (z) ≤ vn; Now assume

p < 0:

∞ X k=n

q

q

Γ( nq + 1)−p rn

−p q (k−n)

k=n

Γ( kq + 1)

x −p

k−n

=

≤

q−1 ∞ X X

q

Γ( nq + 1)−p  |z| k−n Γ( kq + 1)−p rn

Γ( nq + 1)−p Γ( nq + 1 + t + uq )−p

Γ( nq + 1)−p

Γ( nq t=0 u=0 q−1 X

q

p qu

u=0 since

1 −p q (k−n)

then

Γ( nq + 1)−p

1

∞  |z| n X

v˜(j) (rn ) −p qn

xqt+u q −pt−pu/q

∞  X p −pt xqt ≤ h(x) x − q (−pt)! t=0 u

Γ( nq + 1)−p

≤

+1+t+

u −p q)

Γ( nq + 1 + t)−p

≤

1 (−p)−pt ≤ . t!−p (−pt)!

The estimates (22), (23) still hold, hence the tightness of the bound. Once again, the factor

n!p/q



in (24) can be made more explicit if necessary using

Lemma 11.

κ = T = 0 and K ∈ N, the vn; (z) actually admit (αz)n p(n), where p ∈ Q(αz)[n]. Indeed, starting PK K−1 from (18) and writing (for xed K ) (n + k + 1) = i=1 c[i] (n)(k + 1)i−1 , we get   ∞ K 1 (αz)n X (αz)n X (i − 1)! [i] K−1 k = (n+K+1) (αz) = c (n). (1 − αz)K n; (K − 1)! (K − 1)! i=1 (1 − αz)i In the important case where

closed-form expressions of the form

k=0

This is the kind of formula that appears in Example 1(c). Such bounds are easier to read than (24), but they are numerically unstable due to cancellations.

In a

system providing numerical routines for hypergeometric functions, one can use the alternative expression



1 (1 − αz)K



n

= (αz) n;



n+K −1 K −1

which does not suer from this shortcoming.



  1 n + K αz 2F1 n+1

18

MARC MEZZAROBBA AND BRUNO SALVY

Finally, note that it might be worthwhile looking for rened bounds in applications where

n

T

is large and

|z| ' α−1 ,

since (24) becomes tight only for very large

in this case. 5. Applications and Experiments

Implementation.

5.1.

We have implemented the algorithms described in this ar-

ticle (with slight variations) in the computer algebra system Maple v. 12.

Our 2

implementation is part of a development version of the Maple package gfun and will be included in an upcoming release, but the bound computation code is largely self-contained. It provides routines that compute majorant series for rational polynomials (following Ÿ3.2) and D-nite functions (Ÿ3.3, Ÿ4.1), and symbolic bounds for P-recursive sequences specied either using recurrence relations (Ÿ4.1) or as tails of D-nite series (Ÿ4.2). All examples of this article were computed using this implementation

3

.

4

It is also used by the Dynamic Dictionary of Mathematical Functions , an interactive web-based handbook of D-nite functions currently under development. All contents of the Dictionary are automatically generated from a compact description of each function (basically, a dierential equation and initial values) using a mix of symbolic computation algorithms and document templates. The webpages the system produces are interactive, in that they allow the user to trigger more computations, typically by asking for more terms in an asymptotic expansion. This is a situation where being able to display human-readable formulae rather than merely computing numerical bounds represents a signicant benet. Code based on this article provides majorant series for the Taylor expansions of the functions, truncation orders for these expansions to reach a given accuracy over a given disk, and symbolic bounds for their tails involving the truncation order.

Application to the Numerical Evaluation of D-Finite Functions.

5.2.

Guar-

anteed numerical computation with entire classes of functions usually involves the

automatic computation of error bounds relating approximations, e.g., by truncated power series, to the functions they approximate. Elementary results from real and complex analysis commonly used to compute such error bounds include the alternating series criterion, Cauchy's integral formula, and several variants of Taylor's theorem. Du and Yap (2005) provide bounds for the tails of the general hypergeometric series, where the parameters are allowed to vary, based on a detailed analysis of the variations of the coecient sequence. For the more general case of D-nite functions, another

ad hoc

method is given by van der Hoeven (1999).

A further classical tool is the Cauchy-Kovalevskaya majorant series method discussed in Ÿ3.1.

This idea is exploited by van der Hoeven (2001, Ÿ2.4) to bound

the tails of power series expansions of D-nite functions in the neighbourhood of an ordinary point of the equation, and later again in a much more general setting covering a wide range of functional equations (van der Hoeven, 2003). This is the approach we rely on in this article: indeed, the algorithm we described in Ÿ3.3 may actually be seen as a renement of those suggested in Ÿ3.5 and Ÿ5.2 of the latter 2

http://algo.inria.fr/libraries/papers/gfun.html

3The

development snapshot of gfun we used is available from http://algo.inria.fr/

mezzarobba/gfun-bounds-alpha.tgz. 4http://ddmf.msr-inria.inria.fr/

EFFECTIVE BOUNDS FOR P-RECURSIVE SEQUENCES

article.

19

The main originality of our approach is the asymptotic tightness of the

bounds. Finally, it should be noted that in the context of numerical evaluation, instead of using

a priori

bounds, it is often easier to compute successive error bounds in

parallel to successive approximations of the result, until the desired accuracy is reached. The computation of validated numerical enclosures of solutions of ODE, DAE and more general functional equations has been the subject of extensive literature since the sixties (see Rihm, 1994) in the area of interval methods. Of special interest when working with power series is the integration of dierential equations using Taylor models (see Hoefkens, 2001; Neher et al., 2007). Taylor models are one among a fair number of dierent symbolic-numeric representations of functions used in interval arithmetic, several of which have a similar approach of bounds for solutions of functional equations: for more on Taylor models and their relation to other interval methods, see (Makino and Berz, 2003; Neumaier, 2003).

Some of

these methods were imported to computer algebra and revisited by van der Hoeven (2007) in the context of rigorous eective complex analysis. In a nutshell, the common idea is to write the (dierential, say) equation at hand in xed-point form action of

Φ

u = Φ(u),

where

Φ

is an integral operator, and to consider the

on truncated power series augmented with error bounds, using rules

such as

Z

x

Z

2

(a0 + a1 t + a2 t + [α, β])dt ⊆

x

 x3  + [α, β] · B(x). (a0 + a1 t)dt + B a3 3

B(p) is an interval containing the range of p(x) obtained from the range of x. One then computes an approximate solution in the form of a Taylor expansion p(x) = a0 + · · · + an xn and iteratively searches for a tight interval [α, β] such that Φ(p + [α, β]) ⊂ p + [α, β], possibly narrowing the range of x or increasing the expansion order n as necessary. Under mild assumptions, the existence of such p + [α, β] implies that of an actual solution u ∈ p + [α, β] of Φ(u) = u. Here

While this is reported to provide tight numerical enclosures at reasonable cost for computations at machine precision even in the case of nonlinear equations in many variables, we are not aware of any asymptotic tightness result of the kind we are interested in in this paper.

In fact, it is not entirely clear to us under

which conditions methods of this kind are guaranteed to produce arbitrarily tight enclosures. (Note however that van der Hoeven (2007) states initial results in this direction.) Neither do we know how to use them to bound tails of D-nite functions on their whole disk of convergence. And yet, D-nite functions may be evaluated to an absolute precision softly linear time

n(log n)

O(1)

10−n

in

by computing truncations of their Taylor series by bi-

nary splitting. Numerical analytic continuation based on this technique then allows to obtain values of these functions at any point of their Riemann surfaces (Chudnovsky and Chudnovsky, 1987, Ÿ4). Applications include the numerical computation of monodromy matrices of linear dierential equations with polynomial coefcients.

In this context, one benet of the language of majorant series is that a

single majorant encodes both bounds on the values and truncation orders for all elements of a basis of the local solutions of the dierential equations as well as their derivativesall of which are useful to control errors in the numerical analytic continuation process.

20

MARC MEZZAROBBA AND BRUNO SALVY

Excluding degenerated cases, the number of terms of the series to take into

λn+o(n), where λ depends on the position of the evaluation point relative O(n/ log n) in the case of entire functions. tightness result of Theorem 1 translates into the fact that the number N of

account is

to the singularities of the function, or The

terms that get computed is indeed of that order, while most existing methods for computing bounds of seem to ensure only

N = O(n).

This in turn improves the

complexity of the algorithm by a constant factor. The development version of gfun mentioned above contains high-precision numerical evaluation and analytic continuation routines based on this strategy, that rely on the bound computation code for accuracy control. We plan to add similar numerical evaluation facilities to the DDMF in a near future. 5.3.

Experiments.

In Table 1, we report on experiments concerning the tightness

of the bounds for truncating Taylor series expansions of a few common elementary and special functions. Each column label actually stands for a dierential equation that cancels the given function (with suitable initial values), and an evaluation point smaller in absolute value than the dominant singularity of the dierential equation. Each internal cell shows the truncation order computed by the algorithm from this data for a specic accuracy requirement, and compares it to the minimal

2

correct answer, computed by exhaustive search. For instance, the column  erf(1)  corresponds to the evaluation at

z=1

of the function

u(z) = erf(z)

2

represented

as the unique solution of

(2 + 8z 2 ) u0 (z) + 6z u00 (z) + u000 (z),

u(0) = 0,

8 . π |u;189 (1) − u(1)| ≤

u0 (0) = 0,

u00 (0) =

Using a majorant series for u, our algorithm determined that 10−100 , but it happens that only the rst 163 of these 189 terms are really necessary.

It can be seen that the bounds we compute do not stray too far from the optimal values. We consider three cases, corresponding to the three main types of asymptotic behaviours that the coecient sequence of a convergent D-nite series may exhibit, characterized (in generic cases) by the nature of the dominant singularities of the dierential equation: regular singularities (κ

=0=T =0

with the notations of

the previous sections), irregular singularities at nite distance (κ innity (κ

< 0).

(Irregular singularities with

κ>0

= 0, T > 0),

or at

correspond to divergent power

series, and a dierential equation whose only singularity is a regular singular point at innity has only polynomial solutions. The examples of the second set all involve right composition by rational functions because it is unusual to study dierential equations with more than two irregular singular points, and those are usually taken to be

∞

and

0.)

For each of these, the last three columns illustrate how the truncation orders and the bounds vary as

|z|

approaches the radius of convergence of the series. Fortu-

nately, high-order Taylor expansions at values of D-nite functions for such

z:

0 are not the best way to compute numerical it is much more ecient to use several steps

of analytic continuation along a broken-line path from

0

to

z

(Chudnovsky and

Chudnovsky, 1987, Ÿ4). The example of

Si(z)

has an interesting feature: the origin is a regular singular

point of the dierential equation mentioned in Example 1(d), but

Si(z)

may never-

theless be dened by simple initial values at origin, so that our algorithm applies without any adjustment.

EFFECTIVE BOUNDS FOR P-RECURSIVE SEQUENCES

21

Regular dominant singularity

1 (1−z)2

10−10 10−100 10−1000

@ 21

cos z 1−z

cos z 1−z 2

@ 21

cos z (1−z)2

(z+1)2 cos z (z 3 +z+1)2

@ 12

1 @ 10

46/40

46/34

54/33

54/39

24/12

350/342

350/333

364/331

364/341

140/121

3346/3323

3366/3321

ψ(1/2)

arctan 21

3346/3335

arccot(z) (z 2 −1)(z 2 +5)

10−10 10−100 10−1000

@ 12

@ 21

3366/3334

9 arctan 10

1232/1201

99 arctan 100

68/27

40/23

44/28

336/164

4238/1496

386/321

342/313

348/324

2338/2108

25210/21848

3400/3307

3336/3293

3344/3310

22050/21754

231844/227810

z 1 exp (1−z) 2 @2

z 1 exp 1−z 2 @2

erf

Finite irregular dominant singularity

z cos 1−z @ 31

10−10 10−100 10−1000

z sin 1−z @ 13

54/25

52/24

306/224

306/225

2469/2150

exp(1/(1−z)) (1−z)

10−10 10−100 10−1000

118/79 557/497

2469/2149

@ 21

79/54

Bi



1 1−z



4153/4001

@ 21

Ai



1 1−z



@ 12



1+z 2z 2 −1



67/42

28/12

416/364

244/132

3566/3432

Ai



1 1−z



@ 34

@ 19

2384/1292

Ai



1 1−z



@ 78

183/56

175/30

2202/77

35381/215

438/387

729/416

725/345

4199/879

40969/2025

3625/3490

4904/3645

4900/3406

17859/8372

Si(1)

cos(1)

84162/18529

Dominant singularity at innity

Ai(4i + 4) −10

10 10−100 10−1000 −10

10 10−100 10−1000

Bi(4i + 4)

sin(1)

102/59

102/59

18/12

20/13

231/200

231/200

73/68

75/69

75/70

1058/1031

1058/1031

453/448

455/449

455/450

e−100

erf 2 (1)

erf(1)

384/291

60/33

38/24

522/450

189/163

1446/1402

1036/1011

erf(10)

18/14

erf(100)

790/574

71924/54388

150/138

1062/894

72248/54800

910/898

2902/2800

75422/58772

Table 1. Computed/minimal required number of terms of the

Taylor expansion of a D-nite function to approximate this func-

ψ is the solution (1 − z 2 ) ψ 00 (z) − 2(b − 1)z ψ 0 (z) +

tion to a given absolute precision. In this table, of the spheroidal wave equation

(c − 4qz 2 )ψ(z) = 0 given by the choice of parameters and initial 0 values b = 1/2, q = 1/3, c = 1, ψ(0) = 1, ψ (0) = 0; Ai and Bi denote the Airy functions; erf stands for the error function and Si for the integral sine.

Finally, here is a nontrivial non-generic example where our method fails to produce a tight bound.

Example. sequences

ζ(3), Apéry (1979) introduces two un = bn − ζ(3)an satises the (minimal-order)

In his proof or the irrationality of

(an )

and

(bn )

such that

22

MARC MEZZAROBBA AND BRUNO SALVY

linear recurrence relation

(n+2)3 un+2 = (2n+3)(17n2 +51n+39) un+1 −(n+1)3 un ,

u0 = −ζ(3), u1 = 6−5ζ(3).

Applied to this recurrence relation, Algorithm 5 determines that

|un | ≤ 1.21 (n2 + 3n + 2) (17 −

√

12)−n

This bound is asymptotically tight for both Apéry's proof is that

Acknowledgements.

bn − ζ(3)an → 0

fast as

(where

1/(17 −

an and bn , n → ∞.

√

12) ' 33.97)

but the whole point of

The authors would like to thank Moulay Barkatou and Nicolas

Le Roux for interesting discussions and useful suggestions. References

Abramowitz, M., Stegun, I. A., 1972. Handbook of Mathematical Functions. Dover Publications, New York. Apéry, R., 1979. Irrationalité de

ζ(2)

et

ζ(3).

Astérisque 61, 1113.

Bronstein, M., 2005. Symbolic integration. I, 2nd Edition. Vol. 1 of Algorithms and Computation in Mathematics. Springer-Verlag, Berlin. Buslaev, V. I., Buslaeva, S. F., 2005. Poincare theorem for dierence equations. Mathematical Notes 78 (5-6), 877882, translated from Matematicheskie Zametki, vol. 78, no. 6, 2005, pp. 943947. Chudnovsky, D. V., Chudnovsky, G. V., 1987. Computer assisted number theory with applications. In: Number theory (New York, 19841985). Vol. 1240 of Lecture Notes in Mathematics. Springer, Berlin, pp. 168. Chudnovsky, D. V., Chudnovsky, G. V., 1988. Approximations and complex multiplication according to Ramanujan. In: Ramanujan revisited. Academic Press, Boston, MA, pp. 375472. Du, Z., Yap, C., 2005. Uniform complexity approximating hypergeometric functions with absolute error. In: Pae, S.-I., Park, H. (Eds.), Proceedings of the 7th Asian Symposium on Computer Mathematics (ASCM 2005). Korea Institute for Advanced Study, pp. 246249. Flajolet, P., Salvy, B., Zimmermann, P., Feb. 1991. Automatic average-case analysis of algorithms. Theoretical Computer Science, Series A 79 (1), 37109. Flajolet, P., Sedgewick, R., 2009. Analytic Combinatorics. Cambridge University Press. URL

http://algo.inria.fr/flajolet/Publications/book.pdf

Gourdon, X., Salvy, B., 1996. Eective asymptotics of linear recurrences with rational coecients. Discrete Mathematics 153 (1-3), 145163. Graham, R. L., Knuth, D. E., Patashnik, O., 1989. Concrete Mathematics. Addison Wesley. Guelfond, A. O., 1963. Calcul des diérences nies. Collection Universitaire de Mathématiques, XII. Dunod, Paris, translated from the Russian. Hoefkens, J., 2001. Rigorous numerical analysis with high-order Taylor models. Ph.D. thesis, Michigan State University. Ince, E. L., 1956. Ordinary Dierential Equations. Dover, New York. Knuth, D. E., 1997. Seminumerical Algorithms, 3rd Edition. Vol. 2 of The Art of Computer Programming. Addison-Wesley, Reading, Massachusetts. Kooman, R., Tijdeman, R., 1990. Convergence properties of linear recurrence sequences. Nieuw Archief voor Wiskunde, Ser. 4 4, 1325.

EFFECTIVE BOUNDS FOR P-RECURSIVE SEQUENCES

23

Kreuser, P., 1914. Über das Verhalten der Integrale homogener linearer Dierenzengleichungen im Unendlichen. Ph.D. thesis, Universität Tübingen, Borna-Leipzig. Lipshitz, L., 1989.

D-nite

power series. Journal of Algebra 122 (2), 353373.

Makino, K., Berz, M., 2003. Taylor models and other validated functional inclusion methods. International Journal of Pure and Applied Mathematics 4 (4), 379456. URL

http://bt.pa.msu.edu/pub/papers/TMIJPAM03/TMIJPAM03.pdf

Malgrange,

B.,

1974.

Sur

les

points

singuliers

des

équations

diérentielles.

L'Enseignement mathématique XX (1-2), 147176. Meschkowski, H., 1959. Dierenzengleichungen. Vandenhoeck & Ruprecht. Mezzino, M., Pinsky, M., Dec. 1998. Leibniz's formula, Cauchy majorants, and linear dierential equations. Mathematics Magazine 71 (5), 360368. Neher, M., Jackson, K. R., Nedialkov, N. S., 2007. On Taylor model based integration of ODEs. SIAM Journal on Numerical Analysis 45 (1), 236262. URL

http://link.aip.org/link/?SNA/45/236/1

Neumaier, A., 2003. Taylor forms  Use and limits. Reliable Computing 9 (1), 4379. Odlyzko, A. M., 1995. Asymptotic enumeration methods. In:

Graham, R. L.,

Groetschel, M., Lovasz, L. (Eds.), Handbook of Combinatorics. Vol. 2. Elsevier, pp. 10631229. URL

http://www.dtc.umn.edu/~odlyzko/doc/asymptotic.enum.pdf

Perron, O., 1909a. Über die Poincarésche lineare Dierenzengleichung. Journal für die reine und angewandte Mathematik 137, 664. Perron, O., 1909b. Über ein Satz des Herrn Poincaré. Journal für die reine und angewandte Mathematik 136, 1737. Perron, O., 1910. Über lineare Dierenzengleichungen. Acta Mathematica 34, 109 137. Perron, O., 1921. Über Summengleichungen und Poincarésche Dierenzengleichungen. Mathematische Annalen 84, 115. Pituk, M., 1997. Asymptotic behavior of a Poincaré recurrence system. Journal of Approximation Theory 91 (2), 226243. Poincaré, H., 1885. Sur les équations linéaires aux diérentielles ordinaires et aux diérences nies. American Journal of Mathematics 7 (3), 203258. Rihm,

R.,

1994.

Interval

methods

for

initial

value

problems

Herzberger, J. (Ed.), Topics in Validated Computations:

in

ODEs.

In:

Proceedings of the

IMACS-GAMM International Workshop on Validated Computations, University of Oldenburg. Elsevier Studies in Computational Mathematics. Elsevier, pp. 173 207. Salvy, B., Zimmermann, P., 1994. Gfun: A Maple package for the manipulation of generating and holonomic functions in one variable. ACM Transactions on Mathematical Software 20 (2), 163177. Schäfke, F. W., Feb. 1965. Lösungstypen von Dierenzengleichungen und Summengleichungen in normierten abelschen Gruppen. Mathematische Zeitschrift 88 (1), 61104. Schönhage, A., 1982. The fundamental theorem of algebra in terms of computational complexity. Tech. rep., Mathematisches Institut der Universität Tübingen. Stanley, R. P., 1980. Dierentiably nite power series. European Journal of Combinatorics 1 (2), 175188.

24

MARC MEZZAROBBA AND BRUNO SALVY

Stanley, R. P., 1999. Enumerative combinatorics. Vol. 2. Cambridge University Press. van der Hoeven, J., 1999. Fast evaluation of holonomic functions. Theoretical Computer Science 210 (1), 199216. URL

http://www.math.u-psud.fr/~vdhoeven/Publs/1997/TCS.ps.gz

van der Hoeven, J., 2001. Fast evaluation of holonomic functions near and in regular singularities. Journal of Symbolic Computation 31 (6), 717743. URL

http://www.math.u-psud.fr/~vdhoeven/Publs/2000/singhol.ps.gz

van der Hoeven, J., 2003. Majorants for formal power series. Tech. Rep. 2003-15, Université Paris-Sud, Orsay, France. URL

http://www.math.u-psud.fr/~vdhoeven/Publs/2003/maj.ps.gz

van der Hoeven, J., 2007. On eective analytic continuation. Mathematics in Computer Science 1 (1), 111175. Wimp, J., Zeilberger, D., 1985. Resurrecting the asymptotics of linear recurrences. Journal of Mathematical Analysis and Applications 111, 162176. Zeilberger, D., 1990. A holonomic systems approach to special functions identities. Journal of Computational and Applied Mathematics 32 (3), 321368. Zeilberger, D., Apr. 2008. Asyrec: A Maple package for computing the asymptotics of solutions of linear recurrence equations with polynomial coecients. URL

html

http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/asy.

: [email protected] : http://algo.inria.fr/mezzarobba

E-mail address URL

: [email protected] : http://algo.inria.fr/salvy

E-mail address URL

Algorithms Project, Inria Paris-Rocquencourt, France