(* Mixed Gosper's algorithm for Mma 3.0, by Marko Petkovsek *) (* Mixed Poly for Mma 3.0, by Marko Petkovsek *) (* October 6, 1996 *) (* October 23, 1997: corr. treatment of Product in FactorialSimplify[expr] *) (* October 30, 1997: impr. treatment of Product in FactorialSimplify[expr] *) (* May 18, 1998: incl. treatment of qFac in FactorialSimplify[expr] *) (* May 29, 1998: fixed bugs found by W. Gosper (added PowerExpand and Together to ToY *) BeginPackage["GosperSum`"]; Off[General::spell1]; FactorialSimplify::usage = "FactorialSimplify[expr] is a simple-minded simplification routine for expressions containing rational functions, Factorial, Factorial2, Binomial, Gamma, Pochhammer, and Product. The result is expressed by means of Binomial and Factorial functions. FactorialSimplify[expr, n] simplifies linear combinations of terms which are hypergeometric in n. Alias: FS."; FS::usage = "FS is alias for FactorialSimplify."; GetRatio::usage = "GetRatio[h, n] returns the rational consecutive-term ratio h[n+1]/h[n] for a hypergeometric term h[n]."; GosperFunction::usage = "GosperFunction[r, n, qList:{}] uses Gosper's algorithm to find for a given rational function r[n, q1^n, .., qk^n] another rational function R[n, q1^n, .., qk^n] such that r NR - R == 1. If no such function exists then GosperFunction returns $Failed."; GosperSum::usage = "GosperSum[term, {var, n0, n1}, qList:{}] finds a hypergeometric closed form for Sum[term, {var, n0, n1}] if possible. GosperSum[term, var, qList:{}] is the corresponding indefinite sum. Term can contain rational functions, Factorial, Factorial2, Binomial, Gamma, Pochhammer, Product, and qFac. GosperSum uses the generalization of Gosper's algorithm to the mixed multibasic case by A.Bauer and M.Petkovsek."; MixedNormalForm::usage = "MixedNormalForm[r, n, y] returns the normal form {a, b, c} for the rational function r[n, y[q1], .., y[qk]]."; MixedPoly::usage = "MixedPoly[eqn, x[n], qList, paramList:{}] finds polynomial solutions of an inhomogeneous mixed recurrence with parameters."; Zeil::usage = "Zeil[F, s[n], {k, lo, hi}] returns a recurrence rec of the form Ls[n] == rhs[n] satisfied by the sum s[n] = Sum[F[n,k], {k, lo, hi}]. Here s, n and k must be symbols. Options of Zeil: Rhs -> False suppresses the right hand side of the recurrence. Certificate -> True returns {R, rec} where R is a rational function such that G = RF satisfies LF[n,k] == G[n,k+1] - G[n,k]. Bases -> {q, ... } specifies the bases in case of a (multi)basic summand. Start -> d starts the search for L at order d rather than 0." Start::usage = "An option for Zeil specifying the starting value for the order of the sought recurrence. Default value: 0."; Bases::usage = "An option for Zeil specifying the list of bases in the (multi)basic case. Default value: {}."; Rhs::usage = "An option for Zeil specifying whether or not the right-hand side of the recurrence is required. Default value: True."; Certificate::usage = " An option of Zeil specifying whether or not a certificate of the recurrence is required. Default value: False."; qFac::usage = "The q-factorial function: qFac[a,q,n] = Product[1 - a q^i, {i,0,n-1}]."; WZ::usage = "WZ[F, n, k, qList:{}] returns a rational function R(n,k) such that F(n,k) and G(n,k) := R(n,k) F(n,k) are a WZ pair: F(n+1,k) - F(n,k) = G(n,k+1) - G(n,k)." ; ToY::usage = "ToY[expr, n, qList, y] replaces all appearances of q^n in expr with y[q], for all q in qList." ; Begin["Private`"]; Share[]; Unprotect[Power]; 0^0 = 1; Protect[Power]; (* WZ METHOD *) WZ[F_, n_, k_, qList_:{}] := Module[{rn, rk, rr}, rn = GetRatio[F, n]; rk = GetRatio[F, k]; rr = Together[rn - 1]; out = GosperFunction[ Factor[rk (rr /. k->k+1) / rr], k, qList]; If[out === $Failed, Return[$Failed], Return[Together[rr out]] ] ]; (* FACTORIAL SIMPLIFY *) rf[a_, 0] = 1; rf[a_, n_Integer?Positive] := Product[a+i, {i,0,n-1}]; rf[a_, n_Integer?Negative] := 1 / Product[a-i, {i,1,-n}]; FS = FactorialSimplify; FactorialSimplify[expr_] := Module[{k, ex = expr /. {Product[f_, {j_, n_}] :> Product[Factor[f], {j,1,n}], Product[f_, {n_}] :> f^n, Product[f_, {j_, a_, b_, s_}] :> Product[Factor[f] /. j -> a + k s, {k, 0, (b-a)/s}]}, x}, ex = ex //. {Binomial[n_,k_] /; n-k =!= 0 -> rf[n-k,k]/rf[1,k] n/(n-k), Binomial[n_,k_] -> rf[n-k+1,k]/rf[1,k], Factorial[n_] -> rf[1,n], Gamma[x_] -> rf[1,x-1], Pochhammer[a_, x_] -> rf[a,x], Factorial2[a_?EvenQ n_ + b_.] /; EvenQ[b] -> 2^(a/2 n + b/2) (a/2 n + b/2)!, Factorial2[a_?EvenQ n_ + b_] /; OddQ[b] -> 2^(a/2 n + (b+1)/2) rf[1/2, a/2 n + (b+1)/2], Product[f_ g_, {j_, a_, b_}] -> Product[f, {j,a,b}] Product[g, {j,a,b}], Product[f_^k_, {j_, a_, b_}] /; FreeQ[k, j] -> Product[f, {j,a,b}]^k, Product[f_^k_, {j_, a_, b_}] /; FreeQ[f, j] :> f^GosperSum[k, {j,a,b}], Product[u_, {j_, a_, b_}] /; FreeQ[u, j] -> u^(b-a+1), Product[u_. j_ + v_., {j_, a_, b_}] /; FreeQ[u, j] && FreeQ[v, j] -> u^(b-a+1) rf[v/u + a, b-a+1], Product[(u_. j_ + v_.)^(-1), {j_, a_, b_}] /; FreeQ[u, j] && FreeQ[v, j] -> 1/(u^(b-a+1) rf[v/u + a, b-a+1])}; ex = ex /. rf[a_, b_] :> rf[Expand[a], Expand[b]]; ex = ex //. {rf[a_, b_ + n_Integer?Positive] -> rf[a, b] rf[a+b, n], rf[a_, b_ + n_Integer?Negative] -> rf[a, b] / rf[a+b+n, -n], rf[a_ + n_Integer?Positive, b_] -> rf[a, b]rf[a+b,n]/rf[a,n], rf[a_ + n_Integer?Negative, b_] -> rf[a,b]rf[a+n,-n]/rf[a+b+n,-n], rf[a_, b_. + r_Rational?(# > 1 &)] -> rf[a, b+r-Floor[r]] rf[a+b+r-Floor[r], Floor[r]], rf[a_, b_. + r_Rational?Negative] -> rf[a, b+r-Ceiling[r]+1] / rf[a+b+r, 1-Ceiling[r]], rf[a_. + r_Rational?(# > 1 &), b_] -> rf[a+r-Floor[r], b]rf[a+r-Floor[r]+b,Floor[r]]/ rf[a+r-Floor[r],Floor[r]], rf[a_. + r_Rational?Negative, b_] -> rf[a+r, 1-Ceiling[r]]rf[a+r-Ceiling[r]+1,b]/ rf[a+r+b,1-Ceiling[r]], qFac[a_, q_, b_ + n_Integer?Positive] -> qFac[a, q, b] qFac[a q^b, q, n], qFac[a_, q_, b_ + n_Integer?Negative] -> qFac[a, q, b] / qFac[a q^(b+n), q, -n]}; ex = ex /. rf[a_, b_] /; a+b =!= 0 -> Binomial[a+b,b] b! a/(a+b); ex = ex /. rf[a_, b_] -> Binomial[a+b-1,b] b!; ex = ex /. Power[a_, b_] :> Power[a, Expand[b]]; ex = ex /. (-1)^(b_+Optional[u_Integer]) Binomial[a_Integer?Negative, b_] c___ :> (-1)^u Binomial[b-a-1,-a-1] c; ex = ex /. Binomial[a_Integer?Negative, b_] :> (-1)^b Binomial[b-a-1,-a-1]; ex = ex /. Binomial[a_ + b_, b_] :> (-1)^b Binomial[-a-1,b]; ex = ex /. Binomial[Rational[m_?Positive,2], n_] :> (-1)^((m+1)/2) m!! Binomial[2n,n] / ((-1)^n 4^n Product[2n-k, {k,1,m,2}]); ex = ex /. Binomial[Rational[m_?Negative,2], n_] :> Product[2n+k, {k,1,-m-2,2}] Binomial[2n,n] / ((-1)^n 4^n (-m-2)!!); Return[PostSimplify[MyFactor[ex, {Factorial, Binomial, Factorial2, Gamma, Pochhammer, Product, qFac}]]] ]; (* FactorialSimplify[expr_] := Module[{ex = expr, x}, ex = ex //. {Binomial[n_,k_] -> rf[n-k+1,k]/rf[1,k], Factorial[n_] -> rf[1,n], Gamma[x_] -> rf[1,x-1], Pochhammer[a_, x_] -> rf[a,x], Factorial2[a_?EvenQ n_ c_.+ b_] -> (a c)^n Product[ rf[(b+k)/(a c),n], {k,2,a c,2}], Product[f_ g_, {j_, a_, b_}] -> Product[f, {j,a,b}] Product[g, {j,a,b}], Product[f_^k_, {j_, a_, b_}] /; FreeQ[k, j] -> Product[f, {j,a,b}]^k, Product[u_, {j_, a_, b_}] /; FreeQ[u, j] -> u^(b-a+1), Product[u_. j_ + v_., {j_, a_, b_}] /; FreeQ[u, j] && FreeQ[v, j] -> u^(b-a+1) rf[v/u + a, b-a+1], Product[(u_. j_ + v_.)^(-1), {j_, a_, b_}] /; FreeQ[u, j] && FreeQ[v, j] -> 1/(u^(b-a+1) rf[v/u + a, b-a+1])}; ex = ex /. rf[a_, b_] :> rf[Expand[a], Expand[b]]; ex = ex //. {rf[a_, b_ + n_Integer?Positive] -> rf[a, b] rf[a+b, n], rf[a_, b_ + n_Integer?Negative] -> rf[a, b] / rf[a+b+n, -n], rf[a_ + n_Integer?Positive, b_] -> rf[a, b]rf[a+b,n]/rf[a,n], rf[a_ + n_Integer?Negative, b_] -> rf[a,b]rf[a+n,-n]/rf[a+b+n,-n], rf[a_, b_. + r_Rational?(# > 1 &)] -> rf[a, b+r-Floor[r]] rf[a+b+r-Floor[r], Floor[r]], rf[a_, b_. + r_Rational?Negative] -> rf[a, b+r-Ceiling[r]+1] / rf[a+b+r, 1-Ceiling[r]], rf[a_. + r_Rational?(# > 1 &), b_] -> rf[a+r-Floor[r], b]rf[a+r-Floor[r]+b,Floor[r]]/ rf[a+r-Floor[r],Floor[r]], rf[a_. + r_Rational?Negative, b_] -> rf[a+r, 1-Ceiling[r]]rf[a+r-Ceiling[r]+1,b]/ rf[a+r+b,1-Ceiling[r]]}; ex = ex /. rf[a_, b_] -> Binomial[a+b-1,b] b!; ex = ex /. Power[a_, b_] :> Power[a, Expand[b]]; Return[MyFactor[ex, {Factorial, Binomial, Factorial2, Gamma, Pochhammer, Product, qFac}]] ]; *) MyFactor[ex_, heads_] := Module[{flist, vars, exp}, flist = Union[ Cases[ex, _?(MemberQ[heads, Head[#]]&), Infinity]]; vars = Table[Unique[f], {Length[flist]}]; exp = ex /. Thread[flist -> vars]; Return[Factor[exp] /. Thread[vars -> flist]] ]; Similar[r1_, r2_, n_] := Module[{r, a, b, c, a1, b1, c1, y, num, den, k}, r = Cancel[r2/r1]; {num, den} = {Numerator[r], Denominator[r]}; If[Coefficient[num, n, Exponent[num, n]] / Coefficient[den, n, Exponent[den, n]] =!= 1, Return[False]]; {a,b,c} = MixedNormalForm[num/den, n, y]; {a1,b1,c1} = MixedNormalForm[b/a, n, y]; If[Together[a1/b1] =!= 1, Return[False]]; Return[Together[c/c1]] ]; FindSimilar[diss_List, term_, n_] := Module[{r, i, out, rat1, h1, rat2, h2, val, newrat, k, mult}, out = Scan[( {{rat1, h1}, i} = #; {rat2, h2} = term; r = Similar[rat1, rat2, n]; If[r =!= False, k = 1; Off[Power::infy, Infinity::indet]; val = h2/(h1 r) /. n -> k; While[!FreeQ[val, Indeterminate] || !FreeQ[val, DirectedInfinity], k++; val = h2/(h1 r) /. n -> k]; On[Power::infy, Infinity::indet]; mult = Together[1 + FS[val] r]; If[mult === 0, Return[Delete[diss, i]]]; newrat = Factor[rat1 (mult /. n->n+1) / mult]; Return[ReplacePart[diss, {{newrat, h1 mult}, i}, i]] ] )&, diss]; If[out === Null, Return[Append[diss, {term, Length[diss] + 1}]], Return[Thread[{First /@ out, Range[Length[out]]}]]] ]; FactorialSimplify[expr_, n_] /; RationalQ[expr, n] := Together[expr]; FactorialSimplify[expr_, n_] := Module[{den, ex, rat, ratios, diss}, den = Denominator[expr]; If[den =!= 1, FS[Numerator[expr], n] / FS[den, n]]; ex = Together /@ MakeList[expr, Plus]; rat = Select[ex, RationalQ[#, n]&]; ex = Complement[ex, rat]; rat = Together[Plus @@ rat]; ratios = GetRatio[#, n]& /@ ex; ratios = Factor /@ ratios; diss = Fold[FindSimilar[#1, #2, n]&, {}, Thread[{ratios, ex}]]; diss = First /@ diss; Return[rat + PostSimplify[Plus @@ (FormHg[#, n]& /@ diss)]]]; MaxRoot[p_, x_] := (* max nonnegative integer root, or -1 *) Module[{all}, all = x /. Select[Solve[p == 0, x], FreeQ[#, Roots]&]; If[all === x, all = {}]; all = Select[all, IntegerQ[#] && # >= 0&]; If[all === {}, all = {-1}]; Return[Max[all]] ]; Transform[(n_ + a_Integer?Negative)^k_., n_] /; FreeQ[k, n] := (n+a-1)!^k; Transform[n_^k_., n_] /; FreeQ[k,n] := (n-1)!^k; Transform[(n_ + a_)^k_., n_] /; FreeQ[{a,k}, n] := Binomial[n+a-1,n]^k n!^k; (* (-1)^n Binomial[-a,n]^k n!^k; *) Transform[a_, n_] := Product[a /. n->j, {j, n-1}]; ToHg[{r_, h_}, n_] := (* constructs a canonical hg. term for h from its c.t.r. r *) Module[{rr, lc, k, val}, rr = Factor[r] /. {a_?(FreeQ[#,n]&) n + b_?(FreeQ[#,n]&) :> (n + b/a) a}; rr = MakeList[rr, Times]; lc = Select[rr, FreeQ[#, n]&]; rr = Complement[rr, lc]; lc = Times @@ lc; Scan[If[MatchQ[#, a_?(FreeQ[#,n]&) - n], lc = -lc]&, rr]; rr = rr /. {a_?(FreeQ[#,n]&) - n -> n - a}; rr = Transform[#, n]& /@ rr; rr = rr /. Binomial[m_, k_] -> Binomial[k-m-1,k]/(-1)^k; rr = lc^n (Times @@ rr); k = 1; Off[Power::infy, Infinity::indet]; val = (h/rr) /. n -> k; While[!FreeQ[val, Indeterminate] || !FreeQ[val, DirectedInfinity], k++; val = (h/rr) /. n -> k]; On[Power::infy, Infinity::indet]; Return[rr val] ]; FormHg[{rat_, h_}, n_] := Module[{a, b, c, a1, b1, c1, y}, {a,b,c} = MixedNormalForm[rat, n, y]; {a1,b1,c1} = MixedNormalForm[b/a, n, y]; Return[Together[c/c1] ToHg[{b1/a1, h c1/c}, n]] ]; PostSimplify[expr_] := Module[{ex = expr /. n_! :> (Expand[n])!, i}, ex //. {n_!^k_. m_!^l_. :> Product[m + i, {i, n - m}]^k /; IntegerQ[Expand[n - m]] && Expand[n - m] > 0 && Expand[k + l] == 0, n_!^k_. m_^l_. :> Expand[n + 1]!^k Expand[n + 1]^(l - k) /; Expand[m - n] == 1 && Sign[k] == Sign[l], n_!^k_. m_^l_. :> (-1)^k Expand[n + 1]!^k Expand[m]^(l - k) /; Expand[m + n] == -1 && Sign[k] == Sign[l], n_!^k_. n_^l_. :> Expand[n - 1]!^k /; Expand[k + l] == 0, n_!^k_. m_^l_. :> (-1)^k Expand[n - 1]!^k /; Expand[k + l] == 0 == Expand[n + m] } ]; (* MIXED NORMAL FORM *) FindShift[a_, b_, x_] := Module[{d = Exponent[a, x], a0, a1, b0, b1, h}, If[Exponent[b, x] =!= d, Return[{}]]; If[d === 0, Return[{0, {a, b}, a/b}]]; a0 = Coefficient[a, x, d]; a1 = Coefficient[a, x, d-1]; b0 = Coefficient[b, x, d]; b1 = Coefficient[b, x, d-1]; h = Together[Det[{{a1,b1},{a0,b0}}]/(d a0 b0)]; If[!IntegerQ[h] || h < 0 || Expand[b0 a - a0 (b /. x -> x+h)] =!= 0, {}, {h, {a, b}, a0/b0}] ]; FindExponent[flist_, af_] := Scan[If[#[[1]] === af, Return[#[[2]]]]&, flist] /. Null -> 0; FactorOut[{alist_, blist_, a_, b_, c_}, {h_, {af_, bf_}, const_}, x_] := Module[{d = Min[FindExponent[alist, af], FindExponent[blist, bf]]}, {alist/.{af,n_}->{af,n-d}, blist/.{bf,n_}->{bf,n-d}, a const^d / af^d, b / bf^d, c Product[bf /. x -> x+s, {s,0,h-1}]^d} ]; NormalForm[r_, x_] := Module[{fr = Factor[r], a, b, alist, blist, shifts}, a = Numerator[fr]; b = Denominator[fr]; alist = Select[FactorList[a], !FreeQ[#, x]&]; blist = Select[FactorList[b], !FreeQ[#, x]&]; shifts = Select[Union @@ Outer[FindShift[#1, #2, x]&, First /@ alist, First /@ blist], # =!= {}&]; Drop[Fold[FactorOut[#1,#2,x]&, {alist, blist, a, b, 1}, shifts], 2] ]; ToY[expr_, n_, qList_, y_] := Together[PowerExpand[expr] /. q_?(MemberQ[qList,#]&)^(a_. n + b_.) :> q^b y[q]^a]; GetMonomial[term_, y_] := Which[term === 0, 0, FreeQ[term, y], 1, Head[term] =!= Times, term, True, Select[term, Not[FreeQ[#, y]]&] ]; GetCoefficient[expr_, mon_, y_] := If[mon === 1, Plus @@ Select[MakeList[expr, Plus], FreeQ[#, y]&], GetCoefficient[Coefficient[expr, mon], 1, y] ]; MixedFindShift[a_, b_, n_, y_] := (* returns {h, {a, b}, const} where h >= 0 and a(n) = const b(n+h), or {} if no such h and const exist; it is assumed that gcd(a,b) = 1 *) Module[{aa = MakeList[a, Plus], bb = MakeList[b, Plus], e, max, min, amax, amin, bmax, bmin, q, qn, qp, h, exp, quot}, (* compare lists of y-monomials *) aa = Union[GetMonomial[#, y]& /@ aa]; bb = Union[GetMonomial[#, y]& /@ bb]; If[aa =!= bb, Return[{}]]; (* if there are no y's use FindShift *) If[aa === {1}, Return[FindShift[a, b, n]]]; (* use max and min monomials to determine h and const (since gcd(a,b) = 1, there are at least two distinct monomials in a and b) *) max = Last[aa]; min = First[aa]; (* get the exponent of the first base appearing in max/min *) qn = First[MakeList[max/min, Times]]; {q, exp} = Replace[qn, y[q_]^exp_. -> {q, exp}]; (* compare leading coeffs w.r.t. n of the coeffs of max and min *) amax = GetCoefficient[a, max, y]; bmax = GetCoefficient[b, max, y]; amax = Coefficient[amax, n, Exponent[amax, n]]; bmax = Coefficient[bmax, n, Exponent[bmax, n]]; amin = GetCoefficient[a, min, y]; bmin = GetCoefficient[b, min, y]; amin = Coefficient[amin, n, Exponent[amin, n]]; bmin = Coefficient[bmin, n, Exponent[bmin, n]]; quot = MakeList[Cancel[(amax/bmax)/(amin/bmin)], Times]; qp = Select[quot, !FreeQ[#, q]&]; (* there is no shift unless qp = {q^e} and h = e/exp is a nonnegative integer *) h = Together[Replace[qp, {q^e_.} -> e] / exp]; If[!IntegerQ[h] || h < 0, Return[{}]]; (* compute overall constant factor *) const = (amin/bmin) / (min /. y[q_] -> q)^h; (* check shift and const *) If[Together[a - (const b /. {n -> n+h, y[q_] -> q^h y[q]})] =!= 0, {}, {h, {a, b}, const}] ]; MixedFactorOut[{alist_, blist_, a_, b_, c_}, {h_, {af_, bf_}, const_}, n_, y_] := Module[{d = Min[FindExponent[alist, af], FindExponent[blist, bf]]}, {alist /.{af, k_} -> {af, k - d}, blist /.{bf, k_} -> {bf, k - d}, Cancel[a const^d / af^d], Cancel[b / bf^d], c Product[bf /. {n -> n+s, y[q_] -> q^s y[q]}, {s,0,h-1}]^d } ]; MixedNormalForm[r_, n_, y_] := Module[{fr = Cancel[r], a, b, alist, blist, shifts}, a = Numerator[fr]; b = Denominator[fr]; alist = Select[MyFactorList[a], !FreeL[#, {n, y}]&]; blist = Select[MyFactorList[b], !FreeL[#, {n, y}]&]; shifts = Select[Union @@ Outer[MixedFindShift[#1, #2, n, y]&, First /@ alist, First /@ blist], # =!= {}&]; Drop[Fold[MixedFactorOut[#1,#2,n,y]&, {alist, blist, a, b, 1}, shifts], 2] ]; MyFactorList[x__Times] := Join @@ (MyFactorList /@ (List @@ x)); MyFactorList[x_^n_] := {{x, n}}; MyFactorList[x_] := {{x, 1}}; (* MIXED POLY *) MixedPoly[eqn_, f_[n_], qList_, cList_:{}] := Module[{y, eq, lhs, rhs, args, min, max, ord, terms, mons, maxd, ltcoeffs, chp, nmax, roots, coeffs, qmon, nt, sol, sub, newc, genh, part, cc, x, constants}, (* put everything to the left *) If[Head[eqn] === Equal, eq = eqn[[1]] - eqn[[2]], eq = eqn]; (* shift if necessary *) args = First /@ Cases[eq, _f, Infinity] /. n -> 0; {min, max} = {Min[args], Max[args]}; eq = eq /. n -> n - min; ord = max - min; (* separate left-hand side from right-hand side *) rhs = - Select[eq, FreeQ[#, f]&]; lhs = eq + rhs; (* get polynomial coefficients of L *) coeffs = Coefficient[lhs, f[n+#]]& /@ Range[0, ord]; (* replace q^n by y[q] *) {coeffs, rhs} = Expand[{coeffs, rhs}] /. q_?(MemberQ[qList,#]&)^(a_. n + b_.) :> q^b y[q]^a; (* get q-monomials in L *) terms = MakeList[#, Plus]& /@ coeffs; terms = Join @@ terms; mons = Function[x, If[FreeQ[x, y], 1, If[Head[x] =!= Times, x, Select[x, Not[FreeQ[#, y]]&] ]]] /@ terms; mons = Union[mons]; maxd = Last[mons]; (* characteristic polynomial and its roots *) ltcoeffs = If[maxd === 1, coeffs, Expand[Coefficient[#, maxd]& /@ coeffs]]; nmax = Max[Exponent[#, n]& /@ ltcoeffs]; chp = (Coefficient[#, n, nmax]& /@ ltcoeffs).(x^Range[0, ord]); roots = x /. Select[Solve[chp == 0, x], FreeQ[#, Roots]&]; If[roots === x, roots = {}]; qmon[a_ b_] := qmon[a] && qmon[b]; qmon[a_^b_] := MemberQ[qList, a] && IntegerQ[b] && b >= 0; qmon[a_] := MemberQ[Append[qList, 1], a]; roots = Select[Union[roots], qmon] /. q_?(MemberQ[qList,#]&)^a_. :> y[q]^a; roots = Prepend[roots, 0]; (* call the main routine *) nt = NextTerm[coeffs, ltcoeffs, maxd, rhs, n, roots, qList, cList, y]; If[nt === $Failed, Return[{}]]; {sol, sub, newc} = nt; (* spruce up the solution *) sol = MakeList[sol /. sub, Plus]; sol = sol /. y[q_] -> q^n; constants = Join[cList, newc]; part = Select[sol, FreeL[#, constants]&]; genh = Complement[sol, part]; genh = Collect[Plus @@ genh, constants]; genh = MakeList[genh, Plus]; genh = Plus @@ (Together /@ genh); (* absorb constant factors into arbitrary constants *) genh = genh /. {a_?(FreeQ[#, n]&) b_?(MemberQ[newc,#]&) c_. -> b c}; sol = Factor[Plus @@ part] + genh; (* rename and renumber arbitrary constants *) cc = Cases[{sol}, _?(MemberQ[newc,#]&), Infinity]; Return[{sol, Select[Thread[cList -> (cList //. sub)], #[[1]] =!= #[[2]]&]} /. Thread[cc->Array[C, Length[cc]]] ] ]; NextTerm[coeffs_, ltcoeffs_, maxd_, rhs_, n_, roots_, qList_, cList_, y_] := (* finds polynomial solutions of Lf == rhs *) Module[{maxroot, terms, mons, maxr, bound, lt, a, q, pl, sol, x, range = Range[0, Length[coeffs]-1], newroots, newc, nt}, maxroot = Last[roots]; If[maxroot === 0 && rhs === 0, Return[{0,{},{}}]]; (* rhs degree *) terms = MakeList[Expand[rhs], Plus]; mons = Function[x, Which[x === 0, 0, FreeQ[x, y], 1, Head[x] =!= Times, x, True, Select[x, Not[FreeQ[#, y]]&] ]] /@ terms; mons = Union[mons]; maxr = Last[mons]; (* main switch *) If[OrderedQ[{maxr, maxd maxroot}], bound = maxroot; newroots = Drop[roots, -1], If[Denominator[maxr/maxd] === 1, bound = maxr/maxd; newroots = roots, bound = None; If[cList === {}, sol = {}, Off[Solve::svars]; sol = Solve[CoefficientList[Expand[ If[maxr === 1, rhs, Coefficient[rhs, maxr]]], n] == 0, cList] ]; If[sol === {} || sol === {{}}, Return[$Failed], sol = sol[[1]]; nt = NextTerm[coeffs, ltcoeffs, maxd, Expand[Together[rhs /. sol]], n, roots, qList, cList, y]; If[nt === $Failed, Return[$Failed]]; Return[{nt[[1]], Join[sol /. nt[[2]], nt[[2]]], nt[[3]]}] ]]]; (* leading term: call Poly to get the (q^n)-free part *) If[bound =!= None, pl = Poly[ltcoeffs . (a[n+#](bound /. y[q_] -> q)^# & /@ range) == If[maxd bound === 1, rhs, Coefficient[rhs, maxd bound]], a[n], cList]; If[pl === {}, Return[$Failed], {lt, sol, newc} = pl]; lt = lt bound; (* recursive call *) nt = NextTerm[coeffs, ltcoeffs, maxd, Expand[Together[(rhs - coeffs . ((lt /. {n -> n+#, y[q_] -> q^# y[q]})& /@ range)) /. sol]], n, newroots, qList, Join[cList, newc], y]; If[nt === $Failed, Return[$Failed]]; Return[{lt + nt[[1]], Join[sol /. nt[[2]], nt[[2]]], Join[newc, nt[[3]]]}] ]]; MakeList[expr_, head_: List] := If[Head[expr] === head, List @@ expr, List @ expr]; FreeL[eq_, un_] := FreeQ[eq, _?(MemberQ[un, #]&)]; Coef[_,_,_,_?Negative] = 0; Coef[p_,i_,n_,j_] := Coefficient[p[[i]], n, j]; TriSolve[eqns_, unknowns_] := Module[{sub = {}, sub1, x, eq, eqn, term, un}, sub1 = Scan[(eq = Collect[Numerator[Together[# /. sub]], unknowns]; If[eq == 0, Null, (* this equation is satisfied *) Return[Fail], (* lhs != 0 over Q *) If[FreeL[eq, unknowns], Return[Fail], (* lhs != 0 over Q(params) *) eqn = Reverse[MakeList[eq, Plus]]; (* find a term containing an unknown *) term = Scan[If[Or @@ (Function[x, !FreeQ[#, x]] /@ unknowns), Return[#]]&, eqn]; (* extract the unknown *) un = First[Select[MakeList[term, Times], MemberQ[unknowns, #]&]]; sub1 = {un -> (term - eq)/(term/un)}; sub = Join[sub /. sub1, sub1] ] ]) & , eqns]; Return[If[sub1 === Fail, Fail, sub]]]; Poly[eqn_, a_[n_], cList_, verbose_: False] := Module[{min, max, dd, polyCoeffs, i, eq, lhs, rhs}, eq = Numerator[Together[ If[Head[eqn] === Equal, eqn[[1]]-eqn[[2]], eqn] ]]; {min, max} = {Min[#], Max[#]}& @ (First /@ Cases[{eq}, a[_], Infinity]/.n->0); dd = max - min; eq = Collect[Expand[eq /. n -> n - min], Table[a[n+i], {i,0,dd}]]; rhs = - Select[eq, FreeQ[#, a]&]; lhs = eq + rhs; polyCoeffs = Module[{i}, Coefficient[lhs, #]& /@ Table[a[n+i], {i,0,dd}]]; Return[PolyK[polyCoeffs, rhs, a[n], cList, verbose]] ]; PolyK[polyCoeffs_, rhs_, a_[n_], cList_, verbose_: False] := Module[{dd = Length[polyCoeffs] - 1, m, degreePoly, s = 0, d, degrees, deg, unknowns, lhs, i, sol, poly, sub, k, cc, genh, part, verb = MemberQ[{True, Automatic}, verbose], C}, m = Max[Exponent[#, n]& /@ polyCoeffs]; For[degreePoly = 0, degreePoly === 0, s++, degreePoly = Module[{i,j}, Factor[Sum[Binomial[d,j] Sum[i^j Coef[polyCoeffs, i+1, n, m-s+j], {i,0,dd}], {j,0,s}]]]]; s--; If[verb, Print[""]; Print["Degree polynomial: ", degreePoly/.d->n]]; degrees = Flatten[Union[Factor /@ (d/.Solve[degreePoly==0, d]) /. d -> {}]]; If[verb, Print["Roots found: ", degrees]]; degrees = Append[degrees, Exponent[rhs, n] - m + s]; degrees = Select[degrees, (IntegerQ[#] && NonNegative[#])&]; If[verb, Print["Possible degrees: ", degrees]]; deg = Max[degrees]; If[deg === -Infinity, If[rhs === 0, Return[{0,{},{}}], Return[{}] ]]; unknowns = Array[C, deg+1, 0]; poly = unknowns.Table[n^k, {k,0,deg}]; lhs = polyCoeffs . Table[a[n+i], {i,0,dd}]; sub = lhs /. a[x_] :> (poly /. n -> x); sub = Reverse[CoefficientList[Expand[sub - rhs], n]]; If[verbose, Print["Triangular system: "]; Print[""]; Print[MatrixForm[Thread[sub == 0]]]]; sub = TriSolve[sub, Join[unknowns, cList]]; If[sub === Fail, Return[{}], sol = Collect[poly /. sub, unknowns]]; If[sol === 0, Return[{0,{},{}}]]; sol = MakeList[sol, Plus]; part = Select[sol, FreeQ[#, C]&]; genh = Complement[sol, part]; sol = Factor[Plus @@ part] + Plus @@ (Factor[Numerator[Together[#]]]& /@ genh); cc = Cases[{sol}, C[_], Infinity]; sub = Simplify[sub]; Return[{sol, Select[sub, !FreeL[First[#], cList]&], cc} /. Thread[cc->Array[C, Length[cc]]] ]]; (* GET RATIO *) (* q-factorials *) qFac[a_, q_, n_Integer] := Which[ n > 0, Product[1 - a q^i, {i,0,n-1}], n === 0, 1, n < 0, 1 / Product[1 - a/q^i, {i,-n}] ]; (** GetRatio **) (* distribute over factors *) GetRatio[a_Times, n_] := GetRatio[#, n]& /@ a; GetRatio[a_^b_, n_] := GetRatio[a, n]^b /; FreeQ[b, n]; (* built-in hypergeometrics *) GetRatio[Gamma[x_], n_] := GetRatio[Pochhammer[1, x-1], n]; GetRatio[a_!, n_] := GetRatio[Pochhammer[1, a], n]; GetRatio[(a_. n_ + b_.)!!, n_] := Product[a n + b + 2 i, {i, a/2}] /; IntegerQ[a/2] && a > 0; GetRatio[(a_. n_ + b_.)!!, n_] := 1 / Product[a n + b + 2 i, {i, 0, a/2 + 1, -1}] /; IntegerQ[a/2] && a < 0; GetRatio[Binomial[a_,b_], n_] := GetRatio[Pochhammer[a-b+1,b]/b!, n]; GetRatio[Pochhammer[a_. n_ + b_., c_. n_ + d_.], n_] := (If[a < 0, Product[a n + b + i, {i, a, -1}], 1/Product[a n + b + i, {i, 0, a-1}]] * If[a+c > 0, Product[(a+c) n + b + d + i, {i, 0, a+c-1}], 1/Product[(a+c) n + b + d + i, {i, a+c, -1}]]) /; IntegerQ[a] && IntegerQ[c] && FreeQ[b, n] && FreeQ[d, n]; GetRatio[Pochhammer[a_. n_ + b_., d_], n_] := (If[a < 0, Product[a n + b + i, {i, a, -1}], 1/Product[a n + b + i, {i, 0, a-1}]] * If[a > 0, Product[a n + b + d + i, {i, 0, a-1}], 1/Product[a n + b + d + i, {i, a, -1}]]) /; IntegerQ[a] && FreeQ[d, n] && FreeQ[b, n]; GetRatio[Pochhammer[b_, c_. n_ + d_.], n_] := (If[c > 0, Product[c n + b + d + i, {i, 0, c-1}], 1/Product[c n + b + d + i, {i, c, -1}]]) /; IntegerQ[c] && FreeQ[b, n] && FreeQ[d, n]; (* products *) GetRatio[Product[a_, {i_, b_, c_. n_ + d_.}], n_] := (Product[a, {i, c n + d + 1, c n + d + c}]) /; IntegerQ[c] && c >= 0 && FreeQ[b, n] && RationalQ[a, i]; GetRatio[Product[a_, {i_, b_}], n_] := GetRatio[Product[a, {i,1,b}], n]; RationalQ[p_, x_] := PolynomialQ[Numerator[p], x] && PolynomialQ[Denominator[p], x]; (* q-hypergeometrics *) GetRatio[qFac[a_, q_, b_. n_ + c_.], n_] := qFac[a q^(b n + c), q, b] /; IntegerQ[b] && FreeQ[a, n] && FreeQ[q, n] && FreeQ[c, n]; GetRatio[qFac[A_. q_^(a_. n_ + b_.), q_^i_., c_. n_ + d_.], n_] := qFac[A q^((a + i c) n + b + i d), q^i, a/i + c] / qFac[A q^(a n + b), q^i, a/i] /; IntegerQ[a/i] && FreeQ[b, n] && FreeQ[q, n] && FreeQ[c, n] && FreeQ[d, n]; GetRatio[qFac[A_. q_^(a_. n_ + b_.), q_^i_., d_], n_] := qFac[A q^(a n + b + i d), q^i, a/i] / qFac[A q^(a n + b), q^i, a/i] /; IntegerQ[a/i] && FreeQ[b, n] && FreeQ[q, n] && FreeQ[d, n]; (* the rest *) GetRatio[a_, n_] := FS[(a /. n->n+1) / a]; (* MIXED GOSPER *) MyCancel[f_Times] := Module[{f1, f2}, f1 = Select[f, FreeL[#, {Binomial, Factorial, Factorial2, Pochhammer, Gamma, qFac}]&]; f2 = f/f1; Return[Cancel[f1] f2] ]; MyCancel[f_] := f; GosperSum[t_, k_Symbol, qList_:{}] := (* indefinite sum *) Module[{r, f}, r = Together /@ GetRatio[t, k]; f = GosperFunction[r, k, qList]; If[f === $Failed, Return[$Failed], Return[MyCancel[f t] // PostSimplify]]]; GosperSum[t_, {k_, a_, b_}, qList_:{}] := (* definite sum *) Module[{r, f}, r = Together /@ GetRatio[t, k]; f = GosperFunction[r, k, qList]; If[f === $Failed, Return[Sum[t, {k, a, b}]], Return[(MyCancel[(t /. k -> b) Factor[1 + (f /. k -> b)]] - ((f t) /. k -> a)) // PostSimplify] ]]; GosperFunction[r_, k_Symbol, qList_:{}] := (* Gosper's algoritem *) Module[{y, num, den, a, b, c, b1, z, x, T, qmon, f, w, out, a0, b0, g, aa, bb, cc, wList, rr, vars, a2, b2, c2}, rr = ToY[r, k, qList, y]; vars = Prepend[y /@ qList, k]; num = Numerator[rr]; den = Denominator[rr]; If[!PolynomialQ[num, vars] || !PolynomialQ[den, vars], Return[$Failed]]; {a,b,c} = MixedNormalForm[num/den, k, y]; b1 = b /. {k -> k-1, y[q_] -> y[q]/q}; (* try with g = 1 first *) {a2,b2,c2} = {a,b1,c} /. {y[q_]^m_. -> q^(m k)}; z = MixedPoly[a2 f[k+1] - b2 f[k] == c2, f[k], qList] /. _C -> 0; If[z =!= {}, Return[Together[(b1 z[[1]] / c) /. y[q_]^n_. -> q^(n k)]] ]; (* find g != 1 *) T = Select[y /@ qList, FreeQ[{a, b} /. # -> 0, 0]&]; qmon[a_ b_] := qmon[a] && qmon[b]; qmon[a_^b_] := MemberQ[qList, a] && IntegerQ[b]; qmon[a_] := MemberQ[Append[qList, 1], a]; exponents = Map[Function[yi, a0 = a /. yi -> 0; b0 = b1 /. yi -> 0; {aa,bb,cc} = MixedNormalForm[b0/a0, k, y]; If[!qmon[aa/bb], 0, Exponent[Denominator[aa/bb], yi[[1]] ] ]], T]; w = Times @@ ((First /@ T) ^ exponents); g = Times @@ (T ^ exponents); {w2, g2} = {w, g} /. {y[q_]^m_. -> q^(m k)}; If[g === 1, Return[$Failed]]; z = MixedPoly[a2 f[k+1] - b2 w2 f[k] == c2 w2 g2, f[k], qList] /. _C -> 0; If[z === {}, Return[$Failed], Return[Together[(b1 z[[1]] / (g c)) /. y[q_]^n_. -> q^(n k)]] ] ]; Options[Zeil] = {Bases -> {}, Rhs -> True, Certificate -> False, Start -> 0}; Zeil[F_, s_[n_], k_Symbol, opts___] := Zeil[F, s[n], {k, -Infinity, Infinity}, opts]; Zeil[F_, s_[n_], {k_, lo_, hi_}, opts___] := Module[{d, f, g, u, uu, v, vv, vp, out, alpha, pars, p, x, a, b, c, b1, y, z, vars, qList, rhsQ, certQ, rec, cert, dens, den, lclo, lchi, lolim, hilim, rhs, i, FF}, d = Start /. {opts} /. Options[Zeil]; qList = Bases /. {opts} /. Options[Zeil]; rhsQ = Rhs /. {opts} /. Options[Zeil]; certQ = Certificate /. {opts} /. Options[Zeil]; If[!PolynomialQ[lo, n] || !PolynomialQ[hi, n], Print[""]; Print["Non-polynomial limits."]; Return[]]; lclo = Coefficient[lo, n, Exponent[lo, n]]; lchi = Coefficient[hi, n, Exponent[hi, n]]; lolim = Expand[If[lclo > 0, lo /. n->n+d, lo]]; hilim = Expand[If[lchi < 0, hi /. n->n+d, hi] - 1]; alpha[0] = 1; f = Together /@ GetRatio[F, n]; g = Together /@ GetRatio[F, k]; f = ToY[ToY[f, k, qList, z], n, qList, y]; g = ToY[ToY[g, k, qList, z], n, qList, y]; u = Numerator[f]; v = Denominator[f]; vars = Join[{k, n}, y /@ qList, z /@ qList]; If[!PolynomialQ[u, vars] || !PolynomialQ[v, vars] || !PolynomialQ[Numerator[g], vars] || !PolynomialQ[Denominator[g], vars], Print[""]; Print["Non-hypergeometric input."]; Return[]]; out = While[True, Print[""]; Print["trying d = ", d, " ..."]; uu = FoldList[Times, 1, Table[u /. {n -> n+i, y[q_] -> q^i y[q]}, {i,0,d-1}]]; vv = FoldList[Times, 1, Table[v /. {n -> n+i, y[q_] -> q^i y[q]}, {i,d-1,0,-1}]]; vp = Last[vv]; {a,b,c} = MixedNormalForm[g vp / (vp /. {k -> k+1, z[q_] -> q z[q]}), k, z]; b1 = b /. {k -> k-1, z[q_] -> z[q]/q}; pars = alpha /@ Range[0, d]; p = Plus @@ (pars uu Reverse[vv]); {a,b1,c,p} = {a,b1,c,p} /. {y[q_] -> q^n, z[q_] -> q^k}; sol = MixedPoly[a x[k+1] - b1 x[k] == c p, x[k], qList, Rest[pars]]; If[sol === {}, d++, Print["Solution found."]; rec = Together[Reverse[Expand[(pars /. sol[[2]]) /. _alpha -> 0]]]; dens = Denominator /@ rec; den = PolynomialLCM @@ dens; rec = Together[den rec].(s /@ Range[n+d,n,-1]); If[Not[certQ || rhsQ], Return[rec]]; cert = (den sol[[1]]b1/c/vp) /. {y[q_] -> q^n, z[q_] -> q^k, _alpha -> 0}; If[Not[rhsQ], Return[{cert, rec}]]; FF = FS[cert F, k]; rhs = (rec /. s[n+i_.] :> Sum[Evaluate[F /. n->n+i], {k, Expand[lo /. n->n+i], lolim - 1}] + Sum[Evaluate[F /. n->n+i], {k, hilim + 1, Expand[hi /. n->n+i]}]) + FS[FF /. k -> hilim + 1, n] - FS[FF /. k -> lolim, n]; If[FreeQ[rhs, Sum], rhs = FS[rhs, n]]; If[Not[certQ], Return[rec == rhs], Return[{cert, rec == rhs}]]; (* Return[{sol[[1]]b1/c/vp /. {y[q_] -> q^n, z[q_] -> q^k}, Reverse[Expand[(pars /. sol[[2]]) /. _alpha -> 0]]}] *) ] ] ]; Print["N.B.: Besides GosperSum and GosperFunction, this package also contains FactorialSimplify (alias FS), and WZ."]; End[]; EndPackage[]; On[General::spell1];