(**************************************************************************)
(*                                                                        *)
(*                   Algorithms Poly, RatSol and Hyper                    *)
(*                                                                        *)
(*                    Version of January 6, 1998                          *)
(*                 programmed by Marko Petkov"sek                         *)
(*                   IMFM, Ljubljana, Slovenia                            *)
(*                   marko.petkovsek@fmf.uni-lj.si                        *)
(*                                                                        *)
(* This file contains Mathematica implementations of algorithms Poly,     *)
(* RatSol and Hyper which find polynomial resp. rational resp.            *)
(* hypergeometric solutions of linear recurrences with rational           *)
(* coefficients. For details of Poly and Hyper, see the paper             *)
(*                                                                        *)
(*    M. Petkov"sek: Hypergeometric solutions of linear recurrences       *)
(*    with polynomial coefficients, J. Symb. Comp. 14 (1992) 243 - 264.   *)
(**************************************************************************)

BeginPackage["Hyper`"]

Verbose::usage = "Verbose is an option of Poly, RatSol and Hyper specifying 
the printout level. Default value: False. Other values: Automatic, True."

Solutions::usage = "Solutions is an option of Hyper specifying the number of
solutions returned. Default value: 1. Other values: All."

Quadratics::usage = "Quadratics is an option of Hyper specifying whether 
quadratic factors of the leading and trailing coeficients of the recurrence
irreducible over Q should be completely factored. Default value: False. Other 
values: True."

Poly::usage = "Poly[eqn, y[n]] finds a basis for the (linear or affine) 
space of polynomial solutions of eqn which can be nonhomogeneous."

RatSol::usage = "RatSol[eqn, y[n]] finds a basis for the (linear or affine) 
space of rational solutions of eqn which can be nonhomogeneous."

Hyper::usage = "Hyper[eqn, y[n]] finds at least one hypergeometric 
solution of the homogeneous equation eqn over the field of rational numbers Q 
(provided any such solution exists). Hyper[eqn, y[n], Solutions -> All] finds 
a generating set (not necessarily linearly independent) for the space of 
solutions generated by hypergeometric terms over Q. Hyper[eqn, y[n], 
Quadratics -> True] finds solutions over quadratic extensions of Q. Solutions 
y[n] are given by their rational representations y[n+1]/y[n]. 
Warning: The worst-case time complexity of Hyper is exponential 
in the degrees of the leading and trailing coefficients of eqn."

HyperTerm::usage = "HyperTerm[eqn, y[n], y[n0] == y0] works like Hyper except that instead
of the rational representations y[n+1]/y[n] it returns the solutions y[n]
themselves, normalized so that y[n0] = y0. Here y and n must be symbols." 



Begin["`Private`"];

Unprotect[Power];
  0^0 = 1;
Protect[Power];

Options[Hyper] = {Verbose -> False, Solutions -> 1, Quadratics -> False};

Options[Poly] = {Verbose -> False};

Options[RatSol] = {Verbose -> False};

HyperTerm[eqn_, a_[n_], a_[n0_]==a0__, opts___] :=
   Module[{sol},
      sol = Hyper[eqn, a[n], opts];
      sol = ToHg[#, n, n0, a0]& /@ sol;
      Return[sol]
   ];

ToHg[f_, n_, n0_, a0_] :=
  Module[{ff = Factor[f]},
    ff = ff /. (a_. n + b_.)^k_. -> a^k (n + b/a)^k;
    ff = If[Head[ff] === Times, List @@ ff, {ff}];
    ff = Replace[#, (n + b_.)^k_. :> (Factorial[b+n-1])^k]& /@ ff;
    ff = Replace[#, c_ /; FreeQ[c, n] -> c^n]& /@ ff;
	ff = Times @@ ff;
	ff = a0 ff / (ff /. n -> n0);
    Return[GosperSum`FS[ff, n]] ];

Poly[eqn_, a_[n_], opts___] :=
   Module[{min, max, dd, polyCoeffs, i, eq, lhs, rhs, verbose},
   verbose = Verbose /. {opts} /. Options[Poly];
   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], verbose]] ];         


Coeff[_, _, _?Negative] = 0;
 
Coeff[p_, n_, i_] := Coefficient[p, n, i];
 
Bags[{}, {}] = Bags[{0}, {0}];
 
Bags[exps_, degs_] := 
  Module[{tuples, values, i}, tuples = 
     Flatten[Apply[Outer, Join[{List}, (Range[0, #1] & ) /@ exps]], 
      Length[exps] - 1]; values = (Apply[Plus, #1 . degs] & ) /@ tuples; 
    Return[Table[tuples[[Apply[Join, Position[values, i]]]], 
      {i, 0, exps . degs}]]];
 
PolyK[polyCoeffs_, rhs_, (a_)[n_], 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]}, 
   m = Max[(Exponent[#1, 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[#1] && NonNegative[#1] & ]; 
    If[verb, Print["Possible degrees: ", degrees]]; deg = Max[degrees]; 
    If[deg === -Infinity, 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, unknowns]; 
    If[sub === Fail, Return[{}], sol = Collect[poly /. sub, unknowns]]; 
    If[sol === 0, Return[{}]]; sol = MakeList[sol, Plus]; 
    part = Select[sol, FreeQ[#1, C] & ]; genh = Complement[sol, part]; 
    sol = Factor[Apply[Plus, part]] + 
      Apply[Plus, (Factor[Numerator[Together[#1]]] & ) /@ genh]; 
    cc = Cases[{sol}, C[_], Infinity]; 
    Return[sol /. Thread[cc -> Array[C, Length[cc]]]]];
 
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[#1 /. sub]], unknowns]; 
         If[eq == 0, Null, Return[Fail], 
          If[FreeL[eq, unknowns], Return[Fail], 
           eqn = MakeList[eq, Plus]; 
            term = 
             Scan[If[Apply[Or, Function[x, !FreeQ[#1, x]] /@ unknowns], 
                Return[#1]] & , eqn]; 
            un = 
             First[Select[MakeList[term, Times], 
               MemberQ[unknowns, #1] & ]]; 
            sub1 = {un -> (term - eq)/(term/un)}; 
            sub = Join[sub /. sub1, sub1]]]) & , eqns]; 
    Return[If[sub1 === Fail, Fail, sub]]];
 
FreeL[expr_, l_] := 
  Apply[And, (FreeQ[expr, #] & ) /@ l];
 
MakeList[expr_, head_:List] := 
  If[Head[expr] === head, 
   Apply[List, expr], {expr}];
 
Hyper[eqn_, a_[n_], opts___] := 
  Module[{min, max, dd, polyCoeffs, m, degreePoly, s = 0, d, degrees, deg, 
    sol, poly, sub, k, eq, degseq, lc, tc, ld, td, Alist, Blist, Zlist, seq, 
    mm, pos, z, exps, Atempl, Btempl, lcoefs, powers, lhs, rhs, unknowns, i, 
    degsa, degsb, expa, expb, airr, birr, aa, bb, found = {}, llc, ltc, at, 
    bt, Z, A, B, polyeq, lca, lcb, verb, verbp, verbose, solutions, 
    quadratics}, 
    verbose = Verbose /. {opts} /. Options[Hyper];
    solutions = Solutions /. {opts} /. Options[Hyper];
    quadratics = Quadratics /. {opts} /. Options[Hyper];
    eq = Numerator[Together[If[Head[eqn] === Equal, 
      eqn[[1]] - eqn[[2]], eqn]]]; 
    verb = MemberQ[{True, Automatic}, verbose]; 
    {min, max} = 
     ({Min[#1], Max[#1]} & )[First /@ Cases[eq, a[_], Infinity] /. n -> 0]; 
    dd = max - min; unknowns = Table[a[n + i], {i, 0, dd}]; 
    eq = Collect[eq, unknowns]; eq = Factor /@ (eq /. n -> n - min); 
    rhs = -Select[eq, FreeQ[#1, a] & ]; lhs = eq + rhs; 
    polyCoeffs = (Expand[Coefficient[lhs, #1]] & ) /@ unknowns; 
    m = Max[(Exponent[#1, n] & ) /@ polyCoeffs]; lc = Last[polyCoeffs]; 
    tc = First[polyCoeffs]; degseq = (Exponent[#1, n] & ) /@ polyCoeffs; 
    ld = Last[degseq]; td = First[degseq]; powers = Table[z^i, {i, 0, dd}]; 
    lcoefs = (Coeff[#1[[1]], n, #1[[2]]] & ) /@ Thread[{polyCoeffs, degseq}]; 
    llc = Last[lcoefs]; ltc = First[lcoefs]; 
    If[verbose, Print[""]; Print["order = ", dd]; Print["lc = ", lc]; 
      Print["tc = ", tc]; Print["ld = ", ld]; Print["td = ", td]; 
      Print["degree sequence = ", degseq]; 
      Print["leading coeffs = ", lcoefs]]; 
    For[d = -ld, d <= td, d++, 
      seq = degseq + Table[i*d, {i, 0, dd}]; 
      mm = Max[seq]; 
      pos = First /@ Position[seq, mm]; 
      If[verb, Print[""]; Print["d = ", d, "; seq = ", seq, "; pos = ", 
         pos]]; 
      If[Length[pos] <= 1, Null, 
        zeqn = lcoefs[[pos]] . powers[[pos]] == 0; 
        sub = Solve[zeqn, z]; 
        If[!FreeQ[sub, ToRules], 
          Print["\n Warning: algebraic roots of degree > 4 in characteristic \
            polynomial;\n some solutions may not be found"]]; 
        sub = Select[sub, FreeQ[#1, ToRules] & ]; 
        Zlist = Union[If[sub === {}, {}, z /. sub]]; 
        Zlist = Complement[Zlist, {0}]; 
        If[verbose, Print[""]; Print["Zlist = ", Zlist]]; 
        For[ii = 1, ii <= Length[Zlist], ii++, 
          Z = Zlist[[ii]]; 
          If[!FreeQ[Blist[], Blist], 
            bl = If[quadratics === True, 
                FactorList2[lc /. n -> n - dd + 1, n],
                FactorList[lc /. n -> n - dd + 1]]; 
            Blist[] = Select[bl, !FreeQ[#1, n] & ]; 
            If[Blist[] === {}, Blist[] = {{1, 0}}]; 
            birr = First /@ Blist[]; degsb = (Exponent[#1, n] & ) /@ birr; 
            If[verb && Length[Complement[Union[degsb], {0, 1}]] > 0, 
             Print["\n Warning: irreducible factors of degree > 1 in leading \
                    coefficient;\n some solutions may not be found"]]; 	
            expb = Last /@ Blist[]; 
            lcb = 
             (Coefficient[#1[[1]], n, #1[[2]]] & ) /@ Thread[{birr, degsb}]; 
            birr = Expand[birr/lcb]; 
            Blist[] = Thread[{birr, expb}]; 
            Btempl = Bags[expb, degsb]; 
            If[verbose, Print[""]; Print["      Blist = ", Blist[]]]]; 
          If[!FreeQ[Alist[], Alist], 
            al = If[quadratics === True, 
                FactorList2[tc, n], FactorList[tc]]; 
            Alist[] = Select[al, !FreeQ[#1, n] & ]; 
            If[Alist[] === {}, Alist[] = {{1, 0}}]; 
            airr = First /@ Alist[]; 
            degsa = (Exponent[#1, n] & ) /@ airr; 
            If[verb && Length[Complement[Union[degsa], {0, 1}]] > 0, 
             Print["\n Warning: irreducible factors of degree > 1 in trailing \
                    coefficient;\n some solutions may not be found"]]; 
            expa = Last /@ Alist[]; 
            lca = 
             (Coefficient[#1[[1]], n, #1[[2]]] & ) /@ Thread[{airr, degsa}]; 
            airr = Expand[airr/lca]; 
            Alist[] = Thread[{airr, expa}]; 
            Atempl = Bags[expa, degsa]; 
            If[verbose, Print["      Alist = ", Alist[]]]]; 
          For[da = Max[0, d], da <= Min[td, ld + d], da++, 
            db = da - d; 
            If[verb, 
             Print[""]; Print["  da = ", da, ", ", "db = ", db]]; 
            at = Atempl[[da + 1]]; bt = Btempl[[db + 1]]; 
            For[bb = 1, bb <= Length[bt], bb++, 
              B = Apply[Times, birr^bt[[bb]]]; 
              bc = Apply[Times, birr^(expb - bt[[bb]])] /. n -> n + dd - 1; 
              For[aa = 1, aa <= Length[at], aa++, 
                A = Apply[Times, airr^at[[aa]]]; 
                ac = Apply[Times, airr^(expa - at[[aa]])]; 
                If[verb, Print[""]; 
                  Print["     A = ", A, "    B = ", B, "   Z = ", Z]]; 
                polyeq = 
                 Module[{i, j, k}, 
                  Join[{ltc*ac*Product[B /. n -> n + k, {k, 0, dd - 2}]}, 
                   Table[polyCoeffs[[i + 1]]*
                     Product[A /. n -> n + j, {j, i - 1}]*
                     Product[B /. n -> n + k, {k, i, dd - 2}]*Z^i, 
                    {i, dd - 1}], 
                   {llc*bc*Product[A /. n -> n + j, {j, dd - 1}]*Z^dd}]]; 
                polyeq = Expand /@ polyeq; 
                If[verbose, Print[""]; Print["    ", polyeq]]; 
                verbp = Switch[verbose, True, Automatic, _, False]; 
                sol = PolyK[polyeq, 0, a[n], verbp]; 
                If[sol === {}, 
                  If[verbose, Print["     No non-zero polynomial \
                     solution."]; Print[""]], 
                  If[verb, Print[""]; Print[" ===> C = ", sol]]; 
                  sol = 
                   Module[{i}, 
                    Table[Coefficient[sol, C[i]], 
                     {i, Count[{sol}, C[_], Infinity]}]]; 
                  sol = (sol /. n -> n + 1)/sol; 
                  sol = Factor /@ ((Z*A*sol)/B); 
                  If[solutions === 1, Return[sol]];
                  found = Union[found, sol]]]]]]]]; Return[found]];
 
FactorList2[p_, x_] := 
  Module[{pp, lc, lcoefs, kvadr}, 
   pp = FactorList[p]; lc = Select[pp, FreeQ[#1, x] & ]; 
    If[lc === {}, lc = {{1, 1}}]; pp = Complement[pp, lc]; 
    lcoefs = ({LC[First[#1], x], Last[#1]} & ) /@ pp; 
    lc = {Apply[Times, (First[#1]^Last[#1] & ) /@ Join[lc, lcoefs]], 1}; 
    pp = Thread[{Expand[First /@ pp/First /@ lcoefs], Last /@ pp}]; 
    kvadr = Select[pp, Exponent[First[#1], x] === 2 & ]; 
    pp = Complement[pp, kvadr]; kvadr = Apply[Join, (F2[#1, x] & ) /@ kvadr]; 
    Return[Join[{lc}, pp, kvadr]]];
 
LC[p_, x_] := Coefficient[p, x, Exponent[p, x]];
 
F2[{p_, d_}, x_] := ({#1, d} & ) /@ (x - (x /. Solve[p == 0, x]));
 
F2[p_, x_] := Apply[Times, x - (x /. Solve[p == 0, x])];

NaturalZeros[f_, x_] :=
   Module[{fl},
      fl = Select[FactorList[f], !FreeQ[#, x]&];
      fl = Select[First /@ fl, Exponent[#, x] === 1&];
      Select[Join @@ ((x /. Solve[# == 0, x])& /@ fl), 
             IntegerQ[#] && (# >= 0)&]
   ];

dena[num_, den_, x_, ord_] :=      (* Abramov's 94 denominator *)
   Module[{y, A, B, u = 1, res, fl, h, s, d1},
      {A, B} = {num /. x -> x-ord, den};
      res = Resultant[A, B /. x -> x + y, x];
      h = Reverse[Union[NaturalZeros[res, y]]];
      d = PolynomialGCD[A, B /. x -> x+#]& /@ h;
      m = Length[h];
      While[m > 0,
         d1 = First[d];
         s = Product[d1 /. x -> x - j, {j,0,First[h]}];
         u = u s;
         j = 0;
         For[i = 2, i <= m, i++,
           d[[i]] = d[[i]]/PolynomialGCD[d[[i]], s];
           If[Exponent[d[[i]], x] > 0,
             j++; d[[j]] = d[[i]]; h[[j]] = h[[i]]
           ]
         ];
         m = j
      ];
      Return[u]
] ;


RatSol[eqn_, f_[n_], opts___] :=
  Module[{eq, lhs, rhs, args, min, max, ord, coeffs, numer, denom, verbose},

   verbose = Verbose /. {opts} /. Options[RatSol];

(* put everything to the left *)

    eq = Numerator[Together[
      If[Head[eqn] === Equal, 
        eqn[[1]] - eqn[[2]],
        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];

(* find denominator *)

    denom = dena[Last[coeffs], First[coeffs], n, ord];
    If[verbose, Print["denominator = ", denom]];

(* find numerator *)

    coeffs = coeffs/Table[denom /. n->n+i, {i,0,ord}];
    numer = Poly[coeffs.Table[f[n+i], {i,0,ord}] == rhs, f[n], 
      Verbose -> False];
    If[numer === {}, 
      Return[{}],
      Return[Factor /@ Collect[numer/denom, Table[C[i], {i,ord}]]]
    ]
];                



End[];

EndPackage[];