(* Fully automated verification of the ODE from Section 4: Proof of Lemma 4 *)
f[x_]:=((1-x)^(2/p)-x^(2/p))^(q/2)LegendreP[q/2,((1-x)^(2/p)+x^(2/p))/((1-x)^(2/p)-x^(2/p))];
a[x_]:=p^2(1-x)^2x^2((1-x)^(2/p)-x^(2/p));
b[x_]:=p (1-x)x((1-x)^(2/p)(p(1-2x)+2q x)+x^(2/p)(p(2x-1)+2q(1-x)));
c[x_]:=q(1-x)^(2/p)x(p(1-x)+q x)-q(1-x)x^(2/p)(p x+q(1-x));
df[x_]:=Simplify[D[f[x],{x,1}]];
ddf[x_]:=Simplify[D[f[x],{x,2}]];
FullSimplify[a[x] ddf[x]+b[x]df[x] +c[x]f[x]]

(* Half-automated verification of the ODE: using only symbolic algebraic simplification*)
(* g stands for the Legendre function P_{q/2} *)
f[x_]:=((1-x)^(2/p)-x^(2/p))^(q/2)g[((1-x)^(2/p)+x^(2/p))/((1-x)^(2/p)-x^(2/p))];
(* h stands for P_{q/2}, dh stands for P'_{q/2}, ddh stands for P''_{q/2} *)
ff[x_]=f[x]/.{g[(((((1-x))^(2/p))+(x^(2/p)))/((((1-x))^(2/p))-(x^(2/p))))]->h};
dff[x_]=f'[x]/.{g[(((((1-x))^(2/p))+(x^(2/p)))/((((1-x))^(2/p))-(x^(2/p))))]->h,g'[(((((1-x))^(2/p))+(x^(2/p)))/((((1-x))^(2/p))-(x^(2/p))))]->dh};
ddff[x_]=f''[x]/.{g[(((((1-x))^(2/p))+(x^(2/p)))/((((1-x))^(2/p))-(x^(2/p))))]->h,g'[(((((1-x))^(2/p))+(x^(2/p)))/((((1-x))^(2/p))-(x^(2/p))))]->dh, g''[(((((1-x))^(2/p))+(x^(2/p)))/((((1-x))^(2/p))-(x^(2/p))))]->ddh};
ode[x_]=a[x] ddff[x]+b[x]dff[x] +c[x]ff[x];
phi[x_]=(((((1-x))^(2/p))+(x^(2/p)))/((((1-x))^(2/p))-(x^(2/p))));
(* Recall the Legendre ODE: (1-phi[x]^2)ddh-2phi[x]dh+(q/2)(q/2+1)h=0 *)
Simplify[ode[x]/.{ddh->(2phi[x]dh-(q/2)(q/2+1)h)/(1-phi[x]^2)}]

(* Illustrations of exponent regions *)

(* Figure 1, left *)
RegionPlot[x>=If[y<1/2,1-y,1/2],{x,0,1},{y,0,1},PlotStyle->Gray,BoundaryStyle->{Thick,Black},FrameTicks->{{0,1/2,1},{0,1/2,1}},FrameTicksStyle->Directive[FontSize->24,FontFamily->"Times",FontColor->Black],FrameStyle->Directive[Black,Thick]]

(* Figure 1, right *)
RegionPlot[x>=If[y<1/2,y Log[2,((Gamma[1/y+1])/(Gamma[1/(2y)+1]^2))],1/2],{x,0,1},{y,0,1},PlotStyle->Gray,BoundaryStyle->{Thick,Black},FrameTicks->{{0,1/2,1},{0,1/2,1}},FrameTicksStyle->Directive[FontSize->24,FontFamily->"Times",FontColor->Black],FrameStyle->Directive[Black,Thick]]

(* Figure 2, left *)
RegionPlot[2x>=y+1,{x,0,1},{y,0,1},PlotStyle->Gray,BoundaryStyle->{Thick,Black},FrameTicks->{{0,1/2,1},{0,1}},FrameTicksStyle->Directive[FontSize->24,FontFamily->"Times",FontColor->Black],FrameStyle->Directive[Black,Thick],ImageSize->Medium]

(* Figure 2, right *)
RegionPlot[x>=(y/2)Log[2,2^(1/y)+2],{x,0,1},{y,0,1},PlotStyle->Gray,BoundaryStyle->{Thick,Black},FrameTicks->{{0,1/2,1},{0,1}},FrameTicksStyle->Directive[FontSize->24,FontFamily->"Times",FontColor->Black],FrameStyle->Directive[Black,Thick],ImageSize->Medium]

(* Figure 3. Illustration of F_q from Section 4: Proof of Lemma 4 *)
q=4.;p=q/Log[2,Gamma[q+1]/Gamma[q/2+1]^2];
f[x_]=NIntegrate[((1-x)^(2/p)+x^(2/p)+2(1-x)^(1/p) x^(1/p) Cos[2Pi t])^(q/2),{t,0,1}];
Plot[f[x],{x,0,1},Ticks->{{0,1/4,1/2,3/4,1},{19/20,1}},TicksStyle->Directive[FontSize->14,FontFamily->"Times",FontColor->Black],PlotStyle->Directive[Black,Thick]]

(* Figure 5. Illustration of G_q from Section 9: Proof of Lemma 8 *)
q=2.;p=2q/Log[2,2^q+2];
g[x_,y_]=(1-x)^(q/p) (1-y)^(q/p)+((1-x)^(1/p) y^(1/p)+x^(1/p) (1-y)^(1/p))^q+x^(q/p) y^(q/p);
Plot3D[g[x,y],{x,0,1},{y,0,1},Ticks->{{0,1/2,1},{0,1/2,1},{9/10,1}},TicksStyle->Directive[FontSize->14,FontFamily->"Times",FontColor->Black],PlotStyle->Directive[Gray,Thick]]

(* Assistance with the proof of (9.5) of Lemma 12 *)

(* Verification of the formulae for p' and p'' *)
Clear[p,q,s,u];
p[q_]=2q/Log[2,2^q+2];

FullSimplify[p'[q]-(((p[q])/q)(1-((p[q] 2^q)/(2(2^q+2)))))]

FullSimplify[p'[1]]

(* Inequality (9.10) for p' *)
Print[NMaximize[{p'[q]/(((p[q])/q)),1<=q<=4},q,AccuracyGoal->20,PrecisionGoal->18,WorkingPrecision->40]];
Plot[p'[q]/(((p[q])/q)),{q,1,4}]

(* Inequality (9.10) for p'' *)
Print[NMinimize[{-p''[q]/(((p[q])/(q^2))),1<=q<=4},q,AccuracyGoal->20,PrecisionGoal->18,WorkingPrecision->40]];
Print[NMaximize[{-p''[q]/(((p[q])/(q^2))),1<=q<=4},q,AccuracyGoal->20,PrecisionGoal->18,WorkingPrecision->40]];
Plot[-p''[q]/(((p[q])/(q^2))),{q,1,4}]

(* Verification of (∂_q)phi[q,u]>0 when u>=3 *)
Print[NMinimize[{-((p[q] 2^q)/((2q-p[q])(2^q+2)))+Log[2],1<=q<=4},q,AccuracyGoal->20,PrecisionGoal->18,WorkingPrecision->40]];
Plot[-((p[q] 2^q)/((2q-p[q])(2^q+2)))+Log[2],{q,1,4}]

(* Verification of the second derivatives formulae for phi[q,u] *)
Clear[p,q,s,u];
phi[q_,u_]=q Log[Cosh[u]]-Log[Cosh[p[q] u]]+Log[((p[q] 2^(q-1))/(2q-p[q]))];

FullSimplify[(∂_q)((∂_q)phi[q,u])-(-p''[q] u Tanh [p[q] u]-p'[q]^2 ((u^2)/(Cosh[p[q] u]^2))+p''[q] ((2q)/(p[q](2q-p[q])))-p'[q]^2 ((4q(q-p[q]))/(p[q]^2(2q-p[q])^2))-p'[q] (4/((2q-p[q])^2))+(4/((2q-p[q])^2)))]

FullSimplify[(∂_u)((∂_q)phi[q,u])-(-p'[q]Tanh[p[q] u]-p'[q] ((p[q] u)/(Cosh[p[q] u]^2))+Tanh[u])]

(* A function (∂_q)phi[q,u] from the proof of Lemma 12 *)
Clear[p,q,s,u];
p[q_]=2q/Log[2,2^q+2];
phi[q_,u_]=q Log[Cosh[u]]-Log[Cosh[p[q] u]]+Log[((p[q] 2^(q-1))/(2q-p[q]))];
phid[q_,u_]=(∂_q)phi[q,u];
Plot3D[phid[q,u],{q,1,4},{u,0,4},Ticks->{{1,2,3,4},{0,1,2,3,4},Automatic},TicksStyle->Directive[FontSize->14,FontFamily->"Times",FontColor->Black],PlotStyle->Directive[Gray,Thick]]

NMinimize[{phid[q,u],1<=q<=4,0<=u<=4},{q,u},AccuracyGoal->20,PrecisionGoal->18,WorkingPrecision->40]

(* Verification of the single-number inequalities at the end of the proof of Lemma 12 *)
And@@Flatten[Table[phid[1+(i/700),(j/300)]>1/50,{i,1,2100},{j,0,900}]]

(* Assistance with (9.13) *)
Clear[p,q,x,y,g,h,a,b];
g[x_,y_]=(1-x)^(q/p) (1-y)^(q/p)+((1-x)^(1/p) y^(1/p)+x^(1/p) (1-y)^(1/p))^q+x^(q/p) y^(q/p);
h[x_]=g[x,((a x)/(1-x+a x))];
FullSimplify[h''[((1-(Sqrt[a]))/(1-a))],Assumptions->{a>0,a<1,p>1,q>1,q>p}]
(1/(a (p^2)))2 (((1+(Sqrt[a])))^(2-((2 q)/p))) q (-2 (a^((p+q)/(2 p))) p-(Sqrt[a]) (((1+(a^(1/p))))^q) p+2 (a^((p+q)/(2 p))) q+(a^(q/(2 p))) (1+a) q)

(* Verification of the PDE from Section 9, Remark 14 *)
Clear[p,q,x,y,g,a,b,c];
a[x_,y_]=p^2x^(1/p)(1-y)^(1/p);
b[x_,y_]=p (1-x) x (x^(1/p)(1-y)^(1/p) (p (1-2 x)+q (2x+2y-1))-(1-x)^(1/p) y^(1/p)(p(1-2 x)+q (2x+2y-3)));
c[x_,y_]=q x^(1/p)(1-y)^(1/p)(p ((1-x)x+(1-y)y)+q ((((x+y))^2)-x-3y));
g[x_,y_]=(1-x)^(q/p) (1-y)^(q/p)+((1-x)^(1/p) y^(1/p)+x^(1/p) (1-y)^(1/p))^q+x^(q/p) y^(q/p);
FullSimplify[(a[x,y]-a[y,x])((1-x)^2 x^2 Derivative[2,0][g][x,y]+(1-y)^2 y^2 Derivative[0,2][g][x,y]+2(1-x) x(1-y) yDerivative[1,1][g][x,y])+(b[x,y] Derivative[1,0][g][x,y]-b[y,x] Derivative[0,1][g][x,y])+(c[x,y]-c[y,x])g[x,y],Assumptions->{x>0,x<1,y>0,y<1,p>1,q>1,s>0,t>0}]