isl_tab.c

   1 /*
   2  * Copyright 2008-2009 Katholieke Universiteit Leuven
   3  *
   4  * Use of this software is governed by the GNU LGPLv2.1 license
   5  *
   6  * Written by Sven Verdoolaege, K.U.Leuven, Departement
   7  * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
   8  */
   9
  10 #include <isl_ctx_private.h>
  11 #include <isl_mat_private.h>
  12 #include "isl_map_private.h"
  13 #include "isl_tab.h"
  14 #include <isl/seq.h>
  15
  16 /*
  17  * The implementation of tableaus in this file was inspired by Section 8
  18  * of David Detlefs, Greg Nelson and James B. Saxe, "Simplify: a theorem
  19  * prover for program checking".
  20  */
  21
  22 struct isl_tab *isl_tab_alloc(struct isl_ctx *ctx,
  23         unsigned n_row, unsigned n_var, unsigned M)
  24 {
  25         int i;
  26         struct isl_tab *tab;
  27         unsigned off = 2 + M;
  28
  29         tab = isl_calloc_type(ctx, struct isl_tab);
  30         if (!tab)
  31                 return NULL;
  32         tab->mat = isl_mat_alloc(ctx, n_row, off + n_var);
  33         if (!tab->mat)
  34                 goto error;
  35         tab->var = isl_alloc_array(ctx, struct isl_tab_var, n_var);
  36         if (!tab->var)
  37                 goto error;
  38         tab->con = isl_alloc_array(ctx, struct isl_tab_var, n_row);
  39         if (!tab->con)
  40                 goto error;
  41         tab->col_var = isl_alloc_array(ctx, int, n_var);
  42         if (!tab->col_var)
  43                 goto error;
  44         tab->row_var = isl_alloc_array(ctx, int, n_row);
  45         if (!tab->row_var)
  46                 goto error;
  47         for (i = 0; i < n_var; ++i) {
  48                 tab->var[i].index = i;
  49                 tab->var[i].is_row = 0;
  50                 tab->var[i].is_nonneg = 0;
  51                 tab->var[i].is_zero = 0;
  52                 tab->var[i].is_redundant = 0;
  53                 tab->var[i].frozen = 0;
  54                 tab->var[i].negated = 0;
  55                 tab->col_var[i] = i;
  56         }
  57         tab->n_row = 0;
  58         tab->n_con = 0;
  59         tab->n_eq = 0;
  60         tab->max_con = n_row;
  61         tab->n_col = n_var;
  62         tab->n_var = n_var;
  63         tab->max_var = n_var;
  64         tab->n_param = 0;
  65         tab->n_div = 0;
  66         tab->n_dead = 0;
  67         tab->n_redundant = 0;
  68         tab->strict_redundant = 0;
  69         tab->need_undo = 0;
  70         tab->rational = 0;
  71         tab->empty = 0;
  72         tab->in_undo = 0;
  73         tab->M = M;
  74         tab->cone = 0;
  75         tab->bottom.type = isl_tab_undo_bottom;
  76         tab->bottom.next = NULL;
  77         tab->top = &tab->bottom;
  78
  79         tab->n_zero = 0;
  80         tab->n_unbounded = 0;
  81         tab->basis = NULL;
  82
  83         return tab;
  84 error:
  85         isl_tab_free(tab);
  86         return NULL;
  87 }
  88
  89 int isl_tab_extend_cons(struct isl_tab *tab, unsigned n_new)
  90 {
  91         unsigned off;
  92
  93         if (!tab)
  94                 return -1;
  95
  96         off = 2 + tab->M;
  97
  98         if (tab->max_con < tab->n_con + n_new) {
  99                 struct isl_tab_var *con;
 100
 101                 con = isl_realloc_array(tab->mat->ctx, tab->con,
 102                                     struct isl_tab_var, tab->max_con + n_new);
 103                 if (!con)
 104                         return -1;
 105                 tab->con = con;
 106                 tab->max_con += n_new;
 107         }
 108         if (tab->mat->n_row < tab->n_row + n_new) {
 109                 int *row_var;
 110
 111                 tab->mat = isl_mat_extend(tab->mat,
 112                                         tab->n_row + n_new, off + tab->n_col);
 113                 if (!tab->mat)
 114                         return -1;
 115                 row_var = isl_realloc_array(tab->mat->ctx, tab->row_var,
 116                                             int, tab->mat->n_row);
 117                 if (!row_var)
 118                         return -1;
 119                 tab->row_var = row_var;
 120                 if (tab->row_sign) {
 121                         enum isl_tab_row_sign *s;
 122                         s = isl_realloc_array(tab->mat->ctx, tab->row_sign,
 123                                         enum isl_tab_row_sign, tab->mat->n_row);
 124                         if (!s)
 125                                 return -1;
 126                         tab->row_sign = s;
 127                 }
 128         }
 129         return 0;
 130 }
 131
 132 /* Make room for at least n_new extra variables.
 133  * Return -1 if anything went wrong.
 134  */
 135 int isl_tab_extend_vars(struct isl_tab *tab, unsigned n_new)
 136 {
 137         struct isl_tab_var *var;
 138         unsigned off = 2 + tab->M;
 139
 140         if (tab->max_var < tab->n_var + n_new) {
 141                 var = isl_realloc_array(tab->mat->ctx, tab->var,
 142                                     struct isl_tab_var, tab->n_var + n_new);
 143                 if (!var)
 144                         return -1;
 145                 tab->var = var;
 146                 tab->max_var += n_new;
 147         }
 148
 149         if (tab->mat->n_col < off + tab->n_col + n_new) {
 150                 int *p;
 151
 152                 tab->mat = isl_mat_extend(tab->mat,
 153                                     tab->mat->n_row, off + tab->n_col + n_new);
 154                 if (!tab->mat)
 155                         return -1;
 156                 p = isl_realloc_array(tab->mat->ctx, tab->col_var,
 157                                             int, tab->n_col + n_new);
 158                 if (!p)
 159                         return -1;
 160                 tab->col_var = p;
 161         }
 162
 163         return 0;
 164 }
 165
 166 struct isl_tab *isl_tab_extend(struct isl_tab *tab, unsigned n_new)
 167 {
 168         if (isl_tab_extend_cons(tab, n_new) >= 0)
 169                 return tab;
 170
 171         isl_tab_free(tab);
 172         return NULL;
 173 }
 174
 175 static void free_undo_record(struct isl_tab_undo *undo)
 176 {
 177         switch (undo->type) {
 178         case isl_tab_undo_saved_basis:
 179                 free(undo->u.col_var);
 180                 break;
 181         default:;
 182         }
 183         free(undo);
 184 }
 185
 186 static void free_undo(struct isl_tab *tab)
 187 {
 188         struct isl_tab_undo *undo, *next;
 189
 190         for (undo = tab->top; undo && undo != &tab->bottom; undo = next) {
 191                 next = undo->next;
 192                 free_undo_record(undo);
 193         }
 194         tab->top = undo;
 195 }
 196
 197 void isl_tab_free(struct isl_tab *tab)
 198 {
 199         if (!tab)
 200                 return;
 201         free_undo(tab);
 202         isl_mat_free(tab->mat);
 203         isl_vec_free(tab->dual);
 204         isl_basic_map_free(tab->bmap);
 205         free(tab->var);
 206         free(tab->con);
 207         free(tab->row_var);
 208         free(tab->col_var);
 209         free(tab->row_sign);
 210         isl_mat_free(tab->samples);
 211         free(tab->sample_index);
 212         isl_mat_free(tab->basis);
 213         free(tab);
 214 }
 215
 216 struct isl_tab *isl_tab_dup(struct isl_tab *tab)
 217 {
 218         int i;
 219         struct isl_tab *dup;
 220         unsigned off;
 221
 222         if (!tab)
 223                 return NULL;
 224
 225         off = 2 + tab->M;
 226         dup = isl_calloc_type(tab->mat->ctx, struct isl_tab);
 227         if (!dup)
 228                 return NULL;
 229         dup->mat = isl_mat_dup(tab->mat);
 230         if (!dup->mat)
 231                 goto error;
 232         dup->var = isl_alloc_array(tab->mat->ctx, struct isl_tab_var, tab->max_var);
 233         if (!dup->var)
 234                 goto error;
 235         for (i = 0; i < tab->n_var; ++i)
 236                 dup->var[i] = tab->var[i];
 237         dup->con = isl_alloc_array(tab->mat->ctx, struct isl_tab_var, tab->max_con);
 238         if (!dup->con)
 239                 goto error;
 240         for (i = 0; i < tab->n_con; ++i)
 241                 dup->con[i] = tab->con[i];
 242         dup->col_var = isl_alloc_array(tab->mat->ctx, int, tab->mat->n_col - off);
 243         if (!dup->col_var)
 244                 goto error;
 245         for (i = 0; i < tab->n_col; ++i)
 246                 dup->col_var[i] = tab->col_var[i];
 247         dup->row_var = isl_alloc_array(tab->mat->ctx, int, tab->mat->n_row);
 248         if (!dup->row_var)
 249                 goto error;
 250         for (i = 0; i < tab->n_row; ++i)
 251                 dup->row_var[i] = tab->row_var[i];
 252         if (tab->row_sign) {
 253                 dup->row_sign = isl_alloc_array(tab->mat->ctx, enum isl_tab_row_sign,
 254                                                 tab->mat->n_row);
 255                 if (!dup->row_sign)
 256                         goto error;
 257                 for (i = 0; i < tab->n_row; ++i)
 258                         dup->row_sign[i] = tab->row_sign[i];
 259         }
 260         if (tab->samples) {
 261                 dup->samples = isl_mat_dup(tab->samples);
 262                 if (!dup->samples)
 263                         goto error;
 264                 dup->sample_index = isl_alloc_array(tab->mat->ctx, int,
 265                                                         tab->samples->n_row);
 266                 if (!dup->sample_index)
 267                         goto error;
 268                 dup->n_sample = tab->n_sample;
 269                 dup->n_outside = tab->n_outside;
 270         }
 271         dup->n_row = tab->n_row;
 272         dup->n_con = tab->n_con;
 273         dup->n_eq = tab->n_eq;
 274         dup->max_con = tab->max_con;
 275         dup->n_col = tab->n_col;
 276         dup->n_var = tab->n_var;
 277         dup->max_var = tab->max_var;
 278         dup->n_param = tab->n_param;
 279         dup->n_div = tab->n_div;
 280         dup->n_dead = tab->n_dead;
 281         dup->n_redundant = tab->n_redundant;
 282         dup->rational = tab->rational;
 283         dup->empty = tab->empty;
 284         dup->strict_redundant = 0;
 285         dup->need_undo = 0;
 286         dup->in_undo = 0;
 287         dup->M = tab->M;
 288         tab->cone = tab->cone;
 289         dup->bottom.type = isl_tab_undo_bottom;
 290         dup->bottom.next = NULL;
 291         dup->top = &dup->bottom;
 292
 293         dup->n_zero = tab->n_zero;
 294         dup->n_unbounded = tab->n_unbounded;
 295         dup->basis = isl_mat_dup(tab->basis);
 296
 297         return dup;
 298 error:
 299         isl_tab_free(dup);
 300         return NULL;
 301 }
 302
 303 /* Construct the coefficient matrix of the product tableau
 304  * of two tableaus.
 305  * mat{1,2} is the coefficient matrix of tableau {1,2}
 306  * row{1,2} is the number of rows in tableau {1,2}
 307  * col{1,2} is the number of columns in tableau {1,2}
 308  * off is the offset to the coefficient column (skipping the
 309  *      denominator, the constant term and the big parameter if any)
 310  * r{1,2} is the number of redundant rows in tableau {1,2}
 311  * d{1,2} is the number of dead columns in tableau {1,2}
 312  *
 313  * The order of the rows and columns in the result is as explained
 314  * in isl_tab_product.
 315  */
 316 static struct isl_mat *tab_mat_product(struct isl_mat *mat1,
 317         struct isl_mat *mat2, unsigned row1, unsigned row2,
 318         unsigned col1, unsigned col2,
 319         unsigned off, unsigned r1, unsigned r2, unsigned d1, unsigned d2)
 320 {
 321         int i;
 322         struct isl_mat *prod;
 323         unsigned n;
 324
 325         prod = isl_mat_alloc(mat1->ctx, mat1->n_row + mat2->n_row,
 326                                         off + col1 + col2);
 327         if (!prod)
 328                 return NULL;
 329
 330         n = 0;
 331         for (i = 0; i < r1; ++i) {
 332                 isl_seq_cpy(prod->row[n + i], mat1->row[i], off + d1);
 333                 isl_seq_clr(prod->row[n + i] + off + d1, d2);
 334                 isl_seq_cpy(prod->row[n + i] + off + d1 + d2,
 335                                 mat1->row[i] + off + d1, col1 - d1);
 336                 isl_seq_clr(prod->row[n + i] + off + col1 + d1, col2 - d2);
 337         }
 338
 339         n += r1;
 340         for (i = 0; i < r2; ++i) {
 341                 isl_seq_cpy(prod->row[n + i], mat2->row[i], off);
 342                 isl_seq_clr(prod->row[n + i] + off, d1);
 343                 isl_seq_cpy(prod->row[n + i] + off + d1,
 344                             mat2->row[i] + off, d2);
 345                 isl_seq_clr(prod->row[n + i] + off + d1 + d2, col1 - d1);
 346                 isl_seq_cpy(prod->row[n + i] + off + col1 + d1,
 347                             mat2->row[i] + off + d2, col2 - d2);
 348         }
 349
 350         n += r2;
 351         for (i = 0; i < row1 - r1; ++i) {
 352                 isl_seq_cpy(prod->row[n + i], mat1->row[r1 + i], off + d1);
 353                 isl_seq_clr(prod->row[n + i] + off + d1, d2);
 354                 isl_seq_cpy(prod->row[n + i] + off + d1 + d2,
 355                                 mat1->row[r1 + i] + off + d1, col1 - d1);
 356                 isl_seq_clr(prod->row[n + i] + off + col1 + d1, col2 - d2);
 357         }
 358
 359         n += row1 - r1;
 360         for (i = 0; i < row2 - r2; ++i) {
 361                 isl_seq_cpy(prod->row[n + i], mat2->row[r2 + i], off);
 362                 isl_seq_clr(prod->row[n + i] + off, d1);
 363                 isl_seq_cpy(prod->row[n + i] + off + d1,
 364                             mat2->row[r2 + i] + off, d2);
 365                 isl_seq_clr(prod->row[n + i] + off + d1 + d2, col1 - d1);
 366                 isl_seq_cpy(prod->row[n + i] + off + col1 + d1,
 367                             mat2->row[r2 + i] + off + d2, col2 - d2);
 368         }
 369
 370         return prod;
 371 }
 372
 373 /* Update the row or column index of a variable that corresponds
 374  * to a variable in the first input tableau.
 375  */
 376 static void update_index1(struct isl_tab_var *var,
 377         unsigned r1, unsigned r2, unsigned d1, unsigned d2)
 378 {
 379         if (var->index == -1)
 380                 return;
 381         if (var->is_row && var->index >= r1)
 382                 var->index += r2;
 383         if (!var->is_row && var->index >= d1)
 384                 var->index += d2;
 385 }
 386
 387 /* Update the row or column index of a variable that corresponds
 388  * to a variable in the second input tableau.
 389  */
 390 static void update_index2(struct isl_tab_var *var,
 391         unsigned row1, unsigned col1,
 392         unsigned r1, unsigned r2, unsigned d1, unsigned d2)
 393 {
 394         if (var->index == -1)
 395                 return;
 396         if (var->is_row) {
 397                 if (var->index < r2)
 398                         var->index += r1;
 399                 else
 400                         var->index += row1;
 401         } else {
 402                 if (var->index < d2)
 403                         var->index += d1;
 404                 else
 405                         var->index += col1;
 406         }
 407 }
 408
 409 /* Create a tableau that represents the Cartesian product of the sets
 410  * represented by tableaus tab1 and tab2.
 411  * The order of the rows in the product is
 412  *      - redundant rows of tab1
 413  *      - redundant rows of tab2
 414  *      - non-redundant rows of tab1
 415  *      - non-redundant rows of tab2
 416  * The order of the columns is
 417  *      - denominator
 418  *      - constant term
 419  *      - coefficient of big parameter, if any
 420  *      - dead columns of tab1
 421  *      - dead columns of tab2
 422  *      - live columns of tab1
 423  *      - live columns of tab2
 424  * The order of the variables and the constraints is a concatenation
 425  * of order in the two input tableaus.
 426  */
 427 struct isl_tab *isl_tab_product(struct isl_tab *tab1, struct isl_tab *tab2)
 428 {
 429         int i;
 430         struct isl_tab *prod;
 431         unsigned off;
 432         unsigned r1, r2, d1, d2;
 433
 434         if (!tab1 || !tab2)
 435                 return NULL;
 436
 437         isl_assert(tab1->mat->ctx, tab1->M == tab2->M, return NULL);
 438         isl_assert(tab1->mat->ctx, tab1->rational == tab2->rational, return NULL);
 439         isl_assert(tab1->mat->ctx, tab1->cone == tab2->cone, return NULL);
 440         isl_assert(tab1->mat->ctx, !tab1->row_sign, return NULL);
 441         isl_assert(tab1->mat->ctx, !tab2->row_sign, return NULL);
 442         isl_assert(tab1->mat->ctx, tab1->n_param == 0, return NULL);
 443         isl_assert(tab1->mat->ctx, tab2->n_param == 0, return NULL);
 444         isl_assert(tab1->mat->ctx, tab1->n_div == 0, return NULL);
 445         isl_assert(tab1->mat->ctx, tab2->n_div == 0, return NULL);
 446
 447         off = 2 + tab1->M;
 448         r1 = tab1->n_redundant;
 449         r2 = tab2->n_redundant;
 450         d1 = tab1->n_dead;
 451         d2 = tab2->n_dead;
 452         prod = isl_calloc_type(tab1->mat->ctx, struct isl_tab);
 453         if (!prod)
 454                 return NULL;
 455         prod->mat = tab_mat_product(tab1->mat, tab2->mat,
 456                                 tab1->n_row, tab2->n_row,
 457                                 tab1->n_col, tab2->n_col, off, r1, r2, d1, d2);
 458         if (!prod->mat)
 459                 goto error;
 460         prod->var = isl_alloc_array(tab1->mat->ctx, struct isl_tab_var,
 461                                         tab1->max_var + tab2->max_var);
 462         if (!prod->var)
 463                 goto error;
 464         for (i = 0; i < tab1->n_var; ++i) {
 465                 prod->var[i] = tab1->var[i];
 466                 update_index1(&prod->var[i], r1, r2, d1, d2);
 467         }
 468         for (i = 0; i < tab2->n_var; ++i) {
 469                 prod->var[tab1->n_var + i] = tab2->var[i];
 470                 update_index2(&prod->var[tab1->n_var + i],
 471                                 tab1->n_row, tab1->n_col,
 472                                 r1, r2, d1, d2);
 473         }
 474         prod->con = isl_alloc_array(tab1->mat->ctx, struct isl_tab_var,
 475                                         tab1->max_con +  tab2->max_con);
 476         if (!prod->con)
 477                 goto error;
 478         for (i = 0; i < tab1->n_con; ++i) {
 479                 prod->con[i] = tab1->con[i];
 480                 update_index1(&prod->con[i], r1, r2, d1, d2);
 481         }
 482         for (i = 0; i < tab2->n_con; ++i) {
 483                 prod->con[tab1->n_con + i] = tab2->con[i];
 484                 update_index2(&prod->con[tab1->n_con + i],
 485                                 tab1->n_row, tab1->n_col,
 486                                 r1, r2, d1, d2);
 487         }
 488         prod->col_var = isl_alloc_array(tab1->mat->ctx, int,
 489                                         tab1->n_col + tab2->n_col);
 490         if (!prod->col_var)
 491                 goto error;
 492         for (i = 0; i < tab1->n_col; ++i) {
 493                 int pos = i < d1 ? i : i + d2;
 494                 prod->col_var[pos] = tab1->col_var[i];
 495         }
 496         for (i = 0; i < tab2->n_col; ++i) {
 497                 int pos = i < d2 ? d1 + i : tab1->n_col + i;
 498                 int t = tab2->col_var[i];
 499                 if (t >= 0)
 500                         t += tab1->n_var;
 501                 else
 502                         t -= tab1->n_con;
 503                 prod->col_var[pos] = t;
 504         }
 505         prod->row_var = isl_alloc_array(tab1->mat->ctx, int,
 506                                         tab1->mat->n_row + tab2->mat->n_row);
 507         if (!prod->row_var)
 508                 goto error;
 509         for (i = 0; i < tab1->n_row; ++i) {
 510                 int pos = i < r1 ? i : i + r2;
 511                 prod->row_var[pos] = tab1->row_var[i];
 512         }
 513         for (i = 0; i < tab2->n_row; ++i) {
 514                 int pos = i < r2 ? r1 + i : tab1->n_row + i;
 515                 int t = tab2->row_var[i];
 516                 if (t >= 0)
 517                         t += tab1->n_var;
 518                 else
 519                         t -= tab1->n_con;
 520                 prod->row_var[pos] = t;
 521         }
 522         prod->samples = NULL;
 523         prod->sample_index = NULL;
 524         prod->n_row = tab1->n_row + tab2->n_row;
 525         prod->n_con = tab1->n_con + tab2->n_con;
 526         prod->n_eq = 0;
 527         prod->max_con = tab1->max_con + tab2->max_con;
 528         prod->n_col = tab1->n_col + tab2->n_col;
 529         prod->n_var = tab1->n_var + tab2->n_var;
 530         prod->max_var = tab1->max_var + tab2->max_var;
 531         prod->n_param = 0;
 532         prod->n_div = 0;
 533         prod->n_dead = tab1->n_dead + tab2->n_dead;
 534         prod->n_redundant = tab1->n_redundant + tab2->n_redundant;
 535         prod->rational = tab1->rational;
 536         prod->empty = tab1->empty || tab2->empty;
 537         prod->strict_redundant = tab1->strict_redundant || tab2->strict_redundant;
 538         prod->need_undo = 0;
 539         prod->in_undo = 0;
 540         prod->M = tab1->M;
 541         prod->cone = tab1->cone;
 542         prod->bottom.type = isl_tab_undo_bottom;
 543         prod->bottom.next = NULL;
 544         prod->top = &prod->bottom;
 545
 546         prod->n_zero = 0;
 547         prod->n_unbounded = 0;
 548         prod->basis = NULL;
 549
 550         return prod;
 551 error:
 552         isl_tab_free(prod);
 553         return NULL;
 554 }
 555
 556 static struct isl_tab_var *var_from_index(struct isl_tab *tab, int i)
 557 {
 558         if (i >= 0)
 559                 return &tab->var[i];
 560         else
 561                 return &tab->con[~i];
 562 }
 563
 564 struct isl_tab_var *isl_tab_var_from_row(struct isl_tab *tab, int i)
 565 {
 566         return var_from_index(tab, tab->row_var[i]);
 567 }
 568
 569 static struct isl_tab_var *var_from_col(struct isl_tab *tab, int i)
 570 {
 571         return var_from_index(tab, tab->col_var[i]);
 572 }
 573
 574 /* Check if there are any upper bounds on column variable "var",
 575  * i.e., non-negative rows where var appears with a negative coefficient.
 576  * Return 1 if there are no such bounds.
 577  */
 578 static int max_is_manifestly_unbounded(struct isl_tab *tab,
 579         struct isl_tab_var *var)
 580 {
 581         int i;
 582         unsigned off = 2 + tab->M;
 583
 584         if (var->is_row)
 585                 return 0;
 586         for (i = tab->n_redundant; i < tab->n_row; ++i) {
 587                 if (!isl_int_is_neg(tab->mat->row[i][off + var->index]))
 588                         continue;
 589                 if (isl_tab_var_from_row(tab, i)->is_nonneg)
 590                         return 0;
 591         }
 592         return 1;
 593 }
 594
 595 /* Check if there are any lower bounds on column variable "var",
 596  * i.e., non-negative rows where var appears with a positive coefficient.
 597  * Return 1 if there are no such bounds.
 598  */
 599 static int min_is_manifestly_unbounded(struct isl_tab *tab,
 600         struct isl_tab_var *var)
 601 {
 602         int i;
 603         unsigned off = 2 + tab->M;
 604
 605         if (var->is_row)
 606                 return 0;
 607         for (i = tab->n_redundant; i < tab->n_row; ++i) {
 608                 if (!isl_int_is_pos(tab->mat->row[i][off + var->index]))
 609                         continue;
 610                 if (isl_tab_var_from_row(tab, i)->is_nonneg)
 611                         return 0;
 612         }
 613         return 1;
 614 }
 615
 616 static int row_cmp(struct isl_tab *tab, int r1, int r2, int c, isl_int t)
 617 {
 618         unsigned off = 2 + tab->M;
 619
 620         if (tab->M) {
 621                 int s;
 622                 isl_int_mul(t, tab->mat->row[r1][2], tab->mat->row[r2][off+c]);
 623                 isl_int_submul(t, tab->mat->row[r2][2], tab->mat->row[r1][off+c]);
 624                 s = isl_int_sgn(t);
 625                 if (s)
 626                         return s;
 627         }
 628         isl_int_mul(t, tab->mat->row[r1][1], tab->mat->row[r2][off + c]);
 629         isl_int_submul(t, tab->mat->row[r2][1], tab->mat->row[r1][off + c]);
 630         return isl_int_sgn(t);
 631 }
 632
 633 /* Given the index of a column "c", return the index of a row
 634  * that can be used to pivot the column in, with either an increase
 635  * (sgn > 0) or a decrease (sgn < 0) of the corresponding variable.
 636  * If "var" is not NULL, then the row returned will be different from
 637  * the one associated with "var".
 638  *
 639  * Each row in the tableau is of the form
 640  *
 641  *      x_r = a_r0 + \sum_i a_ri x_i
 642  *
 643  * Only rows with x_r >= 0 and with the sign of a_ri opposite to "sgn"
 644  * impose any limit on the increase or decrease in the value of x_c
 645  * and this bound is equal to a_r0 / |a_rc|.  We are therefore looking
 646  * for the row with the smallest (most stringent) such bound.
 647  * Note that the common denominator of each row drops out of the fraction.
 648  * To check if row j has a smaller bound than row r, i.e.,
 649  * a_j0 / |a_jc| < a_r0 / |a_rc| or a_j0 |a_rc| < a_r0 |a_jc|,
 650  * we check if -sign(a_jc) (a_j0 a_rc - a_r0 a_jc) < 0,
 651  * where -sign(a_jc) is equal to "sgn".
 652  */
 653 static int pivot_row(struct isl_tab *tab,
 654         struct isl_tab_var *var, int sgn, int c)
 655 {
 656         int j, r, tsgn;
 657         isl_int t;
 658         unsigned off = 2 + tab->M;
 659
 660         isl_int_init(t);
 661         r = -1;
 662         for (j = tab->n_redundant; j < tab->n_row; ++j) {
 663                 if (var && j == var->index)
 664                         continue;
 665                 if (!isl_tab_var_from_row(tab, j)->is_nonneg)
 666                         continue;
 667                 if (sgn * isl_int_sgn(tab->mat->row[j][off + c]) >= 0)
 668                         continue;
 669                 if (r < 0) {
 670                         r = j;
 671                         continue;
 672                 }
 673                 tsgn = sgn * row_cmp(tab, r, j, c, t);
 674                 if (tsgn < 0 || (tsgn == 0 &&
 675                                             tab->row_var[j] < tab->row_var[r]))
 676                         r = j;
 677         }
 678         isl_int_clear(t);
 679         return r;
 680 }
 681
 682 /* Find a pivot (row and col) that will increase (sgn > 0) or decrease
 683  * (sgn < 0) the value of row variable var.
 684  * If not NULL, then skip_var is a row variable that should be ignored
 685  * while looking for a pivot row.  It is usually equal to var.
 686  *
 687  * As the given row in the tableau is of the form
 688  *
 689  *      x_r = a_r0 + \sum_i a_ri x_i
 690  *
 691  * we need to find a column such that the sign of a_ri is equal to "sgn"
 692  * (such that an increase in x_i will have the desired effect) or a
 693  * column with a variable that may attain negative values.
 694  * If a_ri is positive, then we need to move x_i in the same direction
 695  * to obtain the desired effect.  Otherwise, x_i has to move in the
 696  * opposite direction.
 697  */
 698 static void find_pivot(struct isl_tab *tab,
 699         struct isl_tab_var *var, struct isl_tab_var *skip_var,
 700         int sgn, int *row, int *col)
 701 {
 702         int j, r, c;
 703         isl_int *tr;
 704
 705         *row = *col = -1;
 706
 707         isl_assert(tab->mat->ctx, var->is_row, return);
 708         tr = tab->mat->row[var->index] + 2 + tab->M;
 709
 710         c = -1;
 711         for (j = tab->n_dead; j < tab->n_col; ++j) {
 712                 if (isl_int_is_zero(tr[j]))
 713                         continue;
 714                 if (isl_int_sgn(tr[j]) != sgn &&
 715                     var_from_col(tab, j)->is_nonneg)
 716                         continue;
 717                 if (c < 0 || tab->col_var[j] < tab->col_var[c])
 718                         c = j;
 719         }
 720         if (c < 0)
 721                 return;
 722
 723         sgn *= isl_int_sgn(tr[c]);
 724         r = pivot_row(tab, skip_var, sgn, c);
 725         *row = r < 0 ? var->index : r;
 726         *col = c;
 727 }
 728
 729 /* Return 1 if row "row" represents an obviously redundant inequality.
 730  * This means
 731  *      - it represents an inequality or a variable
 732  *      - that is the sum of a non-negative sample value and a positive
 733  *        combination of zero or more non-negative constraints.
 734  */
 735 int isl_tab_row_is_redundant(struct isl_tab *tab, int row)
 736 {
 737         int i;
 738         unsigned off = 2 + tab->M;
 739
 740         if (tab->row_var[row] < 0 && !isl_tab_var_from_row(tab, row)->is_nonneg)
 741                 return 0;
 742
 743         if (isl_int_is_neg(tab->mat->row[row][1]))
 744                 return 0;
 745         if (tab->strict_redundant && isl_int_is_zero(tab->mat->row[row][1]))
 746                 return 0;
 747         if (tab->M && isl_int_is_neg(tab->mat->row[row][2]))
 748                 return 0;
 749
 750         for (i = tab->n_dead; i < tab->n_col; ++i) {
 751                 if (isl_int_is_zero(tab->mat->row[row][off + i]))
 752                         continue;
 753                 if (tab->col_var[i] >= 0)
 754                         return 0;
 755                 if (isl_int_is_neg(tab->mat->row[row][off + i]))
 756                         return 0;
 757                 if (!var_from_col(tab, i)->is_nonneg)
 758                         return 0;
 759         }
 760         return 1;
 761 }
 762
 763 static void swap_rows(struct isl_tab *tab, int row1, int row2)
 764 {
 765         int t;
 766         enum isl_tab_row_sign s;
 767
 768         t = tab->row_var[row1];
 769         tab->row_var[row1] = tab->row_var[row2];
 770         tab->row_var[row2] = t;
 771         isl_tab_var_from_row(tab, row1)->index = row1;
 772         isl_tab_var_from_row(tab, row2)->index = row2;
 773         tab->mat = isl_mat_swap_rows(tab->mat, row1, row2);
 774
 775         if (!tab->row_sign)
 776                 return;
 777         s = tab->row_sign[row1];
 778         tab->row_sign[row1] = tab->row_sign[row2];
 779         tab->row_sign[row2] = s;
 780 }
 781
 782 static int push_union(struct isl_tab *tab,
 783         enum isl_tab_undo_type type, union isl_tab_undo_val u) WARN_UNUSED;
 784 static int push_union(struct isl_tab *tab,
 785         enum isl_tab_undo_type type, union isl_tab_undo_val u)
 786 {
 787         struct isl_tab_undo *undo;
 788
 789         if (!tab->need_undo)
 790                 return 0;
 791
 792         undo = isl_alloc_type(tab->mat->ctx, struct isl_tab_undo);
 793         if (!undo)
 794                 return -1;
 795         undo->type = type;
 796         undo->u = u;
 797         undo->next = tab->top;
 798         tab->top = undo;
 799
 800         return 0;
 801 }
 802
 803 int isl_tab_push_var(struct isl_tab *tab,
 804         enum isl_tab_undo_type type, struct isl_tab_var *var)
 805 {
 806         union isl_tab_undo_val u;
 807         if (var->is_row)
 808                 u.var_index = tab->row_var[var->index];
 809         else
 810                 u.var_index = tab->col_var[var->index];
 811         return push_union(tab, type, u);
 812 }
 813
 814 int isl_tab_push(struct isl_tab *tab, enum isl_tab_undo_type type)
 815 {
 816         union isl_tab_undo_val u = { 0 };
 817         return push_union(tab, type, u);
 818 }
 819
 820 /* Push a record on the undo stack describing the current basic
 821  * variables, so that the this state can be restored during rollback.
 822  */
 823 int isl_tab_push_basis(struct isl_tab *tab)
 824 {
 825         int i;
 826         union isl_tab_undo_val u;
 827
 828         u.col_var = isl_alloc_array(tab->mat->ctx, int, tab->n_col);
 829         if (!u.col_var)
 830                 return -1;
 831         for (i = 0; i < tab->n_col; ++i)
 832                 u.col_var[i] = tab->col_var[i];
 833         return push_union(tab, isl_tab_undo_saved_basis, u);
 834 }
 835
 836 int isl_tab_push_callback(struct isl_tab *tab, struct isl_tab_callback *callback)
 837 {
 838         union isl_tab_undo_val u;
 839         u.callback = callback;
 840         return push_union(tab, isl_tab_undo_callback, u);
 841 }
 842
 843 struct isl_tab *isl_tab_init_samples(struct isl_tab *tab)
 844 {
 845         if (!tab)
 846                 return NULL;
 847
 848         tab->n_sample = 0;
 849         tab->n_outside = 0;
 850         tab->samples = isl_mat_alloc(tab->mat->ctx, 1, 1 + tab->n_var);
 851         if (!tab->samples)
 852                 goto error;
 853         tab->sample_index = isl_alloc_array(tab->mat->ctx, int, 1);
 854         if (!tab->sample_index)
 855                 goto error;
 856         return tab;
 857 error:
 858         isl_tab_free(tab);
 859         return NULL;
 860 }
 861
 862 struct isl_tab *isl_tab_add_sample(struct isl_tab *tab,
 863         __isl_take isl_vec *sample)
 864 {
 865         if (!tab || !sample)
 866                 goto error;
 867
 868         if (tab->n_sample + 1 > tab->samples->n_row) {
 869                 int *t = isl_realloc_array(tab->mat->ctx,
 870                             tab->sample_index, int, tab->n_sample + 1);
 871                 if (!t)
 872                         goto error;
 873                 tab->sample_index = t;
 874         }
 875
 876         tab->samples = isl_mat_extend(tab->samples,
 877                                 tab->n_sample + 1, tab->samples->n_col);
 878         if (!tab->samples)
 879                 goto error;
 880
 881         isl_seq_cpy(tab->samples->row[tab->n_sample], sample->el, sample->size);
 882         isl_vec_free(sample);
 883         tab->sample_index[tab->n_sample] = tab->n_sample;
 884         tab->n_sample++;
 885
 886         return tab;
 887 error:
 888         isl_vec_free(sample);
 889         isl_tab_free(tab);
 890         return NULL;
 891 }
 892
 893 struct isl_tab *isl_tab_drop_sample(struct isl_tab *tab, int s)
 894 {
 895         if (s != tab->n_outside) {
 896                 int t = tab->sample_index[tab->n_outside];
 897                 tab->sample_index[tab->n_outside] = tab->sample_index[s];
 898                 tab->sample_index[s] = t;
 899                 isl_mat_swap_rows(tab->samples, tab->n_outside, s);
 900         }
 901         tab->n_outside++;
 902         if (isl_tab_push(tab, isl_tab_undo_drop_sample) < 0) {
 903                 isl_tab_free(tab);
 904                 return NULL;
 905         }
 906
 907         return tab;
 908 }
 909
 910 /* Record the current number of samples so that we can remove newer
 911  * samples during a rollback.
 912  */
 913 int isl_tab_save_samples(struct isl_tab *tab)
 914 {
 915         union isl_tab_undo_val u;
 916
 917         if (!tab)
 918                 return -1;
 919
 920         u.n = tab->n_sample;
 921         return push_union(tab, isl_tab_undo_saved_samples, u);
 922 }
 923
 924 /* Mark row with index "row" as being redundant.
 925  * If we may need to undo the operation or if the row represents
 926  * a variable of the original problem, the row is kept,
 927  * but no longer considered when looking for a pivot row.
 928  * Otherwise, the row is simply removed.
 929  *
 930  * The row may be interchanged with some other row.  If it
 931  * is interchanged with a later row, return 1.  Otherwise return 0.
 932  * If the rows are checked in order in the calling function,
 933  * then a return value of 1 means that the row with the given
 934  * row number may now contain a different row that hasn't been checked yet.
 935  */
 936 int isl_tab_mark_redundant(struct isl_tab *tab, int row)
 937 {
 938         struct isl_tab_var *var = isl_tab_var_from_row(tab, row);
 939         var->is_redundant = 1;
 940         isl_assert(tab->mat->ctx, row >= tab->n_redundant, return -1);
 941         if (tab->need_undo || tab->row_var[row] >= 0) {
 942                 if (tab->row_var[row] >= 0 && !var->is_nonneg) {
 943                         var->is_nonneg = 1;
 944                         if (isl_tab_push_var(tab, isl_tab_undo_nonneg, var) < 0)
 945                                 return -1;
 946                 }
 947                 if (row != tab->n_redundant)
 948                         swap_rows(tab, row, tab->n_redundant);
 949                 tab->n_redundant++;
 950                 return isl_tab_push_var(tab, isl_tab_undo_redundant, var);
 951         } else {
 952                 if (row != tab->n_row - 1)
 953                         swap_rows(tab, row, tab->n_row - 1);
 954                 isl_tab_var_from_row(tab, tab->n_row - 1)->index = -1;
 955                 tab->n_row--;
 956                 return 1;
 957         }
 958 }
 959
 960 int isl_tab_mark_empty(struct isl_tab *tab)
 961 {
 962         if (!tab)
 963                 return -1;
 964         if (!tab->empty && tab->need_undo)
 965                 if (isl_tab_push(tab, isl_tab_undo_empty) < 0)
 966                         return -1;
 967         tab->empty = 1;
 968         return 0;
 969 }
 970
 971 int isl_tab_freeze_constraint(struct isl_tab *tab, int con)
 972 {
 973         struct isl_tab_var *var;
 974
 975         if (!tab)
 976                 return -1;
 977
 978         var = &tab->con[con];
 979         if (var->frozen)
 980                 return 0;
 981         if (var->index < 0)
 982                 return 0;
 983         var->frozen = 1;
 984
 985         if (tab->need_undo)
 986                 return isl_tab_push_var(tab, isl_tab_undo_freeze, var);
 987
 988         return 0;
 989 }
 990
 991 /* Update the rows signs after a pivot of "row" and "col", with "row_sgn"
 992  * the original sign of the pivot element.
 993  * We only keep track of row signs during PILP solving and in this case
 994  * we only pivot a row with negative sign (meaning the value is always
 995  * non-positive) using a positive pivot element.
 996  *
 997  * For each row j, the new value of the parametric constant is equal to
 998  *
 999  *      a_j0 - a_jc a_r0/a_rc
1000  *
1001  * where a_j0 is the original parametric constant, a_rc is the pivot element,
1002  * a_r0 is the parametric constant of the pivot row and a_jc is the
1003  * pivot column entry of the row j.
1004  * Since a_r0 is non-positive and a_rc is positive, the sign of row j
1005  * remains the same if a_jc has the same sign as the row j or if
1006  * a_jc is zero.  In all other cases, we reset the sign to "unknown".
1007  */
1008 static void update_row_sign(struct isl_tab *tab, int row, int col, int row_sgn)
1009 {
1010         int i;
1011         struct isl_mat *mat = tab->mat;
1012         unsigned off = 2 + tab->M;
1013
1014         if (!tab->row_sign)
1015                 return;
1016
1017         if (tab->row_sign[row] == 0)
1018                 return;
1019         isl_assert(mat->ctx, row_sgn > 0, return);
1020         isl_assert(mat->ctx, tab->row_sign[row] == isl_tab_row_neg, return);
1021         tab->row_sign[row] = isl_tab_row_pos;
1022         for (i = 0; i < tab->n_row; ++i) {
1023                 int s;
1024                 if (i == row)
1025                         continue;
1026                 s = isl_int_sgn(mat->row[i][off + col]);
1027                 if (!s)
1028                         continue;
1029                 if (!tab->row_sign[i])
1030                         continue;
1031                 if (s < 0 && tab->row_sign[i] == isl_tab_row_neg)
1032                         continue;
1033                 if (s > 0 && tab->row_sign[i] == isl_tab_row_pos)
1034                         continue;
1035                 tab->row_sign[i] = isl_tab_row_unknown;
1036         }
1037 }
1038
1039 /* Given a row number "row" and a column number "col", pivot the tableau
1040  * such that the associated variables are interchanged.
1041  * The given row in the tableau expresses
1042  *
1043  *      x_r = a_r0 + \sum_i a_ri x_i
1044  *
1045  * or
1046  *
1047  *      x_c = 1/a_rc x_r - a_r0/a_rc + sum_{i \ne r} -a_ri/a_rc
1048  *
1049  * Substituting this equality into the other rows
1050  *
1051  *      x_j = a_j0 + \sum_i a_ji x_i
1052  *
1053  * with a_jc \ne 0, we obtain
1054  *
1055  *      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
1056  *
1057  * The tableau
1058  *
1059  *      n_rc/d_r                n_ri/d_r
1060  *      n_jc/d_j                n_ji/d_j
1061  *
1062  * where i is any other column and j is any other row,
1063  * is therefore transformed into
1064  *
1065  * s(n_rc)d_r/|n_rc|            -s(n_rc)n_ri/|n_rc|
1066  * 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)
1067  *
1068  * The transformation is performed along the following steps
1069  *
1070  *      d_r/n_rc                n_ri/n_rc
1071  *      n_jc/d_j                n_ji/d_j
1072  *
1073  *      s(n_rc)d_r/|n_rc|       -s(n_rc)n_ri/|n_rc|
1074  *      n_jc/d_j                n_ji/d_j
1075  *
1076  *      s(n_rc)d_r/|n_rc|       -s(n_rc)n_ri/|n_rc|
1077  *      n_jc/(|n_rc| d_j)       n_ji/(|n_rc| d_j)
1078  *
1079  *      s(n_rc)d_r/|n_rc|       -s(n_rc)n_ri/|n_rc|
1080  *      n_jc/(|n_rc| d_j)       (n_ji |n_rc|)/(|n_rc| d_j)
1081  *
1082  *      s(n_rc)d_r/|n_rc|       -s(n_rc)n_ri/|n_rc|
1083  *      n_jc/(|n_rc| d_j)       (n_ji |n_rc| - s(n_rc)n_jc n_ri)/(|n_rc| d_j)
1084  *
1085  * s(n_rc)d_r/|n_rc|            -s(n_rc)n_ri/|n_rc|
1086  * 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)
1087  *
1088  */
1089 int isl_tab_pivot(struct isl_tab *tab, int row, int col)
1090 {
1091         int i, j;
1092         int sgn;
1093         int t;
1094         struct isl_mat *mat = tab->mat;
1095         struct isl_tab_var *var;
1096         unsigned off = 2 + tab->M;
1097
1098         if (tab->mat->ctx->abort) {
1099                 isl_ctx_set_error(tab->mat->ctx, isl_error_abort);
1100                 return -1;
1101         }
1102
1103         isl_int_swap(mat->row[row][0], mat->row[row][off + col]);
1104         sgn = isl_int_sgn(mat->row[row][0]);
1105         if (sgn < 0) {
1106                 isl_int_neg(mat->row[row][0], mat->row[row][0]);
1107                 isl_int_neg(mat->row[row][off + col], mat->row[row][off + col]);
1108         } else
1109                 for (j = 0; j < off - 1 + tab->n_col; ++j) {
1110                         if (j == off - 1 + col)
1111                                 continue;
1112                         isl_int_neg(mat->row[row][1 + j], mat->row[row][1 + j]);
1113                 }
1114         if (!isl_int_is_one(mat->row[row][0]))
1115                 isl_seq_normalize(mat->ctx, mat->row[row], off + tab->n_col);
1116         for (i = 0; i < tab->n_row; ++i) {
1117                 if (i == row)
1118                         continue;
1119                 if (isl_int_is_zero(mat->row[i][off + col]))
1120                         continue;
1121                 isl_int_mul(mat->row[i][0], mat->row[i][0], mat->row[row][0]);
1122                 for (j = 0; j < off - 1 + tab->n_col; ++j) {
1123                         if (j == off - 1 + col)
1124                                 continue;
1125                         isl_int_mul(mat->row[i][1 + j],
1126                                     mat->row[i][1 + j], mat->row[row][0]);
1127                         isl_int_addmul(mat->row[i][1 + j],
1128                                     mat->row[i][off + col], mat->row[row][1 + j]);
1129                 }
1130                 isl_int_mul(mat->row[i][off + col],
1131                             mat->row[i][off + col], mat->row[row][off + col]);
1132                 if (!isl_int_is_one(mat->row[i][0]))
1133                         isl_seq_normalize(mat->ctx, mat->row[i], off + tab->n_col);
1134         }
1135         t = tab->row_var[row];
1136         tab->row_var[row] = tab->col_var[col];
1137         tab->col_var[col] = t;
1138         var = isl_tab_var_from_row(tab, row);
1139         var->is_row = 1;
1140         var->index = row;
1141         var = var_from_col(tab, col);
1142         var->is_row = 0;
1143         var->index = col;
1144         update_row_sign(tab, row, col, sgn);
1145         if (tab->in_undo)
1146                 return 0;
1147         for (i = tab->n_redundant; i < tab->n_row; ++i) {
1148                 if (isl_int_is_zero(mat->row[i][off + col]))
1149                         continue;
1150                 if (!isl_tab_var_from_row(tab, i)->frozen &&
1151                     isl_tab_row_is_redundant(tab, i)) {
1152                         int redo = isl_tab_mark_redundant(tab, i);
1153                         if (redo < 0)
1154                                 return -1;
1155                         if (redo)
1156                                 --i;
1157                 }
1158         }
1159         return 0;
1160 }
1161
1162 /* If "var" represents a column variable, then pivot is up (sgn > 0)
1163  * or down (sgn < 0) to a row.  The variable is assumed not to be
1164  * unbounded in the specified direction.
1165  * If sgn = 0, then the variable is unbounded in both directions,
1166  * and we pivot with any row we can find.
1167  */
1168 static int to_row(struct isl_tab *tab, struct isl_tab_var *var, int sign) WARN_UNUSED;
1169 static int to_row(struct isl_tab *tab, struct isl_tab_var *var, int sign)
1170 {
1171         int r;
1172         unsigned off = 2 + tab->M;
1173
1174         if (var->is_row)
1175                 return 0;
1176
1177         if (sign == 0) {
1178                 for (r = tab->n_redundant; r < tab->n_row; ++r)
1179                         if (!isl_int_is_zero(tab->mat->row[r][off+var->index]))
1180                                 break;
1181                 isl_assert(tab->mat->ctx, r < tab->n_row, return -1);
1182         } else {
1183                 r = pivot_row(tab, NULL, sign, var->index);
1184                 isl_assert(tab->mat->ctx, r >= 0, return -1);
1185         }
1186
1187         return isl_tab_pivot(tab, r, var->index);
1188 }
1189
1190 static void check_table(struct isl_tab *tab)
1191 {
1192         int i;
1193
1194         if (tab->empty)
1195                 return;
1196         for (i = tab->n_redundant; i < tab->n_row; ++i) {
1197                 struct isl_tab_var *var;
1198                 var = isl_tab_var_from_row(tab, i);
1199                 if (!var->is_nonneg)
1200                         continue;
1201                 if (tab->M) {
1202                         isl_assert(tab->mat->ctx,
1203                                 !isl_int_is_neg(tab->mat->row[i][2]), abort());
1204                         if (isl_int_is_pos(tab->mat->row[i][2]))
1205                                 continue;
1206                 }
1207                 isl_assert(tab->mat->ctx, !isl_int_is_neg(tab->mat->row[i][1]),
1208                                 abort());
1209         }
1210 }
1211
1212 /* Return the sign of the maximal value of "var".
1213  * If the sign is not negative, then on return from this function,
1214  * the sample value will also be non-negative.
1215  *
1216  * If "var" is manifestly unbounded wrt positive values, we are done.
1217  * Otherwise, we pivot the variable up to a row if needed
1218  * Then we continue pivoting down until either
1219  *      - no more down pivots can be performed
1220  *      - the sample value is positive
1221  *      - the variable is pivoted into a manifestly unbounded column
1222  */
1223 static int sign_of_max(struct isl_tab *tab, struct isl_tab_var *var)
1224 {
1225         int row, col;
1226
1227         if (max_is_manifestly_unbounded(tab, var))
1228                 return 1;
1229         if (to_row(tab, var, 1) < 0)
1230                 return -2;
1231         while (!isl_int_is_pos(tab->mat->row[var->index][1])) {
1232                 find_pivot(tab, var, var, 1, &row, &col);
1233                 if (row == -1)
1234                         return isl_int_sgn(tab->mat->row[var->index][1]);
1235                 if (isl_tab_pivot(tab, row, col) < 0)
1236                         return -2;
1237                 if (!var->is_row) /* manifestly unbounded */
1238                         return 1;
1239         }
1240         return 1;
1241 }
1242
1243 int isl_tab_sign_of_max(struct isl_tab *tab, int con)
1244 {
1245         struct isl_tab_var *var;
1246
1247         if (!tab)
1248                 return -2;
1249
1250         var = &tab->con[con];
1251         isl_assert(tab->mat->ctx, !var->is_redundant, return -2);
1252         isl_assert(tab->mat->ctx, !var->is_zero, return -2);
1253
1254         return sign_of_max(tab, var);
1255 }
1256
1257 static int row_is_neg(struct isl_tab *tab, int row)
1258 {
1259         if (!tab->M)
1260                 return isl_int_is_neg(tab->mat->row[row][1]);
1261         if (isl_int_is_pos(tab->mat->row[row][2]))
1262                 return 0;
1263         if (isl_int_is_neg(tab->mat->row[row][2]))
1264                 return 1;
1265         return isl_int_is_neg(tab->mat->row[row][1]);
1266 }
1267
1268 static int row_sgn(struct isl_tab *tab, int row)
1269 {
1270         if (!tab->M)
1271                 return isl_int_sgn(tab->mat->row[row][1]);
1272         if (!isl_int_is_zero(tab->mat->row[row][2]))
1273                 return isl_int_sgn(tab->mat->row[row][2]);
1274         else
1275                 return isl_int_sgn(tab->mat->row[row][1]);
1276 }
1277
1278 /* Perform pivots until the row variable "var" has a non-negative
1279  * sample value or until no more upward pivots can be performed.
1280  * Return the sign of the sample value after the pivots have been
1281  * performed.
1282  */
1283 static int restore_row(struct isl_tab *tab, struct isl_tab_var *var)
1284 {
1285         int row, col;
1286
1287         while (row_is_neg(tab, var->index)) {
1288                 find_pivot(tab, var, var, 1, &row, &col);
1289                 if (row == -1)
1290                         break;
1291                 if (isl_tab_pivot(tab, row, col) < 0)
1292                         return -2;
1293                 if (!var->is_row) /* manifestly unbounded */
1294                         return 1;
1295         }
1296         return row_sgn(tab, var->index);
1297 }
1298
1299 /* Perform pivots until we are sure that the row variable "var"
1300  * can attain non-negative values.  After return from this
1301  * function, "var" is still a row variable, but its sample
1302  * value may not be non-negative, even if the function returns 1.
1303  */
1304 static int at_least_zero(struct isl_tab *tab, struct isl_tab_var *var)
1305 {
1306         int row, col;
1307
1308         while (isl_int_is_neg(tab->mat->row[var->index][1])) {
1309                 find_pivot(tab, var, var, 1, &row, &col);
1310                 if (row == -1)
1311                         break;
1312                 if (row == var->index) /* manifestly unbounded */
1313                         return 1;
1314                 if (isl_tab_pivot(tab, row, col) < 0)
1315                         return -1;
1316         }
1317         return !isl_int_is_neg(tab->mat->row[var->index][1]);
1318 }
1319
1320 /* Return a negative value if "var" can attain negative values.
1321  * Return a non-negative value otherwise.
1322  *
1323  * If "var" is manifestly unbounded wrt negative values, we are done.
1324  * Otherwise, if var is in a column, we can pivot it down to a row.
1325  * Then we continue pivoting down until either
1326  *      - the pivot would result in a manifestly unbounded column
1327  *        => we don't perform the pivot, but simply return -1
1328  *      - no more down pivots can be performed
1329  *      - the sample value is negative
1330  * If the sample value becomes negative and the variable is supposed
1331  * to be nonnegative, then we undo the last pivot.
1332  * However, if the last pivot has made the pivoting variable
1333  * obviously redundant, then it may have moved to another row.
1334  * In that case we look for upward pivots until we reach a non-negative
1335  * value again.
1336  */
1337 static int sign_of_min(struct isl_tab *tab, struct isl_tab_var *var)
1338 {
1339         int row, col;
1340         struct isl_tab_var *pivot_var = NULL;
1341
1342         if (min_is_manifestly_unbounded(tab, var))
1343                 return -1;
1344         if (!var->is_row) {
1345                 col = var->index;
1346                 row = pivot_row(tab, NULL, -1, col);
1347                 pivot_var = var_from_col(tab, col);
1348                 if (isl_tab_pivot(tab, row, col) < 0)
1349                         return -2;
1350                 if (var->is_redundant)
1351                         return 0;
1352                 if (isl_int_is_neg(tab->mat->row[var->index][1])) {
1353                         if (var->is_nonneg) {
1354                                 if (!pivot_var->is_redundant &&
1355                                     pivot_var->index == row) {
1356                                         if (isl_tab_pivot(tab, row, col) < 0)
1357                                                 return -2;
1358                                 } else
1359                                         if (restore_row(tab, var) < -1)
1360                                                 return -2;
1361                         }
1362                         return -1;
1363                 }
1364         }
1365         if (var->is_redundant)
1366                 return 0;
1367         while (!isl_int_is_neg(tab->mat->row[var->index][1])) {
1368                 find_pivot(tab, var, var, -1, &row, &col);
1369                 if (row == var->index)
1370                         return -1;
1371                 if (row == -1)
1372                         return isl_int_sgn(tab->mat->row[var->index][1]);
1373                 pivot_var = var_from_col(tab, col);
1374                 if (isl_tab_pivot(tab, row, col) < 0)
1375                         return -2;
1376                 if (var->is_redundant)
1377                         return 0;
1378         }
1379         if (pivot_var && var->is_nonneg) {
1380                 /* pivot back to non-negative value */
1381                 if (!pivot_var->is_redundant && pivot_var->index == row) {
1382                         if (isl_tab_pivot(tab, row, col) < 0)
1383                                 return -2;
1384                 } else
1385                         if (restore_row(tab, var) < -1)
1386                                 return -2;
1387         }
1388         return -1;
1389 }
1390
1391 static int row_at_most_neg_one(struct isl_tab *tab, int row)
1392 {
1393         if (tab->M) {
1394                 if (isl_int_is_pos(tab->mat->row[row][2]))
1395                         return 0;
1396                 if (isl_int_is_neg(tab->mat->row[row][2]))
1397                         return 1;
1398         }
1399         return isl_int_is_neg(tab->mat->row[row][1]) &&
1400                isl_int_abs_ge(tab->mat->row[row][1],
1401                               tab->mat->row[row][0]);
1402 }
1403
1404 /* Return 1 if "var" can attain values <= -1.
1405  * Return 0 otherwise.
1406  *
1407  * The sample value of "var" is assumed to be non-negative when the
1408  * the function is called.  If 1 is returned then the constraint
1409  * is not redundant and the sample value is made non-negative again before
1410  * the function returns.
1411  */
1412 int isl_tab_min_at_most_neg_one(struct isl_tab *tab, struct isl_tab_var *var)
1413 {
1414         int row, col;
1415         struct isl_tab_var *pivot_var;
1416
1417         if (min_is_manifestly_unbounded(tab, var))
1418                 return 1;
1419         if (!var->is_row) {
1420                 col = var->index;
1421                 row = pivot_row(tab, NULL, -1, col);
1422                 pivot_var = var_from_col(tab, col);
1423                 if (isl_tab_pivot(tab, row, col) < 0)
1424                         return -1;
1425                 if (var->is_redundant)
1426                         return 0;
1427                 if (row_at_most_neg_one(tab, var->index)) {
1428                         if (var->is_nonneg) {
1429                                 if (!pivot_var->is_redundant &&
1430                                     pivot_var->index == row) {
1431                                         if (isl_tab_pivot(tab, row, col) < 0)
1432                                                 return -1;
1433                                 } else
1434                                         if (restore_row(tab, var) < -1)
1435                                                 return -1;
1436                         }
1437                         return 1;
1438                 }
1439         }
1440         if (var->is_redundant)
1441                 return 0;
1442         do {
1443                 find_pivot(tab, var, var, -1, &row, &col);
1444                 if (row == var->index) {
1445                         if (restore_row(tab, var) < -1)
1446                                 return -1;
1447                         return 1;
1448                 }
1449                 if (row == -1)
1450                         return 0;
1451                 pivot_var = var_from_col(tab, col);
1452                 if (isl_tab_pivot(tab, row, col) < 0)
1453                         return -1;
1454                 if (var->is_redundant)
1455                         return 0;
1456         } while (!row_at_most_neg_one(tab, var->index));
1457         if (var->is_nonneg) {
1458                 /* pivot back to non-negative value */
1459                 if (!pivot_var->is_redundant && pivot_var->index == row)
1460                         if (isl_tab_pivot(tab, row, col) < 0)
1461                                 return -1;
1462                 if (restore_row(tab, var) < -1)
1463                         return -1;
1464         }
1465         return 1;
1466 }
1467
1468 /* Return 1 if "var" can attain values >= 1.
1469  * Return 0 otherwise.
1470  */
1471 static int at_least_one(struct isl_tab *tab, struct isl_tab_var *var)
1472 {
1473         int row, col;
1474         isl_int *r;
1475
1476         if (max_is_manifestly_unbounded(tab, var))
1477                 return 1;
1478         if (to_row(tab, var, 1) < 0)
1479                 return -1;
1480         r = tab->mat->row[var->index];
1481         while (isl_int_lt(r[1], r[0])) {
1482                 find_pivot(tab, var, var, 1, &row, &col);
1483                 if (row == -1)
1484                         return isl_int_ge(r[1], r[0]);
1485                 if (row == var->index) /* manifestly unbounded */
1486                         return 1;
1487                 if (isl_tab_pivot(tab, row, col) < 0)
1488                         return -1;
1489         }
1490         return 1;
1491 }
1492
1493 static void swap_cols(struct isl_tab *tab, int col1, int col2)
1494 {
1495         int t;
1496         unsigned off = 2 + tab->M;
1497         t = tab->col_var[col1];
1498         tab->col_var[col1] = tab->col_var[col2];
1499         tab->col_var[col2] = t;
1500         var_from_col(tab, col1)->index = col1;
1501         var_from_col(tab, col2)->index = col2;
1502         tab->mat = isl_mat_swap_cols(tab->mat, off + col1, off + col2);
1503 }
1504
1505 /* Mark column with index "col" as representing a zero variable.
1506  * If we may need to undo the operation the column is kept,
1507  * but no longer considered.
1508  * Otherwise, the column is simply removed.
1509  *
1510  * The column may be interchanged with some other column.  If it
1511  * is interchanged with a later column, return 1.  Otherwise return 0.
1512  * If the columns are checked in order in the calling function,
1513  * then a return value of 1 means that the column with the given
1514  * column number may now contain a different column that
1515  * hasn't been checked yet.
1516  */
1517 int isl_tab_kill_col(struct isl_tab *tab, int col)
1518 {
1519         var_from_col(tab, col)->is_zero = 1;
1520         if (tab->need_undo) {
1521                 if (isl_tab_push_var(tab, isl_tab_undo_zero,
1522                                             var_from_col(tab, col)) < 0)
1523                         return -1;
1524                 if (col != tab->n_dead)
1525                         swap_cols(tab, col, tab->n_dead);
1526                 tab->n_dead++;
1527                 return 0;
1528         } else {
1529                 if (col != tab->n_col - 1)
1530                         swap_cols(tab, col, tab->n_col - 1);
1531                 var_from_col(tab, tab->n_col - 1)->index = -1;
1532                 tab->n_col--;
1533                 return 1;
1534         }
1535 }
1536
1537 static int row_is_manifestly_non_integral(struct isl_tab *tab, int row)
1538 {
1539         unsigned off = 2 + tab->M;
1540
1541         if (tab->M && !isl_int_eq(tab->mat->row[row][2],
1542                                   tab->mat->row[row][0]))
1543                 return 0;
1544         if (isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
1545                                     tab->n_col - tab->n_dead) != -1)
1546                 return 0;
1547
1548         return !isl_int_is_divisible_by(tab->mat->row[row][1],
1549                                         tab->mat->row[row][0]);
1550 }
1551
1552 /* For integer tableaus, check if any of the coordinates are stuck
1553  * at a non-integral value.
1554  */
1555 static int tab_is_manifestly_empty(struct isl_tab *tab)
1556 {
1557         int i;
1558
1559         if (tab->empty)
1560                 return 1;
1561         if (tab->rational)
1562                 return 0;
1563
1564         for (i = 0; i < tab->n_var; ++i) {
1565                 if (!tab->var[i].is_row)
1566                         continue;
1567                 if (row_is_manifestly_non_integral(tab, tab->var[i].index))
1568                         return 1;
1569         }
1570
1571         return 0;
1572 }
1573
1574 /* Row variable "var" is non-negative and cannot attain any values
1575  * larger than zero.  This means that the coefficients of the unrestricted
1576  * column variables are zero and that the coefficients of the non-negative
1577  * column variables are zero or negative.
1578  * Each of the non-negative variables with a negative coefficient can
1579  * then also be written as the negative sum of non-negative variables
1580  * and must therefore also be zero.
1581  */
1582 static int close_row(struct isl_tab *tab, struct isl_tab_var *var) WARN_UNUSED;
1583 static int close_row(struct isl_tab *tab, struct isl_tab_var *var)
1584 {
1585         int j;
1586         struct isl_mat *mat = tab->mat;
1587         unsigned off = 2 + tab->M;
1588
1589         isl_assert(tab->mat->ctx, var->is_nonneg, return -1);
1590         var->is_zero = 1;
1591         if (tab->need_undo)
1592                 if (isl_tab_push_var(tab, isl_tab_undo_zero, var) < 0)
1593                         return -1;
1594         for (j = tab->n_dead; j < tab->n_col; ++j) {
1595                 int recheck;
1596                 if (isl_int_is_zero(mat->row[var->index][off + j]))
1597                         continue;
1598                 isl_assert(tab->mat->ctx,
1599                     isl_int_is_neg(mat->row[var->index][off + j]), return -1);
1600                 recheck = isl_tab_kill_col(tab, j);
1601                 if (recheck < 0)
1602                         return -1;
1603                 if (recheck)
1604                         --j;
1605         }
1606         if (isl_tab_mark_redundant(tab, var->index) < 0)
1607                 return -1;
1608         if (tab_is_manifestly_empty(tab) && isl_tab_mark_empty(tab) < 0)
1609                 return -1;
1610         return 0;
1611 }
1612
1613 /* Add a constraint to the tableau and allocate a row for it.
1614  * Return the index into the constraint array "con".
1615  */
1616 int isl_tab_allocate_con(struct isl_tab *tab)
1617 {
1618         int r;
1619
1620         isl_assert(tab->mat->ctx, tab->n_row < tab->mat->n_row, return -1);
1621         isl_assert(tab->mat->ctx, tab->n_con < tab->max_con, return -1);
1622
1623         r = tab->n_con;
1624         tab->con[r].index = tab->n_row;
1625         tab->con[r].is_row = 1;
1626         tab->con[r].is_nonneg = 0;
1627         tab->con[r].is_zero = 0;
1628         tab->con[r].is_redundant = 0;
1629         tab->con[r].frozen = 0;
1630         tab->con[r].negated = 0;
1631         tab->row_var[tab->n_row] = ~r;
1632
1633         tab->n_row++;
1634         tab->n_con++;
1635         if (isl_tab_push_var(tab, isl_tab_undo_allocate, &tab->con[r]) < 0)
1636                 return -1;
1637
1638         return r;
1639 }
1640
1641 /* Add a variable to the tableau and allocate a column for it.
1642  * Return the index into the variable array "var".
1643  */
1644 int isl_tab_allocate_var(struct isl_tab *tab)
1645 {
1646         int r;
1647         int i;
1648         unsigned off = 2 + tab->M;
1649
1650         isl_assert(tab->mat->ctx, tab->n_col < tab->mat->n_col, return -1);
1651         isl_assert(tab->mat->ctx, tab->n_var < tab->max_var, return -1);
1652
1653         r = tab->n_var;
1654         tab->var[r].index = tab->n_col;
1655         tab->var[r].is_row = 0;
1656         tab->var[r].is_nonneg = 0;
1657         tab->var[r].is_zero = 0;
1658         tab->var[r].is_redundant = 0;
1659         tab->var[r].frozen = 0;
1660         tab->var[r].negated = 0;
1661         tab->col_var[tab->n_col] = r;
1662
1663         for (i = 0; i < tab->n_row; ++i)
1664                 isl_int_set_si(tab->mat->row[i][off + tab->n_col], 0);
1665
1666         tab->n_var++;
1667         tab->n_col++;
1668         if (isl_tab_push_var(tab, isl_tab_undo_allocate, &tab->var[r]) < 0)
1669                 return -1;
1670
1671         return r;
1672 }
1673
1674 /* Add a row to the tableau.  The row is given as an affine combination
1675  * of the original variables and needs to be expressed in terms of the
1676  * column variables.
1677  *
1678  * We add each term in turn.
1679  * If r = n/d_r is the current sum and we need to add k x, then
1680  *      if x is a column variable, we increase the numerator of
1681  *              this column by k d_r
1682  *      if x = f/d_x is a row variable, then the new representation of r is
1683  *
1684  *               n    k f   d_x/g n + d_r/g k f   m/d_r n + m/d_g k f
1685  *              --- + --- = ------------------- = -------------------
1686  *              d_r   d_r        d_r d_x/g                m
1687  *
1688  *      with g the gcd of d_r and d_x and m the lcm of d_r and d_x.
1689  *
1690  * If tab->M is set, then, internally, each variable x is represented
1691  * as x' - M.  We then also need no subtract k d_r from the coefficient of M.
1692  */
1693 int isl_tab_add_row(struct isl_tab *tab, isl_int *line)
1694 {
1695         int i;
1696         int r;
1697         isl_int *row;
1698         isl_int a, b;
1699         unsigned off = 2 + tab->M;
1700
1701         r = isl_tab_allocate_con(tab);
1702         if (r < 0)
1703                 return -1;
1704
1705         isl_int_init(a);
1706         isl_int_init(b);
1707         row = tab->mat->row[tab->con[r].index];
1708         isl_int_set_si(row[0], 1);
1709         isl_int_set(row[1], line[0]);
1710         isl_seq_clr(row + 2, tab->M + tab->n_col);
1711         for (i = 0; i < tab->n_var; ++i) {
1712                 if (tab->var[i].is_zero)
1713                         continue;
1714                 if (tab->var[i].is_row) {
1715                         isl_int_lcm(a,
1716                                 row[0], tab->mat->row[tab->var[i].index][0]);
1717                         isl_int_swap(a, row[0]);
1718                         isl_int_divexact(a, row[0], a);
1719                         isl_int_divexact(b,
1720                                 row[0], tab->mat->row[tab->var[i].index][0]);
1721                         isl_int_mul(b, b, line[1 + i]);
1722                         isl_seq_combine(row + 1, a, row + 1,
1723                             b, tab->mat->row[tab->var[i].index] + 1,
1724                             1 + tab->M + tab->n_col);
1725                 } else
1726                         isl_int_addmul(row[off + tab->var[i].index],
1727                                                         line[1 + i], row[0]);
1728                 if (tab->M && i >= tab->n_param && i < tab->n_var - tab->n_div)
1729                         isl_int_submul(row[2], line[1 + i], row[0]);
1730         }
1731         isl_seq_normalize(tab->mat->ctx, row, off + tab->n_col);
1732         isl_int_clear(a);
1733         isl_int_clear(b);
1734
1735         if (tab->row_sign)
1736                 tab->row_sign[tab->con[r].index] = isl_tab_row_unknown;
1737
1738         return r;
1739 }
1740
1741 static int drop_row(struct isl_tab *tab, int row)
1742 {
1743         isl_assert(tab->mat->ctx, ~tab->row_var[row] == tab->n_con - 1, return -1);
1744         if (row != tab->n_row - 1)
1745                 swap_rows(tab, row, tab->n_row - 1);
1746         tab->n_row--;
1747         tab->n_con--;
1748         return 0;
1749 }
1750
1751 static int drop_col(struct isl_tab *tab, int col)
1752 {
1753         isl_assert(tab->mat->ctx, tab->col_var[col] == tab->n_var - 1, return -1);
1754         if (col != tab->n_col - 1)
1755                 swap_cols(tab, col, tab->n_col - 1);
1756         tab->n_col--;
1757         tab->n_var--;
1758         return 0;
1759 }
1760
1761 /* Add inequality "ineq" and check if it conflicts with the
1762  * previously added constraints or if it is obviously redundant.
1763  */
1764 int isl_tab_add_ineq(struct isl_tab *tab, isl_int *ineq)
1765 {
1766         int r;
1767         int sgn;
1768         isl_int cst;
1769
1770         if (!tab)
1771                 return -1;
1772         if (tab->bmap) {
1773                 struct isl_basic_map *bmap = tab->bmap;
1774
1775                 isl_assert(tab->mat->ctx, tab->n_eq == bmap->n_eq, return -1);
1776                 isl_assert(tab->mat->ctx,
1777                             tab->n_con == bmap->n_eq + bmap->n_ineq, return -1);
1778                 tab->bmap = isl_basic_map_add_ineq(tab->bmap, ineq);
1779                 if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1780                         return -1;
1781                 if (!tab->bmap)
1782                         return -1;
1783         }
1784         if (tab->cone) {
1785                 isl_int_init(cst);
1786                 isl_int_swap(ineq[0], cst);
1787         }
1788         r = isl_tab_add_row(tab, ineq);
1789         if (tab->cone) {
1790                 isl_int_swap(ineq[0], cst);
1791                 isl_int_clear(cst);
1792         }
1793         if (r < 0)
1794                 return -1;
1795         tab->con[r].is_nonneg = 1;
1796         if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1797                 return -1;
1798         if (isl_tab_row_is_redundant(tab, tab->con[r].index)) {
1799                 if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0)
1800                         return -1;
1801                 return 0;
1802         }
1803
1804         sgn = restore_row(tab, &tab->con[r]);
1805         if (sgn < -1)
1806                 return -1;
1807         if (sgn < 0)
1808                 return isl_tab_mark_empty(tab);
1809         if (tab->con[r].is_row && isl_tab_row_is_redundant(tab, tab->con[r].index))
1810                 if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0)
1811                         return -1;
1812         return 0;
1813 }
1814
1815 /* Pivot a non-negative variable down until it reaches the value zero
1816  * and then pivot the variable into a column position.
1817  */
1818 static int to_col(struct isl_tab *tab, struct isl_tab_var *var) WARN_UNUSED;
1819 static int to_col(struct isl_tab *tab, struct isl_tab_var *var)
1820 {
1821         int i;
1822         int row, col;
1823         unsigned off = 2 + tab->M;
1824
1825         if (!var->is_row)
1826                 return 0;
1827
1828         while (isl_int_is_pos(tab->mat->row[var->index][1])) {
1829                 find_pivot(tab, var, NULL, -1, &row, &col);
1830                 isl_assert(tab->mat->ctx, row != -1, return -1);
1831                 if (isl_tab_pivot(tab, row, col) < 0)
1832                         return -1;
1833                 if (!var->is_row)
1834                         return 0;
1835         }
1836
1837         for (i = tab->n_dead; i < tab->n_col; ++i)
1838                 if (!isl_int_is_zero(tab->mat->row[var->index][off + i]))
1839                         break;
1840
1841         isl_assert(tab->mat->ctx, i < tab->n_col, return -1);
1842         if (isl_tab_pivot(tab, var->index, i) < 0)
1843                 return -1;
1844
1845         return 0;
1846 }
1847
1848 /* We assume Gaussian elimination has been performed on the equalities.
1849  * The equalities can therefore never conflict.
1850  * Adding the equalities is currently only really useful for a later call
1851  * to isl_tab_ineq_type.
1852  */
1853 static struct isl_tab *add_eq(struct isl_tab *tab, isl_int *eq)
1854 {
1855         int i;
1856         int r;
1857
1858         if (!tab)
1859                 return NULL;
1860         r = isl_tab_add_row(tab, eq);
1861         if (r < 0)
1862                 goto error;
1863
1864         r = tab->con[r].index;
1865         i = isl_seq_first_non_zero(tab->mat->row[r] + 2 + tab->M + tab->n_dead,
1866                                         tab->n_col - tab->n_dead);
1867         isl_assert(tab->mat->ctx, i >= 0, goto error);
1868         i += tab->n_dead;
1869         if (isl_tab_pivot(tab, r, i) < 0)
1870                 goto error;
1871         if (isl_tab_kill_col(tab, i) < 0)
1872                 goto error;
1873         tab->n_eq++;
1874
1875         return tab;
1876 error:
1877         isl_tab_free(tab);
1878         return NULL;
1879 }
1880
1881 static int row_is_manifestly_zero(struct isl_tab *tab, int row)
1882 {
1883         unsigned off = 2 + tab->M;
1884
1885         if (!isl_int_is_zero(tab->mat->row[row][1]))
1886                 return 0;
1887         if (tab->M && !isl_int_is_zero(tab->mat->row[row][2]))
1888                 return 0;
1889         return isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
1890                                         tab->n_col - tab->n_dead) == -1;
1891 }
1892
1893 /* Add an equality that is known to be valid for the given tableau.
1894  */
1895 int isl_tab_add_valid_eq(struct isl_tab *tab, isl_int *eq)
1896 {
1897         struct isl_tab_var *var;
1898         int r;
1899
1900         if (!tab)
1901                 return -1;
1902         r = isl_tab_add_row(tab, eq);
1903         if (r < 0)
1904                 return -1;
1905
1906         var = &tab->con[r];
1907         r = var->index;
1908         if (row_is_manifestly_zero(tab, r)) {
1909                 var->is_zero = 1;
1910                 if (isl_tab_mark_redundant(tab, r) < 0)
1911                         return -1;
1912                 return 0;
1913         }
1914
1915         if (isl_int_is_neg(tab->mat->row[r][1])) {
1916                 isl_seq_neg(tab->mat->row[r] + 1, tab->mat->row[r] + 1,
1917                             1 + tab->n_col);
1918                 var->negated = 1;
1919         }
1920         var->is_nonneg = 1;
1921         if (to_col(tab, var) < 0)
1922                 return -1;
1923         var->is_nonneg = 0;
1924         if (isl_tab_kill_col(tab, var->index) < 0)
1925                 return -1;
1926
1927         return 0;
1928 }
1929
1930 static int add_zero_row(struct isl_tab *tab)
1931 {
1932         int r;
1933         isl_int *row;
1934
1935         r = isl_tab_allocate_con(tab);
1936         if (r < 0)
1937                 return -1;
1938
1939         row = tab->mat->row[tab->con[r].index];
1940         isl_seq_clr(row + 1, 1 + tab->M + tab->n_col);
1941         isl_int_set_si(row[0], 1);
1942
1943         return r;
1944 }
1945
1946 /* Add equality "eq" and check if it conflicts with the
1947  * previously added constraints or if it is obviously redundant.
1948  */
1949 int isl_tab_add_eq(struct isl_tab *tab, isl_int *eq)
1950 {
1951         struct isl_tab_undo *snap = NULL;
1952         struct isl_tab_var *var;
1953         int r;
1954         int row;
1955         int sgn;
1956         isl_int cst;
1957
1958         if (!tab)
1959                 return -1;
1960         isl_assert(tab->mat->ctx, !tab->M, return -1);
1961
1962         if (tab->need_undo)
1963                 snap = isl_tab_snap(tab);
1964
1965         if (tab->cone) {
1966                 isl_int_init(cst);
1967                 isl_int_swap(eq[0], cst);
1968         }
1969         r = isl_tab_add_row(tab, eq);
1970         if (tab->cone) {
1971                 isl_int_swap(eq[0], cst);
1972                 isl_int_clear(cst);
1973         }
1974         if (r < 0)
1975                 return -1;
1976
1977         var = &tab->con[r];
1978         row = var->index;
1979         if (row_is_manifestly_zero(tab, row)) {
1980                 if (snap) {
1981                         if (isl_tab_rollback(tab, snap) < 0)
1982                                 return -1;
1983                 } else
1984                         drop_row(tab, row);
1985                 return 0;
1986         }
1987
1988         if (tab->bmap) {
1989                 tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);
1990                 if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1991                         return -1;
1992                 isl_seq_neg(eq, eq, 1 + tab->n_var);
1993                 tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);
1994                 isl_seq_neg(eq, eq, 1 + tab->n_var);
1995                 if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1996                         return -1;
1997                 if (!tab->bmap)
1998                         return -1;
1999                 if (add_zero_row(tab) < 0)
2000                         return -1;
2001         }
2002
2003         sgn = isl_int_sgn(tab->mat->row[row][1]);
2004
2005         if (sgn > 0) {
2006                 isl_seq_neg(tab->mat->row[row] + 1, tab->mat->row[row] + 1,
2007                             1 + tab->n_col);
2008                 var->negated = 1;
2009                 sgn = -1;
2010         }
2011
2012         if (sgn < 0) {
2013                 sgn = sign_of_max(tab, var);
2014                 if (sgn < -1)
2015                         return -1;
2016                 if (sgn < 0) {
2017                         if (isl_tab_mark_empty(tab) < 0)
2018                                 return -1;
2019                         return 0;
2020                 }
2021         }
2022
2023         var->is_nonneg = 1;
2024         if (to_col(tab, var) < 0)
2025                 return -1;
2026         var->is_nonneg = 0;
2027         if (isl_tab_kill_col(tab, var->index) < 0)
2028                 return -1;
2029
2030         return 0;
2031 }
2032
2033 /* Construct and return an inequality that expresses an upper bound
2034  * on the given div.
2035  * In particular, if the div is given by
2036  *
2037  *      d = floor(e/m)
2038  *
2039  * then the inequality expresses
2040  *
2041  *      m d <= e
2042  */
2043 static struct isl_vec *ineq_for_div(struct isl_basic_map *bmap, unsigned div)
2044 {
2045         unsigned total;
2046         unsigned div_pos;
2047         struct isl_vec *ineq;
2048
2049         if (!bmap)
2050                 return NULL;
2051
2052         total = isl_basic_map_total_dim(bmap);
2053         div_pos = 1 + total - bmap->n_div + div;
2054
2055         ineq = isl_vec_alloc(bmap->ctx, 1 + total);
2056         if (!ineq)
2057                 return NULL;
2058
2059         isl_seq_cpy(ineq->el, bmap->div[div] + 1, 1 + total);
2060         isl_int_neg(ineq->el[div_pos], bmap->div[div][0]);
2061         return ineq;
2062 }
2063
2064 /* For a div d = floor(f/m), add the constraints
2065  *
2066  *              f - m d >= 0
2067  *              -(f-(m-1)) + m d >= 0
2068  *
2069  * Note that the second constraint is the negation of
2070  *
2071  *              f - m d >= m
2072  *
2073  * If add_ineq is not NULL, then this function is used
2074  * instead of isl_tab_add_ineq to effectively add the inequalities.
2075  */
2076 static int add_div_constraints(struct isl_tab *tab, unsigned div,
2077         int (*add_ineq)(void *user, isl_int *), void *user)
2078 {
2079         unsigned total;
2080         unsigned div_pos;
2081         struct isl_vec *ineq;
2082
2083         total = isl_basic_map_total_dim(tab->bmap);
2084         div_pos = 1 + total - tab->bmap->n_div + div;
2085
2086         ineq = ineq_for_div(tab->bmap, div);
2087         if (!ineq)
2088                 goto error;
2089
2090         if (add_ineq) {
2091                 if (add_ineq(user, ineq->el) < 0)
2092                         goto error;
2093         } else {
2094                 if (isl_tab_add_ineq(tab, ineq->el) < 0)
2095                         goto error;
2096         }
2097
2098         isl_seq_neg(ineq->el, tab->bmap->div[div] + 1, 1 + total);
2099         isl_int_set(ineq->el[div_pos], tab->bmap->div[div][0]);
2100         isl_int_add(ineq->el[0], ineq->el[0], ineq->el[div_pos]);
2101         isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
2102
2103         if (add_ineq) {
2104                 if (add_ineq(user, ineq->el) < 0)
2105                         goto error;
2106         } else {
2107                 if (isl_tab_add_ineq(tab, ineq->el) < 0)
2108                         goto error;
2109         }
2110
2111         isl_vec_free(ineq);
2112
2113         return 0;
2114 error:
2115         isl_vec_free(ineq);
2116         return -1;
2117 }
2118
2119 /* Check whether the div described by "div" is obviously non-negative.
2120  * If we are using a big parameter, then we will encode the div
2121  * as div' = M + div, which is always non-negative.
2122  * Otherwise, we check whether div is a non-negative affine combination
2123  * of non-negative variables.
2124  */
2125 static int div_is_nonneg(struct isl_tab *tab, __isl_keep isl_vec *div)
2126 {
2127         int i;
2128
2129         if (tab->M)
2130                 return 1;
2131
2132         if (isl_int_is_neg(div->el[1]))
2133                 return 0;
2134
2135         for (i = 0; i < tab->n_var; ++i) {
2136                 if (isl_int_is_neg(div->el[2 + i]))
2137                         return 0;
2138                 if (isl_int_is_zero(div->el[2 + i]))
2139                         continue;
2140                 if (!tab->var[i].is_nonneg)
2141                         return 0;
2142         }
2143
2144         return 1;
2145 }
2146
2147 /* Add an extra div, prescribed by "div" to the tableau and
2148  * the associated bmap (which is assumed to be non-NULL).
2149  *
2150  * If add_ineq is not NULL, then this function is used instead
2151  * of isl_tab_add_ineq to add the div constraints.
2152  * This complication is needed because the code in isl_tab_pip
2153  * wants to perform some extra processing when an inequality
2154  * is added to the tableau.
2155  */
2156 int isl_tab_add_div(struct isl_tab *tab, __isl_keep isl_vec *div,
2157         int (*add_ineq)(void *user, isl_int *), void *user)
2158 {
2159         int r;
2160         int k;
2161         int nonneg;
2162
2163         if (!tab || !div)
2164                 return -1;
2165
2166         isl_assert(tab->mat->ctx, tab->bmap, return -1);
2167
2168         nonneg = div_is_nonneg(tab, div);
2169
2170         if (isl_tab_extend_cons(tab, 3) < 0)
2171                 return -1;
2172         if (isl_tab_extend_vars(tab, 1) < 0)
2173                 return -1;
2174         r = isl_tab_allocate_var(tab);
2175         if (r < 0)
2176                 return -1;
2177
2178         if (nonneg)
2179                 tab->var[r].is_nonneg = 1;
2180
2181         tab->bmap = isl_basic_map_extend_dim(tab->bmap,
2182                 isl_basic_map_get_dim(tab->bmap), 1, 0, 2);
2183         k = isl_basic_map_alloc_div(tab->bmap);
2184         if (k < 0)
2185                 return -1;
2186         isl_seq_cpy(tab->bmap->div[k], div->el, div->size);
2187         if (isl_tab_push(tab, isl_tab_undo_bmap_div) < 0)
2188                 return -1;
2189
2190         if (add_div_constraints(tab, k, add_ineq, user) < 0)
2191                 return -1;
2192
2193         return r;
2194 }
2195
2196 struct isl_tab *isl_tab_from_basic_map(struct isl_basic_map *bmap)
2197 {
2198         int i;
2199         struct isl_tab *tab;
2200
2201         if (!bmap)
2202                 return NULL;
2203         tab = isl_tab_alloc(bmap->ctx,
2204                             isl_basic_map_total_dim(bmap) + bmap->n_ineq + 1,
2205                             isl_basic_map_total_dim(bmap), 0);
2206         if (!tab)
2207                 return NULL;
2208         tab->rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
2209         if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY)) {
2210                 if (isl_tab_mark_empty(tab) < 0)
2211                         goto error;
2212                 return tab;
2213         }
2214         for (i = 0; i < bmap->n_eq; ++i) {
2215                 tab = add_eq(tab, bmap->eq[i]);
2216                 if (!tab)
2217                         return tab;
2218         }
2219         for (i = 0; i < bmap->n_ineq; ++i) {
2220                 if (isl_tab_add_ineq(tab, bmap->ineq[i]) < 0)
2221                         goto error;
2222                 if (tab->empty)
2223                         return tab;
2224         }
2225         return tab;
2226 error:
2227         isl_tab_free(tab);
2228         return NULL;
2229 }
2230
2231 struct isl_tab *isl_tab_from_basic_set(struct isl_basic_set *bset)
2232 {
2233         return isl_tab_from_basic_map((struct isl_basic_map *)bset);
2234 }
2235
2236 /* Construct a tableau corresponding to the recession cone of "bset".
2237  */
2238 struct isl_tab *isl_tab_from_recession_cone(__isl_keep isl_basic_set *bset,
2239         int parametric)
2240 {
2241         isl_int cst;
2242         int i;
2243         struct isl_tab *tab;
2244         unsigned offset = 0;
2245
2246         if (!bset)
2247                 return NULL;
2248         if (parametric)
2249                 offset = isl_basic_set_dim(bset, isl_dim_param);
2250         tab = isl_tab_alloc(bset->ctx, bset->n_eq + bset->n_ineq,
2251                                 isl_basic_set_total_dim(bset) - offset, 0);
2252         if (!tab)
2253                 return NULL;
2254         tab->rational = ISL_F_ISSET(bset, ISL_BASIC_SET_RATIONAL);
2255         tab->cone = 1;
2256
2257         isl_int_init(cst);
2258         for (i = 0; i < bset->n_eq; ++i) {
2259                 isl_int_swap(bset->eq[i][offset], cst);
2260                 if (offset > 0) {
2261                         if (isl_tab_add_eq(tab, bset->eq[i] + offset) < 0)
2262                                 goto error;
2263                 } else
2264                         tab = add_eq(tab, bset->eq[i]);
2265                 isl_int_swap(bset->eq[i][offset], cst);
2266                 if (!tab)
2267                         goto done;
2268         }
2269         for (i = 0; i < bset->n_ineq; ++i) {
2270                 int r;
2271                 isl_int_swap(bset->ineq[i][offset], cst);
2272                 r = isl_tab_add_row(tab, bset->ineq[i] + offset);
2273                 isl_int_swap(bset->ineq[i][offset], cst);
2274                 if (r < 0)
2275                         goto error;
2276                 tab->con[r].is_nonneg = 1;
2277                 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
2278                         goto error;
2279         }
2280 done:
2281         isl_int_clear(cst);
2282         return tab;
2283 error:
2284         isl_int_clear(cst);
2285         isl_tab_free(tab);
2286         return NULL;
2287 }
2288
2289 /* Assuming "tab" is the tableau of a cone, check if the cone is
2290  * bounded, i.e., if it is empty or only contains the origin.
2291  */
2292 int isl_tab_cone_is_bounded(struct isl_tab *tab)
2293 {
2294         int i;
2295
2296         if (!tab)
2297                 return -1;
2298         if (tab->empty)
2299                 return 1;
2300         if (tab->n_dead == tab->n_col)
2301                 return 1;
2302
2303         for (;;) {
2304                 for (i = tab->n_redundant; i < tab->n_row; ++i) {
2305                         struct isl_tab_var *var;
2306                         int sgn;
2307                         var = isl_tab_var_from_row(tab, i);
2308                         if (!var->is_nonneg)
2309                                 continue;
2310                         sgn = sign_of_max(tab, var);
2311                         if (sgn < -1)
2312                                 return -1;
2313                         if (sgn != 0)
2314                                 return 0;
2315                         if (close_row(tab, var) < 0)
2316                                 return -1;
2317                         break;
2318                 }
2319                 if (tab->n_dead == tab->n_col)
2320                         return 1;
2321                 if (i == tab->n_row)
2322                         return 0;
2323         }
2324 }
2325
2326 int isl_tab_sample_is_integer(struct isl_tab *tab)
2327 {
2328         int i;
2329
2330         if (!tab)
2331                 return -1;
2332
2333         for (i = 0; i < tab->n_var; ++i) {
2334                 int row;
2335                 if (!tab->var[i].is_row)
2336                         continue;
2337                 row = tab->var[i].index;
2338                 if (!isl_int_is_divisible_by(tab->mat->row[row][1],
2339                                                 tab->mat->row[row][0]))
2340                         return 0;
2341         }
2342         return 1;
2343 }
2344
2345 static struct isl_vec *extract_integer_sample(struct isl_tab *tab)
2346 {
2347         int i;
2348         struct isl_vec *vec;
2349
2350         vec = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
2351         if (!vec)
2352                 return NULL;
2353
2354         isl_int_set_si(vec->block.data[0], 1);
2355         for (i = 0; i < tab->n_var; ++i) {
2356                 if (!tab->var[i].is_row)
2357                         isl_int_set_si(vec->block.data[1 + i], 0);
2358                 else {
2359                         int row = tab->var[i].index;
2360                         isl_int_divexact(vec->block.data[1 + i],
2361                                 tab->mat->row[row][1], tab->mat->row[row][0]);
2362                 }
2363         }
2364
2365         return vec;
2366 }
2367
2368 struct isl_vec *isl_tab_get_sample_value(struct isl_tab *tab)
2369 {
2370         int i;
2371         struct isl_vec *vec;
2372         isl_int m;
2373
2374         if (!tab)
2375                 return NULL;
2376
2377         vec = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
2378         if (!vec)
2379                 return NULL;
2380
2381         isl_int_init(m);
2382
2383         isl_int_set_si(vec->block.data[0], 1);
2384         for (i = 0; i < tab->n_var; ++i) {
2385                 int row;
2386                 if (!tab->var[i].is_row) {
2387                         isl_int_set_si(vec->block.data[1 + i], 0);
2388                         continue;
2389                 }
2390                 row = tab->var[i].index;
2391                 isl_int_gcd(m, vec->block.data[0], tab->mat->row[row][0]);
2392                 isl_int_divexact(m, tab->mat->row[row][0], m);
2393                 isl_seq_scale(vec->block.data, vec->block.data, m, 1 + i);
2394                 isl_int_divexact(m, vec->block.data[0], tab->mat->row[row][0]);
2395                 isl_int_mul(vec->block.data[1 + i], m, tab->mat->row[row][1]);
2396         }
2397         vec = isl_vec_normalize(vec);
2398
2399         isl_int_clear(m);
2400         return vec;
2401 }
2402
2403 /* Update "bmap" based on the results of the tableau "tab".
2404  * In particular, implicit equalities are made explicit, redundant constraints
2405  * are removed and if the sample value happens to be integer, it is stored
2406  * in "bmap" (unless "bmap" already had an integer sample).
2407  *
2408  * The tableau is assumed to have been created from "bmap" using
2409  * isl_tab_from_basic_map.
2410  */
2411 struct isl_basic_map *isl_basic_map_update_from_tab(struct isl_basic_map *bmap,
2412         struct isl_tab *tab)
2413 {
2414         int i;
2415         unsigned n_eq;
2416
2417         if (!bmap)
2418                 return NULL;
2419         if (!tab)
2420                 return bmap;
2421
2422         n_eq = tab->n_eq;
2423         if (tab->empty)
2424                 bmap = isl_basic_map_set_to_empty(bmap);
2425         else
2426                 for (i = bmap->n_ineq - 1; i >= 0; --i) {
2427                         if (isl_tab_is_equality(tab, n_eq + i))
2428                                 isl_basic_map_inequality_to_equality(bmap, i);
2429                         else if (isl_tab_is_redundant(tab, n_eq + i))
2430                                 isl_basic_map_drop_inequality(bmap, i);
2431                 }
2432         if (bmap->n_eq != n_eq)
2433                 isl_basic_map_gauss(bmap, NULL);
2434         if (!tab->rational &&
2435             !bmap->sample && isl_tab_sample_is_integer(tab))
2436                 bmap->sample = extract_integer_sample(tab);
2437         return bmap;
2438 }
2439
2440 struct isl_basic_set *isl_basic_set_update_from_tab(struct isl_basic_set *bset,
2441         struct isl_tab *tab)
2442 {
2443         return (struct isl_basic_set *)isl_basic_map_update_from_tab(
2444                 (struct isl_basic_map *)bset, tab);
2445 }
2446
2447 /* Given a non-negative variable "var", add a new non-negative variable
2448  * that is the opposite of "var", ensuring that var can only attain the
2449  * value zero.
2450  * If var = n/d is a row variable, then the new variable = -n/d.
2451  * If var is a column variables, then the new variable = -var.
2452  * If the new variable cannot attain non-negative values, then
2453  * the resulting tableau is empty.
2454  * Otherwise, we know the value will be zero and we close the row.
2455  */
2456 static int cut_to_hyperplane(struct isl_tab *tab, struct isl_tab_var *var)
2457 {
2458         unsigned r;
2459         isl_int *row;
2460         int sgn;
2461         unsigned off = 2 + tab->M;
2462
2463         if (var->is_zero)
2464                 return 0;
2465         isl_assert(tab->mat->ctx, !var->is_redundant, return -1);
2466         isl_assert(tab->mat->ctx, var->is_nonneg, return -1);
2467
2468         if (isl_tab_extend_cons(tab, 1) < 0)
2469                 return -1;
2470
2471         r = tab->n_con;
2472         tab->con[r].index = tab->n_row;
2473         tab->con[r].is_row = 1;
2474         tab->con[r].is_nonneg = 0;
2475         tab->con[r].is_zero = 0;
2476         tab->con[r].is_redundant = 0;
2477         tab->con[r].frozen = 0;
2478         tab->con[r].negated = 0;
2479         tab->row_var[tab->n_row] = ~r;
2480         row = tab->mat->row[tab->n_row];
2481
2482         if (var->is_row) {
2483                 isl_int_set(row[0], tab->mat->row[var->index][0]);
2484                 isl_seq_neg(row + 1,
2485                             tab->mat->row[var->index] + 1, 1 + tab->n_col);
2486         } else {
2487                 isl_int_set_si(row[0], 1);
2488                 isl_seq_clr(row + 1, 1 + tab->n_col);
2489                 isl_int_set_si(row[off + var->index], -1);
2490         }
2491
2492         tab->n_row++;
2493         tab->n_con++;
2494         if (isl_tab_push_var(tab, isl_tab_undo_allocate, &tab->con[r]) < 0)
2495                 return -1;
2496
2497         sgn = sign_of_max(tab, &tab->con[r]);
2498         if (sgn < -1)
2499                 return -1;
2500         if (sgn < 0) {
2501                 if (isl_tab_mark_empty(tab) < 0)
2502                         return -1;
2503                 return 0;
2504         }
2505         tab->con[r].is_nonneg = 1;
2506         if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
2507                 return -1;
2508         /* sgn == 0 */
2509         if (close_row(tab, &tab->con[r]) < 0)
2510                 return -1;
2511
2512         return 0;
2513 }
2514
2515 /* Given a tableau "tab" and an inequality constraint "con" of the tableau,
2516  * relax the inequality by one.  That is, the inequality r >= 0 is replaced
2517  * by r' = r + 1 >= 0.
2518  * If r is a row variable, we simply increase the constant term by one
2519  * (taking into account the denominator).
2520  * If r is a column variable, then we need to modify each row that
2521  * refers to r = r' - 1 by substituting this equality, effectively
2522  * subtracting the coefficient of the column from the constant.
2523  * We should only do this if the minimum is manifestly unbounded,
2524  * however.  Otherwise, we may end up with negative sample values
2525  * for non-negative variables.
2526  * So, if r is a column variable with a minimum that is not
2527  * manifestly unbounded, then we need to move it to a row.
2528  * However, the sample value of this row may be negative,
2529  * even after the relaxation, so we need to restore it.
2530  * We therefore prefer to pivot a column up to a row, if possible.
2531  */
2532 struct isl_tab *isl_tab_relax(struct isl_tab *tab, int con)
2533 {
2534         struct isl_tab_var *var;
2535         unsigned off = 2 + tab->M;
2536
2537         if (!tab)
2538                 return NULL;
2539
2540         var = &tab->con[con];
2541
2542         if (!var->is_row && !max_is_manifestly_unbounded(tab, var))
2543                 if (to_row(tab, var, 1) < 0)
2544                         goto error;
2545         if (!var->is_row && !min_is_manifestly_unbounded(tab, var))
2546                 if (to_row(tab, var, -1) < 0)
2547                         goto error;
2548
2549         if (var->is_row) {
2550                 isl_int_add(tab->mat->row[var->index][1],
2551                     tab->mat->row[var->index][1], tab->mat->row[var->index][0]);
2552                 if (restore_row(tab, var) < 0)
2553                         goto error;
2554         } else {
2555                 int i;
2556
2557                 for (i = 0; i < tab->n_row; ++i) {
2558                         if (isl_int_is_zero(tab->mat->row[i][off + var->index]))
2559                                 continue;
2560                         isl_int_sub(tab->mat->row[i][1], tab->mat->row[i][1],
2561                             tab->mat->row[i][off + var->index]);
2562                 }
2563
2564         }
2565
2566         if (isl_tab_push_var(tab, isl_tab_undo_relax, var) < 0)
2567                 goto error;
2568
2569         return tab;
2570 error:
2571         isl_tab_free(tab);
2572         return NULL;
2573 }
2574
2575 int isl_tab_select_facet(struct isl_tab *tab, int con)
2576 {
2577         if (!tab)
2578                 return -1;
2579
2580         return cut_to_hyperplane(tab, &tab->con[con]);
2581 }
2582
2583 static int may_be_equality(struct isl_tab *tab, int row)
2584 {
2585         unsigned off = 2 + tab->M;
2586         return tab->rational ? isl_int_is_zero(tab->mat->row[row][1])
2587                              : isl_int_lt(tab->mat->row[row][1],
2588                                             tab->mat->row[row][0]);
2589 }
2590
2591 /* Check for (near) equalities among the constraints.
2592  * A constraint is an equality if it is non-negative and if
2593  * its maximal value is either
2594  *      - zero (in case of rational tableaus), or
2595  *      - strictly less than 1 (in case of integer tableaus)
2596  *
2597  * We first mark all non-redundant and non-dead variables that
2598  * are not frozen and not obviously not an equality.
2599  * Then we iterate over all marked variables if they can attain
2600  * any values larger than zero or at least one.
2601  * If the maximal value is zero, we mark any column variables
2602  * that appear in the row as being zero and mark the row as being redundant.
2603  * Otherwise, if the maximal value is strictly less than one (and the
2604  * tableau is integer), then we restrict the value to being zero
2605  * by adding an opposite non-negative variable.
2606  */
2607 int isl_tab_detect_implicit_equalities(struct isl_tab *tab)
2608 {
2609         int i;
2610         unsigned n_marked;
2611
2612         if (!tab)
2613                 return -1;
2614         if (tab->empty)
2615                 return 0;
2616         if (tab->n_dead == tab->n_col)
2617                 return 0;
2618
2619         n_marked = 0;
2620         for (i = tab->n_redundant; i < tab->n_row; ++i) {
2621                 struct isl_tab_var *var = isl_tab_var_from_row(tab, i);
2622                 var->marked = !var->frozen && var->is_nonneg &&
2623                         may_be_equality(tab, i);
2624                 if (var->marked)
2625                         n_marked++;
2626         }
2627         for (i = tab->n_dead; i < tab->n_col; ++i) {
2628                 struct isl_tab_var *var = var_from_col(tab, i);
2629                 var->marked = !var->frozen && var->is_nonneg;
2630                 if (var->marked)
2631                         n_marked++;
2632         }
2633         while (n_marked) {
2634                 struct isl_tab_var *var;
2635                 int sgn;
2636                 for (i = tab->n_redundant; i < tab->n_row; ++i) {
2637                         var = isl_tab_var_from_row(tab, i);
2638                         if (var->marked)
2639                                 break;
2640                 }
2641                 if (i == tab->n_row) {
2642                         for (i = tab->n_dead; i < tab->n_col; ++i) {
2643                                 var = var_from_col(tab, i);
2644                                 if (var->marked)
2645                                         break;
2646                         }
2647                         if (i == tab->n_col)
2648                                 break;
2649                 }
2650                 var->marked = 0;
2651                 n_marked--;
2652                 sgn = sign_of_max(tab, var);
2653                 if (sgn < 0)
2654                         return -1;
2655                 if (sgn == 0) {
2656                         if (close_row(tab, var) < 0)
2657                                 return -1;
2658                 } else if (!tab->rational && !at_least_one(tab, var)) {
2659                         if (cut_to_hyperplane(tab, var) < 0)
2660                                 return -1;
2661                         return isl_tab_detect_implicit_equalities(tab);
2662                 }
2663                 for (i = tab->n_redundant; i < tab->n_row; ++i) {
2664                         var = isl_tab_var_from_row(tab, i);
2665                         if (!var->marked)
2666                                 continue;
2667                         if (may_be_equality(tab, i))
2668                                 continue;
2669                         var->marked = 0;
2670                         n_marked--;
2671                 }
2672         }
2673
2674         return 0;
2675 }
2676
2677 static int con_is_redundant(struct isl_tab *tab, struct isl_tab_var *var)
2678 {
2679         if (!tab)
2680                 return -1;
2681         if (tab->rational) {
2682                 int sgn = sign_of_min(tab, var);
2683                 if (sgn < -1)
2684                         return -1;
2685                 return sgn >= 0;
2686         } else {
2687                 int irred = isl_tab_min_at_most_neg_one(tab, var);
2688                 if (irred < 0)
2689                         return -1;
2690                 return !irred;
2691         }
2692 }
2693
2694 /* Check for (near) redundant constraints.
2695  * A constraint is redundant if it is non-negative and if
2696  * its minimal value (temporarily ignoring the non-negativity) is either
2697  *      - zero (in case of rational tableaus), or
2698  *      - strictly larger than -1 (in case of integer tableaus)
2699  *
2700  * We first mark all non-redundant and non-dead variables that
2701  * are not frozen and not obviously negatively unbounded.
2702  * Then we iterate over all marked variables if they can attain
2703  * any values smaller than zero or at most negative one.
2704  * If not, we mark the row as being redundant (assuming it hasn't
2705  * been detected as being obviously redundant in the mean time).
2706  */
2707 int isl_tab_detect_redundant(struct isl_tab *tab)
2708 {
2709         int i;
2710         unsigned n_marked;
2711
2712         if (!tab)
2713                 return -1;
2714         if (tab->empty)
2715                 return 0;
2716         if (tab->n_redundant == tab->n_row)
2717                 return 0;
2718
2719         n_marked = 0;
2720         for (i = tab->n_redundant; i < tab->n_row; ++i) {
2721                 struct isl_tab_var *var = isl_tab_var_from_row(tab, i);
2722                 var->marked = !var->frozen && var->is_nonneg;
2723                 if (var->marked)
2724                         n_marked++;
2725         }
2726         for (i = tab->n_dead; i < tab->n_col; ++i) {
2727                 struct isl_tab_var *var = var_from_col(tab, i);
2728                 var->marked = !var->frozen && var->is_nonneg &&
2729                         !min_is_manifestly_unbounded(tab, var);
2730                 if (var->marked)
2731                         n_marked++;
2732         }
2733         while (n_marked) {
2734                 struct isl_tab_var *var;
2735                 int red;
2736                 for (i = tab->n_redundant; i < tab->n_row; ++i) {
2737                         var = isl_tab_var_from_row(tab, i);
2738                         if (var->marked)
2739                                 break;
2740                 }
2741                 if (i == tab->n_row) {
2742                         for (i = tab->n_dead; i < tab->n_col; ++i) {
2743                                 var = var_from_col(tab, i);
2744                                 if (var->marked)
2745                                         break;
2746                         }
2747                         if (i == tab->n_col)
2748                                 break;
2749                 }
2750                 var->marked = 0;
2751                 n_marked--;
2752                 red = con_is_redundant(tab, var);
2753                 if (red < 0)
2754                         return -1;
2755                 if (red && !var->is_redundant)
2756                         if (isl_tab_mark_redundant(tab, var->index) < 0)
2757                                 return -1;
2758                 for (i = tab->n_dead; i < tab->n_col; ++i) {
2759                         var = var_from_col(tab, i);
2760                         if (!var->marked)
2761                                 continue;
2762                         if (!min_is_manifestly_unbounded(tab, var))
2763                                 continue;
2764                         var->marked = 0;
2765                         n_marked--;
2766                 }
2767         }
2768
2769         return 0;
2770 }
2771
2772 int isl_tab_is_equality(struct isl_tab *tab, int con)
2773 {
2774         int row;
2775         unsigned off;
2776
2777         if (!tab)
2778                 return -1;
2779         if (tab->con[con].is_zero)
2780                 return 1;
2781         if (tab->con[con].is_redundant)
2782                 return 0;
2783         if (!tab->con[con].is_row)
2784                 return tab->con[con].index < tab->n_dead;
2785
2786         row = tab->con[con].index;
2787
2788         off = 2 + tab->M;
2789         return isl_int_is_zero(tab->mat->row[row][1]) &&
2790                 (!tab->M || isl_int_is_zero(tab->mat->row[row][2])) &&
2791                 isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
2792                                         tab->n_col - tab->n_dead) == -1;
2793 }
2794
2795 /* Return the minimal value of the affine expression "f" with denominator
2796  * "denom" in *opt, *opt_denom, assuming the tableau is not empty and
2797  * the expression cannot attain arbitrarily small values.
2798  * If opt_denom is NULL, then *opt is rounded up to the nearest integer.
2799  * The return value reflects the nature of the result (empty, unbounded,
2800  * minimal value returned in *opt).
2801  */
2802 enum isl_lp_result isl_tab_min(struct isl_tab *tab,
2803         isl_int *f, isl_int denom, isl_int *opt, isl_int *opt_denom,
2804         unsigned flags)
2805 {
2806         int r;
2807         enum isl_lp_result res = isl_lp_ok;
2808         struct isl_tab_var *var;
2809         struct isl_tab_undo *snap;
2810
2811         if (!tab)
2812                 return isl_lp_error;
2813
2814         if (tab->empty)
2815                 return isl_lp_empty;
2816
2817         snap = isl_tab_snap(tab);
2818         r = isl_tab_add_row(tab, f);
2819         if (r < 0)
2820                 return isl_lp_error;
2821         var = &tab->con[r];
2822         for (;;) {
2823                 int row, col;
2824                 find_pivot(tab, var, var, -1, &row, &col);
2825                 if (row == var->index) {
2826                         res = isl_lp_unbounded;
2827                         break;
2828                 }
2829                 if (row == -1)
2830                         break;
2831                 if (isl_tab_pivot(tab, row, col) < 0)
2832                         return isl_lp_error;
2833         }
2834         isl_int_mul(tab->mat->row[var->index][0],
2835                     tab->mat->row[var->index][0], denom);
2836         if (ISL_FL_ISSET(flags, ISL_TAB_SAVE_DUAL)) {
2837                 int i;
2838
2839                 isl_vec_free(tab->dual);
2840                 tab->dual = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_con);
2841                 if (!tab->dual)
2842                         return isl_lp_error;
2843                 isl_int_set(tab->dual->el[0], tab->mat->row[var->index][0]);
2844                 for (i = 0; i < tab->n_con; ++i) {
2845                         int pos;
2846                         if (tab->con[i].is_row) {
2847                                 isl_int_set_si(tab->dual->el[1 + i], 0);
2848                                 continue;
2849                         }
2850                         pos = 2 + tab->M + tab->con[i].index;
2851                         if (tab->con[i].negated)
2852                                 isl_int_neg(tab->dual->el[1 + i],
2853                                             tab->mat->row[var->index][pos]);
2854                         else
2855                                 isl_int_set(tab->dual->el[1 + i],
2856                                             tab->mat->row[var->index][pos]);
2857                 }
2858         }
2859         if (opt && res == isl_lp_ok) {
2860                 if (opt_denom) {
2861                         isl_int_set(*opt, tab->mat->row[var->index][1]);
2862                         isl_int_set(*opt_denom, tab->mat->row[var->index][0]);
2863                 } else
2864                         isl_int_cdiv_q(*opt, tab->mat->row[var->index][1],
2865                                              tab->mat->row[var->index][0]);
2866         }
2867         if (isl_tab_rollback(tab, snap) < 0)
2868                 return isl_lp_error;
2869         return res;
2870 }
2871
2872 int isl_tab_is_redundant(struct isl_tab *tab, int con)
2873 {
2874         if (!tab)
2875                 return -1;
2876         if (tab->con[con].is_zero)
2877                 return 0;
2878         if (tab->con[con].is_redundant)
2879                 return 1;
2880         return tab->con[con].is_row && tab->con[con].index < tab->n_redundant;
2881 }
2882
2883 /* Take a snapshot of the tableau that can be restored by s call to
2884  * isl_tab_rollback.
2885  */
2886 struct isl_tab_undo *isl_tab_snap(struct isl_tab *tab)
2887 {
2888         if (!tab)
2889                 return NULL;
2890         tab->need_undo = 1;
2891         return tab->top;
2892 }
2893
2894 /* Undo the operation performed by isl_tab_relax.
2895  */
2896 static int unrelax(struct isl_tab *tab, struct isl_tab_var *var) WARN_UNUSED;
2897 static int unrelax(struct isl_tab *tab, struct isl_tab_var *var)
2898 {
2899         unsigned off = 2 + tab->M;
2900
2901         if (!var->is_row && !max_is_manifestly_unbounded(tab, var))
2902                 if (to_row(tab, var, 1) < 0)
2903                         return -1;
2904
2905         if (var->is_row) {
2906                 isl_int_sub(tab->mat->row[var->index][1],
2907                     tab->mat->row[var->index][1], tab->mat->row[var->index][0]);
2908                 if (var->is_nonneg) {
2909                         int sgn = restore_row(tab, var);
2910                         isl_assert(tab->mat->ctx, sgn >= 0, return -1);
2911                 }
2912         } else {
2913                 int i;
2914
2915                 for (i = 0; i < tab->n_row; ++i) {
2916                         if (isl_int_is_zero(tab->mat->row[i][off + var->index]))
2917                                 continue;
2918                         isl_int_add(tab->mat->row[i][1], tab->mat->row[i][1],
2919                             tab->mat->row[i][off + var->index]);
2920                 }
2921
2922         }
2923
2924         return 0;
2925 }
2926
2927 static int perform_undo_var(struct isl_tab *tab, struct isl_tab_undo *undo) WARN_UNUSED;
2928 static int perform_undo_var(struct isl_tab *tab, struct isl_tab_undo *undo)
2929 {
2930         struct isl_tab_var *var = var_from_index(tab, undo->u.var_index);
2931         switch(undo->type) {
2932         case isl_tab_undo_nonneg:
2933                 var->is_nonneg = 0;
2934                 break;
2935         case isl_tab_undo_redundant:
2936                 var->is_redundant = 0;
2937                 tab->n_redundant--;
2938                 restore_row(tab, isl_tab_var_from_row(tab, tab->n_redundant));
2939                 break;
2940         case isl_tab_undo_freeze:
2941                 var->frozen = 0;
2942                 break;
2943         case isl_tab_undo_zero:
2944                 var->is_zero = 0;
2945                 if (!var->is_row)
2946                         tab->n_dead--;
2947                 break;
2948         case isl_tab_undo_allocate:
2949                 if (undo->u.var_index >= 0) {
2950                         isl_assert(tab->mat->ctx, !var->is_row, return -1);
2951                         drop_col(tab, var->index);
2952                         break;
2953                 }
2954                 if (!var->is_row) {
2955                         if (!max_is_manifestly_unbounded(tab, var)) {
2956                                 if (to_row(tab, var, 1) < 0)
2957                                         return -1;
2958                         } else if (!min_is_manifestly_unbounded(tab, var)) {
2959                                 if (to_row(tab, var, -1) < 0)
2960                                         return -1;
2961                         } else
2962                                 if (to_row(tab, var, 0) < 0)
2963                                         return -1;
2964                 }
2965                 drop_row(tab, var->index);
2966                 break;
2967         case isl_tab_undo_relax:
2968                 return unrelax(tab, var);
2969         }
2970
2971         return 0;
2972 }
2973
2974 /* Restore the tableau to the state where the basic variables
2975  * are those in "col_var".
2976  * We first construct a list of variables that are currently in
2977  * the basis, but shouldn't.  Then we iterate over all variables
2978  * that should be in the basis and for each one that is currently
2979  * not in the basis, we exchange it with one of the elements of the
2980  * list constructed before.
2981  * We can always find an appropriate variable to pivot with because
2982  * the current basis is mapped to the old basis by a non-singular
2983  * matrix and so we can never end up with a zero row.
2984  */
2985 static int restore_basis(struct isl_tab *tab, int *col_var)
2986 {
2987         int i, j;
2988         int n_extra = 0;
2989         int *extra = NULL;      /* current columns that contain bad stuff */
2990         unsigned off = 2 + tab->M;
2991
2992         extra = isl_alloc_array(tab->mat->ctx, int, tab->n_col);
2993         if (!extra)
2994                 goto error;
2995         for (i = 0; i < tab->n_col; ++i) {
2996                 for (j = 0; j < tab->n_col; ++j)
2997                         if (tab->col_var[i] == col_var[j])
2998                                 break;
2999                 if (j < tab->n_col)
3000                         continue;
3001                 extra[n_extra++] = i;
3002         }
3003         for (i = 0; i < tab->n_col && n_extra > 0; ++i) {
3004                 struct isl_tab_var *var;
3005                 int row;
3006
3007                 for (j = 0; j < tab->n_col; ++j)
3008                         if (col_var[i] == tab->col_var[j])
3009                                 break;
3010                 if (j < tab->n_col)
3011                         continue;
3012                 var = var_from_index(tab, col_var[i]);
3013                 row = var->index;
3014                 for (j = 0; j < n_extra; ++j)
3015                         if (!isl_int_is_zero(tab->mat->row[row][off+extra[j]]))
3016                                 break;
3017                 isl_assert(tab->mat->ctx, j < n_extra, goto error);
3018                 if (isl_tab_pivot(tab, row, extra[j]) < 0)
3019                         goto error;
3020                 extra[j] = extra[--n_extra];
3021         }
3022
3023         free(extra);
3024         return 0;
3025 error:
3026         free(extra);
3027         return -1;
3028 }
3029
3030 /* Remove all samples with index n or greater, i.e., those samples
3031  * that were added since we saved this number of samples in
3032  * isl_tab_save_samples.
3033  */
3034 static void drop_samples_since(struct isl_tab *tab, int n)
3035 {
3036         int i;
3037
3038         for (i = tab->n_sample - 1; i >= 0 && tab->n_sample > n; --i) {
3039                 if (tab->sample_index[i] < n)
3040                         continue;
3041
3042                 if (i != tab->n_sample - 1) {
3043                         int t = tab->sample_index[tab->n_sample-1];
3044                         tab->sample_index[tab->n_sample-1] = tab->sample_index[i];
3045                         tab->sample_index[i] = t;
3046                         isl_mat_swap_rows(tab->samples, tab->n_sample-1, i);
3047                 }
3048                 tab->n_sample--;
3049         }
3050 }
3051
3052 static int perform_undo(struct isl_tab *tab, struct isl_tab_undo *undo) WARN_UNUSED;
3053 static int perform_undo(struct isl_tab *tab, struct isl_tab_undo *undo)
3054 {
3055         switch (undo->type) {
3056         case isl_tab_undo_empty:
3057                 tab->empty = 0;
3058                 break;
3059         case isl_tab_undo_nonneg:
3060         case isl_tab_undo_redundant:
3061         case isl_tab_undo_freeze:
3062         case isl_tab_undo_zero:
3063         case isl_tab_undo_allocate:
3064         case isl_tab_undo_relax:
3065                 return perform_undo_var(tab, undo);
3066         case isl_tab_undo_bmap_eq:
3067                 return isl_basic_map_free_equality(tab->bmap, 1);
3068         case isl_tab_undo_bmap_ineq:
3069                 return isl_basic_map_free_inequality(tab->bmap, 1);
3070         case isl_tab_undo_bmap_div:
3071                 if (isl_basic_map_free_div(tab->bmap, 1) < 0)
3072                         return -1;
3073                 if (tab->samples)
3074                         tab->samples->n_col--;
3075                 break;
3076         case isl_tab_undo_saved_basis:
3077                 if (restore_basis(tab, undo->u.col_var) < 0)
3078                         return -1;
3079                 break;
3080         case isl_tab_undo_drop_sample:
3081                 tab->n_outside--;
3082                 break;
3083         case isl_tab_undo_saved_samples:
3084                 drop_samples_since(tab, undo->u.n);
3085                 break;
3086         case isl_tab_undo_callback:
3087                 return undo->u.callback->run(undo->u.callback);
3088         default:
3089                 isl_assert(tab->mat->ctx, 0, return -1);
3090         }
3091         return 0;
3092 }
3093
3094 /* Return the tableau to the state it was in when the snapshot "snap"
3095  * was taken.
3096  */
3097 int isl_tab_rollback(struct isl_tab *tab, struct isl_tab_undo *snap)
3098 {
3099         struct isl_tab_undo *undo, *next;
3100
3101         if (!tab)
3102                 return -1;
3103
3104         tab->in_undo = 1;
3105         for (undo = tab->top; undo && undo != &tab->bottom; undo = next) {
3106                 next = undo->next;
3107                 if (undo == snap)
3108                         break;
3109                 if (perform_undo(tab, undo) < 0) {
3110                         tab->top = undo;
3111                         free_undo(tab);
3112                         tab->in_undo = 0;
3113                         return -1;
3114                 }
3115                 free_undo_record(undo);
3116         }
3117         tab->in_undo = 0;
3118         tab->top = undo;
3119         if (!undo)
3120                 return -1;
3121         return 0;
3122 }
3123
3124 /* The given row "row" represents an inequality violated by all
3125  * points in the tableau.  Check for some special cases of such
3126  * separating constraints.
3127  * In particular, if the row has been reduced to the constant -1,
3128  * then we know the inequality is adjacent (but opposite) to
3129  * an equality in the tableau.
3130  * If the row has been reduced to r = c*(-1 -r'), with r' an inequality
3131  * of the tableau and c a positive constant, then the inequality
3132  * is adjacent (but opposite) to the inequality r'.
3133  */
3134 static enum isl_ineq_type separation_type(struct isl_tab *tab, unsigned row)
3135 {
3136         int pos;
3137         unsigned off = 2 + tab->M;
3138
3139         if (tab->rational)
3140                 return isl_ineq_separate;
3141
3142         if (!isl_int_is_one(tab->mat->row[row][0]))
3143                 return isl_ineq_separate;
3144
3145         pos = isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
3146                                         tab->n_col - tab->n_dead);
3147         if (pos == -1) {
3148                 if (isl_int_is_negone(tab->mat->row[row][1]))
3149                         return isl_ineq_adj_eq;
3150                 else
3151                         return isl_ineq_separate;
3152         }
3153
3154         if (!isl_int_eq(tab->mat->row[row][1],
3155                         tab->mat->row[row][off + tab->n_dead + pos]))
3156                 return isl_ineq_separate;
3157
3158         pos = isl_seq_first_non_zero(
3159                         tab->mat->row[row] + off + tab->n_dead + pos + 1,
3160                         tab->n_col - tab->n_dead - pos - 1);
3161
3162         return pos == -1 ? isl_ineq_adj_ineq : isl_ineq_separate;
3163 }
3164
3165 /* Check the effect of inequality "ineq" on the tableau "tab".
3166  * The result may be
3167  *      isl_ineq_redundant:     satisfied by all points in the tableau
3168  *      isl_ineq_separate:      satisfied by no point in the tableau
3169  *      isl_ineq_cut:           satisfied by some by not all points
3170  *      isl_ineq_adj_eq:        adjacent to an equality
3171  *      isl_ineq_adj_ineq:      adjacent to an inequality.
3172  */
3173 enum isl_ineq_type isl_tab_ineq_type(struct isl_tab *tab, isl_int *ineq)
3174 {
3175         enum isl_ineq_type type = isl_ineq_error;
3176         struct isl_tab_undo *snap = NULL;
3177         int con;
3178         int row;
3179
3180         if (!tab)
3181                 return isl_ineq_error;
3182
3183         if (isl_tab_extend_cons(tab, 1) < 0)
3184                 return isl_ineq_error;
3185
3186         snap = isl_tab_snap(tab);
3187
3188         con = isl_tab_add_row(tab, ineq);
3189         if (con < 0)
3190                 goto error;
3191
3192         row = tab->con[con].index;
3193         if (isl_tab_row_is_redundant(tab, row))
3194                 type = isl_ineq_redundant;
3195         else if (isl_int_is_neg(tab->mat->row[row][1]) &&
3196                  (tab->rational ||
3197                     isl_int_abs_ge(tab->mat->row[row][1],
3198                                    tab->mat->row[row][0]))) {
3199                 int nonneg = at_least_zero(tab, &tab->con[con]);
3200                 if (nonneg < 0)
3201                         goto error;
3202                 if (nonneg)
3203                         type = isl_ineq_cut;
3204                 else
3205                         type = separation_type(tab, row);
3206         } else {
3207                 int red = con_is_redundant(tab, &tab->con[con]);
3208                 if (red < 0)
3209                         goto error;
3210                 if (!red)
3211                         type = isl_ineq_cut;
3212                 else
3213                         type = isl_ineq_redundant;
3214         }
3215
3216         if (isl_tab_rollback(tab, snap))
3217                 return isl_ineq_error;
3218         return type;
3219 error:
3220         return isl_ineq_error;
3221 }
3222
3223 int isl_tab_track_bmap(struct isl_tab *tab, __isl_take isl_basic_map *bmap)
3224 {
3225         if (!tab || !bmap)
3226                 goto error;
3227
3228         isl_assert(tab->mat->ctx, tab->n_eq == bmap->n_eq, return -1);
3229         isl_assert(tab->mat->ctx,
3230                     tab->n_con == bmap->n_eq + bmap->n_ineq, return -1);
3231
3232         tab->bmap = bmap;
3233
3234         return 0;
3235 error:
3236         isl_basic_map_free(bmap);
3237         return -1;
3238 }
3239
3240 int isl_tab_track_bset(struct isl_tab *tab, __isl_take isl_basic_set *bset)
3241 {
3242         return isl_tab_track_bmap(tab, (isl_basic_map *)bset);
3243 }
3244
3245 __isl_keep isl_basic_set *isl_tab_peek_bset(struct isl_tab *tab)
3246 {
3247         if (!tab)
3248                 return NULL;
3249
3250         return (isl_basic_set *)tab->bmap;
3251 }
3252
3253 static void isl_tab_print_internal(__isl_keep struct isl_tab *tab,
3254         FILE *out, int indent)
3255 {
3256         unsigned r, c;
3257         int i;
3258
3259         if (!tab) {
3260                 fprintf(out, "%*snull tab\n", indent, "");
3261                 return;
3262         }
3263         fprintf(out, "%*sn_redundant: %d, n_dead: %d", indent, "",
3264                 tab->n_redundant, tab->n_dead);
3265         if (tab->rational)
3266                 fprintf(out, ", rational");
3267         if (tab->empty)
3268                 fprintf(out, ", empty");
3269         fprintf(out, "\n");
3270         fprintf(out, "%*s[", indent, "");
3271         for (i = 0; i < tab->n_var; ++i) {
3272                 if (i)
3273                         fprintf(out, (i == tab->n_param ||
3274                                       i == tab->n_var - tab->n_div) ? "; "
3275                                                                     : ", ");
3276                 fprintf(out, "%c%d%s", tab->var[i].is_row ? 'r' : 'c',
3277                                         tab->var[i].index,
3278                                         tab->var[i].is_zero ? " [=0]" :
3279                                         tab->var[i].is_redundant ? " [R]" : "");
3280         }
3281         fprintf(out, "]\n");
3282         fprintf(out, "%*s[", indent, "");
3283         for (i = 0; i < tab->n_con; ++i) {
3284                 if (i)
3285                         fprintf(out, ", ");
3286                 fprintf(out, "%c%d%s", tab->con[i].is_row ? 'r' : 'c',
3287                                         tab->con[i].index,
3288                                         tab->con[i].is_zero ? " [=0]" :
3289                                         tab->con[i].is_redundant ? " [R]" : "");
3290         }
3291         fprintf(out, "]\n");
3292         fprintf(out, "%*s[", indent, "");
3293         for (i = 0; i < tab->n_row; ++i) {
3294                 const char *sign = "";
3295                 if (i)
3296                         fprintf(out, ", ");
3297                 if (tab->row_sign) {
3298                         if (tab->row_sign[i] == isl_tab_row_unknown)
3299                                 sign = "?";
3300                         else if (tab->row_sign[i] == isl_tab_row_neg)
3301                                 sign = "-";
3302                         else if (tab->row_sign[i] == isl_tab_row_pos)
3303                                 sign = "+";
3304                         else
3305                                 sign = "+-";
3306                 }
3307                 fprintf(out, "r%d: %d%s%s", i, tab->row_var[i],
3308                     isl_tab_var_from_row(tab, i)->is_nonneg ? " [>=0]" : "", sign);
3309         }
3310         fprintf(out, "]\n");
3311         fprintf(out, "%*s[", indent, "");
3312         for (i = 0; i < tab->n_col; ++i) {
3313                 if (i)
3314                         fprintf(out, ", ");
3315                 fprintf(out, "c%d: %d%s", i, tab->col_var[i],
3316                     var_from_col(tab, i)->is_nonneg ? " [>=0]" : "");
3317         }
3318         fprintf(out, "]\n");
3319         r = tab->mat->n_row;
3320         tab->mat->n_row = tab->n_row;
3321         c = tab->mat->n_col;
3322         tab->mat->n_col = 2 + tab->M + tab->n_col;
3323         isl_mat_print_internal(tab->mat, out, indent);
3324         tab->mat->n_row = r;
3325         tab->mat->n_col = c;
3326         if (tab->bmap)
3327                 isl_basic_map_print_internal(tab->bmap, out, indent);
3328 }
3329
3330 void isl_tab_dump(__isl_keep struct isl_tab *tab)
3331 {
3332         isl_tab_print_internal(tab, stderr, 0);
3333 }