lexmin.cc

   1 #include <iostream>
   2 #include <vector>
   3 #include <map>
   4 #include <set>
   5 #include <gmp.h>
   6 #include <NTL/vec_ZZ.h>
   7 #include <NTL/mat_ZZ.h>
   8 extern "C" {
   9 #include <polylib/polylibgmp.h>
  10 }
  11 #include <barvinok/barvinok.h>
  12 #include <barvinok/evalue.h>
  13 #include <barvinok/options.h>
  14 #include <barvinok/util.h>
  15 #include "conversion.h"
  16 #include "decomposer.h"
  17 #include "lattice_point.h"
  18 #include "reduce_domain.h"
  19 #include "mat_util.h"
  20 #include "combine.h"
  21 #include "edomain.h"
  22 #include "evalue_util.h"
  23 #include "sample.h"
  24 #include "config.h"
  25
  26 #ifdef NTL_STD_CXX
  27 using namespace NTL;
  28 #endif
  29
  30 using std::vector;
  31 using std::map;
  32 using std::cerr;
  33 using std::cout;
  34 using std::endl;
  35 using std::ostream;
  36
  37 /* RANGE : normal range for evalutations (-RANGE -> RANGE) */
  38 #define RANGE 50
  39
  40 /* SRANGE : small range for evalutations */
  41 #define SRANGE 15
  42
  43 /* if dimension >= BIDDIM, use SRANGE */
  44 #define BIGDIM 5
  45
  46 /* VSRANGE : very small range for evalutations */
  47 #define VSRANGE 5
  48
  49 /* if dimension >= VBIDDIM, use VSRANGE */
  50 #define VBIGDIM 8
  51
  52 #ifndef HAVE_GETOPT_H
  53 #define getopt_long(a,b,c,d,e) getopt(a,b,c)
  54 #else
  55 #include <getopt.h>
  56 #define NO_EMPTINESS_CHECK  256
  57 struct option lexmin_options[] = {
  58     { "verify",     no_argument,  0,  'T' },
  59     { "print-all",  no_argument,  0,  'A' },
  60     { "no-emptiness-check", no_argument, 0, NO_EMPTINESS_CHECK },
  61     { "min",        required_argument,  0,  'm' },
  62     { "max",        required_argument,  0,  'M' },
  63     { "range",      required_argument,  0,  'r' },
  64     { "version",    no_argument,  0,  'V' },
  65     { 0, 0, 0, 0 }
  66 };
  67 #endif
  68
  69 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
  70
  71 static int type_offset(enode *p)
  72 {
  73    return p->type == fractional ? 1 :
  74           p->type == flooring ? 1 : 0;
  75 }
  76
  77 struct indicator_term {
  78     int sign;
  79     int pos;
  80     mat_ZZ den;
  81     evalue **vertex;
  82
  83     indicator_term(unsigned dim, int pos) {
  84         den.SetDims(0, dim);
  85         vertex = new evalue* [dim];
  86         this->pos = pos;
  87         sign = 0;
  88     }
  89     indicator_term(unsigned dim) {
  90         den.SetDims(dim, dim);
  91         vertex = new evalue* [dim];
  92         pos = -1;
  93     }
  94     indicator_term(const indicator_term& src) {
  95         sign = src.sign;
  96         pos = src.pos;
  97         den = src.den;
  98         unsigned dim = den.NumCols();
  99         vertex = new evalue* [dim];
 100         for (int i = 0; i < dim; ++i) {
 101             vertex[i] = new evalue();
 102             value_init(vertex[i]->d);
 103             evalue_copy(vertex[i], src.vertex[i]);
 104         }
 105     }
 106     ~indicator_term() {
 107         unsigned dim = den.NumCols();
 108         for (int i = 0; i < dim; ++i) {
 109             free_evalue_refs(vertex[i]);
 110             delete vertex[i];
 111         }
 112         delete [] vertex;
 113     }
 114     void print(ostream& os, char **p) const;
 115     void substitute(Matrix *T);
 116     void normalize();
 117     void substitute(evalue *fract, evalue *val);
 118     void substitute(int pos, evalue *val);
 119     void reduce_in_domain(Polyhedron *D);
 120     bool is_opposite(indicator_term *neg);
 121 };
 122
 123 bool indicator_term::is_opposite(indicator_term *neg)
 124 {
 125     if (sign + neg->sign != 0)
 126         return false;
 127     if (den != neg->den)
 128         return false;
 129     for (int k = 0; k < den.NumCols(); ++k)
 130         if (!eequal(vertex[k], neg->vertex[k]))
 131             return false;
 132     return true;
 133 }
 134
 135 void indicator_term::reduce_in_domain(Polyhedron *D)
 136 {
 137     for (int k = 0; k < den.NumCols(); ++k) {
 138         reduce_evalue_in_domain(vertex[k], D);
 139         if (evalue_range_reduction_in_domain(vertex[k], D))
 140             reduce_evalue(vertex[k]);
 141     }
 142 }
 143
 144 void indicator_term::print(ostream& os, char **p) const
 145 {
 146     unsigned dim = den.NumCols();
 147     unsigned factors = den.NumRows();
 148     if (sign == 0)
 149         os << " s ";
 150     else if (sign > 0)
 151         os << " + ";
 152     else
 153         os << " - ";
 154     os << '[';
 155     for (int i = 0; i < dim; ++i) {
 156         if (i)
 157             os << ',';
 158         evalue_print(os, vertex[i], p);
 159     }
 160     os << ']';
 161     for (int i = 0; i < factors; ++i) {
 162         os << " + t" << i << "*[";
 163         for (int j = 0; j < dim; ++j) {
 164             if (j)
 165                 os << ',';
 166             os << den[i][j];
 167         }
 168         os << "]";
 169     }
 170     if (sign == 0)
 171         os << " (" << pos << ")";
 172 }
 173
 174 /* Perform the substitution specified by T on the variables.
 175  * T has dimension (newdim+nparam+1) x (olddim + nparam + 1).
 176  * The substitution is performed as in gen_fun::substitute
 177  */
 178 void indicator_term::substitute(Matrix *T)
 179 {
 180     unsigned dim = den.NumCols();
 181     unsigned nparam = T->NbColumns - dim - 1;
 182     unsigned newdim = T->NbRows - nparam - 1;
 183     evalue **newvertex;
 184     mat_ZZ trans;
 185     matrix2zz(T, trans, newdim, dim);
 186     trans = transpose(trans);
 187     den *= trans;
 188     newvertex = new evalue* [newdim];
 189
 190     vec_ZZ v;
 191     v.SetLength(nparam+1);
 192
 193     evalue factor;
 194     value_init(factor.d);
 195     value_set_si(factor.d, 1);
 196     value_init(factor.x.n);
 197     for (int i = 0; i < newdim; ++i) {
 198         values2zz(T->p[i]+dim, v, nparam+1);
 199         newvertex[i] = multi_monom(v);
 200
 201         for (int j = 0; j < dim; ++j) {
 202             if (value_zero_p(T->p[i][j]))
 203                 continue;
 204             evalue term;
 205             value_init(term.d);
 206             evalue_copy(&term, vertex[j]);
 207             value_assign(factor.x.n, T->p[i][j]);
 208             emul(&factor, &term);
 209             eadd(&term, newvertex[i]);
 210             free_evalue_refs(&term);
 211         }
 212     }
 213     free_evalue_refs(&factor);
 214     for (int i = 0; i < dim; ++i) {
 215         free_evalue_refs(vertex[i]);
 216         delete vertex[i];
 217     }
 218     delete [] vertex;
 219     vertex = newvertex;
 220 }
 221
 222 static void evalue_add_constant(evalue *e, ZZ v)
 223 {
 224     Value tmp;
 225     value_init(tmp);
 226
 227     /* go down to constant term */
 228     while (value_zero_p(e->d))
 229         e = &e->x.p->arr[type_offset(e->x.p)];
 230     /* and add v */
 231     zz2value(v, tmp);
 232     value_multiply(tmp, tmp, e->d);
 233     value_addto(e->x.n, e->x.n, tmp);
 234
 235     value_clear(tmp);
 236 }
 237
 238 /* Make all powers in denominator lexico-positive */
 239 void indicator_term::normalize()
 240 {
 241     vec_ZZ extra_vertex;
 242     extra_vertex.SetLength(den.NumCols());
 243     for (int r = 0; r < den.NumRows(); ++r) {
 244         for (int k = 0; k < den.NumCols(); ++k) {
 245             if (den[r][k] == 0)
 246                 continue;
 247             if (den[r][k] > 0)
 248                 break;
 249             sign = -sign;
 250             den[r] = -den[r];
 251             extra_vertex += den[r];
 252             break;
 253         }
 254     }
 255     for (int k = 0; k < extra_vertex.length(); ++k)
 256         if (extra_vertex[k] != 0)
 257             evalue_add_constant(vertex[k], extra_vertex[k]);
 258 }
 259
 260 static void substitute(evalue *e, evalue *fract, evalue *val)
 261 {
 262     evalue *t;
 263
 264     for (t = e; value_zero_p(t->d); t = &t->x.p->arr[type_offset(t->x.p)]) {
 265         if (t->x.p->type == fractional && eequal(&t->x.p->arr[0], fract))
 266             break;
 267     }
 268     if (value_notzero_p(t->d))
 269         return;
 270
 271     free_evalue_refs(&t->x.p->arr[0]);
 272     evalue *term = &t->x.p->arr[2];
 273     enode *p = t->x.p;
 274     value_clear(t->d);
 275     *t = t->x.p->arr[1];
 276
 277     emul(val, term);
 278     eadd(term, e);
 279     free_evalue_refs(term);
 280     free(p);
 281
 282     reduce_evalue(e);
 283 }
 284
 285 void indicator_term::substitute(evalue *fract, evalue *val)
 286 {
 287     unsigned dim = den.NumCols();
 288     for (int i = 0; i < dim; ++i) {
 289         ::substitute(vertex[i], fract, val);
 290     }
 291 }
 292
 293 static void substitute(evalue *e, int pos, evalue *val)
 294 {
 295     evalue *t;
 296
 297     for (t = e; value_zero_p(t->d); t = &t->x.p->arr[type_offset(t->x.p)]) {
 298         if (t->x.p->type == polynomial && t->x.p->pos == pos)
 299             break;
 300     }
 301     if (value_notzero_p(t->d))
 302         return;
 303
 304     evalue *term = &t->x.p->arr[1];
 305     enode *p = t->x.p;
 306     value_clear(t->d);
 307     *t = t->x.p->arr[0];
 308
 309     emul(val, term);
 310     eadd(term, e);
 311     free_evalue_refs(term);
 312     free(p);
 313
 314     reduce_evalue(e);
 315 }
 316
 317 void indicator_term::substitute(int pos, evalue *val)
 318 {
 319     unsigned dim = den.NumCols();
 320     for (int i = 0; i < dim; ++i) {
 321         ::substitute(vertex[i], pos, val);
 322     }
 323 }
 324
 325 struct indicator_constructor : public polar_decomposer, public vertex_decomposer {
 326     vec_ZZ vertex;
 327     vector<indicator_term*> *terms;
 328     Matrix *T;  /* Transformation to original space */
 329     Param_Polyhedron *PP;
 330
 331     indicator_constructor(Polyhedron *P, unsigned dim, Param_Polyhedron *PP,
 332                           Matrix *T) :
 333                 vertex_decomposer(P, PP->nbV, this), T(T), PP(PP) {
 334         vertex.SetLength(dim);
 335         terms = new vector<indicator_term*>[nbV];
 336     }
 337     ~indicator_constructor() {
 338         for (int i = 0; i < nbV; ++i)
 339             for (int j = 0; j < terms[i].size(); ++j)
 340                 delete terms[i][j];
 341         delete [] terms;
 342     }
 343     void substitute(Matrix *T);
 344     void normalize();
 345     void print(ostream& os, char **p);
 346
 347     virtual void handle_polar(Polyhedron *P, int sign);
 348 };
 349
 350 void indicator_constructor::handle_polar(Polyhedron *C, int s)
 351 {
 352     unsigned dim = vertex.length();
 353
 354     assert(C->NbRays-1 == dim);
 355
 356     indicator_term *term = new indicator_term(dim);
 357     term->sign = s;
 358     terms[vert].push_back(term);
 359
 360     lattice_point(V, C, vertex, term->vertex);
 361
 362     for (int r = 0; r < dim; ++r) {
 363         values2zz(C->Ray[r]+1, term->den[r], dim);
 364         for (int j = 0; j < dim; ++j) {
 365             if (term->den[r][j] == 0)
 366                 continue;
 367             if (term->den[r][j] > 0)
 368                 break;
 369             term->sign = -term->sign;
 370             term->den[r] = -term->den[r];
 371             vertex += term->den[r];
 372             break;
 373         }
 374     }
 375
 376     for (int i = 0; i < dim; ++i) {
 377         if (!term->vertex[i]) {
 378             term->vertex[i] = new evalue();
 379             value_init(term->vertex[i]->d);
 380             value_init(term->vertex[i]->x.n);
 381             zz2value(vertex[i], term->vertex[i]->x.n);
 382             value_set_si(term->vertex[i]->d, 1);
 383             continue;
 384         }
 385         if (vertex[i] == 0)
 386             continue;
 387         evalue_add_constant(term->vertex[i], vertex[i]);
 388     }
 389
 390     if (T) {
 391         term->substitute(T);
 392         term->normalize();
 393     }
 394
 395     lex_order_rows(term->den);
 396 }
 397
 398 void indicator_constructor::print(ostream& os, char **p)
 399 {
 400     for (int i = 0; i < nbV; ++i)
 401         for (int j = 0; j < terms[i].size(); ++j) {
 402             os << "i: " << i << ", j: " << j << endl;
 403             terms[i][j]->print(os, p);
 404             os << endl;
 405         }
 406 }
 407
 408 void indicator_constructor::normalize()
 409 {
 410     for (int i = 0; i < nbV; ++i)
 411         for (int j = 0; j < terms[i].size(); ++j) {
 412             vec_ZZ vertex;
 413             vertex.SetLength(terms[i][j]->den.NumCols());
 414             for (int r = 0; r < terms[i][j]->den.NumRows(); ++r) {
 415                 for (int k = 0; k < terms[i][j]->den.NumCols(); ++k) {
 416                     if (terms[i][j]->den[r][k] == 0)
 417                         continue;
 418                     if (terms[i][j]->den[r][k] > 0)
 419                         break;
 420                     terms[i][j]->sign = -terms[i][j]->sign;
 421                     terms[i][j]->den[r] = -terms[i][j]->den[r];
 422                     vertex += terms[i][j]->den[r];
 423                     break;
 424                 }
 425             }
 426             lex_order_rows(terms[i][j]->den);
 427             for (int k = 0; k < vertex.length(); ++k)
 428                 if (vertex[k] != 0)
 429                     evalue_add_constant(terms[i][j]->vertex[k], vertex[k]);
 430         }
 431 }
 432
 433 struct indicator;
 434
 435 enum order_sign { order_lt, order_le, order_eq, order_ge, order_gt, order_unknown };
 436
 437 struct partial_order {
 438     indicator *ind;
 439
 440     map<indicator_term *, int > pred;
 441     map<indicator_term *, vector<indicator_term * > > lt;
 442     map<indicator_term *, vector<indicator_term * > > le;
 443     map<indicator_term *, vector<indicator_term * > > eq;
 444
 445     map<indicator_term *, vector<indicator_term * > > pending;
 446
 447     partial_order(indicator *ind) : ind(ind) {}
 448     void copy(const partial_order& order,
 449               map< indicator_term *, indicator_term * > old2new);
 450
 451     order_sign compare(indicator_term *a, indicator_term *b);
 452     void set_equal(indicator_term *a, indicator_term *b);
 453     void unset_le(indicator_term *a, indicator_term *b);
 454
 455     bool compared(indicator_term* a, indicator_term* b);
 456     void add(indicator_term* it, std::set<indicator_term *> *filter);
 457     void remove(indicator_term* it);
 458
 459     void print(ostream& os, char **p);
 460 };
 461
 462 void partial_order::unset_le(indicator_term *a, indicator_term *b)
 463 {
 464     vector<indicator_term *>::iterator i;
 465     i = find(le[a].begin(), le[a].end(), b);
 466     le[a].erase(i);
 467     pred[b]--;
 468     i = find(pending[a].begin(), pending[a].end(), b);
 469     if (i != pending[a].end())
 470         pending[a].erase(i);
 471 }
 472
 473 void partial_order::set_equal(indicator_term *a, indicator_term *b)
 474 {
 475     if (eq[a].size() == 0)
 476         eq[a].push_back(a);
 477     if (eq[b].size() == 0)
 478         eq[b].push_back(b);
 479     a = eq[a][0];
 480     b = eq[b][0];
 481     assert(a != b);
 482     if (b < a) {
 483         indicator_term *c = a;
 484         a = b;
 485         b = c;
 486     }
 487
 488     indicator_term *base = a;
 489
 490     map<indicator_term *, vector<indicator_term * > >::iterator i;
 491
 492     for (int j = 0; j < eq[b].size(); ++j) {
 493         eq[base].push_back(eq[b][j]);
 494         eq[eq[b][j]][0] = base;
 495     }
 496     eq[b].resize(1);
 497
 498     i = lt.find(b);
 499     if (i != lt.end()) {
 500         for (int j = 0; j < lt[b].size(); ++j) {
 501             if (find(eq[base].begin(), eq[base].end(), lt[b][j]) != eq[base].end())
 502                 pred[lt[b][j]]--;
 503             else if (find(lt[base].begin(), lt[base].end(), lt[b][j])
 504                         != lt[base].end())
 505                 pred[lt[b][j]]--;
 506             else
 507                 lt[base].push_back(lt[b][j]);
 508         }
 509         lt.erase(i);
 510     }
 511
 512     i = le.find(b);
 513     if (i != le.end()) {
 514         for (int j = 0; j < le[b].size(); ++j) {
 515             if (find(eq[base].begin(), eq[base].end(), le[b][j]) != eq[base].end())
 516                 pred[le[b][j]]--;
 517             else if (find(le[base].begin(), le[base].end(), le[b][j])
 518                         != le[base].end())
 519                 pred[le[b][j]]--;
 520             else
 521                 le[base].push_back(le[b][j]);
 522         }
 523         le.erase(i);
 524     }
 525
 526     i = pending.find(base);
 527     if (i != pending.end()) {
 528         vector<indicator_term * > old = pending[base];
 529         pending[base].clear();
 530         for (int j = 0; j < old.size(); ++j) {
 531             if (find(eq[base].begin(), eq[base].end(), old[j]) == eq[base].end())
 532                 pending[base].push_back(old[j]);
 533         }
 534     }
 535
 536     i = pending.find(b);
 537     if (i != pending.end()) {
 538         for (int j = 0; j < pending[b].size(); ++j) {
 539             if (find(eq[base].begin(), eq[base].end(), pending[b][j]) == eq[base].end())
 540                 pending[base].push_back(pending[b][j]);
 541         }
 542         pending.erase(i);
 543     }
 544 }
 545
 546 void partial_order::copy(const partial_order& order,
 547                          map< indicator_term *, indicator_term * > old2new)
 548 {
 549     map<indicator_term *, vector<indicator_term * > >::const_iterator i;
 550     map<indicator_term *, int >::const_iterator j;
 551
 552     for (j = order.pred.begin(); j != order.pred.end(); ++j)
 553         pred[old2new[(*j).first]] = (*j).second;
 554
 555     for (i = order.lt.begin(); i != order.lt.end(); ++i) {
 556         for (int j = 0; j < (*i).second.size(); ++j)
 557             lt[old2new[(*i).first]].push_back(old2new[(*i).second[j]]);
 558     }
 559     for (i = order.le.begin(); i != order.le.end(); ++i) {
 560         for (int j = 0; j < (*i).second.size(); ++j)
 561             le[old2new[(*i).first]].push_back(old2new[(*i).second[j]]);
 562     }
 563     for (i = order.eq.begin(); i != order.eq.end(); ++i) {
 564         for (int j = 0; j < (*i).second.size(); ++j)
 565             eq[old2new[(*i).first]].push_back(old2new[(*i).second[j]]);
 566     }
 567     for (i = order.pending.begin(); i != order.pending.end(); ++i) {
 568         for (int j = 0; j < (*i).second.size(); ++j)
 569             pending[old2new[(*i).first]].push_back(old2new[(*i).second[j]]);
 570     }
 571 }
 572
 573 static void interval_minmax(Polyhedron *I, int *min, int *max)
 574 {
 575     assert(I->Dimension == 1);
 576     *min = 1;
 577     *max = -1;
 578     POL_ENSURE_VERTICES(I);
 579     for (int i = 0; i < I->NbRays; ++i) {
 580         if (value_pos_p(I->Ray[i][1]))
 581             *max = 1;
 582         else if (value_neg_p(I->Ray[i][1]))
 583             *min = -1;
 584         else {
 585             if (*max < 0)
 586                 *max = 0;
 587             if (*min > 0)
 588                 *min = 0;
 589         }
 590     }
 591 }
 592
 593 static void interval_minmax(Polyhedron *D, Matrix *T, int *min, int *max,
 594                             unsigned MaxRays)
 595 {
 596     Polyhedron *I = Polyhedron_Image(D, T, MaxRays);
 597     if (MaxRays & POL_INTEGER)
 598         I = DomainConstraintSimplify(I, MaxRays);
 599     if (emptyQ2(I)) {
 600         Polyhedron_Free(I);
 601         I = Polyhedron_Image(D, T, MaxRays);
 602     }
 603     interval_minmax(I, min, max);
 604     Polyhedron_Free(I);
 605 }
 606
 607 struct max_term {
 608     unsigned dim;
 609     EDomain *domain;
 610     vector<evalue *> max;
 611
 612     void print(ostream& os, char **p, barvinok_options *options) const;
 613     void substitute(Matrix *T, unsigned MaxRays);
 614     Vector *eval(Value *val, unsigned MaxRays) const;
 615
 616     ~max_term() {
 617         for (int i = 0; i < max.size(); ++i) {
 618             free_evalue_refs(max[i]);
 619             delete max[i];
 620         }
 621         delete domain;
 622     }
 623 };
 624
 625 /*
 626  * Project on first dim dimensions
 627  */
 628 Polyhedron* Polyhedron_Project_Initial(Polyhedron *P, int dim)
 629 {
 630     int i;
 631     Matrix *T;
 632     Polyhedron *I;
 633
 634     if (P->Dimension == dim)
 635         return Polyhedron_Copy(P);
 636
 637     T = Matrix_Alloc(dim+1, P->Dimension+1);
 638     for (i = 0; i < dim; ++i)
 639         value_set_si(T->p[i][i], 1);
 640     value_set_si(T->p[dim][P->Dimension], 1);
 641     I = Polyhedron_Image(P, T, P->NbConstraints);
 642     Matrix_Free(T);
 643     return I;
 644 }
 645
 646 struct indicator {
 647     vector<indicator_term*> term;
 648     indicator_constructor& ic;
 649     partial_order order;
 650     EDomain      *D;
 651     Polyhedron   *P;
 652     Param_Domain *PD;
 653     barvinok_options    *options;
 654
 655     indicator(indicator_constructor& ic, Param_Domain *PD, EDomain *D,
 656               barvinok_options *options) :
 657         ic(ic), PD(PD), D(D), order(this), options(options), P(NULL) {}
 658     indicator(const indicator& ind, EDomain *D) :
 659             ic(ind.ic), PD(ind.PD), D(NULL), order(this), options(ind.options),
 660             P(Polyhedron_Copy(ind.P)) {
 661         map< indicator_term *, indicator_term * > old2new;
 662         for (int i = 0; i < ind.term.size(); ++i) {
 663             indicator_term *it = new indicator_term(*ind.term[i]);
 664             old2new[ind.term[i]] = it;
 665             term.push_back(it);
 666         }
 667         order.copy(ind.order, old2new);
 668         set_domain(D);
 669     }
 670     ~indicator() {
 671         for (int i = 0; i < term.size(); ++i)
 672             delete term[i];
 673         if (D)
 674             delete D;
 675         if (P)
 676             Polyhedron_Free(P);
 677     }
 678
 679     void set_domain(EDomain *D) {
 680         if (this->D)
 681             delete this->D;
 682         this->D = D;
 683         int nparam = ic.P->Dimension - ic.vertex.length();
 684         Polyhedron *Q = Polyhedron_Project_Initial(D->D, nparam);
 685         Q = DomainConstraintSimplify(Q, options->MaxRays);
 686         if (!P || !PolyhedronIncludes(Q, P))
 687             reduce_in_domain(Q);
 688         if (P)
 689             Polyhedron_Free(P);
 690         P = Q;
 691     }
 692
 693     void add(const indicator_term* it);
 694     void remove(indicator_term* it);
 695     void remove_initial_rational_vertices();
 696     void expand_rational_vertex(indicator_term *initial);
 697
 698     void print(ostream& os, char **p);
 699     void simplify();
 700     void peel(int i, int j);
 701     void combine(indicator_term *a, indicator_term *b);
 702     void substitute(evalue *equation);
 703     void reduce_in_domain(Polyhedron *D);
 704     bool handle_equal_numerators(indicator_term *base);
 705
 706     max_term* create_max_term(indicator_term *it);
 707 };
 708
 709 max_term* indicator::create_max_term(indicator_term *it)
 710 {
 711     int dim = it->den.NumCols();
 712     int nparam = ic.P->Dimension - ic.vertex.length();
 713     max_term *maximum = new max_term;
 714     maximum->dim = nparam;
 715     maximum->domain = new EDomain(D);
 716     for (int j = 0; j < dim; ++j) {
 717         evalue *E = new evalue;
 718         value_init(E->d);
 719         evalue_copy(E, it->vertex[j]);
 720         if (evalue_frac2floor_in_domain(E, D->D))
 721             reduce_evalue(E);
 722         maximum->max.push_back(E);
 723     }
 724     return maximum;
 725 }
 726
 727 /* compute a - b */
 728 static evalue *ediff(const evalue *a, const evalue *b)
 729 {
 730     evalue mone;
 731     value_init(mone.d);
 732     evalue_set_si(&mone, -1, 1);
 733     evalue *diff = new evalue;
 734     value_init(diff->d);
 735     evalue_copy(diff, b);
 736     emul(&mone, diff);
 737     eadd(a, diff);
 738     reduce_evalue(diff);
 739     free_evalue_refs(&mone);
 740     return diff;
 741 }
 742
 743 static order_sign evalue_sign(evalue *diff, EDomain *D, unsigned MaxRays)
 744 {
 745     order_sign sign = order_eq;
 746     evalue mone;
 747     value_init(mone.d);
 748     evalue_set_si(&mone, -1, 1);
 749     int len = 1 + D->D->Dimension + 1;
 750     Vector *c = Vector_Alloc(len);
 751     Matrix *T = Matrix_Alloc(2, len-1);
 752
 753     int fract = evalue2constraint(D, diff, c->p, len);
 754     Vector_Copy(c->p+1, T->p[0], len-1);
 755     value_assign(T->p[1][len-2], c->p[0]);
 756
 757     int min, max;
 758     interval_minmax(D->D, T, &min, &max, MaxRays);
 759     if (max < 0)
 760         sign = order_lt;
 761     else {
 762         if (fract) {
 763             emul(&mone, diff);
 764             evalue2constraint(D, diff, c->p, len);
 765             emul(&mone, diff);
 766             Vector_Copy(c->p+1, T->p[0], len-1);
 767             value_assign(T->p[1][len-2], c->p[0]);
 768
 769             int negmin, negmax;
 770             interval_minmax(D->D, T, &negmin, &negmax, MaxRays);
 771             min = -negmax;
 772         }
 773         if (min > 0)
 774             sign = order_gt;
 775         else if (max == 0 && min == 0)
 776             sign = order_eq;
 777         else if (min < 0 && max > 0)
 778             sign = order_unknown;
 779         else if (min < 0)
 780             sign = order_le;
 781         else
 782             sign = order_ge;
 783     }
 784
 785     Matrix_Free(T);
 786     Vector_Free(c);
 787     free_evalue_refs(&mone);
 788
 789     return sign;
 790 }
 791
 792 order_sign partial_order::compare(indicator_term *a, indicator_term *b)
 793 {
 794     unsigned dim = a->den.NumCols();
 795     order_sign sign = order_eq;
 796     EDomain *D = ind->D;
 797     unsigned MaxRays = ind->options->MaxRays;
 798     if (MaxRays & POL_INTEGER && (a->sign == 0 || b->sign == 0))
 799         MaxRays = 0;
 800
 801     for (int k = 0; k < dim; ++k) {
 802         /* compute a->vertex[k] - b->vertex[k] */
 803         evalue *diff = ediff(a->vertex[k], b->vertex[k]);
 804         order_sign diff_sign = evalue_sign(diff, D, MaxRays);
 805
 806         if (diff_sign == order_lt) {
 807             if (sign == order_eq || sign == order_le)
 808                 sign = order_lt;
 809             else
 810                 sign = order_unknown;
 811             free_evalue_refs(diff);
 812             delete diff;
 813             break;
 814         }
 815         if (diff_sign == order_gt) {
 816             if (sign == order_eq || sign == order_ge)
 817                 sign = order_gt;
 818             else
 819                 sign = order_unknown;
 820             free_evalue_refs(diff);
 821             delete diff;
 822             break;
 823         }
 824         if (diff_sign == order_eq) {
 825             if (D == ind->D && !EVALUE_IS_ZERO(*diff))
 826                 ind->substitute(diff);
 827             free_evalue_refs(diff);
 828             delete diff;
 829             continue;
 830         }
 831         if ((diff_sign == order_unknown) ||
 832             ((diff_sign == order_lt || diff_sign == order_le) && sign == order_ge) ||
 833             ((diff_sign == order_gt || diff_sign == order_ge) && sign == order_le)) {
 834             free_evalue_refs(diff);
 835             delete diff;
 836             sign = order_unknown;
 837             break;
 838         }
 839
 840         sign = diff_sign;
 841
 842         Matrix *M;
 843         vector<EDomain_floor *> new_floors;
 844         M = D->add_ge_constraint(diff, new_floors);
 845         value_set_si(M->p[M->NbRows-1][0], 0);
 846         Polyhedron *D2 = Constraints2Polyhedron(M, MaxRays);
 847         EDomain *EDeq = new EDomain(D2, D, new_floors);
 848         Polyhedron_Free(D2);
 849         Matrix_Free(M);
 850         for (int i = 0; i < new_floors.size(); ++i)
 851             EDomain_floor::unref(new_floors[i]);
 852
 853         if (D != ind->D)
 854             delete D;
 855         D = EDeq;
 856
 857         free_evalue_refs(diff);
 858         delete diff;
 859     }
 860
 861     if (D != ind->D)
 862         delete D;
 863
 864     return sign;
 865 }
 866
 867 bool partial_order::compared(indicator_term* a, indicator_term* b)
 868 {
 869     map<indicator_term *, vector<indicator_term * > >::iterator j;
 870
 871     j = lt.find(a);
 872     if (j != lt.end() && find(lt[a].begin(), lt[a].end(), b) != lt[a].end())
 873         return true;
 874
 875     j = le.find(a);
 876     if (j != le.end() && find(le[a].begin(), le[a].end(), b) != le[a].end())
 877         return true;
 878
 879     return false;
 880 }
 881
 882 void partial_order::add(indicator_term* it, std::set<indicator_term *> *filter)
 883 {
 884     if (eq.find(it) != eq.end() && eq[it].size() == 1)
 885         return;
 886
 887     int it_pred = filter ? pred[it] : 0;
 888
 889     map<indicator_term *, int >::iterator i;
 890     for (i = pred.begin(); i != pred.end(); ++i) {
 891         if ((*i).second != 0)
 892             continue;
 893         if (eq.find((*i).first) != eq.end() && eq[(*i).first].size() == 1)
 894             continue;
 895         if (filter) {
 896             if ((*i).first == it)
 897                 continue;
 898             if (filter->find((*i).first) == filter->end())
 899                 continue;
 900             if (compared((*i).first, it))
 901                 continue;
 902         }
 903         order_sign sign = compare(it, (*i).first);
 904         if (sign == order_lt) {
 905             lt[it].push_back((*i).first);
 906             (*i).second++;
 907         } else if (sign == order_le) {
 908             le[it].push_back((*i).first);
 909             (*i).second++;
 910         } else if (sign == order_eq) {
 911             pred[it] = it_pred;
 912             set_equal(it, (*i).first);
 913             return;
 914         } else if (sign == order_gt) {
 915             pending[(*i).first].push_back(it);
 916             lt[(*i).first].push_back(it);
 917             ++it_pred;
 918         } else if (sign == order_ge) {
 919             pending[(*i).first].push_back(it);
 920             le[(*i).first].push_back(it);
 921             ++it_pred;
 922         }
 923     }
 924     pred[it] = it_pred;
 925 }
 926
 927 void partial_order::remove(indicator_term* it)
 928 {
 929     std::set<indicator_term *> filter;
 930     map<indicator_term *, vector<indicator_term * > >::iterator i;
 931
 932     assert(pred[it] == 0);
 933
 934     i = eq.find(it);
 935     if (i != eq.end()) {
 936         assert(eq[it].size() >= 1);
 937         indicator_term *base;
 938         if (eq[it].size() == 1) {
 939             base = eq[it][0];
 940             eq.erase(i);
 941
 942             vector<indicator_term * >::iterator j;
 943             j = find(eq[base].begin(), eq[base].end(), it);
 944             assert(j != eq[base].end());
 945             eq[base].erase(j);
 946         } else {
 947             /* "it" may no longer be the smallest, since the order
 948              * structure may have been copied from another one.
 949              */
 950             sort(eq[it].begin()+1, eq[it].end());
 951             assert(eq[it][0] == it);
 952             eq[it].erase(eq[it].begin());
 953             base = eq[it][0];
 954             eq[base] = eq[it];
 955             eq.erase(i);
 956
 957             for (int j = 1; j < eq[base].size(); ++j)
 958                 eq[eq[base][j]][0] = base;
 959
 960             i = lt.find(it);
 961             if (i != lt.end()) {
 962                 lt[base] = lt[it];
 963                 lt.erase(i);
 964             }
 965
 966             i = le.find(it);
 967             if (i != le.end()) {
 968                 le[base] = le[it];
 969                 le.erase(i);
 970             }
 971
 972             i = pending.find(it);
 973             if (i != pending.end()) {
 974                 pending[base] = pending[it];
 975                 pending.erase(i);
 976             }
 977         }
 978
 979         if (eq[base].size() == 1)
 980             eq.erase(base);
 981
 982         map<indicator_term *, int >::iterator j;
 983         j = pred.find(it);
 984         pred.erase(j);
 985
 986         return;
 987     }
 988
 989     i = lt.find(it);
 990     if (i != lt.end()) {
 991         for (int j = 0; j < lt[it].size(); ++j) {
 992             filter.insert(lt[it][j]);
 993             pred[lt[it][j]]--;
 994         }
 995         lt.erase(i);
 996     }
 997
 998     i = le.find(it);
 999     if (i != le.end()) {
1000         for (int j = 0; j < le[it].size(); ++j) {
1001             filter.insert(le[it][j]);
1002             pred[le[it][j]]--;
1003         }
1004         le.erase(i);
1005     }
1006
1007     map<indicator_term *, int >::iterator j;
1008     j = pred.find(it);
1009     pred.erase(j);
1010
1011     i = pending.find(it);
1012     if (i != pending.end()) {
1013         for (int j = 0; j < pending[it].size(); ++j) {
1014             filter.erase(pending[it][j]);
1015             add(pending[it][j], &filter);
1016         }
1017         pending.erase(i);
1018     }
1019 }
1020
1021 void partial_order::print(ostream& os, char **p)
1022 {
1023     map<indicator_term *, vector<indicator_term * > >::iterator i;
1024     for (i = lt.begin(); i != lt.end(); ++i) {
1025         (*i).first->print(os, p);
1026         assert(pred.find((*i).first) != pred.end());
1027         os << "(" << pred[(*i).first] << ")";
1028         os << " < ";
1029         for (int j = 0; j < (*i).second.size(); ++j) {
1030             if (j)
1031                 os << ", ";
1032             (*i).second[j]->print(os, p);
1033             assert(pred.find((*i).second[j]) != pred.end());
1034             os << "(" << pred[(*i).second[j]] << ")";
1035         }
1036         os << endl;
1037     }
1038     for (i = le.begin(); i != le.end(); ++i) {
1039         (*i).first->print(os, p);
1040         assert(pred.find((*i).first) != pred.end());
1041         os << "(" << pred[(*i).first] << ")";
1042         os << " <= ";
1043         for (int j = 0; j < (*i).second.size(); ++j) {
1044             if (j)
1045                 os << ", ";
1046             (*i).second[j]->print(os, p);
1047             assert(pred.find((*i).second[j]) != pred.end());
1048             os << "(" << pred[(*i).second[j]] << ")";
1049         }
1050         os << endl;
1051     }
1052     for (i = eq.begin(); i != eq.end(); ++i) {
1053         if ((*i).second.size() <= 1)
1054             continue;
1055         (*i).first->print(os, p);
1056         assert(pred.find((*i).first) != pred.end());
1057         os << "(" << pred[(*i).first] << ")";
1058         for (int j = 1; j < (*i).second.size(); ++j) {
1059             if (j)
1060                 os << " = ";
1061             (*i).second[j]->print(os, p);
1062             assert(pred.find((*i).second[j]) != pred.end());
1063             os << "(" << pred[(*i).second[j]] << ")";
1064         }
1065         os << endl;
1066     }
1067     for (i = pending.begin(); i != pending.end(); ++i) {
1068         os << "pending on ";
1069         (*i).first->print(os, p);
1070         assert(pred.find((*i).first) != pred.end());
1071         os << "(" << pred[(*i).first] << ")";
1072         os << ": ";
1073         for (int j = 0; j < (*i).second.size(); ++j) {
1074             if (j)
1075                 os << ", ";
1076             (*i).second[j]->print(os, p);
1077             assert(pred.find((*i).second[j]) != pred.end());
1078             os << "(" << pred[(*i).second[j]] << ")";
1079         }
1080         os << endl;
1081     }
1082 }
1083
1084 void indicator::add(const indicator_term* it)
1085 {
1086     indicator_term *nt = new indicator_term(*it);
1087     nt->reduce_in_domain(P ? P : D->D);
1088     term.push_back(nt);
1089     order.add(nt, NULL);
1090     assert(term.size() == order.pred.size());
1091 }
1092
1093 void indicator::remove(indicator_term* it)
1094 {
1095     vector<indicator_term *>::iterator i;
1096     i = find(term.begin(), term.end(), it);
1097     assert(i!= term.end());
1098     order.remove(it);
1099     term.erase(i);
1100     assert(term.size() == order.pred.size());
1101     delete it;
1102 }
1103
1104 void indicator::expand_rational_vertex(indicator_term *initial)
1105 {
1106     int pos = initial->pos;
1107     remove(initial);
1108     if (ic.terms[pos].size() == 0) {
1109         Param_Vertices *V;
1110         FORALL_PVertex_in_ParamPolyhedron(V, PD, ic.PP) // _i is internal counter
1111             if (_i == pos) {
1112                 ic.decompose_at_vertex(V, pos, options->MaxRays);
1113                 break;
1114             }
1115         END_FORALL_PVertex_in_ParamPolyhedron;
1116     }
1117     for (int j = 0; j < ic.terms[pos].size(); ++j)
1118         add(ic.terms[pos][j]);
1119 }
1120
1121 void indicator::remove_initial_rational_vertices()
1122 {
1123     do {
1124         indicator_term *initial = NULL;
1125         map<indicator_term *, int >::iterator i;
1126         for (i = order.pred.begin(); i != order.pred.end(); ++i) {
1127             if ((*i).second != 0)
1128                 continue;
1129             if ((*i).first->sign != 0)
1130                 continue;
1131             if (order.eq.find((*i).first) != order.eq.end() &&
1132                 order.eq[(*i).first].size() <= 1)
1133                 continue;
1134             initial = (*i).first;
1135             break;
1136         }
1137         if (!initial)
1138             break;
1139         expand_rational_vertex(initial);
1140     } while(1);
1141 }
1142
1143 void indicator::reduce_in_domain(Polyhedron *D)
1144 {
1145     for (int i = 0; i < term.size(); ++i)
1146         term[i]->reduce_in_domain(D);
1147 }
1148
1149 void indicator::print(ostream& os, char **p)
1150 {
1151     assert(term.size() == order.pred.size());
1152     for (int i = 0; i < term.size(); ++i) {
1153         term[i]->print(os, p);
1154         os << endl;
1155     }
1156     order.print(os, p);
1157 }
1158
1159 /* Remove pairs of opposite terms */
1160 void indicator::simplify()
1161 {
1162     for (int i = 0; i < term.size(); ++i) {
1163         for (int j = i+1; j < term.size(); ++j) {
1164             if (term[i]->sign + term[j]->sign != 0)
1165                 continue;
1166             if (term[i]->den != term[j]->den)
1167                 continue;
1168             int k;
1169             for (k = 0; k < term[i]->den.NumCols(); ++k)
1170                 if (!eequal(term[i]->vertex[k], term[j]->vertex[k]))
1171                     break;
1172             if (k < term[i]->den.NumCols())
1173                 continue;
1174             delete term[i];
1175             delete term[j];
1176             term.erase(term.begin()+j);
1177             term.erase(term.begin()+i);
1178             --i;
1179             break;
1180         }
1181     }
1182 }
1183
1184 void indicator::peel(int i, int j)
1185 {
1186     if (j < i) {
1187         int tmp = i;
1188         i = j;
1189         j = tmp;
1190     }
1191     assert (i < j);
1192     int dim = term[i]->den.NumCols();
1193
1194     mat_ZZ common;
1195     mat_ZZ rest_i;
1196     mat_ZZ rest_j;
1197     int n_common = 0, n_i = 0, n_j = 0;
1198
1199     common.SetDims(min(term[i]->den.NumRows(), term[j]->den.NumRows()), dim);
1200     rest_i.SetDims(term[i]->den.NumRows(), dim);
1201     rest_j.SetDims(term[j]->den.NumRows(), dim);
1202
1203     int k, l;
1204     for (k = 0, l = 0; k < term[i]->den.NumRows() && l < term[j]->den.NumRows(); ) {
1205         int s = lex_cmp(term[i]->den[k], term[j]->den[l]);
1206         if (s == 0) {
1207             common[n_common++] = term[i]->den[k];
1208             ++k;
1209             ++l;
1210         } else if (s < 0)
1211             rest_i[n_i++] = term[i]->den[k++];
1212         else
1213             rest_j[n_j++] = term[j]->den[l++];
1214     }
1215     while (k < term[i]->den.NumRows())
1216         rest_i[n_i++] = term[i]->den[k++];
1217     while (l < term[j]->den.NumRows())
1218         rest_j[n_j++] = term[j]->den[l++];
1219     common.SetDims(n_common, dim);
1220     rest_i.SetDims(n_i, dim);
1221     rest_j.SetDims(n_j, dim);
1222
1223     for (k = 0; k <= n_i; ++k) {
1224         indicator_term *it = new indicator_term(*term[i]);
1225         it->den.SetDims(n_common + k, dim);
1226         for (l = 0; l < n_common; ++l)
1227             it->den[l] = common[l];
1228         for (l = 0; l < k; ++l)
1229             it->den[n_common+l] = rest_i[l];
1230         lex_order_rows(it->den);
1231         if (k)
1232             for (l = 0; l < dim; ++l)
1233                 evalue_add_constant(it->vertex[l], rest_i[k-1][l]);
1234         term.push_back(it);
1235     }
1236
1237     for (k = 0; k <= n_j; ++k) {
1238         indicator_term *it = new indicator_term(*term[j]);
1239         it->den.SetDims(n_common + k, dim);
1240         for (l = 0; l < n_common; ++l)
1241             it->den[l] = common[l];
1242         for (l = 0; l < k; ++l)
1243             it->den[n_common+l] = rest_j[l];
1244         lex_order_rows(it->den);
1245         if (k)
1246             for (l = 0; l < dim; ++l)
1247                 evalue_add_constant(it->vertex[l], rest_j[k-1][l]);
1248         term.push_back(it);
1249     }
1250     term.erase(term.begin()+j);
1251     term.erase(term.begin()+i);
1252 }
1253
1254 void indicator::combine(indicator_term *a, indicator_term *b)
1255 {
1256     int dim = a->den.NumCols();
1257
1258     mat_ZZ common;
1259     mat_ZZ rest_i;
1260     mat_ZZ rest_j;
1261     int n_common = 0, n_i = 0, n_j = 0;
1262
1263     common.SetDims(min(a->den.NumRows(), b->den.NumRows()), dim);
1264     rest_i.SetDims(a->den.NumRows(), dim);
1265     rest_j.SetDims(b->den.NumRows(), dim);
1266
1267     int k, l;
1268     for (k = 0, l = 0; k < a->den.NumRows() && l < b->den.NumRows(); ) {
1269         int s = lex_cmp(a->den[k], b->den[l]);
1270         if (s == 0) {
1271             common[n_common++] = a->den[k];
1272             ++k;
1273             ++l;
1274         } else if (s < 0)
1275             rest_i[n_i++] = a->den[k++];
1276         else
1277             rest_j[n_j++] = b->den[l++];
1278     }
1279     while (k < a->den.NumRows())
1280         rest_i[n_i++] = a->den[k++];
1281     while (l < b->den.NumRows())
1282         rest_j[n_j++] = b->den[l++];
1283     common.SetDims(n_common, dim);
1284     rest_i.SetDims(n_i, dim);
1285     rest_j.SetDims(n_j, dim);
1286
1287     assert(n_i < 30);
1288     assert(n_j < 30);
1289
1290     for (k = (1 << n_i)-1; k >= 0; --k) {
1291         indicator_term *it = k ? new indicator_term(*b) : b;
1292         it->den.SetDims(n_common + n_i + n_j, dim);
1293         for (l = 0; l < n_common; ++l)
1294             it->den[l] = common[l];
1295         for (l = 0; l < n_i; ++l)
1296             it->den[n_common+l] = rest_i[l];
1297         for (l = 0; l < n_j; ++l)
1298             it->den[n_common+n_i+l] = rest_j[l];
1299         lex_order_rows(it->den);
1300         int change = 0;
1301         for (l = 0; l < n_i; ++l) {
1302             if (!(k & (1 <<l)))
1303                 continue;
1304             change ^= 1;
1305             for (int m = 0; m < dim; ++m)
1306                 evalue_add_constant(it->vertex[m], rest_i[l][m]);
1307         }
1308         if (change)
1309             it->sign = -it->sign;
1310         if (it != b) {
1311             term.push_back(it);
1312             order.add(it, NULL);
1313         }
1314     }
1315
1316     for (k = (1 << n_j)-1; k >= 0; --k) {
1317         indicator_term *it = k ? new indicator_term(*a) : a;
1318         it->den.SetDims(n_common + n_i + n_j, dim);
1319         for (l = 0; l < n_common; ++l)
1320             it->den[l] = common[l];
1321         for (l = 0; l < n_i; ++l)
1322             it->den[n_common+l] = rest_i[l];
1323         for (l = 0; l < n_j; ++l)
1324             it->den[n_common+n_i+l] = rest_j[l];
1325         lex_order_rows(it->den);
1326         int change = 0;
1327         for (l = 0; l < n_j; ++l) {
1328             if (!(k & (1 <<l)))
1329                 continue;
1330             change ^= 1;
1331             for (int m = 0; m < dim; ++m)
1332                 evalue_add_constant(it->vertex[m], rest_j[l][m]);
1333         }
1334         if (change)
1335             it->sign = -it->sign;
1336         if (it != a) {
1337             term.push_back(it);
1338             order.add(it, NULL);
1339         }
1340     }
1341 }
1342
1343 bool indicator::handle_equal_numerators(indicator_term *base)
1344 {
1345     for (int i = 0; i < order.eq[base].size(); ++i) {
1346         for (int j = i+1; j < order.eq[base].size(); ++j) {
1347             if (order.eq[base][i]->is_opposite(order.eq[base][j])) {
1348                 remove(order.eq[base][j]);
1349                 remove(i ? order.eq[base][i] : base);
1350                 return true;
1351             }
1352         }
1353     }
1354     for (int j = 1; j < order.eq[base].size(); ++j)
1355         if (order.eq[base][j]->sign != base->sign) {
1356             combine(base, order.eq[base][j]);
1357             return true;
1358         }
1359     return false;
1360 }
1361
1362 void indicator::substitute(evalue *equation)
1363 {
1364     evalue *fract = NULL;
1365     evalue *val = new evalue;
1366     value_init(val->d);
1367     evalue_copy(val, equation);
1368
1369     evalue factor;
1370     value_init(factor.d);
1371     value_init(factor.x.n);
1372
1373     evalue *e;
1374     for (e = val; value_zero_p(e->d) && e->x.p->type != fractional;
1375          e = &e->x.p->arr[type_offset(e->x.p)])
1376         ;
1377
1378     if (value_zero_p(e->d) && e->x.p->type == fractional)
1379         fract = &e->x.p->arr[0];
1380     else {
1381         for (e = val; value_zero_p(e->d) && e->x.p->type != polynomial;
1382              e = &e->x.p->arr[type_offset(e->x.p)])
1383             ;
1384         assert(value_zero_p(e->d) && e->x.p->type == polynomial);
1385     }
1386
1387     int offset = type_offset(e->x.p);
1388
1389     assert(value_notzero_p(e->x.p->arr[offset+1].d));
1390     assert(value_notzero_p(e->x.p->arr[offset+1].x.n));
1391     if (value_neg_p(e->x.p->arr[offset+1].x.n)) {
1392         value_oppose(factor.d, e->x.p->arr[offset+1].x.n);
1393         value_assign(factor.x.n, e->x.p->arr[offset+1].d);
1394     } else {
1395         value_assign(factor.d, e->x.p->arr[offset+1].x.n);
1396         value_oppose(factor.x.n, e->x.p->arr[offset+1].d);
1397     }
1398
1399     free_evalue_refs(&e->x.p->arr[offset+1]);
1400     enode *p = e->x.p;
1401     value_clear(e->d);
1402     *e = e->x.p->arr[offset];
1403
1404     emul(&factor, val);
1405
1406     if (fract) {
1407         for (int i = 0; i < term.size(); ++i)
1408             term[i]->substitute(fract, val);
1409
1410         free_evalue_refs(&p->arr[0]);
1411     } else {
1412         for (int i = 0; i < term.size(); ++i)
1413             term[i]->substitute(p->pos, val);
1414     }
1415
1416     free_evalue_refs(&factor);
1417     free_evalue_refs(val);
1418     delete val;
1419
1420     free(p);
1421 }
1422
1423 static void print_varlist(ostream& os, int n, char **names)
1424 {
1425     int i;
1426     os << "[";
1427     for (i = 0; i < n; ++i) {
1428         if (i)
1429             os << ",";
1430         os << names[i];
1431     }
1432     os << "]";
1433 }
1434
1435 void max_term::print(ostream& os, char **p, barvinok_options *options) const
1436 {
1437     os << "{ ";
1438     print_varlist(os, dim, p);
1439     os << " -> ";
1440     os << "[";
1441     for (int i = 0; i < max.size(); ++i) {
1442         if (i)
1443             os << ",";
1444         evalue_print(os, max[i], p);
1445     }
1446     os << "]";
1447     os << " : ";
1448     domain->print_constraints(os, p, options);
1449     os << " }" << endl;
1450 }
1451
1452 static void evalue_substitute(evalue *e, evalue **subs)
1453 {
1454     evalue *v;
1455
1456     if (value_notzero_p(e->d))
1457         return;
1458
1459     enode *p = e->x.p;
1460     for (int i = 0; i < p->size; ++i)
1461         evalue_substitute(&p->arr[i], subs);
1462
1463     if (p->type == polynomial)
1464         v = subs[p->pos-1];
1465     else {
1466         v = new evalue;
1467         value_init(v->d);
1468         value_set_si(v->d, 0);
1469         v->x.p = new_enode(p->type, 3, -1);
1470         value_clear(v->x.p->arr[0].d);
1471         v->x.p->arr[0] = p->arr[0];
1472         evalue_set_si(&v->x.p->arr[1], 0, 1);
1473         evalue_set_si(&v->x.p->arr[2], 1, 1);
1474     }
1475
1476     int offset = type_offset(p);
1477
1478     for (int i = p->size-1; i >= offset+1; i--) {
1479         emul(v, &p->arr[i]);
1480         eadd(&p->arr[i], &p->arr[i-1]);
1481         free_evalue_refs(&(p->arr[i]));
1482     }
1483
1484     if (p->type != polynomial) {
1485         free_evalue_refs(v);
1486         delete v;
1487     }
1488
1489     value_clear(e->d);
1490     *e = p->arr[offset];
1491     free(p);
1492 }
1493
1494 /* "align" matrix to have nrows by inserting
1495  * the necessary number of rows and an equal number of columns at the end
1496  * right before the constant row/column
1497  */
1498 static Matrix *align_matrix_initial(Matrix *M, int nrows)
1499 {
1500     int i;
1501     int newrows = nrows - M->NbRows;
1502     Matrix *M2 = Matrix_Alloc(nrows, newrows + M->NbColumns);
1503     for (i = 0; i < newrows; ++i)
1504         value_set_si(M2->p[M->NbRows-1+i][M->NbColumns-1+i], 1);
1505     for (i = 0; i < M->NbRows-1; ++i) {
1506         Vector_Copy(M->p[i], M2->p[i], M->NbColumns-1);
1507         value_assign(M2->p[i][M2->NbColumns-1], M->p[i][M->NbColumns-1]);
1508     }
1509     value_assign(M2->p[M2->NbRows-1][M2->NbColumns-1],
1510                  M->p[M->NbRows-1][M->NbColumns-1]);
1511     return M2;
1512 }
1513
1514 /* T maps the compressed parameters to the original parameters,
1515  * while this max_term is based on the compressed parameters
1516  * and we want get the original parameters back.
1517  */
1518 void max_term::substitute(Matrix *T, unsigned MaxRays)
1519 {
1520     int nexist = domain->D->Dimension - (T->NbColumns-1);
1521     Matrix *M = align_matrix_initial(T, T->NbRows+nexist);
1522
1523     Polyhedron *D = DomainImage(domain->D, M, MaxRays);
1524     Polyhedron_Free(domain->D);
1525     domain->D = D;
1526     Matrix_Free(M);
1527
1528     assert(T->NbRows == T->NbColumns);
1529     Matrix *T2 = Matrix_Copy(T);
1530     Matrix *inv = Matrix_Alloc(T->NbColumns, T->NbRows);
1531     int ok = Matrix_Inverse(T2, inv);
1532     Matrix_Free(T2);
1533     assert(ok);
1534
1535     evalue denom;
1536     value_init(denom.d);
1537     value_init(denom.x.n);
1538     value_set_si(denom.x.n, 1);
1539     value_assign(denom.d, inv->p[inv->NbRows-1][inv->NbColumns-1]);
1540
1541     vec_ZZ v;
1542     v.SetLength(inv->NbColumns);
1543     evalue* subs[inv->NbRows-1];
1544     for (int i = 0; i < inv->NbRows-1; ++i) {
1545         values2zz(inv->p[i], v, v.length());
1546         subs[i] = multi_monom(v);
1547         emul(&denom, subs[i]);
1548     }
1549     free_evalue_refs(&denom);
1550
1551     for (int i = 0; i < max.size(); ++i) {
1552         evalue_substitute(max[i], subs);
1553         reduce_evalue(max[i]);
1554     }
1555
1556     for (int i = 0; i < inv->NbRows-1; ++i) {
1557         free_evalue_refs(subs[i]);
1558         delete subs[i];
1559     }
1560     Matrix_Free(inv);
1561 }
1562
1563 int Last_Non_Zero(Value *p, unsigned len)
1564 {
1565     for (int i = len-1; i >= 0; --i)
1566         if (value_notzero_p(p[i]))
1567             return i;
1568     return -1;
1569 }
1570
1571 static void SwapColumns(Value **V, int n, int i, int j)
1572 {
1573     for (int r = 0; r < n; ++r)
1574         value_swap(V[r][i], V[r][j]);
1575 }
1576
1577 static void SwapColumns(Polyhedron *P, int i, int j)
1578 {
1579     SwapColumns(P->Constraint, P->NbConstraints, i, j);
1580     SwapColumns(P->Ray, P->NbRays, i, j);
1581 }
1582
1583 void compute_evalue(evalue *e, Value *val, Value *res)
1584 {
1585     double d = compute_evalue(e, val);
1586     if (d > 0)
1587         d += .25;
1588     else
1589         d -= .25;
1590     value_set_double(*res, d);
1591 }
1592
1593 Vector *max_term::eval(Value *val, unsigned MaxRays) const
1594 {
1595     if (!domain->contains(val, dim))
1596         return NULL;
1597     Vector *res = Vector_Alloc(max.size());
1598     for (int i = 0; i < max.size(); ++i) {
1599         compute_evalue(max[i], val, &res->p[i]);
1600     }
1601     return res;
1602 }
1603
1604 static Matrix *remove_equalities(Polyhedron **P, unsigned nparam, unsigned MaxRays);
1605
1606 Vector *Polyhedron_not_empty(Polyhedron *P, unsigned MaxRays)
1607 {
1608     Polyhedron *Porig = P;
1609     Vector *sample = NULL;
1610
1611     POL_ENSURE_VERTICES(P);
1612     if (emptyQ2(P))
1613         return NULL;
1614
1615     for (int i = 0; i < P->NbRays; ++i)
1616         if (value_one_p(P->Ray[i][1+P->Dimension])) {
1617             sample = Vector_Alloc(P->Dimension + 1);
1618             Vector_Copy(P->Ray[i]+1, sample->p, P->Dimension+1);
1619             return sample;
1620         }
1621
1622     Matrix *T = remove_equalities(&P, 0, MaxRays);
1623     if (P)
1624         sample = Polyhedron_Sample(P, MaxRays);
1625     if (sample) {
1626         if (T) {
1627             Vector *P_sample = Vector_Alloc(Porig->Dimension + 1);
1628             Matrix_Vector_Product(T, sample->p, P_sample->p);
1629             Vector_Free(sample);
1630             sample = P_sample;
1631         }
1632     }
1633     if (T) {
1634         Polyhedron_Free(P);
1635         Matrix_Free(T);
1636     }
1637
1638     return sample;
1639 }
1640
1641 struct split {
1642     evalue *constraint;
1643     enum sign { le, ge, lge } sign;
1644
1645     split (evalue *c, enum sign s) : constraint(c), sign(s) {}
1646 };
1647
1648 static void split_on(const split& sp, EDomain *D,
1649                      EDomain **Dlt, EDomain **Deq, EDomain **Dgt,
1650                      barvinok_options *options)
1651 {
1652     Matrix *M, *M2;
1653     EDomain *ED[3];
1654     Polyhedron *D2;
1655     Value mone;
1656     value_init(mone);
1657     value_set_si(mone, -1);
1658     *Dlt = NULL;
1659     *Deq = NULL;
1660     *Dgt = NULL;
1661     vector<EDomain_floor *> new_floors;
1662     M = D->add_ge_constraint(sp.constraint, new_floors);
1663     if (sp.sign == split::lge || sp.sign == split::ge) {
1664         M2 = Matrix_Copy(M);
1665         value_decrement(M2->p[M2->NbRows-1][M2->NbColumns-1],
1666                         M2->p[M2->NbRows-1][M2->NbColumns-1]);
1667         D2 = Constraints2Polyhedron(M2, options->MaxRays);
1668         ED[2] = new EDomain(D2, D, new_floors);
1669         Polyhedron_Free(D2);
1670         Matrix_Free(M2);
1671     } else
1672         ED[2] = NULL;
1673     if (sp.sign == split::lge || sp.sign == split::le) {
1674         M2 = Matrix_Copy(M);
1675         Vector_Scale(M2->p[M2->NbRows-1]+1, M2->p[M2->NbRows-1]+1,
1676                      mone, M2->NbColumns-1);
1677         value_decrement(M2->p[M2->NbRows-1][M2->NbColumns-1],
1678                         M2->p[M2->NbRows-1][M2->NbColumns-1]);
1679         D2 = Constraints2Polyhedron(M2, options->MaxRays);
1680         ED[0] = new EDomain(D2, D, new_floors);
1681         Polyhedron_Free(D2);
1682         Matrix_Free(M2);
1683     } else
1684         ED[0] = NULL;
1685
1686     assert(sp.sign == split::lge || sp.sign == split::ge || sp.sign == split::le);
1687     value_set_si(M->p[M->NbRows-1][0], 0);
1688     D2 = Constraints2Polyhedron(M, options->MaxRays);
1689     ED[1] = new EDomain(D2, D, new_floors);
1690     Polyhedron_Free(D2);
1691     Matrix_Free(M);
1692
1693     Vector *sample = D->sample;
1694     if (sample && new_floors.size() > 0) {
1695         assert(sample->Size == D->D->Dimension+1);
1696         sample = Vector_Alloc(D->D->Dimension+new_floors.size()+1);
1697         Vector_Copy(D->sample->p, sample->p, D->D->Dimension);
1698         value_set_si(sample->p[D->D->Dimension+new_floors.size()], 1);
1699         for (int i = 0; i < new_floors.size(); ++i)
1700             new_floors[i]->eval(sample->p, &sample->p[D->D->Dimension+i]);
1701     }
1702
1703     for (int i = 0; i < new_floors.size(); ++i)
1704         EDomain_floor::unref(new_floors[i]);
1705
1706     for (int i = 0; i < 3; ++i) {
1707         if (!ED[i])
1708             continue;
1709         if (sample && ED[i]->contains(sample->p, sample->Size-1)) {
1710             ED[i]->sample = Vector_Alloc(sample->Size);
1711             Vector_Copy(sample->p, ED[i]->sample->p, sample->Size);
1712         } else if (emptyQ2(ED[i]->D) ||
1713                     (options->emptiness_check == 1 &&
1714                      !(ED[i]->sample = Polyhedron_not_empty(ED[i]->D,
1715                                                             options->MaxRays)))) {
1716             delete ED[i];
1717             ED[i] = NULL;
1718         }
1719     }
1720     *Dlt = ED[0];
1721     *Deq = ED[1];
1722     *Dgt = ED[2];
1723     value_clear(mone);
1724     if (sample != D->sample)
1725         Vector_Free(sample);
1726 }
1727
1728 ostream & operator<< (ostream & os, const vector<int> & v)
1729 {
1730     os << "[";
1731     for (int i = 0; i < v.size(); ++i) {
1732         if (i)
1733             os << ", ";
1734         os << v[i];
1735     }
1736     os << "]";
1737     return os;
1738 }
1739
1740 static bool isTranslation(Matrix *M)
1741 {
1742     unsigned i, j;
1743     if (M->NbRows != M->NbColumns)
1744         return False;
1745
1746     for (i = 0;i < M->NbRows; i ++)
1747         for (j = 0; j < M->NbColumns-1; j ++)
1748             if (i == j) {
1749                 if(value_notone_p(M->p[i][j]))
1750                     return False;
1751             } else {
1752                 if(value_notzero_p(M->p[i][j]))
1753                     return False;
1754             }
1755     return value_one_p(M->p[M->NbRows-1][M->NbColumns-1]);
1756 }
1757
1758 static Matrix *compress_parameters(Polyhedron **P, Polyhedron **C,
1759                                    unsigned nparam, unsigned MaxRays)
1760 {
1761     Matrix *M, *T, *CP;
1762
1763     /* compress_parms doesn't like equalities that only involve parameters */
1764     for (int i = 0; i < (*P)->NbEq; ++i)
1765         assert(First_Non_Zero((*P)->Constraint[i]+1, (*P)->Dimension-nparam) != -1);
1766
1767     M = Matrix_Alloc((*P)->NbEq, (*P)->Dimension+2);
1768     Vector_Copy((*P)->Constraint[0], M->p[0], (*P)->NbEq * ((*P)->Dimension+2));
1769     CP = compress_parms(M, nparam);
1770     Matrix_Free(M);
1771
1772     if (isTranslation(CP)) {
1773         Matrix_Free(CP);
1774         return NULL;
1775     }
1776
1777     T = align_matrix(CP, (*P)->Dimension+1);
1778     *P = Polyhedron_Preimage(*P, T, MaxRays);
1779     Matrix_Free(T);
1780
1781     *C = Polyhedron_Preimage(*C, CP, MaxRays);
1782
1783     return CP;
1784 }
1785
1786 static Matrix *remove_equalities(Polyhedron **P, unsigned nparam, unsigned MaxRays)
1787 {
1788     /* Matrix "view" of equalities */
1789     Matrix M;
1790     M.NbRows = (*P)->NbEq;
1791     M.NbColumns = (*P)->Dimension+2;
1792     M.p_Init = (*P)->p_Init;
1793     M.p = (*P)->Constraint;
1794
1795     Matrix *T = compress_variables(&M, nparam);
1796
1797     if (!T) {
1798         *P = NULL;
1799         return NULL;
1800     }
1801     if (isIdentity(T)) {
1802         Matrix_Free(T);
1803         T = NULL;
1804     } else
1805         *P = Polyhedron_Preimage(*P, T, MaxRays);
1806
1807     return T;
1808 }
1809
1810 void construct_rational_vertices(Param_Polyhedron *PP, Matrix *T, unsigned dim,
1811                                  int nparam, vector<indicator_term *>& vertices)
1812 {
1813     int i;
1814     Param_Vertices *PV;
1815     Value lcm, tmp;
1816     value_init(lcm);
1817     value_init(tmp);
1818
1819     vec_ZZ v;
1820     v.SetLength(nparam+1);
1821
1822     evalue factor;
1823     value_init(factor.d);
1824     value_init(factor.x.n);
1825     value_set_si(factor.x.n, 1);
1826     value_set_si(factor.d, 1);
1827
1828     for (i = 0, PV = PP->V; PV; ++i, PV = PV->next) {
1829         indicator_term *term = new indicator_term(dim, i);
1830         vertices.push_back(term);
1831         Matrix *M = Matrix_Alloc(PV->Vertex->NbRows+nparam+1, nparam+1);
1832         value_set_si(lcm, 1);
1833         for (int j = 0; j < PV->Vertex->NbRows; ++j)
1834             value_lcm(lcm, PV->Vertex->p[j][nparam+1], &lcm);
1835         value_assign(M->p[M->NbRows-1][M->NbColumns-1], lcm);
1836         for (int j = 0; j < PV->Vertex->NbRows; ++j) {
1837             value_division(tmp, lcm, PV->Vertex->p[j][nparam+1]);
1838             Vector_Scale(PV->Vertex->p[j], M->p[j], tmp, nparam+1);
1839         }
1840         for (int j = 0; j < nparam; ++j)
1841             value_assign(M->p[PV->Vertex->NbRows+j][j], lcm);
1842         if (T) {
1843             Matrix *M2 = Matrix_Alloc(T->NbRows, M->NbColumns);
1844             Matrix_Product(T, M, M2);
1845             Matrix_Free(M);
1846             M = M2;
1847         }
1848         for (int j = 0; j < dim; ++j) {
1849             values2zz(M->p[j], v, nparam+1);
1850             term->vertex[j] = multi_monom(v);
1851             value_assign(factor.d, lcm);
1852             emul(&factor, term->vertex[j]);
1853         }
1854         Matrix_Free(M);
1855     }
1856     assert(i == PP->nbV);
1857     free_evalue_refs(&factor);
1858     value_clear(lcm);
1859     value_clear(tmp);
1860 }
1861
1862 /* An auxiliary class that keeps a reference to an evalue
1863  * and frees it when it goes out of scope.
1864  */
1865 struct temp_evalue {
1866     evalue *E;
1867     temp_evalue() : E(NULL) {}
1868     temp_evalue(evalue *e) : E(e) {}
1869     operator evalue* () const { return E; }
1870     evalue *operator=(evalue *e) {
1871         if (E) {
1872             free_evalue_refs(E);
1873             delete E;
1874         }
1875         E = e;
1876         return E;
1877     }
1878     ~temp_evalue() {
1879         if (E) {
1880             free_evalue_refs(E);
1881             delete E;
1882         }
1883     }
1884 };
1885
1886 static vector<max_term*> lexmin(indicator& ind, unsigned nparam,
1887                                     vector<int> loc)
1888 {
1889     vector<max_term*> maxima;
1890     map<indicator_term *, int >::iterator i;
1891     vector<int> best_score;
1892     vector<int> second_score;
1893     vector<int> neg_score;
1894
1895     do {
1896         indicator_term *best = NULL, *second = NULL, *neg = NULL,
1897                        *neg_eq = NULL, *neg_le = NULL;
1898         for (i = ind.order.pred.begin(); i != ind.order.pred.end(); ++i) {
1899             vector<int> score;
1900             if ((*i).second != 0)
1901                 continue;
1902             indicator_term *term = (*i).first;
1903             if (term->sign == 0) {
1904                 ind.expand_rational_vertex(term);
1905                 break;
1906             }
1907
1908             if (ind.order.eq.find(term) != ind.order.eq.end()) {
1909                 int j;
1910                 if (ind.order.eq[term].size() <= 1)
1911                     continue;
1912                 for (j = 1; j < ind.order.eq[term].size(); ++j)
1913                     if (ind.order.pred[ind.order.eq[term][j]] != 0)
1914                         break;
1915                 if (j < ind.order.eq[term].size())
1916                     continue;
1917                 score.push_back(ind.order.eq[term].size());
1918             } else
1919                 score.push_back(0);
1920             if (ind.order.le.find(term) != ind.order.le.end())
1921                 score.push_back(ind.order.le[term].size());
1922             else
1923                 score.push_back(0);
1924             if (ind.order.lt.find(term) != ind.order.lt.end())
1925                 score.push_back(-ind.order.lt[term].size());
1926             else
1927                 score.push_back(0);
1928
1929             if (term->sign > 0) {
1930                 if (!best || score < best_score) {
1931                     second = best;
1932                     second_score = best_score;
1933                     best = term;
1934                     best_score = score;
1935                 } else if (!second || score < second_score) {
1936                     second = term;
1937                     second_score = score;
1938                 }
1939             } else {
1940                 if (!neg_eq && ind.order.eq.find(term) != ind.order.eq.end()) {
1941                     for (int j = 1; j < ind.order.eq[term].size(); ++j)
1942                         if (ind.order.eq[term][j]->sign != term->sign) {
1943                             neg_eq = term;
1944                             break;
1945                         }
1946                 }
1947                 if (!neg_le && ind.order.le.find(term) != ind.order.le.end())
1948                     neg_le = term;
1949                 if (!neg || score < neg_score) {
1950                     neg = term;
1951                     neg_score = score;
1952                 }
1953             }
1954         }
1955         if (i != ind.order.pred.end())
1956             continue;
1957
1958         if (!best && neg_eq) {
1959             assert(ind.order.eq[neg_eq].size() != 0);
1960             bool handled = ind.handle_equal_numerators(neg_eq);
1961             assert(handled);
1962             continue;
1963         }
1964
1965         if (!best && neg_le) {
1966             /* The smallest term is negative and <= some positive term */
1967             best = neg_le;
1968             neg = NULL;
1969         }
1970
1971         if (!best) {
1972             /* apparently there can be negative initial term on empty domains */
1973             if (ind.options->emptiness_check == 1)
1974                 assert(!neg);
1975             break;
1976         }
1977
1978         if (!second && !neg) {
1979             indicator_term *rat = NULL;
1980             assert(best);
1981             if (ind.order.le[best].size() == 0) {
1982                 if (ind.order.eq[best].size() != 0) {
1983                     bool handled = ind.handle_equal_numerators(best);
1984                     if (ind.options->emptiness_check == 1)
1985                         assert(handled);
1986                     /* If !handled then the leading coefficient is bigger than one;
1987                      * must be an empty domain
1988                      */
1989                     if (handled)
1990                         continue;
1991                     else
1992                         break;
1993                 }
1994                 maxima.push_back(ind.create_max_term(best));
1995                 break;
1996             }
1997             for (int j = 0; j < ind.order.le[best].size(); ++j) {
1998                 if (ind.order.le[best][j]->sign == 0) {
1999                     if (!rat && ind.order.pred[ind.order.le[best][j]] == 1)
2000                         rat = ind.order.le[best][j];
2001                 } else if (ind.order.le[best][j]->sign > 0) {
2002                     second = ind.order.le[best][j];
2003                     break;
2004                 } else if (!neg)
2005                     neg = ind.order.le[best][j];
2006             }
2007
2008             if (!second && !neg) {
2009                 assert(rat);
2010                 ind.order.unset_le(best, rat);
2011                 ind.expand_rational_vertex(rat);
2012                 continue;
2013             }
2014
2015             if (!second)
2016                 second = neg;
2017
2018             ind.order.unset_le(best, second);
2019         }
2020
2021         if (!second)
2022             second = neg;
2023
2024         unsigned dim = best->den.NumCols();
2025         temp_evalue diff;
2026         order_sign sign;
2027         for (int k = 0; k < dim; ++k) {
2028             diff = ediff(best->vertex[k], second->vertex[k]);
2029             sign = evalue_sign(diff, ind.D, ind.options->MaxRays);
2030
2031             /* neg can never be smaller than best, unless it may still cancel.
2032              * This can happen if positive terms have been determined to
2033              * be equal or less than or equal to some negative term.
2034              */
2035             if (second == neg && !neg_eq && !neg_le) {
2036                 if (sign == order_ge)
2037                     sign = order_eq;
2038                 if (sign == order_unknown)
2039                     sign = order_le;
2040             }
2041
2042             if (sign != order_eq)
2043                 break;
2044             if (!EVALUE_IS_ZERO(*diff))
2045                 ind.substitute(diff);
2046         }
2047         if (sign == order_eq) {
2048             ind.order.set_equal(best, second);
2049             continue;
2050         }
2051         if (sign == order_lt) {
2052             ind.order.lt[best].push_back(second);
2053             ind.order.pred[second]++;
2054             continue;
2055         }
2056         if (sign == order_gt) {
2057             ind.order.lt[second].push_back(best);
2058             ind.order.pred[best]++;
2059             continue;
2060         }
2061
2062         split sp(diff, sign == order_le ? split::le :
2063                        sign == order_ge ? split::ge : split::lge);
2064
2065         EDomain *Dlt, *Deq, *Dgt;
2066         split_on(sp, ind.D, &Dlt, &Deq, &Dgt, ind.options);
2067         if (ind.options->emptiness_check == 1)
2068             assert(Dlt || Deq || Dgt);
2069         else if (!(Dlt || Deq || Dgt))
2070             /* Must have been empty all along */
2071             break;
2072
2073         if (Deq && (Dlt || Dgt)) {
2074             int locsize = loc.size();
2075             loc.push_back(0);
2076             indicator indeq(ind, Deq);
2077             Deq = NULL;
2078             indeq.substitute(diff);
2079             vector<max_term*> maxeq = lexmin(indeq, nparam, loc);
2080             maxima.insert(maxima.end(), maxeq.begin(), maxeq.end());
2081             loc.resize(locsize);
2082         }
2083         if (Dgt && Dlt) {
2084             int locsize = loc.size();
2085             loc.push_back(1);
2086             indicator indgt(ind, Dgt);
2087             Dgt = NULL;
2088             /* we don't know the new location of these terms in indgt */
2089             /*
2090             indgt.order.lt[second].push_back(best);
2091             indgt.order.pred[best]++;
2092             */
2093             vector<max_term*> maxgt = lexmin(indgt, nparam, loc);
2094             maxima.insert(maxima.end(), maxgt.begin(), maxgt.end());
2095             loc.resize(locsize);
2096         }
2097
2098         if (Deq) {
2099             loc.push_back(0);
2100             ind.substitute(diff);
2101             ind.set_domain(Deq);
2102         }
2103         if (Dlt) {
2104             loc.push_back(-1);
2105             ind.order.lt[best].push_back(second);
2106             ind.order.pred[second]++;
2107             ind.set_domain(Dlt);
2108         }
2109         if (Dgt) {
2110             loc.push_back(1);
2111             ind.order.lt[second].push_back(best);
2112             ind.order.pred[best]++;
2113             ind.set_domain(Dgt);
2114         }
2115     } while(1);
2116
2117     return maxima;
2118 }
2119
2120 static vector<max_term*> lexmin(Polyhedron *P, Polyhedron *C,
2121                                 barvinok_options *options)
2122 {
2123     unsigned nparam = C->Dimension;
2124     Param_Polyhedron *PP = NULL;
2125     Polyhedron *CEq = NULL, *pVD;
2126     Matrix *CT = NULL;
2127     Matrix *T = NULL, *CP = NULL;
2128     Param_Domain *D, *next;
2129     Param_Vertices *V;
2130     Polyhedron *Porig = P;
2131     Polyhedron *Corig = C;
2132     vector<max_term*> all_max;
2133     Polyhedron *Q;
2134     unsigned P2PSD_MaxRays;
2135
2136     if (emptyQ2(P))
2137         return all_max;
2138
2139     POL_ENSURE_VERTICES(P);
2140
2141     if (emptyQ2(P))
2142         return all_max;
2143
2144     assert(P->NbBid == 0);
2145
2146     if (P->NbEq > 0) {
2147         if (nparam > 0)
2148             CP = compress_parameters(&P, &C, nparam, options->MaxRays);
2149         Q = P;
2150         T = remove_equalities(&P, nparam, options->MaxRays);
2151         if (P != Q && Q != Porig)
2152             Polyhedron_Free(Q);
2153         if (!P) {
2154             if (C != Corig)
2155                 Polyhedron_Free(C);
2156             return all_max;
2157         }
2158     }
2159
2160     if (options->MaxRays & POL_NO_DUAL)
2161         P2PSD_MaxRays = 0;
2162     else
2163         P2PSD_MaxRays = options->MaxRays;
2164
2165     Q = P;
2166     PP = Polyhedron2Param_SimplifiedDomain(&P, C, P2PSD_MaxRays, &CEq, &CT);
2167     if (P != Q && Q != Porig)
2168         Polyhedron_Free(Q);
2169
2170     if (CT) {
2171         if (isIdentity(CT)) {
2172             Matrix_Free(CT);
2173             CT = NULL;
2174         } else
2175             nparam = CT->NbRows - 1;
2176     }
2177     assert(!CT);
2178
2179     unsigned dim = P->Dimension - nparam;
2180
2181     int nd;
2182     for (nd = 0, D=PP->D; D; ++nd, D=D->next);
2183     Polyhedron **fVD = new Polyhedron*[nd];
2184
2185     indicator_constructor ic(P, dim, PP, T);
2186
2187     vector<indicator_term *> all_vertices;
2188     construct_rational_vertices(PP, T, T ? T->NbRows-nparam-1 : dim,
2189                                 nparam, all_vertices);
2190
2191     for (nd = 0, D=PP->D; D; D=next) {
2192         next = D->next;
2193
2194         Polyhedron *rVD = reduce_domain(D->Domain, CT, CEq,
2195                                         fVD, nd, options->MaxRays);
2196         if (!rVD)
2197             continue;
2198
2199         pVD = CT ? DomainImage(rVD,CT,options->MaxRays) : rVD;
2200
2201         EDomain *epVD = new EDomain(pVD);
2202         indicator ind(ic, D, epVD, options);
2203
2204         FORALL_PVertex_in_ParamPolyhedron(V,D,PP) // _i is internal counter
2205             ind.add(all_vertices[_i]);
2206         END_FORALL_PVertex_in_ParamPolyhedron;
2207
2208         ind.remove_initial_rational_vertices();
2209
2210         vector<int> loc;
2211         vector<max_term*> maxima = lexmin(ind, nparam, loc);
2212         if (CP)
2213             for (int j = 0; j < maxima.size(); ++j)
2214                 maxima[j]->substitute(CP, options->MaxRays);
2215         all_max.insert(all_max.end(), maxima.begin(), maxima.end());
2216
2217         ++nd;
2218         if (rVD != pVD)
2219             Domain_Free(pVD);
2220         Domain_Free(rVD);
2221     }
2222     for (int i = 0; i < all_vertices.size(); ++i)
2223         delete all_vertices[i];
2224     if (CP)
2225         Matrix_Free(CP);
2226     if (T)
2227         Matrix_Free(T);
2228     Param_Polyhedron_Free(PP);
2229     if (CEq)
2230         Polyhedron_Free(CEq);
2231     for (--nd; nd >= 0; --nd) {
2232         Domain_Free(fVD[nd]);
2233     }
2234     delete [] fVD;
2235     if (C != Corig)
2236         Polyhedron_Free(C);
2237     if (P != Porig)
2238         Polyhedron_Free(P);
2239
2240     return all_max;
2241 }
2242
2243 static void verify_results(Polyhedron *A, Polyhedron *C,
2244                            vector<max_term*>& maxima, int m, int M,
2245                            int print_all, unsigned MaxRays);
2246
2247 int main(int argc, char **argv)
2248 {
2249     Polyhedron *A, *C;
2250     Matrix *MA;
2251     int bignum;
2252     char **iter_names, **param_names;
2253     int c, ind = 0;
2254     int range = 0;
2255     int verify = 0;
2256     int print_all = 0;
2257     int m = INT_MAX, M = INT_MIN, r;
2258     int print_solution = 1;
2259     struct barvinok_options *options;
2260
2261     options = barvinok_options_new_with_defaults();
2262
2263     while ((c = getopt_long(argc, argv, "TAm:M:r:V", lexmin_options, &ind)) != -1) {
2264         switch (c) {
2265         case NO_EMPTINESS_CHECK:
2266             options->emptiness_check = 0;
2267             break;
2268         case 'T':
2269             verify = 1;
2270             break;
2271         case 'A':
2272             print_all = 1;
2273             break;
2274         case 'm':
2275             m = atoi(optarg);
2276             verify = 1;
2277             break;
2278         case 'M':
2279             M = atoi(optarg);
2280             verify = 1;
2281             break;
2282         case 'r':
2283             M = atoi(optarg);
2284             m = -M;
2285             verify = 1;
2286             break;
2287         case 'V':
2288             printf(barvinok_version());
2289             exit(0);
2290             break;
2291         }
2292     }
2293
2294     MA = Matrix_Read();
2295     C = Constraints2Polyhedron(MA, options->MaxRays);
2296     Matrix_Free(MA);
2297     fscanf(stdin, " %d", &bignum);
2298     assert(bignum == -1);
2299     MA = Matrix_Read();
2300     A = Constraints2Polyhedron(MA, options->MaxRays);
2301     Matrix_Free(MA);
2302
2303     if (A->Dimension >= VBIGDIM)
2304         r = VSRANGE;
2305     else if (A->Dimension >= BIGDIM)
2306         r = SRANGE;
2307     else
2308         r = RANGE;
2309     if (M == INT_MIN)
2310         M = r;
2311     if (m == INT_MAX)
2312         m = -r;
2313
2314     if (verify && m > M) {
2315         fprintf(stderr,"Nothing to do: min > max !\n");
2316         return(0);
2317     }
2318     if (verify)
2319         print_solution = 0;
2320
2321     iter_names = util_generate_names(A->Dimension - C->Dimension, "i");
2322     param_names = util_generate_names(C->Dimension, "p");
2323     if (print_solution) {
2324         Polyhedron_Print(stdout, P_VALUE_FMT, A);
2325         Polyhedron_Print(stdout, P_VALUE_FMT, C);
2326     }
2327     vector<max_term*> maxima = lexmin(A, C, options);
2328     if (print_solution)
2329         for (int i = 0; i < maxima.size(); ++i)
2330             maxima[i]->print(cout, param_names, options);
2331
2332     if (verify)
2333         verify_results(A, C, maxima, m, M, print_all, options->MaxRays);
2334
2335     for (int i = 0; i < maxima.size(); ++i)
2336         delete maxima[i];
2337
2338     util_free_names(A->Dimension - C->Dimension, iter_names);
2339     util_free_names(C->Dimension, param_names);
2340     Polyhedron_Free(A);
2341     Polyhedron_Free(C);
2342
2343     free(options);
2344
2345     return 0;
2346 }
2347
2348 static bool lexmin(int pos, Polyhedron *P, Value *context)
2349 {
2350     Value LB, UB, k;
2351     int lu_flags;
2352     bool found = false;
2353
2354     if (emptyQ(P))
2355         return false;
2356
2357     value_init(LB); value_init(UB); value_init(k);
2358     value_set_si(LB,0);
2359     value_set_si(UB,0);
2360     lu_flags = lower_upper_bounds(pos,P,context,&LB,&UB);
2361     assert(!(lu_flags & LB_INFINITY));
2362
2363     value_set_si(context[pos],0);
2364     if (!lu_flags && value_lt(UB,LB)) {
2365         value_clear(LB); value_clear(UB); value_clear(k);
2366         return false;
2367     }
2368     if (!P->next) {
2369         value_assign(context[pos], LB);
2370         value_clear(LB); value_clear(UB); value_clear(k);
2371         return true;
2372     }
2373     for (value_assign(k,LB); lu_flags || value_le(k,UB); value_increment(k,k)) {
2374         value_assign(context[pos],k);
2375         if ((found = lexmin(pos+1, P->next, context)))
2376             break;
2377     }
2378     if (!found)
2379         value_set_si(context[pos],0);
2380     value_clear(LB); value_clear(UB); value_clear(k);
2381     return found;
2382 }
2383
2384 static void print_list(FILE *out, Value *z, char* brackets, int len)
2385 {
2386     fprintf(out, "%c", brackets[0]);
2387     value_print(out, VALUE_FMT,z[0]);
2388     for (int k = 1; k < len; ++k) {
2389         fprintf(out, ", ");
2390         value_print(out, VALUE_FMT,z[k]);
2391     }
2392     fprintf(out, "%c", brackets[1]);
2393 }
2394
2395 static int check_poly(Polyhedron *S, Polyhedron *CS, vector<max_term*>& maxima,
2396                       int nparam, int pos, Value *z, int print_all, int st,
2397                       unsigned MaxRays)
2398 {
2399     if (pos == nparam) {
2400         int k;
2401         bool found = lexmin(1, S, z);
2402
2403         if (print_all) {
2404             printf("lexmin");
2405             print_list(stdout, z+S->Dimension-nparam+1, "()", nparam);
2406             printf(" = ");
2407             if (found)
2408                 print_list(stdout, z+1, "[]", S->Dimension-nparam);
2409             printf(" ");
2410         }
2411
2412         Vector *min = NULL;
2413         for (int i = 0; i < maxima.size(); ++i)
2414             if ((min = maxima[i]->eval(z+S->Dimension-nparam+1, MaxRays)))
2415                 break;
2416
2417         int ok = !(found ^ !!min);
2418         if (found && min)
2419             for (int i = 0; i < S->Dimension-nparam; ++i)
2420                 if (value_ne(z[1+i], min->p[i])) {
2421                     ok = 0;
2422                     break;
2423                 }
2424         if (!ok) {
2425             printf("\n");
2426             fflush(stdout);
2427             fprintf(stderr, "Error !\n");
2428             fprintf(stderr, "lexmin");
2429             print_list(stderr, z+S->Dimension-nparam+1, "()", nparam);
2430             fprintf(stderr, " should be ");
2431             if (found)
2432                 print_list(stderr, z+1, "[]", S->Dimension-nparam);
2433             fprintf(stderr, " while digging gives ");
2434             if (min)
2435                 print_list(stderr, min->p, "[]", S->Dimension-nparam);
2436             fprintf(stderr, ".\n");
2437             return 0;
2438         } else if (print_all)
2439             printf("OK.\n");
2440         if (min)
2441             Vector_Free(min);
2442
2443         for (k = 1; k <= S->Dimension-nparam; ++k)
2444             value_set_si(z[k], 0);
2445     } else {
2446         Value tmp;
2447         Value LB, UB;
2448         value_init(tmp);
2449         value_init(LB);
2450         value_init(UB);
2451         int ok =
2452             !(lower_upper_bounds(1+pos, CS, &z[S->Dimension-nparam], &LB, &UB));
2453         for (value_assign(tmp,LB); value_le(tmp,UB); value_increment(tmp,tmp)) {
2454             if (!print_all) {
2455                 int k = VALUE_TO_INT(tmp);
2456                 if (!pos && !(k%st)) {
2457                     printf("o");
2458                     fflush(stdout);
2459                 }
2460             }
2461             value_assign(z[pos+S->Dimension-nparam+1],tmp);
2462             if (!check_poly(S, CS->next, maxima, nparam, pos+1, z, print_all, st,
2463                             MaxRays)) {
2464                 value_clear(tmp);
2465                 value_clear(LB);
2466                 value_clear(UB);
2467                 return 0;
2468             }
2469         }
2470         value_set_si(z[pos+S->Dimension-nparam+1],0);
2471         value_clear(tmp);
2472         value_clear(LB);
2473         value_clear(UB);
2474     }
2475     return 1;
2476 }
2477
2478 void verify_results(Polyhedron *A, Polyhedron *C, vector<max_term*>& maxima,
2479                     int m, int M, int print_all, unsigned MaxRays)
2480 {
2481     Polyhedron *CC, *CC2, *CS, *S;
2482     unsigned nparam = C->Dimension;
2483     Value *p;
2484     int i;
2485     int st;
2486
2487     CC = Polyhedron_Project(A, nparam);
2488     CC2 = DomainIntersection(C, CC, MaxRays);
2489     Domain_Free(CC);
2490     CC = CC2;
2491
2492     /* Intersect context with range */
2493     if (nparam > 0) {
2494         Matrix *MM;
2495         Polyhedron *U;
2496
2497         MM = Matrix_Alloc(2*C->Dimension, C->Dimension+2);
2498         for (int i = 0; i < C->Dimension; ++i) {
2499             value_set_si(MM->p[2*i][0], 1);
2500             value_set_si(MM->p[2*i][1+i], 1);
2501             value_set_si(MM->p[2*i][1+C->Dimension], -m);
2502             value_set_si(MM->p[2*i+1][0], 1);
2503             value_set_si(MM->p[2*i+1][1+i], -1);
2504             value_set_si(MM->p[2*i+1][1+C->Dimension], M);
2505         }
2506         CC2 = AddConstraints(MM->p[0], 2*CC->Dimension, CC, MaxRays);
2507         U = Universe_Polyhedron(0);
2508         CS = Polyhedron_Scan(CC2, U, MaxRays & POL_NO_DUAL ? 0 : MaxRays);
2509         Polyhedron_Free(U);
2510         Polyhedron_Free(CC2);
2511         Matrix_Free(MM);
2512     } else
2513         CS = NULL;
2514
2515     p = ALLOCN(Value, A->Dimension+2);
2516     for (i=0; i <= A->Dimension; i++) {
2517         value_init(p[i]);
2518         value_set_si(p[i],0);
2519     }
2520     value_init(p[i]);
2521     value_set_si(p[i], 1);
2522
2523     S = Polyhedron_Scan(A, C, MaxRays & POL_NO_DUAL ? 0 : MaxRays);
2524
2525     if (!print_all && C->Dimension > 0) {
2526         if (M-m > 80)
2527             st = 1+(M-m)/80;
2528         else
2529             st = 1;
2530         for (int i = m; i <= M; i += st)
2531             printf(".");
2532         printf( "\r" );
2533         fflush(stdout);
2534     }
2535
2536     if (S) {
2537         if (!(CS && emptyQ2(CS)))
2538             check_poly(S, CS, maxima, nparam, 0, p, print_all, st, MaxRays);
2539         Domain_Free(S);
2540     }
2541
2542     if (!print_all)
2543         printf("\n");
2544
2545     for (i=0; i <= (A->Dimension+1); i++)
2546         value_clear(p[i]);
2547     free(p);
2548     if (CS)
2549         Domain_Free(CS);
2550     Polyhedron_Free(CC);
2551 }