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 8317048188 : 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 8317048188 : unsigned off = 2 + M;
39 :
40 8317048188 : tab = isl_calloc_type(ctx, struct isl_tab);
41 8317048188 : if (!tab)
42 0 : return NULL;
43 8317048188 : tab->mat = isl_mat_alloc(ctx, n_row, off + n_var);
44 8317048188 : if (!tab->mat)
45 0 : goto error;
46 8317048188 : tab->var = isl_alloc_array(ctx, struct isl_tab_var, n_var);
47 8317048188 : if (n_var && !tab->var)
48 0 : goto error;
49 8317048188 : tab->con = isl_alloc_array(ctx, struct isl_tab_var, n_row);
50 8317048188 : if (n_row && !tab->con)
51 0 : goto error;
52 8317048188 : tab->col_var = isl_alloc_array(ctx, int, n_var);
53 8317048188 : if (n_var && !tab->col_var)
54 0 : goto error;
55 8317048188 : tab->row_var = isl_alloc_array(ctx, int, n_row);
56 8317048188 : if (n_row && !tab->row_var)
57 0 : goto error;
58 43526066917 : for (i = 0; i < n_var; ++i) {
59 35209018729 : tab->var[i].index = i;
60 35209018729 : tab->var[i].is_row = 0;
61 35209018729 : tab->var[i].is_nonneg = 0;
62 35209018729 : tab->var[i].is_zero = 0;
63 35209018729 : tab->var[i].is_redundant = 0;
64 35209018729 : tab->var[i].frozen = 0;
65 35209018729 : tab->var[i].negated = 0;
66 35209018729 : tab->col_var[i] = i;
67 : }
68 8317048188 : tab->n_row = 0;
69 8317048188 : tab->n_con = 0;
70 8317048188 : tab->n_eq = 0;
71 8317048188 : tab->max_con = n_row;
72 8317048188 : tab->n_col = n_var;
73 8317048188 : tab->n_var = n_var;
74 8317048188 : tab->max_var = n_var;
75 8317048188 : tab->n_param = 0;
76 8317048188 : tab->n_div = 0;
77 8317048188 : tab->n_dead = 0;
78 8317048188 : tab->n_redundant = 0;
79 8317048188 : tab->strict_redundant = 0;
80 8317048188 : tab->need_undo = 0;
81 8317048188 : tab->rational = 0;
82 8317048188 : tab->empty = 0;
83 8317048188 : tab->in_undo = 0;
84 8317048188 : tab->M = M;
85 8317048188 : tab->cone = 0;
86 8317048188 : tab->bottom.type = isl_tab_undo_bottom;
87 8317048188 : tab->bottom.next = NULL;
88 8317048188 : tab->top = &tab->bottom;
89 :
90 8317048188 : tab->n_zero = 0;
91 8317048188 : tab->n_unbounded = 0;
92 8317048188 : tab->basis = NULL;
93 :
94 8317048188 : return tab;
95 : error:
96 0 : isl_tab_free(tab);
97 0 : return NULL;
98 : }
99 :
100 44493435284 : isl_ctx *isl_tab_get_ctx(struct isl_tab *tab)
101 : {
102 44493435284 : return tab ? isl_mat_get_ctx(tab->mat) : NULL;
103 : }
104 :
105 25496805442 : int isl_tab_extend_cons(struct isl_tab *tab, unsigned n_new)
106 : {
107 : unsigned off;
108 :
109 25496805442 : if (!tab)
110 0 : return -1;
111 :
112 25496805442 : off = 2 + tab->M;
113 :
114 25496805442 : if (tab->max_con < tab->n_con + n_new) {
115 : struct isl_tab_var *con;
116 :
117 36534935 : con = isl_realloc_array(tab->mat->ctx, tab->con,
118 : struct isl_tab_var, tab->max_con + n_new);
119 36534935 : if (!con)
120 0 : return -1;
121 36534935 : tab->con = con;
122 36534935 : tab->max_con += n_new;
123 : }
124 25496805442 : if (tab->mat->n_row < tab->n_row + n_new) {
125 : int *row_var;
126 :
127 204864502 : tab->mat = isl_mat_extend(tab->mat,
128 204864502 : tab->n_row + n_new, off + tab->n_col);
129 102432251 : if (!tab->mat)
130 0 : return -1;
131 102432251 : row_var = isl_realloc_array(tab->mat->ctx, tab->row_var,
132 : int, tab->mat->n_row);
133 102432251 : if (!row_var)
134 0 : return -1;
135 102432251 : tab->row_var = row_var;
136 102432251 : 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 25496805442 : 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 >17771*10^7 : static void free_undo_record(struct isl_tab_undo *undo)
183 : {
184 >17771*10^7 : switch (undo->type) {
185 : case isl_tab_undo_saved_basis:
186 0 : free(undo->u.col_var);
187 0 : break;
188 : default:;
189 : }
190 >17771*10^7 : free(undo);
191 >17771*10^7 : }
192 :
193 8320351023 : static void free_undo(struct isl_tab *tab)
194 : {
195 : struct isl_tab_undo *undo, *next;
196 :
197 9652611126 : for (undo = tab->top; undo && undo != &tab->bottom; undo = next) {
198 1332260103 : next = undo->next;
199 1332260103 : free_undo_record(undo);
200 : }
201 8320351023 : tab->top = undo;
202 8320351023 : }
203 :
204 8321005721 : void isl_tab_free(struct isl_tab *tab)
205 : {
206 8321005721 : if (!tab)
207 654698 : return;
208 8320351023 : free_undo(tab);
209 8320351023 : isl_mat_free(tab->mat);
210 8320351023 : isl_vec_free(tab->dual);
211 8320351023 : isl_basic_map_free(tab->bmap);
212 8320351023 : free(tab->var);
213 8320351023 : free(tab->con);
214 8320351023 : free(tab->row_var);
215 8320351023 : free(tab->col_var);
216 8320351023 : free(tab->row_sign);
217 8320351023 : isl_mat_free(tab->samples);
218 8320351023 : free(tab->sample_index);
219 8320351023 : isl_mat_free(tab->basis);
220 8320351023 : 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 3302835 : 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 3302835 : prod = isl_mat_alloc(mat1->ctx, mat1->n_row + mat2->n_row,
333 3302835 : off + col1 + col2);
334 3302835 : if (!prod)
335 0 : return NULL;
336 :
337 3302835 : n = 0;
338 27671870 : for (i = 0; i < r1; ++i) {
339 24369035 : isl_seq_cpy(prod->row[n + i], mat1->row[i], off + d1);
340 24369035 : isl_seq_clr(prod->row[n + i] + off + d1, d2);
341 48738070 : isl_seq_cpy(prod->row[n + i] + off + d1 + d2,
342 24369035 : mat1->row[i] + off + d1, col1 - d1);
343 24369035 : isl_seq_clr(prod->row[n + i] + off + col1 + d1, col2 - d2);
344 : }
345 :
346 3302835 : n += r1;
347 27671870 : for (i = 0; i < r2; ++i) {
348 24369035 : isl_seq_cpy(prod->row[n + i], mat2->row[i], off);
349 24369035 : isl_seq_clr(prod->row[n + i] + off, d1);
350 48738070 : isl_seq_cpy(prod->row[n + i] + off + d1,
351 48738070 : mat2->row[i] + off, d2);
352 24369035 : isl_seq_clr(prod->row[n + i] + off + d1 + d2, col1 - d1);
353 48738070 : isl_seq_cpy(prod->row[n + i] + off + col1 + d1,
354 24369035 : mat2->row[i] + off + d2, col2 - d2);
355 : }
356 :
357 3302835 : n += r2;
358 42566725 : for (i = 0; i < row1 - r1; ++i) {
359 39263890 : isl_seq_cpy(prod->row[n + i], mat1->row[r1 + i], off + d1);
360 39263890 : isl_seq_clr(prod->row[n + i] + off + d1, d2);
361 78527780 : isl_seq_cpy(prod->row[n + i] + off + d1 + d2,
362 39263890 : mat1->row[r1 + i] + off + d1, col1 - d1);
363 39263890 : isl_seq_clr(prod->row[n + i] + off + col1 + d1, col2 - d2);
364 : }
365 :
366 3302835 : n += row1 - r1;
367 42566725 : for (i = 0; i < row2 - r2; ++i) {
368 39263890 : isl_seq_cpy(prod->row[n + i], mat2->row[r2 + i], off);
369 39263890 : isl_seq_clr(prod->row[n + i] + off, d1);
370 78527780 : isl_seq_cpy(prod->row[n + i] + off + d1,
371 78527780 : mat2->row[r2 + i] + off, d2);
372 39263890 : isl_seq_clr(prod->row[n + i] + off + d1 + d2, col1 - d1);
373 78527780 : isl_seq_cpy(prod->row[n + i] + off + col1 + d1,
374 39263890 : mat2->row[r2 + i] + off + d2, col2 - d2);
375 : }
376 :
377 3302835 : 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 80491241 : static void update_index1(struct isl_tab_var *var,
384 : unsigned r1, unsigned r2, unsigned d1, unsigned d2)
385 : {
386 80491241 : if (var->index == -1)
387 440441 : return;
388 80050800 : if (var->is_row && var->index >= r1)
389 39263890 : var->index += r2;
390 80050800 : if (!var->is_row && var->index >= d1)
391 15226941 : 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 80491241 : 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 80491241 : if (var->index == -1)
402 440441 : return;
403 80050800 : if (var->is_row) {
404 63632925 : if (var->index < r2)
405 24369035 : var->index += r1;
406 : else
407 39263890 : var->index += row1;
408 : } else {
409 16417875 : if (var->index < d2)
410 1190934 : var->index += d1;
411 : else
412 15226941 : 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 3302835 : 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 3302835 : if (!tab1 || !tab2)
442 0 : return NULL;
443 :
444 3302835 : isl_assert(tab1->mat->ctx, tab1->M == tab2->M, return NULL);
445 3302835 : isl_assert(tab1->mat->ctx, tab1->rational == tab2->rational, return NULL);
446 3302835 : isl_assert(tab1->mat->ctx, tab1->cone == tab2->cone, return NULL);
447 3302835 : isl_assert(tab1->mat->ctx, !tab1->row_sign, return NULL);
448 3302835 : isl_assert(tab1->mat->ctx, !tab2->row_sign, return NULL);
449 3302835 : isl_assert(tab1->mat->ctx, tab1->n_param == 0, return NULL);
450 3302835 : isl_assert(tab1->mat->ctx, tab2->n_param == 0, return NULL);
451 3302835 : isl_assert(tab1->mat->ctx, tab1->n_div == 0, return NULL);
452 3302835 : isl_assert(tab1->mat->ctx, tab2->n_div == 0, return NULL);
453 :
454 3302835 : off = 2 + tab1->M;
455 3302835 : r1 = tab1->n_redundant;
456 3302835 : r2 = tab2->n_redundant;
457 3302835 : d1 = tab1->n_dead;
458 3302835 : d2 = tab2->n_dead;
459 3302835 : prod = isl_calloc_type(tab1->mat->ctx, struct isl_tab);
460 3302835 : if (!prod)
461 0 : return NULL;
462 3302835 : 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 3302835 : if (!prod->mat)
466 0 : goto error;
467 3302835 : prod->var = isl_alloc_array(tab1->mat->ctx, struct isl_tab_var,
468 : tab1->max_var + tab2->max_var);
469 3302835 : if ((tab1->max_var + tab2->max_var) && !prod->var)
470 0 : goto error;
471 20161151 : for (i = 0; i < tab1->n_var; ++i) {
472 16858316 : prod->var[i] = tab1->var[i];
473 16858316 : update_index1(&prod->var[i], r1, r2, d1, d2);
474 : }
475 20161151 : for (i = 0; i < tab2->n_var; ++i) {
476 16858316 : prod->var[tab1->n_var + i] = tab2->var[i];
477 16858316 : update_index2(&prod->var[tab1->n_var + i],
478 : tab1->n_row, tab1->n_col,
479 : r1, r2, d1, d2);
480 : }
481 3302835 : prod->con = isl_alloc_array(tab1->mat->ctx, struct isl_tab_var,
482 : tab1->max_con + tab2->max_con);
483 3302835 : if ((tab1->max_con + tab2->max_con) && !prod->con)
484 0 : goto error;
485 66935760 : for (i = 0; i < tab1->n_con; ++i) {
486 63632925 : prod->con[i] = tab1->con[i];
487 63632925 : update_index1(&prod->con[i], r1, r2, d1, d2);
488 : }
489 66935760 : for (i = 0; i < tab2->n_con; ++i) {
490 63632925 : prod->con[tab1->n_con + i] = tab2->con[i];
491 63632925 : update_index2(&prod->con[tab1->n_con + i],
492 : tab1->n_row, tab1->n_col,
493 : r1, r2, d1, d2);
494 : }
495 3302835 : prod->col_var = isl_alloc_array(tab1->mat->ctx, int,
496 : tab1->n_col + tab2->n_col);
497 3302835 : if ((tab1->n_col + tab2->n_col) && !prod->col_var)
498 0 : goto error;
499 19720710 : for (i = 0; i < tab1->n_col; ++i) {
500 16417875 : int pos = i < d1 ? i : i + d2;
501 16417875 : prod->col_var[pos] = tab1->col_var[i];
502 : }
503 19720710 : for (i = 0; i < tab2->n_col; ++i) {
504 16417875 : int pos = i < d2 ? d1 + i : tab1->n_col + i;
505 16417875 : int t = tab2->col_var[i];
506 16417875 : if (t >= 0)
507 5699 : t += tab1->n_var;
508 : else
509 16412176 : t -= tab1->n_con;
510 16417875 : prod->col_var[pos] = t;
511 : }
512 3302835 : prod->row_var = isl_alloc_array(tab1->mat->ctx, int,
513 : tab1->mat->n_row + tab2->mat->n_row);
514 3302835 : if ((tab1->mat->n_row + tab2->mat->n_row) && !prod->row_var)
515 0 : goto error;
516 66935760 : for (i = 0; i < tab1->n_row; ++i) {
517 63632925 : int pos = i < r1 ? i : i + r2;
518 63632925 : prod->row_var[pos] = tab1->row_var[i];
519 : }
520 66935760 : for (i = 0; i < tab2->n_row; ++i) {
521 63632925 : int pos = i < r2 ? r1 + i : tab1->n_row + i;
522 63632925 : int t = tab2->row_var[i];
523 63632925 : if (t >= 0)
524 16852617 : t += tab1->n_var;
525 : else
526 46780308 : t -= tab1->n_con;
527 63632925 : prod->row_var[pos] = t;
528 : }
529 3302835 : prod->samples = NULL;
530 3302835 : prod->sample_index = NULL;
531 3302835 : prod->n_row = tab1->n_row + tab2->n_row;
532 3302835 : prod->n_con = tab1->n_con + tab2->n_con;
533 3302835 : prod->n_eq = 0;
534 3302835 : prod->max_con = tab1->max_con + tab2->max_con;
535 3302835 : prod->n_col = tab1->n_col + tab2->n_col;
536 3302835 : prod->n_var = tab1->n_var + tab2->n_var;
537 3302835 : prod->max_var = tab1->max_var + tab2->max_var;
538 3302835 : prod->n_param = 0;
539 3302835 : prod->n_div = 0;
540 3302835 : prod->n_dead = tab1->n_dead + tab2->n_dead;
541 3302835 : prod->n_redundant = tab1->n_redundant + tab2->n_redundant;
542 3302835 : prod->rational = tab1->rational;
543 3302835 : prod->empty = tab1->empty || tab2->empty;
544 3302835 : prod->strict_redundant = tab1->strict_redundant || tab2->strict_redundant;
545 3302835 : prod->need_undo = 0;
546 3302835 : prod->in_undo = 0;
547 3302835 : prod->M = tab1->M;
548 3302835 : prod->cone = tab1->cone;
549 3302835 : prod->bottom.type = isl_tab_undo_bottom;
550 3302835 : prod->bottom.next = NULL;
551 3302835 : prod->top = &prod->bottom;
552 :
553 3302835 : prod->n_zero = 0;
554 3302835 : prod->n_unbounded = 0;
555 3302835 : prod->basis = NULL;
556 :
557 3302835 : return prod;
558 : error:
559 0 : isl_tab_free(prod);
560 0 : return NULL;
561 : }
562 :
563 >12866*10^8 : static struct isl_tab_var *var_from_index(struct isl_tab *tab, int i)
564 : {
565 >12866*10^8 : if (i >= 0)
566 >30420*10^7 : return &tab->var[i];
567 : else
568 >98239*10^7 : return &tab->con[~i];
569 : }
570 :
571 >79627*10^7 : struct isl_tab_var *isl_tab_var_from_row(struct isl_tab *tab, int i)
572 : {
573 >79627*10^7 : return var_from_index(tab, tab->row_var[i]);
574 : }
575 :
576 >34128*10^7 : static struct isl_tab_var *var_from_col(struct isl_tab *tab, int i)
577 : {
578 >34128*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 12273619118 : static int max_is_manifestly_unbounded(struct isl_tab *tab,
586 : struct isl_tab_var *var)
587 : {
588 : int i;
589 12273619118 : unsigned off = 2 + tab->M;
590 :
591 12273619118 : if (var->is_row)
592 6497819905 : return 0;
593 35368146986 : for (i = tab->n_redundant; i < tab->n_row; ++i) {
594 31271305798 : if (!isl_int_is_neg(tab->mat->row[i][off + var->index]))
595 23647230016 : continue;
596 7624075782 : if (isl_tab_var_from_row(tab, i)->is_nonneg)
597 1678958025 : return 0;
598 : }
599 4096841188 : 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 15824112525 : static int min_is_manifestly_unbounded(struct isl_tab *tab,
607 : struct isl_tab_var *var)
608 : {
609 : int i;
610 15824112525 : unsigned off = 2 + tab->M;
611 :
612 15824112525 : if (var->is_row)
613 12034938905 : return 0;
614 17566302782 : for (i = tab->n_redundant; i < tab->n_row; ++i) {
615 16858033899 : if (!isl_int_is_pos(tab->mat->row[i][off + var->index]))
616 10868393483 : continue;
617 5989640416 : if (isl_tab_var_from_row(tab, i)->is_nonneg)
618 3080904737 : return 0;
619 : }
620 708268883 : return 1;
621 : }
622 :
623 36934495252 : static int row_cmp(struct isl_tab *tab, int r1, int r2, int c, isl_int *t)
624 : {
625 36934495252 : unsigned off = 2 + tab->M;
626 :
627 36934495252 : 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 36934495252 : isl_int_mul(*t, tab->mat->row[r1][1], tab->mat->row[r2][off + c]);
636 36934495252 : isl_int_submul(*t, tab->mat->row[r2][1], tab->mat->row[r1][off + c]);
637 36934495252 : 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 25556208138 : 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 25556208138 : unsigned off = 2 + tab->M;
666 :
667 25556208138 : isl_int_init(t);
668 25556208138 : r = -1;
669 >25769*10^7 : for (j = tab->n_redundant; j < tab->n_row; ++j) {
670 >23214*10^7 : if (var && j == var->index)
671 22386517434 : continue;
672 >20975*10^7 : if (!isl_tab_var_from_row(tab, j)->is_nonneg)
673 76897339959 : continue;
674 >13285*10^7 : if (sgn * isl_int_sgn(tab->mat->row[j][off + c]) >= 0)
675 77601456905 : continue;
676 55255447064 : if (r < 0) {
677 18320951812 : r = j;
678 18320951812 : continue;
679 : }
680 36934495252 : tsgn = sgn * row_cmp(tab, r, j, c, &t);
681 51789099973 : if (tsgn < 0 || (tsgn == 0 &&
682 14854604721 : tab->row_var[j] < tab->row_var[r]))
683 19213956017 : r = j;
684 : }
685 25556208138 : isl_int_clear(t);
686 25556208138 : 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 57586706830 : 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 57586706830 : *row = *col = -1;
713 :
714 57586706830 : isl_assert(tab->mat->ctx, var->is_row, return);
715 57586706830 : tr = tab->mat->row[var->index] + 2 + tab->M;
716 :
717 57586706830 : c = -1;
718 >22710*10^7 : for (j = tab->n_dead; j < tab->n_col; ++j) {
719 >16951*10^7 : if (isl_int_is_zero(tr[j]))
720 69578793417 : continue;
721 >17343*10^7 : if (isl_int_sgn(tr[j]) != sgn &&
722 73495490935 : var_from_col(tab, j)->is_nonneg)
723 60024995852 : continue;
724 39912154634 : if (c < 0 || tab->col_var[j] < tab->col_var[c])
725 25448939684 : c = j;
726 : }
727 57586706830 : if (c < 0)
728 35040166089 : return;
729 :
730 22546540741 : sgn *= isl_int_sgn(tr[c]);
731 22546540741 : r = pivot_row(tab, skip_var, sgn, c);
732 22546540741 : *row = r < 0 ? var->index : r;
733 22546540741 : *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 >25076*10^7 : int isl_tab_row_is_redundant(struct isl_tab *tab, int row)
743 : {
744 : int i;
745 >25076*10^7 : unsigned off = 2 + tab->M;
746 :
747 >25076*10^7 : if (tab->row_var[row] < 0 && !isl_tab_var_from_row(tab, row)->is_nonneg)
748 20322982079 : return 0;
749 :
750 >23044*10^7 : if (isl_int_is_neg(tab->mat->row[row][1]))
751 37738236787 : return 0;
752 >19270*10^7 : if (tab->strict_redundant && isl_int_is_zero(tab->mat->row[row][1]))
753 0 : return 0;
754 >19270*10^7 : if (tab->M && isl_int_is_neg(tab->mat->row[row][2]))
755 0 : return 0;
756 :
757 >55196*10^7 : for (i = tab->n_dead; i < tab->n_col; ++i) {
758 >53026*10^7 : if (isl_int_is_zero(tab->mat->row[row][off + i]))
759 >24176*10^7 : continue;
760 >28849*10^7 : if (tab->col_var[i] >= 0)
761 61370961373 : return 0;
762 >22712*10^7 : if (isl_int_is_neg(tab->mat->row[row][off + i]))
763 94816570947 : return 0;
764 >13231*10^7 : if (!var_from_col(tab, i)->is_nonneg)
765 14812386638 : return 0;
766 : }
767 21705458331 : return 1;
768 : }
769 :
770 30472455619 : 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 30472455619 : t = tab->row_var[row1];
776 30472455619 : tab->row_var[row1] = tab->row_var[row2];
777 30472455619 : tab->row_var[row2] = t;
778 30472455619 : isl_tab_var_from_row(tab, row1)->index = row1;
779 30472455619 : isl_tab_var_from_row(tab, row2)->index = row2;
780 30472455619 : tab->mat = isl_mat_swap_rows(tab->mat, row1, row2);
781 :
782 30472455619 : if (!tab->row_sign)
783 30472455619 : 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 >26352*10^7 : 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 >26352*10^7 : if (!tab)
805 0 : return isl_stat_error;
806 >26352*10^7 : if (!tab->need_undo)
807 85808834205 : return isl_stat_ok;
808 :
809 >17771*10^7 : undo = isl_alloc_type(tab->mat->ctx, struct isl_tab_undo);
810 >17771*10^7 : if (!undo)
811 0 : goto error;
812 >17771*10^7 : undo->type = type;
813 >17771*10^7 : undo->u = u;
814 >17771*10^7 : undo->next = tab->top;
815 >17771*10^7 : tab->top = undo;
816 :
817 >17771*10^7 : 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 >23595*10^7 : 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 >23595*10^7 : if (var->is_row)
829 >19900*10^7 : u.var_index = tab->row_var[var->index];
830 : else
831 36955139141 : u.var_index = tab->col_var[var->index];
832 >23595*10^7 : return push_union(tab, type, u);
833 : }
834 :
835 27568555880 : isl_stat isl_tab_push(struct isl_tab *tab, enum isl_tab_undo_type type)
836 : {
837 27568555880 : union isl_tab_undo_val u = { 0 };
838 27568555880 : 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 27048514827 : int isl_tab_mark_redundant(struct isl_tab *tab, int row)
957 : {
958 27048514827 : struct isl_tab_var *var = isl_tab_var_from_row(tab, row);
959 27048514827 : var->is_redundant = 1;
960 27048514827 : isl_assert(tab->mat->ctx, row >= tab->n_redundant, return -1);
961 27048514827 : if (tab->preserve || tab->need_undo || tab->row_var[row] >= 0) {
962 19010150760 : if (tab->row_var[row] >= 0 && !var->is_nonneg) {
963 9246361619 : var->is_nonneg = 1;
964 9246361619 : if (isl_tab_push_var(tab, isl_tab_undo_nonneg, var) < 0)
965 0 : return -1;
966 : }
967 19010150760 : if (row != tab->n_redundant)
968 16319986643 : swap_rows(tab, row, tab->n_redundant);
969 19010150760 : tab->n_redundant++;
970 19010150760 : return isl_tab_push_var(tab, isl_tab_undo_redundant, var);
971 : } else {
972 8038364067 : if (row != tab->n_row - 1)
973 5447318451 : swap_rows(tab, row, tab->n_row - 1);
974 8038364067 : isl_tab_var_from_row(tab, tab->n_row - 1)->index = -1;
975 8038364067 : tab->n_row--;
976 8038364067 : 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 146948 : int isl_tab_mark_rational(struct isl_tab *tab)
986 : {
987 146948 : if (!tab)
988 0 : return -1;
989 146948 : if (!tab->rational && tab->need_undo)
990 146436 : if (isl_tab_push(tab, isl_tab_undo_rational) < 0)
991 0 : return -1;
992 146948 : tab->rational = 1;
993 146948 : return 0;
994 : }
995 :
996 8741078486 : isl_stat isl_tab_mark_empty(struct isl_tab *tab)
997 : {
998 8741078486 : if (!tab)
999 0 : return isl_stat_error;
1000 8741078486 : if (!tab->empty && tab->need_undo)
1001 6822784797 : if (isl_tab_push(tab, isl_tab_undo_empty) < 0)
1002 0 : return isl_stat_error;
1003 8741078486 : tab->empty = 1;
1004 8741078486 : return isl_stat_ok;
1005 : }
1006 :
1007 71753564463 : int isl_tab_freeze_constraint(struct isl_tab *tab, int con)
1008 : {
1009 : struct isl_tab_var *var;
1010 :
1011 71753564463 : if (!tab)
1012 0 : return -1;
1013 :
1014 71753564463 : var = &tab->con[con];
1015 71753564463 : if (var->frozen)
1016 0 : return 0;
1017 71753564463 : if (var->index < 0)
1018 1483569404 : return 0;
1019 70269995059 : var->frozen = 1;
1020 :
1021 70269995059 : if (tab->need_undo)
1022 70269995059 : 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 44493432228 : static void update_row_sign(struct isl_tab *tab, int row, int col, int row_sgn)
1045 : {
1046 : int i;
1047 44493432228 : struct isl_mat *mat = tab->mat;
1048 44493432228 : unsigned off = 2 + tab->M;
1049 :
1050 44493432228 : if (!tab->row_sign)
1051 44493432228 : 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 44493432228 : 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 44493432228 : struct isl_mat *mat = tab->mat;
1132 : struct isl_tab_var *var;
1133 44493432228 : unsigned off = 2 + tab->M;
1134 :
1135 44493432228 : ctx = isl_tab_get_ctx(tab);
1136 44493432228 : if (isl_ctx_next_operation(ctx) < 0)
1137 0 : return -1;
1138 :
1139 44493432228 : isl_int_swap(mat->row[row][0], mat->row[row][off + col]);
1140 44493432228 : sgn = isl_int_sgn(mat->row[row][0]);
1141 44493432228 : if (sgn < 0) {
1142 16405290621 : isl_int_neg(mat->row[row][0], mat->row[row][0]);
1143 16405290621 : isl_int_neg(mat->row[row][off + col], mat->row[row][off + col]);
1144 : } else
1145 >21273*10^7 : for (j = 0; j < off - 1 + tab->n_col; ++j) {
1146 >18465*10^7 : if (j == off - 1 + col)
1147 28088141607 : continue;
1148 >15656*10^7 : isl_int_neg(mat->row[row][1 + j], mat->row[row][1 + j]);
1149 : }
1150 44493432228 : if (!isl_int_is_one(mat->row[row][0]))
1151 13443961054 : isl_seq_normalize(mat->ctx, mat->row[row], off + tab->n_col);
1152 >41406*10^7 : for (i = 0; i < tab->n_row; ++i) {
1153 >36957*10^7 : if (i == row)
1154 44493432228 : continue;
1155 >32507*10^7 : if (isl_int_is_zero(mat->row[i][off + col]))
1156 >15855*10^7 : continue;
1157 >16652*10^7 : isl_int_mul(mat->row[i][0], mat->row[i][0], mat->row[row][0]);
1158 >12779*10^8 : for (j = 0; j < off - 1 + tab->n_col; ++j) {
1159 >11113*10^8 : if (j == off - 1 + col)
1160 >16652*10^7 : continue;
1161 >94487*10^7 : isl_int_mul(mat->row[i][1 + j],
1162 : mat->row[i][1 + j], mat->row[row][0]);
1163 >94487*10^7 : isl_int_addmul(mat->row[i][1 + j],
1164 : mat->row[i][off + col], mat->row[row][1 + j]);
1165 : }
1166 >16652*10^7 : isl_int_mul(mat->row[i][off + col],
1167 : mat->row[i][off + col], mat->row[row][off + col]);
1168 >16652*10^7 : if (!isl_int_is_one(mat->row[i][0]))
1169 >11986*10^7 : isl_seq_normalize(mat->ctx, mat->row[i], off + tab->n_col);
1170 : }
1171 44493432228 : t = tab->row_var[row];
1172 44493432228 : tab->row_var[row] = tab->col_var[col];
1173 44493432228 : tab->col_var[col] = t;
1174 44493432228 : var = isl_tab_var_from_row(tab, row);
1175 44493432228 : var->is_row = 1;
1176 44493432228 : var->index = row;
1177 44493432228 : var = var_from_col(tab, col);
1178 44493432228 : var->is_row = 0;
1179 44493432228 : var->index = col;
1180 44493432228 : update_row_sign(tab, row, col, sgn);
1181 44493432228 : if (tab->in_undo)
1182 2042571711 : return 0;
1183 >35097*10^7 : for (i = tab->n_redundant; i < tab->n_row; ++i) {
1184 >30852*10^7 : if (isl_int_is_zero(mat->row[i][off + col]))
1185 >13027*10^7 : continue;
1186 >33580*10^7 : if (!isl_tab_var_from_row(tab, i)->frozen &&
1187 >15756*10^7 : isl_tab_row_is_redundant(tab, i)) {
1188 12971790352 : int redo = isl_tab_mark_redundant(tab, i);
1189 12971790352 : if (redo < 0)
1190 0 : return -1;
1191 12971790352 : if (redo)
1192 1837332280 : --i;
1193 : }
1194 : }
1195 42450860517 : 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 9218829984 : static int to_row(struct isl_tab *tab, struct isl_tab_var *var, int sign)
1206 : {
1207 : int r;
1208 9218829984 : unsigned off = 2 + tab->M;
1209 :
1210 9218829984 : if (var->is_row)
1211 6503028973 : return 0;
1212 :
1213 2715801011 : if (sign == 0) {
1214 764467271 : for (r = tab->n_redundant; r < tab->n_row; ++r)
1215 764467271 : if (!isl_int_is_zero(tab->mat->row[r][off+var->index]))
1216 : break;
1217 237733649 : isl_assert(tab->mat->ctx, r < tab->n_row, return -1);
1218 : } else {
1219 2478067362 : r = pivot_row(tab, NULL, sign, var->index);
1220 2478067362 : isl_assert(tab->mat->ctx, r >= 0, return -1);
1221 : }
1222 :
1223 2715801011 : 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 7799083033 : static int sign_of_max(struct isl_tab *tab, struct isl_tab_var *var)
1265 : {
1266 : int row, col;
1267 :
1268 7799083033 : if (max_is_manifestly_unbounded(tab, var))
1269 859444004 : return 1;
1270 6939639029 : if (to_row(tab, var, 1) < 0)
1271 0 : return -2;
1272 20620672118 : while (!isl_int_is_pos(tab->mat->row[var->index][1])) {
1273 13408018163 : find_pivot(tab, var, var, 1, &row, &col);
1274 13408018163 : if (row == -1)
1275 5326264502 : return isl_int_sgn(tab->mat->row[var->index][1]);
1276 8081753661 : if (isl_tab_pivot(tab, row, col) < 0)
1277 0 : return -2;
1278 8081753661 : if (!var->is_row) /* manifestly unbounded */
1279 1340359601 : return 1;
1280 : }
1281 273014926 : 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 59027264137 : static int row_is_neg(struct isl_tab *tab, int row)
1299 : {
1300 59027264137 : if (!tab->M)
1301 59027264137 : 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 48122372165 : static int row_sgn(struct isl_tab *tab, int row)
1310 : {
1311 48122372165 : if (!tab->M)
1312 48122372165 : 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 53645324865 : static int restore_row(struct isl_tab *tab, struct isl_tab_var *var)
1325 : {
1326 : int row, col;
1327 :
1328 >11267*10^7 : while (row_is_neg(tab, var->index)) {
1329 19649467553 : find_pivot(tab, var, var, 1, &row, &col);
1330 19649467553 : if (row == -1)
1331 8744575581 : break;
1332 10904891972 : if (isl_tab_pivot(tab, row, col) < 0)
1333 0 : return -2;
1334 10904891972 : if (!var->is_row) /* manifestly unbounded */
1335 5522952700 : return 1;
1336 : }
1337 48122372165 : 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 7122673533 : static int at_least_zero(struct isl_tab *tab, struct isl_tab_var *var)
1346 : {
1347 : int row, col;
1348 :
1349 14271086468 : while (isl_int_is_neg(tab->mat->row[var->index][1])) {
1350 7138958804 : find_pivot(tab, var, var, 1, &row, &col);
1351 7138958804 : if (row == -1)
1352 7099033178 : break;
1353 39925626 : if (row == var->index) /* manifestly unbounded */
1354 14186224 : return 1;
1355 25739402 : if (isl_tab_pivot(tab, row, col) < 0)
1356 0 : return -1;
1357 : }
1358 7108487309 : 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 11221856972 : static int sign_of_min(struct isl_tab *tab, struct isl_tab_var *var)
1379 : {
1380 : int row, col;
1381 11221856972 : struct isl_tab_var *pivot_var = NULL;
1382 :
1383 11221856972 : if (min_is_manifestly_unbounded(tab, var))
1384 0 : return -1;
1385 11221856972 : if (!var->is_row) {
1386 900220 : col = var->index;
1387 900220 : row = pivot_row(tab, NULL, -1, col);
1388 900220 : pivot_var = var_from_col(tab, col);
1389 900220 : if (isl_tab_pivot(tab, row, col) < 0)
1390 0 : return -2;
1391 900220 : if (var->is_redundant)
1392 144351 : return 0;
1393 755869 : if (isl_int_is_neg(tab->mat->row[var->index][1])) {
1394 449887 : if (var->is_nonneg) {
1395 899774 : if (!pivot_var->is_redundant &&
1396 449887 : pivot_var->index == row) {
1397 838602 : if (isl_tab_pivot(tab, row, col) < 0)
1398 0 : return -2;
1399 : } else
1400 30586 : if (restore_row(tab, var) < -1)
1401 0 : return -2;
1402 : }
1403 449887 : return -1;
1404 : }
1405 : }
1406 11221262734 : if (var->is_redundant)
1407 0 : return 0;
1408 22467094411 : while (!isl_int_is_neg(tab->mat->row[var->index][1])) {
1409 11232056486 : find_pivot(tab, var, var, -1, &row, &col);
1410 11232056486 : if (row == var->index)
1411 25780111 : return -1;
1412 11206276375 : if (row == -1)
1413 11181226503 : return isl_int_sgn(tab->mat->row[var->index][1]);
1414 25049872 : pivot_var = var_from_col(tab, col);
1415 25049872 : if (isl_tab_pivot(tab, row, col) < 0)
1416 0 : return -2;
1417 25049872 : if (var->is_redundant)
1418 480929 : return 0;
1419 : }
1420 13775191 : if (pivot_var && var->is_nonneg) {
1421 : /* pivot back to non-negative value */
1422 512269 : if (!pivot_var->is_redundant && pivot_var->index == row) {
1423 973728 : if (isl_tab_pivot(tab, row, col) < 0)
1424 0 : return -2;
1425 : } else
1426 25405 : if (restore_row(tab, var) < -1)
1427 0 : return -2;
1428 : }
1429 13775191 : return -1;
1430 : }
1431 :
1432 1424913571 : static int row_at_most_neg_one(struct isl_tab *tab, int row)
1433 : {
1434 1424913571 : 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 2044876997 : return isl_int_is_neg(tab->mat->row[row][1]) &&
1441 619963426 : 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 1339472900 : 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 1339472900 : if (min_is_manifestly_unbounded(tab, var))
1460 0 : return 1;
1461 1339472900 : if (!var->is_row) {
1462 530699815 : col = var->index;
1463 530699815 : row = pivot_row(tab, NULL, -1, col);
1464 530699815 : pivot_var = var_from_col(tab, col);
1465 530699815 : if (isl_tab_pivot(tab, row, col) < 0)
1466 0 : return -1;
1467 530699815 : if (var->is_redundant)
1468 61346054 : return 0;
1469 469353761 : if (row_at_most_neg_one(tab, var->index)) {
1470 300696466 : if (var->is_nonneg) {
1471 601392932 : if (!pivot_var->is_redundant &&
1472 300696466 : pivot_var->index == row) {
1473 596397464 : if (isl_tab_pivot(tab, row, col) < 0)
1474 0 : return -1;
1475 : } else
1476 2497734 : if (restore_row(tab, var) < -1)
1477 0 : return -1;
1478 : }
1479 300696466 : return 1;
1480 : }
1481 : }
1482 977430380 : if (var->is_redundant)
1483 1 : return 0;
1484 : do {
1485 1714456894 : find_pivot(tab, var, var, -1, &row, &col);
1486 1714456894 : if (row == var->index) {
1487 304892108 : if (var->is_nonneg && restore_row(tab, var) < -1)
1488 0 : return -1;
1489 304892108 : return 1;
1490 : }
1491 1409564786 : if (row == -1)
1492 168617047 : return 0;
1493 1240947739 : pivot_var = var_from_col(tab, col);
1494 1240947739 : if (isl_tab_pivot(tab, row, col) < 0)
1495 0 : return -1;
1496 1240947739 : if (var->is_redundant)
1497 285387929 : return 0;
1498 955559810 : } while (!row_at_most_neg_one(tab, var->index));
1499 218533295 : if (var->is_nonneg) {
1500 : /* pivot back to non-negative value */
1501 211065489 : if (!pivot_var->is_redundant && pivot_var->index == row)
1502 209392545 : if (isl_tab_pivot(tab, row, col) < 0)
1503 0 : return -1;
1504 211065489 : if (restore_row(tab, var) < -1)
1505 0 : return -1;
1506 : }
1507 218533295 : return 1;
1508 : }
1509 :
1510 : /* Return 1 if "var" can attain values >= 1.
1511 : * Return 0 otherwise.
1512 : */
1513 2436146997 : 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 2436146997 : if (max_is_manifestly_unbounded(tab, var))
1519 2199803605 : return 1;
1520 236343392 : if (to_row(tab, var, 1) < 0)
1521 0 : return -1;
1522 236343392 : r = tab->mat->row[var->index];
1523 489974745 : while (isl_int_lt(r[1], r[0])) {
1524 43305818 : find_pivot(tab, var, var, 1, &row, &col);
1525 43305818 : if (row == -1)
1526 2835708 : return isl_int_ge(r[1], r[0]);
1527 40470110 : if (row == var->index) /* manifestly unbounded */
1528 23182149 : return 1;
1529 17287961 : if (isl_tab_pivot(tab, row, col) < 0)
1530 0 : return -1;
1531 : }
1532 210325535 : return 1;
1533 : }
1534 :
1535 12568376930 : static void swap_cols(struct isl_tab *tab, int col1, int col2)
1536 : {
1537 : int t;
1538 12568376930 : unsigned off = 2 + tab->M;
1539 12568376930 : t = tab->col_var[col1];
1540 12568376930 : tab->col_var[col1] = tab->col_var[col2];
1541 12568376930 : tab->col_var[col2] = t;
1542 12568376930 : var_from_col(tab, col1)->index = col1;
1543 12568376930 : var_from_col(tab, col2)->index = col2;
1544 12568376930 : tab->mat = isl_mat_swap_cols(tab->mat, off + col1, off + col2);
1545 12568376930 : }
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 23843744497 : int isl_tab_kill_col(struct isl_tab *tab, int col)
1560 : {
1561 23843744497 : var_from_col(tab, col)->is_zero = 1;
1562 23843744497 : if (tab->need_undo) {
1563 101549971 : if (isl_tab_push_var(tab, isl_tab_undo_zero,
1564 : var_from_col(tab, col)) < 0)
1565 0 : return -1;
1566 101549971 : if (col != tab->n_dead)
1567 64087872 : swap_cols(tab, col, tab->n_dead);
1568 101549971 : tab->n_dead++;
1569 101549971 : return 0;
1570 : } else {
1571 23742194526 : if (col != tab->n_col - 1)
1572 12504289058 : swap_cols(tab, col, tab->n_col - 1);
1573 23742194526 : var_from_col(tab, tab->n_col - 1)->index = -1;
1574 23742194526 : tab->n_col--;
1575 23742194526 : return 1;
1576 : }
1577 : }
1578 :
1579 23426833698 : static int row_is_manifestly_non_integral(struct isl_tab *tab, int row)
1580 : {
1581 23426833698 : unsigned off = 2 + tab->M;
1582 :
1583 23426833698 : if (tab->M && !isl_int_eq(tab->mat->row[row][2],
1584 : tab->mat->row[row][0]))
1585 0 : return 0;
1586 23426833698 : if (isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
1587 23426833698 : tab->n_col - tab->n_dead) != -1)
1588 3385309585 : return 0;
1589 :
1590 20041524113 : 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 5326213219 : static int tab_is_manifestly_empty(struct isl_tab *tab)
1598 : {
1599 : int i;
1600 :
1601 5326213219 : if (tab->empty)
1602 3255 : return 1;
1603 5326209964 : if (tab->rational)
1604 1147599 : return 0;
1605 :
1606 31554060641 : for (i = 0; i < tab->n_var; ++i) {
1607 26229080103 : if (!tab->var[i].is_row)
1608 2802246405 : continue;
1609 23426833698 : if (row_is_manifestly_non_integral(tab, tab->var[i].index))
1610 81827 : return 1;
1611 : }
1612 :
1613 5324980538 : 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 5326213219 : static isl_stat close_row(struct isl_tab *tab, struct isl_tab_var *var,
1636 : int temp_var)
1637 : {
1638 : int j;
1639 5326213219 : struct isl_mat *mat = tab->mat;
1640 5326213219 : unsigned off = 2 + tab->M;
1641 :
1642 5326213219 : 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 5326213219 : var->is_zero = 1;
1647 5326213219 : 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 15314659418 : for (j = tab->n_dead; j < tab->n_col; ++j) {
1651 : int recheck;
1652 9988446199 : if (isl_int_is_zero(mat->row[var->index][off + j]))
1653 4798749095 : continue;
1654 5189697104 : 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 5189697104 : recheck = isl_tab_kill_col(tab, j);
1659 5189697104 : if (recheck < 0)
1660 0 : return isl_stat_error;
1661 5189697104 : if (recheck)
1662 5189003373 : --j;
1663 : }
1664 5326213219 : if (!temp_var && isl_tab_mark_redundant(tab, var->index) < 0)
1665 0 : return isl_stat_error;
1666 5326213219 : if (tab_is_manifestly_empty(tab) && isl_tab_mark_empty(tab) < 0)
1667 0 : return isl_stat_error;
1668 5326213219 : 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 88495561809 : int isl_tab_allocate_con(struct isl_tab *tab)
1678 : {
1679 : int r;
1680 :
1681 88495561809 : isl_assert(tab->mat->ctx, tab->n_row < tab->mat->n_row, return -1);
1682 88495561809 : isl_assert(tab->mat->ctx, tab->n_con < tab->max_con, return -1);
1683 :
1684 88495561809 : r = tab->n_con;
1685 88495561809 : tab->con[r].index = tab->n_row;
1686 88495561809 : tab->con[r].is_row = 1;
1687 88495561809 : tab->con[r].is_nonneg = 0;
1688 88495561809 : tab->con[r].is_zero = 0;
1689 88495561809 : tab->con[r].is_redundant = 0;
1690 88495561809 : tab->con[r].frozen = 0;
1691 88495561809 : tab->con[r].negated = 0;
1692 88495561809 : tab->row_var[tab->n_row] = ~r;
1693 :
1694 88495561809 : tab->n_row++;
1695 88495561809 : tab->n_con++;
1696 88495561809 : if (isl_tab_push_var(tab, isl_tab_undo_allocate, &tab->con[r]) < 0)
1697 0 : return -1;
1698 :
1699 88495561809 : 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 88495508940 : 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 88495508940 : unsigned off = 2 + tab->M;
1830 :
1831 88495508940 : r = isl_tab_allocate_con(tab);
1832 88495508940 : if (r < 0)
1833 0 : return -1;
1834 :
1835 88495508940 : isl_int_init(a);
1836 88495508940 : isl_int_init(b);
1837 88495508940 : row = tab->mat->row[tab->con[r].index];
1838 88495508940 : isl_int_set_si(row[0], 1);
1839 88495508940 : isl_int_set(row[1], line[0]);
1840 88495508940 : isl_seq_clr(row + 2, tab->M + tab->n_col);
1841 >64781*10^7 : for (i = 0; i < tab->n_var; ++i) {
1842 >55931*10^7 : if (tab->var[i].is_zero)
1843 0 : continue;
1844 >55931*10^7 : if (tab->var[i].is_row) {
1845 >31607*10^7 : isl_int_lcm(a,
1846 : row[0], tab->mat->row[tab->var[i].index][0]);
1847 >31607*10^7 : isl_int_swap(a, row[0]);
1848 >31607*10^7 : isl_int_divexact(a, row[0], a);
1849 >31607*10^7 : isl_int_divexact(b,
1850 : row[0], tab->mat->row[tab->var[i].index][0]);
1851 >31607*10^7 : isl_int_mul(b, b, line[1 + i]);
1852 >63215*10^7 : isl_seq_combine(row + 1, a, row + 1,
1853 >31607*10^7 : b, tab->mat->row[tab->var[i].index] + 1,
1854 >31607*10^7 : 1 + tab->M + tab->n_col);
1855 : } else
1856 >24324*10^7 : isl_int_addmul(row[off + tab->var[i].index],
1857 : line[1 + i], row[0]);
1858 >55931*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 88495508940 : isl_seq_normalize(tab->mat->ctx, row, off + tab->n_col);
1862 88495508940 : isl_int_clear(a);
1863 88495508940 : isl_int_clear(b);
1864 :
1865 88495508940 : if (tab->row_sign)
1866 0 : tab->row_sign[tab->con[r].index] = isl_tab_row_unknown;
1867 :
1868 88495508940 : return r;
1869 : }
1870 :
1871 41579296890 : static isl_stat drop_row(struct isl_tab *tab, int row)
1872 : {
1873 41579296890 : isl_assert(tab->mat->ctx, ~tab->row_var[row] == tab->n_con - 1,
1874 : return isl_stat_error);
1875 41579296890 : if (row != tab->n_row - 1)
1876 8705150525 : swap_rows(tab, row, tab->n_row - 1);
1877 41579296890 : tab->n_row--;
1878 41579296890 : tab->n_con--;
1879 41579296890 : 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 48831805934 : 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 48831805934 : if (!tab)
1914 0 : return isl_stat_error;
1915 48831805934 : if (tab->bmap) {
1916 20745518909 : struct isl_basic_map *bmap = tab->bmap;
1917 :
1918 20745518909 : isl_assert(tab->mat->ctx, tab->n_eq == bmap->n_eq,
1919 : return isl_stat_error);
1920 20745518909 : isl_assert(tab->mat->ctx,
1921 : tab->n_con == bmap->n_eq + bmap->n_ineq,
1922 : return isl_stat_error);
1923 20745518909 : tab->bmap = isl_basic_map_add_ineq(tab->bmap, ineq);
1924 20745518909 : if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1925 0 : return isl_stat_error;
1926 20745518909 : if (!tab->bmap)
1927 0 : return isl_stat_error;
1928 : }
1929 48831805934 : 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 48831805934 : r = isl_tab_add_row(tab, ineq);
1935 48831805934 : if (tab->cone) {
1936 0 : isl_int_swap(ineq[0], cst);
1937 0 : isl_int_clear(cst);
1938 : }
1939 48831805934 : if (r < 0)
1940 0 : return isl_stat_error;
1941 48831805934 : tab->con[r].is_nonneg = 1;
1942 48831805934 : if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1943 0 : return isl_stat_error;
1944 48831805934 : if (isl_tab_row_is_redundant(tab, tab->con[r].index)) {
1945 8733667979 : if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0)
1946 0 : return isl_stat_error;
1947 8733667979 : return isl_stat_ok;
1948 : }
1949 :
1950 40098137955 : sgn = restore_row(tab, &tab->con[r]);
1951 40098137955 : if (sgn < -1)
1952 0 : return isl_stat_error;
1953 40098137955 : if (sgn < 0)
1954 8740948149 : return isl_tab_mark_empty(tab);
1955 31357189806 : 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 31357189806 : 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 100856240 : static int to_col(struct isl_tab *tab, struct isl_tab_var *var)
1966 : {
1967 : int i;
1968 : int row, col;
1969 100856240 : unsigned off = 2 + tab->M;
1970 :
1971 100856240 : if (!var->is_row)
1972 0 : return 0;
1973 :
1974 268561453 : while (isl_int_is_pos(tab->mat->row[var->index][1])) {
1975 157959342 : find_pivot(tab, var, NULL, -1, &row, &col);
1976 157959342 : isl_assert(tab->mat->ctx, row != -1, return -1);
1977 157959342 : if (isl_tab_pivot(tab, row, col) < 0)
1978 0 : return -1;
1979 157959342 : if (!var->is_row)
1980 91110369 : return 0;
1981 : }
1982 :
1983 16498861 : for (i = tab->n_dead; i < tab->n_col; ++i)
1984 16498861 : if (!isl_int_is_zero(tab->mat->row[var->index][off + i]))
1985 : break;
1986 :
1987 9745871 : isl_assert(tab->mat->ctx, i < tab->n_col, return -1);
1988 9745871 : if (isl_tab_pivot(tab, var->index, i) < 0)
1989 0 : return -1;
1990 :
1991 9745871 : 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 18553191153 : static struct isl_tab *add_eq(struct isl_tab *tab, isl_int *eq)
2003 : {
2004 : int i;
2005 : int r;
2006 :
2007 18553191153 : if (!tab)
2008 0 : return NULL;
2009 18553191153 : r = isl_tab_add_row(tab, eq);
2010 18553191153 : if (r < 0)
2011 0 : goto error;
2012 :
2013 18553191153 : r = tab->con[r].index;
2014 18553191153 : i = isl_seq_first_non_zero(tab->mat->row[r] + 2 + tab->M + tab->n_dead,
2015 18553191153 : tab->n_col - tab->n_dead);
2016 18553191153 : isl_assert(tab->mat->ctx, i >= 0, goto error);
2017 18553191153 : i += tab->n_dead;
2018 18553191153 : if (isl_tab_pivot(tab, r, i) < 0)
2019 0 : goto error;
2020 18553191153 : if (isl_tab_kill_col(tab, i) < 0)
2021 0 : goto error;
2022 18553191153 : tab->n_eq++;
2023 :
2024 18553191153 : 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 392851442 : static int row_is_big(struct isl_tab *tab, int row)
2034 : {
2035 392851442 : return tab->M && !isl_int_is_zero(tab->mat->row[row][2]);
2036 : }
2037 :
2038 102174081 : static int row_is_manifestly_zero(struct isl_tab *tab, int row)
2039 : {
2040 102174081 : unsigned off = 2 + tab->M;
2041 :
2042 102174081 : if (!isl_int_is_zero(tab->mat->row[row][1]))
2043 91116397 : return 0;
2044 11057684 : if (row_is_big(tab, row))
2045 0 : return 0;
2046 22115368 : return isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
2047 22115368 : 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 102110736 : 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 102110736 : if (!tab)
2061 0 : return -1;
2062 102110736 : r = isl_tab_add_row(tab, eq);
2063 102110736 : if (r < 0)
2064 0 : return -1;
2065 :
2066 102110736 : var = &tab->con[r];
2067 102110736 : r = var->index;
2068 102110736 : if (row_is_manifestly_zero(tab, r)) {
2069 1313876 : var->is_zero = 1;
2070 1313876 : if (isl_tab_mark_redundant(tab, r) < 0)
2071 0 : return -1;
2072 1313876 : return 0;
2073 : }
2074 :
2075 100796860 : if (isl_int_is_neg(tab->mat->row[r][1])) {
2076 28634946 : isl_seq_neg(tab->mat->row[r] + 1, tab->mat->row[r] + 1,
2077 28634946 : 1 + tab->n_col);
2078 28634946 : var->negated = 1;
2079 : }
2080 100796860 : var->is_nonneg = 1;
2081 100796860 : if (to_col(tab, var) < 0)
2082 0 : return -1;
2083 100796860 : var->is_nonneg = 0;
2084 100796860 : if (isl_tab_kill_col(tab, var->index) < 0)
2085 0 : return -1;
2086 :
2087 100796860 : 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 52869 : static int add_zero_row(struct isl_tab *tab)
2097 : {
2098 : int r;
2099 : isl_int *row;
2100 :
2101 52869 : r = isl_tab_allocate_con(tab);
2102 52869 : if (r < 0)
2103 0 : return -1;
2104 :
2105 52869 : row = tab->mat->row[tab->con[r].index];
2106 52869 : isl_seq_clr(row + 1, 1 + tab->M + tab->n_col);
2107 52869 : isl_int_set_si(row[0], 1);
2108 :
2109 52869 : 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 63345 : int isl_tab_add_eq(struct isl_tab *tab, isl_int *eq)
2120 : {
2121 63345 : 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 63345 : if (!tab)
2129 0 : return -1;
2130 63345 : isl_assert(tab->mat->ctx, !tab->M, return -1);
2131 :
2132 63345 : if (tab->need_undo)
2133 63345 : snap = isl_tab_snap(tab);
2134 :
2135 63345 : 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 63345 : r = isl_tab_add_row(tab, eq);
2141 63345 : if (tab->cone) {
2142 0 : isl_int_swap(eq[0], cst);
2143 0 : isl_int_clear(cst);
2144 : }
2145 63345 : if (r < 0)
2146 0 : return -1;
2147 :
2148 63345 : var = &tab->con[r];
2149 63345 : row = var->index;
2150 63345 : if (row_is_manifestly_zero(tab, row)) {
2151 3965 : if (snap)
2152 3965 : return isl_tab_rollback(tab, snap);
2153 0 : return drop_row(tab, row);
2154 : }
2155 :
2156 59380 : if (tab->bmap) {
2157 52869 : tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);
2158 52869 : if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
2159 0 : return -1;
2160 52869 : isl_seq_neg(eq, eq, 1 + tab->n_var);
2161 52869 : tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);
2162 52869 : isl_seq_neg(eq, eq, 1 + tab->n_var);
2163 52869 : if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
2164 0 : return -1;
2165 52869 : if (!tab->bmap)
2166 0 : return -1;
2167 52869 : if (add_zero_row(tab) < 0)
2168 0 : return -1;
2169 : }
2170 :
2171 59380 : sgn = isl_int_sgn(tab->mat->row[row][1]);
2172 :
2173 59380 : if (sgn > 0) {
2174 1163 : isl_seq_neg(tab->mat->row[row] + 1, tab->mat->row[row] + 1,
2175 1163 : 1 + tab->n_col);
2176 1163 : var->negated = 1;
2177 1163 : sgn = -1;
2178 : }
2179 :
2180 59380 : if (sgn < 0) {
2181 16814 : sgn = sign_of_max(tab, var);
2182 16814 : if (sgn < -1)
2183 0 : return -1;
2184 16814 : if (sgn < 0) {
2185 0 : if (isl_tab_mark_empty(tab) < 0)
2186 0 : return -1;
2187 0 : return 0;
2188 : }
2189 : }
2190 :
2191 59380 : var->is_nonneg = 1;
2192 59380 : if (to_col(tab, var) < 0)
2193 0 : return -1;
2194 59380 : var->is_nonneg = 0;
2195 59380 : if (isl_tab_kill_col(tab, var->index) < 0)
2196 0 : return -1;
2197 :
2198 59380 : 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 8316155043 : __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 8316155043 : if (!bmap)
2397 0 : return NULL;
2398 16632310086 : tab = isl_tab_alloc(bmap->ctx,
2399 8316155043 : isl_basic_map_total_dim(bmap) + bmap->n_ineq + 1,
2400 : isl_basic_map_total_dim(bmap), 0);
2401 8316155043 : if (!tab)
2402 0 : return NULL;
2403 8316155043 : tab->preserve = track;
2404 8316155043 : tab->rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
2405 8316155043 : 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 26864963716 : for (i = 0; i < bmap->n_eq; ++i) {
2411 18548808673 : tab = add_eq(tab, bmap->eq[i]);
2412 18548808673 : if (!tab)
2413 0 : return tab;
2414 : }
2415 35071274743 : for (i = 0; i < bmap->n_ineq; ++i) {
2416 28080554861 : if (isl_tab_add_ineq(tab, bmap->ineq[i]) < 0)
2417 0 : goto error;
2418 28080554861 : if (tab->empty)
2419 1325435161 : goto done;
2420 : }
2421 : done:
2422 8316155043 : if (track && isl_tab_track_bmap(tab, isl_basic_map_copy(bmap)) < 0)
2423 0 : goto error;
2424 8316155043 : return tab;
2425 : error:
2426 0 : isl_tab_free(tab);
2427 0 : return NULL;
2428 : }
2429 :
2430 1421072486 : __isl_give struct isl_tab *isl_tab_from_basic_set(
2431 : __isl_keep isl_basic_set *bset, int track)
2432 : {
2433 1421072486 : return isl_tab_from_basic_map(bset, track);
2434 : }
2435 :
2436 : /* Construct a tableau corresponding to the recession cone of "bset".
2437 : */
2438 893145 : 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 893145 : unsigned offset = 0;
2445 :
2446 893145 : if (!bset)
2447 0 : return NULL;
2448 893145 : if (parametric)
2449 893145 : offset = isl_basic_set_dim(bset, isl_dim_param);
2450 893145 : tab = isl_tab_alloc(bset->ctx, bset->n_eq + bset->n_ineq,
2451 893145 : isl_basic_set_total_dim(bset) - offset, 0);
2452 893145 : if (!tab)
2453 0 : return NULL;
2454 893145 : tab->rational = ISL_F_ISSET(bset, ISL_BASIC_SET_RATIONAL);
2455 893145 : tab->cone = 1;
2456 :
2457 893145 : isl_int_init(cst);
2458 893145 : isl_int_set_si(cst, 0);
2459 5275625 : for (i = 0; i < bset->n_eq; ++i) {
2460 4382480 : isl_int_swap(bset->eq[i][offset], cst);
2461 4382480 : if (offset > 0) {
2462 0 : if (isl_tab_add_eq(tab, bset->eq[i] + offset) < 0)
2463 0 : goto error;
2464 : } else
2465 4382480 : tab = add_eq(tab, bset->eq[i]);
2466 4382480 : isl_int_swap(bset->eq[i][offset], cst);
2467 4382480 : if (!tab)
2468 0 : goto done;
2469 : }
2470 962814 : for (i = 0; i < bset->n_ineq; ++i) {
2471 : int r;
2472 69669 : isl_int_swap(bset->ineq[i][offset], cst);
2473 69669 : r = isl_tab_add_row(tab, bset->ineq[i] + offset);
2474 69669 : isl_int_swap(bset->ineq[i][offset], cst);
2475 69669 : if (r < 0)
2476 0 : goto error;
2477 69669 : tab->con[r].is_nonneg = 1;
2478 69669 : if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
2479 0 : goto error;
2480 : }
2481 : done:
2482 893145 : isl_int_clear(cst);
2483 893145 : 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 893145 : isl_bool isl_tab_cone_is_bounded(struct isl_tab *tab)
2494 : {
2495 : int i;
2496 :
2497 893145 : if (!tab)
2498 0 : return isl_bool_error;
2499 893145 : if (tab->empty)
2500 0 : return isl_bool_true;
2501 893145 : if (tab->n_dead == tab->n_col)
2502 867902 : return isl_bool_true;
2503 :
2504 : for (;;) {
2505 229416 : for (i = tab->n_redundant; i < tab->n_row; ++i) {
2506 : struct isl_tab_var *var;
2507 : int sgn;
2508 114708 : var = isl_tab_var_from_row(tab, i);
2509 114708 : if (!var->is_nonneg)
2510 86523 : continue;
2511 28185 : sgn = sign_of_max(tab, var);
2512 28185 : if (sgn < -1)
2513 0 : return isl_bool_error;
2514 28185 : if (sgn != 0)
2515 0 : return isl_bool_false;
2516 28185 : if (close_row(tab, var, 0) < 0)
2517 0 : return isl_bool_error;
2518 28185 : break;
2519 : }
2520 28185 : if (tab->n_dead == tab->n_col)
2521 25243 : return isl_bool_true;
2522 2942 : if (i == tab->n_row)
2523 0 : return isl_bool_false;
2524 2942 : }
2525 : }
2526 :
2527 1632099111 : int isl_tab_sample_is_integer(struct isl_tab *tab)
2528 : {
2529 : int i;
2530 :
2531 1632099111 : if (!tab)
2532 0 : return -1;
2533 :
2534 9071895402 : for (i = 0; i < tab->n_var; ++i) {
2535 : int row;
2536 7570284094 : if (!tab->var[i].is_row)
2537 713829992 : continue;
2538 6856454102 : row = tab->var[i].index;
2539 6856454102 : if (!isl_int_is_divisible_by(tab->mat->row[row][1],
2540 : tab->mat->row[row][0]))
2541 130487803 : return 0;
2542 : }
2543 1501611308 : return 1;
2544 : }
2545 :
2546 1494243450 : static struct isl_vec *extract_integer_sample(struct isl_tab *tab)
2547 : {
2548 : int i;
2549 : struct isl_vec *vec;
2550 :
2551 1494243450 : vec = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
2552 1494243450 : if (!vec)
2553 0 : return NULL;
2554 :
2555 1494243450 : isl_int_set_si(vec->block.data[0], 1);
2556 8732120513 : for (i = 0; i < tab->n_var; ++i) {
2557 7237877063 : if (!tab->var[i].is_row)
2558 706113021 : isl_int_set_si(vec->block.data[1 + i], 0);
2559 : else {
2560 6531764042 : int row = tab->var[i].index;
2561 6531764042 : isl_int_divexact(vec->block.data[1 + i],
2562 : tab->mat->row[row][1], tab->mat->row[row][0]);
2563 : }
2564 : }
2565 :
2566 1494243450 : return vec;
2567 : }
2568 :
2569 105541573 : 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 105541573 : if (!tab)
2576 0 : return NULL;
2577 :
2578 105541573 : vec = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
2579 105541573 : if (!vec)
2580 0 : return NULL;
2581 :
2582 105541573 : isl_int_init(m);
2583 :
2584 105541573 : isl_int_set_si(vec->block.data[0], 1);
2585 669677594 : for (i = 0; i < tab->n_var; ++i) {
2586 : int row;
2587 564136021 : if (!tab->var[i].is_row) {
2588 127960600 : isl_int_set_si(vec->block.data[1 + i], 0);
2589 127960600 : continue;
2590 : }
2591 436175421 : row = tab->var[i].index;
2592 436175421 : isl_int_gcd(m, vec->block.data[0], tab->mat->row[row][0]);
2593 436175421 : isl_int_divexact(m, tab->mat->row[row][0], m);
2594 436175421 : isl_seq_scale(vec->block.data, vec->block.data, m, 1 + i);
2595 436175421 : isl_int_divexact(m, vec->block.data[0], tab->mat->row[row][0]);
2596 436175421 : isl_int_mul(vec->block.data[1 + i], m, tab->mat->row[row][1]);
2597 : }
2598 105541573 : vec = isl_vec_normalize(vec);
2599 :
2600 105541573 : isl_int_clear(m);
2601 105541573 : 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 148929893 : static void get_rounded_sample_value(struct isl_tab *tab,
2608 : struct isl_tab_var *var, int sgn, isl_int *v)
2609 : {
2610 148929893 : if (!var->is_row)
2611 1278282 : isl_int_set_si(*v, 0);
2612 147651611 : else if (sgn > 0)
2613 142439487 : isl_int_cdiv_q(*v, tab->mat->row[var->index][1],
2614 : tab->mat->row[var->index][0]);
2615 : else
2616 5212124 : isl_int_fdiv_q(*v, tab->mat->row[var->index][1],
2617 : tab->mat->row[var->index][0]);
2618 148929893 : }
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 4628748914 : 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 4628748914 : if (!bmap)
2635 0 : return NULL;
2636 4628748914 : if (!tab)
2637 0 : return bmap;
2638 :
2639 4628748914 : n_eq = tab->n_eq;
2640 4628748914 : if (tab->empty)
2641 2296485 : bmap = isl_basic_map_set_to_empty(bmap);
2642 : else
2643 20347481252 : for (i = bmap->n_ineq - 1; i >= 0; --i) {
2644 15721028823 : if (isl_tab_is_equality(tab, n_eq + i))
2645 10507856375 : isl_basic_map_inequality_to_equality(bmap, i);
2646 5213172448 : else if (isl_tab_is_redundant(tab, n_eq + i))
2647 1954620804 : isl_basic_map_drop_inequality(bmap, i);
2648 : }
2649 4628748914 : if (bmap->n_eq != n_eq)
2650 1346162311 : bmap = isl_basic_map_gauss(bmap, NULL);
2651 4628748914 : if (!tab->rational &&
2652 1742127323 : bmap && !bmap->sample && isl_tab_sample_is_integer(tab))
2653 1494243450 : bmap->sample = extract_integer_sample(tab);
2654 4628748914 : return bmap;
2655 : }
2656 :
2657 1402391 : struct isl_basic_set *isl_basic_set_update_from_tab(struct isl_basic_set *bset,
2658 : struct isl_tab *tab)
2659 : {
2660 1402391 : 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 3524788 : static isl_stat drop_last_con_in_row(struct isl_tab *tab, int r)
2668 : {
2669 3524788 : 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 3524788 : 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 3524788 : if (drop_row(tab, tab->con[r].index) < 0)
2677 0 : return isl_stat_error;
2678 :
2679 3524788 : 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 3524788 : 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 3524788 : unsigned off = 2 + tab->M;
2698 :
2699 3524788 : if (var->is_zero)
2700 0 : return isl_stat_ok;
2701 3524788 : 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 3524788 : if (isl_tab_extend_cons(tab, 1) < 0)
2707 0 : return isl_stat_error;
2708 :
2709 3524788 : r = tab->n_con;
2710 3524788 : tab->con[r].index = tab->n_row;
2711 3524788 : tab->con[r].is_row = 1;
2712 3524788 : tab->con[r].is_nonneg = 0;
2713 3524788 : tab->con[r].is_zero = 0;
2714 3524788 : tab->con[r].is_redundant = 0;
2715 3524788 : tab->con[r].frozen = 0;
2716 3524788 : tab->con[r].negated = 0;
2717 3524788 : tab->row_var[tab->n_row] = ~r;
2718 3524788 : row = tab->mat->row[tab->n_row];
2719 :
2720 3524788 : if (var->is_row) {
2721 2929919 : isl_int_set(row[0], tab->mat->row[var->index][0]);
2722 5859838 : isl_seq_neg(row + 1,
2723 5859838 : tab->mat->row[var->index] + 1, 1 + tab->n_col);
2724 : } else {
2725 594869 : isl_int_set_si(row[0], 1);
2726 594869 : isl_seq_clr(row + 1, 1 + tab->n_col);
2727 594869 : isl_int_set_si(row[off + var->index], -1);
2728 : }
2729 :
2730 3524788 : tab->n_row++;
2731 3524788 : tab->n_con++;
2732 :
2733 3524788 : sgn = sign_of_max(tab, &tab->con[r]);
2734 3524788 : if (sgn < -1)
2735 0 : return isl_stat_error;
2736 3524788 : if (sgn < 0) {
2737 45255 : if (drop_last_con_in_row(tab, r) < 0)
2738 0 : return isl_stat_error;
2739 45255 : if (isl_tab_mark_empty(tab) < 0)
2740 0 : return isl_stat_error;
2741 45255 : return isl_stat_ok;
2742 : }
2743 3479533 : tab->con[r].is_nonneg = 1;
2744 : /* sgn == 0 */
2745 3479533 : if (close_row(tab, &tab->con[r], 1) < 0)
2746 0 : return isl_stat_error;
2747 3479533 : if (drop_last_con_in_row(tab, r) < 0)
2748 0 : return isl_stat_error;
2749 :
2750 3479533 : return isl_stat_ok;
2751 : }
2752 :
2753 : /* Check that "con" is a valid constraint position for "tab".
2754 : */
2755 6556637197 : static isl_stat isl_tab_check_con(struct isl_tab *tab, int con)
2756 : {
2757 6556637197 : if (!tab)
2758 0 : return isl_stat_error;
2759 6556637197 : 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 6556637197 : 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 498520 : int isl_tab_relax(struct isl_tab *tab, int con)
2783 : {
2784 : struct isl_tab_var *var;
2785 :
2786 498520 : if (!tab)
2787 0 : return -1;
2788 :
2789 498520 : var = &tab->con[con];
2790 :
2791 498520 : 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 498520 : 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 498520 : if (!var->is_row && !max_is_manifestly_unbounded(tab, var))
2799 55155 : if (to_row(tab, var, 1) < 0)
2800 0 : return -1;
2801 498520 : if (!var->is_row && !min_is_manifestly_unbounded(tab, var))
2802 42904 : if (to_row(tab, var, -1) < 0)
2803 0 : return -1;
2804 :
2805 498520 : if (var->is_row) {
2806 99838 : isl_int_add(tab->mat->row[var->index][1],
2807 : tab->mat->row[var->index][1], tab->mat->row[var->index][0]);
2808 99838 : if (restore_row(tab, var) < 0)
2809 0 : return -1;
2810 : } else {
2811 : int i;
2812 398682 : unsigned off = 2 + tab->M;
2813 :
2814 1641810 : for (i = 0; i < tab->n_row; ++i) {
2815 1243128 : if (isl_int_is_zero(tab->mat->row[i][off + var->index]))
2816 721901 : continue;
2817 521227 : 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 498520 : if (isl_tab_push_var(tab, isl_tab_undo_relax, var) < 0)
2824 0 : return -1;
2825 :
2826 498520 : 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 1192261 : int isl_tab_unrestrict(struct isl_tab *tab, int con)
2893 : {
2894 : struct isl_tab_var *var;
2895 :
2896 1192261 : if (!tab)
2897 0 : return -1;
2898 :
2899 1192261 : var = &tab->con[con];
2900 1192261 : if (!var->is_nonneg)
2901 0 : return 0;
2902 :
2903 1192261 : var->is_nonneg = 0;
2904 1192261 : if (isl_tab_push_var(tab, isl_tab_undo_unrestrict, var) < 0)
2905 0 : return -1;
2906 :
2907 1192261 : return 0;
2908 : }
2909 :
2910 689080 : int isl_tab_select_facet(struct isl_tab *tab, int con)
2911 : {
2912 689080 : if (!tab)
2913 0 : return -1;
2914 :
2915 689080 : return cut_to_hyperplane(tab, &tab->con[con]);
2916 : }
2917 :
2918 31575684801 : static int may_be_equality(struct isl_tab *tab, int row)
2919 : {
2920 63208821659 : return tab->rational ? isl_int_is_zero(tab->mat->row[row][1])
2921 31633136858 : : 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 10424363111 : static struct isl_tab_var *select_marked(struct isl_tab *tab)
2937 : {
2938 : int i;
2939 : struct isl_tab_var *var;
2940 :
2941 61918871319 : for (i = tab->n_con - 1; i >= 0; --i) {
2942 60482539757 : var = &tab->con[i];
2943 60482539757 : if (var->index < 0)
2944 39382802951 : continue;
2945 21099736806 : if (var->is_row && var->index < tab->n_redundant)
2946 1667711494 : continue;
2947 19432025312 : if (!var->is_row && var->index < tab->n_dead)
2948 85385 : continue;
2949 19431939927 : if (var->marked)
2950 8988031549 : return var;
2951 : }
2952 :
2953 1436331562 : 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 4565119413 : int isl_tab_detect_implicit_equalities(struct isl_tab *tab)
2974 : {
2975 : int i;
2976 : unsigned n_marked;
2977 :
2978 4565119413 : if (!tab)
2979 0 : return -1;
2980 4565119413 : if (tab->empty)
2981 2445752 : return 0;
2982 4562673661 : if (tab->n_dead == tab->n_col)
2983 2851251002 : return 0;
2984 :
2985 1711422659 : n_marked = 0;
2986 16629144516 : for (i = tab->n_redundant; i < tab->n_row; ++i) {
2987 14917721857 : struct isl_tab_var *var = isl_tab_var_from_row(tab, i);
2988 29244624354 : var->marked = !var->frozen && var->is_nonneg &&
2989 14326902497 : may_be_equality(tab, i);
2990 14917721857 : if (var->marked)
2991 13649316650 : n_marked++;
2992 : }
2993 9813452965 : for (i = tab->n_dead; i < tab->n_col; ++i) {
2994 8102030306 : struct isl_tab_var *var = var_from_col(tab, i);
2995 8102030306 : var->marked = !var->frozen && var->is_nonneg;
2996 8102030306 : if (var->marked)
2997 710639129 : n_marked++;
2998 : }
2999 11215522856 : while (n_marked) {
3000 : struct isl_tab_var *var;
3001 : int sgn;
3002 9181894076 : var = select_marked(tab);
3003 9181894076 : if (!var)
3004 1386380830 : break;
3005 7795513246 : var->marked = 0;
3006 7795513246 : n_marked--;
3007 7795513246 : sgn = sign_of_max(tab, var);
3008 7795513246 : if (sgn < 0)
3009 0 : return -1;
3010 7795513246 : if (sgn == 0) {
3011 5322705501 : if (close_row(tab, var, 0) < 0)
3012 0 : return -1;
3013 2472807745 : } else if (!tab->rational && !at_least_one(tab, var)) {
3014 2835708 : if (cut_to_hyperplane(tab, var) < 0)
3015 0 : return -1;
3016 2835708 : return isl_tab_detect_implicit_equalities(tab);
3017 : }
3018 50655476354 : for (i = tab->n_redundant; i < tab->n_row; ++i) {
3019 42862798816 : var = isl_tab_var_from_row(tab, i);
3020 42862798816 : if (!var->marked)
3021 25614016512 : continue;
3022 17248782304 : if (may_be_equality(tab, i))
3023 17225889560 : continue;
3024 22892744 : var->marked = 0;
3025 22892744 : n_marked--;
3026 : }
3027 : }
3028 :
3029 1708586951 : 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 1652116 : static int update_con_after_move(struct isl_tab *tab, int i, int old)
3036 : {
3037 : int *p;
3038 : int index;
3039 :
3040 1652116 : index = tab->con[i].index;
3041 1652116 : if (index == -1)
3042 1175711 : return 0;
3043 476405 : p = tab->con[i].is_row ? tab->row_var : tab->col_var;
3044 476405 : if (p[index] != ~old)
3045 0 : isl_die(tab->mat->ctx, isl_error_internal,
3046 : "broken internal state", return -1);
3047 476405 : p[index] = ~i;
3048 :
3049 476405 : 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 65215 : isl_stat isl_tab_swap_constraints(struct isl_tab *tab, int con1, int con2)
3058 : {
3059 : struct isl_tab_var var;
3060 :
3061 130430 : if (isl_tab_check_con(tab, con1) < 0 ||
3062 65215 : isl_tab_check_con(tab, con2) < 0)
3063 0 : return isl_stat_error;
3064 :
3065 65215 : var = tab->con[con1];
3066 65215 : tab->con[con1] = tab->con[con2];
3067 65215 : if (update_con_after_move(tab, con1, con2) < 0)
3068 0 : return isl_stat_error;
3069 65215 : tab->con[con2] = var;
3070 65215 : if (update_con_after_move(tab, con2, con1) < 0)
3071 0 : return isl_stat_error;
3072 :
3073 65215 : 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 364552 : 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 364552 : if (n <= 1)
3085 116780 : return 0;
3086 :
3087 247772 : last = first + n - 1;
3088 247772 : var = tab->con[last];
3089 1362075 : for (i = last; i > first; --i) {
3090 1114303 : tab->con[i] = tab->con[i - 1];
3091 1114303 : if (update_con_after_move(tab, i, i - 1) < 0)
3092 0 : return -1;
3093 : }
3094 247772 : tab->con[first] = var;
3095 247772 : if (update_con_after_move(tab, first, last) < 0)
3096 0 : return -1;
3097 :
3098 247772 : 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 52099 : static isl_stat con_drop_entries(struct isl_tab *tab,
3107 : unsigned first, unsigned n)
3108 : {
3109 : int i;
3110 :
3111 52099 : 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 52099 : tab->n_con -= n;
3116 :
3117 211710 : for (i = first; i < tab->n_con; ++i) {
3118 159611 : tab->con[i] = tab->con[i + n];
3119 159611 : if (update_con_after_move(tab, i, i + n) < 0)
3120 0 : return isl_stat_error;
3121 : }
3122 :
3123 52099 : 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 65215 : static isl_stat swap_eq(unsigned a, unsigned b, void *user)
3131 : {
3132 65215 : struct isl_tab *tab = user;
3133 :
3134 65215 : 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 52099 : static isl_stat drop_eq(unsigned n, void *user)
3146 : {
3147 52099 : struct isl_tab *tab = user;
3148 :
3149 52099 : if (tab->n_eq == n)
3150 0 : return isl_tab_mark_empty(tab);
3151 :
3152 52099 : tab->n_eq -= n;
3153 52099 : 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 55490 : 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 55490 : single = isl_basic_map_has_single_reference(bmap);
3165 55490 : if (single < 0)
3166 0 : return isl_basic_map_free(bmap);
3167 55490 : if (single)
3168 807 : return bmap;
3169 54683 : 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 2855809721 : __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 2855809721 : if (!tab || !bmap)
3212 0 : return isl_basic_map_free(bmap);
3213 2855809721 : if (tab->empty)
3214 38306 : return bmap;
3215 :
3216 2855771415 : n_eq = tab->n_eq;
3217 2896675534 : for (i = bmap->n_ineq - 1; i >= 0; --i) {
3218 40904119 : if (!isl_tab_is_equality(tab, bmap->n_eq + i))
3219 40721843 : continue;
3220 182276 : isl_basic_map_inequality_to_equality(bmap, i);
3221 182276 : if (rotate_constraints(tab, 0, tab->n_eq + i + 1) < 0)
3222 0 : return isl_basic_map_free(bmap);
3223 182276 : if (rotate_constraints(tab, tab->n_eq + i + 1,
3224 182276 : bmap->n_ineq - i) < 0)
3225 0 : return isl_basic_map_free(bmap);
3226 182276 : tab->n_eq++;
3227 : }
3228 :
3229 2855771415 : if (n_eq != tab->n_eq)
3230 55490 : bmap = gauss_if_shared(bmap, tab);
3231 :
3232 2855771415 : return bmap;
3233 : }
3234 :
3235 12561329872 : static int con_is_redundant(struct isl_tab *tab, struct isl_tab_var *var)
3236 : {
3237 12561329872 : if (!tab)
3238 0 : return -1;
3239 12561329872 : if (tab->rational) {
3240 11221856972 : int sgn = sign_of_min(tab, var);
3241 11221856972 : if (sgn < -1)
3242 0 : return -1;
3243 11221856972 : return sgn >= 0;
3244 : } else {
3245 1339472900 : int irred = isl_tab_min_at_most_neg_one(tab, var);
3246 1339472900 : if (irred < 0)
3247 0 : return -1;
3248 1339472900 : 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 3218490093 : int isl_tab_detect_redundant(struct isl_tab *tab)
3266 : {
3267 : int i;
3268 : unsigned n_marked;
3269 :
3270 3218490093 : if (!tab)
3271 0 : return -1;
3272 3218490093 : if (tab->empty)
3273 445789 : return 0;
3274 3218044304 : if (tab->n_redundant == tab->n_row)
3275 998202 : return 0;
3276 :
3277 3217046102 : n_marked = 0;
3278 10723115629 : for (i = tab->n_redundant; i < tab->n_row; ++i) {
3279 7506069527 : struct isl_tab_var *var = isl_tab_var_from_row(tab, i);
3280 7506069527 : var->marked = !var->frozen && var->is_nonneg;
3281 7506069527 : if (var->marked)
3282 970422210 : n_marked++;
3283 : }
3284 5079608521 : for (i = tab->n_dead; i < tab->n_col; ++i) {
3285 1862562419 : struct isl_tab_var *var = var_from_col(tab, i);
3286 2704152933 : var->marked = !var->frozen && var->is_nonneg &&
3287 841590514 : !min_is_manifestly_unbounded(tab, var);
3288 1862562419 : if (var->marked)
3289 480823849 : n_marked++;
3290 : }
3291 7626610507 : while (n_marked) {
3292 : struct isl_tab_var *var;
3293 : int red;
3294 1242469035 : var = select_marked(tab);
3295 1242469035 : if (!var)
3296 49950732 : break;
3297 1192518303 : var->marked = 0;
3298 1192518303 : n_marked--;
3299 1192518303 : red = con_is_redundant(tab, var);
3300 1192518303 : if (red < 0)
3301 0 : return -1;
3302 1192518303 : if (red && !var->is_redundant)
3303 19008934 : if (isl_tab_mark_redundant(tab, var->index) < 0)
3304 0 : return -1;
3305 7588001692 : for (i = tab->n_dead; i < tab->n_col; ++i) {
3306 6395483389 : var = var_from_col(tab, i);
3307 6395483389 : if (!var->marked)
3308 5017990586 : continue;
3309 1377492803 : if (!min_is_manifestly_unbounded(tab, var))
3310 1269371516 : continue;
3311 108121287 : var->marked = 0;
3312 108121287 : n_marked--;
3313 : }
3314 : }
3315 :
3316 3217046102 : return 0;
3317 : }
3318 :
3319 15762431462 : int isl_tab_is_equality(struct isl_tab *tab, int con)
3320 : {
3321 : int row;
3322 : unsigned off;
3323 :
3324 15762431462 : if (!tab)
3325 0 : return -1;
3326 15762431462 : if (tab->con[con].is_zero)
3327 10508038647 : return 1;
3328 5254392815 : if (tab->con[con].is_redundant)
3329 1955719441 : return 0;
3330 3298673374 : if (!tab->con[con].is_row)
3331 2361878517 : return tab->con[con].index < tab->n_dead;
3332 :
3333 936794857 : row = tab->con[con].index;
3334 :
3335 936794857 : off = 2 + tab->M;
3336 3170548785 : return isl_int_is_zero(tab->mat->row[row][1]) &&
3337 1658157197 : !row_is_big(tab, row) &&
3338 360681170 : isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
3339 360681170 : 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 2516783001 : 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 2516783001 : enum isl_lp_result res = isl_lp_ok;
3358 : struct isl_tab_var *var;
3359 : struct isl_tab_undo *snap;
3360 :
3361 2516783001 : if (!tab)
3362 0 : return isl_lp_error;
3363 :
3364 2516783001 : if (tab->empty)
3365 0 : return isl_lp_empty;
3366 :
3367 2516783001 : snap = isl_tab_snap(tab);
3368 2516783001 : r = isl_tab_add_row(tab, f);
3369 2516783001 : if (r < 0)
3370 0 : return isl_lp_error;
3371 2516783001 : var = &tab->con[r];
3372 : for (;;) {
3373 : int row, col;
3374 4234557161 : find_pivot(tab, var, var, -1, &row, &col);
3375 4234557161 : if (row == var->index) {
3376 570940 : res = isl_lp_unbounded;
3377 2517353941 : break;
3378 : }
3379 4233986221 : if (row == -1)
3380 2516212061 : break;
3381 1717774160 : if (isl_tab_pivot(tab, row, col) < 0)
3382 0 : return isl_lp_error;
3383 1717774160 : }
3384 2516783001 : isl_int_mul(tab->mat->row[var->index][0],
3385 : tab->mat->row[var->index][0], denom);
3386 2516783001 : if (ISL_FL_ISSET(flags, ISL_TAB_SAVE_DUAL)) {
3387 : int i;
3388 :
3389 54971278 : isl_vec_free(tab->dual);
3390 54971278 : tab->dual = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_con);
3391 54971278 : if (!tab->dual)
3392 0 : return isl_lp_error;
3393 54971278 : isl_int_set(tab->dual->el[0], tab->mat->row[var->index][0]);
3394 2625081295 : for (i = 0; i < tab->n_con; ++i) {
3395 : int pos;
3396 2570110017 : if (tab->con[i].is_row) {
3397 1967546236 : isl_int_set_si(tab->dual->el[1 + i], 0);
3398 1967546236 : continue;
3399 : }
3400 602563781 : pos = 2 + tab->M + tab->con[i].index;
3401 602563781 : if (tab->con[i].negated)
3402 96619504 : isl_int_neg(tab->dual->el[1 + i],
3403 : tab->mat->row[var->index][pos]);
3404 : else
3405 505944277 : isl_int_set(tab->dual->el[1 + i],
3406 : tab->mat->row[var->index][pos]);
3407 : }
3408 : }
3409 2516783001 : if (opt && res == isl_lp_ok) {
3410 2516212061 : if (opt_denom) {
3411 2375171027 : isl_int_set(*opt, tab->mat->row[var->index][1]);
3412 2375171027 : isl_int_set(*opt_denom, tab->mat->row[var->index][0]);
3413 : } else
3414 141041034 : get_rounded_sample_value(tab, var, 1, opt);
3415 : }
3416 2516783001 : if (isl_tab_rollback(tab, snap) < 0)
3417 0 : return isl_lp_error;
3418 2516783001 : 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 6556506767 : int isl_tab_is_redundant(struct isl_tab *tab, int con)
3430 : {
3431 6556506767 : if (isl_tab_check_con(tab, con) < 0)
3432 0 : return -1;
3433 6556506767 : if (tab->con[con].is_zero)
3434 0 : return 0;
3435 6556506767 : if (tab->con[con].is_redundant)
3436 2912018833 : return 1;
3437 3644487934 : 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 6919671 : static isl_bool is_constant(struct isl_tab *tab, struct isl_tab_var *var,
3449 : isl_int *value)
3450 : {
3451 6919671 : unsigned off = 2 + tab->M;
3452 6919671 : isl_mat *mat = tab->mat;
3453 : int n;
3454 : int row;
3455 : int pos;
3456 :
3457 6919671 : if (!var->is_row)
3458 1278282 : return isl_bool_false;
3459 5641389 : row = var->index;
3460 5641389 : if (row_is_big(tab, row))
3461 0 : return isl_bool_false;
3462 5641389 : n = tab->n_col - tab->n_dead;
3463 5641389 : pos = isl_seq_first_non_zero(mat->row[row] + off + tab->n_dead, n);
3464 5641389 : if (pos != -1)
3465 5209068 : return isl_bool_false;
3466 432321 : if (value)
3467 0 : isl_int_divexact(*value, mat->row[row][1], mat->row[row][0]);
3468 432321 : 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 9829810 : static int reached(struct isl_tab *tab, struct isl_tab_var *var, int sgn,
3482 : isl_int target, isl_int *tmp)
3483 : {
3484 9829810 : if (row_is_big(tab, var->index))
3485 0 : return 1;
3486 9829810 : isl_int_mul(*tmp, tab->mat->row[var->index][0], target);
3487 9829810 : if (sgn > 0)
3488 2207943 : return isl_int_ge(tab->mat->row[var->index][1], *tmp);
3489 : else
3490 7621867 : 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 7885803 : 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 7885803 : if (sgn < 0 && min_is_manifestly_unbounded(tab, var))
3515 1248600 : return isl_bool_true;
3516 6637203 : if (sgn > 0 && max_is_manifestly_unbounded(tab, var))
3517 0 : return isl_bool_true;
3518 6637203 : if (to_row(tab, var, sgn) < 0)
3519 0 : return isl_bool_error;
3520 16467013 : while (!reached(tab, var, sgn, target, tmp)) {
3521 7926609 : find_pivot(tab, var, var, sgn, &row, &col);
3522 7926609 : if (row == -1)
3523 1401509 : return isl_bool_false;
3524 6525100 : if (row == var->index)
3525 3332493 : return isl_bool_true;
3526 3192607 : if (isl_tab_pivot(tab, row, col) < 0)
3527 0 : return isl_bool_error;
3528 : }
3529 :
3530 1903201 : 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 6487350 : 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 6487350 : get_rounded_sample_value(tab, var, -1, target);
3561 6487350 : isl_int_sub_ui(*target, *target, 1);
3562 6487350 : reached = var_reaches(tab, var, -1, *target, tmp);
3563 6487350 : if (reached < 0 || reached)
3564 5088897 : return isl_bool_not(reached);
3565 1398453 : get_rounded_sample_value(tab, var, 1, target);
3566 1398453 : isl_int_add_ui(*target, *target, 1);
3567 1398453 : reached = var_reaches(tab, var, 1, *target, tmp);
3568 1398453 : if (reached < 0 || reached)
3569 1395397 : return isl_bool_not(reached);
3570 3056 : get_rounded_sample_value(tab, var, -1, tmp);
3571 3056 : isl_int_sub_ui(*target, *target, 1);
3572 3056 : 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 3056 : if (isl_tab_extend_cons(tab, 1) < 0)
3579 0 : return isl_bool_error;
3580 3056 : eq = isl_vec_alloc(isl_tab_get_ctx(tab), 1 + tab->n_var);
3581 3056 : if (!eq)
3582 0 : return isl_bool_error;
3583 3056 : pos = var - tab->var;
3584 3056 : isl_seq_clr(eq->el + 1, tab->n_var);
3585 3056 : isl_int_set_si(eq->el[1 + pos], -1);
3586 3056 : isl_int_set(eq->el[0], *target);
3587 3056 : r = isl_tab_add_eq(tab, eq->el);
3588 3056 : isl_vec_free(eq);
3589 :
3590 3056 : 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 6919671 : 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 6919671 : if (var->is_row && row_is_big(tab, var->index))
3613 0 : return isl_bool_false;
3614 6919671 : is_cst = is_constant(tab, var, value);
3615 6919671 : if (is_cst < 0 || is_cst)
3616 432321 : return is_cst;
3617 :
3618 6487350 : if (!value)
3619 6487350 : isl_int_init(target);
3620 6487350 : isl_int_init(tmp);
3621 :
3622 6487350 : is_cst = detect_constant_with_tmp(tab, var,
3623 : value ? value : &target, &tmp);
3624 :
3625 6487350 : isl_int_clear(tmp);
3626 6487350 : if (!value)
3627 6487350 : isl_int_clear(target);
3628 :
3629 6487350 : 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 1192261 : isl_stat isl_tab_detect_constants(struct isl_tab *tab)
3659 : {
3660 : int i;
3661 :
3662 1192261 : if (!tab)
3663 0 : return isl_stat_error;
3664 1192261 : if (tab->rational)
3665 0 : return isl_stat_ok;
3666 :
3667 8111932 : for (i = 0; i < tab->n_var; ++i) {
3668 6919671 : if (get_constant(tab, &tab->var[i], NULL) < 0)
3669 0 : return isl_stat_error;
3670 : }
3671 :
3672 1192261 : 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 35136170521 : struct isl_tab_undo *isl_tab_snap(struct isl_tab *tab)
3679 : {
3680 35136170521 : if (!tab)
3681 0 : return NULL;
3682 35136170521 : tab->need_undo = 1;
3683 35136170521 : 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 435685 : static isl_stat unrelax(struct isl_tab *tab, struct isl_tab_var *var)
3716 : {
3717 435685 : unsigned off = 2 + tab->M;
3718 :
3719 435685 : if (!var->is_row && !max_is_manifestly_unbounded(tab, var))
3720 48434 : if (to_row(tab, var, 1) < 0)
3721 0 : return isl_stat_error;
3722 :
3723 435685 : if (var->is_row) {
3724 54092 : isl_int_sub(tab->mat->row[var->index][1],
3725 : tab->mat->row[var->index][1], tab->mat->row[var->index][0]);
3726 54092 : if (var->is_nonneg) {
3727 54092 : int sgn = restore_row(tab, var);
3728 54092 : isl_assert(tab->mat->ctx, sgn >= 0,
3729 : return isl_stat_error);
3730 : }
3731 : } else {
3732 : int i;
3733 :
3734 1568344 : for (i = 0; i < tab->n_row; ++i) {
3735 1186751 : if (isl_int_is_zero(tab->mat->row[i][off + var->index]))
3736 674766 : continue;
3737 511985 : 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 435685 : 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 1190442 : static isl_stat ununrestrict(struct isl_tab *tab, struct isl_tab_var *var)
3752 : {
3753 1190442 : var->is_nonneg = 1;
3754 :
3755 1190442 : if (var->is_row && restore_row(tab, var) < -1)
3756 0 : return isl_stat_error;
3757 :
3758 1190442 : 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 13044641571 : static isl_stat restore_last_redundant(struct isl_tab *tab)
3769 : {
3770 : struct isl_tab_var *var;
3771 :
3772 13044641571 : 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 13044641571 : var = isl_tab_var_from_row(tab, tab->n_redundant - 1);
3777 13044641571 : var->is_redundant = 0;
3778 13044641571 : tab->n_redundant--;
3779 13044641571 : restore_row(tab, var);
3780 :
3781 13044641571 : 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 >14904*10^7 : static isl_stat perform_undo_var(struct isl_tab *tab, struct isl_tab_undo *undo)
3787 : {
3788 >14904*10^7 : struct isl_tab_var *var = var_from_index(tab, undo->u.var_index);
3789 >14904*10^7 : switch (undo->type) {
3790 : case isl_tab_undo_nonneg:
3791 24622380636 : var->is_nonneg = 0;
3792 24622380636 : break;
3793 : case isl_tab_undo_redundant:
3794 12709232757 : 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 12709232757 : return restore_last_redundant(tab);
3799 : case isl_tab_undo_freeze:
3800 70090799829 : var->frozen = 0;
3801 70090799829 : break;
3802 : case isl_tab_undo_zero:
3803 46076049 : var->is_zero = 0;
3804 46076049 : if (!var->is_row)
3805 46076049 : tab->n_dead--;
3806 46076049 : break;
3807 : case isl_tab_undo_allocate:
3808 41575772102 : 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 41575772102 : if (!var->is_row) {
3814 2036063867 : if (!max_is_manifestly_unbounded(tab, var)) {
3815 999293467 : if (to_row(tab, var, 1) < 0)
3816 0 : return isl_stat_error;
3817 1036770400 : } else if (!min_is_manifestly_unbounded(tab, var)) {
3818 799036751 : if (to_row(tab, var, -1) < 0)
3819 0 : return isl_stat_error;
3820 : } else
3821 237733649 : if (to_row(tab, var, 0) < 0)
3822 0 : return isl_stat_error;
3823 : }
3824 41575772102 : return drop_row(tab, var->index);
3825 : case isl_tab_undo_relax:
3826 435685 : return unrelax(tab, var);
3827 : case isl_tab_undo_unrestrict:
3828 1190442 : 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 94759256514 : 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 205381885 : isl_stat isl_tab_restore_redundant(struct isl_tab *tab)
3846 : {
3847 205381885 : if (!tab)
3848 0 : return isl_stat_error;
3849 :
3850 205381885 : 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 746172584 : while (tab->n_redundant > 0) {
3857 335408814 : if (tab->row_var[tab->n_redundant - 1] >= 0) {
3858 : struct isl_tab_var *var;
3859 :
3860 264449857 : var = isl_tab_var_from_row(tab, tab->n_redundant - 1);
3861 264449857 : var->is_nonneg = 0;
3862 : }
3863 335408814 : restore_last_redundant(tab);
3864 : }
3865 205381885 : 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 >17638*10^7 : static isl_stat perform_undo(struct isl_tab *tab, struct isl_tab_undo *undo)
3968 : {
3969 >17638*10^7 : switch (undo->type) {
3970 : case isl_tab_undo_rational:
3971 146137 : tab->rational = 0;
3972 146137 : break;
3973 : case isl_tab_undo_empty:
3974 6822784530 : tab->empty = 0;
3975 6822784530 : 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 >14904*10^7 : 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 20515829007 : 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 6822930667 : 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 28076302690 : int isl_tab_rollback(struct isl_tab *tab, struct isl_tab_undo *snap)
4012 : {
4013 : struct isl_tab_undo *undo, *next;
4014 :
4015 28076302690 : if (!tab)
4016 0 : return -1;
4017 :
4018 28076302690 : tab->in_undo = 1;
4019 >20446*10^7 : for (undo = tab->top; undo && undo != &tab->bottom; undo = next) {
4020 >18182*10^7 : next = undo->next;
4021 >18182*10^7 : if (undo == snap)
4022 5442186391 : break;
4023 >17638*10^7 : 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 >17638*10^7 : free_undo_record(undo);
4030 : }
4031 28076302690 : tab->in_undo = 0;
4032 28076302690 : tab->top = undo;
4033 28076302690 : if (!undo)
4034 0 : return -1;
4035 28076302690 : 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 7099033178 : static enum isl_ineq_type separation_type(struct isl_tab *tab, unsigned row)
4049 : {
4050 : int pos;
4051 7099033178 : unsigned off = 2 + tab->M;
4052 :
4053 7099033178 : if (tab->rational)
4054 7037094675 : return isl_ineq_separate;
4055 :
4056 61938503 : if (!isl_int_is_one(tab->mat->row[row][0]))
4057 4262939 : return isl_ineq_separate;
4058 :
4059 57675564 : pos = isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
4060 57675564 : tab->n_col - tab->n_dead);
4061 57675564 : if (pos == -1) {
4062 27605212 : if (isl_int_is_negone(tab->mat->row[row][1]))
4063 2987884 : return isl_ineq_adj_eq;
4064 : else
4065 24617328 : return isl_ineq_separate;
4066 : }
4067 :
4068 30070352 : if (!isl_int_eq(tab->mat->row[row][1],
4069 : tab->mat->row[row][off + tab->n_dead + pos]))
4070 9281613 : return isl_ineq_separate;
4071 :
4072 41577478 : pos = isl_seq_first_non_zero(
4073 20788739 : tab->mat->row[row] + off + tab->n_dead + pos + 1,
4074 20788739 : tab->n_col - tab->n_dead - pos - 1);
4075 :
4076 20788739 : 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 18491485102 : enum isl_ineq_type isl_tab_ineq_type(struct isl_tab *tab, isl_int *ineq)
4088 : {
4089 18491485102 : enum isl_ineq_type type = isl_ineq_error;
4090 18491485102 : struct isl_tab_undo *snap = NULL;
4091 : int con;
4092 : int row;
4093 :
4094 18491485102 : if (!tab)
4095 0 : return isl_ineq_error;
4096 :
4097 18491485102 : if (isl_tab_extend_cons(tab, 1) < 0)
4098 0 : return isl_ineq_error;
4099 :
4100 18491485102 : snap = isl_tab_snap(tab);
4101 :
4102 18491485102 : con = isl_tab_add_row(tab, ineq);
4103 18491485102 : if (con < 0)
4104 0 : goto error;
4105 :
4106 18491485102 : row = tab->con[con].index;
4107 18491485102 : if (isl_tab_row_is_redundant(tab, row))
4108 0 : type = isl_ineq_redundant;
4109 25615972940 : else if (isl_int_is_neg(tab->mat->row[row][1]) &&
4110 7207667426 : (tab->rational ||
4111 83179588 : isl_int_abs_ge(tab->mat->row[row][1],
4112 7122673533 : tab->mat->row[row][0]))) {
4113 7122673533 : int nonneg = at_least_zero(tab, &tab->con[con]);
4114 7122673533 : if (nonneg < 0)
4115 0 : goto error;
4116 7122673533 : if (nonneg)
4117 23640355 : type = isl_ineq_cut;
4118 : else
4119 7099033178 : type = separation_type(tab, row);
4120 : } else {
4121 11368811569 : int red = con_is_redundant(tab, &tab->con[con]);
4122 11368811569 : if (red < 0)
4123 0 : goto error;
4124 11368811569 : if (!red)
4125 37976953 : type = isl_ineq_cut;
4126 : else
4127 11330834616 : type = isl_ineq_redundant;
4128 : }
4129 :
4130 18491485102 : if (isl_tab_rollback(tab, snap))
4131 0 : return isl_ineq_error;
4132 18491485102 : return type;
4133 : error:
4134 0 : return isl_ineq_error;
4135 : }
4136 :
4137 1383503511 : isl_stat isl_tab_track_bmap(struct isl_tab *tab, __isl_take isl_basic_map *bmap)
4138 : {
4139 1383503511 : bmap = isl_basic_map_cow(bmap);
4140 1383503511 : if (!tab || !bmap)
4141 : goto error;
4142 :
4143 1383503511 : if (tab->empty) {
4144 1323103687 : bmap = isl_basic_map_set_to_empty(bmap);
4145 1323103687 : if (!bmap)
4146 0 : goto error;
4147 1323103687 : tab->bmap = bmap;
4148 1323103687 : return isl_stat_ok;
4149 : }
4150 :
4151 60399824 : isl_assert(tab->mat->ctx, tab->n_eq == bmap->n_eq, goto error);
4152 60399824 : isl_assert(tab->mat->ctx,
4153 : tab->n_con == bmap->n_eq + bmap->n_ineq, goto error);
4154 :
4155 60399824 : tab->bmap = bmap;
4156 :
4157 60399824 : 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 1105725 : __isl_keep isl_basic_set *isl_tab_peek_bset(struct isl_tab *tab)
4169 : {
4170 1105725 : if (!tab)
4171 0 : return NULL;
4172 :
4173 1105725 : 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 : }
|
