Symbolic computation of hypergeometric type and non-holonomic power series
aa r X i v : . [ c s . S C ] F e b Symbolic computation of hypergeometric type and non-holonomicpower series
Bertrand Teguia Tabuguia, Wolfram Koepf
University of Kassel, Heinrich-Plett-Str.40. 34132 Kassel, Germany
Abstract
A term a n is m -fold hypergeometric, for a given positive integer m , if the ratio a n + m / a n is a ra-tional function over a field K of characteristic zero. We establish the structure of holonomicrecurrence equation, i.e. linear and homogeneous recurrence equations having polynomial co-efficients, that have m -fold hypergeometric term solutions over K , for any positive integer m .Consequently, we describe an algorithm, say mfoldHyper, that extends van Hoeij’s algorithm(1998) which computes a basis of the subspace of hypergeometric ( m =
1) term solutions ofholonomic recurrence equations to the more general case of m -fold hypergeometric terms.A Laurent-Puiseux series ∞ X n = n a n ( z − z ) n / k ( a n ∈ K , k ∈ N , n ∈ Z ) , (1)where k denotes the corresponding Puiseux number, is mainly characterized by the coefficient a n .We generalize the concept of hypergeometric type power series introduced by Koepf (1992), byconsidering linear combinations of Laurent-Puiseux series whose coefficients are m -fold hyper-geometric terms. Such power series could not be computed before due to the lack of an algorithmto find m -fold hypergeometric term solutions of holonomic recurrence equation which constitutea key step in Koepf’s procedure. Thanks to mfoldHyper, we deduce a complete procedure tocompute these power series; indeed, it turns out that every linear combination of power serieswith m -fold hypergeometric term coefficients, for finitely many values of m , is detected.On the other hand, we investigate an algorithm to represent power series of non-holonomicfunctions like tan( z ) , (1 − tan( z )) / (1 + tan( z )) , z / (exp( z ) − , log(1 + sin( z )) , etc. The algorithmfollows the same steps of Koepf’s algorithm, but instead of seeking holonomic differential equa-tions, quadratic differential equations are computed and the Cauchy product rule is used to de-duce recurrence equations for the power series coefficients. This algorithm defines a normalfunction that yields together with enough initial values normal forms for many power series ofnon-holonomic functions. Therefore, non-trivial identities likeln + tan( z )1 − tan( z ) ! = sin(2 z )1 + cos(2 z ) ! are automatically proved using this approach. This paper is accompanied by implementations inthe Computer Algebra Systems (CAS) Maxima 5.44.0 and Maple 2019. Keywords: m -fold hypergeometric term; holonomic equation; hypergeometric type powerseries; quadratic differential equation; normal form1 . Introduction
The applicability of complex analysis is essentially restricted to analytic functions, since iteasily allows both differentiation and integration. These functions are represented by powerseries with positive radius of convergence. Power series are used to represent orthogonal poly-nomials (Koepf and Schmersau (1998)); in combinatorics, generating functions are power series(Stanley (2011)); in dynamical systems, algebraic properties of power series involve most of theconstructions (see Lubin (1994)); we can also enumerate commutative algebra and algebraic ge-ometry (Brewer (2014)), (Zariski and Samuel, 1960, Chapter VII). It is therefore important toknow the exact general coefficient or formula of a power series. There is no algorithm whichcomputes the power series of every given analytic function. We classify series with a certaincommon property, and build an algorithm which will always find the power series representationfrom such an analytic expression, whenever possible. It is important to notice the word ”expres-sion”, because we are not considering complex functions as abstract objects defined in a certaindomain and its range, but instead as a differentiable object that we can manipulate symbolicallyto characterize its Taylor coefficients by a certain type of linear recurrence equation. Moreover,by the unique power series characterization, this approach does not only lead to the verificationof known identities, but also to the discovery of new ones.Let K be an infinite computable field , and ( a n ) n ∈ Z , a n ∈ K , be an m -fold hypergeometricsequence such that a n + m = r ( n ) a n , ∀ n > n , n ∈ Z , (2)where r ( n ) denotes a rational function in K ( n ), m ∈ N , and n is the first non-zero term index. m -fold hypergeometric sequences are very useful in summation theory (Koepf (2014)). Our firstinterest is to describe an algorithm which computes power series (Puiseux series) of the form ∞ X n = n a n ( z − z ) n / k ( a n ∈ K , k ∈ N , n ∈ Z ) , (3)such that a n is an m -fold hypergeometric term.In 1992, Koepf published an algorithmic approach for computing power series (see Koepf(1992)). The algorithm was implemented in the computer algebra systems (CAS) Maple (Heck(2003)) and Mathematica (Wolfram (2003)). In his original approach, Koepf considered threetypes of functions: two-term recurrence relation type which corresponds to expressions leadingto a linear recurrence equation equivalent to (2). That is Q n a n + m + P n a n = , n ∈ Z , (4)where Q n , P n are polynomials in K [ n ]. The second type called exp-like , corresponding to ex-pressions leading to linear recurrence equations with constant coefficients in K . And the thirdtype with a completely different approach based on partial fraction decomposition (over C ) cor-responding to rational functions in C ( z ). All gathered in the Maple and Mathematica packagesFPS could already recover the power series formulas of a large family of analytic functions. Email addresses: [email protected] (Bertrand Teguia Tabuguia), [email protected] (Wolfram Koepf)
URL: (Bertrand Teguia Tabuguia), (Wolfram Koepf) Mostly K : = Q ( α , . . . , α N ) is the field of rational functions in several variables m -fold hypergeometric terms. Therefore, if we could find all m -fold hypergeometric term solutions of a linear homogeneous recurrence equation, then we couldconsiderably increase the family of power series computed automatically.Marko Petkovˇsek later published an algorithm which finds all hypergeometric ( m = LREtools[hypergeomsols] .An equivalent algorithm is part of our main approach, this is described in (Teguia Tabuguia(2020b)).Note, however, that the Petkovˇsek and van Hoeij algorithms might only find hypergeometricterm solutions in an extension field of Q , which in certain cases, for m >
1, can be equivalentto m -fold hypergeometric term solutions in Q . Indeed, the algorithm is implemented to find allhypergeometric term solutions in Q ( α ), where α ∈ C \ Q ; since α is not always explicitly knownin advance, we will often replace extension fields of Q by C . But this has some disadvantages ofsimplicity. If we consider the power series of the cosine function at z = z ) = ∞ X n = ( − n (2 n )! z n , (5)then we observe that its general coefficient satisfies the recurrence equation(1 + n ) (2 + n ) a n + + a n = . (6)Using Koepf’s algorithm, the type m = a = a = , with the same initial values, we findthe hypergeometric solution i n Γ ( n + + ( − i) n Γ ( n + , i ∈ C , i = − z ) = ∞ X n = ( − i) n + i n Γ ( n + z n . (8) We used Maple 2019 for this paper Q since therecurrence equation obtained is a two-term recurrence relation. In general, an issue occurswith unnecessary algebraic extensions of Q when van Hoeij’s algorithm is used, because itonly looks for hypergeometric term solutions ( m =1). For example, any linear combinationof cos( z ) or sin( z ) with an expression having a hypergeometric general coefficient will havea formula involving (7). > convert(cos(z)+exp(z),FPS); ∞ X k = (cid:16) + i k + ( − i) k (cid:17) z k k ! > convert(log(1+z)+sin(z),FPS); ∞ X k = − ( − k + k + − i · i k + k + + i · ( − i) k + k + ! z k + Note, however, that this paper does not aim to find the power series formula with a simplehypergeometric general coefficient, but to find the formula with the simplest m -fold hyper-geometric general coefficients. Simple here means that the coefficients are not taken in anextension field of Q whenever there exists an m -fold equivalent over Q . We should highlight m -fold hypergeometric, because up to now there is no implemented algorithm able to findsuch solutions of a linear recurrence equation. And it is worth to have such an algorithmsince in many cases, Maple’s convert command fails to find power series of this type. > convert(arcsin(z)+cos(z),FPS); arcsin( z ) + cos( z ) > convert(exp(zˆ2)+log(1+zˆ3),FPS); e z + ln(1 + z )The above Maple failures rely on the incapacity of van Hoeij’s algorithm to detect m -fold( m >
1) hypergeometric term solutions of so-called holonomic recurrence equations, that ishomogeneous linear recurrence equations with polynomial coefficients. Indeed, by using theMaple package
FormalPowerSeries we get the following holonomic recurrence equations. > bind(FormalPowerSeries) > RE1:=SimpleRE(arcsin(z)+cos(z),z,a(n));
RE1 := − n ( n − n + n − a ( n ) + ( n − a ( n − + ( n − n − n + n − a ( n − + n + n + n + n − a ( n + − n + n + n + n + a ( n + = > RE2:=SimpleRE(exp(zˆ2)+log(1+zˆ3),z,a(n)); Maple’s convert command uses Koepf’s original approach followed by an invocation of van Hoeij’s algorithm. − n − a ( n − + n − n − a ( n − − n − n − a ( n − + n − n − a ( n − + n − n − n + a ( n − − n − n − a ( n − + ( n − n − n − a ( n − + n − n − a ( n − − ( n − n − n + a ( n + = > LREtools[hypergeomsols](RE1,a(n), {} ,output=basis); " ( − i) n Γ ( n + , i n Γ ( n + > LREtools[hypergeomsols](RE2,a(n), {} ,output=basis); ( − n n , (cid:18) − i √ (cid:19) n n , (cid:18) + i √ (cid:19) n n showing that the general coefficients of arcsin( z ) in RE1 and the one of exp( z ) in RE2 are missed.Although some algorithms for computing m -fold hypergeometric term solutions of holo-nomic recurrence equations have been described, none of them is implemented. For example, in(Cluzeau and van Hoeij (2006)) and (Van Hoeij (1999)) an algorithm using linear operators isdeveloped, but the described approach needs non-commutative factorization for its implementa-tion. In our approach however, non-commutative algebra is not needed. We will use a differentview on holonomic recurrence equations and develop a new algorithm to detect all their m -foldhypergeometric term solutions. Thus with our Maple and Maxima implementations, the is-sue with m -fold hypergeometric term solutions of holonomic recurrence equations is solved asdemonstrated in Maxima below (see Teguia Tabuguia (2020a)). (%i1) RE1:FindRE(asin(z)+cos(z),z,a[n]); (% o1 ) − · (1 + n ) · (2 + n ) · (3 + n ) · (4 + n ) · a n + + · (1 + n ) · (2 + n ) · (cid:16) − + · n + n (cid:17) · a n + − n · (cid:16) − + · n − · n + n (cid:17) · a n + ( n − · (cid:16) − + · n − · n + n (cid:17) · a n − + ( n − · a n − = (%i2) mfoldHyper(RE1,a[n]); (% o2 ) "" , ( ( − n (2 · n )! , n · n ! n · (2 · n )! ) (%i3) RE2:FindRE(exp(zˆ2)+log(1+zˆ3),z,a[n]); We used Maxima 5.44.0 for this paper. o3 ) − ( n − · ( n − · (1 + n ) · a n + + · ( n − · ( n − · a n − + ( n − · ( n − · ( n − · a n − − · ( n − · ( n − · a n − + · ( n − · (cid:16) − · n + · n (cid:17) · a n − + · ( n − · ( n − · a n − − · ( n − · (2 · n − · a n − + · ( n − · ( n − · a n − − · ( n − · a n − = (%i4) mfoldHyper(RE2,a[n]); (% o4 ) "" , ( ( − n n ) , " , ( n ! ) , " , ( ( − n n ) (%i5) FPS(asin(z)+cos(z),z,n); (% o5 ) ∞ X n = (2 · n )! · z + · n (2 · n + · n · n ! + ∞ X n = ( − n · z · n (2 · n )! (%i6) FPS(exp(zˆ2)+log(1+zˆ3),z,n); (% o6 ) ∞ X n = ( − n · z · (1 + n ) n + + ∞ X n = z · n n !Another important issue that we solve is the step which consists in deducing, when it exists,the correct linear combination needed to find the hypergeometric type representation sought. Let P ( z ) , P ( z ) , . . . , P d ( z ) be d + K ( z ), and f ( z ) , . . . , f d ( z ) some analytic expres-sions that have m -fold hypergeometric term coefficients in their power series expansions. Moregenerally, our algorithm handles formal series formulas of expressions of the form P ( z ) + d X j = P j ( z ) f j ( z ) . (9)The output of such an input is of course a linear combination of hypergeometric type series,plus a polynomial which might be zero. If the correct linear combination of m -fold hypergeomet-ric term solutions of the corresponding linear recurrence equation is not found, then the outputmight be missed. This happens sometimes with Maple for the hypergeometric ( m =
1) case. Forexample, Maple gives > convert((z+zˆ2+1)*exp(z)+(zˆ3+3)*log(1+z),FPS); ( z + z + e z + ( z +
3) ln( z + > convert(1+z+zˆ2+zˆ3*arctan(z),FPS); + z + z + z · arctan( z )whereas our algorithm yields correctly (%i7) FPS((z+zˆ2+1)*exp(z)+(zˆ3+3)*log(1+z),z,n); % o7 ) 8 · z + z + · z + + ∞ X n = − (cid:16) − − · n − · n − · n − n − ( − n · (4 + n )! + · n · ( − n · (4 + n )! (cid:17) · z + n ( n + · ( n + · (4 + n )! (10) (%i8) FPS(1+z+zˆ2+zˆ3*atan(z),z,n); (% o8 ) ∞ X n = ( − n · z · n · n − + z + > LREtools[hypergeomsols](SimpleRE((zˆ2+z+1)*exp(z)+(zˆ3+3) > *log(1+z),z,a(n)),a(n), {} ,output=basis); " ( − n (2 n − n − n , ( n + Γ ( n + > LREtools[hypergeomsols](SimpleRE(1+z+zˆ2+zˆ3*arctan(z), > z,a(n)),a(n), {} ,output=basis); " i n n − , ( − n n − but the power series terms are missed by the Maple command convert . We mention that thisissue is not related to an argument of convert which has to be specified, in particular the orderof the differential equations involved in the computations. Indeed the default value used for theupper bound of the differential equations sought for power series computations is 4. However,using our Maxima procedure HolonomicDE which implements a variant of Koepf’s algorithm tocompute holonomic differential equations, one finds the following differential equations of orderless than 4. (%i9) HolonomicDE((z+zˆ2+1)*exp(z)+(zˆ3+3)*log(1+z),F(z)); (% o9 ) (1 + z ) · (cid:16) + · z − · z − · z − · z + · z + · z − · z + z + z (cid:17) · d d z · F ( z ) ! − (cid:0) + · z − · z − · z + · z + · z + · z + · z + · z + z (cid:1) · d d z · F ( z ) ! + (cid:0) − − · z − · z − · z + · z + · z − · z + · z + · z + · z (cid:1) · dd z · F ( z ) ! − · (cid:16) − − z + z (cid:17) · (cid:0) − − · z + · z + · z + · z + · z + · z (cid:1) · F ( z ) = %i10) HolonomicDE(1+z+zˆ2+zˆ3*atan(z),F(z)); (% o10 ) z · (cid:16) + z (cid:17) · (cid:16) + · z + · z + · z (cid:17) · d d z · F ( z ) ! − · (cid:16) + · z + · z + · z + · z + · z (cid:17) · dd z · F ( z ) ! + · (cid:16) + z + · z + z (cid:17) · F ( z ) = convert cannot find thepower series formulas of (1 + z + z ) exp( z ) + ( z +
3) log(1 + z ) and 1 + z + z + z arctan( z ) isthat the linear combinations of hypergeometric term solutions of the corresponding holonomicrecurrence equations are missed.As observed with the previous computations, our implementation is written in the CAS Max-ima whose internal command powerseries dedicated to power series computations is rather lim-ited. Indeed, this command is based on a pattern matching instead of algorithmic model. Thesyntax is powerseries(expr,z,z ) that calculates the power series formula of expr with respect tothe variable z at the point of development z . Below are some examples showing certain arisingissues with the command powerseries that are solved by our implementation.• Power series written as a square of a power series. (%i11) powerseries(asin(z)ˆ2,z,0); (% o11 ) ∞ X i1 = genfact (2 · i1 − , i1 , · z + · i1 (2 · i1 + · genfact (2 · i1 , i1 , (%i12) FPS(asin(z)ˆ2,z,n); (% o12 ) ∞ X n = n · n ! · z + · n ( n + · (1 + · n )!• Non-classical power series not detected. (%i13) powerseries((1-sqrt(1-4*z))/2,z,0); (% o13 ) powerseries − √ − · z , z , ! (%i14) FPS((1-sqrt(1-4*z))/2,z,n); (% o14 ) ∞ X n = (2 · n )! · z + n ( n + · n ! (%i15) powerseries(asech(z),z,0); (% o15 ) ∞ X i1 = z i1 · (cid:18) d i1 d z i1 · asech ( z ) (cid:12)(cid:12)(cid:12)(cid:12) z = (cid:19) i1 !8 %i16) FPS(asech(z),z,n); (% o16 ) ∞ X n = − − − n · (1 + · n )! · z + · n (1 + n ) · n ! − log ( z ) + log (2)Observe that despite the general rule used for this latter example, the output given by powerseries is wrong since the logarithmic term log( z ) does not allow the computations ofderivatives at 0.• Power series written as multiplication of two power series. (%i17) powerseries(exp(z)*cos(z),z,0); (% o17 ) ∞ X i4 = z i4 i4 ! · ∞ X i4 = ( − i4 · z · i4 (2 · i4 )! (%i18) FPS(exp(z)*cos(z),z,n); (% o18 ) ∞ X n = − ( − n · n · z + · n (cid:16) (cid:17) n · (cid:16) (cid:17) n · (2 · n + · (4 · n + · (4 · n + · n · (2 · n )! + ∞ X n = ( − n · n · z + · n (cid:16) (cid:17) n · (cid:16) (cid:17) n · (4 · n + · n · (2 · n )! + ∞ X n = ( − n · n · z · n (cid:16) (cid:17) n · (cid:16) (cid:17) n · n · (2 · n )!• A bug due to the involvement of complex numbers in the expansion. (%i19) powerseries(log(1+z+zˆ2),z,0); sign: argument cannot be imaginary; found %i– an error. To debug this try: debugmode(true); (%i20) FPS(log(1+z+zˆ2),z,n); (% o20 ) ∞ X n = − · cos (cid:16) · π · (1 + n )3 (cid:17) · z + n n + C (extension field of Q involving i and some irrational numbers) of thecorresponding linear recurrence equation.On the other hand, some expressions like tan( z ) , sec( z ) , csc( z ) , etc . do not lead to linear recur-rence equations, although they are analytic in certain domains. Therefore, we should investigatetheir power series computation. For that purpose, in this paper we consider two approaches.9ur second approach is to follow the same procedure as Koepf, but this time, instead oflooking for a linear differential equation, we look for quadratic ones. For example, for the tangentfunction, one can find the homogeneous differential equation d d z · F ( z ) − · F ( z ) · dd z · F ( z ) ! = , (11)which after the use of the Cauchy product rule, will lead to the recurrence equation(1 + n ) · (2 + n ) · a n + − · n X k = ( k + · a k + · a n − k = z ) case, we will find the recurrence equation − n X k = (cid:16)(cid:16) − · k (cid:17) · a k + + (2 · k + · a k + · n (cid:17) · a n − k + + (cid:16) a k + (cid:16) − k − · k − (cid:17) · a k + (cid:17) · a n − k = z ) = ∞ X n = ( − n − n (cid:16) n − (cid:17) B n z n − (2 n )! , (14)sec( z ) = ∞ X n = ( − n E n z n (2 n )! , (15)are not explicit because of the unknowns B n and E n which represent, respectively, Bernoulli andEuler numbers. Those numbers themselves satisfy rather complicated non-holonomic recurrenceequations.In our third approach, we extend our algorithm of hypergeometric type series. Here we con-sider reciprocals of formal power series and build an algorithm which can compute reciprocalsof power series of some analytic expressions. Using Cauchy’s product rule, some other powerseries are also deduced. This is only available in our Maxima package. (%i21) FPS(tan(z),z,n); (% o21 ) ∞ X n = n X k = A k · ( − n − k (1 − · k + · n )! · z + · n , A k = k X j = − ( − j · A k − j (2 · j )! , A = (%i22) FPS(sec(z),z,n); (% o22 ) ∞ X n = A n · z · n , A n = n X k = − ( − k · A n − k (2 · k )! , A = + tan( z )1 − tan( z ) = exp · arctanh sin(2 z )1 + cos(2 z ) !! , | z | < , (16)which cannot be recognized without using non-trivial transformations (see (Koepf, 2006, Chapter9)). With our Maxima implementations, computing quadratic differential equations for both sidesyields two compatible differential equations as shown below. (%i23) DE1:QDE((1+tan(z))/(1-tan(z)),F(z),Inhomogeneous); (% o23 ) dd z · F ( z ) − F ( z ) − = (%i24) DE2:QDE(exp(2*atanh(sin(2*z)/(1+cos(2*z)))),F(z)); (% o24 ) F ( z ) · d d z · F ( z ) ! − · dd z · F ( z ) ! · d d z · F ( z ) ! + · F ( z ) · dd z · F ( z ) ! = (%i25) CompatibleDE(DE1,DE2,F(z)); The two differential equations are compatible (% o25 ) true Moreover, our FPS algorithm simplifies the difference to zero in a neighborhood of 0. (%i26) FPS((1+tan(z))/(1-tan(z))-exp(2*atanh(sin(2*z)/(1+cos(2*z)))),z,n); (% o26 ) 0Things are a little simpler with Maple, but this is because Maple automatically applies somesimplification on the right-hand side so that the same differential equation is found. > FPS[QDE](exp(2*arctanh(sin(2*z)/(1+cos(2*z)))),F(z)) − dd z F ( z ) ! F ( z ) + d d z F ( z ) = Taylor with the same syntaxas taylor(f,z,z ,d) which computes the Taylor expansion of order d of f ( z ). And it turns out asexpected that our Taylor command is clearly asymptotically faster than taylor for holonomicfunctions. As an example we have: Two differential equations are said to be compatible if every solution of the lower order DE is solution of the other. %i27) taylor(sin(z)ˆ2,z,0,10); (% o27 ) / T / z − z + · z − z + · z + · · · (%i28) Taylor(sin(z)ˆ2,z,0,10); (% o28 ) 2 · z − z + · z − z + z that illustrates the coincidence between both outputs. Testing the efficiency for large order gives: (%i29) taylor(sin(z)ˆ2,z,0,1000)$ Evaluation took 15.8500 seconds (19.6100 elapsed) (%i30) Taylor(sin(z)ˆ2,z,0,1000)$
Evaluation took 1.8300 seconds (1.8900 elapsed)which shows that, asymptotically, our
Taylor command takes just about a fraction of Maxima’sinternal taylor computation timing for sin( z ) .As already observed, most of our computations will be presented with the CAS Maximawhich was the main system used when our results were being developed (see Teguia Tabuguia(2020a)). The further sections are organized as follows.The next section recalls the two first steps in Koepf’s algorithm: computing holonomic differ-ential equations and holonomic recurrence equations. In this section, we add some linear algebratricks in order to gain more efficiency in the process of getting holonomic differential equationsfor hypergeometric type functions. This section ends with the description of our asymptoticallyfast algorithm for computing Taylor expansions of holonomic functions.Section 3 is devoted to our most important result, which is to present a complete algorithmto find all m -fold hypergeometric term solutions of linear recurrence equations with polynomialcoefficients.In Section 4 we complete Koepf’s algorithm with mfoldHyper. We will see in this sectionhow our algorithm handles the starting points, the polynomial part, and the Puiseux numbersinvolved in a given hypergeometric type series expansion.Section 6 describes our approach based on the computation of quadratic differential equationsfor representing non-holonomic power series. This part is another important contribution of ourwork.
2. Computing holonomic equations
This section is about computing holonomic differential equations from given expressions anduse them to deduce holonomic recurrence equations of their Taylor or power series coefficientsprior to computing m -fold hypergeometric term solutions or Taylor expansions. The first interestof this section is our strategy to implement the given algorithm in (Koepf (1992)) for comput-ing holonomic differential equations, which is slightly more efficient than its original version.Moreover, we give an algorithm to check whether two holonomic differential equations are com-patible. Furthermore, we propose a formal algorithm that performs fast computation of largerorder Taylor expansions for holonomic functions as discussed in (Koepf, 2006, Section 10.27).The algorithm to deduce the recurrence equations from holonomic differential equations will begiven just because it contains a rewrite rule that plays an asset role in Section 6.12e recall Koepf’s original approach from an introductory example. Let f ( z ) = exp( z ) + cos( z ).We want to find a holonomic differential equation (DE) with coefficients in Q [ z ] satisfied by f ( z ) and deduce a holonomic recurrence equation (RE) with coefficients in Q [ n ] satisfied by theTaylor coefficients a n of f . Search of a holonomic DE: f ′ ( z ) = exp( z ) − sin( z ), and therefore there is no A ( z ) ∈ Q ( z )such that f ′ ( z ) + A ( z ) f ( z ) = A ( z ) should be − (exp( z ) − sin( z )) / (exp( z ) + cos( z )) whichis not rational. Therefore we move to the second order. We search for A ( z ) , A ( z ) ∈ Q ( z ) suchthat f ′′ ( z ) + A ( z ) f ′ ( z ) + A ( z ) f ( z ) = . Note that at this step f ′ ( z ) is computed again. We write the sum in terms of linearly independentparts and we obtain(1 + A ( z ) + A ( z )) exp( z ) + ( A ( z ) −
1) cos( z ) − A ( z ) sin( z ) = , and we get the linear system A ( z ) − = A ( z ) = A ( z ) + A ( z ) + = , which has no solution. However, for the third order, the relation f (3) ( z ) + A ( z ) f (2) ( z ) + A ( z ) f (1) ( z ) + A ( z ) f ( z ) = , with f (3) ( z ) = exp( z ) + sin( z ), leads to a solvable system. By writing the sum in terms of linearlyindependent parts, we get(1 + A ( z ) + A ( z ) + A ( z )) exp( z ) + ( A ( z ) − A ( z )) cos( z ) + (1 − A ( z )) sin( z ) = , (17)so − A ( z ) = A ( z ) − A ( z ) = + A ( z ) + A ( z ) + A ( z ) = ⇐⇒ A ( z ) = A ( z ) = − A ( z ) = . (18)Hence the corresponding holonomic DE d d z · F ( z ) − d d z · F ( z ) + dd z · F ( z ) − F ( z ) = . (19)for f ( z ) is valid.Notice that if some of the coefficients found in (18) did have polynomials different from 1as their denominators, then a further step would be the multiplication of the resulting holonomicDE with the least common multiple of the denominators. Conversion of (19) into its corresponding RE:
We set f ( z ) = ∞ X n = a n z n , f ′ ( z ) = ∞ X n = na n z n − = ∞ X n = ( n + a n + z n , (20) f ′′ ( z ) = ∞ X n = n ( n − a n z n − = ∞ X n = ( n + n + a n + z n , (21) f (3) ( z ) = ∞ X n = n ( n − n − a n z n − = ∞ X n = ( n + n + n + a n + z n . (22)By substitution of these identities in (19) for f , we get0 = f (3) ( z ) − f ′′ ( z ) + f ′ ( z ) + f ( z ) = ∞ X n = ( n + n + n + a n + z n − ∞ X n = ( n + n + a n + z n + ∞ X n = ( n + a n + z n − ∞ X n = a n z n = ∞ X n = [( n + n + n + a n + − ( n + n + a n + + ( n + a n + − a n ] z n , hence by equating the coefficients we find the holonomic RE( n + n + n + a n + − ( n + n + a n + + ( n + a n + − a n = , n = , , , . . . . (23)for a n . For more details about this procedure see (Koepf (1992); Gruntz and Koepf (1995)). Note that the algorithm we consider here often finds the holonomic DE of lowest order.In this paragraph we look at the algorithm in a different way. Indeed, if we consider f ( z ) = cos( z ) + sin( z ), it is clear that the differential equation that the algorithm will find is a null linearcombination of the derivatives of f ( z ) expanded in the basis (cos( z ) , sin( z )). Thus, we can savethe time spent by computing all the derivatives at each iteration of the original algorithm, bytrying to write each derivative in the same basis. Therefore the computations of the originalalgorithm that are done for all the derivatives at each iteration are reduced to only one derivative. K denotes a field of characteristic zero. Let ( A , A , . . . , A N − ) ∈ K ( z ) N , N ∈ N such that ananalytic expression f satisfies F (cid:16) f , f ′ , . . . , f ( N − , f ( N ) (cid:17) = f ( N ) + A N − · f ( N − + · · · + A · f + A f = . (24)We consider a basis ( e , e , . . . , e l ) of the linear span of all linearly independent summandsover K ( z ) that appear in the complete expansions of the derivatives f , f ′ , . . . , f ( N ) . For example,assume for 0 i , j N , that f ( i ) = e i , + · · · + e i , k i , f ( j ) = e j , + · · · + e j , k j , k i and k j , such that e i , u / e i , v < K ( z ) for all u , v ∈ J , k i K and e j , u / e j , v < K ( z ) for all u , v ∈ J , k j K . Then for f ( i ) and f ( j ) we consider a basis of the linear span of { e i , , . . . , e i , k i , e j , , . . . , e j , k j } which may have less elements since some e i , u , u ∈ J , k i K and e j , v , v ∈ J , k j K can be linearly dependent.Thus each derivative f ( j ) , j ∈ N > (cid:16) f (0) = f (cid:17) can be seen as a vector in the linear space h e , e , . . . , e l i .Since F (cid:16) f , f ′ , . . . , f ( N − , f ( N ) (cid:17) = ⇐⇒ − f ( N ) = A · f + A · f ′ + · · · + A N − · f ( N − , (25)we can write in a matrix representation − f ( N ) = h f , f ′ , . . . , f ( N − i ( e , e ,..., e l ) ( A , A , . . . , A N − ) T . (26)Therefore, one sees that seeking for a holonomic DE of order N satisfied by a given expres-sion f ( z ) is equivalent to find a basis in a K ( z )-linear space where the system (cid:16) f ( N ) ( z ) , f ( N − ( z ) , . . . , f ′ ( z ) , f ( z ) (cid:17) is linearly dependent. The idea of the method described in this section is to construct such a basiswhile computing each derivative of f ( z ) and their components. Thus, in each iteration N , if allthe N + Now let us present the algorithm in general. Consider an expression f ( z ) which is not identi-cally zero with l linearly independent sub-terms over K ( z ). Then we can write f ( z ) = f ( z ) + f ( z ) + · · · + f l ( z ) (27)with f i ( z ) / f j ( z ) < K ( z ) , i , j l . f ( z ) is seen as a vector in the basis E = (cid:0) e , e , . . . , e l (cid:1) where e i = f i . Then we compute the first derivative of f ( z ), and we get the following twopossibilities:• either f ′ ( z ) is expressed in E , which means that there exist α , i = α , i ( z ) ∈ K ( z ) , i = , . . . , l such that f ′ ( z ) = α , e + α , e + . . . + α , l e l . (28)Here in the worst case, f ′ ( z ) and f ( z ) are linearly independent, but then we know that allderivatives can be expanded in E .• Or f ′ ( z ) is not expanded in E , which means that E has to be augmented and there exist α , i ∈ K ( z ) , i = , . . . , l and an integer l > l such that f ′ ( z ) = α , e + α , e + . . . + α , l e l + e l + + . . . + e l . (29)Observe here that the new basis is E = (cid:0) e , . . . , e l (cid:1) with e l + , . . . , e l corresponding toindependent terms brought by f ′ ( z ). And also α , i , i l could be zero. For l , k ∈ Z , l < k we define J l , k K : = { l , l + , . . . , k } f ( z ) satisfies a DE of order N >
1. It follows that the processwill lead to the representation f ( z ) = e + . . . + e l (30) f ′ ( z ) = α , e + · · · + α , l e l + e l + + · · · + e l (31) f ′′ ( z ) = α , e + · · · + α , l e l + α , l + e l + + · · · + α , l e l + e l + + . . . + e l (32) · · · (33) f ( N − ( z ) = α N − , e + · · · + α N − , l N − e l N − + e l N − + + · · · + e l N − (34) f ( N ) = α N , e + · · · + α N , l N − e l N − , (35)with positive integers l l . . . l N − , and α i , j ∈ K ( z ) , i = , . . . , N , j = , . . . , l i − . Note, however, that in this step N only f ( N ) ( z ) is computed by differentiating f ( N − ( z ). Ineach step, the algorithm keeps the coefficients α N , i , the augmented basis and the current deriva-tive.It is straightforward to see that the final basis considered is E N − = (cid:0) e , . . . , e l N − (cid:1) . The algo-rithm keeps information in a matrix form, say H , and at the N th iteration we have H = · · · · · · · · · · · · · · · α , · · · α , l · · · · · · · · · · · · α , · · · α , l · · · · · · α , l · · · · · · · · · ... ... ... ... ... ... ... ... ... ... ... ...α N − , · · · · · · · · · · · · · · · · · · · · · α N − , l N − · · · α N , · · · · · · · · · · · · · · · · · · · · · α N , l N − α N , l N − + · · · α N , l N − . (36) H is a ( N + × l N − matrix in K ( z ), and it contains all information that we need to find the holo-nomic DE sought. Indeed, one can easily show that the coefficients A i ( z ) ∈ K ( z ) , i = , . . . , N − A · v = b , (37)with A = α , α , . . . α N − , ... ... ... ... ... α , l α , l . . . α N − , l α , l + . . . α N − , l + ... ... ... ... ... . . . (38)and b = − α N , α N , ...α N , l N − . (39)16bserve that b is the negative (note the minus in front) of the transpose of the last row of H , and A is the transpose of H deprived of its last row. The above linear system has l N − linear equationsand N unknowns.Let us see how this algorithm works on some examples. Example 1. • f ( z ) = sin( z ) + z cos( z ) . We have two linearly independent terms over Q ( z ) , and we canwrite f ( z ) = e + e , with e = sin( z ) and e = z cos( z ) . Computing the first derivative, we getf ′ ( z ) = − z sin( z ) + z ) = − z · e + z · e . At this step we have H = " − z z , and we get the system " v = " z − z , which has no solution v ∈ Q ( z ) (seen as a one-dimensional vector space). Now we compute the second derivative, and we getf ′′ ( z ) = − z ) − z cos( z ) = − · e − e . H becomes H = − z z − − , which gives the system " − z z v = " , v ∈ Q ( z ) , and we get the solution ( + z + z , − z + z !) . (40) The differential equation sought is therefore (2 + z ) f ′′ ( z ) − z f ′ ( z ) + (6 + z ) f = . (41)• f ( z ) = arctan( z ) . We have only one term so e = arctan( z ) . For the first derivativef ′ ( z ) = + z = · e + e , where e = / (1 + z ) . Since the basis has been augmented there is no system to be solved,and at this step we have H = " . he second derivative givesf ′′ ( z ) = − z (1 + z ) = · e − z + z · e , and we get H = − z + z which produces the system " v = " z + z , v ∈ Q ( z ) . We get v = (cid:16) , z / (1 + z ) (cid:17) , hence the holonomic DE (1 + z ) f ′′ ( z ) + z f ′ ( z ) = . (42)• f ( z ) = exp( z ) + log(1 + z ) = e + e , with e = exp( z ) and e = log(1 + z ) . The first derivativeyields f ′ ( z ) = exp( z ) + + z = e + · e + e , with e = / (1 + z ) . Since a new term is added to the basis, the next step is to compute thesecond derivative f ′′ ( z ) = exp( z ) − + z ) = e + · e − + z ) · e . No term is added to the basis. We try to solve the resulting system. At this stageH = − + z , and we get the system v = − + z , v ∈ Q ( z ) , which has no solution. We move on and compute the third derivativef (3) ( z ) = exp( z ) + + z ) = e + · e + + z ) · e . Thus H = − + z + z ) , nd we obtain the system − + z v = − − + z ) , whose solution in Q ( z ) is ( , − + z (1 + z )(2 + z ) , − − + z + z (1 + z )(2 + z ) !) . (43) Therefore we get the holonomic DE (1 + z )(2 + z ) f (3) ( z ) − ( − + z + z ) f ′′ ( z ) − (3 + z ) f ′ ( z ) = . (44)We implemented this algorithm in Maple and Maxima as HolonomicDE(f,F(z)) to compute aholonomic DE with the indeterminate
F(z) for an expression f of the variable z . In Maxima thepackage contains a global variable Nmax which can be changed in order to look for higher orderdifferential equations. Here are some examples. (%i1) HolonomicDE(asin(z),F(z)); (% o1 ) ( z − · (1 + z ) · d d z · F ( z ) ! + z · dd z · F ( z ) ! = (%i2) HolonomicDE(cos(z)*log(1+z),F(z)); (% o2 ) (1 + z ) · (1 + · z ) · (3 + · z ) · d d z · F ( z ) ! + · (1 + z ) · (cid:16) + · z + · z (cid:17) · d d z · F ( z ) ! + · z · (2 + z ) · (cid:16) + · z + · z (cid:17) · d d z · F ( z ) ! + · (1 + z ) · (cid:16) + · z + · z (cid:17) · dd z · F ( z ) ! + (cid:16) − + · z + · z + · z + · z (cid:17) · F ( z ) = (%i3) HolonomicDE(sin(z)ˆ4*asin(z),F(z)); (% o7 ) false (%i4) Nmax:10$(%i5) HolonomicDE(sin(z)ˆ4*asin(z),F(z))$ Evaluation took 1.7500 seconds (1.7570 elapsed) using 858.342 MB.The latter is a big differential equation of order 10 (default value) > Nmax , that is why thevalue of
Nmax was changed to 10. The following table shows the efficiency gain of this methodon the original one implemented in Maple under the commands
DEtools[FindODE] (Approach2) and
FormalPowerSeries[HolonomicDE] (Approach 3). The difference between these twoMaple commands is that
DEtools[FindODE] uses some simplification in order to consider somespecial functions. Our Maple implementation is accessible through
FPS[HolonomicDE] (Ap-proach 1). 19 able 1: Comparison of timings for computing holonomic DEs f ( z ) CPU timeApproach 1 Approach 2 Approach 3sin( z ) arcsin( z ) 0 .
703 1 .
578 0 . z ) arcsin( z ) .
406 259 .
250 254 . z ) + sin( z ) + log(1 + z ) .
422 11 .
859 11 . z +
3) cos( z ) + log(1 + z ) + cos( z ) sinh( z ) .
968 1466 .
234 1373 . z ) + sinh( z ) .
125 318 .
563 418 . z + z + z ) log(1 + z + z ) + arctanh( z ) cos( z ) .
171 332 .
953 365 . After expanding a holonomic differential equation, one easily establishes the following rewriterule needed to find the corresponding recurrence equation of the underlying power series coeffi-cients (see Koepf (1992)). z l d j dz f −→ ( n + − l ) j · a n + j − l . (45)We used the so-called Pochhammer symbol or shifted factorial defined as ( x ) = x ) n = x · ( x + · · · ( x + n −
1) for a constant x and a non-negative integer n . The final holonomic REsought is obtained after collecting similar terms. To get simpler results, we finally factorize thecoefficients.Our Maxima package contains the function DEtoRE(DE,F(z),a[n]) which converts the holo-nomic differential equation DE depending on the variable z into its corresponding recurrenceequation for the coefficients a[n] . Example 2. (%i1) DE:HolonomicDE(asin(z),F(z))$(%i2) DEtoRE(DE,F(z),a[n]); (% o2 ) n · a n − (1 + n ) · (2 + n ) · a n + = FindRE(f,z,a[n]). (%i3) FindRE(cos(z)+sin(z),z,a[n]); (% o6 ) (1 + n ) · (2 + n ) · a n + + a n = .3. Computing larger-order Taylor expansions of holonomic functions Hypergeometric type functions are strictly contained in the family of holonomic functions.Indeed, it is proved that linear combinations and products of holonomic functions are also holo-nomic (Koepf (1997); Stanley (1980)). Although power series expansions of linear combinationsof hypergeometric type functions remain accessible through the algorithm of Section 4, it is notgenerally the case with their products. This is because algorithm mfoldHyper cannot find ex-plicit formulas for the coefficients of power series expansions of certain holonomic functions.Nevertheless, as we are able to find recurrence equations for the coefficients, the use of enoughinitial values coupled with their corresponding holonomic REs uniquely characterizes their Tay-lor coefficients in a certain neighborhood. It is thanks to this observation that Koepf proceededin computing Taylor polynomials of holonomic expressions by using the outputs of
FindRE (see(Koepf, 2006, Chapter 10)). We are going to use
FindRE to develop an algorithm to computeTaylor polynomials of holonomic functions and compare the result with Maxima’s internal com-mand taylor . First, we give some particular normal forms for holonomic functions.
By an application of the well-known Cauchy-Lipschitz (also called Picard-Lindel¨of) theorem(see (Teschl, 2012, Theorem 2.2)) for uniqueness, the holonomic differential equation of lowestorder and enough initial values corresponding to a holonomic function can be used for identifi-cation purposes. Therefore, such a representation constitutes a normal form (see (Geddes et al.,1992, Chapter 3)). However, our procedure to compute holonomic DEs reduces this normal formdefinition of functions to expressions, because it might happen that equivalent expressions havetwo different representations. Consider the Chebyshev polynomials for example. For | x | <
1, thefollowing two differential equations define the same function. (%i1) DE1:HolonomicDE(cos(4*acos(x)),F(x)); (% o1 ) ( x − · (1 + x ) · d d x · F ( x ) ! + x · dd x · F ( x ) ! − · F ( x ) = (%i2) DE2:HolonomicDE(8*xˆ4-8*xˆ2+1,F(x)); (% o2 ) (cid:16) − · x + · x (cid:17) · dd x · F ( x ) ! − · x · (cid:16) · x − (cid:17) · F ( x ) = HolonomicDE does not use simplifications on its input ex-pressions. However, one can easily prove that these two differential equations are compatibleby substituting the lower order differential equation into the larger one. We implemented such aprocedure in Maxima as
CompatibleDE(DE1,DE2,F(x)) . (%i3) CompatibleDE(DE1,DE2,F(x)); The two differential equations are compatible(% o3 ) trueNext we give the representation that we use as a first step towards computing Taylor ex-pansions. Given an analytic expression f ( z ) at z whose Taylor coefficients satisfy a holonomicrecurrence equation of the form P d ( n ) · a n + d + P d − ( n ) · P n + d − + · · · + P ( n ) · a n = , n ∈ Z , d ∈ N (46)21ith P ( n ) · P d ( n ) , , ∀ n > n , f ( z ) is identified to ∞ X n = a n + n ( z − z ) n + n , with a n + d = P d − ( n ) · a n + d − + . . . + P ( n ) · a n P d ( n ) , n > n a j = lim z → z (cid:16) d j d z j · f (cid:17) ( z ) j ! , j = n , n + , . . . , n + d − . (47)The value of n is deduced by the property P ( n ) · P d ( n ) , , ∀ n > n , that is n = + max { n ∈ Z , P ( n ) · P d ( n ) = } . (48)Notice that n is computed before any cancellation of common factors in (46), which guaranteesthat n does exist in general for any output of FindRE of order d >
1. Indeed, the rewrite rule(45) allows to remark that the differential equation terms with derivative order greater than 1 leadto recurrence equation terms with non-constant polynomial coefficients. The determination of n is crucial to extract parts in the series expansion that are not involved in the summation formula.Much details about the computation of power series extra parts are given in Section 6. Let f ( z ) be a holonomic function. The Taylor expansion of f ( z ) at z is computed as the oneof g ( z ) = f ( z + z ) at 0 if z is a constant, of g ( z ) = f ( − / z ) if z = −∞ , and of g ( z ) = f (1 / z ) if z = ∞ . This algorithm is an immediate use of (47) with the following steps. Algorithm 1
Computing Taylor polynomials of holonomic functions at z ∈ C ∪ {−∞ , ∞} Input:
A holonomic expression f ( z ), a point z , and an integer N . Output:
Taylor polynomial of degree N of f ( z ).1. If z ∈ C , set g ( z ) : = f ( z + z ), else if z = −∞ , set g ( z ) : = f (cid:16) − z (cid:17) , else set g ( z ) : = f (cid:16) z (cid:17) .2. Use FindRE to compute a holonomic recurrence equation satisfied by the Taylor coef-ficients of g ( z ) and write it in the form P d ( n ) · a n + d + P d − ( n ) · a n + d − + · · · + P ( n ) · a n = , n ∈ Z , d ∈ N . (49)3. If d = N with the internal command ofTaylor expansions, say taylor ( f ( z ) , z , N ).4. Compute n = + max { n ∈ Z , P ( n ) · P d ( n ) = } . (50)5. If N n + d − taylor ( f ( z ) , z , N ).6. If N > n + d − T : = taylor ( f ( z ) , z , n + d − . (51)6-1-1 If z ∈ C then compute a j : = coeff ( T , z − z , j ) , j = n , n + , . . . , n + d − , (52)22 lgorithm 1 Computing Taylor polynomials of holonomic functions at z ∈ C ∪ {−∞ , ∞}
6. (resume) where coeff( T , z − z , j ) collects the coefficient of ( z − z ) j in T . The Maximasyntax is adopted.6-1-2 For j = n , n + , . . . , N − d , compute a j + d = P d − ( j ) · a j + d − + . . . + P ( j ) · a j P d ( j ) (53) T = T + a j + d · ( z − z ) j + d (54)6-2-1 Else if | z | = ∞ then compute a j : = coeff ( T , / z , j ) , j = n , n + , . . . , n + d − , (55)6-2-2 For j = n , n + , . . . , N − d , compute a j + d = P d − ( j ) · a j + d − + . . . + P ( j ) · a j P d ( j ) (56) T = T + a j + d · z ! j + d (57)7. Return T . Remark 3. • The relation (53) shows that the coefficients are computed in the same finite number ofoperations. Therefore the complexity is linear. • As we are interested by the asymptotic complexity of this algorithm, there is no issue ofcomparison when the internal command is called in step and . And moreover, this helpsto extract the part of the expansion which cannot be deduced from the recurrence equationused. An example is arcsech( z ) for which Maxima’s command taylor gives the followingexpansion of order at . (%i1) taylor(asech(z),z,0,4); (% o1 ) / T / − log ( z ) + log (2) + · · · − z − · z + · · · • In order to treat certain interesting non-analytic cases like arcsech( z ) in Maxima, insteadof the limit command which can generate errors due to singularities, the internal Maximacommand taylor is used. The initial values are then the coefficients of ( z − z ) j , j = n , . . . , n + d − . We implemented Algorithm 1 as
Taylor with the same syntax as the internal command taylor .Let us present some examples. 23 xample 4. (%i1) Taylor(asech(z),z,0,7); (% o1 ) − log ( z ) − · z − · z − z + log (2) (%i2) Taylor(atan(z),z,inf,7); (% o2 ) − z + · z − · z + · z + π z ) , arctan( z ) , exp( z ) , etc , Maxima seems to use the power series formulaand has very good asymptotic timings. Therefore, for tests we rather use expressions for whichthe internal Maxima command powerseries cannot find the power series formulas. (%i3) Taylor(atan(z)*exp(z),z,0,1000)$ Evaluation took 2.3120 seconds (2.3140 elapsed) using 834.089 MB. (%i4) taylor(atan(z)*exp(z),z,0,1000)$
Evaluation took 34.2350 seconds (34.2220 elapsed) using 5800.949 MB.
The gap between the two computations gets larger as we increase the order.
3. Algorithm mfoldHyper
This section presents a new result in solving holonomic recurrence equations, that is thecomputation of their m -fold hypergeometric term solutions, for positive integers m . Several pro-posals have been given to compute such solutions, among the most recent work in this directionone could cite (Horn et al. (2012)), which is a revisited and improved approach of the one de-scribed in (Petkovˇsek and Salvy (1993)). In the latter, a key step of the proposed algorithm relieson the determination of the linear operator’s right factors of the given holonomic RE. Such a fac-torization is not unique in general because the factors do not commute. In (Horn et al. (2012)),the authors adapted van Hoeij’s approach as explained in (Cluzeau and van Hoeij (2006)) anddefine a concept like the m -Newton polygon for m -fold hypergeometric term solutions of a givenholonomic RE. This approach computes special types of right factors corresponding to m -foldhypergeometric term solutions using the shift operator of order m with the hypothesis that no ra-tional solution exists. With similar approaches, m -fold hypergeometric terms were referred to as m -hypergeometric sequences in (Petkovˇsek and Salvy (1993)), m -interlacings of hypergeometricsequences (see the conclusion of Van Hoeij (1999)), or Liouvillian sequences (see Hendricksand Singer (1999)). Having considered all these developments, we propose to attack the problemfrom another point of view. We will see that contrary to all these previous methods which tryto find solutions in higher dimensions (defined by the solution space) encoded by right factorsof a given RE’s operator, our method only focus on one dimension to deduce a basis of the sub-space of m -fold hypergeometric term solutions. However, before getting in the description ofmfoldHyper, we would like to highlight its importance based in hypergeometric summation.24 .1. Scope of the algorithm Computations of infinite series were connected to holonomic REs by Celine Fasenmyer. Tofind hypergeometric term representations of hypergeometric series, she proposed the followingapproach (see (Koepf, 2014, Chapter 4)).Given a sum P nk = F ( n , k ), write s n = ∞ X k = −∞ F ( n , k ) , (58)and search for polynomials p i , j = p i , j ( n ) , i = , . . . , b , j = , . . . , d with respect to n and do notdepend on k , such that b X i = d X j = p i , j F ( n + j , k + i ) = . (59)If such polynomials are found, then one deduces a holonomic recurrence equation of order atmost d for s n as follows0 = ∞ X k = −∞ b X i = d X j = p i , j F ( n + j , k + i ) = b X i = d X j = p i , j · ∞ X k = −∞ F ( n + j , k + i ) = b X i = d X j = p i , j s n + j = d X j = b X i = p i , j s n + j , where of course we use the fact that p i , j do not depend on k and the advantage of working withbilateral sums: their value is invariant with respect to shifts of the summation variable. In the1940s, Fasenmyer’s method could be used to compute explicit formulas of hypergeometric seriesonly when the obtained holonomic RE was of first order or a two-term recurrence relation. Example 5.
Applied to s n = P nk = k (cid:16) nk (cid:17) , Fasenmyer’s method leads to the recurrence equationns n + − n + s n = , which after use of the initial value s = P k = k (cid:16) k (cid:17) = , yields s n = n n − . For the definite summation case, Gosper (see (Koepf, 2014, Chapter 5)) proposed an algo-rithm which deals with the question of how to find a (forward) anti-difference s k for a given a k ,that is a sequence s k such that a k = ∆ s k = s k + − s k , (60)in the particular case that s k is a hypergeometric term. Thus, once a hypergeometric anti-difference s k of a k is computed, by telescoping definite summation yields n X k = n a k = ( s n + − s n ) + ( s n − s n − ) + · · · + ( s n + − s n ) = s n + − s n , The sum is taken over Z because the summation term vanishes outside a finite set, we say that it has a finite support.
25y an evaluation at the limits of summation. Gosper’s idea is based on the representation a k + a k = p k + p k · q k + r k + , (61)with the property gcd( q k , r k + j ) = , ∀ j ∈ N > (62)that can be algorithmically generated (see (Koepf, 2014, Lemma 5.1 and Algorithm)). Using(61) and (62), one proves that the function f k : = s k + a k + · p k + r k + (63)must be a polynomial for a hypergeometric term anti-difference a k to exist. Thus using (60), (63)and (61) it follows that f k satisfies the inhomogeneous recurrence equation p k = q k + f k + − r k f k − . (64)Gosper gives an upper bound for the degree of f k in terms of p k , q k , and r k which yields a methodfor calculating f k by introducing the appropriate generic polynomial, equating coefficients, andsolving the corresponding linear system so that we finally find s k = r k p k f k − a k . (65)Gosper implemented his algorithm in Maxima as nusum with the same syntax as the Maxima sum command. Example 6. (%i1) nusum(k*k!,k,0,n); solve: dependent equations eliminated: (1) (% o1 ) (1 + n )! − (%i2) nusum(kˆ3,k,0,n); (% o2 ) n · (1 + n ) s n : = ∞ X k = −∞ F ( n , k ) = F ( n , k ) is a hypergeometric term with respect to both n and k with finite support. For thispurpose, one applies Gosper’s algorithm to the expression a k : = F ( n + , k ) − F ( n , k ) (67)26ith respect to the variable k . If successful, this generates G ( n , k ) with a k = F ( n + , k ) − F ( n , k ) = G ( n , k + − G ( n , k ) , (68)and summing over Z yields s n + − s n = ∞ X k = −∞ F ( n + , k ) − F ( n , k ) = ∞ X k = −∞ G ( n , k + − G ( n , k ) = s n is constant, s n = s , and it only remainsto prove that s =
1. In practice, once the function G ( n , k ) is computed one uses the rationalfunction R ( n , k ) : = G ( n , k ) F ( n , k ) (70)called the WZ certificate of F ( n , k ), to establish (66) by proving the rational identity F ( n + , k ) F ( n , k ) − + R ( n , k ) − R ( n , k + F ( n , k + F ( n , k ) = F ( n , k ). Example 7.
For s n : = P nk = F ( n , k ) = P nk = n (cid:16) nk (cid:17) = , the WZ certificate isR ( n , k ) = − k n + − k ) . Therefore the corresponding left hand side of identity (71) isn + n + − k ) − − k n + − k ) + k + n − k )) · n − kk + which trivially yields zero. Although Gosper’s algorithm applies to finite summation, it constitutes a useful tool in dis-covering a method for infinite sums. This is observable in Zeilberger’s algorithm (see (Koepf,2014, Chapter 7)). Zeilberger brings back the computation of a holonomic recurrence equationfor s n : = P ∞ k = −∞ F ( n , k ). The idea is to apply Gosper’s algorithm in the following way: Forsuitable d = , , . . . set a k : = F ( n , k ) + d X j = σ j ( n ) F ( n + j , k ) (72)where σ j is supposed to be a rational function depending on n and not on k . Zeilberger’s mainobservation is that the computation of the polynomial f k defined in (63) yields a linear systemnot only for the unknown coefficients of f k , but also for the rational functions σ j , j = , . . . , d .Thus in a successful case, one obtains an anti-difference G ( n , k ) of a k and rational functions σ j ( n ) , j = , . . . , d such that a k = G ( n , k + − G ( n , k ) = F ( n , k ) + d X j = σ j ( n ) F ( n + j , k ) . (73)27ence, by summation0 = ∞ X k = −∞ G ( n , k + − G ( n , k ) = ∞ X k = −∞ F ( n , k ) + d X j = σ j ( n ) F ( n + j , k ) = s n + d X j = σ j ( n ) s n + j . (74)After multiplication by the common denominator one gets the holonomic recurrence equationsought.Zeilberger’s algorithm gives a much better possibility of computing identities since it com-putes holonomic recurrence equations generally of lowest order (iteration on the order d ) for agiven hypergeometric series. Moreover, this approach is also used to show the coincidence oftwo sums provided their initial values. This constitutes a normal form again, as we have seen inSection 2 . . Example 8.
The sums P nk = (cid:16) nk (cid:17) and P nk = (cid:16) nk (cid:17) (cid:16) kn (cid:17) are proved to coincide by Zeilberger’s algo-rithm as they lead to the same holonomic recurrence equation − ( n + s n + + (7 n + n + s n + + n + s n = , (75) and have the same initial values X k = k ! = X k = k ! k ! = and X k = k ! = X k = k ! k ! = a k is a hypergeometric term with respect to the variable k . When this is not the case,Gosper’s, WZ’s and Zeilberger’s methods cannot be applied. From his algorithm (Koepf, 2014,Algorithm 2.2), Koepf observed that Gosper’s method could miss some results when rational-linear Γ inputs are considered rather than only integer-linear ones. It turns out that this observa-tion constitutes the connection to the general case of m -fold hypergeometric terms (see (Koepf,2014, Chapter 8)). Example 9.
To the Watson’s function F − n b c − n + b + c (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ! , Zeilberger’s algorithm does not directly apply. However, using Koepf’s extended version yieldsthe recurrence equation ( b − c − n − n + a n − ( b − n − c + n + a n + = . (76) Note that the computation of this holonomic RE is made possible after application of (Koepf,2014, Algorithm 8.4) for finding the corresponding m to use. In this example one finds m = . m -fold hypergeometric term solutions of holonomic recurrenceequations could not easily be seen from the known hypergeometric database (see (Koepf, 2014,Chapter 3)). This might be linked to the influence of the combinatorial interpretation oftenpresent in the use of hypergeometric summations in the last century. From this point of view,the m -fold hypergeometric case ( m ∈ N > ) is particularly hidden and Zeilberger’s, Petkovˇsek’sand van Hoeij’s algorithms or their modifications for any hypergeometric situation (multivariatefor example) remain the best approach possible. However, as pointed out in Koepf (1995), thegeneral m -fold hypergeometric case might be the source of a wider family of hypergeometricidentities. This is shown in particular in the computation of power series where some computedholonomic recurrence equations do not have hypergeometric term solutions but only m -fold hy-pergeometric ones ( m ∈ N > ). Note that for the definite summation case the availability ofan algorithm which computes m -fold hypergeometric term solutions of holonomic REs can becombined with Koepf’s extension of Zeilberger’s algorithm to generate new identities. We willnot deal with definite summation, instead as the importance of m -fold hypergeometric terms forhypergeometric series is already shown, we focus on power series computations where the needof such terms is essential for the goal of this paper. There are many limitations of the currentlyused approach to compute power series, we presented some in Section 1; more can be found in(Teguia Tabuguia (2020a)). Let K be a field of characteristic zero. We consider the generic holonomic recurrence equa-tion P d ( n ) a n + d + P d − ( n ) a n + d − + · · · + P ( n ) a n = , (77) P d ( n ) , . . . , P ( n ) ∈ K [ n ] , P d ( n ) · P ( n ) , a n is said to be m -fold hypergeometric, m ∈ N , if there exists afixed rational function r ( n ) ∈ K ( n ) such that r ( n ) = a n + m a n . (78)At first glance, one should remark that m -fold hypergeometric sequences have rational func-tions as the ratio of terms with index difference equal to m . Consequently, if we can find a wayto transform this property to the simple one of hypergeometric sequence then iteratively up tothe order of the given holonomic RE, by van Hoeij’s algorithm we are done.From the characterization (78) one can deduce that for 0 j m − r ( m · n + j ) = a m · ( n + + j a m · n + j . (79)Therefore instead of considering the representation (78) one could rather see an m -fold hy-pergeometric term with m related rational functions as defined in (79). This latter representa-tion is the one used to find a ”simple” formula of an m -fold hypergeometric term as shown in(Teguia Tabuguia (2020b)). Moreover, for fixed m ∈ N and j ∈ N > , if we compute an m -foldhypergeometric term solution of (77) with ratio r ( m · n + j ) for some rational function r , then thisgives the information that there are m − m -fold hypergeometric term solutions of(77). There appears the particularity of our approach. Indeed, r ( m · n + j ) is the term of a sub-sequence of ( r ( n )) n . Unlike previous approaches that try to find r ( n ) directly, we rather compute29he sub-terms r ( m · n + j ) for fixed j ∈ J , m − K so that r ( n ) is constructed as a linear combinationof these sub-terms.Thus for every positive integer m , the computation of an m -fold hypergeometric term solutionof (77) with representation (78) reduces to the computation of an m -fold hypergeometric termsolution of (77) with representation (79) for a fixed j ∈ J , m − K since the other representationscan be similarly computed. When given for a specific value of j , we will say that the solutionis given in an incomplete form. The complete form of the solution is then given with all therepresentations of (79) for all m -fold hypergeometric term solutions. By default in our algorithmwe choose j = K or the index variable subset of Z , holonomic re-currence equations can have m -fold hypergeometric term solutions that cannot be computed over K , but in an extension field instead. This situation occurs with the power series of exp( z ) sin( z )for K = Q . Our implementation HolonomicDE(f,F(z),[destep]) has an optional variable destep whose default value is 1. This number represents the minimum positive difference possible be-tween the derivatives of
F(z) in the holonomic differential equation sought. Our program
FindRE is also adapted for such computations. This particular tool turns out to be important in few cases.Let us examine the situation with exp( z ) sin( z ). (%i1) DE1:HolonomicDE(exp(z)*sin(z),F(z)); (% o1 ) d d z · F ( z ) − · dd z · F ( z ) ! + · F ( z ) = (%i2) DE2:HolonomicDE(exp(z)*sin(z),F(z),2); (% o2 ) d d z · F ( z ) + · F ( z ) = . . (%i3) CompatibleDE(DE1,DE2,F(z)); The two differential equations are compatible(% o3 ) trueLet us now compute the corresponding holonomic REs. (%i4) RE1:FindRE(exp(z)*sin(z),z,a[n]); (% o4 ) (1 + n ) · (2 + n ) · a n + − · (1 + n ) · a n + + · a n = (%i5) RE2: FindRE(exp(z)*sin(z),z,a[n],2); (% o5 ) (1 + n ) · (2 + n ) · (3 + n ) · (4 + n ) · a n + + · a n = RE2 a relationship can be deduced between two of its termswhose index difference is not a multiple of 4. Therefore considering indices in 4 N = { , , , . . . } might be more appropriate. On the other hand, it is trivial that RE2 is a characteristic holonomicrecurrence equation of a 4-fold hypergeometric term. However, this 4-fold hypergeometric termis completely hidden in
RE1 when looking for solutions over Q . Indeed, since the corresponding4-fold symmetric terms are linearly independent over Q , a substitution in the left-hand side of RE1 does not yield zero. Note, however, that
RE1 and
RE2 have hypergeometric term solutions30ver C . We compute these solutions using the first author’s variant of van Hoeij’s algorithm(Teguia Tabuguia (2020b)) implemented in Maxima as HypervanHoeij . (%i6) HypervanHoeij(RE1,a[n],C); (% o6 ) ( (1 − i ) n n ! , (1 + i ) n n ! ) (%i7) HypervanHoeij(RE2,a[n],C); (% o7 ) (cid:16) − ( − · √ (cid:17) n n ! , (cid:16) ( − · √ (cid:17) n n ! , (cid:16) − ( − · √ · i (cid:17) n n ! , (cid:16) ( − · √ · i (cid:17) n n ! The basis of hypergeometric term solutions of
RE1 spans a sub-space of the space of solutionsof
RE2 . The corresponding 4-fold hypergeometric term solutions over Q can be written as alinear combination of the above bases of hypergeometric terms over C . The main thing thatwe point out from this example is that our approach to compute m -fold hypergeometric termsolutions of a given holonomic RE depends on the shifts between indices of the indeterminate( a n ) in that RE and the field considered.Next, we introduce some properties and definitions that clarify this situation and help tocompute m -fold hypergeometric term solutions of (77) in a given field K which we want to bethe smallest algebraic extension field possible of Q in terms of inclusion.The following lemma gives a condition on the order of a given holonomic RE for its m -foldhypergeometric term solutions to be computable over a given field K . An equivalent statementwas given in (Hendricks and Singer, 1999, Theorem 5.1) but with a different perspective asshown by the two proofs. Lemma 10.
Let h n be an m-fold hypergeometric term, m ∈ N . Assume ∀ u ∈ N , u < m , there is no rational function r u ( n ) ∈ K ( n ) : h u + n = r u ( n ) h n . (80) Then there is no holonomic recurrence equation over K of order less than m satisfied by h n .Proof. Let h n be an m -fold hypergeometric term such that h n + m = r ( n ) · h n ⇐⇒ Q m ( n ) · h n + m + Q ( n ) · h n = , (81)where Q m ( n ) , Q ( n ) ∈ K [ n ] and r ( n ) = − Q ( n ) / Q m ( n ) ∈ K ( n ).Suppose that h n satisfies a holonomic recurrence equation of order less than m . Then thereexists an equation of the form P m − a n + m − + P m − a n + m − + · · · + P a n + + P a n = , (82)with polynomials P j = P j ( n ) ∈ K [ n ] , j ∈ J , m − K , and P ( n ) ,
0, satisfied by h n .• If P is the only non-zero polynomial in the equation then h n must be zero, which is acontradiction by definition.• We assume that at least one other polynomial factor in the equation is non-zero. Then h n satisfying (82) yields the following equation after substitution of n by m · nP m − ( mn ) h mn + m − + P m − ( mn ) h mn + m − + · · · + P ( mn ) h mn + + P ( mn ) h mn = . (83)31y assumption (80), we know that ∀ u ∈ N , u < m , h n is not a u -fold hypergeometric term.So the holonomic recurrence equation of lowest order over K satisfied by h n is Q m ( n ) · a n + m + Q ( n ) · a n = , which is a two-term recurrence relation whose subspace of m -fold hypergeometric term ( m is fixed) solutions can be represented by the basis( h mn + m − , h mn + m − , . . . , h mn + , h mn ) . (84)Therefore (83) cannot hold since the left-hand side is a linear combination of linearlyindependent terms with respect to K ( n ), which implies that all the polynomial coefficientsmust be zero. Therefore we get a contradiction.More generally, any shift of a holonomic recurrence equation of order less than m does nothave m -fold hypergeometric term solutions.Remark that checking the hypothesis of Lemma 10 is an important task for the algorithm.Fortunately, this can be done iteratively. Once the field K is fixed, if we have already lookedfor u -fold hypergeometric term solutions for integers u < m , then we can safely proceed to thecomputation of m -fold hypergeometric term solutions knowing that m is less than the order ofthat recurrence equation.Thus, we now know that all the m -fold hypergeometric term solutions of (77) have m d .Furthermore, we can extend this view of m -fold hypergeometric sequences in order to determinewhich type of terms can appear in a holonomic recurrence equation that they satisfy. For thatpurpose, let us first introduce the following definition. Definition 11 ( m -fold holonomic recurrence equation) . A holonomic recurrence equation is saidto be m-fold holonomic, m ∈ N , if it has at least two non-zero polynomial coefficients and thedifference between indices of two appearing terms of the indeterminate sequence in that equationis a multiple of m. Choosing as the trailing term order gives the general formP d ( n ) · a n + md + P d − ( n ) · a n + m ( d − + · · · + P ( n ) · a n + m + P ( n ) · a n = , . (85) so that P d · P , . Assume an m -fold holonomic RE with representation (85) is given. We are going to presenthow to compute a basis of all m -fold hypergeometric term solutions of (85) with representation(79) for j =
0. And by a similar reasoning we will show how to deduce the other bases of m -foldhypergeometric solutions for j ∈ J , m − K .Observe that if we have an m -fold hypergeometric sequence a n starting with a (by shift itis always possible to define the initial term by a ), then using the representation (78) the nextterm which can be computed from a is a m , and afterwards a m . . . , a km , . . . . Thus if we set s n = a mn then all terms computed from s have their indices corresponding to multiples of m for a n . Moreover, since a m · ( n + / a m · n = s n + / s n ∈ K ( n ), s n is a hypergeometric term whoseformula is the same as that of a mn . Therefore we can update (85) accordingly so the algorithm in(Teguia Tabuguia (2020b)) can be applied to compute a basis of all hypergeometric term solution s n of an updated version of (85), which is nothing but the basis of all m -fold hypergeometric termsolutions of (85) with representation (79) for j = m -fold hypergeometric terms is our main idea. As explained, a basis of all m -fold hypergeometric term solutions of the m -fold holonomic RE (85) can be found by van Hoeij’salgorithm provided the following crucial change of variable: m · k = ns k = a m · k . (86)This leads to a 1-fold holonomic RE for s k which has hypergeometric term solutions because s k + s k = a mk + m a mk = r ( mk ) . (87)The resulting RE is P d ( mk ) · s k + d + P d − ( mk ) · s k + ( d − + · · · + P ( mk ) · s k + + P ( mk ) · s k = . (88)In the general case of j ∈ J , m − K , the bases of all m -fold hypergeometric term solutions of(85) with representation (79) are computed by the algorithm in (Teguia Tabuguia (2020b)) afterthe application of the change of variable m · k + j = n , s k = a m · k + j (0 j m − . (89)This is because in (79), the m -fold hypergeometric term indices can always be seen as m · n + j , j ∈ J , m − K .Let us apply this to an example. We consider the two-term recurrence relation of the Taylorcoefficient of exp( z ) sin( z ) which is a 4-fold holonomic RE, so we are going to compute 4-foldhypergeometric term solutions. (%i1) RE:FindRE(exp(z)*sin(z),z,a[n],2); (% o1 ) (1 + n ) · (2 + n ) · (3 + n ) · (4 + n ) · a n + + · a n = (%i2) RE:subst(4*n,n,RE); (% o2 ) (1 + · n ) · (2 + · n ) · (3 + · n ) · (4 + · n ) · a · n + + · a · n = (%i3) RE:subst([a[4*n]=s[n],a[4*n+4]=s[n+1]],RE); (% o3 ) (1 + · n ) · (2 + · n ) · (3 + · n ) · (4 + · n ) · s n + + · s n = (%i4) HypervanHoeij(RE,s[n]); (% o4 ) ( − n · n (cid:16) (cid:17) n · (cid:16) (cid:17) n · n · (2 · n )! This set is a basis of all 4-fold hypergeometric term solutions (90) for j = j = (%i5) RE:FindRE(exp(z)*sin(z),z,a[n],2)$(%i6) RE:subst(4*n+3,n,RE); (% o6 ) (4 + · n ) · (5 + · n ) · (6 + · n ) · (7 + · n ) · a · n + + · a · n + = %i7) RE:subst([a[4*n+3]=s[n],a[4*n+7]=s[n+1]],RE); (% o7 ) (4 + · n ) · (5 + · n ) · (6 + · n ) · (7 + · n ) · s n + + · s n = (%i8) HypervanHoeij(RE,s[n]); (% o8 ) ( − n · n (cid:16) (cid:17) n · (cid:16) (cid:17) n · (cid:0) · n + · n + · n + (cid:1) · n · (2 · n )! On the other hand, note that the m -fold holonomic RE case is the easiest part for the wholealgorithm. Indeed, an arbitrary holonomic recurrence equation is not necessarily m -fold holo-nomic, m ∈ N > . It could have m -fold and m -fold hypergeometric term solutions with positiveintegers m , m . Therefore we should define what to do in the general case.Observe that without the shift that transforms an m -fold holonomic recurrence equation inthe form (85), its general representation is given by P d a n + k + md + P d − a n + k + m ( d − + · · · + P a n + k = , (91)where k ∈ J , m − K .Let us consider the three following 3-fold holonomic REs RE P , · a n + + P , · a n + + P , · a n + = , RE P , · a n + + P , · a n + + P , · a n + + P , · a n + = , RE P , · a n + + P , · a n + + P , · a n + + P , · a n + = . (92)• The difference between an index of the indeterminate ( a n ) in RE RE RE RE a n ) in RE RE RE RE Definition 12.
Let m ∈ N ,RE : P d a n + k + md + P d − a n + k + m ( d − + · · · + P a n + k = , (93) and RE : P d a n + k + md + P d − a n + k + m ( d − + · · · + P a n + k = be two m-fold holonomic recurrence equations. • We say that RE and RE are m-fold distinct holonomic equations if k − k is not divisibleby m. • We say that RE and RE are m-fold equivalent holonomic equations if k − k is divisibleby m.
34n immediate consequence of these definitions is that linear combinations of m -fold equiva-lent holonomic REs always give m -fold holonomic recurrence equations whereas linear combi-nations of m -fold distinct holonomic REs are never m -fold holonomic. For example, let us sum RE RE RE + RE P , · a n + + P , · a n + + (cid:0) P , + P , (cid:1) · a n + + (cid:0) P , + P , (cid:1) · a n + + P , · a n + = . (95)The whole algorithm is based on the following fundamental theorem from which the generalapproach to compute m -fold hypergeometric term solutions of any given holonomic recurrenceequation is deduced. Theorem 13 (Structure of holonomic REs having m -fold hypergeometric term solutions) . Letm ∈ N , K a field of characteristic zero, and h n be an m-fold hypergeometric term which is notu-fold hypergeometric over K for all positive integers u < m. Then h n is a solution of a givenholonomic recurrence equation, if that equation can be written as a linear combination of m-foldholonomic recurrence equations. When this is the case, h n is moreover solution of each of them-fold distinct holonomic recurrence equations of that linear combination.Proof. Let h n be an m -fold hypergeometric term solution of the recurrence equation P d a n + d + P d − a n + d − + · · · + P a n = , d > m , P d · P , . (96)It suffices to show that for any non-zero term P j a n + j in (96), there exists another summand, say P i a n + i , such that m divides j − i . Indeed, by summing m -fold holonomic REs we are sure thatfor each summand appearing on the left-hand side of the sum there must exist another summandwhose index differs from the one of that summand by a multiple of m .We proceed by contradiction. Assume there exists a non-zero term P j a n + j in (96) such thatany other summand P i a n + i , i , j does not verify that m divides j − i . Since h n is a non-zerosolution, we can divide the equation by h n + j and write − P j = d X i = i , j P i · h n + i h n + j . (97)The situation is now two-fold:• For i verifying | i − j | < m , for each corresponding summand P i · h n + i / h n + j on the right-handside of (97), the fact that m does not divide j − i implies that h n + i / h n + j < K ( n ) since byassumption h n is an m -fold hypergeometric term over K that is not u -fold hypergeometricfor all integers u < m . Therefore the whole term P i · h n + i / h n + j < K ( n ).• For i verifying | i − j | > m , for each corresponding summand P i · h n + i / h n + j on the right-handside of (97), we have two possibilities: – either h n + i / h n + j < K ( n ) and we have the same conclusion as in the previous case; – or h n + i / h n + j ∈ K ( n ), but in this case since m does not divide j − i , this implies that h n is not an m -fold hypergeometric term and we get a contradiction.35hus the identity (97) is valid only if all the summands on its right-hand side do not belong to K ( n ). Therefore (97) holds if and only if d X i = i , j P i · h n + i h n + j < K ( n ) , since similarly as in the proof of Lemma 10 we know that the non-rationality is caused by thelinear independency. However, the left-hand side P j ( n ) ∈ K [ n ] ⊂ K ( n ). Hence we obtain acontradiction.Let us now prove the second part of the theorem. Since the multiplication of a holonomicrecurrence equation by a polynomial does not affect the computation of its m -fold hypergeomet-ric term solutions, the linear combination of m -fold holonomic REs can always be considered asa sum of m -fold holonomic REs. Therefore it is enough to show that an m -fold hypergeometricterm solution of a sum of m -fold holonomic recurrence equations is a solution of each of theinvolved m -fold distinct holonomic recurrences.The sum of M m -fold holonomic recurrence equations, M ∈ N , can be written as M X j = RE j ( a n ) = M X j = (cid:16) P d j a n + k j + md j + P d j − a n + k j + m ( d j − + · · · + P j a n + k j (cid:17) = , (98)where k j ∈ J , m − K , and P d j · P j , , j ∈ J , M K .If M =
1, then (98) is an m -fold holonomic recurrence equation and h n is an m -fold hyperge-ometric term solution of it.We assume now that M > m -fold distinct holonomic recur-rence equations in (98). Note that if the M m -fold holonomic REs are m -fold equivalent thenthe situation is similar to the case M = m -fold equivalentholonomic REs is an m -fold holonomic RE.Now suppose that h n is not solution of RE j in (98), j ∈ J , M K , then given that P Mj = RE j ( h n ) =
0, there must be at least one second m -fold holonomic recurrence equation RE j , j ∈ J , M K , m -fold distinct with RE j such that RE j ( h n ) ,
0. Without loss of generality, we consider that RE j is the only second m -fold holonomic RE with these properties. Of course, if RE j ( h n ) , RE j ( h n ) + RE j ( h n ) = RE j ( h n ) ,
0. Thus, we have RE j ( h n ) , RE j ( h n ) , RE j ( h n ) + RE j ( h n ) = . (99)The fact that the m -fold holonomic recurrence equations RE j and RE j are m -fold distinct im-plies that k j − k j is not a multiple of m .Using (99), after substitution of h n in the sum of the equations and division by h n + k j + d j m , wededuce that − P d j = d j − X e j = j P e j h n + k j + e j m h n + k j + d j m + d j X e j = j P e j h n + k j + e j m h n + k j + d j m = S j + S j , j , (100)which is equivalent to − P d j − S j = S j , j . (101)36ll the summands of S j belong to K ( n ) since h n is an m -fold hypergeometric sequence and thecorresponding index differences n + k j + e j m − ( n + k j + d j m ) = m · ( e j − d j )are multiples of m . However, for S j , j the index differences n + k j + e j m − ( n + k j + d j m ) = k j − k j + m · ( e j − d j )are not multiples of m . Therefore by the same argument used in the first part of the proof wededuce that S j , j < K ( n ). Thus (101) holds if and only if − P d j − S j ∈ K ( n ) and S j , j < K ( n ).Therefore we get a contradiction.From this theorem, given m ∈ N , we are now sure to compute a basis of all m -fold hyperge-ometric term solutions of a given holonomic recurrence equation by splitting it into the sum of m -fold distinct holonomic recurrence equations and use van Hoeij’s algorithm or the algorithmin (Teguia Tabuguia (2020b)) to solve these holonomic REs provided the change of variable (89).Note that since we compute m -fold hypergeometric terms as elements of a basis of all m -foldhypergeometric term solutions of holonomic REs, an m -fold hypergeometric term is solutionof two given holonomic REs if it is linearly dependent to an element of the basis of all m -foldhypergeometric term solutions of each of these holonomic REs. Therefore, the solutions soughtare built by all the linearly dependent m -fold hypergeometric term solutions of each involved m -fold distinct holonomic REs. Note that the computation of m -fold hypergeometric term solutionswith representation (79) for j = m -fold holonomic RE P d i a n + k i + md i + P d i − a n + k i + m ( d i − + · · · + P i a n + k i = , (102)is done after writing it in the form (85). Thus (102) is transformed as P d i ( n − k i ) a n + md i + P d i − ( n − k i ) a n + m ( d i − + · · · + P i ( n − k i ) a n = . (103)Let us take as an example the holonomic RE satisfied by the Taylor coefficients ofexp( z ) + cos( z ). (%i1) FindRE(cos(z)+exp(z),z,a[n]); (% o1 ) (1 + n ) · (2 + n ) · (3 + n ) · a n + − (1 + n ) · (2 + n ) · a n + + (1 + n ) · a n + − a n = RE + n ) · (2 + n ) · (3 + n ) · a n + + (1 + n ) · a n + = , and RE − − n ) · (2 + n ) · a n + − a n = . Only RE RE
11 : n · (1 + n ) · (2 + n ) · a n + + n · a n = . From this one easily sees that the given holonomic RE has 2-fold hypergeometric term solutionssince we get two two-term recurrence relations that are linearly dependent: − n · RE = RE . m changes of variable of (89) because aswe explained earlier, once one succeeds in computing a basis of m -fold hypergeometric termsolutions corresponding to the representation (79) for a fixed j ∈ J , m − K , the other ones can becomputed in a similar way. This will be used for power series computations in order to considerall possible linear combinations of hypergeometric type series of type m .This result is a consequence of observing m -fold hypergeometric terms as sequences whoseindices are taken in m · Z + j , j ∈ J , m − K . Commonly this notion is used according to itsdefinition for the set of integers Z which is generally the chosen set of indices. In this case twoterms of a sequence are said to be consecutive if their index difference is 1 or −
1. Such a defi-nition is more useful for hypergeometric terms since it allows van Hoeij’s algorithm to look forhypergeometric term solutions of holonomic REs in such a way that the ratio of two consecutiveterms is a rational function over the considered field. For the m -fold case, however, one ratherneeds to consider m Z as the set of indices so that the computation of m -fold hypergeometric term( m >
2) solutions of holonomic REs is done analogously to the one of hypergeometric terms. Inthis situation one could say that two terms of an m -fold sequence are consecutive if the differenceof their indices is m or − m .To compute the basis of all m -fold hypergeometric term solutions of a given holonomic RE,the algorithm proceeds by iteration up to the order of the RE. Note, however, that more often thenumber of cases to be considered is much smaller than the order of the given RE. For example,the recurrence equation (%i1) RE:FindRE(sin(zˆ3)ˆ3,z,a[n]); (% o1 ) ( n − · ( n − · ( n − · (1 + n ) · a n + + · ( n − · ( n − · a n − + · a n − = m -fold distinct holonomic REs for m < { , , , } .Our algorithm to compute m -fold hypergeometric term solutions of a given holonomic RE,called mfoldHyper, is the following. Algorithm 2 mfoldHyper: m -fold hypergeometric term solutions of holonomic recurrence equa-tion of order d ∈ N Input:
A holonomic recurrence equation P d a n + d + P d − a n + d − + · · · + P a n = , d > m , P d · P , Output:
A basis (incomplete form) of all m -fold hypergeometric term solutions of (104).1. Set H = {} .2. Use the algorithm in (Teguia Tabuguia (2020b)) to find the basis, say H , of all hyper-geometric term solutions of (104). If H , ∅ , then add [1 , H ] to H .3. For 2 m d do:(a) Extract the following m -fold holonomic recurrence equations from (104) and con-struct the system 38 lgorithm 2 mfoldHyper3. (a) P ( n ) · a n + P m ( n ) · a n + m + · · · + P m ·⌊ dm ⌋ ( n ) · a n + m ·⌊ dm ⌋ = P ( n ) · a n + + P m + ( n ) · a n + m + + · · · + P m ·⌊ dm ⌋ + ( n ) · a n + m ·⌊ dm ⌋ + = . . . P m − ( n ) · a n + m − + P m − ( n ) · a n + m − + · · · + P m ·⌊ dm ⌋ + m − ( n ) · a n + m ·⌊ dm ⌋ + m − = , (105)assuming P j ( n ) = j > d .(b) If there exists a holonomic RE with only one non-zero polynomial coefficient in(105), then stop and go back to step 3.(a) for m + m -fold holonomic recurrence equations in (105) so that the order ofthe trailing term equals 0.(d) Apply the change of variable (86) for each m -fold holonomic recurrence equation.(e) Compute a basis of all hypergeometric term solutions s k as defined in (86)for (88) of each resulting holonomic recurrence equation with the algorithm in(Teguia Tabuguia (2020b)).(f) Construct the set H m of hypergeometric terms which are each linearly dependentto one term in each of the m computed bases in step 3.(d).(g) If H m , ∅ then add [ m , H m ] in H .4. Return H .We implemented mfoldHyper in Maxima as mfoldHyper(RE,a[n],[m,j]) , by default [m,j] is an empty list. In that default case each list of m -fold hypergeometric term solutions, say[ m , [ h , m , h , m , . . . ]], contains ”simple formulas” of hypergeometric terms corresponding to j = m -fold hypergeometric term solutions for particular m ∈ N , the algorithm can be called as mfoldHyper(RE,a[n],m,j) for 0 j < m to get the solutionsin their other representations.Let us now apply the algorithm to some examples. We hide the recurrence equations forspace saving purposes. All these computations can be done with our package FPS currentlyavailable as third-party Maxima package on Github. (%i1) RE:FindRE(atan(z)+exp(z),z,a[n])$(%i2) mfoldHyper(RE,a[n]); (% o2 ) "" , ( n ! ) , " , ( ( − n n ) Sometimes computations may involve algebraic extension fields of Q , the syntax is mfold-Hyper(RE,a[n],[K]) for the two possible values K=C or K=Q (default value). To ask for specific m -fold hypergeometric term solutions the syntax is mfoldHyper(RE,a[n],[K,m,j]) . (%i3) RE:FindRE(log(1+z+zˆ2)+cos(z),z,a[n])$(%i4) mfoldHyper(RE,a[n],C); o4 ) , (cid:18) − − √ · i (cid:19) n n , (cid:18) √ · i − (cid:19) n n , ( − i ) n n ! , ( − n n ! , " , ( ( − n (2 · n )! ) where the obtained 2-fold hypergeometric term is the coefficient of the hypergeometric typeseries of cos( z ). (%i5) declare(q1,constant)$(%i6) declare(q2,constant)$(%i7) RE:FindRE(1/((q1-zˆ2)*(q2-zˆ3)),z,a[n])$(%i8) mfoldHyper(RE,a[n],C); (% o8 ) , − s q1 n , q1 ! n , √ · i · (cid:16) q2 (cid:17) − (cid:16) q2 (cid:17) n , q2 ! n , " , ( q1 ! n ) , " , ( q2 ! n ) For these previous examples, the current Maple convert/FormalPowerSeries yields compli-cated power series representations because the above m -fold hypergeometric terms, m >
2, arenot found. Next we compute the power series coefficients of some expressions for which con-vert/FormalPowerSeries does not find representations. (%i9) RE:FindRE(exp(zˆ2)+cos(zˆ2),z,a[n])$(%i10) mfoldHyper(RE,a[n]); (% o10 ) "" , ( n ! ) , " , ( ( − n (2 · n )! ) (%i11) RE:FindRE(cosh(zˆ3)+sin(zˆ2),z,a[n])$(%i12) mfoldHyper(RE,a[n]); (% o12 ) "" , ( n ! , ( − n n ! ) , " , ( ( − n (2 · n )! ) , " , ( · n )! ) (%i13) RE:FindRE(asin(zˆ2)ˆ2+acos(z),z,a[n])$(%i13) mfoldHyper(RE,a[n]); (% o13 ) "" , ( n · n ! n · (2 · n )! ) , " , ( n · n ! n · (2 · n )! ) (%i14) RE:FindRE(sqrt(sqrt(8*zˆ3+1)-1)+sqrt(7+13*zˆ4),z,a[n])$(%i15) mfoldHyper(RE,a[n]); (% o15 ) , (cid:16) (cid:17) n · (cid:16) (cid:17) n · ( − n · n (4 · n − · (2 · n )! , " , ( − − n · ( − n · (2 · n )!(2 · n − · n · n ! ) (%i16) RE:FindRE(sin(zˆ3)ˆ3,z,a[n])$(%i17) mfoldHyper(RE,a[n]); (% o17 ) "" , ( ( − n (2 · n )! , ( − n (2 · n )! ) Let us now use our implementation for the computation of a specific representation of m -foldhypergeometric term solutions. In this case the user has to specify a value for m and j with j ∈ J , m − K . 40 %i18) RE:FindRE(asin(z)ˆ2+log(1+zˆ5),z,a[n])$(%i19) mfoldHyper(RE,a[n],5,0); (% o19 ) ( ( − n · n ) (%i20) mfoldHyper(RE,a[n],5,3); (% o20 ) ( ( − n · (5 · n + ) (%i21) mfoldHyper(RE,a[n],2,1); (% o21 ) ( (2 · n )!(2 · n + · n · n ! ) Eventually, note that the existence of m -fold hypergeometric term solutions of a holonomicrecurrence equation satisfied by the Taylor coefficients of a given expression does not necessarilyguarantee that this expression represents a hypergeometric type function. For example, arctan( z ) · cos( z ) yields a recurrence equations satisfied by the coefficients of cos( z ). (%i22) RE:FindRE(atan(z)*cos(z),z,a[n])$(%i23) mfoldHyper(RE,a[n]); (% o23 ) "" , ( ( − n (2 n ) ! ) However, we know that the coefficient must be different. In the next section, by finding the linearcombination of hypergeometric type power series we will be able to decide using some initialvalues whether a potential coefficient is the correct one.
4. Hypergeometric type power series
The novelty of the results in this section could not be well understood without a precisedefinition of what we consider as hypergeometric type power series.
Definition 14 (Hypergeometric type power series) . Let K be a field of characteristic zero. Foran expansion around z ∈ K , a series s ( z ) is said to be of hypergeometric type if it can be writtenas s ( z ) : = T ( z ) + J X j = s j ( z ) , s j = ∞ X n = n j , a j , n ( z − z ) n / k j (106) where n is the summation variable, T ( z ) ∈ K [ z , / z , ln( z )] , n ∈ Z , J , k j ∈ N , and a j , n is anm j -fold hypergeometric term, m j ∈ N .Thus a hypergeometric type power series is a linear combination of Laurent-Puiseux serieswhose coefficients are m-fold hypergeometric terms. A hypergeometric function is a function thatcan be expanded as a hypergeometric type power series. T is called the Laurent polynomial partof the expansion, and the k j ’s are its Puiseux numbers. z ) in a hypergeometric type expansion is justified by the solution of theunderlying holonomic differential equation (see Kauers and Paule (2011)). The definition in(Koepf (1992)) reduces to the case T = J m , where m is the unique type encounteredin Definition 14. With this new definition, we can define the type of the series (106) as the tuple( m , m , . . . , m J ). Note, however, that we do not compute the coefficients as they appear in (106),but instead for powers of the form z m j · n + i , 0 i < m j which is more suitable for the coefficientscomputed using mfoldHyper.According to the general algorithm described in (Koepf (1992)), we are at the final step ofthe power series computation procedure. Let us first recall what these steps are. For a givenexpression f , we compute its power series in the following way:1. Find a holonomic differential equation for f using the algorithm in Subsection 2.1;2. Convert that holonomic DE into a holonomic recurrence equation satisfied by the powerseries coefficients of f (see (45));3. Solve the obtained holonomic RE which in our case reduces to compute a basis of all the m -fold hypergeometric term solutions of that RE using algorithm mfoldHyper;4. If there are solutions, use initial values to find the linear combination of the resultinghypergeometric type power series that corresponds to the power series expansion of f , ifsuch a linear combination is valid.Regarding Puiseux series, we will give a generalization of an idea in (Gruntz and Koepf,1995, Section 5). We will see that computing Puiseux numbers k j ’s appearing in (106) reducesto finding a number k which can be defined as the Puiseux number of the corresponding hyperge-ometric type series. Once k is found, we use the substitution h ( z ) = f ( z k ) to bring the situation tothe Laurent series one, and finally divide the general power of the indeterminate z in the obtainedpower series representation of h ( z ) by k to get the expansion sought. This is an intermediate stepbetween the second and the third step above.Our goal is to compute a representation of the form f ( z ) = T ( z ) + F ( z ) , (107)where T ( z ) ∈ K [ z , z , log( z )] is a Laurent polynomial in the variable z with coefficients in K [log( z )],and F ( z ) is a linear combination of hypergeometric type series. We mention that T ( z ) is notuniquely determined but its determination will be made more precise by Lemma 16 and Algo-rithm 3.
5. Finding the Puiseux number
Assume we are looking for a representation of the form f ( z ) = T ( z ) + F ( z ) : = T ( z ) + I X i = ∞ X n = s i n z ( m i · n + j i ) / k i (108) Originally the type was used to denote the value of m for an m -fold hypergeometric term coefficient. m i , k i ∈ N , j i ∈ J , m i − K , s i n is an m i -fold hypergeometric term corresponding to j = j i in the representation (79), and T ( z ) is an extra term whose computation will be explained in thenext subsection.It is enough to suppose that F ( z ) in (108) is the sum of two hypergeometric type series oftype m and m since our development works similarly in the general situation. We have F ( z ) : = ∞ X n = s n z ( m · n + j ) / k + ∞ X n = s n z ( m · n + j ) / k , (109)with the same definitions in (108) for I =
2. For simplicity, we also assume that k and k areco-prime. This is to avoid the use of more variables since in particular this assumption impliesthat the least common multiple of k and k is lcm( k , k ) = k · k . Substituting z by z lcm( k , k ) in(109) gives F ( z lcm( k , k ) ) = ∞ X n = s n z ( m · n + j ) · k + ∞ X n = s n z ( m · n + j ) · k (110) = X n ∈ k · ( m · N > + j ) a nk z n + X n ∈ k · ( m · N > + j ) a nk z n , (111)where a i n is obtained from s i n by the change of variable (89), i ∈ { , } .Observe that in (110) the powers of the indeterminate z are integers. In general, the right-hand side of (108) always gives a representation with integer powers when we substitute z by z µ , for any positive multiple µ of lcm( k , k ). Power series with integer powers are dealt with inother sections of this chapter. Thus our aim of determining the positive integers k i , i ∈ J , I K in(108) can be reduced in finding a positive multiple µ of lcm( k , . . . , k I ) so that we can computethe power series of f ( z µ ) and substitute z by z /µ in the obtained representation to get the one of f ( z ).By the general representation (78) of an m -fold hypergeometric term, we know that thereexist rational functions r ( n ) and r ( n ) such that a n + m = r ( n ) · a n and a n + m = r ( n ) · a n , for the coefficients in (111). Therefore we can write a nk + m = r nk ! · a nk and a nk + m = r nk ! · a nk . (112)where nk and nk are not necessarily integers.To compute the holonomic recurrence equation of smallest order for the m -fold and the m -fold hypergeometric terms a nk and a nk , we need to use the smallest integer k such that k · nk ∈ N > and k · nk ∈ N > . Thus k = lcm( k , k ) and the obtained holonomic RE is of coursecompatible with the one computed using FindRE for the input expression F ( z ). From (112),substituting n by lcm( k , k ) · n = k · k · n yields a k · n + m = r ( k · n ) · a k · n and a k · n + m = r ( k · n ) · a k · n . (113)Since a k · n + m and a k · n + m are, respectively, m -fold and m -fold hypergeometric term solu-tions of a holonomic recurrence equation satisfied by the power series coefficients of f ( z ), by43lgorithm mfoldHyper we know how such terms are computed using an algorithm to computethe equivalent hypergeometric terms s i n such that s i n + s i n = a i n + mi a i n = r i ( k i · n ) , i ∈ { , } . By Petkovˇsek’s algorithm we know that ratios of hypergeometric term solutions of holonomicREs are built from monic factors of the corresponding trailing and leading polynomial coef-ficients. This implies in particular that some zeros and poles of r i ( k i · n ) are the roots of theshifted trailing and leading polynomial coefficient of the holonomic recurrence equation com-puted by FindRE for the power series coefficients of f ( z ), i ∈ { , } . Therefore by computingthe least common multiple of all the trailing and leading polynomial coefficient rational rootdenominators of that RE we must obtain a multiple of lcm( k , k ). Example 15.
We consider the expression f ( z ) = exp( z / ) + sin( √ z ) . (%i1) f:exp(zˆ(3/4)) + sin(sqrt(z))$(%i2) RE:FindRE(f,z,a[n])$ We collect the coefficients with our Maxima function
REcoeff (%i3) CoeffsRE: REcoeff(RE,a[n])$
The corresponding leading polynomial coefficient is (%i4) last(CoeffsRE); (% o4 ) − · (8 + n ) · (9 + n ) · (15 + · n ) · (17 + · n ) · (27 + · n ) · (33 + · n )and the trailing one is (%i5) first(CoeffsRE); (% o5 ) − , , , =
4. Indeed the factors (8 + n ) and(9 + n ) have both denominator roots equal to 1, (15 + · n ) and (17 + · n ) have both denominatorroots equal to 2, and (27 + · n ) and (33 + · n ) have both denominator roots equal to 4. Aftersubstitution the new holonomic RE is free of Puiseux numbers. (%i6) RE:FindRE(subst(zˆ4,z,f),z,a[n])$ with leading term (%i7) CoeffsRE:REcoeff(RE,a[n])$(%i8) last(CoeffsRE); (% o8 ) − n +
14) ( n +
15) ( n +
16) ( n + Integer shift used in Petkovˇsek’s algorithm, see also Lemma 16 .1. Computing starting points We have observed that Maple’s current convert/FormalPowerSeries command wrongly rep-resents the power series of arctan( z ) + exp( z ) due to the constant term missing. In this section weshow how to avoid such a situation by explaining how to deduce exact starting points of hyper-geometric type series from the holonomic recurrence equations of their coefficients. By tryingto compute the representations of many examples of sums of polynomials and hypergeometricseries in Maple, one realizes that such a computation is not well-managed in the implementedalgorithm. A simple example is the following.Maple’s FPS gives > FPS(z + zˆ2 * exp(z),z,n);
FPS ( z + z e z , z , n )whereas our Maxima FPS implementation yields correctly (%i1) FPS(z+zˆ2*exp(z),z,n); (% o1 ) ∞ X n = z + n n ! + z which is the sum of the hypergeometric series of z · exp( z ) whose starting point is n = z . The complete information is detected from the corresponding holonomic re-currence equation. Let us now explain how this can be done.Again, we consider the general representation (assuming Puiseux numbers all equal to 1) f ( z ) : = T ( z ) + F ( z ) (114)where F ( z ) is a sum of hypergeometric type series and T ( z ) ∈ K [ z , z , log( z )] is an extra term tobe determined while computing the starting point for F ( z ). Note that T ( z ) can be given explicitlyin the input expression, but also implicitly like for the expressions arcsech( z ), arccosh( z ) andexp( z ) + log(1 + z ).First, we focus on the case where T ( z ) is a Laurent polynomial in K [ z , z ]. For this purposewe need to understand what it means for a Laurent polynomial that its coefficients are solutionof a holonomic recurrence equation. Let us compute the holonomic RE for an unknown Laurentpolynomial and figure out some properties of its coefficients from that RE. (%i1) FindRE(randompoly(z),z,a[n]); (% o1 ) − + n ) a n −
10 (3 + n ) a n − + n − a n − − n − a n − − ( n − a n − +
26 ( n − a n − = . (115)Our code randompoly is used to generate an arbitrary Laurent polynomial. Since all polyno-mials are rational functions, HolonomicDE always computes a holonomic differential equation offirst order for a given polynomial. Thus we can find a general representation of their holonomicrecurrence equations. T ( z ) : = N X i = M c i z i = X i ∈ Z c i z i ∈ K ( z ) , (116)45or M , N ∈ Z , M N where c i = i ∈ Z \ J M , N K , and c − M · c N ,
0, then the differentialequation found is N X i = − M c i z i · F ′ ( z ) − N X i = − M c i iz i − · F ( z ) = . (117)Therefore using the rewrite rule (45) we obtain the recurrence equation N X i = M c i ( n + − i ) · a n + − i − N X i = M c i i · a n − ( i − = N X i = M c i ( n + − i ) · a n + − i = . Hence the holonomic RE found by
FindRE of a Laurent polynomial with representation (116) isgiven by N X i = M c i ( n + − i ) · a n + − i = , (118)or equivalently N X i = M c i ( n + N − i ) · a n + N − i = , (119)after substitution of n by n + N for normalization.Thus, without even using initial values a polynomial whose coefficients satisfy the holonomicRE (115) can easily be found by equating the terms of (118) and (115) to find the unknowncoefficients c i , using FindRE to compute a holonomic RE for the resulting polynomial and checkwhether the REs are identical. We obtain the Laurent polynomial − · z + z − z + z + z + . (120)Of course this is not enough because there might be other solutions. And moreover, when theinput expression is of the form (114), the situation is more complicated since F ( z ) needs initialvalues in order to be computed. Therefore, we have to find the maximum degree N ∈ Z of T ( z )so that F ( z ) starts at N + T ( z ) is computed by a generalized Taylor expansion of order N of f ( z ).Observe for each non-zero coefficient c i , i ∈ J M , N K of T ( z ), that 2 i − N is the root of onepolynomial coefficient in (119). In particular, N is the trailing polynomial coefficient root and M is the root of the leading polynomial coefficient shifted by N − M . These two properties ofthe degrees of a potential Laurent polynomial whose coefficients satisfy a holonomic recurrenceequation is preserved in the general case. This is stated by the following lemma. Lemma 16.
Let K be a field of characteristic zero, N , M ∈ Z , N > M, T ( z ) ∈ K [ z , z ] be aLaurent polynomial of degree N and lowest non-zero monomial degree M. The coefficients ofT ( z ) satisfy the holonomic recurrence equationP d a n + d + P d − a n + d − + . . . + P a n = , (121) d ∈ N , P j ∈ K [ n ] , j ∈ J , d K , P d · P , , if N is a root of P and M is a root of P d ( n − d ) . roof. Suppose that the coefficients of T ( z ) satisfy (121). Since T ( z ) has finitely many non-zerocoefficients we can write T ( z ) = X i ∈ Z c n z n , where c n = n ∈ Z \ J M , N K . Saying that the coefficients of T ( z ) satisfy (121) is equivalentto say that the sequence ( c n ) n ∈ Z is a sequence solution of (121). Given that (121) is valid forall integers, observe that substituting a n by c n in (121) for sufficiently large positive or negativeintegers all the summands on the left-hand side of (121) vanish.Furthermore, we can make a substitution such that either the trailing or the leading term doesnot necessarily give zero. Indeed, since c n = n ∈ Z \ J M , N K , substituting a n by c n in (121)for n = N yields P ( N ) c N = , and therefore using the assumption c N , P ( N ) =
0. Similarly, substituting a n by c n in (121) for n = M − d gives P d ( M − d ) c M = , and therefore as c M , P d ( M − d ) = Remark 17.
Note that generally when T ( z ) = and F ( z ) starts at , N = M = and isnot necessarily a root of the trailing polynomial coefficient. This might be interpreted from thefact that the zero function is always solution of any holonomic RE and moreover takes at .Therefore as T ( z ) does not play a disturbing role, we rather say in this case that it does not exist.This is the case with (%i2) FindRE(exp(z),z,a[n]); (% o2 ) (1 + n ) · a n + − a n = whose trailing polynomial coefficient does not have any root. T ( z ) in this case is or is said tonot exist. However, note that in this example M = is a root of the leading polynomial coefficientwhich represents the starting point of the series expansion of exp( z ) . In general, the computationof M is always possible from all the REs computed by FindRE, and represents moreover thestarting point for the series expansion of the given f ( z ) . Indeed, the fact that FindRE doesnot cancel the common factors after application of the rewrite rule (45) is essential for ourcomputations of starting points. These factors contain necessary information to determine thefirst non-zero coefficient of the series expansion sought. Let for example (%i3) FindRE(z/(1-z),z,a[n]); (% o3 ) ( n − · a n − − ( n − · a n = for which the cancellation of the common factor ( n − (or n after normalization) would hide thestarting point N + = (N = is the root of the trailing polynomial coefficient). Note that using this lemma we can now confirm that any Laurent polynomial whose sequenceof coefficients satisfies the holonomic RE (115) is a constant multiple of the polynomial (120).Indeed, the leading and the trailing polynomial coefficients of (115) have only one integer rooteach which are the degree bounds of (120). Algorithmically, we proceed as follows.47 lgorithm 3
Computing T ( z ) and the starting point of F ( z ) for f = T ( z ) + F ( z ) as in (114) Input:
An expression f whose series coefficients satisfy the holonomic recurrence equation P d a n + d + P d − a n + d − + . . . + P a n = , (122) d ∈ N , P j ∈ K [ n ] , j ∈ J , d K , P d · P , , Output: T ( z ) and a starting point N for F ( z ) for the representation (114) of f .1. Compute the minimum integer roots M of P d ( n − d ) and the maximum integer root N of P ( n ).2. If N does not exist then set T ( z ) : = N : = M .3. If N does exist then set T ( z ) : = T aylor ( f ( z ) , z , , N ) and set N : = N + T ( z ) , N ].Note that Lemma 16 extends to Laurent polynomials in K [log( z )][ z , z ] as connected to gener-alized Taylor expansions or power series in (Kauers and Paule, 2011, Section 7.3). Our Maximapackage has the code LPolyPart(f,z) that implements Algorithm 3. This can be used as follows.
Example 18. (%i1) LPolyPart(asech(z),z); (% o1 ) [log (2) − log ( z ) , (%i2) LPolyPart(exp(z)+log(1+z),z); (% o2 ) [1 , (%i3) LPolyPart(sin(z)/zˆ5,z); (% o3 ) [0 , − z )). (%i4) LPolyPart(cos(4*acos(z)),z); (% o4 ) [8 · z − · z + , F ( z ) is 5. This example leadsto a two-term holonomic recurrence equation. We will see in the next subsection that the linearcombination of the corresponding hypergeometric type series yields 0 so that one finally gets theknown result cos(4 arccos( z )) = T ( z ) = · z − · z + One may ask why we need such an algorithm since it is already generalized by using mfold-Hyper. This could be answered by the following example whose corresponding two-term holo-nomic recurrence relation has hypergeometric and 2-fold hypergeometric term solutions over Q . 48 %i1) RE:FindRE(cosh(z),z,a[n]); (% o1 ) (1 + n ) · (2 + n ) · a n + − a n = (%i2) mfoldHyper(RE,a[n]); (% o2 ) "" , ( n ! , ( − n n ! ) , " , ( · n )! ) (%i3) mfoldHyper(RE,a[n],2,1); (% o3 ) 1(2 · n + · (2 · n )!Thus using mfoldHyper imposes to decide between the representationscosh( z ) : = ∞ X n = + ( − n n ! z n , (123)and cosh( z ) = ∞ X n = z · n (2 · n )! , (124)which are both correct. We will see in the next subsection that this is decided by the linearsystem to be solved; arbitrary constants appearing in the solution are set to zero to ease thedecision and this may lead to a little more complicated representation. Therefore we propose torevisit the algorithm as described in (Koepf (1992)) with our formalism of starting points andPuiseux numbers of series expansions.Let us sketch the two-term holonomic RE algorithm for our introductory example. The re-currence equation found for cosh( z ) is(1 + n ) · (2 + n ) · a n + − a n = . We immediately get the symmetry number m =
2. Therefore the corresponding 2-fold symmetricratios are a n + a n = · n + · (2 · n +
2) and a n + + a n + = · n + · (2 · n + . (125)We get the coefficients 1(2 · n )! and 1(1 + · n )! , respectively. We now write I ( z ) : = α · X n = z · n (2 · n )! + α · X n = z · n + (2 · n + , and use 2 initial values n = = · n = = · + α and α . We have 49 %i4) taylor(cosh(z),z,0,1); (% o4 ) / T / + . . . therefore α + α · z =
1, hence α = α =
0. And finally we obtain the power seriesrepresentation cosh( z ) = ∞ X n = z · n (2 · n )! , (126)as expected.For cos(4 · arccos( z )), FindRE gives the recurrence equation (%i5) RE:FindRE(cos(4*acos(z)),z,a[n]); (% o5 ) ( n − · (4 + n ) · a n − (1 + n ) · (2 + n ) · a n + = − (1 + n ) · (9 + n )( n + · ( n + . Therefore the corresponding 2-fold symmetric ratios are r = − (1 + · n ) · (9 + · n )(2 · n + · (2 · n +
7) and r = − (2 + · n ) · (10 + · n )(2 · n + · (2 · n +
8) (127)which lead to the coefficients2 · (7 + · n ) · ( − n · (2 · n )!7 · ( n + · ( n + · n · n ! and 15 · (1 + n ) · (2 + n ) · (4 + n ) · ( − n · n · n ! (5 + · n )! , respectively. We have to use 2 initial values corresponding to n = · + n = · + + (%i6) I:alpha[0]*subst(0,n,h0)*zˆ5 + alpha[1]*subst(0,n,h1)*zˆ6+ ratdisrep(taylor(cos(4*acos(z)),z,0,4)); (% o6 ) α · z + α · z + · z − · z + α and α . Remember that taylor (cos(4 arccos( z )) , z , ,
4) is T ( z ) forthe representation (114) of cos(4 arccos( z )). To find the values of α and α , we just have to solvethe trivial identity (%i7) I-Taylor(cos(4*acos(z)),z,0,6)=0; (% o7 ) α · z + α · z = α = α = · arccos( z )) = − · z + · z In a nutshell the algorithm is presented as follows.50 lgorithm 4
Power series representations of hypergeometric type functions whose
FindRE com-putes two-term holonomic REs for their coefficients.
Input:
A holonomic expression f and a holonomic RE Q ( n ) a n + m − P ( n ) a n = , (128) P , Q ∈ K [ n ] , P · Q ,
0, computed by
FindRE for the series coefficient of f . Output:
The power series representation of f at 0.1. Use Algorithm 3 to compute the corresponding T ( z ) and a starting point N for therepresentation (114) of f .2. Set r ( n ) = P ( n + N ) / Q ( n + N ) .
3. Compute the m symmetric ratios r j ( n ) : = r ( m · n + j ) , j = , . . . , m − . (129)4. Compute formulas (see ”simple formulas” in (Teguia Tabuguia (2020b)) or hypergeo-metric formula in (Koepf, 2014, Chapter 2)) for h j ( n ) : = Q n − i = r j ( i ), for j = , . . . , m − .
5. Set I ( z ) : = T ( z ) + P m − j = α j · h j (0) · z j + N , where α j , j = , . . . , m − α j by equating the coefficient of the Laurent polynomial T aylor ( f ( z ) , z , , N + m − − I ( z ) (130)to zero.7. Return T ( z ) + P m − j = α ′ j · P ∞ n = h j ( n ) · z m · n + j + N , where α ′ j is the value found for α j .Most results of the category of Algorithm 4 are well handled by the current Maple con-vert/FormalPowerSeries . However, since Puiseux numbers and starting points are not sym-bolically computed in this previous implementation, we give two examples missed by con-vert/FormalPowerSeries but accessible to our implementation. Many more can be constructed. > FPS[FPS](1/(sqrt(1-4*z))*((1-sqrt(1-4*z))/2*z)ˆ2,z,n) ∞ X n = (2 n + n +
2) ( n + z n + (( n + > FPS(1/(z+zˆ2),z,n) ∞ X n = ( − n z n − Although Algorithm 4 can be applied to summands of linear combinations of hypergeometrictype functions to get good results in certain cases, we mention that this is not equivalent with thedirect approach to be described next. Emphasis on that fact can be found in (Teguia Tabuguiaand Koepf (2020)). 51 .3. Hypergeometric type series with arbitrary holonomic REs
As previously we consider an expression f ( z ) related to a hypergeometric type function andwe want to compute the representation f ( z ) : = T ( z ) + F ( z ) , (131)where T ( z ) is generally a Laurent polynomial in K [ z , z , log( z )], and F ( z ) is a hypergeometrictype series. We have already shown how to compute T ( z ) and a starting point for F ( z ).We have to avoid negative arguments for the evaluation of m -fold hypergeometric termswhich are supposed to start at least at 0 according to the algorithm in (Teguia Tabuguia (2020b)).When there are many different types involved in F ( z ) as pointed out with sin( z ) / z + exp( z ), theexact starting point is in fact the minimal one among those of the hypergeometric type seriesin F ( z ). We only have to make sure that this value is positive in order to avoid inappropriateoperations like a division by 0 or factorials of negative integers.In practice, we use N = max { , N } > , (132)from which all necessary m -fold hypergeometric terms can be evaluated for initial conditions.Once the linear combination sought is found, if possible we subtract terms from the correspond-ing F ( z ) for the indices n ∈ J N , N K . Note that T ( z ) should also be modified accordingly. Weset T ( z ) : = T aylor ( f ( z ) , z , , N − . (133)By these measures, the final representation gives a normal form after shifting the power of theindeterminate by the first non-zero term index in each hypergeometric type series.Let us now find a representation (131) of f ( z ) for a computed T ( z ) and a starting point N knowing that we can make further computations to subtract terms in T ( z ) that can be deducedfrom F ( z ) and vice versa.To ease the understanding of the general case, we first describe the details for the situationwhere H = [[2 , { h n } ] , [3 , { h n } ]] , (134)represent the obtained basis of m -fold hypergeometric term solutions of the holonomic RE givenby FindRE for the series coefficients of f ( z ). Let us also assume for simplicity that N = T ( z ) =
0, then the general form for the corresponding F ( z ) can be written as F ( z ) = α , · ∞ X n = h n z n + α , · ∞ X n = h n + z n + + α , · ∞ X n = h n z n + α , · ∞ X n = h n + z n + + α , · ∞ X n = h n + z n + , (135) α , i , α , j ∈ K , i = , , j = , ,
2. Hence we have five unknowns to determine. Observe thatcomputing a series expansion of order 4 of f ( z ) might not be enough. Indeed, the Taylor expan-sion of order 4 would give five linear equations for the unknown constants but it turns out thatthe obtained linear system is not sufficient to determine these. Assume T aylor ( f ( z ) , z , , = t + t z + t z + t z + t z , (136)52hen equating the coefficients with their corresponding terms in (135) yields the linear system α , · h n (0) + α , · h n (0) = t α , · h n + (0) + α , · h n + (0) = t α , · h n (1) + α , · h n + (0) = t α , · h n + (1) + α , · h n (1) = t α , · h n + (2) + α , · h n + (1) = t (137)from which a value for α , cannot be deduced because it appears in three equations with threeother different unknown constants.What we need is to use the series expansion of order p ∈ N of f ( z ) in such a way that thereexists q ∈ N , q p , so that there are at least q linear equations with q unknowns each in theresulting linear system. The minimal value of such a p in this particular example is 2 · x = · x where x and x are the minimal positive integers verifying 2 · x = · x , hence x = , x = p = = lcm(2 , α , and α , and this allows to find their values and deduce those of the other constants. If moreoverthere were a hypergeometric term in (134), then 6 linear equations could not be enough. In thiscase the minimal value for p would be 2 · =
12 in order to have at least three equations for α , ,α , and the unknown constant related to the hypergeometric term.We now move to the general case.Let H : = (cid:20) (cid:2) , (cid:8) h n , , . . . , h n , l (cid:9)(cid:3) , h m , n h m n , , . . . , h m n , l m oi , . . . , h m µ , n h m µ n , , . . . , h m µ n , l m µ oi (cid:21) = h(cid:2) , S , (cid:3) , (cid:2) m , S m , (cid:3) , . . . , h m µ , S m µ , ii (138)for integers 1 < m < · · · < m µ be the non-empty generator of all m -fold hypergeometricterm solutions of a holonomic recurrence equation satisfied by the series coefficients of f ( z ). m µ is the maximum symmetry number, l m is the number of m -fold hypergeometric terms in H m ∈ { , m , . . . , m µ } . The representation (131) for f ( z ) is computed as follows. Algorithm 5
Computing hypergeometric type series
Input: f ( z ), the recurrence equation, say RE computed by FindRE , the span of all m -fold hyper-geometric term solutions of RE , say H , computed by mfoldHyper, T ( z ) and N computed byAlgorithm 3. Output:
The representation (131) of f .1. Find the other m -fold symmetric terms associated to each m -fold hypergeometric termin H for m ∈ { m , . . . , m µ } . For that purpose one calls Algorithm 2 as mfoldHy-per(RE,a[n],m,j) for j = , . . . , m − , m ∈ { m , . . . , m µ } . This allows to build thesets S m : = (cid:8) S m , , S m , . . . , S m , m − (cid:9) , (139)for m ∈ { , m , . . . , m µ } , where S m , j : = n h mn + j , , h mn + j , , . . . , h mn + j , l m o , j m − . (140)53 lgorithm 5 Computing hypergeometric type series2. Set N = max { , N } and T ( z ) : = T aylor ( f ( z ) , z , , N − i m , j = l N − jm m for j = , . . . , m − m ∈ { m , . . . , m µ } .4. Set N = N + X m ∈{ , m ,..., m µ } l m − · lcm(1 , m , . . . , m µ ) + m µ − p m , j = j N− jm k , j = , . . . , m − m ∈ { m , . . . , m µ } .6. Let α m , j , k ∈ K , m ∈ { , m , . . . , m µ } , j = , . . . , m − k = , . . . , l m be some unknownconstants and define I ( z ) : = X m ∈{ , m ,..., m µ } m − X j = l m X k = α m , j , k p m , j X n = i m , j h mn + j , k z mn + j . (142)7. Solve the linear system resulting from the equation I ( z ) + T ( z ) − T aylor ( f ( z ) , z , , N ) = , (143)for the unknown (cid:16) α m , j , k (cid:17) Tm ∈{ , m ,..., m µ } , j m − , k l m ∈ K P m ∈{ , m ,..., m µ } l m · m .8. If there is no solution then stop and return FALSE. No linear combination exists in thiscase.9. If there is a solution then set all parameters of dependency to 0 (if there are some).This gives the choice of the linear combination. We denote by α ′ m , j , k the resulting valuefound for α m , j , k , m ∈ { , m , . . . , m µ } , j = , . . . , m − k = , . . . , l m .10. For each S m , m ∈ { , m , . . . , m µ } construct the term S ′ m : = X S m , j ∈ S m X h mn + j , k ∈ S m , j α ′ m , j , k h mn + j , k z mn + j − i m , j (144): = m − X j = l m X k = α ′ m , j , k h mn + j , k z mn + j − i m , j (145)11. For each S ′ m , m ∈ { , m , . . . , m µ } , make evaluations for n ∈ J N , N K to subtract termsin T ( z ) that can be computed from S ′ m and shift the initial index i m , j accordingly. (Thisstep could also be done before step 10 to get more suitable values for starting points).12. Return T ( z ) + P m ∈{ , m ,..., m µ } P ∞ n = S ′ m .The correctness of this algorithm depends on whether the solution of the linear system in step7 has enough equations to determine the possible coefficients of the linear combination sought.Indeed, we saw for (134) that we need a linear system for which each unknown has enoughequations to be determined. This is established by the following lemma.54 emma 19. In Algorithm , N given in (141) is a valid integer for which the series expansionof order N of f ( z ) allows to determine the linear combination sought.Proof. The computation is similar for any integer N , therefore we assume that N =
0. Thenumber of unknowns in each equation is q = P m ∈{ , m ,..., m µ } l m . The aim is to find N such that T aylor ( f ( z ) , z , , N ) in Algorithm 7 step 7 yields a linear system with at least q equations with q unknowns each. Of course, the minimal value of N is an integer that verifies N = m · x = m · x = · · · = m µ · x µ , for some positive integers x , x , . . . , x µ , since we have to find q equations that correspond to the q first coincidences of z m · n , z m · n , . . . , z m µ · n . The second coincidence is reached at the expansion of order lcm( m , . . . , m µ ), therefore by in-duction we deduce that for any positive integer p , the p th coincidence is reached at the expansionof order ( p − · lcm( m , . . . , m µ ). Hence we finally get N = ( q − · lcm( m , . . . , m µ ) = X m ∈{ , m ,..., m µ } l m − · lcm( m , . . . , · m µ ) + m µ − m µ − z m · n + j , z m · n + j , . . . , z m µ + j · n , with j ∈ J , m K . Remark 20. In (141) we use lcm(1 , m , . . . , m µ ) because it allows to recover the order ( l − · when there are only hypergeometric terms ( µ = ). As already used many times, the command
FPS(f(z),z,n,[z 0]) of our Maxima FPS packagecomputes the power series representation of f(z) at the point of expansion z ∈ C (if given or 0otherwise) with the index variable n by combining FindRE , mfoldHyper , Puiseuxnbrfun and ourimplementations of Algorithms 3 and 4 if the computed holonomic RE is a two-term holonomicRE or Algorithm 5, if f(z) leads to a different type of holonomic RE.
Example 21. (%i1) FPS(sin(z)ˆ2+cos(z)ˆ3,z,n); (% o1 ) ∞ X n = − ( − ( − n − · ( − n + · ( − n · n ) · z · n · (2 · n )! + (%i2) FPS(1/(q1-zˆ2)/(q1-zˆ3),z,n); (% o2 ) ∞ X n = q1 − n − z n + q1 − + ∞ X n = q1 − n − z n + q1 − + ∞ X n = − q1 − n − z n + q1 − + ∞ X n = q1 − n − z n q1 − + ∞ X n = − q1 − n − z n q1 − %i3) FPS(sin(zˆ(1/3))+ cos(zˆ(1/2)),z,n); (% o3 ) ∞ X n = ( − n · z + · n (2 · n + · (2 · n )! + ∞ X n = ( − n · z n (2 · n )! (%i4) FPS(acos(zˆ(1/2))+exp(zˆ2),z,n); (% o4 ) ∞ X n = − (2 · n )! · z + · n (2 · n + · n · n ! + ∞ X n = z · n n ! + + π − (%i5) FPS(log(1+sqrt(z)+z+zˆ(3/2)),z,n); (% o5 ) ∞ X n = ( − n z n + n + + ∞ X n = ( − n z n + n + (%i6) FPS(sin(2*z)+cos(z),z,n,%pi/2); (% o6 ) − ∞ X n = ( − n · (1 + · n ) · (cid:16) z − π (cid:17) + · n (2 · n + · (2 · n )! (%i7) FPS(exp(z)+log(1+z),z,n,%e); (% o7 ) ∞ X n = (cid:16) e e + e + e + e e · n + e + e · n + ( − n · (1 + n )!(1 + e ) n (cid:17) · ( z − e ) + n ( e + · ( n + · (1 + n )! + log ( e + + e e The latter series is wrongly represented by Maple’s current convert/FormalPowerSeries
Maple’scommand due to the missing term log ( e + + e e . Algorithm 19 ends our direct method for hypergeometric type series given in Definition 14.In the next section we present a general approach for expressions that are not holonomic.
6. Non-holonomic power series
We described in sub-subsection 2.3.1 that holonomic recurrence equations together with suf-ficient initial values can be used to identify holonomic functions. In this section we give anapproach to represent the power series of expressions that satisfy homogeneous quadratic differ-ential equations. For a given expression f , the procedure follows the following steps:1. compute a quadratic differential equation satisfied by f ;2. use the Cauchy product rule to convert that quadratic differential equation to a non-holonomicrecurrence equation satisfied by the series coefficients of f ;3. use the obtained recurrence equation to define a recursive formula for the power seriescoefficients of f .Gathering all these steps together we are able to define a normal form representation of non-holonomic power series. 56 .1. Computing quadratic differential equations Let f be an expression, our algorithm in subsection 2.1 searches for a holonomic differentialequation for f by iteration on the order of derivatives of f . Depending on the algebraic simplifi-cation of ratios of these derivatives, the algorithm finds a holonomic DE of lowest order satisfiedby f . Our idea in computing homogeneous quadratic differential equation is to define a ’natural’ordering between products of derivatives of f so that by iteration on this ordering a quadraticdifferential equation satisfied by f is sought.What we refer to as quadratic differential equation is an algebraic ordinary differential equa-tion (AODE) of degree 2. This is a DE with at most one product of two derivatives in at least oneof its summands (differential monomials). We recall that any product of derivatives increasesthe degree of a differential monomial by 1 and therefore the degree of an AODE is the great-est number of products of derivatives among its differential monomials plus 1 (see (Eremenko,1982, Section 6)). In particular, linear differential equations are of degree 1. However, for a gooddefinition of the ordering we are looking for, we first study the general AODE case.Let f ( z ) be a given differentiable function. By observing the product rule when applying thederivative operator ddz (a) ddz (cid:16) f ( z ) (cid:17) = · f ( z ) · ddz f ( z ) , (b) ddz f ( z ) ! = − ddz f ( z ) f ( z ) , (146)one can assume that the maximum degree for f ( z ) in an AODE involving f ( z ) and ddz f ( z ) is 2.Now, differentiating the right-hand sides of ( a ) and ( b ) in (146) yields ( c ) and ( d ), respec-tively, as below.(c) 2 · f ( z ) · ddz f ( z ) + · ddz f ( z ) ! , (d) 2 · ddz f ( z ) ! f ( z ) − d dz f ( z ) f ( z ) . (147)Thus one can assume that the maximum degrees for f ( z ) and ddz f ( z ) in an AODE involving f ( z ), ddz f ( z ), and d dz f ( z ) are, respectively, 3 and 2.Using this process recursively, we can state that• f ( z ) , ddz f ( z ) ! and d dz f ( z ) ! gives the maximum degrees for an AODE of order f ( z ) , ddz f ( z ) ! , d dz f ( z ) ! and d dz f ( z ) ! gives the maximum degrees for an AODE oforder 4;• . . . As in the linear case, the order is taken as the greatest derivative order. f as follows(1) 1 , (2) f , (3) f , (4) f ′ , (5) f ′ f , (6) f ′ , (7) f ′′ , (8) f ′′ f , (9) f ′′ f ′ , (10) f ′′ , (11) f ′′′ , (12) f ′′′ f , (13) f ′′′ f ′ , (14) f ′′′ f ′′ , (15) f ′′′ ,. . . (148)We define d − dz = f ( − = , and d dz = f (0) = f . (149)Observe in (148) that for every derivative of order n , n ∈ N > , we compute the product of f ( n ) and all the derivatives of order less than or equal to n before computing the next derivative.We are going to define a derivative operator, say δ k , z f , k ∈ N so that in (148) the numbers inparenthesis represent the derivative orders. This operator computes the product of two derivativesof f according to the ordering given in (148).Looking at (148) as an infinite lower triangular matrix reduces the definition of δ , z to theone of a bijective map ν between positive integers and the corresponding subspace of N × N :( i , j ) i , j ∈ N , i j . This can be done by counting the couple ( i , j ) in (148) from up to down, fromthe left to the right. We obtain ν ( k ) = ( i , j ) = ( l , l ) if N = k ( l + , k − N ) oherwise , where l = r k + − , and N = l ( l + . (150)It remains to define a correspondence between the couple ( i , j ) = ν ( k ) , k ∈ N and the quadraticproducts in (148). This is straightforward since we have defined (149). We get δ k , z ( f ) = d i − dz j − f · d j − dz i − , where ( i , j ) = ν ( k ) . (151)We implemented this operator in our packages as delta2diff(f,z,k) . One can use it to recoversome products of derivatives in (148). Example 22. (%i1) delta2diff(F(z),z,3); (% o1 ) F ( z ) (%i2) delta2diff(F(z),z,4); (% o2 ) dd z · F ( z ) (%i3) delta2diff(F(z),z,5); (% o3 ) F ( z ) · dd z · F ( z ) ! %i4) delta2diff(F(z),z,6); (% o4 ) dd z · F ( z ) ! (%i5) delta2diff(F(z),z,14); (% o5 ) d d z · F ( z ) ! · d d z · F ( z ) ! Using δ , z instead of ddz in Koepf’s original approach for holonomic functions yields a proce-dure to compute homogeneous quadratic differential equations generally of lowest order satisfiedby a given expression f ( z ). We therefore obtain the following algorithm. Algorithm 6
Computing a quadratic DE satisfied by an expression f Input:
An expression f ( z ). Output:
A quadratic differential equation over K of least order satisfied by f ( z ).1. If f = and we stop.2. f ,
0, compute A ( z ) = δ , z f ( z ) / f ( z ) , (1-a) if A ( z ) ∈ K ( z ) i.e A ( z ) = P ( z ) / Q ( z ) where P and Q are polynomials, then we havefound a quadratic DE satisfied by f : Q ( z ) F ( z ) − P ( z ) F ( z ) = . (1-b) If A ( z ) < K ( z ), then go to 3.3. Fix a number QN max ∈ N , the maximal order of the DE sought; a suitable value is QN max : =
19 which corresponds to the maximum δ , z -order for having a quadratic dif-ferential equation of forth order.(3-a) set N : = δ N + , z f ;(3-c) expand the ansatz δ N + , z f ( z ) + A N − δ N + , z f ( z ) + · · · + A f ( z ) = E X i = S i , (152)in elementary summands with A N , A N − , . . . , A as unknowns. E > N is the totalnumber of summands S i obtained after expansion. A differential equation of order zero: F =
0, since we assumed that d dz f = f . lgorithm 6 Computing a quadratic DE satisfied by an expression f (3-d) For each pair of summands S i and S j (0 i , j E ), group them additively togetherif there exists r ( z ) = S i ( z ) / S j ( z ) ∈ K ( z ). If the number of groups is N then we have N linearly independent expressions. In that case, there exists a solution which can befound by equating each group to zero. The resulting system is linear for the unknowns A , A , . . . , A N − . Solving this system gives rational functions in z , and the solution isunique since we normalized A N =
1. After multiplication by the common denominatorof the values found for A ( z ) , A ( z ) , . . . A N − ( z ) we get the holonomic DE sought. Ifotherwise the number of groups is larger than N , then there is no solution and the step isnot successful.(3-e) If (3-d) is not successful, then increment N , and go back to (3-b), until N = QN max . Note that an algorithm for the general AODE case can be defined similarly provided anappropriate replacement of δ , z for derivation. Our Maxima package has an implementation ofAlgorithm 6 with the syntax QDE(f(z),F(z),[Type]) . The argument
Type is either
Inhomogeneous to allow the search for inhomogeneous quadratic DEs, or
Homogeneous by default to look forhomogeneous ones. We have also implemented the general AODE case as
NLDE(f(z),F(z)) (NLfor non-linear).
Example 23. (%i1) QDE(tan(z),F(z)); (% o1 ) d d z · F ( z ) − · F ( z ) · dd z · F ( z ) ! = (%i2) QDE(tan(z),F(z),Inhomogeneous); (% o2 ) dd z · F ( z ) − F ( z ) − = (%i3) QDE(sec(z)ˆk,F(z)); (% o3 ) k · F ( z ) · d d z · F ( z ) ! + ( − − k ) · dd z · F ( z ) ! − k · F ( z ) = (%i4) QDE(z/(exp(z)-1),F(z)); (% o4 ) z · dd z · F ( z ) ! + F ( z ) + ( z − · F ( z ) = (%i5) QDE(log(1+sin(z)),F(z)); (% o5 ) d d z · F ( z ) + dd z · F ( z ) ! · d d z · F ( z ) ! = (%i6) QDE(tan(z)ˆk,F(z)); (% o6 ) k · F ( z ) · d d z · F ( z ) ! − · (cid:16) · k − (cid:17) · dd z · F ( z ) ! · d d z · F ( z ) ! + · ( k − · (2 + k ) · d d z · F ( z ) ! + · k · F ( z ) · d d z · F ( z ) ! + · (cid:16) · k − (cid:17) · dd z · F ( z ) ! = NLDE in the latter example yields an AODE that does not depend on the exponent k ,but is not quadratic . (%i7) NLDE(tan(z)ˆk,F(z)); (% o7 ) F ( z ) · dd z · F ( z ) ! · d d z · F ( z ) ! − · F ( z ) · d d z · F ( z ) ! + dd z · F ( z ) ! · d d z · F ( z ) ! − · F ( z ) · dd z · F ( z ) ! = QDE , generally
NLDE generates differential equations of lower order but ofhigher degree. However, in terms of timings
QDE is faster and the computed differential equa-tions give much simpler recurrence equations than the outputs of
NLDE . Remark 24.
Observe that unlike the holonomic case where the existence and uniqueness of asolution to the Cauchy problem is quite immediate, in the algebraic case one needs to take intoaccount other important facts since the computed differential equations are not always explicit.However, using the implicit function theorem (see Krantz and Parks (2012)) on underlying alge-braic polynomials, the classical existence and uniqueness theorem (Teschl, 2012, Theorem 2.2)can be applied locally to uniquely determine the solution of outputs of QDE as a well choseninitial value problem.6.2. Converting quadratic differential equations to recurrence equations
We need a rewrite rule similar to (45) for every differential monomial in the expansion of aquadratic differential equation. Let f ( z ) be a power series with representation f ( z ) = ∞ X n = a n z n . It suffices to find a recurrence equation term that corresponds to the quadratic differential equa-tion term z p · f ( z ) ( i ) · f ( z ) ( j ) , i , j , p ∈ N > . (153)This is done in a similar manner as for (45). We have f ( z ) ( k ) = ∞ X n = ( n + k · a n + k · z n , ∀ k ∈ N > , therefore f ( z ) ( i ) · f ( z ) ( j ) = ∞ X n = ( n + i · a n + i · z n · ∞ X n = ( n + j · a n + j · z n = ∞ X n = n X k = ( k + i · a k + i · ( n − k + j · a n − k + j · z n , (154)by application of the Cauchy product rule. Finally multiplying (154) by z p yields the formula z p · f ( z ) ( i ) · f ( z ) ( j ) = ∞ X n = n − p X k = ( k + i · ( n − p − k + j · a k + i · a n − p − k + j · z n , (155)61nd the corresponding rewrite rule z p · f ( z ) ( i ) · f ( z ) ( j ) −→ n − p X k = ( k + i · ( n − p − k + j · a k + i · a n − p − k + j . (156)Observe that (156) is a rewrite rule for the power series coefficients of the given expression f . When dealing with inhomogeneous DEs, the constant term must be considered differently.This is the main reason why we prefer to work with homogeneous DEs.Thus a procedure to convert quadratic differential equations to recurrence equations followsimmediately. Our packages contains the function FindQRE(f,z,a[n]) as analogue of
FindRE forthe quadratic case.
Example 25. (%i1) FindQRE(tan(z),z,a[n]); (% o1 ) (1 + n ) · (2 + n ) · a n + − · n X k = ( k + · a k + · a n − k = (%i2) FindQRE(z/(exp(z)-1),z,a[n]); (% o2 ) n X k = a k · a n − k + ( n − · a n + a n − = (%i3) FindQRE(log(1+sin(z)),z,a[n]); (% o3 ) n X k = ( k + · ( k + · a k + · ( n − k + · a n − k + + (1 + n ) · (2 + n ) · (3 + n ) · a n + = One can use the recurrence equation computed by
FindQRE to define a power series repre-sentation of a given holonomic or non-holonomic expression. For a given recurrence equationcomputed using
FindQRE , we write the highest order term in terms of the others. And evaluatingthe recurrence equation at some integers allows to determine the necessary initial values of therepresentation.Note that to get the highest order term when there are terms with symbolic sums in therecurrence equation, we remove parts corresponding to the minimum and the maximum value ofthe summation variable and substitute the initial conditions until we get a non-zero expressionfrom which the highest order term can be obtained. For example, to get the highest order term ofthe recurrence equation of z / (exp( z ) − n X k = a k · a n − k + ( n − · a n + a n − = , (157)we remove those parts of the symbolic sum corresponding to k = k = n . This gives n − X k = a k · a n − k + · a · a n + ( n − · a n + a n − = . (158)62hen we substitute the value of a and write the resulting highest order term in terms of the othersummands of the equation. In this example a n is necessarily the highest order term to be usedsince the part of the equation that has no sum always depends on a n after substitution of the valueof a .Observe that this process of determining the highest order term can quite easily be managedin the quadratic case because every term in the computed recurrence equation has at most onesum. In the general AODE case we generally have a more complicated situation where somesummands of the equation have many summation symbols.The FPS command of our Maxima package combines all the procedures of this section as thelast method to determine the power series representation of a given expression. Example 26. (%i1) FPS(log(1+sin(z)),z,n); (% o1 ) " ∞ X n = A n · z n , A n + = n + · ( n + · ( n + ( n + · A n + − (2 + n ) · (3 + n ) · A n + − n X k = ( k + · ( k + · A k + · ( n − k + · A n − k + ! , n > = , " A = , A = , A = − , A = (%i2) FPS(1/(1+sin(z)),z,n); (% o2 ) " ∞ X n = A n · z n , A n + = · A n + · P n − k = A k · A n − k ( n + · ( n + , n > = , [ A = , A = − (159)Note that these outputs can also be used to compute Taylor polynomials. We implementedit as QTaylor(f(z),z,z0,N) . However, due to the presence of summation terms the quadratic timecomplexity cannot be avoided and hence the code is generally slower than the built-in Maximacommand taylor . Nevertheless one can use a remembering program so that many calls of closeorders of the same function require a timing only for the first call.
Example 27. (%i1) taylor(log(1+sin(z)),z,0,10); (% o1 ) / T / z − z + z − z + z − z + · z − · z + · z − · z + ... (%i2) QTaylor(log(1+sin(z)),z,0,10); (% o2 ) − · z + · z − · z + · z − z + z − z + z − z + z (%i3) taylor(log(1+sin(z)),z,0,300)$ Evaluation took 1.2500 seconds (1.2510 elapsed) using 287.200 MB. (%i4) QTaylor(log(1+sin(z)),z,0,300)$
Evaluation took 33.5160 seconds (33.5400 elapsed) using 2345.257 MB. %i5) taylor(log(1+sin(z)),z,0,200)$ Evaluation took 0.3280 seconds (0.3410 elapsed) using 88.334 MB (%i6) QTaylor(log(1+sin(z)),z,0,200)$
Evaluation took 0.0000 seconds (0.0130 elapsed) using 1.062 MB.
Notice that our computations of recurrence equations and Taylor polynomials from quadraticdifferential equations sketch another proof on the existence and uniqueness of the solutions ofthese differential equations without a use of the implicit function theorem. Therefore, our normalforms are well defined. Note, however, that we do not consider Puiseux series in these computa-tions because our approach to determine Puiseux numbers from recurrence equations generatedby
FindQRE cannot apply. This could be a topic for further studies on non-holonomic Puiseuxseries. The algebraic geometry approach described in (Falkensteiner (2020)) could be of greathelp for this purpose.
In this section, we present some new results occurring as consequences using our algorithmof the previous subsection. This is an improvement towards identity proving for non-holonomicexpressions. This advantage of symbolic computing was well discussed in the case of hypergeo-metric identities in (Petkovˇsek et al. (1996)). Up to now our algorithms were meant to representthe power series of a given expression A; now we would ask our algorithm to answer the ques-tion whether two given expressions A and B have the same power series representations, andtherefore whether they are identical in a certain neighborhood. This is the conclusion when twoexpressions have the same output using our algorithm or when the representation of their differ-ence is zero. We present two examples where this can be observed. In this subsection we use ourMaple implementation in order to compare our computations with its built-in command simplify which apply many internal simplification rules on its given inputs.As first example, the expressionlog (cid:18) tan (cid:18) z (cid:19) + sec (cid:18) z (cid:19)(cid:19) − arcsinh sin( z )1 + cos( z ) ! (160)from (Geddes et al., 1992, Section 3.3) (see also (Koepf, 2006, Exercise 9.8)) is known to bedifficult to prove equal to zero. One needs non-trivial transformations to simplify this to zero.However, using our algorithm based on the computation of quadratic differential equations yieldsthe same power series representation for > f:=ln(tan(z/2)+sec(z/2)) f : = ln (tan ( z / + sec ( z / > g:=arcsinh(sin(z)/(cos(z)+1)) g : = arcsinh sin ( z )cos ( z ) + ! as shown below. 64 FPS[FPS](f,z,n,fpstype=Quadratic) ∞ X n = A ( n ) z n , A ( n + = −
12 ( n +
2) ( n +
3) ( n + − ( n + A ( n + + n X k = k +
1) ( k +
2) ( k + A ( k +
3) ( n − k + A ( n − k + + n X k = ( − ( k + A ( k +
1) ( n − k + A ( n − k + + n X k = − k +
1) ( k + A ( k +
2) ( n − k +
2) ( n + − k ) A ( n + − k ) ! , ≤ n , " A (0) = , A (1) = , A (2) = , A (3) = > FPS[FPS](g,z,n,fpstype=Quadratic) ∞ X n = A ( n ) z n , A ( n + = −
12 ( n +
2) ( n +
3) ( n + − ( n + A ( n + + n X k = k +
1) ( k +
2) ( k + A ( k +
3) ( n − k + A ( n − k + + n X k = ( − ( k + A ( k +
1) ( n − k + A ( n − k + + n X k = − k +
1) ( k + A ( k +
2) ( n − k +
2) ( n + − k ) A ( n + − k ) ! , ≤ n , " A (0) = , A (1) = , A (2) = , A (3) = The argument fpstype=Quadratic is used to apply our method to non-holonomic functionsdirectly. In fact, our algorithm
QDE computes the same differential equation for f , g , and f − g . That is > FPS[QDE](f-g,F(z)) − dd z F ( z ) ! − d d z F ( z ) ! + d d z F ( z ) ! dd z F ( z ) = . (161)Moreover our implementation identifies their difference to zero, which is inaccessible usingthe built-in command simplify . > CPUTime(FPS[FPS](f-g,z,n,fpstype=Quadratic)) . , .
297 indicates the CPU time used. 65 simplify(f-g,trig) ln + sin ( z / z / ! − arcsinh sin ( z )cos ( z ) + ! Next, we consider log + tan( z )1 − tan( z ) ! − sin(2 z )1 + cos(2 z ) ! . (162)Similarly we get the following computations. > f:=ln((1+tan(z))/(1-tan(z))) f : = ln + tan ( z )1 − tan ( z ) ! > g:=2*arctanh(sin(2*z)/(1+cos(2*z))) g : = sin (2 z )cos (2 z ) + ! f and g satisfy the same differential equation > FPS[QDE](f,F(z)) − dd z F ( z ) ! − d d z F ( z ) ! + d d z F ( z ) ! dd z F ( z ) = . Therefore we get the same power series representation > FPS[FPS](g,z,n,fpstype=Quadratic) ∞ X n = A ( n ) z n , A ( n + = −
12 ( n +
2) ( n +
3) ( n + − n + A ( n + + n X k = ( k +
1) ( k +
2) ( k + A ( k +
3) ( n − k + A ( n − k + + n X k = − k + A ( k +
1) ( n − k + A ( n − k + + n X k = − k +
1) ( k + A ( k +
2) ( n − k +
2) ( n + − k ) A ( n + − k ) ! , ≤ n , " A (0) = , A (1) = , A (2) = , A (3) = (163) However in this case f − g satisfies a trivial differential equation, which makes the computa-tions much faster. > FPS[QDE](f-g,F(z)) dd z F ( z ) = CPUTime(FPS[FPS](f-g,z,n,fpstype=Quadratic)) . , simplify . > simplify(f-g) ln cos ( z ) + sin ( z )cos ( z ) − sin ( z ) ! − sin (2 z )cos (2 z ) + ! Acknowledgment 1.
This work summarizes the results of the first author’s Ph.D. thesis underthe supervision of the second. The first author received a DAAD STIBET and two DAAD ErasmusPlus scholarships to finish his work. Therefore, both authors would like to thank DAAD and theInstitute of Mathematics of the University of Kassel for their important and valuable support.
References
Brewer, T., 2014. Algebraic properties of formal power series composition. Ph.D. thesis, University of Kentucky, https://uknowledge.uky.edu/cgi/viewcontent.cgi?article=1021&context=math_etds .Cluzeau, T., van Hoeij, M., 2006. Computing hypergeometric solutions of linear recurrence equations. Appl. AlgebraEngrg. Comm. Comput. 17 (2), 83–115.Eremenko, A. E., 1982. Meromorphic solutions of algebraic differential equations. Russian Mathematical Surveys 37 (4),61–95.Falkensteiner, S., 2020. Power series solutions of AODEs – existence, uniqueness, convergence and computation. Ph.D.thesis, RISC Hagenberg, Johannes Kepler University Linz.Geddes, K. O., Czapor, S. R., Labahn, G., 1992. Algorithms for Computer Algebra. Kluwer Academic Publishers,Massachusetts.Gruntz, D., Koepf, W., 1995. Maple package on formal power series. Maple Technical Newsletter 2 (2), 22–28.Heck, A., 2003. Introduction to Maple, 3rd Edition. Springer-Verlag, New York.Hendricks, P. A., Singer, M. F., 1999. Solving difference equations in finite terms. J. Symb. Comput. 27 (3), 239–259.Horn, P., Koepf, W., Sprenger, T., 2012. m -fold hypergeometric solutions of linear recurrence equations revisited. Math.Comput. Sci. 6 (1), 61–77.Kauers, M., Paule, P., 2011. The Concrete Tetrahedron. Symbolic Sums, Recurrence Equations, Generating Functions,Asymptotic Estimates. Springer-Verlag, Wien.Koepf, W., 1992. Power series in computer algebra. J. Symb. Comput. 13 (6), 581–603.Koepf, W., 1995. Algorithms for m -fold hypergeometric summation. J. Symb. Comput. 20 (4), 399–417.Koepf, W., 1997. The algebra of holonomic equations. Math. Semesterber. 44, 173–194.Koepf, W., 2006. Computeralgebra: eine algorithmisch orientierte Einf¨uhrung. Springer-Verlag, Berlin Heidelberg, NewYork.Koepf, W., 2014. Hypergeometric Summation, An Algorithmic Approach to Summation and Special Function Identities,2nd Edition. Springer-Verlag, London.Koepf, W., Schmersau, D., 1998. Representations of orthogonal polynomials. J. Comput. Appl. Math. 90 (1), 57–94.Krantz, S. G., Parks, H. R., 2012. The Implicit Function Theorem: History, Theory, and Applications. Birkh¨auser,Boston.Lubin, J., 1994. Nonarchimedean dynamical systems. Compos. Math. 94 (3), 321–346.Petkovˇsek, M., 1992. Hypergeometric solutions of linear recurrences with polynomial coefficients. J. Symb. Comput.14 (2-3), 243–264.Petkovˇsek, M., Salvy, B., 1993. Finding all hypergeometric solutions of linear differential equations. In: ISSAC. Editor:Bronstein, Manuel, Association for Computing Machinery, New York, pp. 27–33.Petkovˇsek, M., Wilf, H. S., Zeilberger, D., 1996. A=B. Vol. 30. AK Peters Ltd, Wellesley, MA.Stanley, R. P., 1980. Differentiably finite power series. European J. Combin. 1 (2), 175–188.Stanley, R. P., 2011. Enumerative Combinatorics, 2nd Edition. Vol. 1. Cambridge Studies in Advanced Mathematics,Cambridge. eguia Tabuguia, B., 2020a. Power series representations of hypergeometric type andnon-holonomic functions in computer algebra. Ph.D. thesis, University of Kassel, https://kobra.uni-kassel.de/handle/123456789/11598 .Teguia Tabuguia, B., 2020b. A variant of van Hoeij’s algorithm to compute hypergeometric term solutions of holonomicrecurrence equations. arXiv:2012.11513 [cs.SC].Teguia Tabuguia, B., Koepf, W., 2020. Power series representations of hypergeometric functions (to appear). In: MapleConference 2020 Proceedings. Editors: Corless, Rob and J¨urgen, Gerhard, Springer.Teschl, G., 2012. Ordinary Differential Equations and Dynamical Systems. Vol. 140. American Mathematical Society,Providence.Van Hoeij, M., 1999. Finite singularities and hypergeometric solutions of linear recurrence equations. J. Pure Appl.Algebra 139 (1-3), 109–131.Wolfram, S., 2003. The Mathematica Book, version 4, 5th Edition. Wolfram Media, Cambridge University Press.Zariski, O., Samuel, P., 1960. Commutative Algebra. Vol. 2. Springer-Verlag, New York..Teguia Tabuguia, B., 2020b. A variant of van Hoeij’s algorithm to compute hypergeometric term solutions of holonomicrecurrence equations. arXiv:2012.11513 [cs.SC].Teguia Tabuguia, B., Koepf, W., 2020. Power series representations of hypergeometric functions (to appear). In: MapleConference 2020 Proceedings. Editors: Corless, Rob and J¨urgen, Gerhard, Springer.Teschl, G., 2012. Ordinary Differential Equations and Dynamical Systems. Vol. 140. American Mathematical Society,Providence.Van Hoeij, M., 1999. Finite singularities and hypergeometric solutions of linear recurrence equations. J. Pure Appl.Algebra 139 (1-3), 109–131.Wolfram, S., 2003. The Mathematica Book, version 4, 5th Edition. Wolfram Media, Cambridge University Press.Zariski, O., Samuel, P., 1960. Commutative Algebra. Vol. 2. Springer-Verlag, New York.