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 /* Add an equality that is known to be valid for the given tableau.
1333  */
1334 struct isl_tab *isl_tab_add_valid_eq(struct isl_tab *tab, isl_int *eq)
1335 {
1336         struct isl_tab_var *var;
1337         int r;
1338
1339         if (!tab)
1340                 return NULL;
1341         r = isl_tab_add_row(tab, eq);
1342         if (r < 0)
1343                 goto error;
1344
1345         var = &tab->con[r];
1346         r = var->index;
1347         if (isl_int_is_neg(tab->mat->row[r][1])) {
1348                 isl_seq_neg(tab->mat->row[r] + 1, tab->mat->row[r] + 1,
1349                             1 + tab->n_col);
1350                 var->negated = 1;
1351         }
1352         var->is_nonneg = 1;
1353         if (to_col(tab, var) < 0)
1354                 goto error;
1355         var->is_nonneg = 0;
1356         isl_tab_kill_col(tab, var->index);
1357
1358         return tab;
1359 error:
1360         isl_tab_free(tab);
1361         return NULL;
1362 }
1363
1364 struct isl_tab *isl_tab_from_basic_map(struct isl_basic_map *bmap)
1365 {
1366         int i;
1367         struct isl_tab *tab;
1368
1369         if (!bmap)
1370                 return NULL;
1371         tab = isl_tab_alloc(bmap->ctx,
1372                             isl_basic_map_total_dim(bmap) + bmap->n_ineq + 1,
1373                             isl_basic_map_total_dim(bmap), 0);
1374         if (!tab)
1375                 return NULL;
1376         tab->rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
1377         if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
1378                 return isl_tab_mark_empty(tab);
1379         for (i = 0; i < bmap->n_eq; ++i) {
1380                 tab = add_eq(tab, bmap->eq[i]);
1381                 if (!tab)
1382                         return tab;
1383         }
1384         for (i = 0; i < bmap->n_ineq; ++i) {
1385                 tab = isl_tab_add_ineq(tab, bmap->ineq[i]);
1386                 if (!tab || tab->empty)
1387                         return tab;
1388         }
1389         return tab;
1390 }
1391
1392 struct isl_tab *isl_tab_from_basic_set(struct isl_basic_set *bset)
1393 {
1394         return isl_tab_from_basic_map((struct isl_basic_map *)bset);
1395 }
1396
1397 /* Construct a tableau corresponding to the recession cone of "bset".
1398  */
1399 struct isl_tab *isl_tab_from_recession_cone(struct isl_basic_set *bset)
1400 {
1401         isl_int cst;
1402         int i;
1403         struct isl_tab *tab;
1404
1405         if (!bset)
1406                 return NULL;
1407         tab = isl_tab_alloc(bset->ctx, bset->n_eq + bset->n_ineq,
1408                                 isl_basic_set_total_dim(bset), 0);
1409         if (!tab)
1410                 return NULL;
1411         tab->rational = ISL_F_ISSET(bset, ISL_BASIC_SET_RATIONAL);
1412
1413         isl_int_init(cst);
1414         for (i = 0; i < bset->n_eq; ++i) {
1415                 isl_int_swap(bset->eq[i][0], cst);
1416                 tab = add_eq(tab, bset->eq[i]);
1417                 isl_int_swap(bset->eq[i][0], cst);
1418                 if (!tab)
1419                         goto done;
1420         }
1421         for (i = 0; i < bset->n_ineq; ++i) {
1422                 int r;
1423                 isl_int_swap(bset->ineq[i][0], cst);
1424                 r = isl_tab_add_row(tab, bset->ineq[i]);
1425                 isl_int_swap(bset->ineq[i][0], cst);
1426                 if (r < 0)
1427                         goto error;
1428                 tab->con[r].is_nonneg = 1;
1429                 isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]);
1430         }
1431 done:
1432         isl_int_clear(cst);
1433         return tab;
1434 error:
1435         isl_int_clear(cst);
1436         isl_tab_free(tab);
1437         return NULL;
1438 }
1439
1440 /* Assuming "tab" is the tableau of a cone, check if the cone is
1441  * bounded, i.e., if it is empty or only contains the origin.
1442  */
1443 int isl_tab_cone_is_bounded(struct isl_tab *tab)
1444 {
1445         int i;
1446
1447         if (!tab)
1448                 return -1;
1449         if (tab->empty)
1450                 return 1;
1451         if (tab->n_dead == tab->n_col)
1452                 return 1;
1453
1454         for (;;) {
1455                 for (i = tab->n_redundant; i < tab->n_row; ++i) {
1456                         struct isl_tab_var *var;
1457                         var = isl_tab_var_from_row(tab, i);
1458                         if (!var->is_nonneg)
1459                                 continue;
1460                         if (sign_of_max(tab, var) != 0)
1461                                 return 0;
1462                         close_row(tab, var);
1463                         break;
1464                 }
1465                 if (tab->n_dead == tab->n_col)
1466                         return 1;
1467                 if (i == tab->n_row)
1468                         return 0;
1469         }
1470 }
1471
1472 int isl_tab_sample_is_integer(struct isl_tab *tab)
1473 {
1474         int i;
1475
1476         if (!tab)
1477                 return -1;
1478
1479         for (i = 0; i < tab->n_var; ++i) {
1480                 int row;
1481                 if (!tab->var[i].is_row)
1482                         continue;
1483                 row = tab->var[i].index;
1484                 if (!isl_int_is_divisible_by(tab->mat->row[row][1],
1485                                                 tab->mat->row[row][0]))
1486                         return 0;
1487         }
1488         return 1;
1489 }
1490
1491 static struct isl_vec *extract_integer_sample(struct isl_tab *tab)
1492 {
1493         int i;
1494         struct isl_vec *vec;
1495
1496         vec = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
1497         if (!vec)
1498                 return NULL;
1499
1500         isl_int_set_si(vec->block.data[0], 1);
1501         for (i = 0; i < tab->n_var; ++i) {
1502                 if (!tab->var[i].is_row)
1503                         isl_int_set_si(vec->block.data[1 + i], 0);
1504                 else {
1505                         int row = tab->var[i].index;
1506                         isl_int_divexact(vec->block.data[1 + i],
1507                                 tab->mat->row[row][1], tab->mat->row[row][0]);
1508                 }
1509         }
1510
1511         return vec;
1512 }
1513
1514 struct isl_vec *isl_tab_get_sample_value(struct isl_tab *tab)
1515 {
1516         int i;
1517         struct isl_vec *vec;
1518         isl_int m;
1519
1520         if (!tab)
1521                 return NULL;
1522
1523         vec = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
1524         if (!vec)
1525                 return NULL;
1526
1527         isl_int_init(m);
1528
1529         isl_int_set_si(vec->block.data[0], 1);
1530         for (i = 0; i < tab->n_var; ++i) {
1531                 int row;
1532                 if (!tab->var[i].is_row) {
1533                         isl_int_set_si(vec->block.data[1 + i], 0);
1534                         continue;
1535                 }
1536                 row = tab->var[i].index;
1537                 isl_int_gcd(m, vec->block.data[0], tab->mat->row[row][0]);
1538                 isl_int_divexact(m, tab->mat->row[row][0], m);
1539                 isl_seq_scale(vec->block.data, vec->block.data, m, 1 + i);
1540                 isl_int_divexact(m, vec->block.data[0], tab->mat->row[row][0]);
1541                 isl_int_mul(vec->block.data[1 + i], m, tab->mat->row[row][1]);
1542         }
1543         vec = isl_vec_normalize(vec);
1544
1545         isl_int_clear(m);
1546         return vec;
1547 }
1548
1549 /* Update "bmap" based on the results of the tableau "tab".
1550  * In particular, implicit equalities are made explicit, redundant constraints
1551  * are removed and if the sample value happens to be integer, it is stored
1552  * in "bmap" (unless "bmap" already had an integer sample).
1553  *
1554  * The tableau is assumed to have been created from "bmap" using
1555  * isl_tab_from_basic_map.
1556  */
1557 struct isl_basic_map *isl_basic_map_update_from_tab(struct isl_basic_map *bmap,
1558         struct isl_tab *tab)
1559 {
1560         int i;
1561         unsigned n_eq;
1562
1563         if (!bmap)
1564                 return NULL;
1565         if (!tab)
1566                 return bmap;
1567
1568         n_eq = tab->n_eq;
1569         if (tab->empty)
1570                 bmap = isl_basic_map_set_to_empty(bmap);
1571         else
1572                 for (i = bmap->n_ineq - 1; i >= 0; --i) {
1573                         if (isl_tab_is_equality(tab, n_eq + i))
1574                                 isl_basic_map_inequality_to_equality(bmap, i);
1575                         else if (isl_tab_is_redundant(tab, n_eq + i))
1576                                 isl_basic_map_drop_inequality(bmap, i);
1577                 }
1578         if (!tab->rational &&
1579             !bmap->sample && isl_tab_sample_is_integer(tab))
1580                 bmap->sample = extract_integer_sample(tab);
1581         return bmap;
1582 }
1583
1584 struct isl_basic_set *isl_basic_set_update_from_tab(struct isl_basic_set *bset,
1585         struct isl_tab *tab)
1586 {
1587         return (struct isl_basic_set *)isl_basic_map_update_from_tab(
1588                 (struct isl_basic_map *)bset, tab);
1589 }
1590
1591 /* Given a non-negative variable "var", add a new non-negative variable
1592  * that is the opposite of "var", ensuring that var can only attain the
1593  * value zero.
1594  * If var = n/d is a row variable, then the new variable = -n/d.
1595  * If var is a column variables, then the new variable = -var.
1596  * If the new variable cannot attain non-negative values, then
1597  * the resulting tableau is empty.
1598  * Otherwise, we know the value will be zero and we close the row.
1599  */
1600 static struct isl_tab *cut_to_hyperplane(struct isl_tab *tab,
1601         struct isl_tab_var *var)
1602 {
1603         unsigned r;
1604         isl_int *row;
1605         int sgn;
1606         unsigned off = 2 + tab->M;
1607
1608         if (var->is_zero)
1609                 return tab;
1610         isl_assert(tab->mat->ctx, !var->is_redundant, goto error);
1611
1612         if (isl_tab_extend_cons(tab, 1) < 0)
1613                 goto error;
1614
1615         r = tab->n_con;
1616         tab->con[r].index = tab->n_row;
1617         tab->con[r].is_row = 1;
1618         tab->con[r].is_nonneg = 0;
1619         tab->con[r].is_zero = 0;
1620         tab->con[r].is_redundant = 0;
1621         tab->con[r].frozen = 0;
1622         tab->con[r].negated = 0;
1623         tab->row_var[tab->n_row] = ~r;
1624         row = tab->mat->row[tab->n_row];
1625
1626         if (var->is_row) {
1627                 isl_int_set(row[0], tab->mat->row[var->index][0]);
1628                 isl_seq_neg(row + 1,
1629                             tab->mat->row[var->index] + 1, 1 + tab->n_col);
1630         } else {
1631                 isl_int_set_si(row[0], 1);
1632                 isl_seq_clr(row + 1, 1 + tab->n_col);
1633                 isl_int_set_si(row[off + var->index], -1);
1634         }
1635
1636         tab->n_row++;
1637         tab->n_con++;
1638         isl_tab_push_var(tab, isl_tab_undo_allocate, &tab->con[r]);
1639
1640         sgn = sign_of_max(tab, &tab->con[r]);
1641         if (sgn < 0)
1642                 return isl_tab_mark_empty(tab);
1643         tab->con[r].is_nonneg = 1;
1644         isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]);
1645         /* sgn == 0 */
1646         close_row(tab, &tab->con[r]);
1647
1648         return tab;
1649 error:
1650         isl_tab_free(tab);
1651         return NULL;
1652 }
1653
1654 /* Given a tableau "tab" and an inequality constraint "con" of the tableau,
1655  * relax the inequality by one.  That is, the inequality r >= 0 is replaced
1656  * by r' = r + 1 >= 0.
1657  * If r is a row variable, we simply increase the constant term by one
1658  * (taking into account the denominator).
1659  * If r is a column variable, then we need to modify each row that
1660  * refers to r = r' - 1 by substituting this equality, effectively
1661  * subtracting the coefficient of the column from the constant.
1662  */
1663 struct isl_tab *isl_tab_relax(struct isl_tab *tab, int con)
1664 {
1665         struct isl_tab_var *var;
1666         unsigned off = 2 + tab->M;
1667
1668         if (!tab)
1669                 return NULL;
1670
1671         var = &tab->con[con];
1672
1673         if (!var->is_row && !max_is_manifestly_unbounded(tab, var))
1674                 to_row(tab, var, 1);
1675
1676         if (var->is_row)
1677                 isl_int_add(tab->mat->row[var->index][1],
1678                     tab->mat->row[var->index][1], tab->mat->row[var->index][0]);
1679         else {
1680                 int i;
1681
1682                 for (i = 0; i < tab->n_row; ++i) {
1683                         if (isl_int_is_zero(tab->mat->row[i][off + var->index]))
1684                                 continue;
1685                         isl_int_sub(tab->mat->row[i][1], tab->mat->row[i][1],
1686                             tab->mat->row[i][off + var->index]);
1687                 }
1688
1689         }
1690
1691         isl_tab_push_var(tab, isl_tab_undo_relax, var);
1692
1693         return tab;
1694 }
1695
1696 struct isl_tab *isl_tab_select_facet(struct isl_tab *tab, int con)
1697 {
1698         if (!tab)
1699                 return NULL;
1700
1701         return cut_to_hyperplane(tab, &tab->con[con]);
1702 }
1703
1704 static int may_be_equality(struct isl_tab *tab, int row)
1705 {
1706         unsigned off = 2 + tab->M;
1707         return (tab->rational ? isl_int_is_zero(tab->mat->row[row][1])
1708                               : isl_int_lt(tab->mat->row[row][1],
1709                                             tab->mat->row[row][0])) &&
1710                 isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
1711                                         tab->n_col - tab->n_dead) != -1;
1712 }
1713
1714 /* Check for (near) equalities among the constraints.
1715  * A constraint is an equality if it is non-negative and if
1716  * its maximal value is either
1717  *      - zero (in case of rational tableaus), or
1718  *      - strictly less than 1 (in case of integer tableaus)
1719  *
1720  * We first mark all non-redundant and non-dead variables that
1721  * are not frozen and not obviously not an equality.
1722  * Then we iterate over all marked variables if they can attain
1723  * any values larger than zero or at least one.
1724  * If the maximal value is zero, we mark any column variables
1725  * that appear in the row as being zero and mark the row as being redundant.
1726  * Otherwise, if the maximal value is strictly less than one (and the
1727  * tableau is integer), then we restrict the value to being zero
1728  * by adding an opposite non-negative variable.
1729  */
1730 struct isl_tab *isl_tab_detect_implicit_equalities(struct isl_tab *tab)
1731 {
1732         int i;
1733         unsigned n_marked;
1734
1735         if (!tab)
1736                 return NULL;
1737         if (tab->empty)
1738                 return tab;
1739         if (tab->n_dead == tab->n_col)
1740                 return tab;
1741
1742         n_marked = 0;
1743         for (i = tab->n_redundant; i < tab->n_row; ++i) {
1744                 struct isl_tab_var *var = isl_tab_var_from_row(tab, i);
1745                 var->marked = !var->frozen && var->is_nonneg &&
1746                         may_be_equality(tab, i);
1747                 if (var->marked)
1748                         n_marked++;
1749         }
1750         for (i = tab->n_dead; i < tab->n_col; ++i) {
1751                 struct isl_tab_var *var = var_from_col(tab, i);
1752                 var->marked = !var->frozen && var->is_nonneg;
1753                 if (var->marked)
1754                         n_marked++;
1755         }
1756         while (n_marked) {
1757                 struct isl_tab_var *var;
1758                 for (i = tab->n_redundant; i < tab->n_row; ++i) {
1759                         var = isl_tab_var_from_row(tab, i);
1760                         if (var->marked)
1761                                 break;
1762                 }
1763                 if (i == tab->n_row) {
1764                         for (i = tab->n_dead; i < tab->n_col; ++i) {
1765                                 var = var_from_col(tab, i);
1766                                 if (var->marked)
1767                                         break;
1768                         }
1769                         if (i == tab->n_col)
1770                                 break;
1771                 }
1772                 var->marked = 0;
1773                 n_marked--;
1774                 if (sign_of_max(tab, var) == 0)
1775                         close_row(tab, var);
1776                 else if (!tab->rational && !at_least_one(tab, var)) {
1777                         tab = cut_to_hyperplane(tab, var);
1778                         return isl_tab_detect_implicit_equalities(tab);
1779                 }
1780                 for (i = tab->n_redundant; i < tab->n_row; ++i) {
1781                         var = isl_tab_var_from_row(tab, i);
1782                         if (!var->marked)
1783                                 continue;
1784                         if (may_be_equality(tab, i))
1785                                 continue;
1786                         var->marked = 0;
1787                         n_marked--;
1788                 }
1789         }
1790
1791         return tab;
1792 }
1793
1794 /* Check for (near) redundant constraints.
1795  * A constraint is redundant if it is non-negative and if
1796  * its minimal value (temporarily ignoring the non-negativity) is either
1797  *      - zero (in case of rational tableaus), or
1798  *      - strictly larger than -1 (in case of integer tableaus)
1799  *
1800  * We first mark all non-redundant and non-dead variables that
1801  * are not frozen and not obviously negatively unbounded.
1802  * Then we iterate over all marked variables if they can attain
1803  * any values smaller than zero or at most negative one.
1804  * If not, we mark the row as being redundant (assuming it hasn't
1805  * been detected as being obviously redundant in the mean time).
1806  */
1807 struct isl_tab *isl_tab_detect_redundant(struct isl_tab *tab)
1808 {
1809         int i;
1810         unsigned n_marked;
1811
1812         if (!tab)
1813                 return NULL;
1814         if (tab->empty)
1815                 return tab;
1816         if (tab->n_redundant == tab->n_row)
1817                 return tab;
1818
1819         n_marked = 0;
1820         for (i = tab->n_redundant; i < tab->n_row; ++i) {
1821                 struct isl_tab_var *var = isl_tab_var_from_row(tab, i);
1822                 var->marked = !var->frozen && var->is_nonneg;
1823                 if (var->marked)
1824                         n_marked++;
1825         }
1826         for (i = tab->n_dead; i < tab->n_col; ++i) {
1827                 struct isl_tab_var *var = var_from_col(tab, i);
1828                 var->marked = !var->frozen && var->is_nonneg &&
1829                         !min_is_manifestly_unbounded(tab, var);
1830                 if (var->marked)
1831                         n_marked++;
1832         }
1833         while (n_marked) {
1834                 struct isl_tab_var *var;
1835                 for (i = tab->n_redundant; i < tab->n_row; ++i) {
1836                         var = isl_tab_var_from_row(tab, i);
1837                         if (var->marked)
1838                                 break;
1839                 }
1840                 if (i == tab->n_row) {
1841                         for (i = tab->n_dead; i < tab->n_col; ++i) {
1842                                 var = var_from_col(tab, i);
1843                                 if (var->marked)
1844                                         break;
1845                         }
1846                         if (i == tab->n_col)
1847                                 break;
1848                 }
1849                 var->marked = 0;
1850                 n_marked--;
1851                 if ((tab->rational ? (sign_of_min(tab, var) >= 0)
1852                                    : !isl_tab_min_at_most_neg_one(tab, var)) &&
1853                     !var->is_redundant)
1854                         isl_tab_mark_redundant(tab, var->index);
1855                 for (i = tab->n_dead; i < tab->n_col; ++i) {
1856                         var = var_from_col(tab, i);
1857                         if (!var->marked)
1858                                 continue;
1859                         if (!min_is_manifestly_unbounded(tab, var))
1860                                 continue;
1861                         var->marked = 0;
1862                         n_marked--;
1863                 }
1864         }
1865
1866         return tab;
1867 }
1868
1869 int isl_tab_is_equality(struct isl_tab *tab, int con)
1870 {
1871         int row;
1872         unsigned off;
1873
1874         if (!tab)
1875                 return -1;
1876         if (tab->con[con].is_zero)
1877                 return 1;
1878         if (tab->con[con].is_redundant)
1879                 return 0;
1880         if (!tab->con[con].is_row)
1881                 return tab->con[con].index < tab->n_dead;
1882
1883         row = tab->con[con].index;
1884
1885         off = 2 + tab->M;
1886         return isl_int_is_zero(tab->mat->row[row][1]) &&
1887                 isl_seq_first_non_zero(tab->mat->row[row] + 2 + tab->n_dead,
1888                                         tab->n_col - tab->n_dead) == -1;
1889 }
1890
1891 /* Return the minimial value of the affine expression "f" with denominator
1892  * "denom" in *opt, *opt_denom, assuming the tableau is not empty and
1893  * the expression cannot attain arbitrarily small values.
1894  * If opt_denom is NULL, then *opt is rounded up to the nearest integer.
1895  * The return value reflects the nature of the result (empty, unbounded,
1896  * minmimal value returned in *opt).
1897  */
1898 enum isl_lp_result isl_tab_min(struct isl_tab *tab,
1899         isl_int *f, isl_int denom, isl_int *opt, isl_int *opt_denom,
1900         unsigned flags)
1901 {
1902         int r;
1903         enum isl_lp_result res = isl_lp_ok;
1904         struct isl_tab_var *var;
1905         struct isl_tab_undo *snap;
1906
1907         if (tab->empty)
1908                 return isl_lp_empty;
1909
1910         snap = isl_tab_snap(tab);
1911         r = isl_tab_add_row(tab, f);
1912         if (r < 0)
1913                 return isl_lp_error;
1914         var = &tab->con[r];
1915         isl_int_mul(tab->mat->row[var->index][0],
1916                     tab->mat->row[var->index][0], denom);
1917         for (;;) {
1918                 int row, col;
1919                 find_pivot(tab, var, var, -1, &row, &col);
1920                 if (row == var->index) {
1921                         res = isl_lp_unbounded;
1922                         break;
1923                 }
1924                 if (row == -1)
1925                         break;
1926                 isl_tab_pivot(tab, row, col);
1927         }
1928         if (ISL_FL_ISSET(flags, ISL_TAB_SAVE_DUAL)) {
1929                 int i;
1930
1931                 isl_vec_free(tab->dual);
1932                 tab->dual = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_con);
1933                 if (!tab->dual)
1934                         return isl_lp_error;
1935                 isl_int_set(tab->dual->el[0], tab->mat->row[var->index][0]);
1936                 for (i = 0; i < tab->n_con; ++i) {
1937                         int pos;
1938                         if (tab->con[i].is_row) {
1939                                 isl_int_set_si(tab->dual->el[1 + i], 0);
1940                                 continue;
1941                         }
1942                         pos = 2 + tab->M + tab->con[i].index;
1943                         if (tab->con[i].negated)
1944                                 isl_int_neg(tab->dual->el[1 + i],
1945                                             tab->mat->row[var->index][pos]);
1946                         else
1947                                 isl_int_set(tab->dual->el[1 + i],
1948                                             tab->mat->row[var->index][pos]);
1949                 }
1950         }
1951         if (opt && res == isl_lp_ok) {
1952                 if (opt_denom) {
1953                         isl_int_set(*opt, tab->mat->row[var->index][1]);
1954                         isl_int_set(*opt_denom, tab->mat->row[var->index][0]);
1955                 } else
1956                         isl_int_cdiv_q(*opt, tab->mat->row[var->index][1],
1957                                              tab->mat->row[var->index][0]);
1958         }
1959         if (isl_tab_rollback(tab, snap) < 0)
1960                 return isl_lp_error;
1961         return res;
1962 }
1963
1964 int isl_tab_is_redundant(struct isl_tab *tab, int con)
1965 {
1966         if (!tab)
1967                 return -1;
1968         if (tab->con[con].is_zero)
1969                 return 0;
1970         if (tab->con[con].is_redundant)
1971                 return 1;
1972         return tab->con[con].is_row && tab->con[con].index < tab->n_redundant;
1973 }
1974
1975 /* Take a snapshot of the tableau that can be restored by s call to
1976  * isl_tab_rollback.
1977  */
1978 struct isl_tab_undo *isl_tab_snap(struct isl_tab *tab)
1979 {
1980         if (!tab)
1981                 return NULL;
1982         tab->need_undo = 1;
1983         return tab->top;
1984 }
1985
1986 /* Undo the operation performed by isl_tab_relax.
1987  */
1988 static void unrelax(struct isl_tab *tab, struct isl_tab_var *var)
1989 {
1990         unsigned off = 2 + tab->M;
1991
1992         if (!var->is_row && !max_is_manifestly_unbounded(tab, var))
1993                 to_row(tab, var, 1);
1994
1995         if (var->is_row)
1996                 isl_int_sub(tab->mat->row[var->index][1],
1997                     tab->mat->row[var->index][1], tab->mat->row[var->index][0]);
1998         else {
1999                 int i;
2000
2001                 for (i = 0; i < tab->n_row; ++i) {
2002                         if (isl_int_is_zero(tab->mat->row[i][off + var->index]))
2003                                 continue;
2004                         isl_int_add(tab->mat->row[i][1], tab->mat->row[i][1],
2005                             tab->mat->row[i][off + var->index]);
2006                 }
2007
2008         }
2009 }
2010
2011 static void perform_undo_var(struct isl_tab *tab, struct isl_tab_undo *undo)
2012 {
2013         struct isl_tab_var *var = var_from_index(tab, undo->u.var_index);
2014         switch(undo->type) {
2015         case isl_tab_undo_nonneg:
2016                 var->is_nonneg = 0;
2017                 break;
2018         case isl_tab_undo_redundant:
2019                 var->is_redundant = 0;
2020                 tab->n_redundant--;
2021                 break;
2022         case isl_tab_undo_zero:
2023                 var->is_zero = 0;
2024                 if (!var->is_row)
2025                         tab->n_dead--;
2026                 break;
2027         case isl_tab_undo_allocate:
2028                 if (undo->u.var_index >= 0) {
2029                         isl_assert(tab->mat->ctx, !var->is_row, return);
2030                         drop_col(tab, var->index);
2031                         break;
2032                 }
2033                 if (!var->is_row) {
2034                         if (!max_is_manifestly_unbounded(tab, var))
2035                                 to_row(tab, var, 1);
2036                         else if (!min_is_manifestly_unbounded(tab, var))
2037                                 to_row(tab, var, -1);
2038                         else
2039                                 to_row(tab, var, 0);
2040                 }
2041                 drop_row(tab, var->index);
2042                 break;
2043         case isl_tab_undo_relax:
2044                 unrelax(tab, var);
2045                 break;
2046         }
2047 }
2048
2049 /* Restore the tableau to the state where the basic variables
2050  * are those in "col_var".
2051  * We first construct a list of variables that are currently in
2052  * the basis, but shouldn't.  Then we iterate over all variables
2053  * that should be in the basis and for each one that is currently
2054  * not in the basis, we exchange it with one of the elements of the
2055  * list constructed before.
2056  * We can always find an appropriate variable to pivot with because
2057  * the current basis is mapped to the old basis by a non-singular
2058  * matrix and so we can never end up with a zero row.
2059  */
2060 static int restore_basis(struct isl_tab *tab, int *col_var)
2061 {
2062         int i, j;
2063         int n_extra = 0;
2064         int *extra = NULL;      /* current columns that contain bad stuff */
2065         unsigned off = 2 + tab->M;
2066
2067         extra = isl_alloc_array(tab->mat->ctx, int, tab->n_col);
2068         if (!extra)
2069                 goto error;
2070         for (i = 0; i < tab->n_col; ++i) {
2071                 for (j = 0; j < tab->n_col; ++j)
2072                         if (tab->col_var[i] == col_var[j])
2073                                 break;
2074                 if (j < tab->n_col)
2075                         continue;
2076                 extra[n_extra++] = i;
2077         }
2078         for (i = 0; i < tab->n_col && n_extra > 0; ++i) {
2079                 struct isl_tab_var *var;
2080                 int row;
2081
2082                 for (j = 0; j < tab->n_col; ++j)
2083                         if (col_var[i] == tab->col_var[j])
2084                                 break;
2085                 if (j < tab->n_col)
2086                         continue;
2087                 var = var_from_index(tab, col_var[i]);
2088                 row = var->index;
2089                 for (j = 0; j < n_extra; ++j)
2090                         if (!isl_int_is_zero(tab->mat->row[row][off+extra[j]]))
2091                                 break;
2092                 isl_assert(tab->mat->ctx, j < n_extra, goto error);
2093                 isl_tab_pivot(tab, row, extra[j]);
2094                 extra[j] = extra[--n_extra];
2095         }
2096
2097         free(extra);
2098         free(col_var);
2099         return 0;
2100 error:
2101         free(extra);
2102         free(col_var);
2103         return -1;
2104 }
2105
2106 static int perform_undo(struct isl_tab *tab, struct isl_tab_undo *undo)
2107 {
2108         switch (undo->type) {
2109         case isl_tab_undo_empty:
2110                 tab->empty = 0;
2111                 break;
2112         case isl_tab_undo_nonneg:
2113         case isl_tab_undo_redundant:
2114         case isl_tab_undo_zero:
2115         case isl_tab_undo_allocate:
2116         case isl_tab_undo_relax:
2117                 perform_undo_var(tab, undo);
2118                 break;
2119         case isl_tab_undo_bset_eq:
2120                 isl_basic_set_free_equality(tab->bset, 1);
2121                 break;
2122         case isl_tab_undo_bset_ineq:
2123                 isl_basic_set_free_inequality(tab->bset, 1);
2124                 break;
2125         case isl_tab_undo_bset_div:
2126                 isl_basic_set_free_div(tab->bset, 1);
2127                 if (tab->samples)
2128                         tab->samples->n_col--;
2129                 break;
2130         case isl_tab_undo_saved_basis:
2131                 if (restore_basis(tab, undo->u.col_var) < 0)
2132                         return -1;
2133                 break;
2134         case isl_tab_undo_drop_sample:
2135                 tab->n_outside--;
2136                 break;
2137         default:
2138                 isl_assert(tab->mat->ctx, 0, return -1);
2139         }
2140         return 0;
2141 }
2142
2143 /* Return the tableau to the state it was in when the snapshot "snap"
2144  * was taken.
2145  */
2146 int isl_tab_rollback(struct isl_tab *tab, struct isl_tab_undo *snap)
2147 {
2148         struct isl_tab_undo *undo, *next;
2149
2150         if (!tab)
2151                 return -1;
2152
2153         tab->in_undo = 1;
2154         for (undo = tab->top; undo && undo != &tab->bottom; undo = next) {
2155                 next = undo->next;
2156                 if (undo == snap)
2157                         break;
2158                 if (perform_undo(tab, undo) < 0) {
2159                         free_undo(tab);
2160                         tab->in_undo = 0;
2161                         return -1;
2162                 }
2163                 free(undo);
2164         }
2165         tab->in_undo = 0;
2166         tab->top = undo;
2167         if (!undo)
2168                 return -1;
2169         return 0;
2170 }
2171
2172 /* The given row "row" represents an inequality violated by all
2173  * points in the tableau.  Check for some special cases of such
2174  * separating constraints.
2175  * In particular, if the row has been reduced to the constant -1,
2176  * then we know the inequality is adjacent (but opposite) to
2177  * an equality in the tableau.
2178  * If the row has been reduced to r = -1 -r', with r' an inequality
2179  * of the tableau, then the inequality is adjacent (but opposite)
2180  * to the inequality r'.
2181  */
2182 static enum isl_ineq_type separation_type(struct isl_tab *tab, unsigned row)
2183 {
2184         int pos;
2185         unsigned off = 2 + tab->M;
2186
2187         if (tab->rational)
2188                 return isl_ineq_separate;
2189
2190         if (!isl_int_is_one(tab->mat->row[row][0]))
2191                 return isl_ineq_separate;
2192         if (!isl_int_is_negone(tab->mat->row[row][1]))
2193                 return isl_ineq_separate;
2194
2195         pos = isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
2196                                         tab->n_col - tab->n_dead);
2197         if (pos == -1)
2198                 return isl_ineq_adj_eq;
2199
2200         if (!isl_int_is_negone(tab->mat->row[row][off + tab->n_dead + pos]))
2201                 return isl_ineq_separate;
2202
2203         pos = isl_seq_first_non_zero(
2204                         tab->mat->row[row] + off + tab->n_dead + pos + 1,
2205                         tab->n_col - tab->n_dead - pos - 1);
2206
2207         return pos == -1 ? isl_ineq_adj_ineq : isl_ineq_separate;
2208 }
2209
2210 /* Check the effect of inequality "ineq" on the tableau "tab".
2211  * The result may be
2212  *      isl_ineq_redundant:     satisfied by all points in the tableau
2213  *      isl_ineq_separate:      satisfied by no point in the tableau
2214  *      isl_ineq_cut:           satisfied by some by not all points
2215  *      isl_ineq_adj_eq:        adjacent to an equality
2216  *      isl_ineq_adj_ineq:      adjacent to an inequality.
2217  */
2218 enum isl_ineq_type isl_tab_ineq_type(struct isl_tab *tab, isl_int *ineq)
2219 {
2220         enum isl_ineq_type type = isl_ineq_error;
2221         struct isl_tab_undo *snap = NULL;
2222         int con;
2223         int row;
2224
2225         if (!tab)
2226                 return isl_ineq_error;
2227
2228         if (isl_tab_extend_cons(tab, 1) < 0)
2229                 return isl_ineq_error;
2230
2231         snap = isl_tab_snap(tab);
2232
2233         con = isl_tab_add_row(tab, ineq);
2234         if (con < 0)
2235                 goto error;
2236
2237         row = tab->con[con].index;
2238         if (isl_tab_row_is_redundant(tab, row))
2239                 type = isl_ineq_redundant;
2240         else if (isl_int_is_neg(tab->mat->row[row][1]) &&
2241                  (tab->rational ||
2242                     isl_int_abs_ge(tab->mat->row[row][1],
2243                                    tab->mat->row[row][0]))) {
2244                 if (at_least_zero(tab, &tab->con[con]))
2245                         type = isl_ineq_cut;
2246                 else
2247                         type = separation_type(tab, row);
2248         } else if (tab->rational ? (sign_of_min(tab, &tab->con[con]) < 0)
2249                              : isl_tab_min_at_most_neg_one(tab, &tab->con[con]))
2250                 type = isl_ineq_cut;
2251         else
2252                 type = isl_ineq_redundant;
2253
2254         if (isl_tab_rollback(tab, snap))
2255                 return isl_ineq_error;
2256         return type;
2257 error:
2258         isl_tab_rollback(tab, snap);
2259         return isl_ineq_error;
2260 }
2261
2262 void isl_tab_dump(struct isl_tab *tab, FILE *out, int indent)
2263 {
2264         unsigned r, c;
2265         int i;
2266
2267         if (!tab) {
2268                 fprintf(out, "%*snull tab\n", indent, "");
2269                 return;
2270         }
2271         fprintf(out, "%*sn_redundant: %d, n_dead: %d", indent, "",
2272                 tab->n_redundant, tab->n_dead);
2273         if (tab->rational)
2274                 fprintf(out, ", rational");
2275         if (tab->empty)
2276                 fprintf(out, ", empty");
2277         fprintf(out, "\n");
2278         fprintf(out, "%*s[", indent, "");
2279         for (i = 0; i < tab->n_var; ++i) {
2280                 if (i)
2281                         fprintf(out, (i == tab->n_param ||
2282                                       i == tab->n_var - tab->n_div) ? "; "
2283                                                                     : ", ");
2284                 fprintf(out, "%c%d%s", tab->var[i].is_row ? 'r' : 'c',
2285                                         tab->var[i].index,
2286                                         tab->var[i].is_zero ? " [=0]" :
2287                                         tab->var[i].is_redundant ? " [R]" : "");
2288         }
2289         fprintf(out, "]\n");
2290         fprintf(out, "%*s[", indent, "");
2291         for (i = 0; i < tab->n_con; ++i) {
2292                 if (i)
2293                         fprintf(out, ", ");
2294                 fprintf(out, "%c%d%s", tab->con[i].is_row ? 'r' : 'c',
2295                                         tab->con[i].index,
2296                                         tab->con[i].is_zero ? " [=0]" :
2297                                         tab->con[i].is_redundant ? " [R]" : "");
2298         }
2299         fprintf(out, "]\n");
2300         fprintf(out, "%*s[", indent, "");
2301         for (i = 0; i < tab->n_row; ++i) {
2302                 const char *sign = "";
2303                 if (i)
2304                         fprintf(out, ", ");
2305                 if (tab->row_sign) {
2306                         if (tab->row_sign[i] == isl_tab_row_unknown)
2307                                 sign = "?";
2308                         else if (tab->row_sign[i] == isl_tab_row_neg)
2309                                 sign = "-";
2310                         else if (tab->row_sign[i] == isl_tab_row_pos)
2311                                 sign = "+";
2312                         else
2313                                 sign = "+-";
2314                 }
2315                 fprintf(out, "r%d: %d%s%s", i, tab->row_var[i],
2316                     isl_tab_var_from_row(tab, i)->is_nonneg ? " [>=0]" : "", sign);
2317         }
2318         fprintf(out, "]\n");
2319         fprintf(out, "%*s[", indent, "");
2320         for (i = 0; i < tab->n_col; ++i) {
2321                 if (i)
2322                         fprintf(out, ", ");
2323                 fprintf(out, "c%d: %d%s", i, tab->col_var[i],
2324                     var_from_col(tab, i)->is_nonneg ? " [>=0]" : "");
2325         }
2326         fprintf(out, "]\n");
2327         r = tab->mat->n_row;
2328         tab->mat->n_row = tab->n_row;
2329         c = tab->mat->n_col;
2330         tab->mat->n_col = 2 + tab->M + tab->n_col;
2331         isl_mat_dump(tab->mat, out, indent);
2332         tab->mat->n_row = r;
2333         tab->mat->n_col = c;
2334         if (tab->bset)
2335                 isl_basic_set_dump(tab->bset, out, indent);
2336 }