File: [local] / OpenXM_contrib / PHC / Ada / Schubert / curves_into_grassmannian.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_Natural_Vectors;
with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
package body Curves_into_Grassmannian is
-- CREATOR :
function Symbolic_Create ( m,p,q : natural; top,bottom : Bracket )
return Standard_Complex_Poly_Matrices.Matrix is
res : Standard_Complex_Poly_Matrices.Matrix(1..(m+p),1..p);
rws : constant natural := (m+p)*(q+1);
n : constant natural := Number_of_Variables(top,bottom) + 2;
row,ind,s_deg,t_deg : natural;
t : Term;
begin
for i in res'range(1) loop -- initialization
for j in res'range(2) loop
res(i,j) := Null_Poly;
end loop;
end loop;
t.dg := new Standard_Natural_Vectors.Vector'(1..n => 0);
t.cf := Create(1.0);
ind := 0; -- ind counts #variables
for j in 1..p loop -- assign columnwise
t_deg := (bottom(j)-1)/(m+p); -- degree in t to homogenize
row := 0; s_deg := 0;
for i in 1..rws loop
row := row + 1;
if i >= top(j) and i <= bottom(j)
then ind := ind+1;
t.dg(n-1) := s_deg; t.dg(n) := t_deg;
t.dg(ind) := 1;
Add(res(row,j),t);
t.dg(n-1) := 0; t.dg(n) := 0;
t.dg(ind) := 0;
end if;
if i mod (m+p) = 0
then row := 0; s_deg := s_deg+1;
if t_deg > 0
then t_deg := t_deg-1;
end if;
end if;
end loop;
end loop;
Clear(t);
return res;
end Symbolic_Create;
-- SELECTORS :
function Number_of_Variables ( top,bottom : Bracket ) return natural is
cnt : natural := 0;
begin
for j in top'range loop
cnt := cnt + (bottom(j) - top(j) + 1);
end loop;
return cnt;
end Number_of_Variables;
function Standard_Coordinate_Frame
( m,p,q : natural; top,bottom : Bracket;
coeff : Standard_Complex_Matrices.Matrix )
return Standard_Natural_Matrices.Matrix is
rws : constant natural := (m+p)*(q+1);
res : Standard_Natural_Matrices.Matrix(1..rws,1..p);
tol : constant double_float := 10.0**(-10);
first : boolean;
begin
for j in 1..p loop
first := true;
for i in 1..rws loop
if i < top(j) or i > bottom(j)
then res(i,j) := 0;
elsif (first and (AbsVal(coeff(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 Eval ( c : Term; s,t : Complex_Number ) return Term is
res : Term;
begin
Copy(c,res);
for i in 1..res.dg(res.dg'last-1) loop -- evaluate s
res.cf := res.cf*s;
end loop;
res.dg(res.dg'last-1) := 0;
for i in 1..res.dg(res.dg'last) loop -- evaluate t
res.cf := res.cf*t;
end loop;
res.dg(res.dg'last) := 0;
return res;
end Eval;
function Eval ( p : Poly; s,t : Complex_Number ) return Poly is
res : Poly := Null_Poly;
procedure Eval_Term ( ct : in Term; continue : out boolean ) is
et : Term := Eval(ct,s,t);
begin
Add(res,et);
Clear(et);
continue := true;
end Eval_Term;
procedure Eval_Terms is new Visiting_Iterator(Eval_Term);
begin
Eval_Terms(p);
return res;
end Eval;
function Eval ( c : Standard_Complex_Poly_Matrices.Matrix;
s,t : Complex_Number )
return Standard_Complex_Poly_Matrices.Matrix is
res : Standard_Complex_Poly_Matrices.Matrix(c'range(1),c'range(2));
begin
for i in c'range(1) loop
for j in c'range(2) loop
if c(i,j) = Null_Poly
then res(i,j) := Null_Poly;
else res(i,j) := Eval(c(i,j),s,t);
end if;
end loop;
end loop;
return res;
end Eval;
function Elim ( c : Term; s,t : Complex_Number ) return Term is
res : Term;
begin
res.dg := new Standard_Natural_Vectors.Vector'
(c.dg(c.dg'first..c.dg'last-2));
res.cf := c.cf;
for i in 1..c.dg(c.dg'last-1) loop -- evaluate s
res.cf := res.cf*s;
end loop;
for i in 1..c.dg(c.dg'last) loop -- evaluate t
res.cf := res.cf*t;
end loop;
return res;
end Elim;
function Elim ( p : Poly; s,t : Complex_Number ) return Poly is
res : Poly := Null_Poly;
procedure Elim_Term ( ct : in Term; continue : out boolean ) is
et : Term := Elim(ct,s,t);
begin
Add(res,et);
Clear(et);
continue := true;
end Elim_Term;
procedure Elim_Terms is new Visiting_Iterator(Elim_Term);
begin
Elim_Terms(p);
return res;
end Elim;
function Elim ( c : Standard_Complex_Poly_matrices.Matrix;
s,t : Complex_Number )
return Standard_Complex_Poly_Matrices.Matrix is
res : Standard_Complex_Poly_Matrices.Matrix(c'range(1),c'range(2));
begin
for i in c'range(1) loop
for j in c'range(2) loop
if c(i,j) = Null_Poly
then res(i,j) := Null_Poly;
else res(i,j) := Elim(c(i,j),s,t);
end if;
end loop;
end loop;
return res;
end Elim;
function Column_Localize ( top,bottom : Bracket;
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.
-- The localization pattern is traversed columnwise.
res : Term;
ndg : Standard_Natural_Vectors.Vector(t.dg'range);
cnt : natural := t.dg'first-1;
ind : natural := cnt;
begin
for j in locmap'range(2) loop -- columnwise order of the variables
for i in top(j)..bottom(j) loop -- restricted range skips the zeros
ind := ind+1;
if locmap(i,j) = 2 -- skip the ones
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 Column_Localize;
function Column_Localize ( top,bottom : Bracket;
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 Column_Localize_Term ( t : in Term; continue : out boolean ) is
lt : Term := Column_Localize(top,bottom,locmap,t);
begin
Add(res,lt);
Clear(lt.dg);
continue := true;
end Column_Localize_Term;
procedure Column_Localize_Terms is
new Visiting_Iterator(Column_Localize_Term);
begin
Column_Localize_Terms(p);
return res;
end Column_Localize;
function Column_Localize ( top,bottom : Bracket;
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) := Column_Localize(top,bottom,locmap,p(i));
end loop;
return res;
end Column_Localize;
function Column_Vector_Rep
( top,bottom : Bracket;
cffmat : Standard_Complex_Matrices.Matrix )
return Standard_Complex_Vectors.Vector is
dim : constant natural := Number_of_Variables(top,bottom);
res : Standard_Complex_Vectors.Vector(1..dim);
cnt : natural := 0;
begin
for j in cffmat'range(2) loop
for i in top(j)..bottom(j) loop
cnt := cnt + 1;
res(cnt) := cffmat(i,j);
end loop;
end loop;
return res;
end Column_Vector_Rep;
function Column_Vector_Rep ( locmap : Standard_Natural_Matrices.Matrix;
cffmat : Standard_Complex_Matrices.Matrix )
return Standard_Complex_Vectors.Vector is
dim : constant natural := cffmat'length(1)*cffmat'length(2);
res : Standard_Complex_Vectors.Vector(1..dim);
cnt : natural := 0;
begin
for j in cffmat'range(2) loop
for i in cffmat'range(1) loop
if locmap(i,j) = 2
then cnt := cnt + 1;
res(cnt) := cffmat(i,j);
end if;
end loop;
end loop;
return res(1..cnt);
end Column_Vector_Rep;
function Column_Matrix_Rep
( locmap : Standard_Natural_Matrices.Matrix;
cffvec : Standard_Complex_Vectors.Vector )
return Standard_Complex_Matrices.Matrix is
res : Standard_Complex_Matrices.Matrix(locmap'range(1),locmap'range(2));
cnt : natural := 0;
begin
for j in locmap'range(2) loop
for i in locmap'range(1) loop
if locmap(i,j) = 0
then res(i,j) := Create(0.0);
elsif locmap(i,j) = 1
then res(i,j) := Create(1.0);
else cnt := cnt + 1;
res(i,j) := cffvec(cnt);
end if;
end loop;
end loop;
return res;
end Column_Matrix_Rep;
procedure Swap ( p : in out Poly; k,l : in natural ) is
procedure Swap_Term ( t : in out Term; continue : out boolean ) is
lval : natural := t.dg(l);
begin
t.dg(l) := t.dg(k);
t.dg(k) := lval;
continue := true;
end Swap_Term;
procedure Swap_Terms is new Changing_Iterator(Swap_Term);
begin
Swap_Terms(p);
end Swap;
procedure Swap ( c : in out Standard_Complex_Poly_Matrices.Matrix;
k,l : in natural ) is
begin
for i in c'range(1) loop
for j in c'range(2) loop
if c(i,j) /= Null_Poly
then Swap(c(i,j),k,l);
end if;
end loop;
end loop;
end Swap;
function Insert ( p : Poly; k : natural ) return Poly is
res : Poly := Null_Poly;
procedure Insert_Term ( t : in Term; continue : out boolean ) is
rt : Term;
begin
rt.cf := t.cf;
rt.dg := new Standard_Natural_Vectors.Vector(t.dg'first..t.dg'last+1);
rt.dg(t.dg'first..k-1) := t.dg(t.dg'first..k-1);
rt.dg(k) := 0;
rt.dg(k+1..rt.dg'last) := t.dg(k..t.dg'last);
Add(res,rt);
Clear(rt);
continue := true;
end Insert_Term;
procedure Insert_Terms is new Visiting_Iterator(Insert_Term);
begin
Insert_Terms(p);
return res;
end Insert;
function Insert ( c : Standard_Complex_Poly_Matrices.Matrix; k : natural )
return Standard_Complex_Poly_Matrices.Matrix is
res : Standard_Complex_Poly_Matrices.Matrix(c'range(1),c'range(2));
begin
for i in c'range(1) loop
for j in c'range(2) loop
if c(i,j) = Null_Poly
then res(i,j) := Null_Poly;
else res(i,j) := Insert(c(i,j),k);
end if;
end loop;
end loop;
return res;
end Insert;
procedure Duplicate ( p : in out Poly; k,l : in natural ) is
procedure Duplicate_Term ( t : in out Term; continue : out boolean ) is
begin
t.dg(l) := t.dg(k);
continue := true;
end Duplicate_Term;
procedure Duplicate_Terms is new Changing_Iterator(Duplicate_Term);
begin
Duplicate_Terms(p);
end Duplicate;
end Curves_into_Grassmannian;