(* These global variables define the computation. kg is the dimension being explored; eg is the epsilon value; dg is the maximum degree of polynomials to be explored. Mgoal is the value of M_kg,eg that we hope to achieve; 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. Tset is the subset of [2,dg] that we allow in signatures; most computations have been done restricting to even numbers, which is particularly efficient in reducing the complexity of the calculations without reducing the eventual M much. *)

kg=50;
dg=12;
eg=1/26;
Mgoal=4;
Mrange=1/3;
Tset=Sort[Join[Table[2*j,{j,1,dg}],{11}]];

(* PossibleDegrees finds all possible signatures up to total degree d0. *)

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

(* The CountRepeats function just 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]]]];

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

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

subslist=Table[Union[Subsets[InitialDegrees[[a]],{counter}]],{a,1,Length[InitialDegrees]},{counter,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}];

faclist=Table[Factorial[n],{n,0,kg+100}];

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

Table[777,{z,1,300},{i,1,dg+2},{j,1,2*dg+2},{m,1,Length[InitialDegrees]}];
TI1vals=%;
TI1valsb=%;
TI2vals=%;

NewDegreeAdder[aindex_,bindex_]:=
Module[{counter,OutSets,asubs,perms,bsubs,newperms,T,newtuple,newstep,j,filla,fillb,count,starter,SetToCompile,countprime,NewPosition,tt,OSc1,a,b},
If[(Mod[aindex,10]==0)&&(bindex==1),Print["NDA ",aindex," ",bindex," ",DateList[]]];
a=InitialDegrees[[aindex]];b=InitialDegrees[[bindex]];
OutSets={};counter=0;
While[counter≤Min[Length[a],Length[b]],
asubs=subslist[[aindex,counter+1]];
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]];
While[T≠{},
T1=First[First[Position[InitialDegrees,Reverse[Sort[T[[1,All,1]]]]]]];
T2=First[First[Position[InitialDegrees,Reverse[Sort[T[[1,All,2]]]]]]];
filla=InitialDegrees[[IDdiffs[[aindex,T1]]]];
fillb=InitialDegrees[[IDdiffs[[bindex,T2]]]];
newtuple=Reverse[Sort[Join[Table[T[[1,i,1]]+T[[1,i,2]],{i,1,counter}],filla,fillb]]];
TTT=Join[T[[1,1;;counter]],Table[{filla[[i]],0},{i,1,Length[filla]}],Table[{0,fillb[[i]]},{i,1,Length[fillb]}]];
TTTG=GatherBy[TTT,First[#]+Last[#]&];
newstep=Product[Apply[Multinomial,Tally[TTTG[[i]]][[All,2]]],{i,1,Length[TTTG]}]; (* 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]
];T=Drop[T,1]
];counter=counter+1
];
count=1;
While[count≤Length[OutSets],
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];
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]]]};
]; count=count+1
];
SSvals=Append[SSvals,OutSets]; (* seems they never repeat anyway *)
Length[SSvals]
];

Table[Table[NewDegreeAdder[p,q],{q,1,p}],{p,1,Length[InitialDegrees]}];

SuperSumSet=%;Print["SSS done at ",DateList[]];

TI3vals=Table[777,{z,1,300},{n,1,Length[SSvals]}];

binlist=Table[Table[1/(m*Binomial[n,m]),{m,1,n}],{n,1,50}];

(* The TermIntegral function tells us the integral over coeff*R_{k0} of (coeff-P1)^z P_{T}. *)
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_]:=Module[{i},
i=TI1vals[[z+1,T[[1]],T[[2]],T[[3]]]];
If[i==777,i=TermIntegral[kg-1,T,z,1-eg];TI1vals[[z+1,T[[1]],T[[2]],T[[3]]]]=i];
i
];

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

(* The TermIntegral2 function tells us the integral over (1-eps)*R_{k0} of (1+eps-P1)^z P_{T}, for signature T. *)
TermIntegral2[T_,z_]:=Module[{a,b,t,i},
a=1;b=1;t=0;i=0;
While[i≤z,
t=t+(a*b*TI1[T,z-i]);
a=a*2*eg;b=b*(z-i)/(i+1);
i=i+1
];
t
];

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

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

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

pdl=Table[Length[PossibleDegrees[d]],{d,0,30}];

pdls=Table[Sum[pdl[[i+1]],{i,j,k}],{j,0,30},{k,0,30}];

(*The zdegFinder helps us figure out what degree z and signature alpha we consider next for (1-P1)^z * P_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[ID_]:=
{
idc=1;cl2={};
While[idc≤Length[ID],
clist={};
tally1=Tally[ID[[idc]]];
tc=1;
While[tc≤Length[tally1],
zinta=ReplacePart[tally1,{tc,2}->Extract[tally1,{tc,2}]-1];
zintaprime={};
While[zinta≠{},
zintaprime=Join[zintaprime,Table[zinta[[1,1]],{m,1,zinta[[1,2]]}]];
zinta=Drop[zinta,1]
];
clist=Append[clist,First[First[Position[ID,zintaprime]]]];
tc=tc+1
];
cl2=Append[cl2,clist];idc=idc+1
];
cl2
}

(* I build matrices M1, M2, as in James' paper, and find the largest eigenvalue of Inverse[M1]M2. *)
Output:=
{
IDkids=First[Children[InitialDegrees]];
varsnum=Sum[Length[PossibleDegrees[m]],{m,0,dg}];
(* varsnum counts the total number of basis elements I will consider.  Remember that we are doubling the signatures, and then using powers of (1-P1) to "replace" those basis elements which had 1's in the signature. *)
M1=Table[Table[0,{pp,1,varsnum}],{qq,1,varsnum}];
M2=M1;
(* We have now constructed empty matrices of the appropriate sizes.  We will fill in each entry successively. *)
ii=1;Print[varsnum];
Print[DateList[]];
While[ii≤varsnum,
If[Mod[ii,10]==0,Print[ii," ",DateList[]]];
{zadeg,nicei}=zdegFinder[ii];
jj=1;
While[jj≤ii,
{zbdeg,nicej}=zdegFinder[jj];
SumSetVal=SSvals[[SuperSumSet[[Max[nicei,nicej],Min[nicei,nicej]]]]];
M1=ReplacePart[M1,{ii,jj}->Sum[SumSetVal[[m,2]]*TI1b[SumSetVal[[m,1]],zadeg+zbdeg-2],{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]];
tally1count=0;
While[tally1count≤Length[tally1],
tally2count=0;
While[tally2count≤Length[tally2],
(* There are 4 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 formulas below follow James' 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]/(zadeg*zbdeg),
M2cell=M2cell+binlist[[zbdeg+t21,zbdeg]]*TI3[SumSet,zadeg+zbdeg+t21]/zadeg
],
If[tally2count==0,
M2cell=M2cell+binlist[[zadeg+t11,zadeg]]*TI3[SumSet,zadeg+zbdeg+t11]/zbdeg,
M2cell=M2cell+binlist[[zadeg+t11,zadeg]]*binlist[[zbdeg+t21,zbdeg]]*TI3[SumSet,zadeg+zbdeg+t11+t21]
]
];tally2count=tally2count+1
];tally1count=tally1count+1
];
M2=ReplacePart[M2,{ii,jj}->M2cell];
If[jj<ii,M1=ReplacePart[M1,{jj,ii}->Extract[M1,{ii,jj}]];M2=ReplacePart[M2,{jj,ii}->Extract[M2,{ii,jj}]]];
jj=jj+1
];
ii=ii+1
];
(* We are now ready to compute the end result, by an iterative positive-definite test. *)
Print["Matrices filled at ",DateList[]];
If[PositiveDefiniteMatrixQ[Mgoal/kg*N[M1,40]-N[M2,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]]>Mrange/100000,
pdvals[[2,1]]=(pdvals[[3,1]]+pdvals[[1,1]])/2;
pdvals[[2,2]]=PositiveDefiniteMatrixQ[pdvals[[2,1]]/kg*N[M1,40]-N[M2,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[]]
]

(* Previously, M was determined in one step by calculating a generalised eigenvector, as shown below, but this is computationally more intensive than the iterative positive-definite test implemented above. *)

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

Output