(* Section 2 *)
Print["Section 2"];
(* isosceles triangle with vertices (0,0), (1,0) and (Sqrt[3],0) *)
g[x_,y_]=Sin [Sqrt[3]Pi y]Sin[Pi x/3] + Sin[Pi y/Sqrt[3]]Sin[5Pi x/3] + Sin[2Pi y/Sqrt[3]]Sin[4Pi x/3];
(* other right triangles *)
g2[x_,y_]=g[1-x,y];          
g3[x_,y_]=g[Sqrt[3]y,Sqrt[3]x];
(* test functions obtained from right triangles *)
Psi1=g[x-(a y /b),Sqrt[3]y/b];
Psi2=g2[x-((a-1) y /b),Sqrt[3]y/b];
Psi3=g3[x-(a y /b),y/(Sqrt[3]b)];
(* equilateral triangle after linear transformation *)
Phi:=Sin[2Pi y/b]-Sin[2Pi(x+(1-a)y/b)]+Sin[2Pi(x-a y/b)]; 
(* final test function *)
test=A Psi1 + B Psi2 + G Psi3 + E Phi;
test // InputForm // Print;

(* Fast integration *)
(* Templates for integrals *)
int[a_,b_]=Integrate[Sin[A+B x],{x,a,b}];
s1 = (X_. Sin[A_.+B_ x]->X int[a y/b,(a-1)y/b+1]);
s2 = (X_. Cos[A_.+B_ x]->(s1[[2]]/.Cos[A_]->-Sin[A]));
s3 = (X_. Sin[A_.+B_ y]->X (int[0,b]/.x->y));
s4 = (X_. y Sin[A_.+B_ y]->X Integrate[y Sin[A+B y],{y,0,b}]);
s5 = (X_. Cos[A_.+B_ y]->(s3[[2]]/.Cos[A_]->-Sin[A]));
s6 = (X_. y Cos[A_.+B_ y]->(s4[[2]]/.{Cos[A_]->-Sin[A],Sin[A_]->Cos[A]}));
(* Expand an integral into a sum of trigonometric functions *)
(* and apply templates to each term *)
TrigInt[fff_]:=(
ff=ExpandAll[fff];
ff=TrigReduce[ff];
ff=Expand[ff];
ff=ff/.Sin[A_]:>Sin[Collect[A,{x},Simplify]]/.Cos[A_]:>Cos[Collect[A,{x},Simplify]];
ff=Replace[ff,{s1,s2,x_->(1 - y/b)x},1];
ff=ff/.Sin[A_]:>Sin[Collect[A,{y},Simplify]]/.Cos[A_]:>Cos[Collect[A,{y},Simplify]];
ff=Expand[ff];
ff=Replace[ff,{s4,s6,s3,s5,A_. y->A b^2/2,x_->b x},1]//Simplify;
ff);

int=TrigInt[test^2];
grad=TrigInt[D[test,x]^2+D[test,y]^2];

(* we have to prove that this is <= 0 *)
in=9b^2grad-4Pi^2(1+Sqrt[a^2+b^2]+Sqrt[(a-1)^2+b^2])^2int;
(* change from (a, b) to (M, N) and cancel b *)
in2=Simplify[in/b /. b^2 -> M^2 - a^2 /. a -> (M^2 - N^2 + 1)/2, (N > 0) && (M > 0)];
(* inequality (2.9) *)
Simplify[308788467187200in/b] // InputForm // Print;

(* Section 3 *)
Print["Section 3"];
W=in2/. E -> 0 /. G -> -1/6 /. B -> 0 /. A -> 1;
(* Inequality (3.1) *)
Apart[1383782400W] // InputForm // Print;
(* Critical point *)
Reduce[(D[W, M] == 0) && (D[W, N] == 0), {M, N}] // N // Print;
(* Boundary : roots and endpoints *)
Reduce[W == 0 /. N -> 2] // N // Print;
Reduce[W == 0 /. M -> N - 1] // N // Print;
Reduce[W == 0 /. M -> N] // N //Print;
Reduce[W == 0 /. M -> 15] // N // Print;
W /. M -> {1, 2} /. N -> 2 // N // Print;
W /. M -> 15 /. N -> {15, 16} // N //Print;

(* Section 4 *)
Print["Section 4"];
W=in2/. E -> 1 /. G -> 0 /. B -> A /.A -> (N + M - 2)/2;
pol = W/. M -> U - V /. N -> U + V /. U -> U + 1;
(* inequality (4.1) *)
Apart[22056319084800pol, V] // InputForm // Print;
(* Critical point *) 
Reduce[D[pol, V] == 0, V] // N // Print;
(* Boundary : roots and endpoints *)
Reduce[pol == 0 /. V -> 0] // N // Print;
Reduce[pol == 0 /. U -> 1 - V] // N // Print;
Reduce[pol == 0 /. U -> 3V] // N // Print;
pol /. V -> 0 /. U -> {0, 1} // N // Print;
pol /. V -> 1/4 /. U -> 3/4 // N // Print;

(* Section 5 *)
Print["Section 5"]
W=in2/. E -> 1 /. B -> 0 /. G -> A /.A -> (M + N - 2)/Sqrt[2];
pol = W /. M -> U - V /. N -> U + V /. U -> U + 1;
(* inequality (5.1) *)
Apart[9609600pol, V] // InputForm // Print;
(* Critical points *)
Vs = Solve[D[pol, V] == 0, V];
Reduce[D[pol, V] == 0, V, Reals] // InputForm // Print;
(* denominator with complx roots only*)
Reduce[Denominator[Together[D[pol, U] /. Vs]] == 0] // N  // Print;
(* polynomial of degree 7 in U *)
Reduce[Numerator[Together[D[pol, U] /. Vs]] == 0] // N // Print;
(* Boundary : roots and endpoints *)
Reduce[pol == 0 /. U -> 3V] // N // Print;
Reduce[pol == 0 /. U -> V] // N // Print;
Reduce[pol == 0 /. V -> 1 - U] // N // Print;
pol /. U -> 1 - V /. V -> {1/4, 1/2} // N // Print;
pol /. U -> 0 /. V -> 0 // N // Print;

Exit[]