isl_equalities.c

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