isl_tab.c

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