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