(* This code computes a lower bound on the quantity M_k,epsilon that has been the focus of the Polymath8b project.
The parameters and computation shown are those which achieve the key threshold M=4 for the smallest k we have found to work, and with the most efficient computation to do so that we have discovered. *)

(* This code can only explore the k-dimensional simplex defined as the region of the positive orthant in which the sum of all coordinates does not exceed (1+epsilon). Other work, exploring larger regions but only for k=3, uses different code posted elsewhere. *)

(* The following global variables, together with the function "Needed" and the argument "eps" of the main function "Output", define the computation. kg is the dimensionality of the polytope being explored, eps is the value of epsilon, and dg is the maximum degree of polynomials to be explored. Mgoal is the value of M_k,epsilon that we hope to achieve, and Mrange is the range of values around Mgoal within which we are confident M will lie and which we will iteratively bisect, in the final step, until we have narrowed it to a range specified by acc. Tset is the subset of [2,dg] allowed in signatures; most of our computations have been done restricting this to even numbers, which (for undiscovered reasons) is particularly efficient in reducing the complexity of the calculations without greatly reducing the eventual lower bound on M. *)

kg=50;
dg=25;
Mgoal=4;
Mrange=1/200;
acc=1/1000000;
Tset=Sort[Join[Table[2*j,{j,1,dg}],{}]];

(* The function Needed further constrains the terms to be included in the eventual polynomial. Needed can be whatever you like, including "True". Its argument is a pair whose first term is the signature and second is the power of (1-P_1). The condition shown arises because it turns out that signatures whose power of (1-P_1) is even can be eliminated at very little cost to the eventual M we derive, and also the smaller a signature’s first (and thus largest) element is the smaller the cost of
omitting terms using it (the signature {} being cheapest of all). More elaborate criteria could probably shrink the calculation quite a bit further while still keeping M above 4. *)

Needed[term_]:=(OddQ[term[[2]]])||((term[[1]]≠{})&&(term[[1,1]]≥12));

(* PossibleDegrees computes all possible signatures up to total degree dstart, sorted by their total. *)

PossibleDegrees[dstart_]:=SortBy[Flatten[Union[Table[IntegerPartitions[ii,All,Intersection[Table[j,{j,2,dg}],Tset]],{ii,0,dstart}]],1],Total];

InitialDegrees=PossibleDegrees[dg];Print["ID ",Length[InitialDegrees]," done at ",DateList[]," ",InitialDegrees];

(* CountRepeats counts the total number of variables in m(a). *)

CountRepeats[k0_,SetInput_]:=Apply[Multinomial,Join[{k0-Total[Tally[SetInput][[All,2]]]},Tally[SetInput][[All,2]]]];

(* The difference of two lists is much more laborious to compute than their combination, so we precompute it here for all our pairs of signatures. *)

IDdiffs=Module[{a,b,ida,idb,c,idd},
idd=Table[0,{i,1,Length[InitialDegrees]},{j,1,Length[InitialDegrees]}];
For[a=1,a≤Length[InitialDegrees],a++,
ida=InitialDegrees[[a]];
For[b=1,b≤a,b++,
idb=InitialDegrees[[b]];
If[MemberQ[InitialDegrees,Reverse[Sort[Join[ida,idb]]]],
c=First[Position[InitialDegrees,Reverse[Sort[Join[ida,idb]]]]];
idd[[c,a]]=b;idd[[c,b]]=a
]
]
];
idd
];

(* It turns out to save a lot of time to precompute these really simple things *)

faclist=Table[Factorial[n],{n,0,kg+100}];
binlist=Table[Table[1/(m*Binomial[n,m]),{m,1,n}],{n,1,50}];
subslist=Table[Union[Subsets[InitialDegrees[[a]],{i}]],{a,1,Length[InitialDegrees]},{i,0,dg}];
permslist=Table[Apply[Union,Table[Permutations[subslist[[j,k,i]]],{i,1,Length[subslist[[j,k]]]}]],{j,1,Length[InitialDegrees]},{k,1,dg}];
pdl=Table[Length[PossibleDegrees[d]],{d,0,dg+2}];
pdls=Table[Sum[pdl[[i+1]],{i,j,k}],{j,0,dg+2},{k,0,dg+2}];

(* WARNING: the required sizes of the following arrays have been determined empirically, Beware of overflow! *)

SSvals={};
Table[{},{i,1,Max[400,Length[InitialDegrees]]},{j,1,Max[400,Length[InitialDegrees]]}];
SSsigvals=%;
SigReps=%;

zmax=300;
TI3vals=Table[777,{z,1,zmax},{n,1,Length[SSvals]}];
TI1vals=Table[{},{z,1,zmax},{i,1,Max[20,dg+2]},{j,1,Max[40,2*dg+2]}];
TI1valsb=TI1vals;
TI2vals=TI1vals;

(* Here we precompute the expressions for which we will later compute the integral, for each pair of signatures to be considered in the main computation. We need to compute these without reference to Needed, because signatures that Needed entirely excludes may still be required for some calculations. *)

NewDegreeAdder[aindex_,bindex_]:=
Module[{counter,OutSets,perms,bsubs,T,newtuple,newstep,j,filla,fillb,count,NewPosition,tt,OSc1,T1,T2,TG,oi1,oi2,fp},
If[(Mod[aindex,10]==0)&&(bindex==1),Print["NDA ",aindex," ",bindex," ",DateList[]]];
OutSets={};
For[counter=0,counter≤Min[Length[InitialDegrees[[aindex]]],Length[InitialDegrees[[bindex]]]],counter++,
perms=permslist[[aindex,counter+1]];
bsubs=subslist[[bindex,counter+1]];
T=Union[Flatten[Table[Sort[Table[{perms[[j,i]],bsubs[[k,i]]},{i,1,counter}]],{j,1,Length[perms]},{k,1,Length[bsubs]}],1]];
For[ti=1,ti≤Length[T],ti++,
T1=First[First[Position[InitialDegrees,Reverse[Sort[T[[ti,All,1]]]]]]];
T2=First[First[Position[InitialDegrees,Reverse[Sort[T[[ti,All,2]]]]]]];
filla=InitialDegrees[[IDdiffs[[aindex,T1]]]];
fillb=InitialDegrees[[IDdiffs[[bindex,T2]]]];
newtuple=Reverse[Sort[Join[Table[T[[ti,i,1]]+T[[ti,i,2]],{i,1,counter}],filla,fillb]]];
TG=GatherBy[Join[T[[ti,1;;counter]],Table[{filla[[i]],0},{i,1,Length[filla]}],Table[{0,fillb[[i]]},{i,1,Length[fillb]}]],First[#]+Last[#]&];
newstep=Product[Apply[Multinomial,Tally[TG[[i]]][[All,2]]],{i,1,Length[TG]}]; (* seems sums are reliably members of U[newtuple] *)
NewPosition=Position[OutSets[[All,1]],newtuple];
If[NewPosition=={}||NewPosition=={{}},
OutSets=Append[OutSets,{newtuple,newstep}],
OutSets=ReplacePart[OutSets,{NewPosition[[1,1]],2}->OutSets[[NewPosition[[1,1]],2]]+newstep]
]
]
];
For[count=1,count≤Length[OutSets],count++,
OSc1=OutSets[[count,1]];oi1=Length[OSc1]+1;oi2=Total[OSc1]+1;
If[MemberQ[SSsigvals[[oi1,oi2]],OSc1],
OutSets[[count,1]]={oi1,oi2,First[First[Position[SSsigvals[[oi1,oi2]],OSc1]]]},
SSsigvals[[oi1,oi2]]=Append[SSsigvals[[oi1,oi2]],OSc1];
(* We make the TI arrays ragged in their 4th dimension for memory reasons. The first dimension, the power of (1-P_1), could also be made ragged but that gains little and complicates the code. *)
For[counter=1,counter≤zmax,counter++,
TI1vals[[counter,oi1,oi2]]=Append[TI1vals[[counter,oi1,oi2]],777];
TI1valsb[[counter,oi1,oi2]]=Append[TI1valsb[[counter,oi1,oi2]],777];
TI2vals[[counter,oi1,oi2]]=Append[TI2vals[[counter,oi1,oi2]],777]
];
fp=Product[faclist[[1+OSc1[[i]]]],{i,1,oi1-1}];
SigReps[[oi1,oi2]]=Append[SigReps[[oi1,oi2]],{fp*CountRepeats[kg,OSc1],fp*CountRepeats[kg-1,OSc1]}];
OutSets[[count,1]]={oi1,oi2,Length[SSsigvals[[oi1,oi2]]]};
]
];
SSvals=Append[SSvals,OutSets]; (* Originally this was saved for lookup and reuse, but it seems they never repeat anyway *)
Length[SSvals]
];

SuperSumSet=Table[Table[NewDegreeAdder[p,q],{q,1,p}],{p,1,Length[InitialDegrees]}];Print["SSS done at ",DateList[]];

(* The TermIntegral function tells us the integral over eps*R_{k0} of (eps-P1)^z m_{T}. Then,the TI functions store and retrieve TermIntegral computations and their combinations, greatly reducing execution time. *)

TermIntegral[k0_,T2_,z_,eps_]:=Power[eps,z+T2[[2]]+k0-1]*faclist[[1+z]]*SigReps[[T2[[1]],T2[[2]],T2[[3]],kg-k0+1]]/faclist[[z+T2[[2]]+k0]];

TI1[T_,z_,eps_]:=Module[{i},
i=TI1vals[[z+1,T[[1]],T[[2]],T[[3]]]];
If[i==777,i=TermIntegral[kg-1,T,z,1-eps];TI1vals[[z+1,T[[1]],T[[2]],T[[3]]]]=i];
i
];

TI1b[T_,z_,eps_]:=Module[{i},
i=TI1valsb[[z+1,T[[1]],T[[2]],T[[3]]]];
If[i==777,i=TermIntegral[kg,T,z,1+eps];TI1valsb[[z+1,T[[1]],T[[2]],T[[3]]]]=i];
i
];

(* TermIntegral2 tells us the integral over (1-eps)*R_{k0} of (1+eps-P1)^z m_{T}, for signature T. *)

TermIntegral2[T_,z_,eps_]:=Module[{a,b,t,i},
a=1;b=1;t=0;
For[i=0,i≤z,i++,
t=t+(a*b*TI1[T,z-i,eps]);
a=a*2*eps;b=b*(z-i)/(i+1)
];
t
];

TI2[T_,z_,eps_]:=Module[{i},
i=TI2vals[[z+1,T[[1]],T[[2]],T[[3]]]];
If[i==777,i=TermIntegral2[T,z,eps];TI2vals[[z+1,T[[1]],T[[2]],T[[3]]]]=i];
i
];

TermIntegral3[SS_,z_,eps_]:=Sum[SS[[m,2]]*TI2[SS[[m,1]],z,eps],{m,1,Length[SS]}];

TI3[T_,z_,eps_]:=Module[{i},
i=TI3vals[[z+1,T]];
If[i==777,i=TermIntegral3[SSvals[[T]],z,eps];TI3vals[[z+1,T]]=i];
i
];

(* The zdegFinder function helps us figure out what degree z and signature alpha we consider next for (1-P1)^z * m_alpha, in our loop in the Output function below. *)

zdegFinder[iplug_]:=Module[{startmeasure,newz,outplace,outzdeg,yy,d0},
d0=dg;
startmeasure=pdls[[1,dg+1]];
If[iplug>startmeasure,
newz=Ceiling[(iplug-startmeasure)/pdl[[dg+1]]];
outplace=iplug-startmeasure-(newz-1)*pdl[[dg+1]];
outzdeg=newz+(d0-Total[PossibleDegrees[d0][[outplace]]]),
yy=0;While[iplug>pdls[[dg-yy+1,dg+1]],yy=yy+1];
outzdeg=yy;outplace=iplug-pdls[[d0-yy+2,d0+1]]
];
{outzdeg+1,outplace}
]

(* Children lists the signatures resulting from removing one element from each of the signatures being explored. *)

Children[ID_]:=
{
cl2={};
For[idc=1,idc≤Length[ID],idc++,
clist={};
tally1=Tally[ID[[idc]]];
For[tc=1,tc≤Length[tally1],tc++,
zinta=ReplacePart[tally1,{tc,2}->Extract[tally1,{tc,2}]-1];
zintaprime={};
For[i=1,i≤Length[zinta],i++,
zintaprime=Join[zintaprime,Table[zinta[[i,1]],{m,1,zinta[[i,2]]}]];
];
clist=Append[clist,First[First[Position[ID,zintaprime]]]]
];
cl2=Append[cl2,clist]
];
cl2
}

(* Now the main routine. We build matrices M1 and M2 as in James's paper, and combine them to derive a lower bound on M. *)

Output[eps_]:=Module[{},
IDkids=First[Children[InitialDegrees]];
TI3vals=Table[777,{z,1,zmax},{n,1,Length[SSvals]}];
(* We store various computations for reuse in 4D arrays that need to be ragged to keep size under control. The needed length of each list in the fourth dimension has already been determined during the construction of SuperSumSet. We reset the values here so that Output can be executed multiple times (within a Table, for instnnce) in the same run, with different epsilon. *)
 For[z=1,z≤zmax,z++,
For[i=1,i<=Max[20,dg+2],i++,
 For[j=1,j<=Max[40,2*dg+2],j++,
For[k=1,k≤Length[TI1vals[[z,i,j]]],k++,
TI1vals[[z,i,j,k]]=777;
TI1valsb[[z,i,j,k]]=777;
TI2vals[[z,i,j,k]]=777
]
]
]
];
(* varsnum is the total number of basis elements we POTENTIALLY consider, i.e. BEFORE applying the constraint specified by Needed. *)
varsnum=Sum[Length[PossibleDegrees[m]],{m,0,dg}];
M1=Table[Table[0,{pp,1,varsnum}],{qq,1,varsnum}];M2=M1;
Mloss=Table[{{},0},{i,1,varsnum}];
(* We have now constructed empty matrices of the appropriate size. We will fill in each entry successively. *)
Print["starting with ",varsnum," terms at ",DateList[]];
For[ii=1,ii≤varsnum,ii++,
If[Mod[ii,10]==0,Print[ii," ",DateList[]]];
{zadeg,nicei}=zdegFinder[ii];
Mloss[[ii]]={InitialDegrees[[nicei]],zadeg-1};
If[Needed[Mloss[[ii]]],
For[jj=1,jj≤ii,jj++,
{zbdeg,nicej}=zdegFinder[jj];
If[Needed[Mloss[[jj]]],
SumSetVal=SSvals[[SuperSumSet[[Max[nicei,nicej],Min[nicei,nicej]]]]];
M1[[ii,jj]]=Sum[SumSetVal[[m,2]]*TI1b[SumSetVal[[m,1]],zadeg+zbdeg-2,eps],{m,1,Length[SumSetVal]}];
(* We have performed the squaring of F, and then put the appropriate entries into the I integral. Our next job is to compute the J entry. To do this, given a function m(a) we have to run over all possible integrations before we do the initial squaring. We will run over the possibilities using "tally1" and "tally2". *)
tally1=Tally[InitialDegrees[[nicei]]];tally2=Tally[InitialDegrees[[nicej]]];
M2cell=M2[[ii,jj]];
For[tally1count=0,tally1count≤Length[tally1],tally1count++,
For[tally2count=0,tally2count≤Length[tally2],tally2count++,
(* There are now four possibilities, according to whether or not tally1=0 or tally2=0. This corresponds to whether or not our initial integration (before squaring) messed with the signature or not. The formulae below follow James's formula (7.5). *)
If[tally1count==0,
newi=nicei,
t11=tally1[[tally1count,1]];
newi=IDkids[[nicei,tally1count]]
];
If[tally2count==0,
newj=nicej,
t21=tally2[[tally2count,1]];
newj=IDkids[[nicej,tally2count]]
];
SumSet=SuperSumSet[[Max[newi,newj],Min[newi,newj]]];
If[tally1count==0,
If[tally2count==0,
M2cell=M2cell+TI3[SumSet,zadeg+zbdeg,eps]/(zadeg*zbdeg),
M2cell=M2cell+binlist[[zbdeg+t21,zbdeg]]*TI3[SumSet,zadeg+zbdeg+t21,eps]/zadeg
],
If[tally2count==0,
M2cell=M2cell+binlist[[zadeg+t11,zadeg]]*TI3[SumSet,zadeg+zbdeg+t11,eps]/zbdeg,
M2cell=M2cell+binlist[[zadeg+t11,zadeg]]*binlist[[zbdeg+t21,zbdeg]]*TI3[SumSet,zadeg+zbdeg+t11+t21,eps]
]
]
]
];
M2[[ii,jj]]=M2cell;
If[jj<ii,M1[[jj,ii]]=M1[[ii,jj]];M2[[jj,ii]]=M2[[ii,jj]]]
]
]
]
];
Print["Matrices filled at ",DateList[]];

(* We now apply Needed, removing rows/columns from M1 and M2 that specify [signature,z] pairs which we guess minorly impact the eventual M. Actually we already applied it in the main routine, by not calculating the cells we're now discarding. *)

M1r={};
(* We could shrink M1 and M2 more succinctly using Drop or Join, but that turns out to run far more slowly *)
For[v=1,v≤varsnum,v++,If[Needed[Mloss[[v]]],M1r=Append[M1r,M1[[v]]]]];
M1r2=Transpose[M1r];M1r={};
For[v=1,v≤varsnum,v++,If[Needed[Mloss[[v]]],M1r=Append[M1r,M1r2[[v]]]]];
M1r=Transpose[M1r];
M2r={};
For[v=1,v≤varsnum,v++,If[Needed[Mloss[[v]]],M2r=Append[M2r,M2[[v]]]]];
M2r2=Transpose[M2r];M2r={};
For[v=1,v≤varsnum,v++,If[Needed[Mloss[[v]]],M2r=Append[M2r,M2r2[[v]]]]];
M2r=Transpose[M2r];
Print["reduced matrix size to ",Length[M1r]];

(* We are now ready to lower-bound M by iterative positive definite tests. M was previously determined in one step by calculating a generalised eigenvector, as shown below, but for large d that is computationally too intensive. *)

(* EigVec=Eigenvectors[{kg*N[M2,100],N[M1,100]},1][[1]];
RatVec=Rationalize[EigVec,10^(-40)];
Ratio=kg*(RatVec.M2.RatVec)/(RatVec.M1.RatVec);
Print["M for ",eps," is at least ",N[Ratio,20]]; *)

If[PositiveDefiniteMatrixQ[Mgoal/kg*N[M1r,40]-N[M2r,40]],
Print["bah! ",DateList[]];pdvals={{Mgoal-Mrange,False},{0,True},{Mgoal,True}},
Print["victory! ",DateList[]];pdvals={{Mgoal,False},{0,True},{Mgoal+Mrange,True}}
];
While[pdvals[[3,1]]-pdvals[[1,1]]>acc,
pdvals[[2,1]]=(pdvals[[3,1]]+pdvals[[1,1]])/2;
pdvals[[2,2]]=PositiveDefiniteMatrixQ[pdvals[[2,1]]/kg*N[M1r,40]-N[M2r,40]];
If[pdvals[[2,2]],pdvals[[3]]=pdvals[[2]],pdvals[[1]]=pdvals[[2]]];
Print["M is between ",N[pdvals[[1,1]],20]," and ",N[pdvals[[3,1]],20]," ",DateList[]]
];
Print["M for ",eps," is at least ",N[pdvals[[1,1]],20]];
N[pdvals[[1,1]],20]
];

(* OK, let's go! To optimise epsilon for a given dimensionality and maximum degree you can call Output within a Table here. Other parameters currently implemented as global variables can of course also be moved into parameters of Output but must then be carried through to the functions it calls. *)

Output[1/27]