with text_io;                            use text_io;
with Communications_with_User;           use Communications_with_User;
with Standard_Natural_Numbers;           use Standard_Natural_Numbers;
with Standard_Integer_Numbers;           use Standard_Integer_Numbers;
with Standard_Integer_Numbers_io;        use Standard_Integer_Numbers_io;
with Standard_Complex_Numbers;           use Standard_Complex_Numbers;
with Standard_Random_Numbers;            use Standard_Random_Numbers;
with Standard_Natural_Vectors;
with Standard_Integer_Vectors;           use Standard_Integer_Vectors;  
with Standard_Integer_VecVecs;           use Standard_Integer_VecVecs;
with Standard_Integer_VecVecs_io;        use Standard_Integer_VecVecs_io;
with Lists_of_Integer_Vectors;           use Lists_of_Integer_Vectors;
with Lists_of_Integer_Vectors_io;        use Lists_of_Integer_Vectors_io;
with Standard_Complex_Polynomials;       use Standard_Complex_Polynomials;
with Standard_Complex_Poly_Systems;      use Standard_Complex_Poly_Systems;
with Standard_Complex_Poly_Systems_io;   use Standard_Complex_Poly_Systems_io;

procedure ts_stair is

-- DESCRIPTION :
--   The user of this program is prompted for a dimension and multiplicity.
--   A monomial system with prescribed multiplicity is generated by
--   choosing m standard monomials under a staircase at random.

  function Is_In ( v : VecVec; x : Vector ) return boolean is

  -- DESCRIPTION :
  --   Returns true if x belongs to v.

  begin
    for i in v'range loop
      if v(i) /= null then
        if Equal(v(i).all,x)
         then return true;
        end if;
      end if;
    end loop;
    return false;
  end Is_In;

  function Divisors_In ( v : VecVec; x : Vector ) return boolean is

  -- DESCRIPTION :
  --   Returns true if all possible decrements of x belong to v.

    y : Vector(x'range) := x;

  begin
    for i in x'range loop
      if x(i) > 0 then
        y(i) := x(i) - 1;
        if not Is_In(v,y)
         then return false;
        end if;
        y(i) := x(i);
      end if;
    end loop;
    return true;
  end Divisors_In;

  function Is_Corner ( v : VecVec; x : Vector ) return boolean is

  -- DESCRIPTION :
  --   The vector x is a corner to v if x does not belong to v
  --   and if all possible divisors of x belong to v.

  begin
    if Is_In(v,x) then
      return false;
    elsif Divisors_In(v,x) then
      return true;
    else
      return false;
    end if;
  end Is_Corner;

  function Random_Staircase ( n,m : integer32 ) return VecVec is

  -- DESCRIPTION :
  --   Returns m monomials under a random staircase in n variables.
  --   The vector on return is of range 0..m-1.

    res : VecVec(0..m-1);
    zero : constant Vector(1..n) := (1..n => 0);
    wrk : Vector(1..n);
    k,j : integer32;

  begin
    res(0) := new Vector'(zero);
    for i in 1..m-1 loop               -- m-1 monomials to add
      k := Random(1,n);                -- pick a direction to increment
      j := Random(0,i-1);              -- pick a monomail to increment
      loop
        wrk := res(j).all;
        wrk(k) := wrk(k) + 1;          -- increment and add if corner
        if Is_Corner(res,wrk)
         then res(i) := new Vector'(wrk); j := i;
         else j := (j+1) mod i;        -- take next available monomial
        end if;
        exit when (j = i);
      end loop;
    end loop;
    return res;
  end Random_Staircase;

  function Corners ( s : VecVec; n : integer32 ) return List is

  -- DESCRIPTION :
  --   Returns all corners to the monomials under the staircase s.

    res,res_last : List;
    wrk : Vector(1..n);

  begin
    for i in s'range loop              -- for all monomials under staircase
      wrk := s(i).all; 
      for k in 1..n loop               -- test all possible increments
        wrk(k) := wrk(k) + 1;          -- and add increment if corner
        if Is_Corner(s,wrk)
         then Append_Diff(res,res_last,wrk);
        end if;
        wrk(k) := wrk(k) - 1;
      end loop;
    end loop;
    return res;
  end Corners;

  function Create ( e : Standard_Integer_Vectors.Vector ) return Poly is

  -- DESCRIPTION :
  --   Returns the monomial with exponent vectors in e.

    res : Poly;
    t : Term;

  begin
    t.cf := Create(1.0);
    t.dg := new Standard_Natural_Vectors.Vector(e'range);
    for i in e'range loop
      t.dg(i) := natural32(e(i));
    end loop;
    res := Create(t);
    Clear(t);
    return res;
  end Create;

  function Create ( L : List ) return Poly_Sys is

  -- DESCRIPTION :
  --   Returns a monomial system in n variables whose equations are
  --   the monomials with exponent vectors in the given list L.

    res : Poly_Sys(1..integer32(Length_Of(L)));
    tmp : List := L;

  begin
    for i in res'range loop
      res(i) := Create(Head_Of(tmp).all);
      tmp := Tail_Of(tmp);
    end loop;
    return res;
  end Create;

  procedure Write_System ( n,m : in integer32; L : List ) is

  -- DESCRIPTION :
  --   Writes the polynomial system defined by the exponents in L.

     neq : constant integer32 := integer32(Length_Of(L));
     p : constant Poly_Sys(1..neq) := Create(L);
     file : file_type;

  begin
    put(natural32(neq),natural32(n),p);
    new_line;
    put_line("Reading the name of the output file ...");
    Read_Name_and_Create_File(file);
    put(file,natural32(neq),natural32(n),p);
    new_line(file);
    put(file,"TITLE : origin is only root, isolated and with multiplicity ");
    put(file,m,1); new_line(file);
  end Write_System;

  procedure Monomial_System ( n,m : in integer32 ) is

  -- DESCRIPTION :
  --   Generates a monomial system in n variables where the origin
  --   occurs as the only solution, isolated and of multiplicity m.

    s : constant VecVec(0..m-1) := Random_Staircase(n,m);
    c : constant List := Corners(s,n);

  begin
    put_line("Monomials under the staircase : "); put(s);
    put_line("The corners of the staircase : "); put(c);
    put_line("The monomial system : ");
    Write_System(n,m,c);
  end Monomial_System;

  procedure Main is

    n,m : integer32 := 0;

  begin
    new_line;
    put_line("Monomial systems with prescribed multiple isolated root ...");
    new_line;
    put("Give #variables : "); get(n);
    put("Give multiplicity : "); get(m); skip_line;
    Monomial_System(n,m);
  end Main;

begin
  Main;
end ts_stair;
