<< MB.m

<< AMBRE.m

int1 = 
        Fullintegral[{1}, {PR[k1, m1, n1] PR[k1 + p1, m2, n2] PR[k1 - p2, 
                m3, n3]}, {k1}];

invariants = {p1^2 -> M1^2, p2^2 -> M2^2, 
        p1*p2 -> 1/2 s -  1/2 M1^2 - 1/2 M2^2}

IntPart[1]

repr = SubLoop[integral];

Fauto[0]

Simplify[fupc] // Expand

fupc = % /. 
            (m1^2) X[1] X[2] - (M1^2) X[1] X[2] + (m2^2) X[1] X[2] -> 
              MM1 X[1] X[2] /. 
          (m1^2) X[1] X[3] - (M2^2) X[1] X[3] + (m3^2) X[1] X[3] -> 
            MM2 X[1] X[3] /. -(M1^2) X[2] X[3] + (m2^2) X[2] X[3] - 
            (M2^2) X[2] X[3] + (m3^2) X[2] X[3] - s X[2] X[3] -> 
          MM3 X[2] X[3]

repr = SubLoop[integral];

fin = BarnesLemma[repr, 1, Shifts -> True]

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

invariants = {p1^2 -> m^2, p2^2 -> m^2, 
          p1*p2 -> 1/2 s -  1/2 m^2 - 1/2 m^2};

IntPart[1]

repr = SubLoop[integral];

BarnesLemma[repr, 1]

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

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

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

ser = 
        MBexpand[integrals, Exp[eps*EulerGamma], {eps, 0, 2}];

MBintegrate[ser, {m -> 1, s -> -100}]