File: [local] / OpenXM_contrib / PHC / Ada / Schubert / determinantal_systems.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 Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
with Bracket_Monomials; use Bracket_Monomials;
with Bracket_Systems; use Bracket_Systems;
with Evaluated_Minors; use Evaluated_Minors;
with Symbolic_Minor_Equations; use Symbolic_Minor_Equations;
with Numeric_Minor_Equations; use Numeric_Minor_Equations;
package body Determinantal_Systems is
-- AUXILIARY TO EVALUATION OF MAXIMAL MINORS :
function Number_of_Maximal_Minors ( nrows,ncols : natural ) return natural is
-- DESCRIPTION :
-- Returns the number of maximal minors of a matrix with nrows and ncols.
-- REQUIRED : nrows > ncols.
begin
if ncols = 1
then return nrows;
else return Number_of_Maximal_Minors(nrows-1,ncols-1)*nrows/ncols;
end if;
end Number_of_Maximal_Minors;
-- LOCALIZATION MAPS :
function Localize ( 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.
res : Term;
ndg : Standard_Natural_Vectors.Vector(t.dg'range);
cnt : natural := t.dg'first-1;
ind : natural := cnt;
begin
for i in locmap'range(1) loop -- lexicographic order of variables
for j in locmap'range(2) loop
ind := ind+1;
if locmap(i,j) = 2
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 Localize;
function Localize ( 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 Localize_Term ( t : in Term; continue : out boolean ) is
lt : Term := Localize(locmap,t);
begin
Add(res,lt);
Clear(lt.dg);
continue := true;
end Localize_Term;
procedure Localize_Terms is new Visiting_Iterator(Localize_Term);
begin
Localize_Terms(p);
return res;
end Localize;
-- TARGET ROUTINES :
function Standard_Coordinate_Frame
( x : Standard_Complex_Poly_Matrices.Matrix;
plane : Standard_Complex_Matrices.Matrix )
return Standard_Natural_Matrices.Matrix is
res : Standard_Natural_Matrices.Matrix(x'range(1),x'range(2));
tol : constant double_float := 10.0**(-10);
first : boolean;
begin
for j in res'range(2) loop
first := true;
for i in res'range(1) loop
if x(i,j) = Null_Poly
then res(i,j) := 0;
elsif (first and (AbsVal(plane(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 Maximal_Coordinate_Frame
( x : Standard_Complex_Poly_Matrices.Matrix;
plane : Standard_Complex_Matrices.Matrix )
return Standard_Natural_Matrices.Matrix is
res : Standard_Natural_Matrices.Matrix(x'range(1),x'range(2));
max,tmp : double_float;
ind : natural;
begin
for j in res'range(2) loop
for i in res'range(1) loop
if x(i,j) = Null_Poly
then res(i,j) := 0;
else res(i,j) := 2;
end if;
end loop;
max := 0.0; ind := 0;
for i in res'range(1) loop
tmp := AbsVal(plane(i,j));
if tmp > max
then max := tmp; ind := i;
end if;
end loop;
res(ind,j) := 1;
end loop;
return res;
end Maximal_Coordinate_Frame;
function Localize ( 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) := Localize(locmap,p(i));
end loop;
return res;
end Localize;
-- CONSTRUCT THE POLYNOMIAL EQUATIONS :
function Polynomial_Equations
( l : Standard_Complex_Matrices.Matrix;
x : Standard_Complex_Poly_Matrices.Matrix ) return Poly_Sys is
n : constant natural := x'length(1);
p : constant natural := x'length(2);
kd : constant natural := p + l'length(2);
bm : Bracket_Monomial := Maximal_Minors(n,kd);
bs : Bracket_System(0..Number_of_Brackets(bm))
:= Minor_Equations(kd,kd-p,bm);
res : Poly_Sys(bs'first+1..bs'last) := Expanded_Minors(l,x,bs);
begin
Clear(bm); Clear(bs);
return res;
end Polynomial_Equations;
procedure Concat ( l : in out Link_to_Poly_Sys; p : Poly_Sys ) is
begin
if l = null
then declare
newsys : Poly_Sys(p'range);
cnt : natural := p'first-1;
begin
for i in p'range loop
if p(i) /= Null_Poly
then cnt := cnt+1;
newsys(cnt) := p(i);
end if;
end loop;
l := new Poly_Sys'(newsys(p'first..cnt));
end;
else declare
newsys : Poly_Sys(l'first..l'last+p'length);
ind : natural := l'last;
begin
newsys(l'range) := l(l'range);
for i in p'range loop
if p(i) /= Null_Poly
then ind := ind+1;
newsys(ind) := p(i);
end if;
end loop;
l := new Poly_Sys'(newsys(l'first..ind));
end;
end if;
end Concat;
function Polynomial_Equations
( l : Standard_Complex_VecMats.VecMat;
x : Standard_Complex_Poly_Matrices.Matrix ) return Poly_Sys is
n : constant natural := x'length(1);
p : constant natural := x'length(2);
kd : constant natural := p + l(l'first)'length(2);
bm : Bracket_Monomial := Maximal_Minors(n,kd);
bs : Bracket_System(0..Number_of_Brackets(bm))
:= Minor_Equations(kd,kd-p,bm);
nb : constant natural := l'length;
res : Poly_Sys(bs'first+1..nb*bs'last);
cnt : natural := bs'first;
begin
for i in 1..nb loop
declare
sys : constant Poly_Sys := Expanded_Minors(l(i).all,x,bs);
begin
for j in sys'range loop
cnt := cnt + 1;
res(cnt) := sys(j);
end loop;
end;
end loop;
Clear(bm); Clear(bs);
return res;
end Polynomial_Equations;
-- EVALUATORS AND DIFFERENTIATORS :
function Eval ( l,x : Standard_Complex_Matrices.Matrix )
return Complex_Number is
wrk : Standard_Complex_Matrices.Matrix(x'range(1),x'range(1));
begin
for j in l'range(2) loop
for i in l'range(1) loop
wrk(i,j) := l(i,j);
end loop;
end loop;
for j in x'range(2) loop
for i in x'range(1) loop
wrk(i,l'last(2)+j) := x(i,j);
end loop;
end loop;
return Determinant(wrk);
end Eval;
function Eval ( l,x : Standard_Complex_Matrices.Matrix ) return Vector is
n : constant natural := l'length(2) + x'length(2);
dim : constant natural := Number_of_Maximal_Minors(l'length(1),n);
res : Standard_Complex_Vectors.Vector(1..dim);
-- := (1..dim => Create(0.0));
cnt : natural := 0;
rws : Standard_Natural_Vectors.Vector(1..n);
wrk : Standard_Complex_Matrices.Matrix(1..n,1..n);
procedure Choose_next_Row ( k,start : in natural ) is
-- DESCRIPTION :
-- Chooses next k-th row within the range start..l'last(1), if k <= n.
-- If k > n, then the minor defined by the current row selection
-- is computed and added to the result.
begin
if k > n
then for j in l'range(2) loop -- select rows
for i in 1..n loop
wrk(i,j) := l(rws(i),j);
end loop;
end loop;
for j in x'range(2) loop
for i in 1..n loop
wrk(i,l'last(2)+j) := x(rws(i),j);
end loop;
end loop;
cnt := cnt + 1;
res(cnt) := Determinant(wrk);
else for i in start..l'last(1) loop
rws(k) := i;
Choose_next_Row(k+1,i+1);
end loop;
end if;
end Choose_next_Row;
begin
Choose_next_Row(1,1);
return res;
end Eval;
function Minor ( l,x : Standard_Complex_Matrices.Matrix; row,col : natural )
return Complex_Number is
-- DESCRIPTION :
-- Returns the minor [l|x] with row and column removed and with
-- the appropriate sign.
wrk : Standard_Complex_Matrices.Matrix
(x'first(1)..x'last(1)-1,x'first(1)..x'last(1)-1);
ind : natural;
begin
for j in l'range(2) loop
for i in l'range(1) loop
if i < row
then wrk(i,j) := l(i,j);
elsif i > row
then wrk(i-1,j) := l(i,j);
end if;
end loop;
end loop;
for j in x'range(2) loop
if j < col
then ind := l'last(2) + j;
elsif j > col
then ind := l'last(2) + j - 1;
end if;
if j /= col
then for i in x'range(1) loop
if i < row
then wrk(i,ind) := x(i,j);
elsif i > row
then wrk(i-1,ind) := x(i,j);
end if;
end loop;
end if;
end loop;
if (row + l'last(2) + col) mod 2 = 0
then return Determinant(wrk);
else return -Determinant(wrk);
end if;
end Minor;
function Diff ( l,x : Standard_Complex_Matrices.Matrix; i : natural )
return Complex_Number is
p : constant natural := x'length(2);
row,col : natural;
begin
row := i/p + 1;
col := i mod p;
if col = 0
then col := p;
row := row - 1;
end if;
return Minor(l,x,row,col);
end Diff;
function Diff ( l,x : Standard_Complex_Matrices.Matrix;
locmap : Standard_Natural_Matrices.Matrix; i : natural )
return Complex_Number is
row,col : natural;
cnt : natural := 0;
begin
for k1 in locmap'range(1) loop
for k2 in locmap'range(2) loop
if locmap(k1,k2) = 2
then cnt := cnt+1;
if cnt = i
then row := k1;
col := k2;
end if;
end if;
exit when (cnt = i);
end loop;
exit when (cnt = i);
end loop;
return Minor(l,x,row,col);
end Diff;
function Eval ( l : Standard_Complex_VecMats.VecMat;
x : Standard_Complex_Matrices.Matrix )
return Standard_Complex_Vectors.Vector is
res : Standard_Complex_Vectors.Vector(l'range);
begin
for i in l'range loop
res(i) := Eval(l(i).all,x);
end loop;
return res;
end Eval;
function Diff ( l : Standard_Complex_VecMats.VecMat;
x : Standard_Complex_Matrices.Matrix )
return Standard_Complex_Matrices.Matrix is
res : Standard_Complex_Matrices.Matrix(l'range,1..x'last(1)*x'last(2));
begin
for i in res'range(1) loop
for j in res'range(2) loop
res(i,j) := Diff(l(i).all,x,j);
end loop;
end loop;
return res;
end Diff;
function Old_Diff ( l : Standard_Complex_VecMats.VecMat;
x : Standard_Complex_Matrices.Matrix; nvars : natural;
locmap : Standard_Natural_Matrices.Matrix )
return Standard_Complex_Matrices.Matrix is
res : Standard_Complex_Matrices.Matrix(l'range,1..nvars);
begin
for i in res'range(1) loop
for j in res'range(2) loop
res(i,j) := Diff(l(i).all,x,locmap,j);
end loop;
end loop;
return res;
end Old_Diff;
function Diff ( l : Standard_Complex_VecMats.VecMat;
x : Standard_Complex_Matrices.Matrix; nvars : natural;
locmap : Standard_Natural_Matrices.Matrix )
return Standard_Complex_Matrices.Matrix is
-- NOTE :
-- This Diff is organized according to the localization map,
-- to avoid multiple searches.
res : Standard_Complex_Matrices.Matrix(l'range,1..nvars);
ind : natural := 0;
begin
for i in locmap'range loop
for j in locmap'range loop
if locmap(i,j) = 2
then ind := ind+1;
for k in l'range loop
res(k,ind) := Minor(l(k).all,x,i,j);
end loop;
end if;
end loop;
end loop;
return res;
end Diff;
end Determinantal_Systems;