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