isl_tab_pip.c

   1 /*
   2  * Copyright 2008-2009 Katholieke Universiteit Leuven
   3  * Copyright 2010      INRIA Saclay
   4  *
   5  * Use of this software is governed by the MIT license
   6  *
   7  * Written by Sven Verdoolaege, K.U.Leuven, Departement
   8  * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
   9  * and INRIA Saclay - Ile-de-France, Parc Club Orsay Universite,
  10  * ZAC des vignes, 4 rue Jacques Monod, 91893 Orsay, France
  11  */
  12
  13 #include <isl_ctx_private.h>
  14 #include "isl_map_private.h"
  15 #include <isl/seq.h>
  16 #include "isl_tab.h"
  17 #include "isl_sample.h"
  18 #include <isl_mat_private.h>
  19 #include <isl_aff_private.h>
  20 #include <isl_options_private.h>
  21 #include <isl_config.h>
  22
  23 /*
  24  * The implementation of parametric integer linear programming in this file
  25  * was inspired by the paper "Parametric Integer Programming" and the
  26  * report "Solving systems of affine (in)equalities" by Paul Feautrier
  27  * (and others).
  28  *
  29  * The strategy used for obtaining a feasible solution is different
  30  * from the one used in isl_tab.c.  In particular, in isl_tab.c,
  31  * upon finding a constraint that is not yet satisfied, we pivot
  32  * in a row that increases the constant term of the row holding the
  33  * constraint, making sure the sample solution remains feasible
  34  * for all the constraints it already satisfied.
  35  * Here, we always pivot in the row holding the constraint,
  36  * choosing a column that induces the lexicographically smallest
  37  * increment to the sample solution.
  38  *
  39  * By starting out from a sample value that is lexicographically
  40  * smaller than any integer point in the problem space, the first
  41  * feasible integer sample point we find will also be the lexicographically
  42  * smallest.  If all variables can be assumed to be non-negative,
  43  * then the initial sample value may be chosen equal to zero.
  44  * However, we will not make this assumption.  Instead, we apply
  45  * the "big parameter" trick.  Any variable x is then not directly
  46  * used in the tableau, but instead it is represented by another
  47  * variable x' = M + x, where M is an arbitrarily large (positive)
  48  * value.  x' is therefore always non-negative, whatever the value of x.
  49  * Taking as initial sample value x' = 0 corresponds to x = -M,
  50  * which is always smaller than any possible value of x.
  51  *
  52  * The big parameter trick is used in the main tableau and
  53  * also in the context tableau if isl_context_lex is used.
  54  * In this case, each tableaus has its own big parameter.
  55  * Before doing any real work, we check if all the parameters
  56  * happen to be non-negative.  If so, we drop the column corresponding
  57  * to M from the initial context tableau.
  58  * If isl_context_gbr is used, then the big parameter trick is only
  59  * used in the main tableau.
  60  */
  61
  62 struct isl_context;
  63 struct isl_context_op {
  64         /* detect nonnegative parameters in context and mark them in tab */
  65         struct isl_tab *(*detect_nonnegative_parameters)(
  66                         struct isl_context *context, struct isl_tab *tab);
  67         /* return temporary reference to basic set representation of context */
  68         struct isl_basic_set *(*peek_basic_set)(struct isl_context *context);
  69         /* return temporary reference to tableau representation of context */
  70         struct isl_tab *(*peek_tab)(struct isl_context *context);
  71         /* add equality; check is 1 if eq may not be valid;
  72          * update is 1 if we may want to call ineq_sign on context later.
  73          */
  74         void (*add_eq)(struct isl_context *context, isl_int *eq,
  75                         int check, int update);
  76         /* add inequality; check is 1 if ineq may not be valid;
  77          * update is 1 if we may want to call ineq_sign on context later.
  78          */
  79         void (*add_ineq)(struct isl_context *context, isl_int *ineq,
  80                         int check, int update);
  81         /* check sign of ineq based on previous information.
  82          * strict is 1 if saturation should be treated as a positive sign.
  83          */
  84         enum isl_tab_row_sign (*ineq_sign)(struct isl_context *context,
  85                         isl_int *ineq, int strict);
  86         /* check if inequality maintains feasibility */
  87         int (*test_ineq)(struct isl_context *context, isl_int *ineq);
  88         /* return index of a div that corresponds to "div" */
  89         int (*get_div)(struct isl_context *context, struct isl_tab *tab,
  90                         struct isl_vec *div);
  91         /* add div "div" to context and return non-negativity */
  92         int (*add_div)(struct isl_context *context, struct isl_vec *div);
  93         int (*detect_equalities)(struct isl_context *context,
  94                         struct isl_tab *tab);
  95         /* return row index of "best" split */
  96         int (*best_split)(struct isl_context *context, struct isl_tab *tab);
  97         /* check if context has already been determined to be empty */
  98         int (*is_empty)(struct isl_context *context);
  99         /* check if context is still usable */
 100         int (*is_ok)(struct isl_context *context);
 101         /* save a copy/snapshot of context */
 102         void *(*save)(struct isl_context *context);
 103         /* restore saved context */
 104         void (*restore)(struct isl_context *context, void *);
 105         /* invalidate context */
 106         void (*invalidate)(struct isl_context *context);
 107         /* free context */
 108         void (*free)(struct isl_context *context);
 109 };
 110
 111 struct isl_context {
 112         struct isl_context_op *op;
 113 };
 114
 115 struct isl_context_lex {
 116         struct isl_context context;
 117         struct isl_tab *tab;
 118 };
 119
 120 struct isl_partial_sol {
 121         int level;
 122         struct isl_basic_set *dom;
 123         struct isl_mat *M;
 124
 125         struct isl_partial_sol *next;
 126 };
 127
 128 struct isl_sol;
 129 struct isl_sol_callback {
 130         struct isl_tab_callback callback;
 131         struct isl_sol *sol;
 132 };
 133
 134 /* isl_sol is an interface for constructing a solution to
 135  * a parametric integer linear programming problem.
 136  * Every time the algorithm reaches a state where a solution
 137  * can be read off from the tableau (including cases where the tableau
 138  * is empty), the function "add" is called on the isl_sol passed
 139  * to find_solutions_main.
 140  *
 141  * The context tableau is owned by isl_sol and is updated incrementally.
 142  *
 143  * There are currently two implementations of this interface,
 144  * isl_sol_map, which simply collects the solutions in an isl_map
 145  * and (optionally) the parts of the context where there is no solution
 146  * in an isl_set, and
 147  * isl_sol_for, which calls a user-defined function for each part of
 148  * the solution.
 149  */
 150 struct isl_sol {
 151         int error;
 152         int rational;
 153         int level;
 154         int max;
 155         int n_out;
 156         struct isl_context *context;
 157         struct isl_partial_sol *partial;
 158         void (*add)(struct isl_sol *sol,
 159                             struct isl_basic_set *dom, struct isl_mat *M);
 160         void (*add_empty)(struct isl_sol *sol, struct isl_basic_set *bset);
 161         void (*free)(struct isl_sol *sol);
 162         struct isl_sol_callback dec_level;
 163 };
 164
 165 static void sol_free(struct isl_sol *sol)
 166 {
 167         struct isl_partial_sol *partial, *next;
 168         if (!sol)
 169                 return;
 170         for (partial = sol->partial; partial; partial = next) {
 171                 next = partial->next;
 172                 isl_basic_set_free(partial->dom);
 173                 isl_mat_free(partial->M);
 174                 free(partial);
 175         }
 176         sol->free(sol);
 177 }
 178
 179 /* Push a partial solution represented by a domain and mapping M
 180  * onto the stack of partial solutions.
 181  */
 182 static void sol_push_sol(struct isl_sol *sol,
 183         struct isl_basic_set *dom, struct isl_mat *M)
 184 {
 185         struct isl_partial_sol *partial;
 186
 187         if (sol->error || !dom)
 188                 goto error;
 189
 190         partial = isl_alloc_type(dom->ctx, struct isl_partial_sol);
 191         if (!partial)
 192                 goto error;
 193
 194         partial->level = sol->level;
 195         partial->dom = dom;
 196         partial->M = M;
 197         partial->next = sol->partial;
 198
 199         sol->partial = partial;
 200
 201         return;
 202 error:
 203         isl_basic_set_free(dom);
 204         isl_mat_free(M);
 205         sol->error = 1;
 206 }
 207
 208 /* Pop one partial solution from the partial solution stack and
 209  * pass it on to sol->add or sol->add_empty.
 210  */
 211 static void sol_pop_one(struct isl_sol *sol)
 212 {
 213         struct isl_partial_sol *partial;
 214
 215         partial = sol->partial;
 216         sol->partial = partial->next;
 217
 218         if (partial->M)
 219                 sol->add(sol, partial->dom, partial->M);
 220         else
 221                 sol->add_empty(sol, partial->dom);
 222         free(partial);
 223 }
 224
 225 /* Return a fresh copy of the domain represented by the context tableau.
 226  */
 227 static struct isl_basic_set *sol_domain(struct isl_sol *sol)
 228 {
 229         struct isl_basic_set *bset;
 230
 231         if (sol->error)
 232                 return NULL;
 233
 234         bset = isl_basic_set_dup(sol->context->op->peek_basic_set(sol->context));
 235         bset = isl_basic_set_update_from_tab(bset,
 236                         sol->context->op->peek_tab(sol->context));
 237
 238         return bset;
 239 }
 240
 241 /* Check whether two partial solutions have the same mapping, where n_div
 242  * is the number of divs that the two partial solutions have in common.
 243  */
 244 static int same_solution(struct isl_partial_sol *s1, struct isl_partial_sol *s2,
 245         unsigned n_div)
 246 {
 247         int i;
 248         unsigned dim;
 249
 250         if (!s1->M != !s2->M)
 251                 return 0;
 252         if (!s1->M)
 253                 return 1;
 254
 255         dim = isl_basic_set_total_dim(s1->dom) - s1->dom->n_div;
 256
 257         for (i = 0; i < s1->M->n_row; ++i) {
 258                 if (isl_seq_first_non_zero(s1->M->row[i]+1+dim+n_div,
 259                                             s1->M->n_col-1-dim-n_div) != -1)
 260                         return 0;
 261                 if (isl_seq_first_non_zero(s2->M->row[i]+1+dim+n_div,
 262                                             s2->M->n_col-1-dim-n_div) != -1)
 263                         return 0;
 264                 if (!isl_seq_eq(s1->M->row[i], s2->M->row[i], 1+dim+n_div))
 265                         return 0;
 266         }
 267         return 1;
 268 }
 269
 270 /* Pop all solutions from the partial solution stack that were pushed onto
 271  * the stack at levels that are deeper than the current level.
 272  * If the two topmost elements on the stack have the same level
 273  * and represent the same solution, then their domains are combined.
 274  * This combined domain is the same as the current context domain
 275  * as sol_pop is called each time we move back to a higher level.
 276  */
 277 static void sol_pop(struct isl_sol *sol)
 278 {
 279         struct isl_partial_sol *partial;
 280         unsigned n_div;
 281
 282         if (sol->error)
 283                 return;
 284
 285         if (sol->level == 0) {
 286                 for (partial = sol->partial; partial; partial = sol->partial)
 287                         sol_pop_one(sol);
 288                 return;
 289         }
 290
 291         partial = sol->partial;
 292         if (!partial)
 293                 return;
 294
 295         if (partial->level <= sol->level)
 296                 return;
 297
 298         if (partial->next && partial->next->level == partial->level) {
 299                 n_div = isl_basic_set_dim(
 300                                 sol->context->op->peek_basic_set(sol->context),
 301                                 isl_dim_div);
 302
 303                 if (!same_solution(partial, partial->next, n_div)) {
 304                         sol_pop_one(sol);
 305                         sol_pop_one(sol);
 306                 } else {
 307                         struct isl_basic_set *bset;
 308
 309                         bset = sol_domain(sol);
 310
 311                         isl_basic_set_free(partial->next->dom);
 312                         partial->next->dom = bset;
 313                         partial->next->level = sol->level;
 314
 315                         sol->partial = partial->next;
 316                         isl_basic_set_free(partial->dom);
 317                         isl_mat_free(partial->M);
 318                         free(partial);
 319                 }
 320         } else
 321                 sol_pop_one(sol);
 322 }
 323
 324 static void sol_dec_level(struct isl_sol *sol)
 325 {
 326         if (sol->error)
 327                 return;
 328
 329         sol->level--;
 330
 331         sol_pop(sol);
 332 }
 333
 334 static int sol_dec_level_wrap(struct isl_tab_callback *cb)
 335 {
 336         struct isl_sol_callback *callback = (struct isl_sol_callback *)cb;
 337
 338         sol_dec_level(callback->sol);
 339
 340         return callback->sol->error ? -1 : 0;
 341 }
 342
 343 /* Move down to next level and push callback onto context tableau
 344  * to decrease the level again when it gets rolled back across
 345  * the current state.  That is, dec_level will be called with
 346  * the context tableau in the same state as it is when inc_level
 347  * is called.
 348  */
 349 static void sol_inc_level(struct isl_sol *sol)
 350 {
 351         struct isl_tab *tab;
 352
 353         if (sol->error)
 354                 return;
 355
 356         sol->level++;
 357         tab = sol->context->op->peek_tab(sol->context);
 358         if (isl_tab_push_callback(tab, &sol->dec_level.callback) < 0)
 359                 sol->error = 1;
 360 }
 361
 362 static void scale_rows(struct isl_mat *mat, isl_int m, int n_row)
 363 {
 364         int i;
 365
 366         if (isl_int_is_one(m))
 367                 return;
 368
 369         for (i = 0; i < n_row; ++i)
 370                 isl_seq_scale(mat->row[i], mat->row[i], m, mat->n_col);
 371 }
 372
 373 /* Add the solution identified by the tableau and the context tableau.
 374  *
 375  * The layout of the variables is as follows.
 376  *      tab->n_var is equal to the total number of variables in the input
 377  *                      map (including divs that were copied from the context)
 378  *                      + the number of extra divs constructed
 379  *      Of these, the first tab->n_param and the last tab->n_div variables
 380  *      correspond to the variables in the context, i.e.,
 381  *              tab->n_param + tab->n_div = context_tab->n_var
 382  *      tab->n_param is equal to the number of parameters and input
 383  *                      dimensions in the input map
 384  *      tab->n_div is equal to the number of divs in the context
 385  *
 386  * If there is no solution, then call add_empty with a basic set
 387  * that corresponds to the context tableau.  (If add_empty is NULL,
 388  * then do nothing).
 389  *
 390  * If there is a solution, then first construct a matrix that maps
 391  * all dimensions of the context to the output variables, i.e.,
 392  * the output dimensions in the input map.
 393  * The divs in the input map (if any) that do not correspond to any
 394  * div in the context do not appear in the solution.
 395  * The algorithm will make sure that they have an integer value,
 396  * but these values themselves are of no interest.
 397  * We have to be careful not to drop or rearrange any divs in the
 398  * context because that would change the meaning of the matrix.
 399  *
 400  * To extract the value of the output variables, it should be noted
 401  * that we always use a big parameter M in the main tableau and so
 402  * the variable stored in this tableau is not an output variable x itself, but
 403  *      x' = M + x (in case of minimization)
 404  * or
 405  *      x' = M - x (in case of maximization)
 406  * If x' appears in a column, then its optimal value is zero,
 407  * which means that the optimal value of x is an unbounded number
 408  * (-M for minimization and M for maximization).
 409  * We currently assume that the output dimensions in the original map
 410  * are bounded, so this cannot occur.
 411  * Similarly, when x' appears in a row, then the coefficient of M in that
 412  * row is necessarily 1.
 413  * If the row in the tableau represents
 414  *      d x' = c + d M + e(y)
 415  * then, in case of minimization, the corresponding row in the matrix
 416  * will be
 417  *      a c + a e(y)
 418  * with a d = m, the (updated) common denominator of the matrix.
 419  * In case of maximization, the row will be
 420  *      -a c - a e(y)
 421  */
 422 static void sol_add(struct isl_sol *sol, struct isl_tab *tab)
 423 {
 424         struct isl_basic_set *bset = NULL;
 425         struct isl_mat *mat = NULL;
 426         unsigned off;
 427         int row;
 428         isl_int m;
 429
 430         if (sol->error || !tab)
 431                 goto error;
 432
 433         if (tab->empty && !sol->add_empty)
 434                 return;
 435         if (sol->context->op->is_empty(sol->context))
 436                 return;
 437
 438         bset = sol_domain(sol);
 439
 440         if (tab->empty) {
 441                 sol_push_sol(sol, bset, NULL);
 442                 return;
 443         }
 444
 445         off = 2 + tab->M;
 446
 447         mat = isl_mat_alloc(tab->mat->ctx, 1 + sol->n_out,
 448                                             1 + tab->n_param + tab->n_div);
 449         if (!mat)
 450                 goto error;
 451
 452         isl_int_init(m);
 453
 454         isl_seq_clr(mat->row[0] + 1, mat->n_col - 1);
 455         isl_int_set_si(mat->row[0][0], 1);
 456         for (row = 0; row < sol->n_out; ++row) {
 457                 int i = tab->n_param + row;
 458                 int r, j;
 459
 460                 isl_seq_clr(mat->row[1 + row], mat->n_col);
 461                 if (!tab->var[i].is_row) {
 462                         if (tab->M)
 463                                 isl_die(mat->ctx, isl_error_invalid,
 464                                         "unbounded optimum", goto error2);
 465                         continue;
 466                 }
 467
 468                 r = tab->var[i].index;
 469                 if (tab->M &&
 470                     isl_int_ne(tab->mat->row[r][2], tab->mat->row[r][0]))
 471                         isl_die(mat->ctx, isl_error_invalid,
 472                                 "unbounded optimum", goto error2);
 473                 isl_int_gcd(m, mat->row[0][0], tab->mat->row[r][0]);
 474                 isl_int_divexact(m, tab->mat->row[r][0], m);
 475                 scale_rows(mat, m, 1 + row);
 476                 isl_int_divexact(m, mat->row[0][0], tab->mat->row[r][0]);
 477                 isl_int_mul(mat->row[1 + row][0], m, tab->mat->row[r][1]);
 478                 for (j = 0; j < tab->n_param; ++j) {
 479                         int col;
 480                         if (tab->var[j].is_row)
 481                                 continue;
 482                         col = tab->var[j].index;
 483                         isl_int_mul(mat->row[1 + row][1 + j], m,
 484                                     tab->mat->row[r][off + col]);
 485                 }
 486                 for (j = 0; j < tab->n_div; ++j) {
 487                         int col;
 488                         if (tab->var[tab->n_var - tab->n_div+j].is_row)
 489                                 continue;
 490                         col = tab->var[tab->n_var - tab->n_div+j].index;
 491                         isl_int_mul(mat->row[1 + row][1 + tab->n_param + j], m,
 492                                     tab->mat->row[r][off + col]);
 493                 }
 494                 if (sol->max)
 495                         isl_seq_neg(mat->row[1 + row], mat->row[1 + row],
 496                                     mat->n_col);
 497         }
 498
 499         isl_int_clear(m);
 500
 501         sol_push_sol(sol, bset, mat);
 502         return;
 503 error2:
 504         isl_int_clear(m);
 505 error:
 506         isl_basic_set_free(bset);
 507         isl_mat_free(mat);
 508         sol->error = 1;
 509 }
 510
 511 struct isl_sol_map {
 512         struct isl_sol  sol;
 513         struct isl_map  *map;
 514         struct isl_set  *empty;
 515 };
 516
 517 static void sol_map_free(struct isl_sol_map *sol_map)
 518 {
 519         if (!sol_map)
 520                 return;
 521         if (sol_map->sol.context)
 522                 sol_map->sol.context->op->free(sol_map->sol.context);
 523         isl_map_free(sol_map->map);
 524         isl_set_free(sol_map->empty);
 525         free(sol_map);
 526 }
 527
 528 static void sol_map_free_wrap(struct isl_sol *sol)
 529 {
 530         sol_map_free((struct isl_sol_map *)sol);
 531 }
 532
 533 /* This function is called for parts of the context where there is
 534  * no solution, with "bset" corresponding to the context tableau.
 535  * Simply add the basic set to the set "empty".
 536  */
 537 static void sol_map_add_empty(struct isl_sol_map *sol,
 538         struct isl_basic_set *bset)
 539 {
 540         if (!bset)
 541                 goto error;
 542         isl_assert(bset->ctx, sol->empty, goto error);
 543
 544         sol->empty = isl_set_grow(sol->empty, 1);
 545         bset = isl_basic_set_simplify(bset);
 546         bset = isl_basic_set_finalize(bset);
 547         sol->empty = isl_set_add_basic_set(sol->empty, isl_basic_set_copy(bset));
 548         if (!sol->empty)
 549                 goto error;
 550         isl_basic_set_free(bset);
 551         return;
 552 error:
 553         isl_basic_set_free(bset);
 554         sol->sol.error = 1;
 555 }
 556
 557 static void sol_map_add_empty_wrap(struct isl_sol *sol,
 558         struct isl_basic_set *bset)
 559 {
 560         sol_map_add_empty((struct isl_sol_map *)sol, bset);
 561 }
 562
 563 /* Given a basic map "dom" that represents the context and an affine
 564  * matrix "M" that maps the dimensions of the context to the
 565  * output variables, construct a basic map with the same parameters
 566  * and divs as the context, the dimensions of the context as input
 567  * dimensions and a number of output dimensions that is equal to
 568  * the number of output dimensions in the input map.
 569  *
 570  * The constraints and divs of the context are simply copied
 571  * from "dom".  For each row
 572  *      x = c + e(y)
 573  * an equality
 574  *      c + e(y) - d x = 0
 575  * is added, with d the common denominator of M.
 576  */
 577 static void sol_map_add(struct isl_sol_map *sol,
 578         struct isl_basic_set *dom, struct isl_mat *M)
 579 {
 580         int i;
 581         struct isl_basic_map *bmap = NULL;
 582         unsigned n_eq;
 583         unsigned n_ineq;
 584         unsigned nparam;
 585         unsigned total;
 586         unsigned n_div;
 587         unsigned n_out;
 588
 589         if (sol->sol.error || !dom || !M)
 590                 goto error;
 591
 592         n_out = sol->sol.n_out;
 593         n_eq = dom->n_eq + n_out;
 594         n_ineq = dom->n_ineq;
 595         n_div = dom->n_div;
 596         nparam = isl_basic_set_total_dim(dom) - n_div;
 597         total = isl_map_dim(sol->map, isl_dim_all);
 598         bmap = isl_basic_map_alloc_space(isl_map_get_space(sol->map),
 599                                         n_div, n_eq, 2 * n_div + n_ineq);
 600         if (!bmap)
 601                 goto error;
 602         if (sol->sol.rational)
 603                 ISL_F_SET(bmap, ISL_BASIC_MAP_RATIONAL);
 604         for (i = 0; i < dom->n_div; ++i) {
 605                 int k = isl_basic_map_alloc_div(bmap);
 606                 if (k < 0)
 607                         goto error;
 608                 isl_seq_cpy(bmap->div[k], dom->div[i], 1 + 1 + nparam);
 609                 isl_seq_clr(bmap->div[k] + 1 + 1 + nparam, total - nparam);
 610                 isl_seq_cpy(bmap->div[k] + 1 + 1 + total,
 611                             dom->div[i] + 1 + 1 + nparam, i);
 612         }
 613         for (i = 0; i < dom->n_eq; ++i) {
 614                 int k = isl_basic_map_alloc_equality(bmap);
 615                 if (k < 0)
 616                         goto error;
 617                 isl_seq_cpy(bmap->eq[k], dom->eq[i], 1 + nparam);
 618                 isl_seq_clr(bmap->eq[k] + 1 + nparam, total - nparam);
 619                 isl_seq_cpy(bmap->eq[k] + 1 + total,
 620                             dom->eq[i] + 1 + nparam, n_div);
 621         }
 622         for (i = 0; i < dom->n_ineq; ++i) {
 623                 int k = isl_basic_map_alloc_inequality(bmap);
 624                 if (k < 0)
 625                         goto error;
 626                 isl_seq_cpy(bmap->ineq[k], dom->ineq[i], 1 + nparam);
 627                 isl_seq_clr(bmap->ineq[k] + 1 + nparam, total - nparam);
 628                 isl_seq_cpy(bmap->ineq[k] + 1 + total,
 629                         dom->ineq[i] + 1 + nparam, n_div);
 630         }
 631         for (i = 0; i < M->n_row - 1; ++i) {
 632                 int k = isl_basic_map_alloc_equality(bmap);
 633                 if (k < 0)
 634                         goto error;
 635                 isl_seq_cpy(bmap->eq[k], M->row[1 + i], 1 + nparam);
 636                 isl_seq_clr(bmap->eq[k] + 1 + nparam, n_out);
 637                 isl_int_neg(bmap->eq[k][1 + nparam + i], M->row[0][0]);
 638                 isl_seq_cpy(bmap->eq[k] + 1 + nparam + n_out,
 639                             M->row[1 + i] + 1 + nparam, n_div);
 640         }
 641         bmap = isl_basic_map_simplify(bmap);
 642         bmap = isl_basic_map_finalize(bmap);
 643         sol->map = isl_map_grow(sol->map, 1);
 644         sol->map = isl_map_add_basic_map(sol->map, bmap);
 645         isl_basic_set_free(dom);
 646         isl_mat_free(M);
 647         if (!sol->map)
 648                 sol->sol.error = 1;
 649         return;
 650 error:
 651         isl_basic_set_free(dom);
 652         isl_mat_free(M);
 653         isl_basic_map_free(bmap);
 654         sol->sol.error = 1;
 655 }
 656
 657 static void sol_map_add_wrap(struct isl_sol *sol,
 658         struct isl_basic_set *dom, struct isl_mat *M)
 659 {
 660         sol_map_add((struct isl_sol_map *)sol, dom, M);
 661 }
 662
 663
 664 /* Store the "parametric constant" of row "row" of tableau "tab" in "line",
 665  * i.e., the constant term and the coefficients of all variables that
 666  * appear in the context tableau.
 667  * Note that the coefficient of the big parameter M is NOT copied.
 668  * The context tableau may not have a big parameter and even when it
 669  * does, it is a different big parameter.
 670  */
 671 static void get_row_parameter_line(struct isl_tab *tab, int row, isl_int *line)
 672 {
 673         int i;
 674         unsigned off = 2 + tab->M;
 675
 676         isl_int_set(line[0], tab->mat->row[row][1]);
 677         for (i = 0; i < tab->n_param; ++i) {
 678                 if (tab->var[i].is_row)
 679                         isl_int_set_si(line[1 + i], 0);
 680                 else {
 681                         int col = tab->var[i].index;
 682                         isl_int_set(line[1 + i], tab->mat->row[row][off + col]);
 683                 }
 684         }
 685         for (i = 0; i < tab->n_div; ++i) {
 686                 if (tab->var[tab->n_var - tab->n_div + i].is_row)
 687                         isl_int_set_si(line[1 + tab->n_param + i], 0);
 688                 else {
 689                         int col = tab->var[tab->n_var - tab->n_div + i].index;
 690                         isl_int_set(line[1 + tab->n_param + i],
 691                                     tab->mat->row[row][off + col]);
 692                 }
 693         }
 694 }
 695
 696 /* Check if rows "row1" and "row2" have identical "parametric constants",
 697  * as explained above.
 698  * In this case, we also insist that the coefficients of the big parameter
 699  * be the same as the values of the constants will only be the same
 700  * if these coefficients are also the same.
 701  */
 702 static int identical_parameter_line(struct isl_tab *tab, int row1, int row2)
 703 {
 704         int i;
 705         unsigned off = 2 + tab->M;
 706
 707         if (isl_int_ne(tab->mat->row[row1][1], tab->mat->row[row2][1]))
 708                 return 0;
 709
 710         if (tab->M && isl_int_ne(tab->mat->row[row1][2],
 711                                  tab->mat->row[row2][2]))
 712                 return 0;
 713
 714         for (i = 0; i < tab->n_param + tab->n_div; ++i) {
 715                 int pos = i < tab->n_param ? i :
 716                         tab->n_var - tab->n_div + i - tab->n_param;
 717                 int col;
 718
 719                 if (tab->var[pos].is_row)
 720                         continue;
 721                 col = tab->var[pos].index;
 722                 if (isl_int_ne(tab->mat->row[row1][off + col],
 723                                tab->mat->row[row2][off + col]))
 724                         return 0;
 725         }
 726         return 1;
 727 }
 728
 729 /* Return an inequality that expresses that the "parametric constant"
 730  * should be non-negative.
 731  * This function is only called when the coefficient of the big parameter
 732  * is equal to zero.
 733  */
 734 static struct isl_vec *get_row_parameter_ineq(struct isl_tab *tab, int row)
 735 {
 736         struct isl_vec *ineq;
 737
 738         ineq = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_param + tab->n_div);
 739         if (!ineq)
 740                 return NULL;
 741
 742         get_row_parameter_line(tab, row, ineq->el);
 743         if (ineq)
 744                 ineq = isl_vec_normalize(ineq);
 745
 746         return ineq;
 747 }
 748
 749 /* Normalize a div expression of the form
 750  *
 751  *      [(g*f(x) + c)/(g * m)]
 752  *
 753  * with c the constant term and f(x) the remaining coefficients, to
 754  *
 755  *      [(f(x) + [c/g])/m]
 756  */
 757 static void normalize_div(__isl_keep isl_vec *div)
 758 {
 759         isl_ctx *ctx = isl_vec_get_ctx(div);
 760         int len = div->size - 2;
 761
 762         isl_seq_gcd(div->el + 2, len, &ctx->normalize_gcd);
 763         isl_int_gcd(ctx->normalize_gcd, ctx->normalize_gcd, div->el[0]);
 764
 765         if (isl_int_is_one(ctx->normalize_gcd))
 766                 return;
 767
 768         isl_int_divexact(div->el[0], div->el[0], ctx->normalize_gcd);
 769         isl_int_fdiv_q(div->el[1], div->el[1], ctx->normalize_gcd);
 770         isl_seq_scale_down(div->el + 2, div->el + 2, ctx->normalize_gcd, len);
 771 }
 772
 773 /* Return a integer division for use in a parametric cut based on the given row.
 774  * In particular, let the parametric constant of the row be
 775  *
 776  *              \sum_i a_i y_i
 777  *
 778  * where y_0 = 1, but none of the y_i corresponds to the big parameter M.
 779  * The div returned is equal to
 780  *
 781  *              floor(\sum_i {-a_i} y_i) = floor((\sum_i (-a_i mod d) y_i)/d)
 782  */
 783 static struct isl_vec *get_row_parameter_div(struct isl_tab *tab, int row)
 784 {
 785         struct isl_vec *div;
 786
 787         div = isl_vec_alloc(tab->mat->ctx, 1 + 1 + tab->n_param + tab->n_div);
 788         if (!div)
 789                 return NULL;
 790
 791         isl_int_set(div->el[0], tab->mat->row[row][0]);
 792         get_row_parameter_line(tab, row, div->el + 1);
 793         isl_seq_neg(div->el + 1, div->el + 1, div->size - 1);
 794         normalize_div(div);
 795         isl_seq_fdiv_r(div->el + 1, div->el + 1, div->el[0], div->size - 1);
 796
 797         return div;
 798 }
 799
 800 /* Return a integer division for use in transferring an integrality constraint
 801  * to the context.
 802  * In particular, let the parametric constant of the row be
 803  *
 804  *              \sum_i a_i y_i
 805  *
 806  * where y_0 = 1, but none of the y_i corresponds to the big parameter M.
 807  * The the returned div is equal to
 808  *
 809  *              floor(\sum_i {a_i} y_i) = floor((\sum_i (a_i mod d) y_i)/d)
 810  */
 811 static struct isl_vec *get_row_split_div(struct isl_tab *tab, int row)
 812 {
 813         struct isl_vec *div;
 814
 815         div = isl_vec_alloc(tab->mat->ctx, 1 + 1 + tab->n_param + tab->n_div);
 816         if (!div)
 817                 return NULL;
 818
 819         isl_int_set(div->el[0], tab->mat->row[row][0]);
 820         get_row_parameter_line(tab, row, div->el + 1);
 821         normalize_div(div);
 822         isl_seq_fdiv_r(div->el + 1, div->el + 1, div->el[0], div->size - 1);
 823
 824         return div;
 825 }
 826
 827 /* Construct and return an inequality that expresses an upper bound
 828  * on the given div.
 829  * In particular, if the div is given by
 830  *
 831  *      d = floor(e/m)
 832  *
 833  * then the inequality expresses
 834  *
 835  *      m d <= e
 836  */
 837 static struct isl_vec *ineq_for_div(struct isl_basic_set *bset, unsigned div)
 838 {
 839         unsigned total;
 840         unsigned div_pos;
 841         struct isl_vec *ineq;
 842
 843         if (!bset)
 844                 return NULL;
 845
 846         total = isl_basic_set_total_dim(bset);
 847         div_pos = 1 + total - bset->n_div + div;
 848
 849         ineq = isl_vec_alloc(bset->ctx, 1 + total);
 850         if (!ineq)
 851                 return NULL;
 852
 853         isl_seq_cpy(ineq->el, bset->div[div] + 1, 1 + total);
 854         isl_int_neg(ineq->el[div_pos], bset->div[div][0]);
 855         return ineq;
 856 }
 857
 858 /* Given a row in the tableau and a div that was created
 859  * using get_row_split_div and that has been constrained to equality, i.e.,
 860  *
 861  *              d = floor(\sum_i {a_i} y_i) = \sum_i {a_i} y_i
 862  *
 863  * replace the expression "\sum_i {a_i} y_i" in the row by d,
 864  * i.e., we subtract "\sum_i {a_i} y_i" and add 1 d.
 865  * The coefficients of the non-parameters in the tableau have been
 866  * verified to be integral.  We can therefore simply replace coefficient b
 867  * by floor(b).  For the coefficients of the parameters we have
 868  * floor(a_i) = a_i - {a_i}, while for the other coefficients, we have
 869  * floor(b) = b.
 870  */
 871 static struct isl_tab *set_row_cst_to_div(struct isl_tab *tab, int row, int div)
 872 {
 873         isl_seq_fdiv_q(tab->mat->row[row] + 1, tab->mat->row[row] + 1,
 874                         tab->mat->row[row][0], 1 + tab->M + tab->n_col);
 875
 876         isl_int_set_si(tab->mat->row[row][0], 1);
 877
 878         if (tab->var[tab->n_var - tab->n_div + div].is_row) {
 879                 int drow = tab->var[tab->n_var - tab->n_div + div].index;
 880
 881                 isl_assert(tab->mat->ctx,
 882                         isl_int_is_one(tab->mat->row[drow][0]), goto error);
 883                 isl_seq_combine(tab->mat->row[row] + 1,
 884                         tab->mat->ctx->one, tab->mat->row[row] + 1,
 885                         tab->mat->ctx->one, tab->mat->row[drow] + 1,
 886                         1 + tab->M + tab->n_col);
 887         } else {
 888                 int dcol = tab->var[tab->n_var - tab->n_div + div].index;
 889
 890                 isl_int_add_ui(tab->mat->row[row][2 + tab->M + dcol],
 891                                 tab->mat->row[row][2 + tab->M + dcol], 1);
 892         }
 893
 894         return tab;
 895 error:
 896         isl_tab_free(tab);
 897         return NULL;
 898 }
 899
 900 /* Check if the (parametric) constant of the given row is obviously
 901  * negative, meaning that we don't need to consult the context tableau.
 902  * If there is a big parameter and its coefficient is non-zero,
 903  * then this coefficient determines the outcome.
 904  * Otherwise, we check whether the constant is negative and
 905  * all non-zero coefficients of parameters are negative and
 906  * belong to non-negative parameters.
 907  */
 908 static int is_obviously_neg(struct isl_tab *tab, int row)
 909 {
 910         int i;
 911         int col;
 912         unsigned off = 2 + tab->M;
 913
 914         if (tab->M) {
 915                 if (isl_int_is_pos(tab->mat->row[row][2]))
 916                         return 0;
 917                 if (isl_int_is_neg(tab->mat->row[row][2]))
 918                         return 1;
 919         }
 920
 921         if (isl_int_is_nonneg(tab->mat->row[row][1]))
 922                 return 0;
 923         for (i = 0; i < tab->n_param; ++i) {
 924                 /* Eliminated parameter */
 925                 if (tab->var[i].is_row)
 926                         continue;
 927                 col = tab->var[i].index;
 928                 if (isl_int_is_zero(tab->mat->row[row][off + col]))
 929                         continue;
 930                 if (!tab->var[i].is_nonneg)
 931                         return 0;
 932                 if (isl_int_is_pos(tab->mat->row[row][off + col]))
 933                         return 0;
 934         }
 935         for (i = 0; i < tab->n_div; ++i) {
 936                 if (tab->var[tab->n_var - tab->n_div + i].is_row)
 937                         continue;
 938                 col = tab->var[tab->n_var - tab->n_div + i].index;
 939                 if (isl_int_is_zero(tab->mat->row[row][off + col]))
 940                         continue;
 941                 if (!tab->var[tab->n_var - tab->n_div + i].is_nonneg)
 942                         return 0;
 943                 if (isl_int_is_pos(tab->mat->row[row][off + col]))
 944                         return 0;
 945         }
 946         return 1;
 947 }
 948
 949 /* Check if the (parametric) constant of the given row is obviously
 950  * non-negative, meaning that we don't need to consult the context tableau.
 951  * If there is a big parameter and its coefficient is non-zero,
 952  * then this coefficient determines the outcome.
 953  * Otherwise, we check whether the constant is non-negative and
 954  * all non-zero coefficients of parameters are positive and
 955  * belong to non-negative parameters.
 956  */
 957 static int is_obviously_nonneg(struct isl_tab *tab, int row)
 958 {
 959         int i;
 960         int col;
 961         unsigned off = 2 + tab->M;
 962
 963         if (tab->M) {
 964                 if (isl_int_is_pos(tab->mat->row[row][2]))
 965                         return 1;
 966                 if (isl_int_is_neg(tab->mat->row[row][2]))
 967                         return 0;
 968         }
 969
 970         if (isl_int_is_neg(tab->mat->row[row][1]))
 971                 return 0;
 972         for (i = 0; i < tab->n_param; ++i) {
 973                 /* Eliminated parameter */
 974                 if (tab->var[i].is_row)
 975                         continue;
 976                 col = tab->var[i].index;
 977                 if (isl_int_is_zero(tab->mat->row[row][off + col]))
 978                         continue;
 979                 if (!tab->var[i].is_nonneg)
 980                         return 0;
 981                 if (isl_int_is_neg(tab->mat->row[row][off + col]))
 982                         return 0;
 983         }
 984         for (i = 0; i < tab->n_div; ++i) {
 985                 if (tab->var[tab->n_var - tab->n_div + i].is_row)
 986                         continue;
 987                 col = tab->var[tab->n_var - tab->n_div + i].index;
 988                 if (isl_int_is_zero(tab->mat->row[row][off + col]))
 989                         continue;
 990                 if (!tab->var[tab->n_var - tab->n_div + i].is_nonneg)
 991                         return 0;
 992                 if (isl_int_is_neg(tab->mat->row[row][off + col]))
 993                         return 0;
 994         }
 995         return 1;
 996 }
 997
 998 /* Given a row r and two columns, return the column that would
 999  * lead to the lexicographically smallest increment in the sample
1000  * solution when leaving the basis in favor of the row.
1001  * Pivoting with column c will increment the sample value by a non-negative
1002  * constant times a_{V,c}/a_{r,c}, with a_{V,c} the elements of column c
1003  * corresponding to the non-parametric variables.
1004  * If variable v appears in a column c_v, the a_{v,c} = 1 iff c = c_v,
1005  * with all other entries in this virtual row equal to zero.
1006  * If variable v appears in a row, then a_{v,c} is the element in column c
1007  * of that row.
1008  *
1009  * Let v be the first variable with a_{v,c1}/a_{r,c1} != a_{v,c2}/a_{r,c2}.
1010  * Then if a_{v,c1}/a_{r,c1} < a_{v,c2}/a_{r,c2}, i.e.,
1011  * a_{v,c2} a_{r,c1} - a_{v,c1} a_{r,c2} > 0, c1 results in the minimal
1012  * increment.  Otherwise, it's c2.
1013  */
1014 static int lexmin_col_pair(struct isl_tab *tab,
1015         int row, int col1, int col2, isl_int tmp)
1016 {
1017         int i;
1018         isl_int *tr;
1019
1020         tr = tab->mat->row[row] + 2 + tab->M;
1021
1022         for (i = tab->n_param; i < tab->n_var - tab->n_div; ++i) {
1023                 int s1, s2;
1024                 isl_int *r;
1025
1026                 if (!tab->var[i].is_row) {
1027                         if (tab->var[i].index == col1)
1028                                 return col2;
1029                         if (tab->var[i].index == col2)
1030                                 return col1;
1031                         continue;
1032                 }
1033
1034                 if (tab->var[i].index == row)
1035                         continue;
1036
1037                 r = tab->mat->row[tab->var[i].index] + 2 + tab->M;
1038                 s1 = isl_int_sgn(r[col1]);
1039                 s2 = isl_int_sgn(r[col2]);
1040                 if (s1 == 0 && s2 == 0)
1041                         continue;
1042                 if (s1 < s2)
1043                         return col1;
1044                 if (s2 < s1)
1045                         return col2;
1046
1047                 isl_int_mul(tmp, r[col2], tr[col1]);
1048                 isl_int_submul(tmp, r[col1], tr[col2]);
1049                 if (isl_int_is_pos(tmp))
1050                         return col1;
1051                 if (isl_int_is_neg(tmp))
1052                         return col2;
1053         }
1054         return -1;
1055 }
1056
1057 /* Given a row in the tableau, find and return the column that would
1058  * result in the lexicographically smallest, but positive, increment
1059  * in the sample point.
1060  * If there is no such column, then return tab->n_col.
1061  * If anything goes wrong, return -1.
1062  */
1063 static int lexmin_pivot_col(struct isl_tab *tab, int row)
1064 {
1065         int j;
1066         int col = tab->n_col;
1067         isl_int *tr;
1068         isl_int tmp;
1069
1070         tr = tab->mat->row[row] + 2 + tab->M;
1071
1072         isl_int_init(tmp);
1073
1074         for (j = tab->n_dead; j < tab->n_col; ++j) {
1075                 if (tab->col_var[j] >= 0 &&
1076                     (tab->col_var[j] < tab->n_param  ||
1077                     tab->col_var[j] >= tab->n_var - tab->n_div))
1078                         continue;
1079
1080                 if (!isl_int_is_pos(tr[j]))
1081                         continue;
1082
1083                 if (col == tab->n_col)
1084                         col = j;
1085                 else
1086                         col = lexmin_col_pair(tab, row, col, j, tmp);
1087                 isl_assert(tab->mat->ctx, col >= 0, goto error);
1088         }
1089
1090         isl_int_clear(tmp);
1091         return col;
1092 error:
1093         isl_int_clear(tmp);
1094         return -1;
1095 }
1096
1097 /* Return the first known violated constraint, i.e., a non-negative
1098  * constraint that currently has an either obviously negative value
1099  * or a previously determined to be negative value.
1100  *
1101  * If any constraint has a negative coefficient for the big parameter,
1102  * if any, then we return one of these first.
1103  */
1104 static int first_neg(struct isl_tab *tab)
1105 {
1106         int row;
1107
1108         if (tab->M)
1109                 for (row = tab->n_redundant; row < tab->n_row; ++row) {
1110                         if (!isl_tab_var_from_row(tab, row)->is_nonneg)
1111                                 continue;
1112                         if (!isl_int_is_neg(tab->mat->row[row][2]))
1113                                 continue;
1114                         if (tab->row_sign)
1115                                 tab->row_sign[row] = isl_tab_row_neg;
1116                         return row;
1117                 }
1118         for (row = tab->n_redundant; row < tab->n_row; ++row) {
1119                 if (!isl_tab_var_from_row(tab, row)->is_nonneg)
1120                         continue;
1121                 if (tab->row_sign) {
1122                         if (tab->row_sign[row] == 0 &&
1123                             is_obviously_neg(tab, row))
1124                                 tab->row_sign[row] = isl_tab_row_neg;
1125                         if (tab->row_sign[row] != isl_tab_row_neg)
1126                                 continue;
1127                 } else if (!is_obviously_neg(tab, row))
1128                         continue;
1129                 return row;
1130         }
1131         return -1;
1132 }
1133
1134 /* Check whether the invariant that all columns are lexico-positive
1135  * is satisfied.  This function is not called from the current code
1136  * but is useful during debugging.
1137  */
1138 static void check_lexpos(struct isl_tab *tab) __attribute__ ((unused));
1139 static void check_lexpos(struct isl_tab *tab)
1140 {
1141         unsigned off = 2 + tab->M;
1142         int col;
1143         int var;
1144         int row;
1145
1146         for (col = tab->n_dead; col < tab->n_col; ++col) {
1147                 if (tab->col_var[col] >= 0 &&
1148                     (tab->col_var[col] < tab->n_param ||
1149                      tab->col_var[col] >= tab->n_var - tab->n_div))
1150                         continue;
1151                 for (var = tab->n_param; var < tab->n_var - tab->n_div; ++var) {
1152                         if (!tab->var[var].is_row) {
1153                                 if (tab->var[var].index == col)
1154                                         break;
1155                                 else
1156                                         continue;
1157                         }
1158                         row = tab->var[var].index;
1159                         if (isl_int_is_zero(tab->mat->row[row][off + col]))
1160                                 continue;
1161                         if (isl_int_is_pos(tab->mat->row[row][off + col]))
1162                                 break;
1163                         fprintf(stderr, "lexneg column %d (row %d)\n",
1164                                 col, row);
1165                 }
1166                 if (var >= tab->n_var - tab->n_div)
1167                         fprintf(stderr, "zero column %d\n", col);
1168         }
1169 }
1170
1171 /* Report to the caller that the given constraint is part of an encountered
1172  * conflict.
1173  */
1174 static int report_conflicting_constraint(struct isl_tab *tab, int con)
1175 {
1176         return tab->conflict(con, tab->conflict_user);
1177 }
1178
1179 /* Given a conflicting row in the tableau, report all constraints
1180  * involved in the row to the caller.  That is, the row itself
1181  * (if it represents a constraint) and all constraint columns with
1182  * non-zero (and therefore negative) coefficients.
1183  */
1184 static int report_conflict(struct isl_tab *tab, int row)
1185 {
1186         int j;
1187         isl_int *tr;
1188
1189         if (!tab->conflict)
1190                 return 0;
1191
1192         if (tab->row_var[row] < 0 &&
1193             report_conflicting_constraint(tab, ~tab->row_var[row]) < 0)
1194                 return -1;
1195
1196         tr = tab->mat->row[row] + 2 + tab->M;
1197
1198         for (j = tab->n_dead; j < tab->n_col; ++j) {
1199                 if (tab->col_var[j] >= 0 &&
1200                     (tab->col_var[j] < tab->n_param  ||
1201                     tab->col_var[j] >= tab->n_var - tab->n_div))
1202                         continue;
1203
1204                 if (!isl_int_is_neg(tr[j]))
1205                         continue;
1206
1207                 if (tab->col_var[j] < 0 &&
1208                     report_conflicting_constraint(tab, ~tab->col_var[j]) < 0)
1209                         return -1;
1210         }
1211
1212         return 0;
1213 }
1214
1215 /* Resolve all known or obviously violated constraints through pivoting.
1216  * In particular, as long as we can find any violated constraint, we
1217  * look for a pivoting column that would result in the lexicographically
1218  * smallest increment in the sample point.  If there is no such column
1219  * then the tableau is infeasible.
1220  */
1221 static int restore_lexmin(struct isl_tab *tab) WARN_UNUSED;
1222 static int restore_lexmin(struct isl_tab *tab)
1223 {
1224         int row, col;
1225
1226         if (!tab)
1227                 return -1;
1228         if (tab->empty)
1229                 return 0;
1230         while ((row = first_neg(tab)) != -1) {
1231                 col = lexmin_pivot_col(tab, row);
1232                 if (col >= tab->n_col) {
1233                         if (report_conflict(tab, row) < 0)
1234                                 return -1;
1235                         if (isl_tab_mark_empty(tab) < 0)
1236                                 return -1;
1237                         return 0;
1238                 }
1239                 if (col < 0)
1240                         return -1;
1241                 if (isl_tab_pivot(tab, row, col) < 0)
1242                         return -1;
1243         }
1244         return 0;
1245 }
1246
1247 /* Given a row that represents an equality, look for an appropriate
1248  * pivoting column.
1249  * In particular, if there are any non-zero coefficients among
1250  * the non-parameter variables, then we take the last of these
1251  * variables.  Eliminating this variable in terms of the other
1252  * variables and/or parameters does not influence the property
1253  * that all column in the initial tableau are lexicographically
1254  * positive.  The row corresponding to the eliminated variable
1255  * will only have non-zero entries below the diagonal of the
1256  * initial tableau.  That is, we transform
1257  *
1258  *              I                               I
1259  *                1             into            a
1260  *                  I                             I
1261  *
1262  * If there is no such non-parameter variable, then we are dealing with
1263  * pure parameter equality and we pick any parameter with coefficient 1 or -1
1264  * for elimination.  This will ensure that the eliminated parameter
1265  * always has an integer value whenever all the other parameters are integral.
1266  * If there is no such parameter then we return -1.
1267  */
1268 static int last_var_col_or_int_par_col(struct isl_tab *tab, int row)
1269 {
1270         unsigned off = 2 + tab->M;
1271         int i;
1272
1273         for (i = tab->n_var - tab->n_div - 1; i >= 0 && i >= tab->n_param; --i) {
1274                 int col;
1275                 if (tab->var[i].is_row)
1276                         continue;
1277                 col = tab->var[i].index;
1278                 if (col <= tab->n_dead)
1279                         continue;
1280                 if (!isl_int_is_zero(tab->mat->row[row][off + col]))
1281                         return col;
1282         }
1283         for (i = tab->n_dead; i < tab->n_col; ++i) {
1284                 if (isl_int_is_one(tab->mat->row[row][off + i]))
1285                         return i;
1286                 if (isl_int_is_negone(tab->mat->row[row][off + i]))
1287                         return i;
1288         }
1289         return -1;
1290 }
1291
1292 /* Add an equality that is known to be valid to the tableau.
1293  * We first check if we can eliminate a variable or a parameter.
1294  * If not, we add the equality as two inequalities.
1295  * In this case, the equality was a pure parameter equality and there
1296  * is no need to resolve any constraint violations.
1297  */
1298 static struct isl_tab *add_lexmin_valid_eq(struct isl_tab *tab, isl_int *eq)
1299 {
1300         int i;
1301         int r;
1302
1303         if (!tab)
1304                 return NULL;
1305         r = isl_tab_add_row(tab, eq);
1306         if (r < 0)
1307                 goto error;
1308
1309         r = tab->con[r].index;
1310         i = last_var_col_or_int_par_col(tab, r);
1311         if (i < 0) {
1312                 tab->con[r].is_nonneg = 1;
1313                 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1314                         goto error;
1315                 isl_seq_neg(eq, eq, 1 + tab->n_var);
1316                 r = isl_tab_add_row(tab, eq);
1317                 if (r < 0)
1318                         goto error;
1319                 tab->con[r].is_nonneg = 1;
1320                 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1321                         goto error;
1322         } else {
1323                 if (isl_tab_pivot(tab, r, i) < 0)
1324                         goto error;
1325                 if (isl_tab_kill_col(tab, i) < 0)
1326                         goto error;
1327                 tab->n_eq++;
1328         }
1329
1330         return tab;
1331 error:
1332         isl_tab_free(tab);
1333         return NULL;
1334 }
1335
1336 /* Check if the given row is a pure constant.
1337  */
1338 static int is_constant(struct isl_tab *tab, int row)
1339 {
1340         unsigned off = 2 + tab->M;
1341
1342         return isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
1343                                         tab->n_col - tab->n_dead) == -1;
1344 }
1345
1346 /* Add an equality that may or may not be valid to the tableau.
1347  * If the resulting row is a pure constant, then it must be zero.
1348  * Otherwise, the resulting tableau is empty.
1349  *
1350  * If the row is not a pure constant, then we add two inequalities,
1351  * each time checking that they can be satisfied.
1352  * In the end we try to use one of the two constraints to eliminate
1353  * a column.
1354  */
1355 static int add_lexmin_eq(struct isl_tab *tab, isl_int *eq) WARN_UNUSED;
1356 static int add_lexmin_eq(struct isl_tab *tab, isl_int *eq)
1357 {
1358         int r1, r2;
1359         int row;
1360         struct isl_tab_undo *snap;
1361
1362         if (!tab)
1363                 return -1;
1364         snap = isl_tab_snap(tab);
1365         r1 = isl_tab_add_row(tab, eq);
1366         if (r1 < 0)
1367                 return -1;
1368         tab->con[r1].is_nonneg = 1;
1369         if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r1]) < 0)
1370                 return -1;
1371
1372         row = tab->con[r1].index;
1373         if (is_constant(tab, row)) {
1374                 if (!isl_int_is_zero(tab->mat->row[row][1]) ||
1375                     (tab->M && !isl_int_is_zero(tab->mat->row[row][2]))) {
1376                         if (isl_tab_mark_empty(tab) < 0)
1377                                 return -1;
1378                         return 0;
1379                 }
1380                 if (isl_tab_rollback(tab, snap) < 0)
1381                         return -1;
1382                 return 0;
1383         }
1384
1385         if (restore_lexmin(tab) < 0)
1386                 return -1;
1387         if (tab->empty)
1388                 return 0;
1389
1390         isl_seq_neg(eq, eq, 1 + tab->n_var);
1391
1392         r2 = isl_tab_add_row(tab, eq);
1393         if (r2 < 0)
1394                 return -1;
1395         tab->con[r2].is_nonneg = 1;
1396         if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r2]) < 0)
1397                 return -1;
1398
1399         if (restore_lexmin(tab) < 0)
1400                 return -1;
1401         if (tab->empty)
1402                 return 0;
1403
1404         if (!tab->con[r1].is_row) {
1405                 if (isl_tab_kill_col(tab, tab->con[r1].index) < 0)
1406                         return -1;
1407         } else if (!tab->con[r2].is_row) {
1408                 if (isl_tab_kill_col(tab, tab->con[r2].index) < 0)
1409                         return -1;
1410         }
1411
1412         if (tab->bmap) {
1413                 tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);
1414                 if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1415                         return -1;
1416                 isl_seq_neg(eq, eq, 1 + tab->n_var);
1417                 tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);
1418                 isl_seq_neg(eq, eq, 1 + tab->n_var);
1419                 if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1420                         return -1;
1421                 if (!tab->bmap)
1422                         return -1;
1423         }
1424
1425         return 0;
1426 }
1427
1428 /* Add an inequality to the tableau, resolving violations using
1429  * restore_lexmin.
1430  */
1431 static struct isl_tab *add_lexmin_ineq(struct isl_tab *tab, isl_int *ineq)
1432 {
1433         int r;
1434
1435         if (!tab)
1436                 return NULL;
1437         if (tab->bmap) {
1438                 tab->bmap = isl_basic_map_add_ineq(tab->bmap, ineq);
1439                 if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1440                         goto error;
1441                 if (!tab->bmap)
1442                         goto error;
1443         }
1444         r = isl_tab_add_row(tab, ineq);
1445         if (r < 0)
1446                 goto error;
1447         tab->con[r].is_nonneg = 1;
1448         if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1449                 goto error;
1450         if (isl_tab_row_is_redundant(tab, tab->con[r].index)) {
1451                 if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0)
1452                         goto error;
1453                 return tab;
1454         }
1455
1456         if (restore_lexmin(tab) < 0)
1457                 goto error;
1458         if (!tab->empty && tab->con[r].is_row &&
1459                  isl_tab_row_is_redundant(tab, tab->con[r].index))
1460                 if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0)
1461                         goto error;
1462         return tab;
1463 error:
1464         isl_tab_free(tab);
1465         return NULL;
1466 }
1467
1468 /* Check if the coefficients of the parameters are all integral.
1469  */
1470 static int integer_parameter(struct isl_tab *tab, int row)
1471 {
1472         int i;
1473         int col;
1474         unsigned off = 2 + tab->M;
1475
1476         for (i = 0; i < tab->n_param; ++i) {
1477                 /* Eliminated parameter */
1478                 if (tab->var[i].is_row)
1479                         continue;
1480                 col = tab->var[i].index;
1481                 if (!isl_int_is_divisible_by(tab->mat->row[row][off + col],
1482                                                 tab->mat->row[row][0]))
1483                         return 0;
1484         }
1485         for (i = 0; i < tab->n_div; ++i) {
1486                 if (tab->var[tab->n_var - tab->n_div + i].is_row)
1487                         continue;
1488                 col = tab->var[tab->n_var - tab->n_div + i].index;
1489                 if (!isl_int_is_divisible_by(tab->mat->row[row][off + col],
1490                                                 tab->mat->row[row][0]))
1491                         return 0;
1492         }
1493         return 1;
1494 }
1495
1496 /* Check if the coefficients of the non-parameter variables are all integral.
1497  */
1498 static int integer_variable(struct isl_tab *tab, int row)
1499 {
1500         int i;
1501         unsigned off = 2 + tab->M;
1502
1503         for (i = tab->n_dead; i < tab->n_col; ++i) {
1504                 if (tab->col_var[i] >= 0 &&
1505                     (tab->col_var[i] < tab->n_param ||
1506                      tab->col_var[i] >= tab->n_var - tab->n_div))
1507                         continue;
1508                 if (!isl_int_is_divisible_by(tab->mat->row[row][off + i],
1509                                                 tab->mat->row[row][0]))
1510                         return 0;
1511         }
1512         return 1;
1513 }
1514
1515 /* Check if the constant term is integral.
1516  */
1517 static int integer_constant(struct isl_tab *tab, int row)
1518 {
1519         return isl_int_is_divisible_by(tab->mat->row[row][1],
1520                                         tab->mat->row[row][0]);
1521 }
1522
1523 #define I_CST   1 << 0
1524 #define I_PAR   1 << 1
1525 #define I_VAR   1 << 2
1526
1527 /* Check for next (non-parameter) variable after "var" (first if var == -1)
1528  * that is non-integer and therefore requires a cut and return
1529  * the index of the variable.
1530  * For parametric tableaus, there are three parts in a row,
1531  * the constant, the coefficients of the parameters and the rest.
1532  * For each part, we check whether the coefficients in that part
1533  * are all integral and if so, set the corresponding flag in *f.
1534  * If the constant and the parameter part are integral, then the
1535  * current sample value is integral and no cut is required
1536  * (irrespective of whether the variable part is integral).
1537  */
1538 static int next_non_integer_var(struct isl_tab *tab, int var, int *f)
1539 {
1540         var = var < 0 ? tab->n_param : var + 1;
1541
1542         for (; var < tab->n_var - tab->n_div; ++var) {
1543                 int flags = 0;
1544                 int row;
1545                 if (!tab->var[var].is_row)
1546                         continue;
1547                 row = tab->var[var].index;
1548                 if (integer_constant(tab, row))
1549                         ISL_FL_SET(flags, I_CST);
1550                 if (integer_parameter(tab, row))
1551                         ISL_FL_SET(flags, I_PAR);
1552                 if (ISL_FL_ISSET(flags, I_CST) && ISL_FL_ISSET(flags, I_PAR))
1553                         continue;
1554                 if (integer_variable(tab, row))
1555                         ISL_FL_SET(flags, I_VAR);
1556                 *f = flags;
1557                 return var;
1558         }
1559         return -1;
1560 }
1561
1562 /* Check for first (non-parameter) variable that is non-integer and
1563  * therefore requires a cut and return the corresponding row.
1564  * For parametric tableaus, there are three parts in a row,
1565  * the constant, the coefficients of the parameters and the rest.
1566  * For each part, we check whether the coefficients in that part
1567  * are all integral and if so, set the corresponding flag in *f.
1568  * If the constant and the parameter part are integral, then the
1569  * current sample value is integral and no cut is required
1570  * (irrespective of whether the variable part is integral).
1571  */
1572 static int first_non_integer_row(struct isl_tab *tab, int *f)
1573 {
1574         int var = next_non_integer_var(tab, -1, f);
1575
1576         return var < 0 ? -1 : tab->var[var].index;
1577 }
1578
1579 /* Add a (non-parametric) cut to cut away the non-integral sample
1580  * value of the given row.
1581  *
1582  * If the row is given by
1583  *
1584  *      m r = f + \sum_i a_i y_i
1585  *
1586  * then the cut is
1587  *
1588  *      c = - {-f/m} + \sum_i {a_i/m} y_i >= 0
1589  *
1590  * The big parameter, if any, is ignored, since it is assumed to be big
1591  * enough to be divisible by any integer.
1592  * If the tableau is actually a parametric tableau, then this function
1593  * is only called when all coefficients of the parameters are integral.
1594  * The cut therefore has zero coefficients for the parameters.
1595  *
1596  * The current value is known to be negative, so row_sign, if it
1597  * exists, is set accordingly.
1598  *
1599  * Return the row of the cut or -1.
1600  */
1601 static int add_cut(struct isl_tab *tab, int row)
1602 {
1603         int i;
1604         int r;
1605         isl_int *r_row;
1606         unsigned off = 2 + tab->M;
1607
1608         if (isl_tab_extend_cons(tab, 1) < 0)
1609                 return -1;
1610         r = isl_tab_allocate_con(tab);
1611         if (r < 0)
1612                 return -1;
1613
1614         r_row = tab->mat->row[tab->con[r].index];
1615         isl_int_set(r_row[0], tab->mat->row[row][0]);
1616         isl_int_neg(r_row[1], tab->mat->row[row][1]);
1617         isl_int_fdiv_r(r_row[1], r_row[1], tab->mat->row[row][0]);
1618         isl_int_neg(r_row[1], r_row[1]);
1619         if (tab->M)
1620                 isl_int_set_si(r_row[2], 0);
1621         for (i = 0; i < tab->n_col; ++i)
1622                 isl_int_fdiv_r(r_row[off + i],
1623                         tab->mat->row[row][off + i], tab->mat->row[row][0]);
1624
1625         tab->con[r].is_nonneg = 1;
1626         if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1627                 return -1;
1628         if (tab->row_sign)
1629                 tab->row_sign[tab->con[r].index] = isl_tab_row_neg;
1630
1631         return tab->con[r].index;
1632 }
1633
1634 #define CUT_ALL 1
1635 #define CUT_ONE 0
1636
1637 /* Given a non-parametric tableau, add cuts until an integer
1638  * sample point is obtained or until the tableau is determined
1639  * to be integer infeasible.
1640  * As long as there is any non-integer value in the sample point,
1641  * we add appropriate cuts, if possible, for each of these
1642  * non-integer values and then resolve the violated
1643  * cut constraints using restore_lexmin.
1644  * If one of the corresponding rows is equal to an integral
1645  * combination of variables/constraints plus a non-integral constant,
1646  * then there is no way to obtain an integer point and we return
1647  * a tableau that is marked empty.
1648  * The parameter cutting_strategy controls the strategy used when adding cuts
1649  * to remove non-integer points. CUT_ALL adds all possible cuts
1650  * before continuing the search. CUT_ONE adds only one cut at a time.
1651  */
1652 static struct isl_tab *cut_to_integer_lexmin(struct isl_tab *tab,
1653         int cutting_strategy)
1654 {
1655         int var;
1656         int row;
1657         int flags;
1658
1659         if (!tab)
1660                 return NULL;
1661         if (tab->empty)
1662                 return tab;
1663
1664         while ((var = next_non_integer_var(tab, -1, &flags)) != -1) {
1665                 do {
1666                         if (ISL_FL_ISSET(flags, I_VAR)) {
1667                                 if (isl_tab_mark_empty(tab) < 0)
1668                                         goto error;
1669                                 return tab;
1670                         }
1671                         row = tab->var[var].index;
1672                         row = add_cut(tab, row);
1673                         if (row < 0)
1674                                 goto error;
1675                         if (cutting_strategy == CUT_ONE)
1676                                 break;
1677                 } while ((var = next_non_integer_var(tab, var, &flags)) != -1);
1678                 if (restore_lexmin(tab) < 0)
1679                         goto error;
1680                 if (tab->empty)
1681                         break;
1682         }
1683         return tab;
1684 error:
1685         isl_tab_free(tab);
1686         return NULL;
1687 }
1688
1689 /* Check whether all the currently active samples also satisfy the inequality
1690  * "ineq" (treated as an equality if eq is set).
1691  * Remove those samples that do not.
1692  */
1693 static struct isl_tab *check_samples(struct isl_tab *tab, isl_int *ineq, int eq)
1694 {
1695         int i;
1696         isl_int v;
1697
1698         if (!tab)
1699                 return NULL;
1700
1701         isl_assert(tab->mat->ctx, tab->bmap, goto error);
1702         isl_assert(tab->mat->ctx, tab->samples, goto error);
1703         isl_assert(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var, goto error);
1704
1705         isl_int_init(v);
1706         for (i = tab->n_outside; i < tab->n_sample; ++i) {
1707                 int sgn;
1708                 isl_seq_inner_product(ineq, tab->samples->row[i],
1709                                         1 + tab->n_var, &v);
1710                 sgn = isl_int_sgn(v);
1711                 if (eq ? (sgn == 0) : (sgn >= 0))
1712                         continue;
1713                 tab = isl_tab_drop_sample(tab, i);
1714                 if (!tab)
1715                         break;
1716         }
1717         isl_int_clear(v);
1718
1719         return tab;
1720 error:
1721         isl_tab_free(tab);
1722         return NULL;
1723 }
1724
1725 /* Check whether the sample value of the tableau is finite,
1726  * i.e., either the tableau does not use a big parameter, or
1727  * all values of the variables are equal to the big parameter plus
1728  * some constant.  This constant is the actual sample value.
1729  */
1730 static int sample_is_finite(struct isl_tab *tab)
1731 {
1732         int i;
1733
1734         if (!tab->M)
1735                 return 1;
1736
1737         for (i = 0; i < tab->n_var; ++i) {
1738                 int row;
1739                 if (!tab->var[i].is_row)
1740                         return 0;
1741                 row = tab->var[i].index;
1742                 if (isl_int_ne(tab->mat->row[row][0], tab->mat->row[row][2]))
1743                         return 0;
1744         }
1745         return 1;
1746 }
1747
1748 /* Check if the context tableau of sol has any integer points.
1749  * Leave tab in empty state if no integer point can be found.
1750  * If an integer point can be found and if moreover it is finite,
1751  * then it is added to the list of sample values.
1752  *
1753  * This function is only called when none of the currently active sample
1754  * values satisfies the most recently added constraint.
1755  */
1756 static struct isl_tab *check_integer_feasible(struct isl_tab *tab)
1757 {
1758         struct isl_tab_undo *snap;
1759
1760         if (!tab)
1761                 return NULL;
1762
1763         snap = isl_tab_snap(tab);
1764         if (isl_tab_push_basis(tab) < 0)
1765                 goto error;
1766
1767         tab = cut_to_integer_lexmin(tab, CUT_ALL);
1768         if (!tab)
1769                 goto error;
1770
1771         if (!tab->empty && sample_is_finite(tab)) {
1772                 struct isl_vec *sample;
1773
1774                 sample = isl_tab_get_sample_value(tab);
1775
1776                 tab = isl_tab_add_sample(tab, sample);
1777         }
1778
1779         if (!tab->empty && isl_tab_rollback(tab, snap) < 0)
1780                 goto error;
1781
1782         return tab;
1783 error:
1784         isl_tab_free(tab);
1785         return NULL;
1786 }
1787
1788 /* Check if any of the currently active sample values satisfies
1789  * the inequality "ineq" (an equality if eq is set).
1790  */
1791 static int tab_has_valid_sample(struct isl_tab *tab, isl_int *ineq, int eq)
1792 {
1793         int i;
1794         isl_int v;
1795
1796         if (!tab)
1797                 return -1;
1798
1799         isl_assert(tab->mat->ctx, tab->bmap, return -1);
1800         isl_assert(tab->mat->ctx, tab->samples, return -1);
1801         isl_assert(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var, return -1);
1802
1803         isl_int_init(v);
1804         for (i = tab->n_outside; i < tab->n_sample; ++i) {
1805                 int sgn;
1806                 isl_seq_inner_product(ineq, tab->samples->row[i],
1807                                         1 + tab->n_var, &v);
1808                 sgn = isl_int_sgn(v);
1809                 if (eq ? (sgn == 0) : (sgn >= 0))
1810                         break;
1811         }
1812         isl_int_clear(v);
1813
1814         return i < tab->n_sample;
1815 }
1816
1817 /* Add a div specified by "div" to the tableau "tab" and return
1818  * 1 if the div is obviously non-negative.
1819  */
1820 static int context_tab_add_div(struct isl_tab *tab, struct isl_vec *div,
1821         int (*add_ineq)(void *user, isl_int *), void *user)
1822 {
1823         int i;
1824         int r;
1825         struct isl_mat *samples;
1826         int nonneg;
1827
1828         r = isl_tab_add_div(tab, div, add_ineq, user);
1829         if (r < 0)
1830                 return -1;
1831         nonneg = tab->var[r].is_nonneg;
1832         tab->var[r].frozen = 1;
1833
1834         samples = isl_mat_extend(tab->samples,
1835                         tab->n_sample, 1 + tab->n_var);
1836         tab->samples = samples;
1837         if (!samples)
1838                 return -1;
1839         for (i = tab->n_outside; i < samples->n_row; ++i) {
1840                 isl_seq_inner_product(div->el + 1, samples->row[i],
1841                         div->size - 1, &samples->row[i][samples->n_col - 1]);
1842                 isl_int_fdiv_q(samples->row[i][samples->n_col - 1],
1843                                samples->row[i][samples->n_col - 1], div->el[0]);
1844         }
1845
1846         return nonneg;
1847 }
1848
1849 /* Add a div specified by "div" to both the main tableau and
1850  * the context tableau.  In case of the main tableau, we only
1851  * need to add an extra div.  In the context tableau, we also
1852  * need to express the meaning of the div.
1853  * Return the index of the div or -1 if anything went wrong.
1854  */
1855 static int add_div(struct isl_tab *tab, struct isl_context *context,
1856         struct isl_vec *div)
1857 {
1858         int r;
1859         int nonneg;
1860
1861         if ((nonneg = context->op->add_div(context, div)) < 0)
1862                 goto error;
1863
1864         if (!context->op->is_ok(context))
1865                 goto error;
1866
1867         if (isl_tab_extend_vars(tab, 1) < 0)
1868                 goto error;
1869         r = isl_tab_allocate_var(tab);
1870         if (r < 0)
1871                 goto error;
1872         if (nonneg)
1873                 tab->var[r].is_nonneg = 1;
1874         tab->var[r].frozen = 1;
1875         tab->n_div++;
1876
1877         return tab->n_div - 1;
1878 error:
1879         context->op->invalidate(context);
1880         return -1;
1881 }
1882
1883 static int find_div(struct isl_tab *tab, isl_int *div, isl_int denom)
1884 {
1885         int i;
1886         unsigned total = isl_basic_map_total_dim(tab->bmap);
1887
1888         for (i = 0; i < tab->bmap->n_div; ++i) {
1889                 if (isl_int_ne(tab->bmap->div[i][0], denom))
1890                         continue;
1891                 if (!isl_seq_eq(tab->bmap->div[i] + 1, div, 1 + total))
1892                         continue;
1893                 return i;
1894         }
1895         return -1;
1896 }
1897
1898 /* Return the index of a div that corresponds to "div".
1899  * We first check if we already have such a div and if not, we create one.
1900  */
1901 static int get_div(struct isl_tab *tab, struct isl_context *context,
1902         struct isl_vec *div)
1903 {
1904         int d;
1905         struct isl_tab *context_tab = context->op->peek_tab(context);
1906
1907         if (!context_tab)
1908                 return -1;
1909
1910         d = find_div(context_tab, div->el + 1, div->el[0]);
1911         if (d != -1)
1912                 return d;
1913
1914         return add_div(tab, context, div);
1915 }
1916
1917 /* Add a parametric cut to cut away the non-integral sample value
1918  * of the give row.
1919  * Let a_i be the coefficients of the constant term and the parameters
1920  * and let b_i be the coefficients of the variables or constraints
1921  * in basis of the tableau.
1922  * Let q be the div q = floor(\sum_i {-a_i} y_i).
1923  *
1924  * The cut is expressed as
1925  *
1926  *      c = \sum_i -{-a_i} y_i + \sum_i {b_i} x_i + q >= 0
1927  *
1928  * If q did not already exist in the context tableau, then it is added first.
1929  * If q is in a column of the main tableau then the "+ q" can be accomplished
1930  * by setting the corresponding entry to the denominator of the constraint.
1931  * If q happens to be in a row of the main tableau, then the corresponding
1932  * row needs to be added instead (taking care of the denominators).
1933  * Note that this is very unlikely, but perhaps not entirely impossible.
1934  *
1935  * The current value of the cut is known to be negative (or at least
1936  * non-positive), so row_sign is set accordingly.
1937  *
1938  * Return the row of the cut or -1.
1939  */
1940 static int add_parametric_cut(struct isl_tab *tab, int row,
1941         struct isl_context *context)
1942 {
1943         struct isl_vec *div;
1944         int d;
1945         int i;
1946         int r;
1947         isl_int *r_row;
1948         int col;
1949         int n;
1950         unsigned off = 2 + tab->M;
1951
1952         if (!context)
1953                 return -1;
1954
1955         div = get_row_parameter_div(tab, row);
1956         if (!div)
1957                 return -1;
1958
1959         n = tab->n_div;
1960         d = context->op->get_div(context, tab, div);
1961         if (d < 0)
1962                 return -1;
1963
1964         if (isl_tab_extend_cons(tab, 1) < 0)
1965                 return -1;
1966         r = isl_tab_allocate_con(tab);
1967         if (r < 0)
1968                 return -1;
1969
1970         r_row = tab->mat->row[tab->con[r].index];
1971         isl_int_set(r_row[0], tab->mat->row[row][0]);
1972         isl_int_neg(r_row[1], tab->mat->row[row][1]);
1973         isl_int_fdiv_r(r_row[1], r_row[1], tab->mat->row[row][0]);
1974         isl_int_neg(r_row[1], r_row[1]);
1975         if (tab->M)
1976                 isl_int_set_si(r_row[2], 0);
1977         for (i = 0; i < tab->n_param; ++i) {
1978                 if (tab->var[i].is_row)
1979                         continue;
1980                 col = tab->var[i].index;
1981                 isl_int_neg(r_row[off + col], tab->mat->row[row][off + col]);
1982                 isl_int_fdiv_r(r_row[off + col], r_row[off + col],
1983                                 tab->mat->row[row][0]);
1984                 isl_int_neg(r_row[off + col], r_row[off + col]);
1985         }
1986         for (i = 0; i < tab->n_div; ++i) {
1987                 if (tab->var[tab->n_var - tab->n_div + i].is_row)
1988                         continue;
1989                 col = tab->var[tab->n_var - tab->n_div + i].index;
1990                 isl_int_neg(r_row[off + col], tab->mat->row[row][off + col]);
1991                 isl_int_fdiv_r(r_row[off + col], r_row[off + col],
1992                                 tab->mat->row[row][0]);
1993                 isl_int_neg(r_row[off + col], r_row[off + col]);
1994         }
1995         for (i = 0; i < tab->n_col; ++i) {
1996                 if (tab->col_var[i] >= 0 &&
1997                     (tab->col_var[i] < tab->n_param ||
1998                      tab->col_var[i] >= tab->n_var - tab->n_div))
1999                         continue;
2000                 isl_int_fdiv_r(r_row[off + i],
2001                         tab->mat->row[row][off + i], tab->mat->row[row][0]);
2002         }
2003         if (tab->var[tab->n_var - tab->n_div + d].is_row) {
2004                 isl_int gcd;
2005                 int d_row = tab->var[tab->n_var - tab->n_div + d].index;
2006                 isl_int_init(gcd);
2007                 isl_int_gcd(gcd, tab->mat->row[d_row][0], r_row[0]);
2008                 isl_int_divexact(r_row[0], r_row[0], gcd);
2009                 isl_int_divexact(gcd, tab->mat->row[d_row][0], gcd);
2010                 isl_seq_combine(r_row + 1, gcd, r_row + 1,
2011                                 r_row[0], tab->mat->row[d_row] + 1,
2012                                 off - 1 + tab->n_col);
2013                 isl_int_mul(r_row[0], r_row[0], tab->mat->row[d_row][0]);
2014                 isl_int_clear(gcd);
2015         } else {
2016                 col = tab->var[tab->n_var - tab->n_div + d].index;
2017                 isl_int_set(r_row[off + col], tab->mat->row[row][0]);
2018         }
2019
2020         tab->con[r].is_nonneg = 1;
2021         if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
2022                 return -1;
2023         if (tab->row_sign)
2024                 tab->row_sign[tab->con[r].index] = isl_tab_row_neg;
2025
2026         isl_vec_free(div);
2027
2028         row = tab->con[r].index;
2029
2030         if (d >= n && context->op->detect_equalities(context, tab) < 0)
2031                 return -1;
2032
2033         return row;
2034 }
2035
2036 /* Construct a tableau for bmap that can be used for computing
2037  * the lexicographic minimum (or maximum) of bmap.
2038  * If not NULL, then dom is the domain where the minimum
2039  * should be computed.  In this case, we set up a parametric
2040  * tableau with row signs (initialized to "unknown").
2041  * If M is set, then the tableau will use a big parameter.
2042  * If max is set, then a maximum should be computed instead of a minimum.
2043  * This means that for each variable x, the tableau will contain the variable
2044  * x' = M - x, rather than x' = M + x.  This in turn means that the coefficient
2045  * of the variables in all constraints are negated prior to adding them
2046  * to the tableau.
2047  */
2048 static struct isl_tab *tab_for_lexmin(struct isl_basic_map *bmap,
2049         struct isl_basic_set *dom, unsigned M, int max)
2050 {
2051         int i;
2052         struct isl_tab *tab;
2053
2054         tab = isl_tab_alloc(bmap->ctx, 2 * bmap->n_eq + bmap->n_ineq + 1,
2055                             isl_basic_map_total_dim(bmap), M);
2056         if (!tab)
2057                 return NULL;
2058
2059         tab->rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
2060         if (dom) {
2061                 tab->n_param = isl_basic_set_total_dim(dom) - dom->n_div;
2062                 tab->n_div = dom->n_div;
2063                 tab->row_sign = isl_calloc_array(bmap->ctx,
2064                                         enum isl_tab_row_sign, tab->mat->n_row);
2065                 if (!tab->row_sign)
2066                         goto error;
2067         }
2068         if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY)) {
2069                 if (isl_tab_mark_empty(tab) < 0)
2070                         goto error;
2071                 return tab;
2072         }
2073
2074         for (i = tab->n_param; i < tab->n_var - tab->n_div; ++i) {
2075                 tab->var[i].is_nonneg = 1;
2076                 tab->var[i].frozen = 1;
2077         }
2078         for (i = 0; i < bmap->n_eq; ++i) {
2079                 if (max)
2080                         isl_seq_neg(bmap->eq[i] + 1 + tab->n_param,
2081                                     bmap->eq[i] + 1 + tab->n_param,
2082                                     tab->n_var - tab->n_param - tab->n_div);
2083                 tab = add_lexmin_valid_eq(tab, bmap->eq[i]);
2084                 if (max)
2085                         isl_seq_neg(bmap->eq[i] + 1 + tab->n_param,
2086                                     bmap->eq[i] + 1 + tab->n_param,
2087                                     tab->n_var - tab->n_param - tab->n_div);
2088                 if (!tab || tab->empty)
2089                         return tab;
2090         }
2091         if (bmap->n_eq && restore_lexmin(tab) < 0)
2092                 goto error;
2093         for (i = 0; i < bmap->n_ineq; ++i) {
2094                 if (max)
2095                         isl_seq_neg(bmap->ineq[i] + 1 + tab->n_param,
2096                                     bmap->ineq[i] + 1 + tab->n_param,
2097                                     tab->n_var - tab->n_param - tab->n_div);
2098                 tab = add_lexmin_ineq(tab, bmap->ineq[i]);
2099                 if (max)
2100                         isl_seq_neg(bmap->ineq[i] + 1 + tab->n_param,
2101                                     bmap->ineq[i] + 1 + tab->n_param,
2102                                     tab->n_var - tab->n_param - tab->n_div);
2103                 if (!tab || tab->empty)
2104                         return tab;
2105         }
2106         return tab;
2107 error:
2108         isl_tab_free(tab);
2109         return NULL;
2110 }
2111
2112 /* Given a main tableau where more than one row requires a split,
2113  * determine and return the "best" row to split on.
2114  *
2115  * Given two rows in the main tableau, if the inequality corresponding
2116  * to the first row is redundant with respect to that of the second row
2117  * in the current tableau, then it is better to split on the second row,
2118  * since in the positive part, both row will be positive.
2119  * (In the negative part a pivot will have to be performed and just about
2120  * anything can happen to the sign of the other row.)
2121  *
2122  * As a simple heuristic, we therefore select the row that makes the most
2123  * of the other rows redundant.
2124  *
2125  * Perhaps it would also be useful to look at the number of constraints
2126  * that conflict with any given constraint.
2127  */
2128 static int best_split(struct isl_tab *tab, struct isl_tab *context_tab)
2129 {
2130         struct isl_tab_undo *snap;
2131         int split;
2132         int row;
2133         int best = -1;
2134         int best_r;
2135
2136         if (isl_tab_extend_cons(context_tab, 2) < 0)
2137                 return -1;
2138
2139         snap = isl_tab_snap(context_tab);
2140
2141         for (split = tab->n_redundant; split < tab->n_row; ++split) {
2142                 struct isl_tab_undo *snap2;
2143                 struct isl_vec *ineq = NULL;
2144                 int r = 0;
2145                 int ok;
2146
2147                 if (!isl_tab_var_from_row(tab, split)->is_nonneg)
2148                         continue;
2149                 if (tab->row_sign[split] != isl_tab_row_any)
2150                         continue;
2151
2152                 ineq = get_row_parameter_ineq(tab, split);
2153                 if (!ineq)
2154                         return -1;
2155                 ok = isl_tab_add_ineq(context_tab, ineq->el) >= 0;
2156                 isl_vec_free(ineq);
2157                 if (!ok)
2158                         return -1;
2159
2160                 snap2 = isl_tab_snap(context_tab);
2161
2162                 for (row = tab->n_redundant; row < tab->n_row; ++row) {
2163                         struct isl_tab_var *var;
2164
2165                         if (row == split)
2166                                 continue;
2167                         if (!isl_tab_var_from_row(tab, row)->is_nonneg)
2168                                 continue;
2169                         if (tab->row_sign[row] != isl_tab_row_any)
2170                                 continue;
2171
2172                         ineq = get_row_parameter_ineq(tab, row);
2173                         if (!ineq)
2174                                 return -1;
2175                         ok = isl_tab_add_ineq(context_tab, ineq->el) >= 0;
2176                         isl_vec_free(ineq);
2177                         if (!ok)
2178                                 return -1;
2179                         var = &context_tab->con[context_tab->n_con - 1];
2180                         if (!context_tab->empty &&
2181                             !isl_tab_min_at_most_neg_one(context_tab, var))
2182                                 r++;
2183                         if (isl_tab_rollback(context_tab, snap2) < 0)
2184                                 return -1;
2185                 }
2186                 if (best == -1 || r > best_r) {
2187                         best = split;
2188                         best_r = r;
2189                 }
2190                 if (isl_tab_rollback(context_tab, snap) < 0)
2191                         return -1;
2192         }
2193
2194         return best;
2195 }
2196
2197 static struct isl_basic_set *context_lex_peek_basic_set(
2198         struct isl_context *context)
2199 {
2200         struct isl_context_lex *clex = (struct isl_context_lex *)context;
2201         if (!clex->tab)
2202                 return NULL;
2203         return isl_tab_peek_bset(clex->tab);
2204 }
2205
2206 static struct isl_tab *context_lex_peek_tab(struct isl_context *context)
2207 {
2208         struct isl_context_lex *clex = (struct isl_context_lex *)context;
2209         return clex->tab;
2210 }
2211
2212 static void context_lex_add_eq(struct isl_context *context, isl_int *eq,
2213                 int check, int update)
2214 {
2215         struct isl_context_lex *clex = (struct isl_context_lex *)context;
2216         if (isl_tab_extend_cons(clex->tab, 2) < 0)
2217                 goto error;
2218         if (add_lexmin_eq(clex->tab, eq) < 0)
2219                 goto error;
2220         if (check) {
2221                 int v = tab_has_valid_sample(clex->tab, eq, 1);
2222                 if (v < 0)
2223                         goto error;
2224                 if (!v)
2225                         clex->tab = check_integer_feasible(clex->tab);
2226         }
2227         if (update)
2228                 clex->tab = check_samples(clex->tab, eq, 1);
2229         return;
2230 error:
2231         isl_tab_free(clex->tab);
2232         clex->tab = NULL;
2233 }
2234
2235 static void context_lex_add_ineq(struct isl_context *context, isl_int *ineq,
2236                 int check, int update)
2237 {
2238         struct isl_context_lex *clex = (struct isl_context_lex *)context;
2239         if (isl_tab_extend_cons(clex->tab, 1) < 0)
2240                 goto error;
2241         clex->tab = add_lexmin_ineq(clex->tab, ineq);
2242         if (check) {
2243                 int v = tab_has_valid_sample(clex->tab, ineq, 0);
2244                 if (v < 0)
2245                         goto error;
2246                 if (!v)
2247                         clex->tab = check_integer_feasible(clex->tab);
2248         }
2249         if (update)
2250                 clex->tab = check_samples(clex->tab, ineq, 0);
2251         return;
2252 error:
2253         isl_tab_free(clex->tab);
2254         clex->tab = NULL;
2255 }
2256
2257 static int context_lex_add_ineq_wrap(void *user, isl_int *ineq)
2258 {
2259         struct isl_context *context = (struct isl_context *)user;
2260         context_lex_add_ineq(context, ineq, 0, 0);
2261         return context->op->is_ok(context) ? 0 : -1;
2262 }
2263
2264 /* Check which signs can be obtained by "ineq" on all the currently
2265  * active sample values.  See row_sign for more information.
2266  */
2267 static enum isl_tab_row_sign tab_ineq_sign(struct isl_tab *tab, isl_int *ineq,
2268         int strict)
2269 {
2270         int i;
2271         int sgn;
2272         isl_int tmp;
2273         enum isl_tab_row_sign res = isl_tab_row_unknown;
2274
2275         isl_assert(tab->mat->ctx, tab->samples, return isl_tab_row_unknown);
2276         isl_assert(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var,
2277                         return isl_tab_row_unknown);
2278
2279         isl_int_init(tmp);
2280         for (i = tab->n_outside; i < tab->n_sample; ++i) {
2281                 isl_seq_inner_product(tab->samples->row[i], ineq,
2282                                         1 + tab->n_var, &tmp);
2283                 sgn = isl_int_sgn(tmp);
2284                 if (sgn > 0 || (sgn == 0 && strict)) {
2285                         if (res == isl_tab_row_unknown)
2286                                 res = isl_tab_row_pos;
2287                         if (res == isl_tab_row_neg)
2288                                 res = isl_tab_row_any;
2289                 }
2290                 if (sgn < 0) {
2291                         if (res == isl_tab_row_unknown)
2292                                 res = isl_tab_row_neg;
2293                         if (res == isl_tab_row_pos)
2294                                 res = isl_tab_row_any;
2295                 }
2296                 if (res == isl_tab_row_any)
2297                         break;
2298         }
2299         isl_int_clear(tmp);
2300
2301         return res;
2302 }
2303
2304 static enum isl_tab_row_sign context_lex_ineq_sign(struct isl_context *context,
2305                         isl_int *ineq, int strict)
2306 {
2307         struct isl_context_lex *clex = (struct isl_context_lex *)context;
2308         return tab_ineq_sign(clex->tab, ineq, strict);
2309 }
2310
2311 /* Check whether "ineq" can be added to the tableau without rendering
2312  * it infeasible.
2313  */
2314 static int context_lex_test_ineq(struct isl_context *context, isl_int *ineq)
2315 {
2316         struct isl_context_lex *clex = (struct isl_context_lex *)context;
2317         struct isl_tab_undo *snap;
2318         int feasible;
2319
2320         if (!clex->tab)
2321                 return -1;
2322
2323         if (isl_tab_extend_cons(clex->tab, 1) < 0)
2324                 return -1;
2325
2326         snap = isl_tab_snap(clex->tab);
2327         if (isl_tab_push_basis(clex->tab) < 0)
2328                 return -1;
2329         clex->tab = add_lexmin_ineq(clex->tab, ineq);
2330         clex->tab = check_integer_feasible(clex->tab);
2331         if (!clex->tab)
2332                 return -1;
2333         feasible = !clex->tab->empty;
2334         if (isl_tab_rollback(clex->tab, snap) < 0)
2335                 return -1;
2336
2337         return feasible;
2338 }
2339
2340 static int context_lex_get_div(struct isl_context *context, struct isl_tab *tab,
2341                 struct isl_vec *div)
2342 {
2343         return get_div(tab, context, div);
2344 }
2345
2346 /* Add a div specified by "div" to the context tableau and return
2347  * 1 if the div is obviously non-negative.
2348  * context_tab_add_div will always return 1, because all variables
2349  * in a isl_context_lex tableau are non-negative.
2350  * However, if we are using a big parameter in the context, then this only
2351  * reflects the non-negativity of the variable used to _encode_ the
2352  * div, i.e., div' = M + div, so we can't draw any conclusions.
2353  */
2354 static int context_lex_add_div(struct isl_context *context, struct isl_vec *div)
2355 {
2356         struct isl_context_lex *clex = (struct isl_context_lex *)context;
2357         int nonneg;
2358         nonneg = context_tab_add_div(clex->tab, div,
2359                                         context_lex_add_ineq_wrap, context);
2360         if (nonneg < 0)
2361                 return -1;
2362         if (clex->tab->M)
2363                 return 0;
2364         return nonneg;
2365 }
2366
2367 static int context_lex_detect_equalities(struct isl_context *context,
2368                 struct isl_tab *tab)
2369 {
2370         return 0;
2371 }
2372
2373 static int context_lex_best_split(struct isl_context *context,
2374                 struct isl_tab *tab)
2375 {
2376         struct isl_context_lex *clex = (struct isl_context_lex *)context;
2377         struct isl_tab_undo *snap;
2378         int r;
2379
2380         snap = isl_tab_snap(clex->tab);
2381         if (isl_tab_push_basis(clex->tab) < 0)
2382                 return -1;
2383         r = best_split(tab, clex->tab);
2384
2385         if (r >= 0 && isl_tab_rollback(clex->tab, snap) < 0)
2386                 return -1;
2387
2388         return r;
2389 }
2390
2391 static int context_lex_is_empty(struct isl_context *context)
2392 {
2393         struct isl_context_lex *clex = (struct isl_context_lex *)context;
2394         if (!clex->tab)
2395                 return -1;
2396         return clex->tab->empty;
2397 }
2398
2399 static void *context_lex_save(struct isl_context *context)
2400 {
2401         struct isl_context_lex *clex = (struct isl_context_lex *)context;
2402         struct isl_tab_undo *snap;
2403
2404         snap = isl_tab_snap(clex->tab);
2405         if (isl_tab_push_basis(clex->tab) < 0)
2406                 return NULL;
2407         if (isl_tab_save_samples(clex->tab) < 0)
2408                 return NULL;
2409
2410         return snap;
2411 }
2412
2413 static void context_lex_restore(struct isl_context *context, void *save)
2414 {
2415         struct isl_context_lex *clex = (struct isl_context_lex *)context;
2416         if (isl_tab_rollback(clex->tab, (struct isl_tab_undo *)save) < 0) {
2417                 isl_tab_free(clex->tab);
2418                 clex->tab = NULL;
2419         }
2420 }
2421
2422 static int context_lex_is_ok(struct isl_context *context)
2423 {
2424         struct isl_context_lex *clex = (struct isl_context_lex *)context;
2425         return !!clex->tab;
2426 }
2427
2428 /* For each variable in the context tableau, check if the variable can
2429  * only attain non-negative values.  If so, mark the parameter as non-negative
2430  * in the main tableau.  This allows for a more direct identification of some
2431  * cases of violated constraints.
2432  */
2433 static struct isl_tab *tab_detect_nonnegative_parameters(struct isl_tab *tab,
2434         struct isl_tab *context_tab)
2435 {
2436         int i;
2437         struct isl_tab_undo *snap;
2438         struct isl_vec *ineq = NULL;
2439         struct isl_tab_var *var;
2440         int n;
2441
2442         if (context_tab->n_var == 0)
2443                 return tab;
2444
2445         ineq = isl_vec_alloc(tab->mat->ctx, 1 + context_tab->n_var);
2446         if (!ineq)
2447                 goto error;
2448
2449         if (isl_tab_extend_cons(context_tab, 1) < 0)
2450                 goto error;
2451
2452         snap = isl_tab_snap(context_tab);
2453
2454         n = 0;
2455         isl_seq_clr(ineq->el, ineq->size);
2456         for (i = 0; i < context_tab->n_var; ++i) {
2457                 isl_int_set_si(ineq->el[1 + i], 1);
2458                 if (isl_tab_add_ineq(context_tab, ineq->el) < 0)
2459                         goto error;
2460                 var = &context_tab->con[context_tab->n_con - 1];
2461                 if (!context_tab->empty &&
2462                     !isl_tab_min_at_most_neg_one(context_tab, var)) {
2463                         int j = i;
2464                         if (i >= tab->n_param)
2465                                 j = i - tab->n_param + tab->n_var - tab->n_div;
2466                         tab->var[j].is_nonneg = 1;
2467                         n++;
2468                 }
2469                 isl_int_set_si(ineq->el[1 + i], 0);
2470                 if (isl_tab_rollback(context_tab, snap) < 0)
2471                         goto error;
2472         }
2473
2474         if (context_tab->M && n == context_tab->n_var) {
2475                 context_tab->mat = isl_mat_drop_cols(context_tab->mat, 2, 1);
2476                 context_tab->M = 0;
2477         }
2478
2479         isl_vec_free(ineq);
2480         return tab;
2481 error:
2482         isl_vec_free(ineq);
2483         isl_tab_free(tab);
2484         return NULL;
2485 }
2486
2487 static struct isl_tab *context_lex_detect_nonnegative_parameters(
2488         struct isl_context *context, struct isl_tab *tab)
2489 {
2490         struct isl_context_lex *clex = (struct isl_context_lex *)context;
2491         struct isl_tab_undo *snap;
2492
2493         if (!tab)
2494                 return NULL;
2495
2496         snap = isl_tab_snap(clex->tab);
2497         if (isl_tab_push_basis(clex->tab) < 0)
2498                 goto error;
2499
2500         tab = tab_detect_nonnegative_parameters(tab, clex->tab);
2501
2502         if (isl_tab_rollback(clex->tab, snap) < 0)
2503                 goto error;
2504
2505         return tab;
2506 error:
2507         isl_tab_free(tab);
2508         return NULL;
2509 }
2510
2511 static void context_lex_invalidate(struct isl_context *context)
2512 {
2513         struct isl_context_lex *clex = (struct isl_context_lex *)context;
2514         isl_tab_free(clex->tab);
2515         clex->tab = NULL;
2516 }
2517
2518 static void context_lex_free(struct isl_context *context)
2519 {
2520         struct isl_context_lex *clex = (struct isl_context_lex *)context;
2521         isl_tab_free(clex->tab);
2522         free(clex);
2523 }
2524
2525 struct isl_context_op isl_context_lex_op = {
2526         context_lex_detect_nonnegative_parameters,
2527         context_lex_peek_basic_set,
2528         context_lex_peek_tab,
2529         context_lex_add_eq,
2530         context_lex_add_ineq,
2531         context_lex_ineq_sign,
2532         context_lex_test_ineq,
2533         context_lex_get_div,
2534         context_lex_add_div,
2535         context_lex_detect_equalities,
2536         context_lex_best_split,
2537         context_lex_is_empty,
2538         context_lex_is_ok,
2539         context_lex_save,
2540         context_lex_restore,
2541         context_lex_invalidate,
2542         context_lex_free,
2543 };
2544
2545 static struct isl_tab *context_tab_for_lexmin(struct isl_basic_set *bset)
2546 {
2547         struct isl_tab *tab;
2548
2549         if (!bset)
2550                 return NULL;
2551         tab = tab_for_lexmin((struct isl_basic_map *)bset, NULL, 1, 0);
2552         if (!tab)
2553                 goto error;
2554         if (isl_tab_track_bset(tab, bset) < 0)
2555                 goto error;
2556         tab = isl_tab_init_samples(tab);
2557         return tab;
2558 error:
2559         isl_basic_set_free(bset);
2560         return NULL;
2561 }
2562
2563 static struct isl_context *isl_context_lex_alloc(struct isl_basic_set *dom)
2564 {
2565         struct isl_context_lex *clex;
2566
2567         if (!dom)
2568                 return NULL;
2569
2570         clex = isl_alloc_type(dom->ctx, struct isl_context_lex);
2571         if (!clex)
2572                 return NULL;
2573
2574         clex->context.op = &isl_context_lex_op;
2575
2576         clex->tab = context_tab_for_lexmin(isl_basic_set_copy(dom));
2577         if (restore_lexmin(clex->tab) < 0)
2578                 goto error;
2579         clex->tab = check_integer_feasible(clex->tab);
2580         if (!clex->tab)
2581                 goto error;
2582
2583         return &clex->context;
2584 error:
2585         clex->context.op->free(&clex->context);
2586         return NULL;
2587 }
2588
2589 struct isl_context_gbr {
2590         struct isl_context context;
2591         struct isl_tab *tab;
2592         struct isl_tab *shifted;
2593         struct isl_tab *cone;
2594 };
2595
2596 static struct isl_tab *context_gbr_detect_nonnegative_parameters(
2597         struct isl_context *context, struct isl_tab *tab)
2598 {
2599         struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2600         if (!tab)
2601                 return NULL;
2602         return tab_detect_nonnegative_parameters(tab, cgbr->tab);
2603 }
2604
2605 static struct isl_basic_set *context_gbr_peek_basic_set(
2606         struct isl_context *context)
2607 {
2608         struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2609         if (!cgbr->tab)
2610                 return NULL;
2611         return isl_tab_peek_bset(cgbr->tab);
2612 }
2613
2614 static struct isl_tab *context_gbr_peek_tab(struct isl_context *context)
2615 {
2616         struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2617         return cgbr->tab;
2618 }
2619
2620 /* Initialize the "shifted" tableau of the context, which
2621  * contains the constraints of the original tableau shifted
2622  * by the sum of all negative coefficients.  This ensures
2623  * that any rational point in the shifted tableau can
2624  * be rounded up to yield an integer point in the original tableau.
2625  */
2626 static void gbr_init_shifted(struct isl_context_gbr *cgbr)
2627 {
2628         int i, j;
2629         struct isl_vec *cst;
2630         struct isl_basic_set *bset = isl_tab_peek_bset(cgbr->tab);
2631         unsigned dim = isl_basic_set_total_dim(bset);
2632
2633         cst = isl_vec_alloc(cgbr->tab->mat->ctx, bset->n_ineq);
2634         if (!cst)
2635                 return;
2636
2637         for (i = 0; i < bset->n_ineq; ++i) {
2638                 isl_int_set(cst->el[i], bset->ineq[i][0]);
2639                 for (j = 0; j < dim; ++j) {
2640                         if (!isl_int_is_neg(bset->ineq[i][1 + j]))
2641                                 continue;
2642                         isl_int_add(bset->ineq[i][0], bset->ineq[i][0],
2643                                     bset->ineq[i][1 + j]);
2644                 }
2645         }
2646
2647         cgbr->shifted = isl_tab_from_basic_set(bset, 0);
2648
2649         for (i = 0; i < bset->n_ineq; ++i)
2650                 isl_int_set(bset->ineq[i][0], cst->el[i]);
2651
2652         isl_vec_free(cst);
2653 }
2654
2655 /* Check if the shifted tableau is non-empty, and if so
2656  * use the sample point to construct an integer point
2657  * of the context tableau.
2658  */
2659 static struct isl_vec *gbr_get_shifted_sample(struct isl_context_gbr *cgbr)
2660 {
2661         struct isl_vec *sample;
2662
2663         if (!cgbr->shifted)
2664                 gbr_init_shifted(cgbr);
2665         if (!cgbr->shifted)
2666                 return NULL;
2667         if (cgbr->shifted->empty)
2668                 return isl_vec_alloc(cgbr->tab->mat->ctx, 0);
2669
2670         sample = isl_tab_get_sample_value(cgbr->shifted);
2671         sample = isl_vec_ceil(sample);
2672
2673         return sample;
2674 }
2675
2676 static struct isl_basic_set *drop_constant_terms(struct isl_basic_set *bset)
2677 {
2678         int i;
2679
2680         if (!bset)
2681                 return NULL;
2682
2683         for (i = 0; i < bset->n_eq; ++i)
2684                 isl_int_set_si(bset->eq[i][0], 0);
2685
2686         for (i = 0; i < bset->n_ineq; ++i)
2687                 isl_int_set_si(bset->ineq[i][0], 0);
2688
2689         return bset;
2690 }
2691
2692 static int use_shifted(struct isl_context_gbr *cgbr)
2693 {
2694         return cgbr->tab->bmap->n_eq == 0 && cgbr->tab->bmap->n_div == 0;
2695 }
2696
2697 static struct isl_vec *gbr_get_sample(struct isl_context_gbr *cgbr)
2698 {
2699         struct isl_basic_set *bset;
2700         struct isl_basic_set *cone;
2701
2702         if (isl_tab_sample_is_integer(cgbr->tab))
2703                 return isl_tab_get_sample_value(cgbr->tab);
2704
2705         if (use_shifted(cgbr)) {
2706                 struct isl_vec *sample;
2707
2708                 sample = gbr_get_shifted_sample(cgbr);
2709                 if (!sample || sample->size > 0)
2710                         return sample;
2711
2712                 isl_vec_free(sample);
2713         }
2714
2715         if (!cgbr->cone) {
2716                 bset = isl_tab_peek_bset(cgbr->tab);
2717                 cgbr->cone = isl_tab_from_recession_cone(bset, 0);
2718                 if (!cgbr->cone)
2719                         return NULL;
2720                 if (isl_tab_track_bset(cgbr->cone,
2721                                         isl_basic_set_copy(bset)) < 0)
2722                         return NULL;
2723         }
2724         if (isl_tab_detect_implicit_equalities(cgbr->cone) < 0)
2725                 return NULL;
2726
2727         if (cgbr->cone->n_dead == cgbr->cone->n_col) {
2728                 struct isl_vec *sample;
2729                 struct isl_tab_undo *snap;
2730
2731                 if (cgbr->tab->basis) {
2732                         if (cgbr->tab->basis->n_col != 1 + cgbr->tab->n_var) {
2733                                 isl_mat_free(cgbr->tab->basis);
2734                                 cgbr->tab->basis = NULL;
2735                         }
2736                         cgbr->tab->n_zero = 0;
2737                         cgbr->tab->n_unbounded = 0;
2738                 }
2739
2740                 snap = isl_tab_snap(cgbr->tab);
2741
2742                 sample = isl_tab_sample(cgbr->tab);
2743
2744                 if (isl_tab_rollback(cgbr->tab, snap) < 0) {
2745                         isl_vec_free(sample);
2746                         return NULL;
2747                 }
2748
2749                 return sample;
2750         }
2751
2752         cone = isl_basic_set_dup(isl_tab_peek_bset(cgbr->cone));
2753         cone = drop_constant_terms(cone);
2754         cone = isl_basic_set_update_from_tab(cone, cgbr->cone);
2755         cone = isl_basic_set_underlying_set(cone);
2756         cone = isl_basic_set_gauss(cone, NULL);
2757
2758         bset = isl_basic_set_dup(isl_tab_peek_bset(cgbr->tab));
2759         bset = isl_basic_set_update_from_tab(bset, cgbr->tab);
2760         bset = isl_basic_set_underlying_set(bset);
2761         bset = isl_basic_set_gauss(bset, NULL);
2762
2763         return isl_basic_set_sample_with_cone(bset, cone);
2764 }
2765
2766 static void check_gbr_integer_feasible(struct isl_context_gbr *cgbr)
2767 {
2768         struct isl_vec *sample;
2769
2770         if (!cgbr->tab)
2771                 return;
2772
2773         if (cgbr->tab->empty)
2774                 return;
2775
2776         sample = gbr_get_sample(cgbr);
2777         if (!sample)
2778                 goto error;
2779
2780         if (sample->size == 0) {
2781                 isl_vec_free(sample);
2782                 if (isl_tab_mark_empty(cgbr->tab) < 0)
2783                         goto error;
2784                 return;
2785         }
2786
2787         cgbr->tab = isl_tab_add_sample(cgbr->tab, sample);
2788
2789         return;
2790 error:
2791         isl_tab_free(cgbr->tab);
2792         cgbr->tab = NULL;
2793 }
2794
2795 static struct isl_tab *add_gbr_eq(struct isl_tab *tab, isl_int *eq)
2796 {
2797         if (!tab)
2798                 return NULL;
2799
2800         if (isl_tab_extend_cons(tab, 2) < 0)
2801                 goto error;
2802
2803         if (isl_tab_add_eq(tab, eq) < 0)
2804                 goto error;
2805
2806         return tab;
2807 error:
2808         isl_tab_free(tab);
2809         return NULL;
2810 }
2811
2812 static void context_gbr_add_eq(struct isl_context *context, isl_int *eq,
2813                 int check, int update)
2814 {
2815         struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2816
2817         cgbr->tab = add_gbr_eq(cgbr->tab, eq);
2818
2819         if (cgbr->cone && cgbr->cone->n_col != cgbr->cone->n_dead) {
2820                 if (isl_tab_extend_cons(cgbr->cone, 2) < 0)
2821                         goto error;
2822                 if (isl_tab_add_eq(cgbr->cone, eq) < 0)
2823                         goto error;
2824         }
2825
2826         if (check) {
2827                 int v = tab_has_valid_sample(cgbr->tab, eq, 1);
2828                 if (v < 0)
2829                         goto error;
2830                 if (!v)
2831                         check_gbr_integer_feasible(cgbr);
2832         }
2833         if (update)
2834                 cgbr->tab = check_samples(cgbr->tab, eq, 1);
2835         return;
2836 error:
2837         isl_tab_free(cgbr->tab);
2838         cgbr->tab = NULL;
2839 }
2840
2841 static void add_gbr_ineq(struct isl_context_gbr *cgbr, isl_int *ineq)
2842 {
2843         if (!cgbr->tab)
2844                 return;
2845
2846         if (isl_tab_extend_cons(cgbr->tab, 1) < 0)
2847                 goto error;
2848
2849         if (isl_tab_add_ineq(cgbr->tab, ineq) < 0)
2850                 goto error;
2851
2852         if (cgbr->shifted && !cgbr->shifted->empty && use_shifted(cgbr)) {
2853                 int i;
2854                 unsigned dim;
2855                 dim = isl_basic_map_total_dim(cgbr->tab->bmap);
2856
2857                 if (isl_tab_extend_cons(cgbr->shifted, 1) < 0)
2858                         goto error;
2859
2860                 for (i = 0; i < dim; ++i) {
2861                         if (!isl_int_is_neg(ineq[1 + i]))
2862                                 continue;
2863                         isl_int_add(ineq[0], ineq[0], ineq[1 + i]);
2864                 }
2865
2866                 if (isl_tab_add_ineq(cgbr->shifted, ineq) < 0)
2867                         goto error;
2868
2869                 for (i = 0; i < dim; ++i) {
2870                         if (!isl_int_is_neg(ineq[1 + i]))
2871                                 continue;
2872                         isl_int_sub(ineq[0], ineq[0], ineq[1 + i]);
2873                 }
2874         }
2875
2876         if (cgbr->cone && cgbr->cone->n_col != cgbr->cone->n_dead) {
2877                 if (isl_tab_extend_cons(cgbr->cone, 1) < 0)
2878                         goto error;
2879                 if (isl_tab_add_ineq(cgbr->cone, ineq) < 0)
2880                         goto error;
2881         }
2882
2883         return;
2884 error:
2885         isl_tab_free(cgbr->tab);
2886         cgbr->tab = NULL;
2887 }
2888
2889 static void context_gbr_add_ineq(struct isl_context *context, isl_int *ineq,
2890                 int check, int update)
2891 {
2892         struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2893
2894         add_gbr_ineq(cgbr, ineq);
2895         if (!cgbr->tab)
2896                 return;
2897
2898         if (check) {
2899                 int v = tab_has_valid_sample(cgbr->tab, ineq, 0);
2900                 if (v < 0)
2901                         goto error;
2902                 if (!v)
2903                         check_gbr_integer_feasible(cgbr);
2904         }
2905         if (update)
2906                 cgbr->tab = check_samples(cgbr->tab, ineq, 0);
2907         return;
2908 error:
2909         isl_tab_free(cgbr->tab);
2910         cgbr->tab = NULL;
2911 }
2912
2913 static int context_gbr_add_ineq_wrap(void *user, isl_int *ineq)
2914 {
2915         struct isl_context *context = (struct isl_context *)user;
2916         context_gbr_add_ineq(context, ineq, 0, 0);
2917         return context->op->is_ok(context) ? 0 : -1;
2918 }
2919
2920 static enum isl_tab_row_sign context_gbr_ineq_sign(struct isl_context *context,
2921                         isl_int *ineq, int strict)
2922 {
2923         struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2924         return tab_ineq_sign(cgbr->tab, ineq, strict);
2925 }
2926
2927 /* Check whether "ineq" can be added to the tableau without rendering
2928  * it infeasible.
2929  */
2930 static int context_gbr_test_ineq(struct isl_context *context, isl_int *ineq)
2931 {
2932         struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2933         struct isl_tab_undo *snap;
2934         struct isl_tab_undo *shifted_snap = NULL;
2935         struct isl_tab_undo *cone_snap = NULL;
2936         int feasible;
2937
2938         if (!cgbr->tab)
2939                 return -1;
2940
2941         if (isl_tab_extend_cons(cgbr->tab, 1) < 0)
2942                 return -1;
2943
2944         snap = isl_tab_snap(cgbr->tab);
2945         if (cgbr->shifted)
2946                 shifted_snap = isl_tab_snap(cgbr->shifted);
2947         if (cgbr->cone)
2948                 cone_snap = isl_tab_snap(cgbr->cone);
2949         add_gbr_ineq(cgbr, ineq);
2950         check_gbr_integer_feasible(cgbr);
2951         if (!cgbr->tab)
2952                 return -1;
2953         feasible = !cgbr->tab->empty;
2954         if (isl_tab_rollback(cgbr->tab, snap) < 0)
2955                 return -1;
2956         if (shifted_snap) {
2957                 if (isl_tab_rollback(cgbr->shifted, shifted_snap))
2958                         return -1;
2959         } else if (cgbr->shifted) {
2960                 isl_tab_free(cgbr->shifted);
2961                 cgbr->shifted = NULL;
2962         }
2963         if (cone_snap) {
2964                 if (isl_tab_rollback(cgbr->cone, cone_snap))
2965                         return -1;
2966         } else if (cgbr->cone) {
2967                 isl_tab_free(cgbr->cone);
2968                 cgbr->cone = NULL;
2969         }
2970
2971         return feasible;
2972 }
2973
2974 /* Return the column of the last of the variables associated to
2975  * a column that has a non-zero coefficient.
2976  * This function is called in a context where only coefficients
2977  * of parameters or divs can be non-zero.
2978  */
2979 static int last_non_zero_var_col(struct isl_tab *tab, isl_int *p)
2980 {
2981         int i;
2982         int col;
2983
2984         if (tab->n_var == 0)
2985                 return -1;
2986
2987         for (i = tab->n_var - 1; i >= 0; --i) {
2988                 if (i >= tab->n_param && i < tab->n_var - tab->n_div)
2989                         continue;
2990                 if (tab->var[i].is_row)
2991                         continue;
2992                 col = tab->var[i].index;
2993                 if (!isl_int_is_zero(p[col]))
2994                         return col;
2995         }
2996
2997         return -1;
2998 }
2999
3000 /* Look through all the recently added equalities in the context
3001  * to see if we can propagate any of them to the main tableau.
3002  *
3003  * The newly added equalities in the context are encoded as pairs
3004  * of inequalities starting at inequality "first".
3005  *
3006  * We tentatively add each of these equalities to the main tableau
3007  * and if this happens to result in a row with a final coefficient
3008  * that is one or negative one, we use it to kill a column
3009  * in the main tableau.  Otherwise, we discard the tentatively
3010  * added row.
3011  */
3012 static void propagate_equalities(struct isl_context_gbr *cgbr,
3013         struct isl_tab *tab, unsigned first)
3014 {
3015         int i;
3016         struct isl_vec *eq = NULL;
3017
3018         eq = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
3019         if (!eq)
3020                 goto error;
3021
3022         if (isl_tab_extend_cons(tab, (cgbr->tab->bmap->n_ineq - first)/2) < 0)
3023                 goto error;
3024
3025         isl_seq_clr(eq->el + 1 + tab->n_param,
3026                     tab->n_var - tab->n_param - tab->n_div);
3027         for (i = first; i < cgbr->tab->bmap->n_ineq; i += 2) {
3028                 int j;
3029                 int r;
3030                 struct isl_tab_undo *snap;
3031                 snap = isl_tab_snap(tab);
3032
3033                 isl_seq_cpy(eq->el, cgbr->tab->bmap->ineq[i], 1 + tab->n_param);
3034                 isl_seq_cpy(eq->el + 1 + tab->n_var - tab->n_div,
3035                             cgbr->tab->bmap->ineq[i] + 1 + tab->n_param,
3036                             tab->n_div);
3037
3038                 r = isl_tab_add_row(tab, eq->el);
3039                 if (r < 0)
3040                         goto error;
3041                 r = tab->con[r].index;
3042                 j = last_non_zero_var_col(tab, tab->mat->row[r] + 2 + tab->M);
3043                 if (j < 0 || j < tab->n_dead ||
3044                     !isl_int_is_one(tab->mat->row[r][0]) ||
3045                     (!isl_int_is_one(tab->mat->row[r][2 + tab->M + j]) &&
3046                      !isl_int_is_negone(tab->mat->row[r][2 + tab->M + j]))) {
3047                         if (isl_tab_rollback(tab, snap) < 0)
3048                                 goto error;
3049                         continue;
3050                 }
3051                 if (isl_tab_pivot(tab, r, j) < 0)
3052                         goto error;
3053                 if (isl_tab_kill_col(tab, j) < 0)
3054                         goto error;
3055
3056                 if (restore_lexmin(tab) < 0)
3057                         goto error;
3058         }
3059
3060         isl_vec_free(eq);
3061
3062         return;
3063 error:
3064         isl_vec_free(eq);
3065         isl_tab_free(cgbr->tab);
3066         cgbr->tab = NULL;
3067 }
3068
3069 static int context_gbr_detect_equalities(struct isl_context *context,
3070         struct isl_tab *tab)
3071 {
3072         struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3073         struct isl_ctx *ctx;
3074         unsigned n_ineq;
3075
3076         ctx = cgbr->tab->mat->ctx;
3077
3078         if (!cgbr->cone) {
3079                 struct isl_basic_set *bset = isl_tab_peek_bset(cgbr->tab);
3080                 cgbr->cone = isl_tab_from_recession_cone(bset, 0);
3081                 if (!cgbr->cone)
3082                         goto error;
3083                 if (isl_tab_track_bset(cgbr->cone,
3084                                         isl_basic_set_copy(bset)) < 0)
3085                         goto error;
3086         }
3087         if (isl_tab_detect_implicit_equalities(cgbr->cone) < 0)
3088                 goto error;
3089
3090         n_ineq = cgbr->tab->bmap->n_ineq;
3091         cgbr->tab = isl_tab_detect_equalities(cgbr->tab, cgbr->cone);
3092         if (cgbr->tab && cgbr->tab->bmap->n_ineq > n_ineq)
3093                 propagate_equalities(cgbr, tab, n_ineq);
3094
3095         return 0;
3096 error:
3097         isl_tab_free(cgbr->tab);
3098         cgbr->tab = NULL;
3099         return -1;
3100 }
3101
3102 static int context_gbr_get_div(struct isl_context *context, struct isl_tab *tab,
3103                 struct isl_vec *div)
3104 {
3105         return get_div(tab, context, div);
3106 }
3107
3108 static int context_gbr_add_div(struct isl_context *context, struct isl_vec *div)
3109 {
3110         struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3111         if (cgbr->cone) {
3112                 int k;
3113
3114                 if (isl_tab_extend_cons(cgbr->cone, 3) < 0)
3115                         return -1;
3116                 if (isl_tab_extend_vars(cgbr->cone, 1) < 0)
3117                         return -1;
3118                 if (isl_tab_allocate_var(cgbr->cone) <0)
3119                         return -1;
3120
3121                 cgbr->cone->bmap = isl_basic_map_extend_space(cgbr->cone->bmap,
3122                         isl_basic_map_get_space(cgbr->cone->bmap), 1, 0, 2);
3123                 k = isl_basic_map_alloc_div(cgbr->cone->bmap);
3124                 if (k < 0)
3125                         return -1;
3126                 isl_seq_cpy(cgbr->cone->bmap->div[k], div->el, div->size);
3127                 if (isl_tab_push(cgbr->cone, isl_tab_undo_bmap_div) < 0)
3128                         return -1;
3129         }
3130         return context_tab_add_div(cgbr->tab, div,
3131                                         context_gbr_add_ineq_wrap, context);
3132 }
3133
3134 static int context_gbr_best_split(struct isl_context *context,
3135                 struct isl_tab *tab)
3136 {
3137         struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3138         struct isl_tab_undo *snap;
3139         int r;
3140
3141         snap = isl_tab_snap(cgbr->tab);
3142         r = best_split(tab, cgbr->tab);
3143
3144         if (r >= 0 && isl_tab_rollback(cgbr->tab, snap) < 0)
3145                 return -1;
3146
3147         return r;
3148 }
3149
3150 static int context_gbr_is_empty(struct isl_context *context)
3151 {
3152         struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3153         if (!cgbr->tab)
3154                 return -1;
3155         return cgbr->tab->empty;
3156 }
3157
3158 struct isl_gbr_tab_undo {
3159         struct isl_tab_undo *tab_snap;
3160         struct isl_tab_undo *shifted_snap;
3161         struct isl_tab_undo *cone_snap;
3162 };
3163
3164 static void *context_gbr_save(struct isl_context *context)
3165 {
3166         struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3167         struct isl_gbr_tab_undo *snap;
3168
3169         snap = isl_alloc_type(cgbr->tab->mat->ctx, struct isl_gbr_tab_undo);
3170         if (!snap)
3171                 return NULL;
3172
3173         snap->tab_snap = isl_tab_snap(cgbr->tab);
3174         if (isl_tab_save_samples(cgbr->tab) < 0)
3175                 goto error;
3176
3177         if (cgbr->shifted)
3178                 snap->shifted_snap = isl_tab_snap(cgbr->shifted);
3179         else
3180                 snap->shifted_snap = NULL;
3181
3182         if (cgbr->cone)
3183                 snap->cone_snap = isl_tab_snap(cgbr->cone);
3184         else
3185                 snap->cone_snap = NULL;
3186
3187         return snap;
3188 error:
3189         free(snap);
3190         return NULL;
3191 }
3192
3193 static void context_gbr_restore(struct isl_context *context, void *save)
3194 {
3195         struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3196         struct isl_gbr_tab_undo *snap = (struct isl_gbr_tab_undo *)save;
3197         if (!snap)
3198                 goto error;
3199         if (isl_tab_rollback(cgbr->tab, snap->tab_snap) < 0) {
3200                 isl_tab_free(cgbr->tab);
3201                 cgbr->tab = NULL;
3202         }
3203
3204         if (snap->shifted_snap) {
3205                 if (isl_tab_rollback(cgbr->shifted, snap->shifted_snap) < 0)
3206                         goto error;
3207         } else if (cgbr->shifted) {
3208                 isl_tab_free(cgbr->shifted);
3209                 cgbr->shifted = NULL;
3210         }
3211
3212         if (snap->cone_snap) {
3213                 if (isl_tab_rollback(cgbr->cone, snap->cone_snap) < 0)
3214                         goto error;
3215         } else if (cgbr->cone) {
3216                 isl_tab_free(cgbr->cone);
3217                 cgbr->cone = NULL;
3218         }
3219
3220         free(snap);
3221
3222         return;
3223 error:
3224         free(snap);
3225         isl_tab_free(cgbr->tab);
3226         cgbr->tab = NULL;
3227 }
3228
3229 static int context_gbr_is_ok(struct isl_context *context)
3230 {
3231         struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3232         return !!cgbr->tab;
3233 }
3234
3235 static void context_gbr_invalidate(struct isl_context *context)
3236 {
3237         struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3238         isl_tab_free(cgbr->tab);
3239         cgbr->tab = NULL;
3240 }
3241
3242 static void context_gbr_free(struct isl_context *context)
3243 {
3244         struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3245         isl_tab_free(cgbr->tab);
3246         isl_tab_free(cgbr->shifted);
3247         isl_tab_free(cgbr->cone);
3248         free(cgbr);
3249 }
3250
3251 struct isl_context_op isl_context_gbr_op = {
3252         context_gbr_detect_nonnegative_parameters,
3253         context_gbr_peek_basic_set,
3254         context_gbr_peek_tab,
3255         context_gbr_add_eq,
3256         context_gbr_add_ineq,
3257         context_gbr_ineq_sign,
3258         context_gbr_test_ineq,
3259         context_gbr_get_div,
3260         context_gbr_add_div,
3261         context_gbr_detect_equalities,
3262         context_gbr_best_split,
3263         context_gbr_is_empty,
3264         context_gbr_is_ok,
3265         context_gbr_save,
3266         context_gbr_restore,
3267         context_gbr_invalidate,
3268         context_gbr_free,
3269 };
3270
3271 static struct isl_context *isl_context_gbr_alloc(struct isl_basic_set *dom)
3272 {
3273         struct isl_context_gbr *cgbr;
3274
3275         if (!dom)
3276                 return NULL;
3277
3278         cgbr = isl_calloc_type(dom->ctx, struct isl_context_gbr);
3279         if (!cgbr)
3280                 return NULL;
3281
3282         cgbr->context.op = &isl_context_gbr_op;
3283
3284         cgbr->shifted = NULL;
3285         cgbr->cone = NULL;
3286         cgbr->tab = isl_tab_from_basic_set(dom, 1);
3287         cgbr->tab = isl_tab_init_samples(cgbr->tab);
3288         if (!cgbr->tab)
3289                 goto error;
3290         check_gbr_integer_feasible(cgbr);
3291
3292         return &cgbr->context;
3293 error:
3294         cgbr->context.op->free(&cgbr->context);
3295         return NULL;
3296 }
3297
3298 static struct isl_context *isl_context_alloc(struct isl_basic_set *dom)
3299 {
3300         if (!dom)
3301                 return NULL;
3302
3303         if (dom->ctx->opt->context == ISL_CONTEXT_LEXMIN)
3304                 return isl_context_lex_alloc(dom);
3305         else
3306                 return isl_context_gbr_alloc(dom);
3307 }
3308
3309 /* Construct an isl_sol_map structure for accumulating the solution.
3310  * If track_empty is set, then we also keep track of the parts
3311  * of the context where there is no solution.
3312  * If max is set, then we are solving a maximization, rather than
3313  * a minimization problem, which means that the variables in the
3314  * tableau have value "M - x" rather than "M + x".
3315  */
3316 static struct isl_sol *sol_map_init(struct isl_basic_map *bmap,
3317         struct isl_basic_set *dom, int track_empty, int max)
3318 {
3319         struct isl_sol_map *sol_map = NULL;
3320
3321         if (!bmap)
3322                 goto error;
3323
3324         sol_map = isl_calloc_type(bmap->ctx, struct isl_sol_map);
3325         if (!sol_map)
3326                 goto error;
3327
3328         sol_map->sol.rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
3329         sol_map->sol.dec_level.callback.run = &sol_dec_level_wrap;
3330         sol_map->sol.dec_level.sol = &sol_map->sol;
3331         sol_map->sol.max = max;
3332         sol_map->sol.n_out = isl_basic_map_dim(bmap, isl_dim_out);
3333         sol_map->sol.add = &sol_map_add_wrap;
3334         sol_map->sol.add_empty = track_empty ? &sol_map_add_empty_wrap : NULL;
3335         sol_map->sol.free = &sol_map_free_wrap;
3336         sol_map->map = isl_map_alloc_space(isl_basic_map_get_space(bmap), 1,
3337                                             ISL_MAP_DISJOINT);
3338         if (!sol_map->map)
3339                 goto error;
3340
3341         sol_map->sol.context = isl_context_alloc(dom);
3342         if (!sol_map->sol.context)
3343                 goto error;
3344
3345         if (track_empty) {
3346                 sol_map->empty = isl_set_alloc_space(isl_basic_set_get_space(dom),
3347                                                         1, ISL_SET_DISJOINT);
3348                 if (!sol_map->empty)
3349                         goto error;
3350         }
3351
3352         isl_basic_set_free(dom);
3353         return &sol_map->sol;
3354 error:
3355         isl_basic_set_free(dom);
3356         sol_map_free(sol_map);
3357         return NULL;
3358 }
3359
3360 /* Check whether all coefficients of (non-parameter) variables
3361  * are non-positive, meaning that no pivots can be performed on the row.
3362  */
3363 static int is_critical(struct isl_tab *tab, int row)
3364 {
3365         int j;
3366         unsigned off = 2 + tab->M;
3367
3368         for (j = tab->n_dead; j < tab->n_col; ++j) {
3369                 if (tab->col_var[j] >= 0 &&
3370                     (tab->col_var[j] < tab->n_param  ||
3371                     tab->col_var[j] >= tab->n_var - tab->n_div))
3372                         continue;
3373
3374                 if (isl_int_is_pos(tab->mat->row[row][off + j]))
3375                         return 0;
3376         }
3377
3378         return 1;
3379 }
3380
3381 /* Check whether the inequality represented by vec is strict over the integers,
3382  * i.e., there are no integer values satisfying the constraint with
3383  * equality.  This happens if the gcd of the coefficients is not a divisor
3384  * of the constant term.  If so, scale the constraint down by the gcd
3385  * of the coefficients.
3386  */
3387 static int is_strict(struct isl_vec *vec)
3388 {
3389         isl_int gcd;
3390         int strict = 0;
3391
3392         isl_int_init(gcd);
3393         isl_seq_gcd(vec->el + 1, vec->size - 1, &gcd);
3394         if (!isl_int_is_one(gcd)) {
3395                 strict = !isl_int_is_divisible_by(vec->el[0], gcd);
3396                 isl_int_fdiv_q(vec->el[0], vec->el[0], gcd);
3397                 isl_seq_scale_down(vec->el + 1, vec->el + 1, gcd, vec->size-1);
3398         }
3399         isl_int_clear(gcd);
3400
3401         return strict;
3402 }
3403
3404 /* Determine the sign of the given row of the main tableau.
3405  * The result is one of
3406  *      isl_tab_row_pos: always non-negative; no pivot needed
3407  *      isl_tab_row_neg: always non-positive; pivot
3408  *      isl_tab_row_any: can be both positive and negative; split
3409  *
3410  * We first handle some simple cases
3411  *      - the row sign may be known already
3412  *      - the row may be obviously non-negative
3413  *      - the parametric constant may be equal to that of another row
3414  *        for which we know the sign.  This sign will be either "pos" or
3415  *        "any".  If it had been "neg" then we would have pivoted before.
3416  *
3417  * If none of these cases hold, we check the value of the row for each
3418  * of the currently active samples.  Based on the signs of these values
3419  * we make an initial determination of the sign of the row.
3420  *
3421  *      all zero                        ->      unk(nown)
3422  *      all non-negative                ->      pos
3423  *      all non-positive                ->      neg
3424  *      both negative and positive      ->      all
3425  *
3426  * If we end up with "all", we are done.
3427  * Otherwise, we perform a check for positive and/or negative
3428  * values as follows.
3429  *
3430  *      samples        neg             unk             pos
3431  *      <0 ?                        Y        N      Y        N
3432  *                                          pos    any      pos
3433  *      >0 ?         Y      N    Y     N
3434  *                  any    neg  any   neg
3435  *
3436  * There is no special sign for "zero", because we can usually treat zero
3437  * as either non-negative or non-positive, whatever works out best.
3438  * However, if the row is "critical", meaning that pivoting is impossible
3439  * then we don't want to limp zero with the non-positive case, because
3440  * then we we would lose the solution for those values of the parameters
3441  * where the value of the row is zero.  Instead, we treat 0 as non-negative
3442  * ensuring a split if the row can attain both zero and negative values.
3443  * The same happens when the original constraint was one that could not
3444  * be satisfied with equality by any integer values of the parameters.
3445  * In this case, we normalize the constraint, but then a value of zero
3446  * for the normalized constraint is actually a positive value for the
3447  * original constraint, so again we need to treat zero as non-negative.
3448  * In both these cases, we have the following decision tree instead:
3449  *
3450  *      all non-negative                ->      pos
3451  *      all negative                    ->      neg
3452  *      both negative and non-negative  ->      all
3453  *
3454  *      samples        neg                             pos
3455  *      <0 ?                                        Y        N
3456  *                                                 any      pos
3457  *      >=0 ?        Y      N
3458  *                  any    neg
3459  */
3460 static enum isl_tab_row_sign row_sign(struct isl_tab *tab,
3461         struct isl_sol *sol, int row)
3462 {
3463         struct isl_vec *ineq = NULL;
3464         enum isl_tab_row_sign res = isl_tab_row_unknown;
3465         int critical;
3466         int strict;
3467         int row2;
3468
3469         if (tab->row_sign[row] != isl_tab_row_unknown)
3470                 return tab->row_sign[row];
3471         if (is_obviously_nonneg(tab, row))
3472                 return isl_tab_row_pos;
3473         for (row2 = tab->n_redundant; row2 < tab->n_row; ++row2) {
3474                 if (tab->row_sign[row2] == isl_tab_row_unknown)
3475                         continue;
3476                 if (identical_parameter_line(tab, row, row2))
3477                         return tab->row_sign[row2];
3478         }
3479
3480         critical = is_critical(tab, row);
3481
3482         ineq = get_row_parameter_ineq(tab, row);
3483         if (!ineq)
3484                 goto error;
3485
3486         strict = is_strict(ineq);
3487
3488         res = sol->context->op->ineq_sign(sol->context, ineq->el,
3489                                           critical || strict);
3490
3491         if (res == isl_tab_row_unknown || res == isl_tab_row_pos) {
3492                 /* test for negative values */
3493                 int feasible;
3494                 isl_seq_neg(ineq->el, ineq->el, ineq->size);
3495                 isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3496
3497                 feasible = sol->context->op->test_ineq(sol->context, ineq->el);
3498                 if (feasible < 0)
3499                         goto error;
3500                 if (!feasible)
3501                         res = isl_tab_row_pos;
3502                 else
3503                         res = (res == isl_tab_row_unknown) ? isl_tab_row_neg
3504                                                            : isl_tab_row_any;
3505                 if (res == isl_tab_row_neg) {
3506                         isl_seq_neg(ineq->el, ineq->el, ineq->size);
3507                         isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3508                 }
3509         }
3510
3511         if (res == isl_tab_row_neg) {
3512                 /* test for positive values */
3513                 int feasible;
3514                 if (!critical && !strict)
3515                         isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3516
3517                 feasible = sol->context->op->test_ineq(sol->context, ineq->el);
3518                 if (feasible < 0)
3519                         goto error;
3520                 if (feasible)
3521                         res = isl_tab_row_any;
3522         }
3523
3524         isl_vec_free(ineq);
3525         return res;
3526 error:
3527         isl_vec_free(ineq);
3528         return isl_tab_row_unknown;
3529 }
3530
3531 static void find_solutions(struct isl_sol *sol, struct isl_tab *tab);
3532
3533 /* Find solutions for values of the parameters that satisfy the given
3534  * inequality.
3535  *
3536  * We currently take a snapshot of the context tableau that is reset
3537  * when we return from this function, while we make a copy of the main
3538  * tableau, leaving the original main tableau untouched.
3539  * These are fairly arbitrary choices.  Making a copy also of the context
3540  * tableau would obviate the need to undo any changes made to it later,
3541  * while taking a snapshot of the main tableau could reduce memory usage.
3542  * If we were to switch to taking a snapshot of the main tableau,
3543  * we would have to keep in mind that we need to save the row signs
3544  * and that we need to do this before saving the current basis
3545  * such that the basis has been restore before we restore the row signs.
3546  */
3547 static void find_in_pos(struct isl_sol *sol, struct isl_tab *tab, isl_int *ineq)
3548 {
3549         void *saved;
3550
3551         if (!sol->context)
3552                 goto error;
3553         saved = sol->context->op->save(sol->context);
3554
3555         tab = isl_tab_dup(tab);
3556         if (!tab)
3557                 goto error;
3558
3559         sol->context->op->add_ineq(sol->context, ineq, 0, 1);
3560
3561         find_solutions(sol, tab);
3562
3563         if (!sol->error)
3564                 sol->context->op->restore(sol->context, saved);
3565         return;
3566 error:
3567         sol->error = 1;
3568 }
3569
3570 /* Record the absence of solutions for those values of the parameters
3571  * that do not satisfy the given inequality with equality.
3572  */
3573 static void no_sol_in_strict(struct isl_sol *sol,
3574         struct isl_tab *tab, struct isl_vec *ineq)
3575 {
3576         int empty;
3577         void *saved;
3578
3579         if (!sol->context || sol->error)
3580                 goto error;
3581         saved = sol->context->op->save(sol->context);
3582
3583         isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3584
3585         sol->context->op->add_ineq(sol->context, ineq->el, 1, 0);
3586         if (!sol->context)
3587                 goto error;
3588
3589         empty = tab->empty;
3590         tab->empty = 1;
3591         sol_add(sol, tab);
3592         tab->empty = empty;
3593
3594         isl_int_add_ui(ineq->el[0], ineq->el[0], 1);
3595
3596         sol->context->op->restore(sol->context, saved);
3597         return;
3598 error:
3599         sol->error = 1;
3600 }
3601
3602 /* Compute the lexicographic minimum of the set represented by the main
3603  * tableau "tab" within the context "sol->context_tab".
3604  * On entry the sample value of the main tableau is lexicographically
3605  * less than or equal to this lexicographic minimum.
3606  * Pivots are performed until a feasible point is found, which is then
3607  * necessarily equal to the minimum, or until the tableau is found to
3608  * be infeasible.  Some pivots may need to be performed for only some
3609  * feasible values of the context tableau.  If so, the context tableau
3610  * is split into a part where the pivot is needed and a part where it is not.
3611  *
3612  * Whenever we enter the main loop, the main tableau is such that no
3613  * "obvious" pivots need to be performed on it, where "obvious" means
3614  * that the given row can be seen to be negative without looking at
3615  * the context tableau.  In particular, for non-parametric problems,
3616  * no pivots need to be performed on the main tableau.
3617  * The caller of find_solutions is responsible for making this property
3618  * hold prior to the first iteration of the loop, while restore_lexmin
3619  * is called before every other iteration.
3620  *
3621  * Inside the main loop, we first examine the signs of the rows of
3622  * the main tableau within the context of the context tableau.
3623  * If we find a row that is always non-positive for all values of
3624  * the parameters satisfying the context tableau and negative for at
3625  * least one value of the parameters, we perform the appropriate pivot
3626  * and start over.  An exception is the case where no pivot can be
3627  * performed on the row.  In this case, we require that the sign of
3628  * the row is negative for all values of the parameters (rather than just
3629  * non-positive).  This special case is handled inside row_sign, which
3630  * will say that the row can have any sign if it determines that it can
3631  * attain both negative and zero values.
3632  *
3633  * If we can't find a row that always requires a pivot, but we can find
3634  * one or more rows that require a pivot for some values of the parameters
3635  * (i.e., the row can attain both positive and negative signs), then we split
3636  * the context tableau into two parts, one where we force the sign to be
3637  * non-negative and one where we force is to be negative.
3638  * The non-negative part is handled by a recursive call (through find_in_pos).
3639  * Upon returning from this call, we continue with the negative part and
3640  * perform the required pivot.
3641  *
3642  * If no such rows can be found, all rows are non-negative and we have
3643  * found a (rational) feasible point.  If we only wanted a rational point
3644  * then we are done.
3645  * Otherwise, we check if all values of the sample point of the tableau
3646  * are integral for the variables.  If so, we have found the minimal
3647  * integral point and we are done.
3648  * If the sample point is not integral, then we need to make a distinction
3649  * based on whether the constant term is non-integral or the coefficients
3650  * of the parameters.  Furthermore, in order to decide how to handle
3651  * the non-integrality, we also need to know whether the coefficients
3652  * of the other columns in the tableau are integral.  This leads
3653  * to the following table.  The first two rows do not correspond
3654  * to a non-integral sample point and are only mentioned for completeness.
3655  *
3656  *      constant        parameters      other
3657  *
3658  *      int             int             int     |
3659  *      int             int             rat     | -> no problem
3660  *
3661  *      rat             int             int       -> fail
3662  *
3663  *      rat             int             rat       -> cut
3664  *
3665  *      int             rat             rat     |
3666  *      rat             rat             rat     | -> parametric cut
3667  *
3668  *      int             rat             int     |
3669  *      rat             rat             int     | -> split context
3670  *
3671  * If the parametric constant is completely integral, then there is nothing
3672  * to be done.  If the constant term is non-integral, but all the other
3673  * coefficient are integral, then there is nothing that can be done
3674  * and the tableau has no integral solution.
3675  * If, on the other hand, one or more of the other columns have rational
3676  * coefficients, but the parameter coefficients are all integral, then
3677  * we can perform a regular (non-parametric) cut.
3678  * Finally, if there is any parameter coefficient that is non-integral,
3679  * then we need to involve the context tableau.  There are two cases here.
3680  * If at least one other column has a rational coefficient, then we
3681  * can perform a parametric cut in the main tableau by adding a new
3682  * integer division in the context tableau.
3683  * If all other columns have integral coefficients, then we need to
3684  * enforce that the rational combination of parameters (c + \sum a_i y_i)/m
3685  * is always integral.  We do this by introducing an integer division
3686  * q = floor((c + \sum a_i y_i)/m) and stipulating that its argument should
3687  * always be integral in the context tableau, i.e., m q = c + \sum a_i y_i.
3688  * Since q is expressed in the tableau as
3689  *      c + \sum a_i y_i - m q >= 0
3690  *      -c - \sum a_i y_i + m q + m - 1 >= 0
3691  * it is sufficient to add the inequality
3692  *      -c - \sum a_i y_i + m q >= 0
3693  * In the part of the context where this inequality does not hold, the
3694  * main tableau is marked as being empty.
3695  */
3696 static void find_solutions(struct isl_sol *sol, struct isl_tab *tab)
3697 {
3698         struct isl_context *context;
3699         int r;
3700
3701         if (!tab || sol->error)
3702                 goto error;
3703
3704         context = sol->context;
3705
3706         if (tab->empty)
3707                 goto done;
3708         if (context->op->is_empty(context))
3709                 goto done;
3710
3711         for (r = 0; r >= 0 && tab && !tab->empty; r = restore_lexmin(tab)) {
3712                 int flags;
3713                 int row;
3714                 enum isl_tab_row_sign sgn;
3715                 int split = -1;
3716                 int n_split = 0;
3717
3718                 for (row = tab->n_redundant; row < tab->n_row; ++row) {
3719                         if (!isl_tab_var_from_row(tab, row)->is_nonneg)
3720                                 continue;
3721                         sgn = row_sign(tab, sol, row);
3722                         if (!sgn)
3723                                 goto error;
3724                         tab->row_sign[row] = sgn;
3725                         if (sgn == isl_tab_row_any)
3726                                 n_split++;
3727                         if (sgn == isl_tab_row_any && split == -1)
3728                                 split = row;
3729                         if (sgn == isl_tab_row_neg)
3730                                 break;
3731                 }
3732                 if (row < tab->n_row)
3733                         continue;
3734                 if (split != -1) {
3735                         struct isl_vec *ineq;
3736                         if (n_split != 1)
3737                                 split = context->op->best_split(context, tab);
3738                         if (split < 0)
3739                                 goto error;
3740                         ineq = get_row_parameter_ineq(tab, split);
3741                         if (!ineq)
3742                                 goto error;
3743                         is_strict(ineq);
3744                         for (row = tab->n_redundant; row < tab->n_row; ++row) {
3745                                 if (!isl_tab_var_from_row(tab, row)->is_nonneg)
3746                                         continue;
3747                                 if (tab->row_sign[row] == isl_tab_row_any)
3748                                         tab->row_sign[row] = isl_tab_row_unknown;
3749                         }
3750                         tab->row_sign[split] = isl_tab_row_pos;
3751                         sol_inc_level(sol);
3752                         find_in_pos(sol, tab, ineq->el);
3753                         tab->row_sign[split] = isl_tab_row_neg;
3754                         row = split;
3755                         isl_seq_neg(ineq->el, ineq->el, ineq->size);
3756                         isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3757                         if (!sol->error)
3758                                 context->op->add_ineq(context, ineq->el, 0, 1);
3759                         isl_vec_free(ineq);
3760                         if (sol->error)
3761                                 goto error;
3762                         continue;
3763                 }
3764                 if (tab->rational)
3765                         break;
3766                 row = first_non_integer_row(tab, &flags);
3767                 if (row < 0)
3768                         break;
3769                 if (ISL_FL_ISSET(flags, I_PAR)) {
3770                         if (ISL_FL_ISSET(flags, I_VAR)) {
3771                                 if (isl_tab_mark_empty(tab) < 0)
3772                                         goto error;
3773                                 break;
3774                         }
3775                         row = add_cut(tab, row);
3776                 } else if (ISL_FL_ISSET(flags, I_VAR)) {
3777                         struct isl_vec *div;
3778                         struct isl_vec *ineq;
3779                         int d;
3780                         div = get_row_split_div(tab, row);
3781                         if (!div)
3782                                 goto error;
3783                         d = context->op->get_div(context, tab, div);
3784                         isl_vec_free(div);
3785                         if (d < 0)
3786                                 goto error;
3787                         ineq = ineq_for_div(context->op->peek_basic_set(context), d);
3788                         if (!ineq)
3789                                 goto error;
3790                         sol_inc_level(sol);
3791                         no_sol_in_strict(sol, tab, ineq);
3792                         isl_seq_neg(ineq->el, ineq->el, ineq->size);
3793                         context->op->add_ineq(context, ineq->el, 1, 1);
3794                         isl_vec_free(ineq);
3795                         if (sol->error || !context->op->is_ok(context))
3796                                 goto error;
3797                         tab = set_row_cst_to_div(tab, row, d);
3798                         if (context->op->is_empty(context))
3799                                 break;
3800                 } else
3801                         row = add_parametric_cut(tab, row, context);
3802                 if (row < 0)
3803                         goto error;
3804         }
3805         if (r < 0)
3806                 goto error;
3807 done:
3808         sol_add(sol, tab);
3809         isl_tab_free(tab);
3810         return;
3811 error:
3812         isl_tab_free(tab);
3813         sol->error = 1;
3814 }
3815
3816 /* Compute the lexicographic minimum of the set represented by the main
3817  * tableau "tab" within the context "sol->context_tab".
3818  *
3819  * As a preprocessing step, we first transfer all the purely parametric
3820  * equalities from the main tableau to the context tableau, i.e.,
3821  * parameters that have been pivoted to a row.
3822  * These equalities are ignored by the main algorithm, because the
3823  * corresponding rows may not be marked as being non-negative.
3824  * In parts of the context where the added equality does not hold,
3825  * the main tableau is marked as being empty.
3826  */
3827 static void find_solutions_main(struct isl_sol *sol, struct isl_tab *tab)
3828 {
3829         int row;
3830
3831         if (!tab)
3832                 goto error;
3833
3834         sol->level = 0;
3835
3836         for (row = tab->n_redundant; row < tab->n_row; ++row) {
3837                 int p;
3838                 struct isl_vec *eq;
3839
3840                 if (tab->row_var[row] < 0)
3841                         continue;
3842                 if (tab->row_var[row] >= tab->n_param &&
3843                     tab->row_var[row] < tab->n_var - tab->n_div)
3844                         continue;
3845                 if (tab->row_var[row] < tab->n_param)
3846                         p = tab->row_var[row];
3847                 else
3848                         p = tab->row_var[row]
3849                                 + tab->n_param - (tab->n_var - tab->n_div);
3850
3851                 eq = isl_vec_alloc(tab->mat->ctx, 1+tab->n_param+tab->n_div);
3852                 if (!eq)
3853                         goto error;
3854                 get_row_parameter_line(tab, row, eq->el);
3855                 isl_int_neg(eq->el[1 + p], tab->mat->row[row][0]);
3856                 eq = isl_vec_normalize(eq);
3857
3858                 sol_inc_level(sol);
3859                 no_sol_in_strict(sol, tab, eq);
3860
3861                 isl_seq_neg(eq->el, eq->el, eq->size);
3862                 sol_inc_level(sol);
3863                 no_sol_in_strict(sol, tab, eq);
3864                 isl_seq_neg(eq->el, eq->el, eq->size);
3865
3866                 sol->context->op->add_eq(sol->context, eq->el, 1, 1);
3867
3868                 isl_vec_free(eq);
3869
3870                 if (isl_tab_mark_redundant(tab, row) < 0)
3871                         goto error;
3872
3873                 if (sol->context->op->is_empty(sol->context))
3874                         break;
3875
3876                 row = tab->n_redundant - 1;
3877         }
3878
3879         find_solutions(sol, tab);
3880
3881         sol->level = 0;
3882         sol_pop(sol);
3883
3884         return;
3885 error:
3886         isl_tab_free(tab);
3887         sol->error = 1;
3888 }
3889
3890 /* Check if integer division "div" of "dom" also occurs in "bmap".
3891  * If so, return its position within the divs.
3892  * If not, return -1.
3893  */
3894 static int find_context_div(struct isl_basic_map *bmap,
3895         struct isl_basic_set *dom, unsigned div)
3896 {
3897         int i;
3898         unsigned b_dim = isl_space_dim(bmap->dim, isl_dim_all);
3899         unsigned d_dim = isl_space_dim(dom->dim, isl_dim_all);
3900
3901         if (isl_int_is_zero(dom->div[div][0]))
3902                 return -1;
3903         if (isl_seq_first_non_zero(dom->div[div] + 2 + d_dim, dom->n_div) != -1)
3904                 return -1;
3905
3906         for (i = 0; i < bmap->n_div; ++i) {
3907                 if (isl_int_is_zero(bmap->div[i][0]))
3908                         continue;
3909                 if (isl_seq_first_non_zero(bmap->div[i] + 2 + d_dim,
3910                                            (b_dim - d_dim) + bmap->n_div) != -1)
3911                         continue;
3912                 if (isl_seq_eq(bmap->div[i], dom->div[div], 2 + d_dim))
3913                         return i;
3914         }
3915         return -1;
3916 }
3917
3918 /* The correspondence between the variables in the main tableau,
3919  * the context tableau, and the input map and domain is as follows.
3920  * The first n_param and the last n_div variables of the main tableau
3921  * form the variables of the context tableau.
3922  * In the basic map, these n_param variables correspond to the
3923  * parameters and the input dimensions.  In the domain, they correspond
3924  * to the parameters and the set dimensions.
3925  * The n_div variables correspond to the integer divisions in the domain.
3926  * To ensure that everything lines up, we may need to copy some of the
3927  * integer divisions of the domain to the map.  These have to be placed
3928  * in the same order as those in the context and they have to be placed
3929  * after any other integer divisions that the map may have.
3930  * This function performs the required reordering.
3931  */
3932 static struct isl_basic_map *align_context_divs(struct isl_basic_map *bmap,
3933         struct isl_basic_set *dom)
3934 {
3935         int i;
3936         int common = 0;
3937         int other;
3938
3939         for (i = 0; i < dom->n_div; ++i)
3940                 if (find_context_div(bmap, dom, i) != -1)
3941                         common++;
3942         other = bmap->n_div - common;
3943         if (dom->n_div - common > 0) {
3944                 bmap = isl_basic_map_extend_space(bmap, isl_space_copy(bmap->dim),
3945                                 dom->n_div - common, 0, 0);
3946                 if (!bmap)
3947                         return NULL;
3948         }
3949         for (i = 0; i < dom->n_div; ++i) {
3950                 int pos = find_context_div(bmap, dom, i);
3951                 if (pos < 0) {
3952                         pos = isl_basic_map_alloc_div(bmap);
3953                         if (pos < 0)
3954                                 goto error;
3955                         isl_int_set_si(bmap->div[pos][0], 0);
3956                 }
3957                 if (pos != other + i)
3958                         isl_basic_map_swap_div(bmap, pos, other + i);
3959         }
3960         return bmap;
3961 error:
3962         isl_basic_map_free(bmap);
3963         return NULL;
3964 }
3965
3966 /* Base case of isl_tab_basic_map_partial_lexopt, after removing
3967  * some obvious symmetries.
3968  *
3969  * We make sure the divs in the domain are properly ordered,
3970  * because they will be added one by one in the given order
3971  * during the construction of the solution map.
3972  */
3973 static struct isl_sol *basic_map_partial_lexopt_base(
3974         __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
3975         __isl_give isl_set **empty, int max,
3976         struct isl_sol *(*init)(__isl_keep isl_basic_map *bmap,
3977                     __isl_take isl_basic_set *dom, int track_empty, int max))
3978 {
3979         struct isl_tab *tab;
3980         struct isl_sol *sol = NULL;
3981         struct isl_context *context;
3982
3983         if (dom->n_div) {
3984                 dom = isl_basic_set_order_divs(dom);
3985                 bmap = align_context_divs(bmap, dom);
3986         }
3987         sol = init(bmap, dom, !!empty, max);
3988         if (!sol)
3989                 goto error;
3990
3991         context = sol->context;
3992         if (isl_basic_set_plain_is_empty(context->op->peek_basic_set(context)))
3993                 /* nothing */;
3994         else if (isl_basic_map_plain_is_empty(bmap)) {
3995                 if (sol->add_empty)
3996                         sol->add_empty(sol,
3997                     isl_basic_set_copy(context->op->peek_basic_set(context)));
3998         } else {
3999                 tab = tab_for_lexmin(bmap,
4000                                     context->op->peek_basic_set(context), 1, max);
4001                 tab = context->op->detect_nonnegative_parameters(context, tab);
4002                 find_solutions_main(sol, tab);
4003         }
4004         if (sol->error)
4005                 goto error;
4006
4007         isl_basic_map_free(bmap);
4008         return sol;
4009 error:
4010         sol_free(sol);
4011         isl_basic_map_free(bmap);
4012         return NULL;
4013 }
4014
4015 /* Base case of isl_tab_basic_map_partial_lexopt, after removing
4016  * some obvious symmetries.
4017  *
4018  * We call basic_map_partial_lexopt_base and extract the results.
4019  */
4020 static __isl_give isl_map *basic_map_partial_lexopt_base_map(
4021         __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
4022         __isl_give isl_set **empty, int max)
4023 {
4024         isl_map *result = NULL;
4025         struct isl_sol *sol;
4026         struct isl_sol_map *sol_map;
4027
4028         sol = basic_map_partial_lexopt_base(bmap, dom, empty, max,
4029                                             &sol_map_init);
4030         if (!sol)
4031                 return NULL;
4032         sol_map = (struct isl_sol_map *) sol;
4033
4034         result = isl_map_copy(sol_map->map);
4035         if (empty)
4036                 *empty = isl_set_copy(sol_map->empty);
4037         sol_free(&sol_map->sol);
4038         return result;
4039 }
4040
4041 /* Structure used during detection of parallel constraints.
4042  * n_in: number of "input" variables: isl_dim_param + isl_dim_in
4043  * n_out: number of "output" variables: isl_dim_out + isl_dim_div
4044  * val: the coefficients of the output variables
4045  */
4046 struct isl_constraint_equal_info {
4047         isl_basic_map *bmap;
4048         unsigned n_in;
4049         unsigned n_out;
4050         isl_int *val;
4051 };
4052
4053 /* Check whether the coefficients of the output variables
4054  * of the constraint in "entry" are equal to info->val.
4055  */
4056 static int constraint_equal(const void *entry, const void *val)
4057 {
4058         isl_int **row = (isl_int **)entry;
4059         const struct isl_constraint_equal_info *info = val;
4060
4061         return isl_seq_eq((*row) + 1 + info->n_in, info->val, info->n_out);
4062 }
4063
4064 /* Check whether "bmap" has a pair of constraints that have
4065  * the same coefficients for the output variables.
4066  * Note that the coefficients of the existentially quantified
4067  * variables need to be zero since the existentially quantified
4068  * of the result are usually not the same as those of the input.
4069  * the isl_dim_out and isl_dim_div dimensions.
4070  * If so, return 1 and return the row indices of the two constraints
4071  * in *first and *second.
4072  */
4073 static int parallel_constraints(__isl_keep isl_basic_map *bmap,
4074         int *first, int *second)
4075 {
4076         int i;
4077         isl_ctx *ctx = isl_basic_map_get_ctx(bmap);
4078         struct isl_hash_table *table = NULL;
4079         struct isl_hash_table_entry *entry;
4080         struct isl_constraint_equal_info info;
4081         unsigned n_out;
4082         unsigned n_div;
4083
4084         ctx = isl_basic_map_get_ctx(bmap);
4085         table = isl_hash_table_alloc(ctx, bmap->n_ineq);
4086         if (!table)
4087                 goto error;
4088
4089         info.n_in = isl_basic_map_dim(bmap, isl_dim_param) +
4090                     isl_basic_map_dim(bmap, isl_dim_in);
4091         info.bmap = bmap;
4092         n_out = isl_basic_map_dim(bmap, isl_dim_out);
4093         n_div = isl_basic_map_dim(bmap, isl_dim_div);
4094         info.n_out = n_out + n_div;
4095         for (i = 0; i < bmap->n_ineq; ++i) {
4096                 uint32_t hash;
4097
4098                 info.val = bmap->ineq[i] + 1 + info.n_in;
4099                 if (isl_seq_first_non_zero(info.val, n_out) < 0)
4100                         continue;
4101                 if (isl_seq_first_non_zero(info.val + n_out, n_div) >= 0)
4102                         continue;
4103                 hash = isl_seq_get_hash(info.val, info.n_out);
4104                 entry = isl_hash_table_find(ctx, table, hash,
4105                                             constraint_equal, &info, 1);
4106                 if (!entry)
4107                         goto error;
4108                 if (entry->data)
4109                         break;
4110                 entry->data = &bmap->ineq[i];
4111         }
4112
4113         if (i < bmap->n_ineq) {
4114                 *first = ((isl_int **)entry->data) - bmap->ineq;
4115                 *second = i;
4116         }
4117
4118         isl_hash_table_free(ctx, table);
4119
4120         return i < bmap->n_ineq;
4121 error:
4122         isl_hash_table_free(ctx, table);
4123         return -1;
4124 }
4125
4126 /* Given a set of upper bounds in "var", add constraints to "bset"
4127  * that make the i-th bound smallest.
4128  *
4129  * In particular, if there are n bounds b_i, then add the constraints
4130  *
4131  *      b_i <= b_j      for j > i
4132  *      b_i <  b_j      for j < i
4133  */
4134 static __isl_give isl_basic_set *select_minimum(__isl_take isl_basic_set *bset,
4135         __isl_keep isl_mat *var, int i)
4136 {
4137         isl_ctx *ctx;
4138         int j, k;
4139
4140         ctx = isl_mat_get_ctx(var);
4141
4142         for (j = 0; j < var->n_row; ++j) {
4143                 if (j == i)
4144                         continue;
4145                 k = isl_basic_set_alloc_inequality(bset);
4146                 if (k < 0)
4147                         goto error;
4148                 isl_seq_combine(bset->ineq[k], ctx->one, var->row[j],
4149                                 ctx->negone, var->row[i], var->n_col);
4150                 isl_int_set_si(bset->ineq[k][var->n_col], 0);
4151                 if (j < i)
4152                         isl_int_sub_ui(bset->ineq[k][0], bset->ineq[k][0], 1);
4153         }
4154
4155         bset = isl_basic_set_finalize(bset);
4156
4157         return bset;
4158 error:
4159         isl_basic_set_free(bset);
4160         return NULL;
4161 }
4162
4163 /* Given a set of upper bounds on the last "input" variable m,
4164  * construct a set that assigns the minimal upper bound to m, i.e.,
4165  * construct a set that divides the space into cells where one
4166  * of the upper bounds is smaller than all the others and assign
4167  * this upper bound to m.
4168  *
4169  * In particular, if there are n bounds b_i, then the result
4170  * consists of n basic sets, each one of the form
4171  *
4172  *      m = b_i
4173  *      b_i <= b_j      for j > i
4174  *      b_i <  b_j      for j < i
4175  */
4176 static __isl_give isl_set *set_minimum(__isl_take isl_space *dim,
4177         __isl_take isl_mat *var)
4178 {
4179         int i, k;
4180         isl_basic_set *bset = NULL;
4181         isl_ctx *ctx;
4182         isl_set *set = NULL;
4183
4184         if (!dim || !var)
4185                 goto error;
4186
4187         ctx = isl_space_get_ctx(dim);
4188         set = isl_set_alloc_space(isl_space_copy(dim),
4189                                 var->n_row, ISL_SET_DISJOINT);
4190
4191         for (i = 0; i < var->n_row; ++i) {
4192                 bset = isl_basic_set_alloc_space(isl_space_copy(dim), 0,
4193                                                1, var->n_row - 1);
4194                 k = isl_basic_set_alloc_equality(bset);
4195                 if (k < 0)
4196                         goto error;
4197                 isl_seq_cpy(bset->eq[k], var->row[i], var->n_col);
4198                 isl_int_set_si(bset->eq[k][var->n_col], -1);
4199                 bset = select_minimum(bset, var, i);
4200                 set = isl_set_add_basic_set(set, bset);
4201         }
4202
4203         isl_space_free(dim);
4204         isl_mat_free(var);
4205         return set;
4206 error:
4207         isl_basic_set_free(bset);
4208         isl_set_free(set);
4209         isl_space_free(dim);
4210         isl_mat_free(var);
4211         return NULL;
4212 }
4213
4214 /* Given that the last input variable of "bmap" represents the minimum
4215  * of the bounds in "cst", check whether we need to split the domain
4216  * based on which bound attains the minimum.
4217  *
4218  * A split is needed when the minimum appears in an integer division
4219  * or in an equality.  Otherwise, it is only needed if it appears in
4220  * an upper bound that is different from the upper bounds on which it
4221  * is defined.
4222  */
4223 static int need_split_basic_map(__isl_keep isl_basic_map *bmap,
4224         __isl_keep isl_mat *cst)
4225 {
4226         int i, j;
4227         unsigned total;
4228         unsigned pos;
4229
4230         pos = cst->n_col - 1;
4231         total = isl_basic_map_dim(bmap, isl_dim_all);
4232
4233         for (i = 0; i < bmap->n_div; ++i)
4234                 if (!isl_int_is_zero(bmap->div[i][2 + pos]))
4235                         return 1;
4236
4237         for (i = 0; i < bmap->n_eq; ++i)
4238                 if (!isl_int_is_zero(bmap->eq[i][1 + pos]))
4239                         return 1;
4240
4241         for (i = 0; i < bmap->n_ineq; ++i) {
4242                 if (isl_int_is_nonneg(bmap->ineq[i][1 + pos]))
4243                         continue;
4244                 if (!isl_int_is_negone(bmap->ineq[i][1 + pos]))
4245                         return 1;
4246                 if (isl_seq_first_non_zero(bmap->ineq[i] + 1 + pos + 1,
4247                                            total - pos - 1) >= 0)
4248                         return 1;
4249
4250                 for (j = 0; j < cst->n_row; ++j)
4251                         if (isl_seq_eq(bmap->ineq[i], cst->row[j], cst->n_col))
4252                                 break;
4253                 if (j >= cst->n_row)
4254                         return 1;
4255         }
4256
4257         return 0;
4258 }
4259
4260 /* Given that the last set variable of "bset" represents the minimum
4261  * of the bounds in "cst", check whether we need to split the domain
4262  * based on which bound attains the minimum.
4263  *
4264  * We simply call need_split_basic_map here.  This is safe because
4265  * the position of the minimum is computed from "cst" and not
4266  * from "bmap".
4267  */
4268 static int need_split_basic_set(__isl_keep isl_basic_set *bset,
4269         __isl_keep isl_mat *cst)
4270 {
4271         return need_split_basic_map((isl_basic_map *)bset, cst);
4272 }
4273
4274 /* Given that the last set variable of "set" represents the minimum
4275  * of the bounds in "cst", check whether we need to split the domain
4276  * based on which bound attains the minimum.
4277  */
4278 static int need_split_set(__isl_keep isl_set *set, __isl_keep isl_mat *cst)
4279 {
4280         int i;
4281
4282         for (i = 0; i < set->n; ++i)
4283                 if (need_split_basic_set(set->p[i], cst))
4284                         return 1;
4285
4286         return 0;
4287 }
4288
4289 /* Given a set of which the last set variable is the minimum
4290  * of the bounds in "cst", split each basic set in the set
4291  * in pieces where one of the bounds is (strictly) smaller than the others.
4292  * This subdivision is given in "min_expr".
4293  * The variable is subsequently projected out.
4294  *
4295  * We only do the split when it is needed.
4296  * For example if the last input variable m = min(a,b) and the only
4297  * constraints in the given basic set are lower bounds on m,
4298  * i.e., l <= m = min(a,b), then we can simply project out m
4299  * to obtain l <= a and l <= b, without having to split on whether
4300  * m is equal to a or b.
4301  */
4302 static __isl_give isl_set *split(__isl_take isl_set *empty,
4303         __isl_take isl_set *min_expr, __isl_take isl_mat *cst)
4304 {
4305         int n_in;
4306         int i;
4307         isl_space *dim;
4308         isl_set *res;
4309
4310         if (!empty || !min_expr || !cst)
4311                 goto error;
4312
4313         n_in = isl_set_dim(empty, isl_dim_set);
4314         dim = isl_set_get_space(empty);
4315         dim = isl_space_drop_dims(dim, isl_dim_set, n_in - 1, 1);
4316         res = isl_set_empty(dim);
4317
4318         for (i = 0; i < empty->n; ++i) {
4319                 isl_set *set;
4320
4321                 set = isl_set_from_basic_set(isl_basic_set_copy(empty->p[i]));
4322                 if (need_split_basic_set(empty->p[i], cst))
4323                         set = isl_set_intersect(set, isl_set_copy(min_expr));
4324                 set = isl_set_remove_dims(set, isl_dim_set, n_in - 1, 1);
4325
4326                 res = isl_set_union_disjoint(res, set);
4327         }
4328
4329         isl_set_free(empty);
4330         isl_set_free(min_expr);
4331         isl_mat_free(cst);
4332         return res;
4333 error:
4334         isl_set_free(empty);
4335         isl_set_free(min_expr);
4336         isl_mat_free(cst);
4337         return NULL;
4338 }
4339
4340 /* Given a map of which the last input variable is the minimum
4341  * of the bounds in "cst", split each basic set in the set
4342  * in pieces where one of the bounds is (strictly) smaller than the others.
4343  * This subdivision is given in "min_expr".
4344  * The variable is subsequently projected out.
4345  *
4346  * The implementation is essentially the same as that of "split".
4347  */
4348 static __isl_give isl_map *split_domain(__isl_take isl_map *opt,
4349         __isl_take isl_set *min_expr, __isl_take isl_mat *cst)
4350 {
4351         int n_in;
4352         int i;
4353         isl_space *dim;
4354         isl_map *res;
4355
4356         if (!opt || !min_expr || !cst)
4357                 goto error;
4358
4359         n_in = isl_map_dim(opt, isl_dim_in);
4360         dim = isl_map_get_space(opt);
4361         dim = isl_space_drop_dims(dim, isl_dim_in, n_in - 1, 1);
4362         res = isl_map_empty(dim);
4363
4364         for (i = 0; i < opt->n; ++i) {
4365                 isl_map *map;
4366
4367                 map = isl_map_from_basic_map(isl_basic_map_copy(opt->p[i]));
4368                 if (need_split_basic_map(opt->p[i], cst))
4369                         map = isl_map_intersect_domain(map,
4370                                                        isl_set_copy(min_expr));
4371                 map = isl_map_remove_dims(map, isl_dim_in, n_in - 1, 1);
4372
4373                 res = isl_map_union_disjoint(res, map);
4374         }
4375
4376         isl_map_free(opt);
4377         isl_set_free(min_expr);
4378         isl_mat_free(cst);
4379         return res;
4380 error:
4381         isl_map_free(opt);
4382         isl_set_free(min_expr);
4383         isl_mat_free(cst);
4384         return NULL;
4385 }
4386
4387 static __isl_give isl_map *basic_map_partial_lexopt(
4388         __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
4389         __isl_give isl_set **empty, int max);
4390
4391 union isl_lex_res {
4392         void *p;
4393         isl_map *map;
4394         isl_pw_multi_aff *pma;
4395 };
4396
4397 /* This function is called from basic_map_partial_lexopt_symm.
4398  * The last variable of "bmap" and "dom" corresponds to the minimum
4399  * of the bounds in "cst".  "map_space" is the space of the original
4400  * input relation (of basic_map_partial_lexopt_symm) and "set_space"
4401  * is the space of the original domain.
4402  *
4403  * We recursively call basic_map_partial_lexopt and then plug in
4404  * the definition of the minimum in the result.
4405  */
4406 static __isl_give union isl_lex_res basic_map_partial_lexopt_symm_map_core(
4407         __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
4408         __isl_give isl_set **empty, int max, __isl_take isl_mat *cst,
4409         __isl_take isl_space *map_space, __isl_take isl_space *set_space)
4410 {
4411         isl_map *opt;
4412         isl_set *min_expr;
4413         union isl_lex_res res;
4414
4415         min_expr = set_minimum(isl_basic_set_get_space(dom), isl_mat_copy(cst));
4416
4417         opt = basic_map_partial_lexopt(bmap, dom, empty, max);
4418
4419         if (empty) {
4420                 *empty = split(*empty,
4421                                isl_set_copy(min_expr), isl_mat_copy(cst));
4422                 *empty = isl_set_reset_space(*empty, set_space);
4423         }
4424
4425         opt = split_domain(opt, min_expr, cst);
4426         opt = isl_map_reset_space(opt, map_space);
4427
4428         res.map = opt;
4429         return res;
4430 }
4431
4432 /* Given a basic map with at least two parallel constraints (as found
4433  * by the function parallel_constraints), first look for more constraints
4434  * parallel to the two constraint and replace the found list of parallel
4435  * constraints by a single constraint with as "input" part the minimum
4436  * of the input parts of the list of constraints.  Then, recursively call
4437  * basic_map_partial_lexopt (possibly finding more parallel constraints)
4438  * and plug in the definition of the minimum in the result.
4439  *
4440  * More specifically, given a set of constraints
4441  *
4442  *      a x + b_i(p) >= 0
4443  *
4444  * Replace this set by a single constraint
4445  *
4446  *      a x + u >= 0
4447  *
4448  * with u a new parameter with constraints
4449  *
4450  *      u <= b_i(p)
4451  *
4452  * Any solution to the new system is also a solution for the original system
4453  * since
4454  *
4455  *      a x >= -u >= -b_i(p)
4456  *
4457  * Moreover, m = min_i(b_i(p)) satisfies the constraints on u and can
4458  * therefore be plugged into the solution.
4459  */
4460 static union isl_lex_res basic_map_partial_lexopt_symm(
4461         __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
4462         __isl_give isl_set **empty, int max, int first, int second,
4463         __isl_give union isl_lex_res (*core)(__isl_take isl_basic_map *bmap,
4464                                             __isl_take isl_basic_set *dom,
4465                                             __isl_give isl_set **empty,
4466                                             int max, __isl_take isl_mat *cst,
4467                                             __isl_take isl_space *map_space,
4468                                             __isl_take isl_space *set_space))
4469 {
4470         int i, n, k;
4471         int *list = NULL;
4472         unsigned n_in, n_out, n_div;
4473         isl_ctx *ctx;
4474         isl_vec *var = NULL;
4475         isl_mat *cst = NULL;
4476         isl_space *map_space, *set_space;
4477         union isl_lex_res res;
4478
4479         map_space = isl_basic_map_get_space(bmap);
4480         set_space = empty ? isl_basic_set_get_space(dom) : NULL;
4481
4482         n_in = isl_basic_map_dim(bmap, isl_dim_param) +
4483                isl_basic_map_dim(bmap, isl_dim_in);
4484         n_out = isl_basic_map_dim(bmap, isl_dim_all) - n_in;
4485
4486         ctx = isl_basic_map_get_ctx(bmap);
4487         list = isl_alloc_array(ctx, int, bmap->n_ineq);
4488         var = isl_vec_alloc(ctx, n_out);
4489         if (!list || !var)
4490                 goto error;
4491
4492         list[0] = first;
4493         list[1] = second;
4494         isl_seq_cpy(var->el, bmap->ineq[first] + 1 + n_in, n_out);
4495         for (i = second + 1, n = 2; i < bmap->n_ineq; ++i) {
4496                 if (isl_seq_eq(var->el, bmap->ineq[i] + 1 + n_in, n_out))
4497                         list[n++] = i;
4498         }
4499
4500         cst = isl_mat_alloc(ctx, n, 1 + n_in);
4501         if (!cst)
4502                 goto error;
4503
4504         for (i = 0; i < n; ++i)
4505                 isl_seq_cpy(cst->row[i], bmap->ineq[list[i]], 1 + n_in);
4506
4507         bmap = isl_basic_map_cow(bmap);
4508         if (!bmap)
4509                 goto error;
4510         for (i = n - 1; i >= 0; --i)
4511                 if (isl_basic_map_drop_inequality(bmap, list[i]) < 0)
4512                         goto error;
4513
4514         bmap = isl_basic_map_add(bmap, isl_dim_in, 1);
4515         bmap = isl_basic_map_extend_constraints(bmap, 0, 1);
4516         k = isl_basic_map_alloc_inequality(bmap);
4517         if (k < 0)
4518                 goto error;
4519         isl_seq_clr(bmap->ineq[k], 1 + n_in);
4520         isl_int_set_si(bmap->ineq[k][1 + n_in], 1);
4521         isl_seq_cpy(bmap->ineq[k] + 1 + n_in + 1, var->el, n_out);
4522         bmap = isl_basic_map_finalize(bmap);
4523
4524         n_div = isl_basic_set_dim(dom, isl_dim_div);
4525         dom = isl_basic_set_add_dims(dom, isl_dim_set, 1);
4526         dom = isl_basic_set_extend_constraints(dom, 0, n);
4527         for (i = 0; i < n; ++i) {
4528                 k = isl_basic_set_alloc_inequality(dom);
4529                 if (k < 0)
4530                         goto error;
4531                 isl_seq_cpy(dom->ineq[k], cst->row[i], 1 + n_in);
4532                 isl_int_set_si(dom->ineq[k][1 + n_in], -1);
4533                 isl_seq_clr(dom->ineq[k] + 1 + n_in + 1, n_div);
4534         }
4535
4536         isl_vec_free(var);
4537         free(list);
4538
4539         return core(bmap, dom, empty, max, cst, map_space, set_space);
4540 error:
4541         isl_space_free(map_space);
4542         isl_space_free(set_space);
4543         isl_mat_free(cst);
4544         isl_vec_free(var);
4545         free(list);
4546         isl_basic_set_free(dom);
4547         isl_basic_map_free(bmap);
4548         res.p = NULL;
4549         return res;
4550 }
4551
4552 static __isl_give isl_map *basic_map_partial_lexopt_symm_map(
4553         __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
4554         __isl_give isl_set **empty, int max, int first, int second)
4555 {
4556         return basic_map_partial_lexopt_symm(bmap, dom, empty, max,
4557                     first, second, &basic_map_partial_lexopt_symm_map_core).map;
4558 }
4559
4560 /* Recursive part of isl_tab_basic_map_partial_lexopt, after detecting
4561  * equalities and removing redundant constraints.
4562  *
4563  * We first check if there are any parallel constraints (left).
4564  * If not, we are in the base case.
4565  * If there are parallel constraints, we replace them by a single
4566  * constraint in basic_map_partial_lexopt_symm and then call
4567  * this function recursively to look for more parallel constraints.
4568  */
4569 static __isl_give isl_map *basic_map_partial_lexopt(
4570         __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
4571         __isl_give isl_set **empty, int max)
4572 {
4573         int par = 0;
4574         int first, second;
4575
4576         if (!bmap)
4577                 goto error;
4578
4579         if (bmap->ctx->opt->pip_symmetry)
4580                 par = parallel_constraints(bmap, &first, &second);
4581         if (par < 0)
4582                 goto error;
4583         if (!par)
4584                 return basic_map_partial_lexopt_base_map(bmap, dom, empty, max);
4585
4586         return basic_map_partial_lexopt_symm_map(bmap, dom, empty, max,
4587                                                  first, second);
4588 error:
4589         isl_basic_set_free(dom);
4590         isl_basic_map_free(bmap);
4591         return NULL;
4592 }
4593
4594 /* Compute the lexicographic minimum (or maximum if "max" is set)
4595  * of "bmap" over the domain "dom" and return the result as a map.
4596  * If "empty" is not NULL, then *empty is assigned a set that
4597  * contains those parts of the domain where there is no solution.
4598  * If "bmap" is marked as rational (ISL_BASIC_MAP_RATIONAL),
4599  * then we compute the rational optimum.  Otherwise, we compute
4600  * the integral optimum.
4601  *
4602  * We perform some preprocessing.  As the PILP solver does not
4603  * handle implicit equalities very well, we first make sure all
4604  * the equalities are explicitly available.
4605  *
4606  * We also add context constraints to the basic map and remove
4607  * redundant constraints.  This is only needed because of the
4608  * way we handle simple symmetries.  In particular, we currently look
4609  * for symmetries on the constraints, before we set up the main tableau.
4610  * It is then no good to look for symmetries on possibly redundant constraints.
4611  */
4612 struct isl_map *isl_tab_basic_map_partial_lexopt(
4613                 struct isl_basic_map *bmap, struct isl_basic_set *dom,
4614                 struct isl_set **empty, int max)
4615 {
4616         if (empty)
4617                 *empty = NULL;
4618         if (!bmap || !dom)
4619                 goto error;
4620
4621         isl_assert(bmap->ctx,
4622             isl_basic_map_compatible_domain(bmap, dom), goto error);
4623
4624         if (isl_basic_set_dim(dom, isl_dim_all) == 0)
4625                 return basic_map_partial_lexopt(bmap, dom, empty, max);
4626
4627         bmap = isl_basic_map_intersect_domain(bmap, isl_basic_set_copy(dom));
4628         bmap = isl_basic_map_detect_equalities(bmap);
4629         bmap = isl_basic_map_remove_redundancies(bmap);
4630
4631         return basic_map_partial_lexopt(bmap, dom, empty, max);
4632 error:
4633         isl_basic_set_free(dom);
4634         isl_basic_map_free(bmap);
4635         return NULL;
4636 }
4637
4638 struct isl_sol_for {
4639         struct isl_sol  sol;
4640         int             (*fn)(__isl_take isl_basic_set *dom,
4641                                 __isl_take isl_aff_list *list, void *user);
4642         void            *user;
4643 };
4644
4645 static void sol_for_free(struct isl_sol_for *sol_for)
4646 {
4647         if (sol_for->sol.context)
4648                 sol_for->sol.context->op->free(sol_for->sol.context);
4649         free(sol_for);
4650 }
4651
4652 static void sol_for_free_wrap(struct isl_sol *sol)
4653 {
4654         sol_for_free((struct isl_sol_for *)sol);
4655 }
4656
4657 /* Add the solution identified by the tableau and the context tableau.
4658  *
4659  * See documentation of sol_add for more details.
4660  *
4661  * Instead of constructing a basic map, this function calls a user
4662  * defined function with the current context as a basic set and
4663  * a list of affine expressions representing the relation between
4664  * the input and output.  The space over which the affine expressions
4665  * are defined is the same as that of the domain.  The number of
4666  * affine expressions in the list is equal to the number of output variables.
4667  */
4668 static void sol_for_add(struct isl_sol_for *sol,
4669         struct isl_basic_set *dom, struct isl_mat *M)
4670 {
4671         int i;
4672         isl_ctx *ctx;
4673         isl_local_space *ls;
4674         isl_aff *aff;
4675         isl_aff_list *list;
4676
4677         if (sol->sol.error || !dom || !M)
4678                 goto error;
4679
4680         ctx = isl_basic_set_get_ctx(dom);
4681         ls = isl_basic_set_get_local_space(dom);
4682         list = isl_aff_list_alloc(ctx, M->n_row - 1);
4683         for (i = 1; i < M->n_row; ++i) {
4684                 aff = isl_aff_alloc(isl_local_space_copy(ls));
4685                 if (aff) {
4686                         isl_int_set(aff->v->el[0], M->row[0][0]);
4687                         isl_seq_cpy(aff->v->el + 1, M->row[i], M->n_col);
4688                 }
4689                 aff = isl_aff_normalize(aff);
4690                 list = isl_aff_list_add(list, aff);
4691         }
4692         isl_local_space_free(ls);
4693
4694         dom = isl_basic_set_finalize(dom);
4695
4696         if (sol->fn(isl_basic_set_copy(dom), list, sol->user) < 0)
4697                 goto error;
4698
4699         isl_basic_set_free(dom);
4700         isl_mat_free(M);
4701         return;
4702 error:
4703         isl_basic_set_free(dom);
4704         isl_mat_free(M);
4705         sol->sol.error = 1;
4706 }
4707
4708 static void sol_for_add_wrap(struct isl_sol *sol,
4709         struct isl_basic_set *dom, struct isl_mat *M)
4710 {
4711         sol_for_add((struct isl_sol_for *)sol, dom, M);
4712 }
4713
4714 static struct isl_sol_for *sol_for_init(struct isl_basic_map *bmap, int max,
4715         int (*fn)(__isl_take isl_basic_set *dom, __isl_take isl_aff_list *list,
4716                   void *user),
4717         void *user)
4718 {
4719         struct isl_sol_for *sol_for = NULL;
4720         isl_space *dom_dim;
4721         struct isl_basic_set *dom = NULL;
4722
4723         sol_for = isl_calloc_type(bmap->ctx, struct isl_sol_for);
4724         if (!sol_for)
4725                 goto error;
4726
4727         dom_dim = isl_space_domain(isl_space_copy(bmap->dim));
4728         dom = isl_basic_set_universe(dom_dim);
4729
4730         sol_for->sol.rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
4731         sol_for->sol.dec_level.callback.run = &sol_dec_level_wrap;
4732         sol_for->sol.dec_level.sol = &sol_for->sol;
4733         sol_for->fn = fn;
4734         sol_for->user = user;
4735         sol_for->sol.max = max;
4736         sol_for->sol.n_out = isl_basic_map_dim(bmap, isl_dim_out);
4737         sol_for->sol.add = &sol_for_add_wrap;
4738         sol_for->sol.add_empty = NULL;
4739         sol_for->sol.free = &sol_for_free_wrap;
4740
4741         sol_for->sol.context = isl_context_alloc(dom);
4742         if (!sol_for->sol.context)
4743                 goto error;
4744
4745         isl_basic_set_free(dom);
4746         return sol_for;
4747 error:
4748         isl_basic_set_free(dom);
4749         sol_for_free(sol_for);
4750         return NULL;
4751 }
4752
4753 static void sol_for_find_solutions(struct isl_sol_for *sol_for,
4754         struct isl_tab *tab)
4755 {
4756         find_solutions_main(&sol_for->sol, tab);
4757 }
4758
4759 int isl_basic_map_foreach_lexopt(__isl_keep isl_basic_map *bmap, int max,
4760         int (*fn)(__isl_take isl_basic_set *dom, __isl_take isl_aff_list *list,
4761                   void *user),
4762         void *user)
4763 {
4764         struct isl_sol_for *sol_for = NULL;
4765
4766         bmap = isl_basic_map_copy(bmap);
4767         if (!bmap)
4768                 return -1;
4769
4770         bmap = isl_basic_map_detect_equalities(bmap);
4771         sol_for = sol_for_init(bmap, max, fn, user);
4772
4773         if (isl_basic_map_plain_is_empty(bmap))
4774                 /* nothing */;
4775         else {
4776                 struct isl_tab *tab;
4777                 struct isl_context *context = sol_for->sol.context;
4778                 tab = tab_for_lexmin(bmap,
4779                                 context->op->peek_basic_set(context), 1, max);
4780                 tab = context->op->detect_nonnegative_parameters(context, tab);
4781                 sol_for_find_solutions(sol_for, tab);
4782                 if (sol_for->sol.error)
4783                         goto error;
4784         }
4785
4786         sol_free(&sol_for->sol);
4787         isl_basic_map_free(bmap);
4788         return 0;
4789 error:
4790         sol_free(&sol_for->sol);
4791         isl_basic_map_free(bmap);
4792         return -1;
4793 }
4794
4795 int isl_basic_set_foreach_lexopt(__isl_keep isl_basic_set *bset, int max,
4796         int (*fn)(__isl_take isl_basic_set *dom, __isl_take isl_aff_list *list,
4797                   void *user),
4798         void *user)
4799 {
4800         return isl_basic_map_foreach_lexopt(bset, max, fn, user);
4801 }
4802
4803 /* Check if the given sequence of len variables starting at pos
4804  * represents a trivial (i.e., zero) solution.
4805  * The variables are assumed to be non-negative and to come in pairs,
4806  * with each pair representing a variable of unrestricted sign.
4807  * The solution is trivial if each such pair in the sequence consists
4808  * of two identical values, meaning that the variable being represented
4809  * has value zero.
4810  */
4811 static int region_is_trivial(struct isl_tab *tab, int pos, int len)
4812 {
4813         int i;
4814
4815         if (len == 0)
4816                 return 0;
4817
4818         for (i = 0; i < len; i +=  2) {
4819                 int neg_row;
4820                 int pos_row;
4821
4822                 neg_row = tab->var[pos + i].is_row ?
4823                                 tab->var[pos + i].index : -1;
4824                 pos_row = tab->var[pos + i + 1].is_row ?
4825                                 tab->var[pos + i + 1].index : -1;
4826
4827                 if ((neg_row < 0 ||
4828                      isl_int_is_zero(tab->mat->row[neg_row][1])) &&
4829                     (pos_row < 0 ||
4830                      isl_int_is_zero(tab->mat->row[pos_row][1])))
4831                         continue;
4832
4833                 if (neg_row < 0 || pos_row < 0)
4834                         return 0;
4835                 if (isl_int_ne(tab->mat->row[neg_row][1],
4836                                tab->mat->row[pos_row][1]))
4837                         return 0;
4838         }
4839
4840         return 1;
4841 }
4842
4843 /* Return the index of the first trivial region or -1 if all regions
4844  * are non-trivial.
4845  */
4846 static int first_trivial_region(struct isl_tab *tab,
4847         int n_region, struct isl_region *region)
4848 {
4849         int i;
4850
4851         for (i = 0; i < n_region; ++i) {
4852                 if (region_is_trivial(tab, region[i].pos, region[i].len))
4853                         return i;
4854         }
4855
4856         return -1;
4857 }
4858
4859 /* Check if the solution is optimal, i.e., whether the first
4860  * n_op entries are zero.
4861  */
4862 static int is_optimal(__isl_keep isl_vec *sol, int n_op)
4863 {
4864         int i;
4865
4866         for (i = 0; i < n_op; ++i)
4867                 if (!isl_int_is_zero(sol->el[1 + i]))
4868                         return 0;
4869         return 1;
4870 }
4871
4872 /* Add constraints to "tab" that ensure that any solution is significantly
4873  * better that that represented by "sol".  That is, find the first
4874  * relevant (within first n_op) non-zero coefficient and force it (along
4875  * with all previous coefficients) to be zero.
4876  * If the solution is already optimal (all relevant coefficients are zero),
4877  * then just mark the table as empty.
4878  */
4879 static int force_better_solution(struct isl_tab *tab,
4880         __isl_keep isl_vec *sol, int n_op)
4881 {
4882         int i;
4883         isl_ctx *ctx;
4884         isl_vec *v = NULL;
4885
4886         if (!sol)
4887                 return -1;
4888
4889         for (i = 0; i < n_op; ++i)
4890                 if (!isl_int_is_zero(sol->el[1 + i]))
4891                         break;
4892
4893         if (i == n_op) {
4894                 if (isl_tab_mark_empty(tab) < 0)
4895                         return -1;
4896                 return 0;
4897         }
4898
4899         ctx = isl_vec_get_ctx(sol);
4900         v = isl_vec_alloc(ctx, 1 + tab->n_var);
4901         if (!v)
4902                 return -1;
4903
4904         for (; i >= 0; --i) {
4905                 v = isl_vec_clr(v);
4906                 isl_int_set_si(v->el[1 + i], -1);
4907                 if (add_lexmin_eq(tab, v->el) < 0)
4908                         goto error;
4909         }
4910
4911         isl_vec_free(v);
4912         return 0;
4913 error:
4914         isl_vec_free(v);
4915         return -1;
4916 }
4917
4918 struct isl_trivial {
4919         int update;
4920         int region;
4921         int side;
4922         struct isl_tab_undo *snap;
4923 };
4924
4925 /* Return the lexicographically smallest non-trivial solution of the
4926  * given ILP problem.
4927  *
4928  * All variables are assumed to be non-negative.
4929  *
4930  * n_op is the number of initial coordinates to optimize.
4931  * That is, once a solution has been found, we will only continue looking
4932  * for solution that result in significantly better values for those
4933  * initial coordinates.  That is, we only continue looking for solutions
4934  * that increase the number of initial zeros in this sequence.
4935  *
4936  * A solution is non-trivial, if it is non-trivial on each of the
4937  * specified regions.  Each region represents a sequence of pairs
4938  * of variables.  A solution is non-trivial on such a region if
4939  * at least one of these pairs consists of different values, i.e.,
4940  * such that the non-negative variable represented by the pair is non-zero.
4941  *
4942  * Whenever a conflict is encountered, all constraints involved are
4943  * reported to the caller through a call to "conflict".
4944  *
4945  * We perform a simple branch-and-bound backtracking search.
4946  * Each level in the search represents initially trivial region that is forced
4947  * to be non-trivial.
4948  * At each level we consider n cases, where n is the length of the region.
4949  * In terms of the n/2 variables of unrestricted signs being encoded by
4950  * the region, we consider the cases
4951  *      x_0 >= 1
4952  *      x_0 <= -1
4953  *      x_0 = 0 and x_1 >= 1
4954  *      x_0 = 0 and x_1 <= -1
4955  *      x_0 = 0 and x_1 = 0 and x_2 >= 1
4956  *      x_0 = 0 and x_1 = 0 and x_2 <= -1
4957  *      ...
4958  * The cases are considered in this order, assuming that each pair
4959  * x_i_a x_i_b represents the value x_i_b - x_i_a.
4960  * That is, x_0 >= 1 is enforced by adding the constraint
4961  *      x_0_b - x_0_a >= 1
4962  */
4963 __isl_give isl_vec *isl_tab_basic_set_non_trivial_lexmin(
4964         __isl_take isl_basic_set *bset, int n_op, int n_region,
4965         struct isl_region *region,
4966         int (*conflict)(int con, void *user), void *user)
4967 {
4968         int i, j;
4969         int r;
4970         isl_ctx *ctx = isl_basic_set_get_ctx(bset);
4971         isl_vec *v = NULL;
4972         isl_vec *sol = isl_vec_alloc(ctx, 0);
4973         struct isl_tab *tab;
4974         struct isl_trivial *triv = NULL;
4975         int level, init;
4976
4977         tab = tab_for_lexmin(bset, NULL, 0, 0);
4978         if (!tab)
4979                 goto error;
4980         tab->conflict = conflict;
4981         tab->conflict_user = user;
4982
4983         v = isl_vec_alloc(ctx, 1 + tab->n_var);
4984         triv = isl_calloc_array(ctx, struct isl_trivial, n_region);
4985         if (!v || !triv)
4986                 goto error;
4987
4988         level = 0;
4989         init = 1;
4990
4991         while (level >= 0) {
4992                 int side, base;
4993
4994                 if (init) {
4995                         tab = cut_to_integer_lexmin(tab, CUT_ONE);
4996                         if (!tab)
4997                                 goto error;
4998                         if (tab->empty)
4999                                 goto backtrack;
5000                         r = first_trivial_region(tab, n_region, region);
5001                         if (r < 0) {
5002                                 for (i = 0; i < level; ++i)
5003                                         triv[i].update = 1;
5004                                 isl_vec_free(sol);
5005                                 sol = isl_tab_get_sample_value(tab);
5006                                 if (!sol)
5007                                         goto error;
5008                                 if (is_optimal(sol, n_op))
5009                                         break;
5010                                 goto backtrack;
5011                         }
5012                         if (level >= n_region)
5013                                 isl_die(ctx, isl_error_internal,
5014                                         "nesting level too deep", goto error);
5015                         if (isl_tab_extend_cons(tab,
5016                                             2 * region[r].len + 2 * n_op) < 0)
5017                                 goto error;
5018                         triv[level].region = r;
5019                         triv[level].side = 0;
5020                 }
5021
5022                 r = triv[level].region;
5023                 side = triv[level].side;
5024                 base = 2 * (side/2);
5025
5026                 if (side >= region[r].len) {
5027 backtrack:
5028                         level--;
5029                         init = 0;
5030                         if (level >= 0)
5031                                 if (isl_tab_rollback(tab, triv[level].snap) < 0)
5032                                         goto error;
5033                         continue;
5034                 }
5035
5036                 if (triv[level].update) {
5037                         if (force_better_solution(tab, sol, n_op) < 0)
5038                                 goto error;
5039                         triv[level].update = 0;
5040                 }
5041
5042                 if (side == base && base >= 2) {
5043                         for (j = base - 2; j < base; ++j) {
5044                                 v = isl_vec_clr(v);
5045                                 isl_int_set_si(v->el[1 + region[r].pos + j], 1);
5046                                 if (add_lexmin_eq(tab, v->el) < 0)
5047                                         goto error;
5048                         }
5049                 }
5050
5051                 triv[level].snap = isl_tab_snap(tab);
5052                 if (isl_tab_push_basis(tab) < 0)
5053                         goto error;
5054
5055                 v = isl_vec_clr(v);
5056                 isl_int_set_si(v->el[0], -1);
5057                 isl_int_set_si(v->el[1 + region[r].pos + side], -1);
5058                 isl_int_set_si(v->el[1 + region[r].pos + (side ^ 1)], 1);
5059                 tab = add_lexmin_ineq(tab, v->el);
5060
5061                 triv[level].side++;
5062                 level++;
5063                 init = 1;
5064         }
5065
5066         free(triv);
5067         isl_vec_free(v);
5068         isl_tab_free(tab);
5069         isl_basic_set_free(bset);
5070
5071         return sol;
5072 error:
5073         free(triv);
5074         isl_vec_free(v);
5075         isl_tab_free(tab);
5076         isl_basic_set_free(bset);
5077         isl_vec_free(sol);
5078         return NULL;
5079 }
5080
5081 /* Return the lexicographically smallest rational point in "bset",
5082  * assuming that all variables are non-negative.
5083  * If "bset" is empty, then return a zero-length vector.
5084  */
5085 __isl_give isl_vec *isl_tab_basic_set_non_neg_lexmin(
5086         __isl_take isl_basic_set *bset)
5087 {
5088         struct isl_tab *tab;
5089         isl_ctx *ctx = isl_basic_set_get_ctx(bset);
5090         isl_vec *sol;
5091
5092         tab = tab_for_lexmin(bset, NULL, 0, 0);
5093         if (!tab)
5094                 goto error;
5095         if (tab->empty)
5096                 sol = isl_vec_alloc(ctx, 0);
5097         else
5098                 sol = isl_tab_get_sample_value(tab);
5099         isl_tab_free(tab);
5100         isl_basic_set_free(bset);
5101         return sol;
5102 error:
5103         isl_tab_free(tab);
5104         isl_basic_set_free(bset);
5105         return NULL;
5106 }
5107
5108 struct isl_sol_pma {
5109         struct isl_sol  sol;
5110         isl_pw_multi_aff *pma;
5111         isl_set *empty;
5112 };
5113
5114 static void sol_pma_free(struct isl_sol_pma *sol_pma)
5115 {
5116         if (!sol_pma)
5117                 return;
5118         if (sol_pma->sol.context)
5119                 sol_pma->sol.context->op->free(sol_pma->sol.context);
5120         isl_pw_multi_aff_free(sol_pma->pma);
5121         isl_set_free(sol_pma->empty);
5122         free(sol_pma);
5123 }
5124
5125 /* This function is called for parts of the context where there is
5126  * no solution, with "bset" corresponding to the context tableau.
5127  * Simply add the basic set to the set "empty".
5128  */
5129 static void sol_pma_add_empty(struct isl_sol_pma *sol,
5130         __isl_take isl_basic_set *bset)
5131 {
5132         if (!bset)
5133                 goto error;
5134         isl_assert(bset->ctx, sol->empty, goto error);
5135
5136         sol->empty = isl_set_grow(sol->empty, 1);
5137         bset = isl_basic_set_simplify(bset);
5138         bset = isl_basic_set_finalize(bset);
5139         sol->empty = isl_set_add_basic_set(sol->empty, bset);
5140         if (!sol->empty)
5141                 sol->sol.error = 1;
5142         return;
5143 error:
5144         isl_basic_set_free(bset);
5145         sol->sol.error = 1;
5146 }
5147
5148 /* Given a basic map "dom" that represents the context and an affine
5149  * matrix "M" that maps the dimensions of the context to the
5150  * output variables, construct an isl_pw_multi_aff with a single
5151  * cell corresponding to "dom" and affine expressions copied from "M".
5152  */
5153 static void sol_pma_add(struct isl_sol_pma *sol,
5154         __isl_take isl_basic_set *dom, __isl_take isl_mat *M)
5155 {
5156         int i;
5157         isl_local_space *ls;
5158         isl_aff *aff;
5159         isl_multi_aff *maff;
5160         isl_pw_multi_aff *pma;
5161
5162         maff = isl_multi_aff_alloc(isl_pw_multi_aff_get_space(sol->pma));
5163         ls = isl_basic_set_get_local_space(dom);
5164         for (i = 1; i < M->n_row; ++i) {
5165                 aff = isl_aff_alloc(isl_local_space_copy(ls));
5166                 if (aff) {
5167                         isl_int_set(aff->v->el[0], M->row[0][0]);
5168                         isl_seq_cpy(aff->v->el + 1, M->row[i], M->n_col);
5169                 }
5170                 aff = isl_aff_normalize(aff);
5171                 maff = isl_multi_aff_set_aff(maff, i - 1, aff);
5172         }
5173         isl_local_space_free(ls);
5174         isl_mat_free(M);
5175         dom = isl_basic_set_simplify(dom);
5176         dom = isl_basic_set_finalize(dom);
5177         pma = isl_pw_multi_aff_alloc(isl_set_from_basic_set(dom), maff);
5178         sol->pma = isl_pw_multi_aff_add_disjoint(sol->pma, pma);
5179         if (!sol->pma)
5180                 sol->sol.error = 1;
5181 }
5182
5183 static void sol_pma_free_wrap(struct isl_sol *sol)
5184 {
5185         sol_pma_free((struct isl_sol_pma *)sol);
5186 }
5187
5188 static void sol_pma_add_empty_wrap(struct isl_sol *sol,
5189         __isl_take isl_basic_set *bset)
5190 {
5191         sol_pma_add_empty((struct isl_sol_pma *)sol, bset);
5192 }
5193
5194 static void sol_pma_add_wrap(struct isl_sol *sol,
5195         __isl_take isl_basic_set *dom, __isl_take isl_mat *M)
5196 {
5197         sol_pma_add((struct isl_sol_pma *)sol, dom, M);
5198 }
5199
5200 /* Construct an isl_sol_pma structure for accumulating the solution.
5201  * If track_empty is set, then we also keep track of the parts
5202  * of the context where there is no solution.
5203  * If max is set, then we are solving a maximization, rather than
5204  * a minimization problem, which means that the variables in the
5205  * tableau have value "M - x" rather than "M + x".
5206  */
5207 static struct isl_sol *sol_pma_init(__isl_keep isl_basic_map *bmap,
5208         __isl_take isl_basic_set *dom, int track_empty, int max)
5209 {
5210         struct isl_sol_pma *sol_pma = NULL;
5211
5212         if (!bmap)
5213                 goto error;
5214
5215         sol_pma = isl_calloc_type(bmap->ctx, struct isl_sol_pma);
5216         if (!sol_pma)
5217                 goto error;
5218
5219         sol_pma->sol.rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
5220         sol_pma->sol.dec_level.callback.run = &sol_dec_level_wrap;
5221         sol_pma->sol.dec_level.sol = &sol_pma->sol;
5222         sol_pma->sol.max = max;
5223         sol_pma->sol.n_out = isl_basic_map_dim(bmap, isl_dim_out);
5224         sol_pma->sol.add = &sol_pma_add_wrap;
5225         sol_pma->sol.add_empty = track_empty ? &sol_pma_add_empty_wrap : NULL;
5226         sol_pma->sol.free = &sol_pma_free_wrap;
5227         sol_pma->pma = isl_pw_multi_aff_empty(isl_basic_map_get_space(bmap));
5228         if (!sol_pma->pma)
5229                 goto error;
5230
5231         sol_pma->sol.context = isl_context_alloc(dom);
5232         if (!sol_pma->sol.context)
5233                 goto error;
5234
5235         if (track_empty) {
5236                 sol_pma->empty = isl_set_alloc_space(isl_basic_set_get_space(dom),
5237                                                         1, ISL_SET_DISJOINT);
5238                 if (!sol_pma->empty)
5239                         goto error;
5240         }
5241
5242         isl_basic_set_free(dom);
5243         return &sol_pma->sol;
5244 error:
5245         isl_basic_set_free(dom);
5246         sol_pma_free(sol_pma);
5247         return NULL;
5248 }
5249
5250 /* Base case of isl_tab_basic_map_partial_lexopt, after removing
5251  * some obvious symmetries.
5252  *
5253  * We call basic_map_partial_lexopt_base and extract the results.
5254  */
5255 static __isl_give isl_pw_multi_aff *basic_map_partial_lexopt_base_pma(
5256         __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
5257         __isl_give isl_set **empty, int max)
5258 {
5259         isl_pw_multi_aff *result = NULL;
5260         struct isl_sol *sol;
5261         struct isl_sol_pma *sol_pma;
5262
5263         sol = basic_map_partial_lexopt_base(bmap, dom, empty, max,
5264                                             &sol_pma_init);
5265         if (!sol)
5266                 return NULL;
5267         sol_pma = (struct isl_sol_pma *) sol;
5268
5269         result = isl_pw_multi_aff_copy(sol_pma->pma);
5270         if (empty)
5271                 *empty = isl_set_copy(sol_pma->empty);
5272         sol_free(&sol_pma->sol);
5273         return result;
5274 }
5275
5276 /* Given that the last input variable of "maff" represents the minimum
5277  * of some bounds, check whether we need to plug in the expression
5278  * of the minimum.
5279  *
5280  * In particular, check if the last input variable appears in any
5281  * of the expressions in "maff".
5282  */
5283 static int need_substitution(__isl_keep isl_multi_aff *maff)
5284 {
5285         int i;
5286         unsigned pos;
5287
5288         pos = isl_multi_aff_dim(maff, isl_dim_in) - 1;
5289
5290         for (i = 0; i < maff->n; ++i)
5291                 if (isl_aff_involves_dims(maff->p[i], isl_dim_in, pos, 1))
5292                         return 1;
5293
5294         return 0;
5295 }
5296
5297 /* Given a set of upper bounds on the last "input" variable m,
5298  * construct a piecewise affine expression that selects
5299  * the minimal upper bound to m, i.e.,
5300  * divide the space into cells where one
5301  * of the upper bounds is smaller than all the others and select
5302  * this upper bound on that cell.
5303  *
5304  * In particular, if there are n bounds b_i, then the result
5305  * consists of n cell, each one of the form
5306  *
5307  *      b_i <= b_j      for j > i
5308  *      b_i <  b_j      for j < i
5309  *
5310  * The affine expression on this cell is
5311  *
5312  *      b_i
5313  */
5314 static __isl_give isl_pw_aff *set_minimum_pa(__isl_take isl_space *space,
5315         __isl_take isl_mat *var)
5316 {
5317         int i;
5318         isl_aff *aff = NULL;
5319         isl_basic_set *bset = NULL;
5320         isl_ctx *ctx;
5321         isl_pw_aff *paff = NULL;
5322         isl_space *pw_space;
5323         isl_local_space *ls = NULL;
5324
5325         if (!space || !var)
5326                 goto error;
5327
5328         ctx = isl_space_get_ctx(space);
5329         ls = isl_local_space_from_space(isl_space_copy(space));
5330         pw_space = isl_space_copy(space);
5331         pw_space = isl_space_from_domain(pw_space);
5332         pw_space = isl_space_add_dims(pw_space, isl_dim_out, 1);
5333         paff = isl_pw_aff_alloc_size(pw_space, var->n_row);
5334
5335         for (i = 0; i < var->n_row; ++i) {
5336                 isl_pw_aff *paff_i;
5337
5338                 aff = isl_aff_alloc(isl_local_space_copy(ls));
5339                 bset = isl_basic_set_alloc_space(isl_space_copy(space), 0,
5340                                                0, var->n_row - 1);
5341                 if (!aff || !bset)
5342                         goto error;
5343                 isl_int_set_si(aff->v->el[0], 1);
5344                 isl_seq_cpy(aff->v->el + 1, var->row[i], var->n_col);
5345                 isl_int_set_si(aff->v->el[1 + var->n_col], 0);
5346                 bset = select_minimum(bset, var, i);
5347                 paff_i = isl_pw_aff_alloc(isl_set_from_basic_set(bset), aff);
5348                 paff = isl_pw_aff_add_disjoint(paff, paff_i);
5349         }
5350
5351         isl_local_space_free(ls);
5352         isl_space_free(space);
5353         isl_mat_free(var);
5354         return paff;
5355 error:
5356         isl_aff_free(aff);
5357         isl_basic_set_free(bset);
5358         isl_pw_aff_free(paff);
5359         isl_local_space_free(ls);
5360         isl_space_free(space);
5361         isl_mat_free(var);
5362         return NULL;
5363 }
5364
5365 /* Given a piecewise multi-affine expression of which the last input variable
5366  * is the minimum of the bounds in "cst", plug in the value of the minimum.
5367  * This minimum expression is given in "min_expr_pa".
5368  * The set "min_expr" contains the same information, but in the form of a set.
5369  * The variable is subsequently projected out.
5370  *
5371  * The implementation is similar to those of "split" and "split_domain".
5372  * If the variable appears in a given expression, then minimum expression
5373  * is plugged in.  Otherwise, if the variable appears in the constraints
5374  * and a split is required, then the domain is split.  Otherwise, no split
5375  * is performed.
5376  */
5377 static __isl_give isl_pw_multi_aff *split_domain_pma(
5378         __isl_take isl_pw_multi_aff *opt, __isl_take isl_pw_aff *min_expr_pa,
5379         __isl_take isl_set *min_expr, __isl_take isl_mat *cst)
5380 {
5381         int n_in;
5382         int i;
5383         isl_space *space;
5384         isl_pw_multi_aff *res;
5385
5386         if (!opt || !min_expr || !cst)
5387                 goto error;
5388
5389         n_in = isl_pw_multi_aff_dim(opt, isl_dim_in);
5390         space = isl_pw_multi_aff_get_space(opt);
5391         space = isl_space_drop_dims(space, isl_dim_in, n_in - 1, 1);
5392         res = isl_pw_multi_aff_empty(space);
5393
5394         for (i = 0; i < opt->n; ++i) {
5395                 isl_pw_multi_aff *pma;
5396
5397                 pma = isl_pw_multi_aff_alloc(isl_set_copy(opt->p[i].set),
5398                                          isl_multi_aff_copy(opt->p[i].maff));
5399                 if (need_substitution(opt->p[i].maff))
5400                         pma = isl_pw_multi_aff_substitute(pma,
5401                                         isl_dim_in, n_in - 1, min_expr_pa);
5402                 else if (need_split_set(opt->p[i].set, cst))
5403                         pma = isl_pw_multi_aff_intersect_domain(pma,
5404                                                        isl_set_copy(min_expr));
5405                 pma = isl_pw_multi_aff_project_out(pma,
5406                                                     isl_dim_in, n_in - 1, 1);
5407
5408                 res = isl_pw_multi_aff_add_disjoint(res, pma);
5409         }
5410
5411         isl_pw_multi_aff_free(opt);
5412         isl_pw_aff_free(min_expr_pa);
5413         isl_set_free(min_expr);
5414         isl_mat_free(cst);
5415         return res;
5416 error:
5417         isl_pw_multi_aff_free(opt);
5418         isl_pw_aff_free(min_expr_pa);
5419         isl_set_free(min_expr);
5420         isl_mat_free(cst);
5421         return NULL;
5422 }
5423
5424 static __isl_give isl_pw_multi_aff *basic_map_partial_lexopt_pma(
5425         __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
5426         __isl_give isl_set **empty, int max);
5427
5428 /* This function is called from basic_map_partial_lexopt_symm.
5429  * The last variable of "bmap" and "dom" corresponds to the minimum
5430  * of the bounds in "cst".  "map_space" is the space of the original
5431  * input relation (of basic_map_partial_lexopt_symm) and "set_space"
5432  * is the space of the original domain.
5433  *
5434  * We recursively call basic_map_partial_lexopt and then plug in
5435  * the definition of the minimum in the result.
5436  */
5437 static __isl_give union isl_lex_res basic_map_partial_lexopt_symm_pma_core(
5438         __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
5439         __isl_give isl_set **empty, int max, __isl_take isl_mat *cst,
5440         __isl_take isl_space *map_space, __isl_take isl_space *set_space)
5441 {
5442         isl_pw_multi_aff *opt;
5443         isl_pw_aff *min_expr_pa;
5444         isl_set *min_expr;
5445         union isl_lex_res res;
5446
5447         min_expr = set_minimum(isl_basic_set_get_space(dom), isl_mat_copy(cst));
5448         min_expr_pa = set_minimum_pa(isl_basic_set_get_space(dom),
5449                                         isl_mat_copy(cst));
5450
5451         opt = basic_map_partial_lexopt_pma(bmap, dom, empty, max);
5452
5453         if (empty) {
5454                 *empty = split(*empty,
5455                                isl_set_copy(min_expr), isl_mat_copy(cst));
5456                 *empty = isl_set_reset_space(*empty, set_space);
5457         }
5458
5459         opt = split_domain_pma(opt, min_expr_pa, min_expr, cst);
5460         opt = isl_pw_multi_aff_reset_space(opt, map_space);
5461
5462         res.pma = opt;
5463         return res;
5464 }
5465
5466 static __isl_give isl_pw_multi_aff *basic_map_partial_lexopt_symm_pma(
5467         __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
5468         __isl_give isl_set **empty, int max, int first, int second)
5469 {
5470         return basic_map_partial_lexopt_symm(bmap, dom, empty, max,
5471                     first, second, &basic_map_partial_lexopt_symm_pma_core).pma;
5472 }
5473
5474 /* Recursive part of isl_basic_map_partial_lexopt_pw_multi_aff, after detecting
5475  * equalities and removing redundant constraints.
5476  *
5477  * We first check if there are any parallel constraints (left).
5478  * If not, we are in the base case.
5479  * If there are parallel constraints, we replace them by a single
5480  * constraint in basic_map_partial_lexopt_symm_pma and then call
5481  * this function recursively to look for more parallel constraints.
5482  */
5483 static __isl_give isl_pw_multi_aff *basic_map_partial_lexopt_pma(
5484         __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
5485         __isl_give isl_set **empty, int max)
5486 {
5487         int par = 0;
5488         int first, second;
5489
5490         if (!bmap)
5491                 goto error;
5492
5493         if (bmap->ctx->opt->pip_symmetry)
5494                 par = parallel_constraints(bmap, &first, &second);
5495         if (par < 0)
5496                 goto error;
5497         if (!par)
5498                 return basic_map_partial_lexopt_base_pma(bmap, dom, empty, max);
5499
5500         return basic_map_partial_lexopt_symm_pma(bmap, dom, empty, max,
5501                                                  first, second);
5502 error:
5503         isl_basic_set_free(dom);
5504         isl_basic_map_free(bmap);
5505         return NULL;
5506 }
5507
5508 /* Compute the lexicographic minimum (or maximum if "max" is set)
5509  * of "bmap" over the domain "dom" and return the result as a piecewise
5510  * multi-affine expression.
5511  * If "empty" is not NULL, then *empty is assigned a set that
5512  * contains those parts of the domain where there is no solution.
5513  * If "bmap" is marked as rational (ISL_BASIC_MAP_RATIONAL),
5514  * then we compute the rational optimum.  Otherwise, we compute
5515  * the integral optimum.
5516  *
5517  * We perform some preprocessing.  As the PILP solver does not
5518  * handle implicit equalities very well, we first make sure all
5519  * the equalities are explicitly available.
5520  *
5521  * We also add context constraints to the basic map and remove
5522  * redundant constraints.  This is only needed because of the
5523  * way we handle simple symmetries.  In particular, we currently look
5524  * for symmetries on the constraints, before we set up the main tableau.
5525  * It is then no good to look for symmetries on possibly redundant constraints.
5526  */
5527 __isl_give isl_pw_multi_aff *isl_basic_map_partial_lexopt_pw_multi_aff(
5528         __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
5529         __isl_give isl_set **empty, int max)
5530 {
5531         if (empty)
5532                 *empty = NULL;
5533         if (!bmap || !dom)
5534                 goto error;
5535
5536         isl_assert(bmap->ctx,
5537             isl_basic_map_compatible_domain(bmap, dom), goto error);
5538
5539         if (isl_basic_set_dim(dom, isl_dim_all) == 0)
5540                 return basic_map_partial_lexopt_pma(bmap, dom, empty, max);
5541
5542         bmap = isl_basic_map_intersect_domain(bmap, isl_basic_set_copy(dom));
5543         bmap = isl_basic_map_detect_equalities(bmap);
5544         bmap = isl_basic_map_remove_redundancies(bmap);
5545
5546         return basic_map_partial_lexopt_pma(bmap, dom, empty, max);
5547 error:
5548         isl_basic_set_free(dom);
5549         isl_basic_map_free(bmap);
5550         return NULL;
5551 }