## Coverage Report

#### Created: 2019-07-24 05:18

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