Return to sagbi_homotopies.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;           use Standard_Natural_Vectors;
with Standard_Complex_Poly_Functions;    use Standard_Complex_Poly_Functions;
with Brackets;                           use Brackets;
with Bracket_Monomials;                  use Bracket_Monomials;
with Bracket_Polynomials;                use Bracket_Polynomials;
with Straightening_Syzygies;             use Straightening_Syzygies;
with Bracket_Expansions;                 use Bracket_Expansions;
with Evaluated_Minors;                   use Evaluated_Minors;

package body SAGBI_Homotopies is

  function Coordinatize_Hexadecimal ( b : Bracket ) return natural is

  -- DESCRIPTION :
  --   Returns the hexadecimal expansion of the entries in the bracket.

	res : natural := 0;

  begin
    for i in b'range loop
      res := res*16 + b(i);
    end loop;
    return res;
  end Coordinatize_Hexadecimal;

  function Unsigned ( i : integer ) return natural is

  -- DESCRIPTION :
  --   Returns the unsigned integer.

    n : natural;

  begin
    if i < 0
     then n := -i;
     else n := i;
    end if;
    return n;
  end Unsigned;

  function Bracketize_Hexadecimal ( n,d : natural ) return Bracket is

  -- DESCRIPTION :
  --   Returns the d-bracket from the hexadecimal expansion n.

	res : Bracket(1..d);
    nn : natural := n;

  begin
    for i in reverse 1..d loop
      res(i) := nn mod 16;
      nn := nn/16;
    end loop;
    return res;
  end Bracketize_Hexadecimal;

  function Coordinatize ( p : Bracket_Polynomial ) return Bracket_Polynomial is

  -- DESCRIPTION :
  --   Replaces the first bracket in every monomial by the decimal expansion.

    res : Bracket_Polynomial;

    procedure Coordinatize_Term ( t : in Bracket_Term; cont1 : out boolean ) is

      first,second : boolean;
      bm : Bracket_Monomial;
      bt : Bracket_Term;

      procedure Coordinatize_Bracket ( b : in Bracket; cont2 : out boolean ) is
      begin
        if first
         then bt.coeff := Create(double_float(Coordinatize_Hexadecimal(b)));
              first := false;
              second := true;
         elsif second
             then bm := Create(b);
             else Multiply(bm,b);
        end if;
        cont2 := true;
      end Coordinatize_Bracket;
      procedure Coordinatize_Brackets is
        new Enumerate_Brackets(Coordinatize_Bracket);

    begin
      first := true; second := false;
      Coordinatize_Brackets(t.monom);
      bt.monom := bm;
      if REAL_PART(t.coeff) < 0.0
       then Min(res,bt);
       else Add(res,bt);
      end if;
      cont1 := true;
    end Coordinatize_Term;
    procedure Coordinatize_Terms is new Enumerate_Terms(Coordinatize_Term);

  begin
    Coordinatize_Terms(p);
    return res;
  end Coordinatize;

  procedure Divide ( p : in out Poly; w : in natural ) is

  -- DESCRIPTION :
  --   Divides the polynomial by t^w.

    procedure Divide_Term ( t : in out Term; continue : out boolean ) is
    begin
      t.dg(t.dg'last) := t.dg(t.dg'last) - w;
      continue := true;
    end Divide_Term;
    procedure Divide_Terms is new Changing_Iterator(Divide_Term);

  begin
    Divide_Terms(p);
  end Divide;

  function Weight ( e : Standard_Natural_Vectors.Vector; n,d : natural )
                  return natural is

  -- DESCRIPTION :
  --   Returns the weight of the exponent vector for the localization that
  --   takes the d-by-d identitity matrix in the lower-right of the d-plane.
  --   The lifting recipe is xij*t^((i-1)*(d-j)).

    res : natural := 0;
    jmp,ind : natural;

  begin
    for j in 1..d loop
      jmp := d-j;
      for i in 1..n-d loop
        ind := (i-1)*d + j;
        if e(ind) > 0
         then res := res + (i-1)*jmp;
        end if;
      end loop;
    end loop;
    return res;
  end Weight;

  function Weight ( locmap : Standard_Natural_Matrices.Matrix;
                    e : Standard_Natural_Vectors.Vector ) return natural is

  -- DESCRIPTION :
  --   Returns the weight of the exponent vector as xij*t^((i-1)*(d-j))
  --   for the localization pattern in locmap.

    res : natural := 0;
    d : constant natural := locmap'length(2);
    jmp,ind : natural;

  begin
    ind := 0;
    for i in locmap'range(1) loop
      for j in locmap'range(2) loop
        jmp := d-j;
        if locmap(i,j) = 2
         then ind := ind+1;
              if e(ind) > 0
               then res := res + (i-1)*jmp;
              end if;
        end if;
      end loop;
    end loop;
    return res;
  end Weight;

  function Lift ( p : Poly; n,d : natural ) return Poly is

  -- DESCRIPTION :
  --   Returns the lifted polynomial, where the xij is lifted according
  --   to xij*t^((i-1)*(d-j)).  The lowest powers of t are divided out.
  --   The d-by-d identity matrix is the lower-right of the d-plane.

    res : Poly := Null_Poly;
    first : boolean := true;
    minwei : natural;

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

      tt : Term;
      wei : natural;

    begin
      tt.cf := t.cf;
      tt.dg := new Standard_Natural_Vectors.Vector(1..t.dg'last+1);
      tt.dg(t.dg'range) := t.dg.all;
      wei := Weight(t.dg.all,n,d);
      tt.dg(tt.dg'last) := wei;
      Add(res,tt);
      Clear(tt.dg);
      if first
       then minwei := wei;
            first := false;
       elsif wei < minwei
           then minwei := wei;
      end if;
      continue := true;
    end Lift_Term;
    procedure Lift_Terms is new Visiting_Iterator(Lift_Term);

  begin
    Lift_Terms(p);
    if minwei /= 0
     then Divide(res,minwei);
    end if;
    return res;
  end Lift;

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

  -- DESCRIPTION :
  --   Lifts p as to xij*t^((i-1)*(d-j)) and divides by the lowest powers
  --   of t, respecting the localization pattern in locmap.

    res : Poly := Null_Poly;
    first : boolean := true;
    minwei : natural;

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

      tt : Term;
      wei : natural;

    begin
      tt.cf := t.cf;
      tt.dg := new Standard_Natural_Vectors.Vector(1..t.dg'last+1);
      tt.dg(t.dg'range) := t.dg.all;
      wei := Weight(locmap,t.dg.all);
      tt.dg(tt.dg'last) := wei;
      Add(res,tt);
      Clear(tt.dg);
      if first
       then minwei := wei;
            first := false;
       elsif wei < minwei
           then minwei := wei;
      end if;
      continue := true;
    end Lift_Term;
    procedure Lift_Terms is new Visiting_Iterator(Lift_Term);

  begin
    Lift_Terms(p);
    if minwei /= 0
     then Divide(res,minwei);
    end if;
    return res;
  end Lift;

-- TARGET ROUTINES :

  function Lifted_Localized_Laplace_Expansion ( n,d : natural ) return Poly is

    res : Poly := Null_Poly;
    p : Bracket_Polynomial := Laplace_Expansion(n,n-d);

    procedure Visit_Term ( t : in Bracket_Term; continue : out boolean ) is

      first : boolean := true;
      cf : integer;

      procedure Visit_Bracket ( b : in Bracket; cont : out boolean ) is

        pb,lp : Poly;

      begin
        if first
         then cf := Coordinatize_Hexadecimal(b);
              first := false;
         else pb := Localized_Expand(n,d,b);
              lp := Lift(pb,n,d);  Clear(pb);
              Mul(lp,Create(double_float(cf)));
              Add(res,lp);
              Clear(lp);
        end if;
        cont := true;
      end Visit_Bracket;
      procedure Visit_Brackets is new Enumerate_Brackets(Visit_Bracket);

    begin
      Visit_Brackets(t.monom);
      continue := true;
    end Visit_Term;
    procedure Visit_Terms is new Enumerate_Terms(Visit_Term);

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

  function Lifted_Localized_Laplace_Expansion
             ( locmap : Standard_Natural_Matrices.Matrix ) return Poly is

    res : Poly := Null_Poly;
    n : constant natural := locmap'length(1);
    d : constant natural := locmap'length(2);
    p : Bracket_Polynomial := Laplace_Expansion(n,n-d);

    procedure Visit_Term ( t : in Bracket_Term; continue : out boolean ) is

      first : boolean := true;
      cf : integer;

      procedure Visit_Bracket ( b : in Bracket; cont : out boolean ) is

        pb,lp : Poly;

      begin
        if first
         then cf := Coordinatize_Hexadecimal(b);
              first := false;
         else pb := Expand(locmap,b);
              Reduce_Variables(locmap,pb);
              lp := Lift(locmap,pb);  Clear(pb);
              Mul(lp,Create(double_float(cf)));
              Add(res,lp);
              Clear(lp);
        end if;
        cont := true;
      end Visit_Bracket;
      procedure Visit_Brackets is new Enumerate_Brackets(Visit_Bracket);

    begin
      Visit_Brackets(t.monom);
      continue := true;
    end Visit_Term;
    procedure Visit_Terms is new Enumerate_Terms(Visit_Term);

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

  function Intersection_Coefficients
              ( m : Standard_Floating_Matrices.Matrix;
                c : Standard_Complex_Vectors.Vector )
              return Standard_Complex_Vectors.Vector is

    res : Standard_Complex_Vectors.Vector(c'range);
    nmd : constant natural := m'last(2);
    ind : integer;
    b : Bracket(1..nmd);

  begin
    for i in c'range loop
      ind := integer(REAL_PART(c(i)));
      b := Bracketize_Hexadecimal(Unsigned(ind),nmd);
      if ind > 0
       then res(i) := Create(Determinant(m,b));
       else res(i) := Create(-Determinant(m,b));
      end if;
    end loop;
    return res;
  end Intersection_Coefficients;

  function Intersection_Coefficients
              ( m : Standard_Complex_Matrices.Matrix;
                c : Standard_Complex_Vectors.Vector )
              return Standard_Complex_Vectors.Vector is

    res : Standard_Complex_Vectors.Vector(c'range);
    nmd : constant natural := m'last(2);
    ind : integer;
    b : Bracket(1..nmd);

  begin
    for i in c'range loop
      ind := integer(REAL_PART(c(i)));
      b := Bracketize_Hexadecimal(Unsigned(ind),nmd);
      if ind > 0
       then res(i) := Determinant(m,b);
       else res(i) := -Determinant(m,b);
      end if;
    end loop;
    return res;
  end Intersection_Coefficients;

  function Intersection_Condition
             ( m : Standard_Floating_Matrices.Matrix; p : Poly ) return Poly is

    res : Poly := Null_Poly;
    nmd : constant natural := m'last(2);

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

      c : integer := integer(REAL_PART(t.cf));
      b : Bracket(1..nmd) := Bracketize_Hexadecimal(Unsigned(c),nmd);
      det : double_float := Determinant(m,b);
      rt : Term;

    begin
      if c > 0
       then rt.cf := Create(det);
       else rt.cf := Create(-det);
      end if;
      rt.dg := t.dg;
      Add(res,rt);
      continue := true;
    end Visit_Term;
    procedure Visit_Terms is new Visiting_Iterator(Visit_Term);

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

  function Intersection_Condition
             ( m : Standard_Complex_Matrices.Matrix; p : Poly ) return Poly is

    res : Poly := Null_Poly;
    nmd : constant natural := m'last(2);

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

      c : integer := integer(REAL_PART(t.cf));
      b : Bracket(1..nmd) := Bracketize_Hexadecimal(Unsigned(c),nmd);
      det : Complex_Number := Determinant(m,b);
      rt : Term;

    begin
      if c > 0
       then rt.cf :=  det;
       else rt.cf := -det;
      end if;
      rt.dg := t.dg;
      Add(res,rt);
      continue := true;
    end Visit_Term;
    procedure Visit_Terms is new Visiting_Iterator(Visit_Term);

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

end SAGBI_Homotopies;