maximize.cc

   1 #include <iostream>
   2 #include <bernstein/bernstein.h>
   3 #include <barvinok/evalue.h>
   4 #include <barvinok/options.h>
   5 #include <barvinok/util.h>
   6 #include <barvinok/bernstein.h>
   7 #include "argp.h"
   8 #include "progname.h"
   9 #include "evalue_convert.h"
  10 #include "verify.h"
  11
  12 using std::cout;
  13 using std::cerr;
  14 using std::endl;
  15 using namespace GiNaC;
  16 using namespace bernstein;
  17 using namespace barvinok;
  18
  19 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
  20
  21 #define OPT_VARS            (BV_OPT_LAST+1)
  22 #define OPT_SPLIT           (BV_OPT_LAST+2)
  23 #define OPT_MIN             (BV_OPT_LAST+3)
  24
  25 struct argp_option argp_options[] = {
  26     { "split",              OPT_SPLIT,  "int" },
  27     { "variables",          OPT_VARS,   "list", 0,
  28         "comma separated list of variables over which to maximize" },
  29     { "verbose",            'v',        0,      0, },
  30     { "minimize",           OPT_MIN,    0, 0,   "minimize instead of maximize"},
  31     { 0 }
  32 };
  33
  34 struct options {
  35     struct convert_options   convert;
  36     struct verify_options    verify;
  37     char* var_list;
  38     int verbose;
  39     int split;
  40     int minimize;
  41 };
  42
  43 static error_t parse_opt(int key, char *arg, struct argp_state *state)
  44 {
  45     struct options *options = (struct options*) state->input;
  46
  47     switch (key) {
  48     case ARGP_KEY_INIT:
  49         state->child_inputs[0] = &options->convert;
  50         state->child_inputs[1] = &options->verify;
  51         state->child_inputs[2] = options->verify.barvinok;
  52         options->var_list = NULL;
  53         options->verbose = 0;
  54         options->split = 0;
  55         options->minimize = 0;
  56         break;
  57     case 'v':
  58         options->verbose = 1;
  59         break;
  60     case OPT_VARS:
  61         options->var_list = strdup(arg);
  62         break;
  63     case OPT_SPLIT:
  64         options->split = atoi(arg);
  65         break;
  66     case OPT_MIN:
  67         options->minimize = 1;
  68         break;
  69     default:
  70         return ARGP_ERR_UNKNOWN;
  71     }
  72     return 0;
  73 }
  74
  75
  76 #define ALLOC(type) (type*)malloc(sizeof(type))
  77 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
  78
  79 enum token_type { TOKEN_UNKNOWN = 256, TOKEN_VALUE, TOKEN_IDENT, TOKEN_GE,
  80                   TOKEN_UNION, TOKEN_VARS };
  81
  82 struct token {
  83     enum token_type  type;
  84
  85     int line;
  86     int col;
  87
  88     union {
  89         Value   v;
  90         char    *s;
  91     } u;
  92 };
  93
  94 static struct token *token_new(int line, int col)
  95 {
  96     struct token *tok = ALLOC(struct token);
  97     tok->line = line;
  98     tok->col = col;
  99     return tok;
 100 }
 101
 102 void token_free(struct token *tok)
 103 {
 104     if (tok->type == TOKEN_VALUE)
 105         value_clear(tok->u.v);
 106     else if (tok->type == TOKEN_IDENT)
 107         free(tok->u.s);
 108     free(tok);
 109 }
 110
 111 struct stream {
 112     FILE            *file;
 113     int             line;
 114     int             col;
 115     int             eof;
 116
 117     char            *buffer;
 118     size_t          size;
 119     size_t          len;
 120     int             c;
 121
 122     struct token    *tokens[5];
 123     int             n_token;
 124 };
 125
 126 static struct stream* stream_new(FILE *file)
 127 {
 128     int i;
 129     struct stream *s = ALLOC(struct stream);
 130     s->file = file;
 131     s->size = 256;
 132     s->buffer = (char*)malloc(s->size);
 133     s->len = 0;
 134     s->line = 1;
 135     s->col = 0;
 136     s->eof = 0;
 137     s->c = -1;
 138     for (i = 0; i < 5; ++i)
 139         s->tokens[i] = NULL;
 140     s->n_token = 0;
 141     return s;
 142 }
 143
 144 static void stream_free(struct stream *s)
 145 {
 146     free(s->buffer);
 147     assert(s->n_token == 0);
 148     free(s);
 149 }
 150
 151 static int stream_getc(struct stream *s)
 152 {
 153     int c;
 154     if (s->eof)
 155         return -1;
 156     c = fgetc(s->file);
 157     if (c == -1)
 158         s->eof = 1;
 159     if (s->c != -1) {
 160         if (s->c == '\n') {
 161             s->line++;
 162             s->col = 0;
 163         } else
 164             s->col++;
 165     }
 166     s->c = c;
 167     return c;
 168 }
 169
 170 static void stream_ungetc(struct stream *s, int c)
 171 {
 172     ungetc(c, s->file);
 173     s->c = -1;
 174 }
 175
 176 static void stream_push_char(struct stream *s, int c)
 177 {
 178     if (s->len >= s->size) {
 179         s->size = (3*s->size)/2;
 180         s->buffer = (char*)realloc(s->buffer, s->size);
 181     }
 182     s->buffer[s->len++] = c;
 183 }
 184
 185 static struct token *stream_push_token(struct stream *s, struct token *tok)
 186 {
 187     assert(s->n_token < 5);
 188     s->tokens[s->n_token++] = tok;
 189 }
 190
 191 static struct token *stream_next_token(struct stream *s)
 192 {
 193     int c;
 194     struct token *tok;
 195     int line, col;
 196
 197     if (s->n_token)
 198         return s->tokens[--s->n_token];
 199
 200     s->len = 0;
 201
 202     /* skip spaces */
 203     while ((c = stream_getc(s)) != -1 && isspace(c))
 204         /* nothing */
 205         ;
 206
 207     line = s->line;
 208     col = s->col;
 209
 210     if (c == -1)
 211         return NULL;
 212     if (c == '(' ||
 213         c == ')' ||
 214         c == '+' ||
 215         c == '/' ||
 216         c == '*' ||
 217         c == '^' ||
 218         c == '=' ||
 219         c == ',' ||
 220         c == '_' ||
 221         c == '[' ||
 222         c == ']' ||
 223         c == '{' ||
 224         c == '}') {
 225         tok = token_new(line, col);
 226         tok->type = (token_type)c;
 227         return tok;
 228     }
 229     if (c == '-' || isdigit(c)) {
 230         tok = token_new(line, col);
 231         tok->type = TOKEN_VALUE;
 232         value_init(tok->u.v);
 233         stream_push_char(s, c);
 234         while ((c = stream_getc(s)) != -1 && isdigit(c))
 235             stream_push_char(s, c);
 236         if (c != -1)
 237             stream_ungetc(s, c);
 238         if (s->len == 1 && s->buffer[0] == '-')
 239             value_set_si(tok->u.v, -1);
 240         else {
 241             stream_push_char(s, '\0');
 242             mpz_set_str(tok->u.v, s->buffer, 0);
 243         }
 244         return tok;
 245     }
 246     if (c == '#' || isalpha(c)) {
 247         tok = token_new(line, col);
 248         stream_push_char(s, c);
 249         while ((c = stream_getc(s)) != -1 && isalpha(c))
 250             stream_push_char(s, c);
 251         if (c != -1)
 252             stream_ungetc(s, c);
 253         stream_push_char(s, '\0');
 254         if (!strcmp(s->buffer, "#variables")) {
 255             tok->type = TOKEN_VARS;
 256         } else if (s->buffer[0] == '#') {
 257             tok->type = TOKEN_UNKNOWN;
 258         } else if (!strcmp(s->buffer, "UNION")) {
 259             tok->type = TOKEN_UNION;
 260         } else {
 261             tok->type = TOKEN_IDENT;
 262             tok->u.s = strdup(s->buffer);
 263         }
 264         return tok;
 265     }
 266     if (c == '>') {
 267         if ((c = stream_getc(s)) == '=') {
 268             tok = token_new(line, col);
 269             tok->type = TOKEN_GE;
 270             return tok;
 271         }
 272         if (c != -1)
 273             stream_ungetc(s, c);
 274     }
 275
 276     tok = token_new(line, col);
 277     tok->type = TOKEN_UNKNOWN;
 278     return tok;
 279 }
 280
 281 void stream_error(struct stream *s, struct token *tok, char *msg)
 282 {
 283     int line = tok ? tok->line : s->line;
 284     int col = tok ? tok->col : s->col;
 285     fprintf(stderr, "syntax error (%d, %d): %s\n", line, col, msg);
 286     if (tok) {
 287         if (tok->type < 256)
 288             fprintf(stderr, "got '%c'\n", tok->type);
 289     }
 290 }
 291
 292 struct parameter {
 293     char                *name;
 294     int                  pos;
 295     struct parameter    *next;
 296 };
 297
 298 struct parameter *parameter_new(char *name, int pos, struct parameter *next)
 299 {
 300     struct parameter *p = ALLOC(struct parameter);
 301     p->name = strdup(name);
 302     p->pos = pos;
 303     p->next = next;
 304     return p;
 305 }
 306
 307 static int parameter_pos(struct parameter **p, char *s)
 308 {
 309     int pos = *p ? (*p)->pos+1 : 0;
 310     struct parameter *q;
 311
 312     for (q = *p; q; q = q->next) {
 313         if (strcmp(q->name, s) == 0)
 314             break;
 315     }
 316     if (q)
 317         pos = q->pos;
 318     else
 319         *p = parameter_new(s, pos, *p);
 320     return pos;
 321 }
 322
 323 static int optional_power(struct stream *s)
 324 {
 325     int pow;
 326     struct token *tok;
 327     tok = stream_next_token(s);
 328     if (!tok)
 329         return 1;
 330     if (tok->type != '^') {
 331         stream_push_token(s, tok);
 332         return 1;
 333     }
 334     token_free(tok);
 335     tok = stream_next_token(s);
 336     if (!tok || tok->type != TOKEN_VALUE) {
 337         stream_error(s, tok, "expecting exponent");
 338         if (tok)
 339             stream_push_token(s, tok);
 340         return 1;
 341     }
 342     pow = VALUE_TO_INT(tok->u.v);
 343     token_free(tok);
 344     return pow;
 345 }
 346
 347 static evalue *evalue_read_factor(struct stream *s, struct parameter **p);
 348 static evalue *evalue_read_term(struct stream *s, struct parameter **p);
 349
 350 static evalue *create_fract_like(struct stream *s, evalue *arg, enode_type type,
 351                                  struct parameter **p)
 352 {
 353     evalue *e;
 354     int pow;
 355     pow = optional_power(s);
 356
 357     e = ALLOC(evalue);
 358     value_init(e->d);
 359     e->x.p = new_enode(type, pow+2, -1);
 360     value_clear(e->x.p->arr[0].d);
 361     e->x.p->arr[0] = *arg;
 362     free(arg);
 363     evalue_set_si(&e->x.p->arr[1+pow], 1, 1);
 364     while (--pow >= 0)
 365         evalue_set_si(&e->x.p->arr[1+pow], 0, 1);
 366
 367     return e;
 368 }
 369
 370 static evalue *read_fract(struct stream *s, struct token *tok, struct parameter **p)
 371 {
 372     evalue *arg;
 373
 374     tok = stream_next_token(s);
 375     assert(tok);
 376     assert(tok->type == '{');
 377
 378     token_free(tok);
 379     arg = evalue_read_term(s, p);
 380     tok = stream_next_token(s);
 381     if (!tok || tok->type != '}') {
 382         stream_error(s, tok, "expecting \"}\"");
 383         if (tok)
 384             stream_push_token(s, tok);
 385     } else
 386         token_free(tok);
 387
 388     return create_fract_like(s, arg, fractional, p);
 389 }
 390
 391 static evalue *read_periodic(struct stream *s, struct parameter **p)
 392 {
 393     evalue **list;
 394     int len;
 395     int n;
 396     evalue *e = NULL;
 397
 398     struct token *tok;
 399     tok = stream_next_token(s);
 400     assert(tok && tok->type == '[');
 401     token_free(tok);
 402
 403     len = 100;
 404     list = (evalue **)malloc(len * sizeof(evalue *));
 405     n = 0;
 406
 407     for (;;) {
 408         evalue *e = evalue_read_term(s, p);
 409         if (!e) {
 410             stream_error(s, NULL, "missing argument or list element");
 411             goto out;
 412         }
 413         if (n >= len) {
 414             len = (3*len)/2;
 415             list = (evalue **)realloc(list, len * sizeof(evalue *));
 416         }
 417         list[n++] = e;
 418
 419         tok = stream_next_token(s);
 420         if (!tok) {
 421             stream_error(s, NULL, "unexpected EOF");
 422             goto out;
 423         }
 424         if (tok->type != ',')
 425             break;
 426         token_free(tok);
 427     }
 428
 429     if (tok->type != ']') {
 430         stream_error(s, tok, "expecting \"]\"");
 431         stream_push_token(s, tok);
 432         goto out;
 433     }
 434
 435     token_free(tok);
 436
 437     tok = stream_next_token(s);
 438     if (tok->type == '_') {
 439         int pos;
 440         token_free(tok);
 441         tok = stream_next_token(s);
 442         if (!tok || tok->type != TOKEN_IDENT) {
 443             stream_error(s, tok, "expecting identifier");
 444             if (tok)
 445                 stream_push_token(s, tok);
 446             goto out;
 447         }
 448         e = ALLOC(evalue);
 449         value_init(e->d);
 450         pos = parameter_pos(p, tok->u.s);
 451         token_free(tok);
 452         e->x.p = new_enode(periodic, n, pos+1);
 453         while (--n >= 0) {
 454             value_clear(e->x.p->arr[n].d);
 455             e->x.p->arr[n] = *list[n];
 456             free(list[n]);
 457         }
 458     } else if (n == 1) {
 459         stream_push_token(s, tok);
 460         e = create_fract_like(s, list[0], flooring, p);
 461         n = 0;
 462     } else {
 463         stream_error(s, tok, "unexpected token");
 464         stream_push_token(s, tok);
 465     }
 466
 467 out:
 468     while (--n >= 0) {
 469         free_evalue_refs(list[n]);
 470         free(list[n]);
 471     }
 472     free(list);
 473     return e;
 474 }
 475
 476 static evalue *evalue_read_factor(struct stream *s, struct parameter **p)
 477 {
 478     struct token *tok;
 479     evalue *e = NULL;
 480
 481     tok = stream_next_token(s);
 482     if (!tok)
 483         return NULL;
 484
 485     if (tok->type == '(') {
 486         token_free(tok);
 487         e = evalue_read_term(s, p);
 488         tok = stream_next_token(s);
 489         if (!tok || tok->type != ')') {
 490             stream_error(s, tok, "expecting \")\"");
 491             if (tok)
 492                 stream_push_token(s, tok);
 493         } else
 494             token_free(tok);
 495     } else if (tok->type == TOKEN_VALUE) {
 496         e = ALLOC(evalue);
 497         value_init(e->d);
 498         value_set_si(e->d, 1);
 499         value_init(e->x.n);
 500         value_assign(e->x.n, tok->u.v);
 501         token_free(tok);
 502         tok = stream_next_token(s);
 503         if (tok && tok->type == '/') {
 504             token_free(tok);
 505             tok = stream_next_token(s);
 506             if (!tok || tok->type != TOKEN_VALUE) {
 507                 stream_error(s, tok, "expecting denominator");
 508                 if (tok)
 509                     stream_push_token(s, tok);
 510                 return NULL;
 511             }
 512             value_assign(e->d, tok->u.v);
 513             token_free(tok);
 514         } else if (tok)
 515             stream_push_token(s, tok);
 516     } else if (tok->type == TOKEN_IDENT) {
 517         int pos = parameter_pos(p, tok->u.s);
 518         int pow = optional_power(s);
 519         token_free(tok);
 520         e = ALLOC(evalue);
 521         value_init(e->d);
 522         e->x.p = new_enode(polynomial, pow+1, pos+1);
 523         evalue_set_si(&e->x.p->arr[pow], 1, 1);
 524         while (--pow >= 0)
 525             evalue_set_si(&e->x.p->arr[pow], 0, 1);
 526     } else if (tok->type == '[') {
 527         stream_push_token(s, tok);
 528         e = read_periodic(s, p);
 529     } else if (tok->type == '{') {
 530         stream_push_token(s, tok);
 531         e = read_fract(s, tok, p);
 532     }
 533
 534     tok = stream_next_token(s);
 535     if (tok && tok->type == '*') {
 536         evalue *e2;
 537         token_free(tok);
 538         e2 = evalue_read_factor(s, p);
 539         if (!e2) {
 540             stream_error(s, NULL, "unexpected EOF");
 541             return NULL;
 542         }
 543         emul(e2, e);
 544         free_evalue_refs(e2);
 545         free(e2);
 546     } else if (tok)
 547         stream_push_token(s, tok);
 548
 549     return e;
 550 }
 551
 552 static evalue *evalue_read_term(struct stream *s, struct parameter **p)
 553 {
 554     struct token *tok;
 555     evalue *e = NULL;
 556
 557     e = evalue_read_factor(s, p);
 558     if (!e)
 559         return NULL;
 560
 561     tok = stream_next_token(s);
 562     if (!tok)
 563         return e;
 564
 565     if (tok->type == '+') {
 566         evalue *e2;
 567         token_free(tok);
 568         e2 = evalue_read_term(s, p);
 569         if (!e2) {
 570             stream_error(s, NULL, "unexpected EOF");
 571             return NULL;
 572         }
 573         eadd(e2, e);
 574         free_evalue_refs(e2);
 575         free(e2);
 576     } else
 577         stream_push_token(s, tok);
 578
 579     return e;
 580 }
 581
 582 struct constraint {
 583     int                  type;
 584     Vector              *v;
 585     struct constraint   *next;
 586     struct constraint   *union_next;
 587 };
 588
 589 static struct constraint *constraint_new()
 590 {
 591     struct constraint *c = ALLOC(struct constraint);
 592     c->type = -1;
 593     c->v = Vector_Alloc(16);
 594     c->next = NULL;
 595     c->union_next = NULL;
 596     return c;
 597 }
 598
 599 static void constraint_free(struct constraint *c)
 600 {
 601     while (c) {
 602         struct constraint *next = c->next ? c->next : c->union_next;
 603         Vector_Free(c->v);
 604         free(c);
 605         c = next;
 606     }
 607 }
 608
 609 static void constraint_extend(struct constraint *c, int pos)
 610 {
 611     Vector *v;
 612     if (pos < c->v->Size)
 613         return;
 614
 615     v = Vector_Alloc((3*c->v->Size)/2);
 616     Vector_Copy(c->v->p, v->p, c->v->Size);
 617     Vector_Free(c->v);
 618     c->v = v;
 619 }
 620
 621 static int evalue_read_constraint(struct stream *s, struct parameter **p,
 622                                   struct constraint **constraints,
 623                                   struct constraint *union_next)
 624 {
 625     struct token *tok;
 626     struct constraint *c = NULL;
 627
 628     while ((tok = stream_next_token(s))) {
 629         struct token *tok2;
 630         int pos;
 631         if (tok->type == '+')
 632             token_free(tok);
 633         else if (tok->type == TOKEN_IDENT) {
 634             if (!c)
 635                 c = constraint_new();
 636             pos = parameter_pos(p, tok->u.s);
 637             constraint_extend(c, 1+pos);
 638             value_set_si(c->v->p[1+pos], 1);
 639             token_free(tok);
 640         } else if (tok->type == TOKEN_VALUE) {
 641             if (!c)
 642                 c = constraint_new();
 643             tok2 = stream_next_token(s);
 644             if (tok2 && tok2->type == TOKEN_IDENT) {
 645                 pos = parameter_pos(p, tok2->u.s);
 646                 constraint_extend(c, 1+pos);
 647                 value_assign(c->v->p[1+pos], tok->u.v);
 648                 token_free(tok);
 649                 token_free(tok2);
 650             } else {
 651                 if (tok2)
 652                     stream_push_token(s, tok2);
 653                 value_assign(c->v->p[0], tok->u.v);
 654                 token_free(tok);
 655             }
 656         } else if (tok->type == TOKEN_GE || tok->type == '=') {
 657             int type = tok->type == TOKEN_GE;
 658             token_free(tok);
 659             tok = stream_next_token(s);
 660             if (!tok || tok->type != TOKEN_VALUE || value_notzero_p(tok->u.v)) {
 661                 stream_error(s, tok, "expecting \"0\"");
 662                 if (tok)
 663                     stream_push_token(s, tok);
 664                 *constraints = NULL;
 665             } else {
 666                 c->type = type;
 667                 c->next = *constraints;
 668                 c->union_next = union_next;
 669                 *constraints = c;
 670                 token_free(tok);
 671             }
 672             break;
 673         } else {
 674             if (!c)
 675                 stream_push_token(s, tok);
 676             else {
 677                 stream_error(s, tok, "unexpected token");
 678                 *constraints = NULL;
 679             }
 680             return 0;
 681         }
 682     }
 683     return tok != NULL;
 684 }
 685
 686 static struct constraint *evalue_read_domain(struct stream *s, struct parameter **p,
 687                                              unsigned MaxRays)
 688 {
 689     struct constraint *constraints = NULL;
 690     struct constraint *union_next = NULL;
 691     struct token *tok;
 692     int line;
 693
 694     tok = stream_next_token(s);
 695     if (!tok)
 696         return NULL;
 697     stream_push_token(s, tok);
 698
 699     line = tok->line;
 700     while (evalue_read_constraint(s, p, &constraints, union_next)) {
 701         tok = stream_next_token(s);
 702         if (tok) {
 703             if (tok->type == TOKEN_UNION) {
 704                 token_free(tok);
 705                 tok = stream_next_token(s);
 706                 if (!tok) {
 707                     stream_error(s, NULL, "unexpected EOF");
 708                     return constraints;
 709                 }
 710                 stream_push_token(s, tok);
 711                 union_next = constraints;
 712                 constraints = NULL;
 713             } else {
 714                 union_next = NULL;
 715                 stream_push_token(s, tok);
 716                 /* empty line separates domain from evalue */
 717                 if (tok->line > line+1)
 718                     break;
 719             }
 720             line = tok->line;
 721         }
 722     }
 723     return constraints;
 724 }
 725
 726 struct section {
 727     struct constraint   *constraints;
 728     evalue              *e;
 729     struct section      *next;
 730 };
 731
 732 static char **extract_parameters(struct parameter *p, unsigned *nparam)
 733 {
 734     int i;
 735     char **params;
 736
 737     *nparam = p ? p->pos+1 : 0;
 738     params = ALLOCN(char *, *nparam);
 739     for (i = 0; i < *nparam; ++i) {
 740         struct parameter *next = p->next;
 741         params[p->pos] = p->name;
 742         free(p);
 743         p = next;
 744     }
 745     return params;
 746 }
 747
 748 static Polyhedron *constraints2domain(struct constraint *constraints,
 749                                       unsigned nparam, unsigned MaxRays)
 750 {
 751     Polyhedron *D;
 752     Matrix *M;
 753     int n;
 754     struct constraint *c;
 755     struct constraint *union_next = NULL;
 756
 757     for (n = 0, c = constraints; c; ++n, c = c->next)
 758         ;
 759     M = Matrix_Alloc(n, 1+nparam+1);
 760     while (--n >= 0) {
 761         struct constraint *next = constraints->next;
 762         union_next = constraints->union_next;
 763         Vector_Copy(constraints->v->p+1, M->p[n]+1, nparam);
 764         if (constraints->type)
 765             value_set_si(M->p[n][0], 1);
 766         value_assign(M->p[n][1+nparam], constraints->v->p[0]);
 767         constraints->next = NULL;
 768         constraints->union_next = NULL;
 769         constraint_free(constraints);
 770         constraints = next;
 771     }
 772     D = Constraints2Polyhedron(M, MaxRays);
 773     Matrix_Free(M);
 774
 775     if (union_next)
 776         D = DomainConcat(D, constraints2domain(union_next, nparam, MaxRays));
 777     return D;
 778 }
 779
 780 static evalue *evalue_read_partition(struct stream *s, struct parameter *p,
 781                                      char ***ppp,
 782                                      unsigned *nparam, unsigned MaxRays)
 783 {
 784     struct section *part = NULL;
 785     struct constraint *constraints;
 786     evalue *e = NULL;
 787     int m = 0;
 788
 789     while ((constraints = evalue_read_domain(s, &p, MaxRays))) {
 790         evalue *e = evalue_read_term(s, &p);
 791         if (!e) {
 792             stream_error(s, NULL, "missing evalue");
 793             break;
 794         }
 795         struct section *sect = ALLOC(struct section);
 796         sect->constraints = constraints;
 797         sect->e = e;
 798         sect->next = part;
 799         part = sect;
 800         ++m;
 801     }
 802
 803     if (part) {
 804         Polyhedron *D;
 805         int j;
 806
 807         *ppp = extract_parameters(p, nparam);
 808         e = ALLOC(evalue);
 809         value_init(e->d);
 810         e->x.p = new_enode(partition, 2*m, *nparam);
 811
 812         for (j = 0; j < m; ++j) {
 813             int n;
 814             struct section *next = part->next;
 815             constraints = part->constraints;
 816             D = constraints2domain(part->constraints, *nparam, MaxRays);
 817             EVALUE_SET_DOMAIN(e->x.p->arr[2*j], D);
 818             value_clear(e->x.p->arr[2*j+1].d);
 819             e->x.p->arr[2*j+1] = *part->e;
 820             free(part->e);
 821             free(part);
 822             part = next;
 823         }
 824     }
 825     return e;
 826 }
 827
 828 static evalue *evalue_read(FILE *in, char *var_list, char ***ppp, unsigned *nvar,
 829                            unsigned *nparam, unsigned MaxRays)
 830 {
 831     struct stream *s = stream_new(in);
 832     struct token *tok;
 833     evalue *e;
 834     struct parameter *p = NULL;
 835     char *next;
 836     int nv;
 837
 838     if (var_list) {
 839         while ((next = strchr(var_list, ','))) {
 840             *next = '\0';
 841             if (next > var_list)
 842                 parameter_pos(&p, var_list);
 843             *next = ',';
 844             var_list = next+1;
 845         }
 846         if (strlen(var_list) > 0)
 847             parameter_pos(&p, var_list);
 848         nv = p ? p->pos+1 : 0;
 849     } else
 850         nv = -1;
 851
 852     if (!(tok = stream_next_token(s)))
 853         return NULL;
 854
 855     if (tok->type == TOKEN_VARS) {
 856         token_free(tok);
 857         for (;;) {
 858             tok = stream_next_token(s);
 859             if (!tok || tok->type != TOKEN_IDENT) {
 860                 stream_error(s, tok, "expecting identifier");
 861                 break;
 862             }
 863             if (nv == -1)
 864                 parameter_pos(&p, tok->u.s);
 865             token_free(tok);
 866             tok = stream_next_token(s);
 867             if (!tok || tok->type != ',')
 868                 break;
 869             token_free(tok);
 870         }
 871         if (!tok)
 872             return NULL;
 873         if (nv = -1)
 874             nv = p ? p->pos+1 : 0;
 875     }
 876
 877     if (tok->type == TOKEN_VALUE) {
 878         struct token *tok2 = stream_next_token(s);
 879         if (tok2)
 880             stream_push_token(s, tok2);
 881         stream_push_token(s, tok);
 882         if (tok2 && (tok2->type == TOKEN_IDENT || tok2->type == TOKEN_GE))
 883             e = evalue_read_partition(s, p, ppp, nparam, MaxRays);
 884         else {
 885             e = evalue_read_term(s, &p);
 886             *ppp = extract_parameters(p, nparam);
 887         }
 888     } else if (tok->type == TOKEN_IDENT) {
 889         stream_push_token(s, tok);
 890         e = evalue_read_partition(s, p, ppp, nparam, MaxRays);
 891     } else {
 892         stream_error(s, tok, "unexpected token");
 893         *nparam = nv == -1 ? 0 : nv;
 894         e = NULL;
 895     }
 896     stream_free(s);
 897     if (nv == -1)
 898         *nvar = *nparam;
 899     else
 900         *nvar = nv;
 901     *nparam -= *nvar;
 902     return e;
 903 }
 904
 905 static int check_poly_max(const struct check_poly_data *data,
 906                           int nparam, Value *z,
 907                           const struct verify_options *options);
 908
 909 struct check_poly_max_data : public check_poly_data {
 910     Polyhedron          **S;
 911     evalue               *EP;
 912     piecewise_lst        *pl;
 913
 914     check_poly_max_data(Value *z, evalue *EP, piecewise_lst *pl) :
 915                         EP(EP), pl(pl) {
 916         this->z = z;
 917         this->check = check_poly_max;
 918     }
 919 };
 920
 921 static void optimum(Polyhedron *S, int pos, const check_poly_max_data *data,
 922                     Value *opt, bool& found,
 923                     const struct verify_options *options)
 924 {
 925     if (!S) {
 926         Value c;
 927         value_init(c);
 928         value_set_double(c, compute_evalue(data->EP, data->z+1)+.25);
 929         if (!found) {
 930             value_assign(*opt, c);
 931             found = true;
 932         } else {
 933             if (options->barvinok->bernstein_optimize == BV_BERNSTEIN_MAX) {
 934                 if (value_gt(c, *opt))
 935                     value_assign(*opt, c);
 936             } else {
 937                 if (value_lt(c, *opt))
 938                     value_assign(*opt, c);
 939             }
 940         }
 941         value_clear(c);
 942     } else {
 943         Value LB, UB;
 944         int ok;
 945         value_init(LB);
 946         value_init(UB);
 947         ok = !(lower_upper_bounds(1+pos, S, data->z, &LB, &UB));
 948         assert(ok);
 949         for (; value_le(LB, UB); value_increment(LB, LB)) {
 950             value_assign(data->z[1+pos], LB);
 951             optimum(S->next, pos+1, data, opt, found, options);
 952         }
 953         value_set_si(data->z[1+pos], 0);
 954         value_clear(LB);
 955         value_clear(UB);
 956     }
 957 }
 958
 959 static void optimum(const check_poly_max_data *data, Value *opt,
 960                     const struct verify_options *options)
 961 {
 962     bool found = false;
 963     for (int i = 0; i < data->EP->x.p->size/2; ++i)
 964         if (!emptyQ2(data->S[i]))
 965             optimum(data->S[i], 0, data, opt, found, options);
 966     assert(found);
 967 }
 968
 969 static int check_poly_max(const struct check_poly_data *data,
 970                           int nparam, Value *z,
 971                           const struct verify_options *options)
 972 {
 973     int k;
 974     int ok;
 975     const check_poly_max_data *max_data;
 976     max_data = static_cast<const check_poly_max_data *>(data);
 977     char *minmax;
 978     Value m, n, d;
 979     value_init(m);
 980     value_init(n);
 981     value_init(d);
 982
 983     if (options->barvinok->bernstein_optimize == BV_BERNSTEIN_MAX)
 984         minmax = "max";
 985     else
 986         minmax = "min";
 987
 988     max_data->pl->evaluate(nparam, z, &n, &d);
 989     if (options->barvinok->bernstein_optimize == BV_BERNSTEIN_MAX)
 990         mpz_fdiv_q(m, n, d);
 991     else
 992         mpz_cdiv_q(m, n, d);
 993
 994     if (options->print_all) {
 995         printf("%s(", minmax);
 996         value_print(stdout, VALUE_FMT, z[0]);
 997         for (k = 1; k < nparam; ++k) {
 998             printf(", ");
 999             value_print(stdout, VALUE_FMT, z[k]);
1000         }
1001         printf(") = ");
1002         value_print(stdout, VALUE_FMT, n);
1003         if (value_notone_p(d)) {
1004             printf("/");
1005             value_print(stdout, VALUE_FMT, d);
1006         }
1007         printf(" (");
1008         value_print(stdout, VALUE_FMT, m);
1009         printf(")");
1010     }
1011
1012     optimum(max_data, &n, options);
1013
1014     if (options->barvinok->bernstein_optimize == BV_BERNSTEIN_MAX)
1015         ok = value_ge(m, n);
1016     else
1017         ok = value_le(m, n);
1018
1019     if (options->print_all) {
1020         printf(", %s(EP) = ", minmax);
1021         value_print(stdout, VALUE_FMT, n);
1022         printf(". ");
1023     }
1024
1025     if (!ok) {
1026         printf("\n");
1027         fflush(stdout);
1028         fprintf(stderr, "Error !\n");
1029         fprintf(stderr, "%s(", minmax);
1030         value_print(stderr, VALUE_FMT, z[0]);
1031         for (k = 1; k < nparam; ++k) {
1032             fprintf(stderr,", ");
1033             value_print(stderr, VALUE_FMT, z[k]);
1034         }
1035         fprintf(stderr, ") should be ");
1036         if (options->barvinok->bernstein_optimize == BV_BERNSTEIN_MAX)
1037             fprintf(stderr, "greater");
1038         else
1039             fprintf(stderr, "smaller");
1040         fprintf(stderr, " than or equal to ");
1041         value_print(stderr, VALUE_FMT, n);
1042         fprintf(stderr, ", while pl eval gives ");
1043         value_print(stderr, VALUE_FMT, m);
1044         fprintf(stderr, ".\n");
1045         cerr << *max_data->pl << endl;
1046     } else if (options->print_all)
1047         printf("OK.\n");
1048
1049     value_clear(m);
1050     value_clear(n);
1051     value_clear(d);
1052
1053     return ok;
1054 }
1055
1056 static int verify(Polyhedron *D, piecewise_lst *pl, evalue *EP,
1057                   unsigned nvar, unsigned nparam, Vector *p,
1058                   struct verify_options *options)
1059 {
1060     Polyhedron *CS, *S;
1061     unsigned MaxRays = options->barvinok->MaxRays;
1062     assert(value_zero_p(EP->d));
1063     assert(EP->x.p->type == partition);
1064     int ok = 1;
1065
1066     CS = check_poly_context_scan(NULL, &D, D->Dimension, options);
1067
1068     check_poly_init(D, options);
1069
1070     if (!(CS && emptyQ2(CS))) {
1071         check_poly_max_data data(p->p, EP, pl);
1072         data.S = ALLOCN(Polyhedron *, EP->x.p->size/2);
1073         for (int i = 0; i < EP->x.p->size/2; ++i) {
1074             Polyhedron *A = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
1075             data.S[i] = Polyhedron_Scan(A, D, MaxRays & POL_NO_DUAL ? 0 : MaxRays);
1076         }
1077         ok = check_poly(CS, &data, nparam, 0, p->p+1+nvar, options);
1078         for (int i = 0; i < EP->x.p->size/2; ++i)
1079             Domain_Free(data.S[i]);
1080         free(data.S);
1081     }
1082
1083     if (!options->print_all)
1084         printf("\n");
1085
1086     if (CS) {
1087         Domain_Free(CS);
1088         Domain_Free(D);
1089     }
1090
1091     return ok;
1092 }
1093
1094 static int verify(piecewise_lst *pl, evalue *EP, unsigned nvar, unsigned nparam,
1095                   struct verify_options *options)
1096 {
1097     Vector *p;
1098
1099     p = Vector_Alloc(nvar+nparam+2);
1100     value_set_si(p->p[nvar+nparam+1], 1);
1101
1102     for (int i = 0; i < pl->list.size(); ++i) {
1103         int ok = verify(pl->list[i].first, pl, EP, nvar, nparam, p, options);
1104         if (!ok && !options->continue_on_error)
1105             break;
1106     }
1107
1108     Vector_Free(p);
1109
1110     return 0;
1111 }
1112
1113 static int optimize(evalue *EP, char **all_vars, unsigned nvar, unsigned nparam,
1114                     struct options *options)
1115 {
1116     Polyhedron *U;
1117     piecewise_lst *pl = NULL;
1118     U = Universe_Polyhedron(nparam);
1119     int print_solution = 1;
1120     int result = 0;
1121
1122     exvector params;
1123     params = constructParameterVector(all_vars+nvar, nparam);
1124
1125     if (options->verify.verify) {
1126         verify_options_set_range(&options->verify, nvar+nparam);
1127         if (!options->verbose)
1128             print_solution = 0;
1129     }
1130
1131     if (options->minimize)
1132         options->verify.barvinok->bernstein_optimize = BV_BERNSTEIN_MIN;
1133     else
1134         options->verify.barvinok->bernstein_optimize = BV_BERNSTEIN_MAX;
1135     pl = evalue_bernstein_coefficients(NULL, EP, U, params,
1136                                        options->verify.barvinok);
1137     assert(pl);
1138     if (options->minimize)
1139         pl->minimize();
1140     else
1141         pl->maximize();
1142     if (print_solution)
1143         cout << *pl << endl;
1144     if (options->verify.verify)
1145         result = verify(pl, EP, nvar, nparam, &options->verify);
1146     delete pl;
1147
1148     Polyhedron_Free(U);
1149
1150     return result;
1151 }
1152
1153 int main(int argc, char **argv)
1154 {
1155     evalue *EP;
1156     char **all_vars = NULL;
1157     unsigned nvar;
1158     unsigned nparam;
1159     struct options options;
1160     struct barvinok_options *bv_options = barvinok_options_new_with_defaults();
1161     static struct argp_child argp_children[] = {
1162         { &convert_argp,        0,      "input conversion",     1 },
1163         { &verify_argp,         0,      "verification",         2 },
1164         { &barvinok_argp,       0,      "barvinok options",     3 },
1165         { 0 }
1166     };
1167     static struct argp argp = { argp_options, parse_opt, 0, 0, argp_children };
1168     int result = 0;
1169
1170     options.verify.barvinok = bv_options;
1171     set_program_name(argv[0]);
1172     argp_parse(&argp, argc, argv, 0, 0, &options);
1173
1174     EP = evalue_read(stdin, options.var_list, &all_vars, &nvar, &nparam,
1175                      bv_options->MaxRays);
1176     assert(EP);
1177
1178     if (options.split)
1179         evalue_split_periods(EP, options.split, bv_options->MaxRays);
1180
1181     evalue_convert(EP, &options.convert, options.verbose, nparam, all_vars);
1182
1183     if (EVALUE_IS_ZERO(*EP))
1184         print_evalue(stdout, EP, all_vars);
1185     else
1186         result = optimize(EP, all_vars, nvar, nparam, &options);
1187
1188     free_evalue_refs(EP);
1189     free(EP);
1190
1191     if (options.var_list)
1192         free(options.var_list);
1193     Free_ParamNames(all_vars, nvar+nparam);
1194     barvinok_options_free(bv_options);
1195     return result;
1196 }