update isl for fix in parameter alignment
[barvinok.git] / barvinok.cc
blobb1f19bc6fe5799ea21df865e727a52572182995c
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 <isl_map_polylib.h>
11 #include <barvinok/util.h>
12 #include <barvinok/evalue.h>
13 #include "config.h"
14 #include <barvinok/barvinok.h>
15 #include <barvinok/genfun.h>
16 #include <barvinok/options.h>
17 #include <barvinok/sample.h>
18 #include "bfcounter.h"
19 #include "conversion.h"
20 #include "counter.h"
21 #include "decomposer.h"
22 #include "euler.h"
23 #include "lattice_point.h"
24 #include "laurent.h"
25 #include "reduce_domain.h"
26 #include "remove_equalities.h"
27 #include "scale.h"
28 #include "volume.h"
29 #include "bernoulli.h"
30 #include "param_util.h"
31 #include "summate.h"
33 #ifdef NTL_STD_CXX
34 using namespace NTL;
35 #endif
36 using std::cerr;
37 using std::cout;
38 using std::endl;
39 using std::vector;
40 using std::deque;
41 using std::string;
42 using std::ostringstream;
44 #define ALLOC(t,p) p = (t*)malloc(sizeof(*p))
46 class dpoly_n {
47 public:
48 Matrix *coeff;
49 ~dpoly_n() {
50 Matrix_Free(coeff);
52 dpoly_n(int d) {
53 Value d0, one;
54 value_init(d0);
55 value_init(one);
56 value_set_si(one, 1);
57 coeff = Matrix_Alloc(d+1, d+1+1);
58 value_set_si(coeff->p[0][0], 1);
59 value_set_si(coeff->p[0][d+1], 1);
60 for (int i = 1; i <= d; ++i) {
61 value_multiply(coeff->p[i][0], coeff->p[i-1][0], d0);
62 Vector_Combine(coeff->p[i-1], coeff->p[i-1]+1, coeff->p[i]+1,
63 one, d0, i);
64 value_set_si(coeff->p[i][d+1], i);
65 value_multiply(coeff->p[i][d+1], coeff->p[i][d+1], coeff->p[i-1][d+1]);
66 value_decrement(d0, d0);
68 value_clear(d0);
69 value_clear(one);
71 void div(dpoly& d, Vector *count, int sign) {
72 int len = coeff->NbRows;
73 Matrix * c = Matrix_Alloc(coeff->NbRows, coeff->NbColumns);
74 Value tmp;
75 value_init(tmp);
76 for (int i = 0; i < len; ++i) {
77 Vector_Copy(coeff->p[i], c->p[i], len+1);
78 for (int j = 1; j <= i; ++j) {
79 value_multiply(tmp, d.coeff->p[j], c->p[i][len]);
80 value_oppose(tmp, tmp);
81 Vector_Combine(c->p[i], c->p[i-j], c->p[i],
82 c->p[i-j][len], tmp, len);
83 value_multiply(c->p[i][len], c->p[i][len], c->p[i-j][len]);
85 value_multiply(c->p[i][len], c->p[i][len], d.coeff->p[0]);
87 if (sign == -1) {
88 value_set_si(tmp, -1);
89 Vector_Scale(c->p[len-1], count->p, tmp, len);
90 value_assign(count->p[len], c->p[len-1][len]);
91 } else
92 Vector_Copy(c->p[len-1], count->p, len+1);
93 Vector_Normalize(count->p, len+1);
94 value_clear(tmp);
95 Matrix_Free(c);
99 static void add_rays(mat_ZZ& rays, Polyhedron *i, int *r, int nvar = -1,
100 bool all = false)
102 unsigned dim = i->Dimension;
103 if (nvar == -1)
104 nvar = dim;
105 for (int k = 0; k < i->NbRays; ++k) {
106 if (!value_zero_p(i->Ray[k][dim+1]))
107 continue;
108 if (!all && nvar != dim && First_Non_Zero(i->Ray[k]+1, nvar) == -1)
109 continue;
110 values2zz(i->Ray[k]+1, rays[(*r)++], nvar);
114 struct bfe_term : public bfc_term_base {
115 vector<evalue *> factors;
117 bfe_term(int len) : bfc_term_base(len) {
120 ~bfe_term() {
121 for (int i = 0; i < factors.size(); ++i) {
122 if (!factors[i])
123 continue;
124 free_evalue_refs(factors[i]);
125 delete factors[i];
130 static void print_int_vector(int *v, int len, const char *name)
132 cerr << name << endl;
133 for (int j = 0; j < len; ++j) {
134 cerr << v[j] << " ";
136 cerr << endl;
139 static void print_bfc_terms(mat_ZZ& factors, bfc_vec& v)
141 cerr << endl;
142 cerr << "factors" << endl;
143 cerr << factors << endl;
144 for (int i = 0; i < v.size(); ++i) {
145 cerr << "term: " << i << endl;
146 print_int_vector(v[i]->powers, factors.NumRows(), "powers");
147 cerr << "terms" << endl;
148 cerr << v[i]->terms << endl;
149 bfc_term* bfct = static_cast<bfc_term *>(v[i]);
150 cerr << bfct->c << endl;
154 static void print_bfe_terms(mat_ZZ& factors, bfc_vec& v)
156 cerr << endl;
157 cerr << "factors" << endl;
158 cerr << factors << endl;
159 for (int i = 0; i < v.size(); ++i) {
160 cerr << "term: " << i << endl;
161 print_int_vector(v[i]->powers, factors.NumRows(), "powers");
162 cerr << "terms" << endl;
163 cerr << v[i]->terms << endl;
164 bfe_term* bfet = static_cast<bfe_term *>(v[i]);
165 for (int j = 0; j < v[i]->terms.NumRows(); ++j) {
166 const char * test[] = {"a", "b"};
167 print_evalue(stderr, bfet->factors[j], test);
168 fprintf(stderr, "\n");
173 struct bfcounter : public bfcounter_base {
174 mpq_t count;
175 Value tz;
177 bfcounter(unsigned dim) : bfcounter_base(dim) {
178 mpq_init(count);
179 lower = 1;
180 value_init(tz);
182 ~bfcounter() {
183 mpq_clear(count);
184 value_clear(tz);
186 virtual void base(mat_ZZ& factors, bfc_vec& v);
187 virtual void get_count(Value *result) {
188 assert(value_one_p(&count[0]._mp_den));
189 value_assign(*result, &count[0]._mp_num);
193 void bfcounter::base(mat_ZZ& factors, bfc_vec& v)
195 unsigned nf = factors.NumRows();
197 for (int i = 0; i < v.size(); ++i) {
198 bfc_term* bfct = static_cast<bfc_term *>(v[i]);
199 int total_power = 0;
200 // factor is always positive, so we always
201 // change signs
202 for (int k = 0; k < nf; ++k)
203 total_power += v[i]->powers[k];
205 int j;
206 for (j = 0; j < nf; ++j)
207 if (v[i]->powers[j] > 0)
208 break;
210 zz2value(factors[j][0], tz);
211 dpoly D(total_power, tz, 1);
212 for (int k = 1; k < v[i]->powers[j]; ++k) {
213 zz2value(factors[j][0], tz);
214 dpoly fact(total_power, tz, 1);
215 D *= fact;
217 for ( ; ++j < nf; )
218 for (int k = 0; k < v[i]->powers[j]; ++k) {
219 zz2value(factors[j][0], tz);
220 dpoly fact(total_power, tz, 1);
221 D *= fact;
224 for (int k = 0; k < v[i]->terms.NumRows(); ++k) {
225 zz2value(v[i]->terms[k][0], tz);
226 dpoly n(total_power, tz);
227 mpq_set_si(tcount, 0, 1);
228 n.div(D, tcount, 1);
229 if (total_power % 2)
230 bfct->c[k].n = -bfct->c[k].n;
231 zz2value(bfct->c[k].n, tn);
232 zz2value(bfct->c[k].d, td);
234 mpz_mul(mpq_numref(tcount), mpq_numref(tcount), tn);
235 mpz_mul(mpq_denref(tcount), mpq_denref(tcount), td);
236 mpq_canonicalize(tcount);
237 mpq_add(count, count, tcount);
239 delete v[i];
244 /* Check whether the polyhedron is unbounded and if so,
245 * check whether it has any (and therefore an infinite number of)
246 * integer points.
247 * If one of the vertices is integer, then we are done.
248 * Otherwise, transform the polyhedron such that one of the rays
249 * is the first unit vector and cut it off at a height that ensures
250 * that if the whole polyhedron has any points, then the remaining part
251 * has integer points. In particular we add the largest coefficient
252 * of a ray to the highest vertex (rounded up).
254 static bool Polyhedron_is_infinite(Polyhedron *P, Value* result,
255 barvinok_options *options)
257 int r = 0;
258 Matrix *M, *M2;
259 Value c, tmp;
260 Value g;
261 bool first;
262 Vector *v;
263 Value offset, size;
264 Polyhedron *R;
266 if (P->NbBid == 0)
267 for (; r < P->NbRays; ++r)
268 if (value_zero_p(P->Ray[r][P->Dimension+1]))
269 break;
270 if (P->NbBid == 0 && r == P->NbRays)
271 return false;
273 if (options->count_sample_infinite) {
274 Vector *sample;
276 sample = Polyhedron_Sample(P, options);
277 if (!sample)
278 value_set_si(*result, 0);
279 else {
280 value_set_si(*result, -1);
281 Vector_Free(sample);
283 return true;
286 for (int i = 0; i < P->NbRays; ++i)
287 if (value_one_p(P->Ray[i][1+P->Dimension])) {
288 value_set_si(*result, -1);
289 return true;
292 value_init(g);
293 M = Matrix_Alloc(P->Dimension+1, P->Dimension+1);
294 Vector_Gcd(P->Ray[r]+1, P->Dimension, &g);
295 Vector_AntiScale(P->Ray[r]+1, M->p[0], g, P->Dimension+1);
296 int ok = unimodular_complete(M, 1);
297 assert(ok);
298 value_set_si(M->p[P->Dimension][P->Dimension], 1);
299 M2 = Transpose(M);
300 Matrix_Free(M);
301 P = Polyhedron_Preimage(P, M2, 0);
302 Matrix_Free(M2);
303 value_clear(g);
305 first = true;
306 value_init(offset);
307 value_init(size);
308 value_init(tmp);
309 value_set_si(size, 0);
311 for (int i = 0; i < P->NbBid; ++i) {
312 value_absolute(tmp, P->Ray[i][1]);
313 if (value_gt(tmp, size))
314 value_assign(size, tmp);
316 for (int i = P->NbBid; i < P->NbRays; ++i) {
317 if (value_zero_p(P->Ray[i][P->Dimension+1])) {
318 if (value_gt(P->Ray[i][1], size))
319 value_assign(size, P->Ray[i][1]);
320 continue;
322 mpz_cdiv_q(tmp, P->Ray[i][1], P->Ray[i][P->Dimension+1]);
323 if (first || value_gt(tmp, offset)) {
324 value_assign(offset, tmp);
325 first = false;
328 value_addto(offset, offset, size);
329 value_clear(size);
330 value_clear(tmp);
332 v = Vector_Alloc(P->Dimension+2);
333 value_set_si(v->p[0], 1);
334 value_set_si(v->p[1], -1);
335 value_assign(v->p[1+P->Dimension], offset);
336 R = AddConstraints(v->p, 1, P, options->MaxRays);
337 Polyhedron_Free(P);
338 P = R;
340 value_clear(offset);
341 Vector_Free(v);
343 value_init(c);
344 barvinok_count_with_options(P, &c, options);
345 Polyhedron_Free(P);
346 if (value_zero_p(c))
347 value_set_si(*result, 0);
348 else
349 value_set_si(*result, -1);
350 value_clear(c);
352 return true;
355 static void evalue2value(evalue *e, Value *v)
357 if (EVALUE_IS_ZERO(*e)) {
358 value_set_si(*v, 0);
359 return;
362 if (value_notzero_p(e->d)) {
363 assert(value_one_p(e->d));
364 value_assign(*v, e->x.n);
365 return;
368 assert(e->x.p->type == partition);
369 assert(e->x.p->size == 2);
370 assert(EVALUE_DOMAIN(e->x.p->arr[0])->Dimension == 0);
371 evalue2value(&e->x.p->arr[1], v);
374 static void barvinok_count_f(Polyhedron *P, Value* result,
375 barvinok_options *options);
377 void barvinok_count_with_options(Polyhedron *P, Value* result,
378 struct barvinok_options *options)
380 unsigned dim;
381 int allocated = 0;
382 Polyhedron *Q;
383 bool infinite = false;
385 if (P->next)
386 fprintf(stderr,
387 "barvinok_count: input is a union; only first polyhedron is counted\n");
389 if (emptyQ2(P)) {
390 value_set_si(*result, 0);
391 return;
393 if (P->NbEq != 0) {
394 Q = NULL;
395 do {
396 P = remove_equalities(P, options->MaxRays);
397 if (P)
398 P = DomainConstraintSimplify(P, options->MaxRays);
399 if (Q)
400 Polyhedron_Free(Q);
401 Q = P;
402 } while (P && !emptyQ(P) && P->NbEq != 0);
403 if (!P || emptyQ(P)) {
404 Polyhedron_Free(P);
405 value_set_si(*result, 0);
406 return;
408 allocated = 1;
410 if (Polyhedron_is_infinite(P, result, options)) {
411 if (allocated)
412 Polyhedron_Free(P);
413 return;
415 if (P->Dimension == 0) {
416 /* Test whether the constraints are satisfied */
417 POL_ENSURE_VERTICES(P);
418 value_set_si(*result, !emptyQ(P));
419 if (allocated)
420 Polyhedron_Free(P);
421 return;
423 if (options->summation == BV_SUM_BERNOULLI) {
424 Polyhedron *C = Universe_Polyhedron(0);
425 evalue *sum = barvinok_summate_unweighted(P, C, options);
426 Polyhedron_Free(C);
427 evalue2value(sum, result);
428 evalue_free(sum);
429 return;
431 Q = Polyhedron_Factor(P, 0, NULL, options->MaxRays);
432 if (Q) {
433 if (allocated)
434 Polyhedron_Free(P);
435 P = Q;
436 allocated = 1;
439 barvinok_count_f(P, result, options);
440 if (value_neg_p(*result))
441 infinite = true;
442 if (Q && P->next && value_notzero_p(*result)) {
443 Value factor;
444 value_init(factor);
446 for (Q = P->next; Q; Q = Q->next) {
447 barvinok_count_f(Q, &factor, options);
448 if (value_neg_p(factor)) {
449 infinite = true;
450 continue;
451 } else if (Q->next && value_zero_p(factor)) {
452 value_set_si(*result, 0);
453 break;
455 value_multiply(*result, *result, factor);
458 value_clear(factor);
461 if (allocated)
462 Domain_Free(P);
463 if (infinite)
464 value_set_si(*result, -1);
467 void barvinok_count(Polyhedron *P, Value* result, unsigned NbMaxCons)
469 barvinok_options *options = barvinok_options_new_with_defaults();
470 options->MaxRays = NbMaxCons;
471 barvinok_count_with_options(P, result, options);
472 barvinok_options_free(options);
475 static void barvinok_count_f(Polyhedron *P, Value* result,
476 barvinok_options *options)
478 if (emptyQ2(P)) {
479 value_set_si(*result, 0);
480 return;
483 if (P->Dimension == 1)
484 return Line_Length(P, result);
486 int c = P->NbConstraints;
487 POL_ENSURE_FACETS(P);
488 if (c != P->NbConstraints || P->NbEq != 0) {
489 Polyhedron *next = P->next;
490 P->next = NULL;
491 barvinok_count_with_options(P, result, options);
492 P->next = next;
493 return;
496 POL_ENSURE_VERTICES(P);
498 if (Polyhedron_is_infinite(P, result, options))
499 return;
501 np_base *cnt;
502 if (options->incremental_specialization == BV_SPECIALIZATION_BF)
503 cnt = new bfcounter(P->Dimension);
504 else if (options->incremental_specialization == BV_SPECIALIZATION_DF)
505 cnt = new icounter(P->Dimension);
506 else if (options->incremental_specialization == BV_SPECIALIZATION_TODD)
507 cnt = new tcounter(P->Dimension, options->max_index);
508 else
509 cnt = new counter(P->Dimension, options->max_index);
510 cnt->start(P, options);
512 cnt->get_count(result);
513 delete cnt;
516 static void uni_polynom(int param, Vector *c, evalue *EP)
518 unsigned dim = c->Size-2;
519 value_init(EP->d);
520 value_set_si(EP->d,0);
521 EP->x.p = new_enode(polynomial, dim+1, param+1);
522 for (int j = 0; j <= dim; ++j)
523 evalue_set(&EP->x.p->arr[j], c->p[j], c->p[dim+1]);
526 typedef evalue * evalue_p;
528 struct enumerator_base {
529 unsigned dim;
530 evalue ** vE;
531 evalue mone;
532 vertex_decomposer *vpd;
534 enumerator_base(unsigned dim, vertex_decomposer *vpd)
536 this->dim = dim;
537 this->vpd = vpd;
539 vE = new evalue_p[vpd->PP->nbV];
540 for (int j = 0; j < vpd->PP->nbV; ++j)
541 vE[j] = 0;
543 value_init(mone.d);
544 evalue_set_si(&mone, -1, 1);
547 void decompose_at(Param_Vertices *V, int _i, barvinok_options *options) {
548 //this->pVD = pVD;
550 vE[_i] = new evalue;
551 value_init(vE[_i]->d);
552 evalue_set_si(vE[_i], 0, 1);
554 vpd->decompose_at_vertex(V, _i, options);
557 virtual ~enumerator_base() {
558 for (int j = 0; j < vpd->PP->nbV; ++j)
559 if (vE[j]) {
560 free_evalue_refs(vE[j]);
561 delete vE[j];
563 delete [] vE;
565 free_evalue_refs(&mone);
568 static enumerator_base *create(Polyhedron *P, unsigned dim,
569 Param_Polyhedron *PP,
570 barvinok_options *options);
573 struct enumerator : public signed_cone_consumer, public vertex_decomposer,
574 public enumerator_base {
575 vec_ZZ lambda;
576 vec_ZZ den;
577 term_info num;
578 Vector *c;
579 mpq_t count;
580 Value tz;
582 enumerator(Polyhedron *P, unsigned dim, Param_Polyhedron *PP) :
583 vertex_decomposer(PP, *this), enumerator_base(dim, this) {
584 randomvector(P, lambda, dim, 0);
585 den.SetLength(dim);
586 c = Vector_Alloc(dim+2);
588 mpq_init(count);
589 value_init(tz);
592 ~enumerator() {
593 mpq_clear(count);
594 Vector_Free(c);
595 value_clear(tz);
598 virtual void handle(const signed_cone& sc, barvinok_options *options);
601 void enumerator::handle(const signed_cone& sc, barvinok_options *options)
603 int sign = sc.sign;
604 int r = 0;
605 assert(sc.rays.NumRows() == dim);
606 for (int k = 0; k < dim; ++k) {
607 if (lambda * sc.rays[k] == 0)
608 throw Orthogonal;
611 lattice_point(V, sc.rays, lambda, &num, sc.det, options);
612 den = sc.rays * lambda;
614 if (dim % 2)
615 sign = -sign;
617 zz2value(den[0], tz);
618 dpoly n(dim, tz, 1);
619 for (int k = 1; k < dim; ++k) {
620 zz2value(den[k], tz);
621 dpoly fact(dim, tz, 1);
622 n *= fact;
624 if (num.E != NULL) {
625 dpoly_n d(dim);
626 d.div(n, c, sign);
627 for (unsigned long i = 0; i < sc.det; ++i) {
628 evalue *EV = evalue_polynomial(c, num.E[i]);
629 eadd(EV, vE[vert]);
630 evalue_free(EV);
631 evalue_free(num.E[i]);
633 delete [] num.E;
634 } else {
635 mpq_set_si(count, 0, 1);
636 if (num.constant.length() == 1) {
637 zz2value(num.constant[0], tz);
638 dpoly d(dim, tz);
639 d.div(n, count, sign);
640 } else {
641 dpoly_n d(dim);
642 d.div(n, c, sign);
643 Value x, sum, acc;
644 value_init(x);
645 value_init(acc);
646 for (unsigned long i = 0; i < sc.det; ++i) {
647 value_assign(acc, c->p[dim]);
648 zz2value(num.constant[i], x);
649 for (int j = dim-1; j >= 0; --j) {
650 value_multiply(acc, acc, x);
651 value_addto(acc, acc, c->p[j]);
653 value_addto(mpq_numref(count), mpq_numref(count), acc);
655 mpz_set(mpq_denref(count), c->p[dim+1]);
656 value_clear(acc);
657 value_clear(x);
659 evalue EV;
660 value_init(EV.d);
661 evalue_set(&EV, &count[0]._mp_num, &count[0]._mp_den);
662 eadd(&EV, vE[vert]);
663 free_evalue_refs(&EV);
667 struct ienumerator_base : enumerator_base {
668 evalue ** E_vertex;
670 ienumerator_base(unsigned dim, vertex_decomposer *vpd) :
671 enumerator_base(dim,vpd) {
672 E_vertex = new evalue_p[dim];
675 virtual ~ienumerator_base() {
676 delete [] E_vertex;
679 evalue *E_num(int i, int d) {
680 return E_vertex[i + (dim-d)];
684 struct cumulator {
685 evalue *factor;
686 evalue *v;
687 dpoly_r *r;
689 cumulator(evalue *factor, evalue *v, dpoly_r *r) :
690 factor(factor), v(v), r(r) {}
692 void cumulate(barvinok_options *options);
694 virtual void add_term(const vector<int>& powers, evalue *f2) = 0;
695 virtual ~cumulator() {}
698 void cumulator::cumulate(barvinok_options *options)
700 evalue cum; // factor * 1 * E_num[0]/1 * (E_num[0]-1)/2 *...
701 evalue f;
702 evalue t; // E_num[0] - (m-1)
703 evalue *cst;
704 evalue mone;
706 if (options->lookup_table) {
707 value_init(mone.d);
708 evalue_set_si(&mone, -1, 1);
711 value_init(cum.d);
712 evalue_copy(&cum, factor);
713 value_init(f.d);
714 value_init(f.x.n);
715 value_set_si(f.d, 1);
716 value_set_si(f.x.n, 1);
717 value_init(t.d);
718 evalue_copy(&t, v);
720 if (!options->lookup_table) {
721 for (cst = &t; value_zero_p(cst->d); ) {
722 if (cst->x.p->type == fractional)
723 cst = &cst->x.p->arr[1];
724 else
725 cst = &cst->x.p->arr[0];
729 for (int m = 0; m < r->len; ++m) {
730 if (m > 0) {
731 if (m > 1) {
732 value_set_si(f.d, m);
733 emul(&f, &cum);
734 if (!options->lookup_table)
735 value_subtract(cst->x.n, cst->x.n, cst->d);
736 else
737 eadd(&mone, &t);
739 emul(&t, &cum);
741 dpoly_r_term_list& current = r->c[r->len-1-m];
742 dpoly_r_term_list::iterator j;
743 for (j = current.begin(); j != current.end(); ++j) {
744 if ((*j)->coeff == 0)
745 continue;
746 evalue *f2 = new evalue;
747 value_init(f2->d);
748 value_init(f2->x.n);
749 zz2value((*j)->coeff, f2->x.n);
750 zz2value(r->denom, f2->d);
751 emul(&cum, f2);
753 add_term((*j)->powers, f2);
756 free_evalue_refs(&f);
757 free_evalue_refs(&t);
758 free_evalue_refs(&cum);
759 if (options->lookup_table)
760 free_evalue_refs(&mone);
763 struct E_poly_term {
764 vector<int> powers;
765 evalue *E;
768 struct ie_cum : public cumulator {
769 vector<E_poly_term *> terms;
771 ie_cum(evalue *factor, evalue *v, dpoly_r *r) : cumulator(factor, v, r) {}
773 virtual void add_term(const vector<int>& powers, evalue *f2);
776 void ie_cum::add_term(const vector<int>& powers, evalue *f2)
778 int k;
779 for (k = 0; k < terms.size(); ++k) {
780 if (terms[k]->powers == powers) {
781 eadd(f2, terms[k]->E);
782 free_evalue_refs(f2);
783 delete f2;
784 break;
787 if (k >= terms.size()) {
788 E_poly_term *ET = new E_poly_term;
789 ET->powers = powers;
790 ET->E = f2;
791 terms.push_back(ET);
795 struct ienumerator : public signed_cone_consumer, public vertex_decomposer,
796 public ienumerator_base {
797 //Polyhedron *pVD;
798 mat_ZZ den;
799 mat_ZZ vertex;
800 mpq_t tcount;
801 Value tz;
803 ienumerator(Polyhedron *P, unsigned dim, Param_Polyhedron *PP) :
804 vertex_decomposer(PP, *this), ienumerator_base(dim, this) {
805 vertex.SetDims(1, dim);
807 den.SetDims(dim, dim);
808 mpq_init(tcount);
809 value_init(tz);
812 ~ienumerator() {
813 mpq_clear(tcount);
814 value_clear(tz);
817 virtual void handle(const signed_cone& sc, barvinok_options *options);
818 void reduce(evalue *factor, const mat_ZZ& num, const mat_ZZ& den_f,
819 barvinok_options *options);
822 void ienumerator::reduce(evalue *factor, const mat_ZZ& num, const mat_ZZ& den_f,
823 barvinok_options *options)
825 unsigned len = den_f.NumRows(); // number of factors in den
826 unsigned dim = num.NumCols();
827 assert(num.NumRows() == 1);
829 if (dim == 0) {
830 eadd(factor, vE[vert]);
831 return;
834 vec_ZZ den_s;
835 mat_ZZ den_r;
836 vec_ZZ num_s;
837 mat_ZZ num_p;
839 split_one(num, num_s, num_p, den_f, den_s, den_r);
841 vec_ZZ den_p;
842 den_p.SetLength(len);
844 ZZ one;
845 one = 1;
846 normalize(one, num_s, num_p, den_s, den_p, den_r);
847 if (one != 1)
848 emul(&mone, factor);
850 int only_param = 0;
851 int no_param = 0;
852 for (int k = 0; k < len; ++k) {
853 if (den_p[k] == 0)
854 ++no_param;
855 else if (den_s[k] == 0)
856 ++only_param;
858 if (no_param == 0) {
859 reduce(factor, num_p, den_r, options);
860 } else {
861 int k, l;
862 mat_ZZ pden;
863 pden.SetDims(only_param, dim-1);
865 for (k = 0, l = 0; k < len; ++k)
866 if (den_s[k] == 0)
867 pden[l++] = den_r[k];
869 for (k = 0; k < len; ++k)
870 if (den_p[k] == 0)
871 break;
873 zz2value(num_s[0], tz);
874 dpoly n(no_param, tz);
875 zz2value(den_s[k], tz);
876 dpoly D(no_param, tz, 1);
877 for ( ; ++k < len; )
878 if (den_p[k] == 0) {
879 zz2value(den_s[k], tz);
880 dpoly fact(no_param, tz, 1);
881 D *= fact;
884 dpoly_r * r = 0;
885 // if no_param + only_param == len then all powers
886 // below will be all zero
887 if (no_param + only_param == len) {
888 if (E_num(0, dim) != 0)
889 r = new dpoly_r(n, len);
890 else {
891 mpq_set_si(tcount, 0, 1);
892 one = 1;
893 n.div(D, tcount, 1);
895 if (value_notzero_p(mpq_numref(tcount))) {
896 evalue f;
897 value_init(f.d);
898 value_init(f.x.n);
899 value_assign(f.x.n, mpq_numref(tcount));
900 value_assign(f.d, mpq_denref(tcount));
901 emul(&f, factor);
902 reduce(factor, num_p, pden, options);
903 free_evalue_refs(&f);
905 return;
907 } else {
908 for (k = 0; k < len; ++k) {
909 if (den_s[k] == 0 || den_p[k] == 0)
910 continue;
912 zz2value(den_s[k], tz);
913 dpoly pd(no_param-1, tz, 1);
915 int l;
916 for (l = 0; l < k; ++l)
917 if (den_r[l] == den_r[k])
918 break;
920 if (r == 0)
921 r = new dpoly_r(n, pd, l, len);
922 else {
923 dpoly_r *nr = new dpoly_r(r, pd, l, len);
924 delete r;
925 r = nr;
929 dpoly_r *rc = r->div(D);
930 delete r;
931 r = rc;
932 if (E_num(0, dim) == 0) {
933 int common = pden.NumRows();
934 dpoly_r_term_list& final = r->c[r->len-1];
935 int rows;
936 evalue t;
937 evalue f;
938 value_init(f.d);
939 value_init(f.x.n);
940 zz2value(r->denom, f.d);
941 dpoly_r_term_list::iterator j;
942 for (j = final.begin(); j != final.end(); ++j) {
943 if ((*j)->coeff == 0)
944 continue;
945 rows = common;
946 for (int k = 0; k < r->dim; ++k) {
947 int n = (*j)->powers[k];
948 if (n == 0)
949 continue;
950 pden.SetDims(rows+n, pden.NumCols());
951 for (int l = 0; l < n; ++l)
952 pden[rows+l] = den_r[k];
953 rows += n;
955 value_init(t.d);
956 evalue_copy(&t, factor);
957 zz2value((*j)->coeff, f.x.n);
958 emul(&f, &t);
959 reduce(&t, num_p, pden, options);
960 free_evalue_refs(&t);
962 free_evalue_refs(&f);
963 } else {
964 ie_cum cum(factor, E_num(0, dim), r);
965 cum.cumulate(options);
967 int common = pden.NumRows();
968 int rows;
969 for (int j = 0; j < cum.terms.size(); ++j) {
970 rows = common;
971 pden.SetDims(rows, pden.NumCols());
972 for (int k = 0; k < r->dim; ++k) {
973 int n = cum.terms[j]->powers[k];
974 if (n == 0)
975 continue;
976 pden.SetDims(rows+n, pden.NumCols());
977 for (int l = 0; l < n; ++l)
978 pden[rows+l] = den_r[k];
979 rows += n;
981 reduce(cum.terms[j]->E, num_p, pden, options);
982 free_evalue_refs(cum.terms[j]->E);
983 delete cum.terms[j]->E;
984 delete cum.terms[j];
987 delete r;
991 static int type_offset(enode *p)
993 return p->type == fractional ? 1 :
994 p->type == flooring ? 1 : 0;
997 static int edegree(evalue *e)
999 int d = 0;
1000 enode *p;
1002 if (value_notzero_p(e->d))
1003 return 0;
1005 p = e->x.p;
1006 int i = type_offset(p);
1007 if (p->size-i-1 > d)
1008 d = p->size - i - 1;
1009 for (; i < p->size; i++) {
1010 int d2 = edegree(&p->arr[i]);
1011 if (d2 > d)
1012 d = d2;
1014 return d;
1017 void ienumerator::handle(const signed_cone& sc, barvinok_options *options)
1019 assert(sc.det == 1);
1020 assert(sc.rays.NumRows() == dim);
1022 lattice_point(V, sc.rays, vertex[0], E_vertex, options);
1024 den = sc.rays;
1026 evalue one;
1027 value_init(one.d);
1028 evalue_set_si(&one, sc.sign, 1);
1029 reduce(&one, vertex, den, options);
1030 free_evalue_refs(&one);
1032 for (int i = 0; i < dim; ++i)
1033 if (E_vertex[i])
1034 evalue_free(E_vertex[i]);
1037 struct bfenumerator : public vertex_decomposer, public bf_base,
1038 public ienumerator_base {
1039 evalue *factor;
1041 bfenumerator(Polyhedron *P, unsigned dim, Param_Polyhedron *PP) :
1042 vertex_decomposer(PP, *this),
1043 bf_base(dim), ienumerator_base(dim, this) {
1044 lower = 0;
1045 factor = NULL;
1048 ~bfenumerator() {
1051 virtual void handle(const signed_cone& sc, barvinok_options *options);
1052 virtual void base(mat_ZZ& factors, bfc_vec& v);
1054 bfc_term_base* new_bf_term(int len) {
1055 bfe_term* t = new bfe_term(len);
1056 return t;
1059 virtual void set_factor(bfc_term_base *t, int k, int change) {
1060 bfe_term* bfet = static_cast<bfe_term *>(t);
1061 factor = bfet->factors[k];
1062 assert(factor != NULL);
1063 bfet->factors[k] = NULL;
1064 if (change)
1065 emul(&mone, factor);
1068 virtual void set_factor(bfc_term_base *t, int k, mpq_t &q, int change) {
1069 bfe_term* bfet = static_cast<bfe_term *>(t);
1070 factor = bfet->factors[k];
1071 assert(factor != NULL);
1072 bfet->factors[k] = NULL;
1074 evalue f;
1075 value_init(f.d);
1076 value_init(f.x.n);
1077 if (change)
1078 value_oppose(f.x.n, mpq_numref(q));
1079 else
1080 value_assign(f.x.n, mpq_numref(q));
1081 value_assign(f.d, mpq_denref(q));
1082 emul(&f, factor);
1083 free_evalue_refs(&f);
1086 virtual void set_factor(bfc_term_base *t, int k, const QQ& c, int change) {
1087 bfe_term* bfet = static_cast<bfe_term *>(t);
1089 factor = new evalue;
1091 evalue f;
1092 value_init(f.d);
1093 value_init(f.x.n);
1094 zz2value(c.n, f.x.n);
1095 if (change)
1096 value_oppose(f.x.n, f.x.n);
1097 zz2value(c.d, f.d);
1099 value_init(factor->d);
1100 evalue_copy(factor, bfet->factors[k]);
1101 emul(&f, factor);
1102 free_evalue_refs(&f);
1105 void set_factor(evalue *f, int change) {
1106 if (change)
1107 emul(&mone, f);
1108 factor = f;
1111 virtual void insert_term(bfc_term_base *t, int i) {
1112 bfe_term* bfet = static_cast<bfe_term *>(t);
1113 int len = t->terms.NumRows()-1; // already increased by one
1115 bfet->factors.resize(len+1);
1116 for (int j = len; j > i; --j) {
1117 bfet->factors[j] = bfet->factors[j-1];
1118 t->terms[j] = t->terms[j-1];
1120 bfet->factors[i] = factor;
1121 factor = NULL;
1124 virtual void update_term(bfc_term_base *t, int i) {
1125 bfe_term* bfet = static_cast<bfe_term *>(t);
1127 eadd(factor, bfet->factors[i]);
1128 free_evalue_refs(factor);
1129 delete factor;
1132 virtual bool constant_vertex(int dim) { return E_num(0, dim) == 0; }
1134 virtual void cum(bf_reducer *bfr, bfc_term_base *t, int k, dpoly_r *r,
1135 barvinok_options *options);
1138 enumerator_base *enumerator_base::create(Polyhedron *P, unsigned dim,
1139 Param_Polyhedron *PP,
1140 barvinok_options *options)
1142 enumerator_base *eb;
1144 if (options->incremental_specialization == BV_SPECIALIZATION_BF)
1145 eb = new bfenumerator(P, dim, PP);
1146 else if (options->incremental_specialization == BV_SPECIALIZATION_DF)
1147 eb = new ienumerator(P, dim, PP);
1148 else
1149 eb = new enumerator(P, dim, PP);
1151 return eb;
1154 struct bfe_cum : public cumulator {
1155 bfenumerator *bfe;
1156 bfc_term_base *told;
1157 int k;
1158 bf_reducer *bfr;
1160 bfe_cum(evalue *factor, evalue *v, dpoly_r *r, bf_reducer *bfr,
1161 bfc_term_base *t, int k, bfenumerator *e) :
1162 cumulator(factor, v, r), told(t), k(k),
1163 bfr(bfr), bfe(e) {
1166 virtual void add_term(const vector<int>& powers, evalue *f2);
1169 void bfe_cum::add_term(const vector<int>& powers, evalue *f2)
1171 bfr->update_powers(powers);
1173 bfc_term_base * t = bfe->find_bfc_term(bfr->vn, bfr->npowers, bfr->nnf);
1174 bfe->set_factor(f2, bfr->l_changes % 2);
1175 bfe->add_term(t, told->terms[k], bfr->l_extra_num);
1178 void bfenumerator::cum(bf_reducer *bfr, bfc_term_base *t, int k,
1179 dpoly_r *r, barvinok_options *options)
1181 bfe_term* bfet = static_cast<bfe_term *>(t);
1182 bfe_cum cum(bfet->factors[k], E_num(0, bfr->d), r, bfr, t, k, this);
1183 cum.cumulate(options);
1186 void bfenumerator::base(mat_ZZ& factors, bfc_vec& v)
1188 for (int i = 0; i < v.size(); ++i) {
1189 assert(v[i]->terms.NumRows() == 1);
1190 evalue *factor = static_cast<bfe_term *>(v[i])->factors[0];
1191 eadd(factor, vE[vert]);
1192 delete v[i];
1196 void bfenumerator::handle(const signed_cone& sc, barvinok_options *options)
1198 assert(sc.det == 1);
1199 assert(sc.rays.NumRows() == enumerator_base::dim);
1201 bfe_term* t = new bfe_term(enumerator_base::dim);
1202 vector< bfc_term_base * > v;
1203 v.push_back(t);
1205 t->factors.resize(1);
1207 t->terms.SetDims(1, enumerator_base::dim);
1208 lattice_point(V, sc.rays, t->terms[0], E_vertex, options);
1210 // the elements of factors are always lexpositive
1211 mat_ZZ factors;
1212 int s = setup_factors(sc.rays, factors, t, sc.sign);
1214 t->factors[0] = new evalue;
1215 value_init(t->factors[0]->d);
1216 evalue_set_si(t->factors[0], s, 1);
1217 reduce(factors, v, options);
1219 for (int i = 0; i < enumerator_base::dim; ++i)
1220 if (E_vertex[i])
1221 evalue_free(E_vertex[i]);
1224 static evalue* barvinok_enumerate_ev_f(Polyhedron *P, Polyhedron* C,
1225 barvinok_options *options);
1227 /* Destroys C */
1228 static evalue* barvinok_enumerate_cst(Polyhedron *P, Polyhedron* C,
1229 struct barvinok_options *options)
1231 evalue *eres;
1233 if (emptyQ2(C)) {
1234 Polyhedron_Free(C);
1235 return evalue_zero();
1238 ALLOC(evalue, eres);
1239 value_init(eres->d);
1240 value_set_si(eres->d, 0);
1241 eres->x.p = new_enode(partition, 2, C->Dimension);
1242 EVALUE_SET_DOMAIN(eres->x.p->arr[0],
1243 DomainConstraintSimplify(C, options->MaxRays));
1244 value_set_si(eres->x.p->arr[1].d, 1);
1245 value_init(eres->x.p->arr[1].x.n);
1246 if (emptyQ2(P))
1247 value_set_si(eres->x.p->arr[1].x.n, 0);
1248 else
1249 barvinok_count_with_options(P, &eres->x.p->arr[1].x.n, options);
1250 if (value_mone_p(eres->x.p->arr[1].x.n)) {
1251 value_clear(eres->x.p->arr[1].x.n);
1252 value_set_si(eres->x.p->arr[1].d, -2); /* NaN */
1255 return eres;
1258 static evalue* enumerate(Polyhedron *P, Polyhedron* C,
1259 struct barvinok_options *options)
1261 Polyhedron *next;
1262 Polyhedron *Porig = P;
1263 Polyhedron *Corig = C;
1264 Polyhedron *CEq = NULL, *rVD;
1265 int r = 0;
1266 unsigned nparam = C->Dimension;
1267 evalue *eres;
1268 Matrix *CP = NULL;
1270 evalue factor;
1271 value_init(factor.d);
1272 evalue_set_si(&factor, 1, 1);
1274 /* for now */
1275 POL_ENSURE_FACETS(P);
1276 POL_ENSURE_VERTICES(P);
1277 POL_ENSURE_FACETS(C);
1278 POL_ENSURE_VERTICES(C);
1280 if (C->Dimension == 0 || emptyQ(P) || emptyQ(C)) {
1281 constant:
1282 if (CEq == Porig)
1283 CEq = Polyhedron_Copy(CEq);
1284 eres = barvinok_enumerate_cst(P, CEq ? CEq : Polyhedron_Copy(C), options);
1285 out:
1286 if (CP) {
1287 evalue_backsubstitute(eres, CP, options->MaxRays);
1288 Matrix_Free(CP);
1291 emul(&factor, eres);
1292 if (options->approx->method == BV_APPROX_DROP) {
1293 if (options->approx->approximation == BV_APPROX_SIGN_UPPER)
1294 evalue_frac2polynomial(eres, 1, options->MaxRays);
1295 if (options->approx->approximation == BV_APPROX_SIGN_LOWER)
1296 evalue_frac2polynomial(eres, -1, options->MaxRays);
1297 if (options->approx->approximation == BV_APPROX_SIGN_APPROX)
1298 evalue_frac2polynomial(eres, 0, options->MaxRays);
1300 reduce_evalue(eres);
1301 free_evalue_refs(&factor);
1302 if (P != Porig)
1303 Domain_Free(P);
1304 if (C != Corig)
1305 Polyhedron_Free(C);
1307 return eres;
1309 if (Polyhedron_is_unbounded(P, nparam, options->MaxRays))
1310 goto constant;
1312 if (P->Dimension == nparam) {
1313 CEq = P;
1314 P = Universe_Polyhedron(0);
1315 goto constant;
1317 if (P->NbEq != 0 || C->NbEq != 0) {
1318 Polyhedron *Q = P;
1319 Polyhedron *D = C;
1320 remove_all_equalities(&P, &C, &CP, NULL, nparam, options->MaxRays);
1321 if (C != D && D != Corig)
1322 Polyhedron_Free(D);
1323 if (P != Q && Q != Porig)
1324 Domain_Free(Q);
1325 eres = enumerate(P, C, options);
1326 goto out;
1329 Polyhedron *T = Polyhedron_Factor(P, nparam, NULL, options->MaxRays);
1330 if (T || (P->Dimension == nparam+1)) {
1331 Polyhedron *C2 = C;
1332 Polyhedron *FC = Factor_Context(T ? T : P, nparam, options->MaxRays);
1333 C = DomainIntersection(C, FC, options->MaxRays);
1334 if (C2 != Corig)
1335 Polyhedron_Free(C2);
1336 Polyhedron_Free(FC);
1338 if (T) {
1339 if (P != Porig)
1340 Polyhedron_Free(P);
1341 P = T;
1342 if (T->Dimension == C->Dimension) {
1343 P = T->next;
1344 T->next = NULL;
1345 Polyhedron_Free(T);
1349 next = P->next;
1350 P->next = NULL;
1351 eres = barvinok_enumerate_ev_f(P, C, options);
1352 P->next = next;
1354 if (P->next) {
1355 Polyhedron *Q;
1356 evalue *f;
1358 for (Q = P->next; Q; Q = Q->next) {
1359 Polyhedron *next = Q->next;
1360 Q->next = NULL;
1362 f = barvinok_enumerate_ev_f(Q, C, options);
1363 emul(f, eres);
1364 evalue_free(f);
1366 Q->next = next;
1370 goto out;
1373 evalue* barvinok_enumerate_with_options(Polyhedron *P, Polyhedron* C,
1374 struct barvinok_options *options)
1376 Polyhedron *next, *Cnext, *C1;
1377 Polyhedron *Corig = C;
1378 evalue *eres;
1380 if (P->next)
1381 fprintf(stderr,
1382 "barvinok_enumerate: input is a union; only first polyhedron is enumerated\n");
1384 if (C->next)
1385 fprintf(stderr,
1386 "barvinok_enumerate: context is a union; only first polyhedron is considered\n");
1388 Cnext = C->next;
1389 C->next = NULL;
1390 C1 = Polyhedron_Project(P, C->Dimension);
1391 C = DomainIntersection(C, C1, options->MaxRays);
1392 Polyhedron_Free(C1);
1393 next = P->next;
1394 P->next = NULL;
1396 if (options->approx->method == BV_APPROX_BERNOULLI ||
1397 options->summation == BV_SUM_BERNOULLI) {
1398 int summation = options->summation;
1399 options->summation = BV_SUM_BERNOULLI;
1400 eres = barvinok_summate_unweighted(P, C, options);
1401 options->summation = summation;
1402 } else
1403 eres = enumerate(P, C, options);
1404 Domain_Free(C);
1406 P->next= next;
1407 Corig->next = Cnext;
1409 return eres;
1412 evalue* barvinok_enumerate_ev(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
1414 evalue *E;
1415 barvinok_options *options = barvinok_options_new_with_defaults();
1416 options->MaxRays = MaxRays;
1417 E = barvinok_enumerate_with_options(P, C, options);
1418 barvinok_options_free(options);
1419 return E;
1422 evalue *Param_Polyhedron_Enumerate(Param_Polyhedron *PP, Polyhedron *P,
1423 Polyhedron *C,
1424 struct barvinok_options *options)
1426 evalue *eres;
1427 Param_Domain *D;
1428 unsigned nparam = C->Dimension;
1429 unsigned dim = P->Dimension - nparam;
1431 int nd;
1432 for (nd = 0, D=PP->D; D; ++nd, D=D->next);
1433 evalue_section *s = new evalue_section[nd];
1435 enumerator_base *et = NULL;
1436 try_again:
1437 if (et)
1438 delete et;
1440 et = enumerator_base::create(P, dim, PP, options);
1442 Polyhedron *TC = true_context(P, C, options->MaxRays);
1443 FORALL_REDUCED_DOMAIN(PP, TC, nd, options, i, D, rVD)
1444 Param_Vertices *V;
1446 s[i].E = evalue_zero();
1447 s[i].D = rVD;
1449 FORALL_PVertex_in_ParamPolyhedron(V,D,PP) // _i is internal counter
1450 if (!et->vE[_i])
1451 try {
1452 et->decompose_at(V, _i, options);
1453 } catch (OrthogonalException &e) {
1454 FORALL_REDUCED_DOMAIN_RESET;
1455 for (; i >= 0; --i) {
1456 evalue_free(s[i].E);
1457 Domain_Free(s[i].D);
1459 goto try_again;
1461 eadd(et->vE[_i] , s[i].E);
1462 END_FORALL_PVertex_in_ParamPolyhedron;
1463 evalue_range_reduction_in_domain(s[i].E, rVD);
1464 END_FORALL_REDUCED_DOMAIN
1465 Polyhedron_Free(TC);
1467 delete et;
1468 eres = evalue_from_section_array(s, nd);
1469 delete [] s;
1471 return eres;
1474 static evalue* barvinok_enumerate_ev_f(Polyhedron *P, Polyhedron* C,
1475 barvinok_options *options)
1477 unsigned nparam = C->Dimension;
1478 bool do_scale = options->approx->method == BV_APPROX_SCALE;
1480 if (options->summation == BV_SUM_EULER)
1481 return barvinok_summate_unweighted(P, C, options);
1483 if (options->approx->method == BV_APPROX_VOLUME)
1484 return Param_Polyhedron_Volume(P, C, options);
1486 if (P->Dimension - nparam == 1 && !do_scale)
1487 return ParamLine_Length(P, C, options);
1489 Param_Polyhedron *PP = NULL;
1490 evalue *eres;
1492 if (do_scale) {
1493 eres = scale_bound(P, C, options);
1494 if (eres)
1495 return eres;
1498 PP = Polyhedron2Param_Polyhedron(P, C, options);
1500 if (do_scale)
1501 eres = scale(PP, P, C, options);
1502 else
1503 eres = Param_Polyhedron_Enumerate(PP, P, C, options);
1505 if (PP)
1506 Param_Polyhedron_Free(PP);
1508 return eres;
1511 Enumeration* barvinok_enumerate(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
1513 evalue *EP = barvinok_enumerate_ev(P, C, MaxRays);
1515 return partition2enumeration(EP);
1518 evalue* barvinok_enumerate_union(Polyhedron *D, Polyhedron* C, unsigned MaxRays)
1520 evalue *EP;
1521 gen_fun *gf = barvinok_enumerate_union_series(D, C, MaxRays);
1522 EP = *gf;
1523 delete gf;
1524 return EP;
1527 static int basic_map_card(__isl_take isl_basic_map *bmap, void *user)
1529 isl_ctx *ctx;
1530 isl_pw_qpolynomial **sum = (isl_pw_qpolynomial **)user;
1531 isl_pw_qpolynomial *pwqp;
1532 Polyhedron *P;
1533 evalue *E;
1534 barvinok_options *options;
1535 int options_allocated = 0;
1536 unsigned nparam = isl_basic_map_dim(bmap, isl_dim_param);
1537 unsigned n_in = isl_basic_map_dim(bmap, isl_dim_in);
1538 Polyhedron *U = Universe_Polyhedron(nparam + n_in);
1539 isl_dim *dim, *target_dim;
1541 ctx = isl_basic_map_get_ctx(bmap);
1542 options = isl_ctx_peek_barvinok_options(ctx);
1543 if (!options) {
1544 options = barvinok_options_new_with_defaults();
1545 options_allocated = 1;
1548 target_dim = isl_basic_map_get_dim(bmap);
1549 target_dim = isl_dim_domain(target_dim);
1551 bmap = isl_basic_map_move_dims(bmap, isl_dim_param, nparam,
1552 isl_dim_in, 0, n_in);
1554 dim = isl_basic_map_get_dim(bmap);
1555 dim = isl_dim_domain(dim);
1557 P = isl_basic_map_to_polylib(bmap);
1558 E = barvinok_enumerate_with_options(P, U, options);
1560 pwqp = isl_pw_qpolynomial_from_evalue(dim, E);
1561 pwqp = isl_pw_qpolynomial_move_dims(pwqp, isl_dim_set, 0,
1562 isl_dim_param, nparam, n_in);
1563 pwqp = isl_pw_qpolynomial_reset_dim(pwqp, target_dim);
1564 *sum = isl_pw_qpolynomial_add(*sum, pwqp);
1566 evalue_free(E);
1567 Polyhedron_Free(P);
1568 Polyhedron_Free(U);
1569 if (options_allocated)
1570 barvinok_options_free(options);
1571 isl_basic_map_free(bmap);
1573 return 0;
1576 __isl_give isl_pw_qpolynomial *isl_map_card(__isl_take isl_map *map)
1578 isl_dim *dim;
1579 isl_pw_qpolynomial *sum;
1581 dim = isl_map_get_dim(map);
1582 dim = isl_dim_domain(dim);
1583 sum = isl_pw_qpolynomial_zero(dim);
1585 map = isl_map_make_disjoint(map);
1586 map = isl_map_compute_divs(map);
1588 if (isl_map_foreach_basic_map(map, &basic_map_card, &sum) < 0)
1589 goto error;
1591 isl_map_free(map);
1593 return sum;
1594 error:
1595 isl_map_free(map);
1596 isl_pw_qpolynomial_free(sum);
1597 return NULL;
1600 __isl_give isl_pw_qpolynomial *isl_set_card(__isl_take isl_set *set)
1602 return isl_map_card(isl_map_from_range(set));
1605 static int set_card(__isl_take isl_set *set, void *user)
1607 isl_union_pw_qpolynomial **res = (isl_union_pw_qpolynomial **)user;
1608 isl_pw_qpolynomial *pwqp;
1610 pwqp = isl_set_card(set);
1611 *res = isl_union_pw_qpolynomial_add_pw_qpolynomial(*res, pwqp);
1613 return 0;
1616 __isl_give isl_union_pw_qpolynomial *isl_union_set_card(
1617 __isl_take isl_union_set *uset)
1619 isl_dim *dim;
1620 isl_union_pw_qpolynomial *res;
1622 dim = isl_union_set_get_dim(uset);
1623 res = isl_union_pw_qpolynomial_zero(dim);
1624 if (isl_union_set_foreach_set(uset, &set_card, &res) < 0)
1625 goto error;
1626 isl_union_set_free(uset);
1628 return res;
1629 error:
1630 isl_union_set_free(uset);
1631 isl_union_pw_qpolynomial_free(res);
1632 return NULL;
1635 static int map_card(__isl_take isl_map *map, void *user)
1637 isl_union_pw_qpolynomial **res = (isl_union_pw_qpolynomial **)user;
1638 isl_pw_qpolynomial *pwqp;
1640 pwqp = isl_map_card(map);
1641 *res = isl_union_pw_qpolynomial_add_pw_qpolynomial(*res, pwqp);
1643 return 0;
1646 __isl_give isl_union_pw_qpolynomial *isl_union_map_card(
1647 __isl_take isl_union_map *umap)
1649 isl_dim *dim;
1650 isl_union_pw_qpolynomial *res;
1652 dim = isl_union_map_get_dim(umap);
1653 res = isl_union_pw_qpolynomial_zero(dim);
1654 if (isl_union_map_foreach_map(umap, &map_card, &res) < 0)
1655 goto error;
1656 isl_union_map_free(umap);
1658 return res;
1659 error:
1660 isl_union_map_free(umap);
1661 isl_union_pw_qpolynomial_free(res);
1662 return NULL;