with Standard_Floating_Numbers;          use Standard_Floating_Numbers;
with Standard_Floating_Matrices;
with Standard_Integer_Norms;             use Standard_Integer_Norms;
with Standard_Integer_Matrices;          use Standard_Integer_Matrices;
with Standard_Integer_Linear_Solvers;    use Standard_Integer_Linear_Solvers;
with Standard_Integer_Linear_Equalities; use Standard_Integer_Linear_Equalities;
with Integer_Linear_Inequalities;        use Integer_Linear_Inequalities;
with Floating_Linear_Inequality_Solvers; use Floating_Linear_Inequality_Solvers;

package body Integer_Pruning_Methods is

-- GENERAL AUXILIARIES :

  function Convert ( fa : Face_Array ) return Array_of_Lists is

  -- DESCRIPTION :
  --   Converts the array of faces into an array of lists by
  --   converting the first element of each list of faces.

    res : Array_of_Lists(fa'range);

  begin
    for k in fa'range loop
      res(k) := Shallow_Create(fa(k).all);
    end loop;
    return res;
  end Convert;

-- AUXILIARIES FOR THE PRUNING ALGORITHMS :

  procedure Initialize ( n : in natural; mat : out Matrix;
                         ipvt : out Vector ) is

  -- DESCRIPTION :
  --   Initializes the matrix of the equality constraints on the normals
  --   and sets the pivoting vector ipvt to 1..n+1.

  begin
    for i in 1..n+1 loop
      ipvt(i) := i;
      for j in 1..n+1 loop
        mat(i,j) := 0;
      end loop;
    end loop;
  end Initialize;

  function Number_of_Inequalities
             ( mix : Vector; lifted : Array_of_Lists ) return natural is

  -- DESCRIPTION :
  --   Returns the maximal number of inequalities on the inner normal.

    res : natural := 0;

  begin
    for k in lifted'range loop
      res := res + Length_Of(lifted(k)) - mix(k) - 1;
    end loop;
    return res;
  end Number_of_Inequalities;

  procedure Ordered_Inequalities ( k : in natural; mat : in out matrix ) is

  -- DESCRIPTION :
  --   Defines k inequalities mat(k,k) - mat(k+1,k) >= 0.

  begin
    for i in mat'first(1)..k loop
      for j in mat'range(2) loop
        mat(i,j) := 0;
      end loop;
      mat(i,i) := 1; mat(i,i+1) := -1;
    end loop;
  end Ordered_Inequalities;

  procedure Check_and_Update
                ( mic : in Face_Array; lifted : in Array_of_Lists;
                  m : in matrix; ipvt : in Vector;
                  mixsub,mixsub_last : in out Mixed_Subdivision ) is

  -- DESCRIPTION :
  --   Computes the normal to the points in mic, by solving the 
  --   linear system defined by m and ipvt.  Updates the mixed subdivision
  --   if the computed normal is an inner normal.

    v : Vector(m'range(2)) := Solve(m,ipvt);
    pts : Array_of_Lists(mic'range);

  begin
    if v(v'last) /= 0
     then Normalize(v);
          if v(v'last) < 0
           then Min(v);
          end if;
          pts := Convert(mic);
          Update(pts,v,mixsub,mixsub_last);
    end if;
  end Check_and_Update;

  procedure Create_Equalities
                ( n : in natural; f : in Face; mat,ineq : in matrix;
                  ipvt : in Vector; row,rowineq : in natural;
                  newmat,newineq : in out matrix; newipvt : in out Vector;
                  newrow,newrowineq : in out natural; fail : out boolean ) is

  -- DESCRIPTION :
  --   Creates new equalities and uses them to eliminate unknowns in the
  --   matrices of equalities and inequalities.  Failure is reported when
  --   a zero row is encountered.  On entry, all new* variables must be
  --   initialized with the corresponding *-ones.

    shi : vector(1..n+1) := f(f'first).all;
    fl : boolean := false;
    pivot : natural;

  begin
    for i in f'first+1..f'last loop
      newrow := newrow + 1;
      for j in f(i)'range loop
        newmat(newrow,j) := f(i)(j) - shi(j);
      end loop;
      Triangulate(newrow,newmat,newipvt,pivot);
      fl := (pivot = 0);
      exit when fl;
      Switch(newrow,pivot,ineq'first,rowineq,newineq);
      --Scale(newrow,newmat);
    end loop;
    fail := fl;
  end Create_Equalities;

  function Check_Feasibility
              ( k,m,n : natural; ineq : matrix ) return boolean is
 
  -- DESCRIPTION :
  --   Returns true if -v(n+1) > 0 can be derived, false otherwise.

  -- ON ENTRY :
  --   k      current unknown that has been eliminated,
  --           for all i <= k : ineq(l,i) = 0, for l in ineq'first..m;
  --   m      number of inequalities;
  --   n      dimension;
  --   ineq   matrix of inequalities.

    tableau : matrix(1..n-k+1,1..m+1);
    feasi : boolean;

  begin
    if m = 0
     then feasi := false;
     else for i in k+1..n+1 loop
            for j in 1..m loop
              tableau(i-k,j) := ineq(j,i);
            end loop;
            tableau(i-k,m+1) := 0;
          end loop;
          tableau(n-k+1,m+1) := -1;
          Integer_Complementary_Slackness(tableau,feasi);
    end if;
    return feasi;
  end Check_Feasibility;

  function New_Check_Feasibility
              ( k,m,n : natural; ineq : matrix; tol : double_float;
                solu : Standard_Floating_Vectors.Vector ) return boolean is

  -- DESCRIPTION :
  --   Returns true if -v(n+1) > 0 can be derived, false otherwise.

  -- ON ENTRY :
  --   k      current unknown that has been eliminated,
  --           for all i <= k : ineq(l,i) = 0, for l in ineq'first..m;
  --   m      number of inequalities;
  --   n      dimension;
  --   ineq   matrix of inequalities.

    tableau : Standard_Floating_Matrices.matrix(1..n-k+1,1..m);
    tabsolu : Standard_Floating_Vectors.Vector(1..n-k);
    cols : Vector(tabsolu'range);
    kfail : integer;
    max,tmpabs : double_float;
         -- used for scaling of the last component of the inner normal

  begin
    for i in k+1..n loop
      for j in 1..m loop
        tableau(i-k,j) := double_float(ineq(j,i));
      end loop;
      tabsolu(i-k) := solu(i);
    end loop;
    max := 0.0;
    for j in 1..m loop
      tableau(n-k+1,j) := -double_float(ineq(j,n+1));
      tmpabs := abs(tableau(n-k+1,j));
      if tmpabs > max
       then max := tmpabs;
      end if;
    end loop;
    if max > 1.0
     then for j in 1..m loop
            tableau(n-k+1,j) := tableau(n-k+1,j)/max;
          end loop;
    end if;
    Scale(tableau,tol);
    Solve(tableau,tol,tabsolu,kfail,cols);
    return (kfail <= m);
  end New_Check_Feasibility;

  procedure Update_Inequalities
               ( k,rowmat1,rowmat2,n : in natural;
                 mat : in matrix; ipvt : in vector; 
                 rowineq : in out natural; ineq : in out matrix;
                 lifted : in Array_of_Lists; mic : in out Face_Array ) is

  -- DESCRIPTION :
  --   The inequalities will be updated w.r.t. the equality
  --   constraints on the inner normal.

    tmp : List;
    pt,shi : Link_to_Vector;

  begin
    for i in ineq'first..rowineq loop   -- triangulate old inequalities
      for j in rowmat1..rowmat2 loop
        Triangulate(j,mat,i,ineq);
      end loop;
      --Scale(i,ineq);
    end loop;
    tmp := lifted(k);                         -- create and triangulate 
    shi := mic(k)(mic(k)'first);                    -- new inequalities
    while not Is_Null(tmp) loop
      pt := Head_Of(tmp);
      if not Is_In(mic(k),pt.all)
       then rowineq := rowineq + 1;
            for j in pt'range loop
              ineq(rowineq,j) := pt(j) - shi(j);
            end loop;
            Switch(ipvt,rowineq,ineq);
            for i in 1..rowmat2 loop
              Triangulate(i,mat,rowineq,ineq);
              --Scale(rowineq,ineq);
            end loop;
      end if;
      tmp := Tail_Of(tmp);
    end loop;
  end Update_Inequalities;

  procedure Eliminate ( k : in natural; m : in Matrix; tol : in double_float;
                        x : in out Standard_Floating_Vectors.Vector ) is

  -- DESCRIPTION :
  --   Eliminates the kth component of x by using the matrix.

   fac : double_float;

  begin
    if abs(x(k)) > tol
     then fac := x(k)/double_float(m(k,k));
          for j in k..x'last loop
            x(j) := x(j) - fac*double_float(m(k,j));
          end loop;
    end if;
  end Eliminate;

  procedure Switch ( l,pivot : in integer;
                     x : in out Standard_Floating_Vectors.Vector ) is

  -- DESCRIPTION :
  --   Applies the pivotation information to the vector x.

    tmp : double_float;

  begin
    if l /= pivot
     then tmp := x(l); x(l) := x(pivot); x(pivot) := tmp;
    end if;
  end Switch;

  procedure Eliminate ( rowmat1,rowmat2 : in natural; mat : in Matrix;
                        ipvt : in Vector; tol : in double_float;
                        x : in out Standard_Floating_Vectors.Vector ) is

  -- DESCRIPTION :
  --   Returns an eliminated solution vector.

  begin
    for k in rowmat1..rowmat2 loop
      Switch(k,ipvt(k),x);
      Eliminate(k,mat,tol,x);
    end loop;
  end Eliminate;

-- GENERAL CONSTRUCTOR in several versions :

  -- procedure Compute_Mixed_Cells ( );  -- general specification.

  -- DESCRIPTION :
  --   Backtrack recursive procedure to enumerate face-face combinations.

  -- ON ENTRY :
  --   k         index for current point set;
  --   row       number of rows already in the matrix;
  --   mat       matrix which determines the inner normal;
  --   ipvt      contains the pivoting information;
  --   rowineq   number of inequality constraints already in ineq;
  --   ineq      matrix for the inequality constraints on the
  --             inner normal v, of type <.,v> >= 0;
  --   testsolu  test solution to check current set of inequalilities;
  --   mic       contains the current selected faces, up to k-1.

  -- ON RETURN :
  --   mic       updated selected faces;
  --   issupp    true if current tuple of faces is supported;
  --   continue  indicates whether to continue the creation or not.

-- CONSTRUCTION WITH PRUNING :

  procedure Gen1_Create_CS
               ( n : in natural; mix : in Vector; fa : in Array_of_Faces; 
                 lifted : in Array_of_Lists;
                 nbsucc,nbfail : in out Standard_Floating_Vectors.Vector;
                 mixsub : out Mixed_Subdivision ) is

    res,res_last : Mixed_Subdivision;
    accu : Face_Array(fa'range);
    ma : matrix(1..n+1,1..n+1);
    ipvt : vector(1..n+1);
    ineqrows : natural;

    procedure Compute_Mixed_Cells
                 ( k,row : in natural; mat : in matrix; ipvt : in vector;
                   rowineq : in natural; ineq : in matrix;
                   mic : in out Face_Array; continue : out boolean );

    -- DESCRIPTION :
    --   Backtrack recursive procedure to enumerate face-face combinations.

    procedure Process_Inequalities
                 ( k,rowmat1,rowmat2 : in natural;
                   mat : in matrix; ipvt : in vector;
                   rowineq : in out natural; ineq : in out matrix;
                   mic : in out Face_Array; cont : out boolean ) is

    -- DESCRIPTION :
    --   Updates inequalities and checks feasibility before proceeding.

    begin
      Update_Inequalities(k,rowmat1,rowmat2,n,mat,ipvt,rowineq,ineq,lifted,mic);
      if Check_Feasibility(rowmat2,rowineq,n,ineq)
       then nbfail(k) := nbfail(k) + 1.0;
            cont := true;
       else nbsucc(k) := nbsucc(k) + 1.0;
            Compute_Mixed_Cells(k+1,rowmat2,mat,ipvt,rowineq,ineq,mic,cont);
      end if;
    end Process_Inequalities;

    procedure Compute_Mixed_Cells
                 ( k,row : in natural; mat : in matrix; ipvt : in vector;
                   rowineq : in natural; ineq : in matrix;
                   mic : in out Face_Array; continue : out boolean ) is

      old : Mixed_Subdivision := res_last;
      cont : boolean := true;
      tmpfa : Faces;

    begin
      if k > mic'last
       then Check_and_Update(mic,lifted,mat,ipvt,res,res_last);
            if old /= res_last
             then Process(Head_Of(res_last),continue);
             else continue := true;
            end if;
       else tmpfa := fa(k);
            while not Is_Null(tmpfa) loop  -- enumerate faces of kth polytope
              mic(k) := Head_Of(tmpfa);
              declare                                    -- update matrices
                fl : boolean;
                newipvt : vector(ipvt'range) := ipvt;
                newmat : matrix(mat'range(1),mat'range(2)) := mat;
                newineq : matrix(ineq'range(1),ineq'range(2)) := ineq;
                newrow : natural := row;
                newrowineq : natural := rowineq;
              begin
                Create_Equalities(n,mic(k),mat,ineq,ipvt,row,rowineq,
                                  newmat,newineq,newipvt,newrow,newrowineq,fl);
                if fl
                 then nbfail(k) := nbfail(k) + 1.0;
                 else Process_Inequalities(k,row+1,newrow,newmat,newipvt,
                                           newrowineq,newineq,mic,cont);
                end if;
              end;
              exit when not cont;
              tmpfa := Tail_Of(tmpfa);
           end loop;
           continue := cont;
      end if;
    end Compute_Mixed_Cells;

  begin
    Initialize(n,ma,ipvt);
    ineqrows := Number_of_Inequalities(mix,lifted);
    declare
      ineq : matrix(1..ineqrows,1..n+1);
      cont : boolean;
    begin
      ineq(1,1) := 0;
      Compute_Mixed_Cells(accu'first,0,ma,ipvt,0,ineq,accu,cont);
    end;
    mixsub := res;
  end Gen1_Create_CS;

  procedure Create_CS
              ( n : in natural; mix : in Vector; fa : in Array_of_Faces; 
                lifted : in Array_of_Lists;
                nbsucc,nbfail : in out Standard_Floating_Vectors.Vector;
                mixsub : out Mixed_Subdivision ) is

    res,res_last : Mixed_Subdivision;
    accu : Face_Array(fa'range);
    ma : matrix(1..n+1,1..n+1);
    ipvt : vector(1..n+1);
    ineqrows : natural;

    procedure Compute_Mixed_Cells
                 ( k,row : in natural; mat : in Matrix; ipvt : in Vector;
                   rowineq : in natural; ineq : in Matrix;
                   mic : in out Face_Array );

    -- DESCRIPTION :
    --   Backtrack recursive procedure to enumerate face-face combinations.

    procedure Process_Inequalities
                 ( k,rowmat1,rowmat2 : in natural;
                   mat : in matrix; ipvt : in vector;
                   rowineq : in out natural; ineq : in out matrix;
                   mic : in out Face_Array ) is

    -- DESCRIPTION :
    --   Updates inequalities and checks feasibility before proceeding.

    begin
      Update_Inequalities(k,rowmat1,rowmat2,n,mat,ipvt,rowineq,ineq,lifted,mic);
      if Check_Feasibility(rowmat2,rowineq,n,ineq)
       then nbfail(k) := nbfail(k) + 1.0;
       else nbsucc(k) := nbsucc(k) + 1.0;
            Compute_Mixed_Cells(k+1,rowmat2,mat,ipvt,rowineq,ineq,mic);
      end if;
    end Process_Inequalities;

    procedure Compute_Mixed_Cells 
                 ( k,row : in natural; mat : in Matrix; ipvt : in Vector;
                   rowineq : in natural; ineq : in Matrix;
                   mic : in out Face_Array ) is

      tmpfa : Faces;

    begin
      if k > mic'last
       then Check_and_Update(mic,lifted,mat,ipvt,res,res_last);
       else tmpfa := fa(k);
            while not Is_Null(tmpfa) loop   -- enumerate faces of kth polytype
              mic(k) := Head_Of(tmpfa);
              declare                                     -- update matrices
                fl : boolean;
                newipvt : vector(ipvt'range) := ipvt;
                newmat : matrix(mat'range(1),mat'range(2)) := mat;
                newineq : matrix(ineq'range(1),ineq'range(2)) := ineq;
                newrow : natural := row;
                newrowineq : natural := rowineq;
              begin
                Create_Equalities(n,mic(k),mat,ineq,ipvt,row,rowineq,newmat,
                                  newineq,newipvt,newrow,newrowineq,fl);
                if fl 
                 then nbfail(k) := nbfail(k) + 1.0;
                 else Process_Inequalities
                       (k,row+1,newrow,newmat,newipvt,newrowineq,newineq,mic);
                end if;
              end;
              tmpfa := Tail_Of(tmpfa);
            end loop;
      end if;
    end Compute_Mixed_Cells;

  begin
    Initialize(n,ma,ipvt);
    ineqrows := Number_of_Inequalities(mix,lifted);
    declare
      ineq : matrix(1..ineqrows,1..n+1);
    begin
      ineq(1,1) := 0;
      Compute_Mixed_Cells(accu'first,0,ma,ipvt,0,ineq,accu);
    end;
    mixsub := res;
  end Create_CS;

  procedure New_Create_CS
              ( n : in natural; mix : in Vector; fa : in Array_of_Faces; 
                lifted : in Array_of_Lists;
                nbsucc,nbfail : in out Standard_Floating_Vectors.Vector;
		mixsub : out Mixed_Subdivision ) is

    res,res_last : Mixed_Subdivision;
    accu : Face_Array(fa'range);
    ma : matrix(1..n+1,1..n+1);
    ipvt : vector(1..n+1);
    tol : constant double_float := double_float(n)*10.0**(-11);
    ineqrows : natural;

    procedure Compute_Mixed_Cells
                 ( k,row : in natural; mat : in matrix; ipvt : in vector;
                   rowineq : in natural; ineq : in matrix;
                   testsolu : in Standard_Floating_Vectors.Vector;
                   mic : in out Face_Array );

    -- DESCRIPTION :
    --   Backtrack recursive procedure to enumerate face-face combinations.

    procedure Process_Inequalities
                 ( k,rowmat1,rowmat2 : in natural;
                   mat : in matrix; ipvt : in vector;
                   rowineq : in out natural; ineq : in out matrix;
                   testsolu : in Standard_Floating_Vectors.Vector;
                   mic : in out Face_Array ) is

    -- DESCRIPTION :
    --   Updates the inequalities and checks feasibility before proceeding.

      tmp : List;
      pt,shi : Link_to_Vector;
      newsolu : Standard_Floating_Vectors.Vector(testsolu'range) := testsolu;

    begin
      Update_Inequalities(k,rowmat1,rowmat2,n,mat,ipvt,rowineq,ineq,
                          lifted,mic);
      Eliminate(rowmat1,rowmat2,mat,ipvt,tol,newsolu);
      if New_Check_Feasibility(rowmat2,rowineq,n,ineq,tol,newsolu)
       then nbfail(k) := nbfail(k) + 1.0;
       else nbsucc(k) := nbsucc(k) + 1.0;
            Compute_Mixed_Cells(k+1,rowmat2,mat,ipvt,rowineq,ineq,testsolu,mic);
      end if;
    end Process_Inequalities;

    procedure Compute_Mixed_Cells 
                 ( k,row : in natural; mat : in matrix; ipvt : in vector;
                   rowineq : in natural; ineq : in matrix;
                   testsolu : in Standard_Floating_Vectors.Vector;
                   mic : in out Face_Array ) is

    -- DESCRIPTION :
    --   Backtrack recursive procedure to enumerate face-face combinations.

      tmpfa : Faces;

    begin
      if k > mic'last
       then Check_and_Update(mic,lifted,mat,ipvt,res,res_last);
       else tmpfa := fa(k); 
            while not Is_Null(tmpfa) loop -- enumerate faces of kth polytope
              mic(k) := Head_Of(tmpfa);
              declare                                     -- update matrices
                fl : boolean;
                newipvt : vector(ipvt'range) := ipvt;
                newmat : matrix(mat'range(1),mat'range(2)) := mat;
                newineq : matrix(ineq'range(1),ineq'range(2)) := ineq;
                newrow : natural := row;
                newrowineq : natural := rowineq;
              begin
                Create_Equalities(n,mic(k),mat,ineq,ipvt,row,rowineq,newmat,
                                  newineq,newipvt,newrow,newrowineq,fl);
                if fl 
                 then nbfail(k) := nbfail(k) + 1.0;
                 else Process_Inequalities(k,row+1,newrow,newmat,newipvt,
                                           newrowineq,newineq,testsolu,mic);
                end if;
              end;
              tmpfa := Tail_Of(tmpfa);
            end loop;
      end if;
    end Compute_Mixed_Cells;

  begin
    Initialize(n,ma,ipvt);
    ineqrows := Number_of_Inequalities(mix,lifted);
    declare
      ineq : matrix(1..ineqrows,1..n+1);
      solu : Standard_Floating_Vectors.Vector(1..n+1);
    begin
      ineq(1,1) := 0;
      solu := (solu'range => 0.0);
      Compute_Mixed_Cells(accu'first,0,ma,ipvt,0,ineq,solu,accu);
    end;
    mixsub := res;
  end New_Create_CS;

-- AUXILIARIES FOR THE CONSTRAINT PRUNING :

  function Is_Supported ( f : Face; normal : Vector; supp : integer )
                        return boolean is

    ip : integer;

  begin
    for i in f'range loop
      ip := f(i).all*normal;
      if ip /= supp
       then return false;
      end if;
    end loop;
    return true;
  end Is_Supported;

  function Is_Supported ( f : Face; k : natural; normals,suppvals : List )
                        return boolean is

  -- DESCRIPTION :
  --   Returns true if there exists a normal for which the inner product
  --   with each point in the face equals the corresponding supporting value
  --   of the kth component.

    tmpnor : List := normals;
    tmpsup : List := suppvals;
    support : integer;

  begin
    while not Is_Null(tmpnor) loop
      support := Head_Of(tmpsup)(k);
      if Is_Supported(f,Head_Of(tmpnor).all,support)
       then return true;
       else tmpnor := Tail_Of(tmpnor);
            tmpsup := Tail_Of(tmpsup);
      end if;
    end loop;
    return false;
  end Is_Supported;
 
  procedure Update ( mic : in Mixed_Cell; normals,suppvals : in out List ) is

  -- DESCRIPTION :
  --   Updates the list of normals and supporting values with the information
  --   of a mixed cell.

    normal : Link_to_Vector := new Vector'(mic.nor.all);
    suppval : Link_to_Vector := new Vector(mic.pts'range);

  begin
    Construct(normal,normals);
    for i in suppval'range loop
      suppval(i) := normal*Head_Of(mic.pts(i));
    end loop;
    Construct(suppval,suppvals);
  end Update;

  procedure Create_CCS
                 ( n : in natural; mix : in Vector; fa : in Array_of_Faces; 
                   lifted : in Array_of_Lists;
                   nbsucc,nbfail : in out Standard_Floating_Vectors.Vector;
                   normals,suppvals : in out List;
                   mixsub : out Mixed_Subdivision ) is

    res,res_last : Mixed_Subdivision;
    accu : Face_Array(fa'range);
    ma : matrix(1..n+1,1..n+1);
    ipvt : vector(1..n+1);
    ineqrows : natural;

    procedure Compute_Mixed_Cells
                 ( k,row : in natural; mat : in matrix; ipvt : in vector;
                   rowineq : in natural; ineq : in matrix;
                   mic : in out Face_Array; issupp : in out boolean );

    -- DESCRIPTION :
    --   Backtrack recursive procedure to enumerate face-face combinations.

    procedure Process_Inequalities
                 ( k,rowmat1,rowmat2 : in natural;
                   mat : in matrix; ipvt : in vector;
                   rowineq : in out natural; ineq : in out matrix;
                   mic : in out Face_Array; issupp : in out boolean ) is

    -- DESCRIPTION :
    --   Updates inequalities and checks feasibility before proceeding.

    begin
      Update_Inequalities(k,rowmat1,rowmat2,n,mat,ipvt,rowineq,ineq,lifted,mic);
      if issupp
       then issupp := Is_Supported(mic(k),k,normals,suppvals);
      end if;
      if not issupp and then Check_Feasibility(rowmat2,rowineq,n,ineq)
       then nbfail(k) := nbfail(k) + 1.0;
       else nbsucc(k) := nbsucc(k) + 1.0;
            Compute_Mixed_Cells(k+1,rowmat2,mat,ipvt,rowineq,ineq,mic,issupp);
      end if;
    end Process_Inequalities;

    procedure Compute_Mixed_Cells 
                 ( k,row : in natural; mat : in matrix; ipvt : in vector;
                   rowineq : in natural; ineq : in matrix;
                   mic : in out Face_Array; issupp : in out boolean ) is

    -- DESCRIPTION :
    --   Backtrack recursive procedure to enumerate face-face combinations.

      old : Mixed_Subdivision;
      tmpfa : Faces;

    begin
      if k > mic'last
       then old := res_last;
            Check_and_Update(mic,lifted,mat,ipvt,res,res_last);
            if old /= res_last
             then Update(Head_Of(res_last),normals,suppvals);
            end if;
       else tmpfa := fa(k);
            while not Is_Null(tmpfa) loop  -- enumerate faces of kth polytope
              mic(k) := Head_Of(tmpfa);
              declare                                      -- update matrices
                fl : boolean;
                newipvt : vector(ipvt'range) := ipvt;
                newmat : matrix(mat'range(1),mat'range(2)) := mat;
                newineq : matrix(ineq'range(1),ineq'range(2)) := ineq;
                newrow : natural := row;
                newrowineq : natural := rowineq;
              begin
                Create_Equalities(n,mic(k),mat,ineq,ipvt,row,rowineq,newmat,
                                  newineq,newipvt,newrow,newrowineq,fl);
                if fl 
                 then nbfail(k) := nbfail(k) + 1.0;
                 else Process_Inequalities(k,row+1,newrow,newmat,newipvt,
                                           newrowineq,newineq,mic,issupp);
                end if;
              end;
              tmpfa := Tail_Of(tmpfa);
            end loop;
      end if;
    end Compute_Mixed_Cells;

  begin
    Initialize(n,ma,ipvt);
    ineqrows := Number_of_Inequalities(mix,lifted);
    declare
      ineq : matrix(1..ineqrows,1..n+1);
      supported : boolean := true;
    begin
      ineq(1,1) := 0;
      Compute_Mixed_Cells(accu'first,0,ma,ipvt,0,ineq,accu,supported);
    end;
    mixsub := res;
  end Create_CCS;

end Integer_Pruning_Methods;