?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
define partition type
[barvinok.git] / barvinok.cc
blobaf620faff82a58ef1caa1cb38a2fbaf8f3d5d2b8
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 }
15 #include "config.h"
16 #include <barvinok.h>
17
18 #ifdef NTL_STD_CXX
19 using namespace NTL;
20 #endif
21 using std::cout;
22 using std::endl;
23 using std::vector;
24 using std::deque;
25 using std::string;
26 using std::ostringstream;
27
28 #define ALLOC(p) (((long *) (p))[0])
29 #define SIZE(p) (((long *) (p))[1])
30 #define DATA(p) ((mp_limb_t *) (((long *) (p)) + 2))
31
32 static void value2zz(Value v, ZZ& z)
33 {
34 int sa = v[0]._mp_size;
35 int abs_sa = sa < 0 ? -sa : sa;
36
37 _ntl_gsetlength(&z.rep, abs_sa);
38 mp_limb_t * adata = DATA(z.rep);
39 for (int i = 0; i < abs_sa; ++i)
40 adata[i] = v[0]._mp_d[i];
41 SIZE(z.rep) = sa;
42 }
43
44 static void zz2value(ZZ& z, Value& v)
45 {
46 if (!z.rep) {
47 value_set_si(v, 0);
48 return;
49 }
50
51 int sa = SIZE(z.rep);
52 int abs_sa = sa < 0 ? -sa : sa;
53
54 mp_limb_t * adata = DATA(z.rep);
55 mpz_realloc2(v, __GMP_BITS_PER_MP_LIMB * abs_sa);
56 for (int i = 0; i < abs_sa; ++i)
57 v[0]._mp_d[i] = adata[i];
58 v[0]._mp_size = sa;
59 }
60
61 #undef ALLOC
62 #define ALLOC(p) p = (typeof(p))malloc(sizeof(*p))
63
64 /*
65 * We just ignore the last column and row
66 * If the final element is not equal to one
67 * then the result will actually be a multiple of the input
68 */
69 static void matrix2zz(Matrix *M, mat_ZZ& m, unsigned nr, unsigned nc)
70 {
71 m.SetDims(nr, nc);
72
73 for (int i = 0; i < nr; ++i) {
74 // assert(value_one_p(M->p[i][M->NbColumns - 1]));
75 for (int j = 0; j < nc; ++j) {
76 value2zz(M->p[i][j], m[i][j]);
77 }
78 }
79 }
80
81 static void values2zz(Value *p, vec_ZZ& v, int len)
82 {
83 v.SetLength(len);
84
85 for (int i = 0; i < len; ++i) {
86 value2zz(p[i], v[i]);
87 }
88 }
89
90 /*
91 * We add a 0 at the end, because we need it afterwards
92 */
93 static Vector * zz2vector(vec_ZZ& v)
94 {
95 Vector *vec = Vector_Alloc(v.length()+1);
96 assert(vec);
97 for (int i = 0; i < v.length(); ++i)
98 zz2value(v[i], vec->p[i]);
99
100 value_set_si(vec->p[v.length()], 0);
101
102 return vec;
103 }
104
105 static void rays(mat_ZZ& r, Polyhedron *C)
106 {
107 unsigned dim = C->NbRays - 1; /* don't count zero vertex */
108 assert(C->NbRays - 1 == C->Dimension);
109 r.SetDims(dim, dim);
110 ZZ tmp;
111
112 int i, c;
113 for (i = 0, c = 0; i < dim; ++i)
114 if (value_zero_p(C->Ray[i][dim+1])) {
115 for (int j = 0; j < dim; ++j) {
116 value2zz(C->Ray[i][j+1], tmp);
117 r[j][c] = tmp;
118 }
119 ++c;
120 }
121 }
122
123 static Matrix * rays(Polyhedron *C)
124 {
125 unsigned dim = C->NbRays - 1; /* don't count zero vertex */
126 assert(C->NbRays - 1 == C->Dimension);
127
128 Matrix *M = Matrix_Alloc(dim+1, dim+1);
129 assert(M);
130
131 int i, c;
132 for (i = 0, c = 0; i <= dim && c < dim; ++i)
133 if (value_zero_p(C->Ray[i][dim+1])) {
134 Vector_Copy(C->Ray[i] + 1, M->p[c], dim);
135 value_set_si(M->p[c++][dim], 0);
136 }
137 assert(c == dim);
138 value_set_si(M->p[dim][dim], 1);
139
140 return M;
141 }
142
143 static Matrix * rays2(Polyhedron *C)
144 {
145 unsigned dim = C->NbRays - 1; /* don't count zero vertex */
146 assert(C->NbRays - 1 == C->Dimension);
147
148 Matrix *M = Matrix_Alloc(dim, dim);
149 assert(M);
150
151 int i, c;
152 for (i = 0, c = 0; i <= dim && c < dim; ++i)
153 if (value_zero_p(C->Ray[i][dim+1]))
154 Vector_Copy(C->Ray[i] + 1, M->p[c++], dim);
155 assert(c == dim);
156
157 return M;
158 }
159
160 /*
161 * Returns the largest absolute value in the vector
162 */
163 static ZZ max(vec_ZZ& v)
164 {
165 ZZ max = abs(v[0]);
166 for (int i = 1; i < v.length(); ++i)
167 if (abs(v[i]) > max)
168 max = abs(v[i]);
169 return max;
170 }
171
172 class cone {
173 public:
174 cone(Matrix *M) {
175 Cone = 0;
176 Rays = Matrix_Copy(M);
177 set_det();
178 }
179 cone(Polyhedron *C) {
180 Cone = Polyhedron_Copy(C);
181 Rays = rays(C);
182 set_det();
183 }
184 void set_det() {
185 mat_ZZ A;
186 matrix2zz(Rays, A, Rays->NbRows - 1, Rays->NbColumns - 1);
187 det = determinant(A);
188 Value v;
189 value_init(v);
190 zz2value(det, v);
191 value_clear(v);
192 }
193
194 Vector* short_vector(vec_ZZ& lambda) {
195 Matrix *M = Matrix_Copy(Rays);
196 Matrix *inv = Matrix_Alloc(M->NbRows, M->NbColumns);
197 int ok = Matrix_Inverse(M, inv);
198 assert(ok);
199 Matrix_Free(M);
200
201 ZZ det2;
202 mat_ZZ B;
203 mat_ZZ U;
204 matrix2zz(inv, B, inv->NbRows - 1, inv->NbColumns - 1);
205 long r = LLL(det2, B, U);
206
207 ZZ min = max(B[0]);
208 int index = 0;
209 for (int i = 1; i < B.NumRows(); ++i) {
210 ZZ tmp = max(B[i]);
211 if (tmp < min) {
212 min = tmp;
213 index = i;
214 }
215 }
216
217 Matrix_Free(inv);
218
219 lambda = B[index];
220
221 Vector *z = zz2vector(U[index]);
222 Value tmp;
223 value_init(tmp);
224 Polyhedron *C = poly();
225 int i;
226 for (i = 0; i < C->NbConstraints; ++i) {
227 Inner_Product(z->p, C->Constraint[i]+1, z->Size-1, &tmp);
228 if (value_pos_p(tmp))
229 break;
230 }
231 if (i == C->NbConstraints) {
232 value_set_si(tmp, -1);
233 Vector_Scale(z->p, z->p, tmp, z->Size-1);
234 }
235 value_clear(tmp);
236 return z;
237 }
238
239 ~cone() {
240 Polyhedron_Free(Cone);
241 Matrix_Free(Rays);
242 }
243
244 Polyhedron *poly() {
245 if (!Cone) {
246 Matrix *M = Matrix_Alloc(Rays->NbRows+1, Rays->NbColumns+1);
247 for (int i = 0; i < Rays->NbRows; ++i) {
248 Vector_Copy(Rays->p[i], M->p[i]+1, Rays->NbColumns);
249 value_set_si(M->p[i][0], 1);
250 }
251 Vector_Set(M->p[Rays->NbRows]+1, 0, Rays->NbColumns-1);
252 value_set_si(M->p[Rays->NbRows][0], 1);
253 value_set_si(M->p[Rays->NbRows][Rays->NbColumns], 1);
254 Cone = Rays2Polyhedron(M, M->NbRows+1);
255 assert(Cone->NbConstraints == Cone->NbRays);
256 Matrix_Free(M);
257 }
258 return Cone;
259 }
260
261 ZZ det;
262 Polyhedron *Cone;
263 Matrix *Rays;
264 };
265
266 class dpoly {
267 public:
268 vec_ZZ coeff;
269 dpoly(int d, ZZ& degree, int offset = 0) {
270 coeff.SetLength(d+1);
271
272 int min = d + offset;
273 if (degree < ZZ(INIT_VAL, min))
274 min = to_int(degree);
275
276 ZZ c = ZZ(INIT_VAL, 1);
277 if (!offset)
278 coeff[0] = c;
279 for (int i = 1; i <= min; ++i) {
280 c *= (degree -i + 1);
281 c /= i;
282 coeff[i-offset] = c;
283 }
284 }
285 void operator *= (dpoly& f) {
286 assert(coeff.length() == f.coeff.length());
287 vec_ZZ old = coeff;
288 coeff = f.coeff[0] * coeff;
289 for (int i = 1; i < coeff.length(); ++i)
290 for (int j = 0; i+j < coeff.length(); ++j)
291 coeff[i+j] += f.coeff[i] * old[j];
292 }
293 void div(dpoly& d, mpq_t count, ZZ& sign) {
294 int len = coeff.length();
295 Value tmp;
296 value_init(tmp);
297 mpq_t* c = new mpq_t[coeff.length()];
298 mpq_t qtmp;
299 mpq_init(qtmp);
300 for (int i = 0; i < len; ++i) {
301 mpq_init(c[i]);
302 zz2value(coeff[i], tmp);
303 mpq_set_z(c[i], tmp);
304
305 for (int j = 1; j <= i; ++j) {
306 zz2value(d.coeff[j], tmp);
307 mpq_set_z(qtmp, tmp);
308 mpq_mul(qtmp, qtmp, c[i-j]);
309 mpq_sub(c[i], c[i], qtmp);
310 }
311
312 zz2value(d.coeff[0], tmp);
313 mpq_set_z(qtmp, tmp);
314 mpq_div(c[i], c[i], qtmp);
315 }
316 if (sign == -1)
317 mpq_sub(count, count, c[len-1]);
318 else
319 mpq_add(count, count, c[len-1]);
320
321 value_clear(tmp);
322 mpq_clear(qtmp);
323 for (int i = 0; i < len; ++i)
324 mpq_clear(c[i]);
325 delete [] c;
326 }
327 };
328
329 class dpoly_n {
330 public:
331 Matrix *coeff;
332 ~dpoly_n() {
333 Matrix_Free(coeff);
334 }
335 dpoly_n(int d, ZZ& degree_0, ZZ& degree_1, int offset = 0) {
336 Value d0, d1;
337 value_init(d0);
338 value_init(d1);
339 zz2value(degree_0, d0);
340 zz2value(degree_1, d1);
341 coeff = Matrix_Alloc(d+1, d+1+1);
342 value_set_si(coeff->p[0][0], 1);
343 value_set_si(coeff->p[0][d+1], 1);
344 for (int i = 1; i <= d; ++i) {
345 value_multiply(coeff->p[i][0], coeff->p[i-1][0], d0);
346 Vector_Combine(coeff->p[i-1], coeff->p[i-1]+1, coeff->p[i]+1,
347 d1, d0, i);
348 value_set_si(coeff->p[i][d+1], i);
349 value_multiply(coeff->p[i][d+1], coeff->p[i][d+1], coeff->p[i-1][d+1]);
350 value_decrement(d0, d0);
351 }
352 value_clear(d0);
353 value_clear(d1);
354 }
355 void div(dpoly& d, Vector *count, ZZ& sign) {
356 int len = coeff->NbRows;
357 Matrix * c = Matrix_Alloc(coeff->NbRows, coeff->NbColumns);
358 Value tmp;
359 value_init(tmp);
360 for (int i = 0; i < len; ++i) {
361 Vector_Copy(coeff->p[i], c->p[i], len+1);
362 for (int j = 1; j <= i; ++j) {
363 zz2value(d.coeff[j], tmp);
364 value_multiply(tmp, tmp, c->p[i][len]);
365 value_oppose(tmp, tmp);
366 Vector_Combine(c->p[i], c->p[i-j], c->p[i],
367 c->p[i-j][len], tmp, len);
368 value_multiply(c->p[i][len], c->p[i][len], c->p[i-j][len]);
369 }
370 zz2value(d.coeff[0], tmp);
371 value_multiply(c->p[i][len], c->p[i][len], tmp);
372 }
373 if (sign == -1) {
374 value_set_si(tmp, -1);
375 Vector_Scale(c->p[len-1], count->p, tmp, len);
376 value_assign(count->p[len], c->p[len-1][len]);
377 } else
378 Vector_Copy(c->p[len-1], count->p, len+1);
379 Vector_Normalize(count->p, len+1);
380 value_clear(tmp);
381 Matrix_Free(c);
382 }
383 };
384
385 /*
386 * Barvinok's Decomposition of a simplicial cone
387 *
388 * Returns two lists of polyhedra
389 */
390 void barvinok_decompose(Polyhedron *C, Polyhedron **ppos, Polyhedron **pneg)
391 {
392 Polyhedron *pos = *ppos, *neg = *pneg;
393 vector<cone *> nonuni;
394 cone * c = new cone(C);
395 ZZ det = c->det;
396 int s = sign(det);
397 assert(det != 0);
398 if (abs(det) > 1) {
399 nonuni.push_back(c);
400 } else {
401 Polyhedron *p = Polyhedron_Copy(c->Cone);
402 p->next = pos;
403 pos = p;
404 delete c;
405 }
406 vec_ZZ lambda;
407 while (!nonuni.empty()) {
408 c = nonuni.back();
409 nonuni.pop_back();
410 Vector* v = c->short_vector(lambda);
411 for (int i = 0; i < c->Rays->NbRows - 1; ++i) {
412 if (lambda[i] == 0)
413 continue;
414 Matrix* M = Matrix_Copy(c->Rays);
415 Vector_Copy(v->p, M->p[i], v->Size);
416 cone * pc = new cone(M);
417 assert (pc->det != 0);
418 if (abs(pc->det) > 1) {
419 assert(abs(pc->det) < abs(c->det));
420 nonuni.push_back(pc);
421 } else {
422 Polyhedron *p = pc->poly();
423 pc->Cone = 0;
424 if (sign(pc->det) == s) {
425 p->next = pos;
426 pos = p;
427 } else {
428 p->next = neg;
429 neg = p;
430 }
431 delete pc;
432 }
433 Matrix_Free(M);
434 }
435 Vector_Free(v);
436 delete c;
437 }
438 *ppos = pos;
439 *pneg = neg;
440 }
441
442 /*
443 * Returns a single list of npos "positive" cones followed by nneg
444 * "negative" cones.
445 * The input cone is freed
446 */
447 void decompose(Polyhedron *cone, Polyhedron **parts, int *npos, int *nneg, unsigned MaxRays)
448 {
449 Polyhedron_Polarize(cone);
450 if (cone->NbRays - 1 != cone->Dimension) {
451 Polyhedron *tmp = cone;
452 cone = triangularize_cone(cone, MaxRays);
453 Polyhedron_Free(tmp);
454 }
455 Polyhedron *polpos = NULL, *polneg = NULL;
456 *npos = 0; *nneg = 0;
457 for (Polyhedron *Polar = cone; Polar; Polar = Polar->next)
458 barvinok_decompose(Polar, &polpos, &polneg);
459
460 Polyhedron *last;
461 for (Polyhedron *i = polpos; i; i = i->next) {
462 Polyhedron_Polarize(i);
463 ++*npos;
464 last = i;
465 }
466 for (Polyhedron *i = polneg; i; i = i->next) {
467 Polyhedron_Polarize(i);
468 ++*nneg;
469 }
470 if (last) {
471 last->next = polneg;
472 *parts = polpos;
473 } else
474 *parts = polneg;
475 Domain_Free(cone);
476 }
477
478 const int MAX_TRY=10;
479 /*
480 * Searches for a vector that is not othogonal to any
481 * of the rays in rays.
482 */
483 static void nonorthog(mat_ZZ& rays, vec_ZZ& lambda)
484 {
485 int dim = rays.NumCols();
486 bool found = false;
487 lambda.SetLength(dim);
488 for (int i = 2; !found && i <= 50*dim; i+=4) {
489 for (int j = 0; j < MAX_TRY; ++j) {
490 for (int k = 0; k < dim; ++k) {
491 int r = random_int(i)+2;
492 int v = (2*(r%2)-1) * (r >> 1);
493 lambda[k] = v;
494 }
495 int k = 0;
496 for (; k < rays.NumRows(); ++k)
497 if (lambda * rays[k] == 0)
498 break;
499 if (k == rays.NumRows()) {
500 found = true;
501 break;
502 }
503 }
504 }
505 assert(found);
506 }
507
508 static void add_rays(mat_ZZ& rays, Polyhedron *i, int *r)
509 {
510 unsigned dim = i->Dimension;
511 for (int k = 0; k < i->NbRays; ++k) {
512 if (!value_zero_p(i->Ray[k][dim+1]))
513 continue;
514 values2zz(i->Ray[k]+1, rays[(*r)++], dim);
515 }
516 }
517
518 void lattice_point(Value* values, Polyhedron *i, vec_ZZ& lambda, ZZ& num)
519 {
520 vec_ZZ vertex;
521 unsigned dim = i->Dimension;
522 if(!value_one_p(values[dim])) {
523 Matrix* Rays = rays(i);
524 Matrix *inv = Matrix_Alloc(Rays->NbRows, Rays->NbColumns);
525 int ok = Matrix_Inverse(Rays, inv);
526 assert(ok);
527 Matrix_Free(Rays);
528 Rays = rays(i);
529 Vector *lambda = Vector_Alloc(dim+1);
530 Vector_Matrix_Product(values, inv, lambda->p);
531 Matrix_Free(inv);
532 for (int j = 0; j < dim; ++j)
533 mpz_cdiv_q(lambda->p[j], lambda->p[j], lambda->p[dim]);
534 value_set_si(lambda->p[dim], 1);
535 Vector *A = Vector_Alloc(dim+1);
536 Vector_Matrix_Product(lambda->p, Rays, A->p);
537 Vector_Free(lambda);
538 Matrix_Free(Rays);
539 values2zz(A->p, vertex, dim);
540 Vector_Free(A);
541 } else
542 values2zz(values, vertex, dim);
543
544 num = vertex * lambda;
545 }
546
547 static evalue *term(int param, ZZ& c, Value *den = NULL)
548 {
549 evalue *EP = new evalue();
550 value_init(EP->d);
551 value_set_si(EP->d,0);
552 EP->x.p = new_enode(polynomial, 2, param + 1);
553 evalue_set_si(&EP->x.p->arr[0], 0, 1);
554 value_init(EP->x.p->arr[1].x.n);
555 if (den == NULL)
556 value_set_si(EP->x.p->arr[1].d, 1);
557 else
558 value_assign(EP->x.p->arr[1].d, *den);
559 zz2value(c, EP->x.p->arr[1].x.n);
560 return EP;
561 }
562
563 static void vertex_period(
564 Polyhedron *i, vec_ZZ& lambda, Matrix *T,
565 Value lcm, int p, Vector *val,
566 evalue *E, evalue* ev,
567 ZZ& offset)
568 {
569 unsigned nparam = T->NbRows - 1;
570 unsigned dim = i->Dimension;
571 Value tmp;
572 ZZ nump;
573
574 if (p == nparam) {
575 ZZ num, l;
576 Vector * values = Vector_Alloc(dim + 1);
577 Vector_Matrix_Product(val->p, T, values->p);
578 value_assign(values->p[dim], lcm);
579 lattice_point(values->p, i, lambda, num);
580 value2zz(lcm, l);
581 num *= l;
582 num += offset;
583 value_init(ev->x.n);
584 zz2value(num, ev->x.n);
585 value_assign(ev->d, lcm);
586 Vector_Free(values);
587 return;
588 }
589
590 value_init(tmp);
591 vec_ZZ vertex;
592 values2zz(T->p[p], vertex, dim);
593 nump = vertex * lambda;
594 if (First_Non_Zero(val->p, p) == -1) {
595 value_assign(tmp, lcm);
596 evalue *ET = term(p, nump, &tmp);
597 eadd(ET, E);
598 free_evalue_refs(ET);
599 delete ET;
600 }
601
602 value_assign(tmp, lcm);
603 if (First_Non_Zero(T->p[p], dim) != -1)
604 Vector_Gcd(T->p[p], dim, &tmp);
605 Gcd(tmp, lcm, &tmp);
606 if (value_lt(tmp, lcm)) {
607 ZZ count;
608
609 value_division(tmp, lcm, tmp);
610 value_set_si(ev->d, 0);
611 ev->x.p = new_enode(periodic, VALUE_TO_INT(tmp), p+1);
612 value2zz(tmp, count);
613 do {
614 value_decrement(tmp, tmp);
615 --count;
616 ZZ new_offset = offset - count * nump;
617 value_assign(val->p[p], tmp);
618 vertex_period(i, lambda, T, lcm, p+1, val, E,
619 &ev->x.p->arr[VALUE_TO_INT(tmp)], new_offset);
620 } while (value_pos_p(tmp));
621 } else
622 vertex_period(i, lambda, T, lcm, p+1, val, E, ev, offset);
623 value_clear(tmp);
624 }
625
626 static void mask_r(Matrix *f, int nr, Vector *lcm, int p, Vector *val, evalue *ev)
627 {
628 unsigned nparam = lcm->Size;
629
630 if (p == nparam) {
631 Vector * prod = Vector_Alloc(f->NbRows);
632 Matrix_Vector_Product(f, val->p, prod->p);
633 int isint = 1;
634 for (int i = 0; i < nr; ++i) {
635 value_modulus(prod->p[i], prod->p[i], f->p[i][nparam+1]);
636 isint &= value_zero_p(prod->p[i]);
637 }
638 value_set_si(ev->d, 1);
639 value_init(ev->x.n);
640 value_set_si(ev->x.n, isint);
641 Vector_Free(prod);
642 return;
643 }
644
645 Value tmp;
646 value_init(tmp);
647 if (value_one_p(lcm->p[p]))
648 mask_r(f, nr, lcm, p+1, val, ev);
649 else {
650 value_assign(tmp, lcm->p[p]);
651 value_set_si(ev->d, 0);
652 ev->x.p = new_enode(periodic, VALUE_TO_INT(tmp), p+1);
653 do {
654 value_decrement(tmp, tmp);
655 value_assign(val->p[p], tmp);
656 mask_r(f, nr, lcm, p+1, val, &ev->x.p->arr[VALUE_TO_INT(tmp)]);
657 } while (value_pos_p(tmp));
658 }
659 value_clear(tmp);
660 }
661
662 static evalue *multi_monom(vec_ZZ& p)
663 {
664 evalue *X = new evalue();
665 value_init(X->d);
666 value_init(X->x.n);
667 unsigned nparam = p.length()-1;
668 zz2value(p[nparam], X->x.n);
669 value_set_si(X->d, 1);
670 for (int i = 0; i < nparam; ++i) {
671 if (p[i] == 0)
672 continue;
673 evalue *T = term(i, p[i]);
674 eadd(T, X);
675 free_evalue_refs(T);
676 delete T;
677 }
678 return X;
679 }
680
681 #ifdef USE_MODULO
682 static void mask(Matrix *f, evalue *factor)
683 {
684 int nr = f->NbRows, nc = f->NbColumns;
685 int n;
686 bool found = false;
687 for (n = 0; n < nr && value_notzero_p(f->p[n][nc-1]); ++n)
688 if (value_notone_p(f->p[n][nc-1]) &&
689 value_notmone_p(f->p[n][nc-1]))
690 found = true;
691 if (!found)
692 return;
693
694 evalue EP;
695 nr = n;
696
697 for (n = 0; n < nr; ++n) {
698 if (value_one_p(f->p[n][nc-1]) ||
699 value_mone_p(f->p[n][nc-1]))
700 continue;
701 value_init(EP.d);
702 value_set_si(EP.d, 0);
703 EP.x.p = new_enode(relation, 2, 0);
704 value_clear(EP.x.p->arr[1].d);
705 EP.x.p->arr[1] = *factor;
706 evalue *ev = &EP.x.p->arr[0];
707 value_set_si(ev->d, 0);
708 ev->x.p = new_enode(modulo, 3, VALUE_TO_INT(f->p[n][nc-1]));
709 evalue_set_si(&ev->x.p->arr[1], 0, 1);
710 evalue_set_si(&ev->x.p->arr[2], 1, 1);
711 vec_ZZ row;
712 values2zz(f->p[n], row, nc-1);
713 evalue *E = multi_monom(row);
714 value_clear(ev->x.p->arr[0].d);
715 ev->x.p->arr[0] = *E;
716 delete E;
717 *factor = EP;
718 }
719 }
720 #else
721 /*
722 *
723 */
724 static void mask(Matrix *f, evalue *factor)
725 {
726 int nr = f->NbRows, nc = f->NbColumns;
727 int n;
728 bool found = false;
729 for (n = 0; n < nr && value_notzero_p(f->p[n][nc-1]); ++n)
730 if (value_notone_p(f->p[n][nc-1]) &&
731 value_notmone_p(f->p[n][nc-1]))
732 found = true;
733 if (!found)
734 return;
735
736 Value tmp;
737 value_init(tmp);
738 nr = n;
739 unsigned np = nc - 2;
740 Vector *lcm = Vector_Alloc(np);
741 Vector *val = Vector_Alloc(nc);
742 Vector_Set(val->p, 0, nc);
743 value_set_si(val->p[np], 1);
744 Vector_Set(lcm->p, 1, np);
745 for (n = 0; n < nr; ++n) {
746 if (value_one_p(f->p[n][nc-1]) ||
747 value_mone_p(f->p[n][nc-1]))
748 continue;
749 for (int j = 0; j < np; ++j)
750 if (value_notzero_p(f->p[n][j])) {
751 Gcd(f->p[n][j], f->p[n][nc-1], &tmp);
752 value_division(tmp, f->p[n][nc-1], tmp);
753 value_lcm(tmp, lcm->p[j], &lcm->p[j]);
754 }
755 }
756 evalue EP;
757 value_init(EP.d);
758 mask_r(f, nr, lcm, 0, val, &EP);
759 value_clear(tmp);
760 Vector_Free(val);
761 Vector_Free(lcm);
762 emul(&EP,factor);
763 free_evalue_refs(&EP);
764 }
765 #endif
766
767 struct term_info {
768 evalue *E;
769 ZZ constant;
770 ZZ coeff;
771 int pos;
772 };
773
774 #ifdef USE_MODULO
775 void ceil_mod(Value *coef, int len, Value d, ZZ& f, evalue *EP)
776 {
777 Value gcd;
778 value_init(gcd);
779 Value mone;
780 value_init(mone);
781 value_set_si(mone, -1);
782
783 vec_ZZ num;
784
785 Vector_Gcd(coef, len, &gcd);
786 Gcd(gcd, d, &gcd);
787 Vector_AntiScale(coef, coef, gcd, len);
788 Vector_Scale(coef, coef, mone, len);
789 values2zz(coef, num, len);
790
791 value_division(gcd, d, gcd);
792 ZZ g;
793 value2zz(gcd, g);
794 if (value_one_p(gcd))
795 goto out;
796
797 int j;
798 for (j = 0; j < len; ++j)
799 num[j] = num[j] % g;
800 for (j = 0; j < len; ++j)
801 if (num[j] != 0)
802 break;
803 if (j == len)
804 goto out;
805
806 evalue tmp;
807 value_init(tmp.d);
808 evalue_set_si(&tmp, 0, 1);
809
810 if (j < len-1 && num[j] > g/2) {
811 for (int k = j; k < len-1; ++k)
812 if (num[k] != 0)
813 num[k] = g - num[k];
814 num[len-1] = g - 1 - num[len-1];
815 value_assign(tmp.d, gcd);
816 ZZ t = f*(g-1);
817 zz2value(t, tmp.x.n);
818 eadd(&tmp, EP);
819 f = -f;
820 }
821
822 if (j >= len-1) {
823 ZZ t = num[len-1] * f;
824 zz2value(t, tmp.x.n);
825 value_assign(tmp.d, gcd);
826 eadd(&tmp, EP);
827 } else {
828 evalue *E = multi_monom(num);
829
830 evalue EV;
831 value_init(EV.d);
832 value_set_si(EV.d, 0);
833 EV.x.p = new_enode(modulo, 3, VALUE_TO_INT(gcd));
834 evalue_copy(&EV.x.p->arr[0], E);
835 evalue_set_si(&EV.x.p->arr[1], 0, 1);
836 value_init(EV.x.p->arr[2].x.n);
837 zz2value(f, EV.x.p->arr[2].x.n);
838 value_assign(EV.x.p->arr[2].d, gcd);
839
840 eadd(&EV, EP);
841 free_evalue_refs(&EV);
842 free_evalue_refs(E);
843 delete E;
844 }
845
846 free_evalue_refs(&tmp);
847
848 out:
849 value_clear(gcd);
850 value_clear(mone);
851 }
852
853 evalue* ceil3(Value *coef, int len, Value d)
854 {
855 Vector *val = Vector_Alloc(len);
856
857 Value mone;
858 value_init(mone);
859 value_set_si(mone, -1);
860 value_absolute(d, d);
861 Vector_Scale(coef, val->p, mone, len);
862 value_clear(mone);
863
864 vec_ZZ num;
865 values2zz(val->p, num, len);
866 evalue *EP = multi_monom(num);
867
868 evalue tmp;
869 value_init(tmp.d);
870 value_init(tmp.x.n);
871 value_set_si(tmp.x.n, 1);
872 value_assign(tmp.d, d);
873
874 emul(&tmp, EP);
875
876 ZZ one;
877 one = 1;
878 ceil_mod(val->p, len, d, one, EP);
879
880 /* copy EP to malloc'ed evalue */
881 evalue *E;
882 ALLOC(E);
883 *E = *EP;
884 delete EP;
885
886 free_evalue_refs(&tmp);
887
888 return E;
889 }
890
891 evalue* lattice_point(Polyhedron *i, vec_ZZ& lambda, Matrix *W, Value lcm)
892 {
893 unsigned nparam = W->NbColumns - 1;
894
895 Matrix* Rays = rays2(i);
896 Matrix *T = Transpose(Rays);
897 Matrix *T2 = Matrix_Copy(T);
898 Matrix *inv = Matrix_Alloc(T2->NbRows, T2->NbColumns);
899 int ok = Matrix_Inverse(T2, inv);
900 assert(ok);
901 Matrix_Free(Rays);
902 Matrix_Free(T2);
903 mat_ZZ vertex;
904 matrix2zz(W, vertex, W->NbRows, W->NbColumns);
905
906 vec_ZZ num;
907 num = lambda * vertex;
908
909 evalue *EP = multi_monom(num);
910
911 evalue tmp;
912 value_init(tmp.d);
913 value_init(tmp.x.n);
914 value_set_si(tmp.x.n, 1);
915 value_assign(tmp.d, lcm);
916
917 emul(&tmp, EP);
918
919 Matrix *L = Matrix_Alloc(inv->NbRows, W->NbColumns);
920 Matrix_Product(inv, W, L);
921
922 mat_ZZ RT;
923 matrix2zz(T, RT, T->NbRows, T->NbColumns);
924 Matrix_Free(T);
925
926 vec_ZZ p = lambda * RT;
927
928 for (int i = 0; i < L->NbRows; ++i) {
929 ceil_mod(L->p[i], nparam+1, lcm, p[i], EP);
930 }
931
932 Matrix_Free(L);
933
934 Matrix_Free(inv);
935 free_evalue_refs(&tmp);
936 return EP;
937 }
938 #else
939 evalue* lattice_point(Polyhedron *i, vec_ZZ& lambda, Matrix *W, Value lcm)
940 {
941 Matrix *T = Transpose(W);
942 unsigned nparam = T->NbRows - 1;
943
944 evalue *EP = new evalue();
945 value_init(EP->d);
946 evalue_set_si(EP, 0, 1);
947
948 evalue ev;
949 Vector *val = Vector_Alloc(nparam+1);
950 value_set_si(val->p[nparam], 1);
951 ZZ offset(INIT_VAL, 0);
952 value_init(ev.d);
953 vertex_period(i, lambda, T, lcm, 0, val, EP, &ev, offset);
954 Vector_Free(val);
955 eadd(&ev, EP);
956 free_evalue_refs(&ev);
957
958 Matrix_Free(T);
959
960 reduce_evalue(EP);
961
962 return EP;
963 }
964 #endif
965
966 void lattice_point(
967 Param_Vertices* V, Polyhedron *i, vec_ZZ& lambda, term_info* term)
968 {
969 unsigned nparam = V->Vertex->NbColumns - 2;
970 unsigned dim = i->Dimension;
971 mat_ZZ vertex;
972 vertex.SetDims(V->Vertex->NbRows, nparam+1);
973 Value lcm, tmp;
974 value_init(lcm);
975 value_init(tmp);
976 value_set_si(lcm, 1);
977 for (int j = 0; j < V->Vertex->NbRows; ++j) {
978 value_lcm(lcm, V->Vertex->p[j][nparam+1], &lcm);
979 }
980 if (value_notone_p(lcm)) {
981 Matrix * mv = Matrix_Alloc(dim, nparam+1);
982 for (int j = 0 ; j < dim; ++j) {
983 value_division(tmp, lcm, V->Vertex->p[j][nparam+1]);
984 Vector_Scale(V->Vertex->p[j], mv->p[j], tmp, nparam+1);
985 }
986
987 term->E = lattice_point(i, lambda, mv, lcm);
988 term->constant = 0;
989
990 Matrix_Free(mv);
991 value_clear(lcm);
992 value_clear(tmp);
993 return;
994 }
995 for (int i = 0; i < V->Vertex->NbRows; ++i) {
996 assert(value_one_p(V->Vertex->p[i][nparam+1])); // for now
997 values2zz(V->Vertex->p[i], vertex[i], nparam+1);
998 }
999
1000 vec_ZZ num;
1001 num = lambda * vertex;
1002
1003 int p = -1;
1004 int nn = 0;
1005 for (int j = 0; j < nparam; ++j)
1006 if (num[j] != 0) {
1007 ++nn;
1008 p = j;
1009 }
1010 if (nn >= 2) {
1011 term->E = multi_monom(num);
1012 term->constant = 0;
1013 } else {
1014 term->E = NULL;
1015 term->constant = num[nparam];
1016 term->pos = p;
1017 if (p != -1)
1018 term->coeff = num[p];
1019 }
1020
1021 value_clear(lcm);
1022 value_clear(tmp);
1023 }
1024
1025 void normalize(Polyhedron *i, vec_ZZ& lambda, ZZ& sign, ZZ& num, vec_ZZ& den)
1026 {
1027 unsigned dim = i->Dimension;
1028
1029 int r = 0;
1030 mat_ZZ rays;
1031 rays.SetDims(dim, dim);
1032 add_rays(rays, i, &r);
1033 den = rays * lambda;
1034 int change = 0;
1035
1036 for (int j = 0; j < den.length(); ++j) {
1037 if (den[j] > 0)
1038 change ^= 1;
1039 else {
1040 den[j] = abs(den[j]);
1041 num += den[j];
1042 }
1043 }
1044 if (change)
1045 sign = -sign;
1046 }
1047
1048 void barvinok_count(Polyhedron *P, Value* result, unsigned NbMaxCons)
1049 {
1050 Polyhedron ** vcone;
1051 vec_ZZ sign;
1052 int ncone = 0;
1053 sign.SetLength(ncone);
1054 unsigned dim;
1055 int allocated = 0;
1056 Value factor;
1057 Polyhedron *Q;
1058 int r = 0;
1059
1060 if (emptyQ(P)) {
1061 value_set_si(*result, 0);
1062 return;
1063 }
1064 if (P->NbBid == 0)
1065 for (; r < P->NbRays; ++r)
1066 if (value_zero_p(P->Ray[r][P->Dimension+1]))
1067 break;
1068 if (P->NbBid !=0 || r < P->NbRays) {
1069 value_set_si(*result, -1);
1070 return;
1071 }
1072 if (P->NbEq != 0) {
1073 P = remove_equalities(P);
1074 if (emptyQ(P)) {
1075 Polyhedron_Free(P);
1076 value_set_si(*result, 0);
1077 return;
1078 }
1079 allocated = 1;
1080 }
1081 value_init(factor);
1082 value_set_si(factor, 1);
1083 Q = Polyhedron_Reduce(P, &factor);
1084 if (Q) {
1085 if (allocated)
1086 Polyhedron_Free(P);
1087 P = Q;
1088 allocated = 1;
1089 }
1090 if (P->Dimension == 0) {
1091 value_assign(*result, factor);
1092 if (allocated)
1093 Polyhedron_Free(P);
1094 value_clear(factor);
1095 return;
1096 }
1097
1098 dim = P->Dimension;
1099 vcone = new (Polyhedron *)[P->NbRays];
1100
1101 for (int j = 0; j < P->NbRays; ++j) {
1102 int npos, nneg;
1103 Polyhedron *C = supporting_cone(P, j);
1104 decompose(C, &vcone[j], &npos, &nneg, NbMaxCons);
1105 ncone += npos + nneg;
1106 sign.SetLength(ncone);
1107 for (int k = 0; k < npos; ++k)
1108 sign[ncone-nneg-k-1] = 1;
1109 for (int k = 0; k < nneg; ++k)
1110 sign[ncone-k-1] = -1;
1111 }
1112
1113 mat_ZZ rays;
1114 rays.SetDims(ncone * dim, dim);
1115 r = 0;
1116 for (int j = 0; j < P->NbRays; ++j) {
1117 for (Polyhedron *i = vcone[j]; i; i = i->next) {
1118 assert(i->NbRays-1 == dim);
1119 add_rays(rays, i, &r);
1120 }
1121 }
1122 vec_ZZ lambda;
1123 nonorthog(rays, lambda);
1124
1125 vec_ZZ num;
1126 mat_ZZ den;
1127 num.SetLength(ncone);
1128 den.SetDims(ncone,dim);
1129
1130 int f = 0;
1131 for (int j = 0; j < P->NbRays; ++j) {
1132 for (Polyhedron *i = vcone[j]; i; i = i->next) {
1133 lattice_point(P->Ray[j]+1, i, lambda, num[f]);
1134 normalize(i, lambda, sign[f], num[f], den[f]);
1135 ++f;
1136 }
1137 }
1138 ZZ min = num[0];
1139 for (int j = 1; j < num.length(); ++j)
1140 if (num[j] < min)
1141 min = num[j];
1142 for (int j = 0; j < num.length(); ++j)
1143 num[j] -= min;
1144
1145 f = 0;
1146 mpq_t count;
1147 mpq_init(count);
1148 for (int j = 0; j < P->NbRays; ++j) {
1149 for (Polyhedron *i = vcone[j]; i; i = i->next) {
1150 dpoly d(dim, num[f]);
1151 dpoly n(dim, den[f][0], 1);
1152 for (int k = 1; k < dim; ++k) {
1153 dpoly fact(dim, den[f][k], 1);
1154 n *= fact;
1155 }
1156 d.div(n, count, sign[f]);
1157 ++f;
1158 }
1159 }
1160 assert(value_one_p(&count[0]._mp_den));
1161 value_multiply(*result, &count[0]._mp_num, factor);
1162 mpq_clear(count);
1163
1164 for (int j = 0; j < P->NbRays; ++j)
1165 Domain_Free(vcone[j]);
1166
1167 delete [] vcone;
1168
1169 if (allocated)
1170 Polyhedron_Free(P);
1171 value_clear(factor);
1172 }
1173
1174 static void uni_polynom(int param, Vector *c, evalue *EP)
1175 {
1176 unsigned dim = c->Size-2;
1177 value_init(EP->d);
1178 value_set_si(EP->d,0);
1179 EP->x.p = new_enode(polynomial, dim+1, param+1);
1180 for (int j = 0; j <= dim; ++j)
1181 evalue_set(&EP->x.p->arr[j], c->p[j], c->p[dim+1]);
1182 }
1183
1184 static void multi_polynom(Vector *c, evalue* X, evalue *EP)
1185 {
1186 unsigned dim = c->Size-2;
1187 evalue EC;
1188
1189 value_init(EC.d);
1190 evalue_set(&EC, c->p[dim], c->p[dim+1]);
1191
1192 value_init(EP->d);
1193 evalue_set(EP, c->p[dim], c->p[dim+1]);
1194
1195 for (int i = dim-1; i >= 0; --i) {
1196 emul(X, EP);
1197 value_assign(EC.x.n, c->p[i]);
1198 eadd(&EC, EP);
1199 }
1200 free_evalue_refs(&EC);
1201 }
1202
1203
1204 Enumeration* barvinok_enumerate(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
1205 {
1206 Polyhedron *CEq = NULL, *rVD, *CA;
1207 Matrix *CT = NULL;
1208 Param_Polyhedron *PP = NULL;
1209 Param_Domain *D, *next;
1210 Param_Vertices *V;
1211 Enumeration *en, *res;
1212 int r = 0;
1213 unsigned nparam = C->Dimension;
1214 evalue factor;
1215 value_init(factor.d);
1216 evalue_set_si(&factor, 1, 1);
1217
1218 res = NULL;
1219
1220 CA = align_context(C, P->Dimension, MaxRays);
1221 P = DomainIntersection(P, CA, MaxRays);
1222 Polyhedron_Free(CA);
1223
1224 if (C->Dimension == 0 || emptyQ(P)) {
1225 constant:
1226 res = (Enumeration *)malloc(sizeof(Enumeration));
1227 res->ValidityDomain = CEq ? CEq : Polyhedron_Copy(C);
1228 res->next = NULL;
1229 value_init(res->EP.d);
1230 value_set_si(res->EP.d, 1);
1231 value_init(res->EP.x.n);
1232 if (emptyQ(P))
1233 value_set_si(res->EP.x.n, 0);
1234 else
1235 barvinok_count(P, &res->EP.x.n, MaxRays);
1236 emul(&factor, &res->EP);
1237 out:
1238 free_evalue_refs(&factor);
1239 Polyhedron_Free(P);
1240 if (CT)
1241 Matrix_Free(CT);
1242 if (PP)
1243 Param_Polyhedron_Free(PP);
1244
1245 return res;
1246 }
1247
1248 if (P->NbEq != 0) {
1249 Matrix *f;
1250 P = remove_equalities_p(P, P->Dimension-nparam, &f);
1251 mask(f, &factor);
1252 Matrix_Free(f);
1253 if (P->Dimension == nparam) {
1254 CEq = P;
1255 P = Universe_Polyhedron(0);
1256 goto constant;
1257 }
1258 }
1259
1260 #ifdef USE_MODULO
1261 Polyhedron *Q = ParamPolyhedron_Reduce(P, P->Dimension-nparam, &factor);
1262 if (Q) {
1263 if (Q->Dimension == nparam) {
1264 CEq = Q;
1265 P = Universe_Polyhedron(0);
1266 goto constant;
1267 }
1268 Polyhedron_Free(P);
1269 P = Q;
1270 }
1271 #endif
1272 Polyhedron *oldP = P;
1273 PP = Polyhedron2Param_SimplifiedDomain(&P,C,MaxRays,&CEq,&CT);
1274 if (P != oldP)
1275 Polyhedron_Free(oldP);
1276
1277 if (isIdentity(CT)) {
1278 Matrix_Free(CT);
1279 CT = NULL;
1280 } else {
1281 assert(CT->NbRows != CT->NbColumns);
1282 if (CT->NbRows == 1) // no more parameters
1283 goto constant;
1284 nparam = CT->NbRows - 1;
1285 }
1286
1287 unsigned dim = P->Dimension - nparam;
1288 Polyhedron ** vcone = new (Polyhedron *)[PP->nbV];
1289 int * npos = new int[PP->nbV];
1290 int * nneg = new int[PP->nbV];
1291 vec_ZZ sign;
1292
1293 int i;
1294 for (i = 0, V = PP->V; V; ++i, V = V->next) {
1295 Polyhedron *C = supporting_cone_p(P, V);
1296 decompose(C, &vcone[i], &npos[i], &nneg[i], MaxRays);
1297 }
1298
1299 Vector *c = Vector_Alloc(dim+2);
1300
1301 for(D=PP->D; D; D=next) {
1302 next = D->next;
1303 if (!CEq) {
1304 rVD = D->Domain;
1305 D->Domain = NULL;
1306 } else {
1307 Polyhedron *Dt;
1308 Dt = CT ? Polyhedron_Preimage(D->Domain,CT,MaxRays) : D->Domain;
1309 rVD = DomainIntersection(Dt,CEq,MaxRays);
1310
1311 /* if rVD is empty or too small in geometric dimension */
1312 if(!rVD || emptyQ(rVD) ||
1313 (rVD->Dimension-rVD->NbEq < Dt->Dimension-Dt->NbEq-CEq->NbEq)) {
1314 if(rVD)
1315 Polyhedron_Free(rVD);
1316 if (CT)
1317 Polyhedron_Free(Dt);
1318 continue; /* empty validity domain */
1319 }
1320 if (CT)
1321 Polyhedron_Free(Dt);
1322 }
1323 int ncone = 0;
1324 sign.SetLength(ncone);
1325 FORALL_PVertex_in_ParamPolyhedron(V,D,PP) // _i is internal counter
1326 ncone += npos[_i] + nneg[_i];
1327 sign.SetLength(ncone);
1328 for (int k = 0; k < npos[_i]; ++k)
1329 sign[ncone-nneg[_i]-k-1] = 1;
1330 for (int k = 0; k < nneg[_i]; ++k)
1331 sign[ncone-k-1] = -1;
1332 END_FORALL_PVertex_in_ParamPolyhedron;
1333
1334 mat_ZZ rays;
1335 rays.SetDims(ncone * dim, dim);
1336 r = 0;
1337 FORALL_PVertex_in_ParamPolyhedron(V,D,PP) // _i is internal counter
1338 for (Polyhedron *i = vcone[_i]; i; i = i->next) {
1339 assert(i->NbRays-1 == dim);
1340 add_rays(rays, i, &r);
1341 }
1342 END_FORALL_PVertex_in_ParamPolyhedron;
1343 vec_ZZ lambda;
1344 nonorthog(rays, lambda);
1345
1346 mat_ZZ den;
1347 den.SetDims(ncone,dim);
1348 term_info *num = new term_info[ncone];
1349
1350 int f = 0;
1351 FORALL_PVertex_in_ParamPolyhedron(V,D,PP)
1352 for (Polyhedron *i = vcone[_i]; i; i = i->next) {
1353 lattice_point(V, i, lambda, &num[f]);
1354 normalize(i, lambda, sign[f], num[f].constant, den[f]);
1355 ++f;
1356 }
1357 END_FORALL_PVertex_in_ParamPolyhedron;
1358 ZZ min = num[0].constant;
1359 for (int j = 1; j < ncone; ++j)
1360 if (num[j].constant < min)
1361 min = num[j].constant;
1362 for (int j = 0; j < ncone; ++j)
1363 num[j].constant -= min;
1364 f = 0;
1365 evalue EP;
1366 value_init(EP.d);
1367 evalue_set_si(&EP, 0, 1);
1368 mpq_t count;
1369 mpq_init(count);
1370 FORALL_PVertex_in_ParamPolyhedron(V,D,PP)
1371 for (Polyhedron *i = vcone[_i]; i; i = i->next) {
1372 dpoly n(dim, den[f][0], 1);
1373 for (int k = 1; k < dim; ++k) {
1374 dpoly fact(dim, den[f][k], 1);
1375 n *= fact;
1376 }
1377 if (num[f].E != NULL) {
1378 ZZ one(INIT_VAL, 1);
1379 dpoly_n d(dim, num[f].constant, one);
1380 d.div(n, c, sign[f]);
1381 evalue EV;
1382 multi_polynom(c, num[f].E, &EV);
1383 eadd(&EV , &EP);
1384 free_evalue_refs(&EV);
1385 free_evalue_refs(num[f].E);
1386 delete num[f].E;
1387 } else if (num[f].pos != -1) {
1388 dpoly_n d(dim, num[f].constant, num[f].coeff);
1389 d.div(n, c, sign[f]);
1390 evalue EV;
1391 uni_polynom(num[f].pos, c, &EV);
1392 eadd(&EV , &EP);
1393 free_evalue_refs(&EV);
1394 } else {
1395 mpq_set_si(count, 0, 1);
1396 dpoly d(dim, num[f].constant);
1397 d.div(n, count, sign[f]);
1398 evalue EV;
1399 value_init(EV.d);
1400 evalue_set(&EV, &count[0]._mp_num, &count[0]._mp_den);
1401 eadd(&EV , &EP);
1402 free_evalue_refs(&EV);
1403 }
1404 ++f;
1405 }
1406 END_FORALL_PVertex_in_ParamPolyhedron;
1407
1408 mpq_clear(count);
1409 delete [] num;
1410
1411 en = (Enumeration *)malloc(sizeof(Enumeration));
1412 en->next = res;
1413 res = en;
1414 res->ValidityDomain = rVD;
1415 if (CT)
1416 addeliminatedparams_evalue(&EP, CT);
1417 emul(&factor, &EP);
1418 res->EP = EP;
1419 reduce_evalue(&res->EP);
1420 }
1421
1422 Vector_Free(c);
1423
1424 for (int j = 0; j < PP->nbV; ++j)
1425 Domain_Free(vcone[j]);
1426 delete [] vcone;
1427 delete [] npos;
1428 delete [] nneg;
1429
1430 if (CEq)
1431 Polyhedron_Free(CEq);
1432
1433 goto out;
1434 }
lattice point enumeration library