isl_morph.c

   1 /*
   2  * Copyright 2010-2011 INRIA Saclay
   3  * Copyright 2014      Ecole Normale Superieure
   4  *
   5  * Use of this software is governed by the MIT license
   6  *
   7  * Written by Sven Verdoolaege, INRIA Saclay - Ile-de-France,
   8  * Parc Club Orsay Universite, ZAC des vignes, 4 rue Jacques Monod,
   9  * 91893 Orsay, France
  10  * and Ecole Normale Superieure, 45 rue d'Ulm, 75230 Paris, France
  11  */
  12
  13 #include <isl_map_private.h>
  14 #include <isl_aff_private.h>
  15 #include <isl_morph.h>
  16 #include <isl_seq.h>
  17 #include <isl_mat_private.h>
  18 #include <isl_space_private.h>
  19 #include <isl_equalities.h>
  20
  21 isl_ctx *isl_morph_get_ctx(__isl_keep isl_morph *morph)
  22 {
  23         if (!morph)
  24                 return NULL;
  25         return isl_basic_set_get_ctx(morph->dom);
  26 }
  27
  28 __isl_give isl_morph *isl_morph_alloc(
  29         __isl_take isl_basic_set *dom, __isl_take isl_basic_set *ran,
  30         __isl_take isl_mat *map, __isl_take isl_mat *inv)
  31 {
  32         isl_morph *morph;
  33
  34         if (!dom || !ran || !map || !inv)
  35                 goto error;
  36
  37         morph = isl_alloc_type(dom->ctx, struct isl_morph);
  38         if (!morph)
  39                 goto error;
  40
  41         morph->ref = 1;
  42         morph->dom = dom;
  43         morph->ran = ran;
  44         morph->map = map;
  45         morph->inv = inv;
  46
  47         return morph;
  48 error:
  49         isl_basic_set_free(dom);
  50         isl_basic_set_free(ran);
  51         isl_mat_free(map);
  52         isl_mat_free(inv);
  53         return NULL;
  54 }
  55
  56 __isl_give isl_morph *isl_morph_copy(__isl_keep isl_morph *morph)
  57 {
  58         if (!morph)
  59                 return NULL;
  60
  61         morph->ref++;
  62         return morph;
  63 }
  64
  65 __isl_give isl_morph *isl_morph_dup(__isl_keep isl_morph *morph)
  66 {
  67         if (!morph)
  68                 return NULL;
  69
  70         return isl_morph_alloc(isl_basic_set_copy(morph->dom),
  71                 isl_basic_set_copy(morph->ran),
  72                 isl_mat_copy(morph->map), isl_mat_copy(morph->inv));
  73 }
  74
  75 __isl_give isl_morph *isl_morph_cow(__isl_take isl_morph *morph)
  76 {
  77         if (!morph)
  78                 return NULL;
  79
  80         if (morph->ref == 1)
  81                 return morph;
  82         morph->ref--;
  83         return isl_morph_dup(morph);
  84 }
  85
  86 void isl_morph_free(__isl_take isl_morph *morph)
  87 {
  88         if (!morph)
  89                 return;
  90
  91         if (--morph->ref > 0)
  92                 return;
  93
  94         isl_basic_set_free(morph->dom);
  95         isl_basic_set_free(morph->ran);
  96         isl_mat_free(morph->map);
  97         isl_mat_free(morph->inv);
  98         free(morph);
  99 }
 100
 101 /* Is "morph" an identity on the parameters?
 102  */
 103 static int identity_on_parameters(__isl_keep isl_morph *morph)
 104 {
 105         int is_identity;
 106         unsigned nparam;
 107         isl_mat *sub;
 108
 109         nparam = isl_morph_dom_dim(morph, isl_dim_param);
 110         if (nparam != isl_morph_ran_dim(morph, isl_dim_param))
 111                 return 0;
 112         if (nparam == 0)
 113                 return 1;
 114         sub = isl_mat_sub_alloc(morph->map, 0, 1 + nparam, 0, 1 + nparam);
 115         is_identity = isl_mat_is_scaled_identity(sub);
 116         isl_mat_free(sub);
 117
 118         return is_identity;
 119 }
 120
 121 /* Return an affine expression of the variables of the range of "morph"
 122  * in terms of the parameters and the variables of the domain on "morph".
 123  *
 124  * In order for the space manipulations to make sense, we require
 125  * that the parameters are not modified by "morph".
 126  */
 127 __isl_give isl_multi_aff *isl_morph_get_var_multi_aff(
 128         __isl_keep isl_morph *morph)
 129 {
 130         isl_space *dom, *ran, *space;
 131         isl_local_space *ls;
 132         isl_multi_aff *ma;
 133         unsigned nparam, nvar;
 134         int i;
 135         int is_identity;
 136
 137         if (!morph)
 138                 return NULL;
 139
 140         is_identity = identity_on_parameters(morph);
 141         if (is_identity < 0)
 142                 return NULL;
 143         if (!is_identity)
 144                 isl_die(isl_morph_get_ctx(morph), isl_error_invalid,
 145                         "cannot handle parameter compression", return NULL);
 146
 147         dom = isl_morph_get_dom_space(morph);
 148         ls = isl_local_space_from_space(isl_space_copy(dom));
 149         ran = isl_morph_get_ran_space(morph);
 150         space = isl_space_map_from_domain_and_range(dom, ran);
 151         ma = isl_multi_aff_zero(space);
 152
 153         nparam = isl_multi_aff_dim(ma, isl_dim_param);
 154         nvar = isl_multi_aff_dim(ma, isl_dim_out);
 155         for (i = 0; i < nvar; ++i) {
 156                 isl_val *val;
 157                 isl_vec *v;
 158                 isl_aff *aff;
 159
 160                 v = isl_mat_get_row(morph->map, 1 + nparam + i);
 161                 v = isl_vec_insert_els(v, 0, 1);
 162                 val = isl_mat_get_element_val(morph->map, 0, 0);
 163                 v = isl_vec_set_element_val(v, 0, val);
 164                 aff = isl_aff_alloc_vec(isl_local_space_copy(ls), v);
 165                 ma = isl_multi_aff_set_aff(ma, i, aff);
 166         }
 167
 168         isl_local_space_free(ls);
 169         return ma;
 170 }
 171
 172 /* Return the domain space of "morph".
 173  */
 174 __isl_give isl_space *isl_morph_get_dom_space(__isl_keep isl_morph *morph)
 175 {
 176         if (!morph)
 177                 return NULL;
 178
 179         return isl_basic_set_get_space(morph->dom);
 180 }
 181
 182 __isl_give isl_space *isl_morph_get_ran_space(__isl_keep isl_morph *morph)
 183 {
 184         if (!morph)
 185                 return NULL;
 186
 187         return isl_space_copy(morph->ran->dim);
 188 }
 189
 190 unsigned isl_morph_dom_dim(__isl_keep isl_morph *morph, enum isl_dim_type type)
 191 {
 192         if (!morph)
 193                 return 0;
 194
 195         return isl_basic_set_dim(morph->dom, type);
 196 }
 197
 198 unsigned isl_morph_ran_dim(__isl_keep isl_morph *morph, enum isl_dim_type type)
 199 {
 200         if (!morph)
 201                 return 0;
 202
 203         return isl_basic_set_dim(morph->ran, type);
 204 }
 205
 206 __isl_give isl_morph *isl_morph_remove_dom_dims(__isl_take isl_morph *morph,
 207         enum isl_dim_type type, unsigned first, unsigned n)
 208 {
 209         unsigned dom_offset;
 210
 211         if (n == 0)
 212                 return morph;
 213
 214         morph = isl_morph_cow(morph);
 215         if (!morph)
 216                 return NULL;
 217
 218         dom_offset = 1 + isl_space_offset(morph->dom->dim, type);
 219
 220         morph->dom = isl_basic_set_remove_dims(morph->dom, type, first, n);
 221
 222         morph->map = isl_mat_drop_cols(morph->map, dom_offset + first, n);
 223
 224         morph->inv = isl_mat_drop_rows(morph->inv, dom_offset + first, n);
 225
 226         if (morph->dom && morph->ran && morph->map && morph->inv)
 227                 return morph;
 228
 229         isl_morph_free(morph);
 230         return NULL;
 231 }
 232
 233 __isl_give isl_morph *isl_morph_remove_ran_dims(__isl_take isl_morph *morph,
 234         enum isl_dim_type type, unsigned first, unsigned n)
 235 {
 236         unsigned ran_offset;
 237
 238         if (n == 0)
 239                 return morph;
 240
 241         morph = isl_morph_cow(morph);
 242         if (!morph)
 243                 return NULL;
 244
 245         ran_offset = 1 + isl_space_offset(morph->ran->dim, type);
 246
 247         morph->ran = isl_basic_set_remove_dims(morph->ran, type, first, n);
 248
 249         morph->map = isl_mat_drop_rows(morph->map, ran_offset + first, n);
 250
 251         morph->inv = isl_mat_drop_cols(morph->inv, ran_offset + first, n);
 252
 253         if (morph->dom && morph->ran && morph->map && morph->inv)
 254                 return morph;
 255
 256         isl_morph_free(morph);
 257         return NULL;
 258 }
 259
 260 /* Project domain of morph onto its parameter domain.
 261  */
 262 __isl_give isl_morph *isl_morph_dom_params(__isl_take isl_morph *morph)
 263 {
 264         unsigned n;
 265
 266         morph = isl_morph_cow(morph);
 267         if (!morph)
 268                 return NULL;
 269         n = isl_basic_set_dim(morph->dom, isl_dim_set);
 270         morph = isl_morph_remove_dom_dims(morph, isl_dim_set, 0, n);
 271         if (!morph)
 272                 return NULL;
 273         morph->dom = isl_basic_set_params(morph->dom);
 274         if (morph->dom)
 275                 return morph;
 276
 277         isl_morph_free(morph);
 278         return NULL;
 279 }
 280
 281 /* Project range of morph onto its parameter domain.
 282  */
 283 __isl_give isl_morph *isl_morph_ran_params(__isl_take isl_morph *morph)
 284 {
 285         unsigned n;
 286
 287         morph = isl_morph_cow(morph);
 288         if (!morph)
 289                 return NULL;
 290         n = isl_basic_set_dim(morph->ran, isl_dim_set);
 291         morph = isl_morph_remove_ran_dims(morph, isl_dim_set, 0, n);
 292         if (!morph)
 293                 return NULL;
 294         morph->ran = isl_basic_set_params(morph->ran);
 295         if (morph->ran)
 296                 return morph;
 297
 298         isl_morph_free(morph);
 299         return NULL;
 300 }
 301
 302 void isl_morph_print_internal(__isl_take isl_morph *morph, FILE *out)
 303 {
 304         if (!morph)
 305                 return;
 306
 307         isl_basic_set_print(morph->dom, out, 0, "", "", ISL_FORMAT_ISL);
 308         isl_basic_set_print(morph->ran, out, 0, "", "", ISL_FORMAT_ISL);
 309         isl_mat_print_internal(morph->map, out, 4);
 310         isl_mat_print_internal(morph->inv, out, 4);
 311 }
 312
 313 void isl_morph_dump(__isl_take isl_morph *morph)
 314 {
 315         isl_morph_print_internal(morph, stderr);
 316 }
 317
 318 __isl_give isl_morph *isl_morph_identity(__isl_keep isl_basic_set *bset)
 319 {
 320         isl_mat *id;
 321         isl_basic_set *universe;
 322         unsigned total;
 323
 324         if (!bset)
 325                 return NULL;
 326
 327         total = isl_basic_set_total_dim(bset);
 328         id = isl_mat_identity(bset->ctx, 1 + total);
 329         universe = isl_basic_set_universe(isl_space_copy(bset->dim));
 330
 331         return isl_morph_alloc(universe, isl_basic_set_copy(universe),
 332                 id, isl_mat_copy(id));
 333 }
 334
 335 /* Create a(n identity) morphism between empty sets of the same dimension
 336  * a "bset".
 337  */
 338 __isl_give isl_morph *isl_morph_empty(__isl_keep isl_basic_set *bset)
 339 {
 340         isl_mat *id;
 341         isl_basic_set *empty;
 342         unsigned total;
 343
 344         if (!bset)
 345                 return NULL;
 346
 347         total = isl_basic_set_total_dim(bset);
 348         id = isl_mat_identity(bset->ctx, 1 + total);
 349         empty = isl_basic_set_empty(isl_space_copy(bset->dim));
 350
 351         return isl_morph_alloc(empty, isl_basic_set_copy(empty),
 352                 id, isl_mat_copy(id));
 353 }
 354
 355 /* Given a matrix that maps a (possibly) parametric domain to
 356  * a parametric domain, add in rows that map the "nparam" parameters onto
 357  * themselves.
 358  */
 359 static __isl_give isl_mat *insert_parameter_rows(__isl_take isl_mat *mat,
 360         unsigned nparam)
 361 {
 362         int i;
 363
 364         if (nparam == 0)
 365                 return mat;
 366         if (!mat)
 367                 return NULL;
 368
 369         mat = isl_mat_insert_rows(mat, 1, nparam);
 370         if (!mat)
 371                 return NULL;
 372
 373         for (i = 0; i < nparam; ++i) {
 374                 isl_seq_clr(mat->row[1 + i], mat->n_col);
 375                 isl_int_set(mat->row[1 + i][1 + i], mat->row[0][0]);
 376         }
 377
 378         return mat;
 379 }
 380
 381 /* Construct a basic set described by the "n" equalities of "bset" starting
 382  * at "first".
 383  */
 384 static __isl_give isl_basic_set *copy_equalities(__isl_keep isl_basic_set *bset,
 385         unsigned first, unsigned n)
 386 {
 387         int i, k;
 388         isl_basic_set *eq;
 389         unsigned total;
 390
 391         isl_assert(bset->ctx, bset->n_div == 0, return NULL);
 392
 393         total = isl_basic_set_total_dim(bset);
 394         eq = isl_basic_set_alloc_space(isl_space_copy(bset->dim), 0, n, 0);
 395         if (!eq)
 396                 return NULL;
 397         for (i = 0; i < n; ++i) {
 398                 k = isl_basic_set_alloc_equality(eq);
 399                 if (k < 0)
 400                         goto error;
 401                 isl_seq_cpy(eq->eq[k], bset->eq[first + k], 1 + total);
 402         }
 403
 404         return eq;
 405 error:
 406         isl_basic_set_free(eq);
 407         return NULL;
 408 }
 409
 410 /* Given a basic set, exploit the equalties in the basic set to construct
 411  * a morphishm that maps the basic set to a lower-dimensional space.
 412  * Specifically, the morphism reduces the number of dimensions of type "type".
 413  *
 414  * This function is a slight generalization of isl_mat_variable_compression
 415  * in that it allows the input to be parametric and that it allows for the
 416  * compression of either parameters or set variables.
 417  *
 418  * We first select the equalities of interest, that is those that involve
 419  * variables of type "type" and no later variables.
 420  * Denote those equalities as
 421  *
 422  *              -C(p) + M x = 0
 423  *
 424  * where C(p) depends on the parameters if type == isl_dim_set and
 425  * is a constant if type == isl_dim_param.
 426  *
 427  * First compute the (left) Hermite normal form of M,
 428  *
 429  *              M [U1 U2] = M U = H = [H1 0]
 430  * or
 431  *                            M = H Q = [H1 0] [Q1]
 432  *                                             [Q2]
 433  *
 434  * with U, Q unimodular, Q = U^{-1} (and H lower triangular).
 435  * Define the transformed variables as
 436  *
 437  *              x = [U1 U2] [ x1' ] = [U1 U2] [Q1] x
 438  *                          [ x2' ]           [Q2]
 439  *
 440  * The equalities then become
 441  *
 442  *              -C(p) + H1 x1' = 0   or   x1' = H1^{-1} C(p) = C'(p)
 443  *
 444  * If the denominator of the constant term does not divide the
 445  * the common denominator of the parametric terms, then every
 446  * integer point is mapped to a non-integer point and then the original set has no
 447  * integer solutions (since the x' are a unimodular transformation
 448  * of the x).  In this case, an empty morphism is returned.
 449  * Otherwise, the transformation is given by
 450  *
 451  *              x = U1 H1^{-1} C(p) + U2 x2'
 452  *
 453  * The inverse transformation is simply
 454  *
 455  *              x2' = Q2 x
 456  *
 457  * Both matrices are extended to map the full original space to the full
 458  * compressed space.
 459  */
 460 __isl_give isl_morph *isl_basic_set_variable_compression(
 461         __isl_keep isl_basic_set *bset, enum isl_dim_type type)
 462 {
 463         unsigned otype;
 464         unsigned ntype;
 465         unsigned orest;
 466         unsigned nrest;
 467         int f_eq, n_eq;
 468         isl_space *dim;
 469         isl_mat *H, *U, *Q, *C = NULL, *H1, *U1, *U2;
 470         isl_basic_set *dom, *ran;
 471
 472         if (!bset)
 473                 return NULL;
 474
 475         if (isl_basic_set_plain_is_empty(bset))
 476                 return isl_morph_empty(bset);
 477
 478         isl_assert(bset->ctx, bset->n_div == 0, return NULL);
 479
 480         otype = 1 + isl_space_offset(bset->dim, type);
 481         ntype = isl_basic_set_dim(bset, type);
 482         orest = otype + ntype;
 483         nrest = isl_basic_set_total_dim(bset) - (orest - 1);
 484
 485         for (f_eq = 0; f_eq < bset->n_eq; ++f_eq)
 486                 if (isl_seq_first_non_zero(bset->eq[f_eq] + orest, nrest) == -1)
 487                         break;
 488         for (n_eq = 0; f_eq + n_eq < bset->n_eq; ++n_eq)
 489                 if (isl_seq_first_non_zero(bset->eq[f_eq + n_eq] + otype, ntype) == -1)
 490                         break;
 491         if (n_eq == 0)
 492                 return isl_morph_identity(bset);
 493
 494         H = isl_mat_sub_alloc6(bset->ctx, bset->eq, f_eq, n_eq, otype, ntype);
 495         H = isl_mat_left_hermite(H, 0, &U, &Q);
 496         if (!H || !U || !Q)
 497                 goto error;
 498         Q = isl_mat_drop_rows(Q, 0, n_eq);
 499         Q = isl_mat_diagonal(isl_mat_identity(bset->ctx, otype), Q);
 500         Q = isl_mat_diagonal(Q, isl_mat_identity(bset->ctx, nrest));
 501         C = isl_mat_alloc(bset->ctx, 1 + n_eq, otype);
 502         if (!C)
 503                 goto error;
 504         isl_int_set_si(C->row[0][0], 1);
 505         isl_seq_clr(C->row[0] + 1, otype - 1);
 506         isl_mat_sub_neg(C->ctx, C->row + 1, bset->eq + f_eq, n_eq, 0, 0, otype);
 507         H1 = isl_mat_sub_alloc(H, 0, H->n_row, 0, H->n_row);
 508         H1 = isl_mat_lin_to_aff(H1);
 509         C = isl_mat_inverse_product(H1, C);
 510         if (!C)
 511                 goto error;
 512         isl_mat_free(H);
 513
 514         if (!isl_int_is_one(C->row[0][0])) {
 515                 int i;
 516                 isl_int g;
 517
 518                 isl_int_init(g);
 519                 for (i = 0; i < n_eq; ++i) {
 520                         isl_seq_gcd(C->row[1 + i] + 1, otype - 1, &g);
 521                         isl_int_gcd(g, g, C->row[0][0]);
 522                         if (!isl_int_is_divisible_by(C->row[1 + i][0], g))
 523                                 break;
 524                 }
 525                 isl_int_clear(g);
 526
 527                 if (i < n_eq) {
 528                         isl_mat_free(C);
 529                         isl_mat_free(U);
 530                         isl_mat_free(Q);
 531                         return isl_morph_empty(bset);
 532                 }
 533
 534                 C = isl_mat_normalize(C);
 535         }
 536
 537         U1 = isl_mat_sub_alloc(U, 0, U->n_row, 0, n_eq);
 538         U1 = isl_mat_lin_to_aff(U1);
 539         U2 = isl_mat_sub_alloc(U, 0, U->n_row, n_eq, U->n_row - n_eq);
 540         U2 = isl_mat_lin_to_aff(U2);
 541         isl_mat_free(U);
 542
 543         C = isl_mat_product(U1, C);
 544         C = isl_mat_aff_direct_sum(C, U2);
 545         C = insert_parameter_rows(C, otype - 1);
 546         C = isl_mat_diagonal(C, isl_mat_identity(bset->ctx, nrest));
 547
 548         dim = isl_space_copy(bset->dim);
 549         dim = isl_space_drop_dims(dim, type, 0, ntype);
 550         dim = isl_space_add_dims(dim, type, ntype - n_eq);
 551         ran = isl_basic_set_universe(dim);
 552         dom = copy_equalities(bset, f_eq, n_eq);
 553
 554         return isl_morph_alloc(dom, ran, Q, C);
 555 error:
 556         isl_mat_free(C);
 557         isl_mat_free(H);
 558         isl_mat_free(U);
 559         isl_mat_free(Q);
 560         return NULL;
 561 }
 562
 563 /* Construct a parameter compression for "bset".
 564  * We basically just call isl_mat_parameter_compression with the right input
 565  * and then extend the resulting matrix to include the variables.
 566  *
 567  * The implementation assumes that "bset" does not have any equalities
 568  * that only involve the parameters and that isl_basic_set_gauss has
 569  * been applied to "bset".
 570  *
 571  * Let the equalities be given as
 572  *
 573  *      B(p) + A x = 0.
 574  *
 575  * We use isl_mat_parameter_compression_ext to compute the compression
 576  *
 577  *      p = T p'.
 578  */
 579 __isl_give isl_morph *isl_basic_set_parameter_compression(
 580         __isl_keep isl_basic_set *bset)
 581 {
 582         unsigned nparam;
 583         unsigned nvar;
 584         unsigned n_div;
 585         int n_eq;
 586         isl_mat *H, *B;
 587         isl_mat *map, *inv;
 588         isl_basic_set *dom, *ran;
 589
 590         if (!bset)
 591                 return NULL;
 592
 593         if (isl_basic_set_plain_is_empty(bset))
 594                 return isl_morph_empty(bset);
 595         if (bset->n_eq == 0)
 596                 return isl_morph_identity(bset);
 597
 598         n_eq = bset->n_eq;
 599         nparam = isl_basic_set_dim(bset, isl_dim_param);
 600         nvar = isl_basic_set_dim(bset, isl_dim_set);
 601         n_div = isl_basic_set_dim(bset, isl_dim_div);
 602
 603         if (isl_seq_first_non_zero(bset->eq[bset->n_eq - 1] + 1 + nparam,
 604                                     nvar + n_div) == -1)
 605                 isl_die(isl_basic_set_get_ctx(bset), isl_error_invalid,
 606                         "input not allowed to have parameter equalities",
 607                         return NULL);
 608         if (n_eq > nvar + n_div)
 609                 isl_die(isl_basic_set_get_ctx(bset), isl_error_invalid,
 610                         "input not gaussed", return NULL);
 611
 612         B = isl_mat_sub_alloc6(bset->ctx, bset->eq, 0, n_eq, 0, 1 + nparam);
 613         H = isl_mat_sub_alloc6(bset->ctx, bset->eq,
 614                                 0, n_eq, 1 + nparam, nvar + n_div);
 615         inv = isl_mat_parameter_compression_ext(B, H);
 616         inv = isl_mat_diagonal(inv, isl_mat_identity(bset->ctx, nvar));
 617         map = isl_mat_right_inverse(isl_mat_copy(inv));
 618
 619         dom = isl_basic_set_universe(isl_space_copy(bset->dim));
 620         ran = isl_basic_set_universe(isl_space_copy(bset->dim));
 621
 622         return isl_morph_alloc(dom, ran, map, inv);
 623 }
 624
 625 /* Add stride constraints to "bset" based on the inverse mapping
 626  * that was plugged in.  In particular, if morph maps x' to x,
 627  * the the constraints of the original input
 628  *
 629  *      A x' + b >= 0
 630  *
 631  * have been rewritten to
 632  *
 633  *      A inv x + b >= 0
 634  *
 635  * However, this substitution may loose information on the integrality of x',
 636  * so we need to impose that
 637  *
 638  *      inv x
 639  *
 640  * is integral.  If inv = B/d, this means that we need to impose that
 641  *
 642  *      B x = 0         mod d
 643  *
 644  * or
 645  *
 646  *      exists alpha in Z^m: B x = d alpha
 647  *
 648  * This function is similar to add_strides in isl_affine_hull.c
 649  */
 650 static __isl_give isl_basic_set *add_strides(__isl_take isl_basic_set *bset,
 651         __isl_keep isl_morph *morph)
 652 {
 653         int i, div, k;
 654         isl_int gcd;
 655
 656         if (isl_int_is_one(morph->inv->row[0][0]))
 657                 return bset;
 658
 659         isl_int_init(gcd);
 660
 661         for (i = 0; 1 + i < morph->inv->n_row; ++i) {
 662                 isl_seq_gcd(morph->inv->row[1 + i], morph->inv->n_col, &gcd);
 663                 if (isl_int_is_divisible_by(gcd, morph->inv->row[0][0]))
 664                         continue;
 665                 div = isl_basic_set_alloc_div(bset);
 666                 if (div < 0)
 667                         goto error;
 668                 isl_int_set_si(bset->div[div][0], 0);
 669                 k = isl_basic_set_alloc_equality(bset);
 670                 if (k < 0)
 671                         goto error;
 672                 isl_seq_cpy(bset->eq[k], morph->inv->row[1 + i],
 673                             morph->inv->n_col);
 674                 isl_seq_clr(bset->eq[k] + morph->inv->n_col, bset->n_div);
 675                 isl_int_set(bset->eq[k][morph->inv->n_col + div],
 676                             morph->inv->row[0][0]);
 677         }
 678
 679         isl_int_clear(gcd);
 680
 681         return bset;
 682 error:
 683         isl_int_clear(gcd);
 684         isl_basic_set_free(bset);
 685         return NULL;
 686 }
 687
 688 /* Apply the morphism to the basic set.
 689  * We basically just compute the preimage of "bset" under the inverse mapping
 690  * in morph, add in stride constraints and intersect with the range
 691  * of the morphism.
 692  */
 693 __isl_give isl_basic_set *isl_morph_basic_set(__isl_take isl_morph *morph,
 694         __isl_take isl_basic_set *bset)
 695 {
 696         isl_basic_set *res = NULL;
 697         isl_mat *mat = NULL;
 698         int i, k;
 699         int max_stride;
 700
 701         if (!morph || !bset)
 702                 goto error;
 703
 704         isl_assert(bset->ctx, isl_space_is_equal(bset->dim, morph->dom->dim),
 705                     goto error);
 706
 707         max_stride = morph->inv->n_row - 1;
 708         if (isl_int_is_one(morph->inv->row[0][0]))
 709                 max_stride = 0;
 710         res = isl_basic_set_alloc_space(isl_space_copy(morph->ran->dim),
 711                 bset->n_div + max_stride, bset->n_eq + max_stride, bset->n_ineq);
 712
 713         for (i = 0; i < bset->n_div; ++i)
 714                 if (isl_basic_set_alloc_div(res) < 0)
 715                         goto error;
 716
 717         mat = isl_mat_sub_alloc6(bset->ctx, bset->eq, 0, bset->n_eq,
 718                                         0, morph->inv->n_row);
 719         mat = isl_mat_product(mat, isl_mat_copy(morph->inv));
 720         if (!mat)
 721                 goto error;
 722         for (i = 0; i < bset->n_eq; ++i) {
 723                 k = isl_basic_set_alloc_equality(res);
 724                 if (k < 0)
 725                         goto error;
 726                 isl_seq_cpy(res->eq[k], mat->row[i], mat->n_col);
 727                 isl_seq_scale(res->eq[k] + mat->n_col, bset->eq[i] + mat->n_col,
 728                                 morph->inv->row[0][0], bset->n_div);
 729         }
 730         isl_mat_free(mat);
 731
 732         mat = isl_mat_sub_alloc6(bset->ctx, bset->ineq, 0, bset->n_ineq,
 733                                         0, morph->inv->n_row);
 734         mat = isl_mat_product(mat, isl_mat_copy(morph->inv));
 735         if (!mat)
 736                 goto error;
 737         for (i = 0; i < bset->n_ineq; ++i) {
 738                 k = isl_basic_set_alloc_inequality(res);
 739                 if (k < 0)
 740                         goto error;
 741                 isl_seq_cpy(res->ineq[k], mat->row[i], mat->n_col);
 742                 isl_seq_scale(res->ineq[k] + mat->n_col,
 743                                 bset->ineq[i] + mat->n_col,
 744                                 morph->inv->row[0][0], bset->n_div);
 745         }
 746         isl_mat_free(mat);
 747
 748         mat = isl_mat_sub_alloc6(bset->ctx, bset->div, 0, bset->n_div,
 749                                         1, morph->inv->n_row);
 750         mat = isl_mat_product(mat, isl_mat_copy(morph->inv));
 751         if (!mat)
 752                 goto error;
 753         for (i = 0; i < bset->n_div; ++i) {
 754                 isl_int_mul(res->div[i][0],
 755                                 morph->inv->row[0][0], bset->div[i][0]);
 756                 isl_seq_cpy(res->div[i] + 1, mat->row[i], mat->n_col);
 757                 isl_seq_scale(res->div[i] + 1 + mat->n_col,
 758                                 bset->div[i] + 1 + mat->n_col,
 759                                 morph->inv->row[0][0], bset->n_div);
 760         }
 761         isl_mat_free(mat);
 762
 763         res = add_strides(res, morph);
 764
 765         if (isl_basic_set_is_rational(bset))
 766                 res = isl_basic_set_set_rational(res);
 767
 768         res = isl_basic_set_simplify(res);
 769         res = isl_basic_set_finalize(res);
 770
 771         res = isl_basic_set_intersect(res, isl_basic_set_copy(morph->ran));
 772
 773         isl_morph_free(morph);
 774         isl_basic_set_free(bset);
 775         return res;
 776 error:
 777         isl_mat_free(mat);
 778         isl_morph_free(morph);
 779         isl_basic_set_free(bset);
 780         isl_basic_set_free(res);
 781         return NULL;
 782 }
 783
 784 /* Apply the morphism to the set.
 785  */
 786 __isl_give isl_set *isl_morph_set(__isl_take isl_morph *morph,
 787         __isl_take isl_set *set)
 788 {
 789         int i;
 790
 791         if (!morph || !set)
 792                 goto error;
 793
 794         isl_assert(set->ctx, isl_space_is_equal(set->dim, morph->dom->dim), goto error);
 795
 796         set = isl_set_cow(set);
 797         if (!set)
 798                 goto error;
 799
 800         isl_space_free(set->dim);
 801         set->dim = isl_space_copy(morph->ran->dim);
 802         if (!set->dim)
 803                 goto error;
 804
 805         for (i = 0; i < set->n; ++i) {
 806                 set->p[i] = isl_morph_basic_set(isl_morph_copy(morph), set->p[i]);
 807                 if (!set->p[i])
 808                         goto error;
 809         }
 810
 811         isl_morph_free(morph);
 812
 813         ISL_F_CLR(set, ISL_SET_NORMALIZED);
 814
 815         return set;
 816 error:
 817         isl_set_free(set);
 818         isl_morph_free(morph);
 819         return NULL;
 820 }
 821
 822 /* Construct a morphism that first does morph2 and then morph1.
 823  */
 824 __isl_give isl_morph *isl_morph_compose(__isl_take isl_morph *morph1,
 825         __isl_take isl_morph *morph2)
 826 {
 827         isl_mat *map, *inv;
 828         isl_basic_set *dom, *ran;
 829
 830         if (!morph1 || !morph2)
 831                 goto error;
 832
 833         map = isl_mat_product(isl_mat_copy(morph1->map), isl_mat_copy(morph2->map));
 834         inv = isl_mat_product(isl_mat_copy(morph2->inv), isl_mat_copy(morph1->inv));
 835         dom = isl_morph_basic_set(isl_morph_inverse(isl_morph_copy(morph2)),
 836                                   isl_basic_set_copy(morph1->dom));
 837         dom = isl_basic_set_intersect(dom, isl_basic_set_copy(morph2->dom));
 838         ran = isl_morph_basic_set(isl_morph_copy(morph1),
 839                                   isl_basic_set_copy(morph2->ran));
 840         ran = isl_basic_set_intersect(ran, isl_basic_set_copy(morph1->ran));
 841
 842         isl_morph_free(morph1);
 843         isl_morph_free(morph2);
 844
 845         return isl_morph_alloc(dom, ran, map, inv);
 846 error:
 847         isl_morph_free(morph1);
 848         isl_morph_free(morph2);
 849         return NULL;
 850 }
 851
 852 __isl_give isl_morph *isl_morph_inverse(__isl_take isl_morph *morph)
 853 {
 854         isl_basic_set *bset;
 855         isl_mat *mat;
 856
 857         morph = isl_morph_cow(morph);
 858         if (!morph)
 859                 return NULL;
 860
 861         bset = morph->dom;
 862         morph->dom = morph->ran;
 863         morph->ran = bset;
 864
 865         mat = morph->map;
 866         morph->map = morph->inv;
 867         morph->inv = mat;
 868
 869         return morph;
 870 }
 871
 872 /* We detect all the equalities first to avoid implicit equalties
 873  * being discovered during the computations.  In particular,
 874  * the compression on the variables could expose additional stride
 875  * constraints on the parameters.  This would result in existentially
 876  * quantified variables after applying the resulting morph, which
 877  * in turn could break invariants of the calling functions.
 878  */
 879 __isl_give isl_morph *isl_basic_set_full_compression(
 880         __isl_keep isl_basic_set *bset)
 881 {
 882         isl_morph *morph, *morph2;
 883
 884         bset = isl_basic_set_copy(bset);
 885         bset = isl_basic_set_detect_equalities(bset);
 886
 887         morph = isl_basic_set_variable_compression(bset, isl_dim_param);
 888         bset = isl_morph_basic_set(isl_morph_copy(morph), bset);
 889
 890         morph2 = isl_basic_set_parameter_compression(bset);
 891         bset = isl_morph_basic_set(isl_morph_copy(morph2), bset);
 892
 893         morph = isl_morph_compose(morph2, morph);
 894
 895         morph2 = isl_basic_set_variable_compression(bset, isl_dim_set);
 896         isl_basic_set_free(bset);
 897
 898         morph = isl_morph_compose(morph2, morph);
 899
 900         return morph;
 901 }
 902
 903 __isl_give isl_vec *isl_morph_vec(__isl_take isl_morph *morph,
 904         __isl_take isl_vec *vec)
 905 {
 906         if (!morph)
 907                 goto error;
 908
 909         vec = isl_mat_vec_product(isl_mat_copy(morph->map), vec);
 910
 911         isl_morph_free(morph);
 912         return vec;
 913 error:
 914         isl_morph_free(morph);
 915         isl_vec_free(vec);
 916         return NULL;
 917 }