#####################################################################
#Ising2D.txt: Save this file as  Ising2D.txt                         #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read `IsingD.txt `<Enter>                          #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: Aug. 2016-March 2018
#VERSION OF March 9, 2018

read `Ons48.txt`:
#read CP2to6:


with(numtheory):
with(linalg):


print(`Created:  Aug. 2016-Jan. 2018. VERSION OF March 9, 2018`):
print(` This is Ising2D.txt `):
print(`It is one of the packages that accompany the article `):
print(`"Symbolic Computation to the Aid of Statistical Mechanics `):
print(`by Manuel Kauers and Doron Zeilberger`):
print(` available from Kauers' and Zeilberger's websites`):
print(``):
print(`Please report bugs to zeilberg at math dot rutgers dot edu`):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/  .`):

print(`------------------------------------------------------`):
 print(`For a list of the supporint procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`------------------------------------------------------`):

print(`------------------------------------------------------`):
 print(`For a list of the Checking procedures type ezraC();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`------------------------------------------------------`):

print(`------------------------------------------------------`):
 print(`For a list of the PLOTTING procedures type ezraP();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`------------------------------------------------------`):

print(`------------------------------------------------------`):
 print(`For a list of the Numerical procedures type ezraN();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`------------------------------------------------------`):

print(`------------------------------------------------------`):
 print(`For a list of the  procedures using our NEW, computational approach type ezraNew();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`------------------------------------------------------`):


print(`------------------------------------------------------`):
 print(`For a list of the  procedures from Findrec type ezraF();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`------------------------------------------------------`):


print(`------------------------------------------------------`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`------------------------------------------------------`):

with(combinat):
with(linalg):

ezraNew:=proc()

if args=NULL then
 print(` The  procedures regarding the NEW, computational approach  are: PtorOpeF, vSeriesFEpureNew, vTozNew  `):
 print(``):

else
ezra(args):
fi:

end:

ezraP:=proc()

if args=NULL then
 print(` The  procedures for Plotting are:  PlotAppxOnsFunTable `):
 print(``):

else
ezra(args):
fi:

end:

ezraF:=proc()

if args=NULL then
 print(` The  procedures for guessing recurrences are: findrec, Findrec, SeqFromRec  `):
 print(``):

else
ezra(args):
fi:

end:

ezraC:=proc()

if args=NULL then
 print(` The checking procedures are: CheckAppxOnsFunTable1, Configs, PmnNaive, Wt  `):
 print(``):

else
ezra(args):
fi:

end:

ezraN:=proc()

if args=NULL then
 print(` The numerical procedures are: AppxOnsFun, AppxOnsFunTable, AppxOnsFunTable1, AppxOnsFunTable1Dx, AppxOnsFunTable1Dy, FlToRat, MyMax1, MyMax, Zinn `):
 print(``):

else
ezra(args):
fi:

end:

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are:  Aitken, CheckFunEqEv, EVm, GPw1, KZpc, FunEq, FunEqA, KWt,KWtC, MKa,  OnsFunm01, OnsPolPCt, PnnPC, rectodiff, Tsamek,`):
 print(` V, ZmnPC, ZnnPC `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are:  Cm, GFm, GFmPC, GPw, GuessRecF, IsingPolsPC, IsProdG, LarsDE, LarsFun, LarsRec,`):
 print(` LarsSeries1,  , LarsSeriesNor, LarsSeriesStExact, LarsSeriesSt`):
 print( `Mm,  KZ, KZD, KZf, KZS, KZf, KZfS, IsingPol `):
 print(` Ons, OnsFun, OnsFunDer, OnsFunm0, OnsPolPC, OnsSeries, OnsSeriesLT, `):
 print(` OnsSeriesNor, Pmn , PmnPC, Rm,  vSeriesFE, vSeriesFEpure, vSeriesFEd, vSeriesFEviaLog, `):
 print(`vSeriesST, vSeriesSTf, vSeriesU, vToz, zSeriesFE, zTov, Zmn, ZmnD`):
 print(`  `):


elif nops([args])=1 and op(1,[args])=Aitken then
print(`Aitken(L): The conjectured limit of the sequence L`):
print(`from its three last elements, using Aitken's fromula`):



elif nops([args])=1 and op(1,[args])=AppxOnsFun then
print(`AppxOnsFun(x,y,K): The first K values of the log of the largest eigenvalue of Mm(m,x,x,y) for m from 2 to K. Try:`):
print(`AppxOnsFun(1.5,1,6);`):

elif nops([args])=1 and op(1,[args])=AppxOnsFunTable then
print(`AppxOnsFunTable(xS,xT,yS,yT,res,K):The values of the Approxminate Onsager function AppxOnsFun(x,y,K)[2], `):
print(`from x=xS to x=xT and from y=yS to y=yT, with resolution res. Try:`):
print(`AppxOnsFunTable(1,1.4,1,1.4,0.1,5);`):

elif nops([args])=1 and op(1,[args])=AppxOnsFunTable1 then
print(`AppxOnsFunTable1(xS,xT,res,y0,K):The values of the Approxminate Onsager function AppxOnsFun(x,y,K)[2], `):
print(`for y=y0 from x=xS to x=xT with resolution res. Try:`):
print(`AppxOnsFunTable1(1,2,0.1,1,6);`):

elif nops([args])=1 and op(1,[args])=AppxOnsFunTable1Dx then
print(`AppxOnsFunTable1Dx(xS,xT,res,y0,K,r):The values of the (numerical) r-th (forward, numerical) `):
print(` derivative, w.r.t. x,  of the Approxminate Onsager function AppxOnsFun(x,y,K)[2], `):
print(`for y=y0 from x=xS to x=xT with resolution res. Try:`):
print(`AppxOnsFunTable1Dx(1,2,0.1,1,6,1);`):

elif nops([args])=1 and op(1,[args])=AppxOnsFunTable1Dy then
print(`AppxOnsFunTable1Dy(xS,xT,res,y0,K,r):The values of the (numerical) r-th (forward, numerical) `):
print(` derivative, w.r.t. y,  of the Approxminate Onsager function AppxOnsFun(x,y,K)[2], `):
print(`for y=y0 from x=xS to x=xT with resolution res. Try:`):
print(`AppxOnsFunTable1Dy(1,2,0.1,1,6,1);`):

elif nops([args])=1 and op(1,[args])=CheckAppxOnsFunTable1 then
print(`CheckAppxOnsFunTable1(xS,xT,res,K):Checks AppxOnsFunTable1(xS,xT,res,y0,K) with y0=1 against the known values given by Onsager`):
print(`by taking the ratios, (they should be close to 1).`):
print(`for y=y0 from x=xS to x=xT with resolution res. Try:`):
print(`CheckAppxOnsFunTable1(1,2,0.1,6);`):

elif nops([args])=1 and op(1,[args])=CheckFunEqEv then
print(`CheckFunEqEv(K,m): Checks the functional equation for the largset eigenvalue sof Mm(exp(K),exp(K),1) given in Eq. (30) of the 1941 article`):
print(`by H.A. Kramers and G.H. Wannier, "Statistics of the Two-Dimensional Ferromagent", Part I, Physical Review v. 60 (1941), 252-262. Try: `):
print(`CheckFunEqEv(2.,3);`):

elif nops([args])=1 and op(1,[args])=Cm then
print(`Cm(m,X,x1,x2,y): the characteristic polynomial in X`):
print(`of Mm(m,x1,x2,y): Try:`):
print(`Cm(2,X,x1,x2,y);`):

elif nops([args])=1 and op(1,[args])=Configs then
print(`Configs(m,n): all the m by n {-1,1} matrices, Try:`):
print(`Configs(2,2);`):

elif nops([args])=1 and op(1,[args])=EVm then
print(`EVm(x,y,m): The eigenvalues of Mm(x,x,y,m). Try:`):
print(`EVm(2,1,3);`):

elif nops([args])=1 and op(1,[args])=GFm then
print(`GFm(m,X,x1,x2,y): the generating function, in X, for P_{m,n}(x1,x2,y)`):
print(`Try:`):
print(`GFm(2,X,x1,x2,y);`):

elif nops([args])=1 and op(1,[args])=GFmPC then
print(`GFmPC(m,X,z,w): a pre-computed version of GFm(m,X,z,w), i.e. P_{m,n}(z,z,y)`):
print(`Try:`):
print(`GFmPC(4,X,z,w);`):


elif nops([args])=1 and op(1,[args])=Cm then
print(`Cm(m,X,x1,x2,y): the characteristic polynomial in X`):
print(`of Mm(m,x1,x2,y): Try:`):
print(`Cm(2,X,x1,x2,y);`):

elif nops([args])=1 and op(1,[args])=Configs then
print(`Configs(m,n): all the m by n {-1,1} matrices, Try:`):
print(`Configs(2,2);`):


elif nargs=1 and args[1]=findrec then
print(`findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating`):
print(`the sequence f of degree DEGREE and order ORDER.`):
print(`For example, try: findrec([seq(i,i=1..10)],0,2,n,N);`):

elif nargs=1 and args[1]=Findrec then
print(`Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with`):
print(`poly coffs. ope(n,N), where n is the discrete variable, and N is the shift operator `):
print(`of maximum DEGREE+ORDER<=MaxC`):
print(`e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);`):

elif nops([args])=1 and op(1,[args])=FlToRat then
print(`FlToRat(a): converts the floating point a to a rational number with parameter K (make it big). Try:`):
print(`FlToRat(2.500000000000000000000000001);`):

elif nops([args])=1 and op(1,[args])=FunEq then
print(`FunEq(n,K) verifies the functional equation for Pmn(n,n,x,x,1) as given in Eq. (2.8), chapter 6 of Thompson's book (p. 2.8)`):
print(`for n a multiple of 4 between 4 and 48`):
print(`FunEq(4,10);`):


elif nops([args])=1 and op(1,[args])=FunEqA then
print(`FunEqA(m,n,mu) verifies the functional equation for Pmn(m,n,exp(nu),exp(nu),1) as given in Eq. (2.8), chapter 6 of Thompson's book (p. 2.8).`):
print(`For numeric nu. Try:`):
print(`FunEqA(4,4,0.3);`):

elif nops([args])=1 and op(1,[args])=GFm then
print(`GFm(m,X,x1,x2,y): the generating function, in X, for P_{m,n}(x1,x2,y)`):
print(`Try:`):
print(`GFm(2,X,x1,x2,y);`):

elif nops([args])=1 and op(1,[args])=GPw then
print(`GPw(m0,k,n,d1,K): inputs a positive integer m0, and a positive integer k, and`):
print(`a symbol n, outputs the polynomial in n, of degree <=d1, if it exists, and`):
print(`n0 such that for n>=n0 the coefficient of v^k in ZmnPC(m0,n,v) is given by it `):
print(`it tries up to K terms`):
print(`Try: `):
print(`GPw(3,4,n,1,20):`):


elif nops([args])=1 and op(1,[args])=IsingPol then
print(`IsingPol(k,N): The Ising Polynomial for even-degree polygons with 2*k sides using pre-computed date. Try: `):
print(`IsingPol(2,N); `):

elif nops([args])=1 and op(1,[args])=GuessRecF then
print(`GuessRecF(f,v,P,C): inputs a function f of v and a parameter P, tries to`):
print(`find a recurrence for the  coefficient of complexity<=C`):
print(`Try:`):
print(`GuessRecF( arccosh((1+2*v^2+v^4-2*v*P+2*v^3*P)/(2*v-2*v^3))+log(4*v/(1-v^2)),v,P,50);`):

elif nops([args])=1 and op(1,[args])=IsingPolsPC then
print(`IsingPolsPC(N): `):
print(`the list of (pre-computed (by us!) Ising Polynomials as defined in Colin Thompson's book "Mathematical Statistical Mechanics"`):
print(`and elsewhere, following van der Waerden's combinatorial approach, the k-th term is the number of ways of placing 2*k-bond graphs`):
print(` in a very large torodial lattice with N vertices. Try:`):
print(`IsingPolsPC(N);`):

elif nops([args])=1 and op(1,[args])=IsProdG then
print(`IsProdG(f,t,L,eps): Given a polynomial  t of t`):
print(`decides whether the C-finite sequence of its `):
print(`roots are Cartesian products of sets`):
print(`of cardinalities given by the list L. eps should be small`):
print(`For example, try:`):
print(`IsProdG((1-t)*(1-2*t)*(1-3*t)*(1-6*t),t,[2,2],10^(-10));`):


elif nops([args])=1 and op(1,[args])=KWt then
print(`KWt(m,X,x): the Kramers-Wannier transform the equation whose root is  Lambda(K)/cosh(2*K). Try: `):
print(`KWt(2,X,x); `):

elif nops([args])=1 and op(1,[args])=KWtC then
print(`KWtC(m,x0): Also like CheckFunEqEv(K,m), but does it exactly`):
print(`checks the functional equation for the eigenvaliesof Mm(exp(K),exp(K),1) given in Eq. (30) of the 1941 article`):
print(`by H.A. Kramers and G.H. Wannier, "Statistics of the Two-Dimensional Ferromagent", Part I, Physical Review v. 60 (1941), 252-262. `):
print(`But via the largest root of the characteristic equation where x^2=x0. Try:`):
print(`KWtC(3,3);`):

elif nops([args])=1 and op(1,[args])=KZ then
print(`KZ(k0,m0,N): The Kauers-Zeilberger polynomials in N=n*m0. Inputs  positive integers k0 and m0>2, outputs the`):
print(`polynomial expression in n (valid for n>=k+1) that equals to the coefficient of v^k0 in the series expansion of`):
print(`the partition function of the Ising model for an m0 by n rectangle, after it is divided`):
print(`by ((x+1/x)/2)^(2*n*m0)*2^(n*m0) and v=(x-1/x)/(x+1/x). Try:`):
print(`KZ(4,3,N);`):

elif nops([args])=1 and op(1,[args])=KZD then
print(`KZD(k0,m0,N,d): The generalized Kauers-Zeilberger polynomials in N=m0*n. Inputs  positive integers k0 and m0>2, outputs the`):
print(`polynomial expression in N (valid for N/m0>=k+1) that equals to the coefficient of v^k0 in the series expansion of`):
print(`the d-th derivative of the partition function of the Ising model for an m0 by n rectangle with magnetic field , evaluated at y=1 and divided by d!`):
print(`after it is divided `):
print(`by ((x+1/x)/2)^(2*n*m0)*2^(n*m0) and v=(x-1/x)/(x+1/x). Try:`):
print(` KZD(4,3,N,2); `):


elif nops([args])=1 and op(1,[args])=KZpc then
print(`KZpc(k0,m0,N):  a version of KZ(k0,m0,N) (q.v.) using pre-computed data. Try: `):
print(`KZpc(4,3,N);`):

elif nops([args])=1 and op(1,[args])=KZS then
print(`KZS(K0,m0,N,v): the first K0 terms in the v-expansion for a strip of width m0, via guessing. Try:`):
print(`KZS(10,3,N,v);`):

elif nops([args])=1 and op(1,[args])=KZf then
print(`KZf(k0,m0,N): A direct version of KZ(k0,m0,N), not by guessing.Try:  `):
print(`KZf(4,3,N);`):

elif nops([args])=1 and op(1,[args])=KZfS then
print(`KZfS(K0,m0,N,v): the first K0 terms in the v-expansion for a strip of width m0, via exact generating functions, NOT guessing. Try:`):
print(`KZfS(10,3,N,v);`):



elif nops([args])=1 and op(1,[args])=LarsDE then
print(`LarsDE(u,Du): the Differntial equation of the Onsager function in terms of u=tanh(J/(k*T))^2 alias ((1-x^2)/(1+x^2))^2`):
print(`try:`):
print(`LarsDE(u,Du);`):

elif nops([args])=1 and op(1,[args])=LarsFun then
print(`LarsFun(v,K): The first K terms of the Lars series. Try`):
print(`LarsFun(v,20);`):


elif nops([args])=1 and op(1,[args])=LarsRec then
print(`LarsRec(n,N): the Recurrence equation satisfied by the coefficients of the v-series`):
print(`try: `):
print(`LarsRec(n,N);`):

elif nops([args])=1 and op(1,[args])=MKa then
print(`MKa(x,y,n,v): Manuel Kauers' algorithm for fast computation of Mv. Try: `):
print(`MKa(x,y,2,[1,1,1,1]);`):

elif nops([args])=1 and op(1,[args])=PnnPC then
print(`PnnPC(n,x): a pre-computed Pmn(n,n,x,x,1) when n is a multiple of 4 and <=48. Try: `):
print(` PnnPC(12,x); `):

elif nops([args])=1 and op(1,[args])=PtorOpeF then
print(`PtorOpeF(ope,n,N): solving the recurrence ope(n,N)x(n)=0, x(0)=1`):
print(`in terms of factorials (R is a symbol)`):
print(`For example try PtorOpeF((n+1)*N-1,n,N);`):

elif nops([args])=1 and op(1,[args])=vSeriesU then
print(`vSeriesU(K): the first K terms of the Taylor series in v=tanh(J/(k*T)), of `):
print(`the Onsager power series, for the high-temperature (low nu=J/(k*T)), zero magnetic field, `):
print(` for the internal energy , U, see Eq. (V.5.9) (p. 132) in Colin Thompson's book.`):
print(`Mathematical Statistical Mechanics. USING Onsager's answer. Try:`):
print(`vSeriesU(10); `):

elif nops([args])=1 and op(1,[args])=vSeriesFE then
print(`vSeriesFE(K): the first K terms of the Taylor series in v=tanh(J/(k*T)), of `):
print(`the Onsager power series, for the high-temperature (low nu=J/(k*T)), zero magnetic field, `):
print(` for the negative of the free-energy per spin divided by (kT), i.e. -psi/(k*T), see Eq. (V.5.5) (p. 131) in Colin Thompson's book.`):
print(`Mathematical Statistical Mechanics. USING Onsager's answer. Try:`):
print(`vSeriesFE(10); `):

elif nops([args])=1 and op(1,[args])=vSeriesFEpure then
print(`vSeriesFEpure(K):  Like vSeriesFE(K), but subtracting the trivial part log(1+v^2)-log(1-v^2). Try:`):
print(`vSeriesFEpure(20);`):

elif nops([args])=1 and op(1,[args])=vSeriesFEpureNew then
print(`vSeriesFEpureNew(K):  Like vSeriesFEpure(K), but NOT using Onsager, only doing by pure Symbolic computations. Try: `):
print(`vSeriesFEpureNew(20); `):

elif nops([args])=1 and op(1,[args])=vSeriesFEd then
print(`vSeriesFEd(N,K): the list of lists of first K terms of vSeriesFE for n by n Ising polynomials  for n from 4 to 4*N`):
print(`(N has to be at most 8). Try:`):
print(`vSeriesFEd(2,20);`):



elif nops([args])=1 and op(1,[args])=vSeriesFEviaLog then
print(`vSeriesFEviaLog(K): the first K terms of the Onsager power series but by naively do`):
print(`log(ZnnPC(n,v))/n^2 for the smallest n>=2*K+2 that is a multiple of 4. Try`):
print(`vSeriesFEviaLog(8);`):

elif nops([args])=1 and op(1,[args])=LarsSeries1 then
print(`LarsSeries1(K): the first K terms of the Onsager power series (in v^2), i.e. vSeriesFE(K),  plus log(1-v^2), i.e. the sequence`):
print(`whose k-th term is the  coefficient of N in the polynomial expression for expressing a graph where every vertex has even degree`):
print(`in a sufficiently large torus with N vertices. Try:`):
print(`LarsSeries1(30); `):



elif nops([args])=1 and op(1,[args])=zSeriesFE then
print(`zSeriesFE(K): the first K terms of the z-series, knowing the answer.Try:`):
print(`zSeriesFE(20);`):


elif nops([args])=1 and op(1,[args])=LarsSeriesSt then
print(`LarsSeriesSt(K,m0): the first K terms of the Onsager power series for an m0 by infinity strip. Try: `):
print(`LarsSeriesSt(10,2); `):



elif nops([args])=1 and op(1,[args])=OnsFunSt1 then
print(`OnsFunSt1(m0,v): the m0 by infinity strip Onsager function in terms of v=(x-1/x)/(x+1/x). Try: `):
print(`plot(OnsFunSt1(4,v),v=1.2..3); `):


elif nops([args])=1 and op(1,[args])=LarsSeriesStN then
print(`LarsSeriesStN(K,m0): the first K terms of the Onsager power series for an m0 by infinity strip.`):
print(`Using Numetics. Need to make Digits sufficiently large. To make sure that it workds`):
print(`do it with Digits and 2*Digits and see whether you get the same answer.`):
print(`Try:`):
print(`LarsSeriesStN(10,4); `):

elif nops([args])=1 and op(1,[args])=LarsSeriesNor then
print(`LarsSeriesNor(K): the first K terms of the Onsager power series normalized by multiplying the n-th  term by n. Try:`):
print(`LarsSeriesNor(10); `):

elif nops([args])=1 and op(1,[args])=Mm then
print(`Mm(m,x1,x2,y):The 2^m by 2^m transfer matrix with`):
print(`x1,x2 for vertical and horiz. interaction respectively`):
print(`and y for magnetizaton. Try:`):
print(`Mm(2,x1,x2,y);`):


elif nops([args])=1 and op(1,[args])=MyMax then
print(`MyMax(f,x,x0,x1,K): numerically finds an interval`):
print(`containing the max`):
print(`an expression f in x, between x0 and x1 with`):
print(`resolution 10^(-K), try:`):
print(`MyMax(1-x^2,x,-.1,.1,2);`):

elif nops([args])=1 and op(1,[args])=MyMax1 then
print(`MyMax1(f,x,x0,x1,K): numerically finds an interval`):
print(`containing the max`):
print(`an expression f in x, between x0 and x1 with`):
print(`resolution 10^(-K). It is done the slow way. Try:`):
print(`MyMax1(1-x^2,x,-.1,.1,1);`):

elif nops([args])=1 and op(1,[args])=Ons then
print(`Ons(x): the pre-computed list (thanks to Lundow et. al.)`):
print(`of length 12, whose ith-entry is the Ising polynomial`):
print(`with periodic boundary conditions, of the 4i by 4i torodial`):
print(`square. Try: `):
print(`Ons(x);`):

elif nops([args])=1 and op(1,[args])=OnsFun then
print(`OnsFun(x): the full Onsager solution`):

elif nops([args])=1 and op(1,[args])=OnsFunDer then
print(`OnsFunDer(x0,r): the r-th derivative of the Onsager solution OnsFun(x) w.r.t. to x, at x=x0, done numerically. Try:`):
print(`OnsFunDer(1.2,2);`):

elif nops([args])=1 and op(1,[args])=OnsFunm0 then
print(`OnsFunm0(m0,x): the Onsager solution for an m0 by infinity strip`):

elif nops([args])=1 and op(1,[args])=OnsFun1 then
print(`OnsFun1(x): the Onsager solution, the integral part`):


elif nops([args])=1 and op(1,[args])=OnsPolPC then
print(`OnsPolPC(n,v): the weight-enumerator of n by n {1,-1} matrices under periodic boundary condition in the Ising model `):
print(`divided by (x+1/x)^(2*n^2) when expressed in term of v=(x-1/x)/(x+1/x), using the pre-computed values. Only valid`):
print(`for n a multipe of 4 between 4 and 48. Try:`):
print(`OnsPolPC(4,v);`):

elif nops([args])=1 and op(1,[args])=OnsPolPCt then
print(`OnsPolPCt(n,v,k): the weight-enumerator of n by n {1,-1} matrices under periodic boundary condition in the Ising model `):
print(`divided by (x+1/x)^(2*n^2) when expressed in term of v=(x-1/x)/(x+1/x), truncated up to v^k, using the pre-computed values. Only valid`):
print(`for n a multipe of 4 between 4 and 48. Try: `):
print(`OnsPolPCt(8,v,6); `):

elif nops([args])=1 and op(1,[args])=OnsSeries then
print(`OnsSeries(K): the first K terms of the expansion, in v, of -psi/(k*T) (given by Eq. (5.8), p. 132, in Colin `):
print(` Thompson's book, where v is written nu)`):
print(`Try:`):
print(`OnsSeries(10);`):

elif nops([args])=1 and op(1,[args])=OnsSeriesLT then
print(`OnsSeriesLT(K): the first K terms of the expansion, in v^2 (it is an even function) of the Low-Temperature series`):
print(`see Steven Finch's book "Mathematical Constants", p. 393`):
print(`Try: `):
print(` OnsSeriesLT(10); `):

elif nops([args])=1 and op(1,[args])=OnsSeriesNor then
print(`OnsSeriesNor(K): the first K terms of the expansion, in v^2 (it is an even function)`):
print(`, of -psi/(k*T) times (2*i)!*2^i:`):
print(`Try:`) :
print(`OnsSeriesNor(10);`):


elif nops([args])=1 and op(1,[args])=PlotAppxOnsFunTable then
print(`PlotAppxOnsFunTable(xS,xT,yS,yT,res,K):Plots AppxOnsFunTable(xS,xT,yS,yT,res,K) (q.v.). Try:`):
print(`PlotAppxOnsFunTable(1,1.4,1,1.4,0.1,5);`):

elif nops([args])=1 and op(1,[args])=Pmn then
print(`Pmn(m,n,x1,x2,y): the Ising polynomial for an m by n`):
print(`torodial rectangle.`):
print(`Pmn(2,2,x1,x2,y);`):

elif nops([args])=1 and op(1,[args])=PmnMK then
print(`PmnMK(m,n,x,y): the Ising polynomial for an m by n done via Manuel Kauers' fast algorithm.`):
print(`torodial rectangle.`):
print(`PmnMK(2,2,x,y);`):

elif nops([args])=1 and op(1,[args])=PmnPC then
print(`PmnPC(m,n,z,w): If m<=4 then it returns Pmn(m,n,z,z,w) much faster, otherwise it does Pmn(m,n,z,z,w) (and it takes for ever)`):
print(`Try:`):
print(`PmnPC(4,10,z,w);`):

elif nops([args])=1 and op(1,[args])=PmnBK then
print(`PmnBK(m,n,x): Pmn(m,n,x,x,1) according to Bruria Kaufman's formula for numeric x. Try:`):
print(`PmnBK(3,4,2);`):


elif nops([args])=1 and op(1,[args])=PmnNaive then
print(`PmnNaive(m,n,x1,x2,y): the Ising polynomial for an m by n`):
print(`torodial rectangle, done by direct weight-enumeration.`):
print(`For checking purposes only. Try:`):
print(`PmnNaive(2,2,x1,x2,y);`):

elif nops([args])=1 and op(1,[args])=rectodiff then
print(`rectodiff(ope,n,N,x,D1): inputs a recurrence operator in n and N (the shift in  n) and outputs the differential equation`):
print(`in x and D1 (the diff. w.r.t. x) annihilating the generating function of the sequence satisfying the recurrence. Try:`):
print(`rectodiff(N-n-1,n,N);`):

elif nops([args])=1 and op(1,[args])=Rm then
print(`Rm(m,v,X): GFm(m,X,x,x,1) with the transormation to get the generating function in the v-polynomials. Try:`):
print(`Rm(2,v,X):`):

elif nargs=1 and args[1]=SeqFromRec then
print(`SeqFromRec(ope,n,N,Ini,K): Given the first L-1`):
print(`terms of the sequence Ini=[f(1), ..., f(L-1)]`):
print(`satisfied by the recurrence ope(n,N)f(n)=0`):
print(`extends it to the first K values`):
print(`For example, try:`):
print(`SeqFromRec(N-n-1,n,N,[1],10);`):

elif nops([args])=1 and op(1,[args])=Tsamek then
print(`Tsamek(P,w,m0): inputs a  polynomial P in w, replaces w^m0 by -1. Try:`):
print(`Tsamek(w^7+4*w^2,w,5);`):

elif nops([args])=1 and op(1,[args])=V then
print(`V(i,m):`):
print(`{-1,1} vector corresponding to the integer i between 1 and 2^m`):
print(`Try: V(4,3);`):

elif nops([args])=1 and op(1,[args])=vSeriesST then
print(`vSeriesST(m0,K0): the first K0 terms in the v-series for an m0 by infinity Ising function. VIA GUESSING. Try`):
print(`vSeriesST(3,10);`):

elif nops([args])=1 and op(1,[args])=vSeriesSTf then
print(`vSeriesSTf(m0,K0): the first K0 terms in the v-series for an m0 by infinity Ising function. VIA GENERTING FUNCTION Rm. Try`):
print(`vSeriesSTf(3,10);`):

elif nops([args])=1 and op(1,[args])=Wt then
print(`Wt(C,x1,x2,y): The weight of a configuration C`):
print(`according to the Ising weight, where x1 is for`):
print(`vertical interaction and x2 for horiz. interaction`):
print(`Try: Wt([[1,-1],[-1,1]],x1,x2,y);`):


elif nops([args])=1 and op(1,[args])=vToz then
print(`vToz(K):  Like zSeriesFE(K), but done via the v-series that we can compute from scratch (up to 11 terms that suffice to guess the Onsager solution)`):
print(`Findrec(vToz(11),n,N,4); `):

elif nops([args])=1 and op(1,[args])=vTozNew then
print(`vTozNew(K):  Like zSeriesFE(K), but done via the v-series that we can compute from scratch (up to 11 terms that suffice to guess the Onsager solution)`):
print(`That we acually DID w/o cheating. Try: `):
print(`vTozNew(7); `):


 elif nops([args])=1 and op(1,[args])=Zinn then
print(`Zinn(resh): Zinn-Justin's method to estimate`):
print(`the C,mu, and theta such that`):
print(`resh[i] is appx. Const*mu^i*i^theta`):
print(`For example, try:`):
print(`Zinn([seq(5*i*2^i,i=1..30)]);`):

elif nops([args])=1 and op(1,[args])=Zmn then
print(`Zmn(m,n,v): the expression in v=(x-1/x)/(x+1/x) of Pmn(m,n,x,x,1)/((x+1/x)/2)^N/2^N`):
print(`Try:`):
print(`Zmn(4,4,v); `):

elif nops([args])=1 and op(1,[args])=ZmnD then
print(`ZmnD(m,n,v,d): the expression in v=(x-1/x)/(x+1/x) of diff(Pmn(m,n,x,x,y), y$d)/d!/((x+1/x)/2)^N/2^N `):
print(`Try: `):
print(`ZmnD(4,4,v,d);`):

elif nops([args])=1 and op(1,[args])=ZmnPC then
print(`ZmnPC(m,n,v): the expression in v=(x-1/x)/(x+1/x) of Pmn(m,n,x,x,1)/((x+1/x)/2)^N/2^N, pre-computed`):
print(`Try:`):
print(`ZmnPC(4,4,v); `):

elif nops([args])=1 and op(1,[args])=ZnnPC then
print(`ZnnPC(n,v): the expression in v=(x-1/x)/(x+1/x) of Pmn(n,n,x,x,1)/((x+1/x)/2)^N/2^N (where N=n^2)`):
print(`using the pre-computed Ising polynomial (with zero magnetic field) given by procedure Ons(x)`):
print(`Only works for n multiple of 4 between 4 and 48`):
print(`Try:`):
print(`ZnnPC(4,v);`):

elif nops([args])=1 and op(1,[args])=zTov then
print(`zTov(K):  Like vSeriesFE(K), via the zSeriesFE(K), the reverse of vToz(K). Try: `):
print(`zTov(20);  `):


else
print(`There is no ezra for`,args):
fi:
 
end:

####Start from twoFone
#PtorOpe(ope,n,N,R): solving the recurrence ope(n,N)x(n)=0, x(0)=1
#in terms of rising factorial R(n,k) (R is a symbol)
#For example try PtorOpe((n+1)*N-1,n,N,R);
PtorOpe:=proc(ope,n,N,R)
local rat,c,t1,b1,t1c,b1c,i:

if degree(ope,N)<>1 then
 RETURN(FAIL):
fi:

rat:=normal(-coeff(ope,N,0)/coeff(ope,N,1)):
t1:=expand(numer(rat)):
b1:=expand(denom(rat)):

t1c:=coeff(t1,n,degree(t1,n)):
b1c:=coeff(b1,n,degree(b1,n)):

c:=t1c/b1c:

t1:=expand(t1/t1c):
b1:=expand(b1/b1c):


t1:=factor(t1):
b1:=factor(b1):

t1:=[solve(t1,n)]:
b1:= [solve(b1,n)]:

c^n* convert([seq(R(-t1[i],n),i=1..nops(t1))],`*`)/
convert([seq(R(-b1[i],n),i=1..nops(b1))],`*`):

end:

#PtorOpeF(ope,n,N): solving the recurrence ope(n,N)x(n)=0, x(0)=1
#in terms of factorials (R is a symbol)
#For example try PtorOpeF((n+1)*N-1,n,N);
PtorOpeF:=proc(ope,n,N)
local rat,c,t1,b1,t1c,b1c,i:

if degree(ope,N)<>1 then
 RETURN(FAIL):
fi:

rat:=normal(-coeff(ope,N,0)/coeff(ope,N,1)):
t1:=expand(numer(rat)):
b1:=expand(denom(rat)):

t1c:=coeff(t1,n,degree(t1,n)):
b1c:=coeff(b1,n,degree(b1,n)):

c:=t1c/b1c:

t1:=expand(t1/t1c):
b1:=expand(b1/b1c):


t1:=factor(t1):
b1:=factor(b1):

t1:=[solve(t1,n)]:
b1:= [solve(b1,n)]:

for i from 1 to nops(t1) do
 if type(t1[i],integer) and t1[i]>=0 then
   RETURN(FAIL):
 fi:
od:

for i from 1 to nops(b1) do
 if type(b1[i],integer) and b1[i]>=0 then
   RETURN(FAIL):
 fi:
od:

c^n* 
normal(convert([seq((-t1[i]+n-1)!/(-t1[i]-1)!,i=1..nops(t1))],`*`)/
convert([seq((-b1[i]+n-1)!/(-b1[i]-1)!,i=1..nops(b1))],`*`)):

end:


####end from twoFone
#Tsamek(P,w,m0): inputs a  polynomial P in w, replaces w^m0 by -1. Try:
#Tsamek(w^7+4*w^2,w,5);
Tsamek:=proc(P,w,m0) local i,gu,mu:

gu:=0:

for i from ldegree(P,w) to degree(P,w) do
gu:=gu+coeff(P,w,i)*(-1)^trunc(abs(i)/m0)*w^(i mod m0):
od:


end:

###from SCHUTZENBERGER and also new rectodiff and Zinn
ez:=proc(): print(` xDn(x,D1,n,Ope), rectodiff(ope,n,N,x,D1) `): end:

pashet:=proc(p,n,N)
local i,gu1,gu,p1,ra:
p1:=normal(p):
gu1:=denom(p1):
ra:=degree(gu1,N):
p1:=subs(n=n+ra,numer(p1)):
p1:=expand(p1):
 
gu:=0:
for i from 0 to degree(p1,N) do
gu:=gu+factor(coeff(p1,N,i))*N^i:
od:
RETURN(gu):
end:

difftorec:=proc(ope1,D,x,n,N)
local i,j,gu,ord,ope,r,mu:
 
ope:=expand(ope1):
gu:=0:
ord:=degree(ope,D):
 
for i from 0 to  ord do
 mu:=coeff(ope,D,i):
 
for j from 0 to degree(mu,x) do
  gu:=gu+coeff(mu,x,j)*product(n+r-j,r=1..i)*N^(i-j):
od:
od:
 
 
pashet(gu,n,N):
end:

#MulD(Ope,D1,x): left-multiplies the differentiation operator D1 by the Operator ope in D1 and x. Try
#MulD(x*D1^2,D1,x):
MulD:=proc(Ope,D1,x) expand(D1*Ope+diff(Ope,x)): end:


#xDn(x,D1,n,Ope): the operator (x*D1)^n applied to Ope where D1 is the differentiation operator and n a non-negative integer
xDn:=proc(x,D1,n,Ope) option remember:
if n=0 then
 Ope:
else
expand(x*MulD(xDn(x,D1,n-1,Ope),D1,x)):
fi:
end:



#rectodiff(ope,n,N,x,D1): inputs a recurrence operator in n and N (the shift in  n) and outputs the differential equation
#in x and D1 (the diff. w.r.t. x) annihilating the generating function of the sequence satisfying the recurrence. Try:
#rectodiff(N-n-1,n,N);
rectodiff:=proc(ope,n,N,x,D1) local gu,i,j,lu:
lu:=expand(ope):

gu:=0:
for i from 0 to degree(lu,N) do
for j from 0 to degree(coeff(lu,N,i),n) do
gu:=gu+expand(coeff(coeff(lu,N,i),n,j)*xDn(x,D1,j,x^(-i))):
od:
od:

collect(numer(gu),D1):

end:



sn:=proc(resh,n1):
-1/log(op(n1+1,resh)*op(n1-1,resh)/op(n1,resh)^2):
end:
 
#Zinn(resh): Zinn-Justin's method to estimate
#the C,mu, and theta such that
#resh[i] is appx. Const*mu^i*i^theta
#For example, try:
#Zinn([seq(5*i*2^i,i=1..30)]);
Zinn:=proc(resh)
local s1,s2,theta,mu,n1,i:
if nops({seq(sign(resh[i]),i=1..nops(resh))})<>1 then
 RETURN(FAIL):
fi:

n1:=nops(resh)-1:
s1:=sn(resh,n1):
s2:=sn(resh,n1-1):
theta:=evalf(2*(s1+s2)/(s1-s2)^2):
mu:=evalf(sqrt(op(n1+1,resh)/op(n1-1,resh))*exp(-(s1+s2)/((s1-s2)*s1))):
[theta,mu]:
end:

###End SCHUTZENBERGER and also new rectodiff

###Findrec 
#findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrec([seq(i,i=1..10)],0,2,n,N);
findrecVerbose:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a:
if (1+DEGREE)*(1+ORDER)+3+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+2 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f):
  od:
 
   eq:= eq union {eq1}:
od:
 
var1:=solve(eq,var):
 
kv:={}:
 
for i from 1 to nops(var1) do
 mu:=op(i,var1):
 
if op(1,mu)=op(2,mu) then
   kv:= kv union {op(1,mu)}:
 fi:
 
od:

ope:=subs(var1,ope):

if ope=0 then
  RETURN(FAIL):
fi:

ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}:

if nops(ope)>1 then
print(`There is some slack, there are `, nops(ope)):
print(ope):
RETURN(Yafe(ope[1],N)[2]):
elif nops(ope)=1 then
RETURN(Yafe(ope[1],N)[2]):
else
 RETURN(FAIL):
fi:

end:
 
#Old
#findrecOld(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrec([seq(i,i=1..10)],0,2,n,N);
findrecOld:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a:
option remember:


if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f):
  od:
 
   eq:= eq union {eq1}:
od:
 
var1:=solve(eq,var):
 
kv:={}:
 
for i from 1 to nops(var1) do
 mu:=op(i,var1):
 
if op(1,mu)=op(2,mu) then
   kv:= kv union {op(1,mu)}:
 fi:
 
od:

ope:=subs(var1,ope):

if ope=0 then
  RETURN(FAIL):
fi:

ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}:

if nops(ope)>1 then
RETURN(Yafe(ope[1],N)[2]):
elif nops(ope)=1 then
RETURN(Yafe(ope[1],N)[2]):
else
 RETURN(FAIL):
fi:

end:
 #Old
#findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrec([seq(i,i=1..10)],0,2,n,N);
findrec:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a:
option remember:


if (1+DEGREE)*(1+ORDER)+2+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
RETURN(FAIL):

fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+1 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f):
  od:
 
   eq:= eq union {eq1}:
od:
 
var1:=solve(eq,var):
 
kv:={}:
 
for i from 1 to nops(var1) do
 mu:=op(i,var1):
 
if op(1,mu)=op(2,mu) then
   kv:= kv union {op(1,mu)}:
 fi:
 
od:

ope:=subs(var1,ope):

if ope=0 then
  RETURN(FAIL):
fi:

ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}:

if nops(ope)>1 then
RETURN(Yafe(ope[1],N)[2]):
elif nops(ope)=1 then
RETURN(Yafe(ope[1],N)[2]):
else
 RETURN(FAIL):
fi:

end:
 


Yafe:=proc(ope,N) local i,ope1,coe1,L: 
if ope=0 then
 RETURN(1,0):
fi:
ope1:=expand(ope):
L:=degree(ope1,N):
coe1:=coeff(ope1,N,L):
ope1:=normal(ope1/coe1):
ope1:=normal(ope1):
ope1:=
convert(
[seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`):
factor(coe1),ope1:
end:


#Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with
#poly coffs.
#of maximum DEGREE+ORDER<=MaxC
#e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);
Findrec:=proc(f,n,N,MaxC)
local DEGREE, ORDER,ope,L:

for L from 0 to MaxC do
for ORDER from 0 to L do
 DEGREE:=L-ORDER:
if (1+DEGREE)*(1+ORDER)+2+ORDER>nops(f) then
 print(`Insufficient data for degree`, DEGREE, `and order `,ORDER):
 RETURN(FAIL):
fi:
 ope:=findrec([op(1..(1+DEGREE)*(1+ORDER)+2+ORDER,f)],DEGREE,ORDER,n,N):
     if ope<>FAIL then
       RETURN(ope):
      fi:
 od:
od:
FAIL:

end:





 
#SeqFromRec(ope,n,N,Ini,K): Given the first L-1
#terms of the sequence Ini=[f(1), ..., f(L-1)]
#satisfied by the recurrence ope(n,N)f(n)=0
#extends it to the first K values
SeqFromRec:=proc(ope,n,N,Ini,K)
local ope1,gu,L,n1,j1:
ope1:=Yafe(ope,N)[2]:
L:=degree(ope1,N):
if nops(Ini)<>L then
 ERROR(`Ini should be of length`, L):
fi:

ope1:=expand(subs(n=n-L,ope1)/N^L):

gu:=Ini:

for n1 from nops(Ini)+1 to K do
gu:=[op(gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)),
j1=1..L)]:
od:

gu:

end:
 
#end Findrec



#V(i,m): the {-1,1} vector corresponding to the integer i between 1 and 2^m
V:=proc(i,m)  local gu,i1:
option remember:

if i<1 or i>2^m then
  RETURN(FAIL):
fi:

gu:=V1(i-1,m):
[seq(2*gu[i1]-1,i1=1..m)]:

end:


#V1(i,m): the binary representation
#of i between 0 and 2^m-1 given as a list of 0's and 1's
V1:=proc(i,m)
option remember:

if i<0 or i>=2^m then
  RETURN(FAIL):
fi:

if m=1 then
 if i=0 then
   RETURN([0]):
 elif i=1 then
   RETURN([1]):
 fi:
fi:

if i<2^(m-1) then
 [0,op(V1(i,m-1))]:
else
  [1,op(V1(i-2^(m-1),m-1))]:
fi:
end:

#V2(i,m): the binary representation
#of i between 0 and 2^m-1 given as a list of 0's and 1's
V2:=proc(i,m)
option remember:

if i<0 or i>=2^m then
  RETURN(FAIL):
fi:

if m=1 then
 if i=0 then
   RETURN([-1]):
 elif i=1 then
   RETURN([1]):
 fi:
fi:

if i<2^(m-1) then
 [-1,op(V2(i,m-1))]:
else
  [1,op(V2(i-2^(m-1),m-1))]:
fi:
end:


#Configs(m,n): all the m by n {-1,1} matrices, Try:
#Configs(2,2);
Configs:=proc(m,n) local gu,lu,i,gu1,lu1:
option remember:

if m=0 then
 RETURN({[]}):
fi:

lu:={seq(V(i,m),i=1..2^m)}:

if n=1 then
RETURN({seq([lu[i]],i=1..nops(lu))}):
fi:

gu:=Configs(m,n-1):

{seq(seq([op(gu1),lu1],gu1 in gu),lu1 in lu)}:

end:


#Wt(C,x1,x2,y): The weight of a configuration C
#according to the Ising weight, where x1 is for
#vertical interaction and x2 for horiz. interaction
#Try: Wt([[1,-1],[-1,1]],x1,x2,y);
Wt:=proc(C,x1,x2,y) local m,n,i,j:
m:=nops(C):
n:=nops(C[1]):
x1^(add(add(C[i][j]*C[i][j+1],j=1..n-1)+C[i][n]*C[i][1],i=1..m))*
x2^(add(add(C[i][j]*C[i+1][j],i=1..m-1)+C[m][j]*C[1][j],j=1..n))*
y^add(add(C[i][j],i=1..m),j=1..n):
end:


#PmnNaive(m,n,x1,x2,y): the Ising polynomial for an m by n
#torodial rectangle, done by direct weight-enumeration
#For checking purposes only. Try:
#PmnNaive(2,2,x1,x2,y);
PmnNaive:=proc(m,n,x1,x2,y) local gu,gu1:
gu:=Configs(m,n):
add(Wt(gu1,x1,x2,y), gu1 in gu):
end:


#Mm(m,x1,x2,y):The 2^m by 2^m transfer matrix with
#x1,x2 for vertical and horiz. interaction respectively
#and y for magnetization. Try:
#Mm(2,x1,x2,y);
Mm:=proc(m,x1,x2,y) local i,j,T,v1,v2,gu,i1:
option remember:

for i from 1 to 2^m do
  v1:=V(i,m):
 gu:=y^convert(v1,`+`)*x1^(add(v1[i1]*v1[i1+1],i1=1..m-1)+v1[m]*v1[1]):
 for j from 1 to 2^m do
   v2:=V(j,m):
   T[i,j]:=gu*x2^(add(v1[i1]*v2[i1],i1=1..m)):
 od:
od:

[seq([seq(T[i,j],j=1..2^m)],i=1..2^m)]:

end:

  

#PnnPC(n,x): a pre-computed Pmn(n,n,x,x,1) when n is a multiple of 4 and <=48. Try
#PnnPC(12,x);
PnnPC:=proc(n,x):
if not n mod 4=0 and n>=4 and n<=48 then
 print(`Not available`):
else
 Ons(x)[n/4]:
fi:
end:


#Pmn(m,n,x1,x2,y): the Ising polynomial for an m by n
#torodial rectangle
#Pmn(2,2,x1,x2,y);
Pmn:=proc(m,n,x1,x2,y) local M,i:
M:=Matrix(Mm(m,x1,x2,y)):
M:=M^n:
expand(add(M[i,i],i=1..2^m)):
end:



#PmnPC(m,n,z,w):  Pmn(m,n,z,z,w) 
#using the precomputed generating function for m<=4, otherwise it does Pmn(m,n,z,z,w) (and it takes for ever). Try:
#PmnPC(2,2,z,w,y);
PmnPC:=proc(m,n,z,w) local gu,X:
if m>4 then
 RETURN(Pmn(m,n,z,z,w)):
fi:
gu:=GFmPC(m,X,z,w):
expand(coeff(taylor(gu,X=0,n+1),X,n)):
end:


#Cm(m,X,x1,x2,y): the characteristic polynomial in X
#of Mm(m,x1,x2,y): Try:
#Cm(2,X,x1,x2,y);
Cm:=proc(m,X,x1,x2,y)
LinearAlgebra[CharacteristicPolynomial](Matrix(Mm(m,x1,x2,y)),X):
end:



#GFmPC(m,X,z,w): Pre-computed values of GFm(m,X,z,z,w) (q.v.) for m from 1 to 4
#Try:
#GFmPC(4,X,z,w);
GFmPC:=proc(m,X,z,w) local L:
if not (m>=1 and m<=4) then
 print(`Not yet precomputed, doing it from scratch`):
 RETURN(GFm(m,X,z,z,w)):
fi:

L:=
[-(X*w^2*z^2+X*z^2-2*w)/(X^2*w*z^4-X*w^2*z^2-X^2*w-X*z^2+w), -(2*X^3*w^2*z^16-2
*X^2*w^2*z^16+X^3*w^4*z^12-4*X^2*w^4*z^12-4*X^3*w^2*z^12+3*X*w^4*z^12-2*X^3*w^4
*z^8+X^3*z^12+4*X^2*w^4*z^8-4*X^2*z^12+2*X^3*w^2*z^8+3*X*z^12+X^3*w^4*z^4-2*X^3
*z^8+6*X*w^2*z^8+4*X^2*z^8-4*w^2*z^8+X^3*z^4+2*X^2*w^2)/(X^4*w^2*z^16-2*X^3*w^2
*z^16+X^2*w^2*z^16-X^3*w^4*z^12-4*X^4*w^2*z^12+2*X^2*w^4*z^12+4*X^3*w^2*z^12-X*
w^4*z^12+2*X^3*w^4*z^8-X^3*z^12+6*X^4*w^2*z^8-2*X^2*w^4*z^8+2*X^2*z^12-2*X^3*w^
2*z^8-X*z^12-X^3*w^4*z^4+2*X^3*z^8-2*X*w^2*z^8-4*X^4*w^2*z^4-2*X^2*z^8+w^2*z^8-
X^3*z^4+X^4*w^2-X^2*w^2), -(3*X^5*w^6*z^34+X^5*w^8*z^30+3*X^5*w^4*z^34-3*X^4*w^
7*z^32-18*X^5*w^6*z^30-5*X^4*w^9*z^28-12*X^4*w^5*z^32-6*X^5*w^8*z^26-18*X^5*w^4
*z^30+X^4*w^7*z^28-3*X^4*w^3*z^32+4*X^3*w^10*z^26+9*X^3*w^6*z^30+45*X^5*w^6*z^
26+X^5*w^2*z^30+20*X^4*w^9*z^24+42*X^4*w^5*z^28+16*X^3*w^8*z^26+9*X^3*w^4*z^30+
15*X^5*w^8*z^22+45*X^5*w^4*z^26+20*X^4*w^7*z^24+X^4*w^3*z^28-8*X^3*w^10*z^22-5*
X^3*w^6*z^26-11*X^2*w^9*z^24-6*X^2*w^5*z^28-60*X^5*w^6*z^22-6*X^5*w^2*z^26-30*X
^4*w^9*z^20-56*X^4*w^5*z^24-5*X^4*w*z^28-28*X^3*w^8*z^22-5*X^3*w^4*z^26-11*X^2*
w^7*z^24-20*X^5*w^8*z^18-60*X^5*w^4*z^22-30*X^4*w^7*z^20+20*X^4*w^3*z^24+4*X^3*
w^10*z^18-4*X^3*w^6*z^22+16*X^3*w^2*z^26+11*X^2*w^9*z^20+7*X*w^8*z^22+45*X^5*w^
6*z^18+15*X^5*w^2*z^22+20*X^4*w^9*z^16+42*X^4*w^5*z^20+20*X^4*w*z^24+8*X^3*w^8*
z^18-4*X^3*w^4*z^22+4*X^3*z^26-12*X^2*w^7*z^20-11*X^2*w^3*z^24+15*X^5*w^8*z^14+
45*X^5*w^4*z^18+5*X^4*w^7*z^16-30*X^4*w^3*z^20-6*X^3*w^6*z^18-28*X^3*w^2*z^22-\
20*X^2*w^5*z^20-11*X^2*w*z^24-18*X^5*w^6*z^14-20*X^5*w^2*z^18-5*X^4*w^9*z^12-30
*X^4*w^5*z^16-30*X^4*w*z^20+4*X^3*w^8*z^14-6*X^3*w^4*z^18-8*X^3*z^22+13*X^2*w^7
*z^16-12*X^2*w^3*z^20+13*X*w^6*z^18+7*X*w^2*z^22-6*X^5*w^8*z^10-18*X^5*w^4*z^14
+13*X^4*w^7*z^12+5*X^4*w^3*z^16-3*X^3*w^6*z^14+8*X^3*w^2*z^18+8*X^2*w^5*z^16+11
*X^2*w*z^20+13*X*w^4*z^18+3*X^5*w^6*z^10+15*X^5*w^2*z^14+22*X^4*w^5*z^12+20*X^4
*w*z^16-3*X^3*w^4*z^14+4*X^3*z^18+10*X^2*w^7*z^12+13*X^2*w^3*z^16+8*X*w^6*z^14-\
8*w^5*z^16+X^5*w^8*z^6+3*X^5*w^4*z^10-6*X^4*w^7*z^8+13*X^4*w^3*z^12+11*X^3*w^6*
z^10+4*X^3*w^2*z^14+8*X^2*w^5*z^12+8*X*w^4*z^14-6*X^5*w^2*z^10-12*X^4*w^5*z^8-5
*X^4*w*z^12+11*X^3*w^4*z^10+10*X^2*w^3*z^12-6*X^4*w^3*z^8-2*X^3*w^6*z^6+10*X^2*
w^5*z^8+X^5*w^2*z^6+6*X^4*w^5*z^4-2*X^3*w^4*z^6-2*X^4*w^5)/(X^6*w^5*z^36-2*X^5*
w^6*z^34-9*X^6*w^5*z^32-X^5*w^8*z^30-2*X^5*w^4*z^34+X^4*w^7*z^32+11*X^5*w^6*z^
30+2*X^4*w^9*z^28+4*X^4*w^5*z^32+36*X^6*w^5*z^28+6*X^5*w^8*z^26+11*X^5*w^4*z^30
+X^4*w^3*z^32-X^3*w^10*z^26-2*X^3*w^6*z^30-24*X^5*w^6*z^26-X^5*w^2*z^30-8*X^4*w
^9*z^24-13*X^4*w^5*z^28-4*X^3*w^8*z^26-2*X^3*w^4*z^30-84*X^6*w^5*z^24-15*X^5*w^
8*z^22-24*X^5*w^4*z^26-7*X^4*w^7*z^24+2*X^3*w^10*z^22+X^3*w^6*z^26+2*X^2*w^9*z^
24+X^2*w^5*z^28+25*X^5*w^6*z^22+6*X^5*w^2*z^26+12*X^4*w^9*z^20+15*X^4*w^5*z^24+
2*X^4*w*z^28+7*X^3*w^8*z^22+X^3*w^4*z^26+2*X^2*w^7*z^24+126*X^6*w^5*z^20+20*X^5
*w^8*z^18+25*X^5*w^4*z^22+8*X^4*w^7*z^20-7*X^4*w^3*z^24-X^3*w^10*z^18-4*X^3*w^2
*z^26-2*X^2*w^9*z^20-X*w^8*z^22-10*X^5*w^6*z^18-15*X^5*w^2*z^22-8*X^4*w^9*z^16-\
9*X^4*w^5*z^20-8*X^4*w*z^24-2*X^3*w^8*z^18-X^3*z^26+2*X^2*w^7*z^20+2*X^2*w^3*z^
24-126*X^6*w^5*z^16-15*X^5*w^8*z^14-10*X^5*w^4*z^18+3*X^4*w^7*z^16+8*X^4*w^3*z^
20+3*X^3*w^6*z^18+7*X^3*w^2*z^22+4*X^2*w^5*z^20+2*X^2*w*z^24-3*X^5*w^6*z^14+20*
X^5*w^2*z^18+2*X^4*w^9*z^12+7*X^4*w^5*z^16+12*X^4*w*z^20-X^3*w^8*z^14+3*X^3*w^4
*z^18+2*X^3*z^22-2*X^2*w^7*z^16+2*X^2*w^3*z^20-2*X*w^6*z^18-X*w^2*z^22+84*X^6*w
^5*z^12+6*X^5*w^8*z^10-3*X^5*w^4*z^14-8*X^4*w^7*z^12+3*X^4*w^3*z^16+X^3*w^6*z^
14-2*X^3*w^2*z^18-2*X^2*w^5*z^16-2*X^2*w*z^20-2*X*w^4*z^18+4*X^5*w^6*z^10-15*X^
5*w^2*z^14-7*X^4*w^5*z^12-8*X^4*w*z^16+X^3*w^4*z^14-X^3*z^18-2*X^2*w^7*z^12-2*X
^2*w^3*z^16-X*w^6*z^14+w^5*z^16-36*X^6*w^5*z^8-X^5*w^8*z^6+4*X^5*w^4*z^10+3*X^4
*w^7*z^8-8*X^4*w^3*z^12-4*X^3*w^6*z^10-X^3*w^2*z^14-X^2*w^5*z^12-X*w^4*z^14-X^5
*w^6*z^6+6*X^5*w^2*z^10+5*X^4*w^5*z^8+2*X^4*w*z^12-4*X^3*w^4*z^10-2*X^2*w^3*z^
12+9*X^6*w^5*z^4-X^5*w^4*z^6+3*X^4*w^3*z^8+X^3*w^6*z^6-2*X^2*w^5*z^8-X^5*w^2*z^
6-3*X^4*w^5*z^4+X^3*w^4*z^6-X^6*w^5+X^4*w^5), -(2*X^11*w^10*z^96-2*X^10*w^10*z^
96+4*X^11*w^12*z^92-14*X^10*w^12*z^92+X^11*w^14*z^88-36*X^11*w^10*z^92+10*X^9*w
^12*z^92-12*X^10*w^14*z^88+20*X^10*w^10*z^92-80*X^11*w^12*z^88+4*X^11*w^8*z^92+
25*X^9*w^14*z^88+8*X^9*w^10*z^92-7*X^10*w^16*z^84+210*X^10*w^12*z^88-14*X^10*w^
8*z^92-14*X^8*w^14*z^88-20*X^11*w^14*z^84+300*X^11*w^10*z^88+24*X^9*w^16*z^84-\
78*X^9*w^12*z^88+10*X^9*w^8*z^92+182*X^10*w^14*z^84-80*X^10*w^10*z^88-23*X^8*w^
16*z^84-32*X^8*w^12*z^88+760*X^11*w^12*z^84-80*X^11*w^8*z^88+11*X^9*w^18*z^80-\
272*X^9*w^14*z^84-16*X^9*w^10*z^88+6*X^7*w^16*z^84+119*X^10*w^16*z^80-1459*X^10
*w^12*z^84+210*X^10*w^8*z^88-28*X^8*w^18*z^80+62*X^8*w^14*z^84-48*X^8*w^10*z^88
+190*X^11*w^14*z^80-1520*X^11*w^10*z^84+X^11*w^6*z^88-290*X^9*w^16*z^80+192*X^9
*w^12*z^84-78*X^9*w^8*z^88+17*X^7*w^18*z^80+40*X^7*w^14*z^84-1266*X^10*w^14*z^
80+192*X^10*w^10*z^84-12*X^10*w^6*z^88-5*X^8*w^20*z^76+147*X^8*w^16*z^80+155*X^
8*w^12*z^84-32*X^8*w^8*z^88-4560*X^11*w^12*z^80+760*X^11*w^8*z^84-154*X^9*w^18*
z^76+1308*X^9*w^14*z^80-404*X^9*w^10*z^84+25*X^9*w^6*z^88+12*X^7*w^20*z^76+24*X
^7*w^16*z^80+76*X^7*w^12*z^84-952*X^10*w^16*z^76+6215*X^10*w^12*z^80-1459*X^10*
w^8*z^84+282*X^8*w^18*z^76-4*X^8*w^14*z^80+252*X^8*w^10*z^84-14*X^8*w^6*z^88-7*
X^6*w^20*z^76-16*X^6*w^16*z^80-1140*X^11*w^14*z^76+5130*X^11*w^10*z^80-20*X^11*
w^6*z^84+1558*X^9*w^16*z^76+126*X^9*w^12*z^80+192*X^9*w^8*z^84-66*X^7*w^18*z^76
-142*X^7*w^14*z^80+120*X^7*w^10*z^84+5304*X^10*w^14*z^76-728*X^10*w^10*z^80+182
*X^10*w^6*z^84+55*X^8*w^20*z^72-236*X^8*w^16*z^76-67*X^8*w^12*z^80+155*X^8*w^8*
z^84-46*X^6*w^18*z^76-48*X^6*w^14*z^80+19380*X^11*w^12*z^76-4560*X^11*w^8*z^80+
1001*X^9*w^18*z^72-3684*X^9*w^14*z^76+3064*X^9*w^10*z^80-272*X^9*w^6*z^84-94*X^
7*w^20*z^72-180*X^7*w^16*z^76-198*X^7*w^12*z^80+76*X^7*w^8*z^84+4760*X^10*w^16*
z^72-18096*X^10*w^12*z^76+6215*X^10*w^8*z^80-7*X^10*w^4*z^84-1250*X^8*w^18*z^72
-386*X^8*w^14*z^76-48*X^8*w^10*z^80+62*X^8*w^6*z^84+17*X^6*w^20*z^72-26*X^6*w^
16*z^76-144*X^6*w^12*z^80+4845*X^11*w^14*z^72-11628*X^11*w^10*z^76+190*X^11*w^6
*z^80-4802*X^9*w^16*z^72-1912*X^9*w^12*z^76+126*X^9*w^8*z^80+24*X^9*w^4*z^84-81
*X^7*w^18*z^72+180*X^7*w^14*z^76-566*X^7*w^10*z^80+40*X^7*w^6*z^84+18*X^5*w^20*
z^72+10*X^5*w^16*z^76-14688*X^10*w^14*z^72+4304*X^10*w^10*z^76-1266*X^10*w^6*z^
80-275*X^8*w^20*z^68-550*X^8*w^16*z^72-775*X^8*w^12*z^76-67*X^8*w^8*z^80-23*X^8
*w^4*z^84+166*X^6*w^18*z^72-46*X^6*w^14*z^76-144*X^6*w^10*z^80-62016*X^11*w^12*
z^72+19380*X^11*w^8*z^76-4004*X^9*w^18*z^68+7068*X^9*w^14*z^72-10372*X^9*w^10*z
^76+1308*X^9*w^6*z^80+306*X^7*w^20*z^68-6*X^7*w^16*z^72-624*X^7*w^12*z^76-198*X
^7*w^8*z^80+6*X^7*w^4*z^84+54*X^5*w^18*z^72+56*X^5*w^14*z^76-16660*X^10*w^16*z^
68+37944*X^10*w^12*z^72-18096*X^10*w^8*z^76+119*X^10*w^4*z^80+3146*X^8*w^18*z^
68+698*X^8*w^14*z^72-2992*X^8*w^10*z^76-4*X^8*w^6*z^80+68*X^6*w^20*z^68+88*X^6*
w^16*z^72+417*X^6*w^12*z^76-144*X^6*w^8*z^80-15504*X^11*w^14*z^68+15504*X^11*w^
10*z^72-1140*X^11*w^6*z^76+8918*X^9*w^16*z^68+6008*X^9*w^12*z^72-1912*X^9*w^8*z
^76-290*X^9*w^4*z^80+960*X^7*w^18*z^68-858*X^7*w^14*z^72+996*X^7*w^10*z^76-142*
X^7*w^6*z^80-70*X^5*w^20*z^68+132*X^5*w^16*z^72+116*X^5*w^12*z^76+27336*X^10*w^
14*z^68-17988*X^10*w^10*z^72+5304*X^10*w^6*z^76+825*X^8*w^20*z^64+2670*X^8*w^16
*z^68+435*X^8*w^12*z^72-775*X^8*w^8*z^76+147*X^8*w^4*z^80-84*X^6*w^18*z^68+614*
X^6*w^14*z^72+184*X^6*w^10*z^76-48*X^6*w^6*z^80-11*X^4*w^20*z^68+155040*X^11*w^
12*z^68-62016*X^11*w^8*z^72+11011*X^9*w^18*z^64-11452*X^9*w^14*z^68+19322*X^9*w
^10*z^72-3684*X^9*w^6*z^76+11*X^9*w^2*z^80-504*X^7*w^20*z^64+1454*X^7*w^16*z^68
+3556*X^7*w^12*z^72-624*X^7*w^8*z^76+24*X^7*w^4*z^80-64*X^5*w^18*z^68+60*X^5*w^
14*z^72+168*X^5*w^10*z^76+43316*X^10*w^16*z^64-58548*X^10*w^12*z^68+37944*X^10*
w^8*z^72-952*X^10*w^4*z^76-4730*X^8*w^18*z^64-1010*X^8*w^14*z^68+11764*X^8*w^10
*z^72-386*X^8*w^6*z^76-28*X^8*w^2*z^80-364*X^6*w^20*z^64+787*X^6*w^16*z^68-457*
X^6*w^12*z^72+417*X^6*w^8*z^76-16*X^6*w^4*z^80-62*X^4*w^18*z^68-34*X^4*w^14*z^
72+38760*X^11*w^14*z^64+4845*X^11*w^6*z^72-8554*X^9*w^16*z^64-13878*X^9*w^12*z^
68+6008*X^9*w^8*z^72+1558*X^9*w^4*z^76-2520*X^7*w^18*z^64+4076*X^7*w^14*z^68-\
1542*X^7*w^10*z^72+180*X^7*w^6*z^76+17*X^7*w^2*z^80+80*X^5*w^20*z^64-316*X^5*w^
16*z^68+116*X^5*w^8*z^76-31416*X^10*w^14*z^64+48448*X^10*w^10*z^68-14688*X^10*w
^6*z^72-1650*X^8*w^20*z^60-3652*X^8*w^16*z^64+6489*X^8*w^12*z^68+435*X^8*w^8*z^
72-236*X^8*w^4*z^76-264*X^6*w^18*z^64-610*X^6*w^14*z^68+1000*X^6*w^10*z^72-46*X
^6*w^6*z^76+33*X^4*w^20*z^64-128*X^4*w^16*z^68-64*X^4*w^12*z^72-310080*X^11*w^
12*z^64+155040*X^11*w^8*z^68-22022*X^9*w^18*z^60+20384*X^9*w^14*z^64-15728*X^9*
w^10*z^68+7068*X^9*w^6*z^72-154*X^9*w^2*z^76+336*X^7*w^20*z^60-2630*X^7*w^16*z^
64-7468*X^7*w^12*z^68+3556*X^7*w^8*z^72-180*X^7*w^4*z^76-379*X^5*w^18*z^64-96*X
^5*w^14*z^68-180*X^5*w^10*z^72+56*X^5*w^6*z^76-86632*X^10*w^16*z^60+65892*X^10*
w^12*z^64-58548*X^10*w^8*z^68+4760*X^10*w^4*z^72+3696*X^8*w^18*z^60+5466*X^8*w^
14*z^64-26248*X^8*w^10*z^68+698*X^8*w^6*z^72+282*X^8*w^2*z^76+686*X^6*w^20*z^60
-3003*X^6*w^16*z^64+1097*X^6*w^12*z^68-457*X^6*w^8*z^72-26*X^6*w^4*z^76+114*X^4
*w^18*z^64-118*X^4*w^14*z^68-112*X^4*w^10*z^72-77520*X^11*w^14*z^60-58140*X^11*
w^10*z^64-15504*X^11*w^6*z^68-2002*X^9*w^16*z^60+30330*X^9*w^12*z^64-13878*X^9*
w^8*z^68-4802*X^9*w^4*z^72+3636*X^7*w^18*z^60-8670*X^7*w^14*z^64+4720*X^7*w^10*
z^68-858*X^7*w^6*z^72-66*X^7*w^2*z^76+20*X^5*w^20*z^60-112*X^5*w^16*z^64-322*X^
5*w^12*z^68+10*X^5*w^4*z^76+37*X^3*w^18*z^64+10608*X^10*w^14*z^60-85034*X^10*w^
10*z^64+27336*X^10*w^6*z^68+2310*X^8*w^20*z^56+231*X^8*w^16*z^60-20433*X^8*w^12
*z^64+6489*X^8*w^8*z^68-550*X^8*w^4*z^72-5*X^8*z^76-84*X^6*w^18*z^60-462*X^6*w^
14*z^64-2900*X^6*w^10*z^68+614*X^6*w^6*z^72-46*X^6*w^2*z^76-33*X^4*w^20*z^60+42
*X^4*w^16*z^64-153*X^4*w^12*z^68-64*X^4*w^8*z^72+503880*X^11*w^12*z^60-310080*X
^11*w^8*z^64+33033*X^9*w^18*z^56-38116*X^9*w^14*z^60-17198*X^9*w^10*z^64-11452*
X^9*w^6*z^68+1001*X^9*w^2*z^72+252*X^7*w^20*z^56-558*X^7*w^16*z^60+11138*X^7*w^
12*z^64-7468*X^7*w^8*z^68-6*X^7*w^4*z^72+12*X^7*z^76+974*X^5*w^18*z^60-960*X^5*
w^14*z^64-88*X^5*w^10*z^68+60*X^5*w^6*z^72+70*X^3*w^16*z^64+38*X^3*w^12*z^68+
136136*X^10*w^16*z^56-49980*X^10*w^12*z^60+65892*X^10*w^8*z^64-16660*X^10*w^4*z
^68+396*X^8*w^18*z^56-18898*X^8*w^14*z^60+43912*X^8*w^10*z^64-1010*X^8*w^6*z^68
-1250*X^8*w^2*z^72-658*X^6*w^20*z^56+4244*X^6*w^16*z^60-2425*X^6*w^12*z^64+1097
*X^6*w^8*z^68+88*X^6*w^4*z^72-7*X^6*z^76+54*X^4*w^18*z^60-28*X^4*w^14*z^64+16*X
^4*w^10*z^68-34*X^4*w^6*z^72+125970*X^11*w^14*z^56+167960*X^11*w^10*z^60+38760*
X^11*w^6*z^64+19734*X^9*w^16*z^56-58602*X^9*w^12*z^60+30330*X^9*w^8*z^64+8918*X
^9*w^4*z^68-3738*X^7*w^18*z^56+9768*X^7*w^14*z^60-10752*X^7*w^10*z^64+4076*X^7*
w^6*z^68-81*X^7*w^2*z^72-110*X^5*w^20*z^56+334*X^5*w^16*z^60-188*X^5*w^12*z^64-\
322*X^5*w^8*z^68+132*X^5*w^4*z^72-74*X^3*w^18*z^60+122*X^3*w^14*z^64+24*X^3*w^
10*z^68+37128*X^10*w^14*z^56+87992*X^10*w^10*z^60-31416*X^10*w^6*z^64-2310*X^8*
w^20*z^52+4257*X^8*w^16*z^56+30255*X^8*w^12*z^60-20433*X^8*w^8*z^64+2670*X^8*w^
4*z^68+55*X^8*z^72+1428*X^6*w^18*z^56-172*X^6*w^14*z^60+3980*X^6*w^10*z^64-610*
X^6*w^6*z^68+166*X^6*w^2*z^72+11*X^4*w^20*z^56+474*X^4*w^16*z^60+37*X^4*w^12*z^
64-153*X^4*w^8*z^68-671840*X^11*w^12*z^56+503880*X^11*w^8*z^60-37752*X^9*w^18*z
^52+58630*X^9*w^14*z^56+78628*X^9*w^10*z^60+20384*X^9*w^6*z^64-4004*X^9*w^2*z^
68-756*X^7*w^20*z^52+7410*X^7*w^16*z^56-15330*X^7*w^12*z^60+11138*X^7*w^8*z^64+
1454*X^7*w^4*z^68-94*X^7*z^72-725*X^5*w^18*z^56+1348*X^5*w^14*z^60-142*X^5*w^10
*z^64-96*X^5*w^6*z^68+54*X^5*w^2*z^72+26*X^3*w^16*z^60+32*X^3*w^12*z^64+38*X^3*
w^8*z^68-170170*X^10*w^16*z^52+15028*X^10*w^12*z^56-49980*X^10*w^8*z^60+43316*X
^10*w^4*z^64-4620*X^8*w^18*z^52+35310*X^8*w^14*z^56-60468*X^8*w^10*z^60+5466*X^
8*w^6*z^64+3146*X^8*w^2*z^68+308*X^6*w^20*z^52-3436*X^6*w^16*z^56+2110*X^6*w^12
*z^60-2425*X^6*w^8*z^64+787*X^6*w^4*z^68+17*X^6*z^72-226*X^4*w^18*z^56+490*X^4*
w^14*z^60-252*X^4*w^10*z^64-118*X^4*w^6*z^68-41*X^2*w^16*z^60-167960*X^11*w^14*
z^52-302328*X^11*w^10*z^56-77520*X^11*w^6*z^60-31746*X^9*w^16*z^52+85982*X^9*w^
12*z^56-58602*X^9*w^8*z^60-8554*X^9*w^4*z^64+3672*X^7*w^18*z^52-8210*X^7*w^14*z
^56+14160*X^7*w^10*z^60-8670*X^7*w^6*z^64+960*X^7*w^2*z^68+82*X^5*w^20*z^52+612
*X^5*w^16*z^56-244*X^5*w^12*z^60-188*X^5*w^8*z^64-316*X^5*w^4*z^68+18*X^5*z^72+
37*X^3*w^18*z^56+68*X^3*w^14*z^60+128*X^3*w^10*z^64-87516*X^10*w^14*z^52-12240*
X^10*w^10*z^56+10608*X^10*w^6*z^60+1650*X^8*w^20*z^48-1155*X^8*w^16*z^52-25811*
X^8*w^12*z^56+30255*X^8*w^8*z^60-3652*X^8*w^4*z^64-275*X^8*z^68-2016*X^6*w^18*z
^52+1684*X^6*w^14*z^56-4992*X^6*w^10*z^60-462*X^6*w^6*z^64-84*X^6*w^2*z^68-308*
X^4*w^16*z^56+335*X^4*w^12*z^60+37*X^4*w^8*z^64-128*X^4*w^4*z^68-26*X^2*w^14*z^
60-14*X^2*w^10*z^64+739024*X^11*w^12*z^52-671840*X^11*w^8*z^56+33033*X^9*w^18*z
^48-64636*X^9*w^14*z^52-144274*X^9*w^10*z^56-38116*X^9*w^6*z^60+11011*X^9*w^2*z
^64+744*X^7*w^20*z^48-10758*X^7*w^16*z^52+17402*X^7*w^12*z^56-15330*X^7*w^8*z^
60-2630*X^7*w^4*z^64+306*X^7*z^68-44*X^5*w^18*z^52+256*X^5*w^14*z^56-644*X^5*w^
10*z^60-960*X^5*w^6*z^64-64*X^5*w^2*z^68-202*X^3*w^16*z^56+250*X^3*w^12*z^60+32
*X^3*w^8*z^64+170170*X^10*w^16*z^48+22542*X^10*w^12*z^52+15028*X^10*w^8*z^56-\
86632*X^10*w^4*z^60+5808*X^8*w^18*z^48-40326*X^8*w^14*z^52+67004*X^8*w^10*z^56-\
18898*X^8*w^6*z^60-4730*X^8*w^2*z^64-28*X^6*w^20*z^48+3026*X^6*w^16*z^52-3298*X
^6*w^12*z^56+2110*X^6*w^8*z^60-3003*X^6*w^4*z^64+68*X^6*z^68+144*X^4*w^18*z^52+
130*X^4*w^14*z^56+320*X^4*w^10*z^60-28*X^4*w^6*z^64-62*X^4*w^2*z^68+41*X^2*w^16
*z^56-41*X^2*w^12*z^60+184756*X^11*w^14*z^48+403104*X^11*w^10*z^52+125970*X^11*
w^6*z^56+28314*X^9*w^16*z^48-85358*X^9*w^12*z^52+85982*X^9*w^8*z^56-2002*X^9*w^
4*z^60-3735*X^7*w^18*z^48+9936*X^7*w^14*z^52-13492*X^7*w^10*z^56+9768*X^7*w^6*z
^60-2520*X^7*w^2*z^64-20*X^5*w^20*z^48-772*X^5*w^16*z^52+1888*X^5*w^12*z^56-244
*X^5*w^8*z^60-112*X^5*w^4*z^64-70*X^5*z^68-91*X^3*w^14*z^56+152*X^3*w^10*z^60+
122*X^3*w^6*z^64+106964*X^10*w^14*z^48-140896*X^10*w^10*z^52+37128*X^10*w^6*z^
56-825*X^8*w^20*z^44-10285*X^8*w^16*z^48+11957*X^8*w^12*z^52-25811*X^8*w^8*z^56
+231*X^8*w^4*z^60+825*X^8*z^64+984*X^6*w^18*z^48+1084*X^6*w^14*z^52+4336*X^6*w^
10*z^56-172*X^6*w^6*z^60-264*X^6*w^2*z^64-363*X^4*w^16*z^52+645*X^4*w^12*z^56+
335*X^4*w^8*z^60+42*X^4*w^4*z^64-11*X^4*z^68-98*X^2*w^14*z^56-671840*X^11*w^12*
z^48+739024*X^11*w^8*z^52-22022*X^9*w^18*z^44+46904*X^9*w^14*z^48+184496*X^9*w^
10*z^52+58630*X^9*w^6*z^56-22022*X^9*w^2*z^60-396*X^7*w^20*z^44+7572*X^7*w^16*z
^48-14582*X^7*w^12*z^52+17402*X^7*w^8*z^56-558*X^7*w^4*z^60-504*X^7*z^64+259*X^
5*w^18*z^48-380*X^5*w^14*z^52+1182*X^5*w^10*z^56+1348*X^5*w^6*z^60-379*X^5*w^2*
z^64+46*X^3*w^16*z^52+10*X^3*w^12*z^56+250*X^3*w^8*z^60+70*X^3*w^4*z^64+15*X*w^
14*z^56-136136*X^10*w^16*z^44-45526*X^10*w^12*z^48+22542*X^10*w^8*z^52+136136*X
^10*w^4*z^56-4070*X^8*w^18*z^44+31042*X^8*w^14*z^48-57444*X^8*w^10*z^52+35310*X
^8*w^6*z^56+3696*X^8*w^2*z^60-31*X^6*w^20*z^44-3082*X^6*w^16*z^48+10010*X^6*w^
12*z^52-3298*X^6*w^8*z^56+4244*X^6*w^4*z^60-364*X^6*z^64-24*X^4*w^18*z^48-394*X
^4*w^14*z^52+1300*X^4*w^10*z^56+490*X^4*w^6*z^60+114*X^4*w^2*z^64-91*X^2*w^12*z
^56-41*X^2*w^8*z^60-167960*X^11*w^14*z^44-419900*X^11*w^10*z^48-167960*X^11*w^6
*z^52-14014*X^9*w^16*z^44+45838*X^9*w^12*z^48-85358*X^9*w^8*z^52+19734*X^9*w^4*
z^56+3070*X^7*w^18*z^44-12442*X^7*w^14*z^48+13256*X^7*w^10*z^52-8210*X^7*w^6*z^
56+3636*X^7*w^2*z^60-374*X^5*w^16*z^48-758*X^5*w^12*z^52+1888*X^5*w^8*z^56+334*
X^5*w^4*z^60+80*X^5*z^64-328*X^3*w^14*z^52+26*X^3*w^10*z^56+68*X^3*w^6*z^60+37*
X^3*w^2*z^64-84864*X^10*w^14*z^44+304096*X^10*w^10*z^48-87516*X^10*w^6*z^52+275
*X^8*w^20*z^40+19338*X^8*w^16*z^44+143*X^8*w^12*z^48+11957*X^8*w^8*z^52+4257*X^
8*w^4*z^56-1650*X^8*z^60+114*X^6*w^18*z^44-5708*X^6*w^14*z^48-1016*X^6*w^10*z^
52+1684*X^6*w^6*z^56-84*X^6*w^2*z^60+255*X^4*w^16*z^48-625*X^4*w^12*z^52+645*X^
4*w^8*z^56+474*X^4*w^4*z^60+33*X^4*z^64+28*X^2*w^14*z^52-132*X^2*w^10*z^56-26*X
^2*w^6*z^60+503880*X^11*w^12*z^44-671840*X^11*w^8*z^48+11011*X^9*w^18*z^40-\
18252*X^9*w^14*z^44-183950*X^9*w^10*z^48-64636*X^9*w^6*z^52+33033*X^9*w^2*z^56+
114*X^7*w^20*z^40-2808*X^7*w^16*z^44+13076*X^7*w^12*z^48-14582*X^7*w^8*z^52+
7410*X^7*w^4*z^56+336*X^7*z^60-66*X^5*w^18*z^44-944*X^5*w^14*z^48+2656*X^5*w^10
*z^52+256*X^5*w^6*z^56+974*X^5*w^2*z^60+60*X^3*w^16*z^48-502*X^3*w^12*z^52+10*X
^3*w^8*z^56+26*X^3*w^4*z^60+86632*X^10*w^16*z^40+45968*X^10*w^12*z^44-45526*X^
10*w^8*z^48-170170*X^10*w^4*z^52+1782*X^8*w^18*z^40-19478*X^8*w^14*z^44+42020*X
^8*w^10*z^48-40326*X^8*w^6*z^52+396*X^8*w^2*z^56+9*X^6*w^20*z^40+1578*X^6*w^16*
z^44-13498*X^6*w^12*z^48+10010*X^6*w^8*z^52-3436*X^6*w^4*z^56+686*X^6*z^60-410*
X^4*w^14*z^48-1124*X^4*w^10*z^52+130*X^4*w^6*z^56+54*X^4*w^2*z^60-118*X^2*w^12*
z^52-91*X^2*w^8*z^56-41*X^2*w^4*z^60+125970*X^11*w^14*z^40+348840*X^11*w^10*z^
44+184756*X^11*w^6*z^48+1274*X^9*w^16*z^40+7514*X^9*w^12*z^44+45838*X^9*w^8*z^
48-31746*X^9*w^4*z^52-1653*X^7*w^18*z^40+8740*X^7*w^14*z^44-12378*X^7*w^10*z^48
+9936*X^7*w^6*z^52-3738*X^7*w^2*z^56+612*X^5*w^16*z^44-494*X^5*w^12*z^48-758*X^
5*w^8*z^52+612*X^5*w^4*z^56+20*X^5*z^60+56*X^3*w^14*z^48-416*X^3*w^10*z^52-91*X
^3*w^6*z^56-74*X^3*w^2*z^60+44*X*w^12*z^52+42432*X^10*w^14*z^40-393040*X^10*w^
10*z^44+106964*X^10*w^6*z^48-55*X^8*w^20*z^36-18028*X^8*w^16*z^40-6029*X^8*w^12
*z^44+143*X^8*w^8*z^48-1155*X^8*w^4*z^52+2310*X^8*z^56-266*X^6*w^18*z^40+6106*X
^6*w^14*z^44+1576*X^6*w^10*z^48+1084*X^6*w^6*z^52+1428*X^6*w^2*z^56+85*X^4*w^16
*z^44-595*X^4*w^12*z^48-625*X^4*w^8*z^52-308*X^4*w^4*z^56-33*X^4*z^60+96*X^2*w^
14*z^48-164*X^2*w^10*z^52-98*X^2*w^6*z^56-310080*X^11*w^12*z^40+503880*X^11*w^8
*z^44-4004*X^9*w^18*z^36-1092*X^9*w^14*z^40+147108*X^9*w^10*z^44+46904*X^9*w^6*
z^48-37752*X^9*w^2*z^52-14*X^7*w^20*z^36+534*X^7*w^16*z^40-16884*X^7*w^12*z^44+
13076*X^7*w^8*z^48-10758*X^7*w^4*z^52+252*X^7*z^56-9*X^5*w^18*z^40+700*X^5*w^14
*z^44-4942*X^5*w^10*z^48-380*X^5*w^6*z^52-725*X^5*w^2*z^56-284*X^3*w^12*z^48-\
502*X^3*w^8*z^52-202*X^3*w^4*z^56+28*X*w^10*z^52+15*X*w^6*z^56-43316*X^10*w^16*
z^36-28560*X^10*w^12*z^40+45968*X^10*w^8*z^44+170170*X^10*w^4*z^48-478*X^8*w^18
*z^36+13358*X^8*w^14*z^40-35680*X^8*w^10*z^44+31042*X^8*w^6*z^48-4620*X^8*w^2*z
^52+172*X^6*w^16*z^40+7541*X^6*w^12*z^44-13498*X^6*w^8*z^48+3026*X^6*w^4*z^52-\
658*X^6*z^56+70*X^4*w^14*z^44-708*X^4*w^10*z^48-394*X^4*w^6*z^52-226*X^4*w^2*z^
56+114*X^2*w^12*z^48-118*X^2*w^8*z^52+41*X^2*w^4*z^56-77520*X^11*w^14*z^36-\
232560*X^11*w^10*z^40-167960*X^11*w^6*z^44+3458*X^9*w^16*z^36-39374*X^9*w^12*z^
40+7514*X^9*w^8*z^44+28314*X^9*w^4*z^48+504*X^7*w^18*z^36-1710*X^7*w^14*z^40+
8276*X^7*w^10*z^44-12442*X^7*w^6*z^48+3672*X^7*w^2*z^52-44*X^5*w^16*z^40-1274*X
^5*w^12*z^44-494*X^5*w^8*z^48-772*X^5*w^4*z^52-110*X^5*z^56+164*X^3*w^14*z^44-\
258*X^3*w^10*z^48-328*X^3*w^6*z^52+37*X^3*w^2*z^56+16*X*w^12*z^48+44*X*w^8*z^52
-8568*X^10*w^14*z^36+374952*X^10*w^10*z^40-84864*X^10*w^6*z^44+5*X^8*w^20*z^32+
10204*X^8*w^16*z^36+6181*X^8*w^12*z^40-6029*X^8*w^8*z^44-10285*X^8*w^4*z^48-\
2310*X^8*z^52+68*X^6*w^18*z^36-3730*X^6*w^14*z^40-5528*X^6*w^10*z^44-5708*X^6*w
^6*z^48-2016*X^6*w^2*z^52-57*X^4*w^16*z^40+87*X^4*w^12*z^44-595*X^4*w^8*z^48-\
363*X^4*w^4*z^52+11*X^4*z^56+68*X^2*w^10*z^48+28*X^2*w^6*z^52+155040*X^11*w^12*
z^36-310080*X^11*w^8*z^40+1001*X^9*w^18*z^32+5292*X^9*w^14*z^36-93998*X^9*w^10*
z^40-18252*X^9*w^6*z^44+33033*X^9*w^2*z^48-30*X^7*w^16*z^36+18024*X^7*w^12*z^40
-16884*X^7*w^8*z^44+7572*X^7*w^4*z^48-756*X^7*z^52-268*X^5*w^14*z^40+1188*X^5*w
^10*z^44-944*X^5*w^6*z^48-44*X^5*w^2*z^52+364*X^3*w^12*z^44-284*X^3*w^8*z^48+46
*X^3*w^4*z^52+62*X*w^10*z^48+16660*X^10*w^16*z^32+6732*X^10*w^12*z^36-28560*X^
10*w^8*z^40-136136*X^10*w^4*z^44+70*X^8*w^18*z^32-9302*X^8*w^14*z^36+35080*X^8*
w^10*z^40-19478*X^8*w^6*z^44+5808*X^8*w^2*z^48-413*X^6*w^16*z^36-1237*X^6*w^12*
z^40+7541*X^6*w^8*z^44-3082*X^6*w^4*z^48+308*X^6*z^52+512*X^4*w^14*z^40+172*X^4
*w^10*z^44-410*X^4*w^6*z^48+144*X^4*w^2*z^52+80*X^2*w^12*z^44+114*X^2*w^8*z^48-\
16*w^10*z^48+38760*X^11*w^14*z^32+124032*X^11*w^10*z^36+125970*X^11*w^6*z^40-\
2702*X^9*w^16*z^32+40222*X^9*w^12*z^36-39374*X^9*w^8*z^40-14014*X^9*w^4*z^44-66
*X^7*w^18*z^32-1620*X^7*w^14*z^36-5070*X^7*w^10*z^40+8740*X^7*w^6*z^44-3735*X^7
*w^2*z^48-92*X^5*w^16*z^36+1694*X^5*w^12*z^40-1274*X^5*w^8*z^44-374*X^5*w^4*z^
48+82*X^5*z^52+9*X^3*w^14*z^40-12*X^3*w^10*z^44+56*X^3*w^6*z^48+16*X*w^8*z^48-\
5304*X^10*w^14*z^32-283424*X^10*w^10*z^36+42432*X^10*w^6*z^40-3570*X^8*w^16*z^
32-2793*X^8*w^12*z^36+6181*X^8*w^8*z^40+19338*X^8*w^4*z^44+1650*X^8*z^48+1862*X
^6*w^14*z^36+3768*X^6*w^10*z^40+6106*X^6*w^6*z^44+984*X^6*w^2*z^48+269*X^4*w^12
*z^40+87*X^4*w^8*z^44+255*X^4*w^4*z^48+176*X^2*w^10*z^44+96*X^2*w^6*z^48-62016*
X^11*w^12*z^32+155040*X^11*w^8*z^36-154*X^9*w^18*z^28-2432*X^9*w^14*z^32+48064*
X^9*w^10*z^36-1092*X^9*w^6*z^40-22022*X^9*w^2*z^44-86*X^7*w^16*z^32-11992*X^7*w
^12*z^36+18024*X^7*w^8*z^40-2808*X^7*w^4*z^44+744*X^7*z^48+508*X^5*w^14*z^36+
1474*X^5*w^10*z^40+700*X^5*w^6*z^44+259*X^5*w^2*z^48+40*X^3*w^12*z^40+364*X^3*w
^8*z^44+60*X^3*w^4*z^48-4760*X^10*w^16*z^28+7140*X^10*w^12*z^32+6732*X^10*w^8*z
^36+86632*X^10*w^4*z^40-4*X^8*w^18*z^28+4478*X^8*w^14*z^32-27232*X^8*w^10*z^36+
13358*X^8*w^6*z^40-4070*X^8*w^2*z^44+61*X^6*w^16*z^32+381*X^6*w^12*z^36-1237*X^
6*w^8*z^40+1578*X^6*w^4*z^44-28*X^6*z^48-208*X^4*w^14*z^36+364*X^4*w^10*z^40+70
*X^4*w^6*z^44-24*X^4*w^2*z^48+56*X^2*w^12*z^40+80*X^2*w^8*z^44-15504*X^11*w^14*
z^28-52326*X^11*w^10*z^32-77520*X^11*w^6*z^36+1018*X^9*w^16*z^28-25698*X^9*w^12
*z^32+40222*X^9*w^8*z^36+1274*X^9*w^4*z^40+1222*X^7*w^14*z^32+4368*X^7*w^10*z^
36-1710*X^7*w^6*z^40+3070*X^7*w^2*z^44+10*X^5*w^16*z^32-344*X^5*w^12*z^36+1694*
X^5*w^8*z^40+612*X^5*w^4*z^44-20*X^5*z^48+304*X^3*w^10*z^40+164*X^3*w^6*z^44+
5712*X^10*w^14*z^28+176426*X^10*w^10*z^32-8568*X^10*w^6*z^36+715*X^8*w^16*z^28-\
39*X^8*w^12*z^32-2793*X^8*w^8*z^36-18028*X^8*w^4*z^40-825*X^8*z^44-598*X^6*w^14
*z^32+1452*X^6*w^10*z^36-3730*X^6*w^6*z^40+114*X^6*w^2*z^44+144*X^4*w^12*z^36+
269*X^4*w^8*z^40+85*X^4*w^4*z^44+52*X^2*w^10*z^40+19380*X^11*w^12*z^28-62016*X^
11*w^8*z^32+11*X^9*w^18*z^24+148*X^9*w^14*z^28-20890*X^9*w^10*z^32+5292*X^9*w^6
*z^36+11011*X^9*w^2*z^40+74*X^7*w^16*z^28+5026*X^7*w^12*z^32-11992*X^7*w^8*z^36
+534*X^7*w^4*z^40-396*X^7*z^44-320*X^5*w^14*z^32-880*X^5*w^10*z^36-268*X^5*w^6*
z^40-66*X^5*w^2*z^44+66*X^3*w^12*z^36+40*X^3*w^8*z^40+952*X^10*w^16*z^24-9894*X
^10*w^12*z^28+7140*X^10*w^8*z^32-43316*X^10*w^4*z^36-918*X^8*w^14*z^28+12496*X^
8*w^10*z^32-9302*X^8*w^6*z^36+1782*X^8*w^2*z^40+20*X^6*w^16*z^28-1053*X^6*w^12*
z^32+381*X^6*w^8*z^36+172*X^6*w^4*z^40-31*X^6*z^44-10*X^4*w^14*z^32+112*X^4*w^
10*z^36+512*X^4*w^6*z^40+56*X^2*w^8*z^40+4845*X^11*w^14*z^24+17100*X^11*w^10*z^
28+38760*X^11*w^6*z^32-206*X^9*w^16*z^24+12120*X^9*w^12*z^28-25698*X^9*w^8*z^32
+3458*X^9*w^4*z^36-384*X^7*w^14*z^28-2764*X^7*w^10*z^32-1620*X^7*w^6*z^36-1653*
X^7*w^2*z^40-22*X^5*w^12*z^32-344*X^5*w^8*z^36-44*X^5*w^4*z^40+80*X^3*w^10*z^36
+9*X^3*w^6*z^40-2652*X^10*w^14*z^24-92716*X^10*w^10*z^28-5304*X^10*w^6*z^32-63*
X^8*w^16*z^24+225*X^8*w^12*z^28-39*X^8*w^8*z^32+10204*X^8*w^4*z^36+275*X^8*z^40
-32*X^6*w^14*z^28-2484*X^6*w^10*z^32+1862*X^6*w^6*z^36-266*X^6*w^2*z^40-12*X^4*
w^12*z^32+144*X^4*w^8*z^36-57*X^4*w^4*z^40-4560*X^11*w^12*z^24+19380*X^11*w^8*z
^28+265*X^9*w^14*z^24+8948*X^9*w^10*z^28-2432*X^9*w^6*z^32-4004*X^9*w^2*z^36-18
*X^7*w^16*z^24-1658*X^7*w^12*z^28+5026*X^7*w^8*z^32-30*X^7*w^4*z^36+114*X^7*z^
40+40*X^5*w^14*z^28+364*X^5*w^10*z^32+508*X^5*w^6*z^36-9*X^5*w^2*z^40-14*X^3*w^
12*z^32+66*X^3*w^8*z^36-119*X^10*w^16*z^20+6474*X^10*w^12*z^24-9894*X^10*w^8*z^
28+16660*X^10*w^4*z^32-200*X^8*w^14*z^24-1204*X^8*w^10*z^28+4478*X^8*w^6*z^32-\
478*X^8*w^2*z^36+876*X^6*w^12*z^28-1053*X^6*w^8*z^32-413*X^6*w^4*z^36+9*X^6*z^
40-84*X^4*w^10*z^32-208*X^4*w^6*z^36-1140*X^11*w^14*z^20-4180*X^11*w^10*z^24-\
15504*X^11*w^6*z^28+18*X^9*w^16*z^20-4536*X^9*w^12*z^24+12120*X^9*w^8*z^28-2702
*X^9*w^4*z^32+90*X^7*w^14*z^24+408*X^7*w^10*z^28+1222*X^7*w^6*z^32+504*X^7*w^2*
z^36-54*X^5*w^12*z^28-22*X^5*w^8*z^32-92*X^5*w^4*z^36-28*X^3*w^10*z^32+726*X^10
*w^14*z^20+41664*X^10*w^10*z^24+5712*X^10*w^6*z^28+507*X^8*w^12*z^24+225*X^8*w^
8*z^28-3570*X^8*w^4*z^32-55*X^8*z^36+56*X^6*w^14*z^24+704*X^6*w^10*z^28-598*X^6
*w^6*z^32+68*X^6*w^2*z^36-68*X^4*w^12*z^28-12*X^4*w^8*z^32+14*X^2*w^10*z^32+760
*X^11*w^12*z^20-4560*X^11*w^8*z^24-100*X^9*w^14*z^20-4006*X^9*w^10*z^24+148*X^9
*w^6*z^28+1001*X^9*w^2*z^32+570*X^7*w^12*z^24-1658*X^7*w^8*z^28-86*X^7*w^4*z^32
-14*X^7*z^36-160*X^5*w^10*z^28-320*X^5*w^6*z^32-14*X^3*w^8*z^32+7*X^10*w^16*z^
16-2675*X^10*w^12*z^20+6474*X^10*w^8*z^24-4760*X^10*w^4*z^28+144*X^8*w^14*z^20-\
2484*X^8*w^10*z^24-918*X^8*w^6*z^28+70*X^8*w^2*z^32-416*X^6*w^12*z^24+876*X^6*w
^8*z^28+61*X^6*w^4*z^32+40*X^4*w^10*z^28-10*X^4*w^6*z^32+190*X^11*w^14*z^16+720
*X^11*w^10*z^20+4845*X^11*w^6*z^24+1350*X^9*w^12*z^20-4536*X^9*w^8*z^24+1018*X^
9*w^4*z^28-16*X^7*w^14*z^20+456*X^7*w^10*z^24-384*X^7*w^6*z^28-66*X^7*w^2*z^32+
2*X^5*w^12*z^24-54*X^5*w^8*z^28+10*X^5*w^4*z^32-114*X^10*w^14*z^16-15968*X^10*w
^10*z^20-2652*X^10*w^6*z^24-544*X^8*w^12*z^20+507*X^8*w^8*z^24+715*X^8*w^4*z^28
+5*X^8*z^32+208*X^6*w^10*z^24-32*X^6*w^6*z^28-68*X^4*w^8*z^28-80*X^11*w^12*z^16
+760*X^11*w^8*z^20+12*X^9*w^14*z^16+1556*X^9*w^10*z^20+265*X^9*w^6*z^24-154*X^9
*w^2*z^28-146*X^7*w^12*z^20+570*X^7*w^8*z^24+74*X^7*w^4*z^28+4*X^5*w^10*z^24+40
*X^5*w^6*z^28+711*X^10*w^12*z^16-2675*X^10*w^8*z^20+952*X^10*w^4*z^24-22*X^8*w^
14*z^16+1936*X^8*w^10*z^20-200*X^8*w^6*z^24-4*X^8*w^2*z^28+96*X^6*w^12*z^20-416
*X^6*w^8*z^24+20*X^6*w^4*z^28-44*X^4*w^10*z^24-20*X^11*w^14*z^12-78*X^11*w^10*z
^16-1140*X^11*w^6*z^20-292*X^9*w^12*z^16+1350*X^9*w^8*z^20-206*X^9*w^4*z^24-224
*X^7*w^10*z^20+90*X^7*w^6*z^24+2*X^5*w^8*z^24+8*X^10*w^14*z^12+5112*X^10*w^10*z
^16+726*X^10*w^6*z^20+204*X^8*w^12*z^16-544*X^8*w^8*z^20-63*X^8*w^4*z^24-192*X^
6*w^10*z^20+56*X^6*w^6*z^24+4*X^11*w^12*z^12-80*X^11*w^8*z^16-418*X^9*w^10*z^16
-100*X^9*w^6*z^20+11*X^9*w^2*z^24+14*X^7*w^12*z^16-146*X^7*w^8*z^20-18*X^7*w^4*
z^24-112*X^10*w^12*z^12+711*X^10*w^8*z^16-119*X^10*w^4*z^20-780*X^8*w^10*z^16+
144*X^8*w^6*z^20+96*X^6*w^8*z^20+X^11*w^14*z^8+4*X^11*w^10*z^12+190*X^11*w^6*z^
16+38*X^9*w^12*z^12-292*X^9*w^8*z^16+18*X^9*w^4*z^20+28*X^7*w^10*z^16-16*X^7*w^
6*z^20-1312*X^10*w^10*z^12-114*X^10*w^6*z^16-28*X^8*w^12*z^12+204*X^8*w^8*z^16+
48*X^6*w^10*z^16+4*X^11*w^8*z^12+64*X^9*w^10*z^12+12*X^9*w^6*z^16+14*X^7*w^8*z^
16+8*X^10*w^12*z^8-112*X^10*w^8*z^12+7*X^10*w^4*z^16+184*X^8*w^10*z^12-22*X^8*w
^6*z^16-20*X^11*w^6*z^12-2*X^9*w^12*z^8+38*X^9*w^8*z^12+252*X^10*w^10*z^8+8*X^
10*w^6*z^12-28*X^8*w^8*z^12-4*X^9*w^10*z^8+8*X^10*w^8*z^8-20*X^8*w^10*z^8+X^11*
w^6*z^8-2*X^9*w^8*z^8-32*X^10*w^10*z^4+2*X^10*w^10)/(X^9*w^8*z^72-2*X^8*w^8*z^
72+X^7*w^8*z^72-2*X^8*w^10*z^68-18*X^9*w^8*z^68+4*X^7*w^10*z^68-X^8*w^12*z^64+
27*X^8*w^8*z^68-2*X^6*w^10*z^68+3*X^7*w^12*z^64-8*X^7*w^8*z^68+26*X^8*w^10*z^64
-2*X^8*w^6*z^68-3*X^6*w^12*z^64-X^6*w^8*z^68+153*X^9*w^8*z^64+2*X^7*w^14*z^60-\
38*X^7*w^10*z^64+4*X^7*w^6*z^68+X^5*w^12*z^64+14*X^8*w^12*z^60-170*X^8*w^8*z^64
-4*X^6*w^14*z^60+10*X^6*w^10*z^64-2*X^6*w^6*z^68-27*X^7*w^12*z^60+33*X^7*w^8*z^
64+2*X^5*w^14*z^60+2*X^5*w^10*z^64-154*X^8*w^10*z^60+26*X^8*w^6*z^64-X^6*w^16*z
^56+15*X^6*w^12*z^60-4*X^6*w^8*z^64-816*X^9*w^8*z^60-22*X^7*w^14*z^56+164*X^7*w
^10*z^60-38*X^7*w^6*z^64+2*X^5*w^16*z^56-X^5*w^12*z^60+4*X^5*w^8*z^64-91*X^8*w^
12*z^56+664*X^8*w^8*z^60-X^8*w^4*z^64+30*X^6*w^14*z^56-24*X^6*w^10*z^60+10*X^6*
w^6*z^64-X^4*w^16*z^56-X^4*w^12*z^60+102*X^7*w^12*z^56-108*X^7*w^8*z^60+3*X^7*w
^4*z^64-6*X^5*w^14*z^56+2*X^5*w^6*z^64+546*X^8*w^10*z^56-154*X^8*w^6*z^60+8*X^6
*w^16*z^52-23*X^6*w^12*z^56+48*X^6*w^8*z^60-3*X^6*w^4*z^64-2*X^4*w^14*z^56-2*X^
4*w^10*z^60+3060*X^9*w^8*z^56+110*X^7*w^14*z^52-434*X^7*w^10*z^56+164*X^7*w^6*z
^60-11*X^5*w^16*z^52+2*X^5*w^12*z^56-8*X^5*w^8*z^60+X^5*w^4*z^64+364*X^8*w^12*z
^52-1806*X^8*w^8*z^56+14*X^8*w^4*z^60-98*X^6*w^14*z^52+54*X^6*w^10*z^56-24*X^6*
w^6*z^60+2*X^4*w^16*z^52-2*X^4*w^12*z^56-4*X^4*w^8*z^60-198*X^7*w^12*z^52+304*X
^7*w^8*z^56-27*X^7*w^4*z^60+4*X^5*w^14*z^52-12*X^5*w^10*z^56+X^3*w^16*z^52-1274
*X^8*w^10*z^52+546*X^8*w^6*z^56-28*X^6*w^16*z^48-3*X^6*w^12*z^52-140*X^6*w^8*z^
56+15*X^6*w^4*z^60-2*X^4*w^14*z^52-4*X^4*w^10*z^56-2*X^4*w^6*z^60-8568*X^9*w^8*
z^52-330*X^7*w^14*z^48+818*X^7*w^10*z^52-434*X^7*w^6*z^56+2*X^7*w^2*z^60+24*X^5
*w^16*z^48-18*X^5*w^12*z^52-4*X^5*w^8*z^56-X^5*w^4*z^60+2*X^3*w^14*z^52+2*X^3*w
^10*z^56-1001*X^8*w^12*z^48+3640*X^8*w^8*z^52-91*X^8*w^4*z^56+184*X^6*w^14*z^48
-126*X^6*w^10*z^52+54*X^6*w^6*z^56-4*X^6*w^2*z^60+2*X^4*w^16*z^48-X^4*w^12*z^52
+4*X^4*w^8*z^56-X^4*w^4*z^60+165*X^7*w^12*z^48-646*X^7*w^8*z^52+102*X^7*w^4*z^
56-4*X^5*w^14*z^48+6*X^5*w^10*z^52-12*X^5*w^6*z^56+2*X^5*w^2*z^60-2*X^3*w^16*z^
48+4*X^3*w^12*z^52+X^3*w^8*z^56+2002*X^8*w^10*z^48-1274*X^8*w^6*z^52+56*X^6*w^
16*z^44+44*X^6*w^12*z^48+188*X^6*w^8*z^52-23*X^6*w^4*z^56+18*X^4*w^14*z^48+10*X
^4*w^10*z^52-4*X^4*w^6*z^56+18564*X^9*w^8*z^48+660*X^7*w^14*z^44-1232*X^7*w^10*
z^48+818*X^7*w^6*z^52-22*X^7*w^2*z^56-25*X^5*w^16*z^44+17*X^5*w^12*z^48+16*X^5*
w^8*z^52+2*X^5*w^4*z^56+2*X^3*w^14*z^48+2*X^3*w^6*z^56+2002*X^8*w^12*z^44-5642*
X^8*w^8*z^48+364*X^8*w^4*z^52-224*X^6*w^14*z^44+228*X^6*w^10*z^48-126*X^6*w^6*z
^52+30*X^6*w^2*z^56-8*X^4*w^16*z^44+19*X^4*w^12*z^48-4*X^4*w^8*z^52-2*X^4*w^4*z
^56-2*X^2*w^14*z^48+99*X^7*w^12*z^44+914*X^7*w^8*z^48-198*X^7*w^4*z^52+24*X^5*w
^14*z^44+6*X^5*w^10*z^48+6*X^5*w^6*z^52-6*X^5*w^2*z^56+X^3*w^16*z^44+3*X^3*w^12
*z^48+6*X^3*w^8*z^52-2002*X^8*w^10*z^44+2002*X^8*w^6*z^48-70*X^6*w^16*z^40-58*X
^6*w^12*z^44-122*X^6*w^8*z^48-3*X^6*w^4*z^52-X^6*z^56-14*X^4*w^14*z^44-6*X^4*w^
10*z^48+10*X^4*w^6*z^52-2*X^4*w^2*z^56-X^2*w^12*z^48-X^2*w^8*z^52-31824*X^9*w^8
*z^44-924*X^7*w^14*z^40+1606*X^7*w^10*z^44-1232*X^7*w^6*z^48+110*X^7*w^2*z^52+
10*X^5*w^16*z^40+46*X^5*w^12*z^44-6*X^5*w^8*z^48-18*X^5*w^4*z^52+2*X^5*z^56-10*
X^3*w^14*z^44+8*X^3*w^10*z^48-3003*X^8*w^12*z^40+6864*X^8*w^8*z^44-1001*X^8*w^4
*z^48+196*X^6*w^14*z^40-334*X^6*w^10*z^44+228*X^6*w^6*z^48-98*X^6*w^2*z^52+7*X^
4*w^16*z^40-14*X^4*w^12*z^44-2*X^4*w^8*z^48-X^4*w^4*z^52-X^4*z^56+2*X^2*w^14*z^
44-2*X^2*w^10*z^48-396*X^7*w^12*z^40-736*X^7*w^8*z^44+165*X^7*w^4*z^48-32*X^5*w
^14*z^40+30*X^5*w^10*z^44+6*X^5*w^6*z^48+4*X^5*w^2*z^52-7*X^3*w^12*z^44+X^3*w^8
*z^48+4*X^3*w^4*z^52+858*X^8*w^10*z^40-2002*X^8*w^6*z^44+56*X^6*w^16*z^36+98*X^
6*w^12*z^40+54*X^6*w^8*z^44+44*X^6*w^4*z^48+8*X^6*z^52-14*X^4*w^14*z^40+28*X^4*
w^10*z^44-6*X^4*w^6*z^48-2*X^4*w^2*z^52-5*X^2*w^12*z^44+43758*X^9*w^8*z^40+924*
X^7*w^14*z^36-1848*X^7*w^10*z^40+1606*X^7*w^6*z^44-330*X^7*w^2*z^48+3*X^5*w^16*
z^36-84*X^5*w^12*z^40-38*X^5*w^8*z^44+17*X^5*w^4*z^48-11*X^5*z^52+6*X^3*w^14*z^
40-4*X^3*w^10*z^44+8*X^3*w^6*z^48+2*X^3*w^2*z^52+X*w^12*z^44+3432*X^8*w^12*z^36
-6578*X^8*w^8*z^40+2002*X^8*w^4*z^44-140*X^6*w^14*z^36+478*X^6*w^10*z^40-334*X^
6*w^6*z^44+184*X^6*w^2*z^48-2*X^4*w^16*z^36+X^4*w^12*z^40+50*X^4*w^8*z^44+19*X^
4*w^4*z^48+2*X^4*z^52-4*X^2*w^10*z^44-2*X^2*w^6*z^48+396*X^7*w^12*z^36+110*X^7*
w^8*z^40+99*X^7*w^4*z^44-4*X^5*w^14*z^36-44*X^5*w^10*z^40+30*X^5*w^6*z^44-4*X^5
*w^2*z^48-6*X^3*w^12*z^40-2*X^3*w^8*z^44+3*X^3*w^4*z^48+X^3*z^52+858*X^8*w^10*z
^36+858*X^8*w^6*z^40-28*X^6*w^16*z^32-214*X^6*w^12*z^36-116*X^6*w^8*z^40-58*X^6
*w^4*z^44-28*X^6*z^48+18*X^4*w^14*z^36-24*X^4*w^10*z^40+28*X^4*w^6*z^44+18*X^4*
w^2*z^48+X^2*w^12*z^40-4*X^2*w^8*z^44-X^2*w^4*z^48-48620*X^9*w^8*z^36-660*X^7*w
^14*z^32+1848*X^7*w^10*z^36-1848*X^7*w^6*z^40+660*X^7*w^2*z^44-4*X^5*w^16*z^32+
20*X^5*w^12*z^36+134*X^5*w^8*z^40+46*X^5*w^4*z^44+24*X^5*z^48-6*X^3*w^10*z^40-4
*X^3*w^6*z^44+2*X^3*w^2*z^48-3003*X^8*w^12*z^32+4862*X^8*w^8*z^36-3003*X^8*w^4*
z^40+88*X^6*w^14*z^32-614*X^6*w^10*z^36+478*X^6*w^6*z^40-224*X^6*w^2*z^44-12*X^
4*w^12*z^36-48*X^4*w^8*z^40-14*X^4*w^4*z^44+2*X^4*z^48-6*X^2*w^10*z^40-4*X^2*w^
6*z^44-2*X^2*w^2*z^48-99*X^7*w^12*z^32+462*X^7*w^8*z^36-396*X^7*w^4*z^40+36*X^5
*w^14*z^32-16*X^5*w^10*z^36-44*X^5*w^6*z^40+24*X^5*w^2*z^44-2*X^3*w^12*z^36-10*
X^3*w^8*z^40-7*X^3*w^4*z^44-2*X^3*z^48+2*X*w^10*z^40-2002*X^8*w^10*z^32+858*X^8
*w^6*z^36+8*X^6*w^16*z^28+293*X^6*w^12*z^32+260*X^6*w^8*z^36+98*X^6*w^4*z^40+56
*X^6*z^44-2*X^4*w^14*z^32-32*X^4*w^10*z^36-24*X^4*w^6*z^40-14*X^4*w^2*z^44+5*X^
2*w^12*z^36-10*X^2*w^8*z^40-5*X^2*w^4*z^44+43758*X^9*w^8*z^32+330*X^7*w^14*z^28
-1606*X^7*w^10*z^32+1848*X^7*w^6*z^36-924*X^7*w^2*z^40+X^5*w^16*z^28+43*X^5*w^
12*z^32-146*X^5*w^8*z^36-84*X^5*w^4*z^40-25*X^5*z^44-20*X^3*w^10*z^36-6*X^3*w^6
*z^40-10*X^3*w^2*z^44+X*w^8*z^40+X*w^4*z^44+2002*X^8*w^12*z^28-2574*X^8*w^8*z^
32+3432*X^8*w^4*z^36-44*X^6*w^14*z^28+570*X^6*w^10*z^32-614*X^6*w^6*z^36+196*X^
6*w^2*z^40-3*X^4*w^12*z^32-28*X^4*w^8*z^36+X^4*w^4*z^40-8*X^4*z^44+6*X^2*w^10*z
^36-6*X^2*w^6*z^40+2*X^2*w^2*z^44-165*X^7*w^12*z^28-616*X^7*w^8*z^32+396*X^7*w^
4*z^36-26*X^5*w^14*z^28+42*X^5*w^10*z^32-16*X^5*w^6*z^36-32*X^5*w^2*z^40+11*X^3
*w^12*z^32-2*X^3*w^8*z^36-6*X^3*w^4*z^40+X^3*z^44+2*X*w^10*z^36+2*X*w^6*z^40+
2002*X^8*w^10*z^28-2002*X^8*w^6*z^32-X^6*w^16*z^24-229*X^6*w^12*z^28-336*X^6*w^
8*z^32-214*X^6*w^4*z^36-70*X^6*z^40-2*X^4*w^14*z^28+32*X^4*w^10*z^32-32*X^4*w^6
*z^36-14*X^4*w^2*z^40+2*X^2*w^8*z^36+X^2*w^4*z^40-31824*X^9*w^8*z^28-110*X^7*w^
14*z^24+1232*X^7*w^10*z^28-1606*X^7*w^6*z^32+924*X^7*w^2*z^36-37*X^5*w^12*z^28-\
48*X^5*w^8*z^32+20*X^5*w^4*z^36+10*X^5*z^40+18*X^3*w^10*z^32-20*X^3*w^6*z^36+6*
X^3*w^2*z^40+4*X*w^8*z^36-1001*X^8*w^12*z^24+728*X^8*w^8*z^28-3003*X^8*w^4*z^32
+14*X^6*w^14*z^24-348*X^6*w^10*z^28+570*X^6*w^6*z^32-140*X^6*w^2*z^36+27*X^4*w^
12*z^28+26*X^4*w^8*z^32-12*X^4*w^4*z^36+7*X^4*z^40+2*X^2*w^10*z^32+6*X^2*w^6*z^
36-w^8*z^36+198*X^7*w^12*z^24+516*X^7*w^8*z^28-99*X^7*w^4*z^32+6*X^5*w^14*z^24-\
16*X^5*w^10*z^28+42*X^5*w^6*z^32-4*X^5*w^2*z^36-3*X^3*w^12*z^28-8*X^3*w^8*z^32-\
2*X^3*w^4*z^36+2*X*w^6*z^36-1274*X^8*w^10*z^24+2002*X^8*w^6*z^28+101*X^6*w^12*z
^24+252*X^6*w^8*z^28+293*X^6*w^4*z^32+56*X^6*z^36-2*X^4*w^10*z^28+32*X^4*w^6*z^
32+18*X^4*w^2*z^36+10*X^2*w^8*z^32+5*X^2*w^4*z^36+18564*X^9*w^8*z^24+22*X^7*w^
14*z^20-818*X^7*w^10*z^24+1232*X^7*w^6*z^28-660*X^7*w^2*z^32+18*X^5*w^12*z^24+
214*X^5*w^8*z^28+43*X^5*w^4*z^32+3*X^5*z^36+2*X^3*w^10*z^28+18*X^3*w^6*z^32+X*w
^8*z^32+364*X^8*w^12*z^20+182*X^8*w^8*z^24+2002*X^8*w^4*z^28-2*X^6*w^14*z^20+
154*X^6*w^10*z^24-348*X^6*w^6*z^28+88*X^6*w^2*z^32-15*X^4*w^12*z^24+12*X^4*w^8*
z^28-3*X^4*w^4*z^32-2*X^4*z^36+4*X^2*w^10*z^28+2*X^2*w^6*z^32-102*X^7*w^12*z^20
-432*X^7*w^8*z^24-165*X^7*w^4*z^28+8*X^5*w^10*z^24-16*X^5*w^6*z^28+36*X^5*w^2*z
^32+12*X^3*w^8*z^28+11*X^3*w^4*z^32+546*X^8*w^10*z^20-1274*X^8*w^6*z^24-23*X^6*
w^12*z^20-76*X^6*w^8*z^24-229*X^6*w^4*z^28-28*X^6*z^32+4*X^4*w^10*z^24-2*X^4*w^
6*z^28-2*X^4*w^2*z^32+2*X^2*w^8*z^28-8568*X^9*w^8*z^20-2*X^7*w^14*z^16+434*X^7*
w^10*z^20-818*X^7*w^6*z^24+330*X^7*w^2*z^28-10*X^5*w^12*z^20-156*X^5*w^8*z^24-\
37*X^5*w^4*z^28-4*X^5*z^32+2*X^3*w^10*z^24+2*X^3*w^6*z^28-91*X^8*w^12*z^16-336*
X^8*w^8*z^20-1001*X^8*w^4*z^24-74*X^6*w^10*z^20+154*X^6*w^6*z^24-44*X^6*w^2*z^
28+X^4*w^12*z^20+27*X^4*w^4*z^28+4*X^2*w^6*z^28+27*X^7*w^12*z^16+374*X^7*w^8*z^
20+198*X^7*w^4*z^24-6*X^5*w^10*z^20+8*X^5*w^6*z^24-26*X^5*w^2*z^28+5*X^3*w^8*z^
24-3*X^3*w^4*z^28-154*X^8*w^10*z^16+546*X^8*w^6*z^20+2*X^6*w^12*z^16-32*X^6*w^8
*z^20+101*X^6*w^4*z^24+8*X^6*z^28-2*X^4*w^10*z^20+4*X^4*w^6*z^24-2*X^4*w^2*z^28
+3060*X^9*w^8*z^16-164*X^7*w^10*z^16+434*X^7*w^6*z^20-110*X^7*w^2*z^24+3*X^5*w^
12*z^16+34*X^5*w^8*z^20+18*X^5*w^4*z^24+X^5*z^28-2*X^3*w^10*z^20+2*X^3*w^6*z^24
+14*X^8*w^12*z^12+194*X^8*w^8*z^16+364*X^8*w^4*z^20+40*X^6*w^10*z^16-74*X^6*w^6
*z^20+14*X^6*w^2*z^24-8*X^4*w^8*z^20-15*X^4*w^4*z^24-3*X^7*w^12*z^12-266*X^7*w^
8*z^16-102*X^7*w^4*z^20-2*X^5*w^10*z^16-6*X^5*w^6*z^20+6*X^5*w^2*z^24-2*X^3*w^8
*z^20+26*X^8*w^10*z^12-154*X^8*w^6*z^16+30*X^6*w^8*z^16-23*X^6*w^4*z^20-X^6*z^
24-2*X^4*w^10*z^16-2*X^4*w^6*z^20+X^2*w^8*z^20-816*X^9*w^8*z^12+38*X^7*w^10*z^
12-164*X^7*w^6*z^16+22*X^7*w^2*z^20+10*X^5*w^8*z^16-10*X^5*w^4*z^20-2*X^3*w^6*z
^20-X^8*w^12*z^8-64*X^8*w^8*z^12-91*X^8*w^4*z^16-14*X^6*w^10*z^12+40*X^6*w^6*z^
16-2*X^6*w^2*z^20+4*X^4*w^8*z^16+X^4*w^4*z^20+136*X^7*w^8*z^12+27*X^7*w^4*z^16+
2*X^5*w^10*z^12-2*X^5*w^6*z^16-X^3*w^8*z^16-2*X^8*w^10*z^8+26*X^8*w^6*z^12-2*X^
6*w^8*z^12+2*X^6*w^4*z^16-2*X^4*w^6*z^16+153*X^9*w^8*z^8-4*X^7*w^10*z^8+38*X^7*
w^6*z^12-2*X^7*w^2*z^16-8*X^5*w^8*z^12+3*X^5*w^4*z^16+12*X^8*w^8*z^8+14*X^8*w^4
*z^12+2*X^6*w^10*z^8-14*X^6*w^6*z^12-2*X^4*w^8*z^12-47*X^7*w^8*z^8-3*X^7*w^4*z^
12+2*X^5*w^6*z^12-2*X^8*w^6*z^8-4*X^6*w^8*z^8-18*X^9*w^8*z^4-4*X^7*w^6*z^8+2*X^
5*w^8*z^8-X^8*w^8*z^4-X^8*w^4*z^8+2*X^6*w^6*z^8+10*X^7*w^8*z^4+X^6*w^8*z^4+X^9*
w^8-X^7*w^8)/(X^3*w^2*z^24-X^2*w^4*z^20-6*X^3*w^2*z^20-X^2*w^2*z^20+3*X^2*w^4*z
^16-X^2*z^20+15*X^3*w^2*z^16+X*w^4*z^16+2*X^2*w^2*z^16+X*w^2*z^16-3*X^2*w^4*z^
12+3*X^2*z^16-20*X^3*w^2*z^12-X*w^4*z^12+X*z^16+X^2*w^4*z^8-3*X^2*z^12-w^2*z^12
+15*X^3*w^2*z^8-X*z^12-2*X^2*w^2*z^8-X*w^2*z^8+X^2*z^8-6*X^3*w^2*z^4+X^2*w^2*z^
4+X^3*w^2)]:
L[m]:
end:

#GFm(m,X,x1,x2,y): the generating function, in X, for P_{m,n}(x1,x2,y)
#Try:
#GFm(2,X,x1,x2,y);
GFm:=proc(m,X,x1,x2,y) local i,j,gu,mu,mu1:
gu:=Mm(m,x1,x2,y):
mu:=[]:

for i from 1 to nops(gu) do
 mu1:=[]:
 for j from 1 to nops(gu[i]) do
  if i=j then
   mu1:=[op(mu1),1-X*gu[i][j]]:
  else
     mu1:=[op(mu1),-X*gu[i][j]]:
  fi:
 od:
mu:=[op(mu),mu1]:
od:
mu:=matrix(mu):    
mu:=inverse(mu):

normal(add(mu[i,i],i=1..2^m)):
end:


####From C-finite (adapted)

#CartProd2(L1,L2): The Cartesian product of lists L1 and L2,
#for example, try:
#CartProd2([2,5],[4,9]);
CartProd2:=proc(L1,L2) local i,j:
[seq(seq(L1[i]*L2[j],i=1..nops(L1)),j=1..nops(L2))]:
end:


#CartProd(L): The Cartesian product of the lists in the list
#of lists L
#for example, try:
#CartProd([[2,5],[4,9],[4,7]]);
CartProd:=proc(L) local L1,L2,i:
if not type(L,list) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

if {seq(type(L[i],list),i=1..nops(L))}<>{true} then
 print(`Bad input`):
 RETURN(FAIL):
fi:

if nops(L)=1 then
RETURN(L[1]):
fi:


L1:=CartProd([op(1..nops(L)-1,L)]):
L2:=L[nops(L)]:
CartProd2(L1,L2):

end:

#Krovim(L,i,eps): Given a list L and an index between
#1 and nops(L), finds the set of indices (including i)
#such that abs(L[j]-L[i])<eps. For example, try:
#Krovim([1.0001,1.01,1.00011],1,1/1000);
Krovim:=proc(L,i,eps) local gu,j,a:
gu:={}:
a:=L[i]:

for j from 1 to nops(L) do
 if abs(L[j]-a)<eps then
  gu:=gu union {j}:
 fi:
od:
gu:
end:

#ProdProfile(L,eps): Given a list of complex numbers finds its
#(approximate profile). For example, try:
#ProdProfile([1,2,3,6],10^(-4));
ProdProfile:=proc(L,eps) local i,j,lu,mu,S,gu:

gu:=[seq(seq(L[i]/L[j],j=1..nops(L)),i=1..nops(L))]:
mu:=[]:
S:={seq(i,i=1..nops(gu))}:

while S<>{} do
 lu:=Krovim(gu,S[1],eps):
 mu:=[op(mu),nops(lu)]:
 S:=S minus lu:
od:
sort(mu):
end:

#ProdIndicatorG(L): The profile of multiplicity of
#the convert(L,`*`)^2 ratios of a list of 
#be the Cartesian product of lists of size given by L
#For example, try:
#ProdIndicatorG([2,2]);
ProdIndicatorG:=proc(L) local a,i,j,gu,x,lu,mu,M:
option remember:
M:=[seq([seq(a[i][j],j=1..L[i])],i=1..nops(L))]:
gu:=CartProd(M):

lu:=[seq(seq(gu[i]/gu[j],j=1..nops(gu)),i=1..nops(gu))]:
lu:=add(x[lu[i]],i=1..nops(lu)):

mu:=[]:

for i from 1 to nops(lu) do
 if type(op(1,op(i,lu)),integer) then
   mu:=[op(mu),op(1,op(i,lu))]:
 else:
   mu:=[op(mu),1]:
 fi:
od:

sort(mu):

end:


#IsProdG(f,t,L,eps): Given a polynomial  t of t
#decides whether the C-finite sequence of its 
#roots are Cartesian products of sets
#of cardinalities given by the list L. eps should be small
#For example, try:
#IsProdG((1-t)*(1-2*t)*(1-3*t)*(1-6*t),t,[2,2],10^(-10));
IsProdG:=proc(lu,t,L,eps) local gu,ku1,ku2:


if degree(lu,t)<>convert(L,`*`) then
 print(`The degree of  of`, lu , `should be`, convert(L,`*`)):
 RETURN(FAIL):
fi:

gu:=[fsolve(lu,t)]:


ku1:=ProdProfile(gu,eps):
ku2:=ProdIndicatorG(L):

if ku1=ku2 then
 RETURN(true):
else

 RETURN(false,[ku1,ku2]):
fi:

end:

#####################End from Cfinite


#MyMax1(f,x,x0,x1,K): numerically finds an interval
#containing the max
#an expression f in x, between x0 and x1 with
#resolution 10^(-K), try:
#MayMax1(1-x^2,x,-.1,.1,1);
MyMax1:=proc(f,x,x0,x1,K) local i:

for i from 0 to (x1-x0)*10^K do
if subs(x=x0+i*10.^(-K),f)<subs(x=x0+(i+1)*10.^(-K),f)
and subs(x=x0+(i+1)*10.^(-K),f)>subs(x=x0+(i+2)*10.^(-K),f) then
 RETURN([x0+i*10.^(-K), x0+(i+2)*10.^(-K)]):
fi:
od:
FAIL:
end:


#MyMax(f,x,x0,x1,K): numerically finds an interval
#containing the max
#an expression f in x, between x0 and x1 with
#resolution 10^(-K), try:
#MayMax(1-x^2,x,-.1,.1,1);
MyMax:=proc(f,x,x0,x1,K) local i,gu:
Digits:=20*K:
gu:=MyMax1(f,x,x0,x1,1):

if gu=FAIL then
 RETURN(FAIL):
fi:

for i from 2 to K do
gu:=MyMax1(f,x,gu[1],gu[2],i):
if gu=FAIL then
 RETURN(FAIL):
fi:
od:
evalf((gu[1]+gu[2])/2,K+1):
end:







#OnsFun1(x): the Onsager solution, the integral part
OnsFun1:=proc(x) local t,gu:
gu:=1/2*(x^2+1/x^2)^2/(x^2-1/x^2):
evalf(1/2/Pi*Int(arccosh(gu-cos(t)), t=0..Pi)):
end:


#OnsFun(x): the Onsager solution
OnsFun:=proc(x) local t,gu:
gu:=1/2*(x^2+1/x^2)^2/(x^2-1/x^2):
evalf(1/2/Pi*Int(arccosh(gu-cos(t)), t=0..Pi)+1/2*log(x^2-1/x^2)):
:
end:

#OnsFunDer(x0,r): the r-th derivative of the Onsager solution OnsFun(x) w.r.t. to x, at x=x0, done numerically. Try:
#OnsFunDer(1.2,2);
OnsFunDer:=proc(x0,r) local t,gu,x:

if r=0 then
 RETURN(OnsFun(x0)):
fi:

gu:=1/2*(x^2+1/x^2)^2/(x^2-1/x^2):

gu:=arccosh(gu-cos(t));

gu:=normal(diff(gu,x$r)):

#+1/2*diff(log(x^2-1/x^2),x):

evalf(1/2/Pi*Int(subs(x=x0,gu), t=0..Pi)+1/2*subs(x=x0,diff(log(x^2-1/x^2),x$r)) ):

end:

#OnsFunm0(m0,x): the Onsager solution for an m0 by infinity strip
OnsFunm0:=proc(m0,x) local t,gu,k:
gu:=1/2*(x^2+1/x^2)^2/(x^2-1/x^2):
1/2*log(x^2-1/x^2)+
add(evalf(arccosh(gu-cos(Pi*(2*k+1)/m0)))/2/m0, k=0..m0-1):
end:



#UC(Pol,P): applies the umbra P^(i)->1/2/Pi*int(cos(t)^i,t=0..Pi):
UC:=proc(Pol,P) local i,gu:
gu:=coeff(Pol,P,0)/2:
for i from 1 to trunc(degree(Pol,P)/2) do
gu:=gu+coeff(Pol,P,2*i)*binomial(2*i-1,i)/2^(2*i):
od:
gu:
end:



#GuessRecF(f,v,P,C): inputs a function f of v and a parameter P, tries to
#find a recurrence for the  coefficient of complexity<=C
#Try:
#GuessRecF( arccosh((1+2*v^2+v^4-2*v*P+2*v^3*P)/(2*v-2*v^3))+log(4*v/(1-v^2)),v,P,50);
GuessRecF:=proc(f,v,P,C) local gu,i,gu1,gu2,ope1,ope2,Deg1,Deg2,Ord1,Ord2:
gu1:=series(subs(P=3/2,f),v=0,C+10):
gu1:=[seq(coeff(gu1,v,i),i=1..C+5)]:
ope1:=Findrec(gu1,n,N,C):

gu2:=series(subs(P=15/7,f),v=0,C+10):
gu2:=[seq(coeff(gu2,v,i),i=1..C+5)]:
ope2:=Findrec(gu2,n,N,C):

Deg1:=degree(numer(ope1),n):
Deg2:=degree(numer(ope2),n):

if Deg1<>Deg2 then
print(Deg1,Deg2):
 RETURN(FAIL):
fi:

Ord1:=degree(ope1,N):
Ord2:=degree(ope2,N):

if Ord1<>Ord2 then
print(`Orders are:`, Ord1, Ord2):
 RETURN(FAIL):
fi:

print(`There is hope for a recurrence of order`, Ord1, `and degree `, Deg1):
print(`Please be patient`):
gu:=series(f,v=0,C+10):
gu:=[seq(coeff(gu2,v,i),i=1..C+5)]:
findrec(gu,Deg1,Ord1,n,N):

end:



#ZmnPC(m,n,v): the expression in v=(x-1/x)/(x+1/x) of Pmn(m,n,x,x,1)/((x+1/x)/2)^N/2^N up to order K in v
#Try:
#ZmnPC(4,4,v,10);
ZmnPC:=proc(m,n,v) local gu,x,N,lu,mu:
option remember:
if m<3 or n<3 then
 print(`m and n must >=3 `):
fi:

gu:=PmnPC(m,n,x,1):
N:=m*n:
mu:=gu/((x+1/x)/2)^(2*N)/2^N:
mu:=simplify(subs(x=sqrt((1+v)/(1-v)),mu)):
mu:=expand(mu):


lu:=expand(normal(subs(v=(x-1/x)/(x+1/x),mu)*((x+1/x)/2)^(2*N)))*2^N:

if expand(lu-gu)<>0 then
 print(`something is wrong`):
  print(gu-lu,mu):
 FAIL:
else
 mu:
fi:

end:


#Zmn(m,n,v): the expression in v=(x-1/x)/(x+1/x) of Pmn(m,n,x,x,1)/((x+1/x)/2)^N/2^N 
#Try:
#Zmn(4,4,v);
Zmn:=proc(m,n,v) local gu,x,N,lu,mu:
option remember:

gu:=Pmn(m,n,x,x,1):
N:=m*n:
mu:=gu/((x+1/x)/2)^(2*N)/2^N:
mu:=simplify(subs(x=sqrt((1+v)/(1-v)),mu)):
mu:=expand(mu):


end:



Di:=proc(L) local i:
[seq(L[i+1]-L[i],i=1..nops(L)-1)]:
end:

#DiK(L,K):
DiK:=proc(L,K) local i,L1:
L1:=L:
for i from 1 to K do
L1:=Di(L1):
if L1=[] then
  RETURN(L1):
fi:
od:

L1:
end:




#GenPol1(x,dx,a): a generic polynomial of degree dx in x
#with coeffs. parametrized by a[i] followed by the set
#of coeffs.
GenPol1:=proc(x,dx,a) local i:
add(a[i]*x^i,i=0..dx),{seq(a[i],i=0..dx)}:
end:

#GP1d(S1,x,d): guesses a polynomial in x of degree d
#that fits the data in DS1
GP1d:=proc(S1,x,d)
local gu,eq,var,P,a,var1,i:
P:=GenPol1(x,d,a):
var:=P[2]:
P:=P[1]:

if nops(S1)<=d+2 then
 print(`Insufficient data`):
 RETURN(FAIL):
fi:

eq:={seq(subs(x=S1[i][1],P)-S1[i][2],i=1..d+2)}:
var1:=solve(eq,var):
if var1=NULL then
  RETURN(FAIL):
fi:
P:=subs(var1,P):


end:

#GP1(S1,x): guesses a polynomial from the data-set S1
GP1:=proc(S1,x) local d,gu:

for d from 1 to nops(S1)-3 do
 gu:=GP1d(S1,x,d):
  if gu<>FAIL then
    RETURN(factor(gu)):
  fi:
od:
FAIL:
end:

#CheckGP1(d): checks GP1(S1,x) for degree d 
CheckGP1:=proc(d)  local x,P,gu,t0:
P:=RandPol1(x,d):
t0:=time():
gu:=DS1(P,x):
evalb(normal(GP1(gu,x)-P)=0),time()-t0:
end:

#GPw1(m0,k,n,d1,K): inputs a positive integer m0, and a positive integer k, and
#a symbol n, outputs the polynomial in n, of degree d1, if it exists, and
#n0 such that for n>=n0 the coefficient of v^k in ZmnPC(m0,n,v) is given by it 
#it tries up to K terms
#Try:
#GP1w(3,4,n,1,20):
GPw1:=proc(m0,k,n,d1,K) local gu,n0,c,v:

gu:=[seq([n0,coeff(ZmnPC(m0,n0,v),v,k)],n0=3..d1+9)]:

for c from 1 to K while GP1d(gu,n,d1)=FAIL do
 gu:=[seq([n0,coeff(ZmnPC(m0,n0,v),v,k)],n0=c+3..c+d1+9)]:
od:

GP1d(gu,n,d1):

 
end:




#GPw(m0,k,n,d1,K): inputs a positive integer m0, and a positive integer k, and
#a symbol n, outputs the polynomial in n, of degree <=d1, if it exists, and
#n0 such that for n>=n0 the coefficient of v^k in ZmnPC(m0,n,v) is given by it 
#it tries up to K terms
#Try:
#GPw(3,4,n,2,20):
GPw:=proc(m0,k,n,d1,K) local gu,d11:

for d11 from 1 to d1 do
 gu:=GPw1(m0,k,n,d11,K):
 if gu<>FAIL then
  RETURN(gu):
 fi:
od:

FAIL:

end:


#OnsSeriesOld(K): the first K terms of the expansion, in v^2 (it is an even function)
#, of -psi/(k*T) (given by Eq. (5.8) in Colon  Thompson's book, where v is written nu)
#Try:
#OnsSeriesOld(10);
OnsSeriesOld:=proc(K) local f,v,t1,t2,i,gu:
f:=log(cosh(2*v)^2-sinh(2*v)*(cos(t1)+cos(t2))):
f:=taylor(f,v=0,2*K+1):

gu:=[seq(coeff(f,v,2*i),i=1..K)]:

[seq(1/2/Pi^2*int(int(gu[i],t1=0..Pi),t2=0..Pi),i=1..K)]:

end:

#UC2(Pol,P1,P2): applies the umbra P1^(i)*P2^j->1/(2*Pi)^2/2*int(int(cos(t1)^i*cos(t2)^j,t1=0..Pi),t2=0..Pi):
UC2:=proc(Pol,P1,P2) local i1,gu:


gu:=UC(coeff(Pol,P1,0),P2)/2:
for i1 from 1 to trunc(degree(Pol,P1)/2) do
gu:=gu+ binomial(2*i1-1,i1)/2^(2*i1)*UC(coeff(Pol,P1,2*i1),P2):
od:
gu:
end:

#OnsSeries(K): the first K terms of the expansion, in v^2 (it is an even function)
#, of -psi/(k*T) (given by Eq. (5.8) in Colon  Thompson's book, where v is written nu)
#Try:
#OnsSeries(10);
OnsSeries:=proc(K) local f,v,P1,P2,i,gu:
f:=log(cosh(2*v)^2-sinh(2*v)*(P1+P2)):
f:=taylor(f,v=0,2*K+1):

gu:=[seq(coeff(f,v,2*i),i=1..K)]:

[seq(2*UC2(gu[i],P1,P2),i=1..K)]:

end:


#OnsSeriesNor(K): the first K terms of the expansion, in v^2 (it is an even function)
#, of -psi/(k*T) times (2*i)!/2^i:
#Try:
#OnsSeriesNor(10);
OnsSeriesNor:=proc(K) local gu,i:

gu:=OnsSeries(K):
[seq(gu[i]*(2*i)!/2^i,i=1..nops(gu))]:
end:

#LarsSeriesNor(K): vSeriesFE(K) times n
#LarsSeriesNor(10);
LarsSeriesNor:=proc(K) local gu,i:

gu:=LarsSeries1(K):
[seq(gu[i]*i,i=1..nops(gu))]:
end:



#KZ(k0,m0,N): The Kauers-Zeilberger polynomials in N=m0*n. Inputs  positive integers k0 and m0>2, outputs the
#polynomial expression in N (valid for N/m0>=k+1) that equals to the coefficient of v^k0 in the series expansion of
#the partition function of the Ising model for an m0 by n rectangle, after it is divided
#by ((x+1/x)/2)^(2*n*m0)*2^(n*m0) and v=(x-1/x)/(x+1/x). Try:
#KZ(4,3,N);
KZ:=proc(k0,m0,N) local gu,d1,i,lu,v,n:
option remember:
gu:=[seq([i,coeff(Zmn(m0,i,v),v,k0)],i=k0+1..k0+5)]:

lu:=GP1(gu,n):

if lu<>FAIL then
lu:=factor(subs(n=N/m0,lu)):
 RETURN(lu):
fi:

for d1 from 1 do

gu:=[seq([i,coeff(Zmn(m0,i,v),v,k0)],i=k0+1..k0+5+d1)]:

lu:=GP1(gu,n):

if lu<>FAIL then
lu:=factor(subs(n=N/m0,lu)):
 RETURN(lu):
fi:

od:


end:

#KZpc(k0,m0,N): pre-computed version of KZ(k0,m0,N)
#KZpc(4,3,N);
KZpc:=proc(k0,m0,N) local gu,d1,i,lu,v,n:
option remember:
gu:=[seq([i,coeff(ZmnPC(m0,i,v),v,k0)],i=k0+1..k0+5)]:

lu:=GP1(gu,n):

if lu<>FAIL then
lu:=factor(subs(n=N/m0,lu)):
 RETURN(lu):
fi:

for d1 from 1 do

gu:=[seq([i,coeff(ZmnPC(m0,i,v),v,k0)],i=k0+1..k0+5+d1)]:

lu:=GP1(gu,n):

if lu<>FAIL then
lu:=factor(subs(n=N/m0,lu)):
 RETURN(lu):
fi:

od:


end:

#Rm(m,v,X): GFm(m,X,x,x,1) with the transormation to get the generating function in the v-polynomials. Try:
#Rm(2,v,X):
Rm:=proc(m,v,X) local gu,x:
option remember:
gu:=GFmPC(m,X,x,1):

gu:=normal(subs(X=X*2^m/(x+1/x)^(2*m),gu)):

gu:=normal(subs(x=sqrt(x),gu)):
gu:=normal(subs(x=(1+v)/(1-v),gu)):
factor(gu):
end:



#KZf(k0,m0,N): A fast version of KZ(k0,m0,n) done directly (no guessing). Try:
#KZf(4,3,N);
KZf:=proc(k0,m0,N) local gu,v,X,gu1,gu2,mu2,i,w,n:
option remember:

gu:=Rm(m0,v,X):

gu:=normal(coeff(taylor(gu,v=0,k0+1),v,k0)):

gu1:=numer(gu):

gu2:=denom(gu):



#mu1:=quo(gu1,gu2,X):
mu2:=rem(gu1,gu2,X):

gu:=mu2/gu2:

gu:=convert(gu,parfrac);

gu:=normal(subs(X=(w-1)/w,gu)):

gu:=factor(expand(add(coeff(gu,w,i)*binomial(n+i-1,i-1),i=1..degree(gu,w)))):

gu:=factor(subs(n=N/m0,gu)):
end:


#KZfS(K0,m0,N,v): the first K0 terms in the v-expansion
KZfS:=proc(K0,m0,N,v) local k0:
1+add(KZf(k0,m0,N)*v^k0,k0=1..K0):
end:

#KZS(K0,m0,N,v): the first K0 terms in the v-expansion by guessing
KZS:=proc(K0,m0,N,v) local k0:
1+add(KZ(k0,m0,N)*v^k0,k0=3..K0):
end:

#vSeries1(m0,K0): the first K0 terms in the v-series for an m0 by infinity Ising function. Try
#vSeries1(3,10);
vSeries1:=proc(m0,K0) local v,n,i,mu,x,lu:



#mu:=expand((x+1/x)^2):
#mu:=expand(subs(x=sqrt(x),mu)):
#mu:=normal(subs(x=(1-v)/(1+v),mu)):


lu:=taylor(log(KZfS(K0,m0,n,v) ),v=0,K0+1):

[seq(limit(coeff(lu,v,i)/n,n=infinity), i=1..K0)]:
end:



#vSeriesST(m0,K0): the first K0 terms in the v-series for an m0 by infinity Ising function. Try
#vSeriesST(3,10);
vSeriesST:=proc(m0,K0) local v,N,i,mu,x,lu,lu1,mu1:
mu:=expand((x+1/x)^2):
mu:=expand(subs(x=sqrt(x),mu)):
mu:=normal(subs(x=(1+v)/(1-v),mu)):
mu:=taylor(log(mu),v=0,K0+2):
mu1:=[seq(coeff(mu,v,i),i=1..K0)]:

lu1:=[0,0,seq(coeff(KZ(i,m0,N),N,1),i=3..K0)]:
lu1+mu1:

end:

#vSeriesSTf(m0,K0): the first K0 terms in the v-series for an m0 by infinity Ising function. Try
#vSeriesSTf(3,10);
vSeriesSTf:=proc(m0,K0) local v,N,i,mu,x,lu,lu1,mu1:
mu:=expand((x+1/x)^2):
mu:=expand(subs(x=sqrt(x),mu)):
mu:=normal(subs(x=(1+v)/(1-v),mu)):
mu:=taylor(log(mu),v=0,K0+2):
mu1:=[seq(coeff(mu,v,i),i=1..K0)]:
lu1:=[0,0,seq(coeff(KZf(i,m0,N),N,1),i=3..K0)]:
lu1+mu1:

end:




#vSeriesFEviaLog(K): the first K terms of the Onsager power series but by naively do
#log(ZnnPC(n,v))/n^2 for the smallest n>=2*K+2 that is a multiple of 4. Try
#vSeriesFEviaLog(8);
vSeriesFEviaLog:=proc(K) local n,gu:
for n from 4 by 4 to 48 while n<2*K+2 do od:
if n=52 then
  RETURN(FAIL):
fi:

gu:=taylor(log(ZnnPC(n,v))/n^2,v=0,2*K+2):

[seq(coeff(gu,v,2*i)+1/i,i=1..K)]:

end:



#vSeriesFE(K): the first K terms of the Onsager power series 
#vSeriesFE(20);
vSeriesFE:=proc(K) local gu,v,i,P1,P2:
gu:=log(1+2*((v^2-1)*v/(1+v^2)^2)*(P1+P2))  -2*log(1-v^2)+2*log(1+v^2):
gu:=series(gu,v=0,2*K+4):

gu:=expand([seq(coeff(gu,v,i),i=1..2*K)]):
[seq(2*UC2(gu[2*i],P1,P2),i=1..K)]:

end:

#vSeriesU(K): the first K terms of the Onsager power series for the specfic Head J(d/d(nu) )(psi/(k*T)). Try
#vSeriesU(20);
vSeriesU:=proc(K) local gu:
gu:=vSeriesFE(K):
gu:=[seq(gu[i]*2*i,i=1..nops(gu))]:
[gu[1],seq(gu[i+1]-gu[i],i=1..K-1)]:
end:


#LarsSeries1(K): the first K terms of the Onsager power series plus log(1-v)
LarsSeries1:=proc(K) local gu,v,i,P1,P2:
gu:=2*log(1+v^2)+log(1+2*v*(v^2-1)/(1+v^2)^2*(P1+P2)):
gu:=series(gu,v=0,2*K+4):

gu:=expand([seq(coeff(gu,v,i),i=1..2*K)]):
[seq(2*UC2(gu[2*i],P1,P2),i=1..K)]:
end:


#LarsRec(n,N): the Recurrence equation satisfied by the coefficients of the v-series
#try:
#LarsRec(n,N);
LarsRec:=proc(n,N):
findrec(LarsSeries1(50),4,6,n,N):
end:


#LarsFun(v,K): The first K terms of the Lars series. Try
#LarsFun(v,20);
LarsFun:=proc(v,K) local i,lu:
lu:=LarsSeries1(K):
log(2)+add(lu[i]*v^(2*i),i=1..K):
end:





#PmnBK(m,n,x): Pmn(m,n,x,x,1) according to Bruria Kaufman's formula for numeric x. Try:
#PmnBK(3,4,2);
PmnBK:=proc(m,n,x) local ga,l,j:

for l from 1 to 2*m do
 ga[l]:=evalf(arccosh((x^2+1/x^2)^2/2/(x^2-1/x^2)-cos(l*Pi/m))):
od:

1/2*(x^2-1/x^2)^(m*n/2)*
(
mul(2*cosh(n*ga[2*j-1]/2),j=1..m)+
mul(2*sinh(n*ga[2*j-1]/2),j=1..m)+
mul(2*cosh(n*ga[2*j]/2),j=1..m)+
mul(2*sinh(n*ga[2*j]/2),j=1..m)
):

end:




#LarsDE(u,Du): the Differntial equation of the Onsager function in terms of u=tanh(J/(k*T))^2 alias ((1-x^2)/(1+x^2))^2
#try:
#LarsDE(u,Du);
LarsDE:=proc(u,Du) local n, N,gu,i,ope:
gu:=numer(LarsRec(n,N)):
ope:=rectodiff(gu,n,N,u,Du):
add(factor(coeff(ope,Du,i))*Du^i,i=0..degree(ope,Du)):
end:


#OnsFunm01(m0,v): the m0 by infinity strip Onsager function in terms of v=(x-1/x)/(x+1/x)
#Try:
#plot(OnsFunm01(4,v),x=1.2..3);
OnsFunm01:=proc(m0,v) local k:
evalf(
1/2*log(4*v/(1-v^2))+
1/2/m0*
add( 
arccosh( (1+v^2)^2/2/(1-v^2)/v- cos(Pi*(2*k+1)/m0)) ,k=0..m0-1)
):
end:





#LarsSeriesStExact(K,m0): the first K terms of the Onsager power series for an m0 by infinity strip
LarsSeriesStExact:=proc(K,m0) local gu,v,k,i,w,mu,top1,bot1,lu1,c:

gu:=1/2*log(4*v/(1-v^2))+1/2/m0*add( arccosh((1+v^2)^2/2/(1-v^2)/v- (w^(2*k+1)+1/w^(2*k+1))/2 ) ,k=0..m0-1):
gu:=expand(simplify(gu)):
gu:=series(gu,v=0,2*K+3):
gu:=expand([seq(coeff(gu,v,i),i=1..2*K)]):

gu:=expand(gu):

mu:=[]:

for i from 1 to nops(gu) do
lu1:=gu[i]:
top1:=numer(lu1):
bot1:=denom(lu1):
c:=lcoeff(bot1,w):
lu1:= round(coeff(evalf(evalc(subs(w=exp(Pi*I/m0),top1/bot1*c))),I,0)):
mu:=[op(mu),lu1/c]:
od:

[seq(mu[2*i],i=1..K)]:

end:



#EvalCF(L): evaluates the simple confinued fraction L
EvalCF:=proc(L)
if nops(L)=1 then
 L[1]:
else
L[1]+1/EvalCF([op(2..nops(L),L)]):
fi:

end:

#FlToRat(a): converts the floating point a to a rational number with parameter K (make it big). Try
#FlToRat(2.500000000000000000000000001);
FlToRat:=proc(a) local L,i,gu:
if a<0 then
RETURN(-FlToRat(-a)):
fi:

L:=convert(a,confrac):
if L[1]=0 then
 if nops(L)=1 then
    RETURN(FAIL):
 else
  gu:=L[2]:
 fi:
else
 gu:=L[1]:
fi:

for i from 1 to nops(L) while L[i]<1000*gu do od:

EvalCF([op(1..i-1,L)]):
end:

#LarsSeriesSt(K,m0): the first K terms of the Onsager power series for an m0 by infinity strip. Try
#LarsSeriesSt(2,10);
LarsSeriesSt:=proc(K,m0) local gu,v,k,i:
Digits:=50:
gu:=1/2*log(4*v/(1-v^2))+1/2/m0*add( arccosh((1+v^2)^2/2/(1-v^2)/v- evalf(cos(Pi*(2*k+1)/m0)) ) ,k=0..m0-1):
gu:=series(gu,v=0,2*K+3):
gu:=expand([seq(coeff(gu,v,i),i=1..2*K)]):

gu:=expand(gu):

[seq(FlToRat(gu[2*i]),i=1..K)]:

end:


#zSeriesFE(K): the first K terms of the nice series done knowing Onsager's answer
#zSeriesFE(20);
zSeriesFE:=proc(K) local gu,z,i,P1,P2:
gu:=log(1+z*(P1+P2)):
gu:=series(gu,z=0,2*K+4):

gu:=expand([seq(coeff(gu,z,i),i=1..2*K)]):
[seq(2*UC2(gu[2*i],P1,P2),i=1..K)]:

end:


#vSeriesFEpure(K):  Like vSeriesFE(K), but subtracting the trivial part log(1+v^2)-log(1-v^2). Try:
#vSeriesFEpure(20);
vSeriesFEpure:=proc(K) local gu,gu1:
gu:=vSeriesFE(K):

gu1:=[seq(coeff(taylor(log(1+v^2)-log(1-v^2),v=0,2*K+2),v,2*i),i=1..K)]:
gu-gu1:
end:


#vToz(K):  Like zSeriesFE(K), but using the v-series (that we found from scratch!)
vToz:=proc(K) local gu,f,mu,i,v,z:

gu:=vSeriesFEpure(K):
f:=add(gu[i]*v^(2*i),i=1..K):

mu:=solve(2*v*(1-v)*(1+v)/(1+v^2)^2=z,v):

f:=taylor(subs(v=mu,f), z=0,2*K+1):

[seq(coeff(f,z,2*i),i=1..K)]:

end:



#zTov(K):  Getting the v series from the z-series, the reverse of vToz(K)
zTov:=proc(K) local gu,f,i,v,z,gu1:

gu:=zSeriesFE(K):

f:=add(gu[i]*z^(2*i),i=1..K):

f:=subs(z=2*v*(1-v)*(1+v)/(1+v^2)^2,f):

gu:=[seq(coeff(taylor(f,v=0,2*K+1),v,2*i),i=1..K)]:

gu1:=[seq(coeff(taylor(log(1+v^2)-log(1-v^2),v=0,2*K+2),v,2*i),i=1..K)]:
gu:=gu+gu1:
gu:

end:


#Aitken(L): The conjectured limit of the sequence L
#from its three last elements, using Aitken's fromula
Aitken:=proc(L)
local n,a,b,c:
n:=nops(L):

if n<3 then
 RETURN(FAIL):
fi:

a:=L[n-2]: b:=L[n-1]: c:=L[n]:

c-(c-b)^2/(c-2*b+a):
end:

#AppxOnsFun(x,y,K): The first K values of the log of the largest eigenvalue of Mm(m,x,x,y) for m from 2 to K. 
#Followed by the Aitken limit. Try:
#AppxOnsFun(1.5,1,6);
AppxOnsFun:=proc(x,y,K)  local m,L:
L:=[seq(log(eigenvalues(evalf(Mm(m,x,x,y)))[1])/m,m=2..K)]:
[L,Aitken(L)]:
end:


#AppxOnsFunTable1(xS,xT,res,y0,K):The values of the Approxminate Onsager function AppxOnsFun(x,y,K)[2], 
#for y=y0 from x=xS to x=xT with resolution res. Try:
#AppxOnsFunTable1(1,2,0.1,1,6);
AppxOnsFunTable1:=proc(xS,xT,res,y0,K) local n,i:
n:=trunc((xT-xS)/res):
[seq([xS+res*i, AppxOnsFun(xS+res*i,y0,K)[2] ],i=1..n)]:
end:

#CheckAppxOnsFunTable1(xS,xT,res,K):Checks AppxOnsFunTable1(xS,xT,res,y0,K) with y0=1 against the known values given by Onsager
#by taking the ratios
#for y=y0 from x=xS to x=xT with resolution res. Try:
#CheckAppxOnsFunTable1(1,2,0.1,6);
CheckAppxOnsFunTable1:=proc(xS,xT,res,K) local n,i:
n:=trunc((xT-xS)/res):
[seq([xS+res*i, evalf(AppxOnsFun(xS+res*i,1.0,K)[2]/OnsFun(xS+res*i),10) ],i=1..n)]:
end:

#AppxOnsFunTable(xS,xT,yS,yT,res,K):The values of the Approxminate Onsager function AppxOnsFun(x,y,K)[2], 
#for y=y0 from x=xS to x=xT and from y=yS to y=yT, with resolution res. Try:
#AppxOnsFunTable(1,1.4,1,1.4,0.1,5);
AppxOnsFunTable:=proc(xS,xT,yS,yT,res,K) local n1,n2, i1,i2, gu:
 n1:=trunc((xT-xS)/res):
 n2:=trunc((yT-yS)/res):

gu:=[]:

for i2 from 1 to n2 do
 gu:=[op(gu), [seq([xS+res*i1, yS+res*i2, evalf(AppxOnsFun(xS+res*i1,yS+res*i2,K)[2],10) ],i1=1..n1)]   ]:
od:

gu:

end:

#PlotAppxOnsFunTable(xS,xT,yS,yT,res,K):Plots AppxOnsFunTable(xS,xT,yS,yT,res,K) (q.v.). Try:
#PlotAppxOnsFunTable(1,1.4,1,1.4,0.1,5);
PlotAppxOnsFunTable:=proc(xS,xT,yS,yT,res,K) local gu,i,j:

gu:=AppxOnsFunTable(xS,xT,yS,yT,res,K):


gu:={seq(seq(gu[i][j],j=1..nops(gu[i])),i=1..nops(gu))};

plots[pointplot3d](gu):
end:

#AppxOnsFunTable1Dx(xS,xT,res,y0,K,r):The values of the (numerical) r-th (forward, numerical) 
#derivative, w.r.t. x,  of the Approxminate Onsager function AppxOnsFun(x,y,K)[2], 
#for y=y0 from x=xS to x=xT with resolution res. Try:
#AppxOnsFunTable1Dx(1,2,0.1,1,6,1);
AppxOnsFunTable1Dx:=proc(xS,xT,res,y0,K,r) local i,gu:
option remember:
if r=0 then
 RETURN(AppxOnsFunTable1(xS,xT,res,y0,K)):
fi:


gu:=AppxOnsFunTable1Dx(xS,xT+res,res,y0,K,r-1):
[seq([gu[i][1],(gu[i+1][2]-gu[i][2])/res],i=1..nops(gu)-1)]:
end:

#AppxOnsFunTable1Dy(xS,xT,res,y0,K,r):The values of the (numerical) r-th (forward, numerical) 
#derivative, w.r.t y, of the Approxminate Onsager function AppxOnsFun(x,y,K)[2], 
#for y=y0 from x=xS to x=xT with resolution res. Try:
#AppxOnsFunTable1Dy(1,2,0.1,1,6,1);
AppxOnsFunTable1Dy:=proc(xS,xT,res,y0,K,r) local i,gu,gu1:
option remember:
if r=0 then
 RETURN(AppxOnsFunTable1(xS,xT,res,y0,K)):
fi:


gu:=AppxOnsFunTable1Dy(xS,xT,res,y0,K,r-1):
gu1:=AppxOnsFunTable1Dy(xS,xT,res,y0+res,K,r-1):
[seq([gu[i][1],(gu1[i+1][2]-gu[i][2])/res],i=1..nops(gu)-1)]:
end:

#IsingPolsPC(N): 
#the list of (pre-computed (by us!) Ising Polynomials as defined in Colin Thompson's book "Mathematical Statistical Mechanics"
#and elsewhere, following van der Waerden's combinatorial approach, the k-th term is the number of ways of placing 2*k-bond graphs
# in a very large torodial lattice with N vertices. Try:
#IsingPolsPC(N);

IsingPolsPC:=proc(N)
[ 0,N,2*N,(N*(9 + N))/2,2*N*(6 + N),(N*(7 + N)*(32 + N))/6,N*(130 + 21*N + N^2),N*(11766 + 1715*N + 102*N^2 + N^3)/24,
N*(5876 + 776*N + 49*N^2 + N^3)/3,(N*(980904 + 118830*N + 7415*N^2 + 210*N^3 + N^4))/120,
N*(423624 + 47666*N + 2855*N^2 + 94*N^3 + N^4)/12,N*(112852800 + 11919274*N + 678945*N^2 + 23725*N^3 + 375*N^4 + N^5)/720,
N*(42723120 + 4272044*N + 231260*N^2 + 8175*N^3 + 160*N^4 + N^5)/60,
N*(16620978240 + 1584498216*N + 81728374*N^2 + 2851695*N^3 + 62545*N^4+ 609*N^5 + N^6)/5040,
N*(5589930384 + 510961484*N + 25204804*N^2 + 857825*N^3 + 19891*N^4 +251*N^5 + N^6)/360,
N*(2990184306480 + 263323487916*N + 12469823436*N^2 + 411847009*N^3 +9754920*N^4 + 143794*N^5 + 924*N^6 + N^7)/40320
]:
end:



UCLT:=proc(f,P1,P2) local i,j,gu:
gu:=0:

for i from 0 to degree(f,P1) do
 for j from 0 to degree(f,P2) do
  if i mod 2=0 and j mod 2=0 then
   gu:=gu+coeff(coeff(f,P1,i),P2,j)*binomial(i,i/2)/2^i*binomial(j,j/2)/2^j:
  fi:
 od:
od:
gu:
end:

#OnsSeriesLT(K): the first K terms of the expansion, in v^2 (it is an even function) of the exponential of the Low-Temperature series
#see Steven Finch's book "Mathematical Constants", p. 393
#Try:
#OnsSeriesLT(10);
OnsSeriesLT:=proc(K) local f,y,P1,P2,i:
f:=(1/2)*log((1+y^2)^2-2*y*(1-y^2)*(P1+P2)):
f:=taylor(f,y=0,2*K+1):
f:=add(UCLT(coeff(f,y,i),P1,P2)*y^i,i=0..2*K):
f:=exp(f):
f:=taylor(f,y=0,2*K+1):
[seq(coeff(f,y,2*i),i=1..K)]:

end:






#OnsPolPC(n,v): the weight-enumerator of n by n {1,-1} matrices under periodic boundary condition in the Ising model 
#divided by (x+1/x)^(2*n^2) when expressed in term of v=(x-1/x)/(x+1/x), using the pre-computed values. Only valid
#for n a multipe of 4 between 4 and 48. Try
#OnsPolPC(4,v);
OnsPolPC:=proc(n,v) local gu,mu,N,lu,x:

if not (n mod 4=0 and n>=4 and n<=48) then
 print(`Not yet implemented`):
  RETURN(FAIL):
fi:


gu:=Ons(x)[n/4]:
N:=n^2:
mu:=gu/((x+1/x)/2)^(2*N)/2^N:
mu:=simplify(subs(x=sqrt((1+v)/(1-v)),mu)):
mu:=expand(mu):


lu:=expand(normal(subs(v=(x-1/x)/(x+1/x),mu)*((x+1/x)/2)^(2*N)))*2^N:

if expand(lu-gu)<>0 then
 print(`something is wrong`):
  print(gu-lu,mu):
 FAIL:
else
 mu:
fi:

end:

#OnsPolPCt(n,v,k): the weight-enumerator of n by n {1,-1} matrices under periodic boundary condition in the Ising model 
#divided by (x+1/x)^(2*n^2) when expressed in term of v=(x-1/x)/(x+1/x), truncated up to v^k, using the pre-computed values. Only valid
#for n a multipe of 4 between 4 and 48. Try
#OnsPolPCt(8,v,6);
OnsPolPCt:=proc(n,v,k) local gu,mu,N,x,i:
option remember:

if not (n mod 4=0 and n>=4 and n<=48) then
 print(`Not yet implemented`):
  RETURN(FAIL):
fi:


gu:=Ons(x)[n/4]:
N:=n^2:
mu:=gu/((x+1/x)/2)^(2*N)/2^N:
mu:=subs(x=sqrt((1+v)/(1-v)),mu):

mu:=taylor(mu,v=0,k+1):
add(coeff(mu,v,i)*v^i,i=0..k):

end:

#IsingPol(k,N): The Ising Polynomial for even-degree polygons with 2*k sides using pre-computed date
IsingPol:=proc(k,N) local mu,n1,v:

if trunc(k/2)+k+1>12 then
print(`k too big`):
 RETURN(FAIL):
fi:
 
mu:=[seq( [(4*n1)^2, coeff(OnsPolPCt(4*n1,v,2*k+1),v,2*k) ], n1=trunc(k/2)+1..trunc(k/2)+k+2)]:


GP1(mu,N):

end:





#ZmnD(m,n,v,d): the expression in v=(x-1/x)/(x+1/x) of diff(Pmn(m,n,x,x,y), y$d)/d!/((x+1/x)/2)^N/2^N 
#Try:
#ZmnD(4,4,v,d);
ZmnD:=proc(m,n,v,d) local gu,x,N,lu,mu,y:
option remember:

if d=0 then
 RETURN(Zmn(m,n,v)):
fi:

gu:=subs(y=1,diff(Pmn(m,n,x,x,y),y$d)/d!):
N:=m*n:
mu:=gu/((x+1/x)/2)^(2*N)/2^N:
mu:=simplify(subs(x=sqrt((1+v)/(1-v)),mu)):
mu:=expand(mu):


end:

#KZD(k0,m0,N,d): The generalized Kauers-Zeilberger polynomials in N=m0*n. Inputs  positive integers k0 and m0>2, outputs the
#polynomial expression in N (valid for N/m0>=k+1) that equals to the coefficient of v^k0 in the series expansion of
#the d-th derivative of the partition function of the Ising model for an m0 by n rectangle with magnetic field , evaluated at y=1 and divided by d!
#after it is divided
#by ((x+1/x)/2)^(2*n*m0)*2^(n*m0) and v=(x-1/x)/(x+1/x). Try:
#KZD(4,3,N,2);
KZD:=proc(k0,m0,N,d) local gu,d1,i,lu,v,n:
option remember:
if d=0 then
RETURN( KZ(k0,m0,N) ):
fi:
gu:=[seq([i,coeff(ZmnD(m0,i,v,d),v,k0)],i=k0+1..k0+5)]:

lu:=GP1(gu,n):

if lu<>FAIL then
lu:=factor(subs(n=N/m0,lu)):
 RETURN(lu):
fi:

for d1 from 1 do

gu:=[seq([i,coeff(ZmnD(m0,i,v,d),v,k0)],i=k0+1..k0+5+d1)]:

lu:=GP1(gu,n):

if lu<>FAIL then
lu:=factor(subs(n=N/m0,lu)):
 RETURN(lu):
fi:

od:


end:



#FunEqOld(m,n) verifies the functional equation for Pmn(m,n,x,x,1) as given in Eq. (2.8), chapter 6 of Thompson's book (p. 2.8). Try:
#FunEqOld(4,4);
FunEqOld:=proc(m,n) local x,gu,lu,mu,N:
N:=m*n:
gu:=Pmn(m,n,x,x,1):
lu:=normal(gu*2^(N-1)/(x^2-x^(-2))^N):
mu:=normal(subs(x=sqrt((x^2+1)/(x^2-1)),gu)):
normal(mu/lu):
end:


#MKa(x,y,n,v): Manuel Kauers' algorithm for fast computation of Mv. Try
#MKa(x,y,2,[1,1,1,1]);
MKa:=proc(x,y,n,v) local j,B1,k,Uj,Vj,vN,w,i:

if nops(v)<>2^n then
 RETURN(FAIL):
fi:

 vN:=[]:

for j from 0 to 2^n-1 do
 B1:=V2(j,n):
  Uj:=add(B1[k]*B1[k+1],k=1..n-1)+B1[n]*B1[1]:
  Vj:=convert(B1,`+`):
  vN:=expand([op(vN),v[j+1]*x^Uj*y^Vj]):
od:

for k from 0 to n-1 do
 for i from 0 to 2^(n-1)-1 do
  w[2*i]:=x*vN[2*i+1]+ 1/x*vN[2*i+2]:
  w[2*i+1]:=1/x*vN[2*i+1]+ x*vN[2*i+2]:
 od:
  vN:=expand([seq(w[2*i],i=0..2^(n-1)-1),seq(w[2*i+1],i=0..2^(n-1)-1)]):
od:
vN:
end:


#PmnMK(m,n,x,y): the trace of Mm(n,x,y)^m (alas Pmn(m,n,x,x,y)), using Manuel Kauers' fast algorithm. Try:
#PmnMK(2,2,x,y);
PmnMK:=proc(m,n,x,y) local t,i,v,k:

t:=0:

for i from 1 to 2^n do
 v:=[0$(i-1),1,0$(2^n-i)]:

 for  k from 0 to m-1 do
  v:=MKa(x,y,n,v):
 od:
t:=t+v[i]:
od:

t:

end:


#Amn(m,n,nu): Pmn(m,n,x,x,1) with x=exp(nu)
Amn:=proc(m,n,nu) local x,gu:
x:=exp(nu):
gu:=Pmn(m,n,x,x,1):
end:

#FunEqA(m,n,nu) verifies the functional equation for Pmn(m,n,x,x,1) where x=exp(nu)
#FunEq(4,4,0.3);
FunEqA:=proc(m,n,nu) local gu,N,nustar,L,R:
N:=m*n:

nustar:=arctanh(exp(-2*nu)):

L:=Amn(m,n,nu):
R:=Amn(m,n,nustar)*2^(1-N)*(2*sinh(2*nu))^N:
L/R:

end:








#OnsPolPCtz(n,z,k):  OnsPolPCt(n,z,k) in terms of z, where 2*v*(1-v)*(1+v)/(1+v^2)^2=z. Try:
#OnsPolPCtz(4,z,10);
OnsPolPCtz:=proc(n,z,k) local gu,v,i,mu,f:

if not (n mod 4=0 and n>=4 and n<=48) then
 print(`Not yet implemented`):
  RETURN(FAIL):
fi:

gu:=OnsPolPCt(n,v,2*k):

mu:=solve(2*v*(1-v)*(1+v)/(1+v^2)^2=z,v):

f:=taylor(subs(v=mu,gu), z=0,2*k+1):

[seq(coeff(f,z,2*i),i=1..k)]:

end:


#FunEq(n,K) verifies the functional equation for Pmn(n,n,x,x,1) as given in Eq. (2.8), chapter 6 of Thompson's book (p. 2.8)
#for n a multiple of 4 between 4 and 48
#FunEq(4,10);
FunEq:=proc(n,x) local gu,v,N,mu:

print(`Warning : Does not work!`):

if not (n mod 4=0 and n>=4 and n<=48) then
 print(`Not yet implemented`):
  RETURN(FAIL):
fi:
N:=n^2:

gu:=Ons(x)[n/4]:
mu:=normal(subs(x=sqrt((x^2+1)/(x^2-1)),gu)):

gu:=gu-2^(1-N)*(x^2-1/x^2)^N*mu:

expand(normal(gu)):

end:




#EVm(x,y,m): The eigenvalues of Mm(x,x,y,m). Try:
#EVm(2,1,3);
EVm:=proc(x,y,m) local gu:
gu:=evalf(LinearAlgebra[Eigenvalues](Matrix(Mm(m,x,x,y)))):
gu:=[seq(gu[i],i=1..2^m)]:
gu:=[seq(coeff(gu[i],I,0),i=1..nops(gu))]:
sort(gu):
end:

#CheckFunEqEv(K,m): 
#checks the functional equation for the eigenvaliesof Mm(exp(K),exp(K),1) given in Eq. (30) of the 1941 article
#by H.A. Kramers and G.H. Wannier, "Statistics of the Two-Dimensional Ferromagent", Part I, Physical Review v. 60 (1941), 252-262. Try
#CheckFunEqEv(2.,3);
CheckFunEqEv:=proc(K,m) local K1:
K1:=evalf(arcsinh(1/sinh(2*evalf(K)))/2):
EVm(exp(K),1,m)[2^m]/evalf(coeff(cosh(2*K),I,0)^m), EVm(exp(K1),1,m)[2^m]/evalf(coeff(cosh(2*K1),I,0))^m :
end:




#KWt(m,X,x): the Kramers-Wannier transform the equation whose root is  Lambda(K)/cosh(2*K). Try
#KWt(2,X,x);
KWt:=proc(m,X,x) local gu,mu:
gu:=Cm(m,X,x,x,1):
gu:=subs(x=sqrt(x),gu):
gu:=normal(subs(X=((x+1/x)/2)^m*X,gu)):

mu:=normal(subs(x=(x+1)/(x-1),gu)):

gu,mu:
end:

#KWtC(m,x0): Also like CheckFunEqEv(K,m), but does it exactly
#checks the functional equation for the eigenvaliesof Mm(exp(K),exp(K),1) given in Eq. (30) of the 1941 article
#by H.A. Kramers and G.H. Wannier, "Statistics of the Two-Dimensional Ferromagent", Part I, Physical Review v. 60 (1941), 252-262. 
#But via the largest root of the characteristic equation where x^2=x0
#Try
#KWtC(3,3);
KWtC:=proc(m,x0) local gu,x,X,mu:
gu:=Cm(m,X,x,x,1):
gu:=subs(x=sqrt(x),gu):
gu:=normal(subs(X=((x+1/x)/2)^m*X,gu)):

mu:=normal(subs(x=(x+1)/(x-1),gu)):

gu:=subs(x=x0,gu):

mu:=subs(x=x0,mu):

gu:=solve(gu,X):

mu:=solve(mu,X):



gu:=max(gu):

mu:=max(mu):

gu,mu:

end:



#ZnnPC(n,v): the expression in v=(x-1/x)/(x+1/x) of Pmn(n,n,x,x,1)/((x+1/x)/2)^N/2^N (where N=n^2)
#using the pre-computed Ising polynomial (with zero magnetic field) given by procedure Ons(x)
#Only works for n multiple of 4 between 4 and 48
#Try:
#ZnnPC(4,v);
ZnnPC:=proc(n,v) local gu,x,N,mu:
option remember:
if not (n mod 4=0 and n>=4 and n<=48) then
 print(`Not yet implemented`):
  RETURN(FAIL):
fi:
N:=n^2:

gu:=Ons(x)[n/4]:
N:=n^2:
mu:=gu/((x+1/x)/2)^(2*N)/2^N:
mu:=simplify(subs(x=sqrt((1+v)/(1-v)),mu)):
mu:=expand(mu):


end:


#ZnnPCt(n,v,K): The first K terms of the the expression in v=(x-1/x)/(x+1/x) of Pmn(n,n,x,x,1)/((x+1/x)/2)^N/2^N (where N=n^2)
#using the pre-computed Ising polynomial (with zero magnetic field) given by procedure Ons(x)
#Only works for n multiple of 4 between 4 and 48
#Try:
#ZnnPCt(4,v,10);
ZnnPCt:=proc(n,v,K) local gu,x,N,mu,i:
option remember:
if not (n mod 4=0 and n>=4 and n<=48) then
 print(`Not yet implemented`):
  RETURN(FAIL):
fi:
N:=n^2:

gu:=Ons(x)[n/4]:
N:=n^2:
mu:=gu/((x+1/x)/2)^(2*N)/2^N:
mu:=subs(x=sqrt(x),mu):
mu:=simplify(subs(x=(1+v)/(1-v),mu)):
mu:=taylor(mu,v=0,K+1):
add(coeff(mu,v,i)*v^i,i=0..K):
end:


#vSeriesFEd(N,K): the list of lists of first K terms of vSeriesFE for n by n Ising polynomials  for n from 4 to 4*N
#(N has to be at most 8). Try
#vSeriesFEd(2,20);
vSeriesFEd:=proc(N,K) local gu,mu,i,v,n:
if N>8 then
 print(`Not yet implemented`):
 RETURN(FAIL):
fi:

gu:=[]:
for n from 4 to 4*N by 4 do
 mu:=ZnnPCt(n,v,2*K+1):
 mu:=taylor(log(mu)/(n^2),v=0,2*K+1):
 mu:=[seq(coeff(mu,v,2*i),i=1..K)]:
 gu:=[op(gu),mu]:
print(`so far gu is`):
print(gu):

od:
gu:

end:



#vSeriesFEpureNew(K):  Like vSeriesFEpure(K), but NOT using Onsager, only doing by pure Symbolic computations. Try
#vSeriesFEpureNew(20);
vSeriesFEpureNew:=proc(K) local gu,n,i,v:

n:=2*K+4:
if n mod 4<>0 then
 n:=trunc(n/4)*4:
fi:
 gu:=ZnnPCt(n,v,2*K+1):

gu:=taylor(log(gu)/n^2 -log(1+v^2),v=0,2*K+1):

[seq(coeff(gu,v,2*i),i=1..K)]:

end:

#vTozNew(K):  Like zSeriesFE(K), but using the v-series (that we found from scratch!)
vTozNew:=proc(K) local gu,f,mu,i,v,z:

gu:=vSeriesFEpureNew(K):
f:=add(gu[i]*v^(2*i),i=1..K):

mu:=solve(2*v*(1-v)*(1+v)/(1+v^2)^2=z,v):

f:=taylor(subs(v=mu,f), z=0,2*K+1):

[seq(coeff(f,z,2*i),i=1..K)]:

end: