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.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 /* Row variable "var" is non-negative and cannot attain any values
1521  * larger than zero.  This means that the coefficients of the unrestricted
1522  * column variables are zero and that the coefficients of the non-negative
1523  * column variables are zero or negative.
1524  * Each of the non-negative variables with a negative coefficient can
1525  * then also be written as the negative sum of non-negative variables
1526  * and must therefore also be zero.
1527  */
1528 static int close_row(struct isl_tab *tab, struct isl_tab_var *var) WARN_UNUSED;
1529 static int close_row(struct isl_tab *tab, struct isl_tab_var *var)
1530 {
1531         int j;
1532         struct isl_mat *mat = tab->mat;
1533         unsigned off = 2 + tab->M;
1534
1535         isl_assert(tab->mat->ctx, var->is_nonneg, return -1);
1536         var->is_zero = 1;
1537         if (tab->need_undo)
1538                 if (isl_tab_push_var(tab, isl_tab_undo_zero, var) < 0)
1539                         return -1;
1540         for (j = tab->n_dead; j < tab->n_col; ++j) {
1541                 int recheck;
1542                 if (isl_int_is_zero(mat->row[var->index][off + j]))
1543                         continue;
1544                 isl_assert(tab->mat->ctx,
1545                     isl_int_is_neg(mat->row[var->index][off + j]), return -1);
1546                 recheck = isl_tab_kill_col(tab, j);
1547                 if (recheck < 0)
1548                         return -1;
1549                 if (recheck)
1550                         --j;
1551         }
1552         if (isl_tab_mark_redundant(tab, var->index) < 0)
1553                 return -1;
1554         return 0;
1555 }
1556
1557 /* Add a constraint to the tableau and allocate a row for it.
1558  * Return the index into the constraint array "con".
1559  */
1560 int isl_tab_allocate_con(struct isl_tab *tab)
1561 {
1562         int r;
1563
1564         isl_assert(tab->mat->ctx, tab->n_row < tab->mat->n_row, return -1);
1565         isl_assert(tab->mat->ctx, tab->n_con < tab->max_con, return -1);
1566
1567         r = tab->n_con;
1568         tab->con[r].index = tab->n_row;
1569         tab->con[r].is_row = 1;
1570         tab->con[r].is_nonneg = 0;
1571         tab->con[r].is_zero = 0;
1572         tab->con[r].is_redundant = 0;
1573         tab->con[r].frozen = 0;
1574         tab->con[r].negated = 0;
1575         tab->row_var[tab->n_row] = ~r;
1576
1577         tab->n_row++;
1578         tab->n_con++;
1579         if (isl_tab_push_var(tab, isl_tab_undo_allocate, &tab->con[r]) < 0)
1580                 return -1;
1581
1582         return r;
1583 }
1584
1585 /* Add a variable to the tableau and allocate a column for it.
1586  * Return the index into the variable array "var".
1587  */
1588 int isl_tab_allocate_var(struct isl_tab *tab)
1589 {
1590         int r;
1591         int i;
1592         unsigned off = 2 + tab->M;
1593
1594         isl_assert(tab->mat->ctx, tab->n_col < tab->mat->n_col, return -1);
1595         isl_assert(tab->mat->ctx, tab->n_var < tab->max_var, return -1);
1596
1597         r = tab->n_var;
1598         tab->var[r].index = tab->n_col;
1599         tab->var[r].is_row = 0;
1600         tab->var[r].is_nonneg = 0;
1601         tab->var[r].is_zero = 0;
1602         tab->var[r].is_redundant = 0;
1603         tab->var[r].frozen = 0;
1604         tab->var[r].negated = 0;
1605         tab->col_var[tab->n_col] = r;
1606
1607         for (i = 0; i < tab->n_row; ++i)
1608                 isl_int_set_si(tab->mat->row[i][off + tab->n_col], 0);
1609
1610         tab->n_var++;
1611         tab->n_col++;
1612         if (isl_tab_push_var(tab, isl_tab_undo_allocate, &tab->var[r]) < 0)
1613                 return -1;
1614
1615         return r;
1616 }
1617
1618 /* Add a row to the tableau.  The row is given as an affine combination
1619  * of the original variables and needs to be expressed in terms of the
1620  * column variables.
1621  *
1622  * We add each term in turn.
1623  * If r = n/d_r is the current sum and we need to add k x, then
1624  *      if x is a column variable, we increase the numerator of
1625  *              this column by k d_r
1626  *      if x = f/d_x is a row variable, then the new representation of r is
1627  *
1628  *               n    k f   d_x/g n + d_r/g k f   m/d_r n + m/d_g k f
1629  *              --- + --- = ------------------- = -------------------
1630  *              d_r   d_r        d_r d_x/g                m
1631  *
1632  *      with g the gcd of d_r and d_x and m the lcm of d_r and d_x.
1633  */
1634 int isl_tab_add_row(struct isl_tab *tab, isl_int *line)
1635 {
1636         int i;
1637         int r;
1638         isl_int *row;
1639         isl_int a, b;
1640         unsigned off = 2 + tab->M;
1641
1642         r = isl_tab_allocate_con(tab);
1643         if (r < 0)
1644                 return -1;
1645
1646         isl_int_init(a);
1647         isl_int_init(b);
1648         row = tab->mat->row[tab->con[r].index];
1649         isl_int_set_si(row[0], 1);
1650         isl_int_set(row[1], line[0]);
1651         isl_seq_clr(row + 2, tab->M + tab->n_col);
1652         for (i = 0; i < tab->n_var; ++i) {
1653                 if (tab->var[i].is_zero)
1654                         continue;
1655                 if (tab->var[i].is_row) {
1656                         isl_int_lcm(a,
1657                                 row[0], tab->mat->row[tab->var[i].index][0]);
1658                         isl_int_swap(a, row[0]);
1659                         isl_int_divexact(a, row[0], a);
1660                         isl_int_divexact(b,
1661                                 row[0], tab->mat->row[tab->var[i].index][0]);
1662                         isl_int_mul(b, b, line[1 + i]);
1663                         isl_seq_combine(row + 1, a, row + 1,
1664                             b, tab->mat->row[tab->var[i].index] + 1,
1665                             1 + tab->M + tab->n_col);
1666                 } else
1667                         isl_int_addmul(row[off + tab->var[i].index],
1668                                                         line[1 + i], row[0]);
1669                 if (tab->M && i >= tab->n_param && i < tab->n_var - tab->n_div)
1670                         isl_int_submul(row[2], line[1 + i], row[0]);
1671         }
1672         isl_seq_normalize(tab->mat->ctx, row, off + tab->n_col);
1673         isl_int_clear(a);
1674         isl_int_clear(b);
1675
1676         if (tab->row_sign)
1677                 tab->row_sign[tab->con[r].index] = isl_tab_row_unknown;
1678
1679         return r;
1680 }
1681
1682 static int drop_row(struct isl_tab *tab, int row)
1683 {
1684         isl_assert(tab->mat->ctx, ~tab->row_var[row] == tab->n_con - 1, return -1);
1685         if (row != tab->n_row - 1)
1686                 swap_rows(tab, row, tab->n_row - 1);
1687         tab->n_row--;
1688         tab->n_con--;
1689         return 0;
1690 }
1691
1692 static int drop_col(struct isl_tab *tab, int col)
1693 {
1694         isl_assert(tab->mat->ctx, tab->col_var[col] == tab->n_var - 1, return -1);
1695         if (col != tab->n_col - 1)
1696                 swap_cols(tab, col, tab->n_col - 1);
1697         tab->n_col--;
1698         tab->n_var--;
1699         return 0;
1700 }
1701
1702 /* Add inequality "ineq" and check if it conflicts with the
1703  * previously added constraints or if it is obviously redundant.
1704  */
1705 int isl_tab_add_ineq(struct isl_tab *tab, isl_int *ineq)
1706 {
1707         int r;
1708         int sgn;
1709         isl_int cst;
1710
1711         if (!tab)
1712                 return -1;
1713         if (tab->bmap) {
1714                 struct isl_basic_map *bmap = tab->bmap;
1715
1716                 isl_assert(tab->mat->ctx, tab->n_eq == bmap->n_eq, return -1);
1717                 isl_assert(tab->mat->ctx,
1718                             tab->n_con == bmap->n_eq + bmap->n_ineq, return -1);
1719                 tab->bmap = isl_basic_map_add_ineq(tab->bmap, ineq);
1720                 if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1721                         return -1;
1722                 if (!tab->bmap)
1723                         return -1;
1724         }
1725         if (tab->cone) {
1726                 isl_int_init(cst);
1727                 isl_int_swap(ineq[0], cst);
1728         }
1729         r = isl_tab_add_row(tab, ineq);
1730         if (tab->cone) {
1731                 isl_int_swap(ineq[0], cst);
1732                 isl_int_clear(cst);
1733         }
1734         if (r < 0)
1735                 return -1;
1736         tab->con[r].is_nonneg = 1;
1737         if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1738                 return -1;
1739         if (isl_tab_row_is_redundant(tab, tab->con[r].index)) {
1740                 if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0)
1741                         return -1;
1742                 return 0;
1743         }
1744
1745         sgn = restore_row(tab, &tab->con[r]);
1746         if (sgn < -1)
1747                 return -1;
1748         if (sgn < 0)
1749                 return isl_tab_mark_empty(tab);
1750         if (tab->con[r].is_row && isl_tab_row_is_redundant(tab, tab->con[r].index))
1751                 if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0)
1752                         return -1;
1753         return 0;
1754 }
1755
1756 /* Pivot a non-negative variable down until it reaches the value zero
1757  * and then pivot the variable into a column position.
1758  */
1759 static int to_col(struct isl_tab *tab, struct isl_tab_var *var) WARN_UNUSED;
1760 static int to_col(struct isl_tab *tab, struct isl_tab_var *var)
1761 {
1762         int i;
1763         int row, col;
1764         unsigned off = 2 + tab->M;
1765
1766         if (!var->is_row)
1767                 return 0;
1768
1769         while (isl_int_is_pos(tab->mat->row[var->index][1])) {
1770                 find_pivot(tab, var, NULL, -1, &row, &col);
1771                 isl_assert(tab->mat->ctx, row != -1, return -1);
1772                 if (isl_tab_pivot(tab, row, col) < 0)
1773                         return -1;
1774                 if (!var->is_row)
1775                         return 0;
1776         }
1777
1778         for (i = tab->n_dead; i < tab->n_col; ++i)
1779                 if (!isl_int_is_zero(tab->mat->row[var->index][off + i]))
1780                         break;
1781
1782         isl_assert(tab->mat->ctx, i < tab->n_col, return -1);
1783         if (isl_tab_pivot(tab, var->index, i) < 0)
1784                 return -1;
1785
1786         return 0;
1787 }
1788
1789 /* We assume Gaussian elimination has been performed on the equalities.
1790  * The equalities can therefore never conflict.
1791  * Adding the equalities is currently only really useful for a later call
1792  * to isl_tab_ineq_type.
1793  */
1794 static struct isl_tab *add_eq(struct isl_tab *tab, isl_int *eq)
1795 {
1796         int i;
1797         int r;
1798
1799         if (!tab)
1800                 return NULL;
1801         r = isl_tab_add_row(tab, eq);
1802         if (r < 0)
1803                 goto error;
1804
1805         r = tab->con[r].index;
1806         i = isl_seq_first_non_zero(tab->mat->row[r] + 2 + tab->M + tab->n_dead,
1807                                         tab->n_col - tab->n_dead);
1808         isl_assert(tab->mat->ctx, i >= 0, goto error);
1809         i += tab->n_dead;
1810         if (isl_tab_pivot(tab, r, i) < 0)
1811                 goto error;
1812         if (isl_tab_kill_col(tab, i) < 0)
1813                 goto error;
1814         tab->n_eq++;
1815
1816         return tab;
1817 error:
1818         isl_tab_free(tab);
1819         return NULL;
1820 }
1821
1822 static int row_is_manifestly_zero(struct isl_tab *tab, int row)
1823 {
1824         unsigned off = 2 + tab->M;
1825
1826         if (!isl_int_is_zero(tab->mat->row[row][1]))
1827                 return 0;
1828         if (tab->M && !isl_int_is_zero(tab->mat->row[row][2]))
1829                 return 0;
1830         return isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
1831                                         tab->n_col - tab->n_dead) == -1;
1832 }
1833
1834 /* Add an equality that is known to be valid for the given tableau.
1835  */
1836 int isl_tab_add_valid_eq(struct isl_tab *tab, isl_int *eq)
1837 {
1838         struct isl_tab_var *var;
1839         int r;
1840
1841         if (!tab)
1842                 return -1;
1843         r = isl_tab_add_row(tab, eq);
1844         if (r < 0)
1845                 return -1;
1846
1847         var = &tab->con[r];
1848         r = var->index;
1849         if (row_is_manifestly_zero(tab, r)) {
1850                 var->is_zero = 1;
1851                 if (isl_tab_mark_redundant(tab, r) < 0)
1852                         return -1;
1853                 return 0;
1854         }
1855
1856         if (isl_int_is_neg(tab->mat->row[r][1])) {
1857                 isl_seq_neg(tab->mat->row[r] + 1, tab->mat->row[r] + 1,
1858                             1 + tab->n_col);
1859                 var->negated = 1;
1860         }
1861         var->is_nonneg = 1;
1862         if (to_col(tab, var) < 0)
1863                 return -1;
1864         var->is_nonneg = 0;
1865         if (isl_tab_kill_col(tab, var->index) < 0)
1866                 return -1;
1867
1868         return 0;
1869 }
1870
1871 static int add_zero_row(struct isl_tab *tab)
1872 {
1873         int r;
1874         isl_int *row;
1875
1876         r = isl_tab_allocate_con(tab);
1877         if (r < 0)
1878                 return -1;
1879
1880         row = tab->mat->row[tab->con[r].index];
1881         isl_seq_clr(row + 1, 1 + tab->M + tab->n_col);
1882         isl_int_set_si(row[0], 1);
1883
1884         return r;
1885 }
1886
1887 /* Add equality "eq" and check if it conflicts with the
1888  * previously added constraints or if it is obviously redundant.
1889  */
1890 int isl_tab_add_eq(struct isl_tab *tab, isl_int *eq)
1891 {
1892         struct isl_tab_undo *snap = NULL;
1893         struct isl_tab_var *var;
1894         int r;
1895         int row;
1896         int sgn;
1897         isl_int cst;
1898
1899         if (!tab)
1900                 return -1;
1901         isl_assert(tab->mat->ctx, !tab->M, return -1);
1902
1903         if (tab->need_undo)
1904                 snap = isl_tab_snap(tab);
1905
1906         if (tab->cone) {
1907                 isl_int_init(cst);
1908                 isl_int_swap(eq[0], cst);
1909         }
1910         r = isl_tab_add_row(tab, eq);
1911         if (tab->cone) {
1912                 isl_int_swap(eq[0], cst);
1913                 isl_int_clear(cst);
1914         }
1915         if (r < 0)
1916                 return -1;
1917
1918         var = &tab->con[r];
1919         row = var->index;
1920         if (row_is_manifestly_zero(tab, row)) {
1921                 if (snap) {
1922                         if (isl_tab_rollback(tab, snap) < 0)
1923                                 return -1;
1924                 } else
1925                         drop_row(tab, row);
1926                 return 0;
1927         }
1928
1929         if (tab->bmap) {
1930                 tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);
1931                 if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1932                         return -1;
1933                 isl_seq_neg(eq, eq, 1 + tab->n_var);
1934                 tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);
1935                 isl_seq_neg(eq, eq, 1 + tab->n_var);
1936                 if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1937                         return -1;
1938                 if (!tab->bmap)
1939                         return -1;
1940                 if (add_zero_row(tab) < 0)
1941                         return -1;
1942         }
1943
1944         sgn = isl_int_sgn(tab->mat->row[row][1]);
1945
1946         if (sgn > 0) {
1947                 isl_seq_neg(tab->mat->row[row] + 1, tab->mat->row[row] + 1,
1948                             1 + tab->n_col);
1949                 var->negated = 1;
1950                 sgn = -1;
1951         }
1952
1953         if (sgn < 0) {
1954                 sgn = sign_of_max(tab, var);
1955                 if (sgn < -1)
1956                         return -1;
1957                 if (sgn < 0) {
1958                         if (isl_tab_mark_empty(tab) < 0)
1959                                 return -1;
1960                         return 0;
1961                 }
1962         }
1963
1964         var->is_nonneg = 1;
1965         if (to_col(tab, var) < 0)
1966                 return -1;
1967         var->is_nonneg = 0;
1968         if (isl_tab_kill_col(tab, var->index) < 0)
1969                 return -1;
1970
1971         return 0;
1972 }
1973
1974 /* Construct and return an inequality that expresses an upper bound
1975  * on the given div.
1976  * In particular, if the div is given by
1977  *
1978  *      d = floor(e/m)
1979  *
1980  * then the inequality expresses
1981  *
1982  *      m d <= e
1983  */
1984 static struct isl_vec *ineq_for_div(struct isl_basic_map *bmap, unsigned div)
1985 {
1986         unsigned total;
1987         unsigned div_pos;
1988         struct isl_vec *ineq;
1989
1990         if (!bmap)
1991                 return NULL;
1992
1993         total = isl_basic_map_total_dim(bmap);
1994         div_pos = 1 + total - bmap->n_div + div;
1995
1996         ineq = isl_vec_alloc(bmap->ctx, 1 + total);
1997         if (!ineq)
1998                 return NULL;
1999
2000         isl_seq_cpy(ineq->el, bmap->div[div] + 1, 1 + total);
2001         isl_int_neg(ineq->el[div_pos], bmap->div[div][0]);
2002         return ineq;
2003 }
2004
2005 /* For a div d = floor(f/m), add the constraints
2006  *
2007  *              f - m d >= 0
2008  *              -(f-(m-1)) + m d >= 0
2009  *
2010  * Note that the second constraint is the negation of
2011  *
2012  *              f - m d >= m
2013  *
2014  * If add_ineq is not NULL, then this function is used
2015  * instead of isl_tab_add_ineq to effectively add the inequalities.
2016  */
2017 static int add_div_constraints(struct isl_tab *tab, unsigned div,
2018         int (*add_ineq)(void *user, isl_int *), void *user)
2019 {
2020         unsigned total;
2021         unsigned div_pos;
2022         struct isl_vec *ineq;
2023
2024         total = isl_basic_map_total_dim(tab->bmap);
2025         div_pos = 1 + total - tab->bmap->n_div + div;
2026
2027         ineq = ineq_for_div(tab->bmap, div);
2028         if (!ineq)
2029                 goto error;
2030
2031         if (add_ineq) {
2032                 if (add_ineq(user, ineq->el) < 0)
2033                         goto error;
2034         } else {
2035                 if (isl_tab_add_ineq(tab, ineq->el) < 0)
2036                         goto error;
2037         }
2038
2039         isl_seq_neg(ineq->el, tab->bmap->div[div] + 1, 1 + total);
2040         isl_int_set(ineq->el[div_pos], tab->bmap->div[div][0]);
2041         isl_int_add(ineq->el[0], ineq->el[0], ineq->el[div_pos]);
2042         isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
2043
2044         if (add_ineq) {
2045                 if (add_ineq(user, ineq->el) < 0)
2046                         goto error;
2047         } else {
2048                 if (isl_tab_add_ineq(tab, ineq->el) < 0)
2049                         goto error;
2050         }
2051
2052         isl_vec_free(ineq);
2053
2054         return 0;
2055 error:
2056         isl_vec_free(ineq);
2057         return -1;
2058 }
2059
2060 /* Add an extra div, prescrived by "div" to the tableau and
2061  * the associated bmap (which is assumed to be non-NULL).
2062  *
2063  * If add_ineq is not NULL, then this function is used instead
2064  * of isl_tab_add_ineq to add the div constraints.
2065  * This complication is needed because the code in isl_tab_pip
2066  * wants to perform some extra processing when an inequality
2067  * is added to the tableau.
2068  */
2069 int isl_tab_add_div(struct isl_tab *tab, __isl_keep isl_vec *div,
2070         int (*add_ineq)(void *user, isl_int *), void *user)
2071 {
2072         int i;
2073         int r;
2074         int k;
2075         int nonneg;
2076
2077         if (!tab || !div)
2078                 return -1;
2079
2080         isl_assert(tab->mat->ctx, tab->bmap, return -1);
2081
2082         for (i = 0; i < tab->n_var; ++i) {
2083                 if (isl_int_is_neg(div->el[2 + i]))
2084                         break;
2085                 if (isl_int_is_zero(div->el[2 + i]))
2086                         continue;
2087                 if (!tab->var[i].is_nonneg)
2088                         break;
2089         }
2090         nonneg = i == tab->n_var && !isl_int_is_neg(div->el[1]);
2091
2092         if (isl_tab_extend_cons(tab, 3) < 0)
2093                 return -1;
2094         if (isl_tab_extend_vars(tab, 1) < 0)
2095                 return -1;
2096         r = isl_tab_allocate_var(tab);
2097         if (r < 0)
2098                 return -1;
2099
2100         if (nonneg)
2101                 tab->var[r].is_nonneg = 1;
2102
2103         tab->bmap = isl_basic_map_extend_dim(tab->bmap,
2104                 isl_basic_map_get_dim(tab->bmap), 1, 0, 2);
2105         k = isl_basic_map_alloc_div(tab->bmap);
2106         if (k < 0)
2107                 return -1;
2108         isl_seq_cpy(tab->bmap->div[k], div->el, div->size);
2109         if (isl_tab_push(tab, isl_tab_undo_bmap_div) < 0)
2110                 return -1;
2111
2112         if (add_div_constraints(tab, k, add_ineq, user) < 0)
2113                 return -1;
2114
2115         return r;
2116 }
2117
2118 struct isl_tab *isl_tab_from_basic_map(struct isl_basic_map *bmap)
2119 {
2120         int i;
2121         struct isl_tab *tab;
2122
2123         if (!bmap)
2124                 return NULL;
2125         tab = isl_tab_alloc(bmap->ctx,
2126                             isl_basic_map_total_dim(bmap) + bmap->n_ineq + 1,
2127                             isl_basic_map_total_dim(bmap), 0);
2128         if (!tab)
2129                 return NULL;
2130         tab->rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
2131         if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY)) {
2132                 if (isl_tab_mark_empty(tab) < 0)
2133                         goto error;
2134                 return tab;
2135         }
2136         for (i = 0; i < bmap->n_eq; ++i) {
2137                 tab = add_eq(tab, bmap->eq[i]);
2138                 if (!tab)
2139                         return tab;
2140         }
2141         for (i = 0; i < bmap->n_ineq; ++i) {
2142                 if (isl_tab_add_ineq(tab, bmap->ineq[i]) < 0)
2143                         goto error;
2144                 if (tab->empty)
2145                         return tab;
2146         }
2147         return tab;
2148 error:
2149         isl_tab_free(tab);
2150         return NULL;
2151 }
2152
2153 struct isl_tab *isl_tab_from_basic_set(struct isl_basic_set *bset)
2154 {
2155         return isl_tab_from_basic_map((struct isl_basic_map *)bset);
2156 }
2157
2158 /* Construct a tableau corresponding to the recession cone of "bset".
2159  */
2160 struct isl_tab *isl_tab_from_recession_cone(__isl_keep isl_basic_set *bset,
2161         int parametric)
2162 {
2163         isl_int cst;
2164         int i;
2165         struct isl_tab *tab;
2166         unsigned offset = 0;
2167
2168         if (!bset)
2169                 return NULL;
2170         if (parametric)
2171                 offset = isl_basic_set_dim(bset, isl_dim_param);
2172         tab = isl_tab_alloc(bset->ctx, bset->n_eq + bset->n_ineq,
2173                                 isl_basic_set_total_dim(bset) - offset, 0);
2174         if (!tab)
2175                 return NULL;
2176         tab->rational = ISL_F_ISSET(bset, ISL_BASIC_SET_RATIONAL);
2177         tab->cone = 1;
2178
2179         isl_int_init(cst);
2180         for (i = 0; i < bset->n_eq; ++i) {
2181                 isl_int_swap(bset->eq[i][offset], cst);
2182                 if (offset > 0) {
2183                         if (isl_tab_add_eq(tab, bset->eq[i] + offset) < 0)
2184                                 goto error;
2185                 } else
2186                         tab = add_eq(tab, bset->eq[i]);
2187                 isl_int_swap(bset->eq[i][offset], cst);
2188                 if (!tab)
2189                         goto done;
2190         }
2191         for (i = 0; i < bset->n_ineq; ++i) {
2192                 int r;
2193                 isl_int_swap(bset->ineq[i][offset], cst);
2194                 r = isl_tab_add_row(tab, bset->ineq[i] + offset);
2195                 isl_int_swap(bset->ineq[i][offset], cst);
2196                 if (r < 0)
2197                         goto error;
2198                 tab->con[r].is_nonneg = 1;
2199                 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
2200                         goto error;
2201         }
2202 done:
2203         isl_int_clear(cst);
2204         return tab;
2205 error:
2206         isl_int_clear(cst);
2207         isl_tab_free(tab);
2208         return NULL;
2209 }
2210
2211 /* Assuming "tab" is the tableau of a cone, check if the cone is
2212  * bounded, i.e., if it is empty or only contains the origin.
2213  */
2214 int isl_tab_cone_is_bounded(struct isl_tab *tab)
2215 {
2216         int i;
2217
2218         if (!tab)
2219                 return -1;
2220         if (tab->empty)
2221                 return 1;
2222         if (tab->n_dead == tab->n_col)
2223                 return 1;
2224
2225         for (;;) {
2226                 for (i = tab->n_redundant; i < tab->n_row; ++i) {
2227                         struct isl_tab_var *var;
2228                         int sgn;
2229                         var = isl_tab_var_from_row(tab, i);
2230                         if (!var->is_nonneg)
2231                                 continue;
2232                         sgn = sign_of_max(tab, var);
2233                         if (sgn < -1)
2234                                 return -1;
2235                         if (sgn != 0)
2236                                 return 0;
2237                         if (close_row(tab, var) < 0)
2238                                 return -1;
2239                         break;
2240                 }
2241                 if (tab->n_dead == tab->n_col)
2242                         return 1;
2243                 if (i == tab->n_row)
2244                         return 0;
2245         }
2246 }
2247
2248 int isl_tab_sample_is_integer(struct isl_tab *tab)
2249 {
2250         int i;
2251
2252         if (!tab)
2253                 return -1;
2254
2255         for (i = 0; i < tab->n_var; ++i) {
2256                 int row;
2257                 if (!tab->var[i].is_row)
2258                         continue;
2259                 row = tab->var[i].index;
2260                 if (!isl_int_is_divisible_by(tab->mat->row[row][1],
2261                                                 tab->mat->row[row][0]))
2262                         return 0;
2263         }
2264         return 1;
2265 }
2266
2267 static struct isl_vec *extract_integer_sample(struct isl_tab *tab)
2268 {
2269         int i;
2270         struct isl_vec *vec;
2271
2272         vec = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
2273         if (!vec)
2274                 return NULL;
2275
2276         isl_int_set_si(vec->block.data[0], 1);
2277         for (i = 0; i < tab->n_var; ++i) {
2278                 if (!tab->var[i].is_row)
2279                         isl_int_set_si(vec->block.data[1 + i], 0);
2280                 else {
2281                         int row = tab->var[i].index;
2282                         isl_int_divexact(vec->block.data[1 + i],
2283                                 tab->mat->row[row][1], tab->mat->row[row][0]);
2284                 }
2285         }
2286
2287         return vec;
2288 }
2289
2290 struct isl_vec *isl_tab_get_sample_value(struct isl_tab *tab)
2291 {
2292         int i;
2293         struct isl_vec *vec;
2294         isl_int m;
2295
2296         if (!tab)
2297                 return NULL;
2298
2299         vec = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
2300         if (!vec)
2301                 return NULL;
2302
2303         isl_int_init(m);
2304
2305         isl_int_set_si(vec->block.data[0], 1);
2306         for (i = 0; i < tab->n_var; ++i) {
2307                 int row;
2308                 if (!tab->var[i].is_row) {
2309                         isl_int_set_si(vec->block.data[1 + i], 0);
2310                         continue;
2311                 }
2312                 row = tab->var[i].index;
2313                 isl_int_gcd(m, vec->block.data[0], tab->mat->row[row][0]);
2314                 isl_int_divexact(m, tab->mat->row[row][0], m);
2315                 isl_seq_scale(vec->block.data, vec->block.data, m, 1 + i);
2316                 isl_int_divexact(m, vec->block.data[0], tab->mat->row[row][0]);
2317                 isl_int_mul(vec->block.data[1 + i], m, tab->mat->row[row][1]);
2318         }
2319         vec = isl_vec_normalize(vec);
2320
2321         isl_int_clear(m);
2322         return vec;
2323 }
2324
2325 /* Update "bmap" based on the results of the tableau "tab".
2326  * In particular, implicit equalities are made explicit, redundant constraints
2327  * are removed and if the sample value happens to be integer, it is stored
2328  * in "bmap" (unless "bmap" already had an integer sample).
2329  *
2330  * The tableau is assumed to have been created from "bmap" using
2331  * isl_tab_from_basic_map.
2332  */
2333 struct isl_basic_map *isl_basic_map_update_from_tab(struct isl_basic_map *bmap,
2334         struct isl_tab *tab)
2335 {
2336         int i;
2337         unsigned n_eq;
2338
2339         if (!bmap)
2340                 return NULL;
2341         if (!tab)
2342                 return bmap;
2343
2344         n_eq = tab->n_eq;
2345         if (tab->empty)
2346                 bmap = isl_basic_map_set_to_empty(bmap);
2347         else
2348                 for (i = bmap->n_ineq - 1; i >= 0; --i) {
2349                         if (isl_tab_is_equality(tab, n_eq + i))
2350                                 isl_basic_map_inequality_to_equality(bmap, i);
2351                         else if (isl_tab_is_redundant(tab, n_eq + i))
2352                                 isl_basic_map_drop_inequality(bmap, i);
2353                 }
2354         if (bmap->n_eq != n_eq)
2355                 isl_basic_map_gauss(bmap, NULL);
2356         if (!tab->rational &&
2357             !bmap->sample && isl_tab_sample_is_integer(tab))
2358                 bmap->sample = extract_integer_sample(tab);
2359         return bmap;
2360 }
2361
2362 struct isl_basic_set *isl_basic_set_update_from_tab(struct isl_basic_set *bset,
2363         struct isl_tab *tab)
2364 {
2365         return (struct isl_basic_set *)isl_basic_map_update_from_tab(
2366                 (struct isl_basic_map *)bset, tab);
2367 }
2368
2369 /* Given a non-negative variable "var", add a new non-negative variable
2370  * that is the opposite of "var", ensuring that var can only attain the
2371  * value zero.
2372  * If var = n/d is a row variable, then the new variable = -n/d.
2373  * If var is a column variables, then the new variable = -var.
2374  * If the new variable cannot attain non-negative values, then
2375  * the resulting tableau is empty.
2376  * Otherwise, we know the value will be zero and we close the row.
2377  */
2378 static int cut_to_hyperplane(struct isl_tab *tab, struct isl_tab_var *var)
2379 {
2380         unsigned r;
2381         isl_int *row;
2382         int sgn;
2383         unsigned off = 2 + tab->M;
2384
2385         if (var->is_zero)
2386                 return 0;
2387         isl_assert(tab->mat->ctx, !var->is_redundant, return -1);
2388         isl_assert(tab->mat->ctx, var->is_nonneg, return -1);
2389
2390         if (isl_tab_extend_cons(tab, 1) < 0)
2391                 return -1;
2392
2393         r = tab->n_con;
2394         tab->con[r].index = tab->n_row;
2395         tab->con[r].is_row = 1;
2396         tab->con[r].is_nonneg = 0;
2397         tab->con[r].is_zero = 0;
2398         tab->con[r].is_redundant = 0;
2399         tab->con[r].frozen = 0;
2400         tab->con[r].negated = 0;
2401         tab->row_var[tab->n_row] = ~r;
2402         row = tab->mat->row[tab->n_row];
2403
2404         if (var->is_row) {
2405                 isl_int_set(row[0], tab->mat->row[var->index][0]);
2406                 isl_seq_neg(row + 1,
2407                             tab->mat->row[var->index] + 1, 1 + tab->n_col);
2408         } else {
2409                 isl_int_set_si(row[0], 1);
2410                 isl_seq_clr(row + 1, 1 + tab->n_col);
2411                 isl_int_set_si(row[off + var->index], -1);
2412         }
2413
2414         tab->n_row++;
2415         tab->n_con++;
2416         if (isl_tab_push_var(tab, isl_tab_undo_allocate, &tab->con[r]) < 0)
2417                 return -1;
2418
2419         sgn = sign_of_max(tab, &tab->con[r]);
2420         if (sgn < -1)
2421                 return -1;
2422         if (sgn < 0) {
2423                 if (isl_tab_mark_empty(tab) < 0)
2424                         return -1;
2425                 return 0;
2426         }
2427         tab->con[r].is_nonneg = 1;
2428         if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
2429                 return -1;
2430         /* sgn == 0 */
2431         if (close_row(tab, &tab->con[r]) < 0)
2432                 return -1;
2433
2434         return 0;
2435 }
2436
2437 /* Given a tableau "tab" and an inequality constraint "con" of the tableau,
2438  * relax the inequality by one.  That is, the inequality r >= 0 is replaced
2439  * by r' = r + 1 >= 0.
2440  * If r is a row variable, we simply increase the constant term by one
2441  * (taking into account the denominator).
2442  * If r is a column variable, then we need to modify each row that
2443  * refers to r = r' - 1 by substituting this equality, effectively
2444  * subtracting the coefficient of the column from the constant.
2445  * We should only do this if the minimum is manifestly unbounded,
2446  * however.  Otherwise, we may end up with negative sample values
2447  * for non-negative variables.
2448  * So, if r is a column variable with a minimum that is not
2449  * manifestly unbounded, then we need to move it to a row.
2450  * However, the sample value of this row may be negative,
2451  * even after the relaxation, so we need to restore it.
2452  * We therefore prefer to pivot a column up to a row, if possible.
2453  */
2454 struct isl_tab *isl_tab_relax(struct isl_tab *tab, int con)
2455 {
2456         struct isl_tab_var *var;
2457         unsigned off = 2 + tab->M;
2458
2459         if (!tab)
2460                 return NULL;
2461
2462         var = &tab->con[con];
2463
2464         if (!var->is_row && !max_is_manifestly_unbounded(tab, var))
2465                 if (to_row(tab, var, 1) < 0)
2466                         goto error;
2467         if (!var->is_row && !min_is_manifestly_unbounded(tab, var))
2468                 if (to_row(tab, var, -1) < 0)
2469                         goto error;
2470
2471         if (var->is_row) {
2472                 isl_int_add(tab->mat->row[var->index][1],
2473                     tab->mat->row[var->index][1], tab->mat->row[var->index][0]);
2474                 if (restore_row(tab, var) < 0)
2475                         goto error;
2476         } else {
2477                 int i;
2478
2479                 for (i = 0; i < tab->n_row; ++i) {
2480                         if (isl_int_is_zero(tab->mat->row[i][off + var->index]))
2481                                 continue;
2482                         isl_int_sub(tab->mat->row[i][1], tab->mat->row[i][1],
2483                             tab->mat->row[i][off + var->index]);
2484                 }
2485
2486         }
2487
2488         if (isl_tab_push_var(tab, isl_tab_undo_relax, var) < 0)
2489                 goto error;
2490
2491         return tab;
2492 error:
2493         isl_tab_free(tab);
2494         return NULL;
2495 }
2496
2497 int isl_tab_select_facet(struct isl_tab *tab, int con)
2498 {
2499         if (!tab)
2500                 return -1;
2501
2502         return cut_to_hyperplane(tab, &tab->con[con]);
2503 }
2504
2505 static int may_be_equality(struct isl_tab *tab, int row)
2506 {
2507         unsigned off = 2 + tab->M;
2508         return tab->rational ? isl_int_is_zero(tab->mat->row[row][1])
2509                              : isl_int_lt(tab->mat->row[row][1],
2510                                             tab->mat->row[row][0]);
2511 }
2512
2513 /* Check for (near) equalities among the constraints.
2514  * A constraint is an equality if it is non-negative and if
2515  * its maximal value is either
2516  *      - zero (in case of rational tableaus), or
2517  *      - strictly less than 1 (in case of integer tableaus)
2518  *
2519  * We first mark all non-redundant and non-dead variables that
2520  * are not frozen and not obviously not an equality.
2521  * Then we iterate over all marked variables if they can attain
2522  * any values larger than zero or at least one.
2523  * If the maximal value is zero, we mark any column variables
2524  * that appear in the row as being zero and mark the row as being redundant.
2525  * Otherwise, if the maximal value is strictly less than one (and the
2526  * tableau is integer), then we restrict the value to being zero
2527  * by adding an opposite non-negative variable.
2528  */
2529 int isl_tab_detect_implicit_equalities(struct isl_tab *tab)
2530 {
2531         int i;
2532         unsigned n_marked;
2533
2534         if (!tab)
2535                 return -1;
2536         if (tab->empty)
2537                 return 0;
2538         if (tab->n_dead == tab->n_col)
2539                 return 0;
2540
2541         n_marked = 0;
2542         for (i = tab->n_redundant; i < tab->n_row; ++i) {
2543                 struct isl_tab_var *var = isl_tab_var_from_row(tab, i);
2544                 var->marked = !var->frozen && var->is_nonneg &&
2545                         may_be_equality(tab, i);
2546                 if (var->marked)
2547                         n_marked++;
2548         }
2549         for (i = tab->n_dead; i < tab->n_col; ++i) {
2550                 struct isl_tab_var *var = var_from_col(tab, i);
2551                 var->marked = !var->frozen && var->is_nonneg;
2552                 if (var->marked)
2553                         n_marked++;
2554         }
2555         while (n_marked) {
2556                 struct isl_tab_var *var;
2557                 int sgn;
2558                 for (i = tab->n_redundant; i < tab->n_row; ++i) {
2559                         var = isl_tab_var_from_row(tab, i);
2560                         if (var->marked)
2561                                 break;
2562                 }
2563                 if (i == tab->n_row) {
2564                         for (i = tab->n_dead; i < tab->n_col; ++i) {
2565                                 var = var_from_col(tab, i);
2566                                 if (var->marked)
2567                                         break;
2568                         }
2569                         if (i == tab->n_col)
2570                                 break;
2571                 }
2572                 var->marked = 0;
2573                 n_marked--;
2574                 sgn = sign_of_max(tab, var);
2575                 if (sgn < 0)
2576                         return -1;
2577                 if (sgn == 0) {
2578                         if (close_row(tab, var) < 0)
2579                                 return -1;
2580                 } else if (!tab->rational && !at_least_one(tab, var)) {
2581                         if (cut_to_hyperplane(tab, var) < 0)
2582                                 return -1;
2583                         return isl_tab_detect_implicit_equalities(tab);
2584                 }
2585                 for (i = tab->n_redundant; i < tab->n_row; ++i) {
2586                         var = isl_tab_var_from_row(tab, i);
2587                         if (!var->marked)
2588                                 continue;
2589                         if (may_be_equality(tab, i))
2590                                 continue;
2591                         var->marked = 0;
2592                         n_marked--;
2593                 }
2594         }
2595
2596         return 0;
2597 }
2598
2599 static int con_is_redundant(struct isl_tab *tab, struct isl_tab_var *var)
2600 {
2601         if (!tab)
2602                 return -1;
2603         if (tab->rational) {
2604                 int sgn = sign_of_min(tab, var);
2605                 if (sgn < -1)
2606                         return -1;
2607                 return sgn >= 0;
2608         } else {
2609                 int irred = isl_tab_min_at_most_neg_one(tab, var);
2610                 if (irred < 0)
2611                         return -1;
2612                 return !irred;
2613         }
2614 }
2615
2616 /* Check for (near) redundant constraints.
2617  * A constraint is redundant if it is non-negative and if
2618  * its minimal value (temporarily ignoring the non-negativity) is either
2619  *      - zero (in case of rational tableaus), or
2620  *      - strictly larger than -1 (in case of integer tableaus)
2621  *
2622  * We first mark all non-redundant and non-dead variables that
2623  * are not frozen and not obviously negatively unbounded.
2624  * Then we iterate over all marked variables if they can attain
2625  * any values smaller than zero or at most negative one.
2626  * If not, we mark the row as being redundant (assuming it hasn't
2627  * been detected as being obviously redundant in the mean time).
2628  */
2629 int isl_tab_detect_redundant(struct isl_tab *tab)
2630 {
2631         int i;
2632         unsigned n_marked;
2633
2634         if (!tab)
2635                 return -1;
2636         if (tab->empty)
2637                 return 0;
2638         if (tab->n_redundant == tab->n_row)
2639                 return 0;
2640
2641         n_marked = 0;
2642         for (i = tab->n_redundant; i < tab->n_row; ++i) {
2643                 struct isl_tab_var *var = isl_tab_var_from_row(tab, i);
2644                 var->marked = !var->frozen && var->is_nonneg;
2645                 if (var->marked)
2646                         n_marked++;
2647         }
2648         for (i = tab->n_dead; i < tab->n_col; ++i) {
2649                 struct isl_tab_var *var = var_from_col(tab, i);
2650                 var->marked = !var->frozen && var->is_nonneg &&
2651                         !min_is_manifestly_unbounded(tab, var);
2652                 if (var->marked)
2653                         n_marked++;
2654         }
2655         while (n_marked) {
2656                 struct isl_tab_var *var;
2657                 int red;
2658                 for (i = tab->n_redundant; i < tab->n_row; ++i) {
2659                         var = isl_tab_var_from_row(tab, i);
2660                         if (var->marked)
2661                                 break;
2662                 }
2663                 if (i == tab->n_row) {
2664                         for (i = tab->n_dead; i < tab->n_col; ++i) {
2665                                 var = var_from_col(tab, i);
2666                                 if (var->marked)
2667                                         break;
2668                         }
2669                         if (i == tab->n_col)
2670                                 break;
2671                 }
2672                 var->marked = 0;
2673                 n_marked--;
2674                 red = con_is_redundant(tab, var);
2675                 if (red < 0)
2676                         return -1;
2677                 if (red && !var->is_redundant)
2678                         if (isl_tab_mark_redundant(tab, var->index) < 0)
2679                                 return -1;
2680                 for (i = tab->n_dead; i < tab->n_col; ++i) {
2681                         var = var_from_col(tab, i);
2682                         if (!var->marked)
2683                                 continue;
2684                         if (!min_is_manifestly_unbounded(tab, var))
2685                                 continue;
2686                         var->marked = 0;
2687                         n_marked--;
2688                 }
2689         }
2690
2691         return 0;
2692 }
2693
2694 int isl_tab_is_equality(struct isl_tab *tab, int con)
2695 {
2696         int row;
2697         unsigned off;
2698
2699         if (!tab)
2700                 return -1;
2701         if (tab->con[con].is_zero)
2702                 return 1;
2703         if (tab->con[con].is_redundant)
2704                 return 0;
2705         if (!tab->con[con].is_row)
2706                 return tab->con[con].index < tab->n_dead;
2707
2708         row = tab->con[con].index;
2709
2710         off = 2 + tab->M;
2711         return isl_int_is_zero(tab->mat->row[row][1]) &&
2712                 isl_seq_first_non_zero(tab->mat->row[row] + 2 + tab->n_dead,
2713                                         tab->n_col - tab->n_dead) == -1;
2714 }
2715
2716 /* Return the minimial value of the affine expression "f" with denominator
2717  * "denom" in *opt, *opt_denom, assuming the tableau is not empty and
2718  * the expression cannot attain arbitrarily small values.
2719  * If opt_denom is NULL, then *opt is rounded up to the nearest integer.
2720  * The return value reflects the nature of the result (empty, unbounded,
2721  * minmimal value returned in *opt).
2722  */
2723 enum isl_lp_result isl_tab_min(struct isl_tab *tab,
2724         isl_int *f, isl_int denom, isl_int *opt, isl_int *opt_denom,
2725         unsigned flags)
2726 {
2727         int r;
2728         enum isl_lp_result res = isl_lp_ok;
2729         struct isl_tab_var *var;
2730         struct isl_tab_undo *snap;
2731
2732         if (!tab)
2733                 return isl_lp_error;
2734
2735         if (tab->empty)
2736                 return isl_lp_empty;
2737
2738         snap = isl_tab_snap(tab);
2739         r = isl_tab_add_row(tab, f);
2740         if (r < 0)
2741                 return isl_lp_error;
2742         var = &tab->con[r];
2743         isl_int_mul(tab->mat->row[var->index][0],
2744                     tab->mat->row[var->index][0], denom);
2745         for (;;) {
2746                 int row, col;
2747                 find_pivot(tab, var, var, -1, &row, &col);
2748                 if (row == var->index) {
2749                         res = isl_lp_unbounded;
2750                         break;
2751                 }
2752                 if (row == -1)
2753                         break;
2754                 if (isl_tab_pivot(tab, row, col) < 0)
2755                         return isl_lp_error;
2756         }
2757         if (ISL_FL_ISSET(flags, ISL_TAB_SAVE_DUAL)) {
2758                 int i;
2759
2760                 isl_vec_free(tab->dual);
2761                 tab->dual = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_con);
2762                 if (!tab->dual)
2763                         return isl_lp_error;
2764                 isl_int_set(tab->dual->el[0], tab->mat->row[var->index][0]);
2765                 for (i = 0; i < tab->n_con; ++i) {
2766                         int pos;
2767                         if (tab->con[i].is_row) {
2768                                 isl_int_set_si(tab->dual->el[1 + i], 0);
2769                                 continue;
2770                         }
2771                         pos = 2 + tab->M + tab->con[i].index;
2772                         if (tab->con[i].negated)
2773                                 isl_int_neg(tab->dual->el[1 + i],
2774                                             tab->mat->row[var->index][pos]);
2775                         else
2776                                 isl_int_set(tab->dual->el[1 + i],
2777                                             tab->mat->row[var->index][pos]);
2778                 }
2779         }
2780         if (opt && res == isl_lp_ok) {
2781                 if (opt_denom) {
2782                         isl_int_set(*opt, tab->mat->row[var->index][1]);
2783                         isl_int_set(*opt_denom, tab->mat->row[var->index][0]);
2784                 } else
2785                         isl_int_cdiv_q(*opt, tab->mat->row[var->index][1],
2786                                              tab->mat->row[var->index][0]);
2787         }
2788         if (isl_tab_rollback(tab, snap) < 0)
2789                 return isl_lp_error;
2790         return res;
2791 }
2792
2793 int isl_tab_is_redundant(struct isl_tab *tab, int con)
2794 {
2795         if (!tab)
2796                 return -1;
2797         if (tab->con[con].is_zero)
2798                 return 0;
2799         if (tab->con[con].is_redundant)
2800                 return 1;
2801         return tab->con[con].is_row && tab->con[con].index < tab->n_redundant;
2802 }
2803
2804 /* Take a snapshot of the tableau that can be restored by s call to
2805  * isl_tab_rollback.
2806  */
2807 struct isl_tab_undo *isl_tab_snap(struct isl_tab *tab)
2808 {
2809         if (!tab)
2810                 return NULL;
2811         tab->need_undo = 1;
2812         return tab->top;
2813 }
2814
2815 /* Undo the operation performed by isl_tab_relax.
2816  */
2817 static int unrelax(struct isl_tab *tab, struct isl_tab_var *var) WARN_UNUSED;
2818 static int unrelax(struct isl_tab *tab, struct isl_tab_var *var)
2819 {
2820         unsigned off = 2 + tab->M;
2821
2822         if (!var->is_row && !max_is_manifestly_unbounded(tab, var))
2823                 if (to_row(tab, var, 1) < 0)
2824                         return -1;
2825
2826         if (var->is_row) {
2827                 isl_int_sub(tab->mat->row[var->index][1],
2828                     tab->mat->row[var->index][1], tab->mat->row[var->index][0]);
2829                 if (var->is_nonneg) {
2830                         int sgn = restore_row(tab, var);
2831                         isl_assert(tab->mat->ctx, sgn >= 0, return -1);
2832                 }
2833         } else {
2834                 int i;
2835
2836                 for (i = 0; i < tab->n_row; ++i) {
2837                         if (isl_int_is_zero(tab->mat->row[i][off + var->index]))
2838                                 continue;
2839                         isl_int_add(tab->mat->row[i][1], tab->mat->row[i][1],
2840                             tab->mat->row[i][off + var->index]);
2841                 }
2842
2843         }
2844
2845         return 0;
2846 }
2847
2848 static int perform_undo_var(struct isl_tab *tab, struct isl_tab_undo *undo) WARN_UNUSED;
2849 static int perform_undo_var(struct isl_tab *tab, struct isl_tab_undo *undo)
2850 {
2851         struct isl_tab_var *var = var_from_index(tab, undo->u.var_index);
2852         switch(undo->type) {
2853         case isl_tab_undo_nonneg:
2854                 var->is_nonneg = 0;
2855                 break;
2856         case isl_tab_undo_redundant:
2857                 var->is_redundant = 0;
2858                 tab->n_redundant--;
2859                 restore_row(tab, isl_tab_var_from_row(tab, tab->n_redundant));
2860                 break;
2861         case isl_tab_undo_freeze:
2862                 var->frozen = 0;
2863                 break;
2864         case isl_tab_undo_zero:
2865                 var->is_zero = 0;
2866                 if (!var->is_row)
2867                         tab->n_dead--;
2868                 break;
2869         case isl_tab_undo_allocate:
2870                 if (undo->u.var_index >= 0) {
2871                         isl_assert(tab->mat->ctx, !var->is_row, return -1);
2872                         drop_col(tab, var->index);
2873                         break;
2874                 }
2875                 if (!var->is_row) {
2876                         if (!max_is_manifestly_unbounded(tab, var)) {
2877                                 if (to_row(tab, var, 1) < 0)
2878                                         return -1;
2879                         } else if (!min_is_manifestly_unbounded(tab, var)) {
2880                                 if (to_row(tab, var, -1) < 0)
2881                                         return -1;
2882                         } else
2883                                 if (to_row(tab, var, 0) < 0)
2884                                         return -1;
2885                 }
2886                 drop_row(tab, var->index);
2887                 break;
2888         case isl_tab_undo_relax:
2889                 return unrelax(tab, var);
2890         }
2891
2892         return 0;
2893 }
2894
2895 /* Restore the tableau to the state where the basic variables
2896  * are those in "col_var".
2897  * We first construct a list of variables that are currently in
2898  * the basis, but shouldn't.  Then we iterate over all variables
2899  * that should be in the basis and for each one that is currently
2900  * not in the basis, we exchange it with one of the elements of the
2901  * list constructed before.
2902  * We can always find an appropriate variable to pivot with because
2903  * the current basis is mapped to the old basis by a non-singular
2904  * matrix and so we can never end up with a zero row.
2905  */
2906 static int restore_basis(struct isl_tab *tab, int *col_var)
2907 {
2908         int i, j;
2909         int n_extra = 0;
2910         int *extra = NULL;      /* current columns that contain bad stuff */
2911         unsigned off = 2 + tab->M;
2912
2913         extra = isl_alloc_array(tab->mat->ctx, int, tab->n_col);
2914         if (!extra)
2915                 goto error;
2916         for (i = 0; i < tab->n_col; ++i) {
2917                 for (j = 0; j < tab->n_col; ++j)
2918                         if (tab->col_var[i] == col_var[j])
2919                                 break;
2920                 if (j < tab->n_col)
2921                         continue;
2922                 extra[n_extra++] = i;
2923         }
2924         for (i = 0; i < tab->n_col && n_extra > 0; ++i) {
2925                 struct isl_tab_var *var;
2926                 int row;
2927
2928                 for (j = 0; j < tab->n_col; ++j)
2929                         if (col_var[i] == tab->col_var[j])
2930                                 break;
2931                 if (j < tab->n_col)
2932                         continue;
2933                 var = var_from_index(tab, col_var[i]);
2934                 row = var->index;
2935                 for (j = 0; j < n_extra; ++j)
2936                         if (!isl_int_is_zero(tab->mat->row[row][off+extra[j]]))
2937                                 break;
2938                 isl_assert(tab->mat->ctx, j < n_extra, goto error);
2939                 if (isl_tab_pivot(tab, row, extra[j]) < 0)
2940                         goto error;
2941                 extra[j] = extra[--n_extra];
2942         }
2943
2944         free(extra);
2945         free(col_var);
2946         return 0;
2947 error:
2948         free(extra);
2949         free(col_var);
2950         return -1;
2951 }
2952
2953 /* Remove all samples with index n or greater, i.e., those samples
2954  * that were added since we saved this number of samples in
2955  * isl_tab_save_samples.
2956  */
2957 static void drop_samples_since(struct isl_tab *tab, int n)
2958 {
2959         int i;
2960
2961         for (i = tab->n_sample - 1; i >= 0 && tab->n_sample > n; --i) {
2962                 if (tab->sample_index[i] < n)
2963                         continue;
2964
2965                 if (i != tab->n_sample - 1) {
2966                         int t = tab->sample_index[tab->n_sample-1];
2967                         tab->sample_index[tab->n_sample-1] = tab->sample_index[i];
2968                         tab->sample_index[i] = t;
2969                         isl_mat_swap_rows(tab->samples, tab->n_sample-1, i);
2970                 }
2971                 tab->n_sample--;
2972         }
2973 }
2974
2975 static int perform_undo(struct isl_tab *tab, struct isl_tab_undo *undo) WARN_UNUSED;
2976 static int perform_undo(struct isl_tab *tab, struct isl_tab_undo *undo)
2977 {
2978         switch (undo->type) {
2979         case isl_tab_undo_empty:
2980                 tab->empty = 0;
2981                 break;
2982         case isl_tab_undo_nonneg:
2983         case isl_tab_undo_redundant:
2984         case isl_tab_undo_freeze:
2985         case isl_tab_undo_zero:
2986         case isl_tab_undo_allocate:
2987         case isl_tab_undo_relax:
2988                 return perform_undo_var(tab, undo);
2989         case isl_tab_undo_bmap_eq:
2990                 return isl_basic_map_free_equality(tab->bmap, 1);
2991         case isl_tab_undo_bmap_ineq:
2992                 return isl_basic_map_free_inequality(tab->bmap, 1);
2993         case isl_tab_undo_bmap_div:
2994                 if (isl_basic_map_free_div(tab->bmap, 1) < 0)
2995                         return -1;
2996                 if (tab->samples)
2997                         tab->samples->n_col--;
2998                 break;
2999         case isl_tab_undo_saved_basis:
3000                 if (restore_basis(tab, undo->u.col_var) < 0)
3001                         return -1;
3002                 break;
3003         case isl_tab_undo_drop_sample:
3004                 tab->n_outside--;
3005                 break;
3006         case isl_tab_undo_saved_samples:
3007                 drop_samples_since(tab, undo->u.n);
3008                 break;
3009         case isl_tab_undo_callback:
3010                 return undo->u.callback->run(undo->u.callback);
3011         default:
3012                 isl_assert(tab->mat->ctx, 0, return -1);
3013         }
3014         return 0;
3015 }
3016
3017 /* Return the tableau to the state it was in when the snapshot "snap"
3018  * was taken.
3019  */
3020 int isl_tab_rollback(struct isl_tab *tab, struct isl_tab_undo *snap)
3021 {
3022         struct isl_tab_undo *undo, *next;
3023
3024         if (!tab)
3025                 return -1;
3026
3027         tab->in_undo = 1;
3028         for (undo = tab->top; undo && undo != &tab->bottom; undo = next) {
3029                 next = undo->next;
3030                 if (undo == snap)
3031                         break;
3032                 if (perform_undo(tab, undo) < 0) {
3033                         tab->top = undo;
3034                         free_undo(tab);
3035                         tab->in_undo = 0;
3036                         return -1;
3037                 }
3038                 free(undo);
3039         }
3040         tab->in_undo = 0;
3041         tab->top = undo;
3042         if (!undo)
3043                 return -1;
3044         return 0;
3045 }
3046
3047 /* The given row "row" represents an inequality violated by all
3048  * points in the tableau.  Check for some special cases of such
3049  * separating constraints.
3050  * In particular, if the row has been reduced to the constant -1,
3051  * then we know the inequality is adjacent (but opposite) to
3052  * an equality in the tableau.
3053  * If the row has been reduced to r = -1 -r', with r' an inequality
3054  * of the tableau, then the inequality is adjacent (but opposite)
3055  * to the inequality r'.
3056  */
3057 static enum isl_ineq_type separation_type(struct isl_tab *tab, unsigned row)
3058 {
3059         int pos;
3060         unsigned off = 2 + tab->M;
3061
3062         if (tab->rational)
3063                 return isl_ineq_separate;
3064
3065         if (!isl_int_is_one(tab->mat->row[row][0]))
3066                 return isl_ineq_separate;
3067         if (!isl_int_is_negone(tab->mat->row[row][1]))
3068                 return isl_ineq_separate;
3069
3070         pos = isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
3071                                         tab->n_col - tab->n_dead);
3072         if (pos == -1)
3073                 return isl_ineq_adj_eq;
3074
3075         if (!isl_int_is_negone(tab->mat->row[row][off + tab->n_dead + pos]))
3076                 return isl_ineq_separate;
3077
3078         pos = isl_seq_first_non_zero(
3079                         tab->mat->row[row] + off + tab->n_dead + pos + 1,
3080                         tab->n_col - tab->n_dead - pos - 1);
3081
3082         return pos == -1 ? isl_ineq_adj_ineq : isl_ineq_separate;
3083 }
3084
3085 /* Check the effect of inequality "ineq" on the tableau "tab".
3086  * The result may be
3087  *      isl_ineq_redundant:     satisfied by all points in the tableau
3088  *      isl_ineq_separate:      satisfied by no point in the tableau
3089  *      isl_ineq_cut:           satisfied by some by not all points
3090  *      isl_ineq_adj_eq:        adjacent to an equality
3091  *      isl_ineq_adj_ineq:      adjacent to an inequality.
3092  */
3093 enum isl_ineq_type isl_tab_ineq_type(struct isl_tab *tab, isl_int *ineq)
3094 {
3095         enum isl_ineq_type type = isl_ineq_error;
3096         struct isl_tab_undo *snap = NULL;
3097         int con;
3098         int row;
3099
3100         if (!tab)
3101                 return isl_ineq_error;
3102
3103         if (isl_tab_extend_cons(tab, 1) < 0)
3104                 return isl_ineq_error;
3105
3106         snap = isl_tab_snap(tab);
3107
3108         con = isl_tab_add_row(tab, ineq);
3109         if (con < 0)
3110                 goto error;
3111
3112         row = tab->con[con].index;
3113         if (isl_tab_row_is_redundant(tab, row))
3114                 type = isl_ineq_redundant;
3115         else if (isl_int_is_neg(tab->mat->row[row][1]) &&
3116                  (tab->rational ||
3117                     isl_int_abs_ge(tab->mat->row[row][1],
3118                                    tab->mat->row[row][0]))) {
3119                 int nonneg = at_least_zero(tab, &tab->con[con]);
3120                 if (nonneg < 0)
3121                         goto error;
3122                 if (nonneg)
3123                         type = isl_ineq_cut;
3124                 else
3125                         type = separation_type(tab, row);
3126         } else {
3127                 int red = con_is_redundant(tab, &tab->con[con]);
3128                 if (red < 0)
3129                         goto error;
3130                 if (!red)
3131                         type = isl_ineq_cut;
3132                 else
3133                         type = isl_ineq_redundant;
3134         }
3135
3136         if (isl_tab_rollback(tab, snap))
3137                 return isl_ineq_error;
3138         return type;
3139 error:
3140         return isl_ineq_error;
3141 }
3142
3143 int isl_tab_track_bmap(struct isl_tab *tab, __isl_take isl_basic_map *bmap)
3144 {
3145         if (!tab || !bmap)
3146                 goto error;
3147
3148         isl_assert(tab->mat->ctx, tab->n_eq == bmap->n_eq, return -1);
3149         isl_assert(tab->mat->ctx,
3150                     tab->n_con == bmap->n_eq + bmap->n_ineq, return -1);
3151
3152         tab->bmap = bmap;
3153
3154         return 0;
3155 error:
3156         isl_basic_map_free(bmap);
3157         return -1;
3158 }
3159
3160 int isl_tab_track_bset(struct isl_tab *tab, __isl_take isl_basic_set *bset)
3161 {
3162         return isl_tab_track_bmap(tab, (isl_basic_map *)bset);
3163 }
3164
3165 __isl_keep isl_basic_set *isl_tab_peek_bset(struct isl_tab *tab)
3166 {
3167         if (!tab)
3168                 return NULL;
3169
3170         return (isl_basic_set *)tab->bmap;
3171 }
3172
3173 void isl_tab_dump(struct isl_tab *tab, FILE *out, int indent)
3174 {
3175         unsigned r, c;
3176         int i;
3177
3178         if (!tab) {
3179                 fprintf(out, "%*snull tab\n", indent, "");
3180                 return;
3181         }
3182         fprintf(out, "%*sn_redundant: %d, n_dead: %d", indent, "",
3183                 tab->n_redundant, tab->n_dead);
3184         if (tab->rational)
3185                 fprintf(out, ", rational");
3186         if (tab->empty)
3187                 fprintf(out, ", empty");
3188         fprintf(out, "\n");
3189         fprintf(out, "%*s[", indent, "");
3190         for (i = 0; i < tab->n_var; ++i) {
3191                 if (i)
3192                         fprintf(out, (i == tab->n_param ||
3193                                       i == tab->n_var - tab->n_div) ? "; "
3194                                                                     : ", ");
3195                 fprintf(out, "%c%d%s", tab->var[i].is_row ? 'r' : 'c',
3196                                         tab->var[i].index,
3197                                         tab->var[i].is_zero ? " [=0]" :
3198                                         tab->var[i].is_redundant ? " [R]" : "");
3199         }
3200         fprintf(out, "]\n");
3201         fprintf(out, "%*s[", indent, "");
3202         for (i = 0; i < tab->n_con; ++i) {
3203                 if (i)
3204                         fprintf(out, ", ");
3205                 fprintf(out, "%c%d%s", tab->con[i].is_row ? 'r' : 'c',
3206                                         tab->con[i].index,
3207                                         tab->con[i].is_zero ? " [=0]" :
3208                                         tab->con[i].is_redundant ? " [R]" : "");
3209         }
3210         fprintf(out, "]\n");
3211         fprintf(out, "%*s[", indent, "");
3212         for (i = 0; i < tab->n_row; ++i) {
3213                 const char *sign = "";
3214                 if (i)
3215                         fprintf(out, ", ");
3216                 if (tab->row_sign) {
3217                         if (tab->row_sign[i] == isl_tab_row_unknown)
3218                                 sign = "?";
3219                         else if (tab->row_sign[i] == isl_tab_row_neg)
3220                                 sign = "-";
3221                         else if (tab->row_sign[i] == isl_tab_row_pos)
3222                                 sign = "+";
3223                         else
3224                                 sign = "+-";
3225                 }
3226                 fprintf(out, "r%d: %d%s%s", i, tab->row_var[i],
3227                     isl_tab_var_from_row(tab, i)->is_nonneg ? " [>=0]" : "", sign);
3228         }
3229         fprintf(out, "]\n");
3230         fprintf(out, "%*s[", indent, "");
3231         for (i = 0; i < tab->n_col; ++i) {
3232                 if (i)
3233                         fprintf(out, ", ");
3234                 fprintf(out, "c%d: %d%s", i, tab->col_var[i],
3235                     var_from_col(tab, i)->is_nonneg ? " [>=0]" : "");
3236         }
3237         fprintf(out, "]\n");
3238         r = tab->mat->n_row;
3239         tab->mat->n_row = tab->n_row;
3240         c = tab->mat->n_col;
3241         tab->mat->n_col = 2 + tab->M + tab->n_col;
3242         isl_mat_dump(tab->mat, out, indent);
3243         tab->mat->n_row = r;
3244         tab->mat->n_col = c;
3245         if (tab->bmap)
3246                 isl_basic_map_dump(tab->bmap, out, indent);
3247 }