?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
doc: update documentation of options and some functions
[barvinok.git] / barvinok.cc
blobcfad4c4906598a7340a6a5af15d2011230b3ab19
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 extern "C" {
13 #include "piputil.h"
14 }
15 #include "config.h"
16 #include <barvinok/barvinok.h>
17 #include <barvinok/genfun.h>
18 #include <barvinok/options.h>
19 #include <barvinok/sample.h>
20 #include "conversion.h"
21 #include "counter.h"
22 #include "tcounter.h"
23 #include "decomposer.h"
24 #include "lattice_point.h"
25 #include "reduce_domain.h"
26 #include "genfun_constructor.h"
27 #include "remove_equalities.h"
28 #include "scale.h"
29 #include "volume.h"
30 #include "bernoulli.h"
31
32 #ifdef NTL_STD_CXX
33 using namespace NTL;
34 #endif
35 using std::cerr;
36 using std::cout;
37 using std::endl;
38 using std::vector;
39 using std::deque;
40 using std::string;
41 using std::ostringstream;
42
43 #define ALLOC(t,p) p = (t*)malloc(sizeof(*p))
44
45 class dpoly_n {
46 public:
47 Matrix *coeff;
48 ~dpoly_n() {
49 Matrix_Free(coeff);
50 }
51 dpoly_n(int d) {
52 Value d0, one;
53 value_init(d0);
54 value_init(one);
55 value_set_si(one, 1);
56 coeff = Matrix_Alloc(d+1, d+1+1);
57 value_set_si(coeff->p[0][0], 1);
58 value_set_si(coeff->p[0][d+1], 1);
59 for (int i = 1; i <= d; ++i) {
60 value_multiply(coeff->p[i][0], coeff->p[i-1][0], d0);
61 Vector_Combine(coeff->p[i-1], coeff->p[i-1]+1, coeff->p[i]+1,
62 one, d0, i);
63 value_set_si(coeff->p[i][d+1], i);
64 value_multiply(coeff->p[i][d+1], coeff->p[i][d+1], coeff->p[i-1][d+1]);
65 value_decrement(d0, d0);
66 }
67 value_clear(d0);
68 value_clear(one);
69 }
70 void div(dpoly& d, Vector *count, ZZ& sign) {
71 int len = coeff->NbRows;
72 Matrix * c = Matrix_Alloc(coeff->NbRows, coeff->NbColumns);
73 Value tmp;
74 value_init(tmp);
75 for (int i = 0; i < len; ++i) {
76 Vector_Copy(coeff->p[i], c->p[i], len+1);
77 for (int j = 1; j <= i; ++j) {
78 value_multiply(tmp, d.coeff->p[j], c->p[i][len]);
79 value_oppose(tmp, tmp);
80 Vector_Combine(c->p[i], c->p[i-j], c->p[i],
81 c->p[i-j][len], tmp, len);
82 value_multiply(c->p[i][len], c->p[i][len], c->p[i-j][len]);
83 }
84 value_multiply(c->p[i][len], c->p[i][len], d.coeff->p[0]);
85 }
86 if (sign == -1) {
87 value_set_si(tmp, -1);
88 Vector_Scale(c->p[len-1], count->p, tmp, len);
89 value_assign(count->p[len], c->p[len-1][len]);
90 } else
91 Vector_Copy(c->p[len-1], count->p, len+1);
92 Vector_Normalize(count->p, len+1);
93 value_clear(tmp);
94 Matrix_Free(c);
95 }
96 };
97
98 const int MAX_TRY=10;
99 /*
100 * Searches for a vector that is not orthogonal to any
101 * of the rays in rays.
102 */
103 static void nonorthog(mat_ZZ& rays, vec_ZZ& lambda)
104 {
105 int dim = rays.NumCols();
106 bool found = false;
107 lambda.SetLength(dim);
108 if (dim == 0)
109 return;
110
111 for (int i = 2; !found && i <= 50*dim; i+=4) {
112 for (int j = 0; j < MAX_TRY; ++j) {
113 for (int k = 0; k < dim; ++k) {
114 int r = random_int(i)+2;
115 int v = (2*(r%2)-1) * (r >> 1);
116 lambda[k] = v;
117 }
118 int k = 0;
119 for (; k < rays.NumRows(); ++k)
120 if (lambda * rays[k] == 0)
121 break;
122 if (k == rays.NumRows()) {
123 found = true;
124 break;
125 }
126 }
127 }
128 assert(found);
129 }
130
131 static void add_rays(mat_ZZ& rays, Polyhedron *i, int *r, int nvar = -1,
132 bool all = false)
133 {
134 unsigned dim = i->Dimension;
135 if (nvar == -1)
136 nvar = dim;
137 for (int k = 0; k < i->NbRays; ++k) {
138 if (!value_zero_p(i->Ray[k][dim+1]))
139 continue;
140 if (!all && nvar != dim && First_Non_Zero(i->Ray[k]+1, nvar) == -1)
141 continue;
142 values2zz(i->Ray[k]+1, rays[(*r)++], nvar);
143 }
144 }
145
146 static void mask_r(Matrix *f, int nr, Vector *lcm, int p, Vector *val, evalue *ev)
147 {
148 unsigned nparam = lcm->Size;
149
150 if (p == nparam) {
151 Vector * prod = Vector_Alloc(f->NbRows);
152 Matrix_Vector_Product(f, val->p, prod->p);
153 int isint = 1;
154 for (int i = 0; i < nr; ++i) {
155 value_modulus(prod->p[i], prod->p[i], f->p[i][nparam+1]);
156 isint &= value_zero_p(prod->p[i]);
157 }
158 value_set_si(ev->d, 1);
159 value_init(ev->x.n);
160 value_set_si(ev->x.n, isint);
161 Vector_Free(prod);
162 return;
163 }
164
165 Value tmp;
166 value_init(tmp);
167 if (value_one_p(lcm->p[p]))
168 mask_r(f, nr, lcm, p+1, val, ev);
169 else {
170 value_assign(tmp, lcm->p[p]);
171 value_set_si(ev->d, 0);
172 ev->x.p = new_enode(periodic, VALUE_TO_INT(tmp), p+1);
173 do {
174 value_decrement(tmp, tmp);
175 value_assign(val->p[p], tmp);
176 mask_r(f, nr, lcm, p+1, val, &ev->x.p->arr[VALUE_TO_INT(tmp)]);
177 } while (value_pos_p(tmp));
178 }
179 value_clear(tmp);
180 }
181
182 static void mask_fractional(Matrix *f, evalue *factor)
183 {
184 int nr = f->NbRows, nc = f->NbColumns;
185 int n;
186 bool found = false;
187 for (n = 0; n < nr && value_notzero_p(f->p[n][nc-1]); ++n)
188 if (value_notone_p(f->p[n][nc-1]) &&
189 value_notmone_p(f->p[n][nc-1]))
190 found = true;
191 if (!found)
192 return;
193
194 evalue EP;
195 nr = n;
196
197 Value m;
198 value_init(m);
199
200 evalue EV;
201 value_init(EV.d);
202 value_init(EV.x.n);
203 value_set_si(EV.x.n, 1);
204
205 for (n = 0; n < nr; ++n) {
206 value_assign(m, f->p[n][nc-1]);
207 if (value_one_p(m) || value_mone_p(m))
208 continue;
209
210 int j = normal_mod(f->p[n], nc-1, &m);
211 if (j == nc-1) {
212 free_evalue_refs(factor);
213 value_init(factor->d);
214 evalue_set_si(factor, 0, 1);
215 break;
216 }
217 vec_ZZ row;
218 values2zz(f->p[n], row, nc-1);
219 ZZ g;
220 value2zz(m, g);
221 if (j < (nc-1)-1 && row[j] > g/2) {
222 for (int k = j; k < (nc-1); ++k)
223 if (row[k] != 0)
224 row[k] = g - row[k];
225 }
226
227 value_init(EP.d);
228 value_set_si(EP.d, 0);
229 EP.x.p = new_enode(relation, 2, 0);
230 value_clear(EP.x.p->arr[1].d);
231 EP.x.p->arr[1] = *factor;
232 evalue *ev = &EP.x.p->arr[0];
233 value_set_si(ev->d, 0);
234 ev->x.p = new_enode(fractional, 3, -1);
235 evalue_set_si(&ev->x.p->arr[1], 0, 1);
236 evalue_set_si(&ev->x.p->arr[2], 1, 1);
237 evalue *E = multi_monom(row);
238 value_assign(EV.d, m);
239 emul(&EV, E);
240 value_clear(ev->x.p->arr[0].d);
241 ev->x.p->arr[0] = *E;
242 delete E;
243 *factor = EP;
244 }
245
246 value_clear(m);
247 free_evalue_refs(&EV);
248 }
249
250 /*
251 *
252 */
253 static void mask_table(Matrix *f, evalue *factor)
254 {
255 int nr = f->NbRows, nc = f->NbColumns;
256 int n;
257 bool found = false;
258 for (n = 0; n < nr && value_notzero_p(f->p[n][nc-1]); ++n)
259 if (value_notone_p(f->p[n][nc-1]) &&
260 value_notmone_p(f->p[n][nc-1]))
261 found = true;
262 if (!found)
263 return;
264
265 Value tmp;
266 value_init(tmp);
267 nr = n;
268 unsigned np = nc - 2;
269 Vector *lcm = Vector_Alloc(np);
270 Vector *val = Vector_Alloc(nc);
271 Vector_Set(val->p, 0, nc);
272 value_set_si(val->p[np], 1);
273 Vector_Set(lcm->p, 1, np);
274 for (n = 0; n < nr; ++n) {
275 if (value_one_p(f->p[n][nc-1]) ||
276 value_mone_p(f->p[n][nc-1]))
277 continue;
278 for (int j = 0; j < np; ++j)
279 if (value_notzero_p(f->p[n][j])) {
280 Gcd(f->p[n][j], f->p[n][nc-1], &tmp);
281 value_division(tmp, f->p[n][nc-1], tmp);
282 value_lcm(tmp, lcm->p[j], &lcm->p[j]);
283 }
284 }
285 evalue EP;
286 value_init(EP.d);
287 mask_r(f, nr, lcm, 0, val, &EP);
288 value_clear(tmp);
289 Vector_Free(val);
290 Vector_Free(lcm);
291 emul(&EP,factor);
292 free_evalue_refs(&EP);
293 }
294
295 static void mask(Matrix *f, evalue *factor, barvinok_options *options)
296 {
297 if (options->lookup_table)
298 mask_table(f, factor);
299 else
300 mask_fractional(f, factor);
301 }
302
303 struct bfe_term : public bfc_term_base {
304 vector<evalue *> factors;
305
306 bfe_term(int len) : bfc_term_base(len) {
307 }
308
309 ~bfe_term() {
310 for (int i = 0; i < factors.size(); ++i) {
311 if (!factors[i])
312 continue;
313 free_evalue_refs(factors[i]);
314 delete factors[i];
315 }
316 }
317 };
318
319 static void print_int_vector(int *v, int len, char *name)
320 {
321 cerr << name << endl;
322 for (int j = 0; j < len; ++j) {
323 cerr << v[j] << " ";
324 }
325 cerr << endl;
326 }
327
328 static void print_bfc_terms(mat_ZZ& factors, bfc_vec& v)
329 {
330 cerr << endl;
331 cerr << "factors" << endl;
332 cerr << factors << endl;
333 for (int i = 0; i < v.size(); ++i) {
334 cerr << "term: " << i << endl;
335 print_int_vector(v[i]->powers, factors.NumRows(), "powers");
336 cerr << "terms" << endl;
337 cerr << v[i]->terms << endl;
338 bfc_term* bfct = static_cast<bfc_term *>(v[i]);
339 cerr << bfct->c << endl;
340 }
341 }
342
343 static void print_bfe_terms(mat_ZZ& factors, bfc_vec& v)
344 {
345 cerr << endl;
346 cerr << "factors" << endl;
347 cerr << factors << endl;
348 for (int i = 0; i < v.size(); ++i) {
349 cerr << "term: " << i << endl;
350 print_int_vector(v[i]->powers, factors.NumRows(), "powers");
351 cerr << "terms" << endl;
352 cerr << v[i]->terms << endl;
353 bfe_term* bfet = static_cast<bfe_term *>(v[i]);
354 for (int j = 0; j < v[i]->terms.NumRows(); ++j) {
355 char * test[] = {"a", "b"};
356 print_evalue(stderr, bfet->factors[j], test);
357 fprintf(stderr, "\n");
358 }
359 }
360 }
361
362 struct bfcounter : public bfcounter_base {
363 mpq_t count;
364 Value tz;
365
366 bfcounter(unsigned dim) : bfcounter_base(dim) {
367 mpq_init(count);
368 lower = 1;
369 value_init(tz);
370 }
371 ~bfcounter() {
372 mpq_clear(count);
373 value_clear(tz);
374 }
375 virtual void base(mat_ZZ& factors, bfc_vec& v);
376 virtual void get_count(Value *result) {
377 assert(value_one_p(&count[0]._mp_den));
378 value_assign(*result, &count[0]._mp_num);
379 }
380 };
381
382 void bfcounter::base(mat_ZZ& factors, bfc_vec& v)
383 {
384 unsigned nf = factors.NumRows();
385
386 for (int i = 0; i < v.size(); ++i) {
387 bfc_term* bfct = static_cast<bfc_term *>(v[i]);
388 int total_power = 0;
389 // factor is always positive, so we always
390 // change signs
391 for (int k = 0; k < nf; ++k)
392 total_power += v[i]->powers[k];
393
394 int j;
395 for (j = 0; j < nf; ++j)
396 if (v[i]->powers[j] > 0)
397 break;
398
399 zz2value(factors[j][0], tz);
400 dpoly D(total_power, tz, 1);
401 for (int k = 1; k < v[i]->powers[j]; ++k) {
402 zz2value(factors[j][0], tz);
403 dpoly fact(total_power, tz, 1);
404 D *= fact;
405 }
406 for ( ; ++j < nf; )
407 for (int k = 0; k < v[i]->powers[j]; ++k) {
408 zz2value(factors[j][0], tz);
409 dpoly fact(total_power, tz, 1);
410 D *= fact;
411 }
412
413 for (int k = 0; k < v[i]->terms.NumRows(); ++k) {
414 zz2value(v[i]->terms[k][0], tz);
415 dpoly n(total_power, tz);
416 mpq_set_si(tcount, 0, 1);
417 n.div(D, tcount, one);
418 if (total_power % 2)
419 bfct->c[k].n = -bfct->c[k].n;
420 zz2value(bfct->c[k].n, tn);
421 zz2value(bfct->c[k].d, td);
422
423 mpz_mul(mpq_numref(tcount), mpq_numref(tcount), tn);
424 mpz_mul(mpq_denref(tcount), mpq_denref(tcount), td);
425 mpq_canonicalize(tcount);
426 mpq_add(count, count, tcount);
427 }
428 delete v[i];
429 }
430 }
431
432
433 /* Check whether the polyhedron is unbounded and if so,
434 * check whether it has any (and therefore an infinite number of)
435 * integer points.
436 * If one of the vertices is integer, then we are done.
437 * Otherwise, transform the polyhedron such that one of the rays
438 * is the first unit vector and cut it off at a height that ensures
439 * that if the whole polyhedron has any points, then the remaining part
440 * has integer points. In particular we add the largest coefficient
441 * of a ray to the highest vertex (rounded up).
442 */
443 static bool Polyhedron_is_infinite(Polyhedron *P, Value* result,
444 barvinok_options *options)
445 {
446 int r = 0;
447 Matrix *M, *M2;
448 Value c, tmp;
449 Value g;
450 bool first;
451 Vector *v;
452 Value offset, size;
453 Polyhedron *R;
454
455 if (P->NbBid == 0)
456 for (; r < P->NbRays; ++r)
457 if (value_zero_p(P->Ray[r][P->Dimension+1]))
458 break;
459 if (P->NbBid == 0 && r == P->NbRays)
460 return false;
461
462 if (options->count_sample_infinite) {
463 Vector *sample;
464
465 sample = Polyhedron_Sample(P, options);
466 if (!sample)
467 value_set_si(*result, 0);
468 else {
469 value_set_si(*result, -1);
470 Vector_Free(sample);
471 }
472 return true;
473 }
474
475 for (int i = 0; i < P->NbRays; ++i)
476 if (value_one_p(P->Ray[i][1+P->Dimension])) {
477 value_set_si(*result, -1);
478 return true;
479 }
480
481 value_init(g);
482 M = Matrix_Alloc(P->Dimension+1, P->Dimension+1);
483 Vector_Gcd(P->Ray[r]+1, P->Dimension, &g);
484 Vector_AntiScale(P->Ray[r]+1, M->p[0], g, P->Dimension+1);
485 int ok = unimodular_complete(M, 1);
486 assert(ok);
487 value_set_si(M->p[P->Dimension][P->Dimension], 1);
488 M2 = Transpose(M);
489 Matrix_Free(M);
490 P = Polyhedron_Preimage(P, M2, 0);
491 Matrix_Free(M2);
492 value_clear(g);
493
494 first = true;
495 value_init(offset);
496 value_init(size);
497 value_init(tmp);
498 value_set_si(size, 0);
499
500 for (int i = 0; i < P->NbBid; ++i) {
501 value_absolute(tmp, P->Ray[i][1]);
502 if (value_gt(tmp, size))
503 value_assign(size, tmp);
504 }
505 for (int i = P->NbBid; i < P->NbRays; ++i) {
506 if (value_zero_p(P->Ray[i][P->Dimension+1])) {
507 if (value_gt(P->Ray[i][1], size))
508 value_assign(size, P->Ray[i][1]);
509 continue;
510 }
511 mpz_cdiv_q(tmp, P->Ray[i][1], P->Ray[i][P->Dimension+1]);
512 if (first || value_gt(tmp, offset)) {
513 value_assign(offset, tmp);
514 first = false;
515 }
516 }
517 value_addto(offset, offset, size);
518 value_clear(size);
519 value_clear(tmp);
520
521 v = Vector_Alloc(P->Dimension+2);
522 value_set_si(v->p[0], 1);
523 value_set_si(v->p[1], -1);
524 value_assign(v->p[1+P->Dimension], offset);
525 R = AddConstraints(v->p, 1, P, options->MaxRays);
526 Polyhedron_Free(P);
527 P = R;
528
529 value_clear(offset);
530 Vector_Free(v);
531
532 value_init(c);
533 barvinok_count_with_options(P, &c, options);
534 Polyhedron_Free(P);
535 if (value_zero_p(c))
536 value_set_si(*result, 0);
537 else
538 value_set_si(*result, -1);
539 value_clear(c);
540
541 return true;
542 }
543
544 typedef Polyhedron * Polyhedron_p;
545
546 static void barvinok_count_f(Polyhedron *P, Value* result,
547 barvinok_options *options);
548
549 void barvinok_count_with_options(Polyhedron *P, Value* result,
550 struct barvinok_options *options)
551 {
552 unsigned dim;
553 int allocated = 0;
554 Polyhedron *Q;
555 bool infinite = false;
556
557 if (P->next)
558 fprintf(stderr,
559 "barvinok_count: input is a union; only first polyhedron is counted\n");
560
561 if (emptyQ2(P)) {
562 value_set_si(*result, 0);
563 return;
564 }
565 if (P->NbEq != 0) {
566 Q = NULL;
567 do {
568 P = remove_equalities(P, options->MaxRays);
569 P = DomainConstraintSimplify(P, options->MaxRays);
570 if (Q)
571 Polyhedron_Free(Q);
572 Q = P;
573 } while (!emptyQ(P) && P->NbEq != 0);
574 if (emptyQ(P)) {
575 Polyhedron_Free(P);
576 value_set_si(*result, 0);
577 return;
578 }
579 allocated = 1;
580 }
581 if (Polyhedron_is_infinite(P, result, options)) {
582 if (allocated)
583 Polyhedron_Free(P);
584 return;
585 }
586 if (P->Dimension == 0) {
587 /* Test whether the constraints are satisfied */
588 POL_ENSURE_VERTICES(P);
589 value_set_si(*result, !emptyQ(P));
590 if (allocated)
591 Polyhedron_Free(P);
592 return;
593 }
594 Q = Polyhedron_Factor(P, 0, NULL, options->MaxRays);
595 if (Q) {
596 if (allocated)
597 Polyhedron_Free(P);
598 P = Q;
599 allocated = 1;
600 }
601
602 barvinok_count_f(P, result, options);
603 if (value_neg_p(*result))
604 infinite = true;
605 if (Q && P->next && value_notzero_p(*result)) {
606 Value factor;
607 value_init(factor);
608
609 for (Q = P->next; Q; Q = Q->next) {
610 barvinok_count_f(Q, &factor, options);
611 if (value_neg_p(factor)) {
612 infinite = true;
613 continue;
614 } else if (Q->next && value_zero_p(factor)) {
615 value_set_si(*result, 0);
616 break;
617 }
618 value_multiply(*result, *result, factor);
619 }
620
621 value_clear(factor);
622 }
623
624 if (allocated)
625 Domain_Free(P);
626 if (infinite)
627 value_set_si(*result, -1);
628 }
629
630 void barvinok_count(Polyhedron *P, Value* result, unsigned NbMaxCons)
631 {
632 barvinok_options *options = barvinok_options_new_with_defaults();
633 options->MaxRays = NbMaxCons;
634 barvinok_count_with_options(P, result, options);
635 barvinok_options_free(options);
636 }
637
638 static void barvinok_count_f(Polyhedron *P, Value* result,
639 barvinok_options *options)
640 {
641 if (emptyQ2(P)) {
642 value_set_si(*result, 0);
643 return;
644 }
645
646 if (P->Dimension == 1)
647 return Line_Length(P, result);
648
649 int c = P->NbConstraints;
650 POL_ENSURE_FACETS(P);
651 if (c != P->NbConstraints || P->NbEq != 0) {
652 Polyhedron *next = P->next;
653 P->next = NULL;
654 barvinok_count_with_options(P, result, options);
655 P->next = next;
656 return;
657 }
658
659 POL_ENSURE_VERTICES(P);
660
661 if (Polyhedron_is_infinite(P, result, options))
662 return;
663
664 np_base *cnt;
665 if (options->incremental_specialization == BV_SPECIALIZATION_BF)
666 cnt = new bfcounter(P->Dimension);
667 else if (options->incremental_specialization == BV_SPECIALIZATION_DF)
668 cnt = new icounter(P->Dimension);
669 else if (options->incremental_specialization == BV_SPECIALIZATION_TODD)
670 cnt = new tcounter(P->Dimension, options->max_index);
671 else
672 cnt = new counter(P->Dimension, options->max_index);
673 cnt->start(P, options);
674
675 cnt->get_count(result);
676 delete cnt;
677 }
678
679 static void uni_polynom(int param, Vector *c, evalue *EP)
680 {
681 unsigned dim = c->Size-2;
682 value_init(EP->d);
683 value_set_si(EP->d,0);
684 EP->x.p = new_enode(polynomial, dim+1, param+1);
685 for (int j = 0; j <= dim; ++j)
686 evalue_set(&EP->x.p->arr[j], c->p[j], c->p[dim+1]);
687 }
688
689 Polyhedron *unfringe (Polyhedron *P, unsigned MaxRays)
690 {
691 int len = P->Dimension+2;
692 Polyhedron *T, *R = P;
693 Value g;
694 value_init(g);
695 Vector *row = Vector_Alloc(len);
696 value_set_si(row->p[0], 1);
697
698 R = DomainConstraintSimplify(Polyhedron_Copy(P), MaxRays);
699
700 Matrix *M = Matrix_Alloc(2, len-1);
701 value_set_si(M->p[1][len-2], 1);
702 for (int v = 0; v < P->Dimension; ++v) {
703 value_set_si(M->p[0][v], 1);
704 Polyhedron *I = Polyhedron_Image(R, M, 2+1);
705 value_set_si(M->p[0][v], 0);
706 for (int r = 0; r < I->NbConstraints; ++r) {
707 if (value_zero_p(I->Constraint[r][0]))
708 continue;
709 if (value_zero_p(I->Constraint[r][1]))
710 continue;
711 if (value_one_p(I->Constraint[r][1]))
712 continue;
713 if (value_mone_p(I->Constraint[r][1]))
714 continue;
715 value_absolute(g, I->Constraint[r][1]);
716 Vector_Set(row->p+1, 0, len-2);
717 value_division(row->p[1+v], I->Constraint[r][1], g);
718 mpz_fdiv_q(row->p[len-1], I->Constraint[r][2], g);
719 T = R;
720 R = AddConstraints(row->p, 1, R, MaxRays);
721 if (T != P)
722 Polyhedron_Free(T);
723 }
724 Polyhedron_Free(I);
725 }
726 Matrix_Free(M);
727 Vector_Free(row);
728 value_clear(g);
729 return R;
730 }
731
732 /* Check whether all rays point in the positive directions
733 * for the parameters
734 */
735 static bool Polyhedron_has_positive_rays(Polyhedron *P, unsigned nparam)
736 {
737 int r;
738 for (r = 0; r < P->NbRays; ++r)
739 if (value_zero_p(P->Ray[r][P->Dimension+1])) {
740 int i;
741 for (i = P->Dimension - nparam; i < P->Dimension; ++i)
742 if (value_neg_p(P->Ray[r][i+1]))
743 return false;
744 }
745 return true;
746 }
747
748 typedef evalue * evalue_p;
749
750 struct enumerator_base {
751 unsigned dim;
752 evalue ** vE;
753 evalue mone;
754 vertex_decomposer *vpd;
755
756 enumerator_base(unsigned dim, vertex_decomposer *vpd)
757 {
758 this->dim = dim;
759 this->vpd = vpd;
760
761 vE = new evalue_p[vpd->nbV];
762 for (int j = 0; j < vpd->nbV; ++j)
763 vE[j] = 0;
764
765 value_init(mone.d);
766 evalue_set_si(&mone, -1, 1);
767 }
768
769 void decompose_at(Param_Vertices *V, int _i, barvinok_options *options) {
770 //this->pVD = pVD;
771
772 vE[_i] = new evalue;
773 value_init(vE[_i]->d);
774 evalue_set_si(vE[_i], 0, 1);
775
776 vpd->decompose_at_vertex(V, _i, options);
777 }
778
779 virtual ~enumerator_base() {
780 for (int j = 0; j < vpd->nbV; ++j)
781 if (vE[j]) {
782 free_evalue_refs(vE[j]);
783 delete vE[j];
784 }
785 delete [] vE;
786
787 free_evalue_refs(&mone);
788 }
789
790 static enumerator_base *create(Polyhedron *P, unsigned dim, unsigned nbV,
791 barvinok_options *options);
792 };
793
794 struct enumerator : public signed_cone_consumer, public vertex_decomposer,
795 public enumerator_base {
796 vec_ZZ lambda;
797 vec_ZZ den;
798 ZZ sign;
799 term_info num;
800 Vector *c;
801 mpq_t count;
802 Value tz;
803
804 enumerator(Polyhedron *P, unsigned dim, unsigned nbV) :
805 vertex_decomposer(P, nbV, *this), enumerator_base(dim, this) {
806 this->P = P;
807 this->nbV = nbV;
808 randomvector(P, lambda, dim);
809 den.SetLength(dim);
810 c = Vector_Alloc(dim+2);
811
812 mpq_init(count);
813 value_init(tz);
814 }
815
816 ~enumerator() {
817 mpq_clear(count);
818 Vector_Free(c);
819 value_clear(tz);
820 }
821
822 virtual void handle(const signed_cone& sc, barvinok_options *options);
823 };
824
825 void enumerator::handle(const signed_cone& sc, barvinok_options *options)
826 {
827 int r = 0;
828 assert(sc.rays.NumRows() == dim);
829 for (int k = 0; k < dim; ++k) {
830 if (lambda * sc.rays[k] == 0)
831 throw Orthogonal;
832 }
833
834 sign = sc.sign;
835
836 lattice_point(V, sc.rays, lambda, &num, sc.det, sc.closed, options);
837 den = sc.rays * lambda;
838
839 if (dim % 2)
840 sign = -sign;
841
842 zz2value(den[0], tz);
843 dpoly n(dim, tz, 1);
844 for (int k = 1; k < dim; ++k) {
845 zz2value(den[k], tz);
846 dpoly fact(dim, tz, 1);
847 n *= fact;
848 }
849 if (num.E != NULL) {
850 dpoly_n d(dim);
851 d.div(n, c, sign);
852 for (unsigned long i = 0; i < sc.det; ++i) {
853 evalue *EV = evalue_polynomial(c, num.E[i]);
854 eadd(EV, vE[vert]);
855 free_evalue_refs(EV);
856 free(EV);
857 free_evalue_refs(num.E[i]);
858 delete num.E[i];
859 }
860 delete [] num.E;
861 } else {
862 mpq_set_si(count, 0, 1);
863 if (num.constant.length() == 1) {
864 zz2value(num.constant[0], tz);
865 dpoly d(dim, tz);
866 d.div(n, count, sign);
867 } else {
868 dpoly_n d(dim);
869 d.div(n, c, sign);
870 Value x, sum, acc;
871 value_init(x);
872 value_init(acc);
873 for (unsigned long i = 0; i < sc.det; ++i) {
874 value_assign(acc, c->p[dim]);
875 zz2value(num.constant[i], x);
876 for (int j = dim-1; j >= 0; --j) {
877 value_multiply(acc, acc, x);
878 value_addto(acc, acc, c->p[j]);
879 }
880 value_addto(mpq_numref(count), mpq_numref(count), acc);
881 }
882 mpz_set(mpq_denref(count), c->p[dim+1]);
883 value_clear(acc);
884 value_clear(x);
885 }
886 evalue EV;
887 value_init(EV.d);
888 evalue_set(&EV, &count[0]._mp_num, &count[0]._mp_den);
889 eadd(&EV, vE[vert]);
890 free_evalue_refs(&EV);
891 }
892 }
893
894 struct ienumerator_base : enumerator_base {
895 evalue ** E_vertex;
896
897 ienumerator_base(unsigned dim, vertex_decomposer *vpd) :
898 enumerator_base(dim,vpd) {
899 E_vertex = new evalue_p[dim];
900 }
901
902 virtual ~ienumerator_base() {
903 delete [] E_vertex;
904 }
905
906 evalue *E_num(int i, int d) {
907 return E_vertex[i + (dim-d)];
908 }
909 };
910
911 struct cumulator {
912 evalue *factor;
913 evalue *v;
914 dpoly_r *r;
915
916 cumulator(evalue *factor, evalue *v, dpoly_r *r) :
917 factor(factor), v(v), r(r) {}
918
919 void cumulate(barvinok_options *options);
920
921 virtual void add_term(const vector<int>& powers, evalue *f2) = 0;
922 virtual ~cumulator() {}
923 };
924
925 void cumulator::cumulate(barvinok_options *options)
926 {
927 evalue cum; // factor * 1 * E_num[0]/1 * (E_num[0]-1)/2 *...
928 evalue f;
929 evalue t; // E_num[0] - (m-1)
930 evalue *cst;
931 evalue mone;
932
933 if (options->lookup_table) {
934 value_init(mone.d);
935 evalue_set_si(&mone, -1, 1);
936 }
937
938 value_init(cum.d);
939 evalue_copy(&cum, factor);
940 value_init(f.d);
941 value_init(f.x.n);
942 value_set_si(f.d, 1);
943 value_set_si(f.x.n, 1);
944 value_init(t.d);
945 evalue_copy(&t, v);
946
947 if (!options->lookup_table) {
948 for (cst = &t; value_zero_p(cst->d); ) {
949 if (cst->x.p->type == fractional)
950 cst = &cst->x.p->arr[1];
951 else
952 cst = &cst->x.p->arr[0];
953 }
954 }
955
956 for (int m = 0; m < r->len; ++m) {
957 if (m > 0) {
958 if (m > 1) {
959 value_set_si(f.d, m);
960 emul(&f, &cum);
961 if (!options->lookup_table)
962 value_subtract(cst->x.n, cst->x.n, cst->d);
963 else
964 eadd(&mone, &t);
965 }
966 emul(&t, &cum);
967 }
968 dpoly_r_term_list& current = r->c[r->len-1-m];
969 dpoly_r_term_list::iterator j;
970 for (j = current.begin(); j != current.end(); ++j) {
971 if ((*j)->coeff == 0)
972 continue;
973 evalue *f2 = new evalue;
974 value_init(f2->d);
975 value_init(f2->x.n);
976 zz2value((*j)->coeff, f2->x.n);
977 zz2value(r->denom, f2->d);
978 emul(&cum, f2);
979
980 add_term((*j)->powers, f2);
981 }
982 }
983 free_evalue_refs(&f);
984 free_evalue_refs(&t);
985 free_evalue_refs(&cum);
986 if (options->lookup_table)
987 free_evalue_refs(&mone);
988 }
989
990 struct E_poly_term {
991 vector<int> powers;
992 evalue *E;
993 };
994
995 struct ie_cum : public cumulator {
996 vector<E_poly_term *> terms;
997
998 ie_cum(evalue *factor, evalue *v, dpoly_r *r) : cumulator(factor, v, r) {}
999
1000 virtual void add_term(const vector<int>& powers, evalue *f2);
1001 };
1002
1003 void ie_cum::add_term(const vector<int>& powers, evalue *f2)
1004 {
1005 int k;
1006 for (k = 0; k < terms.size(); ++k) {
1007 if (terms[k]->powers == powers) {
1008 eadd(f2, terms[k]->E);
1009 free_evalue_refs(f2);
1010 delete f2;
1011 break;
1012 }
1013 }
1014 if (k >= terms.size()) {
1015 E_poly_term *ET = new E_poly_term;
1016 ET->powers = powers;
1017 ET->E = f2;
1018 terms.push_back(ET);
1019 }
1020 }
1021
1022 struct ienumerator : public signed_cone_consumer, public vertex_decomposer,
1023 public ienumerator_base {
1024 //Polyhedron *pVD;
1025 mat_ZZ den;
1026 mat_ZZ vertex;
1027 mpq_t tcount;
1028 Value tz;
1029
1030 ienumerator(Polyhedron *P, unsigned dim, unsigned nbV) :
1031 vertex_decomposer(P, nbV, *this), ienumerator_base(dim, this) {
1032 vertex.SetDims(1, dim);
1033
1034 den.SetDims(dim, dim);
1035 mpq_init(tcount);
1036 value_init(tz);
1037 }
1038
1039 ~ienumerator() {
1040 mpq_clear(tcount);
1041 value_clear(tz);
1042 }
1043
1044 virtual void handle(const signed_cone& sc, barvinok_options *options);
1045 void reduce(evalue *factor, const mat_ZZ& num, const mat_ZZ& den_f,
1046 barvinok_options *options);
1047 };
1048
1049 void ienumerator::reduce(evalue *factor, const mat_ZZ& num, const mat_ZZ& den_f,
1050 barvinok_options *options)
1051 {
1052 unsigned len = den_f.NumRows(); // number of factors in den
1053 unsigned dim = num.NumCols();
1054 assert(num.NumRows() == 1);
1055
1056 if (dim == 0) {
1057 eadd(factor, vE[vert]);
1058 return;
1059 }
1060
1061 vec_ZZ den_s;
1062 mat_ZZ den_r;
1063 vec_ZZ num_s;
1064 mat_ZZ num_p;
1065
1066 split_one(num, num_s, num_p, den_f, den_s, den_r);
1067
1068 vec_ZZ den_p;
1069 den_p.SetLength(len);
1070
1071 ZZ one;
1072 one = 1;
1073 normalize(one, num_s, num_p, den_s, den_p, den_r);
1074 if (one != 1)
1075 emul(&mone, factor);
1076
1077 int only_param = 0;
1078 int no_param = 0;
1079 for (int k = 0; k < len; ++k) {
1080 if (den_p[k] == 0)
1081 ++no_param;
1082 else if (den_s[k] == 0)
1083 ++only_param;
1084 }
1085 if (no_param == 0) {
1086 reduce(factor, num_p, den_r, options);
1087 } else {
1088 int k, l;
1089 mat_ZZ pden;
1090 pden.SetDims(only_param, dim-1);
1091
1092 for (k = 0, l = 0; k < len; ++k)
1093 if (den_s[k] == 0)
1094 pden[l++] = den_r[k];
1095
1096 for (k = 0; k < len; ++k)
1097 if (den_p[k] == 0)
1098 break;
1099
1100 zz2value(num_s[0], tz);
1101 dpoly n(no_param, tz);
1102 zz2value(den_s[k], tz);
1103 dpoly D(no_param, tz, 1);
1104 for ( ; ++k < len; )
1105 if (den_p[k] == 0) {
1106 zz2value(den_s[k], tz);
1107 dpoly fact(no_param, tz, 1);
1108 D *= fact;
1109 }
1110
1111 dpoly_r * r = 0;
1112 // if no_param + only_param == len then all powers
1113 // below will be all zero
1114 if (no_param + only_param == len) {
1115 if (E_num(0, dim) != 0)
1116 r = new dpoly_r(n, len);
1117 else {
1118 mpq_set_si(tcount, 0, 1);
1119 one = 1;
1120 n.div(D, tcount, one);
1121
1122 if (value_notzero_p(mpq_numref(tcount))) {
1123 evalue f;
1124 value_init(f.d);
1125 value_init(f.x.n);
1126 value_assign(f.x.n, mpq_numref(tcount));
1127 value_assign(f.d, mpq_denref(tcount));
1128 emul(&f, factor);
1129 reduce(factor, num_p, pden, options);
1130 free_evalue_refs(&f);
1131 }
1132 return;
1133 }
1134 } else {
1135 for (k = 0; k < len; ++k) {
1136 if (den_s[k] == 0 || den_p[k] == 0)
1137 continue;
1138
1139 zz2value(den_s[k], tz);
1140 dpoly pd(no_param-1, tz, 1);
1141
1142 int l;
1143 for (l = 0; l < k; ++l)
1144 if (den_r[l] == den_r[k])
1145 break;
1146
1147 if (r == 0)
1148 r = new dpoly_r(n, pd, l, len);
1149 else {
1150 dpoly_r *nr = new dpoly_r(r, pd, l, len);
1151 delete r;
1152 r = nr;
1153 }
1154 }
1155 }
1156 dpoly_r *rc = r->div(D);
1157 delete r;
1158 r = rc;
1159 if (E_num(0, dim) == 0) {
1160 int common = pden.NumRows();
1161 dpoly_r_term_list& final = r->c[r->len-1];
1162 int rows;
1163 evalue t;
1164 evalue f;
1165 value_init(f.d);
1166 value_init(f.x.n);
1167 zz2value(r->denom, f.d);
1168 dpoly_r_term_list::iterator j;
1169 for (j = final.begin(); j != final.end(); ++j) {
1170 if ((*j)->coeff == 0)
1171 continue;
1172 rows = common;
1173 for (int k = 0; k < r->dim; ++k) {
1174 int n = (*j)->powers[k];
1175 if (n == 0)
1176 continue;
1177 pden.SetDims(rows+n, pden.NumCols());
1178 for (int l = 0; l < n; ++l)
1179 pden[rows+l] = den_r[k];
1180 rows += n;
1181 }
1182 value_init(t.d);
1183 evalue_copy(&t, factor);
1184 zz2value((*j)->coeff, f.x.n);
1185 emul(&f, &t);
1186 reduce(&t, num_p, pden, options);
1187 free_evalue_refs(&t);
1188 }
1189 free_evalue_refs(&f);
1190 } else {
1191 ie_cum cum(factor, E_num(0, dim), r);
1192 cum.cumulate(options);
1193
1194 int common = pden.NumRows();
1195 int rows;
1196 for (int j = 0; j < cum.terms.size(); ++j) {
1197 rows = common;
1198 pden.SetDims(rows, pden.NumCols());
1199 for (int k = 0; k < r->dim; ++k) {
1200 int n = cum.terms[j]->powers[k];
1201 if (n == 0)
1202 continue;
1203 pden.SetDims(rows+n, pden.NumCols());
1204 for (int l = 0; l < n; ++l)
1205 pden[rows+l] = den_r[k];
1206 rows += n;
1207 }
1208 reduce(cum.terms[j]->E, num_p, pden, options);
1209 free_evalue_refs(cum.terms[j]->E);
1210 delete cum.terms[j]->E;
1211 delete cum.terms[j];
1212 }
1213 }
1214 delete r;
1215 }
1216 }
1217
1218 static int type_offset(enode *p)
1219 {
1220 return p->type == fractional ? 1 :
1221 p->type == flooring ? 1 : 0;
1222 }
1223
1224 static int edegree(evalue *e)
1225 {
1226 int d = 0;
1227 enode *p;
1228
1229 if (value_notzero_p(e->d))
1230 return 0;
1231
1232 p = e->x.p;
1233 int i = type_offset(p);
1234 if (p->size-i-1 > d)
1235 d = p->size - i - 1;
1236 for (; i < p->size; i++) {
1237 int d2 = edegree(&p->arr[i]);
1238 if (d2 > d)
1239 d = d2;
1240 }
1241 return d;
1242 }
1243
1244 void ienumerator::handle(const signed_cone& sc, barvinok_options *options)
1245 {
1246 assert(sc.det == 1);
1247 assert(!sc.closed);
1248 assert(sc.rays.NumRows() == dim);
1249
1250 lattice_point(V, sc.rays, vertex[0], E_vertex, options);
1251
1252 den = sc.rays;
1253
1254 evalue one;
1255 value_init(one.d);
1256 evalue_set_si(&one, sc.sign, 1);
1257 reduce(&one, vertex, den, options);
1258 free_evalue_refs(&one);
1259
1260 for (int i = 0; i < dim; ++i)
1261 if (E_vertex[i]) {
1262 free_evalue_refs(E_vertex[i]);
1263 delete E_vertex[i];
1264 }
1265 }
1266
1267 struct bfenumerator : public vertex_decomposer, public bf_base,
1268 public ienumerator_base {
1269 evalue *factor;
1270
1271 bfenumerator(Polyhedron *P, unsigned dim, unsigned nbV) :
1272 vertex_decomposer(P, nbV, *this),
1273 bf_base(dim), ienumerator_base(dim, this) {
1274 lower = 0;
1275 factor = NULL;
1276 }
1277
1278 ~bfenumerator() {
1279 }
1280
1281 virtual void handle(const signed_cone& sc, barvinok_options *options);
1282 virtual void base(mat_ZZ& factors, bfc_vec& v);
1283
1284 bfc_term_base* new_bf_term(int len) {
1285 bfe_term* t = new bfe_term(len);
1286 return t;
1287 }
1288
1289 virtual void set_factor(bfc_term_base *t, int k, int change) {
1290 bfe_term* bfet = static_cast<bfe_term *>(t);
1291 factor = bfet->factors[k];
1292 assert(factor != NULL);
1293 bfet->factors[k] = NULL;
1294 if (change)
1295 emul(&mone, factor);
1296 }
1297
1298 virtual void set_factor(bfc_term_base *t, int k, mpq_t &q, int change) {
1299 bfe_term* bfet = static_cast<bfe_term *>(t);
1300 factor = bfet->factors[k];
1301 assert(factor != NULL);
1302 bfet->factors[k] = NULL;
1303
1304 evalue f;
1305 value_init(f.d);
1306 value_init(f.x.n);
1307 if (change)
1308 value_oppose(f.x.n, mpq_numref(q));
1309 else
1310 value_assign(f.x.n, mpq_numref(q));
1311 value_assign(f.d, mpq_denref(q));
1312 emul(&f, factor);
1313 free_evalue_refs(&f);
1314 }
1315
1316 virtual void set_factor(bfc_term_base *t, int k, const QQ& c, int change) {
1317 bfe_term* bfet = static_cast<bfe_term *>(t);
1318
1319 factor = new evalue;
1320
1321 evalue f;
1322 value_init(f.d);
1323 value_init(f.x.n);
1324 zz2value(c.n, f.x.n);
1325 if (change)
1326 value_oppose(f.x.n, f.x.n);
1327 zz2value(c.d, f.d);
1328
1329 value_init(factor->d);
1330 evalue_copy(factor, bfet->factors[k]);
1331 emul(&f, factor);
1332 free_evalue_refs(&f);
1333 }
1334
1335 void set_factor(evalue *f, int change) {
1336 if (change)
1337 emul(&mone, f);
1338 factor = f;
1339 }
1340
1341 virtual void insert_term(bfc_term_base *t, int i) {
1342 bfe_term* bfet = static_cast<bfe_term *>(t);
1343 int len = t->terms.NumRows()-1; // already increased by one
1344
1345 bfet->factors.resize(len+1);
1346 for (int j = len; j > i; --j) {
1347 bfet->factors[j] = bfet->factors[j-1];
1348 t->terms[j] = t->terms[j-1];
1349 }
1350 bfet->factors[i] = factor;
1351 factor = NULL;
1352 }
1353
1354 virtual void update_term(bfc_term_base *t, int i) {
1355 bfe_term* bfet = static_cast<bfe_term *>(t);
1356
1357 eadd(factor, bfet->factors[i]);
1358 free_evalue_refs(factor);
1359 delete factor;
1360 }
1361
1362 virtual bool constant_vertex(int dim) { return E_num(0, dim) == 0; }
1363
1364 virtual void cum(bf_reducer *bfr, bfc_term_base *t, int k, dpoly_r *r,
1365 barvinok_options *options);
1366 };
1367
1368 enumerator_base *enumerator_base::create(Polyhedron *P, unsigned dim, unsigned nbV,
1369 barvinok_options *options)
1370 {
1371 enumerator_base *eb;
1372
1373 if (options->incremental_specialization == BV_SPECIALIZATION_BF)
1374 eb = new bfenumerator(P, dim, nbV);
1375 else if (options->incremental_specialization == BV_SPECIALIZATION_DF)
1376 eb = new ienumerator(P, dim, nbV);
1377 else
1378 eb = new enumerator(P, dim, nbV);
1379
1380 return eb;
1381 }
1382
1383 struct bfe_cum : public cumulator {
1384 bfenumerator *bfe;
1385 bfc_term_base *told;
1386 int k;
1387 bf_reducer *bfr;
1388
1389 bfe_cum(evalue *factor, evalue *v, dpoly_r *r, bf_reducer *bfr,
1390 bfc_term_base *t, int k, bfenumerator *e) :
1391 cumulator(factor, v, r), told(t), k(k),
1392 bfr(bfr), bfe(e) {
1393 }
1394
1395 virtual void add_term(const vector<int>& powers, evalue *f2);
1396 };
1397
1398 void bfe_cum::add_term(const vector<int>& powers, evalue *f2)
1399 {
1400 bfr->update_powers(powers);
1401
1402 bfc_term_base * t = bfe->find_bfc_term(bfr->vn, bfr->npowers, bfr->nnf);
1403 bfe->set_factor(f2, bfr->l_changes % 2);
1404 bfe->add_term(t, told->terms[k], bfr->l_extra_num);
1405 }
1406
1407 void bfenumerator::cum(bf_reducer *bfr, bfc_term_base *t, int k,
1408 dpoly_r *r, barvinok_options *options)
1409 {
1410 bfe_term* bfet = static_cast<bfe_term *>(t);
1411 bfe_cum cum(bfet->factors[k], E_num(0, bfr->d), r, bfr, t, k, this);
1412 cum.cumulate(options);
1413 }
1414
1415 void bfenumerator::base(mat_ZZ& factors, bfc_vec& v)
1416 {
1417 for (int i = 0; i < v.size(); ++i) {
1418 assert(v[i]->terms.NumRows() == 1);
1419 evalue *factor = static_cast<bfe_term *>(v[i])->factors[0];
1420 eadd(factor, vE[vert]);
1421 delete v[i];
1422 }
1423 }
1424
1425 void bfenumerator::handle(const signed_cone& sc, barvinok_options *options)
1426 {
1427 assert(sc.det == 1);
1428 assert(!sc.closed);
1429 assert(sc.rays.NumRows() == enumerator_base::dim);
1430
1431 bfe_term* t = new bfe_term(enumerator_base::dim);
1432 vector< bfc_term_base * > v;
1433 v.push_back(t);
1434
1435 t->factors.resize(1);
1436
1437 t->terms.SetDims(1, enumerator_base::dim);
1438 lattice_point(V, sc.rays, t->terms[0], E_vertex, options);
1439
1440 // the elements of factors are always lexpositive
1441 mat_ZZ factors;
1442 int s = setup_factors(sc.rays, factors, t, sc.sign);
1443
1444 t->factors[0] = new evalue;
1445 value_init(t->factors[0]->d);
1446 evalue_set_si(t->factors[0], s, 1);
1447 reduce(factors, v, options);
1448
1449 for (int i = 0; i < enumerator_base::dim; ++i)
1450 if (E_vertex[i]) {
1451 free_evalue_refs(E_vertex[i]);
1452 delete E_vertex[i];
1453 }
1454 }
1455
1456 static inline Param_Polyhedron *Polyhedron2Param_MR(Polyhedron *Din,
1457 Polyhedron *Cin, int WS)
1458 {
1459 if (WS & POL_NO_DUAL)
1460 WS = 0;
1461 return Polyhedron2Param_Domain(Din, Cin, WS);
1462 }
1463
1464 static evalue* barvinok_enumerate_ev_f(Polyhedron *P, Polyhedron* C,
1465 barvinok_options *options);
1466
1467 /* Destroys C */
1468 static evalue* barvinok_enumerate_cst(Polyhedron *P, Polyhedron* C,
1469 struct barvinok_options *options)
1470 {
1471 evalue *eres;
1472
1473 ALLOC(evalue, eres);
1474 value_init(eres->d);
1475 value_set_si(eres->d, 0);
1476 eres->x.p = new_enode(partition, 2, C->Dimension);
1477 EVALUE_SET_DOMAIN(eres->x.p->arr[0],
1478 DomainConstraintSimplify(C, options->MaxRays));
1479 value_set_si(eres->x.p->arr[1].d, 1);
1480 value_init(eres->x.p->arr[1].x.n);
1481 if (emptyQ2(P))
1482 value_set_si(eres->x.p->arr[1].x.n, 0);
1483 else
1484 barvinok_count_with_options(P, &eres->x.p->arr[1].x.n, options);
1485
1486 return eres;
1487 }
1488
1489 /* frees P */
1490 static evalue* enumerate(Polyhedron *P, Polyhedron* C,
1491 struct barvinok_options *options)
1492 {
1493 //P = unfringe(P, MaxRays);
1494 Polyhedron *next;
1495 Polyhedron *Corig = C;
1496 Polyhedron *CEq = NULL, *rVD;
1497 int r = 0;
1498 unsigned nparam = C->Dimension;
1499 evalue *eres;
1500 Matrix *CP = NULL;
1501
1502 evalue factor;
1503 value_init(factor.d);
1504 evalue_set_si(&factor, 1, 1);
1505
1506 /* for now */
1507 POL_ENSURE_FACETS(P);
1508 POL_ENSURE_VERTICES(P);
1509 POL_ENSURE_FACETS(C);
1510 POL_ENSURE_VERTICES(C);
1511
1512 if (C->Dimension == 0 || emptyQ(P)) {
1513 constant:
1514 eres = barvinok_enumerate_cst(P, CEq ? CEq : Polyhedron_Copy(C), options);
1515 out:
1516 if (CP) {
1517 evalue_backsubstitute(eres, CP, options->MaxRays);
1518 Matrix_Free(CP);
1519 }
1520
1521 emul(&factor, eres);
1522 if (options->approximation_method == BV_APPROX_DROP) {
1523 if (options->polynomial_approximation == BV_APPROX_SIGN_UPPER)
1524 evalue_frac2polynomial(eres, 1, options->MaxRays);
1525 if (options->polynomial_approximation == BV_APPROX_SIGN_LOWER)
1526 evalue_frac2polynomial(eres, -1, options->MaxRays);
1527 if (options->polynomial_approximation == BV_APPROX_SIGN_APPROX)
1528 evalue_frac2polynomial(eres, 0, options->MaxRays);
1529 }
1530 reduce_evalue(eres);
1531 free_evalue_refs(&factor);
1532 Domain_Free(P);
1533 if (C != Corig)
1534 Polyhedron_Free(C);
1535
1536 return eres;
1537 }
1538 if (Polyhedron_is_unbounded(P, nparam, options->MaxRays))
1539 goto constant;
1540
1541 if (P->NbEq != 0) {
1542 Matrix *f;
1543 P = remove_equalities_p(P, P->Dimension-nparam, &f, options->MaxRays);
1544 mask(f, &factor, options);
1545 Matrix_Free(f);
1546 }
1547 if (P->Dimension == nparam) {
1548 CEq = P;
1549 P = Universe_Polyhedron(0);
1550 goto constant;
1551 }
1552 if (P->NbEq != 0) {
1553 Polyhedron *Q = P;
1554 Polyhedron *D = C;
1555 remove_all_equalities(&Q, &C, &CP, NULL, nparam, options->MaxRays);
1556 if (C != D && D != Corig)
1557 Polyhedron_Free(D);
1558 eres = enumerate(Q, C, options);
1559 goto out;
1560 }
1561
1562 Polyhedron *T = Polyhedron_Factor(P, nparam, NULL, options->MaxRays);
1563 if (T || (P->Dimension == nparam+1)) {
1564 Polyhedron *Q;
1565 Polyhedron *C2;
1566 for (Q = T ? T : P; Q; Q = Q->next) {
1567 Polyhedron *next = Q->next;
1568 Q->next = NULL;
1569
1570 Polyhedron *QC = Q;
1571 if (Q->Dimension != C->Dimension)
1572 QC = Polyhedron_Project(Q, nparam);
1573
1574 C2 = C;
1575 C = DomainIntersection(C, QC, options->MaxRays);
1576 if (C2 != Corig)
1577 Polyhedron_Free(C2);
1578 if (QC != Q)
1579 Polyhedron_Free(QC);
1580
1581 Q->next = next;
1582 }
1583 }
1584 if (T) {
1585 Polyhedron_Free(P);
1586 P = T;
1587 if (T->Dimension == C->Dimension) {
1588 P = T->next;
1589 T->next = NULL;
1590 Polyhedron_Free(T);
1591 }
1592 }
1593
1594 next = P->next;
1595 P->next = NULL;
1596 eres = barvinok_enumerate_ev_f(P, C, options);
1597 P->next = next;
1598
1599 if (P->next) {
1600 Polyhedron *Q;
1601 evalue *f;
1602
1603 for (Q = P->next; Q; Q = Q->next) {
1604 Polyhedron *next = Q->next;
1605 Q->next = NULL;
1606
1607 f = barvinok_enumerate_ev_f(Q, C, options);
1608 emul(f, eres);
1609 free_evalue_refs(f);
1610 free(f);
1611
1612 Q->next = next;
1613 }
1614 }
1615
1616 goto out;
1617 }
1618
1619 evalue* barvinok_enumerate_with_options(Polyhedron *P, Polyhedron* C,
1620 struct barvinok_options *options)
1621 {
1622 Polyhedron *next, *Cnext, *CA;
1623 Polyhedron *Porig = P;
1624 evalue *eres;
1625
1626 if (P->next)
1627 fprintf(stderr,
1628 "barvinok_enumerate: input is a union; only first polyhedron is enumerated\n");
1629
1630 if (C->next)
1631 fprintf(stderr,
1632 "barvinok_enumerate: context is a union; only first polyhedron is considered\n");
1633
1634 Cnext = C->next;
1635 C->next = NULL;
1636 CA = align_context(C, P->Dimension, options->MaxRays);
1637 next = P->next;
1638 P->next = NULL;
1639 P = DomainIntersection(P, CA, options->MaxRays);
1640 Porig->next = next;
1641 Polyhedron_Free(CA);
1642
1643 if (options->approximation_method == BV_APPROX_BERNOULLI) {
1644 eres = Bernoulli_sum(P, C, options);
1645 Domain_Free(P);
1646 } else
1647 eres = enumerate(P, C, options);
1648
1649 C->next = Cnext;
1650
1651 return eres;
1652 }
1653
1654 evalue* barvinok_enumerate_ev(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
1655 {
1656 evalue *E;
1657 barvinok_options *options = barvinok_options_new_with_defaults();
1658 options->MaxRays = MaxRays;
1659 E = barvinok_enumerate_with_options(P, C, options);
1660 barvinok_options_free(options);
1661 return E;
1662 }
1663
1664 evalue *Param_Polyhedron_Enumerate(Param_Polyhedron *PP, Polyhedron *P,
1665 Polyhedron *C,
1666 struct barvinok_options *options)
1667 {
1668 evalue *eres;
1669 Param_Domain *D;
1670 unsigned nparam = C->Dimension;
1671 unsigned dim = P->Dimension - nparam;
1672
1673 ALLOC(evalue, eres);
1674 value_init(eres->d);
1675 value_set_si(eres->d, 0);
1676
1677 int nd;
1678 for (nd = 0, D=PP->D; D; ++nd, D=D->next);
1679 struct section { Polyhedron *D; evalue E; };
1680 section *s = new section[nd];
1681
1682 enumerator_base *et = NULL;
1683 try_again:
1684 if (et)
1685 delete et;
1686
1687 et = enumerator_base::create(P, dim, PP->nbV, options);
1688
1689 Polyhedron *TC = true_context(P, C, options->MaxRays);
1690 FORALL_REDUCED_DOMAIN(PP, TC, nd, options, i, D, rVD)
1691 Param_Vertices *V;
1692
1693 value_init(s[i].E.d);
1694 evalue_set_si(&s[i].E, 0, 1);
1695 s[i].D = rVD;
1696
1697 FORALL_PVertex_in_ParamPolyhedron(V,D,PP) // _i is internal counter
1698 if (!et->vE[_i])
1699 try {
1700 et->decompose_at(V, _i, options);
1701 } catch (OrthogonalException &e) {
1702 FORALL_REDUCED_DOMAIN_RESET;
1703 for (; i >= 0; --i) {
1704 free_evalue_refs(&s[i].E);
1705 Domain_Free(s[i].D);
1706 }
1707 goto try_again;
1708 }
1709 eadd(et->vE[_i] , &s[i].E);
1710 END_FORALL_PVertex_in_ParamPolyhedron;
1711 evalue_range_reduction_in_domain(&s[i].E, rVD);
1712 END_FORALL_REDUCED_DOMAIN
1713 Polyhedron_Free(TC);
1714
1715 delete et;
1716 if (nd == 0)
1717 evalue_set_si(eres, 0, 1);
1718 else {
1719 eres->x.p = new_enode(partition, 2*nd, C->Dimension);
1720 for (int j = 0; j < nd; ++j) {
1721 EVALUE_SET_DOMAIN(eres->x.p->arr[2*j], s[j].D);
1722 value_clear(eres->x.p->arr[2*j+1].d);
1723 eres->x.p->arr[2*j+1] = s[j].E;
1724 }
1725 }
1726 delete [] s;
1727
1728 return eres;
1729 }
1730
1731 static evalue* barvinok_enumerate_ev_f(Polyhedron *P, Polyhedron* C,
1732 barvinok_options *options)
1733 {
1734 unsigned nparam = C->Dimension;
1735 bool do_scale = options->approximation_method == BV_APPROX_SCALE;
1736
1737 if (options->approximation_method == BV_APPROX_VOLUME)
1738 return Param_Polyhedron_Volume(P, C, options);
1739
1740 if (P->Dimension - nparam == 1 && !do_scale)
1741 return ParamLine_Length(P, C, options);
1742
1743 Param_Polyhedron *PP = NULL;
1744 evalue *eres;
1745
1746 if (do_scale) {
1747 eres = scale_bound(P, C, options);
1748 if (eres)
1749 return eres;
1750 }
1751
1752 PP = Polyhedron2Param_MR(P, C, options->MaxRays);
1753
1754 if (do_scale)
1755 eres = scale(PP, P, C, options);
1756 else
1757 eres = Param_Polyhedron_Enumerate(PP, P, C, options);
1758
1759 if (PP)
1760 Param_Polyhedron_Free(PP);
1761
1762 return eres;
1763 }
1764
1765 Enumeration* barvinok_enumerate(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
1766 {
1767 evalue *EP = barvinok_enumerate_ev(P, C, MaxRays);
1768
1769 return partition2enumeration(EP);
1770 }
1771
1772 static void SwapColumns(Value **V, int n, int i, int j)
1773 {
1774 for (int r = 0; r < n; ++r)
1775 value_swap(V[r][i], V[r][j]);
1776 }
1777
1778 static void SwapColumns(Polyhedron *P, int i, int j)
1779 {
1780 SwapColumns(P->Constraint, P->NbConstraints, i, j);
1781 SwapColumns(P->Ray, P->NbRays, i, j);
1782 }
1783
1784 /* Construct a constraint c from constraints l and u such that if
1785 * if constraint c holds then for each value of the other variables
1786 * there is at most one value of variable pos (position pos+1 in the constraints).
1787 *
1788 * Given a lower and an upper bound
1789 * n_l v_i + <c_l,x> + c_l >= 0
1790 * -n_u v_i + <c_u,x> + c_u >= 0
1791 * the constructed constraint is
1792 *
1793 * -(n_l<c_u,x> + n_u<c_l,x>) + (-n_l c_u - n_u c_l + n_l n_u - 1)
1794 *
1795 * which is then simplified to remove the content of the non-constant coefficients
1796 *
1797 * len is the total length of the constraints.
1798 * v is a temporary variable that can be used by this procedure
1799 */
1800 static void negative_test_constraint(Value *l, Value *u, Value *c, int pos,
1801 int len, Value *v)
1802 {
1803 value_oppose(*v, u[pos+1]);
1804 Vector_Combine(l+1, u+1, c+1, *v, l[pos+1], len-1);
1805 value_multiply(*v, *v, l[pos+1]);
1806 value_subtract(c[len-1], c[len-1], *v);
1807 value_set_si(*v, -1);
1808 Vector_Scale(c+1, c+1, *v, len-1);
1809 value_decrement(c[len-1], c[len-1]);
1810 ConstraintSimplify(c, c, len, v);
1811 }
1812
1813 static bool parallel_constraints(Value *l, Value *u, Value *c, int pos,
1814 int len)
1815 {
1816 bool parallel;
1817 Value g1;
1818 Value g2;
1819 value_init(g1);
1820 value_init(g2);
1821
1822 Vector_Gcd(&l[1+pos], len, &g1);
1823 Vector_Gcd(&u[1+pos], len, &g2);
1824 Vector_Combine(l+1+pos, u+1+pos, c+1, g2, g1, len);
1825 parallel = First_Non_Zero(c+1, len) == -1;
1826
1827 value_clear(g1);
1828 value_clear(g2);
1829
1830 return parallel;
1831 }
1832
1833 static void negative_test_constraint7(Value *l, Value *u, Value *c, int pos,
1834 int exist, int len, Value *v)
1835 {
1836 Value g;
1837 value_init(g);
1838
1839 Vector_Gcd(&u[1+pos], exist, v);
1840 Vector_Gcd(&l[1+pos], exist, &g);
1841 Vector_Combine(l+1, u+1, c+1, *v, g, len-1);
1842 value_multiply(*v, *v, g);
1843 value_subtract(c[len-1], c[len-1], *v);
1844 value_set_si(*v, -1);
1845 Vector_Scale(c+1, c+1, *v, len-1);
1846 value_decrement(c[len-1], c[len-1]);
1847 ConstraintSimplify(c, c, len, v);
1848
1849 value_clear(g);
1850 }
1851
1852 /* Turns a x + b >= 0 into a x + b <= -1
1853 *
1854 * len is the total length of the constraint.
1855 * v is a temporary variable that can be used by this procedure
1856 */
1857 static void oppose_constraint(Value *c, int len, Value *v)
1858 {
1859 value_set_si(*v, -1);
1860 Vector_Scale(c+1, c+1, *v, len-1);
1861 value_decrement(c[len-1], c[len-1]);
1862 }
1863
1864 /* Split polyhedron P into two polyhedra *pos and *neg, where
1865 * existential variable i has at most one solution for each
1866 * value of the other variables in *neg.
1867 *
1868 * The splitting is performed using constraints l and u.
1869 *
1870 * nvar: number of set variables
1871 * row: temporary vector that can be used by this procedure
1872 * f: temporary value that can be used by this procedure
1873 */
1874 static bool SplitOnConstraint(Polyhedron *P, int i, int l, int u,
1875 int nvar, int MaxRays, Vector *row, Value& f,
1876 Polyhedron **pos, Polyhedron **neg)
1877 {
1878 negative_test_constraint(P->Constraint[l], P->Constraint[u],
1879 row->p, nvar+i, P->Dimension+2, &f);
1880 *neg = AddConstraints(row->p, 1, P, MaxRays);
1881 POL_ENSURE_VERTICES(*neg);
1882
1883 /* We found an independent, but useless constraint
1884 * Maybe we should detect this earlier and not
1885 * mark the variable as INDEPENDENT
1886 */
1887 if (emptyQ((*neg))) {
1888 Polyhedron_Free(*neg);
1889 return false;
1890 }
1891
1892 oppose_constraint(row->p, P->Dimension+2, &f);
1893 *pos = AddConstraints(row->p, 1, P, MaxRays);
1894 POL_ENSURE_VERTICES(*pos);
1895
1896 if (emptyQ((*pos))) {
1897 Polyhedron_Free(*neg);
1898 Polyhedron_Free(*pos);
1899 return false;
1900 }
1901
1902 return true;
1903 }
1904
1905 /*
1906 * unimodularly transform P such that constraint r is transformed
1907 * into a constraint that involves only a single (the first)
1908 * existential variable
1909 *
1910 */
1911 static Polyhedron *rotate_along(Polyhedron *P, int r, int nvar, int exist,
1912 unsigned MaxRays)
1913 {
1914 Value g;
1915 value_init(g);
1916
1917 Matrix *M = Matrix_Alloc(exist, exist);
1918 Vector_Copy(P->Constraint[r]+1+nvar, M->p[0], exist);
1919 Vector_Gcd(M->p[0], exist, &g);
1920 if (value_notone_p(g))
1921 Vector_AntiScale(M->p[0], M->p[0], g, exist);
1922 value_clear(g);
1923
1924 int ok = unimodular_complete(M, 1);
1925 assert(ok);
1926 Matrix *M2 = Matrix_Alloc(P->Dimension+1, P->Dimension+1);
1927 for (r = 0; r < nvar; ++r)
1928 value_set_si(M2->p[r][r], 1);
1929 for ( ; r < nvar+exist; ++r)
1930 Vector_Copy(M->p[r-nvar], M2->p[r]+nvar, exist);
1931 for ( ; r < P->Dimension+1; ++r)
1932 value_set_si(M2->p[r][r], 1);
1933 Polyhedron *T = Polyhedron_Image(P, M2, MaxRays);
1934
1935 Matrix_Free(M2);
1936 Matrix_Free(M);
1937
1938 return T;
1939 }
1940
1941 /* Split polyhedron P into two polyhedra *pos and *neg, where
1942 * existential variable i has at most one solution for each
1943 * value of the other variables in *neg.
1944 *
1945 * If independent is set, then the two constraints on which the
1946 * split will be performed need to be independent of the other
1947 * existential variables.
1948 *
1949 * Return true if an appropriate split could be performed.
1950 *
1951 * nvar: number of set variables
1952 * exist: number of existential variables
1953 * row: temporary vector that can be used by this procedure
1954 * f: temporary value that can be used by this procedure
1955 */
1956 static bool SplitOnVar(Polyhedron *P, int i,
1957 int nvar, int exist, int MaxRays,
1958 Vector *row, Value& f, bool independent,
1959 Polyhedron **pos, Polyhedron **neg)
1960 {
1961 int j;
1962
1963 for (int l = P->NbEq; l < P->NbConstraints; ++l) {
1964 if (value_negz_p(P->Constraint[l][nvar+i+1]))
1965 continue;
1966
1967 if (independent) {
1968 for (j = 0; j < exist; ++j)
1969 if (j != i && value_notzero_p(P->Constraint[l][nvar+j+1]))
1970 break;
1971 if (j < exist)
1972 continue;
1973 }
1974
1975 for (int u = P->NbEq; u < P->NbConstraints; ++u) {
1976 if (value_posz_p(P->Constraint[u][nvar+i+1]))
1977 continue;
1978
1979 if (independent) {
1980 for (j = 0; j < exist; ++j)
1981 if (j != i && value_notzero_p(P->Constraint[u][nvar+j+1]))
1982 break;
1983 if (j < exist)
1984 continue;
1985 }
1986
1987 if (SplitOnConstraint(P, i, l, u, nvar, MaxRays, row, f, pos, neg)) {
1988 if (independent) {
1989 if (i != 0)
1990 SwapColumns(*neg, nvar+1, nvar+1+i);
1991 }
1992 return true;
1993 }
1994 }
1995 }
1996
1997 return false;
1998 }
1999
2000 static bool double_bound_pair(Polyhedron *P, int nvar, int exist,
2001 int i, int l1, int l2,
2002 Polyhedron **pos, Polyhedron **neg)
2003 {
2004 Value f;
2005 value_init(f);
2006 Vector *row = Vector_Alloc(P->Dimension+2);
2007 value_set_si(row->p[0], 1);
2008 value_oppose(f, P->Constraint[l1][nvar+i+1]);
2009 Vector_Combine(P->Constraint[l1]+1, P->Constraint[l2]+1,
2010 row->p+1,
2011 P->Constraint[l2][nvar+i+1], f,
2012 P->Dimension+1);
2013 ConstraintSimplify(row->p, row->p, P->Dimension+2, &f);
2014 *pos = AddConstraints(row->p, 1, P, 0);
2015 POL_ENSURE_VERTICES(*pos);
2016 value_set_si(f, -1);
2017 Vector_Scale(row->p+1, row->p+1, f, P->Dimension+1);
2018 value_decrement(row->p[P->Dimension+1], row->p[P->Dimension+1]);
2019 *neg = AddConstraints(row->p, 1, P, 0);
2020 POL_ENSURE_VERTICES(*neg);
2021 Vector_Free(row);
2022 value_clear(f);
2023
2024 return !emptyQ((*pos)) && !emptyQ((*neg));
2025 }
2026
2027 static bool double_bound(Polyhedron *P, int nvar, int exist,
2028 Polyhedron **pos, Polyhedron **neg)
2029 {
2030 for (int i = 0; i < exist; ++i) {
2031 int l1, l2;
2032 for (l1 = P->NbEq; l1 < P->NbConstraints; ++l1) {
2033 if (value_negz_p(P->Constraint[l1][nvar+i+1]))
2034 continue;
2035 for (l2 = l1 + 1; l2 < P->NbConstraints; ++l2) {
2036 if (value_negz_p(P->Constraint[l2][nvar+i+1]))
2037 continue;
2038 if (double_bound_pair(P, nvar, exist, i, l1, l2, pos, neg))
2039 return true;
2040 }
2041 }
2042 for (l1 = P->NbEq; l1 < P->NbConstraints; ++l1) {
2043 if (value_posz_p(P->Constraint[l1][nvar+i+1]))
2044 continue;
2045 if (l1 < P->NbConstraints)
2046 for (l2 = l1 + 1; l2 < P->NbConstraints; ++l2) {
2047 if (value_posz_p(P->Constraint[l2][nvar+i+1]))
2048 continue;
2049 if (double_bound_pair(P, nvar, exist, i, l1, l2, pos, neg))
2050 return true;
2051 }
2052 }
2053 return false;
2054 }
2055 return false;
2056 }
2057
2058 enum constraint {
2059 ALL_POS = 1 << 0,
2060 ONE_NEG = 1 << 1,
2061 INDEPENDENT = 1 << 2,
2062 ROT_NEG = 1 << 3
2063 };
2064
2065 static evalue* enumerate_or(Polyhedron *D,
2066 unsigned exist, unsigned nparam, barvinok_options *options)
2067 {
2068 #ifdef DEBUG_ER
2069 fprintf(stderr, "\nER: Or\n");
2070 #endif /* DEBUG_ER */
2071
2072 Polyhedron *N = D->next;
2073 D->next = 0;
2074 evalue *EP =
2075 barvinok_enumerate_e_with_options(D, exist, nparam, options);
2076 Polyhedron_Free(D);
2077
2078 for (D = N; D; D = N) {
2079 N = D->next;
2080 D->next = 0;
2081
2082 evalue *EN =
2083 barvinok_enumerate_e_with_options(D, exist, nparam, options);
2084
2085 eor(EN, EP);
2086 free_evalue_refs(EN);
2087 free(EN);
2088 Polyhedron_Free(D);
2089 }
2090
2091 reduce_evalue(EP);
2092
2093 return EP;
2094 }
2095
2096 static evalue* enumerate_sum(Polyhedron *P,
2097 unsigned exist, unsigned nparam, barvinok_options *options)
2098 {
2099 int nvar = P->Dimension - exist - nparam;
2100 int toswap = nvar < exist ? nvar : exist;
2101 for (int i = 0; i < toswap; ++i)
2102 SwapColumns(P, 1 + i, nvar+exist - i);
2103 nparam += nvar;
2104
2105 #ifdef DEBUG_ER
2106 fprintf(stderr, "\nER: Sum\n");
2107 #endif /* DEBUG_ER */
2108
2109 evalue *EP = barvinok_enumerate_e_with_options(P, exist, nparam, options);
2110
2111 evalue_split_domains_into_orthants(EP, options->MaxRays);
2112 reduce_evalue(EP);
2113 evalue_range_reduction(EP);
2114
2115 evalue_frac2floor(EP);
2116
2117 evalue *sum = evalue_sum(EP, nvar, options->MaxRays);
2118
2119 free_evalue_refs(EP);
2120 free(EP);
2121 EP = sum;
2122
2123 evalue_range_reduction(EP);
2124
2125 return EP;
2126 }
2127
2128 static evalue* split_sure(Polyhedron *P, Polyhedron *S,
2129 unsigned exist, unsigned nparam, barvinok_options *options)
2130 {
2131 int nvar = P->Dimension - exist - nparam;
2132
2133 Matrix *M = Matrix_Alloc(exist, S->Dimension+2);
2134 for (int i = 0; i < exist; ++i)
2135 value_set_si(M->p[i][nvar+i+1], 1);
2136 Polyhedron *O = S;
2137 S = DomainAddRays(S, M, options->MaxRays);
2138 Polyhedron_Free(O);
2139 Polyhedron *F = DomainAddRays(P, M, options->MaxRays);
2140 Polyhedron *D = DomainDifference(F, S, options->MaxRays);
2141 O = D;
2142 D = Disjoint_Domain(D, 0, options->MaxRays);
2143 Polyhedron_Free(F);
2144 Domain_Free(O);
2145 Matrix_Free(M);
2146
2147 M = Matrix_Alloc(P->Dimension+1-exist, P->Dimension+1);
2148 for (int j = 0; j < nvar; ++j)
2149 value_set_si(M->p[j][j], 1);
2150 for (int j = 0; j < nparam+1; ++j)
2151 value_set_si(M->p[nvar+j][nvar+exist+j], 1);
2152 Polyhedron *T = Polyhedron_Image(S, M, options->MaxRays);
2153 evalue *EP = barvinok_enumerate_e_with_options(T, 0, nparam, options);
2154 Polyhedron_Free(S);
2155 Polyhedron_Free(T);
2156 Matrix_Free(M);
2157
2158 for (Polyhedron *Q = D; Q; Q = Q->next) {
2159 Polyhedron *N = Q->next;
2160 Q->next = 0;
2161 T = DomainIntersection(P, Q, options->MaxRays);
2162 evalue *E = barvinok_enumerate_e_with_options(T, exist, nparam, options);
2163 eadd(E, EP);
2164 free_evalue_refs(E);
2165 free(E);
2166 Polyhedron_Free(T);
2167 Q->next = N;
2168 }
2169 Domain_Free(D);
2170 return EP;
2171 }
2172
2173 static evalue* enumerate_sure(Polyhedron *P,
2174 unsigned exist, unsigned nparam, barvinok_options *options)
2175 {
2176 int i;
2177 Polyhedron *S = P;
2178 int nvar = P->Dimension - exist - nparam;
2179 Value lcm;
2180 Value f;
2181 value_init(lcm);
2182 value_init(f);
2183
2184 for (i = 0; i < exist; ++i) {
2185 Matrix *M = Matrix_Alloc(S->NbConstraints, S->Dimension+2);
2186 int c = 0;
2187 value_set_si(lcm, 1);
2188 for (int j = 0; j < S->NbConstraints; ++j) {
2189 if (value_negz_p(S->Constraint[j][1+nvar+i]))
2190 continue;
2191 if (value_one_p(S->Constraint[j][1+nvar+i]))
2192 continue;
2193 value_lcm(lcm, S->Constraint[j][1+nvar+i], &lcm);
2194 }
2195
2196 for (int j = 0; j < S->NbConstraints; ++j) {
2197 if (value_negz_p(S->Constraint[j][1+nvar+i]))
2198 continue;
2199 if (value_one_p(S->Constraint[j][1+nvar+i]))
2200 continue;
2201 value_division(f, lcm, S->Constraint[j][1+nvar+i]);
2202 Vector_Scale(S->Constraint[j], M->p[c], f, S->Dimension+2);
2203 value_subtract(M->p[c][S->Dimension+1],
2204 M->p[c][S->Dimension+1],
2205 lcm);
2206 value_increment(M->p[c][S->Dimension+1],
2207 M->p[c][S->Dimension+1]);
2208 ++c;
2209 }
2210 Polyhedron *O = S;
2211 S = AddConstraints(M->p[0], c, S, options->MaxRays);
2212 if (O != P)
2213 Polyhedron_Free(O);
2214 Matrix_Free(M);
2215 if (emptyQ(S)) {
2216 Polyhedron_Free(S);
2217 value_clear(lcm);
2218 value_clear(f);
2219 return 0;
2220 }
2221 }
2222 value_clear(lcm);
2223 value_clear(f);
2224
2225 #ifdef DEBUG_ER
2226 fprintf(stderr, "\nER: Sure\n");
2227 #endif /* DEBUG_ER */
2228
2229 return split_sure(P, S, exist, nparam, options);
2230 }
2231
2232 static evalue* enumerate_sure2(Polyhedron *P,
2233 unsigned exist, unsigned nparam, barvinok_options *options)
2234 {
2235 int nvar = P->Dimension - exist - nparam;
2236 int r;
2237 for (r = 0; r < P->NbRays; ++r)
2238 if (value_one_p(P->Ray[r][0]) &&
2239 value_one_p(P->Ray[r][P->Dimension+1]))
2240 break;
2241
2242 if (r >= P->NbRays)
2243 return 0;
2244
2245 Matrix *M = Matrix_Alloc(nvar + 1 + nparam, P->Dimension+2);
2246 for (int i = 0; i < nvar; ++i)
2247 value_set_si(M->p[i][1+i], 1);
2248 for (int i = 0; i < nparam; ++i)
2249 value_set_si(M->p[i+nvar][1+nvar+exist+i], 1);
2250 Vector_Copy(P->Ray[r]+1+nvar, M->p[nvar+nparam]+1+nvar, exist);
2251 value_set_si(M->p[nvar+nparam][0], 1);
2252 value_set_si(M->p[nvar+nparam][P->Dimension+1], 1);
2253 Polyhedron * F = Rays2Polyhedron(M, options->MaxRays);
2254 Matrix_Free(M);
2255
2256 Polyhedron *I = DomainIntersection(F, P, options->MaxRays);
2257 Polyhedron_Free(F);
2258
2259 #ifdef DEBUG_ER
2260 fprintf(stderr, "\nER: Sure2\n");
2261 #endif /* DEBUG_ER */
2262
2263 return split_sure(P, I, exist, nparam, options);
2264 }
2265
2266 static evalue* enumerate_cyclic(Polyhedron *P,
2267 unsigned exist, unsigned nparam,
2268 evalue * EP, int r, int p, unsigned MaxRays)
2269 {
2270 int nvar = P->Dimension - exist - nparam;
2271
2272 /* If EP in its fractional maps only contains references
2273 * to the remainder parameter with appropriate coefficients
2274 * then we could in principle avoid adding existentially
2275 * quantified variables to the validity domains.
2276 * We'd have to replace the remainder by m { p/m }
2277 * and multiply with an appropriate factor that is one
2278 * only in the appropriate range.
2279 * This last multiplication can be avoided if EP
2280 * has a single validity domain with no (further)
2281 * constraints on the remainder parameter
2282 */
2283
2284 Matrix *CT = Matrix_Alloc(nparam+1, nparam+3);
2285 Matrix *M = Matrix_Alloc(1, 1+nparam+3);
2286 for (int j = 0; j < nparam; ++j)
2287 if (j != p)
2288 value_set_si(CT->p[j][j], 1);
2289 value_set_si(CT->p[p][nparam+1], 1);
2290 value_set_si(CT->p[nparam][nparam+2], 1);
2291 value_set_si(M->p[0][1+p], -1);
2292 value_absolute(M->p[0][1+nparam], P->Ray[0][1+nvar+exist+p]);
2293 value_set_si(M->p[0][1+nparam+1], 1);
2294 Polyhedron *CEq = Constraints2Polyhedron(M, 1);
2295 Matrix_Free(M);
2296 addeliminatedparams_enum(EP, CT, CEq, MaxRays, nparam);
2297 Polyhedron_Free(CEq);
2298 Matrix_Free(CT);
2299
2300 return EP;
2301 }
2302
2303 static void enumerate_vd_add_ray(evalue *EP, Matrix *Rays, unsigned MaxRays)
2304 {
2305 if (value_notzero_p(EP->d))
2306 return;
2307
2308 assert(EP->x.p->type == partition);
2309 assert(EP->x.p->pos == EVALUE_DOMAIN(EP->x.p->arr[0])->Dimension);
2310 for (int i = 0; i < EP->x.p->size/2; ++i) {
2311 Polyhedron *D = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
2312 Polyhedron *N = DomainAddRays(D, Rays, MaxRays);
2313 EVALUE_SET_DOMAIN(EP->x.p->arr[2*i], N);
2314 Domain_Free(D);
2315 }
2316 }
2317
2318 static evalue* enumerate_line(Polyhedron *P,
2319 unsigned exist, unsigned nparam, barvinok_options *options)
2320 {
2321 if (P->NbBid == 0)
2322 return 0;
2323
2324 #ifdef DEBUG_ER
2325 fprintf(stderr, "\nER: Line\n");
2326 #endif /* DEBUG_ER */
2327
2328 int nvar = P->Dimension - exist - nparam;
2329 int i, j;
2330 for (i = 0; i < nparam; ++i)
2331 if (value_notzero_p(P->Ray[0][1+nvar+exist+i]))
2332 break;
2333 assert(i < nparam);
2334 for (j = i+1; j < nparam; ++j)
2335 if (value_notzero_p(P->Ray[0][1+nvar+exist+i]))
2336 break;
2337 assert(j >= nparam); // for now
2338
2339 Matrix *M = Matrix_Alloc(2, P->Dimension+2);
2340 value_set_si(M->p[0][0], 1);
2341 value_set_si(M->p[0][1+nvar+exist+i], 1);
2342 value_set_si(M->p[1][0], 1);
2343 value_set_si(M->p[1][1+nvar+exist+i], -1);
2344 value_absolute(M->p[1][1+P->Dimension], P->Ray[0][1+nvar+exist+i]);
2345 value_decrement(M->p[1][1+P->Dimension], M->p[1][1+P->Dimension]);
2346 Polyhedron *S = AddConstraints(M->p[0], 2, P, options->MaxRays);
2347 evalue *EP = barvinok_enumerate_e_with_options(S, exist, nparam, options);
2348 Polyhedron_Free(S);
2349 Matrix_Free(M);
2350
2351 return enumerate_cyclic(P, exist, nparam, EP, 0, i, options->MaxRays);
2352 }
2353
2354 static int single_param_pos(Polyhedron*P, unsigned exist, unsigned nparam,
2355 int r)
2356 {
2357 int nvar = P->Dimension - exist - nparam;
2358 if (First_Non_Zero(P->Ray[r]+1, nvar) != -1)
2359 return -1;
2360 int i = First_Non_Zero(P->Ray[r]+1+nvar+exist, nparam);
2361 if (i == -1)
2362 return -1;
2363 if (First_Non_Zero(P->Ray[r]+1+nvar+exist+1, nparam-i-1) != -1)
2364 return -1;
2365 return i;
2366 }
2367
2368 static evalue* enumerate_remove_ray(Polyhedron *P, int r,
2369 unsigned exist, unsigned nparam, barvinok_options *options)
2370 {
2371 #ifdef DEBUG_ER
2372 fprintf(stderr, "\nER: RedundantRay\n");
2373 #endif /* DEBUG_ER */
2374
2375 Value one;
2376 value_init(one);
2377 value_set_si(one, 1);
2378 int len = P->NbRays-1;
2379 Matrix *M = Matrix_Alloc(2 * len, P->Dimension+2);
2380 Vector_Copy(P->Ray[0], M->p[0], r * (P->Dimension+2));
2381 Vector_Copy(P->Ray[r+1], M->p[r], (len-r) * (P->Dimension+2));
2382 for (int j = 0; j < P->NbRays; ++j) {
2383 if (j == r)
2384 continue;
2385 Vector_Combine(P->Ray[j], P->Ray[r], M->p[len+j-(j>r)],
2386 one, P->Ray[j][P->Dimension+1], P->Dimension+2);
2387 }
2388
2389 P = Rays2Polyhedron(M, options->MaxRays);
2390 Matrix_Free(M);
2391 evalue *EP = barvinok_enumerate_e_with_options(P, exist, nparam, options);
2392 Polyhedron_Free(P);
2393 value_clear(one);
2394
2395 return EP;
2396 }
2397
2398 static evalue* enumerate_redundant_ray(Polyhedron *P,
2399 unsigned exist, unsigned nparam, barvinok_options *options)
2400 {
2401 assert(P->NbBid == 0);
2402 int nvar = P->Dimension - exist - nparam;
2403 Value m;
2404 value_init(m);
2405
2406 for (int r = 0; r < P->NbRays; ++r) {
2407 if (value_notzero_p(P->Ray[r][P->Dimension+1]))
2408 continue;
2409 int i1 = single_param_pos(P, exist, nparam, r);
2410 if (i1 == -1)
2411 continue;
2412 for (int r2 = r+1; r2 < P->NbRays; ++r2) {
2413 if (value_notzero_p(P->Ray[r2][P->Dimension+1]))
2414 continue;
2415 int i2 = single_param_pos(P, exist, nparam, r2);
2416 if (i2 == -1)
2417 continue;
2418 if (i1 != i2)
2419 continue;
2420
2421 value_division(m, P->Ray[r][1+nvar+exist+i1],
2422 P->Ray[r2][1+nvar+exist+i1]);
2423 value_multiply(m, m, P->Ray[r2][1+nvar+exist+i1]);
2424 /* r2 divides r => r redundant */
2425 if (value_eq(m, P->Ray[r][1+nvar+exist+i1])) {
2426 value_clear(m);
2427 return enumerate_remove_ray(P, r, exist, nparam, options);
2428 }
2429
2430 value_division(m, P->Ray[r2][1+nvar+exist+i1],
2431 P->Ray[r][1+nvar+exist+i1]);
2432 value_multiply(m, m, P->Ray[r][1+nvar+exist+i1]);
2433 /* r divides r2 => r2 redundant */
2434 if (value_eq(m, P->Ray[r2][1+nvar+exist+i1])) {
2435 value_clear(m);
2436 return enumerate_remove_ray(P, r2, exist, nparam, options);
2437 }
2438 }
2439 }
2440 value_clear(m);
2441 return 0;
2442 }
2443
2444 static Polyhedron *upper_bound(Polyhedron *P,
2445 int pos, Value *max, Polyhedron **R)
2446 {
2447 Value v;
2448 int r;
2449 value_init(v);
2450
2451 *R = 0;
2452 Polyhedron *N;
2453 Polyhedron *B = 0;
2454 for (Polyhedron *Q = P; Q; Q = N) {
2455 N = Q->next;
2456 for (r = 0; r < P->NbRays; ++r) {
2457 if (value_zero_p(P->Ray[r][P->Dimension+1]) &&
2458 value_pos_p(P->Ray[r][1+pos]))
2459 break;
2460 }
2461 if (r < P->NbRays) {
2462 Q->next = *R;
2463 *R = Q;
2464 continue;
2465 } else {
2466 Q->next = B;
2467 B = Q;
2468 }
2469 for (r = 0; r < P->NbRays; ++r) {
2470 if (value_zero_p(P->Ray[r][P->Dimension+1]))
2471 continue;
2472 mpz_fdiv_q(v, P->Ray[r][1+pos], P->Ray[r][1+P->Dimension]);
2473 if ((!Q->next && r == 0) || value_gt(v, *max))
2474 value_assign(*max, v);
2475 }
2476 }
2477 value_clear(v);
2478 return B;
2479 }
2480
2481 static evalue* enumerate_ray(Polyhedron *P,
2482 unsigned exist, unsigned nparam, barvinok_options *options)
2483 {
2484 assert(P->NbBid == 0);
2485 int nvar = P->Dimension - exist - nparam;
2486
2487 int r;
2488 for (r = 0; r < P->NbRays; ++r)
2489 if (value_zero_p(P->Ray[r][P->Dimension+1]))
2490 break;
2491 if (r >= P->NbRays)
2492 return 0;
2493
2494 int r2;
2495 for (r2 = r+1; r2 < P->NbRays; ++r2)
2496 if (value_zero_p(P->Ray[r2][P->Dimension+1]))
2497 break;
2498 if (r2 < P->NbRays) {
2499 if (nvar > 0)
2500 return enumerate_sum(P, exist, nparam, options);
2501 }
2502
2503 #ifdef DEBUG_ER
2504 fprintf(stderr, "\nER: Ray\n");
2505 #endif /* DEBUG_ER */
2506
2507 Value m;
2508 Value one;
2509 value_init(m);
2510 value_init(one);
2511 value_set_si(one, 1);
2512 int i = single_param_pos(P, exist, nparam, r);
2513 assert(i != -1); // for now;
2514
2515 Matrix *M = Matrix_Alloc(P->NbRays, P->Dimension+2);
2516 for (int j = 0; j < P->NbRays; ++j) {
2517 Vector_Combine(P->Ray[j], P->Ray[r], M->p[j],
2518 one, P->Ray[j][P->Dimension+1], P->Dimension+2);
2519 }
2520 Polyhedron *S = Rays2Polyhedron(M, options->MaxRays);
2521 Matrix_Free(M);
2522 Polyhedron *D = DomainDifference(P, S, options->MaxRays);
2523 Polyhedron_Free(S);
2524 // Polyhedron_Print(stderr, P_VALUE_FMT, D);
2525 assert(value_pos_p(P->Ray[r][1+nvar+exist+i])); // for now
2526 Polyhedron *R;
2527 D = upper_bound(D, nvar+exist+i, &m, &R);
2528 assert(D);
2529 Domain_Free(D);
2530
2531 M = Matrix_Alloc(2, P->Dimension+2);
2532 value_set_si(M->p[0][0], 1);
2533 value_set_si(M->p[1][0], 1);
2534 value_set_si(M->p[0][1+nvar+exist+i], -1);
2535 value_set_si(M->p[1][1+nvar+exist+i], 1);
2536 value_assign(M->p[0][1+P->Dimension], m);
2537 value_oppose(M->p[1][1+P->Dimension], m);
2538 value_addto(M->p[1][1+P->Dimension], M->p[1][1+P->Dimension],
2539 P->Ray[r][1+nvar+exist+i]);
2540 value_decrement(M->p[1][1+P->Dimension], M->p[1][1+P->Dimension]);
2541 // Matrix_Print(stderr, P_VALUE_FMT, M);
2542 D = AddConstraints(M->p[0], 2, P, options->MaxRays);
2543 // Polyhedron_Print(stderr, P_VALUE_FMT, D);
2544 value_subtract(M->p[0][1+P->Dimension], M->p[0][1+P->Dimension],
2545 P->Ray[r][1+nvar+exist+i]);
2546 // Matrix_Print(stderr, P_VALUE_FMT, M);
2547 S = AddConstraints(M->p[0], 1, P, options->MaxRays);
2548 // Polyhedron_Print(stderr, P_VALUE_FMT, S);
2549 Matrix_Free(M);
2550
2551 evalue *EP = barvinok_enumerate_e_with_options(D, exist, nparam, options);
2552 Polyhedron_Free(D);
2553 value_clear(one);
2554 value_clear(m);
2555
2556 if (value_notone_p(P->Ray[r][1+nvar+exist+i]))
2557 EP = enumerate_cyclic(P, exist, nparam, EP, r, i, options->MaxRays);
2558 else {
2559 M = Matrix_Alloc(1, nparam+2);
2560 value_set_si(M->p[0][0], 1);
2561 value_set_si(M->p[0][1+i], 1);
2562 enumerate_vd_add_ray(EP, M, options->MaxRays);
2563 Matrix_Free(M);
2564 }
2565
2566 if (!emptyQ(S)) {
2567 evalue *E = barvinok_enumerate_e_with_options(S, exist, nparam, options);
2568 eadd(E, EP);
2569 free_evalue_refs(E);
2570 free(E);
2571 }
2572 Polyhedron_Free(S);
2573
2574 if (R) {
2575 assert(nvar == 0);
2576 evalue *ER = enumerate_or(R, exist, nparam, options);
2577 eor(ER, EP);
2578 free_evalue_refs(ER);
2579 free(ER);
2580 }
2581
2582 return EP;
2583 }
2584
2585 static evalue* enumerate_vd(Polyhedron **PA,
2586 unsigned exist, unsigned nparam, barvinok_options *options)
2587 {
2588 Polyhedron *P = *PA;
2589 int nvar = P->Dimension - exist - nparam;
2590 Param_Polyhedron *PP = NULL;
2591 Polyhedron *C = Universe_Polyhedron(nparam);
2592 Polyhedron *CEq;
2593 Matrix *CT;
2594 Polyhedron *PR = P;
2595 PP = Polyhedron2Param_Domain(PR,C, options->MaxRays);
2596 Polyhedron_Free(C);
2597
2598 int nd;
2599 Param_Domain *D, *last;
2600 Value c;
2601 value_init(c);
2602 for (nd = 0, D=PP->D; D; D=D->next, ++nd)
2603 ;
2604
2605 Polyhedron **VD = new Polyhedron_p[nd];
2606 Polyhedron *TC = true_context(P, C, options->MaxRays);
2607 FORALL_REDUCED_DOMAIN(PP, TC, nd, options, i, D, rVD)
2608 VD[nd++] = rVD;
2609 last = D;
2610 END_FORALL_REDUCED_DOMAIN
2611 Polyhedron_Free(TC);
2612
2613 evalue *EP = 0;
2614
2615 if (nd == 0)
2616 EP = evalue_zero();
2617
2618 /* This doesn't seem to have any effect */
2619 if (nd == 1) {
2620 Polyhedron *CA = align_context(VD[0], P->Dimension, options->MaxRays);
2621 Polyhedron *O = P;
2622 P = DomainIntersection(P, CA, options->MaxRays);
2623 if (O != *PA)
2624 Polyhedron_Free(O);
2625 Polyhedron_Free(CA);
2626 if (emptyQ(P))
2627 EP = evalue_zero();
2628 }
2629
2630 if (PR != *PA)
2631 Polyhedron_Free(PR);
2632 PR = 0;
2633
2634 if (!EP && nd > 1) {
2635 #ifdef DEBUG_ER
2636 fprintf(stderr, "\nER: VD\n");
2637 #endif /* DEBUG_ER */
2638 for (int i = 0; i < nd; ++i) {
2639 Polyhedron *CA = align_context(VD[i], P->Dimension, options->MaxRays);
2640 Polyhedron *I = DomainIntersection(P, CA, options->MaxRays);
2641
2642 if (i == 0)
2643 EP = barvinok_enumerate_e_with_options(I, exist, nparam, options);
2644 else {
2645 evalue *E = barvinok_enumerate_e_with_options(I, exist, nparam,
2646 options);
2647 eadd(E, EP);
2648 free_evalue_refs(E);
2649 free(E);
2650 }
2651 Polyhedron_Free(I);
2652 Polyhedron_Free(CA);
2653 }
2654 }
2655
2656 for (int i = 0; i < nd; ++i)
2657 Polyhedron_Free(VD[i]);
2658 delete [] VD;
2659 value_clear(c);
2660
2661 if (!EP && nvar == 0) {
2662 Value f;
2663 value_init(f);
2664 Param_Vertices *V, *V2;
2665 Matrix* M = Matrix_Alloc(1, P->Dimension+2);
2666
2667 FORALL_PVertex_in_ParamPolyhedron(V, last, PP) {
2668 bool found = false;
2669 FORALL_PVertex_in_ParamPolyhedron(V2, last, PP) {
2670 if (V == V2) {
2671 found = true;
2672 continue;
2673 }
2674 if (!found)
2675 continue;
2676 for (int i = 0; i < exist; ++i) {
2677 value_oppose(f, V->Vertex->p[i][nparam+1]);
2678 Vector_Combine(V->Vertex->p[i],
2679 V2->Vertex->p[i],
2680 M->p[0] + 1 + nvar + exist,
2681 V2->Vertex->p[i][nparam+1],
2682 f,
2683 nparam+1);
2684 int j;
2685 for (j = 0; j < nparam; ++j)
2686 if (value_notzero_p(M->p[0][1+nvar+exist+j]))
2687 break;
2688 if (j >= nparam)
2689 continue;
2690 ConstraintSimplify(M->p[0], M->p[0],
2691 P->Dimension+2, &f);
2692 value_set_si(M->p[0][0], 0);
2693 Polyhedron *para = AddConstraints(M->p[0], 1, P,
2694 options->MaxRays);
2695 POL_ENSURE_VERTICES(para);
2696 if (emptyQ(para)) {
2697 Polyhedron_Free(para);
2698 continue;
2699 }
2700 Polyhedron *pos, *neg;
2701 value_set_si(M->p[0][0], 1);
2702 value_decrement(M->p[0][P->Dimension+1],
2703 M->p[0][P->Dimension+1]);
2704 neg = AddConstraints(M->p[0], 1, P, options->MaxRays);
2705 value_set_si(f, -1);
2706 Vector_Scale(M->p[0]+1, M->p[0]+1, f,
2707 P->Dimension+1);
2708 value_decrement(M->p[0][P->Dimension+1],
2709 M->p[0][P->Dimension+1]);
2710 value_decrement(M->p[0][P->Dimension+1],
2711 M->p[0][P->Dimension+1]);
2712 pos = AddConstraints(M->p[0], 1, P, options->MaxRays);
2713 POL_ENSURE_VERTICES(neg);
2714 POL_ENSURE_VERTICES(pos);
2715 if (emptyQ(neg) && emptyQ(pos)) {
2716 Polyhedron_Free(para);
2717 Polyhedron_Free(pos);
2718 Polyhedron_Free(neg);
2719 continue;
2720 }
2721 #ifdef DEBUG_ER
2722 fprintf(stderr, "\nER: Order\n");
2723 #endif /* DEBUG_ER */
2724 EP = barvinok_enumerate_e_with_options(para, exist, nparam,
2725 options);
2726 evalue *E;
2727 if (!emptyQ(pos)) {
2728 E = barvinok_enumerate_e_with_options(pos, exist, nparam,
2729 options);
2730 eadd(E, EP);
2731 free_evalue_refs(E);
2732 free(E);
2733 }
2734 if (!emptyQ(neg)) {
2735 E = barvinok_enumerate_e_with_options(neg, exist, nparam,
2736 options);
2737 eadd(E, EP);
2738 free_evalue_refs(E);
2739 free(E);
2740 }
2741 Polyhedron_Free(para);
2742 Polyhedron_Free(pos);
2743 Polyhedron_Free(neg);
2744 break;
2745 }
2746 if (EP)
2747 break;
2748 } END_FORALL_PVertex_in_ParamPolyhedron;
2749 if (EP)
2750 break;
2751 } END_FORALL_PVertex_in_ParamPolyhedron;
2752
2753 if (!EP) {
2754 /* Search for vertex coordinate to split on */
2755 /* First look for one independent of the parameters */
2756 FORALL_PVertex_in_ParamPolyhedron(V, last, PP) {
2757 for (int i = 0; i < exist; ++i) {
2758 int j;
2759 for (j = 0; j < nparam; ++j)
2760 if (value_notzero_p(V->Vertex->p[i][j]))
2761 break;
2762 if (j < nparam)
2763 continue;
2764 value_set_si(M->p[0][0], 1);
2765 Vector_Set(M->p[0]+1, 0, nvar+exist);
2766 Vector_Copy(V->Vertex->p[i],
2767 M->p[0] + 1 + nvar + exist, nparam+1);
2768 value_oppose(M->p[0][1+nvar+i],
2769 V->Vertex->p[i][nparam+1]);
2770
2771 Polyhedron *pos, *neg;
2772 value_set_si(M->p[0][0], 1);
2773 value_decrement(M->p[0][P->Dimension+1],
2774 M->p[0][P->Dimension+1]);
2775 neg = AddConstraints(M->p[0], 1, P, options->MaxRays);
2776 value_set_si(f, -1);
2777 Vector_Scale(M->p[0]+1, M->p[0]+1, f,
2778 P->Dimension+1);
2779 value_decrement(M->p[0][P->Dimension+1],
2780 M->p[0][P->Dimension+1]);
2781 value_decrement(M->p[0][P->Dimension+1],
2782 M->p[0][P->Dimension+1]);
2783 pos = AddConstraints(M->p[0], 1, P, options->MaxRays);
2784 POL_ENSURE_VERTICES(neg);
2785 POL_ENSURE_VERTICES(pos);
2786 if (emptyQ(neg) || emptyQ(pos)) {
2787 Polyhedron_Free(pos);
2788 Polyhedron_Free(neg);
2789 continue;
2790 }
2791 Polyhedron_Free(pos);
2792 value_increment(M->p[0][P->Dimension+1],
2793 M->p[0][P->Dimension+1]);
2794 pos = AddConstraints(M->p[0], 1, P, options->MaxRays);
2795 #ifdef DEBUG_ER
2796 fprintf(stderr, "\nER: Vertex\n");
2797 #endif /* DEBUG_ER */
2798 pos->next = neg;
2799 EP = enumerate_or(pos, exist, nparam, options);
2800 break;
2801 }
2802 if (EP)
2803 break;
2804 } END_FORALL_PVertex_in_ParamPolyhedron;
2805 }
2806
2807 if (!EP) {
2808 /* Search for vertex coordinate to split on */
2809 /* Now look for one that depends on the parameters */
2810 FORALL_PVertex_in_ParamPolyhedron(V, last, PP) {
2811 for (int i = 0; i < exist; ++i) {
2812 value_set_si(M->p[0][0], 1);
2813 Vector_Set(M->p[0]+1, 0, nvar+exist);
2814 Vector_Copy(V->Vertex->p[i],
2815 M->p[0] + 1 + nvar + exist, nparam+1);
2816 value_oppose(M->p[0][1+nvar+i],
2817 V->Vertex->p[i][nparam+1]);
2818
2819 Polyhedron *pos, *neg;
2820 value_set_si(M->p[0][0], 1);
2821 value_decrement(M->p[0][P->Dimension+1],
2822 M->p[0][P->Dimension+1]);
2823 neg = AddConstraints(M->p[0], 1, P, options->MaxRays);
2824 value_set_si(f, -1);
2825 Vector_Scale(M->p[0]+1, M->p[0]+1, f,
2826 P->Dimension+1);
2827 value_decrement(M->p[0][P->Dimension+1],
2828 M->p[0][P->Dimension+1]);
2829 value_decrement(M->p[0][P->Dimension+1],
2830 M->p[0][P->Dimension+1]);
2831 pos = AddConstraints(M->p[0], 1, P, options->MaxRays);
2832 POL_ENSURE_VERTICES(neg);
2833 POL_ENSURE_VERTICES(pos);
2834 if (emptyQ(neg) || emptyQ(pos)) {
2835 Polyhedron_Free(pos);
2836 Polyhedron_Free(neg);
2837 continue;
2838 }
2839 Polyhedron_Free(pos);
2840 value_increment(M->p[0][P->Dimension+1],
2841 M->p[0][P->Dimension+1]);
2842 pos = AddConstraints(M->p[0], 1, P, options->MaxRays);
2843 #ifdef DEBUG_ER
2844 fprintf(stderr, "\nER: ParamVertex\n");
2845 #endif /* DEBUG_ER */
2846 pos->next = neg;
2847 EP = enumerate_or(pos, exist, nparam, options);
2848 break;
2849 }
2850 if (EP)
2851 break;
2852 } END_FORALL_PVertex_in_ParamPolyhedron;
2853 }
2854
2855 Matrix_Free(M);
2856 value_clear(f);
2857 }
2858
2859 if (CEq)
2860 Polyhedron_Free(CEq);
2861 if (CT)
2862 Matrix_Free(CT);
2863 if (PP)
2864 Param_Polyhedron_Free(PP);
2865 *PA = P;
2866
2867 return EP;
2868 }
2869
2870 evalue* barvinok_enumerate_pip(Polyhedron *P, unsigned exist, unsigned nparam,
2871 unsigned MaxRays)
2872 {
2873 evalue *E;
2874 barvinok_options *options = barvinok_options_new_with_defaults();
2875 options->MaxRays = MaxRays;
2876 E = barvinok_enumerate_pip_with_options(P, exist, nparam, options);
2877 barvinok_options_free(options);
2878 return E;
2879 }
2880
2881 #ifndef HAVE_PIPLIB
2882 evalue *barvinok_enumerate_pip_with_options(Polyhedron *P,
2883 unsigned exist, unsigned nparam, struct barvinok_options *options)
2884 {
2885 return 0;
2886 }
2887 #else
2888 evalue *barvinok_enumerate_pip_with_options(Polyhedron *P,
2889 unsigned exist, unsigned nparam, struct barvinok_options *options)
2890 {
2891 int nvar = P->Dimension - exist - nparam;
2892 evalue *EP = evalue_zero();
2893 Polyhedron *Q, *N;
2894
2895 #ifdef DEBUG_ER
2896 fprintf(stderr, "\nER: PIP\n");
2897 #endif /* DEBUG_ER */
2898
2899 Polyhedron *D = pip_projectout(P, nvar, exist, nparam);
2900 for (Q = D; Q; Q = N) {
2901 N = Q->next;
2902 Q->next = 0;
2903 evalue *E;
2904 exist = Q->Dimension - nvar - nparam;
2905 E = barvinok_enumerate_e_with_options(Q, exist, nparam, options);
2906 Polyhedron_Free(Q);
2907 eadd(E, EP);
2908 free_evalue_refs(E);
2909 free(E);
2910 }
2911
2912 return EP;
2913 }
2914 #endif
2915
2916
2917 static bool is_single(Value *row, int pos, int len)
2918 {
2919 return First_Non_Zero(row, pos) == -1 &&
2920 First_Non_Zero(row+pos+1, len-pos-1) == -1;
2921 }
2922
2923 static evalue* barvinok_enumerate_e_r(Polyhedron *P,
2924 unsigned exist, unsigned nparam, barvinok_options *options);
2925
2926 #ifdef DEBUG_ER
2927 static int er_level = 0;
2928
2929 evalue* barvinok_enumerate_e_with_options(Polyhedron *P,
2930 unsigned exist, unsigned nparam, barvinok_options *options)
2931 {
2932 fprintf(stderr, "\nER: level %i\n", er_level);
2933
2934 Polyhedron_PrintConstraints(stderr, P_VALUE_FMT, P);
2935 fprintf(stderr, "\nE %d\nP %d\n", exist, nparam);
2936 ++er_level;
2937 P = DomainConstraintSimplify(Polyhedron_Copy(P), options->MaxRays);
2938 evalue *EP = barvinok_enumerate_e_r(P, exist, nparam, options);
2939 Polyhedron_Free(P);
2940 --er_level;
2941 return EP;
2942 }
2943 #else
2944 evalue* barvinok_enumerate_e_with_options(Polyhedron *P,
2945 unsigned exist, unsigned nparam, barvinok_options *options)
2946 {
2947 P = DomainConstraintSimplify(Polyhedron_Copy(P), options->MaxRays);
2948 evalue *EP = barvinok_enumerate_e_r(P, exist, nparam, options);
2949 Polyhedron_Free(P);
2950 return EP;
2951 }
2952 #endif
2953
2954 evalue* barvinok_enumerate_e(Polyhedron *P, unsigned exist, unsigned nparam,
2955 unsigned MaxRays)
2956 {
2957 evalue *E;
2958 barvinok_options *options = barvinok_options_new_with_defaults();
2959 options->MaxRays = MaxRays;
2960 E = barvinok_enumerate_e_with_options(P, exist, nparam, options);
2961 barvinok_options_free(options);
2962 return E;
2963 }
2964
2965 static evalue* barvinok_enumerate_e_r(Polyhedron *P,
2966 unsigned exist, unsigned nparam, barvinok_options *options)
2967 {
2968 if (exist == 0) {
2969 Polyhedron *U = Universe_Polyhedron(nparam);
2970 evalue *EP = barvinok_enumerate_with_options(P, U, options);
2971 //char *param_name[] = {"P", "Q", "R", "S", "T" };
2972 //print_evalue(stdout, EP, param_name);
2973 Polyhedron_Free(U);
2974 return EP;
2975 }
2976
2977 int nvar = P->Dimension - exist - nparam;
2978 int len = P->Dimension + 2;
2979
2980 /* for now */
2981 POL_ENSURE_FACETS(P);
2982 POL_ENSURE_VERTICES(P);
2983
2984 if (emptyQ(P))
2985 return evalue_zero();
2986
2987 if (nvar == 0 && nparam == 0) {
2988 evalue *EP = evalue_zero();
2989 barvinok_count_with_options(P, &EP->x.n, options);
2990 if (value_pos_p(EP->x.n))
2991 value_set_si(EP->x.n, 1);
2992 return EP;
2993 }
2994
2995 int r;
2996 for (r = 0; r < P->NbRays; ++r)
2997 if (value_zero_p(P->Ray[r][0]) ||
2998 value_zero_p(P->Ray[r][P->Dimension+1])) {
2999 int i;
3000 for (i = 0; i < nvar; ++i)
3001 if (value_notzero_p(P->Ray[r][i+1]))
3002 break;
3003 if (i >= nvar)
3004 continue;
3005 for (i = nvar + exist; i < nvar + exist + nparam; ++i)
3006 if (value_notzero_p(P->Ray[r][i+1]))
3007 break;
3008 if (i >= nvar + exist + nparam)
3009 break;
3010 }
3011 if (r < P->NbRays) {
3012 evalue *EP = evalue_zero();
3013 value_set_si(EP->x.n, -1);
3014 return EP;
3015 }
3016
3017 int first;
3018 for (r = 0; r < P->NbEq; ++r)
3019 if ((first = First_Non_Zero(P->Constraint[r]+1+nvar, exist)) != -1)
3020 break;
3021 if (r < P->NbEq) {
3022 if (First_Non_Zero(P->Constraint[r]+1+nvar+first+1,
3023 exist-first-1) != -1) {
3024 Polyhedron *T = rotate_along(P, r, nvar, exist, options->MaxRays);
3025 #ifdef DEBUG_ER
3026 fprintf(stderr, "\nER: Equality\n");
3027 #endif /* DEBUG_ER */
3028 evalue *EP = barvinok_enumerate_e_with_options(T, exist-1, nparam,
3029 options);
3030 Polyhedron_Free(T);
3031 return EP;
3032 } else {
3033 #ifdef DEBUG_ER
3034 fprintf(stderr, "\nER: Fixed\n");
3035 #endif /* DEBUG_ER */
3036 if (first == 0)
3037 return barvinok_enumerate_e_with_options(P, exist-1, nparam,
3038 options);
3039 else {
3040 Polyhedron *T = Polyhedron_Copy(P);
3041 SwapColumns(T, nvar+1, nvar+1+first);
3042 evalue *EP = barvinok_enumerate_e_with_options(T, exist-1, nparam,
3043 options);
3044 Polyhedron_Free(T);
3045 return EP;
3046 }
3047 }
3048 }
3049
3050 Vector *row = Vector_Alloc(len);
3051 value_set_si(row->p[0], 1);
3052
3053 Value f;
3054 value_init(f);
3055
3056 enum constraint* info = new constraint[exist];
3057 for (int i = 0; i < exist; ++i) {
3058 info[i] = ALL_POS;
3059 for (int l = P->NbEq; l < P->NbConstraints; ++l) {
3060 if (value_negz_p(P->Constraint[l][nvar+i+1]))
3061 continue;
3062 bool l_parallel = is_single(P->Constraint[l]+nvar+1, i, exist);
3063 for (int u = P->NbEq; u < P->NbConstraints; ++u) {
3064 if (value_posz_p(P->Constraint[u][nvar+i+1]))
3065 continue;
3066 bool lu_parallel = l_parallel ||
3067 is_single(P->Constraint[u]+nvar+1, i, exist);
3068 value_oppose(f, P->Constraint[u][nvar+i+1]);
3069 Vector_Combine(P->Constraint[l]+1, P->Constraint[u]+1, row->p+1,
3070 f, P->Constraint[l][nvar+i+1], len-1);
3071 if (!(info[i] & INDEPENDENT)) {
3072 int j;
3073 for (j = 0; j < exist; ++j)
3074 if (j != i && value_notzero_p(row->p[nvar+j+1]))
3075 break;
3076 if (j == exist) {
3077 //printf("independent: i: %d, l: %d, u: %d\n", i, l, u);
3078 info[i] = (constraint)(info[i] | INDEPENDENT);
3079 }
3080 }
3081 if (info[i] & ALL_POS) {
3082 value_addto(row->p[len-1], row->p[len-1],
3083 P->Constraint[l][nvar+i+1]);
3084 value_addto(row->p[len-1], row->p[len-1], f);
3085 value_multiply(f, f, P->Constraint[l][nvar+i+1]);
3086 value_subtract(row->p[len-1], row->p[len-1], f);
3087 value_decrement(row->p[len-1], row->p[len-1]);
3088 ConstraintSimplify(row->p, row->p, len, &f);
3089 value_set_si(f, -1);
3090 Vector_Scale(row->p+1, row->p+1, f, len-1);
3091 value_decrement(row->p[len-1], row->p[len-1]);
3092 Polyhedron *T = AddConstraints(row->p, 1, P, options->MaxRays);
3093 POL_ENSURE_VERTICES(T);
3094 if (!emptyQ(T)) {
3095 //printf("not all_pos: i: %d, l: %d, u: %d\n", i, l, u);
3096 info[i] = (constraint)(info[i] ^ ALL_POS);
3097 }
3098 //puts("pos remainder");
3099 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
3100 Polyhedron_Free(T);
3101 }
3102 if (!(info[i] & ONE_NEG)) {
3103 if (lu_parallel) {
3104 negative_test_constraint(P->Constraint[l],
3105 P->Constraint[u],
3106 row->p, nvar+i, len, &f);
3107 oppose_constraint(row->p, len, &f);
3108 Polyhedron *T = AddConstraints(row->p, 1, P,
3109 options->MaxRays);
3110 POL_ENSURE_VERTICES(T);
3111 if (emptyQ(T)) {
3112 //printf("one_neg i: %d, l: %d, u: %d\n", i, l, u);
3113 info[i] = (constraint)(info[i] | ONE_NEG);
3114 }
3115 //puts("neg remainder");
3116 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
3117 Polyhedron_Free(T);
3118 } else if (!(info[i] & ROT_NEG)) {
3119 if (parallel_constraints(P->Constraint[l],
3120 P->Constraint[u],
3121 row->p, nvar, exist)) {
3122 negative_test_constraint7(P->Constraint[l],
3123 P->Constraint[u],
3124 row->p, nvar, exist,
3125 len, &f);
3126 oppose_constraint(row->p, len, &f);
3127 Polyhedron *T = AddConstraints(row->p, 1, P,
3128 options->MaxRays);
3129 POL_ENSURE_VERTICES(T);
3130 if (emptyQ(T)) {
3131 // printf("rot_neg i: %d, l: %d, u: %d\n", i, l, u);
3132 info[i] = (constraint)(info[i] | ROT_NEG);
3133 r = l;
3134 }
3135 //puts("neg remainder");
3136 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
3137 Polyhedron_Free(T);
3138 }
3139 }
3140 }
3141 if (!(info[i] & ALL_POS) && (info[i] & (ONE_NEG | ROT_NEG)))
3142 goto next;
3143 }
3144 }
3145 if (info[i] & ALL_POS)
3146 break;
3147 next:
3148 ;
3149 }
3150
3151 /*
3152 for (int i = 0; i < exist; ++i)
3153 printf("%i: %i\n", i, info[i]);
3154 */
3155 for (int i = 0; i < exist; ++i)
3156 if (info[i] & ALL_POS) {
3157 #ifdef DEBUG_ER
3158 fprintf(stderr, "\nER: Positive\n");
3159 #endif /* DEBUG_ER */
3160 // Eliminate
3161 // Maybe we should chew off some of the fat here
3162 Matrix *M = Matrix_Alloc(P->Dimension, P->Dimension+1);
3163 for (int j = 0; j < P->Dimension; ++j)
3164 value_set_si(M->p[j][j + (j >= i+nvar)], 1);
3165 Polyhedron *T = Polyhedron_Image(P, M, options->MaxRays);
3166 Matrix_Free(M);
3167 evalue *EP = barvinok_enumerate_e_with_options(T, exist-1, nparam,
3168 options);
3169 Polyhedron_Free(T);
3170 value_clear(f);
3171 Vector_Free(row);
3172 delete [] info;
3173 return EP;
3174 }
3175 for (int i = 0; i < exist; ++i)
3176 if (info[i] & ONE_NEG) {
3177 #ifdef DEBUG_ER
3178 fprintf(stderr, "\nER: Negative\n");
3179 #endif /* DEBUG_ER */
3180 Vector_Free(row);
3181 value_clear(f);
3182 delete [] info;
3183 if (i == 0)
3184 return barvinok_enumerate_e_with_options(P, exist-1, nparam,
3185 options);
3186 else {
3187 Polyhedron *T = Polyhedron_Copy(P);
3188 SwapColumns(T, nvar+1, nvar+1+i);
3189 evalue *EP = barvinok_enumerate_e_with_options(T, exist-1, nparam,
3190 options);
3191 Polyhedron_Free(T);
3192 return EP;
3193 }
3194 }
3195 for (int i = 0; i < exist; ++i)
3196 if (info[i] & ROT_NEG) {
3197 #ifdef DEBUG_ER
3198 fprintf(stderr, "\nER: Rotate\n");
3199 #endif /* DEBUG_ER */
3200 Vector_Free(row);
3201 value_clear(f);
3202 delete [] info;
3203 Polyhedron *T = rotate_along(P, r, nvar, exist, options->MaxRays);
3204 evalue *EP = barvinok_enumerate_e_with_options(T, exist-1, nparam,
3205 options);
3206 Polyhedron_Free(T);
3207 return EP;
3208 }
3209 for (int i = 0; i < exist; ++i)
3210 if (info[i] & INDEPENDENT) {
3211 Polyhedron *pos, *neg;
3212
3213 /* Find constraint again and split off negative part */
3214
3215 if (SplitOnVar(P, i, nvar, exist, options->MaxRays,
3216 row, f, true, &pos, &neg)) {
3217 #ifdef DEBUG_ER
3218 fprintf(stderr, "\nER: Split\n");
3219 #endif /* DEBUG_ER */
3220
3221 evalue *EP =
3222 barvinok_enumerate_e_with_options(neg, exist-1, nparam, options);
3223 evalue *E =
3224 barvinok_enumerate_e_with_options(pos, exist, nparam, options);
3225 eadd(E, EP);
3226 free_evalue_refs(E);
3227 free(E);
3228 Polyhedron_Free(neg);
3229 Polyhedron_Free(pos);
3230 value_clear(f);
3231 Vector_Free(row);
3232 delete [] info;
3233 return EP;
3234 }
3235 }
3236 delete [] info;
3237
3238 Polyhedron *O = P;
3239 Polyhedron *F;
3240
3241 evalue *EP;
3242
3243 EP = enumerate_line(P, exist, nparam, options);
3244 if (EP)
3245 goto out;
3246
3247 EP = barvinok_enumerate_pip_with_options(P, exist, nparam, options);
3248 if (EP)
3249 goto out;
3250
3251 EP = enumerate_redundant_ray(P, exist, nparam, options);
3252 if (EP)
3253 goto out;
3254
3255 EP = enumerate_sure(P, exist, nparam, options);
3256 if (EP)
3257 goto out;
3258
3259 EP = enumerate_ray(P, exist, nparam, options);
3260 if (EP)
3261 goto out;
3262
3263 EP = enumerate_sure2(P, exist, nparam, options);
3264 if (EP)
3265 goto out;
3266
3267 F = unfringe(P, options->MaxRays);
3268 if (!PolyhedronIncludes(F, P)) {
3269 #ifdef DEBUG_ER
3270 fprintf(stderr, "\nER: Fringed\n");
3271 #endif /* DEBUG_ER */
3272 EP = barvinok_enumerate_e_with_options(F, exist, nparam, options);
3273 Polyhedron_Free(F);
3274 goto out;
3275 }
3276 Polyhedron_Free(F);
3277
3278 if (nparam)
3279 EP = enumerate_vd(&P, exist, nparam, options);
3280 if (EP)
3281 goto out2;
3282
3283 if (nvar != 0) {
3284 EP = enumerate_sum(P, exist, nparam, options);
3285 goto out2;
3286 }
3287
3288 assert(nvar == 0);
3289
3290 int i;
3291 Polyhedron *pos, *neg;
3292 for (i = 0; i < exist; ++i)
3293 if (SplitOnVar(P, i, nvar, exist, options->MaxRays,
3294 row, f, false, &pos, &neg))
3295 break;
3296
3297 assert (i < exist);
3298
3299 pos->next = neg;
3300 EP = enumerate_or(pos, exist, nparam, options);
3301
3302 out2:
3303 if (O != P)
3304 Polyhedron_Free(P);
3305
3306 out:
3307 value_clear(f);
3308 Vector_Free(row);
3309 return EP;
3310 }
3311
3312 /*
3313 * remove equalities that require a "compression" of the parameters
3314 */
3315 static Polyhedron *remove_more_equalities(Polyhedron *P, unsigned nparam,
3316 Matrix **CP, unsigned MaxRays)
3317 {
3318 Polyhedron *Q = P;
3319 remove_all_equalities(&P, NULL, CP, NULL, nparam, MaxRays);
3320 if (P != Q)
3321 Polyhedron_Free(Q);
3322 return P;
3323 }
3324
3325 /* frees P */
3326 static gen_fun *series(Polyhedron *P, unsigned nparam, barvinok_options *options)
3327 {
3328 Matrix *CP = NULL;
3329 gen_fun *gf;
3330
3331 if (emptyQ2(P)) {
3332 Polyhedron_Free(P);
3333 return new gen_fun;
3334 }
3335
3336 assert(!Polyhedron_is_unbounded(P, nparam, options->MaxRays));
3337 assert(P->NbBid == 0);
3338 assert(Polyhedron_has_revlex_positive_rays(P, nparam));
3339 if (P->NbEq != 0)
3340 P = remove_more_equalities(P, nparam, &CP, options->MaxRays);
3341 assert(P->NbEq == 0);
3342 if (CP)
3343 nparam = CP->NbColumns-1;
3344
3345 if (nparam == 0) {
3346 Value c;
3347 value_init(c);
3348 barvinok_count_with_options(P, &c, options);
3349 gf = new gen_fun(c);
3350 value_clear(c);
3351 } else {
3352 gf_base *red;
3353 red = gf_base::create(Polyhedron_Project(P, nparam),
3354 P->Dimension, nparam, options);
3355 POL_ENSURE_VERTICES(P);
3356 red->start_gf(P, options);
3357 gf = red->gf;
3358 delete red;
3359 }
3360 if (CP) {
3361 gf->substitute(CP);
3362 Matrix_Free(CP);
3363 }
3364 Polyhedron_Free(P);
3365 return gf;
3366 }
3367
3368 gen_fun * barvinok_series_with_options(Polyhedron *P, Polyhedron* C,
3369 barvinok_options *options)
3370 {
3371 Polyhedron *CA;
3372 unsigned nparam = C->Dimension;
3373 gen_fun *gf;
3374
3375 CA = align_context(C, P->Dimension, options->MaxRays);
3376 P = DomainIntersection(P, CA, options->MaxRays);
3377 Polyhedron_Free(CA);
3378
3379 gf = series(P, nparam, options);
3380
3381 return gf;
3382 }
3383
3384 gen_fun * barvinok_series(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
3385 {
3386 gen_fun *gf;
3387 barvinok_options *options = barvinok_options_new_with_defaults();
3388 options->MaxRays = MaxRays;
3389 gf = barvinok_series_with_options(P, C, options);
3390 barvinok_options_free(options);
3391 return gf;
3392 }
3393
3394 static Polyhedron *skew_into_positive_orthant(Polyhedron *D, unsigned nparam,
3395 unsigned MaxRays)
3396 {
3397 Matrix *M = NULL;
3398 Value tmp;
3399 value_init(tmp);
3400 for (Polyhedron *P = D; P; P = P->next) {
3401 POL_ENSURE_VERTICES(P);
3402 assert(!Polyhedron_is_unbounded(P, nparam, MaxRays));
3403 assert(P->NbBid == 0);
3404 assert(Polyhedron_has_positive_rays(P, nparam));
3405
3406 for (int r = 0; r < P->NbRays; ++r) {
3407 if (value_notzero_p(P->Ray[r][P->Dimension+1]))
3408 continue;
3409 for (int i = 0; i < nparam; ++i) {
3410 int j;
3411 if (value_posz_p(P->Ray[r][i+1]))
3412 continue;
3413 if (!M) {
3414 M = Matrix_Alloc(D->Dimension+1, D->Dimension+1);
3415 for (int i = 0; i < D->Dimension+1; ++i)
3416 value_set_si(M->p[i][i], 1);
3417 } else {
3418 Inner_Product(P->Ray[r]+1, M->p[i], D->Dimension+1, &tmp);
3419 if (value_posz_p(tmp))
3420 continue;
3421 }
3422 for (j = P->Dimension - nparam; j < P->Dimension; ++j)
3423 if (value_pos_p(P->Ray[r][j+1]))
3424 break;
3425 assert(j < P->Dimension);
3426 value_pdivision(tmp, P->Ray[r][j+1], P->Ray[r][i+1]);
3427 value_subtract(M->p[i][j], M->p[i][j], tmp);
3428 }
3429 }
3430 }
3431 value_clear(tmp);
3432 if (M) {
3433 D = DomainImage(D, M, MaxRays);
3434 Matrix_Free(M);
3435 }
3436 return D;
3437 }
3438
3439 gen_fun* barvinok_enumerate_union_series_with_options(Polyhedron *D, Polyhedron* C,
3440 barvinok_options *options)
3441 {
3442 Polyhedron *conv, *D2;
3443 Polyhedron *CA;
3444 gen_fun *gf = NULL, *gf2;
3445 unsigned nparam = C->Dimension;
3446 ZZ one, mone;
3447 one = 1;
3448 mone = -1;
3449
3450 CA = align_context(C, D->Dimension, options->MaxRays);
3451 D = DomainIntersection(D, CA, options->MaxRays);
3452 Polyhedron_Free(CA);
3453
3454 D2 = skew_into_positive_orthant(D, nparam, options->MaxRays);
3455 for (Polyhedron *P = D2; P; P = P->next) {
3456 assert(P->Dimension == D2->Dimension);
3457 gen_fun *P_gf;
3458
3459 P_gf = series(Polyhedron_Copy(P), P->Dimension, options);
3460 if (!gf)
3461 gf = P_gf;
3462 else {
3463 gf->add_union(P_gf, options);
3464 delete P_gf;
3465 }
3466 }
3467 /* we actually only need the convex union of the parameter space
3468 * but the reducer classes currently expect a polyhedron in
3469 * the combined space
3470 */
3471 Polyhedron_Free(gf->context);
3472 gf->context = DomainConvex(D2, options->MaxRays);
3473
3474 gf2 = gf->summate(D2->Dimension - nparam, options);
3475
3476 delete gf;
3477 if (D != D2)
3478 Domain_Free(D2);
3479 Domain_Free(D);
3480 return gf2;
3481 }
3482
3483 gen_fun* barvinok_enumerate_union_series(Polyhedron *D, Polyhedron* C,
3484 unsigned MaxRays)
3485 {
3486 gen_fun *gf;
3487 barvinok_options *options = barvinok_options_new_with_defaults();
3488 options->MaxRays = MaxRays;
3489 gf = barvinok_enumerate_union_series_with_options(D, C, options);
3490 barvinok_options_free(options);
3491 return gf;
3492 }
3493
3494 evalue* barvinok_enumerate_union(Polyhedron *D, Polyhedron* C, unsigned MaxRays)
3495 {
3496 evalue *EP;
3497 gen_fun *gf = barvinok_enumerate_union_series(D, C, MaxRays);
3498 EP = *gf;
3499 delete gf;
3500 return EP;
3501 }
lattice point enumeration library