(* Updae 2020-01-20: The variable Up was renamed in order to avoid a name conflict with a new built in definition  *)

(************ Definitions of functions *************)

Remove["Global`*"]
(* Applying the special product rule for α*)
Multiply[p_,c_]:=
Module[{i,j,m,b},
	FixedPoint[
		Expand[#]
		/.{
			Power[c[i_,m_],b_]:>(c[i,m])^Mod[b,2]*(2+c[2i,m])^Floor[b/2]/;b>=2,
			c[i_,m_]*c[j_,m_]:>c[i+j,m]+c[i-j,m]}&,
		p]
	]

(* Reducing by symmetries *)
SymReduce[p_,c_]:=
Module[{i,m},
	p/.c[i_,m_]:>c[-i,m]/;i<0
	//.c[i_,m_]:>c[i-2m,m]/;i>2m
	/.c[i_,m_]:>c[2m-i,m]/;i>m
	/.c[i_,m_]:>-c[m-i,m]/;2i>m
	/.{c[i_,m_]:>0/;2i==m,c[i_,m_]:>1/;3i==m,c[i_,m_]:>-1/;3i==2m,c[m_,m_]->-2,c[0,_]->2}
]

RemoveCos[p_,c_,acc_]:=
Module[{i,j,x},
	p/.c[i_,j_]:>N[2*Cos[i Pi/j],acc]/.x_Real:>Round[x]
]

TakeProduct[polys_,c_]:=
Module[{a,b,prod},
	Fold[SymReduce[Multiply[#1*#2,c],c]&,1,polys]
]

TakeProductII[polys_,c_]:=
Module[{a,b,prod,p},
	p=polys;
	While[Length[p]>1,
		p1=Take[p,Floor[Length[p]/2]];
		p2=Reverse[Take[p,-Floor[Length[p]/2]]];
		If[OddQ[Length[p]],p3={p[[Floor[Length[p]/2]+1]]},p3={}];
		p=MapThread[SymReduce[Multiply[#1*#2,c],c]&,{p1,p2}];
		p=Join[p,p3];
	];
	p[[1]]
]

FunctionPartition[l_List,f_]:=
Map[
	Last,
	Split[Sort[Transpose[{f[l],l}]],First[#1]==First[#2]&],
	{2}
]

U[K_]=FullSimplify[Coth[2K](1+2/Pi EllipticK[z^2](2*Tanh[2K]^2-1))/.z->2*Sinh[2K]/Cosh[2K]^2,Element[K,Reals]];

(************ Input *************)
Date[]
n=16  (** Should be even! **)

(*************** I **************)
p=n/2;
acc=Floor[N[p^2*Log[10,2]]];

Ti=FunctionPartition[Range[1,n,2],GCD[#,n]&];
Ui=2*FunctionPartition[Range[1,p-1],GCD[#,n]&];

Tp=RemoveCos[#,c,acc]&/@
	(TakeProductII[#,c]&/@
	Map[x-c[#,n]&,Ti,{2}]);
	
Upp=RemoveCos[#,c,acc]&/@
	(TakeProductII[#,c]&/@
	Map[x-c[#,n]&,Ui,{2}]);

(************ Verifying that all is OK so far ******)
And[
	(Simplify[Expand[Times@@Tp]-2*ChebyshevT[p,x/2]])==0,
	(Simplify[Expand[Times@@Upp]-ChebyshevU[p-1,x/2]])==0
]

(************* P[i,j] ***************)
P4polys=Map[
	(P[Y_,a_]=SymReduce[Multiply[#/.x->Y+c[a,n],c],c];
	RemoveCos[TakeProductII[#,c]&/@Map[P[Y,#]&,Ti,{2}],c,acc])&,
	Tp
];

P2polys=Map[
	(P[Y_,a_]=SymReduce[Multiply[#/.x->Y+c[a,n],c],c];
	RemoveCos[TakeProductII[#,c]&/@Map[P[Y,#]&,Ti,{2}],c,acc])&,
	Upp
];

P1polys=Map[
	(P[Y_,a_]=SymReduce[Multiply[#/.x->Y+c[a,n],c],c];
	RemoveCos[TakeProductII[#,c]&/@Map[P[Y,#]&,Ui,{2}],c,acc])&,
	Upp
];

(************* A1, A2, A4 *************)
A4=Expand[(Cancel[b^(p^2)Expand[Times@@Flatten[P4polys]]/.Y->a^2/b]/.{a->1+z^2,b->z(1-z^2)})^2];

P2prod=Expand[(Cancel[b^(p(p-1))Expand[Times@@Flatten[P2polys]]/.Y->a^2/b]/.{a->1+z^2,b->z(1-z^2)})^2];

P2fudge=Expand[Cancel[b^(2p)Expand[4*ChebyshevT[p,(Y+2)/2]ChebyshevT[p,(Y-2)/2]]/.Y->a^2/b]/.a->1+z^2/.b->z(1-z^2)];

A2=Expand[P2prod*P2fudge];

P1prod=Expand[(Cancel[b^((p-1)^2)Expand[Times@@Flatten[P1polys]]/.Y->a^2/b]/.{a->1+z^2,b->z(1-z^2)})^2];

P1fudge=Expand[(1+z^2)^2(1-6z^2+z^4)Cancel[b^(4p-4)Expand[ChebyshevU[p-1,(Y+2)/2]^2*ChebyshevU[p-1,(Y-2)/2]^2]/.Y->a^2/b]/.a->1+z^2/.b->z(1-z^2)];

A1=Expand[P1prod*P1fudge];

(**************** Z ***************)
Z=Expand[(A1+2A2+A4)/2];

Date[]

(************** Checking the moments *************)

s1=Simplify/@Integrate[Series[U[K],{K,0,n-1}],K];
s2=Simplify/@Series[Exp[n^2*s1],{K,0,n-1}];
s3=Simplify/@Series[Z/(2z)^(n^2)/.z->Exp[K]^2,{K,0,n-1}];

s2-s3==O[K]^n