File: [local] / OpenXM_contrib / PHC / Ada / Schubert / numeric_minor_equations.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_Poly_Functions; use Standard_Complex_Poly_Functions;
with Standard_Complex_Poly_SysFun; use Standard_Complex_Poly_SysFun;
with Brackets; use Brackets;
with Plane_Representations; use Plane_Representations;
with Evaluated_Minors; use Evaluated_Minors;
with Symbolic_Minor_Equations; use Symbolic_Minor_Equations;
package body Numeric_Minor_Equations is
tol : constant double_float := 10.0**(-10);
-- EXPANDING ACCORDING A BRACKET MONOMIAL :
function Expanded_Minors
( cffmat : Standard_Floating_Matrices.Matrix;
polmat : Standard_Complex_Poly_Matrices.Matrix;
bm : Bracket_Monomial ) return Poly is
res : Poly := Null_Poly;
first : boolean := true;
procedure Visit_Bracket ( b : in Bracket; continue : out boolean ) is
factor : double_float;
minor : Poly;
begin
if first
then declare
bb : Bracket(b'first..b'last-1) := b(b'first+1..b'last);
begin
factor := Determinant(cffmat,bb);
end;
first := false;
else minor := Expanded_Minor(polmat,b);
if (minor /= Null_Poly) and (abs(factor) > tol)
then Mul(minor,Create(factor));
Add(res,minor);
end if;
Clear(minor);
end if;
continue := true;
end Visit_Bracket;
procedure Visit_Brackets is new Enumerate_Brackets(Visit_Bracket);
begin
Visit_Brackets(bm);
return res;
end Expanded_Minors;
function Expanded_Minors
( cffmat : Standard_Complex_Matrices.Matrix;
polmat : Standard_Complex_Poly_Matrices.Matrix;
bm : Bracket_Monomial ) return Poly is
res : Poly := Null_Poly;
first : boolean := true;
factor : Complex_Number;
procedure Visit_Bracket ( b : in Bracket; continue : out boolean ) is
minor : Poly;
begin
if first
then declare
bb : Bracket(b'first..b'last-1) := b(b'first+1..b'last);
begin
factor := Determinant(cffmat,bb);
end;
first := false;
else minor := Expanded_Minor(polmat,b);
if (minor /= Null_Poly) and (AbsVal(factor) > tol)
then Mul(minor,factor);
Add(res,minor);
end if;
Clear(minor);
end if;
continue := true;
end Visit_Bracket;
procedure Visit_Brackets is new Enumerate_Brackets(Visit_Bracket);
begin
Visit_Brackets(bm);
return res;
end Expanded_Minors;
function Expanded_Minors
( cntmat,polmat : Standard_Complex_Poly_Matrices.Matrix;
bm : Bracket_Monomial ) return Poly is
res : Poly := Null_Poly;
first : boolean := true;
factor : Poly;
procedure Visit_Bracket ( b : in Bracket; continue : out boolean ) is
minor : Poly;
begin
if first
then declare
bb : Bracket(b'first..b'last-1) := b(b'first+1..b'last);
begin
factor := Expanded_Minor(cntmat,bb);
end;
first := false;
else minor := Expanded_Minor(polmat,b);
if (minor /= Null_Poly) and (factor /= Null_Poly)
then Mul(minor,factor);
Add(res,minor);
end if;
Clear(factor); Clear(minor);
end if;
continue := true;
end Visit_Bracket;
procedure Visit_Brackets is new Enumerate_Brackets(Visit_Bracket);
begin
Visit_Brackets(bm);
return res;
end Expanded_Minors;
function Lifted_Expanded_Minors
( cntmat,polmat : Standard_Complex_Poly_Matrices.Matrix;
bm : Bracket_Monomial ) return Poly is
res : Poly := Null_Poly;
first : boolean := true;
factor : Poly;
procedure Visit_Bracket ( b : in Bracket; continue : out boolean ) is
minor,extmin : Poly;
begin
if first
then declare
bb : Bracket(b'first..b'last-1) := b(b'first+1..b'last);
begin
factor := Expanded_Minor(cntmat,bb);
end;
first := false;
else minor := Expanded_Minor(polmat,b);
if (minor /= Null_Poly) and (factor /= Null_Poly)
then extmin := Extend_Zero_Lifting(minor);
Mul(extmin,factor);
Add(res,extmin);
end if;
Clear(factor); Clear(minor); Clear(extmin);
end if;
continue := true;
end Visit_Bracket;
procedure Visit_Brackets is new Enumerate_Brackets(Visit_Bracket);
begin
Visit_Brackets(bm);
return res;
end Lifted_Expanded_Minors;
-- EXPANDING ACCORDING A BRACKET POLYNOMIAL :
function Expanded_Minors
( cffmat : Standard_Floating_Matrices.Matrix;
polmat : Standard_Complex_Poly_Matrices.Matrix;
bp : Bracket_Polynomial ) return Poly is
res : Poly := Null_Poly;
procedure Visit_Term ( t : in Bracket_Term; continue : out boolean ) is
minor : Poly := Expanded_Minors(cffmat,polmat,t.monom);
begin
if REAL_PART(t.coeff) > 0.0
then Add(res,minor);
else Sub(res,minor);
end if;
Clear(minor);
continue := true;
end Visit_Term;
procedure Visit_Terms is new Enumerate_Terms(Visit_Term);
begin
Visit_Terms(bp);
return res;
end Expanded_Minors;
function Expanded_Minors
( cffmat : Standard_Complex_Matrices.Matrix;
polmat : Standard_Complex_Poly_Matrices.Matrix;
bp : Bracket_Polynomial ) return Poly is
res : Poly := Null_Poly;
procedure Visit_Term ( t : in Bracket_Term; continue : out boolean ) is
minor : Poly := Expanded_Minors(cffmat,polmat,t.monom);
begin
if REAL_PART(t.coeff) > 0.0
then Add(res,minor);
else Sub(res,minor);
end if;
Clear(minor);
continue := true;
end Visit_Term;
procedure Visit_Terms is new Enumerate_Terms(Visit_Term);
begin
Visit_Terms(bp);
return res;
end Expanded_Minors;
function Expanded_Minors
( cntmat,polmat : Standard_Complex_Poly_Matrices.Matrix;
bp : Bracket_Polynomial ) return Poly is
res : Poly := Null_Poly;
procedure Visit_Term ( t : in Bracket_Term; continue : out boolean ) is
minor : Poly := Expanded_Minors(cntmat,polmat,t.monom);
begin
if REAL_PART(t.coeff) > 0.0
then Add(res,minor);
else Sub(res,minor);
end if;
Clear(minor);
continue := true;
end Visit_Term;
procedure Visit_Terms is new Enumerate_Terms(Visit_Term);
begin
Visit_Terms(bp);
return res;
end Expanded_Minors;
function Lifted_Expanded_Minors
( cntmat,polmat : Standard_Complex_Poly_Matrices.Matrix;
bp : Bracket_Polynomial ) return Poly is
res : Poly := Null_Poly;
procedure Visit_Term ( t : in Bracket_Term; continue : out boolean ) is
minor : Poly := Lifted_Expanded_Minors(cntmat,polmat,t.monom);
begin
if REAL_PART(t.coeff) > 0.0
then Add(res,minor);
else Sub(res,minor);
end if;
Clear(minor);
continue := true;
end Visit_Term;
procedure Visit_Terms is new Enumerate_Terms(Visit_Term);
begin
Visit_Terms(bp);
return res;
end Lifted_Expanded_Minors;
-- EXPANDING TO CONSTRUCT POLYNOMIAL SYSTEMS :
function Expanded_Minors
( cffmat : Standard_Floating_Matrices.Matrix;
polmat : Standard_Complex_Poly_Matrices.Matrix;
bs : Bracket_System ) return Poly_Sys is
res : Poly_Sys(bs'first+1..bs'last);
begin
for i in res'range loop
res(i) := Expanded_Minors(cffmat,polmat,bs(i));
end loop;
return res;
end Expanded_Minors;
function Expanded_Minors
( cffmat : Standard_Complex_Matrices.Matrix;
polmat : Standard_Complex_Poly_Matrices.Matrix;
bs : Bracket_System ) return Poly_Sys is
res : Poly_Sys(bs'first+1..bs'last);
begin
for i in res'range loop
res(i) := Expanded_Minors(cffmat,polmat,bs(i));
end loop;
return res;
end Expanded_Minors;
function Expanded_Minors
( cntmat,polmat : Standard_Complex_Poly_Matrices.Matrix;
bs : Bracket_System ) return Poly_Sys is
res : Poly_Sys(bs'first+1..bs'last);
begin
for i in res'range loop
res(i) := Expanded_Minors(cntmat,polmat,bs(i));
end loop;
return res;
end Expanded_Minors;
function Lifted_Expanded_Minors
( cntmat,polmat : Standard_Complex_Poly_Matrices.Matrix;
bs : Bracket_System ) return Poly_Sys is
res : Poly_Sys(bs'first+1..bs'last);
begin
for i in res'range loop
res(i) := Lifted_Expanded_Minors(cntmat,polmat,bs(i));
end loop;
return res;
end Lifted_Expanded_Minors;
function Evaluate ( p : Poly; x : Standard_Complex_Matrices.Matrix )
return Complex_Number is
xv : constant Standard_Complex_Vectors.Vector := Vector_Rep(x);
begin
return Eval(p,xv);
end Evaluate;
function Evaluate ( p : Poly_Sys; x : Standard_Complex_Matrices.Matrix )
return Standard_Complex_Vectors.Vector is
xv : constant Standard_Complex_Vectors.Vector := Vector_Rep(x);
begin
return Eval(p,xv);
end Evaluate;
procedure Embed ( t : in out Term ) is
dg : Degrees
:= new Standard_Natural_Vectors.Vector(t.dg'first..t.dg'last+1);
begin
dg(t.dg'range) := t.dg.all;
dg(dg'last) := 0;
Clear(t.dg);
t.dg := dg;
end Embed;
function Embed ( t : Term ) return Term is
res : Term;
begin
res.dg := new Standard_Natural_Vectors.Vector(t.dg'first..t.dg'last+1);
res.dg(t.dg'range) := t.dg.all;
res.dg(res.dg'last) := 0;
res.cf := t.cf;
return res;
end Embed;
procedure Embed ( p : in out Poly ) is
res : Poly := Null_Poly;
procedure Embed_Term ( t : in Term; continue : out boolean ) is
et : Term := Embed(t);
begin
Add(res,et);
Clear(et);
continue := true;
end Embed_Term;
procedure Embed_Terms is new Visiting_Iterator(Embed_Term);
begin
Embed_Terms(p);
Clear(p);
p := res;
end Embed;
procedure Embed ( p : in out Poly_Sys ) is
begin
for i in p'range loop
Embed(p(i));
end loop;
end Embed;
procedure Embed ( m : in out Standard_Complex_Poly_Matrices.Matrix ) is
begin
for i in m'range(1) loop
for j in m'range(2) loop
if m(i,j) /= Null_Poly
then Embed(m(i,j));
end if;
end loop;
end loop;
end Embed;
function Linear_Homotopy ( target,start : Poly ) return Poly is
res : Poly := Null_Poly;
procedure Embed_Target_Term ( t : in Term; continue : out boolean ) is
et : Term := Embed(t);
begin
et.dg(et.dg'last) := 1;
Add(res,et);
Clear(et);
continue := true;
end Embed_Target_Term;
procedure Embed_Target_Terms is new Visiting_Iterator(Embed_Target_Term);
procedure Embed_Start_Term ( t : in Term; continue : out boolean ) is
et : Term := Embed(t);
begin
Add(res,et);
et.dg(et.dg'last) := 1;
Sub(res,et);
Clear(et);
continue := true;
end Embed_Start_Term;
procedure Embed_Start_Terms is new Visiting_Iterator(Embed_Start_Term);
begin
Embed_Target_Terms(target);
Embed_Start_Terms(start);
return res;
end Linear_Homotopy;
function Linear_Interpolation
( target,start : Poly; k : natural ) return Poly is
res : Poly := Null_Poly;
procedure Embed_Target_Term ( t : in Term; continue : out boolean ) is
et : Term;
begin
Copy(t,et);
et.dg(k) := et.dg(k) + 1; -- multiply with t
Add(res,et);
Clear(et);
continue := true;
end Embed_Target_Term;
procedure Embed_Target_Terms is new Visiting_Iterator(Embed_Target_Term);
procedure Embed_Start_Term ( t : in Term; continue : out boolean ) is
et : Term;
begin
Copy(t,et);
Add(res,et); -- res := res + et
et.dg(k) := et.dg(k) + 1; -- multiply with t
Sub(res,et); -- res := res + et - t*et
Clear(et);
continue := true;
end Embed_Start_Term;
procedure Embed_Start_Terms is new Visiting_Iterator(Embed_Start_Term);
begin
Embed_Target_Terms(target);
Embed_Start_Terms(start);
return res;
end Linear_Interpolation;
procedure Divide_Common_Factor ( p : in out Poly; k : in natural ) is
first : boolean := true;
min : natural;
procedure Min_Power ( t : in Term; continue : out boolean ) is
begin
if first
then first := false;
min := t.dg(k); -- initialize minimal power
else if t.dg(k) < min
then min := t.dg(k);
end if;
end if;
continue := true;
end Min_Power;
procedure Find_Min_Power is new Visiting_Iterator(Min_Power);
procedure Divide ( t : in out Term; continue : out boolean ) is
begin
t.dg(k) := t.dg(k) - min;
continue := true;
end Divide;
procedure Divide_Min_Power is new Changing_Iterator(Divide);
begin
Find_Min_Power(p);
if min > 0
then Divide_Min_Power(p);
end if;
end Divide_Common_Factor;
end Numeric_Minor_Equations;