isl_tab_pip.c

   1 /*
   2  * Copyright 2008-2009 Katholieke Universiteit Leuven
   3  *
   4  * Use of this software is governed by the GNU LGPLv2.1 license
   5  *
   6  * Written by Sven Verdoolaege, K.U.Leuven, Departement
   7  * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
   8  */
   9
  10 #include "isl_map_private.h"
  11 #include "isl_seq.h"
  12 #include "isl_tab.h"
  13 #include "isl_sample.h"
  14 #include <isl_mat_private.h>
  15
  16 /*
  17  * The implementation of parametric integer linear programming in this file
  18  * was inspired by the paper "Parametric Integer Programming" and the
  19  * report "Solving systems of affine (in)equalities" by Paul Feautrier
  20  * (and others).
  21  *
  22  * The strategy used for obtaining a feasible solution is different
  23  * from the one used in isl_tab.c.  In particular, in isl_tab.c,
  24  * upon finding a constraint that is not yet satisfied, we pivot
  25  * in a row that increases the constant term of row holding the
  26  * constraint, making sure the sample solution remains feasible
  27  * for all the constraints it already satisfied.
  28  * Here, we always pivot in the row holding the constraint,
  29  * choosing a column that induces the lexicographically smallest
  30  * increment to the sample solution.
  31  *
  32  * By starting out from a sample value that is lexicographically
  33  * smaller than any integer point in the problem space, the first
  34  * feasible integer sample point we find will also be the lexicographically
  35  * smallest.  If all variables can be assumed to be non-negative,
  36  * then the initial sample value may be chosen equal to zero.
  37  * However, we will not make this assumption.  Instead, we apply
  38  * the "big parameter" trick.  Any variable x is then not directly
  39  * used in the tableau, but instead it its represented by another
  40  * variable x' = M + x, where M is an arbitrarily large (positive)
  41  * value.  x' is therefore always non-negative, whatever the value of x.
  42  * Taking as initial smaple value x' = 0 corresponds to x = -M,
  43  * which is always smaller than any possible value of x.
  44  *
  45  * The big parameter trick is used in the main tableau and
  46  * also in the context tableau if isl_context_lex is used.
  47  * In this case, each tableaus has its own big parameter.
  48  * Before doing any real work, we check if all the parameters
  49  * happen to be non-negative.  If so, we drop the column corresponding
  50  * to M from the initial context tableau.
  51  * If isl_context_gbr is used, then the big parameter trick is only
  52  * used in the main tableau.
  53  */
  54
  55 struct isl_context;
  56 struct isl_context_op {
  57         /* detect nonnegative parameters in context and mark them in tab */
  58         struct isl_tab *(*detect_nonnegative_parameters)(
  59                         struct isl_context *context, struct isl_tab *tab);
  60         /* return temporary reference to basic set representation of context */
  61         struct isl_basic_set *(*peek_basic_set)(struct isl_context *context);
  62         /* return temporary reference to tableau representation of context */
  63         struct isl_tab *(*peek_tab)(struct isl_context *context);
  64         /* add equality; check is 1 if eq may not be valid;
  65          * update is 1 if we may want to call ineq_sign on context later.
  66          */
  67         void (*add_eq)(struct isl_context *context, isl_int *eq,
  68                         int check, int update);
  69         /* add inequality; check is 1 if ineq may not be valid;
  70          * update is 1 if we may want to call ineq_sign on context later.
  71          */
  72         void (*add_ineq)(struct isl_context *context, isl_int *ineq,
  73                         int check, int update);
  74         /* check sign of ineq based on previous information.
  75          * strict is 1 if saturation should be treated as a positive sign.
  76          */
  77         enum isl_tab_row_sign (*ineq_sign)(struct isl_context *context,
  78                         isl_int *ineq, int strict);
  79         /* check if inequality maintains feasibility */
  80         int (*test_ineq)(struct isl_context *context, isl_int *ineq);
  81         /* return index of a div that corresponds to "div" */
  82         int (*get_div)(struct isl_context *context, struct isl_tab *tab,
  83                         struct isl_vec *div);
  84         /* add div "div" to context and return non-negativity */
  85         int (*add_div)(struct isl_context *context, struct isl_vec *div);
  86         int (*detect_equalities)(struct isl_context *context,
  87                         struct isl_tab *tab);
  88         /* return row index of "best" split */
  89         int (*best_split)(struct isl_context *context, struct isl_tab *tab);
  90         /* check if context has already been determined to be empty */
  91         int (*is_empty)(struct isl_context *context);
  92         /* check if context is still usable */
  93         int (*is_ok)(struct isl_context *context);
  94         /* save a copy/snapshot of context */
  95         void *(*save)(struct isl_context *context);
  96         /* restore saved context */
  97         void (*restore)(struct isl_context *context, void *);
  98         /* invalidate context */
  99         void (*invalidate)(struct isl_context *context);
 100         /* free context */
 101         void (*free)(struct isl_context *context);
 102 };
 103
 104 struct isl_context {
 105         struct isl_context_op *op;
 106 };
 107
 108 struct isl_context_lex {
 109         struct isl_context context;
 110         struct isl_tab *tab;
 111 };
 112
 113 struct isl_partial_sol {
 114         int level;
 115         struct isl_basic_set *dom;
 116         struct isl_mat *M;
 117
 118         struct isl_partial_sol *next;
 119 };
 120
 121 struct isl_sol;
 122 struct isl_sol_callback {
 123         struct isl_tab_callback callback;
 124         struct isl_sol *sol;
 125 };
 126
 127 /* isl_sol is an interface for constructing a solution to
 128  * a parametric integer linear programming problem.
 129  * Every time the algorithm reaches a state where a solution
 130  * can be read off from the tableau (including cases where the tableau
 131  * is empty), the function "add" is called on the isl_sol passed
 132  * to find_solutions_main.
 133  *
 134  * The context tableau is owned by isl_sol and is updated incrementally.
 135  *
 136  * There are currently two implementations of this interface,
 137  * isl_sol_map, which simply collects the solutions in an isl_map
 138  * and (optionally) the parts of the context where there is no solution
 139  * in an isl_set, and
 140  * isl_sol_for, which calls a user-defined function for each part of
 141  * the solution.
 142  */
 143 struct isl_sol {
 144         int error;
 145         int rational;
 146         int level;
 147         int max;
 148         int n_out;
 149         struct isl_context *context;
 150         struct isl_partial_sol *partial;
 151         void (*add)(struct isl_sol *sol,
 152                             struct isl_basic_set *dom, struct isl_mat *M);
 153         void (*add_empty)(struct isl_sol *sol, struct isl_basic_set *bset);
 154         void (*free)(struct isl_sol *sol);
 155         struct isl_sol_callback dec_level;
 156 };
 157
 158 static void sol_free(struct isl_sol *sol)
 159 {
 160         struct isl_partial_sol *partial, *next;
 161         if (!sol)
 162                 return;
 163         for (partial = sol->partial; partial; partial = next) {
 164                 next = partial->next;
 165                 isl_basic_set_free(partial->dom);
 166                 isl_mat_free(partial->M);
 167                 free(partial);
 168         }
 169         sol->free(sol);
 170 }
 171
 172 /* Push a partial solution represented by a domain and mapping M
 173  * onto the stack of partial solutions.
 174  */
 175 static void sol_push_sol(struct isl_sol *sol,
 176         struct isl_basic_set *dom, struct isl_mat *M)
 177 {
 178         struct isl_partial_sol *partial;
 179
 180         if (sol->error || !dom)
 181                 goto error;
 182
 183         partial = isl_alloc_type(dom->ctx, struct isl_partial_sol);
 184         if (!partial)
 185                 goto error;
 186
 187         partial->level = sol->level;
 188         partial->dom = dom;
 189         partial->M = M;
 190         partial->next = sol->partial;
 191
 192         sol->partial = partial;
 193
 194         return;
 195 error:
 196         isl_basic_set_free(dom);
 197         sol->error = 1;
 198 }
 199
 200 /* Pop one partial solution from the partial solution stack and
 201  * pass it on to sol->add or sol->add_empty.
 202  */
 203 static void sol_pop_one(struct isl_sol *sol)
 204 {
 205         struct isl_partial_sol *partial;
 206
 207         partial = sol->partial;
 208         sol->partial = partial->next;
 209
 210         if (partial->M)
 211                 sol->add(sol, partial->dom, partial->M);
 212         else
 213                 sol->add_empty(sol, partial->dom);
 214         free(partial);
 215 }
 216
 217 /* Return a fresh copy of the domain represented by the context tableau.
 218  */
 219 static struct isl_basic_set *sol_domain(struct isl_sol *sol)
 220 {
 221         struct isl_basic_set *bset;
 222
 223         if (sol->error)
 224                 return NULL;
 225
 226         bset = isl_basic_set_dup(sol->context->op->peek_basic_set(sol->context));
 227         bset = isl_basic_set_update_from_tab(bset,
 228                         sol->context->op->peek_tab(sol->context));
 229
 230         return bset;
 231 }
 232
 233 /* Check whether two partial solutions have the same mapping, where n_div
 234  * is the number of divs that the two partial solutions have in common.
 235  */
 236 static int same_solution(struct isl_partial_sol *s1, struct isl_partial_sol *s2,
 237         unsigned n_div)
 238 {
 239         int i;
 240         unsigned dim;
 241
 242         if (!s1->M != !s2->M)
 243                 return 0;
 244         if (!s1->M)
 245                 return 1;
 246
 247         dim = isl_basic_set_total_dim(s1->dom) - s1->dom->n_div;
 248
 249         for (i = 0; i < s1->M->n_row; ++i) {
 250                 if (isl_seq_first_non_zero(s1->M->row[i]+1+dim+n_div,
 251                                             s1->M->n_col-1-dim-n_div) != -1)
 252                         return 0;
 253                 if (isl_seq_first_non_zero(s2->M->row[i]+1+dim+n_div,
 254                                             s2->M->n_col-1-dim-n_div) != -1)
 255                         return 0;
 256                 if (!isl_seq_eq(s1->M->row[i], s2->M->row[i], 1+dim+n_div))
 257                         return 0;
 258         }
 259         return 1;
 260 }
 261
 262 /* Pop all solutions from the partial solution stack that were pushed onto
 263  * the stack at levels that are deeper than the current level.
 264  * If the two topmost elements on the stack have the same level
 265  * and represent the same solution, then their domains are combined.
 266  * This combined domain is the same as the current context domain
 267  * as sol_pop is called each time we move back to a higher level.
 268  */
 269 static void sol_pop(struct isl_sol *sol)
 270 {
 271         struct isl_partial_sol *partial;
 272         unsigned n_div;
 273
 274         if (sol->error)
 275                 return;
 276
 277         if (sol->level == 0) {
 278                 for (partial = sol->partial; partial; partial = sol->partial)
 279                         sol_pop_one(sol);
 280                 return;
 281         }
 282
 283         partial = sol->partial;
 284         if (!partial)
 285                 return;
 286
 287         if (partial->level <= sol->level)
 288                 return;
 289
 290         if (partial->next && partial->next->level == partial->level) {
 291                 n_div = isl_basic_set_dim(
 292                                 sol->context->op->peek_basic_set(sol->context),
 293                                 isl_dim_div);
 294
 295                 if (!same_solution(partial, partial->next, n_div)) {
 296                         sol_pop_one(sol);
 297                         sol_pop_one(sol);
 298                 } else {
 299                         struct isl_basic_set *bset;
 300
 301                         bset = sol_domain(sol);
 302
 303                         isl_basic_set_free(partial->next->dom);
 304                         partial->next->dom = bset;
 305                         partial->next->level = sol->level;
 306
 307                         sol->partial = partial->next;
 308                         isl_basic_set_free(partial->dom);
 309                         isl_mat_free(partial->M);
 310                         free(partial);
 311                 }
 312         } else
 313                 sol_pop_one(sol);
 314 }
 315
 316 static void sol_dec_level(struct isl_sol *sol)
 317 {
 318         if (sol->error)
 319                 return;
 320
 321         sol->level--;
 322
 323         sol_pop(sol);
 324 }
 325
 326 static int sol_dec_level_wrap(struct isl_tab_callback *cb)
 327 {
 328         struct isl_sol_callback *callback = (struct isl_sol_callback *)cb;
 329
 330         sol_dec_level(callback->sol);
 331
 332         return callback->sol->error ? -1 : 0;
 333 }
 334
 335 /* Move down to next level and push callback onto context tableau
 336  * to decrease the level again when it gets rolled back across
 337  * the current state.  That is, dec_level will be called with
 338  * the context tableau in the same state as it is when inc_level
 339  * is called.
 340  */
 341 static void sol_inc_level(struct isl_sol *sol)
 342 {
 343         struct isl_tab *tab;
 344
 345         if (sol->error)
 346                 return;
 347
 348         sol->level++;
 349         tab = sol->context->op->peek_tab(sol->context);
 350         if (isl_tab_push_callback(tab, &sol->dec_level.callback) < 0)
 351                 sol->error = 1;
 352 }
 353
 354 static void scale_rows(struct isl_mat *mat, isl_int m, int n_row)
 355 {
 356         int i;
 357
 358         if (isl_int_is_one(m))
 359                 return;
 360
 361         for (i = 0; i < n_row; ++i)
 362                 isl_seq_scale(mat->row[i], mat->row[i], m, mat->n_col);
 363 }
 364
 365 /* Add the solution identified by the tableau and the context tableau.
 366  *
 367  * The layout of the variables is as follows.
 368  *      tab->n_var is equal to the total number of variables in the input
 369  *                      map (including divs that were copied from the context)
 370  *                      + the number of extra divs constructed
 371  *      Of these, the first tab->n_param and the last tab->n_div variables
 372  *      correspond to the variables in the context, i.e.,
 373  *              tab->n_param + tab->n_div = context_tab->n_var
 374  *      tab->n_param is equal to the number of parameters and input
 375  *                      dimensions in the input map
 376  *      tab->n_div is equal to the number of divs in the context
 377  *
 378  * If there is no solution, then call add_empty with a basic set
 379  * that corresponds to the context tableau.  (If add_empty is NULL,
 380  * then do nothing).
 381  *
 382  * If there is a solution, then first construct a matrix that maps
 383  * all dimensions of the context to the output variables, i.e.,
 384  * the output dimensions in the input map.
 385  * The divs in the input map (if any) that do not correspond to any
 386  * div in the context do not appear in the solution.
 387  * The algorithm will make sure that they have an integer value,
 388  * but these values themselves are of no interest.
 389  * We have to be careful not to drop or rearrange any divs in the
 390  * context because that would change the meaning of the matrix.
 391  *
 392  * To extract the value of the output variables, it should be noted
 393  * that we always use a big parameter M in the main tableau and so
 394  * the variable stored in this tableau is not an output variable x itself, but
 395  *      x' = M + x (in case of minimization)
 396  * or
 397  *      x' = M - x (in case of maximization)
 398  * If x' appears in a column, then its optimal value is zero,
 399  * which means that the optimal value of x is an unbounded number
 400  * (-M for minimization and M for maximization).
 401  * We currently assume that the output dimensions in the original map
 402  * are bounded, so this cannot occur.
 403  * Similarly, when x' appears in a row, then the coefficient of M in that
 404  * row is necessarily 1.
 405  * If the row in the tableau represents
 406  *      d x' = c + d M + e(y)
 407  * then, in case of minimization, the corresponding row in the matrix
 408  * will be
 409  *      a c + a e(y)
 410  * with a d = m, the (updated) common denominator of the matrix.
 411  * In case of maximization, the row will be
 412  *      -a c - a e(y)
 413  */
 414 static void sol_add(struct isl_sol *sol, struct isl_tab *tab)
 415 {
 416         struct isl_basic_set *bset = NULL;
 417         struct isl_mat *mat = NULL;
 418         unsigned off;
 419         int row, i;
 420         isl_int m;
 421
 422         if (sol->error || !tab)
 423                 goto error;
 424
 425         if (tab->empty && !sol->add_empty)
 426                 return;
 427
 428         bset = sol_domain(sol);
 429
 430         if (tab->empty) {
 431                 sol_push_sol(sol, bset, NULL);
 432                 return;
 433         }
 434
 435         off = 2 + tab->M;
 436
 437         mat = isl_mat_alloc(tab->mat->ctx, 1 + sol->n_out,
 438                                             1 + tab->n_param + tab->n_div);
 439         if (!mat)
 440                 goto error;
 441
 442         isl_int_init(m);
 443
 444         isl_seq_clr(mat->row[0] + 1, mat->n_col - 1);
 445         isl_int_set_si(mat->row[0][0], 1);
 446         for (row = 0; row < sol->n_out; ++row) {
 447                 int i = tab->n_param + row;
 448                 int r, j;
 449
 450                 isl_seq_clr(mat->row[1 + row], mat->n_col);
 451                 if (!tab->var[i].is_row) {
 452                         if (tab->M)
 453                                 isl_die(mat->ctx, isl_error_invalid,
 454                                         "unbounded optimum", goto error2);
 455                         continue;
 456                 }
 457
 458                 r = tab->var[i].index;
 459                 if (tab->M &&
 460                     isl_int_ne(tab->mat->row[r][2], tab->mat->row[r][0]))
 461                         isl_die(mat->ctx, isl_error_invalid,
 462                                 "unbounded optimum", goto error2);
 463                 isl_int_gcd(m, mat->row[0][0], tab->mat->row[r][0]);
 464                 isl_int_divexact(m, tab->mat->row[r][0], m);
 465                 scale_rows(mat, m, 1 + row);
 466                 isl_int_divexact(m, mat->row[0][0], tab->mat->row[r][0]);
 467                 isl_int_mul(mat->row[1 + row][0], m, tab->mat->row[r][1]);
 468                 for (j = 0; j < tab->n_param; ++j) {
 469                         int col;
 470                         if (tab->var[j].is_row)
 471                                 continue;
 472                         col = tab->var[j].index;
 473                         isl_int_mul(mat->row[1 + row][1 + j], m,
 474                                     tab->mat->row[r][off + col]);
 475                 }
 476                 for (j = 0; j < tab->n_div; ++j) {
 477                         int col;
 478                         if (tab->var[tab->n_var - tab->n_div+j].is_row)
 479                                 continue;
 480                         col = tab->var[tab->n_var - tab->n_div+j].index;
 481                         isl_int_mul(mat->row[1 + row][1 + tab->n_param + j], m,
 482                                     tab->mat->row[r][off + col]);
 483                 }
 484                 if (sol->max)
 485                         isl_seq_neg(mat->row[1 + row], mat->row[1 + row],
 486                                     mat->n_col);
 487         }
 488
 489         isl_int_clear(m);
 490
 491         sol_push_sol(sol, bset, mat);
 492         return;
 493 error2:
 494         isl_int_clear(m);
 495 error:
 496         isl_basic_set_free(bset);
 497         isl_mat_free(mat);
 498         sol->error = 1;
 499 }
 500
 501 struct isl_sol_map {
 502         struct isl_sol  sol;
 503         struct isl_map  *map;
 504         struct isl_set  *empty;
 505 };
 506
 507 static void sol_map_free(struct isl_sol_map *sol_map)
 508 {
 509         if (!sol_map)
 510                 return;
 511         if (sol_map->sol.context)
 512                 sol_map->sol.context->op->free(sol_map->sol.context);
 513         isl_map_free(sol_map->map);
 514         isl_set_free(sol_map->empty);
 515         free(sol_map);
 516 }
 517
 518 static void sol_map_free_wrap(struct isl_sol *sol)
 519 {
 520         sol_map_free((struct isl_sol_map *)sol);
 521 }
 522
 523 /* This function is called for parts of the context where there is
 524  * no solution, with "bset" corresponding to the context tableau.
 525  * Simply add the basic set to the set "empty".
 526  */
 527 static void sol_map_add_empty(struct isl_sol_map *sol,
 528         struct isl_basic_set *bset)
 529 {
 530         if (!bset)
 531                 goto error;
 532         isl_assert(bset->ctx, sol->empty, goto error);
 533
 534         sol->empty = isl_set_grow(sol->empty, 1);
 535         bset = isl_basic_set_simplify(bset);
 536         bset = isl_basic_set_finalize(bset);
 537         sol->empty = isl_set_add_basic_set(sol->empty, isl_basic_set_copy(bset));
 538         if (!sol->empty)
 539                 goto error;
 540         isl_basic_set_free(bset);
 541         return;
 542 error:
 543         isl_basic_set_free(bset);
 544         sol->sol.error = 1;
 545 }
 546
 547 static void sol_map_add_empty_wrap(struct isl_sol *sol,
 548         struct isl_basic_set *bset)
 549 {
 550         sol_map_add_empty((struct isl_sol_map *)sol, bset);
 551 }
 552
 553 /* Add bset to sol's empty, but only if we are actually collecting
 554  * the empty set.
 555  */
 556 static void sol_map_add_empty_if_needed(struct isl_sol_map *sol,
 557         struct isl_basic_set *bset)
 558 {
 559         if (sol->empty)
 560                 sol_map_add_empty(sol, bset);
 561         else
 562                 isl_basic_set_free(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         isl_basic_set *context_bset;
 585         unsigned n_eq;
 586         unsigned n_ineq;
 587         unsigned nparam;
 588         unsigned total;
 589         unsigned n_div;
 590         unsigned n_out;
 591
 592         if (sol->sol.error || !dom || !M)
 593                 goto error;
 594
 595         n_out = sol->sol.n_out;
 596         n_eq = dom->n_eq + n_out;
 597         n_ineq = dom->n_ineq;
 598         n_div = dom->n_div;
 599         nparam = isl_basic_set_total_dim(dom) - n_div;
 600         total = isl_map_dim(sol->map, isl_dim_all);
 601         bmap = isl_basic_map_alloc_dim(isl_map_get_dim(sol->map),
 602                                         n_div, n_eq, 2 * n_div + n_ineq);
 603         if (!bmap)
 604                 goto error;
 605         if (sol->sol.rational)
 606                 ISL_F_SET(bmap, ISL_BASIC_MAP_RATIONAL);
 607         for (i = 0; i < dom->n_div; ++i) {
 608                 int k = isl_basic_map_alloc_div(bmap);
 609                 if (k < 0)
 610                         goto error;
 611                 isl_seq_cpy(bmap->div[k], dom->div[i], 1 + 1 + nparam);
 612                 isl_seq_clr(bmap->div[k] + 1 + 1 + nparam, total - nparam);
 613                 isl_seq_cpy(bmap->div[k] + 1 + 1 + total,
 614                             dom->div[i] + 1 + 1 + nparam, i);
 615         }
 616         for (i = 0; i < dom->n_eq; ++i) {
 617                 int k = isl_basic_map_alloc_equality(bmap);
 618                 if (k < 0)
 619                         goto error;
 620                 isl_seq_cpy(bmap->eq[k], dom->eq[i], 1 + nparam);
 621                 isl_seq_clr(bmap->eq[k] + 1 + nparam, total - nparam);
 622                 isl_seq_cpy(bmap->eq[k] + 1 + total,
 623                             dom->eq[i] + 1 + nparam, n_div);
 624         }
 625         for (i = 0; i < dom->n_ineq; ++i) {
 626                 int k = isl_basic_map_alloc_inequality(bmap);
 627                 if (k < 0)
 628                         goto error;
 629                 isl_seq_cpy(bmap->ineq[k], dom->ineq[i], 1 + nparam);
 630                 isl_seq_clr(bmap->ineq[k] + 1 + nparam, total - nparam);
 631                 isl_seq_cpy(bmap->ineq[k] + 1 + total,
 632                         dom->ineq[i] + 1 + nparam, n_div);
 633         }
 634         for (i = 0; i < M->n_row - 1; ++i) {
 635                 int k = isl_basic_map_alloc_equality(bmap);
 636                 if (k < 0)
 637                         goto error;
 638                 isl_seq_cpy(bmap->eq[k], M->row[1 + i], 1 + nparam);
 639                 isl_seq_clr(bmap->eq[k] + 1 + nparam, n_out);
 640                 isl_int_neg(bmap->eq[k][1 + nparam + i], M->row[0][0]);
 641                 isl_seq_cpy(bmap->eq[k] + 1 + nparam + n_out,
 642                             M->row[1 + i] + 1 + nparam, n_div);
 643         }
 644         bmap = isl_basic_map_simplify(bmap);
 645         bmap = isl_basic_map_finalize(bmap);
 646         sol->map = isl_map_grow(sol->map, 1);
 647         sol->map = isl_map_add_basic_map(sol->map, bmap);
 648         if (!sol->map)
 649                 goto error;
 650         isl_basic_set_free(dom);
 651         isl_mat_free(M);
 652         return;
 653 error:
 654         isl_basic_set_free(dom);
 655         isl_mat_free(M);
 656         isl_basic_map_free(bmap);
 657         sol->sol.error = 1;
 658 }
 659
 660 static void sol_map_add_wrap(struct isl_sol *sol,
 661         struct isl_basic_set *dom, struct isl_mat *M)
 662 {
 663         sol_map_add((struct isl_sol_map *)sol, dom, M);
 664 }
 665
 666
 667 /* Store the "parametric constant" of row "row" of tableau "tab" in "line",
 668  * i.e., the constant term and the coefficients of all variables that
 669  * appear in the context tableau.
 670  * Note that the coefficient of the big parameter M is NOT copied.
 671  * The context tableau may not have a big parameter and even when it
 672  * does, it is a different big parameter.
 673  */
 674 static void get_row_parameter_line(struct isl_tab *tab, int row, isl_int *line)
 675 {
 676         int i;
 677         unsigned off = 2 + tab->M;
 678
 679         isl_int_set(line[0], tab->mat->row[row][1]);
 680         for (i = 0; i < tab->n_param; ++i) {
 681                 if (tab->var[i].is_row)
 682                         isl_int_set_si(line[1 + i], 0);
 683                 else {
 684                         int col = tab->var[i].index;
 685                         isl_int_set(line[1 + i], tab->mat->row[row][off + col]);
 686                 }
 687         }
 688         for (i = 0; i < tab->n_div; ++i) {
 689                 if (tab->var[tab->n_var - tab->n_div + i].is_row)
 690                         isl_int_set_si(line[1 + tab->n_param + i], 0);
 691                 else {
 692                         int col = tab->var[tab->n_var - tab->n_div + i].index;
 693                         isl_int_set(line[1 + tab->n_param + i],
 694                                     tab->mat->row[row][off + col]);
 695                 }
 696         }
 697 }
 698
 699 /* Check if rows "row1" and "row2" have identical "parametric constants",
 700  * as explained above.
 701  * In this case, we also insist that the coefficients of the big parameter
 702  * be the same as the values of the constants will only be the same
 703  * if these coefficients are also the same.
 704  */
 705 static int identical_parameter_line(struct isl_tab *tab, int row1, int row2)
 706 {
 707         int i;
 708         unsigned off = 2 + tab->M;
 709
 710         if (isl_int_ne(tab->mat->row[row1][1], tab->mat->row[row2][1]))
 711                 return 0;
 712
 713         if (tab->M && isl_int_ne(tab->mat->row[row1][2],
 714                                  tab->mat->row[row2][2]))
 715                 return 0;
 716
 717         for (i = 0; i < tab->n_param + tab->n_div; ++i) {
 718                 int pos = i < tab->n_param ? i :
 719                         tab->n_var - tab->n_div + i - tab->n_param;
 720                 int col;
 721
 722                 if (tab->var[pos].is_row)
 723                         continue;
 724                 col = tab->var[pos].index;
 725                 if (isl_int_ne(tab->mat->row[row1][off + col],
 726                                tab->mat->row[row2][off + col]))
 727                         return 0;
 728         }
 729         return 1;
 730 }
 731
 732 /* Return an inequality that expresses that the "parametric constant"
 733  * should be non-negative.
 734  * This function is only called when the coefficient of the big parameter
 735  * is equal to zero.
 736  */
 737 static struct isl_vec *get_row_parameter_ineq(struct isl_tab *tab, int row)
 738 {
 739         struct isl_vec *ineq;
 740
 741         ineq = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_param + tab->n_div);
 742         if (!ineq)
 743                 return NULL;
 744
 745         get_row_parameter_line(tab, row, ineq->el);
 746         if (ineq)
 747                 ineq = isl_vec_normalize(ineq);
 748
 749         return ineq;
 750 }
 751
 752 /* Return a integer division for use in a parametric cut based on the given row.
 753  * In particular, let the parametric constant of the row be
 754  *
 755  *              \sum_i a_i y_i
 756  *
 757  * where y_0 = 1, but none of the y_i corresponds to the big parameter M.
 758  * The div returned is equal to
 759  *
 760  *              floor(\sum_i {-a_i} y_i) = floor((\sum_i (-a_i mod d) y_i)/d)
 761  */
 762 static struct isl_vec *get_row_parameter_div(struct isl_tab *tab, int row)
 763 {
 764         struct isl_vec *div;
 765
 766         div = isl_vec_alloc(tab->mat->ctx, 1 + 1 + tab->n_param + tab->n_div);
 767         if (!div)
 768                 return NULL;
 769
 770         isl_int_set(div->el[0], tab->mat->row[row][0]);
 771         get_row_parameter_line(tab, row, div->el + 1);
 772         div = isl_vec_normalize(div);
 773         isl_seq_neg(div->el + 1, div->el + 1, div->size - 1);
 774         isl_seq_fdiv_r(div->el + 1, div->el + 1, div->el[0], div->size - 1);
 775
 776         return div;
 777 }
 778
 779 /* Return a integer division for use in transferring an integrality constraint
 780  * to the context.
 781  * In particular, let the parametric constant of the row be
 782  *
 783  *              \sum_i a_i y_i
 784  *
 785  * where y_0 = 1, but none of the y_i corresponds to the big parameter M.
 786  * The the returned div is equal to
 787  *
 788  *              floor(\sum_i {a_i} y_i) = floor((\sum_i (a_i mod d) y_i)/d)
 789  */
 790 static struct isl_vec *get_row_split_div(struct isl_tab *tab, int row)
 791 {
 792         struct isl_vec *div;
 793
 794         div = isl_vec_alloc(tab->mat->ctx, 1 + 1 + tab->n_param + tab->n_div);
 795         if (!div)
 796                 return NULL;
 797
 798         isl_int_set(div->el[0], tab->mat->row[row][0]);
 799         get_row_parameter_line(tab, row, div->el + 1);
 800         div = isl_vec_normalize(div);
 801         isl_seq_fdiv_r(div->el + 1, div->el + 1, div->el[0], div->size - 1);
 802
 803         return div;
 804 }
 805
 806 /* Construct and return an inequality that expresses an upper bound
 807  * on the given div.
 808  * In particular, if the div is given by
 809  *
 810  *      d = floor(e/m)
 811  *
 812  * then the inequality expresses
 813  *
 814  *      m d <= e
 815  */
 816 static struct isl_vec *ineq_for_div(struct isl_basic_set *bset, unsigned div)
 817 {
 818         unsigned total;
 819         unsigned div_pos;
 820         struct isl_vec *ineq;
 821
 822         if (!bset)
 823                 return NULL;
 824
 825         total = isl_basic_set_total_dim(bset);
 826         div_pos = 1 + total - bset->n_div + div;
 827
 828         ineq = isl_vec_alloc(bset->ctx, 1 + total);
 829         if (!ineq)
 830                 return NULL;
 831
 832         isl_seq_cpy(ineq->el, bset->div[div] + 1, 1 + total);
 833         isl_int_neg(ineq->el[div_pos], bset->div[div][0]);
 834         return ineq;
 835 }
 836
 837 /* Given a row in the tableau and a div that was created
 838  * using get_row_split_div and that been constrained to equality, i.e.,
 839  *
 840  *              d = floor(\sum_i {a_i} y_i) = \sum_i {a_i} y_i
 841  *
 842  * replace the expression "\sum_i {a_i} y_i" in the row by d,
 843  * i.e., we subtract "\sum_i {a_i} y_i" and add 1 d.
 844  * The coefficients of the non-parameters in the tableau have been
 845  * verified to be integral.  We can therefore simply replace coefficient b
 846  * by floor(b).  For the coefficients of the parameters we have
 847  * floor(a_i) = a_i - {a_i}, while for the other coefficients, we have
 848  * floor(b) = b.
 849  */
 850 static struct isl_tab *set_row_cst_to_div(struct isl_tab *tab, int row, int div)
 851 {
 852         isl_seq_fdiv_q(tab->mat->row[row] + 1, tab->mat->row[row] + 1,
 853                         tab->mat->row[row][0], 1 + tab->M + tab->n_col);
 854
 855         isl_int_set_si(tab->mat->row[row][0], 1);
 856
 857         if (tab->var[tab->n_var - tab->n_div + div].is_row) {
 858                 int drow = tab->var[tab->n_var - tab->n_div + div].index;
 859
 860                 isl_assert(tab->mat->ctx,
 861                         isl_int_is_one(tab->mat->row[drow][0]), goto error);
 862                 isl_seq_combine(tab->mat->row[row] + 1,
 863                         tab->mat->ctx->one, tab->mat->row[row] + 1,
 864                         tab->mat->ctx->one, tab->mat->row[drow] + 1,
 865                         1 + tab->M + tab->n_col);
 866         } else {
 867                 int dcol = tab->var[tab->n_var - tab->n_div + div].index;
 868
 869                 isl_int_set_si(tab->mat->row[row][2 + tab->M + dcol], 1);
 870         }
 871
 872         return tab;
 873 error:
 874         isl_tab_free(tab);
 875         return NULL;
 876 }
 877
 878 /* Check if the (parametric) constant of the given row is obviously
 879  * negative, meaning that we don't need to consult the context tableau.
 880  * If there is a big parameter and its coefficient is non-zero,
 881  * then this coefficient determines the outcome.
 882  * Otherwise, we check whether the constant is negative and
 883  * all non-zero coefficients of parameters are negative and
 884  * belong to non-negative parameters.
 885  */
 886 static int is_obviously_neg(struct isl_tab *tab, int row)
 887 {
 888         int i;
 889         int col;
 890         unsigned off = 2 + tab->M;
 891
 892         if (tab->M) {
 893                 if (isl_int_is_pos(tab->mat->row[row][2]))
 894                         return 0;
 895                 if (isl_int_is_neg(tab->mat->row[row][2]))
 896                         return 1;
 897         }
 898
 899         if (isl_int_is_nonneg(tab->mat->row[row][1]))
 900                 return 0;
 901         for (i = 0; i < tab->n_param; ++i) {
 902                 /* Eliminated parameter */
 903                 if (tab->var[i].is_row)
 904                         continue;
 905                 col = tab->var[i].index;
 906                 if (isl_int_is_zero(tab->mat->row[row][off + col]))
 907                         continue;
 908                 if (!tab->var[i].is_nonneg)
 909                         return 0;
 910                 if (isl_int_is_pos(tab->mat->row[row][off + col]))
 911                         return 0;
 912         }
 913         for (i = 0; i < tab->n_div; ++i) {
 914                 if (tab->var[tab->n_var - tab->n_div + i].is_row)
 915                         continue;
 916                 col = tab->var[tab->n_var - tab->n_div + i].index;
 917                 if (isl_int_is_zero(tab->mat->row[row][off + col]))
 918                         continue;
 919                 if (!tab->var[tab->n_var - tab->n_div + i].is_nonneg)
 920                         return 0;
 921                 if (isl_int_is_pos(tab->mat->row[row][off + col]))
 922                         return 0;
 923         }
 924         return 1;
 925 }
 926
 927 /* Check if the (parametric) constant of the given row is obviously
 928  * non-negative, meaning that we don't need to consult the context tableau.
 929  * If there is a big parameter and its coefficient is non-zero,
 930  * then this coefficient determines the outcome.
 931  * Otherwise, we check whether the constant is non-negative and
 932  * all non-zero coefficients of parameters are positive and
 933  * belong to non-negative parameters.
 934  */
 935 static int is_obviously_nonneg(struct isl_tab *tab, int row)
 936 {
 937         int i;
 938         int col;
 939         unsigned off = 2 + tab->M;
 940
 941         if (tab->M) {
 942                 if (isl_int_is_pos(tab->mat->row[row][2]))
 943                         return 1;
 944                 if (isl_int_is_neg(tab->mat->row[row][2]))
 945                         return 0;
 946         }
 947
 948         if (isl_int_is_neg(tab->mat->row[row][1]))
 949                 return 0;
 950         for (i = 0; i < tab->n_param; ++i) {
 951                 /* Eliminated parameter */
 952                 if (tab->var[i].is_row)
 953                         continue;
 954                 col = tab->var[i].index;
 955                 if (isl_int_is_zero(tab->mat->row[row][off + col]))
 956                         continue;
 957                 if (!tab->var[i].is_nonneg)
 958                         return 0;
 959                 if (isl_int_is_neg(tab->mat->row[row][off + col]))
 960                         return 0;
 961         }
 962         for (i = 0; i < tab->n_div; ++i) {
 963                 if (tab->var[tab->n_var - tab->n_div + i].is_row)
 964                         continue;
 965                 col = tab->var[tab->n_var - tab->n_div + i].index;
 966                 if (isl_int_is_zero(tab->mat->row[row][off + col]))
 967                         continue;
 968                 if (!tab->var[tab->n_var - tab->n_div + i].is_nonneg)
 969                         return 0;
 970                 if (isl_int_is_neg(tab->mat->row[row][off + col]))
 971                         return 0;
 972         }
 973         return 1;
 974 }
 975
 976 /* Given a row r and two columns, return the column that would
 977  * lead to the lexicographically smallest increment in the sample
 978  * solution when leaving the basis in favor of the row.
 979  * Pivoting with column c will increment the sample value by a non-negative
 980  * constant times a_{V,c}/a_{r,c}, with a_{V,c} the elements of column c
 981  * corresponding to the non-parametric variables.
 982  * If variable v appears in a column c_v, the a_{v,c} = 1 iff c = c_v,
 983  * with all other entries in this virtual row equal to zero.
 984  * If variable v appears in a row, then a_{v,c} is the element in column c
 985  * of that row.
 986  *
 987  * Let v be the first variable with a_{v,c1}/a_{r,c1} != a_{v,c2}/a_{r,c2}.
 988  * Then if a_{v,c1}/a_{r,c1} < a_{v,c2}/a_{r,c2}, i.e.,
 989  * a_{v,c2} a_{r,c1} - a_{v,c1} a_{r,c2} > 0, c1 results in the minimal
 990  * increment.  Otherwise, it's c2.
 991  */
 992 static int lexmin_col_pair(struct isl_tab *tab,
 993         int row, int col1, int col2, isl_int tmp)
 994 {
 995         int i;
 996         isl_int *tr;
 997
 998         tr = tab->mat->row[row] + 2 + tab->M;
 999
1000         for (i = tab->n_param; i < tab->n_var - tab->n_div; ++i) {
1001                 int s1, s2;
1002                 isl_int *r;
1003
1004                 if (!tab->var[i].is_row) {
1005                         if (tab->var[i].index == col1)
1006                                 return col2;
1007                         if (tab->var[i].index == col2)
1008                                 return col1;
1009                         continue;
1010                 }
1011
1012                 if (tab->var[i].index == row)
1013                         continue;
1014
1015                 r = tab->mat->row[tab->var[i].index] + 2 + tab->M;
1016                 s1 = isl_int_sgn(r[col1]);
1017                 s2 = isl_int_sgn(r[col2]);
1018                 if (s1 == 0 && s2 == 0)
1019                         continue;
1020                 if (s1 < s2)
1021                         return col1;
1022                 if (s2 < s1)
1023                         return col2;
1024
1025                 isl_int_mul(tmp, r[col2], tr[col1]);
1026                 isl_int_submul(tmp, r[col1], tr[col2]);
1027                 if (isl_int_is_pos(tmp))
1028                         return col1;
1029                 if (isl_int_is_neg(tmp))
1030                         return col2;
1031         }
1032         return -1;
1033 }
1034
1035 /* Given a row in the tableau, find and return the column that would
1036  * result in the lexicographically smallest, but positive, increment
1037  * in the sample point.
1038  * If there is no such column, then return tab->n_col.
1039  * If anything goes wrong, return -1.
1040  */
1041 static int lexmin_pivot_col(struct isl_tab *tab, int row)
1042 {
1043         int j;
1044         int col = tab->n_col;
1045         isl_int *tr;
1046         isl_int tmp;
1047
1048         tr = tab->mat->row[row] + 2 + tab->M;
1049
1050         isl_int_init(tmp);
1051
1052         for (j = tab->n_dead; j < tab->n_col; ++j) {
1053                 if (tab->col_var[j] >= 0 &&
1054                     (tab->col_var[j] < tab->n_param  ||
1055                     tab->col_var[j] >= tab->n_var - tab->n_div))
1056                         continue;
1057
1058                 if (!isl_int_is_pos(tr[j]))
1059                         continue;
1060
1061                 if (col == tab->n_col)
1062                         col = j;
1063                 else
1064                         col = lexmin_col_pair(tab, row, col, j, tmp);
1065                 isl_assert(tab->mat->ctx, col >= 0, goto error);
1066         }
1067
1068         isl_int_clear(tmp);
1069         return col;
1070 error:
1071         isl_int_clear(tmp);
1072         return -1;
1073 }
1074
1075 /* Return the first known violated constraint, i.e., a non-negative
1076  * contraint that currently has an either obviously negative value
1077  * or a previously determined to be negative value.
1078  *
1079  * If any constraint has a negative coefficient for the big parameter,
1080  * if any, then we return one of these first.
1081  */
1082 static int first_neg(struct isl_tab *tab)
1083 {
1084         int row;
1085
1086         if (tab->M)
1087                 for (row = tab->n_redundant; row < tab->n_row; ++row) {
1088                         if (!isl_tab_var_from_row(tab, row)->is_nonneg)
1089                                 continue;
1090                         if (!isl_int_is_neg(tab->mat->row[row][2]))
1091                                 continue;
1092                         if (tab->row_sign)
1093                                 tab->row_sign[row] = isl_tab_row_neg;
1094                         return row;
1095                 }
1096         for (row = tab->n_redundant; row < tab->n_row; ++row) {
1097                 if (!isl_tab_var_from_row(tab, row)->is_nonneg)
1098                         continue;
1099                 if (tab->row_sign) {
1100                         if (tab->row_sign[row] == 0 &&
1101                             is_obviously_neg(tab, row))
1102                                 tab->row_sign[row] = isl_tab_row_neg;
1103                         if (tab->row_sign[row] != isl_tab_row_neg)
1104                                 continue;
1105                 } else if (!is_obviously_neg(tab, row))
1106                         continue;
1107                 return row;
1108         }
1109         return -1;
1110 }
1111
1112 /* Resolve all known or obviously violated constraints through pivoting.
1113  * In particular, as long as we can find any violated constraint, we
1114  * look for a pivoting column that would result in the lexicographicallly
1115  * smallest increment in the sample point.  If there is no such column
1116  * then the tableau is infeasible.
1117  */
1118 static struct isl_tab *restore_lexmin(struct isl_tab *tab) WARN_UNUSED;
1119 static struct isl_tab *restore_lexmin(struct isl_tab *tab)
1120 {
1121         int row, col;
1122
1123         if (!tab)
1124                 return NULL;
1125         if (tab->empty)
1126                 return tab;
1127         while ((row = first_neg(tab)) != -1) {
1128                 col = lexmin_pivot_col(tab, row);
1129                 if (col >= tab->n_col) {
1130                         if (isl_tab_mark_empty(tab) < 0)
1131                                 goto error;
1132                         return tab;
1133                 }
1134                 if (col < 0)
1135                         goto error;
1136                 if (isl_tab_pivot(tab, row, col) < 0)
1137                         goto error;
1138         }
1139         return tab;
1140 error:
1141         isl_tab_free(tab);
1142         return NULL;
1143 }
1144
1145 /* Given a row that represents an equality, look for an appropriate
1146  * pivoting column.
1147  * In particular, if there are any non-zero coefficients among
1148  * the non-parameter variables, then we take the last of these
1149  * variables.  Eliminating this variable in terms of the other
1150  * variables and/or parameters does not influence the property
1151  * that all column in the initial tableau are lexicographically
1152  * positive.  The row corresponding to the eliminated variable
1153  * will only have non-zero entries below the diagonal of the
1154  * initial tableau.  That is, we transform
1155  *
1156  *              I                               I
1157  *                1             into            a
1158  *                  I                             I
1159  *
1160  * If there is no such non-parameter variable, then we are dealing with
1161  * pure parameter equality and we pick any parameter with coefficient 1 or -1
1162  * for elimination.  This will ensure that the eliminated parameter
1163  * always has an integer value whenever all the other parameters are integral.
1164  * If there is no such parameter then we return -1.
1165  */
1166 static int last_var_col_or_int_par_col(struct isl_tab *tab, int row)
1167 {
1168         unsigned off = 2 + tab->M;
1169         int i;
1170
1171         for (i = tab->n_var - tab->n_div - 1; i >= 0 && i >= tab->n_param; --i) {
1172                 int col;
1173                 if (tab->var[i].is_row)
1174                         continue;
1175                 col = tab->var[i].index;
1176                 if (col <= tab->n_dead)
1177                         continue;
1178                 if (!isl_int_is_zero(tab->mat->row[row][off + col]))
1179                         return col;
1180         }
1181         for (i = tab->n_dead; i < tab->n_col; ++i) {
1182                 if (isl_int_is_one(tab->mat->row[row][off + i]))
1183                         return i;
1184                 if (isl_int_is_negone(tab->mat->row[row][off + i]))
1185                         return i;
1186         }
1187         return -1;
1188 }
1189
1190 /* Add an equality that is known to be valid to the tableau.
1191  * We first check if we can eliminate a variable or a parameter.
1192  * If not, we add the equality as two inequalities.
1193  * In this case, the equality was a pure parameter equality and there
1194  * is no need to resolve any constraint violations.
1195  */
1196 static struct isl_tab *add_lexmin_valid_eq(struct isl_tab *tab, isl_int *eq)
1197 {
1198         int i;
1199         int r;
1200
1201         if (!tab)
1202                 return NULL;
1203         r = isl_tab_add_row(tab, eq);
1204         if (r < 0)
1205                 goto error;
1206
1207         r = tab->con[r].index;
1208         i = last_var_col_or_int_par_col(tab, r);
1209         if (i < 0) {
1210                 tab->con[r].is_nonneg = 1;
1211                 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1212                         goto error;
1213                 isl_seq_neg(eq, eq, 1 + tab->n_var);
1214                 r = isl_tab_add_row(tab, eq);
1215                 if (r < 0)
1216                         goto error;
1217                 tab->con[r].is_nonneg = 1;
1218                 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1219                         goto error;
1220         } else {
1221                 if (isl_tab_pivot(tab, r, i) < 0)
1222                         goto error;
1223                 if (isl_tab_kill_col(tab, i) < 0)
1224                         goto error;
1225                 tab->n_eq++;
1226
1227                 tab = restore_lexmin(tab);
1228         }
1229
1230         return tab;
1231 error:
1232         isl_tab_free(tab);
1233         return NULL;
1234 }
1235
1236 /* Check if the given row is a pure constant.
1237  */
1238 static int is_constant(struct isl_tab *tab, int row)
1239 {
1240         unsigned off = 2 + tab->M;
1241
1242         return isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
1243                                         tab->n_col - tab->n_dead) == -1;
1244 }
1245
1246 /* Add an equality that may or may not be valid to the tableau.
1247  * If the resulting row is a pure constant, then it must be zero.
1248  * Otherwise, the resulting tableau is empty.
1249  *
1250  * If the row is not a pure constant, then we add two inequalities,
1251  * each time checking that they can be satisfied.
1252  * In the end we try to use one of the two constraints to eliminate
1253  * a column.
1254  */
1255 static struct isl_tab *add_lexmin_eq(struct isl_tab *tab, isl_int *eq) WARN_UNUSED;
1256 static struct isl_tab *add_lexmin_eq(struct isl_tab *tab, isl_int *eq)
1257 {
1258         int r1, r2;
1259         int row;
1260         struct isl_tab_undo *snap;
1261
1262         if (!tab)
1263                 return NULL;
1264         snap = isl_tab_snap(tab);
1265         r1 = isl_tab_add_row(tab, eq);
1266         if (r1 < 0)
1267                 goto error;
1268         tab->con[r1].is_nonneg = 1;
1269         if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r1]) < 0)
1270                 goto error;
1271
1272         row = tab->con[r1].index;
1273         if (is_constant(tab, row)) {
1274                 if (!isl_int_is_zero(tab->mat->row[row][1]) ||
1275                     (tab->M && !isl_int_is_zero(tab->mat->row[row][2]))) {
1276                         if (isl_tab_mark_empty(tab) < 0)
1277                                 goto error;
1278                         return tab;
1279                 }
1280                 if (isl_tab_rollback(tab, snap) < 0)
1281                         goto error;
1282                 return tab;
1283         }
1284
1285         tab = restore_lexmin(tab);
1286         if (!tab || tab->empty)
1287                 return tab;
1288
1289         isl_seq_neg(eq, eq, 1 + tab->n_var);
1290
1291         r2 = isl_tab_add_row(tab, eq);
1292         if (r2 < 0)
1293                 goto error;
1294         tab->con[r2].is_nonneg = 1;
1295         if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r2]) < 0)
1296                 goto error;
1297
1298         tab = restore_lexmin(tab);
1299         if (!tab || tab->empty)
1300                 return tab;
1301
1302         if (!tab->con[r1].is_row) {
1303                 if (isl_tab_kill_col(tab, tab->con[r1].index) < 0)
1304                         goto error;
1305         } else if (!tab->con[r2].is_row) {
1306                 if (isl_tab_kill_col(tab, tab->con[r2].index) < 0)
1307                         goto error;
1308         } else if (isl_int_is_zero(tab->mat->row[tab->con[r1].index][1])) {
1309                 unsigned off = 2 + tab->M;
1310                 int i;
1311                 int row = tab->con[r1].index;
1312                 i = isl_seq_first_non_zero(tab->mat->row[row]+off+tab->n_dead,
1313                                                 tab->n_col - tab->n_dead);
1314                 if (i != -1) {
1315                         if (isl_tab_pivot(tab, row, tab->n_dead + i) < 0)
1316                                 goto error;
1317                         if (isl_tab_kill_col(tab, tab->n_dead + i) < 0)
1318                                 goto error;
1319                 }
1320         }
1321
1322         if (tab->bmap) {
1323                 tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);
1324                 if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1325                         goto error;
1326                 isl_seq_neg(eq, eq, 1 + tab->n_var);
1327                 tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);
1328                 isl_seq_neg(eq, eq, 1 + tab->n_var);
1329                 if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1330                         goto error;
1331                 if (!tab->bmap)
1332                         goto error;
1333         }
1334
1335         return tab;
1336 error:
1337         isl_tab_free(tab);
1338         return NULL;
1339 }
1340
1341 /* Add an inequality to the tableau, resolving violations using
1342  * restore_lexmin.
1343  */
1344 static struct isl_tab *add_lexmin_ineq(struct isl_tab *tab, isl_int *ineq)
1345 {
1346         int r;
1347
1348         if (!tab)
1349                 return NULL;
1350         if (tab->bmap) {
1351                 tab->bmap = isl_basic_map_add_ineq(tab->bmap, ineq);
1352                 if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1353                         goto error;
1354                 if (!tab->bmap)
1355                         goto error;
1356         }
1357         r = isl_tab_add_row(tab, ineq);
1358         if (r < 0)
1359                 goto error;
1360         tab->con[r].is_nonneg = 1;
1361         if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1362                 goto error;
1363         if (isl_tab_row_is_redundant(tab, tab->con[r].index)) {
1364                 if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0)
1365                         goto error;
1366                 return tab;
1367         }
1368
1369         tab = restore_lexmin(tab);
1370         if (tab && !tab->empty && tab->con[r].is_row &&
1371                  isl_tab_row_is_redundant(tab, tab->con[r].index))
1372                 if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0)
1373                         goto error;
1374         return tab;
1375 error:
1376         isl_tab_free(tab);
1377         return NULL;
1378 }
1379
1380 /* Check if the coefficients of the parameters are all integral.
1381  */
1382 static int integer_parameter(struct isl_tab *tab, int row)
1383 {
1384         int i;
1385         int col;
1386         unsigned off = 2 + tab->M;
1387
1388         for (i = 0; i < tab->n_param; ++i) {
1389                 /* Eliminated parameter */
1390                 if (tab->var[i].is_row)
1391                         continue;
1392                 col = tab->var[i].index;
1393                 if (!isl_int_is_divisible_by(tab->mat->row[row][off + col],
1394                                                 tab->mat->row[row][0]))
1395                         return 0;
1396         }
1397         for (i = 0; i < tab->n_div; ++i) {
1398                 if (tab->var[tab->n_var - tab->n_div + i].is_row)
1399                         continue;
1400                 col = tab->var[tab->n_var - tab->n_div + i].index;
1401                 if (!isl_int_is_divisible_by(tab->mat->row[row][off + col],
1402                                                 tab->mat->row[row][0]))
1403                         return 0;
1404         }
1405         return 1;
1406 }
1407
1408 /* Check if the coefficients of the non-parameter variables are all integral.
1409  */
1410 static int integer_variable(struct isl_tab *tab, int row)
1411 {
1412         int i;
1413         unsigned off = 2 + tab->M;
1414
1415         for (i = tab->n_dead; i < tab->n_col; ++i) {
1416                 if (tab->col_var[i] >= 0 &&
1417                     (tab->col_var[i] < tab->n_param ||
1418                      tab->col_var[i] >= tab->n_var - tab->n_div))
1419                         continue;
1420                 if (!isl_int_is_divisible_by(tab->mat->row[row][off + i],
1421                                                 tab->mat->row[row][0]))
1422                         return 0;
1423         }
1424         return 1;
1425 }
1426
1427 /* Check if the constant term is integral.
1428  */
1429 static int integer_constant(struct isl_tab *tab, int row)
1430 {
1431         return isl_int_is_divisible_by(tab->mat->row[row][1],
1432                                         tab->mat->row[row][0]);
1433 }
1434
1435 #define I_CST   1 << 0
1436 #define I_PAR   1 << 1
1437 #define I_VAR   1 << 2
1438
1439 /* Check for next (non-parameter) variable after "var" (first if var == -1)
1440  * that is non-integer and therefore requires a cut and return
1441  * the index of the variable.
1442  * For parametric tableaus, there are three parts in a row,
1443  * the constant, the coefficients of the parameters and the rest.
1444  * For each part, we check whether the coefficients in that part
1445  * are all integral and if so, set the corresponding flag in *f.
1446  * If the constant and the parameter part are integral, then the
1447  * current sample value is integral and no cut is required
1448  * (irrespective of whether the variable part is integral).
1449  */
1450 static int next_non_integer_var(struct isl_tab *tab, int var, int *f)
1451 {
1452         var = var < 0 ? tab->n_param : var + 1;
1453
1454         for (; var < tab->n_var - tab->n_div; ++var) {
1455                 int flags = 0;
1456                 int row;
1457                 if (!tab->var[var].is_row)
1458                         continue;
1459                 row = tab->var[var].index;
1460                 if (integer_constant(tab, row))
1461                         ISL_FL_SET(flags, I_CST);
1462                 if (integer_parameter(tab, row))
1463                         ISL_FL_SET(flags, I_PAR);
1464                 if (ISL_FL_ISSET(flags, I_CST) && ISL_FL_ISSET(flags, I_PAR))
1465                         continue;
1466                 if (integer_variable(tab, row))
1467                         ISL_FL_SET(flags, I_VAR);
1468                 *f = flags;
1469                 return var;
1470         }
1471         return -1;
1472 }
1473
1474 /* Check for first (non-parameter) variable that is non-integer and
1475  * therefore requires a cut and return the corresponding row.
1476  * For parametric tableaus, there are three parts in a row,
1477  * the constant, the coefficients of the parameters and the rest.
1478  * For each part, we check whether the coefficients in that part
1479  * are all integral and if so, set the corresponding flag in *f.
1480  * If the constant and the parameter part are integral, then the
1481  * current sample value is integral and no cut is required
1482  * (irrespective of whether the variable part is integral).
1483  */
1484 static int first_non_integer_row(struct isl_tab *tab, int *f)
1485 {
1486         int var = next_non_integer_var(tab, -1, f);
1487
1488         return var < 0 ? -1 : tab->var[var].index;
1489 }
1490
1491 /* Add a (non-parametric) cut to cut away the non-integral sample
1492  * value of the given row.
1493  *
1494  * If the row is given by
1495  *
1496  *      m r = f + \sum_i a_i y_i
1497  *
1498  * then the cut is
1499  *
1500  *      c = - {-f/m} + \sum_i {a_i/m} y_i >= 0
1501  *
1502  * The big parameter, if any, is ignored, since it is assumed to be big
1503  * enough to be divisible by any integer.
1504  * If the tableau is actually a parametric tableau, then this function
1505  * is only called when all coefficients of the parameters are integral.
1506  * The cut therefore has zero coefficients for the parameters.
1507  *
1508  * The current value is known to be negative, so row_sign, if it
1509  * exists, is set accordingly.
1510  *
1511  * Return the row of the cut or -1.
1512  */
1513 static int add_cut(struct isl_tab *tab, int row)
1514 {
1515         int i;
1516         int r;
1517         isl_int *r_row;
1518         unsigned off = 2 + tab->M;
1519
1520         if (isl_tab_extend_cons(tab, 1) < 0)
1521                 return -1;
1522         r = isl_tab_allocate_con(tab);
1523         if (r < 0)
1524                 return -1;
1525
1526         r_row = tab->mat->row[tab->con[r].index];
1527         isl_int_set(r_row[0], tab->mat->row[row][0]);
1528         isl_int_neg(r_row[1], tab->mat->row[row][1]);
1529         isl_int_fdiv_r(r_row[1], r_row[1], tab->mat->row[row][0]);
1530         isl_int_neg(r_row[1], r_row[1]);
1531         if (tab->M)
1532                 isl_int_set_si(r_row[2], 0);
1533         for (i = 0; i < tab->n_col; ++i)
1534                 isl_int_fdiv_r(r_row[off + i],
1535                         tab->mat->row[row][off + i], tab->mat->row[row][0]);
1536
1537         tab->con[r].is_nonneg = 1;
1538         if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1539                 return -1;
1540         if (tab->row_sign)
1541                 tab->row_sign[tab->con[r].index] = isl_tab_row_neg;
1542
1543         return tab->con[r].index;
1544 }
1545
1546 /* Given a non-parametric tableau, add cuts until an integer
1547  * sample point is obtained or until the tableau is determined
1548  * to be integer infeasible.
1549  * As long as there is any non-integer value in the sample point,
1550  * we add appropriate cuts, if possible, for each of these
1551  * non-integer values and then resolve the violated
1552  * cut constraints using restore_lexmin.
1553  * If one of the corresponding rows is equal to an integral
1554  * combination of variables/constraints plus a non-integral constant,
1555  * then there is no way to obtain an integer point and we return
1556  * a tableau that is marked empty.
1557  */
1558 static struct isl_tab *cut_to_integer_lexmin(struct isl_tab *tab)
1559 {
1560         int var;
1561         int row;
1562         int flags;
1563
1564         if (!tab)
1565                 return NULL;
1566         if (tab->empty)
1567                 return tab;
1568
1569         while ((var = next_non_integer_var(tab, -1, &flags)) != -1) {
1570                 do {
1571                         if (ISL_FL_ISSET(flags, I_VAR)) {
1572                                 if (isl_tab_mark_empty(tab) < 0)
1573                                         goto error;
1574                                 return tab;
1575                         }
1576                         row = tab->var[var].index;
1577                         row = add_cut(tab, row);
1578                         if (row < 0)
1579                                 goto error;
1580                 } while ((var = next_non_integer_var(tab, var, &flags)) != -1);
1581                 tab = restore_lexmin(tab);
1582                 if (!tab || tab->empty)
1583                         break;
1584         }
1585         return tab;
1586 error:
1587         isl_tab_free(tab);
1588         return NULL;
1589 }
1590
1591 /* Check whether all the currently active samples also satisfy the inequality
1592  * "ineq" (treated as an equality if eq is set).
1593  * Remove those samples that do not.
1594  */
1595 static struct isl_tab *check_samples(struct isl_tab *tab, isl_int *ineq, int eq)
1596 {
1597         int i;
1598         isl_int v;
1599
1600         if (!tab)
1601                 return NULL;
1602
1603         isl_assert(tab->mat->ctx, tab->bmap, goto error);
1604         isl_assert(tab->mat->ctx, tab->samples, goto error);
1605         isl_assert(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var, goto error);
1606
1607         isl_int_init(v);
1608         for (i = tab->n_outside; i < tab->n_sample; ++i) {
1609                 int sgn;
1610                 isl_seq_inner_product(ineq, tab->samples->row[i],
1611                                         1 + tab->n_var, &v);
1612                 sgn = isl_int_sgn(v);
1613                 if (eq ? (sgn == 0) : (sgn >= 0))
1614                         continue;
1615                 tab = isl_tab_drop_sample(tab, i);
1616                 if (!tab)
1617                         break;
1618         }
1619         isl_int_clear(v);
1620
1621         return tab;
1622 error:
1623         isl_tab_free(tab);
1624         return NULL;
1625 }
1626
1627 /* Check whether the sample value of the tableau is finite,
1628  * i.e., either the tableau does not use a big parameter, or
1629  * all values of the variables are equal to the big parameter plus
1630  * some constant.  This constant is the actual sample value.
1631  */
1632 static int sample_is_finite(struct isl_tab *tab)
1633 {
1634         int i;
1635
1636         if (!tab->M)
1637                 return 1;
1638
1639         for (i = 0; i < tab->n_var; ++i) {
1640                 int row;
1641                 if (!tab->var[i].is_row)
1642                         return 0;
1643                 row = tab->var[i].index;
1644                 if (isl_int_ne(tab->mat->row[row][0], tab->mat->row[row][2]))
1645                         return 0;
1646         }
1647         return 1;
1648 }
1649
1650 /* Check if the context tableau of sol has any integer points.
1651  * Leave tab in empty state if no integer point can be found.
1652  * If an integer point can be found and if moreover it is finite,
1653  * then it is added to the list of sample values.
1654  *
1655  * This function is only called when none of the currently active sample
1656  * values satisfies the most recently added constraint.
1657  */
1658 static struct isl_tab *check_integer_feasible(struct isl_tab *tab)
1659 {
1660         struct isl_tab_undo *snap;
1661         int feasible;
1662
1663         if (!tab)
1664                 return NULL;
1665
1666         snap = isl_tab_snap(tab);
1667         if (isl_tab_push_basis(tab) < 0)
1668                 goto error;
1669
1670         tab = cut_to_integer_lexmin(tab);
1671         if (!tab)
1672                 goto error;
1673
1674         if (!tab->empty && sample_is_finite(tab)) {
1675                 struct isl_vec *sample;
1676
1677                 sample = isl_tab_get_sample_value(tab);
1678
1679                 tab = isl_tab_add_sample(tab, sample);
1680         }
1681
1682         if (!tab->empty && isl_tab_rollback(tab, snap) < 0)
1683                 goto error;
1684
1685         return tab;
1686 error:
1687         isl_tab_free(tab);
1688         return NULL;
1689 }
1690
1691 /* Check if any of the currently active sample values satisfies
1692  * the inequality "ineq" (an equality if eq is set).
1693  */
1694 static int tab_has_valid_sample(struct isl_tab *tab, isl_int *ineq, int eq)
1695 {
1696         int i;
1697         isl_int v;
1698
1699         if (!tab)
1700                 return -1;
1701
1702         isl_assert(tab->mat->ctx, tab->bmap, return -1);
1703         isl_assert(tab->mat->ctx, tab->samples, return -1);
1704         isl_assert(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var, return -1);
1705
1706         isl_int_init(v);
1707         for (i = tab->n_outside; i < tab->n_sample; ++i) {
1708                 int sgn;
1709                 isl_seq_inner_product(ineq, tab->samples->row[i],
1710                                         1 + tab->n_var, &v);
1711                 sgn = isl_int_sgn(v);
1712                 if (eq ? (sgn == 0) : (sgn >= 0))
1713                         break;
1714         }
1715         isl_int_clear(v);
1716
1717         return i < tab->n_sample;
1718 }
1719
1720 /* Add a div specifed by "div" to the tableau "tab" and return
1721  * 1 if the div is obviously non-negative.
1722  */
1723 static int context_tab_add_div(struct isl_tab *tab, struct isl_vec *div,
1724         int (*add_ineq)(void *user, isl_int *), void *user)
1725 {
1726         int i;
1727         int r;
1728         struct isl_mat *samples;
1729         int nonneg;
1730
1731         r = isl_tab_add_div(tab, div, add_ineq, user);
1732         if (r < 0)
1733                 return -1;
1734         nonneg = tab->var[r].is_nonneg;
1735         tab->var[r].frozen = 1;
1736
1737         samples = isl_mat_extend(tab->samples,
1738                         tab->n_sample, 1 + tab->n_var);
1739         tab->samples = samples;
1740         if (!samples)
1741                 return -1;
1742         for (i = tab->n_outside; i < samples->n_row; ++i) {
1743                 isl_seq_inner_product(div->el + 1, samples->row[i],
1744                         div->size - 1, &samples->row[i][samples->n_col - 1]);
1745                 isl_int_fdiv_q(samples->row[i][samples->n_col - 1],
1746                                samples->row[i][samples->n_col - 1], div->el[0]);
1747         }
1748
1749         return nonneg;
1750 }
1751
1752 /* Add a div specified by "div" to both the main tableau and
1753  * the context tableau.  In case of the main tableau, we only
1754  * need to add an extra div.  In the context tableau, we also
1755  * need to express the meaning of the div.
1756  * Return the index of the div or -1 if anything went wrong.
1757  */
1758 static int add_div(struct isl_tab *tab, struct isl_context *context,
1759         struct isl_vec *div)
1760 {
1761         int r;
1762         int nonneg;
1763
1764         if ((nonneg = context->op->add_div(context, div)) < 0)
1765                 goto error;
1766
1767         if (!context->op->is_ok(context))
1768                 goto error;
1769
1770         if (isl_tab_extend_vars(tab, 1) < 0)
1771                 goto error;
1772         r = isl_tab_allocate_var(tab);
1773         if (r < 0)
1774                 goto error;
1775         if (nonneg)
1776                 tab->var[r].is_nonneg = 1;
1777         tab->var[r].frozen = 1;
1778         tab->n_div++;
1779
1780         return tab->n_div - 1;
1781 error:
1782         context->op->invalidate(context);
1783         return -1;
1784 }
1785
1786 static int find_div(struct isl_tab *tab, isl_int *div, isl_int denom)
1787 {
1788         int i;
1789         unsigned total = isl_basic_map_total_dim(tab->bmap);
1790
1791         for (i = 0; i < tab->bmap->n_div; ++i) {
1792                 if (isl_int_ne(tab->bmap->div[i][0], denom))
1793                         continue;
1794                 if (!isl_seq_eq(tab->bmap->div[i] + 1, div, 1 + total))
1795                         continue;
1796                 return i;
1797         }
1798         return -1;
1799 }
1800
1801 /* Return the index of a div that corresponds to "div".
1802  * We first check if we already have such a div and if not, we create one.
1803  */
1804 static int get_div(struct isl_tab *tab, struct isl_context *context,
1805         struct isl_vec *div)
1806 {
1807         int d;
1808         struct isl_tab *context_tab = context->op->peek_tab(context);
1809
1810         if (!context_tab)
1811                 return -1;
1812
1813         d = find_div(context_tab, div->el + 1, div->el[0]);
1814         if (d != -1)
1815                 return d;
1816
1817         return add_div(tab, context, div);
1818 }
1819
1820 /* Add a parametric cut to cut away the non-integral sample value
1821  * of the give row.
1822  * Let a_i be the coefficients of the constant term and the parameters
1823  * and let b_i be the coefficients of the variables or constraints
1824  * in basis of the tableau.
1825  * Let q be the div q = floor(\sum_i {-a_i} y_i).
1826  *
1827  * The cut is expressed as
1828  *
1829  *      c = \sum_i -{-a_i} y_i + \sum_i {b_i} x_i + q >= 0
1830  *
1831  * If q did not already exist in the context tableau, then it is added first.
1832  * If q is in a column of the main tableau then the "+ q" can be accomplished
1833  * by setting the corresponding entry to the denominator of the constraint.
1834  * If q happens to be in a row of the main tableau, then the corresponding
1835  * row needs to be added instead (taking care of the denominators).
1836  * Note that this is very unlikely, but perhaps not entirely impossible.
1837  *
1838  * The current value of the cut is known to be negative (or at least
1839  * non-positive), so row_sign is set accordingly.
1840  *
1841  * Return the row of the cut or -1.
1842  */
1843 static int add_parametric_cut(struct isl_tab *tab, int row,
1844         struct isl_context *context)
1845 {
1846         struct isl_vec *div;
1847         int d;
1848         int i;
1849         int r;
1850         isl_int *r_row;
1851         int col;
1852         int n;
1853         unsigned off = 2 + tab->M;
1854
1855         if (!context)
1856                 return -1;
1857
1858         div = get_row_parameter_div(tab, row);
1859         if (!div)
1860                 return -1;
1861
1862         n = tab->n_div;
1863         d = context->op->get_div(context, tab, div);
1864         if (d < 0)
1865                 return -1;
1866
1867         if (isl_tab_extend_cons(tab, 1) < 0)
1868                 return -1;
1869         r = isl_tab_allocate_con(tab);
1870         if (r < 0)
1871                 return -1;
1872
1873         r_row = tab->mat->row[tab->con[r].index];
1874         isl_int_set(r_row[0], tab->mat->row[row][0]);
1875         isl_int_neg(r_row[1], tab->mat->row[row][1]);
1876         isl_int_fdiv_r(r_row[1], r_row[1], tab->mat->row[row][0]);
1877         isl_int_neg(r_row[1], r_row[1]);
1878         if (tab->M)
1879                 isl_int_set_si(r_row[2], 0);
1880         for (i = 0; i < tab->n_param; ++i) {
1881                 if (tab->var[i].is_row)
1882                         continue;
1883                 col = tab->var[i].index;
1884                 isl_int_neg(r_row[off + col], tab->mat->row[row][off + col]);
1885                 isl_int_fdiv_r(r_row[off + col], r_row[off + col],
1886                                 tab->mat->row[row][0]);
1887                 isl_int_neg(r_row[off + col], r_row[off + col]);
1888         }
1889         for (i = 0; i < tab->n_div; ++i) {
1890                 if (tab->var[tab->n_var - tab->n_div + i].is_row)
1891                         continue;
1892                 col = tab->var[tab->n_var - tab->n_div + i].index;
1893                 isl_int_neg(r_row[off + col], tab->mat->row[row][off + col]);
1894                 isl_int_fdiv_r(r_row[off + col], r_row[off + col],
1895                                 tab->mat->row[row][0]);
1896                 isl_int_neg(r_row[off + col], r_row[off + col]);
1897         }
1898         for (i = 0; i < tab->n_col; ++i) {
1899                 if (tab->col_var[i] >= 0 &&
1900                     (tab->col_var[i] < tab->n_param ||
1901                      tab->col_var[i] >= tab->n_var - tab->n_div))
1902                         continue;
1903                 isl_int_fdiv_r(r_row[off + i],
1904                         tab->mat->row[row][off + i], tab->mat->row[row][0]);
1905         }
1906         if (tab->var[tab->n_var - tab->n_div + d].is_row) {
1907                 isl_int gcd;
1908                 int d_row = tab->var[tab->n_var - tab->n_div + d].index;
1909                 isl_int_init(gcd);
1910                 isl_int_gcd(gcd, tab->mat->row[d_row][0], r_row[0]);
1911                 isl_int_divexact(r_row[0], r_row[0], gcd);
1912                 isl_int_divexact(gcd, tab->mat->row[d_row][0], gcd);
1913                 isl_seq_combine(r_row + 1, gcd, r_row + 1,
1914                                 r_row[0], tab->mat->row[d_row] + 1,
1915                                 off - 1 + tab->n_col);
1916                 isl_int_mul(r_row[0], r_row[0], tab->mat->row[d_row][0]);
1917                 isl_int_clear(gcd);
1918         } else {
1919                 col = tab->var[tab->n_var - tab->n_div + d].index;
1920                 isl_int_set(r_row[off + col], tab->mat->row[row][0]);
1921         }
1922
1923         tab->con[r].is_nonneg = 1;
1924         if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1925                 return -1;
1926         if (tab->row_sign)
1927                 tab->row_sign[tab->con[r].index] = isl_tab_row_neg;
1928
1929         isl_vec_free(div);
1930
1931         row = tab->con[r].index;
1932
1933         if (d >= n && context->op->detect_equalities(context, tab) < 0)
1934                 return -1;
1935
1936         return row;
1937 }
1938
1939 /* Construct a tableau for bmap that can be used for computing
1940  * the lexicographic minimum (or maximum) of bmap.
1941  * If not NULL, then dom is the domain where the minimum
1942  * should be computed.  In this case, we set up a parametric
1943  * tableau with row signs (initialized to "unknown").
1944  * If M is set, then the tableau will use a big parameter.
1945  * If max is set, then a maximum should be computed instead of a minimum.
1946  * This means that for each variable x, the tableau will contain the variable
1947  * x' = M - x, rather than x' = M + x.  This in turn means that the coefficient
1948  * of the variables in all constraints are negated prior to adding them
1949  * to the tableau.
1950  */
1951 static struct isl_tab *tab_for_lexmin(struct isl_basic_map *bmap,
1952         struct isl_basic_set *dom, unsigned M, int max)
1953 {
1954         int i;
1955         struct isl_tab *tab;
1956
1957         tab = isl_tab_alloc(bmap->ctx, 2 * bmap->n_eq + bmap->n_ineq + 1,
1958                             isl_basic_map_total_dim(bmap), M);
1959         if (!tab)
1960                 return NULL;
1961
1962         tab->rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
1963         if (dom) {
1964                 tab->n_param = isl_basic_set_total_dim(dom) - dom->n_div;
1965                 tab->n_div = dom->n_div;
1966                 tab->row_sign = isl_calloc_array(bmap->ctx,
1967                                         enum isl_tab_row_sign, tab->mat->n_row);
1968                 if (!tab->row_sign)
1969                         goto error;
1970         }
1971         if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY)) {
1972                 if (isl_tab_mark_empty(tab) < 0)
1973                         goto error;
1974                 return tab;
1975         }
1976
1977         for (i = tab->n_param; i < tab->n_var - tab->n_div; ++i) {
1978                 tab->var[i].is_nonneg = 1;
1979                 tab->var[i].frozen = 1;
1980         }
1981         for (i = 0; i < bmap->n_eq; ++i) {
1982                 if (max)
1983                         isl_seq_neg(bmap->eq[i] + 1 + tab->n_param,
1984                                     bmap->eq[i] + 1 + tab->n_param,
1985                                     tab->n_var - tab->n_param - tab->n_div);
1986                 tab = add_lexmin_valid_eq(tab, bmap->eq[i]);
1987                 if (max)
1988                         isl_seq_neg(bmap->eq[i] + 1 + tab->n_param,
1989                                     bmap->eq[i] + 1 + tab->n_param,
1990                                     tab->n_var - tab->n_param - tab->n_div);
1991                 if (!tab || tab->empty)
1992                         return tab;
1993         }
1994         for (i = 0; i < bmap->n_ineq; ++i) {
1995                 if (max)
1996                         isl_seq_neg(bmap->ineq[i] + 1 + tab->n_param,
1997                                     bmap->ineq[i] + 1 + tab->n_param,
1998                                     tab->n_var - tab->n_param - tab->n_div);
1999                 tab = add_lexmin_ineq(tab, bmap->ineq[i]);
2000                 if (max)
2001                         isl_seq_neg(bmap->ineq[i] + 1 + tab->n_param,
2002                                     bmap->ineq[i] + 1 + tab->n_param,
2003                                     tab->n_var - tab->n_param - tab->n_div);
2004                 if (!tab || tab->empty)
2005                         return tab;
2006         }
2007         return tab;
2008 error:
2009         isl_tab_free(tab);
2010         return NULL;
2011 }
2012
2013 /* Given a main tableau where more than one row requires a split,
2014  * determine and return the "best" row to split on.
2015  *
2016  * Given two rows in the main tableau, if the inequality corresponding
2017  * to the first row is redundant with respect to that of the second row
2018  * in the current tableau, then it is better to split on the second row,
2019  * since in the positive part, both row will be positive.
2020  * (In the negative part a pivot will have to be performed and just about
2021  * anything can happen to the sign of the other row.)
2022  *
2023  * As a simple heuristic, we therefore select the row that makes the most
2024  * of the other rows redundant.
2025  *
2026  * Perhaps it would also be useful to look at the number of constraints
2027  * that conflict with any given constraint.
2028  */
2029 static int best_split(struct isl_tab *tab, struct isl_tab *context_tab)
2030 {
2031         struct isl_tab_undo *snap;
2032         int split;
2033         int row;
2034         int best = -1;
2035         int best_r;
2036
2037         if (isl_tab_extend_cons(context_tab, 2) < 0)
2038                 return -1;
2039
2040         snap = isl_tab_snap(context_tab);
2041
2042         for (split = tab->n_redundant; split < tab->n_row; ++split) {
2043                 struct isl_tab_undo *snap2;
2044                 struct isl_vec *ineq = NULL;
2045                 int r = 0;
2046                 int ok;
2047
2048                 if (!isl_tab_var_from_row(tab, split)->is_nonneg)
2049                         continue;
2050                 if (tab->row_sign[split] != isl_tab_row_any)
2051                         continue;
2052
2053                 ineq = get_row_parameter_ineq(tab, split);
2054                 if (!ineq)
2055                         return -1;
2056                 ok = isl_tab_add_ineq(context_tab, ineq->el) >= 0;
2057                 isl_vec_free(ineq);
2058                 if (!ok)
2059                         return -1;
2060
2061                 snap2 = isl_tab_snap(context_tab);
2062
2063                 for (row = tab->n_redundant; row < tab->n_row; ++row) {
2064                         struct isl_tab_var *var;
2065
2066                         if (row == split)
2067                                 continue;
2068                         if (!isl_tab_var_from_row(tab, row)->is_nonneg)
2069                                 continue;
2070                         if (tab->row_sign[row] != isl_tab_row_any)
2071                                 continue;
2072
2073                         ineq = get_row_parameter_ineq(tab, row);
2074                         if (!ineq)
2075                                 return -1;
2076                         ok = isl_tab_add_ineq(context_tab, ineq->el) >= 0;
2077                         isl_vec_free(ineq);
2078                         if (!ok)
2079                                 return -1;
2080                         var = &context_tab->con[context_tab->n_con - 1];
2081                         if (!context_tab->empty &&
2082                             !isl_tab_min_at_most_neg_one(context_tab, var))
2083                                 r++;
2084                         if (isl_tab_rollback(context_tab, snap2) < 0)
2085                                 return -1;
2086                 }
2087                 if (best == -1 || r > best_r) {
2088                         best = split;
2089                         best_r = r;
2090                 }
2091                 if (isl_tab_rollback(context_tab, snap) < 0)
2092                         return -1;
2093         }
2094
2095         return best;
2096 }
2097
2098 static struct isl_basic_set *context_lex_peek_basic_set(
2099         struct isl_context *context)
2100 {
2101         struct isl_context_lex *clex = (struct isl_context_lex *)context;
2102         if (!clex->tab)
2103                 return NULL;
2104         return isl_tab_peek_bset(clex->tab);
2105 }
2106
2107 static struct isl_tab *context_lex_peek_tab(struct isl_context *context)
2108 {
2109         struct isl_context_lex *clex = (struct isl_context_lex *)context;
2110         return clex->tab;
2111 }
2112
2113 static void context_lex_extend(struct isl_context *context, int n)
2114 {
2115         struct isl_context_lex *clex = (struct isl_context_lex *)context;
2116         if (!clex->tab)
2117                 return;
2118         if (isl_tab_extend_cons(clex->tab, n) >= 0)
2119                 return;
2120         isl_tab_free(clex->tab);
2121         clex->tab = NULL;
2122 }
2123
2124 static void context_lex_add_eq(struct isl_context *context, isl_int *eq,
2125                 int check, int update)
2126 {
2127         struct isl_context_lex *clex = (struct isl_context_lex *)context;
2128         if (isl_tab_extend_cons(clex->tab, 2) < 0)
2129                 goto error;
2130         clex->tab = add_lexmin_eq(clex->tab, eq);
2131         if (check) {
2132                 int v = tab_has_valid_sample(clex->tab, eq, 1);
2133                 if (v < 0)
2134                         goto error;
2135                 if (!v)
2136                         clex->tab = check_integer_feasible(clex->tab);
2137         }
2138         if (update)
2139                 clex->tab = check_samples(clex->tab, eq, 1);
2140         return;
2141 error:
2142         isl_tab_free(clex->tab);
2143         clex->tab = NULL;
2144 }
2145
2146 static void context_lex_add_ineq(struct isl_context *context, isl_int *ineq,
2147                 int check, int update)
2148 {
2149         struct isl_context_lex *clex = (struct isl_context_lex *)context;
2150         if (isl_tab_extend_cons(clex->tab, 1) < 0)
2151                 goto error;
2152         clex->tab = add_lexmin_ineq(clex->tab, ineq);
2153         if (check) {
2154                 int v = tab_has_valid_sample(clex->tab, ineq, 0);
2155                 if (v < 0)
2156                         goto error;
2157                 if (!v)
2158                         clex->tab = check_integer_feasible(clex->tab);
2159         }
2160         if (update)
2161                 clex->tab = check_samples(clex->tab, ineq, 0);
2162         return;
2163 error:
2164         isl_tab_free(clex->tab);
2165         clex->tab = NULL;
2166 }
2167
2168 static int context_lex_add_ineq_wrap(void *user, isl_int *ineq)
2169 {
2170         struct isl_context *context = (struct isl_context *)user;
2171         context_lex_add_ineq(context, ineq, 0, 0);
2172         return context->op->is_ok(context) ? 0 : -1;
2173 }
2174
2175 /* Check which signs can be obtained by "ineq" on all the currently
2176  * active sample values.  See row_sign for more information.
2177  */
2178 static enum isl_tab_row_sign tab_ineq_sign(struct isl_tab *tab, isl_int *ineq,
2179         int strict)
2180 {
2181         int i;
2182         int sgn;
2183         isl_int tmp;
2184         enum isl_tab_row_sign res = isl_tab_row_unknown;
2185
2186         isl_assert(tab->mat->ctx, tab->samples, return isl_tab_row_unknown);
2187         isl_assert(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var,
2188                         return isl_tab_row_unknown);
2189
2190         isl_int_init(tmp);
2191         for (i = tab->n_outside; i < tab->n_sample; ++i) {
2192                 isl_seq_inner_product(tab->samples->row[i], ineq,
2193                                         1 + tab->n_var, &tmp);
2194                 sgn = isl_int_sgn(tmp);
2195                 if (sgn > 0 || (sgn == 0 && strict)) {
2196                         if (res == isl_tab_row_unknown)
2197                                 res = isl_tab_row_pos;
2198                         if (res == isl_tab_row_neg)
2199                                 res = isl_tab_row_any;
2200                 }
2201                 if (sgn < 0) {
2202                         if (res == isl_tab_row_unknown)
2203                                 res = isl_tab_row_neg;
2204                         if (res == isl_tab_row_pos)
2205                                 res = isl_tab_row_any;
2206                 }
2207                 if (res == isl_tab_row_any)
2208                         break;
2209         }
2210         isl_int_clear(tmp);
2211
2212         return res;
2213 }
2214
2215 static enum isl_tab_row_sign context_lex_ineq_sign(struct isl_context *context,
2216                         isl_int *ineq, int strict)
2217 {
2218         struct isl_context_lex *clex = (struct isl_context_lex *)context;
2219         return tab_ineq_sign(clex->tab, ineq, strict);
2220 }
2221
2222 /* Check whether "ineq" can be added to the tableau without rendering
2223  * it infeasible.
2224  */
2225 static int context_lex_test_ineq(struct isl_context *context, isl_int *ineq)
2226 {
2227         struct isl_context_lex *clex = (struct isl_context_lex *)context;
2228         struct isl_tab_undo *snap;
2229         int feasible;
2230
2231         if (!clex->tab)
2232                 return -1;
2233
2234         if (isl_tab_extend_cons(clex->tab, 1) < 0)
2235                 return -1;
2236
2237         snap = isl_tab_snap(clex->tab);
2238         if (isl_tab_push_basis(clex->tab) < 0)
2239                 return -1;
2240         clex->tab = add_lexmin_ineq(clex->tab, ineq);
2241         clex->tab = check_integer_feasible(clex->tab);
2242         if (!clex->tab)
2243                 return -1;
2244         feasible = !clex->tab->empty;
2245         if (isl_tab_rollback(clex->tab, snap) < 0)
2246                 return -1;
2247
2248         return feasible;
2249 }
2250
2251 static int context_lex_get_div(struct isl_context *context, struct isl_tab *tab,
2252                 struct isl_vec *div)
2253 {
2254         return get_div(tab, context, div);
2255 }
2256
2257 static int context_lex_add_div(struct isl_context *context, struct isl_vec *div)
2258 {
2259         struct isl_context_lex *clex = (struct isl_context_lex *)context;
2260         return context_tab_add_div(clex->tab, div,
2261                                         context_lex_add_ineq_wrap, context);
2262 }
2263
2264 static int context_lex_detect_equalities(struct isl_context *context,
2265                 struct isl_tab *tab)
2266 {
2267         return 0;
2268 }
2269
2270 static int context_lex_best_split(struct isl_context *context,
2271                 struct isl_tab *tab)
2272 {
2273         struct isl_context_lex *clex = (struct isl_context_lex *)context;
2274         struct isl_tab_undo *snap;
2275         int r;
2276
2277         snap = isl_tab_snap(clex->tab);
2278         if (isl_tab_push_basis(clex->tab) < 0)
2279                 return -1;
2280         r = best_split(tab, clex->tab);
2281
2282         if (r >= 0 && isl_tab_rollback(clex->tab, snap) < 0)
2283                 return -1;
2284
2285         return r;
2286 }
2287
2288 static int context_lex_is_empty(struct isl_context *context)
2289 {
2290         struct isl_context_lex *clex = (struct isl_context_lex *)context;
2291         if (!clex->tab)
2292                 return -1;
2293         return clex->tab->empty;
2294 }
2295
2296 static void *context_lex_save(struct isl_context *context)
2297 {
2298         struct isl_context_lex *clex = (struct isl_context_lex *)context;
2299         struct isl_tab_undo *snap;
2300
2301         snap = isl_tab_snap(clex->tab);
2302         if (isl_tab_push_basis(clex->tab) < 0)
2303                 return NULL;
2304         if (isl_tab_save_samples(clex->tab) < 0)
2305                 return NULL;
2306
2307         return snap;
2308 }
2309
2310 static void context_lex_restore(struct isl_context *context, void *save)
2311 {
2312         struct isl_context_lex *clex = (struct isl_context_lex *)context;
2313         if (isl_tab_rollback(clex->tab, (struct isl_tab_undo *)save) < 0) {
2314                 isl_tab_free(clex->tab);
2315                 clex->tab = NULL;
2316         }
2317 }
2318
2319 static int context_lex_is_ok(struct isl_context *context)
2320 {
2321         struct isl_context_lex *clex = (struct isl_context_lex *)context;
2322         return !!clex->tab;
2323 }
2324
2325 /* For each variable in the context tableau, check if the variable can
2326  * only attain non-negative values.  If so, mark the parameter as non-negative
2327  * in the main tableau.  This allows for a more direct identification of some
2328  * cases of violated constraints.
2329  */
2330 static struct isl_tab *tab_detect_nonnegative_parameters(struct isl_tab *tab,
2331         struct isl_tab *context_tab)
2332 {
2333         int i;
2334         struct isl_tab_undo *snap;
2335         struct isl_vec *ineq = NULL;
2336         struct isl_tab_var *var;
2337         int n;
2338
2339         if (context_tab->n_var == 0)
2340                 return tab;
2341
2342         ineq = isl_vec_alloc(tab->mat->ctx, 1 + context_tab->n_var);
2343         if (!ineq)
2344                 goto error;
2345
2346         if (isl_tab_extend_cons(context_tab, 1) < 0)
2347                 goto error;
2348
2349         snap = isl_tab_snap(context_tab);
2350
2351         n = 0;
2352         isl_seq_clr(ineq->el, ineq->size);
2353         for (i = 0; i < context_tab->n_var; ++i) {
2354                 isl_int_set_si(ineq->el[1 + i], 1);
2355                 if (isl_tab_add_ineq(context_tab, ineq->el) < 0)
2356                         goto error;
2357                 var = &context_tab->con[context_tab->n_con - 1];
2358                 if (!context_tab->empty &&
2359                     !isl_tab_min_at_most_neg_one(context_tab, var)) {
2360                         int j = i;
2361                         if (i >= tab->n_param)
2362                                 j = i - tab->n_param + tab->n_var - tab->n_div;
2363                         tab->var[j].is_nonneg = 1;
2364                         n++;
2365                 }
2366                 isl_int_set_si(ineq->el[1 + i], 0);
2367                 if (isl_tab_rollback(context_tab, snap) < 0)
2368                         goto error;
2369         }
2370
2371         if (context_tab->M && n == context_tab->n_var) {
2372                 context_tab->mat = isl_mat_drop_cols(context_tab->mat, 2, 1);
2373                 context_tab->M = 0;
2374         }
2375
2376         isl_vec_free(ineq);
2377         return tab;
2378 error:
2379         isl_vec_free(ineq);
2380         isl_tab_free(tab);
2381         return NULL;
2382 }
2383
2384 static struct isl_tab *context_lex_detect_nonnegative_parameters(
2385         struct isl_context *context, struct isl_tab *tab)
2386 {
2387         struct isl_context_lex *clex = (struct isl_context_lex *)context;
2388         struct isl_tab_undo *snap;
2389
2390         if (!tab)
2391                 return NULL;
2392
2393         snap = isl_tab_snap(clex->tab);
2394         if (isl_tab_push_basis(clex->tab) < 0)
2395                 goto error;
2396
2397         tab = tab_detect_nonnegative_parameters(tab, clex->tab);
2398
2399         if (isl_tab_rollback(clex->tab, snap) < 0)
2400                 goto error;
2401
2402         return tab;
2403 error:
2404         isl_tab_free(tab);
2405         return NULL;
2406 }
2407
2408 static void context_lex_invalidate(struct isl_context *context)
2409 {
2410         struct isl_context_lex *clex = (struct isl_context_lex *)context;
2411         isl_tab_free(clex->tab);
2412         clex->tab = NULL;
2413 }
2414
2415 static void context_lex_free(struct isl_context *context)
2416 {
2417         struct isl_context_lex *clex = (struct isl_context_lex *)context;
2418         isl_tab_free(clex->tab);
2419         free(clex);
2420 }
2421
2422 struct isl_context_op isl_context_lex_op = {
2423         context_lex_detect_nonnegative_parameters,
2424         context_lex_peek_basic_set,
2425         context_lex_peek_tab,
2426         context_lex_add_eq,
2427         context_lex_add_ineq,
2428         context_lex_ineq_sign,
2429         context_lex_test_ineq,
2430         context_lex_get_div,
2431         context_lex_add_div,
2432         context_lex_detect_equalities,
2433         context_lex_best_split,
2434         context_lex_is_empty,
2435         context_lex_is_ok,
2436         context_lex_save,
2437         context_lex_restore,
2438         context_lex_invalidate,
2439         context_lex_free,
2440 };
2441
2442 static struct isl_tab *context_tab_for_lexmin(struct isl_basic_set *bset)
2443 {
2444         struct isl_tab *tab;
2445
2446         bset = isl_basic_set_cow(bset);
2447         if (!bset)
2448                 return NULL;
2449         tab = tab_for_lexmin((struct isl_basic_map *)bset, NULL, 1, 0);
2450         if (!tab)
2451                 goto error;
2452         if (isl_tab_track_bset(tab, bset) < 0)
2453                 goto error;
2454         tab = isl_tab_init_samples(tab);
2455         return tab;
2456 error:
2457         isl_basic_set_free(bset);
2458         return NULL;
2459 }
2460
2461 static struct isl_context *isl_context_lex_alloc(struct isl_basic_set *dom)
2462 {
2463         struct isl_context_lex *clex;
2464
2465         if (!dom)
2466                 return NULL;
2467
2468         clex = isl_alloc_type(dom->ctx, struct isl_context_lex);
2469         if (!clex)
2470                 return NULL;
2471
2472         clex->context.op = &isl_context_lex_op;
2473
2474         clex->tab = context_tab_for_lexmin(isl_basic_set_copy(dom));
2475         clex->tab = restore_lexmin(clex->tab);
2476         clex->tab = check_integer_feasible(clex->tab);
2477         if (!clex->tab)
2478                 goto error;
2479
2480         return &clex->context;
2481 error:
2482         clex->context.op->free(&clex->context);
2483         return NULL;
2484 }
2485
2486 struct isl_context_gbr {
2487         struct isl_context context;
2488         struct isl_tab *tab;
2489         struct isl_tab *shifted;
2490         struct isl_tab *cone;
2491 };
2492
2493 static struct isl_tab *context_gbr_detect_nonnegative_parameters(
2494         struct isl_context *context, struct isl_tab *tab)
2495 {
2496         struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2497         if (!tab)
2498                 return NULL;
2499         return tab_detect_nonnegative_parameters(tab, cgbr->tab);
2500 }
2501
2502 static struct isl_basic_set *context_gbr_peek_basic_set(
2503         struct isl_context *context)
2504 {
2505         struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2506         if (!cgbr->tab)
2507                 return NULL;
2508         return isl_tab_peek_bset(cgbr->tab);
2509 }
2510
2511 static struct isl_tab *context_gbr_peek_tab(struct isl_context *context)
2512 {
2513         struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2514         return cgbr->tab;
2515 }
2516
2517 /* Initialize the "shifted" tableau of the context, which
2518  * contains the constraints of the original tableau shifted
2519  * by the sum of all negative coefficients.  This ensures
2520  * that any rational point in the shifted tableau can
2521  * be rounded up to yield an integer point in the original tableau.
2522  */
2523 static void gbr_init_shifted(struct isl_context_gbr *cgbr)
2524 {
2525         int i, j;
2526         struct isl_vec *cst;
2527         struct isl_basic_set *bset = isl_tab_peek_bset(cgbr->tab);
2528         unsigned dim = isl_basic_set_total_dim(bset);
2529
2530         cst = isl_vec_alloc(cgbr->tab->mat->ctx, bset->n_ineq);
2531         if (!cst)
2532                 return;
2533
2534         for (i = 0; i < bset->n_ineq; ++i) {
2535                 isl_int_set(cst->el[i], bset->ineq[i][0]);
2536                 for (j = 0; j < dim; ++j) {
2537                         if (!isl_int_is_neg(bset->ineq[i][1 + j]))
2538                                 continue;
2539                         isl_int_add(bset->ineq[i][0], bset->ineq[i][0],
2540                                     bset->ineq[i][1 + j]);
2541                 }
2542         }
2543
2544         cgbr->shifted = isl_tab_from_basic_set(bset);
2545
2546         for (i = 0; i < bset->n_ineq; ++i)
2547                 isl_int_set(bset->ineq[i][0], cst->el[i]);
2548
2549         isl_vec_free(cst);
2550 }
2551
2552 /* Check if the shifted tableau is non-empty, and if so
2553  * use the sample point to construct an integer point
2554  * of the context tableau.
2555  */
2556 static struct isl_vec *gbr_get_shifted_sample(struct isl_context_gbr *cgbr)
2557 {
2558         struct isl_vec *sample;
2559
2560         if (!cgbr->shifted)
2561                 gbr_init_shifted(cgbr);
2562         if (!cgbr->shifted)
2563                 return NULL;
2564         if (cgbr->shifted->empty)
2565                 return isl_vec_alloc(cgbr->tab->mat->ctx, 0);
2566
2567         sample = isl_tab_get_sample_value(cgbr->shifted);
2568         sample = isl_vec_ceil(sample);
2569
2570         return sample;
2571 }
2572
2573 static struct isl_basic_set *drop_constant_terms(struct isl_basic_set *bset)
2574 {
2575         int i;
2576
2577         if (!bset)
2578                 return NULL;
2579
2580         for (i = 0; i < bset->n_eq; ++i)
2581                 isl_int_set_si(bset->eq[i][0], 0);
2582
2583         for (i = 0; i < bset->n_ineq; ++i)
2584                 isl_int_set_si(bset->ineq[i][0], 0);
2585
2586         return bset;
2587 }
2588
2589 static int use_shifted(struct isl_context_gbr *cgbr)
2590 {
2591         return cgbr->tab->bmap->n_eq == 0 && cgbr->tab->bmap->n_div == 0;
2592 }
2593
2594 static struct isl_vec *gbr_get_sample(struct isl_context_gbr *cgbr)
2595 {
2596         struct isl_basic_set *bset;
2597         struct isl_basic_set *cone;
2598
2599         if (isl_tab_sample_is_integer(cgbr->tab))
2600                 return isl_tab_get_sample_value(cgbr->tab);
2601
2602         if (use_shifted(cgbr)) {
2603                 struct isl_vec *sample;
2604
2605                 sample = gbr_get_shifted_sample(cgbr);
2606                 if (!sample || sample->size > 0)
2607                         return sample;
2608
2609                 isl_vec_free(sample);
2610         }
2611
2612         if (!cgbr->cone) {
2613                 bset = isl_tab_peek_bset(cgbr->tab);
2614                 cgbr->cone = isl_tab_from_recession_cone(bset, 0);
2615                 if (!cgbr->cone)
2616                         return NULL;
2617                 if (isl_tab_track_bset(cgbr->cone, isl_basic_set_dup(bset)) < 0)
2618                         return NULL;
2619         }
2620         if (isl_tab_detect_implicit_equalities(cgbr->cone) < 0)
2621                 return NULL;
2622
2623         if (cgbr->cone->n_dead == cgbr->cone->n_col) {
2624                 struct isl_vec *sample;
2625                 struct isl_tab_undo *snap;
2626
2627                 if (cgbr->tab->basis) {
2628                         if (cgbr->tab->basis->n_col != 1 + cgbr->tab->n_var) {
2629                                 isl_mat_free(cgbr->tab->basis);
2630                                 cgbr->tab->basis = NULL;
2631                         }
2632                         cgbr->tab->n_zero = 0;
2633                         cgbr->tab->n_unbounded = 0;
2634                 }
2635
2636                 snap = isl_tab_snap(cgbr->tab);
2637
2638                 sample = isl_tab_sample(cgbr->tab);
2639
2640                 if (isl_tab_rollback(cgbr->tab, snap) < 0) {
2641                         isl_vec_free(sample);
2642                         return NULL;
2643                 }
2644
2645                 return sample;
2646         }
2647
2648         cone = isl_basic_set_dup(isl_tab_peek_bset(cgbr->cone));
2649         cone = drop_constant_terms(cone);
2650         cone = isl_basic_set_update_from_tab(cone, cgbr->cone);
2651         cone = isl_basic_set_underlying_set(cone);
2652         cone = isl_basic_set_gauss(cone, NULL);
2653
2654         bset = isl_basic_set_dup(isl_tab_peek_bset(cgbr->tab));
2655         bset = isl_basic_set_update_from_tab(bset, cgbr->tab);
2656         bset = isl_basic_set_underlying_set(bset);
2657         bset = isl_basic_set_gauss(bset, NULL);
2658
2659         return isl_basic_set_sample_with_cone(bset, cone);
2660 }
2661
2662 static void check_gbr_integer_feasible(struct isl_context_gbr *cgbr)
2663 {
2664         struct isl_vec *sample;
2665
2666         if (!cgbr->tab)
2667                 return;
2668
2669         if (cgbr->tab->empty)
2670                 return;
2671
2672         sample = gbr_get_sample(cgbr);
2673         if (!sample)
2674                 goto error;
2675
2676         if (sample->size == 0) {
2677                 isl_vec_free(sample);
2678                 if (isl_tab_mark_empty(cgbr->tab) < 0)
2679                         goto error;
2680                 return;
2681         }
2682
2683         cgbr->tab = isl_tab_add_sample(cgbr->tab, sample);
2684
2685         return;
2686 error:
2687         isl_tab_free(cgbr->tab);
2688         cgbr->tab = NULL;
2689 }
2690
2691 static struct isl_tab *add_gbr_eq(struct isl_tab *tab, isl_int *eq)
2692 {
2693         int r;
2694
2695         if (!tab)
2696                 return NULL;
2697
2698         if (isl_tab_extend_cons(tab, 2) < 0)
2699                 goto error;
2700
2701         if (isl_tab_add_eq(tab, eq) < 0)
2702                 goto error;
2703
2704         return tab;
2705 error:
2706         isl_tab_free(tab);
2707         return NULL;
2708 }
2709
2710 static void context_gbr_add_eq(struct isl_context *context, isl_int *eq,
2711                 int check, int update)
2712 {
2713         struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2714
2715         cgbr->tab = add_gbr_eq(cgbr->tab, eq);
2716
2717         if (cgbr->cone && cgbr->cone->n_col != cgbr->cone->n_dead) {
2718                 if (isl_tab_extend_cons(cgbr->cone, 2) < 0)
2719                         goto error;
2720                 if (isl_tab_add_eq(cgbr->cone, eq) < 0)
2721                         goto error;
2722         }
2723
2724         if (check) {
2725                 int v = tab_has_valid_sample(cgbr->tab, eq, 1);
2726                 if (v < 0)
2727                         goto error;
2728                 if (!v)
2729                         check_gbr_integer_feasible(cgbr);
2730         }
2731         if (update)
2732                 cgbr->tab = check_samples(cgbr->tab, eq, 1);
2733         return;
2734 error:
2735         isl_tab_free(cgbr->tab);
2736         cgbr->tab = NULL;
2737 }
2738
2739 static void add_gbr_ineq(struct isl_context_gbr *cgbr, isl_int *ineq)
2740 {
2741         if (!cgbr->tab)
2742                 return;
2743
2744         if (isl_tab_extend_cons(cgbr->tab, 1) < 0)
2745                 goto error;
2746
2747         if (isl_tab_add_ineq(cgbr->tab, ineq) < 0)
2748                 goto error;
2749
2750         if (cgbr->shifted && !cgbr->shifted->empty && use_shifted(cgbr)) {
2751                 int i;
2752                 unsigned dim;
2753                 dim = isl_basic_map_total_dim(cgbr->tab->bmap);
2754
2755                 if (isl_tab_extend_cons(cgbr->shifted, 1) < 0)
2756                         goto error;
2757
2758                 for (i = 0; i < dim; ++i) {
2759                         if (!isl_int_is_neg(ineq[1 + i]))
2760                                 continue;
2761                         isl_int_add(ineq[0], ineq[0], ineq[1 + i]);
2762                 }
2763
2764                 if (isl_tab_add_ineq(cgbr->shifted, ineq) < 0)
2765                         goto error;
2766
2767                 for (i = 0; i < dim; ++i) {
2768                         if (!isl_int_is_neg(ineq[1 + i]))
2769                                 continue;
2770                         isl_int_sub(ineq[0], ineq[0], ineq[1 + i]);
2771                 }
2772         }
2773
2774         if (cgbr->cone && cgbr->cone->n_col != cgbr->cone->n_dead) {
2775                 if (isl_tab_extend_cons(cgbr->cone, 1) < 0)
2776                         goto error;
2777                 if (isl_tab_add_ineq(cgbr->cone, ineq) < 0)
2778                         goto error;
2779         }
2780
2781         return;
2782 error:
2783         isl_tab_free(cgbr->tab);
2784         cgbr->tab = NULL;
2785 }
2786
2787 static void context_gbr_add_ineq(struct isl_context *context, isl_int *ineq,
2788                 int check, int update)
2789 {
2790         struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2791
2792         add_gbr_ineq(cgbr, ineq);
2793         if (!cgbr->tab)
2794                 return;
2795
2796         if (check) {
2797                 int v = tab_has_valid_sample(cgbr->tab, ineq, 0);
2798                 if (v < 0)
2799                         goto error;
2800                 if (!v)
2801                         check_gbr_integer_feasible(cgbr);
2802         }
2803         if (update)
2804                 cgbr->tab = check_samples(cgbr->tab, ineq, 0);
2805         return;
2806 error:
2807         isl_tab_free(cgbr->tab);
2808         cgbr->tab = NULL;
2809 }
2810
2811 static int context_gbr_add_ineq_wrap(void *user, isl_int *ineq)
2812 {
2813         struct isl_context *context = (struct isl_context *)user;
2814         context_gbr_add_ineq(context, ineq, 0, 0);
2815         return context->op->is_ok(context) ? 0 : -1;
2816 }
2817
2818 static enum isl_tab_row_sign context_gbr_ineq_sign(struct isl_context *context,
2819                         isl_int *ineq, int strict)
2820 {
2821         struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2822         return tab_ineq_sign(cgbr->tab, ineq, strict);
2823 }
2824
2825 /* Check whether "ineq" can be added to the tableau without rendering
2826  * it infeasible.
2827  */
2828 static int context_gbr_test_ineq(struct isl_context *context, isl_int *ineq)
2829 {
2830         struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2831         struct isl_tab_undo *snap;
2832         struct isl_tab_undo *shifted_snap = NULL;
2833         struct isl_tab_undo *cone_snap = NULL;
2834         int feasible;
2835
2836         if (!cgbr->tab)
2837                 return -1;
2838
2839         if (isl_tab_extend_cons(cgbr->tab, 1) < 0)
2840                 return -1;
2841
2842         snap = isl_tab_snap(cgbr->tab);
2843         if (cgbr->shifted)
2844                 shifted_snap = isl_tab_snap(cgbr->shifted);
2845         if (cgbr->cone)
2846                 cone_snap = isl_tab_snap(cgbr->cone);
2847         add_gbr_ineq(cgbr, ineq);
2848         check_gbr_integer_feasible(cgbr);
2849         if (!cgbr->tab)
2850                 return -1;
2851         feasible = !cgbr->tab->empty;
2852         if (isl_tab_rollback(cgbr->tab, snap) < 0)
2853                 return -1;
2854         if (shifted_snap) {
2855                 if (isl_tab_rollback(cgbr->shifted, shifted_snap))
2856                         return -1;
2857         } else if (cgbr->shifted) {
2858                 isl_tab_free(cgbr->shifted);
2859                 cgbr->shifted = NULL;
2860         }
2861         if (cone_snap) {
2862                 if (isl_tab_rollback(cgbr->cone, cone_snap))
2863                         return -1;
2864         } else if (cgbr->cone) {
2865                 isl_tab_free(cgbr->cone);
2866                 cgbr->cone = NULL;
2867         }
2868
2869         return feasible;
2870 }
2871
2872 /* Return the column of the last of the variables associated to
2873  * a column that has a non-zero coefficient.
2874  * This function is called in a context where only coefficients
2875  * of parameters or divs can be non-zero.
2876  */
2877 static int last_non_zero_var_col(struct isl_tab *tab, isl_int *p)
2878 {
2879         int i;
2880         int col;
2881         unsigned dim = tab->n_var - tab->n_param - tab->n_div;
2882
2883         if (tab->n_var == 0)
2884                 return -1;
2885
2886         for (i = tab->n_var - 1; i >= 0; --i) {
2887                 if (i >= tab->n_param && i < tab->n_var - tab->n_div)
2888                         continue;
2889                 if (tab->var[i].is_row)
2890                         continue;
2891                 col = tab->var[i].index;
2892                 if (!isl_int_is_zero(p[col]))
2893                         return col;
2894         }
2895
2896         return -1;
2897 }
2898
2899 /* Look through all the recently added equalities in the context
2900  * to see if we can propagate any of them to the main tableau.
2901  *
2902  * The newly added equalities in the context are encoded as pairs
2903  * of inequalities starting at inequality "first".
2904  *
2905  * We tentatively add each of these equalities to the main tableau
2906  * and if this happens to result in a row with a final coefficient
2907  * that is one or negative one, we use it to kill a column
2908  * in the main tableau.  Otherwise, we discard the tentatively
2909  * added row.
2910  */
2911 static void propagate_equalities(struct isl_context_gbr *cgbr,
2912         struct isl_tab *tab, unsigned first)
2913 {
2914         int i;
2915         struct isl_vec *eq = NULL;
2916
2917         eq = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
2918         if (!eq)
2919                 goto error;
2920
2921         if (isl_tab_extend_cons(tab, (cgbr->tab->bmap->n_ineq - first)/2) < 0)
2922                 goto error;
2923
2924         isl_seq_clr(eq->el + 1 + tab->n_param,
2925                     tab->n_var - tab->n_param - tab->n_div);
2926         for (i = first; i < cgbr->tab->bmap->n_ineq; i += 2) {
2927                 int j;
2928                 int r;
2929                 struct isl_tab_undo *snap;
2930                 snap = isl_tab_snap(tab);
2931
2932                 isl_seq_cpy(eq->el, cgbr->tab->bmap->ineq[i], 1 + tab->n_param);
2933                 isl_seq_cpy(eq->el + 1 + tab->n_var - tab->n_div,
2934                             cgbr->tab->bmap->ineq[i] + 1 + tab->n_param,
2935                             tab->n_div);
2936
2937                 r = isl_tab_add_row(tab, eq->el);
2938                 if (r < 0)
2939                         goto error;
2940                 r = tab->con[r].index;
2941                 j = last_non_zero_var_col(tab, tab->mat->row[r] + 2 + tab->M);
2942                 if (j < 0 || j < tab->n_dead ||
2943                     !isl_int_is_one(tab->mat->row[r][0]) ||
2944                     (!isl_int_is_one(tab->mat->row[r][2 + tab->M + j]) &&
2945                      !isl_int_is_negone(tab->mat->row[r][2 + tab->M + j]))) {
2946                         if (isl_tab_rollback(tab, snap) < 0)
2947                                 goto error;
2948                         continue;
2949                 }
2950                 if (isl_tab_pivot(tab, r, j) < 0)
2951                         goto error;
2952                 if (isl_tab_kill_col(tab, j) < 0)
2953                         goto error;
2954
2955                 tab = restore_lexmin(tab);
2956         }
2957
2958         isl_vec_free(eq);
2959
2960         return;
2961 error:
2962         isl_vec_free(eq);
2963         isl_tab_free(cgbr->tab);
2964         cgbr->tab = NULL;
2965 }
2966
2967 static int context_gbr_detect_equalities(struct isl_context *context,
2968         struct isl_tab *tab)
2969 {
2970         struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
2971         struct isl_ctx *ctx;
2972         int i;
2973         enum isl_lp_result res;
2974         unsigned n_ineq;
2975
2976         ctx = cgbr->tab->mat->ctx;
2977
2978         if (!cgbr->cone) {
2979                 struct isl_basic_set *bset = isl_tab_peek_bset(cgbr->tab);
2980                 cgbr->cone = isl_tab_from_recession_cone(bset, 0);
2981                 if (!cgbr->cone)
2982                         goto error;
2983                 if (isl_tab_track_bset(cgbr->cone, isl_basic_set_dup(bset)) < 0)
2984                         goto error;
2985         }
2986         if (isl_tab_detect_implicit_equalities(cgbr->cone) < 0)
2987                 goto error;
2988
2989         n_ineq = cgbr->tab->bmap->n_ineq;
2990         cgbr->tab = isl_tab_detect_equalities(cgbr->tab, cgbr->cone);
2991         if (cgbr->tab && cgbr->tab->bmap->n_ineq > n_ineq)
2992                 propagate_equalities(cgbr, tab, n_ineq);
2993
2994         return 0;
2995 error:
2996         isl_tab_free(cgbr->tab);
2997         cgbr->tab = NULL;
2998         return -1;
2999 }
3000
3001 static int context_gbr_get_div(struct isl_context *context, struct isl_tab *tab,
3002                 struct isl_vec *div)
3003 {
3004         return get_div(tab, context, div);
3005 }
3006
3007 static int context_gbr_add_div(struct isl_context *context, struct isl_vec *div)
3008 {
3009         struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3010         if (cgbr->cone) {
3011                 int k;
3012
3013                 if (isl_tab_extend_cons(cgbr->cone, 3) < 0)
3014                         return -1;
3015                 if (isl_tab_extend_vars(cgbr->cone, 1) < 0)
3016                         return -1;
3017                 if (isl_tab_allocate_var(cgbr->cone) <0)
3018                         return -1;
3019
3020                 cgbr->cone->bmap = isl_basic_map_extend_dim(cgbr->cone->bmap,
3021                         isl_basic_map_get_dim(cgbr->cone->bmap), 1, 0, 2);
3022                 k = isl_basic_map_alloc_div(cgbr->cone->bmap);
3023                 if (k < 0)
3024                         return -1;
3025                 isl_seq_cpy(cgbr->cone->bmap->div[k], div->el, div->size);
3026                 if (isl_tab_push(cgbr->cone, isl_tab_undo_bmap_div) < 0)
3027                         return -1;
3028         }
3029         return context_tab_add_div(cgbr->tab, div,
3030                                         context_gbr_add_ineq_wrap, context);
3031 }
3032
3033 static int context_gbr_best_split(struct isl_context *context,
3034                 struct isl_tab *tab)
3035 {
3036         struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3037         struct isl_tab_undo *snap;
3038         int r;
3039
3040         snap = isl_tab_snap(cgbr->tab);
3041         r = best_split(tab, cgbr->tab);
3042
3043         if (r >= 0 && isl_tab_rollback(cgbr->tab, snap) < 0)
3044                 return -1;
3045
3046         return r;
3047 }
3048
3049 static int context_gbr_is_empty(struct isl_context *context)
3050 {
3051         struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3052         if (!cgbr->tab)
3053                 return -1;
3054         return cgbr->tab->empty;
3055 }
3056
3057 struct isl_gbr_tab_undo {
3058         struct isl_tab_undo *tab_snap;
3059         struct isl_tab_undo *shifted_snap;
3060         struct isl_tab_undo *cone_snap;
3061 };
3062
3063 static void *context_gbr_save(struct isl_context *context)
3064 {
3065         struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3066         struct isl_gbr_tab_undo *snap;
3067
3068         snap = isl_alloc_type(cgbr->tab->mat->ctx, struct isl_gbr_tab_undo);
3069         if (!snap)
3070                 return NULL;
3071
3072         snap->tab_snap = isl_tab_snap(cgbr->tab);
3073         if (isl_tab_save_samples(cgbr->tab) < 0)
3074                 goto error;
3075
3076         if (cgbr->shifted)
3077                 snap->shifted_snap = isl_tab_snap(cgbr->shifted);
3078         else
3079                 snap->shifted_snap = NULL;
3080
3081         if (cgbr->cone)
3082                 snap->cone_snap = isl_tab_snap(cgbr->cone);
3083         else
3084                 snap->cone_snap = NULL;
3085
3086         return snap;
3087 error:
3088         free(snap);
3089         return NULL;
3090 }
3091
3092 static void context_gbr_restore(struct isl_context *context, void *save)
3093 {
3094         struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3095         struct isl_gbr_tab_undo *snap = (struct isl_gbr_tab_undo *)save;
3096         if (!snap)
3097                 goto error;
3098         if (isl_tab_rollback(cgbr->tab, snap->tab_snap) < 0) {
3099                 isl_tab_free(cgbr->tab);
3100                 cgbr->tab = NULL;
3101         }
3102
3103         if (snap->shifted_snap) {
3104                 if (isl_tab_rollback(cgbr->shifted, snap->shifted_snap) < 0)
3105                         goto error;
3106         } else if (cgbr->shifted) {
3107                 isl_tab_free(cgbr->shifted);
3108                 cgbr->shifted = NULL;
3109         }
3110
3111         if (snap->cone_snap) {
3112                 if (isl_tab_rollback(cgbr->cone, snap->cone_snap) < 0)
3113                         goto error;
3114         } else if (cgbr->cone) {
3115                 isl_tab_free(cgbr->cone);
3116                 cgbr->cone = NULL;
3117         }
3118
3119         free(snap);
3120
3121         return;
3122 error:
3123         free(snap);
3124         isl_tab_free(cgbr->tab);
3125         cgbr->tab = NULL;
3126 }
3127
3128 static int context_gbr_is_ok(struct isl_context *context)
3129 {
3130         struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3131         return !!cgbr->tab;
3132 }
3133
3134 static void context_gbr_invalidate(struct isl_context *context)
3135 {
3136         struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3137         isl_tab_free(cgbr->tab);
3138         cgbr->tab = NULL;
3139 }
3140
3141 static void context_gbr_free(struct isl_context *context)
3142 {
3143         struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3144         isl_tab_free(cgbr->tab);
3145         isl_tab_free(cgbr->shifted);
3146         isl_tab_free(cgbr->cone);
3147         free(cgbr);
3148 }
3149
3150 struct isl_context_op isl_context_gbr_op = {
3151         context_gbr_detect_nonnegative_parameters,
3152         context_gbr_peek_basic_set,
3153         context_gbr_peek_tab,
3154         context_gbr_add_eq,
3155         context_gbr_add_ineq,
3156         context_gbr_ineq_sign,
3157         context_gbr_test_ineq,
3158         context_gbr_get_div,
3159         context_gbr_add_div,
3160         context_gbr_detect_equalities,
3161         context_gbr_best_split,
3162         context_gbr_is_empty,
3163         context_gbr_is_ok,
3164         context_gbr_save,
3165         context_gbr_restore,
3166         context_gbr_invalidate,
3167         context_gbr_free,
3168 };
3169
3170 static struct isl_context *isl_context_gbr_alloc(struct isl_basic_set *dom)
3171 {
3172         struct isl_context_gbr *cgbr;
3173
3174         if (!dom)
3175                 return NULL;
3176
3177         cgbr = isl_calloc_type(dom->ctx, struct isl_context_gbr);
3178         if (!cgbr)
3179                 return NULL;
3180
3181         cgbr->context.op = &isl_context_gbr_op;
3182
3183         cgbr->shifted = NULL;
3184         cgbr->cone = NULL;
3185         cgbr->tab = isl_tab_from_basic_set(dom);
3186         cgbr->tab = isl_tab_init_samples(cgbr->tab);
3187         if (!cgbr->tab)
3188                 goto error;
3189         if (isl_tab_track_bset(cgbr->tab,
3190                                 isl_basic_set_cow(isl_basic_set_copy(dom))) < 0)
3191                 goto error;
3192         check_gbr_integer_feasible(cgbr);
3193
3194         return &cgbr->context;
3195 error:
3196         cgbr->context.op->free(&cgbr->context);
3197         return NULL;
3198 }
3199
3200 static struct isl_context *isl_context_alloc(struct isl_basic_set *dom)
3201 {
3202         if (!dom)
3203                 return NULL;
3204
3205         if (dom->ctx->opt->context == ISL_CONTEXT_LEXMIN)
3206                 return isl_context_lex_alloc(dom);
3207         else
3208                 return isl_context_gbr_alloc(dom);
3209 }
3210
3211 /* Construct an isl_sol_map structure for accumulating the solution.
3212  * If track_empty is set, then we also keep track of the parts
3213  * of the context where there is no solution.
3214  * If max is set, then we are solving a maximization, rather than
3215  * a minimization problem, which means that the variables in the
3216  * tableau have value "M - x" rather than "M + x".
3217  */
3218 static struct isl_sol_map *sol_map_init(struct isl_basic_map *bmap,
3219         struct isl_basic_set *dom, int track_empty, int max)
3220 {
3221         struct isl_sol_map *sol_map = NULL;
3222
3223         if (!bmap)
3224                 goto error;
3225
3226         sol_map = isl_calloc_type(bmap->ctx, struct isl_sol_map);
3227         if (!sol_map)
3228                 goto error;
3229
3230         sol_map->sol.rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
3231         sol_map->sol.dec_level.callback.run = &sol_dec_level_wrap;
3232         sol_map->sol.dec_level.sol = &sol_map->sol;
3233         sol_map->sol.max = max;
3234         sol_map->sol.n_out = isl_basic_map_dim(bmap, isl_dim_out);
3235         sol_map->sol.add = &sol_map_add_wrap;
3236         sol_map->sol.add_empty = track_empty ? &sol_map_add_empty_wrap : NULL;
3237         sol_map->sol.free = &sol_map_free_wrap;
3238         sol_map->map = isl_map_alloc_dim(isl_basic_map_get_dim(bmap), 1,
3239                                             ISL_MAP_DISJOINT);
3240         if (!sol_map->map)
3241                 goto error;
3242
3243         sol_map->sol.context = isl_context_alloc(dom);
3244         if (!sol_map->sol.context)
3245                 goto error;
3246
3247         if (track_empty) {
3248                 sol_map->empty = isl_set_alloc_dim(isl_basic_set_get_dim(dom),
3249                                                         1, ISL_SET_DISJOINT);
3250                 if (!sol_map->empty)
3251                         goto error;
3252         }
3253
3254         isl_basic_set_free(dom);
3255         return sol_map;
3256 error:
3257         isl_basic_set_free(dom);
3258         sol_map_free(sol_map);
3259         return NULL;
3260 }
3261
3262 /* Check whether all coefficients of (non-parameter) variables
3263  * are non-positive, meaning that no pivots can be performed on the row.
3264  */
3265 static int is_critical(struct isl_tab *tab, int row)
3266 {
3267         int j;
3268         unsigned off = 2 + tab->M;
3269
3270         for (j = tab->n_dead; j < tab->n_col; ++j) {
3271                 if (tab->col_var[j] >= 0 &&
3272                     (tab->col_var[j] < tab->n_param  ||
3273                     tab->col_var[j] >= tab->n_var - tab->n_div))
3274                         continue;
3275
3276                 if (isl_int_is_pos(tab->mat->row[row][off + j]))
3277                         return 0;
3278         }
3279
3280         return 1;
3281 }
3282
3283 /* Check whether the inequality represented by vec is strict over the integers,
3284  * i.e., there are no integer values satisfying the constraint with
3285  * equality.  This happens if the gcd of the coefficients is not a divisor
3286  * of the constant term.  If so, scale the constraint down by the gcd
3287  * of the coefficients.
3288  */
3289 static int is_strict(struct isl_vec *vec)
3290 {
3291         isl_int gcd;
3292         int strict = 0;
3293
3294         isl_int_init(gcd);
3295         isl_seq_gcd(vec->el + 1, vec->size - 1, &gcd);
3296         if (!isl_int_is_one(gcd)) {
3297                 strict = !isl_int_is_divisible_by(vec->el[0], gcd);
3298                 isl_int_fdiv_q(vec->el[0], vec->el[0], gcd);
3299                 isl_seq_scale_down(vec->el + 1, vec->el + 1, gcd, vec->size-1);
3300         }
3301         isl_int_clear(gcd);
3302
3303         return strict;
3304 }
3305
3306 /* Determine the sign of the given row of the main tableau.
3307  * The result is one of
3308  *      isl_tab_row_pos: always non-negative; no pivot needed
3309  *      isl_tab_row_neg: always non-positive; pivot
3310  *      isl_tab_row_any: can be both positive and negative; split
3311  *
3312  * We first handle some simple cases
3313  *      - the row sign may be known already
3314  *      - the row may be obviously non-negative
3315  *      - the parametric constant may be equal to that of another row
3316  *        for which we know the sign.  This sign will be either "pos" or
3317  *        "any".  If it had been "neg" then we would have pivoted before.
3318  *
3319  * If none of these cases hold, we check the value of the row for each
3320  * of the currently active samples.  Based on the signs of these values
3321  * we make an initial determination of the sign of the row.
3322  *
3323  *      all zero                        ->      unk(nown)
3324  *      all non-negative                ->      pos
3325  *      all non-positive                ->      neg
3326  *      both negative and positive      ->      all
3327  *
3328  * If we end up with "all", we are done.
3329  * Otherwise, we perform a check for positive and/or negative
3330  * values as follows.
3331  *
3332  *      samples        neg             unk             pos
3333  *      <0 ?                        Y        N      Y        N
3334  *                                          pos    any      pos
3335  *      >0 ?         Y      N    Y     N
3336  *                  any    neg  any   neg
3337  *
3338  * There is no special sign for "zero", because we can usually treat zero
3339  * as either non-negative or non-positive, whatever works out best.
3340  * However, if the row is "critical", meaning that pivoting is impossible
3341  * then we don't want to limp zero with the non-positive case, because
3342  * then we we would lose the solution for those values of the parameters
3343  * where the value of the row is zero.  Instead, we treat 0 as non-negative
3344  * ensuring a split if the row can attain both zero and negative values.
3345  * The same happens when the original constraint was one that could not
3346  * be satisfied with equality by any integer values of the parameters.
3347  * In this case, we normalize the constraint, but then a value of zero
3348  * for the normalized constraint is actually a positive value for the
3349  * original constraint, so again we need to treat zero as non-negative.
3350  * In both these cases, we have the following decision tree instead:
3351  *
3352  *      all non-negative                ->      pos
3353  *      all negative                    ->      neg
3354  *      both negative and non-negative  ->      all
3355  *
3356  *      samples        neg                             pos
3357  *      <0 ?                                        Y        N
3358  *                                                 any      pos
3359  *      >=0 ?        Y      N
3360  *                  any    neg
3361  */
3362 static enum isl_tab_row_sign row_sign(struct isl_tab *tab,
3363         struct isl_sol *sol, int row)
3364 {
3365         struct isl_vec *ineq = NULL;
3366         enum isl_tab_row_sign res = isl_tab_row_unknown;
3367         int critical;
3368         int strict;
3369         int row2;
3370
3371         if (tab->row_sign[row] != isl_tab_row_unknown)
3372                 return tab->row_sign[row];
3373         if (is_obviously_nonneg(tab, row))
3374                 return isl_tab_row_pos;
3375         for (row2 = tab->n_redundant; row2 < tab->n_row; ++row2) {
3376                 if (tab->row_sign[row2] == isl_tab_row_unknown)
3377                         continue;
3378                 if (identical_parameter_line(tab, row, row2))
3379                         return tab->row_sign[row2];
3380         }
3381
3382         critical = is_critical(tab, row);
3383
3384         ineq = get_row_parameter_ineq(tab, row);
3385         if (!ineq)
3386                 goto error;
3387
3388         strict = is_strict(ineq);
3389
3390         res = sol->context->op->ineq_sign(sol->context, ineq->el,
3391                                           critical || strict);
3392
3393         if (res == isl_tab_row_unknown || res == isl_tab_row_pos) {
3394                 /* test for negative values */
3395                 int feasible;
3396                 isl_seq_neg(ineq->el, ineq->el, ineq->size);
3397                 isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3398
3399                 feasible = sol->context->op->test_ineq(sol->context, ineq->el);
3400                 if (feasible < 0)
3401                         goto error;
3402                 if (!feasible)
3403                         res = isl_tab_row_pos;
3404                 else
3405                         res = (res == isl_tab_row_unknown) ? isl_tab_row_neg
3406                                                            : isl_tab_row_any;
3407                 if (res == isl_tab_row_neg) {
3408                         isl_seq_neg(ineq->el, ineq->el, ineq->size);
3409                         isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3410                 }
3411         }
3412
3413         if (res == isl_tab_row_neg) {
3414                 /* test for positive values */
3415                 int feasible;
3416                 if (!critical && !strict)
3417                         isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3418
3419                 feasible = sol->context->op->test_ineq(sol->context, ineq->el);
3420                 if (feasible < 0)
3421                         goto error;
3422                 if (feasible)
3423                         res = isl_tab_row_any;
3424         }
3425
3426         isl_vec_free(ineq);
3427         return res;
3428 error:
3429         isl_vec_free(ineq);
3430         return isl_tab_row_unknown;
3431 }
3432
3433 static void find_solutions(struct isl_sol *sol, struct isl_tab *tab);
3434
3435 /* Find solutions for values of the parameters that satisfy the given
3436  * inequality.
3437  *
3438  * We currently take a snapshot of the context tableau that is reset
3439  * when we return from this function, while we make a copy of the main
3440  * tableau, leaving the original main tableau untouched.
3441  * These are fairly arbitrary choices.  Making a copy also of the context
3442  * tableau would obviate the need to undo any changes made to it later,
3443  * while taking a snapshot of the main tableau could reduce memory usage.
3444  * If we were to switch to taking a snapshot of the main tableau,
3445  * we would have to keep in mind that we need to save the row signs
3446  * and that we need to do this before saving the current basis
3447  * such that the basis has been restore before we restore the row signs.
3448  */
3449 static void find_in_pos(struct isl_sol *sol, struct isl_tab *tab, isl_int *ineq)
3450 {
3451         void *saved;
3452
3453         if (!sol->context)
3454                 goto error;
3455         saved = sol->context->op->save(sol->context);
3456
3457         tab = isl_tab_dup(tab);
3458         if (!tab)
3459                 goto error;
3460
3461         sol->context->op->add_ineq(sol->context, ineq, 0, 1);
3462
3463         find_solutions(sol, tab);
3464
3465         if (!sol->error)
3466                 sol->context->op->restore(sol->context, saved);
3467         return;
3468 error:
3469         sol->error = 1;
3470 }
3471
3472 /* Record the absence of solutions for those values of the parameters
3473  * that do not satisfy the given inequality with equality.
3474  */
3475 static void no_sol_in_strict(struct isl_sol *sol,
3476         struct isl_tab *tab, struct isl_vec *ineq)
3477 {
3478         int empty;
3479         void *saved;
3480
3481         if (!sol->context || sol->error)
3482                 goto error;
3483         saved = sol->context->op->save(sol->context);
3484
3485         isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3486
3487         sol->context->op->add_ineq(sol->context, ineq->el, 1, 0);
3488         if (!sol->context)
3489                 goto error;
3490
3491         empty = tab->empty;
3492         tab->empty = 1;
3493         sol_add(sol, tab);
3494         tab->empty = empty;
3495
3496         isl_int_add_ui(ineq->el[0], ineq->el[0], 1);
3497
3498         sol->context->op->restore(sol->context, saved);
3499         return;
3500 error:
3501         sol->error = 1;
3502 }
3503
3504 /* Compute the lexicographic minimum of the set represented by the main
3505  * tableau "tab" within the context "sol->context_tab".
3506  * On entry the sample value of the main tableau is lexicographically
3507  * less than or equal to this lexicographic minimum.
3508  * Pivots are performed until a feasible point is found, which is then
3509  * necessarily equal to the minimum, or until the tableau is found to
3510  * be infeasible.  Some pivots may need to be performed for only some
3511  * feasible values of the context tableau.  If so, the context tableau
3512  * is split into a part where the pivot is needed and a part where it is not.
3513  *
3514  * Whenever we enter the main loop, the main tableau is such that no
3515  * "obvious" pivots need to be performed on it, where "obvious" means
3516  * that the given row can be seen to be negative without looking at
3517  * the context tableau.  In particular, for non-parametric problems,
3518  * no pivots need to be performed on the main tableau.
3519  * The caller of find_solutions is responsible for making this property
3520  * hold prior to the first iteration of the loop, while restore_lexmin
3521  * is called before every other iteration.
3522  *
3523  * Inside the main loop, we first examine the signs of the rows of
3524  * the main tableau within the context of the context tableau.
3525  * If we find a row that is always non-positive for all values of
3526  * the parameters satisfying the context tableau and negative for at
3527  * least one value of the parameters, we perform the appropriate pivot
3528  * and start over.  An exception is the case where no pivot can be
3529  * performed on the row.  In this case, we require that the sign of
3530  * the row is negative for all values of the parameters (rather than just
3531  * non-positive).  This special case is handled inside row_sign, which
3532  * will say that the row can have any sign if it determines that it can
3533  * attain both negative and zero values.
3534  *
3535  * If we can't find a row that always requires a pivot, but we can find
3536  * one or more rows that require a pivot for some values of the parameters
3537  * (i.e., the row can attain both positive and negative signs), then we split
3538  * the context tableau into two parts, one where we force the sign to be
3539  * non-negative and one where we force is to be negative.
3540  * The non-negative part is handled by a recursive call (through find_in_pos).
3541  * Upon returning from this call, we continue with the negative part and
3542  * perform the required pivot.
3543  *
3544  * If no such rows can be found, all rows are non-negative and we have
3545  * found a (rational) feasible point.  If we only wanted a rational point
3546  * then we are done.
3547  * Otherwise, we check if all values of the sample point of the tableau
3548  * are integral for the variables.  If so, we have found the minimal
3549  * integral point and we are done.
3550  * If the sample point is not integral, then we need to make a distinction
3551  * based on whether the constant term is non-integral or the coefficients
3552  * of the parameters.  Furthermore, in order to decide how to handle
3553  * the non-integrality, we also need to know whether the coefficients
3554  * of the other columns in the tableau are integral.  This leads
3555  * to the following table.  The first two rows do not correspond
3556  * to a non-integral sample point and are only mentioned for completeness.
3557  *
3558  *      constant        parameters      other
3559  *
3560  *      int             int             int     |
3561  *      int             int             rat     | -> no problem
3562  *
3563  *      rat             int             int       -> fail
3564  *
3565  *      rat             int             rat       -> cut
3566  *
3567  *      int             rat             rat     |
3568  *      rat             rat             rat     | -> parametric cut
3569  *
3570  *      int             rat             int     |
3571  *      rat             rat             int     | -> split context
3572  *
3573  * If the parametric constant is completely integral, then there is nothing
3574  * to be done.  If the constant term is non-integral, but all the other
3575  * coefficient are integral, then there is nothing that can be done
3576  * and the tableau has no integral solution.
3577  * If, on the other hand, one or more of the other columns have rational
3578  * coeffcients, but the parameter coefficients are all integral, then
3579  * we can perform a regular (non-parametric) cut.
3580  * Finally, if there is any parameter coefficient that is non-integral,
3581  * then we need to involve the context tableau.  There are two cases here.
3582  * If at least one other column has a rational coefficient, then we
3583  * can perform a parametric cut in the main tableau by adding a new
3584  * integer division in the context tableau.
3585  * If all other columns have integral coefficients, then we need to
3586  * enforce that the rational combination of parameters (c + \sum a_i y_i)/m
3587  * is always integral.  We do this by introducing an integer division
3588  * q = floor((c + \sum a_i y_i)/m) and stipulating that its argument should
3589  * always be integral in the context tableau, i.e., m q = c + \sum a_i y_i.
3590  * Since q is expressed in the tableau as
3591  *      c + \sum a_i y_i - m q >= 0
3592  *      -c - \sum a_i y_i + m q + m - 1 >= 0
3593  * it is sufficient to add the inequality
3594  *      -c - \sum a_i y_i + m q >= 0
3595  * In the part of the context where this inequality does not hold, the
3596  * main tableau is marked as being empty.
3597  */
3598 static void find_solutions(struct isl_sol *sol, struct isl_tab *tab)
3599 {
3600         struct isl_context *context;
3601
3602         if (!tab || sol->error)
3603                 goto error;
3604
3605         context = sol->context;
3606
3607         if (tab->empty)
3608                 goto done;
3609         if (context->op->is_empty(context))
3610                 goto done;
3611
3612         for (; tab && !tab->empty; tab = restore_lexmin(tab)) {
3613                 int flags;
3614                 int row;
3615                 enum isl_tab_row_sign sgn;
3616                 int split = -1;
3617                 int n_split = 0;
3618
3619                 for (row = tab->n_redundant; row < tab->n_row; ++row) {
3620                         if (!isl_tab_var_from_row(tab, row)->is_nonneg)
3621                                 continue;
3622                         sgn = row_sign(tab, sol, row);
3623                         if (!sgn)
3624                                 goto error;
3625                         tab->row_sign[row] = sgn;
3626                         if (sgn == isl_tab_row_any)
3627                                 n_split++;
3628                         if (sgn == isl_tab_row_any && split == -1)
3629                                 split = row;
3630                         if (sgn == isl_tab_row_neg)
3631                                 break;
3632                 }
3633                 if (row < tab->n_row)
3634                         continue;
3635                 if (split != -1) {
3636                         struct isl_vec *ineq;
3637                         if (n_split != 1)
3638                                 split = context->op->best_split(context, tab);
3639                         if (split < 0)
3640                                 goto error;
3641                         ineq = get_row_parameter_ineq(tab, split);
3642                         if (!ineq)
3643                                 goto error;
3644                         is_strict(ineq);
3645                         for (row = tab->n_redundant; row < tab->n_row; ++row) {
3646                                 if (!isl_tab_var_from_row(tab, row)->is_nonneg)
3647                                         continue;
3648                                 if (tab->row_sign[row] == isl_tab_row_any)
3649                                         tab->row_sign[row] = isl_tab_row_unknown;
3650                         }
3651                         tab->row_sign[split] = isl_tab_row_pos;
3652                         sol_inc_level(sol);
3653                         find_in_pos(sol, tab, ineq->el);
3654                         tab->row_sign[split] = isl_tab_row_neg;
3655                         row = split;
3656                         isl_seq_neg(ineq->el, ineq->el, ineq->size);
3657                         isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
3658                         if (!sol->error)
3659                                 context->op->add_ineq(context, ineq->el, 0, 1);
3660                         isl_vec_free(ineq);
3661                         if (sol->error)
3662                                 goto error;
3663                         continue;
3664                 }
3665                 if (tab->rational)
3666                         break;
3667                 row = first_non_integer_row(tab, &flags);
3668                 if (row < 0)
3669                         break;
3670                 if (ISL_FL_ISSET(flags, I_PAR)) {
3671                         if (ISL_FL_ISSET(flags, I_VAR)) {
3672                                 if (isl_tab_mark_empty(tab) < 0)
3673                                         goto error;
3674                                 break;
3675                         }
3676                         row = add_cut(tab, row);
3677                 } else if (ISL_FL_ISSET(flags, I_VAR)) {
3678                         struct isl_vec *div;
3679                         struct isl_vec *ineq;
3680                         int d;
3681                         div = get_row_split_div(tab, row);
3682                         if (!div)
3683                                 goto error;
3684                         d = context->op->get_div(context, tab, div);
3685                         isl_vec_free(div);
3686                         if (d < 0)
3687                                 goto error;
3688                         ineq = ineq_for_div(context->op->peek_basic_set(context), d);
3689                         if (!ineq)
3690                                 goto error;
3691                         sol_inc_level(sol);
3692                         no_sol_in_strict(sol, tab, ineq);
3693                         isl_seq_neg(ineq->el, ineq->el, ineq->size);
3694                         context->op->add_ineq(context, ineq->el, 1, 1);
3695                         isl_vec_free(ineq);
3696                         if (sol->error || !context->op->is_ok(context))
3697                                 goto error;
3698                         tab = set_row_cst_to_div(tab, row, d);
3699                         if (context->op->is_empty(context))
3700                                 break;
3701                 } else
3702                         row = add_parametric_cut(tab, row, context);
3703                 if (row < 0)
3704                         goto error;
3705         }
3706 done:
3707         sol_add(sol, tab);
3708         isl_tab_free(tab);
3709         return;
3710 error:
3711         isl_tab_free(tab);
3712         sol->error = 1;
3713 }
3714
3715 /* Compute the lexicographic minimum of the set represented by the main
3716  * tableau "tab" within the context "sol->context_tab".
3717  *
3718  * As a preprocessing step, we first transfer all the purely parametric
3719  * equalities from the main tableau to the context tableau, i.e.,
3720  * parameters that have been pivoted to a row.
3721  * These equalities are ignored by the main algorithm, because the
3722  * corresponding rows may not be marked as being non-negative.
3723  * In parts of the context where the added equality does not hold,
3724  * the main tableau is marked as being empty.
3725  */
3726 static void find_solutions_main(struct isl_sol *sol, struct isl_tab *tab)
3727 {
3728         int row;
3729
3730         if (!tab)
3731                 goto error;
3732
3733         sol->level = 0;
3734
3735         for (row = tab->n_redundant; row < tab->n_row; ++row) {
3736                 int p;
3737                 struct isl_vec *eq;
3738
3739                 if (tab->row_var[row] < 0)
3740                         continue;
3741                 if (tab->row_var[row] >= tab->n_param &&
3742                     tab->row_var[row] < tab->n_var - tab->n_div)
3743                         continue;
3744                 if (tab->row_var[row] < tab->n_param)
3745                         p = tab->row_var[row];
3746                 else
3747                         p = tab->row_var[row]
3748                                 + tab->n_param - (tab->n_var - tab->n_div);
3749
3750                 eq = isl_vec_alloc(tab->mat->ctx, 1+tab->n_param+tab->n_div);
3751                 if (!eq)
3752                         goto error;
3753                 get_row_parameter_line(tab, row, eq->el);
3754                 isl_int_neg(eq->el[1 + p], tab->mat->row[row][0]);
3755                 eq = isl_vec_normalize(eq);
3756
3757                 sol_inc_level(sol);
3758                 no_sol_in_strict(sol, tab, eq);
3759
3760                 isl_seq_neg(eq->el, eq->el, eq->size);
3761                 sol_inc_level(sol);
3762                 no_sol_in_strict(sol, tab, eq);
3763                 isl_seq_neg(eq->el, eq->el, eq->size);
3764
3765                 sol->context->op->add_eq(sol->context, eq->el, 1, 1);
3766
3767                 isl_vec_free(eq);
3768
3769                 if (isl_tab_mark_redundant(tab, row) < 0)
3770                         goto error;
3771
3772                 if (sol->context->op->is_empty(sol->context))
3773                         break;
3774
3775                 row = tab->n_redundant - 1;
3776         }
3777
3778         find_solutions(sol, tab);
3779
3780         sol->level = 0;
3781         sol_pop(sol);
3782
3783         return;
3784 error:
3785         isl_tab_free(tab);
3786         sol->error = 1;
3787 }
3788
3789 static void sol_map_find_solutions(struct isl_sol_map *sol_map,
3790         struct isl_tab *tab)
3791 {
3792         find_solutions_main(&sol_map->sol, tab);
3793 }
3794
3795 /* Check if integer division "div" of "dom" also occurs in "bmap".
3796  * If so, return its position within the divs.
3797  * If not, return -1.
3798  */
3799 static int find_context_div(struct isl_basic_map *bmap,
3800         struct isl_basic_set *dom, unsigned div)
3801 {
3802         int i;
3803         unsigned b_dim = isl_dim_total(bmap->dim);
3804         unsigned d_dim = isl_dim_total(dom->dim);
3805
3806         if (isl_int_is_zero(dom->div[div][0]))
3807                 return -1;
3808         if (isl_seq_first_non_zero(dom->div[div] + 2 + d_dim, dom->n_div) != -1)
3809                 return -1;
3810
3811         for (i = 0; i < bmap->n_div; ++i) {
3812                 if (isl_int_is_zero(bmap->div[i][0]))
3813                         continue;
3814                 if (isl_seq_first_non_zero(bmap->div[i] + 2 + d_dim,
3815                                            (b_dim - d_dim) + bmap->n_div) != -1)
3816                         continue;
3817                 if (isl_seq_eq(bmap->div[i], dom->div[div], 2 + d_dim))
3818                         return i;
3819         }
3820         return -1;
3821 }
3822
3823 /* The correspondence between the variables in the main tableau,
3824  * the context tableau, and the input map and domain is as follows.
3825  * The first n_param and the last n_div variables of the main tableau
3826  * form the variables of the context tableau.
3827  * In the basic map, these n_param variables correspond to the
3828  * parameters and the input dimensions.  In the domain, they correspond
3829  * to the parameters and the set dimensions.
3830  * The n_div variables correspond to the integer divisions in the domain.
3831  * To ensure that everything lines up, we may need to copy some of the
3832  * integer divisions of the domain to the map.  These have to be placed
3833  * in the same order as those in the context and they have to be placed
3834  * after any other integer divisions that the map may have.
3835  * This function performs the required reordering.
3836  */
3837 static struct isl_basic_map *align_context_divs(struct isl_basic_map *bmap,
3838         struct isl_basic_set *dom)
3839 {
3840         int i;
3841         int common = 0;
3842         int other;
3843
3844         for (i = 0; i < dom->n_div; ++i)
3845                 if (find_context_div(bmap, dom, i) != -1)
3846                         common++;
3847         other = bmap->n_div - common;
3848         if (dom->n_div - common > 0) {
3849                 bmap = isl_basic_map_extend_dim(bmap, isl_dim_copy(bmap->dim),
3850                                 dom->n_div - common, 0, 0);
3851                 if (!bmap)
3852                         return NULL;
3853         }
3854         for (i = 0; i < dom->n_div; ++i) {
3855                 int pos = find_context_div(bmap, dom, i);
3856                 if (pos < 0) {
3857                         pos = isl_basic_map_alloc_div(bmap);
3858                         if (pos < 0)
3859                                 goto error;
3860                         isl_int_set_si(bmap->div[pos][0], 0);
3861                 }
3862                 if (pos != other + i)
3863                         isl_basic_map_swap_div(bmap, pos, other + i);
3864         }
3865         return bmap;
3866 error:
3867         isl_basic_map_free(bmap);
3868         return NULL;
3869 }
3870
3871 /* Compute the lexicographic minimum (or maximum if "max" is set)
3872  * of "bmap" over the domain "dom" and return the result as a map.
3873  * If "empty" is not NULL, then *empty is assigned a set that
3874  * contains those parts of the domain where there is no solution.
3875  * If "bmap" is marked as rational (ISL_BASIC_MAP_RATIONAL),
3876  * then we compute the rational optimum.  Otherwise, we compute
3877  * the integral optimum.
3878  *
3879  * We perform some preprocessing.  As the PILP solver does not
3880  * handle implicit equalities very well, we first make sure all
3881  * the equalities are explicitly available.
3882  * We also make sure the divs in the domain are properly order,
3883  * because they will be added one by one in the given order
3884  * during the construction of the solution map.
3885  */
3886 struct isl_map *isl_tab_basic_map_partial_lexopt(
3887                 struct isl_basic_map *bmap, struct isl_basic_set *dom,
3888                 struct isl_set **empty, int max)
3889 {
3890         struct isl_tab *tab;
3891         struct isl_map *result = NULL;
3892         struct isl_sol_map *sol_map = NULL;
3893         struct isl_context *context;
3894         struct isl_basic_map *eq;
3895
3896         if (empty)
3897                 *empty = NULL;
3898         if (!bmap || !dom)
3899                 goto error;
3900
3901         isl_assert(bmap->ctx,
3902             isl_basic_map_compatible_domain(bmap, dom), goto error);
3903
3904         eq = isl_basic_map_copy(bmap);
3905         eq = isl_basic_map_intersect_domain(eq, isl_basic_set_copy(dom));
3906         eq = isl_basic_map_affine_hull(eq);
3907         bmap = isl_basic_map_intersect(bmap, eq);
3908
3909         if (dom->n_div) {
3910                 dom = isl_basic_set_order_divs(dom);
3911                 bmap = align_context_divs(bmap, dom);
3912         }
3913         sol_map = sol_map_init(bmap, dom, !!empty, max);
3914         if (!sol_map)
3915                 goto error;
3916
3917         context = sol_map->sol.context;
3918         if (isl_basic_set_fast_is_empty(context->op->peek_basic_set(context)))
3919                 /* nothing */;
3920         else if (isl_basic_map_fast_is_empty(bmap))
3921                 sol_map_add_empty_if_needed(sol_map,
3922                     isl_basic_set_copy(context->op->peek_basic_set(context)));
3923         else {
3924                 tab = tab_for_lexmin(bmap,
3925                                     context->op->peek_basic_set(context), 1, max);
3926                 tab = context->op->detect_nonnegative_parameters(context, tab);
3927                 sol_map_find_solutions(sol_map, tab);
3928         }
3929         if (sol_map->sol.error)
3930                 goto error;
3931
3932         result = isl_map_copy(sol_map->map);
3933         if (empty)
3934                 *empty = isl_set_copy(sol_map->empty);
3935         sol_free(&sol_map->sol);
3936         isl_basic_map_free(bmap);
3937         return result;
3938 error:
3939         sol_free(&sol_map->sol);
3940         isl_basic_map_free(bmap);
3941         return NULL;
3942 }
3943
3944 struct isl_sol_for {
3945         struct isl_sol  sol;
3946         int             (*fn)(__isl_take isl_basic_set *dom,
3947                                 __isl_take isl_mat *map, void *user);
3948         void            *user;
3949 };
3950
3951 static void sol_for_free(struct isl_sol_for *sol_for)
3952 {
3953         if (sol_for->sol.context)
3954                 sol_for->sol.context->op->free(sol_for->sol.context);
3955         free(sol_for);
3956 }
3957
3958 static void sol_for_free_wrap(struct isl_sol *sol)
3959 {
3960         sol_for_free((struct isl_sol_for *)sol);
3961 }
3962
3963 /* Add the solution identified by the tableau and the context tableau.
3964  *
3965  * See documentation of sol_add for more details.
3966  *
3967  * Instead of constructing a basic map, this function calls a user
3968  * defined function with the current context as a basic set and
3969  * an affine matrix reprenting the relation between the input and output.
3970  * The number of rows in this matrix is equal to one plus the number
3971  * of output variables.  The number of columns is equal to one plus
3972  * the total dimension of the context, i.e., the number of parameters,
3973  * input variables and divs.  Since some of the columns in the matrix
3974  * may refer to the divs, the basic set is not simplified.
3975  * (Simplification may reorder or remove divs.)
3976  */
3977 static void sol_for_add(struct isl_sol_for *sol,
3978         struct isl_basic_set *dom, struct isl_mat *M)
3979 {
3980         if (sol->sol.error || !dom || !M)
3981                 goto error;
3982
3983         dom = isl_basic_set_simplify(dom);
3984         dom = isl_basic_set_finalize(dom);
3985
3986         if (sol->fn(isl_basic_set_copy(dom), isl_mat_copy(M), sol->user) < 0)
3987                 goto error;
3988
3989         isl_basic_set_free(dom);
3990         isl_mat_free(M);
3991         return;
3992 error:
3993         isl_basic_set_free(dom);
3994         isl_mat_free(M);
3995         sol->sol.error = 1;
3996 }
3997
3998 static void sol_for_add_wrap(struct isl_sol *sol,
3999         struct isl_basic_set *dom, struct isl_mat *M)
4000 {
4001         sol_for_add((struct isl_sol_for *)sol, dom, M);
4002 }
4003
4004 static struct isl_sol_for *sol_for_init(struct isl_basic_map *bmap, int max,
4005         int (*fn)(__isl_take isl_basic_set *dom, __isl_take isl_mat *map,
4006                   void *user),
4007         void *user)
4008 {
4009         struct isl_sol_for *sol_for = NULL;
4010         struct isl_dim *dom_dim;
4011         struct isl_basic_set *dom = NULL;
4012
4013         sol_for = isl_calloc_type(bmap->ctx, struct isl_sol_for);
4014         if (!sol_for)
4015                 goto error;
4016
4017         dom_dim = isl_dim_domain(isl_dim_copy(bmap->dim));
4018         dom = isl_basic_set_universe(dom_dim);
4019
4020         sol_for->sol.rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
4021         sol_for->sol.dec_level.callback.run = &sol_dec_level_wrap;
4022         sol_for->sol.dec_level.sol = &sol_for->sol;
4023         sol_for->fn = fn;
4024         sol_for->user = user;
4025         sol_for->sol.max = max;
4026         sol_for->sol.n_out = isl_basic_map_dim(bmap, isl_dim_out);
4027         sol_for->sol.add = &sol_for_add_wrap;
4028         sol_for->sol.add_empty = NULL;
4029         sol_for->sol.free = &sol_for_free_wrap;
4030
4031         sol_for->sol.context = isl_context_alloc(dom);
4032         if (!sol_for->sol.context)
4033                 goto error;
4034
4035         isl_basic_set_free(dom);
4036         return sol_for;
4037 error:
4038         isl_basic_set_free(dom);
4039         sol_for_free(sol_for);
4040         return NULL;
4041 }
4042
4043 static void sol_for_find_solutions(struct isl_sol_for *sol_for,
4044         struct isl_tab *tab)
4045 {
4046         find_solutions_main(&sol_for->sol, tab);
4047 }
4048
4049 int isl_basic_map_foreach_lexopt(__isl_keep isl_basic_map *bmap, int max,
4050         int (*fn)(__isl_take isl_basic_set *dom, __isl_take isl_mat *map,
4051                   void *user),
4052         void *user)
4053 {
4054         struct isl_sol_for *sol_for = NULL;
4055
4056         bmap = isl_basic_map_copy(bmap);
4057         if (!bmap)
4058                 return -1;
4059
4060         bmap = isl_basic_map_detect_equalities(bmap);
4061         sol_for = sol_for_init(bmap, max, fn, user);
4062
4063         if (isl_basic_map_fast_is_empty(bmap))
4064                 /* nothing */;
4065         else {
4066                 struct isl_tab *tab;
4067                 struct isl_context *context = sol_for->sol.context;
4068                 tab = tab_for_lexmin(bmap,
4069                                 context->op->peek_basic_set(context), 1, max);
4070                 tab = context->op->detect_nonnegative_parameters(context, tab);
4071                 sol_for_find_solutions(sol_for, tab);
4072                 if (sol_for->sol.error)
4073                         goto error;
4074         }
4075
4076         sol_free(&sol_for->sol);
4077         isl_basic_map_free(bmap);
4078         return 0;
4079 error:
4080         sol_free(&sol_for->sol);
4081         isl_basic_map_free(bmap);
4082         return -1;
4083 }
4084
4085 int isl_basic_map_foreach_lexmin(__isl_keep isl_basic_map *bmap,
4086         int (*fn)(__isl_take isl_basic_set *dom, __isl_take isl_mat *map,
4087                   void *user),
4088         void *user)
4089 {
4090         return isl_basic_map_foreach_lexopt(bmap, 0, fn, user);
4091 }
4092
4093 int isl_basic_map_foreach_lexmax(__isl_keep isl_basic_map *bmap,
4094         int (*fn)(__isl_take isl_basic_set *dom, __isl_take isl_mat *map,
4095                   void *user),
4096         void *user)
4097 {
4098         return isl_basic_map_foreach_lexopt(bmap, 1, fn, user);
4099 }