<< AMBREv1.2.m

B3lNP = PR[k1, 0, n1]*PR[k1 + p1, 0, n2]*PR[k1 + p1 + p2, 0, n3]*PR[k1 - k2, 0, n4]*PR[k2, 0, n5]*
PR[k2 + p1 + p2, 0, n6]*PR[p1 + p2 + p4 + k2 - k3, 0, n7]*PR[k2 - k3, 0, n8]*PR[k3, 0, n9]*PR[p4 - k3, 0, n10]

invariants = {p1^2 -> 0, p2^2 -> 0, p3^2 -> 0, p4^2 -> 0, p1*p2 -> 1/2*s, p3*p4 -> 1/2 s, p1*p3 -> 1/2*t, 
p2*p4 -> 1/2*t, p2*p3 -> 1/2*u, p1*p4 -> 1/2*u};

Fullintegral[{1}, {B3lNP}, {k1, k2, k3}]

(* integration over planar subloop k1 *)

IntPart[1]

loop1res = SubLoop[integral]

(* preparation of the 1-loop result *)

loop1coef = loop1res /. SubLoop1[c_, _] -> c /. StyleForm[arg_, _] -> arg

loop1PRs = loop1res /. SubLoop1[_, pr_] -> pr /. StyleForm[arg_, _] -> arg /. PR[q_, m_, n_] -> PR[q, m, -n]

(* merging with the non-planar part *)

integral = PR[k1, 0, n1]*PR[k1 - k2, 0, n4]*PR[k1 + p1, 0, n2]*PR[k1 + p1 + p2, 0, n3];

loop23PRs = B3lNP / integral

B7NP = loop1PRs*loop23PRs //. {PR[q_, m_, n1_] PR[q_, m_, n2_] -> PR[q, m, n1 + n2]}

(* powers must be a symbol *)

PRlist = List @@ (B7NP);

powrul = {};

Do[{powrul = Append[powrul, PRlist[[i]] /. PR[_, _, n_] -> Rule[ToExpression["a" <> ToString[i]], n]],
   PRlist[[i]] = PRlist[[i]] /. PR[q_, m_, _] -> PR[q, m, ToExpression["a" <> ToString[i]]]}, {i, 1, Length[PRlist]}];

powrul

B7NP = Times @@ PRlist

<< AMBREv3.0.m

invariants = {p1^2 -> 0, p2^2 -> 0, p3^2 -> 0, p4^2 -> 0, p1*p2 -> 1/2*s, p3*p4 -> 1/2 s, 
p1*p3 -> 1/2*t, p2*p4 -> 1/2*t, p2*p3 -> 1/2*u, p1*p4 -> 1/2*u};

SimplifyFpoly[{s + t + u -> 0}];

(* integration over non-planar subgraph k2,k3 *)

NPrepr = MBreprNP[{1}, 
{PR[k2, 0, a1]*PR[k2 - k3, 0, a2]*PR[k3, 0, a3]*PR[k2 + p1, 0, a4]*PR[k2 + p1 + p2, 0, a5]*
PR[- k3 + p4, 0, a6]*PR[k2 - k3 + p1 + p2 + p4, 0, a7]}, {k2, k3}]

(* shifting of z's (k1-representation is 3-dim) *)

NPrepr = NPrepr /. {z1 -> zz1, z2 -> zz2, z3 -> zz3, z4 -> zz4} /. powrul /. {zz1 -> z4, zz2 -> z5, zz3 -> z6, zz4 -> z7}

(* finalizing the output *)

repr = loop1coef*NPrepr

finres = repr[[1]] /. {n1 -> 1, n2 -> 1, n3 -> 1, n4 -> 1, n5 -> 1, n6 -> 1, n7 -> 1, n8 -> 1, n9 -> 1, n10 -> 1}

<< MB.m

<< MBresolve.m

step1 = MBresolve[finres, eps];

step2 = MBexpand[step1, Exp[2*eps*EulerGamma], {eps, 0, 0}];
MBanalytic = MBmerge[step2];

MBintegrate[MBanalytic, {s -> -1, t -> -1, u -> -1}]