3 #include "isl_map_private.h"
4 #include "isl_equalities.h"
6 /* Use the n equalities of bset to unimodularly transform the
7 * variables x such that n transformed variables x1' have a constant value
8 * and rewrite the constraints of bset in terms of the remaining
9 * transformed variables x2'. The matrix pointed to by T maps
10 * the new variables x2' back to the original variables x, while T2
11 * maps the original variables to the new variables.
13 * Let the equalities of bset be
17 * Compute the (left) Hermite normal form of M,
19 * M [U1 U2] = M U = H = [H1 0]
21 * M = H Q = [H1 0] [Q1]
24 * with U, Q unimodular, Q = U^{-1} (and H lower triangular).
25 * Define the transformed variables as
27 * x = [U1 U2] [ x1' ] = [U1 U2] [Q1] x
30 * The equalities then become
32 * H1 x1' - c = 0 or x1' = H1^{-1} c = c'
34 * If any of the c' is non-integer, then the original set has no
35 * integer solutions (since the x' are a unimodular transformation
37 * Otherwise, the transformation is given by
39 * x = U1 H1^{-1} c + U2 x2'
41 * The inverse transformation is simply
45 static struct isl_basic_set *compress_variables(struct isl_ctx *ctx,
46 struct isl_basic_set *bset, struct isl_mat **T, struct isl_mat **T2)
49 struct isl_mat *H = NULL, *C = NULL, *H1, *U = NULL, *U1, *U2, *TC;
58 isl_assert(ctx, isl_basic_set_n_param(bset) == 0, goto error);
59 isl_assert(ctx, bset->n_div == 0, goto error);
60 dim = isl_basic_set_n_dim(bset);
61 isl_assert(ctx, bset->n_eq <= dim, goto error);
65 H = isl_mat_sub_alloc(ctx, bset->eq, 0, bset->n_eq, 1, dim);
66 H = isl_mat_left_hermite(ctx, H, 0, &U, T2);
67 if (!H || !U || (T2 && !*T2))
70 *T2 = isl_mat_drop_rows(ctx, *T2, 0, bset->n_eq);
71 *T2 = isl_mat_lin_to_aff(ctx, *T2);
75 C = isl_mat_alloc(ctx, 1+bset->n_eq, 1);
78 isl_int_set_si(C->row[0][0], 1);
79 isl_mat_sub_neg(ctx, C->row+1, bset->eq, bset->n_eq, 0, 0, 1);
80 H1 = isl_mat_sub_alloc(ctx, H->row, 0, H->n_row, 0, H->n_row);
81 H1 = isl_mat_lin_to_aff(ctx, H1);
82 TC = isl_mat_inverse_product(ctx, H1, C);
86 if (!isl_int_is_one(TC->row[0][0])) {
87 for (i = 0; i < bset->n_eq; ++i) {
88 if (!isl_int_is_divisible_by(TC->row[1+i][0], TC->row[0][0])) {
89 isl_mat_free(ctx, TC);
92 isl_mat_free(ctx, *T2);
95 return isl_basic_set_set_to_empty(bset);
97 isl_seq_scale_down(TC->row[1+i], TC->row[1+i], TC->row[0][0], 1);
99 isl_int_set_si(TC->row[0][0], 1);
101 U1 = isl_mat_sub_alloc(ctx, U->row, 0, U->n_row, 0, bset->n_eq);
102 U1 = isl_mat_lin_to_aff(ctx, U1);
103 U2 = isl_mat_sub_alloc(ctx, U->row, 0, U->n_row,
104 bset->n_eq, U->n_row - bset->n_eq);
105 U2 = isl_mat_lin_to_aff(ctx, U2);
106 isl_mat_free(ctx, U);
107 TC = isl_mat_product(ctx, U1, TC);
108 TC = isl_mat_aff_direct_sum(ctx, TC, U2);
109 bset = isl_basic_set_preimage(ctx, bset, T ? isl_mat_copy(ctx, TC) : TC);
114 isl_mat_free(ctx, H);
115 isl_mat_free(ctx, U);
117 isl_mat_free(ctx, *T2);
118 isl_basic_set_free(bset);
126 struct isl_basic_set *isl_basic_set_remove_equalities(
127 struct isl_basic_set *bset, struct isl_mat **T, struct isl_mat **T2)
135 isl_assert(bset->ctx, isl_basic_set_n_param(bset) == 0, goto error);
136 bset = isl_basic_set_gauss(bset, NULL);
137 if (F_ISSET(bset, ISL_BASIC_SET_EMPTY))
139 bset = compress_variables(bset->ctx, bset, T, T2);
142 isl_basic_set_free(bset);
147 /* Check if dimension dim belongs to a residue class
148 * i_dim \equiv r mod m
149 * with m != 1 and if so return m in *modulo and r in *residue.
151 int isl_basic_set_dim_residue_class(struct isl_basic_set *bset,
152 int pos, isl_int *modulo, isl_int *residue)
155 struct isl_mat *H = NULL, *U = NULL, *C, *H1, *U1;
159 if (!bset || !modulo || !residue)
163 total = isl_basic_set_total_dim(bset);
164 nparam = isl_basic_set_n_param(bset);
165 H = isl_mat_sub_alloc(ctx, bset->eq, 0, bset->n_eq, 1, total);
166 H = isl_mat_left_hermite(ctx, H, 0, &U, NULL);
170 isl_seq_gcd(U->row[nparam + pos]+bset->n_eq,
171 total-bset->n_eq, modulo);
172 if (isl_int_is_zero(*modulo) || isl_int_is_one(*modulo)) {
173 isl_int_set_si(*residue, 0);
174 isl_mat_free(ctx, H);
175 isl_mat_free(ctx, U);
179 C = isl_mat_alloc(ctx, 1+bset->n_eq, 1);
182 isl_int_set_si(C->row[0][0], 1);
183 isl_mat_sub_neg(ctx, C->row+1, bset->eq, bset->n_eq, 0, 0, 1);
184 H1 = isl_mat_sub_alloc(ctx, H->row, 0, H->n_row, 0, H->n_row);
185 H1 = isl_mat_lin_to_aff(ctx, H1);
186 C = isl_mat_inverse_product(ctx, H1, C);
187 isl_mat_free(ctx, H);
188 U1 = isl_mat_sub_alloc(ctx, U->row, nparam+pos, 1, 0, bset->n_eq);
189 U1 = isl_mat_lin_to_aff(ctx, U1);
190 isl_mat_free(ctx, U);
191 C = isl_mat_product(ctx, U1, C);
194 if (!isl_int_is_divisible_by(C->row[1][0], C->row[0][0])) {
195 bset = isl_basic_set_copy(bset);
196 bset = isl_basic_set_set_to_empty(bset);
197 isl_basic_set_free(bset);
198 isl_int_set_si(*modulo, 0);
199 isl_int_set_si(*residue, 0);
202 isl_int_divexact(*residue, C->row[1][0], C->row[0][0]);
203 isl_int_fdiv_r(*residue, *residue, *modulo);
204 isl_mat_free(ctx, C);
207 isl_mat_free(ctx, H);
208 isl_mat_free(ctx, U);