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

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 the associated Polya number of the lattice in question.

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

Command: SeqLimit
Compute the limit of a convergent sequence (Koutschan’s implementation).

(* Given the first values {f[0],...,f[m]} of a sequence f[n] and a basis of its asymptotic solutions, compute the limit Limit[f[n], n->Infinity]. *)
Clear[SeqLimit];
SeqLimit[data_List,asym_,n_]:=
Module[{c,d=Length[asym],pos,ansatz,sol},
pos=Length[data]+Range[-d,-1];
ansatz=Array[c,d].asym;
sol=Solve[((ansatz/.n->#)==data[[#+1]])&/@pos,Array[c,d]][[1]];
Return[N[c[d]/.sol,200]];
];

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

ClearAll[Seq];

Load in advance the REC for (Subscript[(OverTilde[r]), 3,4])(n) in Theorem 4.3 at the end of this file!

Translate the recurrence in term of Ore Polynomials.

RECinS=ToOrePolynomial[REC/.{Seq[k_]->((S[α])^(k-α))}];

Compute the recurrence for the partial Green function: (∑_(0≤n≤(Subscript[n, 0])))(Subscript[(OverTilde[r]), M,N])(n)((((1/((2^M)(({{N}, {M}}))))))^(2n)).

RECPartialGreeninS=DFiniteTimes[{RECinS},Annihilator[((((1/((2^MM)Binomial[NN,MM]))))^(2α)),S[α]]][[1]]**(S[α]-1);

OrePolynomialDegree[RECPartialGreeninS,S[α]]

RECPartialGreen=ApplyOreOperator[RECPartialGreeninS,Seq[α]];

Compute the initial values of the partial Green function by the values of (OverTilde[r]) and then generate a list.

RIni={1,32,6048,2451200,1391236000,921422380032,663895856219904};
PartialGreenIni=Table[Sum[RIni[[i]]*((((1/((2^MM)Binomial[NN,MM]))))^(2(i-1))),{i,1,m}],{m,0,Length@RIni}]

Bound=1000;

PartialGreenList=UnrollRecurrence[RECPartialGreen,Seq[α],PartialGreenIni,Bound];

Analyze the asymptotic behavior of the sequence of partial Green function values.

Asymptotics[RECPartialGreen,Seq[α]]

Compute the limit of partial Green function sequence and the associated Polya number.

lim1=SeqLimit[PartialGreenList,Asymptotics[RECPartialGreen,Seq[α],Order->30],α]

lim2=SeqLimit[PartialGreenList,Asymptotics[RECPartialGreen,Seq[α],Order->32],α]

lim1-lim2

1-1/lim2

Load the REC for (Subscript[(OverTilde[r]), 3,4])(n) in Theorem 4.3.

REC=(221086792032258663383040+3002581182281579476549632 α+18896284453973181469818880 (α^2)+73337056136834742984114176 (α^3)+197017275538043925583364096 (α^4)+389745626428476129286291456 (α^5)+589529476016351811509157888 (α^6)+698690177713813455561031680 (α^7)+659396154092196671988432896 (α^8)+500766687956261350615810048 (α^9)+307887490552535839569608704 (α^10)+153616793330862792246296576 (α^11)+62125104506185984379977728 (α^12)+20265270278609884774662144 (α^13)+5282843409745454510899200 (α^14)+1084193901809507676192768 (α^15)+171154981038855165050880 (α^16)+20040031539432857272320 (α^17)+1638003152561664688128 (α^18)+83373097696100352000 (α^19)+1988330027074191360 (α^20)) Seq[α]+(-123596648884357621088256-1387410081329207115251712 α-7308010505383031273947136 (α^2)-24020604752075269740691456 (α^3)-55262591055735725773815808 (α^4)-94607549345038165436006400 (α^5)-125070786847359746869821440 (α^6)-130760992638503780446109696 (α^7)-109819712522499293630693376 (α^8)-74830049897678615099736064 (α^9)-41599115200046517939601408 (α^10)-18902277196351684209803264 (α^11)-7008965526989775347122176 (α^12)-2109519207312665281560576 (α^13)-510375764108304797663232 (α^14)-97744104267386959429632 (α^15)-14472279363085494386688 (α^16)-1596811738769963089920 (α^17)-123530156260699668480 (α^18)-5975058303292538880 (α^19)-135920997944524800 (α^20)) Seq[1+α]+(2413729498666800513024+25435086835865925058560 α+125542481225411227975680 (α^2)+386097946352750392590336 (α^3)+830183396028360968208384 (α^4)+1327255653860270011465728 (α^5)+1637850112836596110688256 (α^6)+1598197760043557807628288 (α^7)+1252980911862994173739008 (α^8)+797358770338813407952896 (α^9)+414276959391975941603328 (α^10)+176103421096866815410176 (α^11)+61159515859482838548480 (α^12)+17263930413062410149888 (α^13)+3923295133237310914560 (α^14)+706924713366338125824 (α^15)+98652029401005981696 (α^16)+10278087291823325184 (α^17)+752234327699226624 (α^18)+34490272274841600 (α^19)+745214176788480 (α^20)) Seq[2+α]+(-9569617440812835840-97443791378162009856 α-463583339186644316800 (α^2)-1370837922368778354176 (α^3)-2827452328200593850560 (α^4)-4326575055112730856640 (α^5)-5099519612920329528000 (α^6)-4743666552937883189952 (α^7)-3539068890050114722112 (α^8)-2139750587880300657856 (α^9)-1054730779373468537920 (α^10)-424824967934147228480 (α^11)-139643546214642867648 (α^12)-37274084807088072384 (α^13)-8003802897605020608 (α^14)-1361866764260304576 (α^15)-179386646751384192 (α^16)-17635678788631680 (α^17)-1217772669657600 (α^18)-52679537809920 (α^19)-1074030451200 (α^20)) Seq[3+α]+(9051531325562880+90332029095081984 α+420333410362428416 (α^2)+1213206945955473664 (α^3)+2437377188874087136 (α^4)+3625291113645770712 (α^5)+4144688219837114384 (α^6)+3731957019300871994 (α^7)+2689507840271682912 (α^8)+1567534832320365967 (α^9)+743334125295350476 (α^10)+287455002784035524 (α^11)+90539774552500272 (α^12)+23112095925472389 (α^13)+4737102973509780 (α^14)+767930664461310 (α^15)+96195146877576 (α^16)+8977485504456 (α^17)+587451930408 (α^18)+24041253600 (α^19)+462944160 (α^20)) Seq[4+α];