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