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