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