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