File: [local] / OpenXM_contrib / PHC / Ada / Schubert / sagbi_homotopies.adb (download)
Revision 1.1.1.1 (vendor branch), Sun Oct 29 17:45:32 2000 UTC (23 years, 8 months ago) by maekawa
Branch: PHC, MAIN
CVS Tags: v2, maekawa-ipv6, RELEASE_1_2_3, RELEASE_1_2_2_KNOPPIX_b, RELEASE_1_2_2_KNOPPIX, RELEASE_1_2_2, RELEASE_1_2_1, HEAD Changes since 1.1: +0 -0
lines
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;