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