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