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