?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
eliminate a variable if we can (not tested)
[barvinok.git] / barvinok.cc
blob318fb069f606369b056f0e9b167c52d4cf183bd1
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 */
92 static void zz2values(vec_ZZ& v, Value *p)
93 {
94 for (int i = 0; i < v.length(); ++i)
95 zz2value(v[i], p[i]);
96 }
97
98 static void rays(mat_ZZ& r, Polyhedron *C)
99 {
100 unsigned dim = C->NbRays - 1; /* don't count zero vertex */
101 assert(C->NbRays - 1 == C->Dimension);
102 r.SetDims(dim, dim);
103 ZZ tmp;
104
105 int i, c;
106 for (i = 0, c = 0; i < dim; ++i)
107 if (value_zero_p(C->Ray[i][dim+1])) {
108 for (int j = 0; j < dim; ++j) {
109 value2zz(C->Ray[i][j+1], tmp);
110 r[j][c] = tmp;
111 }
112 ++c;
113 }
114 }
115
116 static Matrix * rays(Polyhedron *C)
117 {
118 unsigned dim = C->NbRays - 1; /* don't count zero vertex */
119 assert(C->NbRays - 1 == C->Dimension);
120
121 Matrix *M = Matrix_Alloc(dim+1, dim+1);
122 assert(M);
123
124 int i, c;
125 for (i = 0, c = 0; i <= dim && c < dim; ++i)
126 if (value_zero_p(C->Ray[i][dim+1])) {
127 Vector_Copy(C->Ray[i] + 1, M->p[c], dim);
128 value_set_si(M->p[c++][dim], 0);
129 }
130 assert(c == dim);
131 value_set_si(M->p[dim][dim], 1);
132
133 return M;
134 }
135
136 static Matrix * rays2(Polyhedron *C)
137 {
138 unsigned dim = C->NbRays - 1; /* don't count zero vertex */
139 assert(C->NbRays - 1 == C->Dimension);
140
141 Matrix *M = Matrix_Alloc(dim, dim);
142 assert(M);
143
144 int i, c;
145 for (i = 0, c = 0; i <= dim && c < dim; ++i)
146 if (value_zero_p(C->Ray[i][dim+1]))
147 Vector_Copy(C->Ray[i] + 1, M->p[c++], dim);
148 assert(c == dim);
149
150 return M;
151 }
152
153 /*
154 * Returns the largest absolute value in the vector
155 */
156 static ZZ max(vec_ZZ& v)
157 {
158 ZZ max = abs(v[0]);
159 for (int i = 1; i < v.length(); ++i)
160 if (abs(v[i]) > max)
161 max = abs(v[i]);
162 return max;
163 }
164
165 class cone {
166 public:
167 cone(Matrix *M) {
168 Cone = 0;
169 Rays = Matrix_Copy(M);
170 set_det();
171 }
172 cone(Polyhedron *C) {
173 Cone = Polyhedron_Copy(C);
174 Rays = rays(C);
175 set_det();
176 }
177 void set_det() {
178 mat_ZZ A;
179 matrix2zz(Rays, A, Rays->NbRows - 1, Rays->NbColumns - 1);
180 det = determinant(A);
181 Value v;
182 value_init(v);
183 zz2value(det, v);
184 value_clear(v);
185 }
186
187 Vector* short_vector(vec_ZZ& lambda) {
188 Matrix *M = Matrix_Copy(Rays);
189 Matrix *inv = Matrix_Alloc(M->NbRows, M->NbColumns);
190 int ok = Matrix_Inverse(M, inv);
191 assert(ok);
192 Matrix_Free(M);
193
194 ZZ det2;
195 mat_ZZ B;
196 mat_ZZ U;
197 matrix2zz(inv, B, inv->NbRows - 1, inv->NbColumns - 1);
198 long r = LLL(det2, B, U);
199
200 ZZ min = max(B[0]);
201 int index = 0;
202 for (int i = 1; i < B.NumRows(); ++i) {
203 ZZ tmp = max(B[i]);
204 if (tmp < min) {
205 min = tmp;
206 index = i;
207 }
208 }
209
210 Matrix_Free(inv);
211
212 lambda = B[index];
213
214 Vector *z = Vector_Alloc(U[index].length()+1);
215 assert(z);
216 zz2values(U[index], z->p);
217 value_set_si(z->p[U[index].length()], 0);
218
219 Value tmp;
220 value_init(tmp);
221 Polyhedron *C = poly();
222 int i;
223 for (i = 0; i < C->NbConstraints; ++i) {
224 Inner_Product(z->p, C->Constraint[i]+1, z->Size-1, &tmp);
225 if (value_pos_p(tmp))
226 break;
227 }
228 if (i == C->NbConstraints) {
229 value_set_si(tmp, -1);
230 Vector_Scale(z->p, z->p, tmp, z->Size-1);
231 }
232 value_clear(tmp);
233 return z;
234 }
235
236 ~cone() {
237 Polyhedron_Free(Cone);
238 Matrix_Free(Rays);
239 }
240
241 Polyhedron *poly() {
242 if (!Cone) {
243 Matrix *M = Matrix_Alloc(Rays->NbRows+1, Rays->NbColumns+1);
244 for (int i = 0; i < Rays->NbRows; ++i) {
245 Vector_Copy(Rays->p[i], M->p[i]+1, Rays->NbColumns);
246 value_set_si(M->p[i][0], 1);
247 }
248 Vector_Set(M->p[Rays->NbRows]+1, 0, Rays->NbColumns-1);
249 value_set_si(M->p[Rays->NbRows][0], 1);
250 value_set_si(M->p[Rays->NbRows][Rays->NbColumns], 1);
251 Cone = Rays2Polyhedron(M, M->NbRows+1);
252 assert(Cone->NbConstraints == Cone->NbRays);
253 Matrix_Free(M);
254 }
255 return Cone;
256 }
257
258 ZZ det;
259 Polyhedron *Cone;
260 Matrix *Rays;
261 };
262
263 class dpoly {
264 public:
265 vec_ZZ coeff;
266 dpoly(int d, ZZ& degree, int offset = 0) {
267 coeff.SetLength(d+1);
268
269 int min = d + offset;
270 if (degree < ZZ(INIT_VAL, min))
271 min = to_int(degree);
272
273 ZZ c = ZZ(INIT_VAL, 1);
274 if (!offset)
275 coeff[0] = c;
276 for (int i = 1; i <= min; ++i) {
277 c *= (degree -i + 1);
278 c /= i;
279 coeff[i-offset] = c;
280 }
281 }
282 void operator *= (dpoly& f) {
283 assert(coeff.length() == f.coeff.length());
284 vec_ZZ old = coeff;
285 coeff = f.coeff[0] * coeff;
286 for (int i = 1; i < coeff.length(); ++i)
287 for (int j = 0; i+j < coeff.length(); ++j)
288 coeff[i+j] += f.coeff[i] * old[j];
289 }
290 void div(dpoly& d, mpq_t count, ZZ& sign) {
291 int len = coeff.length();
292 Value tmp;
293 value_init(tmp);
294 mpq_t* c = new mpq_t[coeff.length()];
295 mpq_t qtmp;
296 mpq_init(qtmp);
297 for (int i = 0; i < len; ++i) {
298 mpq_init(c[i]);
299 zz2value(coeff[i], tmp);
300 mpq_set_z(c[i], tmp);
301
302 for (int j = 1; j <= i; ++j) {
303 zz2value(d.coeff[j], tmp);
304 mpq_set_z(qtmp, tmp);
305 mpq_mul(qtmp, qtmp, c[i-j]);
306 mpq_sub(c[i], c[i], qtmp);
307 }
308
309 zz2value(d.coeff[0], tmp);
310 mpq_set_z(qtmp, tmp);
311 mpq_div(c[i], c[i], qtmp);
312 }
313 if (sign == -1)
314 mpq_sub(count, count, c[len-1]);
315 else
316 mpq_add(count, count, c[len-1]);
317
318 value_clear(tmp);
319 mpq_clear(qtmp);
320 for (int i = 0; i < len; ++i)
321 mpq_clear(c[i]);
322 delete [] c;
323 }
324 };
325
326 class dpoly_n {
327 public:
328 Matrix *coeff;
329 ~dpoly_n() {
330 Matrix_Free(coeff);
331 }
332 dpoly_n(int d, ZZ& degree_0, ZZ& degree_1, int offset = 0) {
333 Value d0, d1;
334 value_init(d0);
335 value_init(d1);
336 zz2value(degree_0, d0);
337 zz2value(degree_1, d1);
338 coeff = Matrix_Alloc(d+1, d+1+1);
339 value_set_si(coeff->p[0][0], 1);
340 value_set_si(coeff->p[0][d+1], 1);
341 for (int i = 1; i <= d; ++i) {
342 value_multiply(coeff->p[i][0], coeff->p[i-1][0], d0);
343 Vector_Combine(coeff->p[i-1], coeff->p[i-1]+1, coeff->p[i]+1,
344 d1, d0, i);
345 value_set_si(coeff->p[i][d+1], i);
346 value_multiply(coeff->p[i][d+1], coeff->p[i][d+1], coeff->p[i-1][d+1]);
347 value_decrement(d0, d0);
348 }
349 value_clear(d0);
350 value_clear(d1);
351 }
352 void div(dpoly& d, Vector *count, ZZ& sign) {
353 int len = coeff->NbRows;
354 Matrix * c = Matrix_Alloc(coeff->NbRows, coeff->NbColumns);
355 Value tmp;
356 value_init(tmp);
357 for (int i = 0; i < len; ++i) {
358 Vector_Copy(coeff->p[i], c->p[i], len+1);
359 for (int j = 1; j <= i; ++j) {
360 zz2value(d.coeff[j], tmp);
361 value_multiply(tmp, tmp, c->p[i][len]);
362 value_oppose(tmp, tmp);
363 Vector_Combine(c->p[i], c->p[i-j], c->p[i],
364 c->p[i-j][len], tmp, len);
365 value_multiply(c->p[i][len], c->p[i][len], c->p[i-j][len]);
366 }
367 zz2value(d.coeff[0], tmp);
368 value_multiply(c->p[i][len], c->p[i][len], tmp);
369 }
370 if (sign == -1) {
371 value_set_si(tmp, -1);
372 Vector_Scale(c->p[len-1], count->p, tmp, len);
373 value_assign(count->p[len], c->p[len-1][len]);
374 } else
375 Vector_Copy(c->p[len-1], count->p, len+1);
376 Vector_Normalize(count->p, len+1);
377 value_clear(tmp);
378 Matrix_Free(c);
379 }
380 };
381
382 /*
383 * Barvinok's Decomposition of a simplicial cone
384 *
385 * Returns two lists of polyhedra
386 */
387 void barvinok_decompose(Polyhedron *C, Polyhedron **ppos, Polyhedron **pneg)
388 {
389 Polyhedron *pos = *ppos, *neg = *pneg;
390 vector<cone *> nonuni;
391 cone * c = new cone(C);
392 ZZ det = c->det;
393 int s = sign(det);
394 assert(det != 0);
395 if (abs(det) > 1) {
396 nonuni.push_back(c);
397 } else {
398 Polyhedron *p = Polyhedron_Copy(c->Cone);
399 p->next = pos;
400 pos = p;
401 delete c;
402 }
403 vec_ZZ lambda;
404 while (!nonuni.empty()) {
405 c = nonuni.back();
406 nonuni.pop_back();
407 Vector* v = c->short_vector(lambda);
408 for (int i = 0; i < c->Rays->NbRows - 1; ++i) {
409 if (lambda[i] == 0)
410 continue;
411 Matrix* M = Matrix_Copy(c->Rays);
412 Vector_Copy(v->p, M->p[i], v->Size);
413 cone * pc = new cone(M);
414 assert (pc->det != 0);
415 if (abs(pc->det) > 1) {
416 assert(abs(pc->det) < abs(c->det));
417 nonuni.push_back(pc);
418 } else {
419 Polyhedron *p = pc->poly();
420 pc->Cone = 0;
421 if (sign(pc->det) == s) {
422 p->next = pos;
423 pos = p;
424 } else {
425 p->next = neg;
426 neg = p;
427 }
428 delete pc;
429 }
430 Matrix_Free(M);
431 }
432 Vector_Free(v);
433 delete c;
434 }
435 *ppos = pos;
436 *pneg = neg;
437 }
438
439 /*
440 * Returns a single list of npos "positive" cones followed by nneg
441 * "negative" cones.
442 * The input cone is freed
443 */
444 void decompose(Polyhedron *cone, Polyhedron **parts, int *npos, int *nneg, unsigned MaxRays)
445 {
446 Polyhedron_Polarize(cone);
447 if (cone->NbRays - 1 != cone->Dimension) {
448 Polyhedron *tmp = cone;
449 cone = triangularize_cone(cone, MaxRays);
450 Polyhedron_Free(tmp);
451 }
452 Polyhedron *polpos = NULL, *polneg = NULL;
453 *npos = 0; *nneg = 0;
454 for (Polyhedron *Polar = cone; Polar; Polar = Polar->next)
455 barvinok_decompose(Polar, &polpos, &polneg);
456
457 Polyhedron *last;
458 for (Polyhedron *i = polpos; i; i = i->next) {
459 Polyhedron_Polarize(i);
460 ++*npos;
461 last = i;
462 }
463 for (Polyhedron *i = polneg; i; i = i->next) {
464 Polyhedron_Polarize(i);
465 ++*nneg;
466 }
467 if (last) {
468 last->next = polneg;
469 *parts = polpos;
470 } else
471 *parts = polneg;
472 Domain_Free(cone);
473 }
474
475 const int MAX_TRY=10;
476 /*
477 * Searches for a vector that is not othogonal to any
478 * of the rays in rays.
479 */
480 static void nonorthog(mat_ZZ& rays, vec_ZZ& lambda)
481 {
482 int dim = rays.NumCols();
483 bool found = false;
484 lambda.SetLength(dim);
485 for (int i = 2; !found && i <= 50*dim; i+=4) {
486 for (int j = 0; j < MAX_TRY; ++j) {
487 for (int k = 0; k < dim; ++k) {
488 int r = random_int(i)+2;
489 int v = (2*(r%2)-1) * (r >> 1);
490 lambda[k] = v;
491 }
492 int k = 0;
493 for (; k < rays.NumRows(); ++k)
494 if (lambda * rays[k] == 0)
495 break;
496 if (k == rays.NumRows()) {
497 found = true;
498 break;
499 }
500 }
501 }
502 assert(found);
503 }
504
505 static void add_rays(mat_ZZ& rays, Polyhedron *i, int *r)
506 {
507 unsigned dim = i->Dimension;
508 for (int k = 0; k < i->NbRays; ++k) {
509 if (!value_zero_p(i->Ray[k][dim+1]))
510 continue;
511 values2zz(i->Ray[k]+1, rays[(*r)++], dim);
512 }
513 }
514
515 void lattice_point(Value* values, Polyhedron *i, vec_ZZ& lambda, ZZ& num)
516 {
517 vec_ZZ vertex;
518 unsigned dim = i->Dimension;
519 if(!value_one_p(values[dim])) {
520 Matrix* Rays = rays(i);
521 Matrix *inv = Matrix_Alloc(Rays->NbRows, Rays->NbColumns);
522 int ok = Matrix_Inverse(Rays, inv);
523 assert(ok);
524 Matrix_Free(Rays);
525 Rays = rays(i);
526 Vector *lambda = Vector_Alloc(dim+1);
527 Vector_Matrix_Product(values, inv, lambda->p);
528 Matrix_Free(inv);
529 for (int j = 0; j < dim; ++j)
530 mpz_cdiv_q(lambda->p[j], lambda->p[j], lambda->p[dim]);
531 value_set_si(lambda->p[dim], 1);
532 Vector *A = Vector_Alloc(dim+1);
533 Vector_Matrix_Product(lambda->p, Rays, A->p);
534 Vector_Free(lambda);
535 Matrix_Free(Rays);
536 values2zz(A->p, vertex, dim);
537 Vector_Free(A);
538 } else
539 values2zz(values, vertex, dim);
540
541 num = vertex * lambda;
542 }
543
544 static evalue *term(int param, ZZ& c, Value *den = NULL)
545 {
546 evalue *EP = new evalue();
547 value_init(EP->d);
548 value_set_si(EP->d,0);
549 EP->x.p = new_enode(polynomial, 2, param + 1);
550 evalue_set_si(&EP->x.p->arr[0], 0, 1);
551 value_init(EP->x.p->arr[1].x.n);
552 if (den == NULL)
553 value_set_si(EP->x.p->arr[1].d, 1);
554 else
555 value_assign(EP->x.p->arr[1].d, *den);
556 zz2value(c, EP->x.p->arr[1].x.n);
557 return EP;
558 }
559
560 static void vertex_period(
561 Polyhedron *i, vec_ZZ& lambda, Matrix *T,
562 Value lcm, int p, Vector *val,
563 evalue *E, evalue* ev,
564 ZZ& offset)
565 {
566 unsigned nparam = T->NbRows - 1;
567 unsigned dim = i->Dimension;
568 Value tmp;
569 ZZ nump;
570
571 if (p == nparam) {
572 ZZ num, l;
573 Vector * values = Vector_Alloc(dim + 1);
574 Vector_Matrix_Product(val->p, T, values->p);
575 value_assign(values->p[dim], lcm);
576 lattice_point(values->p, i, lambda, num);
577 value2zz(lcm, l);
578 num *= l;
579 num += offset;
580 value_init(ev->x.n);
581 zz2value(num, ev->x.n);
582 value_assign(ev->d, lcm);
583 Vector_Free(values);
584 return;
585 }
586
587 value_init(tmp);
588 vec_ZZ vertex;
589 values2zz(T->p[p], vertex, dim);
590 nump = vertex * lambda;
591 if (First_Non_Zero(val->p, p) == -1) {
592 value_assign(tmp, lcm);
593 evalue *ET = term(p, nump, &tmp);
594 eadd(ET, E);
595 free_evalue_refs(ET);
596 delete ET;
597 }
598
599 value_assign(tmp, lcm);
600 if (First_Non_Zero(T->p[p], dim) != -1)
601 Vector_Gcd(T->p[p], dim, &tmp);
602 Gcd(tmp, lcm, &tmp);
603 if (value_lt(tmp, lcm)) {
604 ZZ count;
605
606 value_division(tmp, lcm, tmp);
607 value_set_si(ev->d, 0);
608 ev->x.p = new_enode(periodic, VALUE_TO_INT(tmp), p+1);
609 value2zz(tmp, count);
610 do {
611 value_decrement(tmp, tmp);
612 --count;
613 ZZ new_offset = offset - count * nump;
614 value_assign(val->p[p], tmp);
615 vertex_period(i, lambda, T, lcm, p+1, val, E,
616 &ev->x.p->arr[VALUE_TO_INT(tmp)], new_offset);
617 } while (value_pos_p(tmp));
618 } else
619 vertex_period(i, lambda, T, lcm, p+1, val, E, ev, offset);
620 value_clear(tmp);
621 }
622
623 static void mask_r(Matrix *f, int nr, Vector *lcm, int p, Vector *val, evalue *ev)
624 {
625 unsigned nparam = lcm->Size;
626
627 if (p == nparam) {
628 Vector * prod = Vector_Alloc(f->NbRows);
629 Matrix_Vector_Product(f, val->p, prod->p);
630 int isint = 1;
631 for (int i = 0; i < nr; ++i) {
632 value_modulus(prod->p[i], prod->p[i], f->p[i][nparam+1]);
633 isint &= value_zero_p(prod->p[i]);
634 }
635 value_set_si(ev->d, 1);
636 value_init(ev->x.n);
637 value_set_si(ev->x.n, isint);
638 Vector_Free(prod);
639 return;
640 }
641
642 Value tmp;
643 value_init(tmp);
644 if (value_one_p(lcm->p[p]))
645 mask_r(f, nr, lcm, p+1, val, ev);
646 else {
647 value_assign(tmp, lcm->p[p]);
648 value_set_si(ev->d, 0);
649 ev->x.p = new_enode(periodic, VALUE_TO_INT(tmp), p+1);
650 do {
651 value_decrement(tmp, tmp);
652 value_assign(val->p[p], tmp);
653 mask_r(f, nr, lcm, p+1, val, &ev->x.p->arr[VALUE_TO_INT(tmp)]);
654 } while (value_pos_p(tmp));
655 }
656 value_clear(tmp);
657 }
658
659 static evalue *multi_monom(vec_ZZ& p)
660 {
661 evalue *X = new evalue();
662 value_init(X->d);
663 value_init(X->x.n);
664 unsigned nparam = p.length()-1;
665 zz2value(p[nparam], X->x.n);
666 value_set_si(X->d, 1);
667 for (int i = 0; i < nparam; ++i) {
668 if (p[i] == 0)
669 continue;
670 evalue *T = term(i, p[i]);
671 eadd(T, X);
672 free_evalue_refs(T);
673 delete T;
674 }
675 return X;
676 }
677
678 #ifdef USE_MODULO
679 static void mask(Matrix *f, evalue *factor)
680 {
681 int nr = f->NbRows, nc = f->NbColumns;
682 int n;
683 bool found = false;
684 for (n = 0; n < nr && value_notzero_p(f->p[n][nc-1]); ++n)
685 if (value_notone_p(f->p[n][nc-1]) &&
686 value_notmone_p(f->p[n][nc-1]))
687 found = true;
688 if (!found)
689 return;
690
691 evalue EP;
692 nr = n;
693
694 for (n = 0; n < nr; ++n) {
695 if (value_one_p(f->p[n][nc-1]) ||
696 value_mone_p(f->p[n][nc-1]))
697 continue;
698 value_init(EP.d);
699 value_set_si(EP.d, 0);
700 EP.x.p = new_enode(relation, 2, 0);
701 value_clear(EP.x.p->arr[1].d);
702 EP.x.p->arr[1] = *factor;
703 evalue *ev = &EP.x.p->arr[0];
704 value_set_si(ev->d, 0);
705 ev->x.p = new_enode(modulo, 3, VALUE_TO_INT(f->p[n][nc-1]));
706 evalue_set_si(&ev->x.p->arr[1], 0, 1);
707 evalue_set_si(&ev->x.p->arr[2], 1, 1);
708 vec_ZZ row;
709 values2zz(f->p[n], row, nc-1);
710 evalue *E = multi_monom(row);
711 value_clear(ev->x.p->arr[0].d);
712 ev->x.p->arr[0] = *E;
713 delete E;
714 *factor = EP;
715 }
716 }
717 #else
718 /*
719 *
720 */
721 static void mask(Matrix *f, evalue *factor)
722 {
723 int nr = f->NbRows, nc = f->NbColumns;
724 int n;
725 bool found = false;
726 for (n = 0; n < nr && value_notzero_p(f->p[n][nc-1]); ++n)
727 if (value_notone_p(f->p[n][nc-1]) &&
728 value_notmone_p(f->p[n][nc-1]))
729 found = true;
730 if (!found)
731 return;
732
733 Value tmp;
734 value_init(tmp);
735 nr = n;
736 unsigned np = nc - 2;
737 Vector *lcm = Vector_Alloc(np);
738 Vector *val = Vector_Alloc(nc);
739 Vector_Set(val->p, 0, nc);
740 value_set_si(val->p[np], 1);
741 Vector_Set(lcm->p, 1, np);
742 for (n = 0; n < nr; ++n) {
743 if (value_one_p(f->p[n][nc-1]) ||
744 value_mone_p(f->p[n][nc-1]))
745 continue;
746 for (int j = 0; j < np; ++j)
747 if (value_notzero_p(f->p[n][j])) {
748 Gcd(f->p[n][j], f->p[n][nc-1], &tmp);
749 value_division(tmp, f->p[n][nc-1], tmp);
750 value_lcm(tmp, lcm->p[j], &lcm->p[j]);
751 }
752 }
753 evalue EP;
754 value_init(EP.d);
755 mask_r(f, nr, lcm, 0, val, &EP);
756 value_clear(tmp);
757 Vector_Free(val);
758 Vector_Free(lcm);
759 emul(&EP,factor);
760 free_evalue_refs(&EP);
761 }
762 #endif
763
764 struct term_info {
765 evalue *E;
766 ZZ constant;
767 ZZ coeff;
768 int pos;
769 };
770
771 #ifdef USE_MODULO
772 bool mod_needed(Polyhedron *PD, vec_ZZ& num, Value d, evalue *E)
773 {
774 bool ret = true;
775 int len = num.length();
776 Matrix *T = Matrix_Alloc(2, len);
777 zz2values(num, T->p[0]);
778 value_set_si(T->p[1][len-1], 1);
779 Polyhedron *I = Polyhedron_Image(PD, T, PD->NbConstraints);
780 Matrix_Free(T);
781
782 int i;
783 for (i = 0; i < I->NbRays; ++i)
784 if (value_zero_p(I->Ray[i][2])) {
785 Polyhedron_Free(I);
786 return true;
787 }
788
789 Value min, max;
790 value_init(min);
791 value_init(max);
792 for (i = 0; i < I->NbConstraints; ++i) {
793 value_oppose(I->Constraint[i][2], I->Constraint[i][2]);
794 if (value_pos_p(I->Constraint[i][1]))
795 mpz_cdiv_q(min, I->Constraint[i][2], I->Constraint[i][1]);
796 else
797 mpz_fdiv_q(max, I->Constraint[i][2], I->Constraint[i][1]);
798 }
799 mpz_fdiv_q(min, min, d);
800 mpz_fdiv_q(max, max, d);
801 if (value_eq(min, max)) {
802 ret = false;
803 value_oppose(min, min);
804 evalue EV;
805 value_init(EV.d);
806 value_set_si(EV.d, 1);
807 value_init(EV.x.n);
808 value_multiply(EV.x.n, min, d);
809 eadd(&EV, E);
810 free_evalue_refs(&EV);
811 }
812 value_clear(min);
813 value_clear(max);
814 Polyhedron_Free(I);
815 return ret;
816 }
817
818 void ceil_mod(Value *coef, int len, Value d, ZZ& f, evalue *EP, Polyhedron *PD)
819 {
820 Value gcd;
821 value_init(gcd);
822 Value mone;
823 value_init(mone);
824 value_set_si(mone, -1);
825
826 vec_ZZ num;
827
828 Vector_Gcd(coef, len, &gcd);
829 Gcd(gcd, d, &gcd);
830 Vector_AntiScale(coef, coef, gcd, len);
831 Vector_Scale(coef, coef, mone, len);
832
833 value_division(gcd, d, gcd);
834 ZZ g;
835 value2zz(gcd, g);
836 if (value_one_p(gcd))
837 goto out;
838
839 values2zz(coef, num, len);
840
841 int j;
842 for (j = 0; j < len; ++j)
843 num[j] = num[j] % g;
844 for (j = 0; j < len; ++j)
845 if (num[j] != 0)
846 break;
847 if (j == len)
848 goto out;
849
850 evalue tmp;
851 value_init(tmp.d);
852 evalue_set_si(&tmp, 0, 1);
853
854 if (j < len-1 && num[j] > g/2) {
855 for (int k = j; k < len-1; ++k)
856 if (num[k] != 0)
857 num[k] = g - num[k];
858 num[len-1] = g - 1 - num[len-1];
859 value_assign(tmp.d, gcd);
860 ZZ t = f*(g-1);
861 zz2value(t, tmp.x.n);
862 eadd(&tmp, EP);
863 f = -f;
864 }
865
866 if (j >= len-1) {
867 ZZ t = num[len-1] * f;
868 zz2value(t, tmp.x.n);
869 value_assign(tmp.d, gcd);
870 eadd(&tmp, EP);
871 } else {
872 evalue *E = multi_monom(num);
873 evalue EV;
874 value_init(EV.d);
875
876 if (PD && !mod_needed(PD, num, gcd, E)) {
877 value_init(EV.x.n);
878 zz2value(f, EV.x.n);
879 value_assign(EV.d, gcd);
880 emul(&EV, E);
881 eadd(E, EP);
882 } else {
883 value_set_si(EV.d, 0);
884 EV.x.p = new_enode(modulo, 3, VALUE_TO_INT(gcd));
885 evalue_copy(&EV.x.p->arr[0], E);
886 evalue_set_si(&EV.x.p->arr[1], 0, 1);
887 value_init(EV.x.p->arr[2].x.n);
888 zz2value(f, EV.x.p->arr[2].x.n);
889 value_assign(EV.x.p->arr[2].d, gcd);
890
891 eadd(&EV, EP);
892 }
893
894 free_evalue_refs(&EV);
895 free_evalue_refs(E);
896 delete E;
897 }
898
899 free_evalue_refs(&tmp);
900
901 out:
902 value_clear(gcd);
903 value_clear(mone);
904 }
905
906 evalue* ceil3(Value *coef, int len, Value d)
907 {
908 Vector *val = Vector_Alloc(len);
909
910 Value mone;
911 value_init(mone);
912 value_set_si(mone, -1);
913 value_absolute(d, d);
914 Vector_Scale(coef, val->p, mone, len);
915 value_clear(mone);
916
917 vec_ZZ num;
918 values2zz(val->p, num, len);
919 evalue *EP = multi_monom(num);
920
921 evalue tmp;
922 value_init(tmp.d);
923 value_init(tmp.x.n);
924 value_set_si(tmp.x.n, 1);
925 value_assign(tmp.d, d);
926
927 emul(&tmp, EP);
928
929 ZZ one;
930 one = 1;
931 ceil_mod(val->p, len, d, one, EP, NULL);
932
933 /* copy EP to malloc'ed evalue */
934 evalue *E;
935 ALLOC(E);
936 *E = *EP;
937 delete EP;
938
939 free_evalue_refs(&tmp);
940
941 return E;
942 }
943
944 evalue* lattice_point(
945 Polyhedron *i, vec_ZZ& lambda, Matrix *W, Value lcm, Polyhedron *PD)
946 {
947 unsigned nparam = W->NbColumns - 1;
948
949 Matrix* Rays = rays2(i);
950 Matrix *T = Transpose(Rays);
951 Matrix *T2 = Matrix_Copy(T);
952 Matrix *inv = Matrix_Alloc(T2->NbRows, T2->NbColumns);
953 int ok = Matrix_Inverse(T2, inv);
954 assert(ok);
955 Matrix_Free(Rays);
956 Matrix_Free(T2);
957 mat_ZZ vertex;
958 matrix2zz(W, vertex, W->NbRows, W->NbColumns);
959
960 vec_ZZ num;
961 num = lambda * vertex;
962
963 evalue *EP = multi_monom(num);
964
965 evalue tmp;
966 value_init(tmp.d);
967 value_init(tmp.x.n);
968 value_set_si(tmp.x.n, 1);
969 value_assign(tmp.d, lcm);
970
971 emul(&tmp, EP);
972
973 Matrix *L = Matrix_Alloc(inv->NbRows, W->NbColumns);
974 Matrix_Product(inv, W, L);
975
976 mat_ZZ RT;
977 matrix2zz(T, RT, T->NbRows, T->NbColumns);
978 Matrix_Free(T);
979
980 vec_ZZ p = lambda * RT;
981
982 for (int i = 0; i < L->NbRows; ++i) {
983 ceil_mod(L->p[i], nparam+1, lcm, p[i], EP, PD);
984 }
985
986 Matrix_Free(L);
987
988 Matrix_Free(inv);
989 free_evalue_refs(&tmp);
990 return EP;
991 }
992 #else
993 evalue* lattice_point(
994 Polyhedron *i, vec_ZZ& lambda, Matrix *W, Value lcm, Polyhedron *PD)
995 {
996 Matrix *T = Transpose(W);
997 unsigned nparam = T->NbRows - 1;
998
999 evalue *EP = new evalue();
1000 value_init(EP->d);
1001 evalue_set_si(EP, 0, 1);
1002
1003 evalue ev;
1004 Vector *val = Vector_Alloc(nparam+1);
1005 value_set_si(val->p[nparam], 1);
1006 ZZ offset(INIT_VAL, 0);
1007 value_init(ev.d);
1008 vertex_period(i, lambda, T, lcm, 0, val, EP, &ev, offset);
1009 Vector_Free(val);
1010 eadd(&ev, EP);
1011 free_evalue_refs(&ev);
1012
1013 Matrix_Free(T);
1014
1015 reduce_evalue(EP);
1016
1017 return EP;
1018 }
1019 #endif
1020
1021 void lattice_point(
1022 Param_Vertices* V, Polyhedron *i, vec_ZZ& lambda, term_info* term,
1023 Polyhedron *PD)
1024 {
1025 unsigned nparam = V->Vertex->NbColumns - 2;
1026 unsigned dim = i->Dimension;
1027 mat_ZZ vertex;
1028 vertex.SetDims(V->Vertex->NbRows, nparam+1);
1029 Value lcm, tmp;
1030 value_init(lcm);
1031 value_init(tmp);
1032 value_set_si(lcm, 1);
1033 for (int j = 0; j < V->Vertex->NbRows; ++j) {
1034 value_lcm(lcm, V->Vertex->p[j][nparam+1], &lcm);
1035 }
1036 if (value_notone_p(lcm)) {
1037 Matrix * mv = Matrix_Alloc(dim, nparam+1);
1038 for (int j = 0 ; j < dim; ++j) {
1039 value_division(tmp, lcm, V->Vertex->p[j][nparam+1]);
1040 Vector_Scale(V->Vertex->p[j], mv->p[j], tmp, nparam+1);
1041 }
1042
1043 term->E = lattice_point(i, lambda, mv, lcm, PD);
1044 term->constant = 0;
1045
1046 Matrix_Free(mv);
1047 value_clear(lcm);
1048 value_clear(tmp);
1049 return;
1050 }
1051 for (int i = 0; i < V->Vertex->NbRows; ++i) {
1052 assert(value_one_p(V->Vertex->p[i][nparam+1])); // for now
1053 values2zz(V->Vertex->p[i], vertex[i], nparam+1);
1054 }
1055
1056 vec_ZZ num;
1057 num = lambda * vertex;
1058
1059 int p = -1;
1060 int nn = 0;
1061 for (int j = 0; j < nparam; ++j)
1062 if (num[j] != 0) {
1063 ++nn;
1064 p = j;
1065 }
1066 if (nn >= 2) {
1067 term->E = multi_monom(num);
1068 term->constant = 0;
1069 } else {
1070 term->E = NULL;
1071 term->constant = num[nparam];
1072 term->pos = p;
1073 if (p != -1)
1074 term->coeff = num[p];
1075 }
1076
1077 value_clear(lcm);
1078 value_clear(tmp);
1079 }
1080
1081 void normalize(Polyhedron *i, vec_ZZ& lambda, ZZ& sign, ZZ& num, vec_ZZ& den)
1082 {
1083 unsigned dim = i->Dimension;
1084
1085 int r = 0;
1086 mat_ZZ rays;
1087 rays.SetDims(dim, dim);
1088 add_rays(rays, i, &r);
1089 den = rays * lambda;
1090 int change = 0;
1091
1092 for (int j = 0; j < den.length(); ++j) {
1093 if (den[j] > 0)
1094 change ^= 1;
1095 else {
1096 den[j] = abs(den[j]);
1097 num += den[j];
1098 }
1099 }
1100 if (change)
1101 sign = -sign;
1102 }
1103
1104 void barvinok_count(Polyhedron *P, Value* result, unsigned NbMaxCons)
1105 {
1106 Polyhedron ** vcone;
1107 vec_ZZ sign;
1108 int ncone = 0;
1109 sign.SetLength(ncone);
1110 unsigned dim;
1111 int allocated = 0;
1112 Value factor;
1113 Polyhedron *Q;
1114 int r = 0;
1115
1116 if (emptyQ(P)) {
1117 value_set_si(*result, 0);
1118 return;
1119 }
1120 if (P->NbBid == 0)
1121 for (; r < P->NbRays; ++r)
1122 if (value_zero_p(P->Ray[r][P->Dimension+1]))
1123 break;
1124 if (P->NbBid !=0 || r < P->NbRays) {
1125 value_set_si(*result, -1);
1126 return;
1127 }
1128 if (P->NbEq != 0) {
1129 P = remove_equalities(P);
1130 if (emptyQ(P)) {
1131 Polyhedron_Free(P);
1132 value_set_si(*result, 0);
1133 return;
1134 }
1135 allocated = 1;
1136 }
1137 value_init(factor);
1138 value_set_si(factor, 1);
1139 Q = Polyhedron_Reduce(P, &factor);
1140 if (Q) {
1141 if (allocated)
1142 Polyhedron_Free(P);
1143 P = Q;
1144 allocated = 1;
1145 }
1146 if (P->Dimension == 0) {
1147 value_assign(*result, factor);
1148 if (allocated)
1149 Polyhedron_Free(P);
1150 value_clear(factor);
1151 return;
1152 }
1153
1154 dim = P->Dimension;
1155 vcone = new (Polyhedron *)[P->NbRays];
1156
1157 for (int j = 0; j < P->NbRays; ++j) {
1158 int npos, nneg;
1159 Polyhedron *C = supporting_cone(P, j);
1160 decompose(C, &vcone[j], &npos, &nneg, NbMaxCons);
1161 ncone += npos + nneg;
1162 sign.SetLength(ncone);
1163 for (int k = 0; k < npos; ++k)
1164 sign[ncone-nneg-k-1] = 1;
1165 for (int k = 0; k < nneg; ++k)
1166 sign[ncone-k-1] = -1;
1167 }
1168
1169 mat_ZZ rays;
1170 rays.SetDims(ncone * dim, dim);
1171 r = 0;
1172 for (int j = 0; j < P->NbRays; ++j) {
1173 for (Polyhedron *i = vcone[j]; i; i = i->next) {
1174 assert(i->NbRays-1 == dim);
1175 add_rays(rays, i, &r);
1176 }
1177 }
1178 vec_ZZ lambda;
1179 nonorthog(rays, lambda);
1180
1181 vec_ZZ num;
1182 mat_ZZ den;
1183 num.SetLength(ncone);
1184 den.SetDims(ncone,dim);
1185
1186 int f = 0;
1187 for (int j = 0; j < P->NbRays; ++j) {
1188 for (Polyhedron *i = vcone[j]; i; i = i->next) {
1189 lattice_point(P->Ray[j]+1, i, lambda, num[f]);
1190 normalize(i, lambda, sign[f], num[f], den[f]);
1191 ++f;
1192 }
1193 }
1194 ZZ min = num[0];
1195 for (int j = 1; j < num.length(); ++j)
1196 if (num[j] < min)
1197 min = num[j];
1198 for (int j = 0; j < num.length(); ++j)
1199 num[j] -= min;
1200
1201 f = 0;
1202 mpq_t count;
1203 mpq_init(count);
1204 for (int j = 0; j < P->NbRays; ++j) {
1205 for (Polyhedron *i = vcone[j]; i; i = i->next) {
1206 dpoly d(dim, num[f]);
1207 dpoly n(dim, den[f][0], 1);
1208 for (int k = 1; k < dim; ++k) {
1209 dpoly fact(dim, den[f][k], 1);
1210 n *= fact;
1211 }
1212 d.div(n, count, sign[f]);
1213 ++f;
1214 }
1215 }
1216 assert(value_one_p(&count[0]._mp_den));
1217 value_multiply(*result, &count[0]._mp_num, factor);
1218 mpq_clear(count);
1219
1220 for (int j = 0; j < P->NbRays; ++j)
1221 Domain_Free(vcone[j]);
1222
1223 delete [] vcone;
1224
1225 if (allocated)
1226 Polyhedron_Free(P);
1227 value_clear(factor);
1228 }
1229
1230 static void uni_polynom(int param, Vector *c, evalue *EP)
1231 {
1232 unsigned dim = c->Size-2;
1233 value_init(EP->d);
1234 value_set_si(EP->d,0);
1235 EP->x.p = new_enode(polynomial, dim+1, param+1);
1236 for (int j = 0; j <= dim; ++j)
1237 evalue_set(&EP->x.p->arr[j], c->p[j], c->p[dim+1]);
1238 }
1239
1240 static void multi_polynom(Vector *c, evalue* X, evalue *EP)
1241 {
1242 unsigned dim = c->Size-2;
1243 evalue EC;
1244
1245 value_init(EC.d);
1246 evalue_set(&EC, c->p[dim], c->p[dim+1]);
1247
1248 value_init(EP->d);
1249 evalue_set(EP, c->p[dim], c->p[dim+1]);
1250
1251 for (int i = dim-1; i >= 0; --i) {
1252 emul(X, EP);
1253 value_assign(EC.x.n, c->p[i]);
1254 eadd(&EC, EP);
1255 }
1256 free_evalue_refs(&EC);
1257 }
1258
1259
1260 evalue* barvinok_enumerate_ev(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
1261 {
1262 Polyhedron *CEq = NULL, *rVD, *pVD, *CA;
1263 Matrix *CT = NULL;
1264 Param_Polyhedron *PP = NULL;
1265 Param_Domain *D, *next;
1266 Param_Vertices *V;
1267 int r = 0;
1268 unsigned nparam = C->Dimension;
1269 evalue *eres = new evalue;
1270 value_init(eres->d);
1271 value_set_si(eres->d, 0);
1272
1273 evalue factor;
1274 value_init(factor.d);
1275 evalue_set_si(&factor, 1, 1);
1276
1277 CA = align_context(C, P->Dimension, MaxRays);
1278 P = DomainIntersection(P, CA, MaxRays);
1279 Polyhedron_Free(CA);
1280
1281 if (C->Dimension == 0 || emptyQ(P)) {
1282 constant:
1283 eres->x.p = new_enode(partition, 2, -1);
1284 EVALUE_SET_DOMAIN(eres->x.p->arr[0], CEq ? CEq : Polyhedron_Copy(C));
1285 value_set_si(eres->x.p->arr[1].d, 1);
1286 value_init(eres->x.p->arr[1].x.n);
1287 if (emptyQ(P))
1288 value_set_si(eres->x.p->arr[1].x.n, 0);
1289 else
1290 barvinok_count(P, &eres->x.p->arr[1].x.n, MaxRays);
1291 emul(&factor, &eres->x.p->arr[1]);
1292 out:
1293 free_evalue_refs(&factor);
1294 Polyhedron_Free(P);
1295 if (CT)
1296 Matrix_Free(CT);
1297 if (PP)
1298 Param_Polyhedron_Free(PP);
1299
1300 return eres;
1301 }
1302 for (r = 0; r < P->NbRays; ++r)
1303 if (value_zero_p(P->Ray[r][0]) ||
1304 value_zero_p(P->Ray[r][P->Dimension+1])) {
1305 int i;
1306 for (i = P->Dimension - nparam; i < P->Dimension; ++i)
1307 if (value_notzero_p(P->Ray[r][i+1]))
1308 break;
1309 if (i >= P->Dimension)
1310 break;
1311 }
1312 if (r < P->NbRays)
1313 goto constant;
1314
1315 if (P->NbEq != 0) {
1316 Matrix *f;
1317 P = remove_equalities_p(P, P->Dimension-nparam, &f);
1318 mask(f, &factor);
1319 Matrix_Free(f);
1320 if (P->Dimension == nparam) {
1321 CEq = P;
1322 P = Universe_Polyhedron(0);
1323 goto constant;
1324 }
1325 }
1326
1327 #ifdef USE_MODULO
1328 Polyhedron *Q = ParamPolyhedron_Reduce(P, P->Dimension-nparam, &factor);
1329 if (Q) {
1330 if (Q->Dimension == nparam) {
1331 CEq = Q;
1332 P = Universe_Polyhedron(0);
1333 goto constant;
1334 }
1335 Polyhedron_Free(P);
1336 P = Q;
1337 }
1338 #endif
1339 Polyhedron *oldP = P;
1340 PP = Polyhedron2Param_SimplifiedDomain(&P,C,MaxRays,&CEq,&CT);
1341 if (P != oldP)
1342 Polyhedron_Free(oldP);
1343
1344 if (isIdentity(CT)) {
1345 Matrix_Free(CT);
1346 CT = NULL;
1347 } else {
1348 assert(CT->NbRows != CT->NbColumns);
1349 if (CT->NbRows == 1) // no more parameters
1350 goto constant;
1351 nparam = CT->NbRows - 1;
1352 }
1353
1354 unsigned dim = P->Dimension - nparam;
1355 Polyhedron ** vcone = new (Polyhedron *)[PP->nbV];
1356 int * npos = new int[PP->nbV];
1357 int * nneg = new int[PP->nbV];
1358 vec_ZZ sign;
1359
1360 int i;
1361 for (i = 0, V = PP->V; V; ++i, V = V->next) {
1362 Polyhedron *C = supporting_cone_p(P, V);
1363 decompose(C, &vcone[i], &npos[i], &nneg[i], MaxRays);
1364 }
1365
1366 Vector *c = Vector_Alloc(dim+2);
1367
1368 int nd;
1369 for (nd = 0, D=PP->D; D; ++nd, D=D->next);
1370 struct section { Polyhedron * D; evalue E; };
1371 section *s = new section[nd];
1372
1373 for(nd = 0, D=PP->D; D; D=next) {
1374 next = D->next;
1375 if (!CEq) {
1376 pVD = rVD = D->Domain;
1377 D->Domain = NULL;
1378 } else {
1379 Polyhedron *Dt;
1380 Dt = CT ? DomainPreimage(D->Domain,CT,MaxRays) : D->Domain;
1381 rVD = DomainIntersection(Dt,CEq,MaxRays);
1382
1383 /* if rVD is empty or too small in geometric dimension */
1384 if(!rVD || emptyQ(rVD) ||
1385 (rVD->Dimension-rVD->NbEq < Dt->Dimension-Dt->NbEq-CEq->NbEq)) {
1386 if(rVD)
1387 Domain_Free(rVD);
1388 if (CT)
1389 Domain_Free(Dt);
1390 continue; /* empty validity domain */
1391 }
1392 if (CT)
1393 Domain_Free(Dt);
1394 pVD = CT ? DomainImage(rVD,CT,MaxRays) : rVD;
1395 for (Param_Domain *D2 = D->next; D2; D2=D2->next) {
1396 Polyhedron *T = D2->Domain;
1397 D2->Domain = DomainDifference(D2->Domain, D->Domain, MaxRays);
1398 Domain_Free(T);
1399 }
1400 }
1401 int ncone = 0;
1402 sign.SetLength(ncone);
1403 FORALL_PVertex_in_ParamPolyhedron(V,D,PP) // _i is internal counter
1404 ncone += npos[_i] + nneg[_i];
1405 sign.SetLength(ncone);
1406 for (int k = 0; k < npos[_i]; ++k)
1407 sign[ncone-nneg[_i]-k-1] = 1;
1408 for (int k = 0; k < nneg[_i]; ++k)
1409 sign[ncone-k-1] = -1;
1410 END_FORALL_PVertex_in_ParamPolyhedron;
1411
1412 mat_ZZ rays;
1413 rays.SetDims(ncone * dim, dim);
1414 r = 0;
1415 FORALL_PVertex_in_ParamPolyhedron(V,D,PP) // _i is internal counter
1416 for (Polyhedron *i = vcone[_i]; i; i = i->next) {
1417 assert(i->NbRays-1 == dim);
1418 add_rays(rays, i, &r);
1419 }
1420 END_FORALL_PVertex_in_ParamPolyhedron;
1421 vec_ZZ lambda;
1422 nonorthog(rays, lambda);
1423
1424 mat_ZZ den;
1425 den.SetDims(ncone,dim);
1426 term_info *num = new term_info[ncone];
1427
1428 int f = 0;
1429 FORALL_PVertex_in_ParamPolyhedron(V,D,PP)
1430 for (Polyhedron *i = vcone[_i]; i; i = i->next) {
1431 lattice_point(V, i, lambda, &num[f], pVD);
1432 normalize(i, lambda, sign[f], num[f].constant, den[f]);
1433 ++f;
1434 }
1435 END_FORALL_PVertex_in_ParamPolyhedron;
1436 ZZ min = num[0].constant;
1437 for (int j = 1; j < ncone; ++j)
1438 if (num[j].constant < min)
1439 min = num[j].constant;
1440 for (int j = 0; j < ncone; ++j)
1441 num[j].constant -= min;
1442 f = 0;
1443 value_init(s[nd].E.d);
1444 evalue_set_si(&s[nd].E, 0, 1);
1445 mpq_t count;
1446 mpq_init(count);
1447 FORALL_PVertex_in_ParamPolyhedron(V,D,PP)
1448 for (Polyhedron *i = vcone[_i]; i; i = i->next) {
1449 dpoly n(dim, den[f][0], 1);
1450 for (int k = 1; k < dim; ++k) {
1451 dpoly fact(dim, den[f][k], 1);
1452 n *= fact;
1453 }
1454 if (num[f].E != NULL) {
1455 ZZ one(INIT_VAL, 1);
1456 dpoly_n d(dim, num[f].constant, one);
1457 d.div(n, c, sign[f]);
1458 evalue EV;
1459 multi_polynom(c, num[f].E, &EV);
1460 eadd(&EV , &s[nd].E);
1461 free_evalue_refs(&EV);
1462 free_evalue_refs(num[f].E);
1463 delete num[f].E;
1464 } else if (num[f].pos != -1) {
1465 dpoly_n d(dim, num[f].constant, num[f].coeff);
1466 d.div(n, c, sign[f]);
1467 evalue EV;
1468 uni_polynom(num[f].pos, c, &EV);
1469 eadd(&EV , &s[nd].E);
1470 free_evalue_refs(&EV);
1471 } else {
1472 mpq_set_si(count, 0, 1);
1473 dpoly d(dim, num[f].constant);
1474 d.div(n, count, sign[f]);
1475 evalue EV;
1476 value_init(EV.d);
1477 evalue_set(&EV, &count[0]._mp_num, &count[0]._mp_den);
1478 eadd(&EV , &s[nd].E);
1479 free_evalue_refs(&EV);
1480 }
1481 ++f;
1482 }
1483 END_FORALL_PVertex_in_ParamPolyhedron;
1484
1485 mpq_clear(count);
1486 delete [] num;
1487
1488 if (CT)
1489 addeliminatedparams_evalue(&s[nd].E, CT);
1490 emul(&factor, &s[nd].E);
1491 reduce_evalue(&s[nd].E);
1492 s[nd].D = rVD;
1493 ++nd;
1494 if (rVD != pVD)
1495 Domain_Free(pVD);
1496 }
1497
1498 eres->x.p = new_enode(partition, 2*nd, -1);
1499 for (int j = 0; j < nd; ++j) {
1500 EVALUE_SET_DOMAIN(eres->x.p->arr[2*j], s[j].D);
1501 eres->x.p->arr[2*j+1] = s[j].E;
1502 }
1503 delete [] s;
1504
1505 Vector_Free(c);
1506
1507 for (int j = 0; j < PP->nbV; ++j)
1508 Domain_Free(vcone[j]);
1509 delete [] vcone;
1510 delete [] npos;
1511 delete [] nneg;
1512
1513 if (CEq)
1514 Polyhedron_Free(CEq);
1515
1516 goto out;
1517 }
1518
1519 Enumeration* barvinok_enumerate(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
1520 {
1521 Enumeration *en, *res = NULL;
1522 evalue *EP = barvinok_enumerate_ev(P, C, MaxRays);
1523 for (int i = 0; i < EP->x.p->size/2; ++i) {
1524 en = (Enumeration *)malloc(sizeof(Enumeration));
1525 en->next = res;
1526 res = en;
1527 res->ValidityDomain = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
1528 res->EP = EP->x.p->arr[2*i+1];
1529 }
1530 free(EP->x.p);
1531 value_clear(EP->d);
1532 delete EP;
1533 return res;
1534 }
1535
1536 static void SwapColumns(Value **V, int n, int i, int j)
1537 {
1538 for (int r = 0; r < n; ++r)
1539 value_swap(V[r][i], V[r][j]);
1540 }
1541
1542 static void SwapColumns(Polyhedron *P, int i, int j)
1543 {
1544 SwapColumns(P->Constraint, P->NbConstraints, i, j);
1545 SwapColumns(P->Ray, P->NbRays, i, j);
1546 }
1547
1548 enum constraint {
1549 ALL_POS = 1 << 0,
1550 ONE_NEG = 1 << 1,
1551 INDEPENDENT = 1 << 2,
1552 };
1553
1554 evalue* barvinok_enumerate_e(Polyhedron *P,
1555 unsigned exist, unsigned nparam, unsigned MaxRays)
1556 {
1557 if (exist == 0) {
1558 Polyhedron *U = Universe_Polyhedron(nparam);
1559 evalue *EP = barvinok_enumerate_ev(P, U, MaxRays);
1560 //char *param_name[] = {"P", "Q", "R", "S", "T" };
1561 //print_evalue(stdout, EP, param_name);
1562 Polyhedron_Free(U);
1563 return EP;
1564 }
1565
1566 int nvar = P->Dimension - exist - nparam;
1567 int len = P->Dimension + 2;
1568
1569 //printf("%d %d %d\n", nvar, exist, nparam);
1570
1571 int r;
1572 for (r = 0; r < P->NbEq; ++r)
1573 if (First_Non_Zero(P->Constraint[r]+1+nvar, exist) != -1)
1574 break;
1575 if (r < P->NbEq) {
1576 Vector *row = Vector_Alloc(exist);
1577 Vector_Copy(P->Constraint[r]+1+nvar, row->p, exist);
1578 Matrix *M = unimodular_complete(row);
1579 Matrix *M2 = Matrix_Alloc(P->Dimension+1, P->Dimension+1);
1580 for (r = 0; r < nvar; ++r)
1581 value_set_si(M2->p[r][r], 1);
1582 for ( ; r < nvar+exist; ++r)
1583 Vector_Copy(M->p[r-nvar], M2->p[r]+nvar, exist);
1584 for ( ; r < P->Dimension+1; ++r)
1585 value_set_si(M2->p[r][r], 1);
1586 Polyhedron *T = Polyhedron_Image(P, M2, MaxRays);
1587 evalue *EP = barvinok_enumerate_e(T, exist-1, nparam, MaxRays);
1588 Polyhedron_Free(T);
1589 Matrix_Free(M2);
1590 Matrix_Free(M);
1591 Vector_Free(row);
1592 return EP;
1593 }
1594
1595 Vector *row = Vector_Alloc(len);
1596 value_set_si(row->p[0], 1);
1597
1598 Value f;
1599 value_init(f);
1600
1601 enum constraint info[exist];
1602 for (int i = 0; i < exist; ++i) {
1603 info[i] = ALL_POS;
1604 for (int l = P->NbEq; l < P->NbConstraints; ++l) {
1605 if (value_negz_p(P->Constraint[l][nvar+i+1]))
1606 continue;
1607 for (int u = P->NbEq; u < P->NbConstraints; ++u) {
1608 if (value_posz_p(P->Constraint[u][nvar+i+1]))
1609 continue;
1610 value_oppose(f, P->Constraint[u][nvar+i+1]);
1611 Vector_Combine(P->Constraint[l]+1, P->Constraint[u]+1, row->p+1,
1612 f, P->Constraint[l][nvar+i+1], len-1);
1613 if (!(info[i] & INDEPENDENT)) {
1614 int j;
1615 for (j = 0; j < exist; ++j)
1616 if (j != i && value_notzero_p(row->p[nvar+j+1]))
1617 break;
1618 if (j == exist) {
1619 //printf("independent: i: %d, l: %d, u: %d\n", i, l, u);
1620 info[i] = (constraint)(info[i] | INDEPENDENT);
1621 }
1622 }
1623 if (info[i] & ALL_POS) {
1624 value_addto(row->p[len-1], row->p[len-1],
1625 P->Constraint[l][nvar+i+1]);
1626 value_addto(row->p[len-1], row->p[len-1], f);
1627 value_multiply(f, f, P->Constraint[l][nvar+i+1]);
1628 value_substract(row->p[len-1], row->p[len-1], f);
1629 value_decrement(row->p[len-1], row->p[len-1]);
1630 Vector_Gcd(row->p+1, len - 2, &f);
1631 if (value_notone_p(f)) {
1632 Vector_AntiScale(row->p+1, row->p+1, f, len-2);
1633 mpz_fdiv_q(row->p[len-1], row->p[len-1], f);
1634 }
1635 value_set_si(f, -1);
1636 Vector_Scale(row->p+1, row->p+1, f, len-1);
1637 value_decrement(row->p[len-1], row->p[len-1]);
1638 Polyhedron *T = AddConstraints(row->p, 1, P, MaxRays);
1639 if (!emptyQ(T)) {
1640 //printf("not all_pos: i: %d, l: %d, u: %d\n", i, l, u);
1641 info[i] = (constraint)(info[i] ^ ALL_POS);
1642 }
1643 //puts("pos remainder");
1644 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
1645 Polyhedron_Free(T);
1646 }
1647 if (!(info[i] & ONE_NEG)) {
1648 int j;
1649 for (j = 0; j < exist; ++j)
1650 if (j != i &&
1651 (value_notzero_p(P->Constraint[l][nvar+j+1]) ||
1652 value_notzero_p(P->Constraint[u][nvar+j+1])))
1653 break;
1654 if (j == exist) {
1655 /* recalculate constant */
1656 value_oppose(f, P->Constraint[u][nvar+i+1]);
1657 Vector_Combine(P->Constraint[l]+len-1,
1658 P->Constraint[u]+len-1, row->p+len-1,
1659 f, P->Constraint[l][nvar+i+1], 1);
1660 value_multiply(f, f, P->Constraint[l][nvar+i+1]);
1661 value_substract(row->p[len-1], row->p[len-1], f);
1662 value_set_si(f, -1);
1663 Vector_Scale(row->p+1, row->p+1, f, len-1);
1664 value_decrement(row->p[len-1], row->p[len-1]);
1665 Vector_Gcd(row->p+1, len - 2, &f);
1666 if (value_notone_p(f)) {
1667 Vector_AntiScale(row->p+1, row->p+1, f, len-2);
1668 mpz_fdiv_q(row->p[len-1], row->p[len-1], f);
1669 }
1670 value_set_si(f, -1);
1671 Vector_Scale(row->p+1, row->p+1, f, len-1);
1672 value_decrement(row->p[len-1], row->p[len-1]);
1673 //puts("row");
1674 //Vector_Print(stdout, P_VALUE_FMT, row);
1675 Polyhedron *T = AddConstraints(row->p, 1, P, MaxRays);
1676 if (emptyQ(T)) {
1677 //printf("one_neg i: %d, l: %d, u: %d\n", i, l, u);
1678 info[i] = (constraint)(info[i] | ONE_NEG);
1679 }
1680 //puts("neg remainder");
1681 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
1682 Polyhedron_Free(T);
1683 }
1684 }
1685 if (!(info[i] & ALL_POS) && (info[i] & ONE_NEG))
1686 goto next;
1687 }
1688 }
1689 if (info[i] & ALL_POS)
1690 break;
1691 next:
1692 ;
1693 }
1694
1695 /*
1696 for (int i = 0; i < exist; ++i)
1697 printf("%i: %i\n", i, info[i]);
1698 */
1699 for (int i = 0; i < exist; ++i)
1700 if (info[i] & ALL_POS) {
1701 // Eliminate
1702 // Maybe we should chew off some of the fat here
1703 Matrix *M = Matrix_Alloc(P->Dimension, P->Dimension+1);
1704 for (int j = 0; j < P->Dimension; ++j)
1705 value_set_si(M->p[j][j + (j >= i+nvar)], 1);
1706 Polyhedron *T = Polyhedron_Image(P, M, MaxRays);
1707 Matrix_Free(M);
1708 evalue *EP = barvinok_enumerate_e(T, exist-1, nparam, MaxRays);
1709 Polyhedron_Free(T);
1710 return EP;
1711 }
1712 for (int i = 0; i < exist; ++i)
1713 if (info[i] & ONE_NEG) {
1714 value_clear(f);
1715 if (i == 0)
1716 return barvinok_enumerate_e(P, exist-1, nparam, MaxRays);
1717 else {
1718 Polyhedron *T = Polyhedron_Copy(P);
1719 SwapColumns(T, nvar+1, nvar+1+i);
1720 evalue *EP = barvinok_enumerate_e(T, exist-1, nparam, MaxRays);
1721 Polyhedron_Free(T);
1722 return EP;
1723 }
1724 }
1725 for (int i = 0; i < exist; ++i)
1726 if (info[i] & INDEPENDENT) {
1727 /* Find constraint again and split off negative part */
1728
1729 for (int l = P->NbEq; l < P->NbConstraints; ++l) {
1730 if (value_negz_p(P->Constraint[l][nvar+i+1]))
1731 continue;
1732 for (int u = P->NbEq; u < P->NbConstraints; ++u) {
1733 if (value_posz_p(P->Constraint[u][nvar+i+1]))
1734 continue;
1735 value_oppose(f, P->Constraint[u][nvar+i+1]);
1736 Vector_Combine(P->Constraint[l]+1, P->Constraint[u]+1,
1737 row->p+1,
1738 f, P->Constraint[l][nvar+i+1], len-1);
1739
1740 int j;
1741 for (j = 0; j < exist; ++j)
1742 if (j != i && value_notzero_p(row->p[nvar+j+1]))
1743 break;
1744 if (j != exist)
1745 continue;
1746
1747 //printf("l: %d, u: %d\n", l, u);
1748 value_multiply(f, f, P->Constraint[l][nvar+i+1]);
1749 value_substract(row->p[len-1], row->p[len-1], f);
1750 value_set_si(f, -1);
1751 Vector_Scale(row->p+1, row->p+1, f, len-1);
1752 value_decrement(row->p[len-1], row->p[len-1]);
1753 Vector_Gcd(row->p+1, len - 2, &f);
1754 if (value_notone_p(f)) {
1755 Vector_AntiScale(row->p+1, row->p+1, f, len-2);
1756 mpz_fdiv_q(row->p[len-1], row->p[len-1], f);
1757 }
1758 Polyhedron *neg = AddConstraints(row->p, 1, P, MaxRays);
1759 value_set_si(f, -1);
1760 Vector_Scale(row->p+1, row->p+1, f, len-1);
1761 value_decrement(row->p[len-1], row->p[len-1]);
1762 Polyhedron *pos = AddConstraints(row->p, 1, P, MaxRays);
1763
1764 assert(i == 0); // for now
1765 evalue *EP =
1766 barvinok_enumerate_e(neg, exist-1, nparam, MaxRays);
1767 evalue *E =
1768 barvinok_enumerate_e(pos, exist, nparam, MaxRays);
1769 eadd(E, EP);
1770 free_evalue_refs(E);
1771 free(E);
1772 Polyhedron_Free(neg);
1773 Polyhedron_Free(pos);
1774 value_clear(f);
1775 return EP;
1776 }
1777 }
1778 assert(0); // can't happen
1779 }
1780
1781 assert(0);
1782 }
lattice point enumeration library