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