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