Multi-headed Lattice Green Function (N = 4, M = 3)

NN=4;
MM=3;

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:
- recurrences (REC) for r(n) or differential equations (ODE) for R(z).

Case 2. M odd or M=N:
- recurrences (REC) for (OverTilde[r])(n) or differential equations (ODE) for (OverTilde[R])(z).

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

Apply creative telescoping to (OverTilde[r])(n).

ClearAll[k1,k2,k3,k4,z,w,α,β];

k4=α-k1-k2-k3;
summand=Binomial[2α,2k1]Binomial[2α-2k1,2k2]Binomial[2α-2k1-2k2,2k3]Binomial[2(α-k1),α-k1]Binomial[2(α-k2),α-k2]Binomial[2(α-k3),α-k3]Binomial[2(α-k4),α-k4];

Timing[ann0=Annihilator[summand,{S[k1],S[k2],S[k3],S[α]}];]

Timing[ann1=FindCreativeTelescoping[ann0,S[k1]-1][[1]];]

Timing[ann2=FindCreativeTelescoping[ann1,S[k2]-1][[1]];]

Timing[ann3=FindCreativeTelescoping[ann2,S[k3]-1][[1]];]

Alternatively, you may import the value of ann3 from an external file.

ann3=ToExpression[Import[NotebookDirectory[]<>"Data-N4M3-Sum.txt"]];

ann3 gives a REC for (OverTilde[r])(n).

Compute the REC for (OverTilde[r])(n).
Order 4

RECNormalizedinS=NormalizeCoefficients[ann3[[1]]];

RECNormalizedinSOrder=OrePolynomialDegree[RECNormalizedinS,S[α]]

We also write this REC explicitly.
SeqNormalized gives the REC in Theorem 4.3!  (To be displayed at the end of this notebook)

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

The initial values of (OverTilde[r])(n) are as follows.

SeqListIni={};

MAX=20;

For[n=0,n≤MAX,n++,
coord=Select[Tuples[Table[i,{i,0,n}],NN],Total[#]==n&];
size=Length@coord;
p=Sum[Multinomial[Sequence@@(2coord[[i]])]*Product[Binomial[2n-2coord[[i,j]],n-coord[[i,j]]],{j,1,NN}],{i,1,size}];
SeqListIni=Append[SeqListIni,p];
];

SeqListIni

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

Now we may numerically verify our REC.

Table[SeqNormalized/.{Seq->seq,α->n},{n,0,MAX-RECNormalizedinSOrder}]

Let us the generate a list of (OverTilde[r])(n).

Bound=5000;

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

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

Get the ODE for (OverTilde[R])(z).
ODENormalizedinD - in terms of the derivation operator D
ODENormalizedinTheta - in terms of the derivation operator θ - Order 24, Degree 4

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

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

ODENormalizedinThetaOrder=OrePolynomialDegree[ODENormalizedinTheta,Euler[w]]

ODENormalizedinThetaDegree=Max[Exponent[OrePolynomialListCoefficients[ODENormalizedinTheta],w]]

Get the ODE for R(z).
ODEinD - in terms of the derivation operator D
ODEinTheta - in terms of the derivation operator θ - Order 24, Degree 8 (Refer to Table 1)

ODEinD=NormalizeCoefficients[DFiniteSubstitute[{ODENormalizedinD},{w->(z^2)},Algebra->OreAlgebra[Der[z]]][[1]]];

ODEinTheta=NormalizeCoefficients[ChangeOreAlgebra[z**ODEinD,OreAlgebra[Euler[z]]]];

ODEinThetaOrder=OrePolynomialDegree[ODEinTheta,Euler[z]]

ODEinThetaDegree=Max[Exponent[OrePolynomialListCoefficients[ODEinTheta],z]]

Guess a Minimal REC for (OverTilde[r])(n).
Its order is 4, and is identical to that of the REC in Theorem 4.3 (RECNormalizedinS).

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

RECGuessinSOrder=OrePolynomialDegree[RECGuessinS,S[α]]

Compute the asymptotics for (OverTilde[r])(n).

AsyList=Asymptotics[SeqNormalized,Seq[α]];
N[AsyList]

Ind=Reverse[Table[Floor[Bound/i],{i,1,3}]]
Table[N[((seq[Ind[[i]]])/(AsyList[[1]]/.{α->Ind[[i]]}))],{i,1,Length@Ind}]
Table[N[((seq[Ind[[i]]])/(AsyList[[2]]/.{α->Ind[[i]]}))],{i,1,Length@Ind}]
Table[N[((seq[Ind[[i]]])/(AsyList[[3]]/.{α->Ind[[i]]}))],{i,1,Length@Ind}]
Table[N[((seq[Ind[[i]]])/(AsyList[[4]]/.{α->Ind[[i]]}))],{i,1,Length@Ind}]

Approximate the Polya number.

AtOne=N[Sum[seq[n]*((((1/((2^MM)Binomial[NN,MM]))))^(2n)),{n,0,Bound}],11]
N[1-(1/AtOne),10]

Display the REC for (OverTilde[r])(n) in Theorem 4.3

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