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