[AMDGPU] Test codegen'ing True16 additions.
[llvm-project.git] / polly / lib / External / isl / isl_equalities.c
blob90dfb7709ce7ec46942dec864cb2d9de846228b1
1 /*
2 * Copyright 2008-2009 Katholieke Universiteit Leuven
3 * Copyright 2010 INRIA Saclay
5 * Use of this software is governed by the MIT license
7 * Written by Sven Verdoolaege, K.U.Leuven, Departement
8 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
9 * and INRIA Saclay - Ile-de-France, Parc Club Orsay Universite,
10 * ZAC des vignes, 4 rue Jacques Monod, 91893 Orsay, France
13 #include <isl_mat_private.h>
14 #include <isl_vec_private.h>
15 #include <isl_seq.h>
16 #include "isl_map_private.h"
17 #include "isl_equalities.h"
18 #include <isl_val_private.h>
20 /* Given a set of modulo constraints
22 * c + A y = 0 mod d
24 * this function computes a particular solution y_0
26 * The input is given as a matrix B = [ c A ] and a vector d.
28 * The output is matrix containing the solution y_0 or
29 * a zero-column matrix if the constraints admit no integer solution.
31 * The given set of constrains is equivalent to
33 * c + A y = -D x
35 * with D = diag d and x a fresh set of variables.
36 * Reducing both c and A modulo d does not change the
37 * value of y in the solution and may lead to smaller coefficients.
38 * Let M = [ D A ] and [ H 0 ] = M U, the Hermite normal form of M.
39 * Then
40 * [ x ]
41 * M [ y ] = - c
42 * and so
43 * [ x ]
44 * [ H 0 ] U^{-1} [ y ] = - c
45 * Let
46 * [ A ] [ x ]
47 * [ B ] = U^{-1} [ y ]
48 * then
49 * H A + 0 B = -c
51 * so B may be chosen arbitrarily, e.g., B = 0, and then
53 * [ x ] = [ -c ]
54 * U^{-1} [ y ] = [ 0 ]
55 * or
56 * [ x ] [ -c ]
57 * [ y ] = U [ 0 ]
58 * specifically,
60 * y = U_{2,1} (-c)
62 * If any of the coordinates of this y are non-integer
63 * then the constraints admit no integer solution and
64 * a zero-column matrix is returned.
66 static __isl_give isl_mat *particular_solution(__isl_keep isl_mat *B,
67 __isl_keep isl_vec *d)
69 int i, j;
70 struct isl_mat *M = NULL;
71 struct isl_mat *C = NULL;
72 struct isl_mat *U = NULL;
73 struct isl_mat *H = NULL;
74 struct isl_mat *cst = NULL;
75 struct isl_mat *T = NULL;
77 M = isl_mat_alloc(B->ctx, B->n_row, B->n_row + B->n_col - 1);
78 C = isl_mat_alloc(B->ctx, 1 + B->n_row, 1);
79 if (!M || !C)
80 goto error;
81 isl_int_set_si(C->row[0][0], 1);
82 for (i = 0; i < B->n_row; ++i) {
83 isl_seq_clr(M->row[i], B->n_row);
84 isl_int_set(M->row[i][i], d->block.data[i]);
85 isl_int_neg(C->row[1 + i][0], B->row[i][0]);
86 isl_int_fdiv_r(C->row[1+i][0], C->row[1+i][0], M->row[i][i]);
87 for (j = 0; j < B->n_col - 1; ++j)
88 isl_int_fdiv_r(M->row[i][B->n_row + j],
89 B->row[i][1 + j], M->row[i][i]);
91 M = isl_mat_left_hermite(M, 0, &U, NULL);
92 if (!M || !U)
93 goto error;
94 H = isl_mat_sub_alloc(M, 0, B->n_row, 0, B->n_row);
95 H = isl_mat_lin_to_aff(H);
96 C = isl_mat_inverse_product(H, C);
97 if (!C)
98 goto error;
99 for (i = 0; i < B->n_row; ++i) {
100 if (!isl_int_is_divisible_by(C->row[1+i][0], C->row[0][0]))
101 break;
102 isl_int_divexact(C->row[1+i][0], C->row[1+i][0], C->row[0][0]);
104 if (i < B->n_row)
105 cst = isl_mat_alloc(B->ctx, B->n_row, 0);
106 else
107 cst = isl_mat_sub_alloc(C, 1, B->n_row, 0, 1);
108 T = isl_mat_sub_alloc(U, B->n_row, B->n_col - 1, 0, B->n_row);
109 cst = isl_mat_product(T, cst);
110 isl_mat_free(M);
111 isl_mat_free(C);
112 isl_mat_free(U);
113 return cst;
114 error:
115 isl_mat_free(M);
116 isl_mat_free(C);
117 isl_mat_free(U);
118 return NULL;
121 /* Compute and return the matrix
123 * U_1^{-1} diag(d_1, 1, ..., 1)
125 * with U_1 the unimodular completion of the first (and only) row of B.
126 * The columns of this matrix generate the lattice that satisfies
127 * the single (linear) modulo constraint.
129 static __isl_take isl_mat *parameter_compression_1(__isl_keep isl_mat *B,
130 __isl_keep isl_vec *d)
132 struct isl_mat *U;
134 U = isl_mat_alloc(B->ctx, B->n_col - 1, B->n_col - 1);
135 if (!U)
136 return NULL;
137 isl_seq_cpy(U->row[0], B->row[0] + 1, B->n_col - 1);
138 U = isl_mat_unimodular_complete(U, 1);
139 U = isl_mat_right_inverse(U);
140 if (!U)
141 return NULL;
142 isl_mat_col_mul(U, 0, d->block.data[0], 0);
143 U = isl_mat_lin_to_aff(U);
144 return U;
147 /* Compute a common lattice of solutions to the linear modulo
148 * constraints specified by B and d.
149 * See also the documentation of isl_mat_parameter_compression.
150 * We put the matrix
152 * A = [ L_1^{-T} L_2^{-T} ... L_k^{-T} ]
154 * on a common denominator. This denominator D is the lcm of modulos d.
155 * Since L_i = U_i^{-1} diag(d_i, 1, ... 1), we have
156 * L_i^{-T} = U_i^T diag(d_i, 1, ... 1)^{-T} = U_i^T diag(1/d_i, 1, ..., 1).
157 * Putting this on the common denominator, we have
158 * D * L_i^{-T} = U_i^T diag(D/d_i, D, ..., D).
160 static __isl_give isl_mat *parameter_compression_multi(__isl_keep isl_mat *B,
161 __isl_keep isl_vec *d)
163 int i, j, k;
164 isl_int D;
165 struct isl_mat *A = NULL, *U = NULL;
166 struct isl_mat *T;
167 unsigned size;
169 isl_int_init(D);
171 isl_vec_lcm(d, &D);
173 size = B->n_col - 1;
174 A = isl_mat_alloc(B->ctx, size, B->n_row * size);
175 U = isl_mat_alloc(B->ctx, size, size);
176 if (!U || !A)
177 goto error;
178 for (i = 0; i < B->n_row; ++i) {
179 isl_seq_cpy(U->row[0], B->row[i] + 1, size);
180 U = isl_mat_unimodular_complete(U, 1);
181 if (!U)
182 goto error;
183 isl_int_divexact(D, D, d->block.data[i]);
184 for (k = 0; k < U->n_col; ++k)
185 isl_int_mul(A->row[k][i*size+0], D, U->row[0][k]);
186 isl_int_mul(D, D, d->block.data[i]);
187 for (j = 1; j < U->n_row; ++j)
188 for (k = 0; k < U->n_col; ++k)
189 isl_int_mul(A->row[k][i*size+j],
190 D, U->row[j][k]);
192 A = isl_mat_left_hermite(A, 0, NULL, NULL);
193 T = isl_mat_sub_alloc(A, 0, A->n_row, 0, A->n_row);
194 T = isl_mat_lin_to_aff(T);
195 if (!T)
196 goto error;
197 isl_int_set(T->row[0][0], D);
198 T = isl_mat_right_inverse(T);
199 if (!T)
200 goto error;
201 isl_assert(T->ctx, isl_int_is_one(T->row[0][0]), goto error);
202 T = isl_mat_transpose(T);
203 isl_mat_free(A);
204 isl_mat_free(U);
206 isl_int_clear(D);
207 return T;
208 error:
209 isl_mat_free(A);
210 isl_mat_free(U);
211 isl_int_clear(D);
212 return NULL;
215 /* Given a set of modulo constraints
217 * c + A y = 0 mod d
219 * this function returns an affine transformation T,
221 * y = T y'
223 * that bijectively maps the integer vectors y' to integer
224 * vectors y that satisfy the modulo constraints.
226 * This function is inspired by Section 2.5.3
227 * of B. Meister, "Stating and Manipulating Periodicity in the Polytope
228 * Model. Applications to Program Analysis and Optimization".
229 * However, the implementation only follows the algorithm of that
230 * section for computing a particular solution and not for computing
231 * a general homogeneous solution. The latter is incomplete and
232 * may remove some valid solutions.
233 * Instead, we use an adaptation of the algorithm in Section 7 of
234 * B. Meister, S. Verdoolaege, "Polynomial Approximations in the Polytope
235 * Model: Bringing the Power of Quasi-Polynomials to the Masses".
237 * The input is given as a matrix B = [ c A ] and a vector d.
238 * Each element of the vector d corresponds to a row in B.
239 * The output is a lower triangular matrix.
240 * If no integer vector y satisfies the given constraints then
241 * a matrix with zero columns is returned.
243 * We first compute a particular solution y_0 to the given set of
244 * modulo constraints in particular_solution. If no such solution
245 * exists, then we return a zero-columned transformation matrix.
246 * Otherwise, we compute the generic solution to
248 * A y = 0 mod d
250 * That is we want to compute G such that
252 * y = G y''
254 * with y'' integer, describes the set of solutions.
256 * We first remove the common factors of each row.
257 * In particular if gcd(A_i,d_i) != 1, then we divide the whole
258 * row i (including d_i) by this common factor. If afterwards gcd(A_i) != 1,
259 * then we divide this row of A by the common factor, unless gcd(A_i) = 0.
260 * In the later case, we simply drop the row (in both A and d).
262 * If there are no rows left in A, then G is the identity matrix. Otherwise,
263 * for each row i, we now determine the lattice of integer vectors
264 * that satisfies this row. Let U_i be the unimodular extension of the
265 * row A_i. This unimodular extension exists because gcd(A_i) = 1.
266 * The first component of
268 * y' = U_i y
270 * needs to be a multiple of d_i. Let y' = diag(d_i, 1, ..., 1) y''.
271 * Then,
273 * y = U_i^{-1} diag(d_i, 1, ..., 1) y''
275 * for arbitrary integer vectors y''. That is, y belongs to the lattice
276 * generated by the columns of L_i = U_i^{-1} diag(d_i, 1, ..., 1).
277 * If there is only one row, then G = L_1.
279 * If there is more than one row left, we need to compute the intersection
280 * of the lattices. That is, we need to compute an L such that
282 * L = L_i L_i' for all i
284 * with L_i' some integer matrices. Let A be constructed as follows
286 * A = [ L_1^{-T} L_2^{-T} ... L_k^{-T} ]
288 * and computed the Hermite Normal Form of A = [ H 0 ] U
289 * Then,
291 * L_i^{-T} = H U_{1,i}
293 * or
295 * H^{-T} = L_i U_{1,i}^T
297 * In other words G = L = H^{-T}.
298 * To ensure that G is lower triangular, we compute and use its Hermite
299 * normal form.
301 * The affine transformation matrix returned is then
303 * [ 1 0 ]
304 * [ y_0 G ]
306 * as any y = y_0 + G y' with y' integer is a solution to the original
307 * modulo constraints.
309 __isl_give isl_mat *isl_mat_parameter_compression(__isl_take isl_mat *B,
310 __isl_take isl_vec *d)
312 int i;
313 struct isl_mat *cst = NULL;
314 struct isl_mat *T = NULL;
315 isl_int D;
317 if (!B || !d)
318 goto error;
319 isl_assert(B->ctx, B->n_row == d->size, goto error);
320 cst = particular_solution(B, d);
321 if (!cst)
322 goto error;
323 if (cst->n_col == 0) {
324 T = isl_mat_alloc(B->ctx, B->n_col, 0);
325 isl_mat_free(cst);
326 isl_mat_free(B);
327 isl_vec_free(d);
328 return T;
330 isl_int_init(D);
331 /* Replace a*g*row = 0 mod g*m by row = 0 mod m */
332 for (i = 0; i < B->n_row; ++i) {
333 isl_seq_gcd(B->row[i] + 1, B->n_col - 1, &D);
334 if (isl_int_is_one(D))
335 continue;
336 if (isl_int_is_zero(D)) {
337 B = isl_mat_drop_rows(B, i, 1);
338 d = isl_vec_cow(d);
339 if (!B || !d)
340 goto error2;
341 isl_seq_cpy(d->block.data+i, d->block.data+i+1,
342 d->size - (i+1));
343 d->size--;
344 i--;
345 continue;
347 B = isl_mat_cow(B);
348 if (!B)
349 goto error2;
350 isl_seq_scale_down(B->row[i] + 1, B->row[i] + 1, D, B->n_col-1);
351 isl_int_gcd(D, D, d->block.data[i]);
352 d = isl_vec_cow(d);
353 if (!d)
354 goto error2;
355 isl_int_divexact(d->block.data[i], d->block.data[i], D);
357 isl_int_clear(D);
358 if (B->n_row == 0)
359 T = isl_mat_identity(B->ctx, B->n_col);
360 else if (B->n_row == 1)
361 T = parameter_compression_1(B, d);
362 else
363 T = parameter_compression_multi(B, d);
364 T = isl_mat_left_hermite(T, 0, NULL, NULL);
365 if (!T)
366 goto error;
367 isl_mat_sub_copy(T->ctx, T->row + 1, cst->row, cst->n_row, 0, 0, 1);
368 isl_mat_free(cst);
369 isl_mat_free(B);
370 isl_vec_free(d);
371 return T;
372 error2:
373 isl_int_clear(D);
374 error:
375 isl_mat_free(cst);
376 isl_mat_free(B);
377 isl_vec_free(d);
378 return NULL;
381 /* Given a set of equalities
383 * B(y) + A x = 0 (*)
385 * compute and return an affine transformation T,
387 * y = T y'
389 * that bijectively maps the integer vectors y' to integer
390 * vectors y that satisfy the modulo constraints for some value of x.
392 * Let [H 0] be the Hermite Normal Form of A, i.e.,
394 * A = [H 0] Q
396 * Then y is a solution of (*) iff
398 * H^-1 B(y) (= - [I 0] Q x)
400 * is an integer vector. Let d be the common denominator of H^-1.
401 * We impose
403 * d H^-1 B(y) = 0 mod d
405 * and compute the solution using isl_mat_parameter_compression.
407 __isl_give isl_mat *isl_mat_parameter_compression_ext(__isl_take isl_mat *B,
408 __isl_take isl_mat *A)
410 isl_ctx *ctx;
411 isl_vec *d;
412 int n_row, n_col;
414 if (!A)
415 return isl_mat_free(B);
417 ctx = isl_mat_get_ctx(A);
418 n_row = A->n_row;
419 n_col = A->n_col;
420 A = isl_mat_left_hermite(A, 0, NULL, NULL);
421 A = isl_mat_drop_cols(A, n_row, n_col - n_row);
422 A = isl_mat_lin_to_aff(A);
423 A = isl_mat_right_inverse(A);
424 d = isl_vec_alloc(ctx, n_row);
425 if (A)
426 d = isl_vec_set(d, A->row[0][0]);
427 A = isl_mat_drop_rows(A, 0, 1);
428 A = isl_mat_drop_cols(A, 0, 1);
429 B = isl_mat_product(A, B);
431 return isl_mat_parameter_compression(B, d);
434 /* Return a compression matrix that indicates that there are no solutions
435 * to the original constraints. In particular, return a zero-column
436 * matrix with 1 + dim rows. If "T2" is not NULL, then assign *T2
437 * the inverse of this matrix. *T2 may already have been assigned
438 * matrix, so free it first.
439 * "free1", "free2" and "free3" are temporary matrices that are
440 * not useful when an empty compression is returned. They are
441 * simply freed.
443 static __isl_give isl_mat *empty_compression(isl_ctx *ctx, unsigned dim,
444 __isl_give isl_mat **T2, __isl_take isl_mat *free1,
445 __isl_take isl_mat *free2, __isl_take isl_mat *free3)
447 isl_mat_free(free1);
448 isl_mat_free(free2);
449 isl_mat_free(free3);
450 if (T2) {
451 isl_mat_free(*T2);
452 *T2 = isl_mat_alloc(ctx, 0, 1 + dim);
454 return isl_mat_alloc(ctx, 1 + dim, 0);
457 /* Given a matrix that maps a (possibly) parametric domain to
458 * a parametric domain, add in rows that map the "nparam" parameters onto
459 * themselves.
461 static __isl_give isl_mat *insert_parameter_rows(__isl_take isl_mat *mat,
462 unsigned nparam)
464 int i;
466 if (nparam == 0)
467 return mat;
468 if (!mat)
469 return NULL;
471 mat = isl_mat_insert_rows(mat, 1, nparam);
472 if (!mat)
473 return NULL;
475 for (i = 0; i < nparam; ++i) {
476 isl_seq_clr(mat->row[1 + i], mat->n_col);
477 isl_int_set(mat->row[1 + i][1 + i], mat->row[0][0]);
480 return mat;
483 /* Given a set of equalities
485 * -C(y) + M x = 0
487 * this function computes a unimodular transformation from a lower-dimensional
488 * space to the original space that bijectively maps the integer points x'
489 * in the lower-dimensional space to the integer points x in the original
490 * space that satisfy the equalities.
492 * The input is given as a matrix B = [ -C M ] and the output is a
493 * matrix that maps [1 x'] to [1 x].
494 * The number of equality constraints in B is assumed to be smaller than
495 * or equal to the number of variables x.
496 * "first" is the position of the first x variable.
497 * The preceding variables are considered to be y-variables.
498 * If T2 is not NULL, then *T2 is set to a matrix mapping [1 x] to [1 x'].
500 * First compute the (left) Hermite normal form of M,
502 * M [U1 U2] = M U = H = [H1 0]
503 * or
504 * M = H Q = [H1 0] [Q1]
505 * [Q2]
507 * with U, Q unimodular, Q = U^{-1} (and H lower triangular).
508 * Define the transformed variables as
510 * x = [U1 U2] [ x1' ] = [U1 U2] [Q1] x
511 * [ x2' ] [Q2]
513 * The equalities then become
515 * -C(y) + H1 x1' = 0 or x1' = H1^{-1} C(y) = C'(y)
517 * If the denominator of the constant term does not divide the
518 * the common denominator of the coefficients of y, then every
519 * integer point is mapped to a non-integer point and then the original set
520 * has no integer solutions (since the x' are a unimodular transformation
521 * of the x). In this case, a zero-column matrix is returned.
522 * Otherwise, the transformation is given by
524 * x = U1 H1^{-1} C(y) + U2 x2'
526 * The inverse transformation is simply
528 * x2' = Q2 x
530 __isl_give isl_mat *isl_mat_final_variable_compression(__isl_take isl_mat *B,
531 int first, __isl_give isl_mat **T2)
533 int i, n;
534 isl_ctx *ctx;
535 isl_mat *H = NULL, *C, *H1, *U = NULL, *U1, *U2;
536 unsigned dim;
538 if (T2)
539 *T2 = NULL;
540 if (!B)
541 goto error;
543 ctx = isl_mat_get_ctx(B);
544 dim = B->n_col - 1;
545 n = dim - first;
546 if (n < B->n_row)
547 isl_die(ctx, isl_error_invalid, "too many equality constraints",
548 goto error);
549 H = isl_mat_sub_alloc(B, 0, B->n_row, 1 + first, n);
550 H = isl_mat_left_hermite(H, 0, &U, T2);
551 if (!H || !U || (T2 && !*T2))
552 goto error;
553 if (T2) {
554 *T2 = isl_mat_drop_rows(*T2, 0, B->n_row);
555 *T2 = isl_mat_diagonal(isl_mat_identity(ctx, 1 + first), *T2);
556 if (!*T2)
557 goto error;
559 C = isl_mat_alloc(ctx, 1 + B->n_row, 1 + first);
560 if (!C)
561 goto error;
562 isl_int_set_si(C->row[0][0], 1);
563 isl_seq_clr(C->row[0] + 1, first);
564 isl_mat_sub_neg(ctx, C->row + 1, B->row, B->n_row, 0, 0, 1 + first);
565 H1 = isl_mat_sub_alloc(H, 0, H->n_row, 0, H->n_row);
566 H1 = isl_mat_lin_to_aff(H1);
567 C = isl_mat_inverse_product(H1, C);
568 if (!C)
569 goto error;
570 isl_mat_free(H);
571 if (!isl_int_is_one(C->row[0][0])) {
572 isl_int g;
574 isl_int_init(g);
575 for (i = 0; i < B->n_row; ++i) {
576 isl_seq_gcd(C->row[1 + i] + 1, first, &g);
577 isl_int_gcd(g, g, C->row[0][0]);
578 if (!isl_int_is_divisible_by(C->row[1 + i][0], g))
579 break;
581 isl_int_clear(g);
583 if (i < B->n_row)
584 return empty_compression(ctx, dim, T2, B, C, U);
585 C = isl_mat_normalize(C);
587 U1 = isl_mat_sub_alloc(U, 0, U->n_row, 0, B->n_row);
588 U1 = isl_mat_lin_to_aff(U1);
589 U2 = isl_mat_sub_alloc(U, 0, U->n_row, B->n_row, U->n_row - B->n_row);
590 U2 = isl_mat_lin_to_aff(U2);
591 isl_mat_free(U);
592 C = isl_mat_product(U1, C);
593 C = isl_mat_aff_direct_sum(C, U2);
594 C = insert_parameter_rows(C, first);
596 isl_mat_free(B);
598 return C;
599 error:
600 isl_mat_free(B);
601 isl_mat_free(H);
602 isl_mat_free(U);
603 if (T2) {
604 isl_mat_free(*T2);
605 *T2 = NULL;
607 return NULL;
610 /* Given a set of equalities
612 * M x - c = 0
614 * this function computes a unimodular transformation from a lower-dimensional
615 * space to the original space that bijectively maps the integer points x'
616 * in the lower-dimensional space to the integer points x in the original
617 * space that satisfy the equalities.
619 * The input is given as a matrix B = [ -c M ] and the output is a
620 * matrix that maps [1 x'] to [1 x].
621 * The number of equality constraints in B is assumed to be smaller than
622 * or equal to the number of variables x.
623 * If T2 is not NULL, then *T2 is set to a matrix mapping [1 x] to [1 x'].
625 __isl_give isl_mat *isl_mat_variable_compression(__isl_take isl_mat *B,
626 __isl_give isl_mat **T2)
628 return isl_mat_final_variable_compression(B, 0, T2);
631 /* Return "bset" and set *T and *T2 to the identity transformation
632 * on "bset" (provided T and T2 are not NULL).
634 static __isl_give isl_basic_set *return_with_identity(
635 __isl_take isl_basic_set *bset, __isl_give isl_mat **T,
636 __isl_give isl_mat **T2)
638 isl_size dim;
639 isl_mat *id;
641 dim = isl_basic_set_dim(bset, isl_dim_set);
642 if (dim < 0)
643 return isl_basic_set_free(bset);
644 if (!T && !T2)
645 return bset;
647 id = isl_mat_identity(isl_basic_map_get_ctx(bset), 1 + dim);
648 if (T)
649 *T = isl_mat_copy(id);
650 if (T2)
651 *T2 = isl_mat_copy(id);
652 isl_mat_free(id);
654 return bset;
657 /* Use the n equalities of bset to unimodularly transform the
658 * variables x such that n transformed variables x1' have a constant value
659 * and rewrite the constraints of bset in terms of the remaining
660 * transformed variables x2'. The matrix pointed to by T maps
661 * the new variables x2' back to the original variables x, while T2
662 * maps the original variables to the new variables.
664 static __isl_give isl_basic_set *compress_variables(
665 __isl_take isl_basic_set *bset,
666 __isl_give isl_mat **T, __isl_give isl_mat **T2)
668 struct isl_mat *B, *TC;
669 isl_size dim;
671 if (T)
672 *T = NULL;
673 if (T2)
674 *T2 = NULL;
675 if (isl_basic_set_check_no_params(bset) < 0 ||
676 isl_basic_set_check_no_locals(bset) < 0)
677 return isl_basic_set_free(bset);
678 dim = isl_basic_set_dim(bset, isl_dim_set);
679 if (dim < 0)
680 return isl_basic_set_free(bset);
681 isl_assert(bset->ctx, bset->n_eq <= dim, goto error);
682 if (bset->n_eq == 0)
683 return return_with_identity(bset, T, T2);
685 B = isl_mat_sub_alloc6(bset->ctx, bset->eq, 0, bset->n_eq, 0, 1 + dim);
686 TC = isl_mat_variable_compression(B, T2);
687 if (!TC)
688 goto error;
689 if (TC->n_col == 0) {
690 isl_mat_free(TC);
691 if (T2) {
692 isl_mat_free(*T2);
693 *T2 = NULL;
695 bset = isl_basic_set_set_to_empty(bset);
696 return return_with_identity(bset, T, T2);
699 bset = isl_basic_set_preimage(bset, T ? isl_mat_copy(TC) : TC);
700 if (T)
701 *T = TC;
702 return bset;
703 error:
704 isl_basic_set_free(bset);
705 return NULL;
708 __isl_give isl_basic_set *isl_basic_set_remove_equalities(
709 __isl_take isl_basic_set *bset, __isl_give isl_mat **T,
710 __isl_give isl_mat **T2)
712 if (T)
713 *T = NULL;
714 if (T2)
715 *T2 = NULL;
716 if (isl_basic_set_check_no_params(bset) < 0)
717 return isl_basic_set_free(bset);
718 bset = isl_basic_set_gauss(bset, NULL);
719 if (ISL_F_ISSET(bset, ISL_BASIC_SET_EMPTY))
720 return return_with_identity(bset, T, T2);
721 bset = compress_variables(bset, T, T2);
722 return bset;
725 /* Check if dimension dim belongs to a residue class
726 * i_dim \equiv r mod m
727 * with m != 1 and if so return m in *modulo and r in *residue.
728 * As a special case, when i_dim has a fixed value v, then
729 * *modulo is set to 0 and *residue to v.
731 * If i_dim does not belong to such a residue class, then *modulo
732 * is set to 1 and *residue is set to 0.
734 isl_stat isl_basic_set_dim_residue_class(__isl_keep isl_basic_set *bset,
735 int pos, isl_int *modulo, isl_int *residue)
737 isl_bool fixed;
738 struct isl_ctx *ctx;
739 struct isl_mat *H = NULL, *U = NULL, *C, *H1, *U1;
740 isl_size total;
741 isl_size nparam;
743 if (!bset || !modulo || !residue)
744 return isl_stat_error;
746 fixed = isl_basic_set_plain_dim_is_fixed(bset, pos, residue);
747 if (fixed < 0)
748 return isl_stat_error;
749 if (fixed) {
750 isl_int_set_si(*modulo, 0);
751 return isl_stat_ok;
754 ctx = isl_basic_set_get_ctx(bset);
755 total = isl_basic_set_dim(bset, isl_dim_all);
756 nparam = isl_basic_set_dim(bset, isl_dim_param);
757 if (total < 0 || nparam < 0)
758 return isl_stat_error;
759 H = isl_mat_sub_alloc6(ctx, bset->eq, 0, bset->n_eq, 1, total);
760 H = isl_mat_left_hermite(H, 0, &U, NULL);
761 if (!H)
762 return isl_stat_error;
764 isl_seq_gcd(U->row[nparam + pos]+bset->n_eq,
765 total-bset->n_eq, modulo);
766 if (isl_int_is_zero(*modulo))
767 isl_int_set_si(*modulo, 1);
768 if (isl_int_is_one(*modulo)) {
769 isl_int_set_si(*residue, 0);
770 isl_mat_free(H);
771 isl_mat_free(U);
772 return isl_stat_ok;
775 C = isl_mat_alloc(ctx, 1 + bset->n_eq, 1);
776 if (!C)
777 goto error;
778 isl_int_set_si(C->row[0][0], 1);
779 isl_mat_sub_neg(ctx, C->row + 1, bset->eq, bset->n_eq, 0, 0, 1);
780 H1 = isl_mat_sub_alloc(H, 0, H->n_row, 0, H->n_row);
781 H1 = isl_mat_lin_to_aff(H1);
782 C = isl_mat_inverse_product(H1, C);
783 isl_mat_free(H);
784 U1 = isl_mat_sub_alloc(U, nparam+pos, 1, 0, bset->n_eq);
785 U1 = isl_mat_lin_to_aff(U1);
786 isl_mat_free(U);
787 C = isl_mat_product(U1, C);
788 if (!C)
789 return isl_stat_error;
790 if (!isl_int_is_divisible_by(C->row[1][0], C->row[0][0])) {
791 bset = isl_basic_set_copy(bset);
792 bset = isl_basic_set_set_to_empty(bset);
793 isl_basic_set_free(bset);
794 isl_int_set_si(*modulo, 1);
795 isl_int_set_si(*residue, 0);
796 return isl_stat_ok;
798 isl_int_divexact(*residue, C->row[1][0], C->row[0][0]);
799 isl_int_fdiv_r(*residue, *residue, *modulo);
800 isl_mat_free(C);
801 return isl_stat_ok;
802 error:
803 isl_mat_free(H);
804 isl_mat_free(U);
805 return isl_stat_error;
808 /* Check if dimension dim belongs to a residue class
809 * i_dim \equiv r mod m
810 * with m != 1 and if so return m in *modulo and r in *residue.
811 * As a special case, when i_dim has a fixed value v, then
812 * *modulo is set to 0 and *residue to v.
814 * If i_dim does not belong to such a residue class, then *modulo
815 * is set to 1 and *residue is set to 0.
817 isl_stat isl_set_dim_residue_class(__isl_keep isl_set *set,
818 int pos, isl_int *modulo, isl_int *residue)
820 isl_int m;
821 isl_int r;
822 int i;
824 if (!set || !modulo || !residue)
825 return isl_stat_error;
827 if (set->n == 0) {
828 isl_int_set_si(*modulo, 0);
829 isl_int_set_si(*residue, 0);
830 return isl_stat_ok;
833 if (isl_basic_set_dim_residue_class(set->p[0], pos, modulo, residue)<0)
834 return isl_stat_error;
836 if (set->n == 1)
837 return isl_stat_ok;
839 if (isl_int_is_one(*modulo))
840 return isl_stat_ok;
842 isl_int_init(m);
843 isl_int_init(r);
845 for (i = 1; i < set->n; ++i) {
846 if (isl_basic_set_dim_residue_class(set->p[i], pos, &m, &r) < 0)
847 goto error;
848 isl_int_gcd(*modulo, *modulo, m);
849 isl_int_sub(m, *residue, r);
850 isl_int_gcd(*modulo, *modulo, m);
851 if (!isl_int_is_zero(*modulo))
852 isl_int_fdiv_r(*residue, *residue, *modulo);
853 if (isl_int_is_one(*modulo))
854 break;
857 isl_int_clear(m);
858 isl_int_clear(r);
860 return isl_stat_ok;
861 error:
862 isl_int_clear(m);
863 isl_int_clear(r);
864 return isl_stat_error;
867 /* Check if dimension "dim" belongs to a residue class
868 * i_dim \equiv r mod m
869 * with m != 1 and if so return m in *modulo and r in *residue.
870 * As a special case, when i_dim has a fixed value v, then
871 * *modulo is set to 0 and *residue to v.
873 * If i_dim does not belong to such a residue class, then *modulo
874 * is set to 1 and *residue is set to 0.
876 isl_stat isl_set_dim_residue_class_val(__isl_keep isl_set *set,
877 int pos, __isl_give isl_val **modulo, __isl_give isl_val **residue)
879 *modulo = NULL;
880 *residue = NULL;
881 if (!set)
882 return isl_stat_error;
883 *modulo = isl_val_alloc(isl_set_get_ctx(set));
884 *residue = isl_val_alloc(isl_set_get_ctx(set));
885 if (!*modulo || !*residue)
886 goto error;
887 if (isl_set_dim_residue_class(set, pos,
888 &(*modulo)->n, &(*residue)->n) < 0)
889 goto error;
890 isl_int_set_si((*modulo)->d, 1);
891 isl_int_set_si((*residue)->d, 1);
892 return isl_stat_ok;
893 error:
894 isl_val_free(*modulo);
895 isl_val_free(*residue);
896 return isl_stat_error;