?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
barvinok.cc: remove_more_equalities: handle equalities only involving params
[barvinok.git] / barvinok.cc
blob763ef2565eab0195e482e98c345366acbe090b38
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 #include "genfun_constructor.h"
24 #include "sample.h"
25
26 #ifdef NTL_STD_CXX
27 using namespace NTL;
28 #endif
29 using std::cerr;
30 using std::cout;
31 using std::endl;
32 using std::vector;
33 using std::deque;
34 using std::string;
35 using std::ostringstream;
36
37 #define ALLOC(t,p) p = (t*)malloc(sizeof(*p))
38
39 static void rays(mat_ZZ& r, Polyhedron *C)
40 {
41 unsigned dim = C->NbRays - 1; /* don't count zero vertex */
42 assert(C->NbRays - 1 == C->Dimension);
43 r.SetDims(dim, dim);
44 ZZ tmp;
45
46 int i, c;
47 for (i = 0, c = 0; i < dim; ++i)
48 if (value_zero_p(C->Ray[i][dim+1])) {
49 for (int j = 0; j < dim; ++j) {
50 value2zz(C->Ray[i][j+1], tmp);
51 r[j][c] = tmp;
52 }
53 ++c;
54 }
55 }
56
57 class dpoly_n {
58 public:
59 Matrix *coeff;
60 ~dpoly_n() {
61 Matrix_Free(coeff);
62 }
63 dpoly_n(int d, ZZ& degree_0, ZZ& degree_1, int offset = 0) {
64 Value d0, d1;
65 value_init(d0);
66 value_init(d1);
67 zz2value(degree_0, d0);
68 zz2value(degree_1, d1);
69 coeff = Matrix_Alloc(d+1, d+1+1);
70 value_set_si(coeff->p[0][0], 1);
71 value_set_si(coeff->p[0][d+1], 1);
72 for (int i = 1; i <= d; ++i) {
73 value_multiply(coeff->p[i][0], coeff->p[i-1][0], d0);
74 Vector_Combine(coeff->p[i-1], coeff->p[i-1]+1, coeff->p[i]+1,
75 d1, d0, i);
76 value_set_si(coeff->p[i][d+1], i);
77 value_multiply(coeff->p[i][d+1], coeff->p[i][d+1], coeff->p[i-1][d+1]);
78 value_decrement(d0, d0);
79 }
80 value_clear(d0);
81 value_clear(d1);
82 }
83 void div(dpoly& d, Vector *count, ZZ& sign) {
84 int len = coeff->NbRows;
85 Matrix * c = Matrix_Alloc(coeff->NbRows, coeff->NbColumns);
86 Value tmp;
87 value_init(tmp);
88 for (int i = 0; i < len; ++i) {
89 Vector_Copy(coeff->p[i], c->p[i], len+1);
90 for (int j = 1; j <= i; ++j) {
91 zz2value(d.coeff[j], tmp);
92 value_multiply(tmp, tmp, c->p[i][len]);
93 value_oppose(tmp, tmp);
94 Vector_Combine(c->p[i], c->p[i-j], c->p[i],
95 c->p[i-j][len], tmp, len);
96 value_multiply(c->p[i][len], c->p[i][len], c->p[i-j][len]);
97 }
98 zz2value(d.coeff[0], tmp);
99 value_multiply(c->p[i][len], c->p[i][len], tmp);
100 }
101 if (sign == -1) {
102 value_set_si(tmp, -1);
103 Vector_Scale(c->p[len-1], count->p, tmp, len);
104 value_assign(count->p[len], c->p[len-1][len]);
105 } else
106 Vector_Copy(c->p[len-1], count->p, len+1);
107 Vector_Normalize(count->p, len+1);
108 value_clear(tmp);
109 Matrix_Free(c);
110 }
111 };
112
113 const int MAX_TRY=10;
114 /*
115 * Searches for a vector that is not orthogonal to any
116 * of the rays in rays.
117 */
118 static void nonorthog(mat_ZZ& rays, vec_ZZ& lambda)
119 {
120 int dim = rays.NumCols();
121 bool found = false;
122 lambda.SetLength(dim);
123 if (dim == 0)
124 return;
125
126 for (int i = 2; !found && i <= 50*dim; i+=4) {
127 for (int j = 0; j < MAX_TRY; ++j) {
128 for (int k = 0; k < dim; ++k) {
129 int r = random_int(i)+2;
130 int v = (2*(r%2)-1) * (r >> 1);
131 lambda[k] = v;
132 }
133 int k = 0;
134 for (; k < rays.NumRows(); ++k)
135 if (lambda * rays[k] == 0)
136 break;
137 if (k == rays.NumRows()) {
138 found = true;
139 break;
140 }
141 }
142 }
143 assert(found);
144 }
145
146 static void add_rays(mat_ZZ& rays, Polyhedron *i, int *r, int nvar = -1,
147 bool all = false)
148 {
149 unsigned dim = i->Dimension;
150 if (nvar == -1)
151 nvar = dim;
152 for (int k = 0; k < i->NbRays; ++k) {
153 if (!value_zero_p(i->Ray[k][dim+1]))
154 continue;
155 if (!all && nvar != dim && First_Non_Zero(i->Ray[k]+1, nvar) == -1)
156 continue;
157 values2zz(i->Ray[k]+1, rays[(*r)++], nvar);
158 }
159 }
160
161 static void mask_r(Matrix *f, int nr, Vector *lcm, int p, Vector *val, evalue *ev)
162 {
163 unsigned nparam = lcm->Size;
164
165 if (p == nparam) {
166 Vector * prod = Vector_Alloc(f->NbRows);
167 Matrix_Vector_Product(f, val->p, prod->p);
168 int isint = 1;
169 for (int i = 0; i < nr; ++i) {
170 value_modulus(prod->p[i], prod->p[i], f->p[i][nparam+1]);
171 isint &= value_zero_p(prod->p[i]);
172 }
173 value_set_si(ev->d, 1);
174 value_init(ev->x.n);
175 value_set_si(ev->x.n, isint);
176 Vector_Free(prod);
177 return;
178 }
179
180 Value tmp;
181 value_init(tmp);
182 if (value_one_p(lcm->p[p]))
183 mask_r(f, nr, lcm, p+1, val, ev);
184 else {
185 value_assign(tmp, lcm->p[p]);
186 value_set_si(ev->d, 0);
187 ev->x.p = new_enode(periodic, VALUE_TO_INT(tmp), p+1);
188 do {
189 value_decrement(tmp, tmp);
190 value_assign(val->p[p], tmp);
191 mask_r(f, nr, lcm, p+1, val, &ev->x.p->arr[VALUE_TO_INT(tmp)]);
192 } while (value_pos_p(tmp));
193 }
194 value_clear(tmp);
195 }
196
197 #ifdef USE_MODULO
198 static void mask(Matrix *f, evalue *factor)
199 {
200 int nr = f->NbRows, nc = f->NbColumns;
201 int n;
202 bool found = false;
203 for (n = 0; n < nr && value_notzero_p(f->p[n][nc-1]); ++n)
204 if (value_notone_p(f->p[n][nc-1]) &&
205 value_notmone_p(f->p[n][nc-1]))
206 found = true;
207 if (!found)
208 return;
209
210 evalue EP;
211 nr = n;
212
213 Value m;
214 value_init(m);
215
216 evalue EV;
217 value_init(EV.d);
218 value_init(EV.x.n);
219 value_set_si(EV.x.n, 1);
220
221 for (n = 0; n < nr; ++n) {
222 value_assign(m, f->p[n][nc-1]);
223 if (value_one_p(m) || value_mone_p(m))
224 continue;
225
226 int j = normal_mod(f->p[n], nc-1, &m);
227 if (j == nc-1) {
228 free_evalue_refs(factor);
229 value_init(factor->d);
230 evalue_set_si(factor, 0, 1);
231 break;
232 }
233 vec_ZZ row;
234 values2zz(f->p[n], row, nc-1);
235 ZZ g;
236 value2zz(m, g);
237 if (j < (nc-1)-1 && row[j] > g/2) {
238 for (int k = j; k < (nc-1); ++k)
239 if (row[k] != 0)
240 row[k] = g - row[k];
241 }
242
243 value_init(EP.d);
244 value_set_si(EP.d, 0);
245 EP.x.p = new_enode(relation, 2, 0);
246 value_clear(EP.x.p->arr[1].d);
247 EP.x.p->arr[1] = *factor;
248 evalue *ev = &EP.x.p->arr[0];
249 value_set_si(ev->d, 0);
250 ev->x.p = new_enode(fractional, 3, -1);
251 evalue_set_si(&ev->x.p->arr[1], 0, 1);
252 evalue_set_si(&ev->x.p->arr[2], 1, 1);
253 evalue *E = multi_monom(row);
254 value_assign(EV.d, m);
255 emul(&EV, E);
256 value_clear(ev->x.p->arr[0].d);
257 ev->x.p->arr[0] = *E;
258 delete E;
259 *factor = EP;
260 }
261
262 value_clear(m);
263 free_evalue_refs(&EV);
264 }
265 #else
266 /*
267 *
268 */
269 static void mask(Matrix *f, evalue *factor)
270 {
271 int nr = f->NbRows, nc = f->NbColumns;
272 int n;
273 bool found = false;
274 for (n = 0; n < nr && value_notzero_p(f->p[n][nc-1]); ++n)
275 if (value_notone_p(f->p[n][nc-1]) &&
276 value_notmone_p(f->p[n][nc-1]))
277 found = true;
278 if (!found)
279 return;
280
281 Value tmp;
282 value_init(tmp);
283 nr = n;
284 unsigned np = nc - 2;
285 Vector *lcm = Vector_Alloc(np);
286 Vector *val = Vector_Alloc(nc);
287 Vector_Set(val->p, 0, nc);
288 value_set_si(val->p[np], 1);
289 Vector_Set(lcm->p, 1, np);
290 for (n = 0; n < nr; ++n) {
291 if (value_one_p(f->p[n][nc-1]) ||
292 value_mone_p(f->p[n][nc-1]))
293 continue;
294 for (int j = 0; j < np; ++j)
295 if (value_notzero_p(f->p[n][j])) {
296 Gcd(f->p[n][j], f->p[n][nc-1], &tmp);
297 value_division(tmp, f->p[n][nc-1], tmp);
298 value_lcm(tmp, lcm->p[j], &lcm->p[j]);
299 }
300 }
301 evalue EP;
302 value_init(EP.d);
303 mask_r(f, nr, lcm, 0, val, &EP);
304 value_clear(tmp);
305 Vector_Free(val);
306 Vector_Free(lcm);
307 emul(&EP,factor);
308 free_evalue_refs(&EP);
309 }
310 #endif
311
312 /* This structure encodes the power of the term in a rational generating function.
313 *
314 * Either E == NULL or constant = 0
315 * If E != NULL, then the power is E
316 * If E == NULL, then the power is coeff * param[pos] + constant
317 */
318 struct term_info {
319 evalue *E;
320 ZZ constant;
321 ZZ coeff;
322 int pos;
323 };
324
325 /* Returns the power of (t+1) in the term of a rational generating function,
326 * i.e., the scalar product of the actual lattice point and lambda.
327 * The lattice point is the unique lattice point in the fundamental parallelepiped
328 * of the unimodual cone i shifted to the parametric vertex V.
329 *
330 * PD is the parameter domain, which, if != NULL, may be used to simply the
331 * resulting expression.
332 *
333 * The result is returned in term.
334 */
335 void lattice_point(
336 Param_Vertices* V, Polyhedron *i, vec_ZZ& lambda, term_info* term,
337 Polyhedron *PD)
338 {
339 unsigned nparam = V->Vertex->NbColumns - 2;
340 unsigned dim = i->Dimension;
341 mat_ZZ vertex;
342 vertex.SetDims(V->Vertex->NbRows, nparam+1);
343 Value lcm, tmp;
344 value_init(lcm);
345 value_init(tmp);
346 value_set_si(lcm, 1);
347 for (int j = 0; j < V->Vertex->NbRows; ++j) {
348 value_lcm(lcm, V->Vertex->p[j][nparam+1], &lcm);
349 }
350 if (value_notone_p(lcm)) {
351 Matrix * mv = Matrix_Alloc(dim, nparam+1);
352 for (int j = 0 ; j < dim; ++j) {
353 value_division(tmp, lcm, V->Vertex->p[j][nparam+1]);
354 Vector_Scale(V->Vertex->p[j], mv->p[j], tmp, nparam+1);
355 }
356
357 term->E = lattice_point(i, lambda, mv, lcm, PD);
358 term->constant = 0;
359
360 Matrix_Free(mv);
361 value_clear(lcm);
362 value_clear(tmp);
363 return;
364 }
365 for (int i = 0; i < V->Vertex->NbRows; ++i) {
366 assert(value_one_p(V->Vertex->p[i][nparam+1])); // for now
367 values2zz(V->Vertex->p[i], vertex[i], nparam+1);
368 }
369
370 vec_ZZ num;
371 num = lambda * vertex;
372
373 int p = -1;
374 int nn = 0;
375 for (int j = 0; j < nparam; ++j)
376 if (num[j] != 0) {
377 ++nn;
378 p = j;
379 }
380 if (nn >= 2) {
381 term->E = multi_monom(num);
382 term->constant = 0;
383 } else {
384 term->E = NULL;
385 term->constant = num[nparam];
386 term->pos = p;
387 if (p != -1)
388 term->coeff = num[p];
389 }
390
391 value_clear(lcm);
392 value_clear(tmp);
393 }
394
395
396 struct counter : public np_base {
397 vec_ZZ lambda;
398 mat_ZZ rays;
399 vec_ZZ vertex;
400 vec_ZZ den;
401 ZZ sign;
402 ZZ num;
403 int j;
404 mpq_t count;
405
406 counter(unsigned dim) : np_base(dim) {
407 rays.SetDims(dim, dim);
408 den.SetLength(dim);
409 mpq_init(count);
410 }
411
412 void start(Polyhedron *P, unsigned MaxRays);
413
414 ~counter() {
415 mpq_clear(count);
416 }
417
418 virtual void handle_polar(Polyhedron *C, Value *vertex, QQ c);
419 };
420
421 struct OrthogonalException {} Orthogonal;
422
423 void counter::handle_polar(Polyhedron *C, Value *V, QQ c)
424 {
425 int r = 0;
426 add_rays(rays, C, &r);
427 for (int k = 0; k < dim; ++k) {
428 if (lambda * rays[k] == 0)
429 throw Orthogonal;
430 }
431
432 assert(c.d == 1);
433 assert(c.n == 1 || c.n == -1);
434 sign = c.n;
435
436 lattice_point(V, C, vertex);
437 num = vertex * lambda;
438 den = rays * lambda;
439 normalize(sign, num, den);
440
441 dpoly d(dim, num);
442 dpoly n(dim, den[0], 1);
443 for (int k = 1; k < dim; ++k) {
444 dpoly fact(dim, den[k], 1);
445 n *= fact;
446 }
447 d.div(n, count, sign);
448 }
449
450 void counter::start(Polyhedron *P, unsigned MaxRays)
451 {
452 for (;;) {
453 try {
454 randomvector(P, lambda, dim);
455 np_base::start(P, MaxRays);
456 break;
457 } catch (OrthogonalException &e) {
458 mpq_set_si(count, 0, 0);
459 }
460 }
461 }
462
463 struct bfe_term : public bfc_term_base {
464 vector<evalue *> factors;
465
466 bfe_term(int len) : bfc_term_base(len) {
467 }
468
469 ~bfe_term() {
470 for (int i = 0; i < factors.size(); ++i) {
471 if (!factors[i])
472 continue;
473 free_evalue_refs(factors[i]);
474 delete factors[i];
475 }
476 }
477 };
478
479 static void print_int_vector(int *v, int len, char *name)
480 {
481 cerr << name << endl;
482 for (int j = 0; j < len; ++j) {
483 cerr << v[j] << " ";
484 }
485 cerr << endl;
486 }
487
488 static void print_bfc_terms(mat_ZZ& factors, bfc_vec& v)
489 {
490 cerr << endl;
491 cerr << "factors" << endl;
492 cerr << factors << endl;
493 for (int i = 0; i < v.size(); ++i) {
494 cerr << "term: " << i << endl;
495 print_int_vector(v[i]->powers, factors.NumRows(), "powers");
496 cerr << "terms" << endl;
497 cerr << v[i]->terms << endl;
498 bfc_term* bfct = static_cast<bfc_term *>(v[i]);
499 cerr << bfct->c << endl;
500 }
501 }
502
503 static void print_bfe_terms(mat_ZZ& factors, bfc_vec& v)
504 {
505 cerr << endl;
506 cerr << "factors" << endl;
507 cerr << factors << endl;
508 for (int i = 0; i < v.size(); ++i) {
509 cerr << "term: " << i << endl;
510 print_int_vector(v[i]->powers, factors.NumRows(), "powers");
511 cerr << "terms" << endl;
512 cerr << v[i]->terms << endl;
513 bfe_term* bfet = static_cast<bfe_term *>(v[i]);
514 for (int j = 0; j < v[i]->terms.NumRows(); ++j) {
515 char * test[] = {"a", "b"};
516 print_evalue(stderr, bfet->factors[j], test);
517 fprintf(stderr, "\n");
518 }
519 }
520 }
521
522 struct bfcounter : public bfcounter_base {
523 mpq_t count;
524
525 bfcounter(unsigned dim) : bfcounter_base(dim) {
526 mpq_init(count);
527 lower = 1;
528 }
529 ~bfcounter() {
530 mpq_clear(count);
531 }
532 virtual void base(mat_ZZ& factors, bfc_vec& v);
533 };
534
535 void bfcounter::base(mat_ZZ& factors, bfc_vec& v)
536 {
537 unsigned nf = factors.NumRows();
538
539 for (int i = 0; i < v.size(); ++i) {
540 bfc_term* bfct = static_cast<bfc_term *>(v[i]);
541 int total_power = 0;
542 // factor is always positive, so we always
543 // change signs
544 for (int k = 0; k < nf; ++k)
545 total_power += v[i]->powers[k];
546
547 int j;
548 for (j = 0; j < nf; ++j)
549 if (v[i]->powers[j] > 0)
550 break;
551
552 dpoly D(total_power, factors[j][0], 1);
553 for (int k = 1; k < v[i]->powers[j]; ++k) {
554 dpoly fact(total_power, factors[j][0], 1);
555 D *= fact;
556 }
557 for ( ; ++j < nf; )
558 for (int k = 0; k < v[i]->powers[j]; ++k) {
559 dpoly fact(total_power, factors[j][0], 1);
560 D *= fact;
561 }
562
563 for (int k = 0; k < v[i]->terms.NumRows(); ++k) {
564 dpoly n(total_power, v[i]->terms[k][0]);
565 mpq_set_si(tcount, 0, 1);
566 n.div(D, tcount, one);
567 if (total_power % 2)
568 bfct->c[k].n = -bfct->c[k].n;
569 zz2value(bfct->c[k].n, tn);
570 zz2value(bfct->c[k].d, td);
571
572 mpz_mul(mpq_numref(tcount), mpq_numref(tcount), tn);
573 mpz_mul(mpq_denref(tcount), mpq_denref(tcount), td);
574 mpq_canonicalize(tcount);
575 mpq_add(count, count, tcount);
576 }
577 delete v[i];
578 }
579 }
580
581
582 /* Check whether the polyhedron is unbounded and if so,
583 * check whether it has any (and therefore an infinite number of)
584 * integer points.
585 * If one of the vertices is integer, then we are done.
586 * Otherwise, transform the polyhedron such that one of the rays
587 * is the first unit vector and cut it off at a height that ensures
588 * that if the whole polyhedron has any points, then the remaining part
589 * has integer points. In particular we add the largest coefficient
590 * of a ray to the highest vertex (rounded up).
591 */
592 static bool Polyhedron_is_infinite(Polyhedron *P, Value* result, unsigned MaxRays)
593 {
594 int r = 0;
595 Matrix *M, *M2;
596 Value c, tmp;
597 Value g;
598 bool first;
599 Vector *v;
600 Value offset, size;
601 Polyhedron *R;
602
603 if (P->NbBid == 0)
604 for (; r < P->NbRays; ++r)
605 if (value_zero_p(P->Ray[r][P->Dimension+1]))
606 break;
607 if (P->NbBid == 0 && r == P->NbRays)
608 return false;
609
610 #ifdef HAVE_LIBGLPK
611 Vector *sample;
612
613 sample = Polyhedron_Sample(P, MaxRays);
614 if (!sample)
615 value_set_si(*result, 0);
616 else {
617 value_set_si(*result, -1);
618 Vector_Free(sample);
619 }
620 return true;
621 #endif
622
623 for (int i = 0; i < P->NbRays; ++i)
624 if (value_one_p(P->Ray[i][1+P->Dimension])) {
625 value_set_si(*result, -1);
626 return true;
627 }
628
629 value_init(g);
630 v = Vector_Alloc(P->Dimension+1);
631 Vector_Gcd(P->Ray[r]+1, P->Dimension, &g);
632 Vector_AntiScale(P->Ray[r]+1, v->p, g, P->Dimension+1);
633 M = unimodular_complete(v);
634 value_set_si(M->p[P->Dimension][P->Dimension], 1);
635 M2 = Transpose(M);
636 Matrix_Free(M);
637 P = Polyhedron_Preimage(P, M2, 0);
638 Matrix_Free(M2);
639 value_clear(g);
640 Vector_Free(v);
641
642 first = true;
643 value_init(offset);
644 value_init(size);
645 value_init(tmp);
646 value_set_si(size, 0);
647
648 for (int i = 0; i < P->NbBid; ++i) {
649 value_absolute(tmp, P->Ray[i][1]);
650 if (value_gt(tmp, size))
651 value_assign(size, tmp);
652 }
653 for (int i = P->NbBid; i < P->NbRays; ++i) {
654 if (value_zero_p(P->Ray[i][P->Dimension+1])) {
655 if (value_gt(P->Ray[i][1], size))
656 value_assign(size, P->Ray[i][1]);
657 continue;
658 }
659 mpz_cdiv_q(tmp, P->Ray[i][1], P->Ray[i][P->Dimension+1]);
660 if (first || value_gt(tmp, offset)) {
661 value_assign(offset, tmp);
662 first = false;
663 }
664 }
665 value_addto(offset, offset, size);
666 value_clear(size);
667 value_clear(tmp);
668
669 v = Vector_Alloc(P->Dimension+2);
670 value_set_si(v->p[0], 1);
671 value_set_si(v->p[1], -1);
672 value_assign(v->p[1+P->Dimension], offset);
673 R = AddConstraints(v->p, 1, P, MaxRays);
674 Polyhedron_Free(P);
675 P = R;
676
677 value_clear(offset);
678 Vector_Free(v);
679
680 value_init(c);
681 barvinok_count(P, &c, MaxRays);
682 Polyhedron_Free(P);
683 if (value_zero_p(c))
684 value_set_si(*result, 0);
685 else
686 value_set_si(*result, -1);
687 value_clear(c);
688
689 return true;
690 }
691
692 typedef Polyhedron * Polyhedron_p;
693
694 static void barvinok_count_f(Polyhedron *P, Value* result, unsigned NbMaxCons);
695
696 void barvinok_count(Polyhedron *P, Value* result, unsigned NbMaxCons)
697 {
698 unsigned dim;
699 int allocated = 0;
700 Polyhedron *Q;
701 bool infinite = false;
702
703 if (emptyQ2(P)) {
704 value_set_si(*result, 0);
705 return;
706 }
707 if (P->NbEq != 0) {
708 do {
709 P = remove_equalities(P);
710 P = DomainConstraintSimplify(P, NbMaxCons);
711 } while (!emptyQ(P) && P->NbEq != 0);
712 if (emptyQ(P)) {
713 Polyhedron_Free(P);
714 value_set_si(*result, 0);
715 return;
716 }
717 allocated = 1;
718 }
719 if (Polyhedron_is_infinite(P, result, NbMaxCons)) {
720 if (allocated)
721 Polyhedron_Free(P);
722 return;
723 }
724 if (P->Dimension == 0) {
725 /* Test whether the constraints are satisfied */
726 POL_ENSURE_VERTICES(P);
727 value_set_si(*result, !emptyQ(P));
728 if (allocated)
729 Polyhedron_Free(P);
730 return;
731 }
732 Q = Polyhedron_Factor(P, 0, NbMaxCons);
733 if (Q) {
734 if (allocated)
735 Polyhedron_Free(P);
736 P = Q;
737 allocated = 1;
738 }
739
740 barvinok_count_f(P, result, NbMaxCons);
741 if (value_neg_p(*result))
742 infinite = true;
743 if (Q && P->next && value_notzero_p(*result)) {
744 Value factor;
745 value_init(factor);
746
747 for (Q = P->next; Q; Q = Q->next) {
748 barvinok_count_f(Q, &factor, NbMaxCons);
749 if (value_neg_p(factor)) {
750 infinite = true;
751 continue;
752 } else if (Q->next && value_zero_p(factor)) {
753 value_set_si(*result, 0);
754 break;
755 }
756 value_multiply(*result, *result, factor);
757 }
758
759 value_clear(factor);
760 }
761
762 if (allocated)
763 Domain_Free(P);
764 if (infinite)
765 value_set_si(*result, -1);
766 }
767
768 static void barvinok_count_f(Polyhedron *P, Value* result, unsigned NbMaxCons)
769 {
770 if (emptyQ2(P)) {
771 value_set_si(*result, 0);
772 return;
773 }
774
775 if (P->Dimension == 1)
776 return Line_Length(P, result);
777
778 int c = P->NbConstraints;
779 POL_ENSURE_FACETS(P);
780 if (c != P->NbConstraints || P->NbEq != 0)
781 return barvinok_count(P, result, NbMaxCons);
782
783 POL_ENSURE_VERTICES(P);
784
785 if (Polyhedron_is_infinite(P, result, NbMaxCons))
786 return;
787
788 #ifdef USE_INCREMENTAL_BF
789 bfcounter cnt(P->Dimension);
790 #elif defined USE_INCREMENTAL_DF
791 icounter cnt(P->Dimension);
792 #else
793 counter cnt(P->Dimension);
794 #endif
795 cnt.start(P, NbMaxCons);
796
797 assert(value_one_p(&cnt.count[0]._mp_den));
798 value_assign(*result, &cnt.count[0]._mp_num);
799 }
800
801 static void uni_polynom(int param, Vector *c, evalue *EP)
802 {
803 unsigned dim = c->Size-2;
804 value_init(EP->d);
805 value_set_si(EP->d,0);
806 EP->x.p = new_enode(polynomial, dim+1, param+1);
807 for (int j = 0; j <= dim; ++j)
808 evalue_set(&EP->x.p->arr[j], c->p[j], c->p[dim+1]);
809 }
810
811 static void multi_polynom(Vector *c, evalue* X, evalue *EP)
812 {
813 unsigned dim = c->Size-2;
814 evalue EC;
815
816 value_init(EC.d);
817 evalue_set(&EC, c->p[dim], c->p[dim+1]);
818
819 value_init(EP->d);
820 evalue_set(EP, c->p[dim], c->p[dim+1]);
821
822 for (int i = dim-1; i >= 0; --i) {
823 emul(X, EP);
824 value_assign(EC.x.n, c->p[i]);
825 eadd(&EC, EP);
826 }
827 free_evalue_refs(&EC);
828 }
829
830 Polyhedron *unfringe (Polyhedron *P, unsigned MaxRays)
831 {
832 int len = P->Dimension+2;
833 Polyhedron *T, *R = P;
834 Value g;
835 value_init(g);
836 Vector *row = Vector_Alloc(len);
837 value_set_si(row->p[0], 1);
838
839 R = DomainConstraintSimplify(Polyhedron_Copy(P), MaxRays);
840
841 Matrix *M = Matrix_Alloc(2, len-1);
842 value_set_si(M->p[1][len-2], 1);
843 for (int v = 0; v < P->Dimension; ++v) {
844 value_set_si(M->p[0][v], 1);
845 Polyhedron *I = Polyhedron_Image(R, M, 2+1);
846 value_set_si(M->p[0][v], 0);
847 for (int r = 0; r < I->NbConstraints; ++r) {
848 if (value_zero_p(I->Constraint[r][0]))
849 continue;
850 if (value_zero_p(I->Constraint[r][1]))
851 continue;
852 if (value_one_p(I->Constraint[r][1]))
853 continue;
854 if (value_mone_p(I->Constraint[r][1]))
855 continue;
856 value_absolute(g, I->Constraint[r][1]);
857 Vector_Set(row->p+1, 0, len-2);
858 value_division(row->p[1+v], I->Constraint[r][1], g);
859 mpz_fdiv_q(row->p[len-1], I->Constraint[r][2], g);
860 T = R;
861 R = AddConstraints(row->p, 1, R, MaxRays);
862 if (T != P)
863 Polyhedron_Free(T);
864 }
865 Polyhedron_Free(I);
866 }
867 Matrix_Free(M);
868 Vector_Free(row);
869 value_clear(g);
870 return R;
871 }
872
873 /* this procedure may have false negatives */
874 static bool Polyhedron_is_infinite_param(Polyhedron *P, unsigned nparam)
875 {
876 int r;
877 for (r = 0; r < P->NbRays; ++r) {
878 if (!value_zero_p(P->Ray[r][0]) &&
879 !value_zero_p(P->Ray[r][P->Dimension+1]))
880 continue;
881 if (First_Non_Zero(P->Ray[r]+1+P->Dimension-nparam, nparam) == -1)
882 return true;
883 }
884 return false;
885 }
886
887 /* Check whether all rays point in the positive directions
888 * for the parameters
889 */
890 static bool Polyhedron_has_positive_rays(Polyhedron *P, unsigned nparam)
891 {
892 int r;
893 for (r = 0; r < P->NbRays; ++r)
894 if (value_zero_p(P->Ray[r][P->Dimension+1])) {
895 int i;
896 for (i = P->Dimension - nparam; i < P->Dimension; ++i)
897 if (value_neg_p(P->Ray[r][i+1]))
898 return false;
899 }
900 return true;
901 }
902
903 typedef evalue * evalue_p;
904
905 struct enumerator : public polar_decomposer {
906 vec_ZZ lambda;
907 unsigned dim, nbV;
908 evalue ** vE;
909 int _i;
910 mat_ZZ rays;
911 vec_ZZ den;
912 ZZ sign;
913 Polyhedron *P;
914 Param_Vertices *V;
915 term_info num;
916 Vector *c;
917 mpq_t count;
918
919 enumerator(Polyhedron *P, unsigned dim, unsigned nbV) {
920 this->P = P;
921 this->dim = dim;
922 this->nbV = nbV;
923 randomvector(P, lambda, dim);
924 rays.SetDims(dim, dim);
925 den.SetLength(dim);
926 c = Vector_Alloc(dim+2);
927
928 vE = new evalue_p[nbV];
929 for (int j = 0; j < nbV; ++j)
930 vE[j] = 0;
931
932 mpq_init(count);
933 }
934
935 void decompose_at(Param_Vertices *V, int _i, unsigned MaxRays) {
936 Polyhedron *C = supporting_cone_p(P, V);
937 this->_i = _i;
938 this->V = V;
939
940 vE[_i] = new evalue;
941 value_init(vE[_i]->d);
942 evalue_set_si(vE[_i], 0, 1);
943
944 decompose(C, MaxRays);
945 }
946
947 ~enumerator() {
948 mpq_clear(count);
949 Vector_Free(c);
950
951 for (int j = 0; j < nbV; ++j)
952 if (vE[j]) {
953 free_evalue_refs(vE[j]);
954 delete vE[j];
955 }
956 delete [] vE;
957 }
958
959 virtual void handle_polar(Polyhedron *P, int sign);
960 };
961
962 void enumerator::handle_polar(Polyhedron *C, int s)
963 {
964 int r = 0;
965 assert(C->NbRays-1 == dim);
966 add_rays(rays, C, &r);
967 for (int k = 0; k < dim; ++k) {
968 if (lambda * rays[k] == 0)
969 throw Orthogonal;
970 }
971
972 sign = s;
973
974 lattice_point(V, C, lambda, &num, 0);
975 den = rays * lambda;
976 normalize(sign, num.constant, den);
977
978 dpoly n(dim, den[0], 1);
979 for (int k = 1; k < dim; ++k) {
980 dpoly fact(dim, den[k], 1);
981 n *= fact;
982 }
983 if (num.E != NULL) {
984 ZZ one(INIT_VAL, 1);
985 dpoly_n d(dim, num.constant, one);
986 d.div(n, c, sign);
987 evalue EV;
988 multi_polynom(c, num.E, &EV);
989 eadd(&EV , vE[_i]);
990 free_evalue_refs(&EV);
991 free_evalue_refs(num.E);
992 delete num.E;
993 } else if (num.pos != -1) {
994 dpoly_n d(dim, num.constant, num.coeff);
995 d.div(n, c, sign);
996 evalue EV;
997 uni_polynom(num.pos, c, &EV);
998 eadd(&EV , vE[_i]);
999 free_evalue_refs(&EV);
1000 } else {
1001 mpq_set_si(count, 0, 1);
1002 dpoly d(dim, num.constant);
1003 d.div(n, count, sign);
1004 evalue EV;
1005 value_init(EV.d);
1006 evalue_set(&EV, &count[0]._mp_num, &count[0]._mp_den);
1007 eadd(&EV , vE[_i]);
1008 free_evalue_refs(&EV);
1009 }
1010 }
1011
1012 struct enumerator_base {
1013 unsigned dim;
1014 evalue ** vE;
1015 evalue ** E_vertex;
1016 evalue mone;
1017 vertex_decomposer *vpd;
1018
1019 enumerator_base(unsigned dim, vertex_decomposer *vpd)
1020 {
1021 this->dim = dim;
1022 this->vpd = vpd;
1023
1024 vE = new evalue_p[vpd->nbV];
1025 for (int j = 0; j < vpd->nbV; ++j)
1026 vE[j] = 0;
1027
1028 E_vertex = new evalue_p[dim];
1029
1030 value_init(mone.d);
1031 evalue_set_si(&mone, -1, 1);
1032 }
1033
1034 void decompose_at(Param_Vertices *V, int _i, unsigned MaxRays/*, Polyhedron *pVD*/) {
1035 //this->pVD = pVD;
1036
1037 vE[_i] = new evalue;
1038 value_init(vE[_i]->d);
1039 evalue_set_si(vE[_i], 0, 1);
1040
1041 vpd->decompose_at_vertex(V, _i, MaxRays);
1042 }
1043
1044 ~enumerator_base() {
1045 for (int j = 0; j < vpd->nbV; ++j)
1046 if (vE[j]) {
1047 free_evalue_refs(vE[j]);
1048 delete vE[j];
1049 }
1050 delete [] vE;
1051
1052 delete [] E_vertex;
1053
1054 free_evalue_refs(&mone);
1055 }
1056
1057 evalue *E_num(int i, int d) {
1058 return E_vertex[i + (dim-d)];
1059 }
1060 };
1061
1062 struct cumulator {
1063 evalue *factor;
1064 evalue *v;
1065 dpoly_r *r;
1066
1067 cumulator(evalue *factor, evalue *v, dpoly_r *r) :
1068 factor(factor), v(v), r(r) {}
1069
1070 void cumulate();
1071
1072 virtual void add_term(int *powers, int len, evalue *f2) = 0;
1073 };
1074
1075 void cumulator::cumulate()
1076 {
1077 evalue cum; // factor * 1 * E_num[0]/1 * (E_num[0]-1)/2 *...
1078 evalue f;
1079 evalue t; // E_num[0] - (m-1)
1080 #ifdef USE_MODULO
1081 evalue *cst;
1082 #else
1083 evalue mone;
1084 value_init(mone.d);
1085 evalue_set_si(&mone, -1, 1);
1086 #endif
1087
1088 value_init(cum.d);
1089 evalue_copy(&cum, factor);
1090 value_init(f.d);
1091 value_init(f.x.n);
1092 value_set_si(f.d, 1);
1093 value_set_si(f.x.n, 1);
1094 value_init(t.d);
1095 evalue_copy(&t, v);
1096
1097 #ifdef USE_MODULO
1098 for (cst = &t; value_zero_p(cst->d); ) {
1099 if (cst->x.p->type == fractional)
1100 cst = &cst->x.p->arr[1];
1101 else
1102 cst = &cst->x.p->arr[0];
1103 }
1104 #endif
1105
1106 for (int m = 0; m < r->len; ++m) {
1107 if (m > 0) {
1108 if (m > 1) {
1109 value_set_si(f.d, m);
1110 emul(&f, &cum);
1111 #ifdef USE_MODULO
1112 value_subtract(cst->x.n, cst->x.n, cst->d);
1113 #else
1114 eadd(&mone, &t);
1115 #endif
1116 }
1117 emul(&t, &cum);
1118 }
1119 vector< dpoly_r_term * >& current = r->c[r->len-1-m];
1120 for (int j = 0; j < current.size(); ++j) {
1121 if (current[j]->coeff == 0)
1122 continue;
1123 evalue *f2 = new evalue;
1124 value_init(f2->d);
1125 value_init(f2->x.n);
1126 zz2value(current[j]->coeff, f2->x.n);
1127 zz2value(r->denom, f2->d);
1128 emul(&cum, f2);
1129
1130 add_term(current[j]->powers, r->dim, f2);
1131 }
1132 }
1133 free_evalue_refs(&f);
1134 free_evalue_refs(&t);
1135 free_evalue_refs(&cum);
1136 #ifndef USE_MODULO
1137 free_evalue_refs(&mone);
1138 #endif
1139 }
1140
1141 struct E_poly_term {
1142 int *powers;
1143 evalue *E;
1144 };
1145
1146 struct ie_cum : public cumulator {
1147 vector<E_poly_term *> terms;
1148
1149 ie_cum(evalue *factor, evalue *v, dpoly_r *r) : cumulator(factor, v, r) {}
1150
1151 virtual void add_term(int *powers, int len, evalue *f2);
1152 };
1153
1154 void ie_cum::add_term(int *powers, int len, evalue *f2)
1155 {
1156 int k;
1157 for (k = 0; k < terms.size(); ++k) {
1158 if (memcmp(terms[k]->powers, powers, len * sizeof(int)) == 0) {
1159 eadd(f2, terms[k]->E);
1160 free_evalue_refs(f2);
1161 delete f2;
1162 break;
1163 }
1164 }
1165 if (k >= terms.size()) {
1166 E_poly_term *ET = new E_poly_term;
1167 ET->powers = new int[len];
1168 memcpy(ET->powers, powers, len * sizeof(int));
1169 ET->E = f2;
1170 terms.push_back(ET);
1171 }
1172 }
1173
1174 struct ienumerator : public polar_decomposer, public vertex_decomposer,
1175 public enumerator_base {
1176 //Polyhedron *pVD;
1177 mat_ZZ den;
1178 vec_ZZ vertex;
1179 mpq_t tcount;
1180
1181 ienumerator(Polyhedron *P, unsigned dim, unsigned nbV) :
1182 vertex_decomposer(P, nbV, this), enumerator_base(dim, this) {
1183 vertex.SetLength(dim);
1184
1185 den.SetDims(dim, dim);
1186 mpq_init(tcount);
1187 }
1188
1189 ~ienumerator() {
1190 mpq_clear(tcount);
1191 }
1192
1193 virtual void handle_polar(Polyhedron *P, int sign);
1194 void reduce(evalue *factor, vec_ZZ& num, mat_ZZ& den_f);
1195 };
1196
1197 void ienumerator::reduce(
1198 evalue *factor, vec_ZZ& num, mat_ZZ& den_f)
1199 {
1200 unsigned len = den_f.NumRows(); // number of factors in den
1201 unsigned dim = num.length();
1202
1203 if (dim == 0) {
1204 eadd(factor, vE[vert]);
1205 return;
1206 }
1207
1208 vec_ZZ den_s;
1209 den_s.SetLength(len);
1210 mat_ZZ den_r;
1211 den_r.SetDims(len, dim-1);
1212
1213 int r, k;
1214
1215 for (r = 0; r < len; ++r) {
1216 den_s[r] = den_f[r][0];
1217 for (k = 0; k <= dim-1; ++k)
1218 if (k != 0)
1219 den_r[r][k-(k>0)] = den_f[r][k];
1220 }
1221
1222 ZZ num_s = num[0];
1223 vec_ZZ num_p;
1224 num_p.SetLength(dim-1);
1225 for (k = 0 ; k <= dim-1; ++k)
1226 if (k != 0)
1227 num_p[k-(k>0)] = num[k];
1228
1229 vec_ZZ den_p;
1230 den_p.SetLength(len);
1231
1232 ZZ one;
1233 one = 1;
1234 normalize(one, num_s, num_p, den_s, den_p, den_r);
1235 if (one != 1)
1236 emul(&mone, factor);
1237
1238 int only_param = 0;
1239 int no_param = 0;
1240 for (int k = 0; k < len; ++k) {
1241 if (den_p[k] == 0)
1242 ++no_param;
1243 else if (den_s[k] == 0)
1244 ++only_param;
1245 }
1246 if (no_param == 0) {
1247 reduce(factor, num_p, den_r);
1248 } else {
1249 int k, l;
1250 mat_ZZ pden;
1251 pden.SetDims(only_param, dim-1);
1252
1253 for (k = 0, l = 0; k < len; ++k)
1254 if (den_s[k] == 0)
1255 pden[l++] = den_r[k];
1256
1257 for (k = 0; k < len; ++k)
1258 if (den_p[k] == 0)
1259 break;
1260
1261 dpoly n(no_param, num_s);
1262 dpoly D(no_param, den_s[k], 1);
1263 for ( ; ++k < len; )
1264 if (den_p[k] == 0) {
1265 dpoly fact(no_param, den_s[k], 1);
1266 D *= fact;
1267 }
1268
1269 dpoly_r * r = 0;
1270 // if no_param + only_param == len then all powers
1271 // below will be all zero
1272 if (no_param + only_param == len) {
1273 if (E_num(0, dim) != 0)
1274 r = new dpoly_r(n, len);
1275 else {
1276 mpq_set_si(tcount, 0, 1);
1277 one = 1;
1278 n.div(D, tcount, one);
1279
1280 if (value_notzero_p(mpq_numref(tcount))) {
1281 evalue f;
1282 value_init(f.d);
1283 value_init(f.x.n);
1284 value_assign(f.x.n, mpq_numref(tcount));
1285 value_assign(f.d, mpq_denref(tcount));
1286 emul(&f, factor);
1287 reduce(factor, num_p, pden);
1288 free_evalue_refs(&f);
1289 }
1290 return;
1291 }
1292 } else {
1293 for (k = 0; k < len; ++k) {
1294 if (den_s[k] == 0 || den_p[k] == 0)
1295 continue;
1296
1297 dpoly pd(no_param-1, den_s[k], 1);
1298
1299 int l;
1300 for (l = 0; l < k; ++l)
1301 if (den_r[l] == den_r[k])
1302 break;
1303
1304 if (r == 0)
1305 r = new dpoly_r(n, pd, l, len);
1306 else {
1307 dpoly_r *nr = new dpoly_r(r, pd, l, len);
1308 delete r;
1309 r = nr;
1310 }
1311 }
1312 }
1313 dpoly_r *rc = r->div(D);
1314 delete r;
1315 r = rc;
1316 if (E_num(0, dim) == 0) {
1317 int common = pden.NumRows();
1318 vector< dpoly_r_term * >& final = r->c[r->len-1];
1319 int rows;
1320 evalue t;
1321 evalue f;
1322 value_init(f.d);
1323 value_init(f.x.n);
1324 zz2value(r->denom, f.d);
1325 for (int j = 0; j < final.size(); ++j) {
1326 if (final[j]->coeff == 0)
1327 continue;
1328 rows = common;
1329 for (int k = 0; k < r->dim; ++k) {
1330 int n = final[j]->powers[k];
1331 if (n == 0)
1332 continue;
1333 pden.SetDims(rows+n, pden.NumCols());
1334 for (int l = 0; l < n; ++l)
1335 pden[rows+l] = den_r[k];
1336 rows += n;
1337 }
1338 value_init(t.d);
1339 evalue_copy(&t, factor);
1340 zz2value(final[j]->coeff, f.x.n);
1341 emul(&f, &t);
1342 reduce(&t, num_p, pden);
1343 free_evalue_refs(&t);
1344 }
1345 free_evalue_refs(&f);
1346 } else {
1347 ie_cum cum(factor, E_num(0, dim), r);
1348 cum.cumulate();
1349
1350 int common = pden.NumRows();
1351 int rows;
1352 for (int j = 0; j < cum.terms.size(); ++j) {
1353 rows = common;
1354 pden.SetDims(rows, pden.NumCols());
1355 for (int k = 0; k < r->dim; ++k) {
1356 int n = cum.terms[j]->powers[k];
1357 if (n == 0)
1358 continue;
1359 pden.SetDims(rows+n, pden.NumCols());
1360 for (int l = 0; l < n; ++l)
1361 pden[rows+l] = den_r[k];
1362 rows += n;
1363 }
1364 reduce(cum.terms[j]->E, num_p, pden);
1365 free_evalue_refs(cum.terms[j]->E);
1366 delete cum.terms[j]->E;
1367 delete [] cum.terms[j]->powers;
1368 delete cum.terms[j];
1369 }
1370 }
1371 delete r;
1372 }
1373 }
1374
1375 static int type_offset(enode *p)
1376 {
1377 return p->type == fractional ? 1 :
1378 p->type == flooring ? 1 : 0;
1379 }
1380
1381 static int edegree(evalue *e)
1382 {
1383 int d = 0;
1384 enode *p;
1385
1386 if (value_notzero_p(e->d))
1387 return 0;
1388
1389 p = e->x.p;
1390 int i = type_offset(p);
1391 if (p->size-i-1 > d)
1392 d = p->size - i - 1;
1393 for (; i < p->size; i++) {
1394 int d2 = edegree(&p->arr[i]);
1395 if (d2 > d)
1396 d = d2;
1397 }
1398 return d;
1399 }
1400
1401 void ienumerator::handle_polar(Polyhedron *C, int s)
1402 {
1403 assert(C->NbRays-1 == dim);
1404
1405 lattice_point(V, C, vertex, E_vertex);
1406
1407 int r;
1408 for (r = 0; r < dim; ++r)
1409 values2zz(C->Ray[r]+1, den[r], dim);
1410
1411 evalue one;
1412 value_init(one.d);
1413 evalue_set_si(&one, s, 1);
1414 reduce(&one, vertex, den);
1415 free_evalue_refs(&one);
1416
1417 for (int i = 0; i < dim; ++i)
1418 if (E_vertex[i]) {
1419 free_evalue_refs(E_vertex[i]);
1420 delete E_vertex[i];
1421 }
1422 }
1423
1424 struct bfenumerator : public vertex_decomposer, public bf_base,
1425 public enumerator_base {
1426 evalue *factor;
1427
1428 bfenumerator(Polyhedron *P, unsigned dim, unsigned nbV) :
1429 vertex_decomposer(P, nbV, this),
1430 bf_base(dim), enumerator_base(dim, this) {
1431 lower = 0;
1432 factor = NULL;
1433 }
1434
1435 ~bfenumerator() {
1436 }
1437
1438 virtual void handle_polar(Polyhedron *P, int sign);
1439 virtual void base(mat_ZZ& factors, bfc_vec& v);
1440
1441 bfc_term_base* new_bf_term(int len) {
1442 bfe_term* t = new bfe_term(len);
1443 return t;
1444 }
1445
1446 virtual void set_factor(bfc_term_base *t, int k, int change) {
1447 bfe_term* bfet = static_cast<bfe_term *>(t);
1448 factor = bfet->factors[k];
1449 assert(factor != NULL);
1450 bfet->factors[k] = NULL;
1451 if (change)
1452 emul(&mone, factor);
1453 }
1454
1455 virtual void set_factor(bfc_term_base *t, int k, mpq_t &q, int change) {
1456 bfe_term* bfet = static_cast<bfe_term *>(t);
1457 factor = bfet->factors[k];
1458 assert(factor != NULL);
1459 bfet->factors[k] = NULL;
1460
1461 evalue f;
1462 value_init(f.d);
1463 value_init(f.x.n);
1464 if (change)
1465 value_oppose(f.x.n, mpq_numref(q));
1466 else
1467 value_assign(f.x.n, mpq_numref(q));
1468 value_assign(f.d, mpq_denref(q));
1469 emul(&f, factor);
1470 free_evalue_refs(&f);
1471 }
1472
1473 virtual void set_factor(bfc_term_base *t, int k, const QQ& c, int change) {
1474 bfe_term* bfet = static_cast<bfe_term *>(t);
1475
1476 factor = new evalue;
1477
1478 evalue f;
1479 value_init(f.d);
1480 value_init(f.x.n);
1481 zz2value(c.n, f.x.n);
1482 if (change)
1483 value_oppose(f.x.n, f.x.n);
1484 zz2value(c.d, f.d);
1485
1486 value_init(factor->d);
1487 evalue_copy(factor, bfet->factors[k]);
1488 emul(&f, factor);
1489 free_evalue_refs(&f);
1490 }
1491
1492 void set_factor(evalue *f, int change) {
1493 if (change)
1494 emul(&mone, f);
1495 factor = f;
1496 }
1497
1498 virtual void insert_term(bfc_term_base *t, int i) {
1499 bfe_term* bfet = static_cast<bfe_term *>(t);
1500 int len = t->terms.NumRows()-1; // already increased by one
1501
1502 bfet->factors.resize(len+1);
1503 for (int j = len; j > i; --j) {
1504 bfet->factors[j] = bfet->factors[j-1];
1505 t->terms[j] = t->terms[j-1];
1506 }
1507 bfet->factors[i] = factor;
1508 factor = NULL;
1509 }
1510
1511 virtual void update_term(bfc_term_base *t, int i) {
1512 bfe_term* bfet = static_cast<bfe_term *>(t);
1513
1514 eadd(factor, bfet->factors[i]);
1515 free_evalue_refs(factor);
1516 delete factor;
1517 }
1518
1519 virtual bool constant_vertex(int dim) { return E_num(0, dim) == 0; }
1520
1521 virtual void cum(bf_reducer *bfr, bfc_term_base *t, int k, dpoly_r *r);
1522 };
1523
1524 struct bfe_cum : public cumulator {
1525 bfenumerator *bfe;
1526 bfc_term_base *told;
1527 int k;
1528 bf_reducer *bfr;
1529
1530 bfe_cum(evalue *factor, evalue *v, dpoly_r *r, bf_reducer *bfr,
1531 bfc_term_base *t, int k, bfenumerator *e) :
1532 cumulator(factor, v, r), told(t), k(k),
1533 bfr(bfr), bfe(e) {
1534 }
1535
1536 virtual void add_term(int *powers, int len, evalue *f2);
1537 };
1538
1539 void bfe_cum::add_term(int *powers, int len, evalue *f2)
1540 {
1541 bfr->update_powers(powers, len);
1542
1543 bfc_term_base * t = bfe->find_bfc_term(bfr->vn, bfr->npowers, bfr->nnf);
1544 bfe->set_factor(f2, bfr->l_changes % 2);
1545 bfe->add_term(t, told->terms[k], bfr->l_extra_num);
1546 }
1547
1548 void bfenumerator::cum(bf_reducer *bfr, bfc_term_base *t, int k,
1549 dpoly_r *r)
1550 {
1551 bfe_term* bfet = static_cast<bfe_term *>(t);
1552 bfe_cum cum(bfet->factors[k], E_num(0, bfr->d), r, bfr, t, k, this);
1553 cum.cumulate();
1554 }
1555
1556 void bfenumerator::base(mat_ZZ& factors, bfc_vec& v)
1557 {
1558 for (int i = 0; i < v.size(); ++i) {
1559 assert(v[i]->terms.NumRows() == 1);
1560 evalue *factor = static_cast<bfe_term *>(v[i])->factors[0];
1561 eadd(factor, vE[vert]);
1562 delete v[i];
1563 }
1564 }
1565
1566 void bfenumerator::handle_polar(Polyhedron *C, int s)
1567 {
1568 assert(C->NbRays-1 == enumerator_base::dim);
1569
1570 bfe_term* t = new bfe_term(enumerator_base::dim);
1571 vector< bfc_term_base * > v;
1572 v.push_back(t);
1573
1574 t->factors.resize(1);
1575
1576 t->terms.SetDims(1, enumerator_base::dim);
1577 lattice_point(V, C, t->terms[0], E_vertex);
1578
1579 // the elements of factors are always lexpositive
1580 mat_ZZ factors;
1581 s = setup_factors(C, factors, t, s);
1582
1583 t->factors[0] = new evalue;
1584 value_init(t->factors[0]->d);
1585 evalue_set_si(t->factors[0], s, 1);
1586 reduce(factors, v);
1587
1588 for (int i = 0; i < enumerator_base::dim; ++i)
1589 if (E_vertex[i]) {
1590 free_evalue_refs(E_vertex[i]);
1591 delete E_vertex[i];
1592 }
1593 }
1594
1595 #ifdef HAVE_CORRECT_VERTICES
1596 static inline Param_Polyhedron *Polyhedron2Param_SD(Polyhedron **Din,
1597 Polyhedron *Cin,int WS,Polyhedron **CEq,Matrix **CT)
1598 {
1599 if (WS & POL_NO_DUAL)
1600 WS = 0;
1601 return Polyhedron2Param_SimplifiedDomain(Din, Cin, WS, CEq, CT);
1602 }
1603 #else
1604 static Param_Polyhedron *Polyhedron2Param_SD(Polyhedron **Din,
1605 Polyhedron *Cin,int WS,Polyhedron **CEq,Matrix **CT)
1606 {
1607 static char data[] = " 1 0 0 0 0 1 -18 "
1608 " 1 0 0 -20 0 19 1 "
1609 " 1 0 1 20 0 -20 16 "
1610 " 1 0 0 0 0 -1 19 "
1611 " 1 0 -1 0 0 0 4 "
1612 " 1 4 -20 0 0 -1 23 "
1613 " 1 -4 20 0 0 1 -22 "
1614 " 1 0 1 0 20 -20 16 "
1615 " 1 0 0 0 -20 19 1 ";
1616 static int checked = 0;
1617 if (!checked) {
1618 checked = 1;
1619 char *p = data;
1620 int n, v, i;
1621 Matrix *M = Matrix_Alloc(9, 7);
1622 for (i = 0; i < 9; ++i)
1623 for (int j = 0; j < 7; ++j) {
1624 sscanf(p, "%d%n", &v, &n);
1625 p += n;
1626 value_set_si(M->p[i][j], v);
1627 }
1628 Polyhedron *P = Constraints2Polyhedron(M, 1024);
1629 Matrix_Free(M);
1630 Polyhedron *U = Universe_Polyhedron(1);
1631 Param_Polyhedron *PP = Polyhedron2Param_Domain(P, U, 1024);
1632 Polyhedron_Free(P);
1633 Polyhedron_Free(U);
1634 Param_Vertices *V;
1635 for (i = 0, V = PP->V; V; ++i, V = V->next)
1636 ;
1637 if (PP)
1638 Param_Polyhedron_Free(PP);
1639 if (i != 10) {
1640 fprintf(stderr, "WARNING: results may be incorrect\n");
1641 fprintf(stderr,
1642 "WARNING: use latest version of PolyLib to remove this warning\n");
1643 }
1644 }
1645
1646 return Polyhedron2Param_SimplifiedDomain(Din, Cin, WS, CEq, CT);
1647 }
1648 #endif
1649
1650 static evalue* barvinok_enumerate_ev_f(Polyhedron *P, Polyhedron* C,
1651 unsigned MaxRays);
1652
1653 /* Destroys C */
1654 static evalue* barvinok_enumerate_cst(Polyhedron *P, Polyhedron* C,
1655 unsigned MaxRays)
1656 {
1657 evalue *eres;
1658
1659 ALLOC(evalue, eres);
1660 value_init(eres->d);
1661 value_set_si(eres->d, 0);
1662 eres->x.p = new_enode(partition, 2, C->Dimension);
1663 EVALUE_SET_DOMAIN(eres->x.p->arr[0], DomainConstraintSimplify(C, MaxRays));
1664 value_set_si(eres->x.p->arr[1].d, 1);
1665 value_init(eres->x.p->arr[1].x.n);
1666 if (emptyQ(P))
1667 value_set_si(eres->x.p->arr[1].x.n, 0);
1668 else
1669 barvinok_count(P, &eres->x.p->arr[1].x.n, MaxRays);
1670
1671 return eres;
1672 }
1673
1674 evalue* barvinok_enumerate_ev(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
1675 {
1676 //P = unfringe(P, MaxRays);
1677 Polyhedron *Corig = C;
1678 Polyhedron *CEq = NULL, *rVD, *CA;
1679 int r = 0;
1680 unsigned nparam = C->Dimension;
1681 evalue *eres;
1682
1683 evalue factor;
1684 value_init(factor.d);
1685 evalue_set_si(&factor, 1, 1);
1686
1687 CA = align_context(C, P->Dimension, MaxRays);
1688 P = DomainIntersection(P, CA, MaxRays);
1689 Polyhedron_Free(CA);
1690
1691 /* for now */
1692 POL_ENSURE_FACETS(P);
1693 POL_ENSURE_VERTICES(P);
1694 POL_ENSURE_FACETS(C);
1695 POL_ENSURE_VERTICES(C);
1696
1697 if (C->Dimension == 0 || emptyQ(P)) {
1698 constant:
1699 eres = barvinok_enumerate_cst(P, CEq ? CEq : Polyhedron_Copy(C),
1700 MaxRays);
1701 out:
1702 emul(&factor, eres);
1703 reduce_evalue(eres);
1704 free_evalue_refs(&factor);
1705 Domain_Free(P);
1706 if (C != Corig)
1707 Polyhedron_Free(C);
1708
1709 return eres;
1710 }
1711 if (Polyhedron_is_infinite_param(P, nparam))
1712 goto constant;
1713
1714 if (P->NbEq != 0) {
1715 Matrix *f;
1716 P = remove_equalities_p(P, P->Dimension-nparam, &f);
1717 mask(f, &factor);
1718 Matrix_Free(f);
1719 }
1720 if (P->Dimension == nparam) {
1721 CEq = P;
1722 P = Universe_Polyhedron(0);
1723 goto constant;
1724 }
1725
1726 Polyhedron *T = Polyhedron_Factor(P, nparam, MaxRays);
1727 if (T || (P->Dimension == nparam+1)) {
1728 Polyhedron *Q;
1729 Polyhedron *C2;
1730 for (Q = T ? T : P; Q; Q = Q->next) {
1731 Polyhedron *next = Q->next;
1732 Q->next = NULL;
1733
1734 Polyhedron *QC = Q;
1735 if (Q->Dimension != C->Dimension)
1736 QC = Polyhedron_Project(Q, nparam);
1737
1738 C2 = C;
1739 C = DomainIntersection(C, QC, MaxRays);
1740 if (C2 != Corig)
1741 Polyhedron_Free(C2);
1742 if (QC != Q)
1743 Polyhedron_Free(QC);
1744
1745 Q->next = next;
1746 }
1747 }
1748 if (T) {
1749 Polyhedron_Free(P);
1750 P = T;
1751 if (T->Dimension == C->Dimension) {
1752 P = T->next;
1753 T->next = NULL;
1754 Polyhedron_Free(T);
1755 }
1756 }
1757
1758 Polyhedron *next = P->next;
1759 P->next = NULL;
1760 eres = barvinok_enumerate_ev_f(P, C, MaxRays);
1761 P->next = next;
1762
1763 if (P->next) {
1764 Polyhedron *Q;
1765 evalue *f;
1766
1767 for (Q = P->next; Q; Q = Q->next) {
1768 Polyhedron *next = Q->next;
1769 Q->next = NULL;
1770
1771 f = barvinok_enumerate_ev_f(Q, C, MaxRays);
1772 emul(f, eres);
1773 free_evalue_refs(f);
1774 free(f);
1775
1776 Q->next = next;
1777 }
1778 }
1779
1780 goto out;
1781 }
1782
1783 static evalue* barvinok_enumerate_ev_f(Polyhedron *P, Polyhedron* C,
1784 unsigned MaxRays)
1785 {
1786 unsigned nparam = C->Dimension;
1787
1788 if (P->Dimension - nparam == 1)
1789 return ParamLine_Length(P, C, MaxRays);
1790
1791 Param_Polyhedron *PP = NULL;
1792 Polyhedron *CEq = NULL, *pVD;
1793 Matrix *CT = NULL;
1794 Param_Domain *D, *next;
1795 Param_Vertices *V;
1796 evalue *eres;
1797 Polyhedron *Porig = P;
1798
1799 PP = Polyhedron2Param_SD(&P,C,MaxRays,&CEq,&CT);
1800
1801 if (isIdentity(CT)) {
1802 Matrix_Free(CT);
1803 CT = NULL;
1804 } else {
1805 assert(CT->NbRows != CT->NbColumns);
1806 if (CT->NbRows == 1) { // no more parameters
1807 eres = barvinok_enumerate_cst(P, CEq, MaxRays);
1808 out:
1809 if (CT)
1810 Matrix_Free(CT);
1811 if (PP)
1812 Param_Polyhedron_Free(PP);
1813 if (P != Porig)
1814 Polyhedron_Free(P);
1815
1816 return eres;
1817 }
1818 nparam = CT->NbRows - 1;
1819 }
1820
1821 unsigned dim = P->Dimension - nparam;
1822
1823 ALLOC(evalue, eres);
1824 value_init(eres->d);
1825 value_set_si(eres->d, 0);
1826
1827 int nd;
1828 for (nd = 0, D=PP->D; D; ++nd, D=D->next);
1829 struct section { Polyhedron *D; evalue E; };
1830 section *s = new section[nd];
1831 Polyhedron **fVD = new Polyhedron_p[nd];
1832
1833 try_again:
1834 #ifdef USE_INCREMENTAL_BF
1835 bfenumerator et(P, dim, PP->nbV);
1836 #elif defined USE_INCREMENTAL_DF
1837 ienumerator et(P, dim, PP->nbV);
1838 #else
1839 enumerator et(P, dim, PP->nbV);
1840 #endif
1841
1842 for(nd = 0, D=PP->D; D; D=next) {
1843 next = D->next;
1844
1845 Polyhedron *rVD = reduce_domain(D->Domain, CT, CEq,
1846 fVD, nd, MaxRays);
1847 if (!rVD)
1848 continue;
1849
1850 pVD = CT ? DomainImage(rVD,CT,MaxRays) : rVD;
1851
1852 value_init(s[nd].E.d);
1853 evalue_set_si(&s[nd].E, 0, 1);
1854 s[nd].D = rVD;
1855
1856 FORALL_PVertex_in_ParamPolyhedron(V,D,PP) // _i is internal counter
1857 if (!et.vE[_i])
1858 try {
1859 et.decompose_at(V, _i, MaxRays);
1860 } catch (OrthogonalException &e) {
1861 if (rVD != pVD)
1862 Domain_Free(pVD);
1863 for (; nd >= 0; --nd) {
1864 free_evalue_refs(&s[nd].E);
1865 Domain_Free(s[nd].D);
1866 Domain_Free(fVD[nd]);
1867 }
1868 goto try_again;
1869 }
1870 eadd(et.vE[_i] , &s[nd].E);
1871 END_FORALL_PVertex_in_ParamPolyhedron;
1872 evalue_range_reduction_in_domain(&s[nd].E, pVD);
1873
1874 if (CT)
1875 addeliminatedparams_evalue(&s[nd].E, CT);
1876 ++nd;
1877 if (rVD != pVD)
1878 Domain_Free(pVD);
1879 }
1880
1881 if (nd == 0)
1882 evalue_set_si(eres, 0, 1);
1883 else {
1884 eres->x.p = new_enode(partition, 2*nd, C->Dimension);
1885 for (int j = 0; j < nd; ++j) {
1886 EVALUE_SET_DOMAIN(eres->x.p->arr[2*j], s[j].D);
1887 value_clear(eres->x.p->arr[2*j+1].d);
1888 eres->x.p->arr[2*j+1] = s[j].E;
1889 Domain_Free(fVD[j]);
1890 }
1891 }
1892 delete [] s;
1893 delete [] fVD;
1894
1895 if (CEq)
1896 Polyhedron_Free(CEq);
1897 goto out;
1898 }
1899
1900 Enumeration* barvinok_enumerate(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
1901 {
1902 evalue *EP = barvinok_enumerate_ev(P, C, MaxRays);
1903
1904 return partition2enumeration(EP);
1905 }
1906
1907 static void SwapColumns(Value **V, int n, int i, int j)
1908 {
1909 for (int r = 0; r < n; ++r)
1910 value_swap(V[r][i], V[r][j]);
1911 }
1912
1913 static void SwapColumns(Polyhedron *P, int i, int j)
1914 {
1915 SwapColumns(P->Constraint, P->NbConstraints, i, j);
1916 SwapColumns(P->Ray, P->NbRays, i, j);
1917 }
1918
1919 /* Construct a constraint c from constraints l and u such that if
1920 * if constraint c holds then for each value of the other variables
1921 * there is at most one value of variable pos (position pos+1 in the constraints).
1922 *
1923 * Given a lower and an upper bound
1924 * n_l v_i + <c_l,x> + c_l >= 0
1925 * -n_u v_i + <c_u,x> + c_u >= 0
1926 * the constructed constraint is
1927 *
1928 * -(n_l<c_u,x> + n_u<c_l,x>) + (-n_l c_u - n_u c_l + n_l n_u - 1)
1929 *
1930 * which is then simplified to remove the content of the non-constant coefficients
1931 *
1932 * len is the total length of the constraints.
1933 * v is a temporary variable that can be used by this procedure
1934 */
1935 static void negative_test_constraint(Value *l, Value *u, Value *c, int pos,
1936 int len, Value *v)
1937 {
1938 value_oppose(*v, u[pos+1]);
1939 Vector_Combine(l+1, u+1, c+1, *v, l[pos+1], len-1);
1940 value_multiply(*v, *v, l[pos+1]);
1941 value_subtract(c[len-1], c[len-1], *v);
1942 value_set_si(*v, -1);
1943 Vector_Scale(c+1, c+1, *v, len-1);
1944 value_decrement(c[len-1], c[len-1]);
1945 ConstraintSimplify(c, c, len, v);
1946 }
1947
1948 static bool parallel_constraints(Value *l, Value *u, Value *c, int pos,
1949 int len)
1950 {
1951 bool parallel;
1952 Value g1;
1953 Value g2;
1954 value_init(g1);
1955 value_init(g2);
1956
1957 Vector_Gcd(&l[1+pos], len, &g1);
1958 Vector_Gcd(&u[1+pos], len, &g2);
1959 Vector_Combine(l+1+pos, u+1+pos, c+1, g2, g1, len);
1960 parallel = First_Non_Zero(c+1, len) == -1;
1961
1962 value_clear(g1);
1963 value_clear(g2);
1964
1965 return parallel;
1966 }
1967
1968 static void negative_test_constraint7(Value *l, Value *u, Value *c, int pos,
1969 int exist, int len, Value *v)
1970 {
1971 Value g;
1972 value_init(g);
1973
1974 Vector_Gcd(&u[1+pos], exist, v);
1975 Vector_Gcd(&l[1+pos], exist, &g);
1976 Vector_Combine(l+1, u+1, c+1, *v, g, len-1);
1977 value_multiply(*v, *v, g);
1978 value_subtract(c[len-1], c[len-1], *v);
1979 value_set_si(*v, -1);
1980 Vector_Scale(c+1, c+1, *v, len-1);
1981 value_decrement(c[len-1], c[len-1]);
1982 ConstraintSimplify(c, c, len, v);
1983
1984 value_clear(g);
1985 }
1986
1987 /* Turns a x + b >= 0 into a x + b <= -1
1988 *
1989 * len is the total length of the constraint.
1990 * v is a temporary variable that can be used by this procedure
1991 */
1992 static void oppose_constraint(Value *c, int len, Value *v)
1993 {
1994 value_set_si(*v, -1);
1995 Vector_Scale(c+1, c+1, *v, len-1);
1996 value_decrement(c[len-1], c[len-1]);
1997 }
1998
1999 /* Split polyhedron P into two polyhedra *pos and *neg, where
2000 * existential variable i has at most one solution for each
2001 * value of the other variables in *neg.
2002 *
2003 * The splitting is performed using constraints l and u.
2004 *
2005 * nvar: number of set variables
2006 * row: temporary vector that can be used by this procedure
2007 * f: temporary value that can be used by this procedure
2008 */
2009 static bool SplitOnConstraint(Polyhedron *P, int i, int l, int u,
2010 int nvar, int MaxRays, Vector *row, Value& f,
2011 Polyhedron **pos, Polyhedron **neg)
2012 {
2013 negative_test_constraint(P->Constraint[l], P->Constraint[u],
2014 row->p, nvar+i, P->Dimension+2, &f);
2015 *neg = AddConstraints(row->p, 1, P, MaxRays);
2016
2017 /* We found an independent, but useless constraint
2018 * Maybe we should detect this earlier and not
2019 * mark the variable as INDEPENDENT
2020 */
2021 if (emptyQ((*neg))) {
2022 Polyhedron_Free(*neg);
2023 return false;
2024 }
2025
2026 oppose_constraint(row->p, P->Dimension+2, &f);
2027 *pos = AddConstraints(row->p, 1, P, MaxRays);
2028
2029 if (emptyQ((*pos))) {
2030 Polyhedron_Free(*neg);
2031 Polyhedron_Free(*pos);
2032 return false;
2033 }
2034
2035 return true;
2036 }
2037
2038 /*
2039 * unimodularly transform P such that constraint r is transformed
2040 * into a constraint that involves only a single (the first)
2041 * existential variable
2042 *
2043 */
2044 static Polyhedron *rotate_along(Polyhedron *P, int r, int nvar, int exist,
2045 unsigned MaxRays)
2046 {
2047 Value g;
2048 value_init(g);
2049
2050 Vector *row = Vector_Alloc(exist);
2051 Vector_Copy(P->Constraint[r]+1+nvar, row->p, exist);
2052 Vector_Gcd(row->p, exist, &g);
2053 if (value_notone_p(g))
2054 Vector_AntiScale(row->p, row->p, g, exist);
2055 value_clear(g);
2056
2057 Matrix *M = unimodular_complete(row);
2058 Matrix *M2 = Matrix_Alloc(P->Dimension+1, P->Dimension+1);
2059 for (r = 0; r < nvar; ++r)
2060 value_set_si(M2->p[r][r], 1);
2061 for ( ; r < nvar+exist; ++r)
2062 Vector_Copy(M->p[r-nvar], M2->p[r]+nvar, exist);
2063 for ( ; r < P->Dimension+1; ++r)
2064 value_set_si(M2->p[r][r], 1);
2065 Polyhedron *T = Polyhedron_Image(P, M2, MaxRays);
2066
2067 Matrix_Free(M2);
2068 Matrix_Free(M);
2069 Vector_Free(row);
2070
2071 return T;
2072 }
2073
2074 /* Split polyhedron P into two polyhedra *pos and *neg, where
2075 * existential variable i has at most one solution for each
2076 * value of the other variables in *neg.
2077 *
2078 * If independent is set, then the two constraints on which the
2079 * split will be performed need to be independent of the other
2080 * existential variables.
2081 *
2082 * Return true if an appropriate split could be performed.
2083 *
2084 * nvar: number of set variables
2085 * exist: number of existential variables
2086 * row: temporary vector that can be used by this procedure
2087 * f: temporary value that can be used by this procedure
2088 */
2089 static bool SplitOnVar(Polyhedron *P, int i,
2090 int nvar, int exist, int MaxRays,
2091 Vector *row, Value& f, bool independent,
2092 Polyhedron **pos, Polyhedron **neg)
2093 {
2094 int j;
2095
2096 for (int l = P->NbEq; l < P->NbConstraints; ++l) {
2097 if (value_negz_p(P->Constraint[l][nvar+i+1]))
2098 continue;
2099
2100 if (independent) {
2101 for (j = 0; j < exist; ++j)
2102 if (j != i && value_notzero_p(P->Constraint[l][nvar+j+1]))
2103 break;
2104 if (j < exist)
2105 continue;
2106 }
2107
2108 for (int u = P->NbEq; u < P->NbConstraints; ++u) {
2109 if (value_posz_p(P->Constraint[u][nvar+i+1]))
2110 continue;
2111
2112 if (independent) {
2113 for (j = 0; j < exist; ++j)
2114 if (j != i && value_notzero_p(P->Constraint[u][nvar+j+1]))
2115 break;
2116 if (j < exist)
2117 continue;
2118 }
2119
2120 if (SplitOnConstraint(P, i, l, u, nvar, MaxRays, row, f, pos, neg)) {
2121 if (independent) {
2122 if (i != 0)
2123 SwapColumns(*neg, nvar+1, nvar+1+i);
2124 }
2125 return true;
2126 }
2127 }
2128 }
2129
2130 return false;
2131 }
2132
2133 static bool double_bound_pair(Polyhedron *P, int nvar, int exist,
2134 int i, int l1, int l2,
2135 Polyhedron **pos, Polyhedron **neg)
2136 {
2137 Value f;
2138 value_init(f);
2139 Vector *row = Vector_Alloc(P->Dimension+2);
2140 value_set_si(row->p[0], 1);
2141 value_oppose(f, P->Constraint[l1][nvar+i+1]);
2142 Vector_Combine(P->Constraint[l1]+1, P->Constraint[l2]+1,
2143 row->p+1,
2144 P->Constraint[l2][nvar+i+1], f,
2145 P->Dimension+1);
2146 ConstraintSimplify(row->p, row->p, P->Dimension+2, &f);
2147 *pos = AddConstraints(row->p, 1, P, 0);
2148 value_set_si(f, -1);
2149 Vector_Scale(row->p+1, row->p+1, f, P->Dimension+1);
2150 value_decrement(row->p[P->Dimension+1], row->p[P->Dimension+1]);
2151 *neg = AddConstraints(row->p, 1, P, 0);
2152 Vector_Free(row);
2153 value_clear(f);
2154
2155 return !emptyQ((*pos)) && !emptyQ((*neg));
2156 }
2157
2158 static bool double_bound(Polyhedron *P, int nvar, int exist,
2159 Polyhedron **pos, Polyhedron **neg)
2160 {
2161 for (int i = 0; i < exist; ++i) {
2162 int l1, l2;
2163 for (l1 = P->NbEq; l1 < P->NbConstraints; ++l1) {
2164 if (value_negz_p(P->Constraint[l1][nvar+i+1]))
2165 continue;
2166 for (l2 = l1 + 1; l2 < P->NbConstraints; ++l2) {
2167 if (value_negz_p(P->Constraint[l2][nvar+i+1]))
2168 continue;
2169 if (double_bound_pair(P, nvar, exist, i, l1, l2, pos, neg))
2170 return true;
2171 }
2172 }
2173 for (l1 = P->NbEq; l1 < P->NbConstraints; ++l1) {
2174 if (value_posz_p(P->Constraint[l1][nvar+i+1]))
2175 continue;
2176 if (l1 < P->NbConstraints)
2177 for (l2 = l1 + 1; l2 < P->NbConstraints; ++l2) {
2178 if (value_posz_p(P->Constraint[l2][nvar+i+1]))
2179 continue;
2180 if (double_bound_pair(P, nvar, exist, i, l1, l2, pos, neg))
2181 return true;
2182 }
2183 }
2184 return false;
2185 }
2186 return false;
2187 }
2188
2189 enum constraint {
2190 ALL_POS = 1 << 0,
2191 ONE_NEG = 1 << 1,
2192 INDEPENDENT = 1 << 2,
2193 ROT_NEG = 1 << 3
2194 };
2195
2196 static evalue* enumerate_or(Polyhedron *D,
2197 unsigned exist, unsigned nparam, unsigned MaxRays)
2198 {
2199 #ifdef DEBUG_ER
2200 fprintf(stderr, "\nER: Or\n");
2201 #endif /* DEBUG_ER */
2202
2203 Polyhedron *N = D->next;
2204 D->next = 0;
2205 evalue *EP =
2206 barvinok_enumerate_e(D, exist, nparam, MaxRays);
2207 Polyhedron_Free(D);
2208
2209 for (D = N; D; D = N) {
2210 N = D->next;
2211 D->next = 0;
2212
2213 evalue *EN =
2214 barvinok_enumerate_e(D, exist, nparam, MaxRays);
2215
2216 eor(EN, EP);
2217 free_evalue_refs(EN);
2218 free(EN);
2219 Polyhedron_Free(D);
2220 }
2221
2222 reduce_evalue(EP);
2223
2224 return EP;
2225 }
2226
2227 static evalue* enumerate_sum(Polyhedron *P,
2228 unsigned exist, unsigned nparam, unsigned MaxRays)
2229 {
2230 int nvar = P->Dimension - exist - nparam;
2231 int toswap = nvar < exist ? nvar : exist;
2232 for (int i = 0; i < toswap; ++i)
2233 SwapColumns(P, 1 + i, nvar+exist - i);
2234 nparam += nvar;
2235
2236 #ifdef DEBUG_ER
2237 fprintf(stderr, "\nER: Sum\n");
2238 #endif /* DEBUG_ER */
2239
2240 evalue *EP = barvinok_enumerate_e(P, exist, nparam, MaxRays);
2241
2242 for (int i = 0; i < /* nvar */ nparam; ++i) {
2243 Matrix *C = Matrix_Alloc(1, 1 + nparam + 1);
2244 value_set_si(C->p[0][0], 1);
2245 evalue split;
2246 value_init(split.d);
2247 value_set_si(split.d, 0);
2248 split.x.p = new_enode(partition, 4, nparam);
2249 value_set_si(C->p[0][1+i], 1);
2250 Matrix *C2 = Matrix_Copy(C);
2251 EVALUE_SET_DOMAIN(split.x.p->arr[0],
2252 Constraints2Polyhedron(C2, MaxRays));
2253 Matrix_Free(C2);
2254 evalue_set_si(&split.x.p->arr[1], 1, 1);
2255 value_set_si(C->p[0][1+i], -1);
2256 value_set_si(C->p[0][1+nparam], -1);
2257 EVALUE_SET_DOMAIN(split.x.p->arr[2],
2258 Constraints2Polyhedron(C, MaxRays));
2259 evalue_set_si(&split.x.p->arr[3], 1, 1);
2260 emul(&split, EP);
2261 free_evalue_refs(&split);
2262 Matrix_Free(C);
2263 }
2264 reduce_evalue(EP);
2265 evalue_range_reduction(EP);
2266
2267 evalue_frac2floor(EP);
2268
2269 evalue *sum = esum(EP, nvar);
2270
2271 free_evalue_refs(EP);
2272 free(EP);
2273 EP = sum;
2274
2275 evalue_range_reduction(EP);
2276
2277 return EP;
2278 }
2279
2280 static evalue* split_sure(Polyhedron *P, Polyhedron *S,
2281 unsigned exist, unsigned nparam, unsigned MaxRays)
2282 {
2283 int nvar = P->Dimension - exist - nparam;
2284
2285 Matrix *M = Matrix_Alloc(exist, S->Dimension+2);
2286 for (int i = 0; i < exist; ++i)
2287 value_set_si(M->p[i][nvar+i+1], 1);
2288 Polyhedron *O = S;
2289 S = DomainAddRays(S, M, MaxRays);
2290 Polyhedron_Free(O);
2291 Polyhedron *F = DomainAddRays(P, M, MaxRays);
2292 Polyhedron *D = DomainDifference(F, S, MaxRays);
2293 O = D;
2294 D = Disjoint_Domain(D, 0, MaxRays);
2295 Polyhedron_Free(F);
2296 Domain_Free(O);
2297 Matrix_Free(M);
2298
2299 M = Matrix_Alloc(P->Dimension+1-exist, P->Dimension+1);
2300 for (int j = 0; j < nvar; ++j)
2301 value_set_si(M->p[j][j], 1);
2302 for (int j = 0; j < nparam+1; ++j)
2303 value_set_si(M->p[nvar+j][nvar+exist+j], 1);
2304 Polyhedron *T = Polyhedron_Image(S, M, MaxRays);
2305 evalue *EP = barvinok_enumerate_e(T, 0, nparam, MaxRays);
2306 Polyhedron_Free(S);
2307 Polyhedron_Free(T);
2308 Matrix_Free(M);
2309
2310 for (Polyhedron *Q = D; Q; Q = Q->next) {
2311 Polyhedron *N = Q->next;
2312 Q->next = 0;
2313 T = DomainIntersection(P, Q, MaxRays);
2314 evalue *E = barvinok_enumerate_e(T, exist, nparam, MaxRays);
2315 eadd(E, EP);
2316 free_evalue_refs(E);
2317 free(E);
2318 Polyhedron_Free(T);
2319 Q->next = N;
2320 }
2321 Domain_Free(D);
2322 return EP;
2323 }
2324
2325 static evalue* enumerate_sure(Polyhedron *P,
2326 unsigned exist, unsigned nparam, unsigned MaxRays)
2327 {
2328 int i;
2329 Polyhedron *S = P;
2330 int nvar = P->Dimension - exist - nparam;
2331 Value lcm;
2332 Value f;
2333 value_init(lcm);
2334 value_init(f);
2335
2336 for (i = 0; i < exist; ++i) {
2337 Matrix *M = Matrix_Alloc(S->NbConstraints, S->Dimension+2);
2338 int c = 0;
2339 value_set_si(lcm, 1);
2340 for (int j = 0; j < S->NbConstraints; ++j) {
2341 if (value_negz_p(S->Constraint[j][1+nvar+i]))
2342 continue;
2343 if (value_one_p(S->Constraint[j][1+nvar+i]))
2344 continue;
2345 value_lcm(lcm, S->Constraint[j][1+nvar+i], &lcm);
2346 }
2347
2348 for (int j = 0; j < S->NbConstraints; ++j) {
2349 if (value_negz_p(S->Constraint[j][1+nvar+i]))
2350 continue;
2351 if (value_one_p(S->Constraint[j][1+nvar+i]))
2352 continue;
2353 value_division(f, lcm, S->Constraint[j][1+nvar+i]);
2354 Vector_Scale(S->Constraint[j], M->p[c], f, S->Dimension+2);
2355 value_subtract(M->p[c][S->Dimension+1],
2356 M->p[c][S->Dimension+1],
2357 lcm);
2358 value_increment(M->p[c][S->Dimension+1],
2359 M->p[c][S->Dimension+1]);
2360 ++c;
2361 }
2362 Polyhedron *O = S;
2363 S = AddConstraints(M->p[0], c, S, MaxRays);
2364 if (O != P)
2365 Polyhedron_Free(O);
2366 Matrix_Free(M);
2367 if (emptyQ(S)) {
2368 Polyhedron_Free(S);
2369 value_clear(lcm);
2370 value_clear(f);
2371 return 0;
2372 }
2373 }
2374 value_clear(lcm);
2375 value_clear(f);
2376
2377 #ifdef DEBUG_ER
2378 fprintf(stderr, "\nER: Sure\n");
2379 #endif /* DEBUG_ER */
2380
2381 return split_sure(P, S, exist, nparam, MaxRays);
2382 }
2383
2384 static evalue* enumerate_sure2(Polyhedron *P,
2385 unsigned exist, unsigned nparam, unsigned MaxRays)
2386 {
2387 int nvar = P->Dimension - exist - nparam;
2388 int r;
2389 for (r = 0; r < P->NbRays; ++r)
2390 if (value_one_p(P->Ray[r][0]) &&
2391 value_one_p(P->Ray[r][P->Dimension+1]))
2392 break;
2393
2394 if (r >= P->NbRays)
2395 return 0;
2396
2397 Matrix *M = Matrix_Alloc(nvar + 1 + nparam, P->Dimension+2);
2398 for (int i = 0; i < nvar; ++i)
2399 value_set_si(M->p[i][1+i], 1);
2400 for (int i = 0; i < nparam; ++i)
2401 value_set_si(M->p[i+nvar][1+nvar+exist+i], 1);
2402 Vector_Copy(P->Ray[r]+1+nvar, M->p[nvar+nparam]+1+nvar, exist);
2403 value_set_si(M->p[nvar+nparam][0], 1);
2404 value_set_si(M->p[nvar+nparam][P->Dimension+1], 1);
2405 Polyhedron * F = Rays2Polyhedron(M, MaxRays);
2406 Matrix_Free(M);
2407
2408 Polyhedron *I = DomainIntersection(F, P, MaxRays);
2409 Polyhedron_Free(F);
2410
2411 #ifdef DEBUG_ER
2412 fprintf(stderr, "\nER: Sure2\n");
2413 #endif /* DEBUG_ER */
2414
2415 return split_sure(P, I, exist, nparam, MaxRays);
2416 }
2417
2418 static evalue* enumerate_cyclic(Polyhedron *P,
2419 unsigned exist, unsigned nparam,
2420 evalue * EP, int r, int p, unsigned MaxRays)
2421 {
2422 int nvar = P->Dimension - exist - nparam;
2423
2424 /* If EP in its fractional maps only contains references
2425 * to the remainder parameter with appropriate coefficients
2426 * then we could in principle avoid adding existentially
2427 * quantified variables to the validity domains.
2428 * We'd have to replace the remainder by m { p/m }
2429 * and multiply with an appropriate factor that is one
2430 * only in the appropriate range.
2431 * This last multiplication can be avoided if EP
2432 * has a single validity domain with no (further)
2433 * constraints on the remainder parameter
2434 */
2435
2436 Matrix *CT = Matrix_Alloc(nparam+1, nparam+3);
2437 Matrix *M = Matrix_Alloc(1, 1+nparam+3);
2438 for (int j = 0; j < nparam; ++j)
2439 if (j != p)
2440 value_set_si(CT->p[j][j], 1);
2441 value_set_si(CT->p[p][nparam+1], 1);
2442 value_set_si(CT->p[nparam][nparam+2], 1);
2443 value_set_si(M->p[0][1+p], -1);
2444 value_absolute(M->p[0][1+nparam], P->Ray[0][1+nvar+exist+p]);
2445 value_set_si(M->p[0][1+nparam+1], 1);
2446 Polyhedron *CEq = Constraints2Polyhedron(M, 1);
2447 Matrix_Free(M);
2448 addeliminatedparams_enum(EP, CT, CEq, MaxRays, nparam);
2449 Polyhedron_Free(CEq);
2450 Matrix_Free(CT);
2451
2452 return EP;
2453 }
2454
2455 static void enumerate_vd_add_ray(evalue *EP, Matrix *Rays, unsigned MaxRays)
2456 {
2457 if (value_notzero_p(EP->d))
2458 return;
2459
2460 assert(EP->x.p->type == partition);
2461 assert(EP->x.p->pos == EVALUE_DOMAIN(EP->x.p->arr[0])->Dimension);
2462 for (int i = 0; i < EP->x.p->size/2; ++i) {
2463 Polyhedron *D = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
2464 Polyhedron *N = DomainAddRays(D, Rays, MaxRays);
2465 EVALUE_SET_DOMAIN(EP->x.p->arr[2*i], N);
2466 Domain_Free(D);
2467 }
2468 }
2469
2470 static evalue* enumerate_line(Polyhedron *P,
2471 unsigned exist, unsigned nparam, unsigned MaxRays)
2472 {
2473 if (P->NbBid == 0)
2474 return 0;
2475
2476 #ifdef DEBUG_ER
2477 fprintf(stderr, "\nER: Line\n");
2478 #endif /* DEBUG_ER */
2479
2480 int nvar = P->Dimension - exist - nparam;
2481 int i, j;
2482 for (i = 0; i < nparam; ++i)
2483 if (value_notzero_p(P->Ray[0][1+nvar+exist+i]))
2484 break;
2485 assert(i < nparam);
2486 for (j = i+1; j < nparam; ++j)
2487 if (value_notzero_p(P->Ray[0][1+nvar+exist+i]))
2488 break;
2489 assert(j >= nparam); // for now
2490
2491 Matrix *M = Matrix_Alloc(2, P->Dimension+2);
2492 value_set_si(M->p[0][0], 1);
2493 value_set_si(M->p[0][1+nvar+exist+i], 1);
2494 value_set_si(M->p[1][0], 1);
2495 value_set_si(M->p[1][1+nvar+exist+i], -1);
2496 value_absolute(M->p[1][1+P->Dimension], P->Ray[0][1+nvar+exist+i]);
2497 value_decrement(M->p[1][1+P->Dimension], M->p[1][1+P->Dimension]);
2498 Polyhedron *S = AddConstraints(M->p[0], 2, P, MaxRays);
2499 evalue *EP = barvinok_enumerate_e(S, exist, nparam, MaxRays);
2500 Polyhedron_Free(S);
2501 Matrix_Free(M);
2502
2503 return enumerate_cyclic(P, exist, nparam, EP, 0, i, MaxRays);
2504 }
2505
2506 static int single_param_pos(Polyhedron*P, unsigned exist, unsigned nparam,
2507 int r)
2508 {
2509 int nvar = P->Dimension - exist - nparam;
2510 if (First_Non_Zero(P->Ray[r]+1, nvar) != -1)
2511 return -1;
2512 int i = First_Non_Zero(P->Ray[r]+1+nvar+exist, nparam);
2513 if (i == -1)
2514 return -1;
2515 if (First_Non_Zero(P->Ray[r]+1+nvar+exist+1, nparam-i-1) != -1)
2516 return -1;
2517 return i;
2518 }
2519
2520 static evalue* enumerate_remove_ray(Polyhedron *P, int r,
2521 unsigned exist, unsigned nparam, unsigned MaxRays)
2522 {
2523 #ifdef DEBUG_ER
2524 fprintf(stderr, "\nER: RedundantRay\n");
2525 #endif /* DEBUG_ER */
2526
2527 Value one;
2528 value_init(one);
2529 value_set_si(one, 1);
2530 int len = P->NbRays-1;
2531 Matrix *M = Matrix_Alloc(2 * len, P->Dimension+2);
2532 Vector_Copy(P->Ray[0], M->p[0], r * (P->Dimension+2));
2533 Vector_Copy(P->Ray[r+1], M->p[r], (len-r) * (P->Dimension+2));
2534 for (int j = 0; j < P->NbRays; ++j) {
2535 if (j == r)
2536 continue;
2537 Vector_Combine(P->Ray[j], P->Ray[r], M->p[len+j-(j>r)],
2538 one, P->Ray[j][P->Dimension+1], P->Dimension+2);
2539 }
2540
2541 P = Rays2Polyhedron(M, MaxRays);
2542 Matrix_Free(M);
2543 evalue *EP = barvinok_enumerate_e(P, exist, nparam, MaxRays);
2544 Polyhedron_Free(P);
2545 value_clear(one);
2546
2547 return EP;
2548 }
2549
2550 static evalue* enumerate_redundant_ray(Polyhedron *P,
2551 unsigned exist, unsigned nparam, unsigned MaxRays)
2552 {
2553 assert(P->NbBid == 0);
2554 int nvar = P->Dimension - exist - nparam;
2555 Value m;
2556 value_init(m);
2557
2558 for (int r = 0; r < P->NbRays; ++r) {
2559 if (value_notzero_p(P->Ray[r][P->Dimension+1]))
2560 continue;
2561 int i1 = single_param_pos(P, exist, nparam, r);
2562 if (i1 == -1)
2563 continue;
2564 for (int r2 = r+1; r2 < P->NbRays; ++r2) {
2565 if (value_notzero_p(P->Ray[r2][P->Dimension+1]))
2566 continue;
2567 int i2 = single_param_pos(P, exist, nparam, r2);
2568 if (i2 == -1)
2569 continue;
2570 if (i1 != i2)
2571 continue;
2572
2573 value_division(m, P->Ray[r][1+nvar+exist+i1],
2574 P->Ray[r2][1+nvar+exist+i1]);
2575 value_multiply(m, m, P->Ray[r2][1+nvar+exist+i1]);
2576 /* r2 divides r => r redundant */
2577 if (value_eq(m, P->Ray[r][1+nvar+exist+i1])) {
2578 value_clear(m);
2579 return enumerate_remove_ray(P, r, exist, nparam, MaxRays);
2580 }
2581
2582 value_division(m, P->Ray[r2][1+nvar+exist+i1],
2583 P->Ray[r][1+nvar+exist+i1]);
2584 value_multiply(m, m, P->Ray[r][1+nvar+exist+i1]);
2585 /* r divides r2 => r2 redundant */
2586 if (value_eq(m, P->Ray[r2][1+nvar+exist+i1])) {
2587 value_clear(m);
2588 return enumerate_remove_ray(P, r2, exist, nparam, MaxRays);
2589 }
2590 }
2591 }
2592 value_clear(m);
2593 return 0;
2594 }
2595
2596 static Polyhedron *upper_bound(Polyhedron *P,
2597 int pos, Value *max, Polyhedron **R)
2598 {
2599 Value v;
2600 int r;
2601 value_init(v);
2602
2603 *R = 0;
2604 Polyhedron *N;
2605 Polyhedron *B = 0;
2606 for (Polyhedron *Q = P; Q; Q = N) {
2607 N = Q->next;
2608 for (r = 0; r < P->NbRays; ++r) {
2609 if (value_zero_p(P->Ray[r][P->Dimension+1]) &&
2610 value_pos_p(P->Ray[r][1+pos]))
2611 break;
2612 }
2613 if (r < P->NbRays) {
2614 Q->next = *R;
2615 *R = Q;
2616 continue;
2617 } else {
2618 Q->next = B;
2619 B = Q;
2620 }
2621 for (r = 0; r < P->NbRays; ++r) {
2622 if (value_zero_p(P->Ray[r][P->Dimension+1]))
2623 continue;
2624 mpz_fdiv_q(v, P->Ray[r][1+pos], P->Ray[r][1+P->Dimension]);
2625 if ((!Q->next && r == 0) || value_gt(v, *max))
2626 value_assign(*max, v);
2627 }
2628 }
2629 value_clear(v);
2630 return B;
2631 }
2632
2633 static evalue* enumerate_ray(Polyhedron *P,
2634 unsigned exist, unsigned nparam, unsigned MaxRays)
2635 {
2636 assert(P->NbBid == 0);
2637 int nvar = P->Dimension - exist - nparam;
2638
2639 int r;
2640 for (r = 0; r < P->NbRays; ++r)
2641 if (value_zero_p(P->Ray[r][P->Dimension+1]))
2642 break;
2643 if (r >= P->NbRays)
2644 return 0;
2645
2646 int r2;
2647 for (r2 = r+1; r2 < P->NbRays; ++r2)
2648 if (value_zero_p(P->Ray[r2][P->Dimension+1]))
2649 break;
2650 if (r2 < P->NbRays) {
2651 if (nvar > 0)
2652 return enumerate_sum(P, exist, nparam, MaxRays);
2653 }
2654
2655 #ifdef DEBUG_ER
2656 fprintf(stderr, "\nER: Ray\n");
2657 #endif /* DEBUG_ER */
2658
2659 Value m;
2660 Value one;
2661 value_init(m);
2662 value_init(one);
2663 value_set_si(one, 1);
2664 int i = single_param_pos(P, exist, nparam, r);
2665 assert(i != -1); // for now;
2666
2667 Matrix *M = Matrix_Alloc(P->NbRays, P->Dimension+2);
2668 for (int j = 0; j < P->NbRays; ++j) {
2669 Vector_Combine(P->Ray[j], P->Ray[r], M->p[j],
2670 one, P->Ray[j][P->Dimension+1], P->Dimension+2);
2671 }
2672 Polyhedron *S = Rays2Polyhedron(M, MaxRays);
2673 Matrix_Free(M);
2674 Polyhedron *D = DomainDifference(P, S, MaxRays);
2675 Polyhedron_Free(S);
2676 // Polyhedron_Print(stderr, P_VALUE_FMT, D);
2677 assert(value_pos_p(P->Ray[r][1+nvar+exist+i])); // for now
2678 Polyhedron *R;
2679 D = upper_bound(D, nvar+exist+i, &m, &R);
2680 assert(D);
2681 Domain_Free(D);
2682
2683 M = Matrix_Alloc(2, P->Dimension+2);
2684 value_set_si(M->p[0][0], 1);
2685 value_set_si(M->p[1][0], 1);
2686 value_set_si(M->p[0][1+nvar+exist+i], -1);
2687 value_set_si(M->p[1][1+nvar+exist+i], 1);
2688 value_assign(M->p[0][1+P->Dimension], m);
2689 value_oppose(M->p[1][1+P->Dimension], m);
2690 value_addto(M->p[1][1+P->Dimension], M->p[1][1+P->Dimension],
2691 P->Ray[r][1+nvar+exist+i]);
2692 value_decrement(M->p[1][1+P->Dimension], M->p[1][1+P->Dimension]);
2693 // Matrix_Print(stderr, P_VALUE_FMT, M);
2694 D = AddConstraints(M->p[0], 2, P, MaxRays);
2695 // Polyhedron_Print(stderr, P_VALUE_FMT, D);
2696 value_subtract(M->p[0][1+P->Dimension], M->p[0][1+P->Dimension],
2697 P->Ray[r][1+nvar+exist+i]);
2698 // Matrix_Print(stderr, P_VALUE_FMT, M);
2699 S = AddConstraints(M->p[0], 1, P, MaxRays);
2700 // Polyhedron_Print(stderr, P_VALUE_FMT, S);
2701 Matrix_Free(M);
2702
2703 evalue *EP = barvinok_enumerate_e(D, exist, nparam, MaxRays);
2704 Polyhedron_Free(D);
2705 value_clear(one);
2706 value_clear(m);
2707
2708 if (value_notone_p(P->Ray[r][1+nvar+exist+i]))
2709 EP = enumerate_cyclic(P, exist, nparam, EP, r, i, MaxRays);
2710 else {
2711 M = Matrix_Alloc(1, nparam+2);
2712 value_set_si(M->p[0][0], 1);
2713 value_set_si(M->p[0][1+i], 1);
2714 enumerate_vd_add_ray(EP, M, MaxRays);
2715 Matrix_Free(M);
2716 }
2717
2718 if (!emptyQ(S)) {
2719 evalue *E = barvinok_enumerate_e(S, exist, nparam, MaxRays);
2720 eadd(E, EP);
2721 free_evalue_refs(E);
2722 free(E);
2723 }
2724 Polyhedron_Free(S);
2725
2726 if (R) {
2727 assert(nvar == 0);
2728 evalue *ER = enumerate_or(R, exist, nparam, MaxRays);
2729 eor(ER, EP);
2730 free_evalue_refs(ER);
2731 free(ER);
2732 }
2733
2734 return EP;
2735 }
2736
2737 static evalue* enumerate_vd(Polyhedron **PA,
2738 unsigned exist, unsigned nparam, unsigned MaxRays)
2739 {
2740 Polyhedron *P = *PA;
2741 int nvar = P->Dimension - exist - nparam;
2742 Param_Polyhedron *PP = NULL;
2743 Polyhedron *C = Universe_Polyhedron(nparam);
2744 Polyhedron *CEq;
2745 Matrix *CT;
2746 Polyhedron *PR = P;
2747 PP = Polyhedron2Param_SimplifiedDomain(&PR,C,MaxRays,&CEq,&CT);
2748 Polyhedron_Free(C);
2749
2750 int nd;
2751 Param_Domain *D, *last;
2752 Value c;
2753 value_init(c);
2754 for (nd = 0, D=PP->D; D; D=D->next, ++nd)
2755 ;
2756
2757 Polyhedron **VD = new Polyhedron_p[nd];
2758 Polyhedron **fVD = new Polyhedron_p[nd];
2759 for(nd = 0, D=PP->D; D; D=D->next) {
2760 Polyhedron *rVD = reduce_domain(D->Domain, CT, CEq,
2761 fVD, nd, MaxRays);
2762 if (!rVD)
2763 continue;
2764
2765 VD[nd++] = rVD;
2766 last = D;
2767 }
2768
2769 evalue *EP = 0;
2770
2771 if (nd == 0)
2772 EP = evalue_zero();
2773
2774 /* This doesn't seem to have any effect */
2775 if (nd == 1) {
2776 Polyhedron *CA = align_context(VD[0], P->Dimension, MaxRays);
2777 Polyhedron *O = P;
2778 P = DomainIntersection(P, CA, MaxRays);
2779 if (O != *PA)
2780 Polyhedron_Free(O);
2781 Polyhedron_Free(CA);
2782 if (emptyQ(P))
2783 EP = evalue_zero();
2784 }
2785
2786 if (!EP && CT->NbColumns != CT->NbRows) {
2787 Polyhedron *CEqr = DomainImage(CEq, CT, MaxRays);
2788 Polyhedron *CA = align_context(CEqr, PR->Dimension, MaxRays);
2789 Polyhedron *I = DomainIntersection(PR, CA, MaxRays);
2790 Polyhedron_Free(CEqr);
2791 Polyhedron_Free(CA);
2792 #ifdef DEBUG_ER
2793 fprintf(stderr, "\nER: Eliminate\n");
2794 #endif /* DEBUG_ER */
2795 nparam -= CT->NbColumns - CT->NbRows;
2796 EP = barvinok_enumerate_e(I, exist, nparam, MaxRays);
2797 nparam += CT->NbColumns - CT->NbRows;
2798 addeliminatedparams_enum(EP, CT, CEq, MaxRays, nparam);
2799 Polyhedron_Free(I);
2800 }
2801 if (PR != *PA)
2802 Polyhedron_Free(PR);
2803 PR = 0;
2804
2805 if (!EP && nd > 1) {
2806 #ifdef DEBUG_ER
2807 fprintf(stderr, "\nER: VD\n");
2808 #endif /* DEBUG_ER */
2809 for (int i = 0; i < nd; ++i) {
2810 Polyhedron *CA = align_context(VD[i], P->Dimension, MaxRays);
2811 Polyhedron *I = DomainIntersection(P, CA, MaxRays);
2812
2813 if (i == 0)
2814 EP = barvinok_enumerate_e(I, exist, nparam, MaxRays);
2815 else {
2816 evalue *E = barvinok_enumerate_e(I, exist, nparam, MaxRays);
2817 eadd(E, EP);
2818 free_evalue_refs(E);
2819 free(E);
2820 }
2821 Polyhedron_Free(I);
2822 Polyhedron_Free(CA);
2823 }
2824 }
2825
2826 for (int i = 0; i < nd; ++i) {
2827 Polyhedron_Free(VD[i]);
2828 Polyhedron_Free(fVD[i]);
2829 }
2830 delete [] VD;
2831 delete [] fVD;
2832 value_clear(c);
2833
2834 if (!EP && nvar == 0) {
2835 Value f;
2836 value_init(f);
2837 Param_Vertices *V, *V2;
2838 Matrix* M = Matrix_Alloc(1, P->Dimension+2);
2839
2840 FORALL_PVertex_in_ParamPolyhedron(V, last, PP) {
2841 bool found = false;
2842 FORALL_PVertex_in_ParamPolyhedron(V2, last, PP) {
2843 if (V == V2) {
2844 found = true;
2845 continue;
2846 }
2847 if (!found)
2848 continue;
2849 for (int i = 0; i < exist; ++i) {
2850 value_oppose(f, V->Vertex->p[i][nparam+1]);
2851 Vector_Combine(V->Vertex->p[i],
2852 V2->Vertex->p[i],
2853 M->p[0] + 1 + nvar + exist,
2854 V2->Vertex->p[i][nparam+1],
2855 f,
2856 nparam+1);
2857 int j;
2858 for (j = 0; j < nparam; ++j)
2859 if (value_notzero_p(M->p[0][1+nvar+exist+j]))
2860 break;
2861 if (j >= nparam)
2862 continue;
2863 ConstraintSimplify(M->p[0], M->p[0],
2864 P->Dimension+2, &f);
2865 value_set_si(M->p[0][0], 0);
2866 Polyhedron *para = AddConstraints(M->p[0], 1, P,
2867 MaxRays);
2868 if (emptyQ(para)) {
2869 Polyhedron_Free(para);
2870 continue;
2871 }
2872 Polyhedron *pos, *neg;
2873 value_set_si(M->p[0][0], 1);
2874 value_decrement(M->p[0][P->Dimension+1],
2875 M->p[0][P->Dimension+1]);
2876 neg = AddConstraints(M->p[0], 1, P, MaxRays);
2877 value_set_si(f, -1);
2878 Vector_Scale(M->p[0]+1, M->p[0]+1, f,
2879 P->Dimension+1);
2880 value_decrement(M->p[0][P->Dimension+1],
2881 M->p[0][P->Dimension+1]);
2882 value_decrement(M->p[0][P->Dimension+1],
2883 M->p[0][P->Dimension+1]);
2884 pos = AddConstraints(M->p[0], 1, P, MaxRays);
2885 if (emptyQ(neg) && emptyQ(pos)) {
2886 Polyhedron_Free(para);
2887 Polyhedron_Free(pos);
2888 Polyhedron_Free(neg);
2889 continue;
2890 }
2891 #ifdef DEBUG_ER
2892 fprintf(stderr, "\nER: Order\n");
2893 #endif /* DEBUG_ER */
2894 EP = barvinok_enumerate_e(para, exist, nparam, MaxRays);
2895 evalue *E;
2896 if (!emptyQ(pos)) {
2897 E = barvinok_enumerate_e(pos, exist, nparam, MaxRays);
2898 eadd(E, EP);
2899 free_evalue_refs(E);
2900 free(E);
2901 }
2902 if (!emptyQ(neg)) {
2903 E = barvinok_enumerate_e(neg, exist, nparam, MaxRays);
2904 eadd(E, EP);
2905 free_evalue_refs(E);
2906 free(E);
2907 }
2908 Polyhedron_Free(para);
2909 Polyhedron_Free(pos);
2910 Polyhedron_Free(neg);
2911 break;
2912 }
2913 if (EP)
2914 break;
2915 } END_FORALL_PVertex_in_ParamPolyhedron;
2916 if (EP)
2917 break;
2918 } END_FORALL_PVertex_in_ParamPolyhedron;
2919
2920 if (!EP) {
2921 /* Search for vertex coordinate to split on */
2922 /* First look for one independent of the parameters */
2923 FORALL_PVertex_in_ParamPolyhedron(V, last, PP) {
2924 for (int i = 0; i < exist; ++i) {
2925 int j;
2926 for (j = 0; j < nparam; ++j)
2927 if (value_notzero_p(V->Vertex->p[i][j]))
2928 break;
2929 if (j < nparam)
2930 continue;
2931 value_set_si(M->p[0][0], 1);
2932 Vector_Set(M->p[0]+1, 0, nvar+exist);
2933 Vector_Copy(V->Vertex->p[i],
2934 M->p[0] + 1 + nvar + exist, nparam+1);
2935 value_oppose(M->p[0][1+nvar+i],
2936 V->Vertex->p[i][nparam+1]);
2937
2938 Polyhedron *pos, *neg;
2939 value_set_si(M->p[0][0], 1);
2940 value_decrement(M->p[0][P->Dimension+1],
2941 M->p[0][P->Dimension+1]);
2942 neg = AddConstraints(M->p[0], 1, P, MaxRays);
2943 value_set_si(f, -1);
2944 Vector_Scale(M->p[0]+1, M->p[0]+1, f,
2945 P->Dimension+1);
2946 value_decrement(M->p[0][P->Dimension+1],
2947 M->p[0][P->Dimension+1]);
2948 value_decrement(M->p[0][P->Dimension+1],
2949 M->p[0][P->Dimension+1]);
2950 pos = AddConstraints(M->p[0], 1, P, MaxRays);
2951 if (emptyQ(neg) || emptyQ(pos)) {
2952 Polyhedron_Free(pos);
2953 Polyhedron_Free(neg);
2954 continue;
2955 }
2956 Polyhedron_Free(pos);
2957 value_increment(M->p[0][P->Dimension+1],
2958 M->p[0][P->Dimension+1]);
2959 pos = AddConstraints(M->p[0], 1, P, MaxRays);
2960 #ifdef DEBUG_ER
2961 fprintf(stderr, "\nER: Vertex\n");
2962 #endif /* DEBUG_ER */
2963 pos->next = neg;
2964 EP = enumerate_or(pos, exist, nparam, MaxRays);
2965 break;
2966 }
2967 if (EP)
2968 break;
2969 } END_FORALL_PVertex_in_ParamPolyhedron;
2970 }
2971
2972 if (!EP) {
2973 /* Search for vertex coordinate to split on */
2974 /* Now look for one that depends on the parameters */
2975 FORALL_PVertex_in_ParamPolyhedron(V, last, PP) {
2976 for (int i = 0; i < exist; ++i) {
2977 value_set_si(M->p[0][0], 1);
2978 Vector_Set(M->p[0]+1, 0, nvar+exist);
2979 Vector_Copy(V->Vertex->p[i],
2980 M->p[0] + 1 + nvar + exist, nparam+1);
2981 value_oppose(M->p[0][1+nvar+i],
2982 V->Vertex->p[i][nparam+1]);
2983
2984 Polyhedron *pos, *neg;
2985 value_set_si(M->p[0][0], 1);
2986 value_decrement(M->p[0][P->Dimension+1],
2987 M->p[0][P->Dimension+1]);
2988 neg = AddConstraints(M->p[0], 1, P, MaxRays);
2989 value_set_si(f, -1);
2990 Vector_Scale(M->p[0]+1, M->p[0]+1, f,
2991 P->Dimension+1);
2992 value_decrement(M->p[0][P->Dimension+1],
2993 M->p[0][P->Dimension+1]);
2994 value_decrement(M->p[0][P->Dimension+1],
2995 M->p[0][P->Dimension+1]);
2996 pos = AddConstraints(M->p[0], 1, P, MaxRays);
2997 if (emptyQ(neg) || emptyQ(pos)) {
2998 Polyhedron_Free(pos);
2999 Polyhedron_Free(neg);
3000 continue;
3001 }
3002 Polyhedron_Free(pos);
3003 value_increment(M->p[0][P->Dimension+1],
3004 M->p[0][P->Dimension+1]);
3005 pos = AddConstraints(M->p[0], 1, P, MaxRays);
3006 #ifdef DEBUG_ER
3007 fprintf(stderr, "\nER: ParamVertex\n");
3008 #endif /* DEBUG_ER */
3009 pos->next = neg;
3010 EP = enumerate_or(pos, exist, nparam, MaxRays);
3011 break;
3012 }
3013 if (EP)
3014 break;
3015 } END_FORALL_PVertex_in_ParamPolyhedron;
3016 }
3017
3018 Matrix_Free(M);
3019 value_clear(f);
3020 }
3021
3022 if (CEq)
3023 Polyhedron_Free(CEq);
3024 if (CT)
3025 Matrix_Free(CT);
3026 if (PP)
3027 Param_Polyhedron_Free(PP);
3028 *PA = P;
3029
3030 return EP;
3031 }
3032
3033 #ifndef HAVE_PIPLIB
3034 evalue *barvinok_enumerate_pip(Polyhedron *P,
3035 unsigned exist, unsigned nparam, unsigned MaxRays)
3036 {
3037 return 0;
3038 }
3039 #else
3040 evalue *barvinok_enumerate_pip(Polyhedron *P,
3041 unsigned exist, unsigned nparam, unsigned MaxRays)
3042 {
3043 int nvar = P->Dimension - exist - nparam;
3044 evalue *EP = evalue_zero();
3045 Polyhedron *Q, *N;
3046
3047 #ifdef DEBUG_ER
3048 fprintf(stderr, "\nER: PIP\n");
3049 #endif /* DEBUG_ER */
3050
3051 Polyhedron *D = pip_projectout(P, nvar, exist, nparam);
3052 for (Q = D; Q; Q = N) {
3053 N = Q->next;
3054 Q->next = 0;
3055 evalue *E;
3056 exist = Q->Dimension - nvar - nparam;
3057 E = barvinok_enumerate_e(Q, exist, nparam, MaxRays);
3058 Polyhedron_Free(Q);
3059 eadd(E, EP);
3060 free_evalue_refs(E);
3061 free(E);
3062 }
3063
3064 return EP;
3065 }
3066 #endif
3067
3068
3069 static bool is_single(Value *row, int pos, int len)
3070 {
3071 return First_Non_Zero(row, pos) == -1 &&
3072 First_Non_Zero(row+pos+1, len-pos-1) == -1;
3073 }
3074
3075 static evalue* barvinok_enumerate_e_r(Polyhedron *P,
3076 unsigned exist, unsigned nparam, unsigned MaxRays);
3077
3078 #ifdef DEBUG_ER
3079 static int er_level = 0;
3080
3081 evalue* barvinok_enumerate_e(Polyhedron *P,
3082 unsigned exist, unsigned nparam, unsigned MaxRays)
3083 {
3084 fprintf(stderr, "\nER: level %i\n", er_level);
3085
3086 Polyhedron_PrintConstraints(stderr, P_VALUE_FMT, P);
3087 fprintf(stderr, "\nE %d\nP %d\n", exist, nparam);
3088 ++er_level;
3089 P = DomainConstraintSimplify(Polyhedron_Copy(P), MaxRays);
3090 evalue *EP = barvinok_enumerate_e_r(P, exist, nparam, MaxRays);
3091 Polyhedron_Free(P);
3092 --er_level;
3093 return EP;
3094 }
3095 #else
3096 evalue* barvinok_enumerate_e(Polyhedron *P,
3097 unsigned exist, unsigned nparam, unsigned MaxRays)
3098 {
3099 P = DomainConstraintSimplify(Polyhedron_Copy(P), MaxRays);
3100 evalue *EP = barvinok_enumerate_e_r(P, exist, nparam, MaxRays);
3101 Polyhedron_Free(P);
3102 return EP;
3103 }
3104 #endif
3105
3106 static evalue* barvinok_enumerate_e_r(Polyhedron *P,
3107 unsigned exist, unsigned nparam, unsigned MaxRays)
3108 {
3109 if (exist == 0) {
3110 Polyhedron *U = Universe_Polyhedron(nparam);
3111 evalue *EP = barvinok_enumerate_ev(P, U, MaxRays);
3112 //char *param_name[] = {"P", "Q", "R", "S", "T" };
3113 //print_evalue(stdout, EP, param_name);
3114 Polyhedron_Free(U);
3115 return EP;
3116 }
3117
3118 int nvar = P->Dimension - exist - nparam;
3119 int len = P->Dimension + 2;
3120
3121 /* for now */
3122 POL_ENSURE_FACETS(P);
3123 POL_ENSURE_VERTICES(P);
3124
3125 if (emptyQ(P))
3126 return evalue_zero();
3127
3128 if (nvar == 0 && nparam == 0) {
3129 evalue *EP = evalue_zero();
3130 barvinok_count(P, &EP->x.n, MaxRays);
3131 if (value_pos_p(EP->x.n))
3132 value_set_si(EP->x.n, 1);
3133 return EP;
3134 }
3135
3136 int r;
3137 for (r = 0; r < P->NbRays; ++r)
3138 if (value_zero_p(P->Ray[r][0]) ||
3139 value_zero_p(P->Ray[r][P->Dimension+1])) {
3140 int i;
3141 for (i = 0; i < nvar; ++i)
3142 if (value_notzero_p(P->Ray[r][i+1]))
3143 break;
3144 if (i >= nvar)
3145 continue;
3146 for (i = nvar + exist; i < nvar + exist + nparam; ++i)
3147 if (value_notzero_p(P->Ray[r][i+1]))
3148 break;
3149 if (i >= nvar + exist + nparam)
3150 break;
3151 }
3152 if (r < P->NbRays) {
3153 evalue *EP = evalue_zero();
3154 value_set_si(EP->x.n, -1);
3155 return EP;
3156 }
3157
3158 int first;
3159 for (r = 0; r < P->NbEq; ++r)
3160 if ((first = First_Non_Zero(P->Constraint[r]+1+nvar, exist)) != -1)
3161 break;
3162 if (r < P->NbEq) {
3163 if (First_Non_Zero(P->Constraint[r]+1+nvar+first+1,
3164 exist-first-1) != -1) {
3165 Polyhedron *T = rotate_along(P, r, nvar, exist, MaxRays);
3166 #ifdef DEBUG_ER
3167 fprintf(stderr, "\nER: Equality\n");
3168 #endif /* DEBUG_ER */
3169 evalue *EP = barvinok_enumerate_e(T, exist-1, nparam, MaxRays);
3170 Polyhedron_Free(T);
3171 return EP;
3172 } else {
3173 #ifdef DEBUG_ER
3174 fprintf(stderr, "\nER: Fixed\n");
3175 #endif /* DEBUG_ER */
3176 if (first == 0)
3177 return barvinok_enumerate_e(P, exist-1, nparam, MaxRays);
3178 else {
3179 Polyhedron *T = Polyhedron_Copy(P);
3180 SwapColumns(T, nvar+1, nvar+1+first);
3181 evalue *EP = barvinok_enumerate_e(T, exist-1, nparam, MaxRays);
3182 Polyhedron_Free(T);
3183 return EP;
3184 }
3185 }
3186 }
3187
3188 Vector *row = Vector_Alloc(len);
3189 value_set_si(row->p[0], 1);
3190
3191 Value f;
3192 value_init(f);
3193
3194 enum constraint* info = new constraint[exist];
3195 for (int i = 0; i < exist; ++i) {
3196 info[i] = ALL_POS;
3197 for (int l = P->NbEq; l < P->NbConstraints; ++l) {
3198 if (value_negz_p(P->Constraint[l][nvar+i+1]))
3199 continue;
3200 bool l_parallel = is_single(P->Constraint[l]+nvar+1, i, exist);
3201 for (int u = P->NbEq; u < P->NbConstraints; ++u) {
3202 if (value_posz_p(P->Constraint[u][nvar+i+1]))
3203 continue;
3204 bool lu_parallel = l_parallel ||
3205 is_single(P->Constraint[u]+nvar+1, i, exist);
3206 value_oppose(f, P->Constraint[u][nvar+i+1]);
3207 Vector_Combine(P->Constraint[l]+1, P->Constraint[u]+1, row->p+1,
3208 f, P->Constraint[l][nvar+i+1], len-1);
3209 if (!(info[i] & INDEPENDENT)) {
3210 int j;
3211 for (j = 0; j < exist; ++j)
3212 if (j != i && value_notzero_p(row->p[nvar+j+1]))
3213 break;
3214 if (j == exist) {
3215 //printf("independent: i: %d, l: %d, u: %d\n", i, l, u);
3216 info[i] = (constraint)(info[i] | INDEPENDENT);
3217 }
3218 }
3219 if (info[i] & ALL_POS) {
3220 value_addto(row->p[len-1], row->p[len-1],
3221 P->Constraint[l][nvar+i+1]);
3222 value_addto(row->p[len-1], row->p[len-1], f);
3223 value_multiply(f, f, P->Constraint[l][nvar+i+1]);
3224 value_subtract(row->p[len-1], row->p[len-1], f);
3225 value_decrement(row->p[len-1], row->p[len-1]);
3226 ConstraintSimplify(row->p, row->p, len, &f);
3227 value_set_si(f, -1);
3228 Vector_Scale(row->p+1, row->p+1, f, len-1);
3229 value_decrement(row->p[len-1], row->p[len-1]);
3230 Polyhedron *T = AddConstraints(row->p, 1, P, MaxRays);
3231 if (!emptyQ(T)) {
3232 //printf("not all_pos: i: %d, l: %d, u: %d\n", i, l, u);
3233 info[i] = (constraint)(info[i] ^ ALL_POS);
3234 }
3235 //puts("pos remainder");
3236 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
3237 Polyhedron_Free(T);
3238 }
3239 if (!(info[i] & ONE_NEG)) {
3240 if (lu_parallel) {
3241 negative_test_constraint(P->Constraint[l],
3242 P->Constraint[u],
3243 row->p, nvar+i, len, &f);
3244 oppose_constraint(row->p, len, &f);
3245 Polyhedron *T = AddConstraints(row->p, 1, P, MaxRays);
3246 if (emptyQ(T)) {
3247 //printf("one_neg i: %d, l: %d, u: %d\n", i, l, u);
3248 info[i] = (constraint)(info[i] | ONE_NEG);
3249 }
3250 //puts("neg remainder");
3251 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
3252 Polyhedron_Free(T);
3253 } else if (!(info[i] & ROT_NEG)) {
3254 if (parallel_constraints(P->Constraint[l],
3255 P->Constraint[u],
3256 row->p, nvar, exist)) {
3257 negative_test_constraint7(P->Constraint[l],
3258 P->Constraint[u],
3259 row->p, nvar, exist,
3260 len, &f);
3261 oppose_constraint(row->p, len, &f);
3262 Polyhedron *T = AddConstraints(row->p, 1, P, MaxRays);
3263 if (emptyQ(T)) {
3264 // printf("rot_neg i: %d, l: %d, u: %d\n", i, l, u);
3265 info[i] = (constraint)(info[i] | ROT_NEG);
3266 r = l;
3267 }
3268 //puts("neg remainder");
3269 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
3270 Polyhedron_Free(T);
3271 }
3272 }
3273 }
3274 if (!(info[i] & ALL_POS) && (info[i] & (ONE_NEG | ROT_NEG)))
3275 goto next;
3276 }
3277 }
3278 if (info[i] & ALL_POS)
3279 break;
3280 next:
3281 ;
3282 }
3283
3284 /*
3285 for (int i = 0; i < exist; ++i)
3286 printf("%i: %i\n", i, info[i]);
3287 */
3288 for (int i = 0; i < exist; ++i)
3289 if (info[i] & ALL_POS) {
3290 #ifdef DEBUG_ER
3291 fprintf(stderr, "\nER: Positive\n");
3292 #endif /* DEBUG_ER */
3293 // Eliminate
3294 // Maybe we should chew off some of the fat here
3295 Matrix *M = Matrix_Alloc(P->Dimension, P->Dimension+1);
3296 for (int j = 0; j < P->Dimension; ++j)
3297 value_set_si(M->p[j][j + (j >= i+nvar)], 1);
3298 Polyhedron *T = Polyhedron_Image(P, M, MaxRays);
3299 Matrix_Free(M);
3300 evalue *EP = barvinok_enumerate_e(T, exist-1, nparam, MaxRays);
3301 Polyhedron_Free(T);
3302 value_clear(f);
3303 Vector_Free(row);
3304 delete [] info;
3305 return EP;
3306 }
3307 for (int i = 0; i < exist; ++i)
3308 if (info[i] & ONE_NEG) {
3309 #ifdef DEBUG_ER
3310 fprintf(stderr, "\nER: Negative\n");
3311 #endif /* DEBUG_ER */
3312 Vector_Free(row);
3313 value_clear(f);
3314 delete [] info;
3315 if (i == 0)
3316 return barvinok_enumerate_e(P, exist-1, nparam, MaxRays);
3317 else {
3318 Polyhedron *T = Polyhedron_Copy(P);
3319 SwapColumns(T, nvar+1, nvar+1+i);
3320 evalue *EP = barvinok_enumerate_e(T, exist-1, nparam, MaxRays);
3321 Polyhedron_Free(T);
3322 return EP;
3323 }
3324 }
3325 for (int i = 0; i < exist; ++i)
3326 if (info[i] & ROT_NEG) {
3327 #ifdef DEBUG_ER
3328 fprintf(stderr, "\nER: Rotate\n");
3329 #endif /* DEBUG_ER */
3330 Vector_Free(row);
3331 value_clear(f);
3332 delete [] info;
3333 Polyhedron *T = rotate_along(P, r, nvar, exist, MaxRays);
3334 evalue *EP = barvinok_enumerate_e(T, exist-1, nparam, MaxRays);
3335 Polyhedron_Free(T);
3336 return EP;
3337 }
3338 for (int i = 0; i < exist; ++i)
3339 if (info[i] & INDEPENDENT) {
3340 Polyhedron *pos, *neg;
3341
3342 /* Find constraint again and split off negative part */
3343
3344 if (SplitOnVar(P, i, nvar, exist, MaxRays,
3345 row, f, true, &pos, &neg)) {
3346 #ifdef DEBUG_ER
3347 fprintf(stderr, "\nER: Split\n");
3348 #endif /* DEBUG_ER */
3349
3350 evalue *EP =
3351 barvinok_enumerate_e(neg, exist-1, nparam, MaxRays);
3352 evalue *E =
3353 barvinok_enumerate_e(pos, exist, nparam, MaxRays);
3354 eadd(E, EP);
3355 free_evalue_refs(E);
3356 free(E);
3357 Polyhedron_Free(neg);
3358 Polyhedron_Free(pos);
3359 value_clear(f);
3360 Vector_Free(row);
3361 delete [] info;
3362 return EP;
3363 }
3364 }
3365 delete [] info;
3366
3367 Polyhedron *O = P;
3368 Polyhedron *F;
3369
3370 evalue *EP;
3371
3372 EP = enumerate_line(P, exist, nparam, MaxRays);
3373 if (EP)
3374 goto out;
3375
3376 EP = barvinok_enumerate_pip(P, exist, nparam, MaxRays);
3377 if (EP)
3378 goto out;
3379
3380 EP = enumerate_redundant_ray(P, exist, nparam, MaxRays);
3381 if (EP)
3382 goto out;
3383
3384 EP = enumerate_sure(P, exist, nparam, MaxRays);
3385 if (EP)
3386 goto out;
3387
3388 EP = enumerate_ray(P, exist, nparam, MaxRays);
3389 if (EP)
3390 goto out;
3391
3392 EP = enumerate_sure2(P, exist, nparam, MaxRays);
3393 if (EP)
3394 goto out;
3395
3396 F = unfringe(P, MaxRays);
3397 if (!PolyhedronIncludes(F, P)) {
3398 #ifdef DEBUG_ER
3399 fprintf(stderr, "\nER: Fringed\n");
3400 #endif /* DEBUG_ER */
3401 EP = barvinok_enumerate_e(F, exist, nparam, MaxRays);
3402 Polyhedron_Free(F);
3403 goto out;
3404 }
3405 Polyhedron_Free(F);
3406
3407 if (nparam)
3408 EP = enumerate_vd(&P, exist, nparam, MaxRays);
3409 if (EP)
3410 goto out2;
3411
3412 if (nvar != 0) {
3413 EP = enumerate_sum(P, exist, nparam, MaxRays);
3414 goto out2;
3415 }
3416
3417 assert(nvar == 0);
3418
3419 int i;
3420 Polyhedron *pos, *neg;
3421 for (i = 0; i < exist; ++i)
3422 if (SplitOnVar(P, i, nvar, exist, MaxRays,
3423 row, f, false, &pos, &neg))
3424 break;
3425
3426 assert (i < exist);
3427
3428 pos->next = neg;
3429 EP = enumerate_or(pos, exist, nparam, MaxRays);
3430
3431 out2:
3432 if (O != P)
3433 Polyhedron_Free(P);
3434
3435 out:
3436 value_clear(f);
3437 Vector_Free(row);
3438 return EP;
3439 }
3440
3441 /*
3442 * remove equalities that require a "compression" of the parameters
3443 */
3444 #ifndef HAVE_COMPRESS_PARMS
3445 static Polyhedron *remove_more_equalities(Polyhedron *P, unsigned nparam,
3446 Matrix **CP, unsigned MaxRays)
3447 {
3448 return P;
3449 }
3450 #else
3451 static Polyhedron *remove_more_equalities(Polyhedron *P, unsigned nparam,
3452 Matrix **CP, unsigned MaxRays)
3453 {
3454 Matrix *M, *T;
3455 Polyhedron *Q;
3456 Matrix *CV = NULL;
3457 int i;
3458
3459 /* compress_parms doesn't like equalities that only involve parameters */
3460 for (i = 0; i < P->NbEq; ++i)
3461 if (First_Non_Zero(P->Constraint[i]+1, P->Dimension-nparam) == -1)
3462 break;
3463
3464 if (i < P->NbEq) {
3465 Matrix *M = Matrix_Alloc(P->NbEq, 1+nparam+1);
3466 int n = 0;
3467 for (; i < P->NbEq; ++i) {
3468 if (First_Non_Zero(P->Constraint[i]+1, P->Dimension-nparam) == -1)
3469 Vector_Copy(P->Constraint[i]+1+P->Dimension-nparam,
3470 M->p[n++]+1, nparam+1);
3471 }
3472 M->NbRows = n;
3473 CV = compress_variables(M, 0);
3474 T = align_matrix(CV, P->Dimension+1);
3475 Q = Polyhedron_Preimage(P, T, MaxRays);
3476 Matrix_Free(T);
3477 Polyhedron_Free(P);
3478 P = Q;
3479 Matrix_Free(M);
3480 nparam = CV->NbColumns-1;
3481 }
3482
3483 if (P->NbEq == 0) {
3484 *CP = CV;
3485 return P;
3486 }
3487
3488 M = Matrix_Alloc(P->NbEq, P->Dimension+2);
3489 Vector_Copy(P->Constraint[0], M->p[0], P->NbEq * (P->Dimension+2));
3490 *CP = compress_parms(M, nparam);
3491 T = align_matrix(*CP, P->Dimension+1);
3492 Q = Polyhedron_Preimage(P, T, MaxRays);
3493 Polyhedron_Free(P);
3494 P = Q;
3495 P = remove_equalities_p(P, P->Dimension-nparam, NULL);
3496 Matrix_Free(T);
3497 Matrix_Free(M);
3498
3499 if (CV) {
3500 T = *CP;
3501 *CP = Matrix_Alloc(CV->NbRows, T->NbColumns);
3502 Matrix_Product(CV, T, *CP);
3503 Matrix_Free(T);
3504 Matrix_Free(CV);
3505 }
3506
3507 return P;
3508 }
3509 #endif
3510
3511 gen_fun * barvinok_series(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
3512 {
3513 Matrix *CP = NULL;
3514 Polyhedron *CA;
3515 unsigned nparam = C->Dimension;
3516 gen_fun *gf;
3517
3518 CA = align_context(C, P->Dimension, MaxRays);
3519 P = DomainIntersection(P, CA, MaxRays);
3520 Polyhedron_Free(CA);
3521
3522 if (emptyQ2(P)) {
3523 Polyhedron_Free(P);
3524 return new gen_fun;
3525 }
3526
3527 assert(!Polyhedron_is_infinite_param(P, nparam));
3528 assert(P->NbBid == 0);
3529 assert(Polyhedron_has_positive_rays(P, nparam));
3530 if (P->NbEq != 0)
3531 P = remove_equalities_p(P, P->Dimension-nparam, NULL);
3532 if (P->NbEq != 0)
3533 P = remove_more_equalities(P, nparam, &CP, MaxRays);
3534 assert(P->NbEq == 0);
3535 if (CP)
3536 nparam = CP->NbColumns-1;
3537
3538 gf_base *red;
3539 red = gf_base::create(Polyhedron_Project(P, nparam), P->Dimension, nparam);
3540 red->start_gf(P, MaxRays);
3541 Polyhedron_Free(P);
3542 if (CP) {
3543 red->gf->substitute(CP);
3544 Matrix_Free(CP);
3545 }
3546 gf = red->gf;
3547 delete red;
3548 return gf;
3549 }
3550
3551 static Polyhedron *skew_into_positive_orthant(Polyhedron *D, unsigned nparam,
3552 unsigned MaxRays)
3553 {
3554 Matrix *M = NULL;
3555 Value tmp;
3556 value_init(tmp);
3557 for (Polyhedron *P = D; P; P = P->next) {
3558 POL_ENSURE_VERTICES(P);
3559 assert(!Polyhedron_is_infinite_param(P, nparam));
3560 assert(P->NbBid == 0);
3561 assert(Polyhedron_has_positive_rays(P, nparam));
3562
3563 for (int r = 0; r < P->NbRays; ++r) {
3564 if (value_notzero_p(P->Ray[r][P->Dimension+1]))
3565 continue;
3566 for (int i = 0; i < nparam; ++i) {
3567 int j;
3568 if (value_posz_p(P->Ray[r][i+1]))
3569 continue;
3570 if (!M) {
3571 M = Matrix_Alloc(D->Dimension+1, D->Dimension+1);
3572 for (int i = 0; i < D->Dimension+1; ++i)
3573 value_set_si(M->p[i][i], 1);
3574 } else {
3575 Inner_Product(P->Ray[r]+1, M->p[i], D->Dimension+1, &tmp);
3576 if (value_posz_p(tmp))
3577 continue;
3578 }
3579 for (j = P->Dimension - nparam; j < P->Dimension; ++j)
3580 if (value_pos_p(P->Ray[r][j+1]))
3581 break;
3582 assert(j < P->Dimension);
3583 value_pdivision(tmp, P->Ray[r][j+1], P->Ray[r][i+1]);
3584 value_subtract(M->p[i][j], M->p[i][j], tmp);
3585 }
3586 }
3587 }
3588 value_clear(tmp);
3589 if (M) {
3590 D = DomainImage(D, M, MaxRays);
3591 Matrix_Free(M);
3592 }
3593 return D;
3594 }
3595
3596 gen_fun* barvinok_enumerate_union_series(Polyhedron *D, Polyhedron* C,
3597 unsigned MaxRays)
3598 {
3599 Polyhedron *conv, *D2;
3600 gen_fun *gf = NULL, *gf2;
3601 unsigned nparam = C->Dimension;
3602 ZZ one, mone;
3603 one = 1;
3604 mone = -1;
3605 D2 = skew_into_positive_orthant(D, nparam, MaxRays);
3606 for (Polyhedron *P = D2; P; P = P->next) {
3607 assert(P->Dimension == D2->Dimension);
3608 POL_ENSURE_VERTICES(P);
3609 /* it doesn't matter which reducer we use, since we don't actually
3610 * reduce anything here
3611 */
3612 partial_reducer red(Polyhedron_Project(P, P->Dimension), P->Dimension,
3613 P->Dimension);
3614 red.start(P, MaxRays);
3615 if (!gf)
3616 gf = red.gf;
3617 else {
3618 gf->add_union(red.gf, MaxRays);
3619 delete red.gf;
3620 }
3621 }
3622 /* we actually only need the convex union of the parameter space
3623 * but the reducer classes currently expect a polyhedron in
3624 * the combined space
3625 */
3626 Polyhedron_Free(gf->context);
3627 gf->context = DomainConvex(D2, MaxRays);
3628
3629 gf2 = gf->summate(D2->Dimension - nparam);
3630
3631 delete gf;
3632 if (D != D2)
3633 Domain_Free(D2);
3634 return gf2;
3635 }
3636
3637 evalue* barvinok_enumerate_union(Polyhedron *D, Polyhedron* C, unsigned MaxRays)
3638 {
3639 evalue *EP;
3640 gen_fun *gf = barvinok_enumerate_union_series(D, C, MaxRays);
3641 EP = *gf;
3642 delete gf;
3643 return EP;
3644 }
lattice point enumeration library