lexmin.cc

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