update polylib
[barvinok.git] / decomposer.cc
blob6b9344a94877e3b29553d724a66e0e1c82673222
1 #include <vector>
2 #include <assert.h>
3 #include <gmp.h>
4 #include <NTL/mat_ZZ.h>
5 #include <NTL/LLL.h>
6 #include <barvinok/barvinok.h>
7 #include <barvinok/util.h>
8 #include "conversion.h"
9 #include "decomposer.h"
10 #include "param_util.h"
11 #include "reduce_domain.h"
13 #ifdef NTL_STD_CXX
14 using namespace NTL;
15 #endif
16 using std::vector;
17 using std::cerr;
18 using std::endl;
21 * Returns the largest absolute value in the vector
23 static ZZ max(vec_ZZ& v)
25 ZZ max = abs(v[0]);
26 for (int i = 1; i < v.length(); ++i)
27 if (abs(v[i]) > max)
28 max = abs(v[i]);
29 return max;
32 /* Remove common divisor of elements of cols of B */
33 static void normalize_cols(mat_ZZ& B)
35 ZZ gcd;
36 for (int i = 0; i < B.NumCols(); ++i) {
37 gcd = B[0][i];
38 for (int j = 1 ; gcd != 1 && j < B.NumRows(); ++j)
39 GCD(gcd, gcd, B[j][i]);
40 if (gcd != 1)
41 for (int j = 0; j < B.NumRows(); ++j)
42 B[j][i] /= gcd;
46 /* Remove common divisor of elements of B */
47 static void normalize_matrix(mat_ZZ& B)
49 ZZ gcd;
50 for (int i = 0; i < B.NumCols(); ++i)
51 for (int j = 0 ; j < B.NumRows(); ++j) {
52 GCD(gcd, gcd, B[i][j]);
53 if (IsOne(gcd))
54 return;
56 for (int i = 0; i < B.NumCols(); ++i)
57 for (int j = 0; j < B.NumRows(); ++j)
58 B[i][j] /= gcd;
61 class cone {
62 public:
63 cone(const mat_ZZ& r, int row, const vec_ZZ& w, int s) {
64 sgn = s;
65 rays = r;
66 rays[row] = w;
67 set_det();
69 cone(const signed_cone& sc) {
70 rays = sc.rays;
71 sgn = sc.sign;
72 set_det();
74 void set_det() {
75 det = determinant(rays);
76 assert(!IsZero(det));
78 bool needs_split(barvinok_options *options) {
79 index = abs(det);
80 if (IsOne(index))
81 return false;
82 if (options->primal && index <= options->max_index)
83 return false;
85 inv(det, B, rays);
86 normalize_matrix(B);
87 if (sign(det) < 0)
88 negate(B, B);
90 if (!options->primal && options->max_index > 1) {
91 mat_ZZ B2 = B;
92 normalize_cols(B2);
93 index = abs(determinant(B2));
94 if (index <= options->max_index)
95 return false;
98 return true;
101 void short_vector(vec_ZZ& v, vec_ZZ& lambda, barvinok_options *options) {
102 ZZ det2;
103 mat_ZZ U;
104 long r = LLL(det2, B, U, options->LLL_a, options->LLL_b);
106 ZZ min = max(B[0]);
107 int index = 0;
108 for (int i = 1; i < B.NumRows(); ++i) {
109 ZZ tmp = max(B[i]);
110 if (tmp < min) {
111 min = tmp;
112 index = i;
116 lambda = B[index];
118 v = U[index];
120 int i;
121 for (i = 0; i < lambda.length(); ++i)
122 if (lambda[i] > 0)
123 break;
124 if (i == lambda.length()) {
125 v = -v;
126 lambda = -lambda;
130 ZZ det;
131 ZZ index;
132 mat_ZZ rays;
133 mat_ZZ B;
134 int sgn;
137 std::ostream & operator<<(std::ostream & os, const cone& c)
139 os << c.rays << endl;
140 os << "det: " << c.det << endl;
141 os << "sign: " << c.sgn << endl;
142 return os;
145 static void decompose(const signed_cone& sc, signed_cone_consumer& scc,
146 bool primal, barvinok_options *options)
148 vector<cone *> nonuni;
149 cone *c = new cone(sc);
150 if (c->needs_split(options)) {
151 nonuni.push_back(c);
152 } else {
153 try {
154 options->stats->base_cones++;
155 scc.handle(signed_cone(sc.C, sc.rays, sc.sign, to_ulong(c->index)),
156 options);
157 delete c;
158 } catch (...) {
159 delete c;
160 throw;
162 return;
164 vec_ZZ lambda;
165 vec_ZZ v;;
166 while (!nonuni.empty()) {
167 c = nonuni.back();
168 nonuni.pop_back();
169 c->short_vector(v, lambda, options);
170 for (int i = 0; i < c->rays.NumRows(); ++i) {
171 if (lambda[i] == 0)
172 continue;
173 cone *pc = new cone(c->rays, i, v, sign(lambda[i]) * c->sgn);
174 if (primal) {
175 for (int j = 0; j <= i; ++j) {
176 if ((j == i && sign(lambda[i]) < 0) ||
177 (j < i && sign(lambda[i]) == sign(lambda[j]))) {
178 pc->rays[j] = -pc->rays[j];
179 pc->sgn = -pc->sgn;
183 if (pc->needs_split(options)) {
184 assert(abs(pc->det) < abs(c->det));
185 nonuni.push_back(pc);
186 } else {
187 try {
188 options->stats->base_cones++;
189 scc.handle(signed_cone(pc->rays, pc->sgn,
190 to_ulong(pc->index)),
191 options);
192 delete pc;
193 } catch (...) {
194 delete c;
195 delete pc;
196 while (!nonuni.empty()) {
197 c = nonuni.back();
198 nonuni.pop_back();
199 delete c;
201 throw;
205 delete c;
209 struct polar_signed_cone_consumer : public signed_cone_consumer {
210 signed_cone_consumer& scc;
211 mat_ZZ r;
212 polar_signed_cone_consumer(signed_cone_consumer& scc) : scc(scc) {}
213 virtual void handle(const signed_cone& sc, barvinok_options *options) {
214 Polyhedron *C = sc.C;
215 if (!sc.C) {
216 Matrix *M = Matrix_Alloc(sc.rays.NumRows()+1, sc.rays.NumCols()+2);
217 for (int i = 0; i < sc.rays.NumRows(); ++i) {
218 zz2values(sc.rays[i], M->p[i]+1);
219 value_set_si(M->p[i][0], 1);
221 value_set_si(M->p[sc.rays.NumRows()][0], 1);
222 value_set_si(M->p[sc.rays.NumRows()][1+sc.rays.NumCols()], 1);
223 C = Rays2Polyhedron(M, M->NbRows+1);
224 assert(C->NbConstraints == C->NbRays);
225 Matrix_Free(M);
227 Polyhedron_Polarize(C);
228 rays(C, r);
229 try {
230 scc.handle(signed_cone(C, r, sc.sign, sc.det), options);
231 } catch (...) {
232 if (!sc.C)
233 Polyhedron_Free(C);
234 throw;
236 if (!sc.C)
237 Polyhedron_Free(C);
241 /* Remove common divisor of elements of rows of B */
242 static void normalize_rows(mat_ZZ& B)
244 ZZ gcd;
245 for (int i = 0; i < B.NumRows(); ++i) {
246 gcd = B[i][0];
247 for (int j = 1 ; gcd != 1 && j < B.NumCols(); ++j)
248 GCD(gcd, gcd, B[i][j]);
249 if (gcd != 1)
250 for (int j = 0; j < B.NumCols(); ++j)
251 B[i][j] /= gcd;
255 static void polar_decompose(Polyhedron *cone, signed_cone_consumer& scc,
256 barvinok_options *options)
258 POL_ENSURE_VERTICES(cone);
259 Polyhedron_Polarize(cone);
260 if (cone->NbRays - 1 != cone->Dimension) {
261 Polyhedron *tmp = cone;
262 cone = triangulate_cone_with_options(cone, options);
263 Polyhedron_Free(tmp);
265 polar_signed_cone_consumer pssc(scc);
266 mat_ZZ r;
267 try {
268 for (Polyhedron *Polar = cone; Polar; Polar = Polar->next) {
269 rays(Polar, r);
270 normalize_rows(r);
271 decompose(signed_cone(Polar, r, 1), pssc, false, options);
273 Domain_Free(cone);
274 } catch (...) {
275 Domain_Free(cone);
276 throw;
280 static void primal_decompose(Polyhedron *cone, signed_cone_consumer& scc,
281 barvinok_options *options)
283 POL_ENSURE_VERTICES(cone);
284 Polyhedron *parts;
285 if (cone->NbRays - 1 == cone->Dimension)
286 parts = cone;
287 else
288 parts = triangulate_cone_with_options(cone, options);
289 Vector *average = NULL;
290 Value tmp;
291 if (parts != cone) {
292 value_init(tmp);
293 average = inner_point(cone);
295 mat_ZZ ray;
296 try {
297 for (Polyhedron *simple = parts; simple; simple = simple->next) {
298 int sign = 1;
299 Matrix *Rays = rays2(simple);
300 for (int i = 0; i < Rays->NbRows; ++i) {
301 if (simple == cone) {
302 continue;
303 } else {
304 int f;
305 for (f = 0; f < simple->NbConstraints; ++f) {
306 Inner_Product(Rays->p[i], simple->Constraint[f]+1,
307 simple->Dimension, &tmp);
308 if (value_notzero_p(tmp))
309 break;
311 assert(f < simple->NbConstraints);
312 if (!is_internal(average, simple->Constraint[f])) {
313 Vector_Oppose(Rays->p[i], Rays->p[i], Rays->NbColumns);
314 sign = -sign;
318 matrix2zz(Rays, ray, Rays->NbRows, Rays->NbColumns);
319 Matrix_Free(Rays);
320 decompose(signed_cone(simple, ray, sign), scc, true, options);
322 Domain_Free(parts);
323 if (parts != cone) {
324 Domain_Free(cone);
325 value_clear(tmp);
326 Vector_Free(average);
328 } catch (...) {
329 Domain_Free(parts);
330 if (parts != cone) {
331 Domain_Free(cone);
332 value_clear(tmp);
333 Vector_Free(average);
335 throw;
339 void barvinok_decompose(Polyhedron *C, signed_cone_consumer& scc,
340 barvinok_options *options)
342 POL_ENSURE_VERTICES(C);
343 if (options->primal)
344 primal_decompose(C, scc, options);
345 else
346 polar_decompose(C, scc, options);
349 void vertex_decomposer::decompose_at_vertex(Param_Vertices *V, int _i,
350 barvinok_options *options)
352 Polyhedron *C = Param_Vertex_Cone(PP, V, options);
353 vert = _i;
354 this->V = V;
356 barvinok_decompose(C, scc, options);
359 struct posneg_collector : public signed_cone_consumer {
360 posneg_collector(Polyhedron *pos, Polyhedron *neg) : pos(pos), neg(neg) {}
361 virtual void handle(const signed_cone& sc, barvinok_options *options) {
362 Polyhedron *p = Polyhedron_Copy(sc.C);
363 if (sc.sign > 0) {
364 p->next = pos;
365 pos = p;
366 } else {
367 p->next = neg;
368 neg = p;
371 Polyhedron *pos;
372 Polyhedron *neg;
376 * Barvinok's Decomposition of a simplicial cone
378 * Returns two lists of polyhedra
380 void barvinok_decompose(Polyhedron *C, Polyhedron **ppos, Polyhedron **pneg)
382 barvinok_options *options = barvinok_options_new_with_defaults();
383 posneg_collector pc(*ppos, *pneg);
384 POL_ENSURE_VERTICES(C);
385 mat_ZZ r;
386 rays(C, r);
387 decompose(signed_cone(C, r, 1), pc, false, options);
388 *ppos = pc.pos;
389 *pneg = pc.neg;
390 barvinok_options_free(options);