A variant of van Hoeij's algorithm to compute hypergeometric term solutions of holonomic recurrence equations
aa r X i v : . [ c s . S C ] D ec A VARIANT OF VAN HOEIJ’S ALGORITHM TO COMPUTEHYPERGEOMETRIC TERM SOLUTIONS OF HOLONOMIC RECURRENCEEQUATIONS
BERTRAND TEGUIA TABUGUIAA
BSTRACT . Linear homogeneous recurrence equations with polynomial coefficients aresaid to be holonomic. Such equations have been introduced in the last century for provingand discovering combinatorial and hypergeometric identities. Given a field K of character-istic zero, a term a n is called hypergeometric with respect to K , if the ratio a n +1 /a n isa rational function over K . The solutions space of holonomic recurrence equations gainedmore interest in the 1990s from the well known Zeilberger’s algorithm. In particular, al-gorithms computing the subspace of hypergeometric term solutions which covers polyno-mial, rational, and some algebraic solutions of these equations were investigated by MarkoPetkovšek (1993) and Mark van Hoeij (1999). The algorithm proposed by the latter ischaracterized by a much better efficiency than that of the other; it computes, in Gammarepresentations, a basis of the subspace of hypergeometric term solutions of any givenholonomic recurrence equation, and is considered as the current state of the art in this area.Mark van Hoeij implemented his algorithm in the Computer Algebra System (CAS) Maplethrough the command LREtools[hypergeomsols] .We propose a variant of van Hoeij’s algorithm that performs the same efficiency andgives outputs in terms of factorials and shifted factorials, without considering certain rec-ommendations of the original version. We have implementations of our algorithm for theCASs Maxima and Maple. Such an implementation is new for Maxima which is thereforeused for general-purpose examples. Our Maxima code is currently available as a third-party package for Maxima. A comparison between van Hoeij’s implementation and oursis presented for Maple 2020. It appears that both have the same efficiency, and moreover,for some particular cases, our code finds results where
LREtools[hypergeomsols] fails.
1. I
NTRODUCTION
Let K be a field of characteristic zero. K is mostly a finite extension field of the rationals.A hypergeometric term can always be written in the form(1.1) C n · R ( n ) · h ( n ) , where C ∈ K , R ( n ) ∈ K ( n ) , and h ( n ) is a hypergeometric term expressed in terms offactorials and shifted factorials (Pochhammer symbols ) such that h ( n + 1) /h ( n ) ∈ K ( n ) is monic (see [5], [16, Chapter 6]). Notice that the representation (1 . is unique if wechoose to write Pochhammer terms as ( p ) n with the real part of p in a fixed half-open realinterval of unit amplitude. This rewriting creates a multiplicative rational function factorthat is taken into account when computing R ( n ) . We will consider the interval I = (0 , Mathematics Subject Classification.
Primary 33F10, 39A06; Secondary 33C20, 68W30.
Key words and phrases.
Holonomic recurrence equations, hypergeometric terms, van Hoeij’s algorithm,Petkovšek’s algorithm.This work gather results from chapters and of author’s Ph.D. thesis, [16]. The author is grateful to hisadvisor Wolfram Koepf for invaluable guidance. For a given constant p , the Pochhammer symbol ( p ) n is if n = 0 and p · ( p + 1) · · · ( p + n − if n is apositive integer. in our algorithm, and say that Pochhammer parts, corresponding to h ( n ) , is taken modulothe integers ( Z ) with respect to I . Example 1.1. n n ! and n n +1 n +2 (1 / n (3 / n are two hypergeometric terms of the form (1 . .A homogeneous linear recurrence equation with rational function coefficients is equiv-alent to a linear recurrence equation with polynomial coefficients, considering multipli-cation of the equation by the corresponding rational coefficients common denominator.Therefore the algorithms we are dealing with easily adapt to such cases. However, for thispaper we will consider recurrence equations of the form(1.2) d X i =0 P i ( n ) · a n + i = 0 , d ∈ N , for the indeterminate sequence a n , and the coefficients P i ( n ) ∈ K [ n ] , i = 0 , . . . , d . Twomain algorithms were proposed to find all its hypergeometric term solutions.Petkovšek’s approach presented in [12], focuses on the computation of ratios of hyper-geometric term solutions of (1 . and look for formulas afterward. However, his algorithmhas an exponential worst-case complexity on the degree of P and P d . Thus this approachcould not be considered as conclusive for such computations. The commands solve_rec ofthe CAS Maxima and originally RSolve of Mathematica implement Petkovšek’s algorithm.van Hoeij’s approach first described in [18] is much more efficient, and moreover, findsa basis of the subspace of hypergeometric term solutions of (1 . . Note that the finaloutput in Petkovšek’s algorithm is not necessarily a basis. Therefore one could figure outthe gain of efficiency by observing that many computational cases in Petkovšek’s approachare reduced to one case in van Hoeij’s approach. Computational details of van Hoeij’salgorithm were more explained and complemented in [5]. Another important point tonotice in this algorithm is how unnecessary splitting fields that increase the running timeduring computations are avoided. Definition 1.2. (see [18, Definition 9], [5, Definition 8]) A point p + Z , p ∈ K is calledfinite singularity of (1 . if there exists τ ∈ Z such that p + τ is a root of P d ( n − d ) · P ( n ) .From this definition where invariance modulo the integers is clearly put forward, onecan see the connection between finite singularities and our used rewriting of hypergeomet-ric term Pochhammer parts.The power of van Hoeij’s algorithm comes from the following main concepts:(1) local types at infinity of hypergeometric term solutions of (1 . ,(2) local types or valuation growths of hypergeometric term solutions at finite singu-larities of (1 . .The computations of (1) and (2) constitute the key steps of van Hoeij’s algorithm and thatis where our approach proceeds differently. • For (1), van Hoeij’s algorithm uses the Newton polygon algorithm whereas we usea method based on asymptotic expansion inspired by Petkovšek’s algorithm Poly(see [12]). • Computing (2) is inherent in van Hoeij’s algorithm, but in our approach this isautomatically considered in the way we construct h ( n ) in (1 . by taking monicfactors modulo the integers of P and P d ( n − d ) . This consideration is valid thanksto Petkovšek’s approach. a n +1 /a n for a hypergeometric term a n . This notion was introduced to study the local properties of difference operators at infinity (see [4, 6])
VARIANT OF VAN HOEIJ’S ALGORITHM TO COMPUTE HYPERGEOMETRIC TERMS 3
Apart from these essential differences, it is not trivial to notice that both algorithms do thesame thing, because their step orderings do not coincide either.Our first motivation in implementing van Hoeij’s algorithm is to use the outputs forpower series representations as described in [9, Section 10.26]. Indeed, one can representsome holonomic functions as linear combinations of hypergeometric series whose coeffi-cients are hypergeometric term solutions of underlying holonomic recurrence equations.For this purpose, we want outputs ready for evaluation for non-negative integers. Thisis settled in our approach by considering Pochhammer parts modulo the integers with re-spect to I . Contrary to van Hoeij’s Maple implementation LREtools[hypergeomsols] , thischoice is independent of the recurrence equation considered.
LREtools[hypergeomsols] uses Gamma representations for h ( n ) , and similarly as Pochhammer symbols, evaluationis not defined for negative integers. Since such normalization is not taken into account in LREtools[hypergeomsols] , for use in power series computations or combinatory (see [10]),one may have to shift the initialization which might even be different for each hypergeo-metric term appearing in the given output. In other cases it may lead to inconvenient resultswith evaluation like
Γ(1 /
2) = √ π . As we mentioned earlier, we want h ( n ) in terms offactorials and Pochhammer symbols modulo the integers with respect to I ; we will saythat such a formula is "simple".Let us give an example to illustrate all these. Example 1.3.
Consider the following holonomic recurrence equation(1.3) n (cid:0) − n (cid:1)(cid:0) n + 56592 n + 566784 n + 3438888 n + 14040866 n + 40413165 n + 83014167 n + 118689722 n + 105269208 n + 24761376 n − n − n − n − n − n − (cid:1)(cid:0) − n (cid:1) (cid:0) n (cid:1) a n − (cid:0) − n (cid:1)(cid:0) n + 165722112 n + 1895913216 n + 13287379968 n + 63281637504 n + 213327813888 n + 505402785504 n + 757111794432 n + 271146179476 n − n − n − n − n − n − n − n − n − n − n − n − n − n − n − n − (cid:1)(cid:0) n + 1 (cid:1) a n +1 + 32 (cid:0) n + 17712 n + 46656 n + 41208 n − n − n − n − n − n − n − n − n − n − n − n − (cid:1)(cid:0) n + 2 (cid:1) (cid:0) n + 4 (cid:1) a n +2 = 0 Our Maple and Maxima implementation finds the following output with CPU times . and . second respectively.(1.4) ( n ! (cid:0) (cid:1) n ( n − n , n (2 n ) ! ( n −
2) ( n −
1) 4 n n ! ) . Observe that evaluation can be made for non-negative integers. As one can observe, furthersimplification are made in our implementation to transform Pochhammer symbols intofactorials.
BERTRAND TEGUIA TABUGUIA
Maple 2019
LREtools[hypergeomsols] finds(1.5) " Γ ( n −
2) (Γ ( n + 1 / n , (Γ ( n − Γ ( n −
2) ( n − n + 1 / n , with CPU time . second. In this output the Gamma terms Γ ( n − and Γ ( n − cannot be evaluated at . Moreover their arguments differs by , which shows that shiftsmost be considered before initialization. Note, however, that this remark is only for furtheruse of hypergeometric terms, one can easily show that (1 . and (1 . are equivalent.The Maxima (version 5.44) command solve_rec finds(1.6) a n = Γ (cid:0) (cid:1) % k ( n −
2) ! − n − n n Γ (cid:0) n +13 (cid:1) + % k ( n −
3) ! n n +15 Γ (cid:0) n +12 (cid:1) π n , with CPU time . seconds. As one may expect, finding formulas from the obtainedratios of hypergeometric terms with Petkovšek’s algorithm is more complicated. Never-theless, the large computation time here is due to the degrees of the leading and trailingpolynomial coefficients of (1 . .This paper goes as follows. In the next section we derive an algorithm to computeholonomic recurrence equations satisfied by a list of linearly independent hypergeometricterms. This algorithm is useful to generate examples and observe some properties of hy-pergeometric terms by doing forward and backward computations. This can also be doneusing the Maple package gfun (see [14]), but note that this code is limited to two recurrenceequations as input and the strategy used is slightly different.In Section 3 we give more details on how we normalize the Pochhammer parts of hy-pergeometric terms. This will help to make the description of the main algorithm of thispaper self-contained.Section 4 describes our variant of van Hoeij’s algorithm which efficiently computesa basis of the subspace hypergeometric term solutions of (1 . , using the representation (1 . . The paper ends with some comparisons with existing implementations.2. F ROM HYPERGEOMETRIC TERMS TO HOLONOMIC RECURRENCE EQUATIONS
It is well known that linear combinations of holonomic functions are holonomic (see [9,Section 10.9, Section 10.16], [15]). Since hypergeometric terms are holonomic, there existalgorithms to compute a holonomic recurrence equation of least order satisfied by a givenlinear combination of hypergeometric terms. Throughout this section we assume thereexists an algorithm for finding the rational function defined by the ratio of a hypergeometricterm (see [10, Algorithm 2.2]). The algorithm of this section is a generalization of the caseof two given linearly independent hypergeometric terms. Thus, we treat this particular caseand by simple analogy we give the general approach for a given list of linearly independenthypergeometric terms.2.1. Case of two linearly independent hypergeometric terms.
Let a n and b n be twolinearly independent hypergeometric terms over K such that(2.1) a n +1 = r ( n ) a n and b n +1 = r ( n ) b n , Note that this generally reduces to a sum of hypergeometric terms. Therefore the main information here isthat the given hypergeometric terms are distinct.
VARIANT OF VAN HOEIJ’S ALGORITHM TO COMPUTE HYPERGEOMETRIC TERMS 5 where r and r are rational functions in K ( n ) . As we consider two terms, the order of therecurrence equation sought is , so we are looking for a recurrence equation of the form(2.2) P ( n ) s n +2 + P ( n ) s n +1 + P ( n ) s n = 0 , where P , P , P are polynomials over K , satisfied by a n and b n . We must assume that P · P = 0 , otherwise the recurrence equation can be reduced to a first order recurrencerelation. Thus finding (2 . is equivalent to searching for rational functions R and R such that(2.3) R ( n ) s n +2 + R ( n ) s n +1 + s n = 0 . Using (2 . , we have(2.4) a n +2 = r ( n + 1) a n +1 and b n +2 = r ( n + 1) b n +1 . By substitution, a n and b n satisfy (2 . if and only if(2.5) r ( n + 1) R + R = − r ( n ) r ( n + 1) R + R = − r ( n ) , which is a linear system of two equations with two unknowns in K ( n ) . Furthermore, asolution exists and is unique since the determinant of the system(2.6) r a ( n + 1) − r b ( n + 1) = 0 by assumption. As a linear system of two equations, the exact solution is easy to compute,that is R ( n ) = r ( n + 1) r ( n ) − r ( n + 1) r ( n ) r ( n ) r ( n )( r ( n + 1) − r ( n + 1)) , (2.7) R ( n ) = r ( n ) − r ( n ) r ( n ) r ( n )( r ( n + 1) − r ( n + 1)) . (2.8)Finally, the holonomic recurrence equation sought is found by multiplying the equation (2 . by the common denominator of R ( n ) and R ( n ) .2.2. General case.
Now we want to generalize the above approach for finitely many lin-early independent hypergeometric terms. Let a [ i ] n , i = 1 , . . . , d ( d > be d given linearlyindependent hypergeometric terms over K such that(2.9) a [ i ] n +1 = r i ( n ) a [ i ] n , i = 1 , . . . , d, for some rational functions r i . The vector ( R ( n ) , R ( n ) , . . . , R d ( n )) T ∈ K ( n ) d of ratio-nal coefficients of the recurrence equation(2.10) R d ( n ) s n + d + R d − ( n ) s n + d − + . . . + R ( n ) s n +1 + s n = 0 satisfied by each hypergeometric term a [ i ] n , is the unique vector solution v ∈ K ( n ) d of thematrix system(2.11) " j − Y k =1 r i ( n + k ) i,j =1 ,...,d · v = − (cid:18) r i ( n ) (cid:19) Ti =1 ,...,d . In case there are linearly dependent hypergeometric terms, one can still use this process byreplacing the arbitrary constants appearing in the solution of (2 . by zero. This is howwe implemented this method. The steps of the algorithm can be summarized as follows. BERTRAND TEGUIA TABUGUIA
Algorithm 1
Compute the holonomic recurrence equation of least order for a given list L of hypergeometric terms Input:
A list L := [ h , . . . , h d ] of hypergeometric terms in the variable n and a symbol a . Output:
A holonomic recurrence equation in a n of least order satisfied by the elementsin L .(1) Let R := [ r i ( n ) , . . . , r d ( n )] be the ratios of the elements in L .(2) If there are some irrational functions in R then stop and return FALSE. No holo-nomic recurrence equation can be found.(3) Let M := " j − Y k =1 r i ( n + k ) i,j =1 ,...,d . (4) Let b := (cid:18) r i ( n ) (cid:19) Ti =1 ,...,d . (5) Let V be the solution of the matrix system M · v = b .(6) If there are arbitrary constants in V then substitute those by zero.(7) Let RE := a n + P di =1 V [ i ] · a n + i , where V [ i ] denotes the i th component in V .(8) Multiply RE by the common denominator of the components of V and return theresult with equality to 0, after factoring the coefficients. Example 2.1.
We implemented this algorithm as sumhyperRE . Let us consider Example4.1 in [12] and make a backward computation to find the recurrence equation for n + 1)( n + 2) , and ( − n (2 n + 3)( n + 1)( n + 2) . Our Maxima code gets (% i1) sumhyperRE([1/((n+1)*(n+2)), (-1)ˆn*(2*n+3)/((n+1)*(n+2))],a[n]);(% o1) − ( n + 4) a n +2 − a n +1 + ( n + 1) a n = 0 , which is the expected result. Next we recover the Fibonacci recurrence from the goldennumber and its conjugate. (% i2) sumhyperRE([(1-sqrt(5))ˆn/2ˆn, (1+sqrt(5))ˆn/2ˆn], a[n]);(% o2) − a n +2 + a n +1 + a n = 0 The latter example illustrates an important point of the algorithm. Indeed, when con-sidering extension fields to determine hypergeometric term solutions, the conjugates ofalgebraic numbers involved are also part of the solution basis. We will give more detailsabout this in Section 4.3. "S
IMPLE " FORMULAS FOR HYPERGEOMETRIC TERMS
Let a n be a hypergeometric term over a field K of characteristic zero. Then by defini-tion r ( n ) := a n +1 /a n ∈ K ( n ) . K is taken as the minimal extension field of Q where thenumerator and the denominator of r ( n ) split. We wish to write the formula of a n by meansof only numbers appearing in the splitting field of r ( n ) using factorials, and when not triv-ially possible, Pochhammer symbols. We want, moreover, that evaluations of the formulafor non-negative integers do not produce other improper operation (usually at ) than the VARIANT OF VAN HOEIJ’S ALGORITHM TO COMPUTE HYPERGEOMETRIC TERMS 7 division by zero. Nevertheless, the representation we are tagging is mostly valid for evalu-ations with for non-negative integers. This is what we call a "simple" formula. One couldsay that a formula is considered to be "simple" when it presents more familiar objects frommathematical dictionaries in a reduced form. In the sense of computing formulas of hyper-geometric terms, this consists of simplifying as much as possible, Pochhammer symbolsto rational multiples of factorials with positive integer-linear arguments. In this section,we present preliminary steps to recover the representation (1 . and gather some classicalrules as an algorithm to simplify its Pochhammer part. Similar computations can be foundin [8]; what is worth to notice is the consideration we make to get "simple" formulas inSection 4.Consider a rational function(3.1) r ( k ) := P ( k ) Q ( k ) , P ( k ) , Q ( k ) ∈ K [ k ] , Q ( k ) = 0 for integers k > such that P and Q do not have non-negative integer roots. We also assume that rootsof P and Q are all distinct. We will see in the next section how our approach preparesall rational functions used to compute hypergeometric term Pochhammer part "simple"formulas to satisfy these assumptions. For example, a rational function r ( k ) with non-negative integer zeros and poles would implicitly be replaced by r ( k + m ) , where m =max { j ∈ N > : Q ( j ) · P ( j ) = 0 } .We consider a hypergeometric term defined with the property(3.2) a k +1 = r ( k ) a k , for integer k > . Computing a "simple" formula of such a term is to find its general expression a n for apositive integer n provided that the corresponding initial values a is given. That is theresult of the product(3.3) n − Y k =0 r ( k ) . For that purpose, the first step is to split the polynomials of r as follows(3.4) r ( k ) = C ( k + a )( k + a ) · · · ( k + a p )( k + b )( k + b ) · · · ( k + b q ) , where p and q are, respectively, the degrees of P and Q ; C is a constant representing theratio of the leading coefficients of P and Q ; and − a i ’s, i p , and − b j ’s, j q are the zeros and poles of r , respectively. From the Pochhammer symbol definition, using (3 . it follows that(3.5) n − Y k =0 r ( k ) = C n ( a ) n ( a ) n · · · ( a p ) n ( b ) n ( b ) n · · · ( b q ) n . Thus we are called to try simplifications of ratios and products of Pochhammer symbols,and some isolated ones. Many such computations can be found in books or undergraduatecourses, see for example [10, Exercises 1.1 - 1.5]. Bellow we recall some classical ones. x and y denote some numbers, and j an integer. Isolated Rule:
Assume x is rational, then one can simplify ( x ) n , according to thefollowing cases. N > = { , , , . . . } BERTRAND TEGUIA TABUGUIA (1) if x ∈ N , then ( x ) n = x · ( x + 1) · · · ( x + n − x + n − x − (3.6) (2) Else if x has a denominator equal to , then let s ∈ N such that x = s . s isnecessarily an odd integer since x / ∈ N . We set s = 2 t + 1 , t ∈ N > , then itfollows that ( x ) n = s · (cid:16) s (cid:17) · · · (cid:16) s n − (cid:17) = s · ( s + 2) · · · ( s + 2 · ( n − n = (2 t + 1) · (2( t + 1) + 1) · · · (2( t + n −
1) + 1)2 n = (2 ( t + n ))!(2 t )! · (2 t + 2) · · · (2( t + n −
1) + 2) · n = (2 ( t + n ))!(2 t )! · ( t + 1) · · · ( t + n ) · n = (2( t + n ))!(2 t )!4 n (cid:0) t + nn (cid:1) n ! . (3.7) (3) Otherwise no simplification is done for ( x ) n . Ratio Rule:
Assume x − y = j > . Then we have ( y ) n ( x ) n = ( y ) j · ( y + j ) · · · ( y + n − y + j ) · · · ( y + j + n − y ) j ( y + n ) · · · ( y + n + j − y ) j ( y + n ) j . (3.8) Therefore for x − y = j ∈ Z ,(3.9) ( y ) n ( x ) n = ( y ) j ( y + n ) j if j > ( x + n ) − j ( x ) − j if j < . This shows that differences between zeros and poles of r in (3 . should be checkedbefore applying the Isolated Rule in order to apply (3 . which can simplify twoPochhammer symbols at the same time. Fortunately, these nice computations canbe done in Maxima by combining makegamma() , makefact() , minfactorial() and factor() as below. (%i1) r:pochhammer(7/3,n)/pochhammer(1/3,n); (% o1 ) (cid:0) (cid:1) n (cid:0) (cid:1) n (%i2) factor(minfactorial(makefact(makegamma(r)))); (% o2 ) (1 + 3 · n ) · (4 + 3 · n )4 Product Rule:
We consider the following two rules to simplify ( x ) n ( y ) n . VARIANT OF VAN HOEIJ’S ALGORITHM TO COMPUTE HYPERGEOMETRIC TERMS 9 • Assume y − x = 1 / , then multiplying ( x ) n and ( y ) n = ( x + 1 / n by n leads to the relation(3.10) ( x ) n · (cid:18) x + 12 (cid:19) n = (2 x ) n n . • Assume y − x = j > , it is easy to see that(3.11) ( x ) n · ( y ) n = ( x ) n · ( x + j ) n = ( x ) n + j ( x ) j ( x + n ) j . More generally, one can find a hypergeometric term Pochhammer part "simple" formulahaving ratio with non-negative integer zeros and poles by applying the following algorithm.
Algorithm 2
Compute Q n − k =1 r ( k ) Input:
A rational function r := r ( n ) and a variable n . Output:
A formula of Q n − k =1 r ( k ) in terms of factorial and Pochhammer symbols.(1) Factorize r and write it in terms of linear factors and set h := r = C ( n + a )( n + a ) · · · ( n + a p )( n + b )( n + b ) · · · ( n + b q ) . (2) Substitute C by C n in h .(3) For each a i , i = 1 , . . . , p do(a) if there is b j in h such that a i − b j ∈ Z then substitute n + a i n + b i by applyingthe Ratio Rule accordingly.(4) For the remaining a i ’s (resp. b j ’s) do(a) if there is a i ′ (resp. b j ′ ) such that a i − a i ′ = ± / or a i − a i ′ ∈ Z (resp. a j − a j ′ = ± / or a j − a j ′ ∈ Z ) then substitute ( n + a i )( n + a i ′ ) (resp. ( n + b j )( n + b j ′ ) ) by applying the Product Rule accordingly.(5) Substitute the remaining n + a i ’s and n + b j ’s by the result of the Isolated Rule applied to ( a i ) n and ( b j ) n respectively.(6) Return h .For a "fair" comparison of efficiency between our implementation and the current Maple LREtools[hypergeomsols] , we did not considered the
Product Rule in our Maxima andMaple codes (but it will certainly be the case in future versions) since the internal Maplecommand
LREtools[hypergeomsols] does not apply simplifications as we do. So thesecomputations are supplementary steps that we use in the algorithm of Section 4. Algorithm2 is implemented in our Maxima package as pochfactorsimp(r,n) . Below we give someexamples.
Example 3.1. (%i1) pochfactorsimp(-1/(2*(n+1)*(2*n+1)),n); (% o1 ) ( − n (2 · n )! (%i2) pochfactorsimp((2*n+3)^2/((n+1)*(2*n+1)),n); (% o2 ) (1 + 2 · n ) · n − · (2 · (1 + n ))!( n + 1) · n · n ! Now that we have described an algorithm for computing "simple" formulas of hyperge-ometric terms given their ratios with no non-negative integer zeros and poles, let us moveto the main algorithm of this paper where such rational functions of the Pochhammer partsare computed. 4. B
ASIS OF HYPERGEOMETRIC TERM SOLUTIONS
Let us rewrite (1 . as follows.(4.1) P d ( n ) a n + d + P d − ( n ) a n + d − + · · · P ( n ) a n +1 + P ( n ) a n = 0 , with polynomials P i ( n ) ∈ K [ n ] , i = 0 , . . . , d such that P ( n ) · P d ( n ) = 0 . K is a field ofcharacteristic zero.We have seen how to compute a holonomic recurrence equation of lowest order satisfiedby a given number of linearly independent hypergeometric terms. Any computed hyperge-ometric term solution of such a holonomic recurrence equation is a linear combination ofthese linearly independent hypergeometric terms considered. The algorithm of this sectionis a kind of reverse process which for a given holonomic recurrence equation (4 . com-putes a basis of at most d hypergeometric terms of the space of all hypergeometric termsolutions of (4 . .In the first place, we establish (1 . to see hypergeometric terms in normal forms (see [7,Chapter 3]). Let a n , n ∈ N > , be a hypergeometric sequence such that r ( n ) = a n +1 /a n ∈ K ( n ) . Then we have a a = r (0) , a a = r (1) , . . . , a n a n − = r ( n − , n > , and therefore(4.2) a n a = n − Y k =0 a k +1 a k = n − Y k =0 r ( k ) ⇒ a n = a n − Y k =0 r ( k ) . Factorizing r ( n ) over K gives(4.3) r ( n ) = C Q Ii =1 ( n − α i ) Q Jj =1 ( n − β j ) , where C is a constant. Note that contrary to (3 . , in (4 . r ( n ) is considered in a moregeneral setting; α i and β j are not uniquely determined and may have negative or positivereal parts, which is more general than avoiding non-negative integer values.Combining (4 . and (4 . leads to(4.4) a n = a · C n · ( − α ) n · · · ( − α I ) n ( − β ) n · · · ( − β J ) n . VARIANT OF VAN HOEIJ’S ALGORITHM TO COMPUTE HYPERGEOMETRIC TERMS 11
Now we want to write each Pochhammer symbol modulo Z in a certain real interval.That is to say that the real parts of the arguments of Pochhammer terms can be chosenbelonging to an interval of amplitude . This is an interesting observation made by vanHoeij. In our case, we choose to rewrite the Pochhammer symbols modulo Z so that α i , β j ∈ I := [ − , . Each Pochhammer symbol is then substituted by a polynomialtimes another Pochhammer term whose argument differs by an integer u . Precisely, let y be a real number (for the case of complex numbers, the computations are applied on theirreal parts), then its corresponding value in I is u = y − ⌊ y ⌋ − and we have ( y ) n = ( u ) n · ( u + n ) · · · ( y + n − u · ( u + 1) · · · ( y − u ) n · ( u + n ) y − u ( u ) y − u (4.5) = ( y − ⌊ y ⌋ − n · ( n + y − ⌊ y ⌋ − ⌊ y ⌋ +1 ( y − ⌊ y ⌋ − ⌊ y ⌋ +1 . (4.6)After applying (4 . to each Pochhammer symbol in (4 . , the remaining expression willhave Pochhammer terms having arguments with real parts in (0 , . These terms mayhave more coincidence than the ( − α i ) n and ( − β j ) n in (4 . since all Pochhammer termsin (4 . whose arguments differ by an integer give the same Pochhammer term modulo Z after substitution. Therefore there exists a rational function R ( n ) ∈ K ( n ) and someconstant numbers ˜ α , . . . , ˜ α I , ˜ β , . . . , ˜ β J , with real parts in I , such that(4.7) a n = R ( n ) · C n · ( − ˜ α ) n · · · ( − ˜ α I ) n ( − ˜ β ) n · · · ( − ˜ β J ) n . The constant a is neglected by linearity since we will look for a basis of hypergeometricterm solutions of (4 . .Considering multiplicities e k over Z \ { } and replacing − ˜ α i and − ˜ β j , by θ k , we getthe normal form(4.8) a n = C n · R ( n ) · h ( n ) := C n · R ( n ) · K Y k =1 ( θ k ) e k n ( θ k ∈ K , with real part in I ) ,K I + J . This time all the involved data are uniquely determined. The ratio r ( n ) canbe rewritten as(4.9) r ( n ) = a n +1 a n = R ( n + 1) R ( n ) · C · h ( n + 1) /h ( n ) = R ( n + 1) R ( n ) · C · K Y k =1 ( n + θ k ) e k ∈ K ( n ) . Mark van Hoeij uses Gamma representations in (4 . and denotes it singularity structureof a n (see [18, 5]). This representation can be seen as the end point of our algorithm whenit computes an element of the basis of hypergeometric terms sought. In fact, the goal ofcomputing a basis of all hypergeometric term solutions of (4 . is equivalent to findingsolutions of (4 . with the structure (4 . .4.1. Monic factors modulo Z .Lemma 4.1. ( [12, Algorithm Hyper] , [5, Left and Right solutions] ) The ratio of thePochhammer part h ( n + 1) /h ( n ) of hypergeometric term solutions of (4 . are built frommonic factors of P d ( n − d ) for the numerators and P ( n − for the denominators. This lemma allows us to apply factorization modulo Z on P ( n ) and P d ( n ) . In factthe ratio of the Pochhammer part h ( n + 1) /h ( n ) in (4 . is obtained from factorization of P ( n ) and P d ( n ) modulo Z . Let us generate a recurrence equation that will be used whiledescribing the steps of our algorithm. (% i1) RE:sumhyperRE([binomial(n+3,n),1/n!,(-1)\^{}n/n,(-1)\^{}n/pochhammer(1/2,n)\^{}2],a[n])\$
We do not display the output to save space. We will refer to this recurrence equation as ( RE ) . The leading term is (% i2) first(lhs(RE));(% 2) ( n + 2) ( n + 3) ( n + 4) (2 n + 7) ( 64 n +1536 n +16176 n +98080 n +377372 n +955200 n +1584741 n +1631354 n +852544 n − n − n − a n +4 , and the trailing term (% i3) last(lhs(RE));(% o3) n ( n + 4) ( 64 n + 2240 n + 35056 n + 323344 n + 1949788 n + 8053956 n +23188049 n +46338535 n +62583534 n +53821965 n +26011175 n +5133154 ) a n . For simplicity of explanation, we will present computations over the rationals. The caseof extension fields works similarly, this choice is just to avoid lengthy notations for rootslabeling.For the leading polynomial coefficient, the monic factors to be considered after factor-ization in Q modulo Z with roots real parts in I are ( n + 1) e (cid:18) n + 12 (cid:19) e , for e , e . For the trailing term we have ( n + 1) e for e . Therefore ratios of Pochhammer parts of hypergeometric term solutions are among thefollowing , n + 1) , n + 1) , n + 1) , (cid:0) n + (cid:1) , (cid:0) n + (cid:1) , ( n + 1)( n + 1) (cid:0) n + (cid:1) , ( n + 1) (cid:0) n + (cid:1) , ( n + 1) (cid:0) n + (cid:1) , ( n + 1) (cid:0) n + (cid:1) (4.10)Observe that none of these ratios has a non-negative integer zero or pole, hence the typeof rational function that we treat with Algorithm 2. This is made possible by the fact thatwe consider factorization modulo Z with roots real parts in I .Moreover, not all ratios in (4 . should be considered because the exponents of eachlinear factor appearing in the possible ratios of hypergeometric term solutions can bebounded from the given holonomic recurrence equation. For this purpose van Hoeij’s al-gorithm uses the notion of valuation growth or local types of difference operators at finite VARIANT OF VAN HOEIJ’S ALGORITHM TO COMPUTE HYPERGEOMETRIC TERMS 13 singularities [18, Definition 9]. Such a point is simply a root modulo Z of the trailing orthe leading polynomial coefficient of (4 . as we considered.Since we are already computing ratios of Pochhammer parts of hypergeometric termsolutions, we proceed in a slightly different way than what is described in ([18, 5]) tocompute exponent bounds at finite singularities. We observed that this reduces to takeminimum exponents (or valuations) taken by the corresponding factors modulo Z in thetrailing and the leading polynomial coefficients of the initial recurrence equation as lowerbounds. Upper bounds are automatically found while computing ratios. Coming backto our example, it follows that ratios with denominator ( n + 1 / should be removed.Therefore the remaining ratios are , n + 1) , n + 1) , n + 1) , (cid:0) n + (cid:1) , ( n + 1)( n + 1) (cid:0) n + (cid:1) , ( n + 1) (cid:0) n + (cid:1) . (4.11)Note that these considerations on monic factors of leading and trailing polynomial coef-ficients can already be seen as an important efficiency gain when comparing computationswith those in [12] which have to consider more cases. We are going to see more tools inthe remaining part of the algorithm that will also filter the ratios in (4 . .4.2. Local type at infinity.
Without ambiguity, we will more often use the terminology"local type" instead of "local type at infinity" since we only consider computations at infin-ity. This is about a characteristic property of hypergeometric term solutions of holonomicrecurrence equations.We study the behavior of a hypergeometric term ( a n ) ratio r ( n ) at infinity. Indeed, at ∞ we can write(4.12) r ( n ) = c · n ν · (cid:18) bn + O (cid:18) n (cid:19)(cid:19) , with the unique triple ( ν, c, b ) called the local type of a n at ∞ . Theorem 4.2 (Fuchs Relations) . Let R ( n ) = N ( n ) U ( n ) with N ( n ) , U ( n ) ∈ K [ n ] . The follow-ing relations between the local type of a hypergeometric term a n given by (4 . hold: i. ν = P Kk =1 e k , ii. b = P Kk =1 θ k e k + deg( N ( n )) − deg( U ( n )) , iii. c = C, where ( ν, c, b ) denotes the local type of a n at ∞ . Proof.
From (4 . we know that(4.13) r ( n ) = a n +1 a n = C · R ( n + 1) R ( n ) · K Y k =1 ( n + θ k ) e k ! . We would like to compute a truncated asymptotic expansion of (4 . . This can be seenas the result of the product of asymptotic expansions of the form (4 . of R ( n +1) R ( n ) and Q Kk =1 ( n + θ k ) e k times C . Since R ( n ) = N ( n ) U ( n ) , the highest degree of n in its asymptoticexpansion is δ = deg( N ( n )) − deg( U ( n )) . Hence we read as(4.14) R ( n ) = c R · n δ · (cid:18) b R n + O (cid:18) n (cid:19)(cid:19) , for some constants c R , b R . Let us now deduce a truncated asymptotic expansion of R ( n +1) . R ( n + 1) = c R · ( n + 1) δ · (cid:18) b R n + 1 + O (cid:18) n (cid:19)(cid:19) = c R · n δ (cid:18) n (cid:19) δ · b R n (cid:18) n (cid:19) + O (cid:18) n (cid:19) = c R · n δ δn + δ X j =2 (cid:18) δj (cid:19) (cid:18) n (cid:19) j · (cid:18) b R n + O (cid:18) n (cid:19)(cid:19) = c R · n δ (cid:18) b R + δn + O (cid:18) n (cid:19)(cid:19) . (4.15)Thus from (4 . and (4 . the first order asymptotic expansion of R ( n +1) R ( n ) yields R ( n + 1) R ( n ) = 1 + b R + δn + O (cid:18) n (cid:19) b R n + O (cid:18) n (cid:19) = (cid:18) b R + δn + O (cid:18) n (cid:19)(cid:19) · (cid:18) − b R n + O (cid:18) n (cid:19)(cid:19) = 1 + δn + O (cid:18) n (cid:19) = 1 + deg( N ( n )) − deg( U ( n )) n + O (cid:18) n (cid:19) . (4.16)On the other hand ( n + θ k ) e k = n e k · (cid:18) θ k n (cid:19) e k = n e k · θ k e k n + e k X j =2 (cid:18) e k j (cid:19) (cid:18) θ k n (cid:19) j = n e k · (cid:18) θ k e k n + O (cid:18) n (cid:19)(cid:19) , (4.17)therefore(4.18) K Y k =1 ( n + θ k ) e k = n P Kk =1 e k · P Kk =1 θ k e k n + O (cid:18) n (cid:19)! . VARIANT OF VAN HOEIJ’S ALGORITHM TO COMPUTE HYPERGEOMETRIC TERMS 15
Finally according to (4 . , the expansion sought is obtained by the product of (4 . and (4 . times C . That is r ( n ) = C · n P Kk =1 e k · P Kk =1 θ k e k n + O (cid:18) n (cid:19)! · (cid:18) N ( n )) − deg( U ( n )) n + O (cid:18) n (cid:19)(cid:19) = C · n P Kk =1 e k P Kk =1 θ k e k + deg( N ( n )) − deg( U ( n )) n + O (cid:18) n (cid:19)! , (4.19)from which one easily read off the data of the theorem. (cid:3) The first two relations in this theorem tell us that for the local type ( ν, c, b ) of a hyperge-ometric term a n , ν and b can be found directly from a ratio representing its Pochhammerpart. Indeed, observe that modulo Z , the second relation of the theorem reads as(4.20) b = K X k =1 θ k e k . The third relation will be considered later in this subsection. It is straightforward to find ν and b corresponding to a hypergeometric term local type from the ratio of its Pochham-mer part. Thus this constitute the next step after obtaining ratios as in (4 . . For example ( n + 1) / ( n + 1 / = n − (1 − / (4 n ) + O (1 /n )) and therefore ν = − and b = − (modulo Z ). As mentioned earlier, the map y y − ⌊ y ⌋ − is used to find the correspon-dence of y modulo Z in I .Next, we explain how the local types of hypergeometric term solutions of (4 . arecomputed. This step is considered with the highest priority in our algorithm, because ifthe set of local types of hypergeometric term solutions of a given holonomic recurrenceequation is empty, then there is no hypergeometric term solution over the considered field.For this step, van Hoeij’s algorithm uses the Newton polygon of the difference operator(see [18, Section 3]). However, we proceed differently. Our idea is to rewrite (4 . forratios of hypergeometric term solutions, substitute (4 . inside, and compute the asymp-totic expansion of the nonzero side to find equations for the local types by equating theresult to . This process is the same Petkovšek used to develop its algorithm Poly (see [12,Algorithm Poly]).Let a n be a hypergeometric term solution of this equation such that a n +1 = r ( n ) a n fora rational function r . (4 . can then be written for r ( n ) as(4.21) d X i =0 P i i − Y j =0 r ( n + i ) = 0 . We assume(4.22) r ( n ) = c · n ν · (cid:18) O (cid:18) n (cid:19)(cid:19) and we substitute this in (4 . . Similarly as we did in the proof of Theorem . , wemake computations that yield the possible values of ν and c . If such values are found, say ( ν cand , c cand ) , then we rewrite r ( n ) as(4.23) c cand · n ν cand · (cid:18) bn + O (cid:18) n (cid:19)(cid:19) and we make new computations to find b .Summarized, our procedure to find local types ( ν, c, b ) of hypergeometric term solutionsof (4 . consists in the following items:(1) we compute the possible values for ν ;(2) for each value of ν ,2-a we compute possible values for c ,2-b for each value found for c , we use ν and c to compute the possible values for b ;2-c for each value found for b , ( ν, c, b ) constitutes a local type of a hypergeomet-ric term solution of (4 . .Let us now explain how each value is computed. • Computing ν :Substitute (4 . in (4 . gives the following terms on the left-hand side(4.24) c i · n i · ν · P i · (cid:18) O (cid:18) n (cid:19)(cid:19) , (0 i d ) which is equivalent to(4.25) l i · c i · n i · ν +deg( P i ) · (cid:18) O (cid:18) n (cid:19)(cid:19) , (0 i d ) where l i denotes the leading coefficient of P i . Since we are dealing with an equal-ity with right-hand side , the terms having the highest power of n in the asymp-totic expansion of the equation left-hand must be zero. However this is only pos-sible if a term of the form (4 . has the same power of n with some other termsso that they add to . Therefore we deduce that possible candidates for ν are in-teger solutions of linear equations coming from equalities of powers of n for twodifferent terms of the form (4 . . That is for i = j d , we have the equation(4.26) i · ν + deg( P i ) = j · ν + deg( P j ) and therefore a possible value for ν is(4.27) ν i,j = deg( P j ) − deg( P i ) i − j , if the computed value is an integer.We then compute (cid:0) d (cid:1) such values for (4 . and keep the integers. Note thattwo different couples of terms may give the same value for ν , meaning that thecorresponding addition to zero involves all the underlying terms, which is the pointof the next item. • Computing c :Assume that we have found a value ν i,j ∈ Z corresponding to k terms in theequation (4 . with indices u = u = . . . = u k d . Then from (4 . one easily see that a candidate for c is a solution of the polynomial equation(4.28) l u · c u + l u · c u + · · · + l u k · c u k = 0 . In fact, since the corresponding terms must add to zero in the asymptotic expan-sion, their leading coefficients must equal zero. Note that (4 . is solved over theconsidered field K .Thus,for each value c i,j ∈ K which is a zero of (4 . for a given ν i,j , ( ν i,j , c i,j ) is already a possible couple to be completed for the local type of a hypergeometricterm solution of (4 . . VARIANT OF VAN HOEIJ’S ALGORITHM TO COMPUTE HYPERGEOMETRIC TERMS 17 • Computing b :For a computed couple ( ν i,j , c i,j ) as explained above, we rewrite r ( n ) as(4.29) c i,j · n ν i,j · (cid:18) bn + O (cid:18) n (cid:19)(cid:19) , with unknown b .After substituting (4 . in (4 . and computing again the asymptotic expan-sion, terms with highest powers of n add to zero, and therefore the left-hand sideof the resulting equation must have a leading term with coefficient as a polynomialin the variable b . Since that polynomial must be zero, the possible values for b areits roots. This can be done by computing asymptotic expansion and solve the co-efficients equal to zero for the unknown b . Finally if we find values for b ∈ K thenwe have found for each b a local type ( ν, c, b ) of a possible hypergeometric termsolution of (4 . over K .Hence we get the following algorithm. Theorem 4.3.
Algorithm finds all the local types ( ν, c, b ) of all hypergeometric termsolutions of (4 . .Remark . • When extension fields are allowed, Algorithm is used to bound the degree ofsuch extensions. Indeed, the computation of c and b for local types defines thedegree bound of extension fields to carry throughout the remaining part of the al-gorithm. This is sometimes important to avoid computations with algebraic num-bers whose degrees are larger. This is basically how our algorithm avoid splittingfields from the polynomials of degree in the trailing and leading coefficientsof ( RE ) . More theoretical details on dealing with algebraic extensions can beseen in [5, Section 8]. Regarding implementations, this step is better handled withMaple than Maxima: by its command RootOf , Maple is able to manage compu-tations with non-explicit values of algebraic numbers; in Maxima however, sucha command is not yet available and so solutions over some algebraic extensionfields might be missed with our implementation. Nevertheless, this is not to pre-tend that our Maple implementation could work on every extension fields, becausesome computations remains not easily handled. We just want to point out imple-mentation limits that are well considered theoretically. This might be the reasonwhy Maple current
LREtools[hypergeomsols] seems to be an algebraic-numbersbased implementation; we will present examples where
LREtools[hypergeomsols] gets correct results only when polynomial coefficients are converted to symbolicalgebraic numbers using
RootOf . • A further notice about extension fields is a property similar to the conjugate roottheorem (see [5, Lemma 3]). Indeed, when a local type ( ν, c, b ) where c is definedover an extension field leads to a basis of hypergeometric term solutions, say B c ,then local types corresponding to conjugates of c lead to bases of same dimensionsas B c that, only differ from B c by conjugations of c . Therefore an importantefficiency can be gained by using this property. At the current time, this is notconsidered in our implementations but will certainly be the case in future releases. • We mention that computations of Algorithm can sometimes be used to reducethe number of iterations in the implementation. The important point to notice isthat when two linearly independent hypergeometric term solutions have the samelocal type, Algorithm computes it at least twice. Therefore collecting local types Algorithm 3
Compute local types of all hypergeometric term solutions of (4 . Input:
Polynomials P i ( n ) ∈ K [ n ] , i = 0 , . . . , d | P d ( n ) · P ( n ) = 0 Output:
The set of all local types of hypergeometric term solutions of the holonomic RE d X i =0 P i ( n ) a n + i = 0 . (1) Set L = {} .(2) For all pairs { i, j } ∈ { , , . . . , d } , compute(4.30) ν i,j = deg( P j ) − deg( P i ) i − j . (3) For each integer ν i,j computed in (4 . , compute the set of solutions in K , say S c,i,j , of the polynomial equation(4.31) l u · c u + l u · c u + · · · + l u j · c u k = 0 , where l u , l u , . . . , l u k are the leading coefficients of the polynomials P u , P u , . . . , P u k , u = u = . . . = u k d satisfying (4 . for thesame integer ν i,j .(a) For each element c i,j of S c,i,j set(4.32) r ( n ) = c i,j · n ν i,j · (cid:18) bn (cid:19) . (4) (b) Compute the coefficient T i,j ( b ) of the first non-zero term of the asymptoticexpansion of(4.33) d X i =0 P i i − Y j =0 r ( n + i ) . (c) Solve T i,j ( b ) = 0 in K for the unknown b and define S b,i,j to be the set ofsolutions.(d) For each element b i,j ∈ S b,i,j , add the triple ( ν i,j , c i,j , b ) to L .(5) Return L .in a list might be advantageous, so that when a basis of hypergeometric termscorresponding to a particular local type is found, the latter is discarded from thelist of local types. This is a useful tool when the number of computed local types(with repeated values) is less than the order of the given holonomic recurrenceequation. Otherwise there is no reason to save local types in a list.Thus any ratio candidates whose local type is not in the list of local types (deprivedof values for c ) should not be used in further steps. We implemented a Maxima function localtype(L,n) which takes the polynomial coefficients of a holonomic recurrence equationin L in the variable n and returns a list of triples [ ν, c, b ] . Applying it for ( RE ) yields (% i1) localtype(expand(REcoeff(RE,a[n])),n);(% o1) [[ − , − , − , [ − , , − , [0 , − , − , [0 , , − . REcoeff is our code to collect coefficients. These are expanded using the Maxima com-mand expand . From the obtained output it follows that / ( n + 1) and ( n + 1) should VARIANT OF VAN HOEIJ’S ALGORITHM TO COMPUTE HYPERGEOMETRIC TERMS 19 also be removed from potential ratios of hypergeometric term solution Pochhammer partsin (4 . . Note, however, that it is from the computed local types that we get the possiblevalues for C in (4 . according to the third relation in Theorem 4.2. These will be used inthe next step together with their corresponding Pochhammer part ratios.4.3. Rational part of hypergeometric terms.
The algorithm goes further in filtering theset of Pochhammer part ratios. Indeed, once we have found all those better candidatesfor ratios of hypergeometric term solution Pochhammer parts, we need to use again thesecond Fuchs relation from Theorem 4.2 in order to find δ = deg( N ( n )) − deg( U ( n )) ,where N ( n ) and U ( n ) are the numerator and the denominator of R in (4 . . In fact, sincewe have found values for b and its possible ratio candidates, which means that we cancompute P Kk =1 θ k · e k , we therefore deduce that these candidates are valid if and only ifthey satisfy(4.34) δ = b − K X k =1 θ k · e k ∈ Z . In this case the verification of ratios of Pochhammer parts of hypergeometric term solutionsfor the value of b should consist in checking if the difference b − P Kk =1 θ k · e k is an integer.However, (4 . can be used in the algorithm only if b is not computed modulo Z .Another approach is to use again asymptotic expansion. Since now we have the ratios withtheir corresponding values of c , according to (4 . we can write(4.35) r ( n ) = R ( n + 1) R ( n ) · c · h ( n + 1) h ( n ) . Moreover(4.36) R ( n + 1) R ( n ) = 1 + δn + O (cid:18) n (cid:19) , where δ is as in (4 . . Thus the asymptotic expansion of (4 . (left-hand side) with r ( n ) used by combining (4 . and (4 . must have a leading term as a polynomial co-efficient in the variable δ . Therefore values of δ are integer roots (if there are some) ofthat polynomial. If there are not such roots, then the rational function c · h ( n + 1) /h ( n ) isremoved from the potential hypergeometric term solution Pochhammer parts of (4 . . Ourimplementation uses this second approach.Mostly after this step the number of Pochhammer part ratios of hypergeometric termsolution of (4 . is considerably reduced or equal to the exact number of hypergeometricterm solutions.Now, it only remains to find the rational function R in (4 . whose a holonomic recur-rence equation can be easily computed. Let c · h ( n + 1) /h ( n ) be one of the remainingratios times its corresponding c for the local type. Then the recurrence equation(4.37) d X i =0 P i · R ( n + i ) · c n + i · h ( n + 1 + i ) h ( n + i ) = 0 , is an equation for the unknown rational function R ( n ) that we can easily modify to aholonomic recurrence equation. Note, however, that there is no need to use a complete al-gorithm for computing rational solutions of holonomic recurrence equations. Indeed, sincewe already have the difference between the degrees of the numerators and the denominatorsof rational solutions of (4 . , it is enough to use an algorithm that computes a universal denominator U ( n ) of all rational solutions of (4 . , and use δ or its maximum value (forthe second approach we proposed) (see (4 . ) to compute a degree bound δ + deg( U ( n )) for the degrees of the corresponding numerators. Substituting N ( n ) /U ( n ) in (4 . where N ( n ) is an arbitrary polynomial of degree δ + deg( U ( n )) results in a linear system in thecoefficients of the arbitrary polynomial N ( n ) . Finally solving that system gives a basis ofall the rational functions R ( n ) = N ( n ) /U ( n ) sought.Let us then say a few words on the computation of a universal denominator of rationalsolutions of holonomic recurrence equations. Abramov has proposed most key results forthat purpose (see [2, 1, 3]). A crucial step in Abramov’s original algorithm is to computethe dispersion set of two polynomials . The dispersion set can efficiently be obtained fromfull factorization as described in [11] (see also [10, Algorithm 5.2]). We use this method inour implementation of Abramov’s algorithm to compute universal denominators. As a littlestory, note that we could not find neither in Maxima, nor in Maple a satisfactory (in termsof efficiency) implementation for computing dispersion sets. Therefore we implementedthe algorithm in [11] and added it as a by product of our Maxima and Maple packages.We think this last step of computing the rational function R ( n ) might sometimes makea difference of efficiency between our algorithm and van Hoeij’s original version. In hisapproach, van Hoeij uses a special algorithm from his idea of finite singularities to deter-mine R ( n ) (see [17]). Though we mentioned that a complete algorithm for that purposeis not necessary, the algorithm in [17] is sometimes suitable with the computations in [18].However, comparisons in [3] shows that using our approach or van Hoeij’s one at this stepdoes guarantee efficiency gain of one over the other.4.4. Our algorithm.
We can now present the complete algorithm of this paper.
Algorithm 4
Compute hypergeometric term solutions of (4 . Input:
Polynomials P i ( n ) ∈ K ( n ) , i = 0 , . . . , d, | P d ( n ) · P ( n ) = 0 . Output:
A basis of all hypergeometric term solutions of the holonomic recurrence equa-tion(4.38) d X i =0 P i ( n ) a n + i = 0 over K .(1) Set H = {} .(2) Use Algorithm to compute the set L of all local types at infinity of hypergeo-metric term solutions of (4 . .(3) If L = ∅ , then stop and return H . A universal denominator of rational solutions of a holonomic recurrence equation is a polynomial that isdivisible by all the denominators of rational solutions of that holonomic equation [1]. The dispersion set of A ( n ) and B ( n ) can be defined as the set of all non-negative integer roots of theresultant polynomial of A ( n ) and B ( n + h ) in the variable h . VARIANT OF VAN HOEIJ’S ALGORITHM TO COMPUTE HYPERGEOMETRIC TERMS 21
Algorithm 4
Compute hypergeometric term solutions of (4 . (4) Construct the set of couple (numerator, denominator)(4.39) P := (cid:26) ( p ( n ) , q ( n )) ∈ K [ n ] : p ( n ) and q ( n ) are monic factors modulo Z with roots real parts in [ − , of P ( n − and P d ( n − d ) respectively (cid:27) , for ratio candidates of hypergeometric term solution Pochhammer parts.(5) Remove from P all couple whose p ( n ) exponents are less than the minimum mul-tiplicity of the corresponding root modulo Z in the trailing polynomial coefficient P ( n ) . Similarly, clear P by the same consideration for q ( n ) exponents and theleading polynomial coefficient P d ( n ) . Finally substitute each remaining couple ( p ( n ) , q ( n )) in P by p ( n ) q ( n ) .(6) Remove elements in P that have a larger algebraic degree than the bound givenby the local types in L .(7) Construct the set F of c · r , r ∈ P such that c · r has its local type at infinity asan element of L .(4.40) F := (cid:26) c · r : r = n ν r (cid:18) b r n + O (cid:18) n (cid:19)(cid:19) ∈ P and ( ν r , c, b n ) ∈ L (cid:27) . (8) Set F := {} . For each element f ( n ) of F (a) Compute a recurrence equation, say E f with the coefficients(4.41) P i · i Y j =0 f ( n + i ) , i = 0 , . . . , d, for the rational function R ( n ) in (4 . of the possible hypergeometric termsolutions.(b) Substitute the terms R ( n + i ) by (1+ δn + i ) , i = 0 , . . . , d, in E f and computethe coefficient of the leading term of the asymptotic expansion of the lefthand side of E f , say Q f ( δ ) .(c) Compute the set S δ f of integer roots of Q f ( δ ) .(d) If S δ f = ∅ then f ( n ) is discarded.(e) Else set δ f := max( S f ) , rewrite E f in a holonomic form and add ( f ( n ) , δ f , E f ) in F .(9) If F = ∅ then stop and return H .(10) For each ( f ( n ) , δ f , E f ) ∈ F (a) Compute the universal denominator U f ( n ) of rational solutions of E f byusing the approach in [11] to find the needed dispersion set.(b) Update E f as E ′ f with U f ( n ) to get a holonomic recurrence for numeratorsof rational solutions of E f .(c) Set d N f := deg( U f ( n )) + δ f , and find a basis of all polynomial solutionsof degree at most d N f of E ′ f .(d) Use Algorithm to compute h f ( n ) = Q n − k =0 f ( k ) .(e) For each N f ( n ) ∈ S N f add N f ( n ) U f ( n ) · h f ( n ) to H .(11) Return H We implemented Algorithm in Maxima as HypervanHoeij(RE,a[n],[K]) , with the de-fault value Q (for rationals ) for K representing the field where solutions are computed.One must specify C for K to allow computations over extension fields of Q . Using thisimplementation to solve ( RE ) yields (% i1) HypervanHoeij(RE,a[n]);Evaluation took . seconds (0 . elapsed ) using . MB . (% o1) ( ( n + 1) ( n + 2) ( n + 3) , ( − n n , n ! , ( − n n n ! (2 n ) ! ) , with timing (allowed with the Maxima command showtime ) . second and . MB memory used. In Maple, we implemented our algorithm as rectohyperterm in ourpackage
FPS . Le us do the same computation with our Maple implementation and LRE-tools[hypergeomsols] . > RE:=FPS[sumhyperRE]([binomial(n+3,n),1/n!, > (-1)^n/n,(-1)^n/pochhammer(1/2,n)^2],a(n)): > Usage(FPS[rectohyperterm](RE,a(n))) memory used=15.94MiB, alloc change=4.00MiB, cpu time=141.00ms, real time=135.00ms,gc time=0ns ( ( n !) − , ( − n n , ( n + 3) ( n + 2) ( n + 1) , ( − n ( n !) n ((2 n )!) ) > Usage(LREtools[hypergeomsols](RE,a(n),{},output=basis)) memory used=21.48MiB, alloc change=25.99MiB, cpu time=343.00ms, real time=277.00ms,gc time=125.00ms " n + 6 n + 11 n + 6 , ( − n n , (Γ ( n + 1)) − , ( − n (Γ (1 / n )) The
Usage command from the
CodeTools package is used to display timings and mem-ory used. Here one can see the advantage of having an implementation that uses extensionfields from user specifications. Allowing extension fields for this example unnecessarilyincreases the timing (will be the same as for
LREtools[hypergeomsols] ) of computationssince ( RE ) is of order and we already have hypergeometric term solutions in the out-put. 5. S OME COMPARISONS
Our Maple implementation of the given algorithm was tested on many recurrence equa-tions. Regarding efficiency, the difference between our implementation and the internalMaple 2020
LREtools[hypergeomsols] is in the order of milliseconds: for solutions overthe rationals, our implementation generally gives a better efficiency; and for solutions overextension fields
LREtools[hypergeomsols] is generally faster, although the timing are veryclosed. Nevertheless, there are examples where this comparison does not hold but the tim-ings remains closer. Therefore it is more suitable to say that both implementations have Rationally valued, not symbolically rational: a parameter declared as rational is not considered as such inthe implementation. Denoting formal power series (FPS), an implementation of main algorithms from [16].
VARIANT OF VAN HOEIJ’S ALGORITHM TO COMPUTE HYPERGEOMETRIC TERMS 23 same or similar efficiency, and that it may differ depending on internal commands used. Forexample, one particular internal computations that slows down our Maple implementationin certain cases is the asymptotic expansion of Algorithm 3, for which Maple commands series or asympt does not always yield an expansion of the required order. For this reasonwe repeat the computation (of course, by increasing the order) until the desired order isobtained. Furthermore, to facilitate the steps which follow, the coefficients obtained mustoften be simplified since at times they are equivalent to zero. This issue does not occur withour Maxima implementation although it usually comes third when we compare efficiencies.But again, there are some examples where this is not verified, sometimes our Maxima codeis the fastest; indeed, we think for comparisons involving implementations in Maple andMaxima or two different CAS in general, a first look must be taken at kernels and datastructures of both systems. Without going into details on this computer science ’buildings’we only mention that Maple has a C-base kernel whereas Maxima has a Common Lisp-base one, and this could make some differences in the speed of both systems. The readercan visit the programming website https://open.kattis.com to see how often Cprograms are the fastest. The author particularly recommends to check the 0-1 Sequencesproblem.On the other hand, we have been able to find bugs from Maple 2020 internal command.We give here three unexpected examples where our implementation finds results highlight-ing the issues encountered. > RE1:=FPS[sumhyperRE]([(1-sqrt(7))^n/n!, > (1+sqrt(7))^n/(n+1)!],a(n)) RE1 := (cid:16) √ − n − (cid:17) a ( n ) + 4 ( n + 2) (cid:16) − n + 4 √ − (cid:17) a ( n + 1) − ( n + 2) ( n + 3) (cid:16) √ − n − (cid:17) a ( n + 2) = 0 > LREtools[hypergeomsols](RE1,a(n),{},output=basis) which means no solution found; > RE2:=FPS[sumhyperRE]([3^(n/5),3^(n/2)],a(n))
RE2 := 3 a ( n ) + (cid:16) − / − √ (cid:17) a ( n + 1) + 3 / a ( n + 2) = 0 > LREtools[hypergeomsols](RE2,a(n),{},output=basis)
Error, (in mod/Primfield/ReduceField/sort_B) too many levels of recursion.However, the above issues can be avoided by rewriting algebraic numbers in their Maplestandard representation using convert/RootOf . This is a temporary workaround proposedby Mark van Heoij, the bug will be corrected in future Maple releases. Hence the reasonwhy we mentioned that
LREtools[hypergeomsols] is an algebraic-number based implemen-tation. Using our implementation, one gets the expected results > FPS[rectohyperterm](RE1,a(n),C) ( (cid:0) − √ (cid:1) n n ! , (cid:0) √ (cid:1) n ( n + 1) n ! ) > FPS[rectohyperterm](RE2,a(n),C) n(cid:16) √ (cid:17) n , (cid:16) √ (cid:17) n o . It is important that the hypergeometric property is kept in implementations so that com-putations extend to symbolic use of expressions that may not represent numbers. The nextexample type might be of interest for sequence of functions. > RE3:=FPS[sumhyperRE]([ln(x)^n,ln(x*y)^n],a(n))
RE3 := a ( n ) ln ( x ) ln ( xy ) + ( − ln ( x ) − ln ( xy )) a ( n + 1) + a ( n + 2) = 0 > LREtools[hypergeomsols](RE3,a(n),{},output=basis)
Error, (in mod/Normal/Factored) not implemented > FPS[rectohyperterm](RE3,a(n),C) { (ln ( x )) n , (ln ( xy )) n } As with the two first examples, we believe the latter bug can be fixed as well. Thepoint here is to show how having our algorithm, a variant of van Hoeij’s algorithm, con-tributes to ensure and present equivalences of theoretical arguments in [12, 18, 5], and im-prove cutting edge implementations. Although van Hoeij’s algorithm is mentioned at thisfootnote link to a Mathematica webpage, the implementation in Rsolve remains quiteslow: we computed the solutions of ( RE ) and got a much complicated result after about minutes of computations. It is difficult to decide which algorithm is used as the codeis hidden from users. Petkovšek’s algorithm is the most popular implementation encoun-tered in many CASs. Our result and its implementation bring Maxima to the top level incomputing hypergeometric term solutions of holonomic recurrence equations. A softwaredemonstration was recently presented at the International Congress of Mathematical Soft-ware (ICMS) 2020. The algorithm of this paper can be used as an important black box forthe general case of m -interlacings of hypergeometric sequences or m -fold hypergeometricterms, m ∈ N (see [16, 13]). R EFERENCES1. Sergei A Abramov,
Rational solutions of linear difference and q-difference equations with polynomial coeffi-cients , Program. Comput. Softw. (1999), no. 6, 369–374.2. Sergei A Abramov and Moulay A Barkatou, Rational solutions of first order linear difference systems , ISSAC,vol. 98, Editors: Weispfenning, Volker and Trager, Barry, Association for Computing Machinery, New York,1998, pp. 124–131.3. Sergei A Abramov, Amel Gheffar, and DE Khmelnov,
Rational solutions of linear difference equations:Universal denominators and denominator bounds , Program. Comput. Software (2011), no. 2, 78–86.4. Moulay Abdelfattah Barkatou and Mathieu Duval, Sur les séries formelles solutions d’équations aux dif-férences polynomiales , Ann. Inst. Fourier. (1994), no. 2, 495–524.5. Thomas Cluzeau and Mark van Hoeij, Computing hypergeometric solutions of linear recurrence equations ,Appl. Algebra Engrg. Comm. Comput. (2006), no. 2, 83–115.6. Mathieu Duval, Lemmes de Hensel et factorisation formelle pour les opérateurs aux différences , Funkcial.Ekvac. (1996), no. 3, 349–368.7. Keith O. Geddes, Stephen R Czapor, and George Labahn, Algorithms for computer algebra , Kluwer Aca-demic Publishers, Massachusetts, 1992.8. Wolfram Koepf,
Symbolic computation of formal power series with Macsyma , Functional Analytic MethodsIn Complex Analysis And Applications To Partial Differential Equations, Editors: Tutschke, Wolfgang andMshimba, Ali, World Scientific, Singapore, New Jersey, London, Hong Kong, 1995, pp. 306–328. https://reference.wolfram.com/language/tutorial/SomeNotesOnInternalImplementation.html VARIANT OF VAN HOEIJ’S ALGORITHM TO COMPUTE HYPERGEOMETRIC TERMS 25
9. ,
Computeralgebra: eine algorithmisch orientierte Einführung , Springer-Verlag, Berlin Heidelberg,New York, 2006.10. ,
Hypergeometric summation, an algorithmic approach to summation and special function identities ,2nd ed., Springer-Verlag, London, 2014.11. Yiu-Kwong Man and Francis J Wright,
Fast polynomial dispersion computation and its application to indef-inite summation , ISSAC, Editor: MacCallum, Malcolm, Association for Computing Machinery, New York,1994, pp. 175–180.12. Marko Petkovšek,
Hypergeometric solutions of linear recurrences with polynomial coefficients , J. Symb.Comput. (1992), no. 2-3, 243–264.13. Bertrand Teguia Tabuguia and Wolfram Koepf Power series representations of hypergeometric types and non-holonomic functions in computer algebra , in preparation to be submitted in J. Symb. Comput.14. Bruno Salvy and Paul Zimmermann,
Gfun: a maple package for the manipulation of generating and holo-nomic functions in one variable , ACM Trans. Math. Softw. (1994), 163–177.15. Richard P Stanley, Differentiably finite power series , European J. Combin. (1980), no. 2, 175–188.16. Bertrand Teguia Tabuguia, Power series representations of hypergeometric types andnon-holonomic functions in computer algebra , Ph.D. thesis, University of Kassel, https://kobra.uni-kassel.de/handle/123456789/11598 , 2020.17. Mark Van Hoeij,
Rational solutions of linear difference equations , ISSAC, Editors: Weispfenning, Volkerand Trager, Barry, Association for Computing Machinery, New York, 1998, pp. 120–123.18. ,
Finite singularities and hypergeometric solutions of linear recurrence equations , J. Pure Appl. Al-gebra (1999), no. 1-3, 109–131.D
EPARTMENT OF M ATHEMATICS , U
NIVERSITY OF K ASSEL , H
EINRICH -P LETT -S TR .40, 34132 K ASSEL ,G ERMANY
Email address ::