Multi-headed Lattice Green Function (N = 4, M = 2)
	Find Minimal REC

NN=4;
MM=2;

Recall some basic definitions in the paper:

(Subscript[P, M,N])(z):=(1/(((2 π))^N))(∫_(-π)^π)...(∫_(-π)^π)(1/(1 - (z/((({{N}, {M}})))) (Subscript[σ, M])(cos (Subscript[θ, 1]), ..., cos (Subscript[θ, N]))))d(Subscript[θ, 1])...d(Subscript[θ, N])

(Subscript[R, M,N])(z):=(Subscript[P, M,N])((2^M)(({{N}, {M}}))z)   and  (Subscript[R, M,N])(z)=(∑_(n≥0))(Subscript[r, M,N])(n)(z^n)

Also, for M odd or M=N, we always have r(2n+1)=0. Hence, define
(Subscript[(OverTilde[r]), M,N])(n):=(Subscript[r, M,N])(2n)  and  (Subscript[(OverTilde[R]), M,N])(z):=(∑_(n≥0))(Subscript[(OverTilde[r]), M,N])(n)(z^n)=(∑_(n≥0))(Subscript[r, M,N])(2n)(z^n) 

Our goal is to find:

Case 1. M even and M≠N:
- minimal recurrences (REC) for r(n).

Case 2. M odd or M=N:
- minimal recurrences (REC) for (OverTilde[r])(n).

Command: UnrollRecurrence
Generate a sequence from recurrence & initial values (Koutschan’s implementation).

(* given a recurrence rec in f[n], compute the values {f[0],f[1],...,f[bound]} 
where inits are the initial values {f[0],...,f[d-1]} with d being the order of the recurrence *)
Clear[UnrollRecurrence];
UnrollRecurrence[rec1_,f_[n_],inits_,bound_]:=
Module[{i,x,vals=inits,rec=rec1},
If[Head[rec]=!=Equal,rec=(rec==0)];
rec=rec/.n->n-Max[Cases[rec,f[n+a_.]:>a,Infinity]];
Do[
AppendTo[vals,Solve[rec/.n->i/.f[i]->x/.f[a_]:>vals[[a+1]],x][[1,1,2]]];
,{i,Length[inits],bound}];
Return[vals];
];

Load RISC packages.

<<RISC`HolonomicFunctions`
<<RISC`Asymptotics`
<<RISC`Guess`

HolonomicFunctions Package version 1.7.3 (21-Mar-2017)
written by Christoph Koutschan
Copyright Research Institute for Symbolic Computation (RISC),
Johannes Kepler University, Linz, Austria

Asymptotics Package version 0.3
written by Manuel Kauers
Copyright Research Institute for Symbolic Computation (RISC),
Johannes Kepler University, Linz, Austria

Package GeneratingFunctions version 0.9 written by Christian Mallinger
Copyright Research Institute for Symbolic Computation (RISC),
Johannes Kepler University, Linz, Austria

Guess Package version 0.52
written by Manuel Kauers
Copyright Research Institute for Symbolic Computation (RISC),
Johannes Kepler University, Linz, Austria

We start by importing known ODE for R(z).
Note that the ODE in Koutschan (2013, p. 9, Thm 1) is for P(z)=R(z/(({{N}, {M}}))(2^M)).

ODEDiv2=ToOrePolynomial[12*z*(256+632*z+702*z^2+382*z^3+98*z^4+9*z^5)*P[z]+12*(-384+224*z+3716*z^2+7633*z^3+6734*z^4+2939*z^5+604*z^6+45*z^7)*Derivative[1][P][z]+6*z*(-5376-5248*z+11080*z^2+25286*z^3+19898*z^4+7432*z^5+1286*z^6+81*z^7)*Derivative[2][P][z]+2*z^2*(4+3*z)*(-3456-2304*z+3676*z^2+4920*z^3+2079*z^4+356*z^5+21*z^6)*Derivative[3][P][z]+(-1+z)*z^3*(2+z)*(3+z)*(6+z)*(8+z)*(4+3*z)^2*Derivative[4][P][z]/.{Derivative[k_][P][z]->((Der[z])^k)}/.{P[z]->1}];

Process the data.

Write the ODE in terms of the operators D and θ.

ODENormalizedinD=NormalizeCoefficients[DFiniteSubstitute[{ODEDiv2},{z->w*(2^MM)*Binomial[NN,MM]},Algebra->OreAlgebra[Der[w]]][[1]]];

ODENormalizedinTheta=NormalizeCoefficients[ChangeOreAlgebra[w**ODENormalizedinD,OreAlgebra[Euler[w]]]];

Then transform the above to a REC for r(n) and write it explicitly.

RECNormalizedinS=NormalizeCoefficients[DFiniteDE2RE[{ODENormalizedinD},{w},{α}][[1]]];

RECNormalizedinSOrder=OrePolynomialDegree[RECNormalizedinS,S[α]]

ClearAll[Seq];
SeqNormalized=ApplyOreOperator[RECNormalizedinS,Seq[α]];

Compute the initial values of r(n).

MAX=RECNormalizedinSOrder;
ClearAll[a];

SeriesIni=ApplyOreOperator[ODENormalizedinTheta,Sum[a[n](w^n),{n,0,MAX}]];

SeriesIniSol=Solve[Join[Table[Coefficient[SeriesIni,w,i]==0,{i,1,MAX}],{a[0]==1}],Table[a[i],{i,0,MAX}]]

SeqListIni=Table[SeriesIniSol[[1,k,2]],{k,1,Length@SeriesIniSol[[1]]}]

seq[n_]:=SeqListIni[[n+1]];

Generate a list of r(n).

Bound=200;

SeqList=UnrollRecurrence[SeqNormalized,Seq[α],SeqListIni,Bound];

seq[n_]:=SeqList[[n+1]];

Guess a Minimal REC for r(n).
SeqfromRECGuess gives the REC in Theorem 6.1!  (To be displayed at the end of this notebook)
REC: Order 5
ODE: Order 11, Degree 5

ClearAll[Seq];
RECGuess=GuessMinRE[Take[SeqList,200],Seq[α]];
RECGuessinS=NormalizeCoefficients[ToOrePolynomial[RECGuess/.{Seq[k_]->((S[α])^(k-α))}]];

RECGuessinSOrder=OrePolynomialDegree[RECGuessinS,S[α]]

ODEfromRECGuessinD=NormalizeCoefficients[DFiniteRE2DE[{RECGuessinS},{α},{z}][[1]]];

ODEfromRECGuessinTheta=NormalizeCoefficients[ChangeOreAlgebra[z**ODEfromRECGuessinD,OreAlgebra[Euler[z]]]];

ODEfromRECGuessinThetaOrder=OrePolynomialDegree[ODEfromRECGuessinTheta,Euler[z]]

ODEfromRECGuessinThetaDegree=Max[Exponent[OrePolynomialListCoefficients[ODEfromRECGuessinTheta],z]]

We may also write this REC explicitly.

ClearAll[Seq];
SeqfromRECGuess=ApplyOreOperator[RECGuessinS,Seq[α]];

SeqfromRECGuessList=UnrollRecurrence[SeqfromRECGuess,Seq[α],Take[SeqList,RECGuessinSOrder],200];

Prove the minimal REC for r(n).

RECCompare=DFinitePlus[{RECNormalizedinS},{RECGuessinS}][[1]];

Compute the largest positive integral root of the leading coefficient in the recurrence RECCompare.

LeadCoeff=RECCompare[[1,1,1]];
LeadCoeffRoot=Solve[LeadCoeff==0,α][[All,1,2]]

There are no positive integral roots in our case.

Select[Select[LeadCoeffRoot,IntegerQ],#>0&]

RECCompareOrder=OrePolynomialDegree[RECCompare,S[α]]

CheckNum=RECCompareOrder+20;
Take[SeqList,CheckNum]-Take[SeqfromRECGuessList,CheckNum]

Display the REC in Theorem 6.1

Collect[Expand[-SeqfromRECGuess],Seq[_]]