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 /* Check whether the div described by "div" is obviously non-negative.
2064  * If we are using a big parameter, then we will encode the div
2065  * as div' = M + div, which is always non-negative.
2066  * Otherwise, we check whether div is a non-negative affine combination
2067  * of non-negative variables.
2068  */
2069 static int div_is_nonneg(struct isl_tab *tab, __isl_keep isl_vec *div)
2070 {
2071         int i;
2072
2073         if (tab->M)
2074                 return 1;
2075
2076         if (isl_int_is_neg(div->el[1]))
2077                 return 0;
2078
2079         for (i = 0; i < tab->n_var; ++i) {
2080                 if (isl_int_is_neg(div->el[2 + i]))
2081                         return 0;
2082                 if (isl_int_is_zero(div->el[2 + i]))
2083                         continue;
2084                 if (!tab->var[i].is_nonneg)
2085                         return 0;
2086         }
2087
2088         return 1;
2089 }
2090
2091 /* Add an extra div, prescribed by "div" to the tableau and
2092  * the associated bmap (which is assumed to be non-NULL).
2093  *
2094  * If add_ineq is not NULL, then this function is used instead
2095  * of isl_tab_add_ineq to add the div constraints.
2096  * This complication is needed because the code in isl_tab_pip
2097  * wants to perform some extra processing when an inequality
2098  * is added to the tableau.
2099  */
2100 int isl_tab_add_div(struct isl_tab *tab, __isl_keep isl_vec *div,
2101         int (*add_ineq)(void *user, isl_int *), void *user)
2102 {
2103         int r;
2104         int k;
2105         int nonneg;
2106
2107         if (!tab || !div)
2108                 return -1;
2109
2110         isl_assert(tab->mat->ctx, tab->bmap, return -1);
2111
2112         nonneg = div_is_nonneg(tab, div);
2113
2114         if (isl_tab_extend_cons(tab, 3) < 0)
2115                 return -1;
2116         if (isl_tab_extend_vars(tab, 1) < 0)
2117                 return -1;
2118         r = isl_tab_allocate_var(tab);
2119         if (r < 0)
2120                 return -1;
2121
2122         if (nonneg)
2123                 tab->var[r].is_nonneg = 1;
2124
2125         tab->bmap = isl_basic_map_extend_dim(tab->bmap,
2126                 isl_basic_map_get_dim(tab->bmap), 1, 0, 2);
2127         k = isl_basic_map_alloc_div(tab->bmap);
2128         if (k < 0)
2129                 return -1;
2130         isl_seq_cpy(tab->bmap->div[k], div->el, div->size);
2131         if (isl_tab_push(tab, isl_tab_undo_bmap_div) < 0)
2132                 return -1;
2133
2134         if (add_div_constraints(tab, k, add_ineq, user) < 0)
2135                 return -1;
2136
2137         return r;
2138 }
2139
2140 struct isl_tab *isl_tab_from_basic_map(struct isl_basic_map *bmap)
2141 {
2142         int i;
2143         struct isl_tab *tab;
2144
2145         if (!bmap)
2146                 return NULL;
2147         tab = isl_tab_alloc(bmap->ctx,
2148                             isl_basic_map_total_dim(bmap) + bmap->n_ineq + 1,
2149                             isl_basic_map_total_dim(bmap), 0);
2150         if (!tab)
2151                 return NULL;
2152         tab->rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
2153         if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY)) {
2154                 if (isl_tab_mark_empty(tab) < 0)
2155                         goto error;
2156                 return tab;
2157         }
2158         for (i = 0; i < bmap->n_eq; ++i) {
2159                 tab = add_eq(tab, bmap->eq[i]);
2160                 if (!tab)
2161                         return tab;
2162         }
2163         for (i = 0; i < bmap->n_ineq; ++i) {
2164                 if (isl_tab_add_ineq(tab, bmap->ineq[i]) < 0)
2165                         goto error;
2166                 if (tab->empty)
2167                         return tab;
2168         }
2169         return tab;
2170 error:
2171         isl_tab_free(tab);
2172         return NULL;
2173 }
2174
2175 struct isl_tab *isl_tab_from_basic_set(struct isl_basic_set *bset)
2176 {
2177         return isl_tab_from_basic_map((struct isl_basic_map *)bset);
2178 }
2179
2180 /* Construct a tableau corresponding to the recession cone of "bset".
2181  */
2182 struct isl_tab *isl_tab_from_recession_cone(__isl_keep isl_basic_set *bset,
2183         int parametric)
2184 {
2185         isl_int cst;
2186         int i;
2187         struct isl_tab *tab;
2188         unsigned offset = 0;
2189
2190         if (!bset)
2191                 return NULL;
2192         if (parametric)
2193                 offset = isl_basic_set_dim(bset, isl_dim_param);
2194         tab = isl_tab_alloc(bset->ctx, bset->n_eq + bset->n_ineq,
2195                                 isl_basic_set_total_dim(bset) - offset, 0);
2196         if (!tab)
2197                 return NULL;
2198         tab->rational = ISL_F_ISSET(bset, ISL_BASIC_SET_RATIONAL);
2199         tab->cone = 1;
2200
2201         isl_int_init(cst);
2202         for (i = 0; i < bset->n_eq; ++i) {
2203                 isl_int_swap(bset->eq[i][offset], cst);
2204                 if (offset > 0) {
2205                         if (isl_tab_add_eq(tab, bset->eq[i] + offset) < 0)
2206                                 goto error;
2207                 } else
2208                         tab = add_eq(tab, bset->eq[i]);
2209                 isl_int_swap(bset->eq[i][offset], cst);
2210                 if (!tab)
2211                         goto done;
2212         }
2213         for (i = 0; i < bset->n_ineq; ++i) {
2214                 int r;
2215                 isl_int_swap(bset->ineq[i][offset], cst);
2216                 r = isl_tab_add_row(tab, bset->ineq[i] + offset);
2217                 isl_int_swap(bset->ineq[i][offset], cst);
2218                 if (r < 0)
2219                         goto error;
2220                 tab->con[r].is_nonneg = 1;
2221                 if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
2222                         goto error;
2223         }
2224 done:
2225         isl_int_clear(cst);
2226         return tab;
2227 error:
2228         isl_int_clear(cst);
2229         isl_tab_free(tab);
2230         return NULL;
2231 }
2232
2233 /* Assuming "tab" is the tableau of a cone, check if the cone is
2234  * bounded, i.e., if it is empty or only contains the origin.
2235  */
2236 int isl_tab_cone_is_bounded(struct isl_tab *tab)
2237 {
2238         int i;
2239
2240         if (!tab)
2241                 return -1;
2242         if (tab->empty)
2243                 return 1;
2244         if (tab->n_dead == tab->n_col)
2245                 return 1;
2246
2247         for (;;) {
2248                 for (i = tab->n_redundant; i < tab->n_row; ++i) {
2249                         struct isl_tab_var *var;
2250                         int sgn;
2251                         var = isl_tab_var_from_row(tab, i);
2252                         if (!var->is_nonneg)
2253                                 continue;
2254                         sgn = sign_of_max(tab, var);
2255                         if (sgn < -1)
2256                                 return -1;
2257                         if (sgn != 0)
2258                                 return 0;
2259                         if (close_row(tab, var) < 0)
2260                                 return -1;
2261                         break;
2262                 }
2263                 if (tab->n_dead == tab->n_col)
2264                         return 1;
2265                 if (i == tab->n_row)
2266                         return 0;
2267         }
2268 }
2269
2270 int isl_tab_sample_is_integer(struct isl_tab *tab)
2271 {
2272         int i;
2273
2274         if (!tab)
2275                 return -1;
2276
2277         for (i = 0; i < tab->n_var; ++i) {
2278                 int row;
2279                 if (!tab->var[i].is_row)
2280                         continue;
2281                 row = tab->var[i].index;
2282                 if (!isl_int_is_divisible_by(tab->mat->row[row][1],
2283                                                 tab->mat->row[row][0]))
2284                         return 0;
2285         }
2286         return 1;
2287 }
2288
2289 static struct isl_vec *extract_integer_sample(struct isl_tab *tab)
2290 {
2291         int i;
2292         struct isl_vec *vec;
2293
2294         vec = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
2295         if (!vec)
2296                 return NULL;
2297
2298         isl_int_set_si(vec->block.data[0], 1);
2299         for (i = 0; i < tab->n_var; ++i) {
2300                 if (!tab->var[i].is_row)
2301                         isl_int_set_si(vec->block.data[1 + i], 0);
2302                 else {
2303                         int row = tab->var[i].index;
2304                         isl_int_divexact(vec->block.data[1 + i],
2305                                 tab->mat->row[row][1], tab->mat->row[row][0]);
2306                 }
2307         }
2308
2309         return vec;
2310 }
2311
2312 struct isl_vec *isl_tab_get_sample_value(struct isl_tab *tab)
2313 {
2314         int i;
2315         struct isl_vec *vec;
2316         isl_int m;
2317
2318         if (!tab)
2319                 return NULL;
2320
2321         vec = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
2322         if (!vec)
2323                 return NULL;
2324
2325         isl_int_init(m);
2326
2327         isl_int_set_si(vec->block.data[0], 1);
2328         for (i = 0; i < tab->n_var; ++i) {
2329                 int row;
2330                 if (!tab->var[i].is_row) {
2331                         isl_int_set_si(vec->block.data[1 + i], 0);
2332                         continue;
2333                 }
2334                 row = tab->var[i].index;
2335                 isl_int_gcd(m, vec->block.data[0], tab->mat->row[row][0]);
2336                 isl_int_divexact(m, tab->mat->row[row][0], m);
2337                 isl_seq_scale(vec->block.data, vec->block.data, m, 1 + i);
2338                 isl_int_divexact(m, vec->block.data[0], tab->mat->row[row][0]);
2339                 isl_int_mul(vec->block.data[1 + i], m, tab->mat->row[row][1]);
2340         }
2341         vec = isl_vec_normalize(vec);
2342
2343         isl_int_clear(m);
2344         return vec;
2345 }
2346
2347 /* Update "bmap" based on the results of the tableau "tab".
2348  * In particular, implicit equalities are made explicit, redundant constraints
2349  * are removed and if the sample value happens to be integer, it is stored
2350  * in "bmap" (unless "bmap" already had an integer sample).
2351  *
2352  * The tableau is assumed to have been created from "bmap" using
2353  * isl_tab_from_basic_map.
2354  */
2355 struct isl_basic_map *isl_basic_map_update_from_tab(struct isl_basic_map *bmap,
2356         struct isl_tab *tab)
2357 {
2358         int i;
2359         unsigned n_eq;
2360
2361         if (!bmap)
2362                 return NULL;
2363         if (!tab)
2364                 return bmap;
2365
2366         n_eq = tab->n_eq;
2367         if (tab->empty)
2368                 bmap = isl_basic_map_set_to_empty(bmap);
2369         else
2370                 for (i = bmap->n_ineq - 1; i >= 0; --i) {
2371                         if (isl_tab_is_equality(tab, n_eq + i))
2372                                 isl_basic_map_inequality_to_equality(bmap, i);
2373                         else if (isl_tab_is_redundant(tab, n_eq + i))
2374                                 isl_basic_map_drop_inequality(bmap, i);
2375                 }
2376         if (bmap->n_eq != n_eq)
2377                 isl_basic_map_gauss(bmap, NULL);
2378         if (!tab->rational &&
2379             !bmap->sample && isl_tab_sample_is_integer(tab))
2380                 bmap->sample = extract_integer_sample(tab);
2381         return bmap;
2382 }
2383
2384 struct isl_basic_set *isl_basic_set_update_from_tab(struct isl_basic_set *bset,
2385         struct isl_tab *tab)
2386 {
2387         return (struct isl_basic_set *)isl_basic_map_update_from_tab(
2388                 (struct isl_basic_map *)bset, tab);
2389 }
2390
2391 /* Given a non-negative variable "var", add a new non-negative variable
2392  * that is the opposite of "var", ensuring that var can only attain the
2393  * value zero.
2394  * If var = n/d is a row variable, then the new variable = -n/d.
2395  * If var is a column variables, then the new variable = -var.
2396  * If the new variable cannot attain non-negative values, then
2397  * the resulting tableau is empty.
2398  * Otherwise, we know the value will be zero and we close the row.
2399  */
2400 static int cut_to_hyperplane(struct isl_tab *tab, struct isl_tab_var *var)
2401 {
2402         unsigned r;
2403         isl_int *row;
2404         int sgn;
2405         unsigned off = 2 + tab->M;
2406
2407         if (var->is_zero)
2408                 return 0;
2409         isl_assert(tab->mat->ctx, !var->is_redundant, return -1);
2410         isl_assert(tab->mat->ctx, var->is_nonneg, return -1);
2411
2412         if (isl_tab_extend_cons(tab, 1) < 0)
2413                 return -1;
2414
2415         r = tab->n_con;
2416         tab->con[r].index = tab->n_row;
2417         tab->con[r].is_row = 1;
2418         tab->con[r].is_nonneg = 0;
2419         tab->con[r].is_zero = 0;
2420         tab->con[r].is_redundant = 0;
2421         tab->con[r].frozen = 0;
2422         tab->con[r].negated = 0;
2423         tab->row_var[tab->n_row] = ~r;
2424         row = tab->mat->row[tab->n_row];
2425
2426         if (var->is_row) {
2427                 isl_int_set(row[0], tab->mat->row[var->index][0]);
2428                 isl_seq_neg(row + 1,
2429                             tab->mat->row[var->index] + 1, 1 + tab->n_col);
2430         } else {
2431                 isl_int_set_si(row[0], 1);
2432                 isl_seq_clr(row + 1, 1 + tab->n_col);
2433                 isl_int_set_si(row[off + var->index], -1);
2434         }
2435
2436         tab->n_row++;
2437         tab->n_con++;
2438         if (isl_tab_push_var(tab, isl_tab_undo_allocate, &tab->con[r]) < 0)
2439                 return -1;
2440
2441         sgn = sign_of_max(tab, &tab->con[r]);
2442         if (sgn < -1)
2443                 return -1;
2444         if (sgn < 0) {
2445                 if (isl_tab_mark_empty(tab) < 0)
2446                         return -1;
2447                 return 0;
2448         }
2449         tab->con[r].is_nonneg = 1;
2450         if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
2451                 return -1;
2452         /* sgn == 0 */
2453         if (close_row(tab, &tab->con[r]) < 0)
2454                 return -1;
2455
2456         return 0;
2457 }
2458
2459 /* Given a tableau "tab" and an inequality constraint "con" of the tableau,
2460  * relax the inequality by one.  That is, the inequality r >= 0 is replaced
2461  * by r' = r + 1 >= 0.
2462  * If r is a row variable, we simply increase the constant term by one
2463  * (taking into account the denominator).
2464  * If r is a column variable, then we need to modify each row that
2465  * refers to r = r' - 1 by substituting this equality, effectively
2466  * subtracting the coefficient of the column from the constant.
2467  * We should only do this if the minimum is manifestly unbounded,
2468  * however.  Otherwise, we may end up with negative sample values
2469  * for non-negative variables.
2470  * So, if r is a column variable with a minimum that is not
2471  * manifestly unbounded, then we need to move it to a row.
2472  * However, the sample value of this row may be negative,
2473  * even after the relaxation, so we need to restore it.
2474  * We therefore prefer to pivot a column up to a row, if possible.
2475  */
2476 struct isl_tab *isl_tab_relax(struct isl_tab *tab, int con)
2477 {
2478         struct isl_tab_var *var;
2479         unsigned off = 2 + tab->M;
2480
2481         if (!tab)
2482                 return NULL;
2483
2484         var = &tab->con[con];
2485
2486         if (!var->is_row && !max_is_manifestly_unbounded(tab, var))
2487                 if (to_row(tab, var, 1) < 0)
2488                         goto error;
2489         if (!var->is_row && !min_is_manifestly_unbounded(tab, var))
2490                 if (to_row(tab, var, -1) < 0)
2491                         goto error;
2492
2493         if (var->is_row) {
2494                 isl_int_add(tab->mat->row[var->index][1],
2495                     tab->mat->row[var->index][1], tab->mat->row[var->index][0]);
2496                 if (restore_row(tab, var) < 0)
2497                         goto error;
2498         } else {
2499                 int i;
2500
2501                 for (i = 0; i < tab->n_row; ++i) {
2502                         if (isl_int_is_zero(tab->mat->row[i][off + var->index]))
2503                                 continue;
2504                         isl_int_sub(tab->mat->row[i][1], tab->mat->row[i][1],
2505                             tab->mat->row[i][off + var->index]);
2506                 }
2507
2508         }
2509
2510         if (isl_tab_push_var(tab, isl_tab_undo_relax, var) < 0)
2511                 goto error;
2512
2513         return tab;
2514 error:
2515         isl_tab_free(tab);
2516         return NULL;
2517 }
2518
2519 int isl_tab_select_facet(struct isl_tab *tab, int con)
2520 {
2521         if (!tab)
2522                 return -1;
2523
2524         return cut_to_hyperplane(tab, &tab->con[con]);
2525 }
2526
2527 static int may_be_equality(struct isl_tab *tab, int row)
2528 {
2529         unsigned off = 2 + tab->M;
2530         return tab->rational ? isl_int_is_zero(tab->mat->row[row][1])
2531                              : isl_int_lt(tab->mat->row[row][1],
2532                                             tab->mat->row[row][0]);
2533 }
2534
2535 /* Check for (near) equalities among the constraints.
2536  * A constraint is an equality if it is non-negative and if
2537  * its maximal value is either
2538  *      - zero (in case of rational tableaus), or
2539  *      - strictly less than 1 (in case of integer tableaus)
2540  *
2541  * We first mark all non-redundant and non-dead variables that
2542  * are not frozen and not obviously not an equality.
2543  * Then we iterate over all marked variables if they can attain
2544  * any values larger than zero or at least one.
2545  * If the maximal value is zero, we mark any column variables
2546  * that appear in the row as being zero and mark the row as being redundant.
2547  * Otherwise, if the maximal value is strictly less than one (and the
2548  * tableau is integer), then we restrict the value to being zero
2549  * by adding an opposite non-negative variable.
2550  */
2551 int isl_tab_detect_implicit_equalities(struct isl_tab *tab)
2552 {
2553         int i;
2554         unsigned n_marked;
2555
2556         if (!tab)
2557                 return -1;
2558         if (tab->empty)
2559                 return 0;
2560         if (tab->n_dead == tab->n_col)
2561                 return 0;
2562
2563         n_marked = 0;
2564         for (i = tab->n_redundant; i < tab->n_row; ++i) {
2565                 struct isl_tab_var *var = isl_tab_var_from_row(tab, i);
2566                 var->marked = !var->frozen && var->is_nonneg &&
2567                         may_be_equality(tab, i);
2568                 if (var->marked)
2569                         n_marked++;
2570         }
2571         for (i = tab->n_dead; i < tab->n_col; ++i) {
2572                 struct isl_tab_var *var = var_from_col(tab, i);
2573                 var->marked = !var->frozen && var->is_nonneg;
2574                 if (var->marked)
2575                         n_marked++;
2576         }
2577         while (n_marked) {
2578                 struct isl_tab_var *var;
2579                 int sgn;
2580                 for (i = tab->n_redundant; i < tab->n_row; ++i) {
2581                         var = isl_tab_var_from_row(tab, i);
2582                         if (var->marked)
2583                                 break;
2584                 }
2585                 if (i == tab->n_row) {
2586                         for (i = tab->n_dead; i < tab->n_col; ++i) {
2587                                 var = var_from_col(tab, i);
2588                                 if (var->marked)
2589                                         break;
2590                         }
2591                         if (i == tab->n_col)
2592                                 break;
2593                 }
2594                 var->marked = 0;
2595                 n_marked--;
2596                 sgn = sign_of_max(tab, var);
2597                 if (sgn < 0)
2598                         return -1;
2599                 if (sgn == 0) {
2600                         if (close_row(tab, var) < 0)
2601                                 return -1;
2602                 } else if (!tab->rational && !at_least_one(tab, var)) {
2603                         if (cut_to_hyperplane(tab, var) < 0)
2604                                 return -1;
2605                         return isl_tab_detect_implicit_equalities(tab);
2606                 }
2607                 for (i = tab->n_redundant; i < tab->n_row; ++i) {
2608                         var = isl_tab_var_from_row(tab, i);
2609                         if (!var->marked)
2610                                 continue;
2611                         if (may_be_equality(tab, i))
2612                                 continue;
2613                         var->marked = 0;
2614                         n_marked--;
2615                 }
2616         }
2617
2618         return 0;
2619 }
2620
2621 static int con_is_redundant(struct isl_tab *tab, struct isl_tab_var *var)
2622 {
2623         if (!tab)
2624                 return -1;
2625         if (tab->rational) {
2626                 int sgn = sign_of_min(tab, var);
2627                 if (sgn < -1)
2628                         return -1;
2629                 return sgn >= 0;
2630         } else {
2631                 int irred = isl_tab_min_at_most_neg_one(tab, var);
2632                 if (irred < 0)
2633                         return -1;
2634                 return !irred;
2635         }
2636 }
2637
2638 /* Check for (near) redundant constraints.
2639  * A constraint is redundant if it is non-negative and if
2640  * its minimal value (temporarily ignoring the non-negativity) is either
2641  *      - zero (in case of rational tableaus), or
2642  *      - strictly larger than -1 (in case of integer tableaus)
2643  *
2644  * We first mark all non-redundant and non-dead variables that
2645  * are not frozen and not obviously negatively unbounded.
2646  * Then we iterate over all marked variables if they can attain
2647  * any values smaller than zero or at most negative one.
2648  * If not, we mark the row as being redundant (assuming it hasn't
2649  * been detected as being obviously redundant in the mean time).
2650  */
2651 int isl_tab_detect_redundant(struct isl_tab *tab)
2652 {
2653         int i;
2654         unsigned n_marked;
2655
2656         if (!tab)
2657                 return -1;
2658         if (tab->empty)
2659                 return 0;
2660         if (tab->n_redundant == tab->n_row)
2661                 return 0;
2662
2663         n_marked = 0;
2664         for (i = tab->n_redundant; i < tab->n_row; ++i) {
2665                 struct isl_tab_var *var = isl_tab_var_from_row(tab, i);
2666                 var->marked = !var->frozen && var->is_nonneg;
2667                 if (var->marked)
2668                         n_marked++;
2669         }
2670         for (i = tab->n_dead; i < tab->n_col; ++i) {
2671                 struct isl_tab_var *var = var_from_col(tab, i);
2672                 var->marked = !var->frozen && var->is_nonneg &&
2673                         !min_is_manifestly_unbounded(tab, var);
2674                 if (var->marked)
2675                         n_marked++;
2676         }
2677         while (n_marked) {
2678                 struct isl_tab_var *var;
2679                 int red;
2680                 for (i = tab->n_redundant; i < tab->n_row; ++i) {
2681                         var = isl_tab_var_from_row(tab, i);
2682                         if (var->marked)
2683                                 break;
2684                 }
2685                 if (i == tab->n_row) {
2686                         for (i = tab->n_dead; i < tab->n_col; ++i) {
2687                                 var = var_from_col(tab, i);
2688                                 if (var->marked)
2689                                         break;
2690                         }
2691                         if (i == tab->n_col)
2692                                 break;
2693                 }
2694                 var->marked = 0;
2695                 n_marked--;
2696                 red = con_is_redundant(tab, var);
2697                 if (red < 0)
2698                         return -1;
2699                 if (red && !var->is_redundant)
2700                         if (isl_tab_mark_redundant(tab, var->index) < 0)
2701                                 return -1;
2702                 for (i = tab->n_dead; i < tab->n_col; ++i) {
2703                         var = var_from_col(tab, i);
2704                         if (!var->marked)
2705                                 continue;
2706                         if (!min_is_manifestly_unbounded(tab, var))
2707                                 continue;
2708                         var->marked = 0;
2709                         n_marked--;
2710                 }
2711         }
2712
2713         return 0;
2714 }
2715
2716 int isl_tab_is_equality(struct isl_tab *tab, int con)
2717 {
2718         int row;
2719         unsigned off;
2720
2721         if (!tab)
2722                 return -1;
2723         if (tab->con[con].is_zero)
2724                 return 1;
2725         if (tab->con[con].is_redundant)
2726                 return 0;
2727         if (!tab->con[con].is_row)
2728                 return tab->con[con].index < tab->n_dead;
2729
2730         row = tab->con[con].index;
2731
2732         off = 2 + tab->M;
2733         return isl_int_is_zero(tab->mat->row[row][1]) &&
2734                 (!tab->M || isl_int_is_zero(tab->mat->row[row][2])) &&
2735                 isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
2736                                         tab->n_col - tab->n_dead) == -1;
2737 }
2738
2739 /* Return the minimal value of the affine expression "f" with denominator
2740  * "denom" in *opt, *opt_denom, assuming the tableau is not empty and
2741  * the expression cannot attain arbitrarily small values.
2742  * If opt_denom is NULL, then *opt is rounded up to the nearest integer.
2743  * The return value reflects the nature of the result (empty, unbounded,
2744  * minimal value returned in *opt).
2745  */
2746 enum isl_lp_result isl_tab_min(struct isl_tab *tab,
2747         isl_int *f, isl_int denom, isl_int *opt, isl_int *opt_denom,
2748         unsigned flags)
2749 {
2750         int r;
2751         enum isl_lp_result res = isl_lp_ok;
2752         struct isl_tab_var *var;
2753         struct isl_tab_undo *snap;
2754
2755         if (!tab)
2756                 return isl_lp_error;
2757
2758         if (tab->empty)
2759                 return isl_lp_empty;
2760
2761         snap = isl_tab_snap(tab);
2762         r = isl_tab_add_row(tab, f);
2763         if (r < 0)
2764                 return isl_lp_error;
2765         var = &tab->con[r];
2766         isl_int_mul(tab->mat->row[var->index][0],
2767                     tab->mat->row[var->index][0], denom);
2768         for (;;) {
2769                 int row, col;
2770                 find_pivot(tab, var, var, -1, &row, &col);
2771                 if (row == var->index) {
2772                         res = isl_lp_unbounded;
2773                         break;
2774                 }
2775                 if (row == -1)
2776                         break;
2777                 if (isl_tab_pivot(tab, row, col) < 0)
2778                         return isl_lp_error;
2779         }
2780         if (ISL_FL_ISSET(flags, ISL_TAB_SAVE_DUAL)) {
2781                 int i;
2782
2783                 isl_vec_free(tab->dual);
2784                 tab->dual = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_con);
2785                 if (!tab->dual)
2786                         return isl_lp_error;
2787                 isl_int_set(tab->dual->el[0], tab->mat->row[var->index][0]);
2788                 for (i = 0; i < tab->n_con; ++i) {
2789                         int pos;
2790                         if (tab->con[i].is_row) {
2791                                 isl_int_set_si(tab->dual->el[1 + i], 0);
2792                                 continue;
2793                         }
2794                         pos = 2 + tab->M + tab->con[i].index;
2795                         if (tab->con[i].negated)
2796                                 isl_int_neg(tab->dual->el[1 + i],
2797                                             tab->mat->row[var->index][pos]);
2798                         else
2799                                 isl_int_set(tab->dual->el[1 + i],
2800                                             tab->mat->row[var->index][pos]);
2801                 }
2802         }
2803         if (opt && res == isl_lp_ok) {
2804                 if (opt_denom) {
2805                         isl_int_set(*opt, tab->mat->row[var->index][1]);
2806                         isl_int_set(*opt_denom, tab->mat->row[var->index][0]);
2807                 } else
2808                         isl_int_cdiv_q(*opt, tab->mat->row[var->index][1],
2809                                              tab->mat->row[var->index][0]);
2810         }
2811         if (isl_tab_rollback(tab, snap) < 0)
2812                 return isl_lp_error;
2813         return res;
2814 }
2815
2816 int isl_tab_is_redundant(struct isl_tab *tab, int con)
2817 {
2818         if (!tab)
2819                 return -1;
2820         if (tab->con[con].is_zero)
2821                 return 0;
2822         if (tab->con[con].is_redundant)
2823                 return 1;
2824         return tab->con[con].is_row && tab->con[con].index < tab->n_redundant;
2825 }
2826
2827 /* Take a snapshot of the tableau that can be restored by s call to
2828  * isl_tab_rollback.
2829  */
2830 struct isl_tab_undo *isl_tab_snap(struct isl_tab *tab)
2831 {
2832         if (!tab)
2833                 return NULL;
2834         tab->need_undo = 1;
2835         return tab->top;
2836 }
2837
2838 /* Undo the operation performed by isl_tab_relax.
2839  */
2840 static int unrelax(struct isl_tab *tab, struct isl_tab_var *var) WARN_UNUSED;
2841 static int unrelax(struct isl_tab *tab, struct isl_tab_var *var)
2842 {
2843         unsigned off = 2 + tab->M;
2844
2845         if (!var->is_row && !max_is_manifestly_unbounded(tab, var))
2846                 if (to_row(tab, var, 1) < 0)
2847                         return -1;
2848
2849         if (var->is_row) {
2850                 isl_int_sub(tab->mat->row[var->index][1],
2851                     tab->mat->row[var->index][1], tab->mat->row[var->index][0]);
2852                 if (var->is_nonneg) {
2853                         int sgn = restore_row(tab, var);
2854                         isl_assert(tab->mat->ctx, sgn >= 0, return -1);
2855                 }
2856         } else {
2857                 int i;
2858
2859                 for (i = 0; i < tab->n_row; ++i) {
2860                         if (isl_int_is_zero(tab->mat->row[i][off + var->index]))
2861                                 continue;
2862                         isl_int_add(tab->mat->row[i][1], tab->mat->row[i][1],
2863                             tab->mat->row[i][off + var->index]);
2864                 }
2865
2866         }
2867
2868         return 0;
2869 }
2870
2871 static int perform_undo_var(struct isl_tab *tab, struct isl_tab_undo *undo) WARN_UNUSED;
2872 static int perform_undo_var(struct isl_tab *tab, struct isl_tab_undo *undo)
2873 {
2874         struct isl_tab_var *var = var_from_index(tab, undo->u.var_index);
2875         switch(undo->type) {
2876         case isl_tab_undo_nonneg:
2877                 var->is_nonneg = 0;
2878                 break;
2879         case isl_tab_undo_redundant:
2880                 var->is_redundant = 0;
2881                 tab->n_redundant--;
2882                 restore_row(tab, isl_tab_var_from_row(tab, tab->n_redundant));
2883                 break;
2884         case isl_tab_undo_freeze:
2885                 var->frozen = 0;
2886                 break;
2887         case isl_tab_undo_zero:
2888                 var->is_zero = 0;
2889                 if (!var->is_row)
2890                         tab->n_dead--;
2891                 break;
2892         case isl_tab_undo_allocate:
2893                 if (undo->u.var_index >= 0) {
2894                         isl_assert(tab->mat->ctx, !var->is_row, return -1);
2895                         drop_col(tab, var->index);
2896                         break;
2897                 }
2898                 if (!var->is_row) {
2899                         if (!max_is_manifestly_unbounded(tab, var)) {
2900                                 if (to_row(tab, var, 1) < 0)
2901                                         return -1;
2902                         } else if (!min_is_manifestly_unbounded(tab, var)) {
2903                                 if (to_row(tab, var, -1) < 0)
2904                                         return -1;
2905                         } else
2906                                 if (to_row(tab, var, 0) < 0)
2907                                         return -1;
2908                 }
2909                 drop_row(tab, var->index);
2910                 break;
2911         case isl_tab_undo_relax:
2912                 return unrelax(tab, var);
2913         }
2914
2915         return 0;
2916 }
2917
2918 /* Restore the tableau to the state where the basic variables
2919  * are those in "col_var".
2920  * We first construct a list of variables that are currently in
2921  * the basis, but shouldn't.  Then we iterate over all variables
2922  * that should be in the basis and for each one that is currently
2923  * not in the basis, we exchange it with one of the elements of the
2924  * list constructed before.
2925  * We can always find an appropriate variable to pivot with because
2926  * the current basis is mapped to the old basis by a non-singular
2927  * matrix and so we can never end up with a zero row.
2928  */
2929 static int restore_basis(struct isl_tab *tab, int *col_var)
2930 {
2931         int i, j;
2932         int n_extra = 0;
2933         int *extra = NULL;      /* current columns that contain bad stuff */
2934         unsigned off = 2 + tab->M;
2935
2936         extra = isl_alloc_array(tab->mat->ctx, int, tab->n_col);
2937         if (!extra)
2938                 goto error;
2939         for (i = 0; i < tab->n_col; ++i) {
2940                 for (j = 0; j < tab->n_col; ++j)
2941                         if (tab->col_var[i] == col_var[j])
2942                                 break;
2943                 if (j < tab->n_col)
2944                         continue;
2945                 extra[n_extra++] = i;
2946         }
2947         for (i = 0; i < tab->n_col && n_extra > 0; ++i) {
2948                 struct isl_tab_var *var;
2949                 int row;
2950
2951                 for (j = 0; j < tab->n_col; ++j)
2952                         if (col_var[i] == tab->col_var[j])
2953                                 break;
2954                 if (j < tab->n_col)
2955                         continue;
2956                 var = var_from_index(tab, col_var[i]);
2957                 row = var->index;
2958                 for (j = 0; j < n_extra; ++j)
2959                         if (!isl_int_is_zero(tab->mat->row[row][off+extra[j]]))
2960                                 break;
2961                 isl_assert(tab->mat->ctx, j < n_extra, goto error);
2962                 if (isl_tab_pivot(tab, row, extra[j]) < 0)
2963                         goto error;
2964                 extra[j] = extra[--n_extra];
2965         }
2966
2967         free(extra);
2968         free(col_var);
2969         return 0;
2970 error:
2971         free(extra);
2972         free(col_var);
2973         return -1;
2974 }
2975
2976 /* Remove all samples with index n or greater, i.e., those samples
2977  * that were added since we saved this number of samples in
2978  * isl_tab_save_samples.
2979  */
2980 static void drop_samples_since(struct isl_tab *tab, int n)
2981 {
2982         int i;
2983
2984         for (i = tab->n_sample - 1; i >= 0 && tab->n_sample > n; --i) {
2985                 if (tab->sample_index[i] < n)
2986                         continue;
2987
2988                 if (i != tab->n_sample - 1) {
2989                         int t = tab->sample_index[tab->n_sample-1];
2990                         tab->sample_index[tab->n_sample-1] = tab->sample_index[i];
2991                         tab->sample_index[i] = t;
2992                         isl_mat_swap_rows(tab->samples, tab->n_sample-1, i);
2993                 }
2994                 tab->n_sample--;
2995         }
2996 }
2997
2998 static int perform_undo(struct isl_tab *tab, struct isl_tab_undo *undo) WARN_UNUSED;
2999 static int perform_undo(struct isl_tab *tab, struct isl_tab_undo *undo)
3000 {
3001         switch (undo->type) {
3002         case isl_tab_undo_empty:
3003                 tab->empty = 0;
3004                 break;
3005         case isl_tab_undo_nonneg:
3006         case isl_tab_undo_redundant:
3007         case isl_tab_undo_freeze:
3008         case isl_tab_undo_zero:
3009         case isl_tab_undo_allocate:
3010         case isl_tab_undo_relax:
3011                 return perform_undo_var(tab, undo);
3012         case isl_tab_undo_bmap_eq:
3013                 return isl_basic_map_free_equality(tab->bmap, 1);
3014         case isl_tab_undo_bmap_ineq:
3015                 return isl_basic_map_free_inequality(tab->bmap, 1);
3016         case isl_tab_undo_bmap_div:
3017                 if (isl_basic_map_free_div(tab->bmap, 1) < 0)
3018                         return -1;
3019                 if (tab->samples)
3020                         tab->samples->n_col--;
3021                 break;
3022         case isl_tab_undo_saved_basis:
3023                 if (restore_basis(tab, undo->u.col_var) < 0)
3024                         return -1;
3025                 break;
3026         case isl_tab_undo_drop_sample:
3027                 tab->n_outside--;
3028                 break;
3029         case isl_tab_undo_saved_samples:
3030                 drop_samples_since(tab, undo->u.n);
3031                 break;
3032         case isl_tab_undo_callback:
3033                 return undo->u.callback->run(undo->u.callback);
3034         default:
3035                 isl_assert(tab->mat->ctx, 0, return -1);
3036         }
3037         return 0;
3038 }
3039
3040 /* Return the tableau to the state it was in when the snapshot "snap"
3041  * was taken.
3042  */
3043 int isl_tab_rollback(struct isl_tab *tab, struct isl_tab_undo *snap)
3044 {
3045         struct isl_tab_undo *undo, *next;
3046
3047         if (!tab)
3048                 return -1;
3049
3050         tab->in_undo = 1;
3051         for (undo = tab->top; undo && undo != &tab->bottom; undo = next) {
3052                 next = undo->next;
3053                 if (undo == snap)
3054                         break;
3055                 if (perform_undo(tab, undo) < 0) {
3056                         tab->top = undo;
3057                         free_undo(tab);
3058                         tab->in_undo = 0;
3059                         return -1;
3060                 }
3061                 free(undo);
3062         }
3063         tab->in_undo = 0;
3064         tab->top = undo;
3065         if (!undo)
3066                 return -1;
3067         return 0;
3068 }
3069
3070 /* The given row "row" represents an inequality violated by all
3071  * points in the tableau.  Check for some special cases of such
3072  * separating constraints.
3073  * In particular, if the row has been reduced to the constant -1,
3074  * then we know the inequality is adjacent (but opposite) to
3075  * an equality in the tableau.
3076  * If the row has been reduced to r = -1 -r', with r' an inequality
3077  * of the tableau, then the inequality is adjacent (but opposite)
3078  * to the inequality r'.
3079  */
3080 static enum isl_ineq_type separation_type(struct isl_tab *tab, unsigned row)
3081 {
3082         int pos;
3083         unsigned off = 2 + tab->M;
3084
3085         if (tab->rational)
3086                 return isl_ineq_separate;
3087
3088         if (!isl_int_is_one(tab->mat->row[row][0]))
3089                 return isl_ineq_separate;
3090         if (!isl_int_is_negone(tab->mat->row[row][1]))
3091                 return isl_ineq_separate;
3092
3093         pos = isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
3094                                         tab->n_col - tab->n_dead);
3095         if (pos == -1)
3096                 return isl_ineq_adj_eq;
3097
3098         if (!isl_int_is_negone(tab->mat->row[row][off + tab->n_dead + pos]))
3099                 return isl_ineq_separate;
3100
3101         pos = isl_seq_first_non_zero(
3102                         tab->mat->row[row] + off + tab->n_dead + pos + 1,
3103                         tab->n_col - tab->n_dead - pos - 1);
3104
3105         return pos == -1 ? isl_ineq_adj_ineq : isl_ineq_separate;
3106 }
3107
3108 /* Check the effect of inequality "ineq" on the tableau "tab".
3109  * The result may be
3110  *      isl_ineq_redundant:     satisfied by all points in the tableau
3111  *      isl_ineq_separate:      satisfied by no point in the tableau
3112  *      isl_ineq_cut:           satisfied by some by not all points
3113  *      isl_ineq_adj_eq:        adjacent to an equality
3114  *      isl_ineq_adj_ineq:      adjacent to an inequality.
3115  */
3116 enum isl_ineq_type isl_tab_ineq_type(struct isl_tab *tab, isl_int *ineq)
3117 {
3118         enum isl_ineq_type type = isl_ineq_error;
3119         struct isl_tab_undo *snap = NULL;
3120         int con;
3121         int row;
3122
3123         if (!tab)
3124                 return isl_ineq_error;
3125
3126         if (isl_tab_extend_cons(tab, 1) < 0)
3127                 return isl_ineq_error;
3128
3129         snap = isl_tab_snap(tab);
3130
3131         con = isl_tab_add_row(tab, ineq);
3132         if (con < 0)
3133                 goto error;
3134
3135         row = tab->con[con].index;
3136         if (isl_tab_row_is_redundant(tab, row))
3137                 type = isl_ineq_redundant;
3138         else if (isl_int_is_neg(tab->mat->row[row][1]) &&
3139                  (tab->rational ||
3140                     isl_int_abs_ge(tab->mat->row[row][1],
3141                                    tab->mat->row[row][0]))) {
3142                 int nonneg = at_least_zero(tab, &tab->con[con]);
3143                 if (nonneg < 0)
3144                         goto error;
3145                 if (nonneg)
3146                         type = isl_ineq_cut;
3147                 else
3148                         type = separation_type(tab, row);
3149         } else {
3150                 int red = con_is_redundant(tab, &tab->con[con]);
3151                 if (red < 0)
3152                         goto error;
3153                 if (!red)
3154                         type = isl_ineq_cut;
3155                 else
3156                         type = isl_ineq_redundant;
3157         }
3158
3159         if (isl_tab_rollback(tab, snap))
3160                 return isl_ineq_error;
3161         return type;
3162 error:
3163         return isl_ineq_error;
3164 }
3165
3166 int isl_tab_track_bmap(struct isl_tab *tab, __isl_take isl_basic_map *bmap)
3167 {
3168         if (!tab || !bmap)
3169                 goto error;
3170
3171         isl_assert(tab->mat->ctx, tab->n_eq == bmap->n_eq, return -1);
3172         isl_assert(tab->mat->ctx,
3173                     tab->n_con == bmap->n_eq + bmap->n_ineq, return -1);
3174
3175         tab->bmap = bmap;
3176
3177         return 0;
3178 error:
3179         isl_basic_map_free(bmap);
3180         return -1;
3181 }
3182
3183 int isl_tab_track_bset(struct isl_tab *tab, __isl_take isl_basic_set *bset)
3184 {
3185         return isl_tab_track_bmap(tab, (isl_basic_map *)bset);
3186 }
3187
3188 __isl_keep isl_basic_set *isl_tab_peek_bset(struct isl_tab *tab)
3189 {
3190         if (!tab)
3191                 return NULL;
3192
3193         return (isl_basic_set *)tab->bmap;
3194 }
3195
3196 void isl_tab_dump(struct isl_tab *tab, FILE *out, int indent)
3197 {
3198         unsigned r, c;
3199         int i;
3200
3201         if (!tab) {
3202                 fprintf(out, "%*snull tab\n", indent, "");
3203                 return;
3204         }
3205         fprintf(out, "%*sn_redundant: %d, n_dead: %d", indent, "",
3206                 tab->n_redundant, tab->n_dead);
3207         if (tab->rational)
3208                 fprintf(out, ", rational");
3209         if (tab->empty)
3210                 fprintf(out, ", empty");
3211         fprintf(out, "\n");
3212         fprintf(out, "%*s[", indent, "");
3213         for (i = 0; i < tab->n_var; ++i) {
3214                 if (i)
3215                         fprintf(out, (i == tab->n_param ||
3216                                       i == tab->n_var - tab->n_div) ? "; "
3217                                                                     : ", ");
3218                 fprintf(out, "%c%d%s", tab->var[i].is_row ? 'r' : 'c',
3219                                         tab->var[i].index,
3220                                         tab->var[i].is_zero ? " [=0]" :
3221                                         tab->var[i].is_redundant ? " [R]" : "");
3222         }
3223         fprintf(out, "]\n");
3224         fprintf(out, "%*s[", indent, "");
3225         for (i = 0; i < tab->n_con; ++i) {
3226                 if (i)
3227                         fprintf(out, ", ");
3228                 fprintf(out, "%c%d%s", tab->con[i].is_row ? 'r' : 'c',
3229                                         tab->con[i].index,
3230                                         tab->con[i].is_zero ? " [=0]" :
3231                                         tab->con[i].is_redundant ? " [R]" : "");
3232         }
3233         fprintf(out, "]\n");
3234         fprintf(out, "%*s[", indent, "");
3235         for (i = 0; i < tab->n_row; ++i) {
3236                 const char *sign = "";
3237                 if (i)
3238                         fprintf(out, ", ");
3239                 if (tab->row_sign) {
3240                         if (tab->row_sign[i] == isl_tab_row_unknown)
3241                                 sign = "?";
3242                         else if (tab->row_sign[i] == isl_tab_row_neg)
3243                                 sign = "-";
3244                         else if (tab->row_sign[i] == isl_tab_row_pos)
3245                                 sign = "+";
3246                         else
3247                                 sign = "+-";
3248                 }
3249                 fprintf(out, "r%d: %d%s%s", i, tab->row_var[i],
3250                     isl_tab_var_from_row(tab, i)->is_nonneg ? " [>=0]" : "", sign);
3251         }
3252         fprintf(out, "]\n");
3253         fprintf(out, "%*s[", indent, "");
3254         for (i = 0; i < tab->n_col; ++i) {
3255                 if (i)
3256                         fprintf(out, ", ");
3257                 fprintf(out, "c%d: %d%s", i, tab->col_var[i],
3258                     var_from_col(tab, i)->is_nonneg ? " [>=0]" : "");
3259         }
3260         fprintf(out, "]\n");
3261         r = tab->mat->n_row;
3262         tab->mat->n_row = tab->n_row;
3263         c = tab->mat->n_col;
3264         tab->mat->n_col = 2 + tab->M + tab->n_col;
3265         isl_mat_dump(tab->mat, out, indent);
3266         tab->mat->n_row = r;
3267         tab->mat->n_col = c;
3268         if (tab->bmap)
3269                 isl_basic_map_print_internal(tab->bmap, out, indent);
3270 }