?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
extract check for infinite parametric polytopes
[barvinok.git] / barvinok.cc
blob5ff82874a12ad12460843abd897d35af92e2f8f6
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 <util.h>
11 extern "C" {
12 #include <polylib/polylibgmp.h>
13 #include "ev_operations.h"
14 #include "piputil.h"
15 }
16 #include "config.h"
17 #include <barvinok.h>
18
19 #ifdef NTL_STD_CXX
20 using namespace NTL;
21 #endif
22 using std::cout;
23 using std::endl;
24 using std::vector;
25 using std::deque;
26 using std::string;
27 using std::ostringstream;
28
29 #define ALLOC(p) (((long *) (p))[0])
30 #define SIZE(p) (((long *) (p))[1])
31 #define DATA(p) ((mp_limb_t *) (((long *) (p)) + 2))
32
33 static void value2zz(Value v, ZZ& z)
34 {
35 int sa = v[0]._mp_size;
36 int abs_sa = sa < 0 ? -sa : sa;
37
38 _ntl_gsetlength(&z.rep, abs_sa);
39 mp_limb_t * adata = DATA(z.rep);
40 for (int i = 0; i < abs_sa; ++i)
41 adata[i] = v[0]._mp_d[i];
42 SIZE(z.rep) = sa;
43 }
44
45 static void zz2value(ZZ& z, Value& v)
46 {
47 if (!z.rep) {
48 value_set_si(v, 0);
49 return;
50 }
51
52 int sa = SIZE(z.rep);
53 int abs_sa = sa < 0 ? -sa : sa;
54
55 mp_limb_t * adata = DATA(z.rep);
56 _mpz_realloc(v, abs_sa);
57 for (int i = 0; i < abs_sa; ++i)
58 v[0]._mp_d[i] = adata[i];
59 v[0]._mp_size = sa;
60 }
61
62 #undef ALLOC
63 #define ALLOC(t,p) p = (t*)malloc(sizeof(*p))
64
65 /*
66 * We just ignore the last column and row
67 * If the final element is not equal to one
68 * then the result will actually be a multiple of the input
69 */
70 static void matrix2zz(Matrix *M, mat_ZZ& m, unsigned nr, unsigned nc)
71 {
72 m.SetDims(nr, nc);
73
74 for (int i = 0; i < nr; ++i) {
75 // assert(value_one_p(M->p[i][M->NbColumns - 1]));
76 for (int j = 0; j < nc; ++j) {
77 value2zz(M->p[i][j], m[i][j]);
78 }
79 }
80 }
81
82 static void values2zz(Value *p, vec_ZZ& v, int len)
83 {
84 v.SetLength(len);
85
86 for (int i = 0; i < len; ++i) {
87 value2zz(p[i], v[i]);
88 }
89 }
90
91 /*
92 */
93 static void zz2values(vec_ZZ& v, Value *p)
94 {
95 for (int i = 0; i < v.length(); ++i)
96 zz2value(v[i], p[i]);
97 }
98
99 static void rays(mat_ZZ& r, Polyhedron *C)
100 {
101 unsigned dim = C->NbRays - 1; /* don't count zero vertex */
102 assert(C->NbRays - 1 == C->Dimension);
103 r.SetDims(dim, dim);
104 ZZ tmp;
105
106 int i, c;
107 for (i = 0, c = 0; i < dim; ++i)
108 if (value_zero_p(C->Ray[i][dim+1])) {
109 for (int j = 0; j < dim; ++j) {
110 value2zz(C->Ray[i][j+1], tmp);
111 r[j][c] = tmp;
112 }
113 ++c;
114 }
115 }
116
117 static Matrix * rays(Polyhedron *C)
118 {
119 unsigned dim = C->NbRays - 1; /* don't count zero vertex */
120 assert(C->NbRays - 1 == C->Dimension);
121
122 Matrix *M = Matrix_Alloc(dim+1, dim+1);
123 assert(M);
124
125 int i, c;
126 for (i = 0, c = 0; i <= dim && c < dim; ++i)
127 if (value_zero_p(C->Ray[i][dim+1])) {
128 Vector_Copy(C->Ray[i] + 1, M->p[c], dim);
129 value_set_si(M->p[c++][dim], 0);
130 }
131 assert(c == dim);
132 value_set_si(M->p[dim][dim], 1);
133
134 return M;
135 }
136
137 static Matrix * rays2(Polyhedron *C)
138 {
139 unsigned dim = C->NbRays - 1; /* don't count zero vertex */
140 assert(C->NbRays - 1 == C->Dimension);
141
142 Matrix *M = Matrix_Alloc(dim, dim);
143 assert(M);
144
145 int i, c;
146 for (i = 0, c = 0; i <= dim && c < dim; ++i)
147 if (value_zero_p(C->Ray[i][dim+1]))
148 Vector_Copy(C->Ray[i] + 1, M->p[c++], dim);
149 assert(c == dim);
150
151 return M;
152 }
153
154 /*
155 * Returns the largest absolute value in the vector
156 */
157 static ZZ max(vec_ZZ& v)
158 {
159 ZZ max = abs(v[0]);
160 for (int i = 1; i < v.length(); ++i)
161 if (abs(v[i]) > max)
162 max = abs(v[i]);
163 return max;
164 }
165
166 class cone {
167 public:
168 cone(Matrix *M) {
169 Cone = 0;
170 Rays = Matrix_Copy(M);
171 set_det();
172 }
173 cone(Polyhedron *C) {
174 Cone = Polyhedron_Copy(C);
175 Rays = rays(C);
176 set_det();
177 }
178 void set_det() {
179 mat_ZZ A;
180 matrix2zz(Rays, A, Rays->NbRows - 1, Rays->NbColumns - 1);
181 det = determinant(A);
182 Value v;
183 value_init(v);
184 zz2value(det, v);
185 value_clear(v);
186 }
187
188 Vector* short_vector(vec_ZZ& lambda) {
189 Matrix *M = Matrix_Copy(Rays);
190 Matrix *inv = Matrix_Alloc(M->NbRows, M->NbColumns);
191 int ok = Matrix_Inverse(M, inv);
192 assert(ok);
193 Matrix_Free(M);
194
195 ZZ det2;
196 mat_ZZ B;
197 mat_ZZ U;
198 matrix2zz(inv, B, inv->NbRows - 1, inv->NbColumns - 1);
199 long r = LLL(det2, B, U);
200
201 ZZ min = max(B[0]);
202 int index = 0;
203 for (int i = 1; i < B.NumRows(); ++i) {
204 ZZ tmp = max(B[i]);
205 if (tmp < min) {
206 min = tmp;
207 index = i;
208 }
209 }
210
211 Matrix_Free(inv);
212
213 lambda = B[index];
214
215 Vector *z = Vector_Alloc(U[index].length()+1);
216 assert(z);
217 zz2values(U[index], z->p);
218 value_set_si(z->p[U[index].length()], 0);
219
220 Value tmp;
221 value_init(tmp);
222 Polyhedron *C = poly();
223 int i;
224 for (i = 0; i < C->NbConstraints; ++i) {
225 Inner_Product(z->p, C->Constraint[i]+1, z->Size-1, &tmp);
226 if (value_pos_p(tmp))
227 break;
228 }
229 if (i == C->NbConstraints) {
230 value_set_si(tmp, -1);
231 Vector_Scale(z->p, z->p, tmp, z->Size-1);
232 }
233 value_clear(tmp);
234 return z;
235 }
236
237 ~cone() {
238 Polyhedron_Free(Cone);
239 Matrix_Free(Rays);
240 }
241
242 Polyhedron *poly() {
243 if (!Cone) {
244 Matrix *M = Matrix_Alloc(Rays->NbRows+1, Rays->NbColumns+1);
245 for (int i = 0; i < Rays->NbRows; ++i) {
246 Vector_Copy(Rays->p[i], M->p[i]+1, Rays->NbColumns);
247 value_set_si(M->p[i][0], 1);
248 }
249 Vector_Set(M->p[Rays->NbRows]+1, 0, Rays->NbColumns-1);
250 value_set_si(M->p[Rays->NbRows][0], 1);
251 value_set_si(M->p[Rays->NbRows][Rays->NbColumns], 1);
252 Cone = Rays2Polyhedron(M, M->NbRows+1);
253 assert(Cone->NbConstraints == Cone->NbRays);
254 Matrix_Free(M);
255 }
256 return Cone;
257 }
258
259 ZZ det;
260 Polyhedron *Cone;
261 Matrix *Rays;
262 };
263
264 class dpoly {
265 public:
266 vec_ZZ coeff;
267 dpoly(int d, ZZ& degree, int offset = 0) {
268 coeff.SetLength(d+1);
269
270 int min = d + offset;
271 if (degree >= 0 && degree < ZZ(INIT_VAL, min))
272 min = to_int(degree);
273
274 ZZ c = ZZ(INIT_VAL, 1);
275 if (!offset)
276 coeff[0] = c;
277 for (int i = 1; i <= min; ++i) {
278 c *= (degree -i + 1);
279 c /= i;
280 coeff[i-offset] = c;
281 }
282 }
283 void operator *= (dpoly& f) {
284 assert(coeff.length() == f.coeff.length());
285 vec_ZZ old = coeff;
286 coeff = f.coeff[0] * coeff;
287 for (int i = 1; i < coeff.length(); ++i)
288 for (int j = 0; i+j < coeff.length(); ++j)
289 coeff[i+j] += f.coeff[i] * old[j];
290 }
291 void div(dpoly& d, mpq_t count, ZZ& sign) {
292 int len = coeff.length();
293 Value tmp;
294 value_init(tmp);
295 mpq_t* c = new mpq_t[coeff.length()];
296 mpq_t qtmp;
297 mpq_init(qtmp);
298 for (int i = 0; i < len; ++i) {
299 mpq_init(c[i]);
300 zz2value(coeff[i], tmp);
301 mpq_set_z(c[i], tmp);
302
303 for (int j = 1; j <= i; ++j) {
304 zz2value(d.coeff[j], tmp);
305 mpq_set_z(qtmp, tmp);
306 mpq_mul(qtmp, qtmp, c[i-j]);
307 mpq_sub(c[i], c[i], qtmp);
308 }
309
310 zz2value(d.coeff[0], tmp);
311 mpq_set_z(qtmp, tmp);
312 mpq_div(c[i], c[i], qtmp);
313 }
314 if (sign == -1)
315 mpq_sub(count, count, c[len-1]);
316 else
317 mpq_add(count, count, c[len-1]);
318
319 value_clear(tmp);
320 mpq_clear(qtmp);
321 for (int i = 0; i < len; ++i)
322 mpq_clear(c[i]);
323 delete [] c;
324 }
325 };
326
327 class dpoly_n {
328 public:
329 Matrix *coeff;
330 ~dpoly_n() {
331 Matrix_Free(coeff);
332 }
333 dpoly_n(int d, ZZ& degree_0, ZZ& degree_1, int offset = 0) {
334 Value d0, d1;
335 value_init(d0);
336 value_init(d1);
337 zz2value(degree_0, d0);
338 zz2value(degree_1, d1);
339 coeff = Matrix_Alloc(d+1, d+1+1);
340 value_set_si(coeff->p[0][0], 1);
341 value_set_si(coeff->p[0][d+1], 1);
342 for (int i = 1; i <= d; ++i) {
343 value_multiply(coeff->p[i][0], coeff->p[i-1][0], d0);
344 Vector_Combine(coeff->p[i-1], coeff->p[i-1]+1, coeff->p[i]+1,
345 d1, d0, i);
346 value_set_si(coeff->p[i][d+1], i);
347 value_multiply(coeff->p[i][d+1], coeff->p[i][d+1], coeff->p[i-1][d+1]);
348 value_decrement(d0, d0);
349 }
350 value_clear(d0);
351 value_clear(d1);
352 }
353 void div(dpoly& d, Vector *count, ZZ& sign) {
354 int len = coeff->NbRows;
355 Matrix * c = Matrix_Alloc(coeff->NbRows, coeff->NbColumns);
356 Value tmp;
357 value_init(tmp);
358 for (int i = 0; i < len; ++i) {
359 Vector_Copy(coeff->p[i], c->p[i], len+1);
360 for (int j = 1; j <= i; ++j) {
361 zz2value(d.coeff[j], tmp);
362 value_multiply(tmp, tmp, c->p[i][len]);
363 value_oppose(tmp, tmp);
364 Vector_Combine(c->p[i], c->p[i-j], c->p[i],
365 c->p[i-j][len], tmp, len);
366 value_multiply(c->p[i][len], c->p[i][len], c->p[i-j][len]);
367 }
368 zz2value(d.coeff[0], tmp);
369 value_multiply(c->p[i][len], c->p[i][len], tmp);
370 }
371 if (sign == -1) {
372 value_set_si(tmp, -1);
373 Vector_Scale(c->p[len-1], count->p, tmp, len);
374 value_assign(count->p[len], c->p[len-1][len]);
375 } else
376 Vector_Copy(c->p[len-1], count->p, len+1);
377 Vector_Normalize(count->p, len+1);
378 value_clear(tmp);
379 Matrix_Free(c);
380 }
381 };
382
383 /*
384 * Barvinok's Decomposition of a simplicial cone
385 *
386 * Returns two lists of polyhedra
387 */
388 void barvinok_decompose(Polyhedron *C, Polyhedron **ppos, Polyhedron **pneg)
389 {
390 Polyhedron *pos = *ppos, *neg = *pneg;
391 vector<cone *> nonuni;
392 cone * c = new cone(C);
393 ZZ det = c->det;
394 int s = sign(det);
395 assert(det != 0);
396 if (abs(det) > 1) {
397 nonuni.push_back(c);
398 } else {
399 Polyhedron *p = Polyhedron_Copy(c->Cone);
400 p->next = pos;
401 pos = p;
402 delete c;
403 }
404 vec_ZZ lambda;
405 while (!nonuni.empty()) {
406 c = nonuni.back();
407 nonuni.pop_back();
408 Vector* v = c->short_vector(lambda);
409 for (int i = 0; i < c->Rays->NbRows - 1; ++i) {
410 if (lambda[i] == 0)
411 continue;
412 Matrix* M = Matrix_Copy(c->Rays);
413 Vector_Copy(v->p, M->p[i], v->Size);
414 cone * pc = new cone(M);
415 assert (pc->det != 0);
416 if (abs(pc->det) > 1) {
417 assert(abs(pc->det) < abs(c->det));
418 nonuni.push_back(pc);
419 } else {
420 Polyhedron *p = pc->poly();
421 pc->Cone = 0;
422 if (sign(pc->det) == s) {
423 p->next = pos;
424 pos = p;
425 } else {
426 p->next = neg;
427 neg = p;
428 }
429 delete pc;
430 }
431 Matrix_Free(M);
432 }
433 Vector_Free(v);
434 delete c;
435 }
436 *ppos = pos;
437 *pneg = neg;
438 }
439
440 /*
441 * Returns a single list of npos "positive" cones followed by nneg
442 * "negative" cones.
443 * The input cone is freed
444 */
445 void decompose(Polyhedron *cone, Polyhedron **parts, int *npos, int *nneg, unsigned MaxRays)
446 {
447 Polyhedron_Polarize(cone);
448 if (cone->NbRays - 1 != cone->Dimension) {
449 Polyhedron *tmp = cone;
450 cone = triangularize_cone(cone, MaxRays);
451 Polyhedron_Free(tmp);
452 }
453 Polyhedron *polpos = NULL, *polneg = NULL;
454 *npos = 0; *nneg = 0;
455 for (Polyhedron *Polar = cone; Polar; Polar = Polar->next)
456 barvinok_decompose(Polar, &polpos, &polneg);
457
458 Polyhedron *last;
459 for (Polyhedron *i = polpos; i; i = i->next) {
460 Polyhedron_Polarize(i);
461 ++*npos;
462 last = i;
463 }
464 for (Polyhedron *i = polneg; i; i = i->next) {
465 Polyhedron_Polarize(i);
466 ++*nneg;
467 }
468 if (last) {
469 last->next = polneg;
470 *parts = polpos;
471 } else
472 *parts = polneg;
473 Domain_Free(cone);
474 }
475
476 const int MAX_TRY=10;
477 /*
478 * Searches for a vector that is not othogonal to any
479 * of the rays in rays.
480 */
481 static void nonorthog(mat_ZZ& rays, vec_ZZ& lambda)
482 {
483 int dim = rays.NumCols();
484 bool found = false;
485 lambda.SetLength(dim);
486 for (int i = 2; !found && i <= 50*dim; i+=4) {
487 for (int j = 0; j < MAX_TRY; ++j) {
488 for (int k = 0; k < dim; ++k) {
489 int r = random_int(i)+2;
490 int v = (2*(r%2)-1) * (r >> 1);
491 lambda[k] = v;
492 }
493 int k = 0;
494 for (; k < rays.NumRows(); ++k)
495 if (lambda * rays[k] == 0)
496 break;
497 if (k == rays.NumRows()) {
498 found = true;
499 break;
500 }
501 }
502 }
503 assert(found);
504 }
505
506 static void add_rays(mat_ZZ& rays, Polyhedron *i, int *r)
507 {
508 unsigned dim = i->Dimension;
509 for (int k = 0; k < i->NbRays; ++k) {
510 if (!value_zero_p(i->Ray[k][dim+1]))
511 continue;
512 values2zz(i->Ray[k]+1, rays[(*r)++], dim);
513 }
514 }
515
516 void lattice_point(Value* values, Polyhedron *i, vec_ZZ& vertex)
517 {
518 unsigned dim = i->Dimension;
519 if(!value_one_p(values[dim])) {
520 Matrix* Rays = rays(i);
521 Matrix *inv = Matrix_Alloc(Rays->NbRows, Rays->NbColumns);
522 int ok = Matrix_Inverse(Rays, inv);
523 assert(ok);
524 Matrix_Free(Rays);
525 Rays = rays(i);
526 Vector *lambda = Vector_Alloc(dim+1);
527 Vector_Matrix_Product(values, inv, lambda->p);
528 Matrix_Free(inv);
529 for (int j = 0; j < dim; ++j)
530 mpz_cdiv_q(lambda->p[j], lambda->p[j], lambda->p[dim]);
531 value_set_si(lambda->p[dim], 1);
532 Vector *A = Vector_Alloc(dim+1);
533 Vector_Matrix_Product(lambda->p, Rays, A->p);
534 Vector_Free(lambda);
535 Matrix_Free(Rays);
536 values2zz(A->p, vertex, dim);
537 Vector_Free(A);
538 } else
539 values2zz(values, vertex, dim);
540 }
541
542 static evalue *term(int param, ZZ& c, Value *den = NULL)
543 {
544 evalue *EP = new evalue();
545 value_init(EP->d);
546 value_set_si(EP->d,0);
547 EP->x.p = new_enode(polynomial, 2, param + 1);
548 evalue_set_si(&EP->x.p->arr[0], 0, 1);
549 value_init(EP->x.p->arr[1].x.n);
550 if (den == NULL)
551 value_set_si(EP->x.p->arr[1].d, 1);
552 else
553 value_assign(EP->x.p->arr[1].d, *den);
554 zz2value(c, EP->x.p->arr[1].x.n);
555 return EP;
556 }
557
558 static void vertex_period(
559 Polyhedron *i, vec_ZZ& lambda, Matrix *T,
560 Value lcm, int p, Vector *val,
561 evalue *E, evalue* ev,
562 ZZ& offset)
563 {
564 unsigned nparam = T->NbRows - 1;
565 unsigned dim = i->Dimension;
566 Value tmp;
567 ZZ nump;
568
569 if (p == nparam) {
570 vec_ZZ vertex;
571 ZZ num, l;
572 Vector * values = Vector_Alloc(dim + 1);
573 Vector_Matrix_Product(val->p, T, values->p);
574 value_assign(values->p[dim], lcm);
575 lattice_point(values->p, i, vertex);
576 num = vertex * lambda;
577 value2zz(lcm, l);
578 num *= l;
579 num += offset;
580 value_init(ev->x.n);
581 zz2value(num, ev->x.n);
582 value_assign(ev->d, lcm);
583 Vector_Free(values);
584 return;
585 }
586
587 value_init(tmp);
588 vec_ZZ vertex;
589 values2zz(T->p[p], vertex, dim);
590 nump = vertex * lambda;
591 if (First_Non_Zero(val->p, p) == -1) {
592 value_assign(tmp, lcm);
593 evalue *ET = term(p, nump, &tmp);
594 eadd(ET, E);
595 free_evalue_refs(ET);
596 delete ET;
597 }
598
599 value_assign(tmp, lcm);
600 if (First_Non_Zero(T->p[p], dim) != -1)
601 Vector_Gcd(T->p[p], dim, &tmp);
602 Gcd(tmp, lcm, &tmp);
603 if (value_lt(tmp, lcm)) {
604 ZZ count;
605
606 value_division(tmp, lcm, tmp);
607 value_set_si(ev->d, 0);
608 ev->x.p = new_enode(periodic, VALUE_TO_INT(tmp), p+1);
609 value2zz(tmp, count);
610 do {
611 value_decrement(tmp, tmp);
612 --count;
613 ZZ new_offset = offset - count * nump;
614 value_assign(val->p[p], tmp);
615 vertex_period(i, lambda, T, lcm, p+1, val, E,
616 &ev->x.p->arr[VALUE_TO_INT(tmp)], new_offset);
617 } while (value_pos_p(tmp));
618 } else
619 vertex_period(i, lambda, T, lcm, p+1, val, E, ev, offset);
620 value_clear(tmp);
621 }
622
623 static void mask_r(Matrix *f, int nr, Vector *lcm, int p, Vector *val, evalue *ev)
624 {
625 unsigned nparam = lcm->Size;
626
627 if (p == nparam) {
628 Vector * prod = Vector_Alloc(f->NbRows);
629 Matrix_Vector_Product(f, val->p, prod->p);
630 int isint = 1;
631 for (int i = 0; i < nr; ++i) {
632 value_modulus(prod->p[i], prod->p[i], f->p[i][nparam+1]);
633 isint &= value_zero_p(prod->p[i]);
634 }
635 value_set_si(ev->d, 1);
636 value_init(ev->x.n);
637 value_set_si(ev->x.n, isint);
638 Vector_Free(prod);
639 return;
640 }
641
642 Value tmp;
643 value_init(tmp);
644 if (value_one_p(lcm->p[p]))
645 mask_r(f, nr, lcm, p+1, val, ev);
646 else {
647 value_assign(tmp, lcm->p[p]);
648 value_set_si(ev->d, 0);
649 ev->x.p = new_enode(periodic, VALUE_TO_INT(tmp), p+1);
650 do {
651 value_decrement(tmp, tmp);
652 value_assign(val->p[p], tmp);
653 mask_r(f, nr, lcm, p+1, val, &ev->x.p->arr[VALUE_TO_INT(tmp)]);
654 } while (value_pos_p(tmp));
655 }
656 value_clear(tmp);
657 }
658
659 static evalue *multi_monom(vec_ZZ& p)
660 {
661 evalue *X = new evalue();
662 value_init(X->d);
663 value_init(X->x.n);
664 unsigned nparam = p.length()-1;
665 zz2value(p[nparam], X->x.n);
666 value_set_si(X->d, 1);
667 for (int i = 0; i < nparam; ++i) {
668 if (p[i] == 0)
669 continue;
670 evalue *T = term(i, p[i]);
671 eadd(T, X);
672 free_evalue_refs(T);
673 delete T;
674 }
675 return X;
676 }
677
678 /*
679 * Check whether mapping polyhedron P on the affine combination
680 * num yields a range that has a fixed quotient on integer
681 * division by d
682 * If zero is true, then we are only interested in the quotient
683 * for the cases where the remainder is zero.
684 * Returns NULL if false and a newly allocated value if true.
685 */
686 static Value *fixed_quotient(Polyhedron *P, vec_ZZ& num, Value d, bool zero)
687 {
688 Value* ret = NULL;
689 int len = num.length();
690 Matrix *T = Matrix_Alloc(2, len);
691 zz2values(num, T->p[0]);
692 value_set_si(T->p[1][len-1], 1);
693 Polyhedron *I = Polyhedron_Image(P, T, P->NbConstraints);
694 Matrix_Free(T);
695
696 int i;
697 for (i = 0; i < I->NbRays; ++i)
698 if (value_zero_p(I->Ray[i][2])) {
699 Polyhedron_Free(I);
700 return NULL;
701 }
702
703 Value min, max;
704 value_init(min);
705 value_init(max);
706 int bounded = line_minmax(I, &min, &max);
707 assert(bounded);
708
709 if (zero)
710 mpz_cdiv_q(min, min, d);
711 else
712 mpz_fdiv_q(min, min, d);
713 mpz_fdiv_q(max, max, d);
714
715 if (value_eq(min, max)) {
716 ALLOC(Value, ret);
717 value_init(*ret);
718 value_assign(*ret, min);
719 }
720 value_clear(min);
721 value_clear(max);
722 return ret;
723 }
724
725 /*
726 * Normalize linear expression coef modulo m
727 * Removes common factor and reduces coefficients
728 * Returns index of first non-zero coefficient or len
729 */
730 static int normal_mod(Value *coef, int len, Value *m)
731 {
732 Value gcd;
733 value_init(gcd);
734
735 Vector_Gcd(coef, len, &gcd);
736 Gcd(gcd, *m, &gcd);
737 Vector_AntiScale(coef, coef, gcd, len);
738
739 value_division(*m, *m, gcd);
740 value_clear(gcd);
741
742 if (value_one_p(*m))
743 return len;
744
745 int j;
746 for (j = 0; j < len; ++j)
747 mpz_fdiv_r(coef[j], coef[j], *m);
748 for (j = 0; j < len; ++j)
749 if (value_notzero_p(coef[j]))
750 break;
751
752 return j;
753 }
754
755 #ifdef USE_MODULO
756 static void mask(Matrix *f, evalue *factor)
757 {
758 int nr = f->NbRows, nc = f->NbColumns;
759 int n;
760 bool found = false;
761 for (n = 0; n < nr && value_notzero_p(f->p[n][nc-1]); ++n)
762 if (value_notone_p(f->p[n][nc-1]) &&
763 value_notmone_p(f->p[n][nc-1]))
764 found = true;
765 if (!found)
766 return;
767
768 evalue EP;
769 nr = n;
770
771 Value m;
772 value_init(m);
773
774 evalue EV;
775 value_init(EV.d);
776 value_init(EV.x.n);
777 value_set_si(EV.x.n, 1);
778
779 for (n = 0; n < nr; ++n) {
780 value_assign(m, f->p[n][nc-1]);
781 if (value_one_p(m) || value_mone_p(m))
782 continue;
783
784 int j = normal_mod(f->p[n], nc-1, &m);
785 if (j == nc-1) {
786 free_evalue_refs(factor);
787 value_init(factor->d);
788 evalue_set_si(factor, 0, 1);
789 break;
790 }
791 vec_ZZ row;
792 values2zz(f->p[n], row, nc-1);
793 ZZ g;
794 value2zz(m, g);
795 if (j < (nc-1)-1 && row[j] > g/2) {
796 for (int k = j; k < (nc-1); ++k)
797 if (row[k] != 0)
798 row[k] = g - row[k];
799 }
800
801 value_init(EP.d);
802 value_set_si(EP.d, 0);
803 EP.x.p = new_enode(relation, 2, 0);
804 value_clear(EP.x.p->arr[1].d);
805 EP.x.p->arr[1] = *factor;
806 evalue *ev = &EP.x.p->arr[0];
807 value_set_si(ev->d, 0);
808 ev->x.p = new_enode(fractional, 3, -1);
809 evalue_set_si(&ev->x.p->arr[1], 0, 1);
810 evalue_set_si(&ev->x.p->arr[2], 1, 1);
811 evalue *E = multi_monom(row);
812 value_assign(EV.d, m);
813 emul(&EV, E);
814 value_clear(ev->x.p->arr[0].d);
815 ev->x.p->arr[0] = *E;
816 delete E;
817 *factor = EP;
818 }
819
820 value_clear(m);
821 free_evalue_refs(&EV);
822 }
823 #else
824 /*
825 *
826 */
827 static void mask(Matrix *f, evalue *factor)
828 {
829 int nr = f->NbRows, nc = f->NbColumns;
830 int n;
831 bool found = false;
832 for (n = 0; n < nr && value_notzero_p(f->p[n][nc-1]); ++n)
833 if (value_notone_p(f->p[n][nc-1]) &&
834 value_notmone_p(f->p[n][nc-1]))
835 found = true;
836 if (!found)
837 return;
838
839 Value tmp;
840 value_init(tmp);
841 nr = n;
842 unsigned np = nc - 2;
843 Vector *lcm = Vector_Alloc(np);
844 Vector *val = Vector_Alloc(nc);
845 Vector_Set(val->p, 0, nc);
846 value_set_si(val->p[np], 1);
847 Vector_Set(lcm->p, 1, np);
848 for (n = 0; n < nr; ++n) {
849 if (value_one_p(f->p[n][nc-1]) ||
850 value_mone_p(f->p[n][nc-1]))
851 continue;
852 for (int j = 0; j < np; ++j)
853 if (value_notzero_p(f->p[n][j])) {
854 Gcd(f->p[n][j], f->p[n][nc-1], &tmp);
855 value_division(tmp, f->p[n][nc-1], tmp);
856 value_lcm(tmp, lcm->p[j], &lcm->p[j]);
857 }
858 }
859 evalue EP;
860 value_init(EP.d);
861 mask_r(f, nr, lcm, 0, val, &EP);
862 value_clear(tmp);
863 Vector_Free(val);
864 Vector_Free(lcm);
865 emul(&EP,factor);
866 free_evalue_refs(&EP);
867 }
868 #endif
869
870 struct term_info {
871 evalue *E;
872 ZZ constant;
873 ZZ coeff;
874 int pos;
875 };
876
877 static bool mod_needed(Polyhedron *PD, vec_ZZ& num, Value d, evalue *E)
878 {
879 Value *q = fixed_quotient(PD, num, d, false);
880
881 if (!q)
882 return true;
883
884 value_oppose(*q, *q);
885 evalue EV;
886 value_init(EV.d);
887 value_set_si(EV.d, 1);
888 value_init(EV.x.n);
889 value_multiply(EV.x.n, *q, d);
890 eadd(&EV, E);
891 free_evalue_refs(&EV);
892 value_clear(*q);
893 free(q);
894 return false;
895 }
896
897 static void ceil_mod(Value *coef, int len, Value d, ZZ& f, evalue *EP, Polyhedron *PD)
898 {
899 Value m;
900 value_init(m);
901 value_set_si(m, -1);
902
903 Vector_Scale(coef, coef, m, len);
904
905 value_assign(m, d);
906 int j = normal_mod(coef, len, &m);
907
908 if (j == len) {
909 value_clear(m);
910 return;
911 }
912
913 vec_ZZ num;
914 values2zz(coef, num, len);
915
916 ZZ g;
917 value2zz(m, g);
918
919 evalue tmp;
920 value_init(tmp.d);
921 evalue_set_si(&tmp, 0, 1);
922
923 int p = j;
924 if (g % 2 == 0)
925 while (j < len-1 && (num[j] == g/2 || num[j] == 0))
926 ++j;
927 if ((j < len-1 && num[j] > g/2) || (j == len-1 && num[j] >= (g+1)/2)) {
928 for (int k = j; k < len-1; ++k)
929 if (num[k] != 0)
930 num[k] = g - num[k];
931 num[len-1] = g - 1 - num[len-1];
932 value_assign(tmp.d, m);
933 ZZ t = f*(g-1);
934 zz2value(t, tmp.x.n);
935 eadd(&tmp, EP);
936 f = -f;
937 }
938
939 if (p >= len-1) {
940 ZZ t = num[len-1] * f;
941 zz2value(t, tmp.x.n);
942 value_assign(tmp.d, m);
943 eadd(&tmp, EP);
944 } else {
945 evalue *E = multi_monom(num);
946 evalue EV;
947 value_init(EV.d);
948
949 if (PD && !mod_needed(PD, num, m, E)) {
950 value_init(EV.x.n);
951 zz2value(f, EV.x.n);
952 value_assign(EV.d, m);
953 emul(&EV, E);
954 eadd(E, EP);
955 } else {
956 value_init(EV.x.n);
957 value_set_si(EV.x.n, 1);
958 value_assign(EV.d, m);
959 emul(&EV, E);
960 value_clear(EV.x.n);
961 value_set_si(EV.d, 0);
962 EV.x.p = new_enode(fractional, 3, -1);
963 evalue_copy(&EV.x.p->arr[0], E);
964 evalue_set_si(&EV.x.p->arr[1], 0, 1);
965 value_init(EV.x.p->arr[2].x.n);
966 zz2value(f, EV.x.p->arr[2].x.n);
967 value_set_si(EV.x.p->arr[2].d, 1);
968
969 eadd(&EV, EP);
970 }
971
972 free_evalue_refs(&EV);
973 free_evalue_refs(E);
974 delete E;
975 }
976
977 free_evalue_refs(&tmp);
978
979 out:
980 value_clear(m);
981 }
982
983 evalue* bv_ceil3(Value *coef, int len, Value d, Polyhedron *P)
984 {
985 Vector *val = Vector_Alloc(len);
986
987 Value t;
988 value_init(t);
989 value_set_si(t, -1);
990 Vector_Scale(coef, val->p, t, len);
991 value_absolute(t, d);
992
993 vec_ZZ num;
994 values2zz(val->p, num, len);
995 evalue *EP = multi_monom(num);
996
997 evalue tmp;
998 value_init(tmp.d);
999 value_init(tmp.x.n);
1000 value_set_si(tmp.x.n, 1);
1001 value_assign(tmp.d, t);
1002
1003 emul(&tmp, EP);
1004
1005 ZZ one;
1006 one = 1;
1007 ceil_mod(val->p, len, t, one, EP, P);
1008 value_clear(t);
1009
1010 /* copy EP to malloc'ed evalue */
1011 evalue *E;
1012 ALLOC(evalue, E);
1013 *E = *EP;
1014 delete EP;
1015
1016 free_evalue_refs(&tmp);
1017 Vector_Free(val);
1018
1019 return E;
1020 }
1021
1022 #ifdef USE_MODULO
1023 evalue* lattice_point(
1024 Polyhedron *i, vec_ZZ& lambda, Matrix *W, Value lcm, Polyhedron *PD)
1025 {
1026 unsigned nparam = W->NbColumns - 1;
1027
1028 Matrix* Rays = rays2(i);
1029 Matrix *T = Transpose(Rays);
1030 Matrix *T2 = Matrix_Copy(T);
1031 Matrix *inv = Matrix_Alloc(T2->NbRows, T2->NbColumns);
1032 int ok = Matrix_Inverse(T2, inv);
1033 assert(ok);
1034 Matrix_Free(Rays);
1035 Matrix_Free(T2);
1036 mat_ZZ vertex;
1037 matrix2zz(W, vertex, W->NbRows, W->NbColumns);
1038
1039 vec_ZZ num;
1040 num = lambda * vertex;
1041
1042 evalue *EP = multi_monom(num);
1043
1044 evalue tmp;
1045 value_init(tmp.d);
1046 value_init(tmp.x.n);
1047 value_set_si(tmp.x.n, 1);
1048 value_assign(tmp.d, lcm);
1049
1050 emul(&tmp, EP);
1051
1052 Matrix *L = Matrix_Alloc(inv->NbRows, W->NbColumns);
1053 Matrix_Product(inv, W, L);
1054
1055 mat_ZZ RT;
1056 matrix2zz(T, RT, T->NbRows, T->NbColumns);
1057 Matrix_Free(T);
1058
1059 vec_ZZ p = lambda * RT;
1060
1061 for (int i = 0; i < L->NbRows; ++i) {
1062 ceil_mod(L->p[i], nparam+1, lcm, p[i], EP, PD);
1063 }
1064
1065 Matrix_Free(L);
1066
1067 Matrix_Free(inv);
1068 free_evalue_refs(&tmp);
1069 return EP;
1070 }
1071 #else
1072 evalue* lattice_point(
1073 Polyhedron *i, vec_ZZ& lambda, Matrix *W, Value lcm, Polyhedron *PD)
1074 {
1075 Matrix *T = Transpose(W);
1076 unsigned nparam = T->NbRows - 1;
1077
1078 evalue *EP = new evalue();
1079 value_init(EP->d);
1080 evalue_set_si(EP, 0, 1);
1081
1082 evalue ev;
1083 Vector *val = Vector_Alloc(nparam+1);
1084 value_set_si(val->p[nparam], 1);
1085 ZZ offset(INIT_VAL, 0);
1086 value_init(ev.d);
1087 vertex_period(i, lambda, T, lcm, 0, val, EP, &ev, offset);
1088 Vector_Free(val);
1089 eadd(&ev, EP);
1090 free_evalue_refs(&ev);
1091
1092 Matrix_Free(T);
1093
1094 reduce_evalue(EP);
1095
1096 return EP;
1097 }
1098 #endif
1099
1100 void lattice_point(
1101 Param_Vertices* V, Polyhedron *i, vec_ZZ& lambda, term_info* term,
1102 Polyhedron *PD)
1103 {
1104 unsigned nparam = V->Vertex->NbColumns - 2;
1105 unsigned dim = i->Dimension;
1106 mat_ZZ vertex;
1107 vertex.SetDims(V->Vertex->NbRows, nparam+1);
1108 Value lcm, tmp;
1109 value_init(lcm);
1110 value_init(tmp);
1111 value_set_si(lcm, 1);
1112 for (int j = 0; j < V->Vertex->NbRows; ++j) {
1113 value_lcm(lcm, V->Vertex->p[j][nparam+1], &lcm);
1114 }
1115 if (value_notone_p(lcm)) {
1116 Matrix * mv = Matrix_Alloc(dim, nparam+1);
1117 for (int j = 0 ; j < dim; ++j) {
1118 value_division(tmp, lcm, V->Vertex->p[j][nparam+1]);
1119 Vector_Scale(V->Vertex->p[j], mv->p[j], tmp, nparam+1);
1120 }
1121
1122 term->E = lattice_point(i, lambda, mv, lcm, PD);
1123 term->constant = 0;
1124
1125 Matrix_Free(mv);
1126 value_clear(lcm);
1127 value_clear(tmp);
1128 return;
1129 }
1130 for (int i = 0; i < V->Vertex->NbRows; ++i) {
1131 assert(value_one_p(V->Vertex->p[i][nparam+1])); // for now
1132 values2zz(V->Vertex->p[i], vertex[i], nparam+1);
1133 }
1134
1135 vec_ZZ num;
1136 num = lambda * vertex;
1137
1138 int p = -1;
1139 int nn = 0;
1140 for (int j = 0; j < nparam; ++j)
1141 if (num[j] != 0) {
1142 ++nn;
1143 p = j;
1144 }
1145 if (nn >= 2) {
1146 term->E = multi_monom(num);
1147 term->constant = 0;
1148 } else {
1149 term->E = NULL;
1150 term->constant = num[nparam];
1151 term->pos = p;
1152 if (p != -1)
1153 term->coeff = num[p];
1154 }
1155
1156 value_clear(lcm);
1157 value_clear(tmp);
1158 }
1159
1160 void normalize(Polyhedron *i, vec_ZZ& lambda, ZZ& sign, ZZ& num, vec_ZZ& den)
1161 {
1162 unsigned dim = i->Dimension;
1163
1164 int r = 0;
1165 mat_ZZ rays;
1166 rays.SetDims(dim, dim);
1167 add_rays(rays, i, &r);
1168 den = rays * lambda;
1169 int change = 0;
1170
1171 for (int j = 0; j < den.length(); ++j) {
1172 if (den[j] > 0)
1173 change ^= 1;
1174 else {
1175 den[j] = abs(den[j]);
1176 num += den[j];
1177 }
1178 }
1179 if (change)
1180 sign = -sign;
1181 }
1182
1183 typedef Polyhedron * Polyhedron_p;
1184
1185 void barvinok_count(Polyhedron *P, Value* result, unsigned NbMaxCons)
1186 {
1187 Polyhedron ** vcone;
1188 vec_ZZ sign;
1189 int ncone = 0;
1190 sign.SetLength(ncone);
1191 unsigned dim;
1192 int allocated = 0;
1193 Value factor;
1194 Polyhedron *Q;
1195 int r = 0;
1196
1197 if (emptyQ(P)) {
1198 value_set_si(*result, 0);
1199 return;
1200 }
1201 if (P->NbBid == 0)
1202 for (; r < P->NbRays; ++r)
1203 if (value_zero_p(P->Ray[r][P->Dimension+1]))
1204 break;
1205 if (P->NbBid !=0 || r < P->NbRays) {
1206 value_set_si(*result, -1);
1207 return;
1208 }
1209 if (P->NbEq != 0) {
1210 P = remove_equalities(P);
1211 if (emptyQ(P)) {
1212 Polyhedron_Free(P);
1213 value_set_si(*result, 0);
1214 return;
1215 }
1216 allocated = 1;
1217 }
1218 value_init(factor);
1219 value_set_si(factor, 1);
1220 Q = Polyhedron_Reduce(P, &factor);
1221 if (Q) {
1222 if (allocated)
1223 Polyhedron_Free(P);
1224 P = Q;
1225 allocated = 1;
1226 }
1227 if (P->Dimension == 0) {
1228 value_assign(*result, factor);
1229 if (allocated)
1230 Polyhedron_Free(P);
1231 value_clear(factor);
1232 return;
1233 }
1234
1235 dim = P->Dimension;
1236 vcone = new Polyhedron_p[P->NbRays];
1237
1238 for (int j = 0; j < P->NbRays; ++j) {
1239 int npos, nneg;
1240 Polyhedron *C = supporting_cone(P, j);
1241 decompose(C, &vcone[j], &npos, &nneg, NbMaxCons);
1242 ncone += npos + nneg;
1243 sign.SetLength(ncone);
1244 for (int k = 0; k < npos; ++k)
1245 sign[ncone-nneg-k-1] = 1;
1246 for (int k = 0; k < nneg; ++k)
1247 sign[ncone-k-1] = -1;
1248 }
1249
1250 mat_ZZ rays;
1251 rays.SetDims(ncone * dim, dim);
1252 r = 0;
1253 for (int j = 0; j < P->NbRays; ++j) {
1254 for (Polyhedron *i = vcone[j]; i; i = i->next) {
1255 assert(i->NbRays-1 == dim);
1256 add_rays(rays, i, &r);
1257 }
1258 }
1259 vec_ZZ lambda;
1260 nonorthog(rays, lambda);
1261
1262 ZZ num;
1263 vec_ZZ den;
1264 den.SetLength(dim);
1265
1266 int f = 0;
1267 vec_ZZ vertex;
1268 mpq_t count;
1269 mpq_init(count);
1270 for (int j = 0; j < P->NbRays; ++j) {
1271 for (Polyhedron *i = vcone[j]; i; i = i->next) {
1272 lattice_point(P->Ray[j]+1, i, vertex);
1273 num = vertex * lambda;
1274 normalize(i, lambda, sign[f], num, den);
1275
1276 dpoly d(dim, num);
1277 dpoly n(dim, den[0], 1);
1278 for (int k = 1; k < dim; ++k) {
1279 dpoly fact(dim, den[k], 1);
1280 n *= fact;
1281 }
1282 d.div(n, count, sign[f]);
1283
1284 ++f;
1285 }
1286 Domain_Free(vcone[j]);
1287 }
1288
1289 assert(value_one_p(&count[0]._mp_den));
1290 value_multiply(*result, &count[0]._mp_num, factor);
1291 mpq_clear(count);
1292
1293 delete [] vcone;
1294
1295 if (allocated)
1296 Polyhedron_Free(P);
1297 value_clear(factor);
1298 }
1299
1300 static void uni_polynom(int param, Vector *c, evalue *EP)
1301 {
1302 unsigned dim = c->Size-2;
1303 value_init(EP->d);
1304 value_set_si(EP->d,0);
1305 EP->x.p = new_enode(polynomial, dim+1, param+1);
1306 for (int j = 0; j <= dim; ++j)
1307 evalue_set(&EP->x.p->arr[j], c->p[j], c->p[dim+1]);
1308 }
1309
1310 static void multi_polynom(Vector *c, evalue* X, evalue *EP)
1311 {
1312 unsigned dim = c->Size-2;
1313 evalue EC;
1314
1315 value_init(EC.d);
1316 evalue_set(&EC, c->p[dim], c->p[dim+1]);
1317
1318 value_init(EP->d);
1319 evalue_set(EP, c->p[dim], c->p[dim+1]);
1320
1321 for (int i = dim-1; i >= 0; --i) {
1322 emul(X, EP);
1323 value_assign(EC.x.n, c->p[i]);
1324 eadd(&EC, EP);
1325 }
1326 free_evalue_refs(&EC);
1327 }
1328
1329 Polyhedron *unfringe (Polyhedron *P, unsigned MaxRays)
1330 {
1331 int len = P->Dimension+2;
1332 Polyhedron *T, *R = P;
1333 Value g;
1334 value_init(g);
1335 Vector *row = Vector_Alloc(len);
1336 value_set_si(row->p[0], 1);
1337
1338 R = DomainConstraintSimplify(Polyhedron_Copy(P), MaxRays);
1339
1340 Matrix *M = Matrix_Alloc(2, len-1);
1341 value_set_si(M->p[1][len-2], 1);
1342 for (int v = 0; v < P->Dimension; ++v) {
1343 value_set_si(M->p[0][v], 1);
1344 Polyhedron *I = Polyhedron_Image(P, M, 2+1);
1345 value_set_si(M->p[0][v], 0);
1346 for (int r = 0; r < I->NbConstraints; ++r) {
1347 if (value_zero_p(I->Constraint[r][0]))
1348 continue;
1349 if (value_zero_p(I->Constraint[r][1]))
1350 continue;
1351 if (value_one_p(I->Constraint[r][1]))
1352 continue;
1353 if (value_mone_p(I->Constraint[r][1]))
1354 continue;
1355 value_absolute(g, I->Constraint[r][1]);
1356 Vector_Set(row->p+1, 0, len-2);
1357 value_division(row->p[1+v], I->Constraint[r][1], g);
1358 mpz_fdiv_q(row->p[len-1], I->Constraint[r][2], g);
1359 T = R;
1360 R = AddConstraints(row->p, 1, R, MaxRays);
1361 if (T != P)
1362 Polyhedron_Free(T);
1363 }
1364 Polyhedron_Free(I);
1365 }
1366 Matrix_Free(M);
1367 Vector_Free(row);
1368 value_clear(g);
1369 return R;
1370 }
1371
1372 static Polyhedron *reduce_domain(Polyhedron *D, Matrix *CT, Polyhedron *CEq,
1373 Polyhedron **fVD, int nd, unsigned MaxRays)
1374 {
1375 assert(CEq);
1376
1377 Polyhedron *Dt;
1378 Dt = CT ? DomainPreimage(D, CT, MaxRays) : D;
1379 Polyhedron *rVD = DomainIntersection(Dt, CEq, MaxRays);
1380
1381 /* if rVD is empty or too small in geometric dimension */
1382 if(!rVD || emptyQ(rVD) ||
1383 (rVD->Dimension-rVD->NbEq < Dt->Dimension-Dt->NbEq-CEq->NbEq)) {
1384 if(rVD)
1385 Domain_Free(rVD);
1386 if (CT)
1387 Domain_Free(Dt);
1388 return 0; /* empty validity domain */
1389 }
1390
1391 if (CT)
1392 Domain_Free(Dt);
1393
1394 fVD[nd] = Domain_Copy(rVD);
1395 for (int i = 0 ; i < nd; ++i) {
1396 Polyhedron *I = DomainIntersection(fVD[nd], fVD[i], MaxRays);
1397 if (emptyQ(I)) {
1398 Domain_Free(I);
1399 continue;
1400 }
1401 Polyhedron *F = DomainSimplify(I, fVD[nd], MaxRays);
1402 if (F->NbEq == 1) {
1403 Polyhedron *T = rVD;
1404 rVD = DomainDifference(rVD, F, MaxRays);
1405 Domain_Free(T);
1406 }
1407 Domain_Free(F);
1408 Domain_Free(I);
1409 }
1410
1411 rVD = DomainConstraintSimplify(rVD, MaxRays);
1412 if (emptyQ(rVD)) {
1413 Domain_Free(fVD[nd]);
1414 Domain_Free(rVD);
1415 return 0;
1416 }
1417
1418 Value c;
1419 value_init(c);
1420 barvinok_count(rVD, &c, MaxRays);
1421 if (value_zero_p(c)) {
1422 Domain_Free(rVD);
1423 rVD = 0;
1424 }
1425 value_clear(c);
1426
1427 return rVD;
1428 }
1429
1430 static bool Polyhedron_is_infinite(Polyhedron *P, unsigned nparam)
1431 {
1432 int r;
1433 for (r = 0; r < P->NbRays; ++r)
1434 if (value_zero_p(P->Ray[r][0]) ||
1435 value_zero_p(P->Ray[r][P->Dimension+1])) {
1436 int i;
1437 for (i = P->Dimension - nparam; i < P->Dimension; ++i)
1438 if (value_notzero_p(P->Ray[r][i+1]))
1439 break;
1440 if (i >= P->Dimension)
1441 break;
1442 }
1443 return r < P->NbRays;
1444 }
1445
1446 evalue* barvinok_enumerate_ev(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
1447 {
1448 //P = unfringe(P, MaxRays);
1449 Polyhedron *CEq = NULL, *rVD, *pVD, *CA;
1450 Matrix *CT = NULL;
1451 Param_Polyhedron *PP = NULL;
1452 Param_Domain *D, *next;
1453 Param_Vertices *V;
1454 int r = 0;
1455 unsigned nparam = C->Dimension;
1456 evalue *eres;
1457 ALLOC(evalue, eres);
1458 value_init(eres->d);
1459 value_set_si(eres->d, 0);
1460
1461 evalue factor;
1462 value_init(factor.d);
1463 evalue_set_si(&factor, 1, 1);
1464
1465 CA = align_context(C, P->Dimension, MaxRays);
1466 P = DomainIntersection(P, CA, MaxRays);
1467 Polyhedron_Free(CA);
1468
1469 if (C->Dimension == 0 || emptyQ(P)) {
1470 constant:
1471 eres->x.p = new_enode(partition, 2, C->Dimension);
1472 EVALUE_SET_DOMAIN(eres->x.p->arr[0],
1473 DomainConstraintSimplify(CEq ? CEq : Polyhedron_Copy(C), MaxRays));
1474 value_set_si(eres->x.p->arr[1].d, 1);
1475 value_init(eres->x.p->arr[1].x.n);
1476 if (emptyQ(P))
1477 value_set_si(eres->x.p->arr[1].x.n, 0);
1478 else
1479 barvinok_count(P, &eres->x.p->arr[1].x.n, MaxRays);
1480 out:
1481 emul(&factor, eres);
1482 reduce_evalue(eres);
1483 free_evalue_refs(&factor);
1484 Polyhedron_Free(P);
1485 if (CT)
1486 Matrix_Free(CT);
1487 if (PP)
1488 Param_Polyhedron_Free(PP);
1489
1490 return eres;
1491 }
1492 if (Polyhedron_is_infinite(P, nparam))
1493 goto constant;
1494
1495 if (P->NbEq != 0) {
1496 Matrix *f;
1497 P = remove_equalities_p(P, P->Dimension-nparam, &f);
1498 mask(f, &factor);
1499 Matrix_Free(f);
1500 }
1501 if (P->Dimension == nparam) {
1502 CEq = P;
1503 P = Universe_Polyhedron(0);
1504 goto constant;
1505 }
1506
1507 Polyhedron *Q = ParamPolyhedron_Reduce(P, P->Dimension-nparam, &factor);
1508 if (Q) {
1509 Polyhedron_Free(P);
1510 if (Q->Dimension == nparam) {
1511 CEq = Q;
1512 P = Universe_Polyhedron(0);
1513 goto constant;
1514 }
1515 P = Q;
1516 }
1517 Polyhedron *oldP = P;
1518 PP = Polyhedron2Param_SimplifiedDomain(&P,C,MaxRays,&CEq,&CT);
1519 if (P != oldP)
1520 Polyhedron_Free(oldP);
1521
1522 if (isIdentity(CT)) {
1523 Matrix_Free(CT);
1524 CT = NULL;
1525 } else {
1526 assert(CT->NbRows != CT->NbColumns);
1527 if (CT->NbRows == 1) // no more parameters
1528 goto constant;
1529 nparam = CT->NbRows - 1;
1530 }
1531
1532 unsigned dim = P->Dimension - nparam;
1533 Polyhedron ** vcone = new Polyhedron_p[PP->nbV];
1534 int * npos = new int[PP->nbV];
1535 int * nneg = new int[PP->nbV];
1536 ZZ sign;
1537
1538 int i;
1539 for (i = 0, V = PP->V; V; ++i, V = V->next) {
1540 Polyhedron *C = supporting_cone_p(P, V);
1541 decompose(C, &vcone[i], &npos[i], &nneg[i], MaxRays);
1542 }
1543
1544 Vector *c = Vector_Alloc(dim+2);
1545
1546 int nd;
1547 for (nd = 0, D=PP->D; D; ++nd, D=D->next);
1548 struct section { Polyhedron *D; evalue E; };
1549 section *s = new section[nd];
1550 Polyhedron **fVD = new Polyhedron_p[nd];
1551
1552 for(nd = 0, D=PP->D; D; D=next) {
1553 next = D->next;
1554
1555 Polyhedron *rVD = reduce_domain(D->Domain, CT, CEq,
1556 fVD, nd, MaxRays);
1557 if (!rVD)
1558 continue;
1559
1560 pVD = CT ? DomainImage(rVD,CT,MaxRays) : rVD;
1561
1562 int ncone = 0;
1563 FORALL_PVertex_in_ParamPolyhedron(V,D,PP) // _i is internal counter
1564 ncone += npos[_i] + nneg[_i];
1565 END_FORALL_PVertex_in_ParamPolyhedron;
1566
1567 mat_ZZ rays;
1568 rays.SetDims(ncone * dim, dim);
1569 r = 0;
1570 FORALL_PVertex_in_ParamPolyhedron(V,D,PP) // _i is internal counter
1571 for (Polyhedron *i = vcone[_i]; i; i = i->next) {
1572 assert(i->NbRays-1 == dim);
1573 add_rays(rays, i, &r);
1574 }
1575 END_FORALL_PVertex_in_ParamPolyhedron;
1576 vec_ZZ lambda;
1577 nonorthog(rays, lambda);
1578
1579 vec_ZZ den;
1580 den.SetLength(dim);
1581 term_info num;
1582
1583 value_init(s[nd].E.d);
1584 evalue_set_si(&s[nd].E, 0, 1);
1585 mpq_t count;
1586 mpq_init(count);
1587 FORALL_PVertex_in_ParamPolyhedron(V,D,PP)
1588 int f = 0;
1589 for (Polyhedron *i = vcone[_i]; i; i = i->next) {
1590 sign = f < npos[_i] ? 1 : -1;
1591 lattice_point(V, i, lambda, &num, pVD);
1592 normalize(i, lambda, sign, num.constant, den);
1593
1594 dpoly n(dim, den[0], 1);
1595 for (int k = 1; k < dim; ++k) {
1596 dpoly fact(dim, den[k], 1);
1597 n *= fact;
1598 }
1599 if (num.E != NULL) {
1600 ZZ one(INIT_VAL, 1);
1601 dpoly_n d(dim, num.constant, one);
1602 d.div(n, c, sign);
1603 evalue EV;
1604 multi_polynom(c, num.E, &EV);
1605 eadd(&EV , &s[nd].E);
1606 free_evalue_refs(&EV);
1607 free_evalue_refs(num.E);
1608 delete num.E;
1609 } else if (num.pos != -1) {
1610 dpoly_n d(dim, num.constant, num.coeff);
1611 d.div(n, c, sign);
1612 evalue EV;
1613 uni_polynom(num.pos, c, &EV);
1614 eadd(&EV , &s[nd].E);
1615 free_evalue_refs(&EV);
1616 } else {
1617 mpq_set_si(count, 0, 1);
1618 dpoly d(dim, num.constant);
1619 d.div(n, count, sign);
1620 evalue EV;
1621 value_init(EV.d);
1622 evalue_set(&EV, &count[0]._mp_num, &count[0]._mp_den);
1623 eadd(&EV , &s[nd].E);
1624 free_evalue_refs(&EV);
1625 }
1626 ++f;
1627 }
1628 END_FORALL_PVertex_in_ParamPolyhedron;
1629
1630 mpq_clear(count);
1631
1632 if (CT)
1633 addeliminatedparams_evalue(&s[nd].E, CT);
1634 s[nd].D = rVD;
1635 ++nd;
1636 if (rVD != pVD)
1637 Domain_Free(pVD);
1638 }
1639
1640 if (nd == 0)
1641 evalue_set_si(eres, 0, 1);
1642 else {
1643 eres->x.p = new_enode(partition, 2*nd, C->Dimension);
1644 for (int j = 0; j < nd; ++j) {
1645 EVALUE_SET_DOMAIN(eres->x.p->arr[2*j], s[j].D);
1646 value_clear(eres->x.p->arr[2*j+1].d);
1647 eres->x.p->arr[2*j+1] = s[j].E;
1648 Domain_Free(fVD[j]);
1649 }
1650 }
1651 delete [] s;
1652 delete [] fVD;
1653
1654 Vector_Free(c);
1655
1656 for (int j = 0; j < PP->nbV; ++j)
1657 Domain_Free(vcone[j]);
1658 delete [] vcone;
1659 delete [] npos;
1660 delete [] nneg;
1661
1662 if (CEq)
1663 Polyhedron_Free(CEq);
1664
1665 goto out;
1666 }
1667
1668 Enumeration* barvinok_enumerate(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
1669 {
1670 evalue *EP = barvinok_enumerate_ev(P, C, MaxRays);
1671
1672 return partition2enumeration(EP);
1673 }
1674
1675 static void SwapColumns(Value **V, int n, int i, int j)
1676 {
1677 for (int r = 0; r < n; ++r)
1678 value_swap(V[r][i], V[r][j]);
1679 }
1680
1681 static void SwapColumns(Polyhedron *P, int i, int j)
1682 {
1683 SwapColumns(P->Constraint, P->NbConstraints, i, j);
1684 SwapColumns(P->Ray, P->NbRays, i, j);
1685 }
1686
1687 static void negative_test_constraint(Value *l, Value *u, Value *c, int pos,
1688 int len, Value *v)
1689 {
1690 value_oppose(*v, u[pos+1]);
1691 Vector_Combine(l+1, u+1, c+1, *v, l[pos+1], len-1);
1692 value_multiply(*v, *v, l[pos+1]);
1693 value_substract(c[len-1], c[len-1], *v);
1694 value_set_si(*v, -1);
1695 Vector_Scale(c+1, c+1, *v, len-1);
1696 value_decrement(c[len-1], c[len-1]);
1697 ConstraintSimplify(c, c, len, v);
1698 }
1699
1700 static void oppose_constraint(Value *c, int len, Value *v)
1701 {
1702 value_set_si(*v, -1);
1703 Vector_Scale(c+1, c+1, *v, len-1);
1704 value_decrement(c[len-1], c[len-1]);
1705 }
1706
1707 static bool SplitOnConstraint(Polyhedron *P, int i, int l, int u,
1708 int nvar, int len, int exist, int MaxRays,
1709 Vector *row, Value& f, bool independent,
1710 Polyhedron **pos, Polyhedron **neg)
1711 {
1712 negative_test_constraint(P->Constraint[l], P->Constraint[u],
1713 row->p, nvar+i, len, &f);
1714 *neg = AddConstraints(row->p, 1, P, MaxRays);
1715
1716 /* We found an independent, but useless constraint
1717 * Maybe we should detect this earlier and not
1718 * mark the variable as INDEPENDENT
1719 */
1720 if (emptyQ((*neg))) {
1721 Polyhedron_Free(*neg);
1722 return false;
1723 }
1724
1725 oppose_constraint(row->p, len, &f);
1726 *pos = AddConstraints(row->p, 1, P, MaxRays);
1727
1728 if (emptyQ((*pos))) {
1729 Polyhedron_Free(*neg);
1730 Polyhedron_Free(*pos);
1731 return false;
1732 }
1733
1734 return true;
1735 }
1736
1737 /*
1738 * unimodularly transform P such that constraint r is transformed
1739 * into a constraint that involves only a single (the first)
1740 * existential variable
1741 *
1742 */
1743 static Polyhedron *rotate_along(Polyhedron *P, int r, int nvar, int exist,
1744 unsigned MaxRays)
1745 {
1746 Value g;
1747 value_init(g);
1748
1749 Vector *row = Vector_Alloc(exist);
1750 Vector_Copy(P->Constraint[r]+1+nvar, row->p, exist);
1751 Vector_Gcd(row->p, exist, &g);
1752 if (value_notone_p(g))
1753 Vector_AntiScale(row->p, row->p, g, exist);
1754 value_clear(g);
1755
1756 Matrix *M = unimodular_complete(row);
1757 Matrix *M2 = Matrix_Alloc(P->Dimension+1, P->Dimension+1);
1758 for (r = 0; r < nvar; ++r)
1759 value_set_si(M2->p[r][r], 1);
1760 for ( ; r < nvar+exist; ++r)
1761 Vector_Copy(M->p[r-nvar], M2->p[r]+nvar, exist);
1762 for ( ; r < P->Dimension+1; ++r)
1763 value_set_si(M2->p[r][r], 1);
1764 Polyhedron *T = Polyhedron_Image(P, M2, MaxRays);
1765
1766 Matrix_Free(M2);
1767 Matrix_Free(M);
1768 Vector_Free(row);
1769
1770 return T;
1771 }
1772
1773 static bool SplitOnVar(Polyhedron *P, int i,
1774 int nvar, int len, int exist, int MaxRays,
1775 Vector *row, Value& f, bool independent,
1776 Polyhedron **pos, Polyhedron **neg)
1777 {
1778 int j;
1779
1780 for (int l = P->NbEq; l < P->NbConstraints; ++l) {
1781 if (value_negz_p(P->Constraint[l][nvar+i+1]))
1782 continue;
1783
1784 if (independent) {
1785 for (j = 0; j < exist; ++j)
1786 if (j != i && value_notzero_p(P->Constraint[l][nvar+j+1]))
1787 break;
1788 if (j < exist)
1789 continue;
1790 }
1791
1792 for (int u = P->NbEq; u < P->NbConstraints; ++u) {
1793 if (value_posz_p(P->Constraint[u][nvar+i+1]))
1794 continue;
1795
1796 if (independent) {
1797 for (j = 0; j < exist; ++j)
1798 if (j != i && value_notzero_p(P->Constraint[u][nvar+j+1]))
1799 break;
1800 if (j < exist)
1801 continue;
1802 }
1803
1804 if (SplitOnConstraint(P, i, l, u,
1805 nvar, len, exist, MaxRays,
1806 row, f, independent,
1807 pos, neg)) {
1808 if (independent) {
1809 if (i != 0)
1810 SwapColumns(*neg, nvar+1, nvar+1+i);
1811 }
1812 return true;
1813 }
1814 }
1815 }
1816
1817 return false;
1818 }
1819
1820 static bool double_bound_pair(Polyhedron *P, int nvar, int exist,
1821 int i, int l1, int l2,
1822 Polyhedron **pos, Polyhedron **neg)
1823 {
1824 Value f;
1825 value_init(f);
1826 Vector *row = Vector_Alloc(P->Dimension+2);
1827 value_set_si(row->p[0], 1);
1828 value_oppose(f, P->Constraint[l1][nvar+i+1]);
1829 Vector_Combine(P->Constraint[l1]+1, P->Constraint[l2]+1,
1830 row->p+1,
1831 P->Constraint[l2][nvar+i+1], f,
1832 P->Dimension+1);
1833 ConstraintSimplify(row->p, row->p, P->Dimension+2, &f);
1834 *pos = AddConstraints(row->p, 1, P, 0);
1835 value_set_si(f, -1);
1836 Vector_Scale(row->p+1, row->p+1, f, P->Dimension+1);
1837 value_decrement(row->p[P->Dimension+1], row->p[P->Dimension+1]);
1838 *neg = AddConstraints(row->p, 1, P, 0);
1839 Vector_Free(row);
1840 value_clear(f);
1841
1842 return !emptyQ((*pos)) && !emptyQ((*neg));
1843 }
1844
1845 static bool double_bound(Polyhedron *P, int nvar, int exist,
1846 Polyhedron **pos, Polyhedron **neg)
1847 {
1848 for (int i = 0; i < exist; ++i) {
1849 int l1, l2;
1850 for (l1 = P->NbEq; l1 < P->NbConstraints; ++l1) {
1851 if (value_negz_p(P->Constraint[l1][nvar+i+1]))
1852 continue;
1853 for (l2 = l1 + 1; l2 < P->NbConstraints; ++l2) {
1854 if (value_negz_p(P->Constraint[l2][nvar+i+1]))
1855 continue;
1856 if (double_bound_pair(P, nvar, exist, i, l1, l2, pos, neg))
1857 return true;
1858 }
1859 }
1860 for (l1 = P->NbEq; l1 < P->NbConstraints; ++l1) {
1861 if (value_posz_p(P->Constraint[l1][nvar+i+1]))
1862 continue;
1863 if (l1 < P->NbConstraints)
1864 for (l2 = l1 + 1; l2 < P->NbConstraints; ++l2) {
1865 if (value_posz_p(P->Constraint[l2][nvar+i+1]))
1866 continue;
1867 if (double_bound_pair(P, nvar, exist, i, l1, l2, pos, neg))
1868 return true;
1869 }
1870 }
1871 return false;
1872 }
1873 return false;
1874 }
1875
1876 enum constraint {
1877 ALL_POS = 1 << 0,
1878 ONE_NEG = 1 << 1,
1879 INDEPENDENT = 1 << 2
1880 };
1881
1882 static evalue* enumerate_or(Polyhedron *D,
1883 unsigned exist, unsigned nparam, unsigned MaxRays)
1884 {
1885 #ifdef DEBUG_ER
1886 fprintf(stderr, "\nER: Or\n");
1887 #endif /* DEBUG_ER */
1888
1889 Polyhedron *N = D->next;
1890 D->next = 0;
1891 evalue *EP =
1892 barvinok_enumerate_e(D, exist, nparam, MaxRays);
1893 Polyhedron_Free(D);
1894
1895 for (D = N; D; D = N) {
1896 N = D->next;
1897 D->next = 0;
1898
1899 evalue *EN =
1900 barvinok_enumerate_e(D, exist, nparam, MaxRays);
1901
1902 eor(EN, EP);
1903 free_evalue_refs(EN);
1904 free(EN);
1905 Polyhedron_Free(D);
1906 }
1907
1908 reduce_evalue(EP);
1909
1910 return EP;
1911 }
1912
1913 static evalue* enumerate_sum(Polyhedron *P,
1914 unsigned exist, unsigned nparam, unsigned MaxRays)
1915 {
1916 int nvar = P->Dimension - exist - nparam;
1917 int toswap = nvar < exist ? nvar : exist;
1918 for (int i = 0; i < toswap; ++i)
1919 SwapColumns(P, 1 + i, nvar+exist - i);
1920 nparam += nvar;
1921
1922 #ifdef DEBUG_ER
1923 fprintf(stderr, "\nER: Sum\n");
1924 #endif /* DEBUG_ER */
1925
1926 evalue *EP = barvinok_enumerate_e(P, exist, nparam, MaxRays);
1927
1928 for (int i = 0; i < /* nvar */ nparam; ++i) {
1929 Matrix *C = Matrix_Alloc(1, 1 + nparam + 1);
1930 value_set_si(C->p[0][0], 1);
1931 evalue split;
1932 value_init(split.d);
1933 value_set_si(split.d, 0);
1934 split.x.p = new_enode(partition, 4, nparam);
1935 value_set_si(C->p[0][1+i], 1);
1936 Matrix *C2 = Matrix_Copy(C);
1937 EVALUE_SET_DOMAIN(split.x.p->arr[0],
1938 Constraints2Polyhedron(C2, MaxRays));
1939 Matrix_Free(C2);
1940 evalue_set_si(&split.x.p->arr[1], 1, 1);
1941 value_set_si(C->p[0][1+i], -1);
1942 value_set_si(C->p[0][1+nparam], -1);
1943 EVALUE_SET_DOMAIN(split.x.p->arr[2],
1944 Constraints2Polyhedron(C, MaxRays));
1945 evalue_set_si(&split.x.p->arr[3], 1, 1);
1946 emul(&split, EP);
1947 free_evalue_refs(&split);
1948 Matrix_Free(C);
1949 }
1950 reduce_evalue(EP);
1951 evalue_range_reduction(EP);
1952
1953 evalue_frac2floor(EP);
1954
1955 evalue *sum = esum(EP, nvar);
1956
1957 free_evalue_refs(EP);
1958 free(EP);
1959 EP = sum;
1960
1961 evalue_range_reduction(EP);
1962
1963 return EP;
1964 }
1965
1966 static evalue* split_sure(Polyhedron *P, Polyhedron *S,
1967 unsigned exist, unsigned nparam, unsigned MaxRays)
1968 {
1969 int nvar = P->Dimension - exist - nparam;
1970
1971 Matrix *M = Matrix_Alloc(exist, S->Dimension+2);
1972 for (int i = 0; i < exist; ++i)
1973 value_set_si(M->p[i][nvar+i+1], 1);
1974 Polyhedron *O = S;
1975 S = DomainAddRays(S, M, MaxRays);
1976 Polyhedron_Free(O);
1977 Polyhedron *F = DomainAddRays(P, M, MaxRays);
1978 Polyhedron *D = DomainDifference(F, S, MaxRays);
1979 O = D;
1980 D = Disjoint_Domain(D, 0, MaxRays);
1981 Polyhedron_Free(F);
1982 Domain_Free(O);
1983 Matrix_Free(M);
1984
1985 M = Matrix_Alloc(P->Dimension+1-exist, P->Dimension+1);
1986 for (int j = 0; j < nvar; ++j)
1987 value_set_si(M->p[j][j], 1);
1988 for (int j = 0; j < nparam+1; ++j)
1989 value_set_si(M->p[nvar+j][nvar+exist+j], 1);
1990 Polyhedron *T = Polyhedron_Image(S, M, MaxRays);
1991 evalue *EP = barvinok_enumerate_e(T, 0, nparam, MaxRays);
1992 Polyhedron_Free(S);
1993 Polyhedron_Free(T);
1994 Matrix_Free(M);
1995
1996 for (Polyhedron *Q = D; Q; Q = Q->next) {
1997 Polyhedron *N = Q->next;
1998 Q->next = 0;
1999 T = DomainIntersection(P, Q, MaxRays);
2000 evalue *E = barvinok_enumerate_e(T, exist, nparam, MaxRays);
2001 eadd(E, EP);
2002 free_evalue_refs(E);
2003 free(E);
2004 Polyhedron_Free(T);
2005 Q->next = N;
2006 }
2007 Domain_Free(D);
2008 return EP;
2009 }
2010
2011 static evalue* enumerate_sure(Polyhedron *P,
2012 unsigned exist, unsigned nparam, unsigned MaxRays)
2013 {
2014 int i;
2015 Polyhedron *S = P;
2016 int nvar = P->Dimension - exist - nparam;
2017 Value lcm;
2018 Value f;
2019 value_init(lcm);
2020 value_init(f);
2021
2022 for (i = 0; i < exist; ++i) {
2023 Matrix *M = Matrix_Alloc(S->NbConstraints, S->Dimension+2);
2024 int c = 0;
2025 value_set_si(lcm, 1);
2026 for (int j = 0; j < S->NbConstraints; ++j) {
2027 if (value_negz_p(S->Constraint[j][1+nvar+i]))
2028 continue;
2029 if (value_one_p(S->Constraint[j][1+nvar+i]))
2030 continue;
2031 value_lcm(lcm, S->Constraint[j][1+nvar+i], &lcm);
2032 }
2033
2034 for (int j = 0; j < S->NbConstraints; ++j) {
2035 if (value_negz_p(S->Constraint[j][1+nvar+i]))
2036 continue;
2037 if (value_one_p(S->Constraint[j][1+nvar+i]))
2038 continue;
2039 value_division(f, lcm, S->Constraint[j][1+nvar+i]);
2040 Vector_Scale(S->Constraint[j], M->p[c], f, S->Dimension+2);
2041 value_substract(M->p[c][S->Dimension+1],
2042 M->p[c][S->Dimension+1],
2043 lcm);
2044 value_increment(M->p[c][S->Dimension+1],
2045 M->p[c][S->Dimension+1]);
2046 ++c;
2047 }
2048 Polyhedron *O = S;
2049 S = AddConstraints(M->p[0], c, S, MaxRays);
2050 if (O != P)
2051 Polyhedron_Free(O);
2052 Matrix_Free(M);
2053 if (emptyQ(S)) {
2054 Polyhedron_Free(S);
2055 value_clear(lcm);
2056 value_clear(f);
2057 return 0;
2058 }
2059 }
2060 value_clear(lcm);
2061 value_clear(f);
2062
2063 #ifdef DEBUG_ER
2064 fprintf(stderr, "\nER: Sure\n");
2065 #endif /* DEBUG_ER */
2066
2067 return split_sure(P, S, exist, nparam, MaxRays);
2068 }
2069
2070 static evalue* enumerate_sure2(Polyhedron *P,
2071 unsigned exist, unsigned nparam, unsigned MaxRays)
2072 {
2073 int nvar = P->Dimension - exist - nparam;
2074 int r;
2075 for (r = 0; r < P->NbRays; ++r)
2076 if (value_one_p(P->Ray[r][0]) &&
2077 value_one_p(P->Ray[r][P->Dimension+1]))
2078 break;
2079
2080 if (r >= P->NbRays)
2081 return 0;
2082
2083 Matrix *M = Matrix_Alloc(nvar + 1 + nparam, P->Dimension+2);
2084 for (int i = 0; i < nvar; ++i)
2085 value_set_si(M->p[i][1+i], 1);
2086 for (int i = 0; i < nparam; ++i)
2087 value_set_si(M->p[i+nvar][1+nvar+exist+i], 1);
2088 Vector_Copy(P->Ray[r]+1+nvar, M->p[nvar+nparam]+1+nvar, exist);
2089 value_set_si(M->p[nvar+nparam][0], 1);
2090 value_set_si(M->p[nvar+nparam][P->Dimension+1], 1);
2091 Polyhedron * F = Rays2Polyhedron(M, MaxRays);
2092 Matrix_Free(M);
2093
2094 Polyhedron *I = DomainIntersection(F, P, MaxRays);
2095 Polyhedron_Free(F);
2096
2097 #ifdef DEBUG_ER
2098 fprintf(stderr, "\nER: Sure2\n");
2099 #endif /* DEBUG_ER */
2100
2101 return split_sure(P, I, exist, nparam, MaxRays);
2102 }
2103
2104 static evalue* enumerate_cyclic(Polyhedron *P,
2105 unsigned exist, unsigned nparam,
2106 evalue * EP, int r, int p, unsigned MaxRays)
2107 {
2108 int nvar = P->Dimension - exist - nparam;
2109
2110 /* If EP in its fractional maps only contains references
2111 * to the remainder parameter with appropriate coefficients
2112 * then we could in principle avoid adding existentially
2113 * quantified variables to the validity domains.
2114 * We'd have to replace the remainder by m { p/m }
2115 * and multiply with an appropriate factor that is one
2116 * only in the appropriate range.
2117 * This last multiplication can be avoided if EP
2118 * has a single validity domain with no (further)
2119 * constraints on the remainder parameter
2120 */
2121
2122 Matrix *CT = Matrix_Alloc(nparam+1, nparam+3);
2123 Matrix *M = Matrix_Alloc(1, 1+nparam+3);
2124 for (int j = 0; j < nparam; ++j)
2125 if (j != p)
2126 value_set_si(CT->p[j][j], 1);
2127 value_set_si(CT->p[p][nparam+1], 1);
2128 value_set_si(CT->p[nparam][nparam+2], 1);
2129 value_set_si(M->p[0][1+p], -1);
2130 value_absolute(M->p[0][1+nparam], P->Ray[0][1+nvar+exist+p]);
2131 value_set_si(M->p[0][1+nparam+1], 1);
2132 Polyhedron *CEq = Constraints2Polyhedron(M, 1);
2133 Matrix_Free(M);
2134 addeliminatedparams_enum(EP, CT, CEq, MaxRays, nparam);
2135 Polyhedron_Free(CEq);
2136 Matrix_Free(CT);
2137
2138 return EP;
2139 }
2140
2141 static void enumerate_vd_add_ray(evalue *EP, Matrix *Rays, unsigned MaxRays)
2142 {
2143 if (value_notzero_p(EP->d))
2144 return;
2145
2146 assert(EP->x.p->type == partition);
2147 assert(EP->x.p->pos == EVALUE_DOMAIN(EP->x.p->arr[0])->Dimension);
2148 for (int i = 0; i < EP->x.p->size/2; ++i) {
2149 Polyhedron *D = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
2150 Polyhedron *N = DomainAddRays(D, Rays, MaxRays);
2151 EVALUE_SET_DOMAIN(EP->x.p->arr[2*i], N);
2152 Domain_Free(D);
2153 }
2154 }
2155
2156 static evalue* enumerate_line(Polyhedron *P,
2157 unsigned exist, unsigned nparam, unsigned MaxRays)
2158 {
2159 if (P->NbBid == 0)
2160 return 0;
2161
2162 #ifdef DEBUG_ER
2163 fprintf(stderr, "\nER: Line\n");
2164 #endif /* DEBUG_ER */
2165
2166 int nvar = P->Dimension - exist - nparam;
2167 int i, j;
2168 for (i = 0; i < nparam; ++i)
2169 if (value_notzero_p(P->Ray[0][1+nvar+exist+i]))
2170 break;
2171 assert(i < nparam);
2172 for (j = i+1; j < nparam; ++j)
2173 if (value_notzero_p(P->Ray[0][1+nvar+exist+i]))
2174 break;
2175 assert(j >= nparam); // for now
2176
2177 Matrix *M = Matrix_Alloc(2, P->Dimension+2);
2178 value_set_si(M->p[0][0], 1);
2179 value_set_si(M->p[0][1+nvar+exist+i], 1);
2180 value_set_si(M->p[1][0], 1);
2181 value_set_si(M->p[1][1+nvar+exist+i], -1);
2182 value_absolute(M->p[1][1+P->Dimension], P->Ray[0][1+nvar+exist+i]);
2183 value_decrement(M->p[1][1+P->Dimension], M->p[1][1+P->Dimension]);
2184 Polyhedron *S = AddConstraints(M->p[0], 2, P, MaxRays);
2185 evalue *EP = barvinok_enumerate_e(S, exist, nparam, MaxRays);
2186 Polyhedron_Free(S);
2187 Matrix_Free(M);
2188
2189 return enumerate_cyclic(P, exist, nparam, EP, 0, i, MaxRays);
2190 }
2191
2192 static int single_param_pos(Polyhedron*P, unsigned exist, unsigned nparam,
2193 int r)
2194 {
2195 int nvar = P->Dimension - exist - nparam;
2196 if (First_Non_Zero(P->Ray[r]+1, nvar) != -1)
2197 return -1;
2198 int i = First_Non_Zero(P->Ray[r]+1+nvar+exist, nparam);
2199 if (i == -1)
2200 return -1;
2201 if (First_Non_Zero(P->Ray[r]+1+nvar+exist+1, nparam-i-1) != -1)
2202 return -1;
2203 return i;
2204 }
2205
2206 static evalue* enumerate_remove_ray(Polyhedron *P, int r,
2207 unsigned exist, unsigned nparam, unsigned MaxRays)
2208 {
2209 #ifdef DEBUG_ER
2210 fprintf(stderr, "\nER: RedundantRay\n");
2211 #endif /* DEBUG_ER */
2212
2213 Value one;
2214 value_init(one);
2215 value_set_si(one, 1);
2216 int len = P->NbRays-1;
2217 Matrix *M = Matrix_Alloc(2 * len, P->Dimension+2);
2218 Vector_Copy(P->Ray[0], M->p[0], r * (P->Dimension+2));
2219 Vector_Copy(P->Ray[r+1], M->p[r], (len-r) * (P->Dimension+2));
2220 for (int j = 0; j < P->NbRays; ++j) {
2221 if (j == r)
2222 continue;
2223 Vector_Combine(P->Ray[j], P->Ray[r], M->p[len+j-(j>r)],
2224 one, P->Ray[j][P->Dimension+1], P->Dimension+2);
2225 }
2226
2227 P = Rays2Polyhedron(M, MaxRays);
2228 Matrix_Free(M);
2229 evalue *EP = barvinok_enumerate_e(P, exist, nparam, MaxRays);
2230 Polyhedron_Free(P);
2231 value_clear(one);
2232
2233 return EP;
2234 }
2235
2236 static evalue* enumerate_redundant_ray(Polyhedron *P,
2237 unsigned exist, unsigned nparam, unsigned MaxRays)
2238 {
2239 assert(P->NbBid == 0);
2240 int nvar = P->Dimension - exist - nparam;
2241 Value m;
2242 value_init(m);
2243
2244 for (int r = 0; r < P->NbRays; ++r) {
2245 if (value_notzero_p(P->Ray[r][P->Dimension+1]))
2246 continue;
2247 int i1 = single_param_pos(P, exist, nparam, r);
2248 if (i1 == -1)
2249 continue;
2250 for (int r2 = r+1; r2 < P->NbRays; ++r2) {
2251 if (value_notzero_p(P->Ray[r2][P->Dimension+1]))
2252 continue;
2253 int i2 = single_param_pos(P, exist, nparam, r2);
2254 if (i2 == -1)
2255 continue;
2256 if (i1 != i2)
2257 continue;
2258
2259 value_division(m, P->Ray[r][1+nvar+exist+i1],
2260 P->Ray[r2][1+nvar+exist+i1]);
2261 value_multiply(m, m, P->Ray[r2][1+nvar+exist+i1]);
2262 /* r2 divides r => r redundant */
2263 if (value_eq(m, P->Ray[r][1+nvar+exist+i1])) {
2264 value_clear(m);
2265 return enumerate_remove_ray(P, r, exist, nparam, MaxRays);
2266 }
2267
2268 value_division(m, P->Ray[r2][1+nvar+exist+i1],
2269 P->Ray[r][1+nvar+exist+i1]);
2270 value_multiply(m, m, P->Ray[r][1+nvar+exist+i1]);
2271 /* r divides r2 => r2 redundant */
2272 if (value_eq(m, P->Ray[r2][1+nvar+exist+i1])) {
2273 value_clear(m);
2274 return enumerate_remove_ray(P, r2, exist, nparam, MaxRays);
2275 }
2276 }
2277 }
2278 value_clear(m);
2279 return 0;
2280 }
2281
2282 static Polyhedron *upper_bound(Polyhedron *P,
2283 int pos, Value *max, Polyhedron **R)
2284 {
2285 Value v;
2286 int r;
2287 value_init(v);
2288
2289 *R = 0;
2290 Polyhedron *N;
2291 Polyhedron *B = 0;
2292 for (Polyhedron *Q = P; Q; Q = N) {
2293 N = Q->next;
2294 for (r = 0; r < P->NbRays; ++r) {
2295 if (value_zero_p(P->Ray[r][P->Dimension+1]) &&
2296 value_pos_p(P->Ray[r][1+pos]))
2297 break;
2298 }
2299 if (r < P->NbRays) {
2300 Q->next = *R;
2301 *R = Q;
2302 continue;
2303 } else {
2304 Q->next = B;
2305 B = Q;
2306 }
2307 for (r = 0; r < P->NbRays; ++r) {
2308 if (value_zero_p(P->Ray[r][P->Dimension+1]))
2309 continue;
2310 mpz_fdiv_q(v, P->Ray[r][1+pos], P->Ray[r][1+P->Dimension]);
2311 if ((!Q->next && r == 0) || value_gt(v, *max))
2312 value_assign(*max, v);
2313 }
2314 }
2315 value_clear(v);
2316 return B;
2317 }
2318
2319 static evalue* enumerate_ray(Polyhedron *P,
2320 unsigned exist, unsigned nparam, unsigned MaxRays)
2321 {
2322 assert(P->NbBid == 0);
2323 int nvar = P->Dimension - exist - nparam;
2324
2325 int r;
2326 for (r = 0; r < P->NbRays; ++r)
2327 if (value_zero_p(P->Ray[r][P->Dimension+1]))
2328 break;
2329 if (r >= P->NbRays)
2330 return 0;
2331
2332 int r2;
2333 for (r2 = r+1; r2 < P->NbRays; ++r2)
2334 if (value_zero_p(P->Ray[r2][P->Dimension+1]))
2335 break;
2336 if (r2 < P->NbRays) {
2337 if (nvar > 0)
2338 return enumerate_sum(P, exist, nparam, MaxRays);
2339 }
2340
2341 #ifdef DEBUG_ER
2342 fprintf(stderr, "\nER: Ray\n");
2343 #endif /* DEBUG_ER */
2344
2345 Value m;
2346 Value one;
2347 value_init(m);
2348 value_init(one);
2349 value_set_si(one, 1);
2350 int i = single_param_pos(P, exist, nparam, r);
2351 assert(i != -1); // for now;
2352
2353 Matrix *M = Matrix_Alloc(P->NbRays, P->Dimension+2);
2354 for (int j = 0; j < P->NbRays; ++j) {
2355 Vector_Combine(P->Ray[j], P->Ray[r], M->p[j],
2356 one, P->Ray[j][P->Dimension+1], P->Dimension+2);
2357 }
2358 Polyhedron *S = Rays2Polyhedron(M, MaxRays);
2359 Matrix_Free(M);
2360 Polyhedron *D = DomainDifference(P, S, MaxRays);
2361 Polyhedron_Free(S);
2362 // Polyhedron_Print(stderr, P_VALUE_FMT, D);
2363 assert(value_pos_p(P->Ray[r][1+nvar+exist+i])); // for now
2364 Polyhedron *R;
2365 D = upper_bound(D, nvar+exist+i, &m, &R);
2366 assert(D);
2367 Domain_Free(D);
2368
2369 M = Matrix_Alloc(2, P->Dimension+2);
2370 value_set_si(M->p[0][0], 1);
2371 value_set_si(M->p[1][0], 1);
2372 value_set_si(M->p[0][1+nvar+exist+i], -1);
2373 value_set_si(M->p[1][1+nvar+exist+i], 1);
2374 value_assign(M->p[0][1+P->Dimension], m);
2375 value_oppose(M->p[1][1+P->Dimension], m);
2376 value_addto(M->p[1][1+P->Dimension], M->p[1][1+P->Dimension],
2377 P->Ray[r][1+nvar+exist+i]);
2378 value_decrement(M->p[1][1+P->Dimension], M->p[1][1+P->Dimension]);
2379 // Matrix_Print(stderr, P_VALUE_FMT, M);
2380 D = AddConstraints(M->p[0], 2, P, MaxRays);
2381 // Polyhedron_Print(stderr, P_VALUE_FMT, D);
2382 value_substract(M->p[0][1+P->Dimension], M->p[0][1+P->Dimension],
2383 P->Ray[r][1+nvar+exist+i]);
2384 // Matrix_Print(stderr, P_VALUE_FMT, M);
2385 S = AddConstraints(M->p[0], 1, P, MaxRays);
2386 // Polyhedron_Print(stderr, P_VALUE_FMT, S);
2387 Matrix_Free(M);
2388
2389 evalue *EP = barvinok_enumerate_e(D, exist, nparam, MaxRays);
2390 Polyhedron_Free(D);
2391 value_clear(one);
2392 value_clear(m);
2393
2394 if (value_notone_p(P->Ray[r][1+nvar+exist+i]))
2395 EP = enumerate_cyclic(P, exist, nparam, EP, r, i, MaxRays);
2396 else {
2397 M = Matrix_Alloc(1, nparam+2);
2398 value_set_si(M->p[0][0], 1);
2399 value_set_si(M->p[0][1+i], 1);
2400 enumerate_vd_add_ray(EP, M, MaxRays);
2401 Matrix_Free(M);
2402 }
2403
2404 if (!emptyQ(S)) {
2405 evalue *E = barvinok_enumerate_e(S, exist, nparam, MaxRays);
2406 eadd(E, EP);
2407 free_evalue_refs(E);
2408 free(E);
2409 }
2410 Polyhedron_Free(S);
2411
2412 if (R) {
2413 assert(nvar == 0);
2414 evalue *ER = enumerate_or(R, exist, nparam, MaxRays);
2415 eor(ER, EP);
2416 free_evalue_refs(ER);
2417 free(ER);
2418 }
2419
2420 return EP;
2421 }
2422
2423 static evalue* new_zero_ep()
2424 {
2425 evalue *EP;
2426 ALLOC(evalue, EP);
2427 value_init(EP->d);
2428 evalue_set_si(EP, 0, 1);
2429 return EP;
2430 }
2431
2432 static evalue* enumerate_vd(Polyhedron **PA,
2433 unsigned exist, unsigned nparam, unsigned MaxRays)
2434 {
2435 Polyhedron *P = *PA;
2436 int nvar = P->Dimension - exist - nparam;
2437 Param_Polyhedron *PP = NULL;
2438 Polyhedron *C = Universe_Polyhedron(nparam);
2439 Polyhedron *CEq;
2440 Matrix *CT;
2441 Polyhedron *PR = P;
2442 PP = Polyhedron2Param_SimplifiedDomain(&PR,C,MaxRays,&CEq,&CT);
2443 Polyhedron_Free(C);
2444
2445 int nd;
2446 Param_Domain *D, *last;
2447 Value c;
2448 value_init(c);
2449 for (nd = 0, D=PP->D; D; D=D->next, ++nd)
2450 ;
2451
2452 Polyhedron **VD = new Polyhedron_p[nd];
2453 Polyhedron **fVD = new Polyhedron_p[nd];
2454 for(nd = 0, D=PP->D; D; D=D->next) {
2455 Polyhedron *rVD = reduce_domain(D->Domain, CT, CEq,
2456 fVD, nd, MaxRays);
2457 if (!rVD)
2458 continue;
2459
2460 VD[nd++] = rVD;
2461 last = D;
2462 }
2463
2464 evalue *EP = 0;
2465
2466 if (nd == 0)
2467 EP = new_zero_ep();
2468
2469 /* This doesn't seem to have any effect */
2470 if (nd == 1) {
2471 Polyhedron *CA = align_context(VD[0], P->Dimension, MaxRays);
2472 Polyhedron *O = P;
2473 P = DomainIntersection(P, CA, MaxRays);
2474 if (O != *PA)
2475 Polyhedron_Free(O);
2476 Polyhedron_Free(CA);
2477 if (emptyQ(P))
2478 EP = new_zero_ep();
2479 }
2480
2481 if (!EP && CT->NbColumns != CT->NbRows) {
2482 Polyhedron *CEqr = DomainImage(CEq, CT, MaxRays);
2483 Polyhedron *CA = align_context(CEqr, PR->Dimension, MaxRays);
2484 Polyhedron *I = DomainIntersection(PR, CA, MaxRays);
2485 Polyhedron_Free(CEqr);
2486 Polyhedron_Free(CA);
2487 #ifdef DEBUG_ER
2488 fprintf(stderr, "\nER: Eliminate\n");
2489 #endif /* DEBUG_ER */
2490 nparam -= CT->NbColumns - CT->NbRows;
2491 EP = barvinok_enumerate_e(I, exist, nparam, MaxRays);
2492 nparam += CT->NbColumns - CT->NbRows;
2493 addeliminatedparams_enum(EP, CT, CEq, MaxRays, nparam);
2494 Polyhedron_Free(I);
2495 }
2496 if (PR != *PA)
2497 Polyhedron_Free(PR);
2498 PR = 0;
2499
2500 if (!EP && nd > 1) {
2501 #ifdef DEBUG_ER
2502 fprintf(stderr, "\nER: VD\n");
2503 #endif /* DEBUG_ER */
2504 for (int i = 0; i < nd; ++i) {
2505 Polyhedron *CA = align_context(VD[i], P->Dimension, MaxRays);
2506 Polyhedron *I = DomainIntersection(P, CA, MaxRays);
2507
2508 if (i == 0)
2509 EP = barvinok_enumerate_e(I, exist, nparam, MaxRays);
2510 else {
2511 evalue *E = barvinok_enumerate_e(I, exist, nparam, MaxRays);
2512 eadd(E, EP);
2513 free_evalue_refs(E);
2514 free(E);
2515 }
2516 Polyhedron_Free(I);
2517 Polyhedron_Free(CA);
2518 }
2519 }
2520
2521 for (int i = 0; i < nd; ++i) {
2522 Polyhedron_Free(VD[i]);
2523 Polyhedron_Free(fVD[i]);
2524 }
2525 delete [] VD;
2526 delete [] fVD;
2527 value_clear(c);
2528
2529 if (!EP && nvar == 0) {
2530 Value f;
2531 value_init(f);
2532 Param_Vertices *V, *V2;
2533 Matrix* M = Matrix_Alloc(1, P->Dimension+2);
2534
2535 FORALL_PVertex_in_ParamPolyhedron(V, last, PP) {
2536 bool found = false;
2537 FORALL_PVertex_in_ParamPolyhedron(V2, last, PP) {
2538 if (V == V2) {
2539 found = true;
2540 continue;
2541 }
2542 if (!found)
2543 continue;
2544 for (int i = 0; i < exist; ++i) {
2545 value_oppose(f, V->Vertex->p[i][nparam+1]);
2546 Vector_Combine(V->Vertex->p[i],
2547 V2->Vertex->p[i],
2548 M->p[0] + 1 + nvar + exist,
2549 V2->Vertex->p[i][nparam+1],
2550 f,
2551 nparam+1);
2552 int j;
2553 for (j = 0; j < nparam; ++j)
2554 if (value_notzero_p(M->p[0][1+nvar+exist+j]))
2555 break;
2556 if (j >= nparam)
2557 continue;
2558 ConstraintSimplify(M->p[0], M->p[0],
2559 P->Dimension+2, &f);
2560 value_set_si(M->p[0][0], 0);
2561 Polyhedron *para = AddConstraints(M->p[0], 1, P,
2562 MaxRays);
2563 if (emptyQ(para)) {
2564 Polyhedron_Free(para);
2565 continue;
2566 }
2567 Polyhedron *pos, *neg;
2568 value_set_si(M->p[0][0], 1);
2569 value_decrement(M->p[0][P->Dimension+1],
2570 M->p[0][P->Dimension+1]);
2571 neg = AddConstraints(M->p[0], 1, P, MaxRays);
2572 value_set_si(f, -1);
2573 Vector_Scale(M->p[0]+1, M->p[0]+1, f,
2574 P->Dimension+1);
2575 value_decrement(M->p[0][P->Dimension+1],
2576 M->p[0][P->Dimension+1]);
2577 value_decrement(M->p[0][P->Dimension+1],
2578 M->p[0][P->Dimension+1]);
2579 pos = AddConstraints(M->p[0], 1, P, MaxRays);
2580 if (emptyQ(neg) && emptyQ(pos)) {
2581 Polyhedron_Free(para);
2582 Polyhedron_Free(pos);
2583 Polyhedron_Free(neg);
2584 continue;
2585 }
2586 #ifdef DEBUG_ER
2587 fprintf(stderr, "\nER: Order\n");
2588 #endif /* DEBUG_ER */
2589 EP = barvinok_enumerate_e(para, exist, nparam, MaxRays);
2590 evalue *E;
2591 if (!emptyQ(pos)) {
2592 E = barvinok_enumerate_e(pos, exist, nparam, MaxRays);
2593 eadd(E, EP);
2594 free_evalue_refs(E);
2595 free(E);
2596 }
2597 if (!emptyQ(neg)) {
2598 E = barvinok_enumerate_e(neg, exist, nparam, MaxRays);
2599 eadd(E, EP);
2600 free_evalue_refs(E);
2601 free(E);
2602 }
2603 Polyhedron_Free(para);
2604 Polyhedron_Free(pos);
2605 Polyhedron_Free(neg);
2606 break;
2607 }
2608 if (EP)
2609 break;
2610 } END_FORALL_PVertex_in_ParamPolyhedron;
2611 if (EP)
2612 break;
2613 } END_FORALL_PVertex_in_ParamPolyhedron;
2614
2615 if (!EP) {
2616 /* Search for vertex coordinate to split on */
2617 /* First look for one independent of the parameters */
2618 FORALL_PVertex_in_ParamPolyhedron(V, last, PP) {
2619 for (int i = 0; i < exist; ++i) {
2620 int j;
2621 for (j = 0; j < nparam; ++j)
2622 if (value_notzero_p(V->Vertex->p[i][j]))
2623 break;
2624 if (j < nparam)
2625 continue;
2626 value_set_si(M->p[0][0], 1);
2627 Vector_Set(M->p[0]+1, 0, nvar+exist);
2628 Vector_Copy(V->Vertex->p[i],
2629 M->p[0] + 1 + nvar + exist, nparam+1);
2630 value_oppose(M->p[0][1+nvar+i],
2631 V->Vertex->p[i][nparam+1]);
2632
2633 Polyhedron *pos, *neg;
2634 value_set_si(M->p[0][0], 1);
2635 value_decrement(M->p[0][P->Dimension+1],
2636 M->p[0][P->Dimension+1]);
2637 neg = AddConstraints(M->p[0], 1, P, MaxRays);
2638 value_set_si(f, -1);
2639 Vector_Scale(M->p[0]+1, M->p[0]+1, f,
2640 P->Dimension+1);
2641 value_decrement(M->p[0][P->Dimension+1],
2642 M->p[0][P->Dimension+1]);
2643 value_decrement(M->p[0][P->Dimension+1],
2644 M->p[0][P->Dimension+1]);
2645 pos = AddConstraints(M->p[0], 1, P, MaxRays);
2646 if (emptyQ(neg) || emptyQ(pos)) {
2647 Polyhedron_Free(pos);
2648 Polyhedron_Free(neg);
2649 continue;
2650 }
2651 Polyhedron_Free(pos);
2652 value_increment(M->p[0][P->Dimension+1],
2653 M->p[0][P->Dimension+1]);
2654 pos = AddConstraints(M->p[0], 1, P, MaxRays);
2655 #ifdef DEBUG_ER
2656 fprintf(stderr, "\nER: Vertex\n");
2657 #endif /* DEBUG_ER */
2658 pos->next = neg;
2659 EP = enumerate_or(pos, exist, nparam, MaxRays);
2660 break;
2661 }
2662 if (EP)
2663 break;
2664 } END_FORALL_PVertex_in_ParamPolyhedron;
2665 }
2666
2667 if (!EP) {
2668 /* Search for vertex coordinate to split on */
2669 /* Now look for one that depends on the parameters */
2670 FORALL_PVertex_in_ParamPolyhedron(V, last, PP) {
2671 for (int i = 0; i < exist; ++i) {
2672 value_set_si(M->p[0][0], 1);
2673 Vector_Set(M->p[0]+1, 0, nvar+exist);
2674 Vector_Copy(V->Vertex->p[i],
2675 M->p[0] + 1 + nvar + exist, nparam+1);
2676 value_oppose(M->p[0][1+nvar+i],
2677 V->Vertex->p[i][nparam+1]);
2678
2679 Polyhedron *pos, *neg;
2680 value_set_si(M->p[0][0], 1);
2681 value_decrement(M->p[0][P->Dimension+1],
2682 M->p[0][P->Dimension+1]);
2683 neg = AddConstraints(M->p[0], 1, P, MaxRays);
2684 value_set_si(f, -1);
2685 Vector_Scale(M->p[0]+1, M->p[0]+1, f,
2686 P->Dimension+1);
2687 value_decrement(M->p[0][P->Dimension+1],
2688 M->p[0][P->Dimension+1]);
2689 value_decrement(M->p[0][P->Dimension+1],
2690 M->p[0][P->Dimension+1]);
2691 pos = AddConstraints(M->p[0], 1, P, MaxRays);
2692 if (emptyQ(neg) || emptyQ(pos)) {
2693 Polyhedron_Free(pos);
2694 Polyhedron_Free(neg);
2695 continue;
2696 }
2697 Polyhedron_Free(pos);
2698 value_increment(M->p[0][P->Dimension+1],
2699 M->p[0][P->Dimension+1]);
2700 pos = AddConstraints(M->p[0], 1, P, MaxRays);
2701 #ifdef DEBUG_ER
2702 fprintf(stderr, "\nER: ParamVertex\n");
2703 #endif /* DEBUG_ER */
2704 pos->next = neg;
2705 EP = enumerate_or(pos, exist, nparam, MaxRays);
2706 break;
2707 }
2708 if (EP)
2709 break;
2710 } END_FORALL_PVertex_in_ParamPolyhedron;
2711 }
2712
2713 Matrix_Free(M);
2714 value_clear(f);
2715 }
2716
2717 if (CEq)
2718 Polyhedron_Free(CEq);
2719 if (CT)
2720 Matrix_Free(CT);
2721 if (PP)
2722 Param_Polyhedron_Free(PP);
2723 *PA = P;
2724
2725 return EP;
2726 }
2727
2728 #ifndef HAVE_PIPLIB
2729 evalue *barvinok_enumerate_pip(Polyhedron *P,
2730 unsigned exist, unsigned nparam, unsigned MaxRays)
2731 {
2732 return 0;
2733 }
2734 #else
2735 evalue *barvinok_enumerate_pip(Polyhedron *P,
2736 unsigned exist, unsigned nparam, unsigned MaxRays)
2737 {
2738 int nvar = P->Dimension - exist - nparam;
2739 evalue *EP = new_zero_ep();
2740 Polyhedron *Q, *N, *T = 0;
2741 Value min, tmp;
2742 value_init(min);
2743 value_init(tmp);
2744
2745 #ifdef DEBUG_ER
2746 fprintf(stderr, "\nER: PIP\n");
2747 #endif /* DEBUG_ER */
2748
2749 for (int i = 0; i < P->Dimension; ++i) {
2750 bool pos = false;
2751 bool neg = false;
2752 bool posray = false;
2753 bool negray = false;
2754 value_set_si(min, 0);
2755 for (int j = 0; j < P->NbRays; ++j) {
2756 if (value_pos_p(P->Ray[j][1+i])) {
2757 pos = true;
2758 if (value_zero_p(P->Ray[j][1+P->Dimension]))
2759 posray = true;
2760 } else if (value_neg_p(P->Ray[j][1+i])) {
2761 neg = true;
2762 if (value_zero_p(P->Ray[j][1+P->Dimension]))
2763 negray = true;
2764 else {
2765 mpz_fdiv_q(tmp,
2766 P->Ray[j][1+i], P->Ray[j][1+P->Dimension]);
2767 if (value_lt(tmp, min))
2768 value_assign(min, tmp);
2769 }
2770 }
2771 }
2772 if (pos && neg) {
2773 assert(!(posray && negray));
2774 assert(!negray); // for now
2775 Polyhedron *O = T ? T : P;
2776 /* shift by a safe amount */
2777 Matrix *M = Matrix_Alloc(O->NbRays, O->Dimension+2);
2778 Vector_Copy(O->Ray[0], M->p[0], O->NbRays * (O->Dimension+2));
2779 for (int j = 0; j < P->NbRays; ++j) {
2780 if (value_notzero_p(M->p[j][1+P->Dimension])) {
2781 value_multiply(tmp, min, M->p[j][1+P->Dimension]);
2782 value_substract(M->p[j][1+i], M->p[j][1+i], tmp);
2783 }
2784 }
2785 if (T)
2786 Polyhedron_Free(T);
2787 T = Rays2Polyhedron(M, MaxRays);
2788 Matrix_Free(M);
2789 } else if (neg) {
2790 /* negating a parameter requires that we substitute in the
2791 * sign again afterwards.
2792 * Disallow for now.
2793 */
2794 assert(i < nvar+exist);
2795 if (!T)
2796 T = Polyhedron_Copy(P);
2797 for (int j = 0; j < T->NbRays; ++j)
2798 value_oppose(T->Ray[j][1+i], T->Ray[j][1+i]);
2799 for (int j = 0; j < T->NbConstraints; ++j)
2800 value_oppose(T->Constraint[j][1+i], T->Constraint[j][1+i]);
2801 }
2802 }
2803 value_clear(min);
2804 value_clear(tmp);
2805
2806 Polyhedron *D = pip_lexmin(T ? T : P, exist, nparam);
2807 for (Q = D; Q; Q = N) {
2808 N = Q->next;
2809 Q->next = 0;
2810 evalue *E;
2811 exist = Q->Dimension - nvar - nparam;
2812 E = barvinok_enumerate_e(Q, exist, nparam, MaxRays);
2813 Polyhedron_Free(Q);
2814 eadd(E, EP);
2815 free_evalue_refs(E);
2816 free(E);
2817 }
2818
2819 if (T)
2820 Polyhedron_Free(T);
2821
2822 return EP;
2823 }
2824 #endif
2825
2826
2827 static bool is_single(Value *row, int pos, int len)
2828 {
2829 return First_Non_Zero(row, pos) == -1 &&
2830 First_Non_Zero(row+pos+1, len-pos-1) == -1;
2831 }
2832
2833 static evalue* barvinok_enumerate_e_r(Polyhedron *P,
2834 unsigned exist, unsigned nparam, unsigned MaxRays);
2835
2836 #ifdef DEBUG_ER
2837 static int er_level = 0;
2838
2839 evalue* barvinok_enumerate_e(Polyhedron *P,
2840 unsigned exist, unsigned nparam, unsigned MaxRays)
2841 {
2842 fprintf(stderr, "\nER: level %i\n", er_level);
2843 int nvar = P->Dimension - exist - nparam;
2844 fprintf(stderr, "%d %d %d\n", nvar, exist, nparam);
2845
2846 Polyhedron_Print(stderr, P_VALUE_FMT, P);
2847 ++er_level;
2848 P = DomainConstraintSimplify(Polyhedron_Copy(P), MaxRays);
2849 evalue *EP = barvinok_enumerate_e_r(P, exist, nparam, MaxRays);
2850 Polyhedron_Free(P);
2851 --er_level;
2852 return EP;
2853 }
2854 #else
2855 evalue* barvinok_enumerate_e(Polyhedron *P,
2856 unsigned exist, unsigned nparam, unsigned MaxRays)
2857 {
2858 P = DomainConstraintSimplify(Polyhedron_Copy(P), MaxRays);
2859 evalue *EP = barvinok_enumerate_e_r(P, exist, nparam, MaxRays);
2860 Polyhedron_Free(P);
2861 return EP;
2862 }
2863 #endif
2864
2865 static evalue* barvinok_enumerate_e_r(Polyhedron *P,
2866 unsigned exist, unsigned nparam, unsigned MaxRays)
2867 {
2868 if (exist == 0) {
2869 Polyhedron *U = Universe_Polyhedron(nparam);
2870 evalue *EP = barvinok_enumerate_ev(P, U, MaxRays);
2871 //char *param_name[] = {"P", "Q", "R", "S", "T" };
2872 //print_evalue(stdout, EP, param_name);
2873 Polyhedron_Free(U);
2874 return EP;
2875 }
2876
2877 int nvar = P->Dimension - exist - nparam;
2878 int len = P->Dimension + 2;
2879
2880 if (emptyQ(P))
2881 return new_zero_ep();
2882
2883 if (nvar == 0 && nparam == 0) {
2884 evalue *EP = new_zero_ep();
2885 barvinok_count(P, &EP->x.n, MaxRays);
2886 if (value_pos_p(EP->x.n))
2887 value_set_si(EP->x.n, 1);
2888 return EP;
2889 }
2890
2891 int r;
2892 for (r = 0; r < P->NbRays; ++r)
2893 if (value_zero_p(P->Ray[r][0]) ||
2894 value_zero_p(P->Ray[r][P->Dimension+1])) {
2895 int i;
2896 for (i = 0; i < nvar; ++i)
2897 if (value_notzero_p(P->Ray[r][i+1]))
2898 break;
2899 if (i >= nvar)
2900 continue;
2901 for (i = nvar + exist; i < nvar + exist + nparam; ++i)
2902 if (value_notzero_p(P->Ray[r][i+1]))
2903 break;
2904 if (i >= nvar + exist + nparam)
2905 break;
2906 }
2907 if (r < P->NbRays) {
2908 evalue *EP = new_zero_ep();
2909 value_set_si(EP->x.n, -1);
2910 return EP;
2911 }
2912
2913 int first;
2914 for (r = 0; r < P->NbEq; ++r)
2915 if ((first = First_Non_Zero(P->Constraint[r]+1+nvar, exist)) != -1)
2916 break;
2917 if (r < P->NbEq) {
2918 if (First_Non_Zero(P->Constraint[r]+1+nvar+first+1,
2919 exist-first-1) != -1) {
2920 Polyhedron *T = rotate_along(P, r, nvar, exist, MaxRays);
2921 #ifdef DEBUG_ER
2922 fprintf(stderr, "\nER: Equality\n");
2923 #endif /* DEBUG_ER */
2924 evalue *EP = barvinok_enumerate_e(T, exist-1, nparam, MaxRays);
2925 Polyhedron_Free(T);
2926 return EP;
2927 } else {
2928 #ifdef DEBUG_ER
2929 fprintf(stderr, "\nER: Fixed\n");
2930 #endif /* DEBUG_ER */
2931 if (first == 0)
2932 return barvinok_enumerate_e(P, exist-1, nparam, MaxRays);
2933 else {
2934 Polyhedron *T = Polyhedron_Copy(P);
2935 SwapColumns(T, nvar+1, nvar+1+first);
2936 evalue *EP = barvinok_enumerate_e(T, exist-1, nparam, MaxRays);
2937 Polyhedron_Free(T);
2938 return EP;
2939 }
2940 }
2941 }
2942
2943 Vector *row = Vector_Alloc(len);
2944 value_set_si(row->p[0], 1);
2945
2946 Value f;
2947 value_init(f);
2948
2949 enum constraint* info = new constraint[exist];
2950 for (int i = 0; i < exist; ++i) {
2951 info[i] = ALL_POS;
2952 for (int l = P->NbEq; l < P->NbConstraints; ++l) {
2953 if (value_negz_p(P->Constraint[l][nvar+i+1]))
2954 continue;
2955 bool l_parallel = is_single(P->Constraint[l]+nvar+1, i, exist);
2956 for (int u = P->NbEq; u < P->NbConstraints; ++u) {
2957 if (value_posz_p(P->Constraint[u][nvar+i+1]))
2958 continue;
2959 bool lu_parallel = l_parallel ||
2960 is_single(P->Constraint[u]+nvar+1, i, exist);
2961 value_oppose(f, P->Constraint[u][nvar+i+1]);
2962 Vector_Combine(P->Constraint[l]+1, P->Constraint[u]+1, row->p+1,
2963 f, P->Constraint[l][nvar+i+1], len-1);
2964 if (!(info[i] & INDEPENDENT)) {
2965 int j;
2966 for (j = 0; j < exist; ++j)
2967 if (j != i && value_notzero_p(row->p[nvar+j+1]))
2968 break;
2969 if (j == exist) {
2970 //printf("independent: i: %d, l: %d, u: %d\n", i, l, u);
2971 info[i] = (constraint)(info[i] | INDEPENDENT);
2972 }
2973 }
2974 if (info[i] & ALL_POS) {
2975 value_addto(row->p[len-1], row->p[len-1],
2976 P->Constraint[l][nvar+i+1]);
2977 value_addto(row->p[len-1], row->p[len-1], f);
2978 value_multiply(f, f, P->Constraint[l][nvar+i+1]);
2979 value_substract(row->p[len-1], row->p[len-1], f);
2980 value_decrement(row->p[len-1], row->p[len-1]);
2981 ConstraintSimplify(row->p, row->p, len, &f);
2982 value_set_si(f, -1);
2983 Vector_Scale(row->p+1, row->p+1, f, len-1);
2984 value_decrement(row->p[len-1], row->p[len-1]);
2985 Polyhedron *T = AddConstraints(row->p, 1, P, MaxRays);
2986 if (!emptyQ(T)) {
2987 //printf("not all_pos: i: %d, l: %d, u: %d\n", i, l, u);
2988 info[i] = (constraint)(info[i] ^ ALL_POS);
2989 }
2990 //puts("pos remainder");
2991 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
2992 Polyhedron_Free(T);
2993 }
2994 if (!(info[i] & ONE_NEG)) {
2995 if (lu_parallel) {
2996 negative_test_constraint(P->Constraint[l],
2997 P->Constraint[u],
2998 row->p, nvar+i, len, &f);
2999 oppose_constraint(row->p, len, &f);
3000 Polyhedron *T = AddConstraints(row->p, 1, P, MaxRays);
3001 if (emptyQ(T)) {
3002 //printf("one_neg i: %d, l: %d, u: %d\n", i, l, u);
3003 info[i] = (constraint)(info[i] | ONE_NEG);
3004 }
3005 //puts("neg remainder");
3006 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
3007 Polyhedron_Free(T);
3008 }
3009 }
3010 if (!(info[i] & ALL_POS) && (info[i] & ONE_NEG))
3011 goto next;
3012 }
3013 }
3014 if (info[i] & ALL_POS)
3015 break;
3016 next:
3017 ;
3018 }
3019
3020 /*
3021 for (int i = 0; i < exist; ++i)
3022 printf("%i: %i\n", i, info[i]);
3023 */
3024 for (int i = 0; i < exist; ++i)
3025 if (info[i] & ALL_POS) {
3026 #ifdef DEBUG_ER
3027 fprintf(stderr, "\nER: Positive\n");
3028 #endif /* DEBUG_ER */
3029 // Eliminate
3030 // Maybe we should chew off some of the fat here
3031 Matrix *M = Matrix_Alloc(P->Dimension, P->Dimension+1);
3032 for (int j = 0; j < P->Dimension; ++j)
3033 value_set_si(M->p[j][j + (j >= i+nvar)], 1);
3034 Polyhedron *T = Polyhedron_Image(P, M, MaxRays);
3035 Matrix_Free(M);
3036 evalue *EP = barvinok_enumerate_e(T, exist-1, nparam, MaxRays);
3037 Polyhedron_Free(T);
3038 value_clear(f);
3039 Vector_Free(row);
3040 delete [] info;
3041 return EP;
3042 }
3043 for (int i = 0; i < exist; ++i)
3044 if (info[i] & ONE_NEG) {
3045 #ifdef DEBUG_ER
3046 fprintf(stderr, "\nER: Negative\n");
3047 #endif /* DEBUG_ER */
3048 Vector_Free(row);
3049 value_clear(f);
3050 delete [] info;
3051 if (i == 0)
3052 return barvinok_enumerate_e(P, exist-1, nparam, MaxRays);
3053 else {
3054 Polyhedron *T = Polyhedron_Copy(P);
3055 SwapColumns(T, nvar+1, nvar+1+i);
3056 evalue *EP = barvinok_enumerate_e(T, exist-1, nparam, MaxRays);
3057 Polyhedron_Free(T);
3058 return EP;
3059 }
3060 }
3061 for (int i = 0; i < exist; ++i)
3062 if (info[i] & INDEPENDENT) {
3063 Polyhedron *pos, *neg;
3064
3065 /* Find constraint again and split off negative part */
3066
3067 if (SplitOnVar(P, i, nvar, len, exist, MaxRays,
3068 row, f, true, &pos, &neg)) {
3069 #ifdef DEBUG_ER
3070 fprintf(stderr, "\nER: Split\n");
3071 #endif /* DEBUG_ER */
3072
3073 evalue *EP =
3074 barvinok_enumerate_e(neg, exist-1, nparam, MaxRays);
3075 evalue *E =
3076 barvinok_enumerate_e(pos, exist, nparam, MaxRays);
3077 eadd(E, EP);
3078 free_evalue_refs(E);
3079 free(E);
3080 Polyhedron_Free(neg);
3081 Polyhedron_Free(pos);
3082 value_clear(f);
3083 Vector_Free(row);
3084 delete [] info;
3085 return EP;
3086 }
3087 }
3088 delete [] info;
3089
3090 Polyhedron *O = P;
3091 Polyhedron *F;
3092
3093 evalue *EP;
3094
3095 EP = enumerate_line(P, exist, nparam, MaxRays);
3096 if (EP)
3097 goto out;
3098
3099 EP = barvinok_enumerate_pip(P, exist, nparam, MaxRays);
3100 if (EP)
3101 goto out;
3102
3103 EP = enumerate_redundant_ray(P, exist, nparam, MaxRays);
3104 if (EP)
3105 goto out;
3106
3107 EP = enumerate_sure(P, exist, nparam, MaxRays);
3108 if (EP)
3109 goto out;
3110
3111 EP = enumerate_ray(P, exist, nparam, MaxRays);
3112 if (EP)
3113 goto out;
3114
3115 EP = enumerate_sure2(P, exist, nparam, MaxRays);
3116 if (EP)
3117 goto out;
3118
3119 F = unfringe(P, MaxRays);
3120 if (!PolyhedronIncludes(F, P)) {
3121 #ifdef DEBUG_ER
3122 fprintf(stderr, "\nER: Fringed\n");
3123 #endif /* DEBUG_ER */
3124 EP = barvinok_enumerate_e(F, exist, nparam, MaxRays);
3125 Polyhedron_Free(F);
3126 goto out;
3127 }
3128 Polyhedron_Free(F);
3129
3130 if (nparam)
3131 EP = enumerate_vd(&P, exist, nparam, MaxRays);
3132 if (EP)
3133 goto out2;
3134
3135 if (nvar != 0) {
3136 EP = enumerate_sum(P, exist, nparam, MaxRays);
3137 goto out2;
3138 }
3139
3140 assert(nvar == 0);
3141
3142 int i;
3143 Polyhedron *pos, *neg;
3144 for (i = 0; i < exist; ++i)
3145 if (SplitOnVar(P, i, nvar, len, exist, MaxRays,
3146 row, f, false, &pos, &neg))
3147 break;
3148
3149 assert (i < exist);
3150
3151 pos->next = neg;
3152 EP = enumerate_or(pos, exist, nparam, MaxRays);
3153
3154 out2:
3155 if (O != P)
3156 Polyhedron_Free(P);
3157
3158 out:
3159 value_clear(f);
3160 Vector_Free(row);
3161 return EP;
3162 }
lattice point enumeration library