############################### PROCEDURES #############################################

##################
get_explist := proc(mon, vars)

local explist,i:

explist := []:
for i from 1 to nops(vars) do
  explist := [op(explist), degree(mon, vars[i])]:
od:
explist:
end proc:

#########################
get_supp := proc(f, vars)

local mons, supp,i:

coeffs(f, vars, 'mons'):
mons := [mons]:
supp := []:
for i from 1 to nops(mons) do
   supp := [op(supp), get_explist(mons[i], vars)]:
od:
supp:
end proc:

##############################
#colorfacets := proc(cl)
#
#local ctest, ctest1, f, fctclrs, clrassign,fail,c,d:
#
#   fctclrs := []:
#   f := 1: fail := false:
#   while (f <= nops(div) and not(fail)) do
#      fctclrs := [op(fctclrs), []]:
#      c := 0: clrassign := false:
#      while (c < nops(cl)) do
#        c := c + 1:
#        ctest := true: ctest1 := true: d := 1:
#        while (d<=4 and ctest) do
#          if (convert(cl[c][d],set) = convert(div[f][d], set)) then
#                ctest := false : fi:
#          if (`subset`(convert(cl[c][d],set), convert(div[f][d],set))) then
#                ctest1 := false:
#          fi:
#          d := d + 1:
#        od:
#        if (ctest1) then 
#           if (not(clrassign)) then fctclrs[f] := [[c]]: 
#           else fctclrs[f] := [op(fctclrs[f]), [c]]: 
#           fi: clrassign := true: 
#         fi:
#        if (ctest and not(clrassign)) then fctclrs[f] := [op(fctclrs[f]), c]:
#        fi:
#       od:
#       if (fctclrs[f] = []) then fctclrs[f] := [5]: else if (not(clrassign)) then fctclrs[f] := [fctclrs[f]]: fi: fi:
#       f := f+1:
#   od:
#
#   return(fctclrs);
#end proc:
#
##################choose k-th color from the list of colors
#
#choosecolor:=proc(k,clrlist)
# local i,newclrlist:
# newclrlist:=[]:
# for i from 1 to nops(clrlist) do
#  if nops(clrlist[i])>1 then
#   newclrlist:=[op(newclrlist),[clrlist[i][k]]]:
#  else
#   newclrlist:=[op(newclrlist),[clrlist[i][1]]]:
#  fi:
# od:
# return(newclrlist):
#end proc:
#
              
##################### vertex decomposition of vert

vert_dec:=proc(vert)

local vert_coeff,l:

vert_coeff:=coeftayl(f,[x,y,z]=[0,0,0],vert): 
l:=sscanf(convert(vert_coeff,string),"%[a]%d%[*]%[b]%d%[*]%[c]%d%[*]%[d]%d"):

return([l[2],l[5],l[8],l[11]]):

end proc:

################# returns the support of column k of matrix M_f
c_supp:=proc(M_f,k)

local i,cs:

cs:=[]:
for i from 1 to 4 do 
  if M_f[i][k]=1 then cs:=[op(cs),i]: fi:
od:
return(convert(cs,set)):
end proc:


################ returns the restriction of the partition matrix M to the facet number ft 
restriction:=proc(M,ft)
local restr,s,t,r:

               restr:=[]:
               for s from 1 to 4 do
                 r:=[]:
                 for t from 1 to nops(div[ft][s]) do
                   r:=[op(r),M[s][div[ft][s][t]+1]]:
                 od:
                 restr:=[op(restr),r]:
               od:
return(restr):
end proc:

################ returns the color of a coloring matrix with a zero column, otherwise returns {5} 
 facecolor:=proc(M_f)
  local fctcol,i,j,r,s, four:
 
 fctcol:={5}: i:=1:
 four:=permute([1,2,3,4]):
 s:=four[1]:
  
  while nops(c_supp(M_f,s[1]))>0 and i<14 do i:=i+6: s:=four[i]: od:
  if nops(c_supp(M_f,s[1]))=0 then   
    if nops(c_supp(M_f,s[2]))>1 then  
      if nops(c_supp(M_f,s[3]))>1 then  
        if nops(c_supp(M_f,s[4]))>1 then 
          if nops(`union`(c_supp(M_f,s[2]),c_supp(M_f,s[3])))>2 then
 	   if nops(`union`(c_supp(M_f,s[2]),c_supp(M_f,s[4])))>2 then
 	     if nops(`union`(c_supp(M_f,s[3]),c_supp(M_f,s[4])))>2 then
 	       fctcol:={s[1]}: 
 	     else fctcol:={s[1],s[3],s[4]}:
              fi: 
 	   else fctcol:={s[1],s[2],s[4]}:
            fi:
          else fctcol:={s[1],s[2],s[3]}:
          fi:
        else 
          if nops(`union`(c_supp(M_f,s[2]),c_supp(M_f,s[4])))<3 then fctcol:={s[1],s[2],s[4]}: 
          else
            if nops(`union`(c_supp(M_f,s[3]),c_supp(M_f,s[4])))<3 then fctcol:={s[1],s[3],s[4]}: 
            else fctcol:={s[1],s[4]}: 
            fi:
          fi:
        fi:
      else
        if nops(`union`(c_supp(M_f,s[2]),c_supp(M_f,s[3])))<3 then fctcol:={s[1],s[2],s[3]}: 
        else
          if nops(`union`(c_supp(M_f,s[3]),c_supp(M_f,s[4])))<3 then fctcol:={s[1],s[3],s[4]}: 
          else fctcol:={s[1],s[3]}: 
          fi:
        fi:  
      fi:
    else
      if nops(`union`(c_supp(M_f,s[2]),c_supp(M_f,s[4])))<3 then fctcol:={s[1],s[2],s[4]}: 
      else
        if nops(`union`(c_supp(M_f,s[2]),c_supp(M_f,s[3])))<3 then fctcol:={s[1],s[2],s[3]}: 
        else fctcol:={s[1],s[2]}: 
        fi:
      fi:
    fi:
  fi:
 
 return(fctcol):
 end proc:

################compute combinatorial degree with partition matrix M

cdeg:=proc(M)
local facetcolor,edgecolor,vertexcolor,cd,f,v,sl,edgcol,fctcol,facetcolorlist, edgecolorlist, vertexcolorlist:

### returns the color of the facet number f if the first colomn of its coloring matrix is zero, otherwise returns {5} 
facetcolor:=proc(f)
 local i,j,rest,M_fct:

M_fct:=[[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]]:
rest:=restriction(M,f):
for i from 1 to 4 do
  for j from 1 to nops(rest[i]) do                                       ## makes the coloring matrix M_fct of the facet
    M_fct[i][rest[i][j]+1]:=1:
  od:
od:
return(facecolor(M_fct)):
end proc:

### returns the coloring of a vertex v
vertexcolor:=proc(v)
  local vtcolor,vt,i,M_v:

M_v:=[[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]]:
vt:=vert_dec(v):
for i from 1 to 4 do
  M_v[i][M[i][vt[i]+1]+1]:=1:                                                ## makes the coloring matrix M_v of the vertex
od: 

vtcolor:=[]:
for i from 1 to 4 do 
 if nops(c_supp(M_v,i))<2 then vtcolor:=[op(vtcolor),i]: 
 fi:
od:
return(convert(vtcolor,set)):
end proc:

##### returns the color of the edge (vert1, vert2) 
edgecolor:=proc(vert1,vert2)
 local i,vt1,vt2,M_edg:

M_edg:=[[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]]:
vt1:=vert_dec(vert1):
vt2:=vert_dec(vert2):
for i from 1 to 4 do
  M_edg[i][M[i][vt1[i]+1]+1]:=1:                                                ## makes the coloring matrix M_edg of the edge
  M_edg[i][M[i][vt2[i]+1]+1]:=1:
od: 

return(facecolor(M_edg)):
end proc:

##### here the procedure cdeg starts

#fprintf(dfile,"  M=%a, ",M):

cd:=0:

facetcolorlist:=[]:
edgecolorlist:=[]:
vertexcolorlist:=[]:

for f from 1 to nops(newfaces) do 
 fctcol:=facetcolor(f):
 facetcolorlist:=[op(facetcolorlist),fctcol]:             #### remember the color of the facet 
 if fctcol={1} then     #fprintf(dfile,"normal=%a\n",fcts(f)):                              
                        #if (fcts[f][1])^2+(fcts[f][2])^2+(fcts[f][3])^2=1 then fprintf(dfile,"normal=%a",fcts[f]): fi:
  for v from 1 to nops(newfaces[f]) do
   if v<nops(newfaces[f]) then sl:=v+1: else sl:=1:
   fi:
   edgcol:=edgecolor(newfaces[f][v],newfaces[f][sl]):
   if edgcol={1,2} then
    edgecolorlist:=[op(edgecolorlist),[newfaces[f][v],newfaces[f][sl],'blue']]: #### remember the edge and its color
    if vertexcolor(newfaces[f][v])={1,2,3} then
       vertexcolorlist:=[op(vertexcolorlist),[newfaces[f][v],'green']]:  #### remember the vertex and its color
       cd:=cd+1: #fprintf(dfile,"+"):
    fi:
    if vertexcolor(newfaces[f][sl])={1,2,3} then
       vertexcolorlist:=[op(vertexcolorlist),[newfaces[f][sl],'green']]:  #### remember the vertex and its color
       cd:=cd-1: #fprintf(dfile,"-"):
    fi:
   fi:
   if edgcol={1,2,3} then
    edgecolorlist:=[op(edgecolorlist),[newfaces[f][v],newfaces[f][sl],'green']]: #### remember the edge and its color
    if vertexcolor(newfaces[f][v])={1,2} then
       vertexcolorlist:=[op(vertexcolorlist),[newfaces[f][v],'blue']]:  #### remember the vertex and its color
       cd:=cd-1: #fprintf(dfile,"-"):
    fi:
    if vertexcolor(newfaces[f][sl])={1,2} then
       vertexcolorlist:=[op(vertexcolorlist),[newfaces[f][sl],'blue']]:  #### remember the vertex and its color
       cd:=cd+1: #fprintf(dfile,"+"):
    fi:
   fi:
  od:
 fi:
od:
#fprintf(dfile," cdeg=%d\n",cd):
return([cd,facetcolorlist,edgecolorlist,vertexcolorlist]):
end proc:



###################################################plot polytope with multicolored facets
plotpoly := proc(list)
   local i,j,k,l,m,p,s, cmass, newface,colors,newv:
   p:= []:
   colors := ['red', 'blue', 'green', 'yellow', 'gray']:
   for i from 1 to nops(list) do
       s := list[i]:
       p := [op(p), polygon(newfaces[i],color=colors[s[1]],linestyle=1,thickness=1)]:
       if (nops(s) > 1) then
          cmass := [0,0,0]: 
          for j from 1 to nops(newfaces[i]) do
            for k from 1 to 3 do 
            cmass[k] := cmass[k] + newfaces[i][j][k]/nops(newfaces[i]): 
            od:
          od:
       for j from 2 to nops(s) do
          newface := []:
          for k from 1 to nops(newfaces[i]) do 
             newv := [0,0,0]:
             for l from 1 to 3 do 
               newv[l] := cmass[l] + (nops(s)+1-j)/(nops(s))*(newfaces[i][k][l] - cmass[l]): 
             od:
             newface := [op(newface), newv]:
          od:
          p := [op(p), polygon(newface, color=colors[s[j]], linestyle=1, thickness=1)]:
       od:
       fi:
   od:
   return(p):
#   with(plottools):
#   plots[display](p);
end proc:

#################################### plot polytope with colored flags that give the cdeg
#plotpoly:=proc(clrlist) 
#
#  local p,e,v,q,i:
#  
#  p:=[]: e:=[]: v:=[]:
#  for i from 1 to nops(newfaces) do
#    p:=[op(p), polygon(newfaces[i], color=clrlist[1][i], linestyle=1, thickness=1)]:
#  od:
#  for i from 1 to nops(clrlist[2]) do
#    e:=[op(e), line(clrlist[2][i][1], clrlist[2][i][2], color=clrlist[2][i][3],linestyle=1, thickness=2)]:
#  od:
#  for i from 1 to nops(clrlist[3]) do
#    v:=[op(v), tetrahedron(clrlist[3][i][1], .5, color=clrlist[3][i][2])]:
#  od:
#  q:=convert([op(p),op(e),op(v)],'list'):
#   with(plottools):
#   plots[display](q);
#end proc:
#
########################################   

tetr:=proc(v):
  return(polygon([v[1],v[2],v[3]],color='gray', linestyle=1, thickness=1), 
         polygon([v[1],v[2],v[4]],color='gray', linestyle=1, thickness=1),
	 polygon([v[2],v[3],v[4]],color='gray', linestyle=1, thickness=1),
	 polygon([v[3],v[1],v[4]],color='gray', linestyle=1, thickness=1)):
end proc:


######################################### plot 3 simplices with the partition and their sum with flags
plot4poly:=proc(M,clrlist)

   local v0,v1,v2,v3,p0,p1,p2,p3,c,c0,c1,c2,c3,i,p,e,v,q:
   
   v0:=f0supp: v1:=f1supp: v2:=f2supp: v3:=f3supp:
   for i from 1 to 4 do
     v0[i][1]:=v0[i][1]-5: v1[i][1]:=v1[i][1]+5: v2[i][1]:=v2[i][1]+15: v3[i][1]:=v3[i][1]+25:
     v0[i][3]:=v0[i][3]-10: v1[i][3]:=v1[i][3]-10: v2[i][3]:=v2[i][3]-10: v3[i][3]:=v3[i][3]-10:
   od:  
   p0:=[tetr(v0)]: p1:=[tetr(v1)]: p2:=[tetr(v2)]: p3:=[tetr(v3)]: 
   
   c:=['red','blue','green','yellow']:

   c0:=[tetrahedron(v0[1],.3,color=c[M[1][1]+1]),
        tetrahedron(v0[2],.3,color=c[M[1][2]+1]),
	tetrahedron(v0[3],.3,color=c[M[1][3]+1]),
	tetrahedron(v0[4],.3,color=c[M[1][4]+1])]:
   
   c1:=[tetrahedron(v1[1],.3,color=c[M[2][1]+1]),
        tetrahedron(v1[2],.3,color=c[M[2][2]+1]),
	tetrahedron(v1[3],.3,color=c[M[2][3]+1]),
	tetrahedron(v1[4],.3,color=c[M[2][4]+1])]:
   
   c2:=[tetrahedron(v2[1],.3,color=c[M[3][1]+1]),
        tetrahedron(v2[2],.3,color=c[M[3][2]+1]),
	tetrahedron(v2[3],.3,color=c[M[3][3]+1]),
	tetrahedron(v2[4],.3,color=c[M[3][4]+1])]:
    
   c3:=[tetrahedron(v3[1],.3,color=c[M[4][1]+1]),
        tetrahedron(v3[2],.3,color=c[M[4][2]+1]),
	tetrahedron(v3[3],.3,color=c[M[4][3]+1]),
	tetrahedron(v3[4],.3,color=c[M[4][4]+1])]:
 
   e:=[]: v:=[]:
  p:=plotpoly(clrlist[1]): 
  for i from 1 to nops(clrlist[2]) do
    e:=[op(e), line(clrlist[2][i][1], clrlist[2][i][2], color=clrlist[2][i][3],linestyle=1, thickness=2)]:
  od:
  for i from 1 to nops(clrlist[3]) do
    v:=[op(v), tetrahedron(clrlist[3][i][1], .5, color=clrlist[3][i][2])]:
  od:
  q:=convert([op(p),op(e),op(v),op(p0),op(p1),op(p2),op(p3),op(c0),op(c1),op(c2),op(c3)],'list'):
  with(plottools):
  plots[display](q);
  
end proc:
   
#####################################################

######################################################


# 1) pick polynomilas (f0,f1,f2,f3)
# 2) get supports (f0supp,f1supp,f2supp,f3supp,supp) and create filename.poi and filename.site
# 3) run porta and ch (filename.poi.ieq, filename.fct)
# 4) read facet normals (fcts)
# 5) compute Minkowski sum decomposition of facets (div)
# 6) compute facets given by vertices (newfaces)
# 7) compute all partition matrices (Part_Mat)
# 8) compute cdegree for each partition matrix

with(linalg):
with(combinat, permute):

unassign('f0','f1','f2','f3','f0supp','f1supp','f2supp','f3supp','supp','fcts','div','newfaces','colorings','four','good','vertices','v','PartMat','face'):

four := permute([1,2,3,4]):
x := 'x':
y := 'y':
z := 'z':

########################## 1) pick polynomials
   Seed:=randomize();
   r := rand(0..9): 
   if type(numsystem,integer)=true then good := 1: else good:=0 fi:

if numsystem=-2 then 
f0:= a0*x^5*y^2*z^5+a1*x*y^2*z^9+a2*x^4*y^5*z^8+a3*x*y*z:
f1:= b0*x^3*y^2*z+b1*x^7*y^4+b2*x^4*y^8*z^5+b3*x^8*y^8*z:
f2:= c0*x*y^7*z^4+c1*x^6*y^5*z^6+c2*x^9*y^8*z^8+c3*x^4*y^7*z^6:
f3:= d0*x*y^4+d1*x^2*y^3*z^2+d2*x^3*y^3*z^9+d3*y*z^7:
iter:=numsystem:
fi:   
      
if numsystem=-1 then 

f0:= a0*x^8*y^5*z^9+a1*x^2*y*z^3+a2*y*z^6+a3*x^7*y^5*z^9:
f1:= b0*x*y^4*z^2+b1*x^9*y^8*z^3+b2*x^6*y^5*z^2+b3*x^3*y^4*z^5:
f2:= c0*x^8*y^5*z^9+c1*x^7*y*z^2+c2*z^8+c3*x^3*y*z^5:
f3:= d0*y*z^7+d1*x^3*y^2*z^3+d2*x*y^7*z^5+d3*y^9*z^4:

iter:=numsystem:
fi:
if numsystem=0 then 

f0:= a0*x^9*y^8*z^9+a1*y*z^5+a2*x^9*z^7+a3*x*y^5*z^7:
f1:= b0*x^9*z^9+b1*x^7*y*z^4+b2*x^3*y^8*z^9+b3*x*y*z^4:
f2:=c0*x^9*y^8*z^9+c1*y*z^5+c2*x^9*z^7+c3*x*y^5*z^7:
f3:=d0*x^9*z^9+d1*x^7*y*z^4+d2*x^3*y^8*z^9+d3*x*y*z^4:

iter:=numsystem:
fi:

if numsystem=1 then 

f0:= a0*x^9*y*z^5+a1*x^2*y^5*z^5+a2*x^2*y^2*z+a3*x^7*y^6*z^7:
f1:= b0*x*y^2*z^3+b1*x^9*y^7*z^5+b2*x^5*y*z^8+b3*x^6*y^9*z^2:
f2:= c0*x^9*y^9*z^5+c1*x^8*y^3*z^5+c2*x^2*y*z^9+c3*x^6*y^7*z^5:
f3:= d0*x^4*y^5+d1*x^7*y*z^2+d2*x^7*y^9*z^3+d3*x^7*y^8*z^7:

iter:=numsystem:
fi:

if numsystem=2 then 

f0:= a0*y^7*z^2+a1*x^5*y^2*z+a2*x^3*y^3*z^3+a3*x^7*y^6*z^5:
f1:= b0*x^6*y^9*z^2+b1*x*y*z+b2*x^9*y^6+b3*y^2*z:
f2:= c0*x*y^3*z^9+c1*x^6*y^9*z^5+c2*z^3+c3*x^3*y^3*z^2:
f3:= d0*x^9+d1*x*y^9*z^4+d2*x^9*y^2*z+d3*x^5*y^3*z^7:

iter:=numsystem:
fi:


if numsystem=6 then
f0:= a0*x^9*y^6*z^6+a1*x*y^7*z^4+a2*x^6*y^6*z^8+a3*x^2*y^9*z^3:
f1:= b0*x^8*y^8*z^2+b1*x*y^2+b2*x^3*y^6*z+b3*x^9*y^4*z^8:
f2:= c0*x^4*y^4+c1*x^7*y^6*z^9+c2*x^5*y^8*z^2+c3*x*y^3*z^3:
f3:= d0*x^2*y^5*z^6+d1*x^6*y^8*z+d2*x^7*y^9*z^4+d3*x^3*y^3*z^5:

iter:=numsystem:
fi:

if numsystem=11 then
f0:= a0*y^4*z^6+a1*x^6*y^4*z^6+a2*x^8*z^5+a3*x*y^4*z^7:
f1:= b0*x^4*y^2*z+b1*x^2*y^4*z^6+b2*x^2*y^9*z^8+b3*x*y^9*z:
f2:= c0*x^3*y*z^4+c1*x^6*y^3*z+c2*x^9*y^3*z^3+c3*z^4:
f3:= d0*x^4*y^2*z+d1*x^6*y^6*z^6+d2*x^8*y^8*z^8+d3*x^7*y^2*z^8:

iter:=numsystem:
fi:

if numsystem=28 then
f0:= a0*x*y+a1*x^6*y*z^4+a2*x^5*z^5+a3*x^3*y^5*z^4:
f1:= b0*x^7*y^2*z^4+b1*x^5*y^8*z^2+b2*x^8*y^8*z^4+b3*x^5*y^9:
f2:= c0*x^6*z^6+c1*x*y^6*z^6+c2*x^6*y^7*z^2+c3*x^5*y^7*z:
f3:= d0*x^7*y^9*z^2+d1*x*z^9+d2*x^5*y*z+d3*x^9*y^3*z^5:
iter:=numsystem:
fi:


if numsystem=100 then

f0:= a0*x^9+a1*x*y^9*z+a2*y*z^9+a3*x^2*y^2*z^2:
f1:= b0*y*z+b1*y^9+b2*z^9+b3*x*y^3*z^3:
f2:= c0+c1*x^9*z+c2*x*z^9+c3*x^3*y*z^3:
f3:= d0*x+d1*x^9*y+d2*y^9+d3*x^3*y^3*z:
iter:=numsystem:
fi:

fopen(dfile,APPEND):
if type(iter,integer)=true then
fprintf(dfile,"----------------------------------------- System %d-----------------\n",iter):
fi:

while (good=0) do
   exp0 := {}: exp1 := {}: exp2 := {}: exp3 := {}:
   for p from 0 to 3 do
     for i from 1 to 4 do
        exp||p := exp||p union {[r(), r(), r()]}:
     od:
   od:

  f0 := 0: for i from 0 to 3 do f0 := f0 + a||i *x^exp0[i+1][1]*y^exp0[i+1][2]*z^exp0[i+1][3]: od:
  f1 := 0: for i from 0 to 3 do f1 := f1 + b||i *x^exp1[i+1][1]*y^exp1[i+1][2]*z^exp1[i+1][3]: od:
  f2 := 0: for i from 0 to 3 do f2 := f2 + c||i *x^exp2[i+1][1]*y^exp2[i+1][2]*z^exp2[i+1][3]: od:
  f3 := 0: for i from 0 to 3 do f3 := f3 + d||i *x^exp3[i+1][1]*y^exp3[i+1][2]*z^exp3[i+1][3]: od:


A0:=[op(exp0[2]-exp0[1]),op(exp0[3]-exp0[1]),op(exp0[4]-exp0[1])]:
A1:=[op(exp1[2]-exp1[1]),op(exp1[3]-exp1[1]),op(exp1[4]-exp1[1])]:
A2:=[op(exp2[2]-exp2[1]),op(exp2[3]-exp2[1]),op(exp2[4]-exp2[1])]:
A3:=[op(exp3[2]-exp3[1]),op(exp3[3]-exp3[1]),op(exp3[4]-exp3[1])]:

if rank(matrix(3,3,A0))<3 or rank(matrix(3,3,A1))<3 or rank(matrix(3,3,A2))<3 or rank(matrix(3,3,A3))<3
   then fprintf(dfile,"NOT BIG\n"):
else fprintf(dfile,"BIG\n"): good:=1: fi:

od:

   fprintf(dfile,"f0:= %a:\n",f0):
   fprintf(dfile,"f1:= %a:\n",f1):
   fprintf(dfile,"f2:= %a:\n",f2):
   fprintf(dfile,"f3:= %a:\n",f3):

########################################## 2) get supports and product support

   vars := [x,y,z]:
   f := collect(f0*f1*f2*f3, vars, 'distributed'):
   f0supp := [get_explist(coeff(f0, a0), vars), get_explist(coeff(f0, a1), vars), get_explist(coeff(f0, a2), vars), get_explist(coeff(f0, a3), vars)]:
   f1supp := [get_explist(coeff(f1, b0), vars), get_explist(coeff(f1, b1), vars), get_explist(coeff(f1, b2), vars), get_explist(coeff(f1, b3), vars)]:
   f2supp := [get_explist(coeff(f2, c0), vars), get_explist(coeff(f2, c1), vars), get_explist(coeff(f2, c2), vars), get_explist(coeff(f2, c3), vars)]:
   f3supp := [get_explist(coeff(f3, d0), vars), get_explist(coeff(f3, d1), vars), get_explist(coeff(f3, d2), vars), get_explist(coeff(f3, d3), vars)]:
   supp := get_supp(f, vars):

##### write supports to files 1,2

   file1 := open(cat(filename,".sites"), WRITE):
   file2 := open(cat(filename,".poi"), WRITE):
   fprintf(file1, "%d\n", nops(supp)):
   fprintf(file2, "DIM=3\n \n CONV_SECTION\n"):
   for i from 1 to nops(supp) do
      fprintf(file1, "%d %d %d\n", supp[i][1], supp[i][2], supp[i][3]):
      fprintf(file2, "%d %d %d\n", supp[i][1], supp[i][2], supp[i][3]):
   od:
   fprintf(file2, "END\n"):
   fclose(file1):
   fclose(file2):


############################################# 3) run porta and ch

   s1 := cat("ch ",filename):
   s2 := cat("traf -l ", filename, ".poi"):
   s3 := cat("traf -l ", filename, ".poi.ieq"):
   system(s1):
   system(s2):
   system(s3):

############################################# 3') read vertices, write Minkowski sum decomposition of vertices
#   file:=open(cat(filename,".poi.ieq.poi"),READ):
#   for i from 1 to 3 do
#     readline(file):
#   od:
#  vertices:=[]:
#  div_vert:=[]:
#  s:=readline(file):
#  while(s<>"") do:
#    l:=sscanf(s,"%[(]%d%[)]%d%d%d"):
#    vertices:=[op(vertices),[l[4],l[5],l[6]]]:
#    vert_coeff:=coeff(coeff(coeff(f,x^(l[4])),y^(l[5])),z^(l[6])):
#    l:=sscanf(convert(vert_coeff,string),"%[a]%d%[*]%[b]%d%[*]%[c]%d%[*]%[d]%d"): 
#    div_vert:=[op(div_vert),[l[2],l[5],l[8],l[11]]]:
#    s:=readline(file):
#  od:
#  fclose(file):
# printf("div_vert=%a",div_vert):
###########################################  4) read facet inequalities

   file1 := open(cat(filename,".poi.ieq"), READ):
   for i from 1 to 6 do
     readline(file1):
   od:
   s1 := readline(file1):
   v := [];
   fcts :=[];
   while (s1 <> "") do
      l := sscanf(s1, "%[(]%d%[)]%c%[^,<,=]%[<,=]%d"):
      v := [op(v), -l[7]]:
      i := 1:
      s2 := "":
      while(l[5][i] <> "") do
        if (l[5][i] <> " ") then
           s2 := cat(s2,l[5][i]):
        fi:
        i:=i+1:
      od:
      lt := sscanf(s2,"%[+,-]%[^,x]%[x]%d%[+,-]%[^,x]%[x]%d%[+,-]%[^,x]%[x]%d"):
      numterms := nops(lt)/4:
      newfct := [0,0,0]:
      for i from 0 to (numterms-1) do
         if (lt[4*i+1] = "+") then
            term := -1:
         else
            term := 1:
         fi:
         if (lt[4*i+2] <> "") then
            term := term*sscanf(lt[4*i+2], "%d")[1]:
         fi:
         newfct[lt[4*i+4]] := term:
      od:
      fcts := [op(fcts), newfct]:
      s1 := readline(file1):
   od:
   fclose(file1):

######################################## 5) compute minkowski sum decomposition of facets

   dot := proc(a,b) a[1]*b[1] + a[2]*b[2] + a[3]*b[3]: end proc:
   f0div := []:
   f1div := []:
   f2div := []:
   f3div := []:
   f0normals:=[]: f1normals:=[]: f2normals:=[]: f3normals:=[]:
   for i from 1 to nops(fcts) do
      f0min := infinity:
      f1min := infinity:
      f2min := infinity:
      f3min := infinity:
      for j from 1 to 4 do
         newf0 := dot(fcts[i], f0supp[j]):
         if (newf0 < f0min) then
            f0min := newf0:
         fi:
         newf1 := dot(fcts[i], f1supp[j]):
         if (newf1 < f1min) then
            f1min := newf1:
         fi:
         newf2 := dot(fcts[i], f2supp[j]):
         if (newf2 < f2min) then
            f2min := newf2:
         fi:
         newf3 := dot(fcts[i], f3supp[j]):
         if (newf3 < f3min) then
            f3min := newf3:
         fi:
      od:
      newf0div :=[]: newf1div:= []: newf2div:=[]: newf3div:=[]:
      for j from 1 to 4 do
         if (dot(fcts[i], f0supp[j]) = f0min) then
            newf0div := [op(newf0div), j-1]:
	    if (nops(newf0div)=3) then f0normals:=[op(f0normals),fcts[i]]:
            fi:
         fi:
         if (dot(fcts[i], f1supp[j]) = f1min) then
            newf1div := [op(newf1div), j-1]:
	    if (nops(newf1div)=3) then f1normals:=[op(f1normals),fcts[i]]: fi:
         fi:
         if (dot(fcts[i], f2supp[j]) = f2min) then
            newf2div := [op(newf2div), j-1]:
	    if (nops(newf2div)=3) then f2normals:=[op(f2normals),fcts[i]]: fi:
         fi:
         if (dot(fcts[i], f3supp[j]) = f3min) then
            newf3div := [op(newf3div), j-1]:
	    if (nops(newf3div)=3) then f3normals:=[op(f3normals),fcts[i]]: fi:
         fi:
      od:
      f0div := [op(f0div), newf0div]: f1div := [op(f1div), newf1div]:
      f2div := [op(f2div), newf2div]: f3div := [op(f3div), newf3div]:
   od:
   div := []:
   for i from 1 to nops(fcts) do
       div := [op(div), [f0div[i], f1div[i], f2div[i], f3div[i]]]:
       #printf("facet %d decomposes as: %a\n", i, div[i]);
   od:

   
################################################ 6) compute facets, compare to .fct file

for i from 1 to nops(fcts) do
       faces[i] := []:
od:

for j from 1 to nops(fcts) do
   for i from 1 to nops(supp) do
      if (dot(supp[i], fcts[j]) = v[j]) then
         faces[j] := [op(faces[j]), supp[i]]:
      fi:
   od:
od:


   fcs := convert(faces, 'list'):
   file1 := open(cat(filename, ".fct"), READ):
   numlines := fscanf(file1, "%d")[1]:
   oldfaces := []:
   for l from 1 to numlines do
       newface:=[]:
       numvert := fscanf(file1, "%d")[1]:
       for vert from 1 to numvert do
           newface := [op(newface), supp[fscanf(file1, "%d")[1]+1]]:
       od:
       oldfaces := [op(oldfaces), newface]:
   od:
   fclose(file1):
   newfaces :=[]:
   for i from 1 to nops(fcts) do
       s := convert(fcs[i], 'set'):
       j:=1:
       while(not(convert(oldfaces[j], 'set') subset s)) do
            j := j + 1:
       od:
       newfaces := [op(newfaces), oldfaces[j]]:
   od:
   [div,newfaces]:

############################################################# Check if the system is generic

generic:=true:
if (`intersect`(convert(f0normals,'set'),convert(f1normals,'set'),convert(f2normals,'set'),convert(f3normals,'set'))={}) then
  for face in newfaces do 
   if nops(face)>4 then generic:=false: break: 
   else 
    if nops(face)=4 then 
      if face[1][1]+face[3][1]<>face[2][1]+face[4][1] or 
         face[1][2]+face[3][2]<>face[2][2]+face[4][2] or 
	 face[1][3]+face[3][3]<>face[2][3]+face[4][3] 
       then generic:=false: break: 
      fi:
    fi:
   fi:
  od: 
else generic:=false:
fi:
fprintf(dfile,"generic=%a\n",generic):

########################################################### 7) find partition matrices PartMat 

with(combinat,permute):
four := permute([0,1,2,3]):
cycl:=[four[1],four[7],four[9],four[10]]:

PartMat:= []: clrlist:=[]: cdeglist:=[]:
cl := [[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]]:

for i in cycl do
    cl[1][1]:= i[1]: cl[1][2]:= i[2]: cl[1][3]:= i[3]: cl[1][4]:= i[4]:
    for j from 1 to 4 do
      cl[2][j]:= 1:
      for k from 1 to 4 do
        cl[3][k]:= 2:
        for m from 1 to 4 do
	   cl[4][m]:= 3: 
            perm:= permute(cl):  
	    for n from 1 to 24 do
            valid := true:
            permcl:=perm[n]:
            for ft from 1 to nops(fcts) do
               if (not(valid)) then break: fi:
               restr:=restriction(permcl,ft):            # printf("restr[%d]=%a\n",ft,restr):
               for l in four do
                  if (member(l[1], restr[1]) and member(l[2], restr[2]) and member(l[3], restr[3]) and member(l[4], restr[4]))
		   then  valid := false: break: 
                  fi:  #printf("false,%a,l=%a\n",cl,l):
               od:
            od: 
           if (valid) then
              mylist:=cdeg(permcl):
              if mylist[1]<>0 then 
               cdeglist:=[op(cdeglist),mylist[1]]:
               clrlist:=[op(clrlist),[mylist[2],mylist[3],mylist[4]]]: 
               PartMat:=[op(PartMat),permcl]:
               fprintf(dfile,"M=%a, cd=%d\n",permcl,mylist[1]):
              fi:
           fi:
           od:
	   cl[4][m]:=0:   
        od:
        cl[3][k]:=0:
      od:
      cl[2][j]:=0:
    od:
od:

################################################################ 8) compute cdegree for each partition matrix


#clrlist:=cdeg(PartMat[2]):

#plotpoly(clrlist):

#clrlist:=[]:
#for k from 1 to nops(PartMat) do
# fprintf(dfile,"k=%d\n",k):
# clrlist:=[op(clrlist),cdeg(PartMat[k])]:
#od:

fprintf(dfile,"--------------------------------------------------------------------\n"):

fclose(dfile):

###############################################################