Coverage Report

Created: 2018-08-20 19:24

/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
732k
{
36
732k
  int i;
37
732k
  struct isl_tab *tab;
38
732k
  unsigned off = 2 + M;
39
732k
40
732k
  tab = isl_calloc_type(ctx, struct isl_tab);
41
732k
  if (!tab)
42
0
    return NULL;
43
732k
  tab->mat = isl_mat_alloc(ctx, n_row, off + n_var);
44
732k
  if (!tab->mat)
45
0
    goto error;
46
732k
  tab->var = isl_alloc_array(ctx, struct isl_tab_var, n_var);
47
732k
  if (n_var && 
!tab->var731k
)
48
0
    goto error;
49
732k
  tab->con = isl_alloc_array(ctx, struct isl_tab_var, n_row);
50
732k
  if (n_row && 
!tab->con732k
)
51
0
    goto error;
52
732k
  tab->col_var = isl_alloc_array(ctx, int, n_var);
53
732k
  if (n_var && 
!tab->col_var731k
)
54
0
    goto error;
55
732k
  tab->row_var = isl_alloc_array(ctx, int, n_row);
56
732k
  if (n_row && 
!tab->row_var732k
)
57
0
    goto error;
58
4.10M
  
for (i = 0; 732k
i < n_var;
++i3.37M
) {
59
3.37M
    tab->var[i].index = i;
60
3.37M
    tab->var[i].is_row = 0;
61
3.37M
    tab->var[i].is_nonneg = 0;
62
3.37M
    tab->var[i].is_zero = 0;
63
3.37M
    tab->var[i].is_redundant = 0;
64
3.37M
    tab->var[i].frozen = 0;
65
3.37M
    tab->var[i].negated = 0;
66
3.37M
    tab->col_var[i] = i;
67
3.37M
  }
68
732k
  tab->n_row = 0;
69
732k
  tab->n_con = 0;
70
732k
  tab->n_eq = 0;
71
732k
  tab->max_con = n_row;
72
732k
  tab->n_col = n_var;
73
732k
  tab->n_var = n_var;
74
732k
  tab->max_var = n_var;
75
732k
  tab->n_param = 0;
76
732k
  tab->n_div = 0;
77
732k
  tab->n_dead = 0;
78
732k
  tab->n_redundant = 0;
79
732k
  tab->strict_redundant = 0;
80
732k
  tab->need_undo = 0;
81
732k
  tab->rational = 0;
82
732k
  tab->empty = 0;
83
732k
  tab->in_undo = 0;
84
732k
  tab->M = M;
85
732k
  tab->cone = 0;
86
732k
  tab->bottom.type = isl_tab_undo_bottom;
87
732k
  tab->bottom.next = NULL;
88
732k
  tab->top = &tab->bottom;
89
732k
90
732k
  tab->n_zero = 0;
91
732k
  tab->n_unbounded = 0;
92
732k
  tab->basis = NULL;
93
732k
94
732k
  return tab;
95
0
error:
96
0
  isl_tab_free(tab);
97
0
  return NULL;
98
732k
}
99
100
isl_ctx *isl_tab_get_ctx(struct isl_tab *tab)
101
3.37M
{
102
3.37M
  return tab ? isl_mat_get_ctx(tab->mat) : NULL;
103
3.37M
}
104
105
int isl_tab_extend_cons(struct isl_tab *tab, unsigned n_new)
106
808k
{
107
808k
  unsigned off;
108
808k
109
808k
  if (!tab)
110
0
    return -1;
111
808k
112
808k
  off = 2 + tab->M;
113
808k
114
808k
  if (tab->max_con < tab->n_con + n_new) {
115
85.2k
    struct isl_tab_var *con;
116
85.2k
117
85.2k
    con = isl_realloc_array(tab->mat->ctx, tab->con,
118
85.2k
            struct isl_tab_var, tab->max_con + n_new);
119
85.2k
    if (!con)
120
0
      return -1;
121
85.2k
    tab->con = con;
122
85.2k
    tab->max_con += n_new;
123
85.2k
  }
124
808k
  if (tab->mat->n_row < tab->n_row + n_new) {
125
87.5k
    int *row_var;
126
87.5k
127
87.5k
    tab->mat = isl_mat_extend(tab->mat,
128
87.5k
          tab->n_row + n_new, off + tab->n_col);
129
87.5k
    if (!tab->mat)
130
0
      return -1;
131
87.5k
    row_var = isl_realloc_array(tab->mat->ctx, tab->row_var,
132
87.5k
              int, tab->mat->n_row);
133
87.5k
    if (!row_var)
134
0
      return -1;
135
87.5k
    tab->row_var = row_var;
136
87.5k
    if (tab->row_sign) {
137
219
      enum isl_tab_row_sign *s;
138
219
      s = isl_realloc_array(tab->mat->ctx, tab->row_sign,
139
219
          enum isl_tab_row_sign, tab->mat->n_row);
140
219
      if (!s)
141
0
        return -1;
142
219
      tab->row_sign = s;
143
219
    }
144
87.5k
  }
145
808k
  return 0;
146
808k
}
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
7.47k
{
153
7.47k
  struct isl_tab_var *var;
154
7.47k
  unsigned off = 2 + tab->M;
155
7.47k
156
7.47k
  if (tab->max_var < tab->n_var + n_new) {
157
5.00k
    var = isl_realloc_array(tab->mat->ctx, tab->var,
158
5.00k
            struct isl_tab_var, tab->n_var + n_new);
159
5.00k
    if (!var)
160
0
      return -1;
161
5.00k
    tab->var = var;
162
5.00k
    tab->max_var = tab->n_var + n_new;
163
5.00k
  }
164
7.47k
165
7.47k
  if (tab->mat->n_col < off + tab->n_col + n_new) {
166
3.26k
    int *p;
167
3.26k
168
3.26k
    tab->mat = isl_mat_extend(tab->mat,
169
3.26k
            tab->mat->n_row, off + tab->n_col + n_new);
170
3.26k
    if (!tab->mat)
171
0
      return -1;
172
3.26k
    p = isl_realloc_array(tab->mat->ctx, tab->col_var,
173
3.26k
              int, tab->n_col + n_new);
174
3.26k
    if (!p)
175
0
      return -1;
176
3.26k
    tab->col_var = p;
177
3.26k
  }
178
7.47k
179
7.47k
  return 0;
180
7.47k
}
181
182
static void free_undo_record(struct isl_tab_undo *undo)
183
3.34M
{
184
3.34M
  switch (undo->type) {
185
3.34M
  case isl_tab_undo_saved_basis:
186
758
    free(undo->u.col_var);
187
758
    break;
188
3.34M
  
default:;3.34M
189
3.34M
  }
190
3.34M
  free(undo);
191
3.34M
}
192
193
static void free_undo(struct isl_tab *tab)
194
738k
{
195
738k
  struct isl_tab_undo *undo, *next;
196
738k
197
1.29M
  for (undo = tab->top; undo && undo != &tab->bottom; 
undo = next559k
) {
198
559k
    next = undo->next;
199
559k
    free_undo_record(undo);
200
559k
  }
201
738k
  tab->top = undo;
202
738k
}
203
204
void isl_tab_free(struct isl_tab *tab)
205
771k
{
206
771k
  if (!tab)
207
32.9k
    return;
208
738k
  free_undo(tab);
209
738k
  isl_mat_free(tab->mat);
210
738k
  isl_vec_free(tab->dual);
211
738k
  isl_basic_map_free(tab->bmap);
212
738k
  free(tab->var);
213
738k
  free(tab->con);
214
738k
  free(tab->row_var);
215
738k
  free(tab->col_var);
216
738k
  free(tab->row_sign);
217
738k
  isl_mat_free(tab->samples);
218
738k
  free(tab->sample_index);
219
738k
  isl_mat_free(tab->basis);
220
738k
  free(tab);
221
738k
}
222
223
struct isl_tab *isl_tab_dup(struct isl_tab *tab)
224
2.67k
{
225
2.67k
  int i;
226
2.67k
  struct isl_tab *dup;
227
2.67k
  unsigned off;
228
2.67k
229
2.67k
  if (!tab)
230
0
    return NULL;
231
2.67k
232
2.67k
  off = 2 + tab->M;
233
2.67k
  dup = isl_calloc_type(tab->mat->ctx, struct isl_tab);
234
2.67k
  if (!dup)
235
0
    return NULL;
236
2.67k
  dup->mat = isl_mat_dup(tab->mat);
237
2.67k
  if (!dup->mat)
238
0
    goto error;
239
2.67k
  dup->var = isl_alloc_array(tab->mat->ctx, struct isl_tab_var, tab->max_var);
240
2.67k
  if (tab->max_var && !dup->var)
241
0
    goto error;
242
25.3k
  
for (i = 0; 2.67k
i < tab->n_var;
++i22.6k
)
243
22.6k
    dup->var[i] = tab->var[i];
244
2.67k
  dup->con = isl_alloc_array(tab->mat->ctx, struct isl_tab_var, tab->max_con);
245
2.67k
  if (tab->max_con && !dup->con)
246
0
    goto error;
247
28.0k
  
for (i = 0; 2.67k
i < tab->n_con;
++i25.4k
)
248
25.4k
    dup->con[i] = tab->con[i];
249
2.67k
  dup->col_var = isl_alloc_array(tab->mat->ctx, int, tab->mat->n_col - off);
250
2.67k
  if ((tab->mat->n_col - off) && !dup->col_var)
251
0
    goto error;
252
12.9k
  
for (i = 0; 2.67k
i < tab->n_col;
++i10.2k
)
253
10.2k
    dup->col_var[i] = tab->col_var[i];
254
2.67k
  dup->row_var = isl_alloc_array(tab->mat->ctx, int, tab->mat->n_row);
255
2.67k
  if (tab->mat->n_row && !dup->row_var)
256
0
    goto error;
257
27.9k
  
for (i = 0; 2.67k
i < tab->n_row;
++i25.2k
)
258
25.2k
    dup->row_var[i] = tab->row_var[i];
259
2.67k
  if (tab->row_sign) {
260
2.66k
    dup->row_sign = isl_alloc_array(tab->mat->ctx, enum isl_tab_row_sign,
261
2.66k
            tab->mat->n_row);
262
2.66k
    if (tab->mat->n_row && !dup->row_sign)
263
0
      goto error;
264
27.8k
    
for (i = 0; 2.66k
i < tab->n_row;
++i25.2k
)
265
25.2k
      dup->row_sign[i] = tab->row_sign[i];
266
2.66k
  }
267
2.67k
  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.67k
  dup->n_row = tab->n_row;
279
2.67k
  dup->n_con = tab->n_con;
280
2.67k
  dup->n_eq = tab->n_eq;
281
2.67k
  dup->max_con = tab->max_con;
282
2.67k
  dup->n_col = tab->n_col;
283
2.67k
  dup->n_var = tab->n_var;
284
2.67k
  dup->max_var = tab->max_var;
285
2.67k
  dup->n_param = tab->n_param;
286
2.67k
  dup->n_div = tab->n_div;
287
2.67k
  dup->n_dead = tab->n_dead;
288
2.67k
  dup->n_redundant = tab->n_redundant;
289
2.67k
  dup->rational = tab->rational;
290
2.67k
  dup->empty = tab->empty;
291
2.67k
  dup->strict_redundant = 0;
292
2.67k
  dup->need_undo = 0;
293
2.67k
  dup->in_undo = 0;
294
2.67k
  dup->M = tab->M;
295
2.67k
  tab->cone = tab->cone;
296
2.67k
  dup->bottom.type = isl_tab_undo_bottom;
297
2.67k
  dup->bottom.next = NULL;
298
2.67k
  dup->top = &dup->bottom;
299
2.67k
300
2.67k
  dup->n_zero = tab->n_zero;
301
2.67k
  dup->n_unbounded = tab->n_unbounded;
302
2.67k
  dup->basis = isl_mat_dup(tab->basis);
303
2.67k
304
2.67k
  return dup;
305
0
error:
306
0
  isl_tab_free(dup);
307
0
  return NULL;
308
2.67k
}
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
3.15k
{
328
3.15k
  int i;
329
3.15k
  struct isl_mat *prod;
330
3.15k
  unsigned n;
331
3.15k
332
3.15k
  prod = isl_mat_alloc(mat1->ctx, mat1->n_row + mat2->n_row,
333
3.15k
          off + col1 + col2);
334
3.15k
  if (!prod)
335
0
    return NULL;
336
3.15k
337
3.15k
  n = 0;
338
10.9k
  for (i = 0; i < r1; 
++i7.78k
) {
339
7.78k
    isl_seq_cpy(prod->row[n + i], mat1->row[i], off + d1);
340
7.78k
    isl_seq_clr(prod->row[n + i] + off + d1, d2);
341
7.78k
    isl_seq_cpy(prod->row[n + i] + off + d1 + d2,
342
7.78k
        mat1->row[i] + off + d1, col1 - d1);
343
7.78k
    isl_seq_clr(prod->row[n + i] + off + col1 + d1, col2 - d2);
344
7.78k
  }
345
3.15k
346
3.15k
  n += r1;
347
10.9k
  for (i = 0; i < r2; 
++i7.78k
) {
348
7.78k
    isl_seq_cpy(prod->row[n + i], mat2->row[i], off);
349
7.78k
    isl_seq_clr(prod->row[n + i] + off, d1);
350
7.78k
    isl_seq_cpy(prod->row[n + i] + off + d1,
351
7.78k
          mat2->row[i] + off, d2);
352
7.78k
    isl_seq_clr(prod->row[n + i] + off + d1 + d2, col1 - d1);
353
7.78k
    isl_seq_cpy(prod->row[n + i] + off + col1 + d1,
354
7.78k
          mat2->row[i] + off + d2, col2 - d2);
355
7.78k
  }
356
3.15k
357
3.15k
  n += r2;
358
36.4k
  for (i = 0; i < row1 - r1; 
++i33.3k
) {
359
33.3k
    isl_seq_cpy(prod->row[n + i], mat1->row[r1 + i], off + d1);
360
33.3k
    isl_seq_clr(prod->row[n + i] + off + d1, d2);
361
33.3k
    isl_seq_cpy(prod->row[n + i] + off + d1 + d2,
362
33.3k
        mat1->row[r1 + i] + off + d1, col1 - d1);
363
33.3k
    isl_seq_clr(prod->row[n + i] + off + col1 + d1, col2 - d2);
364
33.3k
  }
365
3.15k
366
3.15k
  n += row1 - r1;
367
36.4k
  for (i = 0; i < row2 - r2; 
++i33.3k
) {
368
33.3k
    isl_seq_cpy(prod->row[n + i], mat2->row[r2 + i], off);
369
33.3k
    isl_seq_clr(prod->row[n + i] + off, d1);
370
33.3k
    isl_seq_cpy(prod->row[n + i] + off + d1,
371
33.3k
          mat2->row[r2 + i] + off, d2);
372
33.3k
    isl_seq_clr(prod->row[n + i] + off + d1 + d2, col1 - d1);
373
33.3k
    isl_seq_cpy(prod->row[n + i] + off + col1 + d1,
374
33.3k
          mat2->row[r2 + i] + off + d2, col2 - d2);
375
33.3k
  }
376
3.15k
377
3.15k
  return prod;
378
3.15k
}
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
58.6k
{
386
58.6k
  if (var->index == -1)
387
244
    return;
388
58.4k
  if (var->is_row && 
var->index >= r141.1k
)
389
33.3k
    var->index += r2;
390
58.4k
  if (!var->is_row && 
var->index >= d117.2k
)
391
16.0k
    var->index += d2;
392
58.4k
}
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
58.6k
{
401
58.6k
  if (var->index == -1)
402
244
    return;
403
58.4k
  if (var->is_row) {
404
41.1k
    if (var->index < r2)
405
7.78k
      var->index += r1;
406
33.3k
    else
407
33.3k
      var->index += row1;
408
41.1k
  } else {
409
17.2k
    if (var->index < d2)
410
1.20k
      var->index += d1;
411
16.0k
    else
412
16.0k
      var->index += col1;
413
17.2k
  }
414
58.4k
}
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
3.15k
{
436
3.15k
  int i;
437
3.15k
  struct isl_tab *prod;
438
3.15k
  unsigned off;
439
3.15k
  unsigned r1, r2, d1, d2;
440
3.15k
441
3.15k
  if (!tab1 || !tab2)
442
0
    return NULL;
443
3.15k
444
3.15k
  isl_assert(tab1->mat->ctx, tab1->M == tab2->M, return NULL);
445
3.15k
  isl_assert(tab1->mat->ctx, tab1->rational == tab2->rational, return NULL);
446
3.15k
  isl_assert(tab1->mat->ctx, tab1->cone == tab2->cone, return NULL);
447
3.15k
  isl_assert(tab1->mat->ctx, !tab1->row_sign, return NULL);
448
3.15k
  isl_assert(tab1->mat->ctx, !tab2->row_sign, return NULL);
449
3.15k
  isl_assert(tab1->mat->ctx, tab1->n_param == 0, return NULL);
450
3.15k
  isl_assert(tab1->mat->ctx, tab2->n_param == 0, return NULL);
451
3.15k
  isl_assert(tab1->mat->ctx, tab1->n_div == 0, return NULL);
452
3.15k
  isl_assert(tab1->mat->ctx, tab2->n_div == 0, return NULL);
453
3.15k
454
3.15k
  off = 2 + tab1->M;
455
3.15k
  r1 = tab1->n_redundant;
456
3.15k
  r2 = tab2->n_redundant;
457
3.15k
  d1 = tab1->n_dead;
458
3.15k
  d2 = tab2->n_dead;
459
3.15k
  prod = isl_calloc_type(tab1->mat->ctx, struct isl_tab);
460
3.15k
  if (!prod)
461
0
    return NULL;
462
3.15k
  prod->mat = tab_mat_product(tab1->mat, tab2->mat,
463
3.15k
        tab1->n_row, tab2->n_row,
464
3.15k
        tab1->n_col, tab2->n_col, off, r1, r2, d1, d2);
465
3.15k
  if (!prod->mat)
466
0
    goto error;
467
3.15k
  prod->var = isl_alloc_array(tab1->mat->ctx, struct isl_tab_var,
468
3.15k
          tab1->max_var + tab2->max_var);
469
3.15k
  if ((tab1->max_var + tab2->max_var) && !prod->var)
470
0
    goto error;
471
20.6k
  
for (i = 0; 3.15k
i < tab1->n_var;
++i17.5k
) {
472
17.5k
    prod->var[i] = tab1->var[i];
473
17.5k
    update_index1(&prod->var[i], r1, r2, d1, d2);
474
17.5k
  }
475
20.6k
  for (i = 0; i < tab2->n_var; 
++i17.5k
) {
476
17.5k
    prod->var[tab1->n_var + i] = tab2->var[i];
477
17.5k
    update_index2(&prod->var[tab1->n_var + i],
478
17.5k
        tab1->n_row, tab1->n_col,
479
17.5k
        r1, r2, d1, d2);
480
17.5k
  }
481
3.15k
  prod->con = isl_alloc_array(tab1->mat->ctx, struct isl_tab_var,
482
3.15k
          tab1->max_con +  tab2->max_con);
483
3.15k
  if ((tab1->max_con + tab2->max_con) && !prod->con)
484
0
    goto error;
485
44.2k
  
for (i = 0; 3.15k
i < tab1->n_con;
++i41.1k
) {
486
41.1k
    prod->con[i] = tab1->con[i];
487
41.1k
    update_index1(&prod->con[i], r1, r2, d1, d2);
488
41.1k
  }
489
44.2k
  for (i = 0; i < tab2->n_con; 
++i41.1k
) {
490
41.1k
    prod->con[tab1->n_con + i] = tab2->con[i];
491
41.1k
    update_index2(&prod->con[tab1->n_con + i],
492
41.1k
        tab1->n_row, tab1->n_col,
493
41.1k
        r1, r2, d1, d2);
494
41.1k
  }
495
3.15k
  prod->col_var = isl_alloc_array(tab1->mat->ctx, int,
496
3.15k
          tab1->n_col + tab2->n_col);
497
3.15k
  if ((tab1->n_col + tab2->n_col) && !prod->col_var)
498
0
    goto error;
499
20.4k
  
for (i = 0; 3.15k
i < tab1->n_col;
++i17.2k
) {
500
17.2k
    int pos = i < d1 ? 
i1.20k
:
i + d216.0k
;
501
17.2k
    prod->col_var[pos] = tab1->col_var[i];
502
17.2k
  }
503
20.4k
  for (i = 0; i < tab2->n_col; 
++i17.2k
) {
504
17.2k
    int pos = i < d2 ? 
d1 + i1.20k
:
tab1->n_col + i16.0k
;
505
17.2k
    int t = tab2->col_var[i];
506
17.2k
    if (t >= 0)
507
239
      t += tab1->n_var;
508
17.0k
    else
509
17.0k
      t -= tab1->n_con;
510
17.2k
    prod->col_var[pos] = t;
511
17.2k
  }
512
3.15k
  prod->row_var = isl_alloc_array(tab1->mat->ctx, int,
513
3.15k
          tab1->mat->n_row + tab2->mat->n_row);
514
3.15k
  if ((tab1->mat->n_row + tab2->mat->n_row) && !prod->row_var)
515
0
    goto error;
516
44.2k
  
for (i = 0; 3.15k
i < tab1->n_row;
++i41.1k
) {
517
41.1k
    int pos = i < r1 ? 
i7.78k
:
i + r233.3k
;
518
41.1k
    prod->row_var[pos] = tab1->row_var[i];
519
41.1k
  }
520
44.2k
  for (i = 0; i < tab2->n_row; 
++i41.1k
) {
521
41.1k
    int pos = i < r2 ? 
r1 + i7.78k
:
tab1->n_row + i33.3k
;
522
41.1k
    int t = tab2->row_var[i];
523
41.1k
    if (t >= 0)
524
17.2k
      t += tab1->n_var;
525
23.8k
    else
526
23.8k
      t -= tab1->n_con;
527
41.1k
    prod->row_var[pos] = t;
528
41.1k
  }
529
3.15k
  prod->samples = NULL;
530
3.15k
  prod->sample_index = NULL;
531
3.15k
  prod->n_row = tab1->n_row + tab2->n_row;
532
3.15k
  prod->n_con = tab1->n_con + tab2->n_con;
533
3.15k
  prod->n_eq = 0;
534
3.15k
  prod->max_con = tab1->max_con + tab2->max_con;
535
3.15k
  prod->n_col = tab1->n_col + tab2->n_col;
536
3.15k
  prod->n_var = tab1->n_var + tab2->n_var;
537
3.15k
  prod->max_var = tab1->max_var + tab2->max_var;
538
3.15k
  prod->n_param = 0;
539
3.15k
  prod->n_div = 0;
540
3.15k
  prod->n_dead = tab1->n_dead + tab2->n_dead;
541
3.15k
  prod->n_redundant = tab1->n_redundant + tab2->n_redundant;
542
3.15k
  prod->rational = tab1->rational;
543
3.15k
  prod->empty = tab1->empty || tab2->empty;
544
3.15k
  prod->strict_redundant = tab1->strict_redundant || tab2->strict_redundant;
545
3.15k
  prod->need_undo = 0;
546
3.15k
  prod->in_undo = 0;
547
3.15k
  prod->M = tab1->M;
548
3.15k
  prod->cone = tab1->cone;
549
3.15k
  prod->bottom.type = isl_tab_undo_bottom;
550
3.15k
  prod->bottom.next = NULL;
551
3.15k
  prod->top = &prod->bottom;
552
3.15k
553
3.15k
  prod->n_zero = 0;
554
3.15k
  prod->n_unbounded = 0;
555
3.15k
  prod->basis = NULL;
556
3.15k
557
3.15k
  return prod;
558
0
error:
559
0
  isl_tab_free(prod);
560
0
  return NULL;
561
3.15k
}
562
563
static struct isl_tab_var *var_from_index(struct isl_tab *tab, int i)
564
117M
{
565
117M
  if (i >= 0)
566
33.3M
    return &tab->var[i];
567
83.9M
  else
568
83.9M
    return &tab->con[~i];
569
117M
}
570
571
struct isl_tab_var *isl_tab_var_from_row(struct isl_tab *tab, int i)
572
87.8M
{
573
87.8M
  return var_from_index(tab, tab->row_var[i]);
574
87.8M
}
575
576
static struct isl_tab_var *var_from_col(struct isl_tab *tab, int i)
577
26.9M
{
578
26.9M
  return var_from_index(tab, tab->col_var[i]);
579
26.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
2.16M
{
588
2.16M
  int i;
589
2.16M
  unsigned off = 2 + tab->M;
590
2.16M
591
2.16M
  if (var->is_row)
592
1.46M
    return 0;
593
3.90M
  
for (i = tab->n_redundant; 694k
i < tab->n_row;
++i3.21M
) {
594
3.42M
    if (!isl_int_is_neg(tab->mat->row[i][off + var->index]))
595
3.42M
      
continue2.95M
;
596
463k
    if (isl_tab_var_from_row(tab, i)->is_nonneg)
597
207k
      return 0;
598
463k
  }
599
694k
  
return 1486k
;
600
694k
}
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
1.79M
{
609
1.79M
  int i;
610
1.79M
  unsigned off = 2 + tab->M;
611
1.79M
612
1.79M
  if (var->is_row)
613
849k
    return 0;
614
7.18M
  
for (i = tab->n_redundant; 944k
i < tab->n_row;
++i6.24M
) {
615
6.74M
    if (!isl_int_is_pos(tab->mat->row[i][off + var->index]))
616
6.74M
      
continue5.82M
;
617
918k
    if (isl_tab_var_from_row(tab, i)->is_nonneg)
618
499k
      return 0;
619
918k
  }
620
944k
  
return 1445k
;
621
944k
}
622
623
static int row_cmp(struct isl_tab *tab, int r1, int r2, int c, isl_int *t)
624
1.40M
{
625
1.40M
  unsigned off = 2 + tab->M;
626
1.40M
627
1.40M
  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
1.40M
  }
635
1.40M
  isl_int_mul(*t, tab->mat->row[r1][1], tab->mat->row[r2][off + c]);
636
1.40M
  isl_int_submul(*t, tab->mat->row[r2][1], tab->mat->row[r1][off + c]);
637
1.40M
  return isl_int_sgn(*t);
638
1.40M
}
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
3.38M
{
663
3.38M
  int j, r, tsgn;
664
3.38M
  isl_int t;
665
3.38M
  unsigned off = 2 + tab->M;
666
3.38M
667
3.38M
  isl_int_init(t);
668
3.38M
  r = -1;
669
42.9M
  for (j = tab->n_redundant; j < tab->n_row; 
++j39.5M
) {
670
39.5M
    if (var && 
j == var->index33.7M
)
671
2.93M
      continue;
672
36.5M
    if (!isl_tab_var_from_row(tab, j)->is_nonneg)
673
8.43M
      continue;
674
28.1M
    if (sgn * isl_int_sgn(tab->mat->row[j][off + c]) >= 0)
675
24.4M
      continue;
676
3.65M
    if (r < 0) {
677
2.25M
      r = j;
678
2.25M
      continue;
679
2.25M
    }
680
1.40M
    tsgn = sgn * row_cmp(tab, r, j, c, &t);
681
1.40M
    if (tsgn < 0 || 
(1.03M
tsgn == 01.03M
&&
682
1.03M
              
tab->row_var[j] < tab->row_var[r]421k
))
683
681k
      r = j;
684
1.40M
  }
685
3.38M
  isl_int_clear(t);
686
3.38M
  return r;
687
3.38M
}
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
4.25M
{
709
4.25M
  int j, r, c;
710
4.25M
  isl_int *tr;
711
4.25M
712
4.25M
  *row = *col = -1;
713
4.25M
714
4.25M
  isl_assert(tab->mat->ctx, var->is_row, return);
715
4.25M
  tr = tab->mat->row[var->index] + 2 + tab->M;
716
4.25M
717
4.25M
  c = -1;
718
33.6M
  for (j = tab->n_dead; j < tab->n_col; 
++j29.3M
) {
719
29.3M
    if (isl_int_is_zero(tr[j]))
720
29.3M
      
continue22.9M
;
721
6.41M
    if (isl_int_sgn(tr[j]) != sgn &&
722
6.41M
        
var_from_col(tab, j)->is_nonneg3.52M
)
723
2.34M
      continue;
724
4.06M
    if (c < 0 || 
tab->col_var[j] < tab->col_var[c]1.00M
)
725
3.30M
      c = j;
726
4.06M
  }
727
4.25M
  if (c < 0)
728
1.19M
    return;
729
3.06M
730
3.06M
  sgn *= isl_int_sgn(tr[c]);
731
3.06M
  r = pivot_row(tab, skip_var, sgn, c);
732
3.06M
  *row = r < 0 ? 
var->index1.13M
:
r1.92M
;
733
3.06M
  *col = c;
734
3.06M
}
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
18.2M
{
744
18.2M
  int i;
745
18.2M
  unsigned off = 2 + tab->M;
746
18.2M
747
18.2M
  if (tab->row_var[row] < 0 && 
!isl_tab_var_from_row(tab, row)->is_nonneg13.9M
)
748
908k
    return 0;
749
17.3M
750
17.3M
  if (isl_int_is_neg(tab->mat->row[row][1]))
751
17.3M
    
return 01.55M
;
752
15.7M
  if (tab->strict_redundant && 
isl_int_is_zero43
(tab->mat->row[row][1]))
753
15.7M
    
return 042
;
754
15.7M
  if (tab->M && 
isl_int_is_neg41.8k
(tab->mat->row[row][2]))
755
15.7M
    
return 02.24k
;
756
15.7M
757
65.1M
  
for (i = tab->n_dead; 15.7M
i < tab->n_col;
++i49.4M
) {
758
63.4M
    if (isl_int_is_zero(tab->mat->row[row][off + i]))
759
63.4M
      
continue44.4M
;
760
19.0M
    if (tab->col_var[i] >= 0)
761
7.24M
      return 0;
762
11.7M
    if (isl_int_is_neg(tab->mat->row[row][off + i]))
763
11.7M
      
return 06.61M
;
764
5.17M
    if (!var_from_col(tab, i)->is_nonneg)
765
193k
      return 0;
766
5.17M
  }
767
15.7M
  
return 11.72M
;
768
15.7M
}
769
770
static void swap_rows(struct isl_tab *tab, int row1, int row2)
771
1.61M
{
772
1.61M
  int t;
773
1.61M
  enum isl_tab_row_sign s;
774
1.61M
775
1.61M
  t = tab->row_var[row1];
776
1.61M
  tab->row_var[row1] = tab->row_var[row2];
777
1.61M
  tab->row_var[row2] = t;
778
1.61M
  isl_tab_var_from_row(tab, row1)->index = row1;
779
1.61M
  isl_tab_var_from_row(tab, row2)->index = row2;
780
1.61M
  tab->mat = isl_mat_swap_rows(tab->mat, row1, row2);
781
1.61M
782
1.61M
  if (!tab->row_sign)
783
1.60M
    return;
784
9.39k
  s = tab->row_sign[row1];
785
9.39k
  tab->row_sign[row1] = tab->row_sign[row2];
786
9.39k
  tab->row_sign[row2] = s;
787
9.39k
}
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
12.7M
{
802
12.7M
  struct isl_tab_undo *undo;
803
12.7M
804
12.7M
  if (!tab)
805
0
    return isl_stat_error;
806
12.7M
  if (!tab->need_undo)
807
9.44M
    return isl_stat_ok;
808
3.34M
809
3.34M
  undo = isl_alloc_type(tab->mat->ctx, struct isl_tab_undo);
810
3.34M
  if (!undo)
811
0
    goto error;
812
3.34M
  undo->type = type;
813
3.34M
  undo->u = u;
814
3.34M
  undo->next = tab->top;
815
3.34M
  tab->top = undo;
816
3.34M
817
3.34M
  return isl_stat_ok;
818
0
error:
819
0
  free_undo(tab);
820
0
  tab->top = NULL;
821
0
  return isl_stat_error;
822
3.34M
}
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
12.3M
{
827
12.3M
  union isl_tab_undo_val u;
828
12.3M
  if (var->is_row)
829
12.0M
    u.var_index = tab->row_var[var->index];
830
261k
  else
831
261k
    u.var_index = tab->col_var[var->index];
832
12.3M
  return push_union(tab, type, u);
833
12.3M
}
834
835
isl_stat isl_tab_push(struct isl_tab *tab, enum isl_tab_undo_type type)
836
413k
{
837
413k
  union isl_tab_undo_val u = { 0 };
838
413k
  return push_union(tab, type, u);
839
413k
}
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
758
{
846
758
  int i;
847
758
  union isl_tab_undo_val u;
848
758
849
758
  u.col_var = isl_alloc_array(tab->mat->ctx, int, tab->n_col);
850
758
  if (tab->n_col && !u.col_var)
851
0
    return isl_stat_error;
852
8.59k
  
for (i = 0; 758
i < tab->n_col;
++i7.83k
)
853
7.83k
    u.col_var[i] = tab->col_var[i];
854
758
  return push_union(tab, isl_tab_undo_saved_basis, u);
855
758
}
856
857
isl_stat isl_tab_push_callback(struct isl_tab *tab,
858
  struct isl_tab_callback *callback)
859
21.4k
{
860
21.4k
  union isl_tab_undo_val u;
861
21.4k
  u.callback = callback;
862
21.4k
  return push_union(tab, isl_tab_undo_callback, u);
863
21.4k
}
864
865
struct isl_tab *isl_tab_init_samples(struct isl_tab *tab)
866
7.60k
{
867
7.60k
  if (!tab)
868
0
    return NULL;
869
7.60k
870
7.60k
  tab->n_sample = 0;
871
7.60k
  tab->n_outside = 0;
872
7.60k
  tab->samples = isl_mat_alloc(tab->mat->ctx, 1, 1 + tab->n_var);
873
7.60k
  if (!tab->samples)
874
0
    goto error;
875
7.60k
  tab->sample_index = isl_alloc_array(tab->mat->ctx, int, 1);
876
7.60k
  if (!tab->sample_index)
877
0
    goto error;
878
7.60k
  return tab;
879
0
error:
880
0
  isl_tab_free(tab);
881
0
  return NULL;
882
7.60k
}
883
884
int isl_tab_add_sample(struct isl_tab *tab, __isl_take isl_vec *sample)
885
11.4k
{
886
11.4k
  if (!tab || !sample)
887
0
    goto error;
888
11.4k
889
11.4k
  if (tab->n_sample + 1 > tab->samples->n_row) {
890
3.72k
    int *t = isl_realloc_array(tab->mat->ctx,
891
3.72k
          tab->sample_index, int, tab->n_sample + 1);
892
3.72k
    if (!t)
893
0
      goto error;
894
3.72k
    tab->sample_index = t;
895
3.72k
  }
896
11.4k
897
11.4k
  tab->samples = isl_mat_extend(tab->samples,
898
11.4k
        tab->n_sample + 1, tab->samples->n_col);
899
11.4k
  if (!tab->samples)
900
0
    goto error;
901
11.4k
902
11.4k
  isl_seq_cpy(tab->samples->row[tab->n_sample], sample->el, sample->size);
903
11.4k
  isl_vec_free(sample);
904
11.4k
  tab->sample_index[tab->n_sample] = tab->n_sample;
905
11.4k
  tab->n_sample++;
906
11.4k
907
11.4k
  return 0;
908
0
error:
909
0
  isl_vec_free(sample);
910
0
  return -1;
911
11.4k
}
912
913
struct isl_tab *isl_tab_drop_sample(struct isl_tab *tab, int s)
914
6.06k
{
915
6.06k
  if (s != tab->n_outside) {
916
3.75k
    int t = tab->sample_index[tab->n_outside];
917
3.75k
    tab->sample_index[tab->n_outside] = tab->sample_index[s];
918
3.75k
    tab->sample_index[s] = t;
919
3.75k
    isl_mat_swap_rows(tab->samples, tab->n_outside, s);
920
3.75k
  }
921
6.06k
  tab->n_outside++;
922
6.06k
  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.06k
927
6.06k
  return tab;
928
6.06k
}
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.3k
{
935
28.3k
  union isl_tab_undo_val u;
936
28.3k
937
28.3k
  if (!tab)
938
0
    return isl_stat_error;
939
28.3k
940
28.3k
  u.n = tab->n_sample;
941
28.3k
  return push_union(tab, isl_tab_undo_saved_samples, u);
942
28.3k
}
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
2.24M
{
958
2.24M
  struct isl_tab_var *var = isl_tab_var_from_row(tab, row);
959
2.24M
  var->is_redundant = 1;
960
2.24M
  isl_assert(tab->mat->ctx, row >= tab->n_redundant, return -1);
961
2.24M
  if (tab->preserve || 
tab->need_undo1.57M
||
tab->row_var[row] >= 01.44M
) {
962
1.60M
    if (tab->row_var[row] >= 0 && 
!var->is_nonneg1.23M
) {
963
1.22M
      var->is_nonneg = 1;
964
1.22M
      if (isl_tab_push_var(tab, isl_tab_undo_nonneg, var) < 0)
965
0
        return -1;
966
1.60M
    }
967
1.60M
    if (row != tab->n_redundant)
968
1.10M
      swap_rows(tab, row, tab->n_redundant);
969
1.60M
    tab->n_redundant++;
970
1.60M
    return isl_tab_push_var(tab, isl_tab_undo_redundant, var);
971
1.60M
  } else {
972
646k
    if (row != tab->n_row - 1)
973
392k
      swap_rows(tab, row, tab->n_row - 1);
974
646k
    isl_tab_var_from_row(tab, tab->n_row - 1)->index = -1;
975
646k
    tab->n_row--;
976
646k
    return 1;
977
646k
  }
978
2.24M
}
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.49k
{
987
9.49k
  if (!tab)
988
0
    return -1;
989
9.49k
  if (!tab->rational && 
tab->need_undo9.46k
)
990
9.46k
    if (isl_tab_push(tab, isl_tab_undo_rational) < 0)
991
0
      return -1;
992
9.49k
  tab->rational = 1;
993
9.49k
  return 0;
994
9.49k
}
995
996
isl_stat isl_tab_mark_empty(struct isl_tab *tab)
997
73.9k
{
998
73.9k
  if (!tab)
999
0
    return isl_stat_error;
1000
73.9k
  if (!tab->empty && 
tab->need_undo73.0k
)
1001
62.3k
    if (isl_tab_push(tab, isl_tab_undo_empty) < 0)
1002
0
      return isl_stat_error;
1003
73.9k
  tab->empty = 1;
1004
73.9k
  return isl_stat_ok;
1005
73.9k
}
1006
1007
int isl_tab_freeze_constraint(struct isl_tab *tab, int con)
1008
372k
{
1009
372k
  struct isl_tab_var *var;
1010
372k
1011
372k
  if (!tab)
1012
0
    return -1;
1013
372k
1014
372k
  var = &tab->con[con];
1015
372k
  if (var->frozen)
1016
0
    return 0;
1017
372k
  if (var->index < 0)
1018
25.0k
    return 0;
1019
347k
  var->frozen = 1;
1020
347k
1021
347k
  if (tab->need_undo)
1022
314k
    return isl_tab_push_var(tab, isl_tab_undo_freeze, var);
1023
33.2k
1024
33.2k
  return 0;
1025
33.2k
}
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
3.37M
{
1046
3.37M
  int i;
1047
3.37M
  struct isl_mat *mat = tab->mat;
1048
3.37M
  unsigned off = 2 + tab->M;
1049
3.37M
1050
3.37M
  if (!tab->row_sign)
1051
3.34M
    return;
1052
28.3k
1053
28.3k
  if (tab->row_sign[row] == 0)
1054
22.3k
    return;
1055
6.01k
  isl_assert(mat->ctx, row_sgn > 0, return);
1056
6.01k
  isl_assert(mat->ctx, tab->row_sign[row] == isl_tab_row_neg, return);
1057
6.01k
  tab->row_sign[row] = isl_tab_row_pos;
1058
54.4k
  for (i = 0; i < tab->n_row; 
++i48.4k
) {
1059
48.4k
    int s;
1060
48.4k
    if (i == row)
1061
6.01k
      continue;
1062
42.4k
    s = isl_int_sgn(mat->row[i][off + col]);
1063
42.4k
    if (!s)
1064
24.8k
      continue;
1065
17.6k
    if (!tab->row_sign[i])
1066
6.79k
      continue;
1067
10.8k
    if (s < 0 && 
tab->row_sign[i] == isl_tab_row_neg6.01k
)
1068
0
      continue;
1069
10.8k
    if (s > 0 && 
tab->row_sign[i] == isl_tab_row_pos4.79k
)
1070
4.79k
      continue;
1071
6.01k
    tab->row_sign[i] = isl_tab_row_unknown;
1072
6.01k
  }
1073
6.01k
}
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
3.37M
{
1127
3.37M
  int i, j;
1128
3.37M
  int sgn;
1129
3.37M
  int t;
1130
3.37M
  isl_ctx *ctx;
1131
3.37M
  struct isl_mat *mat = tab->mat;
1132
3.37M
  struct isl_tab_var *var;
1133
3.37M
  unsigned off = 2 + tab->M;
1134
3.37M
1135
3.37M
  ctx = isl_tab_get_ctx(tab);
1136
3.37M
  if (isl_ctx_next_operation(ctx) < 0)
1137
0
    return -1;
1138
3.37M
1139
3.37M
  isl_int_swap(mat->row[row][0], mat->row[row][off + col]);
1140
3.37M
  sgn = isl_int_sgn(mat->row[row][0]);
1141
3.37M
  if (sgn < 0) {
1142
1.93M
    isl_int_neg(mat->row[row][0], mat->row[row][0]);
1143
1.93M
    isl_int_neg(mat->row[row][off + col], mat->row[row][off + col]);
1144
1.93M
  } else
1145
12.0M
    
for (j = 0; 1.43M
j < off - 1 + tab->n_col;
++j10.6M
) {
1146
10.6M
      if (j == off - 1 + col)
1147
1.43M
        continue;
1148
9.22M
      isl_int_neg(mat->row[row][1 + j], mat->row[row][1 + j]);
1149
9.22M
    }
1150
3.37M
  if (!isl_int_is_one(mat->row[row][0]))
1151
3.37M
    
isl_seq_normalize(mat->ctx, mat->row[row], off + tab->n_col)603k
;
1152
40.8M
  for (i = 0; i < tab->n_row; 
++i37.5M
) {
1153
37.5M
    if (i == row)
1154
3.37M
      continue;
1155
34.1M
    if (isl_int_is_zero(mat->row[i][off + col]))
1156
34.1M
      
continue25.2M
;
1157
8.86M
    isl_int_mul(mat->row[i][0], mat->row[i][0], mat->row[row][0]);
1158
95.3M
    for (j = 0; j < off - 1 + tab->n_col; 
++j86.4M
) {
1159
86.4M
      if (j == off - 1 + col)
1160
8.86M
        continue;
1161
77.6M
      isl_int_mul(mat->row[i][1 + j],
1162
77.6M
            mat->row[i][1 + j], mat->row[row][0]);
1163
77.6M
      isl_int_addmul(mat->row[i][1 + j],
1164
77.6M
            mat->row[i][off + col], mat->row[row][1 + j]);
1165
77.6M
    }
1166
8.86M
    isl_int_mul(mat->row[i][off + col],
1167
8.86M
          mat->row[i][off + col], mat->row[row][off + col]);
1168
8.86M
    if (!isl_int_is_one(mat->row[i][0]))
1169
8.86M
      
isl_seq_normalize(mat->ctx, mat->row[i], off + tab->n_col)4.17M
;
1170
8.86M
  }
1171
3.37M
  t = tab->row_var[row];
1172
3.37M
  tab->row_var[row] = tab->col_var[col];
1173
3.37M
  tab->col_var[col] = t;
1174
3.37M
  var = isl_tab_var_from_row(tab, row);
1175
3.37M
  var->is_row = 1;
1176
3.37M
  var->index = row;
1177
3.37M
  var = var_from_col(tab, col);
1178
3.37M
  var->is_row = 0;
1179
3.37M
  var->index = col;
1180
3.37M
  update_row_sign(tab, row, col, sgn);
1181
3.37M
  if (tab->in_undo)
1182
116k
    return 0;
1183
32.3M
  
for (i = tab->n_redundant; 3.25M
i < tab->n_row;
++i29.0M
) {
1184
29.0M
    if (isl_int_is_zero(mat->row[i][off + col]))
1185
29.0M
      
continue18.2M
;
1186
10.8M
    if (!isl_tab_var_from_row(tab, i)->frozen &&
1187
10.8M
        
isl_tab_row_is_redundant(tab, i)10.5M
) {
1188
1.61M
      int redo = isl_tab_mark_redundant(tab, i);
1189
1.61M
      if (redo < 0)
1190
0
        return -1;
1191
1.61M
      if (redo)
1192
102k
        --i;
1193
1.61M
    }
1194
10.8M
  }
1195
3.25M
  return 0;
1196
3.25M
}
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
1.71M
{
1207
1.71M
  int r;
1208
1.71M
  unsigned off = 2 + tab->M;
1209
1.71M
1210
1.71M
  if (var->is_row)
1211
1.47M
    return 0;
1212
239k
1213
239k
  if (sign == 0) {
1214
55.8k
    for (r = tab->n_redundant; r < tab->n_row; 
++r35.4k
)
1215
55.8k
      if (!isl_int_is_zero(tab->mat->row[r][off+var->index]))
1216
55.8k
        
break20.4k
;
1217
20.4k
    isl_assert(tab->mat->ctx, r < tab->n_row, return -1);
1218
219k
  } else {
1219
219k
    r = pivot_row(tab, NULL, sign, var->index);
1220
219k
    isl_assert(tab->mat->ctx, r >= 0, return -1);
1221
219k
  }
1222
239k
1223
239k
  return isl_tab_pivot(tab, r, var->index);
1224
239k
}
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
1.29M
{
1266
1.29M
  int row, col;
1267
1.29M
1268
1.29M
  if (max_is_manifestly_unbounded(tab, var))
1269
118k
    return 1;
1270
1.17M
  if (to_row(tab, var, 1) < 0)
1271
0
    return -2;
1272
2.01M
  
while (1.17M
!isl_int_is_pos(tab->mat->row[var->index][1])) {
1273
1.57M
    find_pivot(tab, var, var, 1, &row, &col);
1274
1.57M
    if (row == -1)
1275
515k
      return isl_int_sgn(tab->mat->row[var->index][1]);
1276
1.06M
    if (isl_tab_pivot(tab, row, col) < 0)
1277
0
      return -2;
1278
1.06M
    if (!var->is_row) /* manifestly unbounded */
1279
216k
      return 1;
1280
1.06M
  }
1281
1.17M
  
return 1442k
;
1282
1.17M
}
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
5.16M
{
1300
5.16M
  if (!tab->M)
1301
5.16M
    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
4.54M
{
1311
4.54M
  if (!tab->M)
1312
4.54M
    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
4.93M
{
1326
4.93M
  int row, col;
1327
4.93M
1328
5.16M
  while (row_is_neg(tab, var->index)) {
1329
690k
    find_pivot(tab, var, var, 1, &row, &col);
1330
690k
    if (row == -1)
1331
70.6k
      break;
1332
619k
    if (isl_tab_pivot(tab, row, col) < 0)
1333
0
      return -2;
1334
619k
    if (!var->is_row) /* manifestly unbounded */
1335
383k
      return 1;
1336
619k
  }
1337
4.93M
  
return row_sgn(tab, var->index)4.54M
;
1338
4.93M
}
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
156k
{
1347
156k
  int row, col;
1348
156k
1349
175k
  while (isl_int_is_neg(tab->mat->row[var->index][1])) {
1350
167k
    find_pivot(tab, var, var, 1, &row, &col);
1351
167k
    if (row == -1)
1352
79.9k
      break;
1353
87.3k
    if (row == var->index) /* manifestly unbounded */
1354
67.9k
      return 1;
1355
19.4k
    if (isl_tab_pivot(tab, row, col) < 0)
1356
0
      return -1;
1357
19.4k
  }
1358
156k
  
return !88.5k
isl_int_is_neg88.5k
(tab->mat->row[var->index][1]);
1359
156k
}
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.4k
{
1380
70.4k
  int row, col;
1381
70.4k
  struct isl_tab_var *pivot_var = NULL;
1382
70.4k
1383
70.4k
  if (min_is_manifestly_unbounded(tab, var))
1384
0
    return -1;
1385
70.4k
  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.4k
  if (var->is_redundant)
1407
0
    return 0;
1408
106k
  
while (69.4k
!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.1k
      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.4k
  
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
290k
{
1434
290k
  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
290k
  }
1440
290k
  return isl_int_is_neg(tab->mat->row[row][1]) &&
1441
290k
         
isl_int_abs_ge166k
(tab->mat->row[row][1],
1442
290k
            tab->mat->row[row][0]);
1443
290k
}
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
882k
{
1456
882k
  int row, col;
1457
882k
  struct isl_tab_var *pivot_var;
1458
882k
1459
882k
  if (min_is_manifestly_unbounded(tab, var))
1460
111
    return 1;
1461
881k
  if (!var->is_row) {
1462
103k
    col = var->index;
1463
103k
    row = pivot_row(tab, NULL, -1, col);
1464
103k
    pivot_var = var_from_col(tab, col);
1465
103k
    if (isl_tab_pivot(tab, row, col) < 0)
1466
0
      return -1;
1467
103k
    if (var->is_redundant)
1468
10.2k
      return 0;
1469
93.4k
    if (row_at_most_neg_one(tab, var->index)) {
1470
73.6k
      if (var->is_nonneg) {
1471
73.6k
        if (!pivot_var->is_redundant &&
1472
73.6k
            pivot_var->index == row) {
1473
68.0k
          if (isl_tab_pivot(tab, row, col) < 0)
1474
0
            return -1;
1475
5.62k
        } else
1476
5.62k
          if (restore_row(tab, var) < -1)
1477
0
            return -1;
1478
73.6k
      }
1479
73.6k
      return 1;
1480
73.6k
    }
1481
93.4k
  }
1482
798k
  if (var->is_redundant)
1483
11.7k
    return 0;
1484
908k
  
do 786k
{
1485
908k
    find_pivot(tab, var, var, -1, &row, &col);
1486
908k
    if (row == var->index) {
1487
439k
      if (var->is_nonneg && 
restore_row(tab, var) < -1398k
)
1488
0
        return -1;
1489
439k
      return 1;
1490
439k
    }
1491
469k
    if (row == -1)
1492
184k
      return 0;
1493
284k
    pivot_var = var_from_col(tab, col);
1494
284k
    if (isl_tab_pivot(tab, row, col) < 0)
1495
0
      return -1;
1496
284k
    if (var->is_redundant)
1497
87.9k
      return 0;
1498
196k
  } while (!row_at_most_neg_one(tab, var->index));
1499
786k
  
if (73.9k
var->is_nonneg73.9k
) {
1500
62.0k
    /* pivot back to non-negative value */
1501
62.0k
    if (!pivot_var->is_redundant && pivot_var->index == row)
1502
57.9k
      if (isl_tab_pivot(tab, row, col) < 0)
1503
0
        return -1;
1504
62.0k
    if (restore_row(tab, var) < -1)
1505
0
      return -1;
1506
73.9k
  }
1507
73.9k
  return 1;
1508
73.9k
}
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
751k
{
1515
751k
  int row, col;
1516
751k
  isl_int *r;
1517
751k
1518
751k
  if (max_is_manifestly_unbounded(tab, var))
1519
334k
    return 1;
1520
417k
  if (to_row(tab, var, 1) < 0)
1521
0
    return -1;
1522
417k
  r = tab->mat->row[var->index];
1523
436k
  while (isl_int_lt(r[1], r[0])) {
1524
19.6k
    find_pivot(tab, var, var, 1, &row, &col);
1525
19.6k
    if (row == -1)
1526
539
      return isl_int_ge(r[1], r[0]);
1527
19.1k
    if (row == var->index) /* manifestly unbounded */
1528
154
      return 1;
1529
18.9k
    if (isl_tab_pivot(tab, row, col) < 0)
1530
0
      return -1;
1531
18.9k
  }
1532
417k
  
return 1416k
;
1533
417k
}
1534
1535
static void swap_cols(struct isl_tab *tab, int col1, int col2)
1536
673k
{
1537
673k
  int t;
1538
673k
  unsigned off = 2 + tab->M;
1539
673k
  t = tab->col_var[col1];
1540
673k
  tab->col_var[col1] = tab->col_var[col2];
1541
673k
  tab->col_var[col2] = t;
1542
673k
  var_from_col(tab, col1)->index = col1;
1543
673k
  var_from_col(tab, col2)->index = col2;
1544
673k
  tab->mat = isl_mat_swap_cols(tab->mat, off + col1, off + col2);
1545
673k
}
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
940k
{
1561
940k
  var_from_col(tab, col)->is_zero = 1;
1562
940k
  if (tab->need_undo) {
1563
114k
    if (isl_tab_push_var(tab, isl_tab_undo_zero,
1564
114k
              var_from_col(tab, col)) < 0)
1565
0
      return -1;
1566
114k
    if (col != tab->n_dead)
1567
43.8k
      swap_cols(tab, col, tab->n_dead);
1568
114k
    tab->n_dead++;
1569
114k
    return 0;
1570
826k
  } else {
1571
826k
    if (col != tab->n_col - 1)
1572
628k
      swap_cols(tab, col, tab->n_col - 1);
1573
826k
    var_from_col(tab, tab->n_col - 1)->index = -1;
1574
826k
    tab->n_col--;
1575
826k
    return 1;
1576
826k
  }
1577
940k
}
1578
1579
static int row_is_manifestly_non_integral(struct isl_tab *tab, int row)
1580
2.93M
{
1581
2.93M
  unsigned off = 2 + tab->M;
1582
2.93M
1583
2.93M
  if (tab->M && 
!0
isl_int_eq0
(tab->mat->row[row][2],
1584
2.93M
          tab->mat->row[row][0]))
1585
2.93M
    
return 00
;
1586
2.93M
  if (isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
1587
2.93M
            tab->n_col - tab->n_dead) != -1)
1588
182k
    return 0;
1589
2.75M
1590
2.75M
  return !isl_int_is_divisible_by(tab->mat->row[row][1],
1591
2.75M
          tab->mat->row[row][0]);
1592
2.75M
}
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
514k
{
1599
514k
  int i;
1600
514k
1601
514k
  if (tab->empty)
1602
1
    return 1;
1603
514k
  if (tab->rational)
1604
11.5k
    return 0;
1605
503k
1606
6.48M
  
for (i = 0; 503k
i < tab->n_var;
++i5.98M
) {
1607
5.98M
    if (!tab->var[i].is_row)
1608
3.04M
      continue;
1609
2.93M
    if (row_is_manifestly_non_integral(tab, tab->var[i].index))
1610
90
      return 1;
1611
2.93M
  }
1612
503k
1613
503k
  
return 0503k
;
1614
503k
}
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
514k
{
1638
514k
  int j;
1639
514k
  struct isl_mat *mat = tab->mat;
1640
514k
  unsigned off = 2 + tab->M;
1641
514k
1642
514k
  if (!var->is_nonneg)
1643
514k
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_internal,
1644
514k
      "expecting non-negative variable",
1645
514k
      return isl_stat_error);
1646
514k
  var->is_zero = 1;
1647
514k
  if (!temp_var && 
tab->need_undo502k
)
1648
454
    if (isl_tab_push_var(tab, isl_tab_undo_zero, var) < 0)
1649
0
      return isl_stat_error;
1650
4.16M
  
for (j = tab->n_dead; 514k
j < tab->n_col;
++j3.65M
) {
1651
3.65M
    int recheck;
1652
3.65M
    if (isl_int_is_zero(mat->row[var->index][off + j]))
1653
3.65M
      
continue3.17M
;
1654
480k
    if (isl_int_is_pos(mat->row[var->index][off + j]))
1655
480k
      
isl_die0
(isl_tab_get_ctx(tab), isl_error_internal,
1656
480k
        "row cannot have positive coefficients",
1657
480k
        return isl_stat_error);
1658
480k
    recheck = isl_tab_kill_col(tab, j);
1659
480k
    if (recheck < 0)
1660
0
      return isl_stat_error;
1661
480k
    if (recheck)
1662
468k
      --j;
1663
480k
  }
1664
514k
  if (!temp_var && 
isl_tab_mark_redundant(tab, var->index) < 0502k
)
1665
0
    return isl_stat_error;
1666
514k
  if (tab_is_manifestly_empty(tab) && 
isl_tab_mark_empty(tab) < 091
)
1667
0
    return isl_stat_error;
1668
514k
  return isl_stat_ok;
1669
514k
}
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
5.13M
{
1679
5.13M
  int r;
1680
5.13M
1681
5.13M
  isl_assert(tab->mat->ctx, tab->n_row < tab->mat->n_row, return -1);
1682
5.13M
  isl_assert(tab->mat->ctx, tab->n_con < tab->max_con, return -1);
1683
5.13M
1684
5.13M
  r = tab->n_con;
1685
5.13M
  tab->con[r].index = tab->n_row;
1686
5.13M
  tab->con[r].is_row = 1;
1687
5.13M
  tab->con[r].is_nonneg = 0;
1688
5.13M
  tab->con[r].is_zero = 0;
1689
5.13M
  tab->con[r].is_redundant = 0;
1690
5.13M
  tab->con[r].frozen = 0;
1691
5.13M
  tab->con[r].negated = 0;
1692
5.13M
  tab->row_var[tab->n_row] = ~r;
1693
5.13M
1694
5.13M
  tab->n_row++;
1695
5.13M
  tab->n_con++;
1696
5.13M
  if (isl_tab_push_var(tab, isl_tab_undo_allocate, &tab->con[r]) < 0)
1697
0
    return -1;
1698
5.13M
1699
5.13M
  return r;
1700
5.13M
}
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
7.62k
{
1709
7.62k
  int i;
1710
7.62k
1711
7.62k
  if (tab->n_var >= tab->max_var)
1712
7.62k
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_internal,
1713
7.62k
      "not enough room for new variable", return -1);
1714
7.62k
  if (first > tab->n_var)
1715
7.62k
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_internal,
1716
7.62k
      "invalid initial position", return -1);
1717
7.62k
1718
8.47k
  
for (i = tab->n_var - 1; 7.62k
i >= first;
--i849
) {
1719
849
    tab->var[i + 1] = tab->var[i];
1720
849
    if (tab->var[i + 1].is_row)
1721
573
      tab->row_var[tab->var[i + 1].index]++;
1722
276
    else
1723
276
      tab->col_var[tab->var[i + 1].index]++;
1724
849
  }
1725
7.62k
1726
7.62k
  tab->n_var++;
1727
7.62k
1728
7.62k
  return 0;
1729
7.62k
}
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
4.98k
{
1738
4.98k
  int i;
1739
4.98k
1740
4.98k
  if (first >= tab->n_var)
1741
4.98k
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_internal,
1742
4.98k
      "invalid initial position", return -1);
1743
4.98k
1744
4.98k
  tab->n_var--;
1745
4.98k
1746
5.57k
  for (i = first; i < tab->n_var; 
++i583
) {
1747
583
    tab->var[i] = tab->var[i + 1];
1748
583
    if (tab->var[i + 1].is_row)
1749
569
      tab->row_var[tab->var[i].index]--;
1750
14
    else
1751
14
      tab->col_var[tab->var[i].index]--;
1752
583
  }
1753
4.98k
1754
4.98k
  return 0;
1755
4.98k
}
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
7.62k
{
1763
7.62k
  int i;
1764
7.62k
  unsigned off = 2 + tab->M;
1765
7.62k
1766
7.62k
  isl_assert(tab->mat->ctx, tab->n_col < tab->mat->n_col, return -1);
1767
7.62k
1768
7.62k
  if (var_insert_entry(tab, r) < 0)
1769
0
    return -1;
1770
7.62k
1771
7.62k
  tab->var[r].index = tab->n_col;
1772
7.62k
  tab->var[r].is_row = 0;
1773
7.62k
  tab->var[r].is_nonneg = 0;
1774
7.62k
  tab->var[r].is_zero = 0;
1775
7.62k
  tab->var[r].is_redundant = 0;
1776
7.62k
  tab->var[r].frozen = 0;
1777
7.62k
  tab->var[r].negated = 0;
1778
7.62k
  tab->col_var[tab->n_col] = r;
1779
7.62k
1780
47.3k
  for (i = 0; i < tab->n_row; 
++i39.6k
)
1781
39.6k
    isl_int_set_si(tab->mat->row[i][off + tab->n_col], 0);
1782
7.62k
1783
7.62k
  tab->n_col++;
1784
7.62k
  if (isl_tab_push_var(tab, isl_tab_undo_allocate, &tab->var[r]) < 0)
1785
0
    return -1;
1786
7.62k
1787
7.62k
  return r;
1788
7.62k
}
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
5.13M
{
1825
5.13M
  int i;
1826
5.13M
  int r;
1827
5.13M
  isl_int *row;
1828
5.13M
  isl_int a, b;
1829
5.13M
  unsigned off = 2 + tab->M;
1830
5.13M
1831
5.13M
  r = isl_tab_allocate_con(tab);
1832
5.13M
  if (r < 0)
1833
0
    return -1;
1834
5.13M
1835
5.13M
  isl_int_init(a);
1836
5.13M
  isl_int_init(b);
1837
5.13M
  row = tab->mat->row[tab->con[r].index];
1838
5.13M
  isl_int_set_si(row[0], 1);
1839
5.13M
  isl_int_set(row[1], line[0]);
1840
5.13M
  isl_seq_clr(row + 2, tab->M + tab->n_col);
1841
48.1M
  for (i = 0; i < tab->n_var; 
++i42.9M
) {
1842
42.9M
    if (tab->var[i].is_zero)
1843
0
      continue;
1844
42.9M
    if (tab->var[i].is_row) {
1845
8.70M
      isl_int_lcm(a,
1846
8.70M
        row[0], tab->mat->row[tab->var[i].index][0]);
1847
8.70M
      isl_int_swap(a, row[0]);
1848
8.70M
      isl_int_divexact(a, row[0], a);
1849
8.70M
      isl_int_divexact(b,
1850
8.70M
        row[0], tab->mat->row[tab->var[i].index][0]);
1851
8.70M
      isl_int_mul(b, b, line[1 + i]);
1852
8.70M
      isl_seq_combine(row + 1, a, row + 1,
1853
8.70M
          b, tab->mat->row[tab->var[i].index] + 1,
1854
8.70M
          1 + tab->M + tab->n_col);
1855
8.70M
    } else
1856
42.9M
      
isl_int_addmul34.2M
(row[off + tab->var[i].index],
1857
42.9M
              line[1 + i], row[0]);
1858
42.9M
    if (tab->M && 
i >= tab->n_param344k
&&
i < tab->n_var - tab->n_div150k
)
1859
42.9M
      
isl_int_submul147k
(row[2], line[1 + i], row[0]);
1860
42.9M
  }
1861
5.13M
  isl_seq_normalize(tab->mat->ctx, row, off + tab->n_col);
1862
5.13M
  isl_int_clear(a);
1863
5.13M
  isl_int_clear(b);
1864
5.13M
1865
5.13M
  if (tab->row_sign)
1866
41.2k
    tab->row_sign[tab->con[r].index] = isl_tab_row_unknown;
1867
5.13M
1868
5.13M
  return r;
1869
5.13M
}
1870
1871
static isl_stat drop_row(struct isl_tab *tab, int row)
1872
1.09M
{
1873
1.09M
  isl_assert(tab->mat->ctx, ~tab->row_var[row] == tab->n_con - 1,
1874
1.09M
    return isl_stat_error);
1875
1.09M
  if (row != tab->n_row - 1)
1876
119k
    swap_rows(tab, row, tab->n_row - 1);
1877
1.09M
  tab->n_row--;
1878
1.09M
  tab->n_con--;
1879
1.09M
  return isl_stat_ok;
1880
1.09M
}
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
4.98k
{
1890
4.98k
  int var;
1891
4.98k
1892
4.98k
  var = tab->col_var[col];
1893
4.98k
  if (col != tab->n_col - 1)
1894
1.07k
    swap_cols(tab, col, tab->n_col - 1);
1895
4.98k
  tab->n_col--;
1896
4.98k
  if (var_drop_entry(tab, var) < 0)
1897
0
    return isl_stat_error;
1898
4.98k
  return isl_stat_ok;
1899
4.98k
}
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
3.88M
{
1909
3.88M
  int r;
1910
3.88M
  int sgn;
1911
3.88M
  isl_int cst;
1912
3.88M
1913
3.88M
  if (!tab)
1914
0
    return isl_stat_error;
1915
3.88M
  if (tab->bmap) {
1916
329k
    struct isl_basic_map *bmap = tab->bmap;
1917
329k
1918
329k
    isl_assert(tab->mat->ctx, tab->n_eq == bmap->n_eq,
1919
329k
      return isl_stat_error);
1920
329k
    isl_assert(tab->mat->ctx,
1921
329k
          tab->n_con == bmap->n_eq + bmap->n_ineq,
1922
329k
          return isl_stat_error);
1923
329k
    tab->bmap = isl_basic_map_add_ineq(tab->bmap, ineq);
1924
329k
    if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1925
0
      return isl_stat_error;
1926
329k
    if (!tab->bmap)
1927
0
      return isl_stat_error;
1928
3.88M
  }
1929
3.88M
  if (tab->cone) {
1930
3.92k
    isl_int_init(cst);
1931
3.92k
    isl_int_set_si(cst, 0);
1932
3.92k
    isl_int_swap(ineq[0], cst);
1933
3.92k
  }
1934
3.88M
  r = isl_tab_add_row(tab, ineq);
1935
3.88M
  if (tab->cone) {
1936
3.92k
    isl_int_swap(ineq[0], cst);
1937
3.92k
    isl_int_clear(cst);
1938
3.92k
  }
1939
3.88M
  if (r < 0)
1940
0
    return isl_stat_error;
1941
3.88M
  tab->con[r].is_nonneg = 1;
1942
3.88M
  if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1943
0
    return isl_stat_error;
1944
3.88M
  if (isl_tab_row_is_redundant(tab, tab->con[r].index)) {
1945
102k
    if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0)
1946
0
      return isl_stat_error;
1947
102k
    return isl_stat_ok;
1948
102k
  }
1949
3.77M
1950
3.77M
  sgn = restore_row(tab, &tab->con[r]);
1951
3.77M
  if (sgn < -1)
1952
0
    return isl_stat_error;
1953
3.77M
  if (sgn < 0)
1954
70.6k
    return isl_tab_mark_empty(tab);
1955
3.70M
  if (tab->con[r].is_row && 
isl_tab_row_is_redundant(tab, tab->con[r].index)3.32M
)
1956
0
    if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0)
1957
0
      return isl_stat_error;
1958
3.70M
  return isl_stat_ok;
1959
3.70M
}
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
102k
{
1967
102k
  int i;
1968
102k
  int row, col;
1969
102k
  unsigned off = 2 + tab->M;
1970
102k
1971
102k
  if (!var->is_row)
1972
215
    return 0;
1973
102k
1974
124k
  
while (102k
isl_int_is_pos(tab->mat->row[var->index][1])) {
1975
121k
    find_pivot(tab, var, NULL, -1, &row, &col);
1976
121k
    isl_assert(tab->mat->ctx, row != -1, return -1);
1977
121k
    if (isl_tab_pivot(tab, row, col) < 0)
1978
0
      return -1;
1979
121k
    if (!var->is_row)
1980
98.8k
      return 0;
1981
121k
  }
1982
102k
1983
102k
  
for (i = tab->n_dead; 3.60k
i < tab->n_col6.46k
;
++i2.85k
)
1984
6.46k
    if (!isl_int_is_zero(tab->mat->row[var->index][off + i]))
1985
6.46k
      
break3.60k
;
1986
3.60k
1987
3.60k
  isl_assert(tab->mat->ctx, i < tab->n_col, return -1);
1988
3.60k
  if (isl_tab_pivot(tab, var->index, i) < 0)
1989
0
    return -1;
1990
3.60k
1991
3.60k
  return 0;
1992
3.60k
}
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
333k
{
2004
333k
  int i;
2005
333k
  int r;
2006
333k
2007
333k
  if (!tab)
2008
0
    return NULL;
2009
333k
  r = isl_tab_add_row(tab, eq);
2010
333k
  if (r < 0)
2011
0
    goto error;
2012
333k
2013
333k
  r = tab->con[r].index;
2014
333k
  i = isl_seq_first_non_zero(tab->mat->row[r] + 2 + tab->M + tab->n_dead,
2015
333k
          tab->n_col - tab->n_dead);
2016
333k
  isl_assert(tab->mat->ctx, i >= 0, goto error);
2017
333k
  i += tab->n_dead;
2018
333k
  if (isl_tab_pivot(tab, r, i) < 0)
2019
0
    goto error;
2020
333k
  if (isl_tab_kill_col(tab, i) < 0)
2021
0
    goto error;
2022
333k
  tab->n_eq++;
2023
333k
2024
333k
  return tab;
2025
0
error:
2026
0
  isl_tab_free(tab);
2027
0
  return NULL;
2028
333k
}
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
86.2k
{
2035
86.2k
  return tab->M && 
!0
isl_int_is_zero0
(tab->mat->row[row][2]);
2036
86.2k
}
2037
2038
static int row_is_manifestly_zero(struct isl_tab *tab, int row)
2039
114k
{
2040
114k
  unsigned off = 2 + tab->M;
2041
114k
2042
114k
  if (!isl_int_is_zero(tab->mat->row[row][1]))
2043
114k
    
return 099.3k
;
2044
14.7k
  if (row_is_big(tab, row))
2045
0
    return 0;
2046
14.7k
  return isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
2047
14.7k
          tab->n_col - tab->n_dead) == -1;
2048
14.7k
}
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
99.2k
{
2057
99.2k
  struct isl_tab_var *var;
2058
99.2k
  int r;
2059
99.2k
2060
99.2k
  if (!tab)
2061
0
    return -1;
2062
99.2k
  r = isl_tab_add_row(tab, eq);
2063
99.2k
  if (r < 0)
2064
0
    return -1;
2065
99.2k
2066
99.2k
  var = &tab->con[r];
2067
99.2k
  r = var->index;
2068
99.2k
  if (row_is_manifestly_zero(tab, r)) {
2069
790
    var->is_zero = 1;
2070
790
    if (isl_tab_mark_redundant(tab, r) < 0)
2071
0
      return -1;
2072
790
    return 0;
2073
790
  }
2074
98.4k
2075
98.4k
  if (isl_int_is_neg(tab->mat->row[r][1])) {
2076
34.9k
    isl_seq_neg(tab->mat->row[r] + 1, tab->mat->row[r] + 1,
2077
34.9k
          1 + tab->n_col);
2078
34.9k
    var->negated = 1;
2079
34.9k
  }
2080
98.4k
  var->is_nonneg = 1;
2081
98.4k
  if (to_col(tab, var) < 0)
2082
0
    return -1;
2083
98.4k
  var->is_nonneg = 0;
2084
98.4k
  if (isl_tab_kill_col(tab, var->index) < 0)
2085
0
    return -1;
2086
98.4k
2087
98.4k
  return 0;
2088
98.4k
}
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
3.10k
{
2098
3.10k
  int r;
2099
3.10k
  isl_int *row;
2100
3.10k
2101
3.10k
  r = isl_tab_allocate_con(tab);
2102
3.10k
  if (r < 0)
2103
0
    return -1;
2104
3.10k
2105
3.10k
  row = tab->mat->row[tab->con[r].index];
2106
3.10k
  isl_seq_clr(row + 1, 1 + tab->M + tab->n_col);
2107
3.10k
  isl_int_set_si(row[0], 1);
2108
3.10k
2109
3.10k
  return r;
2110
3.10k
}
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
14.8k
{
2121
14.8k
  struct isl_tab_undo *snap = NULL;
2122
14.8k
  struct isl_tab_var *var;
2123
14.8k
  int r;
2124
14.8k
  int row;
2125
14.8k
  int sgn;
2126
14.8k
  isl_int cst;
2127
14.8k
2128
14.8k
  if (!tab)
2129
0
    return -1;
2130
14.8k
  isl_assert(tab->mat->ctx, !tab->M, return -1);
2131
14.8k
2132
14.8k
  if (tab->need_undo)
2133
14.2k
    snap = isl_tab_snap(tab);
2134
14.8k
2135
14.8k
  if (tab->cone) {
2136
1.12k
    isl_int_init(cst);
2137
1.12k
    isl_int_set_si(cst, 0);
2138
1.12k
    isl_int_swap(eq[0], cst);
2139
1.12k
  }
2140
14.8k
  r = isl_tab_add_row(tab, eq);
2141
14.8k
  if (tab->cone) {
2142
1.12k
    isl_int_swap(eq[0], cst);
2143
1.12k
    isl_int_clear(cst);
2144
1.12k
  }
2145
14.8k
  if (r < 0)
2146
0
    return -1;
2147
14.8k
2148
14.8k
  var = &tab->con[r];
2149
14.8k
  row = var->index;
2150
14.8k
  if (row_is_manifestly_zero(tab, row)) {
2151
10.6k
    if (snap)
2152
10.5k
      return isl_tab_rollback(tab, snap);
2153
50
    return drop_row(tab, row);
2154
50
  }
2155
4.21k
2156
4.21k
  if (tab->bmap) {
2157
3.10k
    tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);
2158
3.10k
    if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
2159
0
      return -1;
2160
3.10k
    isl_seq_neg(eq, eq, 1 + tab->n_var);
2161
3.10k
    tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);
2162
3.10k
    isl_seq_neg(eq, eq, 1 + tab->n_var);
2163
3.10k
    if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
2164
0
      return -1;
2165
3.10k
    if (!tab->bmap)
2166
0
      return -1;
2167
3.10k
    if (add_zero_row(tab) < 0)
2168
0
      return -1;
2169
4.21k
  }
2170
4.21k
2171
4.21k
  sgn = isl_int_sgn(tab->mat->row[row][1]);
2172
4.21k
2173
4.21k
  if (sgn > 0) {
2174
262
    isl_seq_neg(tab->mat->row[row] + 1, tab->mat->row[row] + 1,
2175
262
          1 + tab->n_col);
2176
262
    var->negated = 1;
2177
262
    sgn = -1;
2178
262
  }
2179
4.21k
2180
4.21k
  if (sgn < 0) {
2181
2.94k
    sgn = sign_of_max(tab, var);
2182
2.94k
    if (sgn < -1)
2183
0
      return -1;
2184
2.94k
    if (sgn < 0) {
2185
0
      if (isl_tab_mark_empty(tab) < 0)
2186
0
        return -1;
2187
0
      return 0;
2188
0
    }
2189
2.94k
  }
2190
4.21k
2191
4.21k
  var->is_nonneg = 1;
2192
4.21k
  if (to_col(tab, var) < 0)
2193
0
    return -1;
2194
4.21k
  var->is_nonneg = 0;
2195
4.21k
  if (isl_tab_kill_col(tab, var->index) < 0)
2196
0
    return -1;
2197
4.21k
2198
4.21k
  return 0;
2199
4.21k
}
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
3.64k
{
2213
3.64k
  unsigned total;
2214
3.64k
  unsigned div_pos;
2215
3.64k
  struct isl_vec *ineq;
2216
3.64k
2217
3.64k
  if (!bmap)
2218
0
    return NULL;
2219
3.64k
2220
3.64k
  total = isl_basic_map_total_dim(bmap);
2221
3.64k
  div_pos = 1 + total - bmap->n_div + div;
2222
3.64k
2223
3.64k
  ineq = isl_vec_alloc(bmap->ctx, 1 + total);
2224
3.64k
  if (!ineq)
2225
0
    return NULL;
2226
3.64k
2227
3.64k
  isl_seq_cpy(ineq->el, bmap->div[div] + 1, 1 + total);
2228
3.64k
  isl_int_neg(ineq->el[div_pos], bmap->div[div][0]);
2229
3.64k
  return ineq;
2230
3.64k
}
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
3.64k
{
2250
3.64k
  unsigned total;
2251
3.64k
  unsigned div_pos;
2252
3.64k
  struct isl_vec *ineq;
2253
3.64k
2254
3.64k
  total = isl_basic_map_total_dim(tab->bmap);
2255
3.64k
  div_pos = 1 + total - tab->bmap->n_div + div;
2256
3.64k
2257
3.64k
  ineq = ineq_for_div(tab->bmap, div);
2258
3.64k
  if (!ineq)
2259
0
    goto error;
2260
3.64k
2261
3.64k
  if (add_ineq) {
2262
624
    if (add_ineq(user, ineq->el) < 0)
2263
0
      goto error;
2264
3.02k
  } else {
2265
3.02k
    if (isl_tab_add_ineq(tab, ineq->el) < 0)
2266
0
      goto error;
2267
3.64k
  }
2268
3.64k
2269
3.64k
  isl_seq_neg(ineq->el, tab->bmap->div[div] + 1, 1 + total);
2270
3.64k
  isl_int_set(ineq->el[div_pos], tab->bmap->div[div][0]);
2271
3.64k
  isl_int_add(ineq->el[0], ineq->el[0], ineq->el[div_pos]);
2272
3.64k
  isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
2273
3.64k
2274
3.64k
  if (add_ineq) {
2275
624
    if (add_ineq(user, ineq->el) < 0)
2276
0
      goto error;
2277
3.02k
  } else {
2278
3.02k
    if (isl_tab_add_ineq(tab, ineq->el) < 0)
2279
0
      goto error;
2280
3.64k
  }
2281
3.64k
2282
3.64k
  isl_vec_free(ineq);
2283
3.64k
2284
3.64k
  return isl_stat_ok;
2285
0
error:
2286
0
  isl_vec_free(ineq);
2287
0
  return isl_stat_error;
2288
3.64k
}
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
3.64k
{
2298
3.64k
  int i;
2299
3.64k
2300
3.64k
  if (tab->M)
2301
0
    return 1;
2302
3.64k
2303
3.64k
  if (isl_int_is_neg(div->el[1]))
2304
3.64k
    
return 0698
;
2305
2.95k
2306
9.12k
  
for (i = 0; 2.95k
i < tab->n_var;
++i6.17k
) {
2307
8.60k
    if (isl_int_is_neg(div->el[2 + i]))
2308
8.60k
      
return 0279
;
2309
8.33k
    if (isl_int_is_zero(div->el[2 + i]))
2310
8.33k
      
continue5.49k
;
2311
2.83k
    if (!tab->var[i].is_nonneg)
2312
2.15k
      return 0;
2313
2.83k
  }
2314
2.95k
2315
2.95k
  
return 1520
;
2316
2.95k
}
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
3.64k
{
2332
3.64k
  int r;
2333
3.64k
  int nonneg;
2334
3.64k
  int n_div, o_div;
2335
3.64k
2336
3.64k
  if (!tab || !div)
2337
0
    return -1;
2338
3.64k
2339
3.64k
  if (div->size != 1 + 1 + tab->n_var)
2340
3.64k
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_invalid,
2341
3.64k
      "unexpected size", return -1);
2342
3.64k
2343
3.64k
  isl_assert(tab->mat->ctx, tab->bmap, return -1);
2344
3.64k
  n_div = isl_basic_map_dim(tab->bmap, isl_dim_div);
2345
3.64k
  o_div = tab->n_var - n_div;
2346
3.64k
  if (pos < o_div || pos > tab->n_var)
2347
3.64k
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_invalid,
2348
3.64k
      "invalid position", return -1);
2349
3.64k
2350
3.64k
  nonneg = div_is_nonneg(tab, div);
2351
3.64k
2352
3.64k
  if (isl_tab_extend_cons(tab, 3) < 0)
2353
0
    return -1;
2354
3.64k
  if (isl_tab_extend_vars(tab, 1) < 0)
2355
0
    return -1;
2356
3.64k
  r = isl_tab_insert_var(tab, pos);
2357
3.64k
  if (r < 0)
2358
0
    return -1;
2359
3.64k
2360
3.64k
  if (nonneg)
2361
520
    tab->var[r].is_nonneg = 1;
2362
3.64k
2363
3.64k
  tab->bmap = isl_basic_map_insert_div(tab->bmap, pos - o_div, div);
2364
3.64k
  if (!tab->bmap)
2365
0
    return -1;
2366
3.64k
  if (isl_tab_push_var(tab, isl_tab_undo_bmap_div, &tab->var[r]) < 0)
2367
0
    return -1;
2368
3.64k
2369
3.64k
  if (add_div_constraints(tab, pos - o_div, add_ineq, user) < 0)
2370
0
    return -1;
2371
3.64k
2372
3.64k
  return r;
2373
3.64k
}
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.02k
{
2380
3.02k
  if (!tab)
2381
0
    return -1;
2382
3.02k
  return isl_tab_insert_div(tab, tab->n_var, div, NULL, NULL);
2383
3.02k
}
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
721k
{
2393
721k
  int i;
2394
721k
  struct isl_tab *tab;
2395
721k
2396
721k
  if (!bmap)
2397
0
    return NULL;
2398
721k
  tab = isl_tab_alloc(bmap->ctx,
2399
721k
          isl_basic_map_total_dim(bmap) + bmap->n_ineq + 1,
2400
721k
          isl_basic_map_total_dim(bmap), 0);
2401
721k
  if (!tab)
2402
0
    return NULL;
2403
721k
  tab->preserve = track;
2404
721k
  tab->rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
2405
721k
  if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY)) {
2406
7
    if (isl_tab_mark_empty(tab) < 0)
2407
0
      goto error;
2408
7
    goto done;
2409
7
  }
2410
1.05M
  
for (i = 0; 721k
i < bmap->n_eq;
++i332k
) {
2411
332k
    tab = add_eq(tab, bmap->eq[i]);
2412
332k
    if (!tab)
2413
0
      return tab;
2414
332k
  }
2415
4.23M
  
for (i = 0; 721k
i < bmap->n_ineq;
++i3.51M
) {
2416
3.52M
    if (isl_tab_add_ineq(tab, bmap->ineq[i]) < 0)
2417
0
      goto error;
2418
3.52M
    if (tab->empty)
2419
5.60k
      goto done;
2420
3.52M
  }
2421
721k
done:
2422
721k
  if (track && 
isl_tab_track_bmap(tab, isl_basic_map_copy(bmap)) < 0177k
)
2423
0
    goto error;
2424
721k
  return tab;
2425
0
error:
2426
0
  isl_tab_free(tab);
2427
0
  return NULL;
2428
721k
}
2429
2430
__isl_give struct isl_tab *isl_tab_from_basic_set(
2431
  __isl_keep isl_basic_set *bset, int track)
2432
256k
{
2433
256k
  return isl_tab_from_basic_map(bset, track);
2434
256k
}
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.62k
{
2441
3.62k
  isl_int cst;
2442
3.62k
  int i;
2443
3.62k
  struct isl_tab *tab;
2444
3.62k
  unsigned offset = 0;
2445
3.62k
2446
3.62k
  if (!bset)
2447
0
    return NULL;
2448
3.62k
  if (parametric)
2449
2.85k
    offset = isl_basic_set_dim(bset, isl_dim_param);
2450
3.62k
  tab = isl_tab_alloc(bset->ctx, bset->n_eq + bset->n_ineq,
2451
3.62k
        isl_basic_set_total_dim(bset) - offset, 0);
2452
3.62k
  if (!tab)
2453
0
    return NULL;
2454
3.62k
  tab->rational = ISL_F_ISSET(bset, ISL_BASIC_SET_RATIONAL);
2455
3.62k
  tab->cone = 1;
2456
3.62k
2457
3.62k
  isl_int_init(cst);
2458
3.62k
  isl_int_set_si(cst, 0);
2459
5.59k
  for (i = 0; i < bset->n_eq; 
++i1.97k
) {
2460
1.97k
    isl_int_swap(bset->eq[i][offset], cst);
2461
1.97k
    if (offset > 0) {
2462
554
      if (isl_tab_add_eq(tab, bset->eq[i] + offset) < 0)
2463
0
        goto error;
2464
1.41k
    } else
2465
1.41k
      tab = add_eq(tab, bset->eq[i]);
2466
1.97k
    isl_int_swap(bset->eq[i][offset], cst);
2467
1.97k
    if (!tab)
2468
0
      goto done;
2469
1.97k
  }
2470
14.9k
  
for (i = 0; 3.62k
i < bset->n_ineq;
++i11.3k
) {
2471
11.3k
    int r;
2472
11.3k
    isl_int_swap(bset->ineq[i][offset], cst);
2473
11.3k
    r = isl_tab_add_row(tab, bset->ineq[i] + offset);
2474
11.3k
    isl_int_swap(bset->ineq[i][offset], cst);
2475
11.3k
    if (r < 0)
2476
0
      goto error;
2477
11.3k
    tab->con[r].is_nonneg = 1;
2478
11.3k
    if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
2479
0
      goto error;
2480
11.3k
  }
2481
3.62k
done:
2482
3.62k
  isl_int_clear(cst);
2483
3.62k
  return tab;
2484
0
error:
2485
0
  isl_int_clear(cst);
2486
0
  isl_tab_free(tab);
2487
0
  return NULL;
2488
3.62k
}
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.85k
{
2495
2.85k
  int i;
2496
2.85k
2497
2.85k
  if (!tab)
2498
0
    return isl_bool_error;
2499
2.85k
  if (tab->empty)
2500
0
    return isl_bool_true;
2501
2.85k
  if (tab->n_dead == tab->n_col)
2502
737
    return isl_bool_true;
2503
2.11k
2504
3.29k
  
for (;;)2.11k
{
2505
3.51k
    for (i = tab->n_redundant; i < tab->n_row; 
++i221
) {
2506
3.51k
      struct isl_tab_var *var;
2507
3.51k
      int sgn;
2508
3.51k
      var = isl_tab_var_from_row(tab, i);
2509
3.51k
      if (!var->is_nonneg)
2510
221
        continue;
2511
3.29k
      sgn = sign_of_max(tab, var);
2512
3.29k
      if (sgn < -1)
2513
0
        return isl_bool_error;
2514
3.29k
      if (sgn != 0)
2515
277
        return isl_bool_false;
2516
3.01k
      if (close_row(tab, var, 0) < 0)
2517
0
        return isl_bool_error;
2518
3.01k
      break;
2519
3.01k
    }
2520
3.29k
    
if (3.01k
tab->n_dead == tab->n_col3.01k
)
2521
1.83k
      return isl_bool_true;
2522
1.18k
    if (i == tab->n_row)
2523
3
      return isl_bool_false;
2524
1.18k
  }
2525
2.11k
}
2526
2527
int isl_tab_sample_is_integer(struct isl_tab *tab)
2528
450k
{
2529
450k
  int i;
2530
450k
2531
450k
  if (!tab)
2532
0
    return -1;
2533
450k
2534
2.20M
  
for (i = 0; 450k
i < tab->n_var;
++i1.75M
) {
2535
1.87M
    int row;
2536
1.87M
    if (!tab->var[i].is_row)
2537
534k
      continue;
2538
1.34M
    row = tab->var[i].index;
2539
1.34M
    if (!isl_int_is_divisible_by(tab->mat->row[row][1],
2540
1.34M
            tab->mat->row[row][0]))
2541
1.34M
      
return 0120k
;
2542
1.34M
  }
2543
450k
  
return 1330k
;
2544
450k
}
2545
2546
static struct isl_vec *extract_integer_sample(struct isl_tab *tab)
2547
214k
{
2548
214k
  int i;
2549
214k
  struct isl_vec *vec;
2550
214k
2551
214k
  vec = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
2552
214k
  if (!vec)
2553
0
    return NULL;
2554
214k
2555
214k
  isl_int_set_si(vec->block.data[0], 1);
2556
1.40M
  for (i = 0; i < tab->n_var; 
++i1.18M
) {
2557
1.18M
    if (!tab->var[i].is_row)
2558
1.18M
      
isl_int_set_si421k
(vec->block.data[1 + i], 0);
2559
1.18M
    else {
2560
767k
      int row = tab->var[i].index;
2561
767k
      isl_int_divexact(vec->block.data[1 + i],
2562
767k
        tab->mat->row[row][1], tab->mat->row[row][0]);
2563
767k
    }
2564
1.18M
  }
2565
214k
2566
214k
  return vec;
2567
214k
}
2568
2569
struct isl_vec *isl_tab_get_sample_value(struct isl_tab *tab)
2570
240k
{
2571
240k
  int i;
2572
240k
  struct isl_vec *vec;
2573
240k
  isl_int m;
2574
240k
2575
240k
  if (!tab)
2576
0
    return NULL;
2577
240k
2578
240k
  vec = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
2579
240k
  if (!vec)
2580
0
    return NULL;
2581
240k
2582
240k
  isl_int_init(m);
2583
240k
2584
240k
  isl_int_set_si(vec->block.data[0], 1);
2585
1.06M
  for (i = 0; i < tab->n_var; 
++i829k
) {
2586
829k
    int row;
2587
829k
    if (!tab->var[i].is_row) {
2588
395k
      isl_int_set_si(vec->block.data[1 + i], 0);
2589
395k
      continue;
2590
395k
    }
2591
433k
    row = tab->var[i].index;
2592
433k
    isl_int_gcd(m, vec->block.data[0], tab->mat->row[row][0]);
2593
433k
    isl_int_divexact(m, tab->mat->row[row][0], m);
2594
433k
    isl_seq_scale(vec->block.data, vec->block.data, m, 1 + i);
2595
433k
    isl_int_divexact(m, vec->block.data[0], tab->mat->row[row][0]);
2596
433k
    isl_int_mul(vec->block.data[1 + i], m, tab->mat->row[row][1]);
2597
433k
  }
2598
240k
  vec = isl_vec_normalize(vec);
2599
240k
2600
240k
  isl_int_clear(m);
2601
240k
  return vec;
2602
240k
}
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
250k
{
2610
250k
  if (!var->is_row)
2611
250k
    
isl_int_set_si3.02k
(*v, 0);
2612
250k
  else 
if (247k
sgn > 0247k
)
2613
247k
    
isl_int_cdiv_q245k
(*v, tab->mat->row[var->index][1],
2614
247k
           tab->mat->row[var->index][0]);
2615
247k
  else
2616
247k
    
isl_int_fdiv_q2.51k
(*v, tab->mat->row[var->index][1],
2617
250k
           tab->mat->row[var->index][0]);
2618
250k
}
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
445k
{
2631
445k
  int i;
2632
445k
  unsigned n_eq;
2633
445k
2634
445k
  if (!bmap)
2635
0
    return NULL;
2636
445k
  if (!tab)
2637
0
    return bmap;
2638
445k
2639
445k
  n_eq = tab->n_eq;
2640
445k
  if (tab->empty)
2641
6.74k
    bmap = isl_basic_map_set_to_empty(bmap);
2642
438k
  else
2643
2.97M
    
for (i = bmap->n_ineq - 1; 438k
i >= 0;
--i2.53M
) {
2644
2.53M
      if (isl_tab_is_equality(tab, n_eq + i))
2645
964k
        isl_basic_map_inequality_to_equality(bmap, i);
2646
1.56M
      else if (isl_tab_is_redundant(tab, n_eq + i))
2647
136k
        isl_basic_map_drop_inequality(bmap, i);
2648
2.53M
    }
2649
445k
  if (bmap->n_eq != n_eq)
2650
149k
    bmap = isl_basic_map_gauss(bmap, NULL);
2651
445k
  if (!tab->rational &&
2652
445k
      
bmap411k
&&
!bmap->sample411k
&&
isl_tab_sample_is_integer(tab)238k
)
2653
214k
    bmap->sample = extract_integer_sample(tab);
2654
445k
  return bmap;
2655
445k
}
2656
2657
struct isl_basic_set *isl_basic_set_update_from_tab(struct isl_basic_set *bset,
2658
  struct isl_tab *tab)
2659
38.6k
{
2660
38.6k
  return bset_from_bmap(isl_basic_map_update_from_tab(bset_to_bmap(bset),
2661
38.6k
                tab));
2662
38.6k
}
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
12.2k
{
2669
12.2k
  if (!tab->con[r].is_row)
2670
12.2k
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_internal,
2671
12.2k
      "row unexpectedly moved to column",
2672
12.2k
      return isl_stat_error);
2673
12.2k
  if (r + 1 != tab->n_con)
2674
12.2k
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_internal,
2675
12.2k
      "additional constraints added", return isl_stat_error);
2676
12.2k
  if (drop_row(tab, tab->con[r].index) < 0)
2677
0
    return isl_stat_error;
2678
12.2k
2679
12.2k
  return isl_stat_ok;
2680
12.2k
}
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
12.2k
{
2694
12.2k
  unsigned r;
2695
12.2k
  isl_int *row;
2696
12.2k
  int sgn;
2697
12.2k
  unsigned off = 2 + tab->M;
2698
12.2k
2699
12.2k
  if (var->is_zero)
2700
0
    return isl_stat_ok;
2701
12.2k
  if (var->is_redundant || !var->is_nonneg)
2702
12.2k
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_invalid,
2703
12.2k
      "expecting non-redundant non-negative variable",
2704
12.2k
      return isl_stat_error);
2705
12.2k
2706
12.2k
  if (isl_tab_extend_cons(tab, 1) < 0)
2707
0
    return isl_stat_error;
2708
12.2k
2709
12.2k
  r = tab->n_con;
2710
12.2k
  tab->con[r].index = tab->n_row;
2711
12.2k
  tab->con[r].is_row = 1;
2712
12.2k
  tab->con[r].is_nonneg = 0;
2713
12.2k
  tab->con[r].is_zero = 0;
2714
12.2k
  tab->con[r].is_redundant = 0;
2715
12.2k
  tab->con[r].frozen = 0;
2716
12.2k
  tab->con[r].negated = 0;
2717
12.2k
  tab->row_var[tab->n_row] = ~r;
2718
12.2k
  row = tab->mat->row[tab->n_row];
2719
12.2k
2720
12.2k
  if (var->is_row) {
2721
1.34k
    isl_int_set(row[0], tab->mat->row[var->index][0]);
2722
1.34k
    isl_seq_neg(row + 1,
2723
1.34k
          tab->mat->row[var->index] + 1, 1 + tab->n_col);
2724
10.8k
  } else {
2725
10.8k
    isl_int_set_si(row[0], 1);
2726
10.8k
    isl_seq_clr(row + 1, 1 + tab->n_col);
2727
10.8k
    isl_int_set_si(row[off + var->index], -1);
2728
10.8k
  }
2729
12.2k
2730
12.2k
  tab->n_row++;
2731
12.2k
  tab->n_con++;
2732
12.2k
2733
12.2k
  sgn = sign_of_max(tab, &tab->con[r]);
2734
12.2k
  if (sgn < -1)
2735
0
    return isl_stat_error;
2736
12.2k
  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
12.1k
  tab->con[r].is_nonneg = 1;
2744
12.1k
  /* sgn == 0 */
2745
12.1k
  if (close_row(tab, &tab->con[r], 1) < 0)
2746
0
    return isl_stat_error;
2747
12.1k
  if (drop_last_con_in_row(tab, r) < 0)
2748
0
    return isl_stat_error;
2749
12.1k
2750
12.1k
  return isl_stat_ok;
2751
12.1k
}
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
1.76k
{
2772
1.76k
  struct isl_tab_var *var;
2773
1.76k
2774
1.76k
  if (!tab)
2775
0
    return -1;
2776
1.76k
2777
1.76k
  var = &tab->con[con];
2778
1.76k
2779
1.76k
  if (var->is_row && 
(57
var->index < 057
||
var->index < tab->n_redundant57
))
2780
1.76k
    
isl_die0
(tab->mat->ctx, isl_error_invalid,
2781
1.76k
      "cannot relax redundant constraint", return -1);
2782
1.76k
  if (!var->is_row && 
(1.70k
var->index < 01.70k
||
var->index < tab->n_dead1.70k
))
2783
1.76k
    
isl_die0
(tab->mat->ctx, isl_error_invalid,
2784
1.76k
      "cannot relax dead constraint", return -1);
2785
1.76k
2786
1.76k
  if (!var->is_row && 
!max_is_manifestly_unbounded(tab, var)1.70k
)
2787
456
    if (to_row(tab, var, 1) < 0)
2788
0
      return -1;
2789
1.76k
  if (!var->is_row && 
!min_is_manifestly_unbounded(tab, var)1.24k
)
2790
18
    if (to_row(tab, var, -1) < 0)
2791
0
      return -1;
2792
1.76k
2793
1.76k
  if (var->is_row) {
2794
531
    isl_int_add(tab->mat->row[var->index][1],
2795
531
        tab->mat->row[var->index][1], tab->mat->row[var->index][0]);
2796
531
    if (restore_row(tab, var) < 0)
2797
0
      return -1;
2798
1.23k
  } else {
2799
1.23k
    int i;
2800
1.23k
    unsigned off = 2 + tab->M;
2801
1.23k
2802
7.16k
    for (i = 0; i < tab->n_row; 
++i5.93k
) {
2803
5.93k
      if (isl_int_is_zero(tab->mat->row[i][off + var->index]))
2804
5.93k
        
continue4.58k
;
2805
1.34k
      isl_int_sub(tab->mat->row[i][1], tab->mat->row[i][1],
2806
1.34k
          tab->mat->row[i][off + var->index]);
2807
1.34k
    }
2808
1.23k
2809
1.23k
  }
2810
1.76k
2811
1.76k
  if (isl_tab_push_var(tab, isl_tab_undo_relax, var) < 0)
2812
0
    return -1;
2813
1.76k
2814
1.76k
  return 0;
2815
1.76k
}
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
158
{
2836
158
  struct isl_tab_var *var;
2837
158
2838
158
  if (!tab)
2839
0
    return -1;
2840
158
  if (isl_int_is_zero(shift))
2841
158
    
return 075
;
2842
83
2843
83
  var = &tab->var[pos];
2844
83
  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
83
    }
2854
15
  }
2855
83
2856
83
  if (var->is_row) {
2857
76
    isl_int_addmul(tab->mat->row[var->index][1],
2858
76
        shift, tab->mat->row[var->index][0]);
2859
76
  } 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
83
2872
83
  return 0;
2873
83
}
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
519
{
2882
519
  struct isl_tab_var *var;
2883
519
2884
519
  if (!tab)
2885
0
    return -1;
2886
519
2887
519
  var = &tab->con[con];
2888
519
  if (!var->is_nonneg)
2889
0
    return 0;
2890
519
2891
519
  var->is_nonneg = 0;
2892
519
  if (isl_tab_push_var(tab, isl_tab_undo_unrestrict, var) < 0)
2893
0
    return -1;
2894
519
2895
519
  return 0;
2896
519
}
2897
2898
int isl_tab_select_facet(struct isl_tab *tab, int con)
2899
11.6k
{
2900
11.6k
  if (!tab)
2901
0
    return -1;
2902
11.6k
2903
11.6k
  return cut_to_hyperplane(tab, &tab->con[con]);
2904
11.6k
}
2905
2906
static int may_be_equality(struct isl_tab *tab, int row)
2907
7.65M
{
2908
7.65M
  return tab->rational ? 
isl_int_is_zero36.3k
(tab->mat->row[row][1])
2909
7.65M
           : 
isl_int_lt7.61M
(tab->mat->row[row][1],
2910
7.65M
              tab->mat->row[row][0]);
2911
7.65M
}
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
2.08M
{
2926
2.08M
  int i;
2927
2.08M
  struct isl_tab_var *var;
2928
2.08M
2929
16.9M
  for (i = tab->n_con - 1; i >= 0; 
--i14.8M
) {
2930
16.7M
    var = &tab->con[i];
2931
16.7M
    if (var->index < 0)
2932
5.88M
      continue;
2933
10.8M
    if (var->is_row && 
var->index < tab->n_redundant8.37M
)
2934
600k
      continue;
2935
10.2M
    if (!var->is_row && 
var->index < tab->n_dead2.49M
)
2936
1.77k
      continue;
2937
10.2M
    if (var->marked)
2938
1.93M
      return var;
2939
10.2M
  }
2940
2.08M
2941
2.08M
  
return NULL155k
;
2942
2.08M
}
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
435k
{
2963
435k
  int i;
2964
435k
  unsigned n_marked;
2965
435k
2966
435k
  if (!tab)
2967
0
    return -1;
2968
435k
  if (tab->empty)
2969
3.95k
    return 0;
2970
431k
  if (tab->n_dead == tab->n_col)
2971
18.1k
    return 0;
2972
413k
2973
413k
  n_marked = 0;
2974
3.00M
  for (i = tab->n_redundant; i < tab->n_row; 
++i2.59M
) {
2975
2.59M
    struct isl_tab_var *var = isl_tab_var_from_row(tab, i);
2976
2.59M
    var->marked = !var->frozen && 
var->is_nonneg2.57M
&&
2977
2.59M
      
may_be_equality(tab, i)2.28M
;
2978
2.59M
    if (var->marked)
2979
1.63M
      n_marked++;
2980
2.59M
  }
2981
2.39M
  for (i = tab->n_dead; i < tab->n_col; 
++i1.97M
) {
2982
1.97M
    struct isl_tab_var *var = var_from_col(tab, i);
2983
1.97M
    var->marked = !var->frozen && 
var->is_nonneg1.97M
;
2984
1.97M
    if (var->marked)
2985
208k
      n_marked++;
2986
1.97M
  }
2987
1.68M
  while (n_marked) {
2988
1.42M
    struct isl_tab_var *var;
2989
1.42M
    int sgn;
2990
1.42M
    var = select_marked(tab);
2991
1.42M
    if (!var)
2992
148k
      break;
2993
1.27M
    var->marked = 0;
2994
1.27M
    n_marked--;
2995
1.27M
    sgn = sign_of_max(tab, var);
2996
1.27M
    if (sgn < 0)
2997
0
      return -1;
2998
1.27M
    if (sgn == 0) {
2999
499k
      if (close_row(tab, var, 0) < 0)
3000
0
        return -1;
3001
774k
    } else if (!tab->rational && 
!at_least_one(tab, var)751k
) {
3002
539
      if (cut_to_hyperplane(tab, var) < 0)
3003
0
        return -1;
3004
539
      return isl_tab_detect_implicit_equalities(tab);
3005
539
    }
3006
11.3M
    
for (i = tab->n_redundant; 1.27M
i < tab->n_row;
++i10.1M
) {
3007
10.1M
      var = isl_tab_var_from_row(tab, i);
3008
10.1M
      if (!var->marked)
3009
4.73M
        continue;
3010
5.36M
      if (may_be_equality(tab, i))
3011
5.31M
        continue;
3012
56.9k
      var->marked = 0;
3013
56.9k
      n_marked--;
3014
56.9k
    }
3015
1.27M
  }
3016
413k
3017
413k
  
return 0412k
;
3018
413k
}
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
57.7k
{
3098
57.7k
  int i;
3099
57.7k
3100
57.7k
  if (!tab || !bmap)
3101
0
    return isl_basic_map_free(bmap);
3102
57.7k
  if (tab->empty)
3103
80
    return bmap;
3104
57.7k
3105
166k
  
for (i = bmap->n_ineq - 1; 57.7k
i >= 0;
--i108k
) {
3106
108k
    if (!isl_tab_is_equality(tab, bmap->n_eq + i))
3107
107k
      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
57.7k
3117
57.7k
  return bmap;
3118
57.7k
}
3119
3120
static int con_is_redundant(struct isl_tab *tab, struct isl_tab_var *var)
3121
930k
{
3122
930k
  if (!tab)
3123
0
    return -1;
3124
930k
  if (tab->rational) {
3125
70.4k
    int sgn = sign_of_min(tab, var);
3126
70.4k
    if (sgn < -1)
3127
0
      return -1;
3128
70.4k
    return sgn >= 0;
3129
860k
  } else {
3130
860k
    int irred = isl_tab_min_at_most_neg_one(tab, var);
3131
860k
    if (irred < 0)
3132
0
      return -1;
3133
860k
    return !irred;
3134
860k
  }
3135
930k
}
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
232k
{
3152
232k
  int i;
3153
232k
  unsigned n_marked;
3154
232k
3155
232k
  if (!tab)
3156
0
    return -1;
3157
232k
  if (tab->empty)
3158
2.95k
    return 0;
3159
229k
  if (tab->n_redundant == tab->n_row)
3160
5.28k
    return 0;
3161
224k
3162
224k
  n_marked = 0;
3163
1.68M
  for (i = tab->n_redundant; i < tab->n_row; 
++i1.46M
) {
3164
1.46M
    struct isl_tab_var *var = isl_tab_var_from_row(tab, i);
3165
1.46M
    var->marked = !var->frozen && 
var->is_nonneg1.32M
;
3166
1.46M
    if (var->marked)
3167
621k
      n_marked++;
3168
1.46M
  }
3169
1.26M
  for (i = tab->n_dead; i < tab->n_col; 
++i1.04M
) {
3170
1.04M
    struct isl_tab_var *var = var_from_col(tab, i);
3171
1.04M
    var->marked = !var->frozen && 
var->is_nonneg985k
&&
3172
1.04M
      
!min_is_manifestly_unbounded(tab, var)465k
;
3173
1.04M
    if (var->marked)
3174
123k
      n_marked++;
3175
1.04M
  }
3176
880k
  while (n_marked) {
3177
663k
    struct isl_tab_var *var;
3178
663k
    int red;
3179
663k
    var = select_marked(tab);
3180
663k
    if (!var)
3181
7.09k
      break;
3182
656k
    var->marked = 0;
3183
656k
    n_marked--;
3184
656k
    red = con_is_redundant(tab, var);
3185
656k
    if (red < 0)
3186
0
      return -1;
3187
656k
    if (red && 
!var->is_redundant111k
)
3188
16.6k
      if (isl_tab_mark_redundant(tab, var->index) < 0)
3189
0
        return -1;
3190
8.89M
    
for (i = tab->n_dead; 656k
i < tab->n_col;
++i8.23M
) {
3191
8.23M
      var = var_from_col(tab, i);
3192
8.23M
      if (!var->marked)
3193
7.89M
        continue;
3194
337k
      if (!min_is_manifestly_unbounded(tab, var))
3195
258k
        continue;
3196
79.2k
      var->marked = 0;
3197
79.2k
      n_marked--;
3198
79.2k
    }
3199
656k
  }
3200
224k
3201
224k
  return 0;
3202
224k
}
3203
3204
int isl_tab_is_equality(struct isl_tab *tab, int con)
3205
2.66M
{
3206
2.66M
  int row;
3207
2.66M
  unsigned off;
3208
2.66M
3209
2.66M
  if (!tab)
3210
0
    return -1;
3211
2.66M
  if (tab->con[con].is_zero)
3212
966k
    return 1;
3213
1.69M
  if (tab->con[con].is_redundant)
3214
138k
    return 0;
3215
1.55M
  if (!tab->con[con].is_row)
3216
862k
    return tab->con[con].index < tab->n_dead;
3217
695k
3218
695k
  row = tab->con[con].index;
3219
695k
3220
695k
  off = 2 + tab->M;
3221
695k
  return isl_int_is_zero(tab->mat->row[row][1]) &&
3222
695k
    
!row_is_big(tab, row)58.7k
&&
3223
695k
    isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
3224
58.7k
          tab->n_col - tab->n_dead) == -1;
3225
695k
}
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
314k
{
3241
314k
  int r;
3242
314k
  enum isl_lp_result res = isl_lp_ok;
3243
314k
  struct isl_tab_var *var;
3244
314k
  struct isl_tab_undo *snap;
3245
314k
3246
314k
  if (!tab)
3247
0
    return isl_lp_error;
3248
314k
3249
314k
  if (tab->empty)
3250
28
    return isl_lp_empty;
3251
314k
3252
314k
  snap = isl_tab_snap(tab);
3253
314k
  r = isl_tab_add_row(tab, f);
3254
314k
  if (r < 0)
3255
0
    return isl_lp_error;
3256
314k
  var = &tab->con[r];
3257
673k
  for (;;) {
3258
673k
    int row, col;
3259
673k
    find_pivot(tab, var, var, -1, &row, &col);
3260
673k
    if (row == var->index) {
3261
9.80k
      res = isl_lp_unbounded;
3262
9.80k
      break;
3263
9.80k
    }
3264
663k
    if (row == -1)
3265
304k
      break;
3266
358k
    if (isl_tab_pivot(tab, row, col) < 0)
3267
0
      return isl_lp_error;
3268
358k
  }
3269
314k
  isl_int_mul(tab->mat->row[var->index][0],
3270
314k
        tab->mat->row[var->index][0], denom);
3271
314k
  if (ISL_FL_ISSET(flags, ISL_TAB_SAVE_DUAL)) {
3272
27.6k
    int i;
3273
27.6k
3274
27.6k
    isl_vec_free(tab->dual);
3275
27.6k
    tab->dual = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_con);
3276
27.6k
    if (!tab->dual)
3277
0
      return isl_lp_error;
3278
27.6k
    isl_int_set(tab->dual->el[0], tab->mat->row[var->index][0]);
3279
1.06M
    for (i = 0; i < tab->n_con; 
++i1.03M
) {
3280
1.03M
      int pos;
3281
1.03M
      if (tab->con[i].is_row) {
3282
629k
        isl_int_set_si(tab->dual->el[1 + i], 0);
3283
629k
        continue;
3284
629k
      }
3285
406k
      pos = 2 + tab->M + tab->con[i].index;
3286
406k
      if (tab->con[i].negated)
3287
406k
        
isl_int_neg72.3k
(tab->dual->el[1 + i],
3288
406k
              tab->mat->row[var->index][pos]);
3289
406k
      else
3290
406k
        
isl_int_set334k
(tab->dual->el[1 + i],
3291
406k
              tab->mat->row[var->index][pos]);
3292
406k
    }
3293
27.6k
  }
3294
314k
  if (opt && 
res == isl_lp_ok314k
) {
3295
304k
    if (opt_denom) {
3296
61.5k
      isl_int_set(*opt, tab->mat->row[var->index][1]);
3297
61.5k
      isl_int_set(*opt_denom, tab->mat->row[var->index][0]);
3298
61.5k
    } else
3299
243k
      get_rounded_sample_value(tab, var, 1, opt);
3300
304k
  }
3301
314k
  if (isl_tab_rollback(tab, snap) < 0)
3302
0
    return isl_lp_error;
3303
314k
  return res;
3304
314k
}
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
1.98M
{
3316
1.98M
  if (!tab)
3317
0
    return -1;
3318
1.98M
  if (con < 0 || con >= tab->n_con)
3319
1.98M
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_invalid,
3320
1.98M
      "position out of bounds", return -1);
3321
1.98M
  if (tab->con[con].is_zero)
3322
157
    return 0;
3323
1.98M
  if (tab->con[con].is_redundant)
3324
246k
    return 1;
3325
1.74M
  return tab->con[con].is_row && 
tab->con[con].index < tab->n_redundant731k
;
3326
1.74M
}
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
5.49k
{
3339
5.49k
  unsigned off = 2 + tab->M;
3340
5.49k
  isl_mat *mat = tab->mat;
3341
5.49k
  int n;
3342
5.49k
  int row;
3343
5.49k
  int pos;
3344
5.49k
3345
5.49k
  if (!var->is_row)
3346
3.02k
    return isl_bool_false;
3347
2.46k
  row = var->index;
3348
2.46k
  if (row_is_big(tab, row))
3349
0
    return isl_bool_false;
3350
2.46k
  n = tab->n_col - tab->n_dead;
3351
2.46k
  pos = isl_seq_first_non_zero(mat->row[row] + off + tab->n_dead, n);
3352
2.46k
  if (pos != -1)
3353
2.28k
    return isl_bool_false;
3354
181
  if (value)
3355
181
    
isl_int_divexact0
(*value, mat->row[row][1], mat->row[row][0]);
3356
181
  return isl_bool_true;
3357
181
}
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
7.80k
{
3372
7.80k
  if (row_is_big(tab, var->index))
3373
0
    return 1;
3374
7.80k
  isl_int_mul(*tmp, tab->mat->row[var->index][0], target);
3375
7.80k
  if (sgn > 0)
3376
3.19k
    return isl_int_ge(tab->mat->row[var->index][1], *tmp);
3377
7.80k
  else
3378
7.80k
    
return 4.61k
isl_int_le4.61k
(tab->mat->row[var->index][1], *tmp);
3379
7.80k
}
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
6.96k
{
3400
6.96k
  int row, col;
3401
6.96k
3402
6.96k
  if (sgn < 0 && 
min_is_manifestly_unbounded(tab, var)5.31k
)
3403
2.10k
    return isl_bool_true;
3404
4.86k
  if (sgn > 0 && 
max_is_manifestly_unbounded(tab, var)1.64k
)
3405
0
    return isl_bool_true;
3406
4.86k
  if (to_row(tab, var, sgn) < 0)
3407
0
    return isl_bool_error;
3408
7.80k
  
while (4.86k
!reached(tab, var, sgn, target, tmp)) {
3409
6.92k
    find_pivot(tab, var, var, sgn, &row, &col);
3410
6.92k
    if (row == -1)
3411
1.87k
      return isl_bool_false;
3412
5.05k
    if (row == var->index)
3413
2.11k
      return isl_bool_true;
3414
2.94k
    if (isl_tab_pivot(tab, row, col) < 0)
3415
0
      return isl_bool_error;
3416
2.94k
  }
3417
4.86k
3418
4.86k
  
return isl_bool_true873
;
3419
4.86k
}
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
5.31k
{
3443
5.31k
  isl_bool reached;
3444
5.31k
  isl_vec *eq;
3445
5.31k
  int pos;
3446
5.31k
  isl_stat r;
3447
5.31k
3448
5.31k
  get_rounded_sample_value(tab, var, -1, target);
3449
5.31k
  isl_int_sub_ui(*target, *target, 1);
3450
5.31k
  reached = var_reaches(tab, var, -1, *target, tmp);
3451
5.31k
  if (reached < 0 || reached)
3452
3.66k
    return isl_bool_not(reached);
3453
1.64k
  get_rounded_sample_value(tab, var, 1, target);
3454
1.64k
  isl_int_add_ui(*target, *target, 1);
3455
1.64k
  reached = var_reaches(tab, var, 1, *target, tmp);
3456
1.64k
  if (reached < 0 || reached)
3457
1.41k
    return isl_bool_not(reached);
3458
230
  get_rounded_sample_value(tab, var, -1, tmp);
3459
230
  isl_int_sub_ui(*target, *target, 1);
3460
230
  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
230
3466
230
  if (isl_tab_extend_cons(tab, 1) < 0)
3467
0
    return isl_bool_error;
3468
230
  eq = isl_vec_alloc(isl_tab_get_ctx(tab), 1 + tab->n_var);
3469
230
  if (!eq)
3470
0
    return isl_bool_error;
3471
230
  pos = var - tab->var;
3472
230
  isl_seq_clr(eq->el + 1, tab->n_var);
3473
230
  isl_int_set_si(eq->el[1 + pos], -1);
3474
230
  isl_int_set(eq->el[0], *target);
3475
230
  r = isl_tab_add_eq(tab, eq->el);
3476
230
  isl_vec_free(eq);
3477
230
3478
230
  return r < 0 ? 
isl_bool_error0
: isl_bool_true;
3479
230
}
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
5.49k
{
3497
5.49k
  isl_int target, tmp;
3498
5.49k
  isl_bool is_cst;
3499
5.49k
3500
5.49k
  if (var->is_row && 
row_is_big(tab, var->index)2.46k
)
3501
0
    return isl_bool_false;
3502
5.49k
  is_cst = is_constant(tab, var, value);
3503
5.49k
  if (is_cst < 0 || is_cst)
3504
181
    return is_cst;
3505
5.31k
3506
5.31k
  if (!value)
3507
5.31k
    
isl_int_init2.55k
(target);
3508
5.31k
  isl_int_init(tmp);
3509
5.31k
3510
5.31k
  is_cst = detect_constant_with_tmp(tab, var,
3511
5.31k
              value ? 
value2.75k
:
&target2.55k
, &tmp);
3512
5.31k
3513
5.31k
  isl_int_clear(tmp);
3514
5.31k
  if (!value)
3515
5.31k
    
isl_int_clear2.55k
(target);
3516
5.31k
3517
5.31k
  return is_cst;
3518
5.31k
}
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.75k
{
3529
2.75k
  if (!tab)
3530
0
    return isl_bool_error;
3531
2.75k
  if (var < 0 || var >= tab->n_var)
3532
2.75k
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_invalid,
3533
2.75k
      "position out of bounds", return isl_bool_error);
3534
2.75k
  if (tab->rational)
3535
0
    return isl_bool_false;
3536
2.75k
3537
2.75k
  return get_constant(tab, &tab->var[var], value);
3538
2.75k
}
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
519
{
3548
519
  int i;
3549
519
3550
519
  if (!tab)
3551
0
    return isl_stat_error;
3552
519
  if (tab->rational)
3553
0
    return isl_stat_ok;
3554
519
3555
3.25k
  
for (i = 0; 519
i < tab->n_var;
++i2.73k
) {
3556
2.73k
    if (get_constant(tab, &tab->var[i], NULL) < 0)
3557
0
      return isl_stat_error;
3558
2.73k
  }
3559
519
3560
519
  return isl_stat_ok;
3561
519
}
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
1.11M
{
3568
1.11M
  if (!tab)
3569
0
    return NULL;
3570
1.11M
  tab->need_undo = 1;
3571
1.11M
  return tab->top;
3572
1.11M
}
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
131
{
3579
131
  if (!tab)
3580
0
    return isl_bool_error;
3581
131
3582
131
  return tab->need_undo;
3583
131
}
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
131
{
3592
131
  if (!tab)
3593
0
    return;
3594
131
3595
131
  free_undo(tab);
3596
131
  tab->need_undo = 0;
3597
131
}
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
1.18k
{
3605
1.18k
  unsigned off = 2 + tab->M;
3606
1.18k
3607
1.18k
  if (!var->is_row && 
!max_is_manifestly_unbounded(tab, var)1.18k
)
3608
116
    if (to_row(tab, var, 1) < 0)
3609
0
      return isl_stat_error;
3610
1.18k
3611
1.18k
  if (var->is_row) {
3612
121
    isl_int_sub(tab->mat->row[var->index][1],
3613
121
        tab->mat->row[var->index][1], tab->mat->row[var->index][0]);
3614
121
    if (var->is_nonneg) {
3615
121
      int sgn = restore_row(tab, var);
3616
121
      isl_assert(tab->mat->ctx, sgn >= 0,
3617
121
        return isl_stat_error);
3618
121
    }
3619
1.06k
  } else {
3620
1.06k
    int i;
3621
1.06k
3622
6.35k
    for (i = 0; i < tab->n_row; 
++i5.29k
) {
3623
5.29k
      if (isl_int_is_zero(tab->mat->row[i][off + var->index]))
3624
5.29k
        
continue4.17k
;
3625
1.11k
      isl_int_add(tab->mat->row[i][1], tab->mat->row[i][1],
3626
1.11k
          tab->mat->row[i][off + var->index]);
3627
1.11k
    }
3628
1.06k
3629
1.06k
  }
3630
1.18k
3631
1.18k
  return isl_stat_ok;
3632
1.18k
}
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
492
{
3641
492
  var->is_nonneg = 1;
3642
492
3643
492
  if (var->is_row && 
restore_row(tab, var) < -1426
)
3644
0
    return isl_stat_error;
3645
492
3646
492
  return isl_stat_ok;
3647
492
}
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
688k
{
3658
688k
  struct isl_tab_var *var;
3659
688k
3660
688k
  if (tab->n_redundant < 1)
3661
688k
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_internal,
3662
688k
      "no redundant rows", return isl_stat_error);
3663
688k
3664
688k
  var = isl_tab_var_from_row(tab, tab->n_redundant - 1);
3665
688k
  var->is_redundant = 0;
3666
688k
  tab->n_redundant--;
3667
688k
  restore_row(tab, var);
3668
688k
3669
688k
  return isl_stat_ok;
3670
688k
}
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
2.43M
{
3676
2.43M
  struct isl_tab_var *var = var_from_index(tab, undo->u.var_index);
3677
2.43M
  switch (undo->type) {
3678
2.43M
  case isl_tab_undo_nonneg:
3679
504k
    var->is_nonneg = 0;
3680
504k
    break;
3681
2.43M
  case isl_tab_undo_redundant:
3682
498k
    if (!var->is_row || var->index != tab->n_redundant - 1)
3683
498k
      
isl_die0
(isl_tab_get_ctx(tab), isl_error_internal,
3684
498k
        "not undoing last redundant row",
3685
498k
        return isl_stat_error);
3686
498k
    return restore_last_redundant(tab);
3687
498k
  case isl_tab_undo_freeze:
3688
254k
    var->frozen = 0;
3689
254k
    break;
3690
498k
  case isl_tab_undo_zero:
3691
82.9k
    var->is_zero = 0;
3692
82.9k
    if (!var->is_row)
3693
82.6k
      tab->n_dead--;
3694
82.9k
    break;
3695
1.08M
  case isl_tab_undo_allocate:
3696
1.08M
    if (undo->u.var_index >= 0) {
3697
4.98k
      isl_assert(tab->mat->ctx, !var->is_row,
3698
4.98k
        return isl_stat_error);
3699
4.98k
      return drop_col(tab, var->index);
3700
1.08M
    }
3701
1.08M
    if (!var->is_row) {
3702
114k
      if (!max_is_manifestly_unbounded(tab, var)) {
3703
83.9k
        if (to_row(tab, var, 1) < 0)
3704
0
          return isl_stat_error;
3705
31.0k
      } else if (!min_is_manifestly_unbounded(tab, var)) {
3706
10.5k
        if (to_row(tab, var, -1) < 0)
3707
0
          return isl_stat_error;
3708
20.4k
      } else
3709
20.4k
        if (to_row(tab, var, 0) < 0)
3710
0
          return isl_stat_error;
3711
1.08M
    }
3712
1.08M
    return drop_row(tab, var->index);
3713
1.08M
  case isl_tab_undo_relax:
3714
1.18k
    return unrelax(tab, var);
3715
1.08M
  case isl_tab_undo_unrestrict:
3716
492
    return ununrestrict(tab, var);
3717
1.08M
  default:
3718
0
    isl_die(tab->mat->ctx, isl_error_internal,
3719
2.43M
      "perform_undo_var called on invalid undo record",
3720
2.43M
      return isl_stat_error);
3721
2.43M
  }
3722
2.43M
3723
2.43M
  
return isl_stat_ok842k
;
3724
2.43M
}
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
141k
{
3735
141k
  if (!tab)
3736
0
    return isl_stat_error;
3737
141k
3738
141k
  if (tab->need_undo)
3739
141k
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_invalid,
3740
141k
      "manually restoring redundant constraints "
3741
141k
      "interferes with undo history",
3742
141k
      return isl_stat_error);
3743
141k
3744
331k
  
while (141k
tab->n_redundant > 0) {
3745
189k
    if (tab->row_var[tab->n_redundant - 1] >= 0) {
3746
174k
      struct isl_tab_var *var;
3747
174k
3748
174k
      var = isl_tab_var_from_row(tab, tab->n_redundant - 1);
3749
174k
      var->is_nonneg = 0;
3750
174k
    }
3751
189k
    restore_last_redundant(tab);
3752
189k
  }
3753
141k
  return isl_stat_ok;
3754
141k
}
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
2.68k
{
3761
2.68k
  int off;
3762
2.68k
3763
2.68k
  off = tab->n_var - isl_basic_map_dim(tab->bmap, isl_dim_div);
3764
2.68k
  if (isl_basic_map_drop_div(tab->bmap, pos - off) < 0)
3765
0
    return isl_stat_error;
3766
2.68k
  if (tab->samples) {
3767
491
    tab->samples = isl_mat_drop_cols(tab->samples, 1 + pos, 1);
3768
491
    if (!tab->samples)
3769
0
      return isl_stat_error;
3770
2.68k
  }
3771
2.68k
3772
2.68k
  return isl_stat_ok;
3773
2.68k
}
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
570
{
3788
570
  int i, j;
3789
570
  int n_extra = 0;
3790
570
  int *extra = NULL;  /* current columns that contain bad stuff */
3791
570
  unsigned off = 2 + tab->M;
3792
570
3793
570
  extra = isl_alloc_array(tab->mat->ctx, int, tab->n_col);
3794
570
  if (tab->n_col && !extra)
3795
0
    goto error;
3796
6.78k
  
for (i = 0; 570
i < tab->n_col;
++i6.21k
) {
3797
50.9k
    for (j = 0; j < tab->n_col; 
++j44.7k
)
3798
49.5k
      if (tab->col_var[i] == col_var[j])
3799
4.75k
        break;
3800
6.21k
    if (j < tab->n_col)
3801
4.75k
      continue;
3802
1.45k
    extra[n_extra++] = i;
3803
1.45k
  }
3804
5.29k
  for (i = 0; i < tab->n_col && 
n_extra > 05.18k
;
++i4.72k
) {
3805
4.72k
    struct isl_tab_var *var;
3806
4.72k
    int row;
3807
4.72k
3808
36.4k
    for (j = 0; j < tab->n_col; 
++j31.7k
)
3809
35.0k
      if (col_var[i] == tab->col_var[j])
3810
3.27k
        break;
3811
4.72k
    if (j < tab->n_col)
3812
3.27k
      continue;
3813
1.45k
    var = var_from_index(tab, col_var[i]);
3814
1.45k
    row = var->index;
3815
1.85k
    for (j = 0; j < n_extra; 
++j396
)
3816
1.85k
      if (!isl_int_is_zero(tab->mat->row[row][off+extra[j]]))
3817
1.85k
        
break1.45k
;
3818
1.45k
    isl_assert(tab->mat->ctx, j < n_extra, goto error);
3819
1.45k
    if (isl_tab_pivot(tab, row, extra[j]) < 0)
3820
0
      goto error;
3821
1.45k
    extra[j] = extra[--n_extra];
3822
1.45k
  }
3823
570
3824
570
  free(extra);
3825
570
  return 0;
3826
0
error:
3827
0
  free(extra);
3828
0
  return -1;
3829
570
}
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.2k
{
3837
23.2k
  int i;
3838
23.2k
3839
28.3k
  for (i = tab->n_sample - 1; i >= 0 && 
tab->n_sample > n26.6k
;
--i5.07k
) {
3840
5.07k
    if (tab->sample_index[i] < n)
3841
1.93k
      continue;
3842
3.13k
3843
3.13k
    if (i != tab->n_sample - 1) {
3844
2.12k
      int t = tab->sample_index[tab->n_sample-1];
3845
2.12k
      tab->sample_index[tab->n_sample-1] = tab->sample_index[i];
3846
2.12k
      tab->sample_index[i] = t;
3847
2.12k
      isl_mat_swap_rows(tab->samples, tab->n_sample-1, i);
3848
2.12k
    }
3849
3.13k
    tab->n_sample--;
3850
3.13k
  }
3851
23.2k
}
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
2.78M
{
3857
2.78M
  switch (undo->type) {
3858
2.78M
  case isl_tab_undo_rational:
3859
9.43k
    tab->rational = 0;
3860
9.43k
    break;
3861
2.78M
  case isl_tab_undo_empty:
3862
62.3k
    tab->empty = 0;
3863
62.3k
    break;
3864
2.78M
  case isl_tab_undo_nonneg:
3865
2.43M
  case isl_tab_undo_redundant:
3866
2.43M
  case isl_tab_undo_freeze:
3867
2.43M
  case isl_tab_undo_zero:
3868
2.43M
  case isl_tab_undo_allocate:
3869
2.43M
  case isl_tab_undo_relax:
3870
2.43M
  case isl_tab_undo_unrestrict:
3871
2.43M
    return perform_undo_var(tab, undo);
3872
2.43M
  case isl_tab_undo_bmap_eq:
3873
0
    return isl_basic_map_free_equality(tab->bmap, 1);
3874
2.43M
  case isl_tab_undo_bmap_ineq:
3875
250k
    return isl_basic_map_free_inequality(tab->bmap, 1);
3876
2.43M
  case isl_tab_undo_bmap_div:
3877
2.68k
    return drop_bmap_div(tab, undo->u.var_index);
3878
2.43M
  case isl_tab_undo_saved_basis:
3879
570
    if (restore_basis(tab, undo->u.col_var) < 0)
3880
0
      return isl_stat_error;
3881
570
    break;
3882
5.47k
  case isl_tab_undo_drop_sample:
3883
5.47k
    tab->n_outside--;
3884
5.47k
    break;
3885
23.2k
  case isl_tab_undo_saved_samples:
3886
23.2k
    drop_samples_since(tab, undo->u.n);
3887
23.2k
    break;
3888
2.29k
  case isl_tab_undo_callback:
3889
2.29k
    return undo->u.callback->run(undo->u.callback);
3890
570
  default:
3891
0
    isl_assert(tab->mat->ctx, 0, return isl_stat_error);
3892
2.78M
  }
3893
2.78M
  
return isl_stat_ok101k
;
3894
2.78M
}
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
1.00M
{
3901
1.00M
  struct isl_tab_undo *undo, *next;
3902
1.00M
3903
1.00M
  if (!tab)
3904
0
    return -1;
3905
1.00M
3906
1.00M
  tab->in_undo = 1;
3907
3.79M
  for (undo = tab->top; undo && undo != &tab->bottom; 
undo = next2.78M
) {
3908
3.08M
    next = undo->next;
3909
3.08M
    if (undo == snap)
3910
294k
      break;
3911
2.78M
    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
2.78M
    free_undo_record(undo);
3918
2.78M
  }
3919
1.00M
  tab->in_undo = 0;
3920
1.00M
  tab->top = undo;
3921
1.00M
  if (!undo)
3922
0
    return -1;
3923
1.00M
  return 0;
3924
1.00M
}
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
79.9k
{
3938
79.9k
  int pos;
3939
79.9k
  unsigned off = 2 + tab->M;
3940
79.9k
3941
79.9k
  if (tab->rational)
3942
9.20k
    return isl_ineq_separate;
3943
70.7k
3944
70.7k
  if (!isl_int_is_one(tab->mat->row[row][0]))
3945
70.7k
    
return isl_ineq_separate265
;
3946
70.5k
3947
70.5k
  pos = isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
3948
70.5k
          tab->n_col - tab->n_dead);
3949
70.5k
  if (pos == -1) {
3950
13.4k
    if (isl_int_is_negone(tab->mat->row[row][1]))
3951
13.4k
      
return isl_ineq_adj_eq11.2k
;
3952
2.21k
    else
3953
2.21k
      return isl_ineq_separate;
3954
57.0k
  }
3955
57.0k
3956
57.0k
  if (!isl_int_eq(tab->mat->row[row][1],
3957
57.0k
      tab->mat->row[row][off + tab->n_dead + pos]))
3958
57.0k
    
return isl_ineq_separate20.9k
;
3959
36.1k
3960
36.1k
  pos = isl_seq_first_non_zero(
3961
36.1k
      tab->mat->row[row] + off + tab->n_dead + pos + 1,
3962
36.1k
      tab->n_col - tab->n_dead - pos - 1);
3963
36.1k
3964
36.1k
  return pos == -1 ? 
isl_ineq_adj_ineq33.3k
:
isl_ineq_separate2.77k
;
3965
36.1k
}
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
430k
{
3977
430k
  enum isl_ineq_type type = isl_ineq_error;
3978
430k
  struct isl_tab_undo *snap = NULL;
3979
430k
  int con;
3980
430k
  int row;
3981
430k
3982
430k
  if (!tab)
3983
0
    return isl_ineq_error;
3984
430k
3985
430k
  if (isl_tab_extend_cons(tab, 1) < 0)
3986
0
    return isl_ineq_error;
3987
430k
3988
430k
  snap = isl_tab_snap(tab);
3989
430k
3990
430k
  con = isl_tab_add_row(tab, ineq);
3991
430k
  if (con < 0)
3992
0
    goto error;
3993
430k
3994
430k
  row = tab->con[con].index;
3995
430k
  if (isl_tab_row_is_redundant(tab, row))
3996
0
    type = isl_ineq_redundant;
3997
430k
  else if (isl_int_is_neg(tab->mat->row[row][1]) &&
3998
430k
     
(157k
tab->rational157k
||
3999
157k
        
isl_int_abs_ge135k
(tab->mat->row[row][1],
4000
157k
           tab->mat->row[row][0]))) {
4001
156k
    int nonneg = at_least_zero(tab, &tab->con[con]);
4002
156k
    if (nonneg < 0)
4003
0
      goto error;
4004
156k
    if (nonneg)
4005
76.5k
      type = isl_ineq_cut;
4006
79.9k
    else
4007
79.9k
      type = separation_type(tab, row);
4008
274k
  } else {
4009
274k
    int red = con_is_redundant(tab, &tab->con[con]);
4010
274k
    if (red < 0)
4011
0
      goto error;
4012
274k
    if (!red)
4013
67.1k
      type = isl_ineq_cut;
4014
207k
    else
4015
207k
      type = isl_ineq_redundant;
4016
274k
  }
4017
430k
4018
430k
  if (isl_tab_rollback(tab, snap))
4019
0
    return isl_ineq_error;
4020
430k
  return type;
4021
0
error:
4022
0
  return isl_ineq_error;
4023
430k
}
4024
4025
isl_stat isl_tab_track_bmap(struct isl_tab *tab, __isl_take isl_basic_map *bmap)
4026
178k
{
4027
178k
  bmap = isl_basic_map_cow(bmap);
4028
178k
  if (!tab || !bmap)
4029
0
    goto error;
4030
178k
4031
178k
  if (tab->empty) {
4032
4.02k
    bmap = isl_basic_map_set_to_empty(bmap);
4033
4.02k
    if (!bmap)
4034
0
      goto error;
4035
4.02k
    tab->bmap = bmap;
4036
4.02k
    return isl_stat_ok;
4037
4.02k
  }
4038
174k
4039
174k
  isl_assert(tab->mat->ctx, tab->n_eq == bmap->n_eq, goto error);
4040
174k
  isl_assert(tab->mat->ctx,
4041
174k
        tab->n_con == bmap->n_eq + bmap->n_ineq, goto error);
4042
174k
4043
174k
  tab->bmap = bmap;
4044
174k
4045
174k
  return isl_stat_ok;
4046
0
error:
4047
0
  isl_basic_map_free(bmap);
4048
0
  return isl_stat_error;
4049
174k
}
4050
4051
isl_stat isl_tab_track_bset(struct isl_tab *tab, __isl_take isl_basic_set *bset)
4052
770
{
4053
770
  return isl_tab_track_bmap(tab, bset_to_bmap(bset));
4054
770
}
4055
4056
__isl_keep isl_basic_set *isl_tab_peek_bset(struct isl_tab *tab)
4057
29.1k
{
4058
29.1k
  if (!tab)
4059
0
    return NULL;
4060
29.1k
4061
29.1k
  return bset_from_bmap(tab->bmap);
4062
29.1k
}
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
}