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