Formulas for Continued Fractions. An Automated Guess and Prove Approach
aa r X i v : . [ c s . S C ] J u l FORMULAS FOR CONTINUED FRACTIONSAN AUTOMATED GUESS AND PROVE APPROACH
S´EBASTIEN MAULAT AND BRUNO SALVY
Abstract.
We describe a simple method that produces automatically closedforms for the coefficients of continued fractions expansions of a large num-ber of special functions. The function is specified by a non-linear differentialequation and initial conditions. This is used to generate the first few coef-ficients and from there a conjectured formula. This formula is then provedautomatically thanks to a linear recurrence satisfied by some remainder terms.Extensive experiments show that this simple approach and its straightforwardgeneralization to difference and q -difference equations capture a large part ofthe formulas in the literature on continued fractions. Introduction
Continued fractions are well known for their approximation properties, their usein acceleration of convergence and analytic continuation, as well as their applica-tion in proofs of irrationality. Any formal power series can be converted into a corresponding continued fraction (C-fraction)(1) a + a ( z )1 + a ( z )1 + a ( z )1 + · · · classically denoted a + K ∞ m =1 a m ( z )1 or [ a , a ( z ) , a ( z ) , · · · ], where a is a constantand a i ( z ) are nonconstant monomials for i > partial numerators .In the frequent case when all the exponents are equal to 1, the C-fraction is called regular . Truncating a continued fraction after its n th term gives a rational functionwhich is called its n th convergent . There is a one-to-one correspondence betweenpower series and C-fractions. It is easily computed by a sequence of extractions ofthe constant coefficient, division by the variable and inversion. This conversion isavailable in the major computer algebra systems.In several isolated cases, the coefficients a m are known to possess a closed form,as in the following formula for exp( z ):exp( z ) = 1 + z − z/ (2 · z/ (2 · − z/ (2 · z/ (2 · · · · , or more compactly(2) a = z, a k = − z/ (2(2 k − , a k +1 = z/ (2(2 k + 1)) . Such formulas are the object of this work. A number of them are listed in theclassical handbook by Abramowitz and Stegun [1], or in its successor [15] and themost extensive list to date is the recent handbook by Cuyt et alii [10]. Our aim isto derive many of these formulas automatically, starting from a description of thefunction to be expanded in continued fraction.We concentrate on functions that are given as solutions of ordinary differentialequations with initial conditions (or difference or q -difference equations, see § § §
6) and proves more suited to the targetted application to the
DynamicDictionary of Mathematical Functions [5]. This is an online encyclopedia of specialfunctions, where the formulas are all generated by computer algebra algorithmsfrom differential equations, in many cases along with a human-readable proof. Inthis context, it makes sense to avoid any table lookup and generate formulas andproofs for continued fractions directly from the differential equation.The work closer to ours is the investigation by Chudnovsky and Chudnovsky [8].They used computer algebra in the study of formulas for continued fractions. Theiraim was to classify all functions possessing continued fractions with explicit formulasof various types and relating them to Painlev´e transcendents. In contrast, we focuson one function that is given as input, and heuristically produce a rational continuedfraction expansion when possible.This article is structured as follows. Section 2 gives an overview of our methodon the example of the tangent function. Next, Section 3 presents a heuristic ofindependent interest that reduces the order of a recurrence given initial conditions.This plays a crucial role in the proving phase of our method. Section 4 is a briefaccount of what guessing means in this context, while Section 5 is the heart of this
ORMULAS FOR CONTINUED FRACTIONSAN AUTOMATED GUESS AND PROVE APPROACH3 work and shows how proofs are achieved automatically. Finally, Section 6 presentsexperiments with this approach.2.
Detailed Example: tan
The tangent function can be defined by the Riccati equation(3) y ′ = 1 + y , y (0) = 0 . The first 15 coefficients of the unique power series solution are easily computed fromthe differential equation (see Proposition 3 below for existence and uniqueness). Aconversion into a continued fraction gives the coefficients[0 , z, − z / , − z / , − z / , − z / , − z / , − z / . The general formula can be deduced from these first terms by rational interpolation,which leads automatically to the (so far conjectural) formula(4) a ( z ) := z ; a n ( z ) := − z / ((2 n − n − , n > . Next, we turn to the automatic proof of this formula. The strategy is to provethat the sequence of rational functions defined by truncating (4) after the n th termfor n = 1 , , . . . converges to the formal power series solution to the differentialequation (3). More precisely, let f n be defined by f n = P n Q n := [0 , z, − z / , . . . , − z / ((2 n − n − , where the rightmost term denotes the finite continued fraction. Then the proof willbe completed by showing thatval(tan − f n ) → ∞ as n → ∞ , where val denotes the valuation of a power series:val (cid:0) X i ≥ c i z i (cid:1) := min { i ≥ | c i = 0 } , with the convention val(0) = ∞ . Proposition 3 below shows that it is sufficient toprove that val( D ( f n )) → ∞ , where D ( f n ) := f ′ n − − f n .It is classical that the numerator and denominator of the convergents of a con-tinued fraction are related to the coefficients a n through a linear recurrence:(5) ( P − , P ) = (1 , , P n = P n − + a n P n − , n ≥ , ( Q − , Q ) = (0 , , Q n = Q n − + a n Q n − , n ≥ . In view of (4), it follows that for all n ≥ Q n (0) = 1. Thus, the valuationof D ( f n ) is that of its numerator(6) H n := P ′ n Q n − Q n − P n − P n Q ′ n . Using (5) to rewrite P n + k and Q n + k in terms of P n , P n +1 , Q n , Q n +1 , it follows thatany shift H n + k ( k ∈ { , , , . . . } ) can be rewritten as a linear combination of P ′ n + i Q n + j , P n + i Q ′ n + j , P n + i P n + j , Q n + i Q n + j , for i and j in { , } . There are finitely many such terms, which implies that alinear dependency between H n , H n +1 , . . . (ie, a linear recurrence for H n ) can be S´EBASTIEN MAULAT AND BRUNO SALVY computed directly from (5) by linear algebra. This computation produces a linearrecurrence of order 4:(7) (2 n + 7) z H n − z (2 n + 7)(2 n + 3) (2 n + 1) H n +1 + 2 z (2 n + 5)(2 n + 3) (2 n + 1) (4 n − z + 20 n + 21) H n +2 − (2 n + 5) (2 n + 1) (2 n + 7) (2 n + 3) H n +3 + (2 n + 5) (2 n + 1) (2 n + 7) (2 n + 3) H n +4 = 0 . This recurrence is satisfied by all sequences defined by (6), with P n and Q n arbitrarysolutions of (5). Using the actual sequences P n and Q n provided by the continuedfraction gives the first values of H n : − , − z , − z , − z , − z , − z . From there, automatic guessing again suggests the following simpler recurrencefor H n :(8) (2 n + 1) H n +1 − z H n = 0 . And again, this recurrence admits of an automatic proof: the right Euclidean divi-sion of the fourth order recurrence operator from (7) by this first order one has aremainder equal to 0. This shows that the solution of (8) with the initial conditionsgiven above coincides with the solution of (7) with the same initial conditions, andthus the numerator of D ( f n ) satisfies (8). On this last recurrence, the increase ofthe valuation with n is clear and this concludes the proof that f n converges to tanand thus that the power series defined by the continued fraction (4) is that of tan.In summary, starting with the differential equation (3), this method producesand proves automatically the general term of the famous continued fractiontan z = z − z / − . . .1 − z / ((2 n − n − π is irrational.3. Reduction of order by Guess and Prove
The transformation of the large recurrence (7) into the shorter one (8) makesit possible to prove automatically that the valuations val H n ( z ) increase with n .This transformation turns out to play a role in most of the examples dealt within our experiments. It is actually of more general interest: the closure propertiesenjoyed by the class of D-finite series or P-recursive sequences give rise to operatorssatisfied by products or sums of zeroes of such operators [17, 16]. These operatorsannihilate all possible cases and are potentially of large size, while operators ofsmaller order may exist for the specific solution of interest. Such an operator maybe a right factor of the large one and could be searched for by factoring, but thisis made difficult by the potentially infinite number of distinct factorizations [18]. ORMULAS FOR CONTINUED FRACTIONSAN AUTOMATED GUESS AND PROVE APPROACH5
Sequences.
Let A be a recurrence operator with polynomial coefficients in n ,of order denoted by ord A and leading coefficient lc( A )( n ). A sequence ( u n ) n ≥ ,abbreviated ( u n ), is said to be defined by the operator A and the initial condi-tions K = ( u i ) i ∈I , when the value u n is given by K for n ∈ I and by the recurrenceoperator evaluated at n − ord A otherwise. Note that the set I must contain { , . . . , ord A − } ∪ { i ∈ N | lc( A )( i − ord A ) = 0 } . Algorithm.
We now detail an efficient heuristic approach finding such rightfactors, whose complexity is controlled with respect to the size of the large operator.The idea is to exploit the initial values of the sequence by a “guess and prove”approach. This is described in Algorithm 1. This algorithm takes as input a linearrecurrence operator A and initial values, as well as an upper bound N on thenumber of coefficients used to find a right factor. It is described here in the case ofrecurrence operators; similar variants apply to differential or q -difference cases.The search for a smaller order operator is performed in two main steps, “guess-ing” and “proving”. First, the input recurrence of order M and its initial values areused to compute the first N terms of the sequence. Next, these N terms are usedto “guess” a linear recurrence. This is done by linear algebra: we search succes-sively for the existence of a linear recurrence operator G of order 1 , , . . . , M withpolynomial coefficients of degrees such that the sum of the numbers of undeter-mined coefficients of the recurrence is smaller than N . The structure of this linearalgebra problem is exploited by computing matrix rational interpolants [4] (in thedifferential case, Hermite-Pad´e approximants are used [3]).When N is sufficiently large, this linear algebra phase is always successful, sinceit can reconstruct A . The next step is to prove that the recurrence G obtained fromthe first N terms of ( u n ) defines the same sequence for all n . The operator G isnot necessarily a right factor of A , but could be merely a left multiple of such aright factor, the factor itself being too large to be found with N terms only. Thisis related to the typical shape of the order-degree curve [7]. Thus the algorithmnext computes the greatest common right divisor of G and A and its numerator R ,obtained by left-multiplication with the least common multiple of the denominatorsof the coefficients.At this stage, the algorithm has produced a right factor R of A . It is thenassociated initial conditions ( u i ) i ∈J , with which R defines a sequence ( v n ) n ≥ . Wenow prove that if v n = u n for n ∈ J + { , . . . , ord A − ord R} , then v n = u n for all n . The induction on n is based on the following. Lemma 1. If u n = v n for n ≤ i +ord A − and i + j / ∈ J for all j ∈ { ord R , . . . , ord A} then u i +ord A = v i +ord A .Proof. The sequences { S j R · v n } ord R≤ j ≤ ord A all cancel at n = i . The applicationof A at ( v n ) n ≥ is a linear combination of them with coefficients that are finite, asshown in the next lemma, so that A · v n is 0 at n = i . (cid:3) Lemma 2.
Let A and C be recurrence operators with polynomial coefficients, sat-isfying A = BC where B has rational coefficients. Then the denominator den B satisfies: den( B ) − (0) ⊆ lc( C ) − (0) + { , − , . . . , − ord B} where addition denotes the sumset. S´EBASTIEN MAULAT AND BRUNO SALVY
Algorithm 1
Reduction of Order
Input: ( A , ( u n ) n ∈I ) defining ( u n ) n ≥ , and N > Output: ( R , ( u n ) n ∈J ) defining ( u n ) n ≥ s.t. ord R ≤ ord A . U ← ( u n ) n =0 ,...,N − , computed using A and ( u n ) n ∈I G ← guessrec( U ) if G 6 = FAIL then
R ← numer (gcd right ( A , G )) J ← I ∪ (cid:16) lc( R ) − (0) − ord R (cid:17) ∩ N V ← ( v n ) n ∈J + { ,..., ord A− ord R} , using ( R , ( u n ) n ∈J ) U ′ ← ( u n ) n ∈J + { ,..., ord A− ord R} , using ( A , ( u n ) n ∈I ) if U ′ = V then return ( R , ( u n ) n ∈J ) end ifend ifreturn ( A , ( u n ) n ∈I ) Proof.
This is seen by following the steps of a right Euclidean division. (cid:3)
In practice, this algorithm is run for increasing values of N = 4 , , , . . . andstopped when either a factor is found or N is larger than the number of coefficientsof A . Note however that if A is not irreducible, then increasing N further is boundto find a nontrivial right factor.4. Guessing Continued Fractions
The first step of our approach to continued fractions is the automatic discoveryof formulas for the partial numerators a k . This section is very short since this partof the computation is straightforward.Starting from the differential equation, a first method would be to produce thefirst terms of the series expansion of the function, convert them into the first termsof the continued fraction and then use the method of the previous section to lookfor a linear recurrence of size bounded by the number of terms that have been com-puted. It turns out that in most of the known examples, explicit formulas are ofrational form (see Section 6). We therefore concentrate on rational coefficients, oron interlacing of rational coefficients as in the example of the exponential function.This means that the “guessing” stage of our approach relies simply on rational in-terpolation, problem for which efficient algorithms are known through its relationto the extended Euclidean algorithm [11, § Proving Continued Fractionsfor solutions of ordinarydifferential equations
The proving phase is the heart of our work. It is also where most of the compu-tational work takes place. We consider first-order non-linear differential equationswith rational coefficients, ie, y ′ = p ( y ), with p ∈ Q ( z )[ Y ] of degree d . In particular,the case d ≤ ORMULAS FOR CONTINUED FRACTIONSAN AUTOMATED GUESS AND PROVE APPROACH7 classes of equations have been provided by Euler and Lagrange and more recentlyby Khovanskii [13].Our procedure goes in the reverse direction. The continued fraction with explicitrational coefficients that was found in the previous stage defines a power series. Theaim is to show that it is a solution of the differential equation.5.1.
Valuations.
The following proposition reduces the proof to that of the ulti-mate increase of a sequence of integers.
Proposition 3.
Let F ∈ K [[ X, Y ]] be a formal power series with coefficients in afield K and let ( f n ( X )) be a sequence of power series in K [[ X ]] . Then the differentialequation Y ′ = F ( X, Y ) with initial condition Y (0) = 0 admits a unique power seriessolution S ( X ) . Moreover, the sequence ( f n ( X )) converges to S ( X ) (ie val( f n − S ) →∞ ) if and only if f n (0) = 0 for n sufficiently large and val( f ′ n ( X ) − F ( X, f n ( X ))) →∞ . Note that equations with an initial condition Y (0) = a = 0 can often be broughtto this setting by changing the unknown function into a + Y . Proof.
Recall that the algebra of power series is a metric space for the distanceinduced by the valuation: if f and g are two power series, then d ( f, g ) = 2 − val( f − g ) ,where val denotes the valuation (this distance does not derive from a norm). It isa simple consequence of the definition that Cauchy sequences for this distanceconverge in K [[ X ]].The first part of the proposition is a variant of Cauchy’s theorem, whose proofis straightforward thanks to Taylor expansions. In detail, the solutions of Y ′ = F ( X, Y ) with initial condition Y (0) = 0 are the fixed points of the operator G : Y R F ( X, Y ); this operator is a contraction:val( G ( Y ) − G ( Y )) = val (cid:18)Z F ( X, Y ) − F ( X, Y ) (cid:19) = val (cid:18)Z ∂F∂Y ( X, Y )( Y − Y ) + O (( Y − Y ) ) (cid:19) > val( Y − Y );this shows both the existence of a solution (start from Y = 0, iterate G and usecompleteness) and its uniqueness.Next, if val( f ′ n − F ( X, f n )) = K , while f n (0) = 0 = S (0), then S − f n = S (0) − f n (0) + Z ( F ( X, S ) − F ( X, f n )) + O ( x K +1 )= ( G ( S ) − G ( f n )) + O ( x K +1 ) . The previous inequality with Y = S and Y = f n shows that the valuation ofthe first term on the right-hand side is larger than that of the left-hand size andthus val( S − f n ) ≥ K + 1, which shows that f n → S . The converse implicationfollows from the continuity of the map Y Y ′ − F ( X, Y ). (cid:3) This proposition extends to more general equations of the type P ( z, y, y ′ , . . . , y ( n ) ) =0, with natural assumptions on the initial and separant of the equation. S´EBASTIEN MAULAT AND BRUNO SALVY
P-recursivity and Convergents.
Recall that a sequence is called P-recursivewhen it satisfies a linear recurrence with coefficients that are polynomial in theindex. P-recursive sequences are closed under sum and product and algorithmscomputing the corresponding recurrences are known and implemented [17, 16].The key to our approach is the following.
Theorem 4.
Let ( P k ( z )) and ( Q k ( z )) be P-recursive sequences of rational functionsin z and let F ∈ K ( z )[ Y ] be a polynomial of degree d > in Y . Then the sequence H k := Q max(2 ,d ) k (cid:18) P k Q k (cid:19) ′ − F (cid:18) z, P k Q k (cid:19)! satisfies a linear recurrence with coefficients in K [ z, k ] . This theorem is used when ( P k ) and ( Q k ) are the sequences of numeratorsand denominators of the continued fraction supposed to converge to a solutionof y ′ = F ( z, y ). Its proof constructs a recurrence for H k from which the increaseof valuation will be obtained using the reduction of order of Section 3. This issufficient thanks to Proposition 3 and the observation that val Q k = 0, which willfollow by induction from Eq. (5) and the fact that val a k > Proof.
Let M be the order of the recurrence satisfied by ( P k ). Using this recurrence,all P k + i , i ∈ N can be rewritten as linear combinations of P k + j for j = 0 , . . . , M − Q ( k, z ), while the polynomials P ′ k + i rewrite as linear combina-tions of those same polynomials complemented by P ′ k + j for j = 0 , . . . , M −
1. Thesame argument applies to the sequence ( Q k ) and we denote by M ′ the order of therecurrence it satisfies.The choice of the exponent of Q k makes H k a polynomial of degree d in P k , Q k , P ′ k and Q ′ k . Thus all the H k + i for i ∈ N can be rewritten as linear combinations ofmonomials of degree d in P k + i , Q k + j , i = 0 , . . . , M − j = 0 , . . . , M ′ − N , whence a linear dependencybetween H k , . . . , H k + N (ie, a linear recurrence of order at most N satisfied by ( H k )).It can be found by linear algebra. (cid:3) As exemplified by the computation in the example of tan in Section 2, the boundon the order on specific examples may not be as large as suggested by this proof.Our implementation thus proceeds by increasing the order one by one and lookingfor a linear dependency until one is found. The algorithm is outlined in Algorithm 2.Its termination is granted by the theorem.We state a simple generalization of this result that could be useful in applicationsto continued fractions: if the partial numerators in the continued fraction expan-sion (1) are of the form r ( k ) z p ( k ) with r rational and p polynomial, then again, thepolynomials H k defined in the theorem satisfy a linear recurrence, this time withcoefficients that are polynomials in k, z and a finite number of z k i , with i ≤ deg p .The proof follows along the same lines.5.3. Riccati Equations.
The case when the polynomial F of Theorem 4 has de-gree 2 gives rise to Riccati equations that are ubiquitous in the theory of continued ORMULAS FOR CONTINUED FRACTIONSAN AUTOMATED GUESS AND PROVE APPROACH9
Algorithm 2
Recurrence for ( H k ) k ≥ Input: linear recurrences L P and L Q of order bounded by M for ( P k ) k ≥ and ( Q k ) k ≥ Output: a linear recurrence L H for ( H k ) k ≥ T ( k ) ← H k for i = 1 , , , . . . do T i ( k ) ← T i − ( k + 1) with P k + M , P ′ k + M , Q k + M , Q ′ k + M rewritten in termsof values of these sequences with smaller indices, using L P , L Q and theirderivatives. if the linear equation P i − j =0 c j T j ( k ) + T i ( k ) in the unknowns c , . . . , c i − hasa solution thenreturn H k + i + c i − H k + i − + · · · + c H k = 0 end ifend for fractions [2, 10.7]. In this case, the computation of a recurrence of the form pre-dicted by the theorem can be made explicit in full generality. Proposition 5.
Let K ∞ k =1 a k ( z )1 be a solution of the Riccati differential equation Y ′ = F ( z, Y ) where F is a polynomial in K ( z )[ Y ] of degree 2 in Y , let ( P k ) and ( Q k ) besequences obeying the linear recurrence u k +2 = u k +1 + a k +2 ( z ) u k , with a ′ k ( z ) = 0 .Let finally H k be defined by H k = Q k (( P k /Q k ) ′ − F ( z, P k /Q k )) . Then the sequence ( H k ) satisfies the following linear recurrence of order 4: a ′ k +1 H k +1 + (cid:18) a k a ′ k − a k +1 + 1 a ′ k +1 (cid:19) H k − (cid:18) a k ( a k + 1) a ′ k + a k +1 ( a k +1 + 1) a ′ k +1 (cid:19) H k − − (cid:18) a k + 1 a ′ k − a k +1 a ′ k +1 (cid:19) a k H k − + a k − a k a ′ k H k − = 0 . The shift of the indices (from H k +1 to H k − ) is only for readability. A niceproperty of this recurrence is that its coefficients do not depend on the differentialequation, but only on the sequence a k . This persists for higher degrees: a differentialequation with a cubic right-hand side leads to a recurrence of order 6 that does notdepend on the equation. Proof.
The formula is obtained automatically by the method from the proof ofTheorem 4, on a differential equation with symbolic coefficients. It could also bederived by hand. However, once it is given, it is a simple matter to produce a proof:inject the definition of H k into the recurrence, rewrite all the P k ’s and Q k ’s usingthe recurrence they satisfy in terms of P k − , Q k − , P k − , Q k − and collect terms toobserve that the left-hand side becomes 0. (cid:3) As an example, setting a k ( z ) = − z / ((2 k − k − Corollary 6.
If the sequences ( P k ) and ( Q k ) satisfy a linear recurrence of theform u k +2 = b k +2 ( z ) u k +1 + a k +2 ( z ) u k , then the sequence ( H k ) satisfies a fourth-order linear recurrence obtained by evaluating that of Prop. 5, replacing a by a /b and a k by a k / ( b k b k − ) for k ≥ .Proof. This is a classical transformation of continued fractions. Setting ˜ P k = P k / ( b · · · b k − b k ) and similarly for ˜ Q k and injecting into the recurrence equationshows that both sequences ( ˜ P k ) and ( ˜ Q k ) satisfy u k +2 = u k +1 + a k +2 b k +2 b k +1 u k . Since ˜ P k / ˜ Q k = P k /Q k , the proposition applies. (cid:3) Nonregularity and Periodicities.
As the example of the continued fractionfor the exponential function in Eq. (2) shows, not all common closed forms forcontinued fractions are given by one rational function. However, most C-fractionsformulas in the literature appear to be “periodic”, in the sense that there exists aperiod ℓ >
0, and ℓ sequences ( a k ) , . . . , ( a ℓ − k ), that alternately define the partialnumerators a k : a k = a ( k mod ℓ ) k . The case ℓ = 2 encountered for exp is the mostcommon, but higher values also happen (e.g., ℓ = 4 for ψ ′′ , where ψ = Γ ′ / Γ is thelogarithmic derivative of the Gamma function).This is not a restriction in our approach, by the following.
Lemma 7.
Given a period ℓ > , a sequence ( u k ) k ≥ is P-recursive if and only ifall its subsequences ( u ℓk + j ) k ≥ are P-recursive, for j = 0 . . . ℓ − .Proof. This lemma is classical. We give a constructive proof for completeness.If the sections ( u ℓk + j ) k ≥ are P-recursive, then their generating series s j ( z ) = P k ≥ u ℓk + j z k are D-finite, then so is s ( z ℓ ) + zs ( z ℓ ) + · · · + z ℓ − s ℓ − ( z ℓ ) andtherefore its sequence of coefficients ( u k ) k ≥ is P-recursive. Conversely, if ( u k ) k ≥ is P-recursive, then its generating series s ( z ) is D-finite and so is its Hadamardproduct with z j / (1 − z ℓ ), then also its quotient by z j evaluated at z /ℓ and this isprecisely the generating series of ( u ℓk + j ) k ≥ . (cid:3) In cases like the exponential function, this lemma implies that the sequenceof partial numerators ( a k ) k ≥ itself satisfies a linear recurrence. With the recur-rences (5), this alone does not imply that ( P k ) and ( Q k ) are also P-recursive forin general, no such closure property exists. For instance, the sequence definedby u n := Q nk =1 k ! satisfies a linear recurrence of order 1 with a coefficient ( k !) thatis P-recursive, but ( u n ) itself is not P-recursive, as can be seen from its asymp-totic growth that is too fast. The crucial property in our application is that thesequences ( a ik ) are rational in k . This allows for the following. Lemma 8.
Let a k , . . . , a ℓ − k be rational functions of k and z , a k be defined for k ≥ by a k := a ( k mod ℓ ) k and let the sequences ( P k ) and ( Q k ) be defined by the recur-rence (5) . Then ( P k ) and ( Q k ) are P-recursive sequences.Proof. The proof is constructive. By the recurrences (5) and the definition of a k ,all P ℓk + j for j = 1 , . . . , ℓ rewrite as linear combinations of P ℓk +1 and P ℓk with coef-ficients that are rational functions of z and k . Thus P ℓk + j , P ℓ ( k +1)+ j and P ℓ ( k +2)+ j are linearly dependent for j = 0 , . . . , d − ORMULAS FOR CONTINUED FRACTIONSAN AUTOMATED GUESS AND PROVE APPROACH11
The same reasoning applies to ( P ℓk + j ) k ≥ for any j ∈ { , . . . , ℓ − } and showsthat it is a P-recursive sequence and thefore that so is ( P k ) by the previous lemma.By construction, ( Q k ) satisfies the same recurrence as ( P k ). (cid:3) Example.
The special case ℓ = 2 is important in applications. Starting from P k = P k − + a k P k − and its first two shifts, the linear combination P k +2 + P k +1 − a k +1 P k gets rid of the terms with odd index, leaving:(9) P k +2 = (1 + a k +1 + a k +2 ) P k − a k a k +1 P k − . A similar computation would give a recurrence between the terms with odd index.The proof of Lemma 8 leads to an algorithm in two steps: compute a recurrencefor ( P k ) and ( Q k ) and then appeal to Lemma 7. A simpler and faster computationproceeds directly from a recurrence for ( P ℓk ), thanks to the following. Proposition 9.
Let F ( X, Y ) be a rational function, that is regular at X = Y =0 . Let a ik , i = 1 , . . . , ℓ be rational functions in X and k with positive valuationin X . Let ( a k ) , ( P k ) and ( Q k ) be defined as in the previous lemma and ( H k ) as inTheorem 4. If val H kℓ → ∞ as k → ∞ , then the continued fraction K ∞ m =1 a k is thesolution of Y ′ = F ( X, Y ) with Y (0) = 0 . Here again, other values for Y (0) are obtained by a change of unknown function. Proof.
Since the denominator of F does not vanish at 0, F admits an expansionin power series and thus by Proposition 3, the differential equation Y ′ = F ( X, Y )with Y (0) = 0 possesses a formal power series solution S ( X ).The condition on the valuations of the sequences ( a ik ) makes the continued frac-tion well-defined, in the sense that the sequence of power series ( P k /Q k ) convergesto a power series G ( X ). Thus if its subsequence f k = P kℓ /Q kℓ converges to S ( X ),then we have G ( X ) = S ( X ). Induction from Eq. (5) shows that Q k (0) = 1 forall k ≥ P k (0) = 0 (by the positive valuation of a ). This gives f k (0) = 0and val H kℓ = val( f ′ k − F ( X, f k )). Thus by Proposition 3, the sequence ( f k ) con-verges to S . (cid:3) Example.
The proof for the continued fraction for exp from the introduction goesas follows. Starting from the recurrences for ( P k ) when k is even and when k isodd and proceeding as for Eq. (9) yields P k +2 = P k + z k − P k − , which is also satisfied by Q k since this computation does not depend on the initialconditions.Next, turn to the numerator of the evaluation of y ′ − y − y = P k /Q k ,namely H k = P ′ k Q k − P k Q ′ k − P k Q k − Q k . Using Proposition 5, or directly as in Section 2, leads to a recurrence of order 4,on which reduction of order yields H k +2 = − z H k / (4(2 k + 1) ), which concludesthe proof. Algorithm 3
Discovery and proof of continued fractions
Input: Y ′ = F ( z, Y ) with F ∈ K ( z )[ Y ];a bound N on the number of coefficients to guess from;a bound L on the period to be found. Output:
In case of success, an explicit expression for the continued fraction ex-pansion of the solution such that Y (0) = 0.Compute the first coefficients a , . . . , a N of the continued fraction expansion ofthe power series solution of Y ′ = F ( z, Y ) with Y (0) = 0. for ℓ = 1 , , , . . . , L do Use rational interpolation to compute a ji interpolating the subse-quences ( a ℓi + j ) i , j = 0 , . . . , ℓ − if this has been successful then compute a recurrence R for H kℓ , with H k defined in Theorem 4compute a new recurrence R ′ from R and the initial conditions for ( H k )using Algorithm 1 if R ′ exhibits the increase of (val H ℓk ) thenreturn the rational functions a ji end ifend ifend forreturn FAIL 6.
Experiments
An overview of the whole approach is given in Figure 3. In practice, N =20 and L = 2 have proved sufficient in our experiments except for one case ofperiod 4. For the computation of the first terms of the continued fraction, one caneither compute a power series expansion first, e.g., by Newton iteration [6], or usetechniques for continued fraction expansions of solutions of Riccati equations [9].Our main experimental result is the following. Empirical Observation.
All the 53 explicit C-fractions formulas of the com-pendium by Cuyt et alii [10] can be guessed and proved by our approach and itsvariants below. Among them the vast majority (44) are solutions of Riccati equa-tions, 2 satisfy q -Riccati equations and the remaining 7 satisfy difference equations. We now give more detail on the calculations in the differential case and thenoutline the variants of our method in the q -difference and difference cases. Animplementation in the differential case is available under the form of a submodule gfun:-ContFrac of the package gfun (for versions ≥ . C-fractions from Differential Equations.
In our experiment, the Riccatiequations were themselves found by a guessing approach on power series expansionsto small order (less than 30). Depending on how one decides to define the powerseries from the computer algebra point of view, these Riccati equations can also beautomatically proved to hold.
ORMULAS FOR CONTINUED FRACTIONSAN AUTOMATED GUESS AND PROVE APPROACH13
Gauss’s continued fraction.
The classical hypergeometric series is F ( a, b ; c ; z ) := X n ≥ ( a ) n ( b ) n ( c ) n z n n ! , where ( a ) n is the Pochhammer symbol ( a ) n = a ( a + 1) · · · ( a + n − F ( a, b ; c ; z ) F ( a, b + 1; c + 1; z ) = 1 + ∞ K m =1 a m z ,a k = − ( k + b )( k + c − a )(2 k + c )(2 k − c ) , a k +1 = − ( k + a )( k + c − b )(2 k + c )(2 k + 1 + c ) , for the quotient of two contiguous hypergeometric series. This is the source ofmany continued fractions for special functions by specialization of the parameters.If y = 1 + F is the function on the left-hand side, then elementary propertiesof the F that can be derived from the first order recurrences satisfied by itscoefficients show that cz ( z − y ′ = a ( c − b ) z + ( c ( a − b ) z + c ) y + c y . This is our starting point. From there, it is easy to compute the first 20 coefficientsand conjecture the formulas for a k and a k +1 by rational interpolation. As inEq. (9), a recurrence for even indices follows. From Corollary 6, a linear recurrenceof order 4 follows for the remainder H k , that can be either obtained by hand fromProposition 5, or by a generic code that searches for linear dependency. Next,reduction of order gives a two-term linear recurrence within a couple of seconds: H n − = z ( n + a ) ( n − a + c ) ( n + b ) ( n − b + c )(2 n + c ) (2 n + c − H n − and this concludes the automatic proof. More parameters.
Khovanskii [13, p. 85] gives an explicit continued fractionwith 5 parameters for the power series solution of the differential equation(1 + αz ) zy ′ + ( β + γz ) y + δy = ǫz, y (0) = 0(an extra parameter k is obtained by changing z into z k and adjusting the coefficientof y ′ ; we have relabeled the parameters). This contains the equation for Gauss’scontinued fraction above as a special case.From there again rational formulas for a k and a k +1 are obtained by guessingon the first 20 values; a recurrence of order 4 can be found for the remainders H k ;Algorithm 1 reduces it to the conclusive recurrence:(2 n + β ) (2 n + β − H n − = − ( αn + ( αβ + γ ) n + βγ + δǫ )( αn + ( αβ − γ ) n + δǫ ) z H n − . Other examples.
We also applied our method to a few functions not men-tioned by Cuyt et alii [10] and in particular found (and proved) experimentally thefollowing nice C-fraction for the Airy function: z Ai ′ Ai (1 /z ) = − − z . (cid:18) ∞ K m =2 a m ( z )1 (cid:19) ,a k = (6 k − z / , a k +1 = (6 k + 1) z / . It also follows from known C-fractions for the divergent F .6.2. q-analogues. The method used in this article also applies to q -analogues. Weoutline the very simple example of the q -exponential: e q ( z ) := X m ≥ (1 − q ) m (1 − q )(1 − q ) · · · (1 − q m ) z m , which satisfies the q -differential equation(10) e q ( qz ) − e q ( z )( q − z − e q ( z ) = 0 . The classical exponential is obtained as the limit when q →
1. The first coefficientsof the continued fraction expansion let one guess a = z , a k = − q k − (1 − q ) z (1 + q k )(1 − q k − ) , a k +1 = q k (1 − q ) z (1 + q k )(1 − q k +1 ) , a clear generalization of the continued fraction (2) for exp. In order to prove thiscontinued fraction, the recurrence for ( P k ) is computed as in Section 5.4, whichgives P k +2 = (cid:18) − (1 − q ) q k z (1 + q k )(1 + q k +1 ) (cid:19) P k + q k − (1 − q ) z (1 + q k ) (1 − q k +1 )(1 − q k − ) P k − . The sequence H k is defined as the numerator of the evaluation of (10), namely P k ( qz ) Q k ( z ) − P k ( z ) Q k ( qz )( q − z − P k ( z ) Q k ( qz ) − Q k ( z ) Q k ( qz ) . Next, we compute a linear dependency between H k , H k +2 , . . . , which is still of or-der 4 (but significantly bigger than its limit as q → q -analogue of reductionof order then proves H k +2 = − q k +2 (1 − q ) z (1 + q k +1 ) (1 − q k +1 ) H k , which concludes the proof, providing with a generalization of the expression for expthat is recovered by letting q → q -analogue of Gauss’scontinued fraction [10, 19.2.1]: the q -hypergeometric series is defined by φ ( a, b ; c ; q ; z ) = X n ≥ ( a ; q ) n ( b ; q ) n ( c ; q ) n z n ( q ; q ) n , where ( a ; q ) n is the q -Pochhammer symbol( a ; q ) n = (1 − a )(1 − aq ) · · · (1 − aq n − ) . Heine’s continued fraction is φ ( a, b ; c ; q ; z ) φ ( a, bq ; cq ; q ; z ) = 1 + ∞ K m =1 a m z ,a k +1 = (1 − aq k )( cq k − b ) q k (1 − cq k )(1 − cq k +1 ) , a k = (1 − bq k )( cq k − a ) q k − (1 − cq k − )(1 − cq k ) . ORMULAS FOR CONTINUED FRACTIONSAN AUTOMATED GUESS AND PROVE APPROACH15
We sketch the main steps of the computation. The q -Riccati equation is(1 − c ) F ( z ) F ( qz ) + (1 − c )( bz − c ) F ( qz )+ (1 − z )(1 − az ) F ( z ) − z ( a − b − c ) = 0 . The sequence H k of interest is therefore the numerator of the evaluation of thisleft-hand size at F ( z ) = P k ( z ) /Q k ( z ). The continued fraction being periodic ofperiod 2, a recurrence for H k (or order 4) is computed. Reduction of order yields H k +2 z H k = (1 − aq k +1 )(1 − bq k +1 )( a − cq k +1 )( b − cq k +1 ) q k +1 (1 − cq k +1 )(1 − cq k +2 ) , which concludes the proof. This automates 2 more of the formulas in [10].6.3. Difference Equations.
The same method applies to difference equations.For instance, it results in one of the classical proofs [14, chap. 3] of Brouncker’scontinued fraction for b ( s ) := Γ (cid:0) s +14 (cid:1) Γ (cid:0) s +34 (cid:1) ! , where Γ is Euler’s Gamma function. Using the functional equation Γ( s +1) = s Γ( s ),it follows that b ( s ) satisfies b ( s ) b ( s + 2) = 16 / ( s + 1) . Looking for a formal power series solution in inverse powers of s (and nonnegativeleading term) leads to a unique solution b ( s ) = 4 /s − /s + · · · This is thenconverted into a continued fraction expansion with coefficients ( a k ) given by4 s , s , s , s , s , s , . . . from which it is easy to conjecture a k = (2 k − / (4 s ) for k ≥
3. The analogueof H k in this context is H k = ( s + 1) P k ( s ) P k ( s + 2) − Q k ( s ) Q k ( s + 2) , for which the same approach as above produces a linear recurrence of order 4 whichis not sufficient to conclude that the valuations increase. From there, reduction oforder with Algorithm 1 yields the shorter H k +1 = − (2 k + 1) s ( s + 2) H k , k ≥ , which exhibits the required increase of valuations.The same technique has been applied to all the explicit C-fractions concerningthe ψ function in [10], thereby completing the experiment on this book.7. Conclusion
In a simple and unified way, our approach to continued fractions recovers anunexpectedly large number of explicit C-fractions from the literature. One miraclethat takes place is that in all cases, the sequence of remainder polynomials turnsout to be hypergeometric or q -hypergeometric. We are currently exploring thisphenomenon in more detail. Acknowledgements
This work was supported in part by the project FastRelax ANR-14-CE25-0018-01.
References [1] M. Abramowitz and I. A. Stegun, editors.
Handbook of mathematical functions with formulas,graphs, and mathematical tables . Dover Publications Inc., New York, 1992.[2] G. A. Baker, Jr. and P. Graves-Morris.
Pad´e approximants , volume 59 of
Encyclopedia ofMathematics and its Applications . Cambridge University Press, Cambridge, second edition,1996.[3] B. Beckermann and G. Labahn. A uniform approach for the fast computation of matrix-typePad´e approximants.
SIAM Journal on Matrix Analysis and Applications , 15(3):804–823, July1994.[4] B. Beckermann and G. Labahn. Recursiveness in matrix rational interpolation problems.
Journal of Computational and Applied Mathematics , 77(1-2):5–34, 1997.[5] A. Benoit, F. Chyzak, A. Darrasse, S. Gerhold, M. Mezzarobba, and B. Salvy. The dynamicdictionary of mathematical functions (DDMF). In K. Fukuda, J. van der Hoeven, M. Joswig,and N. Takayama, editors,
ICMS 2010 , volume 6327 of
LNCS , pages 35–41, 2010.[6] R. P. Brent and H. T. Kung. Fast algorithms for manipulating formal power series.
J. ACM ,25(4):581–595, 1978.[7] S. Chen, M. Jaroschek, M. Kauers, and M. F. Singer. Desingularization explains order-degreecurves for ore operators. In
ISSAC’13 , pages 157–164, New York, NY, USA, 2013. ACM.[8] D. Chudnovsky and G. Chudnovsky. Classical constants and functions: Computationsand continued fraction expansions. In D. Chudnovsky, G. Chudnovsky, H. Cohn, andM. Nathanson, editors,
Number Theory , pages 13–74. Springer New York, 1991.[9] K. D. Cooper, S. C. Cooper, and W. B. Jones. More on C -fraction solutions to Riccati equa-tions. In Proceedings of the U.S.-Western Europe Regional Conference on Pad´e Approximantsand Related Topics (Boulder, CO, 1988) , volume 21, pages 139–158, 1991.[10] A. Cuyt, V. B. Petersen, B. Verdonk, H. Waadeland, and W. B. Jones.
Handbook of continuedfractions for special functions . Springer, New York, 2008.[11] J. v. z. Gathen and J. Gerhard.
Modern Computer Algebra . Cambridge University Press,Cambridge u.a., 3 edition edition, June 2013.[12] W. B. Jones and W. J. Thron.
Continued fractions , volume 11 of
Encyclopedia of Mathematicsand its Applications . Addison-Wesley Publishing Co., Reading, Mass., 1980.[13] A. N. Khovanskii.
The application of continued fractions and their generalizations to problemsin approximation theory . P. Noordhoff N. V., 1963.[14] S. Khrushchev.
Orthogonal polynomials and continued fractions , volume 122 of
Encyclopediaof Mathematics and its Applications . Cambridge University Press, Cambridge, 2008.[15] F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, editors.
NIST Handbook ofMathematical Functions . Cambridge University Press, 2010.[16] B. Salvy and P. Zimmermann. Gfun: a Maple package for the manipulation of generatingand holonomic functions in one variable.
ACM Trans. Math. Softw. , 20(2):163–177, 1994.[17] R. P. Stanley.
Enumerative combinatorics , volume 2. Cambridge University Press, 1999.[18] S. P. Tsarev. An algorithm for complete enumeration of all factorizations of a linear ordinarydifferential operator. In
ISSAC ’96 , pages 226–231. ACM Press, 1996. ´ENS de Lyon, LIP (U. Lyon, CNRS, ENS Lyon, UCBL), France
E-mail address : [email protected]
Inria, LIP (U. Lyon, CNRS, ENS Lyon, UCBL), France
E-mail address ::