/*

 * Copyright 2008-2009 Katholieke Universiteit Leuven

 * Copyright 2013      Ecole Normale Superieure

 *

 * Use of this software is governed by the MIT license

 *

 * Written by Sven Verdoolaege, K.U.Leuven, Departement

 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium

 * and Ecole Normale Superieure, 45 rue d'Ulm, 75230 Paris, France

 */

#include <isl_ctx_private.h>

#include <isl_mat_private.h>

#include <isl_vec_private.h>

#include "isl_map_private.h"

#include "isl_tab.h"

#include <isl_seq.h>

#include <isl_config.h>

/*

 * The implementation of tableaus in this file was inspired by Section 8

 * of David Detlefs, Greg Nelson and James B. Saxe, "Simplify: a theorem

 * prover for program checking".

 */

struct isl_tab *isl_tab_alloc(struct isl_ctx *ctx,

	unsigned n_row, unsigned n_var, unsigned M)

{

	int i;

	struct isl_tab *tab;

	unsigned off = 2 + M;

	tab = isl_calloc_type(ctx, struct isl_tab);

	if (!tab)

		return NULL;

	tab->mat = isl_mat_alloc(ctx, n_row, off + n_var);

	if (!tab->mat)

		goto error;

	tab->var = isl_alloc_array(ctx, struct isl_tab_var, n_var);

	if (n_var && !tab->var)

		goto error;

	tab->con = isl_alloc_array(ctx, struct isl_tab_var, n_row);

	if (n_row && !tab->con)

		goto error;

	tab->col_var = isl_alloc_array(ctx, int, n_var);

	if (n_var && !tab->col_var)

		goto error;

	tab->row_var = isl_alloc_array(ctx, int, n_row);

	if (n_row && !tab->row_var)

		goto error;

	for (i = 0; i < n_var; ++i) {

		tab->var[i].index = i;

		tab->var[i].is_row = 0;

		tab->var[i].is_nonneg = 0;

		tab->var[i].is_zero = 0;

		tab->var[i].is_redundant = 0;

		tab->var[i].frozen = 0;

		tab->var[i].negated = 0;

		tab->col_var[i] = i;

	}

	tab->n_row = 0;

	tab->n_con = 0;

	tab->n_eq = 0;

	tab->max_con = n_row;

	tab->n_col = n_var;

	tab->n_var = n_var;

	tab->max_var = n_var;

	tab->n_param = 0;

	tab->n_div = 0;

	tab->n_dead = 0;

	tab->n_redundant = 0;

	tab->strict_redundant = 0;

	tab->need_undo = 0;

	tab->rational = 0;

	tab->empty = 0;

	tab->in_undo = 0;

	tab->M = M;

	tab->cone = 0;

	tab->bottom.type = isl_tab_undo_bottom;

	tab->bottom.next = NULL;

	tab->top = &tab->bottom;

	tab->n_zero = 0;

	tab->n_unbounded = 0;

	tab->basis = NULL;

	return tab;

error:

	isl_tab_free(tab);

	return NULL;

}

isl_ctx *isl_tab_get_ctx(struct isl_tab *tab)

{

	return tab ? isl_mat_get_ctx(tab->mat) : NULL;

}

int isl_tab_extend_cons(struct isl_tab *tab, unsigned n_new)

{

	unsigned off;

	if (!tab)

		return -1;

	off = 2 + tab->M;

	if (tab->max_con < tab->n_con + n_new) {

		struct isl_tab_var *con;

		con = isl_realloc_array(tab->mat->ctx, tab->con,

				    struct isl_tab_var, tab->max_con + n_new);

		if (!con)

			return -1;

		tab->con = con;

		tab->max_con += n_new;

	}

	if (tab->mat->n_row < tab->n_row + n_new) {

		int *row_var;

		tab->mat = isl_mat_extend(tab->mat,

					tab->n_row + n_new, off + tab->n_col);

		if (!tab->mat)

			return -1;

		row_var = isl_realloc_array(tab->mat->ctx, tab->row_var,

					    int, tab->mat->n_row);

		if (!row_var)

			return -1;

		tab->row_var = row_var;

		if (tab->row_sign) {

			enum isl_tab_row_sign *s;

			s = isl_realloc_array(tab->mat->ctx, tab->row_sign,

					enum isl_tab_row_sign, tab->mat->n_row);

			if (!s)

				return -1;

			tab->row_sign = s;

		}

	}

	return 0;

}

/* Make room for at least n_new extra variables.

 * Return -1 if anything went wrong.

 */

int isl_tab_extend_vars(struct isl_tab *tab, unsigned n_new)

{

	struct isl_tab_var *var;

	unsigned off = 2 + tab->M;

	if (tab->max_var < tab->n_var + n_new) {

		var = isl_realloc_array(tab->mat->ctx, tab->var,

				    struct isl_tab_var, tab->n_var + n_new);

		if (!var)

			return -1;

		tab->var = var;

		tab->max_var += n_new;

	}

	if (tab->mat->n_col < off + tab->n_col + n_new) {

		int *p;

		tab->mat = isl_mat_extend(tab->mat,

				    tab->mat->n_row, off + tab->n_col + n_new);

		if (!tab->mat)

			return -1;

		p = isl_realloc_array(tab->mat->ctx, tab->col_var,

					    int, tab->n_col + n_new);

		if (!p)

			return -1;

		tab->col_var = p;

	}

	return 0;

}

static void free_undo_record(struct isl_tab_undo *undo)

{

	switch (undo->type) {

	case isl_tab_undo_saved_basis:

		free(undo->u.col_var);

		break;

	default:;

	}

	free(undo);

}

static void free_undo(struct isl_tab *tab)

{

	struct isl_tab_undo *undo, *next;

	for (undo = tab->top; undo && undo != &tab->bottom; undo = next) {

		next = undo->next;

		free_undo_record(undo);

	}

	tab->top = undo;

}

void isl_tab_free(struct isl_tab *tab)

{

	if (!tab)

		return;

	free_undo(tab);

	isl_mat_free(tab->mat);

	isl_vec_free(tab->dual);

	isl_basic_map_free(tab->bmap);

	free(tab->var);

	free(tab->con);

	free(tab->row_var);

	free(tab->col_var);

	free(tab->row_sign);

	isl_mat_free(tab->samples);

	free(tab->sample_index);

	isl_mat_free(tab->basis);

	free(tab);

}

struct isl_tab *isl_tab_dup(struct isl_tab *tab)

{

	int i;

	struct isl_tab *dup;

	unsigned off;

	if (!tab)

		return NULL;

	off = 2 + tab->M;

	dup = isl_calloc_type(tab->mat->ctx, struct isl_tab);

	if (!dup)

		return NULL;

	dup->mat = isl_mat_dup(tab->mat);

	if (!dup->mat)

		goto error;

	dup->var = isl_alloc_array(tab->mat->ctx, struct isl_tab_var, tab->max_var);

	if (tab->max_var && !dup->var)

		goto error;

	for (i = 0; i < tab->n_var; ++i)

		dup->var[i] = tab->var[i];

	dup->con = isl_alloc_array(tab->mat->ctx, struct isl_tab_var, tab->max_con);

	if (tab->max_con && !dup->con)

		goto error;

	for (i = 0; i < tab->n_con; ++i)

		dup->con[i] = tab->con[i];

	dup->col_var = isl_alloc_array(tab->mat->ctx, int, tab->mat->n_col - off);

	if ((tab->mat->n_col - off) && !dup->col_var)

		goto error;

	for (i = 0; i < tab->n_col; ++i)

		dup->col_var[i] = tab->col_var[i];

	dup->row_var = isl_alloc_array(tab->mat->ctx, int, tab->mat->n_row);

	if (tab->mat->n_row && !dup->row_var)

		goto error;

	for (i = 0; i < tab->n_row; ++i)

		dup->row_var[i] = tab->row_var[i];

	if (tab->row_sign) {

		dup->row_sign = isl_alloc_array(tab->mat->ctx, enum isl_tab_row_sign,

						tab->mat->n_row);

		if (tab->mat->n_row && !dup->row_sign)

			goto error;

		for (i = 0; i < tab->n_row; ++i)

			dup->row_sign[i] = tab->row_sign[i];

	}

	if (tab->samples) {

		dup->samples = isl_mat_dup(tab->samples);

		if (!dup->samples)

			goto error;

		dup->sample_index = isl_alloc_array(tab->mat->ctx, int,

							tab->samples->n_row);

		if (tab->samples->n_row && !dup->sample_index)

			goto error;

		dup->n_sample = tab->n_sample;

		dup->n_outside = tab->n_outside;

	}

	dup->n_row = tab->n_row;

	dup->n_con = tab->n_con;

	dup->n_eq = tab->n_eq;

	dup->max_con = tab->max_con;

	dup->n_col = tab->n_col;

	dup->n_var = tab->n_var;

	dup->max_var = tab->max_var;

	dup->n_param = tab->n_param;

	dup->n_div = tab->n_div;

	dup->n_dead = tab->n_dead;

	dup->n_redundant = tab->n_redundant;

	dup->rational = tab->rational;

	dup->empty = tab->empty;

	dup->strict_redundant = 0;

	dup->need_undo = 0;

	dup->in_undo = 0;

	dup->M = tab->M;

	tab->cone = tab->cone;

	dup->bottom.type = isl_tab_undo_bottom;

	dup->bottom.next = NULL;

	dup->top = &dup->bottom;

	dup->n_zero = tab->n_zero;

	dup->n_unbounded = tab->n_unbounded;

	dup->basis = isl_mat_dup(tab->basis);

	return dup;

error:

	isl_tab_free(dup);

	return NULL;

}

/* Construct the coefficient matrix of the product tableau

 * of two tableaus.

 * mat{1,2} is the coefficient matrix of tableau {1,2}

 * row{1,2} is the number of rows in tableau {1,2}

 * col{1,2} is the number of columns in tableau {1,2}

 * off is the offset to the coefficient column (skipping the

 *	denominator, the constant term and the big parameter if any)

 * r{1,2} is the number of redundant rows in tableau {1,2}

 * d{1,2} is the number of dead columns in tableau {1,2}

 *

 * The order of the rows and columns in the result is as explained

 * in isl_tab_product.

 */

static struct isl_mat *tab_mat_product(struct isl_mat *mat1,

	struct isl_mat *mat2, unsigned row1, unsigned row2,

	unsigned col1, unsigned col2,

	unsigned off, unsigned r1, unsigned r2, unsigned d1, unsigned d2)

{

	int i;

	struct isl_mat *prod;

	unsigned n;

	prod = isl_mat_alloc(mat1->ctx, mat1->n_row + mat2->n_row,

					off + col1 + col2);

	if (!prod)

		return NULL;

	n = 0;

	for (i = 0; i < r1; ++i) {

		isl_seq_cpy(prod->row[n + i], mat1->row[i], off + d1);

		isl_seq_clr(prod->row[n + i] + off + d1, d2);

		isl_seq_cpy(prod->row[n + i] + off + d1 + d2,

				mat1->row[i] + off + d1, col1 - d1);

		isl_seq_clr(prod->row[n + i] + off + col1 + d1, col2 - d2);

	}

	n += r1;

	for (i = 0; i < r2; ++i) {

		isl_seq_cpy(prod->row[n + i], mat2->row[i], off);

		isl_seq_clr(prod->row[n + i] + off, d1);

		isl_seq_cpy(prod->row[n + i] + off + d1,

			    mat2->row[i] + off, d2);

		isl_seq_clr(prod->row[n + i] + off + d1 + d2, col1 - d1);

		isl_seq_cpy(prod->row[n + i] + off + col1 + d1,

			    mat2->row[i] + off + d2, col2 - d2);

	}

	n += r2;

	for (i = 0; i < row1 - r1; ++i) {

		isl_seq_cpy(prod->row[n + i], mat1->row[r1 + i], off + d1);

		isl_seq_clr(prod->row[n + i] + off + d1, d2);

		isl_seq_cpy(prod->row[n + i] + off + d1 + d2,

				mat1->row[r1 + i] + off + d1, col1 - d1);

		isl_seq_clr(prod->row[n + i] + off + col1 + d1, col2 - d2);

	}

	n += row1 - r1;

	for (i = 0; i < row2 - r2; ++i) {

		isl_seq_cpy(prod->row[n + i], mat2->row[r2 + i], off);

		isl_seq_clr(prod->row[n + i] + off, d1);

		isl_seq_cpy(prod->row[n + i] + off + d1,

			    mat2->row[r2 + i] + off, d2);

		isl_seq_clr(prod->row[n + i] + off + d1 + d2, col1 - d1);

		isl_seq_cpy(prod->row[n + i] + off + col1 + d1,

			    mat2->row[r2 + i] + off + d2, col2 - d2);

	}

	return prod;

}

/* Update the row or column index of a variable that corresponds

 * to a variable in the first input tableau.

 */

static void update_index1(struct isl_tab_var *var,

	unsigned r1, unsigned r2, unsigned d1, unsigned d2)

{

	if (var->index == -1)

		return;

	if (var->is_row && var->index >= r1)

		var->index += r2;

	if (!var->is_row && var->index >= d1)

		var->index += d2;

}

/* Update the row or column index of a variable that corresponds

 * to a variable in the second input tableau.

 */

static void update_index2(struct isl_tab_var *var,

	unsigned row1, unsigned col1,

	unsigned r1, unsigned r2, unsigned d1, unsigned d2)

{

	if (var->index == -1)

		return;

	if (var->is_row) {

		if (var->index < r2)

			var->index += r1;

		else

			var->index += row1;

	} else {

		if (var->index < d2)

			var->index += d1;

		else

			var->index += col1;

	}

}

/* Create a tableau that represents the Cartesian product of the sets

 * represented by tableaus tab1 and tab2.

 * The order of the rows in the product is

 *	- redundant rows of tab1

 *	- redundant rows of tab2

 *	- non-redundant rows of tab1

 *	- non-redundant rows of tab2

 * The order of the columns is

 *	- denominator

 *	- constant term

 *	- coefficient of big parameter, if any

 *	- dead columns of tab1

 *	- dead columns of tab2

 *	- live columns of tab1

 *	- live columns of tab2

 * The order of the variables and the constraints is a concatenation

 * of order in the two input tableaus.

 */

struct isl_tab *isl_tab_product(struct isl_tab *tab1, struct isl_tab *tab2)

{

	int i;

	struct isl_tab *prod;

	unsigned off;

	unsigned r1, r2, d1, d2;

	if (!tab1 || !tab2)

		return NULL;

	isl_assert(tab1->mat->ctx, tab1->M == tab2->M, return NULL);

	isl_assert(tab1->mat->ctx, tab1->rational == tab2->rational, return NULL);

	isl_assert(tab1->mat->ctx, tab1->cone == tab2->cone, return NULL);

	isl_assert(tab1->mat->ctx, !tab1->row_sign, return NULL);

	isl_assert(tab1->mat->ctx, !tab2->row_sign, return NULL);

	isl_assert(tab1->mat->ctx, tab1->n_param == 0, return NULL);

	isl_assert(tab1->mat->ctx, tab2->n_param == 0, return NULL);

	isl_assert(tab1->mat->ctx, tab1->n_div == 0, return NULL);

	isl_assert(tab1->mat->ctx, tab2->n_div == 0, return NULL);

	off = 2 + tab1->M;

	r1 = tab1->n_redundant;

	r2 = tab2->n_redundant;

	d1 = tab1->n_dead;

	d2 = tab2->n_dead;

	prod = isl_calloc_type(tab1->mat->ctx, struct isl_tab);

	if (!prod)

		return NULL;

	prod->mat = tab_mat_product(tab1->mat, tab2->mat,

				tab1->n_row, tab2->n_row,

				tab1->n_col, tab2->n_col, off, r1, r2, d1, d2);

	if (!prod->mat)

		goto error;

	prod->var = isl_alloc_array(tab1->mat->ctx, struct isl_tab_var,

					tab1->max_var + tab2->max_var);

	if ((tab1->max_var + tab2->max_var) && !prod->var)

		goto error;

	for (i = 0; i < tab1->n_var; ++i) {

		prod->var[i] = tab1->var[i];

		update_index1(&prod->var[i], r1, r2, d1, d2);

	}

	for (i = 0; i < tab2->n_var; ++i) {

		prod->var[tab1->n_var + i] = tab2->var[i];

		update_index2(&prod->var[tab1->n_var + i],

				tab1->n_row, tab1->n_col,

				r1, r2, d1, d2);

	}

	prod->con = isl_alloc_array(tab1->mat->ctx, struct isl_tab_var,

					tab1->max_con +  tab2->max_con);

	if ((tab1->max_con + tab2->max_con) && !prod->con)

		goto error;

	for (i = 0; i < tab1->n_con; ++i) {

		prod->con[i] = tab1->con[i];

		update_index1(&prod->con[i], r1, r2, d1, d2);

	}

	for (i = 0; i < tab2->n_con; ++i) {

		prod->con[tab1->n_con + i] = tab2->con[i];

		update_index2(&prod->con[tab1->n_con + i],

				tab1->n_row, tab1->n_col,

				r1, r2, d1, d2);

	}

	prod->col_var = isl_alloc_array(tab1->mat->ctx, int,

					tab1->n_col + tab2->n_col);

	if ((tab1->n_col + tab2->n_col) && !prod->col_var)

		goto error;

	for (i = 0; i < tab1->n_col; ++i) {

		int pos = i < d1 ? i : i + d2;

		prod->col_var[pos] = tab1->col_var[i];

	}

	for (i = 0; i < tab2->n_col; ++i) {

		int pos = i < d2 ? d1 + i : tab1->n_col + i;

		int t = tab2->col_var[i];

		if (t >= 0)

			t += tab1->n_var;

		else

			t -= tab1->n_con;

		prod->col_var[pos] = t;

	}

	prod->row_var = isl_alloc_array(tab1->mat->ctx, int,

					tab1->mat->n_row + tab2->mat->n_row);

	if ((tab1->mat->n_row + tab2->mat->n_row) && !prod->row_var)

		goto error;

	for (i = 0; i < tab1->n_row; ++i) {

		int pos = i < r1 ? i : i + r2;

		prod->row_var[pos] = tab1->row_var[i];

	}

	for (i = 0; i < tab2->n_row; ++i) {

		int pos = i < r2 ? r1 + i : tab1->n_row + i;

		int t = tab2->row_var[i];

		if (t >= 0)

			t += tab1->n_var;

		else

			t -= tab1->n_con;

		prod->row_var[pos] = t;

	}

	prod->samples = NULL;

	prod->sample_index = NULL;

	prod->n_row = tab1->n_row + tab2->n_row;

	prod->n_con = tab1->n_con + tab2->n_con;

	prod->n_eq = 0;

	prod->max_con = tab1->max_con + tab2->max_con;

	prod->n_col = tab1->n_col + tab2->n_col;

	prod->n_var = tab1->n_var + tab2->n_var;

	prod->max_var = tab1->max_var + tab2->max_var;

	prod->n_param = 0;

	prod->n_div = 0;

	prod->n_dead = tab1->n_dead + tab2->n_dead;

	prod->n_redundant = tab1->n_redundant + tab2->n_redundant;

	prod->rational = tab1->rational;

	prod->empty = tab1->empty || tab2->empty;

	prod->strict_redundant = tab1->strict_redundant || tab2->strict_redundant;

	prod->need_undo = 0;

	prod->in_undo = 0;

	prod->M = tab1->M;

	prod->cone = tab1->cone;

	prod->bottom.type = isl_tab_undo_bottom;

	prod->bottom.next = NULL;

	prod->top = &prod->bottom;

	prod->n_zero = 0;

	prod->n_unbounded = 0;

	prod->basis = NULL;

	return prod;

error:

	isl_tab_free(prod);

	return NULL;

}

static struct isl_tab_var *var_from_index(struct isl_tab *tab, int i)

{

	if (i >= 0)

		return &tab->var[i];

	else

		return &tab->con[~i];

}

struct isl_tab_var *isl_tab_var_from_row(struct isl_tab *tab, int i)

{

	return var_from_index(tab, tab->row_var[i]);

}

static struct isl_tab_var *var_from_col(struct isl_tab *tab, int i)

{

	return var_from_index(tab, tab->col_var[i]);

}

/* Check if there are any upper bounds on column variable "var",

 * i.e., non-negative rows where var appears with a negative coefficient.

 * Return 1 if there are no such bounds.

 */

static int max_is_manifestly_unbounded(struct isl_tab *tab,

	struct isl_tab_var *var)

{

	int i;

	unsigned off = 2 + tab->M;

	if (var->is_row)

		return 0;

	for (i = tab->n_redundant; i < tab->n_row; ++i) {

		if (!isl_int_is_neg(tab->mat->row[i][off + var->index]))

			continue;

		if (isl_tab_var_from_row(tab, i)->is_nonneg)

			return 0;

	}

	return 1;

}

/* Check if there are any lower bounds on column variable "var",

 * i.e., non-negative rows where var appears with a positive coefficient.

 * Return 1 if there are no such bounds.

 */

static int min_is_manifestly_unbounded(struct isl_tab *tab,

	struct isl_tab_var *var)

{

	int i;

	unsigned off = 2 + tab->M;

	if (var->is_row)

		return 0;

	for (i = tab->n_redundant; i < tab->n_row; ++i) {

		if (!isl_int_is_pos(tab->mat->row[i][off + var->index]))

			continue;

		if (isl_tab_var_from_row(tab, i)->is_nonneg)

			return 0;

	}

	return 1;

}

static int row_cmp(struct isl_tab *tab, int r1, int r2, int c, isl_int t)

{

	unsigned off = 2 + tab->M;

	if (tab->M) {

		int s;

		isl_int_mul(t, tab->mat->row[r1][2], tab->mat->row[r2][off+c]);

		isl_int_submul(t, tab->mat->row[r2][2], tab->mat->row[r1][off+c]);

		s = isl_int_sgn(t);

		if (s)

			return s;

	}

	isl_int_mul(t, tab->mat->row[r1][1], tab->mat->row[r2][off + c]);

	isl_int_submul(t, tab->mat->row[r2][1], tab->mat->row[r1][off + c]);

	return isl_int_sgn(t);

}

/* Given the index of a column "c", return the index of a row

 * that can be used to pivot the column in, with either an increase

 * (sgn > 0) or a decrease (sgn < 0) of the corresponding variable.

 * If "var" is not NULL, then the row returned will be different from

 * the one associated with "var".

 *

 * Each row in the tableau is of the form

 *

 *	x_r = a_r0 + \sum_i a_ri x_i

 *

 * Only rows with x_r >= 0 and with the sign of a_ri opposite to "sgn"

 * impose any limit on the increase or decrease in the value of x_c

 * and this bound is equal to a_r0 / |a_rc|.  We are therefore looking

 * for the row with the smallest (most stringent) such bound.

 * Note that the common denominator of each row drops out of the fraction.

 * To check if row j has a smaller bound than row r, i.e.,

 * a_j0 / |a_jc| < a_r0 / |a_rc| or a_j0 |a_rc| < a_r0 |a_jc|,

 * we check if -sign(a_jc) (a_j0 a_rc - a_r0 a_jc) < 0,

 * where -sign(a_jc) is equal to "sgn".

 */

static int pivot_row(struct isl_tab *tab,

	struct isl_tab_var *var, int sgn, int c)

{

	int j, r, tsgn;

	isl_int t;

	unsigned off = 2 + tab->M;

	isl_int_init(t);

	r = -1;

	for (j = tab->n_redundant; j < tab->n_row; ++j) {

		if (var && j == var->index)

			continue;

		if (!isl_tab_var_from_row(tab, j)->is_nonneg)

			continue;

		if (sgn * isl_int_sgn(tab->mat->row[j][off + c]) >= 0)

			continue;

		if (r < 0) {

			r = j;

			continue;

		}

		tsgn = sgn * row_cmp(tab, r, j, c, t);

		if (tsgn < 0 || (tsgn == 0 &&

					    tab->row_var[j] < tab->row_var[r]))

			r = j;

	}

	isl_int_clear(t);

	return r;

}

/* Find a pivot (row and col) that will increase (sgn > 0) or decrease

 * (sgn < 0) the value of row variable var.

 * If not NULL, then skip_var is a row variable that should be ignored

 * while looking for a pivot row.  It is usually equal to var.

 *

 * As the given row in the tableau is of the form

 *

 *	x_r = a_r0 + \sum_i a_ri x_i

 *

 * we need to find a column such that the sign of a_ri is equal to "sgn"

 * (such that an increase in x_i will have the desired effect) or a

 * column with a variable that may attain negative values.

 * If a_ri is positive, then we need to move x_i in the same direction

 * to obtain the desired effect.  Otherwise, x_i has to move in the

 * opposite direction.

 */

static void find_pivot(struct isl_tab *tab,

	struct isl_tab_var *var, struct isl_tab_var *skip_var,

	int sgn, int *row, int *col)

{

	int j, r, c;

	isl_int *tr;

	*row = *col = -1;

	isl_assert(tab->mat->ctx, var->is_row, return);

	tr = tab->mat->row[var->index] + 2 + tab->M;

	c = -1;

	for (j = tab->n_dead; j < tab->n_col; ++j) {

		if (isl_int_is_zero(tr[j]))

			continue;

		if (isl_int_sgn(tr[j]) != sgn &&

		    var_from_col(tab, j)->is_nonneg)

			continue;

		if (c < 0 || tab->col_var[j] < tab->col_var[c])

			c = j;

	}

	if (c < 0)

		return;

	sgn *= isl_int_sgn(tr[c]);

	r = pivot_row(tab, skip_var, sgn, c);

	*row = r < 0 ? var->index : r;

	*col = c;

}

/* Return 1 if row "row" represents an obviously redundant inequality.

 * This means

 *	- it represents an inequality or a variable

 *	- that is the sum of a non-negative sample value and a positive

 *	  combination of zero or more non-negative constraints.

 */

int isl_tab_row_is_redundant(struct isl_tab *tab, int row)

{

	int i;

	unsigned off = 2 + tab->M;

	if (tab->row_var[row] < 0 && !isl_tab_var_from_row(tab, row)->is_nonneg)

		return 0;

	if (isl_int_is_neg(tab->mat->row[row][1]))

		return 0;

	if (tab->strict_redundant && isl_int_is_zero(tab->mat->row[row][1]))

		return 0;

	if (tab->M && isl_int_is_neg(tab->mat->row[row][2]))

		return 0;

	for (i = tab->n_dead; i < tab->n_col; ++i) {

		if (isl_int_is_zero(tab->mat->row[row][off + i]))

			continue;

		if (tab->col_var[i] >= 0)

			return 0;

		if (isl_int_is_neg(tab->mat->row[row][off + i]))

			return 0;

		if (!var_from_col(tab, i)->is_nonneg)

			return 0;

	}

	return 1;

}

static void swap_rows(struct isl_tab *tab, int row1, int row2)

{

	int t;

	enum isl_tab_row_sign s;

	t = tab->row_var[row1];

	tab->row_var[row1] = tab->row_var[row2];

	tab->row_var[row2] = t;

	isl_tab_var_from_row(tab, row1)->index = row1;

	isl_tab_var_from_row(tab, row2)->index = row2;

	tab->mat = isl_mat_swap_rows(tab->mat, row1, row2);

	if (!tab->row_sign)

		return;

	s = tab->row_sign[row1];

	tab->row_sign[row1] = tab->row_sign[row2];

	tab->row_sign[row2] = s;

}

static int push_union(struct isl_tab *tab,

	enum isl_tab_undo_type type, union isl_tab_undo_val u) WARN_UNUSED;

static int push_union(struct isl_tab *tab,

	enum isl_tab_undo_type type, union isl_tab_undo_val u)

{

	struct isl_tab_undo *undo;

	if (!tab)

		return -1;

	if (!tab->need_undo)

		return 0;

	undo = isl_alloc_type(tab->mat->ctx, struct isl_tab_undo);

	if (!undo)

		return -1;

	undo->type = type;

	undo->u = u;

	undo->next = tab->top;

	tab->top = undo;

	return 0;

}

int isl_tab_push_var(struct isl_tab *tab,

	enum isl_tab_undo_type type, struct isl_tab_var *var)

{

	union isl_tab_undo_val u;

	if (var->is_row)

		u.var_index = tab->row_var[var->index];

	else

		u.var_index = tab->col_var[var->index];

	return push_union(tab, type, u);

}

int isl_tab_push(struct isl_tab *tab, enum isl_tab_undo_type type)

{

	union isl_tab_undo_val u = { 0 };

	return push_union(tab, type, u);

}

/* Push a record on the undo stack describing the current basic

 * variables, so that the this state can be restored during rollback.

 */

int isl_tab_push_basis(struct isl_tab *tab)

{

	int i;

	union isl_tab_undo_val u;

	u.col_var = isl_alloc_array(tab->mat->ctx, int, tab->n_col);

	if (tab->n_col && !u.col_var)

		return -1;

	for (i = 0; i < tab->n_col; ++i)

		u.col_var[i] = tab->col_var[i];

	return push_union(tab, isl_tab_undo_saved_basis, u);

}

int isl_tab_push_callback(struct isl_tab *tab, struct isl_tab_callback *callback)

{

	union isl_tab_undo_val u;

	u.callback = callback;

	return push_union(tab, isl_tab_undo_callback, u);

}

struct isl_tab *isl_tab_init_samples(struct isl_tab *tab)

{

	if (!tab)

		return NULL;

	tab->n_sample = 0;

	tab->n_outside = 0;

	tab->samples = isl_mat_alloc(tab->mat->ctx, 1, 1 + tab->n_var);

	if (!tab->samples)

		goto error;

	tab->sample_index = isl_alloc_array(tab->mat->ctx, int, 1);

	if (!tab->sample_index)

		goto error;

	return tab;

error:

	isl_tab_free(tab);

	return NULL;

}

int isl_tab_add_sample(struct isl_tab *tab, __isl_take isl_vec *sample)

{

	if (!tab || !sample)

		goto error;

	if (tab->n_sample + 1 > tab->samples->n_row) {

		int *t = isl_realloc_array(tab->mat->ctx,

			    tab->sample_index, int, tab->n_sample + 1);

		if (!t)

			goto error;

		tab->sample_index = t;

	}

	tab->samples = isl_mat_extend(tab->samples,

				tab->n_sample + 1, tab->samples->n_col);

	if (!tab->samples)

		goto error;