?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Bernoulli_sum_evalue: exploit integer bounds to reduce the number of cases
[barvinok.git] / barvinok.cc
blob73366398ca7d9d0ef02cb768613b8ad05c6dbafe
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"
29
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;
40
41 #define ALLOC(t,p) p = (t*)malloc(sizeof(*p))
42
43 class dpoly_n {
44 public:
45 Matrix *coeff;
46 ~dpoly_n() {
47 Matrix_Free(coeff);
48 }
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);
64 }
65 value_clear(d0);
66 value_clear(one);
67 }
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]);
81 }
82 value_multiply(c->p[i][len], c->p[i][len], d.coeff->p[0]);
83 }
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);
93 }
94 };
95
96 const int MAX_TRY=10;
97 /*
98 * Searches for a vector that is not orthogonal to any
99 * of the rays in rays.
100 */
101 static void nonorthog(mat_ZZ& rays, vec_ZZ& lambda)
102 {
103 int dim = rays.NumCols();
104 bool found = false;
105 lambda.SetLength(dim);
106 if (dim == 0)
107 return;
108
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;
115 }
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;
123 }
124 }
125 }
126 assert(found);
127 }
128
129 static void add_rays(mat_ZZ& rays, Polyhedron *i, int *r, int nvar = -1,
130 bool all = false)
131 {
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);
141 }
142 }
143
144 struct bfe_term : public bfc_term_base {
145 vector<evalue *> factors;
146
147 bfe_term(int len) : bfc_term_base(len) {
148 }
149
150 ~bfe_term() {
151 for (int i = 0; i < factors.size(); ++i) {
152 if (!factors[i])
153 continue;
154 free_evalue_refs(factors[i]);
155 delete factors[i];
156 }
157 }
158 };
159
160 static void print_int_vector(int *v, int len, const char *name)
161 {
162 cerr << name << endl;
163 for (int j = 0; j < len; ++j) {
164 cerr << v[j] << " ";
165 }
166 cerr << endl;
167 }
168
169 static void print_bfc_terms(mat_ZZ& factors, bfc_vec& v)
170 {
171 cerr << endl;
172 cerr << "factors" << endl;
173 cerr << factors << endl;
174 for (int i = 0; i < v.size(); ++i) {
175 cerr << "term: " << i << endl;
176 print_int_vector(v[i]->powers, factors.NumRows(), "powers");
177 cerr << "terms" << endl;
178 cerr << v[i]->terms << endl;
179 bfc_term* bfct = static_cast<bfc_term *>(v[i]);
180 cerr << bfct->c << endl;
181 }
182 }
183
184 static void print_bfe_terms(mat_ZZ& factors, bfc_vec& v)
185 {
186 cerr << endl;
187 cerr << "factors" << endl;
188 cerr << factors << endl;
189 for (int i = 0; i < v.size(); ++i) {
190 cerr << "term: " << i << endl;
191 print_int_vector(v[i]->powers, factors.NumRows(), "powers");
192 cerr << "terms" << endl;
193 cerr << v[i]->terms << endl;
194 bfe_term* bfet = static_cast<bfe_term *>(v[i]);
195 for (int j = 0; j < v[i]->terms.NumRows(); ++j) {
196 const char * test[] = {"a", "b"};
197 print_evalue(stderr, bfet->factors[j], test);
198 fprintf(stderr, "\n");
199 }
200 }
201 }
202
203 struct bfcounter : public bfcounter_base {
204 mpq_t count;
205 Value tz;
206
207 bfcounter(unsigned dim) : bfcounter_base(dim) {
208 mpq_init(count);
209 lower = 1;
210 value_init(tz);
211 }
212 ~bfcounter() {
213 mpq_clear(count);
214 value_clear(tz);
215 }
216 virtual void base(mat_ZZ& factors, bfc_vec& v);
217 virtual void get_count(Value *result) {
218 assert(value_one_p(&count[0]._mp_den));
219 value_assign(*result, &count[0]._mp_num);
220 }
221 };
222
223 void bfcounter::base(mat_ZZ& factors, bfc_vec& v)
224 {
225 unsigned nf = factors.NumRows();
226
227 for (int i = 0; i < v.size(); ++i) {
228 bfc_term* bfct = static_cast<bfc_term *>(v[i]);
229 int total_power = 0;
230 // factor is always positive, so we always
231 // change signs
232 for (int k = 0; k < nf; ++k)
233 total_power += v[i]->powers[k];
234
235 int j;
236 for (j = 0; j < nf; ++j)
237 if (v[i]->powers[j] > 0)
238 break;
239
240 zz2value(factors[j][0], tz);
241 dpoly D(total_power, tz, 1);
242 for (int k = 1; k < v[i]->powers[j]; ++k) {
243 zz2value(factors[j][0], tz);
244 dpoly fact(total_power, tz, 1);
245 D *= fact;
246 }
247 for ( ; ++j < nf; )
248 for (int k = 0; k < v[i]->powers[j]; ++k) {
249 zz2value(factors[j][0], tz);
250 dpoly fact(total_power, tz, 1);
251 D *= fact;
252 }
253
254 for (int k = 0; k < v[i]->terms.NumRows(); ++k) {
255 zz2value(v[i]->terms[k][0], tz);
256 dpoly n(total_power, tz);
257 mpq_set_si(tcount, 0, 1);
258 n.div(D, tcount, 1);
259 if (total_power % 2)
260 bfct->c[k].n = -bfct->c[k].n;
261 zz2value(bfct->c[k].n, tn);
262 zz2value(bfct->c[k].d, td);
263
264 mpz_mul(mpq_numref(tcount), mpq_numref(tcount), tn);
265 mpz_mul(mpq_denref(tcount), mpq_denref(tcount), td);
266 mpq_canonicalize(tcount);
267 mpq_add(count, count, tcount);
268 }
269 delete v[i];
270 }
271 }
272
273
274 /* Check whether the polyhedron is unbounded and if so,
275 * check whether it has any (and therefore an infinite number of)
276 * integer points.
277 * If one of the vertices is integer, then we are done.
278 * Otherwise, transform the polyhedron such that one of the rays
279 * is the first unit vector and cut it off at a height that ensures
280 * that if the whole polyhedron has any points, then the remaining part
281 * has integer points. In particular we add the largest coefficient
282 * of a ray to the highest vertex (rounded up).
283 */
284 static bool Polyhedron_is_infinite(Polyhedron *P, Value* result,
285 barvinok_options *options)
286 {
287 int r = 0;
288 Matrix *M, *M2;
289 Value c, tmp;
290 Value g;
291 bool first;
292 Vector *v;
293 Value offset, size;
294 Polyhedron *R;
295
296 if (P->NbBid == 0)
297 for (; r < P->NbRays; ++r)
298 if (value_zero_p(P->Ray[r][P->Dimension+1]))
299 break;
300 if (P->NbBid == 0 && r == P->NbRays)
301 return false;
302
303 if (options->count_sample_infinite) {
304 Vector *sample;
305
306 sample = Polyhedron_Sample(P, options);
307 if (!sample)
308 value_set_si(*result, 0);
309 else {
310 value_set_si(*result, -1);
311 Vector_Free(sample);
312 }
313 return true;
314 }
315
316 for (int i = 0; i < P->NbRays; ++i)
317 if (value_one_p(P->Ray[i][1+P->Dimension])) {
318 value_set_si(*result, -1);
319 return true;
320 }
321
322 value_init(g);
323 M = Matrix_Alloc(P->Dimension+1, P->Dimension+1);
324 Vector_Gcd(P->Ray[r]+1, P->Dimension, &g);
325 Vector_AntiScale(P->Ray[r]+1, M->p[0], g, P->Dimension+1);
326 int ok = unimodular_complete(M, 1);
327 assert(ok);
328 value_set_si(M->p[P->Dimension][P->Dimension], 1);
329 M2 = Transpose(M);
330 Matrix_Free(M);
331 P = Polyhedron_Preimage(P, M2, 0);
332 Matrix_Free(M2);
333 value_clear(g);
334
335 first = true;
336 value_init(offset);
337 value_init(size);
338 value_init(tmp);
339 value_set_si(size, 0);
340
341 for (int i = 0; i < P->NbBid; ++i) {
342 value_absolute(tmp, P->Ray[i][1]);
343 if (value_gt(tmp, size))
344 value_assign(size, tmp);
345 }
346 for (int i = P->NbBid; i < P->NbRays; ++i) {
347 if (value_zero_p(P->Ray[i][P->Dimension+1])) {
348 if (value_gt(P->Ray[i][1], size))
349 value_assign(size, P->Ray[i][1]);
350 continue;
351 }
352 mpz_cdiv_q(tmp, P->Ray[i][1], P->Ray[i][P->Dimension+1]);
353 if (first || value_gt(tmp, offset)) {
354 value_assign(offset, tmp);
355 first = false;
356 }
357 }
358 value_addto(offset, offset, size);
359 value_clear(size);
360 value_clear(tmp);
361
362 v = Vector_Alloc(P->Dimension+2);
363 value_set_si(v->p[0], 1);
364 value_set_si(v->p[1], -1);
365 value_assign(v->p[1+P->Dimension], offset);
366 R = AddConstraints(v->p, 1, P, options->MaxRays);
367 Polyhedron_Free(P);
368 P = R;
369
370 value_clear(offset);
371 Vector_Free(v);
372
373 value_init(c);
374 barvinok_count_with_options(P, &c, options);
375 Polyhedron_Free(P);
376 if (value_zero_p(c))
377 value_set_si(*result, 0);
378 else
379 value_set_si(*result, -1);
380 value_clear(c);
381
382 return true;
383 }
384
385 static void evalue2value(evalue *e, Value *v)
386 {
387 if (EVALUE_IS_ZERO(*e)) {
388 value_set_si(*v, 0);
389 return;
390 }
391
392 if (value_notzero_p(e->d)) {
393 assert(value_one_p(e->d));
394 value_assign(*v, e->x.n);
395 return;
396 }
397
398 assert(e->x.p->type == partition);
399 assert(e->x.p->size == 2);
400 assert(EVALUE_DOMAIN(e->x.p->arr[0])->Dimension == 0);
401 evalue2value(&e->x.p->arr[1], v);
402 }
403
404 static void barvinok_count_f(Polyhedron *P, Value* result,
405 barvinok_options *options);
406
407 void barvinok_count_with_options(Polyhedron *P, Value* result,
408 struct barvinok_options *options)
409 {
410 unsigned dim;
411 int allocated = 0;
412 Polyhedron *Q;
413 bool infinite = false;
414
415 if (P->next)
416 fprintf(stderr,
417 "barvinok_count: input is a union; only first polyhedron is counted\n");
418
419 if (emptyQ2(P)) {
420 value_set_si(*result, 0);
421 return;
422 }
423 if (P->NbEq != 0) {
424 Q = NULL;
425 do {
426 P = remove_equalities(P, options->MaxRays);
427 P = DomainConstraintSimplify(P, options->MaxRays);
428 if (Q)
429 Polyhedron_Free(Q);
430 Q = P;
431 } while (!emptyQ(P) && P->NbEq != 0);
432 if (emptyQ(P)) {
433 Polyhedron_Free(P);
434 value_set_si(*result, 0);
435 return;
436 }
437 allocated = 1;
438 }
439 if (Polyhedron_is_infinite(P, result, options)) {
440 if (allocated)
441 Polyhedron_Free(P);
442 return;
443 }
444 if (P->Dimension == 0) {
445 /* Test whether the constraints are satisfied */
446 POL_ENSURE_VERTICES(P);
447 value_set_si(*result, !emptyQ(P));
448 if (allocated)
449 Polyhedron_Free(P);
450 return;
451 }
452 if (options->summation == BV_SUM_BERNOULLI) {
453 Polyhedron *C = Universe_Polyhedron(0);
454 evalue *sum = Bernoulli_sum(P, C, options);
455 Polyhedron_Free(C);
456 evalue2value(sum, result);
457 evalue_free(sum);
458 return;
459 }
460 Q = Polyhedron_Factor(P, 0, NULL, options->MaxRays);
461 if (Q) {
462 if (allocated)
463 Polyhedron_Free(P);
464 P = Q;
465 allocated = 1;
466 }
467
468 barvinok_count_f(P, result, options);
469 if (value_neg_p(*result))
470 infinite = true;
471 if (Q && P->next && value_notzero_p(*result)) {
472 Value factor;
473 value_init(factor);
474
475 for (Q = P->next; Q; Q = Q->next) {
476 barvinok_count_f(Q, &factor, options);
477 if (value_neg_p(factor)) {
478 infinite = true;
479 continue;
480 } else if (Q->next && value_zero_p(factor)) {
481 value_set_si(*result, 0);
482 break;
483 }
484 value_multiply(*result, *result, factor);
485 }
486
487 value_clear(factor);
488 }
489
490 if (allocated)
491 Domain_Free(P);
492 if (infinite)
493 value_set_si(*result, -1);
494 }
495
496 void barvinok_count(Polyhedron *P, Value* result, unsigned NbMaxCons)
497 {
498 barvinok_options *options = barvinok_options_new_with_defaults();
499 options->MaxRays = NbMaxCons;
500 barvinok_count_with_options(P, result, options);
501 barvinok_options_free(options);
502 }
503
504 static void barvinok_count_f(Polyhedron *P, Value* result,
505 barvinok_options *options)
506 {
507 if (emptyQ2(P)) {
508 value_set_si(*result, 0);
509 return;
510 }
511
512 if (P->Dimension == 1)
513 return Line_Length(P, result);
514
515 int c = P->NbConstraints;
516 POL_ENSURE_FACETS(P);
517 if (c != P->NbConstraints || P->NbEq != 0) {
518 Polyhedron *next = P->next;
519 P->next = NULL;
520 barvinok_count_with_options(P, result, options);
521 P->next = next;
522 return;
523 }
524
525 POL_ENSURE_VERTICES(P);
526
527 if (Polyhedron_is_infinite(P, result, options))
528 return;
529
530 np_base *cnt;
531 if (options->incremental_specialization == BV_SPECIALIZATION_BF)
532 cnt = new bfcounter(P->Dimension);
533 else if (options->incremental_specialization == BV_SPECIALIZATION_DF)
534 cnt = new icounter(P->Dimension);
535 else if (options->incremental_specialization == BV_SPECIALIZATION_TODD)
536 cnt = new tcounter(P->Dimension, options->max_index);
537 else
538 cnt = new counter(P->Dimension, options->max_index);
539 cnt->start(P, options);
540
541 cnt->get_count(result);
542 delete cnt;
543 }
544
545 static void uni_polynom(int param, Vector *c, evalue *EP)
546 {
547 unsigned dim = c->Size-2;
548 value_init(EP->d);
549 value_set_si(EP->d,0);
550 EP->x.p = new_enode(polynomial, dim+1, param+1);
551 for (int j = 0; j <= dim; ++j)
552 evalue_set(&EP->x.p->arr[j], c->p[j], c->p[dim+1]);
553 }
554
555 typedef evalue * evalue_p;
556
557 struct enumerator_base {
558 unsigned dim;
559 evalue ** vE;
560 evalue mone;
561 vertex_decomposer *vpd;
562
563 enumerator_base(unsigned dim, vertex_decomposer *vpd)
564 {
565 this->dim = dim;
566 this->vpd = vpd;
567
568 vE = new evalue_p[vpd->PP->nbV];
569 for (int j = 0; j < vpd->PP->nbV; ++j)
570 vE[j] = 0;
571
572 value_init(mone.d);
573 evalue_set_si(&mone, -1, 1);
574 }
575
576 void decompose_at(Param_Vertices *V, int _i, barvinok_options *options) {
577 //this->pVD = pVD;
578
579 vE[_i] = new evalue;
580 value_init(vE[_i]->d);
581 evalue_set_si(vE[_i], 0, 1);
582
583 vpd->decompose_at_vertex(V, _i, options);
584 }
585
586 virtual ~enumerator_base() {
587 for (int j = 0; j < vpd->PP->nbV; ++j)
588 if (vE[j]) {
589 free_evalue_refs(vE[j]);
590 delete vE[j];
591 }
592 delete [] vE;
593
594 free_evalue_refs(&mone);
595 }
596
597 static enumerator_base *create(Polyhedron *P, unsigned dim,
598 Param_Polyhedron *PP,
599 barvinok_options *options);
600 };
601
602 struct enumerator : public signed_cone_consumer, public vertex_decomposer,
603 public enumerator_base {
604 vec_ZZ lambda;
605 vec_ZZ den;
606 term_info num;
607 Vector *c;
608 mpq_t count;
609 Value tz;
610
611 enumerator(Polyhedron *P, unsigned dim, Param_Polyhedron *PP) :
612 vertex_decomposer(PP, *this), enumerator_base(dim, this) {
613 randomvector(P, lambda, dim);
614 den.SetLength(dim);
615 c = Vector_Alloc(dim+2);
616
617 mpq_init(count);
618 value_init(tz);
619 }
620
621 ~enumerator() {
622 mpq_clear(count);
623 Vector_Free(c);
624 value_clear(tz);
625 }
626
627 virtual void handle(const signed_cone& sc, barvinok_options *options);
628 };
629
630 void enumerator::handle(const signed_cone& sc, barvinok_options *options)
631 {
632 int sign = sc.sign;
633 int r = 0;
634 assert(sc.rays.NumRows() == dim);
635 for (int k = 0; k < dim; ++k) {
636 if (lambda * sc.rays[k] == 0)
637 throw Orthogonal;
638 }
639
640 lattice_point(V, sc.rays, lambda, &num, sc.det, options);
641 den = sc.rays * lambda;
642
643 if (dim % 2)
644 sign = -sign;
645
646 zz2value(den[0], tz);
647 dpoly n(dim, tz, 1);
648 for (int k = 1; k < dim; ++k) {
649 zz2value(den[k], tz);
650 dpoly fact(dim, tz, 1);
651 n *= fact;
652 }
653 if (num.E != NULL) {
654 dpoly_n d(dim);
655 d.div(n, c, sign);
656 for (unsigned long i = 0; i < sc.det; ++i) {
657 evalue *EV = evalue_polynomial(c, num.E[i]);
658 eadd(EV, vE[vert]);
659 evalue_free(EV);
660 free_evalue_refs(num.E[i]);
661 delete num.E[i];
662 }
663 delete [] num.E;
664 } else {
665 mpq_set_si(count, 0, 1);
666 if (num.constant.length() == 1) {
667 zz2value(num.constant[0], tz);
668 dpoly d(dim, tz);
669 d.div(n, count, sign);
670 } else {
671 dpoly_n d(dim);
672 d.div(n, c, sign);
673 Value x, sum, acc;
674 value_init(x);
675 value_init(acc);
676 for (unsigned long i = 0; i < sc.det; ++i) {
677 value_assign(acc, c->p[dim]);
678 zz2value(num.constant[i], x);
679 for (int j = dim-1; j >= 0; --j) {
680 value_multiply(acc, acc, x);
681 value_addto(acc, acc, c->p[j]);
682 }
683 value_addto(mpq_numref(count), mpq_numref(count), acc);
684 }
685 mpz_set(mpq_denref(count), c->p[dim+1]);
686 value_clear(acc);
687 value_clear(x);
688 }
689 evalue EV;
690 value_init(EV.d);
691 evalue_set(&EV, &count[0]._mp_num, &count[0]._mp_den);
692 eadd(&EV, vE[vert]);
693 free_evalue_refs(&EV);
694 }
695 }
696
697 struct ienumerator_base : enumerator_base {
698 evalue ** E_vertex;
699
700 ienumerator_base(unsigned dim, vertex_decomposer *vpd) :
701 enumerator_base(dim,vpd) {
702 E_vertex = new evalue_p[dim];
703 }
704
705 virtual ~ienumerator_base() {
706 delete [] E_vertex;
707 }
708
709 evalue *E_num(int i, int d) {
710 return E_vertex[i + (dim-d)];
711 }
712 };
713
714 struct cumulator {
715 evalue *factor;
716 evalue *v;
717 dpoly_r *r;
718
719 cumulator(evalue *factor, evalue *v, dpoly_r *r) :
720 factor(factor), v(v), r(r) {}
721
722 void cumulate(barvinok_options *options);
723
724 virtual void add_term(const vector<int>& powers, evalue *f2) = 0;
725 virtual ~cumulator() {}
726 };
727
728 void cumulator::cumulate(barvinok_options *options)
729 {
730 evalue cum; // factor * 1 * E_num[0]/1 * (E_num[0]-1)/2 *...
731 evalue f;
732 evalue t; // E_num[0] - (m-1)
733 evalue *cst;
734 evalue mone;
735
736 if (options->lookup_table) {
737 value_init(mone.d);
738 evalue_set_si(&mone, -1, 1);
739 }
740
741 value_init(cum.d);
742 evalue_copy(&cum, factor);
743 value_init(f.d);
744 value_init(f.x.n);
745 value_set_si(f.d, 1);
746 value_set_si(f.x.n, 1);
747 value_init(t.d);
748 evalue_copy(&t, v);
749
750 if (!options->lookup_table) {
751 for (cst = &t; value_zero_p(cst->d); ) {
752 if (cst->x.p->type == fractional)
753 cst = &cst->x.p->arr[1];
754 else
755 cst = &cst->x.p->arr[0];
756 }
757 }
758
759 for (int m = 0; m < r->len; ++m) {
760 if (m > 0) {
761 if (m > 1) {
762 value_set_si(f.d, m);
763 emul(&f, &cum);
764 if (!options->lookup_table)
765 value_subtract(cst->x.n, cst->x.n, cst->d);
766 else
767 eadd(&mone, &t);
768 }
769 emul(&t, &cum);
770 }
771 dpoly_r_term_list& current = r->c[r->len-1-m];
772 dpoly_r_term_list::iterator j;
773 for (j = current.begin(); j != current.end(); ++j) {
774 if ((*j)->coeff == 0)
775 continue;
776 evalue *f2 = new evalue;
777 value_init(f2->d);
778 value_init(f2->x.n);
779 zz2value((*j)->coeff, f2->x.n);
780 zz2value(r->denom, f2->d);
781 emul(&cum, f2);
782
783 add_term((*j)->powers, f2);
784 }
785 }
786 free_evalue_refs(&f);
787 free_evalue_refs(&t);
788 free_evalue_refs(&cum);
789 if (options->lookup_table)
790 free_evalue_refs(&mone);
791 }
792
793 struct E_poly_term {
794 vector<int> powers;
795 evalue *E;
796 };
797
798 struct ie_cum : public cumulator {
799 vector<E_poly_term *> terms;
800
801 ie_cum(evalue *factor, evalue *v, dpoly_r *r) : cumulator(factor, v, r) {}
802
803 virtual void add_term(const vector<int>& powers, evalue *f2);
804 };
805
806 void ie_cum::add_term(const vector<int>& powers, evalue *f2)
807 {
808 int k;
809 for (k = 0; k < terms.size(); ++k) {
810 if (terms[k]->powers == powers) {
811 eadd(f2, terms[k]->E);
812 free_evalue_refs(f2);
813 delete f2;
814 break;
815 }
816 }
817 if (k >= terms.size()) {
818 E_poly_term *ET = new E_poly_term;
819 ET->powers = powers;
820 ET->E = f2;
821 terms.push_back(ET);
822 }
823 }
824
825 struct ienumerator : public signed_cone_consumer, public vertex_decomposer,
826 public ienumerator_base {
827 //Polyhedron *pVD;
828 mat_ZZ den;
829 mat_ZZ vertex;
830 mpq_t tcount;
831 Value tz;
832
833 ienumerator(Polyhedron *P, unsigned dim, Param_Polyhedron *PP) :
834 vertex_decomposer(PP, *this), ienumerator_base(dim, this) {
835 vertex.SetDims(1, dim);
836
837 den.SetDims(dim, dim);
838 mpq_init(tcount);
839 value_init(tz);
840 }
841
842 ~ienumerator() {
843 mpq_clear(tcount);
844 value_clear(tz);
845 }
846
847 virtual void handle(const signed_cone& sc, barvinok_options *options);
848 void reduce(evalue *factor, const mat_ZZ& num, const mat_ZZ& den_f,
849 barvinok_options *options);
850 };
851
852 void ienumerator::reduce(evalue *factor, const mat_ZZ& num, const mat_ZZ& den_f,
853 barvinok_options *options)
854 {
855 unsigned len = den_f.NumRows(); // number of factors in den
856 unsigned dim = num.NumCols();
857 assert(num.NumRows() == 1);
858
859 if (dim == 0) {
860 eadd(factor, vE[vert]);
861 return;
862 }
863
864 vec_ZZ den_s;
865 mat_ZZ den_r;
866 vec_ZZ num_s;
867 mat_ZZ num_p;
868
869 split_one(num, num_s, num_p, den_f, den_s, den_r);
870
871 vec_ZZ den_p;
872 den_p.SetLength(len);
873
874 ZZ one;
875 one = 1;
876 normalize(one, num_s, num_p, den_s, den_p, den_r);
877 if (one != 1)
878 emul(&mone, factor);
879
880 int only_param = 0;
881 int no_param = 0;
882 for (int k = 0; k < len; ++k) {
883 if (den_p[k] == 0)
884 ++no_param;
885 else if (den_s[k] == 0)
886 ++only_param;
887 }
888 if (no_param == 0) {
889 reduce(factor, num_p, den_r, options);
890 } else {
891 int k, l;
892 mat_ZZ pden;
893 pden.SetDims(only_param, dim-1);
894
895 for (k = 0, l = 0; k < len; ++k)
896 if (den_s[k] == 0)
897 pden[l++] = den_r[k];
898
899 for (k = 0; k < len; ++k)
900 if (den_p[k] == 0)
901 break;
902
903 zz2value(num_s[0], tz);
904 dpoly n(no_param, tz);
905 zz2value(den_s[k], tz);
906 dpoly D(no_param, tz, 1);
907 for ( ; ++k < len; )
908 if (den_p[k] == 0) {
909 zz2value(den_s[k], tz);
910 dpoly fact(no_param, tz, 1);
911 D *= fact;
912 }
913
914 dpoly_r * r = 0;
915 // if no_param + only_param == len then all powers
916 // below will be all zero
917 if (no_param + only_param == len) {
918 if (E_num(0, dim) != 0)
919 r = new dpoly_r(n, len);
920 else {
921 mpq_set_si(tcount, 0, 1);
922 one = 1;
923 n.div(D, tcount, 1);
924
925 if (value_notzero_p(mpq_numref(tcount))) {
926 evalue f;
927 value_init(f.d);
928 value_init(f.x.n);
929 value_assign(f.x.n, mpq_numref(tcount));
930 value_assign(f.d, mpq_denref(tcount));
931 emul(&f, factor);
932 reduce(factor, num_p, pden, options);
933 free_evalue_refs(&f);
934 }
935 return;
936 }
937 } else {
938 for (k = 0; k < len; ++k) {
939 if (den_s[k] == 0 || den_p[k] == 0)
940 continue;
941
942 zz2value(den_s[k], tz);
943 dpoly pd(no_param-1, tz, 1);
944
945 int l;
946 for (l = 0; l < k; ++l)
947 if (den_r[l] == den_r[k])
948 break;
949
950 if (r == 0)
951 r = new dpoly_r(n, pd, l, len);
952 else {
953 dpoly_r *nr = new dpoly_r(r, pd, l, len);
954 delete r;
955 r = nr;
956 }
957 }
958 }
959 dpoly_r *rc = r->div(D);
960 delete r;
961 r = rc;
962 if (E_num(0, dim) == 0) {
963 int common = pden.NumRows();
964 dpoly_r_term_list& final = r->c[r->len-1];
965 int rows;
966 evalue t;
967 evalue f;
968 value_init(f.d);
969 value_init(f.x.n);
970 zz2value(r->denom, f.d);
971 dpoly_r_term_list::iterator j;
972 for (j = final.begin(); j != final.end(); ++j) {
973 if ((*j)->coeff == 0)
974 continue;
975 rows = common;
976 for (int k = 0; k < r->dim; ++k) {
977 int n = (*j)->powers[k];
978 if (n == 0)
979 continue;
980 pden.SetDims(rows+n, pden.NumCols());
981 for (int l = 0; l < n; ++l)
982 pden[rows+l] = den_r[k];
983 rows += n;
984 }
985 value_init(t.d);
986 evalue_copy(&t, factor);
987 zz2value((*j)->coeff, f.x.n);
988 emul(&f, &t);
989 reduce(&t, num_p, pden, options);
990 free_evalue_refs(&t);
991 }
992 free_evalue_refs(&f);
993 } else {
994 ie_cum cum(factor, E_num(0, dim), r);
995 cum.cumulate(options);
996
997 int common = pden.NumRows();
998 int rows;
999 for (int j = 0; j < cum.terms.size(); ++j) {
1000 rows = common;
1001 pden.SetDims(rows, pden.NumCols());
1002 for (int k = 0; k < r->dim; ++k) {
1003 int n = cum.terms[j]->powers[k];
1004 if (n == 0)
1005 continue;
1006 pden.SetDims(rows+n, pden.NumCols());
1007 for (int l = 0; l < n; ++l)
1008 pden[rows+l] = den_r[k];
1009 rows += n;
1010 }
1011 reduce(cum.terms[j]->E, num_p, pden, options);
1012 free_evalue_refs(cum.terms[j]->E);
1013 delete cum.terms[j]->E;
1014 delete cum.terms[j];
1015 }
1016 }
1017 delete r;
1018 }
1019 }
1020
1021 static int type_offset(enode *p)
1022 {
1023 return p->type == fractional ? 1 :
1024 p->type == flooring ? 1 : 0;
1025 }
1026
1027 static int edegree(evalue *e)
1028 {
1029 int d = 0;
1030 enode *p;
1031
1032 if (value_notzero_p(e->d))
1033 return 0;
1034
1035 p = e->x.p;
1036 int i = type_offset(p);
1037 if (p->size-i-1 > d)
1038 d = p->size - i - 1;
1039 for (; i < p->size; i++) {
1040 int d2 = edegree(&p->arr[i]);
1041 if (d2 > d)
1042 d = d2;
1043 }
1044 return d;
1045 }
1046
1047 void ienumerator::handle(const signed_cone& sc, barvinok_options *options)
1048 {
1049 assert(sc.det == 1);
1050 assert(sc.rays.NumRows() == dim);
1051
1052 lattice_point(V, sc.rays, vertex[0], E_vertex, options);
1053
1054 den = sc.rays;
1055
1056 evalue one;
1057 value_init(one.d);
1058 evalue_set_si(&one, sc.sign, 1);
1059 reduce(&one, vertex, den, options);
1060 free_evalue_refs(&one);
1061
1062 for (int i = 0; i < dim; ++i)
1063 if (E_vertex[i]) {
1064 free_evalue_refs(E_vertex[i]);
1065 delete E_vertex[i];
1066 }
1067 }
1068
1069 struct bfenumerator : public vertex_decomposer, public bf_base,
1070 public ienumerator_base {
1071 evalue *factor;
1072
1073 bfenumerator(Polyhedron *P, unsigned dim, Param_Polyhedron *PP) :
1074 vertex_decomposer(PP, *this),
1075 bf_base(dim), ienumerator_base(dim, this) {
1076 lower = 0;
1077 factor = NULL;
1078 }
1079
1080 ~bfenumerator() {
1081 }
1082
1083 virtual void handle(const signed_cone& sc, barvinok_options *options);
1084 virtual void base(mat_ZZ& factors, bfc_vec& v);
1085
1086 bfc_term_base* new_bf_term(int len) {
1087 bfe_term* t = new bfe_term(len);
1088 return t;
1089 }
1090
1091 virtual void set_factor(bfc_term_base *t, int k, int change) {
1092 bfe_term* bfet = static_cast<bfe_term *>(t);
1093 factor = bfet->factors[k];
1094 assert(factor != NULL);
1095 bfet->factors[k] = NULL;
1096 if (change)
1097 emul(&mone, factor);
1098 }
1099
1100 virtual void set_factor(bfc_term_base *t, int k, mpq_t &q, int change) {
1101 bfe_term* bfet = static_cast<bfe_term *>(t);
1102 factor = bfet->factors[k];
1103 assert(factor != NULL);
1104 bfet->factors[k] = NULL;
1105
1106 evalue f;
1107 value_init(f.d);
1108 value_init(f.x.n);
1109 if (change)
1110 value_oppose(f.x.n, mpq_numref(q));
1111 else
1112 value_assign(f.x.n, mpq_numref(q));
1113 value_assign(f.d, mpq_denref(q));
1114 emul(&f, factor);
1115 free_evalue_refs(&f);
1116 }
1117
1118 virtual void set_factor(bfc_term_base *t, int k, const QQ& c, int change) {
1119 bfe_term* bfet = static_cast<bfe_term *>(t);
1120
1121 factor = new evalue;
1122
1123 evalue f;
1124 value_init(f.d);
1125 value_init(f.x.n);
1126 zz2value(c.n, f.x.n);
1127 if (change)
1128 value_oppose(f.x.n, f.x.n);
1129 zz2value(c.d, f.d);
1130
1131 value_init(factor->d);
1132 evalue_copy(factor, bfet->factors[k]);
1133 emul(&f, factor);
1134 free_evalue_refs(&f);
1135 }
1136
1137 void set_factor(evalue *f, int change) {
1138 if (change)
1139 emul(&mone, f);
1140 factor = f;
1141 }
1142
1143 virtual void insert_term(bfc_term_base *t, int i) {
1144 bfe_term* bfet = static_cast<bfe_term *>(t);
1145 int len = t->terms.NumRows()-1; // already increased by one
1146
1147 bfet->factors.resize(len+1);
1148 for (int j = len; j > i; --j) {
1149 bfet->factors[j] = bfet->factors[j-1];
1150 t->terms[j] = t->terms[j-1];
1151 }
1152 bfet->factors[i] = factor;
1153 factor = NULL;
1154 }
1155
1156 virtual void update_term(bfc_term_base *t, int i) {
1157 bfe_term* bfet = static_cast<bfe_term *>(t);
1158
1159 eadd(factor, bfet->factors[i]);
1160 free_evalue_refs(factor);
1161 delete factor;
1162 }
1163
1164 virtual bool constant_vertex(int dim) { return E_num(0, dim) == 0; }
1165
1166 virtual void cum(bf_reducer *bfr, bfc_term_base *t, int k, dpoly_r *r,
1167 barvinok_options *options);
1168 };
1169
1170 enumerator_base *enumerator_base::create(Polyhedron *P, unsigned dim,
1171 Param_Polyhedron *PP,
1172 barvinok_options *options)
1173 {
1174 enumerator_base *eb;
1175
1176 if (options->incremental_specialization == BV_SPECIALIZATION_BF)
1177 eb = new bfenumerator(P, dim, PP);
1178 else if (options->incremental_specialization == BV_SPECIALIZATION_DF)
1179 eb = new ienumerator(P, dim, PP);
1180 else
1181 eb = new enumerator(P, dim, PP);
1182
1183 return eb;
1184 }
1185
1186 struct bfe_cum : public cumulator {
1187 bfenumerator *bfe;
1188 bfc_term_base *told;
1189 int k;
1190 bf_reducer *bfr;
1191
1192 bfe_cum(evalue *factor, evalue *v, dpoly_r *r, bf_reducer *bfr,
1193 bfc_term_base *t, int k, bfenumerator *e) :
1194 cumulator(factor, v, r), told(t), k(k),
1195 bfr(bfr), bfe(e) {
1196 }
1197
1198 virtual void add_term(const vector<int>& powers, evalue *f2);
1199 };
1200
1201 void bfe_cum::add_term(const vector<int>& powers, evalue *f2)
1202 {
1203 bfr->update_powers(powers);
1204
1205 bfc_term_base * t = bfe->find_bfc_term(bfr->vn, bfr->npowers, bfr->nnf);
1206 bfe->set_factor(f2, bfr->l_changes % 2);
1207 bfe->add_term(t, told->terms[k], bfr->l_extra_num);
1208 }
1209
1210 void bfenumerator::cum(bf_reducer *bfr, bfc_term_base *t, int k,
1211 dpoly_r *r, barvinok_options *options)
1212 {
1213 bfe_term* bfet = static_cast<bfe_term *>(t);
1214 bfe_cum cum(bfet->factors[k], E_num(0, bfr->d), r, bfr, t, k, this);
1215 cum.cumulate(options);
1216 }
1217
1218 void bfenumerator::base(mat_ZZ& factors, bfc_vec& v)
1219 {
1220 for (int i = 0; i < v.size(); ++i) {
1221 assert(v[i]->terms.NumRows() == 1);
1222 evalue *factor = static_cast<bfe_term *>(v[i])->factors[0];
1223 eadd(factor, vE[vert]);
1224 delete v[i];
1225 }
1226 }
1227
1228 void bfenumerator::handle(const signed_cone& sc, barvinok_options *options)
1229 {
1230 assert(sc.det == 1);
1231 assert(sc.rays.NumRows() == enumerator_base::dim);
1232
1233 bfe_term* t = new bfe_term(enumerator_base::dim);
1234 vector< bfc_term_base * > v;
1235 v.push_back(t);
1236
1237 t->factors.resize(1);
1238
1239 t->terms.SetDims(1, enumerator_base::dim);
1240 lattice_point(V, sc.rays, t->terms[0], E_vertex, options);
1241
1242 // the elements of factors are always lexpositive
1243 mat_ZZ factors;
1244 int s = setup_factors(sc.rays, factors, t, sc.sign);
1245
1246 t->factors[0] = new evalue;
1247 value_init(t->factors[0]->d);
1248 evalue_set_si(t->factors[0], s, 1);
1249 reduce(factors, v, options);
1250
1251 for (int i = 0; i < enumerator_base::dim; ++i)
1252 if (E_vertex[i]) {
1253 free_evalue_refs(E_vertex[i]);
1254 delete E_vertex[i];
1255 }
1256 }
1257
1258 static evalue* barvinok_enumerate_ev_f(Polyhedron *P, Polyhedron* C,
1259 barvinok_options *options);
1260
1261 /* Destroys C */
1262 static evalue* barvinok_enumerate_cst(Polyhedron *P, Polyhedron* C,
1263 struct barvinok_options *options)
1264 {
1265 evalue *eres;
1266
1267 if (emptyQ2(C)) {
1268 Polyhedron_Free(C);
1269 return evalue_zero();
1270 }
1271
1272 ALLOC(evalue, eres);
1273 value_init(eres->d);
1274 value_set_si(eres->d, 0);
1275 eres->x.p = new_enode(partition, 2, C->Dimension);
1276 EVALUE_SET_DOMAIN(eres->x.p->arr[0],
1277 DomainConstraintSimplify(C, options->MaxRays));
1278 value_set_si(eres->x.p->arr[1].d, 1);
1279 value_init(eres->x.p->arr[1].x.n);
1280 if (emptyQ2(P))
1281 value_set_si(eres->x.p->arr[1].x.n, 0);
1282 else
1283 barvinok_count_with_options(P, &eres->x.p->arr[1].x.n, options);
1284
1285 return eres;
1286 }
1287
1288 static evalue* enumerate(Polyhedron *P, Polyhedron* C,
1289 struct barvinok_options *options)
1290 {
1291 Polyhedron *next;
1292 Polyhedron *Porig = P;
1293 Polyhedron *Corig = C;
1294 Polyhedron *CEq = NULL, *rVD;
1295 int r = 0;
1296 unsigned nparam = C->Dimension;
1297 evalue *eres;
1298 Matrix *CP = NULL;
1299
1300 evalue factor;
1301 value_init(factor.d);
1302 evalue_set_si(&factor, 1, 1);
1303
1304 /* for now */
1305 POL_ENSURE_FACETS(P);
1306 POL_ENSURE_VERTICES(P);
1307 POL_ENSURE_FACETS(C);
1308 POL_ENSURE_VERTICES(C);
1309
1310 if (C->Dimension == 0 || emptyQ(P) || emptyQ(C)) {
1311 constant:
1312 if (CEq == Porig)
1313 CEq = Polyhedron_Copy(CEq);
1314 eres = barvinok_enumerate_cst(P, CEq ? CEq : Polyhedron_Copy(C), options);
1315 out:
1316 if (CP) {
1317 evalue_backsubstitute(eres, CP, options->MaxRays);
1318 Matrix_Free(CP);
1319 }
1320
1321 emul(&factor, eres);
1322 if (options->approximation_method == BV_APPROX_DROP) {
1323 if (options->polynomial_approximation == BV_APPROX_SIGN_UPPER)
1324 evalue_frac2polynomial(eres, 1, options->MaxRays);
1325 if (options->polynomial_approximation == BV_APPROX_SIGN_LOWER)
1326 evalue_frac2polynomial(eres, -1, options->MaxRays);
1327 if (options->polynomial_approximation == BV_APPROX_SIGN_APPROX)
1328 evalue_frac2polynomial(eres, 0, options->MaxRays);
1329 }
1330 reduce_evalue(eres);
1331 free_evalue_refs(&factor);
1332 if (P != Porig)
1333 Domain_Free(P);
1334 if (C != Corig)
1335 Polyhedron_Free(C);
1336
1337 return eres;
1338 }
1339 if (Polyhedron_is_unbounded(P, nparam, options->MaxRays))
1340 goto constant;
1341
1342 if (P->Dimension == nparam) {
1343 CEq = P;
1344 P = Universe_Polyhedron(0);
1345 goto constant;
1346 }
1347 if (P->NbEq != 0 || C->NbEq != 0) {
1348 Polyhedron *Q = P;
1349 Polyhedron *D = C;
1350 remove_all_equalities(&P, &C, &CP, NULL, nparam, options->MaxRays);
1351 if (C != D && D != Corig)
1352 Polyhedron_Free(D);
1353 if (P != Q && Q != Porig)
1354 Domain_Free(Q);
1355 eres = enumerate(P, C, options);
1356 goto out;
1357 }
1358
1359 Polyhedron *T = Polyhedron_Factor(P, nparam, NULL, options->MaxRays);
1360 if (T || (P->Dimension == nparam+1)) {
1361 Polyhedron *Q;
1362 Polyhedron *C2;
1363 for (Q = T ? T : P; Q; Q = Q->next) {
1364 Polyhedron *next = Q->next;
1365 Q->next = NULL;
1366
1367 Polyhedron *QC = Q;
1368 if (Q->Dimension != C->Dimension)
1369 QC = Polyhedron_Project(Q, nparam);
1370
1371 C2 = C;
1372 C = DomainIntersection(C, QC, options->MaxRays);
1373 if (C2 != Corig)
1374 Polyhedron_Free(C2);
1375 if (QC != Q)
1376 Polyhedron_Free(QC);
1377
1378 Q->next = next;
1379 }
1380 }
1381 if (T) {
1382 if (P != Porig)
1383 Polyhedron_Free(P);
1384 P = T;
1385 if (T->Dimension == C->Dimension) {
1386 P = T->next;
1387 T->next = NULL;
1388 Polyhedron_Free(T);
1389 }
1390 }
1391
1392 next = P->next;
1393 P->next = NULL;
1394 eres = barvinok_enumerate_ev_f(P, C, options);
1395 P->next = next;
1396
1397 if (P->next) {
1398 Polyhedron *Q;
1399 evalue *f;
1400
1401 for (Q = P->next; Q; Q = Q->next) {
1402 Polyhedron *next = Q->next;
1403 Q->next = NULL;
1404
1405 f = barvinok_enumerate_ev_f(Q, C, options);
1406 emul(f, eres);
1407 evalue_free(f);
1408
1409 Q->next = next;
1410 }
1411 }
1412
1413 goto out;
1414 }
1415
1416 evalue* barvinok_enumerate_with_options(Polyhedron *P, Polyhedron* C,
1417 struct barvinok_options *options)
1418 {
1419 Polyhedron *next, *Cnext, *C1;
1420 Polyhedron *Corig = C;
1421 evalue *eres;
1422
1423 if (P->next)
1424 fprintf(stderr,
1425 "barvinok_enumerate: input is a union; only first polyhedron is enumerated\n");
1426
1427 if (C->next)
1428 fprintf(stderr,
1429 "barvinok_enumerate: context is a union; only first polyhedron is considered\n");
1430
1431 Cnext = C->next;
1432 C->next = NULL;
1433 C1 = Polyhedron_Project(P, C->Dimension);
1434 C = DomainIntersection(C, C1, options->MaxRays);
1435 Polyhedron_Free(C1);
1436 next = P->next;
1437 P->next = NULL;
1438
1439 if (options->approximation_method == BV_APPROX_BERNOULLI ||
1440 options->summation == BV_SUM_BERNOULLI)
1441 eres = Bernoulli_sum(P, C, options);
1442 else
1443 eres = enumerate(P, C, options);
1444 Domain_Free(C);
1445
1446 P->next= next;
1447 Corig->next = Cnext;
1448
1449 return eres;
1450 }
1451
1452 evalue* barvinok_enumerate_ev(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
1453 {
1454 evalue *E;
1455 barvinok_options *options = barvinok_options_new_with_defaults();
1456 options->MaxRays = MaxRays;
1457 E = barvinok_enumerate_with_options(P, C, options);
1458 barvinok_options_free(options);
1459 return E;
1460 }
1461
1462 evalue *Param_Polyhedron_Enumerate(Param_Polyhedron *PP, Polyhedron *P,
1463 Polyhedron *C,
1464 struct barvinok_options *options)
1465 {
1466 evalue *eres;
1467 Param_Domain *D;
1468 unsigned nparam = C->Dimension;
1469 unsigned dim = P->Dimension - nparam;
1470
1471 int nd;
1472 for (nd = 0, D=PP->D; D; ++nd, D=D->next);
1473 evalue_section *s = new evalue_section[nd];
1474
1475 enumerator_base *et = NULL;
1476 try_again:
1477 if (et)
1478 delete et;
1479
1480 et = enumerator_base::create(P, dim, PP, options);
1481
1482 Polyhedron *TC = true_context(P, C, options->MaxRays);
1483 FORALL_REDUCED_DOMAIN(PP, TC, nd, options, i, D, rVD)
1484 Param_Vertices *V;
1485
1486 s[i].E = evalue_zero();
1487 s[i].D = rVD;
1488
1489 FORALL_PVertex_in_ParamPolyhedron(V,D,PP) // _i is internal counter
1490 if (!et->vE[_i])
1491 try {
1492 et->decompose_at(V, _i, options);
1493 } catch (OrthogonalException &e) {
1494 FORALL_REDUCED_DOMAIN_RESET;
1495 for (; i >= 0; --i) {
1496 evalue_free(s[i].E);
1497 Domain_Free(s[i].D);
1498 }
1499 goto try_again;
1500 }
1501 eadd(et->vE[_i] , s[i].E);
1502 END_FORALL_PVertex_in_ParamPolyhedron;
1503 evalue_range_reduction_in_domain(s[i].E, rVD);
1504 END_FORALL_REDUCED_DOMAIN
1505 Polyhedron_Free(TC);
1506
1507 delete et;
1508 eres = evalue_from_section_array(s, nd);
1509 delete [] s;
1510
1511 return eres;
1512 }
1513
1514 static evalue* barvinok_enumerate_ev_f(Polyhedron *P, Polyhedron* C,
1515 barvinok_options *options)
1516 {
1517 unsigned nparam = C->Dimension;
1518 bool do_scale = options->approximation_method == BV_APPROX_SCALE;
1519
1520 if (options->approximation_method == BV_APPROX_VOLUME)
1521 return Param_Polyhedron_Volume(P, C, options);
1522
1523 if (P->Dimension - nparam == 1 && !do_scale)
1524 return ParamLine_Length(P, C, options);
1525
1526 Param_Polyhedron *PP = NULL;
1527 evalue *eres;
1528
1529 if (do_scale) {
1530 eres = scale_bound(P, C, options);
1531 if (eres)
1532 return eres;
1533 }
1534
1535 PP = Polyhedron2Param_Polyhedron(P, C, options);
1536
1537 if (do_scale)
1538 eres = scale(PP, P, C, options);
1539 else
1540 eres = Param_Polyhedron_Enumerate(PP, P, C, options);
1541
1542 if (PP)
1543 Param_Polyhedron_Free(PP);
1544
1545 return eres;
1546 }
1547
1548 Enumeration* barvinok_enumerate(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
1549 {
1550 evalue *EP = barvinok_enumerate_ev(P, C, MaxRays);
1551
1552 return partition2enumeration(EP);
1553 }
1554
1555 evalue* barvinok_enumerate_union(Polyhedron *D, Polyhedron* C, unsigned MaxRays)
1556 {
1557 evalue *EP;
1558 gen_fun *gf = barvinok_enumerate_union_series(D, C, MaxRays);
1559 EP = *gf;
1560 delete gf;
1561 return EP;
1562 }
1563
1564 evalue *barvinok_summate(evalue *e, int nvar, struct barvinok_options *options)
1565 {
1566 if (options->summation == BV_SUM_EULER)
1567 return euler_summate(e, nvar, options);
1568 else if (options->summation == BV_SUM_BERNOULLI)
1569 return Bernoulli_sum_evalue(e, nvar, options);
1570 else
1571 return evalue_sum(e, nvar, options->MaxRays);
1572 }
lattice point enumeration library