R:=(1-r-s-t-u+(2/3)*(r*s+r*t+r*u+s*t+s*u+t*u)):

phi:=1/R:

vars:=[r,s,t,u]:



extract:=proc(a,b,c)
 local  F:

 if a<0 or b<0 or c<0 then RETURN(`values too low`): fi:  

  F:=taylor(phi,r,a+1):
  F:=coeff(F,r,a):
  F:=taylor(F,s,b+1):
  F:=coeff(F,s,b):
  F:=taylor(F,t,c+1):
  F:=coeff(F,t,c):
  F:=factor(normal(expand(F))):

end:



check:=proc(A,B,C)
local  a,b,c,gu,du,lu,i,j,k,l:

a:=A:  b:=B:  c:=C:

lu:=0:
for i from 0 to a do
 for j from i to b do
  for k from j to c do
    gu:=extract(i,j,k):
    gu:=subs(u=1-w,gu):
    gu:=expand(numer(gu)):

    for l from 0 to degree(gu,w) do
      du:=coeff(gu,w,l):
      if du<0 then
	lu:=lu+1: 
 	print(`RATS!  I found a negative coefficient`,du,`when`); 
	print(`i=`,i,`j=`,j,`and k=`,k):
      fi:
    od: 

  od:
 od:
od:

if lu=0 then print(`ALL COEFFICIENTS ARE NON-NEGATIVE!`): fi:

end: