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 /* 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.
  12  *
  13  * Let the equalities of bset be
  14  *
  15  *              M x - c = 0
  16  *
  17  * Compute the (left) Hermite normal form of M,
  18  *
  19  *              M [U1 U2] = M U = H = [H1 0]
  20  * or
  21  *                            M = H Q = [H1 0] [Q1]
  22  *                                             [Q2]
  23  *
  24  * with U, Q unimodular, Q = U^{-1} (and H lower triangular).
  25  * Define the transformed variables as
  26  *
  27  *              x = [U1 U2] [ x1' ] = [U1 U2] [Q1] x
  28  *                          [ x2' ]           [Q2]
  29  *
  30  * The equalities then become
  31  *
  32  *              H1 x1' - c = 0   or   x1' = H1^{-1} c = c'
  33  *
  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
  36  * of the x).
  37  * Otherwise, the transformation is given by
  38  *
  39  *              x = U1 H1^{-1} c + U2 x2'
  40  *
  41  * The inverse transformation is simply
  42  *
  43  *              x2' = Q2 x
  44  */
  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)
  47 {
  48         int i;
  49         struct isl_mat *H = NULL, *C = NULL, *H1, *U = NULL, *U1, *U2, *TC;
  50         unsigned dim;
  51
  52         if (T)
  53                 *T = NULL;
  54         if (T2)
  55                 *T2 = NULL;
  56         if (!bset)
  57                 goto error;
  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);
  62         if (bset->n_eq == 0)
  63                 return bset;
  64
  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))
  68                 goto error;
  69         if (T2) {
  70                 *T2 = isl_mat_drop_rows(ctx, *T2, 0, bset->n_eq);
  71                 *T2 = isl_mat_lin_to_aff(ctx, *T2);
  72                 if (!*T2)
  73                         goto error;
  74         }
  75         C = isl_mat_alloc(ctx, 1+bset->n_eq, 1);
  76         if (!C)
  77                 goto error;
  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);
  83         if (!TC)
  84                 goto error;
  85         isl_mat_free(ctx, H);
  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);
  90                                 isl_mat_free(ctx, U);
  91                                 if (T2) {
  92                                         isl_mat_free(ctx, *T2);
  93                                         *T2 = NULL;
  94                                 }
  95                                 return isl_basic_set_set_to_empty(bset);
  96                         }
  97                         isl_seq_scale_down(TC->row[1+i], TC->row[1+i], TC->row[0][0], 1);
  98                 }
  99                 isl_int_set_si(TC->row[0][0], 1);
 100         }
 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);
 110         if (T)
 111                 *T = TC;
 112         return bset;
 113 error:
 114         isl_mat_free(ctx, H);
 115         isl_mat_free(ctx, U);
 116         if (T2)
 117                 isl_mat_free(ctx, *T2);
 118         isl_basic_set_free(bset);
 119         if (T)
 120                 *T = NULL;
 121         if (T2)
 122                 *T2 = NULL;
 123         return NULL;
 124 }
 125
 126 struct isl_basic_set *isl_basic_set_remove_equalities(
 127         struct isl_basic_set *bset, struct isl_mat **T, struct isl_mat **T2)
 128 {
 129         if (T)
 130                 *T = NULL;
 131         if (T2)
 132                 *T2 = NULL;
 133         if (!bset)
 134                 return NULL;
 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))
 138                 return bset;
 139         bset = compress_variables(bset->ctx, bset, T, T2);
 140         return bset;
 141 error:
 142         isl_basic_set_free(bset);
 143         *T = NULL;
 144         return NULL;
 145 }
 146
 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.
 150  */
 151 int isl_basic_set_dim_residue_class(struct isl_basic_set *bset,
 152         int pos, isl_int *modulo, isl_int *residue)
 153 {
 154         struct isl_ctx *ctx;
 155         struct isl_mat *H = NULL, *U = NULL, *C, *H1, *U1;
 156         unsigned total;
 157         unsigned nparam;
 158
 159         if (!bset || !modulo || !residue)
 160                 return -1;
 161
 162         ctx = bset->ctx;
 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);
 167         if (!H)
 168                 return -1;
 169
 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);
 176                 return 0;
 177         }
 178
 179         C = isl_mat_alloc(ctx, 1+bset->n_eq, 1);
 180         if (!C)
 181                 goto error;
 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);
 192         if (!C)
 193                 goto error;
 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);
 200                 return 0;
 201         }
 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);
 205         return 0;
 206 error:
 207         isl_mat_free(ctx, H);
 208         isl_mat_free(ctx, U);
 209         return -1;
 210 }