3 #include "isl_map_private.h"
4 #include "isl_equalities.h"
6 /* Given a set of modulo constraints
10 * this function computes a particular solution y_0
12 * The input is given as a matrix B = [ c A ] and a vector d.
14 * The output is matrix containing the solution y_0 or
15 * a zero-column matrix if the constraints admit no integer solution.
17 * The given set of constrains is equivalent to
21 * with D = diag d and x a fresh set of variables.
22 * Reducing both c and A modulo d does not change the
23 * value of y in the solution and may lead to smaller coefficients.
24 * Let M = [ D A ] and [ H 0 ] = M U, the Hermite normal form of M.
30 * [ H 0 ] U^{-1} [ y ] = - c
33 * [ B ] = U^{-1} [ y ]
37 * so B may be chosen arbitrarily, e.g., B = 0, and then
40 * U^{-1} [ y ] = [ 0 ]
48 * If any of the coordinates of this y are non-integer
49 * then the constraints admit no integer solution and
50 * a zero-column matrix is returned.
52 static struct isl_mat *particular_solution(struct isl_ctx *ctx,
53 struct isl_mat *B, struct isl_vec *d)
56 struct isl_mat *M = NULL;
57 struct isl_mat *C = NULL;
58 struct isl_mat *U = NULL;
59 struct isl_mat *H = NULL;
60 struct isl_mat *cst = NULL;
61 struct isl_mat *T = NULL;
63 M = isl_mat_alloc(ctx, B->n_row, B->n_row + B->n_col - 1);
64 C = isl_mat_alloc(ctx, 1 + B->n_row, 1);
67 isl_int_set_si(C->row[0][0], 1);
68 for (i = 0; i < B->n_row; ++i) {
69 isl_seq_clr(M->row[i], B->n_row);
70 isl_int_set(M->row[i][i], d->block.data[i]);
71 isl_int_neg(C->row[1 + i][0], B->row[i][0]);
72 isl_int_fdiv_r(C->row[1+i][0], C->row[1+i][0], M->row[i][i]);
73 for (j = 0; j < B->n_col - 1; ++j)
74 isl_int_fdiv_r(M->row[i][B->n_row + j],
75 B->row[i][1 + j], M->row[i][i]);
77 M = isl_mat_left_hermite(ctx, M, 0, &U, NULL);
80 H = isl_mat_sub_alloc(ctx, M->row, 0, B->n_row, 0, B->n_row);
81 H = isl_mat_lin_to_aff(ctx, H);
82 C = isl_mat_inverse_product(ctx, H, C);
85 for (i = 0; i < B->n_row; ++i) {
86 if (!isl_int_is_divisible_by(C->row[1+i][0], C->row[0][0]))
88 isl_int_divexact(C->row[1+i][0], C->row[1+i][0], C->row[0][0]);
91 cst = isl_mat_alloc(ctx, B->n_row, 0);
93 cst = isl_mat_sub_alloc(ctx, C->row, 1, B->n_row, 0, 1);
94 T = isl_mat_sub_alloc(ctx, U->row, B->n_row, B->n_col - 1, 0, B->n_row);
95 cst = isl_mat_product(ctx, T, cst);
101 isl_mat_free(ctx, M);
102 isl_mat_free(ctx, C);
103 isl_mat_free(ctx, U);
107 static struct isl_mat *unimodular_complete(struct isl_ctx *ctx,
108 struct isl_mat *M, int row)
111 struct isl_mat *H = NULL, *Q = NULL;
113 isl_assert(ctx, M->n_row == M->n_col, goto error);
115 H = isl_mat_left_hermite(ctx, isl_mat_copy(ctx, M), 0, NULL, &Q);
119 for (r = 0; r < row; ++r)
120 isl_assert(ctx, isl_int_is_one(H->row[r][r]), goto error);
121 for (r = row; r < M->n_row; ++r)
122 isl_seq_cpy(M->row[r], Q->row[r], M->n_col);
123 isl_mat_free(ctx, H);
124 isl_mat_free(ctx, Q);
127 isl_mat_free(ctx, H);
128 isl_mat_free(ctx, Q);
129 isl_mat_free(ctx, M);
133 /* Compute and return the matrix
135 * U_1^{-1} diag(d_1, 1, ..., 1)
137 * with U_1 the unimodular completion of the first (and only) row of B.
138 * The columns of this matrix generate the lattice that satisfies
139 * the single (linear) modulo constraint.
141 static struct isl_mat *parameter_compression_1(struct isl_ctx *ctx,
142 struct isl_mat *B, struct isl_vec *d)
146 U = isl_mat_alloc(ctx, B->n_col - 1, B->n_col - 1);
149 isl_seq_cpy(U->row[0], B->row[0] + 1, B->n_col - 1);
150 U = unimodular_complete(ctx, U, 1);
151 U = isl_mat_right_inverse(ctx, U);
154 isl_mat_col_mul(U, 0, d->block.data[0], 0);
155 U = isl_mat_lin_to_aff(ctx, U);
158 isl_mat_free(ctx, U);
162 /* Compute a common lattice of solutions to the linear modulo
163 * constraints specified by B and d.
164 * See also the documentation of isl_mat_parameter_compression.
167 * A = [ L_1^{-T} L_2^{-T} ... L_k^{-T} ]
169 * on a common denominator. This denominator D is the lcm of modulos d.
170 * Since L_i = U_i^{-1} diag(d_i, 1, ... 1), we have
171 * L_i^{-T} = U_i^T diag(d_i, 1, ... 1)^{-T} = U_i^T diag(1/d_i, 1, ..., 1).
172 * Putting this on the common denominator, we have
173 * D * L_i^{-T} = U_i^T diag(D/d_i, D, ..., D).
175 static struct isl_mat *parameter_compression_multi(struct isl_ctx *ctx,
176 struct isl_mat *B, struct isl_vec *d)
181 struct isl_mat *A = NULL, *U = NULL;
187 isl_vec_lcm(ctx, d, &D);
190 A = isl_mat_alloc(ctx, size, B->n_row * size);
191 U = isl_mat_alloc(ctx, size, size);
194 for (i = 0; i < B->n_row; ++i) {
195 isl_seq_cpy(U->row[0], B->row[i] + 1, size);
196 U = unimodular_complete(ctx, U, 1);
199 isl_int_divexact(D, D, d->block.data[i]);
200 for (k = 0; k < U->n_col; ++k)
201 isl_int_mul(A->row[k][i*size+0], D, U->row[0][k]);
202 isl_int_mul(D, D, d->block.data[i]);
203 for (j = 1; j < U->n_row; ++j)
204 for (k = 0; k < U->n_col; ++k)
205 isl_int_mul(A->row[k][i*size+j],
208 A = isl_mat_left_hermite(ctx, A, 0, NULL, NULL);
209 T = isl_mat_sub_alloc(ctx, A->row, 0, A->n_row, 0, A->n_row);
210 T = isl_mat_lin_to_aff(ctx, T);
211 isl_int_set(T->row[0][0], D);
212 T = isl_mat_right_inverse(ctx, T);
213 isl_assert(ctx, isl_int_is_one(T->row[0][0]), goto error);
214 T = isl_mat_transpose(ctx, T);
215 isl_mat_free(ctx, A);
216 isl_mat_free(ctx, U);
221 isl_mat_free(ctx, A);
222 isl_mat_free(ctx, U);
227 /* Given a set of modulo constraints
231 * this function returns an affine transformation T,
235 * that bijectively maps the integer vectors y' to integer
236 * vectors y that satisfy the modulo constraints.
238 * This function is inspired by Section 2.5.3
239 * of B. Meister, "Stating and Manipulating Periodicity in the Polytope
240 * Model. Applications to Program Analysis and Optimization".
241 * However, the implementation only follows the algorithm of that
242 * section for computing a particular solution and not for computing
243 * a general homogeneous solution. The latter is incomplete and
244 * may remove some valid solutions.
245 * Instead, we use an adaptation of the algorithm in Section 7 of
246 * B. Meister, S. Verdoolaege, "Polynomial Approximations in the Polytope
247 * Model: Bringing the Power of Quasi-Polynomials to the Masses".
249 * The input is given as a matrix B = [ c A ] and a vector d.
250 * Each element of the vector d corresponds to a row in B.
251 * The output is a lower triangular matrix.
252 * If no integer vector y satisfies the given constraints then
253 * a matrix with zero columns is returned.
255 * We first compute a particular solution y_0 to the given set of
256 * modulo constraints in particular_solution. If no such solution
257 * exists, then we return a zero-columned transformation matrix.
258 * Otherwise, we compute the generic solution to
262 * That is we want to compute G such that
266 * with y'' integer, describes the set of solutions.
268 * We first remove the common factors of each row.
269 * In particular if gcd(A_i,d_i) != 1, then we divide the whole
270 * row i (including d_i) by this common factor. If afterwards gcd(A_i) != 1,
271 * then we divide this row of A by the common factor, unless gcd(A_i) = 0.
272 * In the later case, we simply drop the row (in both A and d).
274 * If there are no rows left in A, the G is the identity matrix. Otherwise,
275 * for each row i, we now determine the lattice of integer vectors
276 * that satisfies this row. Let U_i be the unimodular extension of the
277 * row A_i. This unimodular extension exists because gcd(A_i) = 1.
278 * The first component of
282 * needs to be a multiple of d_i. Let y' = diag(d_i, 1, ..., 1) y''.
285 * y = U_i^{-1} diag(d_i, 1, ..., 1) y''
287 * for arbitrary integer vectors y''. That is, y belongs to the lattice
288 * generated by the columns of L_i = U_i^{-1} diag(d_i, 1, ..., 1).
289 * If there is only one row, then G = L_1.
291 * If there is more than one row left, we need to compute the intersection
292 * of the lattices. That is, we need to compute an L such that
294 * L = L_i L_i' for all i
296 * with L_i' some integer matrices. Let A be constructed as follows
298 * A = [ L_1^{-T} L_2^{-T} ... L_k^{-T} ]
300 * and computed the Hermite Normal Form of A = [ H 0 ] U
303 * L_i^{-T} = H U_{1,i}
307 * H^{-T} = L_i U_{1,i}^T
309 * In other words G = L = H^{-T}.
310 * To ensure that G is lower triangular, we compute and use its Hermite
313 * The affine transformation matrix returned is then
318 * as any y = y_0 + G y' with y' integer is a solution to the original
319 * modulo constraints.
321 struct isl_mat *isl_mat_parameter_compression(struct isl_ctx *ctx,
322 struct isl_mat *B, struct isl_vec *d)
325 struct isl_mat *cst = NULL;
326 struct isl_mat *T = NULL;
331 isl_assert(ctx, B->n_row == d->size, goto error);
332 cst = particular_solution(ctx, B, d);
335 if (cst->n_col == 0) {
336 T = isl_mat_alloc(ctx, B->n_col, 0);
337 isl_mat_free(ctx, cst);
338 isl_mat_free(ctx, B);
339 isl_vec_free(ctx, d);
343 /* Replace a*g*row = 0 mod g*m by row = 0 mod m */
344 for (i = 0; i < B->n_row; ++i) {
345 isl_seq_gcd(B->row[i] + 1, B->n_col - 1, &D);
346 if (isl_int_is_one(D))
348 if (isl_int_is_zero(D)) {
349 B = isl_mat_drop_rows(ctx, B, i, 1);
350 d = isl_vec_cow(ctx, d);
353 isl_seq_cpy(d->block.data+i, d->block.data+i+1,
359 B = isl_mat_cow(ctx, B);
362 isl_seq_scale_down(B->row[i] + 1, B->row[i] + 1, D, B->n_col-1);
363 isl_int_gcd(D, D, d->block.data[i]);
364 d = isl_vec_cow(ctx, d);
367 isl_int_divexact(d->block.data[i], d->block.data[i], D);
371 T = isl_mat_identity(ctx, B->n_col);
372 else if (B->n_row == 1)
373 T = parameter_compression_1(ctx, B, d);
375 T = parameter_compression_multi(ctx, B, d);
376 T = isl_mat_left_hermite(ctx, T, 0, NULL, NULL);
379 isl_mat_sub_copy(ctx, T->row + 1, cst->row, cst->n_row, 0, 0, 1);
380 isl_mat_free(ctx, cst);
381 isl_mat_free(ctx, B);
382 isl_vec_free(ctx, d);
387 isl_mat_free(ctx, cst);
388 isl_mat_free(ctx, B);
389 isl_vec_free(ctx, d);
393 /* Given a set of equalities
397 * this function computes unimodular transformation from a lower-dimensional
398 * space to the original space that bijectively maps the integer points x'
399 * in the lower-dimensional space to the integer points x in the original
400 * space that satisfy the equalities.
402 * The input is given as a matrix B = [ -c M ] and the out is a
403 * matrix that maps [1 x'] to [1 x].
404 * If T2 is not NULL, then *T2 is set to a matrix mapping [1 x] to [1 x'].
406 * First compute the (left) Hermite normal form of M,
408 * M [U1 U2] = M U = H = [H1 0]
410 * M = H Q = [H1 0] [Q1]
413 * with U, Q unimodular, Q = U^{-1} (and H lower triangular).
414 * Define the transformed variables as
416 * x = [U1 U2] [ x1' ] = [U1 U2] [Q1] x
419 * The equalities then become
421 * H1 x1' - c = 0 or x1' = H1^{-1} c = c'
423 * If any of the c' is non-integer, then the original set has no
424 * integer solutions (since the x' are a unimodular transformation
426 * Otherwise, the transformation is given by
428 * x = U1 H1^{-1} c + U2 x2'
430 * The inverse transformation is simply
434 struct isl_mat *isl_mat_variable_compression(struct isl_ctx *ctx,
435 struct isl_mat *B, struct isl_mat **T2)
438 struct isl_mat *H = NULL, *C = NULL, *H1, *U = NULL, *U1, *U2, *TC;
447 H = isl_mat_sub_alloc(ctx, B->row, 0, B->n_row, 1, dim);
448 H = isl_mat_left_hermite(ctx, H, 0, &U, T2);
449 if (!H || !U || (T2 && !*T2))
452 *T2 = isl_mat_drop_rows(ctx, *T2, 0, B->n_row);
453 *T2 = isl_mat_lin_to_aff(ctx, *T2);
457 C = isl_mat_alloc(ctx, 1+B->n_row, 1);
460 isl_int_set_si(C->row[0][0], 1);
461 isl_mat_sub_neg(ctx, C->row+1, B->row, B->n_row, 0, 0, 1);
462 H1 = isl_mat_sub_alloc(ctx, H->row, 0, H->n_row, 0, H->n_row);
463 H1 = isl_mat_lin_to_aff(ctx, H1);
464 TC = isl_mat_inverse_product(ctx, H1, C);
467 isl_mat_free(ctx, H);
468 if (!isl_int_is_one(TC->row[0][0])) {
469 for (i = 0; i < B->n_row; ++i) {
470 if (!isl_int_is_divisible_by(TC->row[1+i][0], TC->row[0][0])) {
471 isl_mat_free(ctx, B);
472 isl_mat_free(ctx, TC);
473 isl_mat_free(ctx, U);
475 isl_mat_free(ctx, *T2);
478 return isl_mat_alloc(ctx, 1 + dim, 0);
480 isl_seq_scale_down(TC->row[1+i], TC->row[1+i], TC->row[0][0], 1);
482 isl_int_set_si(TC->row[0][0], 1);
484 U1 = isl_mat_sub_alloc(ctx, U->row, 0, U->n_row, 0, B->n_row);
485 U1 = isl_mat_lin_to_aff(ctx, U1);
486 U2 = isl_mat_sub_alloc(ctx, U->row, 0, U->n_row,
487 B->n_row, U->n_row - B->n_row);
488 U2 = isl_mat_lin_to_aff(ctx, U2);
489 isl_mat_free(ctx, U);
490 TC = isl_mat_product(ctx, U1, TC);
491 TC = isl_mat_aff_direct_sum(ctx, TC, U2);
493 isl_mat_free(ctx, B);
497 isl_mat_free(ctx, B);
498 isl_mat_free(ctx, H);
499 isl_mat_free(ctx, U);
501 isl_mat_free(ctx, *T2);
507 /* Use the n equalities of bset to unimodularly transform the
508 * variables x such that n transformed variables x1' have a constant value
509 * and rewrite the constraints of bset in terms of the remaining
510 * transformed variables x2'. The matrix pointed to by T maps
511 * the new variables x2' back to the original variables x, while T2
512 * maps the original variables to the new variables.
514 static struct isl_basic_set *compress_variables(struct isl_ctx *ctx,
515 struct isl_basic_set *bset, struct isl_mat **T, struct isl_mat **T2)
517 struct isl_mat *B, *TC;
526 isl_assert(ctx, isl_basic_set_n_param(bset) == 0, goto error);
527 isl_assert(ctx, bset->n_div == 0, goto error);
528 dim = isl_basic_set_n_dim(bset);
529 isl_assert(ctx, bset->n_eq <= dim, goto error);
533 B = isl_mat_sub_alloc(ctx, bset->eq, 0, bset->n_eq, 0, 1 + dim);
534 TC = isl_mat_variable_compression(ctx, B, T2);
537 if (TC->n_col == 0) {
538 isl_mat_free(ctx, TC);
540 isl_mat_free(ctx, *T2);
543 return isl_basic_set_set_to_empty(bset);
546 bset = isl_basic_set_preimage(bset, T ? isl_mat_copy(ctx, TC) : TC);
551 isl_basic_set_free(bset);
555 struct isl_basic_set *isl_basic_set_remove_equalities(
556 struct isl_basic_set *bset, struct isl_mat **T, struct isl_mat **T2)
564 isl_assert(bset->ctx, isl_basic_set_n_param(bset) == 0, goto error);
565 bset = isl_basic_set_gauss(bset, NULL);
566 if (ISL_F_ISSET(bset, ISL_BASIC_SET_EMPTY))
568 bset = compress_variables(bset->ctx, bset, T, T2);
571 isl_basic_set_free(bset);
576 /* Check if dimension dim belongs to a residue class
577 * i_dim \equiv r mod m
578 * with m != 1 and if so return m in *modulo and r in *residue.
580 int isl_basic_set_dim_residue_class(struct isl_basic_set *bset,
581 int pos, isl_int *modulo, isl_int *residue)
584 struct isl_mat *H = NULL, *U = NULL, *C, *H1, *U1;
588 if (!bset || !modulo || !residue)
592 total = isl_basic_set_total_dim(bset);
593 nparam = isl_basic_set_n_param(bset);
594 H = isl_mat_sub_alloc(ctx, bset->eq, 0, bset->n_eq, 1, total);
595 H = isl_mat_left_hermite(ctx, H, 0, &U, NULL);
599 isl_seq_gcd(U->row[nparam + pos]+bset->n_eq,
600 total-bset->n_eq, modulo);
601 if (isl_int_is_zero(*modulo) || isl_int_is_one(*modulo)) {
602 isl_int_set_si(*residue, 0);
603 isl_mat_free(ctx, H);
604 isl_mat_free(ctx, U);
608 C = isl_mat_alloc(ctx, 1+bset->n_eq, 1);
611 isl_int_set_si(C->row[0][0], 1);
612 isl_mat_sub_neg(ctx, C->row+1, bset->eq, bset->n_eq, 0, 0, 1);
613 H1 = isl_mat_sub_alloc(ctx, H->row, 0, H->n_row, 0, H->n_row);
614 H1 = isl_mat_lin_to_aff(ctx, H1);
615 C = isl_mat_inverse_product(ctx, H1, C);
616 isl_mat_free(ctx, H);
617 U1 = isl_mat_sub_alloc(ctx, U->row, nparam+pos, 1, 0, bset->n_eq);
618 U1 = isl_mat_lin_to_aff(ctx, U1);
619 isl_mat_free(ctx, U);
620 C = isl_mat_product(ctx, U1, C);
623 if (!isl_int_is_divisible_by(C->row[1][0], C->row[0][0])) {
624 bset = isl_basic_set_copy(bset);
625 bset = isl_basic_set_set_to_empty(bset);
626 isl_basic_set_free(bset);
627 isl_int_set_si(*modulo, 0);
628 isl_int_set_si(*residue, 0);
631 isl_int_divexact(*residue, C->row[1][0], C->row[0][0]);
632 isl_int_fdiv_r(*residue, *residue, *modulo);
633 isl_mat_free(ctx, C);
636 isl_mat_free(ctx, H);
637 isl_mat_free(ctx, U);