File: [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Supports / face_enumerators_utilities.adb (download)
Revision 1.1.1.1 (vendor branch), Sun Oct 29 17:45:27 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_Common_Divisors; use Standard_Common_Divisors;
with Standard_Integer_Matrices; use Standard_Integer_Matrices;
with Standard_Integer_Linear_Solvers; use Standard_Integer_Linear_Solvers;
package body Face_Enumerators_Utilities is
function Is_Zero ( v : Vector ) return boolean is
begin
for i in v'range loop
if v(i) /= 0
then return false;
end if;
end loop;
return true;
end Is_Zero;
function gcd ( v : Vector ) return integer is
g : integer := v(v'first);
begin
if g < 0
then g := -g;
end if;
if g /= 1
then for i in v'first+1..v'last loop
if v(i) /= 0
then if v(i) < 0
then g := gcd(g,-v(i));
else g := gcd(g,v(i));
end if;
end if;
exit when (g = 1);
end loop;
end if;
return g;
end gcd;
procedure Scale ( v : in out Vector ) is
g : integer := gcd(v);
begin
if (g /= 0) and then (g /= 1)
then for i in v'range loop
v(i) := v(i)/g;
end loop;
end if;
end Scale;
function Is_In ( x : integer; v : Vector ) return boolean is
begin
for k in v'range loop
if v(k) = x
then return true;
end if;
end loop;
return false;
end Is_In;
function In_Edge ( x,a,b : Vector ) return boolean is
ba,xa,piv : integer;
begin
for i in b'range loop -- look for first i: b(i) - a(i) /= 0
ba := b(i) - a(i);
if ba /= 0
then piv := i; exit;
end if;
end loop;
if ba = 0
then return Equal(x,a); -- in this case b = a, so test x = a
else if ba < 0
then ba := -ba;
xa := x(piv) - b(piv);
else xa := x(piv) - a(piv);
end if;
if xa*ba >= 0 and then xa <= ba
then return true;
else return false;
end if;
end if;
end In_Edge;
function Last_Rows_Zero ( mat : matrix; frst : integer ) return boolean is
-- DESCRIPTION :
-- Returns true when the last rows of mat, starting at the index frst,
-- are completely zero.
function Zero_Row ( mat : matrix; i : integer ) return boolean is
begin
for j in mat'range(2) loop
if mat(i,j) /= 0
then return false;
end if;
end loop;
return true;
end Zero_Row;
begin
for i in frst..mat'last(1) loop
if not Zero_Row(mat,i)
then return false;
end if;
end loop;
return true;
end Last_Rows_Zero;
function In_Face ( k : in natural; face,x : Vector; pts : VecVec )
return boolean is
mat : Matrix(x'range,1..k+1);
sol : Vector(1..k+1) := (1..k+1 => 0);
sum : integer := 0;
begin
if k = 0
then return (pts(face(face'first)).all = x);
elsif k = 1
then return In_Edge(x,pts(face(face'first)).all,
pts(face(face'first+1)).all);
end if;
-- put_line("Calling In_Face for ");
-- put("x : "); put(x); new_line;
-- for i in face'range loop
-- put("pts("); put(face(i),1); put(") : "); put(pts(face(i))); new_line;
-- end loop;
-- compute the decomposition of x w.r.t. the face :
for j in 1..k loop
for i in mat'range(1) loop
mat(i,j) := pts(face(j))(i) - pts(face(0))(i);
end loop;
end loop;
for i in mat'range(1) loop
mat(i,k+1) := x(i) - pts(face(0))(i);
end loop;
-- put_line("The matrix which defines the decomposition :");
Upper_Triangulate(mat);
-- put(mat);
if not Last_Rows_Zero(mat,k+2)
then return false;
else declare
matk1 : matrix(1..k+1,1..k+1);
begin
for i in matk1'range(1) loop
for j in matk1'range(2) loop
matk1(i,j) := mat(i,j);
end loop;
end loop;
Scale(matk1); Solve0(matk1,sol);
-- put("The computed solution : "); put(sol); new_line;
for i in 1..k loop
if sol(i)*sol(k+1) > 0
then --put_line("x lies not in the face");
return false;
else sum := sum + sol(i);
end if;
end loop;
-- put("The sum : "); put(sum,1); new_line;
if sum = 0
then --put_line("x lies not in the face");
return false;
end if;
if sol(k+1) < 0
then sol(k+1) := -sol(k+1);
else sum := -sum;
end if;
if sum > sol(k+1)
then --put_line("x lies not in the face");
return false;
else --put_line("x lies in the face");
return true;
end if;
end;
end if;
end In_Face;
end Face_Enumerators_Utilities;