isl_farkas.c

   1 /*
   2  * Copyright 2010      INRIA Saclay
   3  *
   4  * Use of this software is governed by the MIT 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  * 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 "space".
  65  */
  66 static __isl_give isl_space *isl_space_prefix(__isl_take isl_space *space,
  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 (!space)
  75                 return NULL;
  76
  77         ctx = isl_space_get_ctx(space);
  78         nvar = isl_space_dim(space, 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(space, 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                 space = isl_space_set_dim_name(space,
  96                                                 isl_dim_set, i, prefix_name);
  97                 free(prefix_name);
  98         }
  99
 100         return space;
 101 error:
 102         isl_space_free(space);
 103         return NULL;
 104 }
 105
 106 /* Given a dimension specification of the solutions space, construct
 107  * a dimension specification for the space of coefficients.
 108  *
 109  * In particular transform
 110  *
 111  *      [params] -> { S }
 112  *
 113  * to
 114  *
 115  *      { coefficients[[cst, params] -> S] }
 116  *
 117  * and prefix each dimension name with "c_".
 118  */
 119 static __isl_give isl_space *isl_space_coefficients(__isl_take isl_space *space)
 120 {
 121         isl_space *space_param;
 122         unsigned nvar;
 123         unsigned nparam;
 124
 125         nvar = isl_space_dim(space, isl_dim_set);
 126         nparam = isl_space_dim(space, isl_dim_param);
 127         space_param = isl_space_copy(space);
 128         space_param = isl_space_drop_dims(space_param, isl_dim_set, 0, nvar);
 129         space_param = isl_space_move_dims(space_param, isl_dim_set, 0,
 130                                  isl_dim_param, 0, nparam);
 131         space_param = isl_space_prefix(space_param, "c_");
 132         space_param = isl_space_insert_dims(space_param, isl_dim_set, 0, 1);
 133         space_param = isl_space_set_dim_name(space_param,
 134                                 isl_dim_set, 0, "c_cst");
 135         space = isl_space_drop_dims(space, isl_dim_param, 0, nparam);
 136         space = isl_space_prefix(space, "c_");
 137         space = isl_space_join(isl_space_from_domain(space_param),
 138                            isl_space_from_range(space));
 139         space = isl_space_wrap(space);
 140         space = isl_space_set_tuple_name(space, isl_dim_set, "coefficients");
 141
 142         return space;
 143 }
 144
 145 /* Drop the given prefix from all named dimensions of type "type" in "space".
 146  */
 147 static __isl_give isl_space *isl_space_unprefix(__isl_take isl_space *space,
 148         enum isl_dim_type type, const char *prefix)
 149 {
 150         int i;
 151         unsigned n;
 152         size_t prefix_len = strlen(prefix);
 153
 154         n = isl_space_dim(space, type);
 155
 156         for (i = 0; i < n; ++i) {
 157                 const char *name;
 158
 159                 name = isl_space_get_dim_name(space, type, i);
 160                 if (!name)
 161                         continue;
 162                 if (strncmp(name, prefix, prefix_len))
 163                         continue;
 164
 165                 space = isl_space_set_dim_name(space,
 166                                                 type, i, name + prefix_len);
 167         }
 168
 169         return space;
 170 }
 171
 172 /* Given a dimension specification of the space of coefficients, construct
 173  * a dimension specification for the space of solutions.
 174  *
 175  * In particular transform
 176  *
 177  *      { coefficients[[cst, params] -> S] }
 178  *
 179  * to
 180  *
 181  *      [params] -> { S }
 182  *
 183  * and drop the "c_" prefix from the dimension names.
 184  */
 185 static __isl_give isl_space *isl_space_solutions(__isl_take isl_space *space)
 186 {
 187         unsigned nparam;
 188
 189         space = isl_space_unwrap(space);
 190         space = isl_space_drop_dims(space, isl_dim_in, 0, 1);
 191         space = isl_space_unprefix(space, isl_dim_in, "c_");
 192         space = isl_space_unprefix(space, isl_dim_out, "c_");
 193         nparam = isl_space_dim(space, isl_dim_in);
 194         space = isl_space_move_dims(space,
 195                                     isl_dim_param, 0, isl_dim_in, 0, nparam);
 196         space = isl_space_range(space);
 197
 198         return space;
 199 }
 200
 201 /* Return the rational universe basic set in the given space.
 202  */
 203 static __isl_give isl_basic_set *rational_universe(__isl_take isl_space *space)
 204 {
 205         isl_basic_set *bset;
 206
 207         bset = isl_basic_set_universe(space);
 208         bset = isl_basic_set_set_rational(bset);
 209
 210         return bset;
 211 }
 212
 213 /* Compute the dual of "bset" by applying Farkas' lemma.
 214  * As explained above, we add an extra dimension to represent
 215  * the coefficient of the constant term when going from solutions
 216  * to coefficients (shift == 1) and we drop the extra dimension when going
 217  * in the opposite direction (shift == -1).  "dim" is the space in which
 218  * the dual should be created.
 219  *
 220  * If "bset" is (obviously) empty, then the way this emptiness
 221  * is represented by the constraints does not allow for the application
 222  * of the standard farkas algorithm.  We therefore handle this case
 223  * specifically and return the universe basic set.
 224  */
 225 static __isl_give isl_basic_set *farkas(__isl_take isl_space *space,
 226         __isl_take isl_basic_set *bset, int shift)
 227 {
 228         int i, j, k;
 229         isl_basic_set *dual = NULL;
 230         unsigned total;
 231
 232         if (isl_basic_set_plain_is_empty(bset)) {
 233                 isl_basic_set_free(bset);
 234                 return rational_universe(space);
 235         }
 236
 237         total = isl_basic_set_total_dim(bset);
 238
 239         dual = isl_basic_set_alloc_space(space, bset->n_eq + bset->n_ineq,
 240                                         total, bset->n_ineq + (shift > 0));
 241         dual = isl_basic_set_set_rational(dual);
 242
 243         for (i = 0; i < bset->n_eq + bset->n_ineq; ++i) {
 244                 k = isl_basic_set_alloc_div(dual);
 245                 if (k < 0)
 246                         goto error;
 247                 isl_int_set_si(dual->div[k][0], 0);
 248         }
 249
 250         for (i = 0; i < total; ++i) {
 251                 k = isl_basic_set_alloc_equality(dual);
 252                 if (k < 0)
 253                         goto error;
 254                 isl_seq_clr(dual->eq[k], 1 + shift + total);
 255                 isl_int_set_si(dual->eq[k][1 + shift + i], -1);
 256                 for (j = 0; j < bset->n_eq; ++j)
 257                         isl_int_set(dual->eq[k][1 + shift + total + j],
 258                                     bset->eq[j][1 + i]);
 259                 for (j = 0; j < bset->n_ineq; ++j)
 260                         isl_int_set(dual->eq[k][1 + shift + total + bset->n_eq + j],
 261                                     bset->ineq[j][1 + i]);
 262         }
 263
 264         for (i = 0; i < bset->n_ineq; ++i) {
 265                 k = isl_basic_set_alloc_inequality(dual);
 266                 if (k < 0)
 267                         goto error;
 268                 isl_seq_clr(dual->ineq[k],
 269                             1 + shift + total + bset->n_eq + bset->n_ineq);
 270                 isl_int_set_si(dual->ineq[k][1 + shift + total + bset->n_eq + i], 1);
 271         }
 272
 273         if (shift > 0) {
 274                 k = isl_basic_set_alloc_inequality(dual);
 275                 if (k < 0)
 276                         goto error;
 277                 isl_seq_clr(dual->ineq[k], 2 + total);
 278                 isl_int_set_si(dual->ineq[k][1], 1);
 279                 for (j = 0; j < bset->n_eq; ++j)
 280                         isl_int_neg(dual->ineq[k][2 + total + j],
 281                                     bset->eq[j][0]);
 282                 for (j = 0; j < bset->n_ineq; ++j)
 283                         isl_int_neg(dual->ineq[k][2 + total + bset->n_eq + j],
 284                                     bset->ineq[j][0]);
 285         }
 286
 287         dual = isl_basic_set_remove_divs(dual);
 288         dual = isl_basic_set_simplify(dual);
 289         dual = isl_basic_set_finalize(dual);
 290
 291         isl_basic_set_free(bset);
 292         return dual;
 293 error:
 294         isl_basic_set_free(bset);
 295         isl_basic_set_free(dual);
 296         return NULL;
 297 }
 298
 299 /* Construct a basic set containing the tuples of coefficients of all
 300  * valid affine constraints on the given basic set.
 301  */
 302 __isl_give isl_basic_set *isl_basic_set_coefficients(
 303         __isl_take isl_basic_set *bset)
 304 {
 305         isl_space *dim;
 306
 307         if (!bset)
 308                 return NULL;
 309         if (bset->n_div)
 310                 isl_die(bset->ctx, isl_error_invalid,
 311                         "input set not allowed to have local variables",
 312                         goto error);
 313
 314         dim = isl_basic_set_get_space(bset);
 315         dim = isl_space_coefficients(dim);
 316
 317         return farkas(dim, bset, 1);
 318 error:
 319         isl_basic_set_free(bset);
 320         return NULL;
 321 }
 322
 323 /* Construct a basic set containing the elements that satisfy all
 324  * affine constraints whose coefficient tuples are
 325  * contained in the given basic set.
 326  */
 327 __isl_give isl_basic_set *isl_basic_set_solutions(
 328         __isl_take isl_basic_set *bset)
 329 {
 330         isl_space *dim;
 331
 332         if (!bset)
 333                 return NULL;
 334         if (bset->n_div)
 335                 isl_die(bset->ctx, isl_error_invalid,
 336                         "input set not allowed to have local variables",
 337                         goto error);
 338
 339         dim = isl_basic_set_get_space(bset);
 340         dim = isl_space_solutions(dim);
 341
 342         return farkas(dim, bset, -1);
 343 error:
 344         isl_basic_set_free(bset);
 345         return NULL;
 346 }
 347
 348 /* Construct a basic set containing the tuples of coefficients of all
 349  * valid affine constraints on the given set.
 350  */
 351 __isl_give isl_basic_set *isl_set_coefficients(__isl_take isl_set *set)
 352 {
 353         int i;
 354         isl_basic_set *coeff;
 355
 356         if (!set)
 357                 return NULL;
 358         if (set->n == 0) {
 359                 isl_space *space = isl_set_get_space(set);
 360                 space = isl_space_coefficients(space);
 361                 isl_set_free(set);
 362                 return rational_universe(space);
 363         }
 364
 365         coeff = isl_basic_set_coefficients(isl_basic_set_copy(set->p[0]));
 366
 367         for (i = 1; i < set->n; ++i) {
 368                 isl_basic_set *bset, *coeff_i;
 369                 bset = isl_basic_set_copy(set->p[i]);
 370                 coeff_i = isl_basic_set_coefficients(bset);
 371                 coeff = isl_basic_set_intersect(coeff, coeff_i);
 372         }
 373
 374         isl_set_free(set);
 375         return coeff;
 376 }
 377
 378 /* Wrapper around isl_basic_set_coefficients for use
 379  * as a isl_basic_set_list_map callback.
 380  */
 381 static __isl_give isl_basic_set *coefficients_wrap(
 382         __isl_take isl_basic_set *bset, void *user)
 383 {
 384         return isl_basic_set_coefficients(bset);
 385 }
 386
 387 /* Replace the elements of "list" by the result of applying
 388  * isl_basic_set_coefficients to them.
 389  */
 390 __isl_give isl_basic_set_list *isl_basic_set_list_coefficients(
 391         __isl_take isl_basic_set_list *list)
 392 {
 393         return isl_basic_set_list_map(list, &coefficients_wrap, NULL);
 394 }
 395
 396 /* Construct a basic set containing the elements that satisfy all
 397  * affine constraints whose coefficient tuples are
 398  * contained in the given set.
 399  */
 400 __isl_give isl_basic_set *isl_set_solutions(__isl_take isl_set *set)
 401 {
 402         int i;
 403         isl_basic_set *sol;
 404
 405         if (!set)
 406                 return NULL;
 407         if (set->n == 0) {
 408                 isl_space *space = isl_set_get_space(set);
 409                 space = isl_space_solutions(space);
 410                 isl_set_free(set);
 411                 return rational_universe(space);
 412         }
 413
 414         sol = isl_basic_set_solutions(isl_basic_set_copy(set->p[0]));
 415
 416         for (i = 1; i < set->n; ++i) {
 417                 isl_basic_set *bset, *sol_i;
 418                 bset = isl_basic_set_copy(set->p[i]);
 419                 sol_i = isl_basic_set_solutions(bset);
 420                 sol = isl_basic_set_intersect(sol, sol_i);
 421         }
 422
 423         isl_set_free(set);
 424         return sol;
 425 }