Line data Source code
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 5455930109 : struct isl_tab *isl_tab_alloc(struct isl_ctx *ctx,
34 : unsigned n_row, unsigned n_var, unsigned M)
35 : {
36 : int i;
37 : struct isl_tab *tab;
38 5455930109 : unsigned off = 2 + M;
39 :
40 5455930109 : tab = isl_calloc_type(ctx, struct isl_tab);
41 5455930109 : if (!tab)
42 0 : return NULL;
43 5455930109 : tab->mat = isl_mat_alloc(ctx, n_row, off + n_var);
44 5455930109 : if (!tab->mat)
45 0 : goto error;
46 5455930109 : tab->var = isl_alloc_array(ctx, struct isl_tab_var, n_var);
47 5455930109 : if (n_var && !tab->var)
48 0 : goto error;
49 5455930109 : tab->con = isl_alloc_array(ctx, struct isl_tab_var, n_row);
50 5455930109 : if (n_row && !tab->con)
51 0 : goto error;
52 5455930109 : tab->col_var = isl_alloc_array(ctx, int, n_var);
53 5455930109 : if (n_var && !tab->col_var)
54 0 : goto error;
55 5455930109 : tab->row_var = isl_alloc_array(ctx, int, n_row);
56 5455930109 : if (n_row && !tab->row_var)
57 0 : goto error;
58 27999054656 : for (i = 0; i < n_var; ++i) {
59 22543124547 : tab->var[i].index = i;
60 22543124547 : tab->var[i].is_row = 0;
61 22543124547 : tab->var[i].is_nonneg = 0;
62 22543124547 : tab->var[i].is_zero = 0;
63 22543124547 : tab->var[i].is_redundant = 0;
64 22543124547 : tab->var[i].frozen = 0;
65 22543124547 : tab->var[i].negated = 0;
66 22543124547 : tab->col_var[i] = i;
67 : }
68 5455930109 : tab->n_row = 0;
69 5455930109 : tab->n_con = 0;
70 5455930109 : tab->n_eq = 0;
71 5455930109 : tab->max_con = n_row;
72 5455930109 : tab->n_col = n_var;
73 5455930109 : tab->n_var = n_var;
74 5455930109 : tab->max_var = n_var;
75 5455930109 : tab->n_param = 0;
76 5455930109 : tab->n_div = 0;
77 5455930109 : tab->n_dead = 0;
78 5455930109 : tab->n_redundant = 0;
79 5455930109 : tab->strict_redundant = 0;
80 5455930109 : tab->need_undo = 0;
81 5455930109 : tab->rational = 0;
82 5455930109 : tab->empty = 0;
83 5455930109 : tab->in_undo = 0;
84 5455930109 : tab->M = M;
85 5455930109 : tab->cone = 0;
86 5455930109 : tab->bottom.type = isl_tab_undo_bottom;
87 5455930109 : tab->bottom.next = NULL;
88 5455930109 : tab->top = &tab->bottom;
89 :
90 5455930109 : tab->n_zero = 0;
91 5455930109 : tab->n_unbounded = 0;
92 5455930109 : tab->basis = NULL;
93 :
94 5455930109 : return tab;
95 : error:
96 0 : isl_tab_free(tab);
97 0 : return NULL;
98 : }
99 :
100 24716208471 : isl_ctx *isl_tab_get_ctx(struct isl_tab *tab)
101 : {
102 24716208471 : return tab ? isl_mat_get_ctx(tab->mat) : NULL;
103 : }
104 :
105 19534137986 : int isl_tab_extend_cons(struct isl_tab *tab, unsigned n_new)
106 : {
107 : unsigned off;
108 :
109 19534137986 : if (!tab)
110 0 : return -1;
111 :
112 19534137986 : off = 2 + tab->M;
113 :
114 19534137986 : if (tab->max_con < tab->n_con + n_new) {
115 : struct isl_tab_var *con;
116 :
117 18123842 : con = isl_realloc_array(tab->mat->ctx, tab->con,
118 : struct isl_tab_var, tab->max_con + n_new);
119 18123842 : if (!con)
120 0 : return -1;
121 18123842 : tab->con = con;
122 18123842 : tab->max_con += n_new;
123 : }
124 19534137986 : if (tab->mat->n_row < tab->n_row + n_new) {
125 : int *row_var;
126 :
127 98494934 : tab->mat = isl_mat_extend(tab->mat,
128 98494934 : tab->n_row + n_new, off + tab->n_col);
129 49247467 : if (!tab->mat)
130 0 : return -1;
131 49247467 : row_var = isl_realloc_array(tab->mat->ctx, tab->row_var,
132 : int, tab->mat->n_row);
133 49247467 : if (!row_var)
134 0 : return -1;
135 49247467 : tab->row_var = row_var;
136 49247467 : if (tab->row_sign) {
137 : enum isl_tab_row_sign *s;
138 0 : s = isl_realloc_array(tab->mat->ctx, tab->row_sign,
139 : enum isl_tab_row_sign, tab->mat->n_row);
140 0 : if (!s)
141 0 : return -1;
142 0 : tab->row_sign = s;
143 : }
144 : }
145 19534137986 : return 0;
146 : }
147 :
148 : /* Make room for at least n_new extra variables.
149 : * Return -1 if anything went wrong.
150 : */
151 0 : int isl_tab_extend_vars(struct isl_tab *tab, unsigned n_new)
152 : {
153 : struct isl_tab_var *var;
154 0 : unsigned off = 2 + tab->M;
155 :
156 0 : if (tab->max_var < tab->n_var + n_new) {
157 0 : var = isl_realloc_array(tab->mat->ctx, tab->var,
158 : struct isl_tab_var, tab->n_var + n_new);
159 0 : if (!var)
160 0 : return -1;
161 0 : tab->var = var;
162 0 : tab->max_var = tab->n_var + n_new;
163 : }
164 :
165 0 : if (tab->mat->n_col < off + tab->n_col + n_new) {
166 : int *p;
167 :
168 0 : tab->mat = isl_mat_extend(tab->mat,
169 0 : tab->mat->n_row, off + tab->n_col + n_new);
170 0 : if (!tab->mat)
171 0 : return -1;
172 0 : p = isl_realloc_array(tab->mat->ctx, tab->col_var,
173 : int, tab->n_col + n_new);
174 0 : if (!p)
175 0 : return -1;
176 0 : tab->col_var = p;
177 : }
178 :
179 0 : return 0;
180 : }
181 :
182 53023335647 : static void free_undo_record(struct isl_tab_undo *undo)
183 : {
184 53023335647 : switch (undo->type) {
185 : case isl_tab_undo_saved_basis:
186 0 : free(undo->u.col_var);
187 0 : break;
188 : default:;
189 : }
190 53023335647 : free(undo);
191 53023335647 : }
192 :
193 5459292214 : static void free_undo(struct isl_tab *tab)
194 : {
195 : struct isl_tab_undo *undo, *next;
196 :
197 6110050875 : for (undo = tab->top; undo && undo != &tab->bottom; undo = next) {
198 650758661 : next = undo->next;
199 650758661 : free_undo_record(undo);
200 : }
201 5459292214 : tab->top = undo;
202 5459292214 : }
203 :
204 5459467163 : void isl_tab_free(struct isl_tab *tab)
205 : {
206 5459467163 : if (!tab)
207 174949 : return;
208 5459292214 : free_undo(tab);
209 5459292214 : isl_mat_free(tab->mat);
210 5459292214 : isl_vec_free(tab->dual);
211 5459292214 : isl_basic_map_free(tab->bmap);
212 5459292214 : free(tab->var);
213 5459292214 : free(tab->con);
214 5459292214 : free(tab->row_var);
215 5459292214 : free(tab->col_var);
216 5459292214 : free(tab->row_sign);
217 5459292214 : isl_mat_free(tab->samples);
218 5459292214 : free(tab->sample_index);
219 5459292214 : isl_mat_free(tab->basis);
220 5459292214 : free(tab);
221 : }
222 :
223 0 : struct isl_tab *isl_tab_dup(struct isl_tab *tab)
224 : {
225 : int i;
226 : struct isl_tab *dup;
227 : unsigned off;
228 :
229 0 : if (!tab)
230 0 : return NULL;
231 :
232 0 : off = 2 + tab->M;
233 0 : dup = isl_calloc_type(tab->mat->ctx, struct isl_tab);
234 0 : if (!dup)
235 0 : return NULL;
236 0 : dup->mat = isl_mat_dup(tab->mat);
237 0 : if (!dup->mat)
238 0 : goto error;
239 0 : dup->var = isl_alloc_array(tab->mat->ctx, struct isl_tab_var, tab->max_var);
240 0 : if (tab->max_var && !dup->var)
241 0 : goto error;
242 0 : for (i = 0; i < tab->n_var; ++i)
243 0 : dup->var[i] = tab->var[i];
244 0 : dup->con = isl_alloc_array(tab->mat->ctx, struct isl_tab_var, tab->max_con);
245 0 : if (tab->max_con && !dup->con)
246 0 : goto error;
247 0 : for (i = 0; i < tab->n_con; ++i)
248 0 : dup->con[i] = tab->con[i];
249 0 : dup->col_var = isl_alloc_array(tab->mat->ctx, int, tab->mat->n_col - off);
250 0 : if ((tab->mat->n_col - off) && !dup->col_var)
251 0 : goto error;
252 0 : for (i = 0; i < tab->n_col; ++i)
253 0 : dup->col_var[i] = tab->col_var[i];
254 0 : dup->row_var = isl_alloc_array(tab->mat->ctx, int, tab->mat->n_row);
255 0 : if (tab->mat->n_row && !dup->row_var)
256 0 : goto error;
257 0 : for (i = 0; i < tab->n_row; ++i)
258 0 : dup->row_var[i] = tab->row_var[i];
259 0 : if (tab->row_sign) {
260 0 : dup->row_sign = isl_alloc_array(tab->mat->ctx, enum isl_tab_row_sign,
261 : tab->mat->n_row);
262 0 : if (tab->mat->n_row && !dup->row_sign)
263 0 : goto error;
264 0 : for (i = 0; i < tab->n_row; ++i)
265 0 : dup->row_sign[i] = tab->row_sign[i];
266 : }
267 0 : 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 : 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 : }
278 0 : dup->n_row = tab->n_row;
279 0 : dup->n_con = tab->n_con;
280 0 : dup->n_eq = tab->n_eq;
281 0 : dup->max_con = tab->max_con;
282 0 : dup->n_col = tab->n_col;
283 0 : dup->n_var = tab->n_var;
284 0 : dup->max_var = tab->max_var;
285 0 : dup->n_param = tab->n_param;
286 0 : dup->n_div = tab->n_div;
287 0 : dup->n_dead = tab->n_dead;
288 0 : dup->n_redundant = tab->n_redundant;
289 0 : dup->rational = tab->rational;
290 0 : dup->empty = tab->empty;
291 0 : dup->strict_redundant = 0;
292 0 : dup->need_undo = 0;
293 0 : dup->in_undo = 0;
294 0 : dup->M = tab->M;
295 0 : tab->cone = tab->cone;
296 0 : dup->bottom.type = isl_tab_undo_bottom;
297 0 : dup->bottom.next = NULL;
298 0 : dup->top = &dup->bottom;
299 :
300 0 : dup->n_zero = tab->n_zero;
301 0 : dup->n_unbounded = tab->n_unbounded;
302 0 : dup->basis = isl_mat_dup(tab->basis);
303 :
304 0 : return dup;
305 : error:
306 0 : isl_tab_free(dup);
307 0 : return NULL;
308 : }
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 3362105 : 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 : {
328 : int i;
329 : struct isl_mat *prod;
330 : unsigned n;
331 :
332 3362105 : prod = isl_mat_alloc(mat1->ctx, mat1->n_row + mat2->n_row,
333 3362105 : off + col1 + col2);
334 3362105 : if (!prod)
335 0 : return NULL;
336 :
337 3362105 : n = 0;
338 28166891 : for (i = 0; i < r1; ++i) {
339 24804786 : isl_seq_cpy(prod->row[n + i], mat1->row[i], off + d1);
340 24804786 : isl_seq_clr(prod->row[n + i] + off + d1, d2);
341 49609572 : isl_seq_cpy(prod->row[n + i] + off + d1 + d2,
342 24804786 : mat1->row[i] + off + d1, col1 - d1);
343 24804786 : isl_seq_clr(prod->row[n + i] + off + col1 + d1, col2 - d2);
344 : }
345 :
346 3362105 : n += r1;
347 28166891 : for (i = 0; i < r2; ++i) {
348 24804786 : isl_seq_cpy(prod->row[n + i], mat2->row[i], off);
349 24804786 : isl_seq_clr(prod->row[n + i] + off, d1);
350 49609572 : isl_seq_cpy(prod->row[n + i] + off + d1,
351 49609572 : mat2->row[i] + off, d2);
352 24804786 : isl_seq_clr(prod->row[n + i] + off + d1 + d2, col1 - d1);
353 49609572 : isl_seq_cpy(prod->row[n + i] + off + col1 + d1,
354 24804786 : mat2->row[i] + off + d2, col2 - d2);
355 : }
356 :
357 3362105 : n += r2;
358 43292650 : for (i = 0; i < row1 - r1; ++i) {
359 39930545 : isl_seq_cpy(prod->row[n + i], mat1->row[r1 + i], off + d1);
360 39930545 : isl_seq_clr(prod->row[n + i] + off + d1, d2);
361 79861090 : isl_seq_cpy(prod->row[n + i] + off + d1 + d2,
362 39930545 : mat1->row[r1 + i] + off + d1, col1 - d1);
363 39930545 : isl_seq_clr(prod->row[n + i] + off + col1 + d1, col2 - d2);
364 : }
365 :
366 3362105 : n += row1 - r1;
367 43292650 : for (i = 0; i < row2 - r2; ++i) {
368 39930545 : isl_seq_cpy(prod->row[n + i], mat2->row[r2 + i], off);
369 39930545 : isl_seq_clr(prod->row[n + i] + off, d1);
370 79861090 : isl_seq_cpy(prod->row[n + i] + off + d1,
371 79861090 : mat2->row[r2 + i] + off, d2);
372 39930545 : isl_seq_clr(prod->row[n + i] + off + d1 + d2, col1 - d1);
373 79861090 : isl_seq_cpy(prod->row[n + i] + off + col1 + d1,
374 39930545 : mat2->row[r2 + i] + off + d2, col2 - d2);
375 : }
376 :
377 3362105 : return prod;
378 : }
379 :
380 : /* Update the row or column index of a variable that corresponds
381 : * to a variable in the first input tableau.
382 : */
383 81892547 : static void update_index1(struct isl_tab_var *var,
384 : unsigned r1, unsigned r2, unsigned d1, unsigned d2)
385 : {
386 81892547 : if (var->index == -1)
387 451254 : return;
388 81441293 : if (var->is_row && var->index >= r1)
389 39930545 : var->index += r2;
390 81441293 : if (!var->is_row && var->index >= d1)
391 15495213 : var->index += d2;
392 : }
393 :
394 : /* Update the row or column index of a variable that corresponds
395 : * to a variable in the second input tableau.
396 : */
397 81892547 : static void update_index2(struct isl_tab_var *var,
398 : unsigned row1, unsigned col1,
399 : unsigned r1, unsigned r2, unsigned d1, unsigned d2)
400 : {
401 81892547 : if (var->index == -1)
402 451254 : return;
403 81441293 : if (var->is_row) {
404 64735331 : if (var->index < r2)
405 24804786 : var->index += r1;
406 : else
407 39930545 : var->index += row1;
408 : } else {
409 16705962 : if (var->index < d2)
410 1210749 : var->index += d1;
411 : else
412 15495213 : var->index += col1;
413 : }
414 : }
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 3362105 : struct isl_tab *isl_tab_product(struct isl_tab *tab1, struct isl_tab *tab2)
435 : {
436 : int i;
437 : struct isl_tab *prod;
438 : unsigned off;
439 : unsigned r1, r2, d1, d2;
440 :
441 3362105 : if (!tab1 || !tab2)
442 0 : return NULL;
443 :
444 3362105 : isl_assert(tab1->mat->ctx, tab1->M == tab2->M, return NULL);
445 3362105 : isl_assert(tab1->mat->ctx, tab1->rational == tab2->rational, return NULL);
446 3362105 : isl_assert(tab1->mat->ctx, tab1->cone == tab2->cone, return NULL);
447 3362105 : isl_assert(tab1->mat->ctx, !tab1->row_sign, return NULL);
448 3362105 : isl_assert(tab1->mat->ctx, !tab2->row_sign, return NULL);
449 3362105 : isl_assert(tab1->mat->ctx, tab1->n_param == 0, return NULL);
450 3362105 : isl_assert(tab1->mat->ctx, tab2->n_param == 0, return NULL);
451 3362105 : isl_assert(tab1->mat->ctx, tab1->n_div == 0, return NULL);
452 3362105 : isl_assert(tab1->mat->ctx, tab2->n_div == 0, return NULL);
453 :
454 3362105 : off = 2 + tab1->M;
455 3362105 : r1 = tab1->n_redundant;
456 3362105 : r2 = tab2->n_redundant;
457 3362105 : d1 = tab1->n_dead;
458 3362105 : d2 = tab2->n_dead;
459 3362105 : prod = isl_calloc_type(tab1->mat->ctx, struct isl_tab);
460 3362105 : if (!prod)
461 0 : return NULL;
462 3362105 : prod->mat = tab_mat_product(tab1->mat, tab2->mat,
463 : tab1->n_row, tab2->n_row,
464 : tab1->n_col, tab2->n_col, off, r1, r2, d1, d2);
465 3362105 : if (!prod->mat)
466 0 : goto error;
467 3362105 : prod->var = isl_alloc_array(tab1->mat->ctx, struct isl_tab_var,
468 : tab1->max_var + tab2->max_var);
469 3362105 : if ((tab1->max_var + tab2->max_var) && !prod->var)
470 0 : goto error;
471 20519321 : for (i = 0; i < tab1->n_var; ++i) {
472 17157216 : prod->var[i] = tab1->var[i];
473 17157216 : update_index1(&prod->var[i], r1, r2, d1, d2);
474 : }
475 20519321 : for (i = 0; i < tab2->n_var; ++i) {
476 17157216 : prod->var[tab1->n_var + i] = tab2->var[i];
477 17157216 : update_index2(&prod->var[tab1->n_var + i],
478 : tab1->n_row, tab1->n_col,
479 : r1, r2, d1, d2);
480 : }
481 3362105 : prod->con = isl_alloc_array(tab1->mat->ctx, struct isl_tab_var,
482 : tab1->max_con + tab2->max_con);
483 3362105 : if ((tab1->max_con + tab2->max_con) && !prod->con)
484 0 : goto error;
485 68097436 : for (i = 0; i < tab1->n_con; ++i) {
486 64735331 : prod->con[i] = tab1->con[i];
487 64735331 : update_index1(&prod->con[i], r1, r2, d1, d2);
488 : }
489 68097436 : for (i = 0; i < tab2->n_con; ++i) {
490 64735331 : prod->con[tab1->n_con + i] = tab2->con[i];
491 64735331 : update_index2(&prod->con[tab1->n_con + i],
492 : tab1->n_row, tab1->n_col,
493 : r1, r2, d1, d2);
494 : }
495 3362105 : prod->col_var = isl_alloc_array(tab1->mat->ctx, int,
496 : tab1->n_col + tab2->n_col);
497 3362105 : if ((tab1->n_col + tab2->n_col) && !prod->col_var)
498 0 : goto error;
499 20068067 : for (i = 0; i < tab1->n_col; ++i) {
500 16705962 : int pos = i < d1 ? i : i + d2;
501 16705962 : prod->col_var[pos] = tab1->col_var[i];
502 : }
503 20068067 : for (i = 0; i < tab2->n_col; ++i) {
504 16705962 : int pos = i < d2 ? d1 + i : tab1->n_col + i;
505 16705962 : int t = tab2->col_var[i];
506 16705962 : if (t >= 0)
507 5809 : t += tab1->n_var;
508 : else
509 16700153 : t -= tab1->n_con;
510 16705962 : prod->col_var[pos] = t;
511 : }
512 3362105 : prod->row_var = isl_alloc_array(tab1->mat->ctx, int,
513 : tab1->mat->n_row + tab2->mat->n_row);
514 3362105 : if ((tab1->mat->n_row + tab2->mat->n_row) && !prod->row_var)
515 0 : goto error;
516 68097436 : for (i = 0; i < tab1->n_row; ++i) {
517 64735331 : int pos = i < r1 ? i : i + r2;
518 64735331 : prod->row_var[pos] = tab1->row_var[i];
519 : }
520 68097436 : for (i = 0; i < tab2->n_row; ++i) {
521 64735331 : int pos = i < r2 ? r1 + i : tab1->n_row + i;
522 64735331 : int t = tab2->row_var[i];
523 64735331 : if (t >= 0)
524 17151407 : t += tab1->n_var;
525 : else
526 47583924 : t -= tab1->n_con;
527 64735331 : prod->row_var[pos] = t;
528 : }
529 3362105 : prod->samples = NULL;
530 3362105 : prod->sample_index = NULL;
531 3362105 : prod->n_row = tab1->n_row + tab2->n_row;
532 3362105 : prod->n_con = tab1->n_con + tab2->n_con;
533 3362105 : prod->n_eq = 0;
534 3362105 : prod->max_con = tab1->max_con + tab2->max_con;
535 3362105 : prod->n_col = tab1->n_col + tab2->n_col;
536 3362105 : prod->n_var = tab1->n_var + tab2->n_var;
537 3362105 : prod->max_var = tab1->max_var + tab2->max_var;
538 3362105 : prod->n_param = 0;
539 3362105 : prod->n_div = 0;
540 3362105 : prod->n_dead = tab1->n_dead + tab2->n_dead;
541 3362105 : prod->n_redundant = tab1->n_redundant + tab2->n_redundant;
542 3362105 : prod->rational = tab1->rational;
543 3362105 : prod->empty = tab1->empty || tab2->empty;
544 3362105 : prod->strict_redundant = tab1->strict_redundant || tab2->strict_redundant;
545 3362105 : prod->need_undo = 0;
546 3362105 : prod->in_undo = 0;
547 3362105 : prod->M = tab1->M;
548 3362105 : prod->cone = tab1->cone;
549 3362105 : prod->bottom.type = isl_tab_undo_bottom;
550 3362105 : prod->bottom.next = NULL;
551 3362105 : prod->top = &prod->bottom;
552 :
553 3362105 : prod->n_zero = 0;
554 3362105 : prod->n_unbounded = 0;
555 3362105 : prod->basis = NULL;
556 :
557 3362105 : return prod;
558 : error:
559 0 : isl_tab_free(prod);
560 0 : return NULL;
561 : }
562 :
563 >51995*10^7 : static struct isl_tab_var *var_from_index(struct isl_tab *tab, int i)
564 : {
565 >51995*10^7 : if (i >= 0)
566 >12922*10^7 : return &tab->var[i];
567 : else
568 >39072*10^7 : return &tab->con[~i];
569 : }
570 :
571 >31758*10^7 : struct isl_tab_var *isl_tab_var_from_row(struct isl_tab *tab, int i)
572 : {
573 >31758*10^7 : return var_from_index(tab, tab->row_var[i]);
574 : }
575 :
576 >15448*10^7 : static struct isl_tab_var *var_from_col(struct isl_tab *tab, int i)
577 : {
578 >15448*10^7 : return var_from_index(tab, tab->col_var[i]);
579 : }
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 1294877277 : static int max_is_manifestly_unbounded(struct isl_tab *tab,
586 : struct isl_tab_var *var)
587 : {
588 : int i;
589 1294877277 : unsigned off = 2 + tab->M;
590 :
591 1294877277 : if (var->is_row)
592 616703802 : return 0;
593 3849437880 : for (i = tab->n_redundant; i < tab->n_row; ++i) {
594 3637926229 : if (!isl_int_is_neg(tab->mat->row[i][off + var->index]))
595 2722679996 : continue;
596 915246233 : if (isl_tab_var_from_row(tab, i)->is_nonneg)
597 466661824 : return 0;
598 : }
599 211511651 : return 1;
600 : }
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 12817740870 : static int min_is_manifestly_unbounded(struct isl_tab *tab,
607 : struct isl_tab_var *var)
608 : {
609 : int i;
610 12817740870 : unsigned off = 2 + tab->M;
611 :
612 12817740870 : if (var->is_row)
613 11854101727 : return 0;
614 4573174600 : for (i = tab->n_redundant; i < tab->n_row; ++i) {
615 4429753381 : if (!isl_int_is_pos(tab->mat->row[i][off + var->index]))
616 3091662991 : continue;
617 1338090390 : if (isl_tab_var_from_row(tab, i)->is_nonneg)
618 820217924 : return 0;
619 : }
620 143421219 : return 1;
621 : }
622 :
623 21161341164 : static int row_cmp(struct isl_tab *tab, int r1, int r2, int c, isl_int *t)
624 : {
625 21161341164 : unsigned off = 2 + tab->M;
626 :
627 21161341164 : if (tab->M) {
628 : 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 : }
635 21161341164 : isl_int_mul(*t, tab->mat->row[r1][1], tab->mat->row[r2][off + c]);
636 21161341164 : isl_int_submul(*t, tab->mat->row[r2][1], tab->mat->row[r1][off + c]);
637 21161341164 : return isl_int_sgn(*t);
638 : }
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 5707718385 : static int pivot_row(struct isl_tab *tab,
661 : struct isl_tab_var *var, int sgn, int c)
662 : {
663 : int j, r, tsgn;
664 : isl_int t;
665 5707718385 : unsigned off = 2 + tab->M;
666 :
667 5707718385 : isl_int_init(t);
668 5707718385 : r = -1;
669 99518206739 : for (j = tab->n_redundant; j < tab->n_row; ++j) {
670 93810488354 : if (var && j == var->index)
671 4892239600 : continue;
672 88918248754 : if (!isl_tab_var_from_row(tab, j)->is_nonneg)
673 30121924285 : continue;
674 58796324469 : if (sgn * isl_int_sgn(tab->mat->row[j][off + c]) >= 0)
675 32360293081 : continue;
676 26436031388 : if (r < 0) {
677 5274690224 : r = j;
678 5274690224 : continue;
679 : }
680 21161341164 : tsgn = sgn * row_cmp(tab, r, j, c, &t);
681 24558580068 : if (tsgn < 0 || (tsgn == 0 &&
682 3397238904 : tab->row_var[j] < tab->row_var[r]))
683 6486609515 : r = j;
684 : }
685 5707718385 : isl_int_clear(t);
686 5707718385 : return r;
687 : }
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 27492068353 : 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 : {
709 : int j, r, c;
710 : isl_int *tr;
711 :
712 27492068353 : *row = *col = -1;
713 :
714 27492068353 : isl_assert(tab->mat->ctx, var->is_row, return);
715 27492068353 : tr = tab->mat->row[var->index] + 2 + tab->M;
716 :
717 27492068353 : c = -1;
718 59239487936 : for (j = tab->n_dead; j < tab->n_col; ++j) {
719 31747419583 : if (isl_int_is_zero(tr[j]))
720 10216870646 : continue;
721 35463640499 : if (isl_int_sgn(tr[j]) != sgn &&
722 13933091562 : var_from_col(tab, j)->is_nonneg)
723 11845678892 : continue;
724 9684870045 : if (c < 0 || tab->col_var[j] < tab->col_var[c])
725 6490498734 : c = j;
726 : }
727 27492068353 : if (c < 0)
728 22438075420 : return;
729 :
730 5053992933 : sgn *= isl_int_sgn(tr[c]);
731 5053992933 : r = pivot_row(tab, skip_var, sgn, c);
732 5053992933 : *row = r < 0 ? var->index : r;
733 5053992933 : *col = c;
734 : }
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 >10885*10^7 : int isl_tab_row_is_redundant(struct isl_tab *tab, int row)
743 : {
744 : int i;
745 >10885*10^7 : unsigned off = 2 + tab->M;
746 :
747 >10885*10^7 : if (tab->row_var[row] < 0 && !isl_tab_var_from_row(tab, row)->is_nonneg)
748 20541672816 : return 0;
749 :
750 88309122804 : if (isl_int_is_neg(tab->mat->row[row][1]))
751 9873288850 : return 0;
752 78435833954 : if (tab->strict_redundant && isl_int_is_zero(tab->mat->row[row][1]))
753 0 : return 0;
754 78435833954 : if (tab->M && isl_int_is_neg(tab->mat->row[row][2]))
755 0 : return 0;
756 :
757 >24678*10^7 : for (i = tab->n_dead; i < tab->n_col; ++i) {
758 >23988*10^7 : if (isl_int_is_zero(tab->mat->row[row][off + i]))
759 >12850*10^7 : continue;
760 >11137*10^7 : if (tab->col_var[i] >= 0)
761 11925854604 : return 0;
762 99448020500 : if (isl_int_is_neg(tab->mat->row[row][off + i]))
763 44564035637 : return 0;
764 54883984863 : if (!var_from_col(tab, i)->is_nonneg)
765 15044768923 : return 0;
766 : }
767 6901174790 : return 1;
768 : }
769 :
770 7469002579 : static void swap_rows(struct isl_tab *tab, int row1, int row2)
771 : {
772 : int t;
773 : enum isl_tab_row_sign s;
774 :
775 7469002579 : t = tab->row_var[row1];
776 7469002579 : tab->row_var[row1] = tab->row_var[row2];
777 7469002579 : tab->row_var[row2] = t;
778 7469002579 : isl_tab_var_from_row(tab, row1)->index = row1;
779 7469002579 : isl_tab_var_from_row(tab, row2)->index = row2;
780 7469002579 : tab->mat = isl_mat_swap_rows(tab->mat, row1, row2);
781 :
782 7469002579 : if (!tab->row_sign)
783 7469002579 : return;
784 0 : s = tab->row_sign[row1];
785 0 : tab->row_sign[row1] = tab->row_sign[row2];
786 0 : tab->row_sign[row2] = s;
787 : }
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 84629203562 : static isl_stat push_union(struct isl_tab *tab,
800 : enum isl_tab_undo_type type, union isl_tab_undo_val u)
801 : {
802 : struct isl_tab_undo *undo;
803 :
804 84629203562 : if (!tab)
805 0 : return isl_stat_error;
806 84629203562 : if (!tab->need_undo)
807 31605867915 : return isl_stat_ok;
808 :
809 53023335647 : undo = isl_alloc_type(tab->mat->ctx, struct isl_tab_undo);
810 53023335647 : if (!undo)
811 0 : goto error;
812 53023335647 : undo->type = type;
813 53023335647 : undo->u = u;
814 53023335647 : undo->next = tab->top;
815 53023335647 : tab->top = undo;
816 :
817 53023335647 : return isl_stat_ok;
818 : error:
819 0 : free_undo(tab);
820 0 : tab->top = NULL;
821 0 : return isl_stat_error;
822 : }
823 :
824 80049571913 : isl_stat isl_tab_push_var(struct isl_tab *tab,
825 : enum isl_tab_undo_type type, struct isl_tab_var *var)
826 : {
827 : union isl_tab_undo_val u;
828 80049571913 : if (var->is_row)
829 76295722248 : u.var_index = tab->row_var[var->index];
830 : else
831 3753849665 : u.var_index = tab->col_var[var->index];
832 80049571913 : return push_union(tab, type, u);
833 : }
834 :
835 4579631649 : isl_stat isl_tab_push(struct isl_tab *tab, enum isl_tab_undo_type type)
836 : {
837 4579631649 : union isl_tab_undo_val u = { 0 };
838 4579631649 : return push_union(tab, type, u);
839 : }
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 0 : isl_stat isl_tab_push_basis(struct isl_tab *tab)
845 : {
846 : int i;
847 : union isl_tab_undo_val u;
848 :
849 0 : u.col_var = isl_alloc_array(tab->mat->ctx, int, tab->n_col);
850 0 : if (tab->n_col && !u.col_var)
851 0 : return isl_stat_error;
852 0 : for (i = 0; i < tab->n_col; ++i)
853 0 : u.col_var[i] = tab->col_var[i];
854 0 : return push_union(tab, isl_tab_undo_saved_basis, u);
855 : }
856 :
857 0 : isl_stat isl_tab_push_callback(struct isl_tab *tab,
858 : struct isl_tab_callback *callback)
859 : {
860 : union isl_tab_undo_val u;
861 0 : u.callback = callback;
862 0 : return push_union(tab, isl_tab_undo_callback, u);
863 : }
864 :
865 0 : struct isl_tab *isl_tab_init_samples(struct isl_tab *tab)
866 : {
867 0 : if (!tab)
868 0 : return NULL;
869 :
870 0 : tab->n_sample = 0;
871 0 : tab->n_outside = 0;
872 0 : tab->samples = isl_mat_alloc(tab->mat->ctx, 1, 1 + tab->n_var);
873 0 : if (!tab->samples)
874 0 : goto error;
875 0 : tab->sample_index = isl_alloc_array(tab->mat->ctx, int, 1);
876 0 : if (!tab->sample_index)
877 0 : goto error;
878 0 : return tab;
879 : error:
880 0 : isl_tab_free(tab);
881 0 : return NULL;
882 : }
883 :
884 0 : int isl_tab_add_sample(struct isl_tab *tab, __isl_take isl_vec *sample)
885 : {
886 0 : if (!tab || !sample)
887 : goto error;
888 :
889 0 : if (tab->n_sample + 1 > tab->samples->n_row) {
890 0 : int *t = isl_realloc_array(tab->mat->ctx,
891 : tab->sample_index, int, tab->n_sample + 1);
892 0 : if (!t)
893 0 : goto error;
894 0 : tab->sample_index = t;
895 : }
896 :
897 0 : tab->samples = isl_mat_extend(tab->samples,
898 0 : tab->n_sample + 1, tab->samples->n_col);
899 0 : if (!tab->samples)
900 0 : goto error;
901 :
902 0 : isl_seq_cpy(tab->samples->row[tab->n_sample], sample->el, sample->size);
903 0 : isl_vec_free(sample);
904 0 : tab->sample_index[tab->n_sample] = tab->n_sample;
905 0 : tab->n_sample++;
906 :
907 0 : return 0;
908 : error:
909 0 : isl_vec_free(sample);
910 0 : return -1;
911 : }
912 :
913 0 : struct isl_tab *isl_tab_drop_sample(struct isl_tab *tab, int s)
914 : {
915 0 : if (s != tab->n_outside) {
916 0 : int t = tab->sample_index[tab->n_outside];
917 0 : tab->sample_index[tab->n_outside] = tab->sample_index[s];
918 0 : tab->sample_index[s] = t;
919 0 : isl_mat_swap_rows(tab->samples, tab->n_outside, s);
920 : }
921 0 : tab->n_outside++;
922 0 : if (isl_tab_push(tab, isl_tab_undo_drop_sample) < 0) {
923 0 : isl_tab_free(tab);
924 0 : return NULL;
925 : }
926 :
927 0 : return tab;
928 : }
929 :
930 : /* Record the current number of samples so that we can remove newer
931 : * samples during a rollback.
932 : */
933 0 : isl_stat isl_tab_save_samples(struct isl_tab *tab)
934 : {
935 : union isl_tab_undo_val u;
936 :
937 0 : if (!tab)
938 0 : return isl_stat_error;
939 :
940 0 : u.n = tab->n_sample;
941 0 : return push_union(tab, isl_tab_undo_saved_samples, u);
942 : }
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 7292744120 : int isl_tab_mark_redundant(struct isl_tab *tab, int row)
957 : {
958 7292744120 : struct isl_tab_var *var = isl_tab_var_from_row(tab, row);
959 7292744120 : var->is_redundant = 1;
960 7292744120 : isl_assert(tab->mat->ctx, row >= tab->n_redundant, return -1);
961 7292744120 : if (tab->preserve || tab->need_undo || tab->row_var[row] >= 0) {
962 5825065685 : if (tab->row_var[row] >= 0 && !var->is_nonneg) {
963 3579057047 : var->is_nonneg = 1;
964 3579057047 : if (isl_tab_push_var(tab, isl_tab_undo_nonneg, var) < 0)
965 0 : return -1;
966 : }
967 5825065685 : if (row != tab->n_redundant)
968 5184096694 : swap_rows(tab, row, tab->n_redundant);
969 5825065685 : tab->n_redundant++;
970 5825065685 : return isl_tab_push_var(tab, isl_tab_undo_redundant, var);
971 : } else {
972 1467678435 : if (row != tab->n_row - 1)
973 544074291 : swap_rows(tab, row, tab->n_row - 1);
974 1467678435 : isl_tab_var_from_row(tab, tab->n_row - 1)->index = -1;
975 1467678435 : tab->n_row--;
976 1467678435 : return 1;
977 : }
978 : }
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 11179 : int isl_tab_mark_rational(struct isl_tab *tab)
986 : {
987 11179 : if (!tab)
988 0 : return -1;
989 11179 : if (!tab->rational && tab->need_undo)
990 10643 : if (isl_tab_push(tab, isl_tab_undo_rational) < 0)
991 0 : return -1;
992 11179 : tab->rational = 1;
993 11179 : return 0;
994 : }
995 :
996 805460820 : isl_stat isl_tab_mark_empty(struct isl_tab *tab)
997 : {
998 805460820 : if (!tab)
999 0 : return isl_stat_error;
1000 805460820 : if (!tab->empty && tab->need_undo)
1001 749697170 : if (isl_tab_push(tab, isl_tab_undo_empty) < 0)
1002 0 : return isl_stat_error;
1003 805460820 : tab->empty = 1;
1004 805460820 : return isl_stat_ok;
1005 : }
1006 :
1007 10910083355 : int isl_tab_freeze_constraint(struct isl_tab *tab, int con)
1008 : {
1009 : struct isl_tab_var *var;
1010 :
1011 10910083355 : if (!tab)
1012 0 : return -1;
1013 :
1014 10910083355 : var = &tab->con[con];
1015 10910083355 : if (var->frozen)
1016 0 : return 0;
1017 10910083355 : if (var->index < 0)
1018 311580467 : return 0;
1019 10598502888 : var->frozen = 1;
1020 :
1021 10598502888 : if (tab->need_undo)
1022 10598502888 : return isl_tab_push_var(tab, isl_tab_undo_freeze, var);
1023 :
1024 0 : return 0;
1025 : }
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 24716205376 : static void update_row_sign(struct isl_tab *tab, int row, int col, int row_sgn)
1045 : {
1046 : int i;
1047 24716205376 : struct isl_mat *mat = tab->mat;
1048 24716205376 : unsigned off = 2 + tab->M;
1049 :
1050 24716205376 : if (!tab->row_sign)
1051 24716205376 : return;
1052 :
1053 0 : if (tab->row_sign[row] == 0)
1054 0 : return;
1055 0 : isl_assert(mat->ctx, row_sgn > 0, return);
1056 0 : isl_assert(mat->ctx, tab->row_sign[row] == isl_tab_row_neg, return);
1057 0 : tab->row_sign[row] = isl_tab_row_pos;
1058 0 : for (i = 0; i < tab->n_row; ++i) {
1059 : int s;
1060 0 : if (i == row)
1061 0 : continue;
1062 0 : s = isl_int_sgn(mat->row[i][off + col]);
1063 0 : if (!s)
1064 0 : continue;
1065 0 : if (!tab->row_sign[i])
1066 0 : continue;
1067 0 : if (s < 0 && tab->row_sign[i] == isl_tab_row_neg)
1068 0 : continue;
1069 0 : if (s > 0 && tab->row_sign[i] == isl_tab_row_pos)
1070 0 : continue;
1071 0 : tab->row_sign[i] = isl_tab_row_unknown;
1072 : }
1073 : }
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 24716205376 : int isl_tab_pivot(struct isl_tab *tab, int row, int col)
1126 : {
1127 : int i, j;
1128 : int sgn;
1129 : int t;
1130 : isl_ctx *ctx;
1131 24716205376 : struct isl_mat *mat = tab->mat;
1132 : struct isl_tab_var *var;
1133 24716205376 : unsigned off = 2 + tab->M;
1134 :
1135 24716205376 : ctx = isl_tab_get_ctx(tab);
1136 24716205376 : if (isl_ctx_next_operation(ctx) < 0)
1137 0 : return -1;
1138 :
1139 24716205376 : isl_int_swap(mat->row[row][0], mat->row[row][off + col]);
1140 24716205376 : sgn = isl_int_sgn(mat->row[row][0]);
1141 24716205376 : if (sgn < 0) {
1142 5990149498 : isl_int_neg(mat->row[row][0], mat->row[row][0]);
1143 5990149498 : isl_int_neg(mat->row[row][off + col], mat->row[row][off + col]);
1144 : } else
1145 >14923*10^7 : for (j = 0; j < off - 1 + tab->n_col; ++j) {
1146 >13050*10^7 : if (j == off - 1 + col)
1147 18726055878 : continue;
1148 >11178*10^7 : isl_int_neg(mat->row[row][1 + j], mat->row[row][1 + j]);
1149 : }
1150 24716205376 : if (!isl_int_is_one(mat->row[row][0]))
1151 8098515831 : isl_seq_normalize(mat->ctx, mat->row[row], off + tab->n_col);
1152 >23370*10^7 : for (i = 0; i < tab->n_row; ++i) {
1153 >20899*10^7 : if (i == row)
1154 24716205376 : continue;
1155 >18427*10^7 : if (isl_int_is_zero(mat->row[i][off + col]))
1156 >10698*10^7 : continue;
1157 77294032788 : isl_int_mul(mat->row[i][0], mat->row[i][0], mat->row[row][0]);
1158 >62802*10^7 : for (j = 0; j < off - 1 + tab->n_col; ++j) {
1159 >55073*10^7 : if (j == off - 1 + col)
1160 77294032788 : continue;
1161 >47343*10^7 : isl_int_mul(mat->row[i][1 + j],
1162 : mat->row[i][1 + j], mat->row[row][0]);
1163 >47343*10^7 : isl_int_addmul(mat->row[i][1 + j],
1164 : mat->row[i][off + col], mat->row[row][1 + j]);
1165 : }
1166 77294032788 : isl_int_mul(mat->row[i][off + col],
1167 : mat->row[i][off + col], mat->row[row][off + col]);
1168 77294032788 : if (!isl_int_is_one(mat->row[i][0]))
1169 74322745097 : isl_seq_normalize(mat->ctx, mat->row[i], off + tab->n_col);
1170 : }
1171 24716205376 : t = tab->row_var[row];
1172 24716205376 : tab->row_var[row] = tab->col_var[col];
1173 24716205376 : tab->col_var[col] = t;
1174 24716205376 : var = isl_tab_var_from_row(tab, row);
1175 24716205376 : var->is_row = 1;
1176 24716205376 : var->index = row;
1177 24716205376 : var = var_from_col(tab, col);
1178 24716205376 : var->is_row = 0;
1179 24716205376 : var->index = col;
1180 24716205376 : update_row_sign(tab, row, col, sgn);
1181 24716205376 : if (tab->in_undo)
1182 381339560 : return 0;
1183 >20834*10^7 : for (i = tab->n_redundant; i < tab->n_row; ++i) {
1184 >18400*10^7 : if (isl_int_is_zero(mat->row[i][off + col]))
1185 97410802631 : continue;
1186 >16017*10^7 : if (!isl_tab_var_from_row(tab, i)->frozen &&
1187 73575099776 : isl_tab_row_is_redundant(tab, i)) {
1188 4812270692 : int redo = isl_tab_mark_redundant(tab, i);
1189 4812270692 : if (redo < 0)
1190 0 : return -1;
1191 4812270692 : if (redo)
1192 230379432 : --i;
1193 : }
1194 : }
1195 24334865816 : return 0;
1196 : }
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 1171567299 : static int to_row(struct isl_tab *tab, struct isl_tab_var *var, int sign)
1206 : {
1207 : int r;
1208 1171567299 : unsigned off = 2 + tab->M;
1209 :
1210 1171567299 : if (var->is_row)
1211 619370503 : return 0;
1212 :
1213 552196796 : if (sign == 0) {
1214 157771337 : for (r = tab->n_redundant; r < tab->n_row; ++r)
1215 157771337 : if (!isl_int_is_zero(tab->mat->row[r][off+var->index]))
1216 : break;
1217 47261337 : isl_assert(tab->mat->ctx, r < tab->n_row, return -1);
1218 : } else {
1219 504935459 : r = pivot_row(tab, NULL, sign, var->index);
1220 504935459 : isl_assert(tab->mat->ctx, r >= 0, return -1);
1221 : }
1222 :
1223 552196796 : return isl_tab_pivot(tab, r, var->index);
1224 : }
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 0 : static void check_table(struct isl_tab *tab)
1232 : {
1233 : int i;
1234 :
1235 0 : if (tab->empty)
1236 0 : return;
1237 0 : for (i = tab->n_redundant; i < tab->n_row; ++i) {
1238 : 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 : !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 : }
1248 0 : isl_assert(tab->mat->ctx, !isl_int_is_neg(tab->mat->row[i][1]),
1249 : abort());
1250 : }
1251 : }
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 670044224 : static int sign_of_max(struct isl_tab *tab, struct isl_tab_var *var)
1265 : {
1266 : int row, col;
1267 :
1268 670044224 : if (max_is_manifestly_unbounded(tab, var))
1269 34854614 : return 1;
1270 635189610 : if (to_row(tab, var, 1) < 0)
1271 0 : return -2;
1272 1630498137 : while (!isl_int_is_pos(tab->mat->row[var->index][1])) {
1273 802973078 : find_pivot(tab, var, var, 1, &row, &col);
1274 802973078 : if (row == -1)
1275 387244848 : return isl_int_sgn(tab->mat->row[var->index][1]);
1276 415728230 : if (isl_tab_pivot(tab, row, col) < 0)
1277 0 : return -2;
1278 415728230 : if (!var->is_row) /* manifestly unbounded */
1279 55609313 : return 1;
1280 : }
1281 192335449 : return 1;
1282 : }
1283 :
1284 0 : int isl_tab_sign_of_max(struct isl_tab *tab, int con)
1285 : {
1286 : struct isl_tab_var *var;
1287 :
1288 0 : if (!tab)
1289 0 : return -2;
1290 :
1291 0 : var = &tab->con[con];
1292 0 : isl_assert(tab->mat->ctx, !var->is_redundant, return -2);
1293 0 : isl_assert(tab->mat->ctx, !var->is_zero, return -2);
1294 :
1295 0 : return sign_of_max(tab, var);
1296 : }
1297 :
1298 14699681607 : static int row_is_neg(struct isl_tab *tab, int row)
1299 : {
1300 14699681607 : if (!tab->M)
1301 14699681607 : 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 : }
1308 :
1309 12998964181 : static int row_sgn(struct isl_tab *tab, int row)
1310 : {
1311 12998964181 : if (!tab->M)
1312 12998964181 : 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 : else
1316 0 : return isl_int_sgn(tab->mat->row[row][1]);
1317 : }
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 13286314264 : static int restore_row(struct isl_tab *tab, struct isl_tab_var *var)
1325 : {
1326 : int row, col;
1327 :
1328 27985995871 : while (row_is_neg(tab, var->index)) {
1329 2506500974 : find_pivot(tab, var, var, 1, &row, &col);
1330 2506500974 : if (row == -1)
1331 805783548 : break;
1332 1700717426 : if (isl_tab_pivot(tab, row, col) < 0)
1333 0 : return -2;
1334 1700717426 : if (!var->is_row) /* manifestly unbounded */
1335 287350083 : return 1;
1336 : }
1337 12998964181 : return row_sgn(tab, var->index);
1338 : }
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 7208278280 : static int at_least_zero(struct isl_tab *tab, struct isl_tab_var *var)
1346 : {
1347 : int row, col;
1348 :
1349 14431794877 : while (isl_int_is_neg(tab->mat->row[var->index][1])) {
1350 7216359581 : find_pivot(tab, var, var, 1, &row, &col);
1351 7216359581 : if (row == -1)
1352 7196337328 : break;
1353 20022253 : if (row == var->index) /* manifestly unbounded */
1354 4783936 : return 1;
1355 15238317 : if (isl_tab_pivot(tab, row, col) < 0)
1356 0 : return -1;
1357 : }
1358 7203494344 : return !isl_int_is_neg(tab->mat->row[var->index][1]);
1359 : }
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 11417362933 : static int sign_of_min(struct isl_tab *tab, struct isl_tab_var *var)
1379 : {
1380 : int row, col;
1381 11417362933 : struct isl_tab_var *pivot_var = NULL;
1382 :
1383 11417362933 : if (min_is_manifestly_unbounded(tab, var))
1384 0 : return -1;
1385 11417362933 : if (!var->is_row) {
1386 925086 : col = var->index;
1387 925086 : row = pivot_row(tab, NULL, -1, col);
1388 925086 : pivot_var = var_from_col(tab, col);
1389 925086 : if (isl_tab_pivot(tab, row, col) < 0)
1390 0 : return -2;
1391 925086 : if (var->is_redundant)
1392 147649 : return 0;
1393 777437 : if (isl_int_is_neg(tab->mat->row[var->index][1])) {
1394 461312 : if (var->is_nonneg) {
1395 922624 : if (!pivot_var->is_redundant &&
1396 461312 : pivot_var->index == row) {
1397 860244 : if (isl_tab_pivot(tab, row, col) < 0)
1398 0 : return -2;
1399 : } else
1400 31190 : if (restore_row(tab, var) < -1)
1401 0 : return -2;
1402 : }
1403 461312 : return -1;
1404 : }
1405 : }
1406 11416753972 : if (var->is_redundant)
1407 0 : return 0;
1408 22858463627 : while (!isl_int_is_neg(tab->mat->row[var->index][1])) {
1409 11427722368 : find_pivot(tab, var, var, -1, &row, &col);
1410 11427722368 : if (row == var->index)
1411 26222183 : return -1;
1412 11401500185 : if (row == -1)
1413 11376053425 : return isl_int_sgn(tab->mat->row[var->index][1]);
1414 25446760 : pivot_var = var_from_col(tab, col);
1415 25446760 : if (isl_tab_pivot(tab, row, col) < 0)
1416 0 : return -2;
1417 25446760 : if (var->is_redundant)
1418 491077 : return 0;
1419 : }
1420 13987287 : if (pivot_var && var->is_nonneg) {
1421 : /* pivot back to non-negative value */
1422 523082 : if (!pivot_var->is_redundant && pivot_var->index == row) {
1423 994522 : if (isl_tab_pivot(tab, row, col) < 0)
1424 0 : return -2;
1425 : } else
1426 25821 : if (restore_row(tab, var) < -1)
1427 0 : return -2;
1428 : }
1429 13987287 : return -1;
1430 : }
1431 :
1432 836082677 : static int row_at_most_neg_one(struct isl_tab *tab, int row)
1433 : {
1434 836082677 : 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 : }
1440 1092270702 : return isl_int_is_neg(tab->mat->row[row][1]) &&
1441 256188025 : isl_int_abs_ge(tab->mat->row[row][1],
1442 : tab->mat->row[row][0]);
1443 : }
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 582862086 : int isl_tab_min_at_most_neg_one(struct isl_tab *tab, struct isl_tab_var *var)
1455 : {
1456 : int row, col;
1457 : struct isl_tab_var *pivot_var;
1458 :
1459 582862086 : if (min_is_manifestly_unbounded(tab, var))
1460 0 : return 1;
1461 582862086 : if (!var->is_row) {
1462 147864907 : col = var->index;
1463 147864907 : row = pivot_row(tab, NULL, -1, col);
1464 147864907 : pivot_var = var_from_col(tab, col);
1465 147864907 : if (isl_tab_pivot(tab, row, col) < 0)
1466 0 : return -1;
1467 147864907 : if (var->is_redundant)
1468 24661944 : return 0;
1469 123202963 : if (row_at_most_neg_one(tab, var->index)) {
1470 80502838 : if (var->is_nonneg) {
1471 161005676 : if (!pivot_var->is_redundant &&
1472 80502838 : pivot_var->index == row) {
1473 159918060 : if (isl_tab_pivot(tab, row, col) < 0)
1474 0 : return -1;
1475 : } else
1476 543808 : if (restore_row(tab, var) < -1)
1477 0 : return -1;
1478 : }
1479 80502838 : return 1;
1480 : }
1481 : }
1482 477697304 : if (var->is_redundant)
1483 2 : return 0;
1484 : do {
1485 1047951548 : find_pivot(tab, var, var, -1, &row, &col);
1486 1047951548 : if (row == var->index) {
1487 56687786 : if (var->is_nonneg && restore_row(tab, var) < -1)
1488 0 : return -1;
1489 56687786 : return 1;
1490 : }
1491 991263762 : if (row == -1)
1492 109144629 : return 0;
1493 882119133 : pivot_var = var_from_col(tab, col);
1494 882119133 : if (isl_tab_pivot(tab, row, col) < 0)
1495 0 : return -1;
1496 882119133 : if (var->is_redundant)
1497 169239419 : return 0;
1498 712879714 : } while (!row_at_most_neg_one(tab, var->index));
1499 142625468 : if (var->is_nonneg) {
1500 : /* pivot back to non-negative value */
1501 138736701 : if (!pivot_var->is_redundant && pivot_var->index == row)
1502 137853481 : if (isl_tab_pivot(tab, row, col) < 0)
1503 0 : return -1;
1504 138736701 : if (restore_row(tab, var) < -1)
1505 0 : return -1;
1506 : }
1507 142625468 : return 1;
1508 : }
1509 :
1510 : /* Return 1 if "var" can attain values >= 1.
1511 : * Return 0 otherwise.
1512 : */
1513 245410099 : static int at_least_one(struct isl_tab *tab, struct isl_tab_var *var)
1514 : {
1515 : int row, col;
1516 : isl_int *r;
1517 :
1518 245410099 : if (max_is_manifestly_unbounded(tab, var))
1519 90463927 : return 1;
1520 154946172 : if (to_row(tab, var, 1) < 0)
1521 0 : return -1;
1522 154946172 : r = tab->mat->row[var->index];
1523 323523064 : while (isl_int_lt(r[1], r[0])) {
1524 17368142 : find_pivot(tab, var, var, 1, &row, &col);
1525 17368142 : if (row == -1)
1526 2878036 : return isl_int_ge(r[1], r[0]);
1527 14490106 : if (row == var->index) /* manifestly unbounded */
1528 859386 : return 1;
1529 13630720 : if (isl_tab_pivot(tab, row, col) < 0)
1530 0 : return -1;
1531 : }
1532 151208750 : return 1;
1533 : }
1534 :
1535 9032103822 : static void swap_cols(struct isl_tab *tab, int col1, int col2)
1536 : {
1537 : int t;
1538 9032103822 : unsigned off = 2 + tab->M;
1539 9032103822 : t = tab->col_var[col1];
1540 9032103822 : tab->col_var[col1] = tab->col_var[col2];
1541 9032103822 : tab->col_var[col2] = t;
1542 9032103822 : var_from_col(tab, col1)->index = col1;
1543 9032103822 : var_from_col(tab, col2)->index = col2;
1544 9032103822 : tab->mat = isl_mat_swap_cols(tab->mat, off + col1, off + col2);
1545 9032103822 : }
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 19128236403 : int isl_tab_kill_col(struct isl_tab *tab, int col)
1560 : {
1561 19128236403 : var_from_col(tab, col)->is_zero = 1;
1562 19128236403 : if (tab->need_undo) {
1563 103087761 : if (isl_tab_push_var(tab, isl_tab_undo_zero,
1564 : var_from_col(tab, col)) < 0)
1565 0 : return -1;
1566 103087761 : if (col != tab->n_dead)
1567 65094659 : swap_cols(tab, col, tab->n_dead);
1568 103087761 : tab->n_dead++;
1569 103087761 : return 0;
1570 : } else {
1571 19025148642 : if (col != tab->n_col - 1)
1572 8967009163 : swap_cols(tab, col, tab->n_col - 1);
1573 19025148642 : var_from_col(tab, tab->n_col - 1)->index = -1;
1574 19025148642 : tab->n_col--;
1575 19025148642 : return 1;
1576 : }
1577 : }
1578 :
1579 1840570125 : static int row_is_manifestly_non_integral(struct isl_tab *tab, int row)
1580 : {
1581 1840570125 : unsigned off = 2 + tab->M;
1582 :
1583 1840570125 : if (tab->M && !isl_int_eq(tab->mat->row[row][2],
1584 : tab->mat->row[row][0]))
1585 0 : return 0;
1586 1840570125 : if (isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
1587 1840570125 : tab->n_col - tab->n_dead) != -1)
1588 40848126 : return 0;
1589 :
1590 1799721999 : return !isl_int_is_divisible_by(tab->mat->row[row][1],
1591 : tab->mat->row[row][0]);
1592 : }
1593 :
1594 : /* For integer tableaus, check if any of the coordinates are stuck
1595 : * at a non-integral value.
1596 : */
1597 387194228 : static int tab_is_manifestly_empty(struct isl_tab *tab)
1598 : {
1599 : int i;
1600 :
1601 387194228 : if (tab->empty)
1602 2572 : return 1;
1603 387191656 : if (tab->rational)
1604 1034763 : return 0;
1605 :
1606 2246235839 : for (i = 0; i < tab->n_var; ++i) {
1607 1860158993 : if (!tab->var[i].is_row)
1608 19588868 : continue;
1609 1840570125 : if (row_is_manifestly_non_integral(tab, tab->var[i].index))
1610 80047 : return 1;
1611 : }
1612 :
1613 386076846 : return 0;
1614 : }
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 387194228 : static isl_stat close_row(struct isl_tab *tab, struct isl_tab_var *var,
1636 : int temp_var)
1637 : {
1638 : int j;
1639 387194228 : struct isl_mat *mat = tab->mat;
1640 387194228 : unsigned off = 2 + tab->M;
1641 :
1642 387194228 : if (!var->is_nonneg)
1643 0 : isl_die(isl_tab_get_ctx(tab), isl_error_internal,
1644 : "expecting non-negative variable",
1645 : return isl_stat_error);
1646 387194228 : var->is_zero = 1;
1647 387194228 : if (!temp_var && tab->need_undo)
1648 0 : if (isl_tab_push_var(tab, isl_tab_undo_zero, var) < 0)
1649 0 : return isl_stat_error;
1650 638202601 : for (j = tab->n_dead; j < tab->n_col; ++j) {
1651 : int recheck;
1652 251008373 : if (isl_int_is_zero(mat->row[var->index][off + j]))
1653 48132206 : continue;
1654 202876167 : if (isl_int_is_pos(mat->row[var->index][off + j]))
1655 0 : isl_die(isl_tab_get_ctx(tab), isl_error_internal,
1656 : "row cannot have positive coefficients",
1657 : return isl_stat_error);
1658 202876167 : recheck = isl_tab_kill_col(tab, j);
1659 202876167 : if (recheck < 0)
1660 0 : return isl_stat_error;
1661 202876167 : if (recheck)
1662 202490385 : --j;
1663 : }
1664 387194228 : if (!temp_var && isl_tab_mark_redundant(tab, var->index) < 0)
1665 0 : return isl_stat_error;
1666 387194228 : if (tab_is_manifestly_empty(tab) && isl_tab_mark_empty(tab) < 0)
1667 0 : return isl_stat_error;
1668 387194228 : return isl_stat_ok;
1669 : }
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 50070934420 : int isl_tab_allocate_con(struct isl_tab *tab)
1678 : {
1679 : int r;
1680 :
1681 50070934420 : isl_assert(tab->mat->ctx, tab->n_row < tab->mat->n_row, return -1);
1682 50070934420 : isl_assert(tab->mat->ctx, tab->n_con < tab->max_con, return -1);
1683 :
1684 50070934420 : r = tab->n_con;
1685 50070934420 : tab->con[r].index = tab->n_row;
1686 50070934420 : tab->con[r].is_row = 1;
1687 50070934420 : tab->con[r].is_nonneg = 0;
1688 50070934420 : tab->con[r].is_zero = 0;
1689 50070934420 : tab->con[r].is_redundant = 0;
1690 50070934420 : tab->con[r].frozen = 0;
1691 50070934420 : tab->con[r].negated = 0;
1692 50070934420 : tab->row_var[tab->n_row] = ~r;
1693 :
1694 50070934420 : tab->n_row++;
1695 50070934420 : tab->n_con++;
1696 50070934420 : if (isl_tab_push_var(tab, isl_tab_undo_allocate, &tab->con[r]) < 0)
1697 0 : return -1;
1698 :
1699 50070934420 : return r;
1700 : }
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 0 : static int var_insert_entry(struct isl_tab *tab, int first)
1708 : {
1709 : int i;
1710 :
1711 0 : if (tab->n_var >= tab->max_var)
1712 0 : isl_die(isl_tab_get_ctx(tab), isl_error_internal,
1713 : "not enough room for new variable", return -1);
1714 0 : if (first > tab->n_var)
1715 0 : isl_die(isl_tab_get_ctx(tab), isl_error_internal,
1716 : "invalid initial position", return -1);
1717 :
1718 0 : for (i = tab->n_var - 1; i >= first; --i) {
1719 0 : tab->var[i + 1] = tab->var[i];
1720 0 : if (tab->var[i + 1].is_row)
1721 0 : tab->row_var[tab->var[i + 1].index]++;
1722 : else
1723 0 : tab->col_var[tab->var[i + 1].index]++;
1724 : }
1725 :
1726 0 : tab->n_var++;
1727 :
1728 0 : return 0;
1729 : }
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 0 : static int var_drop_entry(struct isl_tab *tab, int first)
1737 : {
1738 : int i;
1739 :
1740 0 : if (first >= tab->n_var)
1741 0 : isl_die(isl_tab_get_ctx(tab), isl_error_internal,
1742 : "invalid initial position", return -1);
1743 :
1744 0 : tab->n_var--;
1745 :
1746 0 : for (i = first; i < tab->n_var; ++i) {
1747 0 : tab->var[i] = tab->var[i + 1];
1748 0 : if (tab->var[i + 1].is_row)
1749 0 : tab->row_var[tab->var[i].index]--;
1750 : else
1751 0 : tab->col_var[tab->var[i].index]--;
1752 : }
1753 :
1754 0 : return 0;
1755 : }
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 0 : int isl_tab_insert_var(struct isl_tab *tab, int r)
1762 : {
1763 : int i;
1764 0 : unsigned off = 2 + tab->M;
1765 :
1766 0 : isl_assert(tab->mat->ctx, tab->n_col < tab->mat->n_col, return -1);
1767 :
1768 0 : if (var_insert_entry(tab, r) < 0)
1769 0 : return -1;
1770 :
1771 0 : tab->var[r].index = tab->n_col;
1772 0 : tab->var[r].is_row = 0;
1773 0 : tab->var[r].is_nonneg = 0;
1774 0 : tab->var[r].is_zero = 0;
1775 0 : tab->var[r].is_redundant = 0;
1776 0 : tab->var[r].frozen = 0;
1777 0 : tab->var[r].negated = 0;
1778 0 : tab->col_var[tab->n_col] = r;
1779 :
1780 0 : for (i = 0; i < tab->n_row; ++i)
1781 0 : isl_int_set_si(tab->mat->row[i][off + tab->n_col], 0);
1782 :
1783 0 : tab->n_col++;
1784 0 : if (isl_tab_push_var(tab, isl_tab_undo_allocate, &tab->var[r]) < 0)
1785 0 : return -1;
1786 :
1787 0 : return r;
1788 : }
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 0 : int isl_tab_allocate_var(struct isl_tab *tab)
1794 : {
1795 0 : if (!tab)
1796 0 : return -1;
1797 :
1798 0 : return isl_tab_insert_var(tab, tab->n_var);
1799 : }
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 50070910073 : int isl_tab_add_row(struct isl_tab *tab, isl_int *line)
1824 : {
1825 : int i;
1826 : int r;
1827 : isl_int *row;
1828 : isl_int a, b;
1829 50070910073 : unsigned off = 2 + tab->M;
1830 :
1831 50070910073 : r = isl_tab_allocate_con(tab);
1832 50070910073 : if (r < 0)
1833 0 : return -1;
1834 :
1835 50070910073 : isl_int_init(a);
1836 50070910073 : isl_int_init(b);
1837 50070910073 : row = tab->mat->row[tab->con[r].index];
1838 50070910073 : isl_int_set_si(row[0], 1);
1839 50070910073 : isl_int_set(row[1], line[0]);
1840 50070910073 : isl_seq_clr(row + 2, tab->M + tab->n_col);
1841 >40675*10^7 : for (i = 0; i < tab->n_var; ++i) {
1842 >35668*10^7 : if (tab->var[i].is_zero)
1843 0 : continue;
1844 >35668*10^7 : if (tab->var[i].is_row) {
1845 >20606*10^7 : isl_int_lcm(a,
1846 : row[0], tab->mat->row[tab->var[i].index][0]);
1847 >20606*10^7 : isl_int_swap(a, row[0]);
1848 >20606*10^7 : isl_int_divexact(a, row[0], a);
1849 >20606*10^7 : isl_int_divexact(b,
1850 : row[0], tab->mat->row[tab->var[i].index][0]);
1851 >20606*10^7 : isl_int_mul(b, b, line[1 + i]);
1852 >41213*10^7 : isl_seq_combine(row + 1, a, row + 1,
1853 >20606*10^7 : b, tab->mat->row[tab->var[i].index] + 1,
1854 >20606*10^7 : 1 + tab->M + tab->n_col);
1855 : } else
1856 >15061*10^7 : isl_int_addmul(row[off + tab->var[i].index],
1857 : line[1 + i], row[0]);
1858 >35668*10^7 : if (tab->M && i >= tab->n_param && i < tab->n_var - tab->n_div)
1859 0 : isl_int_submul(row[2], line[1 + i], row[0]);
1860 : }
1861 50070910073 : isl_seq_normalize(tab->mat->ctx, row, off + tab->n_col);
1862 50070910073 : isl_int_clear(a);
1863 50070910073 : isl_int_clear(b);
1864 :
1865 50070910073 : if (tab->row_sign)
1866 0 : tab->row_sign[tab->con[r].index] = isl_tab_row_unknown;
1867 :
1868 50070910073 : return r;
1869 : }
1870 :
1871 25065482622 : static isl_stat drop_row(struct isl_tab *tab, int row)
1872 : {
1873 25065482622 : isl_assert(tab->mat->ctx, ~tab->row_var[row] == tab->n_con - 1,
1874 : return isl_stat_error);
1875 25065482622 : if (row != tab->n_row - 1)
1876 1740831594 : swap_rows(tab, row, tab->n_row - 1);
1877 25065482622 : tab->n_row--;
1878 25065482622 : tab->n_con--;
1879 25065482622 : return isl_stat_ok;
1880 : }
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 0 : static isl_stat drop_col(struct isl_tab *tab, int col)
1889 : {
1890 : int var;
1891 :
1892 0 : var = tab->col_var[col];
1893 0 : if (col != tab->n_col - 1)
1894 0 : swap_cols(tab, col, tab->n_col - 1);
1895 0 : tab->n_col--;
1896 0 : if (var_drop_entry(tab, var) < 0)
1897 0 : return isl_stat_error;
1898 0 : return isl_stat_ok;
1899 : }
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 9871784566 : isl_stat isl_tab_add_ineq(struct isl_tab *tab, isl_int *ineq)
1908 : {
1909 : int r;
1910 : int sgn;
1911 : isl_int cst;
1912 :
1913 9871784566 : if (!tab)
1914 0 : return isl_stat_error;
1915 9871784566 : if (tab->bmap) {
1916 3829875142 : struct isl_basic_map *bmap = tab->bmap;
1917 :
1918 3829875142 : isl_assert(tab->mat->ctx, tab->n_eq == bmap->n_eq,
1919 : return isl_stat_error);
1920 3829875142 : isl_assert(tab->mat->ctx,
1921 : tab->n_con == bmap->n_eq + bmap->n_ineq,
1922 : return isl_stat_error);
1923 3829875142 : tab->bmap = isl_basic_map_add_ineq(tab->bmap, ineq);
1924 3829875142 : if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1925 0 : return isl_stat_error;
1926 3829875142 : if (!tab->bmap)
1927 0 : return isl_stat_error;
1928 : }
1929 9871784566 : if (tab->cone) {
1930 0 : isl_int_init(cst);
1931 0 : isl_int_set_si(cst, 0);
1932 0 : isl_int_swap(ineq[0], cst);
1933 : }
1934 9871784566 : r = isl_tab_add_row(tab, ineq);
1935 9871784566 : if (tab->cone) {
1936 0 : isl_int_swap(ineq[0], cst);
1937 0 : isl_int_clear(cst);
1938 : }
1939 9871784566 : if (r < 0)
1940 0 : return isl_stat_error;
1941 9871784566 : tab->con[r].is_nonneg = 1;
1942 9871784566 : if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1943 0 : return isl_stat_error;
1944 9871784566 : if (isl_tab_row_is_redundant(tab, tab->con[r].index)) {
1945 2088904098 : if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0)
1946 0 : return isl_stat_error;
1947 2088904098 : return isl_stat_ok;
1948 : }
1949 :
1950 7782880468 : sgn = restore_row(tab, &tab->con[r]);
1951 7782880468 : if (sgn < -1)
1952 0 : return isl_stat_error;
1953 7782880468 : if (sgn < 0)
1954 805333643 : return isl_tab_mark_empty(tab);
1955 6977546825 : if (tab->con[r].is_row && isl_tab_row_is_redundant(tab, tab->con[r].index))
1956 0 : if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0)
1957 0 : return isl_stat_error;
1958 6977546825 : return isl_stat_ok;
1959 : }
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 102701979 : static int to_col(struct isl_tab *tab, struct isl_tab_var *var)
1966 : {
1967 : int i;
1968 : int row, col;
1969 102701979 : unsigned off = 2 + tab->M;
1970 :
1971 102701979 : if (!var->is_row)
1972 0 : return 0;
1973 :
1974 273449658 : while (isl_int_is_pos(tab->mat->row[var->index][1])) {
1975 160851411 : find_pivot(tab, var, NULL, -1, &row, &col);
1976 160851411 : isl_assert(tab->mat->ctx, row != -1, return -1);
1977 160851411 : if (isl_tab_pivot(tab, row, col) < 0)
1978 0 : return -1;
1979 160851411 : if (!var->is_row)
1980 92805711 : return 0;
1981 : }
1982 :
1983 16766505 : for (i = tab->n_dead; i < tab->n_col; ++i)
1984 16766505 : if (!isl_int_is_zero(tab->mat->row[var->index][off + i]))
1985 : break;
1986 :
1987 9896268 : isl_assert(tab->mat->ctx, i < tab->n_col, return -1);
1988 9896268 : if (isl_tab_pivot(tab, var->index, i) < 0)
1989 0 : return -1;
1990 :
1991 9896268 : return 0;
1992 : }
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 18822658257 : static struct isl_tab *add_eq(struct isl_tab *tab, isl_int *eq)
2003 : {
2004 : int i;
2005 : int r;
2006 :
2007 18822658257 : if (!tab)
2008 0 : return NULL;
2009 18822658257 : r = isl_tab_add_row(tab, eq);
2010 18822658257 : if (r < 0)
2011 0 : goto error;
2012 :
2013 18822658257 : r = tab->con[r].index;
2014 18822658257 : i = isl_seq_first_non_zero(tab->mat->row[r] + 2 + tab->M + tab->n_dead,
2015 18822658257 : tab->n_col - tab->n_dead);
2016 18822658257 : isl_assert(tab->mat->ctx, i >= 0, goto error);
2017 18822658257 : i += tab->n_dead;
2018 18822658257 : if (isl_tab_pivot(tab, r, i) < 0)
2019 0 : goto error;
2020 18822658257 : if (isl_tab_kill_col(tab, i) < 0)
2021 0 : goto error;
2022 18822658257 : tab->n_eq++;
2023 :
2024 18822658257 : return tab;
2025 : error:
2026 0 : isl_tab_free(tab);
2027 0 : return NULL;
2028 : }
2029 :
2030 : /* Does the sample value of row "row" of "tab" involve the big parameter,
2031 : * if any?
2032 : */
2033 96541618 : static int row_is_big(struct isl_tab *tab, int row)
2034 : {
2035 96541618 : return tab->M && !isl_int_is_zero(tab->mat->row[row][2]);
2036 : }
2037 :
2038 104041679 : static int row_is_manifestly_zero(struct isl_tab *tab, int row)
2039 : {
2040 104041679 : unsigned off = 2 + tab->M;
2041 :
2042 104041679 : if (!isl_int_is_zero(tab->mat->row[row][1]))
2043 92811773 : return 0;
2044 11229906 : if (row_is_big(tab, row))
2045 0 : return 0;
2046 22459812 : return isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
2047 22459812 : tab->n_col - tab->n_dead) == -1;
2048 : }
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 104008423 : int isl_tab_add_valid_eq(struct isl_tab *tab, isl_int *eq)
2056 : {
2057 : struct isl_tab_var *var;
2058 : int r;
2059 :
2060 104008423 : if (!tab)
2061 0 : return -1;
2062 104008423 : r = isl_tab_add_row(tab, eq);
2063 104008423 : if (r < 0)
2064 0 : return -1;
2065 :
2066 104008423 : var = &tab->con[r];
2067 104008423 : r = var->index;
2068 104008423 : if (row_is_manifestly_zero(tab, r)) {
2069 1335891 : var->is_zero = 1;
2070 1335891 : if (isl_tab_mark_redundant(tab, r) < 0)
2071 0 : return -1;
2072 1335891 : return 0;
2073 : }
2074 :
2075 102672532 : if (isl_int_is_neg(tab->mat->row[r][1])) {
2076 29157967 : isl_seq_neg(tab->mat->row[r] + 1, tab->mat->row[r] + 1,
2077 29157967 : 1 + tab->n_col);
2078 29157967 : var->negated = 1;
2079 : }
2080 102672532 : var->is_nonneg = 1;
2081 102672532 : if (to_col(tab, var) < 0)
2082 0 : return -1;
2083 102672532 : var->is_nonneg = 0;
2084 102672532 : if (isl_tab_kill_col(tab, var->index) < 0)
2085 0 : return -1;
2086 :
2087 102672532 : return 0;
2088 : }
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 24347 : static int add_zero_row(struct isl_tab *tab)
2097 : {
2098 : int r;
2099 : isl_int *row;
2100 :
2101 24347 : r = isl_tab_allocate_con(tab);
2102 24347 : if (r < 0)
2103 0 : return -1;
2104 :
2105 24347 : row = tab->mat->row[tab->con[r].index];
2106 24347 : isl_seq_clr(row + 1, 1 + tab->M + tab->n_col);
2107 24347 : isl_int_set_si(row[0], 1);
2108 :
2109 24347 : return r;
2110 : }
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 33256 : int isl_tab_add_eq(struct isl_tab *tab, isl_int *eq)
2120 : {
2121 33256 : struct isl_tab_undo *snap = NULL;
2122 : struct isl_tab_var *var;
2123 : int r;
2124 : int row;
2125 : int sgn;
2126 : isl_int cst;
2127 :
2128 33256 : if (!tab)
2129 0 : return -1;
2130 33256 : isl_assert(tab->mat->ctx, !tab->M, return -1);
2131 :
2132 33256 : if (tab->need_undo)
2133 33256 : snap = isl_tab_snap(tab);
2134 :
2135 33256 : if (tab->cone) {
2136 0 : isl_int_init(cst);
2137 0 : isl_int_set_si(cst, 0);
2138 0 : isl_int_swap(eq[0], cst);
2139 : }
2140 33256 : r = isl_tab_add_row(tab, eq);
2141 33256 : if (tab->cone) {
2142 0 : isl_int_swap(eq[0], cst);
2143 0 : isl_int_clear(cst);
2144 : }
2145 33256 : if (r < 0)
2146 0 : return -1;
2147 :
2148 33256 : var = &tab->con[r];
2149 33256 : row = var->index;
2150 33256 : if (row_is_manifestly_zero(tab, row)) {
2151 3809 : if (snap)
2152 3809 : return isl_tab_rollback(tab, snap);
2153 0 : return drop_row(tab, row);
2154 : }
2155 :
2156 29447 : if (tab->bmap) {
2157 24347 : tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);
2158 24347 : if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
2159 0 : return -1;
2160 24347 : isl_seq_neg(eq, eq, 1 + tab->n_var);
2161 24347 : tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);
2162 24347 : isl_seq_neg(eq, eq, 1 + tab->n_var);
2163 24347 : if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
2164 0 : return -1;
2165 24347 : if (!tab->bmap)
2166 0 : return -1;
2167 24347 : if (add_zero_row(tab) < 0)
2168 0 : return -1;
2169 : }
2170 :
2171 29447 : sgn = isl_int_sgn(tab->mat->row[row][1]);
2172 :
2173 29447 : if (sgn > 0) {
2174 1154 : isl_seq_neg(tab->mat->row[row] + 1, tab->mat->row[row] + 1,
2175 1154 : 1 + tab->n_col);
2176 1154 : var->negated = 1;
2177 1154 : sgn = -1;
2178 : }
2179 :
2180 29447 : if (sgn < 0) {
2181 16959 : sgn = sign_of_max(tab, var);
2182 16959 : if (sgn < -1)
2183 0 : return -1;
2184 16959 : if (sgn < 0) {
2185 0 : if (isl_tab_mark_empty(tab) < 0)
2186 0 : return -1;
2187 0 : return 0;
2188 : }
2189 : }
2190 :
2191 29447 : var->is_nonneg = 1;
2192 29447 : if (to_col(tab, var) < 0)
2193 0 : return -1;
2194 29447 : var->is_nonneg = 0;
2195 29447 : if (isl_tab_kill_col(tab, var->index) < 0)
2196 0 : return -1;
2197 :
2198 29447 : return 0;
2199 : }
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 0 : static struct isl_vec *ineq_for_div(struct isl_basic_map *bmap, unsigned div)
2212 : {
2213 : unsigned total;
2214 : unsigned div_pos;
2215 : struct isl_vec *ineq;
2216 :
2217 0 : if (!bmap)
2218 0 : return NULL;
2219 :
2220 0 : total = isl_basic_map_total_dim(bmap);
2221 0 : div_pos = 1 + total - bmap->n_div + div;
2222 :
2223 0 : ineq = isl_vec_alloc(bmap->ctx, 1 + total);
2224 0 : if (!ineq)
2225 0 : return NULL;
2226 :
2227 0 : isl_seq_cpy(ineq->el, bmap->div[div] + 1, 1 + total);
2228 0 : isl_int_neg(ineq->el[div_pos], bmap->div[div][0]);
2229 0 : return ineq;
2230 : }
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 0 : static isl_stat add_div_constraints(struct isl_tab *tab, unsigned div,
2248 : isl_stat (*add_ineq)(void *user, isl_int *), void *user)
2249 : {
2250 : unsigned total;
2251 : unsigned div_pos;
2252 : struct isl_vec *ineq;
2253 :
2254 0 : total = isl_basic_map_total_dim(tab->bmap);
2255 0 : div_pos = 1 + total - tab->bmap->n_div + div;
2256 :
2257 0 : ineq = ineq_for_div(tab->bmap, div);
2258 0 : if (!ineq)
2259 0 : goto error;
2260 :
2261 0 : if (add_ineq) {
2262 0 : if (add_ineq(user, ineq->el) < 0)
2263 0 : goto error;
2264 : } else {
2265 0 : if (isl_tab_add_ineq(tab, ineq->el) < 0)
2266 0 : goto error;
2267 : }
2268 :
2269 0 : isl_seq_neg(ineq->el, tab->bmap->div[div] + 1, 1 + total);
2270 0 : isl_int_set(ineq->el[div_pos], tab->bmap->div[div][0]);
2271 0 : isl_int_add(ineq->el[0], ineq->el[0], ineq->el[div_pos]);
2272 0 : isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
2273 :
2274 0 : if (add_ineq) {
2275 0 : if (add_ineq(user, ineq->el) < 0)
2276 0 : goto error;
2277 : } else {
2278 0 : if (isl_tab_add_ineq(tab, ineq->el) < 0)
2279 0 : goto error;
2280 : }
2281 :
2282 0 : isl_vec_free(ineq);
2283 :
2284 0 : return isl_stat_ok;
2285 : error:
2286 0 : isl_vec_free(ineq);
2287 0 : return isl_stat_error;
2288 : }
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 0 : static int div_is_nonneg(struct isl_tab *tab, __isl_keep isl_vec *div)
2297 : {
2298 : int i;
2299 :
2300 0 : if (tab->M)
2301 0 : return 1;
2302 :
2303 0 : if (isl_int_is_neg(div->el[1]))
2304 0 : return 0;
2305 :
2306 0 : for (i = 0; i < tab->n_var; ++i) {
2307 0 : if (isl_int_is_neg(div->el[2 + i]))
2308 0 : return 0;
2309 0 : if (isl_int_is_zero(div->el[2 + i]))
2310 0 : continue;
2311 0 : if (!tab->var[i].is_nonneg)
2312 0 : return 0;
2313 : }
2314 :
2315 0 : return 1;
2316 : }
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 0 : 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 : {
2332 : int r;
2333 : int nonneg;
2334 : int n_div, o_div;
2335 :
2336 0 : if (!tab || !div)
2337 0 : return -1;
2338 :
2339 0 : if (div->size != 1 + 1 + tab->n_var)
2340 0 : isl_die(isl_tab_get_ctx(tab), isl_error_invalid,
2341 : "unexpected size", return -1);
2342 :
2343 0 : isl_assert(tab->mat->ctx, tab->bmap, return -1);
2344 0 : n_div = isl_basic_map_dim(tab->bmap, isl_dim_div);
2345 0 : o_div = tab->n_var - n_div;
2346 0 : if (pos < o_div || pos > tab->n_var)
2347 0 : isl_die(isl_tab_get_ctx(tab), isl_error_invalid,
2348 : "invalid position", return -1);
2349 :
2350 0 : nonneg = div_is_nonneg(tab, div);
2351 :
2352 0 : if (isl_tab_extend_cons(tab, 3) < 0)
2353 0 : return -1;
2354 0 : if (isl_tab_extend_vars(tab, 1) < 0)
2355 0 : return -1;
2356 0 : r = isl_tab_insert_var(tab, pos);
2357 0 : if (r < 0)
2358 0 : return -1;
2359 :
2360 0 : if (nonneg)
2361 0 : tab->var[r].is_nonneg = 1;
2362 :
2363 0 : tab->bmap = isl_basic_map_insert_div(tab->bmap, pos - o_div, div);
2364 0 : if (!tab->bmap)
2365 0 : return -1;
2366 0 : if (isl_tab_push_var(tab, isl_tab_undo_bmap_div, &tab->var[r]) < 0)
2367 0 : return -1;
2368 :
2369 0 : if (add_div_constraints(tab, pos - o_div, add_ineq, user) < 0)
2370 0 : return -1;
2371 :
2372 0 : return r;
2373 : }
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 0 : int isl_tab_add_div(struct isl_tab *tab, __isl_keep isl_vec *div)
2379 : {
2380 0 : if (!tab)
2381 0 : return -1;
2382 0 : return isl_tab_insert_div(tab, tab->n_var, div, NULL, NULL);
2383 : }
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 5455020296 : __isl_give struct isl_tab *isl_tab_from_basic_map(
2391 : __isl_keep isl_basic_map *bmap, int track)
2392 : {
2393 : int i;
2394 : struct isl_tab *tab;
2395 :
2396 5455020296 : if (!bmap)
2397 0 : return NULL;
2398 10910040592 : tab = isl_tab_alloc(bmap->ctx,
2399 5455020296 : isl_basic_map_total_dim(bmap) + bmap->n_ineq + 1,
2400 : isl_basic_map_total_dim(bmap), 0);
2401 5455020296 : if (!tab)
2402 0 : return NULL;
2403 5455020296 : tab->preserve = track;
2404 5455020296 : tab->rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
2405 5455020296 : if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY)) {
2406 0 : if (isl_tab_mark_empty(tab) < 0)
2407 0 : goto error;
2408 0 : goto done;
2409 : }
2410 24273214798 : for (i = 0; i < bmap->n_eq; ++i) {
2411 18818194502 : tab = add_eq(tab, bmap->eq[i]);
2412 18818194502 : if (!tab)
2413 0 : return tab;
2414 : }
2415 11464381692 : for (i = 0; i < bmap->n_ineq; ++i) {
2416 6038590142 : if (isl_tab_add_ineq(tab, bmap->ineq[i]) < 0)
2417 0 : goto error;
2418 6038590142 : if (tab->empty)
2419 29228746 : goto done;
2420 : }
2421 : done:
2422 5455020296 : if (track && isl_tab_track_bmap(tab, isl_basic_map_copy(bmap)) < 0)
2423 0 : goto error;
2424 5455020296 : return tab;
2425 : error:
2426 0 : isl_tab_free(tab);
2427 0 : return NULL;
2428 : }
2429 :
2430 66139487 : __isl_give struct isl_tab *isl_tab_from_basic_set(
2431 : __isl_keep isl_basic_set *bset, int track)
2432 : {
2433 66139487 : return isl_tab_from_basic_map(bset, track);
2434 : }
2435 :
2436 : /* Construct a tableau corresponding to the recession cone of "bset".
2437 : */
2438 909813 : struct isl_tab *isl_tab_from_recession_cone(__isl_keep isl_basic_set *bset,
2439 : int parametric)
2440 : {
2441 : isl_int cst;
2442 : int i;
2443 : struct isl_tab *tab;
2444 909813 : unsigned offset = 0;
2445 :
2446 909813 : if (!bset)
2447 0 : return NULL;
2448 909813 : if (parametric)
2449 909813 : offset = isl_basic_set_dim(bset, isl_dim_param);
2450 909813 : tab = isl_tab_alloc(bset->ctx, bset->n_eq + bset->n_ineq,
2451 909813 : isl_basic_set_total_dim(bset) - offset, 0);
2452 909813 : if (!tab)
2453 0 : return NULL;
2454 909813 : tab->rational = ISL_F_ISSET(bset, ISL_BASIC_SET_RATIONAL);
2455 909813 : tab->cone = 1;
2456 :
2457 909813 : isl_int_init(cst);
2458 909813 : isl_int_set_si(cst, 0);
2459 5373568 : for (i = 0; i < bset->n_eq; ++i) {
2460 4463755 : isl_int_swap(bset->eq[i][offset], cst);
2461 4463755 : if (offset > 0) {
2462 0 : if (isl_tab_add_eq(tab, bset->eq[i] + offset) < 0)
2463 0 : goto error;
2464 : } else
2465 4463755 : tab = add_eq(tab, bset->eq[i]);
2466 4463755 : isl_int_swap(bset->eq[i][offset], cst);
2467 4463755 : if (!tab)
2468 0 : goto done;
2469 : }
2470 980750 : for (i = 0; i < bset->n_ineq; ++i) {
2471 : int r;
2472 70937 : isl_int_swap(bset->ineq[i][offset], cst);
2473 70937 : r = isl_tab_add_row(tab, bset->ineq[i] + offset);
2474 70937 : isl_int_swap(bset->ineq[i][offset], cst);
2475 70937 : if (r < 0)
2476 0 : goto error;
2477 70937 : tab->con[r].is_nonneg = 1;
2478 70937 : if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
2479 0 : goto error;
2480 : }
2481 : done:
2482 909813 : isl_int_clear(cst);
2483 909813 : return tab;
2484 : error:
2485 0 : isl_int_clear(cst);
2486 0 : isl_tab_free(tab);
2487 0 : return NULL;
2488 : }
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 909813 : isl_bool isl_tab_cone_is_bounded(struct isl_tab *tab)
2494 : {
2495 : int i;
2496 :
2497 909813 : if (!tab)
2498 0 : return isl_bool_error;
2499 909813 : if (tab->empty)
2500 0 : return isl_bool_true;
2501 909813 : if (tab->n_dead == tab->n_col)
2502 884106 : return isl_bool_true;
2503 :
2504 : for (;;) {
2505 233624 : for (i = tab->n_redundant; i < tab->n_row; ++i) {
2506 : struct isl_tab_var *var;
2507 : int sgn;
2508 116812 : var = isl_tab_var_from_row(tab, i);
2509 116812 : if (!var->is_nonneg)
2510 88103 : continue;
2511 28709 : sgn = sign_of_max(tab, var);
2512 28709 : if (sgn < -1)
2513 0 : return isl_bool_error;
2514 28709 : if (sgn != 0)
2515 0 : return isl_bool_false;
2516 28709 : if (close_row(tab, var, 0) < 0)
2517 0 : return isl_bool_error;
2518 28709 : break;
2519 : }
2520 28709 : if (tab->n_dead == tab->n_col)
2521 25707 : return isl_bool_true;
2522 3002 : if (i == tab->n_row)
2523 0 : return isl_bool_false;
2524 3002 : }
2525 : }
2526 :
2527 188854843 : int isl_tab_sample_is_integer(struct isl_tab *tab)
2528 : {
2529 : int i;
2530 :
2531 188854843 : if (!tab)
2532 0 : return -1;
2533 :
2534 725708940 : for (i = 0; i < tab->n_var; ++i) {
2535 : int row;
2536 649722286 : if (!tab->var[i].is_row)
2537 65246544 : continue;
2538 584475742 : row = tab->var[i].index;
2539 584475742 : if (!isl_int_is_divisible_by(tab->mat->row[row][1],
2540 : tab->mat->row[row][0]))
2541 112868189 : return 0;
2542 : }
2543 75986654 : return 1;
2544 : }
2545 :
2546 68904148 : static struct isl_vec *extract_integer_sample(struct isl_tab *tab)
2547 : {
2548 : int i;
2549 : struct isl_vec *vec;
2550 :
2551 68904148 : vec = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
2552 68904148 : if (!vec)
2553 0 : return NULL;
2554 :
2555 68904148 : isl_int_set_si(vec->block.data[0], 1);
2556 417960289 : for (i = 0; i < tab->n_var; ++i) {
2557 349056141 : if (!tab->var[i].is_row)
2558 60655936 : isl_int_set_si(vec->block.data[1 + i], 0);
2559 : else {
2560 288400205 : int row = tab->var[i].index;
2561 288400205 : isl_int_divexact(vec->block.data[1 + i],
2562 : tab->mat->row[row][1], tab->mat->row[row][0]);
2563 : }
2564 : }
2565 :
2566 68904148 : return vec;
2567 : }
2568 :
2569 46597212 : struct isl_vec *isl_tab_get_sample_value(struct isl_tab *tab)
2570 : {
2571 : int i;
2572 : struct isl_vec *vec;
2573 : isl_int m;
2574 :
2575 46597212 : if (!tab)
2576 0 : return NULL;
2577 :
2578 46597212 : vec = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
2579 46597212 : if (!vec)
2580 0 : return NULL;
2581 :
2582 46597212 : isl_int_init(m);
2583 :
2584 46597212 : isl_int_set_si(vec->block.data[0], 1);
2585 307386577 : for (i = 0; i < tab->n_var; ++i) {
2586 : int row;
2587 260789365 : if (!tab->var[i].is_row) {
2588 47866252 : isl_int_set_si(vec->block.data[1 + i], 0);
2589 47866252 : continue;
2590 : }
2591 212923113 : row = tab->var[i].index;
2592 212923113 : isl_int_gcd(m, vec->block.data[0], tab->mat->row[row][0]);
2593 212923113 : isl_int_divexact(m, tab->mat->row[row][0], m);
2594 212923113 : isl_seq_scale(vec->block.data, vec->block.data, m, 1 + i);
2595 212923113 : isl_int_divexact(m, vec->block.data[0], tab->mat->row[row][0]);
2596 212923113 : isl_int_mul(vec->block.data[1 + i], m, tab->mat->row[row][1]);
2597 : }
2598 46597212 : vec = isl_vec_normalize(vec);
2599 :
2600 46597212 : isl_int_clear(m);
2601 46597212 : return vec;
2602 : }
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 147512582 : static void get_rounded_sample_value(struct isl_tab *tab,
2608 : struct isl_tab_var *var, int sgn, isl_int *v)
2609 : {
2610 147512582 : if (!var->is_row)
2611 859648 : isl_int_set_si(*v, 0);
2612 146652934 : else if (sgn > 0)
2613 143983138 : isl_int_cdiv_q(*v, tab->mat->row[var->index][1],
2614 : tab->mat->row[var->index][0]);
2615 : else
2616 2669796 : isl_int_fdiv_q(*v, tab->mat->row[var->index][1],
2617 : tab->mat->row[var->index][0]);
2618 147512582 : }
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 3069314595 : struct isl_basic_map *isl_basic_map_update_from_tab(struct isl_basic_map *bmap,
2629 : struct isl_tab *tab)
2630 : {
2631 : int i;
2632 : unsigned n_eq;
2633 :
2634 3069314595 : if (!bmap)
2635 0 : return NULL;
2636 3069314595 : if (!tab)
2637 0 : return bmap;
2638 :
2639 3069314595 : n_eq = tab->n_eq;
2640 3069314595 : if (tab->empty)
2641 1926672 : bmap = isl_basic_map_set_to_empty(bmap);
2642 : else
2643 4724143616 : for (i = bmap->n_ineq - 1; i >= 0; --i) {
2644 1656755693 : if (isl_tab_is_equality(tab, n_eq + i))
2645 583854842 : isl_basic_map_inequality_to_equality(bmap, i);
2646 1072900851 : else if (isl_tab_is_redundant(tab, n_eq + i))
2647 259605480 : isl_basic_map_drop_inequality(bmap, i);
2648 : }
2649 3069314595 : if (bmap->n_eq != n_eq)
2650 48709663 : bmap = isl_basic_map_gauss(bmap, NULL);
2651 3069314595 : if (!tab->rational &&
2652 132209644 : bmap && !bmap->sample && isl_tab_sample_is_integer(tab))
2653 68904148 : bmap->sample = extract_integer_sample(tab);
2654 3069314595 : return bmap;
2655 : }
2656 :
2657 1065652 : struct isl_basic_set *isl_basic_set_update_from_tab(struct isl_basic_set *bset,
2658 : struct isl_tab *tab)
2659 : {
2660 1065652 : return bset_from_bmap(isl_basic_map_update_from_tab(bset_to_bmap(bset),
2661 : tab));
2662 : }
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 3262293 : static isl_stat drop_last_con_in_row(struct isl_tab *tab, int r)
2668 : {
2669 3262293 : if (!tab->con[r].is_row)
2670 0 : isl_die(isl_tab_get_ctx(tab), isl_error_internal,
2671 : "row unexpectedly moved to column",
2672 : return isl_stat_error);
2673 3262293 : if (r + 1 != tab->n_con)
2674 0 : isl_die(isl_tab_get_ctx(tab), isl_error_internal,
2675 : "additional constraints added", return isl_stat_error);
2676 3262293 : if (drop_row(tab, tab->con[r].index) < 0)
2677 0 : return isl_stat_error;
2678 :
2679 3262293 : return isl_stat_ok;
2680 : }
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 3262293 : static isl_stat cut_to_hyperplane(struct isl_tab *tab, struct isl_tab_var *var)
2693 : {
2694 : unsigned r;
2695 : isl_int *row;
2696 : int sgn;
2697 3262293 : unsigned off = 2 + tab->M;
2698 :
2699 3262293 : if (var->is_zero)
2700 0 : return isl_stat_ok;
2701 3262293 : if (var->is_redundant || !var->is_nonneg)
2702 0 : isl_die(isl_tab_get_ctx(tab), isl_error_invalid,
2703 : "expecting non-redundant non-negative variable",
2704 : return isl_stat_error);
2705 :
2706 3262293 : if (isl_tab_extend_cons(tab, 1) < 0)
2707 0 : return isl_stat_error;
2708 :
2709 3262293 : r = tab->n_con;
2710 3262293 : tab->con[r].index = tab->n_row;
2711 3262293 : tab->con[r].is_row = 1;
2712 3262293 : tab->con[r].is_nonneg = 0;
2713 3262293 : tab->con[r].is_zero = 0;
2714 3262293 : tab->con[r].is_redundant = 0;
2715 3262293 : tab->con[r].frozen = 0;
2716 3262293 : tab->con[r].negated = 0;
2717 3262293 : tab->row_var[tab->n_row] = ~r;
2718 3262293 : row = tab->mat->row[tab->n_row];
2719 :
2720 3262293 : if (var->is_row) {
2721 2918079 : isl_int_set(row[0], tab->mat->row[var->index][0]);
2722 5836158 : isl_seq_neg(row + 1,
2723 5836158 : tab->mat->row[var->index] + 1, 1 + tab->n_col);
2724 : } else {
2725 344214 : isl_int_set_si(row[0], 1);
2726 344214 : isl_seq_clr(row + 1, 1 + tab->n_col);
2727 344214 : isl_int_set_si(row[off + var->index], -1);
2728 : }
2729 :
2730 3262293 : tab->n_row++;
2731 3262293 : tab->n_con++;
2732 :
2733 3262293 : sgn = sign_of_max(tab, &tab->con[r]);
2734 3262293 : if (sgn < -1)
2735 0 : return isl_stat_error;
2736 3262293 : if (sgn < 0) {
2737 44558 : if (drop_last_con_in_row(tab, r) < 0)
2738 0 : return isl_stat_error;
2739 44558 : if (isl_tab_mark_empty(tab) < 0)
2740 0 : return isl_stat_error;
2741 44558 : return isl_stat_ok;
2742 : }
2743 3217735 : tab->con[r].is_nonneg = 1;
2744 : /* sgn == 0 */
2745 3217735 : if (close_row(tab, &tab->con[r], 1) < 0)
2746 0 : return isl_stat_error;
2747 3217735 : if (drop_last_con_in_row(tab, r) < 0)
2748 0 : return isl_stat_error;
2749 :
2750 3217735 : return isl_stat_ok;
2751 : }
2752 :
2753 : /* Check that "con" is a valid constraint position for "tab".
2754 : */
2755 1624792070 : static isl_stat isl_tab_check_con(struct isl_tab *tab, int con)
2756 : {
2757 1624792070 : if (!tab)
2758 0 : return isl_stat_error;
2759 1624792070 : if (con < 0 || con >= tab->n_con)
2760 0 : isl_die(isl_tab_get_ctx(tab), isl_error_invalid,
2761 : "position out of bounds", return isl_stat_error);
2762 1624792070 : return isl_stat_ok;
2763 : }
2764 :
2765 : /* Given a tableau "tab" and an inequality constraint "con" of the tableau,
2766 : * relax the inequality by one. That is, the inequality r >= 0 is replaced
2767 : * by r' = r + 1 >= 0.
2768 : * If r is a row variable, we simply increase the constant term by one
2769 : * (taking into account the denominator).
2770 : * If r is a column variable, then we need to modify each row that
2771 : * refers to r = r' - 1 by substituting this equality, effectively
2772 : * subtracting the coefficient of the column from the constant.
2773 : * We should only do this if the minimum is manifestly unbounded,
2774 : * however. Otherwise, we may end up with negative sample values
2775 : * for non-negative variables.
2776 : * So, if r is a column variable with a minimum that is not
2777 : * manifestly unbounded, then we need to move it to a row.
2778 : * However, the sample value of this row may be negative,
2779 : * even after the relaxation, so we need to restore it.
2780 : * We therefore prefer to pivot a column up to a row, if possible.
2781 : */
2782 380231 : int isl_tab_relax(struct isl_tab *tab, int con)
2783 : {
2784 : struct isl_tab_var *var;
2785 :
2786 380231 : if (!tab)
2787 0 : return -1;
2788 :
2789 380231 : var = &tab->con[con];
2790 :
2791 380231 : if (var->is_row && (var->index < 0 || var->index < tab->n_redundant))
2792 0 : isl_die(tab->mat->ctx, isl_error_invalid,
2793 : "cannot relax redundant constraint", return -1);
2794 380231 : if (!var->is_row && (var->index < 0 || var->index < tab->n_dead))
2795 0 : isl_die(tab->mat->ctx, isl_error_invalid,
2796 : "cannot relax dead constraint", return -1);
2797 :
2798 380231 : if (!var->is_row && !max_is_manifestly_unbounded(tab, var))
2799 36170 : if (to_row(tab, var, 1) < 0)
2800 0 : return -1;
2801 380231 : if (!var->is_row && !min_is_manifestly_unbounded(tab, var))
2802 4540 : if (to_row(tab, var, -1) < 0)
2803 0 : return -1;
2804 :
2805 380231 : if (var->is_row) {
2806 41891 : isl_int_add(tab->mat->row[var->index][1],
2807 : tab->mat->row[var->index][1], tab->mat->row[var->index][0]);
2808 41891 : if (restore_row(tab, var) < 0)
2809 0 : return -1;
2810 : } else {
2811 : int i;
2812 338340 : unsigned off = 2 + tab->M;
2813 :
2814 1293049 : for (i = 0; i < tab->n_row; ++i) {
2815 954709 : if (isl_int_is_zero(tab->mat->row[i][off + var->index]))
2816 576556 : continue;
2817 378153 : isl_int_sub(tab->mat->row[i][1], tab->mat->row[i][1],
2818 : tab->mat->row[i][off + var->index]);
2819 : }
2820 :
2821 : }
2822 :
2823 380231 : if (isl_tab_push_var(tab, isl_tab_undo_relax, var) < 0)
2824 0 : return -1;
2825 :
2826 380231 : return 0;
2827 : }
2828 :
2829 : /* Replace the variable v at position "pos" in the tableau "tab"
2830 : * by v' = v + shift.
2831 : *
2832 : * If the variable is in a column, then we first check if we can
2833 : * simply plug in v = v' - shift. The effect on a row with
2834 : * coefficient f/d for variable v is that the constant term c/d
2835 : * is replaced by (c - f * shift)/d. If shift is positive and
2836 : * f is negative for each row that needs to remain non-negative,
2837 : * then this is clearly safe. In other words, if the minimum of v
2838 : * is manifestly unbounded, then we can keep v in a column position.
2839 : * Otherwise, we can pivot it down to a row.
2840 : * Similarly, if shift is negative, we need to check if the maximum
2841 : * of is manifestly unbounded.
2842 : *
2843 : * If the variable is in a row (from the start or after pivoting),
2844 : * then the constant term c/d is replaced by (c + d * shift)/d.
2845 : */
2846 0 : int isl_tab_shift_var(struct isl_tab *tab, int pos, isl_int shift)
2847 : {
2848 : struct isl_tab_var *var;
2849 :
2850 0 : if (!tab)
2851 0 : return -1;
2852 0 : if (isl_int_is_zero(shift))
2853 0 : return 0;
2854 :
2855 0 : var = &tab->var[pos];
2856 0 : if (!var->is_row) {
2857 0 : if (isl_int_is_neg(shift)) {
2858 0 : if (!max_is_manifestly_unbounded(tab, var))
2859 0 : if (to_row(tab, var, 1) < 0)
2860 0 : return -1;
2861 : } else {
2862 0 : if (!min_is_manifestly_unbounded(tab, var))
2863 0 : if (to_row(tab, var, -1) < 0)
2864 0 : return -1;
2865 : }
2866 : }
2867 :
2868 0 : if (var->is_row) {
2869 0 : isl_int_addmul(tab->mat->row[var->index][1],
2870 : shift, tab->mat->row[var->index][0]);
2871 : } else {
2872 : int i;
2873 0 : unsigned off = 2 + tab->M;
2874 :
2875 0 : for (i = 0; i < tab->n_row; ++i) {
2876 0 : if (isl_int_is_zero(tab->mat->row[i][off + var->index]))
2877 0 : continue;
2878 0 : isl_int_submul(tab->mat->row[i][1],
2879 : shift, tab->mat->row[i][off + var->index]);
2880 : }
2881 :
2882 : }
2883 :
2884 0 : return 0;
2885 : }
2886 :
2887 : /* Remove the sign constraint from constraint "con".
2888 : *
2889 : * If the constraint variable was originally marked non-negative,
2890 : * then we make sure we mark it non-negative again during rollback.
2891 : */
2892 688378 : int isl_tab_unrestrict(struct isl_tab *tab, int con)
2893 : {
2894 : struct isl_tab_var *var;
2895 :
2896 688378 : if (!tab)
2897 0 : return -1;
2898 :
2899 688378 : var = &tab->con[con];
2900 688378 : if (!var->is_nonneg)
2901 0 : return 0;
2902 :
2903 688378 : var->is_nonneg = 0;
2904 688378 : if (isl_tab_push_var(tab, isl_tab_undo_unrestrict, var) < 0)
2905 0 : return -1;
2906 :
2907 688378 : return 0;
2908 : }
2909 :
2910 384257 : int isl_tab_select_facet(struct isl_tab *tab, int con)
2911 : {
2912 384257 : if (!tab)
2913 0 : return -1;
2914 :
2915 384257 : return cut_to_hyperplane(tab, &tab->con[con]);
2916 : }
2917 :
2918 4077465203 : static int may_be_equality(struct isl_tab *tab, int row)
2919 : {
2920 8213525812 : return tab->rational ? isl_int_is_zero(tab->mat->row[row][1])
2921 4136060609 : : isl_int_lt(tab->mat->row[row][1],
2922 : tab->mat->row[row][0]);
2923 : }
2924 :
2925 : /* Return an isl_tab_var that has been marked or NULL if no such
2926 : * variable can be found.
2927 : * The marked field has only been set for variables that
2928 : * appear in non-redundant rows or non-dead columns.
2929 : *
2930 : * Pick the last constraint variable that is marked and
2931 : * that appears in either a non-redundant row or a non-dead columns.
2932 : * Since the returned variable is tested for being a redundant constraint or
2933 : * an implicit equality, there is no need to return any tab variable that
2934 : * corresponds to a variable.
2935 : */
2936 1237786113 : static struct isl_tab_var *select_marked(struct isl_tab *tab)
2937 : {
2938 : int i;
2939 : struct isl_tab_var *var;
2940 :
2941 10894147036 : for (i = tab->n_con - 1; i >= 0; --i) {
2942 10819561321 : var = &tab->con[i];
2943 10819561321 : if (var->index < 0)
2944 5373514603 : continue;
2945 5446046718 : if (var->is_row && var->index < tab->n_redundant)
2946 1099995297 : continue;
2947 4346051421 : if (!var->is_row && var->index < tab->n_dead)
2948 15254 : continue;
2949 4346036167 : if (var->marked)
2950 1163200398 : return var;
2951 : }
2952 :
2953 74585715 : return NULL;
2954 : }
2955 :
2956 : /* Check for (near) equalities among the constraints.
2957 : * A constraint is an equality if it is non-negative and if
2958 : * its maximal value is either
2959 : * - zero (in case of rational tableaus), or
2960 : * - strictly less than 1 (in case of integer tableaus)
2961 : *
2962 : * We first mark all non-redundant and non-dead variables that
2963 : * are not frozen and not obviously not an equality.
2964 : * Then we iterate over all marked variables if they can attain
2965 : * any values larger than zero or at least one.
2966 : * If the maximal value is zero, we mark any column variables
2967 : * that appear in the row as being zero and mark the row as being redundant.
2968 : * Otherwise, if the maximal value is strictly less than one (and the
2969 : * tableau is integer), then we restrict the value to being zero
2970 : * by adding an opposite non-negative variable.
2971 : * The order in which the variables are considered is not important.
2972 : */
2973 3051342963 : int isl_tab_detect_implicit_equalities(struct isl_tab *tab)
2974 : {
2975 : int i;
2976 : unsigned n_marked;
2977 :
2978 3051342963 : if (!tab)
2979 0 : return -1;
2980 3051342963 : if (tab->empty)
2981 2078749 : return 0;
2982 3049264214 : if (tab->n_dead == tab->n_col)
2983 2901044469 : return 0;
2984 :
2985 148219745 : n_marked = 0;
2986 1632122324 : for (i = tab->n_redundant; i < tab->n_row; ++i) {
2987 1483902579 : struct isl_tab_var *var = isl_tab_var_from_row(tab, i);
2988 2754253383 : var->marked = !var->frozen && var->is_nonneg &&
2989 1270350804 : may_be_equality(tab, i);
2990 1483902579 : if (var->marked)
2991 772122446 : n_marked++;
2992 : }
2993 718536201 : for (i = tab->n_dead; i < tab->n_col; ++i) {
2994 570316456 : struct isl_tab_var *var = var_from_col(tab, i);
2995 570316456 : var->marked = !var->frozen && var->is_nonneg;
2996 570316456 : if (var->marked)
2997 191221269 : n_marked++;
2998 : }
2999 960297717 : while (n_marked) {
3000 : struct isl_tab_var *var;
3001 : int sgn;
3002 715938933 : var = select_marked(tab);
3003 715938933 : if (!var)
3004 49202670 : break;
3005 666736263 : var->marked = 0;
3006 666736263 : n_marked--;
3007 666736263 : sgn = sign_of_max(tab, var);
3008 666736263 : if (sgn < 0)
3009 0 : return -1;
3010 666736263 : if (sgn == 0) {
3011 383947784 : if (close_row(tab, var, 0) < 0)
3012 0 : return -1;
3013 282788479 : } else if (!tab->rational && !at_least_one(tab, var)) {
3014 2878036 : if (cut_to_hyperplane(tab, var) < 0)
3015 0 : return -1;
3016 2878036 : return isl_tab_detect_implicit_equalities(tab);
3017 : }
3018 7185446198 : for (i = tab->n_redundant; i < tab->n_row; ++i) {
3019 6521587971 : var = isl_tab_var_from_row(tab, i);
3020 6521587971 : if (!var->marked)
3021 3714473572 : continue;
3022 2807114399 : if (may_be_equality(tab, i))
3023 2792379314 : continue;
3024 14735085 : var->marked = 0;
3025 14735085 : n_marked--;
3026 : }
3027 : }
3028 :
3029 145341709 : return 0;
3030 : }
3031 :
3032 : /* Update the element of row_var or col_var that corresponds to
3033 : * constraint tab->con[i] to a move from position "old" to position "i".
3034 : */
3035 1544746 : static int update_con_after_move(struct isl_tab *tab, int i, int old)
3036 : {
3037 : int *p;
3038 : int index;
3039 :
3040 1544746 : index = tab->con[i].index;
3041 1544746 : if (index == -1)
3042 1118901 : return 0;
3043 425845 : p = tab->con[i].is_row ? tab->row_var : tab->col_var;
3044 425845 : if (p[index] != ~old)
3045 0 : isl_die(tab->mat->ctx, isl_error_internal,
3046 : "broken internal state", return -1);
3047 425845 : p[index] = ~i;
3048 :
3049 425845 : return 0;
3050 : }
3051 :
3052 : /* Interchange constraints "con1" and "con2" in "tab".
3053 : * In particular, interchange the contents of these entries in tab->con.
3054 : * Since tab->col_var and tab->row_var point back into this array,
3055 : * they need to be updated accordingly.
3056 : */
3057 63502 : isl_stat isl_tab_swap_constraints(struct isl_tab *tab, int con1, int con2)
3058 : {
3059 : struct isl_tab_var var;
3060 :
3061 127004 : if (isl_tab_check_con(tab, con1) < 0 ||
3062 63502 : isl_tab_check_con(tab, con2) < 0)
3063 0 : return isl_stat_error;
3064 :
3065 63502 : var = tab->con[con1];
3066 63502 : tab->con[con1] = tab->con[con2];
3067 63502 : if (update_con_after_move(tab, con1, con2) < 0)
3068 0 : return isl_stat_error;
3069 63502 : tab->con[con2] = var;
3070 63502 : if (update_con_after_move(tab, con2, con1) < 0)
3071 0 : return isl_stat_error;
3072 :
3073 63502 : return isl_stat_ok;
3074 : }
3075 :
3076 : /* Rotate the "n" constraints starting at "first" to the right,
3077 : * putting the last constraint in the position of the first constraint.
3078 : */
3079 345222 : static int rotate_constraints(struct isl_tab *tab, int first, int n)
3080 : {
3081 : int i, last;
3082 : struct isl_tab_var var;
3083 :
3084 345222 : if (n <= 1)
3085 113290 : return 0;
3086 :
3087 231932 : last = first + n - 1;
3088 231932 : var = tab->con[last];
3089 1274050 : for (i = last; i > first; --i) {
3090 1042118 : tab->con[i] = tab->con[i - 1];
3091 1042118 : if (update_con_after_move(tab, i, i - 1) < 0)
3092 0 : return -1;
3093 : }
3094 231932 : tab->con[first] = var;
3095 231932 : if (update_con_after_move(tab, first, last) < 0)
3096 0 : return -1;
3097 :
3098 231932 : return 0;
3099 : }
3100 :
3101 : /* Drop the "n" entries starting at position "first" in tab->con, moving all
3102 : * subsequent entries down.
3103 : * Since some of the entries of tab->row_var and tab->col_var contain
3104 : * indices into this array, they have to be updated accordingly.
3105 : */
3106 49454 : static isl_stat con_drop_entries(struct isl_tab *tab,
3107 : unsigned first, unsigned n)
3108 : {
3109 : int i;
3110 :
3111 49454 : if (first + n > tab->n_con || first + n < first)
3112 0 : isl_die(isl_tab_get_ctx(tab), isl_error_internal,
3113 : "invalid range", return isl_stat_error);
3114 :
3115 49454 : tab->n_con -= n;
3116 :
3117 193146 : for (i = first; i < tab->n_con; ++i) {
3118 143692 : tab->con[i] = tab->con[i + n];
3119 143692 : if (update_con_after_move(tab, i, i + n) < 0)
3120 0 : return isl_stat_error;
3121 : }
3122 :
3123 49454 : return isl_stat_ok;
3124 : }
3125 :
3126 : /* isl_basic_map_gauss5 callback that gets called when
3127 : * two (equality) constraints "a" and "b" get interchanged
3128 : * in the basic map. Perform the same interchange in "tab".
3129 : */
3130 63502 : static isl_stat swap_eq(unsigned a, unsigned b, void *user)
3131 : {
3132 63502 : struct isl_tab *tab = user;
3133 :
3134 63502 : return isl_tab_swap_constraints(tab, a, b);
3135 : }
3136 :
3137 : /* isl_basic_map_gauss5 callback that gets called when
3138 : * the final "n" equality constraints get removed.
3139 : * As a special case, if "n" is equal to the total number
3140 : * of equality constraints, then this means the basic map
3141 : * turned out to be empty.
3142 : * Drop the same number of equality constraints from "tab" or
3143 : * mark it empty in the special case.
3144 : */
3145 49454 : static isl_stat drop_eq(unsigned n, void *user)
3146 : {
3147 49454 : struct isl_tab *tab = user;
3148 :
3149 49454 : if (tab->n_eq == n)
3150 0 : return isl_tab_mark_empty(tab);
3151 :
3152 49454 : tab->n_eq -= n;
3153 49454 : return con_drop_entries(tab, tab->n_eq, n);
3154 : }
3155 :
3156 : /* If "bmap" has more than a single reference, then call
3157 : * isl_basic_map_gauss on it, updating "tab" accordingly.
3158 : */
3159 52565 : static __isl_give isl_basic_map *gauss_if_shared(__isl_take isl_basic_map *bmap,
3160 : struct isl_tab *tab)
3161 : {
3162 : isl_bool single;
3163 :
3164 52565 : single = isl_basic_map_has_single_reference(bmap);
3165 52565 : if (single < 0)
3166 0 : return isl_basic_map_free(bmap);
3167 52565 : if (single)
3168 499 : return bmap;
3169 52066 : return isl_basic_map_gauss5(bmap, NULL, &swap_eq, &drop_eq, tab);
3170 : }
3171 :
3172 : /* Make the equalities that are implicit in "bmap" but that have been
3173 : * detected in the corresponding "tab" explicit in "bmap" and update
3174 : * "tab" to reflect the new order of the constraints.
3175 : *
3176 : * In particular, if inequality i is an implicit equality then
3177 : * isl_basic_map_inequality_to_equality will move the inequality
3178 : * in front of the other equality and it will move the last inequality
3179 : * in the position of inequality i.
3180 : * In the tableau, the inequalities of "bmap" are stored after the equalities
3181 : * and so the original order
3182 : *
3183 : * E E E E E A A A I B B B B L
3184 : *
3185 : * is changed into
3186 : *
3187 : * I E E E E E A A A L B B B B
3188 : *
3189 : * where I is the implicit equality, the E are equalities,
3190 : * the A inequalities before I, the B inequalities after I and
3191 : * L the last inequality.
3192 : * We therefore need to rotate to the right two sets of constraints,
3193 : * those up to and including I and those after I.
3194 : *
3195 : * If "tab" contains any constraints that are not in "bmap" then they
3196 : * appear after those in "bmap" and they should be left untouched.
3197 : *
3198 : * Note that this function only calls isl_basic_map_gauss
3199 : * (in case some equality constraints got detected)
3200 : * if "bmap" has more than one reference.
3201 : * If it only has a single reference, then it is left in a temporary state,
3202 : * because the caller may require this state.
3203 : * Calling isl_basic_map_gauss is then the responsibility of the caller.
3204 : */
3205 2904988922 : __isl_give isl_basic_map *isl_tab_make_equalities_explicit(struct isl_tab *tab,
3206 : __isl_take isl_basic_map *bmap)
3207 : {
3208 : int i;
3209 : unsigned n_eq;
3210 :
3211 2904988922 : if (!tab || !bmap)
3212 0 : return isl_basic_map_free(bmap);
3213 2904988922 : if (tab->empty)
3214 38835 : return bmap;
3215 :
3216 2904950087 : n_eq = tab->n_eq;
3217 2942094715 : for (i = bmap->n_ineq - 1; i >= 0; --i) {
3218 37144628 : if (!isl_tab_is_equality(tab, bmap->n_eq + i))
3219 36972017 : continue;
3220 172611 : isl_basic_map_inequality_to_equality(bmap, i);
3221 172611 : if (rotate_constraints(tab, 0, tab->n_eq + i + 1) < 0)
3222 0 : return isl_basic_map_free(bmap);
3223 172611 : if (rotate_constraints(tab, tab->n_eq + i + 1,
3224 172611 : bmap->n_ineq - i) < 0)
3225 0 : return isl_basic_map_free(bmap);
3226 172611 : tab->n_eq++;
3227 : }
3228 :
3229 2904950087 : if (n_eq != tab->n_eq)
3230 52565 : bmap = gauss_if_shared(bmap, tab);
3231 :
3232 2904950087 : return bmap;
3233 : }
3234 :
3235 12000225019 : static int con_is_redundant(struct isl_tab *tab, struct isl_tab_var *var)
3236 : {
3237 12000225019 : if (!tab)
3238 0 : return -1;
3239 12000225019 : if (tab->rational) {
3240 11417362933 : int sgn = sign_of_min(tab, var);
3241 11417362933 : if (sgn < -1)
3242 0 : return -1;
3243 11417362933 : return sgn >= 0;
3244 : } else {
3245 582862086 : int irred = isl_tab_min_at_most_neg_one(tab, var);
3246 582862086 : if (irred < 0)
3247 0 : return -1;
3248 582862086 : return !irred;
3249 : }
3250 : }
3251 :
3252 : /* Check for (near) redundant constraints.
3253 : * A constraint is redundant if it is non-negative and if
3254 : * its minimal value (temporarily ignoring the non-negativity) is either
3255 : * - zero (in case of rational tableaus), or
3256 : * - strictly larger than -1 (in case of integer tableaus)
3257 : *
3258 : * We first mark all non-redundant and non-dead variables that
3259 : * are not frozen and not obviously negatively unbounded.
3260 : * Then we iterate over all marked variables if they can attain
3261 : * any values smaller than zero or at most negative one.
3262 : * If not, we mark the row as being redundant (assuming it hasn't
3263 : * been detected as being obviously redundant in the mean time).
3264 : */
3265 3001868136 : int isl_tab_detect_redundant(struct isl_tab *tab)
3266 : {
3267 : int i;
3268 : unsigned n_marked;
3269 :
3270 3001868136 : if (!tab)
3271 0 : return -1;
3272 3001868136 : if (tab->empty)
3273 38835 : return 0;
3274 3001829301 : if (tab->n_redundant == tab->n_row)
3275 260284 : return 0;
3276 :
3277 3001569017 : n_marked = 0;
3278 8830925903 : for (i = tab->n_redundant; i < tab->n_row; ++i) {
3279 5829356886 : struct isl_tab_var *var = isl_tab_var_from_row(tab, i);
3280 5829356886 : var->marked = !var->frozen && var->is_nonneg;
3281 5829356886 : if (var->marked)
3282 528600794 : n_marked++;
3283 : }
3284 3482882398 : for (i = tab->n_dead; i < tab->n_col; ++i) {
3285 481313381 : struct isl_tab_var *var = var_from_col(tab, i);
3286 663027342 : var->marked = !var->frozen && var->is_nonneg &&
3287 181713961 : !min_is_manifestly_unbounded(tab, var);
3288 481313381 : if (var->marked)
3289 114700640 : n_marked++;
3290 : }
3291 6499602169 : while (n_marked) {
3292 : struct isl_tab_var *var;
3293 : int red;
3294 521847180 : var = select_marked(tab);
3295 521847180 : if (!var)
3296 25383045 : break;
3297 496464135 : var->marked = 0;
3298 496464135 : n_marked--;
3299 496464135 : red = con_is_redundant(tab, var);
3300 496464135 : if (red < 0)
3301 0 : return -1;
3302 496464135 : if (red && !var->is_redundant)
3303 6256946 : if (isl_tab_mark_redundant(tab, var->index) < 0)
3304 0 : return -1;
3305 3020247138 : for (i = tab->n_dead; i < tab->n_col; ++i) {
3306 2523783003 : var = var_from_col(tab, i);
3307 2523783003 : if (!var->marked)
3308 1977360783 : continue;
3309 546422220 : if (!min_is_manifestly_unbounded(tab, var))
3310 518453656 : continue;
3311 27968564 : var->marked = 0;
3312 27968564 : n_marked--;
3313 : }
3314 : }
3315 :
3316 3001569017 : return 0;
3317 : }
3318 :
3319 1694280552 : int isl_tab_is_equality(struct isl_tab *tab, int con)
3320 : {
3321 : int row;
3322 : unsigned off;
3323 :
3324 1694280552 : if (!tab)
3325 0 : return -1;
3326 1694280552 : if (tab->con[con].is_zero)
3327 584027445 : return 1;
3328 1110253107 : if (tab->con[con].is_redundant)
3329 260380169 : return 0;
3330 849872938 : if (!tab->con[con].is_row)
3331 396497427 : return tab->con[con].index < tab->n_dead;
3332 :
3333 453375511 : row = tab->con[con].index;
3334 :
3335 453375511 : off = 2 + tab->M;
3336 1433133202 : return isl_int_is_zero(tab->mat->row[row][1]) &&
3337 600067175 : !row_is_big(tab, row) &&
3338 73345832 : isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
3339 73345832 : tab->n_col - tab->n_dead) == -1;
3340 : }
3341 :
3342 : /* Return the minimal value of the affine expression "f" with denominator
3343 : * "denom" in *opt, *opt_denom, assuming the tableau is not empty and
3344 : * the expression cannot attain arbitrarily small values.
3345 : * If opt_denom is NULL, then *opt is rounded up to the nearest integer.
3346 : * The return value reflects the nature of the result (empty, unbounded,
3347 : * minimal value returned in *opt).
3348 : *
3349 : * This function assumes that at least one more row and at least
3350 : * one more element in the constraint array are available in the tableau.
3351 : */
3352 2560315470 : enum isl_lp_result isl_tab_min(struct isl_tab *tab,
3353 : isl_int *f, isl_int denom, isl_int *opt, isl_int *opt_denom,
3354 : unsigned flags)
3355 : {
3356 : int r;
3357 2560315470 : enum isl_lp_result res = isl_lp_ok;
3358 : struct isl_tab_var *var;
3359 : struct isl_tab_undo *snap;
3360 :
3361 2560315470 : if (!tab)
3362 0 : return isl_lp_error;
3363 :
3364 2560315470 : if (tab->empty)
3365 0 : return isl_lp_empty;
3366 :
3367 2560315470 : snap = isl_tab_snap(tab);
3368 2560315470 : r = isl_tab_add_row(tab, f);
3369 2560315470 : if (r < 0)
3370 0 : return isl_lp_error;
3371 2560315470 : var = &tab->con[r];
3372 : for (;;) {
3373 : int row, col;
3374 4308142892 : find_pivot(tab, var, var, -1, &row, &col);
3375 4308142892 : if (row == var->index) {
3376 426448 : res = isl_lp_unbounded;
3377 2560741918 : break;
3378 : }
3379 4307716444 : if (row == -1)
3380 2559889022 : break;
3381 1747827422 : if (isl_tab_pivot(tab, row, col) < 0)
3382 0 : return isl_lp_error;
3383 1747827422 : }
3384 2560315470 : isl_int_mul(tab->mat->row[var->index][0],
3385 : tab->mat->row[var->index][0], denom);
3386 2560315470 : if (ISL_FL_ISSET(flags, ISL_TAB_SAVE_DUAL)) {
3387 : int i;
3388 :
3389 55978629 : isl_vec_free(tab->dual);
3390 55978629 : tab->dual = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_con);
3391 55978629 : if (!tab->dual)
3392 0 : return isl_lp_error;
3393 55978629 : isl_int_set(tab->dual->el[0], tab->mat->row[var->index][0]);
3394 2672641811 : for (i = 0; i < tab->n_con; ++i) {
3395 : int pos;
3396 2616663182 : if (tab->con[i].is_row) {
3397 2003034698 : isl_int_set_si(tab->dual->el[1 + i], 0);
3398 2003034698 : continue;
3399 : }
3400 613628484 : pos = 2 + tab->M + tab->con[i].index;
3401 613628484 : if (tab->con[i].negated)
3402 98396265 : isl_int_neg(tab->dual->el[1 + i],
3403 : tab->mat->row[var->index][pos]);
3404 : else
3405 515232219 : isl_int_set(tab->dual->el[1 + i],
3406 : tab->mat->row[var->index][pos]);
3407 : }
3408 : }
3409 2560315470 : if (opt && res == isl_lp_ok) {
3410 2559889022 : if (opt_denom) {
3411 2416647373 : isl_int_set(*opt, tab->mat->row[var->index][1]);
3412 2416647373 : isl_int_set(*opt_denom, tab->mat->row[var->index][0]);
3413 : } else
3414 143241649 : get_rounded_sample_value(tab, var, 1, opt);
3415 : }
3416 2560315470 : if (isl_tab_rollback(tab, snap) < 0)
3417 0 : return isl_lp_error;
3418 2560315470 : return res;
3419 : }
3420 :
3421 : /* Is the constraint at position "con" marked as being redundant?
3422 : * If it is marked as representing an equality, then it is not
3423 : * considered to be redundant.
3424 : * Note that isl_tab_mark_redundant marks both the isl_tab_var as
3425 : * redundant and moves the corresponding row into the first
3426 : * tab->n_redundant positions (or removes the row, assigning it index -1),
3427 : * so the final test is actually redundant itself.
3428 : */
3429 1624665066 : int isl_tab_is_redundant(struct isl_tab *tab, int con)
3430 : {
3431 1624665066 : if (isl_tab_check_con(tab, con) < 0)
3432 0 : return -1;
3433 1624665066 : if (tab->con[con].is_zero)
3434 0 : return 0;
3435 1624665066 : if (tab->con[con].is_redundant)
3436 645241482 : return 1;
3437 979423584 : return tab->con[con].is_row && tab->con[con].index < tab->n_redundant;
3438 : }
3439 :
3440 : /* Is variable "var" of "tab" fixed to a constant value by its row
3441 : * in the tableau?
3442 : * If so and if "value" is not NULL, then store this constant value
3443 : * in "value".
3444 : *
3445 : * That is, is it a row variable that only has non-zero coefficients
3446 : * for dead columns?
3447 : */
3448 3946123 : static isl_bool is_constant(struct isl_tab *tab, struct isl_tab_var *var,
3449 : isl_int *value)
3450 : {
3451 3946123 : unsigned off = 2 + tab->M;
3452 3946123 : isl_mat *mat = tab->mat;
3453 : int n;
3454 : int row;
3455 : int pos;
3456 :
3457 3946123 : if (!var->is_row)
3458 859648 : return isl_bool_false;
3459 3086475 : row = var->index;
3460 3086475 : if (row_is_big(tab, row))
3461 0 : return isl_bool_false;
3462 3086475 : n = tab->n_col - tab->n_dead;
3463 3086475 : pos = isl_seq_first_non_zero(mat->row[row] + off + tab->n_dead, n);
3464 3086475 : if (pos != -1)
3465 2666701 : return isl_bool_false;
3466 419774 : if (value)
3467 0 : isl_int_divexact(*value, mat->row[row][1], mat->row[row][0]);
3468 419774 : return isl_bool_true;
3469 : }
3470 :
3471 : /* Has the variable "var' of "tab" reached a value that is greater than
3472 : * or equal (if sgn > 0) or smaller than or equal (if sgn < 0) to "target"?
3473 : * "tmp" has been initialized by the caller and can be used
3474 : * to perform local computations.
3475 : *
3476 : * If the sample value involves the big parameter, then any value
3477 : * is reached.
3478 : * Otherwise check if n/d >= t, i.e., n >= d * t (if sgn > 0)
3479 : * or n/d <= t, i.e., n <= d * t (if sgn < 0).
3480 : */
3481 5792930 : static int reached(struct isl_tab *tab, struct isl_tab_var *var, int sgn,
3482 : isl_int target, isl_int *tmp)
3483 : {
3484 5792930 : if (row_is_big(tab, var->index))
3485 0 : return 1;
3486 5792930 : isl_int_mul(*tmp, tab->mat->row[var->index][0], target);
3487 5792930 : if (sgn > 0)
3488 1333495 : return isl_int_ge(tab->mat->row[var->index][1], *tmp);
3489 : else
3490 4459435 : return isl_int_le(tab->mat->row[var->index][1], *tmp);
3491 : }
3492 :
3493 : /* Can variable "var" of "tab" attain the value "target" by
3494 : * pivoting up (if sgn > 0) or down (if sgn < 0)?
3495 : * If not, then pivot up [down] to the greatest [smallest]
3496 : * rational value.
3497 : * "tmp" has been initialized by the caller and can be used
3498 : * to perform local computations.
3499 : *
3500 : * If the variable is manifestly unbounded in the desired direction,
3501 : * then it can attain any value.
3502 : * Otherwise, it can be moved to a row.
3503 : * Continue pivoting until the target is reached.
3504 : * If no more pivoting can be performed, the maximal [minimal]
3505 : * rational value has been reached and the target cannot be reached.
3506 : * If the variable would be pivoted into a manifestly unbounded column,
3507 : * then the target can be reached.
3508 : */
3509 4267838 : static isl_bool var_reaches(struct isl_tab *tab, struct isl_tab_var *var,
3510 : int sgn, isl_int target, isl_int *tmp)
3511 : {
3512 : int row, col;
3513 :
3514 4267838 : if (sgn < 0 && min_is_manifestly_unbounded(tab, var))
3515 839657 : return isl_bool_true;
3516 3428181 : if (sgn > 0 && max_is_manifestly_unbounded(tab, var))
3517 0 : return isl_bool_true;
3518 3428181 : if (to_row(tab, var, sgn) < 0)
3519 0 : return isl_bool_error;
3520 9221111 : while (!reached(tab, var, sgn, target, tmp)) {
3521 4198359 : find_pivot(tab, var, var, sgn, &row, &col);
3522 4198359 : if (row == -1)
3523 744584 : return isl_bool_false;
3524 3453775 : if (row == var->index)
3525 1089026 : return isl_bool_true;
3526 2364749 : if (isl_tab_pivot(tab, row, col) < 0)
3527 0 : return isl_bool_error;
3528 : }
3529 :
3530 1594571 : return isl_bool_true;
3531 : }
3532 :
3533 : /* Check if variable "var" of "tab" can only attain a single (integer)
3534 : * value, and, if so, add an equality constraint to fix the variable
3535 : * to this single value and store the result in "target".
3536 : * "target" and "tmp" have been initialized by the caller.
3537 : *
3538 : * Given the current sample value, round it down and check
3539 : * whether it is possible to attain a strictly smaller integer value.
3540 : * If so, the variable is not restricted to a single integer value.
3541 : * Otherwise, the search stops at the smallest rational value.
3542 : * Round up this value and check whether it is possible to attain
3543 : * a strictly greater integer value.
3544 : * If so, the variable is not restricted to a single integer value.
3545 : * Otherwise, the search stops at the greatest rational value.
3546 : * If rounding down this value yields a value that is different
3547 : * from rounding up the smallest rational value, then the variable
3548 : * cannot attain any integer value. Mark the tableau empty.
3549 : * Otherwise, add an equality constraint that fixes the variable
3550 : * to the single integer value found.
3551 : */
3552 3526349 : static isl_bool detect_constant_with_tmp(struct isl_tab *tab,
3553 : struct isl_tab_var *var, isl_int *target, isl_int *tmp)
3554 : {
3555 : isl_bool reached;
3556 : isl_vec *eq;
3557 : int pos;
3558 : isl_stat r;
3559 :
3560 3526349 : get_rounded_sample_value(tab, var, -1, target);
3561 3526349 : isl_int_sub_ui(*target, *target, 1);
3562 3526349 : reached = var_reaches(tab, var, -1, *target, tmp);
3563 3526349 : if (reached < 0 || reached)
3564 2784860 : return isl_bool_not(reached);
3565 741489 : get_rounded_sample_value(tab, var, 1, target);
3566 741489 : isl_int_add_ui(*target, *target, 1);
3567 741489 : reached = var_reaches(tab, var, 1, *target, tmp);
3568 741489 : if (reached < 0 || reached)
3569 738394 : return isl_bool_not(reached);
3570 3095 : get_rounded_sample_value(tab, var, -1, tmp);
3571 3095 : isl_int_sub_ui(*target, *target, 1);
3572 3095 : if (isl_int_ne(*target, *tmp)) {
3573 0 : if (isl_tab_mark_empty(tab) < 0)
3574 0 : return isl_bool_error;
3575 0 : return isl_bool_false;
3576 : }
3577 :
3578 3095 : if (isl_tab_extend_cons(tab, 1) < 0)
3579 0 : return isl_bool_error;
3580 3095 : eq = isl_vec_alloc(isl_tab_get_ctx(tab), 1 + tab->n_var);
3581 3095 : if (!eq)
3582 0 : return isl_bool_error;
3583 3095 : pos = var - tab->var;
3584 3095 : isl_seq_clr(eq->el + 1, tab->n_var);
3585 3095 : isl_int_set_si(eq->el[1 + pos], -1);
3586 3095 : isl_int_set(eq->el[0], *target);
3587 3095 : r = isl_tab_add_eq(tab, eq->el);
3588 3095 : isl_vec_free(eq);
3589 :
3590 3095 : return r < 0 ? isl_bool_error : isl_bool_true;
3591 : }
3592 :
3593 : /* Check if variable "var" of "tab" can only attain a single (integer)
3594 : * value, and, if so, add an equality constraint to fix the variable
3595 : * to this single value and store the result in "value" (if "value"
3596 : * is not NULL).
3597 : *
3598 : * If the current sample value involves the big parameter,
3599 : * then the variable cannot have a fixed integer value.
3600 : * If the variable is already fixed to a single value by its row, then
3601 : * there is no need to add another equality constraint.
3602 : *
3603 : * Otherwise, allocate some temporary variables and continue
3604 : * with detect_constant_with_tmp.
3605 : */
3606 3946123 : static isl_bool get_constant(struct isl_tab *tab, struct isl_tab_var *var,
3607 : isl_int *value)
3608 : {
3609 : isl_int target, tmp;
3610 : isl_bool is_cst;
3611 :
3612 3946123 : if (var->is_row && row_is_big(tab, var->index))
3613 0 : return isl_bool_false;
3614 3946123 : is_cst = is_constant(tab, var, value);
3615 3946123 : if (is_cst < 0 || is_cst)
3616 419774 : return is_cst;
3617 :
3618 3526349 : if (!value)
3619 3526349 : isl_int_init(target);
3620 3526349 : isl_int_init(tmp);
3621 :
3622 3526349 : is_cst = detect_constant_with_tmp(tab, var,
3623 : value ? value : &target, &tmp);
3624 :
3625 3526349 : isl_int_clear(tmp);
3626 3526349 : if (!value)
3627 3526349 : isl_int_clear(target);
3628 :
3629 3526349 : return is_cst;
3630 : }
3631 :
3632 : /* Check if variable "var" of "tab" can only attain a single (integer)
3633 : * value, and, if so, add an equality constraint to fix the variable
3634 : * to this single value and store the result in "value" (if "value"
3635 : * is not NULL).
3636 : *
3637 : * For rational tableaus, nothing needs to be done.
3638 : */
3639 0 : isl_bool isl_tab_is_constant(struct isl_tab *tab, int var, isl_int *value)
3640 : {
3641 0 : if (!tab)
3642 0 : return isl_bool_error;
3643 0 : if (var < 0 || var >= tab->n_var)
3644 0 : isl_die(isl_tab_get_ctx(tab), isl_error_invalid,
3645 : "position out of bounds", return isl_bool_error);
3646 0 : if (tab->rational)
3647 0 : return isl_bool_false;
3648 :
3649 0 : return get_constant(tab, &tab->var[var], value);
3650 : }
3651 :
3652 : /* Check if any of the variables of "tab" can only attain a single (integer)
3653 : * value, and, if so, add equality constraints to fix those variables
3654 : * to these single values.
3655 : *
3656 : * For rational tableaus, nothing needs to be done.
3657 : */
3658 688378 : isl_stat isl_tab_detect_constants(struct isl_tab *tab)
3659 : {
3660 : int i;
3661 :
3662 688378 : if (!tab)
3663 0 : return isl_stat_error;
3664 688378 : if (tab->rational)
3665 0 : return isl_stat_ok;
3666 :
3667 4634501 : for (i = 0; i < tab->n_var; ++i) {
3668 3946123 : if (get_constant(tab, &tab->var[i], NULL) < 0)
3669 0 : return isl_stat_error;
3670 : }
3671 :
3672 688378 : return isl_stat_ok;
3673 : }
3674 :
3675 : /* Take a snapshot of the tableau that can be restored by a call to
3676 : * isl_tab_rollback.
3677 : */
3678 22975919148 : struct isl_tab_undo *isl_tab_snap(struct isl_tab *tab)
3679 : {
3680 22975919148 : if (!tab)
3681 0 : return NULL;
3682 22975919148 : tab->need_undo = 1;
3683 22975919148 : return tab->top;
3684 : }
3685 :
3686 : /* Does "tab" need to keep track of undo information?
3687 : * That is, was a snapshot taken that may need to be restored?
3688 : */
3689 0 : isl_bool isl_tab_need_undo(struct isl_tab *tab)
3690 : {
3691 0 : if (!tab)
3692 0 : return isl_bool_error;
3693 :
3694 0 : return tab->need_undo;
3695 : }
3696 :
3697 : /* Remove all tracking of undo information from "tab", invalidating
3698 : * any snapshots that may have been taken of the tableau.
3699 : * Since all snapshots have been invalidated, there is also
3700 : * no need to start keeping track of undo information again.
3701 : */
3702 0 : void isl_tab_clear_undo(struct isl_tab *tab)
3703 : {
3704 0 : if (!tab)
3705 0 : return;
3706 :
3707 0 : free_undo(tab);
3708 0 : tab->need_undo = 0;
3709 : }
3710 :
3711 : /* Undo the operation performed by isl_tab_relax.
3712 : */
3713 : static isl_stat unrelax(struct isl_tab *tab, struct isl_tab_var *var)
3714 : WARN_UNUSED;
3715 370821 : static isl_stat unrelax(struct isl_tab *tab, struct isl_tab_var *var)
3716 : {
3717 370821 : unsigned off = 2 + tab->M;
3718 :
3719 370821 : if (!var->is_row && !max_is_manifestly_unbounded(tab, var))
3720 29748 : if (to_row(tab, var, 1) < 0)
3721 0 : return isl_stat_error;
3722 :
3723 370821 : if (var->is_row) {
3724 31032 : isl_int_sub(tab->mat->row[var->index][1],
3725 : tab->mat->row[var->index][1], tab->mat->row[var->index][0]);
3726 31032 : if (var->is_nonneg) {
3727 31032 : int sgn = restore_row(tab, var);
3728 31032 : isl_assert(tab->mat->ctx, sgn >= 0,
3729 : return isl_stat_error);
3730 : }
3731 : } else {
3732 : int i;
3733 :
3734 1308715 : for (i = 0; i < tab->n_row; ++i) {
3735 968926 : if (isl_int_is_zero(tab->mat->row[i][off + var->index]))
3736 579118 : continue;
3737 389808 : isl_int_add(tab->mat->row[i][1], tab->mat->row[i][1],
3738 : tab->mat->row[i][off + var->index]);
3739 : }
3740 :
3741 : }
3742 :
3743 370821 : return isl_stat_ok;
3744 : }
3745 :
3746 : /* Undo the operation performed by isl_tab_unrestrict.
3747 : *
3748 : * In particular, mark the variable as being non-negative and make
3749 : * sure the sample value respects this constraint.
3750 : */
3751 687508 : static isl_stat ununrestrict(struct isl_tab *tab, struct isl_tab_var *var)
3752 : {
3753 687508 : var->is_nonneg = 1;
3754 :
3755 687508 : if (var->is_row && restore_row(tab, var) < -1)
3756 0 : return isl_stat_error;
3757 :
3758 687508 : return isl_stat_ok;
3759 : }
3760 :
3761 : /* Unmark the last redundant row in "tab" as being redundant.
3762 : * This undoes part of the modifications performed by isl_tab_mark_redundant.
3763 : * In particular, remove the redundant mark and make
3764 : * sure the sample value respects the constraint again.
3765 : * A variable that is marked non-negative by isl_tab_mark_redundant
3766 : * is covered by a separate undo record.
3767 : */
3768 5314195918 : static isl_stat restore_last_redundant(struct isl_tab *tab)
3769 : {
3770 : struct isl_tab_var *var;
3771 :
3772 5314195918 : if (tab->n_redundant < 1)
3773 0 : isl_die(isl_tab_get_ctx(tab), isl_error_internal,
3774 : "no redundant rows", return isl_stat_error);
3775 :
3776 5314195918 : var = isl_tab_var_from_row(tab, tab->n_redundant - 1);
3777 5314195918 : var->is_redundant = 0;
3778 5314195918 : tab->n_redundant--;
3779 5314195918 : restore_row(tab, var);
3780 :
3781 5314195918 : return isl_stat_ok;
3782 : }
3783 :
3784 : static isl_stat perform_undo_var(struct isl_tab *tab, struct isl_tab_undo *undo)
3785 : WARN_UNUSED;
3786 47883103905 : static isl_stat perform_undo_var(struct isl_tab *tab, struct isl_tab_undo *undo)
3787 : {
3788 47883103905 : struct isl_tab_var *var = var_from_index(tab, undo->u.var_index);
3789 47883103905 : switch (undo->type) {
3790 : case isl_tab_undo_nonneg:
3791 6996032757 : var->is_nonneg = 0;
3792 6996032757 : break;
3793 : case isl_tab_undo_redundant:
3794 5228617723 : if (!var->is_row || var->index != tab->n_redundant - 1)
3795 0 : isl_die(isl_tab_get_ctx(tab), isl_error_internal,
3796 : "not undoing last redundant row",
3797 : return isl_stat_error);
3798 5228617723 : return restore_last_redundant(tab);
3799 : case isl_tab_undo_freeze:
3800 10548579195 : var->frozen = 0;
3801 10548579195 : break;
3802 : case isl_tab_undo_zero:
3803 46595572 : var->is_zero = 0;
3804 46595572 : if (!var->is_row)
3805 46595572 : tab->n_dead--;
3806 46595572 : break;
3807 : case isl_tab_undo_allocate:
3808 25062220329 : if (undo->u.var_index >= 0) {
3809 0 : isl_assert(tab->mat->ctx, !var->is_row,
3810 : return isl_stat_error);
3811 0 : return drop_col(tab, var->index);
3812 : }
3813 25062220329 : if (!var->is_row) {
3814 377932878 : if (!max_is_manifestly_unbounded(tab, var)) {
3815 292422437 : if (to_row(tab, var, 1) < 0)
3816 0 : return isl_stat_error;
3817 85510441 : } else if (!min_is_manifestly_unbounded(tab, var)) {
3818 38249104 : if (to_row(tab, var, -1) < 0)
3819 0 : return isl_stat_error;
3820 : } else
3821 47261337 : if (to_row(tab, var, 0) < 0)
3822 0 : return isl_stat_error;
3823 : }
3824 25062220329 : return drop_row(tab, var->index);
3825 : case isl_tab_undo_relax:
3826 370821 : return unrelax(tab, var);
3827 : case isl_tab_undo_unrestrict:
3828 687508 : return ununrestrict(tab, var);
3829 : default:
3830 0 : isl_die(tab->mat->ctx, isl_error_internal,
3831 : "perform_undo_var called on invalid undo record",
3832 : return isl_stat_error);
3833 : }
3834 :
3835 17591207524 : return isl_stat_ok;
3836 : }
3837 :
3838 : /* Restore all rows that have been marked redundant by isl_tab_mark_redundant
3839 : * and that have been preserved in the tableau.
3840 : * Note that isl_tab_mark_redundant may also have marked some variables
3841 : * as being non-negative before marking them redundant. These need
3842 : * to be removed as well as otherwise some constraints could end up
3843 : * getting marked redundant with respect to the variable.
3844 : */
3845 48749036 : isl_stat isl_tab_restore_redundant(struct isl_tab *tab)
3846 : {
3847 48749036 : if (!tab)
3848 0 : return isl_stat_error;
3849 :
3850 48749036 : if (tab->need_undo)
3851 0 : isl_die(isl_tab_get_ctx(tab), isl_error_invalid,
3852 : "manually restoring redundant constraints "
3853 : "interferes with undo history",
3854 : return isl_stat_error);
3855 :
3856 183076267 : while (tab->n_redundant > 0) {
3857 85578195 : if (tab->row_var[tab->n_redundant - 1] >= 0) {
3858 : struct isl_tab_var *var;
3859 :
3860 47664391 : var = isl_tab_var_from_row(tab, tab->n_redundant - 1);
3861 47664391 : var->is_nonneg = 0;
3862 : }
3863 85578195 : restore_last_redundant(tab);
3864 : }
3865 48749036 : return isl_stat_ok;
3866 : }
3867 :
3868 : /* Undo the addition of an integer division to the basic map representation
3869 : * of "tab" in position "pos".
3870 : */
3871 0 : static isl_stat drop_bmap_div(struct isl_tab *tab, int pos)
3872 : {
3873 : int off;
3874 :
3875 0 : off = tab->n_var - isl_basic_map_dim(tab->bmap, isl_dim_div);
3876 0 : if (isl_basic_map_drop_div(tab->bmap, pos - off) < 0)
3877 : return isl_stat_error;
3878 0 : if (tab->samples) {
3879 0 : tab->samples = isl_mat_drop_cols(tab->samples, 1 + pos, 1);
3880 0 : if (!tab->samples)
3881 0 : return isl_stat_error;
3882 : }
3883 :
3884 0 : return isl_stat_ok;
3885 : }
3886 :
3887 : /* Restore the tableau to the state where the basic variables
3888 : * are those in "col_var".
3889 : * We first construct a list of variables that are currently in
3890 : * the basis, but shouldn't. Then we iterate over all variables
3891 : * that should be in the basis and for each one that is currently
3892 : * not in the basis, we exchange it with one of the elements of the
3893 : * list constructed before.
3894 : * We can always find an appropriate variable to pivot with because
3895 : * the current basis is mapped to the old basis by a non-singular
3896 : * matrix and so we can never end up with a zero row.
3897 : */
3898 0 : static int restore_basis(struct isl_tab *tab, int *col_var)
3899 : {
3900 : int i, j;
3901 0 : int n_extra = 0;
3902 0 : int *extra = NULL; /* current columns that contain bad stuff */
3903 0 : unsigned off = 2 + tab->M;
3904 :
3905 0 : extra = isl_alloc_array(tab->mat->ctx, int, tab->n_col);
3906 0 : if (tab->n_col && !extra)
3907 0 : goto error;
3908 0 : for (i = 0; i < tab->n_col; ++i) {
3909 0 : for (j = 0; j < tab->n_col; ++j)
3910 0 : if (tab->col_var[i] == col_var[j])
3911 0 : break;
3912 0 : if (j < tab->n_col)
3913 0 : continue;
3914 0 : extra[n_extra++] = i;
3915 : }
3916 0 : for (i = 0; i < tab->n_col && n_extra > 0; ++i) {
3917 : struct isl_tab_var *var;
3918 : int row;
3919 :
3920 0 : for (j = 0; j < tab->n_col; ++j)
3921 0 : if (col_var[i] == tab->col_var[j])
3922 0 : break;
3923 0 : if (j < tab->n_col)
3924 0 : continue;
3925 0 : var = var_from_index(tab, col_var[i]);
3926 0 : row = var->index;
3927 0 : for (j = 0; j < n_extra; ++j)
3928 0 : if (!isl_int_is_zero(tab->mat->row[row][off+extra[j]]))
3929 : break;
3930 0 : isl_assert(tab->mat->ctx, j < n_extra, goto error);
3931 0 : if (isl_tab_pivot(tab, row, extra[j]) < 0)
3932 0 : goto error;
3933 0 : extra[j] = extra[--n_extra];
3934 : }
3935 :
3936 0 : free(extra);
3937 0 : return 0;
3938 : error:
3939 0 : free(extra);
3940 0 : return -1;
3941 : }
3942 :
3943 : /* Remove all samples with index n or greater, i.e., those samples
3944 : * that were added since we saved this number of samples in
3945 : * isl_tab_save_samples.
3946 : */
3947 0 : static void drop_samples_since(struct isl_tab *tab, int n)
3948 : {
3949 : int i;
3950 :
3951 0 : for (i = tab->n_sample - 1; i >= 0 && tab->n_sample > n; --i) {
3952 0 : if (tab->sample_index[i] < n)
3953 0 : continue;
3954 :
3955 0 : if (i != tab->n_sample - 1) {
3956 0 : int t = tab->sample_index[tab->n_sample-1];
3957 0 : tab->sample_index[tab->n_sample-1] = tab->sample_index[i];
3958 0 : tab->sample_index[i] = t;
3959 0 : isl_mat_swap_rows(tab->samples, tab->n_sample-1, i);
3960 : }
3961 0 : tab->n_sample--;
3962 : }
3963 0 : }
3964 :
3965 : static isl_stat perform_undo(struct isl_tab *tab, struct isl_tab_undo *undo)
3966 : WARN_UNUSED;
3967 52372576986 : static isl_stat perform_undo(struct isl_tab *tab, struct isl_tab_undo *undo)
3968 : {
3969 52372576986 : switch (undo->type) {
3970 : case isl_tab_undo_rational:
3971 10639 : tab->rational = 0;
3972 10639 : break;
3973 : case isl_tab_undo_empty:
3974 749697066 : tab->empty = 0;
3975 749697066 : break;
3976 : case isl_tab_undo_nonneg:
3977 : case isl_tab_undo_redundant:
3978 : case isl_tab_undo_freeze:
3979 : case isl_tab_undo_zero:
3980 : case isl_tab_undo_allocate:
3981 : case isl_tab_undo_relax:
3982 : case isl_tab_undo_unrestrict:
3983 47883103905 : return perform_undo_var(tab, undo);
3984 : case isl_tab_undo_bmap_eq:
3985 0 : return isl_basic_map_free_equality(tab->bmap, 1);
3986 : case isl_tab_undo_bmap_ineq:
3987 3739765376 : return isl_basic_map_free_inequality(tab->bmap, 1);
3988 : case isl_tab_undo_bmap_div:
3989 0 : return drop_bmap_div(tab, undo->u.var_index);
3990 : case isl_tab_undo_saved_basis:
3991 0 : if (restore_basis(tab, undo->u.col_var) < 0)
3992 0 : return isl_stat_error;
3993 0 : break;
3994 : case isl_tab_undo_drop_sample:
3995 0 : tab->n_outside--;
3996 0 : break;
3997 : case isl_tab_undo_saved_samples:
3998 0 : drop_samples_since(tab, undo->u.n);
3999 0 : break;
4000 : case isl_tab_undo_callback:
4001 0 : return undo->u.callback->run(undo->u.callback);
4002 : default:
4003 0 : isl_assert(tab->mat->ctx, 0, return isl_stat_error);
4004 : }
4005 749707705 : return isl_stat_ok;
4006 : }
4007 :
4008 : /* Return the tableau to the state it was in when the snapshot "snap"
4009 : * was taken.
4010 : */
4011 22138649812 : int isl_tab_rollback(struct isl_tab *tab, struct isl_tab_undo *snap)
4012 : {
4013 : struct isl_tab_undo *undo, *next;
4014 :
4015 22138649812 : if (!tab)
4016 0 : return -1;
4017 :
4018 22138649812 : tab->in_undo = 1;
4019 74511226798 : for (undo = tab->top; undo && undo != &tab->bottom; undo = next) {
4020 53168171988 : next = undo->next;
4021 53168171988 : if (undo == snap)
4022 795595002 : break;
4023 52372576986 : if (perform_undo(tab, undo) < 0) {
4024 0 : tab->top = undo;
4025 0 : free_undo(tab);
4026 0 : tab->in_undo = 0;
4027 0 : return -1;
4028 : }
4029 52372576986 : free_undo_record(undo);
4030 : }
4031 22138649812 : tab->in_undo = 0;
4032 22138649812 : tab->top = undo;
4033 22138649812 : if (!undo)
4034 0 : return -1;
4035 22138649812 : return 0;
4036 : }
4037 :
4038 : /* The given row "row" represents an inequality violated by all
4039 : * points in the tableau. Check for some special cases of such
4040 : * separating constraints.
4041 : * In particular, if the row has been reduced to the constant -1,
4042 : * then we know the inequality is adjacent (but opposite) to
4043 : * an equality in the tableau.
4044 : * If the row has been reduced to r = c*(-1 -r'), with r' an inequality
4045 : * of the tableau and c a positive constant, then the inequality
4046 : * is adjacent (but opposite) to the inequality r'.
4047 : */
4048 7196337328 : static enum isl_ineq_type separation_type(struct isl_tab *tab, unsigned row)
4049 : {
4050 : int pos;
4051 7196337328 : unsigned off = 2 + tab->M;
4052 :
4053 7196337328 : if (tab->rational)
4054 7160499218 : return isl_ineq_separate;
4055 :
4056 35838110 : if (!isl_int_is_one(tab->mat->row[row][0]))
4057 2292298 : return isl_ineq_separate;
4058 :
4059 33545812 : pos = isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
4060 33545812 : tab->n_col - tab->n_dead);
4061 33545812 : if (pos == -1) {
4062 27028711 : if (isl_int_is_negone(tab->mat->row[row][1]))
4063 1968157 : return isl_ineq_adj_eq;
4064 : else
4065 25060554 : return isl_ineq_separate;
4066 : }
4067 :
4068 6517101 : if (!isl_int_eq(tab->mat->row[row][1],
4069 : tab->mat->row[row][off + tab->n_dead + pos]))
4070 763248 : return isl_ineq_separate;
4071 :
4072 11507706 : pos = isl_seq_first_non_zero(
4073 5753853 : tab->mat->row[row] + off + tab->n_dead + pos + 1,
4074 5753853 : tab->n_col - tab->n_dead - pos - 1);
4075 :
4076 5753853 : return pos == -1 ? isl_ineq_adj_ineq : isl_ineq_separate;
4077 : }
4078 :
4079 : /* Check the effect of inequality "ineq" on the tableau "tab".
4080 : * The result may be
4081 : * isl_ineq_redundant: satisfied by all points in the tableau
4082 : * isl_ineq_separate: satisfied by no point in the tableau
4083 : * isl_ineq_cut: satisfied by some by not all points
4084 : * isl_ineq_adj_eq: adjacent to an equality
4085 : * isl_ineq_adj_ineq: adjacent to an inequality.
4086 : */
4087 18712039164 : enum isl_ineq_type isl_tab_ineq_type(struct isl_tab *tab, isl_int *ineq)
4088 : {
4089 18712039164 : enum isl_ineq_type type = isl_ineq_error;
4090 18712039164 : struct isl_tab_undo *snap = NULL;
4091 : int con;
4092 : int row;
4093 :
4094 18712039164 : if (!tab)
4095 0 : return isl_ineq_error;
4096 :
4097 18712039164 : if (isl_tab_extend_cons(tab, 1) < 0)
4098 0 : return isl_ineq_error;
4099 :
4100 18712039164 : snap = isl_tab_snap(tab);
4101 :
4102 18712039164 : con = isl_tab_add_row(tab, ineq);
4103 18712039164 : if (con < 0)
4104 0 : goto error;
4105 :
4106 18712039164 : row = tab->con[con].index;
4107 18712039164 : if (isl_tab_row_is_redundant(tab, row))
4108 0 : type = isl_ineq_redundant;
4109 25920564982 : else if (isl_int_is_neg(tab->mat->row[row][1]) &&
4110 7252340012 : (tab->rational ||
4111 43814194 : isl_int_abs_ge(tab->mat->row[row][1],
4112 7208278280 : tab->mat->row[row][0]))) {
4113 7208278280 : int nonneg = at_least_zero(tab, &tab->con[con]);
4114 7208278280 : if (nonneg < 0)
4115 0 : goto error;
4116 7208278280 : if (nonneg)
4117 11940952 : type = isl_ineq_cut;
4118 : else
4119 7196337328 : type = separation_type(tab, row);
4120 : } else {
4121 11503760884 : int red = con_is_redundant(tab, &tab->con[con]);
4122 11503760884 : if (red < 0)
4123 0 : goto error;
4124 11503760884 : if (!red)
4125 24819776 : type = isl_ineq_cut;
4126 : else
4127 11478941108 : type = isl_ineq_redundant;
4128 : }
4129 :
4130 18712039164 : if (isl_tab_rollback(tab, snap))
4131 0 : return isl_ineq_error;
4132 18712039164 : return type;
4133 : error:
4134 0 : return isl_ineq_error;
4135 : }
4136 :
4137 55812133 : isl_stat isl_tab_track_bmap(struct isl_tab *tab, __isl_take isl_basic_map *bmap)
4138 : {
4139 55812133 : bmap = isl_basic_map_cow(bmap);
4140 55812133 : if (!tab || !bmap)
4141 : goto error;
4142 :
4143 55812133 : if (tab->empty) {
4144 27263921 : bmap = isl_basic_map_set_to_empty(bmap);
4145 27263921 : if (!bmap)
4146 0 : goto error;
4147 27263921 : tab->bmap = bmap;
4148 27263921 : return isl_stat_ok;
4149 : }
4150 :
4151 28548212 : isl_assert(tab->mat->ctx, tab->n_eq == bmap->n_eq, goto error);
4152 28548212 : isl_assert(tab->mat->ctx,
4153 : tab->n_con == bmap->n_eq + bmap->n_ineq, goto error);
4154 :
4155 28548212 : tab->bmap = bmap;
4156 :
4157 28548212 : return isl_stat_ok;
4158 : error:
4159 0 : isl_basic_map_free(bmap);
4160 0 : return isl_stat_error;
4161 : }
4162 :
4163 0 : isl_stat isl_tab_track_bset(struct isl_tab *tab, __isl_take isl_basic_set *bset)
4164 : {
4165 0 : return isl_tab_track_bmap(tab, bset_to_bmap(bset));
4166 : }
4167 :
4168 1123714 : __isl_keep isl_basic_set *isl_tab_peek_bset(struct isl_tab *tab)
4169 : {
4170 1123714 : if (!tab)
4171 0 : return NULL;
4172 :
4173 1123714 : return bset_from_bmap(tab->bmap);
4174 : }
4175 :
4176 0 : static void isl_tab_print_internal(__isl_keep struct isl_tab *tab,
4177 : FILE *out, int indent)
4178 : {
4179 : unsigned r, c;
4180 : int i;
4181 :
4182 0 : if (!tab) {
4183 0 : fprintf(out, "%*snull tab\n", indent, "");
4184 0 : return;
4185 : }
4186 0 : fprintf(out, "%*sn_redundant: %d, n_dead: %d", indent, "",
4187 : tab->n_redundant, tab->n_dead);
4188 0 : if (tab->rational)
4189 0 : fprintf(out, ", rational");
4190 0 : if (tab->empty)
4191 0 : fprintf(out, ", empty");
4192 0 : fprintf(out, "\n");
4193 0 : fprintf(out, "%*s[", indent, "");
4194 0 : for (i = 0; i < tab->n_var; ++i) {
4195 0 : if (i)
4196 0 : fprintf(out, (i == tab->n_param ||
4197 0 : i == tab->n_var - tab->n_div) ? "; "
4198 : : ", ");
4199 0 : fprintf(out, "%c%d%s", tab->var[i].is_row ? 'r' : 'c',
4200 0 : tab->var[i].index,
4201 0 : tab->var[i].is_zero ? " [=0]" :
4202 0 : tab->var[i].is_redundant ? " [R]" : "");
4203 : }
4204 0 : fprintf(out, "]\n");
4205 0 : fprintf(out, "%*s[", indent, "");
4206 0 : for (i = 0; i < tab->n_con; ++i) {
4207 0 : if (i)
4208 0 : fprintf(out, ", ");
4209 0 : fprintf(out, "%c%d%s", tab->con[i].is_row ? 'r' : 'c',
4210 0 : tab->con[i].index,
4211 0 : tab->con[i].is_zero ? " [=0]" :
4212 0 : tab->con[i].is_redundant ? " [R]" : "");
4213 : }
4214 0 : fprintf(out, "]\n");
4215 0 : fprintf(out, "%*s[", indent, "");
4216 0 : for (i = 0; i < tab->n_row; ++i) {
4217 0 : const char *sign = "";
4218 0 : if (i)
4219 0 : fprintf(out, ", ");
4220 0 : if (tab->row_sign) {
4221 0 : if (tab->row_sign[i] == isl_tab_row_unknown)
4222 0 : sign = "?";
4223 0 : else if (tab->row_sign[i] == isl_tab_row_neg)
4224 0 : sign = "-";
4225 0 : else if (tab->row_sign[i] == isl_tab_row_pos)
4226 0 : sign = "+";
4227 : else
4228 0 : sign = "+-";
4229 : }
4230 0 : fprintf(out, "r%d: %d%s%s", i, tab->row_var[i],
4231 0 : isl_tab_var_from_row(tab, i)->is_nonneg ? " [>=0]" : "", sign);
4232 : }
4233 0 : fprintf(out, "]\n");
4234 0 : fprintf(out, "%*s[", indent, "");
4235 0 : for (i = 0; i < tab->n_col; ++i) {
4236 0 : if (i)
4237 0 : fprintf(out, ", ");
4238 0 : fprintf(out, "c%d: %d%s", i, tab->col_var[i],
4239 0 : var_from_col(tab, i)->is_nonneg ? " [>=0]" : "");
4240 : }
4241 0 : fprintf(out, "]\n");
4242 0 : r = tab->mat->n_row;
4243 0 : tab->mat->n_row = tab->n_row;
4244 0 : c = tab->mat->n_col;
4245 0 : tab->mat->n_col = 2 + tab->M + tab->n_col;
4246 0 : isl_mat_print_internal(tab->mat, out, indent);
4247 0 : tab->mat->n_row = r;
4248 0 : tab->mat->n_col = c;
4249 0 : if (tab->bmap)
4250 0 : isl_basic_map_print_internal(tab->bmap, out, indent);
4251 : }
4252 :
4253 0 : void isl_tab_dump(__isl_keep struct isl_tab *tab)
4254 : {
4255 0 : isl_tab_print_internal(tab, stderr, 0);
4256 0 : }
|
