(* 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 = rr /. Binomial[Rational[m_?Negative,2], u_] :>
         Product[2u+k, {k,1,-m-2,2}] Binomial[2u,u](-1)^u /(4^u (-m-2)!!);
    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},
        Cancel[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[f],
          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];
  