isl_farkas.c

   1 /*
   2  * Copyright 2010      INRIA Saclay
   3  *
   4  * Use of this software is governed by the GNU LGPLv2.1 license
   5  *
   6  * Written by Sven Verdoolaege, INRIA Saclay - Ile-de-France,
   7  * Parc Club Orsay Universite, ZAC des vignes, 4 rue Jacques Monod,
   8  * 91893 Orsay, France
   9  */
  10
  11 #include <isl_map_private.h>
  12 #include <isl/set.h>
  13 #include <isl_space_private.h>
  14 #include <isl/seq.h>
  15
  16 /*
  17  * Let C be a cone and define
  18  *
  19  *      C' := { y | forall x in C : y x >= 0 }
  20  *
  21  * C' contains the coefficients of all linear constraints
  22  * that are valid for C.
  23  * Furthermore, C'' = C.
  24  *
  25  * If C is defined as { x | A x >= 0 }
  26  * then any element in C' must be a non-negative combination
  27  * of the rows of A, i.e., y = t A with t >= 0.  That is,
  28  *
  29  *      C' = { y | exists t >= 0 : y = t A }
  30  *
  31  * If any of the rows in A actually represents an equality, then
  32  * also negative combinations of this row are allowed and so the
  33  * non-negativity constraint on the corresponding element of t
  34  * can be dropped.
  35  *
  36  * A polyhedron P = { x | b + A x >= 0 } can be represented
  37  * in homogeneous coordinates by the cone
  38  * C = { [z,x] | b z + A x >= and z >= 0 }
  39  * The valid linear constraints on C correspond to the valid affine
  40  * constraints on P.
  41  * This is essentially Farkas' lemma.
  42  *
  43  * Let A' = [b A], then, since
  44  *                                [ 1 0 ]
  45  *              [ w y ] = [t_0 t] [ b A ]
  46  *
  47  * we have
  48  *
  49  *      C' = { w, y | exists t_0, t >= 0 : y = t A' and w = t_0 + t b }
  50  * or
  51  *
  52  *      C' = { w, y | exists t >= 0 : y = t A' and w - t b >= 0 }
  53  *
  54  * In practice, we introduce an extra variable (w), shifting all
  55  * other variables to the right, and an extra inequality
  56  * (w - t b >= 0) corresponding to the positivity constraint on
  57  * the homogeneous coordinate.
  58  *
  59  * When going back from coefficients to solutions, we immediately
  60  * plug in 1 for z, which corresponds to shifting all variables
  61  * to the left, with the leftmost ending up in the constant position.
  62  */
  63
  64 /* Add the given prefix to all named isl_dim_set dimensions in "dim".
  65  */
  66 static __isl_give isl_space *isl_space_prefix(__isl_take isl_space *dim,
  67         const char *prefix)
  68 {
  69         int i;
  70         isl_ctx *ctx;
  71         unsigned nvar;
  72         size_t prefix_len = strlen(prefix);
  73
  74         if (!dim)
  75                 return NULL;
  76
  77         ctx = isl_space_get_ctx(dim);
  78         nvar = isl_space_dim(dim, isl_dim_set);
  79
  80         for (i = 0; i < nvar; ++i) {
  81                 const char *name;
  82                 char *prefix_name;
  83
  84                 name = isl_space_get_dim_name(dim, isl_dim_set, i);
  85                 if (!name)
  86                         continue;
  87
  88                 prefix_name = isl_alloc_array(ctx, char,
  89                                               prefix_len + strlen(name) + 1);
  90                 if (!prefix_name)
  91                         goto error;
  92                 memcpy(prefix_name, prefix, prefix_len);
  93                 strcpy(prefix_name + prefix_len, name);
  94
  95                 dim = isl_space_set_dim_name(dim, isl_dim_set, i, prefix_name);
  96                 free(prefix_name);
  97         }
  98
  99         return dim;
 100 error:
 101         isl_space_free(dim);
 102         return NULL;
 103 }
 104
 105 /* Given a dimension specification of the solutions space, construct
 106  * a dimension specification for the space of coefficients.
 107  *
 108  * In particular transform
 109  *
 110  *      [params] -> { S }
 111  *
 112  * to
 113  *
 114  *      { coefficients[[cst, params] -> S] }
 115  *
 116  * and prefix each dimension name with "c_".
 117  */
 118 static __isl_give isl_space *isl_space_coefficients(__isl_take isl_space *dim)
 119 {
 120         isl_space *dim_param;
 121         unsigned nvar;
 122         unsigned nparam;
 123
 124         nvar = isl_space_dim(dim, isl_dim_set);
 125         nparam = isl_space_dim(dim, isl_dim_param);
 126         dim_param = isl_space_copy(dim);
 127         dim_param = isl_space_drop_dims(dim_param, isl_dim_set, 0, nvar);
 128         dim_param = isl_space_move_dims(dim_param, isl_dim_set, 0,
 129                                  isl_dim_param, 0, nparam);
 130         dim_param = isl_space_prefix(dim_param, "c_");
 131         dim_param = isl_space_insert_dims(dim_param, isl_dim_set, 0, 1);
 132         dim_param = isl_space_set_dim_name(dim_param, isl_dim_set, 0, "c_cst");
 133         dim = isl_space_drop_dims(dim, isl_dim_param, 0, nparam);
 134         dim = isl_space_prefix(dim, "c_");
 135         dim = isl_space_join(isl_space_from_domain(dim_param),
 136                            isl_space_from_range(dim));
 137         dim = isl_space_wrap(dim);
 138         dim = isl_space_set_tuple_name(dim, isl_dim_set, "coefficients");
 139
 140         return dim;
 141 }
 142
 143 /* Drop the given prefix from all named dimensions of type "type" in "dim".
 144  */
 145 static __isl_give isl_space *isl_space_unprefix(__isl_take isl_space *dim,
 146         enum isl_dim_type type, const char *prefix)
 147 {
 148         int i;
 149         unsigned n;
 150         size_t prefix_len = strlen(prefix);
 151
 152         n = isl_space_dim(dim, type);
 153
 154         for (i = 0; i < n; ++i) {
 155                 const char *name;
 156
 157                 name = isl_space_get_dim_name(dim, type, i);
 158                 if (!name)
 159                         continue;
 160                 if (strncmp(name, prefix, prefix_len))
 161                         continue;
 162
 163                 dim = isl_space_set_dim_name(dim, type, i, name + prefix_len);
 164         }
 165
 166         return dim;
 167 }
 168
 169 /* Given a dimension specification of the space of coefficients, construct
 170  * a dimension specification for the space of solutions.
 171  *
 172  * In particular transform
 173  *
 174  *      { coefficients[[cst, params] -> S] }
 175  *
 176  * to
 177  *
 178  *      [params] -> { S }
 179  *
 180  * and drop the "c_" prefix from the dimension names.
 181  */
 182 static __isl_give isl_space *isl_space_solutions(__isl_take isl_space *dim)
 183 {
 184         unsigned nparam;
 185
 186         dim = isl_space_unwrap(dim);
 187         dim = isl_space_drop_dims(dim, isl_dim_in, 0, 1);
 188         dim = isl_space_unprefix(dim, isl_dim_in, "c_");
 189         dim = isl_space_unprefix(dim, isl_dim_out, "c_");
 190         nparam = isl_space_dim(dim, isl_dim_in);
 191         dim = isl_space_move_dims(dim, isl_dim_param, 0, isl_dim_in, 0, nparam);
 192         dim = isl_space_range(dim);
 193
 194         return dim;
 195 }
 196
 197 /* Compute the dual of "bset" by applying Farkas' lemma.
 198  * As explained above, we add an extra dimension to represent
 199  * the coefficient of the constant term when going from solutions
 200  * to coefficients (shift == 1) and we drop the extra dimension when going
 201  * in the opposite direction (shift == -1).  "dim" is the space in which
 202  * the dual should be created.
 203  */
 204 static __isl_give isl_basic_set *farkas(__isl_take isl_space *dim,
 205         __isl_take isl_basic_set *bset, int shift)
 206 {
 207         int i, j, k;
 208         isl_basic_set *dual = NULL;
 209         unsigned total;
 210
 211         total = isl_basic_set_total_dim(bset);
 212
 213         dual = isl_basic_set_alloc_space(dim, bset->n_eq + bset->n_ineq,
 214                                         total, bset->n_ineq + (shift > 0));
 215         dual = isl_basic_set_set_rational(dual);
 216
 217         for (i = 0; i < bset->n_eq + bset->n_ineq; ++i) {
 218                 k = isl_basic_set_alloc_div(dual);
 219                 if (k < 0)
 220                         goto error;
 221                 isl_int_set_si(dual->div[k][0], 0);
 222         }
 223
 224         for (i = 0; i < total; ++i) {
 225                 k = isl_basic_set_alloc_equality(dual);
 226                 if (k < 0)
 227                         goto error;
 228                 isl_seq_clr(dual->eq[k], 1 + shift + total);
 229                 isl_int_set_si(dual->eq[k][1 + shift + i], -1);
 230                 for (j = 0; j < bset->n_eq; ++j)
 231                         isl_int_set(dual->eq[k][1 + shift + total + j],
 232                                     bset->eq[j][1 + i]);
 233                 for (j = 0; j < bset->n_ineq; ++j)
 234                         isl_int_set(dual->eq[k][1 + shift + total + bset->n_eq + j],
 235                                     bset->ineq[j][1 + i]);
 236         }
 237
 238         for (i = 0; i < bset->n_ineq; ++i) {
 239                 k = isl_basic_set_alloc_inequality(dual);
 240                 if (k < 0)
 241                         goto error;
 242                 isl_seq_clr(dual->ineq[k],
 243                             1 + shift + total + bset->n_eq + bset->n_ineq);
 244                 isl_int_set_si(dual->ineq[k][1 + shift + total + bset->n_eq + i], 1);
 245         }
 246
 247         if (shift > 0) {
 248                 k = isl_basic_set_alloc_inequality(dual);
 249                 if (k < 0)
 250                         goto error;
 251                 isl_seq_clr(dual->ineq[k], 2 + total);
 252                 isl_int_set_si(dual->ineq[k][1], 1);
 253                 for (j = 0; j < bset->n_eq; ++j)
 254                         isl_int_neg(dual->ineq[k][2 + total + j],
 255                                     bset->eq[j][0]);
 256                 for (j = 0; j < bset->n_ineq; ++j)
 257                         isl_int_neg(dual->ineq[k][2 + total + bset->n_eq + j],
 258                                     bset->ineq[j][0]);
 259         }
 260
 261         dual = isl_basic_set_remove_divs(dual);
 262         isl_basic_set_simplify(dual);
 263         isl_basic_set_finalize(dual);
 264
 265         isl_basic_set_free(bset);
 266         return dual;
 267 error:
 268         isl_basic_set_free(bset);
 269         isl_basic_set_free(dual);
 270         return NULL;
 271 }
 272
 273 /* Construct a basic set containing the tuples of coefficients of all
 274  * valid affine constraints on the given basic set.
 275  */
 276 __isl_give isl_basic_set *isl_basic_set_coefficients(
 277         __isl_take isl_basic_set *bset)
 278 {
 279         isl_space *dim;
 280
 281         if (!bset)
 282                 return NULL;
 283         if (bset->n_div)
 284                 isl_die(bset->ctx, isl_error_invalid,
 285                         "input set not allowed to have local variables",
 286                         goto error);
 287
 288         dim = isl_basic_set_get_space(bset);
 289         dim = isl_space_coefficients(dim);
 290
 291         return farkas(dim, bset, 1);
 292 error:
 293         isl_basic_set_free(bset);
 294         return NULL;
 295 }
 296
 297 /* Construct a basic set containing the elements that satisfy all
 298  * affine constraints whose coefficient tuples are
 299  * contained in the given basic set.
 300  */
 301 __isl_give isl_basic_set *isl_basic_set_solutions(
 302         __isl_take isl_basic_set *bset)
 303 {
 304         isl_space *dim;
 305
 306         if (!bset)
 307                 return NULL;
 308         if (bset->n_div)
 309                 isl_die(bset->ctx, isl_error_invalid,
 310                         "input set not allowed to have local variables",
 311                         goto error);
 312
 313         dim = isl_basic_set_get_space(bset);
 314         dim = isl_space_solutions(dim);
 315
 316         return farkas(dim, bset, -1);
 317 error:
 318         isl_basic_set_free(bset);
 319         return NULL;
 320 }
 321
 322 /* Construct a basic set containing the tuples of coefficients of all
 323  * valid affine constraints on the given set.
 324  */
 325 __isl_give isl_basic_set *isl_set_coefficients(__isl_take isl_set *set)
 326 {
 327         int i;
 328         isl_basic_set *coeff;
 329
 330         if (!set)
 331                 return NULL;
 332         if (set->n == 0) {
 333                 isl_space *dim = isl_set_get_space(set);
 334                 dim = isl_space_coefficients(dim);
 335                 coeff = isl_basic_set_universe(dim);
 336                 coeff = isl_basic_set_set_rational(coeff);
 337                 isl_set_free(set);
 338                 return coeff;
 339         }
 340
 341         coeff = isl_basic_set_coefficients(isl_basic_set_copy(set->p[0]));
 342
 343         for (i = 1; i < set->n; ++i) {
 344                 isl_basic_set *bset, *coeff_i;
 345                 bset = isl_basic_set_copy(set->p[i]);
 346                 coeff_i = isl_basic_set_coefficients(bset);
 347                 coeff = isl_basic_set_intersect(coeff, coeff_i);
 348         }
 349
 350         isl_set_free(set);
 351         return coeff;
 352 }
 353
 354 /* Construct a basic set containing the elements that satisfy all
 355  * affine constraints whose coefficient tuples are
 356  * contained in the given set.
 357  */
 358 __isl_give isl_basic_set *isl_set_solutions(__isl_take isl_set *set)
 359 {
 360         int i;
 361         isl_basic_set *sol;
 362
 363         if (!set)
 364                 return NULL;
 365         if (set->n == 0) {
 366                 isl_space *dim = isl_set_get_space(set);
 367                 dim = isl_space_solutions(dim);
 368                 sol = isl_basic_set_universe(dim);
 369                 sol = isl_basic_set_set_rational(sol);
 370                 isl_set_free(set);
 371                 return sol;
 372         }
 373
 374         sol = isl_basic_set_solutions(isl_basic_set_copy(set->p[0]));
 375
 376         for (i = 1; i < set->n; ++i) {
 377                 isl_basic_set *bset, *sol_i;
 378                 bset = isl_basic_set_copy(set->p[i]);
 379                 sol_i = isl_basic_set_solutions(bset);
 380                 sol = isl_basic_set_intersect(sol, sol_i);
 381         }
 382
 383         isl_set_free(set);
 384         return sol;
 385 }