?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
more partial ray removal (largely untested
[barvinok.git] / barvinok.cc
blob0c2290db7017bc123256cfd7632d2f78fa7322a3
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 }
15 #include "config.h"
16 #include <barvinok.h>
17
18 #ifdef NTL_STD_CXX
19 using namespace NTL;
20 #endif
21 using std::cout;
22 using std::endl;
23 using std::vector;
24 using std::deque;
25 using std::string;
26 using std::ostringstream;
27
28 #define ALLOC(p) (((long *) (p))[0])
29 #define SIZE(p) (((long *) (p))[1])
30 #define DATA(p) ((mp_limb_t *) (((long *) (p)) + 2))
31
32 static void value2zz(Value v, ZZ& z)
33 {
34 int sa = v[0]._mp_size;
35 int abs_sa = sa < 0 ? -sa : sa;
36
37 _ntl_gsetlength(&z.rep, abs_sa);
38 mp_limb_t * adata = DATA(z.rep);
39 for (int i = 0; i < abs_sa; ++i)
40 adata[i] = v[0]._mp_d[i];
41 SIZE(z.rep) = sa;
42 }
43
44 static void zz2value(ZZ& z, Value& v)
45 {
46 if (!z.rep) {
47 value_set_si(v, 0);
48 return;
49 }
50
51 int sa = SIZE(z.rep);
52 int abs_sa = sa < 0 ? -sa : sa;
53
54 mp_limb_t * adata = DATA(z.rep);
55 _mpz_realloc(v, abs_sa);
56 for (int i = 0; i < abs_sa; ++i)
57 v[0]._mp_d[i] = adata[i];
58 v[0]._mp_size = sa;
59 }
60
61 #undef ALLOC
62 #define ALLOC(p) p = (typeof(p))malloc(sizeof(*p))
63
64 /*
65 * We just ignore the last column and row
66 * If the final element is not equal to one
67 * then the result will actually be a multiple of the input
68 */
69 static void matrix2zz(Matrix *M, mat_ZZ& m, unsigned nr, unsigned nc)
70 {
71 m.SetDims(nr, nc);
72
73 for (int i = 0; i < nr; ++i) {
74 // assert(value_one_p(M->p[i][M->NbColumns - 1]));
75 for (int j = 0; j < nc; ++j) {
76 value2zz(M->p[i][j], m[i][j]);
77 }
78 }
79 }
80
81 static void values2zz(Value *p, vec_ZZ& v, int len)
82 {
83 v.SetLength(len);
84
85 for (int i = 0; i < len; ++i) {
86 value2zz(p[i], v[i]);
87 }
88 }
89
90 /*
91 */
92 static void zz2values(vec_ZZ& v, Value *p)
93 {
94 for (int i = 0; i < v.length(); ++i)
95 zz2value(v[i], p[i]);
96 }
97
98 static void rays(mat_ZZ& r, Polyhedron *C)
99 {
100 unsigned dim = C->NbRays - 1; /* don't count zero vertex */
101 assert(C->NbRays - 1 == C->Dimension);
102 r.SetDims(dim, dim);
103 ZZ tmp;
104
105 int i, c;
106 for (i = 0, c = 0; i < dim; ++i)
107 if (value_zero_p(C->Ray[i][dim+1])) {
108 for (int j = 0; j < dim; ++j) {
109 value2zz(C->Ray[i][j+1], tmp);
110 r[j][c] = tmp;
111 }
112 ++c;
113 }
114 }
115
116 static Matrix * rays(Polyhedron *C)
117 {
118 unsigned dim = C->NbRays - 1; /* don't count zero vertex */
119 assert(C->NbRays - 1 == C->Dimension);
120
121 Matrix *M = Matrix_Alloc(dim+1, dim+1);
122 assert(M);
123
124 int i, c;
125 for (i = 0, c = 0; i <= dim && c < dim; ++i)
126 if (value_zero_p(C->Ray[i][dim+1])) {
127 Vector_Copy(C->Ray[i] + 1, M->p[c], dim);
128 value_set_si(M->p[c++][dim], 0);
129 }
130 assert(c == dim);
131 value_set_si(M->p[dim][dim], 1);
132
133 return M;
134 }
135
136 static Matrix * rays2(Polyhedron *C)
137 {
138 unsigned dim = C->NbRays - 1; /* don't count zero vertex */
139 assert(C->NbRays - 1 == C->Dimension);
140
141 Matrix *M = Matrix_Alloc(dim, dim);
142 assert(M);
143
144 int i, c;
145 for (i = 0, c = 0; i <= dim && c < dim; ++i)
146 if (value_zero_p(C->Ray[i][dim+1]))
147 Vector_Copy(C->Ray[i] + 1, M->p[c++], dim);
148 assert(c == dim);
149
150 return M;
151 }
152
153 /*
154 * Returns the largest absolute value in the vector
155 */
156 static ZZ max(vec_ZZ& v)
157 {
158 ZZ max = abs(v[0]);
159 for (int i = 1; i < v.length(); ++i)
160 if (abs(v[i]) > max)
161 max = abs(v[i]);
162 return max;
163 }
164
165 class cone {
166 public:
167 cone(Matrix *M) {
168 Cone = 0;
169 Rays = Matrix_Copy(M);
170 set_det();
171 }
172 cone(Polyhedron *C) {
173 Cone = Polyhedron_Copy(C);
174 Rays = rays(C);
175 set_det();
176 }
177 void set_det() {
178 mat_ZZ A;
179 matrix2zz(Rays, A, Rays->NbRows - 1, Rays->NbColumns - 1);
180 det = determinant(A);
181 Value v;
182 value_init(v);
183 zz2value(det, v);
184 value_clear(v);
185 }
186
187 Vector* short_vector(vec_ZZ& lambda) {
188 Matrix *M = Matrix_Copy(Rays);
189 Matrix *inv = Matrix_Alloc(M->NbRows, M->NbColumns);
190 int ok = Matrix_Inverse(M, inv);
191 assert(ok);
192 Matrix_Free(M);
193
194 ZZ det2;
195 mat_ZZ B;
196 mat_ZZ U;
197 matrix2zz(inv, B, inv->NbRows - 1, inv->NbColumns - 1);
198 long r = LLL(det2, B, U);
199
200 ZZ min = max(B[0]);
201 int index = 0;
202 for (int i = 1; i < B.NumRows(); ++i) {
203 ZZ tmp = max(B[i]);
204 if (tmp < min) {
205 min = tmp;
206 index = i;
207 }
208 }
209
210 Matrix_Free(inv);
211
212 lambda = B[index];
213
214 Vector *z = Vector_Alloc(U[index].length()+1);
215 assert(z);
216 zz2values(U[index], z->p);
217 value_set_si(z->p[U[index].length()], 0);
218
219 Value tmp;
220 value_init(tmp);
221 Polyhedron *C = poly();
222 int i;
223 for (i = 0; i < C->NbConstraints; ++i) {
224 Inner_Product(z->p, C->Constraint[i]+1, z->Size-1, &tmp);
225 if (value_pos_p(tmp))
226 break;
227 }
228 if (i == C->NbConstraints) {
229 value_set_si(tmp, -1);
230 Vector_Scale(z->p, z->p, tmp, z->Size-1);
231 }
232 value_clear(tmp);
233 return z;
234 }
235
236 ~cone() {
237 Polyhedron_Free(Cone);
238 Matrix_Free(Rays);
239 }
240
241 Polyhedron *poly() {
242 if (!Cone) {
243 Matrix *M = Matrix_Alloc(Rays->NbRows+1, Rays->NbColumns+1);
244 for (int i = 0; i < Rays->NbRows; ++i) {
245 Vector_Copy(Rays->p[i], M->p[i]+1, Rays->NbColumns);
246 value_set_si(M->p[i][0], 1);
247 }
248 Vector_Set(M->p[Rays->NbRows]+1, 0, Rays->NbColumns-1);
249 value_set_si(M->p[Rays->NbRows][0], 1);
250 value_set_si(M->p[Rays->NbRows][Rays->NbColumns], 1);
251 Cone = Rays2Polyhedron(M, M->NbRows+1);
252 assert(Cone->NbConstraints == Cone->NbRays);
253 Matrix_Free(M);
254 }
255 return Cone;
256 }
257
258 ZZ det;
259 Polyhedron *Cone;
260 Matrix *Rays;
261 };
262
263 class dpoly {
264 public:
265 vec_ZZ coeff;
266 dpoly(int d, ZZ& degree, int offset = 0) {
267 coeff.SetLength(d+1);
268
269 int min = d + offset;
270 if (degree < ZZ(INIT_VAL, min))
271 min = to_int(degree);
272
273 ZZ c = ZZ(INIT_VAL, 1);
274 if (!offset)
275 coeff[0] = c;
276 for (int i = 1; i <= min; ++i) {
277 c *= (degree -i + 1);
278 c /= i;
279 coeff[i-offset] = c;
280 }
281 }
282 void operator *= (dpoly& f) {
283 assert(coeff.length() == f.coeff.length());
284 vec_ZZ old = coeff;
285 coeff = f.coeff[0] * coeff;
286 for (int i = 1; i < coeff.length(); ++i)
287 for (int j = 0; i+j < coeff.length(); ++j)
288 coeff[i+j] += f.coeff[i] * old[j];
289 }
290 void div(dpoly& d, mpq_t count, ZZ& sign) {
291 int len = coeff.length();
292 Value tmp;
293 value_init(tmp);
294 mpq_t* c = new mpq_t[coeff.length()];
295 mpq_t qtmp;
296 mpq_init(qtmp);
297 for (int i = 0; i < len; ++i) {
298 mpq_init(c[i]);
299 zz2value(coeff[i], tmp);
300 mpq_set_z(c[i], tmp);
301
302 for (int j = 1; j <= i; ++j) {
303 zz2value(d.coeff[j], tmp);
304 mpq_set_z(qtmp, tmp);
305 mpq_mul(qtmp, qtmp, c[i-j]);
306 mpq_sub(c[i], c[i], qtmp);
307 }
308
309 zz2value(d.coeff[0], tmp);
310 mpq_set_z(qtmp, tmp);
311 mpq_div(c[i], c[i], qtmp);
312 }
313 if (sign == -1)
314 mpq_sub(count, count, c[len-1]);
315 else
316 mpq_add(count, count, c[len-1]);
317
318 value_clear(tmp);
319 mpq_clear(qtmp);
320 for (int i = 0; i < len; ++i)
321 mpq_clear(c[i]);
322 delete [] c;
323 }
324 };
325
326 class dpoly_n {
327 public:
328 Matrix *coeff;
329 ~dpoly_n() {
330 Matrix_Free(coeff);
331 }
332 dpoly_n(int d, ZZ& degree_0, ZZ& degree_1, int offset = 0) {
333 Value d0, d1;
334 value_init(d0);
335 value_init(d1);
336 zz2value(degree_0, d0);
337 zz2value(degree_1, d1);
338 coeff = Matrix_Alloc(d+1, d+1+1);
339 value_set_si(coeff->p[0][0], 1);
340 value_set_si(coeff->p[0][d+1], 1);
341 for (int i = 1; i <= d; ++i) {
342 value_multiply(coeff->p[i][0], coeff->p[i-1][0], d0);
343 Vector_Combine(coeff->p[i-1], coeff->p[i-1]+1, coeff->p[i]+1,
344 d1, d0, i);
345 value_set_si(coeff->p[i][d+1], i);
346 value_multiply(coeff->p[i][d+1], coeff->p[i][d+1], coeff->p[i-1][d+1]);
347 value_decrement(d0, d0);
348 }
349 value_clear(d0);
350 value_clear(d1);
351 }
352 void div(dpoly& d, Vector *count, ZZ& sign) {
353 int len = coeff->NbRows;
354 Matrix * c = Matrix_Alloc(coeff->NbRows, coeff->NbColumns);
355 Value tmp;
356 value_init(tmp);
357 for (int i = 0; i < len; ++i) {
358 Vector_Copy(coeff->p[i], c->p[i], len+1);
359 for (int j = 1; j <= i; ++j) {
360 zz2value(d.coeff[j], tmp);
361 value_multiply(tmp, tmp, c->p[i][len]);
362 value_oppose(tmp, tmp);
363 Vector_Combine(c->p[i], c->p[i-j], c->p[i],
364 c->p[i-j][len], tmp, len);
365 value_multiply(c->p[i][len], c->p[i][len], c->p[i-j][len]);
366 }
367 zz2value(d.coeff[0], tmp);
368 value_multiply(c->p[i][len], c->p[i][len], tmp);
369 }
370 if (sign == -1) {
371 value_set_si(tmp, -1);
372 Vector_Scale(c->p[len-1], count->p, tmp, len);
373 value_assign(count->p[len], c->p[len-1][len]);
374 } else
375 Vector_Copy(c->p[len-1], count->p, len+1);
376 Vector_Normalize(count->p, len+1);
377 value_clear(tmp);
378 Matrix_Free(c);
379 }
380 };
381
382 /*
383 * Barvinok's Decomposition of a simplicial cone
384 *
385 * Returns two lists of polyhedra
386 */
387 void barvinok_decompose(Polyhedron *C, Polyhedron **ppos, Polyhedron **pneg)
388 {
389 Polyhedron *pos = *ppos, *neg = *pneg;
390 vector<cone *> nonuni;
391 cone * c = new cone(C);
392 ZZ det = c->det;
393 int s = sign(det);
394 assert(det != 0);
395 if (abs(det) > 1) {
396 nonuni.push_back(c);
397 } else {
398 Polyhedron *p = Polyhedron_Copy(c->Cone);
399 p->next = pos;
400 pos = p;
401 delete c;
402 }
403 vec_ZZ lambda;
404 while (!nonuni.empty()) {
405 c = nonuni.back();
406 nonuni.pop_back();
407 Vector* v = c->short_vector(lambda);
408 for (int i = 0; i < c->Rays->NbRows - 1; ++i) {
409 if (lambda[i] == 0)
410 continue;
411 Matrix* M = Matrix_Copy(c->Rays);
412 Vector_Copy(v->p, M->p[i], v->Size);
413 cone * pc = new cone(M);
414 assert (pc->det != 0);
415 if (abs(pc->det) > 1) {
416 assert(abs(pc->det) < abs(c->det));
417 nonuni.push_back(pc);
418 } else {
419 Polyhedron *p = pc->poly();
420 pc->Cone = 0;
421 if (sign(pc->det) == s) {
422 p->next = pos;
423 pos = p;
424 } else {
425 p->next = neg;
426 neg = p;
427 }
428 delete pc;
429 }
430 Matrix_Free(M);
431 }
432 Vector_Free(v);
433 delete c;
434 }
435 *ppos = pos;
436 *pneg = neg;
437 }
438
439 /*
440 * Returns a single list of npos "positive" cones followed by nneg
441 * "negative" cones.
442 * The input cone is freed
443 */
444 void decompose(Polyhedron *cone, Polyhedron **parts, int *npos, int *nneg, unsigned MaxRays)
445 {
446 Polyhedron_Polarize(cone);
447 if (cone->NbRays - 1 != cone->Dimension) {
448 Polyhedron *tmp = cone;
449 cone = triangularize_cone(cone, MaxRays);
450 Polyhedron_Free(tmp);
451 }
452 Polyhedron *polpos = NULL, *polneg = NULL;
453 *npos = 0; *nneg = 0;
454 for (Polyhedron *Polar = cone; Polar; Polar = Polar->next)
455 barvinok_decompose(Polar, &polpos, &polneg);
456
457 Polyhedron *last;
458 for (Polyhedron *i = polpos; i; i = i->next) {
459 Polyhedron_Polarize(i);
460 ++*npos;
461 last = i;
462 }
463 for (Polyhedron *i = polneg; i; i = i->next) {
464 Polyhedron_Polarize(i);
465 ++*nneg;
466 }
467 if (last) {
468 last->next = polneg;
469 *parts = polpos;
470 } else
471 *parts = polneg;
472 Domain_Free(cone);
473 }
474
475 const int MAX_TRY=10;
476 /*
477 * Searches for a vector that is not othogonal to any
478 * of the rays in rays.
479 */
480 static void nonorthog(mat_ZZ& rays, vec_ZZ& lambda)
481 {
482 int dim = rays.NumCols();
483 bool found = false;
484 lambda.SetLength(dim);
485 for (int i = 2; !found && i <= 50*dim; i+=4) {
486 for (int j = 0; j < MAX_TRY; ++j) {
487 for (int k = 0; k < dim; ++k) {
488 int r = random_int(i)+2;
489 int v = (2*(r%2)-1) * (r >> 1);
490 lambda[k] = v;
491 }
492 int k = 0;
493 for (; k < rays.NumRows(); ++k)
494 if (lambda * rays[k] == 0)
495 break;
496 if (k == rays.NumRows()) {
497 found = true;
498 break;
499 }
500 }
501 }
502 assert(found);
503 }
504
505 static void add_rays(mat_ZZ& rays, Polyhedron *i, int *r)
506 {
507 unsigned dim = i->Dimension;
508 for (int k = 0; k < i->NbRays; ++k) {
509 if (!value_zero_p(i->Ray[k][dim+1]))
510 continue;
511 values2zz(i->Ray[k]+1, rays[(*r)++], dim);
512 }
513 }
514
515 void lattice_point(Value* values, Polyhedron *i, vec_ZZ& lambda, ZZ& num)
516 {
517 vec_ZZ vertex;
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 num = vertex * lambda;
542 }
543
544 static evalue *term(int param, ZZ& c, Value *den = NULL)
545 {
546 evalue *EP = new evalue();
547 value_init(EP->d);
548 value_set_si(EP->d,0);
549 EP->x.p = new_enode(polynomial, 2, param + 1);
550 evalue_set_si(&EP->x.p->arr[0], 0, 1);
551 value_init(EP->x.p->arr[1].x.n);
552 if (den == NULL)
553 value_set_si(EP->x.p->arr[1].d, 1);
554 else
555 value_assign(EP->x.p->arr[1].d, *den);
556 zz2value(c, EP->x.p->arr[1].x.n);
557 return EP;
558 }
559
560 static void vertex_period(
561 Polyhedron *i, vec_ZZ& lambda, Matrix *T,
562 Value lcm, int p, Vector *val,
563 evalue *E, evalue* ev,
564 ZZ& offset)
565 {
566 unsigned nparam = T->NbRows - 1;
567 unsigned dim = i->Dimension;
568 Value tmp;
569 ZZ nump;
570
571 if (p == nparam) {
572 ZZ num, l;
573 Vector * values = Vector_Alloc(dim + 1);
574 Vector_Matrix_Product(val->p, T, values->p);
575 value_assign(values->p[dim], lcm);
576 lattice_point(values->p, i, lambda, num);
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(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(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 void barvinok_count(Polyhedron *P, Value* result, unsigned NbMaxCons)
1184 {
1185 Polyhedron ** vcone;
1186 vec_ZZ sign;
1187 int ncone = 0;
1188 sign.SetLength(ncone);
1189 unsigned dim;
1190 int allocated = 0;
1191 Value factor;
1192 Polyhedron *Q;
1193 int r = 0;
1194
1195 if (emptyQ(P)) {
1196 value_set_si(*result, 0);
1197 return;
1198 }
1199 if (P->NbBid == 0)
1200 for (; r < P->NbRays; ++r)
1201 if (value_zero_p(P->Ray[r][P->Dimension+1]))
1202 break;
1203 if (P->NbBid !=0 || r < P->NbRays) {
1204 value_set_si(*result, -1);
1205 return;
1206 }
1207 if (P->NbEq != 0) {
1208 P = remove_equalities(P);
1209 if (emptyQ(P)) {
1210 Polyhedron_Free(P);
1211 value_set_si(*result, 0);
1212 return;
1213 }
1214 allocated = 1;
1215 }
1216 value_init(factor);
1217 value_set_si(factor, 1);
1218 Q = Polyhedron_Reduce(P, &factor);
1219 if (Q) {
1220 if (allocated)
1221 Polyhedron_Free(P);
1222 P = Q;
1223 allocated = 1;
1224 }
1225 if (P->Dimension == 0) {
1226 value_assign(*result, factor);
1227 if (allocated)
1228 Polyhedron_Free(P);
1229 value_clear(factor);
1230 return;
1231 }
1232
1233 dim = P->Dimension;
1234 vcone = new (Polyhedron *)[P->NbRays];
1235
1236 for (int j = 0; j < P->NbRays; ++j) {
1237 int npos, nneg;
1238 Polyhedron *C = supporting_cone(P, j);
1239 decompose(C, &vcone[j], &npos, &nneg, NbMaxCons);
1240 ncone += npos + nneg;
1241 sign.SetLength(ncone);
1242 for (int k = 0; k < npos; ++k)
1243 sign[ncone-nneg-k-1] = 1;
1244 for (int k = 0; k < nneg; ++k)
1245 sign[ncone-k-1] = -1;
1246 }
1247
1248 mat_ZZ rays;
1249 rays.SetDims(ncone * dim, dim);
1250 r = 0;
1251 for (int j = 0; j < P->NbRays; ++j) {
1252 for (Polyhedron *i = vcone[j]; i; i = i->next) {
1253 assert(i->NbRays-1 == dim);
1254 add_rays(rays, i, &r);
1255 }
1256 }
1257 vec_ZZ lambda;
1258 nonorthog(rays, lambda);
1259
1260 vec_ZZ num;
1261 mat_ZZ den;
1262 num.SetLength(ncone);
1263 den.SetDims(ncone,dim);
1264
1265 int f = 0;
1266 for (int j = 0; j < P->NbRays; ++j) {
1267 for (Polyhedron *i = vcone[j]; i; i = i->next) {
1268 lattice_point(P->Ray[j]+1, i, lambda, num[f]);
1269 normalize(i, lambda, sign[f], num[f], den[f]);
1270 ++f;
1271 }
1272 }
1273 ZZ min = num[0];
1274 for (int j = 1; j < num.length(); ++j)
1275 if (num[j] < min)
1276 min = num[j];
1277 for (int j = 0; j < num.length(); ++j)
1278 num[j] -= min;
1279
1280 f = 0;
1281 mpq_t count;
1282 mpq_init(count);
1283 for (int j = 0; j < P->NbRays; ++j) {
1284 for (Polyhedron *i = vcone[j]; i; i = i->next) {
1285 dpoly d(dim, num[f]);
1286 dpoly n(dim, den[f][0], 1);
1287 for (int k = 1; k < dim; ++k) {
1288 dpoly fact(dim, den[f][k], 1);
1289 n *= fact;
1290 }
1291 d.div(n, count, sign[f]);
1292 ++f;
1293 }
1294 }
1295 assert(value_one_p(&count[0]._mp_den));
1296 value_multiply(*result, &count[0]._mp_num, factor);
1297 mpq_clear(count);
1298
1299 for (int j = 0; j < P->NbRays; ++j)
1300 Domain_Free(vcone[j]);
1301
1302 delete [] vcone;
1303
1304 if (allocated)
1305 Polyhedron_Free(P);
1306 value_clear(factor);
1307 }
1308
1309 static void uni_polynom(int param, Vector *c, evalue *EP)
1310 {
1311 unsigned dim = c->Size-2;
1312 value_init(EP->d);
1313 value_set_si(EP->d,0);
1314 EP->x.p = new_enode(polynomial, dim+1, param+1);
1315 for (int j = 0; j <= dim; ++j)
1316 evalue_set(&EP->x.p->arr[j], c->p[j], c->p[dim+1]);
1317 }
1318
1319 static void multi_polynom(Vector *c, evalue* X, evalue *EP)
1320 {
1321 unsigned dim = c->Size-2;
1322 evalue EC;
1323
1324 value_init(EC.d);
1325 evalue_set(&EC, c->p[dim], c->p[dim+1]);
1326
1327 value_init(EP->d);
1328 evalue_set(EP, c->p[dim], c->p[dim+1]);
1329
1330 for (int i = dim-1; i >= 0; --i) {
1331 emul(X, EP);
1332 value_assign(EC.x.n, c->p[i]);
1333 eadd(&EC, EP);
1334 }
1335 free_evalue_refs(&EC);
1336 }
1337
1338 Polyhedron *unfringe (Polyhedron *P, unsigned MaxRays)
1339 {
1340 int len = P->Dimension+2;
1341 Polyhedron *T, *R = P;
1342 Value g;
1343 value_init(g);
1344 Vector *row = Vector_Alloc(len);
1345 value_set_si(row->p[0], 1);
1346
1347 R = DomainConstraintSimplify(Polyhedron_Copy(P), MaxRays);
1348
1349 Matrix *M = Matrix_Alloc(2, len-1);
1350 value_set_si(M->p[1][len-2], 1);
1351 for (int v = 0; v < P->Dimension; ++v) {
1352 value_set_si(M->p[0][v], 1);
1353 Polyhedron *I = Polyhedron_Image(P, M, 2+1);
1354 value_set_si(M->p[0][v], 0);
1355 for (int r = 0; r < I->NbConstraints; ++r) {
1356 if (value_zero_p(I->Constraint[r][0]))
1357 continue;
1358 if (value_zero_p(I->Constraint[r][1]))
1359 continue;
1360 if (value_one_p(I->Constraint[r][1]))
1361 continue;
1362 if (value_mone_p(I->Constraint[r][1]))
1363 continue;
1364 value_absolute(g, I->Constraint[r][1]);
1365 Vector_Set(row->p+1, 0, len-2);
1366 value_division(row->p[1+v], I->Constraint[r][1], g);
1367 mpz_fdiv_q(row->p[len-1], I->Constraint[r][2], g);
1368 T = R;
1369 R = AddConstraints(row->p, 1, R, MaxRays);
1370 if (T != P)
1371 Polyhedron_Free(T);
1372 }
1373 Polyhedron_Free(I);
1374 }
1375 Matrix_Free(M);
1376 Vector_Free(row);
1377 value_clear(g);
1378 return R;
1379 }
1380
1381 static Polyhedron *reduce_domain(Polyhedron *D, Matrix *CT, Polyhedron *CEq,
1382 Polyhedron **fVD, int nd, unsigned MaxRays)
1383 {
1384 assert(CEq);
1385
1386 Polyhedron *Dt;
1387 Dt = CT ? DomainPreimage(D, CT, MaxRays) : D;
1388 Polyhedron *rVD = DomainIntersection(Dt, CEq, MaxRays);
1389
1390 /* if rVD is empty or too small in geometric dimension */
1391 if(!rVD || emptyQ(rVD) ||
1392 (rVD->Dimension-rVD->NbEq < Dt->Dimension-Dt->NbEq-CEq->NbEq)) {
1393 if(rVD)
1394 Domain_Free(rVD);
1395 if (CT)
1396 Domain_Free(Dt);
1397 return 0; /* empty validity domain */
1398 }
1399
1400 if (CT)
1401 Domain_Free(Dt);
1402
1403 fVD[nd] = Domain_Copy(rVD);
1404 for (int i = 0 ; i < nd; ++i) {
1405 Polyhedron *I = DomainIntersection(fVD[nd], fVD[i], MaxRays);
1406 if (emptyQ(I)) {
1407 Domain_Free(I);
1408 continue;
1409 }
1410 Polyhedron *F = DomainSimplify(I, fVD[nd], MaxRays);
1411 if (F->NbEq == 1) {
1412 Polyhedron *T = rVD;
1413 rVD = DomainDifference(rVD, F, MaxRays);
1414 Domain_Free(T);
1415 }
1416 Domain_Free(F);
1417 Domain_Free(I);
1418 }
1419
1420 rVD = DomainConstraintSimplify(rVD, MaxRays);
1421 if (emptyQ(rVD)) {
1422 Domain_Free(fVD[nd]);
1423 Domain_Free(rVD);
1424 return 0;
1425 }
1426
1427 Value c;
1428 value_init(c);
1429 barvinok_count(rVD, &c, MaxRays);
1430 if (value_zero_p(c)) {
1431 Domain_Free(rVD);
1432 rVD = 0;
1433 }
1434 value_clear(c);
1435
1436 return rVD;
1437 }
1438
1439 evalue* barvinok_enumerate_ev(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
1440 {
1441 //P = unfringe(P, MaxRays);
1442 Polyhedron *CEq = NULL, *rVD, *pVD, *CA;
1443 Matrix *CT = NULL;
1444 Param_Polyhedron *PP = NULL;
1445 Param_Domain *D, *next;
1446 Param_Vertices *V;
1447 int r = 0;
1448 unsigned nparam = C->Dimension;
1449 evalue *eres;
1450 ALLOC(eres);
1451 value_init(eres->d);
1452 value_set_si(eres->d, 0);
1453
1454 evalue factor;
1455 value_init(factor.d);
1456 evalue_set_si(&factor, 1, 1);
1457
1458 CA = align_context(C, P->Dimension, MaxRays);
1459 P = DomainIntersection(P, CA, MaxRays);
1460 Polyhedron_Free(CA);
1461
1462 if (C->Dimension == 0 || emptyQ(P)) {
1463 constant:
1464 eres->x.p = new_enode(partition, 2, C->Dimension);
1465 EVALUE_SET_DOMAIN(eres->x.p->arr[0],
1466 DomainConstraintSimplify(CEq ? CEq : Polyhedron_Copy(C), MaxRays));
1467 value_set_si(eres->x.p->arr[1].d, 1);
1468 value_init(eres->x.p->arr[1].x.n);
1469 if (emptyQ(P))
1470 value_set_si(eres->x.p->arr[1].x.n, 0);
1471 else
1472 barvinok_count(P, &eres->x.p->arr[1].x.n, MaxRays);
1473 out:
1474 emul(&factor, eres);
1475 reduce_evalue(eres);
1476 free_evalue_refs(&factor);
1477 Polyhedron_Free(P);
1478 if (CT)
1479 Matrix_Free(CT);
1480 if (PP)
1481 Param_Polyhedron_Free(PP);
1482
1483 return eres;
1484 }
1485 for (r = 0; r < P->NbRays; ++r)
1486 if (value_zero_p(P->Ray[r][0]) ||
1487 value_zero_p(P->Ray[r][P->Dimension+1])) {
1488 int i;
1489 for (i = P->Dimension - nparam; i < P->Dimension; ++i)
1490 if (value_notzero_p(P->Ray[r][i+1]))
1491 break;
1492 if (i >= P->Dimension)
1493 break;
1494 }
1495 if (r < P->NbRays)
1496 goto constant;
1497
1498 if (P->NbEq != 0) {
1499 Matrix *f;
1500 P = remove_equalities_p(P, P->Dimension-nparam, &f);
1501 mask(f, &factor);
1502 Matrix_Free(f);
1503 }
1504 if (P->Dimension == nparam) {
1505 CEq = P;
1506 P = Universe_Polyhedron(0);
1507 goto constant;
1508 }
1509
1510 Polyhedron *Q = ParamPolyhedron_Reduce(P, P->Dimension-nparam, &factor);
1511 if (Q) {
1512 Polyhedron_Free(P);
1513 if (Q->Dimension == nparam) {
1514 CEq = Q;
1515 P = Universe_Polyhedron(0);
1516 goto constant;
1517 }
1518 P = Q;
1519 }
1520 Polyhedron *oldP = P;
1521 PP = Polyhedron2Param_SimplifiedDomain(&P,C,MaxRays,&CEq,&CT);
1522 if (P != oldP)
1523 Polyhedron_Free(oldP);
1524
1525 if (isIdentity(CT)) {
1526 Matrix_Free(CT);
1527 CT = NULL;
1528 } else {
1529 assert(CT->NbRows != CT->NbColumns);
1530 if (CT->NbRows == 1) // no more parameters
1531 goto constant;
1532 nparam = CT->NbRows - 1;
1533 }
1534
1535 unsigned dim = P->Dimension - nparam;
1536 Polyhedron ** vcone = new (Polyhedron *)[PP->nbV];
1537 int * npos = new int[PP->nbV];
1538 int * nneg = new int[PP->nbV];
1539 vec_ZZ sign;
1540
1541 int i;
1542 for (i = 0, V = PP->V; V; ++i, V = V->next) {
1543 Polyhedron *C = supporting_cone_p(P, V);
1544 decompose(C, &vcone[i], &npos[i], &nneg[i], MaxRays);
1545 }
1546
1547 Vector *c = Vector_Alloc(dim+2);
1548
1549 int nd;
1550 for (nd = 0, D=PP->D; D; ++nd, D=D->next);
1551 struct section { Polyhedron *D; evalue E; };
1552 section *s = new section[nd];
1553 Polyhedron **fVD = new (Polyhedron*)[nd];
1554
1555 for(nd = 0, D=PP->D; D; D=next) {
1556 next = D->next;
1557
1558 Polyhedron *rVD = reduce_domain(D->Domain, CT, CEq,
1559 fVD, nd, MaxRays);
1560 if (!rVD)
1561 continue;
1562
1563 pVD = CT ? DomainImage(rVD,CT,MaxRays) : rVD;
1564
1565 int ncone = 0;
1566 sign.SetLength(ncone);
1567 FORALL_PVertex_in_ParamPolyhedron(V,D,PP) // _i is internal counter
1568 ncone += npos[_i] + nneg[_i];
1569 sign.SetLength(ncone);
1570 for (int k = 0; k < npos[_i]; ++k)
1571 sign[ncone-nneg[_i]-k-1] = 1;
1572 for (int k = 0; k < nneg[_i]; ++k)
1573 sign[ncone-k-1] = -1;
1574 END_FORALL_PVertex_in_ParamPolyhedron;
1575
1576 mat_ZZ rays;
1577 rays.SetDims(ncone * dim, dim);
1578 r = 0;
1579 FORALL_PVertex_in_ParamPolyhedron(V,D,PP) // _i is internal counter
1580 for (Polyhedron *i = vcone[_i]; i; i = i->next) {
1581 assert(i->NbRays-1 == dim);
1582 add_rays(rays, i, &r);
1583 }
1584 END_FORALL_PVertex_in_ParamPolyhedron;
1585 vec_ZZ lambda;
1586 nonorthog(rays, lambda);
1587
1588 mat_ZZ den;
1589 den.SetDims(ncone,dim);
1590 term_info *num = new term_info[ncone];
1591
1592 int f = 0;
1593 FORALL_PVertex_in_ParamPolyhedron(V,D,PP)
1594 for (Polyhedron *i = vcone[_i]; i; i = i->next) {
1595 lattice_point(V, i, lambda, &num[f], pVD);
1596 normalize(i, lambda, sign[f], num[f].constant, den[f]);
1597 ++f;
1598 }
1599 END_FORALL_PVertex_in_ParamPolyhedron;
1600 ZZ min = num[0].constant;
1601 for (int j = 1; j < ncone; ++j)
1602 if (num[j].constant < min)
1603 min = num[j].constant;
1604 for (int j = 0; j < ncone; ++j)
1605 num[j].constant -= min;
1606 f = 0;
1607 value_init(s[nd].E.d);
1608 evalue_set_si(&s[nd].E, 0, 1);
1609 mpq_t count;
1610 mpq_init(count);
1611 FORALL_PVertex_in_ParamPolyhedron(V,D,PP)
1612 for (Polyhedron *i = vcone[_i]; i; i = i->next) {
1613 dpoly n(dim, den[f][0], 1);
1614 for (int k = 1; k < dim; ++k) {
1615 dpoly fact(dim, den[f][k], 1);
1616 n *= fact;
1617 }
1618 if (num[f].E != NULL) {
1619 ZZ one(INIT_VAL, 1);
1620 dpoly_n d(dim, num[f].constant, one);
1621 d.div(n, c, sign[f]);
1622 evalue EV;
1623 multi_polynom(c, num[f].E, &EV);
1624 eadd(&EV , &s[nd].E);
1625 free_evalue_refs(&EV);
1626 free_evalue_refs(num[f].E);
1627 delete num[f].E;
1628 } else if (num[f].pos != -1) {
1629 dpoly_n d(dim, num[f].constant, num[f].coeff);
1630 d.div(n, c, sign[f]);
1631 evalue EV;
1632 uni_polynom(num[f].pos, c, &EV);
1633 eadd(&EV , &s[nd].E);
1634 free_evalue_refs(&EV);
1635 } else {
1636 mpq_set_si(count, 0, 1);
1637 dpoly d(dim, num[f].constant);
1638 d.div(n, count, sign[f]);
1639 evalue EV;
1640 value_init(EV.d);
1641 evalue_set(&EV, &count[0]._mp_num, &count[0]._mp_den);
1642 eadd(&EV , &s[nd].E);
1643 free_evalue_refs(&EV);
1644 }
1645 ++f;
1646 }
1647 END_FORALL_PVertex_in_ParamPolyhedron;
1648
1649 mpq_clear(count);
1650 delete [] num;
1651
1652 if (CT)
1653 addeliminatedparams_evalue(&s[nd].E, CT);
1654 s[nd].D = rVD;
1655 ++nd;
1656 if (rVD != pVD)
1657 Domain_Free(pVD);
1658 }
1659
1660 if (nd == 0)
1661 evalue_set_si(eres, 0, 1);
1662 else {
1663 eres->x.p = new_enode(partition, 2*nd, C->Dimension);
1664 for (int j = 0; j < nd; ++j) {
1665 EVALUE_SET_DOMAIN(eres->x.p->arr[2*j], s[j].D);
1666 value_clear(eres->x.p->arr[2*j+1].d);
1667 eres->x.p->arr[2*j+1] = s[j].E;
1668 Domain_Free(fVD[j]);
1669 }
1670 }
1671 delete [] s;
1672 delete [] fVD;
1673
1674 Vector_Free(c);
1675
1676 for (int j = 0; j < PP->nbV; ++j)
1677 Domain_Free(vcone[j]);
1678 delete [] vcone;
1679 delete [] npos;
1680 delete [] nneg;
1681
1682 if (CEq)
1683 Polyhedron_Free(CEq);
1684
1685 goto out;
1686 }
1687
1688 Enumeration* barvinok_enumerate(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
1689 {
1690 evalue *EP = barvinok_enumerate_ev(P, C, MaxRays);
1691
1692 return partition2enumeration(EP);
1693 }
1694
1695 static void SwapColumns(Value **V, int n, int i, int j)
1696 {
1697 for (int r = 0; r < n; ++r)
1698 value_swap(V[r][i], V[r][j]);
1699 }
1700
1701 static void SwapColumns(Polyhedron *P, int i, int j)
1702 {
1703 SwapColumns(P->Constraint, P->NbConstraints, i, j);
1704 SwapColumns(P->Ray, P->NbRays, i, j);
1705 }
1706
1707 static void negative_test_constraint(Value *l, Value *u, Value *c, int pos,
1708 int len, Value *v)
1709 {
1710 value_oppose(*v, u[pos+1]);
1711 Vector_Combine(l+1, u+1, c+1, *v, l[pos+1], len-1);
1712 value_multiply(*v, *v, l[pos+1]);
1713 value_substract(c[len-1], c[len-1], *v);
1714 value_set_si(*v, -1);
1715 Vector_Scale(c+1, c+1, *v, len-1);
1716 value_decrement(c[len-1], c[len-1]);
1717 ConstraintSimplify(c, c, len, v);
1718 }
1719
1720 static void oppose_constraint(Value *c, int len, Value *v)
1721 {
1722 value_set_si(*v, -1);
1723 Vector_Scale(c+1, c+1, *v, len-1);
1724 value_decrement(c[len-1], c[len-1]);
1725 }
1726
1727 static bool SplitOnConstraint(Polyhedron *P, int i, int l, int u,
1728 int nvar, int len, int exist, int MaxRays,
1729 Vector *row, Value& f, bool independent,
1730 Polyhedron **pos, Polyhedron **neg)
1731 {
1732 negative_test_constraint(P->Constraint[l], P->Constraint[u],
1733 row->p, nvar+i, len, &f);
1734 *neg = AddConstraints(row->p, 1, P, MaxRays);
1735
1736 /* We found an independent, but useless constraint
1737 * Maybe we should detect this earlier and not
1738 * mark the variable as INDEPENDENT
1739 */
1740 if (emptyQ((*neg))) {
1741 Polyhedron_Free(*neg);
1742 return false;
1743 }
1744
1745 oppose_constraint(row->p, len, &f);
1746 *pos = AddConstraints(row->p, 1, P, MaxRays);
1747
1748 if (emptyQ((*pos))) {
1749 Polyhedron_Free(*neg);
1750 Polyhedron_Free(*pos);
1751 return false;
1752 }
1753
1754 return true;
1755 }
1756
1757 /*
1758 * unimodularly transform P such that constraint r is transformed
1759 * into a constraint that involves only a single (the first)
1760 * existential variable
1761 *
1762 */
1763 static Polyhedron *rotate_along(Polyhedron *P, int r, int nvar, int exist,
1764 unsigned MaxRays)
1765 {
1766 Value g;
1767 value_init(g);
1768
1769 Vector *row = Vector_Alloc(exist);
1770 Vector_Copy(P->Constraint[r]+1+nvar, row->p, exist);
1771 Vector_Gcd(row->p, exist, &g);
1772 if (value_notone_p(g))
1773 Vector_AntiScale(row->p, row->p, g, exist);
1774 value_clear(g);
1775
1776 Matrix *M = unimodular_complete(row);
1777 Matrix *M2 = Matrix_Alloc(P->Dimension+1, P->Dimension+1);
1778 for (r = 0; r < nvar; ++r)
1779 value_set_si(M2->p[r][r], 1);
1780 for ( ; r < nvar+exist; ++r)
1781 Vector_Copy(M->p[r-nvar], M2->p[r]+nvar, exist);
1782 for ( ; r < P->Dimension+1; ++r)
1783 value_set_si(M2->p[r][r], 1);
1784 Polyhedron *T = Polyhedron_Image(P, M2, MaxRays);
1785
1786 Matrix_Free(M2);
1787 Matrix_Free(M);
1788 Vector_Free(row);
1789
1790 return T;
1791 }
1792
1793 static bool SplitOnVar(Polyhedron *P, int i,
1794 int nvar, int len, int exist, int MaxRays,
1795 Vector *row, Value& f, bool independent,
1796 Polyhedron **pos, Polyhedron **neg)
1797 {
1798 int j;
1799
1800 for (int l = P->NbEq; l < P->NbConstraints; ++l) {
1801 if (value_negz_p(P->Constraint[l][nvar+i+1]))
1802 continue;
1803
1804 if (independent) {
1805 for (j = 0; j < exist; ++j)
1806 if (j != i && value_notzero_p(P->Constraint[l][nvar+j+1]))
1807 break;
1808 if (j < exist)
1809 continue;
1810 }
1811
1812 for (int u = P->NbEq; u < P->NbConstraints; ++u) {
1813 if (value_posz_p(P->Constraint[u][nvar+i+1]))
1814 continue;
1815
1816 if (independent) {
1817 for (j = 0; j < exist; ++j)
1818 if (j != i && value_notzero_p(P->Constraint[u][nvar+j+1]))
1819 break;
1820 if (j < exist)
1821 continue;
1822 }
1823
1824 if (SplitOnConstraint(P, i, l, u,
1825 nvar, len, exist, MaxRays,
1826 row, f, independent,
1827 pos, neg)) {
1828 if (independent) {
1829 if (i != 0)
1830 SwapColumns(*neg, nvar+1, nvar+1+i);
1831 }
1832 return true;
1833 }
1834 }
1835 }
1836
1837 return false;
1838 }
1839
1840 static bool double_bound_pair(Polyhedron *P, int nvar, int exist,
1841 int i, int l1, int l2,
1842 Polyhedron **pos, Polyhedron **neg)
1843 {
1844 Value f;
1845 value_init(f);
1846 Vector *row = Vector_Alloc(P->Dimension+2);
1847 value_set_si(row->p[0], 1);
1848 value_oppose(f, P->Constraint[l1][nvar+i+1]);
1849 Vector_Combine(P->Constraint[l1]+1, P->Constraint[l2]+1,
1850 row->p+1,
1851 P->Constraint[l2][nvar+i+1], f,
1852 P->Dimension+1);
1853 ConstraintSimplify(row->p, row->p, P->Dimension+2, &f);
1854 *pos = AddConstraints(row->p, 1, P, 0);
1855 value_set_si(f, -1);
1856 Vector_Scale(row->p+1, row->p+1, f, P->Dimension+1);
1857 value_decrement(row->p[P->Dimension+1], row->p[P->Dimension+1]);
1858 *neg = AddConstraints(row->p, 1, P, 0);
1859 Vector_Free(row);
1860 value_clear(f);
1861
1862 return !emptyQ((*pos)) && !emptyQ((*neg));
1863 }
1864
1865 static bool double_bound(Polyhedron *P, int nvar, int exist,
1866 Polyhedron **pos, Polyhedron **neg)
1867 {
1868 for (int i = 0; i < exist; ++i) {
1869 int l1, l2;
1870 for (l1 = P->NbEq; l1 < P->NbConstraints; ++l1) {
1871 if (value_negz_p(P->Constraint[l1][nvar+i+1]))
1872 continue;
1873 for (l2 = l1 + 1; l2 < P->NbConstraints; ++l2) {
1874 if (value_negz_p(P->Constraint[l2][nvar+i+1]))
1875 continue;
1876 if (double_bound_pair(P, nvar, exist, i, l1, l2, pos, neg))
1877 return true;
1878 }
1879 }
1880 for (l1 = P->NbEq; l1 < P->NbConstraints; ++l1) {
1881 if (value_posz_p(P->Constraint[l1][nvar+i+1]))
1882 continue;
1883 if (l1 < P->NbConstraints)
1884 for (l2 = l1 + 1; l2 < P->NbConstraints; ++l2) {
1885 if (value_posz_p(P->Constraint[l2][nvar+i+1]))
1886 continue;
1887 if (double_bound_pair(P, nvar, exist, i, l1, l2, pos, neg))
1888 return true;
1889 }
1890 }
1891 return false;
1892 }
1893 return false;
1894 }
1895
1896 enum constraint {
1897 ALL_POS = 1 << 0,
1898 ONE_NEG = 1 << 1,
1899 INDEPENDENT = 1 << 2,
1900 };
1901
1902 static evalue* enumerate_or(Polyhedron *D,
1903 unsigned exist, unsigned nparam, unsigned MaxRays)
1904 {
1905 #ifdef DEBUG_ER
1906 fprintf(stderr, "\nER: Or\n");
1907 #endif /* DEBUG_ER */
1908
1909 Polyhedron *N = D->next;
1910 D->next = 0;
1911 evalue *EP =
1912 barvinok_enumerate_e(D, exist, nparam, MaxRays);
1913 Polyhedron_Free(D);
1914
1915 for (D = N; D; D = N) {
1916 N = D->next;
1917 D->next = 0;
1918
1919 evalue *EN =
1920 barvinok_enumerate_e(D, exist, nparam, MaxRays);
1921
1922 eor(EN, EP);
1923 free_evalue_refs(EN);
1924 free(EN);
1925 Polyhedron_Free(D);
1926 }
1927
1928 reduce_evalue(EP);
1929
1930 return EP;
1931 }
1932
1933 static evalue* enumerate_sum(Polyhedron *P,
1934 unsigned exist, unsigned nparam, unsigned MaxRays)
1935 {
1936 int nvar = P->Dimension - exist - nparam;
1937 int toswap = nvar < exist ? nvar : exist;
1938 for (int i = 0; i < toswap; ++i)
1939 SwapColumns(P, 1 + i, nvar+exist - i);
1940 nparam += nvar;
1941
1942 #ifdef DEBUG_ER
1943 fprintf(stderr, "\nER: Sum\n");
1944 #endif /* DEBUG_ER */
1945
1946 evalue *EP = barvinok_enumerate_e(P, exist, nparam, MaxRays);
1947
1948 for (int i = 0; i < /* nvar */ nparam; ++i) {
1949 Matrix *C = Matrix_Alloc(1, 1 + nparam + 1);
1950 value_set_si(C->p[0][0], 1);
1951 evalue split;
1952 value_init(split.d);
1953 value_set_si(split.d, 0);
1954 split.x.p = new_enode(partition, 4, nparam);
1955 value_set_si(C->p[0][1+i], 1);
1956 Matrix *C2 = Matrix_Copy(C);
1957 EVALUE_SET_DOMAIN(split.x.p->arr[0],
1958 Constraints2Polyhedron(C2, MaxRays));
1959 Matrix_Free(C2);
1960 evalue_set_si(&split.x.p->arr[1], 1, 1);
1961 value_set_si(C->p[0][1+i], -1);
1962 value_set_si(C->p[0][1+nparam], -1);
1963 EVALUE_SET_DOMAIN(split.x.p->arr[2],
1964 Constraints2Polyhedron(C, MaxRays));
1965 evalue_set_si(&split.x.p->arr[3], 1, 1);
1966 emul(&split, EP);
1967 free_evalue_refs(&split);
1968 Matrix_Free(C);
1969 }
1970 reduce_evalue(EP);
1971 evalue_range_reduction(EP);
1972
1973 evalue_frac2floor(EP);
1974
1975 evalue *sum = esum(EP, nvar);
1976
1977 free_evalue_refs(EP);
1978 free(EP);
1979 EP = sum;
1980
1981 evalue_range_reduction(EP);
1982
1983 return EP;
1984 }
1985
1986 static evalue* split_sure(Polyhedron *P, Polyhedron *S,
1987 unsigned exist, unsigned nparam, unsigned MaxRays)
1988 {
1989 int nvar = P->Dimension - exist - nparam;
1990
1991 Matrix *M = Matrix_Alloc(exist, S->Dimension+2);
1992 for (int i = 0; i < exist; ++i)
1993 value_set_si(M->p[i][nvar+i+1], 1);
1994 Polyhedron *O = S;
1995 S = DomainAddRays(S, M, MaxRays);
1996 Polyhedron_Free(O);
1997 Polyhedron *F = DomainAddRays(P, M, MaxRays);
1998 Polyhedron *D = DomainDifference(F, S, MaxRays);
1999 O = D;
2000 D = Disjoint_Domain(D, 0, MaxRays);
2001 Polyhedron_Free(F);
2002 Domain_Free(O);
2003 Matrix_Free(M);
2004
2005 M = Matrix_Alloc(P->Dimension+1-exist, P->Dimension+1);
2006 for (int j = 0; j < nvar; ++j)
2007 value_set_si(M->p[j][j], 1);
2008 for (int j = 0; j < nparam+1; ++j)
2009 value_set_si(M->p[nvar+j][nvar+exist+j], 1);
2010 Polyhedron *T = Polyhedron_Image(S, M, MaxRays);
2011 evalue *EP = barvinok_enumerate_e(T, 0, nparam, MaxRays);
2012 Polyhedron_Free(S);
2013 Polyhedron_Free(T);
2014 Matrix_Free(M);
2015
2016 for (Polyhedron *Q = D; Q; Q = Q->next) {
2017 Polyhedron *N = Q->next;
2018 Q->next = 0;
2019 T = DomainIntersection(P, Q, MaxRays);
2020 evalue *E = barvinok_enumerate_e(T, exist, nparam, MaxRays);
2021 eadd(E, EP);
2022 free_evalue_refs(E);
2023 free(E);
2024 Polyhedron_Free(T);
2025 Q->next = N;
2026 }
2027 Domain_Free(D);
2028 return EP;
2029 }
2030
2031 static evalue* enumerate_sure(Polyhedron *P,
2032 unsigned exist, unsigned nparam, unsigned MaxRays)
2033 {
2034 int i;
2035 Polyhedron *S = P;
2036 int nvar = P->Dimension - exist - nparam;
2037 Value lcm;
2038 Value f;
2039 value_init(lcm);
2040 value_init(f);
2041
2042 for (i = 0; i < exist; ++i) {
2043 Matrix *M = Matrix_Alloc(S->NbConstraints, S->Dimension+2);
2044 int c = 0;
2045 value_set_si(lcm, 1);
2046 for (int j = 0; j < S->NbConstraints; ++j) {
2047 if (value_negz_p(S->Constraint[j][1+nvar+i]))
2048 continue;
2049 if (value_one_p(S->Constraint[j][1+nvar+i]))
2050 continue;
2051 value_lcm(lcm, S->Constraint[j][1+nvar+i], &lcm);
2052 }
2053
2054 for (int j = 0; j < S->NbConstraints; ++j) {
2055 if (value_negz_p(S->Constraint[j][1+nvar+i]))
2056 continue;
2057 if (value_one_p(S->Constraint[j][1+nvar+i]))
2058 continue;
2059 value_division(f, lcm, S->Constraint[j][1+nvar+i]);
2060 Vector_Scale(S->Constraint[j], M->p[c], f, S->Dimension+2);
2061 value_substract(M->p[c][S->Dimension+1],
2062 M->p[c][S->Dimension+1],
2063 lcm);
2064 value_increment(M->p[c][S->Dimension+1],
2065 M->p[c][S->Dimension+1]);
2066 ++c;
2067 }
2068 Polyhedron *O = S;
2069 S = AddConstraints(M->p[0], c, S, MaxRays);
2070 if (O != P)
2071 Polyhedron_Free(O);
2072 Matrix_Free(M);
2073 if (emptyQ(S)) {
2074 Polyhedron_Free(S);
2075 value_clear(lcm);
2076 value_clear(f);
2077 return 0;
2078 }
2079 }
2080 value_clear(lcm);
2081 value_clear(f);
2082
2083 #ifdef DEBUG_ER
2084 fprintf(stderr, "\nER: Sure\n");
2085 #endif /* DEBUG_ER */
2086
2087 return split_sure(P, S, exist, nparam, MaxRays);
2088 }
2089
2090 static evalue* enumerate_sure2(Polyhedron *P,
2091 unsigned exist, unsigned nparam, unsigned MaxRays)
2092 {
2093 int nvar = P->Dimension - exist - nparam;
2094 int r;
2095 for (r = 0; r < P->NbRays; ++r)
2096 if (value_one_p(P->Ray[r][0]) &&
2097 value_one_p(P->Ray[r][P->Dimension+1]))
2098 break;
2099
2100 if (r >= P->NbRays)
2101 return 0;
2102
2103 Matrix *M = Matrix_Alloc(nvar + 1 + nparam, P->Dimension+2);
2104 for (int i = 0; i < nvar; ++i)
2105 value_set_si(M->p[i][1+i], 1);
2106 for (int i = 0; i < nparam; ++i)
2107 value_set_si(M->p[i+nvar][1+nvar+exist+i], 1);
2108 Vector_Copy(P->Ray[r]+1+nvar, M->p[nvar+nparam]+1+nvar, exist);
2109 value_set_si(M->p[nvar+nparam][0], 1);
2110 value_set_si(M->p[nvar+nparam][P->Dimension+1], 1);
2111 Polyhedron * F = Rays2Polyhedron(M, MaxRays);
2112 Matrix_Free(M);
2113
2114 Polyhedron *I = DomainIntersection(F, P, MaxRays);
2115 Polyhedron_Free(F);
2116
2117 #ifdef DEBUG_ER
2118 fprintf(stderr, "\nER: Sure2\n");
2119 #endif /* DEBUG_ER */
2120
2121 return split_sure(P, I, exist, nparam, MaxRays);
2122 }
2123
2124 static evalue* enumerate_cyclic(Polyhedron *P,
2125 unsigned exist, unsigned nparam,
2126 evalue * EP, int r, int p, unsigned MaxRays)
2127 {
2128 int nvar = P->Dimension - exist - nparam;
2129
2130 /* If EP in its fractional maps only contains references
2131 * to the remainder parameter with appropriate coefficients
2132 * then we could in principle avoid adding existentially
2133 * quantified variables to the validity domains.
2134 * We'd have to replace the remainder by m { p/m }
2135 * and multiply with an appropriate factor that is one
2136 * only in the appropriate range.
2137 * This last multiplication can be avoided if EP
2138 * has a single validity domain with no (further)
2139 * constraints on the remainder parameter
2140 */
2141
2142 Matrix *CT = Matrix_Alloc(nparam+1, nparam+3);
2143 Matrix *M = Matrix_Alloc(1, 1+nparam+3);
2144 for (int j = 0; j < nparam; ++j)
2145 if (j != p)
2146 value_set_si(CT->p[j][j], 1);
2147 value_set_si(CT->p[p][nparam+1], 1);
2148 value_set_si(CT->p[nparam][nparam+2], 1);
2149 value_set_si(M->p[0][1+p], -1);
2150 value_absolute(M->p[0][1+nparam], P->Ray[0][1+nvar+exist+p]);
2151 value_set_si(M->p[0][1+nparam+1], 1);
2152 Polyhedron *CEq = Constraints2Polyhedron(M, 1);
2153 Matrix_Free(M);
2154 addeliminatedparams_enum(EP, CT, CEq, MaxRays, nparam);
2155 Polyhedron_Free(CEq);
2156 Matrix_Free(CT);
2157
2158 return EP;
2159 }
2160
2161 static void enumerate_vd_add_ray(evalue *EP, Matrix *Rays, unsigned MaxRays)
2162 {
2163 if (value_notzero_p(EP->d))
2164 return;
2165
2166 assert(EP->x.p->type == partition);
2167 assert(EP->x.p->pos == EVALUE_DOMAIN(EP->x.p->arr[0])->Dimension);
2168 for (int i = 0; i < EP->x.p->size/2; ++i) {
2169 Polyhedron *D = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
2170 Polyhedron *N = DomainAddRays(D, Rays, MaxRays);
2171 EVALUE_SET_DOMAIN(EP->x.p->arr[2*i], N);
2172 Domain_Free(D);
2173 }
2174 }
2175
2176 static evalue* enumerate_line(Polyhedron *P,
2177 unsigned exist, unsigned nparam, unsigned MaxRays)
2178 {
2179 if (P->NbBid == 0)
2180 return 0;
2181
2182 #ifdef DEBUG_ER
2183 fprintf(stderr, "\nER: Line\n");
2184 #endif /* DEBUG_ER */
2185
2186 int nvar = P->Dimension - exist - nparam;
2187 int i, j;
2188 for (i = 0; i < nparam; ++i)
2189 if (value_notzero_p(P->Ray[0][1+nvar+exist+i]))
2190 break;
2191 assert(i < nparam);
2192 for (j = i+1; j < nparam; ++j)
2193 if (value_notzero_p(P->Ray[0][1+nvar+exist+i]))
2194 break;
2195 assert(j >= nparam); // for now
2196
2197 Matrix *M = Matrix_Alloc(2, P->Dimension+2);
2198 value_set_si(M->p[0][0], 1);
2199 value_set_si(M->p[0][1+nvar+exist+i], 1);
2200 value_set_si(M->p[1][0], 1);
2201 value_set_si(M->p[1][1+nvar+exist+i], -1);
2202 value_absolute(M->p[1][1+P->Dimension], P->Ray[0][1+nvar+exist+i]);
2203 value_decrement(M->p[1][1+P->Dimension], M->p[1][1+P->Dimension]);
2204 Polyhedron *S = AddConstraints(M->p[0], 2, P, MaxRays);
2205 evalue *EP = barvinok_enumerate_e(S, exist, nparam, MaxRays);
2206 Polyhedron_Free(S);
2207 Matrix_Free(M);
2208
2209 return enumerate_cyclic(P, exist, nparam, EP, 0, i, MaxRays);
2210 }
2211
2212 static int single_param_pos(Polyhedron*P, unsigned exist, unsigned nparam,
2213 int r)
2214 {
2215 int nvar = P->Dimension - exist - nparam;
2216 if (First_Non_Zero(P->Ray[r]+1, nvar) != -1)
2217 return -1;
2218 int i = First_Non_Zero(P->Ray[r]+1+nvar+exist, nparam);
2219 if (i == -1)
2220 return -1;
2221 if (First_Non_Zero(P->Ray[r]+1+nvar+exist+1, nparam-i-1) != -1)
2222 return -1;
2223 return i;
2224 }
2225
2226 static evalue* enumerate_remove_ray(Polyhedron *P, int r,
2227 unsigned exist, unsigned nparam, unsigned MaxRays)
2228 {
2229 #ifdef DEBUG_ER
2230 fprintf(stderr, "\nER: RedundantRay\n");
2231 #endif /* DEBUG_ER */
2232
2233 Value one;
2234 value_init(one);
2235 value_set_si(one, 1);
2236 int len = P->NbRays-1;
2237 Matrix *M = Matrix_Alloc(2 * len, P->Dimension+2);
2238 Vector_Copy(P->Ray[0], M->p[0], r * (P->Dimension+2));
2239 Vector_Copy(P->Ray[r+1], M->p[r], (len-r) * (P->Dimension+2));
2240 for (int j = 0; j < P->NbRays; ++j) {
2241 if (j == r)
2242 continue;
2243 Vector_Combine(P->Ray[j], P->Ray[r], M->p[len+j-(j>r)],
2244 one, P->Ray[j][P->Dimension+1], P->Dimension+2);
2245 }
2246
2247 P = Rays2Polyhedron(M, MaxRays);
2248 Matrix_Free(M);
2249 evalue *EP = barvinok_enumerate_e(P, exist, nparam, MaxRays);
2250 Polyhedron_Free(P);
2251 value_clear(one);
2252
2253 return EP;
2254 }
2255
2256 static evalue* enumerate_redundant_ray(Polyhedron *P,
2257 unsigned exist, unsigned nparam, unsigned MaxRays)
2258 {
2259 assert(P->NbBid == 0);
2260 int nvar = P->Dimension - exist - nparam;
2261 Value m;
2262 value_init(m);
2263
2264 for (int r = 0; r < P->NbRays; ++r) {
2265 if (value_notzero_p(P->Ray[r][P->Dimension+1]))
2266 continue;
2267 int i1 = single_param_pos(P, exist, nparam, r);
2268 if (i1 == -1)
2269 continue;
2270 for (int r2 = r+1; r2 < P->NbRays; ++r2) {
2271 if (value_notzero_p(P->Ray[r2][P->Dimension+1]))
2272 continue;
2273 int i2 = single_param_pos(P, exist, nparam, r2);
2274 if (i2 == -1)
2275 continue;
2276 if (i1 != i2)
2277 continue;
2278
2279 value_division(m, P->Ray[r][1+nvar+exist+i1],
2280 P->Ray[r2][1+nvar+exist+i1]);
2281 value_multiply(m, m, P->Ray[r2][1+nvar+exist+i1]);
2282 /* r2 divides r => r redundant */
2283 if (value_eq(m, P->Ray[r][1+nvar+exist+i1])) {
2284 value_clear(m);
2285 return enumerate_remove_ray(P, r, exist, nparam, MaxRays);
2286 }
2287
2288 value_division(m, P->Ray[r2][1+nvar+exist+i1],
2289 P->Ray[r][1+nvar+exist+i1]);
2290 value_multiply(m, m, P->Ray[r][1+nvar+exist+i1]);
2291 /* r divides r2 => r2 redundant */
2292 if (value_eq(m, P->Ray[r2][1+nvar+exist+i1])) {
2293 value_clear(m);
2294 return enumerate_remove_ray(P, r2, exist, nparam, MaxRays);
2295 }
2296 }
2297 }
2298 value_clear(m);
2299 return 0;
2300 }
2301
2302 static Polyhedron *upper_bound(Polyhedron *P,
2303 int pos, Value *max, Polyhedron **R)
2304 {
2305 Value v;
2306 int r;
2307 value_init(v);
2308
2309 *R = 0;
2310 Polyhedron *N;
2311 Polyhedron *B = 0;
2312 for (Polyhedron *Q = P; Q; Q = N) {
2313 N = Q->next;
2314 for (r = 0; r < P->NbRays; ++r) {
2315 if (value_zero_p(P->Ray[r][P->Dimension+1]) &&
2316 value_pos_p(P->Ray[r][1+pos]))
2317 break;
2318 }
2319 if (r < P->NbRays) {
2320 Q->next = *R;
2321 *R = Q;
2322 continue;
2323 } else {
2324 Q->next = B;
2325 B = Q;
2326 }
2327 for (r = 0; r < P->NbRays; ++r) {
2328 if (value_zero_p(P->Ray[r][P->Dimension+1]))
2329 continue;
2330 mpz_fdiv_q(v, P->Ray[r][1+pos], P->Ray[r][1+P->Dimension]);
2331 if ((!Q->next && r == 0) || value_gt(v, *max))
2332 value_assign(*max, v);
2333 }
2334 }
2335 value_clear(v);
2336 return B;
2337 }
2338
2339 static evalue* enumerate_ray(Polyhedron *P,
2340 unsigned exist, unsigned nparam, unsigned MaxRays)
2341 {
2342 assert(P->NbBid == 0);
2343 int nvar = P->Dimension - exist - nparam;
2344
2345 int r;
2346 for (r = 0; r < P->NbRays; ++r)
2347 if (value_zero_p(P->Ray[r][P->Dimension+1]))
2348 break;
2349 if (r >= P->NbRays)
2350 return 0;
2351
2352 int r2;
2353 for (r2 = r+1; r2 < P->NbRays; ++r2)
2354 if (value_zero_p(P->Ray[r2][P->Dimension+1]))
2355 break;
2356 if (r2 < P->NbRays) {
2357 if (nvar > 0)
2358 return enumerate_sum(P, exist, nparam, MaxRays);
2359 }
2360
2361 #ifdef DEBUG_ER
2362 fprintf(stderr, "\nER: Ray\n");
2363 #endif /* DEBUG_ER */
2364
2365 Value m;
2366 Value one;
2367 value_init(m);
2368 value_init(one);
2369 value_set_si(one, 1);
2370 int i = single_param_pos(P, exist, nparam, r);
2371 assert(i != -1); // for now;
2372
2373 Matrix *M = Matrix_Alloc(P->NbRays, P->Dimension+2);
2374 for (int j = 0; j < P->NbRays; ++j) {
2375 Vector_Combine(P->Ray[j], P->Ray[r], M->p[j],
2376 one, P->Ray[j][P->Dimension+1], P->Dimension+2);
2377 }
2378 Polyhedron *S = Rays2Polyhedron(M, MaxRays);
2379 Matrix_Free(M);
2380 Polyhedron *D = DomainDifference(P, S, MaxRays);
2381 Polyhedron_Free(S);
2382 // Polyhedron_Print(stderr, P_VALUE_FMT, D);
2383 assert(value_pos_p(P->Ray[r][1+nvar+exist+i])); // for now
2384 Polyhedron *R;
2385 D = upper_bound(D, nvar+exist+i, &m, &R);
2386 assert(D);
2387 Domain_Free(D);
2388
2389 M = Matrix_Alloc(2, P->Dimension+2);
2390 value_set_si(M->p[0][0], 1);
2391 value_set_si(M->p[1][0], 1);
2392 value_set_si(M->p[0][1+nvar+exist+i], -1);
2393 value_set_si(M->p[1][1+nvar+exist+i], 1);
2394 value_assign(M->p[0][1+P->Dimension], m);
2395 value_oppose(M->p[1][1+P->Dimension], m);
2396 value_addto(M->p[1][1+P->Dimension], M->p[1][1+P->Dimension],
2397 P->Ray[r][1+nvar+exist+i]);
2398 value_decrement(M->p[1][1+P->Dimension], M->p[1][1+P->Dimension]);
2399 // Matrix_Print(stderr, P_VALUE_FMT, M);
2400 D = AddConstraints(M->p[0], 2, P, MaxRays);
2401 // Polyhedron_Print(stderr, P_VALUE_FMT, D);
2402 value_substract(M->p[0][1+P->Dimension], M->p[0][1+P->Dimension],
2403 P->Ray[r][1+nvar+exist+i]);
2404 // Matrix_Print(stderr, P_VALUE_FMT, M);
2405 S = AddConstraints(M->p[0], 1, P, MaxRays);
2406 // Polyhedron_Print(stderr, P_VALUE_FMT, S);
2407 Matrix_Free(M);
2408
2409 evalue *EP = barvinok_enumerate_e(D, exist, nparam, MaxRays);
2410 Polyhedron_Free(D);
2411 value_clear(one);
2412 value_clear(m);
2413
2414 if (value_notone_p(P->Ray[r][1+nvar+exist+i]))
2415 EP = enumerate_cyclic(P, exist, nparam, EP, r, i, MaxRays);
2416 else {
2417 M = Matrix_Alloc(1, nparam+2);
2418 value_set_si(M->p[0][0], 1);
2419 value_set_si(M->p[0][1+i], 1);
2420 enumerate_vd_add_ray(EP, M, MaxRays);
2421 Matrix_Free(M);
2422 }
2423
2424 if (!emptyQ(S)) {
2425 evalue *E = barvinok_enumerate_e(S, exist, nparam, MaxRays);
2426 eadd(E, EP);
2427 free_evalue_refs(E);
2428 free(E);
2429 }
2430 Polyhedron_Free(S);
2431
2432 if (R) {
2433 assert(nvar == 0);
2434 evalue *ER = enumerate_or(R, exist, nparam, MaxRays);
2435 eor(ER, EP);
2436 free_evalue_refs(ER);
2437 free(ER);
2438 }
2439
2440 return EP;
2441 }
2442
2443 static evalue* new_zero_ep()
2444 {
2445 evalue *EP;
2446 ALLOC(EP);
2447 value_init(EP->d);
2448 evalue_set_si(EP, 0, 1);
2449 return EP;
2450 }
2451
2452 static evalue* enumerate_vd(Polyhedron **PA,
2453 unsigned exist, unsigned nparam, unsigned MaxRays)
2454 {
2455 Polyhedron *P = *PA;
2456 int nvar = P->Dimension - exist - nparam;
2457 Param_Polyhedron *PP = NULL;
2458 Polyhedron *C = Universe_Polyhedron(nparam);
2459 Polyhedron *CEq;
2460 Matrix *CT;
2461 Polyhedron *PR = P;
2462 PP = Polyhedron2Param_SimplifiedDomain(&PR,C,MaxRays,&CEq,&CT);
2463 Polyhedron_Free(C);
2464
2465 int nd;
2466 Param_Domain *D, *last;
2467 Value c;
2468 value_init(c);
2469 for (nd = 0, D=PP->D; D; D=D->next, ++nd)
2470 ;
2471
2472 Polyhedron **VD = new (Polyhedron*)[nd];
2473 Polyhedron **fVD = new (Polyhedron*)[nd];
2474 for(nd = 0, D=PP->D; D; D=D->next) {
2475 Polyhedron *rVD = reduce_domain(D->Domain, CT, CEq,
2476 fVD, nd, MaxRays);
2477 if (!rVD)
2478 continue;
2479
2480 VD[nd++] = rVD;
2481 last = D;
2482 }
2483
2484 evalue *EP = 0;
2485
2486 if (nd == 0)
2487 EP = new_zero_ep();
2488
2489 /* This doesn't seem to have any effect */
2490 if (nd == 1) {
2491 Polyhedron *CA = align_context(VD[0], P->Dimension, MaxRays);
2492 Polyhedron *O = P;
2493 P = DomainIntersection(P, CA, MaxRays);
2494 if (O != *PA)
2495 Polyhedron_Free(O);
2496 Polyhedron_Free(CA);
2497 if (emptyQ(P))
2498 EP = new_zero_ep();
2499 }
2500
2501 if (!EP && CT->NbColumns != CT->NbRows) {
2502 Polyhedron *CEqr = DomainImage(CEq, CT, MaxRays);
2503 Polyhedron *CA = align_context(CEqr, PR->Dimension, MaxRays);
2504 Polyhedron *I = DomainIntersection(PR, CA, MaxRays);
2505 Polyhedron_Free(CEqr);
2506 Polyhedron_Free(CA);
2507 #ifdef DEBUG_ER
2508 fprintf(stderr, "\nER: Eliminate\n");
2509 #endif /* DEBUG_ER */
2510 nparam -= CT->NbColumns - CT->NbRows;
2511 EP = barvinok_enumerate_e(I, exist, nparam, MaxRays);
2512 nparam += CT->NbColumns - CT->NbRows;
2513 addeliminatedparams_enum(EP, CT, CEq, MaxRays, nparam);
2514 Polyhedron_Free(I);
2515 }
2516 if (PR != *PA)
2517 Polyhedron_Free(PR);
2518 PR = 0;
2519
2520 if (!EP && nd > 1) {
2521 #ifdef DEBUG_ER
2522 fprintf(stderr, "\nER: VD\n");
2523 #endif /* DEBUG_ER */
2524 for (int i = 0; i < nd; ++i) {
2525 Polyhedron *CA = align_context(VD[i], P->Dimension, MaxRays);
2526 Polyhedron *I = DomainIntersection(P, CA, MaxRays);
2527
2528 if (i == 0)
2529 EP = barvinok_enumerate_e(I, exist, nparam, MaxRays);
2530 else {
2531 evalue *E = barvinok_enumerate_e(I, exist, nparam, MaxRays);
2532 eadd(E, EP);
2533 free_evalue_refs(E);
2534 free(E);
2535 }
2536 Polyhedron_Free(I);
2537 Polyhedron_Free(CA);
2538 }
2539 }
2540
2541 for (int i = 0; i < nd; ++i) {
2542 Polyhedron_Free(VD[i]);
2543 Polyhedron_Free(fVD[i]);
2544 }
2545 delete [] VD;
2546 delete [] fVD;
2547 value_clear(c);
2548
2549 if (!EP && nvar == 0) {
2550 Value f;
2551 value_init(f);
2552 Param_Vertices *V, *V2;
2553 Matrix* M = Matrix_Alloc(1, P->Dimension+2);
2554
2555 FORALL_PVertex_in_ParamPolyhedron(V, last, PP) {
2556 bool found = false;
2557 FORALL_PVertex_in_ParamPolyhedron(V2, last, PP) {
2558 if (V == V2) {
2559 found = true;
2560 continue;
2561 }
2562 if (!found)
2563 continue;
2564 for (int i = 0; i < exist; ++i) {
2565 value_oppose(f, V->Vertex->p[i][nparam+1]);
2566 Vector_Combine(V->Vertex->p[i],
2567 V2->Vertex->p[i],
2568 M->p[0] + 1 + nvar + exist,
2569 V2->Vertex->p[i][nparam+1],
2570 f,
2571 nparam+1);
2572 int j;
2573 for (j = 0; j < nparam; ++j)
2574 if (value_notzero_p(M->p[0][1+nvar+exist+j]))
2575 break;
2576 if (j >= nparam)
2577 continue;
2578 ConstraintSimplify(M->p[0], M->p[0],
2579 P->Dimension+2, &f);
2580 value_set_si(M->p[0][0], 0);
2581 Polyhedron *para = AddConstraints(M->p[0], 1, P,
2582 MaxRays);
2583 if (emptyQ(para)) {
2584 Polyhedron_Free(para);
2585 continue;
2586 }
2587 Polyhedron *pos, *neg;
2588 value_set_si(M->p[0][0], 1);
2589 value_decrement(M->p[0][P->Dimension+1],
2590 M->p[0][P->Dimension+1]);
2591 neg = AddConstraints(M->p[0], 1, P, MaxRays);
2592 value_set_si(f, -1);
2593 Vector_Scale(M->p[0]+1, M->p[0]+1, f,
2594 P->Dimension+1);
2595 value_decrement(M->p[0][P->Dimension+1],
2596 M->p[0][P->Dimension+1]);
2597 value_decrement(M->p[0][P->Dimension+1],
2598 M->p[0][P->Dimension+1]);
2599 pos = AddConstraints(M->p[0], 1, P, MaxRays);
2600 if (emptyQ(neg) && emptyQ(pos)) {
2601 Polyhedron_Free(para);
2602 Polyhedron_Free(pos);
2603 Polyhedron_Free(neg);
2604 continue;
2605 }
2606 #ifdef DEBUG_ER
2607 fprintf(stderr, "\nER: Order\n");
2608 #endif /* DEBUG_ER */
2609 EP = barvinok_enumerate_e(para, exist, nparam, MaxRays);
2610 evalue *E;
2611 if (!emptyQ(pos)) {
2612 E = barvinok_enumerate_e(pos, exist, nparam, MaxRays);
2613 eadd(E, EP);
2614 free_evalue_refs(E);
2615 free(E);
2616 }
2617 if (!emptyQ(neg)) {
2618 E = barvinok_enumerate_e(neg, exist, nparam, MaxRays);
2619 eadd(E, EP);
2620 free_evalue_refs(E);
2621 free(E);
2622 }
2623 Polyhedron_Free(para);
2624 Polyhedron_Free(pos);
2625 Polyhedron_Free(neg);
2626 break;
2627 }
2628 if (EP)
2629 break;
2630 } END_FORALL_PVertex_in_ParamPolyhedron;
2631 if (EP)
2632 break;
2633 } END_FORALL_PVertex_in_ParamPolyhedron;
2634
2635 if (!EP) {
2636 /* Search for vertex coordinate to split on */
2637 /* First look for one independent of the parameters */
2638 FORALL_PVertex_in_ParamPolyhedron(V, last, PP) {
2639 for (int i = 0; i < exist; ++i) {
2640 int j;
2641 for (j = 0; j < nparam; ++j)
2642 if (value_notzero_p(V->Vertex->p[i][j]))
2643 break;
2644 if (j < nparam)
2645 continue;
2646 value_set_si(M->p[0][0], 1);
2647 Vector_Set(M->p[0]+1, 0, nvar+exist);
2648 Vector_Copy(V->Vertex->p[i],
2649 M->p[0] + 1 + nvar + exist, nparam+1);
2650 value_oppose(M->p[0][1+nvar+i],
2651 V->Vertex->p[i][nparam+1]);
2652
2653 Polyhedron *pos, *neg;
2654 value_set_si(M->p[0][0], 1);
2655 value_decrement(M->p[0][P->Dimension+1],
2656 M->p[0][P->Dimension+1]);
2657 neg = AddConstraints(M->p[0], 1, P, MaxRays);
2658 value_set_si(f, -1);
2659 Vector_Scale(M->p[0]+1, M->p[0]+1, f,
2660 P->Dimension+1);
2661 value_decrement(M->p[0][P->Dimension+1],
2662 M->p[0][P->Dimension+1]);
2663 value_decrement(M->p[0][P->Dimension+1],
2664 M->p[0][P->Dimension+1]);
2665 pos = AddConstraints(M->p[0], 1, P, MaxRays);
2666 if (emptyQ(neg) || emptyQ(pos)) {
2667 Polyhedron_Free(pos);
2668 Polyhedron_Free(neg);
2669 continue;
2670 }
2671 Polyhedron_Free(pos);
2672 value_increment(M->p[0][P->Dimension+1],
2673 M->p[0][P->Dimension+1]);
2674 pos = AddConstraints(M->p[0], 1, P, MaxRays);
2675 #ifdef DEBUG_ER
2676 fprintf(stderr, "\nER: Vertex\n");
2677 #endif /* DEBUG_ER */
2678 pos->next = neg;
2679 EP = enumerate_or(pos, exist, nparam, MaxRays);
2680 break;
2681 }
2682 if (EP)
2683 break;
2684 } END_FORALL_PVertex_in_ParamPolyhedron;
2685 }
2686
2687 if (!EP) {
2688 /* Search for vertex coordinate to split on */
2689 /* Now look for one that depends on the parameters */
2690 FORALL_PVertex_in_ParamPolyhedron(V, last, PP) {
2691 for (int i = 0; i < exist; ++i) {
2692 value_set_si(M->p[0][0], 1);
2693 Vector_Set(M->p[0]+1, 0, nvar+exist);
2694 Vector_Copy(V->Vertex->p[i],
2695 M->p[0] + 1 + nvar + exist, nparam+1);
2696 value_oppose(M->p[0][1+nvar+i],
2697 V->Vertex->p[i][nparam+1]);
2698
2699 Polyhedron *pos, *neg;
2700 value_set_si(M->p[0][0], 1);
2701 value_decrement(M->p[0][P->Dimension+1],
2702 M->p[0][P->Dimension+1]);
2703 neg = AddConstraints(M->p[0], 1, P, MaxRays);
2704 value_set_si(f, -1);
2705 Vector_Scale(M->p[0]+1, M->p[0]+1, f,
2706 P->Dimension+1);
2707 value_decrement(M->p[0][P->Dimension+1],
2708 M->p[0][P->Dimension+1]);
2709 value_decrement(M->p[0][P->Dimension+1],
2710 M->p[0][P->Dimension+1]);
2711 pos = AddConstraints(M->p[0], 1, P, MaxRays);
2712 if (emptyQ(neg) || emptyQ(pos)) {
2713 Polyhedron_Free(pos);
2714 Polyhedron_Free(neg);
2715 continue;
2716 }
2717 Polyhedron_Free(pos);
2718 value_increment(M->p[0][P->Dimension+1],
2719 M->p[0][P->Dimension+1]);
2720 pos = AddConstraints(M->p[0], 1, P, MaxRays);
2721 #ifdef DEBUG_ER
2722 fprintf(stderr, "\nER: ParamVertex\n");
2723 #endif /* DEBUG_ER */
2724 pos->next = neg;
2725 EP = enumerate_or(pos, exist, nparam, MaxRays);
2726 break;
2727 }
2728 if (EP)
2729 break;
2730 } END_FORALL_PVertex_in_ParamPolyhedron;
2731 }
2732
2733 Matrix_Free(M);
2734 value_clear(f);
2735 }
2736
2737 if (CEq)
2738 Polyhedron_Free(CEq);
2739 if (CT)
2740 Matrix_Free(CT);
2741 if (PP)
2742 Param_Polyhedron_Free(PP);
2743 *PA = P;
2744
2745 return EP;
2746 }
2747
2748 static bool is_single(Value *row, int pos, int len)
2749 {
2750 return First_Non_Zero(row, pos) == -1 &&
2751 First_Non_Zero(row+pos+1, len-pos-1) == -1;
2752 }
2753
2754 static evalue* barvinok_enumerate_e_r(Polyhedron *P,
2755 unsigned exist, unsigned nparam, unsigned MaxRays);
2756
2757 #ifdef DEBUG_ER
2758 static int er_level = 0;
2759
2760 evalue* barvinok_enumerate_e(Polyhedron *P,
2761 unsigned exist, unsigned nparam, unsigned MaxRays)
2762 {
2763 fprintf(stderr, "\nER: level %i\n", er_level);
2764 int nvar = P->Dimension - exist - nparam;
2765 fprintf(stderr, "%d %d %d\n", nvar, exist, nparam);
2766
2767 Polyhedron_Print(stderr, P_VALUE_FMT, P);
2768 ++er_level;
2769 P = DomainConstraintSimplify(Polyhedron_Copy(P), MaxRays);
2770 evalue *EP = barvinok_enumerate_e_r(P, exist, nparam, MaxRays);
2771 Polyhedron_Free(P);
2772 --er_level;
2773 return EP;
2774 }
2775 #else
2776 evalue* barvinok_enumerate_e(Polyhedron *P,
2777 unsigned exist, unsigned nparam, unsigned MaxRays)
2778 {
2779 P = DomainConstraintSimplify(Polyhedron_Copy(P), MaxRays);
2780 evalue *EP = barvinok_enumerate_e_r(P, exist, nparam, MaxRays);
2781 Polyhedron_Free(P);
2782 return EP;
2783 }
2784 #endif
2785
2786 static evalue* barvinok_enumerate_e_r(Polyhedron *P,
2787 unsigned exist, unsigned nparam, unsigned MaxRays)
2788 {
2789 if (exist == 0) {
2790 Polyhedron *U = Universe_Polyhedron(nparam);
2791 evalue *EP = barvinok_enumerate_ev(P, U, MaxRays);
2792 //char *param_name[] = {"P", "Q", "R", "S", "T" };
2793 //print_evalue(stdout, EP, param_name);
2794 Polyhedron_Free(U);
2795 return EP;
2796 }
2797
2798 int nvar = P->Dimension - exist - nparam;
2799 int len = P->Dimension + 2;
2800
2801 if (emptyQ(P))
2802 return new_zero_ep();
2803
2804 if (nvar == 0 && nparam == 0) {
2805 evalue *EP = new_zero_ep();
2806 barvinok_count(P, &EP->x.n, MaxRays);
2807 if (value_pos_p(EP->x.n))
2808 value_set_si(EP->x.n, 1);
2809 return EP;
2810 }
2811
2812 int r;
2813 for (r = 0; r < P->NbRays; ++r)
2814 if (value_zero_p(P->Ray[r][0]) ||
2815 value_zero_p(P->Ray[r][P->Dimension+1])) {
2816 int i;
2817 for (i = 0; i < nvar; ++i)
2818 if (value_notzero_p(P->Ray[r][i+1]))
2819 break;
2820 if (i >= nvar)
2821 continue;
2822 for (i = nvar + exist; i < nvar + exist + nparam; ++i)
2823 if (value_notzero_p(P->Ray[r][i+1]))
2824 break;
2825 if (i >= nvar + exist + nparam)
2826 break;
2827 }
2828 if (r < P->NbRays) {
2829 evalue *EP = new_zero_ep();
2830 value_set_si(EP->x.n, -1);
2831 return EP;
2832 }
2833
2834 int first;
2835 for (r = 0; r < P->NbEq; ++r)
2836 if ((first = First_Non_Zero(P->Constraint[r]+1+nvar, exist)) != -1)
2837 break;
2838 if (r < P->NbEq) {
2839 if (First_Non_Zero(P->Constraint[r]+1+nvar+first+1,
2840 exist-first-1) != -1) {
2841 Polyhedron *T = rotate_along(P, r, nvar, exist, MaxRays);
2842 #ifdef DEBUG_ER
2843 fprintf(stderr, "\nER: Equality\n");
2844 #endif /* DEBUG_ER */
2845 evalue *EP = barvinok_enumerate_e(T, exist-1, nparam, MaxRays);
2846 Polyhedron_Free(T);
2847 return EP;
2848 } else {
2849 #ifdef DEBUG_ER
2850 fprintf(stderr, "\nER: Fixed\n");
2851 #endif /* DEBUG_ER */
2852 if (first == 0)
2853 return barvinok_enumerate_e(P, exist-1, nparam, MaxRays);
2854 else {
2855 Polyhedron *T = Polyhedron_Copy(P);
2856 SwapColumns(T, nvar+1, nvar+1+first);
2857 evalue *EP = barvinok_enumerate_e(T, exist-1, nparam, MaxRays);
2858 Polyhedron_Free(T);
2859 return EP;
2860 }
2861 }
2862 }
2863
2864 Vector *row = Vector_Alloc(len);
2865 value_set_si(row->p[0], 1);
2866
2867 Value f;
2868 value_init(f);
2869
2870 enum constraint info[exist];
2871 for (int i = 0; i < exist; ++i) {
2872 info[i] = ALL_POS;
2873 for (int l = P->NbEq; l < P->NbConstraints; ++l) {
2874 if (value_negz_p(P->Constraint[l][nvar+i+1]))
2875 continue;
2876 bool l_parallel = is_single(P->Constraint[l]+nvar+1, i, exist);
2877 for (int u = P->NbEq; u < P->NbConstraints; ++u) {
2878 if (value_posz_p(P->Constraint[u][nvar+i+1]))
2879 continue;
2880 bool lu_parallel = l_parallel ||
2881 is_single(P->Constraint[u]+nvar+1, i, exist);
2882 value_oppose(f, P->Constraint[u][nvar+i+1]);
2883 Vector_Combine(P->Constraint[l]+1, P->Constraint[u]+1, row->p+1,
2884 f, P->Constraint[l][nvar+i+1], len-1);
2885 if (!(info[i] & INDEPENDENT)) {
2886 int j;
2887 for (j = 0; j < exist; ++j)
2888 if (j != i && value_notzero_p(row->p[nvar+j+1]))
2889 break;
2890 if (j == exist) {
2891 //printf("independent: i: %d, l: %d, u: %d\n", i, l, u);
2892 info[i] = (constraint)(info[i] | INDEPENDENT);
2893 }
2894 }
2895 if (info[i] & ALL_POS) {
2896 value_addto(row->p[len-1], row->p[len-1],
2897 P->Constraint[l][nvar+i+1]);
2898 value_addto(row->p[len-1], row->p[len-1], f);
2899 value_multiply(f, f, P->Constraint[l][nvar+i+1]);
2900 value_substract(row->p[len-1], row->p[len-1], f);
2901 value_decrement(row->p[len-1], row->p[len-1]);
2902 ConstraintSimplify(row->p, row->p, len, &f);
2903 value_set_si(f, -1);
2904 Vector_Scale(row->p+1, row->p+1, f, len-1);
2905 value_decrement(row->p[len-1], row->p[len-1]);
2906 Polyhedron *T = AddConstraints(row->p, 1, P, MaxRays);
2907 if (!emptyQ(T)) {
2908 //printf("not all_pos: i: %d, l: %d, u: %d\n", i, l, u);
2909 info[i] = (constraint)(info[i] ^ ALL_POS);
2910 }
2911 //puts("pos remainder");
2912 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
2913 Polyhedron_Free(T);
2914 }
2915 if (!(info[i] & ONE_NEG)) {
2916 if (lu_parallel) {
2917 negative_test_constraint(P->Constraint[l],
2918 P->Constraint[u],
2919 row->p, nvar+i, len, &f);
2920 oppose_constraint(row->p, len, &f);
2921 Polyhedron *T = AddConstraints(row->p, 1, P, MaxRays);
2922 if (emptyQ(T)) {
2923 //printf("one_neg i: %d, l: %d, u: %d\n", i, l, u);
2924 info[i] = (constraint)(info[i] | ONE_NEG);
2925 }
2926 //puts("neg remainder");
2927 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
2928 Polyhedron_Free(T);
2929 }
2930 }
2931 if (!(info[i] & ALL_POS) && (info[i] & ONE_NEG))
2932 goto next;
2933 }
2934 }
2935 if (info[i] & ALL_POS)
2936 break;
2937 next:
2938 ;
2939 }
2940
2941 /*
2942 for (int i = 0; i < exist; ++i)
2943 printf("%i: %i\n", i, info[i]);
2944 */
2945 for (int i = 0; i < exist; ++i)
2946 if (info[i] & ALL_POS) {
2947 #ifdef DEBUG_ER
2948 fprintf(stderr, "\nER: Positive\n");
2949 #endif /* DEBUG_ER */
2950 // Eliminate
2951 // Maybe we should chew off some of the fat here
2952 Matrix *M = Matrix_Alloc(P->Dimension, P->Dimension+1);
2953 for (int j = 0; j < P->Dimension; ++j)
2954 value_set_si(M->p[j][j + (j >= i+nvar)], 1);
2955 Polyhedron *T = Polyhedron_Image(P, M, MaxRays);
2956 Matrix_Free(M);
2957 evalue *EP = barvinok_enumerate_e(T, exist-1, nparam, MaxRays);
2958 Polyhedron_Free(T);
2959 value_clear(f);
2960 Vector_Free(row);
2961 return EP;
2962 }
2963 for (int i = 0; i < exist; ++i)
2964 if (info[i] & ONE_NEG) {
2965 #ifdef DEBUG_ER
2966 fprintf(stderr, "\nER: Negative\n");
2967 #endif /* DEBUG_ER */
2968 Vector_Free(row);
2969 value_clear(f);
2970 if (i == 0)
2971 return barvinok_enumerate_e(P, exist-1, nparam, MaxRays);
2972 else {
2973 Polyhedron *T = Polyhedron_Copy(P);
2974 SwapColumns(T, nvar+1, nvar+1+i);
2975 evalue *EP = barvinok_enumerate_e(T, exist-1, nparam, MaxRays);
2976 Polyhedron_Free(T);
2977 return EP;
2978 }
2979 }
2980 for (int i = 0; i < exist; ++i)
2981 if (info[i] & INDEPENDENT) {
2982 Polyhedron *pos, *neg;
2983
2984 /* Find constraint again and split off negative part */
2985
2986 if (SplitOnVar(P, i, nvar, len, exist, MaxRays,
2987 row, f, true, &pos, &neg)) {
2988 #ifdef DEBUG_ER
2989 fprintf(stderr, "\nER: Split\n");
2990 #endif /* DEBUG_ER */
2991
2992 evalue *EP =
2993 barvinok_enumerate_e(neg, exist-1, nparam, MaxRays);
2994 evalue *E =
2995 barvinok_enumerate_e(pos, exist, nparam, MaxRays);
2996 eadd(E, EP);
2997 free_evalue_refs(E);
2998 free(E);
2999 Polyhedron_Free(neg);
3000 Polyhedron_Free(pos);
3001 value_clear(f);
3002 Vector_Free(row);
3003 return EP;
3004 }
3005 }
3006
3007 Polyhedron *O = P;
3008 Polyhedron *F;
3009
3010 evalue *EP;
3011
3012 EP = enumerate_line(P, exist, nparam, MaxRays);
3013 if (EP)
3014 return EP;
3015
3016 EP = enumerate_redundant_ray(P, exist, nparam, MaxRays);
3017 if (EP)
3018 goto out;
3019
3020 EP = enumerate_sure(P, exist, nparam, MaxRays);
3021 if (EP)
3022 goto out;
3023
3024 EP = enumerate_ray(P, exist, nparam, MaxRays);
3025 if (EP)
3026 goto out;
3027
3028 EP = enumerate_sure2(P, exist, nparam, MaxRays);
3029 if (EP)
3030 goto out;
3031
3032 F = unfringe(P, MaxRays);
3033 if (!PolyhedronIncludes(F, P)) {
3034 #ifdef DEBUG_ER
3035 fprintf(stderr, "\nER: Fringed\n");
3036 #endif /* DEBUG_ER */
3037 EP = barvinok_enumerate_e(F, exist, nparam, MaxRays);
3038 Polyhedron_Free(F);
3039 goto out;
3040 }
3041 Polyhedron_Free(F);
3042
3043 if (nparam)
3044 EP = enumerate_vd(&P, exist, nparam, MaxRays);
3045 if (EP)
3046 goto out2;
3047
3048 if (nvar != 0) {
3049 EP = enumerate_sum(P, exist, nparam, MaxRays);
3050 goto out2;
3051 }
3052
3053 assert(nvar == 0);
3054
3055 int i;
3056 Polyhedron *pos, *neg;
3057 for (i = 0; i < exist; ++i)
3058 if (SplitOnVar(P, i, nvar, len, exist, MaxRays,
3059 row, f, false, &pos, &neg))
3060 break;
3061
3062 assert (i < exist);
3063
3064 pos->next = neg;
3065 EP = enumerate_or(pos, exist, nparam, MaxRays);
3066
3067 out2:
3068 if (O != P)
3069 Polyhedron_Free(P);
3070
3071 out:
3072 value_clear(f);
3073 Vector_Free(row);
3074 return EP;
3075 }
lattice point enumeration library