%% Gives back combinatorial factors for a given n=number of replica

function [S,a1,a2]=combinatorial_factors(n)

%% One-loop self-energy
S_2_0=(n.^3-9*n.^2+26*n-22)./(2*(n-1).*(n-2).^2);
S_1_1=(3*n.^2-15*n+16)./(2*(n-1).*(n-2).^2);
S_0_2=(n.^4-8*n.^3+19*n.^2-4*n-16)./(4*(n-1).*(n-2).^2);

S=[S_2_0,S_1_1;S_1_1,S_0_2];


%% One-loop cubic vertices
aI_3_0=(n.^3-11*n.^2+38*n-34)./(2*(n-1).*(n-2).^2);
aII_3_0=(-1)./((n-2).^3);
aI_2_1=(3*n.^2-19*n+20)./(2*(n-1).*(n-2).^2);
aII_2_1=(-n.^3+8*n.^2-17*n+12)./(2*(n-1).*(n-2).^3);
aI_1_2=(-n.^3+5*n.^2+8*n-16)./(4*(n-1).*(n-2).^2);
aII_1_2=(3*n.^3-27*n.^2+64*n-48)./(4*(n-1).*(n-2).^3);
aI_0_3=(-3*n)./(2*(n-2).^2);
aII_0_3=(n.^5-10*n.^4+33*n.^3-8*n.^2-104*n+112)./(8*(n-1).*(n-2).^3);

a1=zeros(2,2,2,2);
a1(1,1,1,1)=aI_3_0;
a1(1,1,1,2)=aI_2_1; a1(1,1,2,1)=aI_2_1; a1(1,2,1,1)=aI_2_1;
a1(1,2,2,1)=aI_1_2; a1(1,2,1,2)=aI_1_2; a1(1,1,2,2)=aI_1_2;
a1(1,2,2,2)=aI_0_3;
a1(2,1,1,1)=aII_3_0;
a1(2,1,1,2)=aII_2_1; a1(2,1,2,1)=aII_2_1; a1(2,2,1,1)=aII_2_1;
a1(2,2,2,1)=aII_1_2; a1(2,2,1,2)=aII_1_2; a1(2,1,2,2)=aII_1_2;
a1(2,2,2,2)=aII_0_3;


%% Two-loop cubic vertices
aI_5_0=(n.^8-26*n.^7+291*n.^6-1816*n.^5+6840*n.^4-15756*n.^3+21586*n.^2-16088*n+5008)./(4*(n-1).^2.*(n-2).^6);
aI_4_1=(3*n.^7-66*n.^6+607*n.^5-2960*n.^4+8132*n.^3-12592*n.^2+10236*n-3392)./(4*(n-1).^2.*(n-2).^6);
bI_4_1=(3*n.^7-66*n.^6+604*n.^5-2930*n.^4+8017*n.^3-12380*n.^2+10048*n-3328)./(4*(n-1).^2.*(n-2).^6);
aI_3_2=(21*n.^6-366*n.^5+2493*n.^4-8316*n.^3+14536*n.^2-12800*n+4480)./(8*(n-1).^2.*(n-2).^6);
bI_3_2=(3*n.^7-27*n.^6-59*n.^5+1471*n.^4-6396*n.^3+12496*n.^2-11664*n+4224)./(8*(n-1).^2.*(n-2).^6);
cI_3_2=(n.^7-7*n.^6-63*n.^5+819*n.^4-3292*n.^3+6262*n.^2-5776*n+2080)./(4*(n-1).^2.*(n-2).^6);
aI_2_3=(n.^9-19*n.^8+145*n.^7-541*n.^6+1018*n.^5-1488*n.^4+4292*n.^3-10192*n.^2+11328*n-4608)./(16*(n-1).^2.*(n-2).^6);
bI_2_3=(-n.^7+20*n.^6-110*n.^5+84*n.^4+871*n.^3-2704*n.^2+3040*n-1216)./(4*(n-1).^2.*(n-2).^6);
cI_2_3=(-7*n.^7+134*n.^6-819*n.^5+1708*n.^4+680*n.^3-7552*n.^2+10144*n-4352)./(16*(n-1).^2.*(n-2).^6);
aI_1_4=(n.^9-15*n.^8+95*n.^7-469*n.^6+2196*n.^5-6368*n.^4+8592*n.^3-2176*n.^2-5376*n+3584)./(32*(n-1).^2.*(n-2).^6);
bI_1_4=(3*n.^8-42*n.^7+169*n.^6+68*n.^5-1750*n.^4+3488*n.^3-1456*n.^2-1984*n+1536)./(16*(n-1).^2.*(n-2).^6);
aI_0_5=n.*(-3*n.^6+54*n.^5-315*n.^4+560*n.^3+376*n.^2-1968*n+1440)./(16*(n-1).*(n-2).^6);

aII_5_0=3*(n.^2-7*n+8)./((n-1).^1.*(n-2).^5);
aII_4_1=(n.^5-15*n.^4+78*n.^3-165*n.^2+159*n-62)./(2*(n-1).^2.*(n-2).^5);
bII_4_1=(3*n.^5-42*n.^4+211*n.^3-448*n.^2+436*n-168)./(4*(n-1).^2.*(n-2).^5);
aII_3_2=(n.^7-18*n.^6+127*n.^5-420*n.^4+574*n.^3-40*n.^2-608*n+416)./(8*(n-1).^2.*(n-2).^5);
bII_3_2=(-n.^5+19*n.^4-118*n.^3+296*n.^2-336*n+148)./(2*(n-1).^2.*(n-2).^5);
cII_3_2=(-2*n.^5+41*n.^4-260*n.^3+659*n.^2-750*n+328)./(4*(n-1).^2.*(n-2).^5);
aII_2_3=(3*n.^5-72*n.^4+531*n.^3-1494*n.^2+1848*n-864)./(8*(n-1).^2.*(n-2).^5);
bII_2_3=(3*n.^6-39*n.^5+151*n.^4-45*n.^3-726*n.^2+1344*n-736)./(8*(n-1).^2.*(n-2).^5);
cII_2_3=(n.^7-14*n.^6+81*n.^5-352*n.^4+1412*n.^3-3384*n.^2+3984*n-1824)./(16*(n-1).^2.*(n-2).^5);
aII_1_4=(3*n.^5-17*n.^4-25*n.^3+243*n.^2-420*n+232)./(2*(n-1).^2.*(n-2).^5);
bII_1_4=(3*n.^6-24*n.^5+147*n.^4-1006*n.^3+3136*n.^2-4240*n+2112)./(16*(n-1).^2.*(n-2).^5);
aII_0_5=(3*n.^8-47*n.^7+315*n.^6-1229*n.^5+3110*n.^4-4088*n.^3+336*n.^2+4928*n-3648)./(32*(n-1).^2.*(n-2).^5);

a2=zeros(2,2,2,2,2,2);
a2(1,1,1,1,1,1)=aI_5_0;
a2(1,2,1,1,1,1)=aI_4_1; a2(1,1,2,1,1,1)=aI_4_1; a2(1,1,1,2,1,1)=aI_4_1;
a2(1,1,1,1,2,1)=bI_4_1; a2(1,1,1,1,1,2)=bI_4_1;
a2(1,1,1,1,2,2)=aI_3_2;
a2(1,2,2,1,1,1)=bI_3_2; a2(1,2,1,2,1,1)=bI_3_2; a2(1,1,2,2,1,1)=bI_3_2;
a2(1,2,1,1,2,1)=cI_3_2; a2(1,1,2,1,2,1)=cI_3_2; a2(1,1,1,2,2,1)=cI_3_2; a2(1,2,1,1,1,2)=cI_3_2; a2(1,1,2,1,1,2)=cI_3_2; a2(1,1,1,2,1,2)=cI_3_2;
a2(1,2,2,2,1,1)=aI_2_3;
a2(1,1,1,2,2,2)=bI_2_3; a2(1,1,2,1,2,2)=bI_2_3; a2(1,2,1,1,2,2)=bI_2_3;
a2(1,1,2,2,1,2)=cI_2_3; a2(1,2,1,2,1,2)=cI_2_3; a2(1,2,2,1,1,2)=cI_2_3; a2(1,1,2,2,2,1)=cI_2_3; a2(1,2,1,2,2,1)=cI_2_3; a2(1,2,2,1,2,1)=cI_2_3;
a2(1,1,2,2,2,2)=aI_1_4; a2(1,2,1,2,2,2)=aI_1_4; a2(1,2,2,1,2,2)=aI_1_4;
a2(1,2,2,2,1,2)=bI_1_4; a2(1,2,2,2,2,1)=bI_1_4;
a2(2,2,2,2,2,2)=aII_0_5;
a2(2,1,1,1,1,1)=aII_5_0;
a2(2,2,1,1,1,1)=aII_4_1; a2(2,1,2,1,1,1)=aII_4_1; a2(2,1,1,2,1,1)=aII_4_1;
a2(2,1,1,1,2,1)=bII_4_1; a2(2,1,1,1,1,2)=bII_4_1;
a2(2,1,1,1,2,2)=aII_3_2;
a2(2,2,2,1,1,1)=bII_3_2; a2(2,2,1,2,1,1)=bII_3_2; a2(2,1,2,2,1,1)=bII_3_2;
a2(2,2,1,1,2,1)=cII_3_2; a2(2,1,2,1,2,1)=cII_3_2; a2(2,1,1,2,2,1)=cII_3_2; a2(2,2,1,1,1,2)=cII_3_2; a2(2,1,2,1,1,2)=cII_3_2; a2(2,1,1,2,1,2)=cII_3_2;
a2(2,2,2,2,1,1)=aII_2_3;
a2(2,1,1,2,2,2)=bII_2_3; a2(2,1,2,1,2,2)=bII_2_3; a2(2,2,1,1,2,2)=bII_2_3;
a2(2,1,2,2,1,2)=cII_2_3; a2(2,2,1,2,1,2)=cII_2_3; a2(2,2,2,1,1,2)=cII_2_3; a2(2,1,2,2,2,1)=cII_2_3; a2(2,2,1,2,2,1)=cII_2_3; a2(2,2,2,1,2,1)=cII_2_3;
a2(2,1,2,2,2,2)=aII_1_4; a2(2,2,1,2,2,2)=aII_1_4; a2(2,2,2,1,2,2)=aII_1_4;
a2(2,2,2,2,1,2)=bII_1_4; a2(2,2,2,2,2,1)=bII_1_4;
a2(2,2,2,2,2,2)=aII_0_5;


end