
Clear[Poly];
Poly::usage = "Poly[a[0,n] y[n]+a[1,n] y[n+1] + ... + a[r,n] y[n+r] == c[n], y[n]] finds general polynomial solutions of the given P-finite recurrences (if exists) following Petkovsek's algorithm"
Poly[rec_,y_[n_]]:=Poly[rec,y[n],False];
Poly[rec_,y_[n_],pl_]:=Module[{acoeffs, bcoeffs, order, beta, deg, yvec, cinhom, i1, j1, d, roots, intpos, d1, u, sol},

   (* -- processing input -- *)
   (* determine the order of the recurrence *)
   yvec = Union[Cases[rec[[1]], y[_], Infinity]];
   order = Max[yvec[[All,1]]-n];

   (* retrieve the recurrence coefficients and the right hand side *)
   acoeffs = Factor[Coefficient[rec[[1]],#]&/@yvec];
   cinhom = rec[[2]];

   (* check if input is of the right type *)
   If[Union[PolynomialQ[#,n]&/@Union[acoeffs,{cinhom}]]=!={True}, Print["Input is not a P-finite recurrence"]; Return[$Fail]];

   (* -- step 1: compute the coefficients in the difference-expansion -- *)
   bcoeffs = Table[Sum[Binomial[i1,j1] acoeffs[[i1+1]],{i1,j1,order}],{j1,0,order}]//Factor;
   If[pl, Print[Sum[bcoeffs[[i1+1]] \[Delta][n]^i1**y[n],{i1,0,order}]]];

   (* -- step 2: determine beta -- *)
   beta = Max[Exponent[bcoeffs,n]-Range[0,order]];
   If[pl, Print["\[Beta] = ",beta]];

   (* -- step 3: compute the integer root d1 -- *)
   (* Note: falling-factorial[a,k] = (-1)^k Pochhammer[-a,k] *)
   roots = d/.Solve[Sum[(-1)^i1 Pochhammer[-d,i1] Coefficient[bcoeffs[[i1+1]],n,beta+i1],{i1,0,order}]==0];
   intpos = Position[IntegerQ/@roots,True]//Flatten;

   If[ intpos === {}, d1 = -Infinity, d1 = Max[roots[[intpos]]]];

   (* degree bound *)
   deg = Max[-beta-1,Exponent[cinhom,n]-beta, d1];
   If[pl, Print["D = Max",{-beta-1,Exponent[cinhom,n]-beta,d1}," = ",deg]];

   If[deg<0, Return[0]];

   (* set up an ansatz and solve for the undetermined coefficients *)
   sol = Solve[Thread[CoefficientList[acoeffs.Table[Sum[u[i1] (n+j1)^i1,{i1,0,deg}],{j1,0,order}]-cinhom,n]==0]];

   If[sol==={}, Print["There is no polynomial solution!"]; Return[$Fail]];

   Return[First[Factor[Sum[u[i1] n^i1,{i1,0,deg}]/.sol/.u->C]]];
];