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