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