with integer_io; use integer_io;
with Communications_with_User; use Communications_with_User;
with Standard_Floating_Numbers; use Standard_Floating_Numbers;
with Standard_Complex_Numbers; use Standard_Complex_Numbers;
with Standard_Complex_Numbers_io; use Standard_Complex_Numbers_io;
with Numbers_io; use Numbers_io;
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 Interpolating_Homotopies; use Interpolating_Homotopies;
procedure Interpolating_Homotopies_Driver
( file : in file_type; p : in Poly_Sys; z : in Partition;
b : in out natural; q : out Poly_Sys;
qsols : in out Solution_List ) is
n : constant natural := p'last;
interpols : Solution_List;
ib,scalind : natural;
ans : character;
function Random_Interpolating ( n,m : natural ) return Solution_List is
-- DESCRIPTION :
-- A list of m random n-dimensional vectors will be returned.
-- The complex numbers will all have modulus one.
res,res_last : Solution_List;
s : Solution(n);
begin
s.t := Create(0.0);
s.m := 1;
s.err := 0.0; s.rco := 1.0; s.res := 0.0;
for i in 1..m loop
for j in 1..n loop
s.v(j) := random1;
end loop;
Append(res,res_last,s);
end loop;
return res;
end Random_Interpolating;
procedure Random_Linear_Scaler ( n : in natural; p : in Poly;
v : out vector; l : out natural ) is
-- DESCRIPTION :
-- Returns a random vector of dimension n+1, with range 0..n.
-- There will be a nonzero entry only for those unknowns that occur in p.
res : Vector(0..n);
last : natural := 0;
begin
for i in res'range loop
if Degree(p,i) > 0
then res(i) := random1;
last := last + 1;
else res(i) := Create(0.0);
end if;
end loop;
v := res; l := last;
end Random_Linear_Scaler;
function Scale ( sc,v : Vector; last : integer ) return Complex_Number is
-- DESCRIPTION :
-- Returns the last component of the vector v, that is v(last),
-- such that sc(0) + sum sc(i)*v(i), i in v'range, holds.
res : Complex_Number := sc(0);
begin
for i in v'first..last-1 loop
res := res + sc(i)*v(i);
end loop;
res := -res/sc(last);
return res;
end Scale;
function Random_Interpolating
( n,m : natural; scaler : Vector; scallast : natural )
return Solution_List is
-- DESCRIPTION :
-- A list of m random n-dimensional vectors will be returned.
-- The complex numbers will all have modulus one, except the last one,
-- indicated by scallast, that has been chosen to satisfy the scaler
-- equation, defined by sum of scaler(i)*x_i = 0, with x_0 = 1.
res,res_last : Solution_List;
s : Solution(n);
begin
s.t := Create(0.0);
s.m := 1;
for i in 1..m loop
for j in 1..(n-1) loop
s.v(j) := random1;
end loop;
s.v(n) := Scale(scaler,s.v,scallast);
Append(res,res_last,s);
end loop;
return res;
end Random_Interpolating;
function Create ( v : Vector ) return Poly is
res : Poly := Null_Poly;
t : Term;
begin
t.dg := new Standard_Natural_Vectors.Vector'(v'first+1..v'last => 0);
for i in v'range loop
t.cf := v(i);
if i > v'first
then t.dg(i) := 1;
end if;
Add(res,t);
if i > v'first
then t.dg(i) := 0;
end if;
end loop;
Clear(t);
return res;
end Create;
function Interpolating_by_User ( n,m : natural ) return Solution_List is
-- DESCRIPTION :
-- A list of m n-dimensional vectors will be read from standard input.
res,res_last : Solution_List;
s : Solution(n);
-- f1,f2 : double_float;
begin
put("Reading "); put(m,1); put(" "); put(n,1);
put_line("-dimensional complex vectors.");
for i in 1..m loop
s.t := Create(0.0);
s.m := 1;
s.err := 0.0; s.rco := 1.0; s.res := 0.0;
put("Give the components of vector "); put(i,1);
put_line(" :");
for j in 1..n loop
-- Read_Double_Float(f1);
-- Read_Double_Float(f2);
-- s.v(j) := Create(f1,f2);
get(s.v(j));
end loop;
Append(res,res_last,s);
end loop;
return res;
end Interpolating_by_User;
procedure Driver_for_Interpolation is
-- DESCRIPTION : interpolation without a scaling equation
dp,ip : Poly_Sys(p'range);
begin
dp := Dense_Representation(p,z);
ip := Independent_Representation(dp);
ib := Independent_Roots(ip);
if ib > b
then ib := b;
end if;
put("The number of independent roots : "); put(ib,1); new_line;
put("Do you want to give interpolation vectors by yourself ? (y/n) ");
Ask_Yes_or_No(ans);
if ans = 'y'
then interpols := Interpolating_by_User(n,ib);
else put_line("Random interpolating vectors will be generated.");
interpols := Random_Interpolating(n,ib);
end if;
q := Interpolate(ip,ib,interpols);
qsols := interpols; b := ib;
Clear(dp); Clear(ip);
end Driver_for_Interpolation;
procedure Driver_for_Scaled_Interpolation is
-- DESCRIPTION : interpolation with a scaling equation, p(scalind).
dp,ip,pp,qq : Poly_Sys(p'range);
scalvec : Vector(0..n);
scalveclast : natural;
begin
Random_Linear_Scaler(n,p(scalind),scalvec,scalveclast);
for i in p'range loop
if i = scalind
then pp(i) := Null_Poly;
else pp(i) := p(i);
end if;
end loop;
dp := Dense_Representation(pp,z);
ip := Independent_Representation(dp);
ib := Independent_Roots(ip,scalveclast);
if ib > b
then ib := b;
end if;
put("The number of independent roots : "); put(ib,1); new_line;
put("Do you want to give interpolation vectors by yourself ? (y/n) ");
Ask_Yes_or_No(ans);
if ans = 'y'
then interpols := Interpolating_by_User(n,ib);
else put_line("Random interpolating vectors will be generated.");
interpols := Random_Interpolating(n,ib,scalvec,scalveclast);
end if;
qq := Interpolate(ip,scalveclast,ib,interpols);
qq(scalind) := Create(scalvec);
qsols := interpols; b := ib; q := qq;
Clear(dp); Clear(ip);
end Driver_for_Scaled_Interpolation;
begin
new_line;
put("Give the number of the linear scaling equation (0 if none) : ");
Read_Natural(scalind);
if scalind = 0
then Driver_for_Interpolation;
else Driver_for_Scaled_Interpolation;
end if;
end Interpolating_Homotopies_Driver;