README: explain how to get submodules
[barvinok.git] / barvinok.cc
blob573b47670f3a62bd4459e037baa776362af5d0e8
1 #include <assert.h>
2 #include <iostream>
3 #include <vector>
4 #include <deque>
5 #include <string>
6 #include <sstream>
7 #include <gmp.h>
8 #include <NTL/mat_ZZ.h>
9 #include <NTL/LLL.h>
10 #include <barvinok/util.h>
11 #include <barvinok/evalue.h>
12 #include "config.h"
13 #include <barvinok/barvinok.h>
14 #include <barvinok/genfun.h>
15 #include <barvinok/options.h>
16 #include <barvinok/sample.h>
17 #include "bfcounter.h"
18 #include "conversion.h"
19 #include "counter.h"
20 #include "decomposer.h"
21 #include "euler.h"
22 #include "lattice_point.h"
23 #include "reduce_domain.h"
24 #include "remove_equalities.h"
25 #include "scale.h"
26 #include "volume.h"
27 #include "bernoulli.h"
28 #include "param_util.h"
30 #ifdef NTL_STD_CXX
31 using namespace NTL;
32 #endif
33 using std::cerr;
34 using std::cout;
35 using std::endl;
36 using std::vector;
37 using std::deque;
38 using std::string;
39 using std::ostringstream;
41 #define ALLOC(t,p) p = (t*)malloc(sizeof(*p))
43 class dpoly_n {
44 public:
45 Matrix *coeff;
46 ~dpoly_n() {
47 Matrix_Free(coeff);
49 dpoly_n(int d) {
50 Value d0, one;
51 value_init(d0);
52 value_init(one);
53 value_set_si(one, 1);
54 coeff = Matrix_Alloc(d+1, d+1+1);
55 value_set_si(coeff->p[0][0], 1);
56 value_set_si(coeff->p[0][d+1], 1);
57 for (int i = 1; i <= d; ++i) {
58 value_multiply(coeff->p[i][0], coeff->p[i-1][0], d0);
59 Vector_Combine(coeff->p[i-1], coeff->p[i-1]+1, coeff->p[i]+1,
60 one, d0, i);
61 value_set_si(coeff->p[i][d+1], i);
62 value_multiply(coeff->p[i][d+1], coeff->p[i][d+1], coeff->p[i-1][d+1]);
63 value_decrement(d0, d0);
65 value_clear(d0);
66 value_clear(one);
68 void div(dpoly& d, Vector *count, int sign) {
69 int len = coeff->NbRows;
70 Matrix * c = Matrix_Alloc(coeff->NbRows, coeff->NbColumns);
71 Value tmp;
72 value_init(tmp);
73 for (int i = 0; i < len; ++i) {
74 Vector_Copy(coeff->p[i], c->p[i], len+1);
75 for (int j = 1; j <= i; ++j) {
76 value_multiply(tmp, d.coeff->p[j], c->p[i][len]);
77 value_oppose(tmp, tmp);
78 Vector_Combine(c->p[i], c->p[i-j], c->p[i],
79 c->p[i-j][len], tmp, len);
80 value_multiply(c->p[i][len], c->p[i][len], c->p[i-j][len]);
82 value_multiply(c->p[i][len], c->p[i][len], d.coeff->p[0]);
84 if (sign == -1) {
85 value_set_si(tmp, -1);
86 Vector_Scale(c->p[len-1], count->p, tmp, len);
87 value_assign(count->p[len], c->p[len-1][len]);
88 } else
89 Vector_Copy(c->p[len-1], count->p, len+1);
90 Vector_Normalize(count->p, len+1);
91 value_clear(tmp);
92 Matrix_Free(c);
96 const int MAX_TRY=10;
98 * Searches for a vector that is not orthogonal to any
99 * of the rays in rays.
101 static void nonorthog(mat_ZZ& rays, vec_ZZ& lambda)
103 int dim = rays.NumCols();
104 bool found = false;
105 lambda.SetLength(dim);
106 if (dim == 0)
107 return;
109 for (int i = 2; !found && i <= 50*dim; i+=4) {
110 for (int j = 0; j < MAX_TRY; ++j) {
111 for (int k = 0; k < dim; ++k) {
112 int r = random_int(i)+2;
113 int v = (2*(r%2)-1) * (r >> 1);
114 lambda[k] = v;
116 int k = 0;
117 for (; k < rays.NumRows(); ++k)
118 if (lambda * rays[k] == 0)
119 break;
120 if (k == rays.NumRows()) {
121 found = true;
122 break;
126 assert(found);
129 static void add_rays(mat_ZZ& rays, Polyhedron *i, int *r, int nvar = -1,
130 bool all = false)
132 unsigned dim = i->Dimension;
133 if (nvar == -1)
134 nvar = dim;
135 for (int k = 0; k < i->NbRays; ++k) {
136 if (!value_zero_p(i->Ray[k][dim+1]))
137 continue;
138 if (!all && nvar != dim && First_Non_Zero(i->Ray[k]+1, nvar) == -1)
139 continue;
140 values2zz(i->Ray[k]+1, rays[(*r)++], nvar);
144 static void mask_r(Matrix *f, int nr, Vector *lcm, int p, Vector *val, evalue *ev)
146 unsigned nparam = lcm->Size;
148 if (p == nparam) {
149 Vector * prod = Vector_Alloc(f->NbRows);
150 Matrix_Vector_Product(f, val->p, prod->p);
151 int isint = 1;
152 for (int i = 0; i < nr; ++i) {
153 value_modulus(prod->p[i], prod->p[i], f->p[i][nparam+1]);
154 isint &= value_zero_p(prod->p[i]);
156 value_set_si(ev->d, 1);
157 value_init(ev->x.n);
158 value_set_si(ev->x.n, isint);
159 Vector_Free(prod);
160 return;
163 Value tmp;
164 value_init(tmp);
165 if (value_one_p(lcm->p[p]))
166 mask_r(f, nr, lcm, p+1, val, ev);
167 else {
168 value_assign(tmp, lcm->p[p]);
169 value_set_si(ev->d, 0);
170 ev->x.p = new_enode(periodic, VALUE_TO_INT(tmp), p+1);
171 do {
172 value_decrement(tmp, tmp);
173 value_assign(val->p[p], tmp);
174 mask_r(f, nr, lcm, p+1, val, &ev->x.p->arr[VALUE_TO_INT(tmp)]);
175 } while (value_pos_p(tmp));
177 value_clear(tmp);
180 static void mask_fractional(Matrix *f, evalue *factor)
182 int nr = f->NbRows, nc = f->NbColumns;
183 int n;
184 bool found = false;
185 for (n = 0; n < nr && value_notzero_p(f->p[n][nc-1]); ++n)
186 if (value_notone_p(f->p[n][nc-1]) &&
187 value_notmone_p(f->p[n][nc-1]))
188 found = true;
189 if (!found)
190 return;
192 evalue EP;
193 nr = n;
195 Value m;
196 value_init(m);
198 evalue EV;
199 value_init(EV.d);
200 value_init(EV.x.n);
201 value_set_si(EV.x.n, 1);
203 for (n = 0; n < nr; ++n) {
204 value_assign(m, f->p[n][nc-1]);
205 if (value_one_p(m) || value_mone_p(m))
206 continue;
208 int j = normal_mod(f->p[n], nc-1, &m);
209 if (j == nc-1) {
210 free_evalue_refs(factor);
211 value_init(factor->d);
212 evalue_set_si(factor, 0, 1);
213 break;
215 vec_ZZ row;
216 values2zz(f->p[n], row, nc-1);
217 ZZ g;
218 value2zz(m, g);
219 if (j < (nc-1)-1 && row[j] > g/2) {
220 for (int k = j; k < (nc-1); ++k)
221 if (row[k] != 0)
222 row[k] = g - row[k];
225 value_init(EP.d);
226 value_set_si(EP.d, 0);
227 EP.x.p = new_enode(relation, 2, 0);
228 value_clear(EP.x.p->arr[1].d);
229 EP.x.p->arr[1] = *factor;
230 evalue *ev = &EP.x.p->arr[0];
231 value_set_si(ev->d, 0);
232 ev->x.p = new_enode(fractional, 3, -1);
233 evalue_set_si(&ev->x.p->arr[1], 0, 1);
234 evalue_set_si(&ev->x.p->arr[2], 1, 1);
235 evalue *E = multi_monom(row);
236 value_assign(EV.d, m);
237 emul(&EV, E);
238 value_clear(ev->x.p->arr[0].d);
239 ev->x.p->arr[0] = *E;
240 delete E;
241 *factor = EP;
244 value_clear(m);
245 free_evalue_refs(&EV);
251 static void mask_table(Matrix *f, evalue *factor)
253 int nr = f->NbRows, nc = f->NbColumns;
254 int n;
255 bool found = false;
256 for (n = 0; n < nr && value_notzero_p(f->p[n][nc-1]); ++n)
257 if (value_notone_p(f->p[n][nc-1]) &&
258 value_notmone_p(f->p[n][nc-1]))
259 found = true;
260 if (!found)
261 return;
263 Value tmp;
264 value_init(tmp);
265 nr = n;
266 unsigned np = nc - 2;
267 Vector *lcm = Vector_Alloc(np);
268 Vector *val = Vector_Alloc(nc);
269 Vector_Set(val->p, 0, nc);
270 value_set_si(val->p[np], 1);
271 Vector_Set(lcm->p, 1, np);
272 for (n = 0; n < nr; ++n) {
273 if (value_one_p(f->p[n][nc-1]) ||
274 value_mone_p(f->p[n][nc-1]))
275 continue;
276 for (int j = 0; j < np; ++j)
277 if (value_notzero_p(f->p[n][j])) {
278 value_gcd(tmp, f->p[n][j], f->p[n][nc-1]);
279 value_division(tmp, f->p[n][nc-1], tmp);
280 value_lcm(lcm->p[j], tmp, lcm->p[j]);
283 evalue EP;
284 value_init(EP.d);
285 mask_r(f, nr, lcm, 0, val, &EP);
286 value_clear(tmp);
287 Vector_Free(val);
288 Vector_Free(lcm);
289 emul(&EP,factor);
290 free_evalue_refs(&EP);
293 static void mask(Matrix *f, evalue *factor, barvinok_options *options)
295 if (options->lookup_table)
296 mask_table(f, factor);
297 else
298 mask_fractional(f, factor);
301 struct bfe_term : public bfc_term_base {
302 vector<evalue *> factors;
304 bfe_term(int len) : bfc_term_base(len) {
307 ~bfe_term() {
308 for (int i = 0; i < factors.size(); ++i) {
309 if (!factors[i])
310 continue;
311 free_evalue_refs(factors[i]);
312 delete factors[i];
317 static void print_int_vector(int *v, int len, const char *name)
319 cerr << name << endl;
320 for (int j = 0; j < len; ++j) {
321 cerr << v[j] << " ";
323 cerr << endl;
326 static void print_bfc_terms(mat_ZZ& factors, bfc_vec& v)
328 cerr << endl;
329 cerr << "factors" << endl;
330 cerr << factors << endl;
331 for (int i = 0; i < v.size(); ++i) {
332 cerr << "term: " << i << endl;
333 print_int_vector(v[i]->powers, factors.NumRows(), "powers");
334 cerr << "terms" << endl;
335 cerr << v[i]->terms << endl;
336 bfc_term* bfct = static_cast<bfc_term *>(v[i]);
337 cerr << bfct->c << endl;
341 static void print_bfe_terms(mat_ZZ& factors, bfc_vec& v)
343 cerr << endl;
344 cerr << "factors" << endl;
345 cerr << factors << endl;
346 for (int i = 0; i < v.size(); ++i) {
347 cerr << "term: " << i << endl;
348 print_int_vector(v[i]->powers, factors.NumRows(), "powers");
349 cerr << "terms" << endl;
350 cerr << v[i]->terms << endl;
351 bfe_term* bfet = static_cast<bfe_term *>(v[i]);
352 for (int j = 0; j < v[i]->terms.NumRows(); ++j) {
353 const char * test[] = {"a", "b"};
354 print_evalue(stderr, bfet->factors[j], test);
355 fprintf(stderr, "\n");
360 struct bfcounter : public bfcounter_base {
361 mpq_t count;
362 Value tz;
364 bfcounter(unsigned dim) : bfcounter_base(dim) {
365 mpq_init(count);
366 lower = 1;
367 value_init(tz);
369 ~bfcounter() {
370 mpq_clear(count);
371 value_clear(tz);
373 virtual void base(mat_ZZ& factors, bfc_vec& v);
374 virtual void get_count(Value *result) {
375 assert(value_one_p(&count[0]._mp_den));
376 value_assign(*result, &count[0]._mp_num);
380 void bfcounter::base(mat_ZZ& factors, bfc_vec& v)
382 unsigned nf = factors.NumRows();
384 for (int i = 0; i < v.size(); ++i) {
385 bfc_term* bfct = static_cast<bfc_term *>(v[i]);
386 int total_power = 0;
387 // factor is always positive, so we always
388 // change signs
389 for (int k = 0; k < nf; ++k)
390 total_power += v[i]->powers[k];
392 int j;
393 for (j = 0; j < nf; ++j)
394 if (v[i]->powers[j] > 0)
395 break;
397 zz2value(factors[j][0], tz);
398 dpoly D(total_power, tz, 1);
399 for (int k = 1; k < v[i]->powers[j]; ++k) {
400 zz2value(factors[j][0], tz);
401 dpoly fact(total_power, tz, 1);
402 D *= fact;
404 for ( ; ++j < nf; )
405 for (int k = 0; k < v[i]->powers[j]; ++k) {
406 zz2value(factors[j][0], tz);
407 dpoly fact(total_power, tz, 1);
408 D *= fact;
411 for (int k = 0; k < v[i]->terms.NumRows(); ++k) {
412 zz2value(v[i]->terms[k][0], tz);
413 dpoly n(total_power, tz);
414 mpq_set_si(tcount, 0, 1);
415 n.div(D, tcount, 1);
416 if (total_power % 2)
417 bfct->c[k].n = -bfct->c[k].n;
418 zz2value(bfct->c[k].n, tn);
419 zz2value(bfct->c[k].d, td);
421 mpz_mul(mpq_numref(tcount), mpq_numref(tcount), tn);
422 mpz_mul(mpq_denref(tcount), mpq_denref(tcount), td);
423 mpq_canonicalize(tcount);
424 mpq_add(count, count, tcount);
426 delete v[i];
431 /* Check whether the polyhedron is unbounded and if so,
432 * check whether it has any (and therefore an infinite number of)
433 * integer points.
434 * If one of the vertices is integer, then we are done.
435 * Otherwise, transform the polyhedron such that one of the rays
436 * is the first unit vector and cut it off at a height that ensures
437 * that if the whole polyhedron has any points, then the remaining part
438 * has integer points. In particular we add the largest coefficient
439 * of a ray to the highest vertex (rounded up).
441 static bool Polyhedron_is_infinite(Polyhedron *P, Value* result,
442 barvinok_options *options)
444 int r = 0;
445 Matrix *M, *M2;
446 Value c, tmp;
447 Value g;
448 bool first;
449 Vector *v;
450 Value offset, size;
451 Polyhedron *R;
453 if (P->NbBid == 0)
454 for (; r < P->NbRays; ++r)
455 if (value_zero_p(P->Ray[r][P->Dimension+1]))
456 break;
457 if (P->NbBid == 0 && r == P->NbRays)
458 return false;
460 if (options->count_sample_infinite) {
461 Vector *sample;
463 sample = Polyhedron_Sample(P, options);
464 if (!sample)
465 value_set_si(*result, 0);
466 else {
467 value_set_si(*result, -1);
468 Vector_Free(sample);
470 return true;
473 for (int i = 0; i < P->NbRays; ++i)
474 if (value_one_p(P->Ray[i][1+P->Dimension])) {
475 value_set_si(*result, -1);
476 return true;
479 value_init(g);
480 M = Matrix_Alloc(P->Dimension+1, P->Dimension+1);
481 Vector_Gcd(P->Ray[r]+1, P->Dimension, &g);
482 Vector_AntiScale(P->Ray[r]+1, M->p[0], g, P->Dimension+1);
483 int ok = unimodular_complete(M, 1);
484 assert(ok);
485 value_set_si(M->p[P->Dimension][P->Dimension], 1);
486 M2 = Transpose(M);
487 Matrix_Free(M);
488 P = Polyhedron_Preimage(P, M2, 0);
489 Matrix_Free(M2);
490 value_clear(g);
492 first = true;
493 value_init(offset);
494 value_init(size);
495 value_init(tmp);
496 value_set_si(size, 0);
498 for (int i = 0; i < P->NbBid; ++i) {
499 value_absolute(tmp, P->Ray[i][1]);
500 if (value_gt(tmp, size))
501 value_assign(size, tmp);
503 for (int i = P->NbBid; i < P->NbRays; ++i) {
504 if (value_zero_p(P->Ray[i][P->Dimension+1])) {
505 if (value_gt(P->Ray[i][1], size))
506 value_assign(size, P->Ray[i][1]);
507 continue;
509 mpz_cdiv_q(tmp, P->Ray[i][1], P->Ray[i][P->Dimension+1]);
510 if (first || value_gt(tmp, offset)) {
511 value_assign(offset, tmp);
512 first = false;
515 value_addto(offset, offset, size);
516 value_clear(size);
517 value_clear(tmp);
519 v = Vector_Alloc(P->Dimension+2);
520 value_set_si(v->p[0], 1);
521 value_set_si(v->p[1], -1);
522 value_assign(v->p[1+P->Dimension], offset);
523 R = AddConstraints(v->p, 1, P, options->MaxRays);
524 Polyhedron_Free(P);
525 P = R;
527 value_clear(offset);
528 Vector_Free(v);
530 value_init(c);
531 barvinok_count_with_options(P, &c, options);
532 Polyhedron_Free(P);
533 if (value_zero_p(c))
534 value_set_si(*result, 0);
535 else
536 value_set_si(*result, -1);
537 value_clear(c);
539 return true;
542 static void barvinok_count_f(Polyhedron *P, Value* result,
543 barvinok_options *options);
545 void barvinok_count_with_options(Polyhedron *P, Value* result,
546 struct barvinok_options *options)
548 unsigned dim;
549 int allocated = 0;
550 Polyhedron *Q;
551 bool infinite = false;
553 if (P->next)
554 fprintf(stderr,
555 "barvinok_count: input is a union; only first polyhedron is counted\n");
557 if (emptyQ2(P)) {
558 value_set_si(*result, 0);
559 return;
561 if (P->NbEq != 0) {
562 Q = NULL;
563 do {
564 P = remove_equalities(P, options->MaxRays);
565 P = DomainConstraintSimplify(P, options->MaxRays);
566 if (Q)
567 Polyhedron_Free(Q);
568 Q = P;
569 } while (!emptyQ(P) && P->NbEq != 0);
570 if (emptyQ(P)) {
571 Polyhedron_Free(P);
572 value_set_si(*result, 0);
573 return;
575 allocated = 1;
577 if (Polyhedron_is_infinite(P, result, options)) {
578 if (allocated)
579 Polyhedron_Free(P);
580 return;
582 if (P->Dimension == 0) {
583 /* Test whether the constraints are satisfied */
584 POL_ENSURE_VERTICES(P);
585 value_set_si(*result, !emptyQ(P));
586 if (allocated)
587 Polyhedron_Free(P);
588 return;
590 Q = Polyhedron_Factor(P, 0, NULL, options->MaxRays);
591 if (Q) {
592 if (allocated)
593 Polyhedron_Free(P);
594 P = Q;
595 allocated = 1;
598 barvinok_count_f(P, result, options);
599 if (value_neg_p(*result))
600 infinite = true;
601 if (Q && P->next && value_notzero_p(*result)) {
602 Value factor;
603 value_init(factor);
605 for (Q = P->next; Q; Q = Q->next) {
606 barvinok_count_f(Q, &factor, options);
607 if (value_neg_p(factor)) {
608 infinite = true;
609 continue;
610 } else if (Q->next && value_zero_p(factor)) {
611 value_set_si(*result, 0);
612 break;
614 value_multiply(*result, *result, factor);
617 value_clear(factor);
620 if (allocated)
621 Domain_Free(P);
622 if (infinite)
623 value_set_si(*result, -1);
626 void barvinok_count(Polyhedron *P, Value* result, unsigned NbMaxCons)
628 barvinok_options *options = barvinok_options_new_with_defaults();
629 options->MaxRays = NbMaxCons;
630 barvinok_count_with_options(P, result, options);
631 barvinok_options_free(options);
634 static void barvinok_count_f(Polyhedron *P, Value* result,
635 barvinok_options *options)
637 if (emptyQ2(P)) {
638 value_set_si(*result, 0);
639 return;
642 if (P->Dimension == 1)
643 return Line_Length(P, result);
645 int c = P->NbConstraints;
646 POL_ENSURE_FACETS(P);
647 if (c != P->NbConstraints || P->NbEq != 0) {
648 Polyhedron *next = P->next;
649 P->next = NULL;
650 barvinok_count_with_options(P, result, options);
651 P->next = next;
652 return;
655 POL_ENSURE_VERTICES(P);
657 if (Polyhedron_is_infinite(P, result, options))
658 return;
660 np_base *cnt;
661 if (options->incremental_specialization == BV_SPECIALIZATION_BF)
662 cnt = new bfcounter(P->Dimension);
663 else if (options->incremental_specialization == BV_SPECIALIZATION_DF)
664 cnt = new icounter(P->Dimension);
665 else if (options->incremental_specialization == BV_SPECIALIZATION_TODD)
666 cnt = new tcounter(P->Dimension, options->max_index);
667 else
668 cnt = new counter(P->Dimension, options->max_index);
669 cnt->start(P, options);
671 cnt->get_count(result);
672 delete cnt;
675 static void uni_polynom(int param, Vector *c, evalue *EP)
677 unsigned dim = c->Size-2;
678 value_init(EP->d);
679 value_set_si(EP->d,0);
680 EP->x.p = new_enode(polynomial, dim+1, param+1);
681 for (int j = 0; j <= dim; ++j)
682 evalue_set(&EP->x.p->arr[j], c->p[j], c->p[dim+1]);
685 typedef evalue * evalue_p;
687 struct enumerator_base {
688 unsigned dim;
689 evalue ** vE;
690 evalue mone;
691 vertex_decomposer *vpd;
693 enumerator_base(unsigned dim, vertex_decomposer *vpd)
695 this->dim = dim;
696 this->vpd = vpd;
698 vE = new evalue_p[vpd->PP->nbV];
699 for (int j = 0; j < vpd->PP->nbV; ++j)
700 vE[j] = 0;
702 value_init(mone.d);
703 evalue_set_si(&mone, -1, 1);
706 void decompose_at(Param_Vertices *V, int _i, barvinok_options *options) {
707 //this->pVD = pVD;
709 vE[_i] = new evalue;
710 value_init(vE[_i]->d);
711 evalue_set_si(vE[_i], 0, 1);
713 vpd->decompose_at_vertex(V, _i, options);
716 virtual ~enumerator_base() {
717 for (int j = 0; j < vpd->PP->nbV; ++j)
718 if (vE[j]) {
719 free_evalue_refs(vE[j]);
720 delete vE[j];
722 delete [] vE;
724 free_evalue_refs(&mone);
727 static enumerator_base *create(Polyhedron *P, unsigned dim,
728 Param_Polyhedron *PP,
729 barvinok_options *options);
732 struct enumerator : public signed_cone_consumer, public vertex_decomposer,
733 public enumerator_base {
734 vec_ZZ lambda;
735 vec_ZZ den;
736 term_info num;
737 Vector *c;
738 mpq_t count;
739 Value tz;
741 enumerator(Polyhedron *P, unsigned dim, Param_Polyhedron *PP) :
742 vertex_decomposer(PP, *this), enumerator_base(dim, this) {
743 randomvector(P, lambda, dim);
744 den.SetLength(dim);
745 c = Vector_Alloc(dim+2);
747 mpq_init(count);
748 value_init(tz);
751 ~enumerator() {
752 mpq_clear(count);
753 Vector_Free(c);
754 value_clear(tz);
757 virtual void handle(const signed_cone& sc, barvinok_options *options);
760 void enumerator::handle(const signed_cone& sc, barvinok_options *options)
762 int sign = sc.sign;
763 int r = 0;
764 assert(sc.rays.NumRows() == dim);
765 for (int k = 0; k < dim; ++k) {
766 if (lambda * sc.rays[k] == 0)
767 throw Orthogonal;
770 lattice_point(V, sc.rays, lambda, &num, sc.det, options);
771 den = sc.rays * lambda;
773 if (dim % 2)
774 sign = -sign;
776 zz2value(den[0], tz);
777 dpoly n(dim, tz, 1);
778 for (int k = 1; k < dim; ++k) {
779 zz2value(den[k], tz);
780 dpoly fact(dim, tz, 1);
781 n *= fact;
783 if (num.E != NULL) {
784 dpoly_n d(dim);
785 d.div(n, c, sign);
786 for (unsigned long i = 0; i < sc.det; ++i) {
787 evalue *EV = evalue_polynomial(c, num.E[i]);
788 eadd(EV, vE[vert]);
789 evalue_free(EV);
790 free_evalue_refs(num.E[i]);
791 delete num.E[i];
793 delete [] num.E;
794 } else {
795 mpq_set_si(count, 0, 1);
796 if (num.constant.length() == 1) {
797 zz2value(num.constant[0], tz);
798 dpoly d(dim, tz);
799 d.div(n, count, sign);
800 } else {
801 dpoly_n d(dim);
802 d.div(n, c, sign);
803 Value x, sum, acc;
804 value_init(x);
805 value_init(acc);
806 for (unsigned long i = 0; i < sc.det; ++i) {
807 value_assign(acc, c->p[dim]);
808 zz2value(num.constant[i], x);
809 for (int j = dim-1; j >= 0; --j) {
810 value_multiply(acc, acc, x);
811 value_addto(acc, acc, c->p[j]);
813 value_addto(mpq_numref(count), mpq_numref(count), acc);
815 mpz_set(mpq_denref(count), c->p[dim+1]);
816 value_clear(acc);
817 value_clear(x);
819 evalue EV;
820 value_init(EV.d);
821 evalue_set(&EV, &count[0]._mp_num, &count[0]._mp_den);
822 eadd(&EV, vE[vert]);
823 free_evalue_refs(&EV);
827 struct ienumerator_base : enumerator_base {
828 evalue ** E_vertex;
830 ienumerator_base(unsigned dim, vertex_decomposer *vpd) :
831 enumerator_base(dim,vpd) {
832 E_vertex = new evalue_p[dim];
835 virtual ~ienumerator_base() {
836 delete [] E_vertex;
839 evalue *E_num(int i, int d) {
840 return E_vertex[i + (dim-d)];
844 struct cumulator {
845 evalue *factor;
846 evalue *v;
847 dpoly_r *r;
849 cumulator(evalue *factor, evalue *v, dpoly_r *r) :
850 factor(factor), v(v), r(r) {}
852 void cumulate(barvinok_options *options);
854 virtual void add_term(const vector<int>& powers, evalue *f2) = 0;
855 virtual ~cumulator() {}
858 void cumulator::cumulate(barvinok_options *options)
860 evalue cum; // factor * 1 * E_num[0]/1 * (E_num[0]-1)/2 *...
861 evalue f;
862 evalue t; // E_num[0] - (m-1)
863 evalue *cst;
864 evalue mone;
866 if (options->lookup_table) {
867 value_init(mone.d);
868 evalue_set_si(&mone, -1, 1);
871 value_init(cum.d);
872 evalue_copy(&cum, factor);
873 value_init(f.d);
874 value_init(f.x.n);
875 value_set_si(f.d, 1);
876 value_set_si(f.x.n, 1);
877 value_init(t.d);
878 evalue_copy(&t, v);
880 if (!options->lookup_table) {
881 for (cst = &t; value_zero_p(cst->d); ) {
882 if (cst->x.p->type == fractional)
883 cst = &cst->x.p->arr[1];
884 else
885 cst = &cst->x.p->arr[0];
889 for (int m = 0; m < r->len; ++m) {
890 if (m > 0) {
891 if (m > 1) {
892 value_set_si(f.d, m);
893 emul(&f, &cum);
894 if (!options->lookup_table)
895 value_subtract(cst->x.n, cst->x.n, cst->d);
896 else
897 eadd(&mone, &t);
899 emul(&t, &cum);
901 dpoly_r_term_list& current = r->c[r->len-1-m];
902 dpoly_r_term_list::iterator j;
903 for (j = current.begin(); j != current.end(); ++j) {
904 if ((*j)->coeff == 0)
905 continue;
906 evalue *f2 = new evalue;
907 value_init(f2->d);
908 value_init(f2->x.n);
909 zz2value((*j)->coeff, f2->x.n);
910 zz2value(r->denom, f2->d);
911 emul(&cum, f2);
913 add_term((*j)->powers, f2);
916 free_evalue_refs(&f);
917 free_evalue_refs(&t);
918 free_evalue_refs(&cum);
919 if (options->lookup_table)
920 free_evalue_refs(&mone);
923 struct E_poly_term {
924 vector<int> powers;
925 evalue *E;
928 struct ie_cum : public cumulator {
929 vector<E_poly_term *> terms;
931 ie_cum(evalue *factor, evalue *v, dpoly_r *r) : cumulator(factor, v, r) {}
933 virtual void add_term(const vector<int>& powers, evalue *f2);
936 void ie_cum::add_term(const vector<int>& powers, evalue *f2)
938 int k;
939 for (k = 0; k < terms.size(); ++k) {
940 if (terms[k]->powers == powers) {
941 eadd(f2, terms[k]->E);
942 free_evalue_refs(f2);
943 delete f2;
944 break;
947 if (k >= terms.size()) {
948 E_poly_term *ET = new E_poly_term;
949 ET->powers = powers;
950 ET->E = f2;
951 terms.push_back(ET);
955 struct ienumerator : public signed_cone_consumer, public vertex_decomposer,
956 public ienumerator_base {
957 //Polyhedron *pVD;
958 mat_ZZ den;
959 mat_ZZ vertex;
960 mpq_t tcount;
961 Value tz;
963 ienumerator(Polyhedron *P, unsigned dim, Param_Polyhedron *PP) :
964 vertex_decomposer(PP, *this), ienumerator_base(dim, this) {
965 vertex.SetDims(1, dim);
967 den.SetDims(dim, dim);
968 mpq_init(tcount);
969 value_init(tz);
972 ~ienumerator() {
973 mpq_clear(tcount);
974 value_clear(tz);
977 virtual void handle(const signed_cone& sc, barvinok_options *options);
978 void reduce(evalue *factor, const mat_ZZ& num, const mat_ZZ& den_f,
979 barvinok_options *options);
982 void ienumerator::reduce(evalue *factor, const mat_ZZ& num, const mat_ZZ& den_f,
983 barvinok_options *options)
985 unsigned len = den_f.NumRows(); // number of factors in den
986 unsigned dim = num.NumCols();
987 assert(num.NumRows() == 1);
989 if (dim == 0) {
990 eadd(factor, vE[vert]);
991 return;
994 vec_ZZ den_s;
995 mat_ZZ den_r;
996 vec_ZZ num_s;
997 mat_ZZ num_p;
999 split_one(num, num_s, num_p, den_f, den_s, den_r);
1001 vec_ZZ den_p;
1002 den_p.SetLength(len);
1004 ZZ one;
1005 one = 1;
1006 normalize(one, num_s, num_p, den_s, den_p, den_r);
1007 if (one != 1)
1008 emul(&mone, factor);
1010 int only_param = 0;
1011 int no_param = 0;
1012 for (int k = 0; k < len; ++k) {
1013 if (den_p[k] == 0)
1014 ++no_param;
1015 else if (den_s[k] == 0)
1016 ++only_param;
1018 if (no_param == 0) {
1019 reduce(factor, num_p, den_r, options);
1020 } else {
1021 int k, l;
1022 mat_ZZ pden;
1023 pden.SetDims(only_param, dim-1);
1025 for (k = 0, l = 0; k < len; ++k)
1026 if (den_s[k] == 0)
1027 pden[l++] = den_r[k];
1029 for (k = 0; k < len; ++k)
1030 if (den_p[k] == 0)
1031 break;
1033 zz2value(num_s[0], tz);
1034 dpoly n(no_param, tz);
1035 zz2value(den_s[k], tz);
1036 dpoly D(no_param, tz, 1);
1037 for ( ; ++k < len; )
1038 if (den_p[k] == 0) {
1039 zz2value(den_s[k], tz);
1040 dpoly fact(no_param, tz, 1);
1041 D *= fact;
1044 dpoly_r * r = 0;
1045 // if no_param + only_param == len then all powers
1046 // below will be all zero
1047 if (no_param + only_param == len) {
1048 if (E_num(0, dim) != 0)
1049 r = new dpoly_r(n, len);
1050 else {
1051 mpq_set_si(tcount, 0, 1);
1052 one = 1;
1053 n.div(D, tcount, 1);
1055 if (value_notzero_p(mpq_numref(tcount))) {
1056 evalue f;
1057 value_init(f.d);
1058 value_init(f.x.n);
1059 value_assign(f.x.n, mpq_numref(tcount));
1060 value_assign(f.d, mpq_denref(tcount));
1061 emul(&f, factor);
1062 reduce(factor, num_p, pden, options);
1063 free_evalue_refs(&f);
1065 return;
1067 } else {
1068 for (k = 0; k < len; ++k) {
1069 if (den_s[k] == 0 || den_p[k] == 0)
1070 continue;
1072 zz2value(den_s[k], tz);
1073 dpoly pd(no_param-1, tz, 1);
1075 int l;
1076 for (l = 0; l < k; ++l)
1077 if (den_r[l] == den_r[k])
1078 break;
1080 if (r == 0)
1081 r = new dpoly_r(n, pd, l, len);
1082 else {
1083 dpoly_r *nr = new dpoly_r(r, pd, l, len);
1084 delete r;
1085 r = nr;
1089 dpoly_r *rc = r->div(D);
1090 delete r;
1091 r = rc;
1092 if (E_num(0, dim) == 0) {
1093 int common = pden.NumRows();
1094 dpoly_r_term_list& final = r->c[r->len-1];
1095 int rows;
1096 evalue t;
1097 evalue f;
1098 value_init(f.d);
1099 value_init(f.x.n);
1100 zz2value(r->denom, f.d);
1101 dpoly_r_term_list::iterator j;
1102 for (j = final.begin(); j != final.end(); ++j) {
1103 if ((*j)->coeff == 0)
1104 continue;
1105 rows = common;
1106 for (int k = 0; k < r->dim; ++k) {
1107 int n = (*j)->powers[k];
1108 if (n == 0)
1109 continue;
1110 pden.SetDims(rows+n, pden.NumCols());
1111 for (int l = 0; l < n; ++l)
1112 pden[rows+l] = den_r[k];
1113 rows += n;
1115 value_init(t.d);
1116 evalue_copy(&t, factor);
1117 zz2value((*j)->coeff, f.x.n);
1118 emul(&f, &t);
1119 reduce(&t, num_p, pden, options);
1120 free_evalue_refs(&t);
1122 free_evalue_refs(&f);
1123 } else {
1124 ie_cum cum(factor, E_num(0, dim), r);
1125 cum.cumulate(options);
1127 int common = pden.NumRows();
1128 int rows;
1129 for (int j = 0; j < cum.terms.size(); ++j) {
1130 rows = common;
1131 pden.SetDims(rows, pden.NumCols());
1132 for (int k = 0; k < r->dim; ++k) {
1133 int n = cum.terms[j]->powers[k];
1134 if (n == 0)
1135 continue;
1136 pden.SetDims(rows+n, pden.NumCols());
1137 for (int l = 0; l < n; ++l)
1138 pden[rows+l] = den_r[k];
1139 rows += n;
1141 reduce(cum.terms[j]->E, num_p, pden, options);
1142 free_evalue_refs(cum.terms[j]->E);
1143 delete cum.terms[j]->E;
1144 delete cum.terms[j];
1147 delete r;
1151 static int type_offset(enode *p)
1153 return p->type == fractional ? 1 :
1154 p->type == flooring ? 1 : 0;
1157 static int edegree(evalue *e)
1159 int d = 0;
1160 enode *p;
1162 if (value_notzero_p(e->d))
1163 return 0;
1165 p = e->x.p;
1166 int i = type_offset(p);
1167 if (p->size-i-1 > d)
1168 d = p->size - i - 1;
1169 for (; i < p->size; i++) {
1170 int d2 = edegree(&p->arr[i]);
1171 if (d2 > d)
1172 d = d2;
1174 return d;
1177 void ienumerator::handle(const signed_cone& sc, barvinok_options *options)
1179 assert(sc.det == 1);
1180 assert(sc.rays.NumRows() == dim);
1182 lattice_point(V, sc.rays, vertex[0], E_vertex, options);
1184 den = sc.rays;
1186 evalue one;
1187 value_init(one.d);
1188 evalue_set_si(&one, sc.sign, 1);
1189 reduce(&one, vertex, den, options);
1190 free_evalue_refs(&one);
1192 for (int i = 0; i < dim; ++i)
1193 if (E_vertex[i]) {
1194 free_evalue_refs(E_vertex[i]);
1195 delete E_vertex[i];
1199 struct bfenumerator : public vertex_decomposer, public bf_base,
1200 public ienumerator_base {
1201 evalue *factor;
1203 bfenumerator(Polyhedron *P, unsigned dim, Param_Polyhedron *PP) :
1204 vertex_decomposer(PP, *this),
1205 bf_base(dim), ienumerator_base(dim, this) {
1206 lower = 0;
1207 factor = NULL;
1210 ~bfenumerator() {
1213 virtual void handle(const signed_cone& sc, barvinok_options *options);
1214 virtual void base(mat_ZZ& factors, bfc_vec& v);
1216 bfc_term_base* new_bf_term(int len) {
1217 bfe_term* t = new bfe_term(len);
1218 return t;
1221 virtual void set_factor(bfc_term_base *t, int k, int change) {
1222 bfe_term* bfet = static_cast<bfe_term *>(t);
1223 factor = bfet->factors[k];
1224 assert(factor != NULL);
1225 bfet->factors[k] = NULL;
1226 if (change)
1227 emul(&mone, factor);
1230 virtual void set_factor(bfc_term_base *t, int k, mpq_t &q, int change) {
1231 bfe_term* bfet = static_cast<bfe_term *>(t);
1232 factor = bfet->factors[k];
1233 assert(factor != NULL);
1234 bfet->factors[k] = NULL;
1236 evalue f;
1237 value_init(f.d);
1238 value_init(f.x.n);
1239 if (change)
1240 value_oppose(f.x.n, mpq_numref(q));
1241 else
1242 value_assign(f.x.n, mpq_numref(q));
1243 value_assign(f.d, mpq_denref(q));
1244 emul(&f, factor);
1245 free_evalue_refs(&f);
1248 virtual void set_factor(bfc_term_base *t, int k, const QQ& c, int change) {
1249 bfe_term* bfet = static_cast<bfe_term *>(t);
1251 factor = new evalue;
1253 evalue f;
1254 value_init(f.d);
1255 value_init(f.x.n);
1256 zz2value(c.n, f.x.n);
1257 if (change)
1258 value_oppose(f.x.n, f.x.n);
1259 zz2value(c.d, f.d);
1261 value_init(factor->d);
1262 evalue_copy(factor, bfet->factors[k]);
1263 emul(&f, factor);
1264 free_evalue_refs(&f);
1267 void set_factor(evalue *f, int change) {
1268 if (change)
1269 emul(&mone, f);
1270 factor = f;
1273 virtual void insert_term(bfc_term_base *t, int i) {
1274 bfe_term* bfet = static_cast<bfe_term *>(t);
1275 int len = t->terms.NumRows()-1; // already increased by one
1277 bfet->factors.resize(len+1);
1278 for (int j = len; j > i; --j) {
1279 bfet->factors[j] = bfet->factors[j-1];
1280 t->terms[j] = t->terms[j-1];
1282 bfet->factors[i] = factor;
1283 factor = NULL;
1286 virtual void update_term(bfc_term_base *t, int i) {
1287 bfe_term* bfet = static_cast<bfe_term *>(t);
1289 eadd(factor, bfet->factors[i]);
1290 free_evalue_refs(factor);
1291 delete factor;
1294 virtual bool constant_vertex(int dim) { return E_num(0, dim) == 0; }
1296 virtual void cum(bf_reducer *bfr, bfc_term_base *t, int k, dpoly_r *r,
1297 barvinok_options *options);
1300 enumerator_base *enumerator_base::create(Polyhedron *P, unsigned dim,
1301 Param_Polyhedron *PP,
1302 barvinok_options *options)
1304 enumerator_base *eb;
1306 if (options->incremental_specialization == BV_SPECIALIZATION_BF)
1307 eb = new bfenumerator(P, dim, PP);
1308 else if (options->incremental_specialization == BV_SPECIALIZATION_DF)
1309 eb = new ienumerator(P, dim, PP);
1310 else
1311 eb = new enumerator(P, dim, PP);
1313 return eb;
1316 struct bfe_cum : public cumulator {
1317 bfenumerator *bfe;
1318 bfc_term_base *told;
1319 int k;
1320 bf_reducer *bfr;
1322 bfe_cum(evalue *factor, evalue *v, dpoly_r *r, bf_reducer *bfr,
1323 bfc_term_base *t, int k, bfenumerator *e) :
1324 cumulator(factor, v, r), told(t), k(k),
1325 bfr(bfr), bfe(e) {
1328 virtual void add_term(const vector<int>& powers, evalue *f2);
1331 void bfe_cum::add_term(const vector<int>& powers, evalue *f2)
1333 bfr->update_powers(powers);
1335 bfc_term_base * t = bfe->find_bfc_term(bfr->vn, bfr->npowers, bfr->nnf);
1336 bfe->set_factor(f2, bfr->l_changes % 2);
1337 bfe->add_term(t, told->terms[k], bfr->l_extra_num);
1340 void bfenumerator::cum(bf_reducer *bfr, bfc_term_base *t, int k,
1341 dpoly_r *r, barvinok_options *options)
1343 bfe_term* bfet = static_cast<bfe_term *>(t);
1344 bfe_cum cum(bfet->factors[k], E_num(0, bfr->d), r, bfr, t, k, this);
1345 cum.cumulate(options);
1348 void bfenumerator::base(mat_ZZ& factors, bfc_vec& v)
1350 for (int i = 0; i < v.size(); ++i) {
1351 assert(v[i]->terms.NumRows() == 1);
1352 evalue *factor = static_cast<bfe_term *>(v[i])->factors[0];
1353 eadd(factor, vE[vert]);
1354 delete v[i];
1358 void bfenumerator::handle(const signed_cone& sc, barvinok_options *options)
1360 assert(sc.det == 1);
1361 assert(sc.rays.NumRows() == enumerator_base::dim);
1363 bfe_term* t = new bfe_term(enumerator_base::dim);
1364 vector< bfc_term_base * > v;
1365 v.push_back(t);
1367 t->factors.resize(1);
1369 t->terms.SetDims(1, enumerator_base::dim);
1370 lattice_point(V, sc.rays, t->terms[0], E_vertex, options);
1372 // the elements of factors are always lexpositive
1373 mat_ZZ factors;
1374 int s = setup_factors(sc.rays, factors, t, sc.sign);
1376 t->factors[0] = new evalue;
1377 value_init(t->factors[0]->d);
1378 evalue_set_si(t->factors[0], s, 1);
1379 reduce(factors, v, options);
1381 for (int i = 0; i < enumerator_base::dim; ++i)
1382 if (E_vertex[i]) {
1383 free_evalue_refs(E_vertex[i]);
1384 delete E_vertex[i];
1388 static evalue* barvinok_enumerate_ev_f(Polyhedron *P, Polyhedron* C,
1389 barvinok_options *options);
1391 /* Destroys C */
1392 static evalue* barvinok_enumerate_cst(Polyhedron *P, Polyhedron* C,
1393 struct barvinok_options *options)
1395 evalue *eres;
1397 if (emptyQ2(C)) {
1398 Polyhedron_Free(C);
1399 return evalue_zero();
1402 ALLOC(evalue, eres);
1403 value_init(eres->d);
1404 value_set_si(eres->d, 0);
1405 eres->x.p = new_enode(partition, 2, C->Dimension);
1406 EVALUE_SET_DOMAIN(eres->x.p->arr[0],
1407 DomainConstraintSimplify(C, options->MaxRays));
1408 value_set_si(eres->x.p->arr[1].d, 1);
1409 value_init(eres->x.p->arr[1].x.n);
1410 if (emptyQ2(P))
1411 value_set_si(eres->x.p->arr[1].x.n, 0);
1412 else
1413 barvinok_count_with_options(P, &eres->x.p->arr[1].x.n, options);
1415 return eres;
1418 static evalue* enumerate(Polyhedron *P, Polyhedron* C,
1419 struct barvinok_options *options)
1421 Polyhedron *next;
1422 Polyhedron *Porig = P;
1423 Polyhedron *Corig = C;
1424 Polyhedron *CEq = NULL, *rVD;
1425 int r = 0;
1426 unsigned nparam = C->Dimension;
1427 evalue *eres;
1428 Matrix *CP = NULL;
1430 evalue factor;
1431 value_init(factor.d);
1432 evalue_set_si(&factor, 1, 1);
1434 /* for now */
1435 POL_ENSURE_FACETS(P);
1436 POL_ENSURE_VERTICES(P);
1437 POL_ENSURE_FACETS(C);
1438 POL_ENSURE_VERTICES(C);
1440 if (C->Dimension == 0 || emptyQ(P) || emptyQ(C)) {
1441 constant:
1442 if (CEq == Porig)
1443 CEq = Polyhedron_Copy(CEq);
1444 eres = barvinok_enumerate_cst(P, CEq ? CEq : Polyhedron_Copy(C), options);
1445 out:
1446 if (CP) {
1447 evalue_backsubstitute(eres, CP, options->MaxRays);
1448 Matrix_Free(CP);
1451 emul(&factor, eres);
1452 if (options->approximation_method == BV_APPROX_DROP) {
1453 if (options->polynomial_approximation == BV_APPROX_SIGN_UPPER)
1454 evalue_frac2polynomial(eres, 1, options->MaxRays);
1455 if (options->polynomial_approximation == BV_APPROX_SIGN_LOWER)
1456 evalue_frac2polynomial(eres, -1, options->MaxRays);
1457 if (options->polynomial_approximation == BV_APPROX_SIGN_APPROX)
1458 evalue_frac2polynomial(eres, 0, options->MaxRays);
1460 reduce_evalue(eres);
1461 free_evalue_refs(&factor);
1462 if (P != Porig)
1463 Domain_Free(P);
1464 if (C != Corig)
1465 Polyhedron_Free(C);
1467 return eres;
1469 if (Polyhedron_is_unbounded(P, nparam, options->MaxRays))
1470 goto constant;
1472 if (P->NbEq != 0) {
1473 Matrix *f;
1474 P = remove_equalities_p(Polyhedron_Copy(P), P->Dimension-nparam, &f,
1475 options->MaxRays);
1476 mask(f, &factor, options);
1477 Matrix_Free(f);
1479 if (P->Dimension == nparam) {
1480 CEq = P;
1481 P = Universe_Polyhedron(0);
1482 goto constant;
1484 if (P->NbEq != 0 || C->NbEq != 0) {
1485 Polyhedron *Q = P;
1486 Polyhedron *D = C;
1487 remove_all_equalities(&P, &C, &CP, NULL, nparam, options->MaxRays);
1488 if (C != D && D != Corig)
1489 Polyhedron_Free(D);
1490 if (P != Q && Q != Porig)
1491 Domain_Free(Q);
1492 eres = enumerate(P, C, options);
1493 goto out;
1496 Polyhedron *T = Polyhedron_Factor(P, nparam, NULL, options->MaxRays);
1497 if (T || (P->Dimension == nparam+1)) {
1498 Polyhedron *Q;
1499 Polyhedron *C2;
1500 for (Q = T ? T : P; Q; Q = Q->next) {
1501 Polyhedron *next = Q->next;
1502 Q->next = NULL;
1504 Polyhedron *QC = Q;
1505 if (Q->Dimension != C->Dimension)
1506 QC = Polyhedron_Project(Q, nparam);
1508 C2 = C;
1509 C = DomainIntersection(C, QC, options->MaxRays);
1510 if (C2 != Corig)
1511 Polyhedron_Free(C2);
1512 if (QC != Q)
1513 Polyhedron_Free(QC);
1515 Q->next = next;
1518 if (T) {
1519 if (P != Porig)
1520 Polyhedron_Free(P);
1521 P = T;
1522 if (T->Dimension == C->Dimension) {
1523 P = T->next;
1524 T->next = NULL;
1525 Polyhedron_Free(T);
1529 next = P->next;
1530 P->next = NULL;
1531 eres = barvinok_enumerate_ev_f(P, C, options);
1532 P->next = next;
1534 if (P->next) {
1535 Polyhedron *Q;
1536 evalue *f;
1538 for (Q = P->next; Q; Q = Q->next) {
1539 Polyhedron *next = Q->next;
1540 Q->next = NULL;
1542 f = barvinok_enumerate_ev_f(Q, C, options);
1543 emul(f, eres);
1544 evalue_free(f);
1546 Q->next = next;
1550 goto out;
1553 evalue* barvinok_enumerate_with_options(Polyhedron *P, Polyhedron* C,
1554 struct barvinok_options *options)
1556 Polyhedron *next, *Cnext, *C1;
1557 Polyhedron *Corig = C;
1558 evalue *eres;
1560 if (P->next)
1561 fprintf(stderr,
1562 "barvinok_enumerate: input is a union; only first polyhedron is enumerated\n");
1564 if (C->next)
1565 fprintf(stderr,
1566 "barvinok_enumerate: context is a union; only first polyhedron is considered\n");
1568 Cnext = C->next;
1569 C->next = NULL;
1570 C1 = Polyhedron_Project(P, C->Dimension);
1571 C = DomainIntersection(C, C1, options->MaxRays);
1572 Polyhedron_Free(C1);
1573 next = P->next;
1574 P->next = NULL;
1576 if (options->approximation_method == BV_APPROX_BERNOULLI)
1577 eres = Bernoulli_sum(P, C, options);
1578 else
1579 eres = enumerate(P, C, options);
1580 Domain_Free(C);
1582 P->next= next;
1583 Corig->next = Cnext;
1585 return eres;
1588 evalue* barvinok_enumerate_ev(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
1590 evalue *E;
1591 barvinok_options *options = barvinok_options_new_with_defaults();
1592 options->MaxRays = MaxRays;
1593 E = barvinok_enumerate_with_options(P, C, options);
1594 barvinok_options_free(options);
1595 return E;
1598 evalue *Param_Polyhedron_Enumerate(Param_Polyhedron *PP, Polyhedron *P,
1599 Polyhedron *C,
1600 struct barvinok_options *options)
1602 evalue *eres;
1603 Param_Domain *D;
1604 unsigned nparam = C->Dimension;
1605 unsigned dim = P->Dimension - nparam;
1607 int nd;
1608 for (nd = 0, D=PP->D; D; ++nd, D=D->next);
1609 evalue_section *s = new evalue_section[nd];
1611 enumerator_base *et = NULL;
1612 try_again:
1613 if (et)
1614 delete et;
1616 et = enumerator_base::create(P, dim, PP, options);
1618 Polyhedron *TC = true_context(P, C, options->MaxRays);
1619 FORALL_REDUCED_DOMAIN(PP, TC, nd, options, i, D, rVD)
1620 Param_Vertices *V;
1622 s[i].E = evalue_zero();
1623 s[i].D = rVD;
1625 FORALL_PVertex_in_ParamPolyhedron(V,D,PP) // _i is internal counter
1626 if (!et->vE[_i])
1627 try {
1628 et->decompose_at(V, _i, options);
1629 } catch (OrthogonalException &e) {
1630 FORALL_REDUCED_DOMAIN_RESET;
1631 for (; i >= 0; --i) {
1632 evalue_free(s[i].E);
1633 Domain_Free(s[i].D);
1635 goto try_again;
1637 eadd(et->vE[_i] , s[i].E);
1638 END_FORALL_PVertex_in_ParamPolyhedron;
1639 evalue_range_reduction_in_domain(s[i].E, rVD);
1640 END_FORALL_REDUCED_DOMAIN
1641 Polyhedron_Free(TC);
1643 delete et;
1644 eres = evalue_from_section_array(s, nd);
1645 delete [] s;
1647 return eres;
1650 static evalue* barvinok_enumerate_ev_f(Polyhedron *P, Polyhedron* C,
1651 barvinok_options *options)
1653 unsigned nparam = C->Dimension;
1654 bool do_scale = options->approximation_method == BV_APPROX_SCALE;
1656 if (options->approximation_method == BV_APPROX_VOLUME)
1657 return Param_Polyhedron_Volume(P, C, options);
1659 if (P->Dimension - nparam == 1 && !do_scale)
1660 return ParamLine_Length(P, C, options);
1662 Param_Polyhedron *PP = NULL;
1663 evalue *eres;
1665 if (do_scale) {
1666 eres = scale_bound(P, C, options);
1667 if (eres)
1668 return eres;
1671 PP = Polyhedron2Param_Polyhedron(P, C, options);
1673 if (do_scale)
1674 eres = scale(PP, P, C, options);
1675 else
1676 eres = Param_Polyhedron_Enumerate(PP, P, C, options);
1678 if (PP)
1679 Param_Polyhedron_Free(PP);
1681 return eres;
1684 Enumeration* barvinok_enumerate(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
1686 evalue *EP = barvinok_enumerate_ev(P, C, MaxRays);
1688 return partition2enumeration(EP);
1691 evalue* barvinok_enumerate_union(Polyhedron *D, Polyhedron* C, unsigned MaxRays)
1693 evalue *EP;
1694 gen_fun *gf = barvinok_enumerate_union_series(D, C, MaxRays);
1695 EP = *gf;
1696 delete gf;
1697 return EP;
1700 evalue *barvinok_summate(evalue *e, int nvar, struct barvinok_options *options)
1702 if (options->summation == BV_SUM_EULER)
1703 return euler_summate(e, nvar, options);
1704 else if (options->summation == BV_SUM_BERNOULLI)
1705 return Bernoulli_sum_evalue(e, nvar, options);
1706 else
1707 return evalue_sum(e, nvar, options->MaxRays);