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 7221391 : 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 7221391 : unsigned off = 2 + M;
39 :
40 7221391 : tab = isl_calloc_type(ctx, struct isl_tab);
41 7221391 : if (!tab)
42 0 : return NULL;
43 7221391 : tab->mat = isl_mat_alloc(ctx, n_row, off + n_var);
44 7221391 : if (!tab->mat)
45 0 : goto error;
46 7221391 : tab->var = isl_alloc_array(ctx, struct isl_tab_var, n_var);
47 7221391 : if (n_var && !tab->var)
48 0 : goto error;
49 7221391 : tab->con = isl_alloc_array(ctx, struct isl_tab_var, n_row);
50 7221391 : if (n_row && !tab->con)
51 0 : goto error;
52 7221391 : tab->col_var = isl_alloc_array(ctx, int, n_var);
53 7221391 : if (n_var && !tab->col_var)
54 0 : goto error;
55 7221391 : tab->row_var = isl_alloc_array(ctx, int, n_row);
56 7221391 : if (n_row && !tab->row_var)
57 0 : goto error;
58 44316938 : for (i = 0; i < n_var; ++i) {
59 37095547 : tab->var[i].index = i;
60 37095547 : tab->var[i].is_row = 0;
61 37095547 : tab->var[i].is_nonneg = 0;
62 37095547 : tab->var[i].is_zero = 0;
63 37095547 : tab->var[i].is_redundant = 0;
64 37095547 : tab->var[i].frozen = 0;
65 37095547 : tab->var[i].negated = 0;
66 37095547 : tab->col_var[i] = i;
67 : }
68 7221391 : tab->n_row = 0;
69 7221391 : tab->n_con = 0;
70 7221391 : tab->n_eq = 0;
71 7221391 : tab->max_con = n_row;
72 7221391 : tab->n_col = n_var;
73 7221391 : tab->n_var = n_var;
74 7221391 : tab->max_var = n_var;
75 7221391 : tab->n_param = 0;
76 7221391 : tab->n_div = 0;
77 7221391 : tab->n_dead = 0;
78 7221391 : tab->n_redundant = 0;
79 7221391 : tab->strict_redundant = 0;
80 7221391 : tab->need_undo = 0;
81 7221391 : tab->rational = 0;
82 7221391 : tab->empty = 0;
83 7221391 : tab->in_undo = 0;
84 7221391 : tab->M = M;
85 7221391 : tab->cone = 0;
86 7221391 : tab->bottom.type = isl_tab_undo_bottom;
87 7221391 : tab->bottom.next = NULL;
88 7221391 : tab->top = &tab->bottom;
89 :
90 7221391 : tab->n_zero = 0;
91 7221391 : tab->n_unbounded = 0;
92 7221391 : tab->basis = NULL;
93 :
94 7221391 : return tab;
95 : error:
96 0 : isl_tab_free(tab);
97 0 : return NULL;
98 : }
99 :
100 58195569 : isl_ctx *isl_tab_get_ctx(struct isl_tab *tab)
101 : {
102 58195569 : return tab ? isl_mat_get_ctx(tab->mat) : NULL;
103 : }
104 :
105 7674611 : int isl_tab_extend_cons(struct isl_tab *tab, unsigned n_new)
106 : {
107 : unsigned off;
108 :
109 7674611 : if (!tab)
110 0 : return -1;
111 :
112 7674611 : off = 2 + tab->M;
113 :
114 7674611 : if (tab->max_con < tab->n_con + n_new) {
115 : struct isl_tab_var *con;
116 :
117 122674 : con = isl_realloc_array(tab->mat->ctx, tab->con,
118 : struct isl_tab_var, tab->max_con + n_new);
119 122674 : if (!con)
120 0 : return -1;
121 122674 : tab->con = con;
122 122674 : tab->max_con += n_new;
123 : }
124 7674611 : if (tab->mat->n_row < tab->n_row + n_new) {
125 : int *row_var;
126 :
127 273436 : tab->mat = isl_mat_extend(tab->mat,
128 273436 : tab->n_row + n_new, off + tab->n_col);
129 136718 : if (!tab->mat)
130 0 : return -1;
131 136718 : row_var = isl_realloc_array(tab->mat->ctx, tab->row_var,
132 : int, tab->mat->n_row);
133 136718 : if (!row_var)
134 0 : return -1;
135 136718 : tab->row_var = row_var;
136 136718 : if (tab->row_sign) {
137 : enum isl_tab_row_sign *s;
138 5542 : s = isl_realloc_array(tab->mat->ctx, tab->row_sign,
139 : enum isl_tab_row_sign, tab->mat->n_row);
140 5542 : if (!s)
141 0 : return -1;
142 5542 : tab->row_sign = s;
143 : }
144 : }
145 7674611 : return 0;
146 : }
147 :
148 : /* Make room for at least n_new extra variables.
149 : * Return -1 if anything went wrong.
150 : */
151 498930 : int isl_tab_extend_vars(struct isl_tab *tab, unsigned n_new)
152 : {
153 : struct isl_tab_var *var;
154 498930 : unsigned off = 2 + tab->M;
155 :
156 498930 : if (tab->max_var < tab->n_var + n_new) {
157 48662 : var = isl_realloc_array(tab->mat->ctx, tab->var,
158 : struct isl_tab_var, tab->n_var + n_new);
159 48662 : if (!var)
160 0 : return -1;
161 48662 : tab->var = var;
162 48662 : tab->max_var = tab->n_var + n_new;
163 : }
164 :
165 498930 : if (tab->mat->n_col < off + tab->n_col + n_new) {
166 : int *p;
167 :
168 90600 : tab->mat = isl_mat_extend(tab->mat,
169 90600 : tab->mat->n_row, off + tab->n_col + n_new);
170 45300 : if (!tab->mat)
171 0 : return -1;
172 45300 : p = isl_realloc_array(tab->mat->ctx, tab->col_var,
173 : int, tab->n_col + n_new);
174 45300 : if (!p)
175 0 : return -1;
176 45300 : tab->col_var = p;
177 : }
178 :
179 498930 : return 0;
180 : }
181 :
182 121034131 : static void free_undo_record(struct isl_tab_undo *undo)
183 : {
184 121034131 : switch (undo->type) {
185 : case isl_tab_undo_saved_basis:
186 0 : free(undo->u.col_var);
187 0 : break;
188 : default:;
189 : }
190 121034131 : free(undo);
191 121034131 : }
192 :
193 7359789 : static void free_undo(struct isl_tab *tab)
194 : {
195 : struct isl_tab_undo *undo, *next;
196 :
197 8399193 : for (undo = tab->top; undo && undo != &tab->bottom; undo = next) {
198 1039404 : next = undo->next;
199 1039404 : free_undo_record(undo);
200 : }
201 7359789 : tab->top = undo;
202 7359789 : }
203 :
204 7385285 : void isl_tab_free(struct isl_tab *tab)
205 : {
206 7385285 : if (!tab)
207 26642 : return;
208 7358643 : free_undo(tab);
209 7358643 : isl_mat_free(tab->mat);
210 7358643 : isl_vec_free(tab->dual);
211 7358643 : isl_basic_map_free(tab->bmap);
212 7358643 : free(tab->var);
213 7358643 : free(tab->con);
214 7358643 : free(tab->row_var);
215 7358643 : free(tab->col_var);
216 7358643 : free(tab->row_sign);
217 7358643 : isl_mat_free(tab->samples);
218 7358643 : free(tab->sample_index);
219 7358643 : isl_mat_free(tab->basis);
220 7358643 : free(tab);
221 : }
222 :
223 136010 : struct isl_tab *isl_tab_dup(struct isl_tab *tab)
224 : {
225 : int i;
226 : struct isl_tab *dup;
227 : unsigned off;
228 :
229 136010 : if (!tab)
230 0 : return NULL;
231 :
232 136010 : off = 2 + tab->M;
233 136010 : dup = isl_calloc_type(tab->mat->ctx, struct isl_tab);
234 136010 : if (!dup)
235 0 : return NULL;
236 136010 : dup->mat = isl_mat_dup(tab->mat);
237 136010 : if (!dup->mat)
238 0 : goto error;
239 136010 : dup->var = isl_alloc_array(tab->mat->ctx, struct isl_tab_var, tab->max_var);
240 136010 : if (tab->max_var && !dup->var)
241 0 : goto error;
242 992184 : for (i = 0; i < tab->n_var; ++i)
243 856174 : dup->var[i] = tab->var[i];
244 136010 : dup->con = isl_alloc_array(tab->mat->ctx, struct isl_tab_var, tab->max_con);
245 136010 : if (tab->max_con && !dup->con)
246 0 : goto error;
247 1162093 : for (i = 0; i < tab->n_con; ++i)
248 1026083 : dup->con[i] = tab->con[i];
249 136010 : dup->col_var = isl_alloc_array(tab->mat->ctx, int, tab->mat->n_col - off);
250 136010 : if ((tab->mat->n_col - off) && !dup->col_var)
251 0 : goto error;
252 992184 : for (i = 0; i < tab->n_col; ++i)
253 856174 : dup->col_var[i] = tab->col_var[i];
254 136010 : dup->row_var = isl_alloc_array(tab->mat->ctx, int, tab->mat->n_row);
255 136010 : if (tab->mat->n_row && !dup->row_var)
256 0 : goto error;
257 1162093 : for (i = 0; i < tab->n_row; ++i)
258 1026083 : dup->row_var[i] = tab->row_var[i];
259 136010 : if (tab->row_sign) {
260 136010 : dup->row_sign = isl_alloc_array(tab->mat->ctx, enum isl_tab_row_sign,
261 : tab->mat->n_row);
262 136010 : if (tab->mat->n_row && !dup->row_sign)
263 0 : goto error;
264 1162093 : for (i = 0; i < tab->n_row; ++i)
265 1026083 : dup->row_sign[i] = tab->row_sign[i];
266 : }
267 136010 : 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 136010 : dup->n_row = tab->n_row;
279 136010 : dup->n_con = tab->n_con;
280 136010 : dup->n_eq = tab->n_eq;
281 136010 : dup->max_con = tab->max_con;
282 136010 : dup->n_col = tab->n_col;
283 136010 : dup->n_var = tab->n_var;
284 136010 : dup->max_var = tab->max_var;
285 136010 : dup->n_param = tab->n_param;
286 136010 : dup->n_div = tab->n_div;
287 136010 : dup->n_dead = tab->n_dead;
288 136010 : dup->n_redundant = tab->n_redundant;
289 136010 : dup->rational = tab->rational;
290 136010 : dup->empty = tab->empty;
291 136010 : dup->strict_redundant = 0;
292 136010 : dup->need_undo = 0;
293 136010 : dup->in_undo = 0;
294 136010 : dup->M = tab->M;
295 136010 : tab->cone = tab->cone;
296 136010 : dup->bottom.type = isl_tab_undo_bottom;
297 136010 : dup->bottom.next = NULL;
298 136010 : dup->top = &dup->bottom;
299 :
300 136010 : dup->n_zero = tab->n_zero;
301 136010 : dup->n_unbounded = tab->n_unbounded;
302 136010 : dup->basis = isl_mat_dup(tab->basis);
303 :
304 136010 : 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 1242 : 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 1242 : prod = isl_mat_alloc(mat1->ctx, mat1->n_row + mat2->n_row,
333 1242 : off + col1 + col2);
334 1242 : if (!prod)
335 0 : return NULL;
336 :
337 1242 : n = 0;
338 12732 : for (i = 0; i < r1; ++i) {
339 11490 : isl_seq_cpy(prod->row[n + i], mat1->row[i], off + d1);
340 11490 : isl_seq_clr(prod->row[n + i] + off + d1, d2);
341 22980 : isl_seq_cpy(prod->row[n + i] + off + d1 + d2,
342 11490 : mat1->row[i] + off + d1, col1 - d1);
343 11490 : isl_seq_clr(prod->row[n + i] + off + col1 + d1, col2 - d2);
344 : }
345 :
346 1242 : n += r1;
347 12732 : for (i = 0; i < r2; ++i) {
348 11490 : isl_seq_cpy(prod->row[n + i], mat2->row[i], off);
349 11490 : isl_seq_clr(prod->row[n + i] + off, d1);
350 22980 : isl_seq_cpy(prod->row[n + i] + off + d1,
351 22980 : mat2->row[i] + off, d2);
352 11490 : isl_seq_clr(prod->row[n + i] + off + d1 + d2, col1 - d1);
353 22980 : isl_seq_cpy(prod->row[n + i] + off + col1 + d1,
354 11490 : mat2->row[i] + off + d2, col2 - d2);
355 : }
356 :
357 1242 : n += r2;
358 12988 : for (i = 0; i < row1 - r1; ++i) {
359 11746 : isl_seq_cpy(prod->row[n + i], mat1->row[r1 + i], off + d1);
360 11746 : isl_seq_clr(prod->row[n + i] + off + d1, d2);
361 23492 : isl_seq_cpy(prod->row[n + i] + off + d1 + d2,
362 11746 : mat1->row[r1 + i] + off + d1, col1 - d1);
363 11746 : isl_seq_clr(prod->row[n + i] + off + col1 + d1, col2 - d2);
364 : }
365 :
366 1242 : n += row1 - r1;
367 12988 : for (i = 0; i < row2 - r2; ++i) {
368 11746 : isl_seq_cpy(prod->row[n + i], mat2->row[r2 + i], off);
369 11746 : isl_seq_clr(prod->row[n + i] + off, d1);
370 23492 : isl_seq_cpy(prod->row[n + i] + off + d1,
371 23492 : mat2->row[r2 + i] + off, d2);
372 11746 : isl_seq_clr(prod->row[n + i] + off + d1 + d2, col1 - d1);
373 23492 : isl_seq_cpy(prod->row[n + i] + off + col1 + d1,
374 11746 : mat2->row[r2 + i] + off + d2, col2 - d2);
375 : }
376 :
377 1242 : 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 30680 : static void update_index1(struct isl_tab_var *var,
384 : unsigned r1, unsigned r2, unsigned d1, unsigned d2)
385 : {
386 30680 : if (var->index == -1)
387 32 : return;
388 30648 : if (var->is_row && var->index >= r1)
389 11746 : var->index += r2;
390 30648 : if (!var->is_row && var->index >= d1)
391 6936 : 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 30680 : 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 30680 : if (var->index == -1)
402 32 : return;
403 30648 : if (var->is_row) {
404 23236 : if (var->index < r2)
405 11490 : var->index += r1;
406 : else
407 11746 : var->index += row1;
408 : } else {
409 7412 : if (var->index < d2)
410 476 : var->index += d1;
411 : else
412 6936 : 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 1242 : 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 1242 : if (!tab1 || !tab2)
442 0 : return NULL;
443 :
444 1242 : isl_assert(tab1->mat->ctx, tab1->M == tab2->M, return NULL);
445 1242 : isl_assert(tab1->mat->ctx, tab1->rational == tab2->rational, return NULL);
446 1242 : isl_assert(tab1->mat->ctx, tab1->cone == tab2->cone, return NULL);
447 1242 : isl_assert(tab1->mat->ctx, !tab1->row_sign, return NULL);
448 1242 : isl_assert(tab1->mat->ctx, !tab2->row_sign, return NULL);
449 1242 : isl_assert(tab1->mat->ctx, tab1->n_param == 0, return NULL);
450 1242 : isl_assert(tab1->mat->ctx, tab2->n_param == 0, return NULL);
451 1242 : isl_assert(tab1->mat->ctx, tab1->n_div == 0, return NULL);
452 1242 : isl_assert(tab1->mat->ctx, tab2->n_div == 0, return NULL);
453 :
454 1242 : off = 2 + tab1->M;
455 1242 : r1 = tab1->n_redundant;
456 1242 : r2 = tab2->n_redundant;
457 1242 : d1 = tab1->n_dead;
458 1242 : d2 = tab2->n_dead;
459 1242 : prod = isl_calloc_type(tab1->mat->ctx, struct isl_tab);
460 1242 : if (!prod)
461 0 : return NULL;
462 1242 : 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 1242 : if (!prod->mat)
466 0 : goto error;
467 1242 : prod->var = isl_alloc_array(tab1->mat->ctx, struct isl_tab_var,
468 : tab1->max_var + tab2->max_var);
469 1242 : if ((tab1->max_var + tab2->max_var) && !prod->var)
470 0 : goto error;
471 8686 : for (i = 0; i < tab1->n_var; ++i) {
472 7444 : prod->var[i] = tab1->var[i];
473 7444 : update_index1(&prod->var[i], r1, r2, d1, d2);
474 : }
475 8686 : for (i = 0; i < tab2->n_var; ++i) {
476 7444 : prod->var[tab1->n_var + i] = tab2->var[i];
477 7444 : update_index2(&prod->var[tab1->n_var + i],
478 : tab1->n_row, tab1->n_col,
479 : r1, r2, d1, d2);
480 : }
481 1242 : prod->con = isl_alloc_array(tab1->mat->ctx, struct isl_tab_var,
482 : tab1->max_con + tab2->max_con);
483 1242 : if ((tab1->max_con + tab2->max_con) && !prod->con)
484 0 : goto error;
485 24478 : for (i = 0; i < tab1->n_con; ++i) {
486 23236 : prod->con[i] = tab1->con[i];
487 23236 : update_index1(&prod->con[i], r1, r2, d1, d2);
488 : }
489 24478 : for (i = 0; i < tab2->n_con; ++i) {
490 23236 : prod->con[tab1->n_con + i] = tab2->con[i];
491 23236 : update_index2(&prod->con[tab1->n_con + i],
492 : tab1->n_row, tab1->n_col,
493 : r1, r2, d1, d2);
494 : }
495 1242 : prod->col_var = isl_alloc_array(tab1->mat->ctx, int,
496 : tab1->n_col + tab2->n_col);
497 1242 : if ((tab1->n_col + tab2->n_col) && !prod->col_var)
498 0 : goto error;
499 8654 : for (i = 0; i < tab1->n_col; ++i) {
500 7412 : int pos = i < d1 ? i : i + d2;
501 7412 : prod->col_var[pos] = tab1->col_var[i];
502 : }
503 8654 : for (i = 0; i < tab2->n_col; ++i) {
504 7412 : int pos = i < d2 ? d1 + i : tab1->n_col + i;
505 7412 : int t = tab2->col_var[i];
506 7412 : if (t >= 0)
507 84 : t += tab1->n_var;
508 : else
509 7328 : t -= tab1->n_con;
510 7412 : prod->col_var[pos] = t;
511 : }
512 1242 : prod->row_var = isl_alloc_array(tab1->mat->ctx, int,
513 : tab1->mat->n_row + tab2->mat->n_row);
514 1242 : if ((tab1->mat->n_row + tab2->mat->n_row) && !prod->row_var)
515 0 : goto error;
516 24478 : for (i = 0; i < tab1->n_row; ++i) {
517 23236 : int pos = i < r1 ? i : i + r2;
518 23236 : prod->row_var[pos] = tab1->row_var[i];
519 : }
520 24478 : for (i = 0; i < tab2->n_row; ++i) {
521 23236 : int pos = i < r2 ? r1 + i : tab1->n_row + i;
522 23236 : int t = tab2->row_var[i];
523 23236 : if (t >= 0)
524 7360 : t += tab1->n_var;
525 : else
526 15876 : t -= tab1->n_con;
527 23236 : prod->row_var[pos] = t;
528 : }
529 1242 : prod->samples = NULL;
530 1242 : prod->sample_index = NULL;
531 1242 : prod->n_row = tab1->n_row + tab2->n_row;
532 1242 : prod->n_con = tab1->n_con + tab2->n_con;
533 1242 : prod->n_eq = 0;
534 1242 : prod->max_con = tab1->max_con + tab2->max_con;
535 1242 : prod->n_col = tab1->n_col + tab2->n_col;
536 1242 : prod->n_var = tab1->n_var + tab2->n_var;
537 1242 : prod->max_var = tab1->max_var + tab2->max_var;
538 1242 : prod->n_param = 0;
539 1242 : prod->n_div = 0;
540 1242 : prod->n_dead = tab1->n_dead + tab2->n_dead;
541 1242 : prod->n_redundant = tab1->n_redundant + tab2->n_redundant;
542 1242 : prod->rational = tab1->rational;
543 1242 : prod->empty = tab1->empty || tab2->empty;
544 1242 : prod->strict_redundant = tab1->strict_redundant || tab2->strict_redundant;
545 1242 : prod->need_undo = 0;
546 1242 : prod->in_undo = 0;
547 1242 : prod->M = tab1->M;
548 1242 : prod->cone = tab1->cone;
549 1242 : prod->bottom.type = isl_tab_undo_bottom;
550 1242 : prod->bottom.next = NULL;
551 1242 : prod->top = &prod->bottom;
552 :
553 1242 : prod->n_zero = 0;
554 1242 : prod->n_unbounded = 0;
555 1242 : prod->basis = NULL;
556 :
557 1242 : return prod;
558 : error:
559 0 : isl_tab_free(prod);
560 0 : return NULL;
561 : }
562 :
563 2004410474 : static struct isl_tab_var *var_from_index(struct isl_tab *tab, int i)
564 : {
565 2004410474 : if (i >= 0)
566 511899261 : return &tab->var[i];
567 : else
568 1492511213 : return &tab->con[~i];
569 : }
570 :
571 1382830212 : struct isl_tab_var *isl_tab_var_from_row(struct isl_tab *tab, int i)
572 : {
573 1382830212 : return var_from_index(tab, tab->row_var[i]);
574 : }
575 :
576 526218978 : static struct isl_tab_var *var_from_col(struct isl_tab *tab, int i)
577 : {
578 526218978 : 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 32209952 : static int max_is_manifestly_unbounded(struct isl_tab *tab,
586 : struct isl_tab_var *var)
587 : {
588 : int i;
589 32209952 : unsigned off = 2 + tab->M;
590 :
591 32209952 : if (var->is_row)
592 15546733 : return 0;
593 113817787 : for (i = tab->n_redundant; i < tab->n_row; ++i) {
594 100387415 : if (!isl_int_is_neg(tab->mat->row[i][off + var->index]))
595 80247164 : continue;
596 20140251 : if (isl_tab_var_from_row(tab, i)->is_nonneg)
597 3232847 : return 0;
598 : }
599 13430372 : 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 8459753 : static int min_is_manifestly_unbounded(struct isl_tab *tab,
607 : struct isl_tab_var *var)
608 : {
609 : int i;
610 8459753 : unsigned off = 2 + tab->M;
611 :
612 8459753 : if (var->is_row)
613 1078124 : return 0;
614 34296190 : for (i = tab->n_redundant; i < tab->n_row; ++i) {
615 32275048 : if (!isl_int_is_pos(tab->mat->row[i][off + var->index]))
616 20183640 : continue;
617 12091408 : if (isl_tab_var_from_row(tab, i)->is_nonneg)
618 5360487 : return 0;
619 : }
620 2021142 : return 1;
621 : }
622 :
623 46226555 : static int row_cmp(struct isl_tab *tab, int r1, int r2, int c, isl_int *t)
624 : {
625 46226555 : unsigned off = 2 + tab->M;
626 :
627 46226555 : 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 46226555 : isl_int_mul(*t, tab->mat->row[r1][1], tab->mat->row[r2][off + c]);
636 46226555 : isl_int_submul(*t, tab->mat->row[r2][1], tab->mat->row[r1][off + c]);
637 46226555 : 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 56963456 : 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 56963456 : unsigned off = 2 + tab->M;
666 :
667 56963456 : isl_int_init(t);
668 56963456 : r = -1;
669 483224449 : for (j = tab->n_redundant; j < tab->n_row; ++j) {
670 426260993 : if (var && j == var->index)
671 51109170 : continue;
672 375151823 : if (!isl_tab_var_from_row(tab, j)->is_nonneg)
673 142320049 : continue;
674 232831774 : if (sgn * isl_int_sgn(tab->mat->row[j][off + c]) >= 0)
675 151360365 : continue;
676 81471409 : if (r < 0) {
677 35244854 : r = j;
678 35244854 : continue;
679 : }
680 46226555 : tsgn = sgn * row_cmp(tab, r, j, c, &t);
681 83164179 : if (tsgn < 0 || (tsgn == 0 &&
682 36937624 : tab->row_var[j] < tab->row_var[r]))
683 38169185 : r = j;
684 : }
685 56963456 : isl_int_clear(t);
686 56963456 : 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 71776998 : 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 71776998 : *row = *col = -1;
713 :
714 71776998 : isl_assert(tab->mat->ctx, var->is_row, return);
715 71776998 : tr = tab->mat->row[var->index] + 2 + tab->M;
716 :
717 71776998 : c = -1;
718 431992366 : for (j = tab->n_dead; j < tab->n_col; ++j) {
719 360215368 : if (isl_int_is_zero(tr[j]))
720 154463226 : continue;
721 348639864 : if (isl_int_sgn(tr[j]) != sgn &&
722 142887722 : var_from_col(tab, j)->is_nonneg)
723 108414082 : continue;
724 97338060 : if (c < 0 || tab->col_var[j] < tab->col_var[c])
725 57546410 : c = j;
726 : }
727 71776998 : if (c < 0)
728 20636653 : return;
729 :
730 51140345 : sgn *= isl_int_sgn(tr[c]);
731 51140345 : r = pivot_row(tab, skip_var, sgn, c);
732 51140345 : *row = r < 0 ? var->index : r;
733 51140345 : *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 419605706 : int isl_tab_row_is_redundant(struct isl_tab *tab, int row)
743 : {
744 : int i;
745 419605706 : unsigned off = 2 + tab->M;
746 :
747 419605706 : if (tab->row_var[row] < 0 && !isl_tab_var_from_row(tab, row)->is_nonneg)
748 402191 : return 0;
749 :
750 419203515 : if (isl_int_is_neg(tab->mat->row[row][1]))
751 60235182 : return 0;
752 358968333 : if (tab->strict_redundant && isl_int_is_zero(tab->mat->row[row][1]))
753 0 : return 0;
754 358968333 : if (tab->M && isl_int_is_neg(tab->mat->row[row][2]))
755 6438 : return 0;
756 :
757 986099572 : for (i = tab->n_dead; i < tab->n_col; ++i) {
758 956311765 : if (isl_int_is_zero(tab->mat->row[row][off + i]))
759 382426994 : continue;
760 573884771 : if (tab->col_var[i] >= 0)
761 154284075 : return 0;
762 419600696 : if (isl_int_is_neg(tab->mat->row[row][off + i]))
763 174843728 : return 0;
764 244756968 : if (!var_from_col(tab, i)->is_nonneg)
765 46285 : return 0;
766 : }
767 29787807 : return 1;
768 : }
769 :
770 39149816 : 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 39149816 : t = tab->row_var[row1];
776 39149816 : tab->row_var[row1] = tab->row_var[row2];
777 39149816 : tab->row_var[row2] = t;
778 39149816 : isl_tab_var_from_row(tab, row1)->index = row1;
779 39149816 : isl_tab_var_from_row(tab, row2)->index = row2;
780 39149816 : tab->mat = isl_mat_swap_rows(tab->mat, row1, row2);
781 :
782 39149816 : if (!tab->row_sign)
783 39149670 : return;
784 146 : s = tab->row_sign[row1];
785 146 : tab->row_sign[row1] = tab->row_sign[row2];
786 146 : 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 274433895 : 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 274433895 : if (!tab)
805 0 : return isl_stat_error;
806 274433895 : if (!tab->need_undo)
807 153399764 : return isl_stat_ok;
808 :
809 121034131 : undo = isl_alloc_type(tab->mat->ctx, struct isl_tab_undo);
810 121034131 : if (!undo)
811 0 : goto error;
812 121034131 : undo->type = type;
813 121034131 : undo->u = u;
814 121034131 : undo->next = tab->top;
815 121034131 : tab->top = undo;
816 :
817 121034131 : 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 250013033 : 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 250013033 : if (var->is_row)
829 222285481 : u.var_index = tab->row_var[var->index];
830 : else
831 27727552 : u.var_index = tab->col_var[var->index];
832 250013033 : return push_union(tab, type, u);
833 : }
834 :
835 24138149 : isl_stat isl_tab_push(struct isl_tab *tab, enum isl_tab_undo_type type)
836 : {
837 24138149 : union isl_tab_undo_val u = { 0 };
838 24138149 : 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 136010 : 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 136010 : u.callback = callback;
862 136010 : return push_union(tab, isl_tab_undo_callback, u);
863 : }
864 :
865 10693 : struct isl_tab *isl_tab_init_samples(struct isl_tab *tab)
866 : {
867 10693 : if (!tab)
868 0 : return NULL;
869 :
870 10693 : tab->n_sample = 0;
871 10693 : tab->n_outside = 0;
872 10693 : tab->samples = isl_mat_alloc(tab->mat->ctx, 1, 1 + tab->n_var);
873 10693 : if (!tab->samples)
874 0 : goto error;
875 10693 : tab->sample_index = isl_alloc_array(tab->mat->ctx, int, 1);
876 10693 : if (!tab->sample_index)
877 0 : goto error;
878 10693 : return tab;
879 : error:
880 0 : isl_tab_free(tab);
881 0 : return NULL;
882 : }
883 :
884 213522 : int isl_tab_add_sample(struct isl_tab *tab, __isl_take isl_vec *sample)
885 : {
886 213522 : if (!tab || !sample)
887 : goto error;
888 :
889 213522 : if (tab->n_sample + 1 > tab->samples->n_row) {
890 118703 : int *t = isl_realloc_array(tab->mat->ctx,
891 : tab->sample_index, int, tab->n_sample + 1);
892 118703 : if (!t)
893 0 : goto error;
894 118703 : tab->sample_index = t;
895 : }
896 :
897 427044 : tab->samples = isl_mat_extend(tab->samples,
898 427044 : tab->n_sample + 1, tab->samples->n_col);
899 213522 : if (!tab->samples)
900 0 : goto error;
901 :
902 213522 : isl_seq_cpy(tab->samples->row[tab->n_sample], sample->el, sample->size);
903 213522 : isl_vec_free(sample);
904 213522 : tab->sample_index[tab->n_sample] = tab->n_sample;
905 213522 : tab->n_sample++;
906 :
907 213522 : return 0;
908 : error:
909 0 : isl_vec_free(sample);
910 0 : return -1;
911 : }
912 :
913 385819 : struct isl_tab *isl_tab_drop_sample(struct isl_tab *tab, int s)
914 : {
915 385819 : if (s != tab->n_outside) {
916 295848 : int t = tab->sample_index[tab->n_outside];
917 295848 : tab->sample_index[tab->n_outside] = tab->sample_index[s];
918 295848 : tab->sample_index[s] = t;
919 295848 : isl_mat_swap_rows(tab->samples, tab->n_outside, s);
920 : }
921 385819 : tab->n_outside++;
922 385819 : if (isl_tab_push(tab, isl_tab_undo_drop_sample) < 0) {
923 0 : isl_tab_free(tab);
924 0 : return NULL;
925 : }
926 :
927 385819 : return tab;
928 : }
929 :
930 : /* Record the current number of samples so that we can remove newer
931 : * samples during a rollback.
932 : */
933 146703 : isl_stat isl_tab_save_samples(struct isl_tab *tab)
934 : {
935 : union isl_tab_undo_val u;
936 :
937 146703 : if (!tab)
938 0 : return isl_stat_error;
939 :
940 146703 : u.n = tab->n_sample;
941 146703 : 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 41404408 : int isl_tab_mark_redundant(struct isl_tab *tab, int row)
957 : {
958 41404408 : struct isl_tab_var *var = isl_tab_var_from_row(tab, row);
959 41404408 : var->is_redundant = 1;
960 41404408 : isl_assert(tab->mat->ctx, row >= tab->n_redundant, return -1);
961 41404408 : if (tab->preserve || tab->need_undo || tab->row_var[row] >= 0) {
962 23953326 : if (tab->row_var[row] >= 0 && !var->is_nonneg) {
963 15193164 : var->is_nonneg = 1;
964 15193164 : if (isl_tab_push_var(tab, isl_tab_undo_nonneg, var) < 0)
965 0 : return -1;
966 : }
967 23953326 : if (row != tab->n_redundant)
968 18690777 : swap_rows(tab, row, tab->n_redundant);
969 23953326 : tab->n_redundant++;
970 23953326 : return isl_tab_push_var(tab, isl_tab_undo_redundant, var);
971 : } else {
972 17451082 : if (row != tab->n_row - 1)
973 13156220 : swap_rows(tab, row, tab->n_row - 1);
974 17451082 : isl_tab_var_from_row(tab, tab->n_row - 1)->index = -1;
975 17451082 : tab->n_row--;
976 17451082 : 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 1234 : int isl_tab_mark_rational(struct isl_tab *tab)
986 : {
987 1234 : if (!tab)
988 0 : return -1;
989 1234 : if (!tab->rational && tab->need_undo)
990 1234 : if (isl_tab_push(tab, isl_tab_undo_rational) < 0)
991 0 : return -1;
992 1234 : tab->rational = 1;
993 1234 : return 0;
994 : }
995 :
996 8771765 : isl_stat isl_tab_mark_empty(struct isl_tab *tab)
997 : {
998 8771765 : if (!tab)
999 0 : return isl_stat_error;
1000 8771765 : if (!tab->empty && tab->need_undo)
1001 5475217 : if (isl_tab_push(tab, isl_tab_undo_empty) < 0)
1002 0 : return isl_stat_error;
1003 8771765 : tab->empty = 1;
1004 8771765 : return isl_stat_ok;
1005 : }
1006 :
1007 50495221 : int isl_tab_freeze_constraint(struct isl_tab *tab, int con)
1008 : {
1009 : struct isl_tab_var *var;
1010 :
1011 50495221 : if (!tab)
1012 0 : return -1;
1013 :
1014 50495221 : var = &tab->con[con];
1015 50495221 : if (var->frozen)
1016 0 : return 0;
1017 50495221 : if (var->index < 0)
1018 1089508 : return 0;
1019 49405713 : var->frozen = 1;
1020 :
1021 49405713 : if (tab->need_undo)
1022 49405713 : 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 58195569 : static void update_row_sign(struct isl_tab *tab, int row, int col, int row_sgn)
1045 : {
1046 : int i;
1047 58195569 : struct isl_mat *mat = tab->mat;
1048 58195569 : unsigned off = 2 + tab->M;
1049 :
1050 58195569 : if (!tab->row_sign)
1051 58121023 : return;
1052 :
1053 74546 : if (tab->row_sign[row] == 0)
1054 1458 : return;
1055 73088 : isl_assert(mat->ctx, row_sgn > 0, return);
1056 73088 : isl_assert(mat->ctx, tab->row_sign[row] == isl_tab_row_neg, return);
1057 73088 : tab->row_sign[row] = isl_tab_row_pos;
1058 605953 : for (i = 0; i < tab->n_row; ++i) {
1059 : int s;
1060 532865 : if (i == row)
1061 73088 : continue;
1062 459777 : s = isl_int_sgn(mat->row[i][off + col]);
1063 459777 : if (!s)
1064 147439 : continue;
1065 312338 : if (!tab->row_sign[i])
1066 43406 : continue;
1067 268932 : if (s < 0 && tab->row_sign[i] == isl_tab_row_neg)
1068 0 : continue;
1069 268932 : if (s > 0 && tab->row_sign[i] == isl_tab_row_pos)
1070 152997 : continue;
1071 115935 : 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 58195569 : 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 58195569 : struct isl_mat *mat = tab->mat;
1132 : struct isl_tab_var *var;
1133 58195569 : unsigned off = 2 + tab->M;
1134 :
1135 58195569 : ctx = isl_tab_get_ctx(tab);
1136 58195569 : if (isl_ctx_next_operation(ctx) < 0)
1137 0 : return -1;
1138 :
1139 58195569 : isl_int_swap(mat->row[row][0], mat->row[row][off + col]);
1140 58195569 : sgn = isl_int_sgn(mat->row[row][0]);
1141 58195569 : if (sgn < 0) {
1142 29626259 : isl_int_neg(mat->row[row][0], mat->row[row][0]);
1143 29626259 : isl_int_neg(mat->row[row][off + col], mat->row[row][off + col]);
1144 : } else
1145 213411493 : for (j = 0; j < off - 1 + tab->n_col; ++j) {
1146 184842183 : if (j == off - 1 + col)
1147 28569310 : continue;
1148 156272873 : isl_int_neg(mat->row[row][1 + j], mat->row[row][1 + j]);
1149 : }
1150 58195569 : if (!isl_int_is_one(mat->row[row][0]))
1151 19533935 : isl_seq_normalize(mat->ctx, mat->row[row], off + tab->n_col);
1152 547356823 : for (i = 0; i < tab->n_row; ++i) {
1153 489161254 : if (i == row)
1154 58195569 : continue;
1155 430965685 : if (isl_int_is_zero(mat->row[i][off + col]))
1156 152025924 : continue;
1157 278939761 : isl_int_mul(mat->row[i][0], mat->row[i][0], mat->row[row][0]);
1158 2232067275 : for (j = 0; j < off - 1 + tab->n_col; ++j) {
1159 1953127514 : if (j == off - 1 + col)
1160 278939761 : continue;
1161 1674187753 : isl_int_mul(mat->row[i][1 + j],
1162 : mat->row[i][1 + j], mat->row[row][0]);
1163 1674187753 : isl_int_addmul(mat->row[i][1 + j],
1164 : mat->row[i][off + col], mat->row[row][1 + j]);
1165 : }
1166 278939761 : isl_int_mul(mat->row[i][off + col],
1167 : mat->row[i][off + col], mat->row[row][off + col]);
1168 278939761 : if (!isl_int_is_one(mat->row[i][0]))
1169 172502236 : isl_seq_normalize(mat->ctx, mat->row[i], off + tab->n_col);
1170 : }
1171 58195569 : t = tab->row_var[row];
1172 58195569 : tab->row_var[row] = tab->col_var[col];
1173 58195569 : tab->col_var[col] = t;
1174 58195569 : var = isl_tab_var_from_row(tab, row);
1175 58195569 : var->is_row = 1;
1176 58195569 : var->index = row;
1177 58195569 : var = var_from_col(tab, col);
1178 58195569 : var->is_row = 0;
1179 58195569 : var->index = col;
1180 58195569 : update_row_sign(tab, row, col, sgn);
1181 58195569 : if (tab->in_undo)
1182 3465228 : return 0;
1183 463150581 : for (i = tab->n_redundant; i < tab->n_row; ++i) {
1184 408420240 : if (isl_int_is_zero(mat->row[i][off + col]))
1185 108097722 : continue;
1186 590996735 : if (!isl_tab_var_from_row(tab, i)->frozen &&
1187 290674217 : isl_tab_row_is_redundant(tab, i)) {
1188 23138002 : int redo = isl_tab_mark_redundant(tab, i);
1189 23138002 : if (redo < 0)
1190 0 : return -1;
1191 23138002 : if (redo)
1192 5541434 : --i;
1193 : }
1194 : }
1195 54730341 : 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 21287503 : static int to_row(struct isl_tab *tab, struct isl_tab_var *var, int sign)
1206 : {
1207 : int r;
1208 21287503 : unsigned off = 2 + tab->M;
1209 :
1210 21287503 : if (var->is_row)
1211 15552019 : return 0;
1212 :
1213 5735484 : if (sign == 0) {
1214 3418208 : for (r = tab->n_redundant; r < tab->n_row; ++r)
1215 3418208 : if (!isl_int_is_zero(tab->mat->row[r][off+var->index]))
1216 : break;
1217 951166 : isl_assert(tab->mat->ctx, r < tab->n_row, return -1);
1218 : } else {
1219 4784318 : r = pivot_row(tab, NULL, sign, var->index);
1220 4784318 : isl_assert(tab->mat->ctx, r >= 0, return -1);
1221 : }
1222 :
1223 5735484 : 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 20183035 : static int sign_of_max(struct isl_tab *tab, struct isl_tab_var *var)
1265 : {
1266 : int row, col;
1267 :
1268 20183035 : if (max_is_manifestly_unbounded(tab, var))
1269 2681605 : return 1;
1270 17501430 : if (to_row(tab, var, 1) < 0)
1271 0 : return -2;
1272 53961064 : while (!isl_int_is_pos(tab->mat->row[var->index][1])) {
1273 36137289 : find_pivot(tab, var, var, 1, &row, &col);
1274 36137289 : if (row == -1)
1275 11614345 : return isl_int_sgn(tab->mat->row[var->index][1]);
1276 24522944 : if (isl_tab_pivot(tab, row, col) < 0)
1277 0 : return -2;
1278 24522944 : if (!var->is_row) /* manifestly unbounded */
1279 5564740 : return 1;
1280 : }
1281 322345 : 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 91788910 : static int row_is_neg(struct isl_tab *tab, int row)
1299 : {
1300 91788910 : if (!tab->M)
1301 91788910 : 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 67147886 : static int row_sgn(struct isl_tab *tab, int row)
1310 : {
1311 67147886 : if (!tab->M)
1312 67147886 : 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 82443238 : static int restore_row(struct isl_tab *tab, struct isl_tab_var *var)
1325 : {
1326 : int row, col;
1327 :
1328 174232148 : while (row_is_neg(tab, var->index)) {
1329 33304877 : find_pivot(tab, var, var, 1, &row, &col);
1330 33304877 : if (row == -1)
1331 8663853 : break;
1332 24641024 : if (isl_tab_pivot(tab, row, col) < 0)
1333 0 : return -2;
1334 24641024 : if (!var->is_row) /* manifestly unbounded */
1335 15295352 : return 1;
1336 : }
1337 67147886 : 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 64231 : static int at_least_zero(struct isl_tab *tab, struct isl_tab_var *var)
1346 : {
1347 : int row, col;
1348 :
1349 148738 : while (isl_int_is_neg(tab->mat->row[var->index][1])) {
1350 79267 : find_pivot(tab, var, var, 1, &row, &col);
1351 79267 : if (row == -1)
1352 41132 : break;
1353 38135 : if (row == var->index) /* manifestly unbounded */
1354 17859 : return 1;
1355 20276 : if (isl_tab_pivot(tab, row, col) < 0)
1356 0 : return -1;
1357 : }
1358 46372 : 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 535 : static int sign_of_min(struct isl_tab *tab, struct isl_tab_var *var)
1379 : {
1380 : int row, col;
1381 535 : struct isl_tab_var *pivot_var = NULL;
1382 :
1383 535 : if (min_is_manifestly_unbounded(tab, var))
1384 0 : return -1;
1385 535 : if (!var->is_row) {
1386 56 : col = var->index;
1387 56 : row = pivot_row(tab, NULL, -1, col);
1388 56 : pivot_var = var_from_col(tab, col);
1389 56 : if (isl_tab_pivot(tab, row, col) < 0)
1390 0 : return -2;
1391 56 : if (var->is_redundant)
1392 7 : return 0;
1393 49 : if (isl_int_is_neg(tab->mat->row[var->index][1])) {
1394 26 : if (var->is_nonneg) {
1395 52 : if (!pivot_var->is_redundant &&
1396 26 : pivot_var->index == row) {
1397 52 : if (isl_tab_pivot(tab, row, col) < 0)
1398 0 : return -2;
1399 : } else
1400 0 : if (restore_row(tab, var) < -1)
1401 0 : return -2;
1402 : }
1403 26 : return -1;
1404 : }
1405 : }
1406 502 : if (var->is_redundant)
1407 0 : return 0;
1408 1206 : while (!isl_int_is_neg(tab->mat->row[var->index][1])) {
1409 601 : find_pivot(tab, var, var, -1, &row, &col);
1410 601 : if (row == var->index)
1411 334 : return -1;
1412 267 : if (row == -1)
1413 62 : return isl_int_sgn(tab->mat->row[var->index][1]);
1414 205 : pivot_var = var_from_col(tab, col);
1415 205 : if (isl_tab_pivot(tab, row, col) < 0)
1416 0 : return -2;
1417 205 : if (var->is_redundant)
1418 3 : return 0;
1419 : }
1420 103 : if (pivot_var && var->is_nonneg) {
1421 : /* pivot back to non-negative value */
1422 15 : if (!pivot_var->is_redundant && pivot_var->index == row) {
1423 30 : if (isl_tab_pivot(tab, row, col) < 0)
1424 0 : return -2;
1425 : } else
1426 0 : if (restore_row(tab, var) < -1)
1427 0 : return -2;
1428 : }
1429 103 : return -1;
1430 : }
1431 :
1432 1716599 : static int row_at_most_neg_one(struct isl_tab *tab, int row)
1433 : {
1434 1716599 : 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 2830036 : return isl_int_is_neg(tab->mat->row[row][1]) &&
1441 1113437 : 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 2285077 : 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 2285077 : if (min_is_manifestly_unbounded(tab, var))
1460 173981 : return 1;
1461 2111096 : if (!var->is_row) {
1462 1038737 : col = var->index;
1463 1038737 : row = pivot_row(tab, NULL, -1, col);
1464 1038737 : pivot_var = var_from_col(tab, col);
1465 1038737 : if (isl_tab_pivot(tab, row, col) < 0)
1466 0 : return -1;
1467 1038737 : if (var->is_redundant)
1468 65771 : return 0;
1469 972966 : if (row_at_most_neg_one(tab, var->index)) {
1470 662452 : if (var->is_nonneg) {
1471 1324904 : if (!pivot_var->is_redundant &&
1472 662452 : pivot_var->index == row) {
1473 1311840 : if (isl_tab_pivot(tab, row, col) < 0)
1474 0 : return -1;
1475 : } else
1476 6532 : if (restore_row(tab, var) < -1)
1477 0 : return -1;
1478 : }
1479 662452 : return 1;
1480 : }
1481 : }
1482 1382873 : if (var->is_redundant)
1483 18150 : return 0;
1484 : do {
1485 1845249 : find_pivot(tab, var, var, -1, &row, &col);
1486 1845249 : if (row == var->index) {
1487 776602 : if (var->is_nonneg && restore_row(tab, var) < -1)
1488 0 : return -1;
1489 776602 : return 1;
1490 : }
1491 1068647 : if (row == -1)
1492 130045 : return 0;
1493 938602 : pivot_var = var_from_col(tab, col);
1494 938602 : if (isl_tab_pivot(tab, row, col) < 0)
1495 0 : return -1;
1496 938602 : if (var->is_redundant)
1497 194969 : return 0;
1498 743633 : } while (!row_at_most_neg_one(tab, var->index));
1499 263107 : if (var->is_nonneg) {
1500 : /* pivot back to non-negative value */
1501 255966 : if (!pivot_var->is_redundant && pivot_var->index == row)
1502 253298 : if (isl_tab_pivot(tab, row, col) < 0)
1503 0 : return -1;
1504 255966 : if (restore_row(tab, var) < -1)
1505 0 : return -1;
1506 : }
1507 263107 : return 1;
1508 : }
1509 :
1510 : /* Return 1 if "var" can attain values >= 1.
1511 : * Return 0 otherwise.
1512 : */
1513 8567565 : 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 8567565 : if (max_is_manifestly_unbounded(tab, var))
1519 8246344 : return 1;
1520 321221 : if (to_row(tab, var, 1) < 0)
1521 0 : return -1;
1522 321221 : r = tab->mat->row[var->index];
1523 660386 : while (isl_int_lt(r[1], r[0])) {
1524 108386 : find_pivot(tab, var, var, 1, &row, &col);
1525 108386 : if (row == -1)
1526 30756 : return isl_int_ge(r[1], r[0]);
1527 77630 : if (row == var->index) /* manifestly unbounded */
1528 59686 : return 1;
1529 17944 : if (isl_tab_pivot(tab, row, col) < 0)
1530 0 : return -1;
1531 : }
1532 230779 : return 1;
1533 : }
1534 :
1535 9622839 : static void swap_cols(struct isl_tab *tab, int col1, int col2)
1536 : {
1537 : int t;
1538 9622839 : unsigned off = 2 + tab->M;
1539 9622839 : t = tab->col_var[col1];
1540 9622839 : tab->col_var[col1] = tab->col_var[col2];
1541 9622839 : tab->col_var[col2] = t;
1542 9622839 : var_from_col(tab, col1)->index = col1;
1543 9622839 : var_from_col(tab, col2)->index = col2;
1544 9622839 : tab->mat = isl_mat_swap_cols(tab->mat, off + col1, off + col2);
1545 9622839 : }
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 12634019 : int isl_tab_kill_col(struct isl_tab *tab, int col)
1560 : {
1561 12634019 : var_from_col(tab, col)->is_zero = 1;
1562 12634019 : if (tab->need_undo) {
1563 82716 : if (isl_tab_push_var(tab, isl_tab_undo_zero,
1564 : var_from_col(tab, col)) < 0)
1565 0 : return -1;
1566 82716 : if (col != tab->n_dead)
1567 59905 : swap_cols(tab, col, tab->n_dead);
1568 82716 : tab->n_dead++;
1569 82716 : return 0;
1570 : } else {
1571 12551303 : if (col != tab->n_col - 1)
1572 9457901 : swap_cols(tab, col, tab->n_col - 1);
1573 12551303 : var_from_col(tab, tab->n_col - 1)->index = -1;
1574 12551303 : tab->n_col--;
1575 12551303 : return 1;
1576 : }
1577 : }
1578 :
1579 56449980 : static int row_is_manifestly_non_integral(struct isl_tab *tab, int row)
1580 : {
1581 56449980 : unsigned off = 2 + tab->M;
1582 :
1583 56449980 : if (tab->M && !isl_int_eq(tab->mat->row[row][2],
1584 : tab->mat->row[row][0]))
1585 0 : return 0;
1586 56449980 : if (isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
1587 56449980 : tab->n_col - tab->n_dead) != -1)
1588 10790786 : return 0;
1589 :
1590 45659194 : 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 11612756 : static int tab_is_manifestly_empty(struct isl_tab *tab)
1598 : {
1599 : int i;
1600 :
1601 11612756 : if (tab->empty)
1602 40 : return 1;
1603 11612716 : if (tab->rational)
1604 1251 : return 0;
1605 :
1606 77962016 : for (i = 0; i < tab->n_var; ++i) {
1607 66351182 : if (!tab->var[i].is_row)
1608 9901202 : continue;
1609 56449980 : if (row_is_manifestly_non_integral(tab, tab->var[i].index))
1610 631 : return 1;
1611 : }
1612 :
1613 11610834 : 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 11612756 : static isl_stat close_row(struct isl_tab *tab, struct isl_tab_var *var,
1636 : int temp_var)
1637 : {
1638 : int j;
1639 11612756 : struct isl_mat *mat = tab->mat;
1640 11612756 : unsigned off = 2 + tab->M;
1641 :
1642 11612756 : 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 11612756 : var->is_zero = 1;
1647 11612756 : if (!temp_var && tab->need_undo)
1648 27157 : if (isl_tab_push_var(tab, isl_tab_undo_zero, var) < 0)
1649 0 : return isl_stat_error;
1650 40461176 : for (j = tab->n_dead; j < tab->n_col; ++j) {
1651 : int recheck;
1652 28848420 : if (isl_int_is_zero(mat->row[var->index][off + j]))
1653 16396378 : continue;
1654 12452042 : 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 12452042 : recheck = isl_tab_kill_col(tab, j);
1659 12452042 : if (recheck < 0)
1660 0 : return isl_stat_error;
1661 12452042 : if (recheck)
1662 12420768 : --j;
1663 : }
1664 11612756 : if (!temp_var && isl_tab_mark_redundant(tab, var->index) < 0)
1665 0 : return isl_stat_error;
1666 11612756 : if (tab_is_manifestly_empty(tab) && isl_tab_mark_empty(tab) < 0)
1667 0 : return isl_stat_error;
1668 11612756 : 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 80456859 : int isl_tab_allocate_con(struct isl_tab *tab)
1678 : {
1679 : int r;
1680 :
1681 80456859 : isl_assert(tab->mat->ctx, tab->n_row < tab->mat->n_row, return -1);
1682 80456859 : isl_assert(tab->mat->ctx, tab->n_con < tab->max_con, return -1);
1683 :
1684 80456859 : r = tab->n_con;
1685 80456859 : tab->con[r].index = tab->n_row;
1686 80456859 : tab->con[r].is_row = 1;
1687 80456859 : tab->con[r].is_nonneg = 0;
1688 80456859 : tab->con[r].is_zero = 0;
1689 80456859 : tab->con[r].is_redundant = 0;
1690 80456859 : tab->con[r].frozen = 0;
1691 80456859 : tab->con[r].negated = 0;
1692 80456859 : tab->row_var[tab->n_row] = ~r;
1693 :
1694 80456859 : tab->n_row++;
1695 80456859 : tab->n_con++;
1696 80456859 : if (isl_tab_push_var(tab, isl_tab_undo_allocate, &tab->con[r]) < 0)
1697 0 : return -1;
1698 :
1699 80456859 : 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 498930 : static int var_insert_entry(struct isl_tab *tab, int first)
1708 : {
1709 : int i;
1710 :
1711 498930 : 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 498930 : 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 498930 : 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 498930 : tab->n_var++;
1727 :
1728 498930 : 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 467087 : static int var_drop_entry(struct isl_tab *tab, int first)
1737 : {
1738 : int i;
1739 :
1740 467087 : 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 467087 : tab->n_var--;
1745 :
1746 467087 : 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 467087 : 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 498930 : int isl_tab_insert_var(struct isl_tab *tab, int r)
1762 : {
1763 : int i;
1764 498930 : unsigned off = 2 + tab->M;
1765 :
1766 498930 : isl_assert(tab->mat->ctx, tab->n_col < tab->mat->n_col, return -1);
1767 :
1768 498930 : if (var_insert_entry(tab, r) < 0)
1769 0 : return -1;
1770 :
1771 498930 : tab->var[r].index = tab->n_col;
1772 498930 : tab->var[r].is_row = 0;
1773 498930 : tab->var[r].is_nonneg = 0;
1774 498930 : tab->var[r].is_zero = 0;
1775 498930 : tab->var[r].is_redundant = 0;
1776 498930 : tab->var[r].frozen = 0;
1777 498930 : tab->var[r].negated = 0;
1778 498930 : tab->col_var[tab->n_col] = r;
1779 :
1780 5107422 : for (i = 0; i < tab->n_row; ++i)
1781 4608492 : isl_int_set_si(tab->mat->row[i][off + tab->n_col], 0);
1782 :
1783 498930 : tab->n_col++;
1784 498930 : if (isl_tab_push_var(tab, isl_tab_undo_allocate, &tab->var[r]) < 0)
1785 0 : return -1;
1786 :
1787 498930 : 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 80418930 : 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 80418930 : unsigned off = 2 + tab->M;
1830 :
1831 80418930 : r = isl_tab_allocate_con(tab);
1832 80418930 : if (r < 0)
1833 0 : return -1;
1834 :
1835 80418930 : isl_int_init(a);
1836 80418930 : isl_int_init(b);
1837 80418930 : row = tab->mat->row[tab->con[r].index];
1838 80418930 : isl_int_set_si(row[0], 1);
1839 80418930 : isl_int_set(row[1], line[0]);
1840 80418930 : isl_seq_clr(row + 2, tab->M + tab->n_col);
1841 542300695 : for (i = 0; i < tab->n_var; ++i) {
1842 461881765 : if (tab->var[i].is_zero)
1843 0 : continue;
1844 461881765 : if (tab->var[i].is_row) {
1845 155014593 : isl_int_lcm(a,
1846 : row[0], tab->mat->row[tab->var[i].index][0]);
1847 155014593 : isl_int_swap(a, row[0]);
1848 155014593 : isl_int_divexact(a, row[0], a);
1849 155014593 : isl_int_divexact(b,
1850 : row[0], tab->mat->row[tab->var[i].index][0]);
1851 155014593 : isl_int_mul(b, b, line[1 + i]);
1852 310029186 : isl_seq_combine(row + 1, a, row + 1,
1853 155014593 : b, tab->mat->row[tab->var[i].index] + 1,
1854 155014593 : 1 + tab->M + tab->n_col);
1855 : } else
1856 306867172 : isl_int_addmul(row[off + tab->var[i].index],
1857 : line[1 + i], row[0]);
1858 461881765 : if (tab->M && i >= tab->n_param && i < tab->n_var - tab->n_div)
1859 90824 : isl_int_submul(row[2], line[1 + i], row[0]);
1860 : }
1861 80418930 : isl_seq_normalize(tab->mat->ctx, row, off + tab->n_col);
1862 80418930 : isl_int_clear(a);
1863 80418930 : isl_int_clear(b);
1864 :
1865 80418930 : if (tab->row_sign)
1866 75172 : tab->row_sign[tab->con[r].index] = isl_tab_row_unknown;
1867 :
1868 80418930 : return r;
1869 : }
1870 :
1871 18526086 : static isl_stat drop_row(struct isl_tab *tab, int row)
1872 : {
1873 18526086 : isl_assert(tab->mat->ctx, ~tab->row_var[row] == tab->n_con - 1,
1874 : return isl_stat_error);
1875 18526086 : if (row != tab->n_row - 1)
1876 7302819 : swap_rows(tab, row, tab->n_row - 1);
1877 18526086 : tab->n_row--;
1878 18526086 : tab->n_con--;
1879 18526086 : 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 467087 : static isl_stat drop_col(struct isl_tab *tab, int col)
1889 : {
1890 : int var;
1891 :
1892 467087 : var = tab->col_var[col];
1893 467087 : if (col != tab->n_col - 1)
1894 105033 : swap_cols(tab, col, tab->n_col - 1);
1895 467087 : tab->n_col--;
1896 467087 : if (var_drop_entry(tab, var) < 0)
1897 0 : return isl_stat_error;
1898 467087 : 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 79557146 : 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 79557146 : if (!tab)
1914 0 : return isl_stat_error;
1915 79557146 : if (tab->bmap) {
1916 18263891 : struct isl_basic_map *bmap = tab->bmap;
1917 :
1918 18263891 : isl_assert(tab->mat->ctx, tab->n_eq == bmap->n_eq,
1919 : return isl_stat_error);
1920 18263891 : isl_assert(tab->mat->ctx,
1921 : tab->n_con == bmap->n_eq + bmap->n_ineq,
1922 : return isl_stat_error);
1923 18263891 : tab->bmap = isl_basic_map_add_ineq(tab->bmap, ineq);
1924 18263891 : if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
1925 0 : return isl_stat_error;
1926 18263891 : if (!tab->bmap)
1927 0 : return isl_stat_error;
1928 : }
1929 79557146 : if (tab->cone) {
1930 136356 : isl_int_init(cst);
1931 136356 : isl_int_set_si(cst, 0);
1932 136356 : isl_int_swap(ineq[0], cst);
1933 : }
1934 79557146 : r = isl_tab_add_row(tab, ineq);
1935 79557146 : if (tab->cone) {
1936 136356 : isl_int_swap(ineq[0], cst);
1937 136356 : isl_int_clear(cst);
1938 : }
1939 79557146 : if (r < 0)
1940 0 : return isl_stat_error;
1941 79557146 : tab->con[r].is_nonneg = 1;
1942 79557146 : if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
1943 0 : return isl_stat_error;
1944 79557146 : if (isl_tab_row_is_redundant(tab, tab->con[r].index)) {
1945 6649805 : if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0)
1946 0 : return isl_stat_error;
1947 6649805 : return isl_stat_ok;
1948 : }
1949 :
1950 72907341 : sgn = restore_row(tab, &tab->con[r]);
1951 72907341 : if (sgn < -1)
1952 0 : return isl_stat_error;
1953 72907341 : if (sgn < 0)
1954 8662170 : return isl_tab_mark_empty(tab);
1955 64245171 : 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 64245171 : 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 49984 : static int to_col(struct isl_tab *tab, struct isl_tab_var *var)
1966 : {
1967 : int i;
1968 : int row, col;
1969 49984 : unsigned off = 2 + tab->M;
1970 :
1971 49984 : if (!var->is_row)
1972 1 : return 0;
1973 :
1974 102340 : while (isl_int_is_pos(tab->mat->row[var->index][1])) {
1975 27240 : find_pivot(tab, var, NULL, -1, &row, &col);
1976 27240 : isl_assert(tab->mat->ctx, row != -1, return -1);
1977 27240 : if (isl_tab_pivot(tab, row, col) < 0)
1978 0 : return -1;
1979 27240 : if (!var->is_row)
1980 24866 : return 0;
1981 : }
1982 :
1983 98131 : for (i = tab->n_dead; i < tab->n_col; ++i)
1984 98131 : if (!isl_int_is_zero(tab->mat->row[var->index][off + i]))
1985 : break;
1986 :
1987 25117 : isl_assert(tab->mat->ctx, i < tab->n_col, return -1);
1988 25117 : if (isl_tab_pivot(tab, var->index, i) < 0)
1989 0 : return -1;
1990 :
1991 25117 : 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 130535 : static struct isl_tab *add_eq(struct isl_tab *tab, isl_int *eq)
2003 : {
2004 : int i;
2005 : int r;
2006 :
2007 130535 : if (!tab)
2008 0 : return NULL;
2009 130535 : r = isl_tab_add_row(tab, eq);
2010 130535 : if (r < 0)
2011 0 : goto error;
2012 :
2013 130535 : r = tab->con[r].index;
2014 130535 : i = isl_seq_first_non_zero(tab->mat->row[r] + 2 + tab->M + tab->n_dead,
2015 130535 : tab->n_col - tab->n_dead);
2016 130535 : isl_assert(tab->mat->ctx, i >= 0, goto error);
2017 130535 : i += tab->n_dead;
2018 130535 : if (isl_tab_pivot(tab, r, i) < 0)
2019 0 : goto error;
2020 130535 : if (isl_tab_kill_col(tab, i) < 0)
2021 0 : goto error;
2022 130535 : tab->n_eq++;
2023 :
2024 130535 : 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 1676974 : static int row_is_big(struct isl_tab *tab, int row)
2034 : {
2035 1676974 : return tab->M && !isl_int_is_zero(tab->mat->row[row][2]);
2036 : }
2037 :
2038 52570 : static int row_is_manifestly_zero(struct isl_tab *tab, int row)
2039 : {
2040 52570 : unsigned off = 2 + tab->M;
2041 :
2042 52570 : if (!isl_int_is_zero(tab->mat->row[row][1]))
2043 26387 : return 0;
2044 26183 : if (row_is_big(tab, row))
2045 0 : return 0;
2046 52366 : return isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
2047 52366 : 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 46127 : 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 46127 : if (!tab)
2061 0 : return -1;
2062 46127 : r = isl_tab_add_row(tab, eq);
2063 46127 : if (r < 0)
2064 0 : return -1;
2065 :
2066 46127 : var = &tab->con[r];
2067 46127 : r = var->index;
2068 46127 : if (row_is_manifestly_zero(tab, r)) {
2069 2156 : var->is_zero = 1;
2070 2156 : if (isl_tab_mark_redundant(tab, r) < 0)
2071 0 : return -1;
2072 2156 : return 0;
2073 : }
2074 :
2075 43971 : if (isl_int_is_neg(tab->mat->row[r][1])) {
2076 2152 : isl_seq_neg(tab->mat->row[r] + 1, tab->mat->row[r] + 1,
2077 2152 : 1 + tab->n_col);
2078 2152 : var->negated = 1;
2079 : }
2080 43971 : var->is_nonneg = 1;
2081 43971 : if (to_col(tab, var) < 0)
2082 0 : return -1;
2083 43971 : var->is_nonneg = 0;
2084 43971 : if (isl_tab_kill_col(tab, var->index) < 0)
2085 0 : return -1;
2086 :
2087 43971 : 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 5994 : static int add_zero_row(struct isl_tab *tab)
2097 : {
2098 : int r;
2099 : isl_int *row;
2100 :
2101 5994 : r = isl_tab_allocate_con(tab);
2102 5994 : if (r < 0)
2103 0 : return -1;
2104 :
2105 5994 : row = tab->mat->row[tab->con[r].index];
2106 5994 : isl_seq_clr(row + 1, 1 + tab->M + tab->n_col);
2107 5994 : isl_int_set_si(row[0], 1);
2108 :
2109 5994 : 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 6443 : int isl_tab_add_eq(struct isl_tab *tab, isl_int *eq)
2120 : {
2121 6443 : 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 6443 : if (!tab)
2129 0 : return -1;
2130 6443 : isl_assert(tab->mat->ctx, !tab->M, return -1);
2131 :
2132 6443 : if (tab->need_undo)
2133 6443 : snap = isl_tab_snap(tab);
2134 :
2135 6443 : 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 6443 : r = isl_tab_add_row(tab, eq);
2141 6443 : if (tab->cone) {
2142 0 : isl_int_swap(eq[0], cst);
2143 0 : isl_int_clear(cst);
2144 : }
2145 6443 : if (r < 0)
2146 0 : return -1;
2147 :
2148 6443 : var = &tab->con[r];
2149 6443 : row = var->index;
2150 6443 : if (row_is_manifestly_zero(tab, row)) {
2151 430 : if (snap)
2152 430 : return isl_tab_rollback(tab, snap);
2153 0 : return drop_row(tab, row);
2154 : }
2155 :
2156 6013 : if (tab->bmap) {
2157 5994 : tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);
2158 5994 : if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
2159 0 : return -1;
2160 5994 : isl_seq_neg(eq, eq, 1 + tab->n_var);
2161 5994 : tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);
2162 5994 : isl_seq_neg(eq, eq, 1 + tab->n_var);
2163 5994 : if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)
2164 0 : return -1;
2165 5994 : if (!tab->bmap)
2166 0 : return -1;
2167 5994 : if (add_zero_row(tab) < 0)
2168 0 : return -1;
2169 : }
2170 :
2171 6013 : sgn = isl_int_sgn(tab->mat->row[row][1]);
2172 :
2173 6013 : if (sgn > 0) {
2174 1294 : isl_seq_neg(tab->mat->row[row] + 1, tab->mat->row[row] + 1,
2175 1294 : 1 + tab->n_col);
2176 1294 : var->negated = 1;
2177 1294 : sgn = -1;
2178 : }
2179 :
2180 6013 : if (sgn < 0) {
2181 2645 : sgn = sign_of_max(tab, var);
2182 2645 : if (sgn < -1)
2183 0 : return -1;
2184 2645 : if (sgn < 0) {
2185 0 : if (isl_tab_mark_empty(tab) < 0)
2186 0 : return -1;
2187 0 : return 0;
2188 : }
2189 : }
2190 :
2191 6013 : var->is_nonneg = 1;
2192 6013 : if (to_col(tab, var) < 0)
2193 0 : return -1;
2194 6013 : var->is_nonneg = 0;
2195 6013 : if (isl_tab_kill_col(tab, var->index) < 0)
2196 0 : return -1;
2197 :
2198 6013 : 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 463343 : 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 463343 : if (!bmap)
2218 0 : return NULL;
2219 :
2220 463343 : total = isl_basic_map_total_dim(bmap);
2221 463343 : div_pos = 1 + total - bmap->n_div + div;
2222 :
2223 463343 : ineq = isl_vec_alloc(bmap->ctx, 1 + total);
2224 463343 : if (!ineq)
2225 0 : return NULL;
2226 :
2227 463343 : isl_seq_cpy(ineq->el, bmap->div[div] + 1, 1 + total);
2228 463343 : isl_int_neg(ineq->el[div_pos], bmap->div[div][0]);
2229 463343 : 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 463343 : 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 463343 : total = isl_basic_map_total_dim(tab->bmap);
2255 463343 : div_pos = 1 + total - tab->bmap->n_div + div;
2256 :
2257 463343 : ineq = ineq_for_div(tab->bmap, div);
2258 463343 : if (!ineq)
2259 0 : goto error;
2260 :
2261 463343 : if (add_ineq) {
2262 31191 : if (add_ineq(user, ineq->el) < 0)
2263 0 : goto error;
2264 : } else {
2265 432152 : if (isl_tab_add_ineq(tab, ineq->el) < 0)
2266 0 : goto error;
2267 : }
2268 :
2269 463343 : isl_seq_neg(ineq->el, tab->bmap->div[div] + 1, 1 + total);
2270 463343 : isl_int_set(ineq->el[div_pos], tab->bmap->div[div][0]);
2271 463343 : isl_int_add(ineq->el[0], ineq->el[0], ineq->el[div_pos]);
2272 463343 : isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);
2273 :
2274 463343 : if (add_ineq) {
2275 31191 : if (add_ineq(user, ineq->el) < 0)
2276 0 : goto error;
2277 : } else {
2278 432152 : if (isl_tab_add_ineq(tab, ineq->el) < 0)
2279 0 : goto error;
2280 : }
2281 :
2282 463343 : isl_vec_free(ineq);
2283 :
2284 463343 : 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 463343 : static int div_is_nonneg(struct isl_tab *tab, __isl_keep isl_vec *div)
2297 : {
2298 : int i;
2299 :
2300 463343 : if (tab->M)
2301 0 : return 1;
2302 :
2303 463343 : if (isl_int_is_neg(div->el[1]))
2304 126433 : return 0;
2305 :
2306 1196111 : for (i = 0; i < tab->n_var; ++i) {
2307 1135318 : if (isl_int_is_neg(div->el[2 + i]))
2308 996 : return 0;
2309 1134322 : if (isl_int_is_zero(div->el[2 + i]))
2310 784581 : continue;
2311 349741 : if (!tab->var[i].is_nonneg)
2312 275121 : return 0;
2313 : }
2314 :
2315 60793 : 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 463343 : 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 463343 : if (!tab || !div)
2337 0 : return -1;
2338 :
2339 463343 : 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 463343 : isl_assert(tab->mat->ctx, tab->bmap, return -1);
2344 463343 : n_div = isl_basic_map_dim(tab->bmap, isl_dim_div);
2345 463343 : o_div = tab->n_var - n_div;
2346 463343 : 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 463343 : nonneg = div_is_nonneg(tab, div);
2351 :
2352 463343 : if (isl_tab_extend_cons(tab, 3) < 0)
2353 0 : return -1;
2354 463343 : if (isl_tab_extend_vars(tab, 1) < 0)
2355 0 : return -1;
2356 463343 : r = isl_tab_insert_var(tab, pos);
2357 463343 : if (r < 0)
2358 0 : return -1;
2359 :
2360 463343 : if (nonneg)
2361 60793 : tab->var[r].is_nonneg = 1;
2362 :
2363 463343 : tab->bmap = isl_basic_map_insert_div(tab->bmap, pos - o_div, div);
2364 463343 : if (!tab->bmap)
2365 0 : return -1;
2366 463343 : if (isl_tab_push_var(tab, isl_tab_undo_bmap_div, &tab->var[r]) < 0)
2367 0 : return -1;
2368 :
2369 463343 : if (add_div_constraints(tab, pos - o_div, add_ineq, user) < 0)
2370 0 : return -1;
2371 :
2372 463343 : 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 432152 : int isl_tab_add_div(struct isl_tab *tab, __isl_keep isl_vec *div)
2379 : {
2380 432152 : if (!tab)
2381 0 : return -1;
2382 432152 : 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 7180445 : __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 7180445 : if (!bmap)
2397 0 : return NULL;
2398 14360890 : tab = isl_tab_alloc(bmap->ctx,
2399 7180445 : isl_basic_map_total_dim(bmap) + bmap->n_ineq + 1,
2400 : isl_basic_map_total_dim(bmap), 0);
2401 7180445 : if (!tab)
2402 0 : return NULL;
2403 7180445 : tab->preserve = track;
2404 7180445 : tab->rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
2405 7180445 : 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 7310980 : for (i = 0; i < bmap->n_eq; ++i) {
2411 130535 : tab = add_eq(tab, bmap->eq[i]);
2412 130535 : if (!tab)
2413 0 : return tab;
2414 : }
2415 65528179 : for (i = 0; i < bmap->n_ineq; ++i) {
2416 61249601 : if (isl_tab_add_ineq(tab, bmap->ineq[i]) < 0)
2417 0 : goto error;
2418 61249601 : if (tab->empty)
2419 2901867 : goto done;
2420 : }
2421 : done:
2422 7180445 : if (track && isl_tab_track_bmap(tab, isl_basic_map_copy(bmap)) < 0)
2423 0 : goto error;
2424 7180445 : return tab;
2425 : error:
2426 0 : isl_tab_free(tab);
2427 0 : return NULL;
2428 : }
2429 :
2430 3421899 : __isl_give struct isl_tab *isl_tab_from_basic_set(
2431 : __isl_keep isl_basic_set *bset, int track)
2432 : {
2433 3421899 : return isl_tab_from_basic_map(bset, track);
2434 : }
2435 :
2436 : /* Construct a tableau corresponding to the recession cone of "bset".
2437 : */
2438 30253 : 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 30253 : unsigned offset = 0;
2445 :
2446 30253 : if (!bset)
2447 0 : return NULL;
2448 30253 : if (parametric)
2449 0 : offset = isl_basic_set_dim(bset, isl_dim_param);
2450 30253 : tab = isl_tab_alloc(bset->ctx, bset->n_eq + bset->n_ineq,
2451 30253 : isl_basic_set_total_dim(bset) - offset, 0);
2452 30253 : if (!tab)
2453 0 : return NULL;
2454 30253 : tab->rational = ISL_F_ISSET(bset, ISL_BASIC_SET_RATIONAL);
2455 30253 : tab->cone = 1;
2456 :
2457 30253 : isl_int_init(cst);
2458 30253 : isl_int_set_si(cst, 0);
2459 30253 : for (i = 0; i < bset->n_eq; ++i) {
2460 0 : isl_int_swap(bset->eq[i][offset], cst);
2461 0 : if (offset > 0) {
2462 0 : if (isl_tab_add_eq(tab, bset->eq[i] + offset) < 0)
2463 0 : goto error;
2464 : } else
2465 0 : tab = add_eq(tab, bset->eq[i]);
2466 0 : isl_int_swap(bset->eq[i][offset], cst);
2467 0 : if (!tab)
2468 0 : goto done;
2469 : }
2470 296238 : for (i = 0; i < bset->n_ineq; ++i) {
2471 : int r;
2472 265985 : isl_int_swap(bset->ineq[i][offset], cst);
2473 265985 : r = isl_tab_add_row(tab, bset->ineq[i] + offset);
2474 265985 : isl_int_swap(bset->ineq[i][offset], cst);
2475 265985 : if (r < 0)
2476 0 : goto error;
2477 265985 : tab->con[r].is_nonneg = 1;
2478 265985 : if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)
2479 0 : goto error;
2480 : }
2481 : done:
2482 30253 : isl_int_clear(cst);
2483 30253 : 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 0 : isl_bool isl_tab_cone_is_bounded(struct isl_tab *tab)
2494 : {
2495 : int i;
2496 :
2497 0 : if (!tab)
2498 0 : return isl_bool_error;
2499 0 : if (tab->empty)
2500 0 : return isl_bool_true;
2501 0 : if (tab->n_dead == tab->n_col)
2502 0 : return isl_bool_true;
2503 :
2504 : for (;;) {
2505 0 : for (i = tab->n_redundant; i < tab->n_row; ++i) {
2506 : struct isl_tab_var *var;
2507 : int sgn;
2508 0 : var = isl_tab_var_from_row(tab, i);
2509 0 : if (!var->is_nonneg)
2510 0 : continue;
2511 0 : sgn = sign_of_max(tab, var);
2512 0 : if (sgn < -1)
2513 0 : return isl_bool_error;
2514 0 : if (sgn != 0)
2515 0 : return isl_bool_false;
2516 0 : if (close_row(tab, var, 0) < 0)
2517 0 : return isl_bool_error;
2518 0 : break;
2519 : }
2520 0 : if (tab->n_dead == tab->n_col)
2521 0 : return isl_bool_true;
2522 0 : if (i == tab->n_row)
2523 0 : return isl_bool_false;
2524 0 : }
2525 : }
2526 :
2527 3729550 : int isl_tab_sample_is_integer(struct isl_tab *tab)
2528 : {
2529 : int i;
2530 :
2531 3729550 : if (!tab)
2532 0 : return -1;
2533 :
2534 23776378 : for (i = 0; i < tab->n_var; ++i) {
2535 : int row;
2536 20277583 : if (!tab->var[i].is_row)
2537 2605207 : continue;
2538 17672376 : row = tab->var[i].index;
2539 17672376 : if (!isl_int_is_divisible_by(tab->mat->row[row][1],
2540 : tab->mat->row[row][0]))
2541 230755 : return 0;
2542 : }
2543 3498795 : return 1;
2544 : }
2545 :
2546 3379675 : static struct isl_vec *extract_integer_sample(struct isl_tab *tab)
2547 : {
2548 : int i;
2549 : struct isl_vec *vec;
2550 :
2551 3379675 : vec = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
2552 3379675 : if (!vec)
2553 0 : return NULL;
2554 :
2555 3379675 : isl_int_set_si(vec->block.data[0], 1);
2556 22651430 : for (i = 0; i < tab->n_var; ++i) {
2557 19271755 : if (!tab->var[i].is_row)
2558 2389100 : isl_int_set_si(vec->block.data[1 + i], 0);
2559 : else {
2560 16882655 : int row = tab->var[i].index;
2561 16882655 : isl_int_divexact(vec->block.data[1 + i],
2562 : tab->mat->row[row][1], tab->mat->row[row][0]);
2563 : }
2564 : }
2565 :
2566 3379675 : return vec;
2567 : }
2568 :
2569 645790 : 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 645790 : if (!tab)
2576 0 : return NULL;
2577 :
2578 645790 : vec = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
2579 645790 : if (!vec)
2580 0 : return NULL;
2581 :
2582 645790 : isl_int_init(m);
2583 :
2584 645790 : isl_int_set_si(vec->block.data[0], 1);
2585 4311583 : for (i = 0; i < tab->n_var; ++i) {
2586 : int row;
2587 3665793 : if (!tab->var[i].is_row) {
2588 763789 : isl_int_set_si(vec->block.data[1 + i], 0);
2589 763789 : continue;
2590 : }
2591 2902004 : row = tab->var[i].index;
2592 2902004 : isl_int_gcd(m, vec->block.data[0], tab->mat->row[row][0]);
2593 2902004 : isl_int_divexact(m, tab->mat->row[row][0], m);
2594 2902004 : isl_seq_scale(vec->block.data, vec->block.data, m, 1 + i);
2595 2902004 : isl_int_divexact(m, vec->block.data[0], tab->mat->row[row][0]);
2596 2902004 : isl_int_mul(vec->block.data[1 + i], m, tab->mat->row[row][1]);
2597 : }
2598 645790 : vec = isl_vec_normalize(vec);
2599 :
2600 645790 : isl_int_clear(m);
2601 645790 : 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 156112 : static void get_rounded_sample_value(struct isl_tab *tab,
2608 : struct isl_tab_var *var, int sgn, isl_int *v)
2609 : {
2610 156112 : if (!var->is_row)
2611 1061 : isl_int_set_si(*v, 0);
2612 155051 : else if (sgn > 0)
2613 149765 : isl_int_cdiv_q(*v, tab->mat->row[var->index][1],
2614 : tab->mat->row[var->index][0]);
2615 : else
2616 5286 : isl_int_fdiv_q(*v, tab->mat->row[var->index][1],
2617 : tab->mat->row[var->index][0]);
2618 156112 : }
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 3828730 : 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 3828730 : if (!bmap)
2635 0 : return NULL;
2636 3828730 : if (!tab)
2637 0 : return bmap;
2638 :
2639 3828730 : n_eq = tab->n_eq;
2640 3828730 : if (tab->empty)
2641 5217 : bmap = isl_basic_map_set_to_empty(bmap);
2642 : else
2643 41727627 : for (i = bmap->n_ineq - 1; i >= 0; --i) {
2644 37904114 : if (isl_tab_is_equality(tab, n_eq + i))
2645 23978381 : isl_basic_map_inequality_to_equality(bmap, i);
2646 13925733 : else if (isl_tab_is_redundant(tab, n_eq + i))
2647 5339253 : isl_basic_map_drop_inequality(bmap, i);
2648 : }
2649 3828730 : if (bmap->n_eq != n_eq)
2650 2931904 : bmap = isl_basic_map_gauss(bmap, NULL);
2651 3828730 : if (!tab->rational &&
2652 3828730 : bmap && !bmap->sample && isl_tab_sample_is_integer(tab))
2653 3379675 : bmap->sample = extract_integer_sample(tab);
2654 3828730 : return bmap;
2655 : }
2656 :
2657 89647 : struct isl_basic_set *isl_basic_set_update_from_tab(struct isl_basic_set *bset,
2658 : struct isl_tab *tab)
2659 : {
2660 89647 : 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 32071 : static isl_stat drop_last_con_in_row(struct isl_tab *tab, int r)
2668 : {
2669 32071 : 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 32071 : 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 32071 : if (drop_row(tab, tab->con[r].index) < 0)
2677 0 : return isl_stat_error;
2678 :
2679 32071 : 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 32071 : 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 32071 : unsigned off = 2 + tab->M;
2698 :
2699 32071 : if (var->is_zero)
2700 0 : return isl_stat_ok;
2701 32071 : 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 32071 : if (isl_tab_extend_cons(tab, 1) < 0)
2707 0 : return isl_stat_error;
2708 :
2709 32071 : r = tab->n_con;
2710 32071 : tab->con[r].index = tab->n_row;
2711 32071 : tab->con[r].is_row = 1;
2712 32071 : tab->con[r].is_nonneg = 0;
2713 32071 : tab->con[r].is_zero = 0;
2714 32071 : tab->con[r].is_redundant = 0;
2715 32071 : tab->con[r].frozen = 0;
2716 32071 : tab->con[r].negated = 0;
2717 32071 : tab->row_var[tab->n_row] = ~r;
2718 32071 : row = tab->mat->row[tab->n_row];
2719 :
2720 32071 : if (var->is_row) {
2721 30897 : isl_int_set(row[0], tab->mat->row[var->index][0]);
2722 61794 : isl_seq_neg(row + 1,
2723 61794 : tab->mat->row[var->index] + 1, 1 + tab->n_col);
2724 : } else {
2725 1174 : isl_int_set_si(row[0], 1);
2726 1174 : isl_seq_clr(row + 1, 1 + tab->n_col);
2727 1174 : isl_int_set_si(row[off + var->index], -1);
2728 : }
2729 :
2730 32071 : tab->n_row++;
2731 32071 : tab->n_con++;
2732 :
2733 32071 : sgn = sign_of_max(tab, &tab->con[r]);
2734 32071 : if (sgn < -1)
2735 0 : return isl_stat_error;
2736 32071 : if (sgn < 0) {
2737 69 : if (drop_last_con_in_row(tab, r) < 0)
2738 0 : return isl_stat_error;
2739 69 : if (isl_tab_mark_empty(tab) < 0)
2740 0 : return isl_stat_error;
2741 69 : return isl_stat_ok;
2742 : }
2743 32002 : tab->con[r].is_nonneg = 1;
2744 : /* sgn == 0 */
2745 32002 : if (close_row(tab, &tab->con[r], 1) < 0)
2746 0 : return isl_stat_error;
2747 32002 : if (drop_last_con_in_row(tab, r) < 0)
2748 0 : return isl_stat_error;
2749 :
2750 32002 : return isl_stat_ok;
2751 : }
2752 :
2753 : /* Check that "con" is a valid constraint position for "tab".
2754 : */
2755 15079105 : static isl_stat isl_tab_check_con(struct isl_tab *tab, int con)
2756 : {
2757 15079105 : if (!tab)
2758 0 : return isl_stat_error;
2759 15079105 : 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 15079105 : 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 65 : int isl_tab_relax(struct isl_tab *tab, int con)
2783 : {
2784 : struct isl_tab_var *var;
2785 :
2786 65 : if (!tab)
2787 0 : return -1;
2788 :
2789 65 : var = &tab->con[con];
2790 :
2791 65 : 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 65 : 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 65 : if (!var->is_row && !max_is_manifestly_unbounded(tab, var))
2799 4 : if (to_row(tab, var, 1) < 0)
2800 0 : return -1;
2801 65 : if (!var->is_row && !min_is_manifestly_unbounded(tab, var))
2802 15 : if (to_row(tab, var, -1) < 0)
2803 0 : return -1;
2804 :
2805 65 : if (var->is_row) {
2806 19 : isl_int_add(tab->mat->row[var->index][1],
2807 : tab->mat->row[var->index][1], tab->mat->row[var->index][0]);
2808 19 : if (restore_row(tab, var) < 0)
2809 0 : return -1;
2810 : } else {
2811 : int i;
2812 46 : unsigned off = 2 + tab->M;
2813 :
2814 213 : for (i = 0; i < tab->n_row; ++i) {
2815 167 : if (isl_int_is_zero(tab->mat->row[i][off + var->index]))
2816 71 : continue;
2817 96 : 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 65 : if (isl_tab_push_var(tab, isl_tab_undo_relax, var) < 0)
2824 0 : return -1;
2825 :
2826 65 : 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 848 : int isl_tab_unrestrict(struct isl_tab *tab, int con)
2893 : {
2894 : struct isl_tab_var *var;
2895 :
2896 848 : if (!tab)
2897 0 : return -1;
2898 :
2899 848 : var = &tab->con[con];
2900 848 : if (!var->is_nonneg)
2901 0 : return 0;
2902 :
2903 848 : var->is_nonneg = 0;
2904 848 : if (isl_tab_push_var(tab, isl_tab_undo_unrestrict, var) < 0)
2905 0 : return -1;
2906 :
2907 848 : return 0;
2908 : }
2909 :
2910 1315 : int isl_tab_select_facet(struct isl_tab *tab, int con)
2911 : {
2912 1315 : if (!tab)
2913 0 : return -1;
2914 :
2915 1315 : return cut_to_hyperplane(tab, &tab->con[con]);
2916 : }
2917 :
2918 83817619 : static int may_be_equality(struct isl_tab *tab, int row)
2919 : {
2920 167635238 : return tab->rational ? isl_int_is_zero(tab->mat->row[row][1])
2921 83817619 : : 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 24997732 : static struct isl_tab_var *select_marked(struct isl_tab *tab)
2937 : {
2938 : int i;
2939 : struct isl_tab_var *var;
2940 :
2941 152235921 : for (i = tab->n_con - 1; i >= 0; --i) {
2942 148957216 : var = &tab->con[i];
2943 148957216 : if (var->index < 0)
2944 90858419 : continue;
2945 58098797 : if (var->is_row && var->index < tab->n_redundant)
2946 1110222 : continue;
2947 56988575 : if (!var->is_row && var->index < tab->n_dead)
2948 132411 : continue;
2949 56856164 : if (var->marked)
2950 21719027 : return var;
2951 : }
2952 :
2953 3278705 : 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 3803028 : int isl_tab_detect_implicit_equalities(struct isl_tab *tab)
2974 : {
2975 : int i;
2976 : unsigned n_marked;
2977 :
2978 3803028 : if (!tab)
2979 0 : return -1;
2980 3803028 : if (tab->empty)
2981 5644 : return 0;
2982 3797384 : if (tab->n_dead == tab->n_col)
2983 11484 : return 0;
2984 :
2985 3785900 : n_marked = 0;
2986 40225512 : for (i = tab->n_redundant; i < tab->n_row; ++i) {
2987 36439612 : struct isl_tab_var *var = isl_tab_var_from_row(tab, i);
2988 71808485 : var->marked = !var->frozen && var->is_nonneg &&
2989 35368873 : may_be_equality(tab, i);
2990 36439612 : if (var->marked)
2991 34741255 : n_marked++;
2992 : }
2993 25124225 : for (i = tab->n_dead; i < tab->n_col; ++i) {
2994 21338325 : struct isl_tab_var *var = var_from_col(tab, i);
2995 21338325 : var->marked = !var->frozen && var->is_nonneg;
2996 21338325 : if (var->marked)
2997 1500630 : n_marked++;
2998 : }
2999 27689363 : while (n_marked) {
3000 : struct isl_tab_var *var;
3001 : int sgn;
3002 23344999 : var = select_marked(tab);
3003 23344999 : if (!var)
3004 3196680 : break;
3005 20148319 : var->marked = 0;
3006 20148319 : n_marked--;
3007 20148319 : sgn = sign_of_max(tab, var);
3008 20148319 : if (sgn < 0)
3009 0 : return -1;
3010 20148319 : if (sgn == 0) {
3011 11580754 : if (close_row(tab, var, 0) < 0)
3012 0 : return -1;
3013 8567565 : } else if (!tab->rational && !at_least_one(tab, var)) {
3014 30756 : if (cut_to_hyperplane(tab, var) < 0)
3015 0 : return -1;
3016 30756 : return isl_tab_detect_implicit_equalities(tab);
3017 : }
3018 146372935 : for (i = tab->n_redundant; i < tab->n_row; ++i) {
3019 126255372 : var = isl_tab_var_from_row(tab, i);
3020 126255372 : if (!var->marked)
3021 77806626 : continue;
3022 48448746 : if (may_be_equality(tab, i))
3023 48401577 : continue;
3024 47169 : var->marked = 0;
3025 47169 : n_marked--;
3026 : }
3027 : }
3028 :
3029 3755144 : 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 0 : static int update_con_after_move(struct isl_tab *tab, int i, int old)
3036 : {
3037 : int *p;
3038 : int index;
3039 :
3040 0 : index = tab->con[i].index;
3041 0 : if (index == -1)
3042 0 : return 0;
3043 0 : p = tab->con[i].is_row ? tab->row_var : tab->col_var;
3044 0 : if (p[index] != ~old)
3045 0 : isl_die(tab->mat->ctx, isl_error_internal,
3046 : "broken internal state", return -1);
3047 0 : p[index] = ~i;
3048 :
3049 0 : 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 0 : isl_stat isl_tab_swap_constraints(struct isl_tab *tab, int con1, int con2)
3058 : {
3059 : struct isl_tab_var var;
3060 :
3061 0 : if (isl_tab_check_con(tab, con1) < 0 ||
3062 0 : isl_tab_check_con(tab, con2) < 0)
3063 0 : return isl_stat_error;
3064 :
3065 0 : var = tab->con[con1];
3066 0 : tab->con[con1] = tab->con[con2];
3067 0 : if (update_con_after_move(tab, con1, con2) < 0)
3068 0 : return isl_stat_error;
3069 0 : tab->con[con2] = var;
3070 0 : if (update_con_after_move(tab, con2, con1) < 0)
3071 0 : return isl_stat_error;
3072 :
3073 0 : 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 0 : 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 0 : if (n <= 1)
3085 0 : return 0;
3086 :
3087 0 : last = first + n - 1;
3088 0 : var = tab->con[last];
3089 0 : for (i = last; i > first; --i) {
3090 0 : tab->con[i] = tab->con[i - 1];
3091 0 : if (update_con_after_move(tab, i, i - 1) < 0)
3092 0 : return -1;
3093 : }
3094 0 : tab->con[first] = var;
3095 0 : if (update_con_after_move(tab, first, last) < 0)
3096 0 : return -1;
3097 :
3098 0 : 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 0 : static isl_stat con_drop_entries(struct isl_tab *tab,
3107 : unsigned first, unsigned n)
3108 : {
3109 : int i;
3110 :
3111 0 : 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 0 : tab->n_con -= n;
3116 :
3117 0 : for (i = first; i < tab->n_con; ++i) {
3118 0 : tab->con[i] = tab->con[i + n];
3119 0 : if (update_con_after_move(tab, i, i + n) < 0)
3120 0 : return isl_stat_error;
3121 : }
3122 :
3123 0 : 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 0 : static isl_stat swap_eq(unsigned a, unsigned b, void *user)
3131 : {
3132 0 : struct isl_tab *tab = user;
3133 :
3134 0 : 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 0 : static isl_stat drop_eq(unsigned n, void *user)
3146 : {
3147 0 : struct isl_tab *tab = user;
3148 :
3149 0 : if (tab->n_eq == n)
3150 0 : return isl_tab_mark_empty(tab);
3151 :
3152 0 : tab->n_eq -= n;
3153 0 : 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 0 : 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 0 : single = isl_basic_map_has_single_reference(bmap);
3165 0 : if (single < 0)
3166 0 : return isl_basic_map_free(bmap);
3167 0 : if (single)
3168 0 : return bmap;
3169 0 : 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 1580 : __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 1580 : if (!tab || !bmap)
3212 0 : return isl_basic_map_free(bmap);
3213 1580 : if (tab->empty)
3214 0 : return bmap;
3215 :
3216 1580 : n_eq = tab->n_eq;
3217 12001 : for (i = bmap->n_ineq - 1; i >= 0; --i) {
3218 10421 : if (!isl_tab_is_equality(tab, bmap->n_eq + i))
3219 10421 : continue;
3220 0 : isl_basic_map_inequality_to_equality(bmap, i);
3221 0 : if (rotate_constraints(tab, 0, tab->n_eq + i + 1) < 0)
3222 0 : return isl_basic_map_free(bmap);
3223 0 : if (rotate_constraints(tab, tab->n_eq + i + 1,
3224 0 : bmap->n_ineq - i) < 0)
3225 0 : return isl_basic_map_free(bmap);
3226 0 : tab->n_eq++;
3227 : }
3228 :
3229 1580 : if (n_eq != tab->n_eq)
3230 0 : bmap = gauss_if_shared(bmap, tab);
3231 :
3232 1580 : return bmap;
3233 : }
3234 :
3235 1689588 : static int con_is_redundant(struct isl_tab *tab, struct isl_tab_var *var)
3236 : {
3237 1689588 : if (!tab)
3238 0 : return -1;
3239 1689588 : if (tab->rational) {
3240 535 : int sgn = sign_of_min(tab, var);
3241 535 : if (sgn < -1)
3242 0 : return -1;
3243 535 : return sgn >= 0;
3244 : } else {
3245 1689053 : int irred = isl_tab_min_at_most_neg_one(tab, var);
3246 1689053 : if (irred < 0)
3247 0 : return -1;
3248 1689053 : 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 544802 : int isl_tab_detect_redundant(struct isl_tab *tab)
3266 : {
3267 : int i;
3268 : unsigned n_marked;
3269 :
3270 544802 : if (!tab)
3271 0 : return -1;
3272 544802 : if (tab->empty)
3273 5217 : return 0;
3274 539585 : if (tab->n_redundant == tab->n_row)
3275 843 : return 0;
3276 :
3277 538742 : n_marked = 0;
3278 4272804 : for (i = tab->n_redundant; i < tab->n_row; ++i) {
3279 3734062 : struct isl_tab_var *var = isl_tab_var_from_row(tab, i);
3280 3734062 : var->marked = !var->frozen && var->is_nonneg;
3281 3734062 : if (var->marked)
3282 1040872 : n_marked++;
3283 : }
3284 3685200 : for (i = tab->n_dead; i < tab->n_col; ++i) {
3285 3146458 : struct isl_tab_var *var = var_from_col(tab, i);
3286 4702835 : var->marked = !var->frozen && var->is_nonneg &&
3287 1556377 : !min_is_manifestly_unbounded(tab, var);
3288 3146458 : if (var->marked)
3289 867987 : n_marked++;
3290 : }
3291 2648192 : while (n_marked) {
3292 : struct isl_tab_var *var;
3293 : int red;
3294 1652733 : var = select_marked(tab);
3295 1652733 : if (!var)
3296 82025 : break;
3297 1570708 : var->marked = 0;
3298 1570708 : n_marked--;
3299 1570708 : red = con_is_redundant(tab, var);
3300 1570708 : if (red < 0)
3301 0 : return -1;
3302 1570708 : if (red && !var->is_redundant)
3303 33691 : if (isl_tab_mark_redundant(tab, var->index) < 0)
3304 0 : return -1;
3305 10973328 : for (i = tab->n_dead; i < tab->n_col; ++i) {
3306 9402620 : var = var_from_col(tab, i);
3307 9402620 : if (!var->marked)
3308 7293596 : continue;
3309 2109024 : if (!min_is_manifestly_unbounded(tab, var))
3310 1902236 : continue;
3311 206788 : var->marked = 0;
3312 206788 : n_marked--;
3313 : }
3314 : }
3315 :
3316 538742 : return 0;
3317 : }
3318 :
3319 38203287 : int isl_tab_is_equality(struct isl_tab *tab, int con)
3320 : {
3321 : int row;
3322 : unsigned off;
3323 :
3324 38203287 : if (!tab)
3325 0 : return -1;
3326 38203287 : if (tab->con[con].is_zero)
3327 24045945 : return 1;
3328 14157342 : if (tab->con[con].is_redundant)
3329 5396171 : return 0;
3330 8761171 : if (!tab->con[con].is_row)
3331 6536110 : return tab->con[con].index < tab->n_dead;
3332 :
3333 2225061 : row = tab->con[con].index;
3334 :
3335 2225061 : off = 2 + tab->M;
3336 8303611 : return isl_int_is_zero(tab->mat->row[row][1]) &&
3337 5482057 : !row_is_big(tab, row) &&
3338 1628498 : isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
3339 1628498 : 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 154411 : 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 154411 : enum isl_lp_result res = isl_lp_ok;
3358 : struct isl_tab_var *var;
3359 : struct isl_tab_undo *snap;
3360 :
3361 154411 : if (!tab)
3362 0 : return isl_lp_error;
3363 :
3364 154411 : if (tab->empty)
3365 0 : return isl_lp_empty;
3366 :
3367 154411 : snap = isl_tab_snap(tab);
3368 154411 : r = isl_tab_add_row(tab, f);
3369 154411 : if (r < 0)
3370 0 : return isl_lp_error;
3371 154411 : var = &tab->con[r];
3372 : for (;;) {
3373 : int row, col;
3374 264060 : find_pivot(tab, var, var, -1, &row, &col);
3375 264060 : if (row == var->index) {
3376 101 : res = isl_lp_unbounded;
3377 154512 : break;
3378 : }
3379 263959 : if (row == -1)
3380 154310 : break;
3381 109649 : if (isl_tab_pivot(tab, row, col) < 0)
3382 0 : return isl_lp_error;
3383 109649 : }
3384 154411 : isl_int_mul(tab->mat->row[var->index][0],
3385 : tab->mat->row[var->index][0], denom);
3386 154411 : if (ISL_FL_ISSET(flags, ISL_TAB_SAVE_DUAL)) {
3387 : int i;
3388 :
3389 3094 : isl_vec_free(tab->dual);
3390 3094 : tab->dual = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_con);
3391 3094 : if (!tab->dual)
3392 0 : return isl_lp_error;
3393 3094 : isl_int_set(tab->dual->el[0], tab->mat->row[var->index][0]);
3394 93764 : for (i = 0; i < tab->n_con; ++i) {
3395 : int pos;
3396 90670 : if (tab->con[i].is_row) {
3397 59622 : isl_int_set_si(tab->dual->el[1 + i], 0);
3398 59622 : continue;
3399 : }
3400 31048 : pos = 2 + tab->M + tab->con[i].index;
3401 31048 : if (tab->con[i].negated)
3402 3496 : isl_int_neg(tab->dual->el[1 + i],
3403 : tab->mat->row[var->index][pos]);
3404 : else
3405 27552 : isl_int_set(tab->dual->el[1 + i],
3406 : tab->mat->row[var->index][pos]);
3407 : }
3408 : }
3409 154411 : if (opt && res == isl_lp_ok) {
3410 154310 : if (opt_denom) {
3411 6695 : isl_int_set(*opt, tab->mat->row[var->index][1]);
3412 6695 : isl_int_set(*opt_denom, tab->mat->row[var->index][0]);
3413 : } else
3414 147615 : get_rounded_sample_value(tab, var, 1, opt);
3415 : }
3416 154411 : if (isl_tab_rollback(tab, snap) < 0)
3417 0 : return isl_lp_error;
3418 154411 : 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 15079105 : int isl_tab_is_redundant(struct isl_tab *tab, int con)
3430 : {
3431 15079105 : if (isl_tab_check_con(tab, con) < 0)
3432 0 : return -1;
3433 15079105 : if (tab->con[con].is_zero)
3434 0 : return 0;
3435 15079105 : if (tab->con[con].is_redundant)
3436 6150378 : return 1;
3437 8928727 : 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 6369 : static isl_bool is_constant(struct isl_tab *tab, struct isl_tab_var *var,
3449 : isl_int *value)
3450 : {
3451 6369 : unsigned off = 2 + tab->M;
3452 6369 : isl_mat *mat = tab->mat;
3453 : int n;
3454 : int row;
3455 : int pos;
3456 :
3457 6369 : if (!var->is_row)
3458 1061 : return isl_bool_false;
3459 5308 : row = var->index;
3460 5308 : if (row_is_big(tab, row))
3461 0 : return isl_bool_false;
3462 5308 : n = tab->n_col - tab->n_dead;
3463 5308 : pos = isl_seq_first_non_zero(mat->row[row] + off + tab->n_dead, n);
3464 5308 : if (pos != -1)
3465 5286 : return isl_bool_false;
3466 22 : if (value)
3467 0 : isl_int_divexact(*value, mat->row[row][1], mat->row[row][0]);
3468 22 : 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 11677 : static int reached(struct isl_tab *tab, struct isl_tab_var *var, int sgn,
3482 : isl_int target, isl_int *tmp)
3483 : {
3484 11677 : if (row_is_big(tab, var->index))
3485 0 : return 1;
3486 11677 : isl_int_mul(*tmp, tab->mat->row[var->index][0], target);
3487 11677 : if (sgn > 0)
3488 3820 : return isl_int_ge(tab->mat->row[var->index][1], *tmp);
3489 : else
3490 7857 : 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 8497 : 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 8497 : if (sgn < 0 && min_is_manifestly_unbounded(tab, var))
3515 771 : return isl_bool_true;
3516 7726 : if (sgn > 0 && max_is_manifestly_unbounded(tab, var))
3517 0 : return isl_bool_true;
3518 7726 : if (to_row(tab, var, sgn) < 0)
3519 0 : return isl_bool_error;
3520 19403 : while (!reached(tab, var, sgn, target, tmp)) {
3521 10029 : find_pivot(tab, var, var, sgn, &row, &col);
3522 10029 : if (row == -1)
3523 2150 : return isl_bool_false;
3524 7879 : if (row == var->index)
3525 3928 : return isl_bool_true;
3526 3951 : if (isl_tab_pivot(tab, row, col) < 0)
3527 0 : return isl_bool_error;
3528 : }
3529 :
3530 1648 : 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 6347 : 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 6347 : get_rounded_sample_value(tab, var, -1, target);
3561 6347 : isl_int_sub_ui(*target, *target, 1);
3562 6347 : reached = var_reaches(tab, var, -1, *target, tmp);
3563 6347 : if (reached < 0 || reached)
3564 4197 : return isl_bool_not(reached);
3565 2150 : get_rounded_sample_value(tab, var, 1, target);
3566 2150 : isl_int_add_ui(*target, *target, 1);
3567 2150 : reached = var_reaches(tab, var, 1, *target, tmp);
3568 2150 : if (reached < 0 || reached)
3569 2150 : return isl_bool_not(reached);
3570 0 : get_rounded_sample_value(tab, var, -1, tmp);
3571 0 : isl_int_sub_ui(*target, *target, 1);
3572 0 : 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 0 : if (isl_tab_extend_cons(tab, 1) < 0)
3579 0 : return isl_bool_error;
3580 0 : eq = isl_vec_alloc(isl_tab_get_ctx(tab), 1 + tab->n_var);
3581 0 : if (!eq)
3582 0 : return isl_bool_error;
3583 0 : pos = var - tab->var;
3584 0 : isl_seq_clr(eq->el + 1, tab->n_var);
3585 0 : isl_int_set_si(eq->el[1 + pos], -1);
3586 0 : isl_int_set(eq->el[0], *target);
3587 0 : r = isl_tab_add_eq(tab, eq->el);
3588 0 : isl_vec_free(eq);
3589 :
3590 0 : 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 6369 : 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 6369 : if (var->is_row && row_is_big(tab, var->index))
3613 0 : return isl_bool_false;
3614 6369 : is_cst = is_constant(tab, var, value);
3615 6369 : if (is_cst < 0 || is_cst)
3616 22 : return is_cst;
3617 :
3618 6347 : if (!value)
3619 5603 : isl_int_init(target);
3620 6347 : isl_int_init(tmp);
3621 :
3622 6347 : is_cst = detect_constant_with_tmp(tab, var,
3623 : value ? value : &target, &tmp);
3624 :
3625 6347 : isl_int_clear(tmp);
3626 6347 : if (!value)
3627 5603 : isl_int_clear(target);
3628 :
3629 6347 : 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 744 : isl_bool isl_tab_is_constant(struct isl_tab *tab, int var, isl_int *value)
3640 : {
3641 744 : if (!tab)
3642 0 : return isl_bool_error;
3643 744 : 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 744 : if (tab->rational)
3647 0 : return isl_bool_false;
3648 :
3649 744 : 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 848 : isl_stat isl_tab_detect_constants(struct isl_tab *tab)
3659 : {
3660 : int i;
3661 :
3662 848 : if (!tab)
3663 0 : return isl_stat_error;
3664 848 : if (tab->rational)
3665 0 : return isl_stat_ok;
3666 :
3667 6473 : for (i = 0; i < tab->n_var; ++i) {
3668 5625 : if (get_constant(tab, &tab->var[i], NULL) < 0)
3669 0 : return isl_stat_error;
3670 : }
3671 :
3672 848 : 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 12570388 : struct isl_tab_undo *isl_tab_snap(struct isl_tab *tab)
3679 : {
3680 12570388 : if (!tab)
3681 0 : return NULL;
3682 12570388 : tab->need_undo = 1;
3683 12570388 : 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 1146 : isl_bool isl_tab_need_undo(struct isl_tab *tab)
3690 : {
3691 1146 : if (!tab)
3692 0 : return isl_bool_error;
3693 :
3694 1146 : 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 1146 : void isl_tab_clear_undo(struct isl_tab *tab)
3703 : {
3704 1146 : if (!tab)
3705 0 : return;
3706 :
3707 1146 : free_undo(tab);
3708 1146 : 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 38 : static isl_stat unrelax(struct isl_tab *tab, struct isl_tab_var *var)
3716 : {
3717 38 : unsigned off = 2 + tab->M;
3718 :
3719 38 : if (!var->is_row && !max_is_manifestly_unbounded(tab, var))
3720 8 : if (to_row(tab, var, 1) < 0)
3721 0 : return isl_stat_error;
3722 :
3723 38 : if (var->is_row) {
3724 8 : isl_int_sub(tab->mat->row[var->index][1],
3725 : tab->mat->row[var->index][1], tab->mat->row[var->index][0]);
3726 8 : if (var->is_nonneg) {
3727 8 : int sgn = restore_row(tab, var);
3728 8 : isl_assert(tab->mat->ctx, sgn >= 0,
3729 : return isl_stat_error);
3730 : }
3731 : } else {
3732 : int i;
3733 :
3734 162 : for (i = 0; i < tab->n_row; ++i) {
3735 132 : if (isl_int_is_zero(tab->mat->row[i][off + var->index]))
3736 47 : continue;
3737 85 : 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 38 : 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 847 : static isl_stat ununrestrict(struct isl_tab *tab, struct isl_tab_var *var)
3752 : {
3753 847 : var->is_nonneg = 1;
3754 :
3755 847 : if (var->is_row && restore_row(tab, var) < -1)
3756 0 : return isl_stat_error;
3757 :
3758 847 : 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 8511602 : static isl_stat restore_last_redundant(struct isl_tab *tab)
3769 : {
3770 : struct isl_tab_var *var;
3771 :
3772 8511602 : 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 8511602 : var = isl_tab_var_from_row(tab, tab->n_redundant - 1);
3777 8511602 : var->is_redundant = 0;
3778 8511602 : tab->n_redundant--;
3779 8511602 : restore_row(tab, var);
3780 :
3781 8511602 : 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 95361284 : static isl_stat perform_undo_var(struct isl_tab *tab, struct isl_tab_undo *undo)
3787 : {
3788 95361284 : struct isl_tab_var *var = var_from_index(tab, undo->u.var_index);
3789 95361284 : switch (undo->type) {
3790 : case isl_tab_undo_nonneg:
3791 19288180 : var->is_nonneg = 0;
3792 19288180 : break;
3793 : case isl_tab_undo_redundant:
3794 7800185 : 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 7800185 : return restore_last_redundant(tab);
3799 : case isl_tab_undo_freeze:
3800 49206564 : var->frozen = 0;
3801 49206564 : break;
3802 : case isl_tab_undo_zero:
3803 104368 : var->is_zero = 0;
3804 104368 : if (!var->is_row)
3805 77911 : tab->n_dead--;
3806 104368 : break;
3807 : case isl_tab_undo_allocate:
3808 18961102 : if (undo->u.var_index >= 0) {
3809 467087 : isl_assert(tab->mat->ctx, !var->is_row,
3810 : return isl_stat_error);
3811 467087 : return drop_col(tab, var->index);
3812 : }
3813 18494015 : if (!var->is_row) {
3814 3457099 : if (!max_is_manifestly_unbounded(tab, var)) {
3815 954767 : if (to_row(tab, var, 1) < 0)
3816 0 : return isl_stat_error;
3817 2502332 : } else if (!min_is_manifestly_unbounded(tab, var)) {
3818 1551166 : if (to_row(tab, var, -1) < 0)
3819 0 : return isl_stat_error;
3820 : } else
3821 951166 : if (to_row(tab, var, 0) < 0)
3822 0 : return isl_stat_error;
3823 : }
3824 18494015 : return drop_row(tab, var->index);
3825 : case isl_tab_undo_relax:
3826 38 : return unrelax(tab, var);
3827 : case isl_tab_undo_unrestrict:
3828 847 : 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 68599112 : 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 366304 : isl_stat isl_tab_restore_redundant(struct isl_tab *tab)
3846 : {
3847 366304 : if (!tab)
3848 0 : return isl_stat_error;
3849 :
3850 366304 : 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 1444025 : while (tab->n_redundant > 0) {
3857 711417 : if (tab->row_var[tab->n_redundant - 1] >= 0) {
3858 : struct isl_tab_var *var;
3859 :
3860 607218 : var = isl_tab_var_from_row(tab, tab->n_redundant - 1);
3861 607218 : var->is_nonneg = 0;
3862 : }
3863 711417 : restore_last_redundant(tab);
3864 : }
3865 366304 : 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 466343 : static isl_stat drop_bmap_div(struct isl_tab *tab, int pos)
3872 : {
3873 : int off;
3874 :
3875 466343 : off = tab->n_var - isl_basic_map_dim(tab->bmap, isl_dim_div);
3876 466343 : if (isl_basic_map_drop_div(tab->bmap, pos - off) < 0)
3877 : return isl_stat_error;
3878 466343 : if (tab->samples) {
3879 31140 : tab->samples = isl_mat_drop_cols(tab->samples, 1 + pos, 1);
3880 31140 : if (!tab->samples)
3881 0 : return isl_stat_error;
3882 : }
3883 :
3884 466343 : 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 136010 : static void drop_samples_since(struct isl_tab *tab, int n)
3948 : {
3949 : int i;
3950 :
3951 375144 : for (i = tab->n_sample - 1; i >= 0 && tab->n_sample > n; --i) {
3952 239134 : if (tab->sample_index[i] < n)
3953 69906 : continue;
3954 :
3955 169228 : if (i != tab->n_sample - 1) {
3956 93403 : int t = tab->sample_index[tab->n_sample-1];
3957 93403 : tab->sample_index[tab->n_sample-1] = tab->sample_index[i];
3958 93403 : tab->sample_index[i] = t;
3959 93403 : isl_mat_swap_rows(tab->samples, tab->n_sample-1, i);
3960 : }
3961 169228 : tab->n_sample--;
3962 : }
3963 136010 : }
3964 :
3965 : static isl_stat perform_undo(struct isl_tab *tab, struct isl_tab_undo *undo)
3966 : WARN_UNUSED;
3967 119994727 : static isl_stat perform_undo(struct isl_tab *tab, struct isl_tab_undo *undo)
3968 : {
3969 119994727 : switch (undo->type) {
3970 : case isl_tab_undo_rational:
3971 1230 : tab->rational = 0;
3972 1230 : break;
3973 : case isl_tab_undo_empty:
3974 5475161 : tab->empty = 0;
3975 5475161 : 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 95361284 : 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 18070665 : return isl_basic_map_free_inequality(tab->bmap, 1);
3988 : case isl_tab_undo_bmap_div:
3989 466343 : 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 363460 : tab->n_outside--;
3996 363460 : break;
3997 : case isl_tab_undo_saved_samples:
3998 136010 : drop_samples_since(tab, undo->u.n);
3999 136010 : break;
4000 : case isl_tab_undo_callback:
4001 120574 : return undo->u.callback->run(undo->u.callback);
4002 : default:
4003 0 : isl_assert(tab->mat->ctx, 0, return isl_stat_error);
4004 : }
4005 5975861 : 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 7400786 : int isl_tab_rollback(struct isl_tab *tab, struct isl_tab_undo *snap)
4012 : {
4013 : struct isl_tab_undo *undo, *next;
4014 :
4015 7400786 : if (!tab)
4016 0 : return -1;
4017 :
4018 7400786 : tab->in_undo = 1;
4019 127395513 : for (undo = tab->top; undo && undo != &tab->bottom; undo = next) {
4020 126044314 : next = undo->next;
4021 126044314 : if (undo == snap)
4022 6049587 : break;
4023 119994727 : 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 119994727 : free_undo_record(undo);
4030 : }
4031 7400786 : tab->in_undo = 0;
4032 7400786 : tab->top = undo;
4033 7400786 : if (!undo)
4034 0 : return -1;
4035 7400786 : 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 41132 : static enum isl_ineq_type separation_type(struct isl_tab *tab, unsigned row)
4049 : {
4050 : int pos;
4051 41132 : unsigned off = 2 + tab->M;
4052 :
4053 41132 : if (tab->rational)
4054 36 : return isl_ineq_separate;
4055 :
4056 41096 : if (!isl_int_is_one(tab->mat->row[row][0]))
4057 2472 : return isl_ineq_separate;
4058 :
4059 38624 : pos = isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,
4060 38624 : tab->n_col - tab->n_dead);
4061 38624 : if (pos == -1) {
4062 959 : if (isl_int_is_negone(tab->mat->row[row][1]))
4063 958 : return isl_ineq_adj_eq;
4064 : else
4065 1 : return isl_ineq_separate;
4066 : }
4067 :
4068 37665 : if (!isl_int_eq(tab->mat->row[row][1],
4069 : tab->mat->row[row][off + tab->n_dead + pos]))
4070 15685 : return isl_ineq_separate;
4071 :
4072 43960 : pos = isl_seq_first_non_zero(
4073 21980 : tab->mat->row[row] + off + tab->n_dead + pos + 1,
4074 21980 : tab->n_col - tab->n_dead - pos - 1);
4075 :
4076 21980 : 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 183111 : enum isl_ineq_type isl_tab_ineq_type(struct isl_tab *tab, isl_int *ineq)
4088 : {
4089 183111 : enum isl_ineq_type type = isl_ineq_error;
4090 183111 : struct isl_tab_undo *snap = NULL;
4091 : int con;
4092 : int row;
4093 :
4094 183111 : if (!tab)
4095 0 : return isl_ineq_error;
4096 :
4097 183111 : if (isl_tab_extend_cons(tab, 1) < 0)
4098 0 : return isl_ineq_error;
4099 :
4100 183111 : snap = isl_tab_snap(tab);
4101 :
4102 183111 : con = isl_tab_add_row(tab, ineq);
4103 183111 : if (con < 0)
4104 0 : goto error;
4105 :
4106 183111 : row = tab->con[con].index;
4107 183111 : if (isl_tab_row_is_redundant(tab, row))
4108 0 : type = isl_ineq_redundant;
4109 249712 : else if (isl_int_is_neg(tab->mat->row[row][1]) &&
4110 132416 : (tab->rational ||
4111 65815 : isl_int_abs_ge(tab->mat->row[row][1],
4112 64231 : tab->mat->row[row][0]))) {
4113 64231 : int nonneg = at_least_zero(tab, &tab->con[con]);
4114 64231 : if (nonneg < 0)
4115 0 : goto error;
4116 64231 : if (nonneg)
4117 23099 : type = isl_ineq_cut;
4118 : else
4119 41132 : type = separation_type(tab, row);
4120 : } else {
4121 118880 : int red = con_is_redundant(tab, &tab->con[con]);
4122 118880 : if (red < 0)
4123 0 : goto error;
4124 118880 : if (!red)
4125 23181 : type = isl_ineq_cut;
4126 : else
4127 95699 : type = isl_ineq_redundant;
4128 : }
4129 :
4130 183111 : if (isl_tab_rollback(tab, snap))
4131 0 : return isl_ineq_error;
4132 183111 : return type;
4133 : error:
4134 0 : return isl_ineq_error;
4135 : }
4136 :
4137 3008625 : isl_stat isl_tab_track_bmap(struct isl_tab *tab, __isl_take isl_basic_map *bmap)
4138 : {
4139 3008625 : bmap = isl_basic_map_cow(bmap);
4140 3008625 : if (!tab || !bmap)
4141 : goto error;
4142 :
4143 3008625 : if (tab->empty) {
4144 2894196 : bmap = isl_basic_map_set_to_empty(bmap);
4145 2894196 : if (!bmap)
4146 0 : goto error;
4147 2894196 : tab->bmap = bmap;
4148 2894196 : return isl_stat_ok;
4149 : }
4150 :
4151 114429 : isl_assert(tab->mat->ctx, tab->n_eq == bmap->n_eq, goto error);
4152 114429 : isl_assert(tab->mat->ctx,
4153 : tab->n_con == bmap->n_eq + bmap->n_ineq, goto error);
4154 :
4155 114429 : tab->bmap = bmap;
4156 :
4157 114429 : return isl_stat_ok;
4158 : error:
4159 0 : isl_basic_map_free(bmap);
4160 0 : return isl_stat_error;
4161 : }
4162 :
4163 30253 : isl_stat isl_tab_track_bset(struct isl_tab *tab, __isl_take isl_basic_set *bset)
4164 : {
4165 30253 : return isl_tab_track_bmap(tab, bset_to_bmap(bset));
4166 : }
4167 :
4168 302001 : __isl_keep isl_basic_set *isl_tab_peek_bset(struct isl_tab *tab)
4169 : {
4170 302001 : if (!tab)
4171 0 : return NULL;
4172 :
4173 302001 : 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 : }
|
