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