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