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