1/* 2 * Copyright 2008-2009 Katholieke Universiteit Leuven 3 * Copyright 2010 INRIA Saclay 4 * 5 * Use of this software is governed by the MIT license 6 * 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 11 */ 12 13#include <isl_mat_private.h> 14#include <isl/seq.h> 15#include "isl_map_private.h" 16#include "isl_equalities.h" 17#include <isl_val_private.h> 18 19/* Given a set of modulo constraints 20 * 21 * c + A y = 0 mod d 22 * 23 * this function computes a particular solution y_0 24 * 25 * The input is given as a matrix B = [ c A ] and a vector d. 26 * 27 * The output is matrix containing the solution y_0 or 28 * a zero-column matrix if the constraints admit no integer solution. 29 * 30 * The given set of constrains is equivalent to 31 * 32 * c + A y = -D x 33 * 34 * with D = diag d and x a fresh set of variables. 35 * Reducing both c and A modulo d does not change the 36 * value of y in the solution and may lead to smaller coefficients. 37 * Let M = [ D A ] and [ H 0 ] = M U, the Hermite normal form of M. 38 * Then 39 * [ x ] 40 * M [ y ] = - c 41 * and so 42 * [ x ] 43 * [ H 0 ] U^{-1} [ y ] = - c 44 * Let 45 * [ A ] [ x ] 46 * [ B ] = U^{-1} [ y ] 47 * then 48 * H A + 0 B = -c 49 * 50 * so B may be chosen arbitrarily, e.g., B = 0, and then 51 * 52 * [ x ] = [ -c ] 53 * U^{-1} [ y ] = [ 0 ] 54 * or 55 * [ x ] [ -c ] 56 * [ y ] = U [ 0 ] 57 * specifically, 58 * 59 * y = U_{2,1} (-c) 60 * 61 * If any of the coordinates of this y are non-integer 62 * then the constraints admit no integer solution and 63 * a zero-column matrix is returned. 64 */ 65static struct isl_mat *particular_solution(struct isl_mat *B, struct isl_vec *d) 66{ 67 int i, j; 68 struct isl_mat *M = NULL; 69 struct isl_mat *C = NULL; 70 struct isl_mat *U = NULL; 71 struct isl_mat *H = NULL; 72 struct isl_mat *cst = NULL; 73 struct isl_mat *T = NULL; 74 75 M = isl_mat_alloc(B->ctx, B->n_row, B->n_row + B->n_col - 1); 76 C = isl_mat_alloc(B->ctx, 1 + B->n_row, 1); 77 if (!M || !C) 78 goto error; 79 isl_int_set_si(C->row[0][0], 1); 80 for (i = 0; i < B->n_row; ++i) { 81 isl_seq_clr(M->row[i], B->n_row); 82 isl_int_set(M->row[i][i], d->block.data[i]); 83 isl_int_neg(C->row[1 + i][0], B->row[i][0]); 84 isl_int_fdiv_r(C->row[1+i][0], C->row[1+i][0], M->row[i][i]); 85 for (j = 0; j < B->n_col - 1; ++j) 86 isl_int_fdiv_r(M->row[i][B->n_row + j], 87 B->row[i][1 + j], M->row[i][i]); 88 } 89 M = isl_mat_left_hermite(M, 0, &U, NULL); 90 if (!M || !U) 91 goto error; 92 H = isl_mat_sub_alloc(M, 0, B->n_row, 0, B->n_row); 93 H = isl_mat_lin_to_aff(H); 94 C = isl_mat_inverse_product(H, C); 95 if (!C) 96 goto error; 97 for (i = 0; i < B->n_row; ++i) { 98 if (!isl_int_is_divisible_by(C->row[1+i][0], C->row[0][0])) 99 break; 100 isl_int_divexact(C->row[1+i][0], C->row[1+i][0], C->row[0][0]); 101 } 102 if (i < B->n_row) 103 cst = isl_mat_alloc(B->ctx, B->n_row, 0); 104 else 105 cst = isl_mat_sub_alloc(C, 1, B->n_row, 0, 1); 106 T = isl_mat_sub_alloc(U, B->n_row, B->n_col - 1, 0, B->n_row); 107 cst = isl_mat_product(T, cst); 108 isl_mat_free(M); 109 isl_mat_free(C); 110 isl_mat_free(U); 111 return cst; 112error: 113 isl_mat_free(M); 114 isl_mat_free(C); 115 isl_mat_free(U); 116 return NULL; 117} 118 119/* Compute and return the matrix 120 * 121 * U_1^{-1} diag(d_1, 1, ..., 1) 122 * 123 * with U_1 the unimodular completion of the first (and only) row of B. 124 * The columns of this matrix generate the lattice that satisfies 125 * the single (linear) modulo constraint. 126 */ 127static struct isl_mat *parameter_compression_1( 128 struct isl_mat *B, struct isl_vec *d) 129{ 130 struct isl_mat *U; 131 132 U = isl_mat_alloc(B->ctx, B->n_col - 1, B->n_col - 1); 133 if (!U) 134 return NULL; 135 isl_seq_cpy(U->row[0], B->row[0] + 1, B->n_col - 1); 136 U = isl_mat_unimodular_complete(U, 1); 137 U = isl_mat_right_inverse(U); 138 if (!U) 139 return NULL; 140 isl_mat_col_mul(U, 0, d->block.data[0], 0); 141 U = isl_mat_lin_to_aff(U); 142 return U; 143} 144 145/* Compute a common lattice of solutions to the linear modulo 146 * constraints specified by B and d. 147 * See also the documentation of isl_mat_parameter_compression. 148 * We put the matrix 149 * 150 * A = [ L_1^{-T} L_2^{-T} ... L_k^{-T} ] 151 * 152 * on a common denominator. This denominator D is the lcm of modulos d. 153 * Since L_i = U_i^{-1} diag(d_i, 1, ... 1), we have 154 * L_i^{-T} = U_i^T diag(d_i, 1, ... 1)^{-T} = U_i^T diag(1/d_i, 1, ..., 1). 155 * Putting this on the common denominator, we have 156 * D * L_i^{-T} = U_i^T diag(D/d_i, D, ..., D). 157 */ 158static struct isl_mat *parameter_compression_multi( 159 struct isl_mat *B, struct isl_vec *d) 160{ 161 int i, j, k; 162 isl_int D; 163 struct isl_mat *A = NULL, *U = NULL; 164 struct isl_mat *T; 165 unsigned size; 166 167 isl_int_init(D); 168 169 isl_vec_lcm(d, &D); 170 171 size = B->n_col - 1; 172 A = isl_mat_alloc(B->ctx, size, B->n_row * size); 173 U = isl_mat_alloc(B->ctx, size, size); 174 if (!U || !A) 175 goto error; 176 for (i = 0; i < B->n_row; ++i) { 177 isl_seq_cpy(U->row[0], B->row[i] + 1, size); 178 U = isl_mat_unimodular_complete(U, 1); 179 if (!U) 180 goto error; 181 isl_int_divexact(D, D, d->block.data[i]); 182 for (k = 0; k < U->n_col; ++k) 183 isl_int_mul(A->row[k][i*size+0], D, U->row[0][k]); 184 isl_int_mul(D, D, d->block.data[i]); 185 for (j = 1; j < U->n_row; ++j) 186 for (k = 0; k < U->n_col; ++k) 187 isl_int_mul(A->row[k][i*size+j], 188 D, U->row[j][k]); 189 } 190 A = isl_mat_left_hermite(A, 0, NULL, NULL); 191 T = isl_mat_sub_alloc(A, 0, A->n_row, 0, A->n_row); 192 T = isl_mat_lin_to_aff(T); 193 if (!T) 194 goto error; 195 isl_int_set(T->row[0][0], D); 196 T = isl_mat_right_inverse(T); 197 if (!T) 198 goto error; 199 isl_assert(T->ctx, isl_int_is_one(T->row[0][0]), goto error); 200 T = isl_mat_transpose(T); 201 isl_mat_free(A); 202 isl_mat_free(U); 203 204 isl_int_clear(D); 205 return T; 206error: 207 isl_mat_free(A); 208 isl_mat_free(U); 209 isl_int_clear(D); 210 return NULL; 211} 212 213/* Given a set of modulo constraints 214 * 215 * c + A y = 0 mod d 216 * 217 * this function returns an affine transformation T, 218 * 219 * y = T y' 220 * 221 * that bijectively maps the integer vectors y' to integer 222 * vectors y that satisfy the modulo constraints. 223 * 224 * This function is inspired by Section 2.5.3 225 * of B. Meister, "Stating and Manipulating Periodicity in the Polytope 226 * Model. Applications to Program Analysis and Optimization". 227 * However, the implementation only follows the algorithm of that 228 * section for computing a particular solution and not for computing 229 * a general homogeneous solution. The latter is incomplete and 230 * may remove some valid solutions. 231 * Instead, we use an adaptation of the algorithm in Section 7 of 232 * B. Meister, S. Verdoolaege, "Polynomial Approximations in the Polytope 233 * Model: Bringing the Power of Quasi-Polynomials to the Masses". 234 * 235 * The input is given as a matrix B = [ c A ] and a vector d. 236 * Each element of the vector d corresponds to a row in B. 237 * The output is a lower triangular matrix. 238 * If no integer vector y satisfies the given constraints then 239 * a matrix with zero columns is returned. 240 * 241 * We first compute a particular solution y_0 to the given set of 242 * modulo constraints in particular_solution. If no such solution 243 * exists, then we return a zero-columned transformation matrix. 244 * Otherwise, we compute the generic solution to 245 * 246 * A y = 0 mod d 247 * 248 * That is we want to compute G such that 249 * 250 * y = G y'' 251 * 252 * with y'' integer, describes the set of solutions. 253 * 254 * We first remove the common factors of each row. 255 * In particular if gcd(A_i,d_i) != 1, then we divide the whole 256 * row i (including d_i) by this common factor. If afterwards gcd(A_i) != 1, 257 * then we divide this row of A by the common factor, unless gcd(A_i) = 0. 258 * In the later case, we simply drop the row (in both A and d). 259 * 260 * If there are no rows left in A, then G is the identity matrix. Otherwise, 261 * for each row i, we now determine the lattice of integer vectors 262 * that satisfies this row. Let U_i be the unimodular extension of the 263 * row A_i. This unimodular extension exists because gcd(A_i) = 1. 264 * The first component of 265 * 266 * y' = U_i y 267 * 268 * needs to be a multiple of d_i. Let y' = diag(d_i, 1, ..., 1) y''. 269 * Then, 270 * 271 * y = U_i^{-1} diag(d_i, 1, ..., 1) y'' 272 * 273 * for arbitrary integer vectors y''. That is, y belongs to the lattice 274 * generated by the columns of L_i = U_i^{-1} diag(d_i, 1, ..., 1). 275 * If there is only one row, then G = L_1. 276 * 277 * If there is more than one row left, we need to compute the intersection 278 * of the lattices. That is, we need to compute an L such that 279 * 280 * L = L_i L_i' for all i 281 * 282 * with L_i' some integer matrices. Let A be constructed as follows 283 * 284 * A = [ L_1^{-T} L_2^{-T} ... L_k^{-T} ] 285 * 286 * and computed the Hermite Normal Form of A = [ H 0 ] U 287 * Then, 288 * 289 * L_i^{-T} = H U_{1,i} 290 * 291 * or 292 * 293 * H^{-T} = L_i U_{1,i}^T 294 * 295 * In other words G = L = H^{-T}. 296 * To ensure that G is lower triangular, we compute and use its Hermite 297 * normal form. 298 * 299 * The affine transformation matrix returned is then 300 * 301 * [ 1 0 ] 302 * [ y_0 G ] 303 * 304 * as any y = y_0 + G y' with y' integer is a solution to the original 305 * modulo constraints. 306 */ 307struct isl_mat *isl_mat_parameter_compression( 308 struct isl_mat *B, struct isl_vec *d) 309{ 310 int i; 311 struct isl_mat *cst = NULL; 312 struct isl_mat *T = NULL; 313 isl_int D; 314 315 if (!B || !d) 316 goto error; 317 isl_assert(B->ctx, B->n_row == d->size, goto error); 318 cst = particular_solution(B, d); 319 if (!cst) 320 goto error; 321 if (cst->n_col == 0) { 322 T = isl_mat_alloc(B->ctx, B->n_col, 0); 323 isl_mat_free(cst); 324 isl_mat_free(B); 325 isl_vec_free(d); 326 return T; 327 } 328 isl_int_init(D); 329 /* Replace a*g*row = 0 mod g*m by row = 0 mod m */ 330 for (i = 0; i < B->n_row; ++i) { 331 isl_seq_gcd(B->row[i] + 1, B->n_col - 1, &D); 332 if (isl_int_is_one(D)) 333 continue; 334 if (isl_int_is_zero(D)) { 335 B = isl_mat_drop_rows(B, i, 1); 336 d = isl_vec_cow(d); 337 if (!B || !d) 338 goto error2; 339 isl_seq_cpy(d->block.data+i, d->block.data+i+1, 340 d->size - (i+1)); 341 d->size--; 342 i--; 343 continue; 344 } 345 B = isl_mat_cow(B); 346 if (!B) 347 goto error2; 348 isl_seq_scale_down(B->row[i] + 1, B->row[i] + 1, D, B->n_col-1); 349 isl_int_gcd(D, D, d->block.data[i]); 350 d = isl_vec_cow(d); 351 if (!d) 352 goto error2; 353 isl_int_divexact(d->block.data[i], d->block.data[i], D); 354 } 355 isl_int_clear(D); 356 if (B->n_row == 0) 357 T = isl_mat_identity(B->ctx, B->n_col); 358 else if (B->n_row == 1) 359 T = parameter_compression_1(B, d); 360 else 361 T = parameter_compression_multi(B, d); 362 T = isl_mat_left_hermite(T, 0, NULL, NULL); 363 if (!T) 364 goto error; 365 isl_mat_sub_copy(T->ctx, T->row + 1, cst->row, cst->n_row, 0, 0, 1); 366 isl_mat_free(cst); 367 isl_mat_free(B); 368 isl_vec_free(d); 369 return T; 370error2: 371 isl_int_clear(D); 372error: 373 isl_mat_free(cst); 374 isl_mat_free(B); 375 isl_vec_free(d); 376 return NULL; 377} 378 379/* Given a set of equalities 380 * 381 * B(y) + A x = 0 (*) 382 * 383 * compute and return an affine transformation T, 384 * 385 * y = T y' 386 * 387 * that bijectively maps the integer vectors y' to integer 388 * vectors y that satisfy the modulo constraints for some value of x. 389 * 390 * Let [H 0] be the Hermite Normal Form of A, i.e., 391 * 392 * A = [H 0] Q 393 * 394 * Then y is a solution of (*) iff 395 * 396 * H^-1 B(y) (= - [I 0] Q x) 397 * 398 * is an integer vector. Let d be the common denominator of H^-1. 399 * We impose 400 * 401 * d H^-1 B(y) = 0 mod d 402 * 403 * and compute the solution using isl_mat_parameter_compression. 404 */ 405__isl_give isl_mat *isl_mat_parameter_compression_ext(__isl_take isl_mat *B, 406 __isl_take isl_mat *A) 407{ 408 isl_ctx *ctx; 409 isl_vec *d; 410 int n_row, n_col; 411 412 if (!A) 413 return isl_mat_free(B); 414 415 ctx = isl_mat_get_ctx(A); 416 n_row = A->n_row; 417 n_col = A->n_col; 418 A = isl_mat_left_hermite(A, 0, NULL, NULL); 419 A = isl_mat_drop_cols(A, n_row, n_col - n_row); 420 A = isl_mat_lin_to_aff(A); 421 A = isl_mat_right_inverse(A); 422 d = isl_vec_alloc(ctx, n_row); 423 if (A) 424 d = isl_vec_set(d, A->row[0][0]); 425 A = isl_mat_drop_rows(A, 0, 1); 426 A = isl_mat_drop_cols(A, 0, 1); 427 B = isl_mat_product(A, B); 428 429 return isl_mat_parameter_compression(B, d); 430} 431 432/* Given a set of equalities 433 * 434 * M x - c = 0 435 * 436 * this function computes a unimodular transformation from a lower-dimensional 437 * space to the original space that bijectively maps the integer points x' 438 * in the lower-dimensional space to the integer points x in the original 439 * space that satisfy the equalities. 440 * 441 * The input is given as a matrix B = [ -c M ] and the output is a 442 * matrix that maps [1 x'] to [1 x]. 443 * If T2 is not NULL, then *T2 is set to a matrix mapping [1 x] to [1 x']. 444 * 445 * First compute the (left) Hermite normal form of M, 446 * 447 * M [U1 U2] = M U = H = [H1 0] 448 * or 449 * M = H Q = [H1 0] [Q1] 450 * [Q2] 451 * 452 * with U, Q unimodular, Q = U^{-1} (and H lower triangular). 453 * Define the transformed variables as 454 * 455 * x = [U1 U2] [ x1' ] = [U1 U2] [Q1] x 456 * [ x2' ] [Q2] 457 * 458 * The equalities then become 459 * 460 * H1 x1' - c = 0 or x1' = H1^{-1} c = c' 461 * 462 * If any of the c' is non-integer, then the original set has no 463 * integer solutions (since the x' are a unimodular transformation 464 * of the x) and a zero-column matrix is returned. 465 * Otherwise, the transformation is given by 466 * 467 * x = U1 H1^{-1} c + U2 x2' 468 * 469 * The inverse transformation is simply 470 * 471 * x2' = Q2 x 472 */ 473__isl_give isl_mat *isl_mat_variable_compression(__isl_take isl_mat *B, 474 __isl_give isl_mat **T2) 475{ 476 int i; 477 struct isl_mat *H = NULL, *C = NULL, *H1, *U = NULL, *U1, *U2, *TC; 478 unsigned dim; 479 480 if (T2) 481 *T2 = NULL; 482 if (!B) 483 goto error; 484 485 dim = B->n_col - 1; 486 H = isl_mat_sub_alloc(B, 0, B->n_row, 1, dim); 487 H = isl_mat_left_hermite(H, 0, &U, T2); 488 if (!H || !U || (T2 && !*T2)) 489 goto error; 490 if (T2) { 491 *T2 = isl_mat_drop_rows(*T2, 0, B->n_row); 492 *T2 = isl_mat_lin_to_aff(*T2); 493 if (!*T2) 494 goto error; 495 } 496 C = isl_mat_alloc(B->ctx, 1+B->n_row, 1); 497 if (!C) 498 goto error; 499 isl_int_set_si(C->row[0][0], 1); 500 isl_mat_sub_neg(C->ctx, C->row+1, B->row, B->n_row, 0, 0, 1); 501 H1 = isl_mat_sub_alloc(H, 0, H->n_row, 0, H->n_row); 502 H1 = isl_mat_lin_to_aff(H1); 503 TC = isl_mat_inverse_product(H1, C); 504 if (!TC) 505 goto error; 506 isl_mat_free(H); 507 if (!isl_int_is_one(TC->row[0][0])) { 508 for (i = 0; i < B->n_row; ++i) { 509 if (!isl_int_is_divisible_by(TC->row[1+i][0], TC->row[0][0])) { 510 struct isl_ctx *ctx = B->ctx; 511 isl_mat_free(B); 512 isl_mat_free(TC); 513 isl_mat_free(U); 514 if (T2) { 515 isl_mat_free(*T2); 516 *T2 = NULL; 517 } 518 return isl_mat_alloc(ctx, 1 + dim, 0); 519 } 520 isl_seq_scale_down(TC->row[1+i], TC->row[1+i], TC->row[0][0], 1); 521 } 522 isl_int_set_si(TC->row[0][0], 1); 523 } 524 U1 = isl_mat_sub_alloc(U, 0, U->n_row, 0, B->n_row); 525 U1 = isl_mat_lin_to_aff(U1); 526 U2 = isl_mat_sub_alloc(U, 0, U->n_row, B->n_row, U->n_row - B->n_row); 527 U2 = isl_mat_lin_to_aff(U2); 528 isl_mat_free(U); 529 TC = isl_mat_product(U1, TC); 530 TC = isl_mat_aff_direct_sum(TC, U2); 531 532 isl_mat_free(B); 533 534 return TC; 535error: 536 isl_mat_free(B); 537 isl_mat_free(H); 538 isl_mat_free(U); 539 if (T2) { 540 isl_mat_free(*T2); 541 *T2 = NULL; 542 } 543 return NULL; 544} 545 546/* Use the n equalities of bset to unimodularly transform the 547 * variables x such that n transformed variables x1' have a constant value 548 * and rewrite the constraints of bset in terms of the remaining 549 * transformed variables x2'. The matrix pointed to by T maps 550 * the new variables x2' back to the original variables x, while T2 551 * maps the original variables to the new variables. 552 */ 553static struct isl_basic_set *compress_variables( 554 struct isl_basic_set *bset, struct isl_mat **T, struct isl_mat **T2) 555{ 556 struct isl_mat *B, *TC; 557 unsigned dim; 558 559 if (T) 560 *T = NULL; 561 if (T2) 562 *T2 = NULL; 563 if (!bset) 564 goto error; 565 isl_assert(bset->ctx, isl_basic_set_n_param(bset) == 0, goto error); 566 isl_assert(bset->ctx, bset->n_div == 0, goto error); 567 dim = isl_basic_set_n_dim(bset); 568 isl_assert(bset->ctx, bset->n_eq <= dim, goto error); 569 if (bset->n_eq == 0) 570 return bset; 571 572 B = isl_mat_sub_alloc6(bset->ctx, bset->eq, 0, bset->n_eq, 0, 1 + dim); 573 TC = isl_mat_variable_compression(B, T2); 574 if (!TC) 575 goto error; 576 if (TC->n_col == 0) { 577 isl_mat_free(TC); 578 if (T2) { 579 isl_mat_free(*T2); 580 *T2 = NULL; 581 } 582 return isl_basic_set_set_to_empty(bset); 583 } 584 585 bset = isl_basic_set_preimage(bset, T ? isl_mat_copy(TC) : TC); 586 if (T) 587 *T = TC; 588 return bset; 589error: 590 isl_basic_set_free(bset); 591 return NULL; 592} 593 594struct isl_basic_set *isl_basic_set_remove_equalities( 595 struct isl_basic_set *bset, struct isl_mat **T, struct isl_mat **T2) 596{ 597 if (T) 598 *T = NULL; 599 if (T2) 600 *T2 = NULL; 601 if (!bset) 602 return NULL; 603 isl_assert(bset->ctx, isl_basic_set_n_param(bset) == 0, goto error); 604 bset = isl_basic_set_gauss(bset, NULL); 605 if (ISL_F_ISSET(bset, ISL_BASIC_SET_EMPTY)) 606 return bset; 607 bset = compress_variables(bset, T, T2); 608 return bset; 609error: 610 isl_basic_set_free(bset); 611 *T = NULL; 612 return NULL; 613} 614 615/* Check if dimension dim belongs to a residue class 616 * i_dim \equiv r mod m 617 * with m != 1 and if so return m in *modulo and r in *residue. 618 * As a special case, when i_dim has a fixed value v, then 619 * *modulo is set to 0 and *residue to v. 620 * 621 * If i_dim does not belong to such a residue class, then *modulo 622 * is set to 1 and *residue is set to 0. 623 */ 624int isl_basic_set_dim_residue_class(struct isl_basic_set *bset, 625 int pos, isl_int *modulo, isl_int *residue) 626{ 627 struct isl_ctx *ctx; 628 struct isl_mat *H = NULL, *U = NULL, *C, *H1, *U1; 629 unsigned total; 630 unsigned nparam; 631 632 if (!bset || !modulo || !residue) 633 return -1; 634 635 if (isl_basic_set_plain_dim_is_fixed(bset, pos, residue)) { 636 isl_int_set_si(*modulo, 0); 637 return 0; 638 } 639 640 ctx = bset->ctx; 641 total = isl_basic_set_total_dim(bset); 642 nparam = isl_basic_set_n_param(bset); 643 H = isl_mat_sub_alloc6(bset->ctx, bset->eq, 0, bset->n_eq, 1, total); 644 H = isl_mat_left_hermite(H, 0, &U, NULL); 645 if (!H) 646 return -1; 647 648 isl_seq_gcd(U->row[nparam + pos]+bset->n_eq, 649 total-bset->n_eq, modulo); 650 if (isl_int_is_zero(*modulo)) 651 isl_int_set_si(*modulo, 1); 652 if (isl_int_is_one(*modulo)) { 653 isl_int_set_si(*residue, 0); 654 isl_mat_free(H); 655 isl_mat_free(U); 656 return 0; 657 } 658 659 C = isl_mat_alloc(bset->ctx, 1+bset->n_eq, 1); 660 if (!C) 661 goto error; 662 isl_int_set_si(C->row[0][0], 1); 663 isl_mat_sub_neg(C->ctx, C->row+1, bset->eq, bset->n_eq, 0, 0, 1); 664 H1 = isl_mat_sub_alloc(H, 0, H->n_row, 0, H->n_row); 665 H1 = isl_mat_lin_to_aff(H1); 666 C = isl_mat_inverse_product(H1, C); 667 isl_mat_free(H); 668 U1 = isl_mat_sub_alloc(U, nparam+pos, 1, 0, bset->n_eq); 669 U1 = isl_mat_lin_to_aff(U1); 670 isl_mat_free(U); 671 C = isl_mat_product(U1, C); 672 if (!C) 673 goto error; 674 if (!isl_int_is_divisible_by(C->row[1][0], C->row[0][0])) { 675 bset = isl_basic_set_copy(bset); 676 bset = isl_basic_set_set_to_empty(bset); 677 isl_basic_set_free(bset); 678 isl_int_set_si(*modulo, 1); 679 isl_int_set_si(*residue, 0); 680 return 0; 681 } 682 isl_int_divexact(*residue, C->row[1][0], C->row[0][0]); 683 isl_int_fdiv_r(*residue, *residue, *modulo); 684 isl_mat_free(C); 685 return 0; 686error: 687 isl_mat_free(H); 688 isl_mat_free(U); 689 return -1; 690} 691 692/* Check if dimension dim belongs to a residue class 693 * i_dim \equiv r mod m 694 * with m != 1 and if so return m in *modulo and r in *residue. 695 * As a special case, when i_dim has a fixed value v, then 696 * *modulo is set to 0 and *residue to v. 697 * 698 * If i_dim does not belong to such a residue class, then *modulo 699 * is set to 1 and *residue is set to 0. 700 */ 701int isl_set_dim_residue_class(struct isl_set *set, 702 int pos, isl_int *modulo, isl_int *residue) 703{ 704 isl_int m; 705 isl_int r; 706 int i; 707 708 if (!set || !modulo || !residue) 709 return -1; 710 711 if (set->n == 0) { 712 isl_int_set_si(*modulo, 0); 713 isl_int_set_si(*residue, 0); 714 return 0; 715 } 716 717 if (isl_basic_set_dim_residue_class(set->p[0], pos, modulo, residue)<0) 718 return -1; 719 720 if (set->n == 1) 721 return 0; 722 723 if (isl_int_is_one(*modulo)) 724 return 0; 725 726 isl_int_init(m); 727 isl_int_init(r); 728 729 for (i = 1; i < set->n; ++i) { 730 if (isl_basic_set_dim_residue_class(set->p[i], pos, &m, &r) < 0) 731 goto error; 732 isl_int_gcd(*modulo, *modulo, m); 733 isl_int_sub(m, *residue, r); 734 isl_int_gcd(*modulo, *modulo, m); 735 if (!isl_int_is_zero(*modulo)) 736 isl_int_fdiv_r(*residue, *residue, *modulo); 737 if (isl_int_is_one(*modulo)) 738 break; 739 } 740 741 isl_int_clear(m); 742 isl_int_clear(r); 743 744 return 0; 745error: 746 isl_int_clear(m); 747 isl_int_clear(r); 748 return -1; 749} 750 751/* Check if dimension "dim" belongs to a residue class 752 * i_dim \equiv r mod m 753 * with m != 1 and if so return m in *modulo and r in *residue. 754 * As a special case, when i_dim has a fixed value v, then 755 * *modulo is set to 0 and *residue to v. 756 * 757 * If i_dim does not belong to such a residue class, then *modulo 758 * is set to 1 and *residue is set to 0. 759 */ 760int isl_set_dim_residue_class_val(__isl_keep isl_set *set, 761 int pos, __isl_give isl_val **modulo, __isl_give isl_val **residue) 762{ 763 *modulo = NULL; 764 *residue = NULL; 765 if (!set) 766 return -1; 767 *modulo = isl_val_alloc(isl_set_get_ctx(set)); 768 *residue = isl_val_alloc(isl_set_get_ctx(set)); 769 if (!*modulo || !*residue) 770 goto error; 771 if (isl_set_dim_residue_class(set, pos, 772 &(*modulo)->n, &(*residue)->n) < 0) 773 goto error; 774 isl_int_set_si((*modulo)->d, 1); 775 isl_int_set_si((*residue)->d, 1); 776 return 0; 777error: 778 isl_val_free(*modulo); 779 isl_val_free(*residue); 780 return -1; 781} 782