update pet for support for recent clangs
[barvinok.git] / barvinok.cc
blob4c8c4cd63d7b9059a54e7223559e4ab5b07131f3
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.h>
11 #include <isl_set_polylib.h>
12 #include <barvinok/util.h>
13 #include <barvinok/evalue.h>
14 #include "config.h"
15 #include <barvinok/barvinok.h>
16 #include <barvinok/genfun.h>
17 #include <barvinok/options.h>
18 #include <barvinok/sample.h>
19 #include "bfcounter.h"
20 #include "conversion.h"
21 #include "counter.h"
22 #include "decomposer.h"
23 #include "euler.h"
24 #include "lattice_point.h"
25 #include "laurent.h"
26 #include "reduce_domain.h"
27 #include "remove_equalities.h"
28 #include "scale.h"
29 #include "volume.h"
30 #include "bernoulli.h"
31 #include "param_util.h"
32 #include "summate.h"
34 #ifdef NTL_STD_CXX
35 using namespace NTL;
36 #endif
37 using std::cerr;
38 using std::cout;
39 using std::endl;
40 using std::vector;
41 using std::deque;
42 using std::string;
43 using std::ostringstream;
45 #define ALLOC(t,p) p = (t*)malloc(sizeof(*p))
47 class dpoly_n {
48 public:
49 Matrix *coeff;
50 ~dpoly_n() {
51 Matrix_Free(coeff);
53 dpoly_n(int d) {
54 Value d0, one;
55 value_init(d0);
56 value_init(one);
57 value_set_si(one, 1);
58 coeff = Matrix_Alloc(d+1, d+1+1);
59 value_set_si(coeff->p[0][0], 1);
60 value_set_si(coeff->p[0][d+1], 1);
61 for (int i = 1; i <= d; ++i) {
62 value_multiply(coeff->p[i][0], coeff->p[i-1][0], d0);
63 Vector_Combine(coeff->p[i-1], coeff->p[i-1]+1, coeff->p[i]+1,
64 one, d0, i);
65 value_set_si(coeff->p[i][d+1], i);
66 value_multiply(coeff->p[i][d+1], coeff->p[i][d+1], coeff->p[i-1][d+1]);
67 value_decrement(d0, d0);
69 value_clear(d0);
70 value_clear(one);
72 void div(dpoly& d, Vector *count, int sign) {
73 int len = coeff->NbRows;
74 Matrix * c = Matrix_Alloc(coeff->NbRows, coeff->NbColumns);
75 Value tmp;
76 value_init(tmp);
77 for (int i = 0; i < len; ++i) {
78 Vector_Copy(coeff->p[i], c->p[i], len+1);
79 for (int j = 1; j <= i; ++j) {
80 value_multiply(tmp, d.coeff->p[j], c->p[i][len]);
81 value_oppose(tmp, tmp);
82 Vector_Combine(c->p[i], c->p[i-j], c->p[i],
83 c->p[i-j][len], tmp, len);
84 value_multiply(c->p[i][len], c->p[i][len], c->p[i-j][len]);
86 value_multiply(c->p[i][len], c->p[i][len], d.coeff->p[0]);
88 if (sign == -1) {
89 value_set_si(tmp, -1);
90 Vector_Scale(c->p[len-1], count->p, tmp, len);
91 value_assign(count->p[len], c->p[len-1][len]);
92 } else
93 Vector_Copy(c->p[len-1], count->p, len+1);
94 Vector_Normalize(count->p, len+1);
95 value_clear(tmp);
96 Matrix_Free(c);
100 static void add_rays(mat_ZZ& rays, Polyhedron *i, int *r, int nvar = -1,
101 bool all = false)
103 unsigned dim = i->Dimension;
104 if (nvar == -1)
105 nvar = dim;
106 for (int k = 0; k < i->NbRays; ++k) {
107 if (!value_zero_p(i->Ray[k][dim+1]))
108 continue;
109 if (!all && nvar != dim && First_Non_Zero(i->Ray[k]+1, nvar) == -1)
110 continue;
111 values2zz(i->Ray[k]+1, rays[(*r)++], nvar);
115 struct bfe_term : public bfc_term_base {
116 vector<evalue *> factors;
118 bfe_term(int len) : bfc_term_base(len) {
121 ~bfe_term() {
122 for (int i = 0; i < factors.size(); ++i) {
123 if (!factors[i])
124 continue;
125 free_evalue_refs(factors[i]);
126 delete factors[i];
131 static void print_int_vector(int *v, int len, const char *name)
133 cerr << name << endl;
134 for (int j = 0; j < len; ++j) {
135 cerr << v[j] << " ";
137 cerr << endl;
140 static void print_bfc_terms(mat_ZZ& factors, bfc_vec& v)
142 cerr << endl;
143 cerr << "factors" << endl;
144 cerr << factors << endl;
145 for (int i = 0; i < v.size(); ++i) {
146 cerr << "term: " << i << endl;
147 print_int_vector(v[i]->powers, factors.NumRows(), "powers");
148 cerr << "terms" << endl;
149 cerr << v[i]->terms << endl;
150 bfc_term* bfct = static_cast<bfc_term *>(v[i]);
151 cerr << bfct->c << endl;
155 static void print_bfe_terms(mat_ZZ& factors, bfc_vec& v)
157 cerr << endl;
158 cerr << "factors" << endl;
159 cerr << factors << endl;
160 for (int i = 0; i < v.size(); ++i) {
161 cerr << "term: " << i << endl;
162 print_int_vector(v[i]->powers, factors.NumRows(), "powers");
163 cerr << "terms" << endl;
164 cerr << v[i]->terms << endl;
165 bfe_term* bfet = static_cast<bfe_term *>(v[i]);
166 for (int j = 0; j < v[i]->terms.NumRows(); ++j) {
167 const char * test[] = {"a", "b"};
168 print_evalue(stderr, bfet->factors[j], test);
169 fprintf(stderr, "\n");
174 struct bfcounter : public bfcounter_base {
175 mpq_t count;
176 Value tz;
178 bfcounter(unsigned dim) : bfcounter_base(dim) {
179 mpq_init(count);
180 lower = 1;
181 value_init(tz);
183 ~bfcounter() {
184 mpq_clear(count);
185 value_clear(tz);
187 virtual void base(mat_ZZ& factors, bfc_vec& v);
188 virtual void get_count(Value *result) {
189 assert(value_one_p(&count[0]._mp_den));
190 value_assign(*result, &count[0]._mp_num);
194 void bfcounter::base(mat_ZZ& factors, bfc_vec& v)
196 unsigned nf = factors.NumRows();
198 for (int i = 0; i < v.size(); ++i) {
199 bfc_term* bfct = static_cast<bfc_term *>(v[i]);
200 int total_power = 0;
201 // factor is always positive, so we always
202 // change signs
203 for (int k = 0; k < nf; ++k)
204 total_power += v[i]->powers[k];
206 int j;
207 for (j = 0; j < nf; ++j)
208 if (v[i]->powers[j] > 0)
209 break;
211 zz2value(factors[j][0], tz);
212 dpoly D(total_power, tz, 1);
213 for (int k = 1; k < v[i]->powers[j]; ++k) {
214 zz2value(factors[j][0], tz);
215 dpoly fact(total_power, tz, 1);
216 D *= fact;
218 for ( ; ++j < nf; )
219 for (int k = 0; k < v[i]->powers[j]; ++k) {
220 zz2value(factors[j][0], tz);
221 dpoly fact(total_power, tz, 1);
222 D *= fact;
225 for (int k = 0; k < v[i]->terms.NumRows(); ++k) {
226 zz2value(v[i]->terms[k][0], tz);
227 dpoly n(total_power, tz);
228 mpq_set_si(tcount, 0, 1);
229 n.div(D, tcount, 1);
230 if (total_power % 2)
231 bfct->c[k].n = -bfct->c[k].n;
232 zz2value(bfct->c[k].n, tn);
233 zz2value(bfct->c[k].d, td);
235 mpz_mul(mpq_numref(tcount), mpq_numref(tcount), tn);
236 mpz_mul(mpq_denref(tcount), mpq_denref(tcount), td);
237 mpq_canonicalize(tcount);
238 mpq_add(count, count, tcount);
240 delete v[i];
245 /* Check whether the polyhedron is unbounded and if so,
246 * check whether it has any (and therefore an infinite number of)
247 * integer points.
248 * If one of the vertices is integer, then we are done.
249 * Otherwise, transform the polyhedron such that one of the rays
250 * is the first unit vector and cut it off at a height that ensures
251 * that if the whole polyhedron has any points, then the remaining part
252 * has integer points. In particular we add the largest coefficient
253 * of a ray to the highest vertex (rounded up).
255 static bool Polyhedron_is_infinite(Polyhedron *P, Value* result,
256 barvinok_options *options)
258 int r = 0;
259 Matrix *M, *M2;
260 Value c, tmp;
261 Value g;
262 bool first;
263 Vector *v;
264 Value offset, size;
265 Polyhedron *R;
267 if (P->NbBid == 0)
268 for (; r < P->NbRays; ++r)
269 if (value_zero_p(P->Ray[r][P->Dimension+1]))
270 break;
271 if (P->NbBid == 0 && r == P->NbRays)
272 return false;
274 if (options->count_sample_infinite) {
275 Vector *sample;
277 sample = Polyhedron_Sample(P, options);
278 if (!sample)
279 value_set_si(*result, 0);
280 else {
281 value_set_si(*result, -1);
282 Vector_Free(sample);
284 return true;
287 for (int i = 0; i < P->NbRays; ++i)
288 if (value_one_p(P->Ray[i][1+P->Dimension])) {
289 value_set_si(*result, -1);
290 return true;
293 value_init(g);
294 M = Matrix_Alloc(P->Dimension+1, P->Dimension+1);
295 Vector_Gcd(P->Ray[r]+1, P->Dimension, &g);
296 Vector_AntiScale(P->Ray[r]+1, M->p[0], g, P->Dimension+1);
297 int ok = unimodular_complete(M, 1);
298 assert(ok);
299 value_set_si(M->p[P->Dimension][P->Dimension], 1);
300 M2 = Transpose(M);
301 Matrix_Free(M);
302 P = Polyhedron_Preimage(P, M2, 0);
303 Matrix_Free(M2);
304 value_clear(g);
306 first = true;
307 value_init(offset);
308 value_init(size);
309 value_init(tmp);
310 value_set_si(size, 0);
312 for (int i = 0; i < P->NbBid; ++i) {
313 value_absolute(tmp, P->Ray[i][1]);
314 if (value_gt(tmp, size))
315 value_assign(size, tmp);
317 for (int i = P->NbBid; i < P->NbRays; ++i) {
318 if (value_zero_p(P->Ray[i][P->Dimension+1])) {
319 if (value_gt(P->Ray[i][1], size))
320 value_assign(size, P->Ray[i][1]);
321 continue;
323 mpz_cdiv_q(tmp, P->Ray[i][1], P->Ray[i][P->Dimension+1]);
324 if (first || value_gt(tmp, offset)) {
325 value_assign(offset, tmp);
326 first = false;
329 value_addto(offset, offset, size);
330 value_clear(size);
331 value_clear(tmp);
333 v = Vector_Alloc(P->Dimension+2);
334 value_set_si(v->p[0], 1);
335 value_set_si(v->p[1], -1);
336 value_assign(v->p[1+P->Dimension], offset);
337 R = AddConstraints(v->p, 1, P, options->MaxRays);
338 Polyhedron_Free(P);
339 P = R;
341 value_clear(offset);
342 Vector_Free(v);
344 value_init(c);
345 barvinok_count_with_options(P, &c, options);
346 Polyhedron_Free(P);
347 if (value_zero_p(c))
348 value_set_si(*result, 0);
349 else
350 value_set_si(*result, -1);
351 value_clear(c);
353 return true;
356 static void evalue2value(evalue *e, Value *v)
358 if (EVALUE_IS_ZERO(*e)) {
359 value_set_si(*v, 0);
360 return;
363 if (value_notzero_p(e->d)) {
364 assert(value_one_p(e->d));
365 value_assign(*v, e->x.n);
366 return;
369 assert(e->x.p->type == partition);
370 assert(e->x.p->size == 2);
371 assert(EVALUE_DOMAIN(e->x.p->arr[0])->Dimension == 0);
372 evalue2value(&e->x.p->arr[1], v);
375 static void barvinok_count_f(Polyhedron *P, Value* result,
376 barvinok_options *options);
378 void barvinok_count_with_options(Polyhedron *P, Value* result,
379 struct barvinok_options *options)
381 unsigned dim;
382 int allocated = 0;
383 Polyhedron *Q;
384 bool infinite = false;
386 if (P->next)
387 fprintf(stderr,
388 "barvinok_count: input is a union; only first polyhedron is counted\n");
390 if (emptyQ2(P)) {
391 value_set_si(*result, 0);
392 return;
394 if (P->NbEq != 0) {
395 Q = NULL;
396 do {
397 P = remove_equalities(P, options->MaxRays);
398 if (P)
399 P = DomainConstraintSimplify(P, options->MaxRays);
400 if (Q)
401 Polyhedron_Free(Q);
402 Q = P;
403 } while (P && !emptyQ(P) && P->NbEq != 0);
404 if (!P || emptyQ(P)) {
405 Polyhedron_Free(P);
406 value_set_si(*result, 0);
407 return;
409 allocated = 1;
411 if (Polyhedron_is_infinite(P, result, options)) {
412 if (allocated)
413 Polyhedron_Free(P);
414 return;
416 if (P->Dimension == 0) {
417 /* Test whether the constraints are satisfied */
418 POL_ENSURE_VERTICES(P);
419 value_set_si(*result, !emptyQ(P));
420 if (allocated)
421 Polyhedron_Free(P);
422 return;
424 if (options->summation == BV_SUM_BERNOULLI) {
425 Polyhedron *C = Universe_Polyhedron(0);
426 evalue *sum = barvinok_summate_unweighted(P, C, options);
427 Polyhedron_Free(C);
428 evalue2value(sum, result);
429 evalue_free(sum);
430 return;
432 Q = Polyhedron_Factor(P, 0, NULL, options->MaxRays);
433 if (Q) {
434 if (allocated)
435 Polyhedron_Free(P);
436 P = Q;
437 allocated = 1;
440 barvinok_count_f(P, result, options);
441 if (value_neg_p(*result))
442 infinite = true;
443 if (Q && P->next && value_notzero_p(*result)) {
444 Value factor;
445 value_init(factor);
447 for (Q = P->next; Q; Q = Q->next) {
448 barvinok_count_f(Q, &factor, options);
449 if (value_neg_p(factor)) {
450 infinite = true;
451 continue;
452 } else if (Q->next && value_zero_p(factor)) {
453 value_set_si(*result, 0);
454 break;
456 value_multiply(*result, *result, factor);
459 value_clear(factor);
462 if (allocated)
463 Domain_Free(P);
464 if (infinite)
465 value_set_si(*result, -1);
468 void barvinok_count(Polyhedron *P, Value* result, unsigned NbMaxCons)
470 barvinok_options *options = barvinok_options_new_with_defaults();
471 options->MaxRays = NbMaxCons;
472 barvinok_count_with_options(P, result, options);
473 barvinok_options_free(options);
476 static void barvinok_count_f(Polyhedron *P, Value* result,
477 barvinok_options *options)
479 if (emptyQ2(P)) {
480 value_set_si(*result, 0);
481 return;
484 if (P->Dimension == 1)
485 return Line_Length(P, result);
487 int c = P->NbConstraints;
488 POL_ENSURE_FACETS(P);
489 if (c != P->NbConstraints || P->NbEq != 0) {
490 Polyhedron *next = P->next;
491 P->next = NULL;
492 barvinok_count_with_options(P, result, options);
493 P->next = next;
494 return;
497 POL_ENSURE_VERTICES(P);
499 if (Polyhedron_is_infinite(P, result, options))
500 return;
502 np_base *cnt;
503 if (options->incremental_specialization == BV_SPECIALIZATION_BF)
504 cnt = new bfcounter(P->Dimension);
505 else if (options->incremental_specialization == BV_SPECIALIZATION_DF)
506 cnt = new icounter(P->Dimension);
507 else if (options->incremental_specialization == BV_SPECIALIZATION_TODD)
508 cnt = new tcounter(P->Dimension, options->max_index);
509 else
510 cnt = new counter(P->Dimension, options->max_index);
511 cnt->start(P, options);
513 cnt->get_count(result);
514 delete cnt;
517 static void uni_polynom(int param, Vector *c, evalue *EP)
519 unsigned dim = c->Size-2;
520 value_init(EP->d);
521 value_set_si(EP->d,0);
522 EP->x.p = new_enode(polynomial, dim+1, param+1);
523 for (int j = 0; j <= dim; ++j)
524 evalue_set(&EP->x.p->arr[j], c->p[j], c->p[dim+1]);
527 typedef evalue * evalue_p;
529 struct enumerator_base {
530 unsigned dim;
531 evalue ** vE;
532 evalue mone;
533 vertex_decomposer *vpd;
535 enumerator_base(unsigned dim, vertex_decomposer *vpd)
537 this->dim = dim;
538 this->vpd = vpd;
540 vE = new evalue_p[vpd->PP->nbV];
541 for (int j = 0; j < vpd->PP->nbV; ++j)
542 vE[j] = 0;
544 value_init(mone.d);
545 evalue_set_si(&mone, -1, 1);
548 void decompose_at(Param_Vertices *V, int _i, barvinok_options *options) {
549 //this->pVD = pVD;
551 vE[_i] = new evalue;
552 value_init(vE[_i]->d);
553 evalue_set_si(vE[_i], 0, 1);
555 vpd->decompose_at_vertex(V, _i, options);
558 virtual ~enumerator_base() {
559 for (int j = 0; j < vpd->PP->nbV; ++j)
560 if (vE[j]) {
561 free_evalue_refs(vE[j]);
562 delete vE[j];
564 delete [] vE;
566 free_evalue_refs(&mone);
569 static enumerator_base *create(Polyhedron *P, unsigned dim,
570 Param_Polyhedron *PP,
571 barvinok_options *options);
574 struct enumerator : public signed_cone_consumer, public vertex_decomposer,
575 public enumerator_base {
576 vec_ZZ lambda;
577 vec_ZZ den;
578 term_info num;
579 Vector *c;
580 mpq_t count;
581 Value tz;
583 enumerator(Polyhedron *P, unsigned dim, Param_Polyhedron *PP) :
584 vertex_decomposer(PP, *this), enumerator_base(dim, this) {
585 randomvector(P, lambda, dim, 0);
586 den.SetLength(dim);
587 c = Vector_Alloc(dim+2);
589 mpq_init(count);
590 value_init(tz);
593 ~enumerator() {
594 mpq_clear(count);
595 Vector_Free(c);
596 value_clear(tz);
599 virtual void handle(const signed_cone& sc, barvinok_options *options);
602 void enumerator::handle(const signed_cone& sc, barvinok_options *options)
604 int sign = sc.sign;
605 int r = 0;
606 assert(sc.rays.NumRows() == dim);
607 for (int k = 0; k < dim; ++k) {
608 if (lambda * sc.rays[k] == 0)
609 throw Orthogonal;
612 lattice_point(V, sc.rays, lambda, &num, sc.det, options);
613 den = sc.rays * lambda;
615 if (dim % 2)
616 sign = -sign;
618 zz2value(den[0], tz);
619 dpoly n(dim, tz, 1);
620 for (int k = 1; k < dim; ++k) {
621 zz2value(den[k], tz);
622 dpoly fact(dim, tz, 1);
623 n *= fact;
625 if (num.E != NULL) {
626 dpoly_n d(dim);
627 d.div(n, c, sign);
628 for (unsigned long i = 0; i < sc.det; ++i) {
629 evalue *EV = evalue_polynomial(c, num.E[i]);
630 eadd(EV, vE[vert]);
631 evalue_free(EV);
632 evalue_free(num.E[i]);
634 delete [] num.E;
635 } else {
636 mpq_set_si(count, 0, 1);
637 if (num.constant.length() == 1) {
638 zz2value(num.constant[0], tz);
639 dpoly d(dim, tz);
640 d.div(n, count, sign);
641 } else {
642 dpoly_n d(dim);
643 d.div(n, c, sign);
644 Value x, sum, acc;
645 value_init(x);
646 value_init(acc);
647 for (unsigned long i = 0; i < sc.det; ++i) {
648 value_assign(acc, c->p[dim]);
649 zz2value(num.constant[i], x);
650 for (int j = dim-1; j >= 0; --j) {
651 value_multiply(acc, acc, x);
652 value_addto(acc, acc, c->p[j]);
654 value_addto(mpq_numref(count), mpq_numref(count), acc);
656 mpz_set(mpq_denref(count), c->p[dim+1]);
657 value_clear(acc);
658 value_clear(x);
660 evalue EV;
661 value_init(EV.d);
662 evalue_set(&EV, &count[0]._mp_num, &count[0]._mp_den);
663 eadd(&EV, vE[vert]);
664 free_evalue_refs(&EV);
668 struct ienumerator_base : enumerator_base {
669 evalue ** E_vertex;
671 ienumerator_base(unsigned dim, vertex_decomposer *vpd) :
672 enumerator_base(dim,vpd) {
673 E_vertex = new evalue_p[dim];
676 virtual ~ienumerator_base() {
677 delete [] E_vertex;
680 evalue *E_num(int i, int d) {
681 return E_vertex[i + (dim-d)];
685 struct cumulator {
686 evalue *factor;
687 evalue *v;
688 dpoly_r *r;
690 cumulator(evalue *factor, evalue *v, dpoly_r *r) :
691 factor(factor), v(v), r(r) {}
693 void cumulate(barvinok_options *options);
695 virtual void add_term(const vector<int>& powers, evalue *f2) = 0;
696 virtual ~cumulator() {}
699 void cumulator::cumulate(barvinok_options *options)
701 evalue cum; // factor * 1 * E_num[0]/1 * (E_num[0]-1)/2 *...
702 evalue f;
703 evalue t; // E_num[0] - (m-1)
704 evalue *cst;
705 evalue mone;
707 if (options->lookup_table) {
708 value_init(mone.d);
709 evalue_set_si(&mone, -1, 1);
712 value_init(cum.d);
713 evalue_copy(&cum, factor);
714 value_init(f.d);
715 value_init(f.x.n);
716 value_set_si(f.d, 1);
717 value_set_si(f.x.n, 1);
718 value_init(t.d);
719 evalue_copy(&t, v);
721 if (!options->lookup_table) {
722 for (cst = &t; value_zero_p(cst->d); ) {
723 if (cst->x.p->type == fractional)
724 cst = &cst->x.p->arr[1];
725 else
726 cst = &cst->x.p->arr[0];
730 for (int m = 0; m < r->len; ++m) {
731 if (m > 0) {
732 if (m > 1) {
733 value_set_si(f.d, m);
734 emul(&f, &cum);
735 if (!options->lookup_table)
736 value_subtract(cst->x.n, cst->x.n, cst->d);
737 else
738 eadd(&mone, &t);
740 emul(&t, &cum);
742 dpoly_r_term_list& current = r->c[r->len-1-m];
743 dpoly_r_term_list::iterator j;
744 for (j = current.begin(); j != current.end(); ++j) {
745 if ((*j)->coeff == 0)
746 continue;
747 evalue *f2 = new evalue;
748 value_init(f2->d);
749 value_init(f2->x.n);
750 zz2value((*j)->coeff, f2->x.n);
751 zz2value(r->denom, f2->d);
752 emul(&cum, f2);
754 add_term((*j)->powers, f2);
757 free_evalue_refs(&f);
758 free_evalue_refs(&t);
759 free_evalue_refs(&cum);
760 if (options->lookup_table)
761 free_evalue_refs(&mone);
764 struct E_poly_term {
765 vector<int> powers;
766 evalue *E;
769 struct ie_cum : public cumulator {
770 vector<E_poly_term *> terms;
772 ie_cum(evalue *factor, evalue *v, dpoly_r *r) : cumulator(factor, v, r) {}
774 virtual void add_term(const vector<int>& powers, evalue *f2);
777 void ie_cum::add_term(const vector<int>& powers, evalue *f2)
779 int k;
780 for (k = 0; k < terms.size(); ++k) {
781 if (terms[k]->powers == powers) {
782 eadd(f2, terms[k]->E);
783 free_evalue_refs(f2);
784 delete f2;
785 break;
788 if (k >= terms.size()) {
789 E_poly_term *ET = new E_poly_term;
790 ET->powers = powers;
791 ET->E = f2;
792 terms.push_back(ET);
796 struct ienumerator : public signed_cone_consumer, public vertex_decomposer,
797 public ienumerator_base {
798 //Polyhedron *pVD;
799 mat_ZZ den;
800 mat_ZZ vertex;
801 mpq_t tcount;
802 Value tz;
804 ienumerator(Polyhedron *P, unsigned dim, Param_Polyhedron *PP) :
805 vertex_decomposer(PP, *this), ienumerator_base(dim, this) {
806 vertex.SetDims(1, dim);
808 den.SetDims(dim, dim);
809 mpq_init(tcount);
810 value_init(tz);
813 ~ienumerator() {
814 mpq_clear(tcount);
815 value_clear(tz);
818 virtual void handle(const signed_cone& sc, barvinok_options *options);
819 void reduce(evalue *factor, const mat_ZZ& num, const mat_ZZ& den_f,
820 barvinok_options *options);
823 void ienumerator::reduce(evalue *factor, const mat_ZZ& num, const mat_ZZ& den_f,
824 barvinok_options *options)
826 unsigned len = den_f.NumRows(); // number of factors in den
827 unsigned dim = num.NumCols();
828 assert(num.NumRows() == 1);
830 if (dim == 0) {
831 eadd(factor, vE[vert]);
832 return;
835 vec_ZZ den_s;
836 mat_ZZ den_r;
837 vec_ZZ num_s;
838 mat_ZZ num_p;
840 split_one(num, num_s, num_p, den_f, den_s, den_r);
842 vec_ZZ den_p;
843 den_p.SetLength(len);
845 ZZ one;
846 one = 1;
847 normalize(one, num_s, num_p, den_s, den_p, den_r);
848 if (one != 1)
849 emul(&mone, factor);
851 int only_param = 0;
852 int no_param = 0;
853 for (int k = 0; k < len; ++k) {
854 if (den_p[k] == 0)
855 ++no_param;
856 else if (den_s[k] == 0)
857 ++only_param;
859 if (no_param == 0) {
860 reduce(factor, num_p, den_r, options);
861 } else {
862 int k, l;
863 mat_ZZ pden;
864 pden.SetDims(only_param, dim-1);
866 for (k = 0, l = 0; k < len; ++k)
867 if (den_s[k] == 0)
868 pden[l++] = den_r[k];
870 for (k = 0; k < len; ++k)
871 if (den_p[k] == 0)
872 break;
874 zz2value(num_s[0], tz);
875 dpoly n(no_param, tz);
876 zz2value(den_s[k], tz);
877 dpoly D(no_param, tz, 1);
878 for ( ; ++k < len; )
879 if (den_p[k] == 0) {
880 zz2value(den_s[k], tz);
881 dpoly fact(no_param, tz, 1);
882 D *= fact;
885 dpoly_r * r = 0;
886 // if no_param + only_param == len then all powers
887 // below will be all zero
888 if (no_param + only_param == len) {
889 if (E_num(0, dim) != 0)
890 r = new dpoly_r(n, len);
891 else {
892 mpq_set_si(tcount, 0, 1);
893 one = 1;
894 n.div(D, tcount, 1);
896 if (value_notzero_p(mpq_numref(tcount))) {
897 evalue f;
898 value_init(f.d);
899 value_init(f.x.n);
900 value_assign(f.x.n, mpq_numref(tcount));
901 value_assign(f.d, mpq_denref(tcount));
902 emul(&f, factor);
903 reduce(factor, num_p, pden, options);
904 free_evalue_refs(&f);
906 return;
908 } else {
909 for (k = 0; k < len; ++k) {
910 if (den_s[k] == 0 || den_p[k] == 0)
911 continue;
913 zz2value(den_s[k], tz);
914 dpoly pd(no_param-1, tz, 1);
916 int l;
917 for (l = 0; l < k; ++l)
918 if (den_r[l] == den_r[k])
919 break;
921 if (r == 0)
922 r = new dpoly_r(n, pd, l, len);
923 else {
924 dpoly_r *nr = new dpoly_r(r, pd, l, len);
925 delete r;
926 r = nr;
930 dpoly_r *rc = r->div(D);
931 delete r;
932 r = rc;
933 if (E_num(0, dim) == 0) {
934 int common = pden.NumRows();
935 dpoly_r_term_list& final = r->c[r->len-1];
936 int rows;
937 evalue t;
938 evalue f;
939 value_init(f.d);
940 value_init(f.x.n);
941 zz2value(r->denom, f.d);
942 dpoly_r_term_list::iterator j;
943 for (j = final.begin(); j != final.end(); ++j) {
944 if ((*j)->coeff == 0)
945 continue;
946 rows = common;
947 for (int k = 0; k < r->dim; ++k) {
948 int n = (*j)->powers[k];
949 if (n == 0)
950 continue;
951 pden.SetDims(rows+n, pden.NumCols());
952 for (int l = 0; l < n; ++l)
953 pden[rows+l] = den_r[k];
954 rows += n;
956 value_init(t.d);
957 evalue_copy(&t, factor);
958 zz2value((*j)->coeff, f.x.n);
959 emul(&f, &t);
960 reduce(&t, num_p, pden, options);
961 free_evalue_refs(&t);
963 free_evalue_refs(&f);
964 } else {
965 ie_cum cum(factor, E_num(0, dim), r);
966 cum.cumulate(options);
968 int common = pden.NumRows();
969 int rows;
970 for (int j = 0; j < cum.terms.size(); ++j) {
971 rows = common;
972 pden.SetDims(rows, pden.NumCols());
973 for (int k = 0; k < r->dim; ++k) {
974 int n = cum.terms[j]->powers[k];
975 if (n == 0)
976 continue;
977 pden.SetDims(rows+n, pden.NumCols());
978 for (int l = 0; l < n; ++l)
979 pden[rows+l] = den_r[k];
980 rows += n;
982 reduce(cum.terms[j]->E, num_p, pden, options);
983 free_evalue_refs(cum.terms[j]->E);
984 delete cum.terms[j]->E;
985 delete cum.terms[j];
988 delete r;
992 static int type_offset(enode *p)
994 return p->type == fractional ? 1 :
995 p->type == flooring ? 1 : 0;
998 static int edegree(evalue *e)
1000 int d = 0;
1001 enode *p;
1003 if (value_notzero_p(e->d))
1004 return 0;
1006 p = e->x.p;
1007 int i = type_offset(p);
1008 if (p->size-i-1 > d)
1009 d = p->size - i - 1;
1010 for (; i < p->size; i++) {
1011 int d2 = edegree(&p->arr[i]);
1012 if (d2 > d)
1013 d = d2;
1015 return d;
1018 void ienumerator::handle(const signed_cone& sc, barvinok_options *options)
1020 assert(sc.det == 1);
1021 assert(sc.rays.NumRows() == dim);
1023 lattice_point(V, sc.rays, vertex[0], E_vertex, options);
1025 den = sc.rays;
1027 evalue one;
1028 value_init(one.d);
1029 evalue_set_si(&one, sc.sign, 1);
1030 reduce(&one, vertex, den, options);
1031 free_evalue_refs(&one);
1033 for (int i = 0; i < dim; ++i)
1034 if (E_vertex[i])
1035 evalue_free(E_vertex[i]);
1038 struct bfenumerator : public vertex_decomposer, public bf_base,
1039 public ienumerator_base {
1040 evalue *factor;
1042 bfenumerator(Polyhedron *P, unsigned dim, Param_Polyhedron *PP) :
1043 vertex_decomposer(PP, *this),
1044 bf_base(dim), ienumerator_base(dim, this) {
1045 lower = 0;
1046 factor = NULL;
1049 ~bfenumerator() {
1052 virtual void handle(const signed_cone& sc, barvinok_options *options);
1053 virtual void base(mat_ZZ& factors, bfc_vec& v);
1055 bfc_term_base* new_bf_term(int len) {
1056 bfe_term* t = new bfe_term(len);
1057 return t;
1060 virtual void set_factor(bfc_term_base *t, int k, int change) {
1061 bfe_term* bfet = static_cast<bfe_term *>(t);
1062 factor = bfet->factors[k];
1063 assert(factor != NULL);
1064 bfet->factors[k] = NULL;
1065 if (change)
1066 emul(&mone, factor);
1069 virtual void set_factor(bfc_term_base *t, int k, mpq_t &q, int change) {
1070 bfe_term* bfet = static_cast<bfe_term *>(t);
1071 factor = bfet->factors[k];
1072 assert(factor != NULL);
1073 bfet->factors[k] = NULL;
1075 evalue f;
1076 value_init(f.d);
1077 value_init(f.x.n);
1078 if (change)
1079 value_oppose(f.x.n, mpq_numref(q));
1080 else
1081 value_assign(f.x.n, mpq_numref(q));
1082 value_assign(f.d, mpq_denref(q));
1083 emul(&f, factor);
1084 free_evalue_refs(&f);
1087 virtual void set_factor(bfc_term_base *t, int k, const QQ& c, int change) {
1088 bfe_term* bfet = static_cast<bfe_term *>(t);
1090 factor = new evalue;
1092 evalue f;
1093 value_init(f.d);
1094 value_init(f.x.n);
1095 zz2value(c.n, f.x.n);
1096 if (change)
1097 value_oppose(f.x.n, f.x.n);
1098 zz2value(c.d, f.d);
1100 value_init(factor->d);
1101 evalue_copy(factor, bfet->factors[k]);
1102 emul(&f, factor);
1103 free_evalue_refs(&f);
1106 void set_factor(evalue *f, int change) {
1107 if (change)
1108 emul(&mone, f);
1109 factor = f;
1112 virtual void insert_term(bfc_term_base *t, int i) {
1113 bfe_term* bfet = static_cast<bfe_term *>(t);
1114 int len = t->terms.NumRows()-1; // already increased by one
1116 bfet->factors.resize(len+1);
1117 for (int j = len; j > i; --j) {
1118 bfet->factors[j] = bfet->factors[j-1];
1119 t->terms[j] = t->terms[j-1];
1121 bfet->factors[i] = factor;
1122 factor = NULL;
1125 virtual void update_term(bfc_term_base *t, int i) {
1126 bfe_term* bfet = static_cast<bfe_term *>(t);
1128 eadd(factor, bfet->factors[i]);
1129 free_evalue_refs(factor);
1130 delete factor;
1133 virtual bool constant_vertex(int dim) { return E_num(0, dim) == 0; }
1135 virtual void cum(bf_reducer *bfr, bfc_term_base *t, int k, dpoly_r *r,
1136 barvinok_options *options);
1139 enumerator_base *enumerator_base::create(Polyhedron *P, unsigned dim,
1140 Param_Polyhedron *PP,
1141 barvinok_options *options)
1143 enumerator_base *eb;
1145 if (options->incremental_specialization == BV_SPECIALIZATION_BF)
1146 eb = new bfenumerator(P, dim, PP);
1147 else if (options->incremental_specialization == BV_SPECIALIZATION_DF)
1148 eb = new ienumerator(P, dim, PP);
1149 else
1150 eb = new enumerator(P, dim, PP);
1152 return eb;
1155 struct bfe_cum : public cumulator {
1156 bfenumerator *bfe;
1157 bfc_term_base *told;
1158 int k;
1159 bf_reducer *bfr;
1161 bfe_cum(evalue *factor, evalue *v, dpoly_r *r, bf_reducer *bfr,
1162 bfc_term_base *t, int k, bfenumerator *e) :
1163 cumulator(factor, v, r), told(t), k(k),
1164 bfr(bfr), bfe(e) {
1167 virtual void add_term(const vector<int>& powers, evalue *f2);
1170 void bfe_cum::add_term(const vector<int>& powers, evalue *f2)
1172 bfr->update_powers(powers);
1174 bfc_term_base * t = bfe->find_bfc_term(bfr->vn, bfr->npowers, bfr->nnf);
1175 bfe->set_factor(f2, bfr->l_changes % 2);
1176 bfe->add_term(t, told->terms[k], bfr->l_extra_num);
1179 void bfenumerator::cum(bf_reducer *bfr, bfc_term_base *t, int k,
1180 dpoly_r *r, barvinok_options *options)
1182 bfe_term* bfet = static_cast<bfe_term *>(t);
1183 bfe_cum cum(bfet->factors[k], E_num(0, bfr->d), r, bfr, t, k, this);
1184 cum.cumulate(options);
1187 void bfenumerator::base(mat_ZZ& factors, bfc_vec& v)
1189 for (int i = 0; i < v.size(); ++i) {
1190 assert(v[i]->terms.NumRows() == 1);
1191 evalue *factor = static_cast<bfe_term *>(v[i])->factors[0];
1192 eadd(factor, vE[vert]);
1193 delete v[i];
1197 void bfenumerator::handle(const signed_cone& sc, barvinok_options *options)
1199 assert(sc.det == 1);
1200 assert(sc.rays.NumRows() == enumerator_base::dim);
1202 bfe_term* t = new bfe_term(enumerator_base::dim);
1203 vector< bfc_term_base * > v;
1204 v.push_back(t);
1206 t->factors.resize(1);
1208 t->terms.SetDims(1, enumerator_base::dim);
1209 lattice_point(V, sc.rays, t->terms[0], E_vertex, options);
1211 // the elements of factors are always lexpositive
1212 mat_ZZ factors;
1213 int s = setup_factors(sc.rays, factors, t, sc.sign);
1215 t->factors[0] = new evalue;
1216 value_init(t->factors[0]->d);
1217 evalue_set_si(t->factors[0], s, 1);
1218 reduce(factors, v, options);
1220 for (int i = 0; i < enumerator_base::dim; ++i)
1221 if (E_vertex[i])
1222 evalue_free(E_vertex[i]);
1225 static evalue* barvinok_enumerate_ev_f(Polyhedron *P, Polyhedron* C,
1226 barvinok_options *options);
1228 /* Destroys C */
1229 static evalue* barvinok_enumerate_cst(Polyhedron *P, Polyhedron* C,
1230 struct barvinok_options *options)
1232 evalue *eres;
1234 if (emptyQ2(C)) {
1235 Polyhedron_Free(C);
1236 return evalue_zero();
1239 ALLOC(evalue, eres);
1240 value_init(eres->d);
1241 value_set_si(eres->d, 0);
1242 eres->x.p = new_enode(partition, 2, C->Dimension);
1243 EVALUE_SET_DOMAIN(eres->x.p->arr[0],
1244 DomainConstraintSimplify(C, options->MaxRays));
1245 value_set_si(eres->x.p->arr[1].d, 1);
1246 value_init(eres->x.p->arr[1].x.n);
1247 if (emptyQ2(P))
1248 value_set_si(eres->x.p->arr[1].x.n, 0);
1249 else
1250 barvinok_count_with_options(P, &eres->x.p->arr[1].x.n, options);
1251 if (value_mone_p(eres->x.p->arr[1].x.n)) {
1252 value_clear(eres->x.p->arr[1].x.n);
1253 value_set_si(eres->x.p->arr[1].d, -2); /* NaN */
1256 return eres;
1259 static evalue* enumerate(Polyhedron *P, Polyhedron* C,
1260 struct barvinok_options *options)
1262 Polyhedron *next;
1263 Polyhedron *Porig = P;
1264 Polyhedron *Corig = C;
1265 Polyhedron *CEq = NULL, *rVD;
1266 int r = 0;
1267 unsigned nparam = C->Dimension;
1268 evalue *eres;
1269 Matrix *CP = NULL;
1271 evalue factor;
1272 value_init(factor.d);
1273 evalue_set_si(&factor, 1, 1);
1275 /* for now */
1276 POL_ENSURE_FACETS(P);
1277 POL_ENSURE_VERTICES(P);
1278 POL_ENSURE_FACETS(C);
1279 POL_ENSURE_VERTICES(C);
1281 if (C->Dimension == 0 || emptyQ(P) || emptyQ(C)) {
1282 constant:
1283 if (CEq == Porig)
1284 CEq = Polyhedron_Copy(CEq);
1285 eres = barvinok_enumerate_cst(P, CEq ? CEq : Polyhedron_Copy(C), options);
1286 out:
1287 if (CP) {
1288 evalue_backsubstitute(eres, CP, options->MaxRays);
1289 Matrix_Free(CP);
1292 emul(&factor, eres);
1293 if (options->approx->method == BV_APPROX_DROP) {
1294 if (options->approx->approximation == BV_APPROX_SIGN_UPPER)
1295 evalue_frac2polynomial(eres, 1, options->MaxRays);
1296 if (options->approx->approximation == BV_APPROX_SIGN_LOWER)
1297 evalue_frac2polynomial(eres, -1, options->MaxRays);
1298 if (options->approx->approximation == BV_APPROX_SIGN_APPROX)
1299 evalue_frac2polynomial(eres, 0, options->MaxRays);
1301 reduce_evalue(eres);
1302 free_evalue_refs(&factor);
1303 if (P != Porig)
1304 Domain_Free(P);
1305 if (C != Corig)
1306 Polyhedron_Free(C);
1308 return eres;
1310 if (Polyhedron_is_unbounded(P, nparam, options->MaxRays))
1311 goto constant;
1313 if (P->Dimension == nparam) {
1314 CEq = DomainIntersection(P, C, options->MaxRays);
1315 P = Universe_Polyhedron(0);
1316 goto constant;
1318 if (P->NbEq != 0 || C->NbEq != 0) {
1319 Polyhedron *Q = P;
1320 Polyhedron *D = C;
1321 remove_all_equalities(&P, &C, &CP, NULL, nparam, options->MaxRays);
1322 if (C != D && D != Corig)
1323 Polyhedron_Free(D);
1324 if (P != Q && Q != Porig)
1325 Domain_Free(Q);
1326 eres = enumerate(P, C, options);
1327 goto out;
1330 Polyhedron *T = Polyhedron_Factor(P, nparam, NULL, options->MaxRays);
1331 if (T || (P->Dimension == nparam+1)) {
1332 Polyhedron *C2 = C;
1333 Polyhedron *FC = Factor_Context(T ? T : P, nparam, options->MaxRays);
1334 C = DomainIntersection(C, FC, options->MaxRays);
1335 if (C2 != Corig)
1336 Polyhedron_Free(C2);
1337 Polyhedron_Free(FC);
1339 if (T) {
1340 if (P != Porig)
1341 Polyhedron_Free(P);
1342 P = T;
1343 if (T->Dimension == C->Dimension) {
1344 P = T->next;
1345 T->next = NULL;
1346 Polyhedron_Free(T);
1350 next = P->next;
1351 P->next = NULL;
1352 eres = barvinok_enumerate_ev_f(P, C, options);
1353 P->next = next;
1355 if (P->next) {
1356 Polyhedron *Q;
1357 evalue *f;
1359 for (Q = P->next; Q; Q = Q->next) {
1360 Polyhedron *next = Q->next;
1361 Q->next = NULL;
1363 f = barvinok_enumerate_ev_f(Q, C, options);
1364 emul(f, eres);
1365 evalue_free(f);
1367 Q->next = next;
1371 goto out;
1374 evalue* barvinok_enumerate_with_options(Polyhedron *P, Polyhedron* C,
1375 struct barvinok_options *options)
1377 Polyhedron *next, *Cnext, *C1;
1378 Polyhedron *Corig = C;
1379 evalue *eres;
1381 if (P->next)
1382 fprintf(stderr,
1383 "barvinok_enumerate: input is a union; only first polyhedron is enumerated\n");
1385 if (C->next)
1386 fprintf(stderr,
1387 "barvinok_enumerate: context is a union; only first polyhedron is considered\n");
1389 Cnext = C->next;
1390 C->next = NULL;
1391 C1 = Polyhedron_Project(P, C->Dimension);
1392 C = DomainIntersection(C, C1, options->MaxRays);
1393 Polyhedron_Free(C1);
1394 next = P->next;
1395 P->next = NULL;
1397 if (options->approx->method == BV_APPROX_BERNOULLI ||
1398 options->summation == BV_SUM_BERNOULLI) {
1399 int summation = options->summation;
1400 options->summation = BV_SUM_BERNOULLI;
1401 eres = barvinok_summate_unweighted(P, C, options);
1402 options->summation = summation;
1403 } else
1404 eres = enumerate(P, C, options);
1405 Domain_Free(C);
1407 P->next= next;
1408 Corig->next = Cnext;
1410 return eres;
1413 evalue* barvinok_enumerate_ev(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
1415 evalue *E;
1416 barvinok_options *options = barvinok_options_new_with_defaults();
1417 options->MaxRays = MaxRays;
1418 E = barvinok_enumerate_with_options(P, C, options);
1419 barvinok_options_free(options);
1420 return E;
1423 evalue *Param_Polyhedron_Enumerate(Param_Polyhedron *PP, Polyhedron *P,
1424 Polyhedron *C,
1425 struct barvinok_options *options)
1427 evalue *eres;
1428 Param_Domain *D;
1429 unsigned nparam = C->Dimension;
1430 unsigned dim = P->Dimension - nparam;
1432 int nd;
1433 for (nd = 0, D=PP->D; D; ++nd, D=D->next);
1434 evalue_section *s = new evalue_section[nd];
1436 enumerator_base *et = NULL;
1437 try_again:
1438 if (et)
1439 delete et;
1441 et = enumerator_base::create(P, dim, PP, options);
1443 Polyhedron *TC = true_context(P, C, options->MaxRays);
1444 FORALL_REDUCED_DOMAIN(PP, TC, nd, options, i, D, rVD)
1445 Param_Vertices *V;
1447 s[i].E = evalue_zero();
1448 s[i].D = rVD;
1450 FORALL_PVertex_in_ParamPolyhedron(V,D,PP) // _i is internal counter
1451 if (!et->vE[_i])
1452 try {
1453 et->decompose_at(V, _i, options);
1454 } catch (OrthogonalException &e) {
1455 FORALL_REDUCED_DOMAIN_RESET;
1456 for (; i >= 0; --i) {
1457 evalue_free(s[i].E);
1458 Domain_Free(s[i].D);
1460 goto try_again;
1462 eadd(et->vE[_i] , s[i].E);
1463 END_FORALL_PVertex_in_ParamPolyhedron;
1464 evalue_range_reduction_in_domain(s[i].E, rVD);
1465 END_FORALL_REDUCED_DOMAIN
1466 Polyhedron_Free(TC);
1468 delete et;
1469 eres = evalue_from_section_array(s, nd);
1470 delete [] s;
1472 return eres;
1475 static evalue* barvinok_enumerate_ev_f(Polyhedron *P, Polyhedron* C,
1476 barvinok_options *options)
1478 unsigned nparam = C->Dimension;
1479 bool do_scale = options->approx->method == BV_APPROX_SCALE;
1481 if (options->summation == BV_SUM_EULER)
1482 return barvinok_summate_unweighted(P, C, options);
1484 if (options->approx->method == BV_APPROX_VOLUME)
1485 return Param_Polyhedron_Volume(P, C, options);
1487 if (P->Dimension - nparam == 1 && !do_scale)
1488 return ParamLine_Length(P, C, options);
1490 Param_Polyhedron *PP = NULL;
1491 evalue *eres;
1493 if (do_scale) {
1494 eres = scale_bound(P, C, options);
1495 if (eres)
1496 return eres;
1499 PP = Polyhedron2Param_Polyhedron(P, C, options);
1501 if (do_scale)
1502 eres = scale(PP, P, C, options);
1503 else
1504 eres = Param_Polyhedron_Enumerate(PP, P, C, options);
1506 if (PP)
1507 Param_Polyhedron_Free(PP);
1509 return eres;
1512 Enumeration* barvinok_enumerate(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
1514 evalue *EP = barvinok_enumerate_ev(P, C, MaxRays);
1516 return partition2enumeration(EP);
1519 evalue* barvinok_enumerate_union(Polyhedron *D, Polyhedron* C, unsigned MaxRays)
1521 evalue *EP;
1522 gen_fun *gf = barvinok_enumerate_union_series(D, C, MaxRays);
1523 EP = *gf;
1524 delete gf;
1525 return EP;
1528 static __isl_give isl_pw_qpolynomial *basic_set_card(
1529 __isl_take isl_basic_set *bset)
1531 isl_ctx *ctx;
1532 isl_space *dim;
1533 isl_pw_qpolynomial *pwqp;
1534 unsigned nparam = isl_basic_set_dim(bset, isl_dim_param);
1535 Polyhedron *U = Universe_Polyhedron(nparam);
1536 Polyhedron *P;
1537 evalue *E;
1538 barvinok_options *options;
1539 int options_allocated = 0;
1541 ctx = isl_basic_set_get_ctx(bset);
1542 options = isl_ctx_peek_barvinok_options(ctx);
1543 if (!options) {
1544 options = barvinok_options_new_with_defaults();
1545 options_allocated = 1;
1548 dim = isl_basic_set_get_space(bset);
1549 dim = isl_space_domain(dim);
1551 P = isl_basic_set_to_polylib(bset);
1552 E = enumerate(P, U, options);
1554 pwqp = isl_pw_qpolynomial_from_evalue(dim, E);
1555 isl_basic_set_free(bset);
1557 evalue_free(E);
1558 Polyhedron_Free(P);
1559 Polyhedron_Free(U);
1560 if (options_allocated)
1561 barvinok_options_free(options);
1563 return pwqp;
1566 static int basic_map_card(__isl_take isl_basic_map *bmap, void *user)
1568 isl_pw_qpolynomial **sum = (isl_pw_qpolynomial **)user;
1569 isl_pw_qpolynomial *pwqp;
1570 unsigned nparam = isl_basic_map_dim(bmap, isl_dim_param);
1571 unsigned n_in = isl_basic_map_dim(bmap, isl_dim_in);
1572 isl_space *target_dim;
1573 isl_basic_set *bset;
1575 target_dim = isl_basic_map_get_space(bmap);
1576 target_dim = isl_space_domain(target_dim);
1578 bmap = isl_basic_map_move_dims(bmap, isl_dim_param, nparam,
1579 isl_dim_in, 0, n_in);
1581 bset = isl_basic_map_range(bmap);
1582 bset = isl_basic_set_lift(bset);
1583 pwqp = isl_basic_set_multiplicative_call(bset, &basic_set_card);
1585 pwqp = isl_pw_qpolynomial_move_dims(pwqp, isl_dim_in, 0,
1586 isl_dim_param, nparam, n_in);
1587 pwqp = isl_pw_qpolynomial_reset_domain_space(pwqp, target_dim);
1588 *sum = isl_pw_qpolynomial_add(*sum, pwqp);
1590 return 0;
1593 static __isl_give isl_pw_qpolynomial *card_as_sum(__isl_take isl_map *map,
1594 barvinok_options *options)
1596 isl_ctx *ctx;
1597 isl_val *one;
1598 isl_space *dim;
1599 isl_set *set;
1600 isl_qpolynomial *qp;
1601 isl_pw_qpolynomial *pwqp;
1602 int summation = options->summation;
1604 if (!map)
1605 return NULL;
1607 options->summation = BV_SUM_BERNOULLI;
1609 set = isl_map_wrap(map);
1610 dim = isl_set_get_space(set);
1611 ctx = isl_map_get_ctx(map);
1612 one = isl_val_one(ctx);
1613 qp = isl_qpolynomial_val_on_domain(dim, one);
1615 pwqp = isl_pw_qpolynomial_alloc(set, qp);
1616 pwqp = isl_pw_qpolynomial_sum(pwqp);
1618 options->summation = summation;
1620 return pwqp;
1623 __isl_give isl_pw_qpolynomial *isl_map_card(__isl_take isl_map *map)
1625 isl_ctx *ctx;
1626 isl_space *dim;
1627 isl_pw_qpolynomial *sum;
1628 barvinok_options *options;
1630 ctx = isl_map_get_ctx(map);
1631 options = isl_ctx_peek_barvinok_options(ctx);
1632 if (options &&
1633 (options->approx->method == BV_APPROX_BERNOULLI ||
1634 options->summation == BV_SUM_BERNOULLI))
1635 return card_as_sum(map, options);
1637 dim = isl_map_get_space(map);
1638 dim = isl_space_domain(dim);
1639 dim = isl_space_from_domain(dim);
1640 dim = isl_space_add_dims(dim, isl_dim_out, 1);
1641 sum = isl_pw_qpolynomial_zero(dim);
1643 map = isl_map_make_disjoint(map);
1644 map = isl_map_compute_divs(map);
1646 if (isl_map_foreach_basic_map(map, &basic_map_card, &sum) < 0)
1647 goto error;
1649 isl_map_free(map);
1651 return sum;
1652 error:
1653 isl_map_free(map);
1654 isl_pw_qpolynomial_free(sum);
1655 return NULL;
1658 __isl_give isl_pw_qpolynomial *isl_set_card(__isl_take isl_set *set)
1660 isl_pw_qpolynomial *pwqp;
1661 pwqp = isl_map_card(isl_map_from_range(set));
1662 pwqp = isl_pw_qpolynomial_project_domain_on_params(pwqp);
1663 return pwqp;
1666 __isl_give isl_pw_qpolynomial *isl_basic_map_card(__isl_take isl_basic_map *bmap)
1668 return isl_map_card(isl_map_from_basic_map(bmap));
1671 __isl_give isl_pw_qpolynomial *isl_basic_set_card(__isl_take isl_basic_set *bset)
1673 isl_pw_qpolynomial *pwqp;
1674 pwqp = isl_basic_map_card(isl_basic_map_from_range(bset));
1675 pwqp = isl_pw_qpolynomial_project_domain_on_params(pwqp);
1676 return pwqp;
1679 static int set_card(__isl_take isl_set *set, void *user)
1681 isl_union_pw_qpolynomial **res = (isl_union_pw_qpolynomial **)user;
1682 isl_pw_qpolynomial *pwqp;
1684 pwqp = isl_set_card(set);
1685 *res = isl_union_pw_qpolynomial_add_pw_qpolynomial(*res, pwqp);
1687 return 0;
1690 __isl_give isl_union_pw_qpolynomial *isl_union_set_card(
1691 __isl_take isl_union_set *uset)
1693 isl_space *dim;
1694 isl_union_pw_qpolynomial *res;
1696 dim = isl_union_set_get_space(uset);
1697 res = isl_union_pw_qpolynomial_zero(dim);
1698 if (isl_union_set_foreach_set(uset, &set_card, &res) < 0)
1699 goto error;
1700 isl_union_set_free(uset);
1702 return res;
1703 error:
1704 isl_union_set_free(uset);
1705 isl_union_pw_qpolynomial_free(res);
1706 return NULL;
1709 static int map_card(__isl_take isl_map *map, void *user)
1711 isl_union_pw_qpolynomial **res = (isl_union_pw_qpolynomial **)user;
1712 isl_pw_qpolynomial *pwqp;
1714 pwqp = isl_map_card(map);
1715 *res = isl_union_pw_qpolynomial_add_pw_qpolynomial(*res, pwqp);
1717 return 0;
1720 __isl_give isl_union_pw_qpolynomial *isl_union_map_card(
1721 __isl_take isl_union_map *umap)
1723 isl_space *dim;
1724 isl_union_pw_qpolynomial *res;
1726 dim = isl_union_map_get_space(umap);
1727 res = isl_union_pw_qpolynomial_zero(dim);
1728 if (isl_union_map_foreach_map(umap, &map_card, &res) < 0)
1729 goto error;
1730 isl_union_map_free(umap);
1732 return res;
1733 error:
1734 isl_union_map_free(umap);
1735 isl_union_pw_qpolynomial_free(res);
1736 return NULL;