?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
barvinok_enumerate_e.cc: memory clean up.
[barvinok.git] / barvinok.cc
blobfd9f8c78225599c1597c061e7115e44ae59b6ec4
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 <barvinok/util.h>
11 extern "C" {
12 #include <polylib/polylibgmp.h>
13 #include <barvinok/evalue.h>
14 #include "piputil.h"
15 }
16 #include "config.h"
17 #include <barvinok/barvinok.h>
18 #include <barvinok/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 = triangulate_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 Domain_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 /* this procedure may have false negatives */
2849 static bool Polyhedron_is_infinite(Polyhedron *P, unsigned nparam)
2850 {
2851 int r;
2852 for (r = 0; r < P->NbRays; ++r) {
2853 if (!value_zero_p(P->Ray[r][0]) &&
2854 !value_zero_p(P->Ray[r][P->Dimension+1]))
2855 continue;
2856 if (First_Non_Zero(P->Ray[r]+1+P->Dimension-nparam, nparam) == -1)
2857 return true;
2858 }
2859 return false;
2860 }
2861
2862 /* Check whether all rays point in the positive directions
2863 * for the parameters
2864 */
2865 static bool Polyhedron_has_positive_rays(Polyhedron *P, unsigned nparam)
2866 {
2867 int r;
2868 for (r = 0; r < P->NbRays; ++r)
2869 if (value_zero_p(P->Ray[r][P->Dimension+1])) {
2870 int i;
2871 for (i = P->Dimension - nparam; i < P->Dimension; ++i)
2872 if (value_neg_p(P->Ray[r][i+1]))
2873 return false;
2874 }
2875 return true;
2876 }
2877
2878 typedef evalue * evalue_p;
2879
2880 struct enumerator : public polar_decomposer {
2881 vec_ZZ lambda;
2882 unsigned dim, nbV;
2883 evalue ** vE;
2884 int _i;
2885 mat_ZZ rays;
2886 vec_ZZ den;
2887 ZZ sign;
2888 Polyhedron *P;
2889 Param_Vertices *V;
2890 term_info num;
2891 Vector *c;
2892 mpq_t count;
2893
2894 enumerator(Polyhedron *P, unsigned dim, unsigned nbV) {
2895 this->P = P;
2896 this->dim = dim;
2897 this->nbV = nbV;
2898 randomvector(P, lambda, dim);
2899 rays.SetDims(dim, dim);
2900 den.SetLength(dim);
2901 c = Vector_Alloc(dim+2);
2902
2903 vE = new evalue_p[nbV];
2904 for (int j = 0; j < nbV; ++j)
2905 vE[j] = 0;
2906
2907 mpq_init(count);
2908 }
2909
2910 void decompose_at(Param_Vertices *V, int _i, unsigned MaxRays) {
2911 Polyhedron *C = supporting_cone_p(P, V);
2912 this->_i = _i;
2913 this->V = V;
2914
2915 vE[_i] = new evalue;
2916 value_init(vE[_i]->d);
2917 evalue_set_si(vE[_i], 0, 1);
2918
2919 decompose(C, MaxRays);
2920 }
2921
2922 ~enumerator() {
2923 mpq_clear(count);
2924 Vector_Free(c);
2925
2926 for (int j = 0; j < nbV; ++j)
2927 if (vE[j]) {
2928 free_evalue_refs(vE[j]);
2929 delete vE[j];
2930 }
2931 delete [] vE;
2932 }
2933
2934 virtual void handle_polar(Polyhedron *P, int sign);
2935 };
2936
2937 void enumerator::handle_polar(Polyhedron *C, int s)
2938 {
2939 int r = 0;
2940 assert(C->NbRays-1 == dim);
2941 add_rays(rays, C, &r);
2942 for (int k = 0; k < dim; ++k) {
2943 if (lambda * rays[k] == 0)
2944 throw Orthogonal;
2945 }
2946
2947 sign = s;
2948
2949 lattice_point(V, C, lambda, &num, 0);
2950 den = rays * lambda;
2951 normalize(sign, num.constant, den);
2952
2953 dpoly n(dim, den[0], 1);
2954 for (int k = 1; k < dim; ++k) {
2955 dpoly fact(dim, den[k], 1);
2956 n *= fact;
2957 }
2958 if (num.E != NULL) {
2959 ZZ one(INIT_VAL, 1);
2960 dpoly_n d(dim, num.constant, one);
2961 d.div(n, c, sign);
2962 evalue EV;
2963 multi_polynom(c, num.E, &EV);
2964 eadd(&EV , vE[_i]);
2965 free_evalue_refs(&EV);
2966 free_evalue_refs(num.E);
2967 delete num.E;
2968 } else if (num.pos != -1) {
2969 dpoly_n d(dim, num.constant, num.coeff);
2970 d.div(n, c, sign);
2971 evalue EV;
2972 uni_polynom(num.pos, c, &EV);
2973 eadd(&EV , vE[_i]);
2974 free_evalue_refs(&EV);
2975 } else {
2976 mpq_set_si(count, 0, 1);
2977 dpoly d(dim, num.constant);
2978 d.div(n, count, sign);
2979 evalue EV;
2980 value_init(EV.d);
2981 evalue_set(&EV, &count[0]._mp_num, &count[0]._mp_den);
2982 eadd(&EV , vE[_i]);
2983 free_evalue_refs(&EV);
2984 }
2985 }
2986
2987 struct enumerator_base : public virtual polar_decomposer {
2988 Polyhedron *P;
2989 unsigned dim, nbV;
2990 evalue ** vE;
2991 int _i;
2992 Param_Vertices *V;
2993 evalue ** E_vertex;
2994 evalue mone;
2995
2996 enumerator_base(Polyhedron *P, unsigned dim, unsigned nbV) {
2997 this->P = P;
2998 this->dim = dim;
2999 this->nbV = nbV;
3000
3001 vE = new evalue_p[nbV];
3002 for (int j = 0; j < nbV; ++j)
3003 vE[j] = 0;
3004
3005 E_vertex = new evalue_p[dim];
3006
3007 value_init(mone.d);
3008 evalue_set_si(&mone, -1, 1);
3009 }
3010
3011 void decompose_at(Param_Vertices *V, int _i, unsigned MaxRays/*, Polyhedron *pVD*/) {
3012 Polyhedron *C = supporting_cone_p(P, V);
3013 this->_i = _i;
3014 this->V = V;
3015 //this->pVD = pVD;
3016
3017 vE[_i] = new evalue;
3018 value_init(vE[_i]->d);
3019 evalue_set_si(vE[_i], 0, 1);
3020
3021 decompose(C, MaxRays);
3022 }
3023
3024 ~enumerator_base() {
3025 for (int j = 0; j < nbV; ++j)
3026 if (vE[j]) {
3027 free_evalue_refs(vE[j]);
3028 delete vE[j];
3029 }
3030 delete [] vE;
3031
3032 delete [] E_vertex;
3033
3034 free_evalue_refs(&mone);
3035 }
3036
3037 evalue *E_num(int i, int d) {
3038 return E_vertex[i + (dim-d)];
3039 }
3040 };
3041
3042 struct cumulator {
3043 evalue *factor;
3044 evalue *v;
3045 dpoly_r *r;
3046
3047 cumulator(evalue *factor, evalue *v, dpoly_r *r) :
3048 factor(factor), v(v), r(r) {}
3049
3050 void cumulate();
3051
3052 virtual void add_term(int *powers, int len, evalue *f2) = 0;
3053 };
3054
3055 void cumulator::cumulate()
3056 {
3057 evalue cum; // factor * 1 * E_num[0]/1 * (E_num[0]-1)/2 *...
3058 evalue f;
3059 evalue t; // E_num[0] - (m-1)
3060 #ifdef USE_MODULO
3061 evalue *cst;
3062 #else
3063 evalue mone;
3064 value_init(mone.d);
3065 evalue_set_si(&mone, -1, 1);
3066 #endif
3067
3068 value_init(cum.d);
3069 evalue_copy(&cum, factor);
3070 value_init(f.d);
3071 value_init(f.x.n);
3072 value_set_si(f.d, 1);
3073 value_set_si(f.x.n, 1);
3074 value_init(t.d);
3075 evalue_copy(&t, v);
3076
3077 #ifdef USE_MODULO
3078 for (cst = &t; value_zero_p(cst->d); ) {
3079 if (cst->x.p->type == fractional)
3080 cst = &cst->x.p->arr[1];
3081 else
3082 cst = &cst->x.p->arr[0];
3083 }
3084 #endif
3085
3086 for (int m = 0; m < r->len; ++m) {
3087 if (m > 0) {
3088 if (m > 1) {
3089 value_set_si(f.d, m);
3090 emul(&f, &cum);
3091 #ifdef USE_MODULO
3092 value_subtract(cst->x.n, cst->x.n, cst->d);
3093 #else
3094 eadd(&mone, &t);
3095 #endif
3096 }
3097 emul(&t, &cum);
3098 }
3099 vector< dpoly_r_term * >& current = r->c[r->len-1-m];
3100 for (int j = 0; j < current.size(); ++j) {
3101 if (current[j]->coeff == 0)
3102 continue;
3103 evalue *f2 = new evalue;
3104 value_init(f2->d);
3105 value_init(f2->x.n);
3106 zz2value(current[j]->coeff, f2->x.n);
3107 zz2value(r->denom, f2->d);
3108 emul(&cum, f2);
3109
3110 add_term(current[j]->powers, r->dim, f2);
3111 }
3112 }
3113 free_evalue_refs(&f);
3114 free_evalue_refs(&t);
3115 free_evalue_refs(&cum);
3116 #ifndef USE_MODULO
3117 free_evalue_refs(&mone);
3118 #endif
3119 }
3120
3121 struct E_poly_term {
3122 int *powers;
3123 evalue *E;
3124 };
3125
3126 struct ie_cum : public cumulator {
3127 vector<E_poly_term *> terms;
3128
3129 ie_cum(evalue *factor, evalue *v, dpoly_r *r) : cumulator(factor, v, r) {}
3130
3131 virtual void add_term(int *powers, int len, evalue *f2);
3132 };
3133
3134 void ie_cum::add_term(int *powers, int len, evalue *f2)
3135 {
3136 int k;
3137 for (k = 0; k < terms.size(); ++k) {
3138 if (memcmp(terms[k]->powers, powers, len * sizeof(int)) == 0) {
3139 eadd(f2, terms[k]->E);
3140 free_evalue_refs(f2);
3141 delete f2;
3142 break;
3143 }
3144 }
3145 if (k >= terms.size()) {
3146 E_poly_term *ET = new E_poly_term;
3147 ET->powers = new int[len];
3148 memcpy(ET->powers, powers, len * sizeof(int));
3149 ET->E = f2;
3150 terms.push_back(ET);
3151 }
3152 }
3153
3154 struct ienumerator : public virtual polar_decomposer, public enumerator_base {
3155 //Polyhedron *pVD;
3156 mat_ZZ den;
3157 vec_ZZ vertex;
3158 mpq_t tcount;
3159
3160 ienumerator(Polyhedron *P, unsigned dim, unsigned nbV) :
3161 enumerator_base(P, dim, nbV) {
3162 vertex.SetLength(dim);
3163
3164 den.SetDims(dim, dim);
3165 mpq_init(tcount);
3166 }
3167
3168 ~ienumerator() {
3169 mpq_clear(tcount);
3170 }
3171
3172 virtual void handle_polar(Polyhedron *P, int sign);
3173 void reduce(evalue *factor, vec_ZZ& num, mat_ZZ& den_f);
3174 };
3175
3176 static evalue* new_zero_ep()
3177 {
3178 evalue *EP;
3179 ALLOC(evalue, EP);
3180 value_init(EP->d);
3181 evalue_set_si(EP, 0, 1);
3182 return EP;
3183 }
3184
3185 void lattice_point(Param_Vertices *V, Polyhedron *C, vec_ZZ& num,
3186 evalue **E_vertex)
3187 {
3188 unsigned nparam = V->Vertex->NbColumns - 2;
3189 unsigned dim = C->Dimension;
3190 vec_ZZ vertex;
3191 vertex.SetLength(nparam+1);
3192
3193 Value lcm, tmp;
3194 value_init(lcm);
3195 value_init(tmp);
3196 value_set_si(lcm, 1);
3197
3198 for (int j = 0; j < V->Vertex->NbRows; ++j) {
3199 value_lcm(lcm, V->Vertex->p[j][nparam+1], &lcm);
3200 }
3201
3202 if (value_notone_p(lcm)) {
3203 Matrix * mv = Matrix_Alloc(dim, nparam+1);
3204 for (int j = 0 ; j < dim; ++j) {
3205 value_division(tmp, lcm, V->Vertex->p[j][nparam+1]);
3206 Vector_Scale(V->Vertex->p[j], mv->p[j], tmp, nparam+1);
3207 }
3208
3209 Matrix* Rays = rays2(C);
3210 Matrix *T = Transpose(Rays);
3211 Matrix *T2 = Matrix_Copy(T);
3212 Matrix *inv = Matrix_Alloc(T2->NbRows, T2->NbColumns);
3213 int ok = Matrix_Inverse(T2, inv);
3214 assert(ok);
3215 Matrix_Free(Rays);
3216 Matrix_Free(T2);
3217 Matrix *L = Matrix_Alloc(inv->NbRows, mv->NbColumns);
3218 Matrix_Product(inv, mv, L);
3219 Matrix_Free(inv);
3220
3221 evalue f;
3222 value_init(f.d);
3223 value_init(f.x.n);
3224
3225 ZZ one;
3226
3227 evalue *remainders[dim];
3228 for (int i = 0; i < dim; ++i) {
3229 remainders[i] = new_zero_ep();
3230 one = 1;
3231 ceil(L->p[i], nparam+1, lcm, one, remainders[i], 0);
3232 }
3233 Matrix_Free(L);
3234
3235
3236 for (int i = 0; i < V->Vertex->NbRows; ++i) {
3237 values2zz(mv->p[i], vertex, nparam+1);
3238 E_vertex[i] = multi_monom(vertex);
3239 num[i] = 0;
3240
3241 value_set_si(f.x.n, 1);
3242 value_assign(f.d, lcm);
3243
3244 emul(&f, E_vertex[i]);
3245
3246 for (int j = 0; j < dim; ++j) {
3247 if (value_zero_p(T->p[i][j]))
3248 continue;
3249 evalue cp;
3250 value_init(cp.d);
3251 evalue_copy(&cp, remainders[j]);
3252 if (value_notone_p(T->p[i][j])) {
3253 value_set_si(f.d, 1);
3254 value_assign(f.x.n, T->p[i][j]);
3255 emul(&f, &cp);
3256 }
3257 eadd(&cp, E_vertex[i]);
3258 free_evalue_refs(&cp);
3259 }
3260 }
3261 for (int i = 0; i < dim; ++i) {
3262 free_evalue_refs(remainders[i]);
3263 free(remainders[i]);
3264 }
3265
3266 free_evalue_refs(&f);
3267
3268 Matrix_Free(T);
3269 Matrix_Free(mv);
3270 value_clear(lcm);
3271 value_clear(tmp);
3272 return;
3273 }
3274 value_clear(lcm);
3275 value_clear(tmp);
3276
3277 for (int i = 0; i < V->Vertex->NbRows; ++i) {
3278 /* fixed value */
3279 if (First_Non_Zero(V->Vertex->p[i], nparam) == -1) {
3280 E_vertex[i] = 0;
3281 value2zz(V->Vertex->p[i][nparam], num[i]);
3282 } else {
3283 values2zz(V->Vertex->p[i], vertex, nparam+1);
3284 E_vertex[i] = multi_monom(vertex);
3285 num[i] = 0;
3286 }
3287 }
3288 }
3289
3290 void ienumerator::reduce(
3291 evalue *factor, vec_ZZ& num, mat_ZZ& den_f)
3292 {
3293 unsigned len = den_f.NumRows(); // number of factors in den
3294 unsigned dim = num.length();
3295
3296 if (dim == 0) {
3297 eadd(factor, vE[_i]);
3298 return;
3299 }
3300
3301 vec_ZZ den_s;
3302 den_s.SetLength(len);
3303 mat_ZZ den_r;
3304 den_r.SetDims(len, dim-1);
3305
3306 int r, k;
3307
3308 for (r = 0; r < len; ++r) {
3309 den_s[r] = den_f[r][0];
3310 for (k = 0; k <= dim-1; ++k)
3311 if (k != 0)
3312 den_r[r][k-(k>0)] = den_f[r][k];
3313 }
3314
3315 ZZ num_s = num[0];
3316 vec_ZZ num_p;
3317 num_p.SetLength(dim-1);
3318 for (k = 0 ; k <= dim-1; ++k)
3319 if (k != 0)
3320 num_p[k-(k>0)] = num[k];
3321
3322 vec_ZZ den_p;
3323 den_p.SetLength(len);
3324
3325 ZZ one;
3326 one = 1;
3327 normalize(one, num_s, num_p, den_s, den_p, den_r);
3328 if (one != 1)
3329 emul(&mone, factor);
3330
3331 int only_param = 0;
3332 int no_param = 0;
3333 for (int k = 0; k < len; ++k) {
3334 if (den_p[k] == 0)
3335 ++no_param;
3336 else if (den_s[k] == 0)
3337 ++only_param;
3338 }
3339 if (no_param == 0) {
3340 reduce(factor, num_p, den_r);
3341 } else {
3342 int k, l;
3343 mat_ZZ pden;
3344 pden.SetDims(only_param, dim-1);
3345
3346 for (k = 0, l = 0; k < len; ++k)
3347 if (den_s[k] == 0)
3348 pden[l++] = den_r[k];
3349
3350 for (k = 0; k < len; ++k)
3351 if (den_p[k] == 0)
3352 break;
3353
3354 dpoly n(no_param, num_s);
3355 dpoly D(no_param, den_s[k], 1);
3356 for ( ; ++k < len; )
3357 if (den_p[k] == 0) {
3358 dpoly fact(no_param, den_s[k], 1);
3359 D *= fact;
3360 }
3361
3362 dpoly_r * r = 0;
3363 // if no_param + only_param == len then all powers
3364 // below will be all zero
3365 if (no_param + only_param == len) {
3366 if (E_num(0, dim) != 0)
3367 r = new dpoly_r(n, len);
3368 else {
3369 mpq_set_si(tcount, 0, 1);
3370 one = 1;
3371 n.div(D, tcount, one);
3372
3373 if (value_notzero_p(mpq_numref(tcount))) {
3374 evalue f;
3375 value_init(f.d);
3376 value_init(f.x.n);
3377 value_assign(f.x.n, mpq_numref(tcount));
3378 value_assign(f.d, mpq_denref(tcount));
3379 emul(&f, factor);
3380 reduce(factor, num_p, pden);
3381 free_evalue_refs(&f);
3382 }
3383 return;
3384 }
3385 } else {
3386 for (k = 0; k < len; ++k) {
3387 if (den_s[k] == 0 || den_p[k] == 0)
3388 continue;
3389
3390 dpoly pd(no_param-1, den_s[k], 1);
3391
3392 int l;
3393 for (l = 0; l < k; ++l)
3394 if (den_r[l] == den_r[k])
3395 break;
3396
3397 if (r == 0)
3398 r = new dpoly_r(n, pd, l, len);
3399 else {
3400 dpoly_r *nr = new dpoly_r(r, pd, l, len);
3401 delete r;
3402 r = nr;
3403 }
3404 }
3405 }
3406 dpoly_r *rc = r->div(D);
3407 delete r;
3408 r = rc;
3409 if (E_num(0, dim) == 0) {
3410 int common = pden.NumRows();
3411 vector< dpoly_r_term * >& final = r->c[r->len-1];
3412 int rows;
3413 evalue t;
3414 evalue f;
3415 value_init(f.d);
3416 value_init(f.x.n);
3417 zz2value(r->denom, f.d);
3418 for (int j = 0; j < final.size(); ++j) {
3419 if (final[j]->coeff == 0)
3420 continue;
3421 rows = common;
3422 for (int k = 0; k < r->dim; ++k) {
3423 int n = final[j]->powers[k];
3424 if (n == 0)
3425 continue;
3426 pden.SetDims(rows+n, pden.NumCols());
3427 for (int l = 0; l < n; ++l)
3428 pden[rows+l] = den_r[k];
3429 rows += n;
3430 }
3431 value_init(t.d);
3432 evalue_copy(&t, factor);
3433 zz2value(final[j]->coeff, f.x.n);
3434 emul(&f, &t);
3435 reduce(&t, num_p, pden);
3436 free_evalue_refs(&t);
3437 }
3438 free_evalue_refs(&f);
3439 } else {
3440 ie_cum cum(factor, E_num(0, dim), r);
3441 cum.cumulate();
3442
3443 int common = pden.NumRows();
3444 int rows;
3445 for (int j = 0; j < cum.terms.size(); ++j) {
3446 rows = common;
3447 pden.SetDims(rows, pden.NumCols());
3448 for (int k = 0; k < r->dim; ++k) {
3449 int n = cum.terms[j]->powers[k];
3450 if (n == 0)
3451 continue;
3452 pden.SetDims(rows+n, pden.NumCols());
3453 for (int l = 0; l < n; ++l)
3454 pden[rows+l] = den_r[k];
3455 rows += n;
3456 }
3457 reduce(cum.terms[j]->E, num_p, pden);
3458 free_evalue_refs(cum.terms[j]->E);
3459 delete cum.terms[j]->E;
3460 delete [] cum.terms[j]->powers;
3461 delete cum.terms[j];
3462 }
3463 }
3464 delete r;
3465 }
3466 }
3467
3468 static int type_offset(enode *p)
3469 {
3470 return p->type == fractional ? 1 :
3471 p->type == flooring ? 1 : 0;
3472 }
3473
3474 static int edegree(evalue *e)
3475 {
3476 int d = 0;
3477 enode *p;
3478
3479 if (value_notzero_p(e->d))
3480 return 0;
3481
3482 p = e->x.p;
3483 int i = type_offset(p);
3484 if (p->size-i-1 > d)
3485 d = p->size - i - 1;
3486 for (; i < p->size; i++) {
3487 int d2 = edegree(&p->arr[i]);
3488 if (d2 > d)
3489 d = d2;
3490 }
3491 return d;
3492 }
3493
3494 void ienumerator::handle_polar(Polyhedron *C, int s)
3495 {
3496 assert(C->NbRays-1 == dim);
3497
3498 lattice_point(V, C, vertex, E_vertex);
3499
3500 int r;
3501 for (r = 0; r < dim; ++r)
3502 values2zz(C->Ray[r]+1, den[r], dim);
3503
3504 evalue one;
3505 value_init(one.d);
3506 evalue_set_si(&one, s, 1);
3507 reduce(&one, vertex, den);
3508 free_evalue_refs(&one);
3509
3510 for (int i = 0; i < dim; ++i)
3511 if (E_vertex[i]) {
3512 free_evalue_refs(E_vertex[i]);
3513 delete E_vertex[i];
3514 }
3515
3516 /*
3517 {
3518 char * test[] = {"a", "b"};
3519 evalue E;
3520 value_init(E.d);
3521 evalue_copy(&E, vE[_i]);
3522 frac2floor_in_domain(&E, pVD);
3523 printf("***** Curr value:");
3524 print_evalue(stdout, &E, test);
3525 fprintf(stdout, "\n");
3526 }
3527 */
3528
3529 }
3530
3531 struct bfenumerator : public bf_base, public enumerator_base {
3532 evalue *factor;
3533
3534 bfenumerator(Polyhedron *P, unsigned dim, unsigned nbV) :
3535 bf_base(P, dim), enumerator_base(P, dim, nbV) {
3536 lower = 0;
3537 factor = NULL;
3538 }
3539
3540 ~bfenumerator() {
3541 }
3542
3543 virtual void handle_polar(Polyhedron *P, int sign);
3544 virtual void base(mat_ZZ& factors, bfc_vec& v);
3545
3546 bfc_term_base* new_bf_term(int len) {
3547 bfe_term* t = new bfe_term(len);
3548 return t;
3549 }
3550
3551 virtual void set_factor(bfc_term_base *t, int k, int change) {
3552 bfe_term* bfet = static_cast<bfe_term *>(t);
3553 factor = bfet->factors[k];
3554 assert(factor != NULL);
3555 bfet->factors[k] = NULL;
3556 if (change)
3557 emul(&mone, factor);
3558 }
3559
3560 virtual void set_factor(bfc_term_base *t, int k, mpq_t &q, int change) {
3561 bfe_term* bfet = static_cast<bfe_term *>(t);
3562 factor = bfet->factors[k];
3563 assert(factor != NULL);
3564 bfet->factors[k] = NULL;
3565
3566 evalue f;
3567 value_init(f.d);
3568 value_init(f.x.n);
3569 if (change)
3570 value_oppose(f.x.n, mpq_numref(q));
3571 else
3572 value_assign(f.x.n, mpq_numref(q));
3573 value_assign(f.d, mpq_denref(q));
3574 emul(&f, factor);
3575 }
3576
3577 virtual void set_factor(bfc_term_base *t, int k, ZZ& n, ZZ& d, int change) {
3578 bfe_term* bfet = static_cast<bfe_term *>(t);
3579
3580 factor = new evalue;
3581
3582 evalue f;
3583 value_init(f.d);
3584 value_init(f.x.n);
3585 zz2value(n, f.x.n);
3586 if (change)
3587 value_oppose(f.x.n, f.x.n);
3588 zz2value(d, f.d);
3589
3590 value_init(factor->d);
3591 evalue_copy(factor, bfet->factors[k]);
3592 emul(&f, factor);
3593 }
3594
3595 void set_factor(evalue *f, int change) {
3596 if (change)
3597 emul(&mone, f);
3598 factor = f;
3599 }
3600
3601 virtual void insert_term(bfc_term_base *t, int i) {
3602 bfe_term* bfet = static_cast<bfe_term *>(t);
3603 int len = t->terms.NumRows()-1; // already increased by one
3604
3605 bfet->factors.resize(len+1);
3606 for (int j = len; j > i; --j) {
3607 bfet->factors[j] = bfet->factors[j-1];
3608 t->terms[j] = t->terms[j-1];
3609 }
3610 bfet->factors[i] = factor;
3611 factor = NULL;
3612 }
3613
3614 virtual void update_term(bfc_term_base *t, int i) {
3615 bfe_term* bfet = static_cast<bfe_term *>(t);
3616
3617 eadd(factor, bfet->factors[i]);
3618 free_evalue_refs(factor);
3619 delete factor;
3620 }
3621
3622 virtual bool constant_vertex(int dim) { return E_num(0, dim) == 0; }
3623
3624 virtual void cum(bf_reducer *bfr, bfc_term_base *t, int k, dpoly_r *r);
3625 };
3626
3627 struct bfe_cum : public cumulator {
3628 bfenumerator *bfe;
3629 bfc_term_base *told;
3630 int k;
3631 bf_reducer *bfr;
3632
3633 bfe_cum(evalue *factor, evalue *v, dpoly_r *r, bf_reducer *bfr,
3634 bfc_term_base *t, int k, bfenumerator *e) :
3635 cumulator(factor, v, r), told(t), k(k),
3636 bfr(bfr), bfe(e) {
3637 }
3638
3639 virtual void add_term(int *powers, int len, evalue *f2);
3640 };
3641
3642 void bfe_cum::add_term(int *powers, int len, evalue *f2)
3643 {
3644 bfr->update_powers(powers, len);
3645
3646 bfc_term_base * t = bfe->find_bfc_term(bfr->vn, bfr->npowers, bfr->nnf);
3647 bfe->set_factor(f2, bfr->l_changes % 2);
3648 bfe->add_term(t, told->terms[k], bfr->l_extra_num);
3649 }
3650
3651 void bfenumerator::cum(bf_reducer *bfr, bfc_term_base *t, int k,
3652 dpoly_r *r)
3653 {
3654 bfe_term* bfet = static_cast<bfe_term *>(t);
3655 bfe_cum cum(bfet->factors[k], E_num(0, bfr->d), r, bfr, t, k, this);
3656 cum.cumulate();
3657 }
3658
3659 void bfenumerator::base(mat_ZZ& factors, bfc_vec& v)
3660 {
3661 for (int i = 0; i < v.size(); ++i) {
3662 assert(v[i]->terms.NumRows() == 1);
3663 evalue *factor = static_cast<bfe_term *>(v[i])->factors[0];
3664 eadd(factor, vE[_i]);
3665 delete v[i];
3666 }
3667 }
3668
3669 void bfenumerator::handle_polar(Polyhedron *C, int s)
3670 {
3671 assert(C->NbRays-1 == enumerator_base::dim);
3672
3673 bfe_term* t = new bfe_term(enumerator_base::dim);
3674 vector< bfc_term_base * > v;
3675 v.push_back(t);
3676
3677 t->factors.resize(1);
3678
3679 t->terms.SetDims(1, enumerator_base::dim);
3680 lattice_point(V, C, t->terms[0], E_vertex);
3681
3682 // the elements of factors are always lexpositive
3683 mat_ZZ factors;
3684 s = setup_factors(C, factors, t, s);
3685
3686 t->factors[0] = new evalue;
3687 value_init(t->factors[0]->d);
3688 evalue_set_si(t->factors[0], s, 1);
3689 reduce(factors, v);
3690
3691 for (int i = 0; i < enumerator_base::dim; ++i)
3692 if (E_vertex[i]) {
3693 free_evalue_refs(E_vertex[i]);
3694 delete E_vertex[i];
3695 }
3696 }
3697
3698 #ifdef HAVE_CORRECT_VERTICES
3699 static inline Param_Polyhedron *Polyhedron2Param_SD(Polyhedron **Din,
3700 Polyhedron *Cin,int WS,Polyhedron **CEq,Matrix **CT)
3701 {
3702 if (WS & POL_NO_DUAL)
3703 WS = 0;
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 CA = align_context(C, P->Dimension, MaxRays);
3796 P = DomainIntersection(P, CA, MaxRays);
3797 Polyhedron_Free(CA);
3798
3799 /* for now */
3800 POL_ENSURE_FACETS(P);
3801 POL_ENSURE_VERTICES(P);
3802 POL_ENSURE_FACETS(C);
3803 POL_ENSURE_VERTICES(C);
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 || (P->Dimension == nparam+1)) {
3836 Polyhedron *Q;
3837 Polyhedron *C2;
3838 for (Q = T ? T : P; 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 }
3856 if (T) {
3857 Polyhedron_Free(P);
3858 P = T;
3859 if (T->Dimension == C->Dimension) {
3860 P = T->next;
3861 T->next = NULL;
3862 Polyhedron_Free(T);
3863 }
3864 }
3865
3866 Polyhedron *next = P->next;
3867 P->next = NULL;
3868 eres = barvinok_enumerate_ev_f(P, C, MaxRays);
3869 P->next = next;
3870
3871 if (P->next) {
3872 Polyhedron *Q;
3873 evalue *f;
3874
3875 for (Q = P->next; Q; Q = Q->next) {
3876 Polyhedron *next = Q->next;
3877 Q->next = NULL;
3878
3879 f = barvinok_enumerate_ev_f(Q, C, MaxRays);
3880 emul(f, eres);
3881 free_evalue_refs(f);
3882 free(f);
3883
3884 Q->next = next;
3885 }
3886 }
3887
3888 goto out;
3889 }
3890
3891 static evalue* barvinok_enumerate_ev_f(Polyhedron *P, Polyhedron* C,
3892 unsigned MaxRays)
3893 {
3894 unsigned nparam = C->Dimension;
3895
3896 if (P->Dimension - nparam == 1)
3897 return ParamLine_Length(P, C, MaxRays);
3898
3899 Param_Polyhedron *PP = NULL;
3900 Polyhedron *CEq = NULL, *pVD;
3901 Matrix *CT = NULL;
3902 Param_Domain *D, *next;
3903 Param_Vertices *V;
3904 evalue *eres;
3905 Polyhedron *Porig = P;
3906
3907 PP = Polyhedron2Param_SD(&P,C,MaxRays,&CEq,&CT);
3908
3909 if (isIdentity(CT)) {
3910 Matrix_Free(CT);
3911 CT = NULL;
3912 } else {
3913 assert(CT->NbRows != CT->NbColumns);
3914 if (CT->NbRows == 1) { // no more parameters
3915 eres = barvinok_enumerate_cst(P, CEq, MaxRays);
3916 out:
3917 if (CT)
3918 Matrix_Free(CT);
3919 if (PP)
3920 Param_Polyhedron_Free(PP);
3921 if (P != Porig)
3922 Polyhedron_Free(P);
3923
3924 return eres;
3925 }
3926 nparam = CT->NbRows - 1;
3927 }
3928
3929 unsigned dim = P->Dimension - nparam;
3930
3931 ALLOC(evalue, eres);
3932 value_init(eres->d);
3933 value_set_si(eres->d, 0);
3934
3935 int nd;
3936 for (nd = 0, D=PP->D; D; ++nd, D=D->next);
3937 struct section { Polyhedron *D; evalue E; };
3938 section *s = new section[nd];
3939 Polyhedron **fVD = new Polyhedron_p[nd];
3940
3941 try_again:
3942 #ifdef USE_INCREMENTAL_BF
3943 bfenumerator et(P, dim, PP->nbV);
3944 #elif defined USE_INCREMENTAL_DF
3945 ienumerator et(P, dim, PP->nbV);
3946 #else
3947 enumerator et(P, dim, PP->nbV);
3948 #endif
3949
3950 for(nd = 0, D=PP->D; D; D=next) {
3951 next = D->next;
3952
3953 Polyhedron *rVD = reduce_domain(D->Domain, CT, CEq,
3954 fVD, nd, MaxRays);
3955 if (!rVD)
3956 continue;
3957
3958 pVD = CT ? DomainImage(rVD,CT,MaxRays) : rVD;
3959
3960 value_init(s[nd].E.d);
3961 evalue_set_si(&s[nd].E, 0, 1);
3962 s[nd].D = rVD;
3963
3964 FORALL_PVertex_in_ParamPolyhedron(V,D,PP) // _i is internal counter
3965 if (!et.vE[_i])
3966 try {
3967 et.decompose_at(V, _i, MaxRays);
3968 } catch (OrthogonalException &e) {
3969 if (rVD != pVD)
3970 Domain_Free(pVD);
3971 for (; nd >= 0; --nd) {
3972 free_evalue_refs(&s[nd].E);
3973 Domain_Free(s[nd].D);
3974 Domain_Free(fVD[nd]);
3975 }
3976 goto try_again;
3977 }
3978 eadd(et.vE[_i] , &s[nd].E);
3979 END_FORALL_PVertex_in_ParamPolyhedron;
3980 reduce_in_domain(&s[nd].E, pVD);
3981
3982 if (CT)
3983 addeliminatedparams_evalue(&s[nd].E, CT);
3984 ++nd;
3985 if (rVD != pVD)
3986 Domain_Free(pVD);
3987 }
3988
3989 if (nd == 0)
3990 evalue_set_si(eres, 0, 1);
3991 else {
3992 eres->x.p = new_enode(partition, 2*nd, C->Dimension);
3993 for (int j = 0; j < nd; ++j) {
3994 EVALUE_SET_DOMAIN(eres->x.p->arr[2*j], s[j].D);
3995 value_clear(eres->x.p->arr[2*j+1].d);
3996 eres->x.p->arr[2*j+1] = s[j].E;
3997 Domain_Free(fVD[j]);
3998 }
3999 }
4000 delete [] s;
4001 delete [] fVD;
4002
4003 if (CEq)
4004 Polyhedron_Free(CEq);
4005 goto out;
4006 }
4007
4008 Enumeration* barvinok_enumerate(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
4009 {
4010 evalue *EP = barvinok_enumerate_ev(P, C, MaxRays);
4011
4012 return partition2enumeration(EP);
4013 }
4014
4015 static void SwapColumns(Value **V, int n, int i, int j)
4016 {
4017 for (int r = 0; r < n; ++r)
4018 value_swap(V[r][i], V[r][j]);
4019 }
4020
4021 static void SwapColumns(Polyhedron *P, int i, int j)
4022 {
4023 SwapColumns(P->Constraint, P->NbConstraints, i, j);
4024 SwapColumns(P->Ray, P->NbRays, i, j);
4025 }
4026
4027 static void negative_test_constraint(Value *l, Value *u, Value *c, int pos,
4028 int len, Value *v)
4029 {
4030 value_oppose(*v, u[pos+1]);
4031 Vector_Combine(l+1, u+1, c+1, *v, l[pos+1], len-1);
4032 value_multiply(*v, *v, l[pos+1]);
4033 value_subtract(c[len-1], c[len-1], *v);
4034 value_set_si(*v, -1);
4035 Vector_Scale(c+1, c+1, *v, len-1);
4036 value_decrement(c[len-1], c[len-1]);
4037 ConstraintSimplify(c, c, len, v);
4038 }
4039
4040 static bool parallel_constraints(Value *l, Value *u, Value *c, int pos,
4041 int len)
4042 {
4043 bool parallel;
4044 Value g1;
4045 Value g2;
4046 value_init(g1);
4047 value_init(g2);
4048
4049 Vector_Gcd(&l[1+pos], len, &g1);
4050 Vector_Gcd(&u[1+pos], len, &g2);
4051 Vector_Combine(l+1+pos, u+1+pos, c+1, g2, g1, len);
4052 parallel = First_Non_Zero(c+1, len) == -1;
4053
4054 value_clear(g1);
4055 value_clear(g2);
4056
4057 return parallel;
4058 }
4059
4060 static void negative_test_constraint7(Value *l, Value *u, Value *c, int pos,
4061 int exist, int len, Value *v)
4062 {
4063 Value g;
4064 value_init(g);
4065
4066 Vector_Gcd(&u[1+pos], exist, v);
4067 Vector_Gcd(&l[1+pos], exist, &g);
4068 Vector_Combine(l+1, u+1, c+1, *v, g, len-1);
4069 value_multiply(*v, *v, g);
4070 value_subtract(c[len-1], c[len-1], *v);
4071 value_set_si(*v, -1);
4072 Vector_Scale(c+1, c+1, *v, len-1);
4073 value_decrement(c[len-1], c[len-1]);
4074 ConstraintSimplify(c, c, len, v);
4075
4076 value_clear(g);
4077 }
4078
4079 static void oppose_constraint(Value *c, int len, Value *v)
4080 {
4081 value_set_si(*v, -1);
4082 Vector_Scale(c+1, c+1, *v, len-1);
4083 value_decrement(c[len-1], c[len-1]);
4084 }
4085
4086 static bool SplitOnConstraint(Polyhedron *P, int i, int l, int u,
4087 int nvar, int len, int exist, int MaxRays,
4088 Vector *row, Value& f, bool independent,
4089 Polyhedron **pos, Polyhedron **neg)
4090 {
4091 negative_test_constraint(P->Constraint[l], P->Constraint[u],
4092 row->p, nvar+i, len, &f);
4093 *neg = AddConstraints(row->p, 1, P, MaxRays);
4094
4095 /* We found an independent, but useless constraint
4096 * Maybe we should detect this earlier and not
4097 * mark the variable as INDEPENDENT
4098 */
4099 if (emptyQ((*neg))) {
4100 Polyhedron_Free(*neg);
4101 return false;
4102 }
4103
4104 oppose_constraint(row->p, len, &f);
4105 *pos = AddConstraints(row->p, 1, P, MaxRays);
4106
4107 if (emptyQ((*pos))) {
4108 Polyhedron_Free(*neg);
4109 Polyhedron_Free(*pos);
4110 return false;
4111 }
4112
4113 return true;
4114 }
4115
4116 /*
4117 * unimodularly transform P such that constraint r is transformed
4118 * into a constraint that involves only a single (the first)
4119 * existential variable
4120 *
4121 */
4122 static Polyhedron *rotate_along(Polyhedron *P, int r, int nvar, int exist,
4123 unsigned MaxRays)
4124 {
4125 Value g;
4126 value_init(g);
4127
4128 Vector *row = Vector_Alloc(exist);
4129 Vector_Copy(P->Constraint[r]+1+nvar, row->p, exist);
4130 Vector_Gcd(row->p, exist, &g);
4131 if (value_notone_p(g))
4132 Vector_AntiScale(row->p, row->p, g, exist);
4133 value_clear(g);
4134
4135 Matrix *M = unimodular_complete(row);
4136 Matrix *M2 = Matrix_Alloc(P->Dimension+1, P->Dimension+1);
4137 for (r = 0; r < nvar; ++r)
4138 value_set_si(M2->p[r][r], 1);
4139 for ( ; r < nvar+exist; ++r)
4140 Vector_Copy(M->p[r-nvar], M2->p[r]+nvar, exist);
4141 for ( ; r < P->Dimension+1; ++r)
4142 value_set_si(M2->p[r][r], 1);
4143 Polyhedron *T = Polyhedron_Image(P, M2, MaxRays);
4144
4145 Matrix_Free(M2);
4146 Matrix_Free(M);
4147 Vector_Free(row);
4148
4149 return T;
4150 }
4151
4152 static bool SplitOnVar(Polyhedron *P, int i,
4153 int nvar, int len, int exist, int MaxRays,
4154 Vector *row, Value& f, bool independent,
4155 Polyhedron **pos, Polyhedron **neg)
4156 {
4157 int j;
4158
4159 for (int l = P->NbEq; l < P->NbConstraints; ++l) {
4160 if (value_negz_p(P->Constraint[l][nvar+i+1]))
4161 continue;
4162
4163 if (independent) {
4164 for (j = 0; j < exist; ++j)
4165 if (j != i && value_notzero_p(P->Constraint[l][nvar+j+1]))
4166 break;
4167 if (j < exist)
4168 continue;
4169 }
4170
4171 for (int u = P->NbEq; u < P->NbConstraints; ++u) {
4172 if (value_posz_p(P->Constraint[u][nvar+i+1]))
4173 continue;
4174
4175 if (independent) {
4176 for (j = 0; j < exist; ++j)
4177 if (j != i && value_notzero_p(P->Constraint[u][nvar+j+1]))
4178 break;
4179 if (j < exist)
4180 continue;
4181 }
4182
4183 if (SplitOnConstraint(P, i, l, u,
4184 nvar, len, exist, MaxRays,
4185 row, f, independent,
4186 pos, neg)) {
4187 if (independent) {
4188 if (i != 0)
4189 SwapColumns(*neg, nvar+1, nvar+1+i);
4190 }
4191 return true;
4192 }
4193 }
4194 }
4195
4196 return false;
4197 }
4198
4199 static bool double_bound_pair(Polyhedron *P, int nvar, int exist,
4200 int i, int l1, int l2,
4201 Polyhedron **pos, Polyhedron **neg)
4202 {
4203 Value f;
4204 value_init(f);
4205 Vector *row = Vector_Alloc(P->Dimension+2);
4206 value_set_si(row->p[0], 1);
4207 value_oppose(f, P->Constraint[l1][nvar+i+1]);
4208 Vector_Combine(P->Constraint[l1]+1, P->Constraint[l2]+1,
4209 row->p+1,
4210 P->Constraint[l2][nvar+i+1], f,
4211 P->Dimension+1);
4212 ConstraintSimplify(row->p, row->p, P->Dimension+2, &f);
4213 *pos = AddConstraints(row->p, 1, P, 0);
4214 value_set_si(f, -1);
4215 Vector_Scale(row->p+1, row->p+1, f, P->Dimension+1);
4216 value_decrement(row->p[P->Dimension+1], row->p[P->Dimension+1]);
4217 *neg = AddConstraints(row->p, 1, P, 0);
4218 Vector_Free(row);
4219 value_clear(f);
4220
4221 return !emptyQ((*pos)) && !emptyQ((*neg));
4222 }
4223
4224 static bool double_bound(Polyhedron *P, int nvar, int exist,
4225 Polyhedron **pos, Polyhedron **neg)
4226 {
4227 for (int i = 0; i < exist; ++i) {
4228 int l1, l2;
4229 for (l1 = P->NbEq; l1 < P->NbConstraints; ++l1) {
4230 if (value_negz_p(P->Constraint[l1][nvar+i+1]))
4231 continue;
4232 for (l2 = l1 + 1; l2 < P->NbConstraints; ++l2) {
4233 if (value_negz_p(P->Constraint[l2][nvar+i+1]))
4234 continue;
4235 if (double_bound_pair(P, nvar, exist, i, l1, l2, pos, neg))
4236 return true;
4237 }
4238 }
4239 for (l1 = P->NbEq; l1 < P->NbConstraints; ++l1) {
4240 if (value_posz_p(P->Constraint[l1][nvar+i+1]))
4241 continue;
4242 if (l1 < P->NbConstraints)
4243 for (l2 = l1 + 1; l2 < P->NbConstraints; ++l2) {
4244 if (value_posz_p(P->Constraint[l2][nvar+i+1]))
4245 continue;
4246 if (double_bound_pair(P, nvar, exist, i, l1, l2, pos, neg))
4247 return true;
4248 }
4249 }
4250 return false;
4251 }
4252 return false;
4253 }
4254
4255 enum constraint {
4256 ALL_POS = 1 << 0,
4257 ONE_NEG = 1 << 1,
4258 INDEPENDENT = 1 << 2,
4259 ROT_NEG = 1 << 3
4260 };
4261
4262 static evalue* enumerate_or(Polyhedron *D,
4263 unsigned exist, unsigned nparam, unsigned MaxRays)
4264 {
4265 #ifdef DEBUG_ER
4266 fprintf(stderr, "\nER: Or\n");
4267 #endif /* DEBUG_ER */
4268
4269 Polyhedron *N = D->next;
4270 D->next = 0;
4271 evalue *EP =
4272 barvinok_enumerate_e(D, exist, nparam, MaxRays);
4273 Polyhedron_Free(D);
4274
4275 for (D = N; D; D = N) {
4276 N = D->next;
4277 D->next = 0;
4278
4279 evalue *EN =
4280 barvinok_enumerate_e(D, exist, nparam, MaxRays);
4281
4282 eor(EN, EP);
4283 free_evalue_refs(EN);
4284 free(EN);
4285 Polyhedron_Free(D);
4286 }
4287
4288 reduce_evalue(EP);
4289
4290 return EP;
4291 }
4292
4293 static evalue* enumerate_sum(Polyhedron *P,
4294 unsigned exist, unsigned nparam, unsigned MaxRays)
4295 {
4296 int nvar = P->Dimension - exist - nparam;
4297 int toswap = nvar < exist ? nvar : exist;
4298 for (int i = 0; i < toswap; ++i)
4299 SwapColumns(P, 1 + i, nvar+exist - i);
4300 nparam += nvar;
4301
4302 #ifdef DEBUG_ER
4303 fprintf(stderr, "\nER: Sum\n");
4304 #endif /* DEBUG_ER */
4305
4306 evalue *EP = barvinok_enumerate_e(P, exist, nparam, MaxRays);
4307
4308 for (int i = 0; i < /* nvar */ nparam; ++i) {
4309 Matrix *C = Matrix_Alloc(1, 1 + nparam + 1);
4310 value_set_si(C->p[0][0], 1);
4311 evalue split;
4312 value_init(split.d);
4313 value_set_si(split.d, 0);
4314 split.x.p = new_enode(partition, 4, nparam);
4315 value_set_si(C->p[0][1+i], 1);
4316 Matrix *C2 = Matrix_Copy(C);
4317 EVALUE_SET_DOMAIN(split.x.p->arr[0],
4318 Constraints2Polyhedron(C2, MaxRays));
4319 Matrix_Free(C2);
4320 evalue_set_si(&split.x.p->arr[1], 1, 1);
4321 value_set_si(C->p[0][1+i], -1);
4322 value_set_si(C->p[0][1+nparam], -1);
4323 EVALUE_SET_DOMAIN(split.x.p->arr[2],
4324 Constraints2Polyhedron(C, MaxRays));
4325 evalue_set_si(&split.x.p->arr[3], 1, 1);
4326 emul(&split, EP);
4327 free_evalue_refs(&split);
4328 Matrix_Free(C);
4329 }
4330 reduce_evalue(EP);
4331 evalue_range_reduction(EP);
4332
4333 evalue_frac2floor(EP);
4334
4335 evalue *sum = esum(EP, nvar);
4336
4337 free_evalue_refs(EP);
4338 free(EP);
4339 EP = sum;
4340
4341 evalue_range_reduction(EP);
4342
4343 return EP;
4344 }
4345
4346 static evalue* split_sure(Polyhedron *P, Polyhedron *S,
4347 unsigned exist, unsigned nparam, unsigned MaxRays)
4348 {
4349 int nvar = P->Dimension - exist - nparam;
4350
4351 Matrix *M = Matrix_Alloc(exist, S->Dimension+2);
4352 for (int i = 0; i < exist; ++i)
4353 value_set_si(M->p[i][nvar+i+1], 1);
4354 Polyhedron *O = S;
4355 S = DomainAddRays(S, M, MaxRays);
4356 Polyhedron_Free(O);
4357 Polyhedron *F = DomainAddRays(P, M, MaxRays);
4358 Polyhedron *D = DomainDifference(F, S, MaxRays);
4359 O = D;
4360 D = Disjoint_Domain(D, 0, MaxRays);
4361 Polyhedron_Free(F);
4362 Domain_Free(O);
4363 Matrix_Free(M);
4364
4365 M = Matrix_Alloc(P->Dimension+1-exist, P->Dimension+1);
4366 for (int j = 0; j < nvar; ++j)
4367 value_set_si(M->p[j][j], 1);
4368 for (int j = 0; j < nparam+1; ++j)
4369 value_set_si(M->p[nvar+j][nvar+exist+j], 1);
4370 Polyhedron *T = Polyhedron_Image(S, M, MaxRays);
4371 evalue *EP = barvinok_enumerate_e(T, 0, nparam, MaxRays);
4372 Polyhedron_Free(S);
4373 Polyhedron_Free(T);
4374 Matrix_Free(M);
4375
4376 for (Polyhedron *Q = D; Q; Q = Q->next) {
4377 Polyhedron *N = Q->next;
4378 Q->next = 0;
4379 T = DomainIntersection(P, Q, MaxRays);
4380 evalue *E = barvinok_enumerate_e(T, exist, nparam, MaxRays);
4381 eadd(E, EP);
4382 free_evalue_refs(E);
4383 free(E);
4384 Polyhedron_Free(T);
4385 Q->next = N;
4386 }
4387 Domain_Free(D);
4388 return EP;
4389 }
4390
4391 static evalue* enumerate_sure(Polyhedron *P,
4392 unsigned exist, unsigned nparam, unsigned MaxRays)
4393 {
4394 int i;
4395 Polyhedron *S = P;
4396 int nvar = P->Dimension - exist - nparam;
4397 Value lcm;
4398 Value f;
4399 value_init(lcm);
4400 value_init(f);
4401
4402 for (i = 0; i < exist; ++i) {
4403 Matrix *M = Matrix_Alloc(S->NbConstraints, S->Dimension+2);
4404 int c = 0;
4405 value_set_si(lcm, 1);
4406 for (int j = 0; j < S->NbConstraints; ++j) {
4407 if (value_negz_p(S->Constraint[j][1+nvar+i]))
4408 continue;
4409 if (value_one_p(S->Constraint[j][1+nvar+i]))
4410 continue;
4411 value_lcm(lcm, S->Constraint[j][1+nvar+i], &lcm);
4412 }
4413
4414 for (int j = 0; j < S->NbConstraints; ++j) {
4415 if (value_negz_p(S->Constraint[j][1+nvar+i]))
4416 continue;
4417 if (value_one_p(S->Constraint[j][1+nvar+i]))
4418 continue;
4419 value_division(f, lcm, S->Constraint[j][1+nvar+i]);
4420 Vector_Scale(S->Constraint[j], M->p[c], f, S->Dimension+2);
4421 value_subtract(M->p[c][S->Dimension+1],
4422 M->p[c][S->Dimension+1],
4423 lcm);
4424 value_increment(M->p[c][S->Dimension+1],
4425 M->p[c][S->Dimension+1]);
4426 ++c;
4427 }
4428 Polyhedron *O = S;
4429 S = AddConstraints(M->p[0], c, S, MaxRays);
4430 if (O != P)
4431 Polyhedron_Free(O);
4432 Matrix_Free(M);
4433 if (emptyQ(S)) {
4434 Polyhedron_Free(S);
4435 value_clear(lcm);
4436 value_clear(f);
4437 return 0;
4438 }
4439 }
4440 value_clear(lcm);
4441 value_clear(f);
4442
4443 #ifdef DEBUG_ER
4444 fprintf(stderr, "\nER: Sure\n");
4445 #endif /* DEBUG_ER */
4446
4447 return split_sure(P, S, exist, nparam, MaxRays);
4448 }
4449
4450 static evalue* enumerate_sure2(Polyhedron *P,
4451 unsigned exist, unsigned nparam, unsigned MaxRays)
4452 {
4453 int nvar = P->Dimension - exist - nparam;
4454 int r;
4455 for (r = 0; r < P->NbRays; ++r)
4456 if (value_one_p(P->Ray[r][0]) &&
4457 value_one_p(P->Ray[r][P->Dimension+1]))
4458 break;
4459
4460 if (r >= P->NbRays)
4461 return 0;
4462
4463 Matrix *M = Matrix_Alloc(nvar + 1 + nparam, P->Dimension+2);
4464 for (int i = 0; i < nvar; ++i)
4465 value_set_si(M->p[i][1+i], 1);
4466 for (int i = 0; i < nparam; ++i)
4467 value_set_si(M->p[i+nvar][1+nvar+exist+i], 1);
4468 Vector_Copy(P->Ray[r]+1+nvar, M->p[nvar+nparam]+1+nvar, exist);
4469 value_set_si(M->p[nvar+nparam][0], 1);
4470 value_set_si(M->p[nvar+nparam][P->Dimension+1], 1);
4471 Polyhedron * F = Rays2Polyhedron(M, MaxRays);
4472 Matrix_Free(M);
4473
4474 Polyhedron *I = DomainIntersection(F, P, MaxRays);
4475 Polyhedron_Free(F);
4476
4477 #ifdef DEBUG_ER
4478 fprintf(stderr, "\nER: Sure2\n");
4479 #endif /* DEBUG_ER */
4480
4481 return split_sure(P, I, exist, nparam, MaxRays);
4482 }
4483
4484 static evalue* enumerate_cyclic(Polyhedron *P,
4485 unsigned exist, unsigned nparam,
4486 evalue * EP, int r, int p, unsigned MaxRays)
4487 {
4488 int nvar = P->Dimension - exist - nparam;
4489
4490 /* If EP in its fractional maps only contains references
4491 * to the remainder parameter with appropriate coefficients
4492 * then we could in principle avoid adding existentially
4493 * quantified variables to the validity domains.
4494 * We'd have to replace the remainder by m { p/m }
4495 * and multiply with an appropriate factor that is one
4496 * only in the appropriate range.
4497 * This last multiplication can be avoided if EP
4498 * has a single validity domain with no (further)
4499 * constraints on the remainder parameter
4500 */
4501
4502 Matrix *CT = Matrix_Alloc(nparam+1, nparam+3);
4503 Matrix *M = Matrix_Alloc(1, 1+nparam+3);
4504 for (int j = 0; j < nparam; ++j)
4505 if (j != p)
4506 value_set_si(CT->p[j][j], 1);
4507 value_set_si(CT->p[p][nparam+1], 1);
4508 value_set_si(CT->p[nparam][nparam+2], 1);
4509 value_set_si(M->p[0][1+p], -1);
4510 value_absolute(M->p[0][1+nparam], P->Ray[0][1+nvar+exist+p]);
4511 value_set_si(M->p[0][1+nparam+1], 1);
4512 Polyhedron *CEq = Constraints2Polyhedron(M, 1);
4513 Matrix_Free(M);
4514 addeliminatedparams_enum(EP, CT, CEq, MaxRays, nparam);
4515 Polyhedron_Free(CEq);
4516 Matrix_Free(CT);
4517
4518 return EP;
4519 }
4520
4521 static void enumerate_vd_add_ray(evalue *EP, Matrix *Rays, unsigned MaxRays)
4522 {
4523 if (value_notzero_p(EP->d))
4524 return;
4525
4526 assert(EP->x.p->type == partition);
4527 assert(EP->x.p->pos == EVALUE_DOMAIN(EP->x.p->arr[0])->Dimension);
4528 for (int i = 0; i < EP->x.p->size/2; ++i) {
4529 Polyhedron *D = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
4530 Polyhedron *N = DomainAddRays(D, Rays, MaxRays);
4531 EVALUE_SET_DOMAIN(EP->x.p->arr[2*i], N);
4532 Domain_Free(D);
4533 }
4534 }
4535
4536 static evalue* enumerate_line(Polyhedron *P,
4537 unsigned exist, unsigned nparam, unsigned MaxRays)
4538 {
4539 if (P->NbBid == 0)
4540 return 0;
4541
4542 #ifdef DEBUG_ER
4543 fprintf(stderr, "\nER: Line\n");
4544 #endif /* DEBUG_ER */
4545
4546 int nvar = P->Dimension - exist - nparam;
4547 int i, j;
4548 for (i = 0; i < nparam; ++i)
4549 if (value_notzero_p(P->Ray[0][1+nvar+exist+i]))
4550 break;
4551 assert(i < nparam);
4552 for (j = i+1; j < nparam; ++j)
4553 if (value_notzero_p(P->Ray[0][1+nvar+exist+i]))
4554 break;
4555 assert(j >= nparam); // for now
4556
4557 Matrix *M = Matrix_Alloc(2, P->Dimension+2);
4558 value_set_si(M->p[0][0], 1);
4559 value_set_si(M->p[0][1+nvar+exist+i], 1);
4560 value_set_si(M->p[1][0], 1);
4561 value_set_si(M->p[1][1+nvar+exist+i], -1);
4562 value_absolute(M->p[1][1+P->Dimension], P->Ray[0][1+nvar+exist+i]);
4563 value_decrement(M->p[1][1+P->Dimension], M->p[1][1+P->Dimension]);
4564 Polyhedron *S = AddConstraints(M->p[0], 2, P, MaxRays);
4565 evalue *EP = barvinok_enumerate_e(S, exist, nparam, MaxRays);
4566 Polyhedron_Free(S);
4567 Matrix_Free(M);
4568
4569 return enumerate_cyclic(P, exist, nparam, EP, 0, i, MaxRays);
4570 }
4571
4572 static int single_param_pos(Polyhedron*P, unsigned exist, unsigned nparam,
4573 int r)
4574 {
4575 int nvar = P->Dimension - exist - nparam;
4576 if (First_Non_Zero(P->Ray[r]+1, nvar) != -1)
4577 return -1;
4578 int i = First_Non_Zero(P->Ray[r]+1+nvar+exist, nparam);
4579 if (i == -1)
4580 return -1;
4581 if (First_Non_Zero(P->Ray[r]+1+nvar+exist+1, nparam-i-1) != -1)
4582 return -1;
4583 return i;
4584 }
4585
4586 static evalue* enumerate_remove_ray(Polyhedron *P, int r,
4587 unsigned exist, unsigned nparam, unsigned MaxRays)
4588 {
4589 #ifdef DEBUG_ER
4590 fprintf(stderr, "\nER: RedundantRay\n");
4591 #endif /* DEBUG_ER */
4592
4593 Value one;
4594 value_init(one);
4595 value_set_si(one, 1);
4596 int len = P->NbRays-1;
4597 Matrix *M = Matrix_Alloc(2 * len, P->Dimension+2);
4598 Vector_Copy(P->Ray[0], M->p[0], r * (P->Dimension+2));
4599 Vector_Copy(P->Ray[r+1], M->p[r], (len-r) * (P->Dimension+2));
4600 for (int j = 0; j < P->NbRays; ++j) {
4601 if (j == r)
4602 continue;
4603 Vector_Combine(P->Ray[j], P->Ray[r], M->p[len+j-(j>r)],
4604 one, P->Ray[j][P->Dimension+1], P->Dimension+2);
4605 }
4606
4607 P = Rays2Polyhedron(M, MaxRays);
4608 Matrix_Free(M);
4609 evalue *EP = barvinok_enumerate_e(P, exist, nparam, MaxRays);
4610 Polyhedron_Free(P);
4611 value_clear(one);
4612
4613 return EP;
4614 }
4615
4616 static evalue* enumerate_redundant_ray(Polyhedron *P,
4617 unsigned exist, unsigned nparam, unsigned MaxRays)
4618 {
4619 assert(P->NbBid == 0);
4620 int nvar = P->Dimension - exist - nparam;
4621 Value m;
4622 value_init(m);
4623
4624 for (int r = 0; r < P->NbRays; ++r) {
4625 if (value_notzero_p(P->Ray[r][P->Dimension+1]))
4626 continue;
4627 int i1 = single_param_pos(P, exist, nparam, r);
4628 if (i1 == -1)
4629 continue;
4630 for (int r2 = r+1; r2 < P->NbRays; ++r2) {
4631 if (value_notzero_p(P->Ray[r2][P->Dimension+1]))
4632 continue;
4633 int i2 = single_param_pos(P, exist, nparam, r2);
4634 if (i2 == -1)
4635 continue;
4636 if (i1 != i2)
4637 continue;
4638
4639 value_division(m, P->Ray[r][1+nvar+exist+i1],
4640 P->Ray[r2][1+nvar+exist+i1]);
4641 value_multiply(m, m, P->Ray[r2][1+nvar+exist+i1]);
4642 /* r2 divides r => r redundant */
4643 if (value_eq(m, P->Ray[r][1+nvar+exist+i1])) {
4644 value_clear(m);
4645 return enumerate_remove_ray(P, r, exist, nparam, MaxRays);
4646 }
4647
4648 value_division(m, P->Ray[r2][1+nvar+exist+i1],
4649 P->Ray[r][1+nvar+exist+i1]);
4650 value_multiply(m, m, P->Ray[r][1+nvar+exist+i1]);
4651 /* r divides r2 => r2 redundant */
4652 if (value_eq(m, P->Ray[r2][1+nvar+exist+i1])) {
4653 value_clear(m);
4654 return enumerate_remove_ray(P, r2, exist, nparam, MaxRays);
4655 }
4656 }
4657 }
4658 value_clear(m);
4659 return 0;
4660 }
4661
4662 static Polyhedron *upper_bound(Polyhedron *P,
4663 int pos, Value *max, Polyhedron **R)
4664 {
4665 Value v;
4666 int r;
4667 value_init(v);
4668
4669 *R = 0;
4670 Polyhedron *N;
4671 Polyhedron *B = 0;
4672 for (Polyhedron *Q = P; Q; Q = N) {
4673 N = Q->next;
4674 for (r = 0; r < P->NbRays; ++r) {
4675 if (value_zero_p(P->Ray[r][P->Dimension+1]) &&
4676 value_pos_p(P->Ray[r][1+pos]))
4677 break;
4678 }
4679 if (r < P->NbRays) {
4680 Q->next = *R;
4681 *R = Q;
4682 continue;
4683 } else {
4684 Q->next = B;
4685 B = Q;
4686 }
4687 for (r = 0; r < P->NbRays; ++r) {
4688 if (value_zero_p(P->Ray[r][P->Dimension+1]))
4689 continue;
4690 mpz_fdiv_q(v, P->Ray[r][1+pos], P->Ray[r][1+P->Dimension]);
4691 if ((!Q->next && r == 0) || value_gt(v, *max))
4692 value_assign(*max, v);
4693 }
4694 }
4695 value_clear(v);
4696 return B;
4697 }
4698
4699 static evalue* enumerate_ray(Polyhedron *P,
4700 unsigned exist, unsigned nparam, unsigned MaxRays)
4701 {
4702 assert(P->NbBid == 0);
4703 int nvar = P->Dimension - exist - nparam;
4704
4705 int r;
4706 for (r = 0; r < P->NbRays; ++r)
4707 if (value_zero_p(P->Ray[r][P->Dimension+1]))
4708 break;
4709 if (r >= P->NbRays)
4710 return 0;
4711
4712 int r2;
4713 for (r2 = r+1; r2 < P->NbRays; ++r2)
4714 if (value_zero_p(P->Ray[r2][P->Dimension+1]))
4715 break;
4716 if (r2 < P->NbRays) {
4717 if (nvar > 0)
4718 return enumerate_sum(P, exist, nparam, MaxRays);
4719 }
4720
4721 #ifdef DEBUG_ER
4722 fprintf(stderr, "\nER: Ray\n");
4723 #endif /* DEBUG_ER */
4724
4725 Value m;
4726 Value one;
4727 value_init(m);
4728 value_init(one);
4729 value_set_si(one, 1);
4730 int i = single_param_pos(P, exist, nparam, r);
4731 assert(i != -1); // for now;
4732
4733 Matrix *M = Matrix_Alloc(P->NbRays, P->Dimension+2);
4734 for (int j = 0; j < P->NbRays; ++j) {
4735 Vector_Combine(P->Ray[j], P->Ray[r], M->p[j],
4736 one, P->Ray[j][P->Dimension+1], P->Dimension+2);
4737 }
4738 Polyhedron *S = Rays2Polyhedron(M, MaxRays);
4739 Matrix_Free(M);
4740 Polyhedron *D = DomainDifference(P, S, MaxRays);
4741 Polyhedron_Free(S);
4742 // Polyhedron_Print(stderr, P_VALUE_FMT, D);
4743 assert(value_pos_p(P->Ray[r][1+nvar+exist+i])); // for now
4744 Polyhedron *R;
4745 D = upper_bound(D, nvar+exist+i, &m, &R);
4746 assert(D);
4747 Domain_Free(D);
4748
4749 M = Matrix_Alloc(2, P->Dimension+2);
4750 value_set_si(M->p[0][0], 1);
4751 value_set_si(M->p[1][0], 1);
4752 value_set_si(M->p[0][1+nvar+exist+i], -1);
4753 value_set_si(M->p[1][1+nvar+exist+i], 1);
4754 value_assign(M->p[0][1+P->Dimension], m);
4755 value_oppose(M->p[1][1+P->Dimension], m);
4756 value_addto(M->p[1][1+P->Dimension], M->p[1][1+P->Dimension],
4757 P->Ray[r][1+nvar+exist+i]);
4758 value_decrement(M->p[1][1+P->Dimension], M->p[1][1+P->Dimension]);
4759 // Matrix_Print(stderr, P_VALUE_FMT, M);
4760 D = AddConstraints(M->p[0], 2, P, MaxRays);
4761 // Polyhedron_Print(stderr, P_VALUE_FMT, D);
4762 value_subtract(M->p[0][1+P->Dimension], M->p[0][1+P->Dimension],
4763 P->Ray[r][1+nvar+exist+i]);
4764 // Matrix_Print(stderr, P_VALUE_FMT, M);
4765 S = AddConstraints(M->p[0], 1, P, MaxRays);
4766 // Polyhedron_Print(stderr, P_VALUE_FMT, S);
4767 Matrix_Free(M);
4768
4769 evalue *EP = barvinok_enumerate_e(D, exist, nparam, MaxRays);
4770 Polyhedron_Free(D);
4771 value_clear(one);
4772 value_clear(m);
4773
4774 if (value_notone_p(P->Ray[r][1+nvar+exist+i]))
4775 EP = enumerate_cyclic(P, exist, nparam, EP, r, i, MaxRays);
4776 else {
4777 M = Matrix_Alloc(1, nparam+2);
4778 value_set_si(M->p[0][0], 1);
4779 value_set_si(M->p[0][1+i], 1);
4780 enumerate_vd_add_ray(EP, M, MaxRays);
4781 Matrix_Free(M);
4782 }
4783
4784 if (!emptyQ(S)) {
4785 evalue *E = barvinok_enumerate_e(S, exist, nparam, MaxRays);
4786 eadd(E, EP);
4787 free_evalue_refs(E);
4788 free(E);
4789 }
4790 Polyhedron_Free(S);
4791
4792 if (R) {
4793 assert(nvar == 0);
4794 evalue *ER = enumerate_or(R, exist, nparam, MaxRays);
4795 eor(ER, EP);
4796 free_evalue_refs(ER);
4797 free(ER);
4798 }
4799
4800 return EP;
4801 }
4802
4803 static evalue* enumerate_vd(Polyhedron **PA,
4804 unsigned exist, unsigned nparam, unsigned MaxRays)
4805 {
4806 Polyhedron *P = *PA;
4807 int nvar = P->Dimension - exist - nparam;
4808 Param_Polyhedron *PP = NULL;
4809 Polyhedron *C = Universe_Polyhedron(nparam);
4810 Polyhedron *CEq;
4811 Matrix *CT;
4812 Polyhedron *PR = P;
4813 PP = Polyhedron2Param_SimplifiedDomain(&PR,C,MaxRays,&CEq,&CT);
4814 Polyhedron_Free(C);
4815
4816 int nd;
4817 Param_Domain *D, *last;
4818 Value c;
4819 value_init(c);
4820 for (nd = 0, D=PP->D; D; D=D->next, ++nd)
4821 ;
4822
4823 Polyhedron **VD = new Polyhedron_p[nd];
4824 Polyhedron **fVD = new Polyhedron_p[nd];
4825 for(nd = 0, D=PP->D; D; D=D->next) {
4826 Polyhedron *rVD = reduce_domain(D->Domain, CT, CEq,
4827 fVD, nd, MaxRays);
4828 if (!rVD)
4829 continue;
4830
4831 VD[nd++] = rVD;
4832 last = D;
4833 }
4834
4835 evalue *EP = 0;
4836
4837 if (nd == 0)
4838 EP = new_zero_ep();
4839
4840 /* This doesn't seem to have any effect */
4841 if (nd == 1) {
4842 Polyhedron *CA = align_context(VD[0], P->Dimension, MaxRays);
4843 Polyhedron *O = P;
4844 P = DomainIntersection(P, CA, MaxRays);
4845 if (O != *PA)
4846 Polyhedron_Free(O);
4847 Polyhedron_Free(CA);
4848 if (emptyQ(P))
4849 EP = new_zero_ep();
4850 }
4851
4852 if (!EP && CT->NbColumns != CT->NbRows) {
4853 Polyhedron *CEqr = DomainImage(CEq, CT, MaxRays);
4854 Polyhedron *CA = align_context(CEqr, PR->Dimension, MaxRays);
4855 Polyhedron *I = DomainIntersection(PR, CA, MaxRays);
4856 Polyhedron_Free(CEqr);
4857 Polyhedron_Free(CA);
4858 #ifdef DEBUG_ER
4859 fprintf(stderr, "\nER: Eliminate\n");
4860 #endif /* DEBUG_ER */
4861 nparam -= CT->NbColumns - CT->NbRows;
4862 EP = barvinok_enumerate_e(I, exist, nparam, MaxRays);
4863 nparam += CT->NbColumns - CT->NbRows;
4864 addeliminatedparams_enum(EP, CT, CEq, MaxRays, nparam);
4865 Polyhedron_Free(I);
4866 }
4867 if (PR != *PA)
4868 Polyhedron_Free(PR);
4869 PR = 0;
4870
4871 if (!EP && nd > 1) {
4872 #ifdef DEBUG_ER
4873 fprintf(stderr, "\nER: VD\n");
4874 #endif /* DEBUG_ER */
4875 for (int i = 0; i < nd; ++i) {
4876 Polyhedron *CA = align_context(VD[i], P->Dimension, MaxRays);
4877 Polyhedron *I = DomainIntersection(P, CA, MaxRays);
4878
4879 if (i == 0)
4880 EP = barvinok_enumerate_e(I, exist, nparam, MaxRays);
4881 else {
4882 evalue *E = barvinok_enumerate_e(I, exist, nparam, MaxRays);
4883 eadd(E, EP);
4884 free_evalue_refs(E);
4885 free(E);
4886 }
4887 Polyhedron_Free(I);
4888 Polyhedron_Free(CA);
4889 }
4890 }
4891
4892 for (int i = 0; i < nd; ++i) {
4893 Polyhedron_Free(VD[i]);
4894 Polyhedron_Free(fVD[i]);
4895 }
4896 delete [] VD;
4897 delete [] fVD;
4898 value_clear(c);
4899
4900 if (!EP && nvar == 0) {
4901 Value f;
4902 value_init(f);
4903 Param_Vertices *V, *V2;
4904 Matrix* M = Matrix_Alloc(1, P->Dimension+2);
4905
4906 FORALL_PVertex_in_ParamPolyhedron(V, last, PP) {
4907 bool found = false;
4908 FORALL_PVertex_in_ParamPolyhedron(V2, last, PP) {
4909 if (V == V2) {
4910 found = true;
4911 continue;
4912 }
4913 if (!found)
4914 continue;
4915 for (int i = 0; i < exist; ++i) {
4916 value_oppose(f, V->Vertex->p[i][nparam+1]);
4917 Vector_Combine(V->Vertex->p[i],
4918 V2->Vertex->p[i],
4919 M->p[0] + 1 + nvar + exist,
4920 V2->Vertex->p[i][nparam+1],
4921 f,
4922 nparam+1);
4923 int j;
4924 for (j = 0; j < nparam; ++j)
4925 if (value_notzero_p(M->p[0][1+nvar+exist+j]))
4926 break;
4927 if (j >= nparam)
4928 continue;
4929 ConstraintSimplify(M->p[0], M->p[0],
4930 P->Dimension+2, &f);
4931 value_set_si(M->p[0][0], 0);
4932 Polyhedron *para = AddConstraints(M->p[0], 1, P,
4933 MaxRays);
4934 if (emptyQ(para)) {
4935 Polyhedron_Free(para);
4936 continue;
4937 }
4938 Polyhedron *pos, *neg;
4939 value_set_si(M->p[0][0], 1);
4940 value_decrement(M->p[0][P->Dimension+1],
4941 M->p[0][P->Dimension+1]);
4942 neg = AddConstraints(M->p[0], 1, P, MaxRays);
4943 value_set_si(f, -1);
4944 Vector_Scale(M->p[0]+1, M->p[0]+1, f,
4945 P->Dimension+1);
4946 value_decrement(M->p[0][P->Dimension+1],
4947 M->p[0][P->Dimension+1]);
4948 value_decrement(M->p[0][P->Dimension+1],
4949 M->p[0][P->Dimension+1]);
4950 pos = AddConstraints(M->p[0], 1, P, MaxRays);
4951 if (emptyQ(neg) && emptyQ(pos)) {
4952 Polyhedron_Free(para);
4953 Polyhedron_Free(pos);
4954 Polyhedron_Free(neg);
4955 continue;
4956 }
4957 #ifdef DEBUG_ER
4958 fprintf(stderr, "\nER: Order\n");
4959 #endif /* DEBUG_ER */
4960 EP = barvinok_enumerate_e(para, exist, nparam, MaxRays);
4961 evalue *E;
4962 if (!emptyQ(pos)) {
4963 E = barvinok_enumerate_e(pos, exist, nparam, MaxRays);
4964 eadd(E, EP);
4965 free_evalue_refs(E);
4966 free(E);
4967 }
4968 if (!emptyQ(neg)) {
4969 E = barvinok_enumerate_e(neg, exist, nparam, MaxRays);
4970 eadd(E, EP);
4971 free_evalue_refs(E);
4972 free(E);
4973 }
4974 Polyhedron_Free(para);
4975 Polyhedron_Free(pos);
4976 Polyhedron_Free(neg);
4977 break;
4978 }
4979 if (EP)
4980 break;
4981 } END_FORALL_PVertex_in_ParamPolyhedron;
4982 if (EP)
4983 break;
4984 } END_FORALL_PVertex_in_ParamPolyhedron;
4985
4986 if (!EP) {
4987 /* Search for vertex coordinate to split on */
4988 /* First look for one independent of the parameters */
4989 FORALL_PVertex_in_ParamPolyhedron(V, last, PP) {
4990 for (int i = 0; i < exist; ++i) {
4991 int j;
4992 for (j = 0; j < nparam; ++j)
4993 if (value_notzero_p(V->Vertex->p[i][j]))
4994 break;
4995 if (j < nparam)
4996 continue;
4997 value_set_si(M->p[0][0], 1);
4998 Vector_Set(M->p[0]+1, 0, nvar+exist);
4999 Vector_Copy(V->Vertex->p[i],
5000 M->p[0] + 1 + nvar + exist, nparam+1);
5001 value_oppose(M->p[0][1+nvar+i],
5002 V->Vertex->p[i][nparam+1]);
5003
5004 Polyhedron *pos, *neg;
5005 value_set_si(M->p[0][0], 1);
5006 value_decrement(M->p[0][P->Dimension+1],
5007 M->p[0][P->Dimension+1]);
5008 neg = AddConstraints(M->p[0], 1, P, MaxRays);
5009 value_set_si(f, -1);
5010 Vector_Scale(M->p[0]+1, M->p[0]+1, f,
5011 P->Dimension+1);
5012 value_decrement(M->p[0][P->Dimension+1],
5013 M->p[0][P->Dimension+1]);
5014 value_decrement(M->p[0][P->Dimension+1],
5015 M->p[0][P->Dimension+1]);
5016 pos = AddConstraints(M->p[0], 1, P, MaxRays);
5017 if (emptyQ(neg) || emptyQ(pos)) {
5018 Polyhedron_Free(pos);
5019 Polyhedron_Free(neg);
5020 continue;
5021 }
5022 Polyhedron_Free(pos);
5023 value_increment(M->p[0][P->Dimension+1],
5024 M->p[0][P->Dimension+1]);
5025 pos = AddConstraints(M->p[0], 1, P, MaxRays);
5026 #ifdef DEBUG_ER
5027 fprintf(stderr, "\nER: Vertex\n");
5028 #endif /* DEBUG_ER */
5029 pos->next = neg;
5030 EP = enumerate_or(pos, exist, nparam, MaxRays);
5031 break;
5032 }
5033 if (EP)
5034 break;
5035 } END_FORALL_PVertex_in_ParamPolyhedron;
5036 }
5037
5038 if (!EP) {
5039 /* Search for vertex coordinate to split on */
5040 /* Now look for one that depends on the parameters */
5041 FORALL_PVertex_in_ParamPolyhedron(V, last, PP) {
5042 for (int i = 0; i < exist; ++i) {
5043 value_set_si(M->p[0][0], 1);
5044 Vector_Set(M->p[0]+1, 0, nvar+exist);
5045 Vector_Copy(V->Vertex->p[i],
5046 M->p[0] + 1 + nvar + exist, nparam+1);
5047 value_oppose(M->p[0][1+nvar+i],
5048 V->Vertex->p[i][nparam+1]);
5049
5050 Polyhedron *pos, *neg;
5051 value_set_si(M->p[0][0], 1);
5052 value_decrement(M->p[0][P->Dimension+1],
5053 M->p[0][P->Dimension+1]);
5054 neg = AddConstraints(M->p[0], 1, P, MaxRays);
5055 value_set_si(f, -1);
5056 Vector_Scale(M->p[0]+1, M->p[0]+1, f,
5057 P->Dimension+1);
5058 value_decrement(M->p[0][P->Dimension+1],
5059 M->p[0][P->Dimension+1]);
5060 value_decrement(M->p[0][P->Dimension+1],
5061 M->p[0][P->Dimension+1]);
5062 pos = AddConstraints(M->p[0], 1, P, MaxRays);
5063 if (emptyQ(neg) || emptyQ(pos)) {
5064 Polyhedron_Free(pos);
5065 Polyhedron_Free(neg);
5066 continue;
5067 }
5068 Polyhedron_Free(pos);
5069 value_increment(M->p[0][P->Dimension+1],
5070 M->p[0][P->Dimension+1]);
5071 pos = AddConstraints(M->p[0], 1, P, MaxRays);
5072 #ifdef DEBUG_ER
5073 fprintf(stderr, "\nER: ParamVertex\n");
5074 #endif /* DEBUG_ER */
5075 pos->next = neg;
5076 EP = enumerate_or(pos, exist, nparam, MaxRays);
5077 break;
5078 }
5079 if (EP)
5080 break;
5081 } END_FORALL_PVertex_in_ParamPolyhedron;
5082 }
5083
5084 Matrix_Free(M);
5085 value_clear(f);
5086 }
5087
5088 if (CEq)
5089 Polyhedron_Free(CEq);
5090 if (CT)
5091 Matrix_Free(CT);
5092 if (PP)
5093 Param_Polyhedron_Free(PP);
5094 *PA = P;
5095
5096 return EP;
5097 }
5098
5099 #ifndef HAVE_PIPLIB
5100 evalue *barvinok_enumerate_pip(Polyhedron *P,
5101 unsigned exist, unsigned nparam, unsigned MaxRays)
5102 {
5103 return 0;
5104 }
5105 #else
5106 evalue *barvinok_enumerate_pip(Polyhedron *P,
5107 unsigned exist, unsigned nparam, unsigned MaxRays)
5108 {
5109 int nvar = P->Dimension - exist - nparam;
5110 evalue *EP = new_zero_ep();
5111 Polyhedron *Q, *N;
5112
5113 #ifdef DEBUG_ER
5114 fprintf(stderr, "\nER: PIP\n");
5115 #endif /* DEBUG_ER */
5116
5117 Polyhedron *D = pip_projectout(P, nvar, exist, nparam);
5118 for (Q = D; Q; Q = N) {
5119 N = Q->next;
5120 Q->next = 0;
5121 evalue *E;
5122 exist = Q->Dimension - nvar - nparam;
5123 E = barvinok_enumerate_e(Q, exist, nparam, MaxRays);
5124 Polyhedron_Free(Q);
5125 eadd(E, EP);
5126 free_evalue_refs(E);
5127 free(E);
5128 }
5129
5130 return EP;
5131 }
5132 #endif
5133
5134
5135 static bool is_single(Value *row, int pos, int len)
5136 {
5137 return First_Non_Zero(row, pos) == -1 &&
5138 First_Non_Zero(row+pos+1, len-pos-1) == -1;
5139 }
5140
5141 static evalue* barvinok_enumerate_e_r(Polyhedron *P,
5142 unsigned exist, unsigned nparam, unsigned MaxRays);
5143
5144 #ifdef DEBUG_ER
5145 static int er_level = 0;
5146
5147 evalue* barvinok_enumerate_e(Polyhedron *P,
5148 unsigned exist, unsigned nparam, unsigned MaxRays)
5149 {
5150 fprintf(stderr, "\nER: level %i\n", er_level);
5151
5152 Polyhedron_PrintConstraints(stderr, P_VALUE_FMT, P);
5153 fprintf(stderr, "\nE %d\nP %d\n", exist, nparam);
5154 ++er_level;
5155 P = DomainConstraintSimplify(Polyhedron_Copy(P), MaxRays);
5156 evalue *EP = barvinok_enumerate_e_r(P, exist, nparam, MaxRays);
5157 Polyhedron_Free(P);
5158 --er_level;
5159 return EP;
5160 }
5161 #else
5162 evalue* barvinok_enumerate_e(Polyhedron *P,
5163 unsigned exist, unsigned nparam, unsigned MaxRays)
5164 {
5165 P = DomainConstraintSimplify(Polyhedron_Copy(P), MaxRays);
5166 evalue *EP = barvinok_enumerate_e_r(P, exist, nparam, MaxRays);
5167 Polyhedron_Free(P);
5168 return EP;
5169 }
5170 #endif
5171
5172 static evalue* barvinok_enumerate_e_r(Polyhedron *P,
5173 unsigned exist, unsigned nparam, unsigned MaxRays)
5174 {
5175 if (exist == 0) {
5176 Polyhedron *U = Universe_Polyhedron(nparam);
5177 evalue *EP = barvinok_enumerate_ev(P, U, MaxRays);
5178 //char *param_name[] = {"P", "Q", "R", "S", "T" };
5179 //print_evalue(stdout, EP, param_name);
5180 Polyhedron_Free(U);
5181 return EP;
5182 }
5183
5184 int nvar = P->Dimension - exist - nparam;
5185 int len = P->Dimension + 2;
5186
5187 /* for now */
5188 POL_ENSURE_FACETS(P);
5189 POL_ENSURE_VERTICES(P);
5190
5191 if (emptyQ(P))
5192 return new_zero_ep();
5193
5194 if (nvar == 0 && nparam == 0) {
5195 evalue *EP = new_zero_ep();
5196 barvinok_count(P, &EP->x.n, MaxRays);
5197 if (value_pos_p(EP->x.n))
5198 value_set_si(EP->x.n, 1);
5199 return EP;
5200 }
5201
5202 int r;
5203 for (r = 0; r < P->NbRays; ++r)
5204 if (value_zero_p(P->Ray[r][0]) ||
5205 value_zero_p(P->Ray[r][P->Dimension+1])) {
5206 int i;
5207 for (i = 0; i < nvar; ++i)
5208 if (value_notzero_p(P->Ray[r][i+1]))
5209 break;
5210 if (i >= nvar)
5211 continue;
5212 for (i = nvar + exist; i < nvar + exist + nparam; ++i)
5213 if (value_notzero_p(P->Ray[r][i+1]))
5214 break;
5215 if (i >= nvar + exist + nparam)
5216 break;
5217 }
5218 if (r < P->NbRays) {
5219 evalue *EP = new_zero_ep();
5220 value_set_si(EP->x.n, -1);
5221 return EP;
5222 }
5223
5224 int first;
5225 for (r = 0; r < P->NbEq; ++r)
5226 if ((first = First_Non_Zero(P->Constraint[r]+1+nvar, exist)) != -1)
5227 break;
5228 if (r < P->NbEq) {
5229 if (First_Non_Zero(P->Constraint[r]+1+nvar+first+1,
5230 exist-first-1) != -1) {
5231 Polyhedron *T = rotate_along(P, r, nvar, exist, MaxRays);
5232 #ifdef DEBUG_ER
5233 fprintf(stderr, "\nER: Equality\n");
5234 #endif /* DEBUG_ER */
5235 evalue *EP = barvinok_enumerate_e(T, exist-1, nparam, MaxRays);
5236 Polyhedron_Free(T);
5237 return EP;
5238 } else {
5239 #ifdef DEBUG_ER
5240 fprintf(stderr, "\nER: Fixed\n");
5241 #endif /* DEBUG_ER */
5242 if (first == 0)
5243 return barvinok_enumerate_e(P, exist-1, nparam, MaxRays);
5244 else {
5245 Polyhedron *T = Polyhedron_Copy(P);
5246 SwapColumns(T, nvar+1, nvar+1+first);
5247 evalue *EP = barvinok_enumerate_e(T, exist-1, nparam, MaxRays);
5248 Polyhedron_Free(T);
5249 return EP;
5250 }
5251 }
5252 }
5253
5254 Vector *row = Vector_Alloc(len);
5255 value_set_si(row->p[0], 1);
5256
5257 Value f;
5258 value_init(f);
5259
5260 enum constraint* info = new constraint[exist];
5261 for (int i = 0; i < exist; ++i) {
5262 info[i] = ALL_POS;
5263 for (int l = P->NbEq; l < P->NbConstraints; ++l) {
5264 if (value_negz_p(P->Constraint[l][nvar+i+1]))
5265 continue;
5266 bool l_parallel = is_single(P->Constraint[l]+nvar+1, i, exist);
5267 for (int u = P->NbEq; u < P->NbConstraints; ++u) {
5268 if (value_posz_p(P->Constraint[u][nvar+i+1]))
5269 continue;
5270 bool lu_parallel = l_parallel ||
5271 is_single(P->Constraint[u]+nvar+1, i, exist);
5272 value_oppose(f, P->Constraint[u][nvar+i+1]);
5273 Vector_Combine(P->Constraint[l]+1, P->Constraint[u]+1, row->p+1,
5274 f, P->Constraint[l][nvar+i+1], len-1);
5275 if (!(info[i] & INDEPENDENT)) {
5276 int j;
5277 for (j = 0; j < exist; ++j)
5278 if (j != i && value_notzero_p(row->p[nvar+j+1]))
5279 break;
5280 if (j == exist) {
5281 //printf("independent: i: %d, l: %d, u: %d\n", i, l, u);
5282 info[i] = (constraint)(info[i] | INDEPENDENT);
5283 }
5284 }
5285 if (info[i] & ALL_POS) {
5286 value_addto(row->p[len-1], row->p[len-1],
5287 P->Constraint[l][nvar+i+1]);
5288 value_addto(row->p[len-1], row->p[len-1], f);
5289 value_multiply(f, f, P->Constraint[l][nvar+i+1]);
5290 value_subtract(row->p[len-1], row->p[len-1], f);
5291 value_decrement(row->p[len-1], row->p[len-1]);
5292 ConstraintSimplify(row->p, row->p, len, &f);
5293 value_set_si(f, -1);
5294 Vector_Scale(row->p+1, row->p+1, f, len-1);
5295 value_decrement(row->p[len-1], row->p[len-1]);
5296 Polyhedron *T = AddConstraints(row->p, 1, P, MaxRays);
5297 if (!emptyQ(T)) {
5298 //printf("not all_pos: i: %d, l: %d, u: %d\n", i, l, u);
5299 info[i] = (constraint)(info[i] ^ ALL_POS);
5300 }
5301 //puts("pos remainder");
5302 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
5303 Polyhedron_Free(T);
5304 }
5305 if (!(info[i] & ONE_NEG)) {
5306 if (lu_parallel) {
5307 negative_test_constraint(P->Constraint[l],
5308 P->Constraint[u],
5309 row->p, nvar+i, len, &f);
5310 oppose_constraint(row->p, len, &f);
5311 Polyhedron *T = AddConstraints(row->p, 1, P, MaxRays);
5312 if (emptyQ(T)) {
5313 //printf("one_neg i: %d, l: %d, u: %d\n", i, l, u);
5314 info[i] = (constraint)(info[i] | ONE_NEG);
5315 }
5316 //puts("neg remainder");
5317 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
5318 Polyhedron_Free(T);
5319 } else if (!(info[i] & ROT_NEG)) {
5320 if (parallel_constraints(P->Constraint[l],
5321 P->Constraint[u],
5322 row->p, nvar, exist)) {
5323 negative_test_constraint7(P->Constraint[l],
5324 P->Constraint[u],
5325 row->p, nvar, exist,
5326 len, &f);
5327 oppose_constraint(row->p, len, &f);
5328 Polyhedron *T = AddConstraints(row->p, 1, P, MaxRays);
5329 if (emptyQ(T)) {
5330 // printf("rot_neg i: %d, l: %d, u: %d\n", i, l, u);
5331 info[i] = (constraint)(info[i] | ROT_NEG);
5332 r = l;
5333 }
5334 //puts("neg remainder");
5335 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
5336 Polyhedron_Free(T);
5337 }
5338 }
5339 }
5340 if (!(info[i] & ALL_POS) && (info[i] & (ONE_NEG | ROT_NEG)))
5341 goto next;
5342 }
5343 }
5344 if (info[i] & ALL_POS)
5345 break;
5346 next:
5347 ;
5348 }
5349
5350 /*
5351 for (int i = 0; i < exist; ++i)
5352 printf("%i: %i\n", i, info[i]);
5353 */
5354 for (int i = 0; i < exist; ++i)
5355 if (info[i] & ALL_POS) {
5356 #ifdef DEBUG_ER
5357 fprintf(stderr, "\nER: Positive\n");
5358 #endif /* DEBUG_ER */
5359 // Eliminate
5360 // Maybe we should chew off some of the fat here
5361 Matrix *M = Matrix_Alloc(P->Dimension, P->Dimension+1);
5362 for (int j = 0; j < P->Dimension; ++j)
5363 value_set_si(M->p[j][j + (j >= i+nvar)], 1);
5364 Polyhedron *T = Polyhedron_Image(P, M, MaxRays);
5365 Matrix_Free(M);
5366 evalue *EP = barvinok_enumerate_e(T, exist-1, nparam, MaxRays);
5367 Polyhedron_Free(T);
5368 value_clear(f);
5369 Vector_Free(row);
5370 delete [] info;
5371 return EP;
5372 }
5373 for (int i = 0; i < exist; ++i)
5374 if (info[i] & ONE_NEG) {
5375 #ifdef DEBUG_ER
5376 fprintf(stderr, "\nER: Negative\n");
5377 #endif /* DEBUG_ER */
5378 Vector_Free(row);
5379 value_clear(f);
5380 delete [] info;
5381 if (i == 0)
5382 return barvinok_enumerate_e(P, exist-1, nparam, MaxRays);
5383 else {
5384 Polyhedron *T = Polyhedron_Copy(P);
5385 SwapColumns(T, nvar+1, nvar+1+i);
5386 evalue *EP = barvinok_enumerate_e(T, exist-1, nparam, MaxRays);
5387 Polyhedron_Free(T);
5388 return EP;
5389 }
5390 }
5391 for (int i = 0; i < exist; ++i)
5392 if (info[i] & ROT_NEG) {
5393 #ifdef DEBUG_ER
5394 fprintf(stderr, "\nER: Rotate\n");
5395 #endif /* DEBUG_ER */
5396 Vector_Free(row);
5397 value_clear(f);
5398 delete [] info;
5399 Polyhedron *T = rotate_along(P, r, nvar, exist, MaxRays);
5400 evalue *EP = barvinok_enumerate_e(T, exist-1, nparam, MaxRays);
5401 Polyhedron_Free(T);
5402 return EP;
5403 }
5404 for (int i = 0; i < exist; ++i)
5405 if (info[i] & INDEPENDENT) {
5406 Polyhedron *pos, *neg;
5407
5408 /* Find constraint again and split off negative part */
5409
5410 if (SplitOnVar(P, i, nvar, len, exist, MaxRays,
5411 row, f, true, &pos, &neg)) {
5412 #ifdef DEBUG_ER
5413 fprintf(stderr, "\nER: Split\n");
5414 #endif /* DEBUG_ER */
5415
5416 evalue *EP =
5417 barvinok_enumerate_e(neg, exist-1, nparam, MaxRays);
5418 evalue *E =
5419 barvinok_enumerate_e(pos, exist, nparam, MaxRays);
5420 eadd(E, EP);
5421 free_evalue_refs(E);
5422 free(E);
5423 Polyhedron_Free(neg);
5424 Polyhedron_Free(pos);
5425 value_clear(f);
5426 Vector_Free(row);
5427 delete [] info;
5428 return EP;
5429 }
5430 }
5431 delete [] info;
5432
5433 Polyhedron *O = P;
5434 Polyhedron *F;
5435
5436 evalue *EP;
5437
5438 EP = enumerate_line(P, exist, nparam, MaxRays);
5439 if (EP)
5440 goto out;
5441
5442 EP = barvinok_enumerate_pip(P, exist, nparam, MaxRays);
5443 if (EP)
5444 goto out;
5445
5446 EP = enumerate_redundant_ray(P, exist, nparam, MaxRays);
5447 if (EP)
5448 goto out;
5449
5450 EP = enumerate_sure(P, exist, nparam, MaxRays);
5451 if (EP)
5452 goto out;
5453
5454 EP = enumerate_ray(P, exist, nparam, MaxRays);
5455 if (EP)
5456 goto out;
5457
5458 EP = enumerate_sure2(P, exist, nparam, MaxRays);
5459 if (EP)
5460 goto out;
5461
5462 F = unfringe(P, MaxRays);
5463 if (!PolyhedronIncludes(F, P)) {
5464 #ifdef DEBUG_ER
5465 fprintf(stderr, "\nER: Fringed\n");
5466 #endif /* DEBUG_ER */
5467 EP = barvinok_enumerate_e(F, exist, nparam, MaxRays);
5468 Polyhedron_Free(F);
5469 goto out;
5470 }
5471 Polyhedron_Free(F);
5472
5473 if (nparam)
5474 EP = enumerate_vd(&P, exist, nparam, MaxRays);
5475 if (EP)
5476 goto out2;
5477
5478 if (nvar != 0) {
5479 EP = enumerate_sum(P, exist, nparam, MaxRays);
5480 goto out2;
5481 }
5482
5483 assert(nvar == 0);
5484
5485 int i;
5486 Polyhedron *pos, *neg;
5487 for (i = 0; i < exist; ++i)
5488 if (SplitOnVar(P, i, nvar, len, exist, MaxRays,
5489 row, f, false, &pos, &neg))
5490 break;
5491
5492 assert (i < exist);
5493
5494 pos->next = neg;
5495 EP = enumerate_or(pos, exist, nparam, MaxRays);
5496
5497 out2:
5498 if (O != P)
5499 Polyhedron_Free(P);
5500
5501 out:
5502 value_clear(f);
5503 Vector_Free(row);
5504 return EP;
5505 }
5506
5507 gen_fun * barvinok_series(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
5508 {
5509 Polyhedron *CA;
5510 unsigned nparam = C->Dimension;
5511 unsigned dim, nvar;
5512 vec_ZZ sign;
5513 int ncone = 0;
5514 sign.SetLength(ncone);
5515
5516 CA = align_context(C, P->Dimension, MaxRays);
5517 P = DomainIntersection(P, CA, MaxRays);
5518 Polyhedron_Free(CA);
5519
5520 assert(!Polyhedron_is_infinite(P, nparam));
5521 assert(P->NbBid == 0);
5522 assert(Polyhedron_has_positive_rays(P, nparam));
5523 if (P->NbEq != 0)
5524 P = remove_equalities_p(P, P->Dimension-nparam, NULL);
5525 assert(P->NbEq == 0);
5526
5527 #ifdef USE_INCREMENTAL_BF
5528 partial_bfcounter red(P, nparam);
5529 #elif defined USE_INCREMENTAL_DF
5530 partial_ireducer red(P, nparam);
5531 #else
5532 partial_reducer red(P, nparam);
5533 #endif
5534 red.start(MaxRays);
5535 Polyhedron_Free(P);
5536 return red.gf;
5537 }
lattice point enumeration library