Coverage Report

Created: 2017-10-03 07:32

/Users/buildslave/jenkins/sharedspace/clang-stage2-coverage-R@2/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 14.8k { 68 14.8k   int i, j; 69 14.8k   struct isl_mat *M = NULL; 70 14.8k   struct isl_mat *C = NULL; 71 14.8k   struct isl_mat *U = NULL; 72 14.8k   struct isl_mat *H = NULL; 73 14.8k   struct isl_mat *cst = NULL; 74 14.8k   struct isl_mat *T = NULL; 75 14.8k 76 14.8k   M = isl_mat_alloc(B->ctx, B->n_row, B->n_row + B->n_col - 1); 77 14.8k   C = isl_mat_alloc(B->ctx, 1 + B->n_row, 1); 78 14.8k   if (!M || 14.8k!C14.8k) 79 0     goto error; 80 14.8k   isl_int_set_si14.8k(C->row[0][0], 1); 81 33.3k   for (i = 0; i < B->n_row33.3k; ++i18.5k) { 82 18.5k     isl_seq_clr(M->row[i], B->n_row); 83 18.5k     isl_int_set(M->row[i][i], d->block.data[i]); 84 18.5k     isl_int_neg(C->row[1 + i][0], B->row[i][0]); 85 18.5k     isl_int_fdiv_r(C->row[1+i][0], C->row[1+i][0], M->row[i][i]); 86 85.2k     for (j = 0; j < B->n_col - 185.2k; ++j66.6k) 87 66.6k       isl_int_fdiv_r(M->row[i][B->n_row + j], 88 18.5k           B->row[i][1 + j], M->row[i][i]); 89 18.5k   } 90 14.8k   M = isl_mat_left_hermite(M, 0, &U, NULL); 91 14.8k   if (!M || 14.8k!U14.8k) 92 0     goto error; 93 14.8k   H = isl_mat_sub_alloc(M, 0, B->n_row, 0, B->n_row); 94 14.8k   H = isl_mat_lin_to_aff(H); 95 14.8k   C = isl_mat_inverse_product(H, C); 96 14.8k   if (!C) 97 0     goto error; 98 31.9k   for (i = 0; 14.8ki < B->n_row31.9k; ++i17.1k) { 99 18.5k     if (!18.5kisl_int_is_divisible_by18.5k(C->row[1+i][0], C->row[0][0])) 100 1.40k       break; 101 17.1k     isl_int_divexact17.1k(C->row[1+i][0], C->row[1+i][0], C->row[0][0]); 102 17.1k   } 103 14.8k   if (i < B->n_row) 104 1.40k     cst = isl_mat_alloc(B->ctx, B->n_row, 0); 105 14.8k   else 106 13.4k     cst = isl_mat_sub_alloc(C, 1, B->n_row, 0, 1); 107 14.8k   T = isl_mat_sub_alloc(U, B->n_row, B->n_col - 1, 0, B->n_row); 108 14.8k   cst = isl_mat_product(T, cst); 109 14.8k   isl_mat_free(M); 110 14.8k   isl_mat_free(C); 111 14.8k   isl_mat_free(U); 112 14.8k   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 14.8k } 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 10.8k { 131 10.8k   struct isl_mat *U; 132 10.8k 133 10.8k   U = isl_mat_alloc(B->ctx, B->n_col - 1, B->n_col - 1); 134 10.8k   if (!U) 135 0     return NULL; 136 10.8k   isl_seq_cpy(U->row[0], B->row[0] + 1, B->n_col - 1); 137 10.8k   U = isl_mat_unimodular_complete(U, 1); 138 10.8k   U = isl_mat_right_inverse(U); 139 10.8k   if (!U) 140 0     return NULL; 141 10.8k   isl_mat_col_mul(U, 0, d->block.data[0], 0); 142 10.8k   U = isl_mat_lin_to_aff(U); 143 10.8k   return U; 144 10.8k } 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 2.46k { 162 2.46k   int i, j, k; 163 2.46k   isl_int D; 164 2.46k   struct isl_mat *A = NULL, *U = NULL; 165 2.46k   struct isl_mat *T; 166 2.46k   unsigned size; 167 2.46k 168 2.46k   isl_int_init(D); 169 2.46k 170 2.46k   isl_vec_lcm(d, &D); 171 2.46k 172 2.46k   size = B->n_col - 1; 173 2.46k   A = isl_mat_alloc(B->ctx, size, B->n_row * size); 174 2.46k   U = isl_mat_alloc(B->ctx, size, size); 175 2.46k   if (!U || 2.46k!A2.46k) 176 0     goto error; 177 7.69k   for (i = 0; 2.46ki < B->n_row7.69k; ++i5.23k) { 178 5.23k     isl_seq_cpy(U->row[0], B->row[i] + 1, size); 179 5.23k     U = isl_mat_unimodular_complete(U, 1); 180 5.23k     if (!U) 181 0       goto error; 182 5.23k     isl_int_divexact5.23k(D, D, d->block.data[i]); 183 29.2k     for (k = 0; k < U->n_col29.2k; ++k24.0k) 184 24.0k       isl_int_mul(A->row[k][i*size+0], D, U->row[0][k]); 185 5.23k     isl_int_mul(D, D, d->block.data[i]); 186 24.0k     for (j = 1; j < U->n_row24.0k; ++j18.7k) 187 140k       for (k = 0; 18.7kk < U->n_col140k; ++k121k) 188 121k         isl_int_mul(A->row[k][i*size+j], 189 5.23k             D, U->row[j][k]); 190 5.23k   } 191 2.46k   A = isl_mat_left_hermite(A, 0, NULL, NULL); 192 2.46k   T = isl_mat_sub_alloc(A, 0, A->n_row, 0, A->n_row); 193 2.46k   T = isl_mat_lin_to_aff(T); 194 2.46k   if (!T) 195 0     goto error; 196 2.46k   isl_int_set2.46k(T->row[0][0], D); 197 2.46k   T = isl_mat_right_inverse(T); 198 2.46k   if (!T) 199 0     goto error; 200 2.46k   isl_assert2.46k(T->ctx, isl_int_is_one(T->row[0][0]), goto error); 201 2.46k   T = isl_mat_transpose(T); 202 2.46k   isl_mat_free(A); 203 2.46k   isl_mat_free(U); 204 2.46k 205 2.46k   isl_int_clear(D); 206 2.46k   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 2.46k } 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 14.8k { 311 14.8k   int i; 312 14.8k   struct isl_mat *cst = NULL; 313 14.8k   struct isl_mat *T = NULL; 314 14.8k   isl_int D; 315 14.8k 316 14.8k   if (!B || 14.8k!d14.8k) 317 0     goto error; 318 14.8k   isl_assert14.8k(B->ctx, B->n_row == d->size, goto error); 319 14.8k   cst = particular_solution(B, d); 320 14.8k   if (!cst) 321 0     goto error; 322 14.8k   if (14.8kcst->n_col == 014.8k) { 323 1.40k     T = isl_mat_alloc(B->ctx, B->n_col, 0); 324 1.40k     isl_mat_free(cst); 325 1.40k     isl_mat_free(B); 326 1.40k     isl_vec_free(d); 327 1.40k     return T; 328 1.40k   } 329 13.4k   isl_int_init13.4k(D); 330 13.4k   /* Replace a*g*row = 0 mod g*m by row = 0 mod m */ 331 29.7k   for (i = 0; i < B->n_row29.7k; ++i16.2k) { 332 16.2k     isl_seq_gcd(B->row[i] + 1, B->n_col - 1, &D); 333 16.2k     if (isl_int_is_one(D)) 334 15.4k       continue; 335 848     if (848isl_int_is_zero848(D)) { 336 228       B = isl_mat_drop_rows(B, i, 1); 337 228       d = isl_vec_cow(d); 338 228       if (!B || 228!d228) 339 0         goto error2; 340 228       isl_seq_cpy(d->block.data+i, d->block.data+i+1, 341 228               d->size - (i+1)); 342 228       d->size--; 343 228       i--; 344 228       continue; 345 228     } 346 620     B = isl_mat_cow(B); 347 620     if (!B) 348 0       goto error2; 349 620     isl_seq_scale_down(B->row[i] + 1, B->row[i] + 1, D, B->n_col-1); 350 620     isl_int_gcd(D, D, d->block.data[i]); 351 620     d = isl_vec_cow(d); 352 620     if (!d) 353 0       goto error2; 354 620     isl_int_divexact620(d->block.data[i], d->block.data[i], D); 355 620   } 356 13.4k   isl_int_clear13.4k(D); 357 13.4k   if (B->n_row == 0) 358 163     T = isl_mat_identity(B->ctx, B->n_col); 359 13.2k   else if (13.2kB->n_row == 113.2k) 360 10.8k     T = parameter_compression_1(B, d); 361 13.2k   else 362 2.46k     T = parameter_compression_multi(B, d); 363 13.4k   T = isl_mat_left_hermite(T, 0, NULL, NULL); 364 13.4k   if (!T) 365 0     goto error; 366 13.4k   isl_mat_sub_copy(T->ctx, T->row + 1, cst->row, cst->n_row, 0, 0, 1); 367 13.4k   isl_mat_free(cst); 368 13.4k   isl_mat_free(B); 369 13.4k   isl_vec_free(d); 370 13.4k   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 14.8k } 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 23 { 446 23   isl_mat_free(free1); 447 23   isl_mat_free(free2); 448 23   isl_mat_free(free3); 449 23   if (T223) { 450 13     isl_mat_free(*T2); 451 13     *T2 = isl_mat_alloc(ctx, 0, 1 + dim); 452 13   } 453 23   return isl_mat_alloc(ctx, 1 + dim, 0); 454 23 } 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 140k { 463 140k   int i; 464 140k 465 140k   if (nparam == 0) 466 140k     return mat; 467 1   if (1!mat1) 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 < nparam2; ++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 140k 479 140k   return mat; 480 140k } 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 140k { 532 140k   int i, n; 533 140k   isl_ctx *ctx; 534 140k   isl_mat *H = NULL, *C, *H1, *U = NULL, *U1, *U2; 535 140k   unsigned dim; 536 140k 537 140k   if (T2) 538 38.8k     *T2 = NULL; 539 140k   if (!B) 540 0     goto error; 541 140k 542 140k   ctx = isl_mat_get_ctx(B); 543 140k   dim = B->n_col - 1; 544 140k   n = dim - first; 545 140k   if (n < B->n_row) 546 0     isl_die(ctx, isl_error_invalid, "too many equality constraints", 547 140k       goto error); 548 140k   H = isl_mat_sub_alloc(B, 0, B->n_row, 1 + first, n); 549 140k   H = isl_mat_left_hermite(H, 0, &U, T2); 550 140k   if (!H || 140k!U140k || (T2 && 140k!*T238.8k)) 551 0     goto error; 552 140k   if (140kT2140k) { 553 38.8k     *T2 = isl_mat_drop_rows(*T2, 0, B->n_row); 554 38.8k     *T2 = isl_mat_diagonal(isl_mat_identity(ctx, 1 + first), *T2); 555 38.8k     if (!*T2) 556 0       goto error; 557 140k   } 558 140k   C = isl_mat_alloc(ctx, 1 + B->n_row, 1 + first); 559 140k   if (!C) 560 0     goto error; 561 140k   isl_int_set_si140k(C->row[0][0], 1); 562 140k   isl_seq_clr(C->row[0] + 1, first); 563 140k   isl_mat_sub_neg(ctx, C->row + 1, B->row, B->n_row, 0, 0, 1 + first); 564 140k   H1 = isl_mat_sub_alloc(H, 0, H->n_row, 0, H->n_row); 565 140k   H1 = isl_mat_lin_to_aff(H1); 566 140k   C = isl_mat_inverse_product(H1, C); 567 140k   if (!C) 568 0     goto error; 569 140k   isl_mat_free(H); 570 140k   if (!140kisl_int_is_one140k(C->row[0][0])) { 571 794     isl_int g; 572 794 573 794     isl_int_init(g); 574 4.00k     for (i = 0; i < B->n_row4.00k; ++i3.21k) { 575 3.23k       isl_seq_gcd(C->row[1 + i] + 1, first, &g); 576 3.23k       isl_int_gcd(g, g, C->row[0][0]); 577 3.23k       if (!3.23kisl_int_is_divisible_by3.23k(C->row[1 + i][0], g)) 578 23         break; 579 3.23k     } 580 794     isl_int_clear(g); 581 794 582 794     if (i < B->n_row) 583 23       return empty_compression(ctx, dim, T2, B, C, U); 584 771     C = isl_mat_normalize(C); 585 771   } 586 140k   U1 = isl_mat_sub_alloc(U, 0, U->n_row, 0, B->n_row); 587 140k   U1 = isl_mat_lin_to_aff(U1); 588 140k   U2 = isl_mat_sub_alloc(U, 0, U->n_row, B->n_row, U->n_row - B->n_row); 589 140k   U2 = isl_mat_lin_to_aff(U2); 590 140k   isl_mat_free(U); 591 140k   C = isl_mat_product(U1, C); 592 140k   C = isl_mat_aff_direct_sum(C, U2); 593 140k   C = insert_parameter_rows(C, first); 594 140k 595 140k   isl_mat_free(B); 596 140k 597 140k   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 (T20) { 603 0     isl_mat_free(*T2); 604 0     *T2 = NULL; 605 0   } 606 0   return NULL; 607 140k } 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 140k { 627 140k   return isl_mat_final_variable_compression(B, 0, T2); 628 140k } 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 10 { 637 10   unsigned dim; 638 10   isl_mat *id; 639 10 640 10   if (!bset) 641 0     return NULL; 642 10   if (10!T && 10!T20) 643 0     return bset; 644 10 645 10   dim = isl_basic_set_dim(bset, isl_dim_set); 646 10   id = isl_mat_identity(isl_basic_map_get_ctx(bset), 1 + dim); 647 10   if (T) 648 10     *T = isl_mat_copy(id); 649 10   if (T2) 650 0     *T2 = isl_mat_copy(id); 651 10   isl_mat_free(id); 652 10 653 10   return bset; 654 10 } 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 120k { 666 120k   struct isl_mat *B, *TC; 667 120k   unsigned dim; 668 120k 669 120k   if (T) 670 120k     *T = NULL; 671 120k   if (T2) 672 29.5k     *T2 = NULL; 673 120k   if (!bset) 674 0     goto error; 675 120k   isl_assert120k(bset->ctx, isl_basic_set_n_param(bset) == 0, goto error); 676 120k   isl_assert120k(bset->ctx, bset->n_div == 0, goto error); 677 120k   dim = isl_basic_set_n_dim(bset); 678 120k   isl_assert(bset->ctx, bset->n_eq <= dim, goto error); 679 120k   if (120kbset->n_eq == 0120k) 680 0     return return_with_identity(bset, T, T2); 681 120k 682 120k   B = isl_mat_sub_alloc6(bset->ctx, bset->eq, 0, bset->n_eq, 0, 1 + dim); 683 120k   TC = isl_mat_variable_compression(B, T2); 684 120k   if (!TC) 685 0     goto error; 686 120k   if (120kTC->n_col == 0120k) { 687 10     isl_mat_free(TC); 688 10     if (T210) { 689 0       isl_mat_free(*T2); 690 0       *T2 = NULL; 691 0     } 692 10     bset = isl_basic_set_set_to_empty(bset); 693 10     return return_with_identity(bset, T, T2); 694 10   } 695 120k 696 120k   bset = isl_basic_set_preimage(bset, T ? 120kisl_mat_copy(TC)120k : TC0); 697 120k   if (T) 698 120k     *T = TC; 699 120k   return bset; 700 0 error: 701 0   isl_basic_set_free(bset); 702 0   return NULL; 703 120k } 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 120k { 708 120k   if (T) 709 120k     *T = NULL; 710 120k   if (T2) 711 29.5k     *T2 = NULL; 712 120k   if (!bset) 713 0     return NULL; 714 120k   isl_assert120k(bset->ctx, isl_basic_set_n_param(bset) == 0, goto error); 715 120k   bset = isl_basic_set_gauss(bset, NULL); 716 120k   if (ISL_F_ISSET(bset, ISL_BASIC_SET_EMPTY)) 717 0     return return_with_identity(bset, T, T2); 718 120k   bset = compress_variables(bset, T, T2); 719 120k   return bset; 720 0 error: 721 0   isl_basic_set_free(bset); 722 0   *T = NULL; 723 0   return NULL; 724 120k } 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 74 { 738 74   isl_bool fixed; 739 74   struct isl_ctx *ctx; 740 74   struct isl_mat *H = NULL, *U = NULL, *C, *H1, *U1; 741 74   unsigned total; 742 74   unsigned nparam; 743 74 744 74   if (!bset || 74!modulo74 || !residue74) 745 0     return isl_stat_error; 746 74 747 74   fixed = isl_basic_set_plain_dim_is_fixed(bset, pos, residue); 748 74   if (fixed < 0) 749 0     return isl_stat_error; 750 74   if (74fixed74) { 751 0     isl_int_set_si(*modulo, 0); 752 0     return isl_stat_ok; 753 0   } 754 74 755 74   ctx = isl_basic_set_get_ctx(bset); 756 74   total = isl_basic_set_total_dim(bset); 757 74   nparam = isl_basic_set_n_param(bset); 758 74   H = isl_mat_sub_alloc6(ctx, bset->eq, 0, bset->n_eq, 1, total); 759 74   H = isl_mat_left_hermite(H, 0, &U, NULL); 760 74   if (!H) 761 0     return isl_stat_error; 762 74 763 74   isl_seq_gcd(U->row[nparam + pos]+bset->n_eq, 764 74       total-bset->n_eq, modulo); 765 74   if (isl_int_is_zero(*modulo)) 766 0     isl_int_set_si(*modulo, 1); 767 74   if (isl_int_is_one74(*modulo)) { 768 59     isl_int_set_si(*residue, 0); 769 59     isl_mat_free(H); 770 59     isl_mat_free(U); 771 59     return isl_stat_ok; 772 59   } 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_si15(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 (15!15isl_int_is_divisible_by15(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_divexact15(*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 74 } 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 73 { 819 73   isl_int m; 820 73   isl_int r; 821 73   int i; 822 73 823 73   if (!set || 73!modulo73 || !residue73) 824 0     return isl_stat_error; 825 73 826 73   if (73set->n == 073) { 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 73 832 73   if (73isl_basic_set_dim_residue_class(set->p[0], pos, modulo, residue)<073) 833 0     return isl_stat_error; 834 73 835 73   if (73set->n == 173) 836 64     return isl_stat_ok; 837 9 838 9   if (9isl_int_is_one9(*modulo)) 839 8     return isl_stat_ok; 840 1 841 1   isl_int_init1(m); 842 1   isl_int_init(r); 843 1 844 2   for (i = 1; i < set->n2; ++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_gcd1(*modulo, *modulo, m); 848 1     isl_int_sub(m, *residue, r); 849 1     isl_int_gcd(*modulo, *modulo, m); 850 1     if (!1isl_int_is_zero1(*modulo)) 851 1       isl_int_fdiv_r(*residue, *residue, *modulo); 852 1     if (isl_int_is_one(*modulo)) 853 0       break; 854 1   } 855 1 856 1   isl_int_clear1(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 73 } 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 72 { 878 72   *modulo = NULL; 879 72   *residue = NULL; 880 72   if (!set) 881 0     return isl_stat_error; 882 72   *modulo = isl_val_alloc(isl_set_get_ctx(set)); 883 72   *residue = isl_val_alloc(isl_set_get_ctx(set)); 884 72   if (!*modulo || 72!*residue72) 885 0     goto error; 886 72   if (72isl_set_dim_residue_class(set, pos, 887 72           &(*modulo)->n, &(*residue)->n) < 0) 888 0     goto error; 889 72   isl_int_set_si72((*modulo)->d, 1); 890 72   isl_int_set_si((*residue)->d, 1); 891 72   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 72 }