isl_tab_pip.c

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