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
  51         if (T)
  52                 *T = NULL;
  53         if (T2)
  54                 *T2 = NULL;
  55         if (!bset)
  56                 goto error;
  57         isl_assert(ctx, bset->nparam == 0, goto error);
  58         isl_assert(ctx, bset->n_div == 0, goto error);
  59         isl_assert(ctx, bset->n_eq <= bset->dim, goto error);
  60         if (bset->n_eq == 0)
  61                 return bset;
  62
  63         H = isl_mat_sub_alloc(ctx, bset->eq, 0, bset->n_eq, 1, bset->dim);
  64         H = isl_mat_left_hermite(ctx, H, 0, &U, T2);
  65         if (!H || !U || (T2 && !*T2))
  66                 goto error;
  67         if (T2) {
  68                 *T2 = isl_mat_drop_rows(ctx, *T2, 0, bset->n_eq);
  69                 *T2 = isl_mat_lin_to_aff(ctx, *T2);
  70                 if (!*T2)
  71                         goto error;
  72         }
  73         C = isl_mat_alloc(ctx, 1+bset->n_eq, 1);
  74         if (!C)
  75                 goto error;
  76         isl_int_set_si(C->row[0][0], 1);
  77         isl_mat_sub_neg(ctx, C->row+1, bset->eq, bset->n_eq, 0, 0, 1);
  78         H1 = isl_mat_sub_alloc(ctx, H->row, 0, H->n_row, 0, H->n_row);
  79         H1 = isl_mat_lin_to_aff(ctx, H1);
  80         TC = isl_mat_inverse_product(ctx, H1, C);
  81         if (!TC)
  82                 goto error;
  83         isl_mat_free(ctx, H);
  84         if (!isl_int_is_one(TC->row[0][0])) {
  85                 for (i = 0; i < bset->n_eq; ++i) {
  86                         if (!isl_int_is_divisible_by(TC->row[1+i][0], TC->row[0][0])) {
  87                                 isl_mat_free(ctx, TC);
  88                                 isl_mat_free(ctx, U);
  89                                 if (T2) {
  90                                         isl_mat_free(ctx, *T2);
  91                                         *T2 = NULL;
  92                                 }
  93                                 return isl_basic_set_set_to_empty(bset);
  94                         }
  95                         isl_seq_scale_down(TC->row[1+i], TC->row[1+i], TC->row[0][0], 1);
  96                 }
  97                 isl_int_set_si(TC->row[0][0], 1);
  98         }
  99         U1 = isl_mat_sub_alloc(ctx, U->row, 0, U->n_row, 0, bset->n_eq);
 100         U1 = isl_mat_lin_to_aff(ctx, U1);
 101         U2 = isl_mat_sub_alloc(ctx, U->row, 0, U->n_row,
 102                                 bset->n_eq, U->n_row - bset->n_eq);
 103         U2 = isl_mat_lin_to_aff(ctx, U2);
 104         isl_mat_free(ctx, U);
 105         TC = isl_mat_product(ctx, U1, TC);
 106         TC = isl_mat_aff_direct_sum(ctx, TC, U2);
 107         bset = isl_basic_set_preimage(ctx, bset, T ? isl_mat_copy(ctx, TC) : TC);
 108         if (T)
 109                 *T = TC;
 110         return bset;
 111 error:
 112         isl_mat_free(ctx, H);
 113         isl_mat_free(ctx, U);
 114         if (T2)
 115                 isl_mat_free(ctx, *T2);
 116         isl_basic_set_free(bset);
 117         if (T)
 118                 *T = NULL;
 119         if (T2)
 120                 *T2 = NULL;
 121         return NULL;
 122 }
 123
 124 struct isl_basic_set *isl_basic_set_remove_equalities(
 125         struct isl_basic_set *bset, struct isl_mat **T, struct isl_mat **T2)
 126 {
 127         if (T)
 128                 *T = NULL;
 129         if (T2)
 130                 *T2 = NULL;
 131         if (!bset)
 132                 return NULL;
 133         isl_assert(bset->ctx, bset->nparam == 0, goto error);
 134         bset = isl_basic_set_gauss(bset, NULL);
 135         if (F_ISSET(bset, ISL_BASIC_SET_EMPTY))
 136                 return bset;
 137         bset = compress_variables(bset->ctx, bset, T, T2);
 138         return bset;
 139 error:
 140         isl_basic_set_free(bset);
 141         *T = NULL;
 142         return NULL;
 143 }
 144
 145 /* Check if dimension dim belongs to a residue class
 146  *              i_dim \equiv r mod m
 147  * with m != 1 and if so return m in *modulo and r in *residue.
 148  */
 149 int isl_basic_set_dim_residue_class(struct isl_basic_set *bset,
 150         int pos, isl_int *modulo, isl_int *residue)
 151 {
 152         struct isl_ctx *ctx;
 153         struct isl_mat *H = NULL, *U = NULL, *C, *H1, *U1;
 154         unsigned total;
 155
 156         if (!bset || !modulo || !residue)
 157                 return -1;
 158
 159         ctx = bset->ctx;
 160         total = bset->nparam + bset->dim + bset->n_div;
 161         H = isl_mat_sub_alloc(ctx, bset->eq, 0, bset->n_eq, 1, total);
 162         H = isl_mat_left_hermite(ctx, H, 0, &U, NULL);
 163         if (!H)
 164                 return -1;
 165
 166         isl_seq_gcd(U->row[bset->nparam + pos]+bset->n_eq,
 167                         total-bset->n_eq, modulo);
 168         if (isl_int_is_zero(*modulo) || isl_int_is_one(*modulo)) {
 169                 isl_int_set_si(*residue, 0);
 170                 isl_mat_free(ctx, H);
 171                 isl_mat_free(ctx, U);
 172                 return 0;
 173         }
 174
 175         C = isl_mat_alloc(ctx, 1+bset->n_eq, 1);
 176         if (!C)
 177                 goto error;
 178         isl_int_set_si(C->row[0][0], 1);
 179         isl_mat_sub_neg(ctx, C->row+1, bset->eq, bset->n_eq, 0, 0, 1);
 180         H1 = isl_mat_sub_alloc(ctx, H->row, 0, H->n_row, 0, H->n_row);
 181         H1 = isl_mat_lin_to_aff(ctx, H1);
 182         C = isl_mat_inverse_product(ctx, H1, C);
 183         isl_mat_free(ctx, H);
 184         U1 = isl_mat_sub_alloc(ctx, U->row, bset->nparam+pos, 1, 0, bset->n_eq);
 185         U1 = isl_mat_lin_to_aff(ctx, U1);
 186         isl_mat_free(ctx, U);
 187         C = isl_mat_product(ctx, U1, C);
 188         if (!C)
 189                 goto error;
 190         if (!isl_int_is_divisible_by(C->row[1][0], C->row[0][0])) {
 191                 bset = isl_basic_set_copy(bset);
 192                 bset = isl_basic_set_set_to_empty(bset);
 193                 isl_basic_set_free(bset);
 194                 isl_int_set_si(*modulo, 0);
 195                 isl_int_set_si(*residue, 0);
 196                 return 0;
 197         }
 198         isl_int_divexact(*residue, C->row[1][0], C->row[0][0]);
 199         isl_int_fdiv_r(*residue, *residue, *modulo);
 200         isl_mat_free(ctx, C);
 201         return 0;
 202 error:
 203         isl_mat_free(ctx, H);
 204         isl_mat_free(ctx, U);
 205         return -1;
 206 }