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