counting.nb

This notebook accompanies the paper M. Levitin, I. Polterovich, D. A. Sher, Pólya’s conjecture for the disk: a computer-assisted proof, available at https://arxiv.org/abs/2203.07696, see also https://michaellevitin.net/polya.html

For simplicity, this notebook differs slightly from the one used in the actual calculations described in the paper, by removing some non-essential embellishments. As a result, the output may also be slightly different from that available for download at https://michaellevitin.net/polya.html

This script reads files smallqminus.csv, smallqplus.csv, capitalqminus.csv, capitalqplus.csv, and capitalqr.csv from https://michaellevitin.net/polya.html

Notebook initialisation

SetDirectory[NotebookDirectory[]];
<<MaTeX`
texStyle={};
SetOptions[MaTeX,"BasePreamble"->{"\\usepackage{amsmath}","\\usepackage{xcolor}","\\usepackage{fourier}","\\usepackage{ebgaramond}"},FontSize->11];

Loading the coefficients of (Subscript[(UnderBar[h]), ℚ]) and (Subscript[(OverBar[h]), ℚ]) and converting them into fractions (compare with Table 4 in the paper)

qdataminus=Import["https://michaellevitin.net/PolyaData/smallqminus.csv"];
qdataplus=Import["https://michaellevitin.net/PolyaData/smallqplus.csv"];
Qdataminus=Import["https://michaellevitin.net/PolyaData/capitalqminus.csv"];
Qdataplus=Import["https://michaellevitin.net/PolyaData/capitalqplus.csv"];
QdataR=Import["https://michaellevitin.net/PolyaData/capitalqr.csv"];

In what follows qminus, qplus, Qminus, Qplus, and QR are the indexed quantities (SubMinus[q]),(SubPlus[q]),(SubMinus[Q]),(SubPlus[Q]) and (Subscript[Q, R]) from the paper, with the second index shifted by one, so that qminus[[i]][[j]]=(Subscript[q, -,i,j-1]), etc.

qminus=Table[qdataminus[[i,j]]/qdataminus[[i,j+1]],{i,14}, {j,1,11,2}];
qminus//TableForm
qplus=Table[qdataplus[[i,j]]/qdataplus[[i,j+1]],{i,14}, {j,1,11,2}];
qplus//TableForm

Qminus=Table[Qdataminus[[i,j]]/Qdataminus[[i,j+1]],{i,13}, {j,1,11,2}];
Qplus=Table[Qdataplus[[i,j]]/Qdataplus[[i,j+1]],{i,13}, {j,1,11,2}];
QR=Table[QdataR[[i,j]]/QdataR[[i,j+1]],{i,10}, {j,1,3,2}];
Qminus//TableForm
Qplus//TableForm
QR//TableForm

Creating partitions of [0,1] and lists of functions composing  (Subscript[(UnderBar[h]), ℚ]) and (Subscript[(OverBar[h]), ℚ])

xs=Join[Table[(i-1)/10, {i,1,7}], Table[6/10+(i-7)/20, {i, 8,15}]] (* points partitioning [0,1] into 14 original intervals (Subscript[E, i]) *)

xcs= Table[1/2(xs[[i]]+xs[[i+1]]),{i,14}] (* centres of the intervals, partitioning them further into subintervals (Subscript[E, ∓]) *)

xRs=Table[95/100+(k-1)/200,{k,1,11}] (* subpartition of the right-most interval (Subscript[E, 14])=[95/100,1]*)

xQs=Union[xs, xcs, xRs]//Sort (* all subinterval partition points for (Subscript[(OverBar[h]), ℚ]); 36 subintervals in total *)
xQs//Length

xqs=Union[xs, xcs]//Sort (* all subinterval partition points for (Subscript[(UnderBar[h]), ℚ]); 28 subintervals in total *)
xqs//Length

The following are the lists of expressions (also defined as functions later) for  (Subscript[(UnderBar[h]), ℚ]) and  (Subscript[(OverBar[h]), ℚ])  on subintervals. The list for (Subscript[(UnderBar[h]), ℚ]) is {(Subscript[(UnderBar[h]), ℚ,1,-]), (Subscript[(UnderBar[h]), ℚ,1,+]), (Subscript[(UnderBar[h]), ℚ,2,-]), (Subscript[(UnderBar[h]), ℚ,2,+]),...,(Subscript[(UnderBar[h]), ℚ,14,-]), (Subscript[(UnderBar[h]), ℚ,14,+])}, and the list for   (Subscript[(OverBar[h]), ℚ]) is  { (Subscript[(OverBar[h]), ℚ,1,-]), (Subscript[(OverBar[h]), ℚ,1,+]), ...,(Subscript[(OverBar[h]), ℚ,13,-]), (Subscript[(OverBar[h]), ℚ,13,+]), (Subscript[(OverBar[h]), ℚ,R,1]),...,(Subscript[(OverBar[h]), ℚ,R,10])}, where
(Subscript[(UnderBar[h]), ℚ,i,±])=(Subscript[q, i,0,±])+(((Subscript[q, i,j,±])(x-(Subscript[x, c,i])))^j),		(Subscript[(OverBar[h]), ℚ,i,±])=(Subscript[Q, i,0,±])+(((Subscript[Q, i,j,±])(x-(Subscript[x, c,i])))^j),		(Subscript[(OverBar[h]), ℚ,R,k])=(Subscript[Q, R,k,0])+(Subscript[Q, R,k,1])(x-(Subscript[x, R,k]))

hDown = Flatten[Table[{qminus[[i]][[1]]+Sum[qminus[[i]][[j]](x-xcs[[i]])^(j-1),{j,2,6}], qplus[[i]][[1]]+Sum[qplus[[i]][[j]](x-xcs[[i]])^(j-1),{j,2,6}]}, {i,1,14}]];
hDownfns=Table[Function[x,hDown[[i]]//Evaluate],{i,1,28}];
hUp = Join[Flatten[Table[{Qminus[[i]][[1]]+Sum[Qminus[[i]][[j]](x-xcs[[i]])^(j-1),{j,2,6}], Qplus[[i]][[1]]+Sum[Qplus[[i]][[j]](x-xcs[[i]])^(j-1),{j,2,6}]}, {i,1,13}]],
Table[QR[[k]][[1]]+QR[[k]][[2]](x-xRs[[k]]), {k,1,10}]
];
hUpfns=Table[Function[x,hUp[[i]]//Evaluate],{i,1,36}];

Checking the inequalities  (Subscript[(UnderBar[h]), ℚ])(x)<=h(x)<=(Subscript[(OverBar[h]), ℚ])(x)

This section is for illustrative purposes only. We demonstrate (numerically) the validity of the inequalities above, cf. Figure 9 in the paper. We omit some embellishments made  in the paper in order to improve the graphics presentation.

h[x_]:=1/Pi(Sqrt[1-x^2]-x ArcCos[x]);
PiecewisehDown=Append[Table[{hDown[[i]], xqs[[i]]<=x<xqs[[i+1]]},{i,1,27}],{hDown[[28]], xqs[[28]]<=x<=xqs[[29]]}];
PiecewisehUp=Append[Table[{hUp[[i]], xQs[[i]]<=x<xQs[[i+1]]},{i,1,35}],{hUp[[36]], xQs[[36]]<=x<=xQs[[37]]}];
xticks1={{0.3,0.6,0.9},MaTeX@{"0.3","0.6","0.9"}}//Transpose;
xticks2={{0.9,0.95,1},MaTeX@{"0.9","0.95","1"}}//Transpose;
yticks1={{1,2,3} 10^(-7),MaTeX@{"10^{-7}","2\\cdot 10^{-7}","3\\cdot 10^{-7}"}}//Transpose;
yticks2={{3,6,9} 10^(-6),MaTeX@{"3\\cdot 10^{-6}","6\\cdot 10^{-6}","9\\cdot 10^{-6}"}}//Transpose;
yticks3={{1,2} 10^(-6),MaTeX@{"10^{-6}","2\\cdot 10^{-6}"}}//Transpose;
xlabel=MaTeX["x"];
ylabel1=MaTeX["\\overline{h}_{\\mathbb{Q}}(x)-h(x)"];
ylabel2=MaTeX["\\underline{h}_{\\mathbb{Q}}(x)-h(x)"];
GraphicsGrid[{{Plot[Piecewise[PiecewisehUp]-h[x],{x,0,0.9},PlotRange->Automatic, PlotStyle->Blue,AxesOrigin->{0,0}, Ticks->{xticks1,yticks1},AxesLabel->{xlabel,ylabel1}],
Plot[Piecewise[PiecewisehUp]-h[x],{x,0.9,1},PlotRange->Automatic, PlotStyle->Blue,AxesOrigin->{0.895,0},Ticks->{xticks2,yticks2},AxesLabel->{xlabel,ylabel1}]},
{Plot[h[x]-Piecewise[PiecewisehDown],{x,0,0.9},PlotRange->Automatic,PlotStyle->Red,AxesOrigin->{0,0},Ticks->{xticks1,yticks1}, AxesLabel->{xlabel,ylabel2}],Plot[h[x]-Piecewise[PiecewisehDown],{x,0.9,1},PlotRange->Automatic,PlotStyle->Red,AxesOrigin->{0.895,0},Ticks->{xticks2,yticks3}, AxesLabel->{xlabel,ylabel2}]}}]

Counting the points and setting up the iteration processes

The next two functions compute ((OverBar[𝒫])^D)(λ) and ((UnderBar[𝒫])^N)(λ) directly

PCountD[λ_]:=Module[{m,i, mP,P, Nf},
Nf=36;
mP=Floor[xQs λ]; mP[[1]]=0;
P= Sum[Floor[λ hUpfns[[i]][m/λ]+1/4],{i,1,Nf},{m,mP[[i]]+1,mP[[i+1]]}];
P=2P+Floor[λ hUpfns[[1]][0]+1/4] 
]
PCountN[λ_]:=Module[{m,i, mP,P, Nf},
Nf=28;
mP=Floor[xqs λ]; mP[[1]]=0;
P= Sum[Floor[λ hDownfns[[i]][m/λ]+3/4],{i,1,Nf},{m,mP[[i]]+1,mP[[i+1]]}];
2P+Floor[λ hDownfns[[1]][0]+3/4] 
]

Example:

Λ1=84123;
pd1=PCountD[Λ1]
Λ1^2/4-pd1
pn1=PCountN[Λ1]
pn1-Λ1^2/4

The following auxiliary function Fori[α0, d] realises Forišek’s algorithm (see the citation in the paper) of finding the smallest possible denominator rational approximation to a rational number α0>0 within the interval (α0-d, α0+d). Its typical uses are Fori[α0+ϵ, δ] with small rational values 0<δ<ϵ to give a rational approximation of α0 from above with a small denominator, and Fori[α0-ϵ, δ] to give a rational approximation of α0 from below with a small denominator. Using Fori is not, strictly speaking, necessary, but ensures that the denominators of the fractions appearing in the iteration processes are kept relatively small, at the cost of slightly slowing the processes.

Fori[α0_,d_]:=Module[{α, denom,αnum,dnum,xnum, xdenom,x, pa, qa, pb,qb,aa,bb,cc, k,knum,kdenom, nsteps,newpa,newpb,newqa,newqb},
α=α0-Floor[α0];
If[α==0, Return[α0]];
denom=LCM[Denominator[α],Denominator[d] ];
αnum=α denom; dnum=d denom;
pa=0;qa=1;pb=1;qb=1; nsteps=0;
While[nsteps<=1000,
xnum=denom pb-αnum qb;xdenom=-denom pa+αnum qa;x=Ceiling[xnum/xdenom];
aa=(α-d<(pb+x pa)/(qb+x qa)<α+d);
bb=(α-d<(pb+(x-1) pa)/(qb+(x-1) qa)<α+d);
If[(aa||bb), 
knum=denom pb-(αnum+dnum)qb; 
kdenom=(αnum+dnum)qa-denom pa; 
k=Ceiling[knum/kdenom];
cc=(pb+k pa)/(qb+k qa);Break[]
];
newpa=pb+(x-1) pa;
newqa=qb+(x-1)qa;
newpb=pb+x pa;
newqb=qb+x qa;
pa=newpa;
pb=newpb;
qa=newqa;
qb=newqb;
(* *************** *)
xnum=-denom pb+αnum qb;
xdenom=denom pa-αnum qa;
x=Ceiling[xnum/xdenom];
aa=(α-d<(pb+x pa)/(qb+x qa)<α+d);
bb=(α-d<(pb+(x-1) pa)/(qb+(x-1) qa)<α+d);
If[(aa||bb), 
knum=-denom pb+(αnum-dnum)qb; 
kdenom=-(αnum-dnum)qa+denom pa; 
k=Ceiling[knum/kdenom];
cc=(pb+k pa)/(qb+k qa);Break[]
];
newpa=pb+(x-1) pa;
newqa=qb+(x-1)qa;
newpb=pb+x pa;
newqb=qb+x qa;
pa=newpa;
pb=newpb;
qa=newqa;
qb=newqb;
nsteps++;
];
cc+Floor[α0]
];

Example of the use of Fori:

α0=(14153214853/168246)
αAbove=Fori[α0+1/100,1/101]
αBelow=Fori[α0-1/100,1/101]
αBelow<α0<αAbove

The following two functions IterD and IterN each perform one step of the iteration process for the Dirichlet and Neumann point counts, respectively, at a point λ, see Figure 7. The output in the Dirichlet case is {λ,  ((OverBar[𝒫])^D)(λ), (e^D)(λ), (shift^D)(λ), (Subscript[λ, new])}, where

(e^D)(λ)=(λ^2)/4- ((OverBar[𝒫])^D)(λ),		(shift^D)(λ)=2(e^D)(λ)/λ,		(Subscript[λ, new])= Fori[λ-(shift^D)(λ)+1/100,1/101].

The output in the Neumann case is  {λ,  ((UnderBar[𝒫])^N)(λ), (e^N)(λ), (shift^N)(λ), (Subscript[λ, new])}, where

(e^N)(λ)= ((UnderBar[𝒫])^N)(λ)-(λ^2)/4,		(shift^N)(λ)=2(e^N)(λ)/(λ+2(e^N)(λ)/λ),	(Subscript[λ, new])= Fori[λ+(shift^N)(λ)-1/100,1/101].

IterD[λ_]:=Module[{err,λnew, shift,pD},
pD=PCountD[λ];
err=λ^2/4-pD;
shift=2err/λ;
λnew=Fori[λ-shift+1/100,1/101];
{λ, pD, err, shift, λnew}
];
IterN[λ_]:=Module[{err,λnew, shift,pN},
pN=PCountN[λ];
err=pN-λ^2/4;
shift=2err/(λ+2err/λ);
λnew=Fori[λ+shift-1/100,1/101];
{λ, pN, err, shift, λnew}
];

Example:

IterD[Λ1]

IterN[5/2]

Actual calculations

The following functions (shown here without commands used for output into external files) have been used to perform the actual calculations. The parameter freq in both functions stands for the frequency of screen output, thus taking freq=1 will output on each step, and freq=100 on every hundredth step.  The parameter id is a string used to identify the output from a particular calculation (useful in parallelised processes). For ease of observation only, some output is converted into floating point numbers.

DoIterationsD[ΛStart_, ΛMin_,id_,freq_]:=Module[{nI,nIfinal, times0,err, Λnew,Λ,pD,shift, time},
times0=TimeUsed[];
Λ=ΛStart;
Do[
{Λ, pD, err, shift, Λnew}=IterD[Λ]; 
time=TimeUsed[]-times0;
If[Mod[nI,freq]==0,Print[id<>": ",nI, " iterations    ", "      λ=",N[Λnew,7], "      time=", time]]; 
Λ=Λnew;
If[shift<=0,Return[id <> ": FAILS at λ="<>ToString[Λ]]];
If[Λ<ΛMin, nIfinal=nI; Break[]];
,
{nI,1,∞}
];
id<>": FINISHED!!!"<>" in time "<>ToString[time]<>" sec  after "<>ToString[nIfinal]<>"  iterations"
];
DoIterationsN[ΛStart_, ΛMax_,id_,freq_]:=Module[{nI,nIfinal,times0,err, Λnew,Λ,pN,shift, time},
times0=TimeUsed[];
Λ=ΛStart;
Do[
{Λ, pN, err, shift, Λnew}=IterN[Λ]; 
time=TimeUsed[]-times0;
If[Mod[nI,freq]==0,Print[id<>": ",nI, " iterations    ", "      λ=",N[Λnew,7], "      time=", time]]; 
Λ=Λnew;
If[shift<=0,Return[id <> ": FAILS at λ="<>ToString[Λ]]];
If[Λ>ΛMax, nIfinal=nI; Break[]];
,
{nI,1,∞}
];
id<>": FINISHED!!!"<>" in time "<>ToString[time]<>" sec  after "<>ToString[nIfinal]<>"  iterations"
];

Examples:

DoIterationsD[84123, 84100,"Dirichlet Test", 5]

DoIterationsN[84100, 84123,"Neumann Test", 5]

The full computer-assisted proof can be verified, for example, by executing the following command (WARNING: may take a very long time, especially on slower computers; ignore the indication of remaining time)

Λ1=84123; Λ0=5/2;
Parallelize[{
DoIterationsD[80000, 70000,"Dirichlet 80K-70K",500],
DoIterationsD[70000, 60000,"Dirichlet 70K-60K",500],
DoIterationsD[60000, 50000,"Dirichlet 60K-50K",500],
DoIterationsD[50000, 40000,"Dirichlet 50K-40K",500],
DoIterationsD[40000, 30000,"Dirichlet 40K-30K",500],
DoIterationsD[30000, 20000,"Dirichlet 30K-20K",500],
DoIterationsD[20000, 10000,"Dirichlet 20K-10K",500],
DoIterationsD[10000, Λ0,"Dirichlet 10K-2",500],
DoIterationsN[Λ0, 10000,"Neumann 2-10K",500],
DoIterationsN[10000, 20000,"Neumann 10K-20K",500],
DoIterationsN[20000, 30000,"Neumann 20K-30K",500],
DoIterationsN[30000, 40000,"Neumann 30K-40K",500],
DoIterationsN[40000, 50000,"Neumann 40K-50K",500],
DoIterationsN[50000, 60000,"Neumann 50K-60K",500],
DoIterationsN[60000, 70000,"Neumann 60K-70K",500],
DoIterationsN[70000, 80000,"Neumann 70K-80K",500],
DoIterationsN[80000, Λ1,"Neumann 80K-84K",500]
}]