File: [local] / OpenXM_contrib / PHC / Ada / Schubert / bracket_expansions.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_Complex_Numbers; use Standard_Complex_Numbers;
with Standard_Natural_Vectors; use Standard_Natural_Vectors;
package body Bracket_Expansions is
-- AUXILIARIES :
function Create_Term ( n,d,i,j : natural ) return Term is
-- DESCRIPTION :
-- Returns the term that corresponds with xij.
res : Term;
begin
res.cf := Create(1.0);
res.dg := new Standard_Natural_Vectors.Vector'(1..d*n => 0);
res.dg((i-1)*d + j) := 1;
return res;
end Create_Term;
function Create_Localized_Term ( n,d,i,j,k : natural ) return Term is
-- DESCRIPTION :
-- Creates a term in localized variables:
-- if i >= k, then xij equals either 1 if i=j, or 0 if i/=j.
-- If i >= k, then the `dg' of the term on return is null.
res : Term;
begin
if i < k
then res.cf := Create(1.0);
res.dg := new Standard_Natural_Vectors.Vector'(1..d*(n-d) => 0);
res.dg((i-1)*d + j) := 1;
else if i = j+k-1
then res.cf := Create(1.0);
else res.cf := Create(0.0);
end if;
return res;
end if;
return res;
end Create_Localized_Term;
function Create_Term ( loc,n,d,i,j : natural ) return Term is
-- DESCRIPTION :
-- Returns the term that corresponds with xij, with respect to
-- the localization information in loc, which is 0,1, or 2.
res : Term;
begin
if loc = 0
then res.cf := Create(0.0);
else res.cf := Create(1.0);
end if;
res.dg := new Standard_Natural_Vectors.Vector'(1..d*n => 0);
if loc = 2
then res.dg((i-1)*d + j) := 1;
end if;
return res;
end Create_Term;
-- TARGET ROUTINES :
function Expand2 ( n,d : natural; b : Bracket ) return Poly is
-- DESCRIPTION :
-- Does the expansion in the two-dimensional case.
res,subtmp : Poly;
b11 : Term := Create_Term(n,d,b(1),1);
b12 : Term := Create_Term(n,d,b(1),2);
b21 : Term := Create_Term(n,d,b(2),1);
b22 : Term := Create_Term(n,d,b(2),2);
d1 : Poly := Create(b11);
d2 : Poly := Create(b21);
begin -- res := b11*b22 - b21*b12;
res := b22*d1;
subtmp := b12*d2;
Sub(res,subtmp);
Clear(subtmp);
Clear(b22); Clear(b12);
Clear(b11); Clear(b21);
Clear(d1); Clear(d2);
return res;
end Expand2;
function Expand ( n,d : natural; b : Bracket ) return Poly is
res : Poly;
begin
if b'last = 2
then res := Expand2(n,d,b);
else res := Null_Poly;
declare
sig : integer := +1;
b1 : Bracket(1..b'last-1);
begin
b1(1..b'last-1) := b(1..b'last-1);
if b'last mod 2 = 0
then sig := -1;
end if;
for i in b'range loop
declare
xt : Term := Create_Term(n,d,b(i),b'last);
expb1,xtexpb1 : Poly;
begin
b1(1..i-1) := b(1..i-1);
b1(i..b1'last) := b(i+1..b'last);
expb1 := Expand(n,d,b1);
xtexpb1 := xt*expb1;
if sig > 0
then Add(res,xtexpb1);
else Sub(res,xtexpb1);
end if;
sig := -sig;
Clear(expb1); Clear(xt); Clear(xtexpb1);
end;
end loop;
end;
end if;
return res;
end Expand;
function Expand ( n,d : natural; b : Bracket_Monomial ) return Poly is
res : Poly := Null_Poly;
fst : boolean := true;
procedure Expand_Bracket ( bb : in Bracket; continue : out boolean ) is
expbb : Poly;
begin
if fst
then res := Expand(n,d,bb);
fst := false;
else expbb := Expand(n,d,bb);
Mul(res,expbb);
Clear(expbb);
end if;
continue := true;
end Expand_Bracket;
procedure Expand_Brackets is new Enumerate_Brackets(Expand_Bracket);
begin
Expand_Brackets(b);
return res;
end Expand;
function Expand ( n,d : natural; b : Bracket_Term ) return Poly is
res : Poly := Expand(n,d,b.monom);
begin
Mul(res,b.coeff);
return res;
end Expand;
function Expand ( n,d : natural; b : Bracket_Polynomial ) return Poly is
res : Poly := Null_Poly;
procedure Expand_Term ( t : in Bracket_Term; continue : out boolean ) is
expt : Poly := Expand(n,d,t);
begin
Add(res,expt);
Clear(expt);
continue := true;
end Expand_Term;
procedure Expand_Terms is new Enumerate_Terms(Expand_Term);
begin
Expand_Terms(b);
return res;
end Expand;
function Localized_Expand2
( n,d : natural; b : Bracket; k : natural ) return Poly is
-- DESCRIPTION :
-- Computes the localized expansion for two dimensions.
-- Exploits the fact that b(1) < b(2).
res : Poly;
b11 : Term := Create_Localized_Term(n,d,b(1),1,k);
b12 : Term := Create_Localized_Term(n,d,b(1),2,k);
b21 : Term := Create_Localized_Term(n,d,b(2),1,k);
b22 : Term := Create_Localized_Term(n,d,b(2),2,k);
begin
if b(2) < k
then declare
subtmp : Poly;
d1 : Poly := Create(b11);
d2 : Poly := Create(b21);
begin -- res := b11*b22 - b21*b12;
res := b22*d1;
subtmp := b12*d2;
Sub(res,subtmp);
Clear(subtmp);
Clear(b22); Clear(b12);
Clear(b11); Clear(b21);
Clear(d1); Clear(d2);
end;
else if b(1) < k
then b11.cf := b11.cf*b22.cf;
b12.cf := b12.cf*b21.cf;
res := Create(b11);
Sub(res,b12);
Clear(b11); Clear(b12);
else declare
rt : Term;
dd : constant natural := d*(n-d);
begin
rt.cf := b11.cf*b22.cf - b21.cf*b12.cf;
rt.dg := new Standard_Natural_Vectors.Vector'(1..dd => 0);
res := Create(rt);
end;
end if;
end if;
return res;
end Localized_Expand2;
function Localized_Expand ( n,d : natural; b : Bracket ) return Poly is
res : Poly;
k : constant natural := n-d+1;
begin
if b'last = 2
then res := Localized_Expand2(n,d,b,k);
else res := Null_Poly;
declare
sig : integer := +1;
b1 : Bracket(1..b'last-1);
begin
b1(1..b'last-1) := b(1..b'last-1);
if b'last mod 2 = 0
then sig := -1;
end if;
for i in b'range loop
declare
xt : Term := Create_Localized_Term(n,d,b(i),b'last,k);
expb1,xtexpb1 : Poly;
begin
if xt.cf /= Create(0.0)
then b1(1..i-1) := b(1..i-1);
b1(i..b1'last) := b(i+1..b'last);
expb1 := Localized_Expand(n,d,b1);
xtexpb1 := xt*expb1;
if sig > 0
then Add(res,xtexpb1);
else Sub(res,xtexpb1);
end if;
Clear(expb1); Clear(xt); Clear(xtexpb1);
end if;
sig := -sig;
end;
end loop;
end;
end if;
return res;
end Localized_Expand;
function Localization_Map ( n,d : natural ) return Matrix is
res : Matrix(1..n,1..d);
row,col : natural;
begin
for i in 1..n loop -- initialization
for j in 1..d loop
res(i,j) := -1; -- means empty space
end loop;
end loop;
col := 0;
for i in 1..n loop -- one free element in every row
col := col+1;
if col > d
then col := 1;
end if;
res(i,col) := 2;
end loop;
row := 0;
for j in 1..d loop -- one 1 in every column
loop
row := row+1;
if row > n
then row := 1;
end if;
exit when (res(row,j) = -1); -- empty space found
end loop;
res(row,j) := 1;
end loop;
for j in 1..d loop -- fill in d-1 zeros in every column
for i in 1..(d-1) loop
row := 0;
loop
row := row+1;
if row > n
then row := 1;
end if;
exit when (res(row,j) = -1); -- empty space found
end loop;
res(row,j) := 0;
end loop;
end loop;
for i in 1..n loop -- fill rest with free elements
for j in 1..d loop
if res(i,j) = -1
then res(i,j) := 2;
end if;
end loop;
end loop;
return res;
end Localization_Map;
function Expand2 ( locmap : Matrix; b : Bracket ) return Poly is
-- DESCRIPTION :
-- Expands a 2-by-2 minor of the matrix, selecting the rows with
-- entries in the 2-bracket b, respecting the localization map.
res,subtmp : Poly;
n : constant natural := locmap'length(1);
d : constant natural := locmap'length(2);
b11 : Term := Create_Term(locmap(b(1),1),n,d,b(1),1);
b12 : Term := Create_Term(locmap(b(1),2),n,d,b(1),2);
b21 : Term := Create_Term(locmap(b(2),1),n,d,b(2),1);
b22 : Term := Create_Term(locmap(b(2),2),n,d,b(2),2);
d1 : Poly := Create(b11);
d2 : Poly := Create(b21);
begin -- res := b11*b22 - b21*b12;
res := b22*d1;
subtmp := b12*d2;
Sub(res,subtmp);
Clear(subtmp);
Clear(b22); Clear(b12);
Clear(b11); Clear(b21);
Clear(d1); Clear(d2);
return res;
end Expand2;
function Expand ( locmap : Matrix; b : Bracket ) return Poly is
-- DESCRIPTION :
-- Expands a d-by-d minor of the matrix, selecting the rows with
-- entries in b, respecting the localization map.
-- The expansion starts at the last column and proceeds recursively.
res : Poly;
n : constant natural := locmap'length(1);
d : constant natural := locmap'length(2);
begin
if b'last = 2
then res := Expand2(locmap,b);
else res := Null_Poly;
declare
sig : integer := +1;
b1 : Bracket(1..b'last-1);
begin
b1(1..b'last-1) := b(1..b'last-1);
if b'last mod 2 = 0
then sig := -1;
end if;
for i in b'range loop
declare
xt : Term := Create_Term(locmap(b(i),b'last),n,d,b(i),b'last);
expb1,xtexpb1 : Poly;
begin
if xt.cf /= Create(0.0)
then b1(1..i-1) := b(1..i-1);
b1(i..b1'last) := b(i+1..b'last);
expb1 := Expand(locmap,b1);
xtexpb1 := xt*expb1;
if sig > 0
then Add(res,xtexpb1);
else Sub(res,xtexpb1);
end if;
Clear(expb1); Clear(xtexpb1);
end if;
Clear(xt);
end;
sig := -sig;
end loop;
end;
end if;
return res;
end Expand;
function Reduce_Variables
( locmap : Matrix; dg : Standard_Natural_Vectors.Vector )
return Standard_Natural_Vectors.Vector is
-- DESCRIPTION :
-- Removes the #variables in the degrees vectors, removing all variables
-- that correspond to zeros in the localization map.
res : Standard_Natural_Vectors.Vector(dg'range) := dg;
cnt : natural := res'last;
d : constant natural := locmap'length(2);
begin
for i in reverse locmap'range(1) loop
for j in reverse locmap'range(2) loop
if locmap(i,j) /= 2
then cnt := cnt - 1;
for k in ((i-1)*d+j)..cnt loop
res(k) := res(k+1);
end loop;
end if;
end loop;
end loop;
return res(res'first..cnt);
end Reduce_Variables;
function Reduce_Variables ( locmap : Matrix; t : Term ) return Term is
-- DESCRIPTION :
-- Reduces the #variables in the term, removing all variables that
-- correspond to zeros in the localization map.
res : Term;
begin
res.cf := t.cf;
res.dg := new Standard_Natural_Vectors.Vector'
(Reduce_Variables(locmap,t.dg.all));
return res;
end Reduce_Variables;
procedure Reduce_Variables ( locmap : in Matrix; p : in out Poly ) is
rp : Poly := Null_Poly;
procedure Reduce_Term ( t : in Term; continue : out boolean ) is
rt : Term := Reduce_Variables(locmap,t);
begin
Add(rp,rt);
continue := true;
end Reduce_Term;
procedure Reduce_Terms is new Visiting_Iterator(Reduce_Term);
begin
Reduce_Terms(p);
Clear(p);
p := rp;
end Reduce_Variables;
end Bracket_Expansions;