File: [local] / OpenXM_contrib / PHC / Ada / Root_Counts / Product / random_product_start_systems.adb (download)
Revision 1.1.1.1 (vendor branch), Sun Oct 29 17:45:29 2000 UTC (23 years, 10 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_Random_Numbers; use Standard_Random_Numbers;
with Standard_Natural_Vectors;
with Standard_Complex_Vectors; use Standard_Complex_Vectors;
with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
with Random_Product_System;
with Set_Structure;
package body Random_Product_Start_Systems is
type Boolean_Array is array(integer range <>) of boolean;
function leq ( d1,d2 : Degrees ) return boolean is
-- DESCRIPTION :
-- returns true if for all k: d1(k) <= d2(k)
begin
for k in d1'range loop
if d1(k) > d2(k)
then return false;
end if;
end loop;
return true;
end leq;
function Dominated ( d : Degrees; p : Poly ) return boolean is
-- DESCRIPTION :
-- Returns true if the term indicated by the degrees d
-- is already in the set structure;
-- the polynomial p contains all terms that are already
-- belonging to the set structure.
res : boolean := false;
procedure Scan ( t : in Term; continue : out boolean ) is
begin
if leq(d,t.dg)
then res := true;
continue := false;
else continue := true;
end if;
end Scan;
procedure Scan_Terms is new Visiting_Iterator (Scan);
begin
Scan_Terms(p);
return res;
end Dominated;
procedure Divide ( dg : in out Standard_Natural_Vectors.Vector;
i,d : in natural; occupied : in out Boolean_Array ) is
-- DESCRIPTION :
-- Divides the monomial by every unknown that occurs already
-- in one set of the ith component of the set structure.
-- Every set that contains already one unknown of the monomial is
-- marked by setting the corresponding entry in occupied on true.
j : natural;
begin
for k in dg'range loop
j := 1;
while (dg(k) /= 0) and (j <= d) loop
if not occupied(j) and Set_Structure.Is_In(i,j,k)
then occupied(j) := true;
dg(k) := dg(k)-1;
end if;
j := j+1;
end loop;
end loop;
end Divide;
procedure Augment ( ba : in Boolean_Array; index : in out integer ) is
-- DESCRIPTION :
-- Augments the index till ba(index) = false.
begin
while ba(index) loop
index := index + 1;
exit when index >= ba'last;
end loop;
end Augment;
procedure Build_Set_Structure ( i,di : in natural; p : in Poly ) is
-- DESCRIPTION :
-- This procedure builds the set structure for the polynomial p
-- ON ENTRY :
-- i the number of the equation
-- di the degree of the polynomial
-- p the polynomial for the i-th equation
temp : Poly := Null_Poly;
procedure Process ( t : in Term; continue : out boolean ) is
ind : natural := 1; -- number of the set
processed : Boolean_Array(1..di) := (1..di => false);
tdg : Standard_Natural_Vectors.Vector(t.dg'range) := t.dg.all;
begin
if not Dominated(t.dg,temp)
then Divide(tdg,i,di,processed);
Augment(processed,ind);
for k in tdg'range loop
for j in 1..tdg(k) loop
if not Set_Structure.Is_In(i,ind,k)
then Set_Structure.Add(i,ind,k);
processed(ind) := true;
end if;
Augment(processed,ind);
end loop;
end loop;
Add(temp,t);
end if;
continue := true;
end Process;
procedure Process_Terms is new Visiting_Iterator (Process);
begin
Process_Terms(p);
Clear(temp);
end Build_Set_Structure;
procedure Build_Set_Structure ( p : in Poly_Sys ) is
n : natural := p'length;
ns : Standard_Natural_Vectors.Vector(1..n);
begin
for i in p'range loop
ns(i) := Degree(p(i));
end loop;
Set_Structure.Init(ns);
for i in p'range loop
Build_Set_Structure(i,Degree(p(i)),p(i));
end loop;
end Build_Set_Structure;
procedure Build_Random_Product_System ( n : in natural ) is
-- DESCRIPTION :
-- Based upon the set structure of p, a random product system
-- will be constructed.
h : Vector(0..n);
first : boolean;
begin
for i in 1..n loop
for j in 1..Set_Structure.Number_Of_Sets(i) loop
h := (0..n => Create(0.0));
first := true;
for k in 1..n loop
if Set_Structure.Is_In(i,j,k)
then if first
then h(k) := Create(1.0); first := false;
else h(k) := Random1;
end if;
end if;
end loop;
h(0) := Random1;
Random_Product_System.Add_Hyperplane(i,h);
end loop;
end loop;
end Build_Random_Product_System;
procedure Construct ( p : in Poly_Sys; q : in out Poly_Sys;
sols : in out Solution_List ) is
nl : natural := 0;
n : natural := p'length;
begin
Build_Set_Structure(p);
Random_Product_System.Init(n);
Build_Random_Product_System(n);
q := Random_Product_System.Polynomial_System;
Random_Product_System.Solve(sols,nl);
Random_Product_System.Clear;
end Construct;
end Random_Product_Start_Systems;