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