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