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 static Matrix *remove_equalities(Polyhedron **P, unsigned nparam, unsigned MaxRays);
1466
1467 Vector *Polyhedron_not_empty(Polyhedron *P, unsigned MaxRays)
1468 {
1469     Polyhedron *Porig = P;
1470     Vector *sample = NULL;
1471
1472     POL_ENSURE_VERTICES(P);
1473     if (emptyQ2(P))
1474         return NULL;
1475
1476     for (int i = 0; i < P->NbRays; ++i)
1477         if (value_one_p(P->Ray[i][1+P->Dimension])) {
1478             sample = Vector_Alloc(P->Dimension + 1);
1479             Vector_Copy(P->Ray[i]+1, sample->p, P->Dimension+1);
1480             return sample;
1481         }
1482
1483     Matrix *T = remove_equalities(&P, 0, MaxRays);
1484     if (P)
1485         sample = Polyhedron_Sample(P, MaxRays);
1486     if (sample) {
1487         if (T) {
1488             Vector *P_sample = Vector_Alloc(Porig->Dimension + 1);
1489             Matrix_Vector_Product(T, sample->p, P_sample->p);
1490             Vector_Free(sample);
1491             sample = P_sample;
1492         }
1493     }
1494     if (T) {
1495         Polyhedron_Free(P);
1496         Matrix_Free(T);
1497     }
1498
1499     return sample;
1500 }
1501
1502 struct split {
1503     evalue *constraint;
1504     enum sign { le, ge, lge } sign;
1505
1506     split (evalue *c, enum sign s) : constraint(c), sign(s) {}
1507 };
1508
1509 static void split_on(const split& sp, EDomain *D,
1510                      EDomain **Dlt, EDomain **Deq, EDomain **Dgt,
1511                      unsigned MaxRays)
1512 {
1513     Matrix *M, *M2;
1514     EDomain *EDlt = NULL, *EDeq = NULL, *EDgt = NULL;
1515     Polyhedron *D2;
1516     Value mone;
1517     value_init(mone);
1518     value_set_si(mone, -1);
1519     *Dlt = NULL;
1520     *Deq = NULL;
1521     *Dgt = NULL;
1522     vector<evalue *> new_floors;
1523     M = add_ge_constraint(D, sp.constraint, new_floors);
1524     if (sp.sign == split::lge || sp.sign == split::ge) {
1525         M2 = Matrix_Copy(M);
1526         value_decrement(M2->p[M2->NbRows-1][M2->NbColumns-1],
1527                         M2->p[M2->NbRows-1][M2->NbColumns-1]);
1528         D2 = Constraints2Polyhedron(M2, MaxRays);
1529         EDgt = new EDomain(D2, D, new_floors);
1530         Polyhedron_Free(D2);
1531         Matrix_Free(M2);
1532     }
1533     if (sp.sign == split::lge || sp.sign == split::le) {
1534         M2 = Matrix_Copy(M);
1535         Vector_Scale(M2->p[M2->NbRows-1]+1, M2->p[M2->NbRows-1]+1,
1536                      mone, M2->NbColumns-1);
1537         value_decrement(M2->p[M2->NbRows-1][M2->NbColumns-1],
1538                         M2->p[M2->NbRows-1][M2->NbColumns-1]);
1539         D2 = Constraints2Polyhedron(M2, MaxRays);
1540         EDlt = new EDomain(D2, D, new_floors);
1541         Polyhedron_Free(D2);
1542         Matrix_Free(M2);
1543     }
1544
1545     assert(sp.sign == split::lge || sp.sign == split::ge || sp.sign == split::le);
1546     value_set_si(M->p[M->NbRows-1][0], 0);
1547     D2 = Constraints2Polyhedron(M, MaxRays);
1548     EDeq = new EDomain(D2, D, new_floors);
1549     Polyhedron_Free(D2);
1550     Matrix_Free(M);
1551
1552     Vector *sample = D->sample;
1553     if (sample && new_floors.size() > 0) {
1554         assert(sample->Size == D->D->Dimension+1);
1555         sample = Vector_Alloc(D->D->Dimension+new_floors.size()+1);
1556         Vector_Copy(D->sample->p, sample->p, D->D->Dimension);
1557         value_set_si(sample->p[D->D->Dimension+new_floors.size()], 1);
1558         for (int i = 0; i < new_floors.size(); ++i)
1559             compute_evalue(new_floors[i], sample->p, sample->p+D->D->Dimension+i);
1560     }
1561
1562     for (int i = 0; i < new_floors.size(); ++i) {
1563         free_evalue_refs(new_floors[i]);
1564         delete new_floors[i];
1565     }
1566
1567     if (EDeq) {
1568         if (sample && in_domain(EDeq->D, sample->p, sample->Size-1, MaxRays)) {
1569             EDeq->sample = Vector_Alloc(sample->Size);
1570             Vector_Copy(sample->p, EDeq->sample->p, sample->Size);
1571         } else if (!(EDeq->sample = Polyhedron_not_empty(EDeq->D, MaxRays))) {
1572             delete EDeq;
1573             EDeq = NULL;
1574         }
1575     }
1576     if (EDgt) {
1577         if (sample && in_domain(EDgt->D, sample->p, sample->Size-1, MaxRays)) {
1578             EDgt->sample = Vector_Alloc(sample->Size);
1579             Vector_Copy(sample->p, EDgt->sample->p, sample->Size);
1580         } else if (!(EDgt->sample = Polyhedron_not_empty(EDgt->D, MaxRays))) {
1581             delete EDgt;
1582             EDgt = NULL;
1583         }
1584     }
1585     if (EDlt) {
1586         if (sample && in_domain(EDlt->D, sample->p, sample->Size-1, MaxRays)) {
1587             EDlt->sample = Vector_Alloc(sample->Size);
1588             Vector_Copy(sample->p, EDlt->sample->p, sample->Size);
1589         } else if (!(EDlt->sample = Polyhedron_not_empty(EDlt->D, MaxRays))) {
1590             delete EDlt;
1591             EDlt = NULL;
1592         }
1593     }
1594     *Dlt = EDlt;
1595     *Deq = EDeq;
1596     *Dgt = EDgt;
1597     value_clear(mone);
1598     if (sample != D->sample)
1599         Vector_Free(sample);
1600 }
1601
1602 ostream & operator<< (ostream & os, const vector<int> & v)
1603 {
1604     os << "[";
1605     for (int i = 0; i < v.size(); ++i) {
1606         if (i)
1607             os << ", ";
1608         os << v[i];
1609     }
1610     os << "]";
1611     return os;
1612 }
1613
1614 /*
1615  * Project on first dim dimensions
1616  */
1617 Polyhedron* Polyhedron_Project_Initial(Polyhedron *P, int dim)
1618 {
1619     int i;
1620     Matrix *T;
1621     Polyhedron *I;
1622
1623     if (P->Dimension == dim)
1624         return Polyhedron_Copy(P);
1625
1626     T = Matrix_Alloc(dim+1, P->Dimension+1);
1627     for (i = 0; i < dim; ++i)
1628         value_set_si(T->p[i][i], 1);
1629     value_set_si(T->p[dim][P->Dimension], 1);
1630     I = Polyhedron_Image(P, T, P->NbConstraints);
1631     Matrix_Free(T);
1632     return I;
1633 }
1634
1635 static vector<max_term*> lexmin(indicator& ind, EDomain *D, unsigned nparam,
1636                                 unsigned MaxRays, vector<int> loc)
1637 {
1638     vector<max_term*> maxima;
1639     int len = 1 + D->D->Dimension + 1;
1640     Value lcm, a, b;
1641     evalue mone;
1642     EDomain *Dorig = D;
1643
1644     value_init(mone.d);
1645     evalue_set_si(&mone, -1, 1);
1646     value_init(lcm);
1647     value_init(a);
1648     value_init(b);
1649     Vector *c = Vector_Alloc(len);
1650     Matrix *T = Matrix_Alloc(2, len-1);
1651     for (int i = 0; i < ind.term.size(); ++i) {
1652         bool restart = false;   /* true if we have modified ind from i up */
1653         bool stop = false;      /* true if i can never be smallest */
1654         int peel = -1;          /* term to peel against */
1655         vector<split> splits;
1656         if (ind.term[i]->sign < 0)
1657             continue;
1658         int dim = ind.term[i]->den.NumCols();
1659         int j;
1660         for (j = 0; j < ind.term.size(); ++j) {
1661             if (i == j)
1662                 continue;
1663             int k;
1664             for (k = 0; k < dim; ++k) {
1665                 /* compute ind.term->[i]->vertex[k] - ind.term->[j]->vertex[k] */
1666                 evalue *diff = new evalue;
1667                 value_init(diff->d);
1668                 evalue_copy(diff, ind.term[j]->vertex[k]);
1669                 emul(&mone, diff);
1670                 eadd(ind.term[i]->vertex[k], diff);
1671                 reduce_evalue(diff);
1672                 int fract = evalue2constraint(D, diff, c->p, len);
1673                 Vector_Copy(c->p+1, T->p[0], len-1);
1674                 value_assign(T->p[1][len-2], c->p[0]);
1675
1676                 int min, max;
1677                 interval_minmax(D->D, T, &min, &max, MaxRays);
1678                 if (max < 0) {
1679                     free_evalue_refs(diff);
1680                     delete diff;
1681                     break;
1682                 }
1683                 if (fract) {
1684                     emul(&mone, diff);
1685                     evalue2constraint(D, diff, c->p, len);
1686                     emul(&mone, diff);
1687                     Vector_Copy(c->p+1, T->p[0], len-1);
1688                     value_assign(T->p[1][len-2], c->p[0]);
1689
1690                     int negmin, negmax;
1691                     interval_minmax(D->D, T, &negmin, &negmax, MaxRays);
1692                     min = -negmax;
1693                 }
1694                 if (min > 0) {
1695                     free_evalue_refs(diff);
1696                     delete diff;
1697                     stop = true;
1698                     break;
1699                 }
1700                 if (max == 0 && min == 0) {
1701                     if (!EVALUE_IS_ZERO(*diff)) {
1702                         ind.substitute(diff);
1703                         ind.simplify();
1704                         restart = true;
1705                     }
1706                     free_evalue_refs(diff);
1707                     delete diff;
1708                     if (restart)
1709                         break;
1710                     continue;
1711                 }
1712                 if (min < 0 && max == 0)
1713                     splits.push_back(split(diff, split::le));
1714                 else if (max > 0 && min == 0)
1715                     splits.push_back(split(diff, split::ge));
1716                 else
1717                     splits.push_back(split(diff, split::lge));
1718                 break;
1719             }
1720             if (k == dim && ind.term[j]->sign < 0)
1721                 peel = j;
1722             if (stop || restart)
1723                 break;
1724         }
1725         if (restart) {
1726             /* The ith entry may have been removed, so we have to consider
1727              * it again.
1728              */
1729             --i;
1730             for (j = 0; j < splits.size(); ++j) {
1731                 free_evalue_refs(splits[j].constraint);
1732                 delete splits[j].constraint;
1733             }
1734             continue;
1735         }
1736         if (stop) {
1737             for (j = 0; j < splits.size(); ++j) {
1738                 free_evalue_refs(splits[j].constraint);
1739                 delete splits[j].constraint;
1740             }
1741             continue;
1742         }
1743         if (peel != -1) {
1744             // ind.peel(i, peel);
1745             ind.combine(i, peel);
1746             ind.simplify();
1747             i = -1;                 /* start over */
1748             for (j = 0; j < splits.size(); ++j) {
1749                 free_evalue_refs(splits[j].constraint);
1750                 delete splits[j].constraint;
1751             }
1752             continue;
1753         }
1754         if (splits.size() != 0) {
1755             for (j = 0; j < splits.size(); ++j)
1756                 if (splits[j].sign == split::le)
1757                     break;
1758             if (j == splits.size())
1759                 j = 0;
1760                 EDomain *Dlt, *Deq, *Dgt;
1761                 split_on(splits[j], D, &Dlt, &Deq, &Dgt, MaxRays);
1762                 assert(Dlt || Deq || Dgt);
1763                 if (Deq) {
1764                     loc.push_back(0);
1765                     indicator indeq(ind);
1766                     indeq.substitute(splits[j].constraint);
1767                     Polyhedron *P = Polyhedron_Project_Initial(Deq->D, nparam);
1768                     P = DomainConstraintSimplify(P, MaxRays);
1769                     indeq.reduce_in_domain(P);
1770                     Polyhedron_Free(P);
1771                     indeq.simplify();
1772                     vector<max_term*> maxeq = lexmin(indeq, Deq, nparam,
1773                                                      MaxRays, loc);
1774                     maxima.insert(maxima.end(), maxeq.begin(), maxeq.end());
1775                     loc.pop_back();
1776                     delete Deq;
1777                 }
1778                 if (Dgt) {
1779                     loc.push_back(1);
1780                     indicator indgt(ind);
1781                     Polyhedron *P = Polyhedron_Project_Initial(Dgt->D, nparam);
1782                     P = DomainConstraintSimplify(P, MaxRays);
1783                     indgt.reduce_in_domain(P);
1784                     Polyhedron_Free(P);
1785                     indgt.simplify();
1786                     vector<max_term*> maxeq = lexmin(indgt, Dgt, nparam,
1787                                                      MaxRays, loc);
1788                     maxima.insert(maxima.end(), maxeq.begin(), maxeq.end());
1789                     loc.pop_back();
1790                     delete Dgt;
1791                 }
1792                 if (Dlt) {
1793                     loc.push_back(-1);
1794                     Polyhedron *P = Polyhedron_Project_Initial(Dlt->D, nparam);
1795                     P = DomainConstraintSimplify(P, MaxRays);
1796                     ind.reduce_in_domain(P);
1797                     Polyhedron_Free(P);
1798                     ind.simplify();
1799                     if (D != Dorig)
1800                         delete D;
1801                     D = Dlt;
1802                     if (splits.size() > 1) {
1803                         vector<max_term*> maxeq = lexmin(ind, Dlt, nparam,
1804                                                          MaxRays, loc);
1805                         maxima.insert(maxima.end(), maxeq.begin(), maxeq.end());
1806                         for (j = 0; j < splits.size(); ++j) {
1807                             free_evalue_refs(splits[j].constraint);
1808                             delete splits[j].constraint;
1809                         }
1810                         break;
1811                     }
1812                 }
1813             /* the vertex turned out not to be minimal */
1814             for (j = 0; j < splits.size(); ++j) {
1815                 free_evalue_refs(splits[j].constraint);
1816                 delete splits[j].constraint;
1817             }
1818             if (!Dlt)
1819                 break;
1820         }
1821             max_term *maximum = new max_term;
1822             maxima.push_back(maximum);
1823             maximum->dim = nparam;
1824             maximum->domain = Polyhedron_Copy(D->D);
1825             for (int j = 0; j < dim; ++j) {
1826                 evalue *E = new evalue;
1827                 value_init(E->d);
1828                 evalue_copy(E, ind.term[i]->vertex[j]);
1829                 if (evalue_frac2floor_in_domain(E, D->D))
1830                     reduce_evalue(E);
1831                 maximum->max.push_back(E);
1832             }
1833             break;
1834     }
1835     Matrix_Free(T);
1836     Vector_Free(c);
1837     value_clear(lcm);
1838     value_clear(a);
1839     value_clear(b);
1840     free_evalue_refs(&mone);
1841     if (D != Dorig)
1842         delete D;
1843     return maxima;
1844 }
1845
1846 static bool isTranslation(Matrix *M)
1847 {
1848     unsigned i, j;
1849     if (M->NbRows != M->NbColumns)
1850         return False;
1851
1852     for (i = 0;i < M->NbRows; i ++)
1853         for (j = 0; j < M->NbColumns-1; j ++)
1854             if (i == j) {
1855                 if(value_notone_p(M->p[i][j]))
1856                     return False;
1857             } else {
1858                 if(value_notzero_p(M->p[i][j]))
1859                     return False;
1860             }
1861     return value_one_p(M->p[M->NbRows-1][M->NbColumns-1]);
1862 }
1863
1864 static Matrix *compress_parameters(Polyhedron **P, Polyhedron **C,
1865                                    unsigned nparam, unsigned MaxRays)
1866 {
1867     Matrix *M, *T, *CP;
1868
1869     /* compress_parms doesn't like equalities that only involve parameters */
1870     for (int i = 0; i < (*P)->NbEq; ++i)
1871         assert(First_Non_Zero((*P)->Constraint[i]+1, (*P)->Dimension-nparam) != -1);
1872
1873     M = Matrix_Alloc((*P)->NbEq, (*P)->Dimension+2);
1874     Vector_Copy((*P)->Constraint[0], M->p[0], (*P)->NbEq * ((*P)->Dimension+2));
1875     CP = compress_parms(M, nparam);
1876     Matrix_Free(M);
1877
1878     if (isTranslation(CP)) {
1879         Matrix_Free(CP);
1880         return NULL;
1881     }
1882
1883     T = align_matrix(CP, (*P)->Dimension+1);
1884     *P = Polyhedron_Preimage(*P, T, MaxRays);
1885     Matrix_Free(T);
1886
1887     *C = Polyhedron_Preimage(*C, CP, MaxRays);
1888
1889     return CP;
1890 }
1891
1892 static Matrix *remove_equalities(Polyhedron **P, unsigned nparam, unsigned MaxRays)
1893 {
1894     /* Matrix "view" of equalities */
1895     Matrix M;
1896     M.NbRows = (*P)->NbEq;
1897     M.NbColumns = (*P)->Dimension+2;
1898     M.p_Init = (*P)->p_Init;
1899     M.p = (*P)->Constraint;
1900
1901     Matrix *T = compress_variables(&M, nparam);
1902
1903     if (!T) {
1904         *P = NULL;
1905         return NULL;
1906     }
1907     if (isIdentity(T)) {
1908         Matrix_Free(T);
1909         T = NULL;
1910     } else
1911         *P = Polyhedron_Preimage(*P, T, MaxRays);
1912
1913     return T;
1914 }
1915
1916 static vector<max_term*> lexmin(Polyhedron *P, Polyhedron *C, unsigned MaxRays)
1917 {
1918     unsigned nparam = C->Dimension;
1919     Param_Polyhedron *PP = NULL;
1920     Polyhedron *CEq = NULL, *pVD;
1921     Matrix *CT = NULL;
1922     Matrix *T = NULL, *CP = NULL;
1923     Param_Domain *D, *next;
1924     Param_Vertices *V;
1925     Polyhedron *Porig = P;
1926     Polyhedron *Corig = C;
1927     int i;
1928     vector<max_term*> all_max;
1929     Polyhedron *Q;
1930
1931     if (emptyQ2(P))
1932         return all_max;
1933
1934     POL_ENSURE_VERTICES(P);
1935
1936     if (emptyQ2(P))
1937         return all_max;
1938
1939     assert(P->NbBid == 0);
1940
1941     if (P->NbEq > 0) {
1942         if (nparam > 0)
1943             CP = compress_parameters(&P, &C, nparam, MaxRays);
1944         Q = P;
1945         T = remove_equalities(&P, nparam, MaxRays);
1946         if (P != Q && Q != Porig)
1947             Polyhedron_Free(Q);
1948         if (!P) {
1949             if (C != Corig)
1950                 Polyhedron_Free(C);
1951             return all_max;
1952         }
1953     }
1954
1955     Q = P;
1956     PP = Polyhedron2Param_SimplifiedDomain(&P,C,
1957                                            (MaxRays & POL_NO_DUAL) ? 0 : MaxRays,
1958                                            &CEq,&CT);
1959     if (P != Q && Q != Porig)
1960         Polyhedron_Free(Q);
1961
1962     if (CT) {
1963         if (isIdentity(CT)) {
1964             Matrix_Free(CT);
1965             CT = NULL;
1966         } else
1967             nparam = CT->NbRows - 1;
1968     }
1969
1970     unsigned dim = P->Dimension - nparam;
1971
1972     int nd;
1973     for (nd = 0, D=PP->D; D; ++nd, D=D->next);
1974     Polyhedron **fVD = new Polyhedron*[nd];
1975
1976     indicator_constructor ic(P, dim, PP->nbV);
1977
1978     for (i = 0, V = PP->V; V; V = V->next, i++) {
1979         ic.decompose_at_vertex(V, i, MaxRays);
1980     }
1981     if (T) {
1982         ic.substitute(T);
1983         ic.normalize();
1984     }
1985
1986     for (nd = 0, D=PP->D; D; D=next) {
1987         next = D->next;
1988
1989         Polyhedron *rVD = reduce_domain(D->Domain, CT, CEq,
1990                                         fVD, nd, MaxRays);
1991         if (!rVD)
1992             continue;
1993
1994         pVD = CT ? DomainImage(rVD,CT,MaxRays) : rVD;
1995
1996         indicator ind;
1997
1998         FORALL_PVertex_in_ParamPolyhedron(V,D,PP) // _i is internal counter
1999             for (int j = 0; j < ic.terms[_i].size(); ++j) {
2000                 indicator_term *term = new indicator_term(*ic.terms[_i][j]);
2001                 term->reduce_in_domain(pVD);
2002                 ind.term.push_back(term);
2003             }
2004         END_FORALL_PVertex_in_ParamPolyhedron;
2005
2006         ind.simplify();
2007
2008         EDomain epVD(pVD);
2009         vector<int> loc;
2010         vector<max_term*> maxima = lexmin(ind, &epVD, nparam, MaxRays, loc);
2011         if (CP)
2012             for (int j = 0; j < maxima.size(); ++j)
2013                 maxima[j]->substitute(CP, MaxRays);
2014         all_max.insert(all_max.end(), maxima.begin(), maxima.end());
2015
2016         ++nd;
2017         if (rVD != pVD)
2018             Domain_Free(pVD);
2019         Domain_Free(rVD);
2020     }
2021     if (CP)
2022         Matrix_Free(CP);
2023     if (T)
2024         Matrix_Free(T);
2025     Param_Polyhedron_Free(PP);
2026     if (CEq)
2027         Polyhedron_Free(CEq);
2028     for (--nd; nd >= 0; --nd) {
2029         Domain_Free(fVD[nd]);
2030     }
2031     delete [] fVD;
2032     if (C != Corig)
2033         Polyhedron_Free(C);
2034     if (P != Porig)
2035         Polyhedron_Free(P);
2036
2037     return all_max;
2038 }
2039
2040 static void verify_results(Polyhedron *A, Polyhedron *C,
2041                            vector<max_term*>& maxima, int m, int M,
2042                            int print_all, unsigned MaxRays);
2043
2044 int main(int argc, char **argv)
2045 {
2046     Polyhedron *A, *C;
2047     Matrix *MA;
2048     int bignum;
2049     char **iter_names, **param_names;
2050     int c, ind = 0;
2051     int range = 0;
2052     int verify = 0;
2053     int print_all = 0;
2054     int m = INT_MAX, M = INT_MIN, r;
2055     int print_solution = 1;
2056
2057     while ((c = getopt_long(argc, argv, "TAm:M:r:V", options, &ind)) != -1) {
2058         switch (c) {
2059         case 'T':
2060             verify = 1;
2061             break;
2062         case 'A':
2063             print_all = 1;
2064             break;
2065         case 'm':
2066             m = atoi(optarg);
2067             verify = 1;
2068             break;
2069         case 'M':
2070             M = atoi(optarg);
2071             verify = 1;
2072             break;
2073         case 'r':
2074             M = atoi(optarg);
2075             m = -M;
2076             verify = 1;
2077             break;
2078         case 'V':
2079             printf(barvinok_version());
2080             exit(0);
2081             break;
2082         }
2083     }
2084
2085     MA = Matrix_Read();
2086     C = Constraints2Polyhedron(MA, MAXRAYS);
2087     Matrix_Free(MA);
2088     fscanf(stdin, " %d", &bignum);
2089     assert(bignum == -1);
2090     MA = Matrix_Read();
2091     A = Constraints2Polyhedron(MA, MAXRAYS);
2092     Matrix_Free(MA);
2093
2094     if (A->Dimension >= VBIGDIM)
2095         r = VSRANGE;
2096     else if (A->Dimension >= BIGDIM)
2097         r = SRANGE;
2098     else
2099         r = RANGE;
2100     if (M == INT_MIN)
2101         M = r;
2102     if (m == INT_MAX)
2103         m = -r;
2104
2105     if (verify && m > M) {
2106         fprintf(stderr,"Nothing to do: min > max !\n");
2107         return(0);
2108     }
2109     if (verify)
2110         print_solution = 0;
2111
2112     iter_names = util_generate_names(A->Dimension - C->Dimension, "i");
2113     param_names = util_generate_names(C->Dimension, "p");
2114     if (print_solution) {
2115         Polyhedron_Print(stdout, P_VALUE_FMT, A);
2116         Polyhedron_Print(stdout, P_VALUE_FMT, C);
2117     }
2118     vector<max_term*> maxima = lexmin(A, C, MAXRAYS);
2119     if (print_solution)
2120         for (int i = 0; i < maxima.size(); ++i)
2121             maxima[i]->print(cout, param_names);
2122
2123     if (verify)
2124         verify_results(A, C, maxima, m, M, print_all, MAXRAYS);
2125
2126     for (int i = 0; i < maxima.size(); ++i)
2127         delete maxima[i];
2128
2129     util_free_names(A->Dimension - C->Dimension, iter_names);
2130     util_free_names(C->Dimension, param_names);
2131     Polyhedron_Free(A);
2132     Polyhedron_Free(C);
2133
2134     return 0;
2135 }
2136
2137 static bool lexmin(int pos, Polyhedron *P, Value *context)
2138 {
2139     Value LB, UB, k;
2140     int lu_flags;
2141     bool found = false;
2142
2143     if (emptyQ(P))
2144         return false;
2145
2146     value_init(LB); value_init(UB); value_init(k);
2147     value_set_si(LB,0);
2148     value_set_si(UB,0);
2149     lu_flags = lower_upper_bounds(pos,P,context,&LB,&UB);
2150     assert(!(lu_flags & LB_INFINITY));
2151
2152     value_set_si(context[pos],0);
2153     if (!lu_flags && value_lt(UB,LB)) {
2154         value_clear(LB); value_clear(UB); value_clear(k);
2155         return false;
2156     }
2157     if (!P->next) {
2158         value_assign(context[pos], LB);
2159         value_clear(LB); value_clear(UB); value_clear(k);
2160         return true;
2161     }
2162     for (value_assign(k,LB); lu_flags || value_le(k,UB); value_increment(k,k)) {
2163         value_assign(context[pos],k);
2164         if ((found = lexmin(pos+1, P->next, context)))
2165             break;
2166     }
2167     if (!found)
2168         value_set_si(context[pos],0);
2169     value_clear(LB); value_clear(UB); value_clear(k);
2170     return found;
2171 }
2172
2173 static void print_list(FILE *out, Value *z, char* brackets, int len)
2174 {
2175     fprintf(out, "%c", brackets[0]);
2176     value_print(out, VALUE_FMT,z[0]);
2177     for (int k = 1; k < len; ++k) {
2178         fprintf(out, ", ");
2179         value_print(out, VALUE_FMT,z[k]);
2180     }
2181     fprintf(out, "%c", brackets[1]);
2182 }
2183
2184 static int check_poly(Polyhedron *S, Polyhedron *CS, vector<max_term*>& maxima,
2185                       int nparam, int pos, Value *z, int print_all, int st,
2186                       unsigned MaxRays)
2187 {
2188     if (pos == nparam) {
2189         int k;
2190         bool found = lexmin(1, S, z);
2191
2192         if (print_all) {
2193             printf("lexmin");
2194             print_list(stdout, z+S->Dimension-nparam+1, "()", nparam);
2195             printf(" = ");
2196             if (found)
2197                 print_list(stdout, z+1, "[]", S->Dimension-nparam);
2198             printf(" ");
2199         }
2200
2201         Vector *min = NULL;
2202         for (int i = 0; i < maxima.size(); ++i)
2203             if ((min = maxima[i]->eval(z+S->Dimension-nparam+1, MaxRays)))
2204                 break;
2205
2206         int ok = !(found ^ !!min);
2207         if (found && min)
2208             for (int i = 0; i < S->Dimension-nparam; ++i)
2209                 if (value_ne(z[1+i], min->p[i])) {
2210                     ok = 0;
2211                     break;
2212                 }
2213         if (!ok) {
2214             printf("\n");
2215             fflush(stdout);
2216             fprintf(stderr, "Error !\n");
2217             fprintf(stderr, "lexmin");
2218             print_list(stderr, z+S->Dimension-nparam+1, "()", nparam);
2219             fprintf(stderr, " should be ");
2220             if (found)
2221                 print_list(stderr, z+1, "[]", S->Dimension-nparam);
2222             fprintf(stderr, " while digging gives ");
2223             if (min)
2224                 print_list(stderr, min->p, "[]", S->Dimension-nparam);
2225             fprintf(stderr, ".\n");
2226             return 0;
2227         } else if (print_all)
2228             printf("OK.\n");
2229         if (min)
2230             Vector_Free(min);
2231
2232         for (k = 1; k <= S->Dimension-nparam; ++k)
2233             value_set_si(z[k], 0);
2234     } else {
2235         Value tmp;
2236         Value LB, UB;
2237         value_init(tmp);
2238         value_init(LB);
2239         value_init(UB);
2240         int ok =
2241             !(lower_upper_bounds(1+pos, CS, &z[S->Dimension-nparam], &LB, &UB));
2242         for (value_assign(tmp,LB); value_le(tmp,UB); value_increment(tmp,tmp)) {
2243             if (!print_all) {
2244                 int k = VALUE_TO_INT(tmp);
2245                 if (!pos && !(k%st)) {
2246                     printf("o");
2247                     fflush(stdout);
2248                 }
2249             }
2250             value_assign(z[pos+S->Dimension-nparam+1],tmp);
2251             if (!check_poly(S, CS->next, maxima, nparam, pos+1, z, print_all, st,
2252                             MaxRays)) {
2253                 value_clear(tmp);
2254                 value_clear(LB);
2255                 value_clear(UB);
2256                 return 0;
2257             }
2258         }
2259         value_set_si(z[pos+S->Dimension-nparam+1],0);
2260         value_clear(tmp);
2261         value_clear(LB);
2262         value_clear(UB);
2263     }
2264     return 1;
2265 }
2266
2267 void verify_results(Polyhedron *A, Polyhedron *C, vector<max_term*>& maxima,
2268                     int m, int M, int print_all, unsigned MaxRays)
2269 {
2270     Polyhedron *CC, *CC2, *CS, *S;
2271     unsigned nparam = C->Dimension;
2272     Value *p;
2273     int i;
2274     int st;
2275
2276     CC = Polyhedron_Project(A, nparam);
2277     CC2 = DomainIntersection(C, CC, MAXRAYS);
2278     Domain_Free(CC);
2279     CC = CC2;
2280
2281     /* Intersect context with range */
2282     if (nparam > 0) {
2283         Matrix *MM;
2284         Polyhedron *U;
2285
2286         MM = Matrix_Alloc(2*C->Dimension, C->Dimension+2);
2287         for (int i = 0; i < C->Dimension; ++i) {
2288             value_set_si(MM->p[2*i][0], 1);
2289             value_set_si(MM->p[2*i][1+i], 1);
2290             value_set_si(MM->p[2*i][1+C->Dimension], -m);
2291             value_set_si(MM->p[2*i+1][0], 1);
2292             value_set_si(MM->p[2*i+1][1+i], -1);
2293             value_set_si(MM->p[2*i+1][1+C->Dimension], M);
2294         }
2295         CC2 = AddConstraints(MM->p[0], 2*CC->Dimension, CC, MAXRAYS);
2296         U = Universe_Polyhedron(0);
2297         CS = Polyhedron_Scan(CC2, U, MAXRAYS & POL_NO_DUAL ? 0 : MAXRAYS);
2298         Polyhedron_Free(U);
2299         Polyhedron_Free(CC2);
2300         Matrix_Free(MM);
2301     } else
2302         CS = NULL;
2303
2304     p = ALLOCN(Value, A->Dimension+2);
2305     for (i=0; i <= A->Dimension; i++) {
2306         value_init(p[i]);
2307         value_set_si(p[i],0);
2308     }
2309     value_init(p[i]);
2310     value_set_si(p[i], 1);
2311
2312     S = Polyhedron_Scan(A, C, MAXRAYS & POL_NO_DUAL ? 0 : MAXRAYS);
2313
2314     if (!print_all && C->Dimension > 0) {
2315         if (M-m > 80)
2316             st = 1+(M-m)/80;
2317         else
2318             st = 1;
2319         for (int i = m; i <= M; i += st)
2320             printf(".");
2321         printf( "\r" );
2322         fflush(stdout);
2323     }
2324
2325     if (S) {
2326         if (!(CS && emptyQ2(CS)))
2327             check_poly(S, CS, maxima, nparam, 0, p, print_all, st, MaxRays);
2328         Domain_Free(S);
2329     }
2330
2331     if (!print_all)
2332         printf("\n");
2333
2334     for (i=0; i <= (A->Dimension+1); i++)
2335         value_clear(p[i]);
2336     free(p);
2337     if (CS)
2338         Domain_Free(CS);
2339     Polyhedron_Free(CC);
2340 }