?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
polymake/configure.in: allow specification of location of barvninok
[barvinok.git] / barvinok.cc
blob239da47e4aca11aeb83a6287955c286c1ce6854e
1 #include <assert.h>
2 #include <iostream>
3 #include <vector>
4 #include <deque>
5 #include <string>
6 #include <sstream>
7 #include <gmp.h>
8 #include <NTL/mat_ZZ.h>
9 #include <NTL/LLL.h>
10 #include <barvinok/util.h>
11 extern "C" {
12 #include <polylib/polylibgmp.h>
13 #include <barvinok/evalue.h>
14 #include "piputil.h"
15 }
16 #include "config.h"
17 #include <barvinok/barvinok.h>
18 #include <barvinok/genfun.h>
19 #include "conversion.h"
20 #include "decomposer.h"
21 #include "lattice_point.h"
22 #include "reduce_domain.h"
23
24 #ifdef NTL_STD_CXX
25 using namespace NTL;
26 #endif
27 using std::cerr;
28 using std::cout;
29 using std::endl;
30 using std::vector;
31 using std::deque;
32 using std::string;
33 using std::ostringstream;
34
35 #define ALLOC(t,p) p = (t*)malloc(sizeof(*p))
36
37 static void rays(mat_ZZ& r, Polyhedron *C)
38 {
39 unsigned dim = C->NbRays - 1; /* don't count zero vertex */
40 assert(C->NbRays - 1 == C->Dimension);
41 r.SetDims(dim, dim);
42 ZZ tmp;
43
44 int i, c;
45 for (i = 0, c = 0; i < dim; ++i)
46 if (value_zero_p(C->Ray[i][dim+1])) {
47 for (int j = 0; j < dim; ++j) {
48 value2zz(C->Ray[i][j+1], tmp);
49 r[j][c] = tmp;
50 }
51 ++c;
52 }
53 }
54
55 class dpoly {
56 public:
57 vec_ZZ coeff;
58 dpoly(int d, ZZ& degree, int offset = 0) {
59 coeff.SetLength(d+1);
60
61 int min = d + offset;
62 if (degree >= 0 && degree < ZZ(INIT_VAL, min))
63 min = to_int(degree);
64
65 ZZ c = ZZ(INIT_VAL, 1);
66 if (!offset)
67 coeff[0] = c;
68 for (int i = 1; i <= min; ++i) {
69 c *= (degree -i + 1);
70 c /= i;
71 coeff[i-offset] = c;
72 }
73 }
74 void operator *= (dpoly& f) {
75 assert(coeff.length() == f.coeff.length());
76 vec_ZZ old = coeff;
77 coeff = f.coeff[0] * coeff;
78 for (int i = 1; i < coeff.length(); ++i)
79 for (int j = 0; i+j < coeff.length(); ++j)
80 coeff[i+j] += f.coeff[i] * old[j];
81 }
82 void div(dpoly& d, mpq_t count, ZZ& sign) {
83 int len = coeff.length();
84 Value tmp;
85 value_init(tmp);
86 mpq_t* c = new mpq_t[coeff.length()];
87 mpq_t qtmp;
88 mpq_init(qtmp);
89 for (int i = 0; i < len; ++i) {
90 mpq_init(c[i]);
91 zz2value(coeff[i], tmp);
92 mpq_set_z(c[i], tmp);
93
94 for (int j = 1; j <= i; ++j) {
95 zz2value(d.coeff[j], tmp);
96 mpq_set_z(qtmp, tmp);
97 mpq_mul(qtmp, qtmp, c[i-j]);
98 mpq_sub(c[i], c[i], qtmp);
99 }
100
101 zz2value(d.coeff[0], tmp);
102 mpq_set_z(qtmp, tmp);
103 mpq_div(c[i], c[i], qtmp);
104 }
105 if (sign == -1)
106 mpq_sub(count, count, c[len-1]);
107 else
108 mpq_add(count, count, c[len-1]);
109
110 value_clear(tmp);
111 mpq_clear(qtmp);
112 for (int i = 0; i < len; ++i)
113 mpq_clear(c[i]);
114 delete [] c;
115 }
116 };
117
118 class dpoly_n {
119 public:
120 Matrix *coeff;
121 ~dpoly_n() {
122 Matrix_Free(coeff);
123 }
124 dpoly_n(int d, ZZ& degree_0, ZZ& degree_1, int offset = 0) {
125 Value d0, d1;
126 value_init(d0);
127 value_init(d1);
128 zz2value(degree_0, d0);
129 zz2value(degree_1, d1);
130 coeff = Matrix_Alloc(d+1, d+1+1);
131 value_set_si(coeff->p[0][0], 1);
132 value_set_si(coeff->p[0][d+1], 1);
133 for (int i = 1; i <= d; ++i) {
134 value_multiply(coeff->p[i][0], coeff->p[i-1][0], d0);
135 Vector_Combine(coeff->p[i-1], coeff->p[i-1]+1, coeff->p[i]+1,
136 d1, d0, i);
137 value_set_si(coeff->p[i][d+1], i);
138 value_multiply(coeff->p[i][d+1], coeff->p[i][d+1], coeff->p[i-1][d+1]);
139 value_decrement(d0, d0);
140 }
141 value_clear(d0);
142 value_clear(d1);
143 }
144 void div(dpoly& d, Vector *count, ZZ& sign) {
145 int len = coeff->NbRows;
146 Matrix * c = Matrix_Alloc(coeff->NbRows, coeff->NbColumns);
147 Value tmp;
148 value_init(tmp);
149 for (int i = 0; i < len; ++i) {
150 Vector_Copy(coeff->p[i], c->p[i], len+1);
151 for (int j = 1; j <= i; ++j) {
152 zz2value(d.coeff[j], tmp);
153 value_multiply(tmp, tmp, c->p[i][len]);
154 value_oppose(tmp, tmp);
155 Vector_Combine(c->p[i], c->p[i-j], c->p[i],
156 c->p[i-j][len], tmp, len);
157 value_multiply(c->p[i][len], c->p[i][len], c->p[i-j][len]);
158 }
159 zz2value(d.coeff[0], tmp);
160 value_multiply(c->p[i][len], c->p[i][len], tmp);
161 }
162 if (sign == -1) {
163 value_set_si(tmp, -1);
164 Vector_Scale(c->p[len-1], count->p, tmp, len);
165 value_assign(count->p[len], c->p[len-1][len]);
166 } else
167 Vector_Copy(c->p[len-1], count->p, len+1);
168 Vector_Normalize(count->p, len+1);
169 value_clear(tmp);
170 Matrix_Free(c);
171 }
172 };
173
174 struct dpoly_r_term {
175 int *powers;
176 ZZ coeff;
177 };
178
179 /* len: number of elements in c
180 * each element in c is the coefficient of a power of t
181 * in the MacLaurin expansion
182 */
183 struct dpoly_r {
184 vector< dpoly_r_term * > *c;
185 int len;
186 int dim;
187 ZZ denom;
188
189 void add_term(int i, int * powers, ZZ& coeff) {
190 if (coeff == 0)
191 return;
192 for (int k = 0; k < c[i].size(); ++k) {
193 if (memcmp(c[i][k]->powers, powers, dim * sizeof(int)) == 0) {
194 c[i][k]->coeff += coeff;
195 return;
196 }
197 }
198 dpoly_r_term *t = new dpoly_r_term;
199 t->powers = new int[dim];
200 memcpy(t->powers, powers, dim * sizeof(int));
201 t->coeff = coeff;
202 c[i].push_back(t);
203 }
204 dpoly_r(int len, int dim) {
205 denom = 1;
206 this->len = len;
207 this->dim = dim;
208 c = new vector< dpoly_r_term * > [len];
209 }
210 dpoly_r(dpoly& num, int dim) {
211 denom = 1;
212 len = num.coeff.length();
213 c = new vector< dpoly_r_term * > [len];
214 this->dim = dim;
215 int powers[dim];
216 memset(powers, 0, dim * sizeof(int));
217
218 for (int i = 0; i < len; ++i) {
219 ZZ coeff = num.coeff[i];
220 add_term(i, powers, coeff);
221 }
222 }
223 dpoly_r(dpoly& num, dpoly& den, int pos, int dim) {
224 denom = 1;
225 len = num.coeff.length();
226 c = new vector< dpoly_r_term * > [len];
227 this->dim = dim;
228 int powers[dim];
229
230 for (int i = 0; i < len; ++i) {
231 ZZ coeff = num.coeff[i];
232 memset(powers, 0, dim * sizeof(int));
233 powers[pos] = 1;
234
235 add_term(i, powers, coeff);
236
237 for (int j = 1; j <= i; ++j) {
238 for (int k = 0; k < c[i-j].size(); ++k) {
239 memcpy(powers, c[i-j][k]->powers, dim*sizeof(int));
240 powers[pos]++;
241 coeff = -den.coeff[j-1] * c[i-j][k]->coeff;
242 add_term(i, powers, coeff);
243 }
244 }
245 }
246 //dump();
247 }
248 dpoly_r(dpoly_r* num, dpoly& den, int pos, int dim) {
249 denom = num->denom;
250 len = num->len;
251 c = new vector< dpoly_r_term * > [len];
252 this->dim = dim;
253 int powers[dim];
254 ZZ coeff;
255
256 for (int i = 0 ; i < len; ++i) {
257 for (int k = 0; k < num->c[i].size(); ++k) {
258 memcpy(powers, num->c[i][k]->powers, dim*sizeof(int));
259 powers[pos]++;
260 add_term(i, powers, num->c[i][k]->coeff);
261 }
262
263 for (int j = 1; j <= i; ++j) {
264 for (int k = 0; k < c[i-j].size(); ++k) {
265 memcpy(powers, c[i-j][k]->powers, dim*sizeof(int));
266 powers[pos]++;
267 coeff = -den.coeff[j-1] * c[i-j][k]->coeff;
268 add_term(i, powers, coeff);
269 }
270 }
271 }
272 }
273 ~dpoly_r() {
274 for (int i = 0 ; i < len; ++i)
275 for (int k = 0; k < c[i].size(); ++k) {
276 delete [] c[i][k]->powers;
277 delete c[i][k];
278 }
279 delete [] c;
280 }
281 dpoly_r *div(dpoly& d) {
282 dpoly_r *rc = new dpoly_r(len, dim);
283 rc->denom = power(d.coeff[0], len);
284 ZZ inv_d = rc->denom / d.coeff[0];
285 ZZ coeff;
286
287 for (int i = 0; i < len; ++i) {
288 for (int k = 0; k < c[i].size(); ++k) {
289 coeff = c[i][k]->coeff * inv_d;
290 rc->add_term(i, c[i][k]->powers, coeff);
291 }
292
293 for (int j = 1; j <= i; ++j) {
294 for (int k = 0; k < rc->c[i-j].size(); ++k) {
295 coeff = - d.coeff[j] * rc->c[i-j][k]->coeff / d.coeff[0];
296 rc->add_term(i, rc->c[i-j][k]->powers, coeff);
297 }
298 }
299 }
300 return rc;
301 }
302 void dump(void) {
303 for (int i = 0; i < len; ++i) {
304 cerr << endl;
305 cerr << i << endl;
306 cerr << c[i].size() << endl;
307 for (int j = 0; j < c[i].size(); ++j) {
308 for (int k = 0; k < dim; ++k) {
309 cerr << c[i][j]->powers[k] << " ";
310 }
311 cerr << ": " << c[i][j]->coeff << "/" << denom << endl;
312 }
313 cerr << endl;
314 }
315 }
316 };
317
318 const int MAX_TRY=10;
319 /*
320 * Searches for a vector that is not orthogonal to any
321 * of the rays in rays.
322 */
323 static void nonorthog(mat_ZZ& rays, vec_ZZ& lambda)
324 {
325 int dim = rays.NumCols();
326 bool found = false;
327 lambda.SetLength(dim);
328 if (dim == 0)
329 return;
330
331 for (int i = 2; !found && i <= 50*dim; i+=4) {
332 for (int j = 0; j < MAX_TRY; ++j) {
333 for (int k = 0; k < dim; ++k) {
334 int r = random_int(i)+2;
335 int v = (2*(r%2)-1) * (r >> 1);
336 lambda[k] = v;
337 }
338 int k = 0;
339 for (; k < rays.NumRows(); ++k)
340 if (lambda * rays[k] == 0)
341 break;
342 if (k == rays.NumRows()) {
343 found = true;
344 break;
345 }
346 }
347 }
348 assert(found);
349 }
350
351 static void randomvector(Polyhedron *P, vec_ZZ& lambda, int nvar)
352 {
353 Value tmp;
354 int max = 10 * 16;
355 unsigned int dim = P->Dimension;
356 value_init(tmp);
357
358 for (int i = 0; i < P->NbRays; ++i) {
359 for (int j = 1; j <= dim; ++j) {
360 value_absolute(tmp, P->Ray[i][j]);
361 int t = VALUE_TO_LONG(tmp) * 16;
362 if (t > max)
363 max = t;
364 }
365 }
366 for (int i = 0; i < P->NbConstraints; ++i) {
367 for (int j = 1; j <= dim; ++j) {
368 value_absolute(tmp, P->Constraint[i][j]);
369 int t = VALUE_TO_LONG(tmp) * 16;
370 if (t > max)
371 max = t;
372 }
373 }
374 value_clear(tmp);
375
376 lambda.SetLength(nvar);
377 for (int k = 0; k < nvar; ++k) {
378 int r = random_int(max*dim)+2;
379 int v = (2*(r%2)-1) * (max/2*dim + (r >> 1));
380 lambda[k] = v;
381 }
382 }
383
384 static void add_rays(mat_ZZ& rays, Polyhedron *i, int *r, int nvar = -1,
385 bool all = false)
386 {
387 unsigned dim = i->Dimension;
388 if (nvar == -1)
389 nvar = dim;
390 for (int k = 0; k < i->NbRays; ++k) {
391 if (!value_zero_p(i->Ray[k][dim+1]))
392 continue;
393 if (!all && nvar != dim && First_Non_Zero(i->Ray[k]+1, nvar) == -1)
394 continue;
395 values2zz(i->Ray[k]+1, rays[(*r)++], nvar);
396 }
397 }
398
399 static void mask_r(Matrix *f, int nr, Vector *lcm, int p, Vector *val, evalue *ev)
400 {
401 unsigned nparam = lcm->Size;
402
403 if (p == nparam) {
404 Vector * prod = Vector_Alloc(f->NbRows);
405 Matrix_Vector_Product(f, val->p, prod->p);
406 int isint = 1;
407 for (int i = 0; i < nr; ++i) {
408 value_modulus(prod->p[i], prod->p[i], f->p[i][nparam+1]);
409 isint &= value_zero_p(prod->p[i]);
410 }
411 value_set_si(ev->d, 1);
412 value_init(ev->x.n);
413 value_set_si(ev->x.n, isint);
414 Vector_Free(prod);
415 return;
416 }
417
418 Value tmp;
419 value_init(tmp);
420 if (value_one_p(lcm->p[p]))
421 mask_r(f, nr, lcm, p+1, val, ev);
422 else {
423 value_assign(tmp, lcm->p[p]);
424 value_set_si(ev->d, 0);
425 ev->x.p = new_enode(periodic, VALUE_TO_INT(tmp), p+1);
426 do {
427 value_decrement(tmp, tmp);
428 value_assign(val->p[p], tmp);
429 mask_r(f, nr, lcm, p+1, val, &ev->x.p->arr[VALUE_TO_INT(tmp)]);
430 } while (value_pos_p(tmp));
431 }
432 value_clear(tmp);
433 }
434
435 #ifdef USE_MODULO
436 static void mask(Matrix *f, evalue *factor)
437 {
438 int nr = f->NbRows, nc = f->NbColumns;
439 int n;
440 bool found = false;
441 for (n = 0; n < nr && value_notzero_p(f->p[n][nc-1]); ++n)
442 if (value_notone_p(f->p[n][nc-1]) &&
443 value_notmone_p(f->p[n][nc-1]))
444 found = true;
445 if (!found)
446 return;
447
448 evalue EP;
449 nr = n;
450
451 Value m;
452 value_init(m);
453
454 evalue EV;
455 value_init(EV.d);
456 value_init(EV.x.n);
457 value_set_si(EV.x.n, 1);
458
459 for (n = 0; n < nr; ++n) {
460 value_assign(m, f->p[n][nc-1]);
461 if (value_one_p(m) || value_mone_p(m))
462 continue;
463
464 int j = normal_mod(f->p[n], nc-1, &m);
465 if (j == nc-1) {
466 free_evalue_refs(factor);
467 value_init(factor->d);
468 evalue_set_si(factor, 0, 1);
469 break;
470 }
471 vec_ZZ row;
472 values2zz(f->p[n], row, nc-1);
473 ZZ g;
474 value2zz(m, g);
475 if (j < (nc-1)-1 && row[j] > g/2) {
476 for (int k = j; k < (nc-1); ++k)
477 if (row[k] != 0)
478 row[k] = g - row[k];
479 }
480
481 value_init(EP.d);
482 value_set_si(EP.d, 0);
483 EP.x.p = new_enode(relation, 2, 0);
484 value_clear(EP.x.p->arr[1].d);
485 EP.x.p->arr[1] = *factor;
486 evalue *ev = &EP.x.p->arr[0];
487 value_set_si(ev->d, 0);
488 ev->x.p = new_enode(fractional, 3, -1);
489 evalue_set_si(&ev->x.p->arr[1], 0, 1);
490 evalue_set_si(&ev->x.p->arr[2], 1, 1);
491 evalue *E = multi_monom(row);
492 value_assign(EV.d, m);
493 emul(&EV, E);
494 value_clear(ev->x.p->arr[0].d);
495 ev->x.p->arr[0] = *E;
496 delete E;
497 *factor = EP;
498 }
499
500 value_clear(m);
501 free_evalue_refs(&EV);
502 }
503 #else
504 /*
505 *
506 */
507 static void mask(Matrix *f, evalue *factor)
508 {
509 int nr = f->NbRows, nc = f->NbColumns;
510 int n;
511 bool found = false;
512 for (n = 0; n < nr && value_notzero_p(f->p[n][nc-1]); ++n)
513 if (value_notone_p(f->p[n][nc-1]) &&
514 value_notmone_p(f->p[n][nc-1]))
515 found = true;
516 if (!found)
517 return;
518
519 Value tmp;
520 value_init(tmp);
521 nr = n;
522 unsigned np = nc - 2;
523 Vector *lcm = Vector_Alloc(np);
524 Vector *val = Vector_Alloc(nc);
525 Vector_Set(val->p, 0, nc);
526 value_set_si(val->p[np], 1);
527 Vector_Set(lcm->p, 1, np);
528 for (n = 0; n < nr; ++n) {
529 if (value_one_p(f->p[n][nc-1]) ||
530 value_mone_p(f->p[n][nc-1]))
531 continue;
532 for (int j = 0; j < np; ++j)
533 if (value_notzero_p(f->p[n][j])) {
534 Gcd(f->p[n][j], f->p[n][nc-1], &tmp);
535 value_division(tmp, f->p[n][nc-1], tmp);
536 value_lcm(tmp, lcm->p[j], &lcm->p[j]);
537 }
538 }
539 evalue EP;
540 value_init(EP.d);
541 mask_r(f, nr, lcm, 0, val, &EP);
542 value_clear(tmp);
543 Vector_Free(val);
544 Vector_Free(lcm);
545 emul(&EP,factor);
546 free_evalue_refs(&EP);
547 }
548 #endif
549
550 /* This structure encodes the power of the term in a rational generating function.
551 *
552 * Either E == NULL or constant = 0
553 * If E != NULL, then the power is E
554 * If E == NULL, then the power is coeff * param[pos] + constant
555 */
556 struct term_info {
557 evalue *E;
558 ZZ constant;
559 ZZ coeff;
560 int pos;
561 };
562
563 /* Returns the power of (t+1) in the term of a rational generating function,
564 * i.e., the scalar product of the actual lattice point and lambda.
565 * The lattice point is the unique lattice point in the fundamental parallelepiped
566 * of the unimodual cone i shifted to the parametric vertex V.
567 *
568 * PD is the parameter domain, which, if != NULL, may be used to simply the
569 * resulting expression.
570 *
571 * The result is returned in term.
572 */
573 void lattice_point(
574 Param_Vertices* V, Polyhedron *i, vec_ZZ& lambda, term_info* term,
575 Polyhedron *PD)
576 {
577 unsigned nparam = V->Vertex->NbColumns - 2;
578 unsigned dim = i->Dimension;
579 mat_ZZ vertex;
580 vertex.SetDims(V->Vertex->NbRows, nparam+1);
581 Value lcm, tmp;
582 value_init(lcm);
583 value_init(tmp);
584 value_set_si(lcm, 1);
585 for (int j = 0; j < V->Vertex->NbRows; ++j) {
586 value_lcm(lcm, V->Vertex->p[j][nparam+1], &lcm);
587 }
588 if (value_notone_p(lcm)) {
589 Matrix * mv = Matrix_Alloc(dim, nparam+1);
590 for (int j = 0 ; j < dim; ++j) {
591 value_division(tmp, lcm, V->Vertex->p[j][nparam+1]);
592 Vector_Scale(V->Vertex->p[j], mv->p[j], tmp, nparam+1);
593 }
594
595 term->E = lattice_point(i, lambda, mv, lcm, PD);
596 term->constant = 0;
597
598 Matrix_Free(mv);
599 value_clear(lcm);
600 value_clear(tmp);
601 return;
602 }
603 for (int i = 0; i < V->Vertex->NbRows; ++i) {
604 assert(value_one_p(V->Vertex->p[i][nparam+1])); // for now
605 values2zz(V->Vertex->p[i], vertex[i], nparam+1);
606 }
607
608 vec_ZZ num;
609 num = lambda * vertex;
610
611 int p = -1;
612 int nn = 0;
613 for (int j = 0; j < nparam; ++j)
614 if (num[j] != 0) {
615 ++nn;
616 p = j;
617 }
618 if (nn >= 2) {
619 term->E = multi_monom(num);
620 term->constant = 0;
621 } else {
622 term->E = NULL;
623 term->constant = num[nparam];
624 term->pos = p;
625 if (p != -1)
626 term->coeff = num[p];
627 }
628
629 value_clear(lcm);
630 value_clear(tmp);
631 }
632
633 static void normalize(ZZ& sign, ZZ& num, vec_ZZ& den)
634 {
635 unsigned dim = den.length();
636
637 int change = 0;
638
639 for (int j = 0; j < den.length(); ++j) {
640 if (den[j] > 0)
641 change ^= 1;
642 else {
643 den[j] = abs(den[j]);
644 num += den[j];
645 }
646 }
647 if (change)
648 sign = -sign;
649 }
650
651 /* input:
652 * f: the powers in the denominator for the remaining vars
653 * each row refers to a factor
654 * den_s: for each factor, the power of (s+1)
655 * sign
656 * num_s: powers in the numerator corresponding to the summed vars
657 * num_p: powers in the numerator corresponding to the remaining vars
658 * number of rays in cone: "dim" = "k"
659 * length of each ray: "dim" = "d"
660 * for now, it is assumed: k == d
661 * output:
662 * den_p: for each factor
663 * 0: independent of remaining vars
664 * 1: power corresponds to corresponding row in f
665 *
666 * all inputs are subject to change
667 */
668 static void normalize(ZZ& sign,
669 ZZ& num_s, vec_ZZ& num_p, vec_ZZ& den_s, vec_ZZ& den_p,
670 mat_ZZ& f)
671 {
672 unsigned dim = f.NumRows();
673 unsigned nparam = num_p.length();
674 unsigned nvar = dim - nparam;
675
676 int change = 0;
677
678 for (int j = 0; j < den_s.length(); ++j) {
679 if (den_s[j] == 0) {
680 den_p[j] = 1;
681 continue;
682 }
683 int k;
684 for (k = 0; k < nparam; ++k)
685 if (f[j][k] != 0)
686 break;
687 if (k < nparam) {
688 den_p[j] = 1;
689 if (den_s[j] > 0) {
690 f[j] = -f[j];
691 num_p += f[j];
692 }
693 } else
694 den_p[j] = 0;
695 if (den_s[j] > 0)
696 change ^= 1;
697 else {
698 den_s[j] = abs(den_s[j]);
699 num_s += den_s[j];
700 }
701 }
702
703 if (change)
704 sign = -sign;
705 }
706
707 struct counter : public polar_decomposer {
708 vec_ZZ lambda;
709 mat_ZZ rays;
710 vec_ZZ vertex;
711 vec_ZZ den;
712 ZZ sign;
713 ZZ num;
714 int j;
715 Polyhedron *P;
716 unsigned dim;
717 mpq_t count;
718
719 counter(Polyhedron *P) {
720 this->P = P;
721 dim = P->Dimension;
722 rays.SetDims(dim, dim);
723 den.SetLength(dim);
724 mpq_init(count);
725 }
726
727 void start(unsigned MaxRays);
728
729 ~counter() {
730 mpq_clear(count);
731 }
732
733 virtual void handle_polar(Polyhedron *P, int sign);
734 };
735
736 struct OrthogonalException {} Orthogonal;
737
738 void counter::handle_polar(Polyhedron *C, int s)
739 {
740 int r = 0;
741 assert(C->NbRays-1 == dim);
742 add_rays(rays, C, &r);
743 for (int k = 0; k < dim; ++k) {
744 if (lambda * rays[k] == 0)
745 throw Orthogonal;
746 }
747
748 sign = s;
749
750 lattice_point(P->Ray[j]+1, C, vertex);
751 num = vertex * lambda;
752 den = rays * lambda;
753 normalize(sign, num, den);
754
755 dpoly d(dim, num);
756 dpoly n(dim, den[0], 1);
757 for (int k = 1; k < dim; ++k) {
758 dpoly fact(dim, den[k], 1);
759 n *= fact;
760 }
761 d.div(n, count, sign);
762 }
763
764 void counter::start(unsigned MaxRays)
765 {
766 for (;;) {
767 try {
768 randomvector(P, lambda, dim);
769 for (j = 0; j < P->NbRays; ++j) {
770 Polyhedron *C = supporting_cone(P, j);
771 decompose(C, MaxRays);
772 }
773 break;
774 } catch (OrthogonalException &e) {
775 mpq_set_si(count, 0, 0);
776 }
777 }
778 }
779
780 /* base for non-parametric counting */
781 struct np_base : public polar_decomposer {
782 int current_vertex;
783 Polyhedron *P;
784 unsigned dim;
785
786 np_base(Polyhedron *P, unsigned dim) {
787 this->P = P;
788 this->dim = dim;
789 }
790 };
791
792 struct reducer : public np_base {
793 vec_ZZ vertex;
794 //vec_ZZ den;
795 ZZ sgn;
796 ZZ num;
797 ZZ one;
798 mpq_t tcount;
799 mpz_t tn;
800 mpz_t td;
801 int lower; // call base when only this many variables is left
802
803 reducer(Polyhedron *P) : np_base(P, P->Dimension) {
804 //den.SetLength(dim);
805 mpq_init(tcount);
806 mpz_init(tn);
807 mpz_init(td);
808 one = 1;
809 }
810
811 void start(unsigned MaxRays);
812
813 ~reducer() {
814 mpq_clear(tcount);
815 mpz_clear(tn);
816 mpz_clear(td);
817 }
818
819 virtual void handle_polar(Polyhedron *P, int sign);
820 void reduce(ZZ c, ZZ cd, vec_ZZ& num, mat_ZZ& den_f);
821 virtual void base(ZZ& c, ZZ& cd, vec_ZZ& num, mat_ZZ& den_f) = 0;
822 virtual void split(vec_ZZ& num, ZZ& num_s, vec_ZZ& num_p,
823 mat_ZZ& den_f, vec_ZZ& den_s, mat_ZZ& den_r) = 0;
824 };
825
826 void reducer::reduce(ZZ c, ZZ cd, vec_ZZ& num, mat_ZZ& den_f)
827 {
828 unsigned len = den_f.NumRows(); // number of factors in den
829
830 if (num.length() == lower) {
831 base(c, cd, num, den_f);
832 return;
833 }
834 assert(num.length() > 1);
835
836 vec_ZZ den_s;
837 mat_ZZ den_r;
838 ZZ num_s;
839 vec_ZZ num_p;
840
841 split(num, num_s, num_p, den_f, den_s, den_r);
842
843 vec_ZZ den_p;
844 den_p.SetLength(len);
845
846 normalize(c, num_s, num_p, den_s, den_p, den_r);
847
848 int only_param = 0; // k-r-s from text
849 int no_param = 0; // r from text
850 for (int k = 0; k < len; ++k) {
851 if (den_p[k] == 0)
852 ++no_param;
853 else if (den_s[k] == 0)
854 ++only_param;
855 }
856 if (no_param == 0) {
857 reduce(c, cd, num_p, den_r);
858 } else {
859 int k, l;
860 mat_ZZ pden;
861 pden.SetDims(only_param, den_r.NumCols());
862
863 for (k = 0, l = 0; k < len; ++k)
864 if (den_s[k] == 0)
865 pden[l++] = den_r[k];
866
867 for (k = 0; k < len; ++k)
868 if (den_p[k] == 0)
869 break;
870
871 dpoly n(no_param, num_s);
872 dpoly D(no_param, den_s[k], 1);
873 for ( ; ++k < len; )
874 if (den_p[k] == 0) {
875 dpoly fact(no_param, den_s[k], 1);
876 D *= fact;
877 }
878
879 if (no_param + only_param == len) {
880 mpq_set_si(tcount, 0, 1);
881 n.div(D, tcount, one);
882
883 ZZ qn, qd;
884 value2zz(mpq_numref(tcount), qn);
885 value2zz(mpq_denref(tcount), qd);
886
887 qn *= c;
888 qd *= cd;
889
890 if (qn != 0)
891 reduce(qn, qd, num_p, pden);
892 } else {
893 dpoly_r * r = 0;
894
895 for (k = 0; k < len; ++k) {
896 if (den_s[k] == 0 || den_p[k] == 0)
897 continue;
898
899 dpoly pd(no_param-1, den_s[k], 1);
900
901 int l;
902 for (l = 0; l < k; ++l)
903 if (den_r[l] == den_r[k])
904 break;
905
906 if (r == 0)
907 r = new dpoly_r(n, pd, l, len);
908 else {
909 dpoly_r *nr = new dpoly_r(r, pd, l, len);
910 delete r;
911 r = nr;
912 }
913 }
914
915 dpoly_r *rc = r->div(D);
916
917 rc->denom *= cd;
918
919 int common = pden.NumRows();
920 vector< dpoly_r_term * >& final = rc->c[rc->len-1];
921 int rows;
922 for (int j = 0; j < final.size(); ++j) {
923 if (final[j]->coeff == 0)
924 continue;
925 rows = common;
926 pden.SetDims(rows, pden.NumCols());
927 for (int k = 0; k < rc->dim; ++k) {
928 int n = final[j]->powers[k];
929 if (n == 0)
930 continue;
931 pden.SetDims(rows+n, pden.NumCols());
932 for (int l = 0; l < n; ++l)
933 pden[rows+l] = den_r[k];
934 rows += n;
935 }
936 final[j]->coeff *= c;
937 reduce(final[j]->coeff, rc->denom, num_p, pden);
938 }
939
940 delete rc;
941 delete r;
942 }
943 }
944 }
945
946 void reducer::handle_polar(Polyhedron *C, int s)
947 {
948 assert(C->NbRays-1 == dim);
949
950 sgn = s;
951
952 lattice_point(P->Ray[current_vertex]+1, C, vertex);
953
954 mat_ZZ den;
955 den.SetDims(dim, dim);
956
957 int r;
958 for (r = 0; r < dim; ++r)
959 values2zz(C->Ray[r]+1, den[r], dim);
960
961 reduce(sgn, one, vertex, den);
962 }
963
964 void reducer::start(unsigned MaxRays)
965 {
966 for (current_vertex = 0; current_vertex < P->NbRays; ++current_vertex) {
967 Polyhedron *C = supporting_cone(P, current_vertex);
968 decompose(C, MaxRays);
969 }
970 }
971
972 struct ireducer : public reducer {
973 ireducer(Polyhedron *P) : reducer(P) {}
974
975 virtual void split(vec_ZZ& num, ZZ& num_s, vec_ZZ& num_p,
976 mat_ZZ& den_f, vec_ZZ& den_s, mat_ZZ& den_r) {
977 unsigned len = den_f.NumRows(); // number of factors in den
978 unsigned d = num.length() - 1;
979
980 den_s.SetLength(len);
981 den_r.SetDims(len, d);
982
983 for (int r = 0; r < len; ++r) {
984 den_s[r] = den_f[r][0];
985 for (int k = 1; k <= d; ++k)
986 den_r[r][k-1] = den_f[r][k];
987 }
988
989 num_s = num[0];
990 num_p.SetLength(d);
991 for (int k = 1 ; k <= d; ++k)
992 num_p[k-1] = num[k];
993 }
994 };
995
996 // incremental counter
997 struct icounter : public ireducer {
998 mpq_t count;
999
1000 icounter(Polyhedron *P) : ireducer(P) {
1001 mpq_init(count);
1002 lower = 1;
1003 }
1004 ~icounter() {
1005 mpq_clear(count);
1006 }
1007 virtual void base(ZZ& c, ZZ& cd, vec_ZZ& num, mat_ZZ& den_f);
1008 };
1009
1010 void icounter::base(ZZ& c, ZZ& cd, vec_ZZ& num, mat_ZZ& den_f)
1011 {
1012 int r;
1013 unsigned len = den_f.NumRows(); // number of factors in den
1014 vec_ZZ den_s;
1015 den_s.SetLength(len);
1016 ZZ num_s = num[0];
1017 for (r = 0; r < len; ++r)
1018 den_s[r] = den_f[r][0];
1019 normalize(c, num_s, den_s);
1020
1021 dpoly n(len, num_s);
1022 dpoly D(len, den_s[0], 1);
1023 for (int k = 1; k < len; ++k) {
1024 dpoly fact(len, den_s[k], 1);
1025 D *= fact;
1026 }
1027 mpq_set_si(tcount, 0, 1);
1028 n.div(D, tcount, one);
1029 zz2value(c, tn);
1030 zz2value(cd, td);
1031 mpz_mul(mpq_numref(tcount), mpq_numref(tcount), tn);
1032 mpz_mul(mpq_denref(tcount), mpq_denref(tcount), td);
1033 mpq_canonicalize(tcount);
1034 mpq_add(count, count, tcount);
1035 }
1036
1037 /* base for generating function counting */
1038 struct gf_base {
1039 np_base *base;
1040 gen_fun *gf;
1041
1042 gf_base(np_base *npb, unsigned nparam) : base(npb) {
1043 gf = new gen_fun(Polyhedron_Project(base->P, nparam));
1044 }
1045 void start(unsigned MaxRays);
1046 };
1047
1048 void gf_base::start(unsigned MaxRays)
1049 {
1050 for (int i = 0; i < base->P->NbRays; ++i) {
1051 if (!value_pos_p(base->P->Ray[i][base->dim+1]))
1052 continue;
1053
1054 Polyhedron *C = supporting_cone(base->P, i);
1055 base->current_vertex = i;
1056 base->decompose(C, MaxRays);
1057 }
1058 }
1059
1060 struct partial_ireducer : public ireducer, public gf_base {
1061 partial_ireducer(Polyhedron *P, unsigned nparam) :
1062 ireducer(P), gf_base(this, nparam) {
1063 lower = nparam;
1064 }
1065 ~partial_ireducer() {
1066 }
1067 virtual void base(ZZ& c, ZZ& cd, vec_ZZ& num, mat_ZZ& den_f);
1068 /* we want to override the start method from reducer with the one from gf_base */
1069 void start(unsigned MaxRays) {
1070 gf_base::start(MaxRays);
1071 }
1072 };
1073
1074 void partial_ireducer::base(ZZ& c, ZZ& cd, vec_ZZ& num, mat_ZZ& den_f)
1075 {
1076 gf->add(c, cd, num, den_f);
1077 }
1078
1079 struct partial_reducer : public reducer, public gf_base {
1080 vec_ZZ lambda;
1081 vec_ZZ tmp;
1082
1083 partial_reducer(Polyhedron *P, unsigned nparam) :
1084 reducer(P), gf_base(this, nparam) {
1085 lower = nparam;
1086
1087 tmp.SetLength(dim - nparam);
1088 randomvector(P, lambda, dim - nparam);
1089 }
1090 ~partial_reducer() {
1091 }
1092 virtual void base(ZZ& c, ZZ& cd, vec_ZZ& num, mat_ZZ& den_f);
1093 /* we want to override the start method from reducer with the one from gf_base */
1094 void start(unsigned MaxRays) {
1095 gf_base::start(MaxRays);
1096 }
1097
1098 virtual void split(vec_ZZ& num, ZZ& num_s, vec_ZZ& num_p,
1099 mat_ZZ& den_f, vec_ZZ& den_s, mat_ZZ& den_r) {
1100 unsigned len = den_f.NumRows(); // number of factors in den
1101 unsigned nvar = tmp.length();
1102
1103 den_s.SetLength(len);
1104 den_r.SetDims(len, lower);
1105
1106 for (int r = 0; r < len; ++r) {
1107 for (int k = 0; k < nvar; ++k)
1108 tmp[k] = den_f[r][k];
1109 den_s[r] = tmp * lambda;
1110
1111 for (int k = nvar; k < dim; ++k)
1112 den_r[r][k-nvar] = den_f[r][k];
1113 }
1114
1115 for (int k = 0; k < nvar; ++k)
1116 tmp[k] = num[k];
1117 num_s = tmp *lambda;
1118 num_p.SetLength(lower);
1119 for (int k = nvar ; k < dim; ++k)
1120 num_p[k-nvar] = num[k];
1121 }
1122 };
1123
1124 void partial_reducer::base(ZZ& c, ZZ& cd, vec_ZZ& num, mat_ZZ& den_f)
1125 {
1126 gf->add(c, cd, num, den_f);
1127 }
1128
1129 struct bfc_term_base {
1130 // the number of times a given factor appears in the denominator
1131 int *powers;
1132 mat_ZZ terms;
1133
1134 bfc_term_base(int len) {
1135 powers = new int[len];
1136 }
1137
1138 virtual ~bfc_term_base() {
1139 delete [] powers;
1140 }
1141 };
1142
1143 struct bfc_term : public bfc_term_base {
1144 vec_ZZ cn;
1145 vec_ZZ cd;
1146
1147 bfc_term(int len) : bfc_term_base(len) {}
1148 };
1149
1150 struct bfe_term : public bfc_term_base {
1151 vector<evalue *> factors;
1152
1153 bfe_term(int len) : bfc_term_base(len) {
1154 }
1155
1156 ~bfe_term() {
1157 for (int i = 0; i < factors.size(); ++i) {
1158 if (!factors[i])
1159 continue;
1160 free_evalue_refs(factors[i]);
1161 delete factors[i];
1162 }
1163 }
1164 };
1165
1166 typedef vector< bfc_term_base * > bfc_vec;
1167
1168 struct bf_reducer;
1169
1170 struct bf_base : public np_base {
1171 ZZ one;
1172 mpq_t tcount;
1173 mpz_t tn;
1174 mpz_t td;
1175 int lower; // call base when only this many variables is left
1176
1177 bf_base(Polyhedron *P, unsigned dim) : np_base(P, dim) {
1178 mpq_init(tcount);
1179 mpz_init(tn);
1180 mpz_init(td);
1181 one = 1;
1182 }
1183
1184 ~bf_base() {
1185 mpq_clear(tcount);
1186 mpz_clear(tn);
1187 mpz_clear(td);
1188 }
1189
1190 void start(unsigned MaxRays);
1191 virtual void handle_polar(Polyhedron *P, int sign);
1192 int setup_factors(Polyhedron *P, mat_ZZ& factors, bfc_term_base* t, int s);
1193
1194 bfc_term_base* find_bfc_term(bfc_vec& v, int *powers, int len);
1195 void add_term(bfc_term_base *t, vec_ZZ& num1, vec_ZZ& num);
1196 void add_term(bfc_term_base *t, vec_ZZ& num);
1197
1198 void reduce(mat_ZZ& factors, bfc_vec& v);
1199 virtual void base(mat_ZZ& factors, bfc_vec& v) = 0;
1200
1201 virtual bfc_term_base* new_bf_term(int len) = 0;
1202 virtual void set_factor(bfc_term_base *t, int k, int change) = 0;
1203 virtual void set_factor(bfc_term_base *t, int k, mpq_t &f, int change) = 0;
1204 virtual void set_factor(bfc_term_base *t, int k, ZZ& n, ZZ& d, int change) = 0;
1205 virtual void update_term(bfc_term_base *t, int i) = 0;
1206 virtual void insert_term(bfc_term_base *t, int i) = 0;
1207 virtual bool constant_vertex(int dim) = 0;
1208 virtual void cum(bf_reducer *bfr, bfc_term_base *t, int k,
1209 dpoly_r *r) = 0;
1210 };
1211
1212 static int lex_cmp(vec_ZZ& a, vec_ZZ& b)
1213 {
1214 assert(a.length() == b.length());
1215
1216 for (int j = 0; j < a.length(); ++j)
1217 if (a[j] != b[j])
1218 return a[j] < b[j] ? -1 : 1;
1219 return 0;
1220 }
1221
1222 void bf_base::add_term(bfc_term_base *t, vec_ZZ& num_orig, vec_ZZ& extra_num)
1223 {
1224 vec_ZZ num;
1225 int d = num_orig.length();
1226 num.SetLength(d-1);
1227 for (int l = 0; l < d-1; ++l)
1228 num[l] = num_orig[l+1] + extra_num[l];
1229
1230 add_term(t, num);
1231 }
1232
1233 void bf_base::add_term(bfc_term_base *t, vec_ZZ& num)
1234 {
1235 int len = t->terms.NumRows();
1236 int i, r;
1237 for (i = 0; i < len; ++i) {
1238 r = lex_cmp(t->terms[i], num);
1239 if (r >= 0)
1240 break;
1241 }
1242 if (i == len || r > 0) {
1243 t->terms.SetDims(len+1, num.length());
1244 insert_term(t, i);
1245 t->terms[i] = num;
1246 } else {
1247 // i < len && r == 0
1248 update_term(t, i);
1249 }
1250 }
1251
1252 static void print_int_vector(int *v, int len, char *name)
1253 {
1254 cerr << name << endl;
1255 for (int j = 0; j < len; ++j) {
1256 cerr << v[j] << " ";
1257 }
1258 cerr << endl;
1259 }
1260
1261 static void print_bfc_terms(mat_ZZ& factors, bfc_vec& v)
1262 {
1263 cerr << endl;
1264 cerr << "factors" << endl;
1265 cerr << factors << endl;
1266 for (int i = 0; i < v.size(); ++i) {
1267 cerr << "term: " << i << endl;
1268 print_int_vector(v[i]->powers, factors.NumRows(), "powers");
1269 cerr << "terms" << endl;
1270 cerr << v[i]->terms << endl;
1271 bfc_term* bfct = static_cast<bfc_term *>(v[i]);
1272 cerr << bfct->cn << endl;
1273 cerr << bfct->cd << endl;
1274 }
1275 }
1276
1277 static void print_bfe_terms(mat_ZZ& factors, bfc_vec& v)
1278 {
1279 cerr << endl;
1280 cerr << "factors" << endl;
1281 cerr << factors << endl;
1282 for (int i = 0; i < v.size(); ++i) {
1283 cerr << "term: " << i << endl;
1284 print_int_vector(v[i]->powers, factors.NumRows(), "powers");
1285 cerr << "terms" << endl;
1286 cerr << v[i]->terms << endl;
1287 bfe_term* bfet = static_cast<bfe_term *>(v[i]);
1288 for (int j = 0; j < v[i]->terms.NumRows(); ++j) {
1289 char * test[] = {"a", "b"};
1290 print_evalue(stderr, bfet->factors[j], test);
1291 fprintf(stderr, "\n");
1292 }
1293 }
1294 }
1295
1296 bfc_term_base* bf_base::find_bfc_term(bfc_vec& v, int *powers, int len)
1297 {
1298 bfc_vec::iterator i;
1299 for (i = v.begin(); i != v.end(); ++i) {
1300 int j;
1301 for (j = 0; j < len; ++j)
1302 if ((*i)->powers[j] != powers[j])
1303 break;
1304 if (j == len)
1305 return (*i);
1306 if ((*i)->powers[j] > powers[j])
1307 break;
1308 }
1309
1310 bfc_term_base* t = new_bf_term(len);
1311 v.insert(i, t);
1312 memcpy(t->powers, powers, len * sizeof(int));
1313
1314 return t;
1315 }
1316
1317 struct bf_reducer {
1318 mat_ZZ& factors;
1319 bfc_vec& v;
1320 bf_base *bf;
1321
1322 unsigned nf;
1323 unsigned d;
1324
1325 mat_ZZ nfactors;
1326 int *old2new;
1327 int *sign;
1328 unsigned int nnf;
1329 bfc_vec vn;
1330
1331 vec_ZZ extra_num;
1332 int changes;
1333 int no_param; // r from text
1334 int only_param; // k-r-s from text
1335 int total_power; // k from text
1336
1337 // created in compute_reduced_factors
1338 int *bpowers;
1339 // set in update_powers
1340 int *npowers;
1341 vec_ZZ l_extra_num;
1342 int l_changes;
1343
1344 bf_reducer(mat_ZZ& factors, bfc_vec& v, bf_base *bf)
1345 : factors(factors), v(v), bf(bf) {
1346 nf = factors.NumRows();
1347 d = factors.NumCols();
1348 old2new = new int[nf];
1349 sign = new int[nf];
1350
1351 extra_num.SetLength(d-1);
1352 }
1353 ~bf_reducer() {
1354 delete [] old2new;
1355 delete [] sign;
1356 delete [] npowers;
1357 delete [] bpowers;
1358 }
1359
1360 void compute_reduced_factors();
1361 void compute_extra_num(int i);
1362
1363 void reduce();
1364
1365 void update_powers(int *powers, int len);
1366 };
1367
1368 void bf_reducer::compute_extra_num(int i)
1369 {
1370 clear(extra_num);
1371 changes = 0;
1372 no_param = 0; // r from text
1373 only_param = 0; // k-r-s from text
1374 total_power = 0; // k from text
1375
1376 for (int j = 0; j < nf; ++j) {
1377 if (v[i]->powers[j] == 0)
1378 continue;
1379
1380 total_power += v[i]->powers[j];
1381 if (factors[j][0] == 0) {
1382 only_param += v[i]->powers[j];
1383 continue;
1384 }
1385
1386 if (old2new[j] == -1)
1387 no_param += v[i]->powers[j];
1388 else
1389 extra_num += -sign[j] * v[i]->powers[j] * nfactors[old2new[j]];
1390 changes += v[i]->powers[j];
1391 }
1392 }
1393
1394 void bf_reducer::update_powers(int *powers, int len)
1395 {
1396 for (int l = 0; l < nnf; ++l)
1397 npowers[l] = bpowers[l];
1398
1399 l_extra_num = extra_num;
1400 l_changes = changes;
1401
1402 for (int l = 0; l < len; ++l) {
1403 int n = powers[l];
1404 if (n == 0)
1405 continue;
1406 assert(old2new[l] != -1);
1407
1408 npowers[old2new[l]] += n;
1409 // interpretation of sign has been inverted
1410 // since we inverted the power for specialization
1411 if (sign[l] == 1) {
1412 l_extra_num += n * nfactors[old2new[l]];
1413 l_changes += n;
1414 }
1415 }
1416 }
1417
1418
1419 void bf_reducer::compute_reduced_factors()
1420 {
1421 unsigned nf = factors.NumRows();
1422 unsigned d = factors.NumCols();
1423 nnf = 0;
1424 nfactors.SetDims(nnf, d-1);
1425
1426 for (int i = 0; i < nf; ++i) {
1427 int j;
1428 int s = 1;
1429 for (j = 0; j < nnf; ++j) {
1430 int k;
1431 for (k = 1; k < d; ++k)
1432 if (factors[i][k] != 0 || nfactors[j][k-1] != 0)
1433 break;
1434 if (k < d && factors[i][k] == -nfactors[j][k-1])
1435 s = -1;
1436 for (; k < d; ++k)
1437 if (factors[i][k] != s * nfactors[j][k-1])
1438 break;
1439 if (k == d)
1440 break;
1441 }
1442 old2new[i] = j;
1443 if (j == nnf) {
1444 int k;
1445 for (k = 1; k < d; ++k)
1446 if (factors[i][k] != 0)
1447 break;
1448 if (k < d) {
1449 if (factors[i][k] < 0)
1450 s = -1;
1451 nfactors.SetDims(++nnf, d-1);
1452 for (int k = 1; k < d; ++k)
1453 nfactors[j][k-1] = s * factors[i][k];
1454 } else
1455 old2new[i] = -1;
1456 }
1457 sign[i] = s;
1458 }
1459 npowers = new int[nnf];
1460 bpowers = new int[nnf];
1461 }
1462
1463 void bf_reducer::reduce()
1464 {
1465 compute_reduced_factors();
1466
1467 for (int i = 0; i < v.size(); ++i) {
1468 compute_extra_num(i);
1469
1470 if (no_param == 0) {
1471 vec_ZZ extra_num;
1472 extra_num.SetLength(d-1);
1473 int changes = 0;
1474 int npowers[nnf];
1475 for (int k = 0; k < nnf; ++k)
1476 npowers[k] = 0;
1477 for (int k = 0; k < nf; ++k) {
1478 assert(old2new[k] != -1);
1479 npowers[old2new[k]] += v[i]->powers[k];
1480 if (sign[k] == -1) {
1481 extra_num += v[i]->powers[k] * nfactors[old2new[k]];
1482 changes += v[i]->powers[k];
1483 }
1484 }
1485
1486 bfc_term_base * t = bf->find_bfc_term(vn, npowers, nnf);
1487 for (int k = 0; k < v[i]->terms.NumRows(); ++k) {
1488 bf->set_factor(v[i], k, changes % 2);
1489 bf->add_term(t, v[i]->terms[k], extra_num);
1490 }
1491 } else {
1492 // powers of "constant" part
1493 for (int k = 0; k < nnf; ++k)
1494 bpowers[k] = 0;
1495 for (int k = 0; k < nf; ++k) {
1496 if (factors[k][0] != 0)
1497 continue;
1498 assert(old2new[k] != -1);
1499 bpowers[old2new[k]] += v[i]->powers[k];
1500 if (sign[k] == -1) {
1501 extra_num += v[i]->powers[k] * nfactors[old2new[k]];
1502 changes += v[i]->powers[k];
1503 }
1504 }
1505
1506 int j;
1507 for (j = 0; j < nf; ++j)
1508 if (old2new[j] == -1 && v[i]->powers[j] > 0)
1509 break;
1510
1511 dpoly D(no_param, factors[j][0], 1);
1512 for (int k = 1; k < v[i]->powers[j]; ++k) {
1513 dpoly fact(no_param, factors[j][0], 1);
1514 D *= fact;
1515 }
1516 for ( ; ++j < nf; )
1517 if (old2new[j] == -1)
1518 for (int k = 0; k < v[i]->powers[j]; ++k) {
1519 dpoly fact(no_param, factors[j][0], 1);
1520 D *= fact;
1521 }
1522
1523 if (no_param + only_param == total_power &&
1524 bf->constant_vertex(d)) {
1525 bfc_term_base * t = NULL;
1526 vec_ZZ num;
1527 num.SetLength(d-1);
1528 ZZ cn;
1529 ZZ cd;
1530 for (int k = 0; k < v[i]->terms.NumRows(); ++k) {
1531 dpoly n(no_param, v[i]->terms[k][0]);
1532 mpq_set_si(bf->tcount, 0, 1);
1533 n.div(D, bf->tcount, bf->one);
1534
1535 if (value_zero_p(mpq_numref(bf->tcount)))
1536 continue;
1537
1538 if (!t)
1539 t = bf->find_bfc_term(vn, bpowers, nnf);
1540 bf->set_factor(v[i], k, bf->tcount, changes % 2);
1541 bf->add_term(t, v[i]->terms[k], extra_num);
1542 }
1543 } else {
1544 for (int j = 0; j < v[i]->terms.NumRows(); ++j) {
1545 dpoly n(no_param, v[i]->terms[j][0]);
1546
1547 dpoly_r * r = 0;
1548 if (no_param + only_param == total_power)
1549 r = new dpoly_r(n, nf);
1550 else
1551 for (int k = 0; k < nf; ++k) {
1552 if (v[i]->powers[k] == 0)
1553 continue;
1554 if (factors[k][0] == 0 || old2new[k] == -1)
1555 continue;
1556
1557 dpoly pd(no_param-1, factors[k][0], 1);
1558
1559 for (int l = 0; l < v[i]->powers[k]; ++l) {
1560 int q;
1561 for (q = 0; q < k; ++q)
1562 if (old2new[q] == old2new[k] &&
1563 sign[q] == sign[k])
1564 break;
1565
1566 if (r == 0)
1567 r = new dpoly_r(n, pd, q, nf);
1568 else {
1569 dpoly_r *nr = new dpoly_r(r, pd, q, nf);
1570 delete r;
1571 r = nr;
1572 }
1573 }
1574 }
1575
1576 dpoly_r *rc = r->div(D);
1577 delete r;
1578
1579 if (bf->constant_vertex(d)) {
1580 vector< dpoly_r_term * >& final = rc->c[rc->len-1];
1581
1582 for (int k = 0; k < final.size(); ++k) {
1583 if (final[k]->coeff == 0)
1584 continue;
1585
1586 update_powers(final[k]->powers, rc->dim);
1587
1588 bfc_term_base * t = bf->find_bfc_term(vn, npowers, nnf);
1589 bf->set_factor(v[i], j, final[k]->coeff, rc->denom, l_changes % 2);
1590 bf->add_term(t, v[i]->terms[j], l_extra_num);
1591 }
1592 } else
1593 bf->cum(this, v[i], j, rc);
1594
1595 delete rc;
1596 }
1597 }
1598 }
1599 delete v[i];
1600 }
1601
1602 }
1603
1604 void bf_base::reduce(mat_ZZ& factors, bfc_vec& v)
1605 {
1606 assert(v.size() > 0);
1607 unsigned nf = factors.NumRows();
1608 unsigned d = factors.NumCols();
1609
1610 if (d == lower)
1611 return base(factors, v);
1612
1613 bf_reducer bfr(factors, v, this);
1614
1615 bfr.reduce();
1616
1617 if (bfr.vn.size() > 0)
1618 reduce(bfr.nfactors, bfr.vn);
1619 }
1620
1621 int bf_base::setup_factors(Polyhedron *C, mat_ZZ& factors,
1622 bfc_term_base* t, int s)
1623 {
1624 factors.SetDims(dim, dim);
1625
1626 int r;
1627
1628 for (r = 0; r < dim; ++r)
1629 t->powers[r] = 1;
1630
1631 for (r = 0; r < dim; ++r) {
1632 values2zz(C->Ray[r]+1, factors[r], dim);
1633 int k;
1634 for (k = 0; k < dim; ++k)
1635 if (factors[r][k] != 0)
1636 break;
1637 if (factors[r][k] < 0) {
1638 factors[r] = -factors[r];
1639 t->terms[0] += factors[r];
1640 s = -s;
1641 }
1642 }
1643
1644 return s;
1645 }
1646
1647 void bf_base::handle_polar(Polyhedron *C, int s)
1648 {
1649 bfc_term* t = new bfc_term(dim);
1650 vector< bfc_term_base * > v;
1651 v.push_back(t);
1652
1653 assert(C->NbRays-1 == dim);
1654
1655 t->cn.SetLength(1);
1656 t->cd.SetLength(1);
1657
1658 t->terms.SetDims(1, dim);
1659 lattice_point(P->Ray[current_vertex]+1, C, t->terms[0]);
1660
1661 // the elements of factors are always lexpositive
1662 mat_ZZ factors;
1663 s = setup_factors(C, factors, t, s);
1664
1665 t->cn[0] = s;
1666 t->cd[0] = 1;
1667
1668 reduce(factors, v);
1669 }
1670
1671 void bf_base::start(unsigned MaxRays)
1672 {
1673 for (current_vertex = 0; current_vertex < P->NbRays; ++current_vertex) {
1674 Polyhedron *C = supporting_cone(P, current_vertex);
1675 decompose(C, MaxRays);
1676 }
1677 }
1678
1679 struct bfcounter_base : public bf_base {
1680 ZZ cn;
1681 ZZ cd;
1682
1683 bfcounter_base(Polyhedron *P) : bf_base(P, P->Dimension) {
1684 }
1685
1686 bfc_term_base* new_bf_term(int len) {
1687 bfc_term* t = new bfc_term(len);
1688 t->cn.SetLength(0);
1689 t->cd.SetLength(0);
1690 return t;
1691 }
1692
1693 virtual void set_factor(bfc_term_base *t, int k, int change) {
1694 bfc_term* bfct = static_cast<bfc_term *>(t);
1695 cn = bfct->cn[k];
1696 if (change)
1697 cn = -cn;
1698 cd = bfct->cd[k];
1699 }
1700
1701 virtual void set_factor(bfc_term_base *t, int k, mpq_t &f, int change) {
1702 bfc_term* bfct = static_cast<bfc_term *>(t);
1703 value2zz(mpq_numref(f), cn);
1704 value2zz(mpq_denref(f), cd);
1705 cn *= bfct->cn[k];
1706 if (change)
1707 cn = -cn;
1708 cd *= bfct->cd[k];
1709 }
1710
1711 virtual void set_factor(bfc_term_base *t, int k, ZZ& n, ZZ& d, int change) {
1712 bfc_term* bfct = static_cast<bfc_term *>(t);
1713 cn = bfct->cn[k] * n;
1714 if (change)
1715 cn = -cn;
1716 cd = bfct->cd[k] * d;
1717 }
1718
1719 virtual void insert_term(bfc_term_base *t, int i) {
1720 bfc_term* bfct = static_cast<bfc_term *>(t);
1721 int len = t->terms.NumRows()-1; // already increased by one
1722
1723 bfct->cn.SetLength(len+1);
1724 bfct->cd.SetLength(len+1);
1725 for (int j = len; j > i; --j) {
1726 bfct->cn[j] = bfct->cn[j-1];
1727 bfct->cd[j] = bfct->cd[j-1];
1728 t->terms[j] = t->terms[j-1];
1729 }
1730 bfct->cn[i] = cn;
1731 bfct->cd[i] = cd;
1732 }
1733
1734 virtual void update_term(bfc_term_base *t, int i) {
1735 bfc_term* bfct = static_cast<bfc_term *>(t);
1736
1737 ZZ g = GCD(bfct->cd[i], cd);
1738 ZZ n = cn * bfct->cd[i]/g + bfct->cn[i] * cd/g;
1739 ZZ d = bfct->cd[i] * cd / g;
1740 bfct->cn[i] = n;
1741 bfct->cd[i] = d;
1742 }
1743
1744 virtual bool constant_vertex(int dim) { return true; }
1745 virtual void cum(bf_reducer *bfr, bfc_term_base *t, int k, dpoly_r *r) {
1746 assert(0);
1747 }
1748 };
1749
1750 struct bfcounter : public bfcounter_base {
1751 mpq_t count;
1752
1753 bfcounter(Polyhedron *P) : bfcounter_base(P) {
1754 mpq_init(count);
1755 lower = 1;
1756 }
1757 ~bfcounter() {
1758 mpq_clear(count);
1759 }
1760 virtual void base(mat_ZZ& factors, bfc_vec& v);
1761 };
1762
1763 void bfcounter::base(mat_ZZ& factors, bfc_vec& v)
1764 {
1765 unsigned nf = factors.NumRows();
1766
1767 for (int i = 0; i < v.size(); ++i) {
1768 bfc_term* bfct = static_cast<bfc_term *>(v[i]);
1769 int total_power = 0;
1770 // factor is always positive, so we always
1771 // change signs
1772 for (int k = 0; k < nf; ++k)
1773 total_power += v[i]->powers[k];
1774
1775 int j;
1776 for (j = 0; j < nf; ++j)
1777 if (v[i]->powers[j] > 0)
1778 break;
1779
1780 dpoly D(total_power, factors[j][0], 1);
1781 for (int k = 1; k < v[i]->powers[j]; ++k) {
1782 dpoly fact(total_power, factors[j][0], 1);
1783 D *= fact;
1784 }
1785 for ( ; ++j < nf; )
1786 for (int k = 0; k < v[i]->powers[j]; ++k) {
1787 dpoly fact(total_power, factors[j][0], 1);
1788 D *= fact;
1789 }
1790
1791 for (int k = 0; k < v[i]->terms.NumRows(); ++k) {
1792 dpoly n(total_power, v[i]->terms[k][0]);
1793 mpq_set_si(tcount, 0, 1);
1794 n.div(D, tcount, one);
1795 if (total_power % 2)
1796 bfct->cn[k] = -bfct->cn[k];
1797 zz2value(bfct->cn[k], tn);
1798 zz2value(bfct->cd[k], td);
1799
1800 mpz_mul(mpq_numref(tcount), mpq_numref(tcount), tn);
1801 mpz_mul(mpq_denref(tcount), mpq_denref(tcount), td);
1802 mpq_canonicalize(tcount);
1803 mpq_add(count, count, tcount);
1804 }
1805 delete v[i];
1806 }
1807 }
1808
1809 struct partial_bfcounter : public bfcounter_base, public gf_base {
1810 partial_bfcounter(Polyhedron *P, unsigned nparam) :
1811 bfcounter_base(P), gf_base(this, nparam) {
1812 lower = nparam;
1813 }
1814 ~partial_bfcounter() {
1815 }
1816 virtual void base(mat_ZZ& factors, bfc_vec& v);
1817 /* we want to override the start method from bf_base with the one from gf_base */
1818 void start(unsigned MaxRays) {
1819 gf_base::start(MaxRays);
1820 }
1821 };
1822
1823 void partial_bfcounter::base(mat_ZZ& factors, bfc_vec& v)
1824 {
1825 mat_ZZ den;
1826 unsigned nf = factors.NumRows();
1827
1828 for (int i = 0; i < v.size(); ++i) {
1829 bfc_term* bfct = static_cast<bfc_term *>(v[i]);
1830 den.SetDims(0, lower);
1831 int total_power = 0;
1832 int p = 0;
1833 for (int j = 0; j < nf; ++j) {
1834 total_power += v[i]->powers[j];
1835 den.SetDims(total_power, lower);
1836 for (int k = 0; k < v[i]->powers[j]; ++k)
1837 den[p++] = factors[j];
1838 }
1839 for (int j = 0; j < v[i]->terms.NumRows(); ++j)
1840 gf->add(bfct->cn[j], bfct->cd[j], v[i]->terms[j], den);
1841 delete v[i];
1842 }
1843 }
1844
1845
1846 typedef Polyhedron * Polyhedron_p;
1847
1848 static void barvinok_count_f(Polyhedron *P, Value* result, unsigned NbMaxCons);
1849
1850 void barvinok_count(Polyhedron *P, Value* result, unsigned NbMaxCons)
1851 {
1852 unsigned dim;
1853 int allocated = 0;
1854 Polyhedron *Q;
1855 int r = 0;
1856 bool infinite = false;
1857
1858 if (emptyQ2(P)) {
1859 value_set_si(*result, 0);
1860 return;
1861 }
1862 if (P->NbEq != 0) {
1863 P = remove_equalities(P);
1864 if (emptyQ(P)) {
1865 Polyhedron_Free(P);
1866 value_set_si(*result, 0);
1867 return;
1868 }
1869 allocated = 1;
1870 }
1871 if (P->NbBid == 0)
1872 for (; r < P->NbRays; ++r)
1873 if (value_zero_p(P->Ray[r][P->Dimension+1]))
1874 break;
1875 if (P->NbBid != 0 || r < P->NbRays) {
1876 value_set_si(*result, -1);
1877 if (allocated)
1878 Polyhedron_Free(P);
1879 return;
1880 }
1881 if (P->Dimension == 0) {
1882 /* Test whether the constraints are satisfied */
1883 POL_ENSURE_VERTICES(P);
1884 value_set_si(*result, !emptyQ(P));
1885 if (allocated)
1886 Polyhedron_Free(P);
1887 return;
1888 }
1889 Q = Polyhedron_Factor(P, 0, NbMaxCons);
1890 if (Q) {
1891 if (allocated)
1892 Polyhedron_Free(P);
1893 P = Q;
1894 allocated = 1;
1895 }
1896
1897 barvinok_count_f(P, result, NbMaxCons);
1898 if (value_neg_p(*result))
1899 infinite = true;
1900 if (Q && P->next && value_notzero_p(*result)) {
1901 Value factor;
1902 value_init(factor);
1903
1904 for (Q = P->next; Q; Q = Q->next) {
1905 barvinok_count_f(Q, &factor, NbMaxCons);
1906 if (value_neg_p(factor)) {
1907 infinite = true;
1908 continue;
1909 } else if (Q->next && value_zero_p(factor)) {
1910 value_set_si(*result, 0);
1911 break;
1912 }
1913 value_multiply(*result, *result, factor);
1914 }
1915
1916 value_clear(factor);
1917 }
1918
1919 if (allocated)
1920 Domain_Free(P);
1921 if (infinite)
1922 value_set_si(*result, -1);
1923 }
1924
1925 static void barvinok_count_f(Polyhedron *P, Value* result, unsigned NbMaxCons)
1926 {
1927 if (emptyQ2(P)) {
1928 value_set_si(*result, 0);
1929 return;
1930 }
1931
1932 if (P->Dimension == 1)
1933 return Line_Length(P, result);
1934
1935 int c = P->NbConstraints;
1936 POL_ENSURE_FACETS(P);
1937 if (c != P->NbConstraints || P->NbEq != 0)
1938 return barvinok_count(P, result, NbMaxCons);
1939
1940 POL_ENSURE_VERTICES(P);
1941
1942 #ifdef USE_INCREMENTAL_BF
1943 bfcounter cnt(P);
1944 #elif defined USE_INCREMENTAL_DF
1945 icounter cnt(P);
1946 #else
1947 counter cnt(P);
1948 #endif
1949 cnt.start(NbMaxCons);
1950
1951 assert(value_one_p(&cnt.count[0]._mp_den));
1952 value_assign(*result, &cnt.count[0]._mp_num);
1953 }
1954
1955 static void uni_polynom(int param, Vector *c, evalue *EP)
1956 {
1957 unsigned dim = c->Size-2;
1958 value_init(EP->d);
1959 value_set_si(EP->d,0);
1960 EP->x.p = new_enode(polynomial, dim+1, param+1);
1961 for (int j = 0; j <= dim; ++j)
1962 evalue_set(&EP->x.p->arr[j], c->p[j], c->p[dim+1]);
1963 }
1964
1965 static void multi_polynom(Vector *c, evalue* X, evalue *EP)
1966 {
1967 unsigned dim = c->Size-2;
1968 evalue EC;
1969
1970 value_init(EC.d);
1971 evalue_set(&EC, c->p[dim], c->p[dim+1]);
1972
1973 value_init(EP->d);
1974 evalue_set(EP, c->p[dim], c->p[dim+1]);
1975
1976 for (int i = dim-1; i >= 0; --i) {
1977 emul(X, EP);
1978 value_assign(EC.x.n, c->p[i]);
1979 eadd(&EC, EP);
1980 }
1981 free_evalue_refs(&EC);
1982 }
1983
1984 Polyhedron *unfringe (Polyhedron *P, unsigned MaxRays)
1985 {
1986 int len = P->Dimension+2;
1987 Polyhedron *T, *R = P;
1988 Value g;
1989 value_init(g);
1990 Vector *row = Vector_Alloc(len);
1991 value_set_si(row->p[0], 1);
1992
1993 R = DomainConstraintSimplify(Polyhedron_Copy(P), MaxRays);
1994
1995 Matrix *M = Matrix_Alloc(2, len-1);
1996 value_set_si(M->p[1][len-2], 1);
1997 for (int v = 0; v < P->Dimension; ++v) {
1998 value_set_si(M->p[0][v], 1);
1999 Polyhedron *I = Polyhedron_Image(P, M, 2+1);
2000 value_set_si(M->p[0][v], 0);
2001 for (int r = 0; r < I->NbConstraints; ++r) {
2002 if (value_zero_p(I->Constraint[r][0]))
2003 continue;
2004 if (value_zero_p(I->Constraint[r][1]))
2005 continue;
2006 if (value_one_p(I->Constraint[r][1]))
2007 continue;
2008 if (value_mone_p(I->Constraint[r][1]))
2009 continue;
2010 value_absolute(g, I->Constraint[r][1]);
2011 Vector_Set(row->p+1, 0, len-2);
2012 value_division(row->p[1+v], I->Constraint[r][1], g);
2013 mpz_fdiv_q(row->p[len-1], I->Constraint[r][2], g);
2014 T = R;
2015 R = AddConstraints(row->p, 1, R, MaxRays);
2016 if (T != P)
2017 Polyhedron_Free(T);
2018 }
2019 Polyhedron_Free(I);
2020 }
2021 Matrix_Free(M);
2022 Vector_Free(row);
2023 value_clear(g);
2024 return R;
2025 }
2026
2027 /* this procedure may have false negatives */
2028 static bool Polyhedron_is_infinite(Polyhedron *P, unsigned nparam)
2029 {
2030 int r;
2031 for (r = 0; r < P->NbRays; ++r) {
2032 if (!value_zero_p(P->Ray[r][0]) &&
2033 !value_zero_p(P->Ray[r][P->Dimension+1]))
2034 continue;
2035 if (First_Non_Zero(P->Ray[r]+1+P->Dimension-nparam, nparam) == -1)
2036 return true;
2037 }
2038 return false;
2039 }
2040
2041 /* Check whether all rays point in the positive directions
2042 * for the parameters
2043 */
2044 static bool Polyhedron_has_positive_rays(Polyhedron *P, unsigned nparam)
2045 {
2046 int r;
2047 for (r = 0; r < P->NbRays; ++r)
2048 if (value_zero_p(P->Ray[r][P->Dimension+1])) {
2049 int i;
2050 for (i = P->Dimension - nparam; i < P->Dimension; ++i)
2051 if (value_neg_p(P->Ray[r][i+1]))
2052 return false;
2053 }
2054 return true;
2055 }
2056
2057 typedef evalue * evalue_p;
2058
2059 struct enumerator : public polar_decomposer {
2060 vec_ZZ lambda;
2061 unsigned dim, nbV;
2062 evalue ** vE;
2063 int _i;
2064 mat_ZZ rays;
2065 vec_ZZ den;
2066 ZZ sign;
2067 Polyhedron *P;
2068 Param_Vertices *V;
2069 term_info num;
2070 Vector *c;
2071 mpq_t count;
2072
2073 enumerator(Polyhedron *P, unsigned dim, unsigned nbV) {
2074 this->P = P;
2075 this->dim = dim;
2076 this->nbV = nbV;
2077 randomvector(P, lambda, dim);
2078 rays.SetDims(dim, dim);
2079 den.SetLength(dim);
2080 c = Vector_Alloc(dim+2);
2081
2082 vE = new evalue_p[nbV];
2083 for (int j = 0; j < nbV; ++j)
2084 vE[j] = 0;
2085
2086 mpq_init(count);
2087 }
2088
2089 void decompose_at(Param_Vertices *V, int _i, unsigned MaxRays) {
2090 Polyhedron *C = supporting_cone_p(P, V);
2091 this->_i = _i;
2092 this->V = V;
2093
2094 vE[_i] = new evalue;
2095 value_init(vE[_i]->d);
2096 evalue_set_si(vE[_i], 0, 1);
2097
2098 decompose(C, MaxRays);
2099 }
2100
2101 ~enumerator() {
2102 mpq_clear(count);
2103 Vector_Free(c);
2104
2105 for (int j = 0; j < nbV; ++j)
2106 if (vE[j]) {
2107 free_evalue_refs(vE[j]);
2108 delete vE[j];
2109 }
2110 delete [] vE;
2111 }
2112
2113 virtual void handle_polar(Polyhedron *P, int sign);
2114 };
2115
2116 void enumerator::handle_polar(Polyhedron *C, int s)
2117 {
2118 int r = 0;
2119 assert(C->NbRays-1 == dim);
2120 add_rays(rays, C, &r);
2121 for (int k = 0; k < dim; ++k) {
2122 if (lambda * rays[k] == 0)
2123 throw Orthogonal;
2124 }
2125
2126 sign = s;
2127
2128 lattice_point(V, C, lambda, &num, 0);
2129 den = rays * lambda;
2130 normalize(sign, num.constant, den);
2131
2132 dpoly n(dim, den[0], 1);
2133 for (int k = 1; k < dim; ++k) {
2134 dpoly fact(dim, den[k], 1);
2135 n *= fact;
2136 }
2137 if (num.E != NULL) {
2138 ZZ one(INIT_VAL, 1);
2139 dpoly_n d(dim, num.constant, one);
2140 d.div(n, c, sign);
2141 evalue EV;
2142 multi_polynom(c, num.E, &EV);
2143 eadd(&EV , vE[_i]);
2144 free_evalue_refs(&EV);
2145 free_evalue_refs(num.E);
2146 delete num.E;
2147 } else if (num.pos != -1) {
2148 dpoly_n d(dim, num.constant, num.coeff);
2149 d.div(n, c, sign);
2150 evalue EV;
2151 uni_polynom(num.pos, c, &EV);
2152 eadd(&EV , vE[_i]);
2153 free_evalue_refs(&EV);
2154 } else {
2155 mpq_set_si(count, 0, 1);
2156 dpoly d(dim, num.constant);
2157 d.div(n, count, sign);
2158 evalue EV;
2159 value_init(EV.d);
2160 evalue_set(&EV, &count[0]._mp_num, &count[0]._mp_den);
2161 eadd(&EV , vE[_i]);
2162 free_evalue_refs(&EV);
2163 }
2164 }
2165
2166 struct enumerator_base {
2167 unsigned dim;
2168 evalue ** vE;
2169 evalue ** E_vertex;
2170 evalue mone;
2171 vertex_decomposer *vpd;
2172
2173 enumerator_base(unsigned dim, vertex_decomposer *vpd)
2174 {
2175 this->dim = dim;
2176 this->vpd = vpd;
2177
2178 vE = new evalue_p[vpd->nbV];
2179 for (int j = 0; j < vpd->nbV; ++j)
2180 vE[j] = 0;
2181
2182 E_vertex = new evalue_p[dim];
2183
2184 value_init(mone.d);
2185 evalue_set_si(&mone, -1, 1);
2186 }
2187
2188 void decompose_at(Param_Vertices *V, int _i, unsigned MaxRays/*, Polyhedron *pVD*/) {
2189 //this->pVD = pVD;
2190
2191 vE[_i] = new evalue;
2192 value_init(vE[_i]->d);
2193 evalue_set_si(vE[_i], 0, 1);
2194
2195 vpd->decompose_at_vertex(V, _i, MaxRays);
2196 }
2197
2198 ~enumerator_base() {
2199 for (int j = 0; j < vpd->nbV; ++j)
2200 if (vE[j]) {
2201 free_evalue_refs(vE[j]);
2202 delete vE[j];
2203 }
2204 delete [] vE;
2205
2206 delete [] E_vertex;
2207
2208 free_evalue_refs(&mone);
2209 }
2210
2211 evalue *E_num(int i, int d) {
2212 return E_vertex[i + (dim-d)];
2213 }
2214 };
2215
2216 struct cumulator {
2217 evalue *factor;
2218 evalue *v;
2219 dpoly_r *r;
2220
2221 cumulator(evalue *factor, evalue *v, dpoly_r *r) :
2222 factor(factor), v(v), r(r) {}
2223
2224 void cumulate();
2225
2226 virtual void add_term(int *powers, int len, evalue *f2) = 0;
2227 };
2228
2229 void cumulator::cumulate()
2230 {
2231 evalue cum; // factor * 1 * E_num[0]/1 * (E_num[0]-1)/2 *...
2232 evalue f;
2233 evalue t; // E_num[0] - (m-1)
2234 #ifdef USE_MODULO
2235 evalue *cst;
2236 #else
2237 evalue mone;
2238 value_init(mone.d);
2239 evalue_set_si(&mone, -1, 1);
2240 #endif
2241
2242 value_init(cum.d);
2243 evalue_copy(&cum, factor);
2244 value_init(f.d);
2245 value_init(f.x.n);
2246 value_set_si(f.d, 1);
2247 value_set_si(f.x.n, 1);
2248 value_init(t.d);
2249 evalue_copy(&t, v);
2250
2251 #ifdef USE_MODULO
2252 for (cst = &t; value_zero_p(cst->d); ) {
2253 if (cst->x.p->type == fractional)
2254 cst = &cst->x.p->arr[1];
2255 else
2256 cst = &cst->x.p->arr[0];
2257 }
2258 #endif
2259
2260 for (int m = 0; m < r->len; ++m) {
2261 if (m > 0) {
2262 if (m > 1) {
2263 value_set_si(f.d, m);
2264 emul(&f, &cum);
2265 #ifdef USE_MODULO
2266 value_subtract(cst->x.n, cst->x.n, cst->d);
2267 #else
2268 eadd(&mone, &t);
2269 #endif
2270 }
2271 emul(&t, &cum);
2272 }
2273 vector< dpoly_r_term * >& current = r->c[r->len-1-m];
2274 for (int j = 0; j < current.size(); ++j) {
2275 if (current[j]->coeff == 0)
2276 continue;
2277 evalue *f2 = new evalue;
2278 value_init(f2->d);
2279 value_init(f2->x.n);
2280 zz2value(current[j]->coeff, f2->x.n);
2281 zz2value(r->denom, f2->d);
2282 emul(&cum, f2);
2283
2284 add_term(current[j]->powers, r->dim, f2);
2285 }
2286 }
2287 free_evalue_refs(&f);
2288 free_evalue_refs(&t);
2289 free_evalue_refs(&cum);
2290 #ifndef USE_MODULO
2291 free_evalue_refs(&mone);
2292 #endif
2293 }
2294
2295 struct E_poly_term {
2296 int *powers;
2297 evalue *E;
2298 };
2299
2300 struct ie_cum : public cumulator {
2301 vector<E_poly_term *> terms;
2302
2303 ie_cum(evalue *factor, evalue *v, dpoly_r *r) : cumulator(factor, v, r) {}
2304
2305 virtual void add_term(int *powers, int len, evalue *f2);
2306 };
2307
2308 void ie_cum::add_term(int *powers, int len, evalue *f2)
2309 {
2310 int k;
2311 for (k = 0; k < terms.size(); ++k) {
2312 if (memcmp(terms[k]->powers, powers, len * sizeof(int)) == 0) {
2313 eadd(f2, terms[k]->E);
2314 free_evalue_refs(f2);
2315 delete f2;
2316 break;
2317 }
2318 }
2319 if (k >= terms.size()) {
2320 E_poly_term *ET = new E_poly_term;
2321 ET->powers = new int[len];
2322 memcpy(ET->powers, powers, len * sizeof(int));
2323 ET->E = f2;
2324 terms.push_back(ET);
2325 }
2326 }
2327
2328 struct ienumerator : public polar_decomposer, public vertex_decomposer,
2329 public enumerator_base {
2330 //Polyhedron *pVD;
2331 mat_ZZ den;
2332 vec_ZZ vertex;
2333 mpq_t tcount;
2334
2335 ienumerator(Polyhedron *P, unsigned dim, unsigned nbV) :
2336 vertex_decomposer(P, nbV, this), enumerator_base(dim, this) {
2337 vertex.SetLength(dim);
2338
2339 den.SetDims(dim, dim);
2340 mpq_init(tcount);
2341 }
2342
2343 ~ienumerator() {
2344 mpq_clear(tcount);
2345 }
2346
2347 virtual void handle_polar(Polyhedron *P, int sign);
2348 void reduce(evalue *factor, vec_ZZ& num, mat_ZZ& den_f);
2349 };
2350
2351 void ienumerator::reduce(
2352 evalue *factor, vec_ZZ& num, mat_ZZ& den_f)
2353 {
2354 unsigned len = den_f.NumRows(); // number of factors in den
2355 unsigned dim = num.length();
2356
2357 if (dim == 0) {
2358 eadd(factor, vE[vert]);
2359 return;
2360 }
2361
2362 vec_ZZ den_s;
2363 den_s.SetLength(len);
2364 mat_ZZ den_r;
2365 den_r.SetDims(len, dim-1);
2366
2367 int r, k;
2368
2369 for (r = 0; r < len; ++r) {
2370 den_s[r] = den_f[r][0];
2371 for (k = 0; k <= dim-1; ++k)
2372 if (k != 0)
2373 den_r[r][k-(k>0)] = den_f[r][k];
2374 }
2375
2376 ZZ num_s = num[0];
2377 vec_ZZ num_p;
2378 num_p.SetLength(dim-1);
2379 for (k = 0 ; k <= dim-1; ++k)
2380 if (k != 0)
2381 num_p[k-(k>0)] = num[k];
2382
2383 vec_ZZ den_p;
2384 den_p.SetLength(len);
2385
2386 ZZ one;
2387 one = 1;
2388 normalize(one, num_s, num_p, den_s, den_p, den_r);
2389 if (one != 1)
2390 emul(&mone, factor);
2391
2392 int only_param = 0;
2393 int no_param = 0;
2394 for (int k = 0; k < len; ++k) {
2395 if (den_p[k] == 0)
2396 ++no_param;
2397 else if (den_s[k] == 0)
2398 ++only_param;
2399 }
2400 if (no_param == 0) {
2401 reduce(factor, num_p, den_r);
2402 } else {
2403 int k, l;
2404 mat_ZZ pden;
2405 pden.SetDims(only_param, dim-1);
2406
2407 for (k = 0, l = 0; k < len; ++k)
2408 if (den_s[k] == 0)
2409 pden[l++] = den_r[k];
2410
2411 for (k = 0; k < len; ++k)
2412 if (den_p[k] == 0)
2413 break;
2414
2415 dpoly n(no_param, num_s);
2416 dpoly D(no_param, den_s[k], 1);
2417 for ( ; ++k < len; )
2418 if (den_p[k] == 0) {
2419 dpoly fact(no_param, den_s[k], 1);
2420 D *= fact;
2421 }
2422
2423 dpoly_r * r = 0;
2424 // if no_param + only_param == len then all powers
2425 // below will be all zero
2426 if (no_param + only_param == len) {
2427 if (E_num(0, dim) != 0)
2428 r = new dpoly_r(n, len);
2429 else {
2430 mpq_set_si(tcount, 0, 1);
2431 one = 1;
2432 n.div(D, tcount, one);
2433
2434 if (value_notzero_p(mpq_numref(tcount))) {
2435 evalue f;
2436 value_init(f.d);
2437 value_init(f.x.n);
2438 value_assign(f.x.n, mpq_numref(tcount));
2439 value_assign(f.d, mpq_denref(tcount));
2440 emul(&f, factor);
2441 reduce(factor, num_p, pden);
2442 free_evalue_refs(&f);
2443 }
2444 return;
2445 }
2446 } else {
2447 for (k = 0; k < len; ++k) {
2448 if (den_s[k] == 0 || den_p[k] == 0)
2449 continue;
2450
2451 dpoly pd(no_param-1, den_s[k], 1);
2452
2453 int l;
2454 for (l = 0; l < k; ++l)
2455 if (den_r[l] == den_r[k])
2456 break;
2457
2458 if (r == 0)
2459 r = new dpoly_r(n, pd, l, len);
2460 else {
2461 dpoly_r *nr = new dpoly_r(r, pd, l, len);
2462 delete r;
2463 r = nr;
2464 }
2465 }
2466 }
2467 dpoly_r *rc = r->div(D);
2468 delete r;
2469 r = rc;
2470 if (E_num(0, dim) == 0) {
2471 int common = pden.NumRows();
2472 vector< dpoly_r_term * >& final = r->c[r->len-1];
2473 int rows;
2474 evalue t;
2475 evalue f;
2476 value_init(f.d);
2477 value_init(f.x.n);
2478 zz2value(r->denom, f.d);
2479 for (int j = 0; j < final.size(); ++j) {
2480 if (final[j]->coeff == 0)
2481 continue;
2482 rows = common;
2483 for (int k = 0; k < r->dim; ++k) {
2484 int n = final[j]->powers[k];
2485 if (n == 0)
2486 continue;
2487 pden.SetDims(rows+n, pden.NumCols());
2488 for (int l = 0; l < n; ++l)
2489 pden[rows+l] = den_r[k];
2490 rows += n;
2491 }
2492 value_init(t.d);
2493 evalue_copy(&t, factor);
2494 zz2value(final[j]->coeff, f.x.n);
2495 emul(&f, &t);
2496 reduce(&t, num_p, pden);
2497 free_evalue_refs(&t);
2498 }
2499 free_evalue_refs(&f);
2500 } else {
2501 ie_cum cum(factor, E_num(0, dim), r);
2502 cum.cumulate();
2503
2504 int common = pden.NumRows();
2505 int rows;
2506 for (int j = 0; j < cum.terms.size(); ++j) {
2507 rows = common;
2508 pden.SetDims(rows, pden.NumCols());
2509 for (int k = 0; k < r->dim; ++k) {
2510 int n = cum.terms[j]->powers[k];
2511 if (n == 0)
2512 continue;
2513 pden.SetDims(rows+n, pden.NumCols());
2514 for (int l = 0; l < n; ++l)
2515 pden[rows+l] = den_r[k];
2516 rows += n;
2517 }
2518 reduce(cum.terms[j]->E, num_p, pden);
2519 free_evalue_refs(cum.terms[j]->E);
2520 delete cum.terms[j]->E;
2521 delete [] cum.terms[j]->powers;
2522 delete cum.terms[j];
2523 }
2524 }
2525 delete r;
2526 }
2527 }
2528
2529 static int type_offset(enode *p)
2530 {
2531 return p->type == fractional ? 1 :
2532 p->type == flooring ? 1 : 0;
2533 }
2534
2535 static int edegree(evalue *e)
2536 {
2537 int d = 0;
2538 enode *p;
2539
2540 if (value_notzero_p(e->d))
2541 return 0;
2542
2543 p = e->x.p;
2544 int i = type_offset(p);
2545 if (p->size-i-1 > d)
2546 d = p->size - i - 1;
2547 for (; i < p->size; i++) {
2548 int d2 = edegree(&p->arr[i]);
2549 if (d2 > d)
2550 d = d2;
2551 }
2552 return d;
2553 }
2554
2555 void ienumerator::handle_polar(Polyhedron *C, int s)
2556 {
2557 assert(C->NbRays-1 == dim);
2558
2559 lattice_point(V, C, vertex, E_vertex);
2560
2561 int r;
2562 for (r = 0; r < dim; ++r)
2563 values2zz(C->Ray[r]+1, den[r], dim);
2564
2565 evalue one;
2566 value_init(one.d);
2567 evalue_set_si(&one, s, 1);
2568 reduce(&one, vertex, den);
2569 free_evalue_refs(&one);
2570
2571 for (int i = 0; i < dim; ++i)
2572 if (E_vertex[i]) {
2573 free_evalue_refs(E_vertex[i]);
2574 delete E_vertex[i];
2575 }
2576 }
2577
2578 struct bfenumerator : public vertex_decomposer, public bf_base,
2579 public enumerator_base {
2580 evalue *factor;
2581
2582 bfenumerator(Polyhedron *P, unsigned dim, unsigned nbV) :
2583 vertex_decomposer(P, nbV, this),
2584 bf_base(P, dim), enumerator_base(dim, this) {
2585 lower = 0;
2586 factor = NULL;
2587 }
2588
2589 ~bfenumerator() {
2590 }
2591
2592 virtual void handle_polar(Polyhedron *P, int sign);
2593 virtual void base(mat_ZZ& factors, bfc_vec& v);
2594
2595 bfc_term_base* new_bf_term(int len) {
2596 bfe_term* t = new bfe_term(len);
2597 return t;
2598 }
2599
2600 virtual void set_factor(bfc_term_base *t, int k, int change) {
2601 bfe_term* bfet = static_cast<bfe_term *>(t);
2602 factor = bfet->factors[k];
2603 assert(factor != NULL);
2604 bfet->factors[k] = NULL;
2605 if (change)
2606 emul(&mone, factor);
2607 }
2608
2609 virtual void set_factor(bfc_term_base *t, int k, mpq_t &q, int change) {
2610 bfe_term* bfet = static_cast<bfe_term *>(t);
2611 factor = bfet->factors[k];
2612 assert(factor != NULL);
2613 bfet->factors[k] = NULL;
2614
2615 evalue f;
2616 value_init(f.d);
2617 value_init(f.x.n);
2618 if (change)
2619 value_oppose(f.x.n, mpq_numref(q));
2620 else
2621 value_assign(f.x.n, mpq_numref(q));
2622 value_assign(f.d, mpq_denref(q));
2623 emul(&f, factor);
2624 }
2625
2626 virtual void set_factor(bfc_term_base *t, int k, ZZ& n, ZZ& d, int change) {
2627 bfe_term* bfet = static_cast<bfe_term *>(t);
2628
2629 factor = new evalue;
2630
2631 evalue f;
2632 value_init(f.d);
2633 value_init(f.x.n);
2634 zz2value(n, f.x.n);
2635 if (change)
2636 value_oppose(f.x.n, f.x.n);
2637 zz2value(d, f.d);
2638
2639 value_init(factor->d);
2640 evalue_copy(factor, bfet->factors[k]);
2641 emul(&f, factor);
2642 }
2643
2644 void set_factor(evalue *f, int change) {
2645 if (change)
2646 emul(&mone, f);
2647 factor = f;
2648 }
2649
2650 virtual void insert_term(bfc_term_base *t, int i) {
2651 bfe_term* bfet = static_cast<bfe_term *>(t);
2652 int len = t->terms.NumRows()-1; // already increased by one
2653
2654 bfet->factors.resize(len+1);
2655 for (int j = len; j > i; --j) {
2656 bfet->factors[j] = bfet->factors[j-1];
2657 t->terms[j] = t->terms[j-1];
2658 }
2659 bfet->factors[i] = factor;
2660 factor = NULL;
2661 }
2662
2663 virtual void update_term(bfc_term_base *t, int i) {
2664 bfe_term* bfet = static_cast<bfe_term *>(t);
2665
2666 eadd(factor, bfet->factors[i]);
2667 free_evalue_refs(factor);
2668 delete factor;
2669 }
2670
2671 virtual bool constant_vertex(int dim) { return E_num(0, dim) == 0; }
2672
2673 virtual void cum(bf_reducer *bfr, bfc_term_base *t, int k, dpoly_r *r);
2674 };
2675
2676 struct bfe_cum : public cumulator {
2677 bfenumerator *bfe;
2678 bfc_term_base *told;
2679 int k;
2680 bf_reducer *bfr;
2681
2682 bfe_cum(evalue *factor, evalue *v, dpoly_r *r, bf_reducer *bfr,
2683 bfc_term_base *t, int k, bfenumerator *e) :
2684 cumulator(factor, v, r), told(t), k(k),
2685 bfr(bfr), bfe(e) {
2686 }
2687
2688 virtual void add_term(int *powers, int len, evalue *f2);
2689 };
2690
2691 void bfe_cum::add_term(int *powers, int len, evalue *f2)
2692 {
2693 bfr->update_powers(powers, len);
2694
2695 bfc_term_base * t = bfe->find_bfc_term(bfr->vn, bfr->npowers, bfr->nnf);
2696 bfe->set_factor(f2, bfr->l_changes % 2);
2697 bfe->add_term(t, told->terms[k], bfr->l_extra_num);
2698 }
2699
2700 void bfenumerator::cum(bf_reducer *bfr, bfc_term_base *t, int k,
2701 dpoly_r *r)
2702 {
2703 bfe_term* bfet = static_cast<bfe_term *>(t);
2704 bfe_cum cum(bfet->factors[k], E_num(0, bfr->d), r, bfr, t, k, this);
2705 cum.cumulate();
2706 }
2707
2708 void bfenumerator::base(mat_ZZ& factors, bfc_vec& v)
2709 {
2710 for (int i = 0; i < v.size(); ++i) {
2711 assert(v[i]->terms.NumRows() == 1);
2712 evalue *factor = static_cast<bfe_term *>(v[i])->factors[0];
2713 eadd(factor, vE[vert]);
2714 delete v[i];
2715 }
2716 }
2717
2718 void bfenumerator::handle_polar(Polyhedron *C, int s)
2719 {
2720 assert(C->NbRays-1 == enumerator_base::dim);
2721
2722 bfe_term* t = new bfe_term(enumerator_base::dim);
2723 vector< bfc_term_base * > v;
2724 v.push_back(t);
2725
2726 t->factors.resize(1);
2727
2728 t->terms.SetDims(1, enumerator_base::dim);
2729 lattice_point(V, C, t->terms[0], E_vertex);
2730
2731 // the elements of factors are always lexpositive
2732 mat_ZZ factors;
2733 s = setup_factors(C, factors, t, s);
2734
2735 t->factors[0] = new evalue;
2736 value_init(t->factors[0]->d);
2737 evalue_set_si(t->factors[0], s, 1);
2738 reduce(factors, v);
2739
2740 for (int i = 0; i < enumerator_base::dim; ++i)
2741 if (E_vertex[i]) {
2742 free_evalue_refs(E_vertex[i]);
2743 delete E_vertex[i];
2744 }
2745 }
2746
2747 #ifdef HAVE_CORRECT_VERTICES
2748 static inline Param_Polyhedron *Polyhedron2Param_SD(Polyhedron **Din,
2749 Polyhedron *Cin,int WS,Polyhedron **CEq,Matrix **CT)
2750 {
2751 if (WS & POL_NO_DUAL)
2752 WS = 0;
2753 return Polyhedron2Param_SimplifiedDomain(Din, Cin, WS, CEq, CT);
2754 }
2755 #else
2756 static Param_Polyhedron *Polyhedron2Param_SD(Polyhedron **Din,
2757 Polyhedron *Cin,int WS,Polyhedron **CEq,Matrix **CT)
2758 {
2759 static char data[] = " 1 0 0 0 0 1 -18 "
2760 " 1 0 0 -20 0 19 1 "
2761 " 1 0 1 20 0 -20 16 "
2762 " 1 0 0 0 0 -1 19 "
2763 " 1 0 -1 0 0 0 4 "
2764 " 1 4 -20 0 0 -1 23 "
2765 " 1 -4 20 0 0 1 -22 "
2766 " 1 0 1 0 20 -20 16 "
2767 " 1 0 0 0 -20 19 1 ";
2768 static int checked = 0;
2769 if (!checked) {
2770 checked = 1;
2771 char *p = data;
2772 int n, v, i;
2773 Matrix *M = Matrix_Alloc(9, 7);
2774 for (i = 0; i < 9; ++i)
2775 for (int j = 0; j < 7; ++j) {
2776 sscanf(p, "%d%n", &v, &n);
2777 p += n;
2778 value_set_si(M->p[i][j], v);
2779 }
2780 Polyhedron *P = Constraints2Polyhedron(M, 1024);
2781 Matrix_Free(M);
2782 Polyhedron *U = Universe_Polyhedron(1);
2783 Param_Polyhedron *PP = Polyhedron2Param_Domain(P, U, 1024);
2784 Polyhedron_Free(P);
2785 Polyhedron_Free(U);
2786 Param_Vertices *V;
2787 for (i = 0, V = PP->V; V; ++i, V = V->next)
2788 ;
2789 if (PP)
2790 Param_Polyhedron_Free(PP);
2791 if (i != 10) {
2792 fprintf(stderr, "WARNING: results may be incorrect\n");
2793 fprintf(stderr,
2794 "WARNING: use latest version of PolyLib to remove this warning\n");
2795 }
2796 }
2797
2798 return Polyhedron2Param_SimplifiedDomain(Din, Cin, WS, CEq, CT);
2799 }
2800 #endif
2801
2802 static evalue* barvinok_enumerate_ev_f(Polyhedron *P, Polyhedron* C,
2803 unsigned MaxRays);
2804
2805 /* Destroys C */
2806 static evalue* barvinok_enumerate_cst(Polyhedron *P, Polyhedron* C,
2807 unsigned MaxRays)
2808 {
2809 evalue *eres;
2810
2811 ALLOC(evalue, eres);
2812 value_init(eres->d);
2813 value_set_si(eres->d, 0);
2814 eres->x.p = new_enode(partition, 2, C->Dimension);
2815 EVALUE_SET_DOMAIN(eres->x.p->arr[0], DomainConstraintSimplify(C, MaxRays));
2816 value_set_si(eres->x.p->arr[1].d, 1);
2817 value_init(eres->x.p->arr[1].x.n);
2818 if (emptyQ(P))
2819 value_set_si(eres->x.p->arr[1].x.n, 0);
2820 else
2821 barvinok_count(P, &eres->x.p->arr[1].x.n, MaxRays);
2822
2823 return eres;
2824 }
2825
2826 evalue* barvinok_enumerate_ev(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
2827 {
2828 //P = unfringe(P, MaxRays);
2829 Polyhedron *Corig = C;
2830 Polyhedron *CEq = NULL, *rVD, *CA;
2831 int r = 0;
2832 unsigned nparam = C->Dimension;
2833 evalue *eres;
2834
2835 evalue factor;
2836 value_init(factor.d);
2837 evalue_set_si(&factor, 1, 1);
2838
2839 CA = align_context(C, P->Dimension, MaxRays);
2840 P = DomainIntersection(P, CA, MaxRays);
2841 Polyhedron_Free(CA);
2842
2843 /* for now */
2844 POL_ENSURE_FACETS(P);
2845 POL_ENSURE_VERTICES(P);
2846 POL_ENSURE_FACETS(C);
2847 POL_ENSURE_VERTICES(C);
2848
2849 if (C->Dimension == 0 || emptyQ(P)) {
2850 constant:
2851 eres = barvinok_enumerate_cst(P, CEq ? CEq : Polyhedron_Copy(C),
2852 MaxRays);
2853 out:
2854 emul(&factor, eres);
2855 reduce_evalue(eres);
2856 free_evalue_refs(&factor);
2857 Domain_Free(P);
2858 if (C != Corig)
2859 Polyhedron_Free(C);
2860
2861 return eres;
2862 }
2863 if (Polyhedron_is_infinite(P, nparam))
2864 goto constant;
2865
2866 if (P->NbEq != 0) {
2867 Matrix *f;
2868 P = remove_equalities_p(P, P->Dimension-nparam, &f);
2869 mask(f, &factor);
2870 Matrix_Free(f);
2871 }
2872 if (P->Dimension == nparam) {
2873 CEq = P;
2874 P = Universe_Polyhedron(0);
2875 goto constant;
2876 }
2877
2878 Polyhedron *T = Polyhedron_Factor(P, nparam, MaxRays);
2879 if (T || (P->Dimension == nparam+1)) {
2880 Polyhedron *Q;
2881 Polyhedron *C2;
2882 for (Q = T ? T : P; Q; Q = Q->next) {
2883 Polyhedron *next = Q->next;
2884 Q->next = NULL;
2885
2886 Polyhedron *QC = Q;
2887 if (Q->Dimension != C->Dimension)
2888 QC = Polyhedron_Project(Q, nparam);
2889
2890 C2 = C;
2891 C = DomainIntersection(C, QC, MaxRays);
2892 if (C2 != Corig)
2893 Polyhedron_Free(C2);
2894 if (QC != Q)
2895 Polyhedron_Free(QC);
2896
2897 Q->next = next;
2898 }
2899 }
2900 if (T) {
2901 Polyhedron_Free(P);
2902 P = T;
2903 if (T->Dimension == C->Dimension) {
2904 P = T->next;
2905 T->next = NULL;
2906 Polyhedron_Free(T);
2907 }
2908 }
2909
2910 Polyhedron *next = P->next;
2911 P->next = NULL;
2912 eres = barvinok_enumerate_ev_f(P, C, MaxRays);
2913 P->next = next;
2914
2915 if (P->next) {
2916 Polyhedron *Q;
2917 evalue *f;
2918
2919 for (Q = P->next; Q; Q = Q->next) {
2920 Polyhedron *next = Q->next;
2921 Q->next = NULL;
2922
2923 f = barvinok_enumerate_ev_f(Q, C, MaxRays);
2924 emul(f, eres);
2925 free_evalue_refs(f);
2926 free(f);
2927
2928 Q->next = next;
2929 }
2930 }
2931
2932 goto out;
2933 }
2934
2935 static evalue* barvinok_enumerate_ev_f(Polyhedron *P, Polyhedron* C,
2936 unsigned MaxRays)
2937 {
2938 unsigned nparam = C->Dimension;
2939
2940 if (P->Dimension - nparam == 1)
2941 return ParamLine_Length(P, C, MaxRays);
2942
2943 Param_Polyhedron *PP = NULL;
2944 Polyhedron *CEq = NULL, *pVD;
2945 Matrix *CT = NULL;
2946 Param_Domain *D, *next;
2947 Param_Vertices *V;
2948 evalue *eres;
2949 Polyhedron *Porig = P;
2950
2951 PP = Polyhedron2Param_SD(&P,C,MaxRays,&CEq,&CT);
2952
2953 if (isIdentity(CT)) {
2954 Matrix_Free(CT);
2955 CT = NULL;
2956 } else {
2957 assert(CT->NbRows != CT->NbColumns);
2958 if (CT->NbRows == 1) { // no more parameters
2959 eres = barvinok_enumerate_cst(P, CEq, MaxRays);
2960 out:
2961 if (CT)
2962 Matrix_Free(CT);
2963 if (PP)
2964 Param_Polyhedron_Free(PP);
2965 if (P != Porig)
2966 Polyhedron_Free(P);
2967
2968 return eres;
2969 }
2970 nparam = CT->NbRows - 1;
2971 }
2972
2973 unsigned dim = P->Dimension - nparam;
2974
2975 ALLOC(evalue, eres);
2976 value_init(eres->d);
2977 value_set_si(eres->d, 0);
2978
2979 int nd;
2980 for (nd = 0, D=PP->D; D; ++nd, D=D->next);
2981 struct section { Polyhedron *D; evalue E; };
2982 section *s = new section[nd];
2983 Polyhedron **fVD = new Polyhedron_p[nd];
2984
2985 try_again:
2986 #ifdef USE_INCREMENTAL_BF
2987 bfenumerator et(P, dim, PP->nbV);
2988 #elif defined USE_INCREMENTAL_DF
2989 ienumerator et(P, dim, PP->nbV);
2990 #else
2991 enumerator et(P, dim, PP->nbV);
2992 #endif
2993
2994 for(nd = 0, D=PP->D; D; D=next) {
2995 next = D->next;
2996
2997 Polyhedron *rVD = reduce_domain(D->Domain, CT, CEq,
2998 fVD, nd, MaxRays);
2999 if (!rVD)
3000 continue;
3001
3002 pVD = CT ? DomainImage(rVD,CT,MaxRays) : rVD;
3003
3004 value_init(s[nd].E.d);
3005 evalue_set_si(&s[nd].E, 0, 1);
3006 s[nd].D = rVD;
3007
3008 FORALL_PVertex_in_ParamPolyhedron(V,D,PP) // _i is internal counter
3009 if (!et.vE[_i])
3010 try {
3011 et.decompose_at(V, _i, MaxRays);
3012 } catch (OrthogonalException &e) {
3013 if (rVD != pVD)
3014 Domain_Free(pVD);
3015 for (; nd >= 0; --nd) {
3016 free_evalue_refs(&s[nd].E);
3017 Domain_Free(s[nd].D);
3018 Domain_Free(fVD[nd]);
3019 }
3020 goto try_again;
3021 }
3022 eadd(et.vE[_i] , &s[nd].E);
3023 END_FORALL_PVertex_in_ParamPolyhedron;
3024 evalue_range_reduction_in_domain(&s[nd].E, pVD);
3025
3026 if (CT)
3027 addeliminatedparams_evalue(&s[nd].E, CT);
3028 ++nd;
3029 if (rVD != pVD)
3030 Domain_Free(pVD);
3031 }
3032
3033 if (nd == 0)
3034 evalue_set_si(eres, 0, 1);
3035 else {
3036 eres->x.p = new_enode(partition, 2*nd, C->Dimension);
3037 for (int j = 0; j < nd; ++j) {
3038 EVALUE_SET_DOMAIN(eres->x.p->arr[2*j], s[j].D);
3039 value_clear(eres->x.p->arr[2*j+1].d);
3040 eres->x.p->arr[2*j+1] = s[j].E;
3041 Domain_Free(fVD[j]);
3042 }
3043 }
3044 delete [] s;
3045 delete [] fVD;
3046
3047 if (CEq)
3048 Polyhedron_Free(CEq);
3049 goto out;
3050 }
3051
3052 Enumeration* barvinok_enumerate(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
3053 {
3054 evalue *EP = barvinok_enumerate_ev(P, C, MaxRays);
3055
3056 return partition2enumeration(EP);
3057 }
3058
3059 static void SwapColumns(Value **V, int n, int i, int j)
3060 {
3061 for (int r = 0; r < n; ++r)
3062 value_swap(V[r][i], V[r][j]);
3063 }
3064
3065 static void SwapColumns(Polyhedron *P, int i, int j)
3066 {
3067 SwapColumns(P->Constraint, P->NbConstraints, i, j);
3068 SwapColumns(P->Ray, P->NbRays, i, j);
3069 }
3070
3071 /* Construct a constraint c from constraints l and u such that if
3072 * if constraint c holds then for each value of the other variables
3073 * there is at most one value of variable pos (position pos+1 in the constraints).
3074 *
3075 * Given a lower and an upper bound
3076 * n_l v_i + <c_l,x> + c_l >= 0
3077 * -n_u v_i + <c_u,x> + c_u >= 0
3078 * the constructed constraint is
3079 *
3080 * -(n_l<c_u,x> + n_u<c_l,x>) + (-n_l c_u - n_u c_l + n_l n_u - 1)
3081 *
3082 * which is then simplified to remove the content of the non-constant coefficients
3083 *
3084 * len is the total length of the constraints.
3085 * v is a temporary variable that can be used by this procedure
3086 */
3087 static void negative_test_constraint(Value *l, Value *u, Value *c, int pos,
3088 int len, Value *v)
3089 {
3090 value_oppose(*v, u[pos+1]);
3091 Vector_Combine(l+1, u+1, c+1, *v, l[pos+1], len-1);
3092 value_multiply(*v, *v, l[pos+1]);
3093 value_subtract(c[len-1], c[len-1], *v);
3094 value_set_si(*v, -1);
3095 Vector_Scale(c+1, c+1, *v, len-1);
3096 value_decrement(c[len-1], c[len-1]);
3097 ConstraintSimplify(c, c, len, v);
3098 }
3099
3100 static bool parallel_constraints(Value *l, Value *u, Value *c, int pos,
3101 int len)
3102 {
3103 bool parallel;
3104 Value g1;
3105 Value g2;
3106 value_init(g1);
3107 value_init(g2);
3108
3109 Vector_Gcd(&l[1+pos], len, &g1);
3110 Vector_Gcd(&u[1+pos], len, &g2);
3111 Vector_Combine(l+1+pos, u+1+pos, c+1, g2, g1, len);
3112 parallel = First_Non_Zero(c+1, len) == -1;
3113
3114 value_clear(g1);
3115 value_clear(g2);
3116
3117 return parallel;
3118 }
3119
3120 static void negative_test_constraint7(Value *l, Value *u, Value *c, int pos,
3121 int exist, int len, Value *v)
3122 {
3123 Value g;
3124 value_init(g);
3125
3126 Vector_Gcd(&u[1+pos], exist, v);
3127 Vector_Gcd(&l[1+pos], exist, &g);
3128 Vector_Combine(l+1, u+1, c+1, *v, g, len-1);
3129 value_multiply(*v, *v, g);
3130 value_subtract(c[len-1], c[len-1], *v);
3131 value_set_si(*v, -1);
3132 Vector_Scale(c+1, c+1, *v, len-1);
3133 value_decrement(c[len-1], c[len-1]);
3134 ConstraintSimplify(c, c, len, v);
3135
3136 value_clear(g);
3137 }
3138
3139 /* Turns a x + b >= 0 into a x + b <= -1
3140 *
3141 * len is the total length of the constraint.
3142 * v is a temporary variable that can be used by this procedure
3143 */
3144 static void oppose_constraint(Value *c, int len, Value *v)
3145 {
3146 value_set_si(*v, -1);
3147 Vector_Scale(c+1, c+1, *v, len-1);
3148 value_decrement(c[len-1], c[len-1]);
3149 }
3150
3151 /* Split polyhedron P into two polyhedra *pos and *neg, where
3152 * existential variable i has at most one solution for each
3153 * value of the other variables in *neg.
3154 *
3155 * The splitting is performed using constraints l and u.
3156 *
3157 * nvar: number of set variables
3158 * row: temporary vector that can be used by this procedure
3159 * f: temporary value that can be used by this procedure
3160 */
3161 static bool SplitOnConstraint(Polyhedron *P, int i, int l, int u,
3162 int nvar, int MaxRays, Vector *row, Value& f,
3163 Polyhedron **pos, Polyhedron **neg)
3164 {
3165 negative_test_constraint(P->Constraint[l], P->Constraint[u],
3166 row->p, nvar+i, P->Dimension+2, &f);
3167 *neg = AddConstraints(row->p, 1, P, MaxRays);
3168
3169 /* We found an independent, but useless constraint
3170 * Maybe we should detect this earlier and not
3171 * mark the variable as INDEPENDENT
3172 */
3173 if (emptyQ((*neg))) {
3174 Polyhedron_Free(*neg);
3175 return false;
3176 }
3177
3178 oppose_constraint(row->p, P->Dimension+2, &f);
3179 *pos = AddConstraints(row->p, 1, P, MaxRays);
3180
3181 if (emptyQ((*pos))) {
3182 Polyhedron_Free(*neg);
3183 Polyhedron_Free(*pos);
3184 return false;
3185 }
3186
3187 return true;
3188 }
3189
3190 /*
3191 * unimodularly transform P such that constraint r is transformed
3192 * into a constraint that involves only a single (the first)
3193 * existential variable
3194 *
3195 */
3196 static Polyhedron *rotate_along(Polyhedron *P, int r, int nvar, int exist,
3197 unsigned MaxRays)
3198 {
3199 Value g;
3200 value_init(g);
3201
3202 Vector *row = Vector_Alloc(exist);
3203 Vector_Copy(P->Constraint[r]+1+nvar, row->p, exist);
3204 Vector_Gcd(row->p, exist, &g);
3205 if (value_notone_p(g))
3206 Vector_AntiScale(row->p, row->p, g, exist);
3207 value_clear(g);
3208
3209 Matrix *M = unimodular_complete(row);
3210 Matrix *M2 = Matrix_Alloc(P->Dimension+1, P->Dimension+1);
3211 for (r = 0; r < nvar; ++r)
3212 value_set_si(M2->p[r][r], 1);
3213 for ( ; r < nvar+exist; ++r)
3214 Vector_Copy(M->p[r-nvar], M2->p[r]+nvar, exist);
3215 for ( ; r < P->Dimension+1; ++r)
3216 value_set_si(M2->p[r][r], 1);
3217 Polyhedron *T = Polyhedron_Image(P, M2, MaxRays);
3218
3219 Matrix_Free(M2);
3220 Matrix_Free(M);
3221 Vector_Free(row);
3222
3223 return T;
3224 }
3225
3226 /* Split polyhedron P into two polyhedra *pos and *neg, where
3227 * existential variable i has at most one solution for each
3228 * value of the other variables in *neg.
3229 *
3230 * If independent is set, then the two constraints on which the
3231 * split will be performed need to be independent of the other
3232 * existential variables.
3233 *
3234 * Return true if an appropriate split could be performed.
3235 *
3236 * nvar: number of set variables
3237 * exist: number of existential variables
3238 * row: temporary vector that can be used by this procedure
3239 * f: temporary value that can be used by this procedure
3240 */
3241 static bool SplitOnVar(Polyhedron *P, int i,
3242 int nvar, int exist, int MaxRays,
3243 Vector *row, Value& f, bool independent,
3244 Polyhedron **pos, Polyhedron **neg)
3245 {
3246 int j;
3247
3248 for (int l = P->NbEq; l < P->NbConstraints; ++l) {
3249 if (value_negz_p(P->Constraint[l][nvar+i+1]))
3250 continue;
3251
3252 if (independent) {
3253 for (j = 0; j < exist; ++j)
3254 if (j != i && value_notzero_p(P->Constraint[l][nvar+j+1]))
3255 break;
3256 if (j < exist)
3257 continue;
3258 }
3259
3260 for (int u = P->NbEq; u < P->NbConstraints; ++u) {
3261 if (value_posz_p(P->Constraint[u][nvar+i+1]))
3262 continue;
3263
3264 if (independent) {
3265 for (j = 0; j < exist; ++j)
3266 if (j != i && value_notzero_p(P->Constraint[u][nvar+j+1]))
3267 break;
3268 if (j < exist)
3269 continue;
3270 }
3271
3272 if (SplitOnConstraint(P, i, l, u, nvar, MaxRays, row, f, pos, neg)) {
3273 if (independent) {
3274 if (i != 0)
3275 SwapColumns(*neg, nvar+1, nvar+1+i);
3276 }
3277 return true;
3278 }
3279 }
3280 }
3281
3282 return false;
3283 }
3284
3285 static bool double_bound_pair(Polyhedron *P, int nvar, int exist,
3286 int i, int l1, int l2,
3287 Polyhedron **pos, Polyhedron **neg)
3288 {
3289 Value f;
3290 value_init(f);
3291 Vector *row = Vector_Alloc(P->Dimension+2);
3292 value_set_si(row->p[0], 1);
3293 value_oppose(f, P->Constraint[l1][nvar+i+1]);
3294 Vector_Combine(P->Constraint[l1]+1, P->Constraint[l2]+1,
3295 row->p+1,
3296 P->Constraint[l2][nvar+i+1], f,
3297 P->Dimension+1);
3298 ConstraintSimplify(row->p, row->p, P->Dimension+2, &f);
3299 *pos = AddConstraints(row->p, 1, P, 0);
3300 value_set_si(f, -1);
3301 Vector_Scale(row->p+1, row->p+1, f, P->Dimension+1);
3302 value_decrement(row->p[P->Dimension+1], row->p[P->Dimension+1]);
3303 *neg = AddConstraints(row->p, 1, P, 0);
3304 Vector_Free(row);
3305 value_clear(f);
3306
3307 return !emptyQ((*pos)) && !emptyQ((*neg));
3308 }
3309
3310 static bool double_bound(Polyhedron *P, int nvar, int exist,
3311 Polyhedron **pos, Polyhedron **neg)
3312 {
3313 for (int i = 0; i < exist; ++i) {
3314 int l1, l2;
3315 for (l1 = P->NbEq; l1 < P->NbConstraints; ++l1) {
3316 if (value_negz_p(P->Constraint[l1][nvar+i+1]))
3317 continue;
3318 for (l2 = l1 + 1; l2 < P->NbConstraints; ++l2) {
3319 if (value_negz_p(P->Constraint[l2][nvar+i+1]))
3320 continue;
3321 if (double_bound_pair(P, nvar, exist, i, l1, l2, pos, neg))
3322 return true;
3323 }
3324 }
3325 for (l1 = P->NbEq; l1 < P->NbConstraints; ++l1) {
3326 if (value_posz_p(P->Constraint[l1][nvar+i+1]))
3327 continue;
3328 if (l1 < P->NbConstraints)
3329 for (l2 = l1 + 1; l2 < P->NbConstraints; ++l2) {
3330 if (value_posz_p(P->Constraint[l2][nvar+i+1]))
3331 continue;
3332 if (double_bound_pair(P, nvar, exist, i, l1, l2, pos, neg))
3333 return true;
3334 }
3335 }
3336 return false;
3337 }
3338 return false;
3339 }
3340
3341 enum constraint {
3342 ALL_POS = 1 << 0,
3343 ONE_NEG = 1 << 1,
3344 INDEPENDENT = 1 << 2,
3345 ROT_NEG = 1 << 3
3346 };
3347
3348 static evalue* enumerate_or(Polyhedron *D,
3349 unsigned exist, unsigned nparam, unsigned MaxRays)
3350 {
3351 #ifdef DEBUG_ER
3352 fprintf(stderr, "\nER: Or\n");
3353 #endif /* DEBUG_ER */
3354
3355 Polyhedron *N = D->next;
3356 D->next = 0;
3357 evalue *EP =
3358 barvinok_enumerate_e(D, exist, nparam, MaxRays);
3359 Polyhedron_Free(D);
3360
3361 for (D = N; D; D = N) {
3362 N = D->next;
3363 D->next = 0;
3364
3365 evalue *EN =
3366 barvinok_enumerate_e(D, exist, nparam, MaxRays);
3367
3368 eor(EN, EP);
3369 free_evalue_refs(EN);
3370 free(EN);
3371 Polyhedron_Free(D);
3372 }
3373
3374 reduce_evalue(EP);
3375
3376 return EP;
3377 }
3378
3379 static evalue* enumerate_sum(Polyhedron *P,
3380 unsigned exist, unsigned nparam, unsigned MaxRays)
3381 {
3382 int nvar = P->Dimension - exist - nparam;
3383 int toswap = nvar < exist ? nvar : exist;
3384 for (int i = 0; i < toswap; ++i)
3385 SwapColumns(P, 1 + i, nvar+exist - i);
3386 nparam += nvar;
3387
3388 #ifdef DEBUG_ER
3389 fprintf(stderr, "\nER: Sum\n");
3390 #endif /* DEBUG_ER */
3391
3392 evalue *EP = barvinok_enumerate_e(P, exist, nparam, MaxRays);
3393
3394 for (int i = 0; i < /* nvar */ nparam; ++i) {
3395 Matrix *C = Matrix_Alloc(1, 1 + nparam + 1);
3396 value_set_si(C->p[0][0], 1);
3397 evalue split;
3398 value_init(split.d);
3399 value_set_si(split.d, 0);
3400 split.x.p = new_enode(partition, 4, nparam);
3401 value_set_si(C->p[0][1+i], 1);
3402 Matrix *C2 = Matrix_Copy(C);
3403 EVALUE_SET_DOMAIN(split.x.p->arr[0],
3404 Constraints2Polyhedron(C2, MaxRays));
3405 Matrix_Free(C2);
3406 evalue_set_si(&split.x.p->arr[1], 1, 1);
3407 value_set_si(C->p[0][1+i], -1);
3408 value_set_si(C->p[0][1+nparam], -1);
3409 EVALUE_SET_DOMAIN(split.x.p->arr[2],
3410 Constraints2Polyhedron(C, MaxRays));
3411 evalue_set_si(&split.x.p->arr[3], 1, 1);
3412 emul(&split, EP);
3413 free_evalue_refs(&split);
3414 Matrix_Free(C);
3415 }
3416 reduce_evalue(EP);
3417 evalue_range_reduction(EP);
3418
3419 evalue_frac2floor(EP);
3420
3421 evalue *sum = esum(EP, nvar);
3422
3423 free_evalue_refs(EP);
3424 free(EP);
3425 EP = sum;
3426
3427 evalue_range_reduction(EP);
3428
3429 return EP;
3430 }
3431
3432 static evalue* split_sure(Polyhedron *P, Polyhedron *S,
3433 unsigned exist, unsigned nparam, unsigned MaxRays)
3434 {
3435 int nvar = P->Dimension - exist - nparam;
3436
3437 Matrix *M = Matrix_Alloc(exist, S->Dimension+2);
3438 for (int i = 0; i < exist; ++i)
3439 value_set_si(M->p[i][nvar+i+1], 1);
3440 Polyhedron *O = S;
3441 S = DomainAddRays(S, M, MaxRays);
3442 Polyhedron_Free(O);
3443 Polyhedron *F = DomainAddRays(P, M, MaxRays);
3444 Polyhedron *D = DomainDifference(F, S, MaxRays);
3445 O = D;
3446 D = Disjoint_Domain(D, 0, MaxRays);
3447 Polyhedron_Free(F);
3448 Domain_Free(O);
3449 Matrix_Free(M);
3450
3451 M = Matrix_Alloc(P->Dimension+1-exist, P->Dimension+1);
3452 for (int j = 0; j < nvar; ++j)
3453 value_set_si(M->p[j][j], 1);
3454 for (int j = 0; j < nparam+1; ++j)
3455 value_set_si(M->p[nvar+j][nvar+exist+j], 1);
3456 Polyhedron *T = Polyhedron_Image(S, M, MaxRays);
3457 evalue *EP = barvinok_enumerate_e(T, 0, nparam, MaxRays);
3458 Polyhedron_Free(S);
3459 Polyhedron_Free(T);
3460 Matrix_Free(M);
3461
3462 for (Polyhedron *Q = D; Q; Q = Q->next) {
3463 Polyhedron *N = Q->next;
3464 Q->next = 0;
3465 T = DomainIntersection(P, Q, MaxRays);
3466 evalue *E = barvinok_enumerate_e(T, exist, nparam, MaxRays);
3467 eadd(E, EP);
3468 free_evalue_refs(E);
3469 free(E);
3470 Polyhedron_Free(T);
3471 Q->next = N;
3472 }
3473 Domain_Free(D);
3474 return EP;
3475 }
3476
3477 static evalue* enumerate_sure(Polyhedron *P,
3478 unsigned exist, unsigned nparam, unsigned MaxRays)
3479 {
3480 int i;
3481 Polyhedron *S = P;
3482 int nvar = P->Dimension - exist - nparam;
3483 Value lcm;
3484 Value f;
3485 value_init(lcm);
3486 value_init(f);
3487
3488 for (i = 0; i < exist; ++i) {
3489 Matrix *M = Matrix_Alloc(S->NbConstraints, S->Dimension+2);
3490 int c = 0;
3491 value_set_si(lcm, 1);
3492 for (int j = 0; j < S->NbConstraints; ++j) {
3493 if (value_negz_p(S->Constraint[j][1+nvar+i]))
3494 continue;
3495 if (value_one_p(S->Constraint[j][1+nvar+i]))
3496 continue;
3497 value_lcm(lcm, S->Constraint[j][1+nvar+i], &lcm);
3498 }
3499
3500 for (int j = 0; j < S->NbConstraints; ++j) {
3501 if (value_negz_p(S->Constraint[j][1+nvar+i]))
3502 continue;
3503 if (value_one_p(S->Constraint[j][1+nvar+i]))
3504 continue;
3505 value_division(f, lcm, S->Constraint[j][1+nvar+i]);
3506 Vector_Scale(S->Constraint[j], M->p[c], f, S->Dimension+2);
3507 value_subtract(M->p[c][S->Dimension+1],
3508 M->p[c][S->Dimension+1],
3509 lcm);
3510 value_increment(M->p[c][S->Dimension+1],
3511 M->p[c][S->Dimension+1]);
3512 ++c;
3513 }
3514 Polyhedron *O = S;
3515 S = AddConstraints(M->p[0], c, S, MaxRays);
3516 if (O != P)
3517 Polyhedron_Free(O);
3518 Matrix_Free(M);
3519 if (emptyQ(S)) {
3520 Polyhedron_Free(S);
3521 value_clear(lcm);
3522 value_clear(f);
3523 return 0;
3524 }
3525 }
3526 value_clear(lcm);
3527 value_clear(f);
3528
3529 #ifdef DEBUG_ER
3530 fprintf(stderr, "\nER: Sure\n");
3531 #endif /* DEBUG_ER */
3532
3533 return split_sure(P, S, exist, nparam, MaxRays);
3534 }
3535
3536 static evalue* enumerate_sure2(Polyhedron *P,
3537 unsigned exist, unsigned nparam, unsigned MaxRays)
3538 {
3539 int nvar = P->Dimension - exist - nparam;
3540 int r;
3541 for (r = 0; r < P->NbRays; ++r)
3542 if (value_one_p(P->Ray[r][0]) &&
3543 value_one_p(P->Ray[r][P->Dimension+1]))
3544 break;
3545
3546 if (r >= P->NbRays)
3547 return 0;
3548
3549 Matrix *M = Matrix_Alloc(nvar + 1 + nparam, P->Dimension+2);
3550 for (int i = 0; i < nvar; ++i)
3551 value_set_si(M->p[i][1+i], 1);
3552 for (int i = 0; i < nparam; ++i)
3553 value_set_si(M->p[i+nvar][1+nvar+exist+i], 1);
3554 Vector_Copy(P->Ray[r]+1+nvar, M->p[nvar+nparam]+1+nvar, exist);
3555 value_set_si(M->p[nvar+nparam][0], 1);
3556 value_set_si(M->p[nvar+nparam][P->Dimension+1], 1);
3557 Polyhedron * F = Rays2Polyhedron(M, MaxRays);
3558 Matrix_Free(M);
3559
3560 Polyhedron *I = DomainIntersection(F, P, MaxRays);
3561 Polyhedron_Free(F);
3562
3563 #ifdef DEBUG_ER
3564 fprintf(stderr, "\nER: Sure2\n");
3565 #endif /* DEBUG_ER */
3566
3567 return split_sure(P, I, exist, nparam, MaxRays);
3568 }
3569
3570 static evalue* enumerate_cyclic(Polyhedron *P,
3571 unsigned exist, unsigned nparam,
3572 evalue * EP, int r, int p, unsigned MaxRays)
3573 {
3574 int nvar = P->Dimension - exist - nparam;
3575
3576 /* If EP in its fractional maps only contains references
3577 * to the remainder parameter with appropriate coefficients
3578 * then we could in principle avoid adding existentially
3579 * quantified variables to the validity domains.
3580 * We'd have to replace the remainder by m { p/m }
3581 * and multiply with an appropriate factor that is one
3582 * only in the appropriate range.
3583 * This last multiplication can be avoided if EP
3584 * has a single validity domain with no (further)
3585 * constraints on the remainder parameter
3586 */
3587
3588 Matrix *CT = Matrix_Alloc(nparam+1, nparam+3);
3589 Matrix *M = Matrix_Alloc(1, 1+nparam+3);
3590 for (int j = 0; j < nparam; ++j)
3591 if (j != p)
3592 value_set_si(CT->p[j][j], 1);
3593 value_set_si(CT->p[p][nparam+1], 1);
3594 value_set_si(CT->p[nparam][nparam+2], 1);
3595 value_set_si(M->p[0][1+p], -1);
3596 value_absolute(M->p[0][1+nparam], P->Ray[0][1+nvar+exist+p]);
3597 value_set_si(M->p[0][1+nparam+1], 1);
3598 Polyhedron *CEq = Constraints2Polyhedron(M, 1);
3599 Matrix_Free(M);
3600 addeliminatedparams_enum(EP, CT, CEq, MaxRays, nparam);
3601 Polyhedron_Free(CEq);
3602 Matrix_Free(CT);
3603
3604 return EP;
3605 }
3606
3607 static void enumerate_vd_add_ray(evalue *EP, Matrix *Rays, unsigned MaxRays)
3608 {
3609 if (value_notzero_p(EP->d))
3610 return;
3611
3612 assert(EP->x.p->type == partition);
3613 assert(EP->x.p->pos == EVALUE_DOMAIN(EP->x.p->arr[0])->Dimension);
3614 for (int i = 0; i < EP->x.p->size/2; ++i) {
3615 Polyhedron *D = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
3616 Polyhedron *N = DomainAddRays(D, Rays, MaxRays);
3617 EVALUE_SET_DOMAIN(EP->x.p->arr[2*i], N);
3618 Domain_Free(D);
3619 }
3620 }
3621
3622 static evalue* enumerate_line(Polyhedron *P,
3623 unsigned exist, unsigned nparam, unsigned MaxRays)
3624 {
3625 if (P->NbBid == 0)
3626 return 0;
3627
3628 #ifdef DEBUG_ER
3629 fprintf(stderr, "\nER: Line\n");
3630 #endif /* DEBUG_ER */
3631
3632 int nvar = P->Dimension - exist - nparam;
3633 int i, j;
3634 for (i = 0; i < nparam; ++i)
3635 if (value_notzero_p(P->Ray[0][1+nvar+exist+i]))
3636 break;
3637 assert(i < nparam);
3638 for (j = i+1; j < nparam; ++j)
3639 if (value_notzero_p(P->Ray[0][1+nvar+exist+i]))
3640 break;
3641 assert(j >= nparam); // for now
3642
3643 Matrix *M = Matrix_Alloc(2, P->Dimension+2);
3644 value_set_si(M->p[0][0], 1);
3645 value_set_si(M->p[0][1+nvar+exist+i], 1);
3646 value_set_si(M->p[1][0], 1);
3647 value_set_si(M->p[1][1+nvar+exist+i], -1);
3648 value_absolute(M->p[1][1+P->Dimension], P->Ray[0][1+nvar+exist+i]);
3649 value_decrement(M->p[1][1+P->Dimension], M->p[1][1+P->Dimension]);
3650 Polyhedron *S = AddConstraints(M->p[0], 2, P, MaxRays);
3651 evalue *EP = barvinok_enumerate_e(S, exist, nparam, MaxRays);
3652 Polyhedron_Free(S);
3653 Matrix_Free(M);
3654
3655 return enumerate_cyclic(P, exist, nparam, EP, 0, i, MaxRays);
3656 }
3657
3658 static int single_param_pos(Polyhedron*P, unsigned exist, unsigned nparam,
3659 int r)
3660 {
3661 int nvar = P->Dimension - exist - nparam;
3662 if (First_Non_Zero(P->Ray[r]+1, nvar) != -1)
3663 return -1;
3664 int i = First_Non_Zero(P->Ray[r]+1+nvar+exist, nparam);
3665 if (i == -1)
3666 return -1;
3667 if (First_Non_Zero(P->Ray[r]+1+nvar+exist+1, nparam-i-1) != -1)
3668 return -1;
3669 return i;
3670 }
3671
3672 static evalue* enumerate_remove_ray(Polyhedron *P, int r,
3673 unsigned exist, unsigned nparam, unsigned MaxRays)
3674 {
3675 #ifdef DEBUG_ER
3676 fprintf(stderr, "\nER: RedundantRay\n");
3677 #endif /* DEBUG_ER */
3678
3679 Value one;
3680 value_init(one);
3681 value_set_si(one, 1);
3682 int len = P->NbRays-1;
3683 Matrix *M = Matrix_Alloc(2 * len, P->Dimension+2);
3684 Vector_Copy(P->Ray[0], M->p[0], r * (P->Dimension+2));
3685 Vector_Copy(P->Ray[r+1], M->p[r], (len-r) * (P->Dimension+2));
3686 for (int j = 0; j < P->NbRays; ++j) {
3687 if (j == r)
3688 continue;
3689 Vector_Combine(P->Ray[j], P->Ray[r], M->p[len+j-(j>r)],
3690 one, P->Ray[j][P->Dimension+1], P->Dimension+2);
3691 }
3692
3693 P = Rays2Polyhedron(M, MaxRays);
3694 Matrix_Free(M);
3695 evalue *EP = barvinok_enumerate_e(P, exist, nparam, MaxRays);
3696 Polyhedron_Free(P);
3697 value_clear(one);
3698
3699 return EP;
3700 }
3701
3702 static evalue* enumerate_redundant_ray(Polyhedron *P,
3703 unsigned exist, unsigned nparam, unsigned MaxRays)
3704 {
3705 assert(P->NbBid == 0);
3706 int nvar = P->Dimension - exist - nparam;
3707 Value m;
3708 value_init(m);
3709
3710 for (int r = 0; r < P->NbRays; ++r) {
3711 if (value_notzero_p(P->Ray[r][P->Dimension+1]))
3712 continue;
3713 int i1 = single_param_pos(P, exist, nparam, r);
3714 if (i1 == -1)
3715 continue;
3716 for (int r2 = r+1; r2 < P->NbRays; ++r2) {
3717 if (value_notzero_p(P->Ray[r2][P->Dimension+1]))
3718 continue;
3719 int i2 = single_param_pos(P, exist, nparam, r2);
3720 if (i2 == -1)
3721 continue;
3722 if (i1 != i2)
3723 continue;
3724
3725 value_division(m, P->Ray[r][1+nvar+exist+i1],
3726 P->Ray[r2][1+nvar+exist+i1]);
3727 value_multiply(m, m, P->Ray[r2][1+nvar+exist+i1]);
3728 /* r2 divides r => r redundant */
3729 if (value_eq(m, P->Ray[r][1+nvar+exist+i1])) {
3730 value_clear(m);
3731 return enumerate_remove_ray(P, r, exist, nparam, MaxRays);
3732 }
3733
3734 value_division(m, P->Ray[r2][1+nvar+exist+i1],
3735 P->Ray[r][1+nvar+exist+i1]);
3736 value_multiply(m, m, P->Ray[r][1+nvar+exist+i1]);
3737 /* r divides r2 => r2 redundant */
3738 if (value_eq(m, P->Ray[r2][1+nvar+exist+i1])) {
3739 value_clear(m);
3740 return enumerate_remove_ray(P, r2, exist, nparam, MaxRays);
3741 }
3742 }
3743 }
3744 value_clear(m);
3745 return 0;
3746 }
3747
3748 static Polyhedron *upper_bound(Polyhedron *P,
3749 int pos, Value *max, Polyhedron **R)
3750 {
3751 Value v;
3752 int r;
3753 value_init(v);
3754
3755 *R = 0;
3756 Polyhedron *N;
3757 Polyhedron *B = 0;
3758 for (Polyhedron *Q = P; Q; Q = N) {
3759 N = Q->next;
3760 for (r = 0; r < P->NbRays; ++r) {
3761 if (value_zero_p(P->Ray[r][P->Dimension+1]) &&
3762 value_pos_p(P->Ray[r][1+pos]))
3763 break;
3764 }
3765 if (r < P->NbRays) {
3766 Q->next = *R;
3767 *R = Q;
3768 continue;
3769 } else {
3770 Q->next = B;
3771 B = Q;
3772 }
3773 for (r = 0; r < P->NbRays; ++r) {
3774 if (value_zero_p(P->Ray[r][P->Dimension+1]))
3775 continue;
3776 mpz_fdiv_q(v, P->Ray[r][1+pos], P->Ray[r][1+P->Dimension]);
3777 if ((!Q->next && r == 0) || value_gt(v, *max))
3778 value_assign(*max, v);
3779 }
3780 }
3781 value_clear(v);
3782 return B;
3783 }
3784
3785 static evalue* enumerate_ray(Polyhedron *P,
3786 unsigned exist, unsigned nparam, unsigned MaxRays)
3787 {
3788 assert(P->NbBid == 0);
3789 int nvar = P->Dimension - exist - nparam;
3790
3791 int r;
3792 for (r = 0; r < P->NbRays; ++r)
3793 if (value_zero_p(P->Ray[r][P->Dimension+1]))
3794 break;
3795 if (r >= P->NbRays)
3796 return 0;
3797
3798 int r2;
3799 for (r2 = r+1; r2 < P->NbRays; ++r2)
3800 if (value_zero_p(P->Ray[r2][P->Dimension+1]))
3801 break;
3802 if (r2 < P->NbRays) {
3803 if (nvar > 0)
3804 return enumerate_sum(P, exist, nparam, MaxRays);
3805 }
3806
3807 #ifdef DEBUG_ER
3808 fprintf(stderr, "\nER: Ray\n");
3809 #endif /* DEBUG_ER */
3810
3811 Value m;
3812 Value one;
3813 value_init(m);
3814 value_init(one);
3815 value_set_si(one, 1);
3816 int i = single_param_pos(P, exist, nparam, r);
3817 assert(i != -1); // for now;
3818
3819 Matrix *M = Matrix_Alloc(P->NbRays, P->Dimension+2);
3820 for (int j = 0; j < P->NbRays; ++j) {
3821 Vector_Combine(P->Ray[j], P->Ray[r], M->p[j],
3822 one, P->Ray[j][P->Dimension+1], P->Dimension+2);
3823 }
3824 Polyhedron *S = Rays2Polyhedron(M, MaxRays);
3825 Matrix_Free(M);
3826 Polyhedron *D = DomainDifference(P, S, MaxRays);
3827 Polyhedron_Free(S);
3828 // Polyhedron_Print(stderr, P_VALUE_FMT, D);
3829 assert(value_pos_p(P->Ray[r][1+nvar+exist+i])); // for now
3830 Polyhedron *R;
3831 D = upper_bound(D, nvar+exist+i, &m, &R);
3832 assert(D);
3833 Domain_Free(D);
3834
3835 M = Matrix_Alloc(2, P->Dimension+2);
3836 value_set_si(M->p[0][0], 1);
3837 value_set_si(M->p[1][0], 1);
3838 value_set_si(M->p[0][1+nvar+exist+i], -1);
3839 value_set_si(M->p[1][1+nvar+exist+i], 1);
3840 value_assign(M->p[0][1+P->Dimension], m);
3841 value_oppose(M->p[1][1+P->Dimension], m);
3842 value_addto(M->p[1][1+P->Dimension], M->p[1][1+P->Dimension],
3843 P->Ray[r][1+nvar+exist+i]);
3844 value_decrement(M->p[1][1+P->Dimension], M->p[1][1+P->Dimension]);
3845 // Matrix_Print(stderr, P_VALUE_FMT, M);
3846 D = AddConstraints(M->p[0], 2, P, MaxRays);
3847 // Polyhedron_Print(stderr, P_VALUE_FMT, D);
3848 value_subtract(M->p[0][1+P->Dimension], M->p[0][1+P->Dimension],
3849 P->Ray[r][1+nvar+exist+i]);
3850 // Matrix_Print(stderr, P_VALUE_FMT, M);
3851 S = AddConstraints(M->p[0], 1, P, MaxRays);
3852 // Polyhedron_Print(stderr, P_VALUE_FMT, S);
3853 Matrix_Free(M);
3854
3855 evalue *EP = barvinok_enumerate_e(D, exist, nparam, MaxRays);
3856 Polyhedron_Free(D);
3857 value_clear(one);
3858 value_clear(m);
3859
3860 if (value_notone_p(P->Ray[r][1+nvar+exist+i]))
3861 EP = enumerate_cyclic(P, exist, nparam, EP, r, i, MaxRays);
3862 else {
3863 M = Matrix_Alloc(1, nparam+2);
3864 value_set_si(M->p[0][0], 1);
3865 value_set_si(M->p[0][1+i], 1);
3866 enumerate_vd_add_ray(EP, M, MaxRays);
3867 Matrix_Free(M);
3868 }
3869
3870 if (!emptyQ(S)) {
3871 evalue *E = barvinok_enumerate_e(S, exist, nparam, MaxRays);
3872 eadd(E, EP);
3873 free_evalue_refs(E);
3874 free(E);
3875 }
3876 Polyhedron_Free(S);
3877
3878 if (R) {
3879 assert(nvar == 0);
3880 evalue *ER = enumerate_or(R, exist, nparam, MaxRays);
3881 eor(ER, EP);
3882 free_evalue_refs(ER);
3883 free(ER);
3884 }
3885
3886 return EP;
3887 }
3888
3889 static evalue* enumerate_vd(Polyhedron **PA,
3890 unsigned exist, unsigned nparam, unsigned MaxRays)
3891 {
3892 Polyhedron *P = *PA;
3893 int nvar = P->Dimension - exist - nparam;
3894 Param_Polyhedron *PP = NULL;
3895 Polyhedron *C = Universe_Polyhedron(nparam);
3896 Polyhedron *CEq;
3897 Matrix *CT;
3898 Polyhedron *PR = P;
3899 PP = Polyhedron2Param_SimplifiedDomain(&PR,C,MaxRays,&CEq,&CT);
3900 Polyhedron_Free(C);
3901
3902 int nd;
3903 Param_Domain *D, *last;
3904 Value c;
3905 value_init(c);
3906 for (nd = 0, D=PP->D; D; D=D->next, ++nd)
3907 ;
3908
3909 Polyhedron **VD = new Polyhedron_p[nd];
3910 Polyhedron **fVD = new Polyhedron_p[nd];
3911 for(nd = 0, D=PP->D; D; D=D->next) {
3912 Polyhedron *rVD = reduce_domain(D->Domain, CT, CEq,
3913 fVD, nd, MaxRays);
3914 if (!rVD)
3915 continue;
3916
3917 VD[nd++] = rVD;
3918 last = D;
3919 }
3920
3921 evalue *EP = 0;
3922
3923 if (nd == 0)
3924 EP = evalue_zero();
3925
3926 /* This doesn't seem to have any effect */
3927 if (nd == 1) {
3928 Polyhedron *CA = align_context(VD[0], P->Dimension, MaxRays);
3929 Polyhedron *O = P;
3930 P = DomainIntersection(P, CA, MaxRays);
3931 if (O != *PA)
3932 Polyhedron_Free(O);
3933 Polyhedron_Free(CA);
3934 if (emptyQ(P))
3935 EP = evalue_zero();
3936 }
3937
3938 if (!EP && CT->NbColumns != CT->NbRows) {
3939 Polyhedron *CEqr = DomainImage(CEq, CT, MaxRays);
3940 Polyhedron *CA = align_context(CEqr, PR->Dimension, MaxRays);
3941 Polyhedron *I = DomainIntersection(PR, CA, MaxRays);
3942 Polyhedron_Free(CEqr);
3943 Polyhedron_Free(CA);
3944 #ifdef DEBUG_ER
3945 fprintf(stderr, "\nER: Eliminate\n");
3946 #endif /* DEBUG_ER */
3947 nparam -= CT->NbColumns - CT->NbRows;
3948 EP = barvinok_enumerate_e(I, exist, nparam, MaxRays);
3949 nparam += CT->NbColumns - CT->NbRows;
3950 addeliminatedparams_enum(EP, CT, CEq, MaxRays, nparam);
3951 Polyhedron_Free(I);
3952 }
3953 if (PR != *PA)
3954 Polyhedron_Free(PR);
3955 PR = 0;
3956
3957 if (!EP && nd > 1) {
3958 #ifdef DEBUG_ER
3959 fprintf(stderr, "\nER: VD\n");
3960 #endif /* DEBUG_ER */
3961 for (int i = 0; i < nd; ++i) {
3962 Polyhedron *CA = align_context(VD[i], P->Dimension, MaxRays);
3963 Polyhedron *I = DomainIntersection(P, CA, MaxRays);
3964
3965 if (i == 0)
3966 EP = barvinok_enumerate_e(I, exist, nparam, MaxRays);
3967 else {
3968 evalue *E = barvinok_enumerate_e(I, exist, nparam, MaxRays);
3969 eadd(E, EP);
3970 free_evalue_refs(E);
3971 free(E);
3972 }
3973 Polyhedron_Free(I);
3974 Polyhedron_Free(CA);
3975 }
3976 }
3977
3978 for (int i = 0; i < nd; ++i) {
3979 Polyhedron_Free(VD[i]);
3980 Polyhedron_Free(fVD[i]);
3981 }
3982 delete [] VD;
3983 delete [] fVD;
3984 value_clear(c);
3985
3986 if (!EP && nvar == 0) {
3987 Value f;
3988 value_init(f);
3989 Param_Vertices *V, *V2;
3990 Matrix* M = Matrix_Alloc(1, P->Dimension+2);
3991
3992 FORALL_PVertex_in_ParamPolyhedron(V, last, PP) {
3993 bool found = false;
3994 FORALL_PVertex_in_ParamPolyhedron(V2, last, PP) {
3995 if (V == V2) {
3996 found = true;
3997 continue;
3998 }
3999 if (!found)
4000 continue;
4001 for (int i = 0; i < exist; ++i) {
4002 value_oppose(f, V->Vertex->p[i][nparam+1]);
4003 Vector_Combine(V->Vertex->p[i],
4004 V2->Vertex->p[i],
4005 M->p[0] + 1 + nvar + exist,
4006 V2->Vertex->p[i][nparam+1],
4007 f,
4008 nparam+1);
4009 int j;
4010 for (j = 0; j < nparam; ++j)
4011 if (value_notzero_p(M->p[0][1+nvar+exist+j]))
4012 break;
4013 if (j >= nparam)
4014 continue;
4015 ConstraintSimplify(M->p[0], M->p[0],
4016 P->Dimension+2, &f);
4017 value_set_si(M->p[0][0], 0);
4018 Polyhedron *para = AddConstraints(M->p[0], 1, P,
4019 MaxRays);
4020 if (emptyQ(para)) {
4021 Polyhedron_Free(para);
4022 continue;
4023 }
4024 Polyhedron *pos, *neg;
4025 value_set_si(M->p[0][0], 1);
4026 value_decrement(M->p[0][P->Dimension+1],
4027 M->p[0][P->Dimension+1]);
4028 neg = AddConstraints(M->p[0], 1, P, MaxRays);
4029 value_set_si(f, -1);
4030 Vector_Scale(M->p[0]+1, M->p[0]+1, f,
4031 P->Dimension+1);
4032 value_decrement(M->p[0][P->Dimension+1],
4033 M->p[0][P->Dimension+1]);
4034 value_decrement(M->p[0][P->Dimension+1],
4035 M->p[0][P->Dimension+1]);
4036 pos = AddConstraints(M->p[0], 1, P, MaxRays);
4037 if (emptyQ(neg) && emptyQ(pos)) {
4038 Polyhedron_Free(para);
4039 Polyhedron_Free(pos);
4040 Polyhedron_Free(neg);
4041 continue;
4042 }
4043 #ifdef DEBUG_ER
4044 fprintf(stderr, "\nER: Order\n");
4045 #endif /* DEBUG_ER */
4046 EP = barvinok_enumerate_e(para, exist, nparam, MaxRays);
4047 evalue *E;
4048 if (!emptyQ(pos)) {
4049 E = barvinok_enumerate_e(pos, exist, nparam, MaxRays);
4050 eadd(E, EP);
4051 free_evalue_refs(E);
4052 free(E);
4053 }
4054 if (!emptyQ(neg)) {
4055 E = barvinok_enumerate_e(neg, exist, nparam, MaxRays);
4056 eadd(E, EP);
4057 free_evalue_refs(E);
4058 free(E);
4059 }
4060 Polyhedron_Free(para);
4061 Polyhedron_Free(pos);
4062 Polyhedron_Free(neg);
4063 break;
4064 }
4065 if (EP)
4066 break;
4067 } END_FORALL_PVertex_in_ParamPolyhedron;
4068 if (EP)
4069 break;
4070 } END_FORALL_PVertex_in_ParamPolyhedron;
4071
4072 if (!EP) {
4073 /* Search for vertex coordinate to split on */
4074 /* First look for one independent of the parameters */
4075 FORALL_PVertex_in_ParamPolyhedron(V, last, PP) {
4076 for (int i = 0; i < exist; ++i) {
4077 int j;
4078 for (j = 0; j < nparam; ++j)
4079 if (value_notzero_p(V->Vertex->p[i][j]))
4080 break;
4081 if (j < nparam)
4082 continue;
4083 value_set_si(M->p[0][0], 1);
4084 Vector_Set(M->p[0]+1, 0, nvar+exist);
4085 Vector_Copy(V->Vertex->p[i],
4086 M->p[0] + 1 + nvar + exist, nparam+1);
4087 value_oppose(M->p[0][1+nvar+i],
4088 V->Vertex->p[i][nparam+1]);
4089
4090 Polyhedron *pos, *neg;
4091 value_set_si(M->p[0][0], 1);
4092 value_decrement(M->p[0][P->Dimension+1],
4093 M->p[0][P->Dimension+1]);
4094 neg = AddConstraints(M->p[0], 1, P, MaxRays);
4095 value_set_si(f, -1);
4096 Vector_Scale(M->p[0]+1, M->p[0]+1, f,
4097 P->Dimension+1);
4098 value_decrement(M->p[0][P->Dimension+1],
4099 M->p[0][P->Dimension+1]);
4100 value_decrement(M->p[0][P->Dimension+1],
4101 M->p[0][P->Dimension+1]);
4102 pos = AddConstraints(M->p[0], 1, P, MaxRays);
4103 if (emptyQ(neg) || emptyQ(pos)) {
4104 Polyhedron_Free(pos);
4105 Polyhedron_Free(neg);
4106 continue;
4107 }
4108 Polyhedron_Free(pos);
4109 value_increment(M->p[0][P->Dimension+1],
4110 M->p[0][P->Dimension+1]);
4111 pos = AddConstraints(M->p[0], 1, P, MaxRays);
4112 #ifdef DEBUG_ER
4113 fprintf(stderr, "\nER: Vertex\n");
4114 #endif /* DEBUG_ER */
4115 pos->next = neg;
4116 EP = enumerate_or(pos, exist, nparam, MaxRays);
4117 break;
4118 }
4119 if (EP)
4120 break;
4121 } END_FORALL_PVertex_in_ParamPolyhedron;
4122 }
4123
4124 if (!EP) {
4125 /* Search for vertex coordinate to split on */
4126 /* Now look for one that depends on the parameters */
4127 FORALL_PVertex_in_ParamPolyhedron(V, last, PP) {
4128 for (int i = 0; i < exist; ++i) {
4129 value_set_si(M->p[0][0], 1);
4130 Vector_Set(M->p[0]+1, 0, nvar+exist);
4131 Vector_Copy(V->Vertex->p[i],
4132 M->p[0] + 1 + nvar + exist, nparam+1);
4133 value_oppose(M->p[0][1+nvar+i],
4134 V->Vertex->p[i][nparam+1]);
4135
4136 Polyhedron *pos, *neg;
4137 value_set_si(M->p[0][0], 1);
4138 value_decrement(M->p[0][P->Dimension+1],
4139 M->p[0][P->Dimension+1]);
4140 neg = AddConstraints(M->p[0], 1, P, MaxRays);
4141 value_set_si(f, -1);
4142 Vector_Scale(M->p[0]+1, M->p[0]+1, f,
4143 P->Dimension+1);
4144 value_decrement(M->p[0][P->Dimension+1],
4145 M->p[0][P->Dimension+1]);
4146 value_decrement(M->p[0][P->Dimension+1],
4147 M->p[0][P->Dimension+1]);
4148 pos = AddConstraints(M->p[0], 1, P, MaxRays);
4149 if (emptyQ(neg) || emptyQ(pos)) {
4150 Polyhedron_Free(pos);
4151 Polyhedron_Free(neg);
4152 continue;
4153 }
4154 Polyhedron_Free(pos);
4155 value_increment(M->p[0][P->Dimension+1],
4156 M->p[0][P->Dimension+1]);
4157 pos = AddConstraints(M->p[0], 1, P, MaxRays);
4158 #ifdef DEBUG_ER
4159 fprintf(stderr, "\nER: ParamVertex\n");
4160 #endif /* DEBUG_ER */
4161 pos->next = neg;
4162 EP = enumerate_or(pos, exist, nparam, MaxRays);
4163 break;
4164 }
4165 if (EP)
4166 break;
4167 } END_FORALL_PVertex_in_ParamPolyhedron;
4168 }
4169
4170 Matrix_Free(M);
4171 value_clear(f);
4172 }
4173
4174 if (CEq)
4175 Polyhedron_Free(CEq);
4176 if (CT)
4177 Matrix_Free(CT);
4178 if (PP)
4179 Param_Polyhedron_Free(PP);
4180 *PA = P;
4181
4182 return EP;
4183 }
4184
4185 #ifndef HAVE_PIPLIB
4186 evalue *barvinok_enumerate_pip(Polyhedron *P,
4187 unsigned exist, unsigned nparam, unsigned MaxRays)
4188 {
4189 return 0;
4190 }
4191 #else
4192 evalue *barvinok_enumerate_pip(Polyhedron *P,
4193 unsigned exist, unsigned nparam, unsigned MaxRays)
4194 {
4195 int nvar = P->Dimension - exist - nparam;
4196 evalue *EP = evalue_zero();
4197 Polyhedron *Q, *N;
4198
4199 #ifdef DEBUG_ER
4200 fprintf(stderr, "\nER: PIP\n");
4201 #endif /* DEBUG_ER */
4202
4203 Polyhedron *D = pip_projectout(P, nvar, exist, nparam);
4204 for (Q = D; Q; Q = N) {
4205 N = Q->next;
4206 Q->next = 0;
4207 evalue *E;
4208 exist = Q->Dimension - nvar - nparam;
4209 E = barvinok_enumerate_e(Q, exist, nparam, MaxRays);
4210 Polyhedron_Free(Q);
4211 eadd(E, EP);
4212 free_evalue_refs(E);
4213 free(E);
4214 }
4215
4216 return EP;
4217 }
4218 #endif
4219
4220
4221 static bool is_single(Value *row, int pos, int len)
4222 {
4223 return First_Non_Zero(row, pos) == -1 &&
4224 First_Non_Zero(row+pos+1, len-pos-1) == -1;
4225 }
4226
4227 static evalue* barvinok_enumerate_e_r(Polyhedron *P,
4228 unsigned exist, unsigned nparam, unsigned MaxRays);
4229
4230 #ifdef DEBUG_ER
4231 static int er_level = 0;
4232
4233 evalue* barvinok_enumerate_e(Polyhedron *P,
4234 unsigned exist, unsigned nparam, unsigned MaxRays)
4235 {
4236 fprintf(stderr, "\nER: level %i\n", er_level);
4237
4238 Polyhedron_PrintConstraints(stderr, P_VALUE_FMT, P);
4239 fprintf(stderr, "\nE %d\nP %d\n", exist, nparam);
4240 ++er_level;
4241 P = DomainConstraintSimplify(Polyhedron_Copy(P), MaxRays);
4242 evalue *EP = barvinok_enumerate_e_r(P, exist, nparam, MaxRays);
4243 Polyhedron_Free(P);
4244 --er_level;
4245 return EP;
4246 }
4247 #else
4248 evalue* barvinok_enumerate_e(Polyhedron *P,
4249 unsigned exist, unsigned nparam, unsigned MaxRays)
4250 {
4251 P = DomainConstraintSimplify(Polyhedron_Copy(P), MaxRays);
4252 evalue *EP = barvinok_enumerate_e_r(P, exist, nparam, MaxRays);
4253 Polyhedron_Free(P);
4254 return EP;
4255 }
4256 #endif
4257
4258 static evalue* barvinok_enumerate_e_r(Polyhedron *P,
4259 unsigned exist, unsigned nparam, unsigned MaxRays)
4260 {
4261 if (exist == 0) {
4262 Polyhedron *U = Universe_Polyhedron(nparam);
4263 evalue *EP = barvinok_enumerate_ev(P, U, MaxRays);
4264 //char *param_name[] = {"P", "Q", "R", "S", "T" };
4265 //print_evalue(stdout, EP, param_name);
4266 Polyhedron_Free(U);
4267 return EP;
4268 }
4269
4270 int nvar = P->Dimension - exist - nparam;
4271 int len = P->Dimension + 2;
4272
4273 /* for now */
4274 POL_ENSURE_FACETS(P);
4275 POL_ENSURE_VERTICES(P);
4276
4277 if (emptyQ(P))
4278 return evalue_zero();
4279
4280 if (nvar == 0 && nparam == 0) {
4281 evalue *EP = evalue_zero();
4282 barvinok_count(P, &EP->x.n, MaxRays);
4283 if (value_pos_p(EP->x.n))
4284 value_set_si(EP->x.n, 1);
4285 return EP;
4286 }
4287
4288 int r;
4289 for (r = 0; r < P->NbRays; ++r)
4290 if (value_zero_p(P->Ray[r][0]) ||
4291 value_zero_p(P->Ray[r][P->Dimension+1])) {
4292 int i;
4293 for (i = 0; i < nvar; ++i)
4294 if (value_notzero_p(P->Ray[r][i+1]))
4295 break;
4296 if (i >= nvar)
4297 continue;
4298 for (i = nvar + exist; i < nvar + exist + nparam; ++i)
4299 if (value_notzero_p(P->Ray[r][i+1]))
4300 break;
4301 if (i >= nvar + exist + nparam)
4302 break;
4303 }
4304 if (r < P->NbRays) {
4305 evalue *EP = evalue_zero();
4306 value_set_si(EP->x.n, -1);
4307 return EP;
4308 }
4309
4310 int first;
4311 for (r = 0; r < P->NbEq; ++r)
4312 if ((first = First_Non_Zero(P->Constraint[r]+1+nvar, exist)) != -1)
4313 break;
4314 if (r < P->NbEq) {
4315 if (First_Non_Zero(P->Constraint[r]+1+nvar+first+1,
4316 exist-first-1) != -1) {
4317 Polyhedron *T = rotate_along(P, r, nvar, exist, MaxRays);
4318 #ifdef DEBUG_ER
4319 fprintf(stderr, "\nER: Equality\n");
4320 #endif /* DEBUG_ER */
4321 evalue *EP = barvinok_enumerate_e(T, exist-1, nparam, MaxRays);
4322 Polyhedron_Free(T);
4323 return EP;
4324 } else {
4325 #ifdef DEBUG_ER
4326 fprintf(stderr, "\nER: Fixed\n");
4327 #endif /* DEBUG_ER */
4328 if (first == 0)
4329 return barvinok_enumerate_e(P, exist-1, nparam, MaxRays);
4330 else {
4331 Polyhedron *T = Polyhedron_Copy(P);
4332 SwapColumns(T, nvar+1, nvar+1+first);
4333 evalue *EP = barvinok_enumerate_e(T, exist-1, nparam, MaxRays);
4334 Polyhedron_Free(T);
4335 return EP;
4336 }
4337 }
4338 }
4339
4340 Vector *row = Vector_Alloc(len);
4341 value_set_si(row->p[0], 1);
4342
4343 Value f;
4344 value_init(f);
4345
4346 enum constraint* info = new constraint[exist];
4347 for (int i = 0; i < exist; ++i) {
4348 info[i] = ALL_POS;
4349 for (int l = P->NbEq; l < P->NbConstraints; ++l) {
4350 if (value_negz_p(P->Constraint[l][nvar+i+1]))
4351 continue;
4352 bool l_parallel = is_single(P->Constraint[l]+nvar+1, i, exist);
4353 for (int u = P->NbEq; u < P->NbConstraints; ++u) {
4354 if (value_posz_p(P->Constraint[u][nvar+i+1]))
4355 continue;
4356 bool lu_parallel = l_parallel ||
4357 is_single(P->Constraint[u]+nvar+1, i, exist);
4358 value_oppose(f, P->Constraint[u][nvar+i+1]);
4359 Vector_Combine(P->Constraint[l]+1, P->Constraint[u]+1, row->p+1,
4360 f, P->Constraint[l][nvar+i+1], len-1);
4361 if (!(info[i] & INDEPENDENT)) {
4362 int j;
4363 for (j = 0; j < exist; ++j)
4364 if (j != i && value_notzero_p(row->p[nvar+j+1]))
4365 break;
4366 if (j == exist) {
4367 //printf("independent: i: %d, l: %d, u: %d\n", i, l, u);
4368 info[i] = (constraint)(info[i] | INDEPENDENT);
4369 }
4370 }
4371 if (info[i] & ALL_POS) {
4372 value_addto(row->p[len-1], row->p[len-1],
4373 P->Constraint[l][nvar+i+1]);
4374 value_addto(row->p[len-1], row->p[len-1], f);
4375 value_multiply(f, f, P->Constraint[l][nvar+i+1]);
4376 value_subtract(row->p[len-1], row->p[len-1], f);
4377 value_decrement(row->p[len-1], row->p[len-1]);
4378 ConstraintSimplify(row->p, row->p, len, &f);
4379 value_set_si(f, -1);
4380 Vector_Scale(row->p+1, row->p+1, f, len-1);
4381 value_decrement(row->p[len-1], row->p[len-1]);
4382 Polyhedron *T = AddConstraints(row->p, 1, P, MaxRays);
4383 if (!emptyQ(T)) {
4384 //printf("not all_pos: i: %d, l: %d, u: %d\n", i, l, u);
4385 info[i] = (constraint)(info[i] ^ ALL_POS);
4386 }
4387 //puts("pos remainder");
4388 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
4389 Polyhedron_Free(T);
4390 }
4391 if (!(info[i] & ONE_NEG)) {
4392 if (lu_parallel) {
4393 negative_test_constraint(P->Constraint[l],
4394 P->Constraint[u],
4395 row->p, nvar+i, len, &f);
4396 oppose_constraint(row->p, len, &f);
4397 Polyhedron *T = AddConstraints(row->p, 1, P, MaxRays);
4398 if (emptyQ(T)) {
4399 //printf("one_neg i: %d, l: %d, u: %d\n", i, l, u);
4400 info[i] = (constraint)(info[i] | ONE_NEG);
4401 }
4402 //puts("neg remainder");
4403 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
4404 Polyhedron_Free(T);
4405 } else if (!(info[i] & ROT_NEG)) {
4406 if (parallel_constraints(P->Constraint[l],
4407 P->Constraint[u],
4408 row->p, nvar, exist)) {
4409 negative_test_constraint7(P->Constraint[l],
4410 P->Constraint[u],
4411 row->p, nvar, exist,
4412 len, &f);
4413 oppose_constraint(row->p, len, &f);
4414 Polyhedron *T = AddConstraints(row->p, 1, P, MaxRays);
4415 if (emptyQ(T)) {
4416 // printf("rot_neg i: %d, l: %d, u: %d\n", i, l, u);
4417 info[i] = (constraint)(info[i] | ROT_NEG);
4418 r = l;
4419 }
4420 //puts("neg remainder");
4421 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
4422 Polyhedron_Free(T);
4423 }
4424 }
4425 }
4426 if (!(info[i] & ALL_POS) && (info[i] & (ONE_NEG | ROT_NEG)))
4427 goto next;
4428 }
4429 }
4430 if (info[i] & ALL_POS)
4431 break;
4432 next:
4433 ;
4434 }
4435
4436 /*
4437 for (int i = 0; i < exist; ++i)
4438 printf("%i: %i\n", i, info[i]);
4439 */
4440 for (int i = 0; i < exist; ++i)
4441 if (info[i] & ALL_POS) {
4442 #ifdef DEBUG_ER
4443 fprintf(stderr, "\nER: Positive\n");
4444 #endif /* DEBUG_ER */
4445 // Eliminate
4446 // Maybe we should chew off some of the fat here
4447 Matrix *M = Matrix_Alloc(P->Dimension, P->Dimension+1);
4448 for (int j = 0; j < P->Dimension; ++j)
4449 value_set_si(M->p[j][j + (j >= i+nvar)], 1);
4450 Polyhedron *T = Polyhedron_Image(P, M, MaxRays);
4451 Matrix_Free(M);
4452 evalue *EP = barvinok_enumerate_e(T, exist-1, nparam, MaxRays);
4453 Polyhedron_Free(T);
4454 value_clear(f);
4455 Vector_Free(row);
4456 delete [] info;
4457 return EP;
4458 }
4459 for (int i = 0; i < exist; ++i)
4460 if (info[i] & ONE_NEG) {
4461 #ifdef DEBUG_ER
4462 fprintf(stderr, "\nER: Negative\n");
4463 #endif /* DEBUG_ER */
4464 Vector_Free(row);
4465 value_clear(f);
4466 delete [] info;
4467 if (i == 0)
4468 return barvinok_enumerate_e(P, exist-1, nparam, MaxRays);
4469 else {
4470 Polyhedron *T = Polyhedron_Copy(P);
4471 SwapColumns(T, nvar+1, nvar+1+i);
4472 evalue *EP = barvinok_enumerate_e(T, exist-1, nparam, MaxRays);
4473 Polyhedron_Free(T);
4474 return EP;
4475 }
4476 }
4477 for (int i = 0; i < exist; ++i)
4478 if (info[i] & ROT_NEG) {
4479 #ifdef DEBUG_ER
4480 fprintf(stderr, "\nER: Rotate\n");
4481 #endif /* DEBUG_ER */
4482 Vector_Free(row);
4483 value_clear(f);
4484 delete [] info;
4485 Polyhedron *T = rotate_along(P, r, nvar, exist, MaxRays);
4486 evalue *EP = barvinok_enumerate_e(T, exist-1, nparam, MaxRays);
4487 Polyhedron_Free(T);
4488 return EP;
4489 }
4490 for (int i = 0; i < exist; ++i)
4491 if (info[i] & INDEPENDENT) {
4492 Polyhedron *pos, *neg;
4493
4494 /* Find constraint again and split off negative part */
4495
4496 if (SplitOnVar(P, i, nvar, exist, MaxRays,
4497 row, f, true, &pos, &neg)) {
4498 #ifdef DEBUG_ER
4499 fprintf(stderr, "\nER: Split\n");
4500 #endif /* DEBUG_ER */
4501
4502 evalue *EP =
4503 barvinok_enumerate_e(neg, exist-1, nparam, MaxRays);
4504 evalue *E =
4505 barvinok_enumerate_e(pos, exist, nparam, MaxRays);
4506 eadd(E, EP);
4507 free_evalue_refs(E);
4508 free(E);
4509 Polyhedron_Free(neg);
4510 Polyhedron_Free(pos);
4511 value_clear(f);
4512 Vector_Free(row);
4513 delete [] info;
4514 return EP;
4515 }
4516 }
4517 delete [] info;
4518
4519 Polyhedron *O = P;
4520 Polyhedron *F;
4521
4522 evalue *EP;
4523
4524 EP = enumerate_line(P, exist, nparam, MaxRays);
4525 if (EP)
4526 goto out;
4527
4528 EP = barvinok_enumerate_pip(P, exist, nparam, MaxRays);
4529 if (EP)
4530 goto out;
4531
4532 EP = enumerate_redundant_ray(P, exist, nparam, MaxRays);
4533 if (EP)
4534 goto out;
4535
4536 EP = enumerate_sure(P, exist, nparam, MaxRays);
4537 if (EP)
4538 goto out;
4539
4540 EP = enumerate_ray(P, exist, nparam, MaxRays);
4541 if (EP)
4542 goto out;
4543
4544 EP = enumerate_sure2(P, exist, nparam, MaxRays);
4545 if (EP)
4546 goto out;
4547
4548 F = unfringe(P, MaxRays);
4549 if (!PolyhedronIncludes(F, P)) {
4550 #ifdef DEBUG_ER
4551 fprintf(stderr, "\nER: Fringed\n");
4552 #endif /* DEBUG_ER */
4553 EP = barvinok_enumerate_e(F, exist, nparam, MaxRays);
4554 Polyhedron_Free(F);
4555 goto out;
4556 }
4557 Polyhedron_Free(F);
4558
4559 if (nparam)
4560 EP = enumerate_vd(&P, exist, nparam, MaxRays);
4561 if (EP)
4562 goto out2;
4563
4564 if (nvar != 0) {
4565 EP = enumerate_sum(P, exist, nparam, MaxRays);
4566 goto out2;
4567 }
4568
4569 assert(nvar == 0);
4570
4571 int i;
4572 Polyhedron *pos, *neg;
4573 for (i = 0; i < exist; ++i)
4574 if (SplitOnVar(P, i, nvar, exist, MaxRays,
4575 row, f, false, &pos, &neg))
4576 break;
4577
4578 assert (i < exist);
4579
4580 pos->next = neg;
4581 EP = enumerate_or(pos, exist, nparam, MaxRays);
4582
4583 out2:
4584 if (O != P)
4585 Polyhedron_Free(P);
4586
4587 out:
4588 value_clear(f);
4589 Vector_Free(row);
4590 return EP;
4591 }
4592
4593 static void split_param_compression(Matrix *CP, mat_ZZ& map, vec_ZZ& offset)
4594 {
4595 Matrix *T = Transpose(CP);
4596 matrix2zz(T, map, T->NbRows-1, T->NbColumns-1);
4597 values2zz(T->p[T->NbRows-1], offset, T->NbColumns-1);
4598 Matrix_Free(T);
4599 }
4600
4601 /*
4602 * remove equalities that require a "compression" of the parameters
4603 */
4604 #ifndef HAVE_COMPRESS_PARMS
4605 static Polyhedron *remove_more_equalities(Polyhedron *P, unsigned nparam,
4606 Matrix **CP, unsigned MaxRays)
4607 {
4608 return P;
4609 }
4610 #else
4611 static Polyhedron *remove_more_equalities(Polyhedron *P, unsigned nparam,
4612 Matrix **CP, unsigned MaxRays)
4613 {
4614 Matrix *M, *T;
4615 Polyhedron *Q;
4616
4617 /* compress_parms doesn't like equalities that only involve parameters */
4618 for (int i = 0; i < P->NbEq; ++i)
4619 assert(First_Non_Zero(P->Constraint[i]+1, P->Dimension-nparam) != -1);
4620
4621 M = Matrix_Alloc(P->NbEq, P->Dimension+2);
4622 Vector_Copy(P->Constraint[0], M->p[0], P->NbEq * (P->Dimension+2));
4623 *CP = compress_parms(M, nparam);
4624 T = align_matrix(*CP, P->Dimension+1);
4625 Q = Polyhedron_Preimage(P, T, MaxRays);
4626 Polyhedron_Free(P);
4627 P = Q;
4628 P = remove_equalities_p(P, P->Dimension-nparam, NULL);
4629 Matrix_Free(T);
4630 Matrix_Free(M);
4631 return P;
4632 }
4633 #endif
4634
4635 gen_fun * barvinok_series(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
4636 {
4637 Matrix *CP = NULL;
4638 Polyhedron *CA;
4639 unsigned nparam = C->Dimension;
4640
4641 CA = align_context(C, P->Dimension, MaxRays);
4642 P = DomainIntersection(P, CA, MaxRays);
4643 Polyhedron_Free(CA);
4644
4645 if (emptyQ2(P)) {
4646 Polyhedron_Free(P);
4647 return new gen_fun;
4648 }
4649
4650 assert(!Polyhedron_is_infinite(P, nparam));
4651 assert(P->NbBid == 0);
4652 assert(Polyhedron_has_positive_rays(P, nparam));
4653 if (P->NbEq != 0)
4654 P = remove_equalities_p(P, P->Dimension-nparam, NULL);
4655 if (P->NbEq != 0)
4656 P = remove_more_equalities(P, nparam, &CP, MaxRays);
4657 assert(P->NbEq == 0);
4658
4659 #ifdef USE_INCREMENTAL_BF
4660 partial_bfcounter red(P, nparam);
4661 #elif defined USE_INCREMENTAL_DF
4662 partial_ireducer red(P, nparam);
4663 #else
4664 partial_reducer red(P, nparam);
4665 #endif
4666 red.start(MaxRays);
4667 Polyhedron_Free(P);
4668 if (CP) {
4669 mat_ZZ map;
4670 vec_ZZ offset;
4671 split_param_compression(CP, map, offset);
4672 red.gf->substitute(CP, map, offset);
4673 Matrix_Free(CP);
4674 }
4675 return red.gf;
4676 }
4677
4678 static Polyhedron *skew_into_positive_orthant(Polyhedron *D, unsigned nparam,
4679 unsigned MaxRays)
4680 {
4681 Matrix *M = NULL;
4682 Value tmp;
4683 value_init(tmp);
4684 for (Polyhedron *P = D; P; P = P->next) {
4685 POL_ENSURE_VERTICES(P);
4686 assert(!Polyhedron_is_infinite(P, nparam));
4687 assert(P->NbBid == 0);
4688 assert(Polyhedron_has_positive_rays(P, nparam));
4689
4690 for (int r = 0; r < P->NbRays; ++r) {
4691 if (value_notzero_p(P->Ray[r][P->Dimension+1]))
4692 continue;
4693 for (int i = 0; i < nparam; ++i) {
4694 int j;
4695 if (value_posz_p(P->Ray[r][i+1]))
4696 continue;
4697 if (!M) {
4698 M = Matrix_Alloc(D->Dimension+1, D->Dimension+1);
4699 for (int i = 0; i < D->Dimension+1; ++i)
4700 value_set_si(M->p[i][i], 1);
4701 } else {
4702 Inner_Product(P->Ray[r]+1, M->p[i], D->Dimension+1, &tmp);
4703 if (value_posz_p(tmp))
4704 continue;
4705 }
4706 for (j = P->Dimension - nparam; j < P->Dimension; ++j)
4707 if (value_pos_p(P->Ray[r][j+1]))
4708 break;
4709 assert(j < P->Dimension);
4710 value_pdivision(tmp, P->Ray[r][j+1], P->Ray[r][i+1]);
4711 value_subtract(M->p[i][j], M->p[i][j], tmp);
4712 }
4713 }
4714 }
4715 value_clear(tmp);
4716 if (M) {
4717 D = DomainImage(D, M, MaxRays);
4718 Matrix_Free(M);
4719 }
4720 return D;
4721 }
4722
4723 gen_fun* barvinok_enumerate_union_series(Polyhedron *D, Polyhedron* C,
4724 unsigned MaxRays)
4725 {
4726 Polyhedron *conv, *D2;
4727 gen_fun *gf = NULL;
4728 unsigned nparam = C->Dimension;
4729 ZZ one, mone;
4730 one = 1;
4731 mone = -1;
4732 D2 = skew_into_positive_orthant(D, nparam, MaxRays);
4733 for (Polyhedron *P = D2; P; P = P->next) {
4734 assert(P->Dimension == D2->Dimension);
4735 POL_ENSURE_VERTICES(P);
4736 /* it doesn't matter which reducer we use, since we don't actually
4737 * reduce anything here
4738 */
4739 partial_reducer red(P, P->Dimension);
4740 red.start(MaxRays);
4741 if (!gf)
4742 gf = red.gf;
4743 else {
4744 gf->add_union(red.gf, MaxRays);
4745 delete red.gf;
4746 }
4747 }
4748 /* we actually only need the convex union of the parameter space
4749 * but the reducer classes currently expect a polyhedron in
4750 * the combined space
4751 */
4752 conv = DomainConvex(D2, MaxRays);
4753 #ifdef USE_INCREMENTAL_DF
4754 partial_ireducer red(conv, nparam);
4755 #else
4756 partial_reducer red(conv, nparam);
4757 #endif
4758 for (int i = 0; i < gf->term.size(); ++i) {
4759 for (int j = 0; j < gf->term[i]->n.power.NumRows(); ++j) {
4760 red.reduce(gf->term[i]->n.coeff[j][0], gf->term[i]->n.coeff[j][1],
4761 gf->term[i]->n.power[j], gf->term[i]->d.power);
4762 }
4763 }
4764 delete gf;
4765 if (D != D2)
4766 Domain_Free(D2);
4767 Polyhedron_Free(conv);
4768 return red.gf;
4769 }
4770
4771 evalue* barvinok_enumerate_union(Polyhedron *D, Polyhedron* C, unsigned MaxRays)
4772 {
4773 evalue *EP;
4774 gen_fun *gf = barvinok_enumerate_union_series(D, C, MaxRays);
4775 EP = *gf;
4776 delete gf;
4777 return EP;
4778 }
lattice point enumeration library