% A.I.Davydychev and A.G.Grozin, hep-ph/0102???

input_case nil$

operator B;
% B(0) = V(0,1,1)
% B(1) = V(1,1,0)
% B(2) = V(1,0,1)
% B(3) = V(1,1,1)

% Gamma(x+n)/Gamma(x)
procedure G(x,n);
if n>0 then for i:=0:(n-1) product (x+i)
else if n<0 then 1/(for i:=1:(-n) product (x-i))
else 1$

procedure V(n0,n1,n2);
if max(n0,n1)<=0 or max(n0,n2)<=0 or max(n1,n2)<=0 then 0
else if n0=0 then
    q2^(2-n1-n2)*(G(-d/2+2,n1+n2-2)*G(d/2-1,1-n1)*G(d/2-1,1-n2))
    /(G(d-2,2-n1-n2)*factorial(n1-1)*factorial(n2-1))*B(0)
else if n1=0 then
    (G(-d+3,n0+2*n2-3)*G(d/2-1,1-n2))/(factorial(n0-1)*factorial(n2-1))
    /((-4)^(n2-1)*w1^(n0+2*n2-3))*B(2)
else if n2=0 then
    (G(-d+3,n0+2*n1-3)*G(d/2-1,1-n1))/(factorial(n0-1)*factorial(n1-1))
    /((-4)^(n1-1)*w^(n0+2*n1-3))*B(1)
else if n1<0 then
    if n2 neq 1 then % 1+ 2- (1)
        q2*V(n0,n1+1,n2)+((d-n0-n2-2*n1-1)*V(n0,n1+1,n2-1)+n0*w*V(n0+1,n1+1,n2-1))/(n2-1)
    else if n0 neq 1 then % n2=1: 0- (0)
        2/(n0-1)*(n1*(w*V(n0-1,n1+1,1)-V(n0-2,n1+1,1))+w1*V(n0-1,n1,2)-V(n0-2,n1,2))
    else % n0=n2=1: (2)
        (n1*(V(1,n1+1,0)-q2*V(1,n1+1,1))-w1*V(2,n1,1))/(d-n1-3)
else if n2<0 then
    if n1 neq 1 then % 2+ 1- (2)
        q2*V(n0,n1,n2+1)+((d-n0-n1-2*n2-1)*V(n0,n1-1,n2+1)+n0*w1*V(n0+1,n1-1,n2+1))/(n1-1)
    else if n0 neq 1 then % n1=1: 0- (0)
        2/(n0-1)*(n2*(w1*V(n0-1,1,n2+1)-V(n0-2,1,n2+1))+w*V(n0-1,2,n2)-V(n0-2,2,n2))
    else % n0=n1=1: (1)
        (n2*(V(1,0,n2+1)-q2*V(1,1,n2+1))-w*V(2,1,n2))/(d-n2-3)
else if n1>1 then % 1- (2)
    (V(n0,n1,n2-1)-(n0*w1*V(n0+1,n1-1,n2)+(d-n0-n1-2*n2+1)*V(n0,n1-1,n2))/(n1-1))/q2
else if n2>1 then % 2- (1)
    (V(n0,n1-1,n2)-(n0*w*V(n0+1,n1,n2-1)+(d-n0-n2-2*n1+1)*V(n0,n1,n2-1))/(n2-1))/q2
% n1=n2=1
else if n0=1 then B(3)
else if n0<0 then % 0+ (3)
    ((d-2*n0-4)*(w+w1)*V(n0+1,1,1)+(n0+1)*(q2+4*w*w1)/2*V(n0+2,1,1)
    +V(n0,2,0)-w*V(n0+1,2,0)+V(n0,0,2)-w1*V(n0+1,0,2))
    /(2*(d-n0-3))
else % n0>1: 0- (3)
    2/((n0-1)*(q2+4*w*w1))*
    (-(d-2*n0)*(w+w1)*V(n0-1,1,1)+2*(d-n0-1)*V(n0-2,1,1)
        +w*V(n0-1,2,0)-V(n0-2,2,0)+w1*V(n0-1,0,2)-V(n0-2,0,2))$

remember V;

end;

% d.v     (0) n0 0+ + 2 n1 1+ (0- - w) + 2 n2 2+ (0- - w1) = 0
procedure T0(n0,n1,n2);
n0*V(n0+1,n1,n2)
+2*n1*(V(n0-1,n1+1,n2)-w*V(n0,n1+1,n2))
+2*n2*(V(n0-1,n1,n2+1)-w1*V(n0,n1,n2+1))$

% d.l     (1) d-n0-n2-2n1 + w n0 0+ + n2 2+ (q2 - 1-) = 0
procedure T1(n0,n1,n2);
(d-n0-n2-2*n1)*V(n0,n1,n2)+w*n0*V(n0+1,n1,n2)
+n2*(q2*V(n0,n1,n2+1)-V(n0,n1-1,n2+1))$

% d.(l-q) (2) d-n0-n1-2n2 + w1 n0 0+ + n1 1+ (q2 - 2-) = 0
procedure T2(n0,n1,n2);
(d-n0-n1-2*n2)*V(n0,n1,n2)+w1*n0*V(n0+1,n1,n2)
+n1*(q2*V(n0,n1+1,n2)-V(n0,n1+1,n2-1))$

% (3): q2/2 (0) + w (2) + w1 (1) - 0- [(1) + (2)]
% (q2 + 4 w w1)/2 n0 0+ + (d-2n0-n1-n2+1)(w+w1) - n1 w1 - n2 w
% - [2(d-n0+1)-3(n1+n2)] 0- - n1 1+ 2- (w - 0-) - n2 2+ 1- (w1 - 0-) = 0
procedure T3(n0,n1,n2);
(q2+4*w*w1)/2*n0*V(n0+1,n1,n2)
+((d-2*n0-n1-n2+1)*(w+w1)-n1*w1-n2*w)*V(n0,n1,n2)
-(2*(d-n0+1)-3*(n1+n2))*V(n0-1,n1,n2)
-n1*(w*V(n0,n1+1,n2-1)-V(n0-1,n1+1,n2-1))
-n2*(w1*V(n0,n1-1,n2+1)-V(n0-1,n1-1,n2+1))$

m:=3$

for n0:=-m:m do for n1:=-m:m do for n2:=-m:m do
if T0(n0,n1,n2) neq 0 then <<write n0,n1,n2; rederr "T0";>>;

for n0:=-m:m do for n1:=-m:m do for n2:=-m:m do
if T1(n0,n1,n2) neq 0 then <<write n0,n1,n2; rederr "T1";>>;

for n0:=-m:m do for n1:=-m:m do for n2:=-m:m do
if T2(n0,n1,n2) neq 0 then <<write n0,n1,n2; rederr "T2";>>;

for n0:=-m:m do for n1:=-m:m do for n2:=-m:m do
if T3(n0,n1,n2) neq 0 then <<write n0,n1,n2; rederr "T3";>>;

end;