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