isl_tab.c

   1 #include "isl_mat.h"
   2 #include "isl_map_private.h"
   3 #include "isl_tab.h"
   4 #include "isl_seq.h"
   5
   6 /*
   7  * The implementation of tableaus in this file was inspired by Section 8
   8  * of David Detlefs, Greg Nelson and James B. Saxe, "Simplify: a theorem
   9  * prover for program checking".
  10  */
  11
  12 struct isl_tab *isl_tab_alloc(struct isl_ctx *ctx,
  13         unsigned n_row, unsigned n_var, unsigned M)
  14 {
  15         int i;
  16         struct isl_tab *tab;
  17         unsigned off = 2 + M;
  18
  19         tab = isl_calloc_type(ctx, struct isl_tab);
  20         if (!tab)
  21                 return NULL;
  22         tab->mat = isl_mat_alloc(ctx, n_row, off + n_var);
  23         if (!tab->mat)
  24                 goto error;
  25         tab->var = isl_alloc_array(ctx, struct isl_tab_var, n_var);
  26         if (!tab->var)
  27                 goto error;
  28         tab->con = isl_alloc_array(ctx, struct isl_tab_var, n_row);
  29         if (!tab->con)
  30                 goto error;
  31         tab->col_var = isl_alloc_array(ctx, int, n_var);
  32         if (!tab->col_var)
  33                 goto error;
  34         tab->row_var = isl_alloc_array(ctx, int, n_row);
  35         if (!tab->row_var)
  36                 goto error;
  37         for (i = 0; i < n_var; ++i) {
  38                 tab->var[i].index = i;
  39                 tab->var[i].is_row = 0;
  40                 tab->var[i].is_nonneg = 0;
  41                 tab->var[i].is_zero = 0;
  42                 tab->var[i].is_redundant = 0;
  43                 tab->var[i].frozen = 0;
  44                 tab->var[i].negated = 0;
  45                 tab->col_var[i] = i;
  46         }
  47         tab->n_row = 0;
  48         tab->n_con = 0;
  49         tab->n_eq = 0;
  50         tab->max_con = n_row;
  51         tab->n_col = n_var;
  52         tab->n_var = n_var;
  53         tab->max_var = n_var;
  54         tab->n_param = 0;
  55         tab->n_div = 0;
  56         tab->n_dead = 0;
  57         tab->n_redundant = 0;
  58         tab->need_undo = 0;
  59         tab->rational = 0;
  60         tab->empty = 0;
  61         tab->in_undo = 0;
  62         tab->M = M;
  63         tab->bottom.type = isl_tab_undo_bottom;
  64         tab->bottom.next = NULL;
  65         tab->top = &tab->bottom;
  66         return tab;
  67 error:
  68         isl_tab_free(tab);
  69         return NULL;
  70 }
  71
  72 int isl_tab_extend_cons(struct isl_tab *tab, unsigned n_new)
  73 {
  74         unsigned off = 2 + tab->M;
  75
  76         if (!tab)
  77                 return -1;
  78
  79         if (tab->max_con < tab->n_con + n_new) {
  80                 struct isl_tab_var *con;
  81
  82                 con = isl_realloc_array(tab->mat->ctx, tab->con,
  83                                     struct isl_tab_var, tab->max_con + n_new);
  84                 if (!con)
  85                         return -1;
  86                 tab->con = con;
  87                 tab->max_con += n_new;
  88         }
  89         if (tab->mat->n_row < tab->n_row + n_new) {
  90                 int *row_var;
  91
  92                 tab->mat = isl_mat_extend(tab->mat,
  93                                         tab->n_row + n_new, off + tab->n_col);
  94                 if (!tab->mat)
  95                         return -1;
  96                 row_var = isl_realloc_array(tab->mat->ctx, tab->row_var,
  97                                             int, tab->mat->n_row);
  98                 if (!row_var)
  99                         return -1;
 100                 tab->row_var = row_var;
 101                 if (tab->row_sign) {
 102                         enum isl_tab_row_sign *s;
 103                         s = isl_realloc_array(tab->mat->ctx, tab->row_sign,
 104                                         enum isl_tab_row_sign, tab->mat->n_row);
 105                         if (!s)
 106                                 return -1;
 107                         tab->row_sign = s;
 108                 }
 109         }
 110         return 0;
 111 }
 112
 113 /* Make room for at least n_new extra variables.
 114  * Return -1 if anything went wrong.
 115  */
 116 int isl_tab_extend_vars(struct isl_tab *tab, unsigned n_new)
 117 {
 118         struct isl_tab_var *var;
 119         unsigned off = 2 + tab->M;
 120
 121         if (tab->max_var < tab->n_var + n_new) {
 122                 var = isl_realloc_array(tab->mat->ctx, tab->var,
 123                                     struct isl_tab_var, tab->n_var + n_new);
 124                 if (!var)
 125                         return -1;
 126                 tab->var = var;
 127                 tab->max_var += n_new;
 128         }
 129
 130         if (tab->mat->n_col < off + tab->n_col + n_new) {
 131                 int *p;
 132
 133                 tab->mat = isl_mat_extend(tab->mat,
 134                                     tab->mat->n_row, off + tab->n_col + n_new);
 135                 if (!tab->mat)
 136                         return -1;
 137                 p = isl_realloc_array(tab->mat->ctx, tab->col_var,
 138                                             int, tab->n_col + n_new);
 139                 if (!p)
 140                         return -1;
 141                 tab->col_var = p;
 142         }
 143
 144         return 0;
 145 }
 146
 147 struct isl_tab *isl_tab_extend(struct isl_tab *tab, unsigned n_new)
 148 {
 149         if (isl_tab_extend_cons(tab, n_new) >= 0)
 150                 return tab;
 151
 152         isl_tab_free(tab);
 153         return NULL;
 154 }
 155
 156 static void free_undo(struct isl_tab *tab)
 157 {
 158         struct isl_tab_undo *undo, *next;
 159
 160         for (undo = tab->top; undo && undo != &tab->bottom; undo = next) {
 161                 next = undo->next;
 162                 free(undo);
 163         }
 164         tab->top = undo;
 165 }
 166
 167 void isl_tab_free(struct isl_tab *tab)
 168 {
 169         if (!tab)
 170                 return;
 171         free_undo(tab);
 172         isl_mat_free(tab->mat);
 173         isl_vec_free(tab->dual);
 174         isl_basic_set_free(tab->bset);
 175         free(tab->var);
 176         free(tab->con);
 177         free(tab->row_var);
 178         free(tab->col_var);
 179         free(tab->row_sign);
 180         isl_mat_free(tab->samples);
 181         free(tab);
 182 }
 183
 184 struct isl_tab *isl_tab_dup(struct isl_tab *tab)
 185 {
 186         int i;
 187         struct isl_tab *dup;
 188         unsigned off;
 189
 190         if (!tab)
 191                 return NULL;
 192
 193         off = 2 + tab->M;
 194         dup = isl_calloc_type(tab->ctx, struct isl_tab);
 195         if (!dup)
 196                 return NULL;
 197         dup->mat = isl_mat_dup(tab->mat);
 198         if (!dup->mat)
 199                 goto error;
 200         dup->var = isl_alloc_array(tab->ctx, struct isl_tab_var, tab->max_var);
 201         if (!dup->var)
 202                 goto error;
 203         for (i = 0; i < tab->n_var; ++i)
 204                 dup->var[i] = tab->var[i];
 205         dup->con = isl_alloc_array(tab->ctx, struct isl_tab_var, tab->max_con);
 206         if (!dup->con)
 207                 goto error;
 208         for (i = 0; i < tab->n_con; ++i)
 209                 dup->con[i] = tab->con[i];
 210         dup->col_var = isl_alloc_array(tab->ctx, int, tab->mat->n_col - off);
 211         if (!dup->col_var)
 212                 goto error;
 213         for (i = 0; i < tab->n_col; ++i)
 214                 dup->col_var[i] = tab->col_var[i];
 215         dup->row_var = isl_alloc_array(tab->ctx, int, tab->mat->n_row);
 216         if (!dup->row_var)
 217                 goto error;
 218         for (i = 0; i < tab->n_row; ++i)
 219                 dup->row_var[i] = tab->row_var[i];
 220         if (tab->row_sign) {
 221                 dup->row_sign = isl_alloc_array(tab->ctx, enum isl_tab_row_sign,
 222                                                 tab->mat->n_row);
 223                 if (!dup->row_sign)
 224                         goto error;
 225                 for (i = 0; i < tab->n_row; ++i)
 226                         dup->row_sign[i] = tab->row_sign[i];
 227         }
 228         if (tab->samples) {
 229                 dup->samples = isl_mat_dup(tab->samples);
 230                 if (!dup->samples)
 231                         goto error;
 232                 dup->n_sample = tab->n_sample;
 233                 dup->n_outside = tab->n_outside;
 234         }
 235         dup->n_row = tab->n_row;
 236         dup->n_con = tab->n_con;
 237         dup->n_eq = tab->n_eq;
 238         dup->max_con = tab->max_con;
 239         dup->n_col = tab->n_col;
 240         dup->n_var = tab->n_var;
 241         dup->max_var = tab->max_var;
 242         dup->n_param = tab->n_param;
 243         dup->n_div = tab->n_div;
 244         dup->n_dead = tab->n_dead;
 245         dup->n_redundant = tab->n_redundant;
 246         dup->rational = tab->rational;
 247         dup->empty = tab->empty;
 248         dup->need_undo = 0;
 249         dup->in_undo = 0;
 250         dup->M = tab->M;
 251         dup->bottom.type = isl_tab_undo_bottom;
 252         dup->bottom.next = NULL;
 253         dup->top = &dup->bottom;
 254         return dup;
 255 error:
 256         isl_tab_free(dup);
 257         return NULL;
 258 }
 259
 260 static struct isl_tab_var *var_from_index(struct isl_tab *tab, int i)
 261 {
 262         if (i >= 0)
 263                 return &tab->var[i];
 264         else
 265                 return &tab->con[~i];
 266 }
 267
 268 struct isl_tab_var *isl_tab_var_from_row(struct isl_tab *tab, int i)
 269 {
 270         return var_from_index(tab, tab->row_var[i]);
 271 }
 272
 273 static struct isl_tab_var *var_from_col(struct isl_tab *tab, int i)
 274 {
 275         return var_from_index(tab, tab->col_var[i]);
 276 }
 277
 278 /* Check if there are any upper bounds on column variable "var",
 279  * i.e., non-negative rows where var appears with a negative coefficient.
 280  * Return 1 if there are no such bounds.
 281  */
 282 static int max_is_manifestly_unbounded(struct isl_tab *tab,
 283         struct isl_tab_var *var)
 284 {
 285         int i;
 286         unsigned off = 2 + tab->M;
 287
 288         if (var->is_row)
 289                 return 0;
 290         for (i = tab->n_redundant; i < tab->n_row; ++i) {
 291                 if (!isl_int_is_neg(tab->mat->row[i][off + var->index]))
 292                         continue;
 293                 if (isl_tab_var_from_row(tab, i)->is_nonneg)
 294                         return 0;
 295         }
 296         return 1;
 297 }
 298
 299 /* Check if there are any lower bounds on column variable "var",
 300  * i.e., non-negative rows where var appears with a positive coefficient.
 301  * Return 1 if there are no such bounds.
 302  */
 303 static int min_is_manifestly_unbounded(struct isl_tab *tab,
 304         struct isl_tab_var *var)
 305 {
 306         int i;
 307         unsigned off = 2 + tab->M;
 308
 309         if (var->is_row)
 310                 return 0;
 311         for (i = tab->n_redundant; i < tab->n_row; ++i) {
 312                 if (!isl_int_is_pos(tab->mat->row[i][off + var->index]))
 313                         continue;
 314                 if (isl_tab_var_from_row(tab, i)->is_nonneg)
 315                         return 0;
 316         }
 317         return 1;
 318 }
 319
 320 static int row_cmp(struct isl_tab *tab, int r1, int r2, int c, isl_int t)
 321 {
 322         unsigned off = 2 + tab->M;
 323
 324         if (tab->M) {
 325                 int s;
 326                 isl_int_mul(t, tab->mat->row[r1][2], tab->mat->row[r2][off+c]);
 327                 isl_int_submul(t, tab->mat->row[r2][2], tab->mat->row[r1][off+c]);
 328                 s = isl_int_sgn(t);
 329                 if (s)
 330                         return s;
 331         }
 332         isl_int_mul(t, tab->mat->row[r1][1], tab->mat->row[r2][off + c]);
 333         isl_int_submul(t, tab->mat->row[r2][1], tab->mat->row[r1][off + c]);
 334         return isl_int_sgn(t);
 335 }
 336
 337 /* Given the index of a column "c", return the index of a row
 338  * that can be used to pivot the column in, with either an increase
 339  * (sgn > 0) or a decrease (sgn < 0) of the corresponding variable.
 340  * If "var" is not NULL, then the row returned will be different from
 341  * the one associated with "var".
 342  *
 343  * Each row in the tableau is of the form
 344  *
 345  *      x_r = a_r0 + \sum_i a_ri x_i
 346  *
 347  * Only rows with x_r >= 0 and with the sign of a_ri opposite to "sgn"
 348  * impose any limit on the increase or decrease in the value of x_c
 349  * and this bound is equal to a_r0 / |a_rc|.  We are therefore looking
 350  * for the row with the smallest (most stringent) such bound.
 351  * Note that the common denominator of each row drops out of the fraction.
 352  * To check if row j has a smaller bound than row r, i.e.,
 353  * a_j0 / |a_jc| < a_r0 / |a_rc| or a_j0 |a_rc| < a_r0 |a_jc|,
 354  * we check if -sign(a_jc) (a_j0 a_rc - a_r0 a_jc) < 0,
 355  * where -sign(a_jc) is equal to "sgn".
 356  */
 357 static int pivot_row(struct isl_tab *tab,
 358         struct isl_tab_var *var, int sgn, int c)
 359 {
 360         int j, r, tsgn;
 361         isl_int t;
 362         unsigned off = 2 + tab->M;
 363
 364         isl_int_init(t);
 365         r = -1;
 366         for (j = tab->n_redundant; j < tab->n_row; ++j) {
 367                 if (var && j == var->index)
 368                         continue;
 369                 if (!isl_tab_var_from_row(tab, j)->is_nonneg)
 370                         continue;
 371                 if (sgn * isl_int_sgn(tab->mat->row[j][off + c]) >= 0)
 372                         continue;
 373                 if (r < 0) {
 374                         r = j;
 375                         continue;
 376                 }
 377                 tsgn = sgn * row_cmp(tab, r, j, c, t);
 378                 if (tsgn < 0 || (tsgn == 0 &&
 379                                             tab->row_var[j] < tab->row_var[r]))
 380                         r = j;
 381         }
 382         isl_int_clear(t);
 383         return r;
 384 }
 385
 386 /* Find a pivot (row and col) that will increase (sgn > 0) or decrease
 387  * (sgn < 0) the value of row variable var.
 388  * If not NULL, then skip_var is a row variable that should be ignored
 389  * while looking for a pivot row.  It is usually equal to var.
 390  *
 391  * As the given row in the tableau is of the form
 392  *
 393  *      x_r = a_r0 + \sum_i a_ri x_i
 394  *
 395  * we need to find a column such that the sign of a_ri is equal to "sgn"
 396  * (such that an increase in x_i will have the desired effect) or a
 397  * column with a variable that may attain negative values.
 398  * If a_ri is positive, then we need to move x_i in the same direction
 399  * to obtain the desired effect.  Otherwise, x_i has to move in the
 400  * opposite direction.
 401  */
 402 static void find_pivot(struct isl_tab *tab,
 403         struct isl_tab_var *var, struct isl_tab_var *skip_var,
 404         int sgn, int *row, int *col)
 405 {
 406         int j, r, c;
 407         isl_int *tr;
 408
 409         *row = *col = -1;
 410
 411         isl_assert(tab->mat->ctx, var->is_row, return);
 412         tr = tab->mat->row[var->index] + 2 + tab->M;
 413
 414         c = -1;
 415         for (j = tab->n_dead; j < tab->n_col; ++j) {
 416                 if (isl_int_is_zero(tr[j]))
 417                         continue;
 418                 if (isl_int_sgn(tr[j]) != sgn &&
 419                     var_from_col(tab, j)->is_nonneg)
 420                         continue;
 421                 if (c < 0 || tab->col_var[j] < tab->col_var[c])
 422                         c = j;
 423         }
 424         if (c < 0)
 425                 return;
 426
 427         sgn *= isl_int_sgn(tr[c]);
 428         r = pivot_row(tab, skip_var, sgn, c);
 429         *row = r < 0 ? var->index : r;
 430         *col = c;
 431 }
 432
 433 /* Return 1 if row "row" represents an obviously redundant inequality.
 434  * This means
 435  *      - it represents an inequality or a variable
 436  *      - that is the sum of a non-negative sample value and a positive
 437  *        combination of zero or more non-negative variables.
 438  */
 439 int isl_tab_row_is_redundant(struct isl_tab *tab, int row)
 440 {
 441         int i;
 442         unsigned off = 2 + tab->M;
 443
 444         if (tab->row_var[row] < 0 && !isl_tab_var_from_row(tab, row)->is_nonneg)
 445                 return 0;
 446
 447         if (isl_int_is_neg(tab->mat->row[row][1]))
 448                 return 0;
 449         if (tab->M && isl_int_is_neg(tab->mat->row[row][2]))
 450                 return 0;
 451
 452         for (i = tab->n_dead; i < tab->n_col; ++i) {
 453                 if (isl_int_is_zero(tab->mat->row[row][off + i]))
 454                         continue;
 455                 if (isl_int_is_neg(tab->mat->row[row][off + i]))
 456                         return 0;
 457                 if (!var_from_col(tab, i)->is_nonneg)
 458                         return 0;
 459         }
 460         return 1;
 461 }
 462
 463 static void swap_rows(struct isl_tab *tab, int row1, int row2)
 464 {
 465         int t;
 466         t = tab->row_var[row1];
 467         tab->row_var[row1] = tab->row_var[row2];
 468         tab->row_var[row2] = t;
 469         isl_tab_var_from_row(tab, row1)->index = row1;
 470         isl_tab_var_from_row(tab, row2)->index = row2;
 471         tab->mat = isl_mat_swap_rows(tab->mat, row1, row2);
 472
 473         if (!tab->row_sign)
 474                 return;
 475         t = tab->row_sign[row1];
 476         tab->row_sign[row1] = tab->row_sign[row2];
 477         tab->row_sign[row2] = t;
 478 }
 479
 480 static void push_union(struct isl_tab *tab,
 481         enum isl_tab_undo_type type, union isl_tab_undo_val u)
 482 {
 483         struct isl_tab_undo *undo;
 484
 485         if (!tab->need_undo)
 486                 return;
 487
 488         undo = isl_alloc_type(tab->mat->ctx, struct isl_tab_undo);
 489         if (!undo) {
 490                 free_undo(tab);
 491                 tab->top = NULL;
 492                 return;
 493         }
 494         undo->type = type;
 495         undo->u = u;
 496         undo->next = tab->top;
 497         tab->top = undo;
 498 }
 499
 500 void isl_tab_push_var(struct isl_tab *tab,
 501         enum isl_tab_undo_type type, struct isl_tab_var *var)
 502 {
 503         union isl_tab_undo_val u;
 504         if (var->is_row)
 505                 u.var_index = tab->row_var[var->index];
 506         else
 507                 u.var_index = tab->col_var[var->index];
 508         push_union(tab, type, u);
 509 }
 510
 511 void isl_tab_push(struct isl_tab *tab, enum isl_tab_undo_type type)
 512 {
 513         union isl_tab_undo_val u = { 0 };
 514         push_union(tab, type, u);
 515 }
 516
 517 /* Push a record on the undo stack describing the current basic
 518  * variables, so that the this state can be restored during rollback.
 519  */
 520 void isl_tab_push_basis(struct isl_tab *tab)
 521 {
 522         int i;
 523         union isl_tab_undo_val u;
 524
 525         u.col_var = isl_alloc_array(tab->mat->ctx, int, tab->n_col);
 526         if (!u.col_var) {
 527                 free_undo(tab);
 528                 tab->top = NULL;
 529                 return;
 530         }
 531         for (i = 0; i < tab->n_col; ++i)
 532                 u.col_var[i] = tab->col_var[i];
 533         push_union(tab, isl_tab_undo_saved_basis, u);
 534 }
 535
 536 /* Mark row with index "row" as being redundant.
 537  * If we may need to undo the operation or if the row represents
 538  * a variable of the original problem, the row is kept,
 539  * but no longer considered when looking for a pivot row.
 540  * Otherwise, the row is simply removed.
 541  *
 542  * The row may be interchanged with some other row.  If it
 543  * is interchanged with a later row, return 1.  Otherwise return 0.
 544  * If the rows are checked in order in the calling function,
 545  * then a return value of 1 means that the row with the given
 546  * row number may now contain a different row that hasn't been checked yet.
 547  */
 548 int isl_tab_mark_redundant(struct isl_tab *tab, int row)
 549 {
 550         struct isl_tab_var *var = isl_tab_var_from_row(tab, row);
 551         var->is_redundant = 1;
 552         isl_assert(tab->mat->ctx, row >= tab->n_redundant, return -1);
 553         if (tab->need_undo || tab->row_var[row] >= 0) {
 554                 if (tab->row_var[row] >= 0 && !var->is_nonneg) {
 555                         var->is_nonneg = 1;
 556                         isl_tab_push_var(tab, isl_tab_undo_nonneg, var);
 557                 }
 558                 if (row != tab->n_redundant)
 559                         swap_rows(tab, row, tab->n_redundant);
 560                 isl_tab_push_var(tab, isl_tab_undo_redundant, var);
 561                 tab->n_redundant++;
 562                 return 0;
 563         } else {
 564                 if (row != tab->n_row - 1)
 565                         swap_rows(tab, row, tab->n_row - 1);
 566                 isl_tab_var_from_row(tab, tab->n_row - 1)->index = -1;
 567                 tab->n_row--;
 568                 return 1;
 569         }
 570 }
 571
 572 struct isl_tab *isl_tab_mark_empty(struct isl_tab *tab)
 573 {
 574         if (!tab->empty && tab->need_undo)
 575                 isl_tab_push(tab, isl_tab_undo_empty);
 576         tab->empty = 1;
 577         return tab;
 578 }
 579
 580 /* Update the rows signs after a pivot of "row" and "col", with "row_sgn"
 581  * the original sign of the pivot element.
 582  * We only keep track of row signs during PILP solving and in this case
 583  * we only pivot a row with negative sign (meaning the value is always
 584  * non-positive) using a positive pivot element.
 585  *
 586  * For each row j, the new value of the parametric constant is equal to
 587  *
 588  *      a_j0 - a_jc a_r0/a_rc
 589  *
 590  * where a_j0 is the original parametric constant, a_rc is the pivot element,
 591  * a_r0 is the parametric constant of the pivot row and a_jc is the
 592  * pivot column entry of the row j.
 593  * Since a_r0 is non-positive and a_rc is positive, the sign of row j
 594  * remains the same if a_jc has the same sign as the row j or if
 595  * a_jc is zero.  In all other cases, we reset the sign to "unknown".
 596  */
 597 static void update_row_sign(struct isl_tab *tab, int row, int col, int row_sgn)
 598 {
 599         int i;
 600         struct isl_mat *mat = tab->mat;
 601         unsigned off = 2 + tab->M;
 602
 603         if (!tab->row_sign)
 604                 return;
 605
 606         if (tab->row_sign[row] == 0)
 607                 return;
 608         isl_assert(mat->ctx, row_sgn > 0, return);
 609         isl_assert(mat->ctx, tab->row_sign[row] == isl_tab_row_neg, return);
 610         tab->row_sign[row] = isl_tab_row_pos;
 611         for (i = 0; i < tab->n_row; ++i) {
 612                 int s;
 613                 if (i == row)
 614                         continue;
 615                 s = isl_int_sgn(mat->row[i][off + col]);
 616                 if (!s)
 617                         continue;
 618                 if (!tab->row_sign[i])
 619                         continue;
 620                 if (s < 0 && tab->row_sign[i] == isl_tab_row_neg)
 621                         continue;
 622                 if (s > 0 && tab->row_sign[i] == isl_tab_row_pos)
 623                         continue;
 624                 tab->row_sign[i] = isl_tab_row_unknown;
 625         }
 626 }
 627
 628 /* Given a row number "row" and a column number "col", pivot the tableau
 629  * such that the associated variables are interchanged.
 630  * The given row in the tableau expresses
 631  *
 632  *      x_r = a_r0 + \sum_i a_ri x_i
 633  *
 634  * or
 635  *
 636  *      x_c = 1/a_rc x_r - a_r0/a_rc + sum_{i \ne r} -a_ri/a_rc
 637  *
 638  * Substituting this equality into the other rows
 639  *
 640  *      x_j = a_j0 + \sum_i a_ji x_i
 641  *
 642  * with a_jc \ne 0, we obtain
 643  *
 644  *      x_j = a_jc/a_rc x_r + a_j0 - a_jc a_r0/a_rc + sum a_ji - a_jc a_ri/a_rc
 645  *
 646  * The tableau
 647  *
 648  *      n_rc/d_r                n_ri/d_r
 649  *      n_jc/d_j                n_ji/d_j
 650  *
 651  * where i is any other column and j is any other row,
 652  * is therefore transformed into
 653  *
 654  * s(n_rc)d_r/|n_rc|            -s(n_rc)n_ri/|n_rc|
 655  * s(n_rc)d_r n_jc/(|n_rc| d_j) (n_ji |n_rc| - s(n_rc)n_jc n_ri)/(|n_rc| d_j)
 656  *
 657  * The transformation is performed along the following steps
 658  *
 659  *      d_r/n_rc                n_ri/n_rc
 660  *      n_jc/d_j                n_ji/d_j
 661  *
 662  *      s(n_rc)d_r/|n_rc|       -s(n_rc)n_ri/|n_rc|
 663  *      n_jc/d_j                n_ji/d_j
 664  *
 665  *      s(n_rc)d_r/|n_rc|       -s(n_rc)n_ri/|n_rc|
 666  *      n_jc/(|n_rc| d_j)       n_ji/(|n_rc| d_j)
 667  *
 668  *      s(n_rc)d_r/|n_rc|       -s(n_rc)n_ri/|n_rc|
 669  *      n_jc/(|n_rc| d_j)       (n_ji |n_rc|)/(|n_rc| d_j)
 670  *
 671  *      s(n_rc)d_r/|n_rc|       -s(n_rc)n_ri/|n_rc|
 672  *      n_jc/(|n_rc| d_j)       (n_ji |n_rc| - s(n_rc)n_jc n_ri)/(|n_rc| d_j)
 673  *
 674  * s(n_rc)d_r/|n_rc|            -s(n_rc)n_ri/|n_rc|
 675  * s(n_rc)d_r n_jc/(|n_rc| d_j) (n_ji |n_rc| - s(n_rc)n_jc n_ri)/(|n_rc| d_j)
 676  *
 677  */
 678 void isl_tab_pivot(struct isl_tab *tab, int row, int col)
 679 {
 680         int i, j;
 681         int sgn;
 682         int t;
 683         struct isl_mat *mat = tab->mat;
 684         struct isl_tab_var *var;
 685         unsigned off = 2 + tab->M;
 686
 687         isl_int_swap(mat->row[row][0], mat->row[row][off + col]);
 688         sgn = isl_int_sgn(mat->row[row][0]);
 689         if (sgn < 0) {
 690                 isl_int_neg(mat->row[row][0], mat->row[row][0]);
 691                 isl_int_neg(mat->row[row][off + col], mat->row[row][off + col]);
 692         } else
 693                 for (j = 0; j < off - 1 + tab->n_col; ++j) {
 694                         if (j == off - 1 + col)
 695                                 continue;
 696                         isl_int_neg(mat->row[row][1 + j], mat->row[row][1 + j]);
 697                 }
 698         if (!isl_int_is_one(mat->row[row][0]))
 699                 isl_seq_normalize(mat->ctx, mat->row[row], off + tab->n_col);
 700         for (i = 0; i < tab->n_row; ++i) {
 701                 if (i == row)
 702                         continue;
 703                 if (isl_int_is_zero(mat->row[i][off + col]))
 704                         continue;
 705                 isl_int_mul(mat->row[i][0], mat->row[i][0], mat->row[row][0]);
 706                 for (j = 0; j < off - 1 + tab->n_col; ++j) {
 707                         if (j == off - 1 + col)
 708                                 continue;
 709                         isl_int_mul(mat->row[i][1 + j],
 710                                     mat->row[i][1 + j], mat->row[row][0]);
 711                         isl_int_addmul(mat->row[i][1 + j],
 712                                     mat->row[i][off + col], mat->row[row][1 + j]);
 713                 }
 714                 isl_int_mul(mat->row[i][off + col],
 715                             mat->row[i][off + col], mat->row[row][off + col]);
 716                 if (!isl_int_is_one(mat->row[i][0]))
 717                         isl_seq_normalize(mat->ctx, mat->row[i], off + tab->n_col);
 718         }
 719         t = tab->row_var[row];
 720         tab->row_var[row] = tab->col_var[col];
 721         tab->col_var[col] = t;
 722         var = isl_tab_var_from_row(tab, row);
 723         var->is_row = 1;
 724         var->index = row;
 725         var = var_from_col(tab, col);
 726         var->is_row = 0;
 727         var->index = col;
 728         update_row_sign(tab, row, col, sgn);
 729         if (tab->in_undo)
 730                 return;
 731         for (i = tab->n_redundant; i < tab->n_row; ++i) {
 732                 if (isl_int_is_zero(mat->row[i][off + col]))
 733                         continue;
 734                 if (!isl_tab_var_from_row(tab, i)->frozen &&
 735                     isl_tab_row_is_redundant(tab, i))
 736                         if (isl_tab_mark_redundant(tab, i))
 737                                 --i;
 738         }
 739 }
 740
 741 /* If "var" represents a column variable, then pivot is up (sgn > 0)
 742  * or down (sgn < 0) to a row.  The variable is assumed not to be
 743  * unbounded in the specified direction.
 744  * If sgn = 0, then the variable is unbounded in both directions,
 745  * and we pivot with any row we can find.
 746  */
 747 static void to_row(struct isl_tab *tab, struct isl_tab_var *var, int sign)
 748 {
 749         int r;
 750         unsigned off = 2 + tab->M;
 751
 752         if (var->is_row)
 753                 return;
 754
 755         if (sign == 0) {
 756                 for (r = tab->n_redundant; r < tab->n_row; ++r)
 757                         if (!isl_int_is_zero(tab->mat->row[r][off+var->index]))
 758                                 break;
 759                 isl_assert(tab->mat->ctx, r < tab->n_row, return);
 760         } else {
 761                 r = pivot_row(tab, NULL, sign, var->index);
 762                 isl_assert(tab->mat->ctx, r >= 0, return);
 763         }
 764
 765         isl_tab_pivot(tab, r, var->index);
 766 }
 767
 768 static void check_table(struct isl_tab *tab)
 769 {
 770         int i;
 771
 772         if (tab->empty)
 773                 return;
 774         for (i = 0; i < tab->n_row; ++i) {
 775                 if (!isl_tab_var_from_row(tab, i)->is_nonneg)
 776                         continue;
 777                 assert(!isl_int_is_neg(tab->mat->row[i][1]));
 778         }
 779 }
 780
 781 /* Return the sign of the maximal value of "var".
 782  * If the sign is not negative, then on return from this function,
 783  * the sample value will also be non-negative.
 784  *
 785  * If "var" is manifestly unbounded wrt positive values, we are done.
 786  * Otherwise, we pivot the variable up to a row if needed
 787  * Then we continue pivoting down until either
 788  *      - no more down pivots can be performed
 789  *      - the sample value is positive
 790  *      - the variable is pivoted into a manifestly unbounded column
 791  */
 792 static int sign_of_max(struct isl_tab *tab, struct isl_tab_var *var)
 793 {
 794         int row, col;
 795
 796         if (max_is_manifestly_unbounded(tab, var))
 797                 return 1;
 798         to_row(tab, var, 1);
 799         while (!isl_int_is_pos(tab->mat->row[var->index][1])) {
 800                 find_pivot(tab, var, var, 1, &row, &col);
 801                 if (row == -1)
 802                         return isl_int_sgn(tab->mat->row[var->index][1]);
 803                 isl_tab_pivot(tab, row, col);
 804                 if (!var->is_row) /* manifestly unbounded */
 805                         return 1;
 806         }
 807         return 1;
 808 }
 809
 810 static int row_is_neg(struct isl_tab *tab, int row)
 811 {
 812         if (!tab->M)
 813                 return isl_int_is_neg(tab->mat->row[row][1]);
 814         if (isl_int_is_pos(tab->mat->row[row][2]))
 815                 return 0;
 816         if (isl_int_is_neg(tab->mat->row[row][2]))
 817                 return 1;
 818         return isl_int_is_neg(tab->mat->row[row][1]);
 819 }
 820
 821 static int row_sgn(struct isl_tab *tab, int row)
 822 {
 823         if (!tab->M)
 824                 return isl_int_sgn(tab->mat->row[row][1]);
 825         if (!isl_int_is_zero(tab->mat->row[row][2]))
 826                 return isl_int_sgn(tab->mat->row[row][2]);
 827         else
 828                 return isl_int_sgn(tab->mat->row[row][1]);
 829 }
 830
 831 /* Perform pivots until the row variable "var" has a non-negative
 832  * sample value or until no more upward pivots can be performed.
 833  * Return the sign of the sample value after the pivots have been
 834  * performed.
 835  */
 836 static int restore_row(struct isl_tab *tab, struct isl_tab_var *var)
 837 {
 838         int row, col;
 839
 840         while (row_is_neg(tab, var->index)) {
 841                 find_pivot(tab, var, var, 1, &row, &col);
 842                 if (row == -1)
 843                         break;
 844                 isl_tab_pivot(tab, row, col);
 845                 if (!var->is_row) /* manifestly unbounded */
 846                         return 1;
 847         }
 848         return row_sgn(tab, var->index);
 849 }
 850
 851 /* Perform pivots until we are sure that the row variable "var"
 852  * can attain non-negative values.  After return from this
 853  * function, "var" is still a row variable, but its sample
 854  * value may not be non-negative, even if the function returns 1.
 855  */
 856 static int at_least_zero(struct isl_tab *tab, struct isl_tab_var *var)
 857 {
 858         int row, col;
 859
 860         while (isl_int_is_neg(tab->mat->row[var->index][1])) {
 861                 find_pivot(tab, var, var, 1, &row, &col);
 862                 if (row == -1)
 863                         break;
 864                 if (row == var->index) /* manifestly unbounded */
 865                         return 1;
 866                 isl_tab_pivot(tab, row, col);
 867         }
 868         return !isl_int_is_neg(tab->mat->row[var->index][1]);
 869 }
 870
 871 /* Return a negative value if "var" can attain negative values.
 872  * Return a non-negative value otherwise.
 873  *
 874  * If "var" is manifestly unbounded wrt negative values, we are done.
 875  * Otherwise, if var is in a column, we can pivot it down to a row.
 876  * Then we continue pivoting down until either
 877  *      - the pivot would result in a manifestly unbounded column
 878  *        => we don't perform the pivot, but simply return -1
 879  *      - no more down pivots can be performed
 880  *      - the sample value is negative
 881  * If the sample value becomes negative and the variable is supposed
 882  * to be nonnegative, then we undo the last pivot.
 883  * However, if the last pivot has made the pivoting variable
 884  * obviously redundant, then it may have moved to another row.
 885  * In that case we look for upward pivots until we reach a non-negative
 886  * value again.
 887  */
 888 static int sign_of_min(struct isl_tab *tab, struct isl_tab_var *var)
 889 {
 890         int row, col;
 891         struct isl_tab_var *pivot_var = NULL;
 892
 893         if (min_is_manifestly_unbounded(tab, var))
 894                 return -1;
 895         if (!var->is_row) {
 896                 col = var->index;
 897                 row = pivot_row(tab, NULL, -1, col);
 898                 pivot_var = var_from_col(tab, col);
 899                 isl_tab_pivot(tab, row, col);
 900                 if (var->is_redundant)
 901                         return 0;
 902                 if (isl_int_is_neg(tab->mat->row[var->index][1])) {
 903                         if (var->is_nonneg) {
 904                                 if (!pivot_var->is_redundant &&
 905                                     pivot_var->index == row)
 906                                         isl_tab_pivot(tab, row, col);
 907                                 else
 908                                         restore_row(tab, var);
 909                         }
 910                         return -1;
 911                 }
 912         }
 913         if (var->is_redundant)
 914                 return 0;
 915         while (!isl_int_is_neg(tab->mat->row[var->index][1])) {
 916                 find_pivot(tab, var, var, -1, &row, &col);
 917                 if (row == var->index)
 918                         return -1;
 919                 if (row == -1)
 920                         return isl_int_sgn(tab->mat->row[var->index][1]);
 921                 pivot_var = var_from_col(tab, col);
 922                 isl_tab_pivot(tab, row, col);
 923                 if (var->is_redundant)
 924                         return 0;
 925         }
 926         if (pivot_var && var->is_nonneg) {
 927                 /* pivot back to non-negative value */
 928                 if (!pivot_var->is_redundant && pivot_var->index == row)
 929                         isl_tab_pivot(tab, row, col);
 930                 else
 931                         restore_row(tab, var);
 932         }
 933         return -1;
 934 }
 935
 936 static int row_at_most_neg_one(struct isl_tab *tab, int row)
 937 {
 938         if (tab->M) {
 939                 if (isl_int_is_pos(tab->mat->row[row][2]))
 940                         return 0;
 941                 if (isl_int_is_neg(tab->mat->row[row][2]))
 942                         return 1;
 943         }
 944         return isl_int_is_neg(tab->mat->row[row][1]) &&
 945                isl_int_abs_ge(tab->mat->row[row][1],
 946                               tab->mat->row[row][0]);
 947 }
 948
 949 /* Return 1 if "var" can attain values <= -1.
 950  * Return 0 otherwise.
 951  *
 952  * The sample value of "var" is assumed to be non-negative when the
 953  * the function is called and will be made non-negative again before
 954  * the function returns.
 955  */
 956 int isl_tab_min_at_most_neg_one(struct isl_tab *tab, struct isl_tab_var *var)
 957 {
 958         int row, col;
 959         struct isl_tab_var *pivot_var;
 960
 961         if (min_is_manifestly_unbounded(tab, var))
 962                 return 1;
 963         if (!var->is_row) {
 964                 col = var->index;
 965                 row = pivot_row(tab, NULL, -1, col);
 966                 pivot_var = var_from_col(tab, col);
 967                 isl_tab_pivot(tab, row, col);
 968                 if (var->is_redundant)
 969                         return 0;
 970                 if (row_at_most_neg_one(tab, var->index)) {
 971                         if (var->is_nonneg) {
 972                                 if (!pivot_var->is_redundant &&
 973                                     pivot_var->index == row)
 974                                         isl_tab_pivot(tab, row, col);
 975                                 else
 976                                         restore_row(tab, var);
 977                         }
 978                         return 1;
 979                 }
 980         }
 981         if (var->is_redundant)
 982                 return 0;
 983         do {
 984                 find_pivot(tab, var, var, -1, &row, &col);
 985                 if (row == var->index)
 986                         return 1;
 987                 if (row == -1)
 988                         return 0;
 989                 pivot_var = var_from_col(tab, col);
 990                 isl_tab_pivot(tab, row, col);
 991                 if (var->is_redundant)
 992                         return 0;
 993         } while (!row_at_most_neg_one(tab, var->index));
 994         if (var->is_nonneg) {
 995                 /* pivot back to non-negative value */
 996                 if (!pivot_var->is_redundant && pivot_var->index == row)
 997                         isl_tab_pivot(tab, row, col);
 998                 restore_row(tab, var);
 999         }
1000         return 1;
1001 }
1002
1003 /* Return 1 if "var" can attain values >= 1.
1004  * Return 0 otherwise.
1005  */
1006 static int at_least_one(struct isl_tab *tab, struct isl_tab_var *var)
1007 {
1008         int row, col;
1009         isl_int *r;
1010
1011         if (max_is_manifestly_unbounded(tab, var))
1012                 return 1;
1013         to_row(tab, var, 1);
1014         r = tab->mat->row[var->index];
1015         while (isl_int_lt(r[1], r[0])) {
1016                 find_pivot(tab, var, var, 1, &row, &col);
1017                 if (row == -1)
1018                         return isl_int_ge(r[1], r[0]);
1019                 if (row == var->index) /* manifestly unbounded */
1020                         return 1;
1021                 isl_tab_pivot(tab, row, col);
1022         }
1023         return 1;
1024 }
1025
1026 static void swap_cols(struct isl_tab *tab, int col1, int col2)
1027 {
1028         int t;
1029         unsigned off = 2 + tab->M;
1030         t = tab->col_var[col1];
1031         tab->col_var[col1] = tab->col_var[col2];
1032         tab->col_var[col2] = t;
1033         var_from_col(tab, col1)->index = col1;
1034         var_from_col(tab, col2)->index = col2;
1035         tab->mat = isl_mat_swap_cols(tab->mat, off + col1, off + col2);
1036 }
1037
1038 /* Mark column with index "col" as representing a zero variable.
1039  * If we may need to undo the operation the column is kept,
1040  * but no longer considered.
1041  * Otherwise, the column is simply removed.
1042  *
1043  * The column may be interchanged with some other column.  If it
1044  * is interchanged with a later column, return 1.  Otherwise return 0.
1045  * If the columns are checked in order in the calling function,
1046  * then a return value of 1 means that the column with the given
1047  * column number may now contain a different column that
1048  * hasn't been checked yet.
1049  */
1050 int isl_tab_kill_col(struct isl_tab *tab, int col)
1051 {
1052         var_from_col(tab, col)->is_zero = 1;
1053         if (tab->need_undo) {
1054                 isl_tab_push_var(tab, isl_tab_undo_zero, var_from_col(tab, col));
1055                 if (col != tab->n_dead)
1056                         swap_cols(tab, col, tab->n_dead);
1057                 tab->n_dead++;
1058                 return 0;
1059         } else {
1060                 if (col != tab->n_col - 1)
1061                         swap_cols(tab, col, tab->n_col - 1);
1062                 var_from_col(tab, tab->n_col - 1)->index = -1;
1063                 tab->n_col--;
1064                 return 1;
1065         }
1066 }
1067
1068 /* Row variable "var" is non-negative and cannot attain any values
1069  * larger than zero.  This means that the coefficients of the unrestricted
1070  * column variables are zero and that the coefficients of the non-negative
1071  * column variables are zero or negative.
1072  * Each of the non-negative variables with a negative coefficient can
1073  * then also be written as the negative sum of non-negative variables
1074  * and must therefore also be zero.
1075  */
1076 static void close_row(struct isl_tab *tab, struct isl_tab_var *var)
1077 {
1078         int j;
1079         struct isl_mat *mat = tab->mat;
1080         unsigned off = 2 + tab->M;
1081
1082         isl_assert(tab->mat->ctx, var->is_nonneg, return);
1083         var->is_zero = 1;
1084         if (tab->need_undo)
1085                 isl_tab_push_var(tab, isl_tab_undo_zero, var);
1086         for (j = tab->n_dead; j < tab->n_col; ++j) {
1087                 if (isl_int_is_zero(mat->row[var->index][off + j]))
1088                         continue;
1089                 isl_assert(tab->mat->ctx,
1090                         isl_int_is_neg(mat->row[var->index][off + j]), return);
1091                 if (isl_tab_kill_col(tab, j))
1092                         --j;
1093         }
1094         isl_tab_mark_redundant(tab, var->index);
1095 }
1096
1097 /* Add a constraint to the tableau and allocate a row for it.
1098  * Return the index into the constraint array "con".
1099  */
1100 int isl_tab_allocate_con(struct isl_tab *tab)
1101 {
1102         int r;
1103
1104         isl_assert(tab->mat->ctx, tab->n_row < tab->mat->n_row, return -1);
1105         isl_assert(tab->mat->ctx, tab->n_con < tab->max_con, return -1);
1106
1107         r = tab->n_con;
1108         tab->con[r].index = tab->n_row;
1109         tab->con[r].is_row = 1;
1110         tab->con[r].is_nonneg = 0;
1111         tab->con[r].is_zero = 0;
1112         tab->con[r].is_redundant = 0;
1113         tab->con[r].frozen = 0;
1114         tab->con[r].negated = 0;
1115         tab->row_var[tab->n_row] = ~r;
1116
1117         tab->n_row++;
1118         tab->n_con++;
1119         isl_tab_push_var(tab, isl_tab_undo_allocate, &tab->con[r]);
1120
1121         return r;
1122 }
1123
1124 /* Add a variable to the tableau and allocate a column for it.
1125  * Return the index into the variable array "var".
1126  */
1127 int isl_tab_allocate_var(struct isl_tab *tab)
1128 {
1129         int r;
1130         int i;
1131         unsigned off = 2 + tab->M;
1132
1133         isl_assert(tab->mat->ctx, tab->n_col < tab->mat->n_col, return -1);
1134         isl_assert(tab->mat->ctx, tab->n_var < tab->max_var, return -1);
1135
1136         r = tab->n_var;
1137         tab->var[r].index = tab->n_col;
1138         tab->var[r].is_row = 0;
1139         tab->var[r].is_nonneg = 0;
1140         tab->var[r].is_zero = 0;
1141         tab->var[r].is_redundant = 0;
1142         tab->var[r].frozen = 0;
1143         tab->var[r].negated = 0;
1144         tab->col_var[tab->n_col] = r;
1145
1146         for (i = 0; i < tab->n_row; ++i)
1147                 isl_int_set_si(tab->mat->row[i][off + tab->n_col], 0);
1148
1149         tab->n_var++;
1150         tab->n_col++;
1151         isl_tab_push_var(tab, isl_tab_undo_allocate, &tab->var[r]);
1152
1153         return r;
1154 }
1155
1156 /* Add a row to the tableau.  The row is given as an affine combination
1157  * of the original variables and needs to be expressed in terms of the
1158  * column variables.
1159  *
1160  * We add each term in turn.
1161  * If r = n/d_r is the current sum and we need to add k x, then
1162  *      if x is a column variable, we increase the numerator of
1163  *              this column by k d_r
1164  *      if x = f/d_x is a row variable, then the new representation of r is
1165  *
1166  *               n    k f   d_x/g n + d_r/g k f   m/d_r n + m/d_g k f
1167  *              --- + --- = ------------------- = -------------------
1168  *              d_r   d_r        d_r d_x/g                m
1169  *
1170  *      with g the gcd of d_r and d_x and m the lcm of d_r and d_x.
1171  */
1172 int isl_tab_add_row(struct isl_tab *tab, isl_int *line)
1173 {
1174         int i;
1175         int r;
1176         isl_int *row;
1177         isl_int a, b;
1178         unsigned off = 2 + tab->M;
1179
1180         r = isl_tab_allocate_con(tab);
1181         if (r < 0)
1182                 return -1;
1183
1184         isl_int_init(a);
1185         isl_int_init(b);
1186         row = tab->mat->row[tab->con[r].index];
1187         isl_int_set_si(row[0], 1);
1188         isl_int_set(row[1], line[0]);
1189         isl_seq_clr(row + 2, tab->M + tab->n_col);
1190         for (i = 0; i < tab->n_var; ++i) {
1191                 if (tab->var[i].is_zero)
1192                         continue;
1193                 if (tab->var[i].is_row) {
1194                         isl_int_lcm(a,
1195                                 row[0], tab->mat->row[tab->var[i].index][0]);
1196                         isl_int_swap(a, row[0]);
1197                         isl_int_divexact(a, row[0], a);
1198                         isl_int_divexact(b,
1199                                 row[0], tab->mat->row[tab->var[i].index][0]);
1200                         isl_int_mul(b, b, line[1 + i]);
1201                         isl_seq_combine(row + 1, a, row + 1,
1202                             b, tab->mat->row[tab->var[i].index] + 1,
1203                             1 + tab->M + tab->n_col);
1204                 } else
1205                         isl_int_addmul(row[off + tab->var[i].index],
1206                                                         line[1 + i], row[0]);
1207                 if (tab->M && i >= tab->n_param && i < tab->n_var - tab->n_div)
1208                         isl_int_submul(row[2], line[1 + i], row[0]);
1209         }
1210         isl_seq_normalize(tab->mat->ctx, row, off + tab->n_col);
1211         isl_int_clear(a);
1212         isl_int_clear(b);
1213
1214         if (tab->row_sign)
1215                 tab->row_sign[tab->con[r].index] = 0;
1216
1217         return r;
1218 }
1219
1220 static int drop_row(struct isl_tab *tab, int row)
1221 {
1222         isl_assert(tab->mat->ctx, ~tab->row_var[row] == tab->n_con - 1, return -1);
1223         if (row != tab->n_row - 1)
1224                 swap_rows(tab, row, tab->n_row - 1);
1225         tab->n_row--;
1226         tab->n_con--;
1227         return 0;
1228 }
1229
1230 static int drop_col(struct isl_tab *tab, int col)
1231 {
1232         isl_assert(tab->mat->ctx, tab->col_var[col] == tab->n_var - 1, return -1);
1233         if (col != tab->n_col - 1)
1234                 swap_cols(tab, col, tab->n_col - 1);
1235         tab->n_col--;
1236         tab->n_var--;
1237         return 0;
1238 }
1239
1240 /* Add inequality "ineq" and check if it conflicts with the
1241  * previously added constraints or if it is obviously redundant.
1242  */
1243 struct isl_tab *isl_tab_add_ineq(struct isl_tab *tab, isl_int *ineq)
1244 {
1245         int r;
1246         int sgn;
1247
1248         if (!tab)
1249                 return NULL;
1250         r = isl_tab_add_row(tab, ineq);
1251         if (r < 0)
1252                 goto error;
1253         tab->con[r].is_nonneg = 1;
1254         isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]);
1255         if (isl_tab_row_is_redundant(tab, tab->con[r].index)) {
1256                 isl_tab_mark_redundant(tab, tab->con[r].index);
1257                 return tab;
1258         }
1259
1260         sgn = restore_row(tab, &tab->con[r]);
1261         if (sgn < 0)
1262                 return isl_tab_mark_empty(tab);
1263         if (tab->con[r].is_row && isl_tab_row_is_redundant(tab, tab->con[r].index))
1264                 isl_tab_mark_redundant(tab, tab->con[r].index);
1265         return tab;
1266 error:
1267         isl_tab_free(tab);
1268         return NULL;
1269 }
1270
1271 /* Pivot a non-negative variable down until it reaches the value zero
1272  * and then pivot the variable into a column position.
1273  */
1274 static int to_col(struct isl_tab *tab, struct isl_tab_var *var)
1275 {
1276         int i;
1277         int row, col;
1278         unsigned off = 2 + tab->M;
1279
1280         if (!var->is_row)
1281                 return 0;
1282
1283         while (isl_int_is_pos(tab->mat->row[var->index][1])) {
1284                 find_pivot(tab, var, NULL, -1, &row, &col);
1285                 isl_assert(tab->mat->ctx, row != -1, return -1);
1286                 isl_tab_pivot(tab, row, col);
1287                 if (!var->is_row)
1288                         return 0;
1289         }
1290
1291         for (i = tab->n_dead; i < tab->n_col; ++i)
1292                 if (!isl_int_is_zero(tab->mat->row[var->index][off + i]))
1293                         break;
1294
1295         isl_assert(tab->mat->ctx, i < tab->n_col, return -1);
1296         isl_tab_pivot(tab, var->index, i);
1297
1298         return 0;
1299 }
1300
1301 /* We assume Gaussian elimination has been performed on the equalities.
1302  * The equalities can therefore never conflict.
1303  * Adding the equalities is currently only really useful for a later call
1304  * to isl_tab_ineq_type.
1305  */
1306 static struct isl_tab *add_eq(struct isl_tab *tab, isl_int *eq)
1307 {
1308         int i;
1309         int r;
1310
1311         if (!tab)
1312                 return NULL;
1313         r = isl_tab_add_row(tab, eq);
1314         if (r < 0)
1315                 goto error;
1316
1317         r = tab->con[r].index;
1318         i = isl_seq_first_non_zero(tab->mat->row[r] + 2 + tab->M + tab->n_dead,
1319                                         tab->n_col - tab->n_dead);
1320         isl_assert(tab->mat->ctx, i >= 0, goto error);
1321         i += tab->n_dead;
1322         isl_tab_pivot(tab, r, i);
1323         isl_tab_kill_col(tab, i);
1324         tab->n_eq++;
1325
1326         return tab;
1327 error:
1328         isl_tab_free(tab);
1329         return NULL;
1330 }
1331
1332 static int row_is_manifestly_zero(struct isl_tab *tab, int row)
1333 {
1334         unsigned off = 2 + tab->M;
1335
1336         if (!isl_int_is_zero(tab->mat->row[row][1]))
1337                 return 0;
1338         if (tab->M && !isl_int_is_zero(tab->mat->row[row][2]))
1339                 return 0;
1340         return isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
1341                                         tab->n_col - tab->n_dead) == -1;
1342 }
1343
1344 /* Add an equality that is known to be valid for the given tableau.
1345  */
1346 struct isl_tab *isl_tab_add_valid_eq(struct isl_tab *tab, isl_int *eq)
1347 {
1348         struct isl_tab_var *var;
1349         int r;
1350
1351         if (!tab)
1352                 return NULL;
1353         r = isl_tab_add_row(tab, eq);
1354         if (r < 0)
1355                 goto error;
1356
1357         var = &tab->con[r];
1358         r = var->index;
1359         if (row_is_manifestly_zero(tab, r)) {
1360                 var->is_zero = 1;
1361                 isl_tab_mark_redundant(tab, r);
1362                 return tab;
1363         }
1364
1365         if (isl_int_is_neg(tab->mat->row[r][1])) {
1366                 isl_seq_neg(tab->mat->row[r] + 1, tab->mat->row[r] + 1,
1367                             1 + tab->n_col);
1368                 var->negated = 1;
1369         }
1370         var->is_nonneg = 1;
1371         if (to_col(tab, var) < 0)
1372                 goto error;
1373         var->is_nonneg = 0;
1374         isl_tab_kill_col(tab, var->index);
1375
1376         return tab;
1377 error:
1378         isl_tab_free(tab);
1379         return NULL;
1380 }
1381
1382 /* Add equality "eq" and check if it conflicts with the
1383  * previously added constraints or if it is obviously redundant.
1384  */
1385 struct isl_tab *isl_tab_add_eq(struct isl_tab *tab, isl_int *eq)
1386 {
1387         struct isl_tab_undo *snap = NULL;
1388         struct isl_tab_var *var;
1389         int r;
1390         int row;
1391         int sgn;
1392
1393         if (!tab)
1394                 return NULL;
1395         isl_assert(tab->mat->ctx, !tab->M, goto error);
1396
1397         if (tab->need_undo)
1398                 snap = isl_tab_snap(tab);
1399
1400         r = isl_tab_add_row(tab, eq);
1401         if (r < 0)
1402                 goto error;
1403
1404         var = &tab->con[r];
1405         row = var->index;
1406         if (row_is_manifestly_zero(tab, row)) {
1407                 if (snap) {
1408                         if (isl_tab_rollback(tab, snap) < 0)
1409                                 goto error;
1410                 } else
1411                         drop_row(tab, row);
1412                 return tab;
1413         }
1414
1415         sgn = isl_int_sgn(tab->mat->row[row][1]);
1416
1417         if (sgn > 0) {
1418                 isl_seq_neg(tab->mat->row[row] + 1, tab->mat->row[row] + 1,
1419                             1 + tab->n_col);
1420                 var->negated = 1;
1421                 sgn = -1;
1422         }
1423
1424         if (sgn < 0 && sign_of_max(tab, var) < 0)
1425                 return isl_tab_mark_empty(tab);
1426
1427         var->is_nonneg = 1;
1428         if (to_col(tab, var) < 0)
1429                 goto error;
1430         var->is_nonneg = 0;
1431         isl_tab_kill_col(tab, var->index);
1432
1433         return tab;
1434 error:
1435         isl_tab_free(tab);
1436         return NULL;
1437 }
1438
1439 struct isl_tab *isl_tab_from_basic_map(struct isl_basic_map *bmap)
1440 {
1441         int i;
1442         struct isl_tab *tab;
1443
1444         if (!bmap)
1445                 return NULL;
1446         tab = isl_tab_alloc(bmap->ctx,
1447                             isl_basic_map_total_dim(bmap) + bmap->n_ineq + 1,
1448                             isl_basic_map_total_dim(bmap), 0);
1449         if (!tab)
1450                 return NULL;
1451         tab->rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
1452         if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
1453                 return isl_tab_mark_empty(tab);
1454         for (i = 0; i < bmap->n_eq; ++i) {
1455                 tab = add_eq(tab, bmap->eq[i]);
1456                 if (!tab)
1457                         return tab;
1458         }
1459         for (i = 0; i < bmap->n_ineq; ++i) {
1460                 tab = isl_tab_add_ineq(tab, bmap->ineq[i]);
1461                 if (!tab || tab->empty)
1462                         return tab;
1463         }
1464         return tab;
1465 }
1466
1467 struct isl_tab *isl_tab_from_basic_set(struct isl_basic_set *bset)
1468 {
1469         return isl_tab_from_basic_map((struct isl_basic_map *)bset);
1470 }
1471
1472 /* Construct a tableau corresponding to the recession cone of "bset".
1473  */
1474 struct isl_tab *isl_tab_from_recession_cone(struct isl_basic_set *bset)
1475 {
1476         isl_int cst;
1477         int i;
1478         struct isl_tab *tab;
1479
1480         if (!bset)
1481                 return NULL;
1482         tab = isl_tab_alloc(bset->ctx, bset->n_eq + bset->n_ineq,
1483                                 isl_basic_set_total_dim(bset), 0);
1484         if (!tab)
1485                 return NULL;
1486         tab->rational = ISL_F_ISSET(bset, ISL_BASIC_SET_RATIONAL);
1487
1488         isl_int_init(cst);
1489         for (i = 0; i < bset->n_eq; ++i) {
1490                 isl_int_swap(bset->eq[i][0], cst);
1491                 tab = add_eq(tab, bset->eq[i]);
1492                 isl_int_swap(bset->eq[i][0], cst);
1493                 if (!tab)
1494                         goto done;
1495         }
1496         for (i = 0; i < bset->n_ineq; ++i) {
1497                 int r;
1498                 isl_int_swap(bset->ineq[i][0], cst);
1499                 r = isl_tab_add_row(tab, bset->ineq[i]);
1500                 isl_int_swap(bset->ineq[i][0], cst);
1501                 if (r < 0)
1502                         goto error;
1503                 tab->con[r].is_nonneg = 1;
1504                 isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]);
1505         }
1506 done:
1507         isl_int_clear(cst);
1508         return tab;
1509 error:
1510         isl_int_clear(cst);
1511         isl_tab_free(tab);
1512         return NULL;
1513 }
1514
1515 /* Assuming "tab" is the tableau of a cone, check if the cone is
1516  * bounded, i.e., if it is empty or only contains the origin.
1517  */
1518 int isl_tab_cone_is_bounded(struct isl_tab *tab)
1519 {
1520         int i;
1521
1522         if (!tab)
1523                 return -1;
1524         if (tab->empty)
1525                 return 1;
1526         if (tab->n_dead == tab->n_col)
1527                 return 1;
1528
1529         for (;;) {
1530                 for (i = tab->n_redundant; i < tab->n_row; ++i) {
1531                         struct isl_tab_var *var;
1532                         var = isl_tab_var_from_row(tab, i);
1533                         if (!var->is_nonneg)
1534                                 continue;
1535                         if (sign_of_max(tab, var) != 0)
1536                                 return 0;
1537                         close_row(tab, var);
1538                         break;
1539                 }
1540                 if (tab->n_dead == tab->n_col)
1541                         return 1;
1542                 if (i == tab->n_row)
1543                         return 0;
1544         }
1545 }
1546
1547 int isl_tab_sample_is_integer(struct isl_tab *tab)
1548 {
1549         int i;
1550
1551         if (!tab)
1552                 return -1;
1553
1554         for (i = 0; i < tab->n_var; ++i) {
1555                 int row;
1556                 if (!tab->var[i].is_row)
1557                         continue;
1558                 row = tab->var[i].index;
1559                 if (!isl_int_is_divisible_by(tab->mat->row[row][1],
1560                                                 tab->mat->row[row][0]))
1561                         return 0;
1562         }
1563         return 1;
1564 }
1565
1566 static struct isl_vec *extract_integer_sample(struct isl_tab *tab)
1567 {
1568         int i;
1569         struct isl_vec *vec;
1570
1571         vec = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
1572         if (!vec)
1573                 return NULL;
1574
1575         isl_int_set_si(vec->block.data[0], 1);
1576         for (i = 0; i < tab->n_var; ++i) {
1577                 if (!tab->var[i].is_row)
1578                         isl_int_set_si(vec->block.data[1 + i], 0);
1579                 else {
1580                         int row = tab->var[i].index;
1581                         isl_int_divexact(vec->block.data[1 + i],
1582                                 tab->mat->row[row][1], tab->mat->row[row][0]);
1583                 }
1584         }
1585
1586         return vec;
1587 }
1588
1589 struct isl_vec *isl_tab_get_sample_value(struct isl_tab *tab)
1590 {
1591         int i;
1592         struct isl_vec *vec;
1593         isl_int m;
1594
1595         if (!tab)
1596                 return NULL;
1597
1598         vec = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
1599         if (!vec)
1600                 return NULL;
1601
1602         isl_int_init(m);
1603
1604         isl_int_set_si(vec->block.data[0], 1);
1605         for (i = 0; i < tab->n_var; ++i) {
1606                 int row;
1607                 if (!tab->var[i].is_row) {
1608                         isl_int_set_si(vec->block.data[1 + i], 0);
1609                         continue;
1610                 }
1611                 row = tab->var[i].index;
1612                 isl_int_gcd(m, vec->block.data[0], tab->mat->row[row][0]);
1613                 isl_int_divexact(m, tab->mat->row[row][0], m);
1614                 isl_seq_scale(vec->block.data, vec->block.data, m, 1 + i);
1615                 isl_int_divexact(m, vec->block.data[0], tab->mat->row[row][0]);
1616                 isl_int_mul(vec->block.data[1 + i], m, tab->mat->row[row][1]);
1617         }
1618         vec = isl_vec_normalize(vec);
1619
1620         isl_int_clear(m);
1621         return vec;
1622 }
1623
1624 /* Update "bmap" based on the results of the tableau "tab".
1625  * In particular, implicit equalities are made explicit, redundant constraints
1626  * are removed and if the sample value happens to be integer, it is stored
1627  * in "bmap" (unless "bmap" already had an integer sample).
1628  *
1629  * The tableau is assumed to have been created from "bmap" using
1630  * isl_tab_from_basic_map.
1631  */
1632 struct isl_basic_map *isl_basic_map_update_from_tab(struct isl_basic_map *bmap,
1633         struct isl_tab *tab)
1634 {
1635         int i;
1636         unsigned n_eq;
1637
1638         if (!bmap)
1639                 return NULL;
1640         if (!tab)
1641                 return bmap;
1642
1643         n_eq = tab->n_eq;
1644         if (tab->empty)
1645                 bmap = isl_basic_map_set_to_empty(bmap);
1646         else
1647                 for (i = bmap->n_ineq - 1; i >= 0; --i) {
1648                         if (isl_tab_is_equality(tab, n_eq + i))
1649                                 isl_basic_map_inequality_to_equality(bmap, i);
1650                         else if (isl_tab_is_redundant(tab, n_eq + i))
1651                                 isl_basic_map_drop_inequality(bmap, i);
1652                 }
1653         if (!tab->rational &&
1654             !bmap->sample && isl_tab_sample_is_integer(tab))
1655                 bmap->sample = extract_integer_sample(tab);
1656         return bmap;
1657 }
1658
1659 struct isl_basic_set *isl_basic_set_update_from_tab(struct isl_basic_set *bset,
1660         struct isl_tab *tab)
1661 {
1662         return (struct isl_basic_set *)isl_basic_map_update_from_tab(
1663                 (struct isl_basic_map *)bset, tab);
1664 }
1665
1666 /* Given a non-negative variable "var", add a new non-negative variable
1667  * that is the opposite of "var", ensuring that var can only attain the
1668  * value zero.
1669  * If var = n/d is a row variable, then the new variable = -n/d.
1670  * If var is a column variables, then the new variable = -var.
1671  * If the new variable cannot attain non-negative values, then
1672  * the resulting tableau is empty.
1673  * Otherwise, we know the value will be zero and we close the row.
1674  */
1675 static struct isl_tab *cut_to_hyperplane(struct isl_tab *tab,
1676         struct isl_tab_var *var)
1677 {
1678         unsigned r;
1679         isl_int *row;
1680         int sgn;
1681         unsigned off = 2 + tab->M;
1682
1683         if (var->is_zero)
1684                 return tab;
1685         isl_assert(tab->mat->ctx, !var->is_redundant, goto error);
1686
1687         if (isl_tab_extend_cons(tab, 1) < 0)
1688                 goto error;
1689
1690         r = tab->n_con;
1691         tab->con[r].index = tab->n_row;
1692         tab->con[r].is_row = 1;
1693         tab->con[r].is_nonneg = 0;
1694         tab->con[r].is_zero = 0;
1695         tab->con[r].is_redundant = 0;
1696         tab->con[r].frozen = 0;
1697         tab->con[r].negated = 0;
1698         tab->row_var[tab->n_row] = ~r;
1699         row = tab->mat->row[tab->n_row];
1700
1701         if (var->is_row) {
1702                 isl_int_set(row[0], tab->mat->row[var->index][0]);
1703                 isl_seq_neg(row + 1,
1704                             tab->mat->row[var->index] + 1, 1 + tab->n_col);
1705         } else {
1706                 isl_int_set_si(row[0], 1);
1707                 isl_seq_clr(row + 1, 1 + tab->n_col);
1708                 isl_int_set_si(row[off + var->index], -1);
1709         }
1710
1711         tab->n_row++;
1712         tab->n_con++;
1713         isl_tab_push_var(tab, isl_tab_undo_allocate, &tab->con[r]);
1714
1715         sgn = sign_of_max(tab, &tab->con[r]);
1716         if (sgn < 0)
1717                 return isl_tab_mark_empty(tab);
1718         tab->con[r].is_nonneg = 1;
1719         isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]);
1720         /* sgn == 0 */
1721         close_row(tab, &tab->con[r]);
1722
1723         return tab;
1724 error:
1725         isl_tab_free(tab);
1726         return NULL;
1727 }
1728
1729 /* Given a tableau "tab" and an inequality constraint "con" of the tableau,
1730  * relax the inequality by one.  That is, the inequality r >= 0 is replaced
1731  * by r' = r + 1 >= 0.
1732  * If r is a row variable, we simply increase the constant term by one
1733  * (taking into account the denominator).
1734  * If r is a column variable, then we need to modify each row that
1735  * refers to r = r' - 1 by substituting this equality, effectively
1736  * subtracting the coefficient of the column from the constant.
1737  */
1738 struct isl_tab *isl_tab_relax(struct isl_tab *tab, int con)
1739 {
1740         struct isl_tab_var *var;
1741         unsigned off = 2 + tab->M;
1742
1743         if (!tab)
1744                 return NULL;
1745
1746         var = &tab->con[con];
1747
1748         if (!var->is_row && !max_is_manifestly_unbounded(tab, var))
1749                 to_row(tab, var, 1);
1750
1751         if (var->is_row)
1752                 isl_int_add(tab->mat->row[var->index][1],
1753                     tab->mat->row[var->index][1], tab->mat->row[var->index][0]);
1754         else {
1755                 int i;
1756
1757                 for (i = 0; i < tab->n_row; ++i) {
1758                         if (isl_int_is_zero(tab->mat->row[i][off + var->index]))
1759                                 continue;
1760                         isl_int_sub(tab->mat->row[i][1], tab->mat->row[i][1],
1761                             tab->mat->row[i][off + var->index]);
1762                 }
1763
1764         }
1765
1766         isl_tab_push_var(tab, isl_tab_undo_relax, var);
1767
1768         return tab;
1769 }
1770
1771 struct isl_tab *isl_tab_select_facet(struct isl_tab *tab, int con)
1772 {
1773         if (!tab)
1774                 return NULL;
1775
1776         return cut_to_hyperplane(tab, &tab->con[con]);
1777 }
1778
1779 static int may_be_equality(struct isl_tab *tab, int row)
1780 {
1781         unsigned off = 2 + tab->M;
1782         return (tab->rational ? isl_int_is_zero(tab->mat->row[row][1])
1783                               : isl_int_lt(tab->mat->row[row][1],
1784                                             tab->mat->row[row][0])) &&
1785                 isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
1786                                         tab->n_col - tab->n_dead) != -1;
1787 }
1788
1789 /* Check for (near) equalities among the constraints.
1790  * A constraint is an equality if it is non-negative and if
1791  * its maximal value is either
1792  *      - zero (in case of rational tableaus), or
1793  *      - strictly less than 1 (in case of integer tableaus)
1794  *
1795  * We first mark all non-redundant and non-dead variables that
1796  * are not frozen and not obviously not an equality.
1797  * Then we iterate over all marked variables if they can attain
1798  * any values larger than zero or at least one.
1799  * If the maximal value is zero, we mark any column variables
1800  * that appear in the row as being zero and mark the row as being redundant.
1801  * Otherwise, if the maximal value is strictly less than one (and the
1802  * tableau is integer), then we restrict the value to being zero
1803  * by adding an opposite non-negative variable.
1804  */
1805 struct isl_tab *isl_tab_detect_implicit_equalities(struct isl_tab *tab)
1806 {
1807         int i;
1808         unsigned n_marked;
1809
1810         if (!tab)
1811                 return NULL;
1812         if (tab->empty)
1813                 return tab;
1814         if (tab->n_dead == tab->n_col)
1815                 return tab;
1816
1817         n_marked = 0;
1818         for (i = tab->n_redundant; i < tab->n_row; ++i) {
1819                 struct isl_tab_var *var = isl_tab_var_from_row(tab, i);
1820                 var->marked = !var->frozen && var->is_nonneg &&
1821                         may_be_equality(tab, i);
1822                 if (var->marked)
1823                         n_marked++;
1824         }
1825         for (i = tab->n_dead; i < tab->n_col; ++i) {
1826                 struct isl_tab_var *var = var_from_col(tab, i);
1827                 var->marked = !var->frozen && var->is_nonneg;
1828                 if (var->marked)
1829                         n_marked++;
1830         }
1831         while (n_marked) {
1832                 struct isl_tab_var *var;
1833                 for (i = tab->n_redundant; i < tab->n_row; ++i) {
1834                         var = isl_tab_var_from_row(tab, i);
1835                         if (var->marked)
1836                                 break;
1837                 }
1838                 if (i == tab->n_row) {
1839                         for (i = tab->n_dead; i < tab->n_col; ++i) {
1840                                 var = var_from_col(tab, i);
1841                                 if (var->marked)
1842                                         break;
1843                         }
1844                         if (i == tab->n_col)
1845                                 break;
1846                 }
1847                 var->marked = 0;
1848                 n_marked--;
1849                 if (sign_of_max(tab, var) == 0)
1850                         close_row(tab, var);
1851                 else if (!tab->rational && !at_least_one(tab, var)) {
1852                         tab = cut_to_hyperplane(tab, var);
1853                         return isl_tab_detect_implicit_equalities(tab);
1854                 }
1855                 for (i = tab->n_redundant; i < tab->n_row; ++i) {
1856                         var = isl_tab_var_from_row(tab, i);
1857                         if (!var->marked)
1858                                 continue;
1859                         if (may_be_equality(tab, i))
1860                                 continue;
1861                         var->marked = 0;
1862                         n_marked--;
1863                 }
1864         }
1865
1866         return tab;
1867 }
1868
1869 /* Check for (near) redundant constraints.
1870  * A constraint is redundant if it is non-negative and if
1871  * its minimal value (temporarily ignoring the non-negativity) is either
1872  *      - zero (in case of rational tableaus), or
1873  *      - strictly larger than -1 (in case of integer tableaus)
1874  *
1875  * We first mark all non-redundant and non-dead variables that
1876  * are not frozen and not obviously negatively unbounded.
1877  * Then we iterate over all marked variables if they can attain
1878  * any values smaller than zero or at most negative one.
1879  * If not, we mark the row as being redundant (assuming it hasn't
1880  * been detected as being obviously redundant in the mean time).
1881  */
1882 struct isl_tab *isl_tab_detect_redundant(struct isl_tab *tab)
1883 {
1884         int i;
1885         unsigned n_marked;
1886
1887         if (!tab)
1888                 return NULL;
1889         if (tab->empty)
1890                 return tab;
1891         if (tab->n_redundant == tab->n_row)
1892                 return tab;
1893
1894         n_marked = 0;
1895         for (i = tab->n_redundant; i < tab->n_row; ++i) {
1896                 struct isl_tab_var *var = isl_tab_var_from_row(tab, i);
1897                 var->marked = !var->frozen && var->is_nonneg;
1898                 if (var->marked)
1899                         n_marked++;
1900         }
1901         for (i = tab->n_dead; i < tab->n_col; ++i) {
1902                 struct isl_tab_var *var = var_from_col(tab, i);
1903                 var->marked = !var->frozen && var->is_nonneg &&
1904                         !min_is_manifestly_unbounded(tab, var);
1905                 if (var->marked)
1906                         n_marked++;
1907         }
1908         while (n_marked) {
1909                 struct isl_tab_var *var;
1910                 for (i = tab->n_redundant; i < tab->n_row; ++i) {
1911                         var = isl_tab_var_from_row(tab, i);
1912                         if (var->marked)
1913                                 break;
1914                 }
1915                 if (i == tab->n_row) {
1916                         for (i = tab->n_dead; i < tab->n_col; ++i) {
1917                                 var = var_from_col(tab, i);
1918                                 if (var->marked)
1919                                         break;
1920                         }
1921                         if (i == tab->n_col)
1922                                 break;
1923                 }
1924                 var->marked = 0;
1925                 n_marked--;
1926                 if ((tab->rational ? (sign_of_min(tab, var) >= 0)
1927                                    : !isl_tab_min_at_most_neg_one(tab, var)) &&
1928                     !var->is_redundant)
1929                         isl_tab_mark_redundant(tab, var->index);
1930                 for (i = tab->n_dead; i < tab->n_col; ++i) {
1931                         var = var_from_col(tab, i);
1932                         if (!var->marked)
1933                                 continue;
1934                         if (!min_is_manifestly_unbounded(tab, var))
1935                                 continue;
1936                         var->marked = 0;
1937                         n_marked--;
1938                 }
1939         }
1940
1941         return tab;
1942 }
1943
1944 int isl_tab_is_equality(struct isl_tab *tab, int con)
1945 {
1946         int row;
1947         unsigned off;
1948
1949         if (!tab)
1950                 return -1;
1951         if (tab->con[con].is_zero)
1952                 return 1;
1953         if (tab->con[con].is_redundant)
1954                 return 0;
1955         if (!tab->con[con].is_row)
1956                 return tab->con[con].index < tab->n_dead;
1957
1958         row = tab->con[con].index;
1959
1960         off = 2 + tab->M;
1961         return isl_int_is_zero(tab->mat->row[row][1]) &&
1962                 isl_seq_first_non_zero(tab->mat->row[row] + 2 + tab->n_dead,
1963                                         tab->n_col - tab->n_dead) == -1;
1964 }
1965
1966 /* Return the minimial value of the affine expression "f" with denominator
1967  * "denom" in *opt, *opt_denom, assuming the tableau is not empty and
1968  * the expression cannot attain arbitrarily small values.
1969  * If opt_denom is NULL, then *opt is rounded up to the nearest integer.
1970  * The return value reflects the nature of the result (empty, unbounded,
1971  * minmimal value returned in *opt).
1972  */
1973 enum isl_lp_result isl_tab_min(struct isl_tab *tab,
1974         isl_int *f, isl_int denom, isl_int *opt, isl_int *opt_denom,
1975         unsigned flags)
1976 {
1977         int r;
1978         enum isl_lp_result res = isl_lp_ok;
1979         struct isl_tab_var *var;
1980         struct isl_tab_undo *snap;
1981
1982         if (tab->empty)
1983                 return isl_lp_empty;
1984
1985         snap = isl_tab_snap(tab);
1986         r = isl_tab_add_row(tab, f);
1987         if (r < 0)
1988                 return isl_lp_error;
1989         var = &tab->con[r];
1990         isl_int_mul(tab->mat->row[var->index][0],
1991                     tab->mat->row[var->index][0], denom);
1992         for (;;) {
1993                 int row, col;
1994                 find_pivot(tab, var, var, -1, &row, &col);
1995                 if (row == var->index) {
1996                         res = isl_lp_unbounded;
1997                         break;
1998                 }
1999                 if (row == -1)
2000                         break;
2001                 isl_tab_pivot(tab, row, col);
2002         }
2003         if (ISL_FL_ISSET(flags, ISL_TAB_SAVE_DUAL)) {
2004                 int i;
2005
2006                 isl_vec_free(tab->dual);
2007                 tab->dual = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_con);
2008                 if (!tab->dual)
2009                         return isl_lp_error;
2010                 isl_int_set(tab->dual->el[0], tab->mat->row[var->index][0]);
2011                 for (i = 0; i < tab->n_con; ++i) {
2012                         int pos;
2013                         if (tab->con[i].is_row) {
2014                                 isl_int_set_si(tab->dual->el[1 + i], 0);
2015                                 continue;
2016                         }
2017                         pos = 2 + tab->M + tab->con[i].index;
2018                         if (tab->con[i].negated)
2019                                 isl_int_neg(tab->dual->el[1 + i],
2020                                             tab->mat->row[var->index][pos]);
2021                         else
2022                                 isl_int_set(tab->dual->el[1 + i],
2023                                             tab->mat->row[var->index][pos]);
2024                 }
2025         }
2026         if (opt && res == isl_lp_ok) {
2027                 if (opt_denom) {
2028                         isl_int_set(*opt, tab->mat->row[var->index][1]);
2029                         isl_int_set(*opt_denom, tab->mat->row[var->index][0]);
2030                 } else
2031                         isl_int_cdiv_q(*opt, tab->mat->row[var->index][1],
2032                                              tab->mat->row[var->index][0]);
2033         }
2034         if (isl_tab_rollback(tab, snap) < 0)
2035                 return isl_lp_error;
2036         return res;
2037 }
2038
2039 int isl_tab_is_redundant(struct isl_tab *tab, int con)
2040 {
2041         if (!tab)
2042                 return -1;
2043         if (tab->con[con].is_zero)
2044                 return 0;
2045         if (tab->con[con].is_redundant)
2046                 return 1;
2047         return tab->con[con].is_row && tab->con[con].index < tab->n_redundant;
2048 }
2049
2050 /* Take a snapshot of the tableau that can be restored by s call to
2051  * isl_tab_rollback.
2052  */
2053 struct isl_tab_undo *isl_tab_snap(struct isl_tab *tab)
2054 {
2055         if (!tab)
2056                 return NULL;
2057         tab->need_undo = 1;
2058         return tab->top;
2059 }
2060
2061 /* Undo the operation performed by isl_tab_relax.
2062  */
2063 static void unrelax(struct isl_tab *tab, struct isl_tab_var *var)
2064 {
2065         unsigned off = 2 + tab->M;
2066
2067         if (!var->is_row && !max_is_manifestly_unbounded(tab, var))
2068                 to_row(tab, var, 1);
2069
2070         if (var->is_row)
2071                 isl_int_sub(tab->mat->row[var->index][1],
2072                     tab->mat->row[var->index][1], tab->mat->row[var->index][0]);
2073         else {
2074                 int i;
2075
2076                 for (i = 0; i < tab->n_row; ++i) {
2077                         if (isl_int_is_zero(tab->mat->row[i][off + var->index]))
2078                                 continue;
2079                         isl_int_add(tab->mat->row[i][1], tab->mat->row[i][1],
2080                             tab->mat->row[i][off + var->index]);
2081                 }
2082
2083         }
2084 }
2085
2086 static void perform_undo_var(struct isl_tab *tab, struct isl_tab_undo *undo)
2087 {
2088         struct isl_tab_var *var = var_from_index(tab, undo->u.var_index);
2089         switch(undo->type) {
2090         case isl_tab_undo_nonneg:
2091                 var->is_nonneg = 0;
2092                 break;
2093         case isl_tab_undo_redundant:
2094                 var->is_redundant = 0;
2095                 tab->n_redundant--;
2096                 break;
2097         case isl_tab_undo_zero:
2098                 var->is_zero = 0;
2099                 if (!var->is_row)
2100                         tab->n_dead--;
2101                 break;
2102         case isl_tab_undo_allocate:
2103                 if (undo->u.var_index >= 0) {
2104                         isl_assert(tab->mat->ctx, !var->is_row, return);
2105                         drop_col(tab, var->index);
2106                         break;
2107                 }
2108                 if (!var->is_row) {
2109                         if (!max_is_manifestly_unbounded(tab, var))
2110                                 to_row(tab, var, 1);
2111                         else if (!min_is_manifestly_unbounded(tab, var))
2112                                 to_row(tab, var, -1);
2113                         else
2114                                 to_row(tab, var, 0);
2115                 }
2116                 drop_row(tab, var->index);
2117                 break;
2118         case isl_tab_undo_relax:
2119                 unrelax(tab, var);
2120                 break;
2121         }
2122 }
2123
2124 /* Restore the tableau to the state where the basic variables
2125  * are those in "col_var".
2126  * We first construct a list of variables that are currently in
2127  * the basis, but shouldn't.  Then we iterate over all variables
2128  * that should be in the basis and for each one that is currently
2129  * not in the basis, we exchange it with one of the elements of the
2130  * list constructed before.
2131  * We can always find an appropriate variable to pivot with because
2132  * the current basis is mapped to the old basis by a non-singular
2133  * matrix and so we can never end up with a zero row.
2134  */
2135 static int restore_basis(struct isl_tab *tab, int *col_var)
2136 {
2137         int i, j;
2138         int n_extra = 0;
2139         int *extra = NULL;      /* current columns that contain bad stuff */
2140         unsigned off = 2 + tab->M;
2141
2142         extra = isl_alloc_array(tab->mat->ctx, int, tab->n_col);
2143         if (!extra)
2144                 goto error;
2145         for (i = 0; i < tab->n_col; ++i) {
2146                 for (j = 0; j < tab->n_col; ++j)
2147                         if (tab->col_var[i] == col_var[j])
2148                                 break;
2149                 if (j < tab->n_col)
2150                         continue;
2151                 extra[n_extra++] = i;
2152         }
2153         for (i = 0; i < tab->n_col && n_extra > 0; ++i) {
2154                 struct isl_tab_var *var;
2155                 int row;
2156
2157                 for (j = 0; j < tab->n_col; ++j)
2158                         if (col_var[i] == tab->col_var[j])
2159                                 break;
2160                 if (j < tab->n_col)
2161                         continue;
2162                 var = var_from_index(tab, col_var[i]);
2163                 row = var->index;
2164                 for (j = 0; j < n_extra; ++j)
2165                         if (!isl_int_is_zero(tab->mat->row[row][off+extra[j]]))
2166                                 break;
2167                 isl_assert(tab->mat->ctx, j < n_extra, goto error);
2168                 isl_tab_pivot(tab, row, extra[j]);
2169                 extra[j] = extra[--n_extra];
2170         }
2171
2172         free(extra);
2173         free(col_var);
2174         return 0;
2175 error:
2176         free(extra);
2177         free(col_var);
2178         return -1;
2179 }
2180
2181 static int perform_undo(struct isl_tab *tab, struct isl_tab_undo *undo)
2182 {
2183         switch (undo->type) {
2184         case isl_tab_undo_empty:
2185                 tab->empty = 0;
2186                 break;
2187         case isl_tab_undo_nonneg:
2188         case isl_tab_undo_redundant:
2189         case isl_tab_undo_zero:
2190         case isl_tab_undo_allocate:
2191         case isl_tab_undo_relax:
2192                 perform_undo_var(tab, undo);
2193                 break;
2194         case isl_tab_undo_bset_eq:
2195                 isl_basic_set_free_equality(tab->bset, 1);
2196                 break;
2197         case isl_tab_undo_bset_ineq:
2198                 isl_basic_set_free_inequality(tab->bset, 1);
2199                 break;
2200         case isl_tab_undo_bset_div:
2201                 isl_basic_set_free_div(tab->bset, 1);
2202                 if (tab->samples)
2203                         tab->samples->n_col--;
2204                 break;
2205         case isl_tab_undo_saved_basis:
2206                 if (restore_basis(tab, undo->u.col_var) < 0)
2207                         return -1;
2208                 break;
2209         case isl_tab_undo_drop_sample:
2210                 tab->n_outside--;
2211                 break;
2212         default:
2213                 isl_assert(tab->mat->ctx, 0, return -1);
2214         }
2215         return 0;
2216 }
2217
2218 /* Return the tableau to the state it was in when the snapshot "snap"
2219  * was taken.
2220  */
2221 int isl_tab_rollback(struct isl_tab *tab, struct isl_tab_undo *snap)
2222 {
2223         struct isl_tab_undo *undo, *next;
2224
2225         if (!tab)
2226                 return -1;
2227
2228         tab->in_undo = 1;
2229         for (undo = tab->top; undo && undo != &tab->bottom; undo = next) {
2230                 next = undo->next;
2231                 if (undo == snap)
2232                         break;
2233                 if (perform_undo(tab, undo) < 0) {
2234                         free_undo(tab);
2235                         tab->in_undo = 0;
2236                         return -1;
2237                 }
2238                 free(undo);
2239         }
2240         tab->in_undo = 0;
2241         tab->top = undo;
2242         if (!undo)
2243                 return -1;
2244         return 0;
2245 }
2246
2247 /* The given row "row" represents an inequality violated by all
2248  * points in the tableau.  Check for some special cases of such
2249  * separating constraints.
2250  * In particular, if the row has been reduced to the constant -1,
2251  * then we know the inequality is adjacent (but opposite) to
2252  * an equality in the tableau.
2253  * If the row has been reduced to r = -1 -r', with r' an inequality
2254  * of the tableau, then the inequality is adjacent (but opposite)
2255  * to the inequality r'.
2256  */
2257 static enum isl_ineq_type separation_type(struct isl_tab *tab, unsigned row)
2258 {
2259         int pos;
2260         unsigned off = 2 + tab->M;
2261
2262         if (tab->rational)
2263                 return isl_ineq_separate;
2264
2265         if (!isl_int_is_one(tab->mat->row[row][0]))
2266                 return isl_ineq_separate;
2267         if (!isl_int_is_negone(tab->mat->row[row][1]))
2268                 return isl_ineq_separate;
2269
2270         pos = isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
2271                                         tab->n_col - tab->n_dead);
2272         if (pos == -1)
2273                 return isl_ineq_adj_eq;
2274
2275         if (!isl_int_is_negone(tab->mat->row[row][off + tab->n_dead + pos]))
2276                 return isl_ineq_separate;
2277
2278         pos = isl_seq_first_non_zero(
2279                         tab->mat->row[row] + off + tab->n_dead + pos + 1,
2280                         tab->n_col - tab->n_dead - pos - 1);
2281
2282         return pos == -1 ? isl_ineq_adj_ineq : isl_ineq_separate;
2283 }
2284
2285 /* Check the effect of inequality "ineq" on the tableau "tab".
2286  * The result may be
2287  *      isl_ineq_redundant:     satisfied by all points in the tableau
2288  *      isl_ineq_separate:      satisfied by no point in the tableau
2289  *      isl_ineq_cut:           satisfied by some by not all points
2290  *      isl_ineq_adj_eq:        adjacent to an equality
2291  *      isl_ineq_adj_ineq:      adjacent to an inequality.
2292  */
2293 enum isl_ineq_type isl_tab_ineq_type(struct isl_tab *tab, isl_int *ineq)
2294 {
2295         enum isl_ineq_type type = isl_ineq_error;
2296         struct isl_tab_undo *snap = NULL;
2297         int con;
2298         int row;
2299
2300         if (!tab)
2301                 return isl_ineq_error;
2302
2303         if (isl_tab_extend_cons(tab, 1) < 0)
2304                 return isl_ineq_error;
2305
2306         snap = isl_tab_snap(tab);
2307
2308         con = isl_tab_add_row(tab, ineq);
2309         if (con < 0)
2310                 goto error;
2311
2312         row = tab->con[con].index;
2313         if (isl_tab_row_is_redundant(tab, row))
2314                 type = isl_ineq_redundant;
2315         else if (isl_int_is_neg(tab->mat->row[row][1]) &&
2316                  (tab->rational ||
2317                     isl_int_abs_ge(tab->mat->row[row][1],
2318                                    tab->mat->row[row][0]))) {
2319                 if (at_least_zero(tab, &tab->con[con]))
2320                         type = isl_ineq_cut;
2321                 else
2322                         type = separation_type(tab, row);
2323         } else if (tab->rational ? (sign_of_min(tab, &tab->con[con]) < 0)
2324                              : isl_tab_min_at_most_neg_one(tab, &tab->con[con]))
2325                 type = isl_ineq_cut;
2326         else
2327                 type = isl_ineq_redundant;
2328
2329         if (isl_tab_rollback(tab, snap))
2330                 return isl_ineq_error;
2331         return type;
2332 error:
2333         isl_tab_rollback(tab, snap);
2334         return isl_ineq_error;
2335 }
2336
2337 void isl_tab_dump(struct isl_tab *tab, FILE *out, int indent)
2338 {
2339         unsigned r, c;
2340         int i;
2341
2342         if (!tab) {
2343                 fprintf(out, "%*snull tab\n", indent, "");
2344                 return;
2345         }
2346         fprintf(out, "%*sn_redundant: %d, n_dead: %d", indent, "",
2347                 tab->n_redundant, tab->n_dead);
2348         if (tab->rational)
2349                 fprintf(out, ", rational");
2350         if (tab->empty)
2351                 fprintf(out, ", empty");
2352         fprintf(out, "\n");
2353         fprintf(out, "%*s[", indent, "");
2354         for (i = 0; i < tab->n_var; ++i) {
2355                 if (i)
2356                         fprintf(out, (i == tab->n_param ||
2357                                       i == tab->n_var - tab->n_div) ? "; "
2358                                                                     : ", ");
2359                 fprintf(out, "%c%d%s", tab->var[i].is_row ? 'r' : 'c',
2360                                         tab->var[i].index,
2361                                         tab->var[i].is_zero ? " [=0]" :
2362                                         tab->var[i].is_redundant ? " [R]" : "");
2363         }
2364         fprintf(out, "]\n");
2365         fprintf(out, "%*s[", indent, "");
2366         for (i = 0; i < tab->n_con; ++i) {
2367                 if (i)
2368                         fprintf(out, ", ");
2369                 fprintf(out, "%c%d%s", tab->con[i].is_row ? 'r' : 'c',
2370                                         tab->con[i].index,
2371                                         tab->con[i].is_zero ? " [=0]" :
2372                                         tab->con[i].is_redundant ? " [R]" : "");
2373         }
2374         fprintf(out, "]\n");
2375         fprintf(out, "%*s[", indent, "");
2376         for (i = 0; i < tab->n_row; ++i) {
2377                 const char *sign = "";
2378                 if (i)
2379                         fprintf(out, ", ");
2380                 if (tab->row_sign) {
2381                         if (tab->row_sign[i] == isl_tab_row_unknown)
2382                                 sign = "?";
2383                         else if (tab->row_sign[i] == isl_tab_row_neg)
2384                                 sign = "-";
2385                         else if (tab->row_sign[i] == isl_tab_row_pos)
2386                                 sign = "+";
2387                         else
2388                                 sign = "+-";
2389                 }
2390                 fprintf(out, "r%d: %d%s%s", i, tab->row_var[i],
2391                     isl_tab_var_from_row(tab, i)->is_nonneg ? " [>=0]" : "", sign);
2392         }
2393         fprintf(out, "]\n");
2394         fprintf(out, "%*s[", indent, "");
2395         for (i = 0; i < tab->n_col; ++i) {
2396                 if (i)
2397                         fprintf(out, ", ");
2398                 fprintf(out, "c%d: %d%s", i, tab->col_var[i],
2399                     var_from_col(tab, i)->is_nonneg ? " [>=0]" : "");
2400         }
2401         fprintf(out, "]\n");
2402         r = tab->mat->n_row;
2403         tab->mat->n_row = tab->n_row;
2404         c = tab->mat->n_col;
2405         tab->mat->n_col = 2 + tab->M + tab->n_col;
2406         isl_mat_dump(tab->mat, out, indent);
2407         tab->mat->n_row = r;
2408         tab->mat->n_col = c;
2409         if (tab->bset)
2410                 isl_basic_set_dump(tab->bset, out, indent);
2411 }