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