?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
add and use incremental counter
[barvinok.git] / barvinok.cc
blobf1cc835bda3f097086b7748e7e6e634292b181a7
1 #include <assert.h>
2 #include <iostream>
3 #include <vector>
4 #include <deque>
5 #include <string>
6 #include <sstream>
7 #include <gmp.h>
8 #include <NTL/mat_ZZ.h>
9 #include <NTL/LLL.h>
10 #include <util.h>
11 extern "C" {
12 #include <polylib/polylibgmp.h>
13 #include "ev_operations.h"
14 #include "piputil.h"
15 }
16 #include "config.h"
17 #include <barvinok.h>
18 #include <genfun.h>
19
20 #ifdef NTL_STD_CXX
21 using namespace NTL;
22 #endif
23 using std::cout;
24 using std::endl;
25 using std::vector;
26 using std::deque;
27 using std::string;
28 using std::ostringstream;
29
30 #define ALLOC(p) (((long *) (p))[0])
31 #define SIZE(p) (((long *) (p))[1])
32 #define DATA(p) ((mp_limb_t *) (((long *) (p)) + 2))
33
34 static void value2zz(Value v, ZZ& z)
35 {
36 int sa = v[0]._mp_size;
37 int abs_sa = sa < 0 ? -sa : sa;
38
39 _ntl_gsetlength(&z.rep, abs_sa);
40 mp_limb_t * adata = DATA(z.rep);
41 for (int i = 0; i < abs_sa; ++i)
42 adata[i] = v[0]._mp_d[i];
43 SIZE(z.rep) = sa;
44 }
45
46 void zz2value(ZZ& z, Value& v)
47 {
48 if (!z.rep) {
49 value_set_si(v, 0);
50 return;
51 }
52
53 int sa = SIZE(z.rep);
54 int abs_sa = sa < 0 ? -sa : sa;
55
56 mp_limb_t * adata = DATA(z.rep);
57 _mpz_realloc(v, abs_sa);
58 for (int i = 0; i < abs_sa; ++i)
59 v[0]._mp_d[i] = adata[i];
60 v[0]._mp_size = sa;
61 }
62
63 #undef ALLOC
64 #define ALLOC(t,p) p = (t*)malloc(sizeof(*p))
65
66 /*
67 * We just ignore the last column and row
68 * If the final element is not equal to one
69 * then the result will actually be a multiple of the input
70 */
71 static void matrix2zz(Matrix *M, mat_ZZ& m, unsigned nr, unsigned nc)
72 {
73 m.SetDims(nr, nc);
74
75 for (int i = 0; i < nr; ++i) {
76 // assert(value_one_p(M->p[i][M->NbColumns - 1]));
77 for (int j = 0; j < nc; ++j) {
78 value2zz(M->p[i][j], m[i][j]);
79 }
80 }
81 }
82
83 static void values2zz(Value *p, vec_ZZ& v, int len)
84 {
85 v.SetLength(len);
86
87 for (int i = 0; i < len; ++i) {
88 value2zz(p[i], v[i]);
89 }
90 }
91
92 /*
93 */
94 static void zz2values(vec_ZZ& v, Value *p)
95 {
96 for (int i = 0; i < v.length(); ++i)
97 zz2value(v[i], p[i]);
98 }
99
100 static void rays(mat_ZZ& r, Polyhedron *C)
101 {
102 unsigned dim = C->NbRays - 1; /* don't count zero vertex */
103 assert(C->NbRays - 1 == C->Dimension);
104 r.SetDims(dim, dim);
105 ZZ tmp;
106
107 int i, c;
108 for (i = 0, c = 0; i < dim; ++i)
109 if (value_zero_p(C->Ray[i][dim+1])) {
110 for (int j = 0; j < dim; ++j) {
111 value2zz(C->Ray[i][j+1], tmp);
112 r[j][c] = tmp;
113 }
114 ++c;
115 }
116 }
117
118 static Matrix * rays(Polyhedron *C)
119 {
120 unsigned dim = C->NbRays - 1; /* don't count zero vertex */
121 assert(C->NbRays - 1 == C->Dimension);
122
123 Matrix *M = Matrix_Alloc(dim+1, dim+1);
124 assert(M);
125
126 int i, c;
127 for (i = 0, c = 0; i <= dim && c < dim; ++i)
128 if (value_zero_p(C->Ray[i][dim+1])) {
129 Vector_Copy(C->Ray[i] + 1, M->p[c], dim);
130 value_set_si(M->p[c++][dim], 0);
131 }
132 assert(c == dim);
133 value_set_si(M->p[dim][dim], 1);
134
135 return M;
136 }
137
138 static Matrix * rays2(Polyhedron *C)
139 {
140 unsigned dim = C->NbRays - 1; /* don't count zero vertex */
141 assert(C->NbRays - 1 == C->Dimension);
142
143 Matrix *M = Matrix_Alloc(dim, dim);
144 assert(M);
145
146 int i, c;
147 for (i = 0, c = 0; i <= dim && c < dim; ++i)
148 if (value_zero_p(C->Ray[i][dim+1]))
149 Vector_Copy(C->Ray[i] + 1, M->p[c++], dim);
150 assert(c == dim);
151
152 return M;
153 }
154
155 /*
156 * Returns the largest absolute value in the vector
157 */
158 static ZZ max(vec_ZZ& v)
159 {
160 ZZ max = abs(v[0]);
161 for (int i = 1; i < v.length(); ++i)
162 if (abs(v[i]) > max)
163 max = abs(v[i]);
164 return max;
165 }
166
167 class cone {
168 public:
169 cone(Matrix *M) {
170 Cone = 0;
171 Rays = Matrix_Copy(M);
172 set_det();
173 }
174 cone(Polyhedron *C) {
175 Cone = Polyhedron_Copy(C);
176 Rays = rays(C);
177 set_det();
178 }
179 void set_det() {
180 mat_ZZ A;
181 matrix2zz(Rays, A, Rays->NbRows - 1, Rays->NbColumns - 1);
182 det = determinant(A);
183 }
184
185 Vector* short_vector(vec_ZZ& lambda) {
186 Matrix *M = Matrix_Copy(Rays);
187 Matrix *inv = Matrix_Alloc(M->NbRows, M->NbColumns);
188 int ok = Matrix_Inverse(M, inv);
189 assert(ok);
190 Matrix_Free(M);
191
192 ZZ det2;
193 mat_ZZ B;
194 mat_ZZ U;
195 matrix2zz(inv, B, inv->NbRows - 1, inv->NbColumns - 1);
196 long r = LLL(det2, B, U);
197
198 ZZ min = max(B[0]);
199 int index = 0;
200 for (int i = 1; i < B.NumRows(); ++i) {
201 ZZ tmp = max(B[i]);
202 if (tmp < min) {
203 min = tmp;
204 index = i;
205 }
206 }
207
208 Matrix_Free(inv);
209
210 lambda = B[index];
211
212 Vector *z = Vector_Alloc(U[index].length()+1);
213 assert(z);
214 zz2values(U[index], z->p);
215 value_set_si(z->p[U[index].length()], 0);
216
217 Value tmp;
218 value_init(tmp);
219 Polyhedron *C = poly();
220 int i;
221 for (i = 0; i < C->NbConstraints; ++i) {
222 Inner_Product(z->p, C->Constraint[i]+1, z->Size-1, &tmp);
223 if (value_pos_p(tmp))
224 break;
225 }
226 if (i == C->NbConstraints) {
227 value_set_si(tmp, -1);
228 Vector_Scale(z->p, z->p, tmp, z->Size-1);
229 }
230 value_clear(tmp);
231 return z;
232 }
233
234 ~cone() {
235 Polyhedron_Free(Cone);
236 Matrix_Free(Rays);
237 }
238
239 Polyhedron *poly() {
240 if (!Cone) {
241 Matrix *M = Matrix_Alloc(Rays->NbRows+1, Rays->NbColumns+1);
242 for (int i = 0; i < Rays->NbRows; ++i) {
243 Vector_Copy(Rays->p[i], M->p[i]+1, Rays->NbColumns);
244 value_set_si(M->p[i][0], 1);
245 }
246 Vector_Set(M->p[Rays->NbRows]+1, 0, Rays->NbColumns-1);
247 value_set_si(M->p[Rays->NbRows][0], 1);
248 value_set_si(M->p[Rays->NbRows][Rays->NbColumns], 1);
249 Cone = Rays2Polyhedron(M, M->NbRows+1);
250 assert(Cone->NbConstraints == Cone->NbRays);
251 Matrix_Free(M);
252 }
253 return Cone;
254 }
255
256 ZZ det;
257 Polyhedron *Cone;
258 Matrix *Rays;
259 };
260
261 class dpoly {
262 public:
263 vec_ZZ coeff;
264 dpoly(int d, ZZ& degree, int offset = 0) {
265 coeff.SetLength(d+1);
266
267 int min = d + offset;
268 if (degree >= 0 && degree < ZZ(INIT_VAL, min))
269 min = to_int(degree);
270
271 ZZ c = ZZ(INIT_VAL, 1);
272 if (!offset)
273 coeff[0] = c;
274 for (int i = 1; i <= min; ++i) {
275 c *= (degree -i + 1);
276 c /= i;
277 coeff[i-offset] = c;
278 }
279 }
280 void operator *= (dpoly& f) {
281 assert(coeff.length() == f.coeff.length());
282 vec_ZZ old = coeff;
283 coeff = f.coeff[0] * coeff;
284 for (int i = 1; i < coeff.length(); ++i)
285 for (int j = 0; i+j < coeff.length(); ++j)
286 coeff[i+j] += f.coeff[i] * old[j];
287 }
288 void div(dpoly& d, mpq_t count, ZZ& sign) {
289 int len = coeff.length();
290 Value tmp;
291 value_init(tmp);
292 mpq_t* c = new mpq_t[coeff.length()];
293 mpq_t qtmp;
294 mpq_init(qtmp);
295 for (int i = 0; i < len; ++i) {
296 mpq_init(c[i]);
297 zz2value(coeff[i], tmp);
298 mpq_set_z(c[i], tmp);
299
300 for (int j = 1; j <= i; ++j) {
301 zz2value(d.coeff[j], tmp);
302 mpq_set_z(qtmp, tmp);
303 mpq_mul(qtmp, qtmp, c[i-j]);
304 mpq_sub(c[i], c[i], qtmp);
305 }
306
307 zz2value(d.coeff[0], tmp);
308 mpq_set_z(qtmp, tmp);
309 mpq_div(c[i], c[i], qtmp);
310 }
311 if (sign == -1)
312 mpq_sub(count, count, c[len-1]);
313 else
314 mpq_add(count, count, c[len-1]);
315
316 value_clear(tmp);
317 mpq_clear(qtmp);
318 for (int i = 0; i < len; ++i)
319 mpq_clear(c[i]);
320 delete [] c;
321 }
322 };
323
324 class dpoly_n {
325 public:
326 Matrix *coeff;
327 ~dpoly_n() {
328 Matrix_Free(coeff);
329 }
330 dpoly_n(int d, ZZ& degree_0, ZZ& degree_1, int offset = 0) {
331 Value d0, d1;
332 value_init(d0);
333 value_init(d1);
334 zz2value(degree_0, d0);
335 zz2value(degree_1, d1);
336 coeff = Matrix_Alloc(d+1, d+1+1);
337 value_set_si(coeff->p[0][0], 1);
338 value_set_si(coeff->p[0][d+1], 1);
339 for (int i = 1; i <= d; ++i) {
340 value_multiply(coeff->p[i][0], coeff->p[i-1][0], d0);
341 Vector_Combine(coeff->p[i-1], coeff->p[i-1]+1, coeff->p[i]+1,
342 d1, d0, i);
343 value_set_si(coeff->p[i][d+1], i);
344 value_multiply(coeff->p[i][d+1], coeff->p[i][d+1], coeff->p[i-1][d+1]);
345 value_decrement(d0, d0);
346 }
347 value_clear(d0);
348 value_clear(d1);
349 }
350 void div(dpoly& d, Vector *count, ZZ& sign) {
351 int len = coeff->NbRows;
352 Matrix * c = Matrix_Alloc(coeff->NbRows, coeff->NbColumns);
353 Value tmp;
354 value_init(tmp);
355 for (int i = 0; i < len; ++i) {
356 Vector_Copy(coeff->p[i], c->p[i], len+1);
357 for (int j = 1; j <= i; ++j) {
358 zz2value(d.coeff[j], tmp);
359 value_multiply(tmp, tmp, c->p[i][len]);
360 value_oppose(tmp, tmp);
361 Vector_Combine(c->p[i], c->p[i-j], c->p[i],
362 c->p[i-j][len], tmp, len);
363 value_multiply(c->p[i][len], c->p[i][len], c->p[i-j][len]);
364 }
365 zz2value(d.coeff[0], tmp);
366 value_multiply(c->p[i][len], c->p[i][len], tmp);
367 }
368 if (sign == -1) {
369 value_set_si(tmp, -1);
370 Vector_Scale(c->p[len-1], count->p, tmp, len);
371 value_assign(count->p[len], c->p[len-1][len]);
372 } else
373 Vector_Copy(c->p[len-1], count->p, len+1);
374 Vector_Normalize(count->p, len+1);
375 value_clear(tmp);
376 Matrix_Free(c);
377 }
378 };
379
380 struct dpoly_r_term {
381 int *powers;
382 ZZ coeff;
383 };
384
385 /* len: number of elements in c
386 * each element in c is the coefficient of a power of t
387 * in the MacLaurin expansion
388 */
389 struct dpoly_r {
390 vector< dpoly_r_term * > *c;
391 int len;
392 int dim;
393 ZZ denom;
394
395 void add_term(int i, int * powers, ZZ& coeff) {
396 for (int k = 0; k < c[i].size(); ++k) {
397 if (memcmp(c[i][k]->powers, powers, dim * sizeof(int)) == 0) {
398 c[i][k]->coeff += coeff;
399 return;
400 }
401 }
402 dpoly_r_term *t = new dpoly_r_term;
403 t->powers = new int[dim];
404 memcpy(t->powers, powers, dim * sizeof(int));
405 t->coeff = coeff;
406 c[i].push_back(t);
407 }
408 dpoly_r(int len, int dim) {
409 denom = 1;
410 this->len = len;
411 this->dim = dim;
412 c = new vector< dpoly_r_term * > [len];
413 }
414 dpoly_r(dpoly& num, dpoly& den, int pos, int sign, int dim) {
415 denom = 1;
416 len = num.coeff.length();
417 c = new vector< dpoly_r_term * > [len];
418 this->dim = dim;
419 int powers[dim];
420
421 for (int i = 0; i < len; ++i) {
422 ZZ coeff = num.coeff[i];
423 memset(powers, 0, dim * sizeof(int));
424 powers[pos] = sign;
425
426 add_term(i, powers, coeff);
427
428 for (int j = 1; j <= i; ++j) {
429 for (int k = 0; k < c[i-j].size(); ++k) {
430 memcpy(powers, c[i-j][k]->powers, dim*sizeof(int));
431 powers[pos] += sign;
432 coeff = -den.coeff[j-1] * c[i-j][k]->coeff;
433 add_term(i, powers, coeff);
434 }
435 }
436 }
437 //dump();
438 }
439 dpoly_r(dpoly_r* num, dpoly& den, int pos, int sign, int dim) {
440 denom = num->denom;
441 len = num->len;
442 c = new vector< dpoly_r_term * > [len];
443 this->dim = dim;
444 int powers[dim];
445 ZZ coeff;
446
447 for (int i = 0 ; i < len; ++i) {
448 for (int k = 0; k < num->c[i].size(); ++k) {
449 memcpy(powers, num->c[i][k]->powers, dim*sizeof(int));
450 powers[pos] += sign;
451 add_term(i, powers, num->c[i][k]->coeff);
452 }
453
454 for (int j = 1; j <= i; ++j) {
455 for (int k = 0; k < c[i-j].size(); ++k) {
456 memcpy(powers, c[i-j][k]->powers, dim*sizeof(int));
457 powers[pos] += sign;
458 coeff = -den.coeff[j-1] * c[i-j][k]->coeff;
459 add_term(i, powers, coeff);
460 }
461 }
462 }
463 }
464 ~dpoly_r() {
465 for (int i = 0 ; i < len; ++i)
466 for (int k = 0; k < c[i].size(); ++k) {
467 delete [] c[i][k]->powers;
468 delete c[i][k];
469 }
470 delete [] c;
471 }
472 dpoly_r *div(dpoly& d) {
473 dpoly_r *rc = new dpoly_r(len, dim);
474 rc->denom = power(d.coeff[0], len);
475 ZZ inv_d = rc->denom / d.coeff[0];
476 ZZ coeff;
477
478 for (int i = 0; i < len; ++i) {
479 for (int k = 0; k < c[i].size(); ++k) {
480 coeff = c[i][k]->coeff * inv_d;
481 rc->add_term(i, c[i][k]->powers, coeff);
482 }
483
484 for (int j = 1; j <= i; ++j) {
485 for (int k = 0; k < rc->c[i-j].size(); ++k) {
486 coeff = - d.coeff[j] * rc->c[i-j][k]->coeff / d.coeff[0];
487 rc->add_term(i, rc->c[i-j][k]->powers, coeff);
488 }
489 }
490 }
491 return rc;
492 }
493 void div(dpoly& d, ZZ& sign, gen_fun *gf, mat_ZZ& pden, mat_ZZ& den,
494 vec_ZZ& num_p) {
495 dpoly_r * rc = div(d);
496 //rc.dump();
497 int common = pden.NumRows();
498
499 vector< dpoly_r_term * >& final = rc->c[len-1];
500 int rows;
501 for (int j = 0; j < final.size(); ++j) {
502 rows = common;
503 pden.SetDims(rows, pden.NumCols());
504 for (int k = 0; k < dim; ++k) {
505 int n = final[j]->powers[k];
506 if (n == 0)
507 continue;
508 int abs_n = n < 0 ? -n : n;
509 pden.SetDims(rows+abs_n, pden.NumCols());
510 for (int l = 0; l < abs_n; ++l) {
511 if (n > 0)
512 pden[rows+l] = den[k];
513 else
514 pden[rows+l] = -den[k];
515 }
516 rows += abs_n;
517 }
518 final[j]->coeff *= sign;
519 gf->add(final[j]->coeff, rc->denom, num_p, pden);
520 }
521 delete rc;
522 }
523 void dump(void) {
524 for (int i = 0; i < len; ++i) {
525 cout << endl;
526 cout << i << endl;
527 cout << c[i].size() << endl;
528 for (int j = 0; j < c[i].size(); ++j) {
529 for (int k = 0; k < dim; ++k) {
530 cout << c[i][j]->powers[k] << " ";
531 }
532 cout << ": " << c[i][j]->coeff << "/" << denom << endl;
533 }
534 cout << endl;
535 }
536 }
537 };
538
539 struct decomposer {
540 void decompose(Polyhedron *C);
541 virtual void handle(Polyhedron *P, int sign) = 0;
542 };
543
544 struct polar_decomposer : public decomposer {
545 void decompose(Polyhedron *C, unsigned MaxRays);
546 virtual void handle(Polyhedron *P, int sign);
547 virtual void handle_polar(Polyhedron *P, int sign) = 0;
548 };
549
550 void decomposer::decompose(Polyhedron *C)
551 {
552 vector<cone *> nonuni;
553 cone * c = new cone(C);
554 ZZ det = c->det;
555 int s = sign(det);
556 assert(det != 0);
557 if (abs(det) > 1) {
558 nonuni.push_back(c);
559 } else {
560 handle(C, 1);
561 delete c;
562 }
563 vec_ZZ lambda;
564 while (!nonuni.empty()) {
565 c = nonuni.back();
566 nonuni.pop_back();
567 Vector* v = c->short_vector(lambda);
568 for (int i = 0; i < c->Rays->NbRows - 1; ++i) {
569 if (lambda[i] == 0)
570 continue;
571 Matrix* M = Matrix_Copy(c->Rays);
572 Vector_Copy(v->p, M->p[i], v->Size);
573 cone * pc = new cone(M);
574 assert (pc->det != 0);
575 if (abs(pc->det) > 1) {
576 assert(abs(pc->det) < abs(c->det));
577 nonuni.push_back(pc);
578 } else {
579 handle(pc->poly(), sign(pc->det) * s);
580 delete pc;
581 }
582 Matrix_Free(M);
583 }
584 Vector_Free(v);
585 delete c;
586 }
587 }
588
589 void polar_decomposer::decompose(Polyhedron *cone, unsigned MaxRays)
590 {
591 Polyhedron_Polarize(cone);
592 if (cone->NbRays - 1 != cone->Dimension) {
593 Polyhedron *tmp = cone;
594 cone = triangularize_cone(cone, MaxRays);
595 Polyhedron_Free(tmp);
596 }
597 for (Polyhedron *Polar = cone; Polar; Polar = Polar->next)
598 decomposer::decompose(Polar);
599 Domain_Free(cone);
600 }
601
602 void polar_decomposer::handle(Polyhedron *P, int sign)
603 {
604 Polyhedron_Polarize(P);
605 handle_polar(P, sign);
606 }
607
608 /*
609 * Barvinok's Decomposition of a simplicial cone
610 *
611 * Returns two lists of polyhedra
612 */
613 void barvinok_decompose(Polyhedron *C, Polyhedron **ppos, Polyhedron **pneg)
614 {
615 Polyhedron *pos = *ppos, *neg = *pneg;
616 vector<cone *> nonuni;
617 cone * c = new cone(C);
618 ZZ det = c->det;
619 int s = sign(det);
620 assert(det != 0);
621 if (abs(det) > 1) {
622 nonuni.push_back(c);
623 } else {
624 Polyhedron *p = Polyhedron_Copy(c->Cone);
625 p->next = pos;
626 pos = p;
627 delete c;
628 }
629 vec_ZZ lambda;
630 while (!nonuni.empty()) {
631 c = nonuni.back();
632 nonuni.pop_back();
633 Vector* v = c->short_vector(lambda);
634 for (int i = 0; i < c->Rays->NbRows - 1; ++i) {
635 if (lambda[i] == 0)
636 continue;
637 Matrix* M = Matrix_Copy(c->Rays);
638 Vector_Copy(v->p, M->p[i], v->Size);
639 cone * pc = new cone(M);
640 assert (pc->det != 0);
641 if (abs(pc->det) > 1) {
642 assert(abs(pc->det) < abs(c->det));
643 nonuni.push_back(pc);
644 } else {
645 Polyhedron *p = pc->poly();
646 pc->Cone = 0;
647 if (sign(pc->det) == s) {
648 p->next = pos;
649 pos = p;
650 } else {
651 p->next = neg;
652 neg = p;
653 }
654 delete pc;
655 }
656 Matrix_Free(M);
657 }
658 Vector_Free(v);
659 delete c;
660 }
661 *ppos = pos;
662 *pneg = neg;
663 }
664
665 /*
666 * Returns a single list of npos "positive" cones followed by nneg
667 * "negative" cones.
668 * The input cone is freed
669 */
670 void decompose(Polyhedron *cone, Polyhedron **parts, int *npos, int *nneg, unsigned MaxRays)
671 {
672 Polyhedron_Polarize(cone);
673 if (cone->NbRays - 1 != cone->Dimension) {
674 Polyhedron *tmp = cone;
675 cone = triangularize_cone(cone, MaxRays);
676 Polyhedron_Free(tmp);
677 }
678 Polyhedron *polpos = NULL, *polneg = NULL;
679 *npos = 0; *nneg = 0;
680 for (Polyhedron *Polar = cone; Polar; Polar = Polar->next)
681 barvinok_decompose(Polar, &polpos, &polneg);
682
683 Polyhedron *last;
684 for (Polyhedron *i = polpos; i; i = i->next) {
685 Polyhedron_Polarize(i);
686 ++*npos;
687 last = i;
688 }
689 for (Polyhedron *i = polneg; i; i = i->next) {
690 Polyhedron_Polarize(i);
691 ++*nneg;
692 }
693 if (last) {
694 last->next = polneg;
695 *parts = polpos;
696 } else
697 *parts = polneg;
698 Domain_Free(cone);
699 }
700
701 const int MAX_TRY=10;
702 /*
703 * Searches for a vector that is not orthogonal to any
704 * of the rays in rays.
705 */
706 static void nonorthog(mat_ZZ& rays, vec_ZZ& lambda)
707 {
708 int dim = rays.NumCols();
709 bool found = false;
710 lambda.SetLength(dim);
711 if (dim == 0)
712 return;
713
714 for (int i = 2; !found && i <= 50*dim; i+=4) {
715 for (int j = 0; j < MAX_TRY; ++j) {
716 for (int k = 0; k < dim; ++k) {
717 int r = random_int(i)+2;
718 int v = (2*(r%2)-1) * (r >> 1);
719 lambda[k] = v;
720 }
721 int k = 0;
722 for (; k < rays.NumRows(); ++k)
723 if (lambda * rays[k] == 0)
724 break;
725 if (k == rays.NumRows()) {
726 found = true;
727 break;
728 }
729 }
730 }
731 assert(found);
732 }
733
734 static void randomvector(Polyhedron *P, vec_ZZ& lambda, int nvar)
735 {
736 Value tmp;
737 int max = 10 * 16;
738 unsigned int dim = P->Dimension;
739 value_init(tmp);
740
741 for (int i = 0; i < P->NbRays; ++i) {
742 for (int j = 1; j <= dim; ++j) {
743 value_absolute(tmp, P->Ray[i][j]);
744 int t = VALUE_TO_LONG(tmp) * 16;
745 if (t > max)
746 max = t;
747 }
748 }
749 for (int i = 0; i < P->NbConstraints; ++i) {
750 for (int j = 1; j <= dim; ++j) {
751 value_absolute(tmp, P->Constraint[i][j]);
752 int t = VALUE_TO_LONG(tmp) * 16;
753 if (t > max)
754 max = t;
755 }
756 }
757 value_clear(tmp);
758
759 lambda.SetLength(nvar);
760 for (int k = 0; k < nvar; ++k) {
761 int r = random_int(max*dim)+2;
762 int v = (2*(r%2)-1) * (max/2*dim + (r >> 1));
763 lambda[k] = v;
764 }
765 }
766
767 static void add_rays(mat_ZZ& rays, Polyhedron *i, int *r, int nvar = -1,
768 bool all = false)
769 {
770 unsigned dim = i->Dimension;
771 if (nvar == -1)
772 nvar = dim;
773 for (int k = 0; k < i->NbRays; ++k) {
774 if (!value_zero_p(i->Ray[k][dim+1]))
775 continue;
776 if (!all && nvar != dim && First_Non_Zero(i->Ray[k]+1, nvar) == -1)
777 continue;
778 values2zz(i->Ray[k]+1, rays[(*r)++], nvar);
779 }
780 }
781
782 void lattice_point(Value* values, Polyhedron *i, vec_ZZ& vertex)
783 {
784 unsigned dim = i->Dimension;
785 if(!value_one_p(values[dim])) {
786 Matrix* Rays = rays(i);
787 Matrix *inv = Matrix_Alloc(Rays->NbRows, Rays->NbColumns);
788 int ok = Matrix_Inverse(Rays, inv);
789 assert(ok);
790 Matrix_Free(Rays);
791 Rays = rays(i);
792 Vector *lambda = Vector_Alloc(dim+1);
793 Vector_Matrix_Product(values, inv, lambda->p);
794 Matrix_Free(inv);
795 for (int j = 0; j < dim; ++j)
796 mpz_cdiv_q(lambda->p[j], lambda->p[j], lambda->p[dim]);
797 value_set_si(lambda->p[dim], 1);
798 Vector *A = Vector_Alloc(dim+1);
799 Vector_Matrix_Product(lambda->p, Rays, A->p);
800 Vector_Free(lambda);
801 Matrix_Free(Rays);
802 values2zz(A->p, vertex, dim);
803 Vector_Free(A);
804 } else
805 values2zz(values, vertex, dim);
806 }
807
808 static evalue *term(int param, ZZ& c, Value *den = NULL)
809 {
810 evalue *EP = new evalue();
811 value_init(EP->d);
812 value_set_si(EP->d,0);
813 EP->x.p = new_enode(polynomial, 2, param + 1);
814 evalue_set_si(&EP->x.p->arr[0], 0, 1);
815 value_init(EP->x.p->arr[1].x.n);
816 if (den == NULL)
817 value_set_si(EP->x.p->arr[1].d, 1);
818 else
819 value_assign(EP->x.p->arr[1].d, *den);
820 zz2value(c, EP->x.p->arr[1].x.n);
821 return EP;
822 }
823
824 static void vertex_period(
825 Polyhedron *i, vec_ZZ& lambda, Matrix *T,
826 Value lcm, int p, Vector *val,
827 evalue *E, evalue* ev,
828 ZZ& offset)
829 {
830 unsigned nparam = T->NbRows - 1;
831 unsigned dim = i->Dimension;
832 Value tmp;
833 ZZ nump;
834
835 if (p == nparam) {
836 vec_ZZ vertex;
837 ZZ num, l;
838 Vector * values = Vector_Alloc(dim + 1);
839 Vector_Matrix_Product(val->p, T, values->p);
840 value_assign(values->p[dim], lcm);
841 lattice_point(values->p, i, vertex);
842 num = vertex * lambda;
843 value2zz(lcm, l);
844 num *= l;
845 num += offset;
846 value_init(ev->x.n);
847 zz2value(num, ev->x.n);
848 value_assign(ev->d, lcm);
849 Vector_Free(values);
850 return;
851 }
852
853 value_init(tmp);
854 vec_ZZ vertex;
855 values2zz(T->p[p], vertex, dim);
856 nump = vertex * lambda;
857 if (First_Non_Zero(val->p, p) == -1) {
858 value_assign(tmp, lcm);
859 evalue *ET = term(p, nump, &tmp);
860 eadd(ET, E);
861 free_evalue_refs(ET);
862 delete ET;
863 }
864
865 value_assign(tmp, lcm);
866 if (First_Non_Zero(T->p[p], dim) != -1)
867 Vector_Gcd(T->p[p], dim, &tmp);
868 Gcd(tmp, lcm, &tmp);
869 if (value_lt(tmp, lcm)) {
870 ZZ count;
871
872 value_division(tmp, lcm, tmp);
873 value_set_si(ev->d, 0);
874 ev->x.p = new_enode(periodic, VALUE_TO_INT(tmp), p+1);
875 value2zz(tmp, count);
876 do {
877 value_decrement(tmp, tmp);
878 --count;
879 ZZ new_offset = offset - count * nump;
880 value_assign(val->p[p], tmp);
881 vertex_period(i, lambda, T, lcm, p+1, val, E,
882 &ev->x.p->arr[VALUE_TO_INT(tmp)], new_offset);
883 } while (value_pos_p(tmp));
884 } else
885 vertex_period(i, lambda, T, lcm, p+1, val, E, ev, offset);
886 value_clear(tmp);
887 }
888
889 static void mask_r(Matrix *f, int nr, Vector *lcm, int p, Vector *val, evalue *ev)
890 {
891 unsigned nparam = lcm->Size;
892
893 if (p == nparam) {
894 Vector * prod = Vector_Alloc(f->NbRows);
895 Matrix_Vector_Product(f, val->p, prod->p);
896 int isint = 1;
897 for (int i = 0; i < nr; ++i) {
898 value_modulus(prod->p[i], prod->p[i], f->p[i][nparam+1]);
899 isint &= value_zero_p(prod->p[i]);
900 }
901 value_set_si(ev->d, 1);
902 value_init(ev->x.n);
903 value_set_si(ev->x.n, isint);
904 Vector_Free(prod);
905 return;
906 }
907
908 Value tmp;
909 value_init(tmp);
910 if (value_one_p(lcm->p[p]))
911 mask_r(f, nr, lcm, p+1, val, ev);
912 else {
913 value_assign(tmp, lcm->p[p]);
914 value_set_si(ev->d, 0);
915 ev->x.p = new_enode(periodic, VALUE_TO_INT(tmp), p+1);
916 do {
917 value_decrement(tmp, tmp);
918 value_assign(val->p[p], tmp);
919 mask_r(f, nr, lcm, p+1, val, &ev->x.p->arr[VALUE_TO_INT(tmp)]);
920 } while (value_pos_p(tmp));
921 }
922 value_clear(tmp);
923 }
924
925 static evalue *multi_monom(vec_ZZ& p)
926 {
927 evalue *X = new evalue();
928 value_init(X->d);
929 value_init(X->x.n);
930 unsigned nparam = p.length()-1;
931 zz2value(p[nparam], X->x.n);
932 value_set_si(X->d, 1);
933 for (int i = 0; i < nparam; ++i) {
934 if (p[i] == 0)
935 continue;
936 evalue *T = term(i, p[i]);
937 eadd(T, X);
938 free_evalue_refs(T);
939 delete T;
940 }
941 return X;
942 }
943
944 /*
945 * Check whether mapping polyhedron P on the affine combination
946 * num yields a range that has a fixed quotient on integer
947 * division by d
948 * If zero is true, then we are only interested in the quotient
949 * for the cases where the remainder is zero.
950 * Returns NULL if false and a newly allocated value if true.
951 */
952 static Value *fixed_quotient(Polyhedron *P, vec_ZZ& num, Value d, bool zero)
953 {
954 Value* ret = NULL;
955 int len = num.length();
956 Matrix *T = Matrix_Alloc(2, len);
957 zz2values(num, T->p[0]);
958 value_set_si(T->p[1][len-1], 1);
959 Polyhedron *I = Polyhedron_Image(P, T, P->NbConstraints);
960 Matrix_Free(T);
961
962 int i;
963 for (i = 0; i < I->NbRays; ++i)
964 if (value_zero_p(I->Ray[i][2])) {
965 Polyhedron_Free(I);
966 return NULL;
967 }
968
969 Value min, max;
970 value_init(min);
971 value_init(max);
972 int bounded = line_minmax(I, &min, &max);
973 assert(bounded);
974
975 if (zero)
976 mpz_cdiv_q(min, min, d);
977 else
978 mpz_fdiv_q(min, min, d);
979 mpz_fdiv_q(max, max, d);
980
981 if (value_eq(min, max)) {
982 ALLOC(Value, ret);
983 value_init(*ret);
984 value_assign(*ret, min);
985 }
986 value_clear(min);
987 value_clear(max);
988 return ret;
989 }
990
991 /*
992 * Normalize linear expression coef modulo m
993 * Removes common factor and reduces coefficients
994 * Returns index of first non-zero coefficient or len
995 */
996 static int normal_mod(Value *coef, int len, Value *m)
997 {
998 Value gcd;
999 value_init(gcd);
1000
1001 Vector_Gcd(coef, len, &gcd);
1002 Gcd(gcd, *m, &gcd);
1003 Vector_AntiScale(coef, coef, gcd, len);
1004
1005 value_division(*m, *m, gcd);
1006 value_clear(gcd);
1007
1008 if (value_one_p(*m))
1009 return len;
1010
1011 int j;
1012 for (j = 0; j < len; ++j)
1013 mpz_fdiv_r(coef[j], coef[j], *m);
1014 for (j = 0; j < len; ++j)
1015 if (value_notzero_p(coef[j]))
1016 break;
1017
1018 return j;
1019 }
1020
1021 #ifdef USE_MODULO
1022 static void mask(Matrix *f, evalue *factor)
1023 {
1024 int nr = f->NbRows, nc = f->NbColumns;
1025 int n;
1026 bool found = false;
1027 for (n = 0; n < nr && value_notzero_p(f->p[n][nc-1]); ++n)
1028 if (value_notone_p(f->p[n][nc-1]) &&
1029 value_notmone_p(f->p[n][nc-1]))
1030 found = true;
1031 if (!found)
1032 return;
1033
1034 evalue EP;
1035 nr = n;
1036
1037 Value m;
1038 value_init(m);
1039
1040 evalue EV;
1041 value_init(EV.d);
1042 value_init(EV.x.n);
1043 value_set_si(EV.x.n, 1);
1044
1045 for (n = 0; n < nr; ++n) {
1046 value_assign(m, f->p[n][nc-1]);
1047 if (value_one_p(m) || value_mone_p(m))
1048 continue;
1049
1050 int j = normal_mod(f->p[n], nc-1, &m);
1051 if (j == nc-1) {
1052 free_evalue_refs(factor);
1053 value_init(factor->d);
1054 evalue_set_si(factor, 0, 1);
1055 break;
1056 }
1057 vec_ZZ row;
1058 values2zz(f->p[n], row, nc-1);
1059 ZZ g;
1060 value2zz(m, g);
1061 if (j < (nc-1)-1 && row[j] > g/2) {
1062 for (int k = j; k < (nc-1); ++k)
1063 if (row[k] != 0)
1064 row[k] = g - row[k];
1065 }
1066
1067 value_init(EP.d);
1068 value_set_si(EP.d, 0);
1069 EP.x.p = new_enode(relation, 2, 0);
1070 value_clear(EP.x.p->arr[1].d);
1071 EP.x.p->arr[1] = *factor;
1072 evalue *ev = &EP.x.p->arr[0];
1073 value_set_si(ev->d, 0);
1074 ev->x.p = new_enode(fractional, 3, -1);
1075 evalue_set_si(&ev->x.p->arr[1], 0, 1);
1076 evalue_set_si(&ev->x.p->arr[2], 1, 1);
1077 evalue *E = multi_monom(row);
1078 value_assign(EV.d, m);
1079 emul(&EV, E);
1080 value_clear(ev->x.p->arr[0].d);
1081 ev->x.p->arr[0] = *E;
1082 delete E;
1083 *factor = EP;
1084 }
1085
1086 value_clear(m);
1087 free_evalue_refs(&EV);
1088 }
1089 #else
1090 /*
1091 *
1092 */
1093 static void mask(Matrix *f, evalue *factor)
1094 {
1095 int nr = f->NbRows, nc = f->NbColumns;
1096 int n;
1097 bool found = false;
1098 for (n = 0; n < nr && value_notzero_p(f->p[n][nc-1]); ++n)
1099 if (value_notone_p(f->p[n][nc-1]) &&
1100 value_notmone_p(f->p[n][nc-1]))
1101 found = true;
1102 if (!found)
1103 return;
1104
1105 Value tmp;
1106 value_init(tmp);
1107 nr = n;
1108 unsigned np = nc - 2;
1109 Vector *lcm = Vector_Alloc(np);
1110 Vector *val = Vector_Alloc(nc);
1111 Vector_Set(val->p, 0, nc);
1112 value_set_si(val->p[np], 1);
1113 Vector_Set(lcm->p, 1, np);
1114 for (n = 0; n < nr; ++n) {
1115 if (value_one_p(f->p[n][nc-1]) ||
1116 value_mone_p(f->p[n][nc-1]))
1117 continue;
1118 for (int j = 0; j < np; ++j)
1119 if (value_notzero_p(f->p[n][j])) {
1120 Gcd(f->p[n][j], f->p[n][nc-1], &tmp);
1121 value_division(tmp, f->p[n][nc-1], tmp);
1122 value_lcm(tmp, lcm->p[j], &lcm->p[j]);
1123 }
1124 }
1125 evalue EP;
1126 value_init(EP.d);
1127 mask_r(f, nr, lcm, 0, val, &EP);
1128 value_clear(tmp);
1129 Vector_Free(val);
1130 Vector_Free(lcm);
1131 emul(&EP,factor);
1132 free_evalue_refs(&EP);
1133 }
1134 #endif
1135
1136 struct term_info {
1137 evalue *E;
1138 ZZ constant;
1139 ZZ coeff;
1140 int pos;
1141 };
1142
1143 static bool mod_needed(Polyhedron *PD, vec_ZZ& num, Value d, evalue *E)
1144 {
1145 Value *q = fixed_quotient(PD, num, d, false);
1146
1147 if (!q)
1148 return true;
1149
1150 value_oppose(*q, *q);
1151 evalue EV;
1152 value_init(EV.d);
1153 value_set_si(EV.d, 1);
1154 value_init(EV.x.n);
1155 value_multiply(EV.x.n, *q, d);
1156 eadd(&EV, E);
1157 free_evalue_refs(&EV);
1158 value_clear(*q);
1159 free(q);
1160 return false;
1161 }
1162
1163 static void ceil_mod(Value *coef, int len, Value d, ZZ& f, evalue *EP, Polyhedron *PD)
1164 {
1165 Value m;
1166 value_init(m);
1167 value_set_si(m, -1);
1168
1169 Vector_Scale(coef, coef, m, len);
1170
1171 value_assign(m, d);
1172 int j = normal_mod(coef, len, &m);
1173
1174 if (j == len) {
1175 value_clear(m);
1176 return;
1177 }
1178
1179 vec_ZZ num;
1180 values2zz(coef, num, len);
1181
1182 ZZ g;
1183 value2zz(m, g);
1184
1185 evalue tmp;
1186 value_init(tmp.d);
1187 evalue_set_si(&tmp, 0, 1);
1188
1189 int p = j;
1190 if (g % 2 == 0)
1191 while (j < len-1 && (num[j] == g/2 || num[j] == 0))
1192 ++j;
1193 if ((j < len-1 && num[j] > g/2) || (j == len-1 && num[j] >= (g+1)/2)) {
1194 for (int k = j; k < len-1; ++k)
1195 if (num[k] != 0)
1196 num[k] = g - num[k];
1197 num[len-1] = g - 1 - num[len-1];
1198 value_assign(tmp.d, m);
1199 ZZ t = f*(g-1);
1200 zz2value(t, tmp.x.n);
1201 eadd(&tmp, EP);
1202 f = -f;
1203 }
1204
1205 if (p >= len-1) {
1206 ZZ t = num[len-1] * f;
1207 zz2value(t, tmp.x.n);
1208 value_assign(tmp.d, m);
1209 eadd(&tmp, EP);
1210 } else {
1211 evalue *E = multi_monom(num);
1212 evalue EV;
1213 value_init(EV.d);
1214
1215 if (PD && !mod_needed(PD, num, m, E)) {
1216 value_init(EV.x.n);
1217 zz2value(f, EV.x.n);
1218 value_assign(EV.d, m);
1219 emul(&EV, E);
1220 eadd(E, EP);
1221 } else {
1222 value_init(EV.x.n);
1223 value_set_si(EV.x.n, 1);
1224 value_assign(EV.d, m);
1225 emul(&EV, E);
1226 value_clear(EV.x.n);
1227 value_set_si(EV.d, 0);
1228 EV.x.p = new_enode(fractional, 3, -1);
1229 evalue_copy(&EV.x.p->arr[0], E);
1230 evalue_set_si(&EV.x.p->arr[1], 0, 1);
1231 value_init(EV.x.p->arr[2].x.n);
1232 zz2value(f, EV.x.p->arr[2].x.n);
1233 value_set_si(EV.x.p->arr[2].d, 1);
1234
1235 eadd(&EV, EP);
1236 }
1237
1238 free_evalue_refs(&EV);
1239 free_evalue_refs(E);
1240 delete E;
1241 }
1242
1243 free_evalue_refs(&tmp);
1244
1245 out:
1246 value_clear(m);
1247 }
1248
1249 evalue* bv_ceil3(Value *coef, int len, Value d, Polyhedron *P)
1250 {
1251 Vector *val = Vector_Alloc(len);
1252
1253 Value t;
1254 value_init(t);
1255 value_set_si(t, -1);
1256 Vector_Scale(coef, val->p, t, len);
1257 value_absolute(t, d);
1258
1259 vec_ZZ num;
1260 values2zz(val->p, num, len);
1261 evalue *EP = multi_monom(num);
1262
1263 evalue tmp;
1264 value_init(tmp.d);
1265 value_init(tmp.x.n);
1266 value_set_si(tmp.x.n, 1);
1267 value_assign(tmp.d, t);
1268
1269 emul(&tmp, EP);
1270
1271 ZZ one;
1272 one = 1;
1273 ceil_mod(val->p, len, t, one, EP, P);
1274 value_clear(t);
1275
1276 /* copy EP to malloc'ed evalue */
1277 evalue *E;
1278 ALLOC(evalue, E);
1279 *E = *EP;
1280 delete EP;
1281
1282 free_evalue_refs(&tmp);
1283 Vector_Free(val);
1284
1285 return E;
1286 }
1287
1288 #ifdef USE_MODULO
1289 evalue* lattice_point(
1290 Polyhedron *i, vec_ZZ& lambda, Matrix *W, Value lcm, Polyhedron *PD)
1291 {
1292 unsigned nparam = W->NbColumns - 1;
1293
1294 Matrix* Rays = rays2(i);
1295 Matrix *T = Transpose(Rays);
1296 Matrix *T2 = Matrix_Copy(T);
1297 Matrix *inv = Matrix_Alloc(T2->NbRows, T2->NbColumns);
1298 int ok = Matrix_Inverse(T2, inv);
1299 assert(ok);
1300 Matrix_Free(Rays);
1301 Matrix_Free(T2);
1302 mat_ZZ vertex;
1303 matrix2zz(W, vertex, W->NbRows, W->NbColumns);
1304
1305 vec_ZZ num;
1306 num = lambda * vertex;
1307
1308 evalue *EP = multi_monom(num);
1309
1310 evalue tmp;
1311 value_init(tmp.d);
1312 value_init(tmp.x.n);
1313 value_set_si(tmp.x.n, 1);
1314 value_assign(tmp.d, lcm);
1315
1316 emul(&tmp, EP);
1317
1318 Matrix *L = Matrix_Alloc(inv->NbRows, W->NbColumns);
1319 Matrix_Product(inv, W, L);
1320
1321 mat_ZZ RT;
1322 matrix2zz(T, RT, T->NbRows, T->NbColumns);
1323 Matrix_Free(T);
1324
1325 vec_ZZ p = lambda * RT;
1326
1327 for (int i = 0; i < L->NbRows; ++i) {
1328 ceil_mod(L->p[i], nparam+1, lcm, p[i], EP, PD);
1329 }
1330
1331 Matrix_Free(L);
1332
1333 Matrix_Free(inv);
1334 free_evalue_refs(&tmp);
1335 return EP;
1336 }
1337 #else
1338 evalue* lattice_point(
1339 Polyhedron *i, vec_ZZ& lambda, Matrix *W, Value lcm, Polyhedron *PD)
1340 {
1341 Matrix *T = Transpose(W);
1342 unsigned nparam = T->NbRows - 1;
1343
1344 evalue *EP = new evalue();
1345 value_init(EP->d);
1346 evalue_set_si(EP, 0, 1);
1347
1348 evalue ev;
1349 Vector *val = Vector_Alloc(nparam+1);
1350 value_set_si(val->p[nparam], 1);
1351 ZZ offset(INIT_VAL, 0);
1352 value_init(ev.d);
1353 vertex_period(i, lambda, T, lcm, 0, val, EP, &ev, offset);
1354 Vector_Free(val);
1355 eadd(&ev, EP);
1356 free_evalue_refs(&ev);
1357
1358 Matrix_Free(T);
1359
1360 reduce_evalue(EP);
1361
1362 return EP;
1363 }
1364 #endif
1365
1366 void lattice_point(
1367 Param_Vertices* V, Polyhedron *i, vec_ZZ& lambda, term_info* term,
1368 Polyhedron *PD)
1369 {
1370 unsigned nparam = V->Vertex->NbColumns - 2;
1371 unsigned dim = i->Dimension;
1372 mat_ZZ vertex;
1373 vertex.SetDims(V->Vertex->NbRows, nparam+1);
1374 Value lcm, tmp;
1375 value_init(lcm);
1376 value_init(tmp);
1377 value_set_si(lcm, 1);
1378 for (int j = 0; j < V->Vertex->NbRows; ++j) {
1379 value_lcm(lcm, V->Vertex->p[j][nparam+1], &lcm);
1380 }
1381 if (value_notone_p(lcm)) {
1382 Matrix * mv = Matrix_Alloc(dim, nparam+1);
1383 for (int j = 0 ; j < dim; ++j) {
1384 value_division(tmp, lcm, V->Vertex->p[j][nparam+1]);
1385 Vector_Scale(V->Vertex->p[j], mv->p[j], tmp, nparam+1);
1386 }
1387
1388 term->E = lattice_point(i, lambda, mv, lcm, PD);
1389 term->constant = 0;
1390
1391 Matrix_Free(mv);
1392 value_clear(lcm);
1393 value_clear(tmp);
1394 return;
1395 }
1396 for (int i = 0; i < V->Vertex->NbRows; ++i) {
1397 assert(value_one_p(V->Vertex->p[i][nparam+1])); // for now
1398 values2zz(V->Vertex->p[i], vertex[i], nparam+1);
1399 }
1400
1401 vec_ZZ num;
1402 num = lambda * vertex;
1403
1404 int p = -1;
1405 int nn = 0;
1406 for (int j = 0; j < nparam; ++j)
1407 if (num[j] != 0) {
1408 ++nn;
1409 p = j;
1410 }
1411 if (nn >= 2) {
1412 term->E = multi_monom(num);
1413 term->constant = 0;
1414 } else {
1415 term->E = NULL;
1416 term->constant = num[nparam];
1417 term->pos = p;
1418 if (p != -1)
1419 term->coeff = num[p];
1420 }
1421
1422 value_clear(lcm);
1423 value_clear(tmp);
1424 }
1425
1426 static void normalize(ZZ& sign, ZZ& num, vec_ZZ& den)
1427 {
1428 unsigned dim = den.length();
1429
1430 int change = 0;
1431
1432 for (int j = 0; j < den.length(); ++j) {
1433 if (den[j] > 0)
1434 change ^= 1;
1435 else {
1436 den[j] = abs(den[j]);
1437 num += den[j];
1438 }
1439 }
1440 if (change)
1441 sign = -sign;
1442 }
1443
1444 /* input:
1445 * f: the powers in the denominator for the remaining vars
1446 * each row refers to a factor
1447 * den_s: for each factor, the power of (s+1)
1448 * sign
1449 * num_s: powers in the numerator corresponding to the summed vars
1450 * num_p: powers in the numerator corresponidng to the remaining vars
1451 * number of rays in cone: "dim" = "k"
1452 * length of each ray: "dim" = "d"
1453 * for now, it is assume: k == d
1454 * output:
1455 * den_p: for each factor
1456 * 0: independent of remaining vars
1457 * 1: power corresponds to corresponding row in f
1458 * -1: power is inverse of corresponding row in f
1459 */
1460 static void normalize(ZZ& sign,
1461 ZZ& num_s, vec_ZZ& num_p, vec_ZZ& den_s, vec_ZZ& den_p,
1462 mat_ZZ& f)
1463 {
1464 unsigned dim = f.NumRows();
1465 unsigned nparam = num_p.length();
1466 unsigned nvar = dim - nparam;
1467
1468 int change = 0;
1469
1470 for (int j = 0; j < den_s.length(); ++j) {
1471 if (den_s[j] == 0) {
1472 den_p[j] = 1;
1473 continue;
1474 }
1475 int k;
1476 for (k = 0; k < nparam; ++k)
1477 if (f[j][k] != 0)
1478 break;
1479 if (k < nparam) {
1480 if (den_s[j] > 0) {
1481 den_p[j] = -1;
1482 num_p -= f[j];
1483 } else
1484 den_p[j] = 1;
1485 } else
1486 den_p[j] = 0;
1487 if (den_s[j] > 0)
1488 change ^= 1;
1489 else {
1490 den_s[j] = abs(den_s[j]);
1491 num_s += den_s[j];
1492 }
1493 }
1494
1495 if (change)
1496 sign = -sign;
1497 }
1498
1499 struct counter : public polar_decomposer {
1500 vec_ZZ lambda;
1501 mat_ZZ rays;
1502 vec_ZZ vertex;
1503 vec_ZZ den;
1504 ZZ sign;
1505 ZZ num;
1506 int j;
1507 Polyhedron *P;
1508 unsigned dim;
1509 mpq_t count;
1510
1511 counter(Polyhedron *P) {
1512 this->P = P;
1513 dim = P->Dimension;
1514 randomvector(P, lambda, dim);
1515 rays.SetDims(dim, dim);
1516 den.SetLength(dim);
1517 mpq_init(count);
1518 }
1519
1520 void start(unsigned MaxRays);
1521
1522 ~counter() {
1523 mpq_clear(count);
1524 }
1525
1526 virtual void handle_polar(Polyhedron *P, int sign);
1527 };
1528
1529 void counter::handle_polar(Polyhedron *C, int s)
1530 {
1531 int r = 0;
1532 assert(C->NbRays-1 == dim);
1533 add_rays(rays, C, &r);
1534 for (int k = 0; k < dim; ++k) {
1535 assert(lambda * rays[k] != 0);
1536 }
1537
1538 sign = s;
1539
1540 lattice_point(P->Ray[j]+1, C, vertex);
1541 num = vertex * lambda;
1542 den = rays * lambda;
1543 normalize(sign, num, den);
1544
1545 dpoly d(dim, num);
1546 dpoly n(dim, den[0], 1);
1547 for (int k = 1; k < dim; ++k) {
1548 dpoly fact(dim, den[k], 1);
1549 n *= fact;
1550 }
1551 d.div(n, count, sign);
1552 }
1553
1554 void counter::start(unsigned MaxRays)
1555 {
1556 for (j = 0; j < P->NbRays; ++j) {
1557 Polyhedron *C = supporting_cone(P, j);
1558 decompose(C, MaxRays);
1559 }
1560 }
1561
1562 // incremental counter
1563 struct icounter : public polar_decomposer {
1564 vec_ZZ lambda;
1565 vec_ZZ vertex;
1566 //vec_ZZ den;
1567 ZZ sgn;
1568 ZZ num;
1569 ZZ one;
1570 int j;
1571 Polyhedron *P;
1572 unsigned dim;
1573 mpq_t count;
1574 mpq_t tcount;
1575 mpz_t tn;
1576 mpz_t td;
1577
1578 icounter(Polyhedron *P) {
1579 this->P = P;
1580 dim = P->Dimension;
1581 lambda.SetLength(1);
1582 lambda[0] = 1;
1583 //den.SetLength(dim);
1584 mpq_init(count);
1585 mpq_init(tcount);
1586 mpz_init(tn);
1587 mpz_init(td);
1588 one = 1;
1589 }
1590
1591 void start(unsigned MaxRays);
1592
1593 ~icounter() {
1594 mpq_clear(count);
1595 mpq_clear(tcount);
1596 mpz_clear(tn);
1597 mpz_clear(td);
1598 }
1599
1600 virtual void handle_polar(Polyhedron *P, int sign);
1601 void reduce(ZZ& c, ZZ& cd, vec_ZZ& num, vec_ZZ& den_s, mat_ZZ& den);
1602 void recurse(ZZ c, ZZ cd, vec_ZZ& num, mat_ZZ& den);
1603 };
1604
1605 void icounter::recurse(ZZ c, ZZ cd, vec_ZZ& num, mat_ZZ& den)
1606 {
1607 unsigned d = num.length();
1608 unsigned len = den.NumRows(); // number of factors in den
1609
1610 vec_ZZ den_s;
1611 den_s.SetLength(len);
1612 mat_ZZ den_r;
1613 den_r.SetDims(len, d-1);
1614
1615 int r, k;
1616 for (r = 0; r < len; ++r) {
1617 den_s[r] = den[r][0];
1618 for (k = 1; k < d; ++k)
1619 den_r[r][k-1] = den[r][k];
1620 }
1621
1622 if (d == 1) {
1623 ZZ num_s = num[0];
1624 normalize(c, num_s, den_s);
1625
1626 dpoly n(len, num_s);
1627 dpoly D(len, den_s[0], 1);
1628 for (int k = 1; k < len; ++k) {
1629 dpoly fact(len, den_s[k], 1);
1630 D *= fact;
1631 }
1632 mpq_set_si(tcount, 0, 1);
1633 n.div(D, tcount, one);
1634 zz2value(c, tn);
1635 zz2value(cd, td);
1636 mpz_mul(mpq_numref(tcount), mpq_numref(tcount), tn);
1637 mpz_mul(mpq_denref(tcount), mpq_denref(tcount), td);
1638 mpq_canonicalize(tcount);
1639 mpq_add(count, count, tcount);
1640 } else {
1641 reduce(c, cd, num, den_s, den_r);
1642 }
1643 }
1644
1645 void icounter::reduce(ZZ& c, ZZ& cd, vec_ZZ& num, vec_ZZ& den_s, mat_ZZ& den)
1646 {
1647 assert(num.length() > 1);
1648 unsigned d = num.length()-1;
1649 unsigned len = den_s.length(); // number of factors in den
1650 ZZ num_s = num[0];
1651 int k = 0;
1652 vec_ZZ num_p;
1653 num_p.SetLength(d);
1654 for (k = 1 ; k <= d; ++k)
1655 num_p[k-1] = num[k];
1656
1657 vec_ZZ den_p;
1658 den_p.SetLength(len);
1659
1660 normalize(c, num_s, num_p, den_s, den_p, den);
1661
1662 /* Since we're working incrementally, we should look
1663 * for the "easiest" parameter first.
1664 * In particular we should handle the parameters such
1665 * that no_param + only_param == len, since that allows
1666 * us to decouple the problem and the split off part
1667 * may very well be zero
1668 */
1669 int only_param = 0;
1670 int no_param = 0;
1671 for (int k = 0; k < len; ++k) {
1672 if (den_p[k] == 0)
1673 ++no_param;
1674 else if (den_s[k] == 0)
1675 ++only_param;
1676 }
1677 if (no_param == 0) {
1678 for (int k = 0; k < len; ++k)
1679 if (den_p[k] == -1)
1680 den[k] = -den[k];
1681 recurse(c, cd, num_p, den);
1682 } else {
1683 int k, l;
1684 mat_ZZ pden;
1685 pden.SetDims(only_param, d);
1686
1687 if (no_param + only_param == len) {
1688 for (k = 0, l = 0; k < len; ++k)
1689 if (den_p[k] != 0)
1690 pden[l++] = den[k];
1691
1692 for (k = 0; k < len; ++k)
1693 if (den_s[k] != 0)
1694 break;
1695
1696 dpoly n(no_param, num_s);
1697 dpoly D(no_param, den_s[k], 1);
1698 for ( ; ++k < len; )
1699 if (den_s[k] != 0) {
1700 dpoly fact(no_param, den_s[k], 1);
1701 D *= fact;
1702 }
1703
1704 mpq_set_si(tcount, 0, 1);
1705 n.div(D, tcount, one);
1706
1707 ZZ qn, qd;
1708 value2zz(mpq_numref(tcount), qn);
1709 value2zz(mpq_denref(tcount), qd);
1710
1711 qn *= c;
1712 qd *= cd;
1713
1714 if (qn != 0)
1715 recurse(qn, qd, num_p, pden);
1716 } else {
1717 dpoly_r * r = 0;
1718
1719 for (k = 0, l = 0; k < len; ++k)
1720 if (den_s[k] == 0)
1721 pden[l++] = den[k];
1722
1723 for (k = 0; k < len; ++k)
1724 if (den_p[k] == 0)
1725 break;
1726
1727 dpoly n(no_param, num_s);
1728 dpoly D(no_param, den_s[k], 1);
1729 for ( ; ++k < len; )
1730 if (den_p[k] == 0) {
1731 dpoly fact(no_param, den_s[k], 1);
1732 D *= fact;
1733 }
1734
1735 for (k = 0; k < len; ++k) {
1736 if (den_s[k] == 0 || den_p[k] == 0)
1737 continue;
1738
1739 dpoly pd(no_param-1, den_s[k], 1);
1740 int s = den_p[k] < 0 ? -1 : 1;
1741
1742 if (r == 0)
1743 r = new dpoly_r(n, pd, k, s, len);
1744 else {
1745 dpoly_r *nr = new dpoly_r(r, pd, k, s, len);
1746 delete r;
1747 r = nr;
1748 }
1749 }
1750
1751 dpoly_r *rc = r->div(D);
1752
1753 rc->denom *= cd;
1754
1755 int common = pden.NumRows();
1756 vector< dpoly_r_term * >& final = rc->c[rc->len-1];
1757 int rows;
1758 for (int j = 0; j < final.size(); ++j) {
1759 if (final[j]->coeff == 0)
1760 continue;
1761 rows = common;
1762 pden.SetDims(rows, pden.NumCols());
1763 for (int k = 0; k < rc->dim; ++k) {
1764 int n = final[j]->powers[k];
1765 if (n == 0)
1766 continue;
1767 int abs_n = n < 0 ? -n : n;
1768 pden.SetDims(rows+abs_n, pden.NumCols());
1769 for (int l = 0; l < abs_n; ++l) {
1770 if (n > 0)
1771 pden[rows+l] = den[k];
1772 else
1773 pden[rows+l] = -den[k];
1774 }
1775 rows += abs_n;
1776 }
1777 final[j]->coeff *= c;
1778 recurse(final[j]->coeff, rc->denom, num_p, pden);
1779 }
1780
1781 delete rc;
1782 delete r;
1783 }
1784 }
1785 }
1786
1787 void icounter::handle_polar(Polyhedron *C, int s)
1788 {
1789 assert(C->NbRays-1 == dim);
1790
1791 sgn = s;
1792
1793 lattice_point(P->Ray[j]+1, C, vertex);
1794
1795 vec_ZZ den_s;
1796 den_s.SetLength(dim);
1797 mat_ZZ den;
1798 den.SetDims(dim, dim-1);
1799
1800 int r;
1801 for (r = 0; r < dim; ++r) {
1802 value2zz(C->Ray[r][1], den_s[r]);
1803 values2zz(C->Ray[r]+1+1, den[r], dim-1);
1804 }
1805
1806 reduce(sgn, one, vertex, den_s, den);
1807 }
1808
1809 void icounter::start(unsigned MaxRays)
1810 {
1811 for (j = 0; j < P->NbRays; ++j) {
1812 Polyhedron *C = supporting_cone(P, j);
1813 decompose(C, MaxRays);
1814 }
1815 }
1816
1817 typedef Polyhedron * Polyhedron_p;
1818
1819 void barvinok_count(Polyhedron *P, Value* result, unsigned NbMaxCons)
1820 {
1821 Polyhedron ** vcone;
1822 ZZ sign;
1823 unsigned dim;
1824 int allocated = 0;
1825 Value factor;
1826 Polyhedron *Q;
1827 int r = 0;
1828
1829 if (emptyQ(P)) {
1830 value_set_si(*result, 0);
1831 return;
1832 }
1833 if (P->NbBid == 0)
1834 for (; r < P->NbRays; ++r)
1835 if (value_zero_p(P->Ray[r][P->Dimension+1]))
1836 break;
1837 if (P->NbBid !=0 || r < P->NbRays) {
1838 value_set_si(*result, -1);
1839 return;
1840 }
1841 if (P->NbEq != 0) {
1842 P = remove_equalities(P);
1843 if (emptyQ(P)) {
1844 Polyhedron_Free(P);
1845 value_set_si(*result, 0);
1846 return;
1847 }
1848 allocated = 1;
1849 }
1850 value_init(factor);
1851 value_set_si(factor, 1);
1852 Q = Polyhedron_Reduce(P, &factor);
1853 if (Q) {
1854 if (allocated)
1855 Polyhedron_Free(P);
1856 P = Q;
1857 allocated = 1;
1858 }
1859 if (P->Dimension == 0) {
1860 value_assign(*result, factor);
1861 if (allocated)
1862 Polyhedron_Free(P);
1863 value_clear(factor);
1864 return;
1865 }
1866
1867 icounter cnt(P);
1868 cnt.start(NbMaxCons);
1869
1870 assert(value_one_p(&cnt.count[0]._mp_den));
1871 value_multiply(*result, &cnt.count[0]._mp_num, factor);
1872
1873 if (allocated)
1874 Polyhedron_Free(P);
1875 value_clear(factor);
1876 }
1877
1878 static void uni_polynom(int param, Vector *c, evalue *EP)
1879 {
1880 unsigned dim = c->Size-2;
1881 value_init(EP->d);
1882 value_set_si(EP->d,0);
1883 EP->x.p = new_enode(polynomial, dim+1, param+1);
1884 for (int j = 0; j <= dim; ++j)
1885 evalue_set(&EP->x.p->arr[j], c->p[j], c->p[dim+1]);
1886 }
1887
1888 static void multi_polynom(Vector *c, evalue* X, evalue *EP)
1889 {
1890 unsigned dim = c->Size-2;
1891 evalue EC;
1892
1893 value_init(EC.d);
1894 evalue_set(&EC, c->p[dim], c->p[dim+1]);
1895
1896 value_init(EP->d);
1897 evalue_set(EP, c->p[dim], c->p[dim+1]);
1898
1899 for (int i = dim-1; i >= 0; --i) {
1900 emul(X, EP);
1901 value_assign(EC.x.n, c->p[i]);
1902 eadd(&EC, EP);
1903 }
1904 free_evalue_refs(&EC);
1905 }
1906
1907 Polyhedron *unfringe (Polyhedron *P, unsigned MaxRays)
1908 {
1909 int len = P->Dimension+2;
1910 Polyhedron *T, *R = P;
1911 Value g;
1912 value_init(g);
1913 Vector *row = Vector_Alloc(len);
1914 value_set_si(row->p[0], 1);
1915
1916 R = DomainConstraintSimplify(Polyhedron_Copy(P), MaxRays);
1917
1918 Matrix *M = Matrix_Alloc(2, len-1);
1919 value_set_si(M->p[1][len-2], 1);
1920 for (int v = 0; v < P->Dimension; ++v) {
1921 value_set_si(M->p[0][v], 1);
1922 Polyhedron *I = Polyhedron_Image(P, M, 2+1);
1923 value_set_si(M->p[0][v], 0);
1924 for (int r = 0; r < I->NbConstraints; ++r) {
1925 if (value_zero_p(I->Constraint[r][0]))
1926 continue;
1927 if (value_zero_p(I->Constraint[r][1]))
1928 continue;
1929 if (value_one_p(I->Constraint[r][1]))
1930 continue;
1931 if (value_mone_p(I->Constraint[r][1]))
1932 continue;
1933 value_absolute(g, I->Constraint[r][1]);
1934 Vector_Set(row->p+1, 0, len-2);
1935 value_division(row->p[1+v], I->Constraint[r][1], g);
1936 mpz_fdiv_q(row->p[len-1], I->Constraint[r][2], g);
1937 T = R;
1938 R = AddConstraints(row->p, 1, R, MaxRays);
1939 if (T != P)
1940 Polyhedron_Free(T);
1941 }
1942 Polyhedron_Free(I);
1943 }
1944 Matrix_Free(M);
1945 Vector_Free(row);
1946 value_clear(g);
1947 return R;
1948 }
1949
1950 static Polyhedron *reduce_domain(Polyhedron *D, Matrix *CT, Polyhedron *CEq,
1951 Polyhedron **fVD, int nd, unsigned MaxRays)
1952 {
1953 assert(CEq);
1954
1955 Polyhedron *Dt;
1956 Dt = CT ? DomainPreimage(D, CT, MaxRays) : D;
1957 Polyhedron *rVD = DomainIntersection(Dt, CEq, MaxRays);
1958
1959 /* if rVD is empty or too small in geometric dimension */
1960 if(!rVD || emptyQ(rVD) ||
1961 (rVD->Dimension-rVD->NbEq < Dt->Dimension-Dt->NbEq-CEq->NbEq)) {
1962 if(rVD)
1963 Domain_Free(rVD);
1964 if (CT)
1965 Domain_Free(Dt);
1966 return 0; /* empty validity domain */
1967 }
1968
1969 if (CT)
1970 Domain_Free(Dt);
1971
1972 fVD[nd] = Domain_Copy(rVD);
1973 for (int i = 0 ; i < nd; ++i) {
1974 Polyhedron *I = DomainIntersection(fVD[nd], fVD[i], MaxRays);
1975 if (emptyQ(I)) {
1976 Domain_Free(I);
1977 continue;
1978 }
1979 Polyhedron *F = DomainSimplify(I, fVD[nd], MaxRays);
1980 if (F->NbEq == 1) {
1981 Polyhedron *T = rVD;
1982 rVD = DomainDifference(rVD, F, MaxRays);
1983 Domain_Free(T);
1984 }
1985 Domain_Free(F);
1986 Domain_Free(I);
1987 }
1988
1989 rVD = DomainConstraintSimplify(rVD, MaxRays);
1990 if (emptyQ(rVD)) {
1991 Domain_Free(fVD[nd]);
1992 Domain_Free(rVD);
1993 return 0;
1994 }
1995
1996 Value c;
1997 value_init(c);
1998 barvinok_count(rVD, &c, MaxRays);
1999 if (value_zero_p(c)) {
2000 Domain_Free(rVD);
2001 rVD = 0;
2002 }
2003 value_clear(c);
2004
2005 return rVD;
2006 }
2007
2008 static bool Polyhedron_is_infinite(Polyhedron *P, unsigned nparam)
2009 {
2010 int r;
2011 for (r = 0; r < P->NbRays; ++r)
2012 if (value_zero_p(P->Ray[r][0]) ||
2013 value_zero_p(P->Ray[r][P->Dimension+1])) {
2014 int i;
2015 for (i = P->Dimension - nparam; i < P->Dimension; ++i)
2016 if (value_notzero_p(P->Ray[r][i+1]))
2017 break;
2018 if (i >= P->Dimension)
2019 break;
2020 }
2021 return r < P->NbRays;
2022 }
2023
2024 /* Check whether all rays point in the positive directions
2025 * for the parameters
2026 */
2027 static bool Polyhedron_has_positive_rays(Polyhedron *P, unsigned nparam)
2028 {
2029 int r;
2030 for (r = 0; r < P->NbRays; ++r)
2031 if (value_zero_p(P->Ray[r][P->Dimension+1])) {
2032 int i;
2033 for (i = P->Dimension - nparam; i < P->Dimension; ++i)
2034 if (value_neg_p(P->Ray[r][i+1]))
2035 return false;
2036 }
2037 return true;
2038 }
2039
2040 typedef evalue * evalue_p;
2041
2042 struct enumerator : public polar_decomposer {
2043 vec_ZZ lambda;
2044 unsigned dim, nbV;
2045 evalue ** vE;
2046 int _i;
2047 mat_ZZ rays;
2048 vec_ZZ den;
2049 ZZ sign;
2050 Polyhedron *P;
2051 Param_Vertices *V;
2052 term_info num;
2053 Vector *c;
2054 mpq_t count;
2055
2056 enumerator(Polyhedron *P, unsigned dim, unsigned nbV) {
2057 this->P = P;
2058 this->dim = dim;
2059 this->nbV = nbV;
2060 randomvector(P, lambda, dim);
2061 rays.SetDims(dim, dim);
2062 den.SetLength(dim);
2063 c = Vector_Alloc(dim+2);
2064
2065 vE = new evalue_p[nbV];
2066 for (int j = 0; j < nbV; ++j)
2067 vE[j] = 0;
2068
2069 mpq_init(count);
2070 }
2071
2072 void decompose_at(Param_Vertices *V, int _i, unsigned MaxRays) {
2073 Polyhedron *C = supporting_cone_p(P, V);
2074 this->_i = _i;
2075 this->V = V;
2076
2077 vE[_i] = new evalue;
2078 value_init(vE[_i]->d);
2079 evalue_set_si(vE[_i], 0, 1);
2080
2081 decompose(C, MaxRays);
2082 }
2083
2084 ~enumerator() {
2085 mpq_clear(count);
2086 Vector_Free(c);
2087
2088 for (int j = 0; j < nbV; ++j)
2089 if (vE[j]) {
2090 free_evalue_refs(vE[j]);
2091 delete vE[j];
2092 }
2093 delete [] vE;
2094 }
2095
2096 virtual void handle_polar(Polyhedron *P, int sign);
2097 };
2098
2099 void enumerator::handle_polar(Polyhedron *C, int s)
2100 {
2101 int r = 0;
2102 assert(C->NbRays-1 == dim);
2103 add_rays(rays, C, &r);
2104 for (int k = 0; k < dim; ++k) {
2105 assert(lambda * rays[k] != 0);
2106 }
2107
2108 sign = s;
2109
2110 lattice_point(V, C, lambda, &num, 0);
2111 den = rays * lambda;
2112 normalize(sign, num.constant, den);
2113
2114 dpoly n(dim, den[0], 1);
2115 for (int k = 1; k < dim; ++k) {
2116 dpoly fact(dim, den[k], 1);
2117 n *= fact;
2118 }
2119 if (num.E != NULL) {
2120 ZZ one(INIT_VAL, 1);
2121 dpoly_n d(dim, num.constant, one);
2122 d.div(n, c, sign);
2123 evalue EV;
2124 multi_polynom(c, num.E, &EV);
2125 eadd(&EV , vE[_i]);
2126 free_evalue_refs(&EV);
2127 free_evalue_refs(num.E);
2128 delete num.E;
2129 } else if (num.pos != -1) {
2130 dpoly_n d(dim, num.constant, num.coeff);
2131 d.div(n, c, sign);
2132 evalue EV;
2133 uni_polynom(num.pos, c, &EV);
2134 eadd(&EV , vE[_i]);
2135 free_evalue_refs(&EV);
2136 } else {
2137 mpq_set_si(count, 0, 1);
2138 dpoly d(dim, num.constant);
2139 d.div(n, count, sign);
2140 evalue EV;
2141 value_init(EV.d);
2142 evalue_set(&EV, &count[0]._mp_num, &count[0]._mp_den);
2143 eadd(&EV , vE[_i]);
2144 free_evalue_refs(&EV);
2145 }
2146 }
2147
2148 evalue* barvinok_enumerate_ev(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
2149 {
2150 //P = unfringe(P, MaxRays);
2151 Polyhedron *CEq = NULL, *rVD, *pVD, *CA;
2152 Matrix *CT = NULL;
2153 Param_Polyhedron *PP = NULL;
2154 Param_Domain *D, *next;
2155 Param_Vertices *V;
2156 int r = 0;
2157 unsigned nparam = C->Dimension;
2158 evalue *eres;
2159 ALLOC(evalue, eres);
2160 value_init(eres->d);
2161 value_set_si(eres->d, 0);
2162
2163 evalue factor;
2164 value_init(factor.d);
2165 evalue_set_si(&factor, 1, 1);
2166
2167 CA = align_context(C, P->Dimension, MaxRays);
2168 P = DomainIntersection(P, CA, MaxRays);
2169 Polyhedron_Free(CA);
2170
2171 if (C->Dimension == 0 || emptyQ(P)) {
2172 constant:
2173 eres->x.p = new_enode(partition, 2, C->Dimension);
2174 EVALUE_SET_DOMAIN(eres->x.p->arr[0],
2175 DomainConstraintSimplify(CEq ? CEq : Polyhedron_Copy(C), MaxRays));
2176 value_set_si(eres->x.p->arr[1].d, 1);
2177 value_init(eres->x.p->arr[1].x.n);
2178 if (emptyQ(P))
2179 value_set_si(eres->x.p->arr[1].x.n, 0);
2180 else
2181 barvinok_count(P, &eres->x.p->arr[1].x.n, MaxRays);
2182 out:
2183 emul(&factor, eres);
2184 reduce_evalue(eres);
2185 free_evalue_refs(&factor);
2186 Polyhedron_Free(P);
2187 if (CT)
2188 Matrix_Free(CT);
2189 if (PP)
2190 Param_Polyhedron_Free(PP);
2191
2192 return eres;
2193 }
2194 if (Polyhedron_is_infinite(P, nparam))
2195 goto constant;
2196
2197 if (P->NbEq != 0) {
2198 Matrix *f;
2199 P = remove_equalities_p(P, P->Dimension-nparam, &f);
2200 mask(f, &factor);
2201 Matrix_Free(f);
2202 }
2203 if (P->Dimension == nparam) {
2204 CEq = P;
2205 P = Universe_Polyhedron(0);
2206 goto constant;
2207 }
2208
2209 Polyhedron *Q = ParamPolyhedron_Reduce(P, P->Dimension-nparam, &factor);
2210 if (Q) {
2211 Polyhedron_Free(P);
2212 if (Q->Dimension == nparam) {
2213 CEq = Q;
2214 P = Universe_Polyhedron(0);
2215 goto constant;
2216 }
2217 P = Q;
2218 }
2219 Polyhedron *oldP = P;
2220 PP = Polyhedron2Param_SimplifiedDomain(&P,C,MaxRays,&CEq,&CT);
2221 if (P != oldP)
2222 Polyhedron_Free(oldP);
2223
2224 if (isIdentity(CT)) {
2225 Matrix_Free(CT);
2226 CT = NULL;
2227 } else {
2228 assert(CT->NbRows != CT->NbColumns);
2229 if (CT->NbRows == 1) // no more parameters
2230 goto constant;
2231 nparam = CT->NbRows - 1;
2232 }
2233
2234 unsigned dim = P->Dimension - nparam;
2235
2236 enumerator et(P, dim, PP->nbV);
2237
2238 int nd;
2239 for (nd = 0, D=PP->D; D; ++nd, D=D->next);
2240 struct section { Polyhedron *D; evalue E; };
2241 section *s = new section[nd];
2242 Polyhedron **fVD = new Polyhedron_p[nd];
2243
2244 for(nd = 0, D=PP->D; D; D=next) {
2245 next = D->next;
2246
2247 Polyhedron *rVD = reduce_domain(D->Domain, CT, CEq,
2248 fVD, nd, MaxRays);
2249 if (!rVD)
2250 continue;
2251
2252 pVD = CT ? DomainImage(rVD,CT,MaxRays) : rVD;
2253
2254 value_init(s[nd].E.d);
2255 evalue_set_si(&s[nd].E, 0, 1);
2256
2257 FORALL_PVertex_in_ParamPolyhedron(V,D,PP) // _i is internal counter
2258 if (!et.vE[_i])
2259 et.decompose_at(V, _i, MaxRays);
2260 eadd(et.vE[_i] , &s[nd].E);
2261 END_FORALL_PVertex_in_ParamPolyhedron;
2262 reduce_in_domain(&s[nd].E, pVD);
2263
2264 if (CT)
2265 addeliminatedparams_evalue(&s[nd].E, CT);
2266 s[nd].D = rVD;
2267 ++nd;
2268 if (rVD != pVD)
2269 Domain_Free(pVD);
2270 }
2271
2272 if (nd == 0)
2273 evalue_set_si(eres, 0, 1);
2274 else {
2275 eres->x.p = new_enode(partition, 2*nd, C->Dimension);
2276 for (int j = 0; j < nd; ++j) {
2277 EVALUE_SET_DOMAIN(eres->x.p->arr[2*j], s[j].D);
2278 value_clear(eres->x.p->arr[2*j+1].d);
2279 eres->x.p->arr[2*j+1] = s[j].E;
2280 Domain_Free(fVD[j]);
2281 }
2282 }
2283 delete [] s;
2284 delete [] fVD;
2285
2286
2287 if (CEq)
2288 Polyhedron_Free(CEq);
2289
2290 goto out;
2291 }
2292
2293 Enumeration* barvinok_enumerate(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
2294 {
2295 evalue *EP = barvinok_enumerate_ev(P, C, MaxRays);
2296
2297 return partition2enumeration(EP);
2298 }
2299
2300 static void SwapColumns(Value **V, int n, int i, int j)
2301 {
2302 for (int r = 0; r < n; ++r)
2303 value_swap(V[r][i], V[r][j]);
2304 }
2305
2306 static void SwapColumns(Polyhedron *P, int i, int j)
2307 {
2308 SwapColumns(P->Constraint, P->NbConstraints, i, j);
2309 SwapColumns(P->Ray, P->NbRays, i, j);
2310 }
2311
2312 static void negative_test_constraint(Value *l, Value *u, Value *c, int pos,
2313 int len, Value *v)
2314 {
2315 value_oppose(*v, u[pos+1]);
2316 Vector_Combine(l+1, u+1, c+1, *v, l[pos+1], len-1);
2317 value_multiply(*v, *v, l[pos+1]);
2318 value_substract(c[len-1], c[len-1], *v);
2319 value_set_si(*v, -1);
2320 Vector_Scale(c+1, c+1, *v, len-1);
2321 value_decrement(c[len-1], c[len-1]);
2322 ConstraintSimplify(c, c, len, v);
2323 }
2324
2325 static bool parallel_constraints(Value *l, Value *u, Value *c, int pos,
2326 int len)
2327 {
2328 bool parallel;
2329 Value g1;
2330 Value g2;
2331 value_init(g1);
2332 value_init(g2);
2333
2334 Vector_Gcd(&l[1+pos], len, &g1);
2335 Vector_Gcd(&u[1+pos], len, &g2);
2336 Vector_Combine(l+1+pos, u+1+pos, c+1, g2, g1, len);
2337 parallel = First_Non_Zero(c+1, len) == -1;
2338
2339 value_clear(g1);
2340 value_clear(g2);
2341
2342 return parallel;
2343 }
2344
2345 static void negative_test_constraint7(Value *l, Value *u, Value *c, int pos,
2346 int exist, int len, Value *v)
2347 {
2348 Value g;
2349 value_init(g);
2350
2351 Vector_Gcd(&u[1+pos], exist, v);
2352 Vector_Gcd(&l[1+pos], exist, &g);
2353 Vector_Combine(l+1, u+1, c+1, *v, g, len-1);
2354 value_multiply(*v, *v, g);
2355 value_substract(c[len-1], c[len-1], *v);
2356 value_set_si(*v, -1);
2357 Vector_Scale(c+1, c+1, *v, len-1);
2358 value_decrement(c[len-1], c[len-1]);
2359 ConstraintSimplify(c, c, len, v);
2360
2361 value_clear(g);
2362 }
2363
2364 static void oppose_constraint(Value *c, int len, Value *v)
2365 {
2366 value_set_si(*v, -1);
2367 Vector_Scale(c+1, c+1, *v, len-1);
2368 value_decrement(c[len-1], c[len-1]);
2369 }
2370
2371 static bool SplitOnConstraint(Polyhedron *P, int i, int l, int u,
2372 int nvar, int len, int exist, int MaxRays,
2373 Vector *row, Value& f, bool independent,
2374 Polyhedron **pos, Polyhedron **neg)
2375 {
2376 negative_test_constraint(P->Constraint[l], P->Constraint[u],
2377 row->p, nvar+i, len, &f);
2378 *neg = AddConstraints(row->p, 1, P, MaxRays);
2379
2380 /* We found an independent, but useless constraint
2381 * Maybe we should detect this earlier and not
2382 * mark the variable as INDEPENDENT
2383 */
2384 if (emptyQ((*neg))) {
2385 Polyhedron_Free(*neg);
2386 return false;
2387 }
2388
2389 oppose_constraint(row->p, len, &f);
2390 *pos = AddConstraints(row->p, 1, P, MaxRays);
2391
2392 if (emptyQ((*pos))) {
2393 Polyhedron_Free(*neg);
2394 Polyhedron_Free(*pos);
2395 return false;
2396 }
2397
2398 return true;
2399 }
2400
2401 /*
2402 * unimodularly transform P such that constraint r is transformed
2403 * into a constraint that involves only a single (the first)
2404 * existential variable
2405 *
2406 */
2407 static Polyhedron *rotate_along(Polyhedron *P, int r, int nvar, int exist,
2408 unsigned MaxRays)
2409 {
2410 Value g;
2411 value_init(g);
2412
2413 Vector *row = Vector_Alloc(exist);
2414 Vector_Copy(P->Constraint[r]+1+nvar, row->p, exist);
2415 Vector_Gcd(row->p, exist, &g);
2416 if (value_notone_p(g))
2417 Vector_AntiScale(row->p, row->p, g, exist);
2418 value_clear(g);
2419
2420 Matrix *M = unimodular_complete(row);
2421 Matrix *M2 = Matrix_Alloc(P->Dimension+1, P->Dimension+1);
2422 for (r = 0; r < nvar; ++r)
2423 value_set_si(M2->p[r][r], 1);
2424 for ( ; r < nvar+exist; ++r)
2425 Vector_Copy(M->p[r-nvar], M2->p[r]+nvar, exist);
2426 for ( ; r < P->Dimension+1; ++r)
2427 value_set_si(M2->p[r][r], 1);
2428 Polyhedron *T = Polyhedron_Image(P, M2, MaxRays);
2429
2430 Matrix_Free(M2);
2431 Matrix_Free(M);
2432 Vector_Free(row);
2433
2434 return T;
2435 }
2436
2437 static bool SplitOnVar(Polyhedron *P, int i,
2438 int nvar, int len, int exist, int MaxRays,
2439 Vector *row, Value& f, bool independent,
2440 Polyhedron **pos, Polyhedron **neg)
2441 {
2442 int j;
2443
2444 for (int l = P->NbEq; l < P->NbConstraints; ++l) {
2445 if (value_negz_p(P->Constraint[l][nvar+i+1]))
2446 continue;
2447
2448 if (independent) {
2449 for (j = 0; j < exist; ++j)
2450 if (j != i && value_notzero_p(P->Constraint[l][nvar+j+1]))
2451 break;
2452 if (j < exist)
2453 continue;
2454 }
2455
2456 for (int u = P->NbEq; u < P->NbConstraints; ++u) {
2457 if (value_posz_p(P->Constraint[u][nvar+i+1]))
2458 continue;
2459
2460 if (independent) {
2461 for (j = 0; j < exist; ++j)
2462 if (j != i && value_notzero_p(P->Constraint[u][nvar+j+1]))
2463 break;
2464 if (j < exist)
2465 continue;
2466 }
2467
2468 if (SplitOnConstraint(P, i, l, u,
2469 nvar, len, exist, MaxRays,
2470 row, f, independent,
2471 pos, neg)) {
2472 if (independent) {
2473 if (i != 0)
2474 SwapColumns(*neg, nvar+1, nvar+1+i);
2475 }
2476 return true;
2477 }
2478 }
2479 }
2480
2481 return false;
2482 }
2483
2484 static bool double_bound_pair(Polyhedron *P, int nvar, int exist,
2485 int i, int l1, int l2,
2486 Polyhedron **pos, Polyhedron **neg)
2487 {
2488 Value f;
2489 value_init(f);
2490 Vector *row = Vector_Alloc(P->Dimension+2);
2491 value_set_si(row->p[0], 1);
2492 value_oppose(f, P->Constraint[l1][nvar+i+1]);
2493 Vector_Combine(P->Constraint[l1]+1, P->Constraint[l2]+1,
2494 row->p+1,
2495 P->Constraint[l2][nvar+i+1], f,
2496 P->Dimension+1);
2497 ConstraintSimplify(row->p, row->p, P->Dimension+2, &f);
2498 *pos = AddConstraints(row->p, 1, P, 0);
2499 value_set_si(f, -1);
2500 Vector_Scale(row->p+1, row->p+1, f, P->Dimension+1);
2501 value_decrement(row->p[P->Dimension+1], row->p[P->Dimension+1]);
2502 *neg = AddConstraints(row->p, 1, P, 0);
2503 Vector_Free(row);
2504 value_clear(f);
2505
2506 return !emptyQ((*pos)) && !emptyQ((*neg));
2507 }
2508
2509 static bool double_bound(Polyhedron *P, int nvar, int exist,
2510 Polyhedron **pos, Polyhedron **neg)
2511 {
2512 for (int i = 0; i < exist; ++i) {
2513 int l1, l2;
2514 for (l1 = P->NbEq; l1 < P->NbConstraints; ++l1) {
2515 if (value_negz_p(P->Constraint[l1][nvar+i+1]))
2516 continue;
2517 for (l2 = l1 + 1; l2 < P->NbConstraints; ++l2) {
2518 if (value_negz_p(P->Constraint[l2][nvar+i+1]))
2519 continue;
2520 if (double_bound_pair(P, nvar, exist, i, l1, l2, pos, neg))
2521 return true;
2522 }
2523 }
2524 for (l1 = P->NbEq; l1 < P->NbConstraints; ++l1) {
2525 if (value_posz_p(P->Constraint[l1][nvar+i+1]))
2526 continue;
2527 if (l1 < P->NbConstraints)
2528 for (l2 = l1 + 1; l2 < P->NbConstraints; ++l2) {
2529 if (value_posz_p(P->Constraint[l2][nvar+i+1]))
2530 continue;
2531 if (double_bound_pair(P, nvar, exist, i, l1, l2, pos, neg))
2532 return true;
2533 }
2534 }
2535 return false;
2536 }
2537 return false;
2538 }
2539
2540 enum constraint {
2541 ALL_POS = 1 << 0,
2542 ONE_NEG = 1 << 1,
2543 INDEPENDENT = 1 << 2,
2544 ROT_NEG = 1 << 3
2545 };
2546
2547 static evalue* enumerate_or(Polyhedron *D,
2548 unsigned exist, unsigned nparam, unsigned MaxRays)
2549 {
2550 #ifdef DEBUG_ER
2551 fprintf(stderr, "\nER: Or\n");
2552 #endif /* DEBUG_ER */
2553
2554 Polyhedron *N = D->next;
2555 D->next = 0;
2556 evalue *EP =
2557 barvinok_enumerate_e(D, exist, nparam, MaxRays);
2558 Polyhedron_Free(D);
2559
2560 for (D = N; D; D = N) {
2561 N = D->next;
2562 D->next = 0;
2563
2564 evalue *EN =
2565 barvinok_enumerate_e(D, exist, nparam, MaxRays);
2566
2567 eor(EN, EP);
2568 free_evalue_refs(EN);
2569 free(EN);
2570 Polyhedron_Free(D);
2571 }
2572
2573 reduce_evalue(EP);
2574
2575 return EP;
2576 }
2577
2578 static evalue* enumerate_sum(Polyhedron *P,
2579 unsigned exist, unsigned nparam, unsigned MaxRays)
2580 {
2581 int nvar = P->Dimension - exist - nparam;
2582 int toswap = nvar < exist ? nvar : exist;
2583 for (int i = 0; i < toswap; ++i)
2584 SwapColumns(P, 1 + i, nvar+exist - i);
2585 nparam += nvar;
2586
2587 #ifdef DEBUG_ER
2588 fprintf(stderr, "\nER: Sum\n");
2589 #endif /* DEBUG_ER */
2590
2591 evalue *EP = barvinok_enumerate_e(P, exist, nparam, MaxRays);
2592
2593 for (int i = 0; i < /* nvar */ nparam; ++i) {
2594 Matrix *C = Matrix_Alloc(1, 1 + nparam + 1);
2595 value_set_si(C->p[0][0], 1);
2596 evalue split;
2597 value_init(split.d);
2598 value_set_si(split.d, 0);
2599 split.x.p = new_enode(partition, 4, nparam);
2600 value_set_si(C->p[0][1+i], 1);
2601 Matrix *C2 = Matrix_Copy(C);
2602 EVALUE_SET_DOMAIN(split.x.p->arr[0],
2603 Constraints2Polyhedron(C2, MaxRays));
2604 Matrix_Free(C2);
2605 evalue_set_si(&split.x.p->arr[1], 1, 1);
2606 value_set_si(C->p[0][1+i], -1);
2607 value_set_si(C->p[0][1+nparam], -1);
2608 EVALUE_SET_DOMAIN(split.x.p->arr[2],
2609 Constraints2Polyhedron(C, MaxRays));
2610 evalue_set_si(&split.x.p->arr[3], 1, 1);
2611 emul(&split, EP);
2612 free_evalue_refs(&split);
2613 Matrix_Free(C);
2614 }
2615 reduce_evalue(EP);
2616 evalue_range_reduction(EP);
2617
2618 evalue_frac2floor(EP);
2619
2620 evalue *sum = esum(EP, nvar);
2621
2622 free_evalue_refs(EP);
2623 free(EP);
2624 EP = sum;
2625
2626 evalue_range_reduction(EP);
2627
2628 return EP;
2629 }
2630
2631 static evalue* split_sure(Polyhedron *P, Polyhedron *S,
2632 unsigned exist, unsigned nparam, unsigned MaxRays)
2633 {
2634 int nvar = P->Dimension - exist - nparam;
2635
2636 Matrix *M = Matrix_Alloc(exist, S->Dimension+2);
2637 for (int i = 0; i < exist; ++i)
2638 value_set_si(M->p[i][nvar+i+1], 1);
2639 Polyhedron *O = S;
2640 S = DomainAddRays(S, M, MaxRays);
2641 Polyhedron_Free(O);
2642 Polyhedron *F = DomainAddRays(P, M, MaxRays);
2643 Polyhedron *D = DomainDifference(F, S, MaxRays);
2644 O = D;
2645 D = Disjoint_Domain(D, 0, MaxRays);
2646 Polyhedron_Free(F);
2647 Domain_Free(O);
2648 Matrix_Free(M);
2649
2650 M = Matrix_Alloc(P->Dimension+1-exist, P->Dimension+1);
2651 for (int j = 0; j < nvar; ++j)
2652 value_set_si(M->p[j][j], 1);
2653 for (int j = 0; j < nparam+1; ++j)
2654 value_set_si(M->p[nvar+j][nvar+exist+j], 1);
2655 Polyhedron *T = Polyhedron_Image(S, M, MaxRays);
2656 evalue *EP = barvinok_enumerate_e(T, 0, nparam, MaxRays);
2657 Polyhedron_Free(S);
2658 Polyhedron_Free(T);
2659 Matrix_Free(M);
2660
2661 for (Polyhedron *Q = D; Q; Q = Q->next) {
2662 Polyhedron *N = Q->next;
2663 Q->next = 0;
2664 T = DomainIntersection(P, Q, MaxRays);
2665 evalue *E = barvinok_enumerate_e(T, exist, nparam, MaxRays);
2666 eadd(E, EP);
2667 free_evalue_refs(E);
2668 free(E);
2669 Polyhedron_Free(T);
2670 Q->next = N;
2671 }
2672 Domain_Free(D);
2673 return EP;
2674 }
2675
2676 static evalue* enumerate_sure(Polyhedron *P,
2677 unsigned exist, unsigned nparam, unsigned MaxRays)
2678 {
2679 int i;
2680 Polyhedron *S = P;
2681 int nvar = P->Dimension - exist - nparam;
2682 Value lcm;
2683 Value f;
2684 value_init(lcm);
2685 value_init(f);
2686
2687 for (i = 0; i < exist; ++i) {
2688 Matrix *M = Matrix_Alloc(S->NbConstraints, S->Dimension+2);
2689 int c = 0;
2690 value_set_si(lcm, 1);
2691 for (int j = 0; j < S->NbConstraints; ++j) {
2692 if (value_negz_p(S->Constraint[j][1+nvar+i]))
2693 continue;
2694 if (value_one_p(S->Constraint[j][1+nvar+i]))
2695 continue;
2696 value_lcm(lcm, S->Constraint[j][1+nvar+i], &lcm);
2697 }
2698
2699 for (int j = 0; j < S->NbConstraints; ++j) {
2700 if (value_negz_p(S->Constraint[j][1+nvar+i]))
2701 continue;
2702 if (value_one_p(S->Constraint[j][1+nvar+i]))
2703 continue;
2704 value_division(f, lcm, S->Constraint[j][1+nvar+i]);
2705 Vector_Scale(S->Constraint[j], M->p[c], f, S->Dimension+2);
2706 value_substract(M->p[c][S->Dimension+1],
2707 M->p[c][S->Dimension+1],
2708 lcm);
2709 value_increment(M->p[c][S->Dimension+1],
2710 M->p[c][S->Dimension+1]);
2711 ++c;
2712 }
2713 Polyhedron *O = S;
2714 S = AddConstraints(M->p[0], c, S, MaxRays);
2715 if (O != P)
2716 Polyhedron_Free(O);
2717 Matrix_Free(M);
2718 if (emptyQ(S)) {
2719 Polyhedron_Free(S);
2720 value_clear(lcm);
2721 value_clear(f);
2722 return 0;
2723 }
2724 }
2725 value_clear(lcm);
2726 value_clear(f);
2727
2728 #ifdef DEBUG_ER
2729 fprintf(stderr, "\nER: Sure\n");
2730 #endif /* DEBUG_ER */
2731
2732 return split_sure(P, S, exist, nparam, MaxRays);
2733 }
2734
2735 static evalue* enumerate_sure2(Polyhedron *P,
2736 unsigned exist, unsigned nparam, unsigned MaxRays)
2737 {
2738 int nvar = P->Dimension - exist - nparam;
2739 int r;
2740 for (r = 0; r < P->NbRays; ++r)
2741 if (value_one_p(P->Ray[r][0]) &&
2742 value_one_p(P->Ray[r][P->Dimension+1]))
2743 break;
2744
2745 if (r >= P->NbRays)
2746 return 0;
2747
2748 Matrix *M = Matrix_Alloc(nvar + 1 + nparam, P->Dimension+2);
2749 for (int i = 0; i < nvar; ++i)
2750 value_set_si(M->p[i][1+i], 1);
2751 for (int i = 0; i < nparam; ++i)
2752 value_set_si(M->p[i+nvar][1+nvar+exist+i], 1);
2753 Vector_Copy(P->Ray[r]+1+nvar, M->p[nvar+nparam]+1+nvar, exist);
2754 value_set_si(M->p[nvar+nparam][0], 1);
2755 value_set_si(M->p[nvar+nparam][P->Dimension+1], 1);
2756 Polyhedron * F = Rays2Polyhedron(M, MaxRays);
2757 Matrix_Free(M);
2758
2759 Polyhedron *I = DomainIntersection(F, P, MaxRays);
2760 Polyhedron_Free(F);
2761
2762 #ifdef DEBUG_ER
2763 fprintf(stderr, "\nER: Sure2\n");
2764 #endif /* DEBUG_ER */
2765
2766 return split_sure(P, I, exist, nparam, MaxRays);
2767 }
2768
2769 static evalue* enumerate_cyclic(Polyhedron *P,
2770 unsigned exist, unsigned nparam,
2771 evalue * EP, int r, int p, unsigned MaxRays)
2772 {
2773 int nvar = P->Dimension - exist - nparam;
2774
2775 /* If EP in its fractional maps only contains references
2776 * to the remainder parameter with appropriate coefficients
2777 * then we could in principle avoid adding existentially
2778 * quantified variables to the validity domains.
2779 * We'd have to replace the remainder by m { p/m }
2780 * and multiply with an appropriate factor that is one
2781 * only in the appropriate range.
2782 * This last multiplication can be avoided if EP
2783 * has a single validity domain with no (further)
2784 * constraints on the remainder parameter
2785 */
2786
2787 Matrix *CT = Matrix_Alloc(nparam+1, nparam+3);
2788 Matrix *M = Matrix_Alloc(1, 1+nparam+3);
2789 for (int j = 0; j < nparam; ++j)
2790 if (j != p)
2791 value_set_si(CT->p[j][j], 1);
2792 value_set_si(CT->p[p][nparam+1], 1);
2793 value_set_si(CT->p[nparam][nparam+2], 1);
2794 value_set_si(M->p[0][1+p], -1);
2795 value_absolute(M->p[0][1+nparam], P->Ray[0][1+nvar+exist+p]);
2796 value_set_si(M->p[0][1+nparam+1], 1);
2797 Polyhedron *CEq = Constraints2Polyhedron(M, 1);
2798 Matrix_Free(M);
2799 addeliminatedparams_enum(EP, CT, CEq, MaxRays, nparam);
2800 Polyhedron_Free(CEq);
2801 Matrix_Free(CT);
2802
2803 return EP;
2804 }
2805
2806 static void enumerate_vd_add_ray(evalue *EP, Matrix *Rays, unsigned MaxRays)
2807 {
2808 if (value_notzero_p(EP->d))
2809 return;
2810
2811 assert(EP->x.p->type == partition);
2812 assert(EP->x.p->pos == EVALUE_DOMAIN(EP->x.p->arr[0])->Dimension);
2813 for (int i = 0; i < EP->x.p->size/2; ++i) {
2814 Polyhedron *D = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
2815 Polyhedron *N = DomainAddRays(D, Rays, MaxRays);
2816 EVALUE_SET_DOMAIN(EP->x.p->arr[2*i], N);
2817 Domain_Free(D);
2818 }
2819 }
2820
2821 static evalue* enumerate_line(Polyhedron *P,
2822 unsigned exist, unsigned nparam, unsigned MaxRays)
2823 {
2824 if (P->NbBid == 0)
2825 return 0;
2826
2827 #ifdef DEBUG_ER
2828 fprintf(stderr, "\nER: Line\n");
2829 #endif /* DEBUG_ER */
2830
2831 int nvar = P->Dimension - exist - nparam;
2832 int i, j;
2833 for (i = 0; i < nparam; ++i)
2834 if (value_notzero_p(P->Ray[0][1+nvar+exist+i]))
2835 break;
2836 assert(i < nparam);
2837 for (j = i+1; j < nparam; ++j)
2838 if (value_notzero_p(P->Ray[0][1+nvar+exist+i]))
2839 break;
2840 assert(j >= nparam); // for now
2841
2842 Matrix *M = Matrix_Alloc(2, P->Dimension+2);
2843 value_set_si(M->p[0][0], 1);
2844 value_set_si(M->p[0][1+nvar+exist+i], 1);
2845 value_set_si(M->p[1][0], 1);
2846 value_set_si(M->p[1][1+nvar+exist+i], -1);
2847 value_absolute(M->p[1][1+P->Dimension], P->Ray[0][1+nvar+exist+i]);
2848 value_decrement(M->p[1][1+P->Dimension], M->p[1][1+P->Dimension]);
2849 Polyhedron *S = AddConstraints(M->p[0], 2, P, MaxRays);
2850 evalue *EP = barvinok_enumerate_e(S, exist, nparam, MaxRays);
2851 Polyhedron_Free(S);
2852 Matrix_Free(M);
2853
2854 return enumerate_cyclic(P, exist, nparam, EP, 0, i, MaxRays);
2855 }
2856
2857 static int single_param_pos(Polyhedron*P, unsigned exist, unsigned nparam,
2858 int r)
2859 {
2860 int nvar = P->Dimension - exist - nparam;
2861 if (First_Non_Zero(P->Ray[r]+1, nvar) != -1)
2862 return -1;
2863 int i = First_Non_Zero(P->Ray[r]+1+nvar+exist, nparam);
2864 if (i == -1)
2865 return -1;
2866 if (First_Non_Zero(P->Ray[r]+1+nvar+exist+1, nparam-i-1) != -1)
2867 return -1;
2868 return i;
2869 }
2870
2871 static evalue* enumerate_remove_ray(Polyhedron *P, int r,
2872 unsigned exist, unsigned nparam, unsigned MaxRays)
2873 {
2874 #ifdef DEBUG_ER
2875 fprintf(stderr, "\nER: RedundantRay\n");
2876 #endif /* DEBUG_ER */
2877
2878 Value one;
2879 value_init(one);
2880 value_set_si(one, 1);
2881 int len = P->NbRays-1;
2882 Matrix *M = Matrix_Alloc(2 * len, P->Dimension+2);
2883 Vector_Copy(P->Ray[0], M->p[0], r * (P->Dimension+2));
2884 Vector_Copy(P->Ray[r+1], M->p[r], (len-r) * (P->Dimension+2));
2885 for (int j = 0; j < P->NbRays; ++j) {
2886 if (j == r)
2887 continue;
2888 Vector_Combine(P->Ray[j], P->Ray[r], M->p[len+j-(j>r)],
2889 one, P->Ray[j][P->Dimension+1], P->Dimension+2);
2890 }
2891
2892 P = Rays2Polyhedron(M, MaxRays);
2893 Matrix_Free(M);
2894 evalue *EP = barvinok_enumerate_e(P, exist, nparam, MaxRays);
2895 Polyhedron_Free(P);
2896 value_clear(one);
2897
2898 return EP;
2899 }
2900
2901 static evalue* enumerate_redundant_ray(Polyhedron *P,
2902 unsigned exist, unsigned nparam, unsigned MaxRays)
2903 {
2904 assert(P->NbBid == 0);
2905 int nvar = P->Dimension - exist - nparam;
2906 Value m;
2907 value_init(m);
2908
2909 for (int r = 0; r < P->NbRays; ++r) {
2910 if (value_notzero_p(P->Ray[r][P->Dimension+1]))
2911 continue;
2912 int i1 = single_param_pos(P, exist, nparam, r);
2913 if (i1 == -1)
2914 continue;
2915 for (int r2 = r+1; r2 < P->NbRays; ++r2) {
2916 if (value_notzero_p(P->Ray[r2][P->Dimension+1]))
2917 continue;
2918 int i2 = single_param_pos(P, exist, nparam, r2);
2919 if (i2 == -1)
2920 continue;
2921 if (i1 != i2)
2922 continue;
2923
2924 value_division(m, P->Ray[r][1+nvar+exist+i1],
2925 P->Ray[r2][1+nvar+exist+i1]);
2926 value_multiply(m, m, P->Ray[r2][1+nvar+exist+i1]);
2927 /* r2 divides r => r redundant */
2928 if (value_eq(m, P->Ray[r][1+nvar+exist+i1])) {
2929 value_clear(m);
2930 return enumerate_remove_ray(P, r, exist, nparam, MaxRays);
2931 }
2932
2933 value_division(m, P->Ray[r2][1+nvar+exist+i1],
2934 P->Ray[r][1+nvar+exist+i1]);
2935 value_multiply(m, m, P->Ray[r][1+nvar+exist+i1]);
2936 /* r divides r2 => r2 redundant */
2937 if (value_eq(m, P->Ray[r2][1+nvar+exist+i1])) {
2938 value_clear(m);
2939 return enumerate_remove_ray(P, r2, exist, nparam, MaxRays);
2940 }
2941 }
2942 }
2943 value_clear(m);
2944 return 0;
2945 }
2946
2947 static Polyhedron *upper_bound(Polyhedron *P,
2948 int pos, Value *max, Polyhedron **R)
2949 {
2950 Value v;
2951 int r;
2952 value_init(v);
2953
2954 *R = 0;
2955 Polyhedron *N;
2956 Polyhedron *B = 0;
2957 for (Polyhedron *Q = P; Q; Q = N) {
2958 N = Q->next;
2959 for (r = 0; r < P->NbRays; ++r) {
2960 if (value_zero_p(P->Ray[r][P->Dimension+1]) &&
2961 value_pos_p(P->Ray[r][1+pos]))
2962 break;
2963 }
2964 if (r < P->NbRays) {
2965 Q->next = *R;
2966 *R = Q;
2967 continue;
2968 } else {
2969 Q->next = B;
2970 B = Q;
2971 }
2972 for (r = 0; r < P->NbRays; ++r) {
2973 if (value_zero_p(P->Ray[r][P->Dimension+1]))
2974 continue;
2975 mpz_fdiv_q(v, P->Ray[r][1+pos], P->Ray[r][1+P->Dimension]);
2976 if ((!Q->next && r == 0) || value_gt(v, *max))
2977 value_assign(*max, v);
2978 }
2979 }
2980 value_clear(v);
2981 return B;
2982 }
2983
2984 static evalue* enumerate_ray(Polyhedron *P,
2985 unsigned exist, unsigned nparam, unsigned MaxRays)
2986 {
2987 assert(P->NbBid == 0);
2988 int nvar = P->Dimension - exist - nparam;
2989
2990 int r;
2991 for (r = 0; r < P->NbRays; ++r)
2992 if (value_zero_p(P->Ray[r][P->Dimension+1]))
2993 break;
2994 if (r >= P->NbRays)
2995 return 0;
2996
2997 int r2;
2998 for (r2 = r+1; r2 < P->NbRays; ++r2)
2999 if (value_zero_p(P->Ray[r2][P->Dimension+1]))
3000 break;
3001 if (r2 < P->NbRays) {
3002 if (nvar > 0)
3003 return enumerate_sum(P, exist, nparam, MaxRays);
3004 }
3005
3006 #ifdef DEBUG_ER
3007 fprintf(stderr, "\nER: Ray\n");
3008 #endif /* DEBUG_ER */
3009
3010 Value m;
3011 Value one;
3012 value_init(m);
3013 value_init(one);
3014 value_set_si(one, 1);
3015 int i = single_param_pos(P, exist, nparam, r);
3016 assert(i != -1); // for now;
3017
3018 Matrix *M = Matrix_Alloc(P->NbRays, P->Dimension+2);
3019 for (int j = 0; j < P->NbRays; ++j) {
3020 Vector_Combine(P->Ray[j], P->Ray[r], M->p[j],
3021 one, P->Ray[j][P->Dimension+1], P->Dimension+2);
3022 }
3023 Polyhedron *S = Rays2Polyhedron(M, MaxRays);
3024 Matrix_Free(M);
3025 Polyhedron *D = DomainDifference(P, S, MaxRays);
3026 Polyhedron_Free(S);
3027 // Polyhedron_Print(stderr, P_VALUE_FMT, D);
3028 assert(value_pos_p(P->Ray[r][1+nvar+exist+i])); // for now
3029 Polyhedron *R;
3030 D = upper_bound(D, nvar+exist+i, &m, &R);
3031 assert(D);
3032 Domain_Free(D);
3033
3034 M = Matrix_Alloc(2, P->Dimension+2);
3035 value_set_si(M->p[0][0], 1);
3036 value_set_si(M->p[1][0], 1);
3037 value_set_si(M->p[0][1+nvar+exist+i], -1);
3038 value_set_si(M->p[1][1+nvar+exist+i], 1);
3039 value_assign(M->p[0][1+P->Dimension], m);
3040 value_oppose(M->p[1][1+P->Dimension], m);
3041 value_addto(M->p[1][1+P->Dimension], M->p[1][1+P->Dimension],
3042 P->Ray[r][1+nvar+exist+i]);
3043 value_decrement(M->p[1][1+P->Dimension], M->p[1][1+P->Dimension]);
3044 // Matrix_Print(stderr, P_VALUE_FMT, M);
3045 D = AddConstraints(M->p[0], 2, P, MaxRays);
3046 // Polyhedron_Print(stderr, P_VALUE_FMT, D);
3047 value_substract(M->p[0][1+P->Dimension], M->p[0][1+P->Dimension],
3048 P->Ray[r][1+nvar+exist+i]);
3049 // Matrix_Print(stderr, P_VALUE_FMT, M);
3050 S = AddConstraints(M->p[0], 1, P, MaxRays);
3051 // Polyhedron_Print(stderr, P_VALUE_FMT, S);
3052 Matrix_Free(M);
3053
3054 evalue *EP = barvinok_enumerate_e(D, exist, nparam, MaxRays);
3055 Polyhedron_Free(D);
3056 value_clear(one);
3057 value_clear(m);
3058
3059 if (value_notone_p(P->Ray[r][1+nvar+exist+i]))
3060 EP = enumerate_cyclic(P, exist, nparam, EP, r, i, MaxRays);
3061 else {
3062 M = Matrix_Alloc(1, nparam+2);
3063 value_set_si(M->p[0][0], 1);
3064 value_set_si(M->p[0][1+i], 1);
3065 enumerate_vd_add_ray(EP, M, MaxRays);
3066 Matrix_Free(M);
3067 }
3068
3069 if (!emptyQ(S)) {
3070 evalue *E = barvinok_enumerate_e(S, exist, nparam, MaxRays);
3071 eadd(E, EP);
3072 free_evalue_refs(E);
3073 free(E);
3074 }
3075 Polyhedron_Free(S);
3076
3077 if (R) {
3078 assert(nvar == 0);
3079 evalue *ER = enumerate_or(R, exist, nparam, MaxRays);
3080 eor(ER, EP);
3081 free_evalue_refs(ER);
3082 free(ER);
3083 }
3084
3085 return EP;
3086 }
3087
3088 static evalue* new_zero_ep()
3089 {
3090 evalue *EP;
3091 ALLOC(evalue, EP);
3092 value_init(EP->d);
3093 evalue_set_si(EP, 0, 1);
3094 return EP;
3095 }
3096
3097 static evalue* enumerate_vd(Polyhedron **PA,
3098 unsigned exist, unsigned nparam, unsigned MaxRays)
3099 {
3100 Polyhedron *P = *PA;
3101 int nvar = P->Dimension - exist - nparam;
3102 Param_Polyhedron *PP = NULL;
3103 Polyhedron *C = Universe_Polyhedron(nparam);
3104 Polyhedron *CEq;
3105 Matrix *CT;
3106 Polyhedron *PR = P;
3107 PP = Polyhedron2Param_SimplifiedDomain(&PR,C,MaxRays,&CEq,&CT);
3108 Polyhedron_Free(C);
3109
3110 int nd;
3111 Param_Domain *D, *last;
3112 Value c;
3113 value_init(c);
3114 for (nd = 0, D=PP->D; D; D=D->next, ++nd)
3115 ;
3116
3117 Polyhedron **VD = new Polyhedron_p[nd];
3118 Polyhedron **fVD = new Polyhedron_p[nd];
3119 for(nd = 0, D=PP->D; D; D=D->next) {
3120 Polyhedron *rVD = reduce_domain(D->Domain, CT, CEq,
3121 fVD, nd, MaxRays);
3122 if (!rVD)
3123 continue;
3124
3125 VD[nd++] = rVD;
3126 last = D;
3127 }
3128
3129 evalue *EP = 0;
3130
3131 if (nd == 0)
3132 EP = new_zero_ep();
3133
3134 /* This doesn't seem to have any effect */
3135 if (nd == 1) {
3136 Polyhedron *CA = align_context(VD[0], P->Dimension, MaxRays);
3137 Polyhedron *O = P;
3138 P = DomainIntersection(P, CA, MaxRays);
3139 if (O != *PA)
3140 Polyhedron_Free(O);
3141 Polyhedron_Free(CA);
3142 if (emptyQ(P))
3143 EP = new_zero_ep();
3144 }
3145
3146 if (!EP && CT->NbColumns != CT->NbRows) {
3147 Polyhedron *CEqr = DomainImage(CEq, CT, MaxRays);
3148 Polyhedron *CA = align_context(CEqr, PR->Dimension, MaxRays);
3149 Polyhedron *I = DomainIntersection(PR, CA, MaxRays);
3150 Polyhedron_Free(CEqr);
3151 Polyhedron_Free(CA);
3152 #ifdef DEBUG_ER
3153 fprintf(stderr, "\nER: Eliminate\n");
3154 #endif /* DEBUG_ER */
3155 nparam -= CT->NbColumns - CT->NbRows;
3156 EP = barvinok_enumerate_e(I, exist, nparam, MaxRays);
3157 nparam += CT->NbColumns - CT->NbRows;
3158 addeliminatedparams_enum(EP, CT, CEq, MaxRays, nparam);
3159 Polyhedron_Free(I);
3160 }
3161 if (PR != *PA)
3162 Polyhedron_Free(PR);
3163 PR = 0;
3164
3165 if (!EP && nd > 1) {
3166 #ifdef DEBUG_ER
3167 fprintf(stderr, "\nER: VD\n");
3168 #endif /* DEBUG_ER */
3169 for (int i = 0; i < nd; ++i) {
3170 Polyhedron *CA = align_context(VD[i], P->Dimension, MaxRays);
3171 Polyhedron *I = DomainIntersection(P, CA, MaxRays);
3172
3173 if (i == 0)
3174 EP = barvinok_enumerate_e(I, exist, nparam, MaxRays);
3175 else {
3176 evalue *E = barvinok_enumerate_e(I, exist, nparam, MaxRays);
3177 eadd(E, EP);
3178 free_evalue_refs(E);
3179 free(E);
3180 }
3181 Polyhedron_Free(I);
3182 Polyhedron_Free(CA);
3183 }
3184 }
3185
3186 for (int i = 0; i < nd; ++i) {
3187 Polyhedron_Free(VD[i]);
3188 Polyhedron_Free(fVD[i]);
3189 }
3190 delete [] VD;
3191 delete [] fVD;
3192 value_clear(c);
3193
3194 if (!EP && nvar == 0) {
3195 Value f;
3196 value_init(f);
3197 Param_Vertices *V, *V2;
3198 Matrix* M = Matrix_Alloc(1, P->Dimension+2);
3199
3200 FORALL_PVertex_in_ParamPolyhedron(V, last, PP) {
3201 bool found = false;
3202 FORALL_PVertex_in_ParamPolyhedron(V2, last, PP) {
3203 if (V == V2) {
3204 found = true;
3205 continue;
3206 }
3207 if (!found)
3208 continue;
3209 for (int i = 0; i < exist; ++i) {
3210 value_oppose(f, V->Vertex->p[i][nparam+1]);
3211 Vector_Combine(V->Vertex->p[i],
3212 V2->Vertex->p[i],
3213 M->p[0] + 1 + nvar + exist,
3214 V2->Vertex->p[i][nparam+1],
3215 f,
3216 nparam+1);
3217 int j;
3218 for (j = 0; j < nparam; ++j)
3219 if (value_notzero_p(M->p[0][1+nvar+exist+j]))
3220 break;
3221 if (j >= nparam)
3222 continue;
3223 ConstraintSimplify(M->p[0], M->p[0],
3224 P->Dimension+2, &f);
3225 value_set_si(M->p[0][0], 0);
3226 Polyhedron *para = AddConstraints(M->p[0], 1, P,
3227 MaxRays);
3228 if (emptyQ(para)) {
3229 Polyhedron_Free(para);
3230 continue;
3231 }
3232 Polyhedron *pos, *neg;
3233 value_set_si(M->p[0][0], 1);
3234 value_decrement(M->p[0][P->Dimension+1],
3235 M->p[0][P->Dimension+1]);
3236 neg = AddConstraints(M->p[0], 1, P, MaxRays);
3237 value_set_si(f, -1);
3238 Vector_Scale(M->p[0]+1, M->p[0]+1, f,
3239 P->Dimension+1);
3240 value_decrement(M->p[0][P->Dimension+1],
3241 M->p[0][P->Dimension+1]);
3242 value_decrement(M->p[0][P->Dimension+1],
3243 M->p[0][P->Dimension+1]);
3244 pos = AddConstraints(M->p[0], 1, P, MaxRays);
3245 if (emptyQ(neg) && emptyQ(pos)) {
3246 Polyhedron_Free(para);
3247 Polyhedron_Free(pos);
3248 Polyhedron_Free(neg);
3249 continue;
3250 }
3251 #ifdef DEBUG_ER
3252 fprintf(stderr, "\nER: Order\n");
3253 #endif /* DEBUG_ER */
3254 EP = barvinok_enumerate_e(para, exist, nparam, MaxRays);
3255 evalue *E;
3256 if (!emptyQ(pos)) {
3257 E = barvinok_enumerate_e(pos, exist, nparam, MaxRays);
3258 eadd(E, EP);
3259 free_evalue_refs(E);
3260 free(E);
3261 }
3262 if (!emptyQ(neg)) {
3263 E = barvinok_enumerate_e(neg, exist, nparam, MaxRays);
3264 eadd(E, EP);
3265 free_evalue_refs(E);
3266 free(E);
3267 }
3268 Polyhedron_Free(para);
3269 Polyhedron_Free(pos);
3270 Polyhedron_Free(neg);
3271 break;
3272 }
3273 if (EP)
3274 break;
3275 } END_FORALL_PVertex_in_ParamPolyhedron;
3276 if (EP)
3277 break;
3278 } END_FORALL_PVertex_in_ParamPolyhedron;
3279
3280 if (!EP) {
3281 /* Search for vertex coordinate to split on */
3282 /* First look for one independent of the parameters */
3283 FORALL_PVertex_in_ParamPolyhedron(V, last, PP) {
3284 for (int i = 0; i < exist; ++i) {
3285 int j;
3286 for (j = 0; j < nparam; ++j)
3287 if (value_notzero_p(V->Vertex->p[i][j]))
3288 break;
3289 if (j < nparam)
3290 continue;
3291 value_set_si(M->p[0][0], 1);
3292 Vector_Set(M->p[0]+1, 0, nvar+exist);
3293 Vector_Copy(V->Vertex->p[i],
3294 M->p[0] + 1 + nvar + exist, nparam+1);
3295 value_oppose(M->p[0][1+nvar+i],
3296 V->Vertex->p[i][nparam+1]);
3297
3298 Polyhedron *pos, *neg;
3299 value_set_si(M->p[0][0], 1);
3300 value_decrement(M->p[0][P->Dimension+1],
3301 M->p[0][P->Dimension+1]);
3302 neg = AddConstraints(M->p[0], 1, P, MaxRays);
3303 value_set_si(f, -1);
3304 Vector_Scale(M->p[0]+1, M->p[0]+1, f,
3305 P->Dimension+1);
3306 value_decrement(M->p[0][P->Dimension+1],
3307 M->p[0][P->Dimension+1]);
3308 value_decrement(M->p[0][P->Dimension+1],
3309 M->p[0][P->Dimension+1]);
3310 pos = AddConstraints(M->p[0], 1, P, MaxRays);
3311 if (emptyQ(neg) || emptyQ(pos)) {
3312 Polyhedron_Free(pos);
3313 Polyhedron_Free(neg);
3314 continue;
3315 }
3316 Polyhedron_Free(pos);
3317 value_increment(M->p[0][P->Dimension+1],
3318 M->p[0][P->Dimension+1]);
3319 pos = AddConstraints(M->p[0], 1, P, MaxRays);
3320 #ifdef DEBUG_ER
3321 fprintf(stderr, "\nER: Vertex\n");
3322 #endif /* DEBUG_ER */
3323 pos->next = neg;
3324 EP = enumerate_or(pos, exist, nparam, MaxRays);
3325 break;
3326 }
3327 if (EP)
3328 break;
3329 } END_FORALL_PVertex_in_ParamPolyhedron;
3330 }
3331
3332 if (!EP) {
3333 /* Search for vertex coordinate to split on */
3334 /* Now look for one that depends on the parameters */
3335 FORALL_PVertex_in_ParamPolyhedron(V, last, PP) {
3336 for (int i = 0; i < exist; ++i) {
3337 value_set_si(M->p[0][0], 1);
3338 Vector_Set(M->p[0]+1, 0, nvar+exist);
3339 Vector_Copy(V->Vertex->p[i],
3340 M->p[0] + 1 + nvar + exist, nparam+1);
3341 value_oppose(M->p[0][1+nvar+i],
3342 V->Vertex->p[i][nparam+1]);
3343
3344 Polyhedron *pos, *neg;
3345 value_set_si(M->p[0][0], 1);
3346 value_decrement(M->p[0][P->Dimension+1],
3347 M->p[0][P->Dimension+1]);
3348 neg = AddConstraints(M->p[0], 1, P, MaxRays);
3349 value_set_si(f, -1);
3350 Vector_Scale(M->p[0]+1, M->p[0]+1, f,
3351 P->Dimension+1);
3352 value_decrement(M->p[0][P->Dimension+1],
3353 M->p[0][P->Dimension+1]);
3354 value_decrement(M->p[0][P->Dimension+1],
3355 M->p[0][P->Dimension+1]);
3356 pos = AddConstraints(M->p[0], 1, P, MaxRays);
3357 if (emptyQ(neg) || emptyQ(pos)) {
3358 Polyhedron_Free(pos);
3359 Polyhedron_Free(neg);
3360 continue;
3361 }
3362 Polyhedron_Free(pos);
3363 value_increment(M->p[0][P->Dimension+1],
3364 M->p[0][P->Dimension+1]);
3365 pos = AddConstraints(M->p[0], 1, P, MaxRays);
3366 #ifdef DEBUG_ER
3367 fprintf(stderr, "\nER: ParamVertex\n");
3368 #endif /* DEBUG_ER */
3369 pos->next = neg;
3370 EP = enumerate_or(pos, exist, nparam, MaxRays);
3371 break;
3372 }
3373 if (EP)
3374 break;
3375 } END_FORALL_PVertex_in_ParamPolyhedron;
3376 }
3377
3378 Matrix_Free(M);
3379 value_clear(f);
3380 }
3381
3382 if (CEq)
3383 Polyhedron_Free(CEq);
3384 if (CT)
3385 Matrix_Free(CT);
3386 if (PP)
3387 Param_Polyhedron_Free(PP);
3388 *PA = P;
3389
3390 return EP;
3391 }
3392
3393 #ifndef HAVE_PIPLIB
3394 evalue *barvinok_enumerate_pip(Polyhedron *P,
3395 unsigned exist, unsigned nparam, unsigned MaxRays)
3396 {
3397 return 0;
3398 }
3399 #else
3400 evalue *barvinok_enumerate_pip(Polyhedron *P,
3401 unsigned exist, unsigned nparam, unsigned MaxRays)
3402 {
3403 int nvar = P->Dimension - exist - nparam;
3404 evalue *EP = new_zero_ep();
3405 Polyhedron *Q, *N, *T = 0;
3406 Value min, tmp;
3407 value_init(min);
3408 value_init(tmp);
3409
3410 #ifdef DEBUG_ER
3411 fprintf(stderr, "\nER: PIP\n");
3412 #endif /* DEBUG_ER */
3413
3414 for (int i = 0; i < P->Dimension; ++i) {
3415 bool pos = false;
3416 bool neg = false;
3417 bool posray = false;
3418 bool negray = false;
3419 value_set_si(min, 0);
3420 for (int j = 0; j < P->NbRays; ++j) {
3421 if (value_pos_p(P->Ray[j][1+i])) {
3422 pos = true;
3423 if (value_zero_p(P->Ray[j][1+P->Dimension]))
3424 posray = true;
3425 } else if (value_neg_p(P->Ray[j][1+i])) {
3426 neg = true;
3427 if (value_zero_p(P->Ray[j][1+P->Dimension]))
3428 negray = true;
3429 else {
3430 mpz_fdiv_q(tmp,
3431 P->Ray[j][1+i], P->Ray[j][1+P->Dimension]);
3432 if (value_lt(tmp, min))
3433 value_assign(min, tmp);
3434 }
3435 }
3436 }
3437 if (pos && neg) {
3438 assert(!(posray && negray));
3439 assert(!negray); // for now
3440 Polyhedron *O = T ? T : P;
3441 /* shift by a safe amount */
3442 Matrix *M = Matrix_Alloc(O->NbRays, O->Dimension+2);
3443 Vector_Copy(O->Ray[0], M->p[0], O->NbRays * (O->Dimension+2));
3444 for (int j = 0; j < P->NbRays; ++j) {
3445 if (value_notzero_p(M->p[j][1+P->Dimension])) {
3446 value_multiply(tmp, min, M->p[j][1+P->Dimension]);
3447 value_substract(M->p[j][1+i], M->p[j][1+i], tmp);
3448 }
3449 }
3450 if (T)
3451 Polyhedron_Free(T);
3452 T = Rays2Polyhedron(M, MaxRays);
3453 Matrix_Free(M);
3454 } else if (neg) {
3455 /* negating a parameter requires that we substitute in the
3456 * sign again afterwards.
3457 * Disallow for now.
3458 */
3459 assert(i < nvar+exist);
3460 if (!T)
3461 T = Polyhedron_Copy(P);
3462 for (int j = 0; j < T->NbRays; ++j)
3463 value_oppose(T->Ray[j][1+i], T->Ray[j][1+i]);
3464 for (int j = 0; j < T->NbConstraints; ++j)
3465 value_oppose(T->Constraint[j][1+i], T->Constraint[j][1+i]);
3466 }
3467 }
3468 value_clear(min);
3469 value_clear(tmp);
3470
3471 Polyhedron *D = pip_lexmin(T ? T : P, exist, nparam);
3472 for (Q = D; Q; Q = N) {
3473 N = Q->next;
3474 Q->next = 0;
3475 evalue *E;
3476 exist = Q->Dimension - nvar - nparam;
3477 E = barvinok_enumerate_e(Q, exist, nparam, MaxRays);
3478 Polyhedron_Free(Q);
3479 eadd(E, EP);
3480 free_evalue_refs(E);
3481 free(E);
3482 }
3483
3484 if (T)
3485 Polyhedron_Free(T);
3486
3487 return EP;
3488 }
3489 #endif
3490
3491
3492 static bool is_single(Value *row, int pos, int len)
3493 {
3494 return First_Non_Zero(row, pos) == -1 &&
3495 First_Non_Zero(row+pos+1, len-pos-1) == -1;
3496 }
3497
3498 static evalue* barvinok_enumerate_e_r(Polyhedron *P,
3499 unsigned exist, unsigned nparam, unsigned MaxRays);
3500
3501 #ifdef DEBUG_ER
3502 static int er_level = 0;
3503
3504 evalue* barvinok_enumerate_e(Polyhedron *P,
3505 unsigned exist, unsigned nparam, unsigned MaxRays)
3506 {
3507 fprintf(stderr, "\nER: level %i\n", er_level);
3508 int nvar = P->Dimension - exist - nparam;
3509 fprintf(stderr, "%d %d %d\n", nvar, exist, nparam);
3510
3511 Polyhedron_Print(stderr, P_VALUE_FMT, P);
3512 ++er_level;
3513 P = DomainConstraintSimplify(Polyhedron_Copy(P), MaxRays);
3514 evalue *EP = barvinok_enumerate_e_r(P, exist, nparam, MaxRays);
3515 Polyhedron_Free(P);
3516 --er_level;
3517 return EP;
3518 }
3519 #else
3520 evalue* barvinok_enumerate_e(Polyhedron *P,
3521 unsigned exist, unsigned nparam, unsigned MaxRays)
3522 {
3523 P = DomainConstraintSimplify(Polyhedron_Copy(P), MaxRays);
3524 evalue *EP = barvinok_enumerate_e_r(P, exist, nparam, MaxRays);
3525 Polyhedron_Free(P);
3526 return EP;
3527 }
3528 #endif
3529
3530 static evalue* barvinok_enumerate_e_r(Polyhedron *P,
3531 unsigned exist, unsigned nparam, unsigned MaxRays)
3532 {
3533 if (exist == 0) {
3534 Polyhedron *U = Universe_Polyhedron(nparam);
3535 evalue *EP = barvinok_enumerate_ev(P, U, MaxRays);
3536 //char *param_name[] = {"P", "Q", "R", "S", "T" };
3537 //print_evalue(stdout, EP, param_name);
3538 Polyhedron_Free(U);
3539 return EP;
3540 }
3541
3542 int nvar = P->Dimension - exist - nparam;
3543 int len = P->Dimension + 2;
3544
3545 if (emptyQ(P))
3546 return new_zero_ep();
3547
3548 if (nvar == 0 && nparam == 0) {
3549 evalue *EP = new_zero_ep();
3550 barvinok_count(P, &EP->x.n, MaxRays);
3551 if (value_pos_p(EP->x.n))
3552 value_set_si(EP->x.n, 1);
3553 return EP;
3554 }
3555
3556 int r;
3557 for (r = 0; r < P->NbRays; ++r)
3558 if (value_zero_p(P->Ray[r][0]) ||
3559 value_zero_p(P->Ray[r][P->Dimension+1])) {
3560 int i;
3561 for (i = 0; i < nvar; ++i)
3562 if (value_notzero_p(P->Ray[r][i+1]))
3563 break;
3564 if (i >= nvar)
3565 continue;
3566 for (i = nvar + exist; i < nvar + exist + nparam; ++i)
3567 if (value_notzero_p(P->Ray[r][i+1]))
3568 break;
3569 if (i >= nvar + exist + nparam)
3570 break;
3571 }
3572 if (r < P->NbRays) {
3573 evalue *EP = new_zero_ep();
3574 value_set_si(EP->x.n, -1);
3575 return EP;
3576 }
3577
3578 int first;
3579 for (r = 0; r < P->NbEq; ++r)
3580 if ((first = First_Non_Zero(P->Constraint[r]+1+nvar, exist)) != -1)
3581 break;
3582 if (r < P->NbEq) {
3583 if (First_Non_Zero(P->Constraint[r]+1+nvar+first+1,
3584 exist-first-1) != -1) {
3585 Polyhedron *T = rotate_along(P, r, nvar, exist, MaxRays);
3586 #ifdef DEBUG_ER
3587 fprintf(stderr, "\nER: Equality\n");
3588 #endif /* DEBUG_ER */
3589 evalue *EP = barvinok_enumerate_e(T, exist-1, nparam, MaxRays);
3590 Polyhedron_Free(T);
3591 return EP;
3592 } else {
3593 #ifdef DEBUG_ER
3594 fprintf(stderr, "\nER: Fixed\n");
3595 #endif /* DEBUG_ER */
3596 if (first == 0)
3597 return barvinok_enumerate_e(P, exist-1, nparam, MaxRays);
3598 else {
3599 Polyhedron *T = Polyhedron_Copy(P);
3600 SwapColumns(T, nvar+1, nvar+1+first);
3601 evalue *EP = barvinok_enumerate_e(T, exist-1, nparam, MaxRays);
3602 Polyhedron_Free(T);
3603 return EP;
3604 }
3605 }
3606 }
3607
3608 Vector *row = Vector_Alloc(len);
3609 value_set_si(row->p[0], 1);
3610
3611 Value f;
3612 value_init(f);
3613
3614 enum constraint* info = new constraint[exist];
3615 for (int i = 0; i < exist; ++i) {
3616 info[i] = ALL_POS;
3617 for (int l = P->NbEq; l < P->NbConstraints; ++l) {
3618 if (value_negz_p(P->Constraint[l][nvar+i+1]))
3619 continue;
3620 bool l_parallel = is_single(P->Constraint[l]+nvar+1, i, exist);
3621 for (int u = P->NbEq; u < P->NbConstraints; ++u) {
3622 if (value_posz_p(P->Constraint[u][nvar+i+1]))
3623 continue;
3624 bool lu_parallel = l_parallel ||
3625 is_single(P->Constraint[u]+nvar+1, i, exist);
3626 value_oppose(f, P->Constraint[u][nvar+i+1]);
3627 Vector_Combine(P->Constraint[l]+1, P->Constraint[u]+1, row->p+1,
3628 f, P->Constraint[l][nvar+i+1], len-1);
3629 if (!(info[i] & INDEPENDENT)) {
3630 int j;
3631 for (j = 0; j < exist; ++j)
3632 if (j != i && value_notzero_p(row->p[nvar+j+1]))
3633 break;
3634 if (j == exist) {
3635 //printf("independent: i: %d, l: %d, u: %d\n", i, l, u);
3636 info[i] = (constraint)(info[i] | INDEPENDENT);
3637 }
3638 }
3639 if (info[i] & ALL_POS) {
3640 value_addto(row->p[len-1], row->p[len-1],
3641 P->Constraint[l][nvar+i+1]);
3642 value_addto(row->p[len-1], row->p[len-1], f);
3643 value_multiply(f, f, P->Constraint[l][nvar+i+1]);
3644 value_substract(row->p[len-1], row->p[len-1], f);
3645 value_decrement(row->p[len-1], row->p[len-1]);
3646 ConstraintSimplify(row->p, row->p, len, &f);
3647 value_set_si(f, -1);
3648 Vector_Scale(row->p+1, row->p+1, f, len-1);
3649 value_decrement(row->p[len-1], row->p[len-1]);
3650 Polyhedron *T = AddConstraints(row->p, 1, P, MaxRays);
3651 if (!emptyQ(T)) {
3652 //printf("not all_pos: i: %d, l: %d, u: %d\n", i, l, u);
3653 info[i] = (constraint)(info[i] ^ ALL_POS);
3654 }
3655 //puts("pos remainder");
3656 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
3657 Polyhedron_Free(T);
3658 }
3659 if (!(info[i] & ONE_NEG)) {
3660 if (lu_parallel) {
3661 negative_test_constraint(P->Constraint[l],
3662 P->Constraint[u],
3663 row->p, nvar+i, len, &f);
3664 oppose_constraint(row->p, len, &f);
3665 Polyhedron *T = AddConstraints(row->p, 1, P, MaxRays);
3666 if (emptyQ(T)) {
3667 //printf("one_neg i: %d, l: %d, u: %d\n", i, l, u);
3668 info[i] = (constraint)(info[i] | ONE_NEG);
3669 }
3670 //puts("neg remainder");
3671 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
3672 Polyhedron_Free(T);
3673 } else if (!(info[i] & ROT_NEG)) {
3674 if (parallel_constraints(P->Constraint[l],
3675 P->Constraint[u],
3676 row->p, nvar, exist)) {
3677 negative_test_constraint7(P->Constraint[l],
3678 P->Constraint[u],
3679 row->p, nvar, exist,
3680 len, &f);
3681 oppose_constraint(row->p, len, &f);
3682 Polyhedron *T = AddConstraints(row->p, 1, P, MaxRays);
3683 if (emptyQ(T)) {
3684 // printf("rot_neg i: %d, l: %d, u: %d\n", i, l, u);
3685 info[i] = (constraint)(info[i] | ROT_NEG);
3686 r = l;
3687 }
3688 //puts("neg remainder");
3689 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
3690 Polyhedron_Free(T);
3691 }
3692 }
3693 }
3694 if (!(info[i] & ALL_POS) && (info[i] & (ONE_NEG | ROT_NEG)))
3695 goto next;
3696 }
3697 }
3698 if (info[i] & ALL_POS)
3699 break;
3700 next:
3701 ;
3702 }
3703
3704 /*
3705 for (int i = 0; i < exist; ++i)
3706 printf("%i: %i\n", i, info[i]);
3707 */
3708 for (int i = 0; i < exist; ++i)
3709 if (info[i] & ALL_POS) {
3710 #ifdef DEBUG_ER
3711 fprintf(stderr, "\nER: Positive\n");
3712 #endif /* DEBUG_ER */
3713 // Eliminate
3714 // Maybe we should chew off some of the fat here
3715 Matrix *M = Matrix_Alloc(P->Dimension, P->Dimension+1);
3716 for (int j = 0; j < P->Dimension; ++j)
3717 value_set_si(M->p[j][j + (j >= i+nvar)], 1);
3718 Polyhedron *T = Polyhedron_Image(P, M, MaxRays);
3719 Matrix_Free(M);
3720 evalue *EP = barvinok_enumerate_e(T, exist-1, nparam, MaxRays);
3721 Polyhedron_Free(T);
3722 value_clear(f);
3723 Vector_Free(row);
3724 delete [] info;
3725 return EP;
3726 }
3727 for (int i = 0; i < exist; ++i)
3728 if (info[i] & ONE_NEG) {
3729 #ifdef DEBUG_ER
3730 fprintf(stderr, "\nER: Negative\n");
3731 #endif /* DEBUG_ER */
3732 Vector_Free(row);
3733 value_clear(f);
3734 delete [] info;
3735 if (i == 0)
3736 return barvinok_enumerate_e(P, exist-1, nparam, MaxRays);
3737 else {
3738 Polyhedron *T = Polyhedron_Copy(P);
3739 SwapColumns(T, nvar+1, nvar+1+i);
3740 evalue *EP = barvinok_enumerate_e(T, exist-1, nparam, MaxRays);
3741 Polyhedron_Free(T);
3742 return EP;
3743 }
3744 }
3745 for (int i = 0; i < exist; ++i)
3746 if (info[i] & ROT_NEG) {
3747 #ifdef DEBUG_ER
3748 fprintf(stderr, "\nER: Rotate\n");
3749 #endif /* DEBUG_ER */
3750 Vector_Free(row);
3751 value_clear(f);
3752 delete [] info;
3753 Polyhedron *T = rotate_along(P, r, nvar, exist, MaxRays);
3754 evalue *EP = barvinok_enumerate_e(T, exist-1, nparam, MaxRays);
3755 Polyhedron_Free(T);
3756 return EP;
3757 }
3758 for (int i = 0; i < exist; ++i)
3759 if (info[i] & INDEPENDENT) {
3760 Polyhedron *pos, *neg;
3761
3762 /* Find constraint again and split off negative part */
3763
3764 if (SplitOnVar(P, i, nvar, len, exist, MaxRays,
3765 row, f, true, &pos, &neg)) {
3766 #ifdef DEBUG_ER
3767 fprintf(stderr, "\nER: Split\n");
3768 #endif /* DEBUG_ER */
3769
3770 evalue *EP =
3771 barvinok_enumerate_e(neg, exist-1, nparam, MaxRays);
3772 evalue *E =
3773 barvinok_enumerate_e(pos, exist, nparam, MaxRays);
3774 eadd(E, EP);
3775 free_evalue_refs(E);
3776 free(E);
3777 Polyhedron_Free(neg);
3778 Polyhedron_Free(pos);
3779 value_clear(f);
3780 Vector_Free(row);
3781 delete [] info;
3782 return EP;
3783 }
3784 }
3785 delete [] info;
3786
3787 Polyhedron *O = P;
3788 Polyhedron *F;
3789
3790 evalue *EP;
3791
3792 EP = enumerate_line(P, exist, nparam, MaxRays);
3793 if (EP)
3794 goto out;
3795
3796 EP = barvinok_enumerate_pip(P, exist, nparam, MaxRays);
3797 if (EP)
3798 goto out;
3799
3800 EP = enumerate_redundant_ray(P, exist, nparam, MaxRays);
3801 if (EP)
3802 goto out;
3803
3804 EP = enumerate_sure(P, exist, nparam, MaxRays);
3805 if (EP)
3806 goto out;
3807
3808 EP = enumerate_ray(P, exist, nparam, MaxRays);
3809 if (EP)
3810 goto out;
3811
3812 EP = enumerate_sure2(P, exist, nparam, MaxRays);
3813 if (EP)
3814 goto out;
3815
3816 F = unfringe(P, MaxRays);
3817 if (!PolyhedronIncludes(F, P)) {
3818 #ifdef DEBUG_ER
3819 fprintf(stderr, "\nER: Fringed\n");
3820 #endif /* DEBUG_ER */
3821 EP = barvinok_enumerate_e(F, exist, nparam, MaxRays);
3822 Polyhedron_Free(F);
3823 goto out;
3824 }
3825 Polyhedron_Free(F);
3826
3827 if (nparam)
3828 EP = enumerate_vd(&P, exist, nparam, MaxRays);
3829 if (EP)
3830 goto out2;
3831
3832 if (nvar != 0) {
3833 EP = enumerate_sum(P, exist, nparam, MaxRays);
3834 goto out2;
3835 }
3836
3837 assert(nvar == 0);
3838
3839 int i;
3840 Polyhedron *pos, *neg;
3841 for (i = 0; i < exist; ++i)
3842 if (SplitOnVar(P, i, nvar, len, exist, MaxRays,
3843 row, f, false, &pos, &neg))
3844 break;
3845
3846 assert (i < exist);
3847
3848 pos->next = neg;
3849 EP = enumerate_or(pos, exist, nparam, MaxRays);
3850
3851 out2:
3852 if (O != P)
3853 Polyhedron_Free(P);
3854
3855 out:
3856 value_clear(f);
3857 Vector_Free(row);
3858 return EP;
3859 }
3860
3861 gen_fun * barvinok_series(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
3862 {
3863 Polyhedron ** vcone;
3864 Polyhedron *CA;
3865 unsigned nparam = C->Dimension;
3866 unsigned dim, nvar;
3867 vec_ZZ sign;
3868 int ncone = 0;
3869 sign.SetLength(ncone);
3870
3871 CA = align_context(C, P->Dimension, MaxRays);
3872 P = DomainIntersection(P, CA, MaxRays);
3873 Polyhedron_Free(CA);
3874
3875 assert(!Polyhedron_is_infinite(P, nparam));
3876 assert(P->NbBid == 0);
3877 assert(Polyhedron_has_positive_rays(P, nparam));
3878 assert(P->NbEq == 0);
3879
3880 dim = P->Dimension;
3881 nvar = dim - nparam;
3882 vcone = new Polyhedron_p[P->NbRays];
3883
3884 for (int j = 0; j < P->NbRays; ++j) {
3885 if (!value_pos_p(P->Ray[j][dim+1]))
3886 continue;
3887
3888 int npos, nneg;
3889 Polyhedron *C = supporting_cone(P, j);
3890 decompose(C, &vcone[j], &npos, &nneg, MaxRays);
3891 ncone += npos + nneg;
3892 sign.SetLength(ncone);
3893 for (int k = 0; k < npos; ++k)
3894 sign[ncone-nneg-k-1] = 1;
3895 for (int k = 0; k < nneg; ++k)
3896 sign[ncone-k-1] = -1;
3897 }
3898
3899 mat_ZZ rays;
3900 rays.SetDims(ncone * dim, nvar);
3901 int r = 0;
3902 for (int j = 0; j < P->NbRays; ++j) {
3903 if (!value_pos_p(P->Ray[j][dim+1]))
3904 continue;
3905
3906 for (Polyhedron *i = vcone[j]; i; i = i->next) {
3907 add_rays(rays, i, &r, nvar);
3908 }
3909 }
3910 rays.SetDims(r, nvar);
3911 vec_ZZ lambda;
3912 nonorthog(rays, lambda);
3913 //randomvector(P, lambda, nvar);
3914
3915 /*
3916 cout << "rays: " << rays;
3917 cout << "lambda: " << lambda;
3918 */
3919
3920 int f = 0;
3921 ZZ num_s;
3922 vec_ZZ num_p;
3923 num_p.SetLength(nparam);
3924 vec_ZZ vertex;
3925 vec_ZZ den_s;
3926 den_s.SetLength(dim);
3927 vec_ZZ den_p;
3928 den_p.SetLength(dim);
3929 mat_ZZ den;
3930 den.SetDims(dim, nparam);
3931 ZZ one;
3932 one = 1;
3933 mpq_t count;
3934 mpq_init(count);
3935
3936 gen_fun * gf = new gen_fun;
3937
3938 rays.SetDims(dim, nvar);
3939
3940 for (int j = 0; j < P->NbRays; ++j) {
3941 if (!value_pos_p(P->Ray[j][dim+1]))
3942 continue;
3943
3944 for (Polyhedron *i = vcone[j]; i; i = i->next, ++f) {
3945 lattice_point(P->Ray[j]+1, i, vertex);
3946 int k = 0;
3947 num_s = 0;
3948 for ( ; k < nvar; ++k)
3949 num_s += vertex[k] * lambda[k];
3950 for ( ; k < dim; ++k)
3951 num_p[k-nvar] = vertex[k];
3952
3953 int r = 0;
3954 add_rays(rays, i, &r, nvar, true);
3955 for (r = 0; r < dim; ++r)
3956 values2zz(i->Ray[r]+1+nvar, den[r], nparam);
3957 den_s = rays * lambda;
3958
3959 normalize(sign[f], num_s, num_p, den_s, den_p, den);
3960
3961 int only_param = 0;
3962 int no_param = 0;
3963 for (int k = 0; k < dim; ++k) {
3964 if (den_p[k] == 0)
3965 ++no_param;
3966 else if (den_s[k] == 0)
3967 ++only_param;
3968 }
3969 if (no_param == 0) {
3970 for (int k = 0; k < dim; ++k)
3971 if (den_p[k] == -1)
3972 den[k] = -den[k];
3973 gf->add(sign[f], one, num_p, den);
3974 } else if (no_param + only_param == dim) {
3975 int k, l;
3976 mat_ZZ pden;
3977 pden.SetDims(only_param, nparam);
3978
3979 for (k = 0, l = 0; k < dim; ++k)
3980 if (den_p[k] != 0)
3981 pden[l++] = den[k];
3982
3983 for (k = 0; k < dim; ++k)
3984 if (den_s[k] != 0)
3985 break;
3986
3987 dpoly n(no_param, num_s);
3988 dpoly d(no_param, den_s[k], 1);
3989 for ( ; ++k < dim; k)
3990 if (den_s[k] != 0) {
3991 dpoly fact(no_param, den_s[k], 1);
3992 d *= fact;
3993 }
3994
3995 mpq_set_si(count, 0, 1);
3996 n.div(d, count, sign[f]);
3997
3998 ZZ qn, qd;
3999 value2zz(mpq_numref(count), qn);
4000 value2zz(mpq_denref(count), qd);
4001
4002 gf->add(qn, qd, num_p, pden);
4003 } else {
4004 int k, l;
4005 dpoly_r * r = 0;
4006 mat_ZZ pden;
4007 pden.SetDims(only_param, nparam);
4008
4009 for (k = 0, l = 0; k < dim; ++k)
4010 if (den_s[k] == 0)
4011 pden[l++] = den[k];
4012
4013 for (k = 0; k < dim; ++k)
4014 if (den_p[k] == 0)
4015 break;
4016
4017 dpoly n(no_param, num_s);
4018 dpoly d(no_param, den_s[k], 1);
4019 for ( ; ++k < dim; )
4020 if (den_p[k] == 0) {
4021 dpoly fact(no_param, den_s[k], 1);
4022 d *= fact;
4023 }
4024
4025 for (k = 0; k < dim; ++k) {
4026 if (den_s[k] == 0 || den_p[k] == 0)
4027 continue;
4028
4029 dpoly pd(no_param-1, den_s[k], 1);
4030 int s = den_p[k] < 0 ? -1 : 1;
4031
4032 if (r == 0)
4033 r = new dpoly_r(n, pd, k, s, dim);
4034 else
4035 assert(0); // for now
4036 }
4037
4038 r->div(d, sign[f], gf, pden, den, num_p);
4039 }
4040
4041 /*
4042 cout << "sign: " << sign[f];
4043 cout << "num_s: " << num_s;
4044 cout << "num_p: " << num_p;
4045 cout << "den_s: " << den_s;
4046 cout << "den_p: " << den_p;
4047 cout << "den: " << den;
4048 cout << "only_param: " << only_param;
4049 cout << "no_param: " << no_param;
4050 cout << endl;
4051 */
4052
4053 }
4054 }
4055
4056 mpq_clear(count);
4057
4058 return gf;
4059 }
lattice point enumeration library