1	/*
2	* Copyright 2008-2009 Katholieke Universiteit Leuven
3	* Copyright 2010 INRIA Saclay
4	* Copyright 2016-2017 Sven Verdoolaege
5	*
6	* Use of this software is governed by the MIT license
7	*
8	* Written by Sven Verdoolaege, K.U.Leuven, Departement
9	* Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
10	* 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 file
30	* was inspired by the paper "Parametric Integer Programming" and the
31	* report "Solving systems of affine (in)equalities" by Paul Feautrier
32	* (and others).
33	*
34	* The strategy used for obtaining a feasible solution is different
35	* 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 pivot
37	* in a row that increases the constant term of the row holding the
38	* constraint, making sure the sample solution remains feasible
39	* 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 smallest
42	* increment to the sample solution.
43	*
44	* By starting out from a sample value that is lexicographically
45	* smaller than any integer point in the problem space, the first
46	* feasible integer sample point we find will also be the lexicographically
47	* 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 apply
50	* the "big parameter" trick. Any variable x is then not directly
51	* used in the tableau, but instead it is represented by another
52	* 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 and
58	* 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 parameters
61	* happen to be non-negative. If so, we drop the column corresponding
62	* to M from the initial context tableau.
63	* If isl_context_gbr is used, then the big parameter trick is only
64	* used in the main tableau.
65	*/
66
67	struct isl_context;
68	struct 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 divisions
122	* in the input domain.
123	*/
124	struct isl_context {
125	struct isl_context_op *op;
126	int n_unknown;
127	};
128
129	struct 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	*/
141	struct 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
149	struct isl_sol;
150	struct 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 to
156	* a parametric integer linear programming problem.
157	* Every time the algorithm reaches a state where a solution
158	* can be read off from the tableau, the function "add" is called
159	* on the isl_sol passed to find_solutions_main. In a state where
160	* 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 invalidates
164	* 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 each
167	* split in the context.
168	* If "max" is set, then a maximization problem is being solved, rather than
169	* a minimization problem, which means that the variables in the
170	* 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_map
178	* and (optionally) the parts of the context where there is no solution
179	* in an isl_set, and
180	* isl_sol_pma, which collects an isl_pw_multi_aff instead.
181	*/
182	struct 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
198	static 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(bset: partial->dom);
206	isl_multi_aff_free(multi: partial->ma);
207	free(ptr: partial);
208	}
209	isl_space_free(space: sol->space);
210	if (sol->context)
211	sol->context->op->free(sol->context);
212	sol->free(sol);
213	free(ptr: 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 domain
219	* with no solution.
220	*/
221	static 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;
241	error:
242	isl_basic_set_free(bset: dom);
243	isl_multi_aff_free(multi: ma);
244	sol->error = 1;
245	}
246
247	/* Check that the final columns of "M", starting at "first", are zero.
248	*/
249	static 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(mat: M);
257	cols = isl_mat_cols(mat: 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(p: M->row[i] + first, len: 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", which
270	* are defined over the local space "ls".
271	* The matrix "M" may have extra (zero) columns beyond the number
272	* of variables in "ls".
273	*/
274	static __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, type: isl_dim_all);
283	if (!ma || dim < 0 || !M)
284	goto error;
285
286	if (check_final_columns_are_zero(M, first: 1 + dim) < 0)
287	goto error;
288	for (i = 1; i < M->n_row; ++i) {
289	aff = isl_aff_alloc(ls: isl_local_space_copy(ls));
290	if (aff) {
291	isl_int_set(aff->v->el[0], M->row[0][0]);
292	isl_seq_cpy(dst: aff->v->el + 1, src: M->row[i], len: 1 + dim);
293	}
294	aff = isl_aff_normalize(aff);
295	ma = isl_multi_aff_set_aff(multi: ma, pos: i - 1, el: aff);
296	}
297	isl_local_space_free(ls);
298	isl_mat_free(mat: M);
299
300	return ma;
301	error:
302	isl_local_space_free(ls);
303	isl_mat_free(mat: M);
304	isl_multi_aff_free(multi: ma);
305	return NULL;
306	}
307
308	/* Push a partial solution represented by a domain and mapping M
309	* onto the stack of partial solutions.
310	*
311	* The affine matrix "M" maps the dimensions of the context
312	* to the output variables. Convert it into an isl_multi_aff and
313	* then call sol_push_sol.
314	*
315	* Note that the description of the initial context may have involved
316	* existentially quantified variables, in which case they also appear
317	* in "dom". These need to be removed before creating the affine
318	* expression because an affine expression cannot be defined in terms
319	* of existentially quantified variables without a known representation.
320	* Since newly added integer divisions are inserted before these
321	* existentially quantified variables, they are still in the final
322	* positions and the corresponding final columns of "M" are zero
323	* because align_context_divs adds the existentially quantified
324	* variables of the context to the main tableau without any constraints and
325	* any equality constraints that are added later on can only serve
326	* to eliminate these existentially quantified variables.
327	*/
328	static 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(bset: dom, type: 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(space: isl_space_copy(space: sol->space));
342	ls = isl_basic_set_get_local_space(bset: dom);
343	ls = isl_local_space_drop_dims(ls, type: isl_dim_div,
344	first: n_known, n: n_div - n_known);
345	ma = set_from_affine_matrix(ma, ls, M);
346
347	if (!ma)
348	dom = isl_basic_set_free(bset: dom);
349	sol_push_sol(sol, dom, ma);
350	return;
351	error:
352	isl_basic_set_free(bset: dom);
353	isl_mat_free(mat: M);
354	sol_push_sol(sol, NULL, NULL);
355	}
356
357	/* Pop one partial solution from the partial solution stack and
358	* pass it on to sol->add or sol->add_empty.
359	*/
360	static 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	else
370	sol->add_empty(sol, partial->dom);
371	free(ptr: partial);
372	}
373
374	/* Return a fresh copy of the domain represented by the context tableau.
375	*/
376	static 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(bset: sol->context->op->peek_basic_set(sol->context));
384	bset = isl_basic_set_update_from_tab(bset,
385	tab: 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	*/
392	static 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(multi1: s1->ma, multi2: s2->ma);
401	}
402
403	/* Swap the initial two partial solutions in "sol".
404	*
405	* That is, go from
406	*
407	* sol->partial = p1; p1->next = p2; p2->next = p3
408	*
409	* to
410	*
411	* sol->partial = p2; p2->next = p1; p1->next = p3
412	*/
413	static 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" into
424	* a partial solution with the current context domain of "sol" and
425	* 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) are
429	* replaced by (D,M2), where D is the domain of "sol", which is assumed
430	* to be the union of D1 and D2, while M1 is assumed to be equal to M2
431	* (at least on D1).
432	*/
433	static 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(bset: 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(bset: partial->dom);
450	isl_multi_aff_free(multi: partial->ma);
451	free(ptr: 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 first
460	* as otherwise they would have to be eliminated twice, in a more complicated
461	* context.
462	*/
463	static 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(bset: isl_basic_set_copy(bset: dom));
471	pma1 = isl_pw_multi_aff_alloc(set: isl_set_copy(set),
472	maff: isl_multi_aff_copy(multi: ma1));
473	pma2 = isl_pw_multi_aff_alloc(set, maff: isl_multi_aff_copy(multi: ma2));
474	equal = isl_pw_multi_aff_is_equal(pma1, pma2);
475	isl_pw_multi_aff_free(pma: pma1);
476	isl_pw_multi_aff_free(pma: pma2);
477
478	return equal;
479	}
480
481	/* The initial two partial solutions of "sol" are known to be at
482	* 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 when
490	* restricted to one of the domains. If so, combine the two
491	* 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 can
494	* be replaced by (D1+D2,M2) and similarly when M1 = M2 on D2.
495	*/
496	static 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(s1: partial, s2: 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(ma1: partial->ma, ma2: partial->next->ma,
510	dom: 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(ma1: partial->ma, ma2: partial->next->ma,
516	dom: 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 onto
530	* the stack at levels that are deeper than the current level.
531	* If the two topmost elements on the stack have the same level
532	* and represent the same solution, then their domains are combined.
533	* This combined domain is the same as the current context domain
534	* 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 solutions
536	* at the current level are also popped off.
537	*/
538	static 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	} else
562	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)
571	error: sol->error = 1;
572	}
573
574	static 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
584	static 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(sol: 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 tableau
594	* to decrease the level again when it gets rolled back across
595	* the current state. That is, dec_level will be called with
596	* the context tableau in the same state as it is when inc_level
597	* is called.
598	*/
599	static 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, callback: &sol->dec_level.callback) < 0)
609	sol->error = 1;
610	}
611
612	static 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(dst: mat->row[i], src: mat->row[i], f: m, len: 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 input
627	* map (including divs that were copied from the context)
628	* + the number of extra divs constructed
629	* Of these, the first tab->n_param and the last tab->n_div variables
630	* correspond to the variables in the context, i.e.,
631	* tab->n_param + tab->n_div = context_tab->n_var
632	* tab->n_param is equal to the number of parameters and input
633	* dimensions in the input map
634	* tab->n_div is equal to the number of divs in the context
635	*
636	* If there is no solution, then call add_empty with a basic set
637	* 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 maps
641	* 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 any
644	* 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 the
648	* context because that would change the meaning of the matrix.
649	*
650	* To extract the value of the output variables, it should be noted
651	* that we always use a big parameter M in the main tableau and so
652	* the variable stored in this tableau is not an output variable x itself, but
653	* x' = M + x (in case of minimization)
654	* or
655	* 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 number
658	* (-M for minimization and M for maximization).
659	* We currently assume that the output dimensions in the original map
660	* are bounded, so this cannot occur.
661	* Similarly, when x' appears in a row, then the coefficient of M in that
662	* row is necessarily 1.
663	* If the row in the tableau represents
664	* d x' = c + d M + e(y)
665	* then, in case of minimization, the corresponding row in the matrix
666	* will be
667	* 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 be
670	* -a c - a e(y)
671	*/
672	static 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, dom: bset, NULL);
692	return;
693	}
694
695	off = 2 + tab->M;
696
697	mat = isl_mat_alloc(ctx: tab->mat->ctx, n_row: 1 + sol->n_out,
698	n_col: 1 + tab->n_param + tab->n_div);
699	if (!mat)
700	goto error;
701
702	isl_int_init(m);
703
704	isl_seq_clr(p: mat->row[0] + 1, len: 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(p: mat->row[1 + row], len: 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, n_row: 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(dst: mat->row[1 + row], src: mat->row[1 + row],
746	len: mat->n_col);
747	}
748
749	isl_int_clear(m);
750
751	sol_push_sol_mat(sol, dom: bset, M: mat);
752	return;
753	error2:
754	isl_int_clear(m);
755	error:
756	isl_basic_set_free(bset);
757	isl_mat_free(mat);
758	sol->error = 1;
759	}
760
761	struct isl_sol_map {
762	struct isl_sol sol;
763	struct isl_map *map;
764	struct isl_set *empty;
765	};
766
767	static 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(map: sol_map->map);
771	isl_set_free(set: sol_map->empty);
772	}
773
774	/* This function is called for parts of the context where there is
775	* no solution, with "bset" corresponding to the context tableau.
776	* Simply add the basic set to the set "empty".
777	*/
778	static 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(set: sol->empty, n: 1);
785	bset = isl_basic_set_simplify(bset);
786	bset = isl_basic_set_finalize(bset);
787	sol->empty = isl_set_add_basic_set(set: sol->empty, bset: isl_basic_set_copy(bset));
788	if (!sol->empty)
789	goto error;
790	isl_basic_set_free(bset);
791	return;
792	error:
793	isl_basic_set_free(bset);
794	sol->sol.error = 1;
795	}
796
797	static void sol_map_add_empty_wrap(struct isl_sol *sol,
798	struct isl_basic_set *bset)
799	{
800	sol_map_add_empty(sol: (struct isl_sol_map *)sol, bset);
801	}
802
803	/* Given a basic set "dom" that represents the context and a tuple of
804	* affine expressions "ma" defined over this domain, construct a basic map
805	* that expresses this function on the domain.
806	*/
807	static 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(maff: ma, rational: sol->sol.rational);
816	bmap = isl_basic_map_intersect_domain(bmap, bset: dom);
817	sol->map = isl_map_grow(map: sol->map, n: 1);
818	sol->map = isl_map_add_basic_map(map: sol->map, bmap);
819	if (!sol->map)
820	sol->sol.error = 1;
821	return;
822	error:
823	isl_basic_set_free(bset: dom);
824	isl_multi_aff_free(multi: ma);
825	sol->sol.error = 1;
826	}
827
828	static 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(sol: (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 that
837	* 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 it
840	* does, it is a different big parameter.
841	*/
842	static 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 parameter
870	* be the same as the values of the constants will only be the same
871	* if these coefficients are also the same.
872	*/
873	static 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 parameter
903	* is equal to zero.
904	*/
905	static 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(ctx: tab->mat->ctx, size: 1 + tab->n_param + tab->n_div);
910	if (!ineq)
911	return NULL;
912
913	get_row_parameter_line(tab, row, line: ineq->el);
914	if (ineq)
915	ineq = isl_vec_normalize(vec: ineq);
916
917	return ineq;
918	}
919
920	/* Normalize a div expression of the form
921	*
922	* [(g*f(x) + c)/(g * m)]
923	*
924	* with c the constant term and f(x) the remaining coefficients, to
925	*
926	* [(f(x) + [c/g])/m]
927	*/
928	static void normalize_div(__isl_keep isl_vec *div)
929	{
930	isl_ctx *ctx = isl_vec_get_ctx(vec: div);
931	int len = div->size - 2;
932
933	isl_seq_gcd(p: div->el + 2, len, gcd: &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(dst: div->el + 2, src: div->el + 2, f: ctx->normalize_gcd, len);
942	}
943
944	/* Return an integer division for use in a parametric cut based
945	* on the given row.
946	* In particular, let the parametric constant of the row be
947	*
948	* \sum_i a_i y_i
949	*
950	* where y_0 = 1, but none of the y_i corresponds to the big parameter M.
951	* The div returned is equal to
952	*
953	* floor(\sum_i {-a_i} y_i) = floor((\sum_i (-a_i mod d) y_i)/d)
954	*/
955	static 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(ctx: tab->mat->ctx, size: 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, line: div->el + 1);
965	isl_seq_neg(dst: div->el + 1, src: div->el + 1, len: div->size - 1);
966	normalize_div(div);
967	isl_seq_fdiv_r(dst: div->el + 1, src: div->el + 1, m: div->el[0], len: div->size - 1);
968
969	return div;
970	}
971
972	/* Return an integer division for use in transferring an integrality constraint
973	* to the context.
974	* In particular, let the parametric constant of the row be
975	*
976	* \sum_i a_i y_i
977	*
978	* where y_0 = 1, but none of the y_i corresponds to the big parameter M.
979	* The the returned div is equal to
980	*
981	* floor(\sum_i {a_i} y_i) = floor((\sum_i (a_i mod d) y_i)/d)
982	*/
983	static 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(ctx: tab->mat->ctx, size: 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, line: div->el + 1);
993	normalize_div(div);
994	isl_seq_fdiv_r(dst: div->el + 1, src: div->el + 1, m: div->el[0], len: div->size - 1);
995
996	return div;
997	}
998
999	/* Construct and return an inequality that expresses an upper bound
1000	* on the given div.
1001	* In particular, if the div is given by
1002	*
1003	* d = floor(e/m)
1004	*
1005	* then the inequality expresses
1006	*
1007	* m d <= e
1008	*/
1009	static __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, type: 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(ctx: bset->ctx, size: 1 + total);
1023	if (!ineq)
1024	return NULL;
1025
1026	isl_seq_cpy(dst: ineq->el, src: bset->div[div] + 1, len: 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 created
1032	* 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_i
1035	*
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 been
1039	* verified to be integral. We can therefore simply replace coefficient b
1040	* by floor(b). For the coefficients of the parameters we have
1041	* floor(a_i) = a_i - {a_i}, while for the other coefficients, we have
1042	* floor(b) = b.
1043	*/
1044	static struct isl_tab *set_row_cst_to_div(struct isl_tab *tab, int row, int div)
1045	{
1046	isl_seq_fdiv_q(dst: tab->mat->row[row] + 1, src: tab->mat->row[row] + 1,
1047	m: tab->mat->row[row][0], len: 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(dst: tab->mat->row[row] + 1,
1057	m1: tab->mat->ctx->one, src1: tab->mat->row[row] + 1,
1058	m2: tab->mat->ctx->one, src2: tab->mat->row[drow] + 1,
1059	len: 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;
1068	error:
1069	isl_tab_free(tab);
1070	return NULL;
1071	}
1072
1073	/* Check if the (parametric) constant of the given row is obviously
1074	* 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 and
1078	* all non-zero coefficients of parameters are negative and
1079	* belong to non-negative parameters.
1080	*/
1081	static 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 obviously
1123	* 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 and
1127	* all non-zero coefficients of parameters are positive and
1128	* belong to non-negative parameters.
1129	*/
1130	static 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 would
1172	* lead to the lexicographically smallest increment in the sample
1173	* solution when leaving the basis in favor of the row.
1174	* Pivoting with column c will increment the sample value by a non-negative
1175	* constant times a_{V,c}/a_{r,c}, with a_{V,c} the elements of column c
1176	* 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 c
1180	* 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 minimal
1185	* increment. Otherwise, it's c2.
1186	*/
1187	static 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 and
1233	* it needs to refer to either one of the first tab->n_param variables or
1234	* one of the last tab->n_div variables.
1235	*/
1236	static 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	*/
1249	static int col_is_parameter_var(struct isl_tab *tab, int col)
1250	{
1251	return is_parameter_var(tab, index: tab->col_var[col]);
1252	}
1253
1254	/* Does row "row" of "tab" refer to a variable in the context tableau?
1255	*/
1256	static int row_is_parameter_var(struct isl_tab *tab, int row)
1257	{
1258	return is_parameter_var(tab, index: tab->row_var[row]);
1259	}
1260
1261	/* Given a row in the tableau, find and return the column that would
1262	* result in the lexicographically smallest, but positive, increment
1263	* in the sample point.
1264	* If there is no such column, then return tab->n_col.
1265	* If anything goes wrong, return -1.
1266	*/
1267	static 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, col: 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	else
1288	col = lexmin_col_pair(tab, row, col1: col, col2: j, tmp);
1289	isl_assert(tab->mat->ctx, col >= 0, goto error);
1290	}
1291
1292	isl_int_clear(tmp);
1293	return col;
1294	error:
1295	isl_int_clear(tmp);
1296	return -1;
1297	}
1298
1299	/* Return the first known violated constraint, i.e., a non-negative
1300	* constraint that currently has an either obviously negative value
1301	* 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	*/
1306	static 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, i: 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, i: 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-positive
1337	* is satisfied. This function is not called from the current code
1338	* but is useful during debugging.
1339	*/
1340	static void check_lexpos(struct isl_tab *tab) __attribute__ ((unused));
1341	static 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	else
1356	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, format: "lexneg column %d (row %d)\n",
1364	col, row);
1365	}
1366	if (var >= tab->n_var - tab->n_div)
1367	fprintf(stderr, format: "zero column %d\n", col);
1368	}
1369	}
1370
1371	/* Report to the caller that the given constraint is part of an encountered
1372	* conflict.
1373	*/
1374	static 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 constraints
1380	* involved in the row to the caller. That is, the row itself
1381	* (if it represents a constraint) and all constraint columns with
1382	* non-zero (and therefore negative) coefficients.
1383	*/
1384	static 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, con: ~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, col: 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, con: ~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, we
1415	* look for a pivoting column that would result in the lexicographically
1416	* smallest increment in the sample point. If there is no such column
1417	* then the tableau is infeasible.
1418	*/
1419	static int restore_lexmin(struct isl_tab *tab) WARN_UNUSED;
1420	static 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 appropriate
1446	* pivoting column.
1447	* In particular, if there are any non-zero coefficients among
1448	* the non-parameter variables, then we take the last of these
1449	* variables. Eliminating this variable in terms of the other
1450	* variables and/or parameters does not influence the property
1451	* that all column in the initial tableau are lexicographically
1452	* positive. The row corresponding to the eliminated variable
1453	* will only have non-zero entries below the diagonal of the
1454	* initial tableau. That is, we transform
1455	*
1456	* I I
1457	* 1 into a
1458	* I I
1459	*
1460	* If there is no such non-parameter variable, then we are dealing with
1461	* pure parameter equality and we pick any parameter with coefficient 1 or -1
1462	* for elimination. This will ensure that the eliminated parameter
1463	* always has an integer value whenever all the other parameters are integral.
1464	* If there is no such parameter then we return -1.
1465	*/
1466	static 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 there
1494	* is no need to resolve any constraint violations.
1495	*
1496	* This function assumes that at least two more rows and at least
1497	* two more elements in the constraint array are available in the tableau.
1498	*/
1499	static 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, line: 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, row: r);
1512	if (i < 0) {
1513	tab->con[r].is_nonneg = 1;
1514	if (isl_tab_push_var(tab, type: isl_tab_undo_nonneg, var: &tab->con[r]) < 0)
1515	goto error;
1516	isl_seq_neg(dst: eq, src: eq, len: 1 + tab->n_var);
1517	r = isl_tab_add_row(tab, line: eq);
1518	if (r < 0)
1519	goto error;
1520	tab->con[r].is_nonneg = 1;
1521	if (isl_tab_push_var(tab, type: isl_tab_undo_nonneg, var: &tab->con[r]) < 0)
1522	goto error;
1523	} else {
1524	if (isl_tab_pivot(tab, row: r, col: i) < 0)
1525	goto error;
1526	if (isl_tab_kill_col(tab, col: i) < 0)
1527	goto error;
1528	tab->n_eq++;
1529	}
1530
1531	return tab;
1532	error:
1533	isl_tab_free(tab);
1534	return NULL;
1535	}
1536
1537	/* Check if the given row is a pure constant.
1538	*/
1539	static int is_constant(struct isl_tab *tab, int row)
1540	{
1541	unsigned off = 2 + tab->M;
1542
1543	return isl_seq_first_non_zero(p: tab->mat->row[row] + off + tab->n_dead,
1544	len: 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	*/
1550	static 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 eliminate
1573	* a column.
1574	*
1575	* This function assumes that at least two more rows and at least
1576	* two more elements in the constraint array are available in the tableau.
1577	*/
1578	static int add_lexmin_eq(struct isl_tab *tab, isl_int *eq) WARN_UNUSED;
1579	static 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, line: eq);
1589	if (r1 < 0)
1590	return -1;
1591	tab->con[r1].is_nonneg = 1;
1592	if (isl_tab_push_var(tab, type: isl_tab_undo_nonneg, var: &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(dst: eq, src: eq, len: 1 + tab->n_var);
1614
1615	r2 = isl_tab_add_row(tab, line: eq);
1616	if (r2 < 0)
1617	return -1;
1618	tab->con[r2].is_nonneg = 1;
1619	if (isl_tab_push_var(tab, type: isl_tab_undo_nonneg, var: &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, col: tab->con[r1].index) < 0)
1629	return -1;
1630	} else if (!tab->con[r2].is_row) {
1631	if (isl_tab_kill_col(tab, col: tab->con[r2].index) < 0)
1632	return -1;
1633	}
1634
1635	if (tab->bmap) {
1636	tab->bmap = isl_basic_map_add_ineq(bmap: tab->bmap, ineq: eq);
1637	if (isl_tab_push(tab, type: isl_tab_undo_bmap_ineq) < 0)
1638	return -1;
1639	isl_seq_neg(dst: eq, src: eq, len: 1 + tab->n_var);
1640	tab->bmap = isl_basic_map_add_ineq(bmap: tab->bmap, ineq: eq);
1641	isl_seq_neg(dst: eq, src: eq, len: 1 + tab->n_var);
1642	if (isl_tab_push(tab, type: 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 using
1652	* restore_lexmin.
1653	*
1654	* This function assumes that at least one more row and at least
1655	* one more element in the constraint array are available in the tableau.
1656	*/
1657	static 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(bmap: tab->bmap, ineq);
1665	if (isl_tab_push(tab, type: isl_tab_undo_bmap_ineq) < 0)
1666	goto error;
1667	if (!tab->bmap)
1668	goto error;
1669	}
1670	r = isl_tab_add_row(tab, line: ineq);
1671	if (r < 0)
1672	goto error;
1673	tab->con[r].is_nonneg = 1;
1674	if (isl_tab_push_var(tab, type: isl_tab_undo_nonneg, var: &tab->con[r]) < 0)
1675	goto error;
1676	if (isl_tab_row_is_redundant(tab, row: tab->con[r].index)) {
1677	if (isl_tab_mark_redundant(tab, row: 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, row: tab->con[r].index))
1686	if (isl_tab_mark_redundant(tab, row: tab->con[r].index) < 0)
1687	goto error;
1688	return tab;
1689	error:
1690	isl_tab_free(tab);
1691	return NULL;
1692	}
1693
1694	/* Check if the coefficients of the parameters are all integral.
1695	*/
1696	static 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	*/
1724	static 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, col: 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	*/
1741	static 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 << 0
1748	#define I_PAR 1 << 1
1749	#define I_VAR 1 << 2
1750
1751	/* Check for next (non-parameter) variable after "var" (first if var == -1)
1752	* that is non-integer and therefore requires a cut and return
1753	* 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 part
1757	* are all integral and if so, set the corresponding flag in *f.
1758	* If the constant and the parameter part are integral, then the
1759	* current sample value is integral and no cut is required
1760	* (irrespective of whether the variable part is integral).
1761	*/
1762	static 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 and
1787	* 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 part
1791	* are all integral and if so, set the corresponding flag in *f.
1792	* If the constant and the parameter part are integral, then the
1793	* current sample value is integral and no cut is required
1794	* (irrespective of whether the variable part is integral).
1795	*/
1796	static int first_non_integer_row(struct isl_tab *tab, int *f)
1797	{
1798	int var = next_non_integer_var(tab, var: -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 sample
1804	* value of the given row.
1805	*
1806	* If the row is given by
1807	*
1808	* m r = f + \sum_i a_i y_i
1809	*
1810	* then the cut is
1811	*
1812	* c = - {-f/m} + \sum_i {a_i/m} y_i >= 0
1813	*
1814	* The big parameter, if any, is ignored, since it is assumed to be big
1815	* enough to be divisible by any integer.
1816	* If the tableau is actually a parametric tableau, then this function
1817	* 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 it
1821	* exists, is set accordingly.
1822	*
1823	* Return the row of the cut or -1.
1824	*/
1825	static 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, n_new: 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, type: isl_tab_undo_nonneg, var: &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 1
1859	#define CUT_ONE 0
1860
1861	/* Given a non-parametric tableau, add cuts until an integer
1862	* sample point is obtained or until the tableau is determined
1863	* 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 these
1866	* non-integer values and then resolve the violated
1867	* cut constraints using restore_lexmin.
1868	* If one of the corresponding rows is equal to an integral
1869	* combination of variables/constraints plus a non-integral constant,
1870	* then there is no way to obtain an integer point and we return
1871	* a tableau that is marked empty.
1872	* The parameter cutting_strategy controls the strategy used when adding cuts
1873	* to remove non-integer points. CUT_ALL adds all possible cuts
1874	* before continuing the search. CUT_ONE adds only one cut at a time.
1875	*/
1876	static 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, var: -1, f: &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, f: &flags)) != -1);
1902	if (restore_lexmin(tab) < 0)
1903	goto error;
1904	if (tab->empty)
1905	break;
1906	}
1907	return tab;
1908	error:
1909	isl_tab_free(tab);
1910	return NULL;
1911	}
1912
1913	/* Check whether all the currently active samples also satisfy the inequality
1914	* "ineq" (treated as an equality if eq is set).
1915	* Remove those samples that do not.
1916	*/
1917	static 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(p1: ineq, p2: tab->samples->row[i],
1933	len: 1 + tab->n_var, prod: &v);
1934	sgn = isl_int_sgn(v);
1935	if (eq ? (sgn == 0) : (sgn >= 0))
1936	continue;
1937	tab = isl_tab_drop_sample(tab, s: i);
1938	if (!tab)
1939	break;
1940	}
1941	isl_int_clear(v);
1942
1943	return tab;
1944	error:
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, or
1951	* all values of the variables are equal to the big parameter plus
1952	* some constant. This constant is the actual sample value.
1953	*/
1954	static 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 sample
1978	* values satisfies the most recently added constraint.
1979	*/
1980	static 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;
2008	error:
2009	isl_tab_free(tab);
2010	return NULL;
2011	}
2012
2013	/* Check if any of the currently active sample values satisfies
2014	* the inequality "ineq" (an equality if eq is set).
2015	*/
2016	static 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(p1: ineq, p2: tab->samples->row[i],
2032	len: 1 + tab->n_var, prod: &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" and
2043	* return isl_bool_true if the div is obviously non-negative.
2044	*/
2045	static 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(mat: tab->samples,
2061	n_row: tab->n_sample, n_col: 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(p1: div->el + 1, p2: samples->row[i],
2067	len: div->size - 1, prod: &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(mat: tab->samples, dst_col: 1 + pos,
2072	src_col: 1 + tab->n_var - 1, n: 1);
2073	if (!tab->samples)
2074	return isl_bool_error;
2075
2076	return isl_bool_ok(b: nonneg);
2077	}
2078
2079	/* Add a div specified by "div" to both the main tableau and
2080	* the context tableau. In case of the main tableau, we only
2081	* need to add an extra div. In the context tableau, we also
2082	* 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 integer
2086	* divisions in the context to ensure that it does not get
2087	* equated to some linear combination involving unknown integer
2088	* divisions.
2089	*/
2090	static 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, n_new: 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;
2120	error:
2121	context->op->invalidate(context);
2122	return -1;
2123	}
2124
2125	/* Return the position of the integer division that is equal to div/denom
2126	* if there is one. Otherwise, return a position beyond the integer divisions.
2127	*/
2128	static 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(bmap: tab->bmap, type: isl_dim_all);
2132	isl_size n_div;
2133
2134	n_div = isl_basic_map_dim(bmap: tab->bmap, type: 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(p1: tab->bmap->div[i] + 1, p2: div, len: 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	*/
2150	static 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(bmap: context_tab->bmap, type: isl_dim_div);
2161	d = find_div(tab: context_tab, div: div->el + 1, denom: 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 value
2171	* of the given row.
2172	* Let a_i be the coefficients of the constant term and the parameters
2173	* and let b_i be the coefficients of the variables or constraints
2174	* 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 as
2178	*
2179	* c = \sum_i -{-a_i} y_i + \sum_i {b_i} x_i + q >= 0
2180	*
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 accomplished
2183	* 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 corresponding
2185	* 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 least
2189	* non-positive), so row_sign is set accordingly.
2190	*
2191	* Return the row of the cut or -1.
2192	*/
2193	static 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(vec: div);
2215	if (d < 0)
2216	return -1;
2217
2218	if (isl_tab_extend_cons(tab, n_new: 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(dst: r_row + 1, m1: gcd, src1: r_row + 1,
2265	m2: r_row[0], src2: tab->mat->row[d_row] + 1,
2266	len: 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, type: isl_tab_undo_nonneg, var: &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 computing
2289	* the lexicographic minimum (or maximum) of bmap.
2290	* If not NULL, then dom is the domain where the minimum
2291	* should be computed. In this case, we set up a parametric
2292	* 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 variable
2296	* x' = M - x, rather than x' = M + x. This in turn means that the coefficient
2297	* of the variables in all constraints are negated prior to adding them
2298	* to the tableau.
2299	*/
2300	static __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, type: isl_dim_all);
2310	if (total < 0)
2311	return NULL;
2312	tab = isl_tab_alloc(ctx: bmap->ctx, n_row: 2 * bmap->n_eq + bmap->n_ineq + 1,
2313	n_var: 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(bset: dom, type: 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(dst: bmap->eq[i] + o_var,
2345	src: bmap->eq[i] + o_var, len: n_var);
2346	tab = add_lexmin_valid_eq(tab, eq: bmap->eq[i]);
2347	if (max)
2348	isl_seq_neg(dst: bmap->eq[i] + o_var,
2349	src: bmap->eq[i] + o_var, len: 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(dst: bmap->ineq[i] + o_var,
2358	src: bmap->ineq[i] + o_var, len: n_var);
2359	tab = add_lexmin_ineq(tab, ineq: bmap->ineq[i]);
2360	if (max)
2361	isl_seq_neg(dst: bmap->ineq[i] + o_var,
2362	src: bmap->ineq[i] + o_var, len: n_var);
2363	if (!tab || tab->empty)
2364	return tab;
2365	}
2366	return tab;
2367	error:
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 involves
2376	* 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 corresponding
2382	* to the first row is redundant with respect to that of the second row
2383	* 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 about
2386	* anything can happen to the sign of the other row.)
2387	*
2388	* As a simple heuristic, we therefore select the row that makes the most
2389	* of the other rows redundant.
2390	*
2391	* Perhaps it would also be useful to look at the number of constraints
2392	* 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 initialized
2397	* to some arbitrary value (0) anyway. Without this redundant initialization
2398	* valgrind may warn about uninitialized memory accesses when isl
2399	* is compiled with some versions of gcc.
2400	*/
2401	static 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(tab: context_tab, n_new: 2) < 0)
2410	return -1;
2411
2412	snap = isl_tab_snap(tab: 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, i: split)->is_nonneg)
2421	continue;
2422	if (tab->row_sign[split] != isl_tab_row_any)
2423	continue;
2424
2425	if (is_parametric_constant(tab, row: split))
2426	return split;
2427
2428	ineq = get_row_parameter_ineq(tab, row: split);
2429	if (!ineq)
2430	return -1;
2431	ok = isl_tab_add_ineq(tab: context_tab, ineq: ineq->el) >= 0;
2432	isl_vec_free(vec: ineq);
2433	if (!ok)
2434	return -1;
2435
2436	snap2 = isl_tab_snap(tab: 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, i: 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(tab: context_tab, ineq: ineq->el) >= 0;
2452	isl_vec_free(vec: 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(tab: context_tab, var))
2458	r++;
2459	if (isl_tab_rollback(tab: context_tab, snap: 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(tab: context_tab, snap) < 0)
2467	return -1;
2468	}
2469
2470	return best;
2471	}
2472
2473	static 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(tab: clex->tab);
2480	}
2481
2482	static 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
2488	static 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(tab: clex->tab, n_new: 2) < 0)
2493	goto error;
2494	if (add_lexmin_eq(tab: clex->tab, eq) < 0)
2495	goto error;
2496	if (check) {
2497	int v = tab_has_valid_sample(tab: clex->tab, ineq: eq, eq: 1);
2498	if (v < 0)
2499	goto error;
2500	if (!v)
2501	clex->tab = check_integer_feasible(tab: clex->tab);
2502	}
2503	if (update)
2504	clex->tab = check_samples(tab: clex->tab, ineq: eq, eq: 1);
2505	return;
2506	error:
2507	isl_tab_free(tab: clex->tab);
2508	clex->tab = NULL;
2509	}
2510
2511	static 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(tab: clex->tab, n_new: 1) < 0)
2516	goto error;
2517	clex->tab = add_lexmin_ineq(tab: clex->tab, ineq);
2518	if (check) {
2519	int v = tab_has_valid_sample(tab: clex->tab, ineq, eq: 0);
2520	if (v < 0)
2521	goto error;
2522	if (!v)
2523	clex->tab = check_integer_feasible(tab: clex->tab);
2524	}
2525	if (update)
2526	clex->tab = check_samples(tab: clex->tab, ineq, eq: 0);
2527	return;
2528	error:
2529	isl_tab_free(tab: clex->tab);
2530	clex->tab = NULL;
2531	}
2532
2533	static 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, check: 0, update: 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 currently
2541	* active sample values. See row_sign for more information.
2542	*/
2543	static 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(p1: tab->samples->row[i], p2: ineq,
2558	len: 1 + tab->n_var, prod: &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
2580	static 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(tab: clex->tab, ineq, strict);
2585	}
2586
2587	/* Check whether "ineq" can be added to the tableau without rendering
2588	* it infeasible.
2589	*/
2590	static 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(tab: clex->tab, n_new: 1) < 0)
2600	return -1;
2601
2602	snap = isl_tab_snap(tab: clex->tab);
2603	if (isl_tab_push_basis(tab: clex->tab) < 0)
2604	return -1;
2605	clex->tab = add_lexmin_ineq(tab: clex->tab, ineq);
2606	clex->tab = check_integer_feasible(tab: clex->tab);
2607	if (!clex->tab)
2608	return -1;
2609	feasible = !clex->tab->empty;
2610	if (isl_tab_rollback(tab: clex->tab, snap) < 0)
2611	return -1;
2612
2613	return feasible;
2614	}
2615
2616	static 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" and
2623	* return isl_bool_true if the div is obviously non-negative.
2624	* context_tab_add_div will always return isl_bool_true, because all variables
2625	* in a isl_context_lex tableau are non-negative.
2626	* However, if we are using a big parameter in the context, then this only
2627	* reflects the non-negativity of the variable used to _encode_ the
2628	* div, i.e., div' = M + div, so we can't draw any conclusions.
2629	*/
2630	static 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(tab: clex->tab, pos, div,
2636	add_ineq: context_lex_add_ineq_wrap, user: 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
2644	static int context_lex_detect_equalities(struct isl_context *context,
2645	struct isl_tab *tab)
2646	{
2647	return 0;
2648	}
2649
2650	static 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(tab: clex->tab);
2658	if (isl_tab_push_basis(tab: clex->tab) < 0)
2659	return -1;
2660	r = best_split(tab, context_tab: clex->tab);
2661
2662	if (r >= 0 && isl_tab_rollback(tab: clex->tab, snap) < 0)
2663	return -1;
2664
2665	return r;
2666	}
2667
2668	static 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
2676	static 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(tab: clex->tab);
2682	if (isl_tab_push_basis(tab: clex->tab) < 0)
2683	return NULL;
2684	if (isl_tab_save_samples(tab: clex->tab) < 0)
2685	return NULL;
2686
2687	return snap;
2688	}
2689
2690	static 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(tab: clex->tab, snap: (struct isl_tab_undo *)save) < 0) {
2694	isl_tab_free(tab: clex->tab);
2695	clex->tab = NULL;
2696	}
2697	}
2698
2699	static void context_lex_discard(void *save)
2700	{
2701	}
2702
2703	static 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 can
2710	* only attain non-negative values. If so, mark the parameter as non-negative
2711	* in the main tableau. This allows for a more direct identification of some
2712	* cases of violated constraints.
2713	*/
2714	static 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(ctx: tab->mat->ctx, size: 1 + context_tab->n_var);
2727	if (!ineq)
2728	goto error;
2729
2730	if (isl_tab_extend_cons(tab: context_tab, n_new: 1) < 0)
2731	goto error;
2732
2733	snap = isl_tab_snap(tab: context_tab);
2734
2735	n = 0;
2736	isl_seq_clr(p: ineq->el, len: 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(tab: context_tab, ineq: 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(tab: 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(tab: 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(mat: context_tab->mat, col: 2, n: 1);
2757	context_tab->M = 0;
2758	}
2759
2760	isl_vec_free(vec: ineq);
2761	return tab;
2762	error:
2763	isl_vec_free(vec: ineq);
2764	isl_tab_free(tab);
2765	return NULL;
2766	}
2767
2768	static 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(tab: clex->tab);
2778	if (isl_tab_push_basis(tab: clex->tab) < 0)
2779	goto error;
2780
2781	tab = tab_detect_nonnegative_parameters(tab, context_tab: clex->tab);
2782
2783	if (isl_tab_rollback(tab: clex->tab, snap) < 0)
2784	goto error;
2785
2786	return tab;
2787	error:
2788	isl_tab_free(tab);
2789	return NULL;
2790	}
2791
2792	static void context_lex_invalidate(struct isl_context *context)
2793	{
2794	struct isl_context_lex *clex = (struct isl_context_lex *)context;
2795	isl_tab_free(tab: clex->tab);
2796	clex->tab = NULL;
2797	}
2798
2799	static __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(tab: clex->tab);
2804	free(ptr: clex);
2805
2806	return NULL;
2807	}
2808
2809	struct 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
2830	static 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(bmap: bset_to_bmap(bset), NULL, M: 1, max: 0);
2837	if (isl_tab_track_bset(tab, bset) < 0)
2838	goto error;
2839	tab = isl_tab_init_samples(tab);
2840	return tab;
2841	error:
2842	isl_tab_free(tab);
2843	return NULL;
2844	}
2845
2846	static 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(bset: isl_basic_set_copy(bset: dom));
2860	if (restore_lexmin(tab: clex->tab) < 0)
2861	goto error;
2862	clex->tab = check_integer_feasible(tab: clex->tab);
2863	if (!clex->tab)
2864	goto error;
2865
2866	return &clex->context;
2867	error:
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 the
2875	* context. Any rational point in "shifted" can therefore be rounded
2876	* up to an integer point in the context.
2877	* If the context is constrained by any equality, then "shifted" is not used
2878	* as it would be empty.
2879	*/
2880	struct 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
2887	static 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, context_tab: cgbr->tab);
2894	}
2895
2896	static 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(tab: cgbr->tab);
2903	}
2904
2905	static 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, which
2912	* contains the constraints of the original tableau shifted
2913	* by the sum of all negative coefficients. This ensures
2914	* that any rational point in the shifted tableau can
2915	* be rounded up to yield an integer point in the original tableau.
2916	*/
2917	static 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(tab: cgbr->tab);
2922	isl_size dim = isl_basic_set_dim(bset, type: isl_dim_all);
2923
2924	if (dim < 0)
2925	return;
2926	cst = isl_vec_alloc(ctx: cgbr->tab->mat->ctx, size: 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, track: 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(vec: cst);
2946	}
2947
2948	/* Check if the shifted tableau is non-empty, and if so
2949	* use the sample point to construct an integer point
2950	* of the context tableau.
2951	*/
2952	static 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(ctx: cgbr->tab->mat->ctx, size: 0);
2962
2963	sample = isl_tab_get_sample_value(tab: cgbr->shifted);
2964	sample = isl_vec_ceil(vec: sample);
2965
2966	return sample;
2967	}
2968
2969	static __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
2986	static 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
2993	static 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(tab: cgbr->tab))
2999	return isl_tab_get_sample_value(tab: 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(vec: sample);
3009	}
3010
3011	if (!cgbr->cone) {
3012	bset = isl_tab_peek_bset(tab: cgbr->tab);
3013	cgbr->cone = isl_tab_from_recession_cone(bset, parametric: 0);
3014	if (!cgbr->cone)
3015	return NULL;
3016	if (isl_tab_track_bset(tab: cgbr->cone,
3017	bset: isl_basic_set_copy(bset)) < 0)
3018	return NULL;
3019	}
3020	if (isl_tab_detect_implicit_equalities(tab: 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(mat: 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(tab: cgbr->tab);
3037
3038	sample = isl_tab_sample(tab: cgbr->tab);
3039
3040	if (!sample || isl_tab_rollback(tab: cgbr->tab, snap) < 0) {
3041	isl_vec_free(vec: sample);
3042	return NULL;
3043	}
3044
3045	return sample;
3046	}
3047
3048	cone = isl_basic_set_dup(bset: isl_tab_peek_bset(tab: cgbr->cone));
3049	cone = drop_constant_terms(bset: cone);
3050	cone = isl_basic_set_update_from_tab(bset: cone, tab: cgbr->cone);
3051	cone = isl_basic_set_underlying_set(bset: cone);
3052	cone = isl_basic_set_gauss(bset: cone, NULL);
3053
3054	bset = isl_basic_set_dup(bset: isl_tab_peek_bset(tab: cgbr->tab));
3055	bset = isl_basic_set_update_from_tab(bset, tab: 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
3062	static 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(vec: sample);
3078	if (isl_tab_mark_empty(tab: cgbr->tab) < 0)
3079	goto error;
3080	return;
3081	}
3082
3083	if (isl_tab_add_sample(tab: cgbr->tab, sample) < 0)
3084	goto error;
3085
3086	return;
3087	error:
3088	isl_tab_free(tab: cgbr->tab);
3089	cgbr->tab = NULL;
3090	}
3091
3092	static 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, n_new: 2) < 0)
3098	goto error;
3099
3100	if (isl_tab_add_eq(tab, eq) < 0)
3101	goto error;
3102
3103	return tab;
3104	error:
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 after
3111	* 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 to
3115	* cgbr->shifted since they would conflict with each other.
3116	* Instead, we directly mark cgbr->shifted empty.
3117	*/
3118	static 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(tab: cgbr->tab, eq);
3124
3125	if (cgbr->shifted && !cgbr->shifted->empty && use_shifted(cgbr)) {
3126	if (isl_tab_mark_empty(tab: 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(tab: cgbr->cone, n_new: 2) < 0)
3132	goto error;
3133	if (isl_tab_add_eq(tab: cgbr->cone, eq) < 0)
3134	goto error;
3135	}
3136
3137	if (check) {
3138	int v = tab_has_valid_sample(tab: cgbr->tab, ineq: eq, 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(tab: cgbr->tab, ineq: eq, eq: 1);
3146	return;
3147	error:
3148	isl_tab_free(tab: cgbr->tab);
3149	cgbr->tab = NULL;
3150	}
3151
3152	static 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(tab: cgbr->tab, n_new: 1) < 0)
3158	goto error;
3159
3160	if (isl_tab_add_ineq(tab: 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(bmap: cgbr->tab->bmap, type: isl_dim_all);
3167	if (dim < 0)
3168	goto error;
3169
3170	if (isl_tab_extend_cons(tab: cgbr->shifted, n_new: 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(tab: 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(tab: cgbr->cone, n_new: 1) < 0)
3191	goto error;
3192	if (isl_tab_add_ineq(tab: cgbr->cone, ineq) < 0)
3193	goto error;
3194	}
3195
3196	return;
3197	error:
3198	isl_tab_free(tab: cgbr->tab);
3199	cgbr->tab = NULL;
3200	}
3201
3202	static 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(tab: cgbr->tab, ineq, eq: 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(tab: cgbr->tab, ineq, eq: 0);
3220	return;
3221	error:
3222	isl_tab_free(tab: cgbr->tab);
3223	cgbr->tab = NULL;
3224	}
3225
3226	static 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, check: 0, update: 0);
3230	return context->op->is_ok(context) ? isl_stat_ok : isl_stat_error;
3231	}
3232
3233	static 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(tab: cgbr->tab, ineq, strict);
3238	}
3239
3240	/* Check whether "ineq" can be added to the tableau without rendering
3241	* it infeasible.
3242	*/
3243	static 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(tab: cgbr->tab, n_new: 1) < 0)
3255	return -1;
3256
3257	snap = isl_tab_snap(tab: cgbr->tab);
3258	if (cgbr->shifted)
3259	shifted_snap = isl_tab_snap(tab: cgbr->shifted);
3260	if (cgbr->cone)
3261	cone_snap = isl_tab_snap(tab: 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(tab: cgbr->tab, snap) < 0)
3268	return -1;
3269	if (shifted_snap) {
3270	if (isl_tab_rollback(tab: cgbr->shifted, snap: shifted_snap))
3271	return -1;
3272	} else if (cgbr->shifted) {
3273	isl_tab_free(tab: cgbr->shifted);
3274	cgbr->shifted = NULL;
3275	}
3276	if (cone_snap) {
3277	if (isl_tab_rollback(tab: cgbr->cone, snap: cone_snap))
3278	return -1;
3279	} else if (cgbr->cone) {
3280	isl_tab_free(tab: cgbr->cone);
3281	cgbr->cone = NULL;
3282	}
3283
3284	return feasible;
3285	}
3286
3287	/* Return the column of the last of the variables associated to
3288	* a column that has a non-zero coefficient.
3289	* This function is called in a context where only coefficients
3290	* of parameters or divs can be non-zero.
3291	*/
3292	static 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 context
3314	* 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 pairs
3317	* of inequalities starting at inequality "first".
3318	*
3319	* We tentatively add each of these equalities to the main tableau
3320	* and if this happens to result in a row with a final coefficient
3321	* that is one or negative one, we use it to kill a column
3322	* in the main tableau. Otherwise, we discard the tentatively
3323	* added row.
3324	* This tentative addition of equality constraints turns
3325	* on the undo facility of the tableau. Turn it off again
3326	* at the end, assuming it was turned off to begin with.
3327	*
3328	* Return 0 on success and -1 on failure.
3329	*/
3330	static 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(ctx: tab->mat->ctx, size: 1 + tab->n_var);
3341	if (!eq)
3342	goto error;
3343
3344	if (isl_tab_extend_cons(tab, n_new: (cgbr->tab->bmap->n_ineq - first)/2) < 0)
3345	goto error;
3346
3347	isl_seq_clr(p: eq->el + 1 + tab->n_param,
3348	len: 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(dst: eq->el, src: cgbr->tab->bmap->ineq[i], len: 1 + tab->n_param);
3356	isl_seq_cpy(dst: eq->el + 1 + tab->n_var - tab->n_div,
3357	src: cgbr->tab->bmap->ineq[i] + 1 + tab->n_param,
3358	len: tab->n_div);
3359
3360	r = isl_tab_add_row(tab, line: eq->el);
3361	if (r < 0)
3362	goto error;
3363	r = tab->con[r].index;
3364	j = last_non_zero_var_col(tab, p: 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, row: r, col: j) < 0)
3374	goto error;
3375	if (isl_tab_kill_col(tab, col: 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(vec: eq);
3385
3386	return 0;
3387	error:
3388	isl_vec_free(vec: eq);
3389	isl_tab_free(tab: cgbr->tab);
3390	cgbr->tab = NULL;
3391	return -1;
3392	}
3393
3394	static 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(tab: cgbr->tab);
3402	cgbr->cone = isl_tab_from_recession_cone(bset, parametric: 0);
3403	if (!cgbr->cone)
3404	goto error;
3405	if (isl_tab_track_bset(tab: cgbr->cone,
3406	bset: isl_basic_set_copy(bset)) < 0)
3407	goto error;
3408	}
3409	if (isl_tab_detect_implicit_equalities(tab: cgbr->cone) < 0)
3410	goto error;
3411
3412	n_ineq = cgbr->tab->bmap->n_ineq;
3413	cgbr->tab = isl_tab_detect_equalities(tab: cgbr->tab, tab_cone: cgbr->cone);
3414	if (!cgbr->tab)
3415	return -1;
3416	if (cgbr->tab->bmap->n_ineq > n_ineq &&
3417	propagate_equalities(cgbr, tab, first: n_ineq) < 0)
3418	return -1;
3419
3420	return 0;
3421	error:
3422	isl_tab_free(tab: cgbr->tab);
3423	cgbr->tab = NULL;
3424	return -1;
3425	}
3426
3427	static 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
3433	static 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(bmap: cgbr->cone->bmap, type: 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(tab: cgbr->cone, n_new: 3) < 0)
3447	return isl_bool_error;
3448	if (isl_tab_extend_vars(tab: cgbr->cone, n_new: 1) < 0)
3449	return isl_bool_error;
3450	if ((r = isl_tab_insert_var(tab: cgbr->cone, pos)) <0)
3451	return isl_bool_error;
3452
3453	cgbr->cone->bmap = isl_basic_map_insert_div(bmap: cgbr->cone->bmap,
3454	pos: r - o_div, div);
3455	if (!cgbr->cone->bmap)
3456	return isl_bool_error;
3457	if (isl_tab_push_var(tab: cgbr->cone, type: isl_tab_undo_bmap_div,
3458	var: &cgbr->cone->var[r]) < 0)
3459	return isl_bool_error;
3460	}
3461	return context_tab_insert_div(tab: cgbr->tab, pos, div,
3462	add_ineq: context_gbr_add_ineq_wrap, user: context);
3463	}
3464
3465	static 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(tab: cgbr->tab);
3473	r = best_split(tab, context_tab: cgbr->tab);
3474
3475	if (r >= 0 && isl_tab_rollback(tab: cgbr->tab, snap) < 0)
3476	return -1;
3477
3478	return r;
3479	}
3480
3481	static 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
3489	struct 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
3495	static 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(tab: cgbr->tab);
3508	if (isl_tab_save_samples(tab: cgbr->tab) < 0)
3509	goto error;
3510
3511	if (cgbr->shifted)
3512	snap->shifted_snap = isl_tab_snap(tab: cgbr->shifted);
3513	else
3514	snap->shifted_snap = NULL;
3515
3516	if (cgbr->cone)
3517	snap->cone_snap = isl_tab_snap(tab: cgbr->cone);
3518	else
3519	snap->cone_snap = NULL;
3520
3521	return snap;
3522	error:
3523	free(ptr: snap);
3524	return NULL;
3525	}
3526
3527	static 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(tab: cgbr->tab, snap: snap->tab_snap) < 0)
3534	goto error;
3535
3536	if (snap->shifted_snap) {
3537	if (isl_tab_rollback(tab: cgbr->shifted, snap: snap->shifted_snap) < 0)
3538	goto error;
3539	} else if (cgbr->shifted) {
3540	isl_tab_free(tab: cgbr->shifted);
3541	cgbr->shifted = NULL;
3542	}
3543
3544	if (snap->cone_snap) {
3545	if (isl_tab_rollback(tab: cgbr->cone, snap: snap->cone_snap) < 0)
3546	goto error;
3547	} else if (cgbr->cone) {
3548	isl_tab_free(tab: cgbr->cone);
3549	cgbr->cone = NULL;
3550	}
3551
3552	free(ptr: snap);
3553
3554	return;
3555	error:
3556	free(ptr: snap);
3557	isl_tab_free(tab: cgbr->tab);
3558	cgbr->tab = NULL;
3559	}
3560
3561	static void context_gbr_discard(void *save)
3562	{
3563	struct isl_gbr_tab_undo *snap = (struct isl_gbr_tab_undo *)save;
3564	free(ptr: snap);
3565	}
3566
3567	static 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
3573	static void context_gbr_invalidate(struct isl_context *context)
3574	{
3575	struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
3576	isl_tab_free(tab: cgbr->tab);
3577	cgbr->tab = NULL;
3578	}
3579
3580	static __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(tab: cgbr->tab);
3585	isl_tab_free(tab: cgbr->shifted);
3586	isl_tab_free(tab: cgbr->cone);
3587	free(ptr: cgbr);
3588
3589	return NULL;
3590	}
3591
3592	struct 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
3613	static 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(bset: dom, track: 1);
3629	cgbr->tab = isl_tab_init_samples(tab: cgbr->tab);
3630	if (!cgbr->tab)
3631	goto error;
3632	check_gbr_integer_feasible(cgbr);
3633
3634	return &cgbr->context;
3635	error:
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 by
3642	* isl_context_lex_alloc or isl_context_gbr_alloc.
3643	* The shared "n_unknown" field is initialized to the number
3644	* of final unknown integer divisions in "dom".
3645	*/
3646	static 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	else
3658	context = isl_context_gbr_alloc(dom);
3659
3660	if (!context)
3661	return NULL;
3662
3663	first = isl_basic_set_first_unknown_div(bset: dom);
3664	n_div = isl_basic_set_dim(bset: dom, type: 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 track
3673	* of the solution of an optimization problem on "bmap" over
3674	* the domain "dom".
3675	* If "max" is set, then a maximization problem is being solved, rather than
3676	* a minimization problem, which means that the variables in the
3677	* tableau have value "M - x" rather than "M + x".
3678	*/
3679	static 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, type: 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 parts
3698	* of the context where there is no solution.
3699	* If max is set, then we are solving a maximization, rather than
3700	* a minimization problem, which means that the variables in the
3701	* tableau have value "M - x" rather than "M + x".
3702	*/
3703	static 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: &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(space: sol_map->sol.space);
3722	sol_map->map = isl_map_alloc_space(space, n: 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(space: isl_basic_set_get_space(bset: dom),
3728	n: 1, ISL_SET_DISJOINT);
3729	if (!sol_map->empty)
3730	goto error;
3731	}
3732
3733	isl_basic_set_free(bset: dom);
3734	return &sol_map->sol;
3735	error:
3736	isl_basic_set_free(bset: dom);
3737	sol_free(sol: &sol_map->sol);
3738	return NULL;
3739	}
3740
3741	/* Check whether all coefficients of (non-parameter) variables
3742	* are non-positive, meaning that no pivots can be performed on the row.
3743	*/
3744	static 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, col: 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 with
3762	* equality. This happens if the gcd of the coefficients is not a divisor
3763	* of the constant term. If so, scale the constraint down by the gcd
3764	* of the coefficients.
3765	*/
3766	static 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(p: vec->el + 1, len: vec->size - 1, gcd: &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(dst: vec->el + 1, src: vec->el + 1, f: gcd, len: 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 of
3785	* isl_tab_row_pos: always non-negative; no pivot needed
3786	* isl_tab_row_neg: always non-positive; pivot
3787	* isl_tab_row_any: can be both positive and negative; split
3788	*
3789	* We first handle some simple cases
3790	* - the row sign may be known already
3791	* - the row may be obviously non-negative
3792	* - the parametric constant may be equal to that of another row
3793	* for which we know the sign. This sign will be either "pos" or
3794	* "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 each
3797	* of the currently active samples. Based on the signs of these values
3798	* we make an initial determination of the sign of the row.
3799	*
3800	* all zero -> unk(nown)
3801	* all non-negative -> pos
3802	* all non-positive -> neg
3803	* both negative and positive -> all
3804	*
3805	* If we end up with "all", we are done.
3806	* Otherwise, we perform a check for positive and/or negative
3807	* values as follows.
3808	*
3809	* samples neg unk pos
3810	* <0 ? Y N Y N
3811	* pos any pos
3812	* >0 ? Y N Y N
3813	* any neg any neg
3814	*
3815	* There is no special sign for "zero", because we can usually treat zero
3816	* as either non-negative or non-positive, whatever works out best.
3817	* However, if the row is "critical", meaning that pivoting is impossible
3818	* then we don't want to limp zero with the non-positive case, because
3819	* then we we would lose the solution for those values of the parameters
3820	* where the value of the row is zero. Instead, we treat 0 as non-negative
3821	* 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 not
3823	* be satisfied with equality by any integer values of the parameters.
3824	* In this case, we normalize the constraint, but then a value of zero
3825	* for the normalized constraint is actually a positive value for the
3826	* 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 -> pos
3830	* all negative -> neg
3831	* both negative and non-negative -> all
3832	*
3833	* samples neg pos
3834	* <0 ? Y N
3835	* any pos
3836	* >=0 ? Y N
3837	* any neg
3838	*/
3839	static 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, row1: 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(vec: 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(dst: ineq->el, src: ineq->el, len: 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	else
3882	res = (res == isl_tab_row_unknown) ? isl_tab_row_neg
3883	: isl_tab_row_any;
3884	if (res == isl_tab_row_neg) {
3885	isl_seq_neg(dst: ineq->el, src: ineq->el, len: 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(vec: ineq);
3904	return res;
3905	error:
3906	isl_vec_free(vec: ineq);
3907	return isl_tab_row_unknown;
3908	}
3909
3910	static void find_solutions(struct isl_sol *sol, struct isl_tab *tab);
3911
3912	/* Find solutions for values of the parameters that satisfy the given
3913	* inequality.
3914	*
3915	* We currently take a snapshot of the context tableau that is reset
3916	* when we return from this function, while we make a copy of the main
3917	* tableau, leaving the original main tableau untouched.
3918	* These are fairly arbitrary choices. Making a copy also of the context
3919	* 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 signs
3923	* and that we need to do this before saving the current basis
3924	* such that the basis has been restore before we restore the row signs.
3925	*/
3926	static 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	else
3945	sol->context->op->discard(saved);
3946	return;
3947	error:
3948	sol->error = 1;
3949	}
3950
3951	/* Record the absence of solutions for those values of the parameters
3952	* that do not satisfy the given inequality with equality.
3953	*/
3954	static 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;
3979	error:
3980	sol->error = 1;
3981	}
3982
3983	/* Reset all row variables that are marked to have a sign that may
3984	* be both positive and negative to have an unknown sign.
3985	*/
3986	static 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, i: 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 main
3999	* tableau "tab" within the context "sol->context_tab".
4000	* On entry the sample value of the main tableau is lexicographically
4001	* less than or equal to this lexicographic minimum.
4002	* Pivots are performed until a feasible point is found, which is then
4003	* necessarily equal to the minimum, or until the tableau is found to
4004	* be infeasible. Some pivots may need to be performed for only some
4005	* feasible values of the context tableau. If so, the context tableau
4006	* 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 no
4009	* "obvious" pivots need to be performed on it, where "obvious" means
4010	* that the given row can be seen to be negative without looking at
4011	* 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 property
4014	* hold prior to the first iteration of the loop, while restore_lexmin
4015	* is called before every other iteration.
4016	*
4017	* Inside the main loop, we first examine the signs of the rows of
4018	* the main tableau within the context of the context tableau.
4019	* If we find a row that is always non-positive for all values of
4020	* the parameters satisfying the context tableau and negative for at
4021	* least one value of the parameters, we perform the appropriate pivot
4022	* and start over. An exception is the case where no pivot can be
4023	* performed on the row. In this case, we require that the sign of
4024	* the row is negative for all values of the parameters (rather than just
4025	* non-positive). This special case is handled inside row_sign, which
4026	* will say that the row can have any sign if it determines that it can
4027	* attain both negative and zero values.
4028	*
4029	* If we can't find a row that always requires a pivot, but we can find
4030	* one or more rows that require a pivot for some values of the parameters
4031	* (i.e., the row can attain both positive and negative signs), then we split
4032	* the context tableau into two parts, one where we force the sign to be
4033	* 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 and
4036	* perform the required pivot.
4037	*
4038	* If no such rows can be found, all rows are non-negative and we have
4039	* found a (rational) feasible point. If we only wanted a rational point
4040	* then we are done.
4041	* Otherwise, we check if all values of the sample point of the tableau
4042	* are integral for the variables. If so, we have found the minimal
4043	* integral point and we are done.
4044	* If the sample point is not integral, then we need to make a distinction
4045	* based on whether the constant term is non-integral or the coefficients
4046	* of the parameters. Furthermore, in order to decide how to handle
4047	* the non-integrality, we also need to know whether the coefficients
4048	* of the other columns in the tableau are integral. This leads
4049	* to the following table. The first two rows do not correspond
4050	* to a non-integral sample point and are only mentioned for completeness.
4051	*
4052	* constant parameters other
4053	*
4054	* int int int |
4055	* int int rat | -> no problem
4056	*
4057	* rat int int -> fail
4058	*
4059	* rat int rat -> cut
4060	*
4061	* int rat rat |
4062	* rat rat rat | -> parametric cut
4063	*
4064	* int rat int |
4065	* rat rat int | -> split context
4066	*
4067	* If the parametric constant is completely integral, then there is nothing
4068	* to be done. If the constant term is non-integral, but all the other
4069	* coefficient are integral, then there is nothing that can be done
4070	* and the tableau has no integral solution.
4071	* If, on the other hand, one or more of the other columns have rational
4072	* coefficients, but the parameter coefficients are all integral, then
4073	* 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 we
4077	* can perform a parametric cut in the main tableau by adding a new
4078	* integer division in the context tableau.
4079	* If all other columns have integral coefficients, then we need to
4080	* enforce that the rational combination of parameters (c + \sum a_i y_i)/m
4081	* is always integral. We do this by introducing an integer division
4082	* q = floor((c + \sum a_i y_i)/m) and stipulating that its argument should
4083	* 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 as
4085	* c + \sum a_i y_i - m q >= 0
4086	* -c - \sum a_i y_i + m q + m - 1 >= 0
4087	* it is sufficient to add the inequality
4088	* -c - \sum a_i y_i + m q >= 0
4089	* In the part of the context where this inequality does not hold, the
4090	* main tableau is marked as being empty.
4091	*/
4092	static 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, i: 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, row: split);
4137	if (!ineq)
4138	goto error;
4139	is_strict(vec: 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: ineq->el);
4144	tab->row_sign[split] = isl_tab_row_neg;
4145	isl_seq_neg(dst: ineq->el, src: ineq->el, len: 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(vec: ineq);
4150	if (sol->error)
4151	goto error;
4152	continue;
4153	}
4154	if (tab->rational)
4155	break;
4156	row = first_non_integer_row(tab, f: &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(vec: div);
4175	if (d < 0)
4176	goto error;
4177	ineq = ineq_for_div(bset: context->op->peek_basic_set(context), div: d);
4178	if (!ineq)
4179	goto error;
4180	sol_inc_level(sol);
4181	no_sol_in_strict(sol, tab, ineq);
4182	isl_seq_neg(dst: ineq->el, src: ineq->el, len: ineq->size);
4183	context->op->add_ineq(context, ineq->el, 1, 1);
4184	isl_vec_free(vec: ineq);
4185	if (sol->error || !context->op->is_ok(context))
4186	goto error;
4187	tab = set_row_cst_to_div(tab, row, div: d);
4188	if (context->op->is_empty(context))
4189	break;
4190	} else
4191	row = add_parametric_cut(tab, row, context);
4192	if (row < 0)
4193	goto error;
4194	}
4195	if (r < 0)
4196	goto error;
4197	done:
4198	sol_add(sol, tab);
4199	isl_tab_free(tab);
4200	return;
4201	error:
4202	isl_tab_free(tab);
4203	sol->error = 1;
4204	}
4205
4206	/* Does "sol" contain a pair of partial solutions that could potentially
4207	* be merged?
4208	*
4209	* We currently only check that "sol" is not in an error state
4210	* and that there are at least two partial solutions of which the final two
4211	* are defined at the same level.
4212	*/
4213	static 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 main
4225	* tableau "tab" within the context "sol->context_tab".
4226	*
4227	* As a preprocessing step, we first transfer all the purely parametric
4228	* 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 the
4231	* 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 copy
4236	* of the context. When we return, we check if there are any
4237	* 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 to
4240	* sol_dec_level that are pushed onto the undo stack of the context.
4241	* If there are no partial solutions that can potentially be merged
4242	* then the rollback is skipped as it would just be wasted effort.
4243	*/
4244	static 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	else
4263	p = tab->row_var[row]
4264	+ tab->n_param - (tab->n_var - tab->n_div);
4265
4266	eq = isl_vec_alloc(ctx: tab->mat->ctx, size: 1+tab->n_param+tab->n_div);
4267	if (!eq)
4268	goto error;
4269	get_row_parameter_line(tab, row, line: eq->el);
4270	isl_int_neg(eq->el[1 + p], tab->mat->row[row][0]);
4271	eq = isl_vec_normalize(vec: eq);
4272
4273	sol_inc_level(sol);
4274	no_sol_in_strict(sol, tab, ineq: eq);
4275
4276	isl_seq_neg(dst: eq->el, src: eq->el, len: eq->size);
4277	sol_inc_level(sol);
4278	no_sol_in_strict(sol, tab, ineq: eq);
4279	isl_seq_neg(dst: eq->el, src: eq->el, len: eq->size);
4280
4281	sol->context->op->add_eq(sol->context, eq->el, 1, 1);
4282
4283	isl_vec_free(vec: 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	else
4301	sol->context->op->discard(saved);
4302
4303	sol->level = 0;
4304	sol_pop(sol);
4305
4306	return;
4307	error:
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	*/
4316	static 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, type: isl_dim_div);
4324	d_v_div = isl_basic_set_var_offset(bset: dom, type: isl_dim_div);
4325	n_div = isl_basic_map_dim(bmap, type: 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(p: dom->div[div] + 2 + d_v_div,
4332	len: 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(p: bmap->div[i] + 2 + d_v_div,
4339	len: (b_v_div - d_v_div) + n_div) != -1)
4340	continue;
4341	if (isl_seq_eq(p1: bmap->div[i], p2: dom->div[div], len: 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 tableau
4350	* form the variables of the context tableau.
4351	* In the basic map, these n_param variables correspond to the
4352	* parameters and the input dimensions. In the domain, they correspond
4353	* 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 the
4356	* integer divisions of the domain to the map. These have to be placed
4357	* in the same order as those in the context and they have to be placed
4358	* after any other integer divisions that the map may have.
4359	* This function performs the required reordering.
4360	*/
4361	static __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, type: isl_dim_div);
4370
4371	for (i = 0; i < dom->n_div; ++i) {
4372	int pos;
4373
4374	pos = find_context_div(bmap, dom, div: 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(base: bmap, extra: dom->n_div - common, n_eq: 0, n_ineq: 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, div: 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, a: pos, b: other + i);
4400	}
4401	return bmap;
4402	error:
4403	isl_basic_map_free(bmap);
4404	return NULL;
4405	}
4406
4407	/* Base case of isl_tab_basic_map_partial_lexopt, after removing
4408	* 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 order
4412	* during the construction of the solution map.
4413	* Furthermore, make sure that the known integer divisions
4414	* appear before any unknown integer division because the solution
4415	* may depend on the known integer divisions, while anything that
4416	* depends on any variable starting from the first unknown integer
4417	* division is ignored in sol_pma_add.
4418	*/
4419	static 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(bset: 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(bset: 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(bset: context->op->peek_basic_set(context)));
4444	} else {
4445	tab = tab_for_lexmin(bmap,
4446	dom: context->op->peek_basic_set(context), M: 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;
4455	error:
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 removing
4462	* some obvious symmetries.
4463	*
4464	* We call basic_map_partial_lexopt_base_sol and extract the results.
4465	*/
4466	static __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	init: &sol_map_init);
4476	if (!sol)
4477	return NULL;
4478	sol_map = (struct isl_sol_map *) sol;
4479
4480	result = isl_map_copy(map: sol_map->map);
4481	if (empty)
4482	*empty = isl_set_copy(set: sol_map->empty);
4483	sol_free(sol: &sol_map->sol);
4484	return result;
4485	}
4486
4487	/* Return a count of the number of occurrences of the "n" first
4488	* variables in the inequality constraints of "bmap".
4489	*/
4490	static __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 number
4517	* of occurrences of each of the variables in the inequality constraints.
4518	*/
4519	static 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 constraint
4534	* "ineq" of "bmap" only occur in that constraint?
4535	*/
4536	static 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_in
4557	* n_out: number of "output" variables: isl_dim_out + isl_dim_div
4558	* val: the coefficients of the output variables
4559	*/
4560	struct 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 variables
4567	* of the constraint in "entry" are equal to info->val.
4568	*/
4569	static 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(p1: (*row) + 1 + info->n_in, p2: info->val, len: info->n_out);
4576	return isl_bool_ok(b: eq);
4577	}
4578
4579	/* Check whether "bmap" has a pair of constraints that have
4580	* the same coefficients for the output variables.
4581	* Note that the coefficients of the existentially quantified
4582	* variables need to be zero since the existentially quantified
4583	* of the result are usually not the same as those of the input.
4584	* Furthermore, check that each of the input variables that occur
4585	* in those constraints does not occur in any other constraint.
4586	* If so, return true and return the row indices of the two constraints
4587	* in *first and *second.
4588	*/
4589	static 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, min_size: bmap->n_ineq);
4602	if (!table)
4603	goto error;
4604
4605	nparam = isl_basic_map_dim(bmap, type: isl_dim_param);
4606	n_in = isl_basic_map_dim(bmap, type: isl_dim_in);
4607	n_out = isl_basic_map_dim(bmap, type: isl_dim_out);
4608	n_div = isl_basic_map_dim(bmap, type: 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, n: 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(p: info.val, len: n_out) < 0)
4621	continue;
4622	if (isl_seq_first_non_zero(p: info.val + n_out, len: n_div) >= 0)
4623	continue;
4624	if (!single_occurrence(n: info.n_in, c: bmap->ineq[i] + 1,
4625	occurrences))
4626	continue;
4627	hash = isl_seq_get_hash(p: info.val, len: info.n_out);
4628	entry = isl_hash_table_find(ctx, table, key_hash: hash,
4629	eq: constraint_equal, val: &info, reserve: 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(ptr: occurrences);
4644
4645	return isl_bool_ok(b: i < bmap->n_ineq);
4646	error:
4647	isl_hash_table_free(ctx, table);
4648	free(ptr: 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 constraints
4656	*
4657	* b_i <= b_j for j > i
4658	* b_i < b_j for j < i
4659	*/
4660	static __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(mat: 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(dst: bset->ineq[k], m1: ctx->one, src1: var->row[j],
4675	m2: ctx->negone, src2: var->row[i], len: 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;
4684	error:
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 one
4692	* of the upper bounds is smaller than all the others and assign
4693	* this upper bound to m.
4694	*
4695	* In particular, if there are n bounds b_i, then the result
4696	* consists of n basic sets, each one of the form
4697	*
4698	* m = b_i
4699	* b_i <= b_j for j > i
4700	* b_i < b_j for j < i
4701	*/
4702	static __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(space: isl_space_copy(space),
4713	n: var->n_row, ISL_SET_DISJOINT);
4714
4715	for (i = 0; i < var->n_row; ++i) {
4716	bset = isl_basic_set_alloc_space(space: isl_space_copy(space), extra: 0,
4717	n_eq: 1, n_ineq: var->n_row - 1);
4718	k = isl_basic_set_alloc_equality(bset);
4719	if (k < 0)
4720	goto error;
4721	isl_seq_cpy(dst: bset->eq[k], src: var->row[i], len: 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(mat: var);
4729	return set;
4730	error:
4731	isl_basic_set_free(bset);
4732	isl_set_free(set);
4733	isl_space_free(space);
4734	isl_mat_free(mat: var);
4735	return NULL;
4736	}
4737
4738	/* Given that the last input variable of "bmap" represents the minimum
4739	* of the bounds in "cst", check whether we need to split the domain
4740	* based on which bound attains the minimum.
4741	*
4742	* A split is needed when the minimum appears in an integer division
4743	* or in an equality. Otherwise, it is only needed if it appears in
4744	* an upper bound that is different from the upper bounds on which it
4745	* is defined.
4746	*/
4747	static 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, type: 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(p: bmap->ineq[i] + 1 + pos + 1,
4773	len: total - pos - 1) >= 0)
4774	return isl_bool_true;
4775
4776	for (j = 0; j < cst->n_row; ++j)
4777	if (isl_seq_eq(p1: bmap->ineq[i], p2: cst->row[j], len: 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 minimum
4787	* of the bounds in "cst", check whether we need to split the domain
4788	* based on which bound attains the minimum.
4789	*
4790	* We simply call need_split_basic_map here. This is safe because
4791	* the position of the minimum is computed from "cst" and not
4792	* from "bmap".
4793	*/
4794	static 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(bmap: bset_to_bmap(bset), cst);
4798	}
4799
4800	/* Given that the last set variable of "set" represents the minimum
4801	* of the bounds in "cst", check whether we need to split the domain
4802	* based on which bound attains the minimum.
4803	*/
4804	static 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(bset: 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 minimum
4820	* of the bounds in "cst", split each basic set in the set
4821	* 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 only
4827	* 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 m
4829	* to obtain l <= a and l <= b, without having to split on whether
4830	* m is equal to a or b.
4831	*/
4832	static __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(map: opt, type: isl_dim_in);
4841	if (n_in < 0 || !min_expr || !cst)
4842	goto error;
4843
4844	space = isl_map_get_space(map: opt);
4845	space = isl_space_drop_dims(space, type: isl_dim_in, first: n_in - 1, num: 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(bmap: isl_basic_map_copy(bmap: opt->p[i]));
4853	split = need_split_basic_map(bmap: 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	set: isl_set_copy(set: min_expr));
4859	map = isl_map_remove_dims(map, type: isl_dim_in, first: n_in - 1, n: 1);
4860
4861	res = isl_map_union_disjoint(map1: res, map2: map);
4862	}
4863
4864	isl_map_free(map: opt);
4865	isl_set_free(set: min_expr);
4866	isl_mat_free(mat: cst);
4867	return res;
4868	error:
4869	isl_map_free(map: opt);
4870	isl_set_free(set: min_expr);
4871	isl_mat_free(mat: cst);
4872	return NULL;
4873	}
4874
4875	/* Given a set of which the last set variable is the minimum
4876	* of the bounds in "cst", split each basic set in the set
4877	* 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	*/
4881	static __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(set: empty);
4887	map = split_domain(opt: map, min_expr, cst);
4888	empty = isl_map_domain(bmap: map);
4889
4890	return empty;
4891	}
4892
4893	static __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 minimum
4899	* of the bounds in "cst". "map_space" is the space of the original
4900	* 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 in
4904	* the definition of the minimum in the result.
4905	*/
4906	static __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(space: isl_basic_set_get_space(bset: dom), var: isl_mat_copy(mat: cst));
4915
4916	opt = basic_map_partial_lexopt(bmap, dom, empty, max);
4917
4918	if (empty) {
4919	*empty = split(empty: *empty,
4920	min_expr: isl_set_copy(set: min_expr), cst: isl_mat_copy(mat: cst));
4921	*empty = isl_set_reset_space(set: *empty, space: set_space);
4922	}
4923
4924	opt = split_domain(opt, min_expr, cst);
4925	opt = isl_map_reset_space(map: opt, space: map_space);
4926
4927	return opt;
4928	}
4929
4930	/* Extract a domain from "bmap" for the purpose of computing
4931	* a lexicographic optimum.
4932	*
4933	* This function is only called when the caller wants to compute a full
4934	* lexicographic optimum, i.e., without specifying a domain. In this case,
4935	* the caller is not interested in the part of the domain space where
4936	* there is no solution and the domain can be initialized to those constraints
4937	* of "bmap" that only involve the parameters and the input dimensions.
4938	* This relieves the parametric programming engine from detecting those
4939	* inequalities and transferring them to the context. More importantly,
4940	* it ensures that those inequalities are transferred first and not
4941	* intermixed with inequalities that actually split the domain.
4942	*
4943	* If the caller does not require the absence of existentially quantified
4944	* 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 that
4946	* the domain does not need to be split at all just to separate out
4947	* 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 involve
4949	* (unknown) existentially quantified variables which will then also
4950	* appear in the solution.
4951	*/
4952	static __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, type: isl_dim_div);
4959	n_out = isl_basic_map_dim(bmap, type: 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	type: isl_dim_div, first: 0, n: n_div);
4966	bmap = isl_basic_map_drop_constraints_involving_dims(bmap,
4967	type: isl_dim_out, first: 0, n: n_out);
4968	}
4969	return isl_basic_map_domain(bmap);
4970	}
4971
4972	#undef TYPE
4973	#define TYPE isl_map
4974	#undef SUFFIX
4975	#define SUFFIX
4976	#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	*/
4981	static __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, size: 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 sequence
5010	* equal to the zero vector?
5011	*/
5012	static 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(mat: 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(mat: trivial);
5027	if (len < 0)
5028	return isl_bool_error;
5029	v = extract_sample_sequence(tab, pos, len);
5030	v = isl_mat_vec_product(mat: isl_mat_copy(mat: trivial), vec: v);
5031	is_trivial = isl_vec_is_zero(vec: v);
5032	isl_vec_free(vec: 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 passed
5042	* 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 each
5046	* (potential) level of the backtracking procedure of
5047	* isl_tab_basic_set_non_trivial_lexmin.
5048	* "v" is a pre-allocated vector that can be used for adding
5049	* 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	*/
5054	struct 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 regions
5067	* are non-trivial or -1 in case of error.
5068	*/
5069	static 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(tab: data->tab, pos: data->region[i].pos,
5076	trivial: 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 first
5087	* n_op entries are zero.
5088	*/
5089	static 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 significantly
5100	* better than that represented by "sol". That is, find the first
5101	* relevant (within first n_op) non-zero coefficient and force it (along
5102	* 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 zero
5106	* 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 and
5110	* at least 2 * (n_op - n_zero) more elements in the constraint array
5111	* are available in the tableau.
5112	*/
5113	static 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(vec: sol);
5134	v = isl_vec_alloc(ctx, size: 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(vec: v);
5141	isl_int_set_si(v->el[1 + i], -1);
5142	if (add_lexmin_eq(tab, eq: v->el) < 0)
5143	goto error;
5144	}
5145
5146	isl_vec_free(vec: v);
5147	return n;
5148	error:
5149	isl_vec_free(vec: 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 least
5156	* two more elements in the constraint array are available in the tableau.
5157	*/
5158	static 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(vec: data->v);
5164	if (!data->v)
5165	return isl_stat_error;
5166	len = isl_mat_cols(mat: region->trivial);
5167	if (len < 0)
5168	return isl_stat_error;
5169	isl_seq_cpy(dst: data->v->el + 1 + region->pos, src: region->trivial->row[dir],
5170	len);
5171	if (add_lexmin_eq(tab, eq: 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 positive
5180	* if side is even and made negative if side is odd.
5181	*
5182	* This function assumes that at least one more row and at least
5183	* one more element in the constraint array are available in the tableau.
5184	*/
5185	static 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(vec: data->v);
5192	if (!data->v)
5193	goto error;
5194	isl_int_set_si(data->v->el[0], -1);
5195	len = isl_mat_cols(mat: region->trivial);
5196	if (len < 0)
5197	goto error;
5198	if (side % 2 == 0)
5199	isl_seq_cpy(dst: data->v->el + 1 + region->pos,
5200	src: region->trivial->row[side / 2], len);
5201	else
5202	isl_seq_neg(dst: data->v->el + 1 + region->pos,
5203	src: region->trivial->row[side / 2], len);
5204	return add_lexmin_ineq(tab, ineq: data->v->el);
5205	error:
5206	isl_tab_free(tab);
5207	return NULL;
5208	}
5209
5210	/* Local data at each level of the backtracking procedure of
5211	* isl_tab_basic_set_non_trivial_lexmin.
5212	*
5213	* "update" is set if a solution has been found in the current case
5214	* of this level, such that a better solution needs to be enforced
5215	* in the next case.
5216	* "n_zero" is the number of initial coordinates that have already
5217	* 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 needs
5222	* to be satisfied by all subsequent cases.
5223	*/
5224	struct 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 solving
5234	* the ILP problem "bset".
5235	*/
5236	static 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(bmap: bset, NULL, M: 0, max: 0);
5244	if (!data->tab)
5245	return isl_stat_error;
5246
5247	data->v = isl_vec_alloc(ctx, size: 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, size: 0);
5256
5257	return isl_stat_ok;
5258	}
5259
5260	/* Mark all outer levels as requiring a better solution
5261	* in the next cases.
5262	*/
5263	static 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" and
5272	* to initiate processing at this level.
5273	*/
5274	static 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(mat: 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; abort
5293	* done: an optimal solution has been found; stop search
5294	* backtrack: backtrack to the previous level
5295	* handle: add the constraints for the current level and
5296	* move to the next level
5297	*/
5298	enum 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 imposed
5309	* in all subsequent cases. This means that if the current
5310	* tableau is empty, then none of those cases should be considered
5311	* anymore and all cases have effectively been considered.
5312	*/
5313	static 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 out
5322	* what to do next. "init" is set if the level was entered
5323	* from a higher level and needs to be initialized.
5324	* Otherwise, the level is entered as a result of backtracking and
5325	* the tableau needs to be restored to a position that can
5326	* 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 initialized
5331	* to point to the first violated region.
5332	* If the constraints of all regions are satisfied by the current
5333	* sample of the tableau, then tell the caller to continue looking
5334	* for a better solution or to stop searching if an optimal solution
5335	* has been found.
5336	*
5337	* If the tableau is empty or if all cases at the current level
5338	* have been considered, then the caller needs to backtrack as well.
5339	*/
5340	static 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(tab: 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(vec: data->sol);
5359	data->sol = isl_tab_get_sample_value(tab: data->tab);
5360	if (!data->sol)
5361	return isl_next_error;
5362	if (is_optimal(sol: data->sol, n_op: 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, region: r, data) < 0)
5371	return isl_next_error;
5372	if (isl_tab_extend_cons(tab: data->tab,
5373	n_new: 2 * local->n + 2 * data->n_op) < 0)
5374	return isl_next_error;
5375	} else {
5376	if (isl_tab_rollback(tab: data->tab, snap: 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 level
5386	* (marked by local->update being set), then add constraints
5387	* that enforce a better solution in the present and all following cases.
5388	* The constraints only need to be imposed once because they are
5389	* included in the snapshot (taken in pick_side) that will be used in
5390	* subsequent cases.
5391	*/
5392	static 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(tab: data->tab,
5399	sol: data->sol, n_op: data->n_op, n_zero: 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 are
5412	* v_0 >= 1
5413	* v_0 <= -1
5414	* v_0 = 0 and v_1 >= 1
5415	* v_0 = 0 and v_1 <= -1
5416	* v_0 = 0 and v_1 = 0 and v_2 >= 1
5417	* v_0 = 0 and v_1 = 0 and v_2 <= -1
5418	* ...
5419	* in this order.
5420	*
5421	* A snapshot is taken after the equality constraint (if any) has been added
5422	* such that the next case can start off from this position.
5423	* The rollback to this position is performed in enter_level.
5424	*/
5425	static 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(tab: data->tab, region, dir: base / 2 - 1, data) < 0)
5437	return isl_stat_error;
5438
5439	local->snap = isl_tab_snap(tab: data->tab);
5440	if (isl_tab_push_basis(tab: data->tab) < 0)
5441	return isl_stat_error;
5442
5443	data->tab = pos_neg(tab: 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	*/
5451	static void clear_lexmin_data(struct isl_lexmin_data *data)
5452	{
5453	free(ptr: data->local);
5454	isl_vec_free(vec: data->v);
5455	isl_tab_free(tab: data->tab);
5456	}
5457
5458	/* Return the lexicographically smallest non-trivial solution of the
5459	* 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 looking
5465	* for solutions that result in significantly better values for those
5466	* initial coordinates. That is, we only continue looking for solutions
5467	* 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 the
5470	* specified regions. Each region represents a sequence of
5471	* triviality directions on a sequence of variables that starts
5472	* at a given position. A solution is non-trivial on such a region if
5473	* at least one of the triviality directions is non-zero
5474	* on that sequence of variables.
5475	*
5476	* Whenever a conflict is encountered, all constraints involved are
5477	* 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 region
5481	* that is forced to be non-trivial.
5482	* At each level we consider 2 * n cases, where n
5483	* is the number of triviality directions.
5484	* In terms of those n directions v_i, we consider the cases
5485	* v_0 >= 1
5486	* v_0 <= -1
5487	* v_0 = 0 and v_1 >= 1
5488	* v_0 = 0 and v_1 <= -1
5489	* v_0 = 0 and v_1 = 0 and v_2 >= 1
5490	* v_0 = 0 and v_1 = 0 and v_2 <= -1
5491	* ...
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: &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: &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: &data) < 0)
5529	goto error;
5530	if (pick_side(local, data: &data) < 0)
5531	goto error;
5532
5533	local->side++;
5534	level++;
5535	init = 1;
5536	}
5537
5538	clear_lexmin_data(data: &data);
5539	isl_basic_set_free(bset);
5540
5541	return data.sol;
5542	error:
5543	clear_lexmin_data(data: &data);
5544	isl_basic_set_free(bset);
5545	isl_vec_free(vec: data.sol);
5546	return NULL;
5547	}
5548
5549	/* Wrapper for a tableau that is used for computing
5550	* 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	*/
5554	struct 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(ctx: tl->ctx);
5566	isl_tab_free(tab: tl->tab);
5567	free(ptr: tl);
5568
5569	return NULL;
5570	}
5571
5572	/* Construct an isl_tab_lexmin for computing
5573	* 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(bmap: bset, NULL, M: 0, max: 0);
5592	isl_basic_set_free(bset);
5593	if (!tl->tab)
5594	return isl_tab_lexmin_free(tl);
5595	return tl;
5596	error:
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	*/
5604	int 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 optimal
5610	* 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(tab: tl->tab, n_new: 2) < 0)
5622	return isl_tab_lexmin_free(tl);
5623	n_var = tl->tab->n_var;
5624	isl_seq_neg(dst: eq, src: eq, len: 1 + n_var);
5625	tl->tab = add_lexmin_ineq(tab: tl->tab, ineq: eq);
5626	isl_seq_neg(dst: eq, src: eq, len: 1 + n_var);
5627	tl->tab = add_lexmin_ineq(tab: tl->tab, ineq: 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 or
5636	* 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(tab: 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 set
5650	* 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(ctx: tl->ctx, size: 0);
5659	else
5660	return isl_tab_get_sample_value(tab: tl->tab);
5661	}
5662
5663	struct isl_sol_pma {
5664	struct isl_sol sol;
5665	isl_pw_multi_aff *pma;
5666	isl_set *empty;
5667	};
5668
5669	static 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(pma: sol_pma->pma);
5673	isl_set_free(set: sol_pma->empty);
5674	}
5675
5676	/* This function is called for parts of the context where there is
5677	* no solution, with "bset" corresponding to the context tableau.
5678	* Simply add the basic set to the set "empty".
5679	*/
5680	static 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(set: sol->empty, n: 1);
5687	bset = isl_basic_set_simplify(bset);
5688	bset = isl_basic_set_finalize(bset);
5689	sol->empty = isl_set_add_basic_set(set: sol->empty, bset);
5690	if (!sol->empty)
5691	sol->sol.error = 1;
5692	return;
5693	error:
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 of
5699	* affine expressions "maff" defined over this domain, construct
5700	* an isl_pw_multi_aff with a single cell corresponding to "dom" and
5701	* the affine expressions in "maff".
5702	*/
5703	static 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(bset: dom);
5709	dom = isl_basic_set_finalize(bset: dom);
5710	pma = isl_pw_multi_aff_alloc(set: isl_set_from_basic_set(bset: dom), maff);
5711	sol->pma = isl_pw_multi_aff_add_disjoint(pma1: sol->pma, pma2: pma);
5712	if (!sol->pma)
5713	sol->sol.error = 1;
5714	}
5715
5716	static void sol_pma_add_empty_wrap(struct isl_sol *sol,
5717	__isl_take isl_basic_set *bset)
5718	{
5719	sol_pma_add_empty(sol: (struct isl_sol_pma *)sol, bset);
5720	}
5721
5722	static 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(sol: (struct isl_sol_pma *)sol, dom, maff: 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 parts
5730	* of the context where there is no solution.
5731	* If max is set, then we are solving a maximization, rather than
5732	* a minimization problem, which means that the variables in the
5733	* tableau have value "M - x" rather than "M + x".
5734	*/
5735	static 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: &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(space: 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(space: isl_basic_set_get_space(bset: dom),
5760	n: 1, ISL_SET_DISJOINT);
5761	if (!sol_pma->empty)
5762	goto error;
5763	}
5764
5765	isl_basic_set_free(bset: dom);
5766	return &sol_pma->sol;
5767	error:
5768	isl_basic_set_free(bset: dom);
5769	sol_free(sol: &sol_pma->sol);
5770	return NULL;
5771	}
5772
5773	/* Base case of isl_tab_basic_map_partial_lexopt, after removing
5774	* some obvious symmetries.
5775	*
5776	* We call basic_map_partial_lexopt_base_sol and extract the results.
5777	*/
5778	static __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	init: &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(pma: sol_pma->pma);
5793	if (empty)
5794	*empty = isl_set_copy(set: sol_pma->empty);
5795	sol_free(sol: &sol_pma->sol);
5796	return result;
5797	}
5798
5799	/* Given that the last input variable of "maff" represents the minimum
5800	* of some bounds, check whether we need to plug in the expression
5801	* of the minimum.
5802	*
5803	* In particular, check if the last input variable appears in any
5804	* of the expressions in "maff".
5805	*/
5806	static 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(multi: maff, type: 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(aff: maff->u.p[i],
5821	type: isl_dim_in, first: pos, n: 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 selects
5831	* the minimal upper bound to m, i.e.,
5832	* divide the space into cells where one
5833	* of the upper bounds is smaller than all the others and select
5834	* this upper bound on that cell.
5835	*
5836	* In particular, if there are n bounds b_i, then the result
5837	* consists of n cell, each one of the form
5838	*
5839	* b_i <= b_j for j > i
5840	* b_i < b_j for j < i
5841	*
5842	* The affine expression on this cell is
5843	*
5844	* b_i
5845	*/
5846	static __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(space: isl_space_copy(space));
5860	pw_space = isl_space_copy(space);
5861	pw_space = isl_space_from_domain(space: pw_space);
5862	pw_space = isl_space_add_dims(space: pw_space, type: isl_dim_out, n: 1);
5863	paff = isl_pw_aff_alloc_size(space: pw_space, n: 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(ls: isl_local_space_copy(ls));
5869	bset = isl_basic_set_alloc_space(space: isl_space_copy(space), extra: 0,
5870	n_eq: 0, n_ineq: 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(dst: aff->v->el + 1, src: var->row[i], len: 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(set: isl_set_from_basic_set(bset), aff);
5878	paff = isl_pw_aff_add_disjoint(pwaff1: paff, pwaff2: paff_i);
5879	}
5880
5881	isl_local_space_free(ls);
5882	isl_space_free(space);
5883	isl_mat_free(mat: var);
5884	return paff;
5885	error:
5886	isl_aff_free(aff);
5887	isl_basic_set_free(bset);
5888	isl_pw_aff_free(pwaff: paff);
5889	isl_local_space_free(ls);
5890	isl_space_free(space);
5891	isl_mat_free(mat: var);
5892	return NULL;
5893	}
5894
5895	/* Given a piecewise multi-affine expression of which the last input variable
5896	* 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 expression
5903	* is plugged in. Otherwise, if the variable appears in the constraints
5904	* and a split is required, then the domain is split. Otherwise, no split
5905	* is performed.
5906	*/
5907	static __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(pma: opt, type: isl_dim_in);
5920	if (n_in < 0)
5921	goto error;
5922	space = isl_pw_multi_aff_get_space(pma: opt);
5923	space = isl_space_drop_dims(space, type: isl_dim_in, first: n_in - 1, num: 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(set: isl_set_copy(set: opt->p[i].set),
5931	maff: isl_multi_aff_copy(multi: opt->p[i].maff));
5932	subs = need_substitution(maff: 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	pos: n_in - 1, subs: min_expr_pa);
5938	} else {
5939	isl_bool split;
5940	split = need_split_set(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	set: isl_set_copy(set: min_expr));
5946	}
5947	pma = isl_pw_multi_aff_project_out(pma,
5948	type: isl_dim_in, first: n_in - 1, n: 1);
5949
5950	res = isl_pw_multi_aff_add_disjoint(pma1: res, pma2: pma);
5951	}
5952
5953	isl_pw_multi_aff_free(pma: opt);
5954	isl_pw_aff_free(pwaff: min_expr_pa);
5955	isl_set_free(set: min_expr);
5956	isl_mat_free(mat: cst);
5957	return res;
5958	error:
5959	isl_pw_multi_aff_free(pma: opt);
5960	isl_pw_aff_free(pwaff: min_expr_pa);
5961	isl_set_free(set: min_expr);
5962	isl_mat_free(mat: cst);
5963	return NULL;
5964	}
5965
5966	static __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 minimum
5972	* of the bounds in "cst". "map_space" is the space of the original
5973	* 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 in
5977	* the definition of the minimum in the result.
5978	*/
5979	static __isl_give isl_pw_multi_aff *
5980	basic_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(space: isl_basic_set_get_space(bset: dom), var: isl_mat_copy(mat: cst));
5990	min_expr_pa = set_minimum_pa(space: isl_basic_set_get_space(bset: dom),
5991	var: isl_mat_copy(mat: cst));
5992
5993	opt = basic_map_partial_lexopt_pw_multi_aff(bmap, dom, empty, max);
5994
5995	if (empty) {
5996	*empty = split(empty: *empty,
5997	min_expr: isl_set_copy(set: min_expr), cst: isl_mat_copy(mat: cst));
5998	*empty = isl_set_reset_space(set: *empty, space: set_space);
5999	}
6000
6001	opt = split_domain_pma(opt, min_expr_pa, min_expr, cst);
6002	opt = isl_pw_multi_aff_reset_space(pwmaff: opt, space: map_space);
6003
6004	return opt;
6005	}
6006
6007	#undef TYPE
6008	#define TYPE isl_pw_multi_aff
6009	#undef SUFFIX
6010	#define SUFFIX _pw_multi_aff
6011	#include "isl_tab_lexopt_templ.c"
6012