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