Return to determinantal_systems.adb CVS log

Up to [local] / OpenXM_contrib / PHC / Ada / Schubert

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_Polynomials;       use Standard_Complex_Polynomials; 
with Bracket_Monomials;                  use Bracket_Monomials;
with Bracket_Systems;                    use Bracket_Systems;
with Evaluated_Minors;                   use Evaluated_Minors;
with Symbolic_Minor_Equations;           use Symbolic_Minor_Equations;
with Numeric_Minor_Equations;            use Numeric_Minor_Equations;

package body Determinantal_Systems is

-- AUXILIARY TO EVALUATION OF MAXIMAL MINORS :

  function Number_of_Maximal_Minors ( nrows,ncols : natural ) return natural is

  -- DESCRIPTION :
  --   Returns the number of maximal minors of a matrix with nrows and ncols.

  -- REQUIRED : nrows > ncols.

  begin
    if ncols = 1
     then return nrows;
     else return Number_of_Maximal_Minors(nrows-1,ncols-1)*nrows/ncols;
    end if;
  end Number_of_Maximal_Minors;

-- LOCALIZATION MAPS :

  function Localize ( locmap : Standard_Natural_Matrices.Matrix;
                      t : Term ) return Term is

  -- DESCRIPTION :
  --   Applies the localization map to the term, eliminating those xij's
  --   xij for which the corresponding entry in locmap is either 0 or 1.

  -- NOTE : 
  --   This localization assumes that t.dg(k) = 0 with k for which the
  --   corresponding (i,j) with locmap(i,j) = 0.

    res : Term;
    ndg : Standard_Natural_Vectors.Vector(t.dg'range);
    cnt : natural := t.dg'first-1;
    ind : natural := cnt;

  begin
    for i in locmap'range(1) loop       -- lexicographic order of variables
      for j in locmap'range(2) loop
        ind := ind+1;
        if locmap(i,j) = 2
         then cnt := cnt + 1;
              ndg(cnt) := t.dg(ind);
        end if;
      end loop;
    end loop;
    for i in ind+1..t.dg'last loop      -- leave the lifting !
      cnt := cnt+1;
      ndg(cnt) := t.dg(i);
    end loop;
    res.cf := t.cf;
    res.dg := new Standard_Natural_Vectors.Vector'(ndg(1..cnt));
    return res;
  end Localize;

  function Localize ( locmap : Standard_Natural_Matrices.Matrix;
                      p : Poly ) return Poly is

  -- DESCRIPTION :
  --   Applies the localization map to the polynomial, eliminating
  --   those xij's for which locmap(i,j) is either 0 or 1.

    res : Poly := Null_Poly;

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

      lt : Term := Localize(locmap,t);

    begin
      Add(res,lt);
      Clear(lt.dg);
      continue := true;
    end Localize_Term;
    procedure Localize_Terms is new Visiting_Iterator(Localize_Term);

  begin
    Localize_Terms(p);
    return res;
  end Localize;

-- TARGET ROUTINES :

  function Standard_Coordinate_Frame
             ( x : Standard_Complex_Poly_Matrices.Matrix;
               plane : Standard_Complex_Matrices.Matrix )
             return Standard_Natural_Matrices.Matrix is

    res : Standard_Natural_Matrices.Matrix(x'range(1),x'range(2));
    tol : constant double_float := 10.0**(-10);
    first : boolean;

  begin
    for j in res'range(2) loop
      first := true;
      for i in res'range(1) loop
        if x(i,j) = Null_Poly
         then res(i,j) := 0;
         elsif (first and (AbsVal(plane(i,j)) > tol))
             then res(i,j) := 1; first := false;
             else res(i,j) := 2;
        end if;
      end loop;
    end loop;
    return res;
  end Standard_Coordinate_Frame;

  function Maximal_Coordinate_Frame
             ( x : Standard_Complex_Poly_Matrices.Matrix;
               plane : Standard_Complex_Matrices.Matrix )
             return Standard_Natural_Matrices.Matrix is

    res : Standard_Natural_Matrices.Matrix(x'range(1),x'range(2));
    max,tmp : double_float;
    ind : natural;

  begin
    for j in res'range(2) loop
      for i in res'range(1) loop
        if x(i,j) = Null_Poly
         then res(i,j) := 0;
         else res(i,j) := 2;
        end if;
      end loop;
      max := 0.0; ind := 0;
      for i in res'range(1) loop
        tmp := AbsVal(plane(i,j));
        if tmp > max
         then max := tmp; ind := i;
        end if;
      end loop;
      res(ind,j) := 1;
    end loop;
    return res;
  end Maximal_Coordinate_Frame;

  function Localize ( locmap : Standard_Natural_Matrices.Matrix;
                      p : Poly_Sys ) return Poly_Sys is

    res : Poly_Sys(p'range);

  begin
    for i in p'range loop
      res(i) := Localize(locmap,p(i));
    end loop;
    return res;
  end Localize;

-- CONSTRUCT THE POLYNOMIAL EQUATIONS :

  function Polynomial_Equations
              ( l : Standard_Complex_Matrices.Matrix;
                x : Standard_Complex_Poly_Matrices.Matrix ) return Poly_Sys is

    n : constant natural := x'length(1);
    p : constant natural := x'length(2);
    kd : constant natural := p + l'length(2);
    bm : Bracket_Monomial := Maximal_Minors(n,kd);
    bs : Bracket_System(0..Number_of_Brackets(bm))
       := Minor_Equations(kd,kd-p,bm);
    res : Poly_Sys(bs'first+1..bs'last) := Expanded_Minors(l,x,bs);

  begin
    Clear(bm); Clear(bs);
    return res;
  end Polynomial_Equations;

  procedure Concat ( l : in out Link_to_Poly_Sys; p : Poly_Sys ) is
  begin
    if l = null
     then declare
            newsys : Poly_Sys(p'range);
            cnt : natural := p'first-1;
          begin
            for i in p'range loop
              if p(i) /= Null_Poly
               then cnt := cnt+1;
                    newsys(cnt) := p(i);
              end if;
            end loop;
            l := new Poly_Sys'(newsys(p'first..cnt));
          end;
     else declare
            newsys : Poly_Sys(l'first..l'last+p'length);
            ind : natural := l'last;
          begin
            newsys(l'range) := l(l'range);
            for i in p'range loop
              if p(i) /= Null_Poly
               then ind := ind+1;
                    newsys(ind) := p(i);
              end if;
            end loop;
            l := new Poly_Sys'(newsys(l'first..ind));
          end;
    end if;
  end Concat;

  function Polynomial_Equations
              ( l : Standard_Complex_VecMats.VecMat;
                x : Standard_Complex_Poly_Matrices.Matrix ) return Poly_Sys is

    n : constant natural := x'length(1);
    p : constant natural := x'length(2);
    kd : constant natural := p + l(l'first)'length(2);
    bm : Bracket_Monomial := Maximal_Minors(n,kd);
    bs : Bracket_System(0..Number_of_Brackets(bm))
       := Minor_Equations(kd,kd-p,bm);
    nb : constant natural := l'length;
    res : Poly_Sys(bs'first+1..nb*bs'last);
    cnt : natural := bs'first;

  begin
    for i in 1..nb loop
      declare 
        sys : constant Poly_Sys := Expanded_Minors(l(i).all,x,bs);
      begin
        for j in sys'range loop
          cnt := cnt + 1;
          res(cnt) := sys(j);
        end loop;
      end;
    end loop;
    Clear(bm); Clear(bs);
    return res;
  end Polynomial_Equations;

-- EVALUATORS AND DIFFERENTIATORS :

  function Eval ( l,x : Standard_Complex_Matrices.Matrix )
                return Complex_Number is

    wrk : Standard_Complex_Matrices.Matrix(x'range(1),x'range(1));

  begin
    for j in l'range(2) loop
      for i in l'range(1) loop
        wrk(i,j) := l(i,j);
      end loop;
    end loop;
    for j in x'range(2) loop
      for i in x'range(1) loop
        wrk(i,l'last(2)+j) := x(i,j);
      end loop;
    end loop;
    return Determinant(wrk);
  end Eval;

  function Eval ( l,x : Standard_Complex_Matrices.Matrix ) return Vector is

    n : constant natural := l'length(2) + x'length(2);
    dim : constant natural := Number_of_Maximal_Minors(l'length(1),n);
    res : Standard_Complex_Vectors.Vector(1..dim);
      -- := (1..dim => Create(0.0));
    cnt : natural := 0;
    rws : Standard_Natural_Vectors.Vector(1..n);
    wrk : Standard_Complex_Matrices.Matrix(1..n,1..n);

    procedure Choose_next_Row ( k,start : in natural ) is

    -- DESCRIPTION :
    --   Chooses next k-th row within the range start..l'last(1), if k <= n.
    --   If k > n, then the minor defined by the current row selection
   --    is computed and added to the result.

    begin
      if k > n
       then for j in l'range(2) loop              -- select rows
              for i in 1..n loop
                wrk(i,j) := l(rws(i),j);
              end loop;
            end loop;
            for j in x'range(2) loop
              for i in 1..n loop
                wrk(i,l'last(2)+j) := x(rws(i),j);
              end loop;
            end loop;
            cnt := cnt + 1;
            res(cnt) := Determinant(wrk);
       else for i in start..l'last(1) loop
              rws(k) := i;
              Choose_next_Row(k+1,i+1);
            end loop;
      end if;
    end Choose_next_Row;

  begin
    Choose_next_Row(1,1);
    return res;
  end Eval;

  function Minor ( l,x : Standard_Complex_Matrices.Matrix; row,col : natural )
                 return Complex_Number is

  -- DESCRIPTION :
  --   Returns the minor [l|x] with row and column removed and with
  --   the appropriate sign.

    wrk : Standard_Complex_Matrices.Matrix
            (x'first(1)..x'last(1)-1,x'first(1)..x'last(1)-1);
    ind : natural;

  begin
    for j in l'range(2) loop
      for i in l'range(1) loop
        if i < row
         then wrk(i,j) := l(i,j);
         elsif i > row
             then wrk(i-1,j) := l(i,j);
        end if;
      end loop;
    end loop;
    for j in x'range(2) loop
      if j < col
       then ind := l'last(2) + j;
       elsif j > col
           then ind := l'last(2) + j - 1;
      end if;
      if j /= col
       then for i in x'range(1) loop
              if i < row
               then wrk(i,ind) := x(i,j);
               elsif i > row
                   then wrk(i-1,ind) := x(i,j);
              end if;
            end loop;
      end if;
    end loop;
    if (row + l'last(2) + col) mod 2 = 0
	 then return Determinant(wrk);
	 else return -Determinant(wrk);
    end if;
  end Minor;

  function Diff ( l,x : Standard_Complex_Matrices.Matrix; i : natural )
                return Complex_Number is

    p : constant natural := x'length(2);
    row,col : natural;

  begin
    row := i/p + 1;
    col := i mod p;
    if col = 0
     then col := p;
          row := row - 1;
    end if;
    return Minor(l,x,row,col);
  end Diff;

  function Diff ( l,x : Standard_Complex_Matrices.Matrix;
                  locmap : Standard_Natural_Matrices.Matrix; i : natural )
                return Complex_Number is

    row,col : natural;
    cnt : natural := 0;

  begin
    for k1 in locmap'range(1) loop
      for k2 in locmap'range(2) loop
        if locmap(k1,k2) = 2
         then cnt := cnt+1;
              if cnt = i
               then row := k1;
                    col := k2;
              end if;
        end if;
        exit when (cnt = i);
      end loop;
      exit when (cnt = i);
    end loop;
    return Minor(l,x,row,col);
  end Diff;

  function Eval ( l : Standard_Complex_VecMats.VecMat;
                  x : Standard_Complex_Matrices.Matrix )
                return Standard_Complex_Vectors.Vector is

    res : Standard_Complex_Vectors.Vector(l'range);

  begin
    for i in l'range loop
      res(i) := Eval(l(i).all,x);
    end loop;
    return res;
  end Eval;

  function Diff ( l : Standard_Complex_VecMats.VecMat;
                  x : Standard_Complex_Matrices.Matrix )
                return Standard_Complex_Matrices.Matrix is

    res : Standard_Complex_Matrices.Matrix(l'range,1..x'last(1)*x'last(2));

  begin
    for i in res'range(1) loop
      for j in res'range(2) loop
        res(i,j) := Diff(l(i).all,x,j);
      end loop;
    end loop;
    return res;
  end Diff;

  function Old_Diff ( l : Standard_Complex_VecMats.VecMat;
                      x : Standard_Complex_Matrices.Matrix; nvars : natural;
                      locmap : Standard_Natural_Matrices.Matrix )
                    return Standard_Complex_Matrices.Matrix is

    res : Standard_Complex_Matrices.Matrix(l'range,1..nvars);

  begin
    for i in res'range(1) loop
      for j in res'range(2) loop
        res(i,j) := Diff(l(i).all,x,locmap,j);
      end loop;
    end loop;
    return res;
  end Old_Diff;

  function Diff ( l : Standard_Complex_VecMats.VecMat;
                  x : Standard_Complex_Matrices.Matrix; nvars : natural;
                  locmap : Standard_Natural_Matrices.Matrix )
                return Standard_Complex_Matrices.Matrix is

  -- NOTE :
  --   This Diff is organized according to the localization map,
  --   to avoid multiple searches.

    res : Standard_Complex_Matrices.Matrix(l'range,1..nvars);
    ind : natural := 0;

  begin
    for i in locmap'range loop
      for j in locmap'range loop
        if locmap(i,j) = 2
         then ind := ind+1;
              for k in l'range loop
                res(k,ind) := Minor(l(k).all,x,i,j);
              end loop;
        end if;
      end loop;
    end loop;
    return res;
  end Diff;

end Determinantal_Systems;