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