![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to face_enumerators_utilities.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Supports|
Return to face_enumerators_utilities.adb CVS log | Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Supports |
1.1 maekawa 1: with Standard_Common_Divisors; use Standard_Common_Divisors;
2: with Standard_Integer_Matrices; use Standard_Integer_Matrices;
3: with Standard_Integer_Linear_Solvers; use Standard_Integer_Linear_Solvers;
4:
5: package body Face_Enumerators_Utilities is
6:
7: function Is_Zero ( v : Vector ) return boolean is
8: begin
9: for i in v'range loop
10: if v(i) /= 0
11: then return false;
12: end if;
13: end loop;
14: return true;
15: end Is_Zero;
16:
17: function gcd ( v : Vector ) return integer is
18:
19: g : integer := v(v'first);
20:
21: begin
22: if g < 0
23: then g := -g;
24: end if;
25: if g /= 1
26: then for i in v'first+1..v'last loop
27: if v(i) /= 0
28: then if v(i) < 0
29: then g := gcd(g,-v(i));
30: else g := gcd(g,v(i));
31: end if;
32: end if;
33: exit when (g = 1);
34: end loop;
35: end if;
36: return g;
37: end gcd;
38:
39: procedure Scale ( v : in out Vector ) is
40:
41: g : integer := gcd(v);
42:
43: begin
44: if (g /= 0) and then (g /= 1)
45: then for i in v'range loop
46: v(i) := v(i)/g;
47: end loop;
48: end if;
49: end Scale;
50:
51: function Is_In ( x : integer; v : Vector ) return boolean is
52: begin
53: for k in v'range loop
54: if v(k) = x
55: then return true;
56: end if;
57: end loop;
58: return false;
59: end Is_In;
60:
61: function In_Edge ( x,a,b : Vector ) return boolean is
62:
63: ba,xa,piv : integer;
64:
65: begin
66: for i in b'range loop -- look for first i: b(i) - a(i) /= 0
67: ba := b(i) - a(i);
68: if ba /= 0
69: then piv := i; exit;
70: end if;
71: end loop;
72: if ba = 0
73: then return Equal(x,a); -- in this case b = a, so test x = a
74: else if ba < 0
75: then ba := -ba;
76: xa := x(piv) - b(piv);
77: else xa := x(piv) - a(piv);
78: end if;
79: if xa*ba >= 0 and then xa <= ba
80: then return true;
81: else return false;
82: end if;
83: end if;
84: end In_Edge;
85:
86: function Last_Rows_Zero ( mat : matrix; frst : integer ) return boolean is
87:
88: -- DESCRIPTION :
89: -- Returns true when the last rows of mat, starting at the index frst,
90: -- are completely zero.
91:
92: function Zero_Row ( mat : matrix; i : integer ) return boolean is
93: begin
94: for j in mat'range(2) loop
95: if mat(i,j) /= 0
96: then return false;
97: end if;
98: end loop;
99: return true;
100: end Zero_Row;
101:
102: begin
103: for i in frst..mat'last(1) loop
104: if not Zero_Row(mat,i)
105: then return false;
106: end if;
107: end loop;
108: return true;
109: end Last_Rows_Zero;
110:
111: function In_Face ( k : in natural; face,x : Vector; pts : VecVec )
112: return boolean is
113:
114: mat : Matrix(x'range,1..k+1);
115: sol : Vector(1..k+1) := (1..k+1 => 0);
116: sum : integer := 0;
117:
118: begin
119: if k = 0
120: then return (pts(face(face'first)).all = x);
121: elsif k = 1
122: then return In_Edge(x,pts(face(face'first)).all,
123: pts(face(face'first+1)).all);
124: end if;
125: -- put_line("Calling In_Face for ");
126: -- put("x : "); put(x); new_line;
127: -- for i in face'range loop
128: -- put("pts("); put(face(i),1); put(") : "); put(pts(face(i))); new_line;
129: -- end loop;
130: -- compute the decomposition of x w.r.t. the face :
131: for j in 1..k loop
132: for i in mat'range(1) loop
133: mat(i,j) := pts(face(j))(i) - pts(face(0))(i);
134: end loop;
135: end loop;
136: for i in mat'range(1) loop
137: mat(i,k+1) := x(i) - pts(face(0))(i);
138: end loop;
139: -- put_line("The matrix which defines the decomposition :");
140: Upper_Triangulate(mat);
141: -- put(mat);
142: if not Last_Rows_Zero(mat,k+2)
143: then return false;
144: else declare
145: matk1 : matrix(1..k+1,1..k+1);
146: begin
147: for i in matk1'range(1) loop
148: for j in matk1'range(2) loop
149: matk1(i,j) := mat(i,j);
150: end loop;
151: end loop;
152: Scale(matk1); Solve0(matk1,sol);
153: -- put("The computed solution : "); put(sol); new_line;
154: for i in 1..k loop
155: if sol(i)*sol(k+1) > 0
156: then --put_line("x lies not in the face");
157: return false;
158: else sum := sum + sol(i);
159: end if;
160: end loop;
161: -- put("The sum : "); put(sum,1); new_line;
162: if sum = 0
163: then --put_line("x lies not in the face");
164: return false;
165: end if;
166: if sol(k+1) < 0
167: then sol(k+1) := -sol(k+1);
168: else sum := -sum;
169: end if;
170: if sum > sol(k+1)
171: then --put_line("x lies not in the face");
172: return false;
173: else --put_line("x lies in the face");
174: return true;
175: end if;
176: end;
177: end if;
178: end In_Face;
179:
180: end Face_Enumerators_Utilities;