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