ClearAll["Global`*"];

Make sure that unzipped folder   qPI_orbit_length_data   is inside the folder containing this Mathematica file.

Data Analysis

Load  data

Load Magma data

Prime powers between 2 and 499

primepowers={2,3,4,5,7,8,9,11,13,16,17,19,23,25,27,29,31,32,37,41,43,47,49,53,59,61,64,67,71,73,79,81,83,89,97,101,103,107,109,113,121,125,127,128,131,137,139,149,151,157,163,167,169,173,179,181,191,193,197,199,211,223,227,229,233,239,241,243,251,256,257,263,269,271,277,281,283,289,293,307,311,313,317,331,337,343,347,349,353,359,361,367,373,379,383,389,397,401,409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499};

(*Set directory to where Magma data should be*)
SetDirectory[StringJoin[NotebookDirectory[],"qPI_orbit_length_data"]];
(*load raw data*)
DATAmagmaraw=Table[Import[StringJoin["qPI_data_q=",ToString[primepowers[[k]]],".txt"]],{k,1,Length[primepowers]}];
(*Turn data into nested Mathematica list*)
DATAmagma=Table[{primepowers[[k]],ToExpression[DATAmagmaraw[[k]]]},{k,1,Length[primepowers]}];
(*Set directory to directoy of notebook*)
SetDirectory[NotebookDirectory[]];

(*Count total number of reduced orbit length data points*)
total=0;
For[i1=1,i1<=Length[DATAmagma],i1++,
sublist1=Last[Part[DATAmagma,i1]];
For[i2=1,i2<=Length[sublist1],i2++,
sublist2=Last[Part[sublist1,i2]];
For[i3=1,i3<=Length[sublist2],i3++,
total = total + Last[Part[sublist2,i3]]]]]
Print["Total nr of reduced orbits lengths: ",total];
Clear[i1,i2,i3,sublist1,sublist2,total];

Example of entry

The fourth prime power is q=5.

DATAmagma[[4]];
First[%]

Here is the full list of data for q=5.

DATAmagma[[4]]

5-1=4 has three divisors, 1, 2 and 4. The data is organised by these divisors.

DATAmagma[[4]];
TableForm[Last[DATAmagma[[4]]],TableDepth->1]

Let’s focus on the entry for the divisor 4

DATAmagma[[4]];
Last[%][[3]]

This means that there are
4 orbits with reduced length 4,
2 orbits with reduced length 5,
2 orbits with reduced length 6,
2 orbits with reduced length 8,
2 orbits with reduced length 9,
of qPI over (Subscript[F, 5]) with s of multiplicative order 4.

Check Conjectures

Hasse upper bound conjecture

Conjecture 1.A



We check this statement against all orbit data.

(*Initialise counters that keep track of number of reduced orbit lengths that satisfy bound and those that do not.*)
ConTrue=0;
ConFalse=0;
(*Next, we iterate over all the data*)
(*The first iterate over prime powers q*)
For[qindex=1,qindex<=Length[DATAmagma],qindex++,
q=First[DATAmagma[[qindex]]];
(*Get corresponding data for this value of q*);
databyorders=Last[DATAmagma[[qindex]]];
(*Iterate over divisors r of q-1*)
For[rindex=1,rindex<=Length[databyorders],rindex++,
r=First[databyorders[[rindex]]];
(*Get reduced orbit lengths and frequencies for values of q and r*);
lengthsandfreqs=Last[databyorders[[rindex]]];
(*Finally, we iterate over all reduced orbit lengths*);
For[lengthindex=1,lengthindex<=Length[lengthsandfreqs],lengthindex++,
{length,freq}=lengthsandfreqs[[lengthindex]];
(*For each reduced orbit length, we check wether it satisfies the Hasse upper bound and add the frequency of this reduced orbit length to the ConTrue counter if the bound holds, and else we add it to do ConFalse counter*);
If[length<=q+2 Sqrt[q]+1,ConTrue =ConTrue+freq,ConFalse = ConFalse +freq]]]]
Print["Statement of conjecture is true for: ",ConTrue," orbits"]
Print["Statement of conjecture is false for:  ",ConFalse," orbits"]
(*Clear all the variables*)
Clear[ConTrue,ConFalse,qindex,rindex,q,r,databyorders,lengthsandfreqs,lengthindex,length,freq]

Check refined conjecture

Conjecture 1.B



We check this statement against all orbit data.

(*Initialise counters that keep track of number of reduced orbit lengths that satisfy statement and those that do not.*)
ConTrue=0;
ConFalse=0;
(*Next, we iterate over all the data*)
(*The first iterate over prime powers q*)
For[qindex=1,qindex<=Length[DATAmagma],qindex++,
q=First[DATAmagma[[qindex]]];
(*Compute integer Mq*)
Mq=Ceiling[(1/4)((Sqrt[q])+(1/(Sqrt[q]))-2)];
(*Definition of mth bin as list of integers between the two bounds*);
Bin[m_]:=If[m<Mq,Range[Ceiling[((q+1-2 (Sqrt[q]))/m)],Floor[((q+1+2 (Sqrt[q]))/m)]],Range[1,Floor[((q+1+2 (Sqrt[q]))/m)]]];
(*Get corresponding data for this value of q*);
databyorders=Last[DATAmagma[[qindex]]];
(*Iterate over divisors r of q-1*)
For[rindex=1,rindex<=Length[databyorders],rindex++,
r=First[databyorders[[rindex]]];
(*Get reduced orbit lengths and frequencies for values of q and r*);
lengthsandfreqs=Last[databyorders[[rindex]]];
(*Finally, we iterate over all reduced orbit lengths*);
For[lengthindex=1,lengthindex<=Length[lengthsandfreqs],lengthindex++,
{length,freq}=lengthsandfreqs[[lengthindex]];
(*For each bin, record whether it sits in there are not*);
withinbins=Table[MemberQ[Bin[m],length],{m,1,Mq}];
(*Next, we check whether it sits in at least one bin*)
inbin=withinbins/.List->Or;
(*Add frequency of this reduced orbit length to the ConTrue counter if reduced orbit length lies in a bin, else add it to do ConFalse counter*);
If[inbin,ConTrue =ConTrue+freq,ConFalse = ConFalse +freq]]]]
Print["Statement of conjecture is true for: ",ConTrue," orbits"]
Print["Statement of conjecture is false for:  ",ConFalse," orbits"]
Clear[ConTrue,ConFalse,qindex,rindex,q,Mq,databyorders,r,lengthsandfreqs,lengthindex,length,freq,withinbins,inbin]

Plotting reduced orbit lengths vs Hasse bound and bins

Plotting code

extractbyorder[data_,r_]:=Module[{output,index,subdata,q,allorbitlengthsreducedfreqlist,k,l,order,orbitlengthsreduced,orderdata},
output={};
For[index=1,index<=Length[data],index++,
subdata=data[[index]];
q=subdata[[1]];
allorbitlengthsreducedfreqlist=subdata[[2]];
For[k=1,k<=Length[allorbitlengthsreducedfreqlist],k++,
orderdata=allorbitlengthsreducedfreqlist[[k]];
order=orderdata[[1]];
If[order==r,orbitlengthsreduced=Table[First[Part[orderdata[[2]],l]],{l,1,Length[orderdata[[2]]]}];
AppendTo[output,{q,orbitlengthsreduced}]]]];
output]

{q1,q2,q3,q4,q5}=Table[Last[q/.Solve[((q+1)/(m+1))+2 ((Sqrt[q])/(m+1))==((q+1)/m)-2 ((Sqrt[q])/m),q]],{m,1,5}];
scatterplotdatavsHWboundadv[data_,r_,max_]:=Module[{orbitlengthsdata,scatterdata,primepowerindex,q,primepower,orbitlengths,lengthindex,HWplot,HWplot1u,HWplot1d,HWplot2u,HWplot2d,HWplot3u,HWplot3d,HWplot4u,HWplot4d,scatterplot},
orbitlengthsdata=extractbyorder[DATAmagma,r];
scatterdata={};
For[primepowerindex=1,primepowerindex<=Length[orbitlengthsdata],primepowerindex++,
q=First[Part[orbitlengthsdata,primepowerindex]];
orbitlengths=Last[Part[orbitlengthsdata,primepowerindex]];
For[lengthindex=1,lengthindex<=Length[orbitlengths],lengthindex++,
AppendTo[scatterdata,{q,Part[orbitlengths,lengthindex]}]]];
HWplot=Plot[qq+2Sqrt[qq]+1,{qq,0,max},PlotStyle->Red,AxesLabel->{Style["q",12],}];
HWplot1u=Plot[qq-2Sqrt[qq]+1,{qq,q1,max},PlotStyle->{Red}];
HWplot1d=Plot[((qq+2Sqrt[qq]+1)/2),{qq,q1,max},PlotStyle->{Red,Dashed}];
HWplot2u=Plot[((qq-2Sqrt[qq]+1)/2),{qq,q2,max},PlotStyle->{Red,Dashed}];
HWplot2d=Plot[((qq+2Sqrt[qq]+1)/3),{qq,q2,max},PlotStyle->{Red}];
HWplot3u=Plot[((qq-2Sqrt[qq]+1)/3),{qq,q3,max},PlotStyle->{Red}];
HWplot3d=Plot[((qq+2Sqrt[qq]+1)/4),{qq,q3,max},PlotStyle->{Red,Dashed}];
HWplot4u=Plot[((qq-2Sqrt[qq]+1)/4),{qq,q4,max},PlotStyle->{Red,Dashed}];
HWplot4d=Plot[((qq+2Sqrt[qq]+1)/5),{qq,q4,max},PlotStyle->{Red}];
scatterplot=ListPlot[scatterdata,PlotRange->{0,1.1 Last[Last[scatterdata]]},PlotMarkers->{"-",1}];
Show[HWplot,HWplot1u,HWplot1d,HWplot2u,HWplot2d,HWplot3u,HWplot3d,HWplot4u,scatterplot]]


scatterplotdatavsHWboundadvorigin[data_,r_,max_]:=Module[{orbitlengthsdata,scatterdata,primepowerindex,q,primepower,orbitlengths,lengthindex,HWplot,HWplot1u,HWplot1d,scatterplot},
orbitlengthsdata=extractbyorder[DATAmagma,r];
scatterdata={};
For[primepowerindex=1,primepowerindex<=Length[orbitlengthsdata],primepowerindex++,
q=First[Part[orbitlengthsdata,primepowerindex]];
orbitlengths=Last[Part[orbitlengthsdata,primepowerindex]];
For[lengthindex=1,lengthindex<=Length[orbitlengths],lengthindex++,
AppendTo[scatterdata,{q,Part[orbitlengths,lengthindex]}]]];
HWplot=Plot[qq+2Sqrt[qq]+1,{qq,0,max},PlotStyle->Red,AxesLabel->{Style["q",12],}];
HWplot1u=Plot[qq-2Sqrt[qq]+1,{qq,q1,max},PlotStyle->{Red}];
HWplot1d=Plot[((qq+2Sqrt[qq]+1)/2),{qq,q1,max},PlotStyle->{Red,Dashed}];
scatterplot=ListPlot[scatterdata,PlotRange->{0,1.1 Last[Last[scatterdata]]},PlotMarkers->{"-",8}];
Show[HWplot,HWplot1u,HWplot1d,scatterplot]]

Plots

Change rchoice to the value you are interested in

rchoice=5;
qmax=First[Last[DATAmagma]];
Print[scatterplotdatavsHWboundadv[DATAmagma,rchoice,qmax]]
Clear[qmax,rchoice]

Plots zoomed in near origin

Change rchoice to the value you are interested in

rchoice=4;
Print[scatterplotdatavsHWboundadvorigin[DATAmagma,rchoice,50]]
Clear[rchoice]