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