Coverage Report

Created: 2018-06-24 14:39

/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
741k
{
36
741k
  int i;
37
741k
  struct isl_tab *tab;
38
741k
  unsigned off = 2 + M;
39
741k
40
741k
  tab = isl_calloc_type(ctx, struct isl_tab);
41
741k
  if (!tab)
42
0
    return NULL;
43
741k
  tab->mat = isl_mat_alloc(ctx, n_row, off + n_var);
44
741k
  if (!tab->mat)
45
0
    goto error;
46
741k
  tab->var = isl_alloc_array(ctx, struct isl_tab_var, n_var);
47
741k
  if (n_var && 
!tab->var740k
)
48
0
    goto error;
49
741k
  tab->con = isl_alloc_array(ctx, struct isl_tab_var, n_row);
50
741k
  if (n_row && 
!tab->con741k
)
51
0
    goto error;
52
741k
  tab->col_var = isl_alloc_array(ctx, int, n_var);
53
741k
  if (n_var && 
!tab->col_var740k
)
54
0
    goto error;
55
741k
  tab->row_var = isl_alloc_array(ctx, int, n_row);
56
741k
  if (n_row && 
!tab->row_var741k
)
57
0
    goto error;
58
4.03M
  
for (i = 0; 741k
i < n_var;
++i3.29M
) {
59
3.29M
    tab->var[i].index = i;
60
3.29M
    tab->var[i].is_row = 0;
61
3.29M
    tab->var[i].is_nonneg = 0;
62
3.29M
    tab->var[i].is_zero = 0;
63
3.29M
    tab->var[i].is_redundant = 0;
64
3.29M
    tab->var[i].frozen = 0;
65
3.29M
    tab->var[i].negated = 0;
66
3.29M
    tab->col_var[i] = i;
67
3.29M
  }
68
741k
  tab->n_row = 0;
69
741k
  tab->n_con = 0;
70
741k
  tab->n_eq = 0;
71
741k
  tab->max_con = n_row;
72
741k
  tab->n_col = n_var;
73
741k
  tab->n_var = n_var;
74
741k
  tab->max_var = n_var;
75
741k
  tab->n_param = 0;
76
741k
  tab->n_div = 0;
77
741k
  tab->n_dead = 0;
78
741k
  tab->n_redundant = 0;
79
741k
  tab->strict_redundant = 0;
80
741k
  tab->need_undo = 0;
81
741k
  tab->rational = 0;
82
741k
  tab->empty = 0;
83
741k
  tab->in_undo = 0;
84
741k
  tab->M = M;
85
741k
  tab->cone = 0;
86
741k
  tab->bottom.type = isl_tab_undo_bottom;
87
741k
  tab->bottom.next = NULL;
88
741k
  tab->top = &tab->bottom;
89
741k
90
741k
  tab->n_zero = 0;
91
741k
  tab->n_unbounded = 0;
92
741k
  tab->basis = NULL;
93
741k
94
741k
  return tab;
95
0
error:
96
0
  isl_tab_free(tab);
97
0
  return NULL;
98
741k
}
99
100
isl_ctx *isl_tab_get_ctx(struct isl_tab *tab)
101
3.28M
{
102
3.28M
  return tab ? isl_mat_get_ctx(tab->mat) : NULL;
103
3.28M
}
104
105
int isl_tab_extend_cons(struct isl_tab *tab, unsigned n_new)
106
847k
{
107
847k
  unsigned off;
108
847k
109
847k
  if (!tab)
110
0
    return -1;
111
847k
112
847k
  off = 2 + tab->M;
113
847k
114
847k
  if (tab->max_con < tab->n_con + n_new) {
115
91.2k
    struct isl_tab_var *con;
116
91.2k
117
91.2k
    con = isl_realloc_array(tab->mat->ctx, tab->con,
118
91.2k
            struct isl_tab_var, tab->max_con + n_new);
119
91.2k
    if (!con)
120
0
      return -1;
121
91.2k
    tab->con = con;
122
91.2k
    tab->max_con += n_new;
123
91.2k
  }
124
847k
  if (tab->mat->n_row < tab->n_row + n_new) {
125
93.4k
    int *row_var;
126
93.4k
127
93.4k
    tab->mat = isl_mat_extend(tab->mat,
128
93.4k
          tab->n_row + n_new, off + tab->n_col);
129
93.4k
    if (!tab->mat)
130
0
      return -1;
131
93.4k
    row_var = isl_realloc_array(tab->mat->ctx, tab->row_var,
132
93.4k
              int, tab->mat->n_row);
133
93.4k
    if (!row_var)
134
0
      return -1;
135
93.4k
    tab->row_var = row_var;
136
93.4k
    if (tab->row_sign) {
137
220
      enum isl_tab_row_sign *s;
138
220
      s = isl_realloc_array(tab->mat->ctx, tab->row_sign,
139
220
          enum isl_tab_row_sign, tab->mat->n_row);
140
220
      if (!s)
141
0
        return -1;
142
220
      tab->row_sign = s;
143
220
    }
144
93.4k
  }
145
847k
  return 0;
146
847k
}
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.92k
{
153
7.92k
  struct isl_tab_var *var;
154
7.92k
  unsigned off = 2 + tab->M;
155
7.92k
156
7.92k
  if (tab->max_var < tab->n_var + n_new) {
157
5.28k
    var = isl_realloc_array(tab->mat->ctx, tab->var,
158
5.28k
            struct isl_tab_var, tab->n_var + n_new);
159
5.28k
    if (!var)
160
0
      return -1;
161
5.28k
    tab->var = var;
162
5.28k
    tab->max_var = tab->n_var + n_new;
163
5.28k
  }
164
7.92k
165
7.92k
  if (tab->mat->n_col < off + tab->n_col + n_new) {
166
3.54k
    int *p;
167
3.54k
168
3.54k
    tab->mat = isl_mat_extend(tab->mat,
169
3.54k
            tab->mat->n_row, off + tab->n_col + n_new);
170
3.54k
    if (!tab->mat)
171
0
      return -1;
172
3.54k
    p = isl_realloc_array(tab->mat->ctx, tab->col_var,
173
3.54k
              int, tab->n_col + n_new);
174
3.54k
    if (!p)
175
0
      return -1;
176
3.54k
    tab->col_var = p;
177
3.54k
  }
178
7.92k
179
7.92k
  return 0;
180
7.92k
}
181
182
static void free_undo_record(struct isl_tab_undo *undo)
183
3.40M
{
184
3.40M
  switch (undo->type) {
185
3.40M
  case isl_tab_undo_saved_basis:
186
778
    free(undo->u.col_var);
187
778
    break;
188
3.40M
  
default:;3.40M
189
3.40M
  }
190
3.40M
  free(undo);
191
3.40M
}
192
193
static void free_undo(struct isl_tab *tab)
194
747k
{
195
747k
  struct isl_tab_undo *undo, *next;
196
747k
197
1.27M
  for (undo = tab->top; undo && undo != &tab->bottom; 
undo = next531k
) {
198
531k
    next = undo->next;
199
531k
    free_undo_record(undo);
200
531k
  }
201
747k
  tab->top = undo;
202
747k
}
203
204
void isl_tab_free(struct isl_tab *tab)
205
781k
{
206
781k
  if (!tab)
207
33.6k
    return;
208
747k
  free_undo(tab);
209
747k
  isl_mat_free(tab->mat);
210
747k
  isl_vec_free(tab->dual);
211
747k
  isl_basic_map_free(tab->bmap);
212
747k
  free(tab->var);
213
747k
  free(tab->con);
214
747k
  free(tab->row_var);
215
747k
  free(tab->col_var);
216
747k
  free(tab->row_sign);
217
747k
  isl_mat_free(tab->samples);
218
747k
  free(tab->sample_index);
219
747k
  isl_mat_free(tab->basis);
220
747k
  free(tab);
221
747k
}
222
223
struct isl_tab *isl_tab_dup(struct isl_tab *tab)
224
2.80k
{
225
2.80k
  int i;
226
2.80k
  struct isl_tab *dup;
227
2.80k
  unsigned off;
228
2.80k
229
2.80k
  if (!tab)
230
0
    return NULL;
231
2.80k
232
2.80k
  off = 2 + tab->M;
233
2.80k
  dup = isl_calloc_type(tab->mat->ctx, struct isl_tab);
234
2.80k
  if (!dup)
235
0
    return NULL;
236
2.80k
  dup->mat = isl_mat_dup(tab->mat);
237
2.80k
  if (!dup->mat)
238
0
    goto error;
239
2.80k
  dup->var = isl_alloc_array(tab->mat->ctx, struct isl_tab_var, tab->max_var);
240
2.80k
  if (tab->max_var && !dup->var)
241
0
    goto error;
242
26.5k
  
for (i = 0; 2.80k
i < tab->n_var;
++i23.7k
)
243
23.7k
    dup->var[i] = tab->var[i];
244
2.80k
  dup->con = isl_alloc_array(tab->mat->ctx, struct isl_tab_var, tab->max_con);
245
2.80k
  if (tab->max_con && !dup->con)
246
0
    goto error;
247
29.2k
  
for (i = 0; 2.80k
i < tab->n_con;
++i26.4k
)
248
26.4k
    dup->con[i] = tab->con[i];
249
2.80k
  dup->col_var = isl_alloc_array(tab->mat->ctx, int, tab->mat->n_col - off);
250
2.80k
  if ((tab->mat->n_col - off) && !dup->col_var)
251
0
    goto error;
252
13.4k
  
for (i = 0; 2.80k
i < tab->n_col;
++i10.6k
)
253
10.6k
    dup->col_var[i] = tab->col_var[i];
254
2.80k
  dup->row_var = isl_alloc_array(tab->mat->ctx, int, tab->mat->n_row);
255
2.80k
  if (tab->mat->n_row && !dup->row_var)
256
0
    goto error;
257
29.0k
  
for (i = 0; 2.80k
i < tab->n_row;
++i26.2k
)
258
26.2k
    dup->row_var[i] = tab->row_var[i];
259
2.80k
  if (tab->row_sign) {
260
2.79k
    dup->row_sign = isl_alloc_array(tab->mat->ctx, enum isl_tab_row_sign,
261
2.79k
            tab->mat->n_row);
262
2.79k
    if (tab->mat->n_row && !dup->row_sign)
263
0
      goto error;
264
29.0k
    
for (i = 0; 2.79k
i < tab->n_row;
++i26.2k
)
265
26.2k
      dup->row_sign[i] = tab->row_sign[i];
266
2.79k
  }
267
2.80k
  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.80k
  dup->n_row = tab->n_row;
279
2.80k
  dup->n_con = tab->n_con;
280
2.80k
  dup->n_eq = tab->n_eq;
281
2.80k
  dup->max_con = tab->max_con;
282
2.80k
  dup->n_col = tab->n_col;
283
2.80k
  dup->n_var = tab->n_var;
284
2.80k
  dup->max_var = tab->max_var;
285
2.80k
  dup->n_param = tab->n_param;
286
2.80k
  dup->n_div = tab->n_div;
287
2.80k
  dup->n_dead = tab->n_dead;
288
2.80k
  dup->n_redundant = tab->n_redundant;
289
2.80k
  dup->rational = tab->rational;
290
2.80k
  dup->empty = tab->empty;
291
2.80k
  dup->strict_redundant = 0;
292
2.80k
  dup->need_undo = 0;
293
2.80k
  dup->in_undo = 0;
294
2.80k
  dup->M = tab->M;
295
2.80k
  tab->cone = tab->cone;
296
2.80k
  dup->bottom.type = isl_tab_undo_bottom;
297
2.80k
  dup->bottom.next = NULL;
298
2.80k
  dup->top = &dup->bottom;
299
2.80k
300
2.80k
  dup->n_zero = tab->n_zero;
301
2.80k
  dup->n_unbounded = tab->n_unbounded;
302
2.80k
  dup->basis = isl_mat_dup(tab->basis);
303
2.80k
304
2.80k
  return dup;
305
0
error:
306
0
  isl_tab_free(dup);
307
0
  return NULL;
308
2.80k
}
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.12k
{
328
3.12k
  int i;
329
3.12k
  struct isl_mat *prod;
330
3.12k
  unsigned n;
331
3.12k
332
3.12k
  prod = isl_mat_alloc(mat1->ctx, mat1->n_row + mat2->n_row,
333
3.12k
          off + col1 + col2);
334
3.12k
  if (!prod)
335
0
    return NULL;
336
3.12k
337
3.12k
  n = 0;
338
11.1k
  for (i = 0; i < r1; 
++i7.99k
) {
339
7.99k
    isl_seq_cpy(prod->row[n + i], mat1->row[i], off + d1);
340
7.99k
    isl_seq_clr(prod->row[n + i] + off + d1, d2);
341
7.99k
    isl_seq_cpy(prod->row[n + i] + off + d1 + d2,
342
7.99k
        mat1->row[i] + off + d1, col1 - d1);
343
7.99k
    isl_seq_clr(prod->row[n + i] + off + col1 + d1, col2 - d2);
344
7.99k
  }
345
3.12k
346
3.12k
  n += r1;
347
11.1k
  for (i = 0; i < r2; 
++i7.99k
) {
348
7.99k
    isl_seq_cpy(prod->row[n + i], mat2->row[i], off);
349
7.99k
    isl_seq_clr(prod->row[n + i] + off, d1);
350
7.99k
    isl_seq_cpy(prod->row[n + i] + off + d1,
351
7.99k
          mat2->row[i] + off, d2);
352
7.99k
    isl_seq_clr(prod->row[n + i] + off + d1 + d2, col1 - d1);
353
7.99k
    isl_seq_cpy(prod->row[n + i] + off + col1 + d1,
354
7.99k
          mat2->row[i] + off + d2, col2 - d2);
355
7.99k
  }
356
3.12k
357
3.12k
  n += r2;
358
36.5k
  for (i = 0; i < row1 - r1; 
++i33.4k
) {
359
33.4k
    isl_seq_cpy(prod->row[n + i], mat1->row[r1 + i], off + d1);
360
33.4k
    isl_seq_clr(prod->row[n + i] + off + d1, d2);
361
33.4k
    isl_seq_cpy(prod->row[n + i] + off + d1 + d2,
362
33.4k
        mat1->row[r1 + i] + off + d1, col1 - d1);
363
33.4k
    isl_seq_clr(prod->row[n + i] + off + col1 + d1, col2 - d2);
364
33.4k
  }
365
3.12k
366
3.12k
  n += row1 - r1;
367
36.5k
  for (i = 0; i < row2 - r2; 
++i33.4k
) {
368
33.4k
    isl_seq_cpy(prod->row[n + i], mat2->row[r2 + i], off);
369
33.4k
    isl_seq_clr(prod->row[n + i] + off, d1);
370
33.4k
    isl_seq_cpy(prod->row[n + i] + off + d1,
371
33.4k
          mat2->row[r2 + i] + off, d2);
372
33.4k
    isl_seq_clr(prod->row[n + i] + off + d1 + d2, col1 - d1);
373
33.4k
    isl_seq_cpy(prod->row[n + i] + off + col1 + d1,
374
33.4k
          mat2->row[r2 + i] + off + d2, col2 - d2);
375
33.4k
  }
376
3.12k
377
3.12k
  return prod;
378
3.12k
}
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.7k
{
386
58.7k
  if (var->index == -1)
387
275
    return;
388
58.4k
  if (var->is_row && 
var->index >= r141.4k
)
389
33.4k
    var->index += r2;
390
58.4k
  if (!var->is_row && 
var->index >= d117.0k
)
391
15.8k
    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.7k
{
401
58.7k
  if (var->index == -1)
402
275
    return;
403
58.4k
  if (var->is_row) {
404
41.4k
    if (var->index < r2)
405
7.99k
      var->index += r1;
406
33.4k
    else
407
33.4k
      var->index += row1;
408
41.4k
  } else {
409
17.0k
    if (var->index < d2)
410
1.21k
      var->index += d1;
411
15.8k
    else
412
15.8k
      var->index += col1;
413
17.0k
  }
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.12k
{
436
3.12k
  int i;
437
3.12k
  struct isl_tab *prod;
438
3.12k
  unsigned off;
439
3.12k
  unsigned r1, r2, d1, d2;
440
3.12k
441
3.12k
  if (!tab1 || !tab2)
442
0
    return NULL;
443
3.12k
444
3.12k
  isl_assert(tab1->mat->ctx, tab1->M == tab2->M, return NULL);
445
3.12k
  isl_assert(tab1->mat->ctx, tab1->rational == tab2->rational, return NULL);
446
3.12k
  isl_assert(tab1->mat->ctx, tab1->cone == tab2->cone, return NULL);
447
3.12k
  isl_assert(tab1->mat->ctx, !tab1->row_sign, return NULL);
448
3.12k
  isl_assert(tab1->mat->ctx, !tab2->row_sign, return NULL);
449
3.12k
  isl_assert(tab1->mat->ctx, tab1->n_param == 0, return NULL);
450
3.12k
  isl_assert(tab1->mat->ctx, tab2->n_param == 0, return NULL);
451
3.12k
  isl_assert(tab1->mat->ctx, tab1->n_div == 0, return NULL);
452
3.12k
  isl_assert(tab1->mat->ctx, tab2->n_div == 0, return NULL);
453
3.12k
454
3.12k
  off = 2 + tab1->M;
455
3.12k
  r1 = tab1->n_redundant;
456
3.12k
  r2 = tab2->n_redundant;
457
3.12k
  d1 = tab1->n_dead;
458
3.12k
  d2 = tab2->n_dead;
459
3.12k
  prod = isl_calloc_type(tab1->mat->ctx, struct isl_tab);
460
3.12k
  if (!prod)
461
0
    return NULL;
462
3.12k
  prod->mat = tab_mat_product(tab1->mat, tab2->mat,
463
3.12k
        tab1->n_row, tab2->n_row,
464
3.12k
        tab1->n_col, tab2->n_col, off, r1, r2, d1, d2);
465
3.12k
  if (!prod->mat)
466
0
    goto error;
467
3.12k
  prod->var = isl_alloc_array(tab1->mat->ctx, struct isl_tab_var,
468
3.12k
          tab1->max_var + tab2->max_var);
469
3.12k
  if ((tab1->max_var + tab2->max_var) && !prod->var)
470
0
    goto error;
471
20.4k
  
for (i = 0; 3.12k
i < tab1->n_var;
++i17.3k
) {
472
17.3k
    prod->var[i] = tab1->var[i];
473
17.3k
    update_index1(&prod->var[i], r1, r2, d1, d2);
474
17.3k
  }
475
20.4k
  for (i = 0; i < tab2->n_var; 
++i17.3k
) {
476
17.3k
    prod->var[tab1->n_var + i] = tab2->var[i];
477
17.3k
    update_index2(&prod->var[tab1->n_var + i],
478
17.3k
        tab1->n_row, tab1->n_col,
479
17.3k
        r1, r2, d1, d2);
480
17.3k
  }
481
3.12k
  prod->con = isl_alloc_array(tab1->mat->ctx, struct isl_tab_var,
482
3.12k
          tab1->max_con +  tab2->max_con);
483
3.12k
  if ((tab1->max_con + tab2->max_con) && !prod->con)
484
0
    goto error;
485
44.5k
  
for (i = 0; 3.12k
i < tab1->n_con;
++i41.4k
) {
486
41.4k
    prod->con[i] = tab1->con[i];
487
41.4k
    update_index1(&prod->con[i], r1, r2, d1, d2);
488
41.4k
  }
489
44.5k
  for (i = 0; i < tab2->n_con; 
++i41.4k
) {
490
41.4k
    prod->con[tab1->n_con + i] = tab2->con[i];
491
41.4k
    update_index2(&prod->con[tab1->n_con + i],
492
41.4k
        tab1->n_row, tab1->n_col,
493
41.4k
        r1, r2, d1, d2);
494
41.4k
  }
495
3.12k
  prod->col_var = isl_alloc_array(tab1->mat->ctx, int,
496
3.12k
          tab1->n_col + tab2->n_col);
497
3.12k
  if ((tab1->n_col + tab2->n_col) && !prod->col_var)
498
0
    goto error;
499
20.1k
  
for (i = 0; 3.12k
i < tab1->n_col;
++i17.0k
) {
500
17.0k
    int pos = i < d1 ? 
i1.21k
:
i + d215.8k
;
501
17.0k
    prod->col_var[pos] = tab1->col_var[i];
502
17.0k
  }
503
20.1k
  for (i = 0; i < tab2->n_col; 
++i17.0k
) {
504
17.0k
    int pos = i < d2 ? 
d1 + i1.21k
:
tab1->n_col + i15.8k
;
505
17.0k
    int t = tab2->col_var[i];
506
17.0k
    if (t >= 0)
507
241
      t += tab1->n_var;
508
16.8k
    else
509
16.8k
      t -= tab1->n_con;
510
17.0k
    prod->col_var[pos] = t;
511
17.0k
  }
512
3.12k
  prod->row_var = isl_alloc_array(tab1->mat->ctx, int,
513
3.12k
          tab1->mat->n_row + tab2->mat->n_row);
514
3.12k
  if ((tab1->mat->n_row + tab2->mat->n_row) && !prod->row_var)
515
0
    goto error;
516
44.5k
  
for (i = 0; 3.12k
i < tab1->n_row;
++i41.4k
) {
517
41.4k
    int pos = i < r1 ? 
i7.99k
:
i + r233.4k
;
518
41.4k
    prod->row_var[pos] = tab1->row_var[i];
519
41.4k
  }
520
44.5k
  for (i = 0; i < tab2->n_row; 
++i41.4k
) {
521
41.4k
    int pos = i < r2 ? 
r1 + i7.99k
:
tab1->n_row + i33.4k
;
522
41.4k
    int t = tab2->row_var[i];
523
41.4k
    if (t >= 0)
524
17.0k
      t += tab1->n_var;
525
24.3k
    else
526
24.3k
      t -= tab1->n_con;
527
41.4k
    prod->row_var[pos] = t;
528
41.4k
  }
529
3.12k
  prod->samples = NULL;
530
3.12k
  prod->sample_index = NULL;
531
3.12k
  prod->n_row = tab1->n_row + tab2->n_row;
532
3.12k
  prod->n_con = tab1->n_con + tab2->n_con;
533
3.12k
  prod->n_eq = 0;
534
3.12k
  prod->max_con = tab1->max_con + tab2->max_con;
535
3.12k
  prod->n_col = tab1->n_col + tab2->n_col;
536
3.12k
  prod->n_var = tab1->n_var + tab2->n_var;
537
3.12k
  prod->max_var = tab1->max_var + tab2->max_var;
538
3.12k
  prod->n_param = 0;
539
3.12k
  prod->n_div = 0;
540
3.12k
  prod->n_dead = tab1->n_dead + tab2->n_dead;
541
3.12k
  prod->n_redundant = tab1->n_redundant + tab2->n_redundant;
542
3.12k
  prod->rational = tab1->rational;
543
3.12k
  prod->empty = tab1->empty || tab2->empty;
544
3.12k
  prod->strict_redundant = tab1->strict_redundant || tab2->strict_redundant;
545
3.12k
  prod->need_undo = 0;
546
3.12k
  prod->in_undo = 0;
547
3.12k
  prod->M = tab1->M;
548
3.12k
  prod->cone = tab1->cone;
549
3.12k
  prod->bottom.type = isl_tab_undo_bottom;
550
3.12k
  prod->bottom.next = NULL;
551
3.12k
  prod->top = &prod->bottom;
552
3.12k
553
3.12k
  prod->n_zero = 0;
554
3.12k
  prod->n_unbounded = 0;
555
3.12k
  prod->basis = NULL;
556
3.12k
557
3.12k
  return prod;
558
0
error:
559
0
  isl_tab_free(prod);
560
0
  return NULL;
561
3.12k
}
562
563
static struct isl_tab_var *var_from_index(struct isl_tab *tab, int i)
564
117M
{
565
117M
  if (i >= 0)
566
31.4M
    return &tab->var[i];
567
85.9M
  else
568
85.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
88.7M
{
573
88.7M
  return var_from_index(tab, tab->row_var[i]);
574
88.7M
}
575
576
static struct isl_tab_var *var_from_col(struct isl_tab *tab, int i)
577
26.1M
{
578
26.1M
  return var_from_index(tab, tab->col_var[i]);
579
26.1M
}
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.11M
{
588
2.11M
  int i;
589
2.11M
  unsigned off = 2 + tab->M;
590
2.11M
591
2.11M
  if (var->is_row)
592
1.45M
    return 0;
593
3.79M
  
for (i = tab->n_redundant; 667k
i < tab->n_row;
++i3.12M
) {
594
3.33M
    if (!isl_int_is_neg(tab->mat->row[i][off + var->index]))
595
3.33M
      
continue2.87M
;
596
463k
    if (isl_tab_var_from_row(tab, i)->is_nonneg)
597
215k
      return 0;
598
463k
  }
599
667k
  
return 1452k
;
600
667k
}
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.64M
{
609
1.64M
  int i;
610
1.64M
  unsigned off = 2 + tab->M;
611
1.64M
612
1.64M
  if (var->is_row)
613
788k
    return 0;
614
6.41M
  
for (i = tab->n_redundant; 856k
i < tab->n_row;
++i5.56M
) {
615
6.04M
    if (!isl_int_is_pos(tab->mat->row[i][off + var->index]))
616
6.04M
      
continue5.26M
;
617
787k
    if (isl_tab_var_from_row(tab, i)->is_nonneg)
618
488k
      return 0;
619
787k
  }
620
856k
  
return 1367k
;
621
856k
}
622
623
static int row_cmp(struct isl_tab *tab, int r1, int r2, int c, isl_int *t)
624
1.57M
{
625
1.57M
  unsigned off = 2 + tab->M;
626
1.57M
627
1.57M
  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.57M
  }
635
1.57M
  isl_int_mul(*t, tab->mat->row[r1][1], tab->mat->row[r2][off + c]);
636
1.57M
  isl_int_submul(*t, tab->mat->row[r2][1], tab->mat->row[r1][off + c]);
637
1.57M
  return isl_int_sgn(*t);
638
1.57M
}
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.29M
{
663
3.29M
  int j, r, tsgn;
664
3.29M
  isl_int t;
665
3.29M
  unsigned off = 2 + tab->M;
666
3.29M
667
3.29M
  isl_int_init(t);
668
3.29M
  r = -1;
669
42.5M
  for (j = tab->n_redundant; j < tab->n_row; 
++j39.2M
) {
670
39.2M
    if (var && 
j == var->index33.5M
)
671
2.86M
      continue;
672
36.3M
    if (!isl_tab_var_from_row(tab, j)->is_nonneg)
673
7.65M
      continue;
674
28.7M
    if (sgn * isl_int_sgn(tab->mat->row[j][off + c]) >= 0)
675
24.9M
      continue;
676
3.82M
    if (r < 0) {
677
2.24M
      r = j;
678
2.24M
      continue;
679
2.24M
    }
680
1.57M
    tsgn = sgn * row_cmp(tab, r, j, c, &t);
681
1.57M
    if (tsgn < 0 || 
(1.19M
tsgn == 01.19M
&&
682
1.19M
              
tab->row_var[j] < tab->row_var[r]535k
))
683
793k
      r = j;
684
1.57M
  }
685
3.29M
  isl_int_clear(t);
686
3.29M
  return r;
687
3.29M
}
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.2M
  for (j = tab->n_dead; j < tab->n_col; 
++j29.0M
) {
719
29.0M
    if (isl_int_is_zero(tr[j]))
720
29.0M
      
continue22.7M
;
721
6.33M
    if (isl_int_sgn(tr[j]) != sgn &&
722
6.33M
        
var_from_col(tab, j)->is_nonneg3.54M
)
723
2.38M
      continue;
724
3.95M
    if (c < 0 || 
tab->col_var[j] < tab->col_var[c]974k
)
725
3.20M
      c = j;
726
3.95M
  }
727
4.25M
  if (c < 0)
728
1.27M
    return;
729
2.98M
730
2.98M
  sgn *= isl_int_sgn(tr[c]);
731
2.98M
  r = pivot_row(tab, skip_var, sgn, c);
732
2.98M
  *row = r < 0 ? 
var->index1.04M
:
r1.93M
;
733
2.98M
  *col = c;
734
2.98M
}
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.8M
{
744
18.8M
  int i;
745
18.8M
  unsigned off = 2 + tab->M;
746
18.8M
747
18.8M
  if (tab->row_var[row] < 0 && 
!isl_tab_var_from_row(tab, row)->is_nonneg14.6M
)
748
936k
    return 0;
749
17.8M
750
17.8M
  if (isl_int_is_neg(tab->mat->row[row][1]))
751
17.8M
    
return 01.52M
;
752
16.3M
  if (tab->strict_redundant && 
isl_int_is_zero43
(tab->mat->row[row][1]))
753
16.3M
    
return 042
;
754
16.3M
  if (tab->M && 
isl_int_is_neg42.5k
(tab->mat->row[row][2]))
755
16.3M
    
return 02.27k
;
756
16.3M
757
65.8M
  
for (i = tab->n_dead; 16.3M
i < tab->n_col;
++i49.5M
) {
758
64.1M
    if (isl_int_is_zero(tab->mat->row[row][off + i]))
759
64.1M
      
continue44.3M
;
760
19.7M
    if (tab->col_var[i] >= 0)
761
7.70M
      return 0;
762
12.0M
    if (isl_int_is_neg(tab->mat->row[row][off + i]))
763
12.0M
      
return 06.73M
;
764
5.30M
    if (!var_from_col(tab, i)->is_nonneg)
765
168k
      return 0;
766
5.30M
  }
767
16.3M
  
return 11.73M
;
768
16.3M
}
769
770
static void swap_rows(struct isl_tab *tab, int row1, int row2)
771
1.72M
{
772
1.72M
  int t;
773
1.72M
  enum isl_tab_row_sign s;
774
1.72M
775
1.72M
  t = tab->row_var[row1];
776
1.72M
  tab->row_var[row1] = tab->row_var[row2];
777
1.72M
  tab->row_var[row2] = t;
778
1.72M
  isl_tab_var_from_row(tab, row1)->index = row1;
779
1.72M
  isl_tab_var_from_row(tab, row2)->index = row2;
780
1.72M
  tab->mat = isl_mat_swap_rows(tab->mat, row1, row2);
781
1.72M
782
1.72M
  if (!tab->row_sign)
783
1.71M
    return;
784
9.79k
  s = tab->row_sign[row1];
785
9.79k
  tab->row_sign[row1] = tab->row_sign[row2];
786
9.79k
  tab->row_sign[row2] = s;
787
9.79k
}
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
13.0M
{
802
13.0M
  struct isl_tab_undo *undo;
803
13.0M
804
13.0M
  if (!tab)
805
0
    return isl_stat_error;
806
13.0M
  if (!tab->need_undo)
807
9.65M
    return isl_stat_ok;
808
3.40M
809
3.40M
  undo = isl_alloc_type(tab->mat->ctx, struct isl_tab_undo);
810
3.40M
  if (!undo)
811
0
    goto error;
812
3.40M
  undo->type = type;
813
3.40M
  undo->u = u;
814
3.40M
  undo->next = tab->top;
815
3.40M
  tab->top = undo;
816
3.40M
817
3.40M
  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.40M
}
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.5M
{
827
12.5M
  union isl_tab_undo_val u;
828
12.5M
  if (var->is_row)
829
12.3M
    u.var_index = tab->row_var[var->index];
830
252k
  else
831
252k
    u.var_index = tab->col_var[var->index];
832
12.5M
  return push_union(tab, type, u);
833
12.5M
}
834
835
isl_stat isl_tab_push(struct isl_tab *tab, enum isl_tab_undo_type type)
836
425k
{
837
425k
  union isl_tab_undo_val u = { 0 };
838
425k
  return push_union(tab, type, u);
839
425k
}
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
778
{
846
778
  int i;
847
778
  union isl_tab_undo_val u;
848
778
849
778
  u.col_var = isl_alloc_array(tab->mat->ctx, int, tab->n_col);
850
778
  if (tab->n_col && !u.col_var)
851
0
    return isl_stat_error;
852
8.83k
  
for (i = 0; 778
i < tab->n_col;
++i8.05k
)
853
8.05k
    u.col_var[i] = tab->col_var[i];
854
778
  return push_union(tab, isl_tab_undo_saved_basis, u);
855
778
}
856
857
isl_stat isl_tab_push_callback(struct isl_tab *tab,
858
  struct isl_tab_callback *callback)
859
22.3k
{
860
22.3k
  union isl_tab_undo_val u;
861
22.3k
  u.callback = callback;
862
22.3k
  return push_union(tab, isl_tab_undo_callback, u);
863
22.3k
}
864
865
struct isl_tab *isl_tab_init_samples(struct isl_tab *tab)
866
7.91k
{
867
7.91k
  if (!tab)
868
0
    return NULL;
869
7.91k
870
7.91k
  tab->n_sample = 0;
871
7.91k
  tab->n_outside = 0;
872
7.91k
  tab->samples = isl_mat_alloc(tab->mat->ctx, 1, 1 + tab->n_var);
873
7.91k
  if (!tab->samples)
874
0
    goto error;
875
7.91k
  tab->sample_index = isl_alloc_array(tab->mat->ctx, int, 1);
876
7.91k
  if (!tab->sample_index)
877
0
    goto error;
878
7.91k
  return tab;
879
0
error:
880
0
  isl_tab_free(tab);
881
0
  return NULL;
882
7.91k
}
883
884
int isl_tab_add_sample(struct isl_tab *tab, __isl_take isl_vec *sample)
885
11.9k
{
886
11.9k
  if (!tab || !sample)
887
0
    goto error;
888
11.9k
889
11.9k
  if (tab->n_sample + 1 > tab->samples->n_row) {
890
3.87k
    int *t = isl_realloc_array(tab->mat->ctx,
891
3.87k
          tab->sample_index, int, tab->n_sample + 1);
892
3.87k
    if (!t)
893
0
      goto error;
894
3.87k
    tab->sample_index = t;
895
3.87k
  }
896
11.9k
897
11.9k
  tab->samples = isl_mat_extend(tab->samples,
898
11.9k
        tab->n_sample + 1, tab->samples->n_col);
899
11.9k
  if (!tab->samples)
900
0
    goto error;
901
11.9k
902
11.9k
  isl_seq_cpy(tab->samples->row[tab->n_sample], sample->el, sample->size);
903
11.9k
  isl_vec_free(sample);
904
11.9k
  tab->sample_index[tab->n_sample] = tab->n_sample;
905
11.9k
  tab->n_sample++;
906
11.9k
907
11.9k
  return 0;
908
0
error:
909
0
  isl_vec_free(sample);
910
0
  return -1;
911
11.9k
}
912
913
struct isl_tab *isl_tab_drop_sample(struct isl_tab *tab, int s)
914
6.33k
{
915
6.33k
  if (s != tab->n_outside) {
916
3.90k
    int t = tab->sample_index[tab->n_outside];
917
3.90k
    tab->sample_index[tab->n_outside] = tab->sample_index[s];
918
3.90k
    tab->sample_index[s] = t;
919
3.90k
    isl_mat_swap_rows(tab->samples, tab->n_outside, s);
920
3.90k
  }
921
6.33k
  tab->n_outside++;
922
6.33k
  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.33k
927
6.33k
  return tab;
928
6.33k
}
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
29.4k
{
935
29.4k
  union isl_tab_undo_val u;
936
29.4k
937
29.4k
  if (!tab)
938
0
    return isl_stat_error;
939
29.4k
940
29.4k
  u.n = tab->n_sample;
941
29.4k
  return push_union(tab, isl_tab_undo_saved_samples, u);
942
29.4k
}
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.33M
{
958
2.33M
  struct isl_tab_var *var = isl_tab_var_from_row(tab, row);
959
2.33M
  var->is_redundant = 1;
960
2.33M
  isl_assert(tab->mat->ctx, row >= tab->n_redundant, return -1);
961
2.33M
  if (tab->preserve || 
tab->need_undo1.65M
||
tab->row_var[row] >= 01.52M
) {
962
1.56M
    if (tab->row_var[row] >= 0 && 
!var->is_nonneg1.18M
) {
963
1.17M
      var->is_nonneg = 1;
964
1.17M
      if (isl_tab_push_var(tab, isl_tab_undo_nonneg, var) < 0)
965
0
        return -1;
966
1.56M
    }
967
1.56M
    if (row != tab->n_redundant)
968
1.08M
      swap_rows(tab, row, tab->n_redundant);
969
1.56M
    tab->n_redundant++;
970
1.56M
    return isl_tab_push_var(tab, isl_tab_undo_redundant, var);
971
1.56M
  } else {
972
768k
    if (row != tab->n_row - 1)
973
514k
      swap_rows(tab, row, tab->n_row - 1);
974
768k
    isl_tab_var_from_row(tab, tab->n_row - 1)->index = -1;
975
768k
    tab->n_row--;
976
768k
    return 1;
977
768k
  }
978
2.33M
}
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
10.6k
{
987
10.6k
  if (!tab)
988
0
    return -1;
989
10.6k
  if (!tab->rational && 
tab->need_undo10.6k
)
990
10.6k
    if (isl_tab_push(tab, isl_tab_undo_rational) < 0)
991
0
      return -1;
992
10.6k
  tab->rational = 1;
993
10.6k
  return 0;
994
10.6k
}
995
996
isl_stat isl_tab_mark_empty(struct isl_tab *tab)
997
76.3k
{
998
76.3k
  if (!tab)
999
0
    return isl_stat_error;
1000
76.3k
  if (!tab->empty && 
tab->need_undo75.4k
)
1001
64.6k
    if (isl_tab_push(tab, isl_tab_undo_empty) < 0)
1002
0
      return isl_stat_error;
1003
76.3k
  tab->empty = 1;
1004
76.3k
  return isl_stat_ok;
1005
76.3k
}
1006
1007
int isl_tab_freeze_constraint(struct isl_tab *tab, int con)
1008
378k
{
1009
378k
  struct isl_tab_var *var;
1010
378k
1011
378k
  if (!tab)
1012
0
    return -1;
1013
378k
1014
378k
  var = &tab->con[con];
1015
378k
  if (var->frozen)
1016
0
    return 0;
1017
378k
  if (var->index < 0)
1018
27.4k
    return 0;
1019
350k
  var->frozen = 1;
1020
350k
1021
350k
  if (tab->need_undo)
1022
317k
    return isl_tab_push_var(tab, isl_tab_undo_freeze, var);
1023
33.8k
1024
33.8k
  return 0;
1025
33.8k
}
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.28M
{
1046
3.28M
  int i;
1047
3.28M
  struct isl_mat *mat = tab->mat;
1048
3.28M
  unsigned off = 2 + tab->M;
1049
3.28M
1050
3.28M
  if (!tab->row_sign)
1051
3.25M
    return;
1052
29.1k
1053
29.1k
  if (tab->row_sign[row] == 0)
1054
23.0k
    return;
1055
6.12k
  isl_assert(mat->ctx, row_sgn > 0, return);
1056
6.12k
  isl_assert(mat->ctx, tab->row_sign[row] == isl_tab_row_neg, return);
1057
6.12k
  tab->row_sign[row] = isl_tab_row_pos;
1058
55.3k
  for (i = 0; i < tab->n_row; 
++i49.1k
) {
1059
49.1k
    int s;
1060
49.1k
    if (i == row)
1061
6.12k
      continue;
1062
43.0k
    s = isl_int_sgn(mat->row[i][off + col]);
1063
43.0k
    if (!s)
1064
25.3k
      continue;
1065
17.7k
    if (!tab->row_sign[i])
1066
6.94k
      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.80k
)
1070
4.80k
      continue;
1071
6.01k
    tab->row_sign[i] = isl_tab_row_unknown;
1072
6.01k
  }
1073
6.12k
}
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.28M
{
1127
3.28M
  int i, j;
1128
3.28M
  int sgn;
1129
3.28M
  int t;
1130
3.28M
  isl_ctx *ctx;
1131
3.28M
  struct isl_mat *mat = tab->mat;
1132
3.28M
  struct isl_tab_var *var;
1133
3.28M
  unsigned off = 2 + tab->M;
1134
3.28M
1135
3.28M
  ctx = isl_tab_get_ctx(tab);
1136
3.28M
  if (isl_ctx_next_operation(ctx) < 0)
1137
0
    return -1;
1138
3.28M
1139
3.28M
  isl_int_swap(mat->row[row][0], mat->row[row][off + col]);
1140
3.28M
  sgn = isl_int_sgn(mat->row[row][0]);
1141
3.28M
  if (sgn < 0) {
1142
1.89M
    isl_int_neg(mat->row[row][0], mat->row[row][0]);
1143
1.89M
    isl_int_neg(mat->row[row][off + col], mat->row[row][off + col]);
1144
1.89M
  } else
1145
11.5M
    
for (j = 0; 1.38M
j < off - 1 + tab->n_col;
++j10.1M
) {
1146
10.1M
      if (j == off - 1 + col)
1147
1.38M
        continue;
1148
8.78M
      isl_int_neg(mat->row[row][1 + j], mat->row[row][1 + j]);
1149
8.78M
    }
1150
3.28M
  if (!isl_int_is_one(mat->row[row][0]))
1151
3.28M
    
isl_seq_normalize(mat->ctx, mat->row[row], off + tab->n_col)586k
;
1152
41.0M
  for (i = 0; i < tab->n_row; 
++i37.7M
) {
1153
37.7M
    if (i == row)
1154
3.28M
      continue;
1155
34.4M
    if (isl_int_is_zero(mat->row[i][off + col]))
1156
34.4M
      
continue25.3M
;
1157
9.07M
    isl_int_mul(mat->row[i][0], mat->row[i][0], mat->row[row][0]);
1158
95.8M
    for (j = 0; j < off - 1 + tab->n_col; 
++j86.7M
) {
1159
86.7M
      if (j == off - 1 + col)
1160
9.07M
        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
9.07M
    isl_int_mul(mat->row[i][off + col],
1167
9.07M
          mat->row[i][off + col], mat->row[row][off + col]);
1168
9.07M
    if (!isl_int_is_one(mat->row[i][0]))
1169
9.07M
      
isl_seq_normalize(mat->ctx, mat->row[i], off + tab->n_col)4.16M
;
1170
9.07M
  }
1171
3.28M
  t = tab->row_var[row];
1172
3.28M
  tab->row_var[row] = tab->col_var[col];
1173
3.28M
  tab->col_var[col] = t;
1174
3.28M
  var = isl_tab_var_from_row(tab, row);
1175
3.28M
  var->is_row = 1;
1176
3.28M
  var->index = row;
1177
3.28M
  var = var_from_col(tab, col);
1178
3.28M
  var->is_row = 0;
1179
3.28M
  var->index = col;
1180
3.28M
  update_row_sign(tab, row, col, sgn);
1181
3.28M
  if (tab->in_undo)
1182
119k
    return 0;
1183
32.4M
  
for (i = tab->n_redundant; 3.16M
i < tab->n_row;
++i29.2M
) {
1184
29.2M
    if (isl_int_is_zero(mat->row[i][off + col]))
1185
29.2M
      
continue18.2M
;
1186
10.9M
    if (!isl_tab_var_from_row(tab, i)->frozen &&
1187
10.9M
        
isl_tab_row_is_redundant(tab, i)10.6M
) {
1188
1.62M
      int redo = isl_tab_mark_redundant(tab, i);
1189
1.62M
      if (redo < 0)
1190
0
        return -1;
1191
1.62M
      if (redo)
1192
150k
        --i;
1193
1.62M
    }
1194
10.9M
  }
1195
3.16M
  return 0;
1196
3.16M
}
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.70M
{
1207
1.70M
  int r;
1208
1.70M
  unsigned off = 2 + tab->M;
1209
1.70M
1210
1.70M
  if (var->is_row)
1211
1.45M
    return 0;
1212
248k
1213
248k
  if (sign == 0) {
1214
56.9k
    for (r = tab->n_redundant; r < tab->n_row; 
++r35.8k
)
1215
56.9k
      if (!isl_int_is_zero(tab->mat->row[r][off+var->index]))
1216
56.9k
        
break21.0k
;
1217
21.0k
    isl_assert(tab->mat->ctx, r < tab->n_row, return -1);
1218
227k
  } else {
1219
227k
    r = pivot_row(tab, NULL, sign, var->index);
1220
227k
    isl_assert(tab->mat->ctx, r >= 0, return -1);
1221
227k
  }
1222
248k
1223
248k
  return isl_tab_pivot(tab, r, var->index);
1224
248k
}
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.30M
{
1266
1.30M
  int row, col;
1267
1.30M
1268
1.30M
  if (max_is_manifestly_unbounded(tab, var))
1269
97.3k
    return 1;
1270
1.20M
  if (to_row(tab, var, 1) < 0)
1271
0
    return -2;
1272
2.04M
  
while (1.20M
!isl_int_is_pos(tab->mat->row[var->index][1])) {
1273
1.65M
    find_pivot(tab, var, var, 1, &row, &col);
1274
1.65M
    if (row == -1)
1275
583k
      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
224k
      return 1;
1280
1.06M
  }
1281
1.20M
  
return 1397k
;
1282
1.20M
}
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.26M
{
1300
5.26M
  if (!tab->M)
1301
5.26M
    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.66M
{
1311
4.66M
  if (!tab->M)
1312
4.66M
    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
5.03M
{
1326
5.03M
  int row, col;
1327
5.03M
1328
5.26M
  while (row_is_neg(tab, var->index)) {
1329
671k
    find_pivot(tab, var, var, 1, &row, &col);
1330
671k
    if (row == -1)
1331
72.9k
      break;
1332
598k
    if (isl_tab_pivot(tab, row, col) < 0)
1333
0
      return -2;
1334
598k
    if (!var->is_row) /* manifestly unbounded */
1335
368k
      return 1;
1336
598k
  }
1337
5.03M
  
return row_sgn(tab, var->index)4.66M
;
1338
5.03M
}
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
163k
{
1347
163k
  int row, col;
1348
163k
1349
183k
  while (isl_int_is_neg(tab->mat->row[var->index][1])) {
1350
174k
    find_pivot(tab, var, var, 1, &row, &col);
1351
174k
    if (row == -1)
1352
82.0k
      break;
1353
92.3k
    if (row == var->index) /* manifestly unbounded */
1354
72.8k
      return 1;
1355
19.5k
    if (isl_tab_pivot(tab, row, col) < 0)
1356
0
      return -1;
1357
19.5k
  }
1358
163k
  
return !90.6k
isl_int_is_neg90.6k
(tab->mat->row[var->index][1]);
1359
163k
}
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.6k
{
1380
70.6k
  int row, col;
1381
70.6k
  struct isl_tab_var *pivot_var = NULL;
1382
70.6k
1383
70.6k
  if (min_is_manifestly_unbounded(tab, var))
1384
0
    return -1;
1385
70.6k
  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.5k
  if (var->is_redundant)
1407
0
    return 0;
1408
107k
  
while (69.5k
!isl_int_is_neg(tab->mat->row[var->index][1])) {
1409
93.7k
    find_pivot(tab, var, var, -1, &row, &col);
1410
93.7k
    if (row == var->index)
1411
16.2k
      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.5k
  
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
284k
{
1434
284k
  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
284k
  }
1440
284k
  return isl_int_is_neg(tab->mat->row[row][1]) &&
1441
284k
         
isl_int_abs_ge159k
(tab->mat->row[row][1],
1442
284k
            tab->mat->row[row][0]);
1443
284k
}
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
806k
{
1456
806k
  int row, col;
1457
806k
  struct isl_tab_var *pivot_var;
1458
806k
1459
806k
  if (min_is_manifestly_unbounded(tab, var))
1460
141
    return 1;
1461
806k
  if (!var->is_row) {
1462
89.1k
    col = var->index;
1463
89.1k
    row = pivot_row(tab, NULL, -1, col);
1464
89.1k
    pivot_var = var_from_col(tab, col);
1465
89.1k
    if (isl_tab_pivot(tab, row, col) < 0)
1466
0
      return -1;
1467
89.1k
    if (var->is_redundant)
1468
12.3k
      return 0;
1469
76.8k
    if (row_at_most_neg_one(tab, var->index)) {
1470
57.9k
      if (var->is_nonneg) {
1471
57.9k
        if (!pivot_var->is_redundant &&
1472
57.9k
            pivot_var->index == row) {
1473
54.3k
          if (isl_tab_pivot(tab, row, col) < 0)
1474
0
            return -1;
1475
3.58k
        } else
1476
3.58k
          if (restore_row(tab, var) < -1)
1477
0
            return -1;
1478
57.9k
      }
1479
57.9k
      return 1;
1480
57.9k
    }
1481
76.8k
  }
1482
736k
  if (var->is_redundant)
1483
12.0k
    return 0;
1484
855k
  
do 724k
{
1485
855k
    find_pivot(tab, var, var, -1, &row, &col);
1486
855k
    if (row == var->index) {
1487
354k
      if (var->is_nonneg && 
restore_row(tab, var) < -1312k
)
1488
0
        return -1;
1489
354k
      return 1;
1490
354k
    }
1491
500k
    if (row == -1)
1492
197k
      return 0;
1493
303k
    pivot_var = var_from_col(tab, col);
1494
303k
    if (isl_tab_pivot(tab, row, col) < 0)
1495
0
      return -1;
1496
303k
    if (var->is_redundant)
1497
95.2k
      return 0;
1498
207k
  } while (!row_at_most_neg_one(tab, var->index));
1499
724k
  
if (76.7k
var->is_nonneg76.7k
) {
1500
64.3k
    /* pivot back to non-negative value */
1501
64.3k
    if (!pivot_var->is_redundant && pivot_var->index == row)
1502
59.6k
      if (isl_tab_pivot(tab, row, col) < 0)
1503
0
        return -1;
1504
64.3k
    if (restore_row(tab, var) < -1)
1505
0
      return -1;
1506
76.7k
  }
1507
76.7k
  return 1;
1508
76.7k
}
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
692k
{
1515
692k
  int row, col;
1516
692k
  isl_int *r;
1517
692k
1518
692k
  if (max_is_manifestly_unbounded(tab, var))
1519
320k
    return 1;
1520
372k
  if (to_row(tab, var, 1) < 0)
1521
0
    return -1;
1522
372k
  r = tab->mat->row[var->index];
1523
386k
  while (isl_int_lt(r[1], r[0])) {
1524
14.7k
    find_pivot(tab, var, var, 1, &row, &col);
1525
14.7k
    if (row == -1)
1526
560
      return isl_int_ge(r[1], r[0]);
1527
14.2k
    if (row == var->index) /* manifestly unbounded */
1528
82
      return 1;
1529
14.1k
    if (isl_tab_pivot(tab, row, col) < 0)
1530
0
      return -1;
1531
14.1k
  }
1532
372k
  
return 1371k
;
1533
372k
}
1534
1535
static void swap_cols(struct isl_tab *tab, int col1, int col2)
1536
632k
{
1537
632k
  int t;
1538
632k
  unsigned off = 2 + tab->M;
1539
632k
  t = tab->col_var[col1];
1540
632k
  tab->col_var[col1] = tab->col_var[col2];
1541
632k
  tab->col_var[col2] = t;
1542
632k
  var_from_col(tab, col1)->index = col1;
1543
632k
  var_from_col(tab, col2)->index = col2;
1544
632k
  tab->mat = isl_mat_swap_cols(tab->mat, off + col1, off + col2);
1545
632k
}
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
898k
{
1561
898k
  var_from_col(tab, col)->is_zero = 1;
1562
898k
  if (tab->need_undo) {
1563
104k
    if (isl_tab_push_var(tab, isl_tab_undo_zero,
1564
104k
              var_from_col(tab, col)) < 0)
1565
0
      return -1;
1566
104k
    if (col != tab->n_dead)
1567
43.9k
      swap_cols(tab, col, tab->n_dead);
1568
104k
    tab->n_dead++;
1569
104k
    return 0;
1570
793k
  } else {
1571
793k
    if (col != tab->n_col - 1)
1572
587k
      swap_cols(tab, col, tab->n_col - 1);
1573
793k
    var_from_col(tab, tab->n_col - 1)->index = -1;
1574
793k
    tab->n_col--;
1575
793k
    return 1;
1576
793k
  }
1577
898k
}
1578
1579
static int row_is_manifestly_non_integral(struct isl_tab *tab, int row)
1580
3.27M
{
1581
3.27M
  unsigned off = 2 + tab->M;
1582
3.27M
1583
3.27M
  if (tab->M && 
!0
isl_int_eq0
(tab->mat->row[row][2],
1584
3.27M
          tab->mat->row[row][0]))
1585
3.27M
    
return 00
;
1586
3.27M
  if (isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
1587
3.27M
            tab->n_col - tab->n_dead) != -1)
1588
273k
    return 0;
1589
3.00M
1590
3.00M
  return !isl_int_is_divisible_by(tab->mat->row[row][1],
1591
3.00M
          tab->mat->row[row][0]);
1592
3.00M
}
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
583k
{
1599
583k
  int i;
1600
583k
1601
583k
  if (tab->empty)
1602
1
    return 1;
1603
583k
  if (tab->rational)
1604
12.6k
    return 0;
1605
570k
1606
7.07M
  
for (i = 0; 570k
i < tab->n_var;
++i6.50M
) {
1607
6.50M
    if (!tab->var[i].is_row)
1608
3.22M
      continue;
1609
3.27M
    if (row_is_manifestly_non_integral(tab, tab->var[i].index))
1610
90
      return 1;
1611
3.27M
  }
1612
570k
1613
570k
  
return 0570k
;
1614
570k
}
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
583k
{
1638
583k
  int j;
1639
583k
  struct isl_mat *mat = tab->mat;
1640
583k
  unsigned off = 2 + tab->M;
1641
583k
1642
583k
  if (!var->is_nonneg)
1643
583k
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_internal,
1644
583k
      "expecting non-negative variable",
1645
583k
      return isl_stat_error);
1646
583k
  var->is_zero = 1;
1647
583k
  if (!temp_var && 
tab->need_undo569k
)
1648
457
    if (isl_tab_push_var(tab, isl_tab_undo_zero, var) < 0)
1649
0
      return isl_stat_error;
1650
4.53M
  
for (j = tab->n_dead; 583k
j < tab->n_col;
++j3.94M
) {
1651
3.94M
    int recheck;
1652
3.94M
    if (isl_int_is_zero(mat->row[var->index][off + j]))
1653
3.94M
      
continue3.43M
;
1654
512k
    if (isl_int_is_pos(mat->row[var->index][off + j]))
1655
512k
      
isl_die0
(isl_tab_get_ctx(tab), isl_error_internal,
1656
512k
        "row cannot have positive coefficients",
1657
512k
        return isl_stat_error);
1658
512k
    recheck = isl_tab_kill_col(tab, j);
1659
512k
    if (recheck < 0)
1660
0
      return isl_stat_error;
1661
512k
    if (recheck)
1662
499k
      --j;
1663
512k
  }
1664
583k
  if (!temp_var && 
isl_tab_mark_redundant(tab, var->index) < 0569k
)
1665
0
    return isl_stat_error;
1666
583k
  if (tab_is_manifestly_empty(tab) && 
isl_tab_mark_empty(tab) < 091
)
1667
0
    return isl_stat_error;
1668
583k
  return isl_stat_ok;
1669
583k
}
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.27M
{
1679
5.27M
  int r;
1680
5.27M
1681
5.27M
  isl_assert(tab->mat->ctx, tab->n_row < tab->mat->n_row, return -1);
1682
5.27M
  isl_assert(tab->mat->ctx, tab->n_con < tab->max_con, return -1);
1683
5.27M
1684
5.27M
  r = tab->n_con;
1685
5.27M
  tab->con[r].index = tab->n_row;
1686
5.27M
  tab->con[r].is_row = 1;
1687
5.27M
  tab->con[r].is_nonneg = 0;
1688
5.27M
  tab->con[r].is_zero = 0;
1689
5.27M
  tab->con[r].is_redundant = 0;
1690
5.27M
  tab->con[r].frozen = 0;
1691
5.27M
  tab->con[r].negated = 0;
1692
5.27M
  tab->row_var[tab->n_row] = ~r;
1693
5.27M
1694
5.27M
  tab->n_row++;
1695
5.27M
  tab->n_con++;
1696
5.27M
  if (isl_tab_push_var(tab, isl_tab_undo_allocate, &tab->con[r]) < 0)
1697
0
    return -1;
1698
5.27M
1699
5.27M
  return r;
1700
5.27M
}
1701
1702
/* Move the entries in tab->var up one position, starting at "first",
1703
 * creating room for an extra entry at position "first".
1704
 * Since some of the entries of tab->row_var and tab->col_var contain
1705
 * indices into this array, they have to be updated accordingly.
1706
 */
1707
static int var_insert_entry(struct isl_tab *tab, int first)
1708
8.08k
{
1709
8.08k
  int i;
1710
8.08k
1711
8.08k
  if (tab->n_var >= tab->max_var)
1712
8.08k
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_internal,
1713
8.08k
      "not enough room for new variable", return -1);
1714
8.08k
  if (first > tab->n_var)
1715
8.08k
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_internal,
1716
8.08k
      "invalid initial position", return -1);
1717
8.08k
1718
8.95k
  
for (i = tab->n_var - 1; 8.08k
i >= first;
--i869
) {
1719
869
    tab->var[i + 1] = tab->var[i];
1720
869
    if (tab->var[i + 1].is_row)
1721
589
      tab->row_var[tab->var[i + 1].index]++;
1722
280
    else
1723
280
      tab->col_var[tab->var[i + 1].index]++;
1724
869
  }
1725
8.08k
1726
8.08k
  tab->n_var++;
1727
8.08k
1728
8.08k
  return 0;
1729
8.08k
}
1730
1731
/* Drop the entry at position "first" in tab->var, moving all
1732
 * subsequent entries down.
1733
 * Since some of the entries of tab->row_var and tab->col_var contain
1734
 * indices into this array, they have to be updated accordingly.
1735
 */
1736
static int var_drop_entry(struct isl_tab *tab, int first)
1737
5.42k
{
1738
5.42k
  int i;
1739
5.42k
1740
5.42k
  if (first >= tab->n_var)
1741
5.42k
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_internal,
1742
5.42k
      "invalid initial position", return -1);
1743
5.42k
1744
5.42k
  tab->n_var--;
1745
5.42k
1746
6.01k
  for (i = first; i < tab->n_var; 
++i585
) {
1747
585
    tab->var[i] = tab->var[i + 1];
1748
585
    if (tab->var[i + 1].is_row)
1749
574
      tab->row_var[tab->var[i].index]--;
1750
11
    else
1751
11
      tab->col_var[tab->var[i].index]--;
1752
585
  }
1753
5.42k
1754
5.42k
  return 0;
1755
5.42k
}
1756
1757
/* Add a variable to the tableau at position "r" and allocate a column for it.
1758
 * Return the index into the variable array "var", i.e., "r",
1759
 * or -1 on error.
1760
 */
1761
int isl_tab_insert_var(struct isl_tab *tab, int r)
1762
8.08k
{
1763
8.08k
  int i;
1764
8.08k
  unsigned off = 2 + tab->M;
1765
8.08k
1766
8.08k
  isl_assert(tab->mat->ctx, tab->n_col < tab->mat->n_col, return -1);
1767
8.08k
1768
8.08k
  if (var_insert_entry(tab, r) < 0)
1769
0
    return -1;
1770
8.08k
1771
8.08k
  tab->var[r].index = tab->n_col;
1772
8.08k
  tab->var[r].is_row = 0;
1773
8.08k
  tab->var[r].is_nonneg = 0;
1774
8.08k
  tab->var[r].is_zero = 0;
1775
8.08k
  tab->var[r].is_redundant = 0;
1776
8.08k
  tab->var[r].frozen = 0;
1777
8.08k
  tab->var[r].negated = 0;
1778
8.08k
  tab->col_var[tab->n_col] = r;
1779
8.08k
1780
48.3k
  for (i = 0; i < tab->n_row; 
++i40.2k
)
1781
40.2k
    isl_int_set_si(tab->mat->row[i][off + tab->n_col], 0);
1782
8.08k
1783
8.08k
  tab->n_col++;
1784
8.08k
  if (isl_tab_push_var(tab, isl_tab_undo_allocate, &tab->var[r]) < 0)
1785
0
    return -1;
1786
8.08k
1787
8.08k
  return r;
1788
8.08k
}
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.27M
{
1825
5.27M
  int i;
1826
5.27M
  int r;
1827
5.27M
  isl_int *row;
1828
5.27M
  isl_int a, b;
1829
5.27M
  unsigned off = 2 + tab->M;
1830
5.27M
1831
5.27M
  r = isl_tab_allocate_con(tab);
1832
5.27M
  if (r < 0)
1833
0
    return -1;
1834
5.27M
1835
5.27M
  isl_int_init(a);
1836
5.27M
  isl_int_init(b);
1837
5.27M
  row = tab->mat->row[tab->con[r].index];
1838
5.27M
  isl_int_set_si(row[0], 1);
1839
5.27M
  isl_int_set(row[1], line[0]);
1840
5.27M
  isl_seq_clr(row + 2, tab->M + tab->n_col);
1841
48.2M
  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.09M
      isl_int_lcm(a,
1846
8.09M
        row[0], tab->mat->row[tab->var[i].index][0]);
1847
8.09M
      isl_int_swap(a, row[0]);
1848
8.09M
      isl_int_divexact(a, row[0], a);
1849
8.09M
      isl_int_divexact(b,
1850
8.09M
        row[0], tab->mat->row[tab->var[i].index][0]);
1851
8.09M
      isl_int_mul(b, b, line[1 + i]);
1852
8.09M
      isl_seq_combine(row + 1, a, row + 1,
1853
8.09M
          b, tab->mat->row[tab->var[i].index] + 1,
1854
8.09M
          1 + tab->M + tab->n_col);
1855
8.09M
    } else
1856
42.9M
      
isl_int_addmul34.8M
(row[off + tab->var[i].index],
1857
42.9M
              line[1 + i], row[0]);
1858
42.9M
    if (tab->M && 
i >= tab->n_param355k
&&
i < tab->n_var - tab->n_div154k
)
1859
42.9M
      
isl_int_submul151k
(row[2], line[1 + i], row[0]);
1860
42.9M
  }
1861
5.27M
  isl_seq_normalize(tab->mat->ctx, row, off + tab->n_col);
1862
5.27M
  isl_int_clear(a);
1863
5.27M
  isl_int_clear(b);
1864
5.27M
1865
5.27M
  if (tab->row_sign)
1866
42.5k
    tab->row_sign[tab->con[r].index] = isl_tab_row_unknown;
1867
5.27M
1868
5.27M
  return r;
1869
5.27M
}
1870
1871
static isl_stat drop_row(struct isl_tab *tab, int row)
1872
1.11M
{
1873
1.11M
  isl_assert(tab->mat->ctx, ~tab->row_var[row] == tab->n_con - 1,
1874
1.11M
    return isl_stat_error);
1875
1.11M
  if (row != tab->n_row - 1)
1876
121k
    swap_rows(tab, row, tab->n_row - 1);
1877
1.11M
  tab->n_row--;
1878
1.11M
  tab->n_con--;
1879
1.11M
  return isl_stat_ok;
1880
1.11M
}
1881
1882
/* Drop the variable in column "col" along with the column.
1883
 * The column is removed first because it may need to be moved
1884
 * into the last position and this process requires
1885
 * the contents of the col_var array in a state
1886
 * before the removal of the variable.
1887
 */
1888
static isl_stat drop_col(struct isl_tab *tab, int col)
1889
5.42k
{
1890
5.42k
  int var;
1891
5.42k
1892
5.42k
  var = tab->col_var[col];
1893
5.42k
  if (col != tab->n_col - 1)
1894
1.08k
    swap_cols(tab, col, tab->n_col - 1);
1895
5.42k
  tab->n_col--;
1896
5.42k
  if (var_drop_entry(tab, var) < 0)
1897
0
    return isl_stat_error;
1898
5.42k
  return isl_stat_ok;
1899
5.42k
}
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
4.08M
{
1909
4.08M
  int r;
1910
4.08M
  int sgn;
1911
4.08M
  isl_int cst;
1912
4.08M
1913
4.08M
  if (!tab)
1914
0
    return isl_stat_error;
1915
4.08M
  if (tab->bmap) {
1916
337k
    struct isl_basic_map *bmap = tab->bmap;
1917
337k
1918
337k
    isl_assert(tab->mat->ctx, tab->n_eq == bmap->n_eq,
1919
337k
      return isl_stat_error);
1920
337k
    isl_assert(tab->mat->ctx,
1921
337k
          tab->n_con == bmap->n_eq + bmap->n_ineq,
1922
337k
          return isl_stat_error);
1923
337k
    tab->bmap = isl_basic_map_add_ineq(tab->bmap, ineq);
1924
337k
    if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1925
0
      return isl_stat_error;
1926
337k
    if (!tab->bmap)
1927
0
      return isl_stat_error;
1928
4.08M
  }
1929
4.08M
  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
4.08M
  r = isl_tab_add_row(tab, ineq);
1935
4.08M
  if (tab->cone) {
1936
3.92k
    isl_int_swap(ineq[0], cst);
1937
3.92k
    isl_int_clear(cst);
1938
3.92k
  }
1939
4.08M
  if (r < 0)
1940
0
    return isl_stat_error;
1941
4.08M
  tab->con[r].is_nonneg = 1;
1942
4.08M
  if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1943
0
    return isl_stat_error;
1944
4.08M
  if (isl_tab_row_is_redundant(tab, tab->con[r].index)) {
1945
106k
    if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0)
1946
0
      return isl_stat_error;
1947
106k
    return isl_stat_ok;
1948
106k
  }
1949
3.97M
1950
3.97M
  sgn = restore_row(tab, &tab->con[r]);
1951
3.97M
  if (sgn < -1)
1952
0
    return isl_stat_error;
1953
3.97M
  if (sgn < 0)
1954
72.9k
    return isl_tab_mark_empty(tab);
1955
3.90M
  if (tab->con[r].is_row && 
isl_tab_row_is_redundant(tab, tab->con[r].index)3.53M
)
1956
0
    if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0)
1957
0
      return isl_stat_error;
1958
3.90M
  return isl_stat_ok;
1959
3.90M
}
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
91.9k
{
1967
91.9k
  int i;
1968
91.9k
  int row, col;
1969
91.9k
  unsigned off = 2 + tab->M;
1970
91.9k
1971
91.9k
  if (!var->is_row)
1972
245
    return 0;
1973
91.6k
1974
112k
  
while (91.6k
isl_int_is_pos(tab->mat->row[var->index][1])) {
1975
108k
    find_pivot(tab, var, NULL, -1, &row, &col);
1976
108k
    isl_assert(tab->mat->ctx, row != -1, return -1);
1977
108k
    if (isl_tab_pivot(tab, row, col) < 0)
1978
0
      return -1;
1979
108k
    if (!var->is_row)
1980
87.6k
      return 0;
1981
108k
  }
1982
91.6k
1983
91.6k
  
for (i = tab->n_dead; 4.01k
i < tab->n_col7.02k
;
++i3.01k
)
1984
7.02k
    if (!isl_int_is_zero(tab->mat->row[var->index][off + i]))
1985
7.02k
      
break4.01k
;
1986
4.01k
1987
4.01k
  isl_assert(tab->mat->ctx, i < tab->n_col, return -1);
1988
4.01k
  if (isl_tab_pivot(tab, var->index, i) < 0)
1989
0
    return -1;
1990
4.01k
1991
4.01k
  return 0;
1992
4.01k
}
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
269k
{
2004
269k
  int i;
2005
269k
  int r;
2006
269k
2007
269k
  if (!tab)
2008
0
    return NULL;
2009
269k
  r = isl_tab_add_row(tab, eq);
2010
269k
  if (r < 0)
2011
0
    goto error;
2012
269k
2013
269k
  r = tab->con[r].index;
2014
269k
  i = isl_seq_first_non_zero(tab->mat->row[r] + 2 + tab->M + tab->n_dead,
2015
269k
          tab->n_col - tab->n_dead);
2016
269k
  isl_assert(tab->mat->ctx, i >= 0, goto error);
2017
269k
  i += tab->n_dead;
2018
269k
  if (isl_tab_pivot(tab, r, i) < 0)
2019
0
    goto error;
2020
269k
  if (isl_tab_kill_col(tab, i) < 0)
2021
0
    goto error;
2022
269k
  tab->n_eq++;
2023
269k
2024
269k
  return tab;
2025
0
error:
2026
0
  isl_tab_free(tab);
2027
0
  return NULL;
2028
269k
}
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
85.5k
{
2035
85.5k
  return tab->M && 
!0
isl_int_is_zero0
(tab->mat->row[row][2]);
2036
85.5k
}
2037
2038
static int row_is_manifestly_zero(struct isl_tab *tab, int row)
2039
103k
{
2040
103k
  unsigned off = 2 + tab->M;
2041
103k
2042
103k
  if (!isl_int_is_zero(tab->mat->row[row][1]))
2043
103k
    
return 088.2k
;
2044
15.4k
  if (row_is_big(tab, row))
2045
0
    return 0;
2046
15.4k
  return isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
2047
15.4k
          tab->n_col - tab->n_dead) == -1;
2048
15.4k
}
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
88.4k
{
2057
88.4k
  struct isl_tab_var *var;
2058
88.4k
  int r;
2059
88.4k
2060
88.4k
  if (!tab)
2061
0
    return -1;
2062
88.4k
  r = isl_tab_add_row(tab, eq);
2063
88.4k
  if (r < 0)
2064
0
    return -1;
2065
88.4k
2066
88.4k
  var = &tab->con[r];
2067
88.4k
  r = var->index;
2068
88.4k
  if (row_is_manifestly_zero(tab, r)) {
2069
802
    var->is_zero = 1;
2070
802
    if (isl_tab_mark_redundant(tab, r) < 0)
2071
0
      return -1;
2072
802
    return 0;
2073
802
  }
2074
87.6k
2075
87.6k
  if (isl_int_is_neg(tab->mat->row[r][1])) {
2076
33.8k
    isl_seq_neg(tab->mat->row[r] + 1, tab->mat->row[r] + 1,
2077
33.8k
          1 + tab->n_col);
2078
33.8k
    var->negated = 1;
2079
33.8k
  }
2080
87.6k
  var->is_nonneg = 1;
2081
87.6k
  if (to_col(tab, var) < 0)
2082
0
    return -1;
2083
87.6k
  var->is_nonneg = 0;
2084
87.6k
  if (isl_tab_kill_col(tab, var->index) < 0)
2085
0
    return -1;
2086
87.6k
2087
87.6k
  return 0;
2088
87.6k
}
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.20k
{
2098
3.20k
  int r;
2099
3.20k
  isl_int *row;
2100
3.20k
2101
3.20k
  r = isl_tab_allocate_con(tab);
2102
3.20k
  if (r < 0)
2103
0
    return -1;
2104
3.20k
2105
3.20k
  row = tab->mat->row[tab->con[r].index];
2106
3.20k
  isl_seq_clr(row + 1, 1 + tab->M + tab->n_col);
2107
3.20k
  isl_int_set_si(row[0], 1);
2108
3.20k
2109
3.20k
  return r;
2110
3.20k
}
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
15.2k
{
2121
15.2k
  struct isl_tab_undo *snap = NULL;
2122
15.2k
  struct isl_tab_var *var;
2123
15.2k
  int r;
2124
15.2k
  int row;
2125
15.2k
  int sgn;
2126
15.2k
  isl_int cst;
2127
15.2k
2128
15.2k
  if (!tab)
2129
0
    return -1;
2130
15.2k
  isl_assert(tab->mat->ctx, !tab->M, return -1);
2131
15.2k
2132
15.2k
  if (tab->need_undo)
2133
14.7k
    snap = isl_tab_snap(tab);
2134
15.2k
2135
15.2k
  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
15.2k
  r = isl_tab_add_row(tab, eq);
2141
15.2k
  if (tab->cone) {
2142
1.12k
    isl_int_swap(eq[0], cst);
2143
1.12k
    isl_int_clear(cst);
2144
1.12k
  }
2145
15.2k
  if (r < 0)
2146
0
    return -1;
2147
15.2k
2148
15.2k
  var = &tab->con[r];
2149
15.2k
  row = var->index;
2150
15.2k
  if (row_is_manifestly_zero(tab, row)) {
2151
10.9k
    if (snap)
2152
10.9k
      return isl_tab_rollback(tab, snap);
2153
50
    return drop_row(tab, row);
2154
50
  }
2155
4.31k
2156
4.31k
  if (tab->bmap) {
2157
3.20k
    tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);
2158
3.20k
    if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
2159
0
      return -1;
2160
3.20k
    isl_seq_neg(eq, eq, 1 + tab->n_var);
2161
3.20k
    tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);
2162
3.20k
    isl_seq_neg(eq, eq, 1 + tab->n_var);
2163
3.20k
    if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
2164
0
      return -1;
2165
3.20k
    if (!tab->bmap)
2166
0
      return -1;
2167
3.20k
    if (add_zero_row(tab) < 0)
2168
0
      return -1;
2169
4.31k
  }
2170
4.31k
2171
4.31k
  sgn = isl_int_sgn(tab->mat->row[row][1]);
2172
4.31k
2173
4.31k
  if (sgn > 0) {
2174
263
    isl_seq_neg(tab->mat->row[row] + 1, tab->mat->row[row] + 1,
2175
263
          1 + tab->n_col);
2176
263
    var->negated = 1;
2177
263
    sgn = -1;
2178
263
  }
2179
4.31k
2180
4.31k
  if (sgn < 0) {
2181
3.01k
    sgn = sign_of_max(tab, var);
2182
3.01k
    if (sgn < -1)
2183
0
      return -1;
2184
3.01k
    if (sgn < 0) {
2185
0
      if (isl_tab_mark_empty(tab) < 0)
2186
0
        return -1;
2187
0
      return 0;
2188
0
    }
2189
3.01k
  }
2190
4.31k
2191
4.31k
  var->is_nonneg = 1;
2192
4.31k
  if (to_col(tab, var) < 0)
2193
0
    return -1;
2194
4.31k
  var->is_nonneg = 0;
2195
4.31k
  if (isl_tab_kill_col(tab, var->index) < 0)
2196
0
    return -1;
2197
4.31k
2198
4.31k
  return 0;
2199
4.31k
}
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.63k
{
2213
3.63k
  unsigned total;
2214
3.63k
  unsigned div_pos;
2215
3.63k
  struct isl_vec *ineq;
2216
3.63k
2217
3.63k
  if (!bmap)
2218
0
    return NULL;
2219
3.63k
2220
3.63k
  total = isl_basic_map_total_dim(bmap);
2221
3.63k
  div_pos = 1 + total - bmap->n_div + div;
2222
3.63k
2223
3.63k
  ineq = isl_vec_alloc(bmap->ctx, 1 + total);
2224
3.63k
  if (!ineq)
2225
0
    return NULL;
2226
3.63k
2227
3.63k
  isl_seq_cpy(ineq->el, bmap->div[div] + 1, 1 + total);
2228
3.63k
  isl_int_neg(ineq->el[div_pos], bmap->div[div][0]);
2229
3.63k
  return ineq;
2230
3.63k
}
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.63k
{
2250
3.63k
  unsigned total;
2251
3.63k
  unsigned div_pos;
2252
3.63k
  struct isl_vec *ineq;
2253
3.63k
2254
3.63k
  total = isl_basic_map_total_dim(tab->bmap);
2255
3.63k
  div_pos = 1 + total - tab->bmap->n_div + div;
2256
3.63k
2257
3.63k
  ineq = ineq_for_div(tab->bmap, div);
2258
3.63k
  if (!ineq)
2259
0
    goto error;
2260
3.63k
2261
3.63k
  if (add_ineq) {
2262
626
    if (add_ineq(user, ineq->el) < 0)
2263
0
      goto error;
2264
3.00k
  } else {
2265
3.00k
    if (isl_tab_add_ineq(tab, ineq->el) < 0)
2266
0
      goto error;
2267
3.63k
  }
2268
3.63k
2269
3.63k
  isl_seq_neg(ineq->el, tab->bmap->div[div] + 1, 1 + total);
2270
3.63k
  isl_int_set(ineq->el[div_pos], tab->bmap->div[div][0]);
2271
3.63k
  isl_int_add(ineq->el[0], ineq->el[0], ineq->el[div_pos]);
2272
3.63k
  isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
2273
3.63k
2274
3.63k
  if (add_ineq) {
2275
626
    if (add_ineq(user, ineq->el) < 0)
2276
0
      goto error;
2277
3.00k
  } else {
2278
3.00k
    if (isl_tab_add_ineq(tab, ineq->el) < 0)
2279
0
      goto error;
2280
3.63k
  }
2281
3.63k
2282
3.63k
  isl_vec_free(ineq);
2283
3.63k
2284
3.63k
  return isl_stat_ok;
2285
0
error:
2286
0
  isl_vec_free(ineq);
2287
0
  return isl_stat_error;
2288
3.63k
}
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.63k
{
2298
3.63k
  int i;
2299
3.63k
2300
3.63k
  if (tab->M)
2301
0
    return 1;
2302
3.63k
2303
3.63k
  if (isl_int_is_neg(div->el[1]))
2304
3.63k
    
return 0683
;
2305
2.94k
2306
9.13k
  
for (i = 0; 2.94k
i < tab->n_var;
++i6.18k
) {
2307
8.60k
    if (isl_int_is_neg(div->el[2 + i]))
2308
8.60k
      
return 0279
;
2309
8.32k
    if (isl_int_is_zero(div->el[2 + i]))
2310
8.32k
      
continue5.49k
;
2311
2.83k
    if (!tab->var[i].is_nonneg)
2312
2.14k
      return 0;
2313
2.83k
  }
2314
2.94k
2315
2.94k
  
return 1525
;
2316
2.94k
}
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.63k
{
2332
3.63k
  int r;
2333
3.63k
  int nonneg;
2334
3.63k
  int n_div, o_div;
2335
3.63k
2336
3.63k
  if (!tab || !div)
2337
0
    return -1;
2338
3.63k
2339
3.63k
  if (div->size != 1 + 1 + tab->n_var)
2340
3.63k
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_invalid,
2341
3.63k
      "unexpected size", return -1);
2342
3.63k
2343
3.63k
  isl_assert(tab->mat->ctx, tab->bmap, return -1);
2344
3.63k
  n_div = isl_basic_map_dim(tab->bmap, isl_dim_div);
2345
3.63k
  o_div = tab->n_var - n_div;
2346
3.63k
  if (pos < o_div || pos > tab->n_var)
2347
3.63k
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_invalid,
2348
3.63k
      "invalid position", return -1);
2349
3.63k
2350
3.63k
  nonneg = div_is_nonneg(tab, div);
2351
3.63k
2352
3.63k
  if (isl_tab_extend_cons(tab, 3) < 0)
2353
0
    return -1;
2354
3.63k
  if (isl_tab_extend_vars(tab, 1) < 0)
2355
0
    return -1;
2356
3.63k
  r = isl_tab_insert_var(tab, pos);
2357
3.63k
  if (r < 0)
2358
0
    return -1;
2359
3.63k
2360
3.63k
  if (nonneg)
2361
525
    tab->var[r].is_nonneg = 1;
2362
3.63k
2363
3.63k
  tab->bmap = isl_basic_map_insert_div(tab->bmap, pos - o_div, div);
2364
3.63k
  if (!tab->bmap)
2365
0
    return -1;
2366
3.63k
  if (isl_tab_push_var(tab, isl_tab_undo_bmap_div, &tab->var[r]) < 0)
2367
0
    return -1;
2368
3.63k
2369
3.63k
  if (add_div_constraints(tab, pos - o_div, add_ineq, user) < 0)
2370
0
    return -1;
2371
3.63k
2372
3.63k
  return r;
2373
3.63k
}
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.00k
{
2380
3.00k
  if (!tab)
2381
0
    return -1;
2382
3.00k
  return isl_tab_insert_div(tab, tab->n_var, div, NULL, NULL);
2383
3.00k
}
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
730k
{
2393
730k
  int i;
2394
730k
  struct isl_tab *tab;
2395
730k
2396
730k
  if (!bmap)
2397
0
    return NULL;
2398
730k
  tab = isl_tab_alloc(bmap->ctx,
2399
730k
          isl_basic_map_total_dim(bmap) + bmap->n_ineq + 1,
2400
730k
          isl_basic_map_total_dim(bmap), 0);
2401
730k
  if (!tab)
2402
0
    return NULL;
2403
730k
  tab->preserve = track;
2404
730k
  tab->rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
2405
730k
  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
998k
  
for (i = 0; 730k
i < bmap->n_eq;
++i267k
) {
2411
267k
    tab = add_eq(tab, bmap->eq[i]);
2412
267k
    if (!tab)
2413
0
      return tab;
2414
267k
  }
2415
4.44M
  
for (i = 0; 730k
i < bmap->n_ineq;
++i3.71M
) {
2416
3.71M
    if (isl_tab_add_ineq(tab, bmap->ineq[i]) < 0)
2417
0
      goto error;
2418
3.71M
    if (tab->empty)
2419
5.47k
      goto done;
2420
3.71M
  }
2421
730k
done:
2422
730k
  if (track && 
isl_tab_track_bmap(tab, isl_basic_map_copy(bmap)) < 0193k
)
2423
0
    goto error;
2424
730k
  return tab;
2425
0
error:
2426
0
  isl_tab_free(tab);
2427
0
  return NULL;
2428
730k
}
2429
2430
__isl_give struct isl_tab *isl_tab_from_basic_set(
2431
  __isl_keep isl_basic_set *bset, int track)
2432
279k
{
2433
279k
  return isl_tab_from_basic_map(bset, track);
2434
279k
}
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.61k
{
2441
3.61k
  isl_int cst;
2442
3.61k
  int i;
2443
3.61k
  struct isl_tab *tab;
2444
3.61k
  unsigned offset = 0;
2445
3.61k
2446
3.61k
  if (!bset)
2447
0
    return NULL;
2448
3.61k
  if (parametric)
2449
2.84k
    offset = isl_basic_set_dim(bset, isl_dim_param);
2450
3.61k
  tab = isl_tab_alloc(bset->ctx, bset->n_eq + bset->n_ineq,
2451
3.61k
        isl_basic_set_total_dim(bset) - offset, 0);
2452
3.61k
  if (!tab)
2453
0
    return NULL;
2454
3.61k
  tab->rational = ISL_F_ISSET(bset, ISL_BASIC_SET_RATIONAL);
2455
3.61k
  tab->cone = 1;
2456
3.61k
2457
3.61k
  isl_int_init(cst);
2458
3.61k
  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
15.0k
  
for (i = 0; 3.61k
i < bset->n_ineq;
++i11.4k
) {
2471
11.4k
    int r;
2472
11.4k
    isl_int_swap(bset->ineq[i][offset], cst);
2473
11.4k
    r = isl_tab_add_row(tab, bset->ineq[i] + offset);
2474
11.4k
    isl_int_swap(bset->ineq[i][offset], cst);
2475
11.4k
    if (r < 0)
2476
0
      goto error;
2477
11.4k
    tab->con[r].is_nonneg = 1;
2478
11.4k
    if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
2479
0
      goto error;
2480
11.4k
  }
2481
3.61k
done:
2482
3.61k
  isl_int_clear(cst);
2483
3.61k
  return tab;
2484
0
error:
2485
0
  isl_int_clear(cst);
2486
0
  isl_tab_free(tab);
2487
0
  return NULL;
2488
3.61k
}
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.84k
{
2495
2.84k
  int i;
2496
2.84k
2497
2.84k
  if (!tab)
2498
0
    return isl_bool_error;
2499
2.84k
  if (tab->empty)
2500
0
    return isl_bool_true;
2501
2.84k
  if (tab->n_dead == tab->n_col)
2502
737
    return isl_bool_true;
2503
2.11k
2504
3.24k
  
for (;;)2.11k
{
2505
3.53k
    for (i = tab->n_redundant; i < tab->n_row; 
++i295
) {
2506
3.53k
      struct isl_tab_var *var;
2507
3.53k
      int sgn;
2508
3.53k
      var = isl_tab_var_from_row(tab, i);
2509
3.53k
      if (!var->is_nonneg)
2510
295
        continue;
2511
3.23k
      sgn = sign_of_max(tab, var);
2512
3.23k
      if (sgn < -1)
2513
0
        return isl_bool_error;
2514
3.23k
      if (sgn != 0)
2515
276
        return isl_bool_false;
2516
2.96k
      if (close_row(tab, var, 0) < 0)
2517
0
        return isl_bool_error;
2518
2.96k
      break;
2519
2.96k
    }
2520
3.24k
    
if (2.96k
tab->n_dead == tab->n_col2.96k
)
2521
1.83k
      return isl_bool_true;
2522
1.13k
    if (i == tab->n_row)
2523
3
      return isl_bool_false;
2524
1.13k
  }
2525
2.11k
}
2526
2527
int isl_tab_sample_is_integer(struct isl_tab *tab)
2528
437k
{
2529
437k
  int i;
2530
437k
2531
437k
  if (!tab)
2532
0
    return -1;
2533
437k
2534
2.11M
  
for (i = 0; 437k
i < tab->n_var;
++i1.67M
) {
2535
1.77M
    int row;
2536
1.77M
    if (!tab->var[i].is_row)
2537
511k
      continue;
2538
1.26M
    row = tab->var[i].index;
2539
1.26M
    if (!isl_int_is_divisible_by(tab->mat->row[row][1],
2540
1.26M
            tab->mat->row[row][0]))
2541
1.26M
      
return 098.1k
;
2542
1.26M
  }
2543
437k
  
return 1339k
;
2544
437k
}
2545
2546
static struct isl_vec *extract_integer_sample(struct isl_tab *tab)
2547
205k
{
2548
205k
  int i;
2549
205k
  struct isl_vec *vec;
2550
205k
2551
205k
  vec = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
2552
205k
  if (!vec)
2553
0
    return NULL;
2554
205k
2555
205k
  isl_int_set_si(vec->block.data[0], 1);
2556
1.30M
  for (i = 0; i < tab->n_var; 
++i1.09M
) {
2557
1.09M
    if (!tab->var[i].is_row)
2558
1.09M
      
isl_int_set_si397k
(vec->block.data[1 + i], 0);
2559
1.09M
    else {
2560
698k
      int row = tab->var[i].index;
2561
698k
      isl_int_divexact(vec->block.data[1 + i],
2562
698k
        tab->mat->row[row][1], tab->mat->row[row][0]);
2563
698k
    }
2564
1.09M
  }
2565
205k
2566
205k
  return vec;
2567
205k
}
2568
2569
struct isl_vec *isl_tab_get_sample_value(struct isl_tab *tab)
2570
259k
{
2571
259k
  int i;
2572
259k
  struct isl_vec *vec;
2573
259k
  isl_int m;
2574
259k
2575
259k
  if (!tab)
2576
0
    return NULL;
2577
259k
2578
259k
  vec = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
2579
259k
  if (!vec)
2580
0
    return NULL;
2581
259k
2582
259k
  isl_int_init(m);
2583
259k
2584
259k
  isl_int_set_si(vec->block.data[0], 1);
2585
1.14M
  for (i = 0; i < tab->n_var; 
++i882k
) {
2586
882k
    int row;
2587
882k
    if (!tab->var[i].is_row) {
2588
417k
      isl_int_set_si(vec->block.data[1 + i], 0);
2589
417k
      continue;
2590
417k
    }
2591
465k
    row = tab->var[i].index;
2592
465k
    isl_int_gcd(m, vec->block.data[0], tab->mat->row[row][0]);
2593
465k
    isl_int_divexact(m, tab->mat->row[row][0], m);
2594
465k
    isl_seq_scale(vec->block.data, vec->block.data, m, 1 + i);
2595
465k
    isl_int_divexact(m, vec->block.data[0], tab->mat->row[row][0]);
2596
465k
    isl_int_mul(vec->block.data[1 + i], m, tab->mat->row[row][1]);
2597
465k
  }
2598
259k
  vec = isl_vec_normalize(vec);
2599
259k
2600
259k
  isl_int_clear(m);
2601
259k
  return vec;
2602
259k
}
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
248k
{
2610
248k
  if (!var->is_row)
2611
248k
    
isl_int_set_si3.49k
(*v, 0);
2612
248k
  else 
if (245k
sgn > 0245k
)
2613
245k
    
isl_int_cdiv_q242k
(*v, tab->mat->row[var->index][1],
2614
245k
           tab->mat->row[var->index][0]);
2615
245k
  else
2616
245k
    
isl_int_fdiv_q2.55k
(*v, tab->mat->row[var->index][1],
2617
248k
           tab->mat->row[var->index][0]);
2618
248k
}
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
431k
{
2631
431k
  int i;
2632
431k
  unsigned n_eq;
2633
431k
2634
431k
  if (!bmap)
2635
0
    return NULL;
2636
431k
  if (!tab)
2637
0
    return bmap;
2638
431k
2639
431k
  n_eq = tab->n_eq;
2640
431k
  if (tab->empty)
2641
6.83k
    bmap = isl_basic_map_set_to_empty(bmap);
2642
424k
  else
2643
2.95M
    
for (i = bmap->n_ineq - 1; 424k
i >= 0;
--i2.53M
) {
2644
2.53M
      if (isl_tab_is_equality(tab, n_eq + i))
2645
1.06M
        isl_basic_map_inequality_to_equality(bmap, i);
2646
1.46M
      else if (isl_tab_is_redundant(tab, n_eq + i))
2647
191k
        isl_basic_map_drop_inequality(bmap, i);
2648
2.53M
    }
2649
431k
  if (bmap->n_eq != n_eq)
2650
160k
    bmap = isl_basic_map_gauss(bmap, NULL);
2651
431k
  if (!tab->rational &&
2652
431k
      
bmap397k
&&
!bmap->sample397k
&&
isl_tab_sample_is_integer(tab)222k
)
2653
205k
    bmap->sample = extract_integer_sample(tab);
2654
431k
  return bmap;
2655
431k
}
2656
2657
struct isl_basic_set *isl_basic_set_update_from_tab(struct isl_basic_set *bset,
2658
  struct isl_tab *tab)
2659
39.7k
{
2660
39.7k
  return bset_from_bmap(isl_basic_map_update_from_tab(bset_to_bmap(bset),
2661
39.7k
                tab));
2662
39.7k
}
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
13.4k
{
2669
13.4k
  if (!tab->con[r].is_row)
2670
13.4k
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_internal,
2671
13.4k
      "row unexpectedly moved to column",
2672
13.4k
      return isl_stat_error);
2673
13.4k
  if (r + 1 != tab->n_con)
2674
13.4k
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_internal,
2675
13.4k
      "additional constraints added", return isl_stat_error);
2676
13.4k
  if (drop_row(tab, tab->con[r].index) < 0)
2677
0
    return isl_stat_error;
2678
13.4k
2679
13.4k
  return isl_stat_ok;
2680
13.4k
}
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
13.4k
{
2694
13.4k
  unsigned r;
2695
13.4k
  isl_int *row;
2696
13.4k
  int sgn;
2697
13.4k
  unsigned off = 2 + tab->M;
2698
13.4k
2699
13.4k
  if (var->is_zero)
2700
0
    return isl_stat_ok;
2701
13.4k
  if (var->is_redundant || !var->is_nonneg)
2702
13.4k
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_invalid,
2703
13.4k
      "expecting non-redundant non-negative variable",
2704
13.4k
      return isl_stat_error);
2705
13.4k
2706
13.4k
  if (isl_tab_extend_cons(tab, 1) < 0)
2707
0
    return isl_stat_error;
2708
13.4k
2709
13.4k
  r = tab->n_con;
2710
13.4k
  tab->con[r].index = tab->n_row;
2711
13.4k
  tab->con[r].is_row = 1;
2712
13.4k
  tab->con[r].is_nonneg = 0;
2713
13.4k
  tab->con[r].is_zero = 0;
2714
13.4k
  tab->con[r].is_redundant = 0;
2715
13.4k
  tab->con[r].frozen = 0;
2716
13.4k
  tab->con[r].negated = 0;
2717
13.4k
  tab->row_var[tab->n_row] = ~r;
2718
13.4k
  row = tab->mat->row[tab->n_row];
2719
13.4k
2720
13.4k
  if (var->is_row) {
2721
1.37k
    isl_int_set(row[0], tab->mat->row[var->index][0]);
2722
1.37k
    isl_seq_neg(row + 1,
2723
1.37k
          tab->mat->row[var->index] + 1, 1 + tab->n_col);
2724
12.1k
  } else {
2725
12.1k
    isl_int_set_si(row[0], 1);
2726
12.1k
    isl_seq_clr(row + 1, 1 + tab->n_col);
2727
12.1k
    isl_int_set_si(row[off + var->index], -1);
2728
12.1k
  }
2729
13.4k
2730
13.4k
  tab->n_row++;
2731
13.4k
  tab->n_con++;
2732
13.4k
2733
13.4k
  sgn = sign_of_max(tab, &tab->con[r]);
2734
13.4k
  if (sgn < -1)
2735
0
    return isl_stat_error;
2736
13.4k
  if (sgn < 0) {
2737
51
    if (drop_last_con_in_row(tab, r) < 0)
2738
0
      return isl_stat_error;
2739
51
    if (isl_tab_mark_empty(tab) < 0)
2740
0
      return isl_stat_error;
2741
51
    return isl_stat_ok;
2742
51
  }
2743
13.4k
  tab->con[r].is_nonneg = 1;
2744
13.4k
  /* sgn == 0 */
2745
13.4k
  if (close_row(tab, &tab->con[r], 1) < 0)
2746
0
    return isl_stat_error;
2747
13.4k
  if (drop_last_con_in_row(tab, r) < 0)
2748
0
    return isl_stat_error;
2749
13.4k
2750
13.4k
  return isl_stat_ok;
2751
13.4k
}
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.80k
{
2772
1.80k
  struct isl_tab_var *var;
2773
1.80k
2774
1.80k
  if (!tab)
2775
0
    return -1;
2776
1.80k
2777
1.80k
  var = &tab->con[con];
2778
1.80k
2779
1.80k
  if (var->is_row && 
(57
var->index < 057
||
var->index < tab->n_redundant57
))
2780
1.80k
    
isl_die0
(tab->mat->ctx, isl_error_invalid,
2781
1.80k
      "cannot relax redundant constraint", return -1);
2782
1.80k
  if (!var->is_row && 
(1.75k
var->index < 01.75k
||
var->index < tab->n_dead1.75k
))
2783
1.80k
    
isl_die0
(tab->mat->ctx, isl_error_invalid,
2784
1.80k
      "cannot relax dead constraint", return -1);
2785
1.80k
2786
1.80k
  if (!var->is_row && 
!max_is_manifestly_unbounded(tab, var)1.75k
)
2787
456
    if (to_row(tab, var, 1) < 0)
2788
0
      return -1;
2789
1.80k
  if (!var->is_row && 
!min_is_manifestly_unbounded(tab, var)1.29k
)
2790
18
    if (to_row(tab, var, -1) < 0)
2791
0
      return -1;
2792
1.80k
2793
1.80k
  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.27k
  } else {
2799
1.27k
    int i;
2800
1.27k
    unsigned off = 2 + tab->M;
2801
1.27k
2802
7.38k
    for (i = 0; i < tab->n_row; 
++i6.10k
) {
2803
6.10k
      if (isl_int_is_zero(tab->mat->row[i][off + var->index]))
2804
6.10k
        
continue4.71k
;
2805
1.38k
      isl_int_sub(tab->mat->row[i][1], tab->mat->row[i][1],
2806
1.38k
          tab->mat->row[i][off + var->index]);
2807
1.38k
    }
2808
1.27k
2809
1.27k
  }
2810
1.80k
2811
1.80k
  if (isl_tab_push_var(tab, isl_tab_undo_relax, var) < 0)
2812
0
    return -1;
2813
1.80k
2814
1.80k
  return 0;
2815
1.80k
}
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
523
{
2882
523
  struct isl_tab_var *var;
2883
523
2884
523
  if (!tab)
2885
0
    return -1;
2886
523
2887
523
  var = &tab->con[con];
2888
523
  if (!var->is_nonneg)
2889
0
    return 0;
2890
523
2891
523
  var->is_nonneg = 0;
2892
523
  if (isl_tab_push_var(tab, isl_tab_undo_unrestrict, var) < 0)
2893
0
    return -1;
2894
523
2895
523
  return 0;
2896
523
}
2897
2898
int isl_tab_select_facet(struct isl_tab *tab, int con)
2899
12.9k
{
2900
12.9k
  if (!tab)
2901
0
    return -1;
2902
12.9k
2903
12.9k
  return cut_to_hyperplane(tab, &tab->con[con]);
2904
12.9k
}
2905
2906
static int may_be_equality(struct isl_tab *tab, int row)
2907
8.23M
{
2908
8.23M
  return tab->rational ? 
isl_int_is_zero36.3k
(tab->mat->row[row][1])
2909
8.23M
           : 
isl_int_lt8.19M
(tab->mat->row[row][1],
2910
8.23M
              tab->mat->row[row][0]);
2911
8.23M
}
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.02M
{
2926
2.02M
  int i;
2927
2.02M
  struct isl_tab_var *var;
2928
2.02M
2929
16.9M
  for (i = tab->n_con - 1; i >= 0; 
--i14.9M
) {
2930
16.7M
    var = &tab->con[i];
2931
16.7M
    if (var->index < 0)
2932
6.55M
      continue;
2933
10.2M
    if (var->is_row && 
var->index < tab->n_redundant8.06M
)
2934
608k
      continue;
2935
9.63M
    if (!var->is_row && 
var->index < tab->n_dead2.17M
)
2936
1.80k
      continue;
2937
9.63M
    if (var->marked)
2938
1.85M
      return var;
2939
9.63M
  }
2940
2.02M
2941
2.02M
  
return NULL166k
;
2942
2.02M
}
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
429k
{
2963
429k
  int i;
2964
429k
  unsigned n_marked;
2965
429k
2966
429k
  if (!tab)
2967
0
    return -1;
2968
429k
  if (tab->empty)
2969
3.84k
    return 0;
2970
425k
  if (tab->n_dead == tab->n_col)
2971
18.1k
    return 0;
2972
407k
2973
407k
  n_marked = 0;
2974
3.03M
  for (i = tab->n_redundant; i < tab->n_row; 
++i2.62M
) {
2975
2.62M
    struct isl_tab_var *var = isl_tab_var_from_row(tab, i);
2976
2.62M
    var->marked = !var->frozen && 
var->is_nonneg2.60M
&&
2977
2.62M
      
may_be_equality(tab, i)2.38M
;
2978
2.62M
    if (var->marked)
2979
1.76M
      n_marked++;
2980
2.62M
  }
2981
2.32M
  for (i = tab->n_dead; i < tab->n_col; 
++i1.91M
) {
2982
1.91M
    struct isl_tab_var *var = var_from_col(tab, i);
2983
1.91M
    var->marked = !var->frozen && 
var->is_nonneg1.90M
;
2984
1.91M
    if (var->marked)
2985
176k
      n_marked++;
2986
1.91M
  }
2987
1.68M
  while (n_marked) {
2988
1.44M
    struct isl_tab_var *var;
2989
1.44M
    int sgn;
2990
1.44M
    var = select_marked(tab);
2991
1.44M
    if (!var)
2992
158k
      break;
2993
1.28M
    var->marked = 0;
2994
1.28M
    n_marked--;
2995
1.28M
    sgn = sign_of_max(tab, var);
2996
1.28M
    if (sgn < 0)
2997
0
      return -1;
2998
1.28M
    if (sgn == 0) {
2999
566k
      if (close_row(tab, var, 0) < 0)
3000
0
        return -1;
3001
715k
    } else if (!tab->rational && 
!at_least_one(tab, var)692k
) {
3002
560
      if (cut_to_hyperplane(tab, var) < 0)
3003
0
        return -1;
3004
560
      return isl_tab_detect_implicit_equalities(tab);
3005
560
    }
3006
11.6M
    
for (i = tab->n_redundant; 1.28M
i < tab->n_row;
++i10.3M
) {
3007
10.3M
      var = isl_tab_var_from_row(tab, i);
3008
10.3M
      if (!var->marked)
3009
4.47M
        continue;
3010
5.84M
      if (may_be_equality(tab, i))
3011
5.76M
        continue;
3012
80.3k
      var->marked = 0;
3013
80.3k
      n_marked--;
3014
80.3k
    }
3015
1.28M
  }
3016
407k
3017
407k
  
return 0406k
;
3018
407k
}
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
58.8k
{
3098
58.8k
  int i;
3099
58.8k
3100
58.8k
  if (!tab || !bmap)
3101
0
    return isl_basic_map_free(bmap);
3102
58.8k
  if (tab->empty)
3103
80
    return bmap;
3104
58.7k
3105
169k
  
for (i = bmap->n_ineq - 1; 58.7k
i >= 0;
--i110k
) {
3106
110k
    if (!isl_tab_is_equality(tab, bmap->n_eq + i))
3107
109k
      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
58.7k
3117
58.7k
  return bmap;
3118
58.7k
}
3119
3120
static int con_is_redundant(struct isl_tab *tab, struct isl_tab_var *var)
3121
854k
{
3122
854k
  if (!tab)
3123
0
    return -1;
3124
854k
  if (tab->rational) {
3125
70.6k
    int sgn = sign_of_min(tab, var);
3126
70.6k
    if (sgn < -1)
3127
0
      return -1;
3128
70.6k
    return sgn >= 0;
3129
784k
  } else {
3130
784k
    int irred = isl_tab_min_at_most_neg_one(tab, var);
3131
784k
    if (irred < 0)
3132
0
      return -1;
3133
784k
    return !irred;
3134
784k
  }
3135
854k
}
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
203k
{
3152
203k
  int i;
3153
203k
  unsigned n_marked;
3154
203k
3155
203k
  if (!tab)
3156
0
    return -1;
3157
203k
  if (tab->empty)
3158
2.84k
    return 0;
3159
200k
  if (tab->n_redundant == tab->n_row)
3160
5.33k
    return 0;
3161
195k
3162
195k
  n_marked = 0;
3163
1.43M
  for (i = tab->n_redundant; i < tab->n_row; 
++i1.23M
) {
3164
1.23M
    struct isl_tab_var *var = isl_tab_var_from_row(tab, i);
3165
1.23M
    var->marked = !var->frozen && 
var->is_nonneg1.09M
;
3166
1.23M
    if (var->marked)
3167
554k
      n_marked++;
3168
1.23M
  }
3169
1.08M
  for (i = tab->n_dead; i < tab->n_col; 
++i885k
) {
3170
885k
    struct isl_tab_var *var = var_from_col(tab, i);
3171
885k
    var->marked = !var->frozen && 
var->is_nonneg828k
&&
3172
885k
      
!min_is_manifestly_unbounded(tab, var)371k
;
3173
885k
    if (var->marked)
3174
116k
      n_marked++;
3175
885k
  }
3176
768k
  while (n_marked) {
3177
580k
    struct isl_tab_var *var;
3178
580k
    int red;
3179
580k
    var = select_marked(tab);
3180
580k
    if (!var)
3181
7.19k
      break;
3182
573k
    var->marked = 0;
3183
573k
    n_marked--;
3184
573k
    red = con_is_redundant(tab, var);
3185
573k
    if (red < 0)
3186
0
      return -1;
3187
573k
    if (red && 
!var->is_redundant127k
)
3188
23.3k
      if (isl_tab_mark_redundant(tab, var->index) < 0)
3189
0
        return -1;
3190
8.28M
    
for (i = tab->n_dead; 573k
i < tab->n_col;
++i7.70M
) {
3191
7.70M
      var = var_from_col(tab, i);
3192
7.70M
      if (!var->marked)
3193
7.34M
        continue;
3194
357k
      if (!min_is_manifestly_unbounded(tab, var))
3195
269k
        continue;
3196
88.3k
      var->marked = 0;
3197
88.3k
      n_marked--;
3198
88.3k
    }
3199
573k
  }
3200
195k
3201
195k
  return 0;
3202
195k
}
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
1.06M
    return 1;
3213
1.59M
  if (tab->con[con].is_redundant)
3214
193k
    return 0;
3215
1.40M
  if (!tab->con[con].is_row)
3216
786k
    return tab->con[con].index < tab->n_dead;
3217
617k
3218
617k
  row = tab->con[con].index;
3219
617k
3220
617k
  off = 2 + tab->M;
3221
617k
  return isl_int_is_zero(tab->mat->row[row][1]) &&
3222
617k
    
!row_is_big(tab, row)56.9k
&&
3223
617k
    isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
3224
56.9k
          tab->n_col - tab->n_dead) == -1;
3225
617k
}
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
311k
{
3241
311k
  int r;
3242
311k
  enum isl_lp_result res = isl_lp_ok;
3243
311k
  struct isl_tab_var *var;
3244
311k
  struct isl_tab_undo *snap;
3245
311k
3246
311k
  if (!tab)
3247
0
    return isl_lp_error;
3248
311k
3249
311k
  if (tab->empty)
3250
20
    return isl_lp_empty;
3251
311k
3252
311k
  snap = isl_tab_snap(tab);
3253
311k
  r = isl_tab_add_row(tab, f);
3254
311k
  if (r < 0)
3255
0
    return isl_lp_error;
3256
311k
  var = &tab->con[r];
3257
680k
  for (;;) {
3258
680k
    int row, col;
3259
680k
    find_pivot(tab, var, var, -1, &row, &col);
3260
680k
    if (row == var->index) {
3261
10.0k
      res = isl_lp_unbounded;
3262
10.0k
      break;
3263
10.0k
    }
3264
670k
    if (row == -1)
3265
301k
      break;
3266
369k
    if (isl_tab_pivot(tab, row, col) < 0)
3267
0
      return isl_lp_error;
3268
369k
  }
3269
311k
  isl_int_mul(tab->mat->row[var->index][0],
3270
311k
        tab->mat->row[var->index][0], denom);
3271
311k
  if (ISL_FL_ISSET(flags, ISL_TAB_SAVE_DUAL)) {
3272
26.8k
    int i;
3273
26.8k
3274
26.8k
    isl_vec_free(tab->dual);
3275
26.8k
    tab->dual = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_con);
3276
26.8k
    if (!tab->dual)
3277
0
      return isl_lp_error;
3278
26.8k
    isl_int_set(tab->dual->el[0], tab->mat->row[var->index][0]);
3279
1.05M
    for (i = 0; i < tab->n_con; 
++i1.02M
) {
3280
1.02M
      int pos;
3281
1.02M
      if (tab->con[i].is_row) {
3282
633k
        isl_int_set_si(tab->dual->el[1 + i], 0);
3283
633k
        continue;
3284
633k
      }
3285
393k
      pos = 2 + tab->M + tab->con[i].index;
3286
393k
      if (tab->con[i].negated)
3287
393k
        
isl_int_neg70.1k
(tab->dual->el[1 + i],
3288
393k
              tab->mat->row[var->index][pos]);
3289
393k
      else
3290
393k
        
isl_int_set322k
(tab->dual->el[1 + i],
3291
393k
              tab->mat->row[var->index][pos]);
3292
393k
    }
3293
26.8k
  }
3294
311k
  if (opt && 
res == isl_lp_ok310k
) {
3295
301k
    if (opt_denom) {
3296
59.9k
      isl_int_set(*opt, tab->mat->row[var->index][1]);
3297
59.9k
      isl_int_set(*opt_denom, tab->mat->row[var->index][0]);
3298
59.9k
    } else
3299
241k
      get_rounded_sample_value(tab, var, 1, opt);
3300
301k
  }
3301
311k
  if (isl_tab_rollback(tab, snap) < 0)
3302
0
    return isl_lp_error;
3303
311k
  return res;
3304
311k
}
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.89M
{
3316
1.89M
  if (!tab)
3317
0
    return -1;
3318
1.89M
  if (con < 0 || con >= tab->n_con)
3319
1.89M
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_invalid,
3320
1.89M
      "position out of bounds", return -1);
3321
1.89M
  if (tab->con[con].is_zero)
3322
164
    return 0;
3323
1.89M
  if (tab->con[con].is_redundant)
3324
305k
    return 1;
3325
1.59M
  return tab->con[con].is_row && 
tab->con[con].index < tab->n_redundant654k
;
3326
1.59M
}
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
6.00k
{
3339
6.00k
  unsigned off = 2 + tab->M;
3340
6.00k
  isl_mat *mat = tab->mat;
3341
6.00k
  int n;
3342
6.00k
  int row;
3343
6.00k
  int pos;
3344
6.00k
3345
6.00k
  if (!var->is_row)
3346
3.49k
    return isl_bool_false;
3347
2.50k
  row = var->index;
3348
2.50k
  if (row_is_big(tab, row))
3349
0
    return isl_bool_false;
3350
2.50k
  n = tab->n_col - tab->n_dead;
3351
2.50k
  pos = isl_seq_first_non_zero(mat->row[row] + off + tab->n_dead, n);
3352
2.50k
  if (pos != -1)
3353
2.32k
    return isl_bool_false;
3354
183
  if (value)
3355
183
    
isl_int_divexact0
(*value, mat->row[row][1], mat->row[row][0]);
3356
183
  return isl_bool_true;
3357
183
}
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
8.04k
{
3372
8.04k
  if (row_is_big(tab, var->index))
3373
0
    return 1;
3374
8.04k
  isl_int_mul(*tmp, tab->mat->row[var->index][0], target);
3375
8.04k
  if (sgn > 0)
3376
3.27k
    return isl_int_ge(tab->mat->row[var->index][1], *tmp);
3377
8.04k
  else
3378
8.04k
    
return 4.76k
isl_int_le4.76k
(tab->mat->row[var->index][1], *tmp);
3379
8.04k
}
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
7.51k
{
3400
7.51k
  int row, col;
3401
7.51k
3402
7.51k
  if (sgn < 0 && 
min_is_manifestly_unbounded(tab, var)5.81k
)
3403
2.45k
    return isl_bool_true;
3404
5.06k
  if (sgn > 0 && 
max_is_manifestly_unbounded(tab, var)1.69k
)
3405
0
    return isl_bool_true;
3406
5.06k
  if (to_row(tab, var, sgn) < 0)
3407
0
    return isl_bool_error;
3408
8.04k
  
while (5.06k
!reached(tab, var, sgn, target, tmp)) {
3409
7.17k
    find_pivot(tab, var, var, sgn, &row, &col);
3410
7.17k
    if (row == -1)
3411
1.92k
      return isl_bool_false;
3412
5.24k
    if (row == var->index)
3413
2.25k
      return isl_bool_true;
3414
2.98k
    if (isl_tab_pivot(tab, row, col) < 0)
3415
0
      return isl_bool_error;
3416
2.98k
  }
3417
5.06k
3418
5.06k
  
return isl_bool_true872
;
3419
5.06k
}
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.81k
{
3443
5.81k
  isl_bool reached;
3444
5.81k
  isl_vec *eq;
3445
5.81k
  int pos;
3446
5.81k
  isl_stat r;
3447
5.81k
3448
5.81k
  get_rounded_sample_value(tab, var, -1, target);
3449
5.81k
  isl_int_sub_ui(*target, *target, 1);
3450
5.81k
  reached = var_reaches(tab, var, -1, *target, tmp);
3451
5.81k
  if (reached < 0 || reached)
3452
4.12k
    return isl_bool_not(reached);
3453
1.69k
  get_rounded_sample_value(tab, var, 1, target);
3454
1.69k
  isl_int_add_ui(*target, *target, 1);
3455
1.69k
  reached = var_reaches(tab, var, 1, *target, tmp);
3456
1.69k
  if (reached < 0 || reached)
3457
1.45k
    return isl_bool_not(reached);
3458
235
  get_rounded_sample_value(tab, var, -1, tmp);
3459
235
  isl_int_sub_ui(*target, *target, 1);
3460
235
  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
235
3466
235
  if (isl_tab_extend_cons(tab, 1) < 0)
3467
0
    return isl_bool_error;
3468
235
  eq = isl_vec_alloc(isl_tab_get_ctx(tab), 1 + tab->n_var);
3469
235
  if (!eq)
3470
0
    return isl_bool_error;
3471
235
  pos = var - tab->var;
3472
235
  isl_seq_clr(eq->el + 1, tab->n_var);
3473
235
  isl_int_set_si(eq->el[1 + pos], -1);
3474
235
  isl_int_set(eq->el[0], *target);
3475
235
  r = isl_tab_add_eq(tab, eq->el);
3476
235
  isl_vec_free(eq);
3477
235
3478
235
  return r < 0 ? 
isl_bool_error0
: isl_bool_true;
3479
235
}
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
6.00k
{
3497
6.00k
  isl_int target, tmp;
3498
6.00k
  isl_bool is_cst;
3499
6.00k
3500
6.00k
  if (var->is_row && 
row_is_big(tab, var->index)2.50k
)
3501
0
    return isl_bool_false;
3502
6.00k
  is_cst = is_constant(tab, var, value);
3503
6.00k
  if (is_cst < 0 || is_cst)
3504
183
    return is_cst;
3505
5.81k
3506
5.81k
  if (!value)
3507
5.81k
    
isl_int_init2.58k
(target);
3508
5.81k
  isl_int_init(tmp);
3509
5.81k
3510
5.81k
  is_cst = detect_constant_with_tmp(tab, var,
3511
5.81k
              value ? 
value3.23k
:
&target2.58k
, &tmp);
3512
5.81k
3513
5.81k
  isl_int_clear(tmp);
3514
5.81k
  if (!value)
3515
5.81k
    
isl_int_clear2.58k
(target);
3516
5.81k
3517
5.81k
  return is_cst;
3518
5.81k
}
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
3.23k
{
3529
3.23k
  if (!tab)
3530
0
    return isl_bool_error;
3531
3.23k
  if (var < 0 || var >= tab->n_var)
3532
3.23k
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_invalid,
3533
3.23k
      "position out of bounds", return isl_bool_error);
3534
3.23k
  if (tab->rational)
3535
0
    return isl_bool_false;
3536
3.23k
3537
3.23k
  return get_constant(tab, &tab->var[var], value);
3538
3.23k
}
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
523
{
3548
523
  int i;
3549
523
3550
523
  if (!tab)
3551
0
    return isl_stat_error;
3552
523
  if (tab->rational)
3553
0
    return isl_stat_ok;
3554
523
3555
3.29k
  
for (i = 0; 523
i < tab->n_var;
++i2.77k
) {
3556
2.77k
    if (get_constant(tab, &tab->var[i], NULL) < 0)
3557
0
      return isl_stat_error;
3558
2.77k
  }
3559
523
3560
523
  return isl_stat_ok;
3561
523
}
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.12M
{
3568
1.12M
  if (!tab)
3569
0
    return NULL;
3570
1.12M
  tab->need_undo = 1;
3571
1.12M
  return tab->top;
3572
1.12M
}
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
132
{
3579
132
  if (!tab)
3580
0
    return isl_bool_error;
3581
132
3582
132
  return tab->need_undo;
3583
132
}
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
132
{
3592
132
  if (!tab)
3593
0
    return;
3594
132
3595
132
  free_undo(tab);
3596
132
  tab->need_undo = 0;
3597
132
}
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.21k
{
3605
1.21k
  unsigned off = 2 + tab->M;
3606
1.21k
3607
1.21k
  if (!var->is_row && 
!max_is_manifestly_unbounded(tab, var)1.20k
)
3608
118
    if (to_row(tab, var, 1) < 0)
3609
0
      return isl_stat_error;
3610
1.21k
3611
1.21k
  if (var->is_row) {
3612
123
    isl_int_sub(tab->mat->row[var->index][1],
3613
123
        tab->mat->row[var->index][1], tab->mat->row[var->index][0]);
3614
123
    if (var->is_nonneg) {
3615
123
      int sgn = restore_row(tab, var);
3616
123
      isl_assert(tab->mat->ctx, sgn >= 0,
3617
123
        return isl_stat_error);
3618
123
    }
3619
1.08k
  } else {
3620
1.08k
    int i;
3621
1.08k
3622
6.48k
    for (i = 0; i < tab->n_row; 
++i5.40k
) {
3623
5.40k
      if (isl_int_is_zero(tab->mat->row[i][off + var->index]))
3624
5.40k
        
continue4.26k
;
3625
1.14k
      isl_int_add(tab->mat->row[i][1], tab->mat->row[i][1],
3626
1.14k
          tab->mat->row[i][off + var->index]);
3627
1.14k
    }
3628
1.08k
3629
1.08k
  }
3630
1.21k
3631
1.21k
  return isl_stat_ok;
3632
1.21k
}
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
496
{
3641
496
  var->is_nonneg = 1;
3642
496
3643
496
  if (var->is_row && 
restore_row(tab, var) < -1426
)
3644
0
    return isl_stat_error;
3645
496
3646
496
  return isl_stat_ok;
3647
496
}
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
672k
{
3658
672k
  struct isl_tab_var *var;
3659
672k
3660
672k
  if (tab->n_redundant < 1)
3661
672k
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_internal,
3662
672k
      "no redundant rows", return isl_stat_error);
3663
672k
3664
672k
  var = isl_tab_var_from_row(tab, tab->n_redundant - 1);
3665
672k
  var->is_redundant = 0;
3666
672k
  tab->n_redundant--;
3667
672k
  restore_row(tab, var);
3668
672k
3669
672k
  return isl_stat_ok;
3670
672k
}
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.50M
{
3676
2.50M
  struct isl_tab_var *var = var_from_index(tab, undo->u.var_index);
3677
2.50M
  switch (undo->type) {
3678
2.50M
  case isl_tab_undo_nonneg:
3679
525k
    var->is_nonneg = 0;
3680
525k
    break;
3681
2.50M
  case isl_tab_undo_redundant:
3682
531k
    if (!var->is_row || var->index != tab->n_redundant - 1)
3683
531k
      
isl_die0
(isl_tab_get_ctx(tab), isl_error_internal,
3684
531k
        "not undoing last redundant row",
3685
531k
        return isl_stat_error);
3686
531k
    return restore_last_redundant(tab);
3687
531k
  case isl_tab_undo_freeze:
3688
256k
    var->frozen = 0;
3689
256k
    break;
3690
531k
  case isl_tab_undo_zero:
3691
84.1k
    var->is_zero = 0;
3692
84.1k
    if (!var->is_row)
3693
83.8k
      tab->n_dead--;
3694
84.1k
    break;
3695
1.10M
  case isl_tab_undo_allocate:
3696
1.10M
    if (undo->u.var_index >= 0) {
3697
5.42k
      isl_assert(tab->mat->ctx, !var->is_row,
3698
5.42k
        return isl_stat_error);
3699
5.42k
      return drop_col(tab, var->index);
3700
1.10M
    }
3701
1.10M
    if (!var->is_row) {
3702
117k
      if (!max_is_manifestly_unbounded(tab, var)) {
3703
85.5k
        if (to_row(tab, var, 1) < 0)
3704
0
          return isl_stat_error;
3705
31.6k
      } else if (!min_is_manifestly_unbounded(tab, var)) {
3706
10.6k
        if (to_row(tab, var, -1) < 0)
3707
0
          return isl_stat_error;
3708
21.0k
      } else
3709
21.0k
        if (to_row(tab, var, 0) < 0)
3710
0
          return isl_stat_error;
3711
1.10M
    }
3712
1.10M
    return drop_row(tab, var->index);
3713
1.10M
  case isl_tab_undo_relax:
3714
1.21k
    return unrelax(tab, var);
3715
1.10M
  case isl_tab_undo_unrestrict:
3716
496
    return ununrestrict(tab, var);
3717
1.10M
  default:
3718
0
    isl_die(tab->mat->ctx, isl_error_internal,
3719
2.50M
      "perform_undo_var called on invalid undo record",
3720
2.50M
      return isl_stat_error);
3721
2.50M
  }
3722
2.50M
3723
2.50M
  
return isl_stat_ok865k
;
3724
2.50M
}
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
111k
{
3735
111k
  if (!tab)
3736
0
    return isl_stat_error;
3737
111k
3738
111k
  if (tab->need_undo)
3739
111k
    
isl_die0
(isl_tab_get_ctx(tab), isl_error_invalid,
3740
111k
      "manually restoring redundant constraints "
3741
111k
      "interferes with undo history",
3742
111k
      return isl_stat_error);
3743
111k
3744
253k
  
while (111k
tab->n_redundant > 0) {
3745
141k
    if (tab->row_var[tab->n_redundant - 1] >= 0) {
3746
125k
      struct isl_tab_var *var;
3747
125k
3748
125k
      var = isl_tab_var_from_row(tab, tab->n_redundant - 1);
3749
125k
      var->is_nonneg = 0;
3750
125k
    }
3751
141k
    restore_last_redundant(tab);
3752
141k
  }
3753
111k
  return isl_stat_ok;
3754
111k
}
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.66k
{
3761
2.66k
  int off;
3762
2.66k
3763
2.66k
  off = tab->n_var - isl_basic_map_dim(tab->bmap, isl_dim_div);
3764
2.66k
  if (isl_basic_map_drop_div(tab->bmap, pos - off) < 0)
3765
0
    return isl_stat_error;
3766
2.66k
  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.66k
  }
3771
2.66k
3772
2.66k
  return isl_stat_ok;
3773
2.66k
}
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
583
{
3788
583
  int i, j;
3789
583
  int n_extra = 0;
3790
583
  int *extra = NULL;  /* current columns that contain bad stuff */
3791
583
  unsigned off = 2 + tab->M;
3792
583
3793
583
  extra = isl_alloc_array(tab->mat->ctx, int, tab->n_col);
3794
583
  if (tab->n_col && !extra)
3795
0
    goto error;
3796
6.88k
  
for (i = 0; 583
i < tab->n_col;
++i6.30k
) {
3797
51.5k
    for (j = 0; j < tab->n_col; 
++j45.2k
)
3798
50.0k
      if (tab->col_var[i] == col_var[j])
3799
4.82k
        break;
3800
6.30k
    if (j < tab->n_col)
3801
4.82k
      continue;
3802
1.48k
    extra[n_extra++] = i;
3803
1.48k
  }
3804
5.38k
  for (i = 0; i < tab->n_col && 
n_extra > 05.27k
;
++i4.80k
) {
3805
4.80k
    struct isl_tab_var *var;
3806
4.80k
    int row;
3807
4.80k
3808
36.9k
    for (j = 0; j < tab->n_col; 
++j32.1k
)
3809
35.4k
      if (col_var[i] == tab->col_var[j])
3810
3.32k
        break;
3811
4.80k
    if (j < tab->n_col)
3812
3.32k
      continue;
3813
1.48k
    var = var_from_index(tab, col_var[i]);
3814
1.48k
    row = var->index;
3815
1.88k
    for (j = 0; j < n_extra; 
++j400
)
3816
1.88k
      if (!isl_int_is_zero(tab->mat->row[row][off+extra[j]]))
3817
1.88k
        
break1.48k
;
3818
1.48k
    isl_assert(tab->mat->ctx, j < n_extra, goto error);
3819
1.48k
    if (isl_tab_pivot(tab, row, extra[j]) < 0)
3820
0
      goto error;
3821
1.48k
    extra[j] = extra[--n_extra];
3822
1.48k
  }
3823
583
3824
583
  free(extra);
3825
583
  return 0;
3826
0
error:
3827
0
  free(extra);
3828
0
  return -1;
3829
583
}
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
24.2k
{
3837
24.2k
  int i;
3838
24.2k
3839
29.5k
  for (i = tab->n_sample - 1; i >= 0 && 
tab->n_sample > n27.8k
;
--i5.29k
) {
3840
5.29k
    if (tab->sample_index[i] < n)
3841
2.03k
      continue;
3842
3.26k
3843
3.26k
    if (i != tab->n_sample - 1) {
3844
2.21k
      int t = tab->sample_index[tab->n_sample-1];
3845
2.21k
      tab->sample_index[tab->n_sample-1] = tab->sample_index[i];
3846
2.21k
      tab->sample_index[i] = t;
3847
2.21k
      isl_mat_swap_rows(tab->samples, tab->n_sample-1, i);
3848
2.21k
    }
3849
3.26k
    tab->n_sample--;
3850
3.26k
  }
3851
24.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.87M
{
3857
2.87M
  switch (undo->type) {
3858
2.87M
  case isl_tab_undo_rational:
3859
10.5k
    tab->rational = 0;
3860
10.5k
    break;
3861
2.87M
  case isl_tab_undo_empty:
3862
64.6k
    tab->empty = 0;
3863
64.6k
    break;
3864
2.87M
  case isl_tab_undo_nonneg:
3865
2.50M
  case isl_tab_undo_redundant:
3866
2.50M
  case isl_tab_undo_freeze:
3867
2.50M
  case isl_tab_undo_zero:
3868
2.50M
  case isl_tab_undo_allocate:
3869
2.50M
  case isl_tab_undo_relax:
3870
2.50M
  case isl_tab_undo_unrestrict:
3871
2.50M
    return perform_undo_var(tab, undo);
3872
2.50M
  case isl_tab_undo_bmap_eq:
3873
0
    return isl_basic_map_free_equality(tab->bmap, 1);
3874
2.50M
  case isl_tab_undo_bmap_ineq:
3875
255k
    return isl_basic_map_free_inequality(tab->bmap, 1);
3876
2.50M
  case isl_tab_undo_bmap_div:
3877
2.66k
    return drop_bmap_div(tab, undo->u.var_index);
3878
2.50M
  case isl_tab_undo_saved_basis:
3879
583
    if (restore_basis(tab, undo->u.col_var) < 0)
3880
0
      return isl_stat_error;
3881
583
    break;
3882
5.67k
  case isl_tab_undo_drop_sample:
3883
5.67k
    tab->n_outside--;
3884
5.67k
    break;
3885
24.2k
  case isl_tab_undo_saved_samples:
3886
24.2k
    drop_samples_since(tab, undo->u.n);
3887
24.2k
    break;
3888
2.38k
  case isl_tab_undo_callback:
3889
2.38k
    return undo->u.callback->run(undo->u.callback);
3890
583
  default:
3891
0
    isl_assert(tab->mat->ctx, 0, return isl_stat_error);
3892
2.87M
  }
3893
2.87M
  
return isl_stat_ok105k
;
3894
2.87M
}
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.02M
{
3901
1.02M
  struct isl_tab_undo *undo, *next;
3902
1.02M
3903
1.02M
  if (!tab)
3904
0
    return -1;
3905
1.02M
3906
1.02M
  tab->in_undo = 1;
3907
3.90M
  for (undo = tab->top; undo && undo != &tab->bottom; 
undo = next2.87M
) {
3908
3.16M
    next = undo->next;
3909
3.16M
    if (undo == snap)
3910
293k
      break;
3911
2.87M
    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.87M
    free_undo_record(undo);
3918
2.87M
  }
3919
1.02M
  tab->in_undo = 0;
3920
1.02M
  tab->top = undo;
3921
1.02M
  if (!undo)
3922
0
    return -1;
3923
1.02M
  return 0;
3924
1.02M
}
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
82.0k
{
3938
82.0k
  int pos;
3939
82.0k
  unsigned off = 2 + tab->M;
3940
82.0k
3941
82.0k
  if (tab->rational)
3942
9.20k
    return isl_ineq_separate;
3943
72.8k
3944
72.8k
  if (!isl_int_is_one(tab->mat->row[row][0]))
3945
72.8k
    
return isl_ineq_separate271
;
3946
72.5k
3947
72.5k
  pos = isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
3948
72.5k
          tab->n_col - tab->n_dead);
3949
72.5k
  if (pos == -1) {
3950
14.1k
    if (isl_int_is_negone(tab->mat->row[row][1]))
3951
14.1k
      
return isl_ineq_adj_eq11.9k
;
3952
2.24k
    else
3953
2.24k
      return isl_ineq_separate;
3954
58.4k
  }
3955
58.4k
3956
58.4k
  if (!isl_int_eq(tab->mat->row[row][1],
3957
58.4k
      tab->mat->row[row][off + tab->n_dead + pos]))
3958
58.4k
    
return isl_ineq_separate21.6k
;
3959
36.8k
3960
36.8k
  pos = isl_seq_first_non_zero(
3961
36.8k
      tab->mat->row[row] + off + tab->n_dead + pos + 1,
3962
36.8k
      tab->n_col - tab->n_dead - pos - 1);
3963
36.8k
3964
36.8k
  return pos == -1 ? 
isl_ineq_adj_ineq34.0k
:
isl_ineq_separate2.80k
;
3965
36.8k
}
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
445k
{
3977
445k
  enum isl_ineq_type type = isl_ineq_error;
3978
445k
  struct isl_tab_undo *snap = NULL;
3979
445k
  int con;
3980
445k
  int row;
3981
445k
3982
445k
  if (!tab)
3983
0
    return isl_ineq_error;
3984
445k
3985
445k
  if (isl_tab_extend_cons(tab, 1) < 0)
3986
0
    return isl_ineq_error;
3987
445k
3988
445k
  snap = isl_tab_snap(tab);
3989
445k
3990
445k
  con = isl_tab_add_row(tab, ineq);
3991
445k
  if (con < 0)
3992
0
    goto error;
3993
445k
3994
445k
  row = tab->con[con].index;
3995
445k
  if (isl_tab_row_is_redundant(tab, row))
3996
0
    type = isl_ineq_redundant;
3997
445k
  else if (isl_int_is_neg(tab->mat->row[row][1]) &&
3998
445k
     
(164k
tab->rational164k
||
3999
164k
        
isl_int_abs_ge141k
(tab->mat->row[row][1],
4000
164k
           tab->mat->row[row][0]))) {
4001
163k
    int nonneg = at_least_zero(tab, &tab->con[con]);
4002
163k
    if (nonneg < 0)
4003
0
      goto error;
4004
163k
    if (nonneg)
4005
81.3k
      type = isl_ineq_cut;
4006
82.0k
    else
4007
82.0k
      type = separation_type(tab, row);
4008
281k
  } else {
4009
281k
    int red = con_is_redundant(tab, &tab->con[con]);
4010
281k
    if (red < 0)
4011
0
      goto error;
4012
281k
    if (!red)
4013
68.9k
      type = isl_ineq_cut;
4014
212k
    else
4015
212k
      type = isl_ineq_redundant;
4016
281k
  }
4017
445k
4018
445k
  if (isl_tab_rollback(tab, snap))
4019
0
    return isl_ineq_error;
4020
445k
  return type;
4021
0
error:
4022
0
  return isl_ineq_error;
4023
445k
}
4024
4025
isl_stat isl_tab_track_bmap(struct isl_tab *tab, __isl_take isl_basic_map *bmap)
4026
194k
{
4027
194k
  bmap = isl_basic_map_cow(bmap);
4028
194k
  if (!tab || !bmap)
4029
0
    goto error;
4030
194k
4031
194k
  if (tab->empty) {
4032
4.18k
    bmap = isl_basic_map_set_to_empty(bmap);
4033
4.18k
    if (!bmap)
4034
0
      goto error;
4035
4.18k
    tab->bmap = bmap;
4036
4.18k
    return isl_stat_ok;
4037
4.18k
  }
4038
190k
4039
190k
  isl_assert(tab->mat->ctx, tab->n_eq == bmap->n_eq, goto error);
4040
190k
  isl_assert(tab->mat->ctx,
4041
190k
        tab->n_con == bmap->n_eq + bmap->n_ineq, goto error);
4042
190k
4043
190k
  tab->bmap = bmap;
4044
190k
4045
190k
  return isl_stat_ok;
4046
0
error:
4047
0
  isl_basic_map_free(bmap);
4048
0
  return isl_stat_error;
4049
190k
}
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
30.1k
{
4058
30.1k
  if (!tab)
4059
0
    return NULL;
4060
30.1k
4061
30.1k
  return bset_from_bmap(tab->bmap);
4062
30.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
}