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