update isl for move of interface/all.h
[barvinok.git] / series.cc
blob155c7eaaa1ed6cde74a14dc8ddefeae8fb840ede
1 #include <assert.h>
2 #include <NTL/vec_ZZ.h>
3 #include <barvinok/polylib.h>
4 #include <barvinok/barvinok.h>
5 #include <barvinok/util.h>
6 #include "genfun_constructor.h"
7 #include "lattice_width.h"
8 #include "remove_equalities.h"
10 using std::cerr;
11 using std::endl;
13 static gen_fun *enumerate_series(Polyhedron *P, unsigned nparam,
14 barvinok_options *options)
16 Matrix *CP = NULL;
17 gen_fun *gf;
18 Polyhedron *P_orig = P;
20 if (emptyQ2(P))
21 return new gen_fun(Empty_Polyhedron(nparam));
23 if (P->NbEq != 0)
24 remove_all_equalities(&P, NULL, &CP, NULL, nparam, options->MaxRays);
25 assert(emptyQ2(P) || P->NbEq == 0);
26 if (CP)
27 nparam = CP->NbColumns-1;
29 if (nparam == 0) {
30 Value c;
31 value_init(c);
32 barvinok_count_with_options(P, &c, options);
33 gf = new gen_fun(c);
34 value_clear(c);
35 } else {
36 POL_ENSURE_VERTICES(P);
37 if (P->NbEq)
38 gf = enumerate_series(P, nparam, options);
39 else {
40 gf_base *red;
41 red = gf_base::create(Polyhedron_Project(P, nparam),
42 P->Dimension, nparam, options);
43 red->start_gf(P, options);
44 gf = red->gf;
45 delete red;
48 if (CP) {
49 gf->substitute(CP);
50 Matrix_Free(CP);
52 if (P != P_orig)
53 Polyhedron_Free(P);
54 return gf;
57 gen_fun *barvinok_enumerate_series(Polyhedron *P, unsigned nparam,
58 barvinok_options *options)
60 if (emptyQ2(P))
61 return new gen_fun(Empty_Polyhedron(nparam));
63 assert(!Polyhedron_is_unbounded(P, nparam, options->MaxRays));
64 assert(P->NbBid == 0);
65 assert(Polyhedron_has_revlex_positive_rays(P, nparam));
66 return enumerate_series(P, nparam, options);
69 gen_fun * barvinok_series_with_options(Polyhedron *P, Polyhedron* C,
70 barvinok_options *options)
72 Polyhedron *CA;
73 unsigned nparam = C->Dimension;
74 gen_fun *gf;
76 CA = align_context(C, P->Dimension, options->MaxRays);
77 P = DomainIntersection(P, CA, options->MaxRays);
78 Polyhedron_Free(CA);
80 gf = barvinok_enumerate_series(P, nparam, options);
81 Polyhedron_Free(P);
83 return gf;
86 gen_fun * barvinok_series(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
88 gen_fun *gf;
89 barvinok_options *options = barvinok_options_new_with_defaults();
90 options->MaxRays = MaxRays;
91 gf = barvinok_series_with_options(P, C, options);
92 barvinok_options_free(options);
93 return gf;
96 gen_fun* barvinok_enumerate_union_series_with_options(Polyhedron *D, Polyhedron* C,
97 barvinok_options *options)
99 Polyhedron *CA;
100 gen_fun *gf = NULL, *gf2;
101 unsigned nparam = C->Dimension;
103 CA = align_context(C, D->Dimension, options->MaxRays);
104 D = DomainIntersection(D, CA, options->MaxRays);
105 Polyhedron_Free(CA);
107 for (Polyhedron *P = D; P; P = P->next) {
108 assert(P->Dimension == D->Dimension);
109 gen_fun *P_gf;
111 P_gf = barvinok_enumerate_series(P, P->Dimension, options);
112 if (!gf)
113 gf = P_gf;
114 else {
115 gf->add_union(P_gf, options);
116 delete P_gf;
119 /* we actually only need the convex union of the parameter space
120 * but the reducer classes currently expect a polyhedron in
121 * the combined space
123 Polyhedron_Free(gf->context);
124 gf->context = DomainConvex(D, options->MaxRays);
126 gf2 = gf->summate(D->Dimension - nparam, options);
128 delete gf;
129 Domain_Free(D);
130 return gf2;
133 gen_fun* barvinok_enumerate_union_series(Polyhedron *D, Polyhedron* C,
134 unsigned MaxRays)
136 gen_fun *gf;
137 barvinok_options *options = barvinok_options_new_with_defaults();
138 options->MaxRays = MaxRays;
139 gf = barvinok_enumerate_union_series_with_options(D, C, options);
140 barvinok_options_free(options);
141 return gf;
144 /* Unimodularly transform the polyhedron P, such that
145 * the direction specified by dir corresponds to the last
146 * variable in the transformed polyhedron.
147 * The number of variables is given by the length of dir.
149 static Polyhedron *put_direction_last(Polyhedron *P, Vector *dir,
150 unsigned MaxRays)
152 Polyhedron *R;
153 Matrix *T;
154 int n = dir->Size;
156 T = Matrix_Alloc(P->Dimension+1, P->Dimension+1);
157 T->NbColumns = T->NbRows = n;
158 Vector_Copy(dir->p, T->p[0], n);
159 unimodular_complete(T, 1);
160 Vector_Exchange(T->p[0], T->p[n-1], n);
161 T->NbColumns = T->NbRows = P->Dimension+1;
162 for (int j = n; j < P->Dimension+1; ++j)
163 value_set_si(T->p[j][j], 1);
165 R = Polyhedron_Image(P, T, MaxRays);
166 Matrix_Free(T);
168 return R;
171 /* Do we need to continue shifting and subtracting?
172 * i is the number of times we shifted so far
173 * n is the number of coordinates projected out
175 static bool more_shifts_needed(int j, int n,
176 gen_fun *S, gen_fun *S_divide, const vec_ZZ& up,
177 barvinok_options *options)
179 bool empty;
180 gen_fun *hp;
182 /* For the 2-dimensional case, we need to subtract at most once */
183 if (n == 2 && j > 0)
184 return false;
186 S_divide->shift(up);
188 /* Assume that we have to subtract at least once */
189 if (j == 0)
190 return true;
192 hp = S->Hadamard_product(S_divide, options);
194 empty = hp->is_zero();
195 delete hp;
197 return !empty;
200 static gen_fun *project(Polyhedron *P, unsigned n, barvinok_options *options,
201 int free_P);
203 /* Return gf of P projected on last dim(P)-n coordinates, i.e.,
204 * project out the first n coordinates.
206 * Assumes P has no equalities.
208 static gen_fun *project_full_dim(Polyhedron *P, unsigned n,
209 barvinok_options *options)
211 QQ one(1, 1);
212 QQ mone(-1, 1);
213 vec_ZZ up;
214 gen_fun *gf = NULL;
215 struct width_direction_array *dirs;
216 Polyhedron *U;
218 if (n == 0)
219 return barvinok_enumerate_series(P, P->Dimension, options);
221 up.SetLength(P->Dimension - (n-1));
222 up[0] = 1;
223 for (int i = 1; i < P->Dimension - (n-1); ++i)
224 up[i] = 0;
226 if (n == 1) {
227 gen_fun *S, *S_shift, *hp;
229 S = barvinok_enumerate_series(P, P->Dimension, options);
230 S_shift = new gen_fun(S);
231 S_shift->shift(up);
232 hp = S->Hadamard_product(S_shift, options);
233 S->add(mone, hp, options);
234 delete S_shift;
235 delete hp;
237 gf = S->summate(1, options);
238 delete S;
240 return gf;
243 U = Universe_Polyhedron(P->Dimension - n);
244 dirs = Polyhedron_Lattice_Width_Directions(P, U, options);
245 Polyhedron_Free(U);
247 for (int i = 0; i < dirs->n; ++i) {
248 Polyhedron *Pi, *R;
249 Polyhedron *CA;
250 gen_fun *S, *S_shift, *S_divide, *sum;
252 CA = align_context(dirs->wd[i].domain, P->Dimension, options->MaxRays);
253 R = DomainIntersection(P, CA, options->MaxRays);
254 Polyhedron_Free(CA);
255 assert(dirs->wd[i].dir->Size == n);
256 Pi = put_direction_last(R, dirs->wd[i].dir, options->MaxRays);
257 Polyhedron_Free(R);
259 S = project(Pi, n-1, options, 1);
261 S_shift = new gen_fun(S);
262 S_divide = new gen_fun(S);
263 S_divide->divide(up);
265 for (int j = 0; more_shifts_needed(j, n, S, S_divide, up, options); ++j) {
266 gen_fun *hp;
268 S_shift->shift(up);
269 hp = S->Hadamard_product(S_shift, options);
270 S->add(mone, hp, options);
272 delete hp;
275 sum = S->summate(1, options);
277 delete S_shift;
278 delete S_divide;
279 delete S;
281 if (!gf)
282 gf = sum;
283 else {
284 gf->add(one, sum, options);
285 delete sum;
288 free_width_direction_array(dirs);
290 return gf;
293 /* Return gf of P projected on last dim(P)-n coordinates, i.e.,
294 * project out the first n coordinates.
296 static gen_fun *project(Polyhedron *P, unsigned n, barvinok_options *options,
297 int free_P)
299 Matrix *CP = NULL;
300 gen_fun *proj;
301 unsigned nparam = P->Dimension - n;
303 if (P->NbEq != 0) {
304 Polyhedron *Q = P;
305 remove_all_equalities(&P, NULL, &CP, NULL, nparam, options->MaxRays);
306 if (CP)
307 nparam = CP->NbColumns - 1;
308 n = P->Dimension - nparam;
309 if (free_P)
310 Polyhedron_Free(Q);
311 free_P = 1;
314 if (emptyQ(P))
315 proj = new gen_fun(Empty_Polyhedron(nparam));
316 else
317 proj = project_full_dim(P, n, options);
318 if (CP) {
319 proj->substitute(CP);
320 Matrix_Free(CP);
323 if (free_P)
324 Polyhedron_Free(P);
326 return proj;
329 gen_fun *barvinok_enumerate_e_series(Polyhedron *P,
330 unsigned exist, unsigned nparam, barvinok_options *options)
332 Polyhedron *P_orig = P;
333 gen_fun *gf, *proj;
334 unsigned nvar = P->Dimension - exist - nparam;
336 if (exist == 0)
337 return barvinok_enumerate_series(P, nparam, options);
339 if (emptyQ(P))
340 return new gen_fun(Empty_Polyhedron(nparam));
342 assert(!Polyhedron_is_unbounded(P, nparam, options->MaxRays));
343 assert(P->NbBid == 0);
344 assert(Polyhedron_has_revlex_positive_rays(P, nparam));
346 /* Move existentially quantified variables to the front.*/
347 if (nvar) {
348 Matrix *T;
349 T = Matrix_Alloc(exist+nvar+nparam+1, nvar+exist+nparam+1);
350 for (int i = 0; i < exist; ++i)
351 value_set_si(T->p[i][nvar+i], 1);
352 for (int i = 0; i < nvar; ++i)
353 value_set_si(T->p[exist+i][i], 1);
354 for (int i = 0; i < nparam+1; ++i)
355 value_set_si(T->p[exist+nvar+i][nvar+exist+i], 1);
356 P = Polyhedron_Image(P, T, options->MaxRays);
357 Matrix_Free(T);
359 proj = project(P, exist, options, P != P_orig);
361 if (!nvar)
362 gf = proj;
363 else {
364 gf = proj->summate(nvar, options);
365 delete proj;
367 return gf;