brintos / llvm-project-archived public Read only

1/*2 * Copyright 2008-2009 Katholieke Universiteit Leuven3 * Copyright 2010      INRIA Saclay4 * Copyright 2016-2017 Sven Verdoolaege5 *6 * Use of this software is governed by the MIT license7 *8 * Written by Sven Verdoolaege, K.U.Leuven, Departement9 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium10 * and INRIA Saclay - Ile-de-France, Parc Club Orsay Universite,11 * ZAC des vignes, 4 rue Jacques Monod, 91893 Orsay, France 12 */13 14#include <isl_ctx_private.h>15#include "isl_map_private.h"16#include <isl_seq.h>17#include "isl_tab.h"18#include "isl_sample.h"19#include <isl_mat_private.h>20#include <isl_vec_private.h>21#include <isl_aff_private.h>22#include <isl_constraint_private.h>23#include <isl_options_private.h>24#include <isl_config.h>25 26#include <bset_to_bmap.c>27 28/*29 * The implementation of parametric integer linear programming in this file30 * was inspired by the paper "Parametric Integer Programming" and the31 * report "Solving systems of affine (in)equalities" by Paul Feautrier32 * (and others).33 *34 * The strategy used for obtaining a feasible solution is different35 * from the one used in isl_tab.c.  In particular, in isl_tab.c,36 * upon finding a constraint that is not yet satisfied, we pivot37 * in a row that increases the constant term of the row holding the38 * constraint, making sure the sample solution remains feasible39 * for all the constraints it already satisfied.40 * Here, we always pivot in the row holding the constraint,41 * choosing a column that induces the lexicographically smallest42 * increment to the sample solution.43 *44 * By starting out from a sample value that is lexicographically45 * smaller than any integer point in the problem space, the first46 * feasible integer sample point we find will also be the lexicographically47 * smallest.  If all variables can be assumed to be non-negative,48 * then the initial sample value may be chosen equal to zero.49 * However, we will not make this assumption.  Instead, we apply50 * the "big parameter" trick.  Any variable x is then not directly51 * used in the tableau, but instead it is represented by another52 * variable x' = M + x, where M is an arbitrarily large (positive)53 * value.  x' is therefore always non-negative, whatever the value of x.54 * Taking as initial sample value x' = 0 corresponds to x = -M,55 * which is always smaller than any possible value of x.56 *57 * The big parameter trick is used in the main tableau and58 * also in the context tableau if isl_context_lex is used.59 * In this case, each tableaus has its own big parameter.60 * Before doing any real work, we check if all the parameters61 * happen to be non-negative.  If so, we drop the column corresponding62 * to M from the initial context tableau.63 * If isl_context_gbr is used, then the big parameter trick is only64 * used in the main tableau.65 */66 67struct isl_context;68struct isl_context_op {69	/* detect nonnegative parameters in context and mark them in tab */70	struct isl_tab *(*detect_nonnegative_parameters)(71			struct isl_context *context, struct isl_tab *tab);72	/* return temporary reference to basic set representation of context */73	struct isl_basic_set *(*peek_basic_set)(struct isl_context *context);74	/* return temporary reference to tableau representation of context */75	struct isl_tab *(*peek_tab)(struct isl_context *context);76	/* add equality; check is 1 if eq may not be valid;77	 * update is 1 if we may want to call ineq_sign on context later.78	 */79	void (*add_eq)(struct isl_context *context, isl_int *eq,80			int check, int update);81	/* add inequality; check is 1 if ineq may not be valid;82	 * update is 1 if we may want to call ineq_sign on context later.83	 */84	void (*add_ineq)(struct isl_context *context, isl_int *ineq,85			int check, int update);86	/* check sign of ineq based on previous information.87	 * strict is 1 if saturation should be treated as a positive sign.88	 */89	enum isl_tab_row_sign (*ineq_sign)(struct isl_context *context,90			isl_int *ineq, int strict);91	/* check if inequality maintains feasibility */92	int (*test_ineq)(struct isl_context *context, isl_int *ineq);93	/* return index of a div that corresponds to "div" */94	int (*get_div)(struct isl_context *context, struct isl_tab *tab,95			struct isl_vec *div);96	/* insert div "div" to context at "pos" and return non-negativity */97	isl_bool (*insert_div)(struct isl_context *context, int pos,98		__isl_keep isl_vec *div);99	int (*detect_equalities)(struct isl_context *context,100			struct isl_tab *tab);101	/* return row index of "best" split */102	int (*best_split)(struct isl_context *context, struct isl_tab *tab);103	/* check if context has already been determined to be empty */104	int (*is_empty)(struct isl_context *context);105	/* check if context is still usable */106	int (*is_ok)(struct isl_context *context);107	/* save a copy/snapshot of context */108	void *(*save)(struct isl_context *context);109	/* restore saved context */110	void (*restore)(struct isl_context *context, void *);111	/* discard saved context */112	void (*discard)(void *);113	/* invalidate context */114	void (*invalidate)(struct isl_context *context);115	/* free context */116	__isl_null struct isl_context *(*free)(struct isl_context *context);117};118 119/* Shared parts of context representation.120 *121 * "n_unknown" is the number of final unknown integer divisions122 * in the input domain.123 */124struct isl_context {125	struct isl_context_op *op;126	int n_unknown;127};128 129struct isl_context_lex {130	struct isl_context context;131	struct isl_tab *tab;132};133 134/* A stack (linked list) of solutions of subtrees of the search space.135 *136 * "ma" describes the solution as a function of "dom".137 * In particular, the domain space of "ma" is equal to the space of "dom".138 *139 * If "ma" is NULL, then there is no solution on "dom".140 */141struct isl_partial_sol {142	int level;143	struct isl_basic_set *dom;144	isl_multi_aff *ma;145 146	struct isl_partial_sol *next;147};148 149struct isl_sol;150struct isl_sol_callback {151	struct isl_tab_callback callback;152	struct isl_sol *sol;153};154 155/* isl_sol is an interface for constructing a solution to156 * a parametric integer linear programming problem.157 * Every time the algorithm reaches a state where a solution158 * can be read off from the tableau, the function "add" is called159 * on the isl_sol passed to find_solutions_main.  In a state where160 * the tableau is empty, "add_empty" is called instead.161 * "free" is called to free the implementation specific fields, if any.162 *163 * "error" is set if some error has occurred.  This flag invalidates164 * the remainder of the data structure.165 * If "rational" is set, then a rational optimization is being performed.166 * "level" is the current level in the tree with nodes for each167 * split in the context.168 * If "max" is set, then a maximization problem is being solved, rather than169 * a minimization problem, which means that the variables in the170 * tableau have value "M - x" rather than "M + x".171 * "n_out" is the number of output dimensions in the input.172 * "space" is the space in which the solution (and also the input) lives.173 *174 * The context tableau is owned by isl_sol and is updated incrementally.175 *176 * There are currently two implementations of this interface,177 * isl_sol_map, which simply collects the solutions in an isl_map178 * and (optionally) the parts of the context where there is no solution179 * in an isl_set, and180 * isl_sol_pma, which collects an isl_pw_multi_aff instead.181 */182struct isl_sol {183	int error;184	int rational;185	int level;186	int max;187	isl_size n_out;188	isl_space *space;189	struct isl_context *context;190	struct isl_partial_sol *partial;191	void (*add)(struct isl_sol *sol,192		__isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma);193	void (*add_empty)(struct isl_sol *sol, struct isl_basic_set *bset);194	void (*free)(struct isl_sol *sol);195	struct isl_sol_callback	dec_level;196};197 198static void sol_free(struct isl_sol *sol)199{200	struct isl_partial_sol *partial, *next;201	if (!sol)202		return;203	for (partial = sol->partial; partial; partial = next) {204		next = partial->next;205		isl_basic_set_free(partial->dom);206		isl_multi_aff_free(partial->ma);207		free(partial);208	}209	isl_space_free(sol->space);210	if (sol->context)211		sol->context->op->free(sol->context);212	sol->free(sol);213	free(sol);214}215 216/* Push a partial solution represented by a domain and function "ma"217 * onto the stack of partial solutions.218 * If "ma" is NULL, then "dom" represents a part of the domain219 * with no solution.220 */221static void sol_push_sol(struct isl_sol *sol,222	__isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma)223{224	struct isl_partial_sol *partial;225 226	if (sol->error || !dom)227		goto error;228 229	partial = isl_alloc_type(dom->ctx, struct isl_partial_sol);230	if (!partial)231		goto error;232 233	partial->level = sol->level;234	partial->dom = dom;235	partial->ma = ma;236	partial->next = sol->partial;237 238	sol->partial = partial;239 240	return;241error:242	isl_basic_set_free(dom);243	isl_multi_aff_free(ma);244	sol->error = 1;245}246 247/* Check that the final columns of "M", starting at "first", are zero.248 */249static isl_stat check_final_columns_are_zero(__isl_keep isl_mat *M,250	unsigned first)251{252	int i;253	isl_size rows, cols;254	unsigned n;255 256	rows = isl_mat_rows(M);257	cols = isl_mat_cols(M);258	if (rows < 0 || cols < 0)259		return isl_stat_error;260	n = cols - first;261	for (i = 0; i < rows; ++i)262		if (isl_seq_first_non_zero(M->row[i] + first, n) != -1)263			isl_die(isl_mat_get_ctx(M), isl_error_internal,264				"final columns should be zero",265				return isl_stat_error);266	return isl_stat_ok;267}268 269/* Set the affine expressions in "ma" according to the rows in "M", which270 * are defined over the local space "ls".271 * The matrix "M" may have extra (zero) columns beyond the number272 * of variables in "ls".273 */274static __isl_give isl_multi_aff *set_from_affine_matrix(275	__isl_take isl_multi_aff *ma, __isl_take isl_local_space *ls,276	__isl_take isl_mat *M)277{278	int i;279	isl_size dim;280	isl_aff *aff;281 282	dim = isl_local_space_dim(ls, isl_dim_all);283	if (!ma || dim < 0 || !M)284		goto error;285 286	if (check_final_columns_are_zero(M, 1 + dim) < 0)287		goto error;288	for (i = 1; i < M->n_row; ++i) {289		aff = isl_aff_alloc(isl_local_space_copy(ls));290		if (aff) {291			isl_int_set(aff->v->el[0], M->row[0][0]);292			isl_seq_cpy(aff->v->el + 1, M->row[i], 1 + dim);293		}294		aff = isl_aff_normalize(aff);295		ma = isl_multi_aff_set_aff(ma, i - 1, aff);296	}297	isl_local_space_free(ls);298	isl_mat_free(M);299 300	return ma;301error:302	isl_local_space_free(ls);303	isl_mat_free(M);304	isl_multi_aff_free(ma);305	return NULL;306}307 308/* Push a partial solution represented by a domain and mapping M309 * onto the stack of partial solutions.310 *311 * The affine matrix "M" maps the dimensions of the context312 * to the output variables.  Convert it into an isl_multi_aff and313 * then call sol_push_sol.314 *315 * Note that the description of the initial context may have involved316 * existentially quantified variables, in which case they also appear317 * in "dom".  These need to be removed before creating the affine318 * expression because an affine expression cannot be defined in terms319 * of existentially quantified variables without a known representation.320 * Since newly added integer divisions are inserted before these321 * existentially quantified variables, they are still in the final322 * positions and the corresponding final columns of "M" are zero323 * because align_context_divs adds the existentially quantified324 * variables of the context to the main tableau without any constraints and325 * any equality constraints that are added later on can only serve326 * to eliminate these existentially quantified variables.327 */328static void sol_push_sol_mat(struct isl_sol *sol,329	__isl_take isl_basic_set *dom, __isl_take isl_mat *M)330{331	isl_local_space *ls;332	isl_multi_aff *ma;333	isl_size n_div;334	int n_known;335 336	n_div = isl_basic_set_dim(dom, isl_dim_div);337	if (n_div < 0)338		goto error;339	n_known = n_div - sol->context->n_unknown;340 341	ma = isl_multi_aff_alloc(isl_space_copy(sol->space));342	ls = isl_basic_set_get_local_space(dom);343	ls = isl_local_space_drop_dims(ls, isl_dim_div,344					n_known, n_div - n_known);345	ma = set_from_affine_matrix(ma, ls, M);346 347	if (!ma)348		dom = isl_basic_set_free(dom);349	sol_push_sol(sol, dom, ma);350	return;351error:352	isl_basic_set_free(dom);353	isl_mat_free(M);354	sol_push_sol(sol, NULL, NULL);355}356 357/* Pop one partial solution from the partial solution stack and358 * pass it on to sol->add or sol->add_empty.359 */360static void sol_pop_one(struct isl_sol *sol)361{362	struct isl_partial_sol *partial;363 364	partial = sol->partial;365	sol->partial = partial->next;366 367	if (partial->ma)368		sol->add(sol, partial->dom, partial->ma);369	else370		sol->add_empty(sol, partial->dom);371	free(partial);372}373 374/* Return a fresh copy of the domain represented by the context tableau.375 */376static struct isl_basic_set *sol_domain(struct isl_sol *sol)377{378	struct isl_basic_set *bset;379 380	if (sol->error)381		return NULL;382 383	bset = isl_basic_set_dup(sol->context->op->peek_basic_set(sol->context));384	bset = isl_basic_set_update_from_tab(bset,385			sol->context->op->peek_tab(sol->context));386 387	return bset;388}389 390/* Check whether two partial solutions have the same affine expressions.391 */392static isl_bool same_solution(struct isl_partial_sol *s1,393	struct isl_partial_sol *s2)394{395	if (!s1->ma != !s2->ma)396		return isl_bool_false;397	if (!s1->ma)398		return isl_bool_true;399 400	return isl_multi_aff_plain_is_equal(s1->ma, s2->ma);401}402 403/* Swap the initial two partial solutions in "sol".404 *405 * That is, go from406 *407 *	sol->partial = p1; p1->next = p2; p2->next = p3408 *409 * to410 *411 *	sol->partial = p2; p2->next = p1; p1->next = p3412 */413static void swap_initial(struct isl_sol *sol)414{415	struct isl_partial_sol *partial;416 417	partial = sol->partial;418	sol->partial = partial->next;419	partial->next = partial->next->next;420	sol->partial->next = partial;421}422 423/* Combine the initial two partial solution of "sol" into424 * a partial solution with the current context domain of "sol" and425 * the function description of the second partial solution in the list.426 * The level of the new partial solution is set to the current level.427 *428 * That is, the first two partial solutions (D1,M1) and (D2,M2) are429 * replaced by (D,M2), where D is the domain of "sol", which is assumed430 * to be the union of D1 and D2, while M1 is assumed to be equal to M2431 * (at least on D1).432 */433static isl_stat combine_initial_into_second(struct isl_sol *sol)434{435	struct isl_partial_sol *partial;436	isl_basic_set *bset;437 438	partial = sol->partial;439 440	bset = sol_domain(sol);441	isl_basic_set_free(partial->next->dom);442	partial->next->dom = bset;443	partial->next->level = sol->level;444 445	if (!bset)446		return isl_stat_error;447 448	sol->partial = partial->next;449	isl_basic_set_free(partial->dom);450	isl_multi_aff_free(partial->ma);451	free(partial);452 453	return isl_stat_ok;454}455 456/* Are "ma1" and "ma2" equal to each other on "dom"?457 *458 * Combine "ma1" and "ma2" with "dom" and check if the results are the same.459 * "dom" may have existentially quantified variables.  Eliminate them first460 * as otherwise they would have to be eliminated twice, in a more complicated461 * context.462 */463static isl_bool equal_on_domain(__isl_keep isl_multi_aff *ma1,464	__isl_keep isl_multi_aff *ma2, __isl_keep isl_basic_set *dom)465{466	isl_set *set;467	isl_pw_multi_aff *pma1, *pma2;468	isl_bool equal;469 470	set = isl_basic_set_compute_divs(isl_basic_set_copy(dom));471	pma1 = isl_pw_multi_aff_alloc(isl_set_copy(set),472					isl_multi_aff_copy(ma1));473	pma2 = isl_pw_multi_aff_alloc(set, isl_multi_aff_copy(ma2));474	equal = isl_pw_multi_aff_is_equal(pma1, pma2);475	isl_pw_multi_aff_free(pma1);476	isl_pw_multi_aff_free(pma2);477 478	return equal;479}480 481/* The initial two partial solutions of "sol" are known to be at482 * the same level.483 * If they represent the same solution (on different parts of the domain),484 * then combine them into a single solution at the current level.485 * Otherwise, pop them both.486 *487 * Even if the two partial solution are not obviously the same,488 * one may still be a simplification of the other over its own domain.489 * Also check if the two sets of affine functions are equal when490 * restricted to one of the domains.  If so, combine the two491 * using the set of affine functions on the other domain.492 * That is, for two partial solutions (D1,M1) and (D2,M2),493 * if M1 = M2 on D1, then the pair of partial solutions can494 * be replaced by (D1+D2,M2) and similarly when M1 = M2 on D2.495 */496static isl_stat combine_initial_if_equal(struct isl_sol *sol)497{498	struct isl_partial_sol *partial;499	isl_bool same;500 501	partial = sol->partial;502 503	same = same_solution(partial, partial->next);504	if (same < 0)505		return isl_stat_error;506	if (same)507		return combine_initial_into_second(sol);508	if (partial->ma && partial->next->ma) {509		same = equal_on_domain(partial->ma, partial->next->ma,510					partial->dom);511		if (same < 0)512			return isl_stat_error;513		if (same)514			return combine_initial_into_second(sol);515		same = equal_on_domain(partial->ma, partial->next->ma,516					partial->next->dom);517		if (same) {518			swap_initial(sol);519			return combine_initial_into_second(sol);520		}521	}522 523	sol_pop_one(sol);524	sol_pop_one(sol);525 526	return isl_stat_ok;527}528 529/* Pop all solutions from the partial solution stack that were pushed onto530 * the stack at levels that are deeper than the current level.531 * If the two topmost elements on the stack have the same level532 * and represent the same solution, then their domains are combined.533 * This combined domain is the same as the current context domain534 * as sol_pop is called each time we move back to a higher level.535 * If the outer level (0) has been reached, then all partial solutions536 * at the current level are also popped off.537 */538static void sol_pop(struct isl_sol *sol)539{540	struct isl_partial_sol *partial;541 542	if (sol->error)543		return;544 545	partial = sol->partial;546	if (!partial)547		return;548 549	if (partial->level == 0 && sol->level == 0) {550		for (partial = sol->partial; partial; partial = sol->partial)551			sol_pop_one(sol);552		return;553	}554 555	if (partial->level <= sol->level)556		return;557 558	if (partial->next && partial->next->level == partial->level) {559		if (combine_initial_if_equal(sol) < 0)560			goto error;561	} else562		sol_pop_one(sol);563 564	if (sol->level == 0) {565		for (partial = sol->partial; partial; partial = sol->partial)566			sol_pop_one(sol);567		return;568	}569 570	if (0)571error:		sol->error = 1;572}573 574static void sol_dec_level(struct isl_sol *sol)575{576	if (sol->error)577		return;578 579	sol->level--;580 581	sol_pop(sol);582}583 584static isl_stat sol_dec_level_wrap(struct isl_tab_callback *cb)585{586	struct isl_sol_callback *callback = (struct isl_sol_callback *)cb;587 588	sol_dec_level(callback->sol);589 590	return callback->sol->error ? isl_stat_error : isl_stat_ok;591}592 593/* Move down to next level and push callback onto context tableau594 * to decrease the level again when it gets rolled back across595 * the current state.  That is, dec_level will be called with596 * the context tableau in the same state as it is when inc_level597 * is called.598 */599static void sol_inc_level(struct isl_sol *sol)600{601	struct isl_tab *tab;602 603	if (sol->error)604		return;605 606	sol->level++;607	tab = sol->context->op->peek_tab(sol->context);608	if (isl_tab_push_callback(tab, &sol->dec_level.callback) < 0)609		sol->error = 1;610}611 612static void scale_rows(struct isl_mat *mat, isl_int m, int n_row)613{614	int i;615 616	if (isl_int_is_one(m))617		return;618 619	for (i = 0; i < n_row; ++i)620		isl_seq_scale(mat->row[i], mat->row[i], m, mat->n_col);621}622 623/* Add the solution identified by the tableau and the context tableau.624 *625 * The layout of the variables is as follows.626 *	tab->n_var is equal to the total number of variables in the input627 *			map (including divs that were copied from the context)628 *			+ the number of extra divs constructed629 *      Of these, the first tab->n_param and the last tab->n_div variables630 *	correspond to the variables in the context, i.e.,631 *		tab->n_param + tab->n_div = context_tab->n_var632 *	tab->n_param is equal to the number of parameters and input633 *			dimensions in the input map634 *	tab->n_div is equal to the number of divs in the context635 *636 * If there is no solution, then call add_empty with a basic set637 * that corresponds to the context tableau.  (If add_empty is NULL,638 * then do nothing).639 *640 * If there is a solution, then first construct a matrix that maps641 * all dimensions of the context to the output variables, i.e.,642 * the output dimensions in the input map.643 * The divs in the input map (if any) that do not correspond to any644 * div in the context do not appear in the solution.645 * The algorithm will make sure that they have an integer value,646 * but these values themselves are of no interest.647 * We have to be careful not to drop or rearrange any divs in the648 * context because that would change the meaning of the matrix.649 *650 * To extract the value of the output variables, it should be noted651 * that we always use a big parameter M in the main tableau and so652 * the variable stored in this tableau is not an output variable x itself, but653 *	x' = M + x (in case of minimization)654 * or655 *	x' = M - x (in case of maximization)656 * If x' appears in a column, then its optimal value is zero,657 * which means that the optimal value of x is an unbounded number658 * (-M for minimization and M for maximization).659 * We currently assume that the output dimensions in the original map660 * are bounded, so this cannot occur.661 * Similarly, when x' appears in a row, then the coefficient of M in that662 * row is necessarily 1.663 * If the row in the tableau represents664 *	d x' = c + d M + e(y)665 * then, in case of minimization, the corresponding row in the matrix666 * will be667 *	a c + a e(y)668 * with a d = m, the (updated) common denominator of the matrix.669 * In case of maximization, the row will be670 *	-a c - a e(y)671 */672static void sol_add(struct isl_sol *sol, struct isl_tab *tab)673{674	struct isl_basic_set *bset = NULL;675	struct isl_mat *mat = NULL;676	unsigned off;677	int row;678	isl_int m;679 680	if (sol->error || !tab)681		goto error;682 683	if (tab->empty && !sol->add_empty)684		return;685	if (sol->context->op->is_empty(sol->context))686		return;687 688	bset = sol_domain(sol);689 690	if (tab->empty) {691		sol_push_sol(sol, bset, NULL);692		return;693	}694 695	off = 2 + tab->M;696 697	mat = isl_mat_alloc(tab->mat->ctx, 1 + sol->n_out,698					    1 + tab->n_param + tab->n_div);699	if (!mat)700		goto error;701 702	isl_int_init(m);703 704	isl_seq_clr(mat->row[0] + 1, mat->n_col - 1);705	isl_int_set_si(mat->row[0][0], 1);706	for (row = 0; row < sol->n_out; ++row) {707		int i = tab->n_param + row;708		int r, j;709 710		isl_seq_clr(mat->row[1 + row], mat->n_col);711		if (!tab->var[i].is_row) {712			if (tab->M)713				isl_die(mat->ctx, isl_error_invalid,714					"unbounded optimum", goto error2);715			continue;716		}717 718		r = tab->var[i].index;719		if (tab->M &&720		    isl_int_ne(tab->mat->row[r][2], tab->mat->row[r][0]))721			isl_die(mat->ctx, isl_error_invalid,722				"unbounded optimum", goto error2);723		isl_int_gcd(m, mat->row[0][0], tab->mat->row[r][0]);724		isl_int_divexact(m, tab->mat->row[r][0], m);725		scale_rows(mat, m, 1 + row);726		isl_int_divexact(m, mat->row[0][0], tab->mat->row[r][0]);727		isl_int_mul(mat->row[1 + row][0], m, tab->mat->row[r][1]);728		for (j = 0; j < tab->n_param; ++j) {729			int col;730			if (tab->var[j].is_row)731				continue;732			col = tab->var[j].index;733			isl_int_mul(mat->row[1 + row][1 + j], m,734				    tab->mat->row[r][off + col]);735		}736		for (j = 0; j < tab->n_div; ++j) {737			int col;738			if (tab->var[tab->n_var - tab->n_div+j].is_row)739				continue;740			col = tab->var[tab->n_var - tab->n_div+j].index;741			isl_int_mul(mat->row[1 + row][1 + tab->n_param + j], m,742				    tab->mat->row[r][off + col]);743		}744		if (sol->max)745			isl_seq_neg(mat->row[1 + row], mat->row[1 + row],746				    mat->n_col);747	}748 749	isl_int_clear(m);750 751	sol_push_sol_mat(sol, bset, mat);752	return;753error2:754	isl_int_clear(m);755error:756	isl_basic_set_free(bset);757	isl_mat_free(mat);758	sol->error = 1;759}760 761struct isl_sol_map {762	struct isl_sol	sol;763	struct isl_map	*map;764	struct isl_set	*empty;765};766 767static void sol_map_free(struct isl_sol *sol)768{769	struct isl_sol_map *sol_map = (struct isl_sol_map *) sol;770	isl_map_free(sol_map->map);771	isl_set_free(sol_map->empty);772}773 774/* This function is called for parts of the context where there is775 * no solution, with "bset" corresponding to the context tableau.776 * Simply add the basic set to the set "empty".777 */778static void sol_map_add_empty(struct isl_sol_map *sol,779	struct isl_basic_set *bset)780{781	if (!bset || !sol->empty)782		goto error;783 784	sol->empty = isl_set_grow(sol->empty, 1);785	bset = isl_basic_set_simplify(bset);786	bset = isl_basic_set_finalize(bset);787	sol->empty = isl_set_add_basic_set(sol->empty, isl_basic_set_copy(bset));788	if (!sol->empty)789		goto error;790	isl_basic_set_free(bset);791	return;792error:793	isl_basic_set_free(bset);794	sol->sol.error = 1;795}796 797static void sol_map_add_empty_wrap(struct isl_sol *sol,798	struct isl_basic_set *bset)799{800	sol_map_add_empty((struct isl_sol_map *)sol, bset);801}802 803/* Given a basic set "dom" that represents the context and a tuple of804 * affine expressions "ma" defined over this domain, construct a basic map805 * that expresses this function on the domain.806 */807static void sol_map_add(struct isl_sol_map *sol,808	__isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma)809{810	isl_basic_map *bmap;811 812	if (sol->sol.error || !dom || !ma)813		goto error;814 815	bmap = isl_basic_map_from_multi_aff2(ma, sol->sol.rational);816	bmap = isl_basic_map_intersect_domain(bmap, dom);817	sol->map = isl_map_grow(sol->map, 1);818	sol->map = isl_map_add_basic_map(sol->map, bmap);819	if (!sol->map)820		sol->sol.error = 1;821	return;822error:823	isl_basic_set_free(dom);824	isl_multi_aff_free(ma);825	sol->sol.error = 1;826}827 828static void sol_map_add_wrap(struct isl_sol *sol,829	__isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma)830{831	sol_map_add((struct isl_sol_map *)sol, dom, ma);832}833 834 835/* Store the "parametric constant" of row "row" of tableau "tab" in "line",836 * i.e., the constant term and the coefficients of all variables that837 * appear in the context tableau.838 * Note that the coefficient of the big parameter M is NOT copied.839 * The context tableau may not have a big parameter and even when it840 * does, it is a different big parameter.841 */842static void get_row_parameter_line(struct isl_tab *tab, int row, isl_int *line)843{844	int i;845	unsigned off = 2 + tab->M;846 847	isl_int_set(line[0], tab->mat->row[row][1]);848	for (i = 0; i < tab->n_param; ++i) {849		if (tab->var[i].is_row)850			isl_int_set_si(line[1 + i], 0);851		else {852			int col = tab->var[i].index;853			isl_int_set(line[1 + i], tab->mat->row[row][off + col]);854		}855	}856	for (i = 0; i < tab->n_div; ++i) {857		if (tab->var[tab->n_var - tab->n_div + i].is_row)858			isl_int_set_si(line[1 + tab->n_param + i], 0);859		else {860			int col = tab->var[tab->n_var - tab->n_div + i].index;861			isl_int_set(line[1 + tab->n_param + i],862				    tab->mat->row[row][off + col]);863		}864	}865}866 867/* Check if rows "row1" and "row2" have identical "parametric constants",868 * as explained above.869 * In this case, we also insist that the coefficients of the big parameter870 * be the same as the values of the constants will only be the same871 * if these coefficients are also the same.872 */873static int identical_parameter_line(struct isl_tab *tab, int row1, int row2)874{875	int i;876	unsigned off = 2 + tab->M;877 878	if (isl_int_ne(tab->mat->row[row1][1], tab->mat->row[row2][1]))879		return 0;880 881	if (tab->M && isl_int_ne(tab->mat->row[row1][2],882				 tab->mat->row[row2][2]))883		return 0;884 885	for (i = 0; i < tab->n_param + tab->n_div; ++i) {886		int pos = i < tab->n_param ? i :887			tab->n_var - tab->n_div + i - tab->n_param;888		int col;889 890		if (tab->var[pos].is_row)891			continue;892		col = tab->var[pos].index;893		if (isl_int_ne(tab->mat->row[row1][off + col],894			       tab->mat->row[row2][off + col]))895			return 0;896	}897	return 1;898}899 900/* Return an inequality that expresses that the "parametric constant"901 * should be non-negative.902 * This function is only called when the coefficient of the big parameter903 * is equal to zero.904 */905static struct isl_vec *get_row_parameter_ineq(struct isl_tab *tab, int row)906{907	struct isl_vec *ineq;908 909	ineq = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_param + tab->n_div);910	if (!ineq)911		return NULL;912 913	get_row_parameter_line(tab, row, ineq->el);914	if (ineq)915		ineq = isl_vec_normalize(ineq);916 917	return ineq;918}919 920/* Normalize a div expression of the form921 *922 *	[(g*f(x) + c)/(g * m)]923 *924 * with c the constant term and f(x) the remaining coefficients, to925 *926 *	[(f(x) + [c/g])/m]927 */928static void normalize_div(__isl_keep isl_vec *div)929{930	isl_ctx *ctx = isl_vec_get_ctx(div);931	int len = div->size - 2;932 933	isl_seq_gcd(div->el + 2, len, &ctx->normalize_gcd);934	isl_int_gcd(ctx->normalize_gcd, ctx->normalize_gcd, div->el[0]);935 936	if (isl_int_is_one(ctx->normalize_gcd))937		return;938 939	isl_int_divexact(div->el[0], div->el[0], ctx->normalize_gcd);940	isl_int_fdiv_q(div->el[1], div->el[1], ctx->normalize_gcd);941	isl_seq_scale_down(div->el + 2, div->el + 2, ctx->normalize_gcd, len);942}943 944/* Return an integer division for use in a parametric cut based945 * on the given row.946 * In particular, let the parametric constant of the row be947 *948 *		\sum_i a_i y_i949 *950 * where y_0 = 1, but none of the y_i corresponds to the big parameter M.951 * The div returned is equal to952 *953 *		floor(\sum_i {-a_i} y_i) = floor((\sum_i (-a_i mod d) y_i)/d)954 */955static struct isl_vec *get_row_parameter_div(struct isl_tab *tab, int row)956{957	struct isl_vec *div;958 959	div = isl_vec_alloc(tab->mat->ctx, 1 + 1 + tab->n_param + tab->n_div);960	if (!div)961		return NULL;962 963	isl_int_set(div->el[0], tab->mat->row[row][0]);964	get_row_parameter_line(tab, row, div->el + 1);965	isl_seq_neg(div->el + 1, div->el + 1, div->size - 1);966	normalize_div(div);967	isl_seq_fdiv_r(div->el + 1, div->el + 1, div->el[0], div->size - 1);968 969	return div;970}971 972/* Return an integer division for use in transferring an integrality constraint973 * to the context.974 * In particular, let the parametric constant of the row be975 *976 *		\sum_i a_i y_i977 *978 * where y_0 = 1, but none of the y_i corresponds to the big parameter M.979 * The the returned div is equal to980 *981 *		floor(\sum_i {a_i} y_i) = floor((\sum_i (a_i mod d) y_i)/d)982 */983static struct isl_vec *get_row_split_div(struct isl_tab *tab, int row)984{985	struct isl_vec *div;986 987	div = isl_vec_alloc(tab->mat->ctx, 1 + 1 + tab->n_param + tab->n_div);988	if (!div)989		return NULL;990 991	isl_int_set(div->el[0], tab->mat->row[row][0]);992	get_row_parameter_line(tab, row, div->el + 1);993	normalize_div(div);994	isl_seq_fdiv_r(div->el + 1, div->el + 1, div->el[0], div->size - 1);995 996	return div;997}998 999/* Construct and return an inequality that expresses an upper bound1000 * on the given div.1001 * In particular, if the div is given by1002 *1003 *	d = floor(e/m)1004 *1005 * then the inequality expresses1006 *1007 *	m d <= e1008 */1009static __isl_give isl_vec *ineq_for_div(__isl_keep isl_basic_set *bset,1010	unsigned div)1011{1012	isl_size total;1013	unsigned div_pos;1014	struct isl_vec *ineq;1015 1016	total = isl_basic_set_dim(bset, isl_dim_all);1017	if (total < 0)1018		return NULL;1019 1020	div_pos = 1 + total - bset->n_div + div;1021 1022	ineq = isl_vec_alloc(bset->ctx, 1 + total);1023	if (!ineq)1024		return NULL;1025 1026	isl_seq_cpy(ineq->el, bset->div[div] + 1, 1 + total);1027	isl_int_neg(ineq->el[div_pos], bset->div[div][0]);1028	return ineq;1029}1030 1031/* Given a row in the tableau and a div that was created1032 * using get_row_split_div and that has been constrained to equality, i.e.,1033 *1034 *		d = floor(\sum_i {a_i} y_i) = \sum_i {a_i} y_i1035 *1036 * replace the expression "\sum_i {a_i} y_i" in the row by d,1037 * i.e., we subtract "\sum_i {a_i} y_i" and add 1 d.1038 * The coefficients of the non-parameters in the tableau have been1039 * verified to be integral.  We can therefore simply replace coefficient b1040 * by floor(b).  For the coefficients of the parameters we have1041 * floor(a_i) = a_i - {a_i}, while for the other coefficients, we have1042 * floor(b) = b.1043 */1044static struct isl_tab *set_row_cst_to_div(struct isl_tab *tab, int row, int div)1045{1046	isl_seq_fdiv_q(tab->mat->row[row] + 1, tab->mat->row[row] + 1,1047			tab->mat->row[row][0], 1 + tab->M + tab->n_col);1048 1049	isl_int_set_si(tab->mat->row[row][0], 1);1050 1051	if (tab->var[tab->n_var - tab->n_div + div].is_row) {1052		int drow = tab->var[tab->n_var - tab->n_div + div].index;1053 1054		isl_assert(tab->mat->ctx,1055			isl_int_is_one(tab->mat->row[drow][0]), goto error);1056		isl_seq_combine(tab->mat->row[row] + 1,1057			tab->mat->ctx->one, tab->mat->row[row] + 1,1058			tab->mat->ctx->one, tab->mat->row[drow] + 1,1059			1 + tab->M + tab->n_col);1060	} else {1061		int dcol = tab->var[tab->n_var - tab->n_div + div].index;1062 1063		isl_int_add_ui(tab->mat->row[row][2 + tab->M + dcol],1064				tab->mat->row[row][2 + tab->M + dcol], 1);1065	}1066 1067	return tab;1068error:1069	isl_tab_free(tab);1070	return NULL;1071}1072 1073/* Check if the (parametric) constant of the given row is obviously1074 * negative, meaning that we don't need to consult the context tableau.1075 * If there is a big parameter and its coefficient is non-zero,1076 * then this coefficient determines the outcome.1077 * Otherwise, we check whether the constant is negative and1078 * all non-zero coefficients of parameters are negative and1079 * belong to non-negative parameters.1080 */1081static int is_obviously_neg(struct isl_tab *tab, int row)1082{1083	int i;1084	int col;1085	unsigned off = 2 + tab->M;1086 1087	if (tab->M) {1088		if (isl_int_is_pos(tab->mat->row[row][2]))1089			return 0;1090		if (isl_int_is_neg(tab->mat->row[row][2]))1091			return 1;1092	}1093 1094	if (isl_int_is_nonneg(tab->mat->row[row][1]))1095		return 0;1096	for (i = 0; i < tab->n_param; ++i) {1097		/* Eliminated parameter */1098		if (tab->var[i].is_row)1099			continue;1100		col = tab->var[i].index;1101		if (isl_int_is_zero(tab->mat->row[row][off + col]))1102			continue;1103		if (!tab->var[i].is_nonneg)1104			return 0;1105		if (isl_int_is_pos(tab->mat->row[row][off + col]))1106			return 0;1107	}1108	for (i = 0; i < tab->n_div; ++i) {1109		if (tab->var[tab->n_var - tab->n_div + i].is_row)1110			continue;1111		col = tab->var[tab->n_var - tab->n_div + i].index;1112		if (isl_int_is_zero(tab->mat->row[row][off + col]))1113			continue;1114		if (!tab->var[tab->n_var - tab->n_div + i].is_nonneg)1115			return 0;1116		if (isl_int_is_pos(tab->mat->row[row][off + col]))1117			return 0;1118	}1119	return 1;1120}1121 1122/* Check if the (parametric) constant of the given row is obviously1123 * non-negative, meaning that we don't need to consult the context tableau.1124 * If there is a big parameter and its coefficient is non-zero,1125 * then this coefficient determines the outcome.1126 * Otherwise, we check whether the constant is non-negative and1127 * all non-zero coefficients of parameters are positive and1128 * belong to non-negative parameters.1129 */1130static int is_obviously_nonneg(struct isl_tab *tab, int row)1131{1132	int i;1133	int col;1134	unsigned off = 2 + tab->M;1135 1136	if (tab->M) {1137		if (isl_int_is_pos(tab->mat->row[row][2]))1138			return 1;1139		if (isl_int_is_neg(tab->mat->row[row][2]))1140			return 0;1141	}1142 1143	if (isl_int_is_neg(tab->mat->row[row][1]))1144		return 0;1145	for (i = 0; i < tab->n_param; ++i) {1146		/* Eliminated parameter */1147		if (tab->var[i].is_row)1148			continue;1149		col = tab->var[i].index;1150		if (isl_int_is_zero(tab->mat->row[row][off + col]))1151			continue;1152		if (!tab->var[i].is_nonneg)1153			return 0;1154		if (isl_int_is_neg(tab->mat->row[row][off + col]))1155			return 0;1156	}1157	for (i = 0; i < tab->n_div; ++i) {1158		if (tab->var[tab->n_var - tab->n_div + i].is_row)1159			continue;1160		col = tab->var[tab->n_var - tab->n_div + i].index;1161		if (isl_int_is_zero(tab->mat->row[row][off + col]))1162			continue;1163		if (!tab->var[tab->n_var - tab->n_div + i].is_nonneg)1164			return 0;1165		if (isl_int_is_neg(tab->mat->row[row][off + col]))1166			return 0;1167	}1168	return 1;1169}1170 1171/* Given a row r and two columns, return the column that would1172 * lead to the lexicographically smallest increment in the sample1173 * solution when leaving the basis in favor of the row.1174 * Pivoting with column c will increment the sample value by a non-negative1175 * constant times a_{V,c}/a_{r,c}, with a_{V,c} the elements of column c1176 * corresponding to the non-parametric variables.1177 * If variable v appears in a column c_v, then a_{v,c} = 1 iff c = c_v,1178 * with all other entries in this virtual row equal to zero.1179 * If variable v appears in a row, then a_{v,c} is the element in column c1180 * of that row.1181 *1182 * Let v be the first variable with a_{v,c1}/a_{r,c1} != a_{v,c2}/a_{r,c2}.1183 * Then if a_{v,c1}/a_{r,c1} < a_{v,c2}/a_{r,c2}, i.e.,1184 * a_{v,c2} a_{r,c1} - a_{v,c1} a_{r,c2} > 0, c1 results in the minimal1185 * increment.  Otherwise, it's c2.1186 */1187static int lexmin_col_pair(struct isl_tab *tab,1188	int row, int col1, int col2, isl_int tmp)1189{1190	int i;1191	isl_int *tr;1192 1193	tr = tab->mat->row[row] + 2 + tab->M;1194 1195	for (i = tab->n_param; i < tab->n_var - tab->n_div; ++i) {1196		int s1, s2;1197		isl_int *r;1198 1199		if (!tab->var[i].is_row) {1200			if (tab->var[i].index == col1)1201				return col2;1202			if (tab->var[i].index == col2)1203				return col1;1204			continue;1205		}1206 1207		if (tab->var[i].index == row)1208			continue;1209 1210		r = tab->mat->row[tab->var[i].index] + 2 + tab->M;1211		s1 = isl_int_sgn(r[col1]);1212		s2 = isl_int_sgn(r[col2]);1213		if (s1 == 0 && s2 == 0)1214			continue;1215		if (s1 < s2)1216			return col1;1217		if (s2 < s1)1218			return col2;1219 1220		isl_int_mul(tmp, r[col2], tr[col1]);1221		isl_int_submul(tmp, r[col1], tr[col2]);1222		if (isl_int_is_pos(tmp))1223			return col1;1224		if (isl_int_is_neg(tmp))1225			return col2;1226	}1227	return -1;1228}1229 1230/* Does the index into the tab->var or tab->con array "index"1231 * correspond to a variable in the context tableau?1232 * In particular, it needs to be an index into the tab->var array and1233 * it needs to refer to either one of the first tab->n_param variables or1234 * one of the last tab->n_div variables.1235 */1236static int is_parameter_var(struct isl_tab *tab, int index)1237{1238	if (index < 0)1239		return 0;1240	if (index < tab->n_param)1241		return 1;1242	if (index >= tab->n_var - tab->n_div)1243		return 1;1244	return 0;1245}1246 1247/* Does column "col" of "tab" refer to a variable in the context tableau?1248 */1249static int col_is_parameter_var(struct isl_tab *tab, int col)1250{1251	return is_parameter_var(tab, tab->col_var[col]);1252}1253 1254/* Does row "row" of "tab" refer to a variable in the context tableau?1255 */1256static int row_is_parameter_var(struct isl_tab *tab, int row)1257{1258	return is_parameter_var(tab, tab->row_var[row]);1259}1260 1261/* Given a row in the tableau, find and return the column that would1262 * result in the lexicographically smallest, but positive, increment1263 * in the sample point.1264 * If there is no such column, then return tab->n_col.1265 * If anything goes wrong, return -1.1266 */1267static int lexmin_pivot_col(struct isl_tab *tab, int row)1268{1269	int j;1270	int col = tab->n_col;1271	isl_int *tr;1272	isl_int tmp;1273 1274	tr = tab->mat->row[row] + 2 + tab->M;1275 1276	isl_int_init(tmp);1277 1278	for (j = tab->n_dead; j < tab->n_col; ++j) {1279		if (col_is_parameter_var(tab, j))1280			continue;1281 1282		if (!isl_int_is_pos(tr[j]))1283			continue;1284 1285		if (col == tab->n_col)1286			col = j;1287		else1288			col = lexmin_col_pair(tab, row, col, j, tmp);1289		isl_assert(tab->mat->ctx, col >= 0, goto error);1290	}1291 1292	isl_int_clear(tmp);1293	return col;1294error:1295	isl_int_clear(tmp);1296	return -1;1297}1298 1299/* Return the first known violated constraint, i.e., a non-negative1300 * constraint that currently has an either obviously negative value1301 * or a previously determined to be negative value.1302 *1303 * If any constraint has a negative coefficient for the big parameter,1304 * if any, then we return one of these first.1305 */1306static int first_neg(struct isl_tab *tab)1307{1308	int row;1309 1310	if (tab->M)1311		for (row = tab->n_redundant; row < tab->n_row; ++row) {1312			if (!isl_tab_var_from_row(tab, row)->is_nonneg)1313				continue;1314			if (!isl_int_is_neg(tab->mat->row[row][2]))1315				continue;1316			if (tab->row_sign)1317				tab->row_sign[row] = isl_tab_row_neg;1318			return row;1319		}1320	for (row = tab->n_redundant; row < tab->n_row; ++row) {1321		if (!isl_tab_var_from_row(tab, row)->is_nonneg)1322			continue;1323		if (tab->row_sign) {1324			if (tab->row_sign[row] == 0 &&1325			    is_obviously_neg(tab, row))1326				tab->row_sign[row] = isl_tab_row_neg;1327			if (tab->row_sign[row] != isl_tab_row_neg)1328				continue;1329		} else if (!is_obviously_neg(tab, row))1330			continue;1331		return row;1332	}1333	return -1;1334}1335 1336/* Check whether the invariant that all columns are lexico-positive1337 * is satisfied.  This function is not called from the current code1338 * but is useful during debugging.1339 */1340static void check_lexpos(struct isl_tab *tab) __attribute__ ((unused));1341static void check_lexpos(struct isl_tab *tab)1342{1343	unsigned off = 2 + tab->M;1344	int col;1345	int var;1346	int row;1347 1348	for (col = tab->n_dead; col < tab->n_col; ++col) {1349		if (col_is_parameter_var(tab, col))1350			continue;1351		for (var = tab->n_param; var < tab->n_var - tab->n_div; ++var) {1352			if (!tab->var[var].is_row) {1353				if (tab->var[var].index == col)1354					break;1355				else1356					continue;1357			}1358			row = tab->var[var].index;1359			if (isl_int_is_zero(tab->mat->row[row][off + col]))1360				continue;1361			if (isl_int_is_pos(tab->mat->row[row][off + col]))1362				break;1363			fprintf(stderr, "lexneg column %d (row %d)\n",1364				col, row);1365		}1366		if (var >= tab->n_var - tab->n_div)1367			fprintf(stderr, "zero column %d\n", col);1368	}1369}1370 1371/* Report to the caller that the given constraint is part of an encountered1372 * conflict.1373 */1374static int report_conflicting_constraint(struct isl_tab *tab, int con)1375{1376	return tab->conflict(con, tab->conflict_user);1377}1378 1379/* Given a conflicting row in the tableau, report all constraints1380 * involved in the row to the caller.  That is, the row itself1381 * (if it represents a constraint) and all constraint columns with1382 * non-zero (and therefore negative) coefficients.1383 */1384static int report_conflict(struct isl_tab *tab, int row)1385{1386	int j;1387	isl_int *tr;1388 1389	if (!tab->conflict)1390		return 0;1391 1392	if (tab->row_var[row] < 0 &&1393	    report_conflicting_constraint(tab, ~tab->row_var[row]) < 0)1394		return -1;1395 1396	tr = tab->mat->row[row] + 2 + tab->M;1397 1398	for (j = tab->n_dead; j < tab->n_col; ++j) {1399		if (col_is_parameter_var(tab, j))1400			continue;1401 1402		if (!isl_int_is_neg(tr[j]))1403			continue;1404 1405		if (tab->col_var[j] < 0 &&1406		    report_conflicting_constraint(tab, ~tab->col_var[j]) < 0)1407			return -1;1408	}1409 1410	return 0;1411}1412 1413/* Resolve all known or obviously violated constraints through pivoting.1414 * In particular, as long as we can find any violated constraint, we1415 * look for a pivoting column that would result in the lexicographically1416 * smallest increment in the sample point.  If there is no such column1417 * then the tableau is infeasible.1418 */1419static int restore_lexmin(struct isl_tab *tab) WARN_UNUSED;1420static int restore_lexmin(struct isl_tab *tab)1421{1422	int row, col;1423 1424	if (!tab)1425		return -1;1426	if (tab->empty)1427		return 0;1428	while ((row = first_neg(tab)) != -1) {1429		col = lexmin_pivot_col(tab, row);1430		if (col >= tab->n_col) {1431			if (report_conflict(tab, row) < 0)1432				return -1;1433			if (isl_tab_mark_empty(tab) < 0)1434				return -1;1435			return 0;1436		}1437		if (col < 0)1438			return -1;1439		if (isl_tab_pivot(tab, row, col) < 0)1440			return -1;1441	}1442	return 0;1443}1444 1445/* Given a row that represents an equality, look for an appropriate1446 * pivoting column.1447 * In particular, if there are any non-zero coefficients among1448 * the non-parameter variables, then we take the last of these1449 * variables.  Eliminating this variable in terms of the other1450 * variables and/or parameters does not influence the property1451 * that all column in the initial tableau are lexicographically1452 * positive.  The row corresponding to the eliminated variable1453 * will only have non-zero entries below the diagonal of the1454 * initial tableau.  That is, we transform1455 *1456 *		I				I1457 *		  1		into		a1458 *		    I				  I1459 *1460 * If there is no such non-parameter variable, then we are dealing with1461 * pure parameter equality and we pick any parameter with coefficient 1 or -11462 * for elimination.  This will ensure that the eliminated parameter1463 * always has an integer value whenever all the other parameters are integral.1464 * If there is no such parameter then we return -1.1465 */1466static int last_var_col_or_int_par_col(struct isl_tab *tab, int row)1467{1468	unsigned off = 2 + tab->M;1469	int i;1470 1471	for (i = tab->n_var - tab->n_div - 1; i >= 0 && i >= tab->n_param; --i) {1472		int col;1473		if (tab->var[i].is_row)1474			continue;1475		col = tab->var[i].index;1476		if (col <= tab->n_dead)1477			continue;1478		if (!isl_int_is_zero(tab->mat->row[row][off + col]))1479			return col;1480	}1481	for (i = tab->n_dead; i < tab->n_col; ++i) {1482		if (isl_int_is_one(tab->mat->row[row][off + i]))1483			return i;1484		if (isl_int_is_negone(tab->mat->row[row][off + i]))1485			return i;1486	}1487	return -1;1488}1489 1490/* Add an equality that is known to be valid to the tableau.1491 * We first check if we can eliminate a variable or a parameter.1492 * If not, we add the equality as two inequalities.1493 * In this case, the equality was a pure parameter equality and there1494 * is no need to resolve any constraint violations.1495 *1496 * This function assumes that at least two more rows and at least1497 * two more elements in the constraint array are available in the tableau.1498 */1499static struct isl_tab *add_lexmin_valid_eq(struct isl_tab *tab, isl_int *eq)1500{1501	int i;1502	int r;1503 1504	if (!tab)1505		return NULL;1506	r = isl_tab_add_row(tab, eq);1507	if (r < 0)1508		goto error;1509 1510	r = tab->con[r].index;1511	i = last_var_col_or_int_par_col(tab, r);1512	if (i < 0) {1513		tab->con[r].is_nonneg = 1;1514		if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)1515			goto error;1516		isl_seq_neg(eq, eq, 1 + tab->n_var);1517		r = isl_tab_add_row(tab, eq);1518		if (r < 0)1519			goto error;1520		tab->con[r].is_nonneg = 1;1521		if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)1522			goto error;1523	} else {1524		if (isl_tab_pivot(tab, r, i) < 0)1525			goto error;1526		if (isl_tab_kill_col(tab, i) < 0)1527			goto error;1528		tab->n_eq++;1529	}1530 1531	return tab;1532error:1533	isl_tab_free(tab);1534	return NULL;1535}1536 1537/* Check if the given row is a pure constant.1538 */1539static int is_constant(struct isl_tab *tab, int row)1540{1541	unsigned off = 2 + tab->M;1542 1543	return isl_seq_first_non_zero(tab->mat->row[row] + off + tab->n_dead,1544					tab->n_col - tab->n_dead) == -1;1545}1546 1547/* Is the given row a parametric constant?1548 * That is, does it only involve variables that also appear in the context?1549 */1550static int is_parametric_constant(struct isl_tab *tab, int row)1551{1552	unsigned off = 2 + tab->M;1553	int col;1554 1555	for (col = tab->n_dead; col < tab->n_col; ++col) {1556		if (col_is_parameter_var(tab, col))1557			continue;1558		if (isl_int_is_zero(tab->mat->row[row][off + col]))1559			continue;1560		return 0;1561	}1562 1563	return 1;1564}1565 1566/* Add an equality that may or may not be valid to the tableau.1567 * If the resulting row is a pure constant, then it must be zero.1568 * Otherwise, the resulting tableau is empty.1569 *1570 * If the row is not a pure constant, then we add two inequalities,1571 * each time checking that they can be satisfied.1572 * In the end we try to use one of the two constraints to eliminate1573 * a column.1574 *1575 * This function assumes that at least two more rows and at least1576 * two more elements in the constraint array are available in the tableau.1577 */1578static int add_lexmin_eq(struct isl_tab *tab, isl_int *eq) WARN_UNUSED;1579static int add_lexmin_eq(struct isl_tab *tab, isl_int *eq)1580{1581	int r1, r2;1582	int row;1583	struct isl_tab_undo *snap;1584 1585	if (!tab)1586		return -1;1587	snap = isl_tab_snap(tab);1588	r1 = isl_tab_add_row(tab, eq);1589	if (r1 < 0)1590		return -1;1591	tab->con[r1].is_nonneg = 1;1592	if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r1]) < 0)1593		return -1;1594 1595	row = tab->con[r1].index;1596	if (is_constant(tab, row)) {1597		if (!isl_int_is_zero(tab->mat->row[row][1]) ||1598		    (tab->M && !isl_int_is_zero(tab->mat->row[row][2]))) {1599			if (isl_tab_mark_empty(tab) < 0)1600				return -1;1601			return 0;1602		}1603		if (isl_tab_rollback(tab, snap) < 0)1604			return -1;1605		return 0;1606	}1607 1608	if (restore_lexmin(tab) < 0)1609		return -1;1610	if (tab->empty)1611		return 0;1612 1613	isl_seq_neg(eq, eq, 1 + tab->n_var);1614 1615	r2 = isl_tab_add_row(tab, eq);1616	if (r2 < 0)1617		return -1;1618	tab->con[r2].is_nonneg = 1;1619	if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r2]) < 0)1620		return -1;1621 1622	if (restore_lexmin(tab) < 0)1623		return -1;1624	if (tab->empty)1625		return 0;1626 1627	if (!tab->con[r1].is_row) {1628		if (isl_tab_kill_col(tab, tab->con[r1].index) < 0)1629			return -1;1630	} else if (!tab->con[r2].is_row) {1631		if (isl_tab_kill_col(tab, tab->con[r2].index) < 0)1632			return -1;1633	}1634 1635	if (tab->bmap) {1636		tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);1637		if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)1638			return -1;1639		isl_seq_neg(eq, eq, 1 + tab->n_var);1640		tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq);1641		isl_seq_neg(eq, eq, 1 + tab->n_var);1642		if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)1643			return -1;1644		if (!tab->bmap)1645			return -1;1646	}1647 1648	return 0;1649}1650 1651/* Add an inequality to the tableau, resolving violations using1652 * restore_lexmin.1653 *1654 * This function assumes that at least one more row and at least1655 * one more element in the constraint array are available in the tableau.1656 */1657static struct isl_tab *add_lexmin_ineq(struct isl_tab *tab, isl_int *ineq)1658{1659	int r;1660 1661	if (!tab)1662		return NULL;1663	if (tab->bmap) {1664		tab->bmap = isl_basic_map_add_ineq(tab->bmap, ineq);1665		if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0)1666			goto error;1667		if (!tab->bmap)1668			goto error;1669	}1670	r = isl_tab_add_row(tab, ineq);1671	if (r < 0)1672		goto error;1673	tab->con[r].is_nonneg = 1;1674	if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)1675		goto error;1676	if (isl_tab_row_is_redundant(tab, tab->con[r].index)) {1677		if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0)1678			goto error;1679		return tab;1680	}1681 1682	if (restore_lexmin(tab) < 0)1683		goto error;1684	if (!tab->empty && tab->con[r].is_row &&1685		 isl_tab_row_is_redundant(tab, tab->con[r].index))1686		if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0)1687			goto error;1688	return tab;1689error:1690	isl_tab_free(tab);1691	return NULL;1692}1693 1694/* Check if the coefficients of the parameters are all integral.1695 */1696static int integer_parameter(struct isl_tab *tab, int row)1697{1698	int i;1699	int col;1700	unsigned off = 2 + tab->M;1701 1702	for (i = 0; i < tab->n_param; ++i) {1703		/* Eliminated parameter */1704		if (tab->var[i].is_row)1705			continue;1706		col = tab->var[i].index;1707		if (!isl_int_is_divisible_by(tab->mat->row[row][off + col],1708						tab->mat->row[row][0]))1709			return 0;1710	}1711	for (i = 0; i < tab->n_div; ++i) {1712		if (tab->var[tab->n_var - tab->n_div + i].is_row)1713			continue;1714		col = tab->var[tab->n_var - tab->n_div + i].index;1715		if (!isl_int_is_divisible_by(tab->mat->row[row][off + col],1716						tab->mat->row[row][0]))1717			return 0;1718	}1719	return 1;1720}1721 1722/* Check if the coefficients of the non-parameter variables are all integral.1723 */1724static int integer_variable(struct isl_tab *tab, int row)1725{1726	int i;1727	unsigned off = 2 + tab->M;1728 1729	for (i = tab->n_dead; i < tab->n_col; ++i) {1730		if (col_is_parameter_var(tab, i))1731			continue;1732		if (!isl_int_is_divisible_by(tab->mat->row[row][off + i],1733						tab->mat->row[row][0]))1734			return 0;1735	}1736	return 1;1737}1738 1739/* Check if the constant term is integral.1740 */1741static int integer_constant(struct isl_tab *tab, int row)1742{1743	return isl_int_is_divisible_by(tab->mat->row[row][1],1744					tab->mat->row[row][0]);1745}1746 1747#define I_CST	1 << 01748#define I_PAR	1 << 11749#define I_VAR	1 << 21750 1751/* Check for next (non-parameter) variable after "var" (first if var == -1)1752 * that is non-integer and therefore requires a cut and return1753 * the index of the variable.1754 * For parametric tableaus, there are three parts in a row,1755 * the constant, the coefficients of the parameters and the rest.1756 * For each part, we check whether the coefficients in that part1757 * are all integral and if so, set the corresponding flag in *f.1758 * If the constant and the parameter part are integral, then the1759 * current sample value is integral and no cut is required1760 * (irrespective of whether the variable part is integral).1761 */1762static int next_non_integer_var(struct isl_tab *tab, int var, int *f)1763{1764	var = var < 0 ? tab->n_param : var + 1;1765 1766	for (; var < tab->n_var - tab->n_div; ++var) {1767		int flags = 0;1768		int row;1769		if (!tab->var[var].is_row)1770			continue;1771		row = tab->var[var].index;1772		if (integer_constant(tab, row))1773			ISL_FL_SET(flags, I_CST);1774		if (integer_parameter(tab, row))1775			ISL_FL_SET(flags, I_PAR);1776		if (ISL_FL_ISSET(flags, I_CST) && ISL_FL_ISSET(flags, I_PAR))1777			continue;1778		if (integer_variable(tab, row))1779			ISL_FL_SET(flags, I_VAR);1780		*f = flags;1781		return var;1782	}1783	return -1;1784}1785 1786/* Check for first (non-parameter) variable that is non-integer and1787 * therefore requires a cut and return the corresponding row.1788 * For parametric tableaus, there are three parts in a row,1789 * the constant, the coefficients of the parameters and the rest.1790 * For each part, we check whether the coefficients in that part1791 * are all integral and if so, set the corresponding flag in *f.1792 * If the constant and the parameter part are integral, then the1793 * current sample value is integral and no cut is required1794 * (irrespective of whether the variable part is integral).1795 */1796static int first_non_integer_row(struct isl_tab *tab, int *f)1797{1798	int var = next_non_integer_var(tab, -1, f);1799 1800	return var < 0 ? -1 : tab->var[var].index;1801}1802 1803/* Add a (non-parametric) cut to cut away the non-integral sample1804 * value of the given row.1805 *1806 * If the row is given by1807 *1808 *	m r = f + \sum_i a_i y_i1809 *1810 * then the cut is1811 *1812 *	c = - {-f/m} + \sum_i {a_i/m} y_i >= 01813 *1814 * The big parameter, if any, is ignored, since it is assumed to be big1815 * enough to be divisible by any integer.1816 * If the tableau is actually a parametric tableau, then this function1817 * is only called when all coefficients of the parameters are integral.1818 * The cut therefore has zero coefficients for the parameters.1819 *1820 * The current value is known to be negative, so row_sign, if it1821 * exists, is set accordingly.1822 *1823 * Return the row of the cut or -1.1824 */1825static int add_cut(struct isl_tab *tab, int row)1826{1827	int i;1828	int r;1829	isl_int *r_row;1830	unsigned off = 2 + tab->M;1831 1832	if (isl_tab_extend_cons(tab, 1) < 0)1833		return -1;1834	r = isl_tab_allocate_con(tab);1835	if (r < 0)1836		return -1;1837 1838	r_row = tab->mat->row[tab->con[r].index];1839	isl_int_set(r_row[0], tab->mat->row[row][0]);1840	isl_int_neg(r_row[1], tab->mat->row[row][1]);1841	isl_int_fdiv_r(r_row[1], r_row[1], tab->mat->row[row][0]);1842	isl_int_neg(r_row[1], r_row[1]);1843	if (tab->M)1844		isl_int_set_si(r_row[2], 0);1845	for (i = 0; i < tab->n_col; ++i)1846		isl_int_fdiv_r(r_row[off + i],1847			tab->mat->row[row][off + i], tab->mat->row[row][0]);1848 1849	tab->con[r].is_nonneg = 1;1850	if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)1851		return -1;1852	if (tab->row_sign)1853		tab->row_sign[tab->con[r].index] = isl_tab_row_neg;1854 1855	return tab->con[r].index;1856}1857 1858#define CUT_ALL 11859#define CUT_ONE 01860 1861/* Given a non-parametric tableau, add cuts until an integer1862 * sample point is obtained or until the tableau is determined1863 * to be integer infeasible.1864 * As long as there is any non-integer value in the sample point,1865 * we add appropriate cuts, if possible, for each of these1866 * non-integer values and then resolve the violated1867 * cut constraints using restore_lexmin.1868 * If one of the corresponding rows is equal to an integral1869 * combination of variables/constraints plus a non-integral constant,1870 * then there is no way to obtain an integer point and we return1871 * a tableau that is marked empty.1872 * The parameter cutting_strategy controls the strategy used when adding cuts1873 * to remove non-integer points. CUT_ALL adds all possible cuts1874 * before continuing the search. CUT_ONE adds only one cut at a time.1875 */1876static struct isl_tab *cut_to_integer_lexmin(struct isl_tab *tab,1877	int cutting_strategy)1878{1879	int var;1880	int row;1881	int flags;1882 1883	if (!tab)1884		return NULL;1885	if (tab->empty)1886		return tab;1887 1888	while ((var = next_non_integer_var(tab, -1, &flags)) != -1) {1889		do {1890			if (ISL_FL_ISSET(flags, I_VAR)) {1891				if (isl_tab_mark_empty(tab) < 0)1892					goto error;1893				return tab;1894			}1895			row = tab->var[var].index;1896			row = add_cut(tab, row);1897			if (row < 0)1898				goto error;1899			if (cutting_strategy == CUT_ONE)1900				break;1901		} while ((var = next_non_integer_var(tab, var, &flags)) != -1);1902		if (restore_lexmin(tab) < 0)1903			goto error;1904		if (tab->empty)1905			break;1906	}1907	return tab;1908error:1909	isl_tab_free(tab);1910	return NULL;1911}1912 1913/* Check whether all the currently active samples also satisfy the inequality1914 * "ineq" (treated as an equality if eq is set).1915 * Remove those samples that do not.1916 */1917static struct isl_tab *check_samples(struct isl_tab *tab, isl_int *ineq, int eq)1918{1919	int i;1920	isl_int v;1921 1922	if (!tab)1923		return NULL;1924 1925	isl_assert(tab->mat->ctx, tab->bmap, goto error);1926	isl_assert(tab->mat->ctx, tab->samples, goto error);1927	isl_assert(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var, goto error);1928 1929	isl_int_init(v);1930	for (i = tab->n_outside; i < tab->n_sample; ++i) {1931		int sgn;1932		isl_seq_inner_product(ineq, tab->samples->row[i],1933					1 + tab->n_var, &v);1934		sgn = isl_int_sgn(v);1935		if (eq ? (sgn == 0) : (sgn >= 0))1936			continue;1937		tab = isl_tab_drop_sample(tab, i);1938		if (!tab)1939			break;1940	}1941	isl_int_clear(v);1942 1943	return tab;1944error:1945	isl_tab_free(tab);1946	return NULL;1947}1948 1949/* Check whether the sample value of the tableau is finite,1950 * i.e., either the tableau does not use a big parameter, or1951 * all values of the variables are equal to the big parameter plus1952 * some constant.  This constant is the actual sample value.1953 */1954static int sample_is_finite(struct isl_tab *tab)1955{1956	int i;1957 1958	if (!tab->M)1959		return 1;1960 1961	for (i = 0; i < tab->n_var; ++i) {1962		int row;1963		if (!tab->var[i].is_row)1964			return 0;1965		row = tab->var[i].index;1966		if (isl_int_ne(tab->mat->row[row][0], tab->mat->row[row][2]))1967			return 0;1968	}1969	return 1;1970}1971 1972/* Check if the context tableau of sol has any integer points.1973 * Leave tab in empty state if no integer point can be found.1974 * If an integer point can be found and if moreover it is finite,1975 * then it is added to the list of sample values.1976 *1977 * This function is only called when none of the currently active sample1978 * values satisfies the most recently added constraint.1979 */1980static struct isl_tab *check_integer_feasible(struct isl_tab *tab)1981{1982	struct isl_tab_undo *snap;1983 1984	if (!tab)1985		return NULL;1986 1987	snap = isl_tab_snap(tab);1988	if (isl_tab_push_basis(tab) < 0)1989		goto error;1990 1991	tab = cut_to_integer_lexmin(tab, CUT_ALL);1992	if (!tab)1993		goto error;1994 1995	if (!tab->empty && sample_is_finite(tab)) {1996		struct isl_vec *sample;1997 1998		sample = isl_tab_get_sample_value(tab);1999 2000		if (isl_tab_add_sample(tab, sample) < 0)2001			goto error;2002	}2003 2004	if (!tab->empty && isl_tab_rollback(tab, snap) < 0)2005		goto error;2006 2007	return tab;2008error:2009	isl_tab_free(tab);2010	return NULL;2011}2012 2013/* Check if any of the currently active sample values satisfies2014 * the inequality "ineq" (an equality if eq is set).2015 */2016static int tab_has_valid_sample(struct isl_tab *tab, isl_int *ineq, int eq)2017{2018	int i;2019	isl_int v;2020 2021	if (!tab)2022		return -1;2023 2024	isl_assert(tab->mat->ctx, tab->bmap, return -1);2025	isl_assert(tab->mat->ctx, tab->samples, return -1);2026	isl_assert(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var, return -1);2027 2028	isl_int_init(v);2029	for (i = tab->n_outside; i < tab->n_sample; ++i) {2030		int sgn;2031		isl_seq_inner_product(ineq, tab->samples->row[i],2032					1 + tab->n_var, &v);2033		sgn = isl_int_sgn(v);2034		if (eq ? (sgn == 0) : (sgn >= 0))2035			break;2036	}2037	isl_int_clear(v);2038 2039	return i < tab->n_sample;2040}2041 2042/* Insert a div specified by "div" to the tableau "tab" at position "pos" and2043 * return isl_bool_true if the div is obviously non-negative.2044 */2045static isl_bool context_tab_insert_div(struct isl_tab *tab, int pos,2046	__isl_keep isl_vec *div,2047	isl_stat (*add_ineq)(void *user, isl_int *), void *user)2048{2049	int i;2050	int r;2051	struct isl_mat *samples;2052	int nonneg;2053 2054	r = isl_tab_insert_div(tab, pos, div, add_ineq, user);2055	if (r < 0)2056		return isl_bool_error;2057	nonneg = tab->var[r].is_nonneg;2058	tab->var[r].frozen = 1;2059 2060	samples = isl_mat_extend(tab->samples,2061			tab->n_sample, 1 + tab->n_var);2062	tab->samples = samples;2063	if (!samples)2064		return isl_bool_error;2065	for (i = tab->n_outside; i < samples->n_row; ++i) {2066		isl_seq_inner_product(div->el + 1, samples->row[i],2067			div->size - 1, &samples->row[i][samples->n_col - 1]);2068		isl_int_fdiv_q(samples->row[i][samples->n_col - 1],2069			       samples->row[i][samples->n_col - 1], div->el[0]);2070	}2071	tab->samples = isl_mat_move_cols(tab->samples, 1 + pos,2072					1 + tab->n_var - 1, 1);2073	if (!tab->samples)2074		return isl_bool_error;2075 2076	return isl_bool_ok(nonneg);2077}2078 2079/* Add a div specified by "div" to both the main tableau and2080 * the context tableau.  In case of the main tableau, we only2081 * need to add an extra div.  In the context tableau, we also2082 * need to express the meaning of the div.2083 * Return the index of the div or -1 if anything went wrong.2084 *2085 * The new integer division is added before any unknown integer2086 * divisions in the context to ensure that it does not get2087 * equated to some linear combination involving unknown integer2088 * divisions.2089 */2090static int add_div(struct isl_tab *tab, struct isl_context *context,2091	__isl_keep isl_vec *div)2092{2093	int r;2094	int pos;2095	isl_bool nonneg;2096	struct isl_tab *context_tab = context->op->peek_tab(context);2097 2098	if (!tab || !context_tab)2099		goto error;2100 2101	pos = context_tab->n_var - context->n_unknown;2102	if ((nonneg = context->op->insert_div(context, pos, div)) < 0)2103		goto error;2104 2105	if (!context->op->is_ok(context))2106		goto error;2107 2108	pos = tab->n_var - context->n_unknown;2109	if (isl_tab_extend_vars(tab, 1) < 0)2110		goto error;2111	r = isl_tab_insert_var(tab, pos);2112	if (r < 0)2113		goto error;2114	if (nonneg)2115		tab->var[r].is_nonneg = 1;2116	tab->var[r].frozen = 1;2117	tab->n_div++;2118 2119	return tab->n_div - 1 - context->n_unknown;2120error:2121	context->op->invalidate(context);2122	return -1;2123}2124 2125/* Return the position of the integer division that is equal to div/denom2126 * if there is one.  Otherwise, return a position beyond the integer divisions.2127 */2128static int find_div(struct isl_tab *tab, isl_int *div, isl_int denom)2129{2130	int i;2131	isl_size total = isl_basic_map_dim(tab->bmap, isl_dim_all);2132	isl_size n_div;2133 2134	n_div = isl_basic_map_dim(tab->bmap, isl_dim_div);2135	if (total < 0 || n_div < 0)2136		return -1;2137	for (i = 0; i < n_div; ++i) {2138		if (isl_int_ne(tab->bmap->div[i][0], denom))2139			continue;2140		if (!isl_seq_eq(tab->bmap->div[i] + 1, div, 1 + total))2141			continue;2142		return i;2143	}2144	return n_div;2145}2146 2147/* Return the index of a div that corresponds to "div".2148 * We first check if we already have such a div and if not, we create one.2149 */2150static int get_div(struct isl_tab *tab, struct isl_context *context,2151	struct isl_vec *div)2152{2153	int d;2154	struct isl_tab *context_tab = context->op->peek_tab(context);2155	unsigned n_div;2156 2157	if (!context_tab)2158		return -1;2159 2160	n_div = isl_basic_map_dim(context_tab->bmap, isl_dim_div);2161	d = find_div(context_tab, div->el + 1, div->el[0]);2162	if (d < 0)2163		return -1;2164	if (d < n_div)2165		return d;2166 2167	return add_div(tab, context, div);2168}2169 2170/* Add a parametric cut to cut away the non-integral sample value2171 * of the given row.2172 * Let a_i be the coefficients of the constant term and the parameters2173 * and let b_i be the coefficients of the variables or constraints2174 * in basis of the tableau.2175 * Let q be the div q = floor(\sum_i {-a_i} y_i).2176 *2177 * The cut is expressed as2178 *2179 *	c = \sum_i -{-a_i} y_i + \sum_i {b_i} x_i + q >= 02180 *2181 * If q did not already exist in the context tableau, then it is added first.2182 * If q is in a column of the main tableau then the "+ q" can be accomplished2183 * by setting the corresponding entry to the denominator of the constraint.2184 * If q happens to be in a row of the main tableau, then the corresponding2185 * row needs to be added instead (taking care of the denominators).2186 * Note that this is very unlikely, but perhaps not entirely impossible.2187 *2188 * The current value of the cut is known to be negative (or at least2189 * non-positive), so row_sign is set accordingly.2190 *2191 * Return the row of the cut or -1.2192 */2193static int add_parametric_cut(struct isl_tab *tab, int row,2194	struct isl_context *context)2195{2196	struct isl_vec *div;2197	int d;2198	int i;2199	int r;2200	isl_int *r_row;2201	int col;2202	int n;2203	unsigned off = 2 + tab->M;2204 2205	if (!context)2206		return -1;2207 2208	div = get_row_parameter_div(tab, row);2209	if (!div)2210		return -1;2211 2212	n = tab->n_div - context->n_unknown;2213	d = context->op->get_div(context, tab, div);2214	isl_vec_free(div);2215	if (d < 0)2216		return -1;2217 2218	if (isl_tab_extend_cons(tab, 1) < 0)2219		return -1;2220	r = isl_tab_allocate_con(tab);2221	if (r < 0)2222		return -1;2223 2224	r_row = tab->mat->row[tab->con[r].index];2225	isl_int_set(r_row[0], tab->mat->row[row][0]);2226	isl_int_neg(r_row[1], tab->mat->row[row][1]);2227	isl_int_fdiv_r(r_row[1], r_row[1], tab->mat->row[row][0]);2228	isl_int_neg(r_row[1], r_row[1]);2229	if (tab->M)2230		isl_int_set_si(r_row[2], 0);2231	for (i = 0; i < tab->n_param; ++i) {2232		if (tab->var[i].is_row)2233			continue;2234		col = tab->var[i].index;2235		isl_int_neg(r_row[off + col], tab->mat->row[row][off + col]);2236		isl_int_fdiv_r(r_row[off + col], r_row[off + col],2237				tab->mat->row[row][0]);2238		isl_int_neg(r_row[off + col], r_row[off + col]);2239	}2240	for (i = 0; i < tab->n_div; ++i) {2241		if (tab->var[tab->n_var - tab->n_div + i].is_row)2242			continue;2243		col = tab->var[tab->n_var - tab->n_div + i].index;2244		isl_int_neg(r_row[off + col], tab->mat->row[row][off + col]);2245		isl_int_fdiv_r(r_row[off + col], r_row[off + col],2246				tab->mat->row[row][0]);2247		isl_int_neg(r_row[off + col], r_row[off + col]);2248	}2249	for (i = 0; i < tab->n_col; ++i) {2250		if (tab->col_var[i] >= 0 &&2251		    (tab->col_var[i] < tab->n_param ||2252		     tab->col_var[i] >= tab->n_var - tab->n_div))2253			continue;2254		isl_int_fdiv_r(r_row[off + i],2255			tab->mat->row[row][off + i], tab->mat->row[row][0]);2256	}2257	if (tab->var[tab->n_var - tab->n_div + d].is_row) {2258		isl_int gcd;2259		int d_row = tab->var[tab->n_var - tab->n_div + d].index;2260		isl_int_init(gcd);2261		isl_int_gcd(gcd, tab->mat->row[d_row][0], r_row[0]);2262		isl_int_divexact(r_row[0], r_row[0], gcd);2263		isl_int_divexact(gcd, tab->mat->row[d_row][0], gcd);2264		isl_seq_combine(r_row + 1, gcd, r_row + 1,2265				r_row[0], tab->mat->row[d_row] + 1,2266				off - 1 + tab->n_col);2267		isl_int_mul(r_row[0], r_row[0], tab->mat->row[d_row][0]);2268		isl_int_clear(gcd);2269	} else {2270		col = tab->var[tab->n_var - tab->n_div + d].index;2271		isl_int_set(r_row[off + col], tab->mat->row[row][0]);2272	}2273 2274	tab->con[r].is_nonneg = 1;2275	if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0)2276		return -1;2277	if (tab->row_sign)2278		tab->row_sign[tab->con[r].index] = isl_tab_row_neg;2279 2280	row = tab->con[r].index;2281 2282	if (d >= n && context->op->detect_equalities(context, tab) < 0)2283		return -1;2284 2285	return row;2286}2287 2288/* Construct a tableau for bmap that can be used for computing2289 * the lexicographic minimum (or maximum) of bmap.2290 * If not NULL, then dom is the domain where the minimum2291 * should be computed.  In this case, we set up a parametric2292 * tableau with row signs (initialized to "unknown").2293 * If M is set, then the tableau will use a big parameter.2294 * If max is set, then a maximum should be computed instead of a minimum.2295 * This means that for each variable x, the tableau will contain the variable2296 * x' = M - x, rather than x' = M + x.  This in turn means that the coefficient2297 * of the variables in all constraints are negated prior to adding them2298 * to the tableau.2299 */2300static __isl_give struct isl_tab *tab_for_lexmin(__isl_keep isl_basic_map *bmap,2301	__isl_keep isl_basic_set *dom, unsigned M, int max)2302{2303	int i;2304	struct isl_tab *tab;2305	unsigned n_var;2306	unsigned o_var;2307	isl_size total;2308 2309	total = isl_basic_map_dim(bmap, isl_dim_all);2310	if (total < 0)2311		return NULL;2312	tab = isl_tab_alloc(bmap->ctx, 2 * bmap->n_eq + bmap->n_ineq + 1,2313			    total, M);2314	if (!tab)2315		return NULL;2316 2317	tab->rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);2318	if (dom) {2319		isl_size dom_total;2320		dom_total = isl_basic_set_dim(dom, isl_dim_all);2321		if (dom_total < 0)2322			goto error;2323		tab->n_param = dom_total - dom->n_div;2324		tab->n_div = dom->n_div;2325		tab->row_sign = isl_calloc_array(bmap->ctx,2326					enum isl_tab_row_sign, tab->mat->n_row);2327		if (tab->mat->n_row && !tab->row_sign)2328			goto error;2329	}2330	if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY)) {2331		if (isl_tab_mark_empty(tab) < 0)2332			goto error;2333		return tab;2334	}2335 2336	for (i = tab->n_param; i < tab->n_var - tab->n_div; ++i) {2337		tab->var[i].is_nonneg = 1;2338		tab->var[i].frozen = 1;2339	}2340	o_var = 1 + tab->n_param;2341	n_var = tab->n_var - tab->n_param - tab->n_div;2342	for (i = 0; i < bmap->n_eq; ++i) {2343		if (max)2344			isl_seq_neg(bmap->eq[i] + o_var,2345				    bmap->eq[i] + o_var, n_var);2346		tab = add_lexmin_valid_eq(tab, bmap->eq[i]);2347		if (max)2348			isl_seq_neg(bmap->eq[i] + o_var,2349				    bmap->eq[i] + o_var, n_var);2350		if (!tab || tab->empty)2351			return tab;2352	}2353	if (bmap->n_eq && restore_lexmin(tab) < 0)2354		goto error;2355	for (i = 0; i < bmap->n_ineq; ++i) {2356		if (max)2357			isl_seq_neg(bmap->ineq[i] + o_var,2358				    bmap->ineq[i] + o_var, n_var);2359		tab = add_lexmin_ineq(tab, bmap->ineq[i]);2360		if (max)2361			isl_seq_neg(bmap->ineq[i] + o_var,2362				    bmap->ineq[i] + o_var, n_var);2363		if (!tab || tab->empty)2364			return tab;2365	}2366	return tab;2367error:2368	isl_tab_free(tab);2369	return NULL;2370}2371 2372/* Given a main tableau where more than one row requires a split,2373 * determine and return the "best" row to split on.2374 *2375 * If any of the rows requiring a split only involves2376 * variables that also appear in the context tableau,2377 * then the negative part is guaranteed not to have a solution.2378 * It is therefore best to split on any of these rows first.2379 *2380 * Otherwise,2381 * given two rows in the main tableau, if the inequality corresponding2382 * to the first row is redundant with respect to that of the second row2383 * in the current tableau, then it is better to split on the second row,2384 * since in the positive part, both rows will be positive.2385 * (In the negative part a pivot will have to be performed and just about2386 * anything can happen to the sign of the other row.)2387 *2388 * As a simple heuristic, we therefore select the row that makes the most2389 * of the other rows redundant.2390 *2391 * Perhaps it would also be useful to look at the number of constraints2392 * that conflict with any given constraint.2393 *2394 * best is the best row so far (-1 when we have not found any row yet).2395 * best_r is the number of other rows made redundant by row best.2396 * When best is still -1, bset_r is meaningless, but it is initialized2397 * to some arbitrary value (0) anyway.  Without this redundant initialization2398 * valgrind may warn about uninitialized memory accesses when isl2399 * is compiled with some versions of gcc.2400 */2401static int best_split(struct isl_tab *tab, struct isl_tab *context_tab)2402{2403	struct isl_tab_undo *snap;2404	int split;2405	int row;2406	int best = -1;2407	int best_r = 0;2408 2409	if (isl_tab_extend_cons(context_tab, 2) < 0)2410		return -1;2411 2412	snap = isl_tab_snap(context_tab);2413 2414	for (split = tab->n_redundant; split < tab->n_row; ++split) {2415		struct isl_tab_undo *snap2;2416		struct isl_vec *ineq = NULL;2417		int r = 0;2418		int ok;2419 2420		if (!isl_tab_var_from_row(tab, split)->is_nonneg)2421			continue;2422		if (tab->row_sign[split] != isl_tab_row_any)2423			continue;2424 2425		if (is_parametric_constant(tab, split))2426			return split;2427 2428		ineq = get_row_parameter_ineq(tab, split);2429		if (!ineq)2430			return -1;2431		ok = isl_tab_add_ineq(context_tab, ineq->el) >= 0;2432		isl_vec_free(ineq);2433		if (!ok)2434			return -1;2435 2436		snap2 = isl_tab_snap(context_tab);2437 2438		for (row = tab->n_redundant; row < tab->n_row; ++row) {2439			struct isl_tab_var *var;2440 2441			if (row == split)2442				continue;2443			if (!isl_tab_var_from_row(tab, row)->is_nonneg)2444				continue;2445			if (tab->row_sign[row] != isl_tab_row_any)2446				continue;2447 2448			ineq = get_row_parameter_ineq(tab, row);2449			if (!ineq)2450				return -1;2451			ok = isl_tab_add_ineq(context_tab, ineq->el) >= 0;2452			isl_vec_free(ineq);2453			if (!ok)2454				return -1;2455			var = &context_tab->con[context_tab->n_con - 1];2456			if (!context_tab->empty &&2457			    !isl_tab_min_at_most_neg_one(context_tab, var))2458				r++;2459			if (isl_tab_rollback(context_tab, snap2) < 0)2460				return -1;2461		}2462		if (best == -1 || r > best_r) {2463			best = split;2464			best_r = r;2465		}2466		if (isl_tab_rollback(context_tab, snap) < 0)2467			return -1;2468	}2469 2470	return best;2471}2472 2473static struct isl_basic_set *context_lex_peek_basic_set(2474	struct isl_context *context)2475{2476	struct isl_context_lex *clex = (struct isl_context_lex *)context;2477	if (!clex->tab)2478		return NULL;2479	return isl_tab_peek_bset(clex->tab);2480}2481 2482static struct isl_tab *context_lex_peek_tab(struct isl_context *context)2483{2484	struct isl_context_lex *clex = (struct isl_context_lex *)context;2485	return clex->tab;2486}2487 2488static void context_lex_add_eq(struct isl_context *context, isl_int *eq,2489		int check, int update)2490{2491	struct isl_context_lex *clex = (struct isl_context_lex *)context;2492	if (isl_tab_extend_cons(clex->tab, 2) < 0)2493		goto error;2494	if (add_lexmin_eq(clex->tab, eq) < 0)2495		goto error;2496	if (check) {2497		int v = tab_has_valid_sample(clex->tab, eq, 1);2498		if (v < 0)2499			goto error;2500		if (!v)2501			clex->tab = check_integer_feasible(clex->tab);2502	}2503	if (update)2504		clex->tab = check_samples(clex->tab, eq, 1);2505	return;2506error:2507	isl_tab_free(clex->tab);2508	clex->tab = NULL;2509}2510 2511static void context_lex_add_ineq(struct isl_context *context, isl_int *ineq,2512		int check, int update)2513{2514	struct isl_context_lex *clex = (struct isl_context_lex *)context;2515	if (isl_tab_extend_cons(clex->tab, 1) < 0)2516		goto error;2517	clex->tab = add_lexmin_ineq(clex->tab, ineq);2518	if (check) {2519		int v = tab_has_valid_sample(clex->tab, ineq, 0);2520		if (v < 0)2521			goto error;2522		if (!v)2523			clex->tab = check_integer_feasible(clex->tab);2524	}2525	if (update)2526		clex->tab = check_samples(clex->tab, ineq, 0);2527	return;2528error:2529	isl_tab_free(clex->tab);2530	clex->tab = NULL;2531}2532 2533static isl_stat context_lex_add_ineq_wrap(void *user, isl_int *ineq)2534{2535	struct isl_context *context = (struct isl_context *)user;2536	context_lex_add_ineq(context, ineq, 0, 0);2537	return context->op->is_ok(context) ? isl_stat_ok : isl_stat_error;2538}2539 2540/* Check which signs can be obtained by "ineq" on all the currently2541 * active sample values.  See row_sign for more information.2542 */2543static enum isl_tab_row_sign tab_ineq_sign(struct isl_tab *tab, isl_int *ineq,2544	int strict)2545{2546	int i;2547	int sgn;2548	isl_int tmp;2549	enum isl_tab_row_sign res = isl_tab_row_unknown;2550 2551	isl_assert(tab->mat->ctx, tab->samples, return isl_tab_row_unknown);2552	isl_assert(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var,2553			return isl_tab_row_unknown);2554 2555	isl_int_init(tmp);2556	for (i = tab->n_outside; i < tab->n_sample; ++i) {2557		isl_seq_inner_product(tab->samples->row[i], ineq,2558					1 + tab->n_var, &tmp);2559		sgn = isl_int_sgn(tmp);2560		if (sgn > 0 || (sgn == 0 && strict)) {2561			if (res == isl_tab_row_unknown)2562				res = isl_tab_row_pos;2563			if (res == isl_tab_row_neg)2564				res = isl_tab_row_any;2565		}2566		if (sgn < 0) {2567			if (res == isl_tab_row_unknown)2568				res = isl_tab_row_neg;2569			if (res == isl_tab_row_pos)2570				res = isl_tab_row_any;2571		}2572		if (res == isl_tab_row_any)2573			break;2574	}2575	isl_int_clear(tmp);2576 2577	return res;2578}2579 2580static enum isl_tab_row_sign context_lex_ineq_sign(struct isl_context *context,2581			isl_int *ineq, int strict)2582{2583	struct isl_context_lex *clex = (struct isl_context_lex *)context;2584	return tab_ineq_sign(clex->tab, ineq, strict);2585}2586 2587/* Check whether "ineq" can be added to the tableau without rendering2588 * it infeasible.2589 */2590static int context_lex_test_ineq(struct isl_context *context, isl_int *ineq)2591{2592	struct isl_context_lex *clex = (struct isl_context_lex *)context;2593	struct isl_tab_undo *snap;2594	int feasible;2595 2596	if (!clex->tab)2597		return -1;2598 2599	if (isl_tab_extend_cons(clex->tab, 1) < 0)2600		return -1;2601 2602	snap = isl_tab_snap(clex->tab);2603	if (isl_tab_push_basis(clex->tab) < 0)2604		return -1;2605	clex->tab = add_lexmin_ineq(clex->tab, ineq);2606	clex->tab = check_integer_feasible(clex->tab);2607	if (!clex->tab)2608		return -1;2609	feasible = !clex->tab->empty;2610	if (isl_tab_rollback(clex->tab, snap) < 0)2611		return -1;2612 2613	return feasible;2614}2615 2616static int context_lex_get_div(struct isl_context *context, struct isl_tab *tab,2617		struct isl_vec *div)2618{2619	return get_div(tab, context, div);2620}2621 2622/* Insert a div specified by "div" to the context tableau at position "pos" and2623 * return isl_bool_true if the div is obviously non-negative.2624 * context_tab_add_div will always return isl_bool_true, because all variables2625 * in a isl_context_lex tableau are non-negative.2626 * However, if we are using a big parameter in the context, then this only2627 * reflects the non-negativity of the variable used to _encode_ the2628 * div, i.e., div' = M + div, so we can't draw any conclusions.2629 */2630static isl_bool context_lex_insert_div(struct isl_context *context, int pos,2631	__isl_keep isl_vec *div)2632{2633	struct isl_context_lex *clex = (struct isl_context_lex *)context;2634	isl_bool nonneg;2635	nonneg = context_tab_insert_div(clex->tab, pos, div,2636					context_lex_add_ineq_wrap, context);2637	if (nonneg < 0)2638		return isl_bool_error;2639	if (clex->tab->M)2640		return isl_bool_false;2641	return nonneg;2642}2643 2644static int context_lex_detect_equalities(struct isl_context *context,2645		struct isl_tab *tab)2646{2647	return 0;2648}2649 2650static int context_lex_best_split(struct isl_context *context,2651		struct isl_tab *tab)2652{2653	struct isl_context_lex *clex = (struct isl_context_lex *)context;2654	struct isl_tab_undo *snap;2655	int r;2656 2657	snap = isl_tab_snap(clex->tab);2658	if (isl_tab_push_basis(clex->tab) < 0)2659		return -1;2660	r = best_split(tab, clex->tab);2661 2662	if (r >= 0 && isl_tab_rollback(clex->tab, snap) < 0)2663		return -1;2664 2665	return r;2666}2667 2668static int context_lex_is_empty(struct isl_context *context)2669{2670	struct isl_context_lex *clex = (struct isl_context_lex *)context;2671	if (!clex->tab)2672		return -1;2673	return clex->tab->empty;2674}2675 2676static void *context_lex_save(struct isl_context *context)2677{2678	struct isl_context_lex *clex = (struct isl_context_lex *)context;2679	struct isl_tab_undo *snap;2680 2681	snap = isl_tab_snap(clex->tab);2682	if (isl_tab_push_basis(clex->tab) < 0)2683		return NULL;2684	if (isl_tab_save_samples(clex->tab) < 0)2685		return NULL;2686 2687	return snap;2688}2689 2690static void context_lex_restore(struct isl_context *context, void *save)2691{2692	struct isl_context_lex *clex = (struct isl_context_lex *)context;2693	if (isl_tab_rollback(clex->tab, (struct isl_tab_undo *)save) < 0) {2694		isl_tab_free(clex->tab);2695		clex->tab = NULL;2696	}2697}2698 2699static void context_lex_discard(void *save)2700{2701}2702 2703static int context_lex_is_ok(struct isl_context *context)2704{2705	struct isl_context_lex *clex = (struct isl_context_lex *)context;2706	return !!clex->tab;2707}2708 2709/* For each variable in the context tableau, check if the variable can2710 * only attain non-negative values.  If so, mark the parameter as non-negative2711 * in the main tableau.  This allows for a more direct identification of some2712 * cases of violated constraints.2713 */2714static struct isl_tab *tab_detect_nonnegative_parameters(struct isl_tab *tab,2715	struct isl_tab *context_tab)2716{2717	int i;2718	struct isl_tab_undo *snap;2719	struct isl_vec *ineq = NULL;2720	struct isl_tab_var *var;2721	int n;2722 2723	if (context_tab->n_var == 0)2724		return tab;2725 2726	ineq = isl_vec_alloc(tab->mat->ctx, 1 + context_tab->n_var);2727	if (!ineq)2728		goto error;2729 2730	if (isl_tab_extend_cons(context_tab, 1) < 0)2731		goto error;2732 2733	snap = isl_tab_snap(context_tab);2734 2735	n = 0;2736	isl_seq_clr(ineq->el, ineq->size);2737	for (i = 0; i < context_tab->n_var; ++i) {2738		isl_int_set_si(ineq->el[1 + i], 1);2739		if (isl_tab_add_ineq(context_tab, ineq->el) < 0)2740			goto error;2741		var = &context_tab->con[context_tab->n_con - 1];2742		if (!context_tab->empty &&2743		    !isl_tab_min_at_most_neg_one(context_tab, var)) {2744			int j = i;2745			if (i >= tab->n_param)2746				j = i - tab->n_param + tab->n_var - tab->n_div;2747			tab->var[j].is_nonneg = 1;2748			n++;2749		}2750		isl_int_set_si(ineq->el[1 + i], 0);2751		if (isl_tab_rollback(context_tab, snap) < 0)2752			goto error;2753	}2754 2755	if (context_tab->M && n == context_tab->n_var) {2756		context_tab->mat = isl_mat_drop_cols(context_tab->mat, 2, 1);2757		context_tab->M = 0;2758	}2759 2760	isl_vec_free(ineq);2761	return tab;2762error:2763	isl_vec_free(ineq);2764	isl_tab_free(tab);2765	return NULL;2766}2767 2768static struct isl_tab *context_lex_detect_nonnegative_parameters(2769	struct isl_context *context, struct isl_tab *tab)2770{2771	struct isl_context_lex *clex = (struct isl_context_lex *)context;2772	struct isl_tab_undo *snap;2773 2774	if (!tab)2775		return NULL;2776 2777	snap = isl_tab_snap(clex->tab);2778	if (isl_tab_push_basis(clex->tab) < 0)2779		goto error;2780 2781	tab = tab_detect_nonnegative_parameters(tab, clex->tab);2782 2783	if (isl_tab_rollback(clex->tab, snap) < 0)2784		goto error;2785 2786	return tab;2787error:2788	isl_tab_free(tab);2789	return NULL;2790}2791 2792static void context_lex_invalidate(struct isl_context *context)2793{2794	struct isl_context_lex *clex = (struct isl_context_lex *)context;2795	isl_tab_free(clex->tab);2796	clex->tab = NULL;2797}2798 2799static __isl_null struct isl_context *context_lex_free(2800	struct isl_context *context)2801{2802	struct isl_context_lex *clex = (struct isl_context_lex *)context;2803	isl_tab_free(clex->tab);2804	free(clex);2805 2806	return NULL;2807}2808 2809struct isl_context_op isl_context_lex_op = {2810	context_lex_detect_nonnegative_parameters,2811	context_lex_peek_basic_set,2812	context_lex_peek_tab,2813	context_lex_add_eq,2814	context_lex_add_ineq,2815	context_lex_ineq_sign,2816	context_lex_test_ineq,2817	context_lex_get_div,2818	context_lex_insert_div,2819	context_lex_detect_equalities,2820	context_lex_best_split,2821	context_lex_is_empty,2822	context_lex_is_ok,2823	context_lex_save,2824	context_lex_restore,2825	context_lex_discard,2826	context_lex_invalidate,2827	context_lex_free,2828};2829 2830static struct isl_tab *context_tab_for_lexmin(__isl_take isl_basic_set *bset)2831{2832	struct isl_tab *tab;2833 2834	if (!bset)2835		return NULL;2836	tab = tab_for_lexmin(bset_to_bmap(bset), NULL, 1, 0);2837	if (isl_tab_track_bset(tab, bset) < 0)2838		goto error;2839	tab = isl_tab_init_samples(tab);2840	return tab;2841error:2842	isl_tab_free(tab);2843	return NULL;2844}2845 2846static struct isl_context *isl_context_lex_alloc(struct isl_basic_set *dom)2847{2848	struct isl_context_lex *clex;2849 2850	if (!dom)2851		return NULL;2852 2853	clex = isl_alloc_type(dom->ctx, struct isl_context_lex);2854	if (!clex)2855		return NULL;2856 2857	clex->context.op = &isl_context_lex_op;2858 2859	clex->tab = context_tab_for_lexmin(isl_basic_set_copy(dom));2860	if (restore_lexmin(clex->tab) < 0)2861		goto error;2862	clex->tab = check_integer_feasible(clex->tab);2863	if (!clex->tab)2864		goto error;2865 2866	return &clex->context;2867error:2868	clex->context.op->free(&clex->context);2869	return NULL;2870}2871 2872/* Representation of the context when using generalized basis reduction.2873 *2874 * "shifted" contains the offsets of the unit hypercubes that lie inside the2875 * context.  Any rational point in "shifted" can therefore be rounded2876 * up to an integer point in the context.2877 * If the context is constrained by any equality, then "shifted" is not used2878 * as it would be empty.2879 */2880struct isl_context_gbr {2881	struct isl_context context;2882	struct isl_tab *tab;2883	struct isl_tab *shifted;2884	struct isl_tab *cone;2885};2886 2887static struct isl_tab *context_gbr_detect_nonnegative_parameters(2888	struct isl_context *context, struct isl_tab *tab)2889{2890	struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;2891	if (!tab)2892		return NULL;2893	return tab_detect_nonnegative_parameters(tab, cgbr->tab);2894}2895 2896static struct isl_basic_set *context_gbr_peek_basic_set(2897	struct isl_context *context)2898{2899	struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;2900	if (!cgbr->tab)2901		return NULL;2902	return isl_tab_peek_bset(cgbr->tab);2903}2904 2905static struct isl_tab *context_gbr_peek_tab(struct isl_context *context)2906{2907	struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;2908	return cgbr->tab;2909}2910 2911/* Initialize the "shifted" tableau of the context, which2912 * contains the constraints of the original tableau shifted2913 * by the sum of all negative coefficients.  This ensures2914 * that any rational point in the shifted tableau can2915 * be rounded up to yield an integer point in the original tableau.2916 */2917static void gbr_init_shifted(struct isl_context_gbr *cgbr)2918{2919	int i, j;2920	struct isl_vec *cst;2921	struct isl_basic_set *bset = isl_tab_peek_bset(cgbr->tab);2922	isl_size dim = isl_basic_set_dim(bset, isl_dim_all);2923 2924	if (dim < 0)2925		return;2926	cst = isl_vec_alloc(cgbr->tab->mat->ctx, bset->n_ineq);2927	if (!cst)2928		return;2929 2930	for (i = 0; i < bset->n_ineq; ++i) {2931		isl_int_set(cst->el[i], bset->ineq[i][0]);2932		for (j = 0; j < dim; ++j) {2933			if (!isl_int_is_neg(bset->ineq[i][1 + j]))2934				continue;2935			isl_int_add(bset->ineq[i][0], bset->ineq[i][0],2936				    bset->ineq[i][1 + j]);2937		}2938	}2939 2940	cgbr->shifted = isl_tab_from_basic_set(bset, 0);2941 2942	for (i = 0; i < bset->n_ineq; ++i)2943		isl_int_set(bset->ineq[i][0], cst->el[i]);2944 2945	isl_vec_free(cst);2946}2947 2948/* Check if the shifted tableau is non-empty, and if so2949 * use the sample point to construct an integer point2950 * of the context tableau.2951 */2952static struct isl_vec *gbr_get_shifted_sample(struct isl_context_gbr *cgbr)2953{2954	struct isl_vec *sample;2955 2956	if (!cgbr->shifted)2957		gbr_init_shifted(cgbr);2958	if (!cgbr->shifted)2959		return NULL;2960	if (cgbr->shifted->empty)2961		return isl_vec_alloc(cgbr->tab->mat->ctx, 0);2962 2963	sample = isl_tab_get_sample_value(cgbr->shifted);2964	sample = isl_vec_ceil(sample);2965 2966	return sample;2967}2968 2969static __isl_give isl_basic_set *drop_constant_terms(2970	__isl_take isl_basic_set *bset)2971{2972	int i;2973 2974	if (!bset)2975		return NULL;2976 2977	for (i = 0; i < bset->n_eq; ++i)2978		isl_int_set_si(bset->eq[i][0], 0);2979 2980	for (i = 0; i < bset->n_ineq; ++i)2981		isl_int_set_si(bset->ineq[i][0], 0);2982 2983	return bset;2984}2985 2986static int use_shifted(struct isl_context_gbr *cgbr)2987{2988	if (!cgbr->tab)2989		return 0;2990	return cgbr->tab->bmap->n_eq == 0 && cgbr->tab->bmap->n_div == 0;2991}2992 2993static struct isl_vec *gbr_get_sample(struct isl_context_gbr *cgbr)2994{2995	struct isl_basic_set *bset;2996	struct isl_basic_set *cone;2997 2998	if (isl_tab_sample_is_integer(cgbr->tab))2999		return isl_tab_get_sample_value(cgbr->tab);3000 3001	if (use_shifted(cgbr)) {3002		struct isl_vec *sample;3003 3004		sample = gbr_get_shifted_sample(cgbr);3005		if (!sample || sample->size > 0)3006			return sample;3007 3008		isl_vec_free(sample);3009	}3010 3011	if (!cgbr->cone) {3012		bset = isl_tab_peek_bset(cgbr->tab);3013		cgbr->cone = isl_tab_from_recession_cone(bset, 0);3014		if (!cgbr->cone)3015			return NULL;3016		if (isl_tab_track_bset(cgbr->cone,3017					isl_basic_set_copy(bset)) < 0)3018			return NULL;3019	}3020	if (isl_tab_detect_implicit_equalities(cgbr->cone) < 0)3021		return NULL;3022 3023	if (cgbr->cone->n_dead == cgbr->cone->n_col) {3024		struct isl_vec *sample;3025		struct isl_tab_undo *snap;3026 3027		if (cgbr->tab->basis) {3028			if (cgbr->tab->basis->n_col != 1 + cgbr->tab->n_var) {3029				isl_mat_free(cgbr->tab->basis);3030				cgbr->tab->basis = NULL;3031			}3032			cgbr->tab->n_zero = 0;3033			cgbr->tab->n_unbounded = 0;3034		}3035 3036		snap = isl_tab_snap(cgbr->tab);3037 3038		sample = isl_tab_sample(cgbr->tab);3039 3040		if (!sample || isl_tab_rollback(cgbr->tab, snap) < 0) {3041			isl_vec_free(sample);3042			return NULL;3043		}3044 3045		return sample;3046	}3047 3048	cone = isl_basic_set_dup(isl_tab_peek_bset(cgbr->cone));3049	cone = drop_constant_terms(cone);3050	cone = isl_basic_set_update_from_tab(cone, cgbr->cone);3051	cone = isl_basic_set_underlying_set(cone);3052	cone = isl_basic_set_gauss(cone, NULL);3053 3054	bset = isl_basic_set_dup(isl_tab_peek_bset(cgbr->tab));3055	bset = isl_basic_set_update_from_tab(bset, cgbr->tab);3056	bset = isl_basic_set_underlying_set(bset);3057	bset = isl_basic_set_gauss(bset, NULL);3058 3059	return isl_basic_set_sample_with_cone(bset, cone);3060}3061 3062static void check_gbr_integer_feasible(struct isl_context_gbr *cgbr)3063{3064	struct isl_vec *sample;3065 3066	if (!cgbr->tab)3067		return;3068 3069	if (cgbr->tab->empty)3070		return;3071 3072	sample = gbr_get_sample(cgbr);3073	if (!sample)3074		goto error;3075 3076	if (sample->size == 0) {3077		isl_vec_free(sample);3078		if (isl_tab_mark_empty(cgbr->tab) < 0)3079			goto error;3080		return;3081	}3082 3083	if (isl_tab_add_sample(cgbr->tab, sample) < 0)3084		goto error;3085 3086	return;3087error:3088	isl_tab_free(cgbr->tab);3089	cgbr->tab = NULL;3090}3091 3092static struct isl_tab *add_gbr_eq(struct isl_tab *tab, isl_int *eq)3093{3094	if (!tab)3095		return NULL;3096 3097	if (isl_tab_extend_cons(tab, 2) < 0)3098		goto error;3099 3100	if (isl_tab_add_eq(tab, eq) < 0)3101		goto error;3102 3103	return tab;3104error:3105	isl_tab_free(tab);3106	return NULL;3107}3108 3109/* Add the equality described by "eq" to the context.3110 * If "check" is set, then we check if the context is empty after3111 * adding the equality.3112 * If "update" is set, then we check if the samples are still valid.3113 *3114 * We do not explicitly add shifted copies of the equality to3115 * cgbr->shifted since they would conflict with each other.3116 * Instead, we directly mark cgbr->shifted empty.3117 */3118static void context_gbr_add_eq(struct isl_context *context, isl_int *eq,3119		int check, int update)3120{3121	struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;3122 3123	cgbr->tab = add_gbr_eq(cgbr->tab, eq);3124 3125	if (cgbr->shifted && !cgbr->shifted->empty && use_shifted(cgbr)) {3126		if (isl_tab_mark_empty(cgbr->shifted) < 0)3127			goto error;3128	}3129 3130	if (cgbr->cone && cgbr->cone->n_col != cgbr->cone->n_dead) {3131		if (isl_tab_extend_cons(cgbr->cone, 2) < 0)3132			goto error;3133		if (isl_tab_add_eq(cgbr->cone, eq) < 0)3134			goto error;3135	}3136 3137	if (check) {3138		int v = tab_has_valid_sample(cgbr->tab, eq, 1);3139		if (v < 0)3140			goto error;3141		if (!v)3142			check_gbr_integer_feasible(cgbr);3143	}3144	if (update)3145		cgbr->tab = check_samples(cgbr->tab, eq, 1);3146	return;3147error:3148	isl_tab_free(cgbr->tab);3149	cgbr->tab = NULL;3150}3151 3152static void add_gbr_ineq(struct isl_context_gbr *cgbr, isl_int *ineq)3153{3154	if (!cgbr->tab)3155		return;3156 3157	if (isl_tab_extend_cons(cgbr->tab, 1) < 0)3158		goto error;3159 3160	if (isl_tab_add_ineq(cgbr->tab, ineq) < 0)3161		goto error;3162 3163	if (cgbr->shifted && !cgbr->shifted->empty && use_shifted(cgbr)) {3164		int i;3165		isl_size dim;3166		dim = isl_basic_map_dim(cgbr->tab->bmap, isl_dim_all);3167		if (dim < 0)3168			goto error;3169 3170		if (isl_tab_extend_cons(cgbr->shifted, 1) < 0)3171			goto error;3172 3173		for (i = 0; i < dim; ++i) {3174			if (!isl_int_is_neg(ineq[1 + i]))3175				continue;3176			isl_int_add(ineq[0], ineq[0], ineq[1 + i]);3177		}3178 3179		if (isl_tab_add_ineq(cgbr->shifted, ineq) < 0)3180			goto error;3181 3182		for (i = 0; i < dim; ++i) {3183			if (!isl_int_is_neg(ineq[1 + i]))3184				continue;3185			isl_int_sub(ineq[0], ineq[0], ineq[1 + i]);3186		}3187	}3188 3189	if (cgbr->cone && cgbr->cone->n_col != cgbr->cone->n_dead) {3190		if (isl_tab_extend_cons(cgbr->cone, 1) < 0)3191			goto error;3192		if (isl_tab_add_ineq(cgbr->cone, ineq) < 0)3193			goto error;3194	}3195 3196	return;3197error:3198	isl_tab_free(cgbr->tab);3199	cgbr->tab = NULL;3200}3201 3202static void context_gbr_add_ineq(struct isl_context *context, isl_int *ineq,3203		int check, int update)3204{3205	struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;3206 3207	add_gbr_ineq(cgbr, ineq);3208	if (!cgbr->tab)3209		return;3210 3211	if (check) {3212		int v = tab_has_valid_sample(cgbr->tab, ineq, 0);3213		if (v < 0)3214			goto error;3215		if (!v)3216			check_gbr_integer_feasible(cgbr);3217	}3218	if (update)3219		cgbr->tab = check_samples(cgbr->tab, ineq, 0);3220	return;3221error:3222	isl_tab_free(cgbr->tab);3223	cgbr->tab = NULL;3224}3225 3226static isl_stat context_gbr_add_ineq_wrap(void *user, isl_int *ineq)3227{3228	struct isl_context *context = (struct isl_context *)user;3229	context_gbr_add_ineq(context, ineq, 0, 0);3230	return context->op->is_ok(context) ? isl_stat_ok : isl_stat_error;3231}3232 3233static enum isl_tab_row_sign context_gbr_ineq_sign(struct isl_context *context,3234			isl_int *ineq, int strict)3235{3236	struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;3237	return tab_ineq_sign(cgbr->tab, ineq, strict);3238}3239 3240/* Check whether "ineq" can be added to the tableau without rendering3241 * it infeasible.3242 */3243static int context_gbr_test_ineq(struct isl_context *context, isl_int *ineq)3244{3245	struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;3246	struct isl_tab_undo *snap;3247	struct isl_tab_undo *shifted_snap = NULL;3248	struct isl_tab_undo *cone_snap = NULL;3249	int feasible;3250 3251	if (!cgbr->tab)3252		return -1;3253 3254	if (isl_tab_extend_cons(cgbr->tab, 1) < 0)3255		return -1;3256 3257	snap = isl_tab_snap(cgbr->tab);3258	if (cgbr->shifted)3259		shifted_snap = isl_tab_snap(cgbr->shifted);3260	if (cgbr->cone)3261		cone_snap = isl_tab_snap(cgbr->cone);3262	add_gbr_ineq(cgbr, ineq);3263	check_gbr_integer_feasible(cgbr);3264	if (!cgbr->tab)3265		return -1;3266	feasible = !cgbr->tab->empty;3267	if (isl_tab_rollback(cgbr->tab, snap) < 0)3268		return -1;3269	if (shifted_snap) {3270		if (isl_tab_rollback(cgbr->shifted, shifted_snap))3271			return -1;3272	} else if (cgbr->shifted) {3273		isl_tab_free(cgbr->shifted);3274		cgbr->shifted = NULL;3275	}3276	if (cone_snap) {3277		if (isl_tab_rollback(cgbr->cone, cone_snap))3278			return -1;3279	} else if (cgbr->cone) {3280		isl_tab_free(cgbr->cone);3281		cgbr->cone = NULL;3282	}3283 3284	return feasible;3285}3286 3287/* Return the column of the last of the variables associated to3288 * a column that has a non-zero coefficient.3289 * This function is called in a context where only coefficients3290 * of parameters or divs can be non-zero.3291 */3292static int last_non_zero_var_col(struct isl_tab *tab, isl_int *p)3293{3294	int i;3295	int col;3296 3297	if (tab->n_var == 0)3298		return -1;3299 3300	for (i = tab->n_var - 1; i >= 0; --i) {3301		if (i >= tab->n_param && i < tab->n_var - tab->n_div)3302			continue;3303		if (tab->var[i].is_row)3304			continue;3305		col = tab->var[i].index;3306		if (!isl_int_is_zero(p[col]))3307			return col;3308	}3309 3310	return -1;3311}3312 3313/* Look through all the recently added equalities in the context3314 * to see if we can propagate any of them to the main tableau.3315 *3316 * The newly added equalities in the context are encoded as pairs3317 * of inequalities starting at inequality "first".3318 *3319 * We tentatively add each of these equalities to the main tableau3320 * and if this happens to result in a row with a final coefficient3321 * that is one or negative one, we use it to kill a column3322 * in the main tableau.  Otherwise, we discard the tentatively3323 * added row.3324 * This tentative addition of equality constraints turns3325 * on the undo facility of the tableau.  Turn it off again3326 * at the end, assuming it was turned off to begin with.3327 *3328 * Return 0 on success and -1 on failure.3329 */3330static int propagate_equalities(struct isl_context_gbr *cgbr,3331	struct isl_tab *tab, unsigned first)3332{3333	int i;3334	struct isl_vec *eq = NULL;3335	isl_bool needs_undo;3336 3337	needs_undo = isl_tab_need_undo(tab);3338	if (needs_undo < 0)3339		goto error;3340	eq = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);3341	if (!eq)3342		goto error;3343 3344	if (isl_tab_extend_cons(tab, (cgbr->tab->bmap->n_ineq - first)/2) < 0)3345		goto error;3346 3347	isl_seq_clr(eq->el + 1 + tab->n_param,3348		    tab->n_var - tab->n_param - tab->n_div);3349	for (i = first; i < cgbr->tab->bmap->n_ineq; i += 2) {3350		int j;3351		int r;3352		struct isl_tab_undo *snap;3353		snap = isl_tab_snap(tab);3354 3355		isl_seq_cpy(eq->el, cgbr->tab->bmap->ineq[i], 1 + tab->n_param);3356		isl_seq_cpy(eq->el + 1 + tab->n_var - tab->n_div,3357			    cgbr->tab->bmap->ineq[i] + 1 + tab->n_param,3358			    tab->n_div);3359 3360		r = isl_tab_add_row(tab, eq->el);3361		if (r < 0)3362			goto error;3363		r = tab->con[r].index;3364		j = last_non_zero_var_col(tab, tab->mat->row[r] + 2 + tab->M);3365		if (j < 0 || j < tab->n_dead ||3366		    !isl_int_is_one(tab->mat->row[r][0]) ||3367		    (!isl_int_is_one(tab->mat->row[r][2 + tab->M + j]) &&3368		     !isl_int_is_negone(tab->mat->row[r][2 + tab->M + j]))) {3369			if (isl_tab_rollback(tab, snap) < 0)3370				goto error;3371			continue;3372		}3373		if (isl_tab_pivot(tab, r, j) < 0)3374			goto error;3375		if (isl_tab_kill_col(tab, j) < 0)3376			goto error;3377 3378		if (restore_lexmin(tab) < 0)3379			goto error;3380	}3381 3382	if (!needs_undo)3383		isl_tab_clear_undo(tab);3384	isl_vec_free(eq);3385 3386	return 0;3387error:3388	isl_vec_free(eq);3389	isl_tab_free(cgbr->tab);3390	cgbr->tab = NULL;3391	return -1;3392}3393 3394static int context_gbr_detect_equalities(struct isl_context *context,3395	struct isl_tab *tab)3396{3397	struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;3398	unsigned n_ineq;3399 3400	if (!cgbr->cone) {3401		struct isl_basic_set *bset = isl_tab_peek_bset(cgbr->tab);3402		cgbr->cone = isl_tab_from_recession_cone(bset, 0);3403		if (!cgbr->cone)3404			goto error;3405		if (isl_tab_track_bset(cgbr->cone,3406					isl_basic_set_copy(bset)) < 0)3407			goto error;3408	}3409	if (isl_tab_detect_implicit_equalities(cgbr->cone) < 0)3410		goto error;3411 3412	n_ineq = cgbr->tab->bmap->n_ineq;3413	cgbr->tab = isl_tab_detect_equalities(cgbr->tab, cgbr->cone);3414	if (!cgbr->tab)3415		return -1;3416	if (cgbr->tab->bmap->n_ineq > n_ineq &&3417	    propagate_equalities(cgbr, tab, n_ineq) < 0)3418		return -1;3419 3420	return 0;3421error:3422	isl_tab_free(cgbr->tab);3423	cgbr->tab = NULL;3424	return -1;3425}3426 3427static int context_gbr_get_div(struct isl_context *context, struct isl_tab *tab,3428		struct isl_vec *div)3429{3430	return get_div(tab, context, div);3431}3432 3433static isl_bool context_gbr_insert_div(struct isl_context *context, int pos,3434	__isl_keep isl_vec *div)3435{3436	struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;3437	if (cgbr->cone) {3438		int r, o_div;3439		isl_size n_div;3440 3441		n_div = isl_basic_map_dim(cgbr->cone->bmap, isl_dim_div);3442		if (n_div < 0)3443			return isl_bool_error;3444		o_div = cgbr->cone->n_var - n_div;3445 3446		if (isl_tab_extend_cons(cgbr->cone, 3) < 0)3447			return isl_bool_error;3448		if (isl_tab_extend_vars(cgbr->cone, 1) < 0)3449			return isl_bool_error;3450		if ((r = isl_tab_insert_var(cgbr->cone, pos)) <0)3451			return isl_bool_error;3452 3453		cgbr->cone->bmap = isl_basic_map_insert_div(cgbr->cone->bmap,3454						    r - o_div, div);3455		if (!cgbr->cone->bmap)3456			return isl_bool_error;3457		if (isl_tab_push_var(cgbr->cone, isl_tab_undo_bmap_div,3458				    &cgbr->cone->var[r]) < 0)3459			return isl_bool_error;3460	}3461	return context_tab_insert_div(cgbr->tab, pos, div,3462					context_gbr_add_ineq_wrap, context);3463}3464 3465static int context_gbr_best_split(struct isl_context *context,3466		struct isl_tab *tab)3467{3468	struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;3469	struct isl_tab_undo *snap;3470	int r;3471 3472	snap = isl_tab_snap(cgbr->tab);3473	r = best_split(tab, cgbr->tab);3474 3475	if (r >= 0 && isl_tab_rollback(cgbr->tab, snap) < 0)3476		return -1;3477 3478	return r;3479}3480 3481static int context_gbr_is_empty(struct isl_context *context)3482{3483	struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;3484	if (!cgbr->tab)3485		return -1;3486	return cgbr->tab->empty;3487}3488 3489struct isl_gbr_tab_undo {3490	struct isl_tab_undo *tab_snap;3491	struct isl_tab_undo *shifted_snap;3492	struct isl_tab_undo *cone_snap;3493};3494 3495static void *context_gbr_save(struct isl_context *context)3496{3497	struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;3498	struct isl_gbr_tab_undo *snap;3499 3500	if (!cgbr->tab)3501		return NULL;3502 3503	snap = isl_alloc_type(cgbr->tab->mat->ctx, struct isl_gbr_tab_undo);3504	if (!snap)3505		return NULL;3506 3507	snap->tab_snap = isl_tab_snap(cgbr->tab);3508	if (isl_tab_save_samples(cgbr->tab) < 0)3509		goto error;3510 3511	if (cgbr->shifted)3512		snap->shifted_snap = isl_tab_snap(cgbr->shifted);3513	else3514		snap->shifted_snap = NULL;3515 3516	if (cgbr->cone)3517		snap->cone_snap = isl_tab_snap(cgbr->cone);3518	else3519		snap->cone_snap = NULL;3520 3521	return snap;3522error:3523	free(snap);3524	return NULL;3525}3526 3527static void context_gbr_restore(struct isl_context *context, void *save)3528{3529	struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;3530	struct isl_gbr_tab_undo *snap = (struct isl_gbr_tab_undo *)save;3531	if (!snap)3532		goto error;3533	if (isl_tab_rollback(cgbr->tab, snap->tab_snap) < 0)3534		goto error;3535 3536	if (snap->shifted_snap) {3537		if (isl_tab_rollback(cgbr->shifted, snap->shifted_snap) < 0)3538			goto error;3539	} else if (cgbr->shifted) {3540		isl_tab_free(cgbr->shifted);3541		cgbr->shifted = NULL;3542	}3543 3544	if (snap->cone_snap) {3545		if (isl_tab_rollback(cgbr->cone, snap->cone_snap) < 0)3546			goto error;3547	} else if (cgbr->cone) {3548		isl_tab_free(cgbr->cone);3549		cgbr->cone = NULL;3550	}3551 3552	free(snap);3553 3554	return;3555error:3556	free(snap);3557	isl_tab_free(cgbr->tab);3558	cgbr->tab = NULL;3559}3560 3561static void context_gbr_discard(void *save)3562{3563	struct isl_gbr_tab_undo *snap = (struct isl_gbr_tab_undo *)save;3564	free(snap);3565}3566 3567static int context_gbr_is_ok(struct isl_context *context)3568{3569	struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;3570	return !!cgbr->tab;3571}3572 3573static void context_gbr_invalidate(struct isl_context *context)3574{3575	struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;3576	isl_tab_free(cgbr->tab);3577	cgbr->tab = NULL;3578}3579 3580static __isl_null struct isl_context *context_gbr_free(3581	struct isl_context *context)3582{3583	struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;3584	isl_tab_free(cgbr->tab);3585	isl_tab_free(cgbr->shifted);3586	isl_tab_free(cgbr->cone);3587	free(cgbr);3588 3589	return NULL;3590}3591 3592struct isl_context_op isl_context_gbr_op = {3593	context_gbr_detect_nonnegative_parameters,3594	context_gbr_peek_basic_set,3595	context_gbr_peek_tab,3596	context_gbr_add_eq,3597	context_gbr_add_ineq,3598	context_gbr_ineq_sign,3599	context_gbr_test_ineq,3600	context_gbr_get_div,3601	context_gbr_insert_div,3602	context_gbr_detect_equalities,3603	context_gbr_best_split,3604	context_gbr_is_empty,3605	context_gbr_is_ok,3606	context_gbr_save,3607	context_gbr_restore,3608	context_gbr_discard,3609	context_gbr_invalidate,3610	context_gbr_free,3611};3612 3613static struct isl_context *isl_context_gbr_alloc(__isl_keep isl_basic_set *dom)3614{3615	struct isl_context_gbr *cgbr;3616 3617	if (!dom)3618		return NULL;3619 3620	cgbr = isl_calloc_type(dom->ctx, struct isl_context_gbr);3621	if (!cgbr)3622		return NULL;3623 3624	cgbr->context.op = &isl_context_gbr_op;3625 3626	cgbr->shifted = NULL;3627	cgbr->cone = NULL;3628	cgbr->tab = isl_tab_from_basic_set(dom, 1);3629	cgbr->tab = isl_tab_init_samples(cgbr->tab);3630	if (!cgbr->tab)3631		goto error;3632	check_gbr_integer_feasible(cgbr);3633 3634	return &cgbr->context;3635error:3636	cgbr->context.op->free(&cgbr->context);3637	return NULL;3638}3639 3640/* Allocate a context corresponding to "dom".3641 * The representation specific fields are initialized by3642 * isl_context_lex_alloc or isl_context_gbr_alloc.3643 * The shared "n_unknown" field is initialized to the number3644 * of final unknown integer divisions in "dom".3645 */3646static struct isl_context *isl_context_alloc(__isl_keep isl_basic_set *dom)3647{3648	struct isl_context *context;3649	int first;3650	isl_size n_div;3651 3652	if (!dom)3653		return NULL;3654 3655	if (dom->ctx->opt->context == ISL_CONTEXT_LEXMIN)3656		context = isl_context_lex_alloc(dom);3657	else3658		context = isl_context_gbr_alloc(dom);3659 3660	if (!context)3661		return NULL;3662 3663	first = isl_basic_set_first_unknown_div(dom);3664	n_div = isl_basic_set_dim(dom, isl_dim_div);3665	if (first < 0 || n_div < 0)3666		return context->op->free(context);3667	context->n_unknown = n_div - first;3668 3669	return context;3670}3671 3672/* Initialize some common fields of "sol", which keeps track3673 * of the solution of an optimization problem on "bmap" over3674 * the domain "dom".3675 * If "max" is set, then a maximization problem is being solved, rather than3676 * a minimization problem, which means that the variables in the3677 * tableau have value "M - x" rather than "M + x".3678 */3679static isl_stat sol_init(struct isl_sol *sol, __isl_keep isl_basic_map *bmap,3680	__isl_keep isl_basic_set *dom, int max)3681{3682	sol->rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);3683	sol->dec_level.callback.run = &sol_dec_level_wrap;3684	sol->dec_level.sol = sol;3685	sol->max = max;3686	sol->n_out = isl_basic_map_dim(bmap, isl_dim_out);3687	sol->space = isl_basic_map_get_space(bmap);3688 3689	sol->context = isl_context_alloc(dom);3690	if (sol->n_out < 0 || !sol->space || !sol->context)3691		return isl_stat_error;3692 3693	return isl_stat_ok;3694}3695 3696/* Construct an isl_sol_map structure for accumulating the solution.3697 * If track_empty is set, then we also keep track of the parts3698 * of the context where there is no solution.3699 * If max is set, then we are solving a maximization, rather than3700 * a minimization problem, which means that the variables in the3701 * tableau have value "M - x" rather than "M + x".3702 */3703static struct isl_sol *sol_map_init(__isl_keep isl_basic_map *bmap,3704	__isl_take isl_basic_set *dom, int track_empty, int max)3705{3706	struct isl_sol_map *sol_map = NULL;3707	isl_space *space;3708 3709	if (!bmap)3710		goto error;3711 3712	sol_map = isl_calloc_type(bmap->ctx, struct isl_sol_map);3713	if (!sol_map)3714		goto error;3715 3716	sol_map->sol.free = &sol_map_free;3717	if (sol_init(&sol_map->sol, bmap, dom, max) < 0)3718		goto error;3719	sol_map->sol.add = &sol_map_add_wrap;3720	sol_map->sol.add_empty = track_empty ? &sol_map_add_empty_wrap : NULL;3721	space = isl_space_copy(sol_map->sol.space);3722	sol_map->map = isl_map_alloc_space(space, 1, ISL_MAP_DISJOINT);3723	if (!sol_map->map)3724		goto error;3725 3726	if (track_empty) {3727		sol_map->empty = isl_set_alloc_space(isl_basic_set_get_space(dom),3728							1, ISL_SET_DISJOINT);3729		if (!sol_map->empty)3730			goto error;3731	}3732 3733	isl_basic_set_free(dom);3734	return &sol_map->sol;3735error:3736	isl_basic_set_free(dom);3737	sol_free(&sol_map->sol);3738	return NULL;3739}3740 3741/* Check whether all coefficients of (non-parameter) variables3742 * are non-positive, meaning that no pivots can be performed on the row.3743 */3744static int is_critical(struct isl_tab *tab, int row)3745{3746	int j;3747	unsigned off = 2 + tab->M;3748 3749	for (j = tab->n_dead; j < tab->n_col; ++j) {3750		if (col_is_parameter_var(tab, j))3751			continue;3752 3753		if (isl_int_is_pos(tab->mat->row[row][off + j]))3754			return 0;3755	}3756 3757	return 1;3758}3759 3760/* Check whether the inequality represented by vec is strict over the integers,3761 * i.e., there are no integer values satisfying the constraint with3762 * equality.  This happens if the gcd of the coefficients is not a divisor3763 * of the constant term.  If so, scale the constraint down by the gcd3764 * of the coefficients.3765 */3766static int is_strict(struct isl_vec *vec)3767{3768	isl_int gcd;3769	int strict = 0;3770 3771	isl_int_init(gcd);3772	isl_seq_gcd(vec->el + 1, vec->size - 1, &gcd);3773	if (!isl_int_is_one(gcd)) {3774		strict = !isl_int_is_divisible_by(vec->el[0], gcd);3775		isl_int_fdiv_q(vec->el[0], vec->el[0], gcd);3776		isl_seq_scale_down(vec->el + 1, vec->el + 1, gcd, vec->size-1);3777	}3778	isl_int_clear(gcd);3779 3780	return strict;3781}3782 3783/* Determine the sign of the given row of the main tableau.3784 * The result is one of3785 *	isl_tab_row_pos: always non-negative; no pivot needed3786 *	isl_tab_row_neg: always non-positive; pivot3787 *	isl_tab_row_any: can be both positive and negative; split3788 *3789 * We first handle some simple cases3790 *	- the row sign may be known already3791 *	- the row may be obviously non-negative3792 *	- the parametric constant may be equal to that of another row3793 *	  for which we know the sign.  This sign will be either "pos" or3794 *	  "any".  If it had been "neg" then we would have pivoted before.3795 *3796 * If none of these cases hold, we check the value of the row for each3797 * of the currently active samples.  Based on the signs of these values3798 * we make an initial determination of the sign of the row.3799 *3800 *	all zero			->	unk(nown)3801 *	all non-negative		->	pos3802 *	all non-positive		->	neg3803 *	both negative and positive	->	all3804 *3805 * If we end up with "all", we are done.3806 * Otherwise, we perform a check for positive and/or negative3807 * values as follows.3808 *3809 *	samples	       neg	       unk	       pos3810 *	<0 ?			    Y        N	    Y        N3811 *					    pos    any      pos3812 *	>0 ?	     Y      N	 Y     N3813 *		    any    neg  any   neg3814 *3815 * There is no special sign for "zero", because we can usually treat zero3816 * as either non-negative or non-positive, whatever works out best.3817 * However, if the row is "critical", meaning that pivoting is impossible3818 * then we don't want to limp zero with the non-positive case, because3819 * then we we would lose the solution for those values of the parameters3820 * where the value of the row is zero.  Instead, we treat 0 as non-negative3821 * ensuring a split if the row can attain both zero and negative values.3822 * The same happens when the original constraint was one that could not3823 * be satisfied with equality by any integer values of the parameters.3824 * In this case, we normalize the constraint, but then a value of zero3825 * for the normalized constraint is actually a positive value for the3826 * original constraint, so again we need to treat zero as non-negative.3827 * In both these cases, we have the following decision tree instead:3828 *3829 *	all non-negative		->	pos3830 *	all negative			->	neg3831 *	both negative and non-negative	->	all3832 *3833 *	samples	       neg	          	       pos3834 *	<0 ?			             	    Y        N3835 *					           any      pos3836 *	>=0 ?	     Y      N3837 *		    any    neg3838 */3839static enum isl_tab_row_sign row_sign(struct isl_tab *tab,3840	struct isl_sol *sol, int row)3841{3842	struct isl_vec *ineq = NULL;3843	enum isl_tab_row_sign res = isl_tab_row_unknown;3844	int critical;3845	int strict;3846	int row2;3847 3848	if (tab->row_sign[row] != isl_tab_row_unknown)3849		return tab->row_sign[row];3850	if (is_obviously_nonneg(tab, row))3851		return isl_tab_row_pos;3852	for (row2 = tab->n_redundant; row2 < tab->n_row; ++row2) {3853		if (tab->row_sign[row2] == isl_tab_row_unknown)3854			continue;3855		if (identical_parameter_line(tab, row, row2))3856			return tab->row_sign[row2];3857	}3858 3859	critical = is_critical(tab, row);3860 3861	ineq = get_row_parameter_ineq(tab, row);3862	if (!ineq)3863		goto error;3864 3865	strict = is_strict(ineq);3866 3867	res = sol->context->op->ineq_sign(sol->context, ineq->el,3868					  critical || strict);3869 3870	if (res == isl_tab_row_unknown || res == isl_tab_row_pos) {3871		/* test for negative values */3872		int feasible;3873		isl_seq_neg(ineq->el, ineq->el, ineq->size);3874		isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);3875 3876		feasible = sol->context->op->test_ineq(sol->context, ineq->el);3877		if (feasible < 0)3878			goto error;3879		if (!feasible)3880			res = isl_tab_row_pos;3881		else3882			res = (res == isl_tab_row_unknown) ? isl_tab_row_neg3883							   : isl_tab_row_any;3884		if (res == isl_tab_row_neg) {3885			isl_seq_neg(ineq->el, ineq->el, ineq->size);3886			isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);3887		}3888	}3889 3890	if (res == isl_tab_row_neg) {3891		/* test for positive values */3892		int feasible;3893		if (!critical && !strict)3894			isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);3895 3896		feasible = sol->context->op->test_ineq(sol->context, ineq->el);3897		if (feasible < 0)3898			goto error;3899		if (feasible)3900			res = isl_tab_row_any;3901	}3902 3903	isl_vec_free(ineq);3904	return res;3905error:3906	isl_vec_free(ineq);3907	return isl_tab_row_unknown;3908}3909 3910static void find_solutions(struct isl_sol *sol, struct isl_tab *tab);3911 3912/* Find solutions for values of the parameters that satisfy the given3913 * inequality.3914 *3915 * We currently take a snapshot of the context tableau that is reset3916 * when we return from this function, while we make a copy of the main3917 * tableau, leaving the original main tableau untouched.3918 * These are fairly arbitrary choices.  Making a copy also of the context3919 * tableau would obviate the need to undo any changes made to it later,3920 * while taking a snapshot of the main tableau could reduce memory usage.3921 * If we were to switch to taking a snapshot of the main tableau,3922 * we would have to keep in mind that we need to save the row signs3923 * and that we need to do this before saving the current basis3924 * such that the basis has been restore before we restore the row signs.3925 */3926static void find_in_pos(struct isl_sol *sol, struct isl_tab *tab, isl_int *ineq)3927{3928	void *saved;3929 3930	if (!sol->context)3931		goto error;3932	saved = sol->context->op->save(sol->context);3933 3934	tab = isl_tab_dup(tab);3935	if (!tab)3936		goto error;3937 3938	sol->context->op->add_ineq(sol->context, ineq, 0, 1);3939 3940	find_solutions(sol, tab);3941 3942	if (!sol->error)3943		sol->context->op->restore(sol->context, saved);3944	else3945		sol->context->op->discard(saved);3946	return;3947error:3948	sol->error = 1;3949}3950 3951/* Record the absence of solutions for those values of the parameters3952 * that do not satisfy the given inequality with equality.3953 */3954static void no_sol_in_strict(struct isl_sol *sol,3955	struct isl_tab *tab, struct isl_vec *ineq)3956{3957	int empty;3958	void *saved;3959 3960	if (!sol->context || sol->error)3961		goto error;3962	saved = sol->context->op->save(sol->context);3963 3964	isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);3965 3966	sol->context->op->add_ineq(sol->context, ineq->el, 1, 0);3967	if (!sol->context)3968		goto error;3969 3970	empty = tab->empty;3971	tab->empty = 1;3972	sol_add(sol, tab);3973	tab->empty = empty;3974 3975	isl_int_add_ui(ineq->el[0], ineq->el[0], 1);3976 3977	sol->context->op->restore(sol->context, saved);3978	return;3979error:3980	sol->error = 1;3981}3982 3983/* Reset all row variables that are marked to have a sign that may3984 * be both positive and negative to have an unknown sign.3985 */3986static void reset_any_to_unknown(struct isl_tab *tab)3987{3988	int row;3989 3990	for (row = tab->n_redundant; row < tab->n_row; ++row) {3991		if (!isl_tab_var_from_row(tab, row)->is_nonneg)3992			continue;3993		if (tab->row_sign[row] == isl_tab_row_any)3994			tab->row_sign[row] = isl_tab_row_unknown;3995	}3996}3997 3998/* Compute the lexicographic minimum of the set represented by the main3999 * tableau "tab" within the context "sol->context_tab".4000 * On entry the sample value of the main tableau is lexicographically4001 * less than or equal to this lexicographic minimum.4002 * Pivots are performed until a feasible point is found, which is then4003 * necessarily equal to the minimum, or until the tableau is found to4004 * be infeasible.  Some pivots may need to be performed for only some4005 * feasible values of the context tableau.  If so, the context tableau4006 * is split into a part where the pivot is needed and a part where it is not.4007 *4008 * Whenever we enter the main loop, the main tableau is such that no4009 * "obvious" pivots need to be performed on it, where "obvious" means4010 * that the given row can be seen to be negative without looking at4011 * the context tableau.  In particular, for non-parametric problems,4012 * no pivots need to be performed on the main tableau.4013 * The caller of find_solutions is responsible for making this property4014 * hold prior to the first iteration of the loop, while restore_lexmin4015 * is called before every other iteration.4016 *4017 * Inside the main loop, we first examine the signs of the rows of4018 * the main tableau within the context of the context tableau.4019 * If we find a row that is always non-positive for all values of4020 * the parameters satisfying the context tableau and negative for at4021 * least one value of the parameters, we perform the appropriate pivot4022 * and start over.  An exception is the case where no pivot can be4023 * performed on the row.  In this case, we require that the sign of4024 * the row is negative for all values of the parameters (rather than just4025 * non-positive).  This special case is handled inside row_sign, which4026 * will say that the row can have any sign if it determines that it can4027 * attain both negative and zero values.4028 *4029 * If we can't find a row that always requires a pivot, but we can find4030 * one or more rows that require a pivot for some values of the parameters4031 * (i.e., the row can attain both positive and negative signs), then we split4032 * the context tableau into two parts, one where we force the sign to be4033 * non-negative and one where we force is to be negative.4034 * The non-negative part is handled by a recursive call (through find_in_pos).4035 * Upon returning from this call, we continue with the negative part and4036 * perform the required pivot.4037 *4038 * If no such rows can be found, all rows are non-negative and we have4039 * found a (rational) feasible point.  If we only wanted a rational point4040 * then we are done.4041 * Otherwise, we check if all values of the sample point of the tableau4042 * are integral for the variables.  If so, we have found the minimal4043 * integral point and we are done.4044 * If the sample point is not integral, then we need to make a distinction4045 * based on whether the constant term is non-integral or the coefficients4046 * of the parameters.  Furthermore, in order to decide how to handle4047 * the non-integrality, we also need to know whether the coefficients4048 * of the other columns in the tableau are integral.  This leads4049 * to the following table.  The first two rows do not correspond4050 * to a non-integral sample point and are only mentioned for completeness.4051 *4052 *	constant	parameters	other4053 *4054 *	int		int		int	|4055 *	int		int		rat	| -> no problem4056 *4057 *	rat		int		int	  -> fail4058 *4059 *	rat		int		rat	  -> cut4060 *4061 *	int		rat		rat	|4062 *	rat		rat		rat	| -> parametric cut4063 *4064 *	int		rat		int	|4065 *	rat		rat		int	| -> split context4066 *4067 * If the parametric constant is completely integral, then there is nothing4068 * to be done.  If the constant term is non-integral, but all the other4069 * coefficient are integral, then there is nothing that can be done4070 * and the tableau has no integral solution.4071 * If, on the other hand, one or more of the other columns have rational4072 * coefficients, but the parameter coefficients are all integral, then4073 * we can perform a regular (non-parametric) cut.4074 * Finally, if there is any parameter coefficient that is non-integral,4075 * then we need to involve the context tableau.  There are two cases here.4076 * If at least one other column has a rational coefficient, then we4077 * can perform a parametric cut in the main tableau by adding a new4078 * integer division in the context tableau.4079 * If all other columns have integral coefficients, then we need to4080 * enforce that the rational combination of parameters (c + \sum a_i y_i)/m4081 * is always integral.  We do this by introducing an integer division4082 * q = floor((c + \sum a_i y_i)/m) and stipulating that its argument should4083 * always be integral in the context tableau, i.e., m q = c + \sum a_i y_i.4084 * Since q is expressed in the tableau as4085 *	c + \sum a_i y_i - m q >= 04086 *	-c - \sum a_i y_i + m q + m - 1 >= 04087 * it is sufficient to add the inequality4088 *	-c - \sum a_i y_i + m q >= 04089 * In the part of the context where this inequality does not hold, the4090 * main tableau is marked as being empty.4091 */4092static void find_solutions(struct isl_sol *sol, struct isl_tab *tab)4093{4094	struct isl_context *context;4095	int r;4096 4097	if (!tab || sol->error)4098		goto error;4099 4100	context = sol->context;4101 4102	if (tab->empty)4103		goto done;4104	if (context->op->is_empty(context))4105		goto done;4106 4107	for (r = 0; r >= 0 && tab && !tab->empty; r = restore_lexmin(tab)) {4108		int flags;4109		int row;4110		enum isl_tab_row_sign sgn;4111		int split = -1;4112		int n_split = 0;4113 4114		for (row = tab->n_redundant; row < tab->n_row; ++row) {4115			if (!isl_tab_var_from_row(tab, row)->is_nonneg)4116				continue;4117			sgn = row_sign(tab, sol, row);4118			if (!sgn)4119				goto error;4120			tab->row_sign[row] = sgn;4121			if (sgn == isl_tab_row_any)4122				n_split++;4123			if (sgn == isl_tab_row_any && split == -1)4124				split = row;4125			if (sgn == isl_tab_row_neg)4126				break;4127		}4128		if (row < tab->n_row)4129			continue;4130		if (split != -1) {4131			struct isl_vec *ineq;4132			if (n_split != 1)4133				split = context->op->best_split(context, tab);4134			if (split < 0)4135				goto error;4136			ineq = get_row_parameter_ineq(tab, split);4137			if (!ineq)4138				goto error;4139			is_strict(ineq);4140			reset_any_to_unknown(tab);4141			tab->row_sign[split] = isl_tab_row_pos;4142			sol_inc_level(sol);4143			find_in_pos(sol, tab, ineq->el);4144			tab->row_sign[split] = isl_tab_row_neg;4145			isl_seq_neg(ineq->el, ineq->el, ineq->size);4146			isl_int_sub_ui(ineq->el[0], ineq->el[0], 1);4147			if (!sol->error)4148				context->op->add_ineq(context, ineq->el, 0, 1);4149			isl_vec_free(ineq);4150			if (sol->error)4151				goto error;4152			continue;4153		}4154		if (tab->rational)4155			break;4156		row = first_non_integer_row(tab, &flags);4157		if (row < 0)4158			break;4159		if (ISL_FL_ISSET(flags, I_PAR)) {4160			if (ISL_FL_ISSET(flags, I_VAR)) {4161				if (isl_tab_mark_empty(tab) < 0)4162					goto error;4163				break;4164			}4165			row = add_cut(tab, row);4166		} else if (ISL_FL_ISSET(flags, I_VAR)) {4167			struct isl_vec *div;4168			struct isl_vec *ineq;4169			int d;4170			div = get_row_split_div(tab, row);4171			if (!div)4172				goto error;4173			d = context->op->get_div(context, tab, div);4174			isl_vec_free(div);4175			if (d < 0)4176				goto error;4177			ineq = ineq_for_div(context->op->peek_basic_set(context), d);4178			if (!ineq)4179				goto error;4180			sol_inc_level(sol);4181			no_sol_in_strict(sol, tab, ineq);4182			isl_seq_neg(ineq->el, ineq->el, ineq->size);4183			context->op->add_ineq(context, ineq->el, 1, 1);4184			isl_vec_free(ineq);4185			if (sol->error || !context->op->is_ok(context))4186				goto error;4187			tab = set_row_cst_to_div(tab, row, d);4188			if (context->op->is_empty(context))4189				break;4190		} else4191			row = add_parametric_cut(tab, row, context);4192		if (row < 0)4193			goto error;4194	}4195	if (r < 0)4196		goto error;4197done:4198	sol_add(sol, tab);4199	isl_tab_free(tab);4200	return;4201error:4202	isl_tab_free(tab);4203	sol->error = 1;4204}4205 4206/* Does "sol" contain a pair of partial solutions that could potentially4207 * be merged?4208 *4209 * We currently only check that "sol" is not in an error state4210 * and that there are at least two partial solutions of which the final two4211 * are defined at the same level.4212 */4213static int sol_has_mergeable_solutions(struct isl_sol *sol)4214{4215	if (sol->error)4216		return 0;4217	if (!sol->partial)4218		return 0;4219	if (!sol->partial->next)4220		return 0;4221	return sol->partial->level == sol->partial->next->level;4222}4223 4224/* Compute the lexicographic minimum of the set represented by the main4225 * tableau "tab" within the context "sol->context_tab".4226 *4227 * As a preprocessing step, we first transfer all the purely parametric4228 * equalities from the main tableau to the context tableau, i.e.,4229 * parameters that have been pivoted to a row.4230 * These equalities are ignored by the main algorithm, because the4231 * corresponding rows may not be marked as being non-negative.4232 * In parts of the context where the added equality does not hold,4233 * the main tableau is marked as being empty.4234 *4235 * Before we embark on the actual computation, we save a copy4236 * of the context.  When we return, we check if there are any4237 * partial solutions that can potentially be merged.  If so,4238 * we perform a rollback to the initial state of the context.4239 * The merging of partial solutions happens inside calls to4240 * sol_dec_level that are pushed onto the undo stack of the context.4241 * If there are no partial solutions that can potentially be merged4242 * then the rollback is skipped as it would just be wasted effort.4243 */4244static void find_solutions_main(struct isl_sol *sol, struct isl_tab *tab)4245{4246	int row;4247	void *saved;4248 4249	if (!tab)4250		goto error;4251 4252	sol->level = 0;4253 4254	for (row = tab->n_redundant; row < tab->n_row; ++row) {4255		int p;4256		struct isl_vec *eq;4257 4258		if (!row_is_parameter_var(tab, row))4259			continue;4260		if (tab->row_var[row] < tab->n_param)4261			p = tab->row_var[row];4262		else4263			p = tab->row_var[row]4264				+ tab->n_param - (tab->n_var - tab->n_div);4265 4266		eq = isl_vec_alloc(tab->mat->ctx, 1+tab->n_param+tab->n_div);4267		if (!eq)4268			goto error;4269		get_row_parameter_line(tab, row, eq->el);4270		isl_int_neg(eq->el[1 + p], tab->mat->row[row][0]);4271		eq = isl_vec_normalize(eq);4272 4273		sol_inc_level(sol);4274		no_sol_in_strict(sol, tab, eq);4275 4276		isl_seq_neg(eq->el, eq->el, eq->size);4277		sol_inc_level(sol);4278		no_sol_in_strict(sol, tab, eq);4279		isl_seq_neg(eq->el, eq->el, eq->size);4280 4281		sol->context->op->add_eq(sol->context, eq->el, 1, 1);4282 4283		isl_vec_free(eq);4284 4285		if (isl_tab_mark_redundant(tab, row) < 0)4286			goto error;4287 4288		if (sol->context->op->is_empty(sol->context))4289			break;4290 4291		row = tab->n_redundant - 1;4292	}4293 4294	saved = sol->context->op->save(sol->context);4295 4296	find_solutions(sol, tab);4297 4298	if (sol_has_mergeable_solutions(sol))4299		sol->context->op->restore(sol->context, saved);4300	else4301		sol->context->op->discard(saved);4302 4303	sol->level = 0;4304	sol_pop(sol);4305 4306	return;4307error:4308	isl_tab_free(tab);4309	sol->error = 1;4310}4311 4312/* Check if integer division "div" of "dom" also occurs in "bmap".4313 * If so, return its position within the divs.4314 * Otherwise, return a position beyond the integer divisions.4315 */4316static int find_context_div(__isl_keep isl_basic_map *bmap,4317	__isl_keep isl_basic_set *dom, unsigned div)4318{4319	int i;4320	isl_size b_v_div, d_v_div;4321	isl_size n_div;4322 4323	b_v_div = isl_basic_map_var_offset(bmap, isl_dim_div);4324	d_v_div = isl_basic_set_var_offset(dom, isl_dim_div);4325	n_div = isl_basic_map_dim(bmap, isl_dim_div);4326	if (b_v_div < 0 || d_v_div < 0 || n_div < 0)4327		return -1;4328 4329	if (isl_int_is_zero(dom->div[div][0]))4330		return n_div;4331	if (isl_seq_first_non_zero(dom->div[div] + 2 + d_v_div,4332				    dom->n_div) != -1)4333		return n_div;4334 4335	for (i = 0; i < n_div; ++i) {4336		if (isl_int_is_zero(bmap->div[i][0]))4337			continue;4338		if (isl_seq_first_non_zero(bmap->div[i] + 2 + d_v_div,4339					   (b_v_div - d_v_div) + n_div) != -1)4340			continue;4341		if (isl_seq_eq(bmap->div[i], dom->div[div], 2 + d_v_div))4342			return i;4343	}4344	return n_div;4345}4346 4347/* The correspondence between the variables in the main tableau,4348 * the context tableau, and the input map and domain is as follows.4349 * The first n_param and the last n_div variables of the main tableau4350 * form the variables of the context tableau.4351 * In the basic map, these n_param variables correspond to the4352 * parameters and the input dimensions.  In the domain, they correspond4353 * to the parameters and the set dimensions.4354 * The n_div variables correspond to the integer divisions in the domain.4355 * To ensure that everything lines up, we may need to copy some of the4356 * integer divisions of the domain to the map.  These have to be placed4357 * in the same order as those in the context and they have to be placed4358 * after any other integer divisions that the map may have.4359 * This function performs the required reordering.4360 */4361static __isl_give isl_basic_map *align_context_divs(4362	__isl_take isl_basic_map *bmap, __isl_keep isl_basic_set *dom)4363{4364	int i;4365	int common = 0;4366	int other;4367	unsigned bmap_n_div;4368 4369	bmap_n_div = isl_basic_map_dim(bmap, isl_dim_div);4370 4371	for (i = 0; i < dom->n_div; ++i) {4372		int pos;4373 4374		pos = find_context_div(bmap, dom, i);4375		if (pos < 0)4376			return isl_basic_map_free(bmap);4377		if (pos < bmap_n_div)4378			common++;4379	}4380	other = bmap_n_div - common;4381	if (dom->n_div - common > 0) {4382		bmap = isl_basic_map_cow(bmap);4383		bmap = isl_basic_map_extend(bmap, dom->n_div - common, 0, 0);4384		if (!bmap)4385			return NULL;4386	}4387	for (i = 0; i < dom->n_div; ++i) {4388		int pos = find_context_div(bmap, dom, i);4389		if (pos < 0)4390			bmap = isl_basic_map_free(bmap);4391		if (pos >= bmap_n_div) {4392			pos = isl_basic_map_alloc_div(bmap);4393			if (pos < 0)4394				goto error;4395			isl_int_set_si(bmap->div[pos][0], 0);4396			bmap_n_div++;4397		}4398		if (pos != other + i)4399			bmap = isl_basic_map_swap_div(bmap, pos, other + i);4400	}4401	return bmap;4402error:4403	isl_basic_map_free(bmap);4404	return NULL;4405}4406 4407/* Base case of isl_tab_basic_map_partial_lexopt, after removing4408 * some obvious symmetries.4409 *4410 * We make sure the divs in the domain are properly ordered,4411 * because they will be added one by one in the given order4412 * during the construction of the solution map.4413 * Furthermore, make sure that the known integer divisions4414 * appear before any unknown integer division because the solution4415 * may depend on the known integer divisions, while anything that4416 * depends on any variable starting from the first unknown integer4417 * division is ignored in sol_pma_add.4418 */4419static struct isl_sol *basic_map_partial_lexopt_base_sol(4420	__isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,4421	__isl_give isl_set **empty, int max,4422	struct isl_sol *(*init)(__isl_keep isl_basic_map *bmap,4423		    __isl_take isl_basic_set *dom, int track_empty, int max))4424{4425	struct isl_tab *tab;4426	struct isl_sol *sol = NULL;4427	struct isl_context *context;4428 4429	if (dom->n_div) {4430		dom = isl_basic_set_sort_divs(dom);4431		bmap = align_context_divs(bmap, dom);4432	}4433	sol = init(bmap, dom, !!empty, max);4434	if (!sol)4435		goto error;4436 4437	context = sol->context;4438	if (isl_basic_set_plain_is_empty(context->op->peek_basic_set(context)))4439		/* nothing */;4440	else if (isl_basic_map_plain_is_empty(bmap)) {4441		if (sol->add_empty)4442			sol->add_empty(sol,4443		    isl_basic_set_copy(context->op->peek_basic_set(context)));4444	} else {4445		tab = tab_for_lexmin(bmap,4446				    context->op->peek_basic_set(context), 1, max);4447		tab = context->op->detect_nonnegative_parameters(context, tab);4448		find_solutions_main(sol, tab);4449	}4450	if (sol->error)4451		goto error;4452 4453	isl_basic_map_free(bmap);4454	return sol;4455error:4456	sol_free(sol);4457	isl_basic_map_free(bmap);4458	return NULL;4459}4460 4461/* Base case of isl_tab_basic_map_partial_lexopt, after removing4462 * some obvious symmetries.4463 *4464 * We call basic_map_partial_lexopt_base_sol and extract the results.4465 */4466static __isl_give isl_map *basic_map_partial_lexopt_base(4467	__isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,4468	__isl_give isl_set **empty, int max)4469{4470	isl_map *result = NULL;4471	struct isl_sol *sol;4472	struct isl_sol_map *sol_map;4473 4474	sol = basic_map_partial_lexopt_base_sol(bmap, dom, empty, max,4475						&sol_map_init);4476	if (!sol)4477		return NULL;4478	sol_map = (struct isl_sol_map *) sol;4479 4480	result = isl_map_copy(sol_map->map);4481	if (empty)4482		*empty = isl_set_copy(sol_map->empty);4483	sol_free(&sol_map->sol);4484	return result;4485}4486 4487/* Return a count of the number of occurrences of the "n" first4488 * variables in the inequality constraints of "bmap".4489 */4490static __isl_give int *count_occurrences(__isl_keep isl_basic_map *bmap,4491	int n)4492{4493	int i, j;4494	isl_ctx *ctx;4495	int *occurrences;4496 4497	if (!bmap)4498		return NULL;4499	ctx = isl_basic_map_get_ctx(bmap);4500	occurrences = isl_calloc_array(ctx, int, n);4501	if (!occurrences)4502		return NULL;4503 4504	for (i = 0; i < bmap->n_ineq; ++i) {4505		for (j = 0; j < n; ++j) {4506			if (!isl_int_is_zero(bmap->ineq[i][1 + j]))4507				occurrences[j]++;4508		}4509	}4510 4511	return occurrences;4512}4513 4514/* Do all of the "n" variables with non-zero coefficients in "c"4515 * occur in exactly a single constraint.4516 * "occurrences" is an array of length "n" containing the number4517 * of occurrences of each of the variables in the inequality constraints.4518 */4519static int single_occurrence(int n, isl_int *c, int *occurrences)4520{4521	int i;4522 4523	for (i = 0; i < n; ++i) {4524		if (isl_int_is_zero(c[i]))4525			continue;4526		if (occurrences[i] != 1)4527			return 0;4528	}4529 4530	return 1;4531}4532 4533/* Do all of the "n" initial variables that occur in inequality constraint4534 * "ineq" of "bmap" only occur in that constraint?4535 */4536static int all_single_occurrence(__isl_keep isl_basic_map *bmap, int ineq,4537	int n)4538{4539	int i, j;4540 4541	for (i = 0; i < n; ++i) {4542		if (isl_int_is_zero(bmap->ineq[ineq][1 + i]))4543			continue;4544		for (j = 0; j < bmap->n_ineq; ++j) {4545			if (j == ineq)4546				continue;4547			if (!isl_int_is_zero(bmap->ineq[j][1 + i]))4548				return 0;4549		}4550	}4551 4552	return 1;4553}4554 4555/* Structure used during detection of parallel constraints.4556 * n_in: number of "input" variables: isl_dim_param + isl_dim_in4557 * n_out: number of "output" variables: isl_dim_out + isl_dim_div4558 * val: the coefficients of the output variables4559 */4560struct isl_constraint_equal_info {4561	unsigned n_in;4562	unsigned n_out;4563	isl_int *val;4564};4565 4566/* Check whether the coefficients of the output variables4567 * of the constraint in "entry" are equal to info->val.4568 */4569static isl_bool constraint_equal(const void *entry, const void *val)4570{4571	isl_int **row = (isl_int **)entry;4572	const struct isl_constraint_equal_info *info = val;4573	int eq;4574 4575	eq = isl_seq_eq((*row) + 1 + info->n_in, info->val, info->n_out);4576	return isl_bool_ok(eq);4577}4578 4579/* Check whether "bmap" has a pair of constraints that have4580 * the same coefficients for the output variables.4581 * Note that the coefficients of the existentially quantified4582 * variables need to be zero since the existentially quantified4583 * of the result are usually not the same as those of the input.4584 * Furthermore, check that each of the input variables that occur4585 * in those constraints does not occur in any other constraint.4586 * If so, return true and return the row indices of the two constraints4587 * in *first and *second.4588 */4589static isl_bool parallel_constraints(__isl_keep isl_basic_map *bmap,4590	int *first, int *second)4591{4592	int i;4593	isl_ctx *ctx;4594	int *occurrences = NULL;4595	struct isl_hash_table *table = NULL;4596	struct isl_hash_table_entry *entry;4597	struct isl_constraint_equal_info info;4598	isl_size nparam, n_in, n_out, n_div;4599 4600	ctx = isl_basic_map_get_ctx(bmap);4601	table = isl_hash_table_alloc(ctx, bmap->n_ineq);4602	if (!table)4603		goto error;4604 4605	nparam = isl_basic_map_dim(bmap, isl_dim_param);4606	n_in = isl_basic_map_dim(bmap, isl_dim_in);4607	n_out = isl_basic_map_dim(bmap, isl_dim_out);4608	n_div = isl_basic_map_dim(bmap, isl_dim_div);4609	if (nparam < 0 || n_in < 0 || n_out < 0 || n_div < 0)4610		goto error;4611	info.n_in = nparam + n_in;4612	occurrences = count_occurrences(bmap, info.n_in);4613	if (info.n_in && !occurrences)4614		goto error;4615	info.n_out = n_out + n_div;4616	for (i = 0; i < bmap->n_ineq; ++i) {4617		uint32_t hash;4618 4619		info.val = bmap->ineq[i] + 1 + info.n_in;4620		if (isl_seq_first_non_zero(info.val, n_out) < 0)4621			continue;4622		if (isl_seq_first_non_zero(info.val + n_out, n_div) >= 0)4623			continue;4624		if (!single_occurrence(info.n_in, bmap->ineq[i] + 1,4625					occurrences))4626			continue;4627		hash = isl_seq_get_hash(info.val, info.n_out);4628		entry = isl_hash_table_find(ctx, table, hash,4629					    constraint_equal, &info, 1);4630		if (!entry)4631			goto error;4632		if (entry->data)4633			break;4634		entry->data = &bmap->ineq[i];4635	}4636 4637	if (i < bmap->n_ineq) {4638		*first = ((isl_int **)entry->data) - bmap->ineq; 4639		*second = i;4640	}4641 4642	isl_hash_table_free(ctx, table);4643	free(occurrences);4644 4645	return isl_bool_ok(i < bmap->n_ineq);4646error:4647	isl_hash_table_free(ctx, table);4648	free(occurrences);4649	return isl_bool_error;4650}4651 4652/* Given a set of upper bounds in "var", add constraints to "bset"4653 * that make the i-th bound smallest.4654 *4655 * In particular, if there are n bounds b_i, then add the constraints4656 *4657 *	b_i <= b_j	for j > i4658 *	b_i <  b_j	for j < i4659 */4660static __isl_give isl_basic_set *select_minimum(__isl_take isl_basic_set *bset,4661	__isl_keep isl_mat *var, int i)4662{4663	isl_ctx *ctx;4664	int j, k;4665 4666	ctx = isl_mat_get_ctx(var);4667 4668	for (j = 0; j < var->n_row; ++j) {4669		if (j == i)4670			continue;4671		k = isl_basic_set_alloc_inequality(bset);4672		if (k < 0)4673			goto error;4674		isl_seq_combine(bset->ineq[k], ctx->one, var->row[j],4675				ctx->negone, var->row[i], var->n_col);4676		isl_int_set_si(bset->ineq[k][var->n_col], 0);4677		if (j < i)4678			isl_int_sub_ui(bset->ineq[k][0], bset->ineq[k][0], 1);4679	}4680 4681	bset = isl_basic_set_finalize(bset);4682 4683	return bset;4684error:4685	isl_basic_set_free(bset);4686	return NULL;4687}4688 4689/* Given a set of upper bounds on the last "input" variable m,4690 * construct a set that assigns the minimal upper bound to m, i.e.,4691 * construct a set that divides the space into cells where one4692 * of the upper bounds is smaller than all the others and assign4693 * this upper bound to m.4694 *4695 * In particular, if there are n bounds b_i, then the result4696 * consists of n basic sets, each one of the form4697 *4698 *	m = b_i4699 *	b_i <= b_j	for j > i4700 *	b_i <  b_j	for j < i4701 */4702static __isl_give isl_set *set_minimum(__isl_take isl_space *space,4703	__isl_take isl_mat *var)4704{4705	int i, k;4706	isl_basic_set *bset = NULL;4707	isl_set *set = NULL;4708 4709	if (!space || !var)4710		goto error;4711 4712	set = isl_set_alloc_space(isl_space_copy(space),4713				var->n_row, ISL_SET_DISJOINT);4714 4715	for (i = 0; i < var->n_row; ++i) {4716		bset = isl_basic_set_alloc_space(isl_space_copy(space), 0,4717					       1, var->n_row - 1);4718		k = isl_basic_set_alloc_equality(bset);4719		if (k < 0)4720			goto error;4721		isl_seq_cpy(bset->eq[k], var->row[i], var->n_col);4722		isl_int_set_si(bset->eq[k][var->n_col], -1);4723		bset = select_minimum(bset, var, i);4724		set = isl_set_add_basic_set(set, bset);4725	}4726 4727	isl_space_free(space);4728	isl_mat_free(var);4729	return set;4730error:4731	isl_basic_set_free(bset);4732	isl_set_free(set);4733	isl_space_free(space);4734	isl_mat_free(var);4735	return NULL;4736}4737 4738/* Given that the last input variable of "bmap" represents the minimum4739 * of the bounds in "cst", check whether we need to split the domain4740 * based on which bound attains the minimum.4741 *4742 * A split is needed when the minimum appears in an integer division4743 * or in an equality.  Otherwise, it is only needed if it appears in4744 * an upper bound that is different from the upper bounds on which it4745 * is defined.4746 */4747static isl_bool need_split_basic_map(__isl_keep isl_basic_map *bmap,4748	__isl_keep isl_mat *cst)4749{4750	int i, j;4751	isl_size total;4752	unsigned pos;4753 4754	pos = cst->n_col - 1;4755	total = isl_basic_map_dim(bmap, isl_dim_all);4756	if (total < 0)4757		return isl_bool_error;4758 4759	for (i = 0; i < bmap->n_div; ++i)4760		if (!isl_int_is_zero(bmap->div[i][2 + pos]))4761			return isl_bool_true;4762 4763	for (i = 0; i < bmap->n_eq; ++i)4764		if (!isl_int_is_zero(bmap->eq[i][1 + pos]))4765			return isl_bool_true;4766 4767	for (i = 0; i < bmap->n_ineq; ++i) {4768		if (isl_int_is_nonneg(bmap->ineq[i][1 + pos]))4769			continue;4770		if (!isl_int_is_negone(bmap->ineq[i][1 + pos]))4771			return isl_bool_true;4772		if (isl_seq_first_non_zero(bmap->ineq[i] + 1 + pos + 1,4773					   total - pos - 1) >= 0)4774			return isl_bool_true;4775 4776		for (j = 0; j < cst->n_row; ++j)4777			if (isl_seq_eq(bmap->ineq[i], cst->row[j], cst->n_col))4778				break;4779		if (j >= cst->n_row)4780			return isl_bool_true;4781	}4782 4783	return isl_bool_false;4784}4785 4786/* Given that the last set variable of "bset" represents the minimum4787 * of the bounds in "cst", check whether we need to split the domain4788 * based on which bound attains the minimum.4789 *4790 * We simply call need_split_basic_map here.  This is safe because4791 * the position of the minimum is computed from "cst" and not4792 * from "bmap".4793 */4794static isl_bool need_split_basic_set(__isl_keep isl_basic_set *bset,4795	__isl_keep isl_mat *cst)4796{4797	return need_split_basic_map(bset_to_bmap(bset), cst);4798}4799 4800/* Given that the last set variable of "set" represents the minimum4801 * of the bounds in "cst", check whether we need to split the domain4802 * based on which bound attains the minimum.4803 */4804static isl_bool need_split_set(__isl_keep isl_set *set, __isl_keep isl_mat *cst)4805{4806	int i;4807 4808	for (i = 0; i < set->n; ++i) {4809		isl_bool split;4810 4811		split = need_split_basic_set(set->p[i], cst);4812		if (split < 0 || split)4813			return split;4814	}4815 4816	return isl_bool_false;4817}4818 4819/* Given a map of which the last input variable is the minimum4820 * of the bounds in "cst", split each basic set in the set4821 * in pieces where one of the bounds is (strictly) smaller than the others.4822 * This subdivision is given in "min_expr".4823 * The variable is subsequently projected out.4824 *4825 * We only do the split when it is needed.4826 * For example if the last input variable m = min(a,b) and the only4827 * constraints in the given basic set are lower bounds on m,4828 * i.e., l <= m = min(a,b), then we can simply project out m4829 * to obtain l <= a and l <= b, without having to split on whether4830 * m is equal to a or b.4831 */4832static __isl_give isl_map *split_domain(__isl_take isl_map *opt,4833	__isl_take isl_set *min_expr, __isl_take isl_mat *cst)4834{4835	isl_size n_in;4836	int i;4837	isl_space *space;4838	isl_map *res;4839 4840	n_in = isl_map_dim(opt, isl_dim_in);4841	if (n_in < 0 || !min_expr || !cst)4842		goto error;4843 4844	space = isl_map_get_space(opt);4845	space = isl_space_drop_dims(space, isl_dim_in, n_in - 1, 1);4846	res = isl_map_empty(space);4847 4848	for (i = 0; i < opt->n; ++i) {4849		isl_map *map;4850		isl_bool split;4851 4852		map = isl_map_from_basic_map(isl_basic_map_copy(opt->p[i]));4853		split = need_split_basic_map(opt->p[i], cst);4854		if (split < 0)4855			map = isl_map_free(map);4856		else if (split)4857			map = isl_map_intersect_domain(map,4858						       isl_set_copy(min_expr));4859		map = isl_map_remove_dims(map, isl_dim_in, n_in - 1, 1);4860 4861		res = isl_map_union_disjoint(res, map);4862	}4863 4864	isl_map_free(opt);4865	isl_set_free(min_expr);4866	isl_mat_free(cst);4867	return res;4868error:4869	isl_map_free(opt);4870	isl_set_free(min_expr);4871	isl_mat_free(cst);4872	return NULL;4873}4874 4875/* Given a set of which the last set variable is the minimum4876 * of the bounds in "cst", split each basic set in the set4877 * in pieces where one of the bounds is (strictly) smaller than the others.4878 * This subdivision is given in "min_expr".4879 * The variable is subsequently projected out.4880 */4881static __isl_give isl_set *split(__isl_take isl_set *empty,4882	__isl_take isl_set *min_expr, __isl_take isl_mat *cst)4883{4884	isl_map *map;4885 4886	map = isl_map_from_domain(empty);4887	map = split_domain(map, min_expr, cst);4888	empty = isl_map_domain(map);4889 4890	return empty;4891}4892 4893static __isl_give isl_map *basic_map_partial_lexopt(4894	__isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,4895	__isl_give isl_set **empty, int max);4896 4897/* This function is called from basic_map_partial_lexopt_symm.4898 * The last variable of "bmap" and "dom" corresponds to the minimum4899 * of the bounds in "cst".  "map_space" is the space of the original4900 * input relation (of basic_map_partial_lexopt_symm) and "set_space"4901 * is the space of the original domain.4902 *4903 * We recursively call basic_map_partial_lexopt and then plug in4904 * the definition of the minimum in the result.4905 */4906static __isl_give isl_map *basic_map_partial_lexopt_symm_core(4907	__isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,4908	__isl_give isl_set **empty, int max, __isl_take isl_mat *cst,4909	__isl_take isl_space *map_space, __isl_take isl_space *set_space)4910{4911	isl_map *opt;4912	isl_set *min_expr;4913 4914	min_expr = set_minimum(isl_basic_set_get_space(dom), isl_mat_copy(cst));4915 4916	opt = basic_map_partial_lexopt(bmap, dom, empty, max);4917 4918	if (empty) {4919		*empty = split(*empty,4920			       isl_set_copy(min_expr), isl_mat_copy(cst));4921		*empty = isl_set_reset_space(*empty, set_space);4922	}4923 4924	opt = split_domain(opt, min_expr, cst);4925	opt = isl_map_reset_space(opt, map_space);4926 4927	return opt;4928}4929 4930/* Extract a domain from "bmap" for the purpose of computing4931 * a lexicographic optimum.4932 *4933 * This function is only called when the caller wants to compute a full4934 * lexicographic optimum, i.e., without specifying a domain.  In this case,4935 * the caller is not interested in the part of the domain space where4936 * there is no solution and the domain can be initialized to those constraints4937 * of "bmap" that only involve the parameters and the input dimensions.4938 * This relieves the parametric programming engine from detecting those4939 * inequalities and transferring them to the context.  More importantly,4940 * it ensures that those inequalities are transferred first and not4941 * intermixed with inequalities that actually split the domain.4942 *4943 * If the caller does not require the absence of existentially quantified4944 * variables in the result (i.e., if ISL_OPT_QE is not set in "flags"),4945 * then the actual domain of "bmap" can be used.  This ensures that4946 * the domain does not need to be split at all just to separate out4947 * pieces of the domain that do not have a solution from piece that do.4948 * This domain cannot be used in general because it may involve4949 * (unknown) existentially quantified variables which will then also4950 * appear in the solution.4951 */4952static __isl_give isl_basic_set *extract_domain(__isl_keep isl_basic_map *bmap,4953	unsigned flags)4954{4955	isl_size n_div;4956	isl_size n_out;4957 4958	n_div = isl_basic_map_dim(bmap, isl_dim_div);4959	n_out = isl_basic_map_dim(bmap, isl_dim_out);4960	if (n_div < 0 || n_out < 0)4961		return NULL;4962	bmap = isl_basic_map_copy(bmap);4963	if (ISL_FL_ISSET(flags, ISL_OPT_QE)) {4964		bmap = isl_basic_map_drop_constraints_involving_dims(bmap,4965							isl_dim_div, 0, n_div);4966		bmap = isl_basic_map_drop_constraints_involving_dims(bmap,4967							isl_dim_out, 0, n_out);4968	}4969	return isl_basic_map_domain(bmap);4970}4971 4972#undef TYPE4973#define TYPE	isl_map4974#undef SUFFIX4975#define SUFFIX4976#include "isl_tab_lexopt_templ.c"4977 4978/* Extract the subsequence of the sample value of "tab"4979 * starting at "pos" and of length "len".4980 */4981static __isl_give isl_vec *extract_sample_sequence(struct isl_tab *tab,4982	int pos, int len)4983{4984	int i;4985	isl_ctx *ctx;4986	isl_vec *v;4987 4988	ctx = isl_tab_get_ctx(tab);4989	v = isl_vec_alloc(ctx, len);4990	if (!v)4991		return NULL;4992	for (i = 0; i < len; ++i) {4993		if (!tab->var[pos + i].is_row) {4994			isl_int_set_si(v->el[i], 0);4995		} else {4996			int row;4997 4998			row = tab->var[pos + i].index;4999			isl_int_divexact(v->el[i], tab->mat->row[row][1],5000					tab->mat->row[row][0]);5001		}5002	}5003 5004	return v;5005}5006 5007/* Check if the sequence of variables starting at "pos"5008 * represents a trivial solution according to "trivial".5009 * That is, is the result of applying "trivial" to this sequence5010 * equal to the zero vector?5011 */5012static isl_bool region_is_trivial(struct isl_tab *tab, int pos,5013	__isl_keep isl_mat *trivial)5014{5015	isl_size n, len;5016	isl_vec *v;5017	isl_bool is_trivial;5018 5019	n = isl_mat_rows(trivial);5020	if (n < 0)5021		return isl_bool_error;5022 5023	if (n == 0)5024		return isl_bool_false;5025 5026	len = isl_mat_cols(trivial);5027	if (len < 0)5028		return isl_bool_error;5029	v = extract_sample_sequence(tab, pos, len);5030	v = isl_mat_vec_product(isl_mat_copy(trivial), v);5031	is_trivial = isl_vec_is_zero(v);5032	isl_vec_free(v);5033 5034	return is_trivial;5035}5036 5037/* Global internal data for isl_tab_basic_set_non_trivial_lexmin.5038 *5039 * "n_op" is the number of initial coordinates to optimize,5040 * as passed to isl_tab_basic_set_non_trivial_lexmin.5041 * "region" is the "n_region"-sized array of regions passed5042 * to isl_tab_basic_set_non_trivial_lexmin.5043 *5044 * "tab" is the tableau that corresponds to the ILP problem.5045 * "local" is an array of local data structure, one for each5046 * (potential) level of the backtracking procedure of5047 * isl_tab_basic_set_non_trivial_lexmin.5048 * "v" is a pre-allocated vector that can be used for adding5049 * constraints to the tableau.5050 *5051 * "sol" contains the best solution found so far.5052 * It is initialized to a vector of size zero.5053 */5054struct isl_lexmin_data {5055	int n_op;5056	int n_region;5057	struct isl_trivial_region *region;5058 5059	struct isl_tab *tab;5060	struct isl_local_region *local;5061	isl_vec *v;5062 5063	isl_vec *sol;5064};5065 5066/* Return the index of the first trivial region, "n_region" if all regions5067 * are non-trivial or -1 in case of error.5068 */5069static int first_trivial_region(struct isl_lexmin_data *data)5070{5071	int i;5072 5073	for (i = 0; i < data->n_region; ++i) {5074		isl_bool trivial;5075		trivial = region_is_trivial(data->tab, data->region[i].pos,5076					data->region[i].trivial);5077		if (trivial < 0)5078			return -1;5079		if (trivial)5080			return i;5081	}5082 5083	return data->n_region;5084}5085 5086/* Check if the solution is optimal, i.e., whether the first5087 * n_op entries are zero.5088 */5089static int is_optimal(__isl_keep isl_vec *sol, int n_op)5090{5091	int i;5092 5093	for (i = 0; i < n_op; ++i)5094		if (!isl_int_is_zero(sol->el[1 + i]))5095			return 0;5096	return 1;5097}5098 5099/* Add constraints to "tab" that ensure that any solution is significantly5100 * better than that represented by "sol".  That is, find the first5101 * relevant (within first n_op) non-zero coefficient and force it (along5102 * with all previous coefficients) to be zero.5103 * If the solution is already optimal (all relevant coefficients are zero),5104 * then just mark the table as empty.5105 * "n_zero" is the number of coefficients that have been forced zero5106 * by previous calls to this function at the same level.5107 * Return the updated number of forced zero coefficients or -1 on error.5108 *5109 * This function assumes that at least 2 * (n_op - n_zero) more rows and5110 * at least 2 * (n_op - n_zero) more elements in the constraint array5111 * are available in the tableau.5112 */5113static int force_better_solution(struct isl_tab *tab,5114	__isl_keep isl_vec *sol, int n_op, int n_zero)5115{5116	int i, n;5117	isl_ctx *ctx;5118	isl_vec *v = NULL;5119 5120	if (!sol)5121		return -1;5122 5123	for (i = n_zero; i < n_op; ++i)5124		if (!isl_int_is_zero(sol->el[1 + i]))5125			break;5126 5127	if (i == n_op) {5128		if (isl_tab_mark_empty(tab) < 0)5129			return -1;5130		return n_op;5131	}5132 5133	ctx = isl_vec_get_ctx(sol);5134	v = isl_vec_alloc(ctx, 1 + tab->n_var);5135	if (!v)5136		return -1;5137 5138	n = i + 1;5139	for (; i >= n_zero; --i) {5140		v = isl_vec_clr(v);5141		isl_int_set_si(v->el[1 + i], -1);5142		if (add_lexmin_eq(tab, v->el) < 0)5143			goto error;5144	}5145 5146	isl_vec_free(v);5147	return n;5148error:5149	isl_vec_free(v);5150	return -1;5151}5152 5153/* Fix triviality direction "dir" of the given region to zero.5154 *5155 * This function assumes that at least two more rows and at least5156 * two more elements in the constraint array are available in the tableau.5157 */5158static isl_stat fix_zero(struct isl_tab *tab, struct isl_trivial_region *region,5159	int dir, struct isl_lexmin_data *data)5160{5161	isl_size len;5162 5163	data->v = isl_vec_clr(data->v);5164	if (!data->v)5165		return isl_stat_error;5166	len = isl_mat_cols(region->trivial);5167	if (len < 0)5168		return isl_stat_error;5169	isl_seq_cpy(data->v->el + 1 + region->pos, region->trivial->row[dir],5170		    len);5171	if (add_lexmin_eq(tab, data->v->el) < 0)5172		return isl_stat_error;5173 5174	return isl_stat_ok;5175}5176 5177/* This function selects case "side" for non-triviality region "region",5178 * assuming all the equality constraints have been imposed already.5179 * In particular, the triviality direction side/2 is made positive5180 * if side is even and made negative if side is odd.5181 *5182 * This function assumes that at least one more row and at least5183 * one more element in the constraint array are available in the tableau.5184 */5185static struct isl_tab *pos_neg(struct isl_tab *tab,5186	struct isl_trivial_region *region,5187	int side, struct isl_lexmin_data *data)5188{5189	isl_size len;5190 5191	data->v = isl_vec_clr(data->v);5192	if (!data->v)5193		goto error;5194	isl_int_set_si(data->v->el[0], -1);5195	len = isl_mat_cols(region->trivial);5196	if (len < 0)5197		goto error;5198	if (side % 2 == 0)5199		isl_seq_cpy(data->v->el + 1 + region->pos,5200			    region->trivial->row[side / 2], len);5201	else5202		isl_seq_neg(data->v->el + 1 + region->pos,5203			    region->trivial->row[side / 2], len);5204	return add_lexmin_ineq(tab, data->v->el);5205error:5206	isl_tab_free(tab);5207	return NULL;5208}5209 5210/* Local data at each level of the backtracking procedure of5211 * isl_tab_basic_set_non_trivial_lexmin.5212 *5213 * "update" is set if a solution has been found in the current case5214 * of this level, such that a better solution needs to be enforced5215 * in the next case.5216 * "n_zero" is the number of initial coordinates that have already5217 * been forced to be zero at this level.5218 * "region" is the non-triviality region considered at this level.5219 * "side" is the index of the current case at this level.5220 * "n" is the number of triviality directions.5221 * "snap" is a snapshot of the tableau holding a state that needs5222 * to be satisfied by all subsequent cases.5223 */5224struct isl_local_region {5225	int update;5226	int n_zero;5227	int region;5228	int side;5229	int n;5230	struct isl_tab_undo *snap;5231};5232 5233/* Initialize the global data structure "data" used while solving5234 * the ILP problem "bset".5235 */5236static isl_stat init_lexmin_data(struct isl_lexmin_data *data,5237	__isl_keep isl_basic_set *bset)5238{5239	isl_ctx *ctx;5240 5241	ctx = isl_basic_set_get_ctx(bset);5242 5243	data->tab = tab_for_lexmin(bset, NULL, 0, 0);5244	if (!data->tab)5245		return isl_stat_error;5246 5247	data->v = isl_vec_alloc(ctx, 1 + data->tab->n_var);5248	if (!data->v)5249		return isl_stat_error;5250	data->local = isl_calloc_array(ctx, struct isl_local_region,5251					data->n_region);5252	if (data->n_region && !data->local)5253		return isl_stat_error;5254 5255	data->sol = isl_vec_alloc(ctx, 0);5256 5257	return isl_stat_ok;5258}5259 5260/* Mark all outer levels as requiring a better solution5261 * in the next cases.5262 */5263static void update_outer_levels(struct isl_lexmin_data *data, int level)5264{5265	int i;5266 5267	for (i = 0; i < level; ++i)5268		data->local[i].update = 1;5269}5270 5271/* Initialize "local" to refer to region "region" and5272 * to initiate processing at this level.5273 */5274static isl_stat init_local_region(struct isl_local_region *local, int region,5275	struct isl_lexmin_data *data)5276{5277	isl_size n = isl_mat_rows(data->region[region].trivial);5278 5279	if (n < 0)5280		return isl_stat_error;5281	local->n = n;5282	local->region = region;5283	local->side = 0;5284	local->update = 0;5285	local->n_zero = 0;5286 5287	return isl_stat_ok;5288}5289 5290/* What to do next after entering a level of the backtracking procedure.5291 *5292 * error: some error has occurred; abort5293 * done: an optimal solution has been found; stop search5294 * backtrack: backtrack to the previous level5295 * handle: add the constraints for the current level and5296 * 	move to the next level5297 */5298enum isl_next {5299	isl_next_error = -1,5300	isl_next_done,5301	isl_next_backtrack,5302	isl_next_handle,5303};5304 5305/* Have all cases of the current region been considered?5306 * If there are n directions, then there are 2n cases.5307 *5308 * The constraints in the current tableau are imposed5309 * in all subsequent cases.  This means that if the current5310 * tableau is empty, then none of those cases should be considered5311 * anymore and all cases have effectively been considered.5312 */5313static int finished_all_cases(struct isl_local_region *local,5314	struct isl_lexmin_data *data)5315{5316	if (data->tab->empty)5317		return 1;5318	return local->side >= 2 * local->n;5319}5320 5321/* Enter level "level" of the backtracking search and figure out5322 * what to do next.  "init" is set if the level was entered5323 * from a higher level and needs to be initialized.5324 * Otherwise, the level is entered as a result of backtracking and5325 * the tableau needs to be restored to a position that can5326 * be used for the next case at this level.5327 * The snapshot is assumed to have been saved in the previous case,5328 * before the constraints specific to that case were added.5329 *5330 * In the initialization case, the local region is initialized5331 * to point to the first violated region.5332 * If the constraints of all regions are satisfied by the current5333 * sample of the tableau, then tell the caller to continue looking5334 * for a better solution or to stop searching if an optimal solution5335 * has been found.5336 *5337 * If the tableau is empty or if all cases at the current level5338 * have been considered, then the caller needs to backtrack as well.5339 */5340static enum isl_next enter_level(int level, int init,5341	struct isl_lexmin_data *data)5342{5343	struct isl_local_region *local = &data->local[level];5344 5345	if (init) {5346		int r;5347 5348		data->tab = cut_to_integer_lexmin(data->tab, CUT_ONE);5349		if (!data->tab)5350			return isl_next_error;5351		if (data->tab->empty)5352			return isl_next_backtrack;5353		r = first_trivial_region(data);5354		if (r < 0)5355			return isl_next_error;5356		if (r == data->n_region) {5357			update_outer_levels(data, level);5358			isl_vec_free(data->sol);5359			data->sol = isl_tab_get_sample_value(data->tab);5360			if (!data->sol)5361				return isl_next_error;5362			if (is_optimal(data->sol, data->n_op))5363				return isl_next_done;5364			return isl_next_backtrack;5365		}5366		if (level >= data->n_region)5367			isl_die(isl_vec_get_ctx(data->v), isl_error_internal,5368				"nesting level too deep",5369				return isl_next_error);5370		if (init_local_region(local, r, data) < 0)5371			return isl_next_error;5372		if (isl_tab_extend_cons(data->tab,5373				    2 * local->n + 2 * data->n_op) < 0)5374			return isl_next_error;5375	} else {5376		if (isl_tab_rollback(data->tab, local->snap) < 0)5377			return isl_next_error;5378	}5379 5380	if (finished_all_cases(local, data))5381		return isl_next_backtrack;5382	return isl_next_handle;5383}5384 5385/* If a solution has been found in the previous case at this level5386 * (marked by local->update being set), then add constraints5387 * that enforce a better solution in the present and all following cases.5388 * The constraints only need to be imposed once because they are5389 * included in the snapshot (taken in pick_side) that will be used in5390 * subsequent cases.5391 */5392static isl_stat better_next_side(struct isl_local_region *local,5393	struct isl_lexmin_data *data)5394{5395	if (!local->update)5396		return isl_stat_ok;5397 5398	local->n_zero = force_better_solution(data->tab,5399				data->sol, data->n_op, local->n_zero);5400	if (local->n_zero < 0)5401		return isl_stat_error;5402 5403	local->update = 0;5404 5405	return isl_stat_ok;5406}5407 5408/* Add constraints to data->tab that select the current case (local->side)5409 * at the current level.5410 *5411 * If the linear combinations v should not be zero, then the cases are5412 *	v_0 >= 15413 *	v_0 <= -15414 *	v_0 = 0 and v_1 >= 15415 *	v_0 = 0 and v_1 <= -15416 *	v_0 = 0 and v_1 = 0 and v_2 >= 15417 *	v_0 = 0 and v_1 = 0 and v_2 <= -15418 *	...5419 * in this order.5420 *5421 * A snapshot is taken after the equality constraint (if any) has been added5422 * such that the next case can start off from this position.5423 * The rollback to this position is performed in enter_level.5424 */5425static isl_stat pick_side(struct isl_local_region *local,5426	struct isl_lexmin_data *data)5427{5428	struct isl_trivial_region *region;5429	int side, base;5430 5431	region = &data->region[local->region];5432	side = local->side;5433	base = 2 * (side/2);5434 5435	if (side == base && base >= 2 &&5436	    fix_zero(data->tab, region, base / 2 - 1, data) < 0)5437		return isl_stat_error;5438 5439	local->snap = isl_tab_snap(data->tab);5440	if (isl_tab_push_basis(data->tab) < 0)5441		return isl_stat_error;5442 5443	data->tab = pos_neg(data->tab, region, side, data);5444	if (!data->tab)5445		return isl_stat_error;5446	return isl_stat_ok;5447}5448 5449/* Free the memory associated to "data".5450 */5451static void clear_lexmin_data(struct isl_lexmin_data *data)5452{5453	free(data->local);5454	isl_vec_free(data->v);5455	isl_tab_free(data->tab);5456}5457 5458/* Return the lexicographically smallest non-trivial solution of the5459 * given ILP problem.5460 *5461 * All variables are assumed to be non-negative.5462 *5463 * n_op is the number of initial coordinates to optimize.5464 * That is, once a solution has been found, we will only continue looking5465 * for solutions that result in significantly better values for those5466 * initial coordinates.  That is, we only continue looking for solutions5467 * that increase the number of initial zeros in this sequence.5468 *5469 * A solution is non-trivial, if it is non-trivial on each of the5470 * specified regions.  Each region represents a sequence of5471 * triviality directions on a sequence of variables that starts5472 * at a given position.  A solution is non-trivial on such a region if5473 * at least one of the triviality directions is non-zero5474 * on that sequence of variables.5475 *5476 * Whenever a conflict is encountered, all constraints involved are5477 * reported to the caller through a call to "conflict".5478 *5479 * We perform a simple branch-and-bound backtracking search.5480 * Each level in the search represents an initially trivial region5481 * that is forced to be non-trivial.5482 * At each level we consider 2 * n cases, where n5483 * is the number of triviality directions.5484 * In terms of those n directions v_i, we consider the cases5485 *	v_0 >= 15486 *	v_0 <= -15487 *	v_0 = 0 and v_1 >= 15488 *	v_0 = 0 and v_1 <= -15489 *	v_0 = 0 and v_1 = 0 and v_2 >= 15490 *	v_0 = 0 and v_1 = 0 and v_2 <= -15491 *	...5492 * in this order.5493 */5494__isl_give isl_vec *isl_tab_basic_set_non_trivial_lexmin(5495	__isl_take isl_basic_set *bset, int n_op, int n_region,5496	struct isl_trivial_region *region,5497	int (*conflict)(int con, void *user), void *user)5498{5499	struct isl_lexmin_data data = { n_op, n_region, region };5500	int level, init;5501 5502	if (!bset)5503		return NULL;5504 5505	if (init_lexmin_data(&data, bset) < 0)5506		goto error;5507	data.tab->conflict = conflict;5508	data.tab->conflict_user = user;5509 5510	level = 0;5511	init = 1;5512 5513	while (level >= 0) {5514		enum isl_next next;5515		struct isl_local_region *local = &data.local[level];5516 5517		next = enter_level(level, init, &data);5518		if (next < 0)5519			goto error;5520		if (next == isl_next_done)5521			break;5522		if (next == isl_next_backtrack) {5523			level--;5524			init = 0;5525			continue;5526		}5527 5528		if (better_next_side(local, &data) < 0)5529			goto error;5530		if (pick_side(local, &data) < 0)5531			goto error;5532 5533		local->side++;5534		level++;5535		init = 1;5536	}5537 5538	clear_lexmin_data(&data);5539	isl_basic_set_free(bset);5540 5541	return data.sol;5542error:5543	clear_lexmin_data(&data);5544	isl_basic_set_free(bset);5545	isl_vec_free(data.sol);5546	return NULL;5547}5548 5549/* Wrapper for a tableau that is used for computing5550 * the lexicographically smallest rational point of a non-negative set.5551 * This point is represented by the sample value of "tab",5552 * unless "tab" is empty.5553 */5554struct isl_tab_lexmin {5555	isl_ctx *ctx;5556	struct isl_tab *tab;5557};5558 5559/* Free "tl" and return NULL.5560 */5561__isl_null isl_tab_lexmin *isl_tab_lexmin_free(__isl_take isl_tab_lexmin *tl)5562{5563	if (!tl)5564		return NULL;5565	isl_ctx_deref(tl->ctx);5566	isl_tab_free(tl->tab);5567	free(tl);5568 5569	return NULL;5570}5571 5572/* Construct an isl_tab_lexmin for computing5573 * the lexicographically smallest rational point in "bset",5574 * assuming that all variables are non-negative.5575 */5576__isl_give isl_tab_lexmin *isl_tab_lexmin_from_basic_set(5577	__isl_take isl_basic_set *bset)5578{5579	isl_ctx *ctx;5580	isl_tab_lexmin *tl;5581 5582	if (!bset)5583		return NULL;5584 5585	ctx = isl_basic_set_get_ctx(bset);5586	tl = isl_calloc_type(ctx, struct isl_tab_lexmin);5587	if (!tl)5588		goto error;5589	tl->ctx = ctx;5590	isl_ctx_ref(ctx);5591	tl->tab = tab_for_lexmin(bset, NULL, 0, 0);5592	isl_basic_set_free(bset);5593	if (!tl->tab)5594		return isl_tab_lexmin_free(tl);5595	return tl;5596error:5597	isl_basic_set_free(bset);5598	isl_tab_lexmin_free(tl);5599	return NULL;5600}5601 5602/* Return the dimension of the set represented by "tl".5603 */5604int isl_tab_lexmin_dim(__isl_keep isl_tab_lexmin *tl)5605{5606	return tl ? tl->tab->n_var : -1;5607}5608 5609/* Add the equality with coefficients "eq" to "tl", updating the optimal5610 * solution if needed.5611 * The equality is added as two opposite inequality constraints.5612 */5613__isl_give isl_tab_lexmin *isl_tab_lexmin_add_eq(__isl_take isl_tab_lexmin *tl,5614	isl_int *eq)5615{5616	unsigned n_var;5617 5618	if (!tl || !eq)5619		return isl_tab_lexmin_free(tl);5620 5621	if (isl_tab_extend_cons(tl->tab, 2) < 0)5622		return isl_tab_lexmin_free(tl);5623	n_var = tl->tab->n_var;5624	isl_seq_neg(eq, eq, 1 + n_var);5625	tl->tab = add_lexmin_ineq(tl->tab, eq);5626	isl_seq_neg(eq, eq, 1 + n_var);5627	tl->tab = add_lexmin_ineq(tl->tab, eq);5628 5629	if (!tl->tab)5630		return isl_tab_lexmin_free(tl);5631 5632	return tl;5633}5634 5635/* Add cuts to "tl" until the sample value reaches an integer value or5636 * until the result becomes empty.5637 */5638__isl_give isl_tab_lexmin *isl_tab_lexmin_cut_to_integer(5639	__isl_take isl_tab_lexmin *tl)5640{5641	if (!tl)5642		return NULL;5643	tl->tab = cut_to_integer_lexmin(tl->tab, CUT_ONE);5644	if (!tl->tab)5645		return isl_tab_lexmin_free(tl);5646	return tl;5647}5648 5649/* Return the lexicographically smallest rational point in the basic set5650 * from which "tl" was constructed.5651 * If the original input was empty, then return a zero-length vector.5652 */5653__isl_give isl_vec *isl_tab_lexmin_get_solution(__isl_keep isl_tab_lexmin *tl)5654{5655	if (!tl)5656		return NULL;5657	if (tl->tab->empty)5658		return isl_vec_alloc(tl->ctx, 0);5659	else5660		return isl_tab_get_sample_value(tl->tab);5661}5662 5663struct isl_sol_pma {5664	struct isl_sol	sol;5665	isl_pw_multi_aff *pma;5666	isl_set *empty;5667};5668 5669static void sol_pma_free(struct isl_sol *sol)5670{5671	struct isl_sol_pma *sol_pma = (struct isl_sol_pma *) sol;5672	isl_pw_multi_aff_free(sol_pma->pma);5673	isl_set_free(sol_pma->empty);5674}5675 5676/* This function is called for parts of the context where there is5677 * no solution, with "bset" corresponding to the context tableau.5678 * Simply add the basic set to the set "empty".5679 */5680static void sol_pma_add_empty(struct isl_sol_pma *sol,5681	__isl_take isl_basic_set *bset)5682{5683	if (!bset || !sol->empty)5684		goto error;5685 5686	sol->empty = isl_set_grow(sol->empty, 1);5687	bset = isl_basic_set_simplify(bset);5688	bset = isl_basic_set_finalize(bset);5689	sol->empty = isl_set_add_basic_set(sol->empty, bset);5690	if (!sol->empty)5691		sol->sol.error = 1;5692	return;5693error:5694	isl_basic_set_free(bset);5695	sol->sol.error = 1;5696}5697 5698/* Given a basic set "dom" that represents the context and a tuple of5699 * affine expressions "maff" defined over this domain, construct5700 * an isl_pw_multi_aff with a single cell corresponding to "dom" and5701 * the affine expressions in "maff".5702 */5703static void sol_pma_add(struct isl_sol_pma *sol,5704	__isl_take isl_basic_set *dom, __isl_take isl_multi_aff *maff)5705{5706	isl_pw_multi_aff *pma;5707 5708	dom = isl_basic_set_simplify(dom);5709	dom = isl_basic_set_finalize(dom);5710	pma = isl_pw_multi_aff_alloc(isl_set_from_basic_set(dom), maff);5711	sol->pma = isl_pw_multi_aff_add_disjoint(sol->pma, pma);5712	if (!sol->pma)5713		sol->sol.error = 1;5714}5715 5716static void sol_pma_add_empty_wrap(struct isl_sol *sol,5717	__isl_take isl_basic_set *bset)5718{5719	sol_pma_add_empty((struct isl_sol_pma *)sol, bset);5720}5721 5722static void sol_pma_add_wrap(struct isl_sol *sol,5723	__isl_take isl_basic_set *dom, __isl_take isl_multi_aff *ma)5724{5725	sol_pma_add((struct isl_sol_pma *)sol, dom, ma);5726}5727 5728/* Construct an isl_sol_pma structure for accumulating the solution.5729 * If track_empty is set, then we also keep track of the parts5730 * of the context where there is no solution.5731 * If max is set, then we are solving a maximization, rather than5732 * a minimization problem, which means that the variables in the5733 * tableau have value "M - x" rather than "M + x".5734 */5735static struct isl_sol *sol_pma_init(__isl_keep isl_basic_map *bmap,5736	__isl_take isl_basic_set *dom, int track_empty, int max)5737{5738	struct isl_sol_pma *sol_pma = NULL;5739	isl_space *space;5740 5741	if (!bmap)5742		goto error;5743 5744	sol_pma = isl_calloc_type(bmap->ctx, struct isl_sol_pma);5745	if (!sol_pma)5746		goto error;5747 5748	sol_pma->sol.free = &sol_pma_free;5749	if (sol_init(&sol_pma->sol, bmap, dom, max) < 0)5750		goto error;5751	sol_pma->sol.add = &sol_pma_add_wrap;5752	sol_pma->sol.add_empty = track_empty ? &sol_pma_add_empty_wrap : NULL;5753	space = isl_space_copy(sol_pma->sol.space);5754	sol_pma->pma = isl_pw_multi_aff_empty(space);5755	if (!sol_pma->pma)5756		goto error;5757 5758	if (track_empty) {5759		sol_pma->empty = isl_set_alloc_space(isl_basic_set_get_space(dom),5760							1, ISL_SET_DISJOINT);5761		if (!sol_pma->empty)5762			goto error;5763	}5764 5765	isl_basic_set_free(dom);5766	return &sol_pma->sol;5767error:5768	isl_basic_set_free(dom);5769	sol_free(&sol_pma->sol);5770	return NULL;5771}5772 5773/* Base case of isl_tab_basic_map_partial_lexopt, after removing5774 * some obvious symmetries.5775 *5776 * We call basic_map_partial_lexopt_base_sol and extract the results.5777 */5778static __isl_give isl_pw_multi_aff *basic_map_partial_lexopt_base_pw_multi_aff(5779	__isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,5780	__isl_give isl_set **empty, int max)5781{5782	isl_pw_multi_aff *result = NULL;5783	struct isl_sol *sol;5784	struct isl_sol_pma *sol_pma;5785 5786	sol = basic_map_partial_lexopt_base_sol(bmap, dom, empty, max,5787						&sol_pma_init);5788	if (!sol)5789		return NULL;5790	sol_pma = (struct isl_sol_pma *) sol;5791 5792	result = isl_pw_multi_aff_copy(sol_pma->pma);5793	if (empty)5794		*empty = isl_set_copy(sol_pma->empty);5795	sol_free(&sol_pma->sol);5796	return result;5797}5798 5799/* Given that the last input variable of "maff" represents the minimum5800 * of some bounds, check whether we need to plug in the expression5801 * of the minimum.5802 *5803 * In particular, check if the last input variable appears in any5804 * of the expressions in "maff".5805 */5806static isl_bool need_substitution(__isl_keep isl_multi_aff *maff)5807{5808	int i;5809	isl_size n_in;5810	unsigned pos;5811 5812	n_in = isl_multi_aff_dim(maff, isl_dim_in);5813	if (n_in < 0)5814		return isl_bool_error;5815	pos = n_in - 1;5816 5817	for (i = 0; i < maff->n; ++i) {5818		isl_bool involves;5819 5820		involves = isl_aff_involves_dims(maff->u.p[i],5821						isl_dim_in, pos, 1);5822		if (involves < 0 || involves)5823			return involves;5824	}5825 5826	return isl_bool_false;5827}5828 5829/* Given a set of upper bounds on the last "input" variable m,5830 * construct a piecewise affine expression that selects5831 * the minimal upper bound to m, i.e.,5832 * divide the space into cells where one5833 * of the upper bounds is smaller than all the others and select5834 * this upper bound on that cell.5835 *5836 * In particular, if there are n bounds b_i, then the result5837 * consists of n cell, each one of the form5838 *5839 *	b_i <= b_j	for j > i5840 *	b_i <  b_j	for j < i5841 *5842 * The affine expression on this cell is5843 *5844 *	b_i5845 */5846static __isl_give isl_pw_aff *set_minimum_pa(__isl_take isl_space *space,5847	__isl_take isl_mat *var)5848{5849	int i;5850	isl_aff *aff = NULL;5851	isl_basic_set *bset = NULL;5852	isl_pw_aff *paff = NULL;5853	isl_space *pw_space;5854	isl_local_space *ls = NULL;5855 5856	if (!space || !var)5857		goto error;5858 5859	ls = isl_local_space_from_space(isl_space_copy(space));5860	pw_space = isl_space_copy(space);5861	pw_space = isl_space_from_domain(pw_space);5862	pw_space = isl_space_add_dims(pw_space, isl_dim_out, 1);5863	paff = isl_pw_aff_alloc_size(pw_space, var->n_row);5864 5865	for (i = 0; i < var->n_row; ++i) {5866		isl_pw_aff *paff_i;5867 5868		aff = isl_aff_alloc(isl_local_space_copy(ls));5869		bset = isl_basic_set_alloc_space(isl_space_copy(space), 0,5870					       0, var->n_row - 1);5871		if (!aff || !bset)5872			goto error;5873		isl_int_set_si(aff->v->el[0], 1);5874		isl_seq_cpy(aff->v->el + 1, var->row[i], var->n_col);5875		isl_int_set_si(aff->v->el[1 + var->n_col], 0);5876		bset = select_minimum(bset, var, i);5877		paff_i = isl_pw_aff_alloc(isl_set_from_basic_set(bset), aff);5878		paff = isl_pw_aff_add_disjoint(paff, paff_i);5879	}5880 5881	isl_local_space_free(ls);5882	isl_space_free(space);5883	isl_mat_free(var);5884	return paff;5885error:5886	isl_aff_free(aff);5887	isl_basic_set_free(bset);5888	isl_pw_aff_free(paff);5889	isl_local_space_free(ls);5890	isl_space_free(space);5891	isl_mat_free(var);5892	return NULL;5893}5894 5895/* Given a piecewise multi-affine expression of which the last input variable5896 * is the minimum of the bounds in "cst", plug in the value of the minimum.5897 * This minimum expression is given in "min_expr_pa".5898 * The set "min_expr" contains the same information, but in the form of a set.5899 * The variable is subsequently projected out.5900 *5901 * The implementation is similar to those of "split" and "split_domain".5902 * If the variable appears in a given expression, then minimum expression5903 * is plugged in.  Otherwise, if the variable appears in the constraints5904 * and a split is required, then the domain is split.  Otherwise, no split5905 * is performed.5906 */5907static __isl_give isl_pw_multi_aff *split_domain_pma(5908	__isl_take isl_pw_multi_aff *opt, __isl_take isl_pw_aff *min_expr_pa,5909	__isl_take isl_set *min_expr, __isl_take isl_mat *cst)5910{5911	isl_size n_in;5912	int i;5913	isl_space *space;5914	isl_pw_multi_aff *res;5915 5916	if (!opt || !min_expr || !cst)5917		goto error;5918 5919	n_in = isl_pw_multi_aff_dim(opt, isl_dim_in);5920	if (n_in < 0)5921		goto error;5922	space = isl_pw_multi_aff_get_space(opt);5923	space = isl_space_drop_dims(space, isl_dim_in, n_in - 1, 1);5924	res = isl_pw_multi_aff_empty(space);5925 5926	for (i = 0; i < opt->n; ++i) {5927		isl_bool subs;5928		isl_pw_multi_aff *pma;5929 5930		pma = isl_pw_multi_aff_alloc(isl_set_copy(opt->p[i].set),5931					 isl_multi_aff_copy(opt->p[i].maff));5932		subs = need_substitution(opt->p[i].maff);5933		if (subs < 0) {5934			pma = isl_pw_multi_aff_free(pma);5935		} else if (subs) {5936			pma = isl_pw_multi_aff_substitute(pma,5937					n_in - 1, min_expr_pa);5938		} else {5939			isl_bool split;5940			split = need_split_set(opt->p[i].set, cst);5941			if (split < 0)5942				pma = isl_pw_multi_aff_free(pma);5943			else if (split)5944				pma = isl_pw_multi_aff_intersect_domain(pma,5945						       isl_set_copy(min_expr));5946		}5947		pma = isl_pw_multi_aff_project_out(pma,5948						    isl_dim_in, n_in - 1, 1);5949 5950		res = isl_pw_multi_aff_add_disjoint(res, pma);5951	}5952 5953	isl_pw_multi_aff_free(opt);5954	isl_pw_aff_free(min_expr_pa);5955	isl_set_free(min_expr);5956	isl_mat_free(cst);5957	return res;5958error:5959	isl_pw_multi_aff_free(opt);5960	isl_pw_aff_free(min_expr_pa);5961	isl_set_free(min_expr);5962	isl_mat_free(cst);5963	return NULL;5964}5965 5966static __isl_give isl_pw_multi_aff *basic_map_partial_lexopt_pw_multi_aff(5967	__isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,5968	__isl_give isl_set **empty, int max);5969 5970/* This function is called from basic_map_partial_lexopt_symm.5971 * The last variable of "bmap" and "dom" corresponds to the minimum5972 * of the bounds in "cst".  "map_space" is the space of the original5973 * input relation (of basic_map_partial_lexopt_symm) and "set_space"5974 * is the space of the original domain.5975 *5976 * We recursively call basic_map_partial_lexopt and then plug in5977 * the definition of the minimum in the result.5978 */5979static __isl_give isl_pw_multi_aff *5980basic_map_partial_lexopt_symm_core_pw_multi_aff(5981	__isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,5982	__isl_give isl_set **empty, int max, __isl_take isl_mat *cst,5983	__isl_take isl_space *map_space, __isl_take isl_space *set_space)5984{5985	isl_pw_multi_aff *opt;5986	isl_pw_aff *min_expr_pa;5987	isl_set *min_expr;5988 5989	min_expr = set_minimum(isl_basic_set_get_space(dom), isl_mat_copy(cst));5990	min_expr_pa = set_minimum_pa(isl_basic_set_get_space(dom),5991					isl_mat_copy(cst));5992 5993	opt = basic_map_partial_lexopt_pw_multi_aff(bmap, dom, empty, max);5994 5995	if (empty) {5996		*empty = split(*empty,5997			       isl_set_copy(min_expr), isl_mat_copy(cst));5998		*empty = isl_set_reset_space(*empty, set_space);5999	}6000 6001	opt = split_domain_pma(opt, min_expr_pa, min_expr, cst);6002	opt = isl_pw_multi_aff_reset_space(opt, map_space);6003 6004	return opt;6005}6006 6007#undef TYPE6008#define TYPE	isl_pw_multi_aff6009#undef SUFFIX6010#define SUFFIX	_pw_multi_aff6011#include "isl_tab_lexopt_templ.c"6012