barvinok.cc

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