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