isl_basic_map_new_from_polylib: add memory management annotations
[barvinok.git] / barvinok.cc
blob9a20cf59e2ba492918cc387bee2790d4c96867e9
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/ZZ.h>
9 #include <NTL/vec_ZZ.h>
10 #include <NTL/mat_ZZ.h>
11 #include <isl/val.h>
12 #include <isl/space.h>
13 #include <isl/set.h>
14 #include <isl/map.h>
15 #include <isl/union_set.h>
16 #include <isl/union_map.h>
17 #include <isl/polynomial.h>
18 #include <isl_set_polylib.h>
19 #include <barvinok/polylib.h>
20 #include <barvinok/util.h>
21 #include <barvinok/evalue.h>
22 #include "config.h"
23 #include <barvinok/barvinok.h>
24 #include <barvinok/genfun.h>
25 #include <barvinok/options.h>
26 #include <barvinok/sample.h>
27 #include "bfcounter.h"
28 #include "conversion.h"
29 #include "counter.h"
30 #include "decomposer.h"
31 #include "euler.h"
32 #include "lattice_point.h"
33 #include "laurent.h"
34 #include "reduce_domain.h"
35 #include "remove_equalities.h"
36 #include "scale.h"
37 #include "volume.h"
38 #include "bernoulli.h"
39 #include "param_util.h"
40 #include "summate.h"
42 using namespace NTL;
43 using std::cerr;
44 using std::cout;
45 using std::endl;
46 using std::vector;
47 using std::deque;
48 using std::string;
49 using std::ostringstream;
51 #define ALLOC(t,p) p = (t*)malloc(sizeof(*p))
53 class dpoly_n {
54 public:
55 Matrix *coeff;
56 ~dpoly_n() {
57 Matrix_Free(coeff);
59 dpoly_n(int d) {
60 Value d0, one;
61 value_init(d0);
62 value_init(one);
63 value_set_si(one, 1);
64 coeff = Matrix_Alloc(d+1, d+1+1);
65 value_set_si(coeff->p[0][0], 1);
66 value_set_si(coeff->p[0][d+1], 1);
67 for (int i = 1; i <= d; ++i) {
68 value_multiply(coeff->p[i][0], coeff->p[i-1][0], d0);
69 Vector_Combine(coeff->p[i-1], coeff->p[i-1]+1, coeff->p[i]+1,
70 one, d0, i);
71 value_set_si(coeff->p[i][d+1], i);
72 value_multiply(coeff->p[i][d+1], coeff->p[i][d+1], coeff->p[i-1][d+1]);
73 value_decrement(d0, d0);
75 value_clear(d0);
76 value_clear(one);
78 void div(dpoly& d, Vector *count, int sign) {
79 int len = coeff->NbRows;
80 Matrix * c = Matrix_Alloc(coeff->NbRows, coeff->NbColumns);
81 Value tmp;
82 value_init(tmp);
83 for (int i = 0; i < len; ++i) {
84 Vector_Copy(coeff->p[i], c->p[i], len+1);
85 for (int j = 1; j <= i; ++j) {
86 value_multiply(tmp, d.coeff->p[j], c->p[i][len]);
87 value_oppose(tmp, tmp);
88 Vector_Combine(c->p[i], c->p[i-j], c->p[i],
89 c->p[i-j][len], tmp, len);
90 value_multiply(c->p[i][len], c->p[i][len], c->p[i-j][len]);
92 value_multiply(c->p[i][len], c->p[i][len], d.coeff->p[0]);
94 if (sign == -1) {
95 value_set_si(tmp, -1);
96 Vector_Scale(c->p[len-1], count->p, tmp, len);
97 value_assign(count->p[len], c->p[len-1][len]);
98 } else
99 Vector_Copy(c->p[len-1], count->p, len+1);
100 Vector_Normalize(count->p, len+1);
101 value_clear(tmp);
102 Matrix_Free(c);
106 struct bfe_term : public bfc_term_base {
107 vector<evalue *> factors;
109 bfe_term(int len) : bfc_term_base(len) {
112 ~bfe_term() {
113 for (int i = 0; i < factors.size(); ++i) {
114 if (!factors[i])
115 continue;
116 free_evalue_refs(factors[i]);
117 delete factors[i];
122 static void print_int_vector(int *v, int len, const char *name)
124 cerr << name << endl;
125 for (int j = 0; j < len; ++j) {
126 cerr << v[j] << " ";
128 cerr << endl;
131 static void print_bfc_terms(mat_ZZ& factors, bfc_vec& v)
132 __attribute__((unused));
133 static void print_bfc_terms(mat_ZZ& factors, bfc_vec& v)
135 cerr << endl;
136 cerr << "factors" << endl;
137 cerr << factors << endl;
138 for (int i = 0; i < v.size(); ++i) {
139 cerr << "term: " << i << endl;
140 print_int_vector(v[i]->powers, factors.NumRows(), "powers");
141 cerr << "terms" << endl;
142 cerr << v[i]->terms << endl;
143 bfc_term* bfct = static_cast<bfc_term *>(v[i]);
144 cerr << bfct->c << endl;
148 static void print_bfe_terms(mat_ZZ& factors, bfc_vec& v)
149 __attribute__((unused));
150 static void print_bfe_terms(mat_ZZ& factors, bfc_vec& v)
152 cerr << endl;
153 cerr << "factors" << endl;
154 cerr << factors << endl;
155 for (int i = 0; i < v.size(); ++i) {
156 cerr << "term: " << i << endl;
157 print_int_vector(v[i]->powers, factors.NumRows(), "powers");
158 cerr << "terms" << endl;
159 cerr << v[i]->terms << endl;
160 bfe_term* bfet = static_cast<bfe_term *>(v[i]);
161 for (int j = 0; j < v[i]->terms.NumRows(); ++j) {
162 const char * test[] = {"a", "b"};
163 print_evalue(stderr, bfet->factors[j], test);
164 fprintf(stderr, "\n");
169 struct bfcounter : public bfcounter_base {
170 mpq_t count;
171 Value tz;
173 bfcounter(unsigned dim) : bfcounter_base(dim) {
174 mpq_init(count);
175 lower = 1;
176 value_init(tz);
178 ~bfcounter() {
179 mpq_clear(count);
180 value_clear(tz);
182 virtual void base(mat_ZZ& factors, bfc_vec& v);
183 virtual void get_count(Value *result) {
184 assert(value_one_p(&count[0]._mp_den));
185 value_assign(*result, &count[0]._mp_num);
189 void bfcounter::base(mat_ZZ& factors, bfc_vec& v)
191 unsigned nf = factors.NumRows();
193 for (int i = 0; i < v.size(); ++i) {
194 bfc_term* bfct = static_cast<bfc_term *>(v[i]);
195 int total_power = 0;
196 // factor is always positive, so we always
197 // change signs
198 for (int k = 0; k < nf; ++k)
199 total_power += v[i]->powers[k];
201 int j;
202 for (j = 0; j < nf; ++j)
203 if (v[i]->powers[j] > 0)
204 break;
206 zz2value(factors[j][0], tz);
207 dpoly D(total_power, tz, 1);
208 for (int k = 1; k < v[i]->powers[j]; ++k) {
209 zz2value(factors[j][0], tz);
210 dpoly fact(total_power, tz, 1);
211 D *= fact;
213 for ( ; ++j < nf; )
214 for (int k = 0; k < v[i]->powers[j]; ++k) {
215 zz2value(factors[j][0], tz);
216 dpoly fact(total_power, tz, 1);
217 D *= fact;
220 for (int k = 0; k < v[i]->terms.NumRows(); ++k) {
221 zz2value(v[i]->terms[k][0], tz);
222 dpoly n(total_power, tz);
223 mpq_set_si(tcount, 0, 1);
224 n.div(D, tcount, 1);
225 if (total_power % 2)
226 bfct->c[k].n = -bfct->c[k].n;
227 zz2value(bfct->c[k].n, tn);
228 zz2value(bfct->c[k].d, td);
230 mpz_mul(mpq_numref(tcount), mpq_numref(tcount), tn);
231 mpz_mul(mpq_denref(tcount), mpq_denref(tcount), td);
232 mpq_canonicalize(tcount);
233 mpq_add(count, count, tcount);
235 delete v[i];
240 /* Check whether the polyhedron is unbounded and if so,
241 * check whether it has any (and therefore an infinite number of)
242 * integer points.
243 * If one of the vertices is integer, then we are done.
244 * Otherwise, transform the polyhedron such that one of the rays
245 * is the first unit vector and cut it off at a height that ensures
246 * that if the whole polyhedron has any points, then the remaining part
247 * has integer points. In particular we add the largest coefficient
248 * of a ray to the highest vertex (rounded up).
250 static bool Polyhedron_is_infinite(Polyhedron *P, Value* result,
251 barvinok_options *options)
253 int r = 0;
254 Matrix *M, *M2;
255 Value c, tmp;
256 Value g;
257 bool first;
258 Vector *v;
259 Value offset, size;
260 Polyhedron *R;
262 if (P->NbBid == 0)
263 for (; r < P->NbRays; ++r)
264 if (value_zero_p(P->Ray[r][P->Dimension+1]))
265 break;
266 if (P->NbBid == 0 && r == P->NbRays)
267 return false;
269 if (options->count_sample_infinite) {
270 Vector *sample;
272 sample = Polyhedron_Sample(P, options);
273 if (!sample)
274 value_set_si(*result, 0);
275 else {
276 value_set_si(*result, -1);
277 Vector_Free(sample);
279 return true;
282 for (int i = 0; i < P->NbRays; ++i)
283 if (value_one_p(P->Ray[i][1+P->Dimension])) {
284 value_set_si(*result, -1);
285 return true;
288 value_init(g);
289 M = Matrix_Alloc(P->Dimension+1, P->Dimension+1);
290 Vector_Gcd(P->Ray[r]+1, P->Dimension, &g);
291 Vector_AntiScale(P->Ray[r]+1, M->p[0], g, P->Dimension+1);
292 int ok = unimodular_complete(M, 1);
293 assert(ok);
294 value_set_si(M->p[P->Dimension][P->Dimension], 1);
295 M2 = Transpose(M);
296 Matrix_Free(M);
297 P = Polyhedron_Preimage(P, M2, 0);
298 Matrix_Free(M2);
299 value_clear(g);
301 first = true;
302 value_init(offset);
303 value_init(size);
304 value_init(tmp);
305 value_set_si(size, 0);
307 for (int i = 0; i < P->NbBid; ++i) {
308 value_absolute(tmp, P->Ray[i][1]);
309 if (value_gt(tmp, size))
310 value_assign(size, tmp);
312 for (int i = P->NbBid; i < P->NbRays; ++i) {
313 if (value_zero_p(P->Ray[i][P->Dimension+1])) {
314 if (value_gt(P->Ray[i][1], size))
315 value_assign(size, P->Ray[i][1]);
316 continue;
318 mpz_cdiv_q(tmp, P->Ray[i][1], P->Ray[i][P->Dimension+1]);
319 if (first || value_gt(tmp, offset)) {
320 value_assign(offset, tmp);
321 first = false;
324 value_addto(offset, offset, size);
325 value_clear(size);
326 value_clear(tmp);
328 v = Vector_Alloc(P->Dimension+2);
329 value_set_si(v->p[0], 1);
330 value_set_si(v->p[1], -1);
331 value_assign(v->p[1+P->Dimension], offset);
332 R = AddConstraints(v->p, 1, P, options->MaxRays);
333 Polyhedron_Free(P);
334 P = R;
336 value_clear(offset);
337 Vector_Free(v);
339 value_init(c);
340 barvinok_count_with_options(P, &c, options);
341 Polyhedron_Free(P);
342 if (value_zero_p(c))
343 value_set_si(*result, 0);
344 else
345 value_set_si(*result, -1);
346 value_clear(c);
348 return true;
351 static void evalue2value(evalue *e, Value *v)
353 if (EVALUE_IS_ZERO(*e)) {
354 value_set_si(*v, 0);
355 return;
358 if (value_notzero_p(e->d)) {
359 assert(value_one_p(e->d));
360 value_assign(*v, e->x.n);
361 return;
364 assert(e->x.p->type == partition);
365 assert(e->x.p->size == 2);
366 assert(EVALUE_DOMAIN(e->x.p->arr[0])->Dimension == 0);
367 evalue2value(&e->x.p->arr[1], v);
370 static void barvinok_count_f(Polyhedron *P, Value* result,
371 barvinok_options *options);
373 void barvinok_count_with_options(Polyhedron *P, Value* result,
374 struct barvinok_options *options)
376 int allocated = 0;
377 Polyhedron *Q;
378 bool infinite = false;
380 if (P->next)
381 fprintf(stderr,
382 "barvinok_count: input is a union; only first polyhedron is counted\n");
384 if (emptyQ2(P)) {
385 value_set_si(*result, 0);
386 return;
388 if (P->NbEq != 0) {
389 Q = NULL;
390 do {
391 P = remove_equalities(P, options->MaxRays);
392 if (P)
393 P = DomainConstraintSimplify(P, options->MaxRays);
394 if (Q)
395 Polyhedron_Free(Q);
396 Q = P;
397 } while (P && !emptyQ(P) && P->NbEq != 0);
398 if (!P || emptyQ(P)) {
399 Polyhedron_Free(P);
400 value_set_si(*result, 0);
401 return;
403 allocated = 1;
405 if (Polyhedron_is_infinite(P, result, options)) {
406 if (allocated)
407 Polyhedron_Free(P);
408 return;
410 if (P->Dimension == 0) {
411 /* Test whether the constraints are satisfied */
412 POL_ENSURE_VERTICES(P);
413 value_set_si(*result, !emptyQ(P));
414 if (allocated)
415 Polyhedron_Free(P);
416 return;
418 if (options->summation == BV_SUM_BERNOULLI) {
419 Polyhedron *C = Universe_Polyhedron(0);
420 evalue *sum = barvinok_summate_unweighted(P, C, options);
421 Polyhedron_Free(C);
422 evalue2value(sum, result);
423 evalue_free(sum);
424 return;
426 Q = Polyhedron_Factor(P, 0, NULL, options->MaxRays);
427 if (Q) {
428 if (allocated)
429 Polyhedron_Free(P);
430 P = Q;
431 allocated = 1;
434 barvinok_count_f(P, result, options);
435 if (value_neg_p(*result))
436 infinite = true;
437 if (Q && P->next && value_notzero_p(*result)) {
438 Value factor;
439 value_init(factor);
441 for (Q = P->next; Q; Q = Q->next) {
442 barvinok_count_f(Q, &factor, options);
443 if (value_neg_p(factor)) {
444 infinite = true;
445 continue;
446 } else if (Q->next && value_zero_p(factor)) {
447 value_set_si(*result, 0);
448 break;
450 value_multiply(*result, *result, factor);
453 value_clear(factor);
456 if (allocated)
457 Domain_Free(P);
458 if (infinite)
459 value_set_si(*result, -1);
462 void barvinok_count(Polyhedron *P, Value* result, unsigned NbMaxCons)
464 barvinok_options *options = barvinok_options_new_with_defaults();
465 options->MaxRays = NbMaxCons;
466 barvinok_count_with_options(P, result, options);
467 barvinok_options_free(options);
470 static void barvinok_count_f(Polyhedron *P, Value* result,
471 barvinok_options *options)
473 if (emptyQ2(P)) {
474 value_set_si(*result, 0);
475 return;
478 if (P->Dimension == 1)
479 return Line_Length(P, result);
481 int c = P->NbConstraints;
482 POL_ENSURE_FACETS(P);
483 if (c != P->NbConstraints || P->NbEq != 0) {
484 Polyhedron *next = P->next;
485 P->next = NULL;
486 barvinok_count_with_options(P, result, options);
487 P->next = next;
488 return;
491 POL_ENSURE_VERTICES(P);
493 if (Polyhedron_is_infinite(P, result, options))
494 return;
496 np_base *cnt;
497 if (options->incremental_specialization == BV_SPECIALIZATION_BF)
498 cnt = new bfcounter(P->Dimension);
499 else if (options->incremental_specialization == BV_SPECIALIZATION_DF)
500 cnt = new icounter(P->Dimension);
501 else if (options->incremental_specialization == BV_SPECIALIZATION_TODD)
502 cnt = new tcounter(P->Dimension, options->max_index);
503 else
504 cnt = new counter(P->Dimension, options->max_index);
505 cnt->start(P, options);
507 cnt->get_count(result);
508 delete cnt;
511 typedef evalue * evalue_p;
513 struct enumerator_base {
514 unsigned dim;
515 evalue ** vE;
516 evalue mone;
517 vertex_decomposer *vpd;
519 enumerator_base(unsigned dim, vertex_decomposer *vpd)
521 this->dim = dim;
522 this->vpd = vpd;
524 vE = new evalue_p[vpd->PP->nbV];
525 for (int j = 0; j < vpd->PP->nbV; ++j)
526 vE[j] = 0;
528 value_init(mone.d);
529 evalue_set_si(&mone, -1, 1);
532 void decompose_at(Param_Vertices *V, int _i, barvinok_options *options) {
533 //this->pVD = pVD;
535 vE[_i] = new evalue;
536 value_init(vE[_i]->d);
537 evalue_set_si(vE[_i], 0, 1);
539 vpd->decompose_at_vertex(V, _i, options);
542 virtual ~enumerator_base() {
543 for (int j = 0; j < vpd->PP->nbV; ++j)
544 if (vE[j]) {
545 free_evalue_refs(vE[j]);
546 delete vE[j];
548 delete [] vE;
550 free_evalue_refs(&mone);
553 static enumerator_base *create(Polyhedron *P, unsigned dim,
554 Param_Polyhedron *PP,
555 barvinok_options *options);
558 struct enumerator : public signed_cone_consumer, public vertex_decomposer,
559 public enumerator_base {
560 vec_ZZ lambda;
561 vec_ZZ den;
562 term_info num;
563 Vector *c;
564 mpq_t count;
565 Value tz;
567 enumerator(Polyhedron *P, unsigned dim, Param_Polyhedron *PP) :
568 vertex_decomposer(PP, *this), enumerator_base(dim, this) {
569 randomvector(P, lambda, dim, 0);
570 den.SetLength(dim);
571 c = Vector_Alloc(dim+2);
573 mpq_init(count);
574 value_init(tz);
577 ~enumerator() {
578 mpq_clear(count);
579 Vector_Free(c);
580 value_clear(tz);
583 virtual void handle(const signed_cone& sc, barvinok_options *options);
586 void enumerator::handle(const signed_cone& sc, barvinok_options *options)
588 int sign = sc.sign;
589 assert(sc.rays.NumRows() == dim);
590 for (int k = 0; k < dim; ++k) {
591 if (lambda * sc.rays[k] == 0)
592 throw Orthogonal;
595 lattice_point(V, sc.rays, lambda, &num, sc.det, options);
596 den = sc.rays * lambda;
598 if (dim % 2)
599 sign = -sign;
601 zz2value(den[0], tz);
602 dpoly n(dim, tz, 1);
603 for (int k = 1; k < dim; ++k) {
604 zz2value(den[k], tz);
605 dpoly fact(dim, tz, 1);
606 n *= fact;
608 if (num.E != NULL) {
609 dpoly_n d(dim);
610 d.div(n, c, sign);
611 for (unsigned long i = 0; i < sc.det; ++i) {
612 evalue *EV = evalue_polynomial(c, num.E[i]);
613 eadd(EV, vE[vert]);
614 evalue_free(EV);
615 evalue_free(num.E[i]);
617 delete [] num.E;
618 } else {
619 mpq_set_si(count, 0, 1);
620 if (num.constant.length() == 1) {
621 zz2value(num.constant[0], tz);
622 dpoly d(dim, tz);
623 d.div(n, count, sign);
624 } else {
625 dpoly_n d(dim);
626 d.div(n, c, sign);
627 Value x, acc;
628 value_init(x);
629 value_init(acc);
630 for (unsigned long i = 0; i < sc.det; ++i) {
631 value_assign(acc, c->p[dim]);
632 zz2value(num.constant[i], x);
633 for (int j = dim-1; j >= 0; --j) {
634 value_multiply(acc, acc, x);
635 value_addto(acc, acc, c->p[j]);
637 value_addto(mpq_numref(count), mpq_numref(count), acc);
639 mpz_set(mpq_denref(count), c->p[dim+1]);
640 value_clear(acc);
641 value_clear(x);
643 evalue EV;
644 value_init(EV.d);
645 evalue_set(&EV, &count[0]._mp_num, &count[0]._mp_den);
646 eadd(&EV, vE[vert]);
647 free_evalue_refs(&EV);
651 struct ienumerator_base : enumerator_base {
652 evalue ** E_vertex;
654 ienumerator_base(unsigned dim, vertex_decomposer *vpd) :
655 enumerator_base(dim,vpd) {
656 E_vertex = new evalue_p[dim];
659 virtual ~ienumerator_base() {
660 delete [] E_vertex;
663 evalue *E_num(int i, int d) {
664 return E_vertex[i + (dim-d)];
668 struct cumulator {
669 evalue *factor;
670 evalue *v;
671 dpoly_r *r;
673 cumulator(evalue *factor, evalue *v, dpoly_r *r) :
674 factor(factor), v(v), r(r) {}
676 void cumulate(barvinok_options *options);
678 virtual void add_term(const vector<int>& powers, evalue *f2) = 0;
679 virtual ~cumulator() {}
682 void cumulator::cumulate(barvinok_options *options)
684 evalue cum; // factor * 1 * E_num[0]/1 * (E_num[0]-1)/2 *...
685 evalue f;
686 evalue t; // E_num[0] - (m-1)
687 evalue *cst;
688 evalue mone;
690 if (options->lookup_table) {
691 value_init(mone.d);
692 evalue_set_si(&mone, -1, 1);
695 value_init(cum.d);
696 evalue_copy(&cum, factor);
697 value_init(f.d);
698 value_init(f.x.n);
699 value_set_si(f.d, 1);
700 value_set_si(f.x.n, 1);
701 value_init(t.d);
702 evalue_copy(&t, v);
704 if (!options->lookup_table) {
705 for (cst = &t; value_zero_p(cst->d); ) {
706 if (cst->x.p->type == fractional)
707 cst = &cst->x.p->arr[1];
708 else
709 cst = &cst->x.p->arr[0];
713 for (int m = 0; m < r->len; ++m) {
714 if (m > 0) {
715 if (m > 1) {
716 value_set_si(f.d, m);
717 emul(&f, &cum);
718 if (!options->lookup_table)
719 value_subtract(cst->x.n, cst->x.n, cst->d);
720 else
721 eadd(&mone, &t);
723 emul(&t, &cum);
725 dpoly_r_term_list& current = r->c[r->len-1-m];
726 dpoly_r_term_list::iterator j;
727 for (j = current.begin(); j != current.end(); ++j) {
728 if ((*j)->coeff == 0)
729 continue;
730 evalue *f2 = new evalue;
731 value_init(f2->d);
732 value_init(f2->x.n);
733 zz2value((*j)->coeff, f2->x.n);
734 zz2value(r->denom, f2->d);
735 emul(&cum, f2);
737 add_term((*j)->powers, f2);
740 free_evalue_refs(&f);
741 free_evalue_refs(&t);
742 free_evalue_refs(&cum);
743 if (options->lookup_table)
744 free_evalue_refs(&mone);
747 struct E_poly_term {
748 vector<int> powers;
749 evalue *E;
752 struct ie_cum : public cumulator {
753 vector<E_poly_term *> terms;
755 ie_cum(evalue *factor, evalue *v, dpoly_r *r) : cumulator(factor, v, r) {}
757 virtual void add_term(const vector<int>& powers, evalue *f2);
760 void ie_cum::add_term(const vector<int>& powers, evalue *f2)
762 int k;
763 for (k = 0; k < terms.size(); ++k) {
764 if (terms[k]->powers == powers) {
765 eadd(f2, terms[k]->E);
766 free_evalue_refs(f2);
767 delete f2;
768 break;
771 if (k >= terms.size()) {
772 E_poly_term *ET = new E_poly_term;
773 ET->powers = powers;
774 ET->E = f2;
775 terms.push_back(ET);
779 struct ienumerator : public signed_cone_consumer, public vertex_decomposer,
780 public ienumerator_base {
781 mat_ZZ den;
782 mat_ZZ vertex;
783 mpq_t tcount;
784 Value tz;
786 ienumerator(unsigned dim, Param_Polyhedron *PP) :
787 vertex_decomposer(PP, *this), ienumerator_base(dim, this) {
788 vertex.SetDims(1, dim);
790 den.SetDims(dim, dim);
791 mpq_init(tcount);
792 value_init(tz);
795 ~ienumerator() {
796 mpq_clear(tcount);
797 value_clear(tz);
800 virtual void handle(const signed_cone& sc, barvinok_options *options);
801 void reduce(evalue *factor, const mat_ZZ& num, const mat_ZZ& den_f,
802 barvinok_options *options);
805 void ienumerator::reduce(evalue *factor, const mat_ZZ& num, const mat_ZZ& den_f,
806 barvinok_options *options)
808 unsigned len = den_f.NumRows(); // number of factors in den
809 unsigned dim = num.NumCols();
810 assert(num.NumRows() == 1);
812 if (dim == 0) {
813 eadd(factor, vE[vert]);
814 return;
817 vec_ZZ den_s;
818 mat_ZZ den_r;
819 vec_ZZ num_s;
820 mat_ZZ num_p;
822 split_one(num, num_s, num_p, den_f, den_s, den_r);
824 vec_ZZ den_p;
825 den_p.SetLength(len);
827 ZZ one;
828 one = 1;
829 normalize(one, num_s, num_p, den_s, den_p, den_r);
830 if (one != 1)
831 emul(&mone, factor);
833 int only_param = 0;
834 int no_param = 0;
835 for (int k = 0; k < len; ++k) {
836 if (den_p[k] == 0)
837 ++no_param;
838 else if (den_s[k] == 0)
839 ++only_param;
841 if (no_param == 0) {
842 reduce(factor, num_p, den_r, options);
843 } else {
844 int k, l;
845 mat_ZZ pden;
846 pden.SetDims(only_param, dim-1);
848 for (k = 0, l = 0; k < len; ++k)
849 if (den_s[k] == 0)
850 pden[l++] = den_r[k];
852 for (k = 0; k < len; ++k)
853 if (den_p[k] == 0)
854 break;
856 zz2value(num_s[0], tz);
857 dpoly n(no_param, tz);
858 zz2value(den_s[k], tz);
859 dpoly D(no_param, tz, 1);
860 for ( ; ++k < len; )
861 if (den_p[k] == 0) {
862 zz2value(den_s[k], tz);
863 dpoly fact(no_param, tz, 1);
864 D *= fact;
867 dpoly_r * r = 0;
868 // if no_param + only_param == len then all powers
869 // below will be all zero
870 if (no_param + only_param == len) {
871 if (E_num(0, dim) != 0)
872 r = new dpoly_r(n, len);
873 else {
874 mpq_set_si(tcount, 0, 1);
875 one = 1;
876 n.div(D, tcount, 1);
878 if (value_notzero_p(mpq_numref(tcount))) {
879 evalue f;
880 value_init(f.d);
881 value_init(f.x.n);
882 value_assign(f.x.n, mpq_numref(tcount));
883 value_assign(f.d, mpq_denref(tcount));
884 emul(&f, factor);
885 reduce(factor, num_p, pden, options);
886 free_evalue_refs(&f);
888 return;
890 } else {
891 for (k = 0; k < len; ++k) {
892 if (den_s[k] == 0 || den_p[k] == 0)
893 continue;
895 zz2value(den_s[k], tz);
896 dpoly pd(no_param-1, tz, 1);
898 int l;
899 for (l = 0; l < k; ++l)
900 if (den_r[l] == den_r[k])
901 break;
903 if (r == 0)
904 r = new dpoly_r(n, pd, l, len);
905 else {
906 dpoly_r *nr = new dpoly_r(r, pd, l, len);
907 delete r;
908 r = nr;
912 dpoly_r *rc = r->div(D);
913 delete r;
914 r = rc;
915 if (E_num(0, dim) == 0) {
916 int common = pden.NumRows();
917 dpoly_r_term_list& final = r->c[r->len-1];
918 int rows;
919 evalue t;
920 evalue f;
921 value_init(f.d);
922 value_init(f.x.n);
923 zz2value(r->denom, f.d);
924 dpoly_r_term_list::iterator j;
925 for (j = final.begin(); j != final.end(); ++j) {
926 if ((*j)->coeff == 0)
927 continue;
928 rows = common;
929 for (int k = 0; k < r->dim; ++k) {
930 int n = (*j)->powers[k];
931 if (n == 0)
932 continue;
933 pden.SetDims(rows+n, pden.NumCols());
934 for (int l = 0; l < n; ++l)
935 pden[rows+l] = den_r[k];
936 rows += n;
938 value_init(t.d);
939 evalue_copy(&t, factor);
940 zz2value((*j)->coeff, f.x.n);
941 emul(&f, &t);
942 reduce(&t, num_p, pden, options);
943 free_evalue_refs(&t);
945 free_evalue_refs(&f);
946 } else {
947 ie_cum cum(factor, E_num(0, dim), r);
948 cum.cumulate(options);
950 int common = pden.NumRows();
951 int rows;
952 for (int j = 0; j < cum.terms.size(); ++j) {
953 rows = common;
954 pden.SetDims(rows, pden.NumCols());
955 for (int k = 0; k < r->dim; ++k) {
956 int n = cum.terms[j]->powers[k];
957 if (n == 0)
958 continue;
959 pden.SetDims(rows+n, pden.NumCols());
960 for (int l = 0; l < n; ++l)
961 pden[rows+l] = den_r[k];
962 rows += n;
964 reduce(cum.terms[j]->E, num_p, pden, options);
965 free_evalue_refs(cum.terms[j]->E);
966 delete cum.terms[j]->E;
967 delete cum.terms[j];
970 delete r;
974 void ienumerator::handle(const signed_cone& sc, barvinok_options *options)
976 assert(sc.det == 1);
977 assert(sc.rays.NumRows() == dim);
979 lattice_point(V, sc.rays, vertex[0], E_vertex, options);
981 den = sc.rays;
983 evalue one;
984 value_init(one.d);
985 evalue_set_si(&one, sc.sign, 1);
986 reduce(&one, vertex, den, options);
987 free_evalue_refs(&one);
989 for (int i = 0; i < dim; ++i)
990 if (E_vertex[i])
991 evalue_free(E_vertex[i]);
994 struct bfenumerator : public vertex_decomposer, public bf_base,
995 public ienumerator_base {
996 evalue *factor;
998 bfenumerator(unsigned dim, Param_Polyhedron *PP) :
999 vertex_decomposer(PP, *this),
1000 bf_base(dim), ienumerator_base(dim, this) {
1001 lower = 0;
1002 factor = NULL;
1005 ~bfenumerator() {
1008 virtual void handle(const signed_cone& sc, barvinok_options *options);
1009 virtual void base(mat_ZZ& factors, bfc_vec& v);
1011 bfc_term_base* new_bf_term(int len) {
1012 bfe_term* t = new bfe_term(len);
1013 return t;
1016 virtual void set_factor(bfc_term_base *t, int k, int change) {
1017 bfe_term* bfet = static_cast<bfe_term *>(t);
1018 factor = bfet->factors[k];
1019 assert(factor != NULL);
1020 bfet->factors[k] = NULL;
1021 if (change)
1022 emul(&mone, factor);
1025 virtual void set_factor(bfc_term_base *t, int k, mpq_t &q, int change) {
1026 bfe_term* bfet = static_cast<bfe_term *>(t);
1027 factor = bfet->factors[k];
1028 assert(factor != NULL);
1029 bfet->factors[k] = NULL;
1031 evalue f;
1032 value_init(f.d);
1033 value_init(f.x.n);
1034 if (change)
1035 value_oppose(f.x.n, mpq_numref(q));
1036 else
1037 value_assign(f.x.n, mpq_numref(q));
1038 value_assign(f.d, mpq_denref(q));
1039 emul(&f, factor);
1040 free_evalue_refs(&f);
1043 virtual void set_factor(bfc_term_base *t, int k, const QQ& c, int change) {
1044 bfe_term* bfet = static_cast<bfe_term *>(t);
1046 factor = new evalue;
1048 evalue f;
1049 value_init(f.d);
1050 value_init(f.x.n);
1051 zz2value(c.n, f.x.n);
1052 if (change)
1053 value_oppose(f.x.n, f.x.n);
1054 zz2value(c.d, f.d);
1056 value_init(factor->d);
1057 evalue_copy(factor, bfet->factors[k]);
1058 emul(&f, factor);
1059 free_evalue_refs(&f);
1062 void set_factor(evalue *f, int change) {
1063 if (change)
1064 emul(&mone, f);
1065 factor = f;
1068 virtual void insert_term(bfc_term_base *t, int i) {
1069 bfe_term* bfet = static_cast<bfe_term *>(t);
1070 int len = t->terms.NumRows()-1; // already increased by one
1072 bfet->factors.resize(len+1);
1073 for (int j = len; j > i; --j) {
1074 bfet->factors[j] = bfet->factors[j-1];
1075 t->terms[j] = t->terms[j-1];
1077 bfet->factors[i] = factor;
1078 factor = NULL;
1081 virtual void update_term(bfc_term_base *t, int i) {
1082 bfe_term* bfet = static_cast<bfe_term *>(t);
1084 eadd(factor, bfet->factors[i]);
1085 free_evalue_refs(factor);
1086 delete factor;
1089 virtual bool constant_vertex(int dim) { return E_num(0, dim) == 0; }
1091 virtual void cum(bf_reducer *bfr, bfc_term_base *t, int k, dpoly_r *r,
1092 barvinok_options *options);
1095 enumerator_base *enumerator_base::create(Polyhedron *P, unsigned dim,
1096 Param_Polyhedron *PP,
1097 barvinok_options *options)
1099 enumerator_base *eb;
1101 if (options->incremental_specialization == BV_SPECIALIZATION_BF)
1102 eb = new bfenumerator(dim, PP);
1103 else if (options->incremental_specialization == BV_SPECIALIZATION_DF)
1104 eb = new ienumerator(dim, PP);
1105 else
1106 eb = new enumerator(P, dim, PP);
1108 return eb;
1111 struct bfe_cum : public cumulator {
1112 bfenumerator *bfe;
1113 bfc_term_base *told;
1114 int k;
1115 bf_reducer *bfr;
1117 bfe_cum(evalue *factor, evalue *v, dpoly_r *r, bf_reducer *bfr,
1118 bfc_term_base *t, int k, bfenumerator *e) :
1119 cumulator(factor, v, r), bfe(e), told(t), k(k),
1120 bfr(bfr) {
1123 virtual void add_term(const vector<int>& powers, evalue *f2);
1126 void bfe_cum::add_term(const vector<int>& powers, evalue *f2)
1128 bfr->update_powers(powers);
1130 bfc_term_base * t = bfe->find_bfc_term(bfr->vn, bfr->npowers, bfr->nnf);
1131 bfe->set_factor(f2, bfr->l_changes % 2);
1132 bfe->add_term(t, told->terms[k], bfr->l_extra_num);
1135 void bfenumerator::cum(bf_reducer *bfr, bfc_term_base *t, int k,
1136 dpoly_r *r, barvinok_options *options)
1138 bfe_term* bfet = static_cast<bfe_term *>(t);
1139 bfe_cum cum(bfet->factors[k], E_num(0, bfr->d), r, bfr, t, k, this);
1140 cum.cumulate(options);
1143 void bfenumerator::base(mat_ZZ& factors, bfc_vec& v)
1145 for (int i = 0; i < v.size(); ++i) {
1146 assert(v[i]->terms.NumRows() == 1);
1147 evalue *factor = static_cast<bfe_term *>(v[i])->factors[0];
1148 eadd(factor, vE[vert]);
1149 delete v[i];
1153 void bfenumerator::handle(const signed_cone& sc, barvinok_options *options)
1155 assert(sc.det == 1);
1156 assert(sc.rays.NumRows() == enumerator_base::dim);
1158 bfe_term* t = new bfe_term(enumerator_base::dim);
1159 vector< bfc_term_base * > v;
1160 v.push_back(t);
1162 t->factors.resize(1);
1164 t->terms.SetDims(1, enumerator_base::dim);
1165 lattice_point(V, sc.rays, t->terms[0], E_vertex, options);
1167 // the elements of factors are always lexpositive
1168 mat_ZZ factors;
1169 int s = setup_factors(sc.rays, factors, t, sc.sign);
1171 t->factors[0] = new evalue;
1172 value_init(t->factors[0]->d);
1173 evalue_set_si(t->factors[0], s, 1);
1174 reduce(factors, v, options);
1176 for (int i = 0; i < enumerator_base::dim; ++i)
1177 if (E_vertex[i])
1178 evalue_free(E_vertex[i]);
1181 static evalue* barvinok_enumerate_ev_f(Polyhedron *P, Polyhedron* C,
1182 barvinok_options *options);
1184 /* Destroys C */
1185 static evalue* barvinok_enumerate_cst(Polyhedron *P, Polyhedron* C,
1186 struct barvinok_options *options)
1188 evalue *eres;
1190 if (emptyQ2(C)) {
1191 Polyhedron_Free(C);
1192 return evalue_zero();
1195 ALLOC(evalue, eres);
1196 value_init(eres->d);
1197 value_set_si(eres->d, 0);
1198 eres->x.p = new_enode(partition, 2, C->Dimension);
1199 EVALUE_SET_DOMAIN(eres->x.p->arr[0],
1200 DomainConstraintSimplify(C, options->MaxRays));
1201 value_set_si(eres->x.p->arr[1].d, 1);
1202 value_init(eres->x.p->arr[1].x.n);
1203 if (emptyQ2(P))
1204 value_set_si(eres->x.p->arr[1].x.n, 0);
1205 else
1206 barvinok_count_with_options(P, &eres->x.p->arr[1].x.n, options);
1207 if (value_mone_p(eres->x.p->arr[1].x.n)) {
1208 value_clear(eres->x.p->arr[1].x.n);
1209 value_set_si(eres->x.p->arr[1].d, -2); /* NaN */
1212 return eres;
1215 static evalue* enumerate(Polyhedron *P, Polyhedron* C,
1216 struct barvinok_options *options)
1218 Polyhedron *next;
1219 Polyhedron *Porig = P;
1220 Polyhedron *Corig = C;
1221 Polyhedron *CEq = NULL;
1222 unsigned nparam = C->Dimension;
1223 evalue *eres;
1224 Matrix *CP = NULL;
1226 evalue factor;
1227 value_init(factor.d);
1228 evalue_set_si(&factor, 1, 1);
1230 /* for now */
1231 POL_ENSURE_FACETS(P);
1232 POL_ENSURE_VERTICES(P);
1233 POL_ENSURE_FACETS(C);
1234 POL_ENSURE_VERTICES(C);
1236 if (C->Dimension == 0 || emptyQ(P) || emptyQ(C)) {
1237 constant:
1238 if (CEq == Porig)
1239 CEq = Polyhedron_Copy(CEq);
1240 eres = barvinok_enumerate_cst(P, CEq ? CEq : Polyhedron_Copy(C), options);
1241 out:
1242 if (CP) {
1243 evalue_backsubstitute(eres, CP, options->MaxRays);
1244 Matrix_Free(CP);
1247 emul(&factor, eres);
1248 if (options->approx->method == BV_APPROX_DROP) {
1249 if (options->approx->approximation == BV_APPROX_SIGN_UPPER)
1250 evalue_frac2polynomial(eres, 1, options->MaxRays);
1251 if (options->approx->approximation == BV_APPROX_SIGN_LOWER)
1252 evalue_frac2polynomial(eres, -1, options->MaxRays);
1253 if (options->approx->approximation == BV_APPROX_SIGN_APPROX)
1254 evalue_frac2polynomial(eres, 0, options->MaxRays);
1256 reduce_evalue(eres);
1257 free_evalue_refs(&factor);
1258 if (P != Porig)
1259 Domain_Free(P);
1260 if (C != Corig)
1261 Polyhedron_Free(C);
1263 return eres;
1265 if (Polyhedron_is_unbounded(P, nparam, options->MaxRays))
1266 goto constant;
1268 if (P->Dimension == nparam) {
1269 CEq = DomainIntersection(P, C, options->MaxRays);
1270 P = Universe_Polyhedron(0);
1271 goto constant;
1273 if (P->NbEq != 0 || C->NbEq != 0) {
1274 Polyhedron *Q = P;
1275 Polyhedron *D = C;
1276 remove_all_equalities(&P, &C, &CP, NULL, nparam, options->MaxRays);
1277 if (C != D && D != Corig)
1278 Polyhedron_Free(D);
1279 if (P != Q && Q != Porig)
1280 Domain_Free(Q);
1281 eres = enumerate(P, C, options);
1282 goto out;
1285 Polyhedron *T = Polyhedron_Factor(P, nparam, NULL, options->MaxRays);
1286 if (T || (P->Dimension == nparam+1)) {
1287 Polyhedron *C2 = C;
1288 Polyhedron *FC = Factor_Context(T ? T : P, nparam, options->MaxRays);
1289 C = DomainIntersection(C, FC, options->MaxRays);
1290 if (C2 != Corig)
1291 Polyhedron_Free(C2);
1292 Polyhedron_Free(FC);
1294 if (T) {
1295 if (P != Porig)
1296 Polyhedron_Free(P);
1297 P = T;
1298 if (T->Dimension == C->Dimension) {
1299 P = T->next;
1300 T->next = NULL;
1301 Polyhedron_Free(T);
1305 next = P->next;
1306 P->next = NULL;
1307 eres = barvinok_enumerate_ev_f(P, C, options);
1308 P->next = next;
1310 if (P->next) {
1311 Polyhedron *Q;
1312 evalue *f;
1314 for (Q = P->next; Q; Q = Q->next) {
1315 Polyhedron *next = Q->next;
1316 Q->next = NULL;
1318 f = barvinok_enumerate_ev_f(Q, C, options);
1319 emul(f, eres);
1320 evalue_free(f);
1322 Q->next = next;
1326 goto out;
1329 evalue* barvinok_enumerate_with_options(Polyhedron *P, Polyhedron* C,
1330 struct barvinok_options *options)
1332 Polyhedron *next, *Cnext, *C1;
1333 Polyhedron *Corig = C;
1334 evalue *eres;
1336 if (P->next)
1337 fprintf(stderr,
1338 "barvinok_enumerate: input is a union; only first polyhedron is enumerated\n");
1340 if (C->next)
1341 fprintf(stderr,
1342 "barvinok_enumerate: context is a union; only first polyhedron is considered\n");
1344 Cnext = C->next;
1345 C->next = NULL;
1346 C1 = Polyhedron_Project(P, C->Dimension);
1347 C = DomainIntersection(C, C1, options->MaxRays);
1348 Polyhedron_Free(C1);
1349 next = P->next;
1350 P->next = NULL;
1352 if (options->approx->method == BV_APPROX_BERNOULLI ||
1353 options->summation == BV_SUM_BERNOULLI) {
1354 int summation = options->summation;
1355 options->summation = BV_SUM_BERNOULLI;
1356 eres = barvinok_summate_unweighted(P, C, options);
1357 options->summation = summation;
1358 } else
1359 eres = enumerate(P, C, options);
1360 Domain_Free(C);
1362 P->next= next;
1363 Corig->next = Cnext;
1365 return eres;
1368 evalue* barvinok_enumerate_ev(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
1370 evalue *E;
1371 barvinok_options *options = barvinok_options_new_with_defaults();
1372 options->MaxRays = MaxRays;
1373 E = barvinok_enumerate_with_options(P, C, options);
1374 barvinok_options_free(options);
1375 return E;
1378 evalue *Param_Polyhedron_Enumerate(Param_Polyhedron *PP, Polyhedron *P,
1379 Polyhedron *C,
1380 struct barvinok_options *options)
1382 evalue *eres;
1383 Param_Domain *D;
1384 unsigned nparam = C->Dimension;
1385 unsigned dim = P->Dimension - nparam;
1387 int nd;
1388 for (nd = 0, D=PP->D; D; ++nd, D=D->next);
1389 evalue_section *s = new evalue_section[nd];
1390 Polyhedron *TC = true_context(P, C, options->MaxRays);
1392 enumerator_base *et = NULL;
1393 try_again:
1394 if (et)
1395 delete et;
1397 et = enumerator_base::create(P, dim, PP, options);
1399 FORALL_REDUCED_DOMAIN(PP, TC, nd, options, i, D, rVD)
1400 Param_Vertices *V;
1402 s[i].E = evalue_zero();
1403 s[i].D = rVD;
1405 FORALL_PVertex_in_ParamPolyhedron(V,D,PP) // _i is internal counter
1406 if (!et->vE[_i])
1407 try {
1408 et->decompose_at(V, _i, options);
1409 } catch (OrthogonalException &e) {
1410 FORALL_REDUCED_DOMAIN_RESET;
1411 for (; i >= 0; --i) {
1412 evalue_free(s[i].E);
1413 Domain_Free(s[i].D);
1415 goto try_again;
1417 eadd(et->vE[_i] , s[i].E);
1418 END_FORALL_PVertex_in_ParamPolyhedron;
1419 evalue_range_reduction_in_domain(s[i].E, rVD);
1420 END_FORALL_REDUCED_DOMAIN
1421 Polyhedron_Free(TC);
1423 delete et;
1424 eres = evalue_from_section_array(s, nd);
1425 delete [] s;
1427 return eres;
1430 static evalue* barvinok_enumerate_ev_f(Polyhedron *P, Polyhedron* C,
1431 barvinok_options *options)
1433 unsigned nparam = C->Dimension;
1434 bool do_scale = options->approx->method == BV_APPROX_SCALE;
1436 if (options->summation == BV_SUM_EULER)
1437 return barvinok_summate_unweighted(P, C, options);
1439 if (options->approx->method == BV_APPROX_VOLUME)
1440 return Param_Polyhedron_Volume(P, C, options);
1442 if (P->Dimension - nparam == 1 && !do_scale)
1443 return ParamLine_Length(P, C, options);
1445 Param_Polyhedron *PP = NULL;
1446 evalue *eres;
1448 if (do_scale) {
1449 eres = scale_bound(P, C, options);
1450 if (eres)
1451 return eres;
1454 PP = Polyhedron2Param_Polyhedron(P, C, options);
1456 if (do_scale)
1457 eres = scale(PP, P, C, options);
1458 else
1459 eres = Param_Polyhedron_Enumerate(PP, P, C, options);
1461 if (PP)
1462 Param_Polyhedron_Free(PP);
1464 return eres;
1467 Enumeration* barvinok_enumerate(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
1469 evalue *EP = barvinok_enumerate_ev(P, C, MaxRays);
1471 return partition2enumeration(EP);
1474 evalue* barvinok_enumerate_union(Polyhedron *D, Polyhedron* C, unsigned MaxRays)
1476 evalue *EP;
1477 gen_fun *gf = barvinok_enumerate_union_series(D, C, MaxRays);
1478 EP = *gf;
1479 delete gf;
1480 return EP;
1483 static __isl_give isl_pw_qpolynomial *basic_set_card(
1484 __isl_take isl_basic_set *bset)
1486 isl_ctx *ctx;
1487 isl_space *dim;
1488 isl_pw_qpolynomial *pwqp;
1489 unsigned nparam = isl_basic_set_dim(bset, isl_dim_param);
1490 Polyhedron *U = Universe_Polyhedron(nparam);
1491 Polyhedron *P;
1492 evalue *E;
1493 barvinok_options *options;
1494 int options_allocated = 0;
1496 ctx = isl_basic_set_get_ctx(bset);
1497 options = isl_ctx_peek_barvinok_options(ctx);
1498 if (!options) {
1499 options = barvinok_options_new_with_defaults();
1500 options_allocated = 1;
1503 dim = isl_basic_set_get_space(bset);
1504 dim = isl_space_domain(dim);
1506 P = isl_basic_set_to_polylib(bset);
1507 E = enumerate(P, U, options);
1509 pwqp = isl_pw_qpolynomial_from_evalue(dim, E);
1510 isl_basic_set_free(bset);
1512 evalue_free(E);
1513 Polyhedron_Free(P);
1514 Polyhedron_Free(U);
1515 if (options_allocated)
1516 barvinok_options_free(options);
1518 return pwqp;
1521 static isl_stat basic_map_card(__isl_take isl_basic_map *bmap, void *user)
1523 isl_pw_qpolynomial **sum = (isl_pw_qpolynomial **)user;
1524 isl_pw_qpolynomial *pwqp;
1525 unsigned nparam = isl_basic_map_dim(bmap, isl_dim_param);
1526 unsigned n_in = isl_basic_map_dim(bmap, isl_dim_in);
1527 isl_space *target_dim;
1528 isl_basic_set *bset;
1530 target_dim = isl_basic_map_get_space(bmap);
1531 target_dim = isl_space_domain(target_dim);
1533 bmap = isl_basic_map_move_dims(bmap, isl_dim_param, nparam,
1534 isl_dim_in, 0, n_in);
1536 bset = isl_basic_map_range(bmap);
1537 bset = isl_basic_set_lift(bset);
1538 pwqp = isl_basic_set_multiplicative_call(bset, &basic_set_card);
1540 pwqp = isl_pw_qpolynomial_move_dims(pwqp, isl_dim_in, 0,
1541 isl_dim_param, nparam, n_in);
1542 pwqp = isl_pw_qpolynomial_reset_domain_space(pwqp, target_dim);
1543 *sum = isl_pw_qpolynomial_add(*sum, pwqp);
1545 return isl_stat_ok;
1548 static __isl_give isl_pw_qpolynomial *card_as_sum(__isl_take isl_map *map,
1549 barvinok_options *options)
1551 isl_ctx *ctx;
1552 isl_val *one;
1553 isl_space *dim;
1554 isl_set *set;
1555 isl_qpolynomial *qp;
1556 isl_pw_qpolynomial *pwqp;
1557 int summation = options->summation;
1559 if (!map)
1560 return NULL;
1562 options->summation = BV_SUM_BERNOULLI;
1564 set = isl_map_wrap(map);
1565 dim = isl_set_get_space(set);
1566 ctx = isl_map_get_ctx(map);
1567 one = isl_val_one(ctx);
1568 qp = isl_qpolynomial_val_on_domain(dim, one);
1570 pwqp = isl_pw_qpolynomial_alloc(set, qp);
1571 pwqp = isl_pw_qpolynomial_sum(pwqp);
1573 options->summation = summation;
1575 return pwqp;
1578 __isl_give isl_pw_qpolynomial *isl_map_card(__isl_take isl_map *map)
1580 isl_ctx *ctx;
1581 isl_space *dim;
1582 isl_pw_qpolynomial *sum;
1583 barvinok_options *options;
1585 ctx = isl_map_get_ctx(map);
1586 options = isl_ctx_peek_barvinok_options(ctx);
1587 if (options &&
1588 (options->approx->method == BV_APPROX_BERNOULLI ||
1589 options->summation == BV_SUM_BERNOULLI))
1590 return card_as_sum(map, options);
1592 dim = isl_map_get_space(map);
1593 dim = isl_space_domain(dim);
1594 dim = isl_space_from_domain(dim);
1595 dim = isl_space_add_dims(dim, isl_dim_out, 1);
1596 sum = isl_pw_qpolynomial_zero(dim);
1598 map = isl_map_make_disjoint(map);
1599 map = isl_map_compute_divs(map);
1601 if (isl_map_foreach_basic_map(map, &basic_map_card, &sum) < 0)
1602 goto error;
1604 isl_map_free(map);
1606 return sum;
1607 error:
1608 isl_map_free(map);
1609 isl_pw_qpolynomial_free(sum);
1610 return NULL;
1613 __isl_give isl_pw_qpolynomial *isl_set_card(__isl_take isl_set *set)
1615 isl_pw_qpolynomial *pwqp;
1616 pwqp = isl_map_card(isl_map_from_range(set));
1617 pwqp = isl_pw_qpolynomial_project_domain_on_params(pwqp);
1618 return pwqp;
1621 __isl_give isl_pw_qpolynomial *isl_basic_map_card(__isl_take isl_basic_map *bmap)
1623 return isl_map_card(isl_map_from_basic_map(bmap));
1626 __isl_give isl_pw_qpolynomial *isl_basic_set_card(__isl_take isl_basic_set *bset)
1628 isl_pw_qpolynomial *pwqp;
1629 pwqp = isl_basic_map_card(isl_basic_map_from_range(bset));
1630 pwqp = isl_pw_qpolynomial_project_domain_on_params(pwqp);
1631 return pwqp;
1634 static isl_stat set_card(__isl_take isl_set *set, void *user)
1636 isl_union_pw_qpolynomial **res = (isl_union_pw_qpolynomial **)user;
1637 isl_pw_qpolynomial *pwqp;
1638 isl_union_pw_qpolynomial *upwqp;
1640 pwqp = isl_set_card(set);
1641 upwqp = isl_union_pw_qpolynomial_from_pw_qpolynomial(pwqp);
1642 *res = isl_union_pw_qpolynomial_add(*res, upwqp);
1644 return isl_stat_ok;
1647 __isl_give isl_union_pw_qpolynomial *isl_union_set_card(
1648 __isl_take isl_union_set *uset)
1650 isl_space *dim;
1651 isl_union_pw_qpolynomial *res;
1653 dim = isl_union_set_get_space(uset);
1654 res = isl_union_pw_qpolynomial_zero(dim);
1655 if (isl_union_set_foreach_set(uset, &set_card, &res) < 0)
1656 goto error;
1657 isl_union_set_free(uset);
1659 return res;
1660 error:
1661 isl_union_set_free(uset);
1662 isl_union_pw_qpolynomial_free(res);
1663 return NULL;
1666 static isl_stat map_card(__isl_take isl_map *map, void *user)
1668 isl_union_pw_qpolynomial **res = (isl_union_pw_qpolynomial **)user;
1669 isl_pw_qpolynomial *pwqp;
1670 isl_union_pw_qpolynomial *upwqp;
1672 pwqp = isl_map_card(map);
1673 upwqp = isl_union_pw_qpolynomial_from_pw_qpolynomial(pwqp);
1674 *res = isl_union_pw_qpolynomial_add(*res, upwqp);
1676 return isl_stat_ok;
1679 __isl_give isl_union_pw_qpolynomial *isl_union_map_card(
1680 __isl_take isl_union_map *umap)
1682 isl_space *dim;
1683 isl_union_pw_qpolynomial *res;
1685 dim = isl_union_map_get_space(umap);
1686 res = isl_union_pw_qpolynomial_zero(dim);
1687 if (isl_union_map_foreach_map(umap, &map_card, &res) < 0)
1688 goto error;
1689 isl_union_map_free(umap);
1691 return res;
1692 error:
1693 isl_union_map_free(umap);
1694 isl_union_pw_qpolynomial_free(res);
1695 return NULL;