with Standard_Integer_Numbers;           use Standard_Integer_Numbers;
with Standard_Integer_Vectors;           use Standard_Integer_Vectors;
with Standard_Integer_Norms;             use Standard_Integer_Norms;
with Standard_Integer_VecVecs;           use Standard_Integer_VecVecs;
with Standard_Integer_Matrices;          use Standard_Integer_Matrices;
with Standard_Integer_Linear_Solvers;    use Standard_Integer_Linear_Solvers;
with Standard_Integer32_Transformations;
 use Standard_Integer32_Transformations;
with Transforming_Integer32_Vector_Lists;
 use Transforming_Integer32_Vector_Lists;
with Integer32_Vectors_Utilities;        use Integer32_Vectors_Utilities;
with Integer_Support_Functions;          use Integer_Support_Functions;
with Integer_Face_Enumerators;           use Integer_Face_Enumerators;
with Lists_of_Vectors32_Utilities;       use Lists_of_Vectors32_Utilities;
with Arrays_of_Lists_Utilities;          use Arrays_of_Lists_Utilities;
--with Face_Enumerator_of_Sum;             use Face_Enumerator_of_Sum;

package body Volumes is

-- INVARIANT CONDITION :
--   In order to get p(v) > 0, the zero vector must be in the first list,
--   so shifting is necessary.
--   The shift vector must be equal to the maximal element in the list,
--   w.r.t. the graded lexicographic ordening.
--   In this way, the shift vector is unique, which allows to `mirror'
--   the same operations for the mixed homotopy continuation.

-- AUXILIARIES :

  function Create_Facet
             ( n : integer32; facet : Vector; pts : VecVec ) return VecVec is

    res : VecVec(1..n);
    cnt : integer32 := 0;

  begin
    for i in facet'range loop
      cnt := cnt+1;
      res(cnt) := pts(facet(i));
    end loop;
    return res;
  end Create_Facet;

  function Determinant ( vv : VecVec ) return integer32 is

    a : Matrix(vv'range,vv'range);

  begin
    for k in a'range(1) loop
      for l in a'range(2) loop
        a(k,l) := vv(k)(l);
      end loop;
    end loop;
   -- Upper_Triangulate(a);
    return Det(a);
  end Determinant;

-- VOLUME COMPUTATIONS :

  function Vol ( n : natural32; L : List ) return natural32 is

  -- DESCRIPTION :
  --   Computes the volume of the simplex spanned by the list
  --   and the origin.

    res : integer32;
    vv : VecVec(1..integer32(n));
    tmp : List;

  begin
    tmp := L;
    for i in vv'range loop
      vv(i) := Head_Of(tmp);
      tmp := Tail_Of(tmp);
    end loop;
    res := Determinant(vv);
    if res < 0
     then return natural32(-res);
     else return natural32(res);
    end if;
  end Vol;

  function Volume ( n,i : integer32; L : List; v : Vector; pv : integer32 )
                  return natural32 is

  -- DESCRIPTION :
  --   The points in l all belong to the same hyperplane:
  --    < x , v > - pv = 0;
  --   this procedure computes the volume of the polytope generated by
  --   the points in l and the origin.

    ll : constant integer32 := integer32(Length_Of(L));

  begin
    if ll < n then
      return 0;
    elsif ll = n then
      return Vol(natural32(n),L);
    else
      declare
        t : Transfo; 
        tl,rl : List;
        vl : integer32;
      begin
        if pv > 0 then
          t := Build_Transfo(v,i);
          tl := t*L;
          rl := Reduce(tl,i);
          Clear(t); Deep_Clear(tl);
         -- Remove_Duplicates(rl);
          if integer32(Length_Of(rl)) >= n-1
           then vl := pv*integer32(Volume(natural32(n)-1,rl));
           else vl := 0;
          end if;
          Deep_Clear(rl);
          return natural32(vl);
        else
          return 0;
        end if;
      end;
    end if;
  end Volume;

  function Sum ( n,m : integer32; L : List ) return natural32 is

  -- DESCRIPTION :
  --   Computes the volume of the polytope generated by the points in l;
  --   where n > 1 and n < m = Length_Of(l).

    res : natural32;
    fix : Link_to_Vector;
    pts : VecVec(1..m-1);
    nulvec : constant Vector(1..n) := (1..n => 0);
    vl,wl,vl_last : List;
    sh : boolean;

    procedure Compute_Facet ( facet : in vector; cont : out boolean ) is

      vv : constant VecVec := Create_Facet(n,facet,pts);
      pl : List;
      v : Link_to_Vector;
      i,pv,sup : integer32;

    begin
      v := Compute_Normal(vv);
      i := Pivot(v);
      if i <= v'last then
        pv := vv(vv'first)*v;
        if pv < 0 then
          for j in v'range loop
            v(j) := -v(j);
          end loop;
          pv := -pv;
        end if;
        if (pv > 0) and then not Is_In(vl,v) then
          sup := Maximal_Support(wl,v.all);
          if sup = pv then
            Append(vl,vl_last,v);
            pl := Face(wl,v.all,pv);
            res := res + Volume(n,i,pl,v.all,pv);
            Deep_Clear(pl);
          end if;
        end if;
      end if;
      cont := true;
    end Compute_Facet;

    procedure Compute_Facets is new Enumerate_Faces(Compute_Facet);
   
  begin
    res := 0;
    if not Is_In(l,nulvec)
     then fix := Graded_Max(l);
          wl := Shift(l,fix); sh := true;
     else wl := l;            sh := false;
    end if;
    Move_to_Front(wl,nulvec);
    pts := Shallow_Create(Tail_Of(wl));
    Compute_Facets(n-1,pts);
    Deep_Clear(vl);
    if sh
     then Deep_Clear(wl); Clear(fix);
    end if;
    return res;
  end Sum;

  function Volume ( n : natural32; L : List ) return natural32 is

    m : natural32;

  begin
    if n = 1 then
      declare
        min,max : integer32 := 0;
        d : constant Link_to_Vector := Head_Of(L);
      begin
        Min_Max(L,d'first,min,max);
        return natural32(max - min);
      end;
    else
      m := Length_Of(L);
      if m <= n
       then return 0;
       else return Sum(integer32(n),integer32(m),L);
      end if;
    end if;
  end Volume;

  function Volume ( n,i : integer32; L : List; v : Vector;
                    pv : integer32; tv : Tree_of_Vectors ) return natural32 is

  -- DESCRIPTION :
  --   The points in l all belong to the same hyperplane:
  --    < x , v > - pv = 0;
  --   this procedure computes the volume of the polytope generated by
  --   the points in l and the origin.

    ll : constant integer32 := integer32(Length_Of(L));

  begin
    if ll < n then
      return 0;
    elsif ll = n then
      return Vol(natural32(n),L);
    else
      declare
        t : Transfo; 
        tl,rl : List;
        vl : integer32;
      begin
        if pv > 0 then
          t := Build_Transfo(v,I);
          tl := t*L;
          rl := Reduce(tl,i);
          Clear(t); Deep_Clear(tl);
         -- Remove_Duplicates(rl);
          if integer32(Length_Of(rl)) >= n-1
           then vl := pv*integer32(Volume(natural32(n-1),rl,tv));
           else vl := 0;
          end if;
          Deep_Clear(rl);
	  return natural32(vl);
        else
          return 0;
        end if;
      end;
    end if;
  end Volume;

  function Sum ( n,m : integer32; L : List; tv : Tree_of_Vectors )
	       return natural32 is

  -- DESCRIPTION :
  --   Computes the volume of the polytope generated by the points in l;
  --   where n > 1 and n < m = Length_Of(l).
  --   The tree of degrees tv is not empty.

    res : natural32;
    fix : Link_to_Vector;
    nulvec : constant Vector(1..n) := (1..n => 0);
    wl : List;
    tmp : Tree_of_Vectors;
    sh : boolean;

  begin
    res := 0;
    if not Is_In(L,nulvec)
     then fix := Graded_Max(L);
          wl := Shift(L,fix); sh := true;
     else wl := L;            sh := false;
    end if;
    Move_to_Front(wl,nulvec);
    tmp := tv;
    while not Is_Null(tmp) loop
      declare
	nd : constant node := Head_Of(tmp);
	v : constant Link_to_Vector := nd.d;
	pv,i : integer32;
	pl : List;
      begin
	i := Pivot(v);
	pv := Maximal_Support(wl,v.all);
	pl := Face(wl,v.all,pv);
	if (nd.ltv = null) or else Is_Null(nd.ltv.all)
	 then res := res + Volume(n,i,pl,v.all,pv);
	 else res := res + Volume(n,i,pl,v.all,pv,nd.ltv.all);
        end if;
	Deep_Clear(pl);
      end;
      tmp := Tail_Of(tmp);
    end loop;
    if sh
     then Deep_Clear(wl); Clear(fix);
    end if;
    return res;
  end Sum;

  function Volume ( n : natural32; L : List; tv : Tree_of_Vectors ) 
                  return natural32 is

    m : natural32;

  begin
    if n = 1 then
      declare
        min,max : integer32 := 0;
        d : constant Link_to_Vector := Head_Of(L);
      begin
        Min_Max(L,d'first,min,max);
        return natural32(max - min);
      end;
    else
      m := Length_Of(L);
      if m <= n then
        return 0;
      elsif not Is_Null(tv) then
        return Sum(integer32(n),integer32(m),L,tv);
      else
        return Sum(integer32(n),integer32(m),L);
      end if;
    end if;
  end Volume;

  procedure Volume ( n,i : in integer32; L : in List; v : in Vector;
                     pv : in integer32; tv : in out Tree_of_Vectors;
                     vlm : out natural32 ) is

  -- DESCRIPTION :
  --   The points in l all belong to the same hyperplane:
  --    < x , v > - pv = 0;
  --   this procedure computes the volume of the polytope generated by
  --   the points in l and the origin.

    ll : constant integer32 := integer32(Length_Of(L));
    vl : natural32;

  begin
    if ll < n then
      vlm := 0;
    elsif ll = n then
      vl := Vol(natural32(n),L);
      if vl > 0 then
        declare
          nd : node;
        begin
          nd.d := new Vector'(v);
          Construct(nd,tv);
        end;
      end if;
      vlm := vl;
    else
      declare
        t : Transfo; 
        tl,rl : List;
      begin
        if pv > 0 then
          t := Build_Transfo(v,i);
          tl := t*L;
          rl := Reduce(tl,i);
          Clear(t); Deep_Clear(tl);
         -- Remove_Duplicates(rl);
          if integer32(Length_Of(rl)) >= n-1 then
            declare
              tmp : Tree_of_Vectors;
            begin
              Volume(natural32(n)-1,rl,tmp,vl);
              if vl = 0 then
                Clear(tmp);
              else
                declare
                  nd : node;
                begin
                  nd.d := new Vector'(v);
                  if not Is_Null(tmp)
                   then nd.ltv := new Tree_of_Vectors'(tmp);
                  end if;
                  Construct(nd,tv);
                end;
              end if;
            end;
            vlm := natural32(pv)*vl;
          else
            vlm := 0;
          end if;
          Deep_Clear(rl);
        else
          vlm := 0;
        end if;
      end;
    end if;
  end Volume;

  procedure Sum ( n,m : in integer32; L : in List;
		  tv : in out Tree_of_Vectors; vlm : out natural32 ) is

  -- DESCRIPTION :
  --   Computes the volume of the polytope generated by the points in l;
  --   where n > 1 and n < m = Length_Of(L).

    res : natural32;
    pts : VecVec(1..m-1);
    fix : Link_to_Vector;
    nulvec : constant Vector(1..n) := (1..n => 0);
    wl : List;
    sh : boolean;

    procedure Compute_Facet ( facet : in vector; cont : out boolean ) is

      vv : constant VecVec := Create_Facet(n,facet,pts);
      pl : List;
      v : Link_to_Vector;
      i,pv,sup : integer32;
      pvol : natural32;

    begin
      v := Compute_Normal(vv);
      i := Pivot(v);
      if i <= v'last then
        pv := vv(vv'first)*v;
        if pv < 0 then
          for j in v'range loop
            v(j) := -v(j);
          end loop;
          pv := -pv;
        end if;
        if (pv > 0) and then not Is_In(tv,v) then
          sup := Maximal_Support(wl,v.all);
          if sup = pv then
            pl := Face(wl,v.all,pv);
            Volume(n,i,pl,v.all,pv,tv,pvol);
            res := res + pvol;
            Deep_Clear(pl);
          end if;
        end if;
      end if;
      cont := true;
    end Compute_Facet;

    procedure Compute_Facets is new Enumerate_Faces(Compute_Facet);
   
  begin
    res := 0;
    if not Is_In(L,nulvec) then
      fix := Graded_Max(L);
      wl := Shift(L,fix); sh := true;
    else
      wl := L;            sh := false;
    end if;
    Move_to_Front(wl,nulvec);
    pts := Shallow_Create(Tail_Of(wl));
    Compute_Facets(n-1,pts);
    if sh
     then Deep_Clear(wl); Clear(fix);
    end if;
    vlm := res;
  end Sum;

  procedure Volume ( n : in natural32; L : in List; 
		     tv : in out Tree_of_Vectors; vlm : out natural32 ) is

    m : natural32;

  begin
    if n = 1 then
      declare
        min,max : integer32 := 0;
        d : constant Link_to_Vector := Head_Of(L);
      begin
        Min_Max(L,d'first,min,max);
        vlm := natural32(max - min);
      end;
    else
      m := Length_Of(L);
      if m <= n
       then vlm := 0;
       else Sum(integer32(n),integer32(m),L,tv,vlm);
      end if;
    end if;
  end Volume;

-- MIXED VOLUME COMPUTATIONS WITHOUT TREE OF USEFUL DIRECTIONS :

  function Two_Terms_Mixed_Vol ( n : natural32; al : Array_of_Lists )
			       return natural32 is

  -- DESCRIPTION :
  --   returns the mixed volume of the polytopes generated by the
  --   points in al, where the first polytope is generated by only
  --   two points.

    first,second : Link_to_Vector;
    L : List renames al(al'first);
    res : natural32;

  begin
    first := Head_Of(L);
    second := Head_Of(Tail_Of(L));
    declare
      d : Vector(first'range);
      piv : integer32 := 0;
    begin
      for i in d'range loop
	d(i) := first(i) - second(i);
	if (piv = 0) and then (d(i) /= 0)
	 then piv := i;
        end if;
      end loop;
      if piv = 0 then
        return 0;
      else 
        if d(piv) < 0
         then Min(d);
        end if;
        declare
          t : Transfo := Rotate(d,piv);
	  tr_al : Array_of_Lists(al'first..(al'last-1));
          degen : boolean := false;
        begin
          Apply(t,d);
          for i in tr_al'range loop
            tr_al(i) := Transform_and_Reduce(t,piv,al(i+1));
            Remove_Duplicates(tr_al(i));
            degen := (Length_Of(tr_al(i)) <= 1);
            exit when degen;
          end loop;
          Clear(t);
          if not degen
           then res := natural32(d(piv))*Mixed_Volume(n-1,tr_al);
           else res := 0;
          end if;
          Deep_Clear(tr_al);
        end;
      end if;
    end;
    return res;
  end Two_Terms_Mixed_Vol;

  function Facet_Normal ( n : integer32; facet,pts : VecVec ) return Vector is

    res,pt1,pt2 : Vector(1..n);
    im : Matrix(1..n,1..n);
    cnt : integer32 := 0;

  begin
    for i in facet'range loop
      if facet(i)'length > 1 then
        pt1 := pts(facet(i)(facet(i)'first)).all;
        for j in facet(i)'first+1..facet(i)'last loop
          pt2 := pts(facet(i)(j)).all;
          cnt := cnt + 1;
          for k in pt1'range loop
            im(cnt,k) := pt2(k) - pt1(k);
          end loop;
        end loop;
      end if;
    end loop;
    for j in 1..n loop
      im(n,j) := 0;
    end loop;
    Upper_Triangulate(im);
    Scale(im);
    res := (1..n => 0);
    Solve0(im,res);
    Normalize(res);
   -- put("The facet normal : "); put(res); new_line;
    return res;
  end Facet_Normal;
     
  function Minkowski_Sum ( n : natural32; al : Array_of_Lists )
                         return natural32 is

  -- DESCRIPTION :
  --   Computes the mixed volume of the polytope generated
  --   by the points in al, where n > 1.

    res,m,ptslen : natural32;
    vl,vl_last,al1 : List;
    typ,ind : Vector(1..integer32(n));
    perm,mix : Link_to_Vector;
    wal : Array_of_Lists(al'range) := Interchange2(al);

    procedure Update ( v : in Vector; i : in integer32;
		       added : in out boolean ) is

    -- DESCRIPTION :
    --   This procedure computes the support of the first list
    --   in the direction v; if this support is not zero,
    --   the projection will be computed.
    --   The parameter added becomes true if v has been added to vl.

      pal : Array_of_Lists(al'first..(al'last-1));
      pv : integer32;
      degen : boolean;

    begin
      if not Is_In(vl,v) then
        pv := Maximal_Support(al1,v);
        if pv > 0 then
          Projection(wal,v,i,pal,degen);
	  if not degen then
            declare
              mv : constant natural32 := Mixed_Volume(n-1,pal);
            begin
              if mv > 0 then
                res := res + natural32(pv)*mv;
                Append(vl,vl_last,v);
                added := true;
              end if;
            end;
            Deep_Clear(pal);
          end if;
        end if;
      end if;
    end Update;

    procedure Enumerate_Facets
                ( lpts : in Array_of_Lists; len : in natural32 ) is

      pts : VecVec(1..integer32(len));
      cnt : integer32;

      procedure Compute_Facet ( facet : in VecVec; cont : out boolean ) is

        v,w : Vector(1..integer32(n));
        i : integer32;
        added : boolean;

      begin
        v := Facet_Normal(integer32(n),facet,pts);
        i := Pivot(v);
        if i <= v'last then
          added := false;
          Update(v,i,added);
          if not added
	   then Min(v); w := v;
	  else w := -v; added := false;
          end if;
	  Update(w,i,added);
        end if;
        cont := true;
      end Compute_Facet;

      procedure Compute_Facets is new Enumerate_Faces_of_Sum(Compute_Facet);

    begin
      pts(ind(1)..ind(2)-1) := Shallow_Create(Tail_Of(al1));
      cnt := lpts'first + mix(mix'first);
      for i in mix'first+1..mix'last loop
        if i < mix'last then
          pts(ind(i)..ind(i+1)-1) := Shallow_Create(lpts(cnt));
          cnt := cnt + mix(i);
        else
          pts(ind(i)..integer32(len)) := Shallow_Create(lpts(cnt));
        end if;
      end loop;
      Compute_Facets(ind(mix'range),typ(mix'range),integer32(n)-1,pts);
    end Enumerate_Facets;
   
  begin
    m := Length_Of(wal(wal'first));
    if m = 2 then
      return Two_Terms_Mixed_Vol(n,wal);
    elsif m > 2 then
      Mixture(al,perm,mix);
     -- put("Type of mixture : "); put(mix); new_line;
     -- put(" with permutation : "); put(perm); new_line;
      wal := Permute(perm.all,al);
      declare
        shiftvec : Link_to_Vector;
        nulvec : constant Vector(1..integer32(n)) := (1..integer32(n) => 0);
        shifted : boolean;
        cnt : integer32;
      begin
       -- SHIFT OF THE FIRST LIST ( then all pv >= 0) :
        if not Is_In(wal(wal'first),nulvec) then
          shiftvec := Graded_Max(wal(wal'first));
          al1 := Shift(wal(wal'first),shiftvec); shifted := true;
        else
          al1 := wal(wal'first);                 shifted := false;
        end if;
       -- ENUMERATE FACES OF SUM OF THE RIGHT TYPE :
        Move_to_Front(al1,nulvec);
        wal(wal'first) := al1;
        res := 0;
        typ(1) := mix(mix'first)-1;
        typ(2) := mix(mix'first+1); ind(1) := 1;
        ind(2) := ind(1) + integer32(Length_Of(al1)) - 1;  -- skip null vector
        cnt := wal'first + mix(mix'first);
        for i in mix'first+2..mix'last loop
          typ(i) := mix(i);
          ind(i) := ind(i-1) + integer32(Length_Of(wal(cnt)));
          cnt := cnt + mix(i-1);
        end loop;
        ptslen := natural32(ind(mix'last)) + Length_Of(wal(cnt)) - 1;
        Enumerate_Facets(wal,ptslen);
       -- CLEANING UP :
        Deep_Clear(vl); Clear(perm); Clear(mix);
        if shifted 
         then Clear(shiftvec); Deep_Clear(al1);
        end if;
        return res;
      end;
    else -- m < 2
      return 0;
    end if;
  end Minkowski_Sum;

  function Mixed_Volume ( n : natural32; al : Array_of_Lists )
			return natural32 is
  begin
    if (n = 0) or Is_Null(al(al'first)) then
      return 0;
    elsif n = 1 then
      declare
        min,max : integer32 := 0;
        d : constant Link_to_Vector := Head_Of(al(al'first));
      begin
        Min_Max(al(al'first),d'first,min,max);
        return natural32(max - min);
      end;
    elsif All_Equal(al) then
      return Volume(n,al(al'first));
    else
      return Minkowski_Sum(n,al);
    end if;
  end Mixed_Volume;

-- MIXED VOLUME COMPUTATIONS WITH TREE OF USEFUL DIRECTIONS :

  function Two_Terms_Mixed_Volume
             ( n : natural32; al : Array_of_Lists; tv : Tree_of_Vectors )
             return natural32 is

  -- DESCRIPTION :
  --   returns the mixed volume of the polytopes generated by the
  --   points in al, where the first polytope is generated by only
  --   two points.

    first,second : Link_to_Vector;
    l : List renames al(al'first);
    res : natural32;

  begin
    first := Head_Of(l);
    second := Head_Of(Tail_Of(l));
    declare
      d : Vector(first'range);
      piv : integer32 := 0;
    begin
      for i in d'range loop
	d(i) := first(i) - second(i);
	if (piv = 0) and then (d(i) /= 0)
	 then piv := i;
        end if;
      end loop;
      if piv = 0 then
        return 0;
      else
        if d(piv) < 0
         then Min(d);
        end if;
        declare
          t : Transfo := Rotate(d,piv);
          tr_al : Array_of_Lists(al'first..(al'last-1));
          degen : boolean := false;
        begin
          Apply(t,d);
	  for i in tr_al'range loop
            tr_al(i) := Transform_and_Reduce(t,piv,al(i+1));
            Remove_Duplicates(tr_al(i));
            degen := (Length_Of(tr_al(i)) <= 1);
            exit when degen;
          end loop;
          Clear(t);
          if not degen
           then res := natural32(d(piv))*Mixed_Volume(n-1,tr_al,tv);
           else res := 0;
          end if;
          Deep_Clear(tr_al);
        end;
      end if;
    end;
    return res;
  end Two_Terms_Mixed_Volume;
     
  function Minkowski_Sum ( n : natural32; al : Array_of_Lists;
                           tv : Tree_of_Vectors ) return natural32 is

  -- DESCRIPTION :
  --   Computes the mixed volume of the polytope generated
  --   by the points in al, where n > 1.
  --   The tree of degrees is not empty.

    res,m : natural32;
    al1 : List;
    wal : Array_of_Lists(al'range) := Interchange2(al);
    tmp : Tree_of_Vectors;
    perm,mix : Link_to_Vector;

  begin
    m := Length_Of(wal(wal'first));
    if m = 2 then
      return Two_Terms_Mixed_Volume(n,wal,tv);
    elsif m > 2 then
      Mixture(al,perm,mix);
      wal := Permute(perm.all,al);
      declare
        shiftvec : Link_to_Vector;
        nulvec : constant Vector(1..integer32(n)) := (1..integer32(n) => 0);
        shifted : boolean;
      begin
       -- SHIFT OF THE FIRST LIST ( then all pv >= 0) :
        if not Is_In(wal(wal'first),nulvec) then
          shiftvec := Graded_Max(wal(wal'first));
      	  al1 := Shift(wal(wal'first),shiftvec); shifted := true;
        else
          al1 := wal(wal'first);                 shifted := false;
        end if;
        Move_to_Front(al1,nulvec);
        wal(wal'first) := al1;
       -- COMPUTING THE MIXED VOLUME :
        tmp := tv;  res := 0;
        while not Is_Null(tmp) loop
          declare
	    nd : constant node := Head_Of(tmp);
	    v : constant Link_to_Vector := nd.d;
	    i : constant integer32 := Pivot(v);
            pv : constant integer32 := Maximal_Support(al1,v.all);
            pal : Array_of_Lists(al'first..(al'last-1));
            degen : boolean;
          begin
            Projection(wal,v.all,i,pal,degen);
            if (nd.ltv = null) or else Is_Null(nd.ltv.all)
             then res := res + natural32(pv)*Mixed_Volume(n-1,pal);
             else res := res + natural32(pv)*Mixed_Volume(n-1,pal,nd.ltv.all);
            end if;
            Deep_Clear(pal);
          end;
          tmp := Tail_Of(tmp);
        end loop;
       -- CLEANING UP :
        Clear(perm); Clear(mix);
        if shifted 
         then Clear(shiftvec); Deep_Clear(al1);
        end if;
        return res;
      end;
    else -- m < 2
      return 0;
    end if;
  end Minkowski_Sum;

  function Mixed_Volume ( n : natural32; al : Array_of_Lists;
			  tv : Tree_of_Vectors ) return natural32 is
  begin
    if (n = 0) or Is_Null(al(al'first)) then
      return 0;
    elsif n = 1 then
      declare
        min,max : integer32 := 0;
        d : constant Link_to_Vector := Head_Of(al(al'first));
      begin
        Min_Max(al(al'first),d'first,min,max);
        return natural32(max - min);
      end;
    elsif All_Equal(al) then
      return Volume(n,al(al'first),tv);
    elsif not Is_Null(tv) then
      return Minkowski_Sum(n,al,tv);
    else
      return Minkowski_Sum(n,al);
    end if;
  end Mixed_Volume;

-- MIXED VOLUME COMPUTATIONS WITH CREATION OF TREE OF USEFUL DIRECTIONS :

  procedure Two_Terms_Mixed_Vol 
             ( n : in natural32; al : in Array_of_Lists;
	       tv : in out Tree_of_Vectors; mv : out natural32 ) is

  -- DESCRIPTION :
  --   returns the mixed volume of the polytopes generated by the
  --   points in al, where the first polytope is generated by only
  --   two points.

    first,second : Link_to_Vector;
    l : List renames al(al'first);

  begin
    first := Head_Of(l);
    second := Head_Of(Tail_Of(l));
    declare
      d : Vector(first'range);
      piv : integer32 := 0;
    begin
      for i in d'range loop
	d(i) := first(i) - second(i);
	if (piv = 0) and then (d(i) /= 0)
	 then piv := i;
        end if;
      end loop;
      if piv = 0 then
        mv := 0;
      else 
        if d(piv) < 0
         then Min(d);
        end if;
        declare
          t : Transfo := Rotate(d,piv);
          tr_al : Array_of_Lists(al'first..(al'last-1));
          mvl : natural32;
          degen : boolean := false;
        begin
          Apply(t,d);
          for i in tr_al'range loop
            tr_al(i) := Transform_and_Reduce(t,piv,al(i+1));
            Remove_Duplicates(tr_al(i));
            degen := (Length_Of(tr_al(i)) <= 1);
            exit when degen;
          end loop;
          Clear(t);
          if not degen then
            Mixed_Volume(n-1,tr_al,tv,mvl);
            mv := natural32(d(piv))*mvl;
          else
            mv := 0;
          end if;
          Deep_Clear(tr_al);
        end;
      end if;
    end;
  end Two_Terms_Mixed_Vol;
     
  procedure Minkowski_Sum
              ( n : in natural32; al : in Array_of_Lists;
                tv : in out Tree_of_Vectors; mv : out natural32 ) is

  -- DESCRIPTION :
  --   Computes the mixed volume of the polytope generated
  --   by the points in al, where n > 1.

    res,m,ptslen : natural32;
    al1 : List;
    typ,ind : Vector(1..integer32(n));
    perm,mix : Link_to_Vector;
    wal : Array_of_Lists(al'range) := Interchange2(al);

    procedure Update ( v : in Vector; i : in integer32;
		       added : in out boolean ) is

    -- DESCRIPTION :
    --   This procedure computes the support of the first list
    --   in the direction v; if this support is not zero,
    --   the projection will be computed.
    --   The parameter added becomes true if v has been added to vl.

      pal : Array_of_Lists(al'first..(al'last-1));
      pv : integer32;
      degen : boolean;

    begin
      if not Is_In(tv,v) then
        pv := Maximal_Support(al1,v);
        if pv > 0 then
          Projection(wal,v,i,pal,degen);
          if not degen then
            declare
              tmp : Tree_of_Vectors;
              mv : natural32;
            begin
              Mixed_Volume(n-1,pal,tmp,mv);
              if mv = 0 then
                Clear(tmp);
              else
                res := res + natural32(pv)*mv;
                declare
                  nd : node;
                begin
                  nd.d := new Vector'(v);
                  if not Is_Null(tmp)
                   then nd.ltv := new Tree_of_Vectors'(tmp);
                  end if;
                  Construct(nd,tv);
                end;
                added := true;
              end if;
            end;
            Deep_Clear(pal);
          end if;
        end if;
      end if;
    end Update;

    procedure Enumerate_Facets
                ( lpts : in Array_of_Lists; len : in natural32 ) is

      pts : VecVec(1..integer32(len));
      cnt : integer32;

      procedure Compute_Facet ( facet : in VecVec; cont : out boolean ) is

        v,w : Vector(1..integer32(n));
        i : integer32;
        added : boolean;

      begin
        v := Facet_Normal(integer32(n),facet,pts);
        i := Pivot(v);
        if i <= v'last then
          added := false;
          Update(v,i,added);
          if not added
           then Min(v); w := v;
           else w := -v; added := false;
          end if;
          Update(w,i,added);
        end if;
        cont := true;
      end Compute_Facet;

      procedure Compute_Facets is
        new Enumerate_Faces_of_Sum(Compute_Facet);

    begin
      pts(ind(1)..ind(2)-1) := Shallow_Create(Tail_Of(al1));
      cnt := lpts'first + mix(mix'first);
      for i in mix'first+1..mix'last loop
        if i < mix'last
         then pts(ind(i)..ind(i+1)-1) := Shallow_Create(lpts(cnt));
              cnt := cnt + mix(i);
         else pts(ind(i)..integer32(len)) := Shallow_Create(lpts(cnt));
        end if;
      end loop;
      Compute_Facets(ind(mix'range),typ(mix'range),integer32(n-1),pts);
    end Enumerate_Facets;

  begin
    m := Length_Of(wal(wal'first));
    if m = 2 then
      Two_Terms_Mixed_Vol(n,wal,tv,mv);
    elsif m > 2 then
      Mixture(al,perm,mix);
      wal := Permute(perm.all,al);
      declare
        shiftvec : Link_to_Vector;
        nulvec : constant Vector(1..integer32(n)) := (1..integer32(n) => 0);
        shifted : boolean;
        cnt : integer32;
      begin
       -- SHIFT OF THE FIRST LIST ( then all pv >= 0) :
        if not Is_In(wal(wal'first),nulvec) then
          shiftvec := Graded_Max(wal(wal'first));
          al1 := Shift(wal(wal'first),shiftvec); shifted := true;
        else
          al1 := wal(wal'first);                 shifted := false;
        end if;
       -- ENUMERATE FACES OF SUM OF THE RIGHT TYPE :
        Move_to_Front(al1,nulvec);
        wal(wal'first) := al1;
        res := 0;
        typ(1) := mix(mix'first)-1;
        typ(2) := mix(mix'first+1); ind(1) := 1;
        ind(2) := ind(1) + integer32(Length_Of(al1)) - 1;  -- skip null vector
        cnt := wal'first + mix(mix'first);
        for i in mix'first+2..mix'last loop
          typ(i) := mix(i);
          ind(i) := ind(i-1) + integer32(Length_Of(wal(cnt)));
          cnt := cnt + mix(i-1);
        end loop;
        ptslen := natural32(ind(mix'last)) + Length_Of(wal(cnt)) - 1;
        Enumerate_Facets(wal,ptslen);
       -- CLEANING UP :
        Clear(perm); Clear(mix);
        if shifted 
         then Clear(shiftvec); Deep_Clear(al1);
        end if;
        mv := res;
      end;
    else -- m < 2
      mv := 0;
    end if;
  end Minkowski_Sum;

  procedure Mixed_Volume
              ( n : in natural32; al : in Array_of_Lists;
                tv : in out Tree_of_Vectors; mv : out natural32 ) is
  begin
    if (n = 0) or Is_Null(al(al'first)) then
      mv := 0;
    elsif n = 1 then
      declare
        min,max : integer32 := 0;
        d : constant Link_to_Vector := Head_Of(al(al'first));
      begin
        Min_Max(al(al'first),d'first,min,max);
        mv := natural32(max - min);
      end;
    elsif All_Equal(al) then
      Volume(n,al(al'first),tv,mv);
    else
      Minkowski_Sum(n,al,tv,mv);
    end if;
  end Mixed_Volume;

end Volumes;
