File: [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Polynomials / generic_laur_poly_functions.adb (download)
Revision 1.1.1.1 (vendor branch), Sun Oct 29 17:45:26 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 unchecked_deallocation;
with Standard_Integer_Vectors;
with Standard_Integer_Matrices; use Standard_Integer_Matrices;
package body Generic_Laur_Poly_Functions is
-- NOTE : Evaluation is not guaranteed to work for negative exponents...
-- DATA STRUCTURES :
MAX_INT : constant natural := 100000;
type kind is (coefficient,polynomial);
type Poly_Rec ( k : kind := coefficient ) is record
case k is
when coefficient => c : number;
when polynomial => p : Eval_Poly;
end case;
end record;
type Coeff_Poly_Rec ( k : kind := coefficient ) is record
case k is
when coefficient => c : integer;
when polynomial => p : Eval_Coeff_Poly;
end case;
end record;
type Eval_Poly_Rep is array(integer range <>) of Poly_Rec;
type Eval_Coeff_Poly_Rep is array(integer range <>) of Coeff_Poly_Rec;
procedure free is new unchecked_deallocation(Eval_Poly_Rep,Eval_Poly);
procedure free is
new unchecked_deallocation(Eval_Coeff_Poly_Rep,Eval_Coeff_Poly);
-- AUXILIARY OPERATIONS :
function Convert ( c : number; n : natural ) return natural is
-- DESCRIPTION :
-- Returns the corresponding value for c, when it lies in 1..n,
-- otherwise 0 is returned.
begin
for i in 1..n loop
if c = Create(i)
then return i;
end if;
end loop;
return 0;
end Convert;
function Head_Term ( p : Poly ) return Term is
-- DESCRIPTION :
-- Returns the leading term of the polynomial.
res : Term;
procedure Scan_Term ( t : in Term; continue : out boolean ) is
begin
res := t;
continue := false;
end Scan_Term;
procedure Scan_Terms is new Visiting_Iterator(Scan_Term);
begin
Scan_Terms(p);
return res;
end Head_Term;
procedure Initialize ( evpr : in out Eval_Poly_Rep ) is
begin
for i in evpr'range loop
declare
nullpr : Poly_Rec(polynomial);
begin
nullpr.p := null;
evpr(i) := nullpr;
end;
end loop;
end Initialize;
procedure Initialize ( evpr : in out Eval_Coeff_Poly_Rep ) is
begin
for i in evpr'range loop
declare
nullpr : Coeff_Poly_Rec(polynomial);
begin
nullpr.p := null;
evpr(i) := nullpr;
end;
end loop;
end Initialize;
procedure Initialize ( p : in Poly; cff : out Vector; deg : out Matrix ) is
-- DESCRIPTION :
-- Returns a vector/matrix representation of the polynomial p.
-- Starts filling in backwards, since the highest degrees are in front
-- of the list. This could reduce the sorting time later.
ind : natural := cff'last+1;
procedure Scan_Term ( t : in Term; continue : out boolean ) is
begin
ind := ind-1;
Copy(t.cf,cff(ind));
for j in t.dg'range loop
deg(ind,j) := t.dg(j);
end loop;
continue := true;
end Scan_Term;
procedure Scan_Terms is new Visiting_Iterator(Scan_Term);
begin
Scan_Terms(p);
end Initialize;
procedure Swap ( cff : in out Vector; deg : in out Matrix;
i,j,k : in natural ) is
-- DESCRIPTION :
-- Swaps the ith row with the jth row, starting from the kth column.
tmpcff : number;
tmpdeg : natural;
begin
Copy(cff(i),tmpcff);
Copy(cff(j),cff(i));
Copy(tmpcff,cff(j));
Clear(tmpcff);
for kk in k..deg'last(2) loop
tmpdeg := deg(i,kk);
deg(i,kk) := deg(j,kk);
deg(j,kk) := tmpdeg;
end loop;
end Swap;
procedure Sort ( cff : in out Vector; deg : in out Matrix;
i1,i2,k : in natural ) is
-- DESCRIPTION :
-- Sorts the elements in the kth column of the degree matrix, in
-- the range i1..i2. The coefficient vector gets sorted along.
ind,min : natural;
begin
for i in i1..i2 loop -- sort by swapping minimal element
min := deg(i,k);
ind := i;
for j in i+1..i2 loop -- search for minimal element
if deg(j,k) < min
then min := deg(j,k);
ind := j;
end if;
end loop;
if ind /= i -- swap cff and deg
then Swap(cff,deg,i,ind,k);
end if;
end loop;
end Sort;
procedure Create ( cff : in out Vector; deg : in out Matrix;
i1,i2,k : in natural; ep : out Eval_Poly ) is
-- DESCRIPTION :
-- Returns in ep a nested Horner scheme to evaluate a polynomial given
-- in vector/matrix representation with coefficients in cff and degrees
-- in deg. The range being considered is i1..i2, from the kth column.
-- REQUIRED :
-- The entries in the kth column of deg in i1..i2 are sorted
-- in increasing order.
evpr : Eval_Poly_Rep(0..deg(i2,k));
ind,j1,j2 : natural;
begin
Initialize(evpr);
if k = deg'last(2) -- polynomial in one unknown
then for i in i1..i2 loop
declare
pr : Poly_Rec(Coefficient);
begin
Copy(cff(i),pr.c);
evpr(deg(i,k)) := pr;
end;
end loop;
else ind := i1; -- recursive call
while ind <= i2 loop
j1 := ind; j2 := ind;
while j2 < i2 and then deg(j1,k) = deg(j2+1,k) loop
j2 := j2 + 1;
end loop;
declare
pr : Poly_Rec(Polynomial);
begin
Sort(cff,deg,j1,j2,k+1);
Create(cff,deg,j1,j2,k+1,pr.p);
evpr(deg(ind,k)) := pr;
end;
ind := j2+1;
end loop;
end if;
ep := new Eval_Poly_Rep'(evpr);
end Create;
-- CONSTRUCTORS :
function Create ( p : Poly; n : natural; d : integer ) return Eval_Poly is
-- DESCRIPTION :
-- An evaluable polynomial is returned for p, with d = Maximal_Degree(p,x1)
-- and n = Number_of_Unknowns(p) >= 1.
res : Eval_Poly;
evpr : Eval_Poly_Rep(0..d);
terms : array(0..d) of Poly := (0..d => Null_Poly);
procedure Add_Term1 ( t : in Term; cont : out boolean ) is
pr : Poly_Rec(coefficient);
begin
copy(t.cf,pr.c);
evpr(t.dg(t.dg'first)) := pr;
cont := true;
end Add_Term1;
procedure Add_Terms1 is new Visiting_Iterator(Add_Term1);
procedure Add_Term ( t : in Term; cont : out boolean ) is
nt : Term;
begin
nt.cf := t.cf;
nt.dg := new Standard_Integer_Vectors.Vector(t.dg'first..t.dg'last-1);
for i in nt.dg'range loop
nt.dg(i) := t.dg(i+1);
end loop;
Clear(nt);
Add(terms(t.dg(t.dg'first)),nt);
cont := true;
end Add_Term;
procedure Add_Terms is new Visiting_Iterator(Add_Term);
begin
Initialize(evpr);
if n = 1
then Add_Terms1(p);
else Add_Terms(p);
for i in terms'range loop
declare
pr : Poly_Rec(polynomial);
begin
pr.p := Create(terms(i));
evpr(i) := pr;
Clear(terms(i));
end;
end loop;
end if;
res := new Eval_Poly_Rep'(evpr);
return res;
end Create;
function Create ( p : Poly; n,nb : natural; d : integer )
return Eval_Coeff_Poly is
-- DESCRIPTION :
-- An evaluable polynomial is returned for p, with d = Maximal_Degree(p,x1),
-- n = Number_of_Unknowns(p) >= 1 and nb = Number_of_Terms(p).
-- The coefficients of p are converted natural numbers.
res : Eval_Coeff_Poly;
evpr : Eval_Coeff_Poly_Rep(0..d);
terms : array(0..d) of Poly := (0..d => Null_Poly);
procedure Add_Term1 ( t : in Term; cont : out boolean ) is
pr : Coeff_Poly_Rec(coefficient);
begin
pr.c := Convert(t.cf,nb);
evpr(t.dg(t.dg'first)) := pr;
cont := true;
end Add_Term1;
procedure Add_Terms1 is new Visiting_Iterator(Add_Term1);
procedure Add_Term ( t : in Term; cont : out boolean ) is
nt : Term;
begin
nt.cf := t.cf;
nt.dg := new Standard_Integer_Vectors.Vector(t.dg'first..t.dg'last-1);
for i in nt.dg'range loop
nt.dg(i) := t.dg(i+1);
end loop;
Add(terms(t.dg(t.dg'first)),nt);
Clear(nt);
cont := true;
end Add_Term;
procedure Add_Terms is new Visiting_Iterator(Add_Term);
begin
Initialize(evpr);
if n = 1
then for i in evpr'range loop -- initialization
declare
nullpr : Coeff_Poly_Rec(polynomial);
begin
nullpr.p := null;
evpr(i) := nullpr;
end;
end loop;
Add_Terms1(p);
else Add_Terms(p);
for i in terms'range loop
declare
pr : Coeff_Poly_Rec(polynomial);
ind,mdi : integer;
begin
if terms(i) = Null_Poly
then pr.p := null;
else ind := Head_Term(terms(i)).dg'first;
mdi := Maximal_Degree(terms(i),ind);
pr.p := Create(terms(i),n-1,nb,mdi);
end if;
evpr(i) := pr;
Clear(terms(i));
end;
end loop;
end if;
res := new Eval_Coeff_Poly_Rep'(evpr);
return res;
end Create;
function Create ( p : Poly ) return Eval_Poly is
res : Eval_Poly;
nbvar : constant natural := Number_of_Unknowns(p);
nbtms : constant natural := Number_of_Terms(p);
cff : Vector(1..nbtms);
deg : Matrix(1..nbtms,1..nbvar);
begin
if (p = Null_Poly) or else (nbvar = 0)
then return null;
else Initialize(p,cff,deg);
Sort(cff,deg,1,nbtms,1);
Create(cff,deg,1,nbtms,1,res);
Clear(cff);
return res;
end if;
end Create;
function Create ( p : Poly ) return Eval_Coeff_Poly is
res : Eval_Coeff_Poly;
lp : Poly := Null_Poly;
n : constant natural := Number_of_Unknowns(p);
nb : constant natural := Number_of_Terms(p);
cnt : natural := 0;
ind : integer;
procedure Label_Term ( t : in Term; cont : out boolean ) is
lt : Term;
begin
cnt := cnt + 1;
lt.cf := Create(cnt);
lt.dg := new Standard_Integer_Vectors.Vector'(t.dg.all);
Add(lp,lt);
Clear(lt);
cont := true;
end Label_Term;
procedure Label_Terms is new Visiting_Iterator(Label_Term);
begin
if (p = Null_Poly) or else (nb = 0)
then return null;
else Label_Terms(p);
ind := Head_Term(p).dg'first;
res := Create(lp,n,nb,Maximal_Degree(p,ind));
Clear(lp);
end if;
return res;
end Create;
procedure Diff ( p : in Poly; i : in integer;
cp : out Eval_Coeff_Poly; m : out Vector ) is
nb : constant natural := Number_of_Terms(p);
n : constant natural := Number_of_Unknowns(p);
ind,cnt : integer;
dp : Poly := Null_Poly;
procedure Diff_Term ( t : in Term; cont : out boolean ) is
dt : Term;
begin
cnt := cnt + 1;
if t.dg(i) > 0
then dt.cf := Create(cnt);
dt.dg := new Standard_Integer_Vectors.Vector'(t.dg.all);
m(cnt) := Create(t.dg(i));
dt.dg(i) := dt.dg(i) - 1;
Add(dp,dt);
Clear(dt);
else m(cnt) := Create(0);
end if;
cont := true;
end Diff_Term;
procedure Diff_Terms is new Visiting_Iterator(Diff_Term);
begin
cnt := 0;
Diff_Terms(p);
if dp /= Null_Poly
then ind := Head_Term(dp).dg'first;
cp := Create(dp,n,nb,Maximal_Degree(dp,ind));
end if;
end Diff;
function Coeff ( p : Poly ) return Vector is
res : Vector(1..Number_of_Terms(p));
cnt : natural := 0;
procedure Collect_Term ( t : in Term; cont : out boolean ) is
begin
cnt := cnt + 1;
copy(t.cf,res(cnt));
cont := true;
end Collect_Term;
procedure Collect_Terms is new Visiting_Iterator(Collect_Term);
begin
Collect_Terms(p);
return res;
end Coeff;
-- EVALUATORS :
function Eval ( p : Poly; x : number; i : integer ) return Poly is
res : Poly := Null_Poly;
procedure Eval_Term ( t : in Term; cont : out boolean ) is
nt : Term;
begin
copy(t.cf,nt.cf);
nt.dg := new Standard_Integer_Vectors.Vector(t.dg'first..t.dg'last-1);
for j in t.dg'range loop
if j < i
then nt.dg(j) := t.dg(j);
elsif j > i
then nt.dg(j-1) := t.dg(j);
else for k in 1..t.dg(i) loop
Mul(nt.cf,x);
end loop;
end if;
end loop;
Add(res,nt);
Clear(nt);
cont := true;
end Eval_Term;
procedure Eval_Terms is new Visiting_Iterator(Eval_Term);
begin
Eval_Terms(p);
return res;
end Eval;
function Eval ( d : Degrees; c : number; x : Vector ) return number is
res : number;
begin
copy(c,res);
for i in d'range loop
for j in 1..(-d(i)) loop
Div(res,x(i));
end loop;
for j in 1..d(i) loop
Mul(res,x(i));
end loop;
end loop;
return res;
end Eval;
function Eval ( t : Term; x : Vector ) return number is
begin
return Eval(t.dg,t.cf,x);
end Eval;
function Eval ( t : Term; c : number; x : Vector ) return number is
begin
return Eval(t.dg,c,x);
end Eval;
function Eval ( p : Poly; x : Vector ) return number is
res : number;
procedure Eval_Term ( t : in Term; cont : out boolean ) is
tmp : number := Eval(t,x);
begin
Add(res,tmp);
Clear(tmp);
cont := true;
end Eval_Term;
procedure Eval_Terms is new Visiting_Iterator(Eval_Term);
begin
Copy(zero,res);
Eval_Terms(p);
return res;
end Eval;
function Eval ( p : Poly; c,x : Vector ) return number is
res : number;
cnt : natural := 0;
procedure Eval_Term ( t : in Term; cont : out boolean ) is
tmp : number := Eval(t,c(cnt),x);
begin
cnt := cnt + 1;
Add(res,tmp);
Clear(tmp);
cont := true;
end Eval_Term;
procedure Eval_Terms is new Visiting_Iterator(Eval_Term);
begin
Copy(zero,res);
Eval_Terms(p);
return res;
end Eval;
function Eval ( vp : Eval_Poly_Rep; x : Vector; i : integer ) return number;
function Eval ( vprec : Poly_Rec; x : Vector; i : integer ) return number is
res : number;
begin
if vprec.k = coefficient
then copy(vprec.c,res);
return res;
elsif vprec.p = null
then copy(zero,res);
return res;
else return Eval(vprec.p.all,x,i);
end if;
end Eval;
function Eval ( vp : Eval_Poly_Rep; x : Vector; i : integer ) return number is
deg : natural := vp'length-1;
res : number;
begin
if deg = 0
then return Eval(vp(0),x,i+1);
else res := Eval(vp(deg),x,i+1);
for j in reverse 0..(deg-1) loop
Mul(res,x(i));
if (vp(j).k = coefficient) or else (vp(j).p /= null)
then declare
temp : number := Eval(vp(j),x,i+1);
begin
Add(res,temp);
Clear(temp);
end;
end if;
end loop;
return res;
end if;
end Eval;
function Eval ( p : Eval_Poly; x : Vector ) return number is
begin
if p = null
then declare
res : number;
begin
Copy(zero,res);
return res;
end;
else return Eval(p.all,x,x'first);
end if;
end Eval;
function Eval ( vp : Eval_Coeff_Poly_Rep; c,x : Vector; i : integer )
return number;
function Eval ( vprec : Coeff_Poly_Rec; c,x : Vector; i : integer )
return number is
res : number;
begin
if vprec.k = coefficient
then copy(c(vprec.c),res);
return res;
elsif vprec.p = null
then copy(zero,res);
return res;
else return Eval(vprec.p.all,c,x,i);
end if;
end Eval;
function Eval ( vp : Eval_Coeff_Poly_Rep; c,x : Vector; i : integer )
return number is
deg : natural := vp'length-1;
res : number;
begin
if deg = 0
then return Eval(vp(0),c,x,i+1);
else res := Eval(vp(deg),c,x,i+1);
for j in reverse 0..(deg-1) loop
Mul(res,x(i));
if (vp(j).k = coefficient) or else (vp(j).p /= null)
then declare
temp : number := Eval(vp(j),c,x,i+1);
begin
Add(res,temp);
Clear(temp);
end;
end if;
end loop;
return res;
end if;
end Eval;
function Eval ( p : Eval_Coeff_Poly; c,x : Vector ) return number is
begin
if p = null
then declare
res : number;
begin
Copy(zero,res);
return res;
end;
else return Eval(p.all,c,x,x'first);
end if;
end Eval;
-- DESTRUCTORS :
procedure Clear ( p : in out Eval_Poly ) is
begin
if p /= null
then declare
vp : Eval_Poly_Rep renames p.all;
begin
for i in vp'range loop
if vp(i).k = coefficient
then Clear(vp(i).c);
else Clear(vp(i).p);
end if;
end loop;
end;
free(p);
end if;
end Clear;
procedure Clear ( p : in out Eval_Coeff_Poly ) is
begin
if p /= null
then declare
vp : Eval_Coeff_Poly_Rep renames p.all;
begin
for i in vp'range loop
if vp(i).k /= coefficient
then Clear(vp(i).p);
end if;
end loop;
end;
free(p);
end if;
end Clear;
end Generic_Laur_Poly_Functions;