?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
barvinok_enumerate: make sure the input polyhedra are fully specified.
[barvinok.git] / barvinok.cc
blobd1c75d7dfc8f15c8416fb2b1e2bf067c5701af83
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::cerr;
24 using std::cout;
25 using std::endl;
26 using std::vector;
27 using std::deque;
28 using std::string;
29 using std::ostringstream;
30
31 #define ALLOC(p) (((long *) (p))[0])
32 #define SIZE(p) (((long *) (p))[1])
33 #define DATA(p) ((mp_limb_t *) (((long *) (p)) + 2))
34
35 static void value2zz(Value v, ZZ& z)
36 {
37 int sa = v[0]._mp_size;
38 int abs_sa = sa < 0 ? -sa : sa;
39
40 _ntl_gsetlength(&z.rep, abs_sa);
41 mp_limb_t * adata = DATA(z.rep);
42 for (int i = 0; i < abs_sa; ++i)
43 adata[i] = v[0]._mp_d[i];
44 SIZE(z.rep) = sa;
45 }
46
47 void zz2value(ZZ& z, Value& v)
48 {
49 if (!z.rep) {
50 value_set_si(v, 0);
51 return;
52 }
53
54 int sa = SIZE(z.rep);
55 int abs_sa = sa < 0 ? -sa : sa;
56
57 mp_limb_t * adata = DATA(z.rep);
58 _mpz_realloc(v, abs_sa);
59 for (int i = 0; i < abs_sa; ++i)
60 v[0]._mp_d[i] = adata[i];
61 v[0]._mp_size = sa;
62 }
63
64 #undef ALLOC
65 #define ALLOC(t,p) p = (t*)malloc(sizeof(*p))
66
67 /*
68 * We just ignore the last column and row
69 * If the final element is not equal to one
70 * then the result will actually be a multiple of the input
71 */
72 static void matrix2zz(Matrix *M, mat_ZZ& m, unsigned nr, unsigned nc)
73 {
74 m.SetDims(nr, nc);
75
76 for (int i = 0; i < nr; ++i) {
77 // assert(value_one_p(M->p[i][M->NbColumns - 1]));
78 for (int j = 0; j < nc; ++j) {
79 value2zz(M->p[i][j], m[i][j]);
80 }
81 }
82 }
83
84 static void values2zz(Value *p, vec_ZZ& v, int len)
85 {
86 v.SetLength(len);
87
88 for (int i = 0; i < len; ++i) {
89 value2zz(p[i], v[i]);
90 }
91 }
92
93 /*
94 */
95 static void zz2values(vec_ZZ& v, Value *p)
96 {
97 for (int i = 0; i < v.length(); ++i)
98 zz2value(v[i], p[i]);
99 }
100
101 static void rays(mat_ZZ& r, Polyhedron *C)
102 {
103 unsigned dim = C->NbRays - 1; /* don't count zero vertex */
104 assert(C->NbRays - 1 == C->Dimension);
105 r.SetDims(dim, dim);
106 ZZ tmp;
107
108 int i, c;
109 for (i = 0, c = 0; i < dim; ++i)
110 if (value_zero_p(C->Ray[i][dim+1])) {
111 for (int j = 0; j < dim; ++j) {
112 value2zz(C->Ray[i][j+1], tmp);
113 r[j][c] = tmp;
114 }
115 ++c;
116 }
117 }
118
119 static Matrix * rays(Polyhedron *C)
120 {
121 unsigned dim = C->NbRays - 1; /* don't count zero vertex */
122 assert(C->NbRays - 1 == C->Dimension);
123
124 Matrix *M = Matrix_Alloc(dim+1, dim+1);
125 assert(M);
126
127 int i, c;
128 for (i = 0, c = 0; i <= dim && c < dim; ++i)
129 if (value_zero_p(C->Ray[i][dim+1])) {
130 Vector_Copy(C->Ray[i] + 1, M->p[c], dim);
131 value_set_si(M->p[c++][dim], 0);
132 }
133 assert(c == dim);
134 value_set_si(M->p[dim][dim], 1);
135
136 return M;
137 }
138
139 static Matrix * rays2(Polyhedron *C)
140 {
141 unsigned dim = C->NbRays - 1; /* don't count zero vertex */
142 assert(C->NbRays - 1 == C->Dimension);
143
144 Matrix *M = Matrix_Alloc(dim, dim);
145 assert(M);
146
147 int i, c;
148 for (i = 0, c = 0; i <= dim && c < dim; ++i)
149 if (value_zero_p(C->Ray[i][dim+1]))
150 Vector_Copy(C->Ray[i] + 1, M->p[c++], dim);
151 assert(c == dim);
152
153 return M;
154 }
155
156 /*
157 * Returns the largest absolute value in the vector
158 */
159 static ZZ max(vec_ZZ& v)
160 {
161 ZZ max = abs(v[0]);
162 for (int i = 1; i < v.length(); ++i)
163 if (abs(v[i]) > max)
164 max = abs(v[i]);
165 return max;
166 }
167
168 class cone {
169 public:
170 cone(Matrix *M) {
171 Cone = 0;
172 Rays = Matrix_Copy(M);
173 set_det();
174 }
175 cone(Polyhedron *C) {
176 Cone = Polyhedron_Copy(C);
177 Rays = rays(C);
178 set_det();
179 }
180 void set_det() {
181 mat_ZZ A;
182 matrix2zz(Rays, A, Rays->NbRows - 1, Rays->NbColumns - 1);
183 det = determinant(A);
184 }
185
186 Vector* short_vector(vec_ZZ& lambda) {
187 Matrix *M = Matrix_Copy(Rays);
188 Matrix *inv = Matrix_Alloc(M->NbRows, M->NbColumns);
189 int ok = Matrix_Inverse(M, inv);
190 assert(ok);
191 Matrix_Free(M);
192
193 ZZ det2;
194 mat_ZZ B;
195 mat_ZZ U;
196 matrix2zz(inv, B, inv->NbRows - 1, inv->NbColumns - 1);
197 long r = LLL(det2, B, U);
198
199 ZZ min = max(B[0]);
200 int index = 0;
201 for (int i = 1; i < B.NumRows(); ++i) {
202 ZZ tmp = max(B[i]);
203 if (tmp < min) {
204 min = tmp;
205 index = i;
206 }
207 }
208
209 Matrix_Free(inv);
210
211 lambda = B[index];
212
213 Vector *z = Vector_Alloc(U[index].length()+1);
214 assert(z);
215 zz2values(U[index], z->p);
216 value_set_si(z->p[U[index].length()], 0);
217
218 Polyhedron *C = poly();
219 int i;
220 for (i = 0; i < lambda.length(); ++i)
221 if (lambda[i] > 0)
222 break;
223 if (i == lambda.length()) {
224 Value tmp;
225 value_init(tmp);
226 value_set_si(tmp, -1);
227 Vector_Scale(z->p, z->p, tmp, z->Size-1);
228 value_clear(tmp);
229 }
230 return z;
231 }
232
233 ~cone() {
234 Polyhedron_Free(Cone);
235 Matrix_Free(Rays);
236 }
237
238 Polyhedron *poly() {
239 if (!Cone) {
240 Matrix *M = Matrix_Alloc(Rays->NbRows+1, Rays->NbColumns+1);
241 for (int i = 0; i < Rays->NbRows; ++i) {
242 Vector_Copy(Rays->p[i], M->p[i]+1, Rays->NbColumns);
243 value_set_si(M->p[i][0], 1);
244 }
245 Vector_Set(M->p[Rays->NbRows]+1, 0, Rays->NbColumns-1);
246 value_set_si(M->p[Rays->NbRows][0], 1);
247 value_set_si(M->p[Rays->NbRows][Rays->NbColumns], 1);
248 Cone = Rays2Polyhedron(M, M->NbRows+1);
249 assert(Cone->NbConstraints == Cone->NbRays);
250 Matrix_Free(M);
251 }
252 return Cone;
253 }
254
255 ZZ det;
256 Polyhedron *Cone;
257 Matrix *Rays;
258 };
259
260 class dpoly {
261 public:
262 vec_ZZ coeff;
263 dpoly(int d, ZZ& degree, int offset = 0) {
264 coeff.SetLength(d+1);
265
266 int min = d + offset;
267 if (degree >= 0 && degree < ZZ(INIT_VAL, min))
268 min = to_int(degree);
269
270 ZZ c = ZZ(INIT_VAL, 1);
271 if (!offset)
272 coeff[0] = c;
273 for (int i = 1; i <= min; ++i) {
274 c *= (degree -i + 1);
275 c /= i;
276 coeff[i-offset] = c;
277 }
278 }
279 void operator *= (dpoly& f) {
280 assert(coeff.length() == f.coeff.length());
281 vec_ZZ old = coeff;
282 coeff = f.coeff[0] * coeff;
283 for (int i = 1; i < coeff.length(); ++i)
284 for (int j = 0; i+j < coeff.length(); ++j)
285 coeff[i+j] += f.coeff[i] * old[j];
286 }
287 void div(dpoly& d, mpq_t count, ZZ& sign) {
288 int len = coeff.length();
289 Value tmp;
290 value_init(tmp);
291 mpq_t* c = new mpq_t[coeff.length()];
292 mpq_t qtmp;
293 mpq_init(qtmp);
294 for (int i = 0; i < len; ++i) {
295 mpq_init(c[i]);
296 zz2value(coeff[i], tmp);
297 mpq_set_z(c[i], tmp);
298
299 for (int j = 1; j <= i; ++j) {
300 zz2value(d.coeff[j], tmp);
301 mpq_set_z(qtmp, tmp);
302 mpq_mul(qtmp, qtmp, c[i-j]);
303 mpq_sub(c[i], c[i], qtmp);
304 }
305
306 zz2value(d.coeff[0], tmp);
307 mpq_set_z(qtmp, tmp);
308 mpq_div(c[i], c[i], qtmp);
309 }
310 if (sign == -1)
311 mpq_sub(count, count, c[len-1]);
312 else
313 mpq_add(count, count, c[len-1]);
314
315 value_clear(tmp);
316 mpq_clear(qtmp);
317 for (int i = 0; i < len; ++i)
318 mpq_clear(c[i]);
319 delete [] c;
320 }
321 };
322
323 class dpoly_n {
324 public:
325 Matrix *coeff;
326 ~dpoly_n() {
327 Matrix_Free(coeff);
328 }
329 dpoly_n(int d, ZZ& degree_0, ZZ& degree_1, int offset = 0) {
330 Value d0, d1;
331 value_init(d0);
332 value_init(d1);
333 zz2value(degree_0, d0);
334 zz2value(degree_1, d1);
335 coeff = Matrix_Alloc(d+1, d+1+1);
336 value_set_si(coeff->p[0][0], 1);
337 value_set_si(coeff->p[0][d+1], 1);
338 for (int i = 1; i <= d; ++i) {
339 value_multiply(coeff->p[i][0], coeff->p[i-1][0], d0);
340 Vector_Combine(coeff->p[i-1], coeff->p[i-1]+1, coeff->p[i]+1,
341 d1, d0, i);
342 value_set_si(coeff->p[i][d+1], i);
343 value_multiply(coeff->p[i][d+1], coeff->p[i][d+1], coeff->p[i-1][d+1]);
344 value_decrement(d0, d0);
345 }
346 value_clear(d0);
347 value_clear(d1);
348 }
349 void div(dpoly& d, Vector *count, ZZ& sign) {
350 int len = coeff->NbRows;
351 Matrix * c = Matrix_Alloc(coeff->NbRows, coeff->NbColumns);
352 Value tmp;
353 value_init(tmp);
354 for (int i = 0; i < len; ++i) {
355 Vector_Copy(coeff->p[i], c->p[i], len+1);
356 for (int j = 1; j <= i; ++j) {
357 zz2value(d.coeff[j], tmp);
358 value_multiply(tmp, tmp, c->p[i][len]);
359 value_oppose(tmp, tmp);
360 Vector_Combine(c->p[i], c->p[i-j], c->p[i],
361 c->p[i-j][len], tmp, len);
362 value_multiply(c->p[i][len], c->p[i][len], c->p[i-j][len]);
363 }
364 zz2value(d.coeff[0], tmp);
365 value_multiply(c->p[i][len], c->p[i][len], tmp);
366 }
367 if (sign == -1) {
368 value_set_si(tmp, -1);
369 Vector_Scale(c->p[len-1], count->p, tmp, len);
370 value_assign(count->p[len], c->p[len-1][len]);
371 } else
372 Vector_Copy(c->p[len-1], count->p, len+1);
373 Vector_Normalize(count->p, len+1);
374 value_clear(tmp);
375 Matrix_Free(c);
376 }
377 };
378
379 struct dpoly_r_term {
380 int *powers;
381 ZZ coeff;
382 };
383
384 /* len: number of elements in c
385 * each element in c is the coefficient of a power of t
386 * in the MacLaurin expansion
387 */
388 struct dpoly_r {
389 vector< dpoly_r_term * > *c;
390 int len;
391 int dim;
392 ZZ denom;
393
394 void add_term(int i, int * powers, ZZ& coeff) {
395 if (coeff == 0)
396 return;
397 for (int k = 0; k < c[i].size(); ++k) {
398 if (memcmp(c[i][k]->powers, powers, dim * sizeof(int)) == 0) {
399 c[i][k]->coeff += coeff;
400 return;
401 }
402 }
403 dpoly_r_term *t = new dpoly_r_term;
404 t->powers = new int[dim];
405 memcpy(t->powers, powers, dim * sizeof(int));
406 t->coeff = coeff;
407 c[i].push_back(t);
408 }
409 dpoly_r(int len, int dim) {
410 denom = 1;
411 this->len = len;
412 this->dim = dim;
413 c = new vector< dpoly_r_term * > [len];
414 }
415 dpoly_r(dpoly& num, int dim) {
416 denom = 1;
417 len = num.coeff.length();
418 c = new vector< dpoly_r_term * > [len];
419 this->dim = dim;
420 int powers[dim];
421 memset(powers, 0, dim * sizeof(int));
422
423 for (int i = 0; i < len; ++i) {
424 ZZ coeff = num.coeff[i];
425 add_term(i, powers, coeff);
426 }
427 }
428 dpoly_r(dpoly& num, dpoly& den, int pos, int dim) {
429 denom = 1;
430 len = num.coeff.length();
431 c = new vector< dpoly_r_term * > [len];
432 this->dim = dim;
433 int powers[dim];
434
435 for (int i = 0; i < len; ++i) {
436 ZZ coeff = num.coeff[i];
437 memset(powers, 0, dim * sizeof(int));
438 powers[pos] = 1;
439
440 add_term(i, powers, coeff);
441
442 for (int j = 1; j <= i; ++j) {
443 for (int k = 0; k < c[i-j].size(); ++k) {
444 memcpy(powers, c[i-j][k]->powers, dim*sizeof(int));
445 powers[pos]++;
446 coeff = -den.coeff[j-1] * c[i-j][k]->coeff;
447 add_term(i, powers, coeff);
448 }
449 }
450 }
451 //dump();
452 }
453 dpoly_r(dpoly_r* num, dpoly& den, int pos, int dim) {
454 denom = num->denom;
455 len = num->len;
456 c = new vector< dpoly_r_term * > [len];
457 this->dim = dim;
458 int powers[dim];
459 ZZ coeff;
460
461 for (int i = 0 ; i < len; ++i) {
462 for (int k = 0; k < num->c[i].size(); ++k) {
463 memcpy(powers, num->c[i][k]->powers, dim*sizeof(int));
464 powers[pos]++;
465 add_term(i, powers, num->c[i][k]->coeff);
466 }
467
468 for (int j = 1; j <= i; ++j) {
469 for (int k = 0; k < c[i-j].size(); ++k) {
470 memcpy(powers, c[i-j][k]->powers, dim*sizeof(int));
471 powers[pos]++;
472 coeff = -den.coeff[j-1] * c[i-j][k]->coeff;
473 add_term(i, powers, coeff);
474 }
475 }
476 }
477 }
478 ~dpoly_r() {
479 for (int i = 0 ; i < len; ++i)
480 for (int k = 0; k < c[i].size(); ++k) {
481 delete [] c[i][k]->powers;
482 delete c[i][k];
483 }
484 delete [] c;
485 }
486 dpoly_r *div(dpoly& d) {
487 dpoly_r *rc = new dpoly_r(len, dim);
488 rc->denom = power(d.coeff[0], len);
489 ZZ inv_d = rc->denom / d.coeff[0];
490 ZZ coeff;
491
492 for (int i = 0; i < len; ++i) {
493 for (int k = 0; k < c[i].size(); ++k) {
494 coeff = c[i][k]->coeff * inv_d;
495 rc->add_term(i, c[i][k]->powers, coeff);
496 }
497
498 for (int j = 1; j <= i; ++j) {
499 for (int k = 0; k < rc->c[i-j].size(); ++k) {
500 coeff = - d.coeff[j] * rc->c[i-j][k]->coeff / d.coeff[0];
501 rc->add_term(i, rc->c[i-j][k]->powers, coeff);
502 }
503 }
504 }
505 return rc;
506 }
507 void dump(void) {
508 for (int i = 0; i < len; ++i) {
509 cerr << endl;
510 cerr << i << endl;
511 cerr << c[i].size() << endl;
512 for (int j = 0; j < c[i].size(); ++j) {
513 for (int k = 0; k < dim; ++k) {
514 cerr << c[i][j]->powers[k] << " ";
515 }
516 cerr << ": " << c[i][j]->coeff << "/" << denom << endl;
517 }
518 cerr << endl;
519 }
520 }
521 };
522
523 struct decomposer {
524 void decompose(Polyhedron *C);
525 virtual void handle(Polyhedron *P, int sign) = 0;
526 };
527
528 struct polar_decomposer : public decomposer {
529 void decompose(Polyhedron *C, unsigned MaxRays);
530 virtual void handle(Polyhedron *P, int sign);
531 virtual void handle_polar(Polyhedron *P, int sign) = 0;
532 };
533
534 void decomposer::decompose(Polyhedron *C)
535 {
536 vector<cone *> nonuni;
537 cone * c = new cone(C);
538 ZZ det = c->det;
539 int s = sign(det);
540 assert(det != 0);
541 if (abs(det) > 1) {
542 nonuni.push_back(c);
543 } else {
544 try {
545 handle(C, 1);
546 delete c;
547 } catch (...) {
548 delete c;
549 throw;
550 }
551 }
552 vec_ZZ lambda;
553 while (!nonuni.empty()) {
554 c = nonuni.back();
555 nonuni.pop_back();
556 Vector* v = c->short_vector(lambda);
557 for (int i = 0; i < c->Rays->NbRows - 1; ++i) {
558 if (lambda[i] == 0)
559 continue;
560 Matrix* M = Matrix_Copy(c->Rays);
561 Vector_Copy(v->p, M->p[i], v->Size);
562 cone * pc = new cone(M);
563 assert (pc->det != 0);
564 if (abs(pc->det) > 1) {
565 assert(abs(pc->det) < abs(c->det));
566 nonuni.push_back(pc);
567 } else {
568 try {
569 handle(pc->poly(), sign(pc->det) * s);
570 delete pc;
571 } catch (...) {
572 delete c;
573 delete pc;
574 while (!nonuni.empty()) {
575 c = nonuni.back();
576 nonuni.pop_back();
577 delete c;
578 }
579 Matrix_Free(M);
580 Vector_Free(v);
581 throw;
582 }
583 }
584 Matrix_Free(M);
585 }
586 Vector_Free(v);
587 delete c;
588 }
589 }
590
591 void polar_decomposer::decompose(Polyhedron *cone, unsigned MaxRays)
592 {
593 Polyhedron_Polarize(cone);
594 if (cone->NbRays - 1 != cone->Dimension) {
595 Polyhedron *tmp = cone;
596 cone = triangularize_cone(cone, MaxRays);
597 Polyhedron_Free(tmp);
598 }
599 try {
600 for (Polyhedron *Polar = cone; Polar; Polar = Polar->next)
601 decomposer::decompose(Polar);
602 Domain_Free(cone);
603 } catch (...) {
604 Domain_Free(cone);
605 throw;
606 }
607 }
608
609 void polar_decomposer::handle(Polyhedron *P, int sign)
610 {
611 Polyhedron_Polarize(P);
612 handle_polar(P, sign);
613 }
614
615 /*
616 * Barvinok's Decomposition of a simplicial cone
617 *
618 * Returns two lists of polyhedra
619 */
620 void barvinok_decompose(Polyhedron *C, Polyhedron **ppos, Polyhedron **pneg)
621 {
622 Polyhedron *pos = *ppos, *neg = *pneg;
623 vector<cone *> nonuni;
624 cone * c = new cone(C);
625 ZZ det = c->det;
626 int s = sign(det);
627 assert(det != 0);
628 if (abs(det) > 1) {
629 nonuni.push_back(c);
630 } else {
631 Polyhedron *p = Polyhedron_Copy(c->Cone);
632 p->next = pos;
633 pos = p;
634 delete c;
635 }
636 vec_ZZ lambda;
637 while (!nonuni.empty()) {
638 c = nonuni.back();
639 nonuni.pop_back();
640 Vector* v = c->short_vector(lambda);
641 for (int i = 0; i < c->Rays->NbRows - 1; ++i) {
642 if (lambda[i] == 0)
643 continue;
644 Matrix* M = Matrix_Copy(c->Rays);
645 Vector_Copy(v->p, M->p[i], v->Size);
646 cone * pc = new cone(M);
647 assert (pc->det != 0);
648 if (abs(pc->det) > 1) {
649 assert(abs(pc->det) < abs(c->det));
650 nonuni.push_back(pc);
651 } else {
652 Polyhedron *p = pc->poly();
653 pc->Cone = 0;
654 if (sign(pc->det) == s) {
655 p->next = pos;
656 pos = p;
657 } else {
658 p->next = neg;
659 neg = p;
660 }
661 delete pc;
662 }
663 Matrix_Free(M);
664 }
665 Vector_Free(v);
666 delete c;
667 }
668 *ppos = pos;
669 *pneg = neg;
670 }
671
672 const int MAX_TRY=10;
673 /*
674 * Searches for a vector that is not orthogonal to any
675 * of the rays in rays.
676 */
677 static void nonorthog(mat_ZZ& rays, vec_ZZ& lambda)
678 {
679 int dim = rays.NumCols();
680 bool found = false;
681 lambda.SetLength(dim);
682 if (dim == 0)
683 return;
684
685 for (int i = 2; !found && i <= 50*dim; i+=4) {
686 for (int j = 0; j < MAX_TRY; ++j) {
687 for (int k = 0; k < dim; ++k) {
688 int r = random_int(i)+2;
689 int v = (2*(r%2)-1) * (r >> 1);
690 lambda[k] = v;
691 }
692 int k = 0;
693 for (; k < rays.NumRows(); ++k)
694 if (lambda * rays[k] == 0)
695 break;
696 if (k == rays.NumRows()) {
697 found = true;
698 break;
699 }
700 }
701 }
702 assert(found);
703 }
704
705 static void randomvector(Polyhedron *P, vec_ZZ& lambda, int nvar)
706 {
707 Value tmp;
708 int max = 10 * 16;
709 unsigned int dim = P->Dimension;
710 value_init(tmp);
711
712 for (int i = 0; i < P->NbRays; ++i) {
713 for (int j = 1; j <= dim; ++j) {
714 value_absolute(tmp, P->Ray[i][j]);
715 int t = VALUE_TO_LONG(tmp) * 16;
716 if (t > max)
717 max = t;
718 }
719 }
720 for (int i = 0; i < P->NbConstraints; ++i) {
721 for (int j = 1; j <= dim; ++j) {
722 value_absolute(tmp, P->Constraint[i][j]);
723 int t = VALUE_TO_LONG(tmp) * 16;
724 if (t > max)
725 max = t;
726 }
727 }
728 value_clear(tmp);
729
730 lambda.SetLength(nvar);
731 for (int k = 0; k < nvar; ++k) {
732 int r = random_int(max*dim)+2;
733 int v = (2*(r%2)-1) * (max/2*dim + (r >> 1));
734 lambda[k] = v;
735 }
736 }
737
738 static void add_rays(mat_ZZ& rays, Polyhedron *i, int *r, int nvar = -1,
739 bool all = false)
740 {
741 unsigned dim = i->Dimension;
742 if (nvar == -1)
743 nvar = dim;
744 for (int k = 0; k < i->NbRays; ++k) {
745 if (!value_zero_p(i->Ray[k][dim+1]))
746 continue;
747 if (!all && nvar != dim && First_Non_Zero(i->Ray[k]+1, nvar) == -1)
748 continue;
749 values2zz(i->Ray[k]+1, rays[(*r)++], nvar);
750 }
751 }
752
753 void lattice_point(Value* values, Polyhedron *i, vec_ZZ& vertex)
754 {
755 unsigned dim = i->Dimension;
756 if(!value_one_p(values[dim])) {
757 Matrix* Rays = rays(i);
758 Matrix *inv = Matrix_Alloc(Rays->NbRows, Rays->NbColumns);
759 int ok = Matrix_Inverse(Rays, inv);
760 assert(ok);
761 Matrix_Free(Rays);
762 Rays = rays(i);
763 Vector *lambda = Vector_Alloc(dim+1);
764 Vector_Matrix_Product(values, inv, lambda->p);
765 Matrix_Free(inv);
766 for (int j = 0; j < dim; ++j)
767 mpz_cdiv_q(lambda->p[j], lambda->p[j], lambda->p[dim]);
768 value_set_si(lambda->p[dim], 1);
769 Vector *A = Vector_Alloc(dim+1);
770 Vector_Matrix_Product(lambda->p, Rays, A->p);
771 Vector_Free(lambda);
772 Matrix_Free(Rays);
773 values2zz(A->p, vertex, dim);
774 Vector_Free(A);
775 } else
776 values2zz(values, vertex, dim);
777 }
778
779 static evalue *term(int param, ZZ& c, Value *den = NULL)
780 {
781 evalue *EP = new evalue();
782 value_init(EP->d);
783 value_set_si(EP->d,0);
784 EP->x.p = new_enode(polynomial, 2, param + 1);
785 evalue_set_si(&EP->x.p->arr[0], 0, 1);
786 value_init(EP->x.p->arr[1].x.n);
787 if (den == NULL)
788 value_set_si(EP->x.p->arr[1].d, 1);
789 else
790 value_assign(EP->x.p->arr[1].d, *den);
791 zz2value(c, EP->x.p->arr[1].x.n);
792 return EP;
793 }
794
795 static void vertex_period(
796 Polyhedron *i, vec_ZZ& lambda, Matrix *T,
797 Value lcm, int p, Vector *val,
798 evalue *E, evalue* ev,
799 ZZ& offset)
800 {
801 unsigned nparam = T->NbRows - 1;
802 unsigned dim = i->Dimension;
803 Value tmp;
804 ZZ nump;
805
806 if (p == nparam) {
807 vec_ZZ vertex;
808 ZZ num, l;
809 Vector * values = Vector_Alloc(dim + 1);
810 Vector_Matrix_Product(val->p, T, values->p);
811 value_assign(values->p[dim], lcm);
812 lattice_point(values->p, i, vertex);
813 num = vertex * lambda;
814 value2zz(lcm, l);
815 num *= l;
816 num += offset;
817 value_init(ev->x.n);
818 zz2value(num, ev->x.n);
819 value_assign(ev->d, lcm);
820 Vector_Free(values);
821 return;
822 }
823
824 value_init(tmp);
825 vec_ZZ vertex;
826 values2zz(T->p[p], vertex, dim);
827 nump = vertex * lambda;
828 if (First_Non_Zero(val->p, p) == -1) {
829 value_assign(tmp, lcm);
830 evalue *ET = term(p, nump, &tmp);
831 eadd(ET, E);
832 free_evalue_refs(ET);
833 delete ET;
834 }
835
836 value_assign(tmp, lcm);
837 if (First_Non_Zero(T->p[p], dim) != -1)
838 Vector_Gcd(T->p[p], dim, &tmp);
839 Gcd(tmp, lcm, &tmp);
840 if (value_lt(tmp, lcm)) {
841 ZZ count;
842
843 value_division(tmp, lcm, tmp);
844 value_set_si(ev->d, 0);
845 ev->x.p = new_enode(periodic, VALUE_TO_INT(tmp), p+1);
846 value2zz(tmp, count);
847 do {
848 value_decrement(tmp, tmp);
849 --count;
850 ZZ new_offset = offset - count * nump;
851 value_assign(val->p[p], tmp);
852 vertex_period(i, lambda, T, lcm, p+1, val, E,
853 &ev->x.p->arr[VALUE_TO_INT(tmp)], new_offset);
854 } while (value_pos_p(tmp));
855 } else
856 vertex_period(i, lambda, T, lcm, p+1, val, E, ev, offset);
857 value_clear(tmp);
858 }
859
860 static void mask_r(Matrix *f, int nr, Vector *lcm, int p, Vector *val, evalue *ev)
861 {
862 unsigned nparam = lcm->Size;
863
864 if (p == nparam) {
865 Vector * prod = Vector_Alloc(f->NbRows);
866 Matrix_Vector_Product(f, val->p, prod->p);
867 int isint = 1;
868 for (int i = 0; i < nr; ++i) {
869 value_modulus(prod->p[i], prod->p[i], f->p[i][nparam+1]);
870 isint &= value_zero_p(prod->p[i]);
871 }
872 value_set_si(ev->d, 1);
873 value_init(ev->x.n);
874 value_set_si(ev->x.n, isint);
875 Vector_Free(prod);
876 return;
877 }
878
879 Value tmp;
880 value_init(tmp);
881 if (value_one_p(lcm->p[p]))
882 mask_r(f, nr, lcm, p+1, val, ev);
883 else {
884 value_assign(tmp, lcm->p[p]);
885 value_set_si(ev->d, 0);
886 ev->x.p = new_enode(periodic, VALUE_TO_INT(tmp), p+1);
887 do {
888 value_decrement(tmp, tmp);
889 value_assign(val->p[p], tmp);
890 mask_r(f, nr, lcm, p+1, val, &ev->x.p->arr[VALUE_TO_INT(tmp)]);
891 } while (value_pos_p(tmp));
892 }
893 value_clear(tmp);
894 }
895
896 static evalue *multi_monom(vec_ZZ& p)
897 {
898 evalue *X = new evalue();
899 value_init(X->d);
900 value_init(X->x.n);
901 unsigned nparam = p.length()-1;
902 zz2value(p[nparam], X->x.n);
903 value_set_si(X->d, 1);
904 for (int i = 0; i < nparam; ++i) {
905 if (p[i] == 0)
906 continue;
907 evalue *T = term(i, p[i]);
908 eadd(T, X);
909 free_evalue_refs(T);
910 delete T;
911 }
912 return X;
913 }
914
915 /*
916 * Check whether mapping polyhedron P on the affine combination
917 * num yields a range that has a fixed quotient on integer
918 * division by d
919 * If zero is true, then we are only interested in the quotient
920 * for the cases where the remainder is zero.
921 * Returns NULL if false and a newly allocated value if true.
922 */
923 static Value *fixed_quotient(Polyhedron *P, vec_ZZ& num, Value d, bool zero)
924 {
925 Value* ret = NULL;
926 int len = num.length();
927 Matrix *T = Matrix_Alloc(2, len);
928 zz2values(num, T->p[0]);
929 value_set_si(T->p[1][len-1], 1);
930 Polyhedron *I = Polyhedron_Image(P, T, P->NbConstraints);
931 Matrix_Free(T);
932
933 int i;
934 for (i = 0; i < I->NbRays; ++i)
935 if (value_zero_p(I->Ray[i][2])) {
936 Polyhedron_Free(I);
937 return NULL;
938 }
939
940 Value min, max;
941 value_init(min);
942 value_init(max);
943 int bounded = line_minmax(I, &min, &max);
944 assert(bounded);
945
946 if (zero)
947 mpz_cdiv_q(min, min, d);
948 else
949 mpz_fdiv_q(min, min, d);
950 mpz_fdiv_q(max, max, d);
951
952 if (value_eq(min, max)) {
953 ALLOC(Value, ret);
954 value_init(*ret);
955 value_assign(*ret, min);
956 }
957 value_clear(min);
958 value_clear(max);
959 return ret;
960 }
961
962 /*
963 * Normalize linear expression coef modulo m
964 * Removes common factor and reduces coefficients
965 * Returns index of first non-zero coefficient or len
966 */
967 static int normal_mod(Value *coef, int len, Value *m)
968 {
969 Value gcd;
970 value_init(gcd);
971
972 Vector_Gcd(coef, len, &gcd);
973 Gcd(gcd, *m, &gcd);
974 Vector_AntiScale(coef, coef, gcd, len);
975
976 value_division(*m, *m, gcd);
977 value_clear(gcd);
978
979 if (value_one_p(*m))
980 return len;
981
982 int j;
983 for (j = 0; j < len; ++j)
984 mpz_fdiv_r(coef[j], coef[j], *m);
985 for (j = 0; j < len; ++j)
986 if (value_notzero_p(coef[j]))
987 break;
988
989 return j;
990 }
991
992 #ifdef USE_MODULO
993 static void mask(Matrix *f, evalue *factor)
994 {
995 int nr = f->NbRows, nc = f->NbColumns;
996 int n;
997 bool found = false;
998 for (n = 0; n < nr && value_notzero_p(f->p[n][nc-1]); ++n)
999 if (value_notone_p(f->p[n][nc-1]) &&
1000 value_notmone_p(f->p[n][nc-1]))
1001 found = true;
1002 if (!found)
1003 return;
1004
1005 evalue EP;
1006 nr = n;
1007
1008 Value m;
1009 value_init(m);
1010
1011 evalue EV;
1012 value_init(EV.d);
1013 value_init(EV.x.n);
1014 value_set_si(EV.x.n, 1);
1015
1016 for (n = 0; n < nr; ++n) {
1017 value_assign(m, f->p[n][nc-1]);
1018 if (value_one_p(m) || value_mone_p(m))
1019 continue;
1020
1021 int j = normal_mod(f->p[n], nc-1, &m);
1022 if (j == nc-1) {
1023 free_evalue_refs(factor);
1024 value_init(factor->d);
1025 evalue_set_si(factor, 0, 1);
1026 break;
1027 }
1028 vec_ZZ row;
1029 values2zz(f->p[n], row, nc-1);
1030 ZZ g;
1031 value2zz(m, g);
1032 if (j < (nc-1)-1 && row[j] > g/2) {
1033 for (int k = j; k < (nc-1); ++k)
1034 if (row[k] != 0)
1035 row[k] = g - row[k];
1036 }
1037
1038 value_init(EP.d);
1039 value_set_si(EP.d, 0);
1040 EP.x.p = new_enode(relation, 2, 0);
1041 value_clear(EP.x.p->arr[1].d);
1042 EP.x.p->arr[1] = *factor;
1043 evalue *ev = &EP.x.p->arr[0];
1044 value_set_si(ev->d, 0);
1045 ev->x.p = new_enode(fractional, 3, -1);
1046 evalue_set_si(&ev->x.p->arr[1], 0, 1);
1047 evalue_set_si(&ev->x.p->arr[2], 1, 1);
1048 evalue *E = multi_monom(row);
1049 value_assign(EV.d, m);
1050 emul(&EV, E);
1051 value_clear(ev->x.p->arr[0].d);
1052 ev->x.p->arr[0] = *E;
1053 delete E;
1054 *factor = EP;
1055 }
1056
1057 value_clear(m);
1058 free_evalue_refs(&EV);
1059 }
1060 #else
1061 /*
1062 *
1063 */
1064 static void mask(Matrix *f, evalue *factor)
1065 {
1066 int nr = f->NbRows, nc = f->NbColumns;
1067 int n;
1068 bool found = false;
1069 for (n = 0; n < nr && value_notzero_p(f->p[n][nc-1]); ++n)
1070 if (value_notone_p(f->p[n][nc-1]) &&
1071 value_notmone_p(f->p[n][nc-1]))
1072 found = true;
1073 if (!found)
1074 return;
1075
1076 Value tmp;
1077 value_init(tmp);
1078 nr = n;
1079 unsigned np = nc - 2;
1080 Vector *lcm = Vector_Alloc(np);
1081 Vector *val = Vector_Alloc(nc);
1082 Vector_Set(val->p, 0, nc);
1083 value_set_si(val->p[np], 1);
1084 Vector_Set(lcm->p, 1, np);
1085 for (n = 0; n < nr; ++n) {
1086 if (value_one_p(f->p[n][nc-1]) ||
1087 value_mone_p(f->p[n][nc-1]))
1088 continue;
1089 for (int j = 0; j < np; ++j)
1090 if (value_notzero_p(f->p[n][j])) {
1091 Gcd(f->p[n][j], f->p[n][nc-1], &tmp);
1092 value_division(tmp, f->p[n][nc-1], tmp);
1093 value_lcm(tmp, lcm->p[j], &lcm->p[j]);
1094 }
1095 }
1096 evalue EP;
1097 value_init(EP.d);
1098 mask_r(f, nr, lcm, 0, val, &EP);
1099 value_clear(tmp);
1100 Vector_Free(val);
1101 Vector_Free(lcm);
1102 emul(&EP,factor);
1103 free_evalue_refs(&EP);
1104 }
1105 #endif
1106
1107 struct term_info {
1108 evalue *E;
1109 ZZ constant;
1110 ZZ coeff;
1111 int pos;
1112 };
1113
1114 static bool mod_needed(Polyhedron *PD, vec_ZZ& num, Value d, evalue *E)
1115 {
1116 Value *q = fixed_quotient(PD, num, d, false);
1117
1118 if (!q)
1119 return true;
1120
1121 value_oppose(*q, *q);
1122 evalue EV;
1123 value_init(EV.d);
1124 value_set_si(EV.d, 1);
1125 value_init(EV.x.n);
1126 value_multiply(EV.x.n, *q, d);
1127 eadd(&EV, E);
1128 free_evalue_refs(&EV);
1129 value_clear(*q);
1130 free(q);
1131 return false;
1132 }
1133
1134 /* modifies f argument ! */
1135 static void ceil_mod(Value *coef, int len, Value d, ZZ& f, evalue *EP, Polyhedron *PD)
1136 {
1137 Value m;
1138 value_init(m);
1139 value_set_si(m, -1);
1140
1141 Vector_Scale(coef, coef, m, len);
1142
1143 value_assign(m, d);
1144 int j = normal_mod(coef, len, &m);
1145
1146 if (j == len) {
1147 value_clear(m);
1148 return;
1149 }
1150
1151 vec_ZZ num;
1152 values2zz(coef, num, len);
1153
1154 ZZ g;
1155 value2zz(m, g);
1156
1157 evalue tmp;
1158 value_init(tmp.d);
1159 evalue_set_si(&tmp, 0, 1);
1160
1161 int p = j;
1162 if (g % 2 == 0)
1163 while (j < len-1 && (num[j] == g/2 || num[j] == 0))
1164 ++j;
1165 if ((j < len-1 && num[j] > g/2) || (j == len-1 && num[j] >= (g+1)/2)) {
1166 for (int k = j; k < len-1; ++k)
1167 if (num[k] != 0)
1168 num[k] = g - num[k];
1169 num[len-1] = g - 1 - num[len-1];
1170 value_assign(tmp.d, m);
1171 ZZ t = f*(g-1);
1172 zz2value(t, tmp.x.n);
1173 eadd(&tmp, EP);
1174 f = -f;
1175 }
1176
1177 if (p >= len-1) {
1178 ZZ t = num[len-1] * f;
1179 zz2value(t, tmp.x.n);
1180 value_assign(tmp.d, m);
1181 eadd(&tmp, EP);
1182 } else {
1183 evalue *E = multi_monom(num);
1184 evalue EV;
1185 value_init(EV.d);
1186
1187 if (PD && !mod_needed(PD, num, m, E)) {
1188 value_init(EV.x.n);
1189 zz2value(f, EV.x.n);
1190 value_assign(EV.d, m);
1191 emul(&EV, E);
1192 eadd(E, EP);
1193 } else {
1194 value_init(EV.x.n);
1195 value_set_si(EV.x.n, 1);
1196 value_assign(EV.d, m);
1197 emul(&EV, E);
1198 value_clear(EV.x.n);
1199 value_set_si(EV.d, 0);
1200 EV.x.p = new_enode(fractional, 3, -1);
1201 evalue_copy(&EV.x.p->arr[0], E);
1202 evalue_set_si(&EV.x.p->arr[1], 0, 1);
1203 value_init(EV.x.p->arr[2].x.n);
1204 zz2value(f, EV.x.p->arr[2].x.n);
1205 value_set_si(EV.x.p->arr[2].d, 1);
1206
1207 eadd(&EV, EP);
1208 }
1209
1210 free_evalue_refs(&EV);
1211 free_evalue_refs(E);
1212 delete E;
1213 }
1214
1215 free_evalue_refs(&tmp);
1216
1217 out:
1218 value_clear(m);
1219 }
1220
1221 #ifdef USE_MODULO
1222 static void ceil(Value *coef, int len, Value d, ZZ& f,
1223 evalue *EP, Polyhedron *PD) {
1224 ceil_mod(coef, len, d, f, EP, PD);
1225 }
1226 #else
1227 static void ceil(Value *coef, int len, Value d, ZZ& f,
1228 evalue *EP, Polyhedron *PD) {
1229 ceil_mod(coef, len, d, f, EP, PD);
1230 evalue_mod2table(EP, len-1);
1231 }
1232 #endif
1233
1234 evalue* bv_ceil3(Value *coef, int len, Value d, Polyhedron *P)
1235 {
1236 Vector *val = Vector_Alloc(len);
1237
1238 Value t;
1239 value_init(t);
1240 value_set_si(t, -1);
1241 Vector_Scale(coef, val->p, t, len);
1242 value_absolute(t, d);
1243
1244 vec_ZZ num;
1245 values2zz(val->p, num, len);
1246 evalue *EP = multi_monom(num);
1247
1248 evalue tmp;
1249 value_init(tmp.d);
1250 value_init(tmp.x.n);
1251 value_set_si(tmp.x.n, 1);
1252 value_assign(tmp.d, t);
1253
1254 emul(&tmp, EP);
1255
1256 ZZ one;
1257 one = 1;
1258 ceil_mod(val->p, len, t, one, EP, P);
1259 value_clear(t);
1260
1261 /* copy EP to malloc'ed evalue */
1262 evalue *E;
1263 ALLOC(evalue, E);
1264 *E = *EP;
1265 delete EP;
1266
1267 free_evalue_refs(&tmp);
1268 Vector_Free(val);
1269
1270 return E;
1271 }
1272
1273 #ifdef USE_MODULO
1274 evalue* lattice_point(
1275 Polyhedron *i, vec_ZZ& lambda, Matrix *W, Value lcm, Polyhedron *PD)
1276 {
1277 unsigned nparam = W->NbColumns - 1;
1278
1279 Matrix* Rays = rays2(i);
1280 Matrix *T = Transpose(Rays);
1281 Matrix *T2 = Matrix_Copy(T);
1282 Matrix *inv = Matrix_Alloc(T2->NbRows, T2->NbColumns);
1283 int ok = Matrix_Inverse(T2, inv);
1284 assert(ok);
1285 Matrix_Free(Rays);
1286 Matrix_Free(T2);
1287 mat_ZZ vertex;
1288 matrix2zz(W, vertex, W->NbRows, W->NbColumns);
1289
1290 vec_ZZ num;
1291 num = lambda * vertex;
1292
1293 evalue *EP = multi_monom(num);
1294
1295 evalue tmp;
1296 value_init(tmp.d);
1297 value_init(tmp.x.n);
1298 value_set_si(tmp.x.n, 1);
1299 value_assign(tmp.d, lcm);
1300
1301 emul(&tmp, EP);
1302
1303 Matrix *L = Matrix_Alloc(inv->NbRows, W->NbColumns);
1304 Matrix_Product(inv, W, L);
1305
1306 mat_ZZ RT;
1307 matrix2zz(T, RT, T->NbRows, T->NbColumns);
1308 Matrix_Free(T);
1309
1310 vec_ZZ p = lambda * RT;
1311
1312 for (int i = 0; i < L->NbRows; ++i) {
1313 ceil_mod(L->p[i], nparam+1, lcm, p[i], EP, PD);
1314 }
1315
1316 Matrix_Free(L);
1317
1318 Matrix_Free(inv);
1319 free_evalue_refs(&tmp);
1320 return EP;
1321 }
1322 #else
1323 evalue* lattice_point(
1324 Polyhedron *i, vec_ZZ& lambda, Matrix *W, Value lcm, Polyhedron *PD)
1325 {
1326 Matrix *T = Transpose(W);
1327 unsigned nparam = T->NbRows - 1;
1328
1329 evalue *EP = new evalue();
1330 value_init(EP->d);
1331 evalue_set_si(EP, 0, 1);
1332
1333 evalue ev;
1334 Vector *val = Vector_Alloc(nparam+1);
1335 value_set_si(val->p[nparam], 1);
1336 ZZ offset(INIT_VAL, 0);
1337 value_init(ev.d);
1338 vertex_period(i, lambda, T, lcm, 0, val, EP, &ev, offset);
1339 Vector_Free(val);
1340 eadd(&ev, EP);
1341 free_evalue_refs(&ev);
1342
1343 Matrix_Free(T);
1344
1345 reduce_evalue(EP);
1346
1347 return EP;
1348 }
1349 #endif
1350
1351 void lattice_point(
1352 Param_Vertices* V, Polyhedron *i, vec_ZZ& lambda, term_info* term,
1353 Polyhedron *PD)
1354 {
1355 unsigned nparam = V->Vertex->NbColumns - 2;
1356 unsigned dim = i->Dimension;
1357 mat_ZZ vertex;
1358 vertex.SetDims(V->Vertex->NbRows, nparam+1);
1359 Value lcm, tmp;
1360 value_init(lcm);
1361 value_init(tmp);
1362 value_set_si(lcm, 1);
1363 for (int j = 0; j < V->Vertex->NbRows; ++j) {
1364 value_lcm(lcm, V->Vertex->p[j][nparam+1], &lcm);
1365 }
1366 if (value_notone_p(lcm)) {
1367 Matrix * mv = Matrix_Alloc(dim, nparam+1);
1368 for (int j = 0 ; j < dim; ++j) {
1369 value_division(tmp, lcm, V->Vertex->p[j][nparam+1]);
1370 Vector_Scale(V->Vertex->p[j], mv->p[j], tmp, nparam+1);
1371 }
1372
1373 term->E = lattice_point(i, lambda, mv, lcm, PD);
1374 term->constant = 0;
1375
1376 Matrix_Free(mv);
1377 value_clear(lcm);
1378 value_clear(tmp);
1379 return;
1380 }
1381 for (int i = 0; i < V->Vertex->NbRows; ++i) {
1382 assert(value_one_p(V->Vertex->p[i][nparam+1])); // for now
1383 values2zz(V->Vertex->p[i], vertex[i], nparam+1);
1384 }
1385
1386 vec_ZZ num;
1387 num = lambda * vertex;
1388
1389 int p = -1;
1390 int nn = 0;
1391 for (int j = 0; j < nparam; ++j)
1392 if (num[j] != 0) {
1393 ++nn;
1394 p = j;
1395 }
1396 if (nn >= 2) {
1397 term->E = multi_monom(num);
1398 term->constant = 0;
1399 } else {
1400 term->E = NULL;
1401 term->constant = num[nparam];
1402 term->pos = p;
1403 if (p != -1)
1404 term->coeff = num[p];
1405 }
1406
1407 value_clear(lcm);
1408 value_clear(tmp);
1409 }
1410
1411 static void normalize(ZZ& sign, ZZ& num, vec_ZZ& den)
1412 {
1413 unsigned dim = den.length();
1414
1415 int change = 0;
1416
1417 for (int j = 0; j < den.length(); ++j) {
1418 if (den[j] > 0)
1419 change ^= 1;
1420 else {
1421 den[j] = abs(den[j]);
1422 num += den[j];
1423 }
1424 }
1425 if (change)
1426 sign = -sign;
1427 }
1428
1429 /* input:
1430 * f: the powers in the denominator for the remaining vars
1431 * each row refers to a factor
1432 * den_s: for each factor, the power of (s+1)
1433 * sign
1434 * num_s: powers in the numerator corresponding to the summed vars
1435 * num_p: powers in the numerator corresponding to the remaining vars
1436 * number of rays in cone: "dim" = "k"
1437 * length of each ray: "dim" = "d"
1438 * for now, it is assumed: k == d
1439 * output:
1440 * den_p: for each factor
1441 * 0: independent of remaining vars
1442 * 1: power corresponds to corresponding row in f
1443 *
1444 * all inputs are subject to change
1445 */
1446 static void normalize(ZZ& sign,
1447 ZZ& num_s, vec_ZZ& num_p, vec_ZZ& den_s, vec_ZZ& den_p,
1448 mat_ZZ& f)
1449 {
1450 unsigned dim = f.NumRows();
1451 unsigned nparam = num_p.length();
1452 unsigned nvar = dim - nparam;
1453
1454 int change = 0;
1455
1456 for (int j = 0; j < den_s.length(); ++j) {
1457 if (den_s[j] == 0) {
1458 den_p[j] = 1;
1459 continue;
1460 }
1461 int k;
1462 for (k = 0; k < nparam; ++k)
1463 if (f[j][k] != 0)
1464 break;
1465 if (k < nparam) {
1466 den_p[j] = 1;
1467 if (den_s[j] > 0) {
1468 f[j] = -f[j];
1469 num_p += f[j];
1470 }
1471 } else
1472 den_p[j] = 0;
1473 if (den_s[j] > 0)
1474 change ^= 1;
1475 else {
1476 den_s[j] = abs(den_s[j]);
1477 num_s += den_s[j];
1478 }
1479 }
1480
1481 if (change)
1482 sign = -sign;
1483 }
1484
1485 struct counter : public polar_decomposer {
1486 vec_ZZ lambda;
1487 mat_ZZ rays;
1488 vec_ZZ vertex;
1489 vec_ZZ den;
1490 ZZ sign;
1491 ZZ num;
1492 int j;
1493 Polyhedron *P;
1494 unsigned dim;
1495 mpq_t count;
1496
1497 counter(Polyhedron *P) {
1498 this->P = P;
1499 dim = P->Dimension;
1500 rays.SetDims(dim, dim);
1501 den.SetLength(dim);
1502 mpq_init(count);
1503 }
1504
1505 void start(unsigned MaxRays);
1506
1507 ~counter() {
1508 mpq_clear(count);
1509 }
1510
1511 virtual void handle_polar(Polyhedron *P, int sign);
1512 };
1513
1514 struct OrthogonalException {} Orthogonal;
1515
1516 void counter::handle_polar(Polyhedron *C, int s)
1517 {
1518 int r = 0;
1519 assert(C->NbRays-1 == dim);
1520 add_rays(rays, C, &r);
1521 for (int k = 0; k < dim; ++k) {
1522 if (lambda * rays[k] == 0)
1523 throw Orthogonal;
1524 }
1525
1526 sign = s;
1527
1528 lattice_point(P->Ray[j]+1, C, vertex);
1529 num = vertex * lambda;
1530 den = rays * lambda;
1531 normalize(sign, num, den);
1532
1533 dpoly d(dim, num);
1534 dpoly n(dim, den[0], 1);
1535 for (int k = 1; k < dim; ++k) {
1536 dpoly fact(dim, den[k], 1);
1537 n *= fact;
1538 }
1539 d.div(n, count, sign);
1540 }
1541
1542 void counter::start(unsigned MaxRays)
1543 {
1544 for (;;) {
1545 try {
1546 randomvector(P, lambda, dim);
1547 for (j = 0; j < P->NbRays; ++j) {
1548 Polyhedron *C = supporting_cone(P, j);
1549 decompose(C, MaxRays);
1550 }
1551 break;
1552 } catch (OrthogonalException &e) {
1553 mpq_set_si(count, 0, 0);
1554 }
1555 }
1556 }
1557
1558 struct reducer : public polar_decomposer {
1559 vec_ZZ vertex;
1560 //vec_ZZ den;
1561 ZZ sgn;
1562 ZZ num;
1563 ZZ one;
1564 int j;
1565 Polyhedron *P;
1566 unsigned dim;
1567 mpq_t tcount;
1568 mpz_t tn;
1569 mpz_t td;
1570 int lower; // call base when only this many variables is left
1571
1572 reducer(Polyhedron *P) {
1573 this->P = P;
1574 dim = P->Dimension;
1575 //den.SetLength(dim);
1576 mpq_init(tcount);
1577 mpz_init(tn);
1578 mpz_init(td);
1579 one = 1;
1580 }
1581
1582 void start(unsigned MaxRays);
1583
1584 ~reducer() {
1585 mpq_clear(tcount);
1586 mpz_clear(tn);
1587 mpz_clear(td);
1588 }
1589
1590 virtual void handle_polar(Polyhedron *P, int sign);
1591 void reduce(ZZ c, ZZ cd, vec_ZZ& num, mat_ZZ& den_f);
1592 virtual void base(ZZ& c, ZZ& cd, vec_ZZ& num, mat_ZZ& den_f) = 0;
1593 virtual void split(vec_ZZ& num, ZZ& num_s, vec_ZZ& num_p,
1594 mat_ZZ& den_f, vec_ZZ& den_s, mat_ZZ& den_r) = 0;
1595 };
1596
1597 void reducer::reduce(ZZ c, ZZ cd, vec_ZZ& num, mat_ZZ& den_f)
1598 {
1599 unsigned len = den_f.NumRows(); // number of factors in den
1600
1601 if (num.length() == lower) {
1602 base(c, cd, num, den_f);
1603 return;
1604 }
1605 assert(num.length() > 1);
1606
1607 vec_ZZ den_s;
1608 mat_ZZ den_r;
1609 ZZ num_s;
1610 vec_ZZ num_p;
1611
1612 split(num, num_s, num_p, den_f, den_s, den_r);
1613
1614 vec_ZZ den_p;
1615 den_p.SetLength(len);
1616
1617 normalize(c, num_s, num_p, den_s, den_p, den_r);
1618
1619 int only_param = 0; // k-r-s from text
1620 int no_param = 0; // r from text
1621 for (int k = 0; k < len; ++k) {
1622 if (den_p[k] == 0)
1623 ++no_param;
1624 else if (den_s[k] == 0)
1625 ++only_param;
1626 }
1627 if (no_param == 0) {
1628 reduce(c, cd, num_p, den_r);
1629 } else {
1630 int k, l;
1631 mat_ZZ pden;
1632 pden.SetDims(only_param, den_r.NumCols());
1633
1634 for (k = 0, l = 0; k < len; ++k)
1635 if (den_s[k] == 0)
1636 pden[l++] = den_r[k];
1637
1638 for (k = 0; k < len; ++k)
1639 if (den_p[k] == 0)
1640 break;
1641
1642 dpoly n(no_param, num_s);
1643 dpoly D(no_param, den_s[k], 1);
1644 for ( ; ++k < len; )
1645 if (den_p[k] == 0) {
1646 dpoly fact(no_param, den_s[k], 1);
1647 D *= fact;
1648 }
1649
1650 if (no_param + only_param == len) {
1651 mpq_set_si(tcount, 0, 1);
1652 n.div(D, tcount, one);
1653
1654 ZZ qn, qd;
1655 value2zz(mpq_numref(tcount), qn);
1656 value2zz(mpq_denref(tcount), qd);
1657
1658 qn *= c;
1659 qd *= cd;
1660
1661 if (qn != 0)
1662 reduce(qn, qd, num_p, pden);
1663 } else {
1664 dpoly_r * r = 0;
1665
1666 for (k = 0; k < len; ++k) {
1667 if (den_s[k] == 0 || den_p[k] == 0)
1668 continue;
1669
1670 dpoly pd(no_param-1, den_s[k], 1);
1671
1672 int l;
1673 for (l = 0; l < k; ++l)
1674 if (den_r[l] == den_r[k])
1675 break;
1676
1677 if (r == 0)
1678 r = new dpoly_r(n, pd, l, len);
1679 else {
1680 dpoly_r *nr = new dpoly_r(r, pd, l, len);
1681 delete r;
1682 r = nr;
1683 }
1684 }
1685
1686 dpoly_r *rc = r->div(D);
1687
1688 rc->denom *= cd;
1689
1690 int common = pden.NumRows();
1691 vector< dpoly_r_term * >& final = rc->c[rc->len-1];
1692 int rows;
1693 for (int j = 0; j < final.size(); ++j) {
1694 if (final[j]->coeff == 0)
1695 continue;
1696 rows = common;
1697 pden.SetDims(rows, pden.NumCols());
1698 for (int k = 0; k < rc->dim; ++k) {
1699 int n = final[j]->powers[k];
1700 if (n == 0)
1701 continue;
1702 pden.SetDims(rows+n, pden.NumCols());
1703 for (int l = 0; l < n; ++l)
1704 pden[rows+l] = den_r[k];
1705 rows += n;
1706 }
1707 final[j]->coeff *= c;
1708 reduce(final[j]->coeff, rc->denom, num_p, pden);
1709 }
1710
1711 delete rc;
1712 delete r;
1713 }
1714 }
1715 }
1716
1717 void reducer::handle_polar(Polyhedron *C, int s)
1718 {
1719 assert(C->NbRays-1 == dim);
1720
1721 sgn = s;
1722
1723 lattice_point(P->Ray[j]+1, C, vertex);
1724
1725 mat_ZZ den;
1726 den.SetDims(dim, dim);
1727
1728 int r;
1729 for (r = 0; r < dim; ++r)
1730 values2zz(C->Ray[r]+1, den[r], dim);
1731
1732 reduce(sgn, one, vertex, den);
1733 }
1734
1735 void reducer::start(unsigned MaxRays)
1736 {
1737 for (j = 0; j < P->NbRays; ++j) {
1738 Polyhedron *C = supporting_cone(P, j);
1739 decompose(C, MaxRays);
1740 }
1741 }
1742
1743 struct ireducer : public reducer {
1744 ireducer(Polyhedron *P) : reducer(P) {}
1745
1746 virtual void split(vec_ZZ& num, ZZ& num_s, vec_ZZ& num_p,
1747 mat_ZZ& den_f, vec_ZZ& den_s, mat_ZZ& den_r) {
1748 unsigned len = den_f.NumRows(); // number of factors in den
1749 unsigned d = num.length() - 1;
1750
1751 den_s.SetLength(len);
1752 den_r.SetDims(len, d);
1753
1754 for (int r = 0; r < len; ++r) {
1755 den_s[r] = den_f[r][0];
1756 for (int k = 1; k <= d; ++k)
1757 den_r[r][k-1] = den_f[r][k];
1758 }
1759
1760 num_s = num[0];
1761 num_p.SetLength(d);
1762 for (int k = 1 ; k <= d; ++k)
1763 num_p[k-1] = num[k];
1764 }
1765 };
1766
1767 // incremental counter
1768 struct icounter : public ireducer {
1769 mpq_t count;
1770
1771 icounter(Polyhedron *P) : ireducer(P) {
1772 mpq_init(count);
1773 lower = 1;
1774 }
1775 ~icounter() {
1776 mpq_clear(count);
1777 }
1778 virtual void base(ZZ& c, ZZ& cd, vec_ZZ& num, mat_ZZ& den_f);
1779 };
1780
1781 void icounter::base(ZZ& c, ZZ& cd, vec_ZZ& num, mat_ZZ& den_f)
1782 {
1783 int r;
1784 unsigned len = den_f.NumRows(); // number of factors in den
1785 vec_ZZ den_s;
1786 den_s.SetLength(len);
1787 ZZ num_s = num[0];
1788 for (r = 0; r < len; ++r)
1789 den_s[r] = den_f[r][0];
1790 normalize(c, num_s, den_s);
1791
1792 dpoly n(len, num_s);
1793 dpoly D(len, den_s[0], 1);
1794 for (int k = 1; k < len; ++k) {
1795 dpoly fact(len, den_s[k], 1);
1796 D *= fact;
1797 }
1798 mpq_set_si(tcount, 0, 1);
1799 n.div(D, tcount, one);
1800 zz2value(c, tn);
1801 zz2value(cd, td);
1802 mpz_mul(mpq_numref(tcount), mpq_numref(tcount), tn);
1803 mpz_mul(mpq_denref(tcount), mpq_denref(tcount), td);
1804 mpq_canonicalize(tcount);
1805 mpq_add(count, count, tcount);
1806 }
1807
1808 struct partial_ireducer : public ireducer {
1809 gen_fun * gf;
1810
1811 partial_ireducer(Polyhedron *P, unsigned nparam) : ireducer(P) {
1812 gf = new gen_fun(Polyhedron_Project(P, nparam));
1813 lower = nparam;
1814 }
1815 ~partial_ireducer() {
1816 }
1817 virtual void base(ZZ& c, ZZ& cd, vec_ZZ& num, mat_ZZ& den_f);
1818 void start(unsigned MaxRays);
1819 };
1820
1821 void partial_ireducer::base(ZZ& c, ZZ& cd, vec_ZZ& num, mat_ZZ& den_f)
1822 {
1823 gf->add(c, cd, num, den_f);
1824 }
1825
1826 void partial_ireducer::start(unsigned MaxRays)
1827 {
1828 for (j = 0; j < P->NbRays; ++j) {
1829 if (!value_pos_p(P->Ray[j][dim+1]))
1830 continue;
1831
1832 Polyhedron *C = supporting_cone(P, j);
1833 decompose(C, MaxRays);
1834 }
1835 }
1836
1837 struct partial_reducer : public reducer {
1838 gen_fun * gf;
1839 vec_ZZ lambda;
1840 vec_ZZ tmp;
1841
1842 partial_reducer(Polyhedron *P, unsigned nparam) : reducer(P) {
1843 gf = new gen_fun(Polyhedron_Project(P, nparam));
1844 lower = nparam;
1845
1846 tmp.SetLength(dim - nparam);
1847 randomvector(P, lambda, dim - nparam);
1848 }
1849 ~partial_reducer() {
1850 }
1851 virtual void base(ZZ& c, ZZ& cd, vec_ZZ& num, mat_ZZ& den_f);
1852 void start(unsigned MaxRays);
1853
1854 virtual void split(vec_ZZ& num, ZZ& num_s, vec_ZZ& num_p,
1855 mat_ZZ& den_f, vec_ZZ& den_s, mat_ZZ& den_r) {
1856 unsigned len = den_f.NumRows(); // number of factors in den
1857 unsigned nvar = tmp.length();
1858
1859 den_s.SetLength(len);
1860 den_r.SetDims(len, lower);
1861
1862 for (int r = 0; r < len; ++r) {
1863 for (int k = 0; k < nvar; ++k)
1864 tmp[k] = den_f[r][k];
1865 den_s[r] = tmp * lambda;
1866
1867 for (int k = nvar; k < dim; ++k)
1868 den_r[r][k-nvar] = den_f[r][k];
1869 }
1870
1871 for (int k = 0; k < nvar; ++k)
1872 tmp[k] = num[k];
1873 num_s = tmp *lambda;
1874 num_p.SetLength(lower);
1875 for (int k = nvar ; k < dim; ++k)
1876 num_p[k-nvar] = num[k];
1877 }
1878 };
1879
1880 void partial_reducer::base(ZZ& c, ZZ& cd, vec_ZZ& num, mat_ZZ& den_f)
1881 {
1882 gf->add(c, cd, num, den_f);
1883 }
1884
1885 void partial_reducer::start(unsigned MaxRays)
1886 {
1887 for (j = 0; j < P->NbRays; ++j) {
1888 if (!value_pos_p(P->Ray[j][dim+1]))
1889 continue;
1890
1891 Polyhedron *C = supporting_cone(P, j);
1892 decompose(C, MaxRays);
1893 }
1894 }
1895
1896 struct bfc_term_base {
1897 // the number of times a given factor appears in the denominator
1898 int *powers;
1899 mat_ZZ terms;
1900
1901 bfc_term_base(int len) {
1902 powers = new int[len];
1903 }
1904
1905 virtual ~bfc_term_base() {
1906 delete [] powers;
1907 }
1908 };
1909
1910 struct bfc_term : public bfc_term_base {
1911 vec_ZZ cn;
1912 vec_ZZ cd;
1913
1914 bfc_term(int len) : bfc_term_base(len) {}
1915 };
1916
1917 struct bfe_term : public bfc_term_base {
1918 vector<evalue *> factors;
1919
1920 bfe_term(int len) : bfc_term_base(len) {
1921 }
1922
1923 ~bfe_term() {
1924 for (int i = 0; i < factors.size(); ++i) {
1925 if (!factors[i])
1926 continue;
1927 free_evalue_refs(factors[i]);
1928 delete factors[i];
1929 }
1930 }
1931 };
1932
1933 typedef vector< bfc_term_base * > bfc_vec;
1934
1935 struct bf_reducer;
1936
1937 struct bf_base : public virtual polar_decomposer {
1938 Polyhedron *P;
1939 unsigned dim;
1940 int j;
1941 ZZ one;
1942 mpq_t tcount;
1943 mpz_t tn;
1944 mpz_t td;
1945 int lower; // call base when only this many variables is left
1946
1947 bf_base(Polyhedron *P, unsigned dim) {
1948 this->P = P;
1949 this->dim = dim;
1950 mpq_init(tcount);
1951 mpz_init(tn);
1952 mpz_init(td);
1953 one = 1;
1954 }
1955
1956 ~bf_base() {
1957 mpq_clear(tcount);
1958 mpz_clear(tn);
1959 mpz_clear(td);
1960 }
1961
1962 void start(unsigned MaxRays);
1963 virtual void handle_polar(Polyhedron *P, int sign);
1964 int setup_factors(Polyhedron *P, mat_ZZ& factors, bfc_term_base* t, int s);
1965
1966 bfc_term_base* find_bfc_term(bfc_vec& v, int *powers, int len);
1967 void add_term(bfc_term_base *t, vec_ZZ& num1, vec_ZZ& num);
1968 void add_term(bfc_term_base *t, vec_ZZ& num);
1969
1970 void reduce(mat_ZZ& factors, bfc_vec& v);
1971 virtual void base(mat_ZZ& factors, bfc_vec& v) = 0;
1972
1973 virtual bfc_term_base* new_bf_term(int len) = 0;
1974 virtual void set_factor(bfc_term_base *t, int k, int change) = 0;
1975 virtual void set_factor(bfc_term_base *t, int k, mpq_t &f, int change) = 0;
1976 virtual void set_factor(bfc_term_base *t, int k, ZZ& n, ZZ& d, int change) = 0;
1977 virtual void update_term(bfc_term_base *t, int i) = 0;
1978 virtual void insert_term(bfc_term_base *t, int i) = 0;
1979 virtual bool constant_vertex(int dim) = 0;
1980 virtual void cum(bf_reducer *bfr, bfc_term_base *t, int k,
1981 dpoly_r *r) = 0;
1982 };
1983
1984 static int lex_cmp(vec_ZZ& a, vec_ZZ& b)
1985 {
1986 assert(a.length() == b.length());
1987
1988 for (int j = 0; j < a.length(); ++j)
1989 if (a[j] != b[j])
1990 return a[j] < b[j] ? -1 : 1;
1991 return 0;
1992 }
1993
1994 void bf_base::add_term(bfc_term_base *t, vec_ZZ& num_orig, vec_ZZ& extra_num)
1995 {
1996 vec_ZZ num;
1997 int d = num_orig.length();
1998 num.SetLength(d-1);
1999 for (int l = 0; l < d-1; ++l)
2000 num[l] = num_orig[l+1] + extra_num[l];
2001
2002 add_term(t, num);
2003 }
2004
2005 void bf_base::add_term(bfc_term_base *t, vec_ZZ& num)
2006 {
2007 int len = t->terms.NumRows();
2008 int i, r;
2009 for (i = 0; i < len; ++i) {
2010 r = lex_cmp(t->terms[i], num);
2011 if (r >= 0)
2012 break;
2013 }
2014 if (i == len || r > 0) {
2015 t->terms.SetDims(len+1, num.length());
2016 insert_term(t, i);
2017 t->terms[i] = num;
2018 } else {
2019 // i < len && r == 0
2020 update_term(t, i);
2021 }
2022 }
2023
2024 static void print_int_vector(int *v, int len, char *name)
2025 {
2026 cerr << name << endl;
2027 for (int j = 0; j < len; ++j) {
2028 cerr << v[j] << " ";
2029 }
2030 cerr << endl;
2031 }
2032
2033 static void print_bfc_terms(mat_ZZ& factors, bfc_vec& v)
2034 {
2035 cerr << endl;
2036 cerr << "factors" << endl;
2037 cerr << factors << endl;
2038 for (int i = 0; i < v.size(); ++i) {
2039 cerr << "term: " << i << endl;
2040 print_int_vector(v[i]->powers, factors.NumRows(), "powers");
2041 cerr << "terms" << endl;
2042 cerr << v[i]->terms << endl;
2043 bfc_term* bfct = static_cast<bfc_term *>(v[i]);
2044 cerr << bfct->cn << endl;
2045 cerr << bfct->cd << endl;
2046 }
2047 }
2048
2049 static void print_bfe_terms(mat_ZZ& factors, bfc_vec& v)
2050 {
2051 cerr << endl;
2052 cerr << "factors" << endl;
2053 cerr << factors << endl;
2054 for (int i = 0; i < v.size(); ++i) {
2055 cerr << "term: " << i << endl;
2056 print_int_vector(v[i]->powers, factors.NumRows(), "powers");
2057 cerr << "terms" << endl;
2058 cerr << v[i]->terms << endl;
2059 bfe_term* bfet = static_cast<bfe_term *>(v[i]);
2060 for (int j = 0; j < v[i]->terms.NumRows(); ++j) {
2061 char * test[] = {"a", "b"};
2062 print_evalue(stderr, bfet->factors[j], test);
2063 fprintf(stderr, "\n");
2064 }
2065 }
2066 }
2067
2068 bfc_term_base* bf_base::find_bfc_term(bfc_vec& v, int *powers, int len)
2069 {
2070 bfc_vec::iterator i;
2071 for (i = v.begin(); i != v.end(); ++i) {
2072 int j;
2073 for (j = 0; j < len; ++j)
2074 if ((*i)->powers[j] != powers[j])
2075 break;
2076 if (j == len)
2077 return (*i);
2078 if ((*i)->powers[j] > powers[j])
2079 break;
2080 }
2081
2082 bfc_term_base* t = new_bf_term(len);
2083 v.insert(i, t);
2084 memcpy(t->powers, powers, len * sizeof(int));
2085
2086 return t;
2087 }
2088
2089 struct bf_reducer {
2090 mat_ZZ& factors;
2091 bfc_vec& v;
2092 bf_base *bf;
2093
2094 unsigned nf;
2095 unsigned d;
2096
2097 mat_ZZ nfactors;
2098 int *old2new;
2099 int *sign;
2100 unsigned int nnf;
2101 bfc_vec vn;
2102
2103 vec_ZZ extra_num;
2104 int changes;
2105 int no_param; // r from text
2106 int only_param; // k-r-s from text
2107 int total_power; // k from text
2108
2109 // created in compute_reduced_factors
2110 int *bpowers;
2111 // set in update_powers
2112 int *npowers;
2113 vec_ZZ l_extra_num;
2114 int l_changes;
2115
2116 bf_reducer(mat_ZZ& factors, bfc_vec& v, bf_base *bf)
2117 : factors(factors), v(v), bf(bf) {
2118 nf = factors.NumRows();
2119 d = factors.NumCols();
2120 old2new = new int[nf];
2121 sign = new int[nf];
2122
2123 extra_num.SetLength(d-1);
2124 }
2125 ~bf_reducer() {
2126 delete [] old2new;
2127 delete [] sign;
2128 delete [] npowers;
2129 delete [] bpowers;
2130 }
2131
2132 void compute_reduced_factors();
2133 void compute_extra_num(int i);
2134
2135 void reduce();
2136
2137 void update_powers(int *powers, int len);
2138 };
2139
2140 void bf_reducer::compute_extra_num(int i)
2141 {
2142 clear(extra_num);
2143 changes = 0;
2144 no_param = 0; // r from text
2145 only_param = 0; // k-r-s from text
2146 total_power = 0; // k from text
2147
2148 for (int j = 0; j < nf; ++j) {
2149 if (v[i]->powers[j] == 0)
2150 continue;
2151
2152 total_power += v[i]->powers[j];
2153 if (factors[j][0] == 0) {
2154 only_param += v[i]->powers[j];
2155 continue;
2156 }
2157
2158 if (old2new[j] == -1)
2159 no_param += v[i]->powers[j];
2160 else
2161 extra_num += -sign[j] * v[i]->powers[j] * nfactors[old2new[j]];
2162 changes += v[i]->powers[j];
2163 }
2164 }
2165
2166 void bf_reducer::update_powers(int *powers, int len)
2167 {
2168 for (int l = 0; l < nnf; ++l)
2169 npowers[l] = bpowers[l];
2170
2171 l_extra_num = extra_num;
2172 l_changes = changes;
2173
2174 for (int l = 0; l < len; ++l) {
2175 int n = powers[l];
2176 if (n == 0)
2177 continue;
2178 assert(old2new[l] != -1);
2179
2180 npowers[old2new[l]] += n;
2181 // interpretation of sign has been inverted
2182 // since we inverted the power for specialization
2183 if (sign[l] == 1) {
2184 l_extra_num += n * nfactors[old2new[l]];
2185 l_changes += n;
2186 }
2187 }
2188 }
2189
2190
2191 void bf_reducer::compute_reduced_factors()
2192 {
2193 unsigned nf = factors.NumRows();
2194 unsigned d = factors.NumCols();
2195 nnf = 0;
2196 nfactors.SetDims(nnf, d-1);
2197
2198 for (int i = 0; i < nf; ++i) {
2199 int j;
2200 int s = 1;
2201 for (j = 0; j < nnf; ++j) {
2202 int k;
2203 for (k = 1; k < d; ++k)
2204 if (factors[i][k] != 0 || nfactors[j][k-1] != 0)
2205 break;
2206 if (k < d && factors[i][k] == -nfactors[j][k-1])
2207 s = -1;
2208 for (; k < d; ++k)
2209 if (factors[i][k] != s * nfactors[j][k-1])
2210 break;
2211 if (k == d)
2212 break;
2213 }
2214 old2new[i] = j;
2215 if (j == nnf) {
2216 int k;
2217 for (k = 1; k < d; ++k)
2218 if (factors[i][k] != 0)
2219 break;
2220 if (k < d) {
2221 if (factors[i][k] < 0)
2222 s = -1;
2223 nfactors.SetDims(++nnf, d-1);
2224 for (int k = 1; k < d; ++k)
2225 nfactors[j][k-1] = s * factors[i][k];
2226 } else
2227 old2new[i] = -1;
2228 }
2229 sign[i] = s;
2230 }
2231 npowers = new int[nnf];
2232 bpowers = new int[nnf];
2233 }
2234
2235 void bf_reducer::reduce()
2236 {
2237 compute_reduced_factors();
2238
2239 for (int i = 0; i < v.size(); ++i) {
2240 compute_extra_num(i);
2241
2242 if (no_param == 0) {
2243 vec_ZZ extra_num;
2244 extra_num.SetLength(d-1);
2245 int changes = 0;
2246 int npowers[nnf];
2247 for (int k = 0; k < nnf; ++k)
2248 npowers[k] = 0;
2249 for (int k = 0; k < nf; ++k) {
2250 assert(old2new[k] != -1);
2251 npowers[old2new[k]] += v[i]->powers[k];
2252 if (sign[k] == -1) {
2253 extra_num += v[i]->powers[k] * nfactors[old2new[k]];
2254 changes += v[i]->powers[k];
2255 }
2256 }
2257
2258 bfc_term_base * t = bf->find_bfc_term(vn, npowers, nnf);
2259 for (int k = 0; k < v[i]->terms.NumRows(); ++k) {
2260 bf->set_factor(v[i], k, changes % 2);
2261 bf->add_term(t, v[i]->terms[k], extra_num);
2262 }
2263 } else {
2264 // powers of "constant" part
2265 for (int k = 0; k < nnf; ++k)
2266 bpowers[k] = 0;
2267 for (int k = 0; k < nf; ++k) {
2268 if (factors[k][0] != 0)
2269 continue;
2270 assert(old2new[k] != -1);
2271 bpowers[old2new[k]] += v[i]->powers[k];
2272 if (sign[k] == -1) {
2273 extra_num += v[i]->powers[k] * nfactors[old2new[k]];
2274 changes += v[i]->powers[k];
2275 }
2276 }
2277
2278 int j;
2279 for (j = 0; j < nf; ++j)
2280 if (old2new[j] == -1 && v[i]->powers[j] > 0)
2281 break;
2282
2283 dpoly D(no_param, factors[j][0], 1);
2284 for (int k = 1; k < v[i]->powers[j]; ++k) {
2285 dpoly fact(no_param, factors[j][0], 1);
2286 D *= fact;
2287 }
2288 for ( ; ++j < nf; )
2289 if (old2new[j] == -1)
2290 for (int k = 0; k < v[i]->powers[j]; ++k) {
2291 dpoly fact(no_param, factors[j][0], 1);
2292 D *= fact;
2293 }
2294
2295 if (no_param + only_param == total_power &&
2296 bf->constant_vertex(d)) {
2297 bfc_term_base * t = NULL;
2298 vec_ZZ num;
2299 num.SetLength(d-1);
2300 ZZ cn;
2301 ZZ cd;
2302 for (int k = 0; k < v[i]->terms.NumRows(); ++k) {
2303 dpoly n(no_param, v[i]->terms[k][0]);
2304 mpq_set_si(bf->tcount, 0, 1);
2305 n.div(D, bf->tcount, bf->one);
2306
2307 if (value_zero_p(mpq_numref(bf->tcount)))
2308 continue;
2309
2310 if (!t)
2311 t = bf->find_bfc_term(vn, bpowers, nnf);
2312 bf->set_factor(v[i], k, bf->tcount, changes % 2);
2313 bf->add_term(t, v[i]->terms[k], extra_num);
2314 }
2315 } else {
2316 for (int j = 0; j < v[i]->terms.NumRows(); ++j) {
2317 dpoly n(no_param, v[i]->terms[j][0]);
2318
2319 dpoly_r * r = 0;
2320 if (no_param + only_param == total_power)
2321 r = new dpoly_r(n, nf);
2322 else
2323 for (int k = 0; k < nf; ++k) {
2324 if (v[i]->powers[k] == 0)
2325 continue;
2326 if (factors[k][0] == 0 || old2new[k] == -1)
2327 continue;
2328
2329 dpoly pd(no_param-1, factors[k][0], 1);
2330
2331 for (int l = 0; l < v[i]->powers[k]; ++l) {
2332 int q;
2333 for (q = 0; q < k; ++q)
2334 if (old2new[q] == old2new[k] &&
2335 sign[q] == sign[k])
2336 break;
2337
2338 if (r == 0)
2339 r = new dpoly_r(n, pd, q, nf);
2340 else {
2341 dpoly_r *nr = new dpoly_r(r, pd, q, nf);
2342 delete r;
2343 r = nr;
2344 }
2345 }
2346 }
2347
2348 dpoly_r *rc = r->div(D);
2349 delete r;
2350
2351 if (bf->constant_vertex(d)) {
2352 vector< dpoly_r_term * >& final = rc->c[rc->len-1];
2353
2354 for (int k = 0; k < final.size(); ++k) {
2355 if (final[k]->coeff == 0)
2356 continue;
2357
2358 update_powers(final[k]->powers, rc->dim);
2359
2360 bfc_term_base * t = bf->find_bfc_term(vn, npowers, nnf);
2361 bf->set_factor(v[i], j, final[k]->coeff, rc->denom, l_changes % 2);
2362 bf->add_term(t, v[i]->terms[j], l_extra_num);
2363 }
2364 } else
2365 bf->cum(this, v[i], j, rc);
2366
2367 delete rc;
2368 }
2369 }
2370 }
2371 delete v[i];
2372 }
2373
2374 }
2375
2376 void bf_base::reduce(mat_ZZ& factors, bfc_vec& v)
2377 {
2378 assert(v.size() > 0);
2379 unsigned nf = factors.NumRows();
2380 unsigned d = factors.NumCols();
2381
2382 if (d == lower)
2383 return base(factors, v);
2384
2385 bf_reducer bfr(factors, v, this);
2386
2387 bfr.reduce();
2388
2389 if (bfr.vn.size() > 0)
2390 reduce(bfr.nfactors, bfr.vn);
2391 }
2392
2393 int bf_base::setup_factors(Polyhedron *C, mat_ZZ& factors,
2394 bfc_term_base* t, int s)
2395 {
2396 factors.SetDims(dim, dim);
2397
2398 int r;
2399
2400 for (r = 0; r < dim; ++r)
2401 t->powers[r] = 1;
2402
2403 for (r = 0; r < dim; ++r) {
2404 values2zz(C->Ray[r]+1, factors[r], dim);
2405 int k;
2406 for (k = 0; k < dim; ++k)
2407 if (factors[r][k] != 0)
2408 break;
2409 if (factors[r][k] < 0) {
2410 factors[r] = -factors[r];
2411 t->terms[0] += factors[r];
2412 s = -s;
2413 }
2414 }
2415
2416 return s;
2417 }
2418
2419 void bf_base::handle_polar(Polyhedron *C, int s)
2420 {
2421 bfc_term* t = new bfc_term(dim);
2422 vector< bfc_term_base * > v;
2423 v.push_back(t);
2424
2425 assert(C->NbRays-1 == dim);
2426
2427 t->cn.SetLength(1);
2428 t->cd.SetLength(1);
2429
2430 t->terms.SetDims(1, dim);
2431 lattice_point(P->Ray[j]+1, C, t->terms[0]);
2432
2433 // the elements of factors are always lexpositive
2434 mat_ZZ factors;
2435 s = setup_factors(C, factors, t, s);
2436
2437 t->cn[0] = s;
2438 t->cd[0] = 1;
2439
2440 reduce(factors, v);
2441 }
2442
2443 void bf_base::start(unsigned MaxRays)
2444 {
2445 for (j = 0; j < P->NbRays; ++j) {
2446 Polyhedron *C = supporting_cone(P, j);
2447 decompose(C, MaxRays);
2448 }
2449 }
2450
2451 struct bfcounter_base : public bf_base {
2452 ZZ cn;
2453 ZZ cd;
2454
2455 bfcounter_base(Polyhedron *P) : bf_base(P, P->Dimension) {
2456 }
2457
2458 bfc_term_base* new_bf_term(int len) {
2459 bfc_term* t = new bfc_term(len);
2460 t->cn.SetLength(0);
2461 t->cd.SetLength(0);
2462 return t;
2463 }
2464
2465 virtual void set_factor(bfc_term_base *t, int k, int change) {
2466 bfc_term* bfct = static_cast<bfc_term *>(t);
2467 cn = bfct->cn[k];
2468 if (change)
2469 cn = -cn;
2470 cd = bfct->cd[k];
2471 }
2472
2473 virtual void set_factor(bfc_term_base *t, int k, mpq_t &f, int change) {
2474 bfc_term* bfct = static_cast<bfc_term *>(t);
2475 value2zz(mpq_numref(f), cn);
2476 value2zz(mpq_denref(f), cd);
2477 cn *= bfct->cn[k];
2478 if (change)
2479 cn = -cn;
2480 cd *= bfct->cd[k];
2481 }
2482
2483 virtual void set_factor(bfc_term_base *t, int k, ZZ& n, ZZ& d, int change) {
2484 bfc_term* bfct = static_cast<bfc_term *>(t);
2485 cn = bfct->cn[k] * n;
2486 if (change)
2487 cn = -cn;
2488 cd = bfct->cd[k] * d;
2489 }
2490
2491 virtual void insert_term(bfc_term_base *t, int i) {
2492 bfc_term* bfct = static_cast<bfc_term *>(t);
2493 int len = t->terms.NumRows()-1; // already increased by one
2494
2495 bfct->cn.SetLength(len+1);
2496 bfct->cd.SetLength(len+1);
2497 for (int j = len; j > i; --j) {
2498 bfct->cn[j] = bfct->cn[j-1];
2499 bfct->cd[j] = bfct->cd[j-1];
2500 t->terms[j] = t->terms[j-1];
2501 }
2502 bfct->cn[i] = cn;
2503 bfct->cd[i] = cd;
2504 }
2505
2506 virtual void update_term(bfc_term_base *t, int i) {
2507 bfc_term* bfct = static_cast<bfc_term *>(t);
2508
2509 ZZ g = GCD(bfct->cd[i], cd);
2510 ZZ n = cn * bfct->cd[i]/g + bfct->cn[i] * cd/g;
2511 ZZ d = bfct->cd[i] * cd / g;
2512 bfct->cn[i] = n;
2513 bfct->cd[i] = d;
2514 }
2515
2516 virtual bool constant_vertex(int dim) { return true; }
2517 virtual void cum(bf_reducer *bfr, bfc_term_base *t, int k, dpoly_r *r) {
2518 assert(0);
2519 }
2520 };
2521
2522 struct bfcounter : public bfcounter_base {
2523 mpq_t count;
2524
2525 bfcounter(Polyhedron *P) : bfcounter_base(P) {
2526 mpq_init(count);
2527 lower = 1;
2528 }
2529 ~bfcounter() {
2530 mpq_clear(count);
2531 }
2532 virtual void base(mat_ZZ& factors, bfc_vec& v);
2533 };
2534
2535 void bfcounter::base(mat_ZZ& factors, bfc_vec& v)
2536 {
2537 unsigned nf = factors.NumRows();
2538
2539 for (int i = 0; i < v.size(); ++i) {
2540 bfc_term* bfct = static_cast<bfc_term *>(v[i]);
2541 int total_power = 0;
2542 // factor is always positive, so we always
2543 // change signs
2544 for (int k = 0; k < nf; ++k)
2545 total_power += v[i]->powers[k];
2546
2547 int j;
2548 for (j = 0; j < nf; ++j)
2549 if (v[i]->powers[j] > 0)
2550 break;
2551
2552 dpoly D(total_power, factors[j][0], 1);
2553 for (int k = 1; k < v[i]->powers[j]; ++k) {
2554 dpoly fact(total_power, factors[j][0], 1);
2555 D *= fact;
2556 }
2557 for ( ; ++j < nf; )
2558 for (int k = 0; k < v[i]->powers[j]; ++k) {
2559 dpoly fact(total_power, factors[j][0], 1);
2560 D *= fact;
2561 }
2562
2563 for (int k = 0; k < v[i]->terms.NumRows(); ++k) {
2564 dpoly n(total_power, v[i]->terms[k][0]);
2565 mpq_set_si(tcount, 0, 1);
2566 n.div(D, tcount, one);
2567 if (total_power % 2)
2568 bfct->cn[k] = -bfct->cn[k];
2569 zz2value(bfct->cn[k], tn);
2570 zz2value(bfct->cd[k], td);
2571
2572 mpz_mul(mpq_numref(tcount), mpq_numref(tcount), tn);
2573 mpz_mul(mpq_denref(tcount), mpq_denref(tcount), td);
2574 mpq_canonicalize(tcount);
2575 mpq_add(count, count, tcount);
2576 }
2577 delete v[i];
2578 }
2579 }
2580
2581 struct partial_bfcounter : public bfcounter_base {
2582 gen_fun * gf;
2583
2584 partial_bfcounter(Polyhedron *P, unsigned nparam) : bfcounter_base(P) {
2585 gf = new gen_fun(Polyhedron_Project(P, nparam));
2586 lower = nparam;
2587 }
2588 ~partial_bfcounter() {
2589 }
2590 virtual void base(mat_ZZ& factors, bfc_vec& v);
2591 void start(unsigned MaxRays);
2592 };
2593
2594 void partial_bfcounter::base(mat_ZZ& factors, bfc_vec& v)
2595 {
2596 mat_ZZ den;
2597 unsigned nf = factors.NumRows();
2598
2599 for (int i = 0; i < v.size(); ++i) {
2600 bfc_term* bfct = static_cast<bfc_term *>(v[i]);
2601 den.SetDims(0, lower);
2602 int total_power = 0;
2603 int p = 0;
2604 for (int j = 0; j < nf; ++j) {
2605 total_power += v[i]->powers[j];
2606 den.SetDims(total_power, lower);
2607 for (int k = 0; k < v[i]->powers[j]; ++k)
2608 den[p++] = factors[j];
2609 }
2610 for (int j = 0; j < v[i]->terms.NumRows(); ++j)
2611 gf->add(bfct->cn[j], bfct->cd[j], v[i]->terms[j], den);
2612 delete v[i];
2613 }
2614 }
2615
2616 void partial_bfcounter::start(unsigned MaxRays)
2617 {
2618 for (j = 0; j < P->NbRays; ++j) {
2619 if (!value_pos_p(P->Ray[j][dim+1]))
2620 continue;
2621
2622 Polyhedron *C = supporting_cone(P, j);
2623 decompose(C, MaxRays);
2624 }
2625 }
2626
2627
2628 typedef Polyhedron * Polyhedron_p;
2629
2630 static void barvinok_count_f(Polyhedron *P, Value* result, unsigned NbMaxCons);
2631
2632 void barvinok_count(Polyhedron *P, Value* result, unsigned NbMaxCons)
2633 {
2634 unsigned dim;
2635 int allocated = 0;
2636 Polyhedron *Q;
2637 int r = 0;
2638
2639 if (emptyQ2(P)) {
2640 value_set_si(*result, 0);
2641 return;
2642 }
2643 if (P->NbBid == 0)
2644 for (; r < P->NbRays; ++r)
2645 if (value_zero_p(P->Ray[r][P->Dimension+1]))
2646 break;
2647 if (P->NbBid !=0 || r < P->NbRays) {
2648 value_set_si(*result, -1);
2649 return;
2650 }
2651 if (P->NbEq != 0) {
2652 P = remove_equalities(P);
2653 if (emptyQ(P)) {
2654 Polyhedron_Free(P);
2655 value_set_si(*result, 0);
2656 return;
2657 }
2658 allocated = 1;
2659 }
2660 if (P->Dimension == 0) {
2661 /* Test whether the constraints are satisfied */
2662 POL_ENSURE_VERTICES(P);
2663 value_set_si(*result, !emptyQ(P));
2664 if (allocated)
2665 Polyhedron_Free(P);
2666 return;
2667 }
2668 Q = Polyhedron_Factor(P, 0, NbMaxCons);
2669 if (Q) {
2670 if (allocated)
2671 Polyhedron_Free(P);
2672 P = Q;
2673 allocated = 1;
2674 }
2675
2676 barvinok_count_f(P, result, NbMaxCons);
2677 if (Q && P->next) {
2678 Value factor;
2679 value_init(factor);
2680
2681 for (Q = P->next; Q; Q = Q->next) {
2682 barvinok_count_f(Q, &factor, NbMaxCons);
2683 value_multiply(*result, *result, factor);
2684 }
2685
2686 value_clear(factor);
2687 }
2688
2689 if (allocated)
2690 Polyhedron_Free(P);
2691 }
2692
2693 static void barvinok_count_f(Polyhedron *P, Value* result, unsigned NbMaxCons)
2694 {
2695 if (P->Dimension == 1)
2696 return Line_Length(P, result);
2697
2698 int c = P->NbConstraints;
2699 POL_ENSURE_FACETS(P);
2700 if (c != P->NbConstraints || P->NbEq != 0)
2701 return barvinok_count(P, result, NbMaxCons);
2702
2703 POL_ENSURE_VERTICES(P);
2704
2705 #ifdef USE_INCREMENTAL_BF
2706 bfcounter cnt(P);
2707 #elif defined USE_INCREMENTAL_DF
2708 icounter cnt(P);
2709 #else
2710 counter cnt(P);
2711 #endif
2712 cnt.start(NbMaxCons);
2713
2714 assert(value_one_p(&cnt.count[0]._mp_den));
2715 value_assign(*result, &cnt.count[0]._mp_num);
2716 }
2717
2718 static void uni_polynom(int param, Vector *c, evalue *EP)
2719 {
2720 unsigned dim = c->Size-2;
2721 value_init(EP->d);
2722 value_set_si(EP->d,0);
2723 EP->x.p = new_enode(polynomial, dim+1, param+1);
2724 for (int j = 0; j <= dim; ++j)
2725 evalue_set(&EP->x.p->arr[j], c->p[j], c->p[dim+1]);
2726 }
2727
2728 static void multi_polynom(Vector *c, evalue* X, evalue *EP)
2729 {
2730 unsigned dim = c->Size-2;
2731 evalue EC;
2732
2733 value_init(EC.d);
2734 evalue_set(&EC, c->p[dim], c->p[dim+1]);
2735
2736 value_init(EP->d);
2737 evalue_set(EP, c->p[dim], c->p[dim+1]);
2738
2739 for (int i = dim-1; i >= 0; --i) {
2740 emul(X, EP);
2741 value_assign(EC.x.n, c->p[i]);
2742 eadd(&EC, EP);
2743 }
2744 free_evalue_refs(&EC);
2745 }
2746
2747 Polyhedron *unfringe (Polyhedron *P, unsigned MaxRays)
2748 {
2749 int len = P->Dimension+2;
2750 Polyhedron *T, *R = P;
2751 Value g;
2752 value_init(g);
2753 Vector *row = Vector_Alloc(len);
2754 value_set_si(row->p[0], 1);
2755
2756 R = DomainConstraintSimplify(Polyhedron_Copy(P), MaxRays);
2757
2758 Matrix *M = Matrix_Alloc(2, len-1);
2759 value_set_si(M->p[1][len-2], 1);
2760 for (int v = 0; v < P->Dimension; ++v) {
2761 value_set_si(M->p[0][v], 1);
2762 Polyhedron *I = Polyhedron_Image(P, M, 2+1);
2763 value_set_si(M->p[0][v], 0);
2764 for (int r = 0; r < I->NbConstraints; ++r) {
2765 if (value_zero_p(I->Constraint[r][0]))
2766 continue;
2767 if (value_zero_p(I->Constraint[r][1]))
2768 continue;
2769 if (value_one_p(I->Constraint[r][1]))
2770 continue;
2771 if (value_mone_p(I->Constraint[r][1]))
2772 continue;
2773 value_absolute(g, I->Constraint[r][1]);
2774 Vector_Set(row->p+1, 0, len-2);
2775 value_division(row->p[1+v], I->Constraint[r][1], g);
2776 mpz_fdiv_q(row->p[len-1], I->Constraint[r][2], g);
2777 T = R;
2778 R = AddConstraints(row->p, 1, R, MaxRays);
2779 if (T != P)
2780 Polyhedron_Free(T);
2781 }
2782 Polyhedron_Free(I);
2783 }
2784 Matrix_Free(M);
2785 Vector_Free(row);
2786 value_clear(g);
2787 return R;
2788 }
2789
2790 static Polyhedron *reduce_domain(Polyhedron *D, Matrix *CT, Polyhedron *CEq,
2791 Polyhedron **fVD, int nd, unsigned MaxRays)
2792 {
2793 assert(CEq);
2794
2795 Polyhedron *Dt;
2796 Dt = CT ? DomainPreimage(D, CT, MaxRays) : D;
2797 Polyhedron *rVD = DomainIntersection(Dt, CEq, MaxRays);
2798
2799 /* if rVD is empty or too small in geometric dimension */
2800 if(!rVD || emptyQ(rVD) ||
2801 (rVD->Dimension-rVD->NbEq < Dt->Dimension-Dt->NbEq-CEq->NbEq)) {
2802 if(rVD)
2803 Domain_Free(rVD);
2804 if (CT)
2805 Domain_Free(Dt);
2806 return 0; /* empty validity domain */
2807 }
2808
2809 if (CT)
2810 Domain_Free(Dt);
2811
2812 fVD[nd] = Domain_Copy(rVD);
2813 for (int i = 0 ; i < nd; ++i) {
2814 Polyhedron *I = DomainIntersection(fVD[nd], fVD[i], MaxRays);
2815 if (emptyQ(I)) {
2816 Domain_Free(I);
2817 continue;
2818 }
2819 Polyhedron *F = DomainSimplify(I, fVD[nd], MaxRays);
2820 if (F->NbEq == 1) {
2821 Polyhedron *T = rVD;
2822 rVD = DomainDifference(rVD, F, MaxRays);
2823 Domain_Free(T);
2824 }
2825 Domain_Free(F);
2826 Domain_Free(I);
2827 }
2828
2829 rVD = DomainConstraintSimplify(rVD, MaxRays);
2830 if (emptyQ(rVD)) {
2831 Domain_Free(fVD[nd]);
2832 Domain_Free(rVD);
2833 return 0;
2834 }
2835
2836 Value c;
2837 value_init(c);
2838 barvinok_count(rVD, &c, MaxRays);
2839 if (value_zero_p(c)) {
2840 Domain_Free(rVD);
2841 rVD = 0;
2842 }
2843 value_clear(c);
2844
2845 return rVD;
2846 }
2847
2848 static bool Polyhedron_is_infinite(Polyhedron *P, unsigned nparam)
2849 {
2850 int r;
2851 for (r = 0; r < P->NbRays; ++r)
2852 if (value_zero_p(P->Ray[r][0]) ||
2853 value_zero_p(P->Ray[r][P->Dimension+1])) {
2854 int i;
2855 for (i = P->Dimension - nparam; i < P->Dimension; ++i)
2856 if (value_notzero_p(P->Ray[r][i+1]))
2857 break;
2858 if (i >= P->Dimension)
2859 break;
2860 }
2861 return r < P->NbRays;
2862 }
2863
2864 /* Check whether all rays point in the positive directions
2865 * for the parameters
2866 */
2867 static bool Polyhedron_has_positive_rays(Polyhedron *P, unsigned nparam)
2868 {
2869 int r;
2870 for (r = 0; r < P->NbRays; ++r)
2871 if (value_zero_p(P->Ray[r][P->Dimension+1])) {
2872 int i;
2873 for (i = P->Dimension - nparam; i < P->Dimension; ++i)
2874 if (value_neg_p(P->Ray[r][i+1]))
2875 return false;
2876 }
2877 return true;
2878 }
2879
2880 typedef evalue * evalue_p;
2881
2882 struct enumerator : public polar_decomposer {
2883 vec_ZZ lambda;
2884 unsigned dim, nbV;
2885 evalue ** vE;
2886 int _i;
2887 mat_ZZ rays;
2888 vec_ZZ den;
2889 ZZ sign;
2890 Polyhedron *P;
2891 Param_Vertices *V;
2892 term_info num;
2893 Vector *c;
2894 mpq_t count;
2895
2896 enumerator(Polyhedron *P, unsigned dim, unsigned nbV) {
2897 this->P = P;
2898 this->dim = dim;
2899 this->nbV = nbV;
2900 randomvector(P, lambda, dim);
2901 rays.SetDims(dim, dim);
2902 den.SetLength(dim);
2903 c = Vector_Alloc(dim+2);
2904
2905 vE = new evalue_p[nbV];
2906 for (int j = 0; j < nbV; ++j)
2907 vE[j] = 0;
2908
2909 mpq_init(count);
2910 }
2911
2912 void decompose_at(Param_Vertices *V, int _i, unsigned MaxRays) {
2913 Polyhedron *C = supporting_cone_p(P, V);
2914 this->_i = _i;
2915 this->V = V;
2916
2917 vE[_i] = new evalue;
2918 value_init(vE[_i]->d);
2919 evalue_set_si(vE[_i], 0, 1);
2920
2921 decompose(C, MaxRays);
2922 }
2923
2924 ~enumerator() {
2925 mpq_clear(count);
2926 Vector_Free(c);
2927
2928 for (int j = 0; j < nbV; ++j)
2929 if (vE[j]) {
2930 free_evalue_refs(vE[j]);
2931 delete vE[j];
2932 }
2933 delete [] vE;
2934 }
2935
2936 virtual void handle_polar(Polyhedron *P, int sign);
2937 };
2938
2939 void enumerator::handle_polar(Polyhedron *C, int s)
2940 {
2941 int r = 0;
2942 assert(C->NbRays-1 == dim);
2943 add_rays(rays, C, &r);
2944 for (int k = 0; k < dim; ++k) {
2945 if (lambda * rays[k] == 0)
2946 throw Orthogonal;
2947 }
2948
2949 sign = s;
2950
2951 lattice_point(V, C, lambda, &num, 0);
2952 den = rays * lambda;
2953 normalize(sign, num.constant, den);
2954
2955 dpoly n(dim, den[0], 1);
2956 for (int k = 1; k < dim; ++k) {
2957 dpoly fact(dim, den[k], 1);
2958 n *= fact;
2959 }
2960 if (num.E != NULL) {
2961 ZZ one(INIT_VAL, 1);
2962 dpoly_n d(dim, num.constant, one);
2963 d.div(n, c, sign);
2964 evalue EV;
2965 multi_polynom(c, num.E, &EV);
2966 eadd(&EV , vE[_i]);
2967 free_evalue_refs(&EV);
2968 free_evalue_refs(num.E);
2969 delete num.E;
2970 } else if (num.pos != -1) {
2971 dpoly_n d(dim, num.constant, num.coeff);
2972 d.div(n, c, sign);
2973 evalue EV;
2974 uni_polynom(num.pos, c, &EV);
2975 eadd(&EV , vE[_i]);
2976 free_evalue_refs(&EV);
2977 } else {
2978 mpq_set_si(count, 0, 1);
2979 dpoly d(dim, num.constant);
2980 d.div(n, count, sign);
2981 evalue EV;
2982 value_init(EV.d);
2983 evalue_set(&EV, &count[0]._mp_num, &count[0]._mp_den);
2984 eadd(&EV , vE[_i]);
2985 free_evalue_refs(&EV);
2986 }
2987 }
2988
2989 struct enumerator_base : public virtual polar_decomposer {
2990 Polyhedron *P;
2991 unsigned dim, nbV;
2992 evalue ** vE;
2993 int _i;
2994 Param_Vertices *V;
2995 evalue ** E_vertex;
2996 evalue mone;
2997
2998 enumerator_base(Polyhedron *P, unsigned dim, unsigned nbV) {
2999 this->P = P;
3000 this->dim = dim;
3001 this->nbV = nbV;
3002
3003 vE = new evalue_p[nbV];
3004 for (int j = 0; j < nbV; ++j)
3005 vE[j] = 0;
3006
3007 E_vertex = new evalue_p[dim];
3008
3009 value_init(mone.d);
3010 evalue_set_si(&mone, -1, 1);
3011 }
3012
3013 void decompose_at(Param_Vertices *V, int _i, unsigned MaxRays/*, Polyhedron *pVD*/) {
3014 Polyhedron *C = supporting_cone_p(P, V);
3015 this->_i = _i;
3016 this->V = V;
3017 //this->pVD = pVD;
3018
3019 vE[_i] = new evalue;
3020 value_init(vE[_i]->d);
3021 evalue_set_si(vE[_i], 0, 1);
3022
3023 decompose(C, MaxRays);
3024 }
3025
3026 ~enumerator_base() {
3027 for (int j = 0; j < nbV; ++j)
3028 if (vE[j]) {
3029 free_evalue_refs(vE[j]);
3030 delete vE[j];
3031 }
3032 delete [] vE;
3033
3034 delete [] E_vertex;
3035
3036 free_evalue_refs(&mone);
3037 }
3038
3039 evalue *E_num(int i, int d) {
3040 return E_vertex[i + (dim-d)];
3041 }
3042 };
3043
3044 struct cumulator {
3045 evalue *factor;
3046 evalue *v;
3047 dpoly_r *r;
3048
3049 cumulator(evalue *factor, evalue *v, dpoly_r *r) :
3050 factor(factor), v(v), r(r) {}
3051
3052 void cumulate();
3053
3054 virtual void add_term(int *powers, int len, evalue *f2) = 0;
3055 };
3056
3057 void cumulator::cumulate()
3058 {
3059 evalue cum; // factor * 1 * E_num[0]/1 * (E_num[0]-1)/2 *...
3060 evalue f;
3061 evalue t; // E_num[0] - (m-1)
3062 #ifdef USE_MODULO
3063 evalue *cst;
3064 #else
3065 evalue mone;
3066 value_init(mone.d);
3067 evalue_set_si(&mone, -1, 1);
3068 #endif
3069
3070 value_init(cum.d);
3071 evalue_copy(&cum, factor);
3072 value_init(f.d);
3073 value_init(f.x.n);
3074 value_set_si(f.d, 1);
3075 value_set_si(f.x.n, 1);
3076 value_init(t.d);
3077 evalue_copy(&t, v);
3078
3079 #ifdef USE_MODULO
3080 for (cst = &t; value_zero_p(cst->d); ) {
3081 if (cst->x.p->type == fractional)
3082 cst = &cst->x.p->arr[1];
3083 else
3084 cst = &cst->x.p->arr[0];
3085 }
3086 #endif
3087
3088 for (int m = 0; m < r->len; ++m) {
3089 if (m > 0) {
3090 if (m > 1) {
3091 value_set_si(f.d, m);
3092 emul(&f, &cum);
3093 #ifdef USE_MODULO
3094 value_subtract(cst->x.n, cst->x.n, cst->d);
3095 #else
3096 eadd(&mone, &t);
3097 #endif
3098 }
3099 emul(&t, &cum);
3100 }
3101 vector< dpoly_r_term * >& current = r->c[r->len-1-m];
3102 for (int j = 0; j < current.size(); ++j) {
3103 if (current[j]->coeff == 0)
3104 continue;
3105 evalue *f2 = new evalue;
3106 value_init(f2->d);
3107 value_init(f2->x.n);
3108 zz2value(current[j]->coeff, f2->x.n);
3109 zz2value(r->denom, f2->d);
3110 emul(&cum, f2);
3111
3112 add_term(current[j]->powers, r->dim, f2);
3113 }
3114 }
3115 free_evalue_refs(&f);
3116 free_evalue_refs(&t);
3117 free_evalue_refs(&cum);
3118 #ifndef USE_MODULO
3119 free_evalue_refs(&mone);
3120 #endif
3121 }
3122
3123 struct E_poly_term {
3124 int *powers;
3125 evalue *E;
3126 };
3127
3128 struct ie_cum : public cumulator {
3129 vector<E_poly_term *> terms;
3130
3131 ie_cum(evalue *factor, evalue *v, dpoly_r *r) : cumulator(factor, v, r) {}
3132
3133 virtual void add_term(int *powers, int len, evalue *f2);
3134 };
3135
3136 void ie_cum::add_term(int *powers, int len, evalue *f2)
3137 {
3138 int k;
3139 for (k = 0; k < terms.size(); ++k) {
3140 if (memcmp(terms[k]->powers, powers, len * sizeof(int)) == 0) {
3141 eadd(f2, terms[k]->E);
3142 free_evalue_refs(f2);
3143 delete f2;
3144 break;
3145 }
3146 }
3147 if (k >= terms.size()) {
3148 E_poly_term *ET = new E_poly_term;
3149 ET->powers = new int[len];
3150 memcpy(ET->powers, powers, len * sizeof(int));
3151 ET->E = f2;
3152 terms.push_back(ET);
3153 }
3154 }
3155
3156 struct ienumerator : public virtual polar_decomposer, public enumerator_base {
3157 //Polyhedron *pVD;
3158 mat_ZZ den;
3159 vec_ZZ vertex;
3160 mpq_t tcount;
3161
3162 ienumerator(Polyhedron *P, unsigned dim, unsigned nbV) :
3163 enumerator_base(P, dim, nbV) {
3164 vertex.SetLength(dim);
3165
3166 den.SetDims(dim, dim);
3167 mpq_init(tcount);
3168 }
3169
3170 ~ienumerator() {
3171 mpq_clear(tcount);
3172 }
3173
3174 virtual void handle_polar(Polyhedron *P, int sign);
3175 void reduce(evalue *factor, vec_ZZ& num, mat_ZZ& den_f);
3176 };
3177
3178 static evalue* new_zero_ep()
3179 {
3180 evalue *EP;
3181 ALLOC(evalue, EP);
3182 value_init(EP->d);
3183 evalue_set_si(EP, 0, 1);
3184 return EP;
3185 }
3186
3187 void lattice_point(Param_Vertices *V, Polyhedron *C, vec_ZZ& num,
3188 evalue **E_vertex)
3189 {
3190 unsigned nparam = V->Vertex->NbColumns - 2;
3191 unsigned dim = C->Dimension;
3192 vec_ZZ vertex;
3193 vertex.SetLength(nparam+1);
3194
3195 Value lcm, tmp;
3196 value_init(lcm);
3197 value_init(tmp);
3198 value_set_si(lcm, 1);
3199
3200 for (int j = 0; j < V->Vertex->NbRows; ++j) {
3201 value_lcm(lcm, V->Vertex->p[j][nparam+1], &lcm);
3202 }
3203
3204 if (value_notone_p(lcm)) {
3205 Matrix * mv = Matrix_Alloc(dim, nparam+1);
3206 for (int j = 0 ; j < dim; ++j) {
3207 value_division(tmp, lcm, V->Vertex->p[j][nparam+1]);
3208 Vector_Scale(V->Vertex->p[j], mv->p[j], tmp, nparam+1);
3209 }
3210
3211 Matrix* Rays = rays2(C);
3212 Matrix *T = Transpose(Rays);
3213 Matrix *T2 = Matrix_Copy(T);
3214 Matrix *inv = Matrix_Alloc(T2->NbRows, T2->NbColumns);
3215 int ok = Matrix_Inverse(T2, inv);
3216 assert(ok);
3217 Matrix_Free(Rays);
3218 Matrix_Free(T2);
3219 Matrix *L = Matrix_Alloc(inv->NbRows, mv->NbColumns);
3220 Matrix_Product(inv, mv, L);
3221 Matrix_Free(inv);
3222
3223 evalue f;
3224 value_init(f.d);
3225 value_init(f.x.n);
3226
3227 ZZ one;
3228
3229 evalue *remainders[dim];
3230 for (int i = 0; i < dim; ++i) {
3231 remainders[i] = new_zero_ep();
3232 one = 1;
3233 ceil(L->p[i], nparam+1, lcm, one, remainders[i], 0);
3234 }
3235 Matrix_Free(L);
3236
3237
3238 for (int i = 0; i < V->Vertex->NbRows; ++i) {
3239 values2zz(mv->p[i], vertex, nparam+1);
3240 E_vertex[i] = multi_monom(vertex);
3241 num[i] = 0;
3242
3243 value_set_si(f.x.n, 1);
3244 value_assign(f.d, lcm);
3245
3246 emul(&f, E_vertex[i]);
3247
3248 for (int j = 0; j < dim; ++j) {
3249 if (value_zero_p(T->p[i][j]))
3250 continue;
3251 evalue cp;
3252 value_init(cp.d);
3253 evalue_copy(&cp, remainders[j]);
3254 if (value_notone_p(T->p[i][j])) {
3255 value_set_si(f.d, 1);
3256 value_assign(f.x.n, T->p[i][j]);
3257 emul(&f, &cp);
3258 }
3259 eadd(&cp, E_vertex[i]);
3260 free_evalue_refs(&cp);
3261 }
3262 }
3263 for (int i = 0; i < dim; ++i) {
3264 free_evalue_refs(remainders[i]);
3265 free(remainders[i]);
3266 }
3267
3268 free_evalue_refs(&f);
3269
3270 Matrix_Free(T);
3271 Matrix_Free(mv);
3272 value_clear(lcm);
3273 value_clear(tmp);
3274 return;
3275 }
3276 value_clear(lcm);
3277 value_clear(tmp);
3278
3279 for (int i = 0; i < V->Vertex->NbRows; ++i) {
3280 /* fixed value */
3281 if (First_Non_Zero(V->Vertex->p[i], nparam) == -1) {
3282 E_vertex[i] = 0;
3283 value2zz(V->Vertex->p[i][nparam], num[i]);
3284 } else {
3285 values2zz(V->Vertex->p[i], vertex, nparam+1);
3286 E_vertex[i] = multi_monom(vertex);
3287 num[i] = 0;
3288 }
3289 }
3290 }
3291
3292 void ienumerator::reduce(
3293 evalue *factor, vec_ZZ& num, mat_ZZ& den_f)
3294 {
3295 unsigned len = den_f.NumRows(); // number of factors in den
3296 unsigned dim = num.length();
3297
3298 if (dim == 0) {
3299 eadd(factor, vE[_i]);
3300 return;
3301 }
3302
3303 vec_ZZ den_s;
3304 den_s.SetLength(len);
3305 mat_ZZ den_r;
3306 den_r.SetDims(len, dim-1);
3307
3308 int r, k;
3309
3310 for (r = 0; r < len; ++r) {
3311 den_s[r] = den_f[r][0];
3312 for (k = 0; k <= dim-1; ++k)
3313 if (k != 0)
3314 den_r[r][k-(k>0)] = den_f[r][k];
3315 }
3316
3317 ZZ num_s = num[0];
3318 vec_ZZ num_p;
3319 num_p.SetLength(dim-1);
3320 for (k = 0 ; k <= dim-1; ++k)
3321 if (k != 0)
3322 num_p[k-(k>0)] = num[k];
3323
3324 vec_ZZ den_p;
3325 den_p.SetLength(len);
3326
3327 ZZ one;
3328 one = 1;
3329 normalize(one, num_s, num_p, den_s, den_p, den_r);
3330 if (one != 1)
3331 emul(&mone, factor);
3332
3333 int only_param = 0;
3334 int no_param = 0;
3335 for (int k = 0; k < len; ++k) {
3336 if (den_p[k] == 0)
3337 ++no_param;
3338 else if (den_s[k] == 0)
3339 ++only_param;
3340 }
3341 if (no_param == 0) {
3342 reduce(factor, num_p, den_r);
3343 } else {
3344 int k, l;
3345 mat_ZZ pden;
3346 pden.SetDims(only_param, dim-1);
3347
3348 for (k = 0, l = 0; k < len; ++k)
3349 if (den_s[k] == 0)
3350 pden[l++] = den_r[k];
3351
3352 for (k = 0; k < len; ++k)
3353 if (den_p[k] == 0)
3354 break;
3355
3356 dpoly n(no_param, num_s);
3357 dpoly D(no_param, den_s[k], 1);
3358 for ( ; ++k < len; )
3359 if (den_p[k] == 0) {
3360 dpoly fact(no_param, den_s[k], 1);
3361 D *= fact;
3362 }
3363
3364 dpoly_r * r = 0;
3365 // if no_param + only_param == len then all powers
3366 // below will be all zero
3367 if (no_param + only_param == len) {
3368 if (E_num(0, dim) != 0)
3369 r = new dpoly_r(n, len);
3370 else {
3371 mpq_set_si(tcount, 0, 1);
3372 one = 1;
3373 n.div(D, tcount, one);
3374
3375 if (value_notzero_p(mpq_numref(tcount))) {
3376 evalue f;
3377 value_init(f.d);
3378 value_init(f.x.n);
3379 value_assign(f.x.n, mpq_numref(tcount));
3380 value_assign(f.d, mpq_denref(tcount));
3381 emul(&f, factor);
3382 reduce(factor, num_p, pden);
3383 free_evalue_refs(&f);
3384 }
3385 return;
3386 }
3387 } else {
3388 for (k = 0; k < len; ++k) {
3389 if (den_s[k] == 0 || den_p[k] == 0)
3390 continue;
3391
3392 dpoly pd(no_param-1, den_s[k], 1);
3393
3394 int l;
3395 for (l = 0; l < k; ++l)
3396 if (den_r[l] == den_r[k])
3397 break;
3398
3399 if (r == 0)
3400 r = new dpoly_r(n, pd, l, len);
3401 else {
3402 dpoly_r *nr = new dpoly_r(r, pd, l, len);
3403 delete r;
3404 r = nr;
3405 }
3406 }
3407 }
3408 dpoly_r *rc = r->div(D);
3409 delete r;
3410 r = rc;
3411 if (E_num(0, dim) == 0) {
3412 int common = pden.NumRows();
3413 vector< dpoly_r_term * >& final = r->c[r->len-1];
3414 int rows;
3415 evalue t;
3416 evalue f;
3417 value_init(f.d);
3418 value_init(f.x.n);
3419 zz2value(r->denom, f.d);
3420 for (int j = 0; j < final.size(); ++j) {
3421 if (final[j]->coeff == 0)
3422 continue;
3423 rows = common;
3424 for (int k = 0; k < r->dim; ++k) {
3425 int n = final[j]->powers[k];
3426 if (n == 0)
3427 continue;
3428 pden.SetDims(rows+n, pden.NumCols());
3429 for (int l = 0; l < n; ++l)
3430 pden[rows+l] = den_r[k];
3431 rows += n;
3432 }
3433 value_init(t.d);
3434 evalue_copy(&t, factor);
3435 zz2value(final[j]->coeff, f.x.n);
3436 emul(&f, &t);
3437 reduce(&t, num_p, pden);
3438 free_evalue_refs(&t);
3439 }
3440 free_evalue_refs(&f);
3441 } else {
3442 ie_cum cum(factor, E_num(0, dim), r);
3443 cum.cumulate();
3444
3445 int common = pden.NumRows();
3446 int rows;
3447 for (int j = 0; j < cum.terms.size(); ++j) {
3448 rows = common;
3449 pden.SetDims(rows, pden.NumCols());
3450 for (int k = 0; k < r->dim; ++k) {
3451 int n = cum.terms[j]->powers[k];
3452 if (n == 0)
3453 continue;
3454 pden.SetDims(rows+n, pden.NumCols());
3455 for (int l = 0; l < n; ++l)
3456 pden[rows+l] = den_r[k];
3457 rows += n;
3458 }
3459 reduce(cum.terms[j]->E, num_p, pden);
3460 free_evalue_refs(cum.terms[j]->E);
3461 delete cum.terms[j]->E;
3462 delete [] cum.terms[j]->powers;
3463 delete cum.terms[j];
3464 }
3465 }
3466 delete r;
3467 }
3468 }
3469
3470 static int type_offset(enode *p)
3471 {
3472 return p->type == fractional ? 1 :
3473 p->type == flooring ? 1 : 0;
3474 }
3475
3476 static int edegree(evalue *e)
3477 {
3478 int d = 0;
3479 enode *p;
3480
3481 if (value_notzero_p(e->d))
3482 return 0;
3483
3484 p = e->x.p;
3485 int i = type_offset(p);
3486 if (p->size-i-1 > d)
3487 d = p->size - i - 1;
3488 for (; i < p->size; i++) {
3489 int d2 = edegree(&p->arr[i]);
3490 if (d2 > d)
3491 d = d2;
3492 }
3493 return d;
3494 }
3495
3496 void ienumerator::handle_polar(Polyhedron *C, int s)
3497 {
3498 assert(C->NbRays-1 == dim);
3499
3500 lattice_point(V, C, vertex, E_vertex);
3501
3502 int r;
3503 for (r = 0; r < dim; ++r)
3504 values2zz(C->Ray[r]+1, den[r], dim);
3505
3506 evalue one;
3507 value_init(one.d);
3508 evalue_set_si(&one, s, 1);
3509 reduce(&one, vertex, den);
3510 free_evalue_refs(&one);
3511
3512 for (int i = 0; i < dim; ++i)
3513 if (E_vertex[i]) {
3514 free_evalue_refs(E_vertex[i]);
3515 delete E_vertex[i];
3516 }
3517
3518 /*
3519 {
3520 char * test[] = {"a", "b"};
3521 evalue E;
3522 value_init(E.d);
3523 evalue_copy(&E, vE[_i]);
3524 frac2floor_in_domain(&E, pVD);
3525 printf("***** Curr value:");
3526 print_evalue(stdout, &E, test);
3527 fprintf(stdout, "\n");
3528 }
3529 */
3530
3531 }
3532
3533 struct bfenumerator : public bf_base, public enumerator_base {
3534 evalue *factor;
3535
3536 bfenumerator(Polyhedron *P, unsigned dim, unsigned nbV) :
3537 bf_base(P, dim), enumerator_base(P, dim, nbV) {
3538 lower = 0;
3539 factor = NULL;
3540 }
3541
3542 ~bfenumerator() {
3543 }
3544
3545 virtual void handle_polar(Polyhedron *P, int sign);
3546 virtual void base(mat_ZZ& factors, bfc_vec& v);
3547
3548 bfc_term_base* new_bf_term(int len) {
3549 bfe_term* t = new bfe_term(len);
3550 return t;
3551 }
3552
3553 virtual void set_factor(bfc_term_base *t, int k, int change) {
3554 bfe_term* bfet = static_cast<bfe_term *>(t);
3555 factor = bfet->factors[k];
3556 assert(factor != NULL);
3557 bfet->factors[k] = NULL;
3558 if (change)
3559 emul(&mone, factor);
3560 }
3561
3562 virtual void set_factor(bfc_term_base *t, int k, mpq_t &q, int change) {
3563 bfe_term* bfet = static_cast<bfe_term *>(t);
3564 factor = bfet->factors[k];
3565 assert(factor != NULL);
3566 bfet->factors[k] = NULL;
3567
3568 evalue f;
3569 value_init(f.d);
3570 value_init(f.x.n);
3571 if (change)
3572 value_oppose(f.x.n, mpq_numref(q));
3573 else
3574 value_assign(f.x.n, mpq_numref(q));
3575 value_assign(f.d, mpq_denref(q));
3576 emul(&f, factor);
3577 }
3578
3579 virtual void set_factor(bfc_term_base *t, int k, ZZ& n, ZZ& d, int change) {
3580 bfe_term* bfet = static_cast<bfe_term *>(t);
3581
3582 factor = new evalue;
3583
3584 evalue f;
3585 value_init(f.d);
3586 value_init(f.x.n);
3587 zz2value(n, f.x.n);
3588 if (change)
3589 value_oppose(f.x.n, f.x.n);
3590 zz2value(d, f.d);
3591
3592 value_init(factor->d);
3593 evalue_copy(factor, bfet->factors[k]);
3594 emul(&f, factor);
3595 }
3596
3597 void set_factor(evalue *f, int change) {
3598 if (change)
3599 emul(&mone, f);
3600 factor = f;
3601 }
3602
3603 virtual void insert_term(bfc_term_base *t, int i) {
3604 bfe_term* bfet = static_cast<bfe_term *>(t);
3605 int len = t->terms.NumRows()-1; // already increased by one
3606
3607 bfet->factors.resize(len+1);
3608 for (int j = len; j > i; --j) {
3609 bfet->factors[j] = bfet->factors[j-1];
3610 t->terms[j] = t->terms[j-1];
3611 }
3612 bfet->factors[i] = factor;
3613 factor = NULL;
3614 }
3615
3616 virtual void update_term(bfc_term_base *t, int i) {
3617 bfe_term* bfet = static_cast<bfe_term *>(t);
3618
3619 eadd(factor, bfet->factors[i]);
3620 free_evalue_refs(factor);
3621 delete factor;
3622 }
3623
3624 virtual bool constant_vertex(int dim) { return E_num(0, dim) == 0; }
3625
3626 virtual void cum(bf_reducer *bfr, bfc_term_base *t, int k, dpoly_r *r);
3627 };
3628
3629 struct bfe_cum : public cumulator {
3630 bfenumerator *bfe;
3631 bfc_term_base *told;
3632 int k;
3633 bf_reducer *bfr;
3634
3635 bfe_cum(evalue *factor, evalue *v, dpoly_r *r, bf_reducer *bfr,
3636 bfc_term_base *t, int k, bfenumerator *e) :
3637 cumulator(factor, v, r), told(t), k(k),
3638 bfr(bfr), bfe(e) {
3639 }
3640
3641 virtual void add_term(int *powers, int len, evalue *f2);
3642 };
3643
3644 void bfe_cum::add_term(int *powers, int len, evalue *f2)
3645 {
3646 bfr->update_powers(powers, len);
3647
3648 bfc_term_base * t = bfe->find_bfc_term(bfr->vn, bfr->npowers, bfr->nnf);
3649 bfe->set_factor(f2, bfr->l_changes % 2);
3650 bfe->add_term(t, told->terms[k], bfr->l_extra_num);
3651 }
3652
3653 void bfenumerator::cum(bf_reducer *bfr, bfc_term_base *t, int k,
3654 dpoly_r *r)
3655 {
3656 bfe_term* bfet = static_cast<bfe_term *>(t);
3657 bfe_cum cum(bfet->factors[k], E_num(0, bfr->d), r, bfr, t, k, this);
3658 cum.cumulate();
3659 }
3660
3661 void bfenumerator::base(mat_ZZ& factors, bfc_vec& v)
3662 {
3663 for (int i = 0; i < v.size(); ++i) {
3664 assert(v[i]->terms.NumRows() == 1);
3665 evalue *factor = static_cast<bfe_term *>(v[i])->factors[0];
3666 eadd(factor, vE[_i]);
3667 delete v[i];
3668 }
3669 }
3670
3671 void bfenumerator::handle_polar(Polyhedron *C, int s)
3672 {
3673 assert(C->NbRays-1 == enumerator_base::dim);
3674
3675 bfe_term* t = new bfe_term(enumerator_base::dim);
3676 vector< bfc_term_base * > v;
3677 v.push_back(t);
3678
3679 t->factors.resize(1);
3680
3681 t->terms.SetDims(1, enumerator_base::dim);
3682 lattice_point(V, C, t->terms[0], E_vertex);
3683
3684 // the elements of factors are always lexpositive
3685 mat_ZZ factors;
3686 s = setup_factors(C, factors, t, s);
3687
3688 t->factors[0] = new evalue;
3689 value_init(t->factors[0]->d);
3690 evalue_set_si(t->factors[0], s, 1);
3691 reduce(factors, v);
3692
3693 for (int i = 0; i < enumerator_base::dim; ++i)
3694 if (E_vertex[i]) {
3695 free_evalue_refs(E_vertex[i]);
3696 delete E_vertex[i];
3697 }
3698 }
3699
3700 #ifdef HAVE_CORRECT_VERTICES
3701 static inline Param_Polyhedron *Polyhedron2Param_SD(Polyhedron **Din,
3702 Polyhedron *Cin,int WS,Polyhedron **CEq,Matrix **CT)
3703 {
3704 return Polyhedron2Param_SimplifiedDomain(Din, Cin, WS, CEq, CT);
3705 }
3706 #else
3707 static Param_Polyhedron *Polyhedron2Param_SD(Polyhedron **Din,
3708 Polyhedron *Cin,int WS,Polyhedron **CEq,Matrix **CT)
3709 {
3710 static char data[] = " 1 0 0 0 0 1 -18 "
3711 " 1 0 0 -20 0 19 1 "
3712 " 1 0 1 20 0 -20 16 "
3713 " 1 0 0 0 0 -1 19 "
3714 " 1 0 -1 0 0 0 4 "
3715 " 1 4 -20 0 0 -1 23 "
3716 " 1 -4 20 0 0 1 -22 "
3717 " 1 0 1 0 20 -20 16 "
3718 " 1 0 0 0 -20 19 1 ";
3719 static int checked = 0;
3720 if (!checked) {
3721 checked = 1;
3722 char *p = data;
3723 int n, v, i;
3724 Matrix *M = Matrix_Alloc(9, 7);
3725 for (i = 0; i < 9; ++i)
3726 for (int j = 0; j < 7; ++j) {
3727 sscanf(p, "%d%n", &v, &n);
3728 p += n;
3729 value_set_si(M->p[i][j], v);
3730 }
3731 Polyhedron *P = Constraints2Polyhedron(M, 1024);
3732 Matrix_Free(M);
3733 Polyhedron *P2;
3734 Polyhedron *U = Universe_Polyhedron(1);
3735 P2 = P;
3736 Param_Polyhedron *PP =
3737 Polyhedron2Param_SimplifiedDomain(&P, U, 1024, NULL, NULL);
3738 Polyhedron_Free(P);
3739 if (P2 != P)
3740 Polyhedron_Free(P2);
3741 Polyhedron_Free(U);
3742 Param_Vertices *V;
3743 for (i = 0, V = PP->V; V; ++i, V = V->next)
3744 ;
3745 if (PP)
3746 Param_Polyhedron_Free(PP);
3747 if (i != 10) {
3748 fprintf(stderr, "WARNING: results may be incorrect\n");
3749 fprintf(stderr,
3750 "WARNING: use latest version of PolyLib to remove this warning\n");
3751 }
3752 }
3753
3754 return Polyhedron2Param_SimplifiedDomain(Din, Cin, WS, CEq, CT);
3755 }
3756 #endif
3757
3758 static evalue* barvinok_enumerate_ev_f(Polyhedron *P, Polyhedron* C,
3759 unsigned MaxRays);
3760
3761 /* Destroys C */
3762 static evalue* barvinok_enumerate_cst(Polyhedron *P, Polyhedron* C,
3763 unsigned MaxRays)
3764 {
3765 evalue *eres;
3766
3767 ALLOC(evalue, eres);
3768 value_init(eres->d);
3769 value_set_si(eres->d, 0);
3770 eres->x.p = new_enode(partition, 2, C->Dimension);
3771 EVALUE_SET_DOMAIN(eres->x.p->arr[0], DomainConstraintSimplify(C, MaxRays));
3772 value_set_si(eres->x.p->arr[1].d, 1);
3773 value_init(eres->x.p->arr[1].x.n);
3774 if (emptyQ(P))
3775 value_set_si(eres->x.p->arr[1].x.n, 0);
3776 else
3777 barvinok_count(P, &eres->x.p->arr[1].x.n, MaxRays);
3778
3779 return eres;
3780 }
3781
3782 evalue* barvinok_enumerate_ev(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
3783 {
3784 //P = unfringe(P, MaxRays);
3785 Polyhedron *Corig = C;
3786 Polyhedron *CEq = NULL, *rVD, *CA;
3787 int r = 0;
3788 unsigned nparam = C->Dimension;
3789 evalue *eres;
3790
3791 evalue factor;
3792 value_init(factor.d);
3793 evalue_set_si(&factor, 1, 1);
3794
3795 /* for now */
3796 POL_ENSURE_FACETS(P);
3797 POL_ENSURE_VERTICES(P);
3798 POL_ENSURE_FACETS(C);
3799 POL_ENSURE_VERTICES(C);
3800
3801 CA = align_context(C, P->Dimension, MaxRays);
3802 P = DomainIntersection(P, CA, MaxRays);
3803 Polyhedron_Free(CA);
3804
3805 if (C->Dimension == 0 || emptyQ(P)) {
3806 constant:
3807 eres = barvinok_enumerate_cst(P, CEq ? CEq : Polyhedron_Copy(C),
3808 MaxRays);
3809 out:
3810 emul(&factor, eres);
3811 reduce_evalue(eres);
3812 free_evalue_refs(&factor);
3813 Domain_Free(P);
3814 if (C != Corig)
3815 Polyhedron_Free(C);
3816
3817 return eres;
3818 }
3819 if (Polyhedron_is_infinite(P, nparam))
3820 goto constant;
3821
3822 if (P->NbEq != 0) {
3823 Matrix *f;
3824 P = remove_equalities_p(P, P->Dimension-nparam, &f);
3825 mask(f, &factor);
3826 Matrix_Free(f);
3827 }
3828 if (P->Dimension == nparam) {
3829 CEq = P;
3830 P = Universe_Polyhedron(0);
3831 goto constant;
3832 }
3833
3834 Polyhedron *T = Polyhedron_Factor(P, nparam, MaxRays);
3835 if (T) {
3836 Polyhedron *Q;
3837 Polyhedron *C2;
3838 for (Q = T; Q; Q = Q->next) {
3839 Polyhedron *next = Q->next;
3840 Q->next = NULL;
3841
3842 Polyhedron *QC = Q;
3843 if (Q->Dimension != C->Dimension)
3844 QC = Polyhedron_Project(Q, nparam);
3845
3846 C2 = C;
3847 C = DomainIntersection(C, QC, MaxRays);
3848 if (C2 != Corig)
3849 Polyhedron_Free(C2);
3850 if (QC != Q)
3851 Polyhedron_Free(QC);
3852
3853 Q->next = next;
3854 }
3855 Polyhedron_Free(P);
3856 P = T;
3857 if (T->Dimension == C->Dimension) {
3858 P = T->next;
3859 T->next = NULL;
3860 Polyhedron_Free(T);
3861 }
3862 }
3863
3864 Polyhedron *next = P->next;
3865 P->next = NULL;
3866 eres = barvinok_enumerate_ev_f(P, C, MaxRays);
3867 P->next = next;
3868
3869 if (P->next) {
3870 Polyhedron *Q;
3871 evalue *f;
3872
3873 for (Q = P->next; Q; Q = Q->next) {
3874 Polyhedron *next = Q->next;
3875 Q->next = NULL;
3876
3877 f = barvinok_enumerate_ev_f(Q, C, MaxRays);
3878 emul(f, eres);
3879 free_evalue_refs(f);
3880 free(f);
3881
3882 Q->next = next;
3883 }
3884 }
3885
3886 goto out;
3887 }
3888
3889 static evalue* barvinok_enumerate_ev_f(Polyhedron *P, Polyhedron* C,
3890 unsigned MaxRays)
3891 {
3892 unsigned nparam = C->Dimension;
3893
3894 if (P->Dimension - nparam == 1)
3895 return ParamLine_Length(P, C, MaxRays);
3896
3897 Param_Polyhedron *PP = NULL;
3898 Polyhedron *CEq = NULL, *pVD;
3899 Matrix *CT = NULL;
3900 Param_Domain *D, *next;
3901 Param_Vertices *V;
3902 evalue *eres;
3903 Polyhedron *Porig = P;
3904
3905 PP = Polyhedron2Param_SD(&P,C,MaxRays,&CEq,&CT);
3906
3907 if (isIdentity(CT)) {
3908 Matrix_Free(CT);
3909 CT = NULL;
3910 } else {
3911 assert(CT->NbRows != CT->NbColumns);
3912 if (CT->NbRows == 1) { // no more parameters
3913 eres = barvinok_enumerate_cst(P, CEq, MaxRays);
3914 out:
3915 if (CT)
3916 Matrix_Free(CT);
3917 if (PP)
3918 Param_Polyhedron_Free(PP);
3919 if (P != Porig)
3920 Polyhedron_Free(P);
3921
3922 return eres;
3923 }
3924 nparam = CT->NbRows - 1;
3925 }
3926
3927 unsigned dim = P->Dimension - nparam;
3928
3929 ALLOC(evalue, eres);
3930 value_init(eres->d);
3931 value_set_si(eres->d, 0);
3932
3933 int nd;
3934 for (nd = 0, D=PP->D; D; ++nd, D=D->next);
3935 struct section { Polyhedron *D; evalue E; };
3936 section *s = new section[nd];
3937 Polyhedron **fVD = new Polyhedron_p[nd];
3938
3939 try_again:
3940 #ifdef USE_INCREMENTAL_BF
3941 bfenumerator et(P, dim, PP->nbV);
3942 #elif defined USE_INCREMENTAL_DF
3943 ienumerator et(P, dim, PP->nbV);
3944 #else
3945 enumerator et(P, dim, PP->nbV);
3946 #endif
3947
3948 for(nd = 0, D=PP->D; D; D=next) {
3949 next = D->next;
3950
3951 Polyhedron *rVD = reduce_domain(D->Domain, CT, CEq,
3952 fVD, nd, MaxRays);
3953 if (!rVD)
3954 continue;
3955
3956 pVD = CT ? DomainImage(rVD,CT,MaxRays) : rVD;
3957
3958 value_init(s[nd].E.d);
3959 evalue_set_si(&s[nd].E, 0, 1);
3960 s[nd].D = rVD;
3961
3962 FORALL_PVertex_in_ParamPolyhedron(V,D,PP) // _i is internal counter
3963 if (!et.vE[_i])
3964 try {
3965 et.decompose_at(V, _i, MaxRays);
3966 } catch (OrthogonalException &e) {
3967 if (rVD != pVD)
3968 Domain_Free(pVD);
3969 for (; nd >= 0; --nd) {
3970 free_evalue_refs(&s[nd].E);
3971 Domain_Free(s[nd].D);
3972 Domain_Free(fVD[nd]);
3973 }
3974 goto try_again;
3975 }
3976 eadd(et.vE[_i] , &s[nd].E);
3977 END_FORALL_PVertex_in_ParamPolyhedron;
3978 reduce_in_domain(&s[nd].E, pVD);
3979
3980 if (CT)
3981 addeliminatedparams_evalue(&s[nd].E, CT);
3982 ++nd;
3983 if (rVD != pVD)
3984 Domain_Free(pVD);
3985 }
3986
3987 if (nd == 0)
3988 evalue_set_si(eres, 0, 1);
3989 else {
3990 eres->x.p = new_enode(partition, 2*nd, C->Dimension);
3991 for (int j = 0; j < nd; ++j) {
3992 EVALUE_SET_DOMAIN(eres->x.p->arr[2*j], s[j].D);
3993 value_clear(eres->x.p->arr[2*j+1].d);
3994 eres->x.p->arr[2*j+1] = s[j].E;
3995 Domain_Free(fVD[j]);
3996 }
3997 }
3998 delete [] s;
3999 delete [] fVD;
4000
4001 if (CEq)
4002 Polyhedron_Free(CEq);
4003 goto out;
4004 }
4005
4006 Enumeration* barvinok_enumerate(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
4007 {
4008 evalue *EP = barvinok_enumerate_ev(P, C, MaxRays);
4009
4010 return partition2enumeration(EP);
4011 }
4012
4013 static void SwapColumns(Value **V, int n, int i, int j)
4014 {
4015 for (int r = 0; r < n; ++r)
4016 value_swap(V[r][i], V[r][j]);
4017 }
4018
4019 static void SwapColumns(Polyhedron *P, int i, int j)
4020 {
4021 SwapColumns(P->Constraint, P->NbConstraints, i, j);
4022 SwapColumns(P->Ray, P->NbRays, i, j);
4023 }
4024
4025 static void negative_test_constraint(Value *l, Value *u, Value *c, int pos,
4026 int len, Value *v)
4027 {
4028 value_oppose(*v, u[pos+1]);
4029 Vector_Combine(l+1, u+1, c+1, *v, l[pos+1], len-1);
4030 value_multiply(*v, *v, l[pos+1]);
4031 value_subtract(c[len-1], c[len-1], *v);
4032 value_set_si(*v, -1);
4033 Vector_Scale(c+1, c+1, *v, len-1);
4034 value_decrement(c[len-1], c[len-1]);
4035 ConstraintSimplify(c, c, len, v);
4036 }
4037
4038 static bool parallel_constraints(Value *l, Value *u, Value *c, int pos,
4039 int len)
4040 {
4041 bool parallel;
4042 Value g1;
4043 Value g2;
4044 value_init(g1);
4045 value_init(g2);
4046
4047 Vector_Gcd(&l[1+pos], len, &g1);
4048 Vector_Gcd(&u[1+pos], len, &g2);
4049 Vector_Combine(l+1+pos, u+1+pos, c+1, g2, g1, len);
4050 parallel = First_Non_Zero(c+1, len) == -1;
4051
4052 value_clear(g1);
4053 value_clear(g2);
4054
4055 return parallel;
4056 }
4057
4058 static void negative_test_constraint7(Value *l, Value *u, Value *c, int pos,
4059 int exist, int len, Value *v)
4060 {
4061 Value g;
4062 value_init(g);
4063
4064 Vector_Gcd(&u[1+pos], exist, v);
4065 Vector_Gcd(&l[1+pos], exist, &g);
4066 Vector_Combine(l+1, u+1, c+1, *v, g, len-1);
4067 value_multiply(*v, *v, g);
4068 value_subtract(c[len-1], c[len-1], *v);
4069 value_set_si(*v, -1);
4070 Vector_Scale(c+1, c+1, *v, len-1);
4071 value_decrement(c[len-1], c[len-1]);
4072 ConstraintSimplify(c, c, len, v);
4073
4074 value_clear(g);
4075 }
4076
4077 static void oppose_constraint(Value *c, int len, Value *v)
4078 {
4079 value_set_si(*v, -1);
4080 Vector_Scale(c+1, c+1, *v, len-1);
4081 value_decrement(c[len-1], c[len-1]);
4082 }
4083
4084 static bool SplitOnConstraint(Polyhedron *P, int i, int l, int u,
4085 int nvar, int len, int exist, int MaxRays,
4086 Vector *row, Value& f, bool independent,
4087 Polyhedron **pos, Polyhedron **neg)
4088 {
4089 negative_test_constraint(P->Constraint[l], P->Constraint[u],
4090 row->p, nvar+i, len, &f);
4091 *neg = AddConstraints(row->p, 1, P, MaxRays);
4092
4093 /* We found an independent, but useless constraint
4094 * Maybe we should detect this earlier and not
4095 * mark the variable as INDEPENDENT
4096 */
4097 if (emptyQ((*neg))) {
4098 Polyhedron_Free(*neg);
4099 return false;
4100 }
4101
4102 oppose_constraint(row->p, len, &f);
4103 *pos = AddConstraints(row->p, 1, P, MaxRays);
4104
4105 if (emptyQ((*pos))) {
4106 Polyhedron_Free(*neg);
4107 Polyhedron_Free(*pos);
4108 return false;
4109 }
4110
4111 return true;
4112 }
4113
4114 /*
4115 * unimodularly transform P such that constraint r is transformed
4116 * into a constraint that involves only a single (the first)
4117 * existential variable
4118 *
4119 */
4120 static Polyhedron *rotate_along(Polyhedron *P, int r, int nvar, int exist,
4121 unsigned MaxRays)
4122 {
4123 Value g;
4124 value_init(g);
4125
4126 Vector *row = Vector_Alloc(exist);
4127 Vector_Copy(P->Constraint[r]+1+nvar, row->p, exist);
4128 Vector_Gcd(row->p, exist, &g);
4129 if (value_notone_p(g))
4130 Vector_AntiScale(row->p, row->p, g, exist);
4131 value_clear(g);
4132
4133 Matrix *M = unimodular_complete(row);
4134 Matrix *M2 = Matrix_Alloc(P->Dimension+1, P->Dimension+1);
4135 for (r = 0; r < nvar; ++r)
4136 value_set_si(M2->p[r][r], 1);
4137 for ( ; r < nvar+exist; ++r)
4138 Vector_Copy(M->p[r-nvar], M2->p[r]+nvar, exist);
4139 for ( ; r < P->Dimension+1; ++r)
4140 value_set_si(M2->p[r][r], 1);
4141 Polyhedron *T = Polyhedron_Image(P, M2, MaxRays);
4142
4143 Matrix_Free(M2);
4144 Matrix_Free(M);
4145 Vector_Free(row);
4146
4147 return T;
4148 }
4149
4150 static bool SplitOnVar(Polyhedron *P, int i,
4151 int nvar, int len, int exist, int MaxRays,
4152 Vector *row, Value& f, bool independent,
4153 Polyhedron **pos, Polyhedron **neg)
4154 {
4155 int j;
4156
4157 for (int l = P->NbEq; l < P->NbConstraints; ++l) {
4158 if (value_negz_p(P->Constraint[l][nvar+i+1]))
4159 continue;
4160
4161 if (independent) {
4162 for (j = 0; j < exist; ++j)
4163 if (j != i && value_notzero_p(P->Constraint[l][nvar+j+1]))
4164 break;
4165 if (j < exist)
4166 continue;
4167 }
4168
4169 for (int u = P->NbEq; u < P->NbConstraints; ++u) {
4170 if (value_posz_p(P->Constraint[u][nvar+i+1]))
4171 continue;
4172
4173 if (independent) {
4174 for (j = 0; j < exist; ++j)
4175 if (j != i && value_notzero_p(P->Constraint[u][nvar+j+1]))
4176 break;
4177 if (j < exist)
4178 continue;
4179 }
4180
4181 if (SplitOnConstraint(P, i, l, u,
4182 nvar, len, exist, MaxRays,
4183 row, f, independent,
4184 pos, neg)) {
4185 if (independent) {
4186 if (i != 0)
4187 SwapColumns(*neg, nvar+1, nvar+1+i);
4188 }
4189 return true;
4190 }
4191 }
4192 }
4193
4194 return false;
4195 }
4196
4197 static bool double_bound_pair(Polyhedron *P, int nvar, int exist,
4198 int i, int l1, int l2,
4199 Polyhedron **pos, Polyhedron **neg)
4200 {
4201 Value f;
4202 value_init(f);
4203 Vector *row = Vector_Alloc(P->Dimension+2);
4204 value_set_si(row->p[0], 1);
4205 value_oppose(f, P->Constraint[l1][nvar+i+1]);
4206 Vector_Combine(P->Constraint[l1]+1, P->Constraint[l2]+1,
4207 row->p+1,
4208 P->Constraint[l2][nvar+i+1], f,
4209 P->Dimension+1);
4210 ConstraintSimplify(row->p, row->p, P->Dimension+2, &f);
4211 *pos = AddConstraints(row->p, 1, P, 0);
4212 value_set_si(f, -1);
4213 Vector_Scale(row->p+1, row->p+1, f, P->Dimension+1);
4214 value_decrement(row->p[P->Dimension+1], row->p[P->Dimension+1]);
4215 *neg = AddConstraints(row->p, 1, P, 0);
4216 Vector_Free(row);
4217 value_clear(f);
4218
4219 return !emptyQ((*pos)) && !emptyQ((*neg));
4220 }
4221
4222 static bool double_bound(Polyhedron *P, int nvar, int exist,
4223 Polyhedron **pos, Polyhedron **neg)
4224 {
4225 for (int i = 0; i < exist; ++i) {
4226 int l1, l2;
4227 for (l1 = P->NbEq; l1 < P->NbConstraints; ++l1) {
4228 if (value_negz_p(P->Constraint[l1][nvar+i+1]))
4229 continue;
4230 for (l2 = l1 + 1; l2 < P->NbConstraints; ++l2) {
4231 if (value_negz_p(P->Constraint[l2][nvar+i+1]))
4232 continue;
4233 if (double_bound_pair(P, nvar, exist, i, l1, l2, pos, neg))
4234 return true;
4235 }
4236 }
4237 for (l1 = P->NbEq; l1 < P->NbConstraints; ++l1) {
4238 if (value_posz_p(P->Constraint[l1][nvar+i+1]))
4239 continue;
4240 if (l1 < P->NbConstraints)
4241 for (l2 = l1 + 1; l2 < P->NbConstraints; ++l2) {
4242 if (value_posz_p(P->Constraint[l2][nvar+i+1]))
4243 continue;
4244 if (double_bound_pair(P, nvar, exist, i, l1, l2, pos, neg))
4245 return true;
4246 }
4247 }
4248 return false;
4249 }
4250 return false;
4251 }
4252
4253 enum constraint {
4254 ALL_POS = 1 << 0,
4255 ONE_NEG = 1 << 1,
4256 INDEPENDENT = 1 << 2,
4257 ROT_NEG = 1 << 3
4258 };
4259
4260 static evalue* enumerate_or(Polyhedron *D,
4261 unsigned exist, unsigned nparam, unsigned MaxRays)
4262 {
4263 #ifdef DEBUG_ER
4264 fprintf(stderr, "\nER: Or\n");
4265 #endif /* DEBUG_ER */
4266
4267 Polyhedron *N = D->next;
4268 D->next = 0;
4269 evalue *EP =
4270 barvinok_enumerate_e(D, exist, nparam, MaxRays);
4271 Polyhedron_Free(D);
4272
4273 for (D = N; D; D = N) {
4274 N = D->next;
4275 D->next = 0;
4276
4277 evalue *EN =
4278 barvinok_enumerate_e(D, exist, nparam, MaxRays);
4279
4280 eor(EN, EP);
4281 free_evalue_refs(EN);
4282 free(EN);
4283 Polyhedron_Free(D);
4284 }
4285
4286 reduce_evalue(EP);
4287
4288 return EP;
4289 }
4290
4291 static evalue* enumerate_sum(Polyhedron *P,
4292 unsigned exist, unsigned nparam, unsigned MaxRays)
4293 {
4294 int nvar = P->Dimension - exist - nparam;
4295 int toswap = nvar < exist ? nvar : exist;
4296 for (int i = 0; i < toswap; ++i)
4297 SwapColumns(P, 1 + i, nvar+exist - i);
4298 nparam += nvar;
4299
4300 #ifdef DEBUG_ER
4301 fprintf(stderr, "\nER: Sum\n");
4302 #endif /* DEBUG_ER */
4303
4304 evalue *EP = barvinok_enumerate_e(P, exist, nparam, MaxRays);
4305
4306 for (int i = 0; i < /* nvar */ nparam; ++i) {
4307 Matrix *C = Matrix_Alloc(1, 1 + nparam + 1);
4308 value_set_si(C->p[0][0], 1);
4309 evalue split;
4310 value_init(split.d);
4311 value_set_si(split.d, 0);
4312 split.x.p = new_enode(partition, 4, nparam);
4313 value_set_si(C->p[0][1+i], 1);
4314 Matrix *C2 = Matrix_Copy(C);
4315 EVALUE_SET_DOMAIN(split.x.p->arr[0],
4316 Constraints2Polyhedron(C2, MaxRays));
4317 Matrix_Free(C2);
4318 evalue_set_si(&split.x.p->arr[1], 1, 1);
4319 value_set_si(C->p[0][1+i], -1);
4320 value_set_si(C->p[0][1+nparam], -1);
4321 EVALUE_SET_DOMAIN(split.x.p->arr[2],
4322 Constraints2Polyhedron(C, MaxRays));
4323 evalue_set_si(&split.x.p->arr[3], 1, 1);
4324 emul(&split, EP);
4325 free_evalue_refs(&split);
4326 Matrix_Free(C);
4327 }
4328 reduce_evalue(EP);
4329 evalue_range_reduction(EP);
4330
4331 evalue_frac2floor(EP);
4332
4333 evalue *sum = esum(EP, nvar);
4334
4335 free_evalue_refs(EP);
4336 free(EP);
4337 EP = sum;
4338
4339 evalue_range_reduction(EP);
4340
4341 return EP;
4342 }
4343
4344 static evalue* split_sure(Polyhedron *P, Polyhedron *S,
4345 unsigned exist, unsigned nparam, unsigned MaxRays)
4346 {
4347 int nvar = P->Dimension - exist - nparam;
4348
4349 Matrix *M = Matrix_Alloc(exist, S->Dimension+2);
4350 for (int i = 0; i < exist; ++i)
4351 value_set_si(M->p[i][nvar+i+1], 1);
4352 Polyhedron *O = S;
4353 S = DomainAddRays(S, M, MaxRays);
4354 Polyhedron_Free(O);
4355 Polyhedron *F = DomainAddRays(P, M, MaxRays);
4356 Polyhedron *D = DomainDifference(F, S, MaxRays);
4357 O = D;
4358 D = Disjoint_Domain(D, 0, MaxRays);
4359 Polyhedron_Free(F);
4360 Domain_Free(O);
4361 Matrix_Free(M);
4362
4363 M = Matrix_Alloc(P->Dimension+1-exist, P->Dimension+1);
4364 for (int j = 0; j < nvar; ++j)
4365 value_set_si(M->p[j][j], 1);
4366 for (int j = 0; j < nparam+1; ++j)
4367 value_set_si(M->p[nvar+j][nvar+exist+j], 1);
4368 Polyhedron *T = Polyhedron_Image(S, M, MaxRays);
4369 evalue *EP = barvinok_enumerate_e(T, 0, nparam, MaxRays);
4370 Polyhedron_Free(S);
4371 Polyhedron_Free(T);
4372 Matrix_Free(M);
4373
4374 for (Polyhedron *Q = D; Q; Q = Q->next) {
4375 Polyhedron *N = Q->next;
4376 Q->next = 0;
4377 T = DomainIntersection(P, Q, MaxRays);
4378 evalue *E = barvinok_enumerate_e(T, exist, nparam, MaxRays);
4379 eadd(E, EP);
4380 free_evalue_refs(E);
4381 free(E);
4382 Polyhedron_Free(T);
4383 Q->next = N;
4384 }
4385 Domain_Free(D);
4386 return EP;
4387 }
4388
4389 static evalue* enumerate_sure(Polyhedron *P,
4390 unsigned exist, unsigned nparam, unsigned MaxRays)
4391 {
4392 int i;
4393 Polyhedron *S = P;
4394 int nvar = P->Dimension - exist - nparam;
4395 Value lcm;
4396 Value f;
4397 value_init(lcm);
4398 value_init(f);
4399
4400 for (i = 0; i < exist; ++i) {
4401 Matrix *M = Matrix_Alloc(S->NbConstraints, S->Dimension+2);
4402 int c = 0;
4403 value_set_si(lcm, 1);
4404 for (int j = 0; j < S->NbConstraints; ++j) {
4405 if (value_negz_p(S->Constraint[j][1+nvar+i]))
4406 continue;
4407 if (value_one_p(S->Constraint[j][1+nvar+i]))
4408 continue;
4409 value_lcm(lcm, S->Constraint[j][1+nvar+i], &lcm);
4410 }
4411
4412 for (int j = 0; j < S->NbConstraints; ++j) {
4413 if (value_negz_p(S->Constraint[j][1+nvar+i]))
4414 continue;
4415 if (value_one_p(S->Constraint[j][1+nvar+i]))
4416 continue;
4417 value_division(f, lcm, S->Constraint[j][1+nvar+i]);
4418 Vector_Scale(S->Constraint[j], M->p[c], f, S->Dimension+2);
4419 value_subtract(M->p[c][S->Dimension+1],
4420 M->p[c][S->Dimension+1],
4421 lcm);
4422 value_increment(M->p[c][S->Dimension+1],
4423 M->p[c][S->Dimension+1]);
4424 ++c;
4425 }
4426 Polyhedron *O = S;
4427 S = AddConstraints(M->p[0], c, S, MaxRays);
4428 if (O != P)
4429 Polyhedron_Free(O);
4430 Matrix_Free(M);
4431 if (emptyQ(S)) {
4432 Polyhedron_Free(S);
4433 value_clear(lcm);
4434 value_clear(f);
4435 return 0;
4436 }
4437 }
4438 value_clear(lcm);
4439 value_clear(f);
4440
4441 #ifdef DEBUG_ER
4442 fprintf(stderr, "\nER: Sure\n");
4443 #endif /* DEBUG_ER */
4444
4445 return split_sure(P, S, exist, nparam, MaxRays);
4446 }
4447
4448 static evalue* enumerate_sure2(Polyhedron *P,
4449 unsigned exist, unsigned nparam, unsigned MaxRays)
4450 {
4451 int nvar = P->Dimension - exist - nparam;
4452 int r;
4453 for (r = 0; r < P->NbRays; ++r)
4454 if (value_one_p(P->Ray[r][0]) &&
4455 value_one_p(P->Ray[r][P->Dimension+1]))
4456 break;
4457
4458 if (r >= P->NbRays)
4459 return 0;
4460
4461 Matrix *M = Matrix_Alloc(nvar + 1 + nparam, P->Dimension+2);
4462 for (int i = 0; i < nvar; ++i)
4463 value_set_si(M->p[i][1+i], 1);
4464 for (int i = 0; i < nparam; ++i)
4465 value_set_si(M->p[i+nvar][1+nvar+exist+i], 1);
4466 Vector_Copy(P->Ray[r]+1+nvar, M->p[nvar+nparam]+1+nvar, exist);
4467 value_set_si(M->p[nvar+nparam][0], 1);
4468 value_set_si(M->p[nvar+nparam][P->Dimension+1], 1);
4469 Polyhedron * F = Rays2Polyhedron(M, MaxRays);
4470 Matrix_Free(M);
4471
4472 Polyhedron *I = DomainIntersection(F, P, MaxRays);
4473 Polyhedron_Free(F);
4474
4475 #ifdef DEBUG_ER
4476 fprintf(stderr, "\nER: Sure2\n");
4477 #endif /* DEBUG_ER */
4478
4479 return split_sure(P, I, exist, nparam, MaxRays);
4480 }
4481
4482 static evalue* enumerate_cyclic(Polyhedron *P,
4483 unsigned exist, unsigned nparam,
4484 evalue * EP, int r, int p, unsigned MaxRays)
4485 {
4486 int nvar = P->Dimension - exist - nparam;
4487
4488 /* If EP in its fractional maps only contains references
4489 * to the remainder parameter with appropriate coefficients
4490 * then we could in principle avoid adding existentially
4491 * quantified variables to the validity domains.
4492 * We'd have to replace the remainder by m { p/m }
4493 * and multiply with an appropriate factor that is one
4494 * only in the appropriate range.
4495 * This last multiplication can be avoided if EP
4496 * has a single validity domain with no (further)
4497 * constraints on the remainder parameter
4498 */
4499
4500 Matrix *CT = Matrix_Alloc(nparam+1, nparam+3);
4501 Matrix *M = Matrix_Alloc(1, 1+nparam+3);
4502 for (int j = 0; j < nparam; ++j)
4503 if (j != p)
4504 value_set_si(CT->p[j][j], 1);
4505 value_set_si(CT->p[p][nparam+1], 1);
4506 value_set_si(CT->p[nparam][nparam+2], 1);
4507 value_set_si(M->p[0][1+p], -1);
4508 value_absolute(M->p[0][1+nparam], P->Ray[0][1+nvar+exist+p]);
4509 value_set_si(M->p[0][1+nparam+1], 1);
4510 Polyhedron *CEq = Constraints2Polyhedron(M, 1);
4511 Matrix_Free(M);
4512 addeliminatedparams_enum(EP, CT, CEq, MaxRays, nparam);
4513 Polyhedron_Free(CEq);
4514 Matrix_Free(CT);
4515
4516 return EP;
4517 }
4518
4519 static void enumerate_vd_add_ray(evalue *EP, Matrix *Rays, unsigned MaxRays)
4520 {
4521 if (value_notzero_p(EP->d))
4522 return;
4523
4524 assert(EP->x.p->type == partition);
4525 assert(EP->x.p->pos == EVALUE_DOMAIN(EP->x.p->arr[0])->Dimension);
4526 for (int i = 0; i < EP->x.p->size/2; ++i) {
4527 Polyhedron *D = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
4528 Polyhedron *N = DomainAddRays(D, Rays, MaxRays);
4529 EVALUE_SET_DOMAIN(EP->x.p->arr[2*i], N);
4530 Domain_Free(D);
4531 }
4532 }
4533
4534 static evalue* enumerate_line(Polyhedron *P,
4535 unsigned exist, unsigned nparam, unsigned MaxRays)
4536 {
4537 if (P->NbBid == 0)
4538 return 0;
4539
4540 #ifdef DEBUG_ER
4541 fprintf(stderr, "\nER: Line\n");
4542 #endif /* DEBUG_ER */
4543
4544 int nvar = P->Dimension - exist - nparam;
4545 int i, j;
4546 for (i = 0; i < nparam; ++i)
4547 if (value_notzero_p(P->Ray[0][1+nvar+exist+i]))
4548 break;
4549 assert(i < nparam);
4550 for (j = i+1; j < nparam; ++j)
4551 if (value_notzero_p(P->Ray[0][1+nvar+exist+i]))
4552 break;
4553 assert(j >= nparam); // for now
4554
4555 Matrix *M = Matrix_Alloc(2, P->Dimension+2);
4556 value_set_si(M->p[0][0], 1);
4557 value_set_si(M->p[0][1+nvar+exist+i], 1);
4558 value_set_si(M->p[1][0], 1);
4559 value_set_si(M->p[1][1+nvar+exist+i], -1);
4560 value_absolute(M->p[1][1+P->Dimension], P->Ray[0][1+nvar+exist+i]);
4561 value_decrement(M->p[1][1+P->Dimension], M->p[1][1+P->Dimension]);
4562 Polyhedron *S = AddConstraints(M->p[0], 2, P, MaxRays);
4563 evalue *EP = barvinok_enumerate_e(S, exist, nparam, MaxRays);
4564 Polyhedron_Free(S);
4565 Matrix_Free(M);
4566
4567 return enumerate_cyclic(P, exist, nparam, EP, 0, i, MaxRays);
4568 }
4569
4570 static int single_param_pos(Polyhedron*P, unsigned exist, unsigned nparam,
4571 int r)
4572 {
4573 int nvar = P->Dimension - exist - nparam;
4574 if (First_Non_Zero(P->Ray[r]+1, nvar) != -1)
4575 return -1;
4576 int i = First_Non_Zero(P->Ray[r]+1+nvar+exist, nparam);
4577 if (i == -1)
4578 return -1;
4579 if (First_Non_Zero(P->Ray[r]+1+nvar+exist+1, nparam-i-1) != -1)
4580 return -1;
4581 return i;
4582 }
4583
4584 static evalue* enumerate_remove_ray(Polyhedron *P, int r,
4585 unsigned exist, unsigned nparam, unsigned MaxRays)
4586 {
4587 #ifdef DEBUG_ER
4588 fprintf(stderr, "\nER: RedundantRay\n");
4589 #endif /* DEBUG_ER */
4590
4591 Value one;
4592 value_init(one);
4593 value_set_si(one, 1);
4594 int len = P->NbRays-1;
4595 Matrix *M = Matrix_Alloc(2 * len, P->Dimension+2);
4596 Vector_Copy(P->Ray[0], M->p[0], r * (P->Dimension+2));
4597 Vector_Copy(P->Ray[r+1], M->p[r], (len-r) * (P->Dimension+2));
4598 for (int j = 0; j < P->NbRays; ++j) {
4599 if (j == r)
4600 continue;
4601 Vector_Combine(P->Ray[j], P->Ray[r], M->p[len+j-(j>r)],
4602 one, P->Ray[j][P->Dimension+1], P->Dimension+2);
4603 }
4604
4605 P = Rays2Polyhedron(M, MaxRays);
4606 Matrix_Free(M);
4607 evalue *EP = barvinok_enumerate_e(P, exist, nparam, MaxRays);
4608 Polyhedron_Free(P);
4609 value_clear(one);
4610
4611 return EP;
4612 }
4613
4614 static evalue* enumerate_redundant_ray(Polyhedron *P,
4615 unsigned exist, unsigned nparam, unsigned MaxRays)
4616 {
4617 assert(P->NbBid == 0);
4618 int nvar = P->Dimension - exist - nparam;
4619 Value m;
4620 value_init(m);
4621
4622 for (int r = 0; r < P->NbRays; ++r) {
4623 if (value_notzero_p(P->Ray[r][P->Dimension+1]))
4624 continue;
4625 int i1 = single_param_pos(P, exist, nparam, r);
4626 if (i1 == -1)
4627 continue;
4628 for (int r2 = r+1; r2 < P->NbRays; ++r2) {
4629 if (value_notzero_p(P->Ray[r2][P->Dimension+1]))
4630 continue;
4631 int i2 = single_param_pos(P, exist, nparam, r2);
4632 if (i2 == -1)
4633 continue;
4634 if (i1 != i2)
4635 continue;
4636
4637 value_division(m, P->Ray[r][1+nvar+exist+i1],
4638 P->Ray[r2][1+nvar+exist+i1]);
4639 value_multiply(m, m, P->Ray[r2][1+nvar+exist+i1]);
4640 /* r2 divides r => r redundant */
4641 if (value_eq(m, P->Ray[r][1+nvar+exist+i1])) {
4642 value_clear(m);
4643 return enumerate_remove_ray(P, r, exist, nparam, MaxRays);
4644 }
4645
4646 value_division(m, P->Ray[r2][1+nvar+exist+i1],
4647 P->Ray[r][1+nvar+exist+i1]);
4648 value_multiply(m, m, P->Ray[r][1+nvar+exist+i1]);
4649 /* r divides r2 => r2 redundant */
4650 if (value_eq(m, P->Ray[r2][1+nvar+exist+i1])) {
4651 value_clear(m);
4652 return enumerate_remove_ray(P, r2, exist, nparam, MaxRays);
4653 }
4654 }
4655 }
4656 value_clear(m);
4657 return 0;
4658 }
4659
4660 static Polyhedron *upper_bound(Polyhedron *P,
4661 int pos, Value *max, Polyhedron **R)
4662 {
4663 Value v;
4664 int r;
4665 value_init(v);
4666
4667 *R = 0;
4668 Polyhedron *N;
4669 Polyhedron *B = 0;
4670 for (Polyhedron *Q = P; Q; Q = N) {
4671 N = Q->next;
4672 for (r = 0; r < P->NbRays; ++r) {
4673 if (value_zero_p(P->Ray[r][P->Dimension+1]) &&
4674 value_pos_p(P->Ray[r][1+pos]))
4675 break;
4676 }
4677 if (r < P->NbRays) {
4678 Q->next = *R;
4679 *R = Q;
4680 continue;
4681 } else {
4682 Q->next = B;
4683 B = Q;
4684 }
4685 for (r = 0; r < P->NbRays; ++r) {
4686 if (value_zero_p(P->Ray[r][P->Dimension+1]))
4687 continue;
4688 mpz_fdiv_q(v, P->Ray[r][1+pos], P->Ray[r][1+P->Dimension]);
4689 if ((!Q->next && r == 0) || value_gt(v, *max))
4690 value_assign(*max, v);
4691 }
4692 }
4693 value_clear(v);
4694 return B;
4695 }
4696
4697 static evalue* enumerate_ray(Polyhedron *P,
4698 unsigned exist, unsigned nparam, unsigned MaxRays)
4699 {
4700 assert(P->NbBid == 0);
4701 int nvar = P->Dimension - exist - nparam;
4702
4703 int r;
4704 for (r = 0; r < P->NbRays; ++r)
4705 if (value_zero_p(P->Ray[r][P->Dimension+1]))
4706 break;
4707 if (r >= P->NbRays)
4708 return 0;
4709
4710 int r2;
4711 for (r2 = r+1; r2 < P->NbRays; ++r2)
4712 if (value_zero_p(P->Ray[r2][P->Dimension+1]))
4713 break;
4714 if (r2 < P->NbRays) {
4715 if (nvar > 0)
4716 return enumerate_sum(P, exist, nparam, MaxRays);
4717 }
4718
4719 #ifdef DEBUG_ER
4720 fprintf(stderr, "\nER: Ray\n");
4721 #endif /* DEBUG_ER */
4722
4723 Value m;
4724 Value one;
4725 value_init(m);
4726 value_init(one);
4727 value_set_si(one, 1);
4728 int i = single_param_pos(P, exist, nparam, r);
4729 assert(i != -1); // for now;
4730
4731 Matrix *M = Matrix_Alloc(P->NbRays, P->Dimension+2);
4732 for (int j = 0; j < P->NbRays; ++j) {
4733 Vector_Combine(P->Ray[j], P->Ray[r], M->p[j],
4734 one, P->Ray[j][P->Dimension+1], P->Dimension+2);
4735 }
4736 Polyhedron *S = Rays2Polyhedron(M, MaxRays);
4737 Matrix_Free(M);
4738 Polyhedron *D = DomainDifference(P, S, MaxRays);
4739 Polyhedron_Free(S);
4740 // Polyhedron_Print(stderr, P_VALUE_FMT, D);
4741 assert(value_pos_p(P->Ray[r][1+nvar+exist+i])); // for now
4742 Polyhedron *R;
4743 D = upper_bound(D, nvar+exist+i, &m, &R);
4744 assert(D);
4745 Domain_Free(D);
4746
4747 M = Matrix_Alloc(2, P->Dimension+2);
4748 value_set_si(M->p[0][0], 1);
4749 value_set_si(M->p[1][0], 1);
4750 value_set_si(M->p[0][1+nvar+exist+i], -1);
4751 value_set_si(M->p[1][1+nvar+exist+i], 1);
4752 value_assign(M->p[0][1+P->Dimension], m);
4753 value_oppose(M->p[1][1+P->Dimension], m);
4754 value_addto(M->p[1][1+P->Dimension], M->p[1][1+P->Dimension],
4755 P->Ray[r][1+nvar+exist+i]);
4756 value_decrement(M->p[1][1+P->Dimension], M->p[1][1+P->Dimension]);
4757 // Matrix_Print(stderr, P_VALUE_FMT, M);
4758 D = AddConstraints(M->p[0], 2, P, MaxRays);
4759 // Polyhedron_Print(stderr, P_VALUE_FMT, D);
4760 value_subtract(M->p[0][1+P->Dimension], M->p[0][1+P->Dimension],
4761 P->Ray[r][1+nvar+exist+i]);
4762 // Matrix_Print(stderr, P_VALUE_FMT, M);
4763 S = AddConstraints(M->p[0], 1, P, MaxRays);
4764 // Polyhedron_Print(stderr, P_VALUE_FMT, S);
4765 Matrix_Free(M);
4766
4767 evalue *EP = barvinok_enumerate_e(D, exist, nparam, MaxRays);
4768 Polyhedron_Free(D);
4769 value_clear(one);
4770 value_clear(m);
4771
4772 if (value_notone_p(P->Ray[r][1+nvar+exist+i]))
4773 EP = enumerate_cyclic(P, exist, nparam, EP, r, i, MaxRays);
4774 else {
4775 M = Matrix_Alloc(1, nparam+2);
4776 value_set_si(M->p[0][0], 1);
4777 value_set_si(M->p[0][1+i], 1);
4778 enumerate_vd_add_ray(EP, M, MaxRays);
4779 Matrix_Free(M);
4780 }
4781
4782 if (!emptyQ(S)) {
4783 evalue *E = barvinok_enumerate_e(S, exist, nparam, MaxRays);
4784 eadd(E, EP);
4785 free_evalue_refs(E);
4786 free(E);
4787 }
4788 Polyhedron_Free(S);
4789
4790 if (R) {
4791 assert(nvar == 0);
4792 evalue *ER = enumerate_or(R, exist, nparam, MaxRays);
4793 eor(ER, EP);
4794 free_evalue_refs(ER);
4795 free(ER);
4796 }
4797
4798 return EP;
4799 }
4800
4801 static evalue* enumerate_vd(Polyhedron **PA,
4802 unsigned exist, unsigned nparam, unsigned MaxRays)
4803 {
4804 Polyhedron *P = *PA;
4805 int nvar = P->Dimension - exist - nparam;
4806 Param_Polyhedron *PP = NULL;
4807 Polyhedron *C = Universe_Polyhedron(nparam);
4808 Polyhedron *CEq;
4809 Matrix *CT;
4810 Polyhedron *PR = P;
4811 PP = Polyhedron2Param_SimplifiedDomain(&PR,C,MaxRays,&CEq,&CT);
4812 Polyhedron_Free(C);
4813
4814 int nd;
4815 Param_Domain *D, *last;
4816 Value c;
4817 value_init(c);
4818 for (nd = 0, D=PP->D; D; D=D->next, ++nd)
4819 ;
4820
4821 Polyhedron **VD = new Polyhedron_p[nd];
4822 Polyhedron **fVD = new Polyhedron_p[nd];
4823 for(nd = 0, D=PP->D; D; D=D->next) {
4824 Polyhedron *rVD = reduce_domain(D->Domain, CT, CEq,
4825 fVD, nd, MaxRays);
4826 if (!rVD)
4827 continue;
4828
4829 VD[nd++] = rVD;
4830 last = D;
4831 }
4832
4833 evalue *EP = 0;
4834
4835 if (nd == 0)
4836 EP = new_zero_ep();
4837
4838 /* This doesn't seem to have any effect */
4839 if (nd == 1) {
4840 Polyhedron *CA = align_context(VD[0], P->Dimension, MaxRays);
4841 Polyhedron *O = P;
4842 P = DomainIntersection(P, CA, MaxRays);
4843 if (O != *PA)
4844 Polyhedron_Free(O);
4845 Polyhedron_Free(CA);
4846 if (emptyQ(P))
4847 EP = new_zero_ep();
4848 }
4849
4850 if (!EP && CT->NbColumns != CT->NbRows) {
4851 Polyhedron *CEqr = DomainImage(CEq, CT, MaxRays);
4852 Polyhedron *CA = align_context(CEqr, PR->Dimension, MaxRays);
4853 Polyhedron *I = DomainIntersection(PR, CA, MaxRays);
4854 Polyhedron_Free(CEqr);
4855 Polyhedron_Free(CA);
4856 #ifdef DEBUG_ER
4857 fprintf(stderr, "\nER: Eliminate\n");
4858 #endif /* DEBUG_ER */
4859 nparam -= CT->NbColumns - CT->NbRows;
4860 EP = barvinok_enumerate_e(I, exist, nparam, MaxRays);
4861 nparam += CT->NbColumns - CT->NbRows;
4862 addeliminatedparams_enum(EP, CT, CEq, MaxRays, nparam);
4863 Polyhedron_Free(I);
4864 }
4865 if (PR != *PA)
4866 Polyhedron_Free(PR);
4867 PR = 0;
4868
4869 if (!EP && nd > 1) {
4870 #ifdef DEBUG_ER
4871 fprintf(stderr, "\nER: VD\n");
4872 #endif /* DEBUG_ER */
4873 for (int i = 0; i < nd; ++i) {
4874 Polyhedron *CA = align_context(VD[i], P->Dimension, MaxRays);
4875 Polyhedron *I = DomainIntersection(P, CA, MaxRays);
4876
4877 if (i == 0)
4878 EP = barvinok_enumerate_e(I, exist, nparam, MaxRays);
4879 else {
4880 evalue *E = barvinok_enumerate_e(I, exist, nparam, MaxRays);
4881 eadd(E, EP);
4882 free_evalue_refs(E);
4883 free(E);
4884 }
4885 Polyhedron_Free(I);
4886 Polyhedron_Free(CA);
4887 }
4888 }
4889
4890 for (int i = 0; i < nd; ++i) {
4891 Polyhedron_Free(VD[i]);
4892 Polyhedron_Free(fVD[i]);
4893 }
4894 delete [] VD;
4895 delete [] fVD;
4896 value_clear(c);
4897
4898 if (!EP && nvar == 0) {
4899 Value f;
4900 value_init(f);
4901 Param_Vertices *V, *V2;
4902 Matrix* M = Matrix_Alloc(1, P->Dimension+2);
4903
4904 FORALL_PVertex_in_ParamPolyhedron(V, last, PP) {
4905 bool found = false;
4906 FORALL_PVertex_in_ParamPolyhedron(V2, last, PP) {
4907 if (V == V2) {
4908 found = true;
4909 continue;
4910 }
4911 if (!found)
4912 continue;
4913 for (int i = 0; i < exist; ++i) {
4914 value_oppose(f, V->Vertex->p[i][nparam+1]);
4915 Vector_Combine(V->Vertex->p[i],
4916 V2->Vertex->p[i],
4917 M->p[0] + 1 + nvar + exist,
4918 V2->Vertex->p[i][nparam+1],
4919 f,
4920 nparam+1);
4921 int j;
4922 for (j = 0; j < nparam; ++j)
4923 if (value_notzero_p(M->p[0][1+nvar+exist+j]))
4924 break;
4925 if (j >= nparam)
4926 continue;
4927 ConstraintSimplify(M->p[0], M->p[0],
4928 P->Dimension+2, &f);
4929 value_set_si(M->p[0][0], 0);
4930 Polyhedron *para = AddConstraints(M->p[0], 1, P,
4931 MaxRays);
4932 if (emptyQ(para)) {
4933 Polyhedron_Free(para);
4934 continue;
4935 }
4936 Polyhedron *pos, *neg;
4937 value_set_si(M->p[0][0], 1);
4938 value_decrement(M->p[0][P->Dimension+1],
4939 M->p[0][P->Dimension+1]);
4940 neg = AddConstraints(M->p[0], 1, P, MaxRays);
4941 value_set_si(f, -1);
4942 Vector_Scale(M->p[0]+1, M->p[0]+1, f,
4943 P->Dimension+1);
4944 value_decrement(M->p[0][P->Dimension+1],
4945 M->p[0][P->Dimension+1]);
4946 value_decrement(M->p[0][P->Dimension+1],
4947 M->p[0][P->Dimension+1]);
4948 pos = AddConstraints(M->p[0], 1, P, MaxRays);
4949 if (emptyQ(neg) && emptyQ(pos)) {
4950 Polyhedron_Free(para);
4951 Polyhedron_Free(pos);
4952 Polyhedron_Free(neg);
4953 continue;
4954 }
4955 #ifdef DEBUG_ER
4956 fprintf(stderr, "\nER: Order\n");
4957 #endif /* DEBUG_ER */
4958 EP = barvinok_enumerate_e(para, exist, nparam, MaxRays);
4959 evalue *E;
4960 if (!emptyQ(pos)) {
4961 E = barvinok_enumerate_e(pos, exist, nparam, MaxRays);
4962 eadd(E, EP);
4963 free_evalue_refs(E);
4964 free(E);
4965 }
4966 if (!emptyQ(neg)) {
4967 E = barvinok_enumerate_e(neg, exist, nparam, MaxRays);
4968 eadd(E, EP);
4969 free_evalue_refs(E);
4970 free(E);
4971 }
4972 Polyhedron_Free(para);
4973 Polyhedron_Free(pos);
4974 Polyhedron_Free(neg);
4975 break;
4976 }
4977 if (EP)
4978 break;
4979 } END_FORALL_PVertex_in_ParamPolyhedron;
4980 if (EP)
4981 break;
4982 } END_FORALL_PVertex_in_ParamPolyhedron;
4983
4984 if (!EP) {
4985 /* Search for vertex coordinate to split on */
4986 /* First look for one independent of the parameters */
4987 FORALL_PVertex_in_ParamPolyhedron(V, last, PP) {
4988 for (int i = 0; i < exist; ++i) {
4989 int j;
4990 for (j = 0; j < nparam; ++j)
4991 if (value_notzero_p(V->Vertex->p[i][j]))
4992 break;
4993 if (j < nparam)
4994 continue;
4995 value_set_si(M->p[0][0], 1);
4996 Vector_Set(M->p[0]+1, 0, nvar+exist);
4997 Vector_Copy(V->Vertex->p[i],
4998 M->p[0] + 1 + nvar + exist, nparam+1);
4999 value_oppose(M->p[0][1+nvar+i],
5000 V->Vertex->p[i][nparam+1]);
5001
5002 Polyhedron *pos, *neg;
5003 value_set_si(M->p[0][0], 1);
5004 value_decrement(M->p[0][P->Dimension+1],
5005 M->p[0][P->Dimension+1]);
5006 neg = AddConstraints(M->p[0], 1, P, MaxRays);
5007 value_set_si(f, -1);
5008 Vector_Scale(M->p[0]+1, M->p[0]+1, f,
5009 P->Dimension+1);
5010 value_decrement(M->p[0][P->Dimension+1],
5011 M->p[0][P->Dimension+1]);
5012 value_decrement(M->p[0][P->Dimension+1],
5013 M->p[0][P->Dimension+1]);
5014 pos = AddConstraints(M->p[0], 1, P, MaxRays);
5015 if (emptyQ(neg) || emptyQ(pos)) {
5016 Polyhedron_Free(pos);
5017 Polyhedron_Free(neg);
5018 continue;
5019 }
5020 Polyhedron_Free(pos);
5021 value_increment(M->p[0][P->Dimension+1],
5022 M->p[0][P->Dimension+1]);
5023 pos = AddConstraints(M->p[0], 1, P, MaxRays);
5024 #ifdef DEBUG_ER
5025 fprintf(stderr, "\nER: Vertex\n");
5026 #endif /* DEBUG_ER */
5027 pos->next = neg;
5028 EP = enumerate_or(pos, exist, nparam, MaxRays);
5029 break;
5030 }
5031 if (EP)
5032 break;
5033 } END_FORALL_PVertex_in_ParamPolyhedron;
5034 }
5035
5036 if (!EP) {
5037 /* Search for vertex coordinate to split on */
5038 /* Now look for one that depends on the parameters */
5039 FORALL_PVertex_in_ParamPolyhedron(V, last, PP) {
5040 for (int i = 0; i < exist; ++i) {
5041 value_set_si(M->p[0][0], 1);
5042 Vector_Set(M->p[0]+1, 0, nvar+exist);
5043 Vector_Copy(V->Vertex->p[i],
5044 M->p[0] + 1 + nvar + exist, nparam+1);
5045 value_oppose(M->p[0][1+nvar+i],
5046 V->Vertex->p[i][nparam+1]);
5047
5048 Polyhedron *pos, *neg;
5049 value_set_si(M->p[0][0], 1);
5050 value_decrement(M->p[0][P->Dimension+1],
5051 M->p[0][P->Dimension+1]);
5052 neg = AddConstraints(M->p[0], 1, P, MaxRays);
5053 value_set_si(f, -1);
5054 Vector_Scale(M->p[0]+1, M->p[0]+1, f,
5055 P->Dimension+1);
5056 value_decrement(M->p[0][P->Dimension+1],
5057 M->p[0][P->Dimension+1]);
5058 value_decrement(M->p[0][P->Dimension+1],
5059 M->p[0][P->Dimension+1]);
5060 pos = AddConstraints(M->p[0], 1, P, MaxRays);
5061 if (emptyQ(neg) || emptyQ(pos)) {
5062 Polyhedron_Free(pos);
5063 Polyhedron_Free(neg);
5064 continue;
5065 }
5066 Polyhedron_Free(pos);
5067 value_increment(M->p[0][P->Dimension+1],
5068 M->p[0][P->Dimension+1]);
5069 pos = AddConstraints(M->p[0], 1, P, MaxRays);
5070 #ifdef DEBUG_ER
5071 fprintf(stderr, "\nER: ParamVertex\n");
5072 #endif /* DEBUG_ER */
5073 pos->next = neg;
5074 EP = enumerate_or(pos, exist, nparam, MaxRays);
5075 break;
5076 }
5077 if (EP)
5078 break;
5079 } END_FORALL_PVertex_in_ParamPolyhedron;
5080 }
5081
5082 Matrix_Free(M);
5083 value_clear(f);
5084 }
5085
5086 if (CEq)
5087 Polyhedron_Free(CEq);
5088 if (CT)
5089 Matrix_Free(CT);
5090 if (PP)
5091 Param_Polyhedron_Free(PP);
5092 *PA = P;
5093
5094 return EP;
5095 }
5096
5097 #ifndef HAVE_PIPLIB
5098 evalue *barvinok_enumerate_pip(Polyhedron *P,
5099 unsigned exist, unsigned nparam, unsigned MaxRays)
5100 {
5101 return 0;
5102 }
5103 #else
5104 evalue *barvinok_enumerate_pip(Polyhedron *P,
5105 unsigned exist, unsigned nparam, unsigned MaxRays)
5106 {
5107 int nvar = P->Dimension - exist - nparam;
5108 evalue *EP = new_zero_ep();
5109 Polyhedron *Q, *N, *T = 0;
5110 Value min, tmp;
5111 value_init(min);
5112 value_init(tmp);
5113
5114 #ifdef DEBUG_ER
5115 fprintf(stderr, "\nER: PIP\n");
5116 #endif /* DEBUG_ER */
5117
5118 for (int i = 0; i < P->Dimension; ++i) {
5119 bool pos = false;
5120 bool neg = false;
5121 bool posray = false;
5122 bool negray = false;
5123 value_set_si(min, 0);
5124 for (int j = 0; j < P->NbRays; ++j) {
5125 if (value_pos_p(P->Ray[j][1+i])) {
5126 pos = true;
5127 if (value_zero_p(P->Ray[j][1+P->Dimension]))
5128 posray = true;
5129 } else if (value_neg_p(P->Ray[j][1+i])) {
5130 neg = true;
5131 if (value_zero_p(P->Ray[j][1+P->Dimension]))
5132 negray = true;
5133 else {
5134 mpz_fdiv_q(tmp,
5135 P->Ray[j][1+i], P->Ray[j][1+P->Dimension]);
5136 if (value_lt(tmp, min))
5137 value_assign(min, tmp);
5138 }
5139 }
5140 }
5141 if (pos && neg) {
5142 assert(!(posray && negray));
5143 assert(!negray); // for now
5144 Polyhedron *O = T ? T : P;
5145 /* shift by a safe amount */
5146 Matrix *M = Matrix_Alloc(O->NbRays, O->Dimension+2);
5147 Vector_Copy(O->Ray[0], M->p[0], O->NbRays * (O->Dimension+2));
5148 for (int j = 0; j < P->NbRays; ++j) {
5149 if (value_notzero_p(M->p[j][1+P->Dimension])) {
5150 value_multiply(tmp, min, M->p[j][1+P->Dimension]);
5151 value_subtract(M->p[j][1+i], M->p[j][1+i], tmp);
5152 }
5153 }
5154 if (T)
5155 Polyhedron_Free(T);
5156 T = Rays2Polyhedron(M, MaxRays);
5157 Matrix_Free(M);
5158 } else if (neg) {
5159 /* negating a parameter requires that we substitute in the
5160 * sign again afterwards.
5161 * Disallow for now.
5162 */
5163 assert(i < nvar+exist);
5164 if (!T)
5165 T = Polyhedron_Copy(P);
5166 for (int j = 0; j < T->NbRays; ++j)
5167 value_oppose(T->Ray[j][1+i], T->Ray[j][1+i]);
5168 for (int j = 0; j < T->NbConstraints; ++j)
5169 value_oppose(T->Constraint[j][1+i], T->Constraint[j][1+i]);
5170 }
5171 }
5172 value_clear(min);
5173 value_clear(tmp);
5174
5175 Polyhedron *D = pip_lexmin(T ? T : P, exist, nparam);
5176 for (Q = D; Q; Q = N) {
5177 N = Q->next;
5178 Q->next = 0;
5179 evalue *E;
5180 exist = Q->Dimension - nvar - nparam;
5181 E = barvinok_enumerate_e(Q, exist, nparam, MaxRays);
5182 Polyhedron_Free(Q);
5183 eadd(E, EP);
5184 free_evalue_refs(E);
5185 free(E);
5186 }
5187
5188 if (T)
5189 Polyhedron_Free(T);
5190
5191 return EP;
5192 }
5193 #endif
5194
5195
5196 static bool is_single(Value *row, int pos, int len)
5197 {
5198 return First_Non_Zero(row, pos) == -1 &&
5199 First_Non_Zero(row+pos+1, len-pos-1) == -1;
5200 }
5201
5202 static evalue* barvinok_enumerate_e_r(Polyhedron *P,
5203 unsigned exist, unsigned nparam, unsigned MaxRays);
5204
5205 #ifdef DEBUG_ER
5206 static int er_level = 0;
5207
5208 evalue* barvinok_enumerate_e(Polyhedron *P,
5209 unsigned exist, unsigned nparam, unsigned MaxRays)
5210 {
5211 fprintf(stderr, "\nER: level %i\n", er_level);
5212 int nvar = P->Dimension - exist - nparam;
5213 fprintf(stderr, "%d %d %d\n", nvar, exist, nparam);
5214
5215 Polyhedron_Print(stderr, P_VALUE_FMT, P);
5216 ++er_level;
5217 P = DomainConstraintSimplify(Polyhedron_Copy(P), MaxRays);
5218 evalue *EP = barvinok_enumerate_e_r(P, exist, nparam, MaxRays);
5219 Polyhedron_Free(P);
5220 --er_level;
5221 return EP;
5222 }
5223 #else
5224 evalue* barvinok_enumerate_e(Polyhedron *P,
5225 unsigned exist, unsigned nparam, unsigned MaxRays)
5226 {
5227 P = DomainConstraintSimplify(Polyhedron_Copy(P), MaxRays);
5228 evalue *EP = barvinok_enumerate_e_r(P, exist, nparam, MaxRays);
5229 Polyhedron_Free(P);
5230 return EP;
5231 }
5232 #endif
5233
5234 static evalue* barvinok_enumerate_e_r(Polyhedron *P,
5235 unsigned exist, unsigned nparam, unsigned MaxRays)
5236 {
5237 if (exist == 0) {
5238 Polyhedron *U = Universe_Polyhedron(nparam);
5239 evalue *EP = barvinok_enumerate_ev(P, U, MaxRays);
5240 //char *param_name[] = {"P", "Q", "R", "S", "T" };
5241 //print_evalue(stdout, EP, param_name);
5242 Polyhedron_Free(U);
5243 return EP;
5244 }
5245
5246 int nvar = P->Dimension - exist - nparam;
5247 int len = P->Dimension + 2;
5248
5249 if (emptyQ(P))
5250 return new_zero_ep();
5251
5252 if (nvar == 0 && nparam == 0) {
5253 evalue *EP = new_zero_ep();
5254 barvinok_count(P, &EP->x.n, MaxRays);
5255 if (value_pos_p(EP->x.n))
5256 value_set_si(EP->x.n, 1);
5257 return EP;
5258 }
5259
5260 int r;
5261 for (r = 0; r < P->NbRays; ++r)
5262 if (value_zero_p(P->Ray[r][0]) ||
5263 value_zero_p(P->Ray[r][P->Dimension+1])) {
5264 int i;
5265 for (i = 0; i < nvar; ++i)
5266 if (value_notzero_p(P->Ray[r][i+1]))
5267 break;
5268 if (i >= nvar)
5269 continue;
5270 for (i = nvar + exist; i < nvar + exist + nparam; ++i)
5271 if (value_notzero_p(P->Ray[r][i+1]))
5272 break;
5273 if (i >= nvar + exist + nparam)
5274 break;
5275 }
5276 if (r < P->NbRays) {
5277 evalue *EP = new_zero_ep();
5278 value_set_si(EP->x.n, -1);
5279 return EP;
5280 }
5281
5282 int first;
5283 for (r = 0; r < P->NbEq; ++r)
5284 if ((first = First_Non_Zero(P->Constraint[r]+1+nvar, exist)) != -1)
5285 break;
5286 if (r < P->NbEq) {
5287 if (First_Non_Zero(P->Constraint[r]+1+nvar+first+1,
5288 exist-first-1) != -1) {
5289 Polyhedron *T = rotate_along(P, r, nvar, exist, MaxRays);
5290 #ifdef DEBUG_ER
5291 fprintf(stderr, "\nER: Equality\n");
5292 #endif /* DEBUG_ER */
5293 evalue *EP = barvinok_enumerate_e(T, exist-1, nparam, MaxRays);
5294 Polyhedron_Free(T);
5295 return EP;
5296 } else {
5297 #ifdef DEBUG_ER
5298 fprintf(stderr, "\nER: Fixed\n");
5299 #endif /* DEBUG_ER */
5300 if (first == 0)
5301 return barvinok_enumerate_e(P, exist-1, nparam, MaxRays);
5302 else {
5303 Polyhedron *T = Polyhedron_Copy(P);
5304 SwapColumns(T, nvar+1, nvar+1+first);
5305 evalue *EP = barvinok_enumerate_e(T, exist-1, nparam, MaxRays);
5306 Polyhedron_Free(T);
5307 return EP;
5308 }
5309 }
5310 }
5311
5312 Vector *row = Vector_Alloc(len);
5313 value_set_si(row->p[0], 1);
5314
5315 Value f;
5316 value_init(f);
5317
5318 enum constraint* info = new constraint[exist];
5319 for (int i = 0; i < exist; ++i) {
5320 info[i] = ALL_POS;
5321 for (int l = P->NbEq; l < P->NbConstraints; ++l) {
5322 if (value_negz_p(P->Constraint[l][nvar+i+1]))
5323 continue;
5324 bool l_parallel = is_single(P->Constraint[l]+nvar+1, i, exist);
5325 for (int u = P->NbEq; u < P->NbConstraints; ++u) {
5326 if (value_posz_p(P->Constraint[u][nvar+i+1]))
5327 continue;
5328 bool lu_parallel = l_parallel ||
5329 is_single(P->Constraint[u]+nvar+1, i, exist);
5330 value_oppose(f, P->Constraint[u][nvar+i+1]);
5331 Vector_Combine(P->Constraint[l]+1, P->Constraint[u]+1, row->p+1,
5332 f, P->Constraint[l][nvar+i+1], len-1);
5333 if (!(info[i] & INDEPENDENT)) {
5334 int j;
5335 for (j = 0; j < exist; ++j)
5336 if (j != i && value_notzero_p(row->p[nvar+j+1]))
5337 break;
5338 if (j == exist) {
5339 //printf("independent: i: %d, l: %d, u: %d\n", i, l, u);
5340 info[i] = (constraint)(info[i] | INDEPENDENT);
5341 }
5342 }
5343 if (info[i] & ALL_POS) {
5344 value_addto(row->p[len-1], row->p[len-1],
5345 P->Constraint[l][nvar+i+1]);
5346 value_addto(row->p[len-1], row->p[len-1], f);
5347 value_multiply(f, f, P->Constraint[l][nvar+i+1]);
5348 value_subtract(row->p[len-1], row->p[len-1], f);
5349 value_decrement(row->p[len-1], row->p[len-1]);
5350 ConstraintSimplify(row->p, row->p, len, &f);
5351 value_set_si(f, -1);
5352 Vector_Scale(row->p+1, row->p+1, f, len-1);
5353 value_decrement(row->p[len-1], row->p[len-1]);
5354 Polyhedron *T = AddConstraints(row->p, 1, P, MaxRays);
5355 if (!emptyQ(T)) {
5356 //printf("not all_pos: i: %d, l: %d, u: %d\n", i, l, u);
5357 info[i] = (constraint)(info[i] ^ ALL_POS);
5358 }
5359 //puts("pos remainder");
5360 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
5361 Polyhedron_Free(T);
5362 }
5363 if (!(info[i] & ONE_NEG)) {
5364 if (lu_parallel) {
5365 negative_test_constraint(P->Constraint[l],
5366 P->Constraint[u],
5367 row->p, nvar+i, len, &f);
5368 oppose_constraint(row->p, len, &f);
5369 Polyhedron *T = AddConstraints(row->p, 1, P, MaxRays);
5370 if (emptyQ(T)) {
5371 //printf("one_neg i: %d, l: %d, u: %d\n", i, l, u);
5372 info[i] = (constraint)(info[i] | ONE_NEG);
5373 }
5374 //puts("neg remainder");
5375 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
5376 Polyhedron_Free(T);
5377 } else if (!(info[i] & ROT_NEG)) {
5378 if (parallel_constraints(P->Constraint[l],
5379 P->Constraint[u],
5380 row->p, nvar, exist)) {
5381 negative_test_constraint7(P->Constraint[l],
5382 P->Constraint[u],
5383 row->p, nvar, exist,
5384 len, &f);
5385 oppose_constraint(row->p, len, &f);
5386 Polyhedron *T = AddConstraints(row->p, 1, P, MaxRays);
5387 if (emptyQ(T)) {
5388 // printf("rot_neg i: %d, l: %d, u: %d\n", i, l, u);
5389 info[i] = (constraint)(info[i] | ROT_NEG);
5390 r = l;
5391 }
5392 //puts("neg remainder");
5393 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
5394 Polyhedron_Free(T);
5395 }
5396 }
5397 }
5398 if (!(info[i] & ALL_POS) && (info[i] & (ONE_NEG | ROT_NEG)))
5399 goto next;
5400 }
5401 }
5402 if (info[i] & ALL_POS)
5403 break;
5404 next:
5405 ;
5406 }
5407
5408 /*
5409 for (int i = 0; i < exist; ++i)
5410 printf("%i: %i\n", i, info[i]);
5411 */
5412 for (int i = 0; i < exist; ++i)
5413 if (info[i] & ALL_POS) {
5414 #ifdef DEBUG_ER
5415 fprintf(stderr, "\nER: Positive\n");
5416 #endif /* DEBUG_ER */
5417 // Eliminate
5418 // Maybe we should chew off some of the fat here
5419 Matrix *M = Matrix_Alloc(P->Dimension, P->Dimension+1);
5420 for (int j = 0; j < P->Dimension; ++j)
5421 value_set_si(M->p[j][j + (j >= i+nvar)], 1);
5422 Polyhedron *T = Polyhedron_Image(P, M, MaxRays);
5423 Matrix_Free(M);
5424 evalue *EP = barvinok_enumerate_e(T, exist-1, nparam, MaxRays);
5425 Polyhedron_Free(T);
5426 value_clear(f);
5427 Vector_Free(row);
5428 delete [] info;
5429 return EP;
5430 }
5431 for (int i = 0; i < exist; ++i)
5432 if (info[i] & ONE_NEG) {
5433 #ifdef DEBUG_ER
5434 fprintf(stderr, "\nER: Negative\n");
5435 #endif /* DEBUG_ER */
5436 Vector_Free(row);
5437 value_clear(f);
5438 delete [] info;
5439 if (i == 0)
5440 return barvinok_enumerate_e(P, exist-1, nparam, MaxRays);
5441 else {
5442 Polyhedron *T = Polyhedron_Copy(P);
5443 SwapColumns(T, nvar+1, nvar+1+i);
5444 evalue *EP = barvinok_enumerate_e(T, exist-1, nparam, MaxRays);
5445 Polyhedron_Free(T);
5446 return EP;
5447 }
5448 }
5449 for (int i = 0; i < exist; ++i)
5450 if (info[i] & ROT_NEG) {
5451 #ifdef DEBUG_ER
5452 fprintf(stderr, "\nER: Rotate\n");
5453 #endif /* DEBUG_ER */
5454 Vector_Free(row);
5455 value_clear(f);
5456 delete [] info;
5457 Polyhedron *T = rotate_along(P, r, nvar, exist, MaxRays);
5458 evalue *EP = barvinok_enumerate_e(T, exist-1, nparam, MaxRays);
5459 Polyhedron_Free(T);
5460 return EP;
5461 }
5462 for (int i = 0; i < exist; ++i)
5463 if (info[i] & INDEPENDENT) {
5464 Polyhedron *pos, *neg;
5465
5466 /* Find constraint again and split off negative part */
5467
5468 if (SplitOnVar(P, i, nvar, len, exist, MaxRays,
5469 row, f, true, &pos, &neg)) {
5470 #ifdef DEBUG_ER
5471 fprintf(stderr, "\nER: Split\n");
5472 #endif /* DEBUG_ER */
5473
5474 evalue *EP =
5475 barvinok_enumerate_e(neg, exist-1, nparam, MaxRays);
5476 evalue *E =
5477 barvinok_enumerate_e(pos, exist, nparam, MaxRays);
5478 eadd(E, EP);
5479 free_evalue_refs(E);
5480 free(E);
5481 Polyhedron_Free(neg);
5482 Polyhedron_Free(pos);
5483 value_clear(f);
5484 Vector_Free(row);
5485 delete [] info;
5486 return EP;
5487 }
5488 }
5489 delete [] info;
5490
5491 Polyhedron *O = P;
5492 Polyhedron *F;
5493
5494 evalue *EP;
5495
5496 EP = enumerate_line(P, exist, nparam, MaxRays);
5497 if (EP)
5498 goto out;
5499
5500 EP = barvinok_enumerate_pip(P, exist, nparam, MaxRays);
5501 if (EP)
5502 goto out;
5503
5504 EP = enumerate_redundant_ray(P, exist, nparam, MaxRays);
5505 if (EP)
5506 goto out;
5507
5508 EP = enumerate_sure(P, exist, nparam, MaxRays);
5509 if (EP)
5510 goto out;
5511
5512 EP = enumerate_ray(P, exist, nparam, MaxRays);
5513 if (EP)
5514 goto out;
5515
5516 EP = enumerate_sure2(P, exist, nparam, MaxRays);
5517 if (EP)
5518 goto out;
5519
5520 F = unfringe(P, MaxRays);
5521 if (!PolyhedronIncludes(F, P)) {
5522 #ifdef DEBUG_ER
5523 fprintf(stderr, "\nER: Fringed\n");
5524 #endif /* DEBUG_ER */
5525 EP = barvinok_enumerate_e(F, exist, nparam, MaxRays);
5526 Polyhedron_Free(F);
5527 goto out;
5528 }
5529 Polyhedron_Free(F);
5530
5531 if (nparam)
5532 EP = enumerate_vd(&P, exist, nparam, MaxRays);
5533 if (EP)
5534 goto out2;
5535
5536 if (nvar != 0) {
5537 EP = enumerate_sum(P, exist, nparam, MaxRays);
5538 goto out2;
5539 }
5540
5541 assert(nvar == 0);
5542
5543 int i;
5544 Polyhedron *pos, *neg;
5545 for (i = 0; i < exist; ++i)
5546 if (SplitOnVar(P, i, nvar, len, exist, MaxRays,
5547 row, f, false, &pos, &neg))
5548 break;
5549
5550 assert (i < exist);
5551
5552 pos->next = neg;
5553 EP = enumerate_or(pos, exist, nparam, MaxRays);
5554
5555 out2:
5556 if (O != P)
5557 Polyhedron_Free(P);
5558
5559 out:
5560 value_clear(f);
5561 Vector_Free(row);
5562 return EP;
5563 }
5564
5565 gen_fun * barvinok_series(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
5566 {
5567 Polyhedron *CA;
5568 unsigned nparam = C->Dimension;
5569 unsigned dim, nvar;
5570 vec_ZZ sign;
5571 int ncone = 0;
5572 sign.SetLength(ncone);
5573
5574 CA = align_context(C, P->Dimension, MaxRays);
5575 P = DomainIntersection(P, CA, MaxRays);
5576 Polyhedron_Free(CA);
5577
5578 assert(!Polyhedron_is_infinite(P, nparam));
5579 assert(P->NbBid == 0);
5580 assert(Polyhedron_has_positive_rays(P, nparam));
5581 if (P->NbEq != 0)
5582 P = remove_equalities_p(P, P->Dimension-nparam, NULL);
5583 assert(P->NbEq == 0);
5584
5585 #ifdef USE_INCREMENTAL_BF
5586 partial_bfcounter red(P, nparam);
5587 #elif defined USE_INCREMENTAL_DF
5588 partial_ireducer red(P, nparam);
5589 #else
5590 partial_reducer red(P, nparam);
5591 #endif
5592 red.start(MaxRays);
5593 Polyhedron_Free(P);
5594 return red.gf;
5595 }
lattice point enumeration library