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