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