Return to interpolating_homotopies.adb CVS log

Up to [local] / OpenXM_contrib / PHC / Ada / Root_Counts / Product

Import the second public release of PHCpack.

OKed by Jan Verschelde.

with Standard_Floating_Numbers;          use Standard_Floating_Numbers;
with Standard_Complex_Numbers;           use Standard_Complex_Numbers;
with Standard_Natural_Vectors;
with Standard_Complex_Vectors;
with Standard_Complex_Matrices;
with Standard_Complex_Linear_Solvers;    use Standard_Complex_Linear_Solvers;
with Standard_Complex_Polynomials;       use Standard_Complex_Polynomials;
with Standard_Complex_Poly_Randomizers;  use Standard_Complex_Poly_Randomizers;
with Sets_of_Unknowns;                   use Sets_of_Unknowns;
with Degrees_in_Sets_of_Unknowns;        use Degrees_in_Sets_of_Unknowns;

package body Interpolating_Homotopies is

-- AUXILIARY OPERATIONS  :

  function Initial_Term ( p : Poly ) return Term is

  -- DESCRIPTION :
  --   Returns the initial term of p.

    res : Term;

    procedure Init_Term ( t : Term; continue : out boolean ) is
    begin
      res := t;
      continue := false;
    end Init_Term;
    procedure Init_Terms is new Visiting_Iterator(Init_Term);

  begin
    Init_Terms(p);
    return res;
  end Initial_Term;

  procedure Eliminate_Term ( p : in out Poly; dg : in Degrees ) is

  -- DESCRIPTION :
  --   Eliminates the monomial in p that has the given exponent vector.

    c : Complex_Number := Coeff(p,dg);

  begin
    if c /= Create(0.0)
     then declare
            t : Term;
          begin
            t.cf := c; t.dg := dg;
            Sub(p,t);
          end;
    end if;
  end Eliminate_Term;

  function Admitted ( t : Term; i : natural; z : partition;
                      d : Standard_Integer_Matrices.Matrix ) return boolean is

  -- DESCRIPTION :
  --   Returns true if the given term does not violate the m-homogeneous
  --   structure, i.e., if Degree(t,z(j)) <= d(i,j), for j in z'range.

  begin
    for j in z'range loop
      if Degree(t,z(j)) > d(i,j)
       then return false; 
      end if;
    end loop;
    return true; 
  end Admitted;

  function Dense_Representation ( p : Poly; i : natural; z : Partition;
                                  d : Matrix ) return Poly is

  -- DESCRIPTION :
  --   Returns the dense representation of the given polynomial p, known as
  --   p(i) of some polynomial system, of the given m-homogeneous structure.
  --   The returned polynomial has all its coefficients equal to one.

    res : Poly := Null_Poly;
    maxdegs,accu : Standard_Natural_Vectors.Vector(d'range(1));

    procedure Generate_Monomials
                  ( k : in natural; max : in Standard_Natural_Vectors.Vector;
                    acc : in out Standard_Natural_Vectors.Vector) is

    -- DESCRIPTION :
    --   Generates all monomials with exponents in a box, with upper bounds
    --   given by the vector max.  Every monomial that respects the 
    --   m-homogeneous structure will be added to the result res.
    --   The current unknown is indicated by the parameter k.

      t : Term;

    begin
      for j in 0..max(k) loop
        acc(k) := j;
        t.cf := Create(1.0);
        t.dg := new Standard_Natural_Vectors.Vector'(acc);
        if Admitted(t,i,z,d)
         then Add(res,t);
        end if;
        Clear(t);
        if k < max'last
         then Generate_Monomials(k+1,max,acc);
        end if;
      end loop;
    end Generate_Monomials;

  begin
    for j in maxdegs'range loop
      maxdegs(j) := Degree(p,j); accu(j) := 0;
    end loop;
    Generate_Monomials(1,maxdegs,accu);
    return res;
  end Dense_Representation;

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

    res : Complex_Number := Create(1.0);

  begin
    for i in x'range loop
      for k in 1..t.dg(i) loop
        res := res*x(i);
      end loop;
    end loop;
    return res;
  end Evaluate;

  procedure Interpolation_System
                ( p : in Poly; b : in natural; sols : in Solution_List;
                  mat : out Standard_Complex_Matrices.Matrix;
                  rhs : out Standard_Complex_Vectors.Vector ) is

  -- DESCRIPTION :
  --   Returns the matrix and right hand side vector of the linear system
  --   that expresses the interpolationg conditions.

    m : Standard_Complex_Matrices.Matrix(1..b,1..b);
    v : Standard_Complex_Vectors.Vector(1..b);
    cnt : natural := 0;

    procedure Scan_Term ( t : in Term; continue : out boolean ) is

      tmp : Solution_List := sols;
      s : Link_to_Solution;

    begin
      for row in m'range(1) loop        -- fill in column indicated by cnt
        s := Head_Of(tmp);
        if cnt = 0
         then v(row) := -t.cf*Evaluate(t,s.v);
         elsif cnt <= b
             then m(row,cnt) := Evaluate(t,s.v);
             else v(row) := v(row) - t.cf*Evaluate(t,s.v);
        end if;
        tmp := Tail_Of(tmp);
      end loop;
      cnt := cnt + 1;
      continue := true;
    end Scan_Term;
    procedure Scan_Poly is new Visiting_Iterator(Scan_Term);

  begin
    Scan_Poly(p);
    mat := m; rhs := v;
  end Interpolation_System;

  procedure Interpolation_System
                ( p : in Poly; i,b : in natural; sols : in Solution_List;
                  mat : out Standard_Complex_Matrices.Matrix;
                  rhs : out Standard_Complex_Vectors.Vector ) is

  -- DESCRIPTION :
  --   Returns the matrix and right hand side vector of the linear system
  --   that expresses the interpolationg conditions.  Monomials with degree
  --   one in x_i will not appear in the interpolation matrix.

    m : Standard_Complex_Matrices.Matrix(1..b,1..b);
    v : Standard_Complex_Vectors.Vector(1..b);
    cnt : natural := 0;

    procedure Scan_Term ( t : in Term; continue : out boolean ) is

      tmp : Solution_List := sols;
      s : Link_to_Solution;

    begin
      for row in m'range(1) loop        -- fill in column indicated by cnt
        s := Head_Of(tmp);
        if cnt = 0
         then v(row) := -t.cf*Evaluate(t,s.v);
         else if cnt <= b and t.dg(i) /= 1
               then m(row,cnt) := Evaluate(t,s.v);
               else v(row) := v(row) - t.cf*Evaluate(t,s.v);
              end if;
        end if;
        tmp := Tail_Of(tmp);
      end loop;
      if cnt = 0 or t.dg(i) /= 1
       then cnt := cnt + 1;
      end if;
      continue := true;
    end Scan_Term;
    procedure Scan_Poly is new Visiting_Iterator(Scan_Term);

  begin
    Scan_Poly(p);
    mat := m; rhs := v;
  end Interpolation_System;

  procedure Construct_Interpolant ( p : in out Poly;
                                    cv : in Standard_Complex_Vectors.Vector ) is

  -- DESCRIPTION :
  --   With the coefficients of its monomials, the interpolant will be
  --   constructed.

    cnt : natural := 0;

    procedure Scan_Term ( t : in out Term; continue : out boolean ) is
    begin
      if cnt > cv'last
       then continue := false;
       else if cnt >= cv'first
             then t.cf := cv(cnt);
            end if;
            cnt := cnt + 1;
            continue := true;
      end if;
    end Scan_Term;
    procedure Scan_Poly is new Changing_Iterator(Scan_Term);

  begin
    Scan_Poly(p);
  end Construct_Interpolant;

  procedure Construct_Interpolant ( p : in out Poly; i : in natural;
                                    cv : in Standard_Complex_Vectors.Vector ) is

  -- DESCRIPTION :
  --   With the coefficients of its monomials, the interpolant will be
  --   constructed.  Monomials that have degree one in x_i will be ignored.

    cnt : natural := 0;

    procedure Scan_Term ( t : in out Term; continue : out boolean ) is
    begin
      if cnt > cv'last
       then continue := false;
       else if cnt = 0 or t.dg(i) /= 1
             then if cnt >= cv'first
                   then t.cf := cv(cnt);
                  end if;
                  cnt := cnt + 1;
            end if;
            continue := true;
      end if;
    end Scan_Term;
    procedure Scan_Poly is new Changing_Iterator(Scan_Term);

  begin
    Scan_Poly(p);
  end Construct_Interpolant;

  function Interpolate ( p : Poly; b : natural; sols : Solution_List )
                       return Poly is

  -- DESCRIPTION :
  --   Constructs the interpolating polynomial with the same monomial
  --   structure as the given polynomial.

    res : Poly := Complex_Randomize(p);
    mat : Standard_Complex_Matrices.Matrix(1..b,1..b);
    rhs : Standard_Complex_Vectors.Vector(1..b);
    ipvt : Standard_Natural_Vectors.Vector(1..b);
    info : integer;

  begin
    Interpolation_System(res,b,sols,mat,rhs);
    for i in ipvt'range loop
      ipvt(i) := i;
    end loop;
    lufac(mat,b,ipvt,info);
    if info = 0
     then lusolve(mat,b,ipvt,rhs);
    end if;
    Construct_Interpolant(res,rhs);
    return res;
  end Interpolate;

  function Interpolate ( p : Poly; i,b : natural; sols : Solution_List )
                       return Poly is

  -- DESCRIPTION :
  --   Constructs the interpolating polynomial with the same monomial
  --   structure as the given polynomial.  The monomials that have degree 
  --   one in x_i will not appear in the interpolation matrix.

    res : Poly := Complex_Randomize(p);
    mat : Standard_Complex_Matrices.Matrix(1..b,1..b);
    rhs : Standard_Complex_Vectors.Vector(1..b);
    ipvt : Standard_Natural_Vectors.Vector(1..b);
    info : integer;

  begin
    Interpolation_System(res,i,b,sols,mat,rhs);
    for j in ipvt'range loop
      ipvt(j) := j;
    end loop;
    lufac(mat,b,ipvt,info);
    if info = 0
     then lusolve(mat,b,ipvt,rhs);
    end if;
    Construct_Interpolant(res,i,rhs);
    return res;
  end Interpolate;

-- TARGET ROUTINES :

  function Dense_Representation
              ( p : Poly_Sys; z : partition ) return Poly_Sys is

    d : constant Standard_Integer_Matrices.Matrix := Degree_Table(p,z);

  begin
    return Dense_Representation(p,z,d);
  end Dense_Representation;

  function Dense_Representation
              ( p : Poly_Sys; z : partition; d : Matrix ) return Poly_Sys is

    res : Poly_Sys(d'range(1));

  begin
    for i in res'range loop
      if p(i) /= Null_Poly
       then res(i) := Dense_Representation(p(i),i,z,d);
      end if;
    end loop;
    return res;
  end Dense_Representation;

  function Independent_Representation ( p : Poly_Sys ) return Poly_Sys is
 
    res : Poly_Sys(p'range);
    it : Term;

  begin
    Copy(p,res);
    for i in res'range loop
      if p(i) /= Null_Poly
       then it := Initial_Term(res(i));
            for j in res'range loop
              if j /= i and then (p(j) /= Null_Poly)
               then Eliminate_Term(res(j),it.dg);
              end if;
            end loop;
      end if;
    end loop;
    return res;
  end Independent_Representation;

  function Independent_Roots ( p : Poly_Sys ) return natural is

    ntp : natural := 0;
    min : natural := ntp;

  begin
    for i in p'first..p'last loop
      if p(i) /= Null_Poly
       then ntp := Number_of_Terms(p(i));
            if min = 0
             then min := ntp;
             elsif ntp < min
                 then min := ntp;
            end if;
      end if;
    end loop;
    if min = 0
     then return 0;
     else return (min-1);
    end if;
  end Independent_Roots;

  function Number_of_Terms ( p : Poly; i : natural ) return natural is

  -- DESCRIPTION :
  --   Returns the number of monomials of p, without those monomials that
  --   have degree one in x_i.

    cnt : natural := 0;

    procedure Count_Term ( t : in Term; continue : out boolean ) is
    begin
      if t.dg(i) /= 1
       then cnt := cnt+1;
      end if;
      continue := true;
    end Count_Term;
    procedure Count_Terms is new Visiting_Iterator(Count_Term);

  begin
    Count_Terms(p);
    return cnt;
  end Number_of_Terms;

  function Independent_Roots ( p : Poly_Sys; i : natural ) return natural is

    ntp : natural := 0;
    min : natural := ntp;

  begin
    for j in p'first..p'last loop
      if p(j) /= Null_Poly
       then ntp := Number_of_Terms(p(j),i);
            if min = 0
             then min := ntp;
             elsif ntp < min
                 then min := ntp;
            end if;
      end if;
    end loop;
    if min = 0
     then return 0;
     else return (min-1);
    end if;
  end Independent_Roots;

  function Interpolate ( p : Poly_Sys; b : natural; sols : Solution_List )
                       return Poly_Sys is

    res : Poly_Sys(p'range);

  begin 
    for i in p'range loop
      if p(i) /= Null_Poly
       then res(i) := Interpolate(p(i),b,sols);
      end if;
    end loop;
    return res;
  end Interpolate;

  function Interpolate ( p : Poly_Sys; i,b : natural; sols : Solution_List )
                       return Poly_Sys is

    res : Poly_Sys(p'range);

  begin 
    for j in p'range loop
      if p(j) /= Null_Poly
       then res(j) := Interpolate(p(j),i,b,sols);
      end if;
    end loop;
    return res;
  end Interpolate;

end Interpolating_Homotopies;