doc: add another paper refering to the library
[barvinok.git] / series.cc
blob4ebf160e368f04d92076c919586b1629b848b593
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 *enumerate_series(Polyhedron *P, unsigned nparam,
28 barvinok_options *options)
30 Matrix *CP = NULL;
31 gen_fun *gf;
32 Polyhedron *P_orig = P;
34 if (emptyQ2(P))
35 return new gen_fun(Empty_Polyhedron(nparam));
37 if (P->NbEq != 0)
38 remove_all_equalities(&P, NULL, &CP, NULL, nparam, options->MaxRays);
39 assert(emptyQ2(P) || P->NbEq == 0);
40 if (CP)
41 nparam = CP->NbColumns-1;
43 if (nparam == 0) {
44 Value c;
45 value_init(c);
46 barvinok_count_with_options(P, &c, options);
47 gf = new gen_fun(c);
48 value_clear(c);
49 } else {
50 POL_ENSURE_VERTICES(P);
51 if (P->NbEq)
52 gf = enumerate_series(P, nparam, options);
53 else {
54 gf_base *red;
55 red = gf_base::create(Polyhedron_Project(P, nparam),
56 P->Dimension, nparam, options);
57 red->start_gf(P, options);
58 gf = red->gf;
59 delete red;
62 if (CP) {
63 gf->substitute(CP);
64 Matrix_Free(CP);
66 if (P != P_orig)
67 Polyhedron_Free(P);
68 return gf;
71 gen_fun *barvinok_enumerate_series(Polyhedron *P, unsigned nparam,
72 barvinok_options *options)
74 if (emptyQ2(P))
75 return new gen_fun(Empty_Polyhedron(nparam));
77 assert(!Polyhedron_is_unbounded(P, nparam, options->MaxRays));
78 assert(P->NbBid == 0);
79 assert(Polyhedron_has_revlex_positive_rays(P, nparam));
80 return enumerate_series(P, nparam, options);
83 gen_fun * barvinok_series_with_options(Polyhedron *P, Polyhedron* C,
84 barvinok_options *options)
86 Polyhedron *CA;
87 unsigned nparam = C->Dimension;
88 gen_fun *gf;
90 CA = align_context(C, P->Dimension, options->MaxRays);
91 P = DomainIntersection(P, CA, options->MaxRays);
92 Polyhedron_Free(CA);
94 gf = barvinok_enumerate_series(P, nparam, options);
95 Polyhedron_Free(P);
97 return gf;
100 gen_fun * barvinok_series(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
102 gen_fun *gf;
103 barvinok_options *options = barvinok_options_new_with_defaults();
104 options->MaxRays = MaxRays;
105 gf = barvinok_series_with_options(P, C, options);
106 barvinok_options_free(options);
107 return gf;
110 static Polyhedron *skew_into_positive_orthant(Polyhedron *D, unsigned nparam,
111 unsigned MaxRays)
113 Matrix *M = NULL;
114 Value tmp;
115 value_init(tmp);
116 for (Polyhedron *P = D; P; P = P->next) {
117 POL_ENSURE_VERTICES(P);
118 assert(!Polyhedron_is_unbounded(P, nparam, MaxRays));
119 assert(P->NbBid == 0);
120 assert(Polyhedron_has_positive_rays(P, nparam));
122 for (int r = 0; r < P->NbRays; ++r) {
123 if (value_notzero_p(P->Ray[r][P->Dimension+1]))
124 continue;
125 for (int i = 0; i < nparam; ++i) {
126 int j;
127 if (value_posz_p(P->Ray[r][i+1]))
128 continue;
129 if (!M) {
130 M = Matrix_Alloc(D->Dimension+1, D->Dimension+1);
131 for (int i = 0; i < D->Dimension+1; ++i)
132 value_set_si(M->p[i][i], 1);
133 } else {
134 Inner_Product(P->Ray[r]+1, M->p[i], D->Dimension+1, &tmp);
135 if (value_posz_p(tmp))
136 continue;
138 for (j = P->Dimension - nparam; j < P->Dimension; ++j)
139 if (value_pos_p(P->Ray[r][j+1]))
140 break;
141 assert(j < P->Dimension);
142 value_pdivision(tmp, P->Ray[r][j+1], P->Ray[r][i+1]);
143 value_subtract(M->p[i][j], M->p[i][j], tmp);
147 value_clear(tmp);
148 if (M) {
149 D = DomainImage(D, M, MaxRays);
150 Matrix_Free(M);
152 return D;
155 gen_fun* barvinok_enumerate_union_series_with_options(Polyhedron *D, Polyhedron* C,
156 barvinok_options *options)
158 Polyhedron *conv, *D2;
159 Polyhedron *CA;
160 gen_fun *gf = NULL, *gf2;
161 unsigned nparam = C->Dimension;
162 ZZ one, mone;
163 one = 1;
164 mone = -1;
166 CA = align_context(C, D->Dimension, options->MaxRays);
167 D = DomainIntersection(D, CA, options->MaxRays);
168 Polyhedron_Free(CA);
170 D2 = skew_into_positive_orthant(D, nparam, options->MaxRays);
171 for (Polyhedron *P = D2; P; P = P->next) {
172 assert(P->Dimension == D2->Dimension);
173 gen_fun *P_gf;
175 P_gf = barvinok_enumerate_series(P, P->Dimension, options);
176 if (!gf)
177 gf = P_gf;
178 else {
179 gf->add_union(P_gf, options);
180 delete P_gf;
183 /* we actually only need the convex union of the parameter space
184 * but the reducer classes currently expect a polyhedron in
185 * the combined space
187 Polyhedron_Free(gf->context);
188 gf->context = DomainConvex(D2, options->MaxRays);
190 gf2 = gf->summate(D2->Dimension - nparam, options);
192 delete gf;
193 if (D != D2)
194 Domain_Free(D2);
195 Domain_Free(D);
196 return gf2;
199 gen_fun* barvinok_enumerate_union_series(Polyhedron *D, Polyhedron* C,
200 unsigned MaxRays)
202 gen_fun *gf;
203 barvinok_options *options = barvinok_options_new_with_defaults();
204 options->MaxRays = MaxRays;
205 gf = barvinok_enumerate_union_series_with_options(D, C, options);
206 barvinok_options_free(options);
207 return gf;
210 /* Unimodularly transform the polyhedron P, such that
211 * the direction specified by dir corresponds to the last
212 * variable in the transformed polyhedron.
213 * The number of variables is given by the length of dir.
215 static Polyhedron *put_direction_last(Polyhedron *P, Vector *dir,
216 unsigned MaxRays)
218 Polyhedron *R;
219 Matrix *T;
220 int n = dir->Size;
222 T = Matrix_Alloc(P->Dimension+1, P->Dimension+1);
223 T->NbColumns = T->NbRows = n;
224 Vector_Copy(dir->p, T->p[0], n);
225 unimodular_complete(T, 1);
226 Vector_Exchange(T->p[0], T->p[n-1], n);
227 T->NbColumns = T->NbRows = P->Dimension+1;
228 for (int j = n; j < P->Dimension+1; ++j)
229 value_set_si(T->p[j][j], 1);
231 R = Polyhedron_Image(P, T, MaxRays);
232 Matrix_Free(T);
234 return R;
237 /* Do we need to continue shifting and subtracting?
238 * i is the number of times we shifted so far
239 * n is the number of coordinates projected out
241 static bool more_shifts_needed(int j, int n,
242 gen_fun *S, gen_fun *S_divide, const vec_ZZ& up,
243 barvinok_options *options)
245 bool empty;
246 gen_fun *hp;
248 /* For the 2-dimensional case, we need to subtract at most once */
249 if (n == 2 && j > 0)
250 return false;
252 S_divide->shift(up);
254 /* Assume that we have to subtract at least once */
255 if (j == 0)
256 return true;
258 hp = S->Hadamard_product(S_divide, options);
260 empty = hp->is_zero();
261 delete hp;
263 return !empty;
266 static gen_fun *project(Polyhedron *P, unsigned n, barvinok_options *options,
267 int free_P);
269 /* Return gf of P projected on last dim(P)-n coordinates, i.e.,
270 * project out the first n coordinates.
272 * Assumes P has no equalities.
274 static gen_fun *project_full_dim(Polyhedron *P, unsigned n,
275 barvinok_options *options)
277 QQ one(1, 1);
278 QQ mone(-1, 1);
279 vec_ZZ up;
280 gen_fun *gf = NULL;
281 struct width_direction_array *dirs;
282 Polyhedron *U;
284 if (n == 0)
285 return barvinok_enumerate_series(P, P->Dimension, options);
287 up.SetLength(P->Dimension - (n-1));
288 up[0] = 1;
289 for (int i = 1; i < P->Dimension - (n-1); ++i)
290 up[i] = 0;
292 if (n == 1) {
293 gen_fun *S, *S_shift, *hp;
295 S = barvinok_enumerate_series(P, P->Dimension, options);
296 S_shift = new gen_fun(S);
297 S_shift->shift(up);
298 hp = S->Hadamard_product(S_shift, options);
299 S->add(mone, hp, options);
300 delete S_shift;
301 delete hp;
303 gf = S->summate(1, options);
304 delete S;
306 return gf;
309 U = Universe_Polyhedron(P->Dimension - n);
310 dirs = Polyhedron_Lattice_Width_Directions(P, U, options);
311 Polyhedron_Free(U);
313 for (int i = 0; i < dirs->n; ++i) {
314 Polyhedron *Pi, *R;
315 Polyhedron *CA;
316 gen_fun *S, *S_shift, *S_divide, *sum;
318 CA = align_context(dirs->wd[i].domain, P->Dimension, options->MaxRays);
319 R = DomainIntersection(P, CA, options->MaxRays);
320 Polyhedron_Free(CA);
321 assert(dirs->wd[i].dir->Size == n);
322 Pi = put_direction_last(R, dirs->wd[i].dir, options->MaxRays);
323 Polyhedron_Free(R);
325 S = project(Pi, n-1, options, 1);
327 S_shift = new gen_fun(S);
328 S_divide = new gen_fun(S);
329 S_divide->divide(up);
331 for (int j = 0; more_shifts_needed(j, n, S, S_divide, up, options); ++j) {
332 gen_fun *hp;
334 S_shift->shift(up);
335 hp = S->Hadamard_product(S_shift, options);
336 S->add(mone, hp, options);
338 delete hp;
341 sum = S->summate(1, options);
343 delete S_shift;
344 delete S_divide;
345 delete S;
347 if (!gf)
348 gf = sum;
349 else {
350 gf->add(one, sum, options);
351 delete sum;
354 free_width_direction_array(dirs);
356 return gf;
359 /* Return gf of P projected on last dim(P)-n coordinates, i.e.,
360 * project out the first n coordinates.
362 static gen_fun *project(Polyhedron *P, unsigned n, barvinok_options *options,
363 int free_P)
365 Matrix *CP = NULL;
366 gen_fun *proj;
367 unsigned nparam = P->Dimension - n;
369 if (P->NbEq != 0) {
370 Polyhedron *Q = P;
371 remove_all_equalities(&P, NULL, &CP, NULL, nparam, options->MaxRays);
372 if (CP)
373 nparam = CP->NbColumns - 1;
374 n = P->Dimension - nparam;
375 if (free_P)
376 Polyhedron_Free(Q);
377 free_P = 1;
380 if (emptyQ(P))
381 proj = new gen_fun(Empty_Polyhedron(nparam));
382 else
383 proj = project_full_dim(P, n, options);
384 if (CP) {
385 proj->substitute(CP);
386 Matrix_Free(CP);
389 if (free_P)
390 Polyhedron_Free(P);
392 return proj;
395 gen_fun *barvinok_enumerate_e_series(Polyhedron *P,
396 unsigned exist, unsigned nparam, barvinok_options *options)
398 Polyhedron *P_orig = P;
399 gen_fun *gf, *proj;
400 unsigned nvar = P->Dimension - exist - nparam;
402 if (exist == 0)
403 return barvinok_enumerate_series(P, nparam, options);
405 if (emptyQ(P))
406 return new gen_fun(Empty_Polyhedron(nparam));
408 assert(!Polyhedron_is_unbounded(P, nparam, options->MaxRays));
409 assert(P->NbBid == 0);
410 assert(Polyhedron_has_revlex_positive_rays(P, nparam));
412 /* Move existentially quantified variables to the front.*/
413 if (nvar) {
414 Matrix *T;
415 T = Matrix_Alloc(exist+nvar+nparam+1, nvar+exist+nparam+1);
416 for (int i = 0; i < exist; ++i)
417 value_set_si(T->p[i][nvar+i], 1);
418 for (int i = 0; i < nvar; ++i)
419 value_set_si(T->p[exist+i][i], 1);
420 for (int i = 0; i < nparam+1; ++i)
421 value_set_si(T->p[exist+nvar+i][nvar+exist+i], 1);
422 P = Polyhedron_Image(P, T, options->MaxRays);
423 Matrix_Free(T);
425 proj = project(P, exist, options, P != P_orig);
427 if (!nvar)
428 gf = proj;
429 else {
430 gf = proj->summate(nvar, options);
431 delete proj;
433 return gf;