(* :Author:
	P.H. Lundow
	Bug reports to phl@kth.se
*)
(* :Summary:
	Implementation of Beale's method for computing Ising partition 
	functions.
*)
(* :History:
	070628: Slight speed-up, added IsingFreeEnergy.
*)
(* :References: 
	P. D. Beale,
	Exact distribution of energies in the two-dimensional Ising model,
	Phys. Rev. Lett. 76 (1996), 78--81
*)


BeginPackage["Beale`"]

IsingFreeEnergy::usage = 
	"IsingFreeEnergy[m, n, K] returns the free energy for the 
	mxn-torus at inverse temperature K."

IsingPoly::usage = 
	"IsingPoly[m, n, x] returns the Ising partition function for the
	mxn-torus in variable x."

Begin["`Private`"]

definitions = {IsingFreeEnergy, IsingPoly}

Scan[Unprotect, definitions]


IsingFreeEnergy[m_Integer, n_Integer, K_] :=
	Module[{sum,a,b,c0,cn,c2,s0,sn,s2,z1,z2,z3,z4,j,k,x},
		x = Exp[-2*K];
		c0 = (1-x)^m + (x*(1+x))^m;
		s0 = (1-x)^m - (x*(1+x))^m;
		cn = (1+x)^m + (x*(1-x))^m;
		sn = (1+x)^m - (x*(1-x))^m;
		b = 2 * x * (1-x^2);
		a[k_] = (1+x^2)^2 - b * Cos[Pi * k/n];
		sum[k_] = Sum[
			Binomial[m, 2j] * (a[k]^2 - b^2)^j * a[k]^(m-2j),
			{j, 0, Floor[m/2]}
		];
		c2[k_] = 2^(1-m) * (sum[k] + b^m);
		s2[k_] = 2^(1-m) * (sum[k] - b^m);
		If[Mod[n, 2] == 0,
			z1 = Product[c2[2k+1],{k,0,n/2-1}];
			z2 = Product[s2[2k+1],{k,0,n/2-1}];
			z3 = c0 * cn * Product[c2[2k],{k,1,n/2-1}];
			z4 = s0 * sn * Product[s2[2k],{k,1,n/2-1}],
		(*Else*)
			z1 = cn * Product[c2[2k+1],{k,0,(n-3)/2}];
			z2 = sn * Product[s2[2k+1],{k,0,(n-3)/2}];
			z3 = c0 * Product[c2[2k],{k,1,(n-1)/2}];
			z4 = s0 * Product[s2[2k],{k,1,(n-1)/2}]	  
		];
		Log[(z1+z2+z3+z4)/2]/m/n+2*K
	]


IsingPoly[m_Integer, n_Integer, x_] :=
	Module[{acc,sum,a,b,c0,cn,c2,s0,sn,s2,z,z1,z2,z3,z4,j,k},
		acc = Floor[N[m * n * Log[10,2]]] + 20; 
		c0 = (1-x)^m + (x*(1+x))^m;
		s0 = (1-x)^m - (x*(1+x))^m;
		cn = (1+x)^m + (x*(1-x))^m;
		sn = (1+x)^m - (x*(1-x))^m;
		b = 2 * x * (1-x^2);
		a[k_] = (1+x^2)^2 - b * Cos[Pi * k/n];
		sum[k_] = Sum[
			Binomial[m, 2j] * (a[k]^2 - b^2)^j * a[k]^(m-2j),
			{j, 0, Floor[m/2]}
		];
		c2[k_] = 2^(1-m) * (sum[k] + b^m);
		s2[k_] = 2^(1-m) * (sum[k] - b^m);
		If[Mod[n, 2] == 0,
			z1 = Product[c2[2k+1],{k,0,n/2-1}];
			z2 = Product[s2[2k+1],{k,0,n/2-1}];
			z3 = c0 * cn * Product[c2[2k],{k,1,n/2-1}];
			z4 = s0 * sn * Product[s2[2k],{k,1,n/2-1}],
		(*Else*)
			z1 = cn * Product[c2[2k+1],{k,0,(n-3)/2}];
			z2 = sn * Product[s2[2k+1],{k,0,(n-3)/2}];
			z3 = c0 * Product[c2[2k],{k,1,(n-1)/2}];
			z4 = s0 * Product[s2[2k],{k,1,(n-1)/2}]	  
		];
		z = N[(z1 + z2 + z3 + z4) / 2, acc];
		Dot[
			Round[CoefficientList[z, x]],
			x^Range[2*m*n, -2*m*n, -2]
		]
	]


Scan[Protect, definitions]

End[]

EndPackage[]