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