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 #include "config.h"
  13 #include <barvinok/barvinok.h>
  14 #include <barvinok/genfun.h>
  15 #include <barvinok/options.h>
  16 #include <barvinok/sample.h>
  17 #include "bfcounter.h"
  18 #include "conversion.h"
  19 #include "counter.h"
  20 #include "tcounter.h"
  21 #include "decomposer.h"
  22 #include "euler.h"
  23 #include "lattice_point.h"
  24 #include "reduce_domain.h"
  25 #include "remove_equalities.h"
  26 #include "scale.h"
  27 #include "volume.h"
  28 #include "bernoulli.h"
  29 #include "param_util.h"
  30
  31 #ifdef NTL_STD_CXX
  32 using namespace NTL;
  33 #endif
  34 using std::cerr;
  35 using std::cout;
  36 using std::endl;
  37 using std::vector;
  38 using std::deque;
  39 using std::string;
  40 using std::ostringstream;
  41
  42 #define ALLOC(t,p) p = (t*)malloc(sizeof(*p))
  43
  44 class dpoly_n {
  45 public:
  46     Matrix *coeff;
  47     ~dpoly_n() {
  48         Matrix_Free(coeff);
  49     }
  50     dpoly_n(int d) {
  51         Value d0, one;
  52         value_init(d0);
  53         value_init(one);
  54         value_set_si(one, 1);
  55         coeff = Matrix_Alloc(d+1, d+1+1);
  56         value_set_si(coeff->p[0][0], 1);
  57         value_set_si(coeff->p[0][d+1], 1);
  58         for (int i = 1; i <= d; ++i) {
  59             value_multiply(coeff->p[i][0], coeff->p[i-1][0], d0);
  60             Vector_Combine(coeff->p[i-1], coeff->p[i-1]+1, coeff->p[i]+1,
  61                            one, d0, i);
  62             value_set_si(coeff->p[i][d+1], i);
  63             value_multiply(coeff->p[i][d+1], coeff->p[i][d+1], coeff->p[i-1][d+1]);
  64             value_decrement(d0, d0);
  65         }
  66         value_clear(d0);
  67         value_clear(one);
  68     }
  69     void div(dpoly& d, Vector *count, int sign) {
  70         int len = coeff->NbRows;
  71         Matrix * c = Matrix_Alloc(coeff->NbRows, coeff->NbColumns);
  72         Value tmp;
  73         value_init(tmp);
  74         for (int i = 0; i < len; ++i) {
  75             Vector_Copy(coeff->p[i], c->p[i], len+1);
  76             for (int j = 1; j <= i; ++j) {
  77                 value_multiply(tmp, d.coeff->p[j], c->p[i][len]);
  78                 value_oppose(tmp, tmp);
  79                 Vector_Combine(c->p[i], c->p[i-j], c->p[i],
  80                                c->p[i-j][len], tmp, len);
  81                 value_multiply(c->p[i][len], c->p[i][len], c->p[i-j][len]);
  82             }
  83             value_multiply(c->p[i][len], c->p[i][len], d.coeff->p[0]);
  84         }
  85         if (sign == -1) {
  86             value_set_si(tmp, -1);
  87             Vector_Scale(c->p[len-1], count->p, tmp, len);
  88             value_assign(count->p[len], c->p[len-1][len]);
  89         } else
  90             Vector_Copy(c->p[len-1], count->p, len+1);
  91         Vector_Normalize(count->p, len+1);
  92         value_clear(tmp);
  93         Matrix_Free(c);
  94     }
  95 };
  96
  97 const int MAX_TRY=10;
  98 /*
  99  * Searches for a vector that is not orthogonal to any
 100  * of the rays in rays.
 101  */
 102 static void nonorthog(mat_ZZ& rays, vec_ZZ& lambda)
 103 {
 104     int dim = rays.NumCols();
 105     bool found = false;
 106     lambda.SetLength(dim);
 107     if (dim == 0)
 108         return;
 109
 110     for (int i = 2; !found && i <= 50*dim; i+=4) {
 111         for (int j = 0; j < MAX_TRY; ++j) {
 112             for (int k = 0; k < dim; ++k) {
 113                 int r = random_int(i)+2;
 114                 int v = (2*(r%2)-1) * (r >> 1);
 115                 lambda[k] = v;
 116             }
 117             int k = 0;
 118             for (; k < rays.NumRows(); ++k)
 119                 if (lambda * rays[k] == 0)
 120                     break;
 121             if (k == rays.NumRows()) {
 122                 found = true;
 123                 break;
 124             }
 125         }
 126     }
 127     assert(found);
 128 }
 129
 130 static void add_rays(mat_ZZ& rays, Polyhedron *i, int *r, int nvar = -1,
 131                      bool all = false)
 132 {
 133     unsigned dim = i->Dimension;
 134     if (nvar == -1)
 135         nvar = dim;
 136     for (int k = 0; k < i->NbRays; ++k) {
 137         if (!value_zero_p(i->Ray[k][dim+1]))
 138             continue;
 139         if (!all && nvar != dim && First_Non_Zero(i->Ray[k]+1, nvar) == -1)
 140             continue;
 141         values2zz(i->Ray[k]+1, rays[(*r)++], nvar);
 142     }
 143 }
 144
 145 static void mask_r(Matrix *f, int nr, Vector *lcm, int p, Vector *val, evalue *ev)
 146 {
 147     unsigned nparam = lcm->Size;
 148
 149     if (p == nparam) {
 150         Vector * prod = Vector_Alloc(f->NbRows);
 151         Matrix_Vector_Product(f, val->p, prod->p);
 152         int isint = 1;
 153         for (int i = 0; i < nr; ++i) {
 154             value_modulus(prod->p[i], prod->p[i], f->p[i][nparam+1]);
 155             isint &= value_zero_p(prod->p[i]);
 156         }
 157         value_set_si(ev->d, 1);
 158         value_init(ev->x.n);
 159         value_set_si(ev->x.n, isint);
 160         Vector_Free(prod);
 161         return;
 162     }
 163
 164     Value tmp;
 165     value_init(tmp);
 166     if (value_one_p(lcm->p[p]))
 167         mask_r(f, nr, lcm, p+1, val, ev);
 168     else {
 169         value_assign(tmp, lcm->p[p]);
 170         value_set_si(ev->d, 0);
 171         ev->x.p = new_enode(periodic, VALUE_TO_INT(tmp), p+1);
 172         do {
 173             value_decrement(tmp, tmp);
 174             value_assign(val->p[p], tmp);
 175             mask_r(f, nr, lcm, p+1, val, &ev->x.p->arr[VALUE_TO_INT(tmp)]);
 176         } while (value_pos_p(tmp));
 177     }
 178     value_clear(tmp);
 179 }
 180
 181 static void mask_fractional(Matrix *f, evalue *factor)
 182 {
 183     int nr = f->NbRows, nc = f->NbColumns;
 184     int n;
 185     bool found = false;
 186     for (n = 0; n < nr && value_notzero_p(f->p[n][nc-1]); ++n)
 187         if (value_notone_p(f->p[n][nc-1]) &&
 188             value_notmone_p(f->p[n][nc-1]))
 189                 found = true;
 190     if (!found)
 191         return;
 192
 193     evalue EP;
 194     nr = n;
 195
 196     Value m;
 197     value_init(m);
 198
 199     evalue EV;
 200     value_init(EV.d);
 201     value_init(EV.x.n);
 202     value_set_si(EV.x.n, 1);
 203
 204     for (n = 0; n < nr; ++n) {
 205         value_assign(m, f->p[n][nc-1]);
 206         if (value_one_p(m) || value_mone_p(m))
 207             continue;
 208
 209         int j = normal_mod(f->p[n], nc-1, &m);
 210         if (j == nc-1) {
 211             free_evalue_refs(factor);
 212             value_init(factor->d);
 213             evalue_set_si(factor, 0, 1);
 214             break;
 215         }
 216         vec_ZZ row;
 217         values2zz(f->p[n], row, nc-1);
 218         ZZ g;
 219         value2zz(m, g);
 220         if (j < (nc-1)-1 && row[j] > g/2) {
 221             for (int k = j; k < (nc-1); ++k)
 222                 if (row[k] != 0)
 223                     row[k] = g - row[k];
 224         }
 225
 226         value_init(EP.d);
 227         value_set_si(EP.d, 0);
 228         EP.x.p = new_enode(relation, 2, 0);
 229         value_clear(EP.x.p->arr[1].d);
 230         EP.x.p->arr[1] = *factor;
 231         evalue *ev = &EP.x.p->arr[0];
 232         value_set_si(ev->d, 0);
 233         ev->x.p = new_enode(fractional, 3, -1);
 234         evalue_set_si(&ev->x.p->arr[1], 0, 1);
 235         evalue_set_si(&ev->x.p->arr[2], 1, 1);
 236         evalue *E = multi_monom(row);
 237         value_assign(EV.d, m);
 238         emul(&EV, E);
 239         value_clear(ev->x.p->arr[0].d);
 240         ev->x.p->arr[0] = *E;
 241         delete E;
 242         *factor = EP;
 243     }
 244
 245     value_clear(m);
 246     free_evalue_refs(&EV);
 247 }
 248
 249 /*
 250  *
 251  */
 252 static void mask_table(Matrix *f, evalue *factor)
 253 {
 254     int nr = f->NbRows, nc = f->NbColumns;
 255     int n;
 256     bool found = false;
 257     for (n = 0; n < nr && value_notzero_p(f->p[n][nc-1]); ++n)
 258         if (value_notone_p(f->p[n][nc-1]) &&
 259             value_notmone_p(f->p[n][nc-1]))
 260                 found = true;
 261     if (!found)
 262         return;
 263
 264     Value tmp;
 265     value_init(tmp);
 266     nr = n;
 267     unsigned np = nc - 2;
 268     Vector *lcm = Vector_Alloc(np);
 269     Vector *val = Vector_Alloc(nc);
 270     Vector_Set(val->p, 0, nc);
 271     value_set_si(val->p[np], 1);
 272     Vector_Set(lcm->p, 1, np);
 273     for (n = 0; n < nr; ++n) {
 274         if (value_one_p(f->p[n][nc-1]) ||
 275             value_mone_p(f->p[n][nc-1]))
 276             continue;
 277         for (int j = 0; j < np; ++j)
 278             if (value_notzero_p(f->p[n][j])) {
 279                 value_gcd(tmp, f->p[n][j], f->p[n][nc-1]);
 280                 value_division(tmp, f->p[n][nc-1], tmp);
 281                 value_lcm(lcm->p[j], tmp, lcm->p[j]);
 282             }
 283     }
 284     evalue EP;
 285     value_init(EP.d);
 286     mask_r(f, nr, lcm, 0, val, &EP);
 287     value_clear(tmp);
 288     Vector_Free(val);
 289     Vector_Free(lcm);
 290     emul(&EP,factor);
 291     free_evalue_refs(&EP);
 292 }
 293
 294 static void mask(Matrix *f, evalue *factor, barvinok_options *options)
 295 {
 296     if (options->lookup_table)
 297         mask_table(f, factor);
 298     else
 299         mask_fractional(f, factor);
 300 }
 301
 302 struct bfe_term : public bfc_term_base {
 303     vector<evalue *> factors;
 304
 305     bfe_term(int len) : bfc_term_base(len) {
 306     }
 307
 308     ~bfe_term() {
 309         for (int i = 0; i < factors.size(); ++i) {
 310             if (!factors[i])
 311                 continue;
 312             free_evalue_refs(factors[i]);
 313             delete factors[i];
 314         }
 315     }
 316 };
 317
 318 static void print_int_vector(int *v, int len, const char *name)
 319 {
 320     cerr << name << endl;
 321     for (int j = 0; j < len; ++j) {
 322         cerr << v[j] << " ";
 323     }
 324     cerr << endl;
 325 }
 326
 327 static void print_bfc_terms(mat_ZZ& factors, bfc_vec& v)
 328 {
 329     cerr << endl;
 330     cerr << "factors" << endl;
 331     cerr << factors << endl;
 332     for (int i = 0; i < v.size(); ++i) {
 333         cerr << "term: " << i << endl;
 334         print_int_vector(v[i]->powers, factors.NumRows(), "powers");
 335         cerr << "terms" << endl;
 336         cerr << v[i]->terms << endl;
 337         bfc_term* bfct = static_cast<bfc_term *>(v[i]);
 338         cerr << bfct->c << endl;
 339     }
 340 }
 341
 342 static void print_bfe_terms(mat_ZZ& factors, bfc_vec& v)
 343 {
 344     cerr << endl;
 345     cerr << "factors" << endl;
 346     cerr << factors << endl;
 347     for (int i = 0; i < v.size(); ++i) {
 348         cerr << "term: " << i << endl;
 349         print_int_vector(v[i]->powers, factors.NumRows(), "powers");
 350         cerr << "terms" << endl;
 351         cerr << v[i]->terms << endl;
 352         bfe_term* bfet = static_cast<bfe_term *>(v[i]);
 353         for (int j = 0; j < v[i]->terms.NumRows(); ++j) {
 354            const char * test[] = {"a", "b"};
 355            print_evalue(stderr, bfet->factors[j], test);
 356            fprintf(stderr, "\n");
 357         }
 358     }
 359 }
 360
 361 struct bfcounter : public bfcounter_base {
 362     mpq_t count;
 363     Value tz;
 364
 365     bfcounter(unsigned dim) : bfcounter_base(dim) {
 366         mpq_init(count);
 367         lower = 1;
 368         value_init(tz);
 369     }
 370     ~bfcounter() {
 371         mpq_clear(count);
 372         value_clear(tz);
 373     }
 374     virtual void base(mat_ZZ& factors, bfc_vec& v);
 375     virtual void get_count(Value *result) {
 376         assert(value_one_p(&count[0]._mp_den));
 377         value_assign(*result, &count[0]._mp_num);
 378     }
 379 };
 380
 381 void bfcounter::base(mat_ZZ& factors, bfc_vec& v)
 382 {
 383     unsigned nf = factors.NumRows();
 384
 385     for (int i = 0; i < v.size(); ++i) {
 386         bfc_term* bfct = static_cast<bfc_term *>(v[i]);
 387         int total_power = 0;
 388         // factor is always positive, so we always
 389         // change signs
 390         for (int k = 0; k < nf; ++k)
 391             total_power += v[i]->powers[k];
 392
 393         int j;
 394         for (j = 0; j < nf; ++j)
 395             if (v[i]->powers[j] > 0)
 396                 break;
 397
 398         zz2value(factors[j][0], tz);
 399         dpoly D(total_power, tz, 1);
 400         for (int k = 1; k < v[i]->powers[j]; ++k) {
 401             zz2value(factors[j][0], tz);
 402             dpoly fact(total_power, tz, 1);
 403             D *= fact;
 404         }
 405         for ( ; ++j < nf; )
 406             for (int k = 0; k < v[i]->powers[j]; ++k) {
 407                 zz2value(factors[j][0], tz);
 408                 dpoly fact(total_power, tz, 1);
 409                 D *= fact;
 410             }
 411
 412         for (int k = 0; k < v[i]->terms.NumRows(); ++k) {
 413             zz2value(v[i]->terms[k][0], tz);
 414             dpoly n(total_power, tz);
 415             mpq_set_si(tcount, 0, 1);
 416             n.div(D, tcount, 1);
 417             if (total_power % 2)
 418                 bfct->c[k].n = -bfct->c[k].n;
 419             zz2value(bfct->c[k].n, tn);
 420             zz2value(bfct->c[k].d, td);
 421
 422             mpz_mul(mpq_numref(tcount), mpq_numref(tcount), tn);
 423             mpz_mul(mpq_denref(tcount), mpq_denref(tcount), td);
 424             mpq_canonicalize(tcount);
 425             mpq_add(count, count, tcount);
 426         }
 427         delete v[i];
 428     }
 429 }
 430
 431
 432 /* Check whether the polyhedron is unbounded and if so,
 433  * check whether it has any (and therefore an infinite number of)
 434  * integer points.
 435  * If one of the vertices is integer, then we are done.
 436  * Otherwise, transform the polyhedron such that one of the rays
 437  * is the first unit vector and cut it off at a height that ensures
 438  * that if the whole polyhedron has any points, then the remaining part
 439  * has integer points.  In particular we add the largest coefficient
 440  * of a ray to the highest vertex (rounded up).
 441  */
 442 static bool Polyhedron_is_infinite(Polyhedron *P, Value* result,
 443                                    barvinok_options *options)
 444 {
 445     int r = 0;
 446     Matrix *M, *M2;
 447     Value c, tmp;
 448     Value g;
 449     bool first;
 450     Vector *v;
 451     Value offset, size;
 452     Polyhedron *R;
 453
 454     if (P->NbBid == 0)
 455         for (; r < P->NbRays; ++r)
 456             if (value_zero_p(P->Ray[r][P->Dimension+1]))
 457                 break;
 458     if (P->NbBid == 0 && r == P->NbRays)
 459         return false;
 460
 461     if (options->count_sample_infinite) {
 462         Vector *sample;
 463
 464         sample = Polyhedron_Sample(P, options);
 465         if (!sample)
 466             value_set_si(*result, 0);
 467         else {
 468             value_set_si(*result, -1);
 469             Vector_Free(sample);
 470         }
 471         return true;
 472     }
 473
 474     for (int i = 0; i < P->NbRays; ++i)
 475         if (value_one_p(P->Ray[i][1+P->Dimension])) {
 476             value_set_si(*result, -1);
 477             return true;
 478         }
 479
 480     value_init(g);
 481     M = Matrix_Alloc(P->Dimension+1, P->Dimension+1);
 482     Vector_Gcd(P->Ray[r]+1, P->Dimension, &g);
 483     Vector_AntiScale(P->Ray[r]+1, M->p[0], g, P->Dimension+1);
 484     int ok = unimodular_complete(M, 1);
 485     assert(ok);
 486     value_set_si(M->p[P->Dimension][P->Dimension], 1);
 487     M2 = Transpose(M);
 488     Matrix_Free(M);
 489     P = Polyhedron_Preimage(P, M2, 0);
 490     Matrix_Free(M2);
 491     value_clear(g);
 492
 493     first = true;
 494     value_init(offset);
 495     value_init(size);
 496     value_init(tmp);
 497     value_set_si(size, 0);
 498
 499     for (int i = 0; i < P->NbBid; ++i) {
 500         value_absolute(tmp, P->Ray[i][1]);
 501         if (value_gt(tmp, size))
 502             value_assign(size, tmp);
 503     }
 504     for (int i = P->NbBid; i < P->NbRays; ++i) {
 505         if (value_zero_p(P->Ray[i][P->Dimension+1])) {
 506             if (value_gt(P->Ray[i][1], size))
 507                 value_assign(size, P->Ray[i][1]);
 508             continue;
 509         }
 510         mpz_cdiv_q(tmp, P->Ray[i][1], P->Ray[i][P->Dimension+1]);
 511         if (first || value_gt(tmp, offset)) {
 512             value_assign(offset, tmp);
 513             first = false;
 514         }
 515     }
 516     value_addto(offset, offset, size);
 517     value_clear(size);
 518     value_clear(tmp);
 519
 520     v = Vector_Alloc(P->Dimension+2);
 521     value_set_si(v->p[0], 1);
 522     value_set_si(v->p[1], -1);
 523     value_assign(v->p[1+P->Dimension], offset);
 524     R = AddConstraints(v->p, 1, P, options->MaxRays);
 525     Polyhedron_Free(P);
 526     P = R;
 527
 528     value_clear(offset);
 529     Vector_Free(v);
 530
 531     value_init(c);
 532     barvinok_count_with_options(P, &c, options);
 533     Polyhedron_Free(P);
 534     if (value_zero_p(c))
 535         value_set_si(*result, 0);
 536     else
 537         value_set_si(*result, -1);
 538     value_clear(c);
 539
 540     return true;
 541 }
 542
 543 static void barvinok_count_f(Polyhedron *P, Value* result,
 544                              barvinok_options *options);
 545
 546 void barvinok_count_with_options(Polyhedron *P, Value* result,
 547                                  struct barvinok_options *options)
 548 {
 549     unsigned dim;
 550     int allocated = 0;
 551     Polyhedron *Q;
 552     bool infinite = false;
 553
 554     if (P->next)
 555         fprintf(stderr,
 556             "barvinok_count: input is a union; only first polyhedron is counted\n");
 557
 558     if (emptyQ2(P)) {
 559         value_set_si(*result, 0);
 560         return;
 561     }
 562     if (P->NbEq != 0) {
 563         Q = NULL;
 564         do {
 565             P = remove_equalities(P, options->MaxRays);
 566             P = DomainConstraintSimplify(P, options->MaxRays);
 567             if (Q)
 568                 Polyhedron_Free(Q);
 569             Q = P;
 570         } while (!emptyQ(P) && P->NbEq != 0);
 571         if (emptyQ(P)) {
 572             Polyhedron_Free(P);
 573             value_set_si(*result, 0);
 574             return;
 575         }
 576         allocated = 1;
 577     }
 578     if (Polyhedron_is_infinite(P, result, options)) {
 579         if (allocated)
 580             Polyhedron_Free(P);
 581         return;
 582     }
 583     if (P->Dimension == 0) {
 584         /* Test whether the constraints are satisfied */
 585         POL_ENSURE_VERTICES(P);
 586         value_set_si(*result, !emptyQ(P));
 587         if (allocated)
 588             Polyhedron_Free(P);
 589         return;
 590     }
 591     Q = Polyhedron_Factor(P, 0, NULL, options->MaxRays);
 592     if (Q) {
 593         if (allocated)
 594             Polyhedron_Free(P);
 595         P = Q;
 596         allocated = 1;
 597     }
 598
 599     barvinok_count_f(P, result, options);
 600     if (value_neg_p(*result))
 601         infinite = true;
 602     if (Q && P->next && value_notzero_p(*result)) {
 603         Value factor;
 604         value_init(factor);
 605
 606         for (Q = P->next; Q; Q = Q->next) {
 607             barvinok_count_f(Q, &factor, options);
 608             if (value_neg_p(factor)) {
 609                 infinite = true;
 610                 continue;
 611             } else if (Q->next && value_zero_p(factor)) {
 612                 value_set_si(*result, 0);
 613                 break;
 614             }
 615             value_multiply(*result, *result, factor);
 616         }
 617
 618         value_clear(factor);
 619     }
 620
 621     if (allocated)
 622         Domain_Free(P);
 623     if (infinite)
 624         value_set_si(*result, -1);
 625 }
 626
 627 void barvinok_count(Polyhedron *P, Value* result, unsigned NbMaxCons)
 628 {
 629     barvinok_options *options = barvinok_options_new_with_defaults();
 630     options->MaxRays = NbMaxCons;
 631     barvinok_count_with_options(P, result, options);
 632     barvinok_options_free(options);
 633 }
 634
 635 static void barvinok_count_f(Polyhedron *P, Value* result,
 636                              barvinok_options *options)
 637 {
 638     if (emptyQ2(P)) {
 639         value_set_si(*result, 0);
 640         return;
 641     }
 642
 643     if (P->Dimension == 1)
 644         return Line_Length(P, result);
 645
 646     int c = P->NbConstraints;
 647     POL_ENSURE_FACETS(P);
 648     if (c != P->NbConstraints || P->NbEq != 0) {
 649         Polyhedron *next = P->next;
 650         P->next = NULL;
 651         barvinok_count_with_options(P, result, options);
 652         P->next = next;
 653         return;
 654     }
 655
 656     POL_ENSURE_VERTICES(P);
 657
 658     if (Polyhedron_is_infinite(P, result, options))
 659         return;
 660
 661     np_base *cnt;
 662     if (options->incremental_specialization == BV_SPECIALIZATION_BF)
 663         cnt = new bfcounter(P->Dimension);
 664     else if (options->incremental_specialization == BV_SPECIALIZATION_DF)
 665         cnt = new icounter(P->Dimension);
 666     else if (options->incremental_specialization == BV_SPECIALIZATION_TODD)
 667         cnt = new tcounter(P->Dimension, options->max_index);
 668     else
 669         cnt = new counter(P->Dimension, options->max_index);
 670     cnt->start(P, options);
 671
 672     cnt->get_count(result);
 673     delete cnt;
 674 }
 675
 676 static void uni_polynom(int param, Vector *c, evalue *EP)
 677 {
 678     unsigned dim = c->Size-2;
 679     value_init(EP->d);
 680     value_set_si(EP->d,0);
 681     EP->x.p = new_enode(polynomial, dim+1, param+1);
 682     for (int j = 0; j <= dim; ++j)
 683         evalue_set(&EP->x.p->arr[j], c->p[j], c->p[dim+1]);
 684 }
 685
 686 typedef evalue * evalue_p;
 687
 688 struct enumerator_base {
 689     unsigned dim;
 690     evalue ** vE;
 691     evalue mone;
 692     vertex_decomposer *vpd;
 693
 694     enumerator_base(unsigned dim, vertex_decomposer *vpd)
 695     {
 696         this->dim = dim;
 697         this->vpd = vpd;
 698
 699         vE = new evalue_p[vpd->PP->nbV];
 700         for (int j = 0; j < vpd->PP->nbV; ++j)
 701             vE[j] = 0;
 702
 703         value_init(mone.d);
 704         evalue_set_si(&mone, -1, 1);
 705     }
 706
 707     void decompose_at(Param_Vertices *V, int _i, barvinok_options *options) {
 708         //this->pVD = pVD;
 709
 710         vE[_i] = new evalue;
 711         value_init(vE[_i]->d);
 712         evalue_set_si(vE[_i], 0, 1);
 713
 714         vpd->decompose_at_vertex(V, _i, options);
 715     }
 716
 717     virtual ~enumerator_base() {
 718         for (int j = 0; j < vpd->PP->nbV; ++j)
 719             if (vE[j]) {
 720                 free_evalue_refs(vE[j]);
 721                 delete vE[j];
 722             }
 723         delete [] vE;
 724
 725         free_evalue_refs(&mone);
 726     }
 727
 728     static enumerator_base *create(Polyhedron *P, unsigned dim,
 729                                      Param_Polyhedron *PP,
 730                                      barvinok_options *options);
 731 };
 732
 733 struct enumerator : public signed_cone_consumer, public vertex_decomposer,
 734                     public enumerator_base {
 735     vec_ZZ lambda;
 736     vec_ZZ den;
 737     term_info num;
 738     Vector *c;
 739     mpq_t count;
 740     Value tz;
 741
 742     enumerator(Polyhedron *P, unsigned dim, Param_Polyhedron *PP) :
 743                 vertex_decomposer(PP, *this), enumerator_base(dim, this) {
 744         randomvector(P, lambda, dim);
 745         den.SetLength(dim);
 746         c = Vector_Alloc(dim+2);
 747
 748         mpq_init(count);
 749         value_init(tz);
 750     }
 751
 752     ~enumerator() {
 753         mpq_clear(count);
 754         Vector_Free(c);
 755         value_clear(tz);
 756     }
 757
 758     virtual void handle(const signed_cone& sc, barvinok_options *options);
 759 };
 760
 761 void enumerator::handle(const signed_cone& sc, barvinok_options *options)
 762 {
 763     int sign = sc.sign;
 764     int r = 0;
 765     assert(sc.rays.NumRows() == dim);
 766     for (int k = 0; k < dim; ++k) {
 767         if (lambda * sc.rays[k] == 0)
 768             throw Orthogonal;
 769     }
 770
 771     lattice_point(V, sc.rays, lambda, &num, sc.det, options);
 772     den = sc.rays * lambda;
 773
 774     if (dim % 2)
 775         sign = -sign;
 776
 777     zz2value(den[0], tz);
 778     dpoly n(dim, tz, 1);
 779     for (int k = 1; k < dim; ++k) {
 780         zz2value(den[k], tz);
 781         dpoly fact(dim, tz, 1);
 782         n *= fact;
 783     }
 784     if (num.E != NULL) {
 785         dpoly_n d(dim);
 786         d.div(n, c, sign);
 787         for (unsigned long i = 0; i < sc.det; ++i) {
 788             evalue *EV = evalue_polynomial(c, num.E[i]);
 789             eadd(EV, vE[vert]);
 790             evalue_free(EV);
 791             free_evalue_refs(num.E[i]);
 792             delete num.E[i];
 793         }
 794         delete [] num.E;
 795     } else {
 796         mpq_set_si(count, 0, 1);
 797         if (num.constant.length() == 1) {
 798             zz2value(num.constant[0], tz);
 799             dpoly d(dim, tz);
 800             d.div(n, count, sign);
 801         } else {
 802             dpoly_n d(dim);
 803             d.div(n, c, sign);
 804             Value x, sum, acc;
 805             value_init(x);
 806             value_init(acc);
 807             for (unsigned long i = 0; i < sc.det; ++i) {
 808                 value_assign(acc, c->p[dim]);
 809                 zz2value(num.constant[i], x);
 810                 for (int j = dim-1; j >= 0; --j) {
 811                     value_multiply(acc, acc, x);
 812                     value_addto(acc, acc, c->p[j]);
 813                 }
 814                 value_addto(mpq_numref(count), mpq_numref(count), acc);
 815             }
 816             mpz_set(mpq_denref(count), c->p[dim+1]);
 817             value_clear(acc);
 818             value_clear(x);
 819         }
 820         evalue EV;
 821         value_init(EV.d);
 822         evalue_set(&EV, &count[0]._mp_num, &count[0]._mp_den);
 823         eadd(&EV, vE[vert]);
 824         free_evalue_refs(&EV);
 825     }
 826 }
 827
 828 struct ienumerator_base : enumerator_base {
 829     evalue ** E_vertex;
 830
 831     ienumerator_base(unsigned dim, vertex_decomposer *vpd) :
 832                         enumerator_base(dim,vpd) {
 833         E_vertex = new evalue_p[dim];
 834     }
 835
 836     virtual ~ienumerator_base() {
 837         delete [] E_vertex;
 838     }
 839
 840     evalue *E_num(int i, int d) {
 841         return E_vertex[i + (dim-d)];
 842     }
 843 };
 844
 845 struct cumulator {
 846     evalue *factor;
 847     evalue *v;
 848     dpoly_r *r;
 849
 850     cumulator(evalue *factor, evalue *v, dpoly_r *r) :
 851         factor(factor), v(v), r(r) {}
 852
 853     void cumulate(barvinok_options *options);
 854
 855     virtual void add_term(const vector<int>& powers, evalue *f2) = 0;
 856     virtual ~cumulator() {}
 857 };
 858
 859 void cumulator::cumulate(barvinok_options *options)
 860 {
 861     evalue cum;  // factor * 1 * E_num[0]/1 * (E_num[0]-1)/2 *...
 862     evalue f;
 863     evalue t;   // E_num[0] - (m-1)
 864     evalue *cst;
 865     evalue mone;
 866
 867     if (options->lookup_table) {
 868         value_init(mone.d);
 869         evalue_set_si(&mone, -1, 1);
 870     }
 871
 872     value_init(cum.d);
 873     evalue_copy(&cum, factor);
 874     value_init(f.d);
 875     value_init(f.x.n);
 876     value_set_si(f.d, 1);
 877     value_set_si(f.x.n, 1);
 878     value_init(t.d);
 879     evalue_copy(&t, v);
 880
 881     if (!options->lookup_table) {
 882         for (cst = &t; value_zero_p(cst->d); ) {
 883             if (cst->x.p->type == fractional)
 884                 cst = &cst->x.p->arr[1];
 885             else
 886                 cst = &cst->x.p->arr[0];
 887         }
 888     }
 889
 890     for (int m = 0; m < r->len; ++m) {
 891         if (m > 0) {
 892             if (m > 1) {
 893                 value_set_si(f.d, m);
 894                 emul(&f, &cum);
 895                 if (!options->lookup_table)
 896                     value_subtract(cst->x.n, cst->x.n, cst->d);
 897                 else
 898                     eadd(&mone, &t);
 899             }
 900             emul(&t, &cum);
 901         }
 902         dpoly_r_term_list& current = r->c[r->len-1-m];
 903         dpoly_r_term_list::iterator j;
 904         for (j = current.begin(); j != current.end(); ++j) {
 905             if ((*j)->coeff == 0)
 906                 continue;
 907             evalue *f2 = new evalue;
 908             value_init(f2->d);
 909             value_init(f2->x.n);
 910             zz2value((*j)->coeff, f2->x.n);
 911             zz2value(r->denom, f2->d);
 912             emul(&cum, f2);
 913
 914             add_term((*j)->powers, f2);
 915         }
 916     }
 917     free_evalue_refs(&f);
 918     free_evalue_refs(&t);
 919     free_evalue_refs(&cum);
 920     if (options->lookup_table)
 921         free_evalue_refs(&mone);
 922 }
 923
 924 struct E_poly_term {
 925     vector<int> powers;
 926     evalue  *E;
 927 };
 928
 929 struct ie_cum : public cumulator {
 930     vector<E_poly_term *> terms;
 931
 932     ie_cum(evalue *factor, evalue *v, dpoly_r *r) : cumulator(factor, v, r) {}
 933
 934     virtual void add_term(const vector<int>& powers, evalue *f2);
 935 };
 936
 937 void ie_cum::add_term(const vector<int>& powers, evalue *f2)
 938 {
 939     int k;
 940     for (k = 0; k < terms.size(); ++k) {
 941         if (terms[k]->powers == powers) {
 942             eadd(f2, terms[k]->E);
 943             free_evalue_refs(f2);
 944             delete f2;
 945             break;
 946         }
 947     }
 948     if (k >= terms.size()) {
 949         E_poly_term *ET = new E_poly_term;
 950         ET->powers = powers;
 951         ET->E = f2;
 952         terms.push_back(ET);
 953     }
 954 }
 955
 956 struct ienumerator : public signed_cone_consumer, public vertex_decomposer,
 957                      public ienumerator_base {
 958     //Polyhedron *pVD;
 959     mat_ZZ den;
 960     mat_ZZ vertex;
 961     mpq_t tcount;
 962     Value tz;
 963
 964     ienumerator(Polyhedron *P, unsigned dim, Param_Polyhedron *PP) :
 965                 vertex_decomposer(PP, *this), ienumerator_base(dim, this) {
 966         vertex.SetDims(1, dim);
 967
 968         den.SetDims(dim, dim);
 969         mpq_init(tcount);
 970         value_init(tz);
 971     }
 972
 973     ~ienumerator() {
 974         mpq_clear(tcount);
 975         value_clear(tz);
 976     }
 977
 978     virtual void handle(const signed_cone& sc, barvinok_options *options);
 979     void reduce(evalue *factor, const mat_ZZ& num, const mat_ZZ& den_f,
 980                 barvinok_options *options);
 981 };
 982
 983 void ienumerator::reduce(evalue *factor, const mat_ZZ& num, const mat_ZZ& den_f,
 984                          barvinok_options *options)
 985 {
 986     unsigned len = den_f.NumRows();  // number of factors in den
 987     unsigned dim = num.NumCols();
 988     assert(num.NumRows() == 1);
 989
 990     if (dim == 0) {
 991         eadd(factor, vE[vert]);
 992         return;
 993     }
 994
 995     vec_ZZ den_s;
 996     mat_ZZ den_r;
 997     vec_ZZ num_s;
 998     mat_ZZ num_p;
 999
1000     split_one(num, num_s, num_p, den_f, den_s, den_r);
1001
1002     vec_ZZ den_p;
1003     den_p.SetLength(len);
1004
1005     ZZ one;
1006     one = 1;
1007     normalize(one, num_s, num_p, den_s, den_p, den_r);
1008     if (one != 1)
1009         emul(&mone, factor);
1010
1011     int only_param = 0;
1012     int no_param = 0;
1013     for (int k = 0; k < len; ++k) {
1014         if (den_p[k] == 0)
1015             ++no_param;
1016         else if (den_s[k] == 0)
1017             ++only_param;
1018     }
1019     if (no_param == 0) {
1020         reduce(factor, num_p, den_r, options);
1021     } else {
1022         int k, l;
1023         mat_ZZ pden;
1024         pden.SetDims(only_param, dim-1);
1025
1026         for (k = 0, l = 0; k < len; ++k)
1027             if (den_s[k] == 0)
1028                 pden[l++] = den_r[k];
1029
1030         for (k = 0; k < len; ++k)
1031             if (den_p[k] == 0)
1032                 break;
1033
1034         zz2value(num_s[0], tz);
1035         dpoly n(no_param, tz);
1036         zz2value(den_s[k], tz);
1037         dpoly D(no_param, tz, 1);
1038         for ( ; ++k < len; )
1039             if (den_p[k] == 0) {
1040                 zz2value(den_s[k], tz);
1041                 dpoly fact(no_param, tz, 1);
1042                 D *= fact;
1043             }
1044
1045         dpoly_r * r = 0;
1046         // if no_param + only_param == len then all powers
1047         // below will be all zero
1048         if (no_param + only_param == len) {
1049             if (E_num(0, dim) != 0)
1050                 r = new dpoly_r(n, len);
1051             else {
1052                 mpq_set_si(tcount, 0, 1);
1053                 one = 1;
1054                 n.div(D, tcount, 1);
1055
1056                 if (value_notzero_p(mpq_numref(tcount))) {
1057                     evalue f;
1058                     value_init(f.d);
1059                     value_init(f.x.n);
1060                     value_assign(f.x.n, mpq_numref(tcount));
1061                     value_assign(f.d, mpq_denref(tcount));
1062                     emul(&f, factor);
1063                     reduce(factor, num_p, pden, options);
1064                     free_evalue_refs(&f);
1065                 }
1066                 return;
1067             }
1068         } else {
1069             for (k = 0; k < len; ++k) {
1070                 if (den_s[k] == 0 || den_p[k] == 0)
1071                     continue;
1072
1073                 zz2value(den_s[k], tz);
1074                 dpoly pd(no_param-1, tz, 1);
1075
1076                 int l;
1077                 for (l = 0; l < k; ++l)
1078                     if (den_r[l] == den_r[k])
1079                         break;
1080
1081                 if (r == 0)
1082                     r = new dpoly_r(n, pd, l, len);
1083                 else {
1084                     dpoly_r *nr = new dpoly_r(r, pd, l, len);
1085                     delete r;
1086                     r = nr;
1087                 }
1088             }
1089         }
1090         dpoly_r *rc = r->div(D);
1091         delete r;
1092         r = rc;
1093         if (E_num(0, dim) == 0) {
1094             int common = pden.NumRows();
1095             dpoly_r_term_list& final = r->c[r->len-1];
1096             int rows;
1097             evalue t;
1098             evalue f;
1099             value_init(f.d);
1100             value_init(f.x.n);
1101             zz2value(r->denom, f.d);
1102             dpoly_r_term_list::iterator j;
1103             for (j = final.begin(); j != final.end(); ++j) {
1104                 if ((*j)->coeff == 0)
1105                     continue;
1106                 rows = common;
1107                 for (int k = 0; k < r->dim; ++k) {
1108                     int n = (*j)->powers[k];
1109                     if (n == 0)
1110                         continue;
1111                     pden.SetDims(rows+n, pden.NumCols());
1112                     for (int l = 0; l < n; ++l)
1113                         pden[rows+l] = den_r[k];
1114                     rows += n;
1115                 }
1116                 value_init(t.d);
1117                 evalue_copy(&t, factor);
1118                 zz2value((*j)->coeff, f.x.n);
1119                 emul(&f, &t);
1120                 reduce(&t, num_p, pden, options);
1121                 free_evalue_refs(&t);
1122             }
1123             free_evalue_refs(&f);
1124         } else {
1125             ie_cum cum(factor, E_num(0, dim), r);
1126             cum.cumulate(options);
1127
1128             int common = pden.NumRows();
1129             int rows;
1130             for (int j = 0; j < cum.terms.size(); ++j) {
1131                 rows = common;
1132                 pden.SetDims(rows, pden.NumCols());
1133                 for (int k = 0; k < r->dim; ++k) {
1134                     int n = cum.terms[j]->powers[k];
1135                     if (n == 0)
1136                         continue;
1137                     pden.SetDims(rows+n, pden.NumCols());
1138                     for (int l = 0; l < n; ++l)
1139                         pden[rows+l] = den_r[k];
1140                     rows += n;
1141                 }
1142                 reduce(cum.terms[j]->E, num_p, pden, options);
1143                 free_evalue_refs(cum.terms[j]->E);
1144                 delete cum.terms[j]->E;
1145                 delete cum.terms[j];
1146             }
1147         }
1148         delete r;
1149     }
1150 }
1151
1152 static int type_offset(enode *p)
1153 {
1154    return p->type == fractional ? 1 :
1155           p->type == flooring ? 1 : 0;
1156 }
1157
1158 static int edegree(evalue *e)
1159 {
1160     int d = 0;
1161     enode *p;
1162
1163     if (value_notzero_p(e->d))
1164         return 0;
1165
1166     p = e->x.p;
1167     int i = type_offset(p);
1168     if (p->size-i-1 > d)
1169         d = p->size - i - 1;
1170     for (; i < p->size; i++) {
1171         int d2 = edegree(&p->arr[i]);
1172         if (d2 > d)
1173             d = d2;
1174     }
1175     return d;
1176 }
1177
1178 void ienumerator::handle(const signed_cone& sc, barvinok_options *options)
1179 {
1180     assert(sc.det == 1);
1181     assert(sc.rays.NumRows() == dim);
1182
1183     lattice_point(V, sc.rays, vertex[0], E_vertex, options);
1184
1185     den = sc.rays;
1186
1187     evalue one;
1188     value_init(one.d);
1189     evalue_set_si(&one, sc.sign, 1);
1190     reduce(&one, vertex, den, options);
1191     free_evalue_refs(&one);
1192
1193     for (int i = 0; i < dim; ++i)
1194         if (E_vertex[i]) {
1195             free_evalue_refs(E_vertex[i]);
1196             delete E_vertex[i];
1197         }
1198 }
1199
1200 struct bfenumerator : public vertex_decomposer, public bf_base,
1201                       public ienumerator_base {
1202     evalue *factor;
1203
1204     bfenumerator(Polyhedron *P, unsigned dim, Param_Polyhedron *PP) :
1205                     vertex_decomposer(PP, *this),
1206                     bf_base(dim), ienumerator_base(dim, this) {
1207         lower = 0;
1208         factor = NULL;
1209     }
1210
1211     ~bfenumerator() {
1212     }
1213
1214     virtual void handle(const signed_cone& sc, barvinok_options *options);
1215     virtual void base(mat_ZZ& factors, bfc_vec& v);
1216
1217     bfc_term_base* new_bf_term(int len) {
1218         bfe_term* t = new bfe_term(len);
1219         return t;
1220     }
1221
1222     virtual void set_factor(bfc_term_base *t, int k, int change) {
1223         bfe_term* bfet = static_cast<bfe_term *>(t);
1224         factor = bfet->factors[k];
1225         assert(factor != NULL);
1226         bfet->factors[k] = NULL;
1227         if (change)
1228             emul(&mone, factor);
1229     }
1230
1231     virtual void set_factor(bfc_term_base *t, int k, mpq_t &q, int change) {
1232         bfe_term* bfet = static_cast<bfe_term *>(t);
1233         factor = bfet->factors[k];
1234         assert(factor != NULL);
1235         bfet->factors[k] = NULL;
1236
1237         evalue f;
1238         value_init(f.d);
1239         value_init(f.x.n);
1240         if (change)
1241             value_oppose(f.x.n, mpq_numref(q));
1242         else
1243             value_assign(f.x.n, mpq_numref(q));
1244         value_assign(f.d, mpq_denref(q));
1245         emul(&f, factor);
1246         free_evalue_refs(&f);
1247     }
1248
1249     virtual void set_factor(bfc_term_base *t, int k, const QQ& c, int change) {
1250         bfe_term* bfet = static_cast<bfe_term *>(t);
1251
1252         factor = new evalue;
1253
1254         evalue f;
1255         value_init(f.d);
1256         value_init(f.x.n);
1257         zz2value(c.n, f.x.n);
1258         if (change)
1259             value_oppose(f.x.n, f.x.n);
1260         zz2value(c.d, f.d);
1261
1262         value_init(factor->d);
1263         evalue_copy(factor, bfet->factors[k]);
1264         emul(&f, factor);
1265         free_evalue_refs(&f);
1266     }
1267
1268     void set_factor(evalue *f, int change) {
1269         if (change)
1270             emul(&mone, f);
1271         factor = f;
1272     }
1273
1274     virtual void insert_term(bfc_term_base *t, int i) {
1275         bfe_term* bfet = static_cast<bfe_term *>(t);
1276         int len = t->terms.NumRows()-1; // already increased by one
1277
1278         bfet->factors.resize(len+1);
1279         for (int j = len; j > i; --j) {
1280             bfet->factors[j] = bfet->factors[j-1];
1281             t->terms[j] = t->terms[j-1];
1282         }
1283         bfet->factors[i] = factor;
1284         factor = NULL;
1285     }
1286
1287     virtual void update_term(bfc_term_base *t, int i) {
1288         bfe_term* bfet = static_cast<bfe_term *>(t);
1289
1290         eadd(factor, bfet->factors[i]);
1291         free_evalue_refs(factor);
1292         delete factor;
1293     }
1294
1295     virtual bool constant_vertex(int dim) { return E_num(0, dim) == 0; }
1296
1297     virtual void cum(bf_reducer *bfr, bfc_term_base *t, int k, dpoly_r *r,
1298                      barvinok_options *options);
1299 };
1300
1301 enumerator_base *enumerator_base::create(Polyhedron *P, unsigned dim,
1302                                          Param_Polyhedron *PP,
1303                                          barvinok_options *options)
1304 {
1305     enumerator_base *eb;
1306
1307     if (options->incremental_specialization == BV_SPECIALIZATION_BF)
1308         eb = new bfenumerator(P, dim, PP);
1309     else if (options->incremental_specialization == BV_SPECIALIZATION_DF)
1310         eb = new ienumerator(P, dim, PP);
1311     else
1312         eb = new enumerator(P, dim, PP);
1313
1314     return eb;
1315 }
1316
1317 struct bfe_cum : public cumulator {
1318     bfenumerator *bfe;
1319     bfc_term_base *told;
1320     int k;
1321     bf_reducer *bfr;
1322
1323     bfe_cum(evalue *factor, evalue *v, dpoly_r *r, bf_reducer *bfr,
1324             bfc_term_base *t, int k, bfenumerator *e) :
1325             cumulator(factor, v, r), told(t), k(k),
1326             bfr(bfr), bfe(e) {
1327     }
1328
1329     virtual void add_term(const vector<int>& powers, evalue *f2);
1330 };
1331
1332 void bfe_cum::add_term(const vector<int>& powers, evalue *f2)
1333 {
1334     bfr->update_powers(powers);
1335
1336     bfc_term_base * t = bfe->find_bfc_term(bfr->vn, bfr->npowers, bfr->nnf);
1337     bfe->set_factor(f2, bfr->l_changes % 2);
1338     bfe->add_term(t, told->terms[k], bfr->l_extra_num);
1339 }
1340
1341 void bfenumerator::cum(bf_reducer *bfr, bfc_term_base *t, int k,
1342                        dpoly_r *r, barvinok_options *options)
1343 {
1344     bfe_term* bfet = static_cast<bfe_term *>(t);
1345     bfe_cum cum(bfet->factors[k], E_num(0, bfr->d), r, bfr, t, k, this);
1346     cum.cumulate(options);
1347 }
1348
1349 void bfenumerator::base(mat_ZZ& factors, bfc_vec& v)
1350 {
1351     for (int i = 0; i < v.size(); ++i) {
1352         assert(v[i]->terms.NumRows() == 1);
1353         evalue *factor = static_cast<bfe_term *>(v[i])->factors[0];
1354         eadd(factor, vE[vert]);
1355         delete v[i];
1356     }
1357 }
1358
1359 void bfenumerator::handle(const signed_cone& sc, barvinok_options *options)
1360 {
1361     assert(sc.det == 1);
1362     assert(sc.rays.NumRows() == enumerator_base::dim);
1363
1364     bfe_term* t = new bfe_term(enumerator_base::dim);
1365     vector< bfc_term_base * > v;
1366     v.push_back(t);
1367
1368     t->factors.resize(1);
1369
1370     t->terms.SetDims(1, enumerator_base::dim);
1371     lattice_point(V, sc.rays, t->terms[0], E_vertex, options);
1372
1373     // the elements of factors are always lexpositive
1374     mat_ZZ   factors;
1375     int s = setup_factors(sc.rays, factors, t, sc.sign);
1376
1377     t->factors[0] = new evalue;
1378     value_init(t->factors[0]->d);
1379     evalue_set_si(t->factors[0], s, 1);
1380     reduce(factors, v, options);
1381
1382     for (int i = 0; i < enumerator_base::dim; ++i)
1383         if (E_vertex[i]) {
1384             free_evalue_refs(E_vertex[i]);
1385             delete E_vertex[i];
1386         }
1387 }
1388
1389 static evalue* barvinok_enumerate_ev_f(Polyhedron *P, Polyhedron* C,
1390                                        barvinok_options *options);
1391
1392 /* Destroys C */
1393 static evalue* barvinok_enumerate_cst(Polyhedron *P, Polyhedron* C,
1394                                       struct barvinok_options *options)
1395 {
1396     evalue *eres;
1397
1398     if (emptyQ2(C)) {
1399         Polyhedron_Free(C);
1400         return evalue_zero();
1401     }
1402
1403     ALLOC(evalue, eres);
1404     value_init(eres->d);
1405     value_set_si(eres->d, 0);
1406     eres->x.p = new_enode(partition, 2, C->Dimension);
1407     EVALUE_SET_DOMAIN(eres->x.p->arr[0],
1408                       DomainConstraintSimplify(C, options->MaxRays));
1409     value_set_si(eres->x.p->arr[1].d, 1);
1410     value_init(eres->x.p->arr[1].x.n);
1411     if (emptyQ2(P))
1412         value_set_si(eres->x.p->arr[1].x.n, 0);
1413     else
1414         barvinok_count_with_options(P, &eres->x.p->arr[1].x.n, options);
1415
1416     return eres;
1417 }
1418
1419 static evalue* enumerate(Polyhedron *P, Polyhedron* C,
1420                                         struct barvinok_options *options)
1421 {
1422     Polyhedron *next;
1423     Polyhedron *Porig = P;
1424     Polyhedron *Corig = C;
1425     Polyhedron *CEq = NULL, *rVD;
1426     int r = 0;
1427     unsigned nparam = C->Dimension;
1428     evalue *eres;
1429     Matrix *CP = NULL;
1430
1431     evalue factor;
1432     value_init(factor.d);
1433     evalue_set_si(&factor, 1, 1);
1434
1435     /* for now */
1436     POL_ENSURE_FACETS(P);
1437     POL_ENSURE_VERTICES(P);
1438     POL_ENSURE_FACETS(C);
1439     POL_ENSURE_VERTICES(C);
1440
1441     if (C->Dimension == 0 || emptyQ(P) || emptyQ(C)) {
1442 constant:
1443         if (CEq == Porig)
1444             CEq = Polyhedron_Copy(CEq);
1445         eres = barvinok_enumerate_cst(P, CEq ? CEq : Polyhedron_Copy(C), options);
1446 out:
1447         if (CP) {
1448             evalue_backsubstitute(eres, CP, options->MaxRays);
1449             Matrix_Free(CP);
1450         }
1451
1452         emul(&factor, eres);
1453         if (options->approximation_method == BV_APPROX_DROP) {
1454             if (options->polynomial_approximation == BV_APPROX_SIGN_UPPER)
1455                 evalue_frac2polynomial(eres, 1, options->MaxRays);
1456             if (options->polynomial_approximation == BV_APPROX_SIGN_LOWER)
1457                 evalue_frac2polynomial(eres, -1, options->MaxRays);
1458             if (options->polynomial_approximation == BV_APPROX_SIGN_APPROX)
1459                 evalue_frac2polynomial(eres, 0, options->MaxRays);
1460         }
1461         reduce_evalue(eres);
1462         free_evalue_refs(&factor);
1463         if (P != Porig)
1464             Domain_Free(P);
1465         if (C != Corig)
1466             Polyhedron_Free(C);
1467
1468         return eres;
1469     }
1470     if (Polyhedron_is_unbounded(P, nparam, options->MaxRays))
1471         goto constant;
1472
1473     if (P->NbEq != 0) {
1474         Matrix *f;
1475         P = remove_equalities_p(Polyhedron_Copy(P), P->Dimension-nparam, &f,
1476                                 options->MaxRays);
1477         mask(f, &factor, options);
1478         Matrix_Free(f);
1479     }
1480     if (P->Dimension == nparam) {
1481         CEq = P;
1482         P = Universe_Polyhedron(0);
1483         goto constant;
1484     }
1485     if (P->NbEq != 0 || C->NbEq != 0) {
1486         Polyhedron *Q = P;
1487         Polyhedron *D = C;
1488         remove_all_equalities(&P, &C, &CP, NULL, nparam, options->MaxRays);
1489         if (C != D && D != Corig)
1490             Polyhedron_Free(D);
1491         if (P != Q && Q != Porig)
1492             Domain_Free(Q);
1493         eres = enumerate(P, C, options);
1494         goto out;
1495     }
1496
1497     Polyhedron *T = Polyhedron_Factor(P, nparam, NULL, options->MaxRays);
1498     if (T || (P->Dimension == nparam+1)) {
1499         Polyhedron *Q;
1500         Polyhedron *C2;
1501         for (Q = T ? T : P; Q; Q = Q->next) {
1502             Polyhedron *next = Q->next;
1503             Q->next = NULL;
1504
1505             Polyhedron *QC = Q;
1506             if (Q->Dimension != C->Dimension)
1507                 QC = Polyhedron_Project(Q, nparam);
1508
1509             C2 = C;
1510             C = DomainIntersection(C, QC, options->MaxRays);
1511             if (C2 != Corig)
1512                 Polyhedron_Free(C2);
1513             if (QC != Q)
1514                 Polyhedron_Free(QC);
1515
1516             Q->next = next;
1517         }
1518     }
1519     if (T) {
1520         if (P != Porig)
1521             Polyhedron_Free(P);
1522         P = T;
1523         if (T->Dimension == C->Dimension) {
1524             P = T->next;
1525             T->next = NULL;
1526             Polyhedron_Free(T);
1527         }
1528     }
1529
1530     next = P->next;
1531     P->next = NULL;
1532     eres = barvinok_enumerate_ev_f(P, C, options);
1533     P->next = next;
1534
1535     if (P->next) {
1536         Polyhedron *Q;
1537         evalue *f;
1538
1539         for (Q = P->next; Q; Q = Q->next) {
1540             Polyhedron *next = Q->next;
1541             Q->next = NULL;
1542
1543             f = barvinok_enumerate_ev_f(Q, C, options);
1544             emul(f, eres);
1545             evalue_free(f);
1546
1547             Q->next = next;
1548         }
1549     }
1550
1551     goto out;
1552 }
1553
1554 evalue* barvinok_enumerate_with_options(Polyhedron *P, Polyhedron* C,
1555                                         struct barvinok_options *options)
1556 {
1557     Polyhedron *next, *Cnext, *C1;
1558     Polyhedron *Corig = C;
1559     evalue *eres;
1560
1561     if (P->next)
1562         fprintf(stderr,
1563     "barvinok_enumerate: input is a union; only first polyhedron is enumerated\n");
1564
1565     if (C->next)
1566         fprintf(stderr,
1567     "barvinok_enumerate: context is a union; only first polyhedron is considered\n");
1568
1569     Cnext = C->next;
1570     C->next = NULL;
1571     C1 = Polyhedron_Project(P, C->Dimension);
1572     C = DomainIntersection(C, C1, options->MaxRays);
1573     Polyhedron_Free(C1);
1574     next = P->next;
1575     P->next = NULL;
1576
1577     if (options->approximation_method == BV_APPROX_BERNOULLI)
1578         eres = Bernoulli_sum(P, C, options);
1579     else
1580         eres = enumerate(P, C, options);
1581     Domain_Free(C);
1582
1583     P->next= next;
1584     Corig->next = Cnext;
1585
1586     return eres;
1587 }
1588
1589 evalue* barvinok_enumerate_ev(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
1590 {
1591     evalue *E;
1592     barvinok_options *options = barvinok_options_new_with_defaults();
1593     options->MaxRays = MaxRays;
1594     E = barvinok_enumerate_with_options(P, C, options);
1595     barvinok_options_free(options);
1596     return E;
1597 }
1598
1599 evalue *Param_Polyhedron_Enumerate(Param_Polyhedron *PP, Polyhedron *P,
1600                                    Polyhedron *C,
1601                                    struct barvinok_options *options)
1602 {
1603     evalue *eres;
1604     Param_Domain *D;
1605     unsigned nparam = C->Dimension;
1606     unsigned dim = P->Dimension - nparam;
1607
1608     int nd;
1609     for (nd = 0, D=PP->D; D; ++nd, D=D->next);
1610     evalue_section *s = new evalue_section[nd];
1611
1612     enumerator_base *et = NULL;
1613 try_again:
1614     if (et)
1615         delete et;
1616
1617     et = enumerator_base::create(P, dim, PP, options);
1618
1619     Polyhedron *TC = true_context(P, C, options->MaxRays);
1620     FORALL_REDUCED_DOMAIN(PP, TC, nd, options, i, D, rVD)
1621         Param_Vertices *V;
1622
1623         s[i].E = evalue_zero();
1624         s[i].D = rVD;
1625
1626         FORALL_PVertex_in_ParamPolyhedron(V,D,PP) // _i is internal counter
1627             if (!et->vE[_i])
1628                 try {
1629                     et->decompose_at(V, _i, options);
1630                 } catch (OrthogonalException &e) {
1631                     FORALL_REDUCED_DOMAIN_RESET;
1632                     for (; i >= 0; --i) {
1633                         evalue_free(s[i].E);
1634                         Domain_Free(s[i].D);
1635                     }
1636                     goto try_again;
1637                 }
1638             eadd(et->vE[_i] , s[i].E);
1639         END_FORALL_PVertex_in_ParamPolyhedron;
1640         evalue_range_reduction_in_domain(s[i].E, rVD);
1641     END_FORALL_REDUCED_DOMAIN
1642     Polyhedron_Free(TC);
1643
1644     delete et;
1645     eres = evalue_from_section_array(s, nd);
1646     delete [] s;
1647
1648     return eres;
1649 }
1650
1651 static evalue* barvinok_enumerate_ev_f(Polyhedron *P, Polyhedron* C,
1652                                        barvinok_options *options)
1653 {
1654     unsigned nparam = C->Dimension;
1655     bool do_scale = options->approximation_method == BV_APPROX_SCALE;
1656
1657     if (options->approximation_method == BV_APPROX_VOLUME)
1658         return Param_Polyhedron_Volume(P, C, options);
1659
1660     if (P->Dimension - nparam == 1 && !do_scale)
1661         return ParamLine_Length(P, C, options);
1662
1663     Param_Polyhedron *PP = NULL;
1664     evalue *eres;
1665
1666     if (do_scale) {
1667         eres = scale_bound(P, C, options);
1668         if (eres)
1669             return eres;
1670     }
1671
1672     PP = Polyhedron2Param_Polyhedron(P, C, options);
1673
1674     if (do_scale)
1675         eres = scale(PP, P, C, options);
1676     else
1677         eres = Param_Polyhedron_Enumerate(PP, P, C, options);
1678
1679     if (PP)
1680         Param_Polyhedron_Free(PP);
1681
1682     return eres;
1683 }
1684
1685 Enumeration* barvinok_enumerate(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
1686 {
1687     evalue *EP = barvinok_enumerate_ev(P, C, MaxRays);
1688
1689     return partition2enumeration(EP);
1690 }
1691
1692 evalue* barvinok_enumerate_union(Polyhedron *D, Polyhedron* C, unsigned MaxRays)
1693 {
1694     evalue *EP;
1695     gen_fun *gf = barvinok_enumerate_union_series(D, C, MaxRays);
1696     EP = *gf;
1697     delete gf;
1698     return EP;
1699 }
1700
1701 evalue *barvinok_summate(evalue *e, int nvar, struct barvinok_options *options)
1702 {
1703     if (options->summation == BV_SUM_EULER)
1704         return euler_summate(e, nvar, options);
1705     else if (options->summation == BV_SUM_BERNOULLI)
1706         return Bernoulli_sum_evalue(e, nvar, options);
1707     else
1708         return evalue_sum(e, nvar, options->MaxRays);
1709 }