File: [local] / OpenXM_contrib / PHC / Ada / Schubert / symbolic_minor_equations.adb (download)
Revision 1.1.1.1 (vendor branch), Sun Oct 29 17:45:33 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;
with Matrix_Indeterminates;
with Bracket_Polynomials; use Bracket_Polynomials;
with Straightening_Syzygies; use Straightening_Syzygies;
package body Symbolic_Minor_Equations is
-- AUXILIARIES TO Minor_Equations :
function Substitute ( b,minor : Bracket ) return Bracket is
-- DESCRIPTION :
-- The ith entry in the bracket b is replaced by minor(b(i)).
res : Bracket(b'range);
start : integer;
begin
if b(b'first) = 0
then start := b'first+1; -- bracket is a coefficient bracket
res(res'first) := 0; -- result is also a coefficient bracket
else start := b'first; -- second bracket in expansion
end if;
for i in start..b'last loop
res(i) := minor(b(i));
end loop;
return res;
end Substitute;
function Substitute ( bm : Bracket_Monomial; minor : Bracket )
return Bracket_Monomial is
-- DESCRIPTION :
-- Substitutes the entries in the brackets according to the minor.
-- REQUIRED : bm is a quadratic monomial.
res : Bracket_Monomial;
lb : Link_to_Bracket;
first : boolean := true;
procedure Substitute_Bracket ( b : in Bracket; continue : out boolean ) is
sb : Bracket(b'range) := Substitute(b,minor);
begin
if first
then lb := new Bracket'(sb); -- save to preserve order
first := false;
else res := sb*lb.all;
end if;
continue := true;
end Substitute_Bracket;
procedure Substitute_Brackets is
new Enumerate_Brackets(Substitute_Bracket);
begin
Substitute_Brackets(bm);
Clear(lb);
return res;
end Substitute;
function Substitute ( bt : Bracket_Term; minor : Bracket )
return Bracket_Term is
-- DESCRIPTION :
-- The ith entry in every bracket of the term according to the minor.
res : Bracket_Term;
begin
res.coeff := bt.coeff;
res.monom := Substitute(bt.monom,minor);
return res;
end Substitute;
function Substitute ( bp : Bracket_Polynomial; minor : Bracket )
return Bracket_Polynomial is
-- DESCRIPTION :
-- The labels in the brackets of bp are replaced by those in the minor.
-- The bracket polynomial represents the Laplace expansion of that minor.
res : Bracket_Polynomial;
procedure Substitute_Term ( t : in Bracket_Term; continue : out boolean ) is
st : Bracket_Term := Substitute(t,minor);
begin
Add(res,st);
Clear(st);
continue := true;
end Substitute_Term;
procedure Substitute_Terms is new Enumerate_Terms(Substitute_Term);
begin
Substitute_Terms(bp);
return res;
end Substitute;
-- AUXILIARIES TO Expanded_Minor :
function Subtract ( b : Bracket; i : natural ) return Bracket is
-- DESCRIPTION :
-- Returns a smaller bracket with the ith entry removed.
res : Bracket(b'first..b'last-1);
begin
res(b'first..(i-1)) := b(b'first..(i-1));
res(i..res'last) := b((i+1)..b'last);
return res;
end Subtract;
function General_Expanded_Minor ( m : Standard_Complex_Poly_Matrices.Matrix;
b : Bracket ) return Poly is
-- DESCRIPTION :
-- This function treats the case for Expanded_Minor when b'length > 2.
res : Poly := Null_Poly;
sig : integer;
submin,acc : Poly;
begin
if b'last mod 2 = 0
then sig := -1;
else sig := +1;
end if;
for i in b'range loop
if m(b(i),b'last) /= Null_Poly
then submin := Expanded_Minor(m,Subtract(b,i));
if submin /= Null_Poly
then acc := m(b(i),b'last)*submin;
if sig > 0
then Add(res,acc);
else Sub(res,acc);
end if;
Clear(acc);
end if;
Clear(submin);
end if;
sig := -sig;
end loop;
return res;
end General_Expanded_Minor;
procedure Purify ( p : in out Poly ) is
-- DESCRIPTION :
-- Eliminates terms that have coefficients less that 10.0**(-10).
tol : constant double_float := 10.0**(-10);
res : Poly := Null_Poly;
procedure Scan_Term ( t : in Term; cont : out boolean ) is
begin
if AbsVal(t.cf) > tol
then Add(res,t);
end if;
cont := true;
end Scan_Term;
procedure Scan_Terms is new Visiting_Iterator(Scan_Term);
begin
Scan_Terms(p);
Clear(p);
p := res;
end Purify;
-- TARGET ROUTINES :
function Schubert_Pattern ( n : natural; b1,b2 : Bracket )
return Standard_Complex_Poly_Matrices.Matrix is
res : Standard_Complex_Poly_Matrices.Matrix(1..n,b1'range);
d : constant natural := b1'last;
begin
for i in 1..n loop
for j in b1'range loop
if ((i < b2(j)) or (i > n+1-b1(b1'last+1-j)))
then res(i,j) := Null_Poly;
else res(i,j) := Matrix_Indeterminates.Monomial(n,d,i,j);
end if;
end loop;
end loop;
return res;
end Schubert_Pattern;
function Localization_Pattern ( n : natural; top,bottom : Bracket )
return Standard_Complex_Poly_Matrices.Matrix is
p : constant natural := top'length;
res : Standard_Complex_Poly_Matrices.Matrix(1..n,1..p);
begin
for i in res'range(1) loop
for j in res'range(2) loop
if i < top(j) or i > bottom(j)
then res(i,j) := Null_Poly;
else res(i,j) := Matrix_Indeterminates.Monomial(n,p,i,j);
end if;
end loop;
end loop;
return res;
end Localization_Pattern;
-- SYMBOLIC REPRESENTATION OF THE EQUATIONS :
function Maximal_Minors ( n,m : natural ) return Bracket_Monomial is
res : Bracket_Monomial;
first : boolean := true;
accu : Bracket(1..m);
procedure Enumerate_Brackets ( k,start : natural ) is
-- DESCRIPTION :
-- Enumerate all m-brackets with entries between 1 and n,
-- beginning at start, and currently filling in the kth entry.
begin
if k > m
then if first
then res := Create(accu); first := false;
else Multiply(res,accu);
end if;
else for i in start..n loop
accu(k) := i;
Enumerate_Brackets(k+1,i+1);
end loop;
end if;
end Enumerate_Brackets;
begin
Enumerate_Brackets(1,1);
return res;
end Maximal_Minors;
function Minor_Equations
( m,d : natural; bm : Bracket_Monomial ) return Bracket_System is
res : Bracket_System(0..Number_of_Brackets(bm));
gen : Bracket_Polynomial := Laplace_Expansion(m,d);
cnt : natural := 0;
procedure Substitute_Bracket ( b : in Bracket; continue : out boolean ) is
begin
cnt := cnt + 1;
res(cnt) := Substitute(gen,b);
continue := true;
end Substitute_Bracket;
procedure Substitute_Brackets is
new Enumerate_Brackets(Substitute_Bracket);
begin
res(0) := gen;
Substitute_Brackets(bm);
return res;
end Minor_Equations;
function Expanded_Minor ( m : Standard_Complex_Poly_Matrices.Matrix;
b : Bracket ) return Poly is
res : Poly := Null_Poly;
acc : Poly;
begin
if b'length = 1
then Copy(m(b(1),1),res);
elsif b'length = 2
then if (m(b(1),1) /= Null_Poly) and (m(b(2),2) /= Null_Poly)
then res := m(b(1),1)*m(b(2),2);
end if;
if (m(b(2),1) /= Null_Poly) and (m(b(1),2) /= Null_Poly)
then acc := m(b(2),1)*m(b(1),2);
Sub(res,acc);
Clear(acc);
end if;
else res := General_Expanded_Minor(m,b);
end if;
Purify(res);
return res;
end Expanded_Minor;
function Extend_Zero_Lifting ( p : Poly ) return Poly is
res : Poly := Null_Poly;
procedure Extend_Term ( t : in Term; continue : out boolean ) is
et : Term;
begin
et.dg := new Standard_Natural_Vectors.Vector(t.dg'first..t.dg'last+1);
et.dg(t.dg'range) := t.dg.all;
et.dg(et.dg'last) := 0;
et.cf := t.cf;
Add(res,et);
Clear(et.dg);
continue := true;
end Extend_Term;
procedure Extend_Terms is new Visiting_Iterator(Extend_Term);
begin
Extend_Terms(p);
return res;
end Extend_Zero_Lifting;
end Symbolic_Minor_Equations;