Print["Math version of OLEC by J.Gluza, ver: 0.9\nlast modified 28 Oct 2012"]
Print["One Loop External Contractions of tensor integrals,  arXiv:1211.3921"]
Print["Some internal objects taken directly from hexagon.m by K.Kajda"]

BeginPackage["Reducer5pt`"]

p1::usage = "External momentum p1."
p2::usage = "External momentum p2."
p3::usage = "External momentum p3."
p4::usage = "External momentum p4."
p5::usage = "External momentum p5."
q::usage = "Chord"
g::usage = "Metric tensor"
mu1::usage = "Index"
mu2::usage = "Index"
mu3::usage = "Index"
mu3::usage = "Index"
mu4::usage = "Index"
  MCM::usage = ""
  YY::usage = ""
  qq::usage = ""
  CH::usage = ""
  I4::usage = ""
  I3::usage = ""
  I2::usage = ""
  SM::usage = ""
  qaQ0::usage = ""
  qaQS::usage = ""
  R::usage = ""
  bracket::usage = ""
  C5::usage = ""

RedE0::usage = "Hexagon function. Reduces scalar five point function.\nUsage: RedE0[p1s,p2s,p3s,p4s,p5s,s12,s23,s34,s45,s15,m1s,m2s,m3s,m4s,m5s]\nwhere p1s=p1^2,...  m1s=m1^2,... and s12=(p1+p2)^2,..."
RedE1::usage = "Hexagon function. Reduces vector five point function and returns result as a sum over chords with coefficients.\nUsage: RedE1[p1s,p2s,p3s,p4s,p5s,s12,s23,s34,s45,s15,m1s,m2s,m3s,m4s,m5s]\nwhere p1s=p1^2,...  m1s=m1^2,... and s12=(p1+p2)^2,..."
RedE2::usage = "Hexagon function. Reduces tensor (rank 2) five point function and returns result as a basis formed from metric tensor g[mu,nu] and chords q[mu,i].\nUsage: RedE2[p1s,p2s,p3s,p4s,p5s,s12,s23,s34,s45,s15,m1s,m2s,m3s,m4s,m5s]\nwhere p1s=p1^2,...  m1s=m1^2,... and s12=(p1+p2)^2,..."
RedE3::usage = "Hexagon function. Reduces tensor (rank 3) five point function and returns result as a basis formed from metric tensor g[mu,nu] and chords q[mu,i].\nUsage: RedE3[p1s,p2s,p3s,p4s,p5s,s12,s23,s34,s45,s15,m1s,m2s,m3s,m4s,m5s]\nwhere p1s=p1^2,...  m1s=m1^2,... and s12=(p1+p2)^2,..."

CE1::usage = "hexagon 2012 chord version of vector"
CE2::usage = "hexagon 2012 chord version of rank 2"
CE2g::usage = "hexagon 2012 chord version of rank 2 contracted with metric tensor"
CE3::usage = "hexagon 2012 chord version of rank 3"
CE3g::usage = "hexagon 2012 chord version of rank 3 contracted with metric tensor"

RedEcoef::usage = "Returns specific tensor coefficient of five point function.\nUsage: RedEcoef[coef_,p1s,p2s,p3s,p4s,p5s,s12,s23,s34,s45,s15,m1s,m2s,m3s,m4s,m5s]\nwhere p1s=p1^2,...  m1s=m1^2,... and s12=(p1+p2)^2,..."
RedEget::usage = "Returns rules of coefficients of five point function.\nUsage: RedEget[rank_,p1s,p2s,p3s,p4s,p5s,s12,s23,s34,s45,s15,m1s,m2s,m3s,m4s,m5s]\nwhere p1s=p1^2,...  m1s=m1^2,... and s12=(p1+p2)^2,..."

D4::usage = "Four point scalar integral."
C3::usage = "Three point scalar integral."
B2::usage = "Two point scalar integral."

chords1::usage = ""
chords2::usage = ""
chords3::usage = ""
chords4::usage = ""
chords5::usage = ""

chords11::usage = ""
chords12::usage = ""
chords13::usage = ""
chords14::usage = ""
chords15::usage = ""
chords22::usage = ""
chords23::usage = ""
chords24::usage = ""
chords25::usage = ""
chords33::usage = ""
chords34::usage = ""
chords35::usage = ""
chords44::usage = ""
chords45::usage = ""
chords55::usage = ""

chords111::usage = ""
chords112::usage = ""
chords113::usage = ""
chords114::usage = ""
chords115::usage = ""
chords122::usage = ""
chords123::usage = ""
chords124::usage = ""
chords125::usage = ""
chords133::usage = ""
chords134::usage = ""
chords135::usage = ""
chords144::usage = ""
chords145::usage = ""
chords155::usage = ""
chords222::usage = ""
chords223::usage = ""
chords224::usage = ""
chords225::usage = ""
chords333::usage = ""
chords334::usage = ""
chords444::usage = ""
chords445::usage = ""
chords455::usage = ""
chords555::usage = ""

Begin["`Private`"]

Y[a_,b_]:=-(qq[a]-qq[b])^2+m[a]^2+m[b]^2;
YY[a_,b_]:=-CH[Expand[(qq[a]-qq[b])^2 /. chords]]+m[a]^2+m[b]^2;
TakeChordNum[i_, chords_] := ToExpression[StringTake[ToString[chords], {6 + i}]];


(******************************************)
(*       Modified Cayley determinant      *)
(******************************************)
MCM[size_] := Module[{matrix,i,j},

	Switch[size,
		5,
	matrix=
	{{0, 1, 1, 1, 1, 1}, 
         {1, Y[1,1], Y[1,2], Y[1,3], Y[1,4], Y[1,5]}, 
	 {1, Y[1,2], Y[2,2], Y[2,3], Y[2,4], Y[2,5]}, 
         {1, Y[1,3], Y[2,3], Y[3,3], Y[3,4], Y[3,5]}, 
         {1, Y[1,4], Y[2,4], Y[3,4], Y[4,4], Y[4,5]}, 
         {1, Y[1,5], Y[2,5], Y[3,5], Y[4,5], Y[5,5]}},
		6,
	matrix=
	{{0, 1, 1, 1, 1, 1, 1}, 
	 {1, Y[1,1], Y[1,2], Y[1,3], Y[1,4], Y[1,5], Y[1,6]}, 
         {1, Y[1,2], Y[2,2], Y[2,3], Y[2,4], Y[2,5], Y[2,6]}, 
         {1, Y[1,3], Y[2,3], Y[3,3], Y[3,4], Y[3,5], Y[3,6]}, 
         {1, Y[1,4], Y[2,4], Y[3,4], Y[4,4], Y[4,5], Y[4,6]}, 
         {1, Y[1,5], Y[2,5], Y[3,5], Y[4,5], Y[5,5], Y[5,6]}, 
         {1, Y[1,6], Y[2,6], Y[3,6], Y[4,6], Y[5,6], Y[6,6]}}
	];
		   

  Return[matrix]]


(******************************************)
(*              Determinants              *)
(******************************************)
SM[matrix_][Sequence[x__]][Sequence[y__]] := 
  Module[{sigx,sigy,rowSEQ,colSEQ,
	  norows,nocolumns,determinant},

	sigx=Signature[List[x]];
        sigy=Signature[List[y]];

        If[sigx==0 || sigy==0,
        {   
           determinant=0
        },
        {
        rowSEQ={#+1}& /@ List[Sequence[x]];
        colSEQ={#+1}& /@ List[Sequence[y]];

        norows=Delete[matrix,rowSEQ];
        nocolumns=Delete[#,colSEQ]& /@ norows;
        
        determinant=
          (-1)^(Plus[x,y])*
          sigx*sigy*Det[nocolumns]
        }];

  Return[determinant]]



(******************************************)
(*        definitions, rules, etc.        *)
(******************************************)
chords={
  qq[1]->0,
  qq[2]->p1,
  qq[3]->p1+p2,
  qq[4]->p1+p2+p3,
  qq[5]->p1+p2+p3+p4,
  qq[6]->p1+p2+p3+p4+p5
};

kinematics={
	  p1^2->ps1,
	  p2^2->ps2,
	  p3^2->ps3,
	  p4^2->ps4,
      p5^2->ps5,
	  m[1]^2->ms1, 
	  m[2]^2->ms2, 
	  m[3]^2->ms3, 
	  m[4]^2->ms4, 
	  m[5]^2->ms5};

scalars5pt={
          p1*p2->(s12-ps1-ps2)/2,
	  p2*p3->(s23-ps2-ps3)/2,
	  p3*p4->(s34-ps3-ps4)/2,
	  p4*p5->(s45-ps4-ps5)/2,
	  p1*p5->(s15-ps1-ps5)/2,
	  p3*p5->(s12-s34-s45+ps4)/2,
	  p1*p4->(s23-s15-s45+ps5)/2,
          p2*p5->(s34-s12-s15+ps1)/2,
          p1*p3->(s45-s12-s23+ps2)/2,
          p2*p4->(s15-s23-s34+ps3)/2};
 	  

(*********************************************)
(* Drops ee for tensor coefficients of  *)
(* 5pt function                        *)
(*********************************************)
DROPstr[x_, dropn_] := ToExpression[#] & /@ Characters[StringDrop[ToString[x], dropn]];
JOINstr[Sequence[x__]] := StringJoin[Map[ToString[#] &, List[x] ]];
CreateRules[list1_,list2_]:=Map[(Rule @@ #)&, Transpose[{list1,list2}]];

Options[RedE1]={Listable->False};
Options[RedE2]={Listable->False};
Options[RedE3]={Listable->False};

OptChange[fun_,new___Rule]:=
  Module[{i},
     Map[If[Length[i=Position[{new},First[#]->_,1,1]]===0, 
            #,{new}[[Sequence @@ First[i]]]]&, Options[fun]] ]
            
            
(*******************************************) 
(*        Chord reduction relations        *)
(*******************************************)           
qaQ0[a_] := -1/2*(YY[a, 1] - YY[1, 1]);
qaQs[a_, s_] := 1/2*(KroneckerDelta[a, s] - KroneckerDelta[1, s]);


(*******************************************)
(*           5pt chord vector              *)
(*******************************************)
CE1[chordsNum_, pq1_, pq2_, pq3_, pq4_, pq5_,
		 pq1pq2_, pq2pq3_, pq3pq4_, pq4pq5_, pq1pq5_,
		 mq1_, mq2_, mq3_, mq4_, mq5_]:=

 Module[{a, ind, E0,kinrules,kinematics5pt,cayley,result},
 
	a = TakeChordNum[1, chordsNum];

	E0 = RedE0[pq1, pq2, pq3, pq4, pq5, pq1pq2, pq2pq3, pq3pq4, pq4pq5, pq1pq5, mq1, mq2, mq3, mq4, mq5];
	
	kinrules={ps1->pq1,ps2->pq2,ps3->pq3,ps4->pq4,ps5->pq5,
		  ms1->mq1,ms2->mq2,ms3->mq3,ms4->mq4,ms5->mq5,
		  s12->pq1pq2,s23->pq2pq3,s34->pq3pq4,
		  s45->pq4pq5,s15->pq1pq5};
	kinematics5pt=Join[kinematics,scalars5pt] /.kinrules;
	cayley=Expand[(MCM[5] /. chords)]/.kinematics5pt;	

	ind = a + 1;
	If[ind == 6, ind = 1];

	(* CHANGED MINUS SIGN *)
	result = -(E0*qaQ0[ind] - Sum[I4[s]*qaQs[ind, s], {s, 1, 5}]) /. CH[n_] -> n;
	
  Return[Expand[result] /. kinematics5pt]]


CE2[chordsNum_,pq1_,pq2_,pq3_,pq4_,pq5_,
      pq1pq2_,pq2pq3_,pq3pq4_,pq4pq5_,pq1pq5_,
      mq1_,mq2_,mq3_,mq4_,mq5_]:=

 Module[{a,b,ind1,ind2,E1,kinrules,kinematics5pt,cayley,result},
 
	a = TakeChordNum[1, chordsNum];
	b = TakeChordNum[2, chordsNum];

	E1 = -CE1[ToExpression[StringTake[ToString[chordsNum], 7]], pq1, pq2, pq3, pq4, pq5, pq1pq2, pq2pq3, pq3pq4, pq4pq5, pq1pq5, mq1, mq2, mq3, mq4, mq5];

	
	kinrules={ps1->pq1,ps2->pq2,ps3->pq3,ps4->pq4,ps5->pq5,
		  ms1->mq1,ms2->mq2,ms3->mq3,ms4->mq4,ms5->mq5,
		  s12->pq1pq2,s23->pq2pq3,s34->pq3pq4,
		  s45->pq4pq5,s15->pq1pq5};
	kinematics5pt=Join[kinematics,scalars5pt] /.kinrules;
	cayley=Expand[(MCM[5] /.chords)]/.kinematics5pt;
				
	ind1 = a + 1;
	ind2 = b + 1;	
	If[ind1 == 6, ind1 = 1];
	If[ind2 == 6, ind2 = 1];				

	qaQs0[i_, s_] := 1/SM[cayley][s][s]*
		(-1/2*(SM[cayley][s][s]*(YY[i, 1] - YY[1, 1]) + SM[cayley][s][0]*(KroneckerDelta[i, s] - KroneckerDelta[1, s])));				
				
	qaQst[i_, s_, t_] := 1/SM[cayley][s][s]*
		(1/2*(1 - KroneckerDelta[s, t])*(SM[cayley][s][s]*(KroneckerDelta[i, t] - KroneckerDelta[1, t]) - SM[cayley][s][t]*(KroneckerDelta[i, s] - KroneckerDelta[1, s])));

	result = Expand[E1*qaQ0[ind2] - 
			 Sum[(qaQs0[ind1, s]*I4[s] - Sum[If[s != t, qaQst[ind1, s, t]*I3[s, t], 0], {t, 1, 5}])*qaQs[ind2, s],{s, 1, 5}]] /. CH[n_] -> n;

  Return[result /.kinematics5pt]]
  

CE2g[pq1_,pq2_,pq3_,pq4_,pq5_,
      pq1pq2_,pq2pq3_,pq3pq4_,pq4pq5_,pq1pq5_,
      mq1_,mq2_,mq3_,mq4_,mq5_]:=

 Module[{E0,kinrules,kinematics5pt,cayley,result},

	E0 = RedE0[pq1, pq2, pq3, pq4, pq5, pq1pq2, pq2pq3, pq3pq4, pq4pq5, pq1pq5, mq1, mq2, mq3, mq4, mq5];
	
	kinrules={ps1->pq1,ps2->pq2,ps3->pq3,ps4->pq4,ps5->pq5,
		  ms1->mq1,ms2->mq2,ms3->mq3,ms4->mq4,ms5->mq5,
		  s12->pq1pq2,s23->pq2pq3,s34->pq3pq4,
		  s45->pq4pq5,s15->pq1pq5};
	kinematics5pt=Join[kinematics,scalars5pt] /.kinrules;
	cayley=Expand[(MCM[5] /.chords)]/.kinematics5pt;
				
	result = YY[1, 1]*E0/2 + I4[1] /. CH[n_] -> n;
	
  Return[result /.kinematics5pt]]  

  
CE3[chordsNum_,pq1_,pq2_,pq3_,pq4_,pq5_,
      pq1pq2_,pq2pq3_,pq3pq4_,pq4pq5_,pq1pq5_,
      mq1_,mq2_,mq3_,mq4_,mq5_]:=

 Module[{E2,kinrules,kinematics5pt,cayley,result},

	a = TakeChordNum[1, chordsNum];
	b = TakeChordNum[2, chordsNum];
	c = TakeChordNum[3, chordsNum];

	E2 = -CE2[ToExpression[StringTake[ToString[chordsNum], 8]], pq1, pq2, pq3, pq4, pq5, pq1pq2, pq2pq3, pq3pq4, pq4pq5, pq1pq5, mq1, mq2, mq3, mq4, mq5];
	
	ind1 = a + 1;
	ind2 = b + 1;
	ind3 = c + 1;
	
	If[ind1 == 6, ind1 = 1];
	If[ind2 == 6, ind2 = 1];
	If[ind3 == 6, ind3 = 1];
	
	kinrules={ps1->pq1,ps2->pq2,ps3->pq3,ps4->pq4,ps5->pq5,
		  ms1->mq1,ms2->mq2,ms3->mq3,ms4->mq4,ms5->mq5,
		  s12->pq1pq2,s23->pq2pq3,s34->pq3pq4,
		  s45->pq4pq5,s15->pq1pq5};
	kinematics5pt=Join[kinematics,scalars5pt] /.kinrules;
	cayley=Expand[(MCM[5] /.chords)]/.kinematics5pt;

	R[s_] := 1/SM[cayley][s][s]*(SM[cayley][s][0]*I4[s] - Sum[If[s != t, SM[cayley][s][t]*I3[s, t], 0], {t, 1, 5}]) /.kinematics5pt;
	R[s_, t_] := 1/SM[cayley][s, t][s, t]*(SM[cayley][s, t][0, t]*I3[s, t] - Sum[If[s != u && t != u && s != t, SM[cayley][s, t][u, t]*I2[s, t, u], 0], {u, 1, 5}]) /.kinematics5pt;

	BRACKET[s_] := 1/SM[cayley][s][s]*(
				 - (SM[cayley][0][0]*I4[s] - Sum[If[s != t, SM[cayley][t][0]*I3[s, t], 0], {t, 1, 5}] - SM[cayley][s][0]*R[s])
				 + SM[cayley][s][0]*R[s]
				 - Sum[If[s != t, SM[cayley][s][t]*R[s, t], 0], {t, 1, 5}]) /.kinematics5pt;

	YS[i_] := YY[i, 1] - YY[1, 1];
	
	C5[i_, j_, k_] := 1/8*(
				- KroneckerDelta[i, j]*KroneckerDelta[i, k]*BRACKET[i] 
				+ BRACKET[1]
				- YS[i]*YS[j]*I4[k] 
				+ YS[i]*YS[j]*I4[1] 
				- If[j != k, YS[i]*I3[j, k], 0] 
				- If[i != k, YS[j]*I3[i, k], 0]
				+ If[j != 1, YS[i]*I3[j, 1], 0]
				+ If[k != 1, YS[i]*I3[k, 1], 0]
				+ If[i != 1, YS[j]*I3[i, 1], 0]
				+ If[k != 1, YS[j]*I3[k, 1], 0]
				+ If[i != j && i != 1 && j != 1, I2[i, j, 1], 0]
				+ If[i != k && i != 1 && k != 1, I2[i, k, 1], 0]
				+ If[j != k && k != 1 && j != 1, I2[j, k, 1], 0]
				- If[i != j && i != k && j != k, I2[i, j, k], 0]
				- KroneckerDelta[i, k]*YS[j]*R[i]
				- KroneckerDelta[j, k]*YS[i]*R[j]
				- YS[i]*R[1]
				- YS[j]*R[1]
				- If[i != 1, R[1, i], 0]
				- If[j != 1, R[1, j], 0]
				- If[k != 1, R[1, k], 0]
				+ If[i != 1, KroneckerDelta[i, j]*R[i, 1], 0]
				+ If[i != 1, KroneckerDelta[i, k]*R[i, 1], 0]
				- If[i != j, KroneckerDelta[i, k]*R[i, j], 0]
				- If[i != k, KroneckerDelta[i, j]*R[i, k], 0]
				- If[j != i, KroneckerDelta[j, k]*R[j, i], 0]
				+ If[j != 1, KroneckerDelta[j, k]*R[j, 1], 0]);

	(* Result is zero when ind1 or ind2 or ind3 points to zero chord *)
	result = If[ind1 == 1 || ind2 == 1 || ind3 == 1, 0, E2*qaQ0[ind3] - C5[ind1, ind2, ind3] /. CH[n_]->n];

  Return[result /.kinematics5pt]]  
  

CE3g[chordsNum_,pq1_,pq2_,pq3_,pq4_,pq5_,
      pq1pq2_,pq2pq3_,pq3pq4_,pq4pq5_,pq1pq5_,
      mq1_,mq2_,mq3_,mq4_,mq5_]:=

 Module[{a, ind, E0,kinrules,kinematics5pt,cayley,result},
	a = TakeChordNum[1, chordsNum];

	E0 = RedE0[pq1, pq2, pq3, pq4, pq5, pq1pq2, pq2pq3, pq3pq4, pq4pq5, pq1pq5, mq1, mq2, mq3, mq4, mq5];
	
	kinrules={ps1->pq1,ps2->pq2,ps3->pq3,ps4->pq4,ps5->pq5,
		  ms1->mq1,ms2->mq2,ms3->mq3,ms4->mq4,ms5->mq5,
		  s12->pq1pq2,s23->pq2pq3,s34->pq3pq4,
		  s45->pq4pq5,s15->pq1pq5};
	kinematics5pt=Join[kinematics,scalars5pt] /.kinrules;
	cayley=Expand[(MCM[5] /.chords)]/.kinematics5pt;
R[s_] := 1/SM[cayley][s][s]*(SM[cayley][s][0]*I4[s] - Sum[If[s != t, SM[cayley][s][t]*I3[s, t], 0], {t, 1, 5}]) /.kinematics5pt;				
	ind = a + 1;
	If[ind == 6, ind = 1];				
	
	qaQ0s[i_, s_] := 1/SM[cayley][0][0]*
		(1/2*(SM[cayley][s][0]*(YY[i, 1] - YY[1, 1]) + SM[cayley][0][0]*(KroneckerDelta[i, s] - KroneckerDelta[1, s])));				
		
	qaQs0[i_, s_] := 1/SM[cayley][s][s]*
		(-1/2*(SM[cayley][s][s]*(YY[i, 1] - YY[1, 1]) + SM[cayley][s][0]*(KroneckerDelta[i, s] - KroneckerDelta[1, s])));				
				
	qaQst[i_, s_, t_] := 1/SM[cayley][s][s]*
		(1/2*(1 - KroneckerDelta[s, t])*(SM[cayley][s][s]*(KroneckerDelta[i, t] - KroneckerDelta[1, t]) - SM[cayley][s][t]*(KroneckerDelta[i, s] - KroneckerDelta[1, s])));
		
	(* CHANGED MINUS SIGN, connected with def. of YY *)
	(*	result = YY[1, 1]/2*Sum[I4[s]*qaQ0s[ind, s], {s, 1, 5}] - I4[1]*qaQs0[ind, 1] + Sum[I3[1, t]*qaQst[ind, 1, t], {t, 2, 5}] /. CH[n_]->n;
	 *)
result = If[ind == 1, 0, 
YY[1, 1]/4 * ((YY[ind, 1] - YY[1, 1]) * E0 + I4[ind] - I4[1]) + 
1/2 * ((YY[ind, 1] - YY[1, 1]) * I4[1] + I3[ind, 1] - R[1]) /. CH[n_] -> n];

  Return[result /.kinematics5pt]]  
    

(******************************************)
(*         Scalar 5pt integral            *)
(******************************************)
RedE0[pq1_,pq2_,pq3_,pq4_,pq5_,
      pq1pq2_,pq2pq3_,pq3pq4_,pq4pq5_,pq1pq5_,
      mq1_,mq2_,mq3_,mq4_,mq5_]:=

 Module[{kinrules,kinematics5pt,cayley,result},

	kinrules={ps1->pq1,ps2->pq2,ps3->pq3,ps4->pq4,ps5->pq5,
		  ms1->mq1,ms2->mq2,ms3->mq3,ms4->mq4,ms5->mq5,
		  s12->pq1pq2,s23->pq2pq3,s34->pq3pq4,
		  s45->pq4pq5,s15->pq1pq5};
	kinematics5pt=Join[kinematics,scalars5pt] /.kinrules;
	cayley=Expand[(MCM[5] /.chords)] /.kinematics5pt;

	result=
	  Sum[SM[cayley][0][s]*I4[s],{s,1,5}];

  Return[result/SM[cayley][0][0] /.kinematics5pt]]


(*******************************************)
(*               5pt vector                *)
(*******************************************)
RedE1[pq1_,pq2_,pq3_,pq4_,pq5_,
      pq1pq2_,pq2pq3_,pq3pq4_,pq4pq5_,pq1pq5_,
      mq1_,mq2_,mq3_,mq4_,mq5_,options___Rule]:=

 Module[{option,kinrules,kinematics5pt,cayley,result},

	option=OptChange[RedE1,options];
	kinrules={ps1->pq1,ps2->pq2,ps3->pq3,ps4->pq4,ps5->pq5,
		  ms1->mq1,ms2->mq2,ms3->mq3,ms4->mq4,ms5->mq5,
		  s12->pq1pq2,s23->pq2pq3,s34->pq3pq4,
		  s45->pq4pq5,s15->pq1pq5};
	kinematics5pt=Join[kinematics,scalars5pt] /.kinrules;
	cayley=Expand[(MCM[5] /.chords)]/.kinematics5pt;
       
		result=
	Table[q[mu1,i-1]*
	  Sum[SM[cayley][0,i][0,s]*I4[s],{s,1,5}], {i,2,5}]/SM[cayley][0][0];
	 
	If[Not[Listable/.option],result=Apply[Plus,result]];


  Return[result /.kinematics5pt]]


(*******************************************)
(*          5pt tensor 2nd order           *)
(*******************************************)
RedE2[pq1_,pq2_,pq3_,pq4_,pq5_,
      pq1pq2_,pq2pq3_,pq3pq4_,pq4pq5_,pq1pq5_,
      mq1_,mq2_,mq3_,mq4_,mq5_,options___Rule]:=

  Module[{option,kinrules,kinematics5pt,rules,cayley,
	  SM2D0dArray,SM2D00Array,SM1DnoArray,
	  SM2D0d,SM2D00,SM1Dno,II5,ID5,ArrayID5,E00,result},

	option=OptChange[RedE2,options];
	kinrules={ps1->pq1,ps2->pq2,ps3->pq3,ps4->pq4,ps5->pq5,
		  ms1->mq1,ms2->mq2,ms3->mq3,ms4->mq4,ms5->mq5,
		  s12->pq1pq2,s23->pq2pq3,s34->pq3pq4,
		  s45->pq4pq5,s15->pq1pq5};
	kinematics5pt=Join[kinematics,scalars5pt] /.kinrules;
	cayley=Expand[(MCM[5] /.chords)]/.kinematics5pt;
	
	SM2D0dArray=Table[SM[cayley][0,i][j,k],{i,1,5},{j,1,5},{k,1,5}];
	SM2D00Array=Table[SM[cayley][0,i][0,j],{i,1,5},{j,1,5}];
	SM1DnoArray=Table[SM[cayley][i][i],{i,1,5}];
	SM2D0d[i_,j_,k_]:=SM2D0dArray[[i,j,k]];
	SM2D00[i_,j_]:=SM2D00Array[[i,j]];
	SM1Dno[i_]:=SM1DnoArray[[i]];

	ID5[i_,j_]:=
	Sum[((-SM2D00[s,j]*SM2D0d[s,i,s]+SM2D0d[i,s,j]*SM2D00[s,s])*I4[s]-
	(Sum[If[s!=t,
            (-SM2D00[s,j]*SM[cayley][t,s][i,s]+SM2D0d[i,s,j]*SM2D0d[s,t,s])*I3[s,t],0],{t,1,5}])
	)/SM1Dno[s],{s,1,5}];

	ArrayID5=Flatten[Table[II5[i,j],{i,2,5},{j,i,5}]];
	E00=
	Sum[SM[cayley][0][s]/SM1Dno[s]*
	(SM2D00[s,s]*I4[s] - Sum[If[s!=t,SM2D0d[s,t,s]*I3[s,t],0],{t,1,5}]), 
		{s,1,5}]/2;
   
        If[Listable/.option,
	 result=Flatten[{-E00,ArrayID5}/SM[cayley][0][0] ]/.II5->ID5, 
	{result={-g[mu1,mu2]*E00,
	    Table[q[mu1,i-1]*q[mu2,j-1]*II5[Sequence @@ Sort[{i,j}] ],
	    	{i,2,5},{j,2,5}]}/SM[cayley][0][0];
	 result=Plus @@ Flatten[result /.II5->ID5]} ];
      

  Return[result/.kinematics5pt]]


(*******************************************)
(*          5pt tensor 3rd order           *)
(*******************************************)
RedE3[pq1_,pq2_,pq3_,pq4_,pq5_,
      pq1pq2_,pq2pq3_,pq3pq4_,pq4pq5_,pq1pq5_,
      mq1_,mq2_,mq3_,mq4_,mq5_,options___Rule]:=

  Module[{option,kinrules,kinematics5pt,cayley,SM3DnoArray,
          SM2DnoArray,SM2D00Array,SM2D0uArray,SM1DnoArray,
	  SM3Dno,SM2Dno,SM2D00,SM2D0u,SM1Dno,
	  S4,S3,S2,E00j,S4Array,S3Array,S2Array,E00jArray,result},

	option=OptChange[RedE3,options];
	kinrules={ps1->pq1,ps2->pq2,ps3->pq3,ps4->pq4,ps5->pq5,
		  ms1->mq1,ms2->mq2,ms3->mq3,ms4->mq4,ms5->mq5,
		  s12->pq1pq2,s23->pq2pq3,s34->pq3pq4,
		  s45->pq4pq5,s15->pq1pq5};
	kinematics5pt=Join[kinematics,scalars5pt] /.kinrules;
	cayley=Expand[(MCM[5] /.chords)]/.kinematics5pt;       
	
	SM3DnoArray=Table[SM[cayley][i,s,t][j,s,t],{i,0,5},{s,1,5},{t,1,5},{j,0,5}];
	SM2DnoArray=Table[SM[cayley][i,s][j,s],{i,1,5},{s,1,5},{j,1,5}];
	SM2D00Array=Table[SM[cayley][0,i][0,j],{i,1,5},{j,1,5}];
	SM2D0uArray=Table[SM[cayley][0,i][j,k],{i,1,5},{j,1,5},{k,1,5}];
	SM1DnoArray=Table[SM[cayley][i][j],{i,0,5},{j,0,5}];
	SM3Dno[i_,s_,t_,j_]:=SM3DnoArray[[i+1,s,t,j+1]];
	SM2Dno[i_,s_,j_]:=SM2DnoArray[[i,s,j]];
	SM2D00[i_,j_]:=SM2D00Array[[i,j]];
	SM2D0u[i_,j_,k_]:=SM2D0uArray[[i,j,k]];
	SM1Dno[i_,j_]:=SM1DnoArray[[i+1,j+1]];

	 S4[i_,j_,k_]:=
	 Sum[(-SM2D00[s,k]*
	  (SM2D0u[s,i,s]*SM2D0u[s,j,s] + SM2Dno[i,s,j]*SM2D00[s,s])+
	  (SM2D0u[i,s,k]*SM2D0u[s,j,s] + SM2D0u[j,s,k]*SM2D0u[s,i,s])*
	   SM2D00[s,s])/(SM1Dno[s,s])^2*I4[s], {s,1,5}];

	 S3[i_,j_,k_]:=
	 Sum[If[s!=t,
	  1/(SM1Dno[s,s])^2*
          (SM2D00[s,k]*
	  (SM2Dno[t,s,i]*SM2D0u[s,j,s] + SM2Dno[i,s,j]*SM2D0u[s,t,s]+
	   SM1Dno[s,s]*SM3Dno[0,s,t,i]*SM2Dno[t,s,j]/SM2Dno[s,t,s])-
          (SM2D0u[i,s,k]*SM2D0u[s,j,s] + SM2D0u[j,s,k]*SM2D0u[s,i,s])*
	   SM2D0u[s,t,s]-
	  (SM2D0u[i,s,k]*SM2Dno[t,s,j] + SM2D0u[j,s,k]*SM2Dno[t,s,i])*
	   SM1Dno[s,s]*SM3Dno[0,s,t,0]/(2*SM2Dno[s,t,s]))*
  		I3[s,t],0], {s,1,5},{t,1,5}];

	S2[i_,j_,k_]:=
	 Sum[If[s!=t && s!=u && t!=u,
	 -1/(SM1Dno[s,s]*SM2Dno[s,t,s])*
	  (SM2D00[s,k]*SM2Dno[t,s,j]*SM3Dno[u,s,t,i]-
	  (SM2D0u[j,s,k]*SM3Dno[u,s,t,i] + SM2D0u[i,s,k]*SM3Dno[u,s,t,j])*
           SM2D0u[s,t,s]/2)*I2[s,t,u], 0],{s,1,5},{t,1,5},{u,1,5}];

	E00j[j_]:=
       -Sum[(3*SM1Dno[s,0]*SM2D0u[s,j,s]-SM1Dno[s,j]*SM2D00[s,s])*
	        SM2D00[s,s]/(SM1Dno[s,s])^2 * I4[s],{s,1,5}]+
	Sum[If[s!=t,
	    (3*SM1Dno[s,0]*SM2D0u[s,j,s]-SM1Dno[s,j]*(SM2D0u[s,t,s])^2/SM2Dno[s,t,s])*
	       SM2D0u[s,t,s]/(SM1Dno[s,s])^2 * I3[s,t], 0],{s,1,5},{t,1,5}]-
	Sum[If[s!=t && s!=u && t!=u,
	    SM1Dno[s,j]*SM1Dno[s,s]*SM3Dno[u,s,t,0]*SM2D0u[s,t,s]/
	   (SM2Dno[s,t,s]*(SM1Dno[s,s])^2) * I2[s,t,u], 0],{s,1,5},{t,1,5},{u,1,5}];


	 S4Array=Table[S4[i,j,k],{i,2,5},{j,2,5},{k,2,5}]/.kinematics5pt;
	 S3Array=Table[S3[i,j,k],{i,2,5},{j,2,5},{k,2,5}]/.kinematics5pt;
	 S2Array=Table[S2[i,j,k],{i,2,5},{j,2,5},{k,2,5}]/.kinematics5pt;
	 

	If[Listable/.option,
	 result=-{Table[E00j[j],{j,2,5}]/2,
                  Table[S4Array[[i,j,k]]+S4Array[[k,j,i]]+S4Array[[i,k,j]]+
                        S3Array[[i,j,k]]+S3Array[[k,j,i]]+S3Array[[i,k,j]]+
		        S2Array[[i,j,k]]+S2Array[[k,j,i]]+S2Array[[i,k,j]],
	           {i,1,4},{j,i,4},{k,j,4}]}/(3*SM1Dno[0,0]) //Flatten, 
	{result={Table[q[mu1,i]*q[mu2,j]*q[mu3,k]*
		       (S4Array[[i,j,k]]+S4Array[[k,j,i]]+S4Array[[i,k,j]]+
                        S3Array[[i,j,k]]+S3Array[[k,j,i]]+S3Array[[i,k,j]]+
		        S2Array[[i,j,k]]+S2Array[[k,j,i]]+S2Array[[i,k,j]]),
	   	   {i,1,4},{j,1,4},{k,1,4}],
		 Table[(g[mu1,mu2]*q[mu3,j]+g[mu3,mu1]*q[mu2,j]+g[mu2,mu3]*q[mu1,j])*
			E00j[j+1], {j,1,4}]/2}/(3*SM1Dno[0,0]);
	 result=Plus @@ Flatten[-result]} 
	];


 Return[result/.kinematics5pt]]


(*******************************************)
(*           5pt coefficients              *)
(*******************************************)
RedEcoef[coef_,pq1_,pq2_,pq3_,pq4_,pq5_,
      pq1pq2_,pq2pq3_,pq3pq4_,pq4pq5_,pq1pq5_,
      mq1_,mq2_,mq3_,mq4_,mq5_]:=

 Module[{kinrules,kinematics5pt,IndList,i1,i2,i3,i,j,k,X,
         S4,S3,S2,SM3DnoArray,SM2DnoArray,SM2D00Array,
         SM2D0uArray,SM2D0dArray,SM1DnoArray,
	 SM3Dno,SM2Dno,SM2D00,SM2D0u,SM2D0d,SM1Dnocayley,result},


	kinrules={ps1->pq1,ps2->pq2,ps3->pq3,ps4->pq4,ps5->pq5,
		  ms1->mq1,ms2->mq2,ms3->mq3,ms4->mq4,ms5->mq5,
		  s12->pq1pq2,s23->pq2pq3,s34->pq3pq4,
		  s45->pq4pq5,s15->pq1pq5};
	kinematics5pt=Join[kinematics,scalars5pt] /.kinrules;

	cayley=Expand[(MCM[5] /.chords)]/.kinematics5pt;

        IndList=DROPstr[coef, 2];


	Switch[Length[IndList],
        	1,
	i1=First[IndList];
	result=
	 If[i1!=0,
	  Sum[SM[cayley][0,i1+1][0,s]*I4[s], {s,1,5}],
	  Sum[SM[cayley][0][s]*I4[s], {s,1,5}] ],

		2,
	i1=IndList[[1]];
	i2=IndList[[2]];
	result=
	 If[i1!=0 && i2!=0,
          X[s_,t_,i_,j_]:=-SM[cayley][0,s][0,i2+1]*SM[cayley][t,s][i1+1,s] + SM[cayley][0,i1+1][s,i2+1]*SM[cayley][t,s][0,s];
	  Sum[(X[s,0,i1+1,i2+1]*I4[s] - Sum[X[s,t,i1+1,i2+1]*I3[s,t], {t,1,5}])/SM[cayley][s][s], {s,1,5}],
	 -Sum[SM[cayley][s][0]/(2*SM[cayley][s][s])*
	     (SM[cayley][0,s][0,s]*I4[s] - Sum[SM[cayley][t,s][0,s]*I3[s,t], {t,1,5}]), {s,1,5}] ],

		3,
	i1=IndList[[1]];
	i2=IndList[[2]];
	i3=IndList[[3]];

	SM3DnoArray=Table[SM[cayley][i,s,t][j,s,t],{i,0,5},{s,1,5},{t,1,5},{j,0,5}];
	SM2DnoArray=Table[SM[cayley][i,s][j,s],{i,1,5},{s,1,5},{j,1,5}];
	SM2D00Array=Table[SM[cayley][0,i][0,j],{i,1,5},{j,1,5}];
	SM2D0uArray=Table[SM[cayley][0,i][j,k],{i,1,5},{j,1,5},{k,1,5}];
	SM2D0dArray=Table[SM[cayley][i,s][0,s],{i,1,5},{s,1,5}];
	SM1DnoArray=Table[SM[cayley][i][j],{i,0,5},{j,0,5}];

	SM3Dno[i_,s_,t_,j_]:=SM3DnoArray[[i+1,s,t,j+1]];
	SM2Dno[i_,s_,j_]:=SM2DnoArray[[i,s,j]];
	SM2D00[i_,j_]:=SM2D00Array[[i,j]];
	SM2D0u[i_,j_,k_]:=SM2D0uArray[[i,j,k]];
	SM2D0d[i_,s_]:=SM2D0dArray[[i,s]];
	SM1Dno[i_,j_]:=SM1DnoArray[[i+1,j+1]];

	 S4[i_,j_,k_]:=
	 Sum[(-SM2D00[s,k]*
	  (SM2D0u[s,i,s]*SM2D0u[s,j,s] + SM2Dno[i,s,j]*SM2D00[s,s])+
	  (SM2D0u[i,s,k]*SM2D0u[s,j,s] + SM2D0u[j,s,k]*SM2D0u[s,i,s])*
	   SM2D00[s,s])/(SM1Dno[s,s])^2*I4[s], {s,1,5}];

	 S3[i_,j_,k_]:=
	 Sum[If[s!=t,
	  1/(SM1Dno[s,s])^2*
          (SM2D00[s,k]*
	  (SM2Dno[t,s,i]*SM2D0u[s,j,s] + SM2Dno[i,s,j]*SM2D0d[t,s]+
	   SM1Dno[s,s]*SM3Dno[0,s,t,i]*SM2Dno[t,s,j]/SM2Dno[s,t,s])-
          (SM2D0u[i,s,k]*SM2D0u[s,j,s] + SM2D0u[j,s,k]*SM2D0u[s,i,s])*
	   SM2D0d[t,s]-
	  (SM2D0u[i,s,k]*SM2Dno[t,s,j] + SM2D0u[j,s,k]*SM2Dno[t,s,i])*
	   SM1Dno[s,s]*SM3Dno[0,s,t,0]/(2*SM2Dno[s,t,s]))*
  		I3[s,t],0], {s,1,5},{t,1,5}];

	 S2[i_,j_,k_]:=
	 Sum[If[s!=t && s!=u && t!=u,
	 -1/(SM1Dno[s,s]*SM2Dno[s,t,s])*
	  (SM2D00[s,k]*SM2Dno[t,s,j]*SM3Dno[u,s,t,i]-
	  (SM2D0u[j,s,k]*SM3Dno[u,s,t,i] + SM2D0u[i,s,k]*SM3Dno[u,s,t,j])*
           SM2D0d[t,s]/2)*I2[s,t,u], 0],
	 {s,1,5},{t,1,5},{u,1,5}];

	i=i1+1;
	j=i2+1;
	k=i3+1;

	result=
	-If[i1!=0 && i2!=0,
	  (S4[i,j,k]+S4[k,j,i]+S4[i,k,j]+
	   S3[i,j,k]+S3[k,j,i]+S3[i,k,j]+
           S2[i,j,k]+S2[k,j,i]+S2[i,k,j])/3,		
	 1/6*
	 (-Sum[(2*SM1Dno[s,0]*SM2D0u[s,k,s] - SM1Dno[s,s]*SM2D00[s,k])*
	          SM2D00[s,s]*I4[s]/(SM1Dno[s,s])^2, {s,1,5}]+
	  Sum[If[s!=t,
         ((SM1Dno[s,0]*SM2D0u[s,k,s] - SM1Dno[s,s]*SM2D00[s,k])*
	   SM2D0d[t,s]/(SM1Dno[s,s])^2+
	   SM1Dno[s,0]*SM2Dno[t,s,k]*SM2D00[s,s]/(SM1Dno[s,s])^2+
	   SM1Dno[s,0]*SM2D0d[t,s]*SM3Dno[0,s,t,k]/(SM1Dno[s,s]*SM2Dno[s,t,s]))*
		 I3[s,t], 0], {s,1,5},{t,1,5}]-
	  Sum[If[s!=t && s!=u && t!=u,
	   SM1Dno[s,0]*SM2D0d[t,s]*
	   SM3Dno[u,s,t,k]/(SM1Dno[s,s]*SM2Dno[s,t,s])*
		I2[s,t,u], 0], {s,1,5},{t,1,5},{u,1,5}]) ] ];


  Return[result/SM[cayley][0][0] /.kinematics5pt]]


(*******************************************)
(*           5pt all coefficients          *)
(*******************************************)
RedEget[rank_,pq1_,pq2_,pq3_,pq4_,pq5_,
      pq1pq2_,pq2pq3_,pq3pq4_,pq4pq5_,pq1pq5_,
      mq1_,mq2_,mq3_,mq4_,mq5_]:=

 Module[{kinrules,kinematics5pt,cayley,FivePT,rankDroped,
         coeficients,eecoef,ctr,fin,result={}},


	kinrules={ps1->pq1,ps2->pq2,ps3->pq3,ps4->pq4,ps5->pq5,
		  ms1->mq1,ms2->mq2,ms3->mq3,ms4->mq4,ms5->mq5,
		  s12->pq1pq2,s23->pq2pq3,s34->pq3pq4,
		  s45->pq4pq5,s15->pq1pq5};
	kinematics5pt=Join[kinematics,scalars5pt] /.kinrules;

	rankDroped=DROPstr[rank,4];
	If[rankDroped===ToExpression["{A,L,L}"],{ctr=4;rankDroped={0}},ctr=1];

	Do[

	Switch[rankDroped[[1]],
		0,
	fin={ToExpression["ee0"]->RedE0[pq1,pq2,pq3,pq4,pq5,
		     pq1pq2,pq2pq3,pq3pq4,pq4pq5,pq1pq5,
                     mq1,mq2,mq3,mq4,mq5]},
		1,
	fin=RedE1[pq1,pq2,pq3,pq4,pq5,
                     pq1pq2,pq2pq3,pq3pq4,pq4pq5,pq1pq5,
                     mq1,mq2,mq3,mq4,mq5,Listable->True]
                   /.q[_,n_]*F_:>ToExpression["ee"<>ToString[n]]->F,
		2,
	coeficients=RedE2[pq1,pq2,pq3,pq4,pq5,
                     pq1pq2,pq2pq3,pq3pq4,pq4pq5,pq1pq5,
                     mq1,mq2,mq3,mq4,mq5,Listable->True];
	eecoef=ToExpression["{ee00,ee11,ee12,ee13,ee14,ee22,ee23,ee24,ee33,ee34,ee44}"];
        fin=CreateRules[eecoef,coeficients],
		3,
	coeficients=RedE3[pq1,pq2,pq3,pq4,pq5,
                     pq1pq2,pq2pq3,pq3pq4,pq4pq5,pq1pq5,
                     mq1,mq2,mq3,mq4,mq5,Listable->True];
	eecoef=ToExpression["{ee001,ee002,ee003,ee004,
	                      ee111, ee112, ee113, ee114, ee122, ee123, ee124, ee133, 
	                      ee134, ee144, ee222, ee223, ee224, ee233, ee234, ee244, 
                              ee333, ee334, ee344, ee444}"];
	fin=CreateRules[eecoef,coeficients] ];

	rankDroped={i};
	result=Join[result,fin],{i,1,ctr}];
	
	
  Return[ Flatten[result] ]]




(*---------------------------------------*)
(*       Reduction to 4pt function       *)
(*---------------------------------------*)
I4[r_]:=Module[{boxes,result},

	boxes={D4[p1^2 + 2*p1*p5 + p5^2, p2^2, p3^2, p4^2, 
  	          p1^2 + 2*p1*p2 + p2^2 + 2*p1*p5 + 2*p2*p5 + p5^2, 
                  p2^2 + 2*p2*p3 + p3^2, 
                  m[5]^2, m[2]^2, m[3]^2, m[4]^2], 
	       D4[p1^2 + 2*p1*p2 + p2^2, p3^2, p4^2, p5^2, 
                  p1^2 + 2*p1*p2 + p2^2 + 2*p1*p3 + 2*p2*p3 + p3^2, 
                  p3^2 + 2*p3*p4 + p4^2, 
                  m[1]^2, m[3]^2, m[4]^2, m[5]^2], 
               D4[p1^2, p2^2 + 2*p2*p3 + p3^2, p4^2, p5^2, 
                  p1^2 + 2*p1*p2 + p2^2 + 2*p1*p3 + 2*p2*p3 + p3^2, 
                  p2^2 + 2*p2*p3 + p3^2 + 2*p2*p4 + 2*p3*p4 + p4^2, 
                  m[1]^2, m[2]^2, m[4]^2, m[5]^2], 
               D4[p1^2, p2^2, p3^2 + 2*p3*p4 + p4^2, p5^2, 
                  p1^2 + 2*p1*p2 + p2^2, 
                  p2^2 + 2*p2*p3 + p3^2 + 2*p2*p4 + 2*p3*p4 + p4^2, 
                  m[1]^2, m[2]^2, m[3]^2, m[5]^2], 
               D4[p1^2, p2^2, p3^2, p4^2 + 2*p4*p5 + p5^2, 
                  p1^2 + 2*p1*p2 + p2^2, 
                  p2^2 + 2*p2*p3 + p3^2, 
                  m[1]^2, m[2]^2, m[3]^2, m[4]^2]};

	result=boxes[[r]];

 Return[result]]

(*---------------------------------------*)
(*       Reduction to 4pt function       *)
(*---------------------------------------*)
I3[Sequence[n__]]:=Module[{rules,TT,result},

     rules={
	TT[{1, 2}] -> 
	C3[p1^2 + 2*p1*p2 + p2^2 + 2*p1*p5 + 2*p2*p5 + p5^2, p3^2, 
	   p1^2 + 2*p1*p2 + p2^2 + 2*p1*p3 + 2*p2*p3 + p3^2 + 2*p1*p5 + 2*p2*p5 + 2*p3*p5 + p5^2, 
           m[5]^2, m[3]^2, m[4]^2], 
 	TT[{1, 3}] -> 
        C3[p1^2 + 2*p1*p5 + p5^2, p2^2 + 2*p2*p3 + p3^2, 
           p1^2 + 2*p1*p2 + p2^2 + 2*p1*p3 + 2*p2*p3 + p3^2 + 2*p1*p5 + 2*p2*p5 + 2*p3*p5 + p5^2, 
           m[5]^2, m[2]^2, m[4]^2], 
        TT[{1, 4}] -> 
        C3[p1^2 + 2*p1*p5 + p5^2, p2^2, 
           p1^2 + 2*p1*p2 + p2^2 + 2*p1*p5 + 2*p2*p5 + p5^2, 
           m[5]^2, m[2]^2, m[3]^2], 
        TT[{1, 5}] -> 
        C3[p1^2 + 2*p1*p4 + p4^2 + 2*p1*p5 + 2*p4*p5 + p5^2, p2^2, 
           p1^2 + 2*p1*p2 + p2^2 + 2*p1*p4 + 2*p2*p4 + p4^2 + 2*p1*p5 + 2*p2*p5 + 2*p4*p5 + p5^2, 
           m[4]^2, m[2]^2, m[3]^2], 
        TT[{2, 3}] -> 
        C3[p1^2 + 2*p1*p2 + p2^2 + 2*p1*p3 + 2*p2*p3 + p3^2, p4^2, 
           p1^2 + 2*p1*p2 + p2^2 + 2*p1*p3 + 2*p2*p3 + p3^2 + 2*p1*p4 + 2*p2*p4 + 2*p3*p4 + p4^2, 
           m[1]^2, m[4]^2, m[5]^2], 
        TT[{2, 4}] -> 
        C3[p1^2 + 2*p1*p2 + p2^2, p3^2 + 2*p3*p4 + p4^2, 
           p1^2 + 2*p1*p2 + p2^2 + 2*p1*p3 + 2*p2*p3 + p3^2 + 2*p1*p4 + 2*p2*p4 + 2*p3*p4 + p4^2, 
           m[1]^2, m[3]^2, m[5]^2], 
        TT[{2, 5}] -> 
        C3[p1^2 + 2*p1*p2 + p2^2, p3^2, 
           p1^2 + 2*p1*p2 + p2^2 + 2*p1*p3 + 2*p2*p3 + p3^2, 
           m[1]^2, m[3]^2, m[4]^2], 
        TT[{3, 4}] -> 
        C3[p1^2, p2^2 + 2*p2*p3 + p3^2 + 2*p2*p4 + 2*p3*p4 + p4^2, 
           p1^2 + 2*p1*p2 + p2^2 + 2*p1*p3 + 2*p2*p3 + p3^2 + 2*p1*p4 + 2*p2*p4 + 2*p3*p4 + p4^2, 
           m[1]^2, m[2]^2, m[5]^2], 
        TT[{3, 5}] -> 
        C3[p1^2, p2^2 + 2*p2*p3 + p3^2, 
           p1^2 + 2*p1*p2 + p2^2 + 2*p1*p3 + 2*p2*p3 + p3^2, 
           m[1]^2, m[2]^2, m[4]^2], 
        TT[{4, 5}] -> 
        C3[p1^2, p2^2, 
           p1^2 + 2*p1*p2 + p2^2, 
           m[1]^2, m[2]^2, m[3]^2]};


	result=TT[Sort[List[n]]] /.rules;

 
 Return[result]]

(*---------------------------------------*)
(*       Reduction to 2pt function       *)
(*---------------------------------------*)
I2[Sequence[n__]]:=Module[{rules,TT,result},

     rules={
	TT[{1, 2, 3}] -> 
	B2[p1^2 + 2*p1*p2 + p2^2 + 2*p1*p3 + 2*p2*p3 + p3^2 + 2*p1*p5 + 2*p2*p5 + 2*p3*p5 + p5^2, 
           m[5]^2, m[4]^2], 
	TT[{1, 2, 4}] -> 
	B2[p1^2 + 2*p1*p2 + p2^2 + 2*p1*p5 + 2*p2*p5 + p5^2, 
           m[5]^2, m[3]^2], 
	TT[{1, 2, 5}] -> 
        B2[p1^2 + 2*p1*p2 + p2^2 + 2*p1*p4 + 2*p2*p4 + p4^2 + 2*p1*p5 + 2*p2*p5 + 2*p4*p5 + p5^2, 
	   m[4]^2, m[3]^2], 
        TT[{1, 3, 4}] -> 
        B2[p1^2 + 2*p1*p5 + p5^2, 
	   m[5]^2, m[2]^2], 
        TT[{1, 3, 5}] -> 
        B2[p1^2 + 2*p1*p4 + p4^2 + 2*p1*p5 + 2*p4*p5 + p5^2, 
           m[4]^2, m[2]^2], 
        TT[{1, 4, 5}] -> 
        B2[p1^2 + 2*p1*p3 + p3^2 + 2*p1*p4 + 2*p3*p4 + p4^2 + 2*p1*p5 + 2*p3*p5 + 2*p4*p5 + p5^2, 
           m[3]^2, m[2]^2], 
        TT[{2, 3, 4}] -> 
        B2[p1^2 + 2*p1*p2 + p2^2 + 2*p1*p3 + 2*p2*p3 + p3^2 + 2*p1*p4 + 2*p2*p4 + 2*p3*p4 + p4^2, 
           m[1]^2, m[5]^2], 
        TT[{2, 3, 5}] -> 
        B2[p1^2 + 2*p1*p2 + p2^2 + 2*p1*p3 + 2*p2*p3 + p3^2, 
           m[1]^2, m[4]^2], 
        TT[{2, 4, 5}] -> 
        B2[p1^2 + 2*p1*p2 + p2^2, 
           m[1]^2, m[3]^2], 
        TT[{3, 4, 5}] -> 
        B2[p1^2, 
           m[1]^2, m[2]^2]};

  
	result=TT[Sort[List[n]]] /.rules;


 Return[result]]



End[]


EndPackage[]

