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