File: [local] / OpenXM_contrib / PHC / Ada / Root_Counts / Implift / set_structures_and_volumes.adb (download)
Revision 1.1.1.1 (vendor branch), Sun Oct 29 17:45:29 2000 UTC (23 years, 7 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_Vectors; use Standard_Complex_Vectors;
with Standard_Complex_Matrices; use Standard_Complex_Matrices;
with Standard_Complex_Linear_Solvers; use Standard_Complex_Linear_Solvers;
with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
with Standard_Complex_Poly_Randomizers; use Standard_Complex_Poly_Randomizers;
with Standard_Complex_Substitutors; use Standard_Complex_Substitutors;
with Standard_Complex_Laur_Systems; use Standard_Complex_Laur_Systems;
with Standard_Poly_Laur_Convertors; use Standard_Poly_Laur_Convertors;
with Arrays_of_Integer_Vector_Lists; use Arrays_of_Integer_Vector_Lists;
with Set_Structure;
with Random_Product_System;
with Random_Product_Start_Systems;
with Power_Lists; use Power_Lists;
with Volumes; use Volumes;
with Mixed_Homotopy_Continuation; use Mixed_Homotopy_Continuation;
package body Set_Structures_and_Volumes is
-- AUXILIARIES :
procedure Build_RPS ( k,n : in natural; p : in Poly_Sys ) is
-- DESCRIPTION :
-- If the set structure is empty, then a set structure will
-- be constructed for the first k polynomials of p;
-- the data in Random_Product_System will be build.
begin
if Set_Structure.Empty
then declare
tmp : Poly_Sys(p'range);
t : Term;
cnst : Poly;
begin
tmp(1..k) := p(1..k);
t.dg := new Standard_Natural_Vectors.Vector'(1..n => 0);
t.cf := Create(1.0);
cnst := Create(t);
tmp(k+1..n) := (k+1..n => cnst);
Random_Product_Start_Systems.Build_Set_Structure(tmp);
Standard_Natural_Vectors.Clear
(Standard_Natural_Vectors.Link_to_Vector(t.dg));
Clear(cnst);
end;
end if;
Random_Product_System.Init(n);
Random_Product_Start_Systems.Build_Random_Product_System(n);
end Build_RPS;
procedure Build_Elimination_Matrix
( k,n : in natural; ind : in Standard_Natural_Vectors.Vector;
m : in out Matrix;
pivots : in out Standard_Natural_Vectors.Vector;
degenerate : out boolean ) is
-- DESCRIPTION :
-- builds the elimination matrix m.
-- ON ENTRY :
-- k the number of unknowns to be eliminated;
-- n the number of polynomials and unknowns;
-- ind entries in Random_Product_System;
-- ind(l) indicates lth hyperplane;
-- ON RETURN :
-- m contains all k hyperplanes to use for elimination;
-- pivots if not degenerate, then m(i,pivots(i)) /= 0;
-- degenerate is true if the first k hyperplanes were inconsistent.
degen : boolean;
begin
for i in ind'range loop -- build the matrix
declare
h : Vector(0..n) := Random_Product_System.Get_Hyperplane(i,ind(i));
begin
for j in 1..n loop
m(i,j) := h(j);
end loop;
m(i,n+1) := h(0);
end;
end loop;
diagonalize(m,k,n+1);
degen := false; -- check degeneracy
declare
eps : constant double_float := 10.0**(-10);
begin
for i in pivots'range loop
for j in 1..n loop
if AbsVal(m(i,j)) > eps
then pivots(i) := j;
end if;
exit when (pivots(i) /= 0);
end loop;
degen := (pivots(i) = 0);
exit when degen;
end loop;
end;
degenerate := degen;
end Build_Elimination_Matrix;
procedure Eliminate ( k,n : in natural; p : in Poly_Sys; m : in Matrix;
pivots : in Standard_Natural_Vectors.Vector;
q : in out Poly_Sys ) is
-- DESCRIPTION :
-- eliminates k unknowns of the polynomial system p.
-- ON ENTRY :
-- k the number of unknowns to be eliminated;
-- n the number of polynomials and unknowns;
-- p a polynomial system with random coefficients;
-- m the diagonalized matrix, to be used for elimination,
-- for i in m'range(1) : m(i,pivots(i)) /= 0;
-- pivots a vector indicating nonzero entries in m.
-- ON RETURN :
-- q the reduced system.
begin
Clear(q); Copy(p,q);
for i in pivots'range loop
declare
h : Vector(0..n-i+1);
piv : natural := pivots(i)-i+1;
begin
h(0) := m(i,n+1);
h(1..n-i+1) := (1..n-i+1 => Create(0.0));
h(piv) := m(i,pivots(i));
for j in piv+1..n-i+1 loop
h(j) := m(i,j+i-1);
end loop;
Substitute(piv,h,q);
end;
end loop;
end Eliminate;
procedure Update ( sols : in out Solution_List; k,n : in natural;
m : in Matrix;
pivots : in Standard_Natural_Vectors.Vector ) is
-- DESCRIPTION :
-- Based on the elimination matrix m, solution components for x1..xk
-- can be computed.
-- ON ENTRY :
-- sols a list of solutions, with values for x(n-k+1)..x(n);
-- k the number of unknowns that have been eliminated;
-- n total number of unknowns;
-- m the elimination matrix;
-- pivots for i in pivots'range : m(i,pivots(i)) /= 0.
-- ON RETURN :
-- sols the solution list with n-dimensional vectors.
tmp,res,res_last : Solution_List;
begin
tmp := sols;
while not Is_Null(tmp) loop
declare
sol : Solution(n-k) := Head_Of(tmp).all;
soln : Solution(n);
ind : natural;
begin
soln.t := Create(0.0);
soln.m := sol.m;
soln.v := (1..n => Create(0.0));
for i in 1..k loop
soln.v(pivots(i)) := -m(i,n+1);
end loop;
ind := 0;
for i in 1..n loop
if soln.v(i) = Create(0.0)
then ind := ind + 1;
soln.v(i) := sol.v(ind);
end if;
end loop;
for i in 1..k loop
ind := pivots(i);
for j in ind+1..n loop
soln.v(ind) := soln.v(ind) - m(i,j)*soln.v(j);
end loop;
soln.v(ind) := soln.v(ind) / m(i,ind);
end loop;
Append(res,res_last,soln);
end;
tmp := Tail_Of(tmp);
end loop;
Clear(sols); Copy(res,sols);
end Update;
procedure Build_Indices ( l,k,n : in natural; p : in Poly_Sys;
acc : in out Standard_Natural_Vectors.Vector;
bb : in out natural ) is
-- DESCRIPTION :
-- Builds the indices in the vector acc.
-- ON ENTRY :
-- l index in the array acc;
-- k number of unknowns to eliminate;
-- n number of polynomials and unknowns;
-- p a polynomial system;
-- acc(1..l-1) contains entries indicating hyperplanes
-- in Random_Product_System.
-- ON RETURN :
-- acc(1..l) contains entries indicating hyperplanes
-- in Random_Product_System;
-- bb the mixed Bezout BKK bound.
begin
if l > k
then declare
degenerate : boolean;
m : Matrix(1..k,1..n+1);
pivots : Standard_Natural_Vectors.Vector(1..k);
begin
pivots := (1..k => 0);
Build_Elimination_Matrix(k,n,acc,m,pivots,degenerate);
if not degenerate
then declare
q : Poly_Sys(p'range);
begin
Eliminate(k,n,p,m,pivots,q);
bb := bb + Bezout_BKK(0,n-k,q);
Clear(q);
end;
end if;
end;
else for j in 1..Random_Product_System.Number_of_Hyperplanes(l) loop
acc(l) := j;
Build_Indices(l+1,k,n,p,acc,bb);
end loop;
end if;
end Build_Indices;
procedure Build_Indices ( l,k,n : in natural; p : in Poly_Sys;
tv : in Tree_of_Vectors;
acc : in out Standard_Natural_Vectors.Vector;
bb : in out natural ) is
-- DESCRIPTION :
-- Builds the indices in the vector acc.
-- ON ENTRY :
-- l index in the array acc;
-- k number of unknowns to eliminate;
-- n number of polynomials and unknowns;
-- p a polynomial system;
-- tv the tree of degrees used to compute the mixed volume;
-- acc(1..l-1) contains entries indicating hyperplanes
-- in Random_Product_System.
-- ON RETURN :
-- acc(1..l) contains entries indicating hyperplanes
-- in Random_Product_System;
-- bb the mixed Bezout BKK bound.
begin
if l > k
then declare
degenerate : boolean;
m : Matrix(1..k,1..n+1);
pivots : Standard_Natural_Vectors.Vector(1..k);
begin
pivots := (1..k => 0);
Build_Elimination_Matrix(k,n,acc,m,pivots,degenerate);
if not degenerate
then declare
q : Poly_Sys(p'range);
begin
Eliminate(k,n,p,m,pivots,q);
bb := bb + Bezout_BKK(0,n-k,q,tv);
Clear(q);
end;
end if;
end;
else for j in 1..Random_Product_System.Number_of_Hyperplanes(l) loop
acc(l) := j;
Build_Indices(l+1,k,n,p,tv,acc,bb);
end loop;
end if;
end Build_Indices;
procedure Build_Indices2 ( l,k,n : in natural; p : in Poly_Sys;
tv : in out Tree_of_Vectors;
acc : in out Standard_Natural_Vectors.Vector;
bb : in out natural ) is
-- DESCRIPTION :
-- Builds the indices in the vector acc.
-- ON ENTRY :
-- l index in the array acc;
-- k number of unknowns to eliminate;
-- n number of polynomials and unknowns;
-- p a polynomial system;
-- acc(1..l-1) contains entries indicating hyperplanes
-- in Random_Product_System.
-- ON RETURN :
-- tv the tree of degrees used to compute the mixed volume;
-- acc(1..l) contains entries indicating hyperplanes
-- in Random_Product_System;
-- bb the mixed Bezout BKK bound.
begin
if l > k
then declare
degenerate : boolean;
m : Matrix(1..k,1..n+1);
pivots : Standard_Natural_Vectors.Vector(1..k);
begin
pivots := (1..k => 0);
Build_Elimination_Matrix(k,n,acc,m,pivots,degenerate);
if not degenerate
then declare
q : Poly_Sys(p'range);
qtv,tmp : Tree_of_Vectors;
mv : natural;
begin
Eliminate(k,n,p,m,pivots,q);
Bezout_BKK(0,n-k,q,qtv,mv);
bb := bb + mv;
tmp := qtv;
while not Is_Null(tmp) loop
declare
nd : node := Head_Of(tmp);
begin
if Is_In(tv,nd.d)
then Clear(nd);
else Construct(nd,tv);
end if;
end;
tmp := Tail_Of(tmp);
end loop;
Clear(q);
end;
end if;
end;
else for j in 1..Random_Product_System.Number_of_Hyperplanes(l) loop
acc(l) := j;
Build_Indices2(l+1,k,n,p,tv,acc,bb);
end loop;
end if;
end Build_Indices2;
procedure Build_Indices
( file : in file_type; p : in Poly_Sys; l,k,n : in natural;
acc : in out Standard_Natural_Vectors.Vector;
bb : in out natural;
sols,sols_last : in out Solution_List ) is
-- DESCRIPTION :
-- Builds the indices in the vector acc.
-- ON ENTRY :
-- file file to write intermediate results on;
-- p a polynomial system;
-- l index in the array acc;
-- k number of unknowns to eliminate;
-- n number of polynomials and unknowns;
-- acc(1..l-1) contains entries indicating hyperplanes
-- in Random_Product_System.
-- ON RETURN :
-- acc(1..l) contains entries indicating hyperplanes
-- in Random_Product_System;
-- bb the mixed Bezout BKK bound;
-- sols the solutions of pi;
-- sols_last points to the last element of the list sols.
begin
if l > k
then declare
degenerate : boolean;
m : Matrix(1..k,1..n+1);
pivots : Standard_Natural_Vectors.Vector(1..k);
begin
pivots := (1..k => 0);
Build_Elimination_Matrix(k,n,acc,m,pivots,degenerate);
if not degenerate
then declare
q : Poly_Sys(p'range);
begin
Eliminate(k,n,p,m,pivots,q);
declare
las : Laur_Sys(q'range);
qsols : Solution_List;
bkk : natural;
begin
las := Polynomial_to_Laurent_System(q);
Solve(file,las,bkk,qsols);
bb := bb + bkk;
Update(qsols,k,n,m,pivots);
Concat(sols,sols_last,qsols);
Clear(las); Shallow_Clear(qsols);
end;
Clear(q);
end;
end if;
end;
else for j in 1..Random_Product_System.Number_of_Hyperplanes(l) loop
acc(l) := j;
Build_Indices(file,p,l+1,k,n,acc,bb,sols,sols_last);
end loop;
end if;
end Build_Indices;
procedure Build_Indices
( file : in file_type; p : in Poly_Sys; l,k,n : in natural;
acc : in out Standard_Natural_Vectors.Vector;
tv : in Tree_of_Vectors; bb : in out natural;
sols,sols_last : in out Solution_List ) is
-- DESCRIPTION :
-- Builds the indices in the vector acc.
-- ON ENTRY :
-- file file to write intermediate results on;
-- p a polynomial system;
-- l index in the array acc;
-- k number of unknowns to eliminate;
-- n number of polynomials and unknowns;
-- acc(1..l-1) contains entries indicating hyperplanes
-- in Random_Product_System;
-- tv the tree containing useful directions;
-- ON RETURN :
-- acc(1..l) contains entries indicating hyperplanes
-- in Random_Product_System;
-- bb the mixed Bezout BKK bound;
-- sols the solutions of pi;
-- sols_last points to the last element of the list sols.
begin
if l > k
then declare
degenerate : boolean;
m : Matrix(1..k,1..n+1);
pivots : Standard_Natural_Vectors.Vector(1..k);
begin
pivots := (1..k => 0);
Build_Elimination_Matrix(k,n,acc,m,pivots,degenerate);
if not degenerate
then declare
q : Poly_Sys(p'range);
begin
Eliminate(k,n,p,m,pivots,q);
declare
las : Laur_Sys(q'range);
qsols : Solution_List;
bkk : natural;
begin
las := Polynomial_to_Laurent_System(q);
Solve(file,las,tv,bkk,qsols);
bb := bb + bkk;
Update(qsols,k,n,m,pivots);
Concat(sols,sols_last,qsols);
Clear(las); Shallow_Clear(qsols);
end;
Clear(q);
end;
end if;
end;
else for j in 1..Random_Product_System.Number_of_Hyperplanes(l) loop
acc(l) := j;
Build_Indices(file,p,l+1,k,n,acc,tv,bb,sols,sols_last);
end loop;
end if;
end Build_Indices;
-- TARGET ROUTINES :
function Bezout_BKK ( k,n : natural; p : Poly_Sys ) return natural is
res : natural;
begin
if k = 0
then declare
adl : Array_of_Lists(p'range) := Create(p);
begin
res := Mixed_Volume(n,adl);
Deep_Clear(adl);
end;
elsif k = n
then if Set_Structure.Empty
then Random_Product_Start_Systems.Build_Set_Structure(p);
res := Set_Structure.B;
Set_Structure.Clear;
else res := Set_Structure.B;
end if;
else Build_RPS(k,n,p);
declare
acc : Standard_Natural_Vectors.Vector(1..k);
q : Poly_Sys(1..n-k);
begin
acc := (1..k => 0); res := 0;
for i in q'range loop
q(i) := Complex_Randomize1(p(i+k));
end loop;
Build_Indices(1,k,n,q,acc,res);
Clear(q);
end;
Random_Product_System.Clear;
end if;
return res;
end Bezout_BKK;
function Bezout_BKK ( k,n : natural; p : Poly_Sys; tv : Tree_of_Vectors )
return natural is
res : natural;
begin
if k = 0
then declare
adl : Array_of_Lists(p'range) := Create(p);
begin
res := Mixed_Volume(n,adl,tv);
Deep_Clear(adl);
end;
elsif k = n
then if Set_Structure.Empty
then Random_Product_Start_Systems.Build_Set_Structure(p);
res := Set_Structure.B;
Set_Structure.Clear;
else res := Set_Structure.B;
end if;
else Build_RPS(k,n,p);
declare
acc : Standard_Natural_Vectors.Vector(1..k);
q : Poly_Sys(1..n-k);
begin
acc := (1..k => 0); res := 0;
for i in q'range loop
q(i) := Complex_Randomize1(p(i+k));
end loop;
Build_Indices(1,k,n,q,tv,acc,res);
Clear(q);
end;
Random_Product_System.Clear;
end if;
return res;
end Bezout_BKK;
procedure Bezout_BKK ( k,n : in natural; p : in Poly_Sys;
tv : in out Tree_of_Vectors; bb : out natural ) is
res : natural;
begin
if k = 0
then declare
adl : Array_of_Lists(p'range) := Create(p);
begin
Mixed_Volume(n,adl,tv,res);
Deep_Clear(adl);
end;
elsif k = n
then if Set_Structure.Empty
then Random_Product_Start_Systems.Build_Set_Structure(p);
res := Set_Structure.B;
Set_Structure.Clear;
else res := Set_Structure.B;
end if;
else Build_RPS(k,n,p);
declare
acc : Standard_Natural_Vectors.Vector(1..k);
q : Poly_Sys(1..n-k);
begin
acc := (1..k => 0); res := 0;
for i in q'range loop
q(i) := Complex_Randomize1(p(i+k));
end loop;
Build_Indices2(1,k,n,q,tv,acc,res);
Clear(q);
end;
Random_Product_System.Clear;
end if;
bb := res;
end Bezout_BKK;
procedure Mixed_Solve
( file : in file_type; k,n : in natural; p : in Poly_Sys;
bb : out natural;
g : in out Poly_Sys; sols : in out Solution_List ) is
begin
if k = 0
then declare
l : Laur_Sys(g'range);
tmp : Solution_List;
begin
l := Polynomial_to_Laurent_System(g);
Solve(file,l,bb,sols);
tmp := sols;
while not Is_Null(tmp) loop
Head_Of(tmp).t := Create(0.0);
tmp := Tail_Of(tmp);
end loop;
Clear(l);
end;
elsif k = n
then Random_Product_Start_Systems.Construct(p,g,sols);
bb := Length_Of(sols);
else Build_RPS(k,n,p);
declare
acc : Standard_Natural_Vectors.Vector(1..k);
q : Poly_Sys(1..n-k);
rq : Poly_Sys(p'range);
sols_last : Solution_List;
res : natural;
begin
acc := (1..k => 0); res := 0;
for i in q'range loop
q(i) := g(i+k);
end loop;
Build_Indices(file,q,1,k,n,acc,res,sols,sols_last);
bb := res;
rq := Random_Product_System.Polynomial_System;
g(1..k) := rq(1..k);
end;
Random_Product_System.Clear;
end if;
end Mixed_Solve;
procedure Mixed_Solve
( file : in file_type; k,n : in natural; p : in Poly_Sys;
tv : in Tree_of_Vectors; bb : out natural;
g : in out Poly_Sys; sols : in out Solution_List ) is
begin
if k = 0
then declare
l : Laur_Sys(g'range);
tmp : Solution_List;
begin
l := Polynomial_to_Laurent_System(g);
Solve(file,l,tv,bb,sols);
tmp := sols;
while not Is_Null(tmp) loop
Head_Of(tmp).t := Create(0.0);
tmp := Tail_Of(tmp);
end loop;
Clear(l);
end;
elsif k = n
then Random_Product_Start_Systems.Construct(p,g,sols);
bb := Length_Of(sols);
else Build_RPS(k,n,p);
declare
acc : Standard_Natural_Vectors.Vector(1..k);
q : Poly_Sys(1..n-k);
rq : Poly_Sys(p'range);
sols_last : Solution_List;
res : natural;
begin
acc := (1..k => 0); res := 0;
for i in q'range loop
q(i) := g(i+k);
end loop;
Build_Indices(file,q,1,k,n,acc,tv,res,sols,sols_last);
bb := res;
rq := Random_Product_System.Polynomial_System;
g(1..k) := rq(1..k);
end;
Random_Product_System.Clear;
end if;
end Mixed_Solve;
end Set_Structures_and_Volumes;
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to set_structures_and_volumes.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Root_Counts / Implift