remove_parameter_equalities: stop "removing" equalities if polyhedron is empty
[barvinok.git] / series.cc
blobff448f1ac0bd96db7f1efd7bcfdfe6ce73038b60
1 #include <assert.h>
2 #include <barvinok/barvinok.h>
3 #include <barvinok/util.h>
4 #include "genfun_constructor.h"
5 #include "lattice_width.h"
6 #include "remove_equalities.h"
8 using std::cerr;
9 using std::endl;
11 /* Check whether all rays point in the positive directions
12 * for the parameters
14 static bool Polyhedron_has_positive_rays(Polyhedron *P, unsigned nparam)
16 int r;
17 for (r = 0; r < P->NbRays; ++r)
18 if (value_zero_p(P->Ray[r][P->Dimension+1])) {
19 int i;
20 for (i = P->Dimension - nparam; i < P->Dimension; ++i)
21 if (value_neg_p(P->Ray[r][i+1]))
22 return false;
24 return true;
27 static gen_fun *series(Polyhedron *P, unsigned nparam, barvinok_options *options)
29 Matrix *CP = NULL;
30 gen_fun *gf;
31 Polyhedron *P_orig = P;
33 if (emptyQ2(P))
34 return new gen_fun(Empty_Polyhedron(nparam));
36 assert(!Polyhedron_is_unbounded(P, nparam, options->MaxRays));
37 assert(P->NbBid == 0);
38 assert(Polyhedron_has_revlex_positive_rays(P, nparam));
39 if (P->NbEq != 0)
40 remove_all_equalities(&P, NULL, &CP, NULL, nparam, options->MaxRays);
41 assert(emptyQ2(P) || P->NbEq == 0);
42 if (CP)
43 nparam = CP->NbColumns-1;
45 if (nparam == 0) {
46 Value c;
47 value_init(c);
48 barvinok_count_with_options(P, &c, options);
49 gf = new gen_fun(c);
50 value_clear(c);
51 } else {
52 POL_ENSURE_VERTICES(P);
53 if (P->NbEq)
54 gf = series(P, nparam, options);
55 else {
56 gf_base *red;
57 red = gf_base::create(Polyhedron_Project(P, nparam),
58 P->Dimension, nparam, options);
59 red->start_gf(P, options);
60 gf = red->gf;
61 delete red;
64 if (CP) {
65 gf->substitute(CP);
66 Matrix_Free(CP);
68 if (P != P_orig)
69 Polyhedron_Free(P);
70 return gf;
73 gen_fun * barvinok_series_with_options(Polyhedron *P, Polyhedron* C,
74 barvinok_options *options)
76 Polyhedron *CA;
77 unsigned nparam = C->Dimension;
78 gen_fun *gf;
80 CA = align_context(C, P->Dimension, options->MaxRays);
81 P = DomainIntersection(P, CA, options->MaxRays);
82 Polyhedron_Free(CA);
84 gf = series(P, nparam, options);
85 Polyhedron_Free(P);
87 return gf;
90 gen_fun * barvinok_series(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
92 gen_fun *gf;
93 barvinok_options *options = barvinok_options_new_with_defaults();
94 options->MaxRays = MaxRays;
95 gf = barvinok_series_with_options(P, C, options);
96 barvinok_options_free(options);
97 return gf;
100 static Polyhedron *skew_into_positive_orthant(Polyhedron *D, unsigned nparam,
101 unsigned MaxRays)
103 Matrix *M = NULL;
104 Value tmp;
105 value_init(tmp);
106 for (Polyhedron *P = D; P; P = P->next) {
107 POL_ENSURE_VERTICES(P);
108 assert(!Polyhedron_is_unbounded(P, nparam, MaxRays));
109 assert(P->NbBid == 0);
110 assert(Polyhedron_has_positive_rays(P, nparam));
112 for (int r = 0; r < P->NbRays; ++r) {
113 if (value_notzero_p(P->Ray[r][P->Dimension+1]))
114 continue;
115 for (int i = 0; i < nparam; ++i) {
116 int j;
117 if (value_posz_p(P->Ray[r][i+1]))
118 continue;
119 if (!M) {
120 M = Matrix_Alloc(D->Dimension+1, D->Dimension+1);
121 for (int i = 0; i < D->Dimension+1; ++i)
122 value_set_si(M->p[i][i], 1);
123 } else {
124 Inner_Product(P->Ray[r]+1, M->p[i], D->Dimension+1, &tmp);
125 if (value_posz_p(tmp))
126 continue;
128 for (j = P->Dimension - nparam; j < P->Dimension; ++j)
129 if (value_pos_p(P->Ray[r][j+1]))
130 break;
131 assert(j < P->Dimension);
132 value_pdivision(tmp, P->Ray[r][j+1], P->Ray[r][i+1]);
133 value_subtract(M->p[i][j], M->p[i][j], tmp);
137 value_clear(tmp);
138 if (M) {
139 D = DomainImage(D, M, MaxRays);
140 Matrix_Free(M);
142 return D;
145 gen_fun* barvinok_enumerate_union_series_with_options(Polyhedron *D, Polyhedron* C,
146 barvinok_options *options)
148 Polyhedron *conv, *D2;
149 Polyhedron *CA;
150 gen_fun *gf = NULL, *gf2;
151 unsigned nparam = C->Dimension;
152 ZZ one, mone;
153 one = 1;
154 mone = -1;
156 CA = align_context(C, D->Dimension, options->MaxRays);
157 D = DomainIntersection(D, CA, options->MaxRays);
158 Polyhedron_Free(CA);
160 D2 = skew_into_positive_orthant(D, nparam, options->MaxRays);
161 for (Polyhedron *P = D2; P; P = P->next) {
162 assert(P->Dimension == D2->Dimension);
163 gen_fun *P_gf;
165 P_gf = series(P, P->Dimension, options);
166 if (!gf)
167 gf = P_gf;
168 else {
169 gf->add_union(P_gf, options);
170 delete P_gf;
173 /* we actually only need the convex union of the parameter space
174 * but the reducer classes currently expect a polyhedron in
175 * the combined space
177 Polyhedron_Free(gf->context);
178 gf->context = DomainConvex(D2, options->MaxRays);
180 gf2 = gf->summate(D2->Dimension - nparam, options);
182 delete gf;
183 if (D != D2)
184 Domain_Free(D2);
185 Domain_Free(D);
186 return gf2;
189 gen_fun* barvinok_enumerate_union_series(Polyhedron *D, Polyhedron* C,
190 unsigned MaxRays)
192 gen_fun *gf;
193 barvinok_options *options = barvinok_options_new_with_defaults();
194 options->MaxRays = MaxRays;
195 gf = barvinok_enumerate_union_series_with_options(D, C, options);
196 barvinok_options_free(options);
197 return gf;
200 /* Unimodularly transform the polyhedron P, such that
201 * the direction specified by dir corresponds to the last
202 * variable in the transformed polyhedron.
203 * The number of variables is given by the length of dir.
205 static Polyhedron *put_direction_last(Polyhedron *P, Vector *dir,
206 unsigned MaxRays)
208 Polyhedron *R;
209 Matrix *T;
210 int n = dir->Size;
212 T = Matrix_Alloc(P->Dimension+1, P->Dimension+1);
213 T->NbColumns = T->NbRows = n;
214 Vector_Copy(dir->p, T->p[0], n);
215 unimodular_complete(T, 1);
216 Vector_Exchange(T->p[0], T->p[n-1], n);
217 T->NbColumns = T->NbRows = P->Dimension+1;
218 for (int j = n; j < P->Dimension+1; ++j)
219 value_set_si(T->p[j][j], 1);
221 R = Polyhedron_Image(P, T, MaxRays);
222 Matrix_Free(T);
224 return R;
227 /* Do we need to continue shifting and subtracting?
228 * i is the number of times we shifted so far
229 * n is the number of coordinates projected out
231 static bool more_shifts_needed(int j, int n,
232 gen_fun *S, gen_fun *S_divide, const vec_ZZ& up,
233 barvinok_options *options)
235 bool empty;
236 gen_fun *hp;
237 Value c;
239 /* For the 2-dimensional case, we need to subtract at most once */
240 if (n == 2 && j > 0)
241 return false;
243 S_divide->shift(up);
245 /* Assume that we have to subtract at least once */
246 if (j == 0)
247 return true;
249 hp = S->Hadamard_product(S_divide, options);
251 value_init(c);
253 empty = hp->summate(&c) && value_zero_p(c);
254 delete hp;
256 value_clear(c);
258 return !empty;
261 /* Return gf of P projected on last dim(P)-n coordinates, i.e.,
262 * project out the first n coordinates.
264 gen_fun *project(Polyhedron *P, unsigned n, barvinok_options *options)
266 QQ one(1, 1);
267 QQ mone(-1, 1);
268 vec_ZZ up;
269 gen_fun *gf = NULL;
270 struct width_direction_array *dirs;
271 Polyhedron *U;
273 up.SetLength(P->Dimension - (n-1));
274 up[0] = 1;
275 for (int i = 1; i < P->Dimension - (n-1); ++i)
276 up[i] = 0;
278 if (n == 1) {
279 gen_fun *S, *S_shift, *hp;
281 S = series(P, P->Dimension, options);
282 S_shift = new gen_fun(S);
283 S_shift->shift(up);
284 hp = S->Hadamard_product(S_shift, options);
285 S->add(mone, hp, options);
286 delete S_shift;
287 delete hp;
289 gf = S->summate(1, options);
290 delete S;
292 return gf;
295 U = Universe_Polyhedron(P->Dimension - n);
296 dirs = Polyhedron_Lattice_Width_Directions(P, U, options);
297 Polyhedron_Free(U);
299 for (int i = 0; i < dirs->n; ++i) {
300 Polyhedron *Pi, *R;
301 Polyhedron *CA;
302 gen_fun *S, *S_shift, *S_divide, *sum;
304 CA = align_context(dirs->wd[i].domain, P->Dimension, options->MaxRays);
305 R = DomainIntersection(P, CA, options->MaxRays);
306 Polyhedron_Free(CA);
307 assert(dirs->wd[i].dir->Size == n);
308 Pi = put_direction_last(R, dirs->wd[i].dir, options->MaxRays);
309 Polyhedron_Free(R);
311 S = project(Pi, n-1, options);
312 Polyhedron_Free(Pi);
314 S_shift = new gen_fun(S);
315 S_divide = new gen_fun(S);
316 S_divide->divide(up);
318 for (int j = 0; more_shifts_needed(j, n, S, S_divide, up, options); ++j) {
319 gen_fun *hp;
321 S_shift->shift(up);
322 hp = S->Hadamard_product(S_shift, options);
323 S->add(mone, hp, options);
325 delete hp;
328 sum = S->summate(1, options);
330 delete S_shift;
331 delete S_divide;
332 delete S;
334 if (!gf)
335 gf = sum;
336 else {
337 gf->add(one, sum, options);
338 delete sum;
341 free_width_direction_array(dirs);
343 return gf;
346 gen_fun *barvinok_enumerate_e_series(Polyhedron *P,
347 unsigned exist, unsigned nparam, barvinok_options *options)
349 Matrix *CP = NULL;
350 Polyhedron *P_orig = P;
351 gen_fun *gf, *proj;
352 unsigned nvar = P->Dimension - exist - nparam;
354 if (exist == 0)
355 return series(P, nparam, options);
357 if (emptyQ(P))
358 return new gen_fun(Empty_Polyhedron(nparam));
360 assert(!Polyhedron_is_unbounded(P, nparam, options->MaxRays));
361 assert(P->NbBid == 0);
362 assert(Polyhedron_has_revlex_positive_rays(P, nparam));
364 /* Move existentially quantified variables to the front.*/
365 if (nvar) {
366 Matrix *T;
367 T = Matrix_Alloc(exist+nvar+nparam+1, nvar+exist+nparam+1);
368 for (int i = 0; i < exist; ++i)
369 value_set_si(T->p[i][nvar+i], 1);
370 for (int i = 0; i < nvar; ++i)
371 value_set_si(T->p[exist+i][i], 1);
372 for (int i = 0; i < nparam+1; ++i)
373 value_set_si(T->p[exist+nvar+i][nvar+exist+i], 1);
374 P = Polyhedron_Image(P, T, options->MaxRays);
375 Matrix_Free(T);
377 if (P->NbEq != 0) {
378 Polyhedron *Q = P;
379 remove_all_equalities(&P, NULL, &CP, NULL, nvar+nparam,
380 options->MaxRays);
381 exist = P->Dimension - (CP->NbColumns-1);
382 if (Q != P_orig)
383 Polyhedron_Free(Q);
386 proj = project(P, exist, options);
387 if (CP) {
388 proj->substitute(CP);
389 Matrix_Free(CP);
392 if (P != P_orig)
393 Polyhedron_Free(P);
395 if (!nvar)
396 gf = proj;
397 else {
398 gf = proj->summate(nvar, options);
399 delete proj;
401 return gf;