Quit[]

<< MB.m

<< AMBRE.m

Fullintegral[{1}, {PR[k1 + k2, m, n1] 
            PR[k2, 0, n2] PR[k1 + p1, 0, n3] }, {k2, k1}];

invariants = {p1^2 -> m^2};

IntPart[1]

SubLoop[integral]

IntPart[2]

repr = SubLoop[integral];

BarnesLemma[repr, 1]

fin = BarnesLemma[repr, 2]

se3l1mOS = fin /. n1 -> 1 /. n2 -> 1 /. 
        n3 -> 1

rules = MBoptimizedRules[se3l1mOS, eps -> 0, {}, {eps}]

integrals = MBcontinue[se3l1mOS, eps -> 0, rules]

ser3 = 
        MBexpand[{integrals}, 
          Exp[2 eps EulerGamma], {eps, 0, 4}];

num1 = MBintegrate[ser3, {m -> 1}]

SE3l1mOnSh = (+ep^(-2)*(-1/2) + 
                ep^(-1)*(-5/4) + 
                ep*(55/16 - 5*zt[2] - 4*zt[3]) + 
                ep^2*(949/32 - 11/2*zt[2] - 10*zt[3] - 32*zt[4]) + 
                ep^3*(8575/64 + 55/4*zt[2] - 16*zt[2]*zt[3] - 11*zt[3] - 
                      80*zt[4] - 72*zt[5]) + 
                ep^4*(64189/128 + 949/8*zt[2] - 40*zt[2]*zt[3] + 
                      55/2*zt[3] - 16*zt[3]*zt[3] - 88*zt[4] - 180*zt[5] - 
                      458*zt[6]) - 11/8 - 2*zt[2]) /. 
            Li[x_, y_] -> PolyLog[x, y] /. zt[x_] -> Zeta[x] /. 
        ln[x_] -> Log[x]

res = 
        Normal[Series[-(SE3l1mOnSh) Gamma[
                      1 + ep]^2 (m^2)^(1 - 2 ep), {ep, 0, 4}]] /. 
            EulerGamma -> 0 // Expand;

Expand[%] /. m -> 1 // N

num1

Expand[res] /. m -> 3 // N