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