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

source code of polly/lib/External/isl/isl_equalities.c