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