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