File: [local] / OpenXM_contrib / PHC / Ada / Root_Counts / Implift / mixed_homotopy_continuation.adb (download)
Revision 1.1.1.1 (vendor branch), Sun Oct 29 17:45:28 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_Floating_Numbers; use Standard_Floating_Numbers;
with Standard_Complex_Numbers; use Standard_Complex_Numbers;
with Standard_Complex_Numbers_Polar; use Standard_Complex_Numbers_Polar;
with Standard_Random_Numbers; use Standard_Random_Numbers;
with Standard_Integer_Vectors;
with Standard_Integer_Vectors_io; use Standard_Integer_Vectors_io;
with Standard_Complex_Vectors;
with Standard_Complex_Norms_Equals; use Standard_Complex_Norms_Equals;
with Standard_Complex_Matrices; use Standard_Complex_Matrices;
with Standard_Complex_Polynomials;
with Standard_Complex_Laur_Polys; use Standard_Complex_Laur_Polys;
with Standard_Complex_Poly_Systems; use Standard_Complex_Poly_Systems;
with Standard_Complex_Poly_Functions; use Standard_Complex_Poly_Functions;
with Standard_Complex_Poly_SysFun; use Standard_Complex_Poly_SysFun;
with Standard_Complex_Jaco_Matrices; use Standard_Complex_Jaco_Matrices;
with Standard_Poly_Laur_Convertors; use Standard_Poly_Laur_Convertors;
with Standard_Laur_Poly_Convertors; use Standard_Laur_Poly_Convertors;
with Lists_of_Integer_Vectors; use Lists_of_Integer_Vectors;
with Integer_Support_Functions; use Integer_Support_Functions;
with Arrays_of_Integer_Vector_Lists; use Arrays_of_Integer_Vector_Lists;
with Increment_and_Fix_Continuation; use Increment_and_Fix_Continuation;
with Standard_Root_Refiners; use Standard_Root_Refiners;
with Integer_Vectors_Utilities; use Integer_Vectors_Utilities;
with Lists_of_Vectors_Utilities; use Lists_of_Vectors_Utilities;
with Arrays_of_Lists_Utilities; use Arrays_of_Lists_Utilities;
with Transformations; use Transformations;
with Transforming_Solutions; use Transforming_Solutions;
with Transforming_Laurent_Systems; use Transforming_Laurent_Systems;
with Transforming_Integer_Vector_Lists; use Transforming_Integer_Vector_Lists;
with Power_Lists; use Power_Lists;
with Volumes;
with Durand_Kerner;
package body Mixed_Homotopy_Continuation is
-- INVARIANT CONDITION :
-- The procedures and functions in this package `mirror' corresponding
-- routines in the package Volumes.
-- AUXILIARIES :
type bar is array ( integer range <> ) of boolean;
function Interchange2 ( p : Laur_Sys; index : integer ) return Laur_Sys is
-- DESCRIPTION :
-- Returns a polynomial system where the first equation is interchanged
-- with the equation given by the index.
res : Laur_Sys(p'range);
begin
if index = p'first
then res := p;
else res(res'first) := p(index);
res(index) := p(p'first);
res(res'first+1..index-1) := p(p'first+1..index-1);
res(index+1..res'last) := p(index+1..p'last);
end if;
return res;
end Interchange2;
procedure Interchange2 ( adl : in out Array_of_Lists; index : integer ) is
-- DESCRIPTION :
-- Interchanges the first list with the list given by the index.
tmp : List;
begin
if index /= adl'first
then tmp := adl(adl'first);
adl(adl'first) := adl(index);
adl(index) := tmp;
end if;
end Interchange2;
function Permute ( perm : Standard_Integer_Vectors.Vector; p : Laur_Sys )
return laur_Sys is
-- DESCRIPTION :
-- Returns a permuted polynomial system.
res : Laur_Sys(p'range);
begin
for i in p'range loop
res(i) := p(perm(i));
end loop;
return res;
end Permute;
function Initial_Degrees ( p : Poly ) return Degrees is
-- DESCRIPTION :
-- Returns the initial degrees of the polynomial p.
init : Degrees;
procedure Init_Term ( t : in Term; cont : out boolean ) is
begin
init := new Standard_Integer_Vectors.Vector'(t.dg.all);
cont := false;
end Init_Term;
procedure Initial_Term is new Visiting_Iterator (Init_Term);
begin
Initial_Term(p);
return init;
end Initial_Degrees;
procedure Binomial ( p : in Poly;
d : out Standard_Integer_Vectors.Link_to_Vector;
k : in out integer; c : in out Complex_Number ) is
-- DESCRIPTION :
-- p consists of two terms, this procedure gets the degrees d and
-- the constant c of the binomial equation.
first : boolean := true;
dd : Degrees;
procedure Scan_Term ( t : in Term; cont : out boolean ) is
begin
if first
then dd := new Standard_Integer_Vectors.Vector'(t.dg.all);
c := t.cf;
first := false;
else k := dd'first - 1;
for i in dd'range loop
dd(i) := dd(i) - t.dg(i);
if k < dd'first and then dd(i) /= 0
then k := i;
end if;
end loop;
c := -t.cf/c;
end if;
cont := true;
end Scan_Term;
procedure Scan_Terms is new Visiting_Iterator (Scan_Term);
begin
Scan_Terms(p);
d := new Standard_Integer_Vectors.Vector'(dd.all);
Standard_Integer_Vectors.Clear(Standard_Integer_Vectors.Link_to_Vector(dd));
end Binomial;
procedure Normalize ( p : in Laur_Sys; dl : in out List;
wp : in out Laur_Sys; shifted : out boolean ) is
-- DESCRIPTION :
-- Makes sure that the first element of dl contains all zeroes.
-- REQUIRED :
-- The list dl is not empty.
-- ON ENTRY :
-- p a Laurent polynomial system;
-- dl the power list of p(p'first).
-- ON RETURN :
-- dl the first element of contains all zeroes,
-- therefore dl has been shifted;
-- wp a Laurent polynomial system where
-- dl is the power list of wp(wp'first);
-- shifted is true if dl has been shifted.
use Standard_Integer_Vectors;
first : Link_to_Vector := Head_Of(dl);
nullvec : Vector(first'range) := (first'range => 0);
shiftvec : Link_to_Vector;
begin
if not Is_In(dl,nullvec)
then shiftvec := Graded_Max(dl);
Shift(dl,shiftvec);
Clear(shiftvec);
Copy(p,wp);
for i in p'range loop
Shift(wp(i));
end loop;
else wp := p;
shifted := false;
end if;
Move_to_Front(dl,nullvec);
end Normalize;
function Evaluate ( p : Poly; x : Complex_Number; k : integer )
return Poly is
-- DESCRIPTION :
-- This function returns a polynomial where the kth unknown has
-- been replaced by x.
-- It is important to use this function as the term order in p
-- must remain the same!
res : Poly;
procedure Evaluate_Term ( t : in out Term; cont : out boolean ) is
fac : Complex_Number;
pow : natural;
begin
if t.dg(k) < 0
then fac := Create(1.0)/x;
pow := -t.dg(k);
else fac := x;
pow := t.dg(k);
end if;
for i in 1..pow loop
t.cf := t.cf * fac;
end loop;
declare
tmp : constant Standard_Integer_Vectors.Vector := t.dg.all;
begin
Clear(t);
t.dg := new Standard_Integer_Vectors.Vector'(Reduce(tmp,k));
end;
cont := true;
end Evaluate_Term;
procedure Evaluate_Terms is new Changing_Iterator (Evaluate_Term);
begin
Copy(p,res);
Evaluate_Terms(res);
return res;
end Evaluate;
function Re_Arrange ( p : poly ) return Poly is
-- DESCRIPTION :
-- Returns a polynomial whose terms are sorted
-- in graded lexicographical ordening.
res : Poly := Null_Poly;
procedure Re_Arrange_Term ( t : in Term; cont : out boolean ) is
begin
Add(res,t);
cont := true;
end Re_Arrange_Term;
procedure Re_Arrange_Terms is new Visiting_Iterator (Re_Arrange_Term);
begin
Re_Arrange_Terms(p);
return res;
end Re_Arrange;
function Substitute ( p : Poly; v : Standard_Complex_Vectors.Vector;
k : integer ) return Poly is
-- DESCRIPTION :
-- Substitutes the values in v into the polynomial p,
-- starting at the last unknowns of p.
init : Degrees := Initial_Degrees(p);
index : integer := init'last;
res,tmp : Poly;
begin
if index = k
then index := index - 1;
res := Evaluate(p,v(v'last),index);
else res := Evaluate(p,v(v'last),index);
end if;
for i in reverse v'first..(v'last-1) loop
index := index - 1;
if index = k
then index := index - 1;
end if;
tmp := Evaluate(res,v(i),index);
Clear(res); Copy(tmp,res); Clear(tmp);
end loop;
Standard_Integer_Vectors.Clear
(Standard_Integer_Vectors.Link_to_Vector(init));
return res;
end Substitute;
procedure Refine_and_Concat
( p : in Laur_Sys;
newsols,sols,sols_last : in out Solution_List ) is
-- DESCRIPTION :
-- This procedure refines the solutions of a given
-- polynomial system and adds them to the solution list.
-- This can be very useful for eliminating rounding errors
-- after transformating the solutions.
pp : Poly_Sys(p'range) := Laurent_to_Polynomial_System(p);
numb : natural := 0;
begin
Silent_Root_Refiner(pp,newsols,10.0**(-12),10.0**(-12),10.0**(-8),numb,5);
Concat(sols,sols_last,newsols);
Clear(pp); Shallow_Clear(newsols);
end Refine_and_Concat;
procedure Refine_and_Concat
( file : in file_type; p : in Laur_Sys;
newsols,sols,sols_last : in out Solution_List ) is
-- DESCRIPTION :
-- This procedure refines the solutions of a given
-- polynomial system and adds them to the solution list.
-- This can be very useful for eliminating rounding errors
-- after transformating the solutions.
pp : Poly_Sys(p'range) := Laurent_to_Polynomial_System(p);
numb : natural := 0;
begin
Reporting_Root_Refiner
(file,pp,newsols,10.0**(-12),10.0**(-12),10.0**(-8),numb,5,false);
Concat(sols,sols_last,newsols);
Clear(pp); Shallow_Clear(newsols);
end Refine_and_Concat;
procedure Write_Direction
( file : in file_type;
v : in Standard_Integer_Vectors.Link_to_Vector ) is
begin
put(file,"++++ considering direction "); put(file,v);
put_line(file," ++++");
end Write_Direction;
-- INTERMEDIATE LAYER :
procedure Mixed_Continue
( file : in file_type; p : in Laur_Sys;
k : in integer; m : in Standard_Integer_Vectors.Vector;
sols : in out Solution_List ) is
-- DESCRIPTION :
-- This continuation routine computes a part of the solution list
-- of a Laurent polynomial system.
-- ON ENTRY :
-- file a file where the intermediate results are written;
-- p the transformed Laurent polynomial system to be solved;
-- k the index;
-- m m(k) = p1(v), m(k/=l) = Maximal_Degree(p(l),k);
-- sols the start solutions.
-- ON RETURN :
-- sols the computed solutions.
h : Laur_Sys(p'range);
hp : Poly_Sys(h'range);
hpe : Eval_Poly_Sys(hp'range);
j : Jaco_Mat(p'range,p'first..p'last+1);
je : Eval_Jaco_Mat(j'range(1),j'range(2));
function Construct_Homotopy
( p : Laur_Sys; k : integer;
m : Standard_Integer_Vectors.Vector ) return Laur_Sys is
res : Laur_Sys(p'range);
ran : Complex_Number;
re : boolean;
zeroes : Degrees := new Standard_Integer_Vectors.Vector'(p'range => 0);
function Construct_First_Polynomial
( pp : Poly; max : integer ) return Poly is
r : Poly := Null_Poly;
procedure Construct_Term ( t : in Term; cont : out boolean ) is
rt : Term;
begin
rt.cf := t.cf;
rt.dg := new Standard_Integer_Vectors.Vector(t.dg'first..t.dg'last+1);
rt.dg(t.dg'range) := t.dg.all;
rt.dg(k) := -t.dg(k) + max;
if Equal(t.dg,zeroes)
then rt.dg(rt.dg'last) := 0;
re := ( IMAG_PART(rt.cf) + 1.0 = 1.0 );
else rt.dg(rt.dg'last) := rt.dg(k);
end if;
Add(r,rt);
Clear(rt);
cont := true;
end Construct_Term;
procedure Construct_Terms is new Visiting_Iterator(Construct_Term);
begin
Construct_Terms(pp);
Standard_Integer_Vectors.Clear
(Standard_Integer_Vectors.Link_to_Vector(zeroes));
return r;
end Construct_First_Polynomial;
function Construct_Polynomial ( pp : Poly; max : integer ) return Poly is
r : Poly := Null_Poly;
procedure Construct_Term ( t : in Term; cont : out boolean ) is
rt : Term;
begin
rt.cf := t.cf;
rt.dg := new Standard_Integer_Vectors.Vector(t.dg'first..t.dg'last+1);
rt.dg(t.dg'range) := t.dg.all;
rt.dg(k) := -t.dg(k) + max;
rt.dg(rt.dg'last) := rt.dg(k);
Add(r,rt);
Clear(rt);
cont := true;
end Construct_Term;
procedure Construct_Terms is new Visiting_Iterator(Construct_Term);
begin
Construct_Terms(pp);
return r;
end Construct_Polynomial;
begin
res(res'first) := Construct_First_Polynomial(p(p'first),m(m'first));
for i in p'first+1..p'last loop
res(i) := Construct_Polynomial(p(i),m(i));
end loop;
if re
then for i in res'range loop
ran := Random;
Mul(res(i),ran);
end loop;
end if;
return res;
end Construct_Homotopy;
function To_Be_Removed ( flag : in natural ) return boolean is
begin
return ( flag /= 1 );
end To_Be_Removed;
procedure Extract_Regular_Solutions is
new Standard_Complex_Solutions.Delete(To_Be_Removed);
begin
h := Construct_Homotopy(p,k,m); -- CONSTRUCTION OF HOMOTOPY
hp := Laurent_to_Polynomial_System(h);
hpe := Create(hp);
j := Create(hp);
je := Create(j);
declare -- CONTINUATION
use Standard_Complex_Vectors;
function Eval ( x : Vector; t : Complex_Number ) return Vector is
xt : Vector(x'first..x'last+1);
begin
xt(x'range) := x;
xt(xt'last) := t;
return Eval(hpe,xt);
end Eval;
function dHt ( x : Vector; t : Complex_Number ) return Vector is
xt : Vector(x'first..x'last+1);
res : Vector(p'range);
begin
xt(x'range) := x;
xt(xt'last) := t;
for i in res'range loop
res(i) := Eval(je(i,xt'last),xt);
end loop;
return res;
end dHt;
function dHx ( x : Vector; t : Complex_Number ) return Matrix is
xt : Vector(x'first..x'last+1);
m : Matrix(x'range,x'range);
begin
xt(x'range) := x;
xt(xt'last) := t;
for i in m'range(1) loop
for j in m'range(1) loop
m(i,j) := Eval(je(i,j),xt);
end loop;
end loop;
return m;
end dHx;
procedure Invert ( k : in integer; sols : in out Solution_List ) is
-- DESCRIPTION :
-- For all solutions s in sols : s.v(k) := 1/s.v(k).
tmp : Solution_List := sols;
begin
while not Is_Null(tmp) loop
declare
l : Link_to_Solution := Head_Of(tmp);
begin
l.v(k) := Create(1.0)/l.v(k);
l.t := Create(0.0);
end;
tmp := Tail_Of(tmp);
end loop;
end Invert;
procedure Cont is new Reporting_Continue(Max_Norm,Eval,dHt,dHx);
begin
Invert(k,sols);
Cont(file,sols,false);
Invert(k,sols);
Extract_Regular_Solutions(sols);
end;
Clear(h); Clear(hp); Clear(hpe); Clear(j); Clear(je);
end Mixed_Continue;
-- THE SOLVERS :
function One_Unknown_Solve ( p : Poly )
return Standard_Complex_Vectors.Vector is
-- DESCRIPTION :
-- Returns the solution vector of p, a polynomial in one unknown.
-- p will be solved by using the method of Durand-Kerner.
p1 : Poly := Re_Arrange(p);
init : Degrees := Initial_Degrees(p1);
index : integer := init'first;
min : integer := -Minimal_Degree(p1,index);
pv : Standard_Complex_Vectors.Vector(0..Degree(p1)+min);
z,res : Standard_Complex_Vectors.Vector(1..pv'last);
maxsteps : constant natural := 10;
eps : constant double_float := 10.0**(-10);
nb : integer := pv'last + 1;
procedure Store_Coeff ( t : in Term; cont : out boolean ) is
begin
nb := nb - 1;
if t.dg(index) = (nb - min)
then pv(nb) := t.cf;
else for i in reverse nb..(t.dg(index)+min+1) loop
pv(i) := Create(0.0);
end loop;
nb := t.dg(index) + min;
pv(nb) := t.cf;
end if;
cont := true;
end Store_Coeff;
procedure Polynomial_To_Vector is new Visiting_Iterator (Store_Coeff);
procedure Write ( step : in natural;
z,res : in Standard_Complex_Vectors.Vector ) is
begin
null; -- no output desired during the iterations
end Write;
procedure DuKe is new Durand_Kerner (Write);
begin
for i in pv'range loop -- initialize coefficient vector
pv(i) := Create(0.0);
end loop;
Standard_Integer_Vectors.Clear
(Standard_Integer_Vectors.Link_to_Vector(init));
Polynomial_To_Vector(p1);
Clear(p1);
for i in z'range loop
z(i) := Random;
res(i) := z(i);
end loop;
DuKe(pv,z,res,maxsteps,eps,nb);
return z;
end One_Unknown_Solve;
procedure One_Unknown_Solve ( p : in Poly; sols : in out Solution_List ) is
-- DESCRIPTION :
-- If p is a polynomial in one unknown,
-- p can be solved efficiently by the application of Durand-Kerner.
init : Degrees := Initial_Degrees(p);
function Make_Solutions ( x : in Standard_Complex_Vectors.Vector )
return Solution_List is
res,res_last : Solution_List;
s : Solution(1);
begin
s.m := 1;
s.t := Create(0.0);
for i in x'range loop
s.v(init'first) := x(i);
Append(res,res_last,s);
end loop;
return res;
end Make_Solutions;
begin
sols := Make_Solutions(One_Unknown_Solve(p));
Standard_Integer_Vectors.Clear
(Standard_Integer_Vectors.Link_to_Vector(init));
end One_Unknown_Solve;
procedure Two_Terms_Solve
( file : in file_type; p : in Laur_Sys;
tv : in Tree_of_Vectors; bkk : out natural;
sols : in out Solution_List ) is
-- DESCRIPTION :
-- The first polynomial of p consists of two terms.
-- A binomial system can be solved efficiently by
-- transforming and using de Moivre's rule.
d : Standard_Integer_Vectors.Link_to_Vector;
c : Complex_Number := Create(0.0);
k : natural := 0;
sols_last : Solution_List;
begin
-- put_line(file,"Applying Two_Terms_Solve on "); Write(file,p);
Binomial(p(p'first),d,k,c);
if k < d'first
then bkk := 0;
Standard_Integer_Vectors.Clear(d); return;
elsif ( c + Create(1.0) = Create(1.0) )
then bkk := 0;
Standard_Integer_Vectors.Clear(d); return;
elsif d(k) < 0
then Standard_Integer_Vectors.Min(d);
c := Create(1.0)/c;
end if;
declare
t : Transfo := Rotate(d,k);
tmp_bkk : natural := 0;
begin
Apply(t,d);
declare
v : Standard_Complex_Vectors.Vector(1..d(k));
tmp : Poly;
rtp : Laur_Sys(p'first..(p'last-1));
rtp_bkk : natural;
rtp_sols : Solution_List;
begin
for i in v'range loop
v(i) := Root(c,d(k),i);
for j in rtp'range loop
tmp := Transform(t,p(j+1));
rtp(j) := Evaluate(tmp,v(i),k);
Clear(tmp);
end loop;
Solve(file,rtp,tv,rtp_bkk,rtp_sols);
Clear(rtp);
tmp_bkk := tmp_bkk + rtp_bkk;
Insert(v(i),k,rtp_sols);
Transform(t,rtp_sols);
--Concat(sols,sols_last,rtp_sols);
Refine_and_Concat(file,p,rtp_sols,sols,sols_last);
end loop;
end;
Clear(t);
bkk := tmp_bkk;
end;
-- end if;
Standard_Integer_Vectors.Clear(d);
-- put_line(file,"The solutions found : "); put(file,sols);
end Two_Terms_Solve;
function Project_on_First_and_Solve
( p : Poly; k : integer; sols : Solution_List )
return Solution_List is
-- ON ENTRY :
-- p a Laurent polynomial in n unknowns x1,..,xk,..,xn;
-- sols contains values for x1,..,xn, except xk.
-- ON RETURN :
-- a solution list for p, obtained after substition of the values
-- for x1,..,xn into the polynomial p.
tmp,res,res_last : Solution_List;
init : Degrees := Initial_Degrees(p);
begin
-- put_line(file,"Calling Project_on_First_and_Solve");
-- put_line(file," with polynomial : ");
-- put(file,Laurent_Polynomial_to_Polynomial(p)); new_line(file);
tmp := sols;
while not Is_Null(tmp) loop
declare
p1 : Poly := Substitute(p,Head_Of(tmp).v,k);
sols1 : Solution_List;
begin
-- put(file,"k : "); put(file,k,1); new_line(file);
-- put(file,"v : "); put(file,Head_Of(tmp).v,3,3,3); new_line(file);
-- put_line(file,"After substitution : "); Write(file,p1);
if Number_of_Terms(p1) < 2
then null;
elsif Number_of_Terms(p1) = 2
then
declare
d : Standard_Integer_Vectors.Link_to_Vector;
l : natural := 0;
c : Complex_Number := Create(0.0);
begin
Binomial(p1,d,l,c);
if l < init'first
then null;
elsif ( c + Create(1.0) = Create(1.0) )
then null;
else
if d(l) < 0
then d(l) := -d(l); c := Create(1.0)/c;
end if;
declare
v : Standard_Complex_Vectors.Vector(1..d(l));
begin
for i in v'range loop
v(i) := Root(c,d(l),i);
end loop;
sols1 := Insert(v,k,Head_Of(tmp).all);
-- put_line(file,"Sols1 after Binomial :");
-- put(file,sols1);
Concat(res,res_last,sols1);
Shallow_Clear(sols1);
end;
end if;
Standard_Integer_Vectors.Clear(d);
end;
else
sols1 := Insert(One_Unknown_Solve(p1),k,Head_Of(tmp).all);
Concat(res,res_last,sols1);
-- put_line(file,"Sols1 :"); put(file,sols1);
Shallow_Clear(sols1);
end if;
Clear(p1);
end;
tmp := Tail_Of(tmp);
end loop;
Standard_Integer_Vectors.Clear
(Standard_Integer_Vectors.Link_to_Vector(init));
return res;
end Project_on_First_and_Solve;
procedure Project_and_Solve
( file : in file_type; p : in Laur_Sys; k : in integer;
m : in out Standard_Integer_Vectors.Vector;
nd : in node; bkk : out natural;
sols : in out Solution_List ) is
-- DESCRIPTION :
-- Solves the projected start system of p along a direction v.
-- ON ENTRY :
-- file a file to write intermediate results on;
-- p a Laurent polynomial system;
-- k entry in the degrees of p;
-- m m(m'first) equals Maximal_Support(p(p'first),v) > 0;
-- nd a node in the tree of useful directions.
-- ON RETURN :
-- m m(m'first+1..m'last) contains the maximal degrees
-- of the last n-1 equations of p in xk;
-- bkk the BKK bound of the projected system;
-- sols the solutions of the projected system.
g_v : Laur_Sys(p'first..(p'last-1));
bkk_g_v : natural;
sols_g_v : Solution_List;
begin
-- put_line(file,"Applying Project_and_Solve on"); Write(file,p);
for i in g_v'range loop
m(i+1) := Maximal_Degree(p(i+1),k);
g_v(i) := Face(k,m(i+1),p(i+1));
Reduce(k,g_v(i));
end loop;
if (nd.ltv = null) or else Is_Null(nd.ltv.all)
then Solve(file,g_v,bkk_g_v,sols_g_v);
else Solve(file,g_v,nd.ltv.all,bkk_g_v,sols_g_v);
end if;
-- put(file,"After Solve (without tv) bkk_g_v = "); put(file,bkk_g_v,1);
-- new_line(file);
declare
p0 : Poly := Re_Arrange(p(p'first));
p1 : Poly := Face(k,m(m'first),p0);
cnst : Term;
begin
cnst.dg := new Standard_Integer_Vectors.Vector'(p'range => 0);
if Coeff(p1,cnst.dg) = Create(0.0)
then cnst.cf := Coeff(p0,cnst.dg);
Add(p1,cnst);
end if;
Standard_Integer_Vectors.Clear
(Standard_Integer_Vectors.Link_to_Vector(cnst.dg));
Clear(p0);
sols := Project_on_First_and_Solve(p1,k,sols_g_v);
-- sols := Project_on_First_and_Solve(file,p1,k,sols_g_v);
Set_Continuation_Parameter(sols,Create(0.0));
Clear(p1);
end;
bkk := m(m'first)*bkk_g_v;
Clear(sols_g_v);
Clear(g_v);
-- put_line(file,"The solutions found : "); put(file,sols);
end Project_and_Solve;
procedure Unmixed_Solve
( file : in file_type; p : in Laur_Sys; dl : in List;
tv : in Tree_of_Vectors;
bkk : out natural; sols : in out Solution_List ) is
-- DESCRIPTION :
-- Solves a Laurent polynomial system where all polytopes are the same.
-- ON ENTRY :
-- file a file to write intermediate results on;
-- p a Laurent polynomial system;
-- dl the list of powers for p;
-- tv the tree of degrees containing the useful directions.
-- ON RETURN :
-- bkk the bkk bound;
-- sols the list of solutions.
sols_last : Solution_List;
tmp_bkk : natural := 0;
tmp : Tree_of_Vectors := tv;
begin
tmp_bkk := 0;
tmp := tv;
while not Is_Null(tmp) loop
declare
nd : node := Head_Of(tmp);
v : Standard_Integer_Vectors.Link_to_Vector := nd.d;
i : integer := Pivot(v);
pv : integer := Maximal_Support(dl,v.all);
t : Transfo := Build_Transfo(v,i);
tp : Laur_Sys(p'range) := Transform(t,p);
bkk_tp : natural;
sols_tp : Solution_List;
max : Standard_Integer_Vectors.Vector(p'range);
begin
Write_Direction(file,v);
max(max'first) := pv;
-- if (nd.ltv = null) or else Is_Null(nd.ltv.all)
-- then Projected_Solve(file,tp,i,max,bkk_tp,sols_tp);
-- else Projected_Solve
-- (file,tp,i,max,nd.ltv.all,bkk_tp,sols_tp);
-- end if;
Project_and_Solve(file,tp,i,max,nd,bkk_tp,sols_tp);
Mixed_Continue(file,tp,i,max,sols_tp);
tmp_bkk := tmp_bkk + bkk_tp;
Transform(t,sols_tp);
--Concat(sols,sols_last,sols_tp);
Refine_and_Concat(file,p,sols_tp,sols,sols_last);
Clear(t); Clear(tp);
end;
tmp := Tail_Of(tmp);
end loop;
bkk := tmp_bkk;
end Unmixed_Solve;
procedure Unmixed_Solve
( file : in file_type; p : in Laur_Sys; dl : in List;
bkk : out natural; sols : in out Solution_List ) is
tv : Tree_of_Vectors;
begin
Volumes.Volume(p'last,dl,tv,bkk);
Unmixed_Solve(file,p,dl,tv,bkk,sols);
Clear(tv);
end Unmixed_Solve;
procedure Mixed_Solve
( file : in file_type; p : in Laur_Sys;
adl : in out Array_of_Lists; tv : in Tree_of_Vectors;
bkk : out natural; sols : in out Solution_List ) is
-- DESCRIPTION :
-- Computes the solutions of the Laurent polynomial system p,
-- where p has more than one equation.
-- NOTE :
-- This procedure mirrors the procedure Minkowski_Sum in the body
-- of the package Volumes.
tmp_bkk,len : natural;
tmp : Tree_of_Vectors;
index : integer := Index2(adl);
wp : Laur_Sys(p'range);
sols_last : Solution_List;
shifted : boolean;
perm,mix : Standard_Integer_Vectors.Link_to_Vector;
begin
Interchange2(adl,index);
len := Length_Of(adl(adl'first));
-- put_line(file,"Applying Mixed_Solve on"); Write(file,p);
if len = 2
then wp := Interchange2(p,index);
Two_Terms_Solve(file,wp,tv,bkk,sols);
else
if len > 2
then -- INITIALIZATION :
Mixture(adl,perm,mix);
wp := Permute(perm.all,p);
declare
zeroes : Degrees
:= new Standard_Integer_Vectors.Vector'(p'range => 0);
tmpwpi : Poly;
begin
if Coeff(wp(wp'first),zeroes) = Create(0.0)
then shifted := true;
-- wp(wp'first) := Shift(p(p'first));
Copy(p(index),tmpwpi); wp(wp'first) := tmpwpi;
Shift(wp(wp'first));
else shifted := false;
end if;
Standard_Integer_Vectors.Clear
(Standard_Integer_Vectors.Link_to_Vector(zeroes));
end;
-- MIXED HOMOTOPY CONTINUATION :
tmp_bkk := 0;
tmp := tv;
while not Is_Null(tmp) loop
declare
nd : node := Head_Of(tmp);
v : Standard_Integer_Vectors.Link_to_Vector := nd.d;
k : integer := Pivot(v);
pv : integer := Maximal_Support(wp(wp'first),v);
t : Transfo := Build_Transfo(v,k);
twp : Laur_Sys(wp'range) := Transform(t,wp);
bkk_twp : natural;
sols_twp : Solution_List;
m : Standard_Integer_Vectors.Vector(wp'range);
begin
Write_Direction(file,v);
m(m'first) := pv;
-- if (nd.ltv = null) or else Is_Null(nd.ltv.all)
-- then Projected_Solve(file,twp,k,m,bkk_twp,sols_twp);
-- else Projected_Solve
-- (file,twp,k,m,nd.ltv.all,bkk_twp,sols_twp);
-- end if;
Project_and_Solve(file,twp,k,m,nd,bkk_twp,sols_twp);
Mixed_Continue(file,twp,k,m,sols_twp);
tmp_bkk := tmp_bkk + bkk_twp;
Transform(t,sols_twp);
--Concat(sols,sols_last,sols_twp);
Refine_and_Concat(file,wp,sols_twp,sols,sols_last);
Clear(t); Clear(twp);
end;
tmp := Tail_Of(tmp);
end loop;
bkk := tmp_bkk;
Standard_Integer_Vectors.Clear(perm);
Standard_Integer_Vectors.Clear(mix);
if shifted
then Clear(wp(wp'first));
end if;
else -- len < 2
bkk := 0;
end if;
end if;
end Mixed_Solve;
procedure Mixed_Solve
( file : in file_type; p : in Laur_Sys;
adl : in out Array_of_Lists; bkk : out natural;
sols : in out Solution_List ) is
tv : Tree_of_Vectors;
begin
Volumes.Mixed_Volume(adl'last,adl,tv,bkk);
Mixed_Solve(file,p,adl,tv,bkk,sols);
Clear(tv);
end Mixed_Solve;
procedure Solve ( file : in file_type; p : in Laur_Sys;
bkk : out natural; sols : in out Solution_List ) is
al : Array_of_Lists(p'range) := Create(p);
tv : Tree_of_Vectors;
begin
Volumes.Mixed_Volume(p'last,al,tv,bkk);
Solve(file,p,tv,bkk,sols);
Deep_Clear(al); Clear(tv);
end Solve;
procedure Solve ( file : in file_type; p : in Laur_Sys;
tv : in Tree_of_Vectors; bkk : out natural;
sols : in out Solution_List ) is
-- NOTE :
-- This procedure mirrors the procedure Volumes.Mixed_Volume,
-- with a tree of useful directions on entry.
begin
if p'first = p'last
then One_Unknown_Solve(p(p'first),sols);
bkk := Length_Of(sols);
else --if Is_Fewnomial_System(p)
-- then
-- declare
-- fail : boolean;
-- begin
-- Fewnomials.Solve(p,sols,fail);
-- if fail
-- then bkk := 0; Clear(sols);
-- else bkk := Length_Of(sols);
-- end if;
-- end;
-- else
declare
adl : Array_of_Lists(p'range) := Create(p);
begin
if All_Equal(adl)
then
for i in (adl'first+1)..adl'last loop
Deep_Clear(adl(i));
end loop;
declare
wp : Laur_Sys(p'range);
shifted : boolean;
begin
Normalize(p,adl(adl'first),wp,shifted);
if Is_Null(tv)
then Unmixed_Solve(file,wp,adl(adl'first),bkk,sols);
else Unmixed_Solve(file,wp,adl(adl'first),tv,bkk,sols);
end if;
if shifted
then Clear(wp);
end if;
end;
elsif Is_Null(tv)
then Mixed_Solve(file,p,adl,bkk,sols);
else Mixed_Solve(file,p,adl,tv,bkk,sols);
end if;
end;
end if;
end Solve;
end Mixed_Homotopy_Continuation;