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