898 lines · c
1/*2 * Copyright 2008-2009 Katholieke Universiteit Leuven3 * Copyright 2010 INRIA Saclay4 *5 * Use of this software is governed by the MIT license6 *7 * Written by Sven Verdoolaege, K.U.Leuven, Departement8 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium9 * and INRIA Saclay - Ile-de-France, Parc Club Orsay Universite,10 * ZAC des vignes, 4 rue Jacques Monod, 91893 Orsay, France11 */12 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>19 20/* Given a set of modulo constraints21 *22 * c + A y = 0 mod d23 *24 * this function computes a particular solution y_025 *26 * The input is given as a matrix B = [ c A ] and a vector d.27 *28 * The output is matrix containing the solution y_0 or29 * a zero-column matrix if the constraints admit no integer solution.30 *31 * The given set of constrains is equivalent to32 *33 * c + A y = -D x34 *35 * with D = diag d and x a fresh set of variables.36 * Reducing both c and A modulo d does not change the37 * 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 * Then40 * [ x ]41 * M [ y ] = - c42 * and so43 * [ x ]44 * [ H 0 ] U^{-1} [ y ] = - c45 * Let46 * [ A ] [ x ]47 * [ B ] = U^{-1} [ y ]48 * then49 * H A + 0 B = -c50 *51 * so B may be chosen arbitrarily, e.g., B = 0, and then52 *53 * [ x ] = [ -c ]54 * U^{-1} [ y ] = [ 0 ]55 * or56 * [ x ] [ -c ]57 * [ y ] = U [ 0 ]58 * specifically,59 *60 * y = U_{2,1} (-c)61 *62 * If any of the coordinates of this y are non-integer63 * then the constraints admit no integer solution and64 * a zero-column matrix is returned.65 */66static __isl_give isl_mat *particular_solution(__isl_keep isl_mat *B,67 __isl_keep isl_vec *d)68{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;76 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]);90 }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]);103 }104 if (i < B->n_row)105 cst = isl_mat_alloc(B->ctx, B->n_row, 0);106 else107 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;114error:115 isl_mat_free(M);116 isl_mat_free(C);117 isl_mat_free(U);118 return NULL;119}120 121/* Compute and return the matrix122 *123 * U_1^{-1} diag(d_1, 1, ..., 1)124 *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 satisfies127 * the single (linear) modulo constraint.128 */129static __isl_take isl_mat *parameter_compression_1(__isl_keep isl_mat *B,130 __isl_keep isl_vec *d)131{132 struct isl_mat *U;133 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;145}146 147/* Compute a common lattice of solutions to the linear modulo148 * constraints specified by B and d.149 * See also the documentation of isl_mat_parameter_compression.150 * We put the matrix151 * 152 * A = [ L_1^{-T} L_2^{-T} ... L_k^{-T} ]153 *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 have156 * 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 have158 * D * L_i^{-T} = U_i^T diag(D/d_i, D, ..., D).159 */160static __isl_give isl_mat *parameter_compression_multi(__isl_keep isl_mat *B,161 __isl_keep isl_vec *d)162{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;168 169 isl_int_init(D);170 171 isl_vec_lcm(d, &D);172 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]);191 }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);205 206 isl_int_clear(D);207 return T;208error:209 isl_mat_free(A);210 isl_mat_free(U);211 isl_int_clear(D);212 return NULL;213}214 215/* Given a set of modulo constraints216 *217 * c + A y = 0 mod d218 *219 * this function returns an affine transformation T,220 *221 * y = T y'222 *223 * that bijectively maps the integer vectors y' to integer224 * vectors y that satisfy the modulo constraints.225 *226 * This function is inspired by Section 2.5.3227 * of B. Meister, "Stating and Manipulating Periodicity in the Polytope228 * Model. Applications to Program Analysis and Optimization".229 * However, the implementation only follows the algorithm of that230 * section for computing a particular solution and not for computing231 * a general homogeneous solution. The latter is incomplete and232 * may remove some valid solutions.233 * Instead, we use an adaptation of the algorithm in Section 7 of234 * B. Meister, S. Verdoolaege, "Polynomial Approximations in the Polytope235 * Model: Bringing the Power of Quasi-Polynomials to the Masses".236 *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 then241 * a matrix with zero columns is returned.242 *243 * We first compute a particular solution y_0 to the given set of244 * modulo constraints in particular_solution. If no such solution245 * exists, then we return a zero-columned transformation matrix.246 * Otherwise, we compute the generic solution to247 *248 * A y = 0 mod d249 *250 * That is we want to compute G such that251 *252 * y = G y''253 *254 * with y'' integer, describes the set of solutions.255 *256 * We first remove the common factors of each row.257 * In particular if gcd(A_i,d_i) != 1, then we divide the whole258 * 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).261 *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 vectors264 * that satisfies this row. Let U_i be the unimodular extension of the265 * row A_i. This unimodular extension exists because gcd(A_i) = 1.266 * The first component of267 *268 * y' = U_i y269 *270 * needs to be a multiple of d_i. Let y' = diag(d_i, 1, ..., 1) y''.271 * Then,272 *273 * y = U_i^{-1} diag(d_i, 1, ..., 1) y''274 *275 * for arbitrary integer vectors y''. That is, y belongs to the lattice276 * 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.278 *279 * If there is more than one row left, we need to compute the intersection280 * of the lattices. That is, we need to compute an L such that281 *282 * L = L_i L_i' for all i283 *284 * with L_i' some integer matrices. Let A be constructed as follows285 *286 * A = [ L_1^{-T} L_2^{-T} ... L_k^{-T} ]287 *288 * and computed the Hermite Normal Form of A = [ H 0 ] U289 * Then,290 *291 * L_i^{-T} = H U_{1,i}292 *293 * or294 *295 * H^{-T} = L_i U_{1,i}^T296 *297 * In other words G = L = H^{-T}.298 * To ensure that G is lower triangular, we compute and use its Hermite299 * normal form.300 *301 * The affine transformation matrix returned is then302 *303 * [ 1 0 ]304 * [ y_0 G ]305 *306 * as any y = y_0 + G y' with y' integer is a solution to the original307 * modulo constraints.308 */309__isl_give isl_mat *isl_mat_parameter_compression(__isl_take isl_mat *B,310 __isl_take isl_vec *d)311{312 int i;313 struct isl_mat *cst = NULL;314 struct isl_mat *T = NULL;315 isl_int D;316 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;329 }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;346 }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);356 }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 else363 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;372error2:373 isl_int_clear(D);374error:375 isl_mat_free(cst);376 isl_mat_free(B);377 isl_vec_free(d);378 return NULL;379}380 381/* Given a set of equalities382 *383 * B(y) + A x = 0 (*)384 *385 * compute and return an affine transformation T,386 *387 * y = T y'388 *389 * that bijectively maps the integer vectors y' to integer390 * vectors y that satisfy the modulo constraints for some value of x.391 *392 * Let [H 0] be the Hermite Normal Form of A, i.e.,393 *394 * A = [H 0] Q395 *396 * Then y is a solution of (*) iff397 *398 * H^-1 B(y) (= - [I 0] Q x)399 *400 * is an integer vector. Let d be the common denominator of H^-1.401 * We impose402 *403 * d H^-1 B(y) = 0 mod d404 *405 * and compute the solution using isl_mat_parameter_compression.406 */407__isl_give isl_mat *isl_mat_parameter_compression_ext(__isl_take isl_mat *B,408 __isl_take isl_mat *A)409{410 isl_ctx *ctx;411 isl_vec *d;412 int n_row, n_col;413 414 if (!A)415 return isl_mat_free(B);416 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);430 431 return isl_mat_parameter_compression(B, d);432}433 434/* Return a compression matrix that indicates that there are no solutions435 * to the original constraints. In particular, return a zero-column436 * matrix with 1 + dim rows. If "T2" is not NULL, then assign *T2437 * the inverse of this matrix. *T2 may already have been assigned438 * matrix, so free it first.439 * "free1", "free2" and "free3" are temporary matrices that are440 * not useful when an empty compression is returned. They are441 * simply freed.442 */443static __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)446{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);453 }454 return isl_mat_alloc(ctx, 1 + dim, 0);455}456 457/* Given a matrix that maps a (possibly) parametric domain to458 * a parametric domain, add in rows that map the "nparam" parameters onto459 * themselves.460 */461static __isl_give isl_mat *insert_parameter_rows(__isl_take isl_mat *mat,462 unsigned nparam)463{464 int i;465 466 if (nparam == 0)467 return mat;468 if (!mat)469 return NULL;470 471 mat = isl_mat_insert_rows(mat, 1, nparam);472 if (!mat)473 return NULL;474 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]);478 }479 480 return mat;481}482 483/* Given a set of equalities484 *485 * -C(y) + M x = 0486 *487 * this function computes a unimodular transformation from a lower-dimensional488 * 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 original490 * space that satisfy the equalities.491 *492 * The input is given as a matrix B = [ -C M ] and the output is a493 * matrix that maps [1 x'] to [1 x].494 * The number of equality constraints in B is assumed to be smaller than495 * 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'].499 *500 * First compute the (left) Hermite normal form of M,501 *502 * M [U1 U2] = M U = H = [H1 0]503 * or504 * M = H Q = [H1 0] [Q1]505 * [Q2]506 *507 * with U, Q unimodular, Q = U^{-1} (and H lower triangular).508 * Define the transformed variables as509 *510 * x = [U1 U2] [ x1' ] = [U1 U2] [Q1] x511 * [ x2' ] [Q2]512 *513 * The equalities then become514 *515 * -C(y) + H1 x1' = 0 or x1' = H1^{-1} C(y) = C'(y)516 *517 * If the denominator of the constant term does not divide the518 * the common denominator of the coefficients of y, then every519 * integer point is mapped to a non-integer point and then the original set520 * has no integer solutions (since the x' are a unimodular transformation521 * of the x). In this case, a zero-column matrix is returned.522 * Otherwise, the transformation is given by523 *524 * x = U1 H1^{-1} C(y) + U2 x2'525 *526 * The inverse transformation is simply527 *528 * x2' = Q2 x529 */530__isl_give isl_mat *isl_mat_final_variable_compression(__isl_take isl_mat *B,531 int first, __isl_give isl_mat **T2)532{533 int i, n;534 isl_ctx *ctx;535 isl_mat *H = NULL, *C, *H1, *U = NULL, *U1, *U2;536 unsigned dim;537 538 if (T2)539 *T2 = NULL;540 if (!B)541 goto error;542 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;558 }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;573 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;580 }581 isl_int_clear(g);582 583 if (i < B->n_row)584 return empty_compression(ctx, dim, T2, B, C, U);585 C = isl_mat_normalize(C);586 }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);595 596 isl_mat_free(B);597 598 return C;599error: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;606 }607 return NULL;608}609 610/* Given a set of equalities611 *612 * M x - c = 0613 *614 * this function computes a unimodular transformation from a lower-dimensional615 * 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 original617 * space that satisfy the equalities.618 *619 * The input is given as a matrix B = [ -c M ] and the output is a620 * matrix that maps [1 x'] to [1 x].621 * The number of equality constraints in B is assumed to be smaller than622 * 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'].624 */625__isl_give isl_mat *isl_mat_variable_compression(__isl_take isl_mat *B,626 __isl_give isl_mat **T2)627{628 return isl_mat_final_variable_compression(B, 0, T2);629}630 631/* Return "bset" and set *T and *T2 to the identity transformation632 * on "bset" (provided T and T2 are not NULL).633 */634static __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)637{638 isl_size dim;639 isl_mat *id;640 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;646 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);653 654 return bset;655}656 657/* Use the n equalities of bset to unimodularly transform the658 * variables x such that n transformed variables x1' have a constant value659 * and rewrite the constraints of bset in terms of the remaining660 * transformed variables x2'. The matrix pointed to by T maps661 * the new variables x2' back to the original variables x, while T2662 * maps the original variables to the new variables.663 */664static __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)667{668 struct isl_mat *B, *TC;669 isl_size dim;670 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);684 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;694 }695 bset = isl_basic_set_set_to_empty(bset);696 return return_with_identity(bset, T, T2);697 }698 699 bset = isl_basic_set_preimage(bset, T ? isl_mat_copy(TC) : TC);700 if (T)701 *T = TC;702 return bset;703error:704 isl_basic_set_free(bset);705 return NULL;706}707 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)711{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;723}724 725/* Check if dimension dim belongs to a residue class726 * i_dim \equiv r mod m727 * 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, then729 * *modulo is set to 0 and *residue to v.730 *731 * If i_dim does not belong to such a residue class, then *modulo732 * is set to 1 and *residue is set to 0.733 */734isl_stat isl_basic_set_dim_residue_class(__isl_keep isl_basic_set *bset,735 int pos, isl_int *modulo, isl_int *residue)736{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;742 743 if (!bset || !modulo || !residue)744 return isl_stat_error;745 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;752 }753 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;763 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;773 }774 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;797 }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;802error:803 isl_mat_free(H);804 isl_mat_free(U);805 return isl_stat_error;806}807 808/* Check if dimension dim belongs to a residue class809 * i_dim \equiv r mod m810 * 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, then812 * *modulo is set to 0 and *residue to v.813 *814 * If i_dim does not belong to such a residue class, then *modulo815 * is set to 1 and *residue is set to 0.816 */817isl_stat isl_set_dim_residue_class(__isl_keep isl_set *set,818 int pos, isl_int *modulo, isl_int *residue)819{820 isl_int m;821 isl_int r;822 int i;823 824 if (!set || !modulo || !residue)825 return isl_stat_error;826 827 if (set->n == 0) {828 isl_int_set_si(*modulo, 0);829 isl_int_set_si(*residue, 0);830 return isl_stat_ok;831 }832 833 if (isl_basic_set_dim_residue_class(set->p[0], pos, modulo, residue)<0)834 return isl_stat_error;835 836 if (set->n == 1)837 return isl_stat_ok;838 839 if (isl_int_is_one(*modulo))840 return isl_stat_ok;841 842 isl_int_init(m);843 isl_int_init(r);844 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;855 }856 857 isl_int_clear(m);858 isl_int_clear(r);859 860 return isl_stat_ok;861error:862 isl_int_clear(m);863 isl_int_clear(r);864 return isl_stat_error;865}866 867/* Check if dimension "dim" belongs to a residue class868 * i_dim \equiv r mod m869 * 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, then871 * *modulo is set to 0 and *residue to v.872 *873 * If i_dim does not belong to such a residue class, then *modulo874 * is set to 1 and *residue is set to 0.875 */876isl_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)878{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;893error:894 isl_val_free(*modulo);895 isl_val_free(*residue);896 return isl_stat_error;897}898