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