Coverage Report

Created: 2019-07-24 05:18

/Users/buildslave/jenkins/workspace/clang-stage2-coverage-R/llvm/tools/polly/lib/External/isl/isl_tab.c
Line
Count
Source (jump to first uncovered line)
1
/*
2
 * Copyright 2008-2009 Katholieke Universiteit Leuven
3
 * Copyright 2013      Ecole Normale Superieure
4
 * Copyright 2014      INRIA Rocquencourt
5
 * Copyright 2016      Sven Verdoolaege
6
 *
7
 * Use of this software is governed by the MIT license
8
 *
9
 * Written by Sven Verdoolaege, K.U.Leuven, Departement
10
 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
11
 * and Ecole Normale Superieure, 45 rue d'Ulm, 75230 Paris, France
12
 * and Inria Paris - Rocquencourt, Domaine de Voluceau - Rocquencourt,
13
 * B.P. 105 - 78153 Le Chesnay, France
14
 */
15
16
#include <isl_ctx_private.h>
17
#include <isl_mat_private.h>
18
#include <isl_vec_private.h>
19
#include "isl_map_private.h"
20
#include "isl_tab.h"
21
#include <isl_seq.h>
22
#include <isl_config.h>
23
24
#include <bset_to_bmap.c>
25
#include <bset_from_bmap.c>
26
27
/*
28
 * The implementation of tableaus in this file was inspired by Section 8
29
 * of David Detlefs, Greg Nelson and James B. Saxe, "Simplify: a theorem
30
 * prover for program checking".
31
 */
32
33
struct isl_tab *isl_tab_alloc(struct isl_ctx *ctx,
34
  unsigned n_row, unsigned n_var, unsigned M)
35
1.22M
{
36
1.22M
  int i;
37
1.22M
  struct isl_tab *tab;
38
1.22M
  unsigned off = 2 + M;
39
1.22M
40
1.22M
  tab = isl_calloc_type(ctx, struct isl_tab);
41
1.22M
  if (!tab)
42
0
    return NULL;
43
1.22M
  tab->mat = isl_mat_alloc(ctx, n_row, off + n_var);
44
1.22M
  if (!tab->mat)
45
0
    goto error;
46
1.22M
  tab->var = isl_alloc_array(ctx, struct isl_tab_var, n_var);
47
1.22M
  if (n_var && 
!tab->var1.22M
)
48
0
    goto error;
49
1.22M
  tab->con = isl_alloc_array(ctx, struct isl_tab_var, n_row);
50
1.22M
  if (n_row && 
!tab->con1.22M
)
51
0
    goto error;
52
1.22M
  tab->col_var = isl_alloc_array(ctx, int, n_var);
53
1.22M
  if (n_var && 
!tab->col_var1.22M
)
54
0
    goto error;
55
1.22M
  tab->row_var = isl_alloc_array(ctx, int, n_row);
56
1.22M
  if (n_row && 
!tab->row_var1.22M
)
57
0
    goto error;
58
7.66M
  
for (i = 0; 1.22M
i < n_var;
++i6.43M
) {
59
6.43M
    tab->var[i].index = i;
60
6.43M
    tab->var[i].is_row = 0;
61
6.43M
    tab->var[i].is_nonneg = 0;
62
6.43M
    tab->var[i].is_zero = 0;
63
6.43M
    tab->var[i].is_redundant = 0;
64
6.43M
    tab->var[i].frozen = 0;
65
6.43M
    tab->var[i].negated = 0;
66
6.43M
    tab->col_var[i] = i;
67
6.43M
  }
68
1.22M
  tab->n_row = 0;
69
1.22M
  tab->n_con = 0;
70
1.22M
  tab->n_eq = 0;
71
1.22M
  tab->max_con = n_row;
72
1.22M
  tab->n_col = n_var;
73
1.22M
  tab->n_var = n_var;
74
1.22M
  tab->max_var = n_var;
75
1.22M
  tab->n_param = 0;
76
1.22M
  tab->n_div = 0;
77
1.22M
  tab->n_dead = 0;
78
1.22M
  tab->n_redundant = 0;
79
1.22M
  tab->strict_redundant = 0;
80
1.22M
  tab->need_undo = 0;
81
1.22M
  tab->rational = 0;
82
1.22M
  tab->empty = 0;
83
1.22M
  tab->in_undo = 0;
84
1.22M
  tab->M = M;
85
1.22M
  tab->cone = 0;
86
1.22M
  tab->bottom.type = isl_tab_undo_bottom;
87
1.22M
  tab->bottom.next = NULL;
88
1.22M
  tab->top = &tab->bottom;
89
1.22M
90
1.22M
  tab->n_zero = 0;
91
1.22M
  tab->n_unbounded = 0;
92
1.22M
  tab->basis = NULL;
93
1.22M
94
1.22M
  return tab;
95
0
error:
96
0
  isl_tab_free(tab);
97
0
  return NULL;
98
1.22M
}
99
100
isl_ctx *isl_tab_get_ctx(struct isl_tab *tab)
101
6.74M
{
102
6.74M
  return tab ? isl_mat_get_ctx(tab->mat) : NULL;
103
6.74M
}
104
105
int isl_tab_extend_cons(struct isl_tab *tab, unsigned n_new)
106
1.45M
{
107
1.45M
  unsigned off;
108
1.45M
109
1.45M
  if (!tab)
110
0
    return -1;
111
1.45M
112
1.45M
  off = 2 + tab->M;
113
1.45M
114
1.45M
  if (tab->max_con < tab->n_con + n_new) {
115
192k
    struct isl_tab_var *con;
116
192k
117
192k
    con = isl_realloc_array(tab->mat->ctx, tab->con,
118
192k
            struct isl_tab_var, tab->max_con + n_new);
119
192k
    if (!con)
120
0
      return -1;
121
192k
    tab->con = con;
122
192k
    tab->max_con += n_new;
123
192k
  }
124
1.45M
  if (tab->mat->n_row < tab->n_row + n_new) {
125
195k
    int *row_var;
126
195k
127
195k
    tab->mat = isl_mat_extend(tab->mat,
128
195k
          tab->n_row + n_new, off + tab->n_col);
129
195k
    if (!tab->mat)
130
0
      return -1;
131
195k
    row_var = isl_realloc_array(tab->mat->ctx, tab->row_var,
132
195k
              int, tab->mat->n_row);
133
195k
    if (!row_var)
134
0
      return -1;
135
195k
    tab->row_var = row_var;
136
195k
    if (tab->row_sign) {
137
408
      enum isl_tab_row_sign *s;
138
408
      s = isl_realloc_array(tab->mat->ctx, tab->row_sign,
139
408
          enum isl_tab_row_sign, tab->mat->n_row);
140
408
      if (!s)
141
0
        return -1;
142
408
      tab->row_sign = s;
143
408
    }
144
195k
  }
145
1.45M
  return 0;
146
1.45M
}
147
148
/* Make room for at least n_new extra variables.
149
 * Return -1 if anything went wrong.
150
 */
151
int isl_tab_extend_vars(struct isl_tab *tab, unsigned n_new)
152
8.31k
{
153
8.31k
  struct isl_tab_var *var;
154
8.31k
  unsigned off = 2 + tab->M;
155
8.31k
156
8.31k
  if (tab->max_var < tab->n_var + n_new) {
157
5.63k
    var = isl_realloc_array(tab->mat->ctx, tab->var,
158
5.63k
            struct isl_tab_var, tab->n_var + n_new);
159
5.63k
    if (!var)
160
0
      return -1;
161
5.63k
    tab->var = var;
162
5.63k
    tab->max_var = tab->n_var + n_new;
163
5.63k
  }
164
8.31k
165
8.31k
  if (tab->mat->n_col < off + tab->n_col + n_new) {
166
3.75k
    int *p;
167
3.75k
168
3.75k
    tab->mat = isl_mat_extend(tab->mat,
169
3.75k
            tab->mat->n_row, off + tab->n_col + n_new);
170
3.75k
    if (!tab->mat)
171
0
      return -1;
172
3.75k
    p = isl_realloc_array(tab->mat->ctx, tab->col_var,
173
3.75k
              int, tab->n_col + n_new);
174
3.75k
    if (!p)
175
0
      return -1;
176
3.75k
    tab->col_var = p;
177
3.75k
  }
178
8.31k
179
8.31k
  return 0;
180
8.31k
}
181
182
static void free_undo_record(struct isl_tab_undo *undo)
183
5.99M
{
184
5.99M
  switch (undo->type) {
185
5.99M
  case isl_tab_undo_saved_basis:
186
766
    free(undo->u.col_var);
187
766
    break;
188
5.99M
  
default:;5.99M
189
5.99M
  }
190
5.99M
  free(undo);
191
5.99M
}
192
193
static void free_undo(struct isl_tab *tab)
194
1.23M
{
195
1.23M
  struct isl_tab_undo *undo, *next;
196
1.23M
197
2.10M
  for (undo = tab->top; undo && undo != &tab->bottom; 
undo = next867k
) {
198
867k
    next = undo->next;
199
867k
    free_undo_record(undo);
200
867k
  }
201
1.23M
  tab->top = undo;
202
1.23M
}
203
204
void isl_tab_free(struct isl_tab *tab)
205
1.28M
{
206
1.28M
  if (!tab)
207
51.1k
    return;
208
1.23M
  free_undo(tab);
209
1.23M
  isl_mat_free(tab->mat);
210
1.23M
  isl_vec_free(tab->dual);
211
1.23M
  isl_basic_map_free(tab->bmap);
212
1.23M
  free(tab->var);
213
1.23M
  free(tab->con);
214
1.23M
  free(tab->row_var);
215
1.23M
  free(tab->col_var);
216
1.23M
  free(tab->row_sign);
217
1.23M
  isl_mat_free(tab->samples);
218
1.23M
  free(tab->sample_index);
219
1.23M
  isl_mat_free(tab->basis);
220
1.23M
  free(tab);
221
1.23M
}
222
223
struct isl_tab *isl_tab_dup(struct isl_tab *tab)
224
2.85k
{
225
2.85k
  int i;
226
2.85k
  struct isl_tab *dup;
227
2.85k
  unsigned off;
228
2.85k
229
2.85k
  if (!tab)
230
0
    return NULL;
231
2.85k
232
2.85k
  off = 2 + tab->M;
233
2.85k
  dup = isl_calloc_type(tab->mat->ctx, struct isl_tab);
234
2.85k
  if (!dup)
235
0
    return NULL;
236
2.85k
  dup->mat = isl_mat_dup(tab->mat);
237
2.85k
  if (!dup->mat)
238
0
    goto error;
239
2.85k
  dup->var = isl_alloc_array(tab->mat->ctx, struct isl_tab_var, tab->max_var);
240
2.85k
  if (tab->max_var && !dup->var)
241
0
    goto error;
242
27.2k
  
for (i = 0; 2.85k
i < tab->n_var;
++i24.4k
)
243
24.4k
    dup->var[i] = tab->var[i];
244
2.85k
  dup->con = isl_alloc_array(tab->mat->ctx, struct isl_tab_var, tab->max_con);
245
2.85k
  if (tab->max_con && !dup->con)
246
0
    goto error;
247
30.4k
  
for (i = 0; 2.85k
i < tab->n_con;
++i27.5k
)
248
27.5k
    dup->con[i] = tab->con[i];
249
2.85k
  dup->col_var = isl_alloc_array(tab->mat->ctx, int, tab->mat->n_col - off);
250
2.85k
  if ((tab->mat->n_col - off) && !dup->col_var)
251
0
    goto error;
252
14.5k
  
for (i = 0; 2.85k
i < tab->n_col;
++i11.7k
)
253
11.7k
    dup->col_var[i] = tab->col_var[i];
254
2.85k
  dup->row_var = isl_alloc_array(tab->mat->ctx, int, tab->mat->n_row);
255
2.85k
  if (tab->mat->n_row && !dup->row_var)
256
0
    goto error;
257
30.2k
  
for (i = 0; 2.85k
i < tab->n_row;
++i27.4k
)
258
27.4k
    dup->row_var[i] = tab->row_var[i];
259
2.85k
  if (tab->row_sign) {
260
2.84k
    dup->row_sign = isl_alloc_array(tab->mat->ctx, enum isl_tab_row_sign,
261
2.84k
            tab->mat->n_row);
262
2.84k
    if (tab->mat->n_row && !dup->row_sign)
263
0
      goto error;
264
30.2k
    
for (i = 0; 2.84k
i < tab->n_row;
++i27.3k
)
265
27.3k
      dup->row_sign[i] = tab->row_sign[i];
266
2.84k
  }
267
2.85k
  if (tab->samples) {
268
0
    dup->samples = isl_mat_dup(tab->samples);
269
0
    if (!dup->samples)
270
0
      goto error;
271
0
    dup->sample_index = isl_alloc_array(tab->mat->ctx, int,
272
0
              tab->samples->n_row);
273
0
    if (tab->samples->n_row && !dup->sample_index)
274
0
      goto error;
275
0
    dup->n_sample = tab->n_sample;
276
0
    dup->n_outside = tab->n_outside;
277
0
  }
278
2.85k
  dup->n_row = tab->n_row;
279
2.85k
  dup->n_con = tab->n_con;
280
2.85k
  dup->n_eq = tab->n_eq;
281
2.85k
  dup->max_con = tab->max_con;
282
2.85k
  dup->n_col = tab->n_col;
283
2.85k
  dup->n_var = tab->n_var;
284
2.85k
  dup->max_var = tab->max_var;
285
2.85k
  dup->n_param = tab->n_param;
286
2.85k
  dup->n_div = tab->n_div;
287
2.85k
  dup->n_dead = tab->n_dead;
288
2.85k
  dup->n_redundant = tab->n_redundant;
289
2.85k
  dup->rational = tab->rational;
290
2.85k
  dup->empty = tab->empty;
291
2.85k
  dup->strict_redundant = 0;
292
2.85k
  dup->need_undo = 0;
293
2.85k
  dup->in_undo = 0;
294
2.85k
  dup->M = tab->M;
295
2.85k
  tab->cone = tab->cone;
296
2.85k
  dup->bottom.type = isl_tab_undo_bottom;
297
2.85k
  dup->bottom.next = NULL;
298
2.85k
  dup->top = &dup->bottom;
299
2.85k
300
2.85k
  dup->n_zero = tab->n_zero;
301
2.85k
  dup->n_unbounded = tab->n_unbounded;
302
2.85k
  dup->basis = isl_mat_dup(tab->basis);
303
2.85k
304
2.85k
  return dup;
305
0
error:
306
0
  isl_tab_free(dup);
307
0
  return NULL;
308
2.85k
}
309
310
/* Construct the coefficient matrix of the product tableau
311
 * of two tableaus.
312
 * mat{1,2} is the coefficient matrix of tableau {1,2}
313
 * row{1,2} is the number of rows in tableau {1,2}
314
 * col{1,2} is the number of columns in tableau {1,2}
315
 * off is the offset to the coefficient column (skipping the
316
 *  denominator, the constant term and the big parameter if any)
317
 * r{1,2} is the number of redundant rows in tableau {1,2}
318
 * d{1,2} is the number of dead columns in tableau {1,2}
319
 *
320
 * The order of the rows and columns in the result is as explained
321
 * in isl_tab_product.
322
 */
323
static struct isl_mat *tab_mat_product(struct isl_mat *mat1,
324
  struct isl_mat *mat2, unsigned row1, unsigned row2,
325
  unsigned col1, unsigned col2,
326
  unsigned off, unsigned r1, unsigned r2, unsigned d1, unsigned d2)
327
6.84k
{
328
6.84k
  int i;
329
6.84k
  struct isl_mat *prod;
330
6.84k
  unsigned n;
331
6.84k
332
6.84k
  prod = isl_mat_alloc(mat1->ctx, mat1->n_row + mat2->n_row,
333
6.84k
          off + col1 + col2);
334
6.84k
  if (!prod)
335
0
    return NULL;
336
6.84k
337
6.84k
  n = 0;
338
27.3k
  for (i = 0; i < r1; 
++i20.4k
) {
339
20.4k
    isl_seq_cpy(prod->row[n + i], mat1->row[i], off + d1);
340
20.4k
    isl_seq_clr(prod->row[n + i] + off + d1, d2);
341
20.4k
    isl_seq_cpy(prod->row[n + i] + off + d1 + d2,
342
20.4k
        mat1->row[i] + off + d1, col1 - d1);
343
20.4k
    isl_seq_clr(prod->row[n + i] + off + col1 + d1, col2 - d2);
344
20.4k
  }
345
6.84k
346
6.84k
  n += r1;
347
27.3k
  for (i = 0; i < r2; 
++i20.4k
) {
348
20.4k
    isl_seq_cpy(prod->row[n + i], mat2->row[i], off);
349
20.4k
    isl_seq_clr(prod->row[n + i] + off, d1);
350
20.4k
    isl_seq_cpy(prod->row[n + i] + off + d1,
351
20.4k
          mat2->row[i] + off, d2);
352
20.4k
    isl_seq_clr(prod->row[n + i] + off + d1 + d2, col1 - d1);
353
20.4k
    isl_seq_cpy(prod->row[n + i] + off + col1 + d1,
354
20.4k
          mat2->row[i] + off + d2, col2 - d2);
355
20.4k
  }
356
6.84k
357
6.84k
  n += r2;
358
77.9k
  for (i = 0; i < row1 - r1; 
++i71.0k
) {
359
71.0k
    isl_seq_cpy(prod->row[n + i], mat1->row[r1 + i], off + d1);
360
71.0k
    isl_seq_clr(prod->row[n + i] + off + d1, d2);
361
71.0k
    isl_seq_cpy(prod->row[n + i] + off + d1 + d2,
362
71.0k
        mat1->row[r1 + i] + off + d1, col1 - d1);
363
71.0k
    isl_seq_clr(prod->row[n + i] + off + col1 + d1, col2 - d2);
364
71.0k
  }
365
6.84k
366
6.84k
  n += row1 - r1;
367
77.9k
  for (i = 0; i < row2 - r2; 
++i71.0k
) {
368
71.0k
    isl_seq_cpy(prod->row[n + i], mat2->row[r2 + i], off);
369
71.0k
    isl_seq_clr(prod->row[n + i] + off, d1);
370
71.0k
    isl_seq_cpy(prod->row[n + i] + off + d1,
371
71.0k
          mat2->row[r2 + i] + off, d2);
372
71.0k
    isl_seq_clr(prod->row[n + i] + off + d1 + d2, col1 - d1);
373
71.0k
    isl_seq_cpy(prod->row[n + i] + off + col1 + d1,
374
71.0k
          mat2->row[r2 + i] + off + d2, col2 - d2);
375
71.0k
  }
376
6.84k
377
6.84k
  return prod;
378
6.84k
}
379
380
/* Update the row or column index of a variable that corresponds
381
 * to a variable in the first input tableau.
382
 */
383
static void update_index1(struct isl_tab_var *var,
384
  unsigned r1, unsigned r2, unsigned d1, unsigned d2)
385
133k
{
386
133k
  if (var->index == -1)
387
327
    return;
388
133k
  if (var->is_row && 
var->index >= r191.5k
)
389
71.0k
    var->index += r2;
390
133k
  if (!var->is_row && 
var->index >= d141.4k
)
391
39.6k
    var->index += d2;
392
133k
}
393
394
/* Update the row or column index of a variable that corresponds
395
 * to a variable in the second input tableau.
396
 */
397
static void update_index2(struct isl_tab_var *var,
398
  unsigned row1, unsigned col1,
399
  unsigned r1, unsigned r2, unsigned d1, unsigned d2)
400
133k
{
401
133k
  if (var->index == -1)
402
327
    return;
403
133k
  if (var->is_row) {
404
91.5k
    if (var->index < r2)
405
20.4k
      var->index += r1;
406
71.0k
    else
407
71.0k
      var->index += row1;
408
91.5k
  } else {
409
41.4k
    if (var->index < d2)
410
1.79k
      var->index += d1;
411
39.6k
    else
412
39.6k
      var->index += col1;
413
41.4k
  }
414
133k
}
415
416
/* Create a tableau that represents the Cartesian product of the sets
417
 * represented by tableaus tab1 and tab2.
418
 * The order of the rows in the product is
419
 *  - redundant rows of tab1
420
 *  - redundant rows of tab2
421
 *  - non-redundant rows of tab1
422
 *  - non-redundant rows of tab2
423
 * The order of the columns is
424
 *  - denominator
425
 *  - constant term
426
 *  - coefficient of big parameter, if any
427
 *  - dead columns of tab1
428
 *  - dead columns of tab2
429
 *  - live columns of tab1
430
 *  - live columns of tab2
431
 * The order of the variables and the constraints is a concatenation
432
 * of order in the two input tableaus.
433
 */
434
struct isl_tab *isl_tab_product(struct isl_tab *tab1, struct isl_tab *tab2)
435
6.84k
{
436
6.84k
  int i;
437
6.84k
  struct isl_tab *prod;
438
6.84k
  unsigned off;
439
6.84k
  unsigned r1, r2, d1, d2;
440
6.84k
441
6.84k
  if (!tab1 || !tab2)
442
0
    return NULL;
443
6.84k
444
6.84k
  isl_assert(tab1->mat->ctx, tab1->M == tab2->M, return NULL);
445
6.84k
  isl_assert(tab1->mat->ctx, tab1->rational == tab2->rational, return NULL);
446
6.84k
  isl_assert(tab1->mat->ctx, tab1->cone == tab2->cone, return NULL);
447
6.84k
  isl_assert(tab1->mat->ctx, !tab1->row_sign, return NULL);
448
6.84k
  isl_assert(tab1->mat->ctx, !tab2->row_sign, return NULL);
449
6.84k
  isl_assert(tab1->mat->ctx, tab1->n_param == 0, return NULL);
450
6.84k
  isl_assert(tab1->mat->ctx, tab2->n_param == 0, return NULL);
451
6.84k
  isl_assert(tab1->mat->ctx, tab1->n_div == 0, return NULL);
452
6.84k
  isl_assert(tab1->mat->ctx, tab2->n_div == 0, return NULL);
453
6.84k
454
6.84k
  off = 2 + tab1->M;
455
6.84k
  r1 = tab1->n_redundant;
456
6.84k
  r2 = tab2->n_redundant;
457
6.84k
  d1 = tab1->n_dead;
458
6.84k
  d2 = tab2->n_dead;
459
6.84k
  prod = isl_calloc_type(tab1->mat->ctx, struct isl_tab);
460
6.84k
  if (!prod)
461
0
    return NULL;
462
6.84k
  prod->mat = tab_mat_product(tab1->mat, tab2->mat,
463
6.84k
        tab1->n_row, tab2->n_row,
464
6.84k
        tab1->n_col, tab2->n_col, off, r1, r2, d1, d2);
465
6.84k
  if (!prod->mat)
466
0
    goto error;
467
6.84k
  prod->var = isl_alloc_array(tab1->mat->ctx, struct isl_tab_var,
468
6.84k
          tab1->max_var + tab2->max_var);
469
6.84k
  if ((tab1->max_var + tab2->max_var) && !prod->var)
470
0
    goto error;
471
48.6k
  
for (i = 0; 6.84k
i < tab1->n_var;
++i41.7k
) {
472
41.7k
    prod->var[i] = tab1->var[i];
473
41.7k
    update_index1(&prod->var[i], r1, r2, d1, d2);
474
41.7k
  }
475
48.6k
  for (i = 0; i < tab2->n_var; 
++i41.7k
) {
476
41.7k
    prod->var[tab1->n_var + i] = tab2->var[i];
477
41.7k
    update_index2(&prod->var[tab1->n_var + i],
478
41.7k
        tab1->n_row, tab1->n_col,
479
41.7k
        r1, r2, d1, d2);
480
41.7k
  }
481
6.84k
  prod->con = isl_alloc_array(tab1->mat->ctx, struct isl_tab_var,
482
6.84k
          tab1->max_con +  tab2->max_con);
483
6.84k
  if ((tab1->max_con + tab2->max_con) && !prod->con)
484
0
    goto error;
485
98.4k
  
for (i = 0; 6.84k
i < tab1->n_con;
++i91.5k
) {
486
91.5k
    prod->con[i] = tab1->con[i];
487
91.5k
    update_index1(&prod->con[i], r1, r2, d1, d2);
488
91.5k
  }
489
98.4k
  for (i = 0; i < tab2->n_con; 
++i91.5k
) {
490
91.5k
    prod->con[tab1->n_con + i] = tab2->con[i];
491
91.5k
    update_index2(&prod->con[tab1->n_con + i],
492
91.5k
        tab1->n_row, tab1->n_col,
493
91.5k
        r1, r2, d1, d2);
494
91.5k
  }
495
6.84k
  prod->col_var = isl_alloc_array(tab1->mat->ctx, int,
496
6.84k
          tab1->n_col + tab2->n_col);
497
6.84k
  if ((tab1->n_col + tab2->n_col) && !prod->col_var)
498
0
    goto error;
499
48.2k
  
for (i = 0; 6.84k
i < tab1->n_col;
++i41.4k
) {
500
41.4k
    int pos = i < d1 ? 
i1.79k
:
i + d239.6k
;
501
41.4k
    prod->col_var[pos] = tab1->col_var[i];
502
41.4k
  }
503
48.2k
  for (i = 0; i < tab2->n_col; 
++i41.4k
) {
504
41.4k
    int pos = i < d2 ? 
d1 + i1.79k
:
tab1->n_col + i39.6k
;
505
41.4k
    int t = tab2->col_var[i];
506
41.4k
    if (t >= 0)
507
597
      t += tab1->n_var;
508
40.8k
    else
509
40.8k
      t -= tab1->n_con;
510
41.4k
    prod->col_var[pos] = t;
511
41.4k
  }
512
6.84k
  prod->row_var = isl_alloc_array(tab1->mat->ctx, int,
513
6.84k
          tab1->mat->n_row + tab2->mat->n_row);
514
6.84k
  if ((tab1->mat->n_row + tab2->mat->n_row) && !prod->row_var)
515
0
    goto error;
516
98.4k
  
for (i = 0; 6.84k
i < tab1->n_row;
++i91.5k
) {
517
91.5k
    int pos = i < r1 ? 
i20.4k
:
i + r271.0k
;
518
91.5k
    prod->row_var[pos] = tab1->row_var[i];
519
91.5k
  }
520
98.4k
  for (i = 0; i < tab2->n_row; 
++i91.5k
) {
521
91.5k
    int pos = i < r2 ? 
r1 + i20.4k
:
tab1->n_row + i71.0k
;
522
91.5k
    int t = tab2->row_var[i];
523
91.5k
    if (t >= 0)
524
41.1k
      t += tab1->n_var;
525
50.3k
    else
526
50.3k
      t -= tab1->n_con;
527
91.5k
    prod->row_var[pos] = t;
528
91.5k
  }
529
6.84k
  prod->samples = NULL;
530
6.84k
  prod->sample_index = NULL;
531
6.84k
  prod->n_row = tab1->n_row + tab2->n_row;
532
6.84k
  prod->n_con = tab1->n_con + tab2->n_con;
533
6.84k
  prod->n_eq = 0;
534
6.84k
  prod->max_con = tab1->max_con + tab2->max_con;
535
6.84k
  prod->n_col = tab1->n_col + tab2->n_col;
536
6.84k
  prod->n_var = tab1->n_var + tab2->n_var;
537
6.84k
  prod->max_var = tab1->max_var + tab2->max_var;
538
6.84k
  prod->n_param = 0;
539
6.84k
  prod->n_div = 0;
540
6.84k
  prod->n_dead = tab1->n_dead + tab2->n_dead;
541
6.84k
  prod->n_redundant = tab1->n_redundant + tab2->n_redundant;
542
6.84k
  prod->rational = tab1->rational;
543
6.84k
  prod->empty = tab1->empty || tab2->empty;
544
6.84k
  prod->strict_redundant = tab1->strict_redundant || tab2->strict_redundant;
545
6.84k
  prod->need_undo = 0;
546
6.84k
  prod->in_undo = 0;
547
6.84k
  prod->M = tab1->M;
548
6.84k
  prod->cone = tab1->cone;
549
6.84k
  prod->bottom.type = isl_tab_undo_bottom;
550
6.84k
  prod->bottom.next = NULL;
551
6.84k
  prod->top = &prod->bottom;
552
6.84k
553
6.84k
  prod->n_zero = 0;
554
6.84k
  prod->n_unbounded = 0;
555
6.84k
  prod->basis = NULL;
556
6.84k
557
6.84k
  return prod;
558
0
error:
559
0
  isl_tab_free(prod);
560
0
  return NULL;
561
6.84k
}
562
563
static struct isl_tab_var *var_from_index(struct isl_tab *tab, int i)
564
210M
{
565
210M
  if (i >= 0)
566
58.7M
    return &tab->var[i];
567
151M
  else
568
151M
    return &tab->con[~i];
569
210M
}
570
571
struct isl_tab_var *isl_tab_var_from_row(struct isl_tab *tab, int i)
572
158M
{
573
158M
  return var_from_index(tab, tab->row_var[i]);
574
158M
}
575
576
static struct isl_tab_var *var_from_col(struct isl_tab *tab, int i)
577
46.9M
{
578
46.9M
  return var_from_index(tab, tab->col_var[i]);
579
46.9M
}
580
581
/* Check if there are any upper bounds on column variable "var",
582
 * i.e., non-negative rows where var appears with a negative coefficient.
583
 * Return 1 if there are no such bounds.
584
 */
585
static int max_is_manifestly_unbounded(struct isl_tab *tab,
586
  struct isl_tab_var *var)
587
3.94M
{
588
3.94M
  int i;
589
3.94M
  unsigned off = 2 + tab->M;
590
3.94M
591
3.94M
  if (var->is_row)
592
2.82M
    return 0;
593
7.75M
  
for (i = tab->n_redundant; 1.12M
i < tab->n_row;
++i6.63M
) {
594
7.13M
    if (!isl_int_is_neg(tab->mat->row[i][off + var->index]))
595
7.13M
      
continue6.25M
;
596
880k
    if (isl_tab_var_from_row(tab, i)->is_nonneg)
597
496k
      return 0;
598
880k
  }
599
1.12M
  
return 1623k
;
600
1.12M
}
601
602
/* Check if there are any lower bounds on column variable "var",
603
 * i.e., non-negative rows where var appears with a positive coefficient.
604
 * Return 1 if there are no such bounds.
605
 */
606
static int min_is_manifestly_unbounded(struct isl_tab *tab,
607
  struct isl_tab_var *var)
608
3.24M
{
609
3.24M
  int i;
610
3.24M
  unsigned off = 2 + tab->M;
611
3.24M
612
3.24M
  if (var->is_row)
613
1.51M
    return 0;
614
13.8M
  
for (i = tab->n_redundant; 1.72M
i < tab->n_row;
++i12.1M
) {
615
13.0M
    if (!isl_int_is_pos(tab->mat->row[i][off + var->index]))
616
13.0M
      
continue11.5M
;
617
1.51M
    if (isl_tab_var_from_row(tab, i)->is_nonneg)
618
956k
      return 0;
619
1.51M
  }
620
1.72M
  
return 1767k
;
621
1.72M
}
622
623
static int row_cmp(struct isl_tab *tab, int r1, int r2, int c, isl_int *t)
624
2.29M
{
625
2.29M
  unsigned off = 2 + tab->M;
626
2.29M
627
2.29M
  if (tab->M) {
628
0
    int s;
629
0
    isl_int_mul(*t, tab->mat->row[r1][2], tab->mat->row[r2][off+c]);
630
0
    isl_int_submul(*t, tab->mat->row[r2][2], tab->mat->row[r1][off+c]);
631
0
    s = isl_int_sgn(*t);
632
0
    if (s)
633
0
      return s;
634
2.29M
  }
635
2.29M
  isl_int_mul(*t, tab->mat->row[r1][1], tab->mat->row[r2][off + c]);
636
2.29M
  isl_int_submul(*t, tab->mat->row[r2][1], tab->mat->row[r1][off + c]);
637
2.29M
  return isl_int_sgn(*t);
638
2.29M
}
639
640
/* Given the index of a column "c", return the index of a row
641
 * that can be used to pivot the column in, with either an increase
642
 * (sgn > 0) or a decrease (sgn < 0) of the corresponding variable.
643
 * If "var" is not NULL, then the row returned will be different from
644
 * the one associated with "var".
645
 *
646
 * Each row in the tableau is of the form
647
 *
648
 *  x_r = a_r0 + \sum_i a_ri x_i
649
 *
650
 * Only rows with x_r >= 0 and with the sign of a_ri opposite to "sgn"
651
 * impose any limit on the increase or decrease in the value of x_c
652
 * and this bound is equal to a_r0 / |a_rc|.  We are therefore looking
653
 * for the row with the smallest (most stringent) such bound.
654
 * Note that the common denominator of each row drops out of the fraction.
655
 * To check if row j has a smaller bound than row r, i.e.,
656
 * a_j0 / |a_jc| < a_r0 / |a_rc| or a_j0 |a_rc| < a_r0 |a_jc|,
657
 * we check if -sign(a_jc) (a_j0 a_rc - a_r0 a_jc) < 0,
658
 * where -sign(a_jc) is equal to "sgn".
659
 */
660
static int pivot_row(struct isl_tab *tab,
661
  struct isl_tab_var *var, int sgn, int c)
662
6.48M
{
663
6.48M
  int j, r, tsgn;
664
6.48M
  isl_int t;
665
6.48M
  unsigned off = 2 + tab->M;
666
6.48M
667
6.48M
  isl_int_init(t);
668
6.48M
  r = -1;
669
75.5M
  for (j = tab->n_redundant; j < tab->n_row; 
++j69.0M
) {
670
69.0M
    if (var && 
j == var->index54.4M
)
671
5.43M
      continue;
672
63.6M
    if (!isl_tab_var_from_row(tab, j)->is_nonneg)
673
16.8M
      continue;
674
46.8M
    if (sgn * isl_int_sgn(tab->mat->row[j][off + c]) >= 0)
675
39.7M
      continue;
676
7.04M
    if (r < 0) {
677
4.75M
      r = j;
678
4.75M
      continue;
679
4.75M
    }
680
2.29M
    tsgn = sgn * row_cmp(tab, r, j, c, &t);
681
2.29M
    if (tsgn < 0 || 
(1.50M
tsgn == 01.50M
&&
682
1.50M
              
tab->row_var[j] < tab->row_var[r]700k
))
683
1.35M
      r = j;
684
2.29M
  }
685
6.48M
  isl_int_clear(t);
686
6.48M
  return r;
687
6.48M
}
688
689
/* Find a pivot (row and col) that will increase (sgn > 0) or decrease
690
 * (sgn < 0) the value of row variable var.
691
 * If not NULL, then skip_var is a row variable that should be ignored
692
 * while looking for a pivot row.  It is usually equal to var.
693
 *
694
 * As the given row in the tableau is of the form
695
 *
696
 *  x_r = a_r0 + \sum_i a_ri x_i
697
 *
698
 * we need to find a column such that the sign of a_ri is equal to "sgn"
699
 * (such that an increase in x_i will have the desired effect) or a
700
 * column with a variable that may attain negative values.
701
 * If a_ri is positive, then we need to move x_i in the same direction
702
 * to obtain the desired effect.  Otherwise, x_i has to move in the
703
 * opposite direction.
704
 */
705
static void find_pivot(struct isl_tab *tab,
706
  struct isl_tab_var *var, struct isl_tab_var *skip_var,
707
  int sgn, int *row, int *col)
708
8.44M
{
709
8.44M
  int j, r, c;
710
8.44M
  isl_int *tr;
711
8.44M
712
8.44M
  *row = *col = -1;
713
8.44M
714
8.44M
  isl_assert(tab->mat->ctx, var->is_row, return);
715
8.44M
  tr = tab->mat->row[var->index] + 2 + tab->M;
716
8.44M
717
8.44M
  c = -1;
718
64.1M
  for (j = tab->n_dead; j < tab->n_col; 
++j55.7M
) {
719
55.7M
    if (isl_int_is_zero(tr[j]))
720
55.7M
      
continue43.8M
;
721
11.8M
    if (isl_int_sgn(tr[j]) != sgn &&
722
11.8M
        
var_from_col(tab, j)->is_nonneg6.55M
)
723
4.44M
      continue;
724
7.43M
    if (c < 0 || 
tab->col_var[j] < tab->col_var[c]1.67M
)
725
6.18M
      c = j;
726
7.43M
  }
727
8.44M
  if (c < 0)
728
2.67M
    return;
729
5.76M
730
5.76M
  sgn *= isl_int_sgn(tr[c]);
731
5.76M
  r = pivot_row(tab, skip_var, sgn, c);
732
5.76M
  *row = r < 0 ? 
var->index1.73M
:
r4.02M
;
733
5.76M
  *col = c;
734
5.76M
}
735
736
/* Return 1 if row "row" represents an obviously redundant inequality.
737
 * This means
738
 *  - it represents an inequality or a variable
739
 *  - that is the sum of a non-negative sample value and a positive
740
 *    combination of zero or more non-negative constraints.
741
 */
742
int isl_tab_row_is_redundant(struct isl_tab *tab, int row)
743
34.3M
{
744
34.3M
  int i;
745
34.3M
  unsigned off = 2 + tab->M;
746
34.3M
747
34.3M
  if (tab->row_var[row] < 0 && 
!isl_tab_var_from_row(tab, row)->is_nonneg26.9M
)
748
1.87M
    return 0;
749
32.4M
750
32.4M
  if (isl_int_is_neg(tab->mat->row[row][1]))
751
32.4M
    
return 02.42M
;
752
30.0M
  if (tab->strict_redundant && 
isl_int_is_zero43
(tab->mat->row[row][1]))
753
30.0M
    
return 042
;
754
30.0M
  if (tab->M && 
isl_int_is_neg43.8k
(tab->mat->row[row][2]))
755
30.0M
    
return 02.37k
;
756
30.0M
757
117M
  
for (i = tab->n_dead; 30.0M
i < tab->n_col;
++i87.0M
) {
758
113M
    if (isl_int_is_zero(tab->mat->row[row][off + i]))
759
113M
      
continue78.8M
;
760
34.6M
    if (tab->col_var[i] >= 0)
761
13.8M
      return 0;
762
20.8M
    if (isl_int_is_neg(tab->mat->row[row][off + i]))
763
20.8M
      
return 012.1M
;
764
8.67M
    if (!var_from_col(tab, i)->is_nonneg)
765
497k
      return 0;
766
8.67M
  }
767
30.0M
  
return 13.58M
;
768
30.0M
}
769
770
static void swap_rows(struct isl_tab *tab, int row1, int row2)
771
3.47M
{
772
3.47M
  int t;
773
3.47M
  enum isl_tab_row_sign s;
774
3.47M
775
3.47M
  t = tab->row_var[row1];
776
3.47M
  tab->row_var[row1] = tab->row_var[row2];
777
3.47M
  tab->row_var[row2] = t;
778
3.47M
  isl_tab_var_from_row(tab, row1)->index = row1;
779
3.47M
  isl_tab_var_from_row(tab, row2)->index = row2;
780
3.47M
  tab->mat = isl_mat_swap_rows(tab->mat, row1, row2);
781
3.47M
782
3.47M
  if (!tab->row_sign)
783
3.46M
    return;
784
9.57k
  s = tab->row_sign[row1];
785
9.57k
  tab->row_sign[row1] = tab->row_sign[row2];
786
9.57k
  tab->row_sign[row2] = s;
787
9.57k
}
788
789
static isl_stat push_union(struct isl_tab *tab,
790
  enum isl_tab_undo_type type, union isl_tab_undo_val u) WARN_UNUSED;
791
792
/* Push record "u" onto the undo stack of "tab", provided "tab"
793
 * keeps track of undo information.
794
 *
795
 * If the record cannot be pushed, then mark the undo stack as invalid
796
 * such that a later rollback attempt will not try to undo earlier
797
 * records without having been able to undo the current record.
798
 */
799
static isl_stat push_union(struct isl_tab *tab,
800
  enum isl_tab_undo_type type, union isl_tab_undo_val u)
801
25.4M
{
802
25.4M
  struct isl_tab_undo *undo;
803
25.4M
804
25.4M
  if (!tab)
805
0
    return isl_stat_error;
806
25.4M
  if (!tab->need_undo)
807
19.4M
    return isl_stat_ok;
808
5.99M
809
5.99M
  undo = isl_alloc_type(tab->mat->ctx, struct isl_tab_undo);
810
5.99M
  if (!undo)
811
0
    goto error;
812
5.99M
  undo->type = type;
813
5.99M
  undo->u = u;
814
5.99M
  undo->next = tab->top;
815
5.99M
  tab->top = undo;
816
5.99M
817
5.99M
  return isl_stat_ok;
818
0
error:
819
0
  free_undo(tab);
820
0
  tab->top = NULL;
821
0
  return isl_stat_error;
822
5.99M
}
823
824
isl_stat isl_tab_push_var(struct isl_tab *tab,
825
  enum isl_tab_undo_type type, struct isl_tab_var *var)
826
24.8M
{
827
24.8M
  union isl_tab_undo_val u;
828
24.8M
  if (var->is_row)
829
24.3M
    u.var_index = tab->row_var[var->index];
830
475k
  else
831
475k
    u.var_index = tab->col_var[var->index];
832
24.8M
  return push_union(tab, type, u);
833
24.8M
}
834
835
isl_stat isl_tab_push(struct isl_tab *tab, enum isl_tab_undo_type type)
836
544k
{
837
544k
  union isl_tab_undo_val u = { 0 };
838
544k
  return push_union(tab, type, u);
839
544k
}
840
841
/* Push a record on the undo stack describing the current basic
842
 * variables, so that the this state can be restored during rollback.
843
 */
844
isl_stat isl_tab_push_basis(struct isl_tab *tab)
845
766
{
846
766
  int i;
847
766
  union isl_tab_undo_val u;
848
766
849
766
  u.col_var = isl_alloc_array(tab->mat->ctx, int, tab->n_col);
850
766
  if (tab->n_col && !u.col_var)
851
0
    return isl_stat_error;
852
8.75k
  
for (i = 0; 766
i < tab->n_col;
++i7.99k
)
853
7.99k
    u.col_var[i] = tab->col_var[i];
854
766
  return push_union(tab, isl_tab_undo_saved_basis, u);
855
766
}
856
857
isl_stat isl_tab_push_callback(struct isl_tab *tab,
858
  struct isl_tab_callback *callback)
859
21.7k
{
860
21.7k
  union isl_tab_undo_val u;
861
21.7k
  u.callback = callback;
862
21.7k
  return push_union(tab, isl_tab_undo_callback, u);
863
21.7k
}
864
865
struct isl_tab *isl_tab_init_samples(struct isl_tab *tab)
866
7.85k
{
867
7.85k
  if (!tab)
868
0
    return NULL;
869
7.85k
870
7.85k
  tab->n_sample = 0;
871
7.85k
  tab->n_outside = 0;
872
7.85k
  tab->samples = isl_mat_alloc(tab->mat->ctx, 1, 1 + tab->n_var);
873
7.85k
  if (!tab->samples)
874
0
    goto error;
875
7.85k
  tab->sample_index = isl_alloc_array(tab->mat->ctx, int, 1);
876
7.85k
  if (!tab->sample_index)
877
0
    goto error;
878
7.85k
  return tab;
879
0
error:
880
0
  isl_tab_free(tab);
881
0
  return NULL;
882
7.85k
}
883
884
int isl_tab_add_sample(struct isl_tab *tab, __isl_take isl_vec *sample)
885
12.1k
{
886
12.1k
  if (!tab || !sample)
887
0
    goto error;
888
12.1k
889
12.1k
  if (tab->n_sample + 1 > tab->samples->n_row) {
890
3.97k
    int *t = isl_realloc_array(tab->mat->ctx,
891
3.97k
          tab->sample_index, int, tab->n_sample + 1);
892
3.97k
    if (!t)
893
0
      goto error;
894
3.97k
    tab->sample_index = t;
895
3.97k
  }
896
12.1k
897
12.1k
  tab->samples = isl_mat_extend(tab->samples,
898
12.1k
        tab->n_sample + 1, tab->samples->n_col);
899
12.1k
  if (!tab->samples)
900
0
    goto error;
901
12.1k
902
12.1k
  isl_seq_cpy(tab->samples->row[tab->n_sample], sample->el, sample->size);
903
12.1k
  isl_vec_free(sample);
904
12.1k
  tab->sample_index[tab->n_sample] = tab->n_sample;
905
12.1k
  tab->n_sample++;
906
12.1k
907
12.1k
  return 0;
908
0
error:
909
0
  isl_vec_free(sample);
910
0
  return -1;
911
12.1k
}
912
913
struct isl_tab *isl_tab_drop_sample(struct isl_tab *tab, int s)
914
6.69k
{
915
6.69k
  if (s != tab->n_outside) {
916
4.16k
    int t = tab->sample_index[tab->n_outside];
917
4.16k
    tab->sample_index[tab->n_outside] = tab->sample_index[s];
918
4.16k
    tab->sample_index[s] = t;
919
4.16k
    isl_mat_swap_rows(tab->samples, tab->n_outside, s);
920
4.16k
  }
921
6.69k
  tab->n_outside++;
922
6.69k
  if (isl_tab_push(tab, isl_tab_undo_drop_sample) < 0) {
923
0
    isl_tab_free(tab);
924
0
    return NULL;
925
0
  }
926
6.69k
927
6.69k
  return tab;
928
6.69k
}
929
930
/* Record the current number of samples so that we can remove newer
931
 * samples during a rollback.
932
 */
933
isl_stat isl_tab_save_samples(struct isl_tab *tab)
934
28.9k
{
935
28.9k
  union isl_tab_undo_val u;
936
28.9k
937
28.9k
  if (!tab)
938
0
    return isl_stat_error;
939
28.9k
940
28.9k
  u.n = tab->n_sample;
941
28.9k
  return push_union(tab, isl_tab_undo_saved_samples, u);
942
28.9k
}
943
944
/* Mark row with index "row" as being redundant.
945
 * If we may need to undo the operation or if the row represents
946
 * a variable of the original problem, the row is kept,
947
 * but no longer considered when looking for a pivot row.
948
 * Otherwise, the row is simply removed.
949
 *
950
 * The row may be interchanged with some other row.  If it
951
 * is interchanged with a later row, return 1.  Otherwise return 0.
952
 * If the rows are checked in order in the calling function,
953
 * then a return value of 1 means that the row with the given
954
 * row number may now contain a different row that hasn't been checked yet.
955
 */
956
int isl_tab_mark_redundant(struct isl_tab *tab, int row)
957
4.75M
{
958
4.75M
  struct isl_tab_var *var = isl_tab_var_from_row(tab, row);
959
4.75M
  var->is_redundant = 1;
960
4.75M
  isl_assert(tab->mat->ctx, row >= tab->n_redundant, return -1);
961
4.75M
  if (tab->preserve || 
tab->need_undo3.41M
||
tab->row_var[row] >= 03.09M
) {
962
3.35M
    if (tab->row_var[row] >= 0 && 
!var->is_nonneg2.63M
) {
963
2.62M
      var->is_nonneg = 1;
964
2.62M
      if (isl_tab_push_var(tab, isl_tab_undo_nonneg, var) < 0)
965
0
        return -1;
966
3.35M
    }
967
3.35M
    if (row != tab->n_redundant)
968
2.45M
      swap_rows(tab, row, tab->n_redundant);
969
3.35M
    tab->n_redundant++;
970
3.35M
    return isl_tab_push_var(tab, isl_tab_undo_redundant, var);
971
3.35M
  } else {
972
1.39M
    if (row != tab->n_row - 1)
973
813k
      swap_rows(tab, row, tab->n_row - 1);
974
1.39M
    isl_tab_var_from_row(tab, tab->n_row - 1)->index = -1;
975
1.39M
    tab->n_row--;
976
1.39M
    return 1;
977
1.39M
  }
978
4.75M
}
979
980
/* Mark "tab" as a rational tableau.
981
 * If it wasn't marked as a rational tableau already and if we may
982
 * need to undo changes, then arrange for the marking to be undone
983
 * during the undo.
984
 */
985
int isl_tab_mark_rational(struct isl_tab *tab)
986
9.26k
{
987
9.26k
  if (!tab)
988
0
    return -1;
989
9.26k
  if (!tab->rational && 
tab->need_undo9.24k
)
990
9.24k
    if (isl_tab_push(tab, isl_tab_undo_rational) < 0)
991
0
      return -1;
992
9.26k
  tab->rational = 1;
993
9.26k
  return 0;
994
9.26k
}
995
996
isl_stat isl_tab_mark_empty(struct isl_tab *tab)
997
82.3k
{
998
82.3k
  if (!tab)
999
0
    return isl_stat_error;
1000
82.3k
  if (!tab->empty && 
tab->need_undo81.1k
)
1001
69.7k
    if (isl_tab_push(tab, isl_tab_undo_empty) < 0)
1002
0
      return isl_stat_error;
1003
82.3k
  tab->empty = 1;
1004
82.3k
  return isl_stat_ok;
1005
82.3k
}
1006
1007
int isl_tab_freeze_constraint(struct isl_tab *tab, int con)
1008
471k
{
1009
471k
  struct isl_tab_var *var;
1010
471k
1011
471k
  if (!tab)
1012
0
    return -1;
1013
471k
1014
471k
  var = &tab->con[con];
1015
471k
  if (var->frozen)
1016
0
    return 0;
1017
471k
  if (var->index < 0)
1018
43.0k
    return 0;
1019
428k
  var->frozen = 1;
1020
428k
1021
428k
  if (tab->need_undo)
1022
381k
    return isl_tab_push_var(tab, isl_tab_undo_freeze, var);
1023
47.0k
1024
47.0k
  return 0;
1025
47.0k
}
1026
1027
/* Update the rows signs after a pivot of "row" and "col", with "row_sgn"
1028
 * the original sign of the pivot element.
1029
 * We only keep track of row signs during PILP solving and in this case
1030
 * we only pivot a row with negative sign (meaning the value is always
1031
 * non-positive) using a positive pivot element.
1032
 *
1033
 * For each row j, the new value of the parametric constant is equal to
1034
 *
1035
 *  a_j0 - a_jc a_r0/a_rc
1036
 *
1037
 * where a_j0 is the original parametric constant, a_rc is the pivot element,
1038
 * a_r0 is the parametric constant of the pivot row and a_jc is the
1039
 * pivot column entry of the row j.
1040
 * Since a_r0 is non-positive and a_rc is positive, the sign of row j
1041
 * remains the same if a_jc has the same sign as the row j or if
1042
 * a_jc is zero.  In all other cases, we reset the sign to "unknown".
1043
 */
1044
static void update_row_sign(struct isl_tab *tab, int row, int col, int row_sgn)
1045
6.73M
{
1046
6.73M
  int i;
1047
6.73M
  struct isl_mat *mat = tab->mat;
1048
6.73M
  unsigned off = 2 + tab->M;
1049
6.73M
1050
6.73M
  if (!tab->row_sign)
1051
6.71M
    return;
1052
29.2k
1053
29.2k
  if (tab->row_sign[row] == 0)
1054
22.7k
    return;
1055
6.49k
  isl_assert(mat->ctx, row_sgn > 0, return);
1056
6.49k
  isl_assert(mat->ctx, tab->row_sign[row] == isl_tab_row_neg, return);
1057
6.49k
  tab->row_sign[row] = isl_tab_row_pos;
1058
60.8k
  for (i = 0; i < tab->n_row; 
++i54.3k
) {
1059
54.3k
    int s;
1060
54.3k
    if (i == row)
1061
6.49k
      continue;
1062
47.8k
    s = isl_int_sgn(mat->row[i][off + col]);
1063
47.8k
    if (!s)
1064
28.7k
      continue;
1065
19.0k
    if (!tab->row_sign[i])
1066
6.95k
      continue;
1067
12.0k
    if (s < 0 && 
tab->row_sign[i] == isl_tab_row_neg6.64k
)
1068
0
      continue;
1069
12.0k
    if (s > 0 && 
tab->row_sign[i] == isl_tab_row_pos5.42k
)
1070
5.42k
      continue;
1071
6.64k
    tab->row_sign[i] = isl_tab_row_unknown;
1072
6.64k
  }
1073
6.49k
}
1074
1075
/* Given a row number "row" and a column number "col", pivot the tableau
1076
 * such that the associated variables are interchanged.
1077
 * The given row in the tableau expresses
1078
 *
1079
 *  x_r = a_r0 + \sum_i a_ri x_i
1080
 *
1081
 * or
1082
 *
1083
 *  x_c = 1/a_rc x_r - a_r0/a_rc + sum_{i \ne r} -a_ri/a_rc
1084
 *
1085
 * Substituting this equality into the other rows
1086
 *
1087
 *  x_j = a_j0 + \sum_i a_ji x_i
1088
 *
1089
 * with a_jc \ne 0, we obtain
1090
 *
1091
 *  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 
1092
 *
1093
 * The tableau
1094
 *
1095
 *  n_rc/d_r    n_ri/d_r
1096
 *  n_jc/d_j    n_ji/d_j
1097
 *
1098
 * where i is any other column and j is any other row,
1099
 * is therefore transformed into
1100
 *
1101
 * s(n_rc)d_r/|n_rc|    -s(n_rc)n_ri/|n_rc|
1102
 * 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)
1103
 *
1104
 * The transformation is performed along the following steps
1105
 *
1106
 *  d_r/n_rc    n_ri/n_rc
1107
 *  n_jc/d_j    n_ji/d_j
1108
 *
1109
 *  s(n_rc)d_r/|n_rc| -s(n_rc)n_ri/|n_rc|
1110
 *  n_jc/d_j    n_ji/d_j
1111
 *
1112
 *  s(n_rc)d_r/|n_rc| -s(n_rc)n_ri/|n_rc|
1113
 *  n_jc/(|n_rc| d_j) n_ji/(|n_rc| d_j)
1114
 *
1115
 *  s(n_rc)d_r/|n_rc| -s(n_rc)n_ri/|n_rc|
1116
 *  n_jc/(|n_rc| d_j) (n_ji |n_rc|)/(|n_rc| d_j)
1117
 *
1118
 *  s(n_rc)d_r/|n_rc| -s(n_rc)n_ri/|n_rc|
1119
 *  n_jc/(|n_rc| d_j) (n_ji |n_rc| - s(n_rc)n_jc n_ri)/(|n_rc| d_j)
1120
 *
1121
 * s(n_rc)d_r/|n_rc|    -s(n_rc)n_ri/|n_rc|
1122
 * 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)
1123
 *
1124
 */
1125
int isl_tab_pivot(struct isl_tab *tab, int row, int col)
1126
6.73M
{
1127
6.73M
  int i, j;
1128
6.73M
  int sgn;
1129
6.73M
  int t;
1130
6.73M
  isl_ctx *ctx;
1131
6.73M
  struct isl_mat *mat = tab->mat;
1132
6.73M
  struct isl_tab_var *var;
1133
6.73M
  unsigned off = 2 + tab->M;
1134
6.73M
1135
6.73M
  ctx = isl_tab_get_ctx(tab);
1136
6.73M
  if (isl_ctx_next_operation(ctx) < 0)
1137
0
    return -1;
1138
6.73M
1139
6.73M
  isl_int_swap(mat->row[row][0], mat->row[row][off + col]);
1140
6.73M
  sgn = isl_int_sgn(mat->row[row][0]);
1141
6.73M
  if (sgn < 0) {
1142
3.82M
    isl_int_neg(mat->row[row][0], mat->row[row][0]);
1143
3.82M
    isl_int_neg(mat->row[row][off + col], mat->row[row][off + col]);
1144
3.82M
  } else
1145
24.5M
    
for (j = 0; 2.91M
j < off - 1 + tab->n_col;
++j21.6M
) {
1146
21.6M
      if (j == off - 1 + col)
1147
2.91M
        continue;
1148
18.7M
      isl_int_neg(mat->row[row][1 + j], mat->row[row][1 + j]);
1149
18.7M
    }
1150
6.73M
  if (!isl_int_is_one(mat->row[row][0]))
1151
6.73M
    
isl_seq_normalize(mat->ctx, mat->row[row], off + tab->n_col)1.09M
;
1152
83.1M
  for (i = 0; i < tab->n_row; 
++i76.3M
) {
1153
76.3M
    if (i == row)
1154
6.73M
      continue;
1155
69.6M
    if (isl_int_is_zero(mat->row[i][off + col]))
1156
69.6M
      
continue53.8M
;
1157
15.8M
    isl_int_mul(mat->row[i][0], mat->row[i][0], mat->row[row][0]);
1158
165M
    for (j = 0; j < off - 1 + tab->n_col; 
++j149M
) {
1159
149M
      if (j == off - 1 + col)
1160
15.8M
        continue;
1161
133M
      isl_int_mul(mat->row[i][1 + j],
1162
133M
            mat->row[i][1 + j], mat->row[row][0]);
1163
133M
      isl_int_addmul(mat->row[i][1 + j],
1164
133M
            mat->row[i][off + col], mat->row[row][1 + j]);
1165
133M
    }
1166
15.8M
    isl_int_mul(mat->row[i][off + col],
1167
15.8M
          mat->row[i][off + col], mat->row[row][off + col]);
1168
15.8M
    if (!isl_int_is_one(mat->row[i][0]))
1169
15.8M
      
isl_seq_normalize(mat->ctx, mat->row[i], off + tab->n_col)7.12M
;
1170
15.8M
  }
1171
6.73M
  t = tab->row_var[row];
1172
6.73M
  tab->row_var[row] = tab->col_var[col];
1173
6.73M
  tab->col_var[col] = t;
1174
6.73M
  var = isl_tab_var_from_row(tab, row);
1175
6.73M
  var->is_row = 1;
1176
6.73M
  var->index = row;
1177
6.73M
  var = var_from_col(tab, col);
1178
6.73M
  var->is_row = 0;
1179
6.73M
  var->index = col;
1180
6.73M
  update_row_sign(tab, row, col, sgn);
1181
6.73M
  if (tab->in_undo)
1182
274k
    return 0;
1183
62.4M
  
for (i = tab->n_redundant; 6.46M
i < tab->n_row;
++i56.0M
) {
1184
56.0M
    if (isl_int_is_zero(mat->row[i][off + col]))
1185
56.0M
      
continue36.7M
;
1186
19.3M
    if (!isl_tab_var_from_row(tab, i)->frozen &&
1187
19.3M
        
isl_tab_row_is_redundant(tab, i)18.9M
) {
1188
3.43M
      int redo = isl_tab_mark_redundant(tab, i);
1189
3.43M
      if (redo < 0)
1190
0
        return -1;
1191
3.43M
      if (redo)
1192
204k
        --i;
1193
3.43M
    }
1194
19.3M
  }
1195
6.46M
  return 0;
1196
6.46M
}
1197
1198
/* If "var" represents a column variable, then pivot is up (sgn > 0)
1199
 * or down (sgn < 0) to a row.  The variable is assumed not to be
1200
 * unbounded in the specified direction.
1201
 * If sgn = 0, then the variable is unbounded in both directions,
1202
 * and we pivot with any row we can find.
1203
 */
1204
static int to_row(struct isl_tab *tab, struct isl_tab_var *var, int sign) WARN_UNUSED;
1205
static int to_row(struct isl_tab *tab, struct isl_tab_var *var, int sign)
1206
3.36M
{
1207
3.36M
  int r;
1208
3.36M
  unsigned off = 2 + tab->M;
1209
3.36M
1210
3.36M
  if (var->is_row)
1211
2.82M
    return 0;
1212
534k
1213
534k
  if (sign == 0) {
1214
100k
    for (r = tab->n_redundant; r < tab->n_row; 
++r74.3k
)
1215
100k
      if (!isl_int_is_zero(tab->mat->row[r][off+var->index]))
1216
100k
        
break25.7k
;
1217
25.7k
    isl_assert(tab->mat->ctx, r < tab->n_row, return -1);
1218
508k
  } else {
1219
508k
    r = pivot_row(tab, NULL, sign, var->index);
1220
508k
    isl_assert(tab->mat->ctx, r >= 0, return -1);
1221
508k
  }
1222
534k
1223
534k
  return isl_tab_pivot(tab, r, var->index);
1224
534k
}
1225
1226
/* Check whether all variables that are marked as non-negative
1227
 * also have a non-negative sample value.  This function is not
1228
 * called from the current code but is useful during debugging.
1229
 */
1230
static void check_table(struct isl_tab *tab) __attribute__ ((unused));
1231
static void check_table(struct isl_tab *tab)
1232
0
{
1233
0
  int i;
1234
0
1235
0
  if (tab->empty)
1236
0
    return;
1237
0
  for (i = tab->n_redundant; i < tab->n_row; ++i) {
1238
0
    struct isl_tab_var *var;
1239
0
    var = isl_tab_var_from_row(tab, i);
1240
0
    if (!var->is_nonneg)
1241
0
      continue;
1242
0
    if (tab->M) {
1243
0
      isl_assert(tab->mat->ctx,
1244
0
        !isl_int_is_neg(tab->mat->row[i][2]), abort());
1245
0
      if (isl_int_is_pos(tab->mat->row[i][2]))
1246
0
        continue;
1247
0
    }
1248
0
    isl_assert(tab->mat->ctx, !isl_int_is_neg(tab->mat->row[i][1]),
1249
0
        abort());
1250
0
  }
1251
0
}
1252
1253
/* Return the sign of the maximal value of "var".
1254
 * If the sign is not negative, then on return from this function,
1255
 * the sample value will also be non-negative.
1256
 *
1257
 * If "var" is manifestly unbounded wrt positive values, we are done.
1258
 * Otherwise, we pivot the variable up to a row if needed
1259
 * Then we continue pivoting down until either
1260
 *  - no more down pivots can be performed
1261
 *  - the sample value is positive
1262
 *  - the variable is pivoted into a manifestly unbounded column
1263
 */
1264
static int sign_of_max(struct isl_tab *tab, struct isl_tab_var *var)
1265
2.40M
{
1266
2.40M
  int row, col;
1267
2.40M
1268
2.40M
  if (max_is_manifestly_unbounded(tab, var))
1269
145k
    return 1;
1270
2.26M
  if (to_row(tab, var, 1) < 0)
1271
0
    return -2;
1272
3.97M
  
while (2.26M
!isl_int_is_pos(tab->mat->row[var->index][1])) {
1273
3.14M
    find_pivot(tab, var, var, 1, &row, &col);
1274
3.14M
    if (row == -1)
1275
1.13M
      return isl_int_sgn(tab->mat->row[var->index][1]);
1276
2.00M
    if (isl_tab_pivot(tab, row, col) < 0)
1277
0
      return -2;
1278
2.00M
    if (!var->is_row) /* manifestly unbounded */
1279
293k
      return 1;
1280
2.00M
  }
1281
2.26M
  
return 1838k
;
1282
2.26M
}
1283
1284
int isl_tab_sign_of_max(struct isl_tab *tab, int con)
1285
130
{
1286
130
  struct isl_tab_var *var;
1287
130
1288
130
  if (!tab)
1289
0
    return -2;
1290
130
1291
130
  var = &tab->con[con];
1292
130
  isl_assert(tab->mat->ctx, !var->is_redundant, return -2);
1293
130
  isl_assert(tab->mat->ctx, !var->is_zero, return -2);
1294
130
1295
130
  return sign_of_max(tab, var);
1296
130
}
1297
1298
static int row_is_neg(struct isl_tab *tab, int row)
1299
10.1M
{
1300
10.1M
  if (!tab->M)
1301
10.1M
    return isl_int_is_neg(tab->mat->row[row][1]);
1302
0
  if (isl_int_is_pos(tab->mat->row[row][2]))
1303
0
    return 0;
1304
0
  if (isl_int_is_neg(tab->mat->row[row][2]))
1305
0
    return 1;
1306
0
  return isl_int_is_neg(tab->mat->row[row][1]);
1307
0
}
1308
1309
static int row_sgn(struct isl_tab *tab, int row)
1310
8.98M
{
1311
8.98M
  if (!tab->M)
1312
8.98M
    return isl_int_sgn(tab->mat->row[row][1]);
1313
0
  if (!isl_int_is_zero(tab->mat->row[row][2]))
1314
0
    return isl_int_sgn(tab->mat->row[row][2]);
1315
0
  else
1316
0
    return isl_int_sgn(tab->mat->row[row][1]);
1317
0
}
1318
1319
/* Perform pivots until the row variable "var" has a non-negative
1320
 * sample value or until no more upward pivots can be performed.
1321
 * Return the sign of the sample value after the pivots have been
1322
 * performed.
1323
 */
1324
static int restore_row(struct isl_tab *tab, struct isl_tab_var *var)
1325
9.69M
{
1326
9.69M
  int row, col;
1327
9.69M
1328
10.1M
  while (row_is_neg(tab, var->index)) {
1329
1.20M
    find_pivot(tab, var, var, 1, &row, &col);
1330
1.20M
    if (row == -1)
1331
78.9k
      break;
1332
1.12M
    if (isl_tab_pivot(tab, row, col) < 0)
1333
0
      return -2;
1334
1.12M
    if (!var->is_row) /* manifestly unbounded */
1335
707k
      return 1;
1336
1.12M
  }
1337
9.69M
  
return row_sgn(tab, var->index)8.98M
;
1338
9.69M
}
1339
1340
/* Perform pivots until we are sure that the row variable "var"
1341
 * can attain non-negative values.  After return from this
1342
 * function, "var" is still a row variable, but its sample
1343
 * value may not be non-negative, even if the function returns 1.
1344
 */
1345
static int at_least_zero(struct isl_tab *tab, struct isl_tab_var *var)
1346
183k
{
1347
183k
  int row, col;
1348
183k
1349
206k
  while (isl_int_is_neg(tab->mat->row[var->index][1])) {
1350
197k
    find_pivot(tab, var, var, 1, &row, &col);
1351
197k
    if (row == -1)
1352
105k
      break;
1353
91.3k
    if (row == var->index) /* manifestly unbounded */
1354
67.9k
      return 1;
1355
23.3k
    if (isl_tab_pivot(tab, row, col) < 0)
1356
0
      return -1;
1357
23.3k
  }
1358
183k
  
return !115k
isl_int_is_neg115k
(tab->mat->row[var->index][1]);
1359
183k
}
1360
1361
/* Return a negative value if "var" can attain negative values.
1362
 * Return a non-negative value otherwise.
1363
 *
1364
 * If "var" is manifestly unbounded wrt negative values, we are done.
1365
 * Otherwise, if var is in a column, we can pivot it down to a row.
1366
 * Then we continue pivoting down until either
1367
 *  - the pivot would result in a manifestly unbounded column
1368
 *    => we don't perform the pivot, but simply return -1
1369
 *  - no more down pivots can be performed
1370
 *  - the sample value is negative
1371
 * If the sample value becomes negative and the variable is supposed
1372
 * to be nonnegative, then we undo the last pivot.
1373
 * However, if the last pivot has made the pivoting variable
1374
 * obviously redundant, then it may have moved to another row.
1375
 * In that case we look for upward pivots until we reach a non-negative
1376
 * value again.
1377
 */
1378
static int sign_of_min(struct isl_tab *tab, struct isl_tab_var *var)
1379
70.3k
{
1380
70.3k
  int row, col;
1381
70.3k
  struct isl_tab_var *pivot_var = NULL;
1382
70.3k
1383
70.3k
  if (min_is_manifestly_unbounded(tab, var))
1384
0
    return -1;
1385
70.3k
  if (!var->is_row) {
1386
1.86k
    col = var->index;
1387
1.86k
    row = pivot_row(tab, NULL, -1, col);
1388
1.86k
    pivot_var = var_from_col(tab, col);
1389
1.86k
    if (isl_tab_pivot(tab, row, col) < 0)
1390
0
      return -2;
1391
1.86k
    if (var->is_redundant)
1392
317
      return 0;
1393
1.54k
    if (isl_int_is_neg(tab->mat->row[var->index][1])) {
1394
693
      if (var->is_nonneg) {
1395
693
        if (!pivot_var->is_redundant &&
1396
693
            pivot_var->index == row) {
1397
680
          if (isl_tab_pivot(tab, row, col) < 0)
1398
0
            return -2;
1399
13
        } else
1400
13
          if (restore_row(tab, var) < -1)
1401
0
            return -2;
1402
693
      }
1403
693
      return -1;
1404
693
    }
1405
1.54k
  }
1406
69.3k
  if (var->is_redundant)
1407
0
    return 0;
1408
106k
  
while (69.3k
!isl_int_is_neg(tab->mat->row[var->index][1])) {
1409
93.5k
    find_pivot(tab, var, var, -1, &row, &col);
1410
93.5k
    if (row == var->index)
1411
16.0k
      return -1;
1412
77.4k
    if (row == -1)
1413
38.8k
      return isl_int_sgn(tab->mat->row[var->index][1]);
1414
38.6k
    pivot_var = var_from_col(tab, col);
1415
38.6k
    if (isl_tab_pivot(tab, row, col) < 0)
1416
0
      return -2;
1417
38.6k
    if (var->is_redundant)
1418
1.14k
      return 0;
1419
38.6k
  }
1420
69.3k
  
if (13.3k
pivot_var13.3k
&&
var->is_nonneg13.3k
) {
1421
968
    /* pivot back to non-negative value */
1422
968
    if (!pivot_var->is_redundant && pivot_var->index == row) {
1423
955
      if (isl_tab_pivot(tab, row, col) < 0)
1424
0
        return -2;
1425
13
    } else
1426
13
      if (restore_row(tab, var) < -1)
1427
0
        return -2;
1428
13.3k
  }
1429
13.3k
  return -1;
1430
13.3k
}
1431
1432
static int row_at_most_neg_one(struct isl_tab *tab, int row)
1433
673k
{
1434
673k
  if (tab->M) {
1435
0
    if (isl_int_is_pos(tab->mat->row[row][2]))
1436
0
      return 0;
1437
0
    if (isl_int_is_neg(tab->mat->row[row][2]))
1438
0
      return 1;
1439
673k
  }
1440
673k
  return isl_int_is_neg(tab->mat->row[row][1]) &&
1441
673k
         
isl_int_abs_ge341k
(tab->mat->row[row][1],
1442
673k
            tab->mat->row[row][0]);
1443
673k
}
1444
1445
/* Return 1 if "var" can attain values <= -1.
1446
 * Return 0 otherwise.
1447
 *
1448
 * If the variable "var" is supposed to be non-negative (is_nonneg is set),
1449
 * then the sample value of "var" is assumed to be non-negative when the
1450
 * the function is called.  If 1 is returned then the constraint
1451
 * is not redundant and the sample value is made non-negative again before
1452
 * the function returns.
1453
 */
1454
int isl_tab_min_at_most_neg_one(struct isl_tab *tab, struct isl_tab_var *var)
1455
1.66M
{
1456
1.66M
  int row, col;
1457
1.66M
  struct isl_tab_var *pivot_var;
1458
1.66M
1459
1.66M
  if (min_is_manifestly_unbounded(tab, var))
1460
117
    return 1;
1461
1.66M
  if (!var->is_row) {
1462
214k
    col = var->index;
1463
214k
    row = pivot_row(tab, NULL, -1, col);
1464
214k
    pivot_var = var_from_col(tab, col);
1465
214k
    if (isl_tab_pivot(tab, row, col) < 0)
1466
0
      return -1;
1467
214k
    if (var->is_redundant)
1468
24.0k
      return 0;
1469
190k
    if (row_at_most_neg_one(tab, var->index)) {
1470
153k
      if (var->is_nonneg) {
1471
153k
        if (!pivot_var->is_redundant &&
1472
153k
            pivot_var->index == row) {
1473
144k
          if (isl_tab_pivot(tab, row, col) < 0)
1474
0
            return -1;
1475
8.43k
        } else
1476
8.43k
          if (restore_row(tab, var) < -1)
1477
0
            return -1;
1478
153k
      }
1479
153k
      return 1;
1480
153k
    }
1481
190k
  }
1482
1.48M
  if (var->is_redundant)
1483
12.1k
    return 0;
1484
1.82M
  
do 1.47M
{
1485
1.82M
    find_pivot(tab, var, var, -1, &row, &col);
1486
1.82M
    if (row == var->index) {
1487
636k
      if (var->is_nonneg && 
restore_row(tab, var) < -1594k
)
1488
0
        return -1;
1489
636k
      return 1;
1490
636k
    }
1491
1.18M
    if (row == -1)
1492
543k
      return 0;
1493
646k
    pivot_var = var_from_col(tab, col);
1494
646k
    if (isl_tab_pivot(tab, row, col) < 0)
1495
0
      return -1;
1496
646k
    if (var->is_redundant)
1497
164k
      return 0;
1498
482k
  } while (!row_at_most_neg_one(tab, var->index));
1499
1.47M
  
if (127k
var->is_nonneg127k
) {
1500
110k
    /* pivot back to non-negative value */
1501
110k
    if (!pivot_var->is_redundant && pivot_var->index == row)
1502
109k
      if (isl_tab_pivot(tab, row, col) < 0)
1503
0
        return -1;
1504
110k
    if (restore_row(tab, var) < -1)
1505
0
      return -1;
1506
127k
  }
1507
127k
  return 1;
1508
127k
}
1509
1510
/* Return 1 if "var" can attain values >= 1.
1511
 * Return 0 otherwise.
1512
 */
1513
static int at_least_one(struct isl_tab *tab, struct isl_tab_var *var)
1514
1.25M
{
1515
1.25M
  int row, col;
1516
1.25M
  isl_int *r;
1517
1.25M
1518
1.25M
  if (max_is_manifestly_unbounded(tab, var))
1519
438k
    return 1;
1520
811k
  if (to_row(tab, var, 1) < 0)
1521
0
    return -1;
1522
811k
  r = tab->mat->row[var->index];
1523
860k
  while (isl_int_lt(r[1], r[0])) {
1524
50.1k
    find_pivot(tab, var, var, 1, &row, &col);
1525
50.1k
    if (row == -1)
1526
856
      return isl_int_ge(r[1], r[0]);
1527
49.2k
    if (row == var->index) /* manifestly unbounded */
1528
132
      return 1;
1529
49.1k
    if (isl_tab_pivot(tab, row, col) < 0)
1530
0
      return -1;
1531
49.1k
  }
1532
811k
  
return 1810k
;
1533
811k
}
1534
1535
static void swap_cols(struct isl_tab *tab, int col1, int col2)
1536
1.53M
{
1537
1.53M
  int t;
1538
1.53M
  unsigned off = 2 + tab->M;
1539
1.53M
  t = tab->col_var[col1];
1540
1.53M
  tab->col_var[col1] = tab->col_var[col2];
1541
1.53M
  tab->col_var[col2] = t;
1542
1.53M
  var_from_col(tab, col1)->index = col1;
1543
1.53M
  var_from_col(tab, col2)->index = col2;
1544
1.53M
  tab->mat = isl_mat_swap_cols(tab->mat, off + col1, off + col2);
1545
1.53M
}
1546
1547
/* Mark column with index "col" as representing a zero variable.
1548
 * If we may need to undo the operation the column is kept,
1549
 * but no longer considered.
1550
 * Otherwise, the column is simply removed.
1551
 *
1552
 * The column may be interchanged with some other column.  If it
1553
 * is interchanged with a later column, return 1.  Otherwise return 0.
1554
 * If the columns are checked in order in the calling function,
1555
 * then a return value of 1 means that the column with the given
1556
 * column number may now contain a different column that
1557
 * hasn't been checked yet.
1558
 */
1559
int isl_tab_kill_col(struct isl_tab *tab, int col)
1560
2.05M
{
1561
2.05M
  var_from_col(tab, col)->is_zero = 1;
1562
2.05M
  if (tab->need_undo) {
1563
308k
    if (isl_tab_push_var(tab, isl_tab_undo_zero,
1564
308k
              var_from_col(tab, col)) < 0)
1565
0
      return -1;
1566
308k
    if (col != tab->n_dead)
1567
113k
      swap_cols(tab, col, tab->n_dead);
1568
308k
    tab->n_dead++;
1569
308k
    return 0;
1570
1.75M
  } else {
1571
1.75M
    if (col != tab->n_col - 1)
1572
1.42M
      swap_cols(tab, col, tab->n_col - 1);
1573
1.75M
    var_from_col(tab, tab->n_col - 1)->index = -1;
1574
1.75M
    tab->n_col--;
1575
1.75M
    return 1;
1576
1.75M
  }
1577
2.05M
}
1578
1579
static int row_is_manifestly_non_integral(struct isl_tab *tab, int row)
1580
4.76M
{
1581
4.76M
  unsigned off = 2 + tab->M;
1582
4.76M
1583
4.76M
  if (tab->M && 
!0
isl_int_eq0
(tab->mat->row[row][2],
1584
4.76M
          tab->mat->row[row][0]))
1585
4.76M
    
return 00
;
1586
4.76M
  if (isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
1587
4.76M
            tab->n_col - tab->n_dead) != -1)
1588
367k
    return 0;
1589
4.39M
1590
4.39M
  return !isl_int_is_divisible_by(tab->mat->row[row][1],
1591
4.39M
          tab->mat->row[row][0]);
1592
4.39M
}
1593
1594
/* For integer tableaus, check if any of the coordinates are stuck
1595
 * at a non-integral value.
1596
 */
1597
static int tab_is_manifestly_empty(struct isl_tab *tab)
1598
1.13M
{
1599
1.13M
  int i;
1600
1.13M
1601
1.13M
  if (tab->empty)
1602
4
    return 1;
1603
1.13M
  if (tab->rational)
1604
11.2k
    return 0;
1605
1.12M
1606
11.4M
  
for (i = 0; 1.12M
i < tab->n_var;
++i10.3M
) {
1607
10.3M
    if (!tab->var[i].is_row)
1608
5.57M
      continue;
1609
4.76M
    if (row_is_manifestly_non_integral(tab, tab->var[i].index))
1610
94
      return 1;
1611
4.76M
  }
1612
1.12M
1613
1.12M
  
return 01.12M
;
1614
1.12M
}
1615
1616
/* Row variable "var" is non-negative and cannot attain any values
1617
 * larger than zero.  This means that the coefficients of the unrestricted
1618
 * column variables are zero and that the coefficients of the non-negative
1619
 * column variables are zero or negative.
1620
 * Each of the non-negative variables with a negative coefficient can
1621
 * then also be written as the negative sum of non-negative variables
1622
 * and must therefore also be zero.
1623
 *
1624
 * If "temp_var" is set, then "var" is a temporary variable that
1625
 * will be removed after this function returns and for which
1626
 * no information is recorded on the undo stack.
1627
 * Do not add any undo records involving this variable in this case
1628
 * since the variable will have been removed before any future undo
1629
 * operations.  Also avoid marking the variable as redundant,
1630
 * since that either adds an undo record or needlessly removes the row
1631
 * (the caller will take care of removing the row).
1632
 */
1633
static isl_stat close_row(struct isl_tab *tab, struct isl_tab_var *var,
1634
  int temp_var) WARN_UNUSED;
1635
static isl_stat close_row(struct isl_tab *tab, struct isl_tab_var *var,
1636
  int temp_var)
1637
1.13M
{
1638
1.13M
  int j;
1639
1.13M
  struct isl_mat *mat = tab->mat;
1640
1.13M
  unsigned off = 2 + tab->M;
1641
1.13M
1642
1.13M
  if (!var->is_nonneg)
1643
1.13M
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_internal,
1644
1.13M
      "expecting non-negative variable",
1645
1.13M
      return isl_stat_error);
1646
1.13M
  var->is_zero = 1;
1647
1.13M
  if (!temp_var && 
tab->need_undo1.11M
)
1648
578
    if (isl_tab_push_var(tab, isl_tab_undo_zero, var) < 0)
1649
0
      return isl_stat_error;
1650
8.01M
  
for (j = tab->n_dead; 1.13M
j < tab->n_col;
++j6.88M
) {
1651
6.88M
    int recheck;
1652
6.88M
    if (isl_int_is_zero(mat->row[var->index][off + j]))
1653
6.88M
      
continue5.80M
;
1654
1.07M
    if (isl_int_is_pos(mat->row[var->index][off + j]))
1655
1.07M
      
isl_die0
(isl_tab_get_ctx(tab), isl_error_internal,
1656
1.07M
        "row cannot have positive coefficients",
1657
1.07M
        return isl_stat_error);
1658
1.07M
    recheck = isl_tab_kill_col(tab, j);
1659
1.07M
    if (recheck < 0)
1660
0
      return isl_stat_error;
1661
1.07M
    if (recheck)
1662
1.05M
      --j;
1663
1.07M
  }
1664
1.13M
  if (!temp_var && 
isl_tab_mark_redundant(tab, var->index) < 01.11M
)
1665
0
    return isl_stat_error;
1666
1.13M
  if (tab_is_manifestly_empty(tab) && 
isl_tab_mark_empty(tab) < 098
)
1667
0
    return isl_stat_error;
1668
1.13M
  return isl_stat_ok;
1669
1.13M
}
1670
1671
/* Add a constraint to the tableau and allocate a row for it.
1672
 * Return the index into the constraint array "con".
1673
 *
1674
 * This function assumes that at least one more row and at least
1675
 * one more element in the constraint array are available in the tableau.
1676
 */
1677
int isl_tab_allocate_con(struct isl_tab *tab)
1678
10.3M
{
1679
10.3M
  int r;
1680
10.3M
1681
10.3M
  isl_assert(tab->mat->ctx, tab->n_row < tab->mat->n_row, return -1);
1682
10.3M
  isl_assert(tab->mat->ctx, tab->n_con < tab->max_con, return -1);
1683
10.3M
1684
10.3M
  r = tab->n_con;
1685
10.3M
  tab->con[r].index = tab->n_row;
1686
10.3M
  tab->con[r].is_row = 1;
1687
10.3M
  tab->con[r].is_nonneg = 0;
1688
10.3M
  tab->con[r].is_zero = 0;
1689
10.3M
  tab->con[r].is_redundant = 0;
1690
10.3M
  tab->con[r].frozen = 0;
1691
10.3M
  tab->con[r].negated = 0;
1692
10.3M
  tab->row_var[tab->n_row] = ~r;
1693
10.3M
1694
10.3M
  tab->n_row++;
1695
10.3M
  tab->n_con++;
1696
10.3M
  if (isl_tab_push_var(tab, isl_tab_undo_allocate, &tab->con[r]) < 0)
1697
0
    return -1;
1698
10.3M
1699
10.3M
  return r;
1700
10.3M
}
1701
1702
/* Move the entries in tab->var up one position, starting at "first",
1703
 * creating room for an extra entry at position "first".
1704
 * Since some of the entries of tab->row_var and tab->col_var contain
1705
 * indices into this array, they have to be updated accordingly.
1706
 */
1707
static int var_insert_entry(struct isl_tab *tab, int first)
1708
8.41k
{
1709
8.41k
  int i;
1710
8.41k
1711
8.41k
  if (tab->n_var >= tab->max_var)
1712
8.41k
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_internal,
1713
8.41k
      "not enough room for new variable", return -1);
1714
8.41k
  if (first > tab->n_var)
1715
8.41k
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_internal,
1716
8.41k
      "invalid initial position", return -1);
1717
8.41k
1718
9.27k
  
for (i = tab->n_var - 1; 8.41k
i >= first;
--i860
) {
1719
860
    tab->var[i + 1] = tab->var[i];
1720
860
    if (tab->var[i + 1].is_row)
1721
580
      tab->row_var[tab->var[i + 1].index]++;
1722
280
    else
1723
280
      tab->col_var[tab->var[i + 1].index]++;
1724
860
  }
1725
8.41k
1726
8.41k
  tab->n_var++;
1727
8.41k
1728
8.41k
  return 0;
1729
8.41k
}
1730
1731
/* Drop the entry at position "first" in tab->var, moving all
1732
 * subsequent entries down.
1733
 * Since some of the entries of tab->row_var and tab->col_var contain
1734
 * indices into this array, they have to be updated accordingly.
1735
 */
1736
static int var_drop_entry(struct isl_tab *tab, int first)
1737
5.51k
{
1738
5.51k
  int i;
1739
5.51k
1740
5.51k
  if (first >= tab->n_var)
1741
5.51k
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_internal,
1742
5.51k
      "invalid initial position", return -1);
1743
5.51k
1744
5.51k
  tab->n_var--;
1745
5.51k
1746
6.09k
  for (i = first; i < tab->n_var; 
++i582
) {
1747
582
    tab->var[i] = tab->var[i + 1];
1748
582
    if (tab->var[i + 1].is_row)
1749
568
      tab->row_var[tab->var[i].index]--;
1750
14
    else
1751
14
      tab->col_var[tab->var[i].index]--;
1752
582
  }
1753
5.51k
1754
5.51k
  return 0;
1755
5.51k
}
1756
1757
/* Add a variable to the tableau at position "r" and allocate a column for it.
1758
 * Return the index into the variable array "var", i.e., "r",
1759
 * or -1 on error.
1760
 */
1761
int isl_tab_insert_var(struct isl_tab *tab, int r)
1762
8.41k
{
1763
8.41k
  int i;
1764
8.41k
  unsigned off = 2 + tab->M;
1765
8.41k
1766
8.41k
  isl_assert(tab->mat->ctx, tab->n_col < tab->mat->n_col, return -1);
1767
8.41k
1768
8.41k
  if (var_insert_entry(tab, r) < 0)
1769
0
    return -1;
1770
8.41k
1771
8.41k
  tab->var[r].index = tab->n_col;
1772
8.41k
  tab->var[r].is_row = 0;
1773
8.41k
  tab->var[r].is_nonneg = 0;
1774
8.41k
  tab->var[r].is_zero = 0;
1775
8.41k
  tab->var[r].is_redundant = 0;
1776
8.41k
  tab->var[r].frozen = 0;
1777
8.41k
  tab->var[r].negated = 0;
1778
8.41k
  tab->col_var[tab->n_col] = r;
1779
8.41k
1780
58.0k
  for (i = 0; i < tab->n_row; 
++i49.6k
)
1781
49.6k
    isl_int_set_si(tab->mat->row[i][off + tab->n_col], 0);
1782
8.41k
1783
8.41k
  tab->n_col++;
1784
8.41k
  if (isl_tab_push_var(tab, isl_tab_undo_allocate, &tab->var[r]) < 0)
1785
0
    return -1;
1786
8.41k
1787
8.41k
  return r;
1788
8.41k
}
1789
1790
/* Add a variable to the tableau and allocate a column for it.
1791
 * Return the index into the variable array "var".
1792
 */
1793
int isl_tab_allocate_var(struct isl_tab *tab)
1794
0
{
1795
0
  if (!tab)
1796
0
    return -1;
1797
0
1798
0
  return isl_tab_insert_var(tab, tab->n_var);
1799
0
}
1800
1801
/* Add a row to the tableau.  The row is given as an affine combination
1802
 * of the original variables and needs to be expressed in terms of the
1803
 * column variables.
1804
 *
1805
 * This function assumes that at least one more row and at least
1806
 * one more element in the constraint array are available in the tableau.
1807
 *
1808
 * We add each term in turn.
1809
 * If r = n/d_r is the current sum and we need to add k x, then
1810
 *  if x is a column variable, we increase the numerator of
1811
 *    this column by k d_r
1812
 *  if x = f/d_x is a row variable, then the new representation of r is
1813
 *
1814
 *     n    k f   d_x/g n + d_r/g k f   m/d_r n + m/d_g k f
1815
 *    --- + --- = ------------------- = -------------------
1816
 *    d_r   d_r        d_r d_x/g                m
1817
 *
1818
 *  with g the gcd of d_r and d_x and m the lcm of d_r and d_x.
1819
 *
1820
 * If tab->M is set, then, internally, each variable x is represented
1821
 * as x' - M.  We then also need no subtract k d_r from the coefficient of M.
1822
 */
1823
int isl_tab_add_row(struct isl_tab *tab, isl_int *line)
1824
10.3M
{
1825
10.3M
  int i;
1826
10.3M
  int r;
1827
10.3M
  isl_int *row;
1828
10.3M
  isl_int a, b;
1829
10.3M
  unsigned off = 2 + tab->M;
1830
10.3M
1831
10.3M
  r = isl_tab_allocate_con(tab);
1832
10.3M
  if (r < 0)
1833
0
    return -1;
1834
10.3M
1835
10.3M
  isl_int_init(a);
1836
10.3M
  isl_int_init(b);
1837
10.3M
  row = tab->mat->row[tab->con[r].index];
1838
10.3M
  isl_int_set_si(row[0], 1);
1839
10.3M
  isl_int_set(row[1], line[0]);
1840
10.3M
  isl_seq_clr(row + 2, tab->M + tab->n_col);
1841
93.7M
  for (i = 0; i < tab->n_var; 
++i83.3M
) {
1842
83.3M
    if (tab->var[i].is_zero)
1843
0
      continue;
1844
83.3M
    if (tab->var[i].is_row) {
1845
22.0M
      isl_int_lcm(a,
1846
22.0M
        row[0], tab->mat->row[tab->var[i].index][0]);
1847
22.0M
      isl_int_swap(a, row[0]);
1848
22.0M
      isl_int_divexact(a, row[0], a);
1849
22.0M
      isl_int_divexact(b,
1850
22.0M
        row[0], tab->mat->row[tab->var[i].index][0]);
1851
22.0M
      isl_int_mul(b, b, line[1 + i]);
1852
22.0M
      isl_seq_combine(row + 1, a, row + 1,
1853
22.0M
          b, tab->mat->row[tab->var[i].index] + 1,
1854
22.0M
          1 + tab->M + tab->n_col);
1855
22.0M
    } else
1856
83.3M
      
isl_int_addmul61.3M
(row[off + tab->var[i].index],
1857
83.3M
              line[1 + i], row[0]);
1858
83.3M
    if (tab->M && 
i >= tab->n_param355k
&&
i < tab->n_var - tab->n_div154k
)
1859
83.3M
      
isl_int_submul151k
(row[2], line[1 + i], row[0]);
1860
83.3M
  }
1861
10.3M
  isl_seq_normalize(tab->mat->ctx, row, off + tab->n_col);
1862
10.3M
  isl_int_clear(a);
1863
10.3M
  isl_int_clear(b);
1864
10.3M
1865
10.3M
  if (tab->row_sign)
1866
42.5k
    tab->row_sign[tab->con[r].index] = isl_tab_row_unknown;
1867
10.3M
1868
10.3M
  return r;
1869
10.3M
}
1870
1871
static isl_stat drop_row(struct isl_tab *tab, int row)
1872
2.17M
{
1873
2.17M
  isl_assert(tab->mat->ctx, ~tab->row_var[row] == tab->n_con - 1,
1874
2.17M
    return isl_stat_error);
1875
2.17M
  if (row != tab->n_row - 1)
1876
204k
    swap_rows(tab, row, tab->n_row - 1);
1877
2.17M
  tab->n_row--;
1878
2.17M
  tab->n_con--;
1879
2.17M
  return isl_stat_ok;
1880
2.17M
}
1881
1882
/* Drop the variable in column "col" along with the column.
1883
 * The column is removed first because it may need to be moved
1884
 * into the last position and this process requires
1885
 * the contents of the col_var array in a state
1886
 * before the removal of the variable.
1887
 */
1888
static isl_stat drop_col(struct isl_tab *tab, int col)
1889
5.51k
{
1890
5.51k
  int var;
1891
5.51k
1892
5.51k
  var = tab->col_var[col];
1893
5.51k
  if (col != tab->n_col - 1)
1894
1.39k
    swap_cols(tab, col, tab->n_col - 1);
1895
5.51k
  tab->n_col--;
1896
5.51k
  if (var_drop_entry(tab, var) < 0)
1897
0
    return isl_stat_error;
1898
5.51k
  return isl_stat_ok;
1899
5.51k
}
1900
1901
/* Add inequality "ineq" and check if it conflicts with the
1902
 * previously added constraints or if it is obviously redundant.
1903
 *
1904
 * This function assumes that at least one more row and at least
1905
 * one more element in the constraint array are available in the tableau.
1906
 */
1907
isl_stat isl_tab_add_ineq(struct isl_tab *tab, isl_int *ineq)
1908
7.74M
{
1909
7.74M
  int r;
1910
7.74M
  int sgn;
1911
7.74M
  isl_int cst;
1912
7.74M
1913
7.74M
  if (!tab)
1914
0
    return isl_stat_error;
1915
7.74M
  if (tab->bmap) {
1916
449k
    struct isl_basic_map *bmap = tab->bmap;
1917
449k
1918
449k
    isl_assert(tab->mat->ctx, tab->n_eq == bmap->n_eq,
1919
449k
      return isl_stat_error);
1920
449k
    isl_assert(tab->mat->ctx,
1921
449k
          tab->n_con == bmap->n_eq + bmap->n_ineq,
1922
449k
          return isl_stat_error);
1923
449k
    tab->bmap = isl_basic_map_add_ineq(tab->bmap, ineq);
1924
449k
    if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1925
0
      return isl_stat_error;
1926
449k
    if (!tab->bmap)
1927
0
      return isl_stat_error;
1928
7.74M
  }
1929
7.74M
  if (tab->cone) {
1930
4.68k
    isl_int_init(cst);
1931
4.68k
    isl_int_set_si(cst, 0);
1932
4.68k
    isl_int_swap(ineq[0], cst);
1933
4.68k
  }
1934
7.74M
  r = isl_tab_add_row(tab, ineq);
1935
7.74M
  if (tab->cone) {
1936
4.68k
    isl_int_swap(ineq[0], cst);
1937
4.68k
    isl_int_clear(cst);
1938
4.68k
  }
1939
7.74M
  if (r < 0)
1940
0
    return isl_stat_error;
1941
7.74M
  tab->con[r].is_nonneg = 1;
1942
7.74M
  if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1943
0
    return isl_stat_error;
1944
7.74M
  if (isl_tab_row_is_redundant(tab, tab->con[r].index)) {
1945
148k
    if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0)
1946
0
      return isl_stat_error;
1947
148k
    return isl_stat_ok;
1948
148k
  }
1949
7.59M
1950
7.59M
  sgn = restore_row(tab, &tab->con[r]);
1951
7.59M
  if (sgn < -1)
1952
0
    return isl_stat_error;
1953
7.59M
  if (sgn < 0)
1954
78.9k
    return isl_tab_mark_empty(tab);
1955
7.51M
  if (tab->con[r].is_row && 
isl_tab_row_is_redundant(tab, tab->con[r].index)6.81M
)
1956
0
    if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0)
1957
0
      return isl_stat_error;
1958
7.51M
  return isl_stat_ok;
1959
7.51M
}
1960
1961
/* Pivot a non-negative variable down until it reaches the value zero
1962
 * and then pivot the variable into a column position.
1963
 */
1964
static int to_col(struct isl_tab *tab, struct isl_tab_var *var) WARN_UNUSED;
1965
static int to_col(struct isl_tab *tab, struct isl_tab_var *var)
1966
289k
{
1967
289k
  int i;
1968
289k
  int row, col;
1969
289k
  unsigned off = 2 + tab->M;
1970
289k
1971
289k
  if (!var->is_row)
1972
214
    return 0;
1973
289k
1974
325k
  
while (289k
isl_int_is_pos(tab->mat->row[var->index][1])) {
1975
320k
    find_pivot(tab, var, NULL, -1, &row, &col);
1976
320k
    isl_assert(tab->mat->ctx, row != -1, return -1);
1977
320k
    if (isl_tab_pivot(tab, row, col) < 0)
1978
0
      return -1;
1979
320k
    if (!var->is_row)
1980
285k
      return 0;
1981
320k
  }
1982
289k
1983
289k
  
for (i = tab->n_dead; 4.46k
i < tab->n_col9.23k
;
++i4.77k
)
1984
9.23k
    if (!isl_int_is_zero(tab->mat->row[var->index][off + i]))
1985
9.23k
      
break4.46k
;
1986
4.46k
1987
4.46k
  isl_assert(tab->mat->ctx, i < tab->n_col, return -1);
1988
4.46k
  if (isl_tab_pivot(tab, var->index, i) < 0)
1989
0
    return -1;
1990
4.46k
1991
4.46k
  return 0;
1992
4.46k
}
1993
1994
/* We assume Gaussian elimination has been performed on the equalities.
1995
 * The equalities can therefore never conflict.
1996
 * Adding the equalities is currently only really useful for a later call
1997
 * to isl_tab_ineq_type.
1998
 *
1999
 * This function assumes that at least one more row and at least
2000
 * one more element in the constraint array are available in the tableau.
2001
 */
2002
static struct isl_tab *add_eq(struct isl_tab *tab, isl_int *eq)
2003
667k
{
2004
667k
  int i;
2005
667k
  int r;
2006
667k
2007
667k
  if (!tab)
2008
0
    return NULL;
2009
667k
  r = isl_tab_add_row(tab, eq);
2010
667k
  if (r < 0)
2011
0
    goto error;
2012
667k
2013
667k
  r = tab->con[r].index;
2014
667k
  i = isl_seq_first_non_zero(tab->mat->row[r] + 2 + tab->M + tab->n_dead,
2015
667k
          tab->n_col - tab->n_dead);
2016
667k
  isl_assert(tab->mat->ctx, i >= 0, goto error);
2017
667k
  i += tab->n_dead;
2018
667k
  if (isl_tab_pivot(tab, r, i) < 0)
2019
0
    goto error;
2020
667k
  if (isl_tab_kill_col(tab, i) < 0)
2021
0
    goto error;
2022
667k
  tab->n_eq++;
2023
667k
2024
667k
  return tab;
2025
0
error:
2026
0
  isl_tab_free(tab);
2027
0
  return NULL;
2028
667k
}
2029
2030
/* Does the sample value of row "row" of "tab" involve the big parameter,
2031
 * if any?
2032
 */
2033
static int row_is_big(struct isl_tab *tab, int row)
2034
103k
{
2035
103k
  return tab->M && 
!0
isl_int_is_zero0
(tab->mat->row[row][2]);
2036
103k
}
2037
2038
static int row_is_manifestly_zero(struct isl_tab *tab, int row)
2039
302k
{
2040
302k
  unsigned off = 2 + tab->M;
2041
302k
2042
302k
  if (!isl_int_is_zero(tab->mat->row[row][1]))
2043
302k
    
return 0285k
;
2044
16.8k
  if (row_is_big(tab, row))
2045
0
    return 0;
2046
16.8k
  return isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
2047
16.8k
          tab->n_col - tab->n_dead) == -1;
2048
16.8k
}
2049
2050
/* Add an equality that is known to be valid for the given tableau.
2051
 *
2052
 * This function assumes that at least one more row and at least
2053
 * one more element in the constraint array are available in the tableau.
2054
 */
2055
int isl_tab_add_valid_eq(struct isl_tab *tab, isl_int *eq)
2056
285k
{
2057
285k
  struct isl_tab_var *var;
2058
285k
  int r;
2059
285k
2060
285k
  if (!tab)
2061
0
    return -1;
2062
285k
  r = isl_tab_add_row(tab, eq);
2063
285k
  if (r < 0)
2064
0
    return -1;
2065
285k
2066
285k
  var = &tab->con[r];
2067
285k
  r = var->index;
2068
285k
  if (row_is_manifestly_zero(tab, r)) {
2069
1.48k
    var->is_zero = 1;
2070
1.48k
    if (isl_tab_mark_redundant(tab, r) < 0)
2071
0
      return -1;
2072
1.48k
    return 0;
2073
1.48k
  }
2074
283k
2075
283k
  if (isl_int_is_neg(tab->mat->row[r][1])) {
2076
120k
    isl_seq_neg(tab->mat->row[r] + 1, tab->mat->row[r] + 1,
2077
120k
          1 + tab->n_col);
2078
120k
    var->negated = 1;
2079
120k
  }
2080
283k
  var->is_nonneg = 1;
2081
283k
  if (to_col(tab, var) < 0)
2082
0
    return -1;
2083
283k
  var->is_nonneg = 0;
2084
283k
  if (isl_tab_kill_col(tab, var->index) < 0)
2085
0
    return -1;
2086
283k
2087
283k
  return 0;
2088
283k
}
2089
2090
/* Add a zero row to "tab" and return the corresponding index
2091
 * in the constraint array.
2092
 *
2093
 * This function assumes that at least one more row and at least
2094
 * one more element in the constraint array are available in the tableau.
2095
 */
2096
static int add_zero_row(struct isl_tab *tab)
2097
4.74k
{
2098
4.74k
  int r;
2099
4.74k
  isl_int *row;
2100
4.74k
2101
4.74k
  r = isl_tab_allocate_con(tab);
2102
4.74k
  if (r < 0)
2103
0
    return -1;
2104
4.74k
2105
4.74k
  row = tab->mat->row[tab->con[r].index];
2106
4.74k
  isl_seq_clr(row + 1, 1 + tab->M + tab->n_col);
2107
4.74k
  isl_int_set_si(row[0], 1);
2108
4.74k
2109
4.74k
  return r;
2110
4.74k
}
2111
2112
/* Add equality "eq" and check if it conflicts with the
2113
 * previously added constraints or if it is obviously redundant.
2114
 *
2115
 * This function assumes that at least one more row and at least
2116
 * one more element in the constraint array are available in the tableau.
2117
 * If tab->bmap is set, then two rows are needed instead of one.
2118
 */
2119
int isl_tab_add_eq(struct isl_tab *tab, isl_int *eq)
2120
17.3k
{
2121
17.3k
  struct isl_tab_undo *snap = NULL;
2122
17.3k
  struct isl_tab_var *var;
2123
17.3k
  int r;
2124
17.3k
  int row;
2125
17.3k
  int sgn;
2126
17.3k
  isl_int cst;
2127
17.3k
2128
17.3k
  if (!tab)
2129
0
    return -1;
2130
17.3k
  isl_assert(tab->mat->ctx, !tab->M, return -1);
2131
17.3k
2132
17.3k
  if (tab->need_undo)
2133
16.7k
    snap = isl_tab_snap(tab);
2134
17.3k
2135
17.3k
  if (tab->cone) {
2136
1.18k
    isl_int_init(cst);
2137
1.18k
    isl_int_set_si(cst, 0);
2138
1.18k
    isl_int_swap(eq[0], cst);
2139
1.18k
  }
2140
17.3k
  r = isl_tab_add_row(tab, eq);
2141
17.3k
  if (tab->cone) {
2142
1.18k
    isl_int_swap(eq[0], cst);
2143
1.18k
    isl_int_clear(cst);
2144
1.18k
  }
2145
17.3k
  if (r < 0)
2146
0
    return -1;
2147
17.3k
2148
17.3k
  var = &tab->con[r];
2149
17.3k
  row = var->index;
2150
17.3k
  if (row_is_manifestly_zero(tab, row)) {
2151
11.3k
    if (snap)
2152
11.2k
      return isl_tab_rollback(tab, snap);
2153
50
    return drop_row(tab, row);
2154
50
  }
2155
6.00k
2156
6.00k
  if (tab->bmap) {
2157
4.74k
    tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);
2158
4.74k
    if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
2159
0
      return -1;
2160
4.74k
    isl_seq_neg(eq, eq, 1 + tab->n_var);
2161
4.74k
    tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);
2162
4.74k
    isl_seq_neg(eq, eq, 1 + tab->n_var);
2163
4.74k
    if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
2164
0
      return -1;
2165
4.74k
    if (!tab->bmap)
2166
0
      return -1;
2167
4.74k
    if (add_zero_row(tab) < 0)
2168
0
      return -1;
2169
6.00k
  }
2170
6.00k
2171
6.00k
  sgn = isl_int_sgn(tab->mat->row[row][1]);
2172
6.00k
2173
6.00k
  if (sgn > 0) {
2174
295
    isl_seq_neg(tab->mat->row[row] + 1, tab->mat->row[row] + 1,
2175
295
          1 + tab->n_col);
2176
295
    var->negated = 1;
2177
295
    sgn = -1;
2178
295
  }
2179
6.00k
2180
6.00k
  if (sgn < 0) {
2181
4.68k
    sgn = sign_of_max(tab, var);
2182
4.68k
    if (sgn < -1)
2183
0
      return -1;
2184
4.68k
    if (sgn < 0) {
2185
0
      if (isl_tab_mark_empty(tab) < 0)
2186
0
        return -1;
2187
0
      return 0;
2188
0
    }
2189
4.68k
  }
2190
6.00k
2191
6.00k
  var->is_nonneg = 1;
2192
6.00k
  if (to_col(tab, var) < 0)
2193
0
    return -1;
2194
6.00k
  var->is_nonneg = 0;
2195
6.00k
  if (isl_tab_kill_col(tab, var->index) < 0)
2196
0
    return -1;
2197
6.00k
2198
6.00k
  return 0;
2199
6.00k
}
2200
2201
/* Construct and return an inequality that expresses an upper bound
2202
 * on the given div.
2203
 * In particular, if the div is given by
2204
 *
2205
 *  d = floor(e/m)
2206
 *
2207
 * then the inequality expresses
2208
 *
2209
 *  m d <= e
2210
 */
2211
static struct isl_vec *ineq_for_div(struct isl_basic_map *bmap, unsigned div)
2212
4.14k
{
2213
4.14k
  unsigned total;
2214
4.14k
  unsigned div_pos;
2215
4.14k
  struct isl_vec *ineq;
2216
4.14k
2217
4.14k
  if (!bmap)
2218
0
    return NULL;
2219
4.14k
2220
4.14k
  total = isl_basic_map_total_dim(bmap);
2221
4.14k
  div_pos = 1 + total - bmap->n_div + div;
2222
4.14k
2223
4.14k
  ineq = isl_vec_alloc(bmap->ctx, 1 + total);
2224
4.14k
  if (!ineq)
2225
0
    return NULL;
2226
4.14k
2227
4.14k
  isl_seq_cpy(ineq->el, bmap->div[div] + 1, 1 + total);
2228
4.14k
  isl_int_neg(ineq->el[div_pos], bmap->div[div][0]);
2229
4.14k
  return ineq;
2230
4.14k
}
2231
2232
/* For a div d = floor(f/m), add the constraints
2233
 *
2234
 *    f - m d >= 0
2235
 *    -(f-(m-1)) + m d >= 0
2236
 *
2237
 * Note that the second constraint is the negation of
2238
 *
2239
 *    f - m d >= m
2240
 *
2241
 * If add_ineq is not NULL, then this function is used
2242
 * instead of isl_tab_add_ineq to effectively add the inequalities.
2243
 *
2244
 * This function assumes that at least two more rows and at least
2245
 * two more elements in the constraint array are available in the tableau.
2246
 */
2247
static isl_stat add_div_constraints(struct isl_tab *tab, unsigned div,
2248
  isl_stat (*add_ineq)(void *user, isl_int *), void *user)
2249
4.14k
{
2250
4.14k
  unsigned total;
2251
4.14k
  unsigned div_pos;
2252
4.14k
  struct isl_vec *ineq;
2253
4.14k
2254
4.14k
  total = isl_basic_map_total_dim(tab->bmap);
2255
4.14k
  div_pos = 1 + total - tab->bmap->n_div + div;
2256
4.14k
2257
4.14k
  ineq = ineq_for_div(tab->bmap, div);
2258
4.14k
  if (!ineq)
2259
0
    goto error;
2260
4.14k
2261
4.14k
  if (add_ineq) {
2262
759
    if (add_ineq(user, ineq->el) < 0)
2263
0
      goto error;
2264
3.38k
  } else {
2265
3.38k
    if (isl_tab_add_ineq(tab, ineq->el) < 0)
2266
0
      goto error;
2267
4.14k
  }
2268
4.14k
2269
4.14k
  isl_seq_neg(ineq->el, tab->bmap->div[div] + 1, 1 + total);
2270
4.14k
  isl_int_set(ineq->el[div_pos], tab->bmap->div[div][0]);
2271
4.14k
  isl_int_add(ineq->el[0], ineq->el[0], ineq->el[div_pos]);
2272
4.14k
  isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
2273
4.14k
2274
4.14k
  if (add_ineq) {
2275
759
    if (add_ineq(user, ineq->el) < 0)
2276
0
      goto error;
2277
3.38k
  } else {
2278
3.38k
    if (isl_tab_add_ineq(tab, ineq->el) < 0)
2279
0
      goto error;
2280
4.14k
  }
2281
4.14k
2282
4.14k
  isl_vec_free(ineq);
2283
4.14k
2284
4.14k
  return isl_stat_ok;
2285
0
error:
2286
0
  isl_vec_free(ineq);
2287
0
  return isl_stat_error;
2288
4.14k
}
2289
2290
/* Check whether the div described by "div" is obviously non-negative.
2291
 * If we are using a big parameter, then we will encode the div
2292
 * as div' = M + div, which is always non-negative.
2293
 * Otherwise, we check whether div is a non-negative affine combination
2294
 * of non-negative variables.
2295
 */
2296
static int div_is_nonneg(struct isl_tab *tab, __isl_keep isl_vec *div)
2297
4.14k
{
2298
4.14k
  int i;
2299
4.14k
2300
4.14k
  if (tab->M)
2301
0
    return 1;
2302
4.14k
2303
4.14k
  if (isl_int_is_neg(div->el[1]))
2304
4.14k
    
return 0700
;
2305
3.44k
2306
11.0k
  
for (i = 0; 3.44k
i < tab->n_var;
++i7.55k
) {
2307
10.3k
    if (isl_int_is_neg(div->el[2 + i]))
2308
10.3k
      
return 0289
;
2309
10.1k
    if (isl_int_is_zero(div->el[2 + i]))
2310
10.1k
      
continue6.77k
;
2311
3.32k
    if (!tab->var[i].is_nonneg)
2312
2.54k
      return 0;
2313
3.32k
  }
2314
3.44k
2315
3.44k
  
return 1612
;
2316
3.44k
}
2317
2318
/* Insert an extra div, prescribed by "div", to the tableau and
2319
 * the associated bmap (which is assumed to be non-NULL).
2320
 * The extra integer division is inserted at (tableau) position "pos".
2321
 * Return "pos" or -1 if an error occurred.
2322
 *
2323
 * If add_ineq is not NULL, then this function is used instead
2324
 * of isl_tab_add_ineq to add the div constraints.
2325
 * This complication is needed because the code in isl_tab_pip
2326
 * wants to perform some extra processing when an inequality
2327
 * is added to the tableau.
2328
 */
2329
int isl_tab_insert_div(struct isl_tab *tab, int pos, __isl_keep isl_vec *div,
2330
  isl_stat (*add_ineq)(void *user, isl_int *), void *user)
2331
4.14k
{
2332
4.14k
  int r;
2333
4.14k
  int nonneg;
2334
4.14k
  int n_div, o_div;
2335
4.14k
2336
4.14k
  if (!tab || !div)
2337
0
    return -1;
2338
4.14k
2339
4.14k
  if (div->size != 1 + 1 + tab->n_var)
2340
4.14k
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_invalid,
2341
4.14k
      "unexpected size", return -1);
2342
4.14k
2343
4.14k
  isl_assert(tab->mat->ctx, tab->bmap, return -1);
2344
4.14k
  n_div = isl_basic_map_dim(tab->bmap, isl_dim_div);
2345
4.14k
  o_div = tab->n_var - n_div;
2346
4.14k
  if (pos < o_div || pos > tab->n_var)
2347
4.14k
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_invalid,
2348
4.14k
      "invalid position", return -1);
2349
4.14k
2350
4.14k
  nonneg = div_is_nonneg(tab, div);
2351
4.14k
2352
4.14k
  if (isl_tab_extend_cons(tab, 3) < 0)
2353
0
    return -1;
2354
4.14k
  if (isl_tab_extend_vars(tab, 1) < 0)
2355
0
    return -1;
2356
4.14k
  r = isl_tab_insert_var(tab, pos);
2357
4.14k
  if (r < 0)
2358
0
    return -1;
2359
4.14k
2360
4.14k
  if (nonneg)
2361
612
    tab->var[r].is_nonneg = 1;
2362
4.14k
2363
4.14k
  tab->bmap = isl_basic_map_insert_div(tab->bmap, pos - o_div, div);
2364
4.14k
  if (!tab->bmap)
2365
0
    return -1;
2366
4.14k
  if (isl_tab_push_var(tab, isl_tab_undo_bmap_div, &tab->var[r]) < 0)
2367
0
    return -1;
2368
4.14k
2369
4.14k
  if (add_div_constraints(tab, pos - o_div, add_ineq, user) < 0)
2370
0
    return -1;
2371
4.14k
2372
4.14k
  return r;
2373
4.14k
}
2374
2375
/* Add an extra div, prescribed by "div", to the tableau and
2376
 * the associated bmap (which is assumed to be non-NULL).
2377
 */
2378
int isl_tab_add_div(struct isl_tab *tab, __isl_keep isl_vec *div)
2379
3.38k
{
2380
3.38k
  if (!tab)
2381
0
    return -1;
2382
3.38k
  return isl_tab_insert_div(tab, tab->n_var, div, NULL, NULL);
2383
3.38k
}
2384
2385
/* If "track" is set, then we want to keep track of all constraints in tab
2386
 * in its bmap field.  This field is initialized from a copy of "bmap",
2387
 * so we need to make sure that all constraints in "bmap" also appear
2388
 * in the constructed tab.
2389
 */
2390
__isl_give struct isl_tab *isl_tab_from_basic_map(
2391
  __isl_keep isl_basic_map *bmap, int track)
2392
1.21M
{
2393
1.21M
  int i;
2394
1.21M
  struct isl_tab *tab;
2395
1.21M
2396
1.21M
  if (!bmap)
2397
0
    return NULL;
2398
1.21M
  tab = isl_tab_alloc(bmap->ctx,
2399
1.21M
          isl_basic_map_total_dim(bmap) + bmap->n_ineq + 1,
2400
1.21M
          isl_basic_map_total_dim(bmap), 0);
2401
1.21M
  if (!tab)
2402
0
    return NULL;
2403
1.21M
  tab->preserve = track;
2404
1.21M
  tab->rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
2405
1.21M
  if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY)) {
2406
3
    if (isl_tab_mark_empty(tab) < 0)
2407
0
      goto error;
2408
3
    goto done;
2409
3
  }
2410
1.88M
  
for (i = 0; 1.21M
i < bmap->n_eq;
++i665k
) {
2411
665k
    tab = add_eq(tab, bmap->eq[i]);
2412
665k
    if (!tab)
2413
0
      return tab;
2414
665k
  }
2415
8.46M
  
for (i = 0; 1.21M
i < bmap->n_ineq;
++i7.24M
) {
2416
7.25M
    if (isl_tab_add_ineq(tab, bmap->ineq[i]) < 0)
2417
0
      goto error;
2418
7.25M
    if (tab->empty)
2419
6.30k
      goto done;
2420
7.25M
  }
2421
1.21M
done:
2422
1.21M
  if (track && 
isl_tab_track_bmap(tab, isl_basic_map_copy(bmap)) < 0325k
)
2423
0
    goto error;
2424
1.21M
  return tab;
2425
0
error:
2426
0
  isl_tab_free(tab);
2427
0
  return NULL;
2428
1.21M
}
2429
2430
__isl_give struct isl_tab *isl_tab_from_basic_set(
2431
  __isl_keep isl_basic_set *bset, int track)
2432
503k
{
2433
503k
  return isl_tab_from_basic_map(bset, track);
2434
503k
}
2435
2436
/* Construct a tableau corresponding to the recession cone of "bset".
2437
 */
2438
struct isl_tab *isl_tab_from_recession_cone(__isl_keep isl_basic_set *bset,
2439
  int parametric)
2440
3.76k
{
2441
3.76k
  isl_int cst;
2442
3.76k
  int i;
2443
3.76k
  struct isl_tab *tab;
2444
3.76k
  unsigned offset = 0;
2445
3.76k
2446
3.76k
  if (!bset)
2447
0
    return NULL;
2448
3.76k
  if (parametric)
2449
2.89k
    offset = isl_basic_set_dim(bset, isl_dim_param);
2450
3.76k
  tab = isl_tab_alloc(bset->ctx, bset->n_eq + bset->n_ineq,
2451
3.76k
        isl_basic_set_total_dim(bset) - offset, 0);
2452
3.76k
  if (!tab)
2453
0
    return NULL;
2454
3.76k
  tab->rational = ISL_F_ISSET(bset, ISL_BASIC_SET_RATIONAL);
2455
3.76k
  tab->cone = 1;
2456
3.76k
2457
3.76k
  isl_int_init(cst);
2458
3.76k
  isl_int_set_si(cst, 0);
2459
6.04k
  for (i = 0; i < bset->n_eq; 
++i2.28k
) {
2460
2.28k
    isl_int_swap(bset->eq[i][offset], cst);
2461
2.28k
    if (offset > 0) {
2462
579
      if (isl_tab_add_eq(tab, bset->eq[i] + offset) < 0)
2463
0
        goto error;
2464
1.70k
    } else
2465
1.70k
      tab = add_eq(tab, bset->eq[i]);
2466
2.28k
    isl_int_swap(bset->eq[i][offset], cst);
2467
2.28k
    if (!tab)
2468
0
      goto done;
2469
2.28k
  }
2470
15.8k
  
for (i = 0; 3.76k
i < bset->n_ineq;
++i12.0k
) {
2471
12.0k
    int r;
2472
12.0k
    isl_int_swap(bset->ineq[i][offset], cst);
2473
12.0k
    r = isl_tab_add_row(tab, bset->ineq[i] + offset);
2474
12.0k
    isl_int_swap(bset->ineq[i][offset], cst);
2475
12.0k
    if (r < 0)
2476
0
      goto error;
2477
12.0k
    tab->con[r].is_nonneg = 1;
2478
12.0k
    if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
2479
0
      goto error;
2480
12.0k
  }
2481
3.76k
done:
2482
3.76k
  isl_int_clear(cst);
2483
3.76k
  return tab;
2484
0
error:
2485
0
  isl_int_clear(cst);
2486
0
  isl_tab_free(tab);
2487
0
  return NULL;
2488
3.76k
}
2489
2490
/* Assuming "tab" is the tableau of a cone, check if the cone is
2491
 * bounded, i.e., if it is empty or only contains the origin.
2492
 */
2493
isl_bool isl_tab_cone_is_bounded(struct isl_tab *tab)
2494
2.89k
{
2495
2.89k
  int i;
2496
2.89k
2497
2.89k
  if (!tab)
2498
0
    return isl_bool_error;
2499
2.89k
  if (tab->empty)
2500
0
    return isl_bool_true;
2501
2.89k
  if (tab->n_dead == tab->n_col)
2502
741
    return isl_bool_true;
2503
2.15k
2504
3.37k
  
for (;;)2.15k
{
2505
3.59k
    for (i = tab->n_redundant; i < tab->n_row; 
++i221
) {
2506
3.59k
      struct isl_tab_var *var;
2507
3.59k
      int sgn;
2508
3.59k
      var = isl_tab_var_from_row(tab, i);
2509
3.59k
      if (!var->is_nonneg)
2510
221
        continue;
2511
3.36k
      sgn = sign_of_max(tab, var);
2512
3.36k
      if (sgn < -1)
2513
0
        return isl_bool_error;
2514
3.36k
      if (sgn != 0)
2515
280
        return isl_bool_false;
2516
3.08k
      if (close_row(tab, var, 0) < 0)
2517
0
        return isl_bool_error;
2518
3.08k
      break;
2519
3.08k
    }
2520
3.37k
    
if (3.09k
tab->n_dead == tab->n_col3.09k
)
2521
1.86k
      return isl_bool_true;
2522
1.22k
    if (i == tab->n_row)
2523
3
      return isl_bool_false;
2524
1.22k
  }
2525
2.15k
}
2526
2527
int isl_tab_sample_is_integer(struct isl_tab *tab)
2528
837k
{
2529
837k
  int i;
2530
837k
2531
837k
  if (!tab)
2532
0
    return -1;
2533
837k
2534
4.16M
  
for (i = 0; 837k
i < tab->n_var;
++i3.32M
) {
2535
3.59M
    int row;
2536
3.59M
    if (!tab->var[i].is_row)
2537
964k
      continue;
2538
2.63M
    row = tab->var[i].index;
2539
2.63M
    if (!isl_int_is_divisible_by(tab->mat->row[row][1],
2540
2.63M
            tab->mat->row[row][0]))
2541
2.63M
      
return 0271k
;
2542
2.63M
  }
2543
837k
  
return 1565k
;
2544
837k
}
2545
2546
static struct isl_vec *extract_integer_sample(struct isl_tab *tab)
2547
317k
{
2548
317k
  int i;
2549
317k
  struct isl_vec *vec;
2550
317k
2551
317k
  vec = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
2552
317k
  if (!vec)
2553
0
    return NULL;
2554
317k
2555
317k
  isl_int_set_si(vec->block.data[0], 1);
2556
2.21M
  for (i = 0; i < tab->n_var; 
++i1.89M
) {
2557
1.89M
    if (!tab->var[i].is_row)
2558
1.89M
      
isl_int_set_si662k
(vec->block.data[1 + i], 0);
2559
1.89M
    else {
2560
1.23M
      int row = tab->var[i].index;
2561
1.23M
      isl_int_divexact(vec->block.data[1 + i],
2562
1.23M
        tab->mat->row[row][1], tab->mat->row[row][0]);
2563
1.23M
    }
2564
1.89M
  }
2565
317k
2566
317k
  return vec;
2567
317k
}
2568
2569
struct isl_vec *isl_tab_get_sample_value(struct isl_tab *tab)
2570
459k
{
2571
459k
  int i;
2572
459k
  struct isl_vec *vec;
2573
459k
  isl_int m;
2574
459k
2575
459k
  if (!tab)
2576
0
    return NULL;
2577
459k
2578
459k
  vec = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
2579
459k
  if (!vec)
2580
0
    return NULL;
2581
459k
2582
459k
  isl_int_init(m);
2583
459k
2584
459k
  isl_int_set_si(vec->block.data[0], 1);
2585
2.14M
  for (i = 0; i < tab->n_var; 
++i1.68M
) {
2586
1.68M
    int row;
2587
1.68M
    if (!tab->var[i].is_row) {
2588
744k
      isl_int_set_si(vec->block.data[1 + i], 0);
2589
744k
      continue;
2590
744k
    }
2591
943k
    row = tab->var[i].index;
2592
943k
    isl_int_gcd(m, vec->block.data[0], tab->mat->row[row][0]);
2593
943k
    isl_int_divexact(m, tab->mat->row[row][0], m);
2594
943k
    isl_seq_scale(vec->block.data, vec->block.data, m, 1 + i);
2595
943k
    isl_int_divexact(m, vec->block.data[0], tab->mat->row[row][0]);
2596
943k
    isl_int_mul(vec->block.data[1 + i], m, tab->mat->row[row][1]);
2597
943k
  }
2598
459k
  vec = isl_vec_normalize(vec);
2599
459k
2600
459k
  isl_int_clear(m);
2601
459k
  return vec;
2602
459k
}
2603
2604
/* Store the sample value of "var" of "tab" rounded up (if sgn > 0)
2605
 * or down (if sgn < 0) to the nearest integer in *v.
2606
 */
2607
static void get_rounded_sample_value(struct isl_tab *tab,
2608
  struct isl_tab_var *var, int sgn, isl_int *v)
2609
661k
{
2610
661k
  if (!var->is_row)
2611
661k
    
isl_int_set_si6.22k
(*v, 0);
2612
661k
  else 
if (655k
sgn > 0655k
)
2613
655k
    
isl_int_cdiv_q650k
(*v, tab->mat->row[var->index][1],
2614
655k
           tab->mat->row[var->index][0]);
2615
655k
  else
2616
655k
    
isl_int_fdiv_q4.44k
(*v, tab->mat->row[var->index][1],
2617
661k
           tab->mat->row[var->index][0]);
2618
661k
}
2619
2620
/* Update "bmap" based on the results of the tableau "tab".
2621
 * In particular, implicit equalities are made explicit, redundant constraints
2622
 * are removed and if the sample value happens to be integer, it is stored
2623
 * in "bmap" (unless "bmap" already had an integer sample).
2624
 *
2625
 * The tableau is assumed to have been created from "bmap" using
2626
 * isl_tab_from_basic_map.
2627
 */
2628
struct isl_basic_map *isl_basic_map_update_from_tab(struct isl_basic_map *bmap,
2629
  struct isl_tab *tab)
2630
682k
{
2631
682k
  int i;
2632
682k
  unsigned n_eq;
2633
682k
2634
682k
  if (!bmap)
2635
0
    return NULL;
2636
682k
  if (!tab)
2637
0
    return bmap;
2638
682k
2639
682k
  n_eq = tab->n_eq;
2640
682k
  if (tab->empty)
2641
7.02k
    bmap = isl_basic_map_set_to_empty(bmap);
2642
675k
  else
2643
5.33M
    
for (i = bmap->n_ineq - 1; 675k
i >= 0;
--i4.66M
) {
2644
4.66M
      if (isl_tab_is_equality(tab, n_eq + i))
2645
2.16M
        isl_basic_map_inequality_to_equality(bmap, i);
2646
2.49M
      else if (isl_tab_is_redundant(tab, n_eq + i))
2647
277k
        isl_basic_map_drop_inequality(bmap, i);
2648
4.66M
    }
2649
682k
  if (bmap->n_eq != n_eq)
2650
299k
    bmap = isl_basic_map_gauss(bmap, NULL);
2651
682k
  if (!tab->rational &&
2652
682k
      
bmap648k
&&
!bmap->sample648k
&&
isl_tab_sample_is_integer(tab)355k
)
2653
317k
    bmap->sample = extract_integer_sample(tab);
2654
682k
  return bmap;
2655
682k
}
2656
2657
struct isl_basic_set *isl_basic_set_update_from_tab(struct isl_basic_set *bset,
2658
  struct isl_tab *tab)
2659
47.3k
{
2660
47.3k
  return bset_from_bmap(isl_basic_map_update_from_tab(bset_to_bmap(bset),
2661
47.3k
                tab));
2662
47.3k
}
2663
2664
/* Drop the last constraint added to "tab" in position "r".
2665
 * The constraint is expected to have remained in a row.
2666
 */
2667
static isl_stat drop_last_con_in_row(struct isl_tab *tab, int r)
2668
19.0k
{
2669
19.0k
  if (!tab->con[r].is_row)
2670
19.0k
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_internal,
2671
19.0k
      "row unexpectedly moved to column",
2672
19.0k
      return isl_stat_error);
2673
19.0k
  if (r + 1 != tab->n_con)
2674
19.0k
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_internal,
2675
19.0k
      "additional constraints added", return isl_stat_error);
2676
19.0k
  if (drop_row(tab, tab->con[r].index) < 0)
2677
0
    return isl_stat_error;
2678
19.0k
2679
19.0k
  return isl_stat_ok;
2680
19.0k
}
2681
2682
/* Given a non-negative variable "var", temporarily add a new non-negative
2683
 * variable that is the opposite of "var", ensuring that "var" can only attain
2684
 * the value zero.  The new variable is removed again before this function
2685
 * returns.  However, the effect of forcing "var" to be zero remains.
2686
 * If var = n/d is a row variable, then the new variable = -n/d.
2687
 * If var is a column variables, then the new variable = -var.
2688
 * If the new variable cannot attain non-negative values, then
2689
 * the resulting tableau is empty.
2690
 * Otherwise, we know the value will be zero and we close the row.
2691
 */
2692
static isl_stat cut_to_hyperplane(struct isl_tab *tab, struct isl_tab_var *var)
2693
19.0k
{
2694
19.0k
  unsigned r;
2695
19.0k
  isl_int *row;
2696
19.0k
  int sgn;
2697
19.0k
  unsigned off = 2 + tab->M;
2698
19.0k
2699
19.0k
  if (var->is_zero)
2700
0
    return isl_stat_ok;
2701
19.0k
  if (var->is_redundant || !var->is_nonneg)
2702
19.0k
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_invalid,
2703
19.0k
      "expecting non-redundant non-negative variable",
2704
19.0k
      return isl_stat_error);
2705
19.0k
2706
19.0k
  if (isl_tab_extend_cons(tab, 1) < 0)
2707
0
    return isl_stat_error;
2708
19.0k
2709
19.0k
  r = tab->n_con;
2710
19.0k
  tab->con[r].index = tab->n_row;
2711
19.0k
  tab->con[r].is_row = 1;
2712
19.0k
  tab->con[r].is_nonneg = 0;
2713
19.0k
  tab->con[r].is_zero = 0;
2714
19.0k
  tab->con[r].is_redundant = 0;
2715
19.0k
  tab->con[r].frozen = 0;
2716
19.0k
  tab->con[r].negated = 0;
2717
19.0k
  tab->row_var[tab->n_row] = ~r;
2718
19.0k
  row = tab->mat->row[tab->n_row];
2719
19.0k
2720
19.0k
  if (var->is_row) {
2721
8.08k
    isl_int_set(row[0], tab->mat->row[var->index][0]);
2722
8.08k
    isl_seq_neg(row + 1,
2723
8.08k
          tab->mat->row[var->index] + 1, 1 + tab->n_col);
2724
11.0k
  } else {
2725
11.0k
    isl_int_set_si(row[0], 1);
2726
11.0k
    isl_seq_clr(row + 1, 1 + tab->n_col);
2727
11.0k
    isl_int_set_si(row[off + var->index], -1);
2728
11.0k
  }
2729
19.0k
2730
19.0k
  tab->n_row++;
2731
19.0k
  tab->n_con++;
2732
19.0k
2733
19.0k
  sgn = sign_of_max(tab, &tab->con[r]);
2734
19.0k
  if (sgn < -1)
2735
0
    return isl_stat_error;
2736
19.0k
  if (sgn < 0) {
2737
47
    if (drop_last_con_in_row(tab, r) < 0)
2738
0
      return isl_stat_error;
2739
47
    if (isl_tab_mark_empty(tab) < 0)
2740
0
      return isl_stat_error;
2741
47
    return isl_stat_ok;
2742
47
  }
2743
19.0k
  tab->con[r].is_nonneg = 1;
2744
19.0k
  /* sgn == 0 */
2745
19.0k
  if (close_row(tab, &tab->con[r], 1) < 0)
2746
0
    return isl_stat_error;
2747
19.0k
  if (drop_last_con_in_row(tab, r) < 0)
2748
0
    return isl_stat_error;
2749
19.0k
2750
19.0k
  return isl_stat_ok;
2751
19.0k
}
2752
2753
/* Given a tableau "tab" and an inequality constraint "con" of the tableau,
2754
 * relax the inequality by one.  That is, the inequality r >= 0 is replaced
2755
 * by r' = r + 1 >= 0.
2756
 * If r is a row variable, we simply increase the constant term by one
2757
 * (taking into account the denominator).
2758
 * If r is a column variable, then we need to modify each row that
2759
 * refers to r = r' - 1 by substituting this equality, effectively
2760
 * subtracting the coefficient of the column from the constant.
2761
 * We should only do this if the minimum is manifestly unbounded,
2762
 * however.  Otherwise, we may end up with negative sample values
2763
 * for non-negative variables.
2764
 * So, if r is a column variable with a minimum that is not
2765
 * manifestly unbounded, then we need to move it to a row.
2766
 * However, the sample value of this row may be negative,
2767
 * even after the relaxation, so we need to restore it.
2768
 * We therefore prefer to pivot a column up to a row, if possible.
2769
 */
2770
int isl_tab_relax(struct isl_tab *tab, int con)
2771
8.44k
{
2772
8.44k
  struct isl_tab_var *var;
2773
8.44k
2774
8.44k
  if (!tab)
2775
0
    return -1;
2776
8.44k
2777
8.44k
  var = &tab->con[con];
2778
8.44k
2779
8.44k
  if (var->is_row && 
(55
var->index < 055
||
var->index < tab->n_redundant55
))
2780
8.44k
    
isl_die0
(tab->mat->ctx, isl_error_invalid,
2781
8.44k
      "cannot relax redundant constraint", return -1);
2782
8.44k
  if (!var->is_row && 
(8.38k
var->index < 08.38k
||
var->index < tab->n_dead8.38k
))
2783
8.44k
    
isl_die0
(tab->mat->ctx, isl_error_invalid,
2784
8.44k
      "cannot relax dead constraint", return -1);
2785
8.44k
2786
8.44k
  if (!var->is_row && 
!max_is_manifestly_unbounded(tab, var)8.38k
)
2787
6.92k
    if (to_row(tab, var, 1) < 0)
2788
0
      return -1;
2789
8.44k
  if (!var->is_row && 
!min_is_manifestly_unbounded(tab, var)1.46k
)
2790
18
    if (to_row(tab, var, -1) < 0)
2791
0
      return -1;
2792
8.44k
2793
8.44k
  if (var->is_row) {
2794
6.99k
    isl_int_add(tab->mat->row[var->index][1],
2795
6.99k
        tab->mat->row[var->index][1], tab->mat->row[var->index][0]);
2796
6.99k
    if (restore_row(tab, var) < 0)
2797
0
      return -1;
2798
1.44k
  } else {
2799
1.44k
    int i;
2800
1.44k
    unsigned off = 2 + tab->M;
2801
1.44k
2802
9.30k
    for (i = 0; i < tab->n_row; 
++i7.86k
) {
2803
7.86k
      if (isl_int_is_zero(tab->mat->row[i][off + var->index]))
2804
7.86k
        
continue6.30k
;
2805
1.55k
      isl_int_sub(tab->mat->row[i][1], tab->mat->row[i][1],
2806
1.55k
          tab->mat->row[i][off + var->index]);
2807
1.55k
    }
2808
1.44k
2809
1.44k
  }
2810
8.44k
2811
8.44k
  if (isl_tab_push_var(tab, isl_tab_undo_relax, var) < 0)
2812
0
    return -1;
2813
8.44k
2814
8.44k
  return 0;
2815
8.44k
}
2816
2817
/* Replace the variable v at position "pos" in the tableau "tab"
2818
 * by v' = v + shift.
2819
 *
2820
 * If the variable is in a column, then we first check if we can
2821
 * simply plug in v = v' - shift.  The effect on a row with
2822
 * coefficient f/d for variable v is that the constant term c/d
2823
 * is replaced by (c - f * shift)/d.  If shift is positive and
2824
 * f is negative for each row that needs to remain non-negative,
2825
 * then this is clearly safe.  In other words, if the minimum of v
2826
 * is manifestly unbounded, then we can keep v in a column position.
2827
 * Otherwise, we can pivot it down to a row.
2828
 * Similarly, if shift is negative, we need to check if the maximum
2829
 * of is manifestly unbounded.
2830
 *
2831
 * If the variable is in a row (from the start or after pivoting),
2832
 * then the constant term c/d is replaced by (c + d * shift)/d.
2833
 */
2834
int isl_tab_shift_var(struct isl_tab *tab, int pos, isl_int shift)
2835
176
{
2836
176
  struct isl_tab_var *var;
2837
176
2838
176
  if (!tab)
2839
0
    return -1;
2840
176
  if (isl_int_is_zero(shift))
2841
176
    
return 091
;
2842
85
2843
85
  var = &tab->var[pos];
2844
85
  if (!var->is_row) {
2845
15
    if (isl_int_is_neg(shift)) {
2846
10
      if (!max_is_manifestly_unbounded(tab, var))
2847
7
        if (to_row(tab, var, 1) < 0)
2848
0
          return -1;
2849
5
    } else {
2850
5
      if (!min_is_manifestly_unbounded(tab, var))
2851
1
        if (to_row(tab, var, -1) < 0)
2852
0
          return -1;
2853
85
    }
2854
15
  }
2855
85
2856
85
  if (var->is_row) {
2857
78
    isl_int_addmul(tab->mat->row[var->index][1],
2858
78
        shift, tab->mat->row[var->index][0]);
2859
78
  } else {
2860
7
    int i;
2861
7
    unsigned off = 2 + tab->M;
2862
7
2863
30
    for (i = 0; i < tab->n_row; 
++i23
) {
2864
23
      if (isl_int_is_zero(tab->mat->row[i][off + var->index]))
2865
23
        
continue15
;
2866
8
      isl_int_submul(tab->mat->row[i][1],
2867
8
            shift, tab->mat->row[i][off + var->index]);
2868
8
    }
2869
7
2870
7
  }
2871
85
2872
85
  return 0;
2873
85
}
2874
2875
/* Remove the sign constraint from constraint "con".
2876
 *
2877
 * If the constraint variable was originally marked non-negative,
2878
 * then we make sure we mark it non-negative again during rollback.
2879
 */
2880
int isl_tab_unrestrict(struct isl_tab *tab, int con)
2881
977
{
2882
977
  struct isl_tab_var *var;
2883
977
2884
977
  if (!tab)
2885
0
    return -1;
2886
977
2887
977
  var = &tab->con[con];
2888
977
  if (!var->is_nonneg)
2889
0
    return 0;
2890
977
2891
977
  var->is_nonneg = 0;
2892
977
  if (isl_tab_push_var(tab, isl_tab_undo_unrestrict, var) < 0)
2893
0
    return -1;
2894
977
2895
977
  return 0;
2896
977
}
2897
2898
int isl_tab_select_facet(struct isl_tab *tab, int con)
2899
18.2k
{
2900
18.2k
  if (!tab)
2901
0
    return -1;
2902
18.2k
2903
18.2k
  return cut_to_hyperplane(tab, &tab->con[con]);
2904
18.2k
}
2905
2906
static int may_be_equality(struct isl_tab *tab, int row)
2907
12.3M
{
2908
12.3M
  return tab->rational ? 
isl_int_is_zero36.3k
(tab->mat->row[row][1])
2909
12.3M
           : 
isl_int_lt12.3M
(tab->mat->row[row][1],
2910
12.3M
              tab->mat->row[row][0]);
2911
12.3M
}
2912
2913
/* Return an isl_tab_var that has been marked or NULL if no such
2914
 * variable can be found.
2915
 * The marked field has only been set for variables that
2916
 * appear in non-redundant rows or non-dead columns.
2917
 *
2918
 * Pick the last constraint variable that is marked and
2919
 * that appears in either a non-redundant row or a non-dead columns.
2920
 * Since the returned variable is tested for being a redundant constraint or
2921
 * an implicit equality, there is no need to return any tab variable that
2922
 * corresponds to a variable.
2923
 */
2924
static struct isl_tab_var *select_marked(struct isl_tab *tab)
2925
3.79M
{
2926
3.79M
  int i;
2927
3.79M
  struct isl_tab_var *var;
2928
3.79M
2929
26.4M
  for (i = tab->n_con - 1; i >= 0; 
--i22.6M
) {
2930
26.0M
    var = &tab->con[i];
2931
26.0M
    if (var->index < 0)
2932
9.59M
      continue;
2933
16.4M
    if (var->is_row && 
var->index < tab->n_redundant11.9M
)
2934
868k
      continue;
2935
15.6M
    if (!var->is_row && 
var->index < tab->n_dead4.52M
)
2936
2.49k
      continue;
2937
15.6M
    if (var->marked)
2938
3.48M
      return var;
2939
15.6M
  }
2940
3.79M
2941
3.79M
  
return NULL309k
;
2942
3.79M
}
2943
2944
/* Check for (near) equalities among the constraints.
2945
 * A constraint is an equality if it is non-negative and if
2946
 * its maximal value is either
2947
 *  - zero (in case of rational tableaus), or
2948
 *  - strictly less than 1 (in case of integer tableaus)
2949
 *
2950
 * We first mark all non-redundant and non-dead variables that
2951
 * are not frozen and not obviously not an equality.
2952
 * Then we iterate over all marked variables if they can attain
2953
 * any values larger than zero or at least one.
2954
 * If the maximal value is zero, we mark any column variables
2955
 * that appear in the row as being zero and mark the row as being redundant.
2956
 * Otherwise, if the maximal value is strictly less than one (and the
2957
 * tableau is integer), then we restrict the value to being zero
2958
 * by adding an opposite non-negative variable.
2959
 * The order in which the variables are considered is not important.
2960
 */
2961
int isl_tab_detect_implicit_equalities(struct isl_tab *tab)
2962
705k
{
2963
705k
  int i;
2964
705k
  unsigned n_marked;
2965
705k
2966
705k
  if (!tab)
2967
0
    return -1;
2968
705k
  if (tab->empty)
2969
4.14k
    return 0;
2970
701k
  if (tab->n_dead == tab->n_col)
2971
18.2k
    return 0;
2972
682k
2973
682k
  n_marked = 0;
2974
5.56M
  for (i = tab->n_redundant; i < tab->n_row; 
++i4.88M
) {
2975
4.88M
    struct isl_tab_var *var = isl_tab_var_from_row(tab, i);
2976
4.88M
    var->marked = !var->frozen && 
var->is_nonneg4.85M
&&
2977
4.88M
      
may_be_equality(tab, i)4.32M
;
2978
4.88M
    if (var->marked)
2979
3.22M
      n_marked++;
2980
4.88M
  }
2981
4.31M
  for (i = tab->n_dead; i < tab->n_col; 
++i3.63M
) {
2982
3.63M
    struct isl_tab_var *var = var_from_col(tab, i);
2983
3.63M
    var->marked = !var->frozen && 
var->is_nonneg3.62M
;
2984
3.63M
    if (var->marked)
2985
370k
      n_marked++;
2986
3.63M
  }
2987
3.06M
  while (n_marked) {
2988
2.68M
    struct isl_tab_var *var;
2989
2.68M
    int sgn;
2990
2.68M
    var = select_marked(tab);
2991
2.68M
    if (!var)
2992
297k
      break;
2993
2.38M
    var->marked = 0;
2994
2.38M
    n_marked--;
2995
2.38M
    sgn = sign_of_max(tab, var);
2996
2.38M
    if (sgn < 0)
2997
0
      return -1;
2998
2.38M
    if (sgn == 0) {
2999
1.10M
      if (close_row(tab, var, 0) < 0)
3000
0
        return -1;
3001
1.27M
    } else if (!tab->rational && 
!at_least_one(tab, var)1.25M
) {
3002
856
      if (cut_to_hyperplane(tab, var) < 0)
3003
0
        return -1;
3004
856
      return isl_tab_detect_implicit_equalities(tab);
3005
856
    }
3006
19.1M
    
for (i = tab->n_redundant; 2.38M
i < tab->n_row;
++i16.7M
) {
3007
16.7M
      var = isl_tab_var_from_row(tab, i);
3008
16.7M
      if (!var->marked)
3009
8.72M
        continue;
3010
8.03M
      if (may_be_equality(tab, i))
3011
7.95M
        continue;
3012
84.1k
      var->marked = 0;
3013
84.1k
      n_marked--;
3014
84.1k
    }
3015
2.38M
  }
3016
682k
3017
682k
  
return 0682k
;
3018
682k
}
3019
3020
/* Update the element of row_var or col_var that corresponds to
3021
 * constraint tab->con[i] to a move from position "old" to position "i".
3022
 */
3023
static int update_con_after_move(struct isl_tab *tab, int i, int old)
3024
6.29k
{
3025
6.29k
  int *p;
3026
6.29k
  int index;
3027
6.29k
3028
6.29k
  index = tab->con[i].index;
3029
6.29k
  if (index == -1)
3030
4.15k
    return 0;
3031
2.14k
  p = tab->con[i].is_row ? 
tab->row_var1.47k
:
tab->col_var667
;
3032
2.14k
  if (p[index] != ~old)
3033
2.14k
    
isl_die0
(tab->mat->ctx, isl_error_internal,
3034
2.14k
      "broken internal state", return -1);
3035
2.14k
  p[index] = ~i;
3036
2.14k
3037
2.14k
  return 0;
3038
2.14k
}
3039
3040
/* Rotate the "n" constraints starting at "first" to the right,
3041
 * putting the last constraint in the position of the first constraint.
3042
 */
3043
static int rotate_constraints(struct isl_tab *tab, int first, int n)
3044
1.80k
{
3045
1.80k
  int i, last;
3046
1.80k
  struct isl_tab_var var;
3047
1.80k
3048
1.80k
  if (n <= 1)
3049
578
    return 0;
3050
1.22k
3051
1.22k
  last = first + n - 1;
3052
1.22k
  var = tab->con[last];
3053
6.29k
  for (i = last; i > first; 
--i5.07k
) {
3054
5.07k
    tab->con[i] = tab->con[i - 1];
3055
5.07k
    if (update_con_after_move(tab, i, i - 1) < 0)
3056
0
      return -1;
3057
5.07k
  }
3058
1.22k
  tab->con[first] = var;
3059
1.22k
  if (update_con_after_move(tab, first, last) < 0)
3060
0
    return -1;
3061
1.22k
3062
1.22k
  return 0;
3063
1.22k
}
3064
3065
/* Make the equalities that are implicit in "bmap" but that have been
3066
 * detected in the corresponding "tab" explicit in "bmap" and update
3067
 * "tab" to reflect the new order of the constraints.
3068
 *
3069
 * In particular, if inequality i is an implicit equality then
3070
 * isl_basic_map_inequality_to_equality will move the inequality
3071
 * in front of the other equality and it will move the last inequality
3072
 * in the position of inequality i.
3073
 * In the tableau, the inequalities of "bmap" are stored after the equalities
3074
 * and so the original order
3075
 *
3076
 *    E E E E E A A A I B B B B L
3077
 *
3078
 * is changed into
3079
 *
3080
 *    I E E E E E A A A L B B B B
3081
 *
3082
 * where I is the implicit equality, the E are equalities,
3083
 * the A inequalities before I, the B inequalities after I and
3084
 * L the last inequality.
3085
 * We therefore need to rotate to the right two sets of constraints,
3086
 * those up to and including I and those after I.
3087
 *
3088
 * If "tab" contains any constraints that are not in "bmap" then they
3089
 * appear after those in "bmap" and they should be left untouched.
3090
 *
3091
 * Note that this function leaves "bmap" in a temporary state
3092
 * as it does not call isl_basic_map_gauss.  Calling this function
3093
 * is the responsibility of the caller.
3094
 */
3095
__isl_give isl_basic_map *isl_tab_make_equalities_explicit(struct isl_tab *tab,
3096
  __isl_take isl_basic_map *bmap)
3097
67.6k
{
3098
67.6k
  int i;
3099
67.6k
3100
67.6k
  if (!tab || !bmap)
3101
0
    return isl_basic_map_free(bmap);
3102
67.6k
  if (tab->empty)
3103
80
    return bmap;
3104
67.6k
3105
271k
  
for (i = bmap->n_ineq - 1; 67.6k
i >= 0;
--i204k
) {
3106
204k
    if (!isl_tab_is_equality(tab, bmap->n_eq + i))
3107
203k
      continue;
3108
900
    isl_basic_map_inequality_to_equality(bmap, i);
3109
900
    if (rotate_constraints(tab, 0, tab->n_eq + i + 1) < 0)
3110
0
      return isl_basic_map_free(bmap);
3111
900
    if (rotate_constraints(tab, tab->n_eq + i + 1,
3112
900
          bmap->n_ineq - i) < 0)
3113
0
      return isl_basic_map_free(bmap);
3114
900
    tab->n_eq++;
3115
900
  }
3116
67.6k
3117
67.6k
  return bmap;
3118
67.6k
}
3119
3120
static int con_is_redundant(struct isl_tab *tab, struct isl_tab_var *var)
3121
1.70M
{
3122
1.70M
  if (!tab)
3123
0
    return -1;
3124
1.70M
  if (tab->rational) {
3125
70.3k
    int sgn = sign_of_min(tab, var);
3126
70.3k
    if (sgn < -1)
3127
0
      return -1;
3128
70.3k
    return sgn >= 0;
3129
1.63M
  } else {
3130
1.63M
    int irred = isl_tab_min_at_most_neg_one(tab, var);
3131
1.63M
    if (irred < 0)
3132
0
      return -1;
3133
1.63M
    return !irred;
3134
1.63M
  }
3135
1.70M
}
3136
3137
/* Check for (near) redundant constraints.
3138
 * A constraint is redundant if it is non-negative and if
3139
 * its minimal value (temporarily ignoring the non-negativity) is either
3140
 *  - zero (in case of rational tableaus), or
3141
 *  - strictly larger than -1 (in case of integer tableaus)
3142
 *
3143
 * We first mark all non-redundant and non-dead variables that
3144
 * are not frozen and not obviously negatively unbounded.
3145
 * Then we iterate over all marked variables if they can attain
3146
 * any values smaller than zero or at most negative one.
3147
 * If not, we mark the row as being redundant (assuming it hasn't
3148
 * been detected as being obviously redundant in the mean time).
3149
 */
3150
int isl_tab_detect_redundant(struct isl_tab *tab)
3151
316k
{
3152
316k
  int i;
3153
316k
  unsigned n_marked;
3154
316k
3155
316k
  if (!tab)
3156
0
    return -1;
3157
316k
  if (tab->empty)
3158
3.14k
    return 0;
3159
313k
  if (tab->n_redundant == tab->n_row)
3160
5.42k
    return 0;
3161
307k
3162
307k
  n_marked = 0;
3163
2.75M
  for (i = tab->n_redundant; i < tab->n_row; 
++i2.44M
) {
3164
2.44M
    struct isl_tab_var *var = isl_tab_var_from_row(tab, i);
3165
2.44M
    var->marked = !var->frozen && 
var->is_nonneg2.25M
;
3166
2.44M
    if (var->marked)
3167
1.06M
      n_marked++;
3168
2.44M
  }
3169
1.94M
  for (i = tab->n_dead; i < tab->n_col; 
++i1.64M
) {
3170
1.64M
    struct isl_tab_var *var = var_from_col(tab, i);
3171
1.64M
    var->marked = !var->frozen && 
var->is_nonneg1.57M
&&
3172
1.64M
      
!min_is_manifestly_unbounded(tab, var)798k
;
3173
1.64M
    if (var->marked)
3174
245k
      n_marked++;
3175
1.64M
  }
3176
1.41M
  while (n_marked) {
3177
1.11M
    struct isl_tab_var *var;
3178
1.11M
    int red;
3179
1.11M
    var = select_marked(tab);
3180
1.11M
    if (!var)
3181
11.6k
      break;
3182
1.10M
    var->marked = 0;
3183
1.10M
    n_marked--;
3184
1.10M
    red = con_is_redundant(tab, var);
3185
1.10M
    if (red < 0)
3186
0
      return -1;
3187
1.10M
    if (red && 
!var->is_redundant237k
)
3188
52.8k
      if (isl_tab_mark_redundant(tab, var->index) < 0)
3189
0
        return -1;
3190
12.7M
    
for (i = tab->n_dead; 1.10M
i < tab->n_col;
++i11.6M
) {
3191
11.6M
      var = var_from_col(tab, i);
3192
11.6M
      if (!var->marked)
3193
10.9M
        continue;
3194
663k
      if (!min_is_manifestly_unbounded(tab, var))
3195
482k
        continue;
3196
181k
      var->marked = 0;
3197
181k
      n_marked--;
3198
181k
    }
3199
1.10M
  }
3200
307k
3201
307k
  return 0;
3202
307k
}
3203
3204
int isl_tab_is_equality(struct isl_tab *tab, int con)
3205
4.90M
{
3206
4.90M
  int row;
3207
4.90M
  unsigned off;
3208
4.90M
3209
4.90M
  if (!tab)
3210
0
    return -1;
3211
4.90M
  if (tab->con[con].is_zero)
3212
2.16M
    return 1;
3213
2.73M
  if (tab->con[con].is_redundant)
3214
286k
    return 0;
3215
2.44M
  if (!tab->con[con].is_row)
3216
1.36M
    return tab->con[con].index < tab->n_dead;
3217
1.08M
3218
1.08M
  row = tab->con[con].index;
3219
1.08M
3220
1.08M
  off = 2 + tab->M;
3221
1.08M
  return isl_int_is_zero(tab->mat->row[row][1]) &&
3222
1.08M
    
!row_is_big(tab, row)65.5k
&&
3223
1.08M
    isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
3224
65.5k
          tab->n_col - tab->n_dead) == -1;
3225
1.08M
}
3226
3227
/* Return the minimal value of the affine expression "f" with denominator
3228
 * "denom" in *opt, *opt_denom, assuming the tableau is not empty and
3229
 * the expression cannot attain arbitrarily small values.
3230
 * If opt_denom is NULL, then *opt is rounded up to the nearest integer.
3231
 * The return value reflects the nature of the result (empty, unbounded,
3232
 * minimal value returned in *opt).
3233
 *
3234
 * This function assumes that at least one more row and at least
3235
 * one more element in the constraint array are available in the tableau.
3236
 */
3237
enum isl_lp_result isl_tab_min(struct isl_tab *tab,
3238
  isl_int *f, isl_int denom, isl_int *opt, isl_int *opt_denom,
3239
  unsigned flags)
3240
788k
{
3241
788k
  int r;
3242
788k
  enum isl_lp_result res = isl_lp_ok;
3243
788k
  struct isl_tab_var *var;
3244
788k
  struct isl_tab_undo *snap;
3245
788k
3246
788k
  if (!tab)
3247
0
    return isl_lp_error;
3248
788k
3249
788k
  if (tab->empty)
3250
148
    return isl_lp_empty;
3251
788k
3252
788k
  snap = isl_tab_snap(tab);
3253
788k
  r = isl_tab_add_row(tab, f);
3254
788k
  if (r < 0)
3255
0
    return isl_lp_error;
3256
788k
  var = &tab->con[r];
3257
1.59M
  for (;;) {
3258
1.59M
    int row, col;
3259
1.59M
    find_pivot(tab, var, var, -1, &row, &col);
3260
1.59M
    if (row == var->index) {
3261
11.1k
      res = isl_lp_unbounded;
3262
11.1k
      break;
3263
11.1k
    }
3264
1.58M
    if (row == -1)
3265
777k
      break;
3266
807k
    if (isl_tab_pivot(tab, row, col) < 0)
3267
0
      return isl_lp_error;
3268
807k
  }
3269
788k
  isl_int_mul(tab->mat->row[var->index][0],
3270
788k
        tab->mat->row[var->index][0], denom);
3271
788k
  if (ISL_FL_ISSET(flags, ISL_TAB_SAVE_DUAL)) {
3272
84.2k
    int i;
3273
84.2k
3274
84.2k
    isl_vec_free(tab->dual);
3275
84.2k
    tab->dual = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_con);
3276
84.2k
    if (!tab->dual)
3277
0
      return isl_lp_error;
3278
84.2k
    isl_int_set(tab->dual->el[0], tab->mat->row[var->index][0]);
3279
3.15M
    for (i = 0; i < tab->n_con; 
++i3.07M
) {
3280
3.07M
      int pos;
3281
3.07M
      if (tab->con[i].is_row) {
3282
1.81M
        isl_int_set_si(tab->dual->el[1 + i], 0);
3283
1.81M
        continue;
3284
1.81M
      }
3285
1.26M
      pos = 2 + tab->M + tab->con[i].index;
3286
1.26M
      if (tab->con[i].negated)
3287
1.26M
        
isl_int_neg244k
(tab->dual->el[1 + i],
3288
1.26M
              tab->mat->row[var->index][pos]);
3289
1.26M
      else
3290
1.26M
        
isl_int_set1.01M
(tab->dual->el[1 + i],
3291
1.26M
              tab->mat->row[var->index][pos]);
3292
1.26M
    }
3293
84.2k
  }
3294
788k
  if (opt && 
res == isl_lp_ok788k
) {
3295
777k
    if (opt_denom) {
3296
129k
      isl_int_set(*opt, tab->mat->row[var->index][1]);
3297
129k
      isl_int_set(*opt_denom, tab->mat->row[var->index][0]);
3298
129k
    } else
3299
648k
      get_rounded_sample_value(tab, var, 1, opt);
3300
777k
  }
3301
788k
  if (isl_tab_rollback(tab, snap) < 0)
3302
0
    return isl_lp_error;
3303
788k
  return res;
3304
788k
}
3305
3306
/* Is the constraint at position "con" marked as being redundant?
3307
 * If it is marked as representing an equality, then it is not
3308
 * considered to be redundant.
3309
 * Note that isl_tab_mark_redundant marks both the isl_tab_var as
3310
 * redundant and moves the corresponding row into the first
3311
 * tab->n_redundant positions (or removes the row, assigning it index -1),
3312
 * so the final test is actually redundant itself.
3313
 */
3314
int isl_tab_is_redundant(struct isl_tab *tab, int con)
3315
3.20M
{
3316
3.20M
  if (!tab)
3317
0
    return -1;
3318
3.20M
  if (con < 0 || con >= tab->n_con)
3319
3.20M
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_invalid,
3320
3.20M
      "position out of bounds", return -1);
3321
3.20M
  if (tab->con[con].is_zero)
3322
160
    return 0;
3323
3.20M
  if (tab->con[con].is_redundant)
3324
500k
    return 1;
3325
2.69M
  return tab->con[con].is_row && 
tab->con[con].index < tab->n_redundant1.15M
;
3326
2.69M
}
3327
3328
/* Is variable "var" of "tab" fixed to a constant value by its row
3329
 * in the tableau?
3330
 * If so and if "value" is not NULL, then store this constant value
3331
 * in "value".
3332
 *
3333
 * That is, is it a row variable that only has non-zero coefficients
3334
 * for dead columns?
3335
 */
3336
static isl_bool is_constant(struct isl_tab *tab, struct isl_tab_var *var,
3337
  isl_int *value)
3338
10.5k
{
3339
10.5k
  unsigned off = 2 + tab->M;
3340
10.5k
  isl_mat *mat = tab->mat;
3341
10.5k
  int n;
3342
10.5k
  int row;
3343
10.5k
  int pos;
3344
10.5k
3345
10.5k
  if (!var->is_row)
3346
6.22k
    return isl_bool_false;
3347
4.36k
  row = var->index;
3348
4.36k
  if (row_is_big(tab, row))
3349
0
    return isl_bool_false;
3350
4.36k
  n = tab->n_col - tab->n_dead;
3351
4.36k
  pos = isl_seq_first_non_zero(mat->row[row] + off + tab->n_dead, n);
3352
4.36k
  if (pos != -1)
3353
4.14k
    return isl_bool_false;
3354
219
  if (value)
3355
219
    
isl_int_divexact0
(*value, mat->row[row][1], mat->row[row][0]);
3356
219
  return isl_bool_true;
3357
219
}
3358
3359
/* Has the variable "var' of "tab" reached a value that is greater than
3360
 * or equal (if sgn > 0) or smaller than or equal (if sgn < 0) to "target"?
3361
 * "tmp" has been initialized by the caller and can be used
3362
 * to perform local computations.
3363
 *
3364
 * If the sample value involves the big parameter, then any value
3365
 * is reached.
3366
 * Otherwise check if n/d >= t, i.e., n >= d * t (if sgn > 0)
3367
 * or n/d <= t, i.e., n <= d * t (if sgn < 0).
3368
 */
3369
static int reached(struct isl_tab *tab, struct isl_tab_var *var, int sgn,
3370
  isl_int target, isl_int *tmp)
3371
12.3k
{
3372
12.3k
  if (row_is_big(tab, var->index))
3373
0
    return 1;
3374
12.3k
  isl_int_mul(*tmp, tab->mat->row[var->index][0], target);
3375
12.3k
  if (sgn > 0)
3376
4.93k
    return isl_int_ge(tab->mat->row[var->index][1], *tmp);
3377
12.3k
  else
3378
12.3k
    
return 7.45k
isl_int_le7.45k
(tab->mat->row[var->index][1], *tmp);
3379
12.3k
}
3380
3381
/* Can variable "var" of "tab" attain the value "target" by
3382
 * pivoting up (if sgn > 0) or down (if sgn < 0)?
3383
 * If not, then pivot up [down] to the greatest [smallest]
3384
 * rational value.
3385
 * "tmp" has been initialized by the caller and can be used
3386
 * to perform local computations.
3387
 *
3388
 * If the variable is manifestly unbounded in the desired direction,
3389
 * then it can attain any value.
3390
 * Otherwise, it can be moved to a row.
3391
 * Continue pivoting until the target is reached.
3392
 * If no more pivoting can be performed, the maximal [minimal]
3393
 * rational value has been reached and the target cannot be reached.
3394
 * If the variable would be pivoted into a manifestly unbounded column,
3395
 * then the target can be reached.
3396
 */
3397
static isl_bool var_reaches(struct isl_tab *tab, struct isl_tab_var *var,
3398
  int sgn, isl_int target, isl_int *tmp)
3399
12.7k
{
3400
12.7k
  int row, col;
3401
12.7k
3402
12.7k
  if (sgn < 0 && 
min_is_manifestly_unbounded(tab, var)10.3k
)
3403
5.27k
    return isl_bool_true;
3404
7.52k
  if (sgn > 0 && 
max_is_manifestly_unbounded(tab, var)2.41k
)
3405
0
    return isl_bool_true;
3406
7.52k
  if (to_row(tab, var, sgn) < 0)
3407
0
    return isl_bool_error;
3408
12.3k
  
while (7.52k
!reached(tab, var, sgn, target, tmp)) {
3409
10.0k
    find_pivot(tab, var, var, sgn, &row, &col);
3410
10.0k
    if (row == -1)
3411
2.71k
      return isl_bool_false;
3412
7.31k
    if (row == var->index)
3413
2.44k
      return isl_bool_true;
3414
4.87k
    if (isl_tab_pivot(tab, row, col) < 0)
3415
0
      return isl_bool_error;
3416
4.87k
  }
3417
7.52k
3418
7.52k
  
return isl_bool_true2.36k
;
3419
7.52k
}
3420
3421
/* Check if variable "var" of "tab" can only attain a single (integer)
3422
 * value, and, if so, add an equality constraint to fix the variable
3423
 * to this single value and store the result in "target".
3424
 * "target" and "tmp" have been initialized by the caller.
3425
 *
3426
 * Given the current sample value, round it down and check
3427
 * whether it is possible to attain a strictly smaller integer value.
3428
 * If so, the variable is not restricted to a single integer value.
3429
 * Otherwise, the search stops at the smallest rational value.
3430
 * Round up this value and check whether it is possible to attain
3431
 * a strictly greater integer value.
3432
 * If so, the variable is not restricted to a single integer value.
3433
 * Otherwise, the search stops at the greatest rational value.
3434
 * If rounding down this value yields a value that is different
3435
 * from rounding up the smallest rational value, then the variable
3436
 * cannot attain any integer value.  Mark the tableau empty.
3437
 * Otherwise, add an equality constraint that fixes the variable
3438
 * to the single integer value found.
3439
 */
3440
static isl_bool detect_constant_with_tmp(struct isl_tab *tab,
3441
  struct isl_tab_var *var, isl_int *target, isl_int *tmp)
3442
10.3k
{
3443
10.3k
  isl_bool reached;
3444
10.3k
  isl_vec *eq;
3445
10.3k
  int pos;
3446
10.3k
  isl_stat r;
3447
10.3k
3448
10.3k
  get_rounded_sample_value(tab, var, -1, target);
3449
10.3k
  isl_int_sub_ui(*target, *target, 1);
3450
10.3k
  reached = var_reaches(tab, var, -1, *target, tmp);
3451
10.3k
  if (reached < 0 || reached)
3452
7.95k
    return isl_bool_not(reached);
3453
2.41k
  get_rounded_sample_value(tab, var, 1, target);
3454
2.41k
  isl_int_add_ui(*target, *target, 1);
3455
2.41k
  reached = var_reaches(tab, var, 1, *target, tmp);
3456
2.41k
  if (reached < 0 || reached)
3457
2.12k
    return isl_bool_not(reached);
3458
295
  get_rounded_sample_value(tab, var, -1, tmp);
3459
295
  isl_int_sub_ui(*target, *target, 1);
3460
295
  if (isl_int_ne(*target, *tmp)) {
3461
0
    if (isl_tab_mark_empty(tab) < 0)
3462
0
      return isl_bool_error;
3463
0
    return isl_bool_false;
3464
0
  }
3465
295
3466
295
  if (isl_tab_extend_cons(tab, 1) < 0)
3467
0
    return isl_bool_error;
3468
295
  eq = isl_vec_alloc(isl_tab_get_ctx(tab), 1 + tab->n_var);
3469
295
  if (!eq)
3470
0
    return isl_bool_error;
3471
295
  pos = var - tab->var;
3472
295
  isl_seq_clr(eq->el + 1, tab->n_var);
3473
295
  isl_int_set_si(eq->el[1 + pos], -1);
3474
295
  isl_int_set(eq->el[0], *target);
3475
295
  r = isl_tab_add_eq(tab, eq->el);
3476
295
  isl_vec_free(eq);
3477
295
3478
295
  return r < 0 ? 
isl_bool_error0
: isl_bool_true;
3479
295
}
3480
3481
/* Check if variable "var" of "tab" can only attain a single (integer)
3482
 * value, and, if so, add an equality constraint to fix the variable
3483
 * to this single value and store the result in "value" (if "value"
3484
 * is not NULL).
3485
 *
3486
 * If the current sample value involves the big parameter,
3487
 * then the variable cannot have a fixed integer value.
3488
 * If the variable is already fixed to a single value by its row, then
3489
 * there is no need to add another equality constraint.
3490
 *
3491
 * Otherwise, allocate some temporary variables and continue
3492
 * with detect_constant_with_tmp.
3493
 */
3494
static isl_bool get_constant(struct isl_tab *tab, struct isl_tab_var *var,
3495
  isl_int *value)
3496
10.5k
{
3497
10.5k
  isl_int target, tmp;
3498
10.5k
  isl_bool is_cst;
3499
10.5k
3500
10.5k
  if (var->is_row && 
row_is_big(tab, var->index)4.36k
)
3501
0
    return isl_bool_false;
3502
10.5k
  is_cst = is_constant(tab, var, value);
3503
10.5k
  if (is_cst < 0 || is_cst)
3504
219
    return is_cst;
3505
10.3k
3506
10.3k
  if (!value)
3507
10.3k
    
isl_int_init7.58k
(target);
3508
10.3k
  isl_int_init(tmp);
3509
10.3k
3510
10.3k
  is_cst = detect_constant_with_tmp(tab, var,
3511
10.3k
              value ? 
value2.78k
:
&target7.58k
, &tmp);
3512
10.3k
3513
10.3k
  isl_int_clear(tmp);
3514
10.3k
  if (!value)
3515
10.3k
    
isl_int_clear7.58k
(target);
3516
10.3k
3517
10.3k
  return is_cst;
3518
10.3k
}
3519
3520
/* Check if variable "var" of "tab" can only attain a single (integer)
3521
 * value, and, if so, add an equality constraint to fix the variable
3522
 * to this single value and store the result in "value" (if "value"
3523
 * is not NULL).
3524
 *
3525
 * For rational tableaus, nothing needs to be done.
3526
 */
3527
isl_bool isl_tab_is_constant(struct isl_tab *tab, int var, isl_int *value)
3528
2.78k
{
3529
2.78k
  if (!tab)
3530
0
    return isl_bool_error;
3531
2.78k
  if (var < 0 || var >= tab->n_var)
3532
2.78k
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_invalid,
3533
2.78k
      "position out of bounds", return isl_bool_error);
3534
2.78k
  if (tab->rational)
3535
0
    return isl_bool_false;
3536
2.78k
3537
2.78k
  return get_constant(tab, &tab->var[var], value);
3538
2.78k
}
3539
3540
/* Check if any of the variables of "tab" can only attain a single (integer)
3541
 * value, and, if so, add equality constraints to fix those variables
3542
 * to these single values.
3543
 *
3544
 * For rational tableaus, nothing needs to be done.
3545
 */
3546
isl_stat isl_tab_detect_constants(struct isl_tab *tab)
3547
977
{
3548
977
  int i;
3549
977
3550
977
  if (!tab)
3551
0
    return isl_stat_error;
3552
977
  if (tab->rational)
3553
0
    return isl_stat_ok;
3554
977
3555
8.78k
  
for (i = 0; 977
i < tab->n_var;
++i7.80k
) {
3556
7.80k
    if (get_constant(tab, &tab->var[i], NULL) < 0)
3557
0
      return isl_stat_error;
3558
7.80k
  }
3559
977
3560
977
  return isl_stat_ok;
3561
977
}
3562
3563
/* Take a snapshot of the tableau that can be restored by a call to
3564
 * isl_tab_rollback.
3565
 */
3566
struct isl_tab_undo *isl_tab_snap(struct isl_tab *tab)
3567
2.21M
{
3568
2.21M
  if (!tab)
3569
0
    return NULL;
3570
2.21M
  tab->need_undo = 1;
3571
2.21M
  return tab->top;
3572
2.21M
}
3573
3574
/* Does "tab" need to keep track of undo information?
3575
 * That is, was a snapshot taken that may need to be restored?
3576
 */
3577
isl_bool isl_tab_need_undo(struct isl_tab *tab)
3578
176
{
3579
176
  if (!tab)
3580
0
    return isl_bool_error;
3581
176
3582
176
  return tab->need_undo;
3583
176
}
3584
3585
/* Remove all tracking of undo information from "tab", invalidating
3586
 * any snapshots that may have been taken of the tableau.
3587
 * Since all snapshots have been invalidated, there is also
3588
 * no need to start keeping track of undo information again.
3589
 */
3590
void isl_tab_clear_undo(struct isl_tab *tab)
3591
176
{
3592
176
  if (!tab)
3593
0
    return;
3594
176
3595
176
  free_undo(tab);
3596
176
  tab->need_undo = 0;
3597
176
}
3598
3599
/* Undo the operation performed by isl_tab_relax.
3600
 */
3601
static isl_stat unrelax(struct isl_tab *tab, struct isl_tab_var *var)
3602
  WARN_UNUSED;
3603
static isl_stat unrelax(struct isl_tab *tab, struct isl_tab_var *var)
3604
3.00k
{
3605
3.00k
  unsigned off = 2 + tab->M;
3606
3.00k
3607
3.00k
  if (!var->is_row && 
!max_is_manifestly_unbounded(tab, var)3.00k
)
3608
1.74k
    if (to_row(tab, var, 1) < 0)
3609
0
      return isl_stat_error;
3610
3.00k
3611
3.00k
  if (var->is_row) {
3612
1.75k
    isl_int_sub(tab->mat->row[var->index][1],
3613
1.75k
        tab->mat->row[var->index][1], tab->mat->row[var->index][0]);
3614
1.75k
    if (var->is_nonneg) {
3615
1.75k
      int sgn = restore_row(tab, var);
3616
1.75k
      isl_assert(tab->mat->ctx, sgn >= 0,
3617
1.75k
        return isl_stat_error);
3618
1.75k
    }
3619
1.75k
  } else {
3620
1.25k
    int i;
3621
1.25k
3622
8.42k
    for (i = 0; i < tab->n_row; 
++i7.17k
) {
3623
7.17k
      if (isl_int_is_zero(tab->mat->row[i][off + var->index]))
3624
7.17k
        
continue5.85k
;
3625
1.31k
      isl_int_add(tab->mat->row[i][1], tab->mat->row[i][1],
3626
1.31k
          tab->mat->row[i][off + var->index]);
3627
1.31k
    }
3628
1.25k
3629
1.25k
  }
3630
3.00k
3631
3.00k
  return isl_stat_ok;
3632
3.00k
}
3633
3634
/* Undo the operation performed by isl_tab_unrestrict.
3635
 *
3636
 * In particular, mark the variable as being non-negative and make
3637
 * sure the sample value respects this constraint.
3638
 */
3639
static isl_stat ununrestrict(struct isl_tab *tab, struct isl_tab_var *var)
3640
851
{
3641
851
  var->is_nonneg = 1;
3642
851
3643
851
  if (var->is_row && 
restore_row(tab, var) < -1788
)
3644
0
    return isl_stat_error;
3645
851
3646
851
  return isl_stat_ok;
3647
851
}
3648
3649
/* Unmark the last redundant row in "tab" as being redundant.
3650
 * This undoes part of the modifications performed by isl_tab_mark_redundant.
3651
 * In particular, remove the redundant mark and make
3652
 * sure the sample value respects the constraint again.
3653
 * A variable that is marked non-negative by isl_tab_mark_redundant
3654
 * is covered by a separate undo record.
3655
 */
3656
static isl_stat restore_last_redundant(struct isl_tab *tab)
3657
1.37M
{
3658
1.37M
  struct isl_tab_var *var;
3659
1.37M
3660
1.37M
  if (tab->n_redundant < 1)
3661
1.37M
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_internal,
3662
1.37M
      "no redundant rows", return isl_stat_error);
3663
1.37M
3664
1.37M
  var = isl_tab_var_from_row(tab, tab->n_redundant - 1);
3665
1.37M
  var->is_redundant = 0;
3666
1.37M
  tab->n_redundant--;
3667
1.37M
  restore_row(tab, var);
3668
1.37M
3669
1.37M
  return isl_stat_ok;
3670
1.37M
}
3671
3672
static isl_stat perform_undo_var(struct isl_tab *tab, struct isl_tab_undo *undo)
3673
  WARN_UNUSED;
3674
static isl_stat perform_undo_var(struct isl_tab *tab, struct isl_tab_undo *undo)
3675
4.66M
{
3676
4.66M
  struct isl_tab_var *var = var_from_index(tab, undo->u.var_index);
3677
4.66M
  switch (undo->type) {
3678
4.66M
  case isl_tab_undo_nonneg:
3679
896k
    var->is_nonneg = 0;
3680
896k
    break;
3681
4.66M
  case isl_tab_undo_redundant:
3682
1.06M
    if (!var->is_row || var->index != tab->n_redundant - 1)
3683
1.06M
      
isl_die0
(isl_tab_get_ctx(tab), isl_error_internal,
3684
1.06M
        "not undoing last redundant row",
3685
1.06M
        return isl_stat_error);
3686
1.06M
    return restore_last_redundant(tab);
3687
1.06M
  case isl_tab_undo_freeze:
3688
311k
    var->frozen = 0;
3689
311k
    break;
3690
1.06M
  case isl_tab_undo_zero:
3691
230k
    var->is_zero = 0;
3692
230k
    if (!var->is_row)
3693
230k
      tab->n_dead--;
3694
230k
    break;
3695
2.16M
  case isl_tab_undo_allocate:
3696
2.16M
    if (undo->u.var_index >= 0) {
3697
5.51k
      isl_assert(tab->mat->ctx, !var->is_row,
3698
5.51k
        return isl_stat_error);
3699
5.51k
      return drop_col(tab, var->index);
3700
2.15M
    }
3701
2.15M
    if (!var->is_row) {
3702
270k
      if (!max_is_manifestly_unbounded(tab, var)) {
3703
234k
        if (to_row(tab, var, 1) < 0)
3704
0
          return isl_stat_error;
3705
36.7k
      } else if (!min_is_manifestly_unbounded(tab, var)) {
3706
10.9k
        if (to_row(tab, var, -1) < 0)
3707
0
          return isl_stat_error;
3708
25.7k
      } else
3709
25.7k
        if (to_row(tab, var, 0) < 0)
3710
0
          return isl_stat_error;
3711
2.15M
    }
3712
2.15M
    return drop_row(tab, var->index);
3713
2.15M
  case isl_tab_undo_relax:
3714
3.00k
    return unrelax(tab, var);
3715
2.15M
  case isl_tab_undo_unrestrict:
3716
851
    return ununrestrict(tab, var);
3717
2.15M
  default:
3718
0
    isl_die(tab->mat->ctx, isl_error_internal,
3719
4.66M
      "perform_undo_var called on invalid undo record",
3720
4.66M
      return isl_stat_error);
3721
4.66M
  }
3722
4.66M
3723
4.66M
  
return isl_stat_ok1.43M
;
3724
4.66M
}
3725
3726
/* Restore all rows that have been marked redundant by isl_tab_mark_redundant
3727
 * and that have been preserved in the tableau.
3728
 * Note that isl_tab_mark_redundant may also have marked some variables
3729
 * as being non-negative before marking them redundant.  These need
3730
 * to be removed as well as otherwise some constraints could end up
3731
 * getting marked redundant with respect to the variable.
3732
 */
3733
isl_stat isl_tab_restore_redundant(struct isl_tab *tab)
3734
207k
{
3735
207k
  if (!tab)
3736
0
    return isl_stat_error;
3737
207k
3738
207k
  if (tab->need_undo)
3739
207k
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_invalid,
3740
207k
      "manually restoring redundant constraints "
3741
207k
      "interferes with undo history",
3742
207k
      return isl_stat_error);
3743
207k
3744
520k
  
while (207k
tab->n_redundant > 0) {
3745
313k
    if (tab->row_var[tab->n_redundant - 1] >= 0) {
3746
284k
      struct isl_tab_var *var;
3747
284k
3748
284k
      var = isl_tab_var_from_row(tab, tab->n_redundant - 1);
3749
284k
      var->is_nonneg = 0;
3750
284k
    }
3751
313k
    restore_last_redundant(tab);
3752
313k
  }
3753
207k
  return isl_stat_ok;
3754
207k
}
3755
3756
/* Undo the addition of an integer division to the basic map representation
3757
 * of "tab" in position "pos".
3758
 */
3759
static isl_stat drop_bmap_div(struct isl_tab *tab, int pos)
3760
3.20k
{
3761
3.20k
  int off;
3762
3.20k
3763
3.20k
  off = tab->n_var - isl_basic_map_dim(tab->bmap, isl_dim_div);
3764
3.20k
  if (isl_basic_map_drop_div(tab->bmap, pos - off) < 0)
3765
0
    return isl_stat_error;
3766
3.20k
  if (tab->samples) {
3767
626
    tab->samples = isl_mat_drop_cols(tab->samples, 1 + pos, 1);
3768
626
    if (!tab->samples)
3769
0
      return isl_stat_error;
3770
3.20k
  }
3771
3.20k
3772
3.20k
  return isl_stat_ok;
3773
3.20k
}
3774
3775
/* Restore the tableau to the state where the basic variables
3776
 * are those in "col_var".
3777
 * We first construct a list of variables that are currently in
3778
 * the basis, but shouldn't.  Then we iterate over all variables
3779
 * that should be in the basis and for each one that is currently
3780
 * not in the basis, we exchange it with one of the elements of the
3781
 * list constructed before.
3782
 * We can always find an appropriate variable to pivot with because
3783
 * the current basis is mapped to the old basis by a non-singular
3784
 * matrix and so we can never end up with a zero row.
3785
 */
3786
static int restore_basis(struct isl_tab *tab, int *col_var)
3787
576
{
3788
576
  int i, j;
3789
576
  int n_extra = 0;
3790
576
  int *extra = NULL;  /* current columns that contain bad stuff */
3791
576
  unsigned off = 2 + tab->M;
3792
576
3793
576
  extra = isl_alloc_array(tab->mat->ctx, int, tab->n_col);
3794
576
  if (tab->n_col && !extra)
3795
0
    goto error;
3796
6.90k
  
for (i = 0; 576
i < tab->n_col;
++i6.32k
) {
3797
52.2k
    for (j = 0; j < tab->n_col; 
++j45.9k
)
3798
50.7k
      if (tab->col_var[i] == col_var[j])
3799
4.85k
        break;
3800
6.32k
    if (j < tab->n_col)
3801
4.85k
      continue;
3802
1.46k
    extra[n_extra++] = i;
3803
1.46k
  }
3804
5.39k
  for (i = 0; i < tab->n_col && 
n_extra > 05.28k
;
++i4.82k
) {
3805
4.82k
    struct isl_tab_var *var;
3806
4.82k
    int row;
3807
4.82k
3808
37.3k
    for (j = 0; j < tab->n_col; 
++j32.5k
)
3809
35.9k
      if (col_var[i] == tab->col_var[j])
3810
3.35k
        break;
3811
4.82k
    if (j < tab->n_col)
3812
3.35k
      continue;
3813
1.46k
    var = var_from_index(tab, col_var[i]);
3814
1.46k
    row = var->index;
3815
1.86k
    for (j = 0; j < n_extra; 
++j396
)
3816
1.86k
      if (!isl_int_is_zero(tab->mat->row[row][off+extra[j]]))
3817
1.86k
        
break1.46k
;
3818
1.46k
    isl_assert(tab->mat->ctx, j < n_extra, goto error);
3819
1.46k
    if (isl_tab_pivot(tab, row, extra[j]) < 0)
3820
0
      goto error;
3821
1.46k
    extra[j] = extra[--n_extra];
3822
1.46k
  }
3823
576
3824
576
  free(extra);
3825
576
  return 0;
3826
0
error:
3827
0
  free(extra);
3828
0
  return -1;
3829
576
}
3830
3831
/* Remove all samples with index n or greater, i.e., those samples
3832
 * that were added since we saved this number of samples in
3833
 * isl_tab_save_samples.
3834
 */
3835
static void drop_samples_since(struct isl_tab *tab, int n)
3836
23.6k
{
3837
23.6k
  int i;
3838
23.6k
3839
29.1k
  for (i = tab->n_sample - 1; i >= 0 && 
tab->n_sample > n27.4k
;
--i5.49k
) {
3840
5.49k
    if (tab->sample_index[i] < n)
3841
1.99k
      continue;
3842
3.49k
3843
3.49k
    if (i != tab->n_sample - 1) {
3844
2.20k
      int t = tab->sample_index[tab->n_sample-1];
3845
2.20k
      tab->sample_index[tab->n_sample-1] = tab->sample_index[i];
3846
2.20k
      tab->sample_index[i] = t;
3847
2.20k
      isl_mat_swap_rows(tab->samples, tab->n_sample-1, i);
3848
2.20k
    }
3849
3.49k
    tab->n_sample--;
3850
3.49k
  }
3851
23.6k
}
3852
3853
static isl_stat perform_undo(struct isl_tab *tab, struct isl_tab_undo *undo)
3854
  WARN_UNUSED;
3855
static isl_stat perform_undo(struct isl_tab *tab, struct isl_tab_undo *undo)
3856
5.12M
{
3857
5.12M
  switch (undo->type) {
3858
5.12M
  case isl_tab_undo_rational:
3859
9.20k
    tab->rational = 0;
3860
9.20k
    break;
3861
5.12M
  case isl_tab_undo_empty:
3862
69.7k
    tab->empty = 0;
3863
69.7k
    break;
3864
5.12M
  case isl_tab_undo_nonneg:
3865
4.66M
  case isl_tab_undo_redundant:
3866
4.66M
  case isl_tab_undo_freeze:
3867
4.66M
  case isl_tab_undo_zero:
3868
4.66M
  case isl_tab_undo_allocate:
3869
4.66M
  case isl_tab_undo_relax:
3870
4.66M
  case isl_tab_undo_unrestrict:
3871
4.66M
    return perform_undo_var(tab, undo);
3872
4.66M
  case isl_tab_undo_bmap_eq:
3873
0
    return isl_basic_map_free_equality(tab->bmap, 1);
3874
4.66M
  case isl_tab_undo_bmap_ineq:
3875
348k
    return isl_basic_map_free_inequality(tab->bmap, 1);
3876
4.66M
  case isl_tab_undo_bmap_div:
3877
3.20k
    return drop_bmap_div(tab, undo->u.var_index);
3878
4.66M
  case isl_tab_undo_saved_basis:
3879
576
    if (restore_basis(tab, undo->u.col_var) < 0)
3880
0
      return isl_stat_error;
3881
576
    break;
3882
6.04k
  case isl_tab_undo_drop_sample:
3883
6.04k
    tab->n_outside--;
3884
6.04k
    break;
3885
23.6k
  case isl_tab_undo_saved_samples:
3886
23.6k
    drop_samples_since(tab, undo->u.n);
3887
23.6k
    break;
3888
2.41k
  case isl_tab_undo_callback:
3889
2.41k
    return undo->u.callback->run(undo->u.callback);
3890
576
  default:
3891
0
    isl_assert(tab->mat->ctx, 0, return isl_stat_error);
3892
5.12M
  }
3893
5.12M
  
return isl_stat_ok109k
;
3894
5.12M
}
3895
3896
/* Return the tableau to the state it was in when the snapshot "snap"
3897
 * was taken.
3898
 */
3899
int isl_tab_rollback(struct isl_tab *tab, struct isl_tab_undo *snap)
3900
2.05M
{
3901
2.05M
  struct isl_tab_undo *undo, *next;
3902
2.05M
3903
2.05M
  if (!tab)
3904
0
    return -1;
3905
2.05M
3906
2.05M
  tab->in_undo = 1;
3907
7.18M
  for (undo = tab->top; undo && undo != &tab->bottom; 
undo = next5.12M
) {
3908
5.89M
    next = undo->next;
3909
5.89M
    if (undo == snap)
3910
769k
      break;
3911
5.12M
    if (perform_undo(tab, undo) < 0) {
3912
0
      tab->top = undo;
3913
0
      free_undo(tab);
3914
0
      tab->in_undo = 0;
3915
0
      return -1;
3916
0
    }
3917
5.12M
    free_undo_record(undo);
3918
5.12M
  }
3919
2.05M
  tab->in_undo = 0;
3920
2.05M
  tab->top = undo;
3921
2.05M
  if (!undo)
3922
0
    return -1;
3923
2.05M
  return 0;
3924
2.05M
}
3925
3926
/* The given row "row" represents an inequality violated by all
3927
 * points in the tableau.  Check for some special cases of such
3928
 * separating constraints.
3929
 * In particular, if the row has been reduced to the constant -1,
3930
 * then we know the inequality is adjacent (but opposite) to
3931
 * an equality in the tableau.
3932
 * If the row has been reduced to r = c*(-1 -r'), with r' an inequality
3933
 * of the tableau and c a positive constant, then the inequality
3934
 * is adjacent (but opposite) to the inequality r'.
3935
 */
3936
static enum isl_ineq_type separation_type(struct isl_tab *tab, unsigned row)
3937
105k
{
3938
105k
  int pos;
3939
105k
  unsigned off = 2 + tab->M;
3940
105k
3941
105k
  if (tab->rational)
3942
9.20k
    return isl_ineq_separate;
3943
96.5k
3944
96.5k
  if (!isl_int_is_one(tab->mat->row[row][0]))
3945
96.5k
    
return isl_ineq_separate242
;
3946
96.3k
3947
96.3k
  pos = isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
3948
96.3k
          tab->n_col - tab->n_dead);
3949
96.3k
  if (pos == -1) {
3950
24.7k
    if (isl_int_is_negone(tab->mat->row[row][1]))
3951
24.7k
      
return isl_ineq_adj_eq21.4k
;
3952
3.34k
    else
3953
3.34k
      return isl_ineq_separate;
3954
71.5k
  }
3955
71.5k
3956
71.5k
  if (!isl_int_eq(tab->mat->row[row][1],
3957
71.5k
      tab->mat->row[row][off + tab->n_dead + pos]))
3958
71.5k
    
return isl_ineq_separate21.5k
;
3959
50.0k
3960
50.0k
  pos = isl_seq_first_non_zero(
3961
50.0k
      tab->mat->row[row] + off + tab->n_dead + pos + 1,
3962
50.0k
      tab->n_col - tab->n_dead - pos - 1);
3963
50.0k
3964
50.0k
  return pos == -1 ? 
isl_ineq_adj_ineq46.8k
:
isl_ineq_separate3.13k
;
3965
50.0k
}
3966
3967
/* Check the effect of inequality "ineq" on the tableau "tab".
3968
 * The result may be
3969
 *  isl_ineq_redundant: satisfied by all points in the tableau
3970
 *  isl_ineq_separate:  satisfied by no point in the tableau
3971
 *  isl_ineq_cut:   satisfied by some by not all points
3972
 *  isl_ineq_adj_eq:  adjacent to an equality
3973
 *  isl_ineq_adj_ineq:  adjacent to an inequality.
3974
 */
3975
enum isl_ineq_type isl_tab_ineq_type(struct isl_tab *tab, isl_int *ineq)
3976
785k
{
3977
785k
  enum isl_ineq_type type = isl_ineq_error;
3978
785k
  struct isl_tab_undo *snap = NULL;
3979
785k
  int con;
3980
785k
  int row;
3981
785k
3982
785k
  if (!tab)
3983
0
    return isl_ineq_error;
3984
785k
3985
785k
  if (isl_tab_extend_cons(tab, 1) < 0)
3986
0
    return isl_ineq_error;
3987
785k
3988
785k
  snap = isl_tab_snap(tab);
3989
785k
3990
785k
  con = isl_tab_add_row(tab, ineq);
3991
785k
  if (con < 0)
3992
0
    goto error;
3993
785k
3994
785k
  row = tab->con[con].index;
3995
785k
  if (isl_tab_row_is_redundant(tab, row))
3996
0
    type = isl_ineq_redundant;
3997
785k
  else if (isl_int_is_neg(tab->mat->row[row][1]) &&
3998
785k
     
(184k
tab->rational184k
||
3999
184k
        
isl_int_abs_ge163k
(tab->mat->row[row][1],
4000
184k
           tab->mat->row[row][0]))) {
4001
183k
    int nonneg = at_least_zero(tab, &tab->con[con]);
4002
183k
    if (nonneg < 0)
4003
0
      goto error;
4004
183k
    if (nonneg)
4005
77.8k
      type = isl_ineq_cut;
4006
105k
    else
4007
105k
      type = separation_type(tab, row);
4008
602k
  } else {
4009
602k
    int red = con_is_redundant(tab, &tab->con[con]);
4010
602k
    if (red < 0)
4011
0
      goto error;
4012
602k
    if (!red)
4013
73.3k
      type = isl_ineq_cut;
4014
528k
    else
4015
528k
      type = isl_ineq_redundant;
4016
602k
  }
4017
785k
4018
785k
  if (isl_tab_rollback(tab, snap))
4019
0
    return isl_ineq_error;
4020
785k
  return type;
4021
0
error:
4022
0
  return isl_ineq_error;
4023
785k
}
4024
4025
isl_stat isl_tab_track_bmap(struct isl_tab *tab, __isl_take isl_basic_map *bmap)
4026
325k
{
4027
325k
  bmap = isl_basic_map_cow(bmap);
4028
325k
  if (!tab || !bmap)
4029
0
    goto error;
4030
325k
4031
325k
  if (tab->empty) {
4032
4.38k
    bmap = isl_basic_map_set_to_empty(bmap);
4033
4.38k
    if (!bmap)
4034
0
      goto error;
4035
4.38k
    tab->bmap = bmap;
4036
4.38k
    return isl_stat_ok;
4037
4.38k
  }
4038
321k
4039
321k
  isl_assert(tab->mat->ctx, tab->n_eq == bmap->n_eq, goto error);
4040
321k
  isl_assert(tab->mat->ctx,
4041
321k
        tab->n_con == bmap->n_eq + bmap->n_ineq, goto error);
4042
321k
4043
321k
  tab->bmap = bmap;
4044
321k
4045
321k
  return isl_stat_ok;
4046
0
error:
4047
0
  isl_basic_map_free(bmap);
4048
0
  return isl_stat_error;
4049
321k
}
4050
4051
isl_stat isl_tab_track_bset(struct isl_tab *tab, __isl_take isl_basic_set *bset)
4052
869
{
4053
869
  return isl_tab_track_bmap(tab, bset_to_bmap(bset));
4054
869
}
4055
4056
__isl_keep isl_basic_set *isl_tab_peek_bset(struct isl_tab *tab)
4057
31.2k
{
4058
31.2k
  if (!tab)
4059
0
    return NULL;
4060
31.2k
4061
31.2k
  return bset_from_bmap(tab->bmap);
4062
31.2k
}
4063
4064
static void isl_tab_print_internal(__isl_keep struct isl_tab *tab,
4065
  FILE *out, int indent)
4066
0
{
4067
0
  unsigned r, c;
4068
0
  int i;
4069
0
4070
0
  if (!tab) {
4071
0
    fprintf(out, "%*snull tab\n", indent, "");
4072
0
    return;
4073
0
  }
4074
0
  fprintf(out, "%*sn_redundant: %d, n_dead: %d", indent, "",
4075
0
    tab->n_redundant, tab->n_dead);
4076
0
  if (tab->rational)
4077
0
    fprintf(out, ", rational");
4078
0
  if (tab->empty)
4079
0
    fprintf(out, ", empty");
4080
0
  fprintf(out, "\n");
4081
0
  fprintf(out, "%*s[", indent, "");
4082
0
  for (i = 0; i < tab->n_var; ++i) {
4083
0
    if (i)
4084
0
      fprintf(out, (i == tab->n_param ||
4085
0
              i == tab->n_var - tab->n_div) ? "; "
4086
0
                    : ", ");
4087
0
    fprintf(out, "%c%d%s", tab->var[i].is_row ? 'r' : 'c',
4088
0
          tab->var[i].index,
4089
0
          tab->var[i].is_zero ? " [=0]" :
4090
0
          tab->var[i].is_redundant ? " [R]" : "");
4091
0
  }
4092
0
  fprintf(out, "]\n");
4093
0
  fprintf(out, "%*s[", indent, "");
4094
0
  for (i = 0; i < tab->n_con; ++i) {
4095
0
    if (i)
4096
0
      fprintf(out, ", ");
4097
0
    fprintf(out, "%c%d%s", tab->con[i].is_row ? 'r' : 'c',
4098
0
          tab->con[i].index,
4099
0
          tab->con[i].is_zero ? " [=0]" :
4100
0
          tab->con[i].is_redundant ? " [R]" : "");
4101
0
  }
4102
0
  fprintf(out, "]\n");
4103
0
  fprintf(out, "%*s[", indent, "");
4104
0
  for (i = 0; i < tab->n_row; ++i) {
4105
0
    const char *sign = "";
4106
0
    if (i)
4107
0
      fprintf(out, ", ");
4108
0
    if (tab->row_sign) {
4109
0
      if (tab->row_sign[i] == isl_tab_row_unknown)
4110
0
        sign = "?";
4111
0
      else if (tab->row_sign[i] == isl_tab_row_neg)
4112
0
        sign = "-";
4113
0
      else if (tab->row_sign[i] == isl_tab_row_pos)
4114
0
        sign = "+";
4115
0
      else
4116
0
        sign = "+-";
4117
0
    }
4118
0
    fprintf(out, "r%d: %d%s%s", i, tab->row_var[i],
4119
0
        isl_tab_var_from_row(tab, i)->is_nonneg ? " [>=0]" : "", sign);
4120
0
  }
4121
0
  fprintf(out, "]\n");
4122
0
  fprintf(out, "%*s[", indent, "");
4123
0
  for (i = 0; i < tab->n_col; ++i) {
4124
0
    if (i)
4125
0
      fprintf(out, ", ");
4126
0
    fprintf(out, "c%d: %d%s", i, tab->col_var[i],
4127
0
        var_from_col(tab, i)->is_nonneg ? " [>=0]" : "");
4128
0
  }
4129
0
  fprintf(out, "]\n");
4130
0
  r = tab->mat->n_row;
4131
0
  tab->mat->n_row = tab->n_row;
4132
0
  c = tab->mat->n_col;
4133
0
  tab->mat->n_col = 2 + tab->M + tab->n_col;
4134
0
  isl_mat_print_internal(tab->mat, out, indent);
4135
0
  tab->mat->n_row = r;
4136
0
  tab->mat->n_col = c;
4137
0
  if (tab->bmap)
4138
0
    isl_basic_map_print_internal(tab->bmap, out, indent);
4139
0
}
4140
4141
void isl_tab_dump(__isl_keep struct isl_tab *tab)
4142
0
{
4143
0
  isl_tab_print_internal(tab, stderr, 0);
4144
0
}