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