<<MB/MB.m;
<<MB/MBresolve.m

(*       auxiliary functions         *)

SortByDimension[l_List]:=Sort[l,Length[#1[[2,2]]]>Length[#2[[2,2]]]&];
CoeffEps[X_,n_]:=(X /.X[[1]]-> Simplify[Coefficient[X[[1]],ep,n]]);
MBDimension[int_MBint]:=Length[int[[2,2]]];

(*   The end of derivation of MB representation for the planar massless double box diagram at 
p1^2=p2^2=p3^2=p4^2=0.

Notation:   S=-s=-(p1+p2)^2, T=-t=-(p1+p3)^2 *)

(*  This is MB representation for the box   *)

Box1[a1_,a2_,a3_,a4_] :=
((S^(2-a1-a2-a3-a4-ep-z1)) (T^z1) Gamma[a1+a2+a3+a4-2+ep+z1]Gamma[a2+z1] Gamma[a4+z1] Gamma[2-a1-a2-a4-ep-z1] Gamma[2-a2-a3-a4-ep-z1] Gamma[-z1] )/(Gamma[a1] Gamma[a2] Gamma[a3] Gamma[a4] Gamma[4-a1-a2-a3-a4-2 ep]);

( 1/S^(a4+a5+a6+a7+ep-2+z2+z3+z4)/(Gamma[a4] Gamma[a5] Gamma[a6] Gamma[a7] Gamma[4-a4-a5-a6-a7-2 ep])Gamma[a4+a5+a6+a7-2+ep+z2+z3+z4]Gamma[a7+z2+z3+z4] Gamma[a5+z4] Gamma[2-a4-a5-a7-ep-z3-z4] Gamma[2-a5-a6-a7-ep-z2-z4] Gamma[-z2] Gamma[-z3] Gamma[-z4] )Box1[a1-z2,a2,a3-z3,a8-z4];

Simplify[%]

(*   Changing variables    *)

%/.  {z2-> z2-z4,z3-> z3-z4};

% /. {z3-> z3+z1};

% /. z4-> z4+z1;

%/. z2->z2+z1;

Simplify[% /. z4->z4+z2+z3]

(*  In addition to the usual factor  (I Pi^(d/2))^2,
let us pull out the factor

(S^(4-a1-a2-a3-a4-a5-a6-a7-a8-2 ep)).

Let us turn to the variable x=T/S = t/s      *)

K2[a1_,a2_,a3_,a4_,a5_,a6_,a7_,a8_] := (  (x^z1) Gamma[-z1] Gamma[a2+z1] Gamma[2-a5-a6-a7-ep-z1-z2] Gamma[2-a1-a2-a8-ep+z2] Gamma[2-a4-a5-a7-ep-z1-z3] Gamma[2-a2-a3-a8-ep+z3] Gamma[a7+z1-z4] Gamma[-2+a4+a5+a6+a7+ep+z1-z4] Gamma[a8-z2-z3-z4] Gamma[-z1-z2-z3-z4] Gamma[-2+a1+a2+a3+a8+ep+z4] Gamma[z2+z4] Gamma[z3+z4] Gamma[a5+z1+z2+z3+z4])/(Gamma[a2] Gamma[a4] Gamma[a5] Gamma[a6] Gamma[a7] Gamma[4-a4-a5-a6-a7-2 ep] Gamma[4-a1-a2-a3-a8-2 ep+z1-z4] Gamma[a8-z1-z2-z3-z4] Gamma[a3+z2+z4] Gamma[a1+z3+z4])

(*  The 2-box with the powers of the propagators equal to one *)

B2=x K2[1,1,1,1,1,1,1,0]

B2rules=MBoptimizedRules[B2,ep->0,{},{ep}]

B2cont=MBcontinue[B2,ep->0,B2rules];

B2select=MBpreselect[B2cont,{ep,0,0}]

B2exp=Simplify[MBexpand[B2select,E^(2EulerGamma ep),{ep,0,0}]]

B2expS=SortByDimension[MBmerge[B2exp]]

Length[B2expS]

MBDimension /@ B2expS

B2expS[[1]]

Barnes2[B2expS[[1]],z4]

res01=(17/4)Zeta[4]

%//N

NIntegrate[Barnes2[B2expS[[1]],z4] [[1]]/(2Pi)  /.{z3->-(3/8)+I*y1},{y1,-Infinity,Infinity} ]

B2expS[[2]]

Barnes2[B2expS[[2]],z4]

res02=(17/4)Zeta[4]

%//N

NIntegrate[Barnes2[B2expS[[2]],z4] [[1]]/(2Pi)  /.{z2->-(5/16)+I*y1},{y1,-Infinity,Infinity} ]

B2expS[[3]]

Barnes1[B2expS[[3]],z3]

B2expS[[4]]

Barnes1[B2expS[[4]],z2]

B2expS[[5]]

CoeffEps[B2expS[[5]],-1]

CoeffEps[B2expS[[5]],-1] /. {Gamma[1-z4]->-z4 Gamma[-z4],Gamma[1+z4]->z4 Gamma[z4]}

res11=Barnes1[%,z4][[1]]

CoeffEps[B2expS[[5]],0]

res03=-(11/2)Zeta[4]

%//N

NIntegrate[CoeffEps[B2expS[[5]],0][[1]]/(2Pi)  /.{z4->(7/16)+I*y1},{y1,-Infinity,Infinity}  ]

B2expS[[6]]

B2expS[[7]]

F67=MBint[Simplify[(B2expS[[6]][[1]] /. z3->z2)+B2expS[[7]][[1]]],{{ep->0},{z2->-(1/2)}}]

res04=-((π^2)/(6ep^2))+(1/ep)(((π^2)/3)Log[x]+Zeta[3])+(-((π^2)/3)Log[x]^2-2Zeta[3]Log[x]+((π^4)/9));

res04/. {ep->0.3,x->0.1}

NIntegrate[F67[[1]]/(2Pi)  /.{ep->0.3,x->0.1,z2->-(1/2)+I*y1},{y1,-Infinity,Infinity}  ]

B2expS[[8]]

B2expS[[9]][[1]]

(* Collecting result for terms with trivial dependence on x  *)

restr=Apart[FullSimplify[res01+res02+res11/ep+res03+res04+B2expS[[9]][[1]]],ep]

(* Collecting terms with nontrivial dependence on x  *)

MBmerge[{Barnes1[B2expS[[3]],z3] ,Barnes1[B2expS[[4]],z2]   ,B2expS[[8]]}][[1]]

NIntegrate[%[[1]]/(2Pi)  /.{ep->0.3,x->0.1,z1->-(1/2)+I*y1},{y1,-Infinity,Infinity}  ]

resnontr=-((2 (π^2) Log[1+x])/ep)+(10/3) (π^2) Log[x] Log[1+x]-((2 ((Log[x])^2) Log[1+x])/ep)+(8/3) ((Log[x])^3) Log[1+x]-(π^2) ((Log[1+x])^2)-((Log[x])^2) ((Log[1+x])^2)-(20/3) (π^2) PolyLog[2,-x]-((4 Log[x] PolyLog[2,-x])/ep)-2 ((Log[x])^2) PolyLog[2,-x]+((4 PolyLog[3,-x])/ep)+24 Log[x] PolyLog[3,-x]-4 (-(1/2) ((PolyLog[2,-x])^2)-Log[1+x] PolyLog[3,-x])-44 PolyLog[4,-x]+4 Log[x] (-Log[1+x] PolyLog[2,-x]-2 PolyLog[1,2,-x])+4 Log[x] PolyLog[1,2,-x]-4 ((1/2) ((PolyLog[2,-x])^2)-2 PolyLog[2,2,-x])-4 PolyLog[2,2,-x]-4 Log[1+x] Zeta[3];

%/.{ep->0.3,x->0.1}

resnontr+restr;

Apart[%,ep]

(*             the result            *)

1/x((4/(ep^4))-((5 Log[x])/(ep^3))+((-5 (π^2)+4 ((Log[x])^2))/(2 (ep^2)))+(1/(6 ep))(33 (π^2) Log[x]+4 ((Log[x])^3)-12 (π^2) Log[1+x]-12 ((Log[x])^2) Log[1+x]-24 Log[x] PolyLog[2,-x]+24 PolyLog[3,-x]-130 Zeta[3])+(1/30) (-29 (π^4)-180 (π^2) ((Log[x])^2)-40 ((Log[x])^4)+100 (π^2) Log[x] Log[1+x]+80 ((Log[x])^3) Log[1+x]-30 (π^2) ((Log[1+x])^2)-30 ((Log[x])^2) ((Log[1+x])^2)-200 (π^2) PolyLog[2,-x]-60 ((Log[x])^2) PolyLog[2,-x]-120 Log[x] Log[1+x] PolyLog[2,-x]+720 Log[x] PolyLog[3,-x]+120 Log[1+x] PolyLog[3,-x]-1320 PolyLog[4,-x]-120 Log[x] PolyLog[1,2,-x]+120 PolyLog[2,2,-x]+880 Log[x] Zeta[3]-120 Log[1+x] Zeta[3]))