isl_equalities.c

   1 #include "isl_mat.h"
   2 #include "isl_seq.h"
   3 #include "isl_map_private.h"
   4 #include "isl_equalities.h"
   5
   6 /* Given a set of equalities
   7  *
   8  *              M x - c = 0
   9  *
  10  * this function computes unimodular transformation from a lower-dimensional
  11  * space to the original space that bijectively maps the integer points x'
  12  * in the lower-dimensional space to the integer points x in the original
  13  * space that satisfy the equalities.
  14  *
  15  * The input is given as a matrix B = [ -c M ] and the out is a
  16  * matrix that maps [1 x'] to [1 x].
  17  * If T2 is not NULL, then *T2 is set to a matrix mapping [1 x] to [1 x'].
  18  *
  19  * First compute the (left) Hermite normal form of M,
  20  *
  21  *              M [U1 U2] = M U = H = [H1 0]
  22  * or
  23  *                            M = H Q = [H1 0] [Q1]
  24  *                                             [Q2]
  25  *
  26  * with U, Q unimodular, Q = U^{-1} (and H lower triangular).
  27  * Define the transformed variables as
  28  *
  29  *              x = [U1 U2] [ x1' ] = [U1 U2] [Q1] x
  30  *                          [ x2' ]           [Q2]
  31  *
  32  * The equalities then become
  33  *
  34  *              H1 x1' - c = 0   or   x1' = H1^{-1} c = c'
  35  *
  36  * If any of the c' is non-integer, then the original set has no
  37  * integer solutions (since the x' are a unimodular transformation
  38  * of the x).
  39  * Otherwise, the transformation is given by
  40  *
  41  *              x = U1 H1^{-1} c + U2 x2'
  42  *
  43  * The inverse transformation is simply
  44  *
  45  *              x2' = Q2 x
  46  */
  47 struct isl_mat *isl_mat_variable_compression(struct isl_ctx *ctx,
  48                         struct isl_mat *B, struct isl_mat **T2)
  49 {
  50         int i;
  51         struct isl_mat *H = NULL, *C = NULL, *H1, *U = NULL, *U1, *U2, *TC;
  52         unsigned dim;
  53
  54         if (T2)
  55                 *T2 = NULL;
  56         if (!B)
  57                 goto error;
  58
  59         dim = B->n_col - 1;
  60         H = isl_mat_sub_alloc(ctx, B->row, 0, B->n_row, 1, dim);
  61         H = isl_mat_left_hermite(ctx, H, 0, &U, T2);
  62         if (!H || !U || (T2 && !*T2))
  63                 goto error;
  64         if (T2) {
  65                 *T2 = isl_mat_drop_rows(ctx, *T2, 0, B->n_row);
  66                 *T2 = isl_mat_lin_to_aff(ctx, *T2);
  67                 if (!*T2)
  68                         goto error;
  69         }
  70         C = isl_mat_alloc(ctx, 1+B->n_row, 1);
  71         if (!C)
  72                 goto error;
  73         isl_int_set_si(C->row[0][0], 1);
  74         isl_mat_sub_neg(ctx, C->row+1, B->row, B->n_row, 0, 0, 1);
  75         H1 = isl_mat_sub_alloc(ctx, H->row, 0, H->n_row, 0, H->n_row);
  76         H1 = isl_mat_lin_to_aff(ctx, H1);
  77         TC = isl_mat_inverse_product(ctx, H1, C);
  78         if (!TC)
  79                 goto error;
  80         isl_mat_free(ctx, H);
  81         if (!isl_int_is_one(TC->row[0][0])) {
  82                 for (i = 0; i < B->n_row; ++i) {
  83                         if (!isl_int_is_divisible_by(TC->row[1+i][0], TC->row[0][0])) {
  84                                 isl_mat_free(ctx, B);
  85                                 isl_mat_free(ctx, TC);
  86                                 isl_mat_free(ctx, U);
  87                                 if (T2) {
  88                                         isl_mat_free(ctx, *T2);
  89                                         *T2 = NULL;
  90                                 }
  91                                 return isl_mat_alloc(ctx, 1 + B->n_col, 0);
  92                         }
  93                         isl_seq_scale_down(TC->row[1+i], TC->row[1+i], TC->row[0][0], 1);
  94                 }
  95                 isl_int_set_si(TC->row[0][0], 1);
  96         }
  97         U1 = isl_mat_sub_alloc(ctx, U->row, 0, U->n_row, 0, B->n_row);
  98         U1 = isl_mat_lin_to_aff(ctx, U1);
  99         U2 = isl_mat_sub_alloc(ctx, U->row, 0, U->n_row,
 100                                 B->n_row, U->n_row - B->n_row);
 101         U2 = isl_mat_lin_to_aff(ctx, U2);
 102         isl_mat_free(ctx, U);
 103         TC = isl_mat_product(ctx, U1, TC);
 104         TC = isl_mat_aff_direct_sum(ctx, TC, U2);
 105
 106         isl_mat_free(ctx, B);
 107
 108         return TC;
 109 error:
 110         isl_mat_free(ctx, B);
 111         isl_mat_free(ctx, H);
 112         isl_mat_free(ctx, U);
 113         if (T2) {
 114                 isl_mat_free(ctx, *T2);
 115                 *T2 = NULL;
 116         }
 117         return NULL;
 118 }
 119
 120 /* Use the n equalities of bset to unimodularly transform the
 121  * variables x such that n transformed variables x1' have a constant value
 122  * and rewrite the constraints of bset in terms of the remaining
 123  * transformed variables x2'.  The matrix pointed to by T maps
 124  * the new variables x2' back to the original variables x, while T2
 125  * maps the original variables to the new variables.
 126  */
 127 static struct isl_basic_set *compress_variables(struct isl_ctx *ctx,
 128         struct isl_basic_set *bset, struct isl_mat **T, struct isl_mat **T2)
 129 {
 130         struct isl_mat *B, *TC;
 131         unsigned dim;
 132
 133         if (T)
 134                 *T = NULL;
 135         if (T2)
 136                 *T2 = NULL;
 137         if (!bset)
 138                 goto error;
 139         isl_assert(ctx, isl_basic_set_n_param(bset) == 0, goto error);
 140         isl_assert(ctx, bset->n_div == 0, goto error);
 141         dim = isl_basic_set_n_dim(bset);
 142         isl_assert(ctx, bset->n_eq <= dim, goto error);
 143         if (bset->n_eq == 0)
 144                 return bset;
 145
 146         B = isl_mat_sub_alloc(ctx, bset->eq, 0, bset->n_eq, 0, 1 + dim);
 147         TC = isl_mat_variable_compression(ctx, B, T2);
 148         if (!TC)
 149                 goto error;
 150         if (TC->n_col == 0) {
 151                 isl_mat_free(ctx, TC);
 152                 if (T2) {
 153                         isl_mat_free(ctx, *T2);
 154                         *T2 = NULL;
 155                 }
 156                 return isl_basic_set_set_to_empty(bset);
 157         }
 158
 159         bset = isl_basic_set_preimage(ctx, bset, T ? isl_mat_copy(ctx, TC) : TC);
 160         if (T)
 161                 *T = TC;
 162         return bset;
 163 error:
 164         isl_basic_set_free(bset);
 165         return NULL;
 166 }
 167
 168 struct isl_basic_set *isl_basic_set_remove_equalities(
 169         struct isl_basic_set *bset, struct isl_mat **T, struct isl_mat **T2)
 170 {
 171         if (T)
 172                 *T = NULL;
 173         if (T2)
 174                 *T2 = NULL;
 175         if (!bset)
 176                 return NULL;
 177         isl_assert(bset->ctx, isl_basic_set_n_param(bset) == 0, goto error);
 178         bset = isl_basic_set_gauss(bset, NULL);
 179         if (F_ISSET(bset, ISL_BASIC_SET_EMPTY))
 180                 return bset;
 181         bset = compress_variables(bset->ctx, bset, T, T2);
 182         return bset;
 183 error:
 184         isl_basic_set_free(bset);
 185         *T = NULL;
 186         return NULL;
 187 }
 188
 189 /* Check if dimension dim belongs to a residue class
 190  *              i_dim \equiv r mod m
 191  * with m != 1 and if so return m in *modulo and r in *residue.
 192  */
 193 int isl_basic_set_dim_residue_class(struct isl_basic_set *bset,
 194         int pos, isl_int *modulo, isl_int *residue)
 195 {
 196         struct isl_ctx *ctx;
 197         struct isl_mat *H = NULL, *U = NULL, *C, *H1, *U1;
 198         unsigned total;
 199         unsigned nparam;
 200
 201         if (!bset || !modulo || !residue)
 202                 return -1;
 203
 204         ctx = bset->ctx;
 205         total = isl_basic_set_total_dim(bset);
 206         nparam = isl_basic_set_n_param(bset);
 207         H = isl_mat_sub_alloc(ctx, bset->eq, 0, bset->n_eq, 1, total);
 208         H = isl_mat_left_hermite(ctx, H, 0, &U, NULL);
 209         if (!H)
 210                 return -1;
 211
 212         isl_seq_gcd(U->row[nparam + pos]+bset->n_eq,
 213                         total-bset->n_eq, modulo);
 214         if (isl_int_is_zero(*modulo) || isl_int_is_one(*modulo)) {
 215                 isl_int_set_si(*residue, 0);
 216                 isl_mat_free(ctx, H);
 217                 isl_mat_free(ctx, U);
 218                 return 0;
 219         }
 220
 221         C = isl_mat_alloc(ctx, 1+bset->n_eq, 1);
 222         if (!C)
 223                 goto error;
 224         isl_int_set_si(C->row[0][0], 1);
 225         isl_mat_sub_neg(ctx, C->row+1, bset->eq, bset->n_eq, 0, 0, 1);
 226         H1 = isl_mat_sub_alloc(ctx, H->row, 0, H->n_row, 0, H->n_row);
 227         H1 = isl_mat_lin_to_aff(ctx, H1);
 228         C = isl_mat_inverse_product(ctx, H1, C);
 229         isl_mat_free(ctx, H);
 230         U1 = isl_mat_sub_alloc(ctx, U->row, nparam+pos, 1, 0, bset->n_eq);
 231         U1 = isl_mat_lin_to_aff(ctx, U1);
 232         isl_mat_free(ctx, U);
 233         C = isl_mat_product(ctx, U1, C);
 234         if (!C)
 235                 goto error;
 236         if (!isl_int_is_divisible_by(C->row[1][0], C->row[0][0])) {
 237                 bset = isl_basic_set_copy(bset);
 238                 bset = isl_basic_set_set_to_empty(bset);
 239                 isl_basic_set_free(bset);
 240                 isl_int_set_si(*modulo, 0);
 241                 isl_int_set_si(*residue, 0);
 242                 return 0;
 243         }
 244         isl_int_divexact(*residue, C->row[1][0], C->row[0][0]);
 245         isl_int_fdiv_r(*residue, *residue, *modulo);
 246         isl_mat_free(ctx, C);
 247         return 0;
 248 error:
 249         isl_mat_free(ctx, H);
 250         isl_mat_free(ctx, U);
 251         return -1;
 252 }