######################################################################
#Lars.txt: Save this file as  Lars.txt                               #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read Lars.txt<Enter>                               #              
##Then follow the instructions given there                           #
##(based on joint work with Manuel Kauers)                           #                                                                    
##Written by Doron Zeilberger Rutgers University ,                   #
#DoronZeil at gmail dot com                                          #
######################################################################
 
#Created:  March-May 2018
 

read `Ons48PC.txt`:

print(`Created:  March-May 2018`):
print(`Version of May 1, 2018`):
print(` This is Lars.txt `):
print(`It  accompanies the article `):
print(`A Simple Re-Derivation of Onsager's Solution of the 2D Ising Model using Experimental Mathematics `):

print(`by Manuel Kauers and Doron Zeilberger`):

print(`and also available from Zeilberger's website`):
print(``):
print(`Please report bugs to DoronZeil at gmail dot com `):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://sites.math.rutgers.edu/~zeilberg/  .`):

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

print(`---------------------------------------`):
 print(`For a list of the direct procedures type ezraD();, 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 Cheating procedures type ezraCh();, 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):

ezraC:=proc()

if args=NULL then
 print(` The Checking procedures are: Check1, CheckDuality, CheckDualityS, CheckT27, `):
 print(``):

else
ezra(args):
fi:

end:

ezraCh:=proc()

if args=NULL then
 print(` The Cheating procedures (using the known Onsager solution) are: OnsPC, vSeriesCh `):
 print(``):

else
ezra(args):
fi:

end:

ezraD:=proc()

if args=NULL then
 print(` The Direct (very slow) procedures are: vSeriesD`):
 print(``):

else
ezra(args):
fi:

end:

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are:   Gfun, GP, GR, IsingPolST, IsingPolSTseq, IsingPolSTseqD, `):
 print(` MKa, Mult1,MMult1, OnsFun1, OnsFun2, PmnMK, Trace1, V, ZmnSeq, ZmnSeqD `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are:  IsingPols,  IsingPolsPC,  Mm, Mm2, OnsAI, PmnMK, Pmn,    OnsFun, OnsFunAI, vSeries, Zmn, zSeries `):
 print(` `):



elif nops([args])=1 and op(1,[args])=Check1 then
print(`Check1(x): Checks that [log(sinh(2*nu)),log(2*cosh(2*nu))-log(2*cosh(2*nu1)) ] are equal :`):
print(`where nu=log(x) and  x1:=sqrt((x^2+1)/(x^2-1)): and nu1=log(x1). Try:`):
print(`Check1(1.3);`):

elif nops([args])=1 and op(1,[args])=CheckDuality then
print(`CheckDuality(x): Checks the duality functional equation for the Onsager function  OnsFun(x), in the Direct Form.`):
print(`OnsFun(x*)=OnsFun(x)+log(2*x^2/(x^4-1)). Try:`):
print(`CheckDuality(1.2);`):

elif nops([args])=1 and op(1,[args])=CheckDualityS then
print(`CheckDualityS(x): Checks the duality functional equation for the Onsager function  OnsFun(x), in the SYMMETRIZED form.`):
print(`OnsFun(x*)- log(x*^2+1/x*^2)= OnsFun(x)-log(x^2+1/x^2):`):
print(`CheckDualityS(1.2);`):


elif nops([args])=1 and op(1,[args])=CheckT27 then
print(`CheckT27(m,n,x) ; checks numerically equation (2.7) for the polynomial Pmn(m,n,x), in Thompson's book, p. 153. Try`):
print(`CheckT27(4,5,2.1);`):

elif nops([args])=1 and op(1,[args])=Gfun then
print(`Gfun(z,eps): An approximation (with error err) of Sum(-(2*r-1)!^2/(r!^3*(r-1)!)*z^(2*r),r=1..infinity).`):
print(`Try: `):
print(` Gfun(2.4,10^(-12)); `):

elif nops([args])=1 and op(1,[args])=GP then
print(`GP(L,n): guesses a polynomial P of n,  such that P(L[i][1])=L[i][2].`):
print(`Try: GP([seq([i,i^2],i=1..4)],n);`):

elif nops([args])=1 and op(1,[args])=GR then
print(`GR(L,n): guesses a rational function R of n,  such that R(L[i][1])=L[i][2].`):
print(`Try: GR([seq([i,(i+2)/(i+3)],i=1..4)],n);`):
print(`Also Try: mu:=zSeries(16): GR([seq([i,mu[i+1]/mu[i]],i=1..nops(mu)-1)],n);`):

elif nops([args])=1 and op(1,[args])=IsingPols then
print(`IsingPols(k,N): Inputs a positive even integer, at least 4, and a symbol N.`):
print(`Outputs the list of the first k Ising polynomials in the variable N. Recall that the k-th Ising polynomial is`):
print(`the polynomial in N describing the number of ways, in a torodial strip with N sites (with N larger than 2k+2), of`):
print(`picking 2k bonds such that every site has an even number of neighbors.`):
print(` Try: `):
print(` IsingPols(4,N); `):

elif nops([args])=1 and op(1,[args])=IsingPolsPC then
print(`IsingPolsPC(N):  the output of IsingPols(16,N). `):
print(`Try: `):
print(` IsingPolsPC(N); `):

elif nops([args])=1 and op(1,[args])=IsingPolST then
print(`IsingPolST(m,k,N): the coefficient of v^k, as a polynomial in N=m*n, in Zmn(m,n) for n>k. Try: `):
print(`IsingPolST(4,4,N);`):

elif nops([args])=1 and op(1,[args])=IsingPolSTseq then
print(`IsingPolSTseq(m,k,N): Fast version. The Ising polynomials for torodial strip of width m from the coeff. of v^1 to the coeff.of v^k.`):
print(` m must be at least 3. `):
print(`Try: `):
print(` IsingPolSTseq(4,6,N); `):

elif nops([args])=1 and op(1,[args])=IsingPolSTseqD then
print(`IsingPolSTseqD(m,k,N): the Ising polynomials for torodial strip of width m from the coeff. of v^1 to the coeff.of v^k.`):
print(`Done directly, for checking purposes only.`):
print(` m must be at least 3. `):
print(`Try: `):
print(` IsingPolSTseqD(4,6,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])=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 magnetization. Try: `):
print(`Mm(2,x1,x2,y);`):

elif nops([args])=1 and op(1,[args])=Mm2 then
print(`Mm2(m,v,y):The 2^m by 2^m transfer matrix with`):
print(`x for vertical and horiz. interaction `):
print(`and y for magnetization, BUT in terms of v. Try:`):
print(`Mm2(2,v,y);`):

elif nops([args])=1 and op(1,[args])=Mult1 then
print(`Mult1(P,Q,x,d): the product of two polynomials P and Q of the variable x trncated to x^d. Try`):
print(`Mult1(1+x+x^2,1+x+x^2,x,3);`):

elif nops([args])=1 and op(1,[args])=MMult1 then
print(`MMult1(P,Q,x,d): given two matrices of polynomials P and Q in the variable x, finds the matrix product P, Q up to the k-th power. Try:`):
print(`MMult1([[1+x+x^2,1+x+x^2], [1+x+x^2,1+x+x^2]],[[1+x+x^2,1+x+x^2], [1+x+x^2,1+x+x^2]], x,3);`):

elif nops([args])=1 and op(1,[args])=OnsPC then
print(`OnsPC(x):  the list of length 12 whose i-th entry is Pmn(4*i,4*i,x,1)`):


elif nops([args])=1 and op(1,[args])=OnsAI then
print(`OnsAI(v,DP): getting the Onsager expression for -psi/(k*T) by pure experimental mathematics, Ab Initio, using DP date points`):
print(`for the so-called a[r]^(1) in the combinatorial approach. Here v=exp(J/(k*T)), J is the coupling constant, k is Avogardro's number`):
print(`and T is the absolute temperature. `):
print(` Try:`):
print(`OnsAI(v,10);`):


elif nops([args])=1 and op(1,[args])=OnsFun then
print(`OnsFun(x): the Onsager solution, found by Onsager evaluated at x (=exp(nu)), using  equation (Chatper 5, (5.8)) (p. 132) in Thompson's book. Try:`):
print(`plot([seq([1+0.01*i,OnsFun(1+0.01*i)],i=1..50)]);`):

elif nops([args])=1 and op(1,[args])=OnsFun1 then
print(`OnsFun1(x): the Onsager solution, found by Onsager evaluated at x (=exp(nu)), using  equation (Chapter 5, (5.4)) p. 131) in Thompson's book. Try:`):
print(`plot([seq([1+0.01*i,OnsFun1(1+0.01*i)],i=1..50)]);`):

elif nops([args])=1 and op(1,[args])=OnsFun2 then
print(`OnsFun2(x,N): the Onsager solution, found by Onsager evaluated at x (=exp(nu)), using  equation (Chatper 5, (5.8)) (p. 132) in Thompson's book.`):
print(`but using the implied Taylor expansion instead of integrating. Try:`):
print(`plot([seq([1+0.01*i,OnsFun2(1+0.01*i,10)],i=1..50)]);`):

elif nops([args])=1 and op(1,[args])=OnsFunAI then
print(`OnsFunAI(x,eps): the Onsager solution (with error eps) via the new (Ab Initio) approach. Try:`):
print(`plot([seq([1+0.01*i,OnsFunAI(1+0.01*i,10^(-12))],i=1..50)]);`):

elif nops([args])=1 and op(1,[args])=PmnMK then
print(`PmnMK(m,n,x,y): the trace of Mm(m,x,x,y)^n, using Manuel Kauers' fast algorithm. Try:`):
print(`PmnMK(2,2,x,y);`):

elif nops([args])=1 and op(1,[args])=Pmn then
print(`Pmn(m,n,x,y): the trace of Mm(n,x,x,y)^n  done directly. Try:`):
print(`Pmn(2,2,x,y);`):

elif nops([args])=1 and op(1,[args])=Trace1 then
print(`Trace1(P): the trace of the square matrix P. Try: `):
print(` Trace1([[1,2],[2,3]]); `):


elif nops([args])=1 and op(1,[args])=V then
print(`V(i,m): the binary representation`):
print(`of i between 0 and 2^m-1 given as a list of 0's and 1's`):

elif nops([args])=1 and op(1,[args])=vSeries then
print(`vSeries(K):  The v Series for the Ising model, via the Ising Polynomials. Try`):
print(`vSeries(16);`):

elif nops([args])=1 and op(1,[args])=vSeriesCh then
print(`vSeriesCh(K):  Like vSeriesD(K),  but using pre-computed data for the Ising Generating functions given in the`):
print(`data file (also available from the front of this article) Ons48PC.txt .`):
print(`K must be at most 23. Try:`):
print(`vSeriesCh(20);`):

elif nops([args])=1 and op(1,[args])=vSeriesD then
print(`vSeriesD(K):  The v Series for the Ising model. Done directly. VERY SLOW. Try`):
print(`vSeriesD(2);`):

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])=ZmnSeq then
print(`ZmnSeq(m,v,n0,d): Inputs a pos. integer m, a variable name v, and positive integers n0 and d, and outputs`):
print(`the list of length n0 whose n-th entry is the polynomial in v, of degree k, that is the truncation (up to v^d)`):
print(`of the polynomial Zmn(m,n,v). Try:`):
print(`ZmnSeq(4,v,8, 6);`):

elif nops([args])=1 and op(1,[args])=ZmnSeqD then
print(`ZmnSeqD(m,v,n0,d):  Like ZmnSeq(m,v,n0,d), but done directly. Try: `):
print(`ZmnSeqD(4,v,8, 6);`):


elif nops([args])=1 and op(1,[args])=zSeriesD then
print(`zSeriesD(K):  The z Series for the Ising model. Try`):
print(`zSeriesD(2);`):

elif nops([args])=1 and op(1,[args])=zSeries then
print(`zSeries(K):  The z Series for the Ising model, via the Ising Polynomials. Try`):
print(`zSeries(16);`):


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

###start guess rational
#GR1(L,n,d): guesses a rational function R of n,  of degree d, such that R(L[i][1])=L[i][2].
#Try: GR1([seq([i,(i+2)/(i+3)],i=1..4)],n,1);
GR1:=proc(L,n,d) local i,eq,var,a,b,R:

if nops(L)<2*d+2 then
 print(L, `should have at least`, 2*d+2, `members. `):
 RETURN(FAIL):
fi:

R:=add(a[i]*n^i,i=0..d)/(1+add(b[i]*n^i,i=1..d)):

eq:={seq(subs(n=L[i][1],R)-L[i][2],i=1..2*d+2)};

var:={seq(a[i],i=0..d),seq(b[i],i=1..d)}:

var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
fi:

R:=factor(subs(var,R)):

if {seq(subs(n=L[i][1],R)-L[i][2],i=2*d+2..nops(L))}<>{0} then
 print(R, `did not work out`):
 RETURN(FAIL):
fi:

R:

end:





#GR(L,n): guesses a rational function R of n,  such that R(L[i][1])=L[i][2].
#Try: GR([seq([i,(i+2)/(i+3)],i=1..4)],n);
GR:=proc(L,n) local d,gu:

for d from 1 to trunc((nops(L)-2)/2) do
 gu:=GR1(L,n,d):
 if gu<>FAIL then
   RETURN(gu):
 fi:
od:

gu:

end:

###end guess rational

#OnsFun1(x): the Onsager solution, found by Onsager, using  equation (5.4) (Chapter 5) in Thompson's book. Try: 
#plot([seq([1+0.01*i,OnsFun1(1+0.01*i)],i=1..100)]);
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)+1/2*log(x^2-1/x^2)):
end:


#OnsFun2(x,N): the Onsager solution, found by Onsager, using  the Taylor expansion in z up to N terms
#plot([seq([1+0.01*i,OnsFun2(1+0.01*i,20)],i=1..100)]);
OnsFun2:=proc(x,N) local t1,t2,gu,z,n:
z:=2*(x^2-1/x^2)/(x^2+1/x^2)^2:
evalf(log(x^2+1/x^2))-add(binomial(2*n,n)^2/n/4^(n+1)*z^(2*n),n=1..N):
end:

###start guess poly
#GP1(L,n,d): guesses a polynomial R of n,  of degree d, such that R(L[i][1])=L[i][2].
#Try: GP1([seq([i,i^3+1],i=1..4)],n,1);
GP1:=proc(L,n,d) local i,eq,var,a,b,R:

if nops(L)<d+2 then
 print(L, `should have at least`, d+2, `members. `):
 RETURN(FAIL):
fi:

R:=add(a[i]*n^i,i=0..d):

eq:={seq(subs(n=L[i][1],R)-L[i][2],i=1..d+2)};

var:={seq(a[i],i=0..d)}:

var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
fi:

R:=factor(subs(var,R)):

if {seq(subs(n=L[i][1],R)-L[i][2],i=d+3..nops(L))}<>{0} then
 print(R, `did not work out`):
 RETURN(FAIL):
fi:

R:

end:





#GP(L,n): guesses a polynomial R of n,  such that R(L[i][1])=L[i][2].
#Try: GP([seq([i, i^3],i=1..5)],n);
GP:=proc(L,n) local d,gu:


for d from 0 to trunc((nops(L)-2)/2) do
 gu:=GP1(L,n,d):
 if gu<>FAIL then
   RETURN(gu):
 fi:
od:

gu:

end:

###end guess poly


#V(i,m): the binary representation
#of i between 0 and 2^m-1 given as a list of 0's and 1's
V:=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(V(i,m-1))]:
else
  [1,op(V(i-2^(m-1),m-1))]:
fi:
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:=V(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:

#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 0 to 2^m-1 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 0 to 2^m-1 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=0..2^m-1)],i=0..2^m-1)]:

end:

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

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

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

end:


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

for i from 0 to 2^m-1 do
  v1:=V(i,m):
 gu:=y^convert(v1,`+`)*x^(add(v1[i1]*v1[i1+1],i1=1..m-1)+v1[m]*v1[1]):
 for j from 0 to 2^m-1 do
   v2:=V(j,m):
   T[i,j]:=gu*x^(add(v1[i1]*v2[i1],i1=1..m))/((x+1/x)^2/2)^m:
 T[i,j]:=simplify(subs(x=sqrt((1+v)/(1-v)),T[i,j])):
 od:
od:

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

end:


#Pmn(m,n,x,y): the trace of Mm(n,x,x,y)^m done directly
#Pmn(2,2,x,y);
Pmn:=proc(m,n,x,y) local t,i,v,k,gu:

gu:=Matrix(Mm(m,x,x,y)):
expand(LinearAlgebra[Trace](gu^n)):
end:

#Pmn1(m,n,x,y): the trace of Mm1(n,x,x,y)^m done directly
#Pmn1(2,2,x,y);
Pmn1:=proc(m,n,x,y) local t,i,v,k,gu:

gu:=Matrix(Mm1(m,x,y)):
expand(LinearAlgebra[Trace](gu^n)):
end:

#PmnMK(m,n,x,y): the trace of Mm(n,x,x,y)^m 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:


#Zmn(m,n,v): the expression in v=(x-1/x)/(x+1/x) of Pmn(m,n,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,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:


#Zmn1(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:
#Zmn1(4,4,v);
Zmn1:=proc(m,n,v) local gu,x,N,lu,mu:
option remember:

gu:=Pmn1(m,n,x,1):
mu:=simplify(subs(x=sqrt((1+v)/(1-v)),gu)):
mu:=expand(mu):


end:

#vSeriesD(K):  vSeries(K) done directly
#vSeriesD(20);
vSeriesD:=proc(K) local gu,n,i,v:

n:=2*K+1:
 gu:=Zmn(n,n,v):

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

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

end:

#vSeriesCh(K):  Like vSeriesD(K),  but using pre-computed data for the Ising Generating functions given in the
#data file (also available from the front of this article) Ons48PC.txt .
#K must be at most 23. Try:
#vSeriesCh(20);
vSeriesCh:=proc(K) local x,gu,mu,n,N,i,v:

if K>22 then
 print(`The input must be<=22`):
 RETURN(FAIL):
fi:

n:=2*K+1:

for i from 1 while 4*i<n do od:



if i>12 then
  RETURN(FAIL):
fi:

 gu:=OnsPC(x)[i]:

n:=4*i:
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):

n:=4*i:

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

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

end:








#Zmn2(m,n,v): Like Zmn(m,n,v) but via Mm2(m,v,1)
#Try:
#Zmn2(4,4,v);
Zmn2:=proc(m,n,v) local gu:
option remember:

gu:=Matrix(Mm2(m,v,1)):

expand(LinearAlgebra[Trace](gu^n)):

end:


#IsingPolST(m,k,N): the coefficient of v^k, as a polynomial in N=m*n, in Zmn(m,n) for n>k. Try
#IsingPolST(4,4);
IsingPolST:=proc(m,k,N) local gu,n,v:
gu:=GP([seq([n,coeff(Zmn(m,n,v),v,k)],n=k+2..2*k+2)],n):
subs(n=N/m,gu):
end:

###start truncated polynomials

#Mult1(P,Q,x,d): the product of two polynomials P and Q of the variable x trncated to x^d. Try
#Mult1(1+x+x^2,1+x+x^2,x,3);
Mult1:=proc(P,Q,x,d) local i,gu:
gu:=P*Q:
gu:=taylor(gu,x=0,d+1):
add(coeff(gu,x,i)*x^i,i=0..d):
end:

#MMult1(P,Q,x,d): given two matrices of polynomials P and Q in the variable x, finds the matrix product P, Q up to the k-th power. Try:
#MMult1([[1+x+x^2,1+x+x^2], [1+x+x^2,1+x+x^2]],x,3);
MMult1:=proc(P,Q,x,d) local i,j,k:
if nops(P[1])<>nops(Q) then
  RETURN(FAIL):
fi:

[seq([seq(add(Mult1(P[i][k],Q[k][j],x,d),k=1..nops(Q)),j=1..nops(Q))],i=1..nops(P))]:
end:

#Trace1(P): the trace of the square matrix P. Try
#Trace1([[1,2],[2,3]]);
Trace1:=proc(P) local i:

if nops(P)<>nops(P[1]) then
 RETURN(FAIL):
fi:

expand(add(P[i][i],i=1..nops(P))):

end:

Trunc1:=proc(P,x,d) local i:
 add(coeff(P,x,i)*x^i,i=0..d):
end:

#ZmnSeq(m,v,n0,d): Inputs a pos. integer m, a variable name v, and positive integers n0 and d, and outputs
#the list of length n0 whose n-th entry is the polynomial in v, of degree k, that is the truncation (up to v^d)
#of the polynomial Zmn(m,n,v). Try
#ZmnSeq(4,v,10,14);
ZmnSeq:=proc(m,v,n0,d) local M1, M,gu,gu1,n1:

M1:=expand(Mm2(m,v,1)):
M:=M1:

gu1:=Trunc1(Trace1(M1),v,d):

gu:=[gu1]:

for n1 from 2 to n0 do
 M:=MMult1(M,M1,v,d):
 gu:=[op(gu),Trace1(M)]:
od:

gu:

end:



#ZmnSeqD(m,v,n0,d): Like ZmnSeq(m,v,n0,d) but done directly. Try:
#ZmnSeqD(4,v,10,14);
ZmnSeqD:=proc(m,v,n0,d) local n1:

[seq(Trunc1(Zmn(m,n1,v),v,d),n1=1..n0)]:

end:



#IsingPolSTseqD(m,k,N): the Ising polynomials for torodial strip of width m from the coeff. of v to the coeff.of v^k.
#Try:
#IsingPolSTseqD(4,6,N);
IsingPolSTseqD:=proc(m,k,N) local k1:
if m<3 then
 print(m ,` must be at least 3`):
 RETURN(FAIL):
fi:

[0,seq(IsingPolST(m,k1,N) ,k1=2..k)]:
end:

#IsingPolSTseq(m,k,N): the Ising polynomials for torodial strip of width m from the coeff. of v to the coeff.of v^k. Fast way
#Try:
#IsingPolSTseq(4,6,N);
IsingPolSTseq:=proc(m,k,N) local k1,n1,mu,gu,gu1,n,v,lu:
option remember:

if m<3 then
 print(m ,` must be at least 3`):
 RETURN(FAIL):
fi:

mu:=ZmnSeq(m,v,2*k+2,k):


gu:=[0]:

for k1 from 2 to k do
 lu:=[seq([n1,coeff(mu[n1],v,k1)],n1=k1+2..2*k+2)]:

 gu1:=GP(lu,n):
  gu1:=subs(n=N/m,gu1):
 gu:=[op(gu),gu1]:

od:

gu:

end:

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:




#IsingPols(k,N): Inputs a positive even integer, at least 4, and a symbol N.
#Outputs the list of the first k Ising polynomials in the variable N. Recall that the k-th Ising polynomial is
#the polynomial in N describing the number of ways, in a torodial strip with N sites (with N larger than 2k+2), of
#picking 2k bonds such that every site has an even number of neighbors.
#Try:
#IsingPols(4,N);
IsingPols:=proc(k,N) local m,gu,lu,i:
option remember:
if k<4 then
 print(k ,` must be at least 4`):
 RETURN(FAIL):
fi:

 if k mod 2<>0 then
  print(k, `must be even`):
 RETURN(FAIL):
 fi:
m:=k-1:
gu:=IsingPolSTseq(m,2*k,N):

lu:=[seq(gu[2*i],i=1..m-1)]:
lu:=[op(lu),gu[2*m]-N/m*(N/m-1)/2]:
lu:=[op(lu),gu[2*m+2]-N*((m-1)/m*N-(2*m-1))]:
factor(lu):
end:



#vSeries(K):  The first K terms of the free energy in terms of w^2=tanh(J/(k*T)), using the Ising Polynomials
#vSeries(16);
vSeries:=proc(K) local gu,N,i:

if K>16 then
 print(`Not yet implemented`):
 RETURN(FAIL):
fi:

gu:=IsingPolsPC(N):
[seq(coeff(gu[i],N,1),i=1..K)]:
end:

#zSeries(K):  The z-series using the Ising Polynomials
zSeries:=proc(K) local gu,f,mu,i,v,z,gu1:

gu:=vSeries(K):

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

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:






#CheckT27(m,n,x) ; checks numerically equation (2.7) for the polynomial Pmn(m,n,x), in Thompson's book, p. 153. Try
#CheckT27(4,5,2.1);
CheckT27:=proc(m,n,nu) local N,x,nu1,gu1,gu2:

 N:=m*n:

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

 gu1:=subs( x=exp(nu), Pmn(m,n,x,1))*exp(-2*N*nu):

 gu2:=2*subs( x=exp(nu1), Pmn(m,n,x,1))/2^N/cosh(nu1)^(2*N):

[gu1,gu2]:

end:



#OnsAI(v,DP): getting the Onsager expression for -psi/(k*T) by pure experimental mathematics, Ab Initio. 
#v=exp(J/(k*T)), and k is the number of terms used for guessing the z-zeries. Try:
#OnsAI(v,9);
OnsAI:=proc(v,DP) local psi,r,w,n,N,T,k,a,gu,r1,lu,F,lu1,lu2,mu,i,z,ku,G,b:
print(`An Experimental Mathematics Derivation of Onsager's famous Explicit formula for the Free Energy of the Ising Model`):
print(``):
print(`By Manuel Kauers' Computer and Doron Zeilberger's Computer [Shalosh B. EKhad] `):
print(``):
print(`Abstract: One of the most celebrated results in Statistical Physics is Lars Onsager's tour-de-force explicit formula`):
print(`for the free energy of the 2-dimensional Ising model with zero magnetic field.  (Physical Reviews v.65 (1944), p. 117ff). `):
print(`His proof, as well as all the subsequence proofs,`):
print(`are extremely complicated and ad hoc. In this note we will rederive his famous formula AB INITIO. While we use an experimental mathematics`):
print(`approach, based on guessing, that is non-rigorous, it is nevertheless absolutely certain, and should satisfy physicists who do not`):
print(`share the mathematicians' dogmatic insistence on rigor. Recall that many people got Physics Nobel prizes for  physical theories`):
print(`based on non-rigorous mathematics. `):
print(``):
print(`Let -psi/(k*T), the Onsager free energy, be denoted by Ons(v),`):
print(`where , v=J/(k*T),  T is the absolute temperature , and J is the coupling constant. `):
print(``):
print(`By the so-called (elementary!) combinatorial approach (due to B.L. van der Waerden, 1941),`):
print(`that is described, inter alia, in Colin Thompson's classic monograph "Mathematical Statistical Physics", MacMillan, 1972, `):
print(`[Note that this method predates  Onsager's 1944 proof]`):
print(``):
print(Ons(v)= log(2) +2*log(cosh(v))+ Sum(a[r]^({1})*w^r,r=1..infinity)):
print(``):
print(`where w=tanh(v). See Eq. (1.18) in Chapter 6 of Thompson's book (p. 150) . `):
print(``):
print( ` Here for each`,r,  a[r]^({1}), `is the coefficient of N in the polynomial expression  `, n[r](N), ` for the number `):
print(`of so-called lattice polygons with r edges that can be placed in a rectangular torus with N lattice points`):
print(`whose sides are larger than r. Lattice polygons are subgraphs of the rectangular lattice where every`):
print(`vertex has an even number of neighbors (i.e. 0,2, or 4) .`):
print(``):
print(`By using the (finite!) transfer matrices for width <=7, one can guess explicit expressions for the first`, DP, ` polynomials `):
gu:=IsingPolsPC(N):
gu:=[op(1..DP,gu)]:

for r1 from 1 to DP do 
 print(n[2*r1]=gu[r1]):
od:

print(``):
print(`of course n[odd] are all zero. `):

print(``):
print(` Taking the coefficients of N, we get that the first `, DP, ` terms of `, a[2*r]^({1}), ` are `):
print(``):
print( seq(a[2*r1]^({1})=coeff(gu[r1],N,1), r1=1..DP) ):
print(``):

print(`By the Kramres-Wannier duality argument, given a combinatorial explanation by van der Waerden in 1941, Ons(v) satisfies the`):
print(` FUNCTIONAL EQUATION`):
print(``):
print(`Ons(v*)=Ons(v)-log(sinh(2v)) `):
print(``):
print(`where v*=arctanh(exp(-2*v)) . `):
print(``):
print(`This is equivalent to`):
print(`Ons(v*)-log(2cosh(2v*))=Ons(v)-log(2cosh(2v)) .`):
print(` Since `):
print(``):
print(z=sinh(2*v)/cosh(2*v)^2):
print(``):
print(`is invariant under v<->v*, it is natural to consider the function`):
print(Ons(v)-log(2*cosh(2*v))  ):
print(``):
print(`and later express it in terms of z, that is the most natural variable.`):
print(``):
print(`This function equals :`):
print(``):
print(-log(2*cosh(2*v))+ log(2) +2*log(cosh(v))+ Sum(a[r]^({1})*w^r,r=1..infinity)):
print(``):
print(`where w=tanh(v) `):
print(``):
print(` Since `):
print(``):
print(-log(2*cosh(2*v))+ log(2) +2*log(cosh(v))=-log(1+w^2)):
print(``):
print(`this function, entirely in terms of w, is `):
print(``):
print(  -log(1+w^2)+ Sum(a[2*r]^({1})*w^(2*r),r=1..infinity)):
print(``):
lu1:=[seq(coeff(taylor(-log(1+w^2),w=0,4*DP+2),w,2*r1),r1=1..DP)]:
lu2:=[seq(coeff(gu[r1],N,1), r1=1..DP)]:
lu:=lu1+lu2:
F:=add(lu[r1]*w^(2*r1),r1=1..DP):
print(`whose beginning, up to, and including, the power`, w^(2*DP), `is `):
print(F+O(w^(2*DP+2)) ):
print(``):
print(`changing to the natural z variable, that in terms of w is `):
print(``):
print(z=2*w*(1-w)*(1+w)/(1+w^2)^2):
print(``):
mu:=solve(2*w*(1-w)*(1+w)/(1+w^2)^2=z,w):

F:=taylor(subs(w=mu,F), z=0,2*DP+1):

print(`we get that the beginning Taylor series is `):
print(``):
print(F):

lu:=[seq(coeff(F,z,2*i),i=1..DP)]:

print(`The first`, DP, `coefficients are `):
print(``):
print(lu):
print(``):
print(`factoring suggests that there is a closed form`):
print(``):
print(ifactor(lu)):
print(``):
print(`The sequence of  ratios is`):
print(``):
lu:=[seq(lu[i+1]/lu[i],i=1..nops(lu)-1)]:
print(``):
print(lu):
print(``):
print(`that obviously fits`):
print(``):
ku:=GR([seq([i,lu[i]],i=1..nops(lu)-1)],r):
print( ku):
print(``):
print(`Hence we discovered the recurrence for the coefficients: `):
print(``):
print(b[2*r+2]=b[2*r]*ku , `, `):
print(``):
print(`that immediately implies the closed-form expression for`, b[2*r] , `r>=1 `):
print(``):
print(b[2*r]=-binomial(2*r,r)^2/(r*4^(r+1))):
print(``):
print(`and indeed the first few terms are `):
print(``):
print(seq(-binomial(2*r1,r1)^2/(r1*4^(r1+1)),r1=1..DP)):
print(``):
print(`Hence we discovered`):
print(``):
print(`Theorem: Let G(z) be`):
print(``):
print(Sum(-binomial(2*r,r)^2/(r*4^(r+1))*z^(2*r),r=1..infinity)):
print(``):
print(`Then `):
print(``):
print(-psi/(k*T)=log(2*cosh(2*v))+G(sinh(2*v)/cosh(2*v)^2)):

print(`----------------------------`):
print(`this ends our ab initio, non-rigorous, but absolutely certain,  derivation of Onsager's famous formula, by purely computational means`):
print(`and elementary combinatorial reasonings due to van der Waerden and Kramers and Wannier. `):
end:







#Gfun(z,eps): An approximation (with error err) of Sum(-(2*r-1)!^2/(r!^3*(r-1)!)*z^(2*r),r=1..infinity).
#Try:
#Gfun(2.4,10^(-12));

Gfun:=proc(z,eps) local gu,gu1,r:

if evalf(abs(z))>0.49 then
 RETURN(FAIL):
fi:


gu1:=evalf(subs(r=1, -binomial(2*r,r)^2/4^(r+1)/r *z^(2*r))):

gu:=gu1:

for r from 2 to 1000 while abs(gu1)>eps/1000 do

gu1:=evalf(-binomial(2*r,r)^2/4^(r+1)/r*z^(2*r)):

gu:=gu+gu1:
od:
if r>1000 then
 RETURN(FAIL):
fi:
gu:

end:

#OnsFun(x): the Onsager solution, found by Onsager, using  equation (5.8) (Chapter 5) in Thompson's book. Try: 
#plot([seq([1+0.01*i,OnsFun(1+0.01*i)],i=1..100)]);
OnsFun:=proc(x) local t1,t2,z:
z:=2*(x^2-1/x^2)/(x^2+1/x^2)^2:
evalf(log(x^2+1/x^2))+ evalf(1/(2*Pi^2))*evalf(Int(Int(log(1-z*(cos(t1)+cos(t2)))  ,t1=0..Pi),t2=0..Pi)):
end:

#OnsFunAI(x,eps): the Onsager solution (with error err) via the new approach. Try:
#plot([seq([1+0.01*i,OnsFunAI(1+0.01*i, 10^(-12))],i=1..50)]);
OnsFunAI:=proc(x,eps):
evalf(log(x^2+1/x^2))+Gfun(2*(x^2-1/x^2)/(x^2+1/x^2)^2,eps):

end:


#CheckDuality0(x): 
#CheckDuality0(1.2);
CheckDuality0:=proc(x) local nu:
nu:=log(x):
 evalf(OnsFun(sqrt((x^2+1)/(x^2-1)))-OnsFun(x)+log(sinh(2*nu)));end:

#CheckDuality(x): Checks the duality functional equation for the Onsager function  OnsFun(x), in the direct form.
#Let x*=sqrt((x^2+1)/(x^2-1)) then:
#OnsFun(x*)=OnsFun(x)+log(2*x^2/(x^4-1)). Try:
#CheckDuality(1.2);
CheckDuality:=proc(x): evalf(OnsFun(sqrt((x^2+1)/(x^2-1)))-OnsFun(x)-log(2*x^2/(x^4-1)));end:



#CheckDualityS(x): Checks the duality functional equation for the Onsager function  OnsFun(x), in the symmetrized form.
#Let x*=sqrt((x^2+1)/(x^2-1)) then:
#OnsFun(x*)- log(x*^2+1/x*^2)= OnsFun(x)-log(x^2+1/x^2):
CheckDualityS:=proc(x) local x1:
x1:=sqrt((x^2+1)/(x^2-1)):
[OnsFun(x1)- log(x1^2+1/x1^2), OnsFun(x)-log(x^2+1/x^2)]:
end:

#CheckDualityS0(x): Checks the duality functional equation for the Onsager function  OnsFun(x), in the symmetrized form.
#Let x*=sqrt((x^2+1)/(x^2-1)) then:
#OnsFun(x*)- log(x*^2+1/x*^2)= OnsFun(x)-log(x^2+1/x^2):
CheckDualityS0:=proc(x) local x1,nu,nu1:
nu:=log(x):
x1:=sqrt((x^2+1)/(x^2-1)):
nu1:=log(x1):
[OnsFun(x1)- log(2*cosh(2*nu1)), OnsFun(x)-log(2*cosh(2*nu))]:
end:



##############
#GfunS(z,N): The beginning of Sum(-(2*r-1)!^2/(r!^3*(r-1)!)*z^(2*r),r=1..infinity).
#Try:
#GfunS(z,N);

GfunS:=proc(z,N) local r:
add( -binomial(2*r,r)^2/4^(r+1)/r*z^(2*r),r=1..N):

end:





#Check1(x): Checks that [log(sinh(2*nu)),log(2*cosh(2*nu))-log(2*cosh(2*nu1)) ] are equal :
#where nu=log(x) and  x1:=sqrt((x^2+1)/(x^2-1)): and nu1=log(x1). Try:
#Check1(1.3);
Check1:=proc(x) local nu, x1,nu1:
nu:=log(x):
x1:=sqrt((x^2+1)/(x^2-1)):
nu1:=log(x1):

[log(sinh(2*nu)),log(2*cosh(2*nu))-log(2*cosh(2*nu1)) ]:

end:

#Check1a(x): Checks that [log(sinh(2*nu)),log(2*cosh(2*nu))-log(2*cosh(2*nu1)) ] are equal :
#where nu=log(x) and  x1:=sqrt((x^2+1)/(x^2-1)): and nu1=log(x1). Try:
#Check1a(1.3);
Check1a:=proc(x) local nu, x1,nu1:
nu:=log(x):
x1:=sqrt((x^2+1)/(x^2-1)):
nu1:=log(x1):

[sinh(2*nu),cosh(2*nu)/cosh(2*nu1) ]:

end:

#Check1b(x): Checks that [log(sinh(2*nu)),log(2*cosh(2*nu))-log(2*cosh(2*nu1)) ] are equal :
#where nu=log(x) and  x1:=sqrt((x^2+1)/(x^2-1)): and nu1=log(x1). Try:
#Check1b(x);
Check1b:=proc(x) local x1:

x1:=sqrt((x^2+1)/(x^2-1)):

normal((x^2-1/x^2)/2-(x^2+1/x^2)/(x1^2+1/x1^2)):


end:

#Check1c(X): Checks that [log(sinh(2*nu)),log(2*cosh(2*nu))-log(2*cosh(2*nu1)) ] are equal :
#where nu=log(x) and  x1:=sqrt((x^2+1)/(x^2-1)): and nu1=log(x1). Try:
#Check1c();
Check1c:=proc() local X,X1:

X1:=(X+1)/(X-1):

evalb(normal((X-1/X)/2-(X+1/X)/(X1+1/X1))=0):


end: