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

NN=5;
MM=4;

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 the even-indexed subsequence (Subscript[(OverTilde[r]), e])(n):=r(2n).

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

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

Timing[ann0EVEN=Annihilator[summandEVEN,{S[k1],S[k2],S[k3],S[k4],S[α]}];]

Timing[ann1EVEN=FindCreativeTelescoping[ann0EVEN,S[k1]-1][[1]];]

Timing[ann2EVEN=FindCreativeTelescoping[ann1EVEN,S[k2]-1][[1]];]

Timing[ann3EVEN=FindCreativeTelescoping[ann2EVEN,S[k3]-1][[1]];]

Timing[ann4EVEN=FindCreativeTelescoping[ann3EVEN,S[k4]-1][[1]];]

Alternatively, you may import the value of {ann1EVEN, ..., ann4EVEN} from an external file.

{ann1EVEN,ann2EVEN,ann3EVEN,ann4EVEN}=ToExpression[Import[NotebookDirectory[]<>"Data-N5M4-Sum-EVEN.txt"]];

ann4EVEN gives a REC for (Subscript[(OverTilde[r]), e])(n).

Apply creative telescoping to the odd-indexed subsequence (Subscript[(OverTilde[r]), o])(n):=r(2n+1).

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

k5=α+((1-NN)/2)-k1-k2-k3-k4;
summandODD=Binomial[2α+1,2k1+1]Binomial[(2α+1)-(2k1+1),2k2+1]Binomial[(2α+1)-(2k1+1)-(2k2+1),2k3+1]Binomial[(2α+1)-(2k1+1)-(2k2+1)-(2k3+1),2k4+1]Binomial[2(α-k1),α-k1]Binomial[2(α-k2),α-k2]Binomial[2(α-k3),α-k3]Binomial[2(α-k4),α-k4]Binomial[2(α-k5),α-k5];

Timing[ann0ODD=Annihilator[summandODD,{S[k1],S[k2],S[k3],S[k4],S[α]}];]

Timing[ann1ODD=FindCreativeTelescoping[ann0ODD,S[k1]-1][[1]];]

Timing[ann2ODD=FindCreativeTelescoping[ann1ODD,S[k2]-1][[1]];]

Timing[ann3ODD=FindCreativeTelescoping[ann2ODD,S[k3]-1][[1]];]

Timing[ann4ODD=FindCreativeTelescoping[ann3ODD,S[k4]-1][[1]];]

Alternatively, you may import the value of {ann1ODD, ..., ann4ODD} from an external file.

{ann1ODD,ann2ODD,ann3ODD,ann4ODD}=ToExpression[Import[NotebookDirectory[]<>"Data-N5M4-Sum-ODD.txt"]];

ann4ODD gives a REC for (Subscript[(OverTilde[r]), o])(n).

Compute the REC for r(n).
REC: Order 12
ODE: Order 83, Degree 12

We first store the RECs for (Subscript[(OverTilde[r]), e])(n) and (Subscript[(OverTilde[r]), o])(n).

RECNormalizedinSEVEN=ann4EVEN[[1]];
RECNormalizedinSOrderEVEN=OrePolynomialDegree[RECNormalizedinSEVEN]

RECNormalizedinSODD=ann4ODD[[1]];
RECNormalizedinSOrderODD=OrePolynomialDegree[RECNormalizedinSODD]

Then we derive the RECs for sequences
{r(0), 0, r(2), 0, ...} and
{0, r(1), 0, r(3), ...},
and compute the REC for their linear combination, including
{r(0), 0, r(2), 0, ...} + {0, r(1), 0, r(3), ...} = {r(0), r(1), r(2), r(3), ...}.

RECNormalizedEVENnew=OrePolynomialSubstitute[{RECNormalizedinSEVEN},{α->(α-0)/2,S[α]->((S[α])^2)}];

RECNormalizedODDnew=OrePolynomialSubstitute[{RECNormalizedinSODD},{α->(α-1)/2,S[α]->((S[α])^2)}];

RECNormalizedinS=DFinitePlus[RECNormalizedEVENnew,RECNormalizedODDnew][[1]];

RECNormalizedinSOrder=OrePolynomialDegree[RECNormalizedinS]

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]]

We also write this REC explicitly.

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

The initial values of 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];

coord=Select[Tuples[Table[i,{i,0,n}],NN],Total[#]==n+(1-NN)/2&];
size=Length@coord;
p=Sum[Multinomial[Sequence@@(2coord[[i]]+1)]*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,2MAX-RECNormalizedinSOrder}]

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

Bound=5000;

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

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

Guess a Minimal ODE for R(z).
ODEGuessinTheta gives the ODE in Theorem 4.6!  (To be displayed at the end of this notebook)
Order 9, Degree 24

ClearAll[Diff];
ODEGuess=GuessMinDE[Take[SeqList,300],Diff[z]];
ODEGuessinD=NormalizeCoefficients[ToOrePolynomial[ODEGuess/.{Derivative[k_][Diff][z]->((Der[z])^k)}/.{Diff[z]->1}]];

ODEGuessinTheta=NormalizeCoefficients[ChangeOreAlgebra[z**ODEGuessinD,OreAlgebra[Euler[z]]]];

ODEGuessinThetaOrder=OrePolynomialDegree[ODEGuessinTheta,Euler[z]]

ODEGuessinThetaDegree=Max[Exponent[OrePolynomialListCoefficients[ODEGuessinTheta],z]]

Get the REC from ODE and write it explicitly.

RECfromODEGuessinS=DFiniteDE2RE[{ODEGuessinD},{z},{α}][[1]];

RECfromODEGuessinSOrder=OrePolynomialDegree[RECfromODEGuessinS,S[α]]

ClearAll[Seq];
SeqfromODEGuess=ApplyOreOperator[RECfromODEGuessinS,Seq[α]];

SeqfromODEGuessList=UnrollRecurrence[SeqfromODEGuess,Seq[α],Take[SeqList,RECfromODEGuessinSOrder],200];

Prove the minimal ODE for R(z).

RECCompare=DFinitePlus[{RECNormalizedinS},{RECfromODEGuessinS}][[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[SeqfromODEGuessList,CheckNum]

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

RECGuess=GuessMinRE[Take[SeqList,300],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=LeadingExponent[RECCompare][[1]]

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

Compute the asymptotics for r(n).

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

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

Approximate the Polya number.

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

Display the ODE in Theorem 4.6

ODEGuessinTheta

Display the REC in Theorem 4.7

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