isl_equalities.c

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