with integer_io;                         use integer_io;
with Standard_Floating_Numbers;          use Standard_Floating_Numbers;
with Standard_Floating_Numbers_io;       use Standard_Floating_Numbers_io;
with Standard_Mathematical_Functions;    use Standard_Mathematical_Functions;
with Standard_Complex_Numbers;           use Standard_Complex_Numbers;
with Standard_Floating_Vectors_io;       use Standard_Floating_Vectors_io;
with Standard_Complex_Vectors;
with Standard_Complex_Vectors_io;        use Standard_Complex_Vectors_io;
with Standard_Complex_Norms_Equals;      use Standard_Complex_Norms_Equals;
with Standard_Integer_VecVecs;
with Standard_Floating_VecVecs;
with Standard_Complex_Matrices;          use Standard_Complex_Matrices;
with Standard_Complex_Polynomials;
with Standard_Complex_Laur_Polys;
with Standard_Complex_Laur_Functions;    use Standard_Complex_Laur_Functions;
with Standard_Complex_Poly_Systems;      use Standard_Complex_Poly_Systems;
with Standard_Laur_Poly_Convertors;      use Standard_Laur_Poly_Convertors;
with Standard_Complex_Substitutors;      use Standard_Complex_Substitutors;
with Power_Lists;                        use Power_Lists;
with Lists_of_Floating_Vectors;          use Lists_of_Floating_Vectors;
with Floating_Lifting_Utilities;         use Floating_Lifting_Utilities;
with Floating_Integer_Convertors;        use Floating_Integer_Convertors;
with Transforming_Integer_Vector_Lists;  use Transforming_Integer_Vector_Lists;
with Transforming_Laurent_Systems;       use Transforming_Laurent_Systems;
with Continuation_Parameters;
with Increment_and_Fix_Continuation;     use Increment_and_Fix_Continuation;
with Fewnomial_System_Solvers;           use Fewnomial_System_Solvers;
with BKK_Bound_Computations;             use BKK_Bound_Computations;

package body Floating_Polyhedral_Continuation is

-- UTILITIES FOR POLYHEDRAL COEFFICIENT HOMOTOPY :

  procedure put ( file : in file_type;
                  v : in Standard_Floating_Vectors.Vector;
                  fore,aft,exp : in natural ) is
  begin
    for i in v'range loop
      put(file,v(i),fore,aft,exp);
    end loop;
  end put;

  procedure put ( file : in file_type;
                  v : in Standard_Complex_Vectors.Vector;
                  fore,aft,exp : in natural ) is
  begin
    for i in v'range loop
      put(file,REAL_PART(v(i)),fore,aft,exp);
      put(file,IMAG_PART(v(i)),fore,aft,exp);
    end loop;
  end put;

  function Power ( x,m : double_float ) return double_float is

  -- DESCRIPTION :
  --   Computes x**m for high powers of m to avoid RANGE_ERROR.

   -- intm : integer := integer(m);
   -- fltm : double_float;

  begin
   -- if m < 10.0
   --  then
    return (x**m);  -- GNAT is better at this
   --  else if double_float(intm) > m
   --        then intm := intm-1;
   --       end if;
   --       fltm := m - double_float(intm);
   --       return ((x**intm)*(x**fltm));
   -- end if;
  end Power;

  function Minimum ( v : Standard_Floating_Vectors.Vector )
                   return double_float is

  -- DESCRIPTION :
  --   Returns the minimal element (/= 0)  in the vector v.

    tol : constant double_float := 10.0**(-7);
    min : double_float := 0.0;
    tmp : double_float;

  begin
    for i in v'range loop
      tmp := abs(v(i));
      if tmp > tol
       then if (min = 0.0) or else (tmp < min)
             then min := tmp;
            end if;
      end if;
    end loop;
    return min;
  end Minimum;

  function Minimum ( v : Standard_Floating_VecVecs.VecVec )
                   return double_float is

  -- DESCRIPTION :
  --   Returns the minimal nonzero element in the vectors of v.

    tol : constant double_float := 10.0**(-7);
    min : double_float := 0.0;
    tmp : double_float;

  begin
    for i in v'range loop
      tmp := Minimum(v(i).all);
      if abs(tmp) > tol
       then if (min = 0.0) or else (tmp < min)
             then min := tmp;
            end if;
      end if;
    end loop;
    return min;
  end Minimum;

  function Scale ( v : Standard_Floating_Vectors.Vector; s : double_float )
                 return Standard_Floating_Vectors.Vector is

  -- DESCRIPTION :
  --   Returns the scaled vector, by dividing every element by s.

    res : Standard_Floating_Vectors.Vector(v'range);

  begin
    for i in res'range loop
      res(i) := v(i)/s;
    end loop;
    return res;
  end Scale;

  function Scale ( v : Standard_Floating_VecVecs.VecVec )
                 return Standard_Floating_VecVecs.VecVec is

  -- DESCRIPTION :
  --   Returns an array of scaled vectors.

    res : Standard_Floating_VecVecs.VecVec(v'range);
    min : constant double_float := Minimum(v);
    tol : constant double_float := 10.0**(-8);

  begin
    if (abs(min) > tol) and (min /= 1.0)
     then
       for i in v'range loop
         res(i) := new Standard_Floating_Vectors.Vector'(Scale(v(i).all,min));
       end loop;
     else 
       for i in v'range loop
         res(i) := new Standard_Floating_Vectors.Vector'(v(i).all);
       end loop;
    end if;
    return res;
  end Scale;

  procedure Shift ( v : in out Standard_Floating_Vectors.Vector ) is

  -- DESCRIPTION :
  --   Shifts the elements in v, such that the minimal element equals zero.

    min : double_float := v(v'first);

  begin
    for i in v'first+1..v'last loop
      if v(i) < min
       then min := v(i);
      end if;
    end loop;
    if min /= 0.0
     then for i in v'range loop
            v(i) := v(i) - min;
          end loop;
    end if;
  end Shift;

  function Create ( e : Standard_Integer_VecVecs.VecVec;
                    l : List; normal :  Standard_Floating_Vectors.Vector ) 
                  return Standard_Floating_Vectors.Vector is

  -- DESCRIPTION :
  --   Returns a vector with all inner products of the normal with
  --   the exponents in the list, such that minimal value equals zero.

    res : Standard_Floating_Vectors.Vector(e'range);
    use Standard_Floating_Vectors;
    found : boolean;
    lif : double_float;

  begin
    for i in e'range loop
      declare
        fei : constant Standard_Floating_Vectors.Vector := Convert(e(i).all);
      begin
        Search_Lifting(l,fei,found,lif);
        if not found
         then res(i) := 1000.0;
         else res(i) := fei*normal(fei'range) + lif*normal(normal'last);
        end if;
      end;
    end loop;
    Shift(res);
    return res;
  end Create;

  function Create ( e : Exponent_Vectors_Array;
                    l : Array_of_Lists; mix : Standard_Integer_Vectors.Vector;
                    normal : Standard_Floating_Vectors.Vector )
                  return Standard_Floating_VecVecs.VecVec is

    res : Standard_Floating_VecVecs.VecVec(normal'first..normal'last-1);
    cnt : natural := res'first; 

  begin
    for i in mix'range loop
      declare
        rescnt : constant Standard_Floating_Vectors.Vector
               := Create(e(cnt).all,l(i),normal);
      begin
        res(cnt) := new Standard_Floating_Vectors.Vector(rescnt'range);
        for j in rescnt'range loop
          res(cnt)(j) := rescnt(j);
        end loop;
      end;
      Shift(res(cnt).all);
      for k in 1..(mix(i)-1) loop
        res(cnt+k) := new Standard_Floating_Vectors.Vector(res(cnt)'range);
        for j in res(cnt)'range loop
          res(cnt+k)(j) := res(cnt)(j);
        end loop;
      end loop;
      cnt := cnt + mix(i);
    end loop;
    return res;
  end Create;

  procedure Eval ( c : in Standard_Complex_Vectors.Vector;
                   t : in double_float; m : in Standard_Floating_Vectors.Vector;
                   ctm : in out Standard_Complex_Vectors.Vector ) is

  -- DESCRIPTION : ctm = c*t**m.

  begin
    for i in ctm'range loop
      -- ctm(i) := c(i)*Create((t**m(i)));
      ctm(i) := c(i)*Create(Power(t,m(i)));
    end loop;
  end Eval;

  procedure Eval ( c : in Standard_Complex_VecVecs.VecVec;
                   t : in double_float; m : in Standard_Floating_VecVecs.VecVec;
                   ctm : in out Standard_Complex_VecVecs.VecVec ) is

  -- DESCRIPTION : ctm = c*t**m.

  begin
    for i in ctm'range loop
      Eval(c(i).all,t,m(i).all,ctm(i).all);
    end loop;
  end Eval;

-- USEFUL PRIMITIVES :

  function Is_In ( l : List;
                   d : Standard_Complex_Laur_Polys.Degrees )
                 return boolean is

  -- DESCRIPTION :
  --   Returns true if the degrees belong to the list l.

    tmp : List := l;
    pt : Standard_Floating_Vectors.Link_to_Vector;
    fld : Standard_Floating_Vectors.Vector(d'range);
    equ : boolean;

  begin
    for i in fld'range loop
      fld(i) := double_float(d(i));
    end loop;
    while not Is_Null(tmp) loop
      pt := Head_Of(tmp);
      equ := true;
      for i in fld'range loop
        if pt(i) /= fld(i)
         then equ := false;
        end if;
        exit when not equ;
      end loop;
      if equ
       then return true;
       else tmp := Tail_Of(tmp);
      end if;
    end loop;
    return false;
  end Is_In;

  function Select_Terms ( p : Standard_Complex_Laur_Polys.Poly; l : List )
                        return Standard_Complex_Laur_Polys.Poly is

    use Standard_Complex_Laur_Polys;
    res : Poly := Null_Poly;

    procedure Visit_Term ( t : in Term; cont : out boolean ) is
    begin
      if Is_In(l,t.dg)
       then Add(res,t);
      end if;
      cont := true;
    end Visit_Term;
    procedure Visit_Terms is new Visiting_Iterator(Visit_Term);

  begin
    Visit_Terms(p);
    return res;
  end Select_Terms;

  function Select_Subsystem
               ( p : Laur_Sys; mix : Standard_Integer_Vectors.Vector;
                 mic : Mixed_Cell ) return Laur_Sys is

  -- DESCRIPTION :
  --   Given a Laurent polynomial system and a mixed cell,
  --   the corresponding subsystem will be returned.

  -- ON ENTRY :
  --   p          a Laurent polynomial system;
  --   mix        type of mixture: occurencies of the supports;
  --   mic        a mixed cell.

  -- REQUIRED :
  --   The polynomials in p must be ordered according to the type of mixture.

    res : Laur_Sys(p'range);
    cnt : natural := 0;

  begin
    for k in mix'range loop
      for l in 1..mix(k) loop
        cnt := cnt + 1;
        res(cnt) := Select_Terms(p(cnt),mic.pts(k));
      end loop;
    end loop;
    return res;
  end Select_Subsystem;

  procedure Extract_Regular ( sols : in out Solution_List ) is

    function To_Be_Removed ( flag : in natural ) return boolean is
    begin
      return ( flag /= 1 );
    end To_Be_Removed;
    procedure Extract_Regular_Solutions is new 
                Standard_Complex_Solutions.Delete(To_Be_Removed);

  begin
    Extract_Regular_Solutions(sols);
  end Extract_Regular;

  procedure Refine ( file : in file_type; p : in Laur_Sys;
                     sols : in out Solution_List ) is

  -- DESCRIPTION :
  --   Given a polyhedral homotopy p and a list of solution for t=1,
  --   this list of solutions will be refined.

    pp : Poly_Sys(p'range) := Laurent_to_Polynomial_System(p);
    n : constant natural := p'length;
    eps : constant double_float := 10.0**(-12);
    tolsing : constant double_float := 10.0**(-8);
    max : constant natural := 3;
    numb : natural := 0;

  begin
    pp := Laurent_to_Polynomial_System(p);
    Substitute(n+1,Create(1.0),pp);
   -- Reporting_Root_Refiner(file,pp,sols,eps,eps,tolsing,numb,max,false);
    Clear(pp); Extract_Regular(sols);
  end Refine;

-- FIRST LAYER OF CONTINUATION ROUTINES :

  procedure Mixed_Continuation
                ( mix : in Standard_Integer_Vectors.Vector;
                  lifted : in Array_of_Lists; h : in Eval_Coeff_Laur_Sys;
                  c : in Standard_Complex_VecVecs.VecVec;
                  e : in Exponent_Vectors_Array;
                  j : in Eval_Coeff_Jaco_Mat; m : in Mult_Factors;
                  normal : in Standard_Floating_Vectors.Vector;
                  sols : in out Solution_List ) is

    pow : Standard_Floating_VecVecs.VecVec(c'range)
        := Create(e,lifted,mix,normal);
    scapow : Standard_Floating_VecVecs.VecVec(c'range) := Scale(pow);
    ctm : Standard_Complex_VecVecs.VecVec(c'range);

    use Standard_Floating_Vectors;
    use Standard_Complex_Laur_Polys;

    function Eval ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
                  return Standard_Complex_Vectors.Vector is
    begin
      Eval(c,REAL_PART(t),scapow,ctm);
      return Eval(h,ctm,x);
    end Eval;

    function dHt ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
                 return Standard_Complex_Vectors.Vector is

      res : Standard_Complex_Vectors.Vector(h'range);
      xtl : constant integer := x'last+1;

    begin
      Eval(c,REAL_PART(t),scapow,ctm);
      for i in res'range loop
        res(i) := Eval(j(i,xtl),m(i,xtl).all,ctm(i).all,x);
      end loop;
      return res;
    end dHt;

    function dHx ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
                 return matrix is

      mt : Matrix(x'range,x'range);

    begin
      Eval(c,REAL_PART(t),scapow,ctm);
      for k in mt'range(1) loop
        for l in mt'range(2) loop
          mt(k,l) := Eval(j(k,l),m(k,l).all,ctm(k).all,x);
        end loop;
      end loop;
      return mt;
    end dHx;

    procedure Laur_Cont is new Silent_Continue(Max_Norm,Eval,dHt,dHx);

  begin
    for i in c'range loop
      ctm(i) 
        := new Standard_Complex_Vectors.Vector'(c(i).all'range => Create(0.0));
    end loop;
    Laur_Cont(sols,false);
    Standard_Floating_VecVecs.Clear(pow);
    Standard_Floating_VecVecs.Clear(scapow);
    Standard_Complex_VecVecs.Clear(ctm);
    Extract_Regular(sols);
  end Mixed_Continuation;

  procedure Mixed_Continuation
                ( file : in file_type;
                  mix : in Standard_Integer_Vectors.Vector;
                  lifted : in Array_of_Lists; h : in Eval_Coeff_Laur_Sys;
                  c : in Standard_Complex_VecVecs.VecVec;
                  e : in Exponent_Vectors_Array;
                  j : in Eval_Coeff_Jaco_Mat; m : in Mult_Factors;
                  normal : in Standard_Floating_Vectors.Vector;
                  sols : in out Solution_List ) is

    pow : Standard_Floating_VecVecs.VecVec(c'range)
        := Create(e,lifted,mix,normal);
    scapow : Standard_Floating_VecVecs.VecVec(c'range) := Scale(pow);
    ctm : Standard_Complex_VecVecs.VecVec(c'range);

    use Standard_Floating_Vectors;
    use Standard_Complex_Laur_Polys;

    function Eval ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
                  return Standard_Complex_Vectors.Vector is
    begin
      Eval(c,REAL_PART(t),scapow,ctm);
      return Eval(h,ctm,x);
    end Eval;

    function dHt ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
                 return Standard_Complex_Vectors.Vector is

      res : Standard_Complex_Vectors.Vector(h'range);
      xtl : constant integer := x'last+1;

    begin
      Eval(c,REAL_PART(t),scapow,ctm);
      for i in res'range loop
        res(i) := Eval(j(i,xtl),m(i,xtl).all,ctm(i).all,x);
      end loop;
      return res;
    end dHt;

    function dHx ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
                 return Matrix is

      mt : Matrix(x'range,x'range);

    begin
      Eval(c,REAL_PART(t),scapow,ctm);
      for k in m'range(1) loop
        for l in m'range(1) loop
          mt(k,l) := Eval(j(k,l),m(k,l).all,ctm(k).all,x);
        end loop;
      end loop;
      return mt;
    end dHx;

    procedure Laur_Cont is new Reporting_Continue(Max_Norm,Eval,dHt,dHx);

  begin
   -- put_line(file,"The coefficient vectors :" );
   -- for i in c'range loop
   --   put(file,c(i).all,3,3,3); new_line(file);
   -- end loop;
    for i in c'range loop
      ctm(i) 
        := new Standard_Complex_Vectors.Vector'(c(i).all'range => Create(0.0));
    end loop;
   -- put(file,"The normal : "); put(file,normal,3,3,3); new_line(file);
   -- put_line(file,"The exponent vector : ");
   -- for i in pow'range loop
   --   put(file,pow(i).all,3,3,3); new_line(file);
   -- end loop;
   -- put_line(file,"The scaled exponent vector : ");
   -- for i in pow'range loop
   --   put(file,scapow(i).all,3,3,3); new_line(file);
   -- end loop;
    Laur_Cont(file,sols,false);
    Standard_Floating_VecVecs.Clear(pow);
    Standard_Floating_VecVecs.Clear(scapow);
    Standard_Complex_VecVecs.Clear(ctm);
    Extract_Regular(sols);
  end Mixed_Continuation;

-- UTILITIES FOR SECOND LAYER :

  function Remove_Lifting ( l : List ) return List is

  -- DESCRIPTION :
  --   Removes the lifting value from the vectors.

    tmp,res,res_last : List;

  begin
    tmp := l;
    while not Is_Null(tmp) loop
      declare
        d1 : constant Standard_Floating_Vectors.Vector := Head_Of(tmp).all;
        d2 : constant Standard_Floating_Vectors.Vector
           := d1(d1'first..d1'last-1);
      begin
        Append(res,res_last,d2);
      end;
      tmp := Tail_Of(tmp);
    end loop;
    return res;
  end Remove_Lifting;

  function Sub_Lifting ( q : Laur_Sys; mix : Standard_Integer_Vectors.Vector;
                         mic : Mixed_Cell ) return Array_of_Lists is

  -- DESCRIPTION :
  --   Returns the lifting used to subdivide the cell.

    res : Array_of_Lists(mix'range);
    sup : Array_of_Lists(q'range);
    n : constant natural := q'last;
    cnt : natural := sup'first;

  begin
    for i in mic.pts'range loop
      sup(cnt) := Remove_Lifting(mic.pts(i));
      for j in 1..(mix(i)-1) loop
        Copy(sup(cnt),sup(cnt+j));
      end loop;
      cnt := cnt + mix(i);
    end loop;
    res := Induced_Lifting(n,mix,sup,mic.sub.all);
    Deep_Clear(sup);
    return res;
  end Sub_Lifting;

  function Sub_Polyhedral_Homotopy
               ( l : List; e : Standard_Integer_VecVecs.VecVec;
                 c : Standard_Complex_Vectors.Vector )
               return Standard_Complex_Vectors.Vector is

  -- DESCRIPTION :
  --   For every vector in e that does not belong to l, the corresponding
  --   index in c will be set to zero, otherwise it is copied to the result.

    res : Standard_Complex_Vectors.Vector(c'range);
    found : boolean;
    lif : double_float;

  begin
    for i in e'range loop
      declare
        fei : constant Standard_Floating_Vectors.Vector := Convert(e(i).all);
      begin
        Search_Lifting(l,fei,found,lif);
        if not found
         then res(i) := Create(0.0);
         else res(i) := c(i);
        end if;
      end;
    end loop;
    return res;
  end Sub_Polyhedral_Homotopy;

  function Sub_Polyhedral_Homotopy
               ( mix : Standard_Integer_Vectors.Vector; mic : Mixed_Cell;
                 e : Exponent_Vectors_Array;
                 c : Standard_Complex_VecVecs.VecVec )
               return Standard_Complex_VecVecs.VecVec is

  -- DESCRIPTION :
  --   Given a subsystem q of p and the coefficient vector of p, the
  --   vector on return will have only nonzero entries for coefficients
  --   that belong to q.

    res : Standard_Complex_VecVecs.VecVec(c'range);

  begin
    for i in mix'range loop
      declare
        cri : constant Standard_Complex_Vectors.Vector
            := Sub_Polyhedral_Homotopy(mic.pts(i),e(i).all,c(i).all);
      begin
        res(i) := new Standard_Complex_Vectors.Vector'(cri);
        for j in 1..(mix(i)-1) loop
          declare
            crj : constant Standard_Complex_Vectors.Vector
                := Sub_Polyhedral_Homotopy(mic.pts(i),e(i+j).all,c(i+j).all);
          begin
            res(i+j) := new Standard_Complex_Vectors.Vector'(crj);
          end;
        end loop;
      end;
    end loop;
    return res;
  end Sub_Polyhedral_Homotopy;

  procedure Refined_Mixed_Solve
                ( q : in Laur_Sys; mix : in Standard_Integer_Vectors.Vector;
                  mic : in Mixed_Cell; h : in Eval_Coeff_Laur_Sys;
                  c : in Standard_Complex_VecVecs.VecVec;
                  e : in Exponent_Vectors_Array;
                  j : in Eval_Coeff_Jaco_Mat; m : in Mult_Factors;
                  qsols : in out Solution_List ) is

  -- DESCRIPTION :
  --   Polyhedral coeffient-homotopy for subsystem q.

  -- REQUIRED : mic.sub /= null.

    lif : Array_of_Lists(mix'range) := Sub_Lifting(q,mix,mic);
    cq : Standard_Complex_VecVecs.VecVec(c'range)
       := Sub_Polyhedral_Homotopy(mix,mic,e,c);

  begin
    Mixed_Solve(q,lif,h,cq,e,j,m,mix,mic.sub.all,qsols);
    Standard_Complex_VecVecs.Clear(cq); Deep_Clear(lif);
  end Refined_Mixed_Solve;

  procedure Refined_Mixed_Solve
                ( file : in file_type; q : in Laur_Sys;
                  mix : in Standard_Integer_Vectors.Vector;
                  mic : in Mixed_Cell; h : in Eval_Coeff_Laur_Sys;
                  c : in Standard_Complex_VecVecs.VecVec;
                  e : in Exponent_Vectors_Array;
                  j : in Eval_Coeff_Jaco_Mat; m : in Mult_Factors;
                  qsols : in out Solution_List ) is

  -- DESCRIPTION :
  --   Polyhedral coeffient-homotopy for subsystem q.

  -- REQUIRED : mic.sub /= null.

    lif : Array_of_Lists(mix'range) := Sub_Lifting(q,mix,mic);
    cq : Standard_Complex_VecVecs.VecVec(c'range)
       := Sub_Polyhedral_Homotopy(mix,mic,e,c);

  begin
    Mixed_Solve(file,q,lif,h,cq,e,j,m,mix,mic.sub.all,qsols);
    Standard_Complex_VecVecs.Clear(cq); Deep_Clear(lif);
  end Refined_Mixed_Solve;

-- SECOND LAYER :

  procedure Mixed_Solve
                ( p : in Laur_Sys; lifted : in Array_of_Lists;
                  h : in Eval_Coeff_Laur_Sys;
                  c : in Standard_Complex_VecVecs.VecVec;
                  e : in Exponent_Vectors_Array;
                  j : in Eval_Coeff_Jaco_Mat; m : in Mult_Factors;
                  mix : in Standard_Integer_Vectors.Vector;
                  mic : in Mixed_Cell;
                  sols,sols_last : in out Solution_List ) is

    q : Laur_Sys(p'range) := Select_Subsystem(p,mix,mic);
    sq : Laur_Sys(q'range) := Shift(q);
    pq : Poly_Sys(q'range);
    qsols : Solution_List;
    len : natural := 0;
    fail : boolean;

  begin
    Solve(sq,qsols,fail);
    if fail
     then if mic.sub = null
           then pq := Laurent_to_Polynomial_System(sq);
                qsols := Solve_by_Static_Lifting(pq);
                Clear(pq);
           else Refined_Mixed_Solve(q,mix,mic,h,c,e,j,m,qsols);
          end if;
          Set_Continuation_Parameter(qsols,Create(0.0));
    end if;
    len := Length_Of(qsols);
    if len > 0
     then Mixed_Continuation(mix,lifted,h,c,e,j,m,mic.nor.all,qsols);
          Concat(sols,sols_last,qsols);
    end if;
    Clear(sq); Clear(q); Clear(qsols);
  end Mixed_Solve;

  procedure Mixed_Solve
                ( file : in file_type;
                  p : in Laur_Sys; lifted : in Array_of_Lists;
                  h : in Eval_Coeff_Laur_Sys;
                  c : in Standard_Complex_VecVecs.VecVec;
                  e : in Exponent_Vectors_Array;
                  j : in Eval_Coeff_Jaco_Mat; m : in Mult_Factors;
                  mix : in Standard_Integer_Vectors.Vector;
                  mic : in Mixed_Cell;
                  sols,sols_last : in out Solution_List ) is

    q : Laur_Sys(p'range) := Select_Subsystem(p,mix,mic);
    sq : Laur_Sys(q'range) := Shift(q);
    pq : Poly_Sys(q'range);
    qsols : Solution_List;
    len : natural := 0;
    fail : boolean;

  begin
    Solve(sq,qsols,fail);
    if not fail
     then put_line(file,"It is a fewnomial system.");
     else put_line(file,"No fewnomial system.");
          if mic.sub = null
           then put_line(file,"Calling the black box solver.");
                pq := Laurent_to_Polynomial_System(sq);
                qsols := Solve_by_Static_Lifting(file,pq);
                Clear(pq);
           else put_line(file,"Using the refinement of the cell.");
                Refined_Mixed_Solve(file,q,mix,mic,h,c,e,j,m,qsols);
          end if;
          Set_Continuation_Parameter(qsols,Create(0.0));
    end if;
    len := Length_Of(qsols);
    put(file,len,1); put_line(file," solutions found.");
    if len > 0
     then 
       Mixed_Continuation(file,mix,lifted,h,c,e,j,m,mic.nor.all,qsols);
       Concat(sols,sols_last,qsols);
    end if;
    Clear(sq); Clear(q); Clear(qsols);
  end Mixed_Solve;

-- THIRD LAYER :

  procedure Mixed_Solve
                ( p : in Laur_Sys; lifted : in Array_of_Lists;
                  h : in Eval_Coeff_Laur_Sys;
                  c : in Standard_Complex_VecVecs.VecVec;
                  e : in Exponent_Vectors_Array;
                  j : in Eval_Coeff_Jaco_Mat; m : in Mult_Factors;
                  mix : in Standard_Integer_Vectors.Vector;
                  mixsub : in Mixed_Subdivision;
                  sols : in out Solution_List ) is

    tmp : Mixed_Subdivision := mixsub;
    mic : Mixed_Cell;
    sols_last : Solution_List;

  begin
    while not Is_Null(tmp) loop
      mic := Head_Of(tmp);
      Mixed_Solve(p,lifted,h,c,e,j,m,mix,mic,sols,sols_last);
      tmp := Tail_Of(tmp);
    end loop;
  end Mixed_Solve;

  procedure Mixed_Solve
                ( file : in file_type;
                  p : in Laur_Sys; lifted : in Array_of_Lists;
                  h : in Eval_Coeff_Laur_Sys;
                  c : in Standard_Complex_VecVecs.VecVec;
                  e : in Exponent_Vectors_Array;
                  j : in Eval_Coeff_Jaco_Mat; m : in Mult_Factors;
                  mix : in Standard_Integer_Vectors.Vector;
                  mixsub : in Mixed_Subdivision;
                  sols : in out Solution_List ) is

    tmp : Mixed_Subdivision := mixsub;
    mic : Mixed_Cell;
    sols_last : Solution_List;
    cnt : natural := 0;

  begin
    while not Is_Null(tmp) loop
      mic := Head_Of(tmp);
      cnt := cnt + 1;
      new_line(file);
      put(file,"*** PROCESSING SUBSYSTEM "); put(file,cnt,1);
      put_line(file," ***");
      new_line(file);
      Mixed_Solve(file,p,lifted,h,c,e,j,m,mix,mic,sols,sols_last);
      tmp := Tail_Of(tmp);
    end loop;
  end Mixed_Solve;

end Floating_Polyhedral_Continuation;