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_NE, 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 && isalnum(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     if (c == '!') {
 276         if ((c = stream_getc(s)) == '=') {
 277             tok = token_new(line, col);
 278             tok->type = TOKEN_NE;
 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 *create_relation(evalue *arg, int ne)
 380 {
 381     evalue *e;
 382
 383     e = ALLOC(evalue);
 384     value_init(e->d);
 385     e->x.p = new_enode(relation, 2+ne, 0);
 386     value_clear(e->x.p->arr[0].d);
 387     e->x.p->arr[0] = *arg;
 388     free(arg);
 389     if (ne)
 390         evalue_set_si(&e->x.p->arr[1], 0, 1);
 391     evalue_set_si(&e->x.p->arr[1+ne], 1, 1);
 392
 393     return e;
 394 }
 395
 396 static evalue *read_fract(struct stream *s, struct token *tok, struct parameter **p)
 397 {
 398     evalue *arg;
 399
 400     tok = stream_next_token(s);
 401     assert(tok);
 402     assert(tok->type == '{');
 403
 404     token_free(tok);
 405     arg = evalue_read_term(s, p);
 406     tok = stream_next_token(s);
 407     if (!tok || tok->type != '}') {
 408         stream_error(s, tok, "expecting \"}\"");
 409         if (tok)
 410             stream_push_token(s, tok);
 411     } else
 412         token_free(tok);
 413
 414     return create_fract_like(s, arg, fractional, p);
 415 }
 416
 417 static evalue *read_periodic(struct stream *s, struct parameter **p)
 418 {
 419     evalue **list;
 420     int len;
 421     int n;
 422     evalue *e = NULL;
 423
 424     struct token *tok;
 425     tok = stream_next_token(s);
 426     assert(tok && tok->type == '[');
 427     token_free(tok);
 428
 429     len = 100;
 430     list = (evalue **)malloc(len * sizeof(evalue *));
 431     n = 0;
 432
 433     for (;;) {
 434         evalue *e = evalue_read_term(s, p);
 435         if (!e) {
 436             stream_error(s, NULL, "missing argument or list element");
 437             goto out;
 438         }
 439         if (n >= len) {
 440             len = (3*len)/2;
 441             list = (evalue **)realloc(list, len * sizeof(evalue *));
 442         }
 443         list[n++] = e;
 444
 445         tok = stream_next_token(s);
 446         if (!tok) {
 447             stream_error(s, NULL, "unexpected EOF");
 448             goto out;
 449         }
 450         if (tok->type != ',')
 451             break;
 452         token_free(tok);
 453     }
 454
 455     if (n == 1 && (tok->type == '=' || tok->type == TOKEN_NE)) {
 456         int ne = tok->type == TOKEN_NE;
 457         token_free(tok);
 458         tok = stream_next_token(s);
 459         if (!tok || tok->type != TOKEN_VALUE) {
 460             stream_error(s, tok, "expecting \"0\"");
 461             if (tok)
 462                 stream_push_token(s, tok);
 463             goto out;
 464         }
 465         token_free(tok);
 466         tok = stream_next_token(s);
 467         if (!tok || tok->type != ']') {
 468             stream_error(s, tok, "expecting \"]\"");
 469             if (tok)
 470                 stream_push_token(s, tok);
 471             goto out;
 472         }
 473         token_free(tok);
 474         e = create_relation(list[0], ne);
 475         n = 0;
 476         goto out;
 477     }
 478
 479     if (tok->type != ']') {
 480         stream_error(s, tok, "expecting \"]\"");
 481         stream_push_token(s, tok);
 482         goto out;
 483     }
 484
 485     token_free(tok);
 486
 487     tok = stream_next_token(s);
 488     if (tok->type == '_') {
 489         int pos;
 490         token_free(tok);
 491         tok = stream_next_token(s);
 492         if (!tok || tok->type != TOKEN_IDENT) {
 493             stream_error(s, tok, "expecting identifier");
 494             if (tok)
 495                 stream_push_token(s, tok);
 496             goto out;
 497         }
 498         e = ALLOC(evalue);
 499         value_init(e->d);
 500         pos = parameter_pos(p, tok->u.s);
 501         token_free(tok);
 502         e->x.p = new_enode(periodic, n, pos+1);
 503         while (--n >= 0) {
 504             value_clear(e->x.p->arr[n].d);
 505             e->x.p->arr[n] = *list[n];
 506             free(list[n]);
 507         }
 508     } else if (n == 1) {
 509         stream_push_token(s, tok);
 510         e = create_fract_like(s, list[0], flooring, p);
 511         n = 0;
 512     } else {
 513         stream_error(s, tok, "unexpected token");
 514         stream_push_token(s, tok);
 515     }
 516
 517 out:
 518     while (--n >= 0) {
 519         free_evalue_refs(list[n]);
 520         free(list[n]);
 521     }
 522     free(list);
 523     return e;
 524 }
 525
 526 static evalue *evalue_read_factor(struct stream *s, struct parameter **p)
 527 {
 528     struct token *tok;
 529     evalue *e = NULL;
 530
 531     tok = stream_next_token(s);
 532     if (!tok)
 533         return NULL;
 534
 535     if (tok->type == '(') {
 536         token_free(tok);
 537         e = evalue_read_term(s, p);
 538         tok = stream_next_token(s);
 539         if (!tok || tok->type != ')') {
 540             stream_error(s, tok, "expecting \")\"");
 541             if (tok)
 542                 stream_push_token(s, tok);
 543         } else
 544             token_free(tok);
 545     } else if (tok->type == TOKEN_VALUE) {
 546         e = ALLOC(evalue);
 547         value_init(e->d);
 548         value_set_si(e->d, 1);
 549         value_init(e->x.n);
 550         value_assign(e->x.n, tok->u.v);
 551         token_free(tok);
 552         tok = stream_next_token(s);
 553         if (tok && tok->type == '/') {
 554             token_free(tok);
 555             tok = stream_next_token(s);
 556             if (!tok || tok->type != TOKEN_VALUE) {
 557                 stream_error(s, tok, "expecting denominator");
 558                 if (tok)
 559                     stream_push_token(s, tok);
 560                 return NULL;
 561             }
 562             value_assign(e->d, tok->u.v);
 563             token_free(tok);
 564         } else if (tok)
 565             stream_push_token(s, tok);
 566     } else if (tok->type == TOKEN_IDENT) {
 567         int pos = parameter_pos(p, tok->u.s);
 568         int pow = optional_power(s);
 569         token_free(tok);
 570         e = ALLOC(evalue);
 571         value_init(e->d);
 572         e->x.p = new_enode(polynomial, pow+1, pos+1);
 573         evalue_set_si(&e->x.p->arr[pow], 1, 1);
 574         while (--pow >= 0)
 575             evalue_set_si(&e->x.p->arr[pow], 0, 1);
 576     } else if (tok->type == '[') {
 577         stream_push_token(s, tok);
 578         e = read_periodic(s, p);
 579     } else if (tok->type == '{') {
 580         stream_push_token(s, tok);
 581         e = read_fract(s, tok, p);
 582     }
 583
 584     tok = stream_next_token(s);
 585     if (tok && tok->type == '*') {
 586         evalue *e2;
 587         token_free(tok);
 588         e2 = evalue_read_factor(s, p);
 589         if (!e2) {
 590             stream_error(s, NULL, "unexpected EOF");
 591             return NULL;
 592         }
 593         emul(e2, e);
 594         free_evalue_refs(e2);
 595         free(e2);
 596     } else if (tok)
 597         stream_push_token(s, tok);
 598
 599     return e;
 600 }
 601
 602 static evalue *evalue_read_term(struct stream *s, struct parameter **p)
 603 {
 604     struct token *tok;
 605     evalue *e = NULL;
 606
 607     e = evalue_read_factor(s, p);
 608     if (!e)
 609         return NULL;
 610
 611     tok = stream_next_token(s);
 612     if (!tok)
 613         return e;
 614
 615     if (tok->type == '+') {
 616         evalue *e2;
 617         token_free(tok);
 618         e2 = evalue_read_term(s, p);
 619         if (!e2) {
 620             stream_error(s, NULL, "unexpected EOF");
 621             return NULL;
 622         }
 623         eadd(e2, e);
 624         free_evalue_refs(e2);
 625         free(e2);
 626     } else
 627         stream_push_token(s, tok);
 628
 629     return e;
 630 }
 631
 632 struct constraint {
 633     int                  type;
 634     Vector              *v;
 635     struct constraint   *next;
 636     struct constraint   *union_next;
 637 };
 638
 639 static struct constraint *constraint_new()
 640 {
 641     struct constraint *c = ALLOC(struct constraint);
 642     c->type = -1;
 643     c->v = Vector_Alloc(16);
 644     c->next = NULL;
 645     c->union_next = NULL;
 646     return c;
 647 }
 648
 649 static void constraint_free(struct constraint *c)
 650 {
 651     while (c) {
 652         struct constraint *next = c->next ? c->next : c->union_next;
 653         Vector_Free(c->v);
 654         free(c);
 655         c = next;
 656     }
 657 }
 658
 659 static void constraint_extend(struct constraint *c, int pos)
 660 {
 661     Vector *v;
 662     if (pos < c->v->Size)
 663         return;
 664
 665     v = Vector_Alloc((3*c->v->Size)/2);
 666     Vector_Copy(c->v->p, v->p, c->v->Size);
 667     Vector_Free(c->v);
 668     c->v = v;
 669 }
 670
 671 static int evalue_read_constraint(struct stream *s, struct parameter **p,
 672                                   struct constraint **constraints,
 673                                   struct constraint *union_next)
 674 {
 675     struct token *tok;
 676     struct constraint *c = NULL;
 677
 678     while ((tok = stream_next_token(s))) {
 679         struct token *tok2;
 680         int pos;
 681         if (tok->type == '+')
 682             token_free(tok);
 683         else if (tok->type == TOKEN_IDENT) {
 684             if (!c)
 685                 c = constraint_new();
 686             pos = parameter_pos(p, tok->u.s);
 687             constraint_extend(c, 1+pos);
 688             value_set_si(c->v->p[1+pos], 1);
 689             token_free(tok);
 690         } else if (tok->type == TOKEN_VALUE) {
 691             if (!c)
 692                 c = constraint_new();
 693             tok2 = stream_next_token(s);
 694             if (tok2 && tok2->type == TOKEN_IDENT) {
 695                 pos = parameter_pos(p, tok2->u.s);
 696                 constraint_extend(c, 1+pos);
 697                 value_assign(c->v->p[1+pos], tok->u.v);
 698                 token_free(tok);
 699                 token_free(tok2);
 700             } else {
 701                 if (tok2)
 702                     stream_push_token(s, tok2);
 703                 value_assign(c->v->p[0], tok->u.v);
 704                 token_free(tok);
 705             }
 706         } else if (tok->type == TOKEN_GE || tok->type == '=') {
 707             int type = tok->type == TOKEN_GE;
 708             token_free(tok);
 709             tok = stream_next_token(s);
 710             if (!tok || tok->type != TOKEN_VALUE || value_notzero_p(tok->u.v)) {
 711                 stream_error(s, tok, "expecting \"0\"");
 712                 if (tok)
 713                     stream_push_token(s, tok);
 714                 *constraints = NULL;
 715             } else {
 716                 c->type = type;
 717                 c->next = *constraints;
 718                 c->union_next = union_next;
 719                 *constraints = c;
 720                 token_free(tok);
 721             }
 722             break;
 723         } else {
 724             if (!c)
 725                 stream_push_token(s, tok);
 726             else {
 727                 stream_error(s, tok, "unexpected token");
 728                 *constraints = NULL;
 729             }
 730             return 0;
 731         }
 732     }
 733     return tok != NULL;
 734 }
 735
 736 static struct constraint *evalue_read_domain(struct stream *s, struct parameter **p,
 737                                              unsigned MaxRays)
 738 {
 739     struct constraint *constraints = NULL;
 740     struct constraint *union_next = NULL;
 741     struct token *tok;
 742     int line;
 743
 744     tok = stream_next_token(s);
 745     if (!tok)
 746         return NULL;
 747     stream_push_token(s, tok);
 748
 749     line = tok->line;
 750     while (evalue_read_constraint(s, p, &constraints, union_next)) {
 751         tok = stream_next_token(s);
 752         if (tok) {
 753             if (tok->type == TOKEN_UNION) {
 754                 token_free(tok);
 755                 tok = stream_next_token(s);
 756                 if (!tok) {
 757                     stream_error(s, NULL, "unexpected EOF");
 758                     return constraints;
 759                 }
 760                 stream_push_token(s, tok);
 761                 union_next = constraints;
 762                 constraints = NULL;
 763             } else {
 764                 union_next = NULL;
 765                 stream_push_token(s, tok);
 766                 /* empty line separates domain from evalue */
 767                 if (tok->line > line+1)
 768                     break;
 769             }
 770             line = tok->line;
 771         }
 772     }
 773     return constraints;
 774 }
 775
 776 struct section {
 777     struct constraint   *constraints;
 778     evalue              *e;
 779     struct section      *next;
 780 };
 781
 782 static char **extract_parameters(struct parameter *p, unsigned *nparam)
 783 {
 784     int i;
 785     char **params;
 786
 787     *nparam = p ? p->pos+1 : 0;
 788     params = ALLOCN(char *, *nparam);
 789     for (i = 0; i < *nparam; ++i) {
 790         struct parameter *next = p->next;
 791         params[p->pos] = p->name;
 792         free(p);
 793         p = next;
 794     }
 795     return params;
 796 }
 797
 798 static Polyhedron *constraints2domain(struct constraint *constraints,
 799                                       unsigned nparam, unsigned MaxRays)
 800 {
 801     Polyhedron *D;
 802     Matrix *M;
 803     int n;
 804     struct constraint *c;
 805     struct constraint *union_next = NULL;
 806
 807     for (n = 0, c = constraints; c; ++n, c = c->next)
 808         ;
 809     M = Matrix_Alloc(n, 1+nparam+1);
 810     while (--n >= 0) {
 811         struct constraint *next = constraints->next;
 812         union_next = constraints->union_next;
 813         Vector_Copy(constraints->v->p+1, M->p[n]+1, nparam);
 814         if (constraints->type)
 815             value_set_si(M->p[n][0], 1);
 816         value_assign(M->p[n][1+nparam], constraints->v->p[0]);
 817         constraints->next = NULL;
 818         constraints->union_next = NULL;
 819         constraint_free(constraints);
 820         constraints = next;
 821     }
 822     D = Constraints2Polyhedron(M, MaxRays);
 823     Matrix_Free(M);
 824
 825     if (union_next)
 826         D = DomainConcat(D, constraints2domain(union_next, nparam, MaxRays));
 827     return D;
 828 }
 829
 830 static evalue *evalue_read_partition(struct stream *s, struct parameter *p,
 831                                      char ***ppp,
 832                                      unsigned *nparam, unsigned MaxRays)
 833 {
 834     struct section *part = NULL;
 835     struct constraint *constraints;
 836     evalue *e = NULL;
 837     int m = 0;
 838
 839     while ((constraints = evalue_read_domain(s, &p, MaxRays))) {
 840         evalue *e = evalue_read_term(s, &p);
 841         if (!e) {
 842             stream_error(s, NULL, "missing evalue");
 843             break;
 844         }
 845         struct section *sect = ALLOC(struct section);
 846         sect->constraints = constraints;
 847         sect->e = e;
 848         sect->next = part;
 849         part = sect;
 850         ++m;
 851     }
 852
 853     if (part) {
 854         Polyhedron *D;
 855         int j;
 856
 857         *ppp = extract_parameters(p, nparam);
 858         e = ALLOC(evalue);
 859         value_init(e->d);
 860         e->x.p = new_enode(partition, 2*m, *nparam);
 861
 862         for (j = 0; j < m; ++j) {
 863             int n;
 864             struct section *next = part->next;
 865             constraints = part->constraints;
 866             D = constraints2domain(part->constraints, *nparam, MaxRays);
 867             EVALUE_SET_DOMAIN(e->x.p->arr[2*j], D);
 868             value_clear(e->x.p->arr[2*j+1].d);
 869             e->x.p->arr[2*j+1] = *part->e;
 870             free(part->e);
 871             free(part);
 872             part = next;
 873         }
 874     }
 875     return e;
 876 }
 877
 878 static evalue *evalue_read(FILE *in, char *var_list, char ***ppp, unsigned *nvar,
 879                            unsigned *nparam, unsigned MaxRays)
 880 {
 881     struct stream *s = stream_new(in);
 882     struct token *tok;
 883     evalue *e;
 884     struct parameter *p = NULL;
 885     char *next;
 886     int nv;
 887
 888     if (var_list) {
 889         while ((next = strchr(var_list, ','))) {
 890             *next = '\0';
 891             if (next > var_list)
 892                 parameter_pos(&p, var_list);
 893             *next = ',';
 894             var_list = next+1;
 895         }
 896         if (strlen(var_list) > 0)
 897             parameter_pos(&p, var_list);
 898         nv = p ? p->pos+1 : 0;
 899     } else
 900         nv = -1;
 901
 902     if (!(tok = stream_next_token(s)))
 903         return NULL;
 904
 905     if (tok->type == TOKEN_VARS) {
 906         token_free(tok);
 907         for (;;) {
 908             tok = stream_next_token(s);
 909             if (!tok || tok->type != TOKEN_IDENT) {
 910                 stream_error(s, tok, "expecting identifier");
 911                 break;
 912             }
 913             if (nv == -1)
 914                 parameter_pos(&p, tok->u.s);
 915             token_free(tok);
 916             tok = stream_next_token(s);
 917             if (!tok || tok->type != ',')
 918                 break;
 919             token_free(tok);
 920         }
 921         if (!tok)
 922             return NULL;
 923         if (nv = -1)
 924             nv = p ? p->pos+1 : 0;
 925     }
 926
 927     if (tok->type == TOKEN_VALUE) {
 928         struct token *tok2 = stream_next_token(s);
 929         if (tok2)
 930             stream_push_token(s, tok2);
 931         stream_push_token(s, tok);
 932         if (tok2 && (tok2->type == TOKEN_IDENT || tok2->type == TOKEN_GE))
 933             e = evalue_read_partition(s, p, ppp, nparam, MaxRays);
 934         else {
 935             e = evalue_read_term(s, &p);
 936             *ppp = extract_parameters(p, nparam);
 937         }
 938     } else if (tok->type == TOKEN_IDENT) {
 939         stream_push_token(s, tok);
 940         e = evalue_read_partition(s, p, ppp, nparam, MaxRays);
 941     } else {
 942         stream_error(s, tok, "unexpected token");
 943         *nparam = nv == -1 ? 0 : nv;
 944         e = NULL;
 945     }
 946     stream_free(s);
 947     if (nv == -1)
 948         *nvar = *nparam;
 949     else
 950         *nvar = nv;
 951     *nparam -= *nvar;
 952     return e;
 953 }
 954
 955 static int check_poly_max(const struct check_poly_data *data,
 956                           int nparam, Value *z,
 957                           const struct verify_options *options);
 958
 959 struct check_poly_max_data : public check_poly_data {
 960     Polyhedron          **S;
 961     evalue               *EP;
 962     piecewise_lst        *pl;
 963
 964     check_poly_max_data(Value *z, evalue *EP, piecewise_lst *pl) :
 965                         EP(EP), pl(pl) {
 966         this->z = z;
 967         this->check = check_poly_max;
 968     }
 969 };
 970
 971 static void optimum(Polyhedron *S, int pos, const check_poly_max_data *data,
 972                     Value *opt, bool& found,
 973                     const struct verify_options *options)
 974 {
 975     if (!S) {
 976         Value c;
 977         value_init(c);
 978         value_set_double(c, compute_evalue(data->EP, data->z+1)+.25);
 979         if (!found) {
 980             value_assign(*opt, c);
 981             found = true;
 982         } else {
 983             if (options->barvinok->bernstein_optimize == BV_BERNSTEIN_MAX) {
 984                 if (value_gt(c, *opt))
 985                     value_assign(*opt, c);
 986             } else {
 987                 if (value_lt(c, *opt))
 988                     value_assign(*opt, c);
 989             }
 990         }
 991         value_clear(c);
 992     } else {
 993         Value LB, UB;
 994         int ok;
 995         value_init(LB);
 996         value_init(UB);
 997         ok = !(lower_upper_bounds(1+pos, S, data->z, &LB, &UB));
 998         assert(ok);
 999         for (; value_le(LB, UB); value_increment(LB, LB)) {
1000             value_assign(data->z[1+pos], LB);
1001             optimum(S->next, pos+1, data, opt, found, options);
1002         }
1003         value_set_si(data->z[1+pos], 0);
1004         value_clear(LB);
1005         value_clear(UB);
1006     }
1007 }
1008
1009 static void optimum(const check_poly_max_data *data, Value *opt,
1010                     const struct verify_options *options)
1011 {
1012     bool found = false;
1013     for (int i = 0; i < data->EP->x.p->size/2; ++i)
1014         if (!emptyQ2(data->S[i]))
1015             optimum(data->S[i], 0, data, opt, found, options);
1016     assert(found);
1017 }
1018
1019 static int check_poly_max(const struct check_poly_data *data,
1020                           int nparam, Value *z,
1021                           const struct verify_options *options)
1022 {
1023     int k;
1024     int ok;
1025     const check_poly_max_data *max_data;
1026     max_data = static_cast<const check_poly_max_data *>(data);
1027     char *minmax;
1028     Value m, n, d;
1029     value_init(m);
1030     value_init(n);
1031     value_init(d);
1032
1033     if (options->barvinok->bernstein_optimize == BV_BERNSTEIN_MAX)
1034         minmax = "max";
1035     else
1036         minmax = "min";
1037
1038     max_data->pl->evaluate(nparam, z, &n, &d);
1039     if (options->barvinok->bernstein_optimize == BV_BERNSTEIN_MAX)
1040         mpz_fdiv_q(m, n, d);
1041     else
1042         mpz_cdiv_q(m, n, d);
1043
1044     if (options->print_all) {
1045         printf("%s(", minmax);
1046         value_print(stdout, VALUE_FMT, z[0]);
1047         for (k = 1; k < nparam; ++k) {
1048             printf(", ");
1049             value_print(stdout, VALUE_FMT, z[k]);
1050         }
1051         printf(") = ");
1052         value_print(stdout, VALUE_FMT, n);
1053         if (value_notone_p(d)) {
1054             printf("/");
1055             value_print(stdout, VALUE_FMT, d);
1056         }
1057         printf(" (");
1058         value_print(stdout, VALUE_FMT, m);
1059         printf(")");
1060     }
1061
1062     optimum(max_data, &n, options);
1063
1064     if (options->barvinok->bernstein_optimize == BV_BERNSTEIN_MAX)
1065         ok = value_ge(m, n);
1066     else
1067         ok = value_le(m, n);
1068
1069     if (options->print_all) {
1070         printf(", %s(EP) = ", minmax);
1071         value_print(stdout, VALUE_FMT, n);
1072         printf(". ");
1073     }
1074
1075     if (!ok) {
1076         printf("\n");
1077         fflush(stdout);
1078         fprintf(stderr, "Error !\n");
1079         fprintf(stderr, "%s(", minmax);
1080         value_print(stderr, VALUE_FMT, z[0]);
1081         for (k = 1; k < nparam; ++k) {
1082             fprintf(stderr,", ");
1083             value_print(stderr, VALUE_FMT, z[k]);
1084         }
1085         fprintf(stderr, ") should be ");
1086         if (options->barvinok->bernstein_optimize == BV_BERNSTEIN_MAX)
1087             fprintf(stderr, "greater");
1088         else
1089             fprintf(stderr, "smaller");
1090         fprintf(stderr, " than or equal to ");
1091         value_print(stderr, VALUE_FMT, n);
1092         fprintf(stderr, ", while pl eval gives ");
1093         value_print(stderr, VALUE_FMT, m);
1094         fprintf(stderr, ".\n");
1095         cerr << *max_data->pl << endl;
1096     } else if (options->print_all)
1097         printf("OK.\n");
1098
1099     value_clear(m);
1100     value_clear(n);
1101     value_clear(d);
1102
1103     return ok;
1104 }
1105
1106 static int verify(Polyhedron *D, piecewise_lst *pl, evalue *EP,
1107                   unsigned nvar, unsigned nparam, Vector *p,
1108                   struct verify_options *options)
1109 {
1110     Polyhedron *CS, *S;
1111     unsigned MaxRays = options->barvinok->MaxRays;
1112     assert(value_zero_p(EP->d));
1113     assert(EP->x.p->type == partition);
1114     int ok = 1;
1115
1116     CS = check_poly_context_scan(NULL, &D, D->Dimension, options);
1117
1118     check_poly_init(D, options);
1119
1120     if (!(CS && emptyQ2(CS))) {
1121         check_poly_max_data data(p->p, EP, pl);
1122         data.S = ALLOCN(Polyhedron *, EP->x.p->size/2);
1123         for (int i = 0; i < EP->x.p->size/2; ++i) {
1124             Polyhedron *A = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
1125             data.S[i] = Polyhedron_Scan(A, D, MaxRays & POL_NO_DUAL ? 0 : MaxRays);
1126         }
1127         ok = check_poly(CS, &data, nparam, 0, p->p+1+nvar, options);
1128         for (int i = 0; i < EP->x.p->size/2; ++i)
1129             Domain_Free(data.S[i]);
1130         free(data.S);
1131     }
1132
1133     if (!options->print_all)
1134         printf("\n");
1135
1136     if (CS) {
1137         Domain_Free(CS);
1138         Domain_Free(D);
1139     }
1140
1141     return ok;
1142 }
1143
1144 static int verify(piecewise_lst *pl, evalue *EP, unsigned nvar, unsigned nparam,
1145                   struct verify_options *options)
1146 {
1147     Vector *p;
1148
1149     p = Vector_Alloc(nvar+nparam+2);
1150     value_set_si(p->p[nvar+nparam+1], 1);
1151
1152     for (int i = 0; i < pl->list.size(); ++i) {
1153         int ok = verify(pl->list[i].first, pl, EP, nvar, nparam, p, options);
1154         if (!ok && !options->continue_on_error)
1155             break;
1156     }
1157
1158     Vector_Free(p);
1159
1160     return 0;
1161 }
1162
1163 static int optimize(evalue *EP, char **all_vars, unsigned nvar, unsigned nparam,
1164                     struct options *options)
1165 {
1166     Polyhedron *U;
1167     piecewise_lst *pl = NULL;
1168     U = Universe_Polyhedron(nparam);
1169     int print_solution = 1;
1170     int result = 0;
1171
1172     exvector params;
1173     params = constructParameterVector(all_vars+nvar, nparam);
1174
1175     if (options->verify.verify) {
1176         verify_options_set_range(&options->verify, nvar+nparam);
1177         if (!options->verbose)
1178             print_solution = 0;
1179     }
1180
1181     if (options->minimize)
1182         options->verify.barvinok->bernstein_optimize = BV_BERNSTEIN_MIN;
1183     else
1184         options->verify.barvinok->bernstein_optimize = BV_BERNSTEIN_MAX;
1185     pl = evalue_bernstein_coefficients(NULL, EP, U, params,
1186                                        options->verify.barvinok);
1187     assert(pl);
1188     if (options->minimize)
1189         pl->minimize();
1190     else
1191         pl->maximize();
1192     if (print_solution)
1193         cout << *pl << endl;
1194     if (options->verify.verify)
1195         result = verify(pl, EP, nvar, nparam, &options->verify);
1196     delete pl;
1197
1198     Polyhedron_Free(U);
1199
1200     return result;
1201 }
1202
1203 int main(int argc, char **argv)
1204 {
1205     evalue *EP;
1206     char **all_vars = NULL;
1207     unsigned nvar;
1208     unsigned nparam;
1209     struct options options;
1210     struct barvinok_options *bv_options = barvinok_options_new_with_defaults();
1211     static struct argp_child argp_children[] = {
1212         { &convert_argp,        0,      "input conversion",     1 },
1213         { &verify_argp,         0,      "verification",         2 },
1214         { &barvinok_argp,       0,      "barvinok options",     3 },
1215         { 0 }
1216     };
1217     static struct argp argp = { argp_options, parse_opt, 0, 0, argp_children };
1218     int result = 0;
1219
1220     options.verify.barvinok = bv_options;
1221     set_program_name(argv[0]);
1222     argp_parse(&argp, argc, argv, 0, 0, &options);
1223
1224     EP = evalue_read(stdin, options.var_list, &all_vars, &nvar, &nparam,
1225                      bv_options->MaxRays);
1226     assert(EP);
1227
1228     if (options.split)
1229         evalue_split_periods(EP, options.split, bv_options->MaxRays);
1230
1231     evalue_convert(EP, &options.convert, options.verbose, nparam, all_vars);
1232
1233     if (EVALUE_IS_ZERO(*EP))
1234         print_evalue(stdout, EP, all_vars);
1235     else
1236         result = optimize(EP, all_vars, nvar, nparam, &options);
1237
1238     free_evalue_refs(EP);
1239     free(EP);
1240
1241     if (options.var_list)
1242         free(options.var_list);
1243     Free_ParamNames(all_vars, nvar+nparam);
1244     barvinok_options_free(bv_options);
1245     return result;
1246 }