(* PossibleDegrees finds all possible signatures up to total degree d0. *)
PossibleDegrees[d0_]:=Prepend[Flatten[Union[Table[IntegerPartitions[i],{i,1,d0}]],1],{}]

{};
SSvals=%;
SSsigvals=%;

NewDegreeAdder[a_,b_]:=
Module[{counter,OutSets,asubs,perms,bsubs,newperms,T,newtuple,newstep,j,filla,fillb,count,starter,SetToCompile,countprime,NewPosition,tt},
counter=0;
OutSets={};
While[counter≤Min[Length[a],Length[b]],
asubs=Union[Subsets[a,{counter}]];
perms=Apply[Union,Table[Permutations[asubs[[i]]],{i,1,Length[asubs]}]];
bsubs=Union[Subsets[b,{counter}]];
T={};
While[bsubs≠{},
newperms=perms;
While[newperms≠{},
T=Append[T,Sort[Table[{newperms[[1,i]],bsubs[[1,i]]},{i,1,counter}]]];
newperms=Drop[newperms,1]
];
bsubs=Drop[bsubs,1]
];
T=Union[T];
While[T≠{},
newtuple=Reverse[Sort[Apply[Join,{Table[T[[1,i,1]]+T[[1,i,2]],{i,1,counter}],
j=1;filla=a;While[j≤counter,filla=Drop[filla,First[Position[filla,T[[1,j,1]]]]];j=j+1];filla,
j=1;fillb=b;While[j≤counter,fillb=Drop[fillb,First[Position[fillb,T[[1,j,2]]]]];j=j+1];fillb}]]];
newstep=1;
count=1;
While[count≤Length[Union[newtuple]],
starter=Union[newtuple][[count]];
SetToCompile={};
countprime=1;
While[countprime≤counter,
If[T[[1,countprime,1]]+T[[1,countprime,2]]==starter,
SetToCompile=Append[SetToCompile,{T[[1,countprime,1]],T[[1,countprime,2]]}],
];
countprime=countprime+1
];
While[countprime≤counter+Length[filla],
If[filla[[countprime-counter]]==starter,
SetToCompile=Append[SetToCompile,{filla[[countprime-counter]],0}]
];
countprime=countprime+1
];
While[countprime≤counter+Length[filla]+Length[fillb],
If[fillb[[countprime-counter-Length[filla]]]==starter,
SetToCompile=Append[SetToCompile,{0,fillb[[countprime-counter-Length[filla]]]}]
];
countprime=countprime+1
];
newstep=newstep*Apply[Multinomial,Tally[SetToCompile][[All,2]]];
count=count+1
];
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;
OutSets=Sort[OutSets];
While[count≤Length[OutSets],
If[MemberQ[SSsigvals,OutSets[[count,1]]],
OutSets[[count,1]]=First[First[Position[SSsigvals,OutSets[[count,1]]]]],
SSsigvals=Append[SSsigvals,OutSets[[count,1]]];
OutSets[[count,1]]=Length[SSsigvals]
]; count=count+1
];
If[MemberQ[SSvals,OutSets],
tt=First[First[Position[SSvals,OutSets]]],
SSvals=Append[SSvals,OutSets];
tt=Length[SSvals]
];
tt
];

InitialDegrees=PossibleDegrees[5];Print["ID done ",DateList[]];

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

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

Table[777,{z,1,300},{n,1,Length[SSsigvals]}];
TI1vals=%;
TI1valsb=%;
TI2vals=%;
Table[777,{z,1,300},{n,1,Length[SSvals]}];
TI3vals=%;

(* 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]]]]

Table[Factorial[n],{n,0,200}];

faclist=%;

Table[Product[faclist[[1+2*SSsigvals[[j,i]]]],{i,1,Length[SSsigvals[[j]]]}],{j,1,Length[SSsigvals]}];

fplist=%;

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

binlist=%;

(* The TermIntegral function tells us the integral over coeff*R_{k0} of (coeff-P1)^z P_{T}, for signature T DOUBLED!!!!!!!!!!!!. *)
TermIntegral[k0_,T2_,z_,coeff_]:=Module[{T,T3,x},
T=2*SSsigvals[[T2]];T3=z+Total[T]+k0;
x=Power[coeff,T3]*faclist[[1+z]]*fplist[[T2]]*CountRepeats[k0,T]/faclist[[1+T3]];
x
];

TI1[k_,T_,z_,eps_]:=Module[{i},
i=TI1vals[[z+1,T]];
If[i==777,i=TermIntegral[k,T,z,eps];TI1vals[[z+1,T]]=i];
i
];

TI1b[k_,T_,z_,eps_]:=Module[{i},
i=TI1valsb[[z+1,T]];
If[i==777,i=TermIntegral[k,T,z,eps];TI1valsb[[z+1,T]]=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[k0_,T_,z_,eps_]:=Module[{a,b,t,i},
a=1;b=1;t=0;i=0;
While[i≤z,
t=t+(a*b*TI1[k0-1,T,z-i,1-eps]);
a=a*2*eps;b=b*(z-i)/(i+1);
i=i+1
];
t
];

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

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

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

(*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_,d0_]:=
{
If[iplug≤Length[PossibleDegrees[d0]],
outplace=iplug;
outzdeg=0,
If[iplug≤2(Length[PossibleDegrees[d0]]),
outplace=iplug-(Length[PossibleDegrees[d0]]);
outzdeg=1,
yy=1;
While[iplug>2(Length[PossibleDegrees[d0]])+2*Sum[Length[PossibleDegrees[m]],{m,d0-yy,d0-1}],
yy=yy+1
];
uu=0;
If[iplug>2(Length[PossibleDegrees[d0]])+2*Sum[Length[PossibleDegrees[m]],{m,d0-yy+1,d0-1}]+Length[PossibleDegrees[d0-yy]],
uu=1
];
outzdeg=2*yy-uu+1;
outplace=iplug-2(Length[PossibleDegrees[d0]])-2*Sum[Length[PossibleDegrees[m]],{m,d0-yy+1,d0-1}]-uu*Length[PossibleDegrees[d0-yy]]
]
];
{outzdeg+1,outplace}
}

(* I build matrices M1, M2, as in James' paper, and find the largest eigenvalue of Inverse[M1]M2. *)
Output[k_,d_,eps_]:=
{
InitialDegrees=PossibleDegrees[d];
varsnum=2*Length[InitialDegrees]+2*Sum[Length[PossibleDegrees[m]],{m,0,d-1}];
(* 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]];
{zadeg,nicei}=First[zdegFinder[ii,d]];
jj=1;
While[jj≤ii,
{zbdeg,nicej}=First[zdegFinder[jj,d]];
SumSetVal=SSvals[[SuperSumSet[[Max[nicei,nicej],Min[nicei,nicej]]]]];
M1=ReplacePart[M1,{ii,jj}->Sum[SumSetVal[[m,2]]*TI1b[k,SumSetVal[[m,1]],zadeg+zbdeg-2,1+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]];
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=2*tally1[[tally1count,1]];
zinta=ReplacePart[tally1,{tally1count,2}->Extract[tally1,{tally1count,2}]-1];
zintaprime={};
While[zinta≠{},
zintaprime=Join[zintaprime,Table[zinta[[1,1]],{m,1,zinta[[1,2]]}]];
zinta=Drop[zinta,1]
];
newi=First[First[Position[InitialDegrees,zintaprime]]]
];
If[tally2count==0,
newj=nicej,
t21=2*tally2[[tally2count,1]];
zintb=ReplacePart[tally2,{tally2count,2}->Extract[tally2,{tally2count,2}]-1];
zintbprime={};
While[zintb≠{},
zintbprime=Join[zintbprime,Table[zintb[[1,1]],{m,1,zintb[[1,2]]}]];
zintb=Drop[zintb,1]
];
newj=First[First[Position[InitialDegrees,zintbprime]]]
];
SumSet=SuperSumSet[[Max[newi,newj],Min[newi,newj]]];
If[tally1count==0,
If[tally2count==0,
M2cell=M2cell+((1/zadeg))*((1/zbdeg))*TI3[k,SumSet,zadeg+zbdeg,eps],
M2cell=M2cell+((1/zadeg))*binlist[[zbdeg+t21+1,zbdeg+1]]*TI3[k,SumSet,zadeg+zbdeg+t21,eps]
],
If[tally2count==0,
M2cell=M2cell+binlist[[zadeg+t11+1,zadeg+1]]*((1/zbdeg))*TI3[k,SumSet,zadeg+zbdeg+t11,eps],
M2cell=M2cell+binlist[[zadeg+t11+1,zadeg+1]]*binlist[[zbdeg+t21+1,zbdeg+1]]*TI3[k,SumSet,zadeg+zbdeg+t11+t21,eps]
]
];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 finding our eigenvector. *)
Print["Matrices filled at ",DateList[]];
RatVec=Rationalize[Eigenvectors[{k*N[M2,40],N[M1,40]},1][[1]],10^(-40)];
Print["RatVec is done at ",DateList[]];
Ratio=k*(RatVec.M2.RatVec)/(RatVec.M1.RatVec);
Print["Ratio done at ",DateList[]];
N[Ratio,20]
}

DateList[]

Output[49,4,1/48]

DateList[]