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