(* ::Package:: *)

(* IBP identities in opearator notation *)
ibp1 = (d-2a1-a2)*op[0,0] - 2a1*m^2*op[1,0] - a2*op[-1,1] + a2*(q2-m^2)*op[0,1];
ibp2 = (a2-a1)*op[0,0] - a1*(q2+m^2)*op[1,0] + a1*op[1,-1] - a2*op[-1,1] + a2*(q2-m^2)*op[0,1];


(* Linear combination of ibp1 and ibp2 where op[1,0] cancels *)
ibpLC = Collect[(q2+m^2)*ibp1 - 2m^2*ibp2,{op[0,1],op[1,0],op[-1,1],op[1,-1],op[0,0],q2},Simplify]


(* Derive identities for Fm[a1,a2] *)
Simplify[Solve[Expand[ibpLC*Fm[a1,a2]==0] /. {op[i1_,i2_]*Fm[a1,a2] -> Fm[a1+i1,a2+i2]},Fm[a1,1+a2]]/.{a2->a2-1}]
Simplify[Solve[Expand[ibp1*Fm[a1,a2]==0] /. {op[i1_,i2_]*Fm[a1,a2] -> Fm[a1+i1,a2+i2]},Fm[1+a1,a2]]/.{a1->a1-1}]


(* Together with the BC for scaleless integrals this are the identities we use for the reduction *)
ReduceA2 = {Fm[a1_,a2_] /; TrueQ[a1>=1 && a2>1] -> (-(-1+a2) (m^2-q2) Fm[-1+a1,a2]+((-3+3 a2-d) m^2+(-1+2 a1+a2-d) q2) Fm[a1,-1+a2]+2 a1 m^2 Fm[1+a1,-2+a2])/((-1+a2) (m^2-q2)^2)};
ShuffleA1 = {Fm[a1_,a2_] /; TrueQ[a1>1] -> (-a2 Fm[-2+a1,1+a2]+(2-2 a1-a2+d) Fm[-1+a1,a2]+a2 (-m^2+q2) Fm[-1+a1,1+a2])/(2 (-1+a1) m^2)};
BC = {Fm[a1_,a2_] /; TrueQ[a1\[Element]Integers && a1<=0] :> 0};
ReduceAll = Union[ShuffleA1,ReduceA2,BC];


(* Apply to Fm[1,2] *)
Simplify[Fm[1,2] /. ReduceA2]
Simplify[% /. ShuffleA1 /. BC]


Collect[Fm[3,2] //. ReduceAll,{Fm[1,0],Fm[1,1]},Simplify]


Collect[Fm[5,3] //. ReduceAll,{Fm[1,0],Fm[1,1]},Simplify]