isl_aff.c

   1 /*
   2  * Copyright 2011      INRIA Saclay
   3  * Copyright 2011      Sven Verdoolaege
   4  * Copyright 2012      Ecole Normale Superieure
   5  *
   6  * Use of this software is governed by the GNU LGPLv2.1 license
   7  *
   8  * Written by Sven Verdoolaege, INRIA Saclay - Ile-de-France,
   9  * Parc Club Orsay Universite, ZAC des vignes, 4 rue Jacques Monod,
  10  * 91893 Orsay, France
  11  * and Ecole Normale Superieure, 45 rue d’Ulm, 75230 Paris, France
  12  */
  13
  14 #include <isl_ctx_private.h>
  15 #include <isl_map_private.h>
  16 #include <isl_aff_private.h>
  17 #include <isl_space_private.h>
  18 #include <isl_local_space_private.h>
  19 #include <isl_mat_private.h>
  20 #include <isl_list_private.h>
  21 #include <isl/constraint.h>
  22 #include <isl/seq.h>
  23 #include <isl/set.h>
  24 #include <isl_config.h>
  25
  26 __isl_give isl_aff *isl_aff_alloc_vec(__isl_take isl_local_space *ls,
  27         __isl_take isl_vec *v)
  28 {
  29         isl_aff *aff;
  30
  31         if (!ls || !v)
  32                 goto error;
  33
  34         aff = isl_calloc_type(v->ctx, struct isl_aff);
  35         if (!aff)
  36                 goto error;
  37
  38         aff->ref = 1;
  39         aff->ls = ls;
  40         aff->v = v;
  41
  42         return aff;
  43 error:
  44         isl_local_space_free(ls);
  45         isl_vec_free(v);
  46         return NULL;
  47 }
  48
  49 __isl_give isl_aff *isl_aff_alloc(__isl_take isl_local_space *ls)
  50 {
  51         isl_ctx *ctx;
  52         isl_vec *v;
  53         unsigned total;
  54
  55         if (!ls)
  56                 return NULL;
  57
  58         ctx = isl_local_space_get_ctx(ls);
  59         if (!isl_local_space_divs_known(ls))
  60                 isl_die(ctx, isl_error_invalid, "local space has unknown divs",
  61                         goto error);
  62         if (!isl_local_space_is_set(ls))
  63                 isl_die(ctx, isl_error_invalid,
  64                         "domain of affine expression should be a set",
  65                         goto error);
  66
  67         total = isl_local_space_dim(ls, isl_dim_all);
  68         v = isl_vec_alloc(ctx, 1 + 1 + total);
  69         return isl_aff_alloc_vec(ls, v);
  70 error:
  71         isl_local_space_free(ls);
  72         return NULL;
  73 }
  74
  75 __isl_give isl_aff *isl_aff_zero_on_domain(__isl_take isl_local_space *ls)
  76 {
  77         isl_aff *aff;
  78
  79         aff = isl_aff_alloc(ls);
  80         if (!aff)
  81                 return NULL;
  82
  83         isl_int_set_si(aff->v->el[0], 1);
  84         isl_seq_clr(aff->v->el + 1, aff->v->size - 1);
  85
  86         return aff;
  87 }
  88
  89 __isl_give isl_aff *isl_aff_copy(__isl_keep isl_aff *aff)
  90 {
  91         if (!aff)
  92                 return NULL;
  93
  94         aff->ref++;
  95         return aff;
  96 }
  97
  98 __isl_give isl_aff *isl_aff_dup(__isl_keep isl_aff *aff)
  99 {
 100         if (!aff)
 101                 return NULL;
 102
 103         return isl_aff_alloc_vec(isl_local_space_copy(aff->ls),
 104                                  isl_vec_copy(aff->v));
 105 }
 106
 107 __isl_give isl_aff *isl_aff_cow(__isl_take isl_aff *aff)
 108 {
 109         if (!aff)
 110                 return NULL;
 111
 112         if (aff->ref == 1)
 113                 return aff;
 114         aff->ref--;
 115         return isl_aff_dup(aff);
 116 }
 117
 118 void *isl_aff_free(__isl_take isl_aff *aff)
 119 {
 120         if (!aff)
 121                 return NULL;
 122
 123         if (--aff->ref > 0)
 124                 return NULL;
 125
 126         isl_local_space_free(aff->ls);
 127         isl_vec_free(aff->v);
 128
 129         free(aff);
 130
 131         return NULL;
 132 }
 133
 134 isl_ctx *isl_aff_get_ctx(__isl_keep isl_aff *aff)
 135 {
 136         return aff ? isl_local_space_get_ctx(aff->ls) : NULL;
 137 }
 138
 139 /* Externally, an isl_aff has a map space, but internally, the
 140  * ls field corresponds to the domain of that space.
 141  */
 142 int isl_aff_dim(__isl_keep isl_aff *aff, enum isl_dim_type type)
 143 {
 144         if (!aff)
 145                 return 0;
 146         if (type == isl_dim_out)
 147                 return 1;
 148         if (type == isl_dim_in)
 149                 type = isl_dim_set;
 150         return isl_local_space_dim(aff->ls, type);
 151 }
 152
 153 __isl_give isl_space *isl_aff_get_domain_space(__isl_keep isl_aff *aff)
 154 {
 155         return aff ? isl_local_space_get_space(aff->ls) : NULL;
 156 }
 157
 158 __isl_give isl_space *isl_aff_get_space(__isl_keep isl_aff *aff)
 159 {
 160         isl_space *space;
 161         if (!aff)
 162                 return NULL;
 163         space = isl_local_space_get_space(aff->ls);
 164         space = isl_space_from_domain(space);
 165         space = isl_space_add_dims(space, isl_dim_out, 1);
 166         return space;
 167 }
 168
 169 __isl_give isl_local_space *isl_aff_get_domain_local_space(
 170         __isl_keep isl_aff *aff)
 171 {
 172         return aff ? isl_local_space_copy(aff->ls) : NULL;
 173 }
 174
 175 __isl_give isl_local_space *isl_aff_get_local_space(__isl_keep isl_aff *aff)
 176 {
 177         isl_local_space *ls;
 178         if (!aff)
 179                 return NULL;
 180         ls = isl_local_space_copy(aff->ls);
 181         ls = isl_local_space_from_domain(ls);
 182         ls = isl_local_space_add_dims(ls, isl_dim_out, 1);
 183         return ls;
 184 }
 185
 186 /* Externally, an isl_aff has a map space, but internally, the
 187  * ls field corresponds to the domain of that space.
 188  */
 189 const char *isl_aff_get_dim_name(__isl_keep isl_aff *aff,
 190         enum isl_dim_type type, unsigned pos)
 191 {
 192         if (!aff)
 193                 return NULL;
 194         if (type == isl_dim_out)
 195                 return NULL;
 196         if (type == isl_dim_in)
 197                 type = isl_dim_set;
 198         return isl_local_space_get_dim_name(aff->ls, type, pos);
 199 }
 200
 201 __isl_give isl_aff *isl_aff_reset_domain_space(__isl_take isl_aff *aff,
 202         __isl_take isl_space *dim)
 203 {
 204         aff = isl_aff_cow(aff);
 205         if (!aff || !dim)
 206                 goto error;
 207
 208         aff->ls = isl_local_space_reset_space(aff->ls, dim);
 209         if (!aff->ls)
 210                 return isl_aff_free(aff);
 211
 212         return aff;
 213 error:
 214         isl_aff_free(aff);
 215         isl_space_free(dim);
 216         return NULL;
 217 }
 218
 219 /* Reset the space of "aff".  This function is called from isl_pw_templ.c
 220  * and doesn't know if the space of an element object is represented
 221  * directly or through its domain.  It therefore passes along both.
 222  */
 223 __isl_give isl_aff *isl_aff_reset_space_and_domain(__isl_take isl_aff *aff,
 224         __isl_take isl_space *space, __isl_take isl_space *domain)
 225 {
 226         isl_space_free(space);
 227         return isl_aff_reset_domain_space(aff, domain);
 228 }
 229
 230 /* Reorder the coefficients of the affine expression based
 231  * on the given reodering.
 232  * The reordering r is assumed to have been extended with the local
 233  * variables.
 234  */
 235 static __isl_give isl_vec *vec_reorder(__isl_take isl_vec *vec,
 236         __isl_take isl_reordering *r, int n_div)
 237 {
 238         isl_vec *res;
 239         int i;
 240
 241         if (!vec || !r)
 242                 goto error;
 243
 244         res = isl_vec_alloc(vec->ctx,
 245                             2 + isl_space_dim(r->dim, isl_dim_all) + n_div);
 246         isl_seq_cpy(res->el, vec->el, 2);
 247         isl_seq_clr(res->el + 2, res->size - 2);
 248         for (i = 0; i < r->len; ++i)
 249                 isl_int_set(res->el[2 + r->pos[i]], vec->el[2 + i]);
 250
 251         isl_reordering_free(r);
 252         isl_vec_free(vec);
 253         return res;
 254 error:
 255         isl_vec_free(vec);
 256         isl_reordering_free(r);
 257         return NULL;
 258 }
 259
 260 /* Reorder the dimensions of the domain of "aff" according
 261  * to the given reordering.
 262  */
 263 __isl_give isl_aff *isl_aff_realign_domain(__isl_take isl_aff *aff,
 264         __isl_take isl_reordering *r)
 265 {
 266         aff = isl_aff_cow(aff);
 267         if (!aff)
 268                 goto error;
 269
 270         r = isl_reordering_extend(r, aff->ls->div->n_row);
 271         aff->v = vec_reorder(aff->v, isl_reordering_copy(r),
 272                                 aff->ls->div->n_row);
 273         aff->ls = isl_local_space_realign(aff->ls, r);
 274
 275         if (!aff->v || !aff->ls)
 276                 return isl_aff_free(aff);
 277
 278         return aff;
 279 error:
 280         isl_aff_free(aff);
 281         isl_reordering_free(r);
 282         return NULL;
 283 }
 284
 285 __isl_give isl_aff *isl_aff_align_params(__isl_take isl_aff *aff,
 286         __isl_take isl_space *model)
 287 {
 288         if (!aff || !model)
 289                 goto error;
 290
 291         if (!isl_space_match(aff->ls->dim, isl_dim_param,
 292                              model, isl_dim_param)) {
 293                 isl_reordering *exp;
 294
 295                 model = isl_space_drop_dims(model, isl_dim_in,
 296                                         0, isl_space_dim(model, isl_dim_in));
 297                 model = isl_space_drop_dims(model, isl_dim_out,
 298                                         0, isl_space_dim(model, isl_dim_out));
 299                 exp = isl_parameter_alignment_reordering(aff->ls->dim, model);
 300                 exp = isl_reordering_extend_space(exp,
 301                                         isl_aff_get_domain_space(aff));
 302                 aff = isl_aff_realign_domain(aff, exp);
 303         }
 304
 305         isl_space_free(model);
 306         return aff;
 307 error:
 308         isl_space_free(model);
 309         isl_aff_free(aff);
 310         return NULL;
 311 }
 312
 313 int isl_aff_plain_is_zero(__isl_keep isl_aff *aff)
 314 {
 315         if (!aff)
 316                 return -1;
 317
 318         return isl_seq_first_non_zero(aff->v->el + 1, aff->v->size - 1) < 0;
 319 }
 320
 321 int isl_aff_plain_is_equal(__isl_keep isl_aff *aff1, __isl_keep isl_aff *aff2)
 322 {
 323         int equal;
 324
 325         if (!aff1 || !aff2)
 326                 return -1;
 327
 328         equal = isl_local_space_is_equal(aff1->ls, aff2->ls);
 329         if (equal < 0 || !equal)
 330                 return equal;
 331
 332         return isl_vec_is_equal(aff1->v, aff2->v);
 333 }
 334
 335 int isl_aff_get_denominator(__isl_keep isl_aff *aff, isl_int *v)
 336 {
 337         if (!aff)
 338                 return -1;
 339         isl_int_set(*v, aff->v->el[0]);
 340         return 0;
 341 }
 342
 343 int isl_aff_get_constant(__isl_keep isl_aff *aff, isl_int *v)
 344 {
 345         if (!aff)
 346                 return -1;
 347         isl_int_set(*v, aff->v->el[1]);
 348         return 0;
 349 }
 350
 351 int isl_aff_get_coefficient(__isl_keep isl_aff *aff,
 352         enum isl_dim_type type, int pos, isl_int *v)
 353 {
 354         if (!aff)
 355                 return -1;
 356
 357         if (type == isl_dim_out)
 358                 isl_die(aff->v->ctx, isl_error_invalid,
 359                         "output/set dimension does not have a coefficient",
 360                         return -1);
 361         if (type == isl_dim_in)
 362                 type = isl_dim_set;
 363
 364         if (pos >= isl_local_space_dim(aff->ls, type))
 365                 isl_die(aff->v->ctx, isl_error_invalid,
 366                         "position out of bounds", return -1);
 367
 368         pos += isl_local_space_offset(aff->ls, type);
 369         isl_int_set(*v, aff->v->el[1 + pos]);
 370
 371         return 0;
 372 }
 373
 374 __isl_give isl_aff *isl_aff_set_denominator(__isl_take isl_aff *aff, isl_int v)
 375 {
 376         aff = isl_aff_cow(aff);
 377         if (!aff)
 378                 return NULL;
 379
 380         aff->v = isl_vec_cow(aff->v);
 381         if (!aff->v)
 382                 return isl_aff_free(aff);
 383
 384         isl_int_set(aff->v->el[0], v);
 385
 386         return aff;
 387 }
 388
 389 __isl_give isl_aff *isl_aff_set_constant(__isl_take isl_aff *aff, isl_int v)
 390 {
 391         aff = isl_aff_cow(aff);
 392         if (!aff)
 393                 return NULL;
 394
 395         aff->v = isl_vec_cow(aff->v);
 396         if (!aff->v)
 397                 return isl_aff_free(aff);
 398
 399         isl_int_set(aff->v->el[1], v);
 400
 401         return aff;
 402 }
 403
 404 __isl_give isl_aff *isl_aff_add_constant(__isl_take isl_aff *aff, isl_int v)
 405 {
 406         if (isl_int_is_zero(v))
 407                 return aff;
 408
 409         aff = isl_aff_cow(aff);
 410         if (!aff)
 411                 return NULL;
 412
 413         aff->v = isl_vec_cow(aff->v);
 414         if (!aff->v)
 415                 return isl_aff_free(aff);
 416
 417         isl_int_addmul(aff->v->el[1], aff->v->el[0], v);
 418
 419         return aff;
 420 }
 421
 422 __isl_give isl_aff *isl_aff_add_constant_si(__isl_take isl_aff *aff, int v)
 423 {
 424         isl_int t;
 425
 426         isl_int_init(t);
 427         isl_int_set_si(t, v);
 428         aff = isl_aff_add_constant(aff, t);
 429         isl_int_clear(t);
 430
 431         return aff;
 432 }
 433
 434 __isl_give isl_aff *isl_aff_set_constant_si(__isl_take isl_aff *aff, int v)
 435 {
 436         aff = isl_aff_cow(aff);
 437         if (!aff)
 438                 return NULL;
 439
 440         aff->v = isl_vec_cow(aff->v);
 441         if (!aff->v)
 442                 return isl_aff_free(aff);
 443
 444         isl_int_set_si(aff->v->el[1], v);
 445
 446         return aff;
 447 }
 448
 449 __isl_give isl_aff *isl_aff_set_coefficient(__isl_take isl_aff *aff,
 450         enum isl_dim_type type, int pos, isl_int v)
 451 {
 452         if (!aff)
 453                 return NULL;
 454
 455         if (type == isl_dim_out)
 456                 isl_die(aff->v->ctx, isl_error_invalid,
 457                         "output/set dimension does not have a coefficient",
 458                         return isl_aff_free(aff));
 459         if (type == isl_dim_in)
 460                 type = isl_dim_set;
 461
 462         if (pos >= isl_local_space_dim(aff->ls, type))
 463                 isl_die(aff->v->ctx, isl_error_invalid,
 464                         "position out of bounds", return isl_aff_free(aff));
 465
 466         aff = isl_aff_cow(aff);
 467         if (!aff)
 468                 return NULL;
 469
 470         aff->v = isl_vec_cow(aff->v);
 471         if (!aff->v)
 472                 return isl_aff_free(aff);
 473
 474         pos += isl_local_space_offset(aff->ls, type);
 475         isl_int_set(aff->v->el[1 + pos], v);
 476
 477         return aff;
 478 }
 479
 480 __isl_give isl_aff *isl_aff_set_coefficient_si(__isl_take isl_aff *aff,
 481         enum isl_dim_type type, int pos, int v)
 482 {
 483         if (!aff)
 484                 return NULL;
 485
 486         if (type == isl_dim_out)
 487                 isl_die(aff->v->ctx, isl_error_invalid,
 488                         "output/set dimension does not have a coefficient",
 489                         return isl_aff_free(aff));
 490         if (type == isl_dim_in)
 491                 type = isl_dim_set;
 492
 493         if (pos >= isl_local_space_dim(aff->ls, type))
 494                 isl_die(aff->v->ctx, isl_error_invalid,
 495                         "position out of bounds", return isl_aff_free(aff));
 496
 497         aff = isl_aff_cow(aff);
 498         if (!aff)
 499                 return NULL;
 500
 501         aff->v = isl_vec_cow(aff->v);
 502         if (!aff->v)
 503                 return isl_aff_free(aff);
 504
 505         pos += isl_local_space_offset(aff->ls, type);
 506         isl_int_set_si(aff->v->el[1 + pos], v);
 507
 508         return aff;
 509 }
 510
 511 __isl_give isl_aff *isl_aff_add_coefficient(__isl_take isl_aff *aff,
 512         enum isl_dim_type type, int pos, isl_int v)
 513 {
 514         if (!aff)
 515                 return NULL;
 516
 517         if (type == isl_dim_out)
 518                 isl_die(aff->v->ctx, isl_error_invalid,
 519                         "output/set dimension does not have a coefficient",
 520                         return isl_aff_free(aff));
 521         if (type == isl_dim_in)
 522                 type = isl_dim_set;
 523
 524         if (pos >= isl_local_space_dim(aff->ls, type))
 525                 isl_die(aff->v->ctx, isl_error_invalid,
 526                         "position out of bounds", return isl_aff_free(aff));
 527
 528         aff = isl_aff_cow(aff);
 529         if (!aff)
 530                 return NULL;
 531
 532         aff->v = isl_vec_cow(aff->v);
 533         if (!aff->v)
 534                 return isl_aff_free(aff);
 535
 536         pos += isl_local_space_offset(aff->ls, type);
 537         isl_int_addmul(aff->v->el[1 + pos], aff->v->el[0], v);
 538
 539         return aff;
 540 }
 541
 542 __isl_give isl_aff *isl_aff_add_coefficient_si(__isl_take isl_aff *aff,
 543         enum isl_dim_type type, int pos, int v)
 544 {
 545         isl_int t;
 546
 547         isl_int_init(t);
 548         isl_int_set_si(t, v);
 549         aff = isl_aff_add_coefficient(aff, type, pos, t);
 550         isl_int_clear(t);
 551
 552         return aff;
 553 }
 554
 555 __isl_give isl_aff *isl_aff_get_div(__isl_keep isl_aff *aff, int pos)
 556 {
 557         if (!aff)
 558                 return NULL;
 559
 560         return isl_local_space_get_div(aff->ls, pos);
 561 }
 562
 563 __isl_give isl_aff *isl_aff_neg(__isl_take isl_aff *aff)
 564 {
 565         aff = isl_aff_cow(aff);
 566         if (!aff)
 567                 return NULL;
 568         aff->v = isl_vec_cow(aff->v);
 569         if (!aff->v)
 570                 return isl_aff_free(aff);
 571
 572         isl_seq_neg(aff->v->el + 1, aff->v->el + 1, aff->v->size - 1);
 573
 574         return aff;
 575 }
 576
 577 /* Remove divs from the local space that do not appear in the affine
 578  * expression.
 579  * We currently only remove divs at the end.
 580  * Some intermediate divs may also not appear directly in the affine
 581  * expression, but we would also need to check that no other divs are
 582  * defined in terms of them.
 583  */
 584 __isl_give isl_aff *isl_aff_remove_unused_divs( __isl_take isl_aff *aff)
 585 {
 586         int pos;
 587         int off;
 588         int n;
 589
 590         if (!aff)
 591                 return NULL;
 592
 593         n = isl_local_space_dim(aff->ls, isl_dim_div);
 594         off = isl_local_space_offset(aff->ls, isl_dim_div);
 595
 596         pos = isl_seq_last_non_zero(aff->v->el + 1 + off, n) + 1;
 597         if (pos == n)
 598                 return aff;
 599
 600         aff = isl_aff_cow(aff);
 601         if (!aff)
 602                 return NULL;
 603
 604         aff->ls = isl_local_space_drop_dims(aff->ls, isl_dim_div, pos, n - pos);
 605         aff->v = isl_vec_drop_els(aff->v, 1 + off + pos, n - pos);
 606         if (!aff->ls || !aff->v)
 607                 return isl_aff_free(aff);
 608
 609         return aff;
 610 }
 611
 612 __isl_give isl_aff *isl_aff_normalize(__isl_take isl_aff *aff)
 613 {
 614         if (!aff)
 615                 return NULL;
 616         aff->v = isl_vec_normalize(aff->v);
 617         if (!aff->v)
 618                 return isl_aff_free(aff);
 619         aff = isl_aff_remove_unused_divs(aff);
 620         return aff;
 621 }
 622
 623 /* Given f, return floor(f).
 624  * If f is an integer expression, then just return f.
 625  * Otherwise, if f = g/m, write g = q m + r,
 626  * create a new div d = [r/m] and return the expression q + d.
 627  * The coefficients in r are taken to lie between -m/2 and m/2.
 628  */
 629 __isl_give isl_aff *isl_aff_floor(__isl_take isl_aff *aff)
 630 {
 631         int i;
 632         int size;
 633         isl_ctx *ctx;
 634         isl_vec *div;
 635
 636         if (!aff)
 637                 return NULL;
 638
 639         if (isl_int_is_one(aff->v->el[0]))
 640                 return aff;
 641
 642         aff = isl_aff_cow(aff);
 643         if (!aff)
 644                 return NULL;
 645
 646         aff->v = isl_vec_cow(aff->v);
 647         div = isl_vec_copy(aff->v);
 648         div = isl_vec_cow(div);
 649         if (!div)
 650                 return isl_aff_free(aff);
 651
 652         ctx = isl_aff_get_ctx(aff);
 653         isl_int_fdiv_q(aff->v->el[0], aff->v->el[0], ctx->two);
 654         for (i = 1; i < aff->v->size; ++i) {
 655                 isl_int_fdiv_r(div->el[i], div->el[i], div->el[0]);
 656                 isl_int_fdiv_q(aff->v->el[i], aff->v->el[i], div->el[0]);
 657                 if (isl_int_gt(div->el[i], aff->v->el[0])) {
 658                         isl_int_sub(div->el[i], div->el[i], div->el[0]);
 659                         isl_int_add_ui(aff->v->el[i], aff->v->el[i], 1);
 660                 }
 661         }
 662
 663         aff->ls = isl_local_space_add_div(aff->ls, div);
 664         if (!aff->ls)
 665                 return isl_aff_free(aff);
 666
 667         size = aff->v->size;
 668         aff->v = isl_vec_extend(aff->v, size + 1);
 669         if (!aff->v)
 670                 return isl_aff_free(aff);
 671         isl_int_set_si(aff->v->el[0], 1);
 672         isl_int_set_si(aff->v->el[size], 1);
 673
 674         return aff;
 675 }
 676
 677 /* Compute
 678  *
 679  *      aff mod m = aff - m * floor(aff/m)
 680  */
 681 __isl_give isl_aff *isl_aff_mod(__isl_take isl_aff *aff, isl_int m)
 682 {
 683         isl_aff *res;
 684
 685         res = isl_aff_copy(aff);
 686         aff = isl_aff_scale_down(aff, m);
 687         aff = isl_aff_floor(aff);
 688         aff = isl_aff_scale(aff, m);
 689         res = isl_aff_sub(res, aff);
 690
 691         return res;
 692 }
 693
 694 /* Compute
 695  *
 696  *      pwaff mod m = pwaff - m * floor(pwaff/m)
 697  */
 698 __isl_give isl_pw_aff *isl_pw_aff_mod(__isl_take isl_pw_aff *pwaff, isl_int m)
 699 {
 700         isl_pw_aff *res;
 701
 702         res = isl_pw_aff_copy(pwaff);
 703         pwaff = isl_pw_aff_scale_down(pwaff, m);
 704         pwaff = isl_pw_aff_floor(pwaff);
 705         pwaff = isl_pw_aff_scale(pwaff, m);
 706         res = isl_pw_aff_sub(res, pwaff);
 707
 708         return res;
 709 }
 710
 711 /* Given f, return ceil(f).
 712  * If f is an integer expression, then just return f.
 713  * Otherwise, create a new div d = [-f] and return the expression -d.
 714  */
 715 __isl_give isl_aff *isl_aff_ceil(__isl_take isl_aff *aff)
 716 {
 717         if (!aff)
 718                 return NULL;
 719
 720         if (isl_int_is_one(aff->v->el[0]))
 721                 return aff;
 722
 723         aff = isl_aff_neg(aff);
 724         aff = isl_aff_floor(aff);
 725         aff = isl_aff_neg(aff);
 726
 727         return aff;
 728 }
 729
 730 /* Apply the expansion computed by isl_merge_divs.
 731  * The expansion itself is given by "exp" while the resulting
 732  * list of divs is given by "div".
 733  */
 734 __isl_give isl_aff *isl_aff_expand_divs( __isl_take isl_aff *aff,
 735         __isl_take isl_mat *div, int *exp)
 736 {
 737         int i, j;
 738         int old_n_div;
 739         int new_n_div;
 740         int offset;
 741
 742         aff = isl_aff_cow(aff);
 743         if (!aff || !div)
 744                 goto error;
 745
 746         old_n_div = isl_local_space_dim(aff->ls, isl_dim_div);
 747         new_n_div = isl_mat_rows(div);
 748         if (new_n_div < old_n_div)
 749                 isl_die(isl_mat_get_ctx(div), isl_error_invalid,
 750                         "not an expansion", goto error);
 751
 752         aff->v = isl_vec_extend(aff->v, aff->v->size + new_n_div - old_n_div);
 753         if (!aff->v)
 754                 goto error;
 755
 756         offset = 1 + isl_local_space_offset(aff->ls, isl_dim_div);
 757         j = old_n_div - 1;
 758         for (i = new_n_div - 1; i >= 0; --i) {
 759                 if (j >= 0 && exp[j] == i) {
 760                         if (i != j)
 761                                 isl_int_swap(aff->v->el[offset + i],
 762                                              aff->v->el[offset + j]);
 763                         j--;
 764                 } else
 765                         isl_int_set_si(aff->v->el[offset + i], 0);
 766         }
 767
 768         aff->ls = isl_local_space_replace_divs(aff->ls, isl_mat_copy(div));
 769         if (!aff->ls)
 770                 goto error;
 771         isl_mat_free(div);
 772         return aff;
 773 error:
 774         isl_aff_free(aff);
 775         isl_mat_free(div);
 776         return NULL;
 777 }
 778
 779 /* Add two affine expressions that live in the same local space.
 780  */
 781 static __isl_give isl_aff *add_expanded(__isl_take isl_aff *aff1,
 782         __isl_take isl_aff *aff2)
 783 {
 784         isl_int gcd, f;
 785
 786         aff1 = isl_aff_cow(aff1);
 787         if (!aff1 || !aff2)
 788                 goto error;
 789
 790         aff1->v = isl_vec_cow(aff1->v);
 791         if (!aff1->v)
 792                 goto error;
 793
 794         isl_int_init(gcd);
 795         isl_int_init(f);
 796         isl_int_gcd(gcd, aff1->v->el[0], aff2->v->el[0]);
 797         isl_int_divexact(f, aff2->v->el[0], gcd);
 798         isl_seq_scale(aff1->v->el + 1, aff1->v->el + 1, f, aff1->v->size - 1);
 799         isl_int_divexact(f, aff1->v->el[0], gcd);
 800         isl_seq_addmul(aff1->v->el + 1, f, aff2->v->el + 1, aff1->v->size - 1);
 801         isl_int_divexact(f, aff2->v->el[0], gcd);
 802         isl_int_mul(aff1->v->el[0], aff1->v->el[0], f);
 803         isl_int_clear(f);
 804         isl_int_clear(gcd);
 805
 806         isl_aff_free(aff2);
 807         return aff1;
 808 error:
 809         isl_aff_free(aff1);
 810         isl_aff_free(aff2);
 811         return NULL;
 812 }
 813
 814 __isl_give isl_aff *isl_aff_add(__isl_take isl_aff *aff1,
 815         __isl_take isl_aff *aff2)
 816 {
 817         isl_ctx *ctx;
 818         int *exp1 = NULL;
 819         int *exp2 = NULL;
 820         isl_mat *div;
 821
 822         if (!aff1 || !aff2)
 823                 goto error;
 824
 825         ctx = isl_aff_get_ctx(aff1);
 826         if (!isl_space_is_equal(aff1->ls->dim, aff2->ls->dim))
 827                 isl_die(ctx, isl_error_invalid,
 828                         "spaces don't match", goto error);
 829
 830         if (aff1->ls->div->n_row == 0 && aff2->ls->div->n_row == 0)
 831                 return add_expanded(aff1, aff2);
 832
 833         exp1 = isl_alloc_array(ctx, int, aff1->ls->div->n_row);
 834         exp2 = isl_alloc_array(ctx, int, aff2->ls->div->n_row);
 835         if (!exp1 || !exp2)
 836                 goto error;
 837
 838         div = isl_merge_divs(aff1->ls->div, aff2->ls->div, exp1, exp2);
 839         aff1 = isl_aff_expand_divs(aff1, isl_mat_copy(div), exp1);
 840         aff2 = isl_aff_expand_divs(aff2, div, exp2);
 841         free(exp1);
 842         free(exp2);
 843
 844         return add_expanded(aff1, aff2);
 845 error:
 846         free(exp1);
 847         free(exp2);
 848         isl_aff_free(aff1);
 849         isl_aff_free(aff2);
 850         return NULL;
 851 }
 852
 853 __isl_give isl_aff *isl_aff_sub(__isl_take isl_aff *aff1,
 854         __isl_take isl_aff *aff2)
 855 {
 856         return isl_aff_add(aff1, isl_aff_neg(aff2));
 857 }
 858
 859 __isl_give isl_aff *isl_aff_scale(__isl_take isl_aff *aff, isl_int f)
 860 {
 861         isl_int gcd;
 862
 863         if (isl_int_is_one(f))
 864                 return aff;
 865
 866         aff = isl_aff_cow(aff);
 867         if (!aff)
 868                 return NULL;
 869         aff->v = isl_vec_cow(aff->v);
 870         if (!aff->v)
 871                 return isl_aff_free(aff);
 872
 873         isl_int_init(gcd);
 874         isl_int_gcd(gcd, aff->v->el[0], f);
 875         isl_int_divexact(aff->v->el[0], aff->v->el[0], gcd);
 876         isl_int_divexact(gcd, f, gcd);
 877         isl_seq_scale(aff->v->el + 1, aff->v->el + 1, gcd, aff->v->size - 1);
 878         isl_int_clear(gcd);
 879
 880         return aff;
 881 }
 882
 883 __isl_give isl_aff *isl_aff_scale_down(__isl_take isl_aff *aff, isl_int f)
 884 {
 885         isl_int gcd;
 886
 887         if (isl_int_is_one(f))
 888                 return aff;
 889
 890         aff = isl_aff_cow(aff);
 891         if (!aff)
 892                 return NULL;
 893         aff->v = isl_vec_cow(aff->v);
 894         if (!aff->v)
 895                 return isl_aff_free(aff);
 896
 897         isl_int_init(gcd);
 898         isl_seq_gcd(aff->v->el + 1, aff->v->size - 1, &gcd);
 899         isl_int_gcd(gcd, gcd, f);
 900         isl_seq_scale_down(aff->v->el + 1, aff->v->el + 1, gcd, aff->v->size - 1);
 901         isl_int_divexact(gcd, f, gcd);
 902         isl_int_mul(aff->v->el[0], aff->v->el[0], gcd);
 903         isl_int_clear(gcd);
 904
 905         return aff;
 906 }
 907
 908 __isl_give isl_aff *isl_aff_scale_down_ui(__isl_take isl_aff *aff, unsigned f)
 909 {
 910         isl_int v;
 911
 912         if (f == 1)
 913                 return aff;
 914
 915         isl_int_init(v);
 916         isl_int_set_ui(v, f);
 917         aff = isl_aff_scale_down(aff, v);
 918         isl_int_clear(v);
 919
 920         return aff;
 921 }
 922
 923 __isl_give isl_aff *isl_aff_set_dim_name(__isl_take isl_aff *aff,
 924         enum isl_dim_type type, unsigned pos, const char *s)
 925 {
 926         aff = isl_aff_cow(aff);
 927         if (!aff)
 928                 return NULL;
 929         if (type == isl_dim_out)
 930                 isl_die(aff->v->ctx, isl_error_invalid,
 931                         "cannot set name of output/set dimension",
 932                         return isl_aff_free(aff));
 933         if (type == isl_dim_in)
 934                 type = isl_dim_set;
 935         aff->ls = isl_local_space_set_dim_name(aff->ls, type, pos, s);
 936         if (!aff->ls)
 937                 return isl_aff_free(aff);
 938
 939         return aff;
 940 }
 941
 942 __isl_give isl_aff *isl_aff_set_dim_id(__isl_take isl_aff *aff,
 943         enum isl_dim_type type, unsigned pos, __isl_take isl_id *id)
 944 {
 945         aff = isl_aff_cow(aff);
 946         if (!aff)
 947                 return isl_id_free(id);
 948         if (type == isl_dim_out)
 949                 isl_die(aff->v->ctx, isl_error_invalid,
 950                         "cannot set name of output/set dimension",
 951                         goto error);
 952         if (type == isl_dim_in)
 953                 type = isl_dim_set;
 954         aff->ls = isl_local_space_set_dim_id(aff->ls, type, pos, id);
 955         if (!aff->ls)
 956                 return isl_aff_free(aff);
 957
 958         return aff;
 959 error:
 960         isl_id_free(id);
 961         isl_aff_free(aff);
 962         return NULL;
 963 }
 964
 965 /* Exploit the equalities in "eq" to simplify the affine expression
 966  * and the expressions of the integer divisions in the local space.
 967  * The integer divisions in this local space are assumed to appear
 968  * as regular dimensions in "eq".
 969  */
 970 static __isl_give isl_aff *isl_aff_substitute_equalities_lifted(
 971         __isl_take isl_aff *aff, __isl_take isl_basic_set *eq)
 972 {
 973         int i, j;
 974         unsigned total;
 975         unsigned n_div;
 976
 977         if (!eq)
 978                 goto error;
 979         if (eq->n_eq == 0) {
 980                 isl_basic_set_free(eq);
 981                 return aff;
 982         }
 983
 984         aff = isl_aff_cow(aff);
 985         if (!aff)
 986                 goto error;
 987
 988         aff->ls = isl_local_space_substitute_equalities(aff->ls,
 989                                                         isl_basic_set_copy(eq));
 990         if (!aff->ls)
 991                 goto error;
 992
 993         total = 1 + isl_space_dim(eq->dim, isl_dim_all);
 994         n_div = eq->n_div;
 995         for (i = 0; i < eq->n_eq; ++i) {
 996                 j = isl_seq_last_non_zero(eq->eq[i], total + n_div);
 997                 if (j < 0 || j == 0 || j >= total)
 998                         continue;
 999
1000                 isl_seq_elim(aff->v->el + 1, eq->eq[i], j, total,
1001                                 &aff->v->el[0]);
1002         }
1003
1004         isl_basic_set_free(eq);
1005         aff = isl_aff_normalize(aff);
1006         return aff;
1007 error:
1008         isl_basic_set_free(eq);
1009         isl_aff_free(aff);
1010         return NULL;
1011 }
1012
1013 /* Exploit the equalities in "eq" to simplify the affine expression
1014  * and the expressions of the integer divisions in the local space.
1015  */
1016 static __isl_give isl_aff *isl_aff_substitute_equalities(
1017         __isl_take isl_aff *aff, __isl_take isl_basic_set *eq)
1018 {
1019         int n_div;
1020
1021         if (!aff || !eq)
1022                 goto error;
1023         n_div = isl_local_space_dim(aff->ls, isl_dim_div);
1024         if (n_div > 0)
1025                 eq = isl_basic_set_add(eq, isl_dim_set, n_div);
1026         return isl_aff_substitute_equalities_lifted(aff, eq);
1027 error:
1028         isl_basic_set_free(eq);
1029         isl_aff_free(aff);
1030         return NULL;
1031 }
1032
1033 /* Look for equalities among the variables shared by context and aff
1034  * and the integer divisions of aff, if any.
1035  * The equalities are then used to eliminate coefficients and/or integer
1036  * divisions from aff.
1037  */
1038 __isl_give isl_aff *isl_aff_gist(__isl_take isl_aff *aff,
1039         __isl_take isl_set *context)
1040 {
1041         isl_basic_set *hull;
1042         int n_div;
1043
1044         if (!aff)
1045                 goto error;
1046         n_div = isl_local_space_dim(aff->ls, isl_dim_div);
1047         if (n_div > 0) {
1048                 isl_basic_set *bset;
1049                 isl_local_space *ls;
1050                 context = isl_set_add_dims(context, isl_dim_set, n_div);
1051                 ls = isl_aff_get_domain_local_space(aff);
1052                 bset = isl_basic_set_from_local_space(ls);
1053                 bset = isl_basic_set_lift(bset);
1054                 bset = isl_basic_set_flatten(bset);
1055                 context = isl_set_intersect(context,
1056                                             isl_set_from_basic_set(bset));
1057         }
1058
1059         hull = isl_set_affine_hull(context);
1060         return isl_aff_substitute_equalities_lifted(aff, hull);
1061 error:
1062         isl_aff_free(aff);
1063         isl_set_free(context);
1064         return NULL;
1065 }
1066
1067 __isl_give isl_aff *isl_aff_gist_params(__isl_take isl_aff *aff,
1068         __isl_take isl_set *context)
1069 {
1070         isl_set *dom_context = isl_set_universe(isl_aff_get_domain_space(aff));
1071         dom_context = isl_set_intersect_params(dom_context, context);
1072         return isl_aff_gist(aff, dom_context);
1073 }
1074
1075 /* Return a basic set containing those elements in the space
1076  * of aff where it is non-negative.
1077  */
1078 __isl_give isl_basic_set *isl_aff_nonneg_basic_set(__isl_take isl_aff *aff)
1079 {
1080         isl_constraint *ineq;
1081         isl_basic_set *bset;
1082
1083         ineq = isl_inequality_from_aff(aff);
1084
1085         bset = isl_basic_set_from_constraint(ineq);
1086         bset = isl_basic_set_simplify(bset);
1087         return bset;
1088 }
1089
1090 /* Return a basic set containing those elements in the space
1091  * of aff where it is zero.
1092  */
1093 __isl_give isl_basic_set *isl_aff_zero_basic_set(__isl_take isl_aff *aff)
1094 {
1095         isl_constraint *ineq;
1096         isl_basic_set *bset;
1097
1098         ineq = isl_equality_from_aff(aff);
1099
1100         bset = isl_basic_set_from_constraint(ineq);
1101         bset = isl_basic_set_simplify(bset);
1102         return bset;
1103 }
1104
1105 /* Return a basic set containing those elements in the shared space
1106  * of aff1 and aff2 where aff1 is greater than or equal to aff2.
1107  */
1108 __isl_give isl_basic_set *isl_aff_ge_basic_set(__isl_take isl_aff *aff1,
1109         __isl_take isl_aff *aff2)
1110 {
1111         aff1 = isl_aff_sub(aff1, aff2);
1112
1113         return isl_aff_nonneg_basic_set(aff1);
1114 }
1115
1116 /* Return a basic set containing those elements in the shared space
1117  * of aff1 and aff2 where aff1 is smaller than or equal to aff2.
1118  */
1119 __isl_give isl_basic_set *isl_aff_le_basic_set(__isl_take isl_aff *aff1,
1120         __isl_take isl_aff *aff2)
1121 {
1122         return isl_aff_ge_basic_set(aff2, aff1);
1123 }
1124
1125 __isl_give isl_aff *isl_aff_add_on_domain(__isl_keep isl_set *dom,
1126         __isl_take isl_aff *aff1, __isl_take isl_aff *aff2)
1127 {
1128         aff1 = isl_aff_add(aff1, aff2);
1129         aff1 = isl_aff_gist(aff1, isl_set_copy(dom));
1130         return aff1;
1131 }
1132
1133 int isl_aff_is_empty(__isl_keep isl_aff *aff)
1134 {
1135         if (!aff)
1136                 return -1;
1137
1138         return 0;
1139 }
1140
1141 /* Check whether the given affine expression has non-zero coefficient
1142  * for any dimension in the given range or if any of these dimensions
1143  * appear with non-zero coefficients in any of the integer divisions
1144  * involved in the affine expression.
1145  */
1146 int isl_aff_involves_dims(__isl_keep isl_aff *aff,
1147         enum isl_dim_type type, unsigned first, unsigned n)
1148 {
1149         int i;
1150         isl_ctx *ctx;
1151         int *active = NULL;
1152         int involves = 0;
1153
1154         if (!aff)
1155                 return -1;
1156         if (n == 0)
1157                 return 0;
1158
1159         ctx = isl_aff_get_ctx(aff);
1160         if (first + n > isl_aff_dim(aff, type))
1161                 isl_die(ctx, isl_error_invalid,
1162                         "range out of bounds", return -1);
1163
1164         active = isl_local_space_get_active(aff->ls, aff->v->el + 2);
1165         if (!active)
1166                 goto error;
1167
1168         first += isl_local_space_offset(aff->ls, type) - 1;
1169         for (i = 0; i < n; ++i)
1170                 if (active[first + i]) {
1171                         involves = 1;
1172                         break;
1173                 }
1174
1175         free(active);
1176
1177         return involves;
1178 error:
1179         free(active);
1180         return -1;
1181 }
1182
1183 __isl_give isl_aff *isl_aff_drop_dims(__isl_take isl_aff *aff,
1184         enum isl_dim_type type, unsigned first, unsigned n)
1185 {
1186         isl_ctx *ctx;
1187
1188         if (!aff)
1189                 return NULL;
1190         if (type == isl_dim_out)
1191                 isl_die(aff->v->ctx, isl_error_invalid,
1192                         "cannot drop output/set dimension",
1193                         return isl_aff_free(aff));
1194         if (type == isl_dim_in)
1195                 type = isl_dim_set;
1196         if (n == 0 && !isl_local_space_is_named_or_nested(aff->ls, type))
1197                 return aff;
1198
1199         ctx = isl_aff_get_ctx(aff);
1200         if (first + n > isl_local_space_dim(aff->ls, type))
1201                 isl_die(ctx, isl_error_invalid, "range out of bounds",
1202                         return isl_aff_free(aff));
1203
1204         aff = isl_aff_cow(aff);
1205         if (!aff)
1206                 return NULL;
1207
1208         aff->ls = isl_local_space_drop_dims(aff->ls, type, first, n);
1209         if (!aff->ls)
1210                 return isl_aff_free(aff);
1211
1212         first += 1 + isl_local_space_offset(aff->ls, type);
1213         aff->v = isl_vec_drop_els(aff->v, first, n);
1214         if (!aff->v)
1215                 return isl_aff_free(aff);
1216
1217         return aff;
1218 }
1219
1220 /* Project the domain of the affine expression onto its parameter space.
1221  * The affine expression may not involve any of the domain dimensions.
1222  */
1223 __isl_give isl_aff *isl_aff_project_domain_on_params(__isl_take isl_aff *aff)
1224 {
1225         isl_space *space;
1226         unsigned n;
1227         int involves;
1228
1229         n = isl_aff_dim(aff, isl_dim_in);
1230         involves = isl_aff_involves_dims(aff, isl_dim_in, 0, n);
1231         if (involves < 0)
1232                 return isl_aff_free(aff);
1233         if (involves)
1234                 isl_die(isl_aff_get_ctx(aff), isl_error_invalid,
1235                     "affine expression involves some of the domain dimensions",
1236                     return isl_aff_free(aff));
1237         aff = isl_aff_drop_dims(aff, isl_dim_in, 0, n);
1238         space = isl_aff_get_domain_space(aff);
1239         space = isl_space_params(space);
1240         aff = isl_aff_reset_domain_space(aff, space);
1241         return aff;
1242 }
1243
1244 __isl_give isl_aff *isl_aff_insert_dims(__isl_take isl_aff *aff,
1245         enum isl_dim_type type, unsigned first, unsigned n)
1246 {
1247         isl_ctx *ctx;
1248
1249         if (!aff)
1250                 return NULL;
1251         if (type == isl_dim_out)
1252                 isl_die(aff->v->ctx, isl_error_invalid,
1253                         "cannot insert output/set dimensions",
1254                         return isl_aff_free(aff));
1255         if (type == isl_dim_in)
1256                 type = isl_dim_set;
1257         if (n == 0 && !isl_local_space_is_named_or_nested(aff->ls, type))
1258                 return aff;
1259
1260         ctx = isl_aff_get_ctx(aff);
1261         if (first > isl_local_space_dim(aff->ls, type))
1262                 isl_die(ctx, isl_error_invalid, "position out of bounds",
1263                         return isl_aff_free(aff));
1264
1265         aff = isl_aff_cow(aff);
1266         if (!aff)
1267                 return NULL;
1268
1269         aff->ls = isl_local_space_insert_dims(aff->ls, type, first, n);
1270         if (!aff->ls)
1271                 return isl_aff_free(aff);
1272
1273         first += 1 + isl_local_space_offset(aff->ls, type);
1274         aff->v = isl_vec_insert_zero_els(aff->v, first, n);
1275         if (!aff->v)
1276                 return isl_aff_free(aff);
1277
1278         return aff;
1279 }
1280
1281 __isl_give isl_aff *isl_aff_add_dims(__isl_take isl_aff *aff,
1282         enum isl_dim_type type, unsigned n)
1283 {
1284         unsigned pos;
1285
1286         pos = isl_aff_dim(aff, type);
1287
1288         return isl_aff_insert_dims(aff, type, pos, n);
1289 }
1290
1291 __isl_give isl_pw_aff *isl_pw_aff_add_dims(__isl_take isl_pw_aff *pwaff,
1292         enum isl_dim_type type, unsigned n)
1293 {
1294         unsigned pos;
1295
1296         pos = isl_pw_aff_dim(pwaff, type);
1297
1298         return isl_pw_aff_insert_dims(pwaff, type, pos, n);
1299 }
1300
1301 __isl_give isl_pw_aff *isl_pw_aff_from_aff(__isl_take isl_aff *aff)
1302 {
1303         isl_set *dom = isl_set_universe(isl_aff_get_domain_space(aff));
1304         return isl_pw_aff_alloc(dom, aff);
1305 }
1306
1307 #undef PW
1308 #define PW isl_pw_aff
1309 #undef EL
1310 #define EL isl_aff
1311 #undef EL_IS_ZERO
1312 #define EL_IS_ZERO is_empty
1313 #undef ZERO
1314 #define ZERO empty
1315 #undef IS_ZERO
1316 #define IS_ZERO is_empty
1317 #undef FIELD
1318 #define FIELD aff
1319 #undef DEFAULT_IS_ZERO
1320 #define DEFAULT_IS_ZERO 0
1321
1322 #define NO_EVAL
1323 #define NO_OPT
1324 #define NO_MOVE_DIMS
1325 #define NO_LIFT
1326 #define NO_MORPH
1327
1328 #include <isl_pw_templ.c>
1329
1330 static __isl_give isl_set *align_params_pw_pw_set_and(
1331         __isl_take isl_pw_aff *pwaff1, __isl_take isl_pw_aff *pwaff2,
1332         __isl_give isl_set *(*fn)(__isl_take isl_pw_aff *pwaff1,
1333                                     __isl_take isl_pw_aff *pwaff2))
1334 {
1335         if (!pwaff1 || !pwaff2)
1336                 goto error;
1337         if (isl_space_match(pwaff1->dim, isl_dim_param,
1338                           pwaff2->dim, isl_dim_param))
1339                 return fn(pwaff1, pwaff2);
1340         if (!isl_space_has_named_params(pwaff1->dim) ||
1341             !isl_space_has_named_params(pwaff2->dim))
1342                 isl_die(isl_pw_aff_get_ctx(pwaff1), isl_error_invalid,
1343                         "unaligned unnamed parameters", goto error);
1344         pwaff1 = isl_pw_aff_align_params(pwaff1, isl_pw_aff_get_space(pwaff2));
1345         pwaff2 = isl_pw_aff_align_params(pwaff2, isl_pw_aff_get_space(pwaff1));
1346         return fn(pwaff1, pwaff2);
1347 error:
1348         isl_pw_aff_free(pwaff1);
1349         isl_pw_aff_free(pwaff2);
1350         return NULL;
1351 }
1352
1353 /* Compute a piecewise quasi-affine expression with a domain that
1354  * is the union of those of pwaff1 and pwaff2 and such that on each
1355  * cell, the quasi-affine expression is the better (according to cmp)
1356  * of those of pwaff1 and pwaff2.  If only one of pwaff1 or pwaff2
1357  * is defined on a given cell, then the associated expression
1358  * is the defined one.
1359  */
1360 static __isl_give isl_pw_aff *pw_aff_union_opt(__isl_take isl_pw_aff *pwaff1,
1361         __isl_take isl_pw_aff *pwaff2,
1362         __isl_give isl_basic_set *(*cmp)(__isl_take isl_aff *aff1,
1363                                         __isl_take isl_aff *aff2))
1364 {
1365         int i, j, n;
1366         isl_pw_aff *res;
1367         isl_ctx *ctx;
1368         isl_set *set;
1369
1370         if (!pwaff1 || !pwaff2)
1371                 goto error;
1372
1373         ctx = isl_space_get_ctx(pwaff1->dim);
1374         if (!isl_space_is_equal(pwaff1->dim, pwaff2->dim))
1375                 isl_die(ctx, isl_error_invalid,
1376                         "arguments should live in same space", goto error);
1377
1378         if (isl_pw_aff_is_empty(pwaff1)) {
1379                 isl_pw_aff_free(pwaff1);
1380                 return pwaff2;
1381         }
1382
1383         if (isl_pw_aff_is_empty(pwaff2)) {
1384                 isl_pw_aff_free(pwaff2);
1385                 return pwaff1;
1386         }
1387
1388         n = 2 * (pwaff1->n + 1) * (pwaff2->n + 1);
1389         res = isl_pw_aff_alloc_size(isl_space_copy(pwaff1->dim), n);
1390
1391         for (i = 0; i < pwaff1->n; ++i) {
1392                 set = isl_set_copy(pwaff1->p[i].set);
1393                 for (j = 0; j < pwaff2->n; ++j) {
1394                         struct isl_set *common;
1395                         isl_set *better;
1396
1397                         common = isl_set_intersect(
1398                                         isl_set_copy(pwaff1->p[i].set),
1399                                         isl_set_copy(pwaff2->p[j].set));
1400                         better = isl_set_from_basic_set(cmp(
1401                                         isl_aff_copy(pwaff2->p[j].aff),
1402                                         isl_aff_copy(pwaff1->p[i].aff)));
1403                         better = isl_set_intersect(common, better);
1404                         if (isl_set_plain_is_empty(better)) {
1405                                 isl_set_free(better);
1406                                 continue;
1407                         }
1408                         set = isl_set_subtract(set, isl_set_copy(better));
1409
1410                         res = isl_pw_aff_add_piece(res, better,
1411                                                 isl_aff_copy(pwaff2->p[j].aff));
1412                 }
1413                 res = isl_pw_aff_add_piece(res, set,
1414                                                 isl_aff_copy(pwaff1->p[i].aff));
1415         }
1416
1417         for (j = 0; j < pwaff2->n; ++j) {
1418                 set = isl_set_copy(pwaff2->p[j].set);
1419                 for (i = 0; i < pwaff1->n; ++i)
1420                         set = isl_set_subtract(set,
1421                                         isl_set_copy(pwaff1->p[i].set));
1422                 res = isl_pw_aff_add_piece(res, set,
1423                                                 isl_aff_copy(pwaff2->p[j].aff));
1424         }
1425
1426         isl_pw_aff_free(pwaff1);
1427         isl_pw_aff_free(pwaff2);
1428
1429         return res;
1430 error:
1431         isl_pw_aff_free(pwaff1);
1432         isl_pw_aff_free(pwaff2);
1433         return NULL;
1434 }
1435
1436 /* Compute a piecewise quasi-affine expression with a domain that
1437  * is the union of those of pwaff1 and pwaff2 and such that on each
1438  * cell, the quasi-affine expression is the maximum of those of pwaff1
1439  * and pwaff2.  If only one of pwaff1 or pwaff2 is defined on a given
1440  * cell, then the associated expression is the defined one.
1441  */
1442 static __isl_give isl_pw_aff *pw_aff_union_max(__isl_take isl_pw_aff *pwaff1,
1443         __isl_take isl_pw_aff *pwaff2)
1444 {
1445         return pw_aff_union_opt(pwaff1, pwaff2, &isl_aff_ge_basic_set);
1446 }
1447
1448 __isl_give isl_pw_aff *isl_pw_aff_union_max(__isl_take isl_pw_aff *pwaff1,
1449         __isl_take isl_pw_aff *pwaff2)
1450 {
1451         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2,
1452                                                         &pw_aff_union_max);
1453 }
1454
1455 /* Compute a piecewise quasi-affine expression with a domain that
1456  * is the union of those of pwaff1 and pwaff2 and such that on each
1457  * cell, the quasi-affine expression is the minimum of those of pwaff1
1458  * and pwaff2.  If only one of pwaff1 or pwaff2 is defined on a given
1459  * cell, then the associated expression is the defined one.
1460  */
1461 static __isl_give isl_pw_aff *pw_aff_union_min(__isl_take isl_pw_aff *pwaff1,
1462         __isl_take isl_pw_aff *pwaff2)
1463 {
1464         return pw_aff_union_opt(pwaff1, pwaff2, &isl_aff_le_basic_set);
1465 }
1466
1467 __isl_give isl_pw_aff *isl_pw_aff_union_min(__isl_take isl_pw_aff *pwaff1,
1468         __isl_take isl_pw_aff *pwaff2)
1469 {
1470         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2,
1471                                                         &pw_aff_union_min);
1472 }
1473
1474 __isl_give isl_pw_aff *isl_pw_aff_union_opt(__isl_take isl_pw_aff *pwaff1,
1475         __isl_take isl_pw_aff *pwaff2, int max)
1476 {
1477         if (max)
1478                 return isl_pw_aff_union_max(pwaff1, pwaff2);
1479         else
1480                 return isl_pw_aff_union_min(pwaff1, pwaff2);
1481 }
1482
1483 /* Construct a map with as domain the domain of pwaff and
1484  * one-dimensional range corresponding to the affine expressions.
1485  */
1486 static __isl_give isl_map *map_from_pw_aff(__isl_take isl_pw_aff *pwaff)
1487 {
1488         int i;
1489         isl_space *dim;
1490         isl_map *map;
1491
1492         if (!pwaff)
1493                 return NULL;
1494
1495         dim = isl_pw_aff_get_space(pwaff);
1496         map = isl_map_empty(dim);
1497
1498         for (i = 0; i < pwaff->n; ++i) {
1499                 isl_basic_map *bmap;
1500                 isl_map *map_i;
1501
1502                 bmap = isl_basic_map_from_aff(isl_aff_copy(pwaff->p[i].aff));
1503                 map_i = isl_map_from_basic_map(bmap);
1504                 map_i = isl_map_intersect_domain(map_i,
1505                                                 isl_set_copy(pwaff->p[i].set));
1506                 map = isl_map_union_disjoint(map, map_i);
1507         }
1508
1509         isl_pw_aff_free(pwaff);
1510
1511         return map;
1512 }
1513
1514 /* Construct a map with as domain the domain of pwaff and
1515  * one-dimensional range corresponding to the affine expressions.
1516  */
1517 __isl_give isl_map *isl_map_from_pw_aff(__isl_take isl_pw_aff *pwaff)
1518 {
1519         if (!pwaff)
1520                 return NULL;
1521         if (isl_space_is_set(pwaff->dim))
1522                 isl_die(isl_pw_aff_get_ctx(pwaff), isl_error_invalid,
1523                         "space of input is not a map",
1524                         return isl_pw_aff_free(pwaff));
1525         return map_from_pw_aff(pwaff);
1526 }
1527
1528 /* Construct a one-dimensional set with as parameter domain
1529  * the domain of pwaff and the single set dimension
1530  * corresponding to the affine expressions.
1531  */
1532 __isl_give isl_set *isl_set_from_pw_aff(__isl_take isl_pw_aff *pwaff)
1533 {
1534         if (!pwaff)
1535                 return NULL;
1536         if (!isl_space_is_set(pwaff->dim))
1537                 isl_die(isl_pw_aff_get_ctx(pwaff), isl_error_invalid,
1538                         "space of input is not a set",
1539                         return isl_pw_aff_free(pwaff));
1540         return map_from_pw_aff(pwaff);
1541 }
1542
1543 /* Return a set containing those elements in the domain
1544  * of pwaff where it is non-negative.
1545  */
1546 __isl_give isl_set *isl_pw_aff_nonneg_set(__isl_take isl_pw_aff *pwaff)
1547 {
1548         int i;
1549         isl_set *set;
1550
1551         if (!pwaff)
1552                 return NULL;
1553
1554         set = isl_set_empty(isl_pw_aff_get_domain_space(pwaff));
1555
1556         for (i = 0; i < pwaff->n; ++i) {
1557                 isl_basic_set *bset;
1558                 isl_set *set_i;
1559
1560                 bset = isl_aff_nonneg_basic_set(isl_aff_copy(pwaff->p[i].aff));
1561                 set_i = isl_set_from_basic_set(bset);
1562                 set_i = isl_set_intersect(set_i, isl_set_copy(pwaff->p[i].set));
1563                 set = isl_set_union_disjoint(set, set_i);
1564         }
1565
1566         isl_pw_aff_free(pwaff);
1567
1568         return set;
1569 }
1570
1571 /* Return a set containing those elements in the domain
1572  * of pwaff where it is zero (if complement is 0) or not zero
1573  * (if complement is 1).
1574  */
1575 static __isl_give isl_set *pw_aff_zero_set(__isl_take isl_pw_aff *pwaff,
1576         int complement)
1577 {
1578         int i;
1579         isl_set *set;
1580
1581         if (!pwaff)
1582                 return NULL;
1583
1584         set = isl_set_empty(isl_pw_aff_get_domain_space(pwaff));
1585
1586         for (i = 0; i < pwaff->n; ++i) {
1587                 isl_basic_set *bset;
1588                 isl_set *set_i, *zero;
1589
1590                 bset = isl_aff_zero_basic_set(isl_aff_copy(pwaff->p[i].aff));
1591                 zero = isl_set_from_basic_set(bset);
1592                 set_i = isl_set_copy(pwaff->p[i].set);
1593                 if (complement)
1594                         set_i = isl_set_subtract(set_i, zero);
1595                 else
1596                         set_i = isl_set_intersect(set_i, zero);
1597                 set = isl_set_union_disjoint(set, set_i);
1598         }
1599
1600         isl_pw_aff_free(pwaff);
1601
1602         return set;
1603 }
1604
1605 /* Return a set containing those elements in the domain
1606  * of pwaff where it is zero.
1607  */
1608 __isl_give isl_set *isl_pw_aff_zero_set(__isl_take isl_pw_aff *pwaff)
1609 {
1610         return pw_aff_zero_set(pwaff, 0);
1611 }
1612
1613 /* Return a set containing those elements in the domain
1614  * of pwaff where it is not zero.
1615  */
1616 __isl_give isl_set *isl_pw_aff_non_zero_set(__isl_take isl_pw_aff *pwaff)
1617 {
1618         return pw_aff_zero_set(pwaff, 1);
1619 }
1620
1621 /* Return a set containing those elements in the shared domain
1622  * of pwaff1 and pwaff2 where pwaff1 is greater than (or equal) to pwaff2.
1623  *
1624  * We compute the difference on the shared domain and then construct
1625  * the set of values where this difference is non-negative.
1626  * If strict is set, we first subtract 1 from the difference.
1627  * If equal is set, we only return the elements where pwaff1 and pwaff2
1628  * are equal.
1629  */
1630 static __isl_give isl_set *pw_aff_gte_set(__isl_take isl_pw_aff *pwaff1,
1631         __isl_take isl_pw_aff *pwaff2, int strict, int equal)
1632 {
1633         isl_set *set1, *set2;
1634
1635         set1 = isl_pw_aff_domain(isl_pw_aff_copy(pwaff1));
1636         set2 = isl_pw_aff_domain(isl_pw_aff_copy(pwaff2));
1637         set1 = isl_set_intersect(set1, set2);
1638         pwaff1 = isl_pw_aff_intersect_domain(pwaff1, isl_set_copy(set1));
1639         pwaff2 = isl_pw_aff_intersect_domain(pwaff2, isl_set_copy(set1));
1640         pwaff1 = isl_pw_aff_add(pwaff1, isl_pw_aff_neg(pwaff2));
1641
1642         if (strict) {
1643                 isl_space *dim = isl_set_get_space(set1);
1644                 isl_aff *aff;
1645                 aff = isl_aff_zero_on_domain(isl_local_space_from_space(dim));
1646                 aff = isl_aff_add_constant_si(aff, -1);
1647                 pwaff1 = isl_pw_aff_add(pwaff1, isl_pw_aff_alloc(set1, aff));
1648         } else
1649                 isl_set_free(set1);
1650
1651         if (equal)
1652                 return isl_pw_aff_zero_set(pwaff1);
1653         return isl_pw_aff_nonneg_set(pwaff1);
1654 }
1655
1656 /* Return a set containing those elements in the shared domain
1657  * of pwaff1 and pwaff2 where pwaff1 is equal to pwaff2.
1658  */
1659 static __isl_give isl_set *pw_aff_eq_set(__isl_take isl_pw_aff *pwaff1,
1660         __isl_take isl_pw_aff *pwaff2)
1661 {
1662         return pw_aff_gte_set(pwaff1, pwaff2, 0, 1);
1663 }
1664
1665 __isl_give isl_set *isl_pw_aff_eq_set(__isl_take isl_pw_aff *pwaff1,
1666         __isl_take isl_pw_aff *pwaff2)
1667 {
1668         return align_params_pw_pw_set_and(pwaff1, pwaff2, &pw_aff_eq_set);
1669 }
1670
1671 /* Return a set containing those elements in the shared domain
1672  * of pwaff1 and pwaff2 where pwaff1 is greater than or equal to pwaff2.
1673  */
1674 static __isl_give isl_set *pw_aff_ge_set(__isl_take isl_pw_aff *pwaff1,
1675         __isl_take isl_pw_aff *pwaff2)
1676 {
1677         return pw_aff_gte_set(pwaff1, pwaff2, 0, 0);
1678 }
1679
1680 __isl_give isl_set *isl_pw_aff_ge_set(__isl_take isl_pw_aff *pwaff1,
1681         __isl_take isl_pw_aff *pwaff2)
1682 {
1683         return align_params_pw_pw_set_and(pwaff1, pwaff2, &pw_aff_ge_set);
1684 }
1685
1686 /* Return a set containing those elements in the shared domain
1687  * of pwaff1 and pwaff2 where pwaff1 is strictly greater than pwaff2.
1688  */
1689 static __isl_give isl_set *pw_aff_gt_set(__isl_take isl_pw_aff *pwaff1,
1690         __isl_take isl_pw_aff *pwaff2)
1691 {
1692         return pw_aff_gte_set(pwaff1, pwaff2, 1, 0);
1693 }
1694
1695 __isl_give isl_set *isl_pw_aff_gt_set(__isl_take isl_pw_aff *pwaff1,
1696         __isl_take isl_pw_aff *pwaff2)
1697 {
1698         return align_params_pw_pw_set_and(pwaff1, pwaff2, &pw_aff_gt_set);
1699 }
1700
1701 __isl_give isl_set *isl_pw_aff_le_set(__isl_take isl_pw_aff *pwaff1,
1702         __isl_take isl_pw_aff *pwaff2)
1703 {
1704         return isl_pw_aff_ge_set(pwaff2, pwaff1);
1705 }
1706
1707 __isl_give isl_set *isl_pw_aff_lt_set(__isl_take isl_pw_aff *pwaff1,
1708         __isl_take isl_pw_aff *pwaff2)
1709 {
1710         return isl_pw_aff_gt_set(pwaff2, pwaff1);
1711 }
1712
1713 /* Return a set containing those elements in the shared domain
1714  * of the elements of list1 and list2 where each element in list1
1715  * has the relation specified by "fn" with each element in list2.
1716  */
1717 static __isl_give isl_set *pw_aff_list_set(__isl_take isl_pw_aff_list *list1,
1718         __isl_take isl_pw_aff_list *list2,
1719         __isl_give isl_set *(*fn)(__isl_take isl_pw_aff *pwaff1,
1720                                     __isl_take isl_pw_aff *pwaff2))
1721 {
1722         int i, j;
1723         isl_ctx *ctx;
1724         isl_set *set;
1725
1726         if (!list1 || !list2)
1727                 goto error;
1728
1729         ctx = isl_pw_aff_list_get_ctx(list1);
1730         if (list1->n < 1 || list2->n < 1)
1731                 isl_die(ctx, isl_error_invalid,
1732                         "list should contain at least one element", goto error);
1733
1734         set = isl_set_universe(isl_pw_aff_get_domain_space(list1->p[0]));
1735         for (i = 0; i < list1->n; ++i)
1736                 for (j = 0; j < list2->n; ++j) {
1737                         isl_set *set_ij;
1738
1739                         set_ij = fn(isl_pw_aff_copy(list1->p[i]),
1740                                     isl_pw_aff_copy(list2->p[j]));
1741                         set = isl_set_intersect(set, set_ij);
1742                 }
1743
1744         isl_pw_aff_list_free(list1);
1745         isl_pw_aff_list_free(list2);
1746         return set;
1747 error:
1748         isl_pw_aff_list_free(list1);
1749         isl_pw_aff_list_free(list2);
1750         return NULL;
1751 }
1752
1753 /* Return a set containing those elements in the shared domain
1754  * of the elements of list1 and list2 where each element in list1
1755  * is equal to each element in list2.
1756  */
1757 __isl_give isl_set *isl_pw_aff_list_eq_set(__isl_take isl_pw_aff_list *list1,
1758         __isl_take isl_pw_aff_list *list2)
1759 {
1760         return pw_aff_list_set(list1, list2, &isl_pw_aff_eq_set);
1761 }
1762
1763 __isl_give isl_set *isl_pw_aff_list_ne_set(__isl_take isl_pw_aff_list *list1,
1764         __isl_take isl_pw_aff_list *list2)
1765 {
1766         return pw_aff_list_set(list1, list2, &isl_pw_aff_ne_set);
1767 }
1768
1769 /* Return a set containing those elements in the shared domain
1770  * of the elements of list1 and list2 where each element in list1
1771  * is less than or equal to each element in list2.
1772  */
1773 __isl_give isl_set *isl_pw_aff_list_le_set(__isl_take isl_pw_aff_list *list1,
1774         __isl_take isl_pw_aff_list *list2)
1775 {
1776         return pw_aff_list_set(list1, list2, &isl_pw_aff_le_set);
1777 }
1778
1779 __isl_give isl_set *isl_pw_aff_list_lt_set(__isl_take isl_pw_aff_list *list1,
1780         __isl_take isl_pw_aff_list *list2)
1781 {
1782         return pw_aff_list_set(list1, list2, &isl_pw_aff_lt_set);
1783 }
1784
1785 __isl_give isl_set *isl_pw_aff_list_ge_set(__isl_take isl_pw_aff_list *list1,
1786         __isl_take isl_pw_aff_list *list2)
1787 {
1788         return pw_aff_list_set(list1, list2, &isl_pw_aff_ge_set);
1789 }
1790
1791 __isl_give isl_set *isl_pw_aff_list_gt_set(__isl_take isl_pw_aff_list *list1,
1792         __isl_take isl_pw_aff_list *list2)
1793 {
1794         return pw_aff_list_set(list1, list2, &isl_pw_aff_gt_set);
1795 }
1796
1797
1798 /* Return a set containing those elements in the shared domain
1799  * of pwaff1 and pwaff2 where pwaff1 is not equal to pwaff2.
1800  */
1801 static __isl_give isl_set *pw_aff_ne_set(__isl_take isl_pw_aff *pwaff1,
1802         __isl_take isl_pw_aff *pwaff2)
1803 {
1804         isl_set *set_lt, *set_gt;
1805
1806         set_lt = isl_pw_aff_lt_set(isl_pw_aff_copy(pwaff1),
1807                                    isl_pw_aff_copy(pwaff2));
1808         set_gt = isl_pw_aff_gt_set(pwaff1, pwaff2);
1809         return isl_set_union_disjoint(set_lt, set_gt);
1810 }
1811
1812 __isl_give isl_set *isl_pw_aff_ne_set(__isl_take isl_pw_aff *pwaff1,
1813         __isl_take isl_pw_aff *pwaff2)
1814 {
1815         return align_params_pw_pw_set_and(pwaff1, pwaff2, &pw_aff_ne_set);
1816 }
1817
1818 __isl_give isl_pw_aff *isl_pw_aff_scale_down(__isl_take isl_pw_aff *pwaff,
1819         isl_int v)
1820 {
1821         int i;
1822
1823         if (isl_int_is_one(v))
1824                 return pwaff;
1825         if (!isl_int_is_pos(v))
1826                 isl_die(isl_pw_aff_get_ctx(pwaff), isl_error_invalid,
1827                         "factor needs to be positive",
1828                         return isl_pw_aff_free(pwaff));
1829         pwaff = isl_pw_aff_cow(pwaff);
1830         if (!pwaff)
1831                 return NULL;
1832         if (pwaff->n == 0)
1833                 return pwaff;
1834
1835         for (i = 0; i < pwaff->n; ++i) {
1836                 pwaff->p[i].aff = isl_aff_scale_down(pwaff->p[i].aff, v);
1837                 if (!pwaff->p[i].aff)
1838                         return isl_pw_aff_free(pwaff);
1839         }
1840
1841         return pwaff;
1842 }
1843
1844 __isl_give isl_pw_aff *isl_pw_aff_floor(__isl_take isl_pw_aff *pwaff)
1845 {
1846         int i;
1847
1848         pwaff = isl_pw_aff_cow(pwaff);
1849         if (!pwaff)
1850                 return NULL;
1851         if (pwaff->n == 0)
1852                 return pwaff;
1853
1854         for (i = 0; i < pwaff->n; ++i) {
1855                 pwaff->p[i].aff = isl_aff_floor(pwaff->p[i].aff);
1856                 if (!pwaff->p[i].aff)
1857                         return isl_pw_aff_free(pwaff);
1858         }
1859
1860         return pwaff;
1861 }
1862
1863 __isl_give isl_pw_aff *isl_pw_aff_ceil(__isl_take isl_pw_aff *pwaff)
1864 {
1865         int i;
1866
1867         pwaff = isl_pw_aff_cow(pwaff);
1868         if (!pwaff)
1869                 return NULL;
1870         if (pwaff->n == 0)
1871                 return pwaff;
1872
1873         for (i = 0; i < pwaff->n; ++i) {
1874                 pwaff->p[i].aff = isl_aff_ceil(pwaff->p[i].aff);
1875                 if (!pwaff->p[i].aff)
1876                         return isl_pw_aff_free(pwaff);
1877         }
1878
1879         return pwaff;
1880 }
1881
1882 /* Return an affine expression that is equal to pwaff_true for elements
1883  * in "cond" and to pwaff_false for elements not in "cond".
1884  * That is, return cond ? pwaff_true : pwaff_false;
1885  */
1886 __isl_give isl_pw_aff *isl_pw_aff_cond(__isl_take isl_set *cond,
1887         __isl_take isl_pw_aff *pwaff_true, __isl_take isl_pw_aff *pwaff_false)
1888 {
1889         isl_set *comp;
1890
1891         comp = isl_set_complement(isl_set_copy(cond));
1892         pwaff_true = isl_pw_aff_intersect_domain(pwaff_true, cond);
1893         pwaff_false = isl_pw_aff_intersect_domain(pwaff_false, comp);
1894
1895         return isl_pw_aff_add_disjoint(pwaff_true, pwaff_false);
1896 }
1897
1898 int isl_aff_is_cst(__isl_keep isl_aff *aff)
1899 {
1900         if (!aff)
1901                 return -1;
1902
1903         return isl_seq_first_non_zero(aff->v->el + 2, aff->v->size - 2) == -1;
1904 }
1905
1906 /* Check whether pwaff is a piecewise constant.
1907  */
1908 int isl_pw_aff_is_cst(__isl_keep isl_pw_aff *pwaff)
1909 {
1910         int i;
1911
1912         if (!pwaff)
1913                 return -1;
1914
1915         for (i = 0; i < pwaff->n; ++i) {
1916                 int is_cst = isl_aff_is_cst(pwaff->p[i].aff);
1917                 if (is_cst < 0 || !is_cst)
1918                         return is_cst;
1919         }
1920
1921         return 1;
1922 }
1923
1924 __isl_give isl_aff *isl_aff_mul(__isl_take isl_aff *aff1,
1925         __isl_take isl_aff *aff2)
1926 {
1927         if (!isl_aff_is_cst(aff2) && isl_aff_is_cst(aff1))
1928                 return isl_aff_mul(aff2, aff1);
1929
1930         if (!isl_aff_is_cst(aff2))
1931                 isl_die(isl_aff_get_ctx(aff1), isl_error_invalid,
1932                         "at least one affine expression should be constant",
1933                         goto error);
1934
1935         aff1 = isl_aff_cow(aff1);
1936         if (!aff1 || !aff2)
1937                 goto error;
1938
1939         aff1 = isl_aff_scale(aff1, aff2->v->el[1]);
1940         aff1 = isl_aff_scale_down(aff1, aff2->v->el[0]);
1941
1942         isl_aff_free(aff2);
1943         return aff1;
1944 error:
1945         isl_aff_free(aff1);
1946         isl_aff_free(aff2);
1947         return NULL;
1948 }
1949
1950 static __isl_give isl_pw_aff *pw_aff_add(__isl_take isl_pw_aff *pwaff1,
1951         __isl_take isl_pw_aff *pwaff2)
1952 {
1953         return isl_pw_aff_on_shared_domain(pwaff1, pwaff2, &isl_aff_add);
1954 }
1955
1956 __isl_give isl_pw_aff *isl_pw_aff_add(__isl_take isl_pw_aff *pwaff1,
1957         __isl_take isl_pw_aff *pwaff2)
1958 {
1959         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2, &pw_aff_add);
1960 }
1961
1962 __isl_give isl_pw_aff *isl_pw_aff_union_add(__isl_take isl_pw_aff *pwaff1,
1963         __isl_take isl_pw_aff *pwaff2)
1964 {
1965         return isl_pw_aff_union_add_(pwaff1, pwaff2);
1966 }
1967
1968 static __isl_give isl_pw_aff *pw_aff_mul(__isl_take isl_pw_aff *pwaff1,
1969         __isl_take isl_pw_aff *pwaff2)
1970 {
1971         return isl_pw_aff_on_shared_domain(pwaff1, pwaff2, &isl_aff_mul);
1972 }
1973
1974 __isl_give isl_pw_aff *isl_pw_aff_mul(__isl_take isl_pw_aff *pwaff1,
1975         __isl_take isl_pw_aff *pwaff2)
1976 {
1977         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2, &pw_aff_mul);
1978 }
1979
1980 static __isl_give isl_pw_aff *pw_aff_min(__isl_take isl_pw_aff *pwaff1,
1981         __isl_take isl_pw_aff *pwaff2)
1982 {
1983         isl_set *le;
1984
1985         le = isl_pw_aff_le_set(isl_pw_aff_copy(pwaff1),
1986                                 isl_pw_aff_copy(pwaff2));
1987         return isl_pw_aff_cond(le, pwaff1, pwaff2);
1988 }
1989
1990 __isl_give isl_pw_aff *isl_pw_aff_min(__isl_take isl_pw_aff *pwaff1,
1991         __isl_take isl_pw_aff *pwaff2)
1992 {
1993         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2, &pw_aff_min);
1994 }
1995
1996 static __isl_give isl_pw_aff *pw_aff_max(__isl_take isl_pw_aff *pwaff1,
1997         __isl_take isl_pw_aff *pwaff2)
1998 {
1999         isl_set *le;
2000
2001         le = isl_pw_aff_ge_set(isl_pw_aff_copy(pwaff1),
2002                                 isl_pw_aff_copy(pwaff2));
2003         return isl_pw_aff_cond(le, pwaff1, pwaff2);
2004 }
2005
2006 __isl_give isl_pw_aff *isl_pw_aff_max(__isl_take isl_pw_aff *pwaff1,
2007         __isl_take isl_pw_aff *pwaff2)
2008 {
2009         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2, &pw_aff_max);
2010 }
2011
2012 static __isl_give isl_pw_aff *pw_aff_list_reduce(
2013         __isl_take isl_pw_aff_list *list,
2014         __isl_give isl_pw_aff *(*fn)(__isl_take isl_pw_aff *pwaff1,
2015                                         __isl_take isl_pw_aff *pwaff2))
2016 {
2017         int i;
2018         isl_ctx *ctx;
2019         isl_pw_aff *res;
2020
2021         if (!list)
2022                 return NULL;
2023
2024         ctx = isl_pw_aff_list_get_ctx(list);
2025         if (list->n < 1)
2026                 isl_die(ctx, isl_error_invalid,
2027                         "list should contain at least one element",
2028                         return isl_pw_aff_list_free(list));
2029
2030         res = isl_pw_aff_copy(list->p[0]);
2031         for (i = 1; i < list->n; ++i)
2032                 res = fn(res, isl_pw_aff_copy(list->p[i]));
2033
2034         isl_pw_aff_list_free(list);
2035         return res;
2036 }
2037
2038 /* Return an isl_pw_aff that maps each element in the intersection of the
2039  * domains of the elements of list to the minimal corresponding affine
2040  * expression.
2041  */
2042 __isl_give isl_pw_aff *isl_pw_aff_list_min(__isl_take isl_pw_aff_list *list)
2043 {
2044         return pw_aff_list_reduce(list, &isl_pw_aff_min);
2045 }
2046
2047 /* Return an isl_pw_aff that maps each element in the intersection of the
2048  * domains of the elements of list to the maximal corresponding affine
2049  * expression.
2050  */
2051 __isl_give isl_pw_aff *isl_pw_aff_list_max(__isl_take isl_pw_aff_list *list)
2052 {
2053         return pw_aff_list_reduce(list, &isl_pw_aff_max);
2054 }
2055
2056 #undef BASE
2057 #define BASE aff
2058
2059 #include <isl_multi_templ.c>
2060
2061 __isl_give isl_multi_aff *isl_multi_aff_add(__isl_take isl_multi_aff *maff1,
2062         __isl_take isl_multi_aff *maff2)
2063 {
2064         int i;
2065         isl_ctx *ctx;
2066
2067         maff1 = isl_multi_aff_cow(maff1);
2068         if (!maff1 || !maff2)
2069                 goto error;
2070
2071         ctx = isl_multi_aff_get_ctx(maff1);
2072         if (!isl_space_is_equal(maff1->space, maff2->space))
2073                 isl_die(ctx, isl_error_invalid,
2074                         "spaces don't match", goto error);
2075
2076         for (i = 0; i < maff1->n; ++i) {
2077                 maff1->p[i] = isl_aff_add(maff1->p[i],
2078                                             isl_aff_copy(maff2->p[i]));
2079                 if (!maff1->p[i])
2080                         goto error;
2081         }
2082
2083         isl_multi_aff_free(maff2);
2084         return maff1;
2085 error:
2086         isl_multi_aff_free(maff1);
2087         isl_multi_aff_free(maff2);
2088         return NULL;
2089 }
2090
2091 /* Exploit the equalities in "eq" to simplify the affine expressions.
2092  */
2093 static __isl_give isl_multi_aff *isl_multi_aff_substitute_equalities(
2094         __isl_take isl_multi_aff *maff, __isl_take isl_basic_set *eq)
2095 {
2096         int i;
2097
2098         maff = isl_multi_aff_cow(maff);
2099         if (!maff || !eq)
2100                 goto error;
2101
2102         for (i = 0; i < maff->n; ++i) {
2103                 maff->p[i] = isl_aff_substitute_equalities(maff->p[i],
2104                                                     isl_basic_set_copy(eq));
2105                 if (!maff->p[i])
2106                         goto error;
2107         }
2108
2109         isl_basic_set_free(eq);
2110         return maff;
2111 error:
2112         isl_basic_set_free(eq);
2113         isl_multi_aff_free(maff);
2114         return NULL;
2115 }
2116
2117 __isl_give isl_multi_aff *isl_multi_aff_scale(__isl_take isl_multi_aff *maff,
2118         isl_int f)
2119 {
2120         int i;
2121
2122         maff = isl_multi_aff_cow(maff);
2123         if (!maff)
2124                 return NULL;
2125
2126         for (i = 0; i < maff->n; ++i) {
2127                 maff->p[i] = isl_aff_scale(maff->p[i], f);
2128                 if (!maff->p[i])
2129                         return isl_multi_aff_free(maff);
2130         }
2131
2132         return maff;
2133 }
2134
2135 __isl_give isl_multi_aff *isl_multi_aff_add_on_domain(__isl_keep isl_set *dom,
2136         __isl_take isl_multi_aff *maff1, __isl_take isl_multi_aff *maff2)
2137 {
2138         maff1 = isl_multi_aff_add(maff1, maff2);
2139         maff1 = isl_multi_aff_gist(maff1, isl_set_copy(dom));
2140         return maff1;
2141 }
2142
2143 int isl_multi_aff_is_empty(__isl_keep isl_multi_aff *maff)
2144 {
2145         if (!maff)
2146                 return -1;
2147
2148         return 0;
2149 }
2150
2151 int isl_multi_aff_plain_is_equal(__isl_keep isl_multi_aff *maff1,
2152         __isl_keep isl_multi_aff *maff2)
2153 {
2154         int i;
2155         int equal;
2156
2157         if (!maff1 || !maff2)
2158                 return -1;
2159         if (maff1->n != maff2->n)
2160                 return 0;
2161         equal = isl_space_is_equal(maff1->space, maff2->space);
2162         if (equal < 0 || !equal)
2163                 return equal;
2164
2165         for (i = 0; i < maff1->n; ++i) {
2166                 equal = isl_aff_plain_is_equal(maff1->p[i], maff2->p[i]);
2167                 if (equal < 0 || !equal)
2168                         return equal;
2169         }
2170
2171         return 1;
2172 }
2173
2174 __isl_give isl_multi_aff *isl_multi_aff_set_dim_name(
2175         __isl_take isl_multi_aff *maff,
2176         enum isl_dim_type type, unsigned pos, const char *s)
2177 {
2178         int i;
2179
2180         maff = isl_multi_aff_cow(maff);
2181         if (!maff)
2182                 return NULL;
2183
2184         maff->space = isl_space_set_dim_name(maff->space, type, pos, s);
2185         if (!maff->space)
2186                 return isl_multi_aff_free(maff);
2187         for (i = 0; i < maff->n; ++i) {
2188                 maff->p[i] = isl_aff_set_dim_name(maff->p[i], type, pos, s);
2189                 if (!maff->p[i])
2190                         return isl_multi_aff_free(maff);
2191         }
2192
2193         return maff;
2194 }
2195
2196 __isl_give isl_multi_aff *isl_multi_aff_drop_dims(__isl_take isl_multi_aff *maff,
2197         enum isl_dim_type type, unsigned first, unsigned n)
2198 {
2199         int i;
2200
2201         maff = isl_multi_aff_cow(maff);
2202         if (!maff)
2203                 return NULL;
2204
2205         maff->space = isl_space_drop_dims(maff->space, type, first, n);
2206         if (!maff->space)
2207                 return isl_multi_aff_free(maff);
2208
2209         if (type == isl_dim_out) {
2210                 for (i = 0; i < n; ++i)
2211                         isl_aff_free(maff->p[first + i]);
2212                 for (i = first; i + n < maff->n; ++i)
2213                         maff->p[i] = maff->p[i + n];
2214                 maff->n -= n;
2215                 return maff;
2216         }
2217
2218         for (i = 0; i < maff->n; ++i) {
2219                 maff->p[i] = isl_aff_drop_dims(maff->p[i], type, first, n);
2220                 if (!maff->p[i])
2221                         return isl_multi_aff_free(maff);
2222         }
2223
2224         return maff;
2225 }
2226
2227 #undef PW
2228 #define PW isl_pw_multi_aff
2229 #undef EL
2230 #define EL isl_multi_aff
2231 #undef EL_IS_ZERO
2232 #define EL_IS_ZERO is_empty
2233 #undef ZERO
2234 #define ZERO empty
2235 #undef IS_ZERO
2236 #define IS_ZERO is_empty
2237 #undef FIELD
2238 #define FIELD maff
2239 #undef DEFAULT_IS_ZERO
2240 #define DEFAULT_IS_ZERO 0
2241
2242 #define NO_NEG
2243 #define NO_EVAL
2244 #define NO_OPT
2245 #define NO_INVOLVES_DIMS
2246 #define NO_MOVE_DIMS
2247 #define NO_INSERT_DIMS
2248 #define NO_LIFT
2249 #define NO_MORPH
2250
2251 #include <isl_pw_templ.c>
2252
2253 static __isl_give isl_pw_multi_aff *pw_multi_aff_add(
2254         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
2255 {
2256         return isl_pw_multi_aff_on_shared_domain(pma1, pma2,
2257                                                 &isl_multi_aff_add);
2258 }
2259
2260 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_add(
2261         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
2262 {
2263         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
2264                                                 &pw_multi_aff_add);
2265 }
2266
2267 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_union_add(
2268         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
2269 {
2270         return isl_pw_multi_aff_union_add_(pma1, pma2);
2271 }
2272
2273 /* Construct a map mapping the domain the piecewise multi-affine expression
2274  * to its range, with each dimension in the range equated to the
2275  * corresponding affine expression on its cell.
2276  */
2277 __isl_give isl_map *isl_map_from_pw_multi_aff(__isl_take isl_pw_multi_aff *pma)
2278 {
2279         int i;
2280         isl_map *map;
2281
2282         if (!pma)
2283                 return NULL;
2284
2285         map = isl_map_empty(isl_pw_multi_aff_get_space(pma));
2286
2287         for (i = 0; i < pma->n; ++i) {
2288                 isl_multi_aff *maff;
2289                 isl_basic_map *bmap;
2290                 isl_map *map_i;
2291
2292                 maff = isl_multi_aff_copy(pma->p[i].maff);
2293                 bmap = isl_basic_map_from_multi_aff(maff);
2294                 map_i = isl_map_from_basic_map(bmap);
2295                 map_i = isl_map_intersect_domain(map_i,
2296                                                 isl_set_copy(pma->p[i].set));
2297                 map = isl_map_union_disjoint(map, map_i);
2298         }
2299
2300         isl_pw_multi_aff_free(pma);
2301         return map;
2302 }
2303
2304 __isl_give isl_set *isl_set_from_pw_multi_aff(__isl_take isl_pw_multi_aff *pma)
2305 {
2306         if (!isl_space_is_set(pma->dim))
2307                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
2308                         "isl_pw_multi_aff cannot be converted into an isl_set",
2309                         return isl_pw_multi_aff_free(pma));
2310
2311         return isl_map_from_pw_multi_aff(pma);
2312 }
2313
2314 /* Given a basic map with a single output dimension that is defined
2315  * in terms of the parameters and input dimensions using an equality,
2316  * extract an isl_aff that expresses the output dimension in terms
2317  * of the parameters and input dimensions.
2318  *
2319  * Since some applications expect the result of isl_pw_multi_aff_from_map
2320  * to only contain integer affine expressions, we compute the floor
2321  * of the expression before returning.
2322  *
2323  * This function shares some similarities with
2324  * isl_basic_map_has_defining_equality and isl_constraint_get_bound.
2325  */
2326 static __isl_give isl_aff *extract_isl_aff_from_basic_map(
2327         __isl_take isl_basic_map *bmap)
2328 {
2329         int i;
2330         unsigned offset;
2331         unsigned total;
2332         isl_local_space *ls;
2333         isl_aff *aff;
2334
2335         if (!bmap)
2336                 return NULL;
2337         if (isl_basic_map_dim(bmap, isl_dim_out) != 1)
2338                 isl_die(isl_basic_map_get_ctx(bmap), isl_error_invalid,
2339                         "basic map should have a single output dimension",
2340                         goto error);
2341         offset = isl_basic_map_offset(bmap, isl_dim_out);
2342         total = isl_basic_map_total_dim(bmap);
2343         for (i = 0; i < bmap->n_eq; ++i) {
2344                 if (isl_int_is_zero(bmap->eq[i][offset]))
2345                         continue;
2346                 if (isl_seq_first_non_zero(bmap->eq[i] + offset + 1,
2347                                            1 + total - (offset + 1)) != -1)
2348                         continue;
2349                 break;
2350         }
2351         if (i >= bmap->n_eq)
2352                 isl_die(isl_basic_map_get_ctx(bmap), isl_error_invalid,
2353                         "unable to find suitable equality", goto error);
2354         ls = isl_basic_map_get_local_space(bmap);
2355         aff = isl_aff_alloc(isl_local_space_domain(ls));
2356         if (!aff)
2357                 goto error;
2358         if (isl_int_is_neg(bmap->eq[i][offset]))
2359                 isl_seq_cpy(aff->v->el + 1, bmap->eq[i], offset);
2360         else
2361                 isl_seq_neg(aff->v->el + 1, bmap->eq[i], offset);
2362         isl_seq_clr(aff->v->el + 1 + offset, aff->v->size - (1 + offset));
2363         isl_int_abs(aff->v->el[0], bmap->eq[i][offset]);
2364         isl_basic_map_free(bmap);
2365
2366         aff = isl_aff_remove_unused_divs(aff);
2367         aff = isl_aff_floor(aff);
2368         return aff;
2369 error:
2370         isl_basic_map_free(bmap);
2371         return NULL;
2372 }
2373
2374 /* Given a basic map where each output dimension is defined
2375  * in terms of the parameters and input dimensions using an equality,
2376  * extract an isl_multi_aff that expresses the output dimensions in terms
2377  * of the parameters and input dimensions.
2378  */
2379 static __isl_give isl_multi_aff *extract_isl_multi_aff_from_basic_map(
2380         __isl_take isl_basic_map *bmap)
2381 {
2382         int i;
2383         unsigned n_out;
2384         isl_multi_aff *ma;
2385
2386         if (!bmap)
2387                 return NULL;
2388
2389         ma = isl_multi_aff_alloc(isl_basic_map_get_space(bmap));
2390         n_out = isl_basic_map_dim(bmap, isl_dim_out);
2391
2392         for (i = 0; i < n_out; ++i) {
2393                 isl_basic_map *bmap_i;
2394                 isl_aff *aff;
2395
2396                 bmap_i = isl_basic_map_copy(bmap);
2397                 bmap_i = isl_basic_map_project_out(bmap_i, isl_dim_out,
2398                                                         i + 1, n_out - (1 + i));
2399                 bmap_i = isl_basic_map_project_out(bmap_i, isl_dim_out, 0, i);
2400                 aff = extract_isl_aff_from_basic_map(bmap_i);
2401                 ma = isl_multi_aff_set_aff(ma, i, aff);
2402         }
2403
2404         isl_basic_map_free(bmap);
2405
2406         return ma;
2407 }
2408
2409 /* Create an isl_pw_multi_aff that is equivalent to
2410  * isl_map_intersect_domain(isl_map_from_basic_map(bmap), domain).
2411  * The given basic map is such that each output dimension is defined
2412  * in terms of the parameters and input dimensions using an equality.
2413  */
2414 static __isl_give isl_pw_multi_aff *plain_pw_multi_aff_from_map(
2415         __isl_take isl_set *domain, __isl_take isl_basic_map *bmap)
2416 {
2417         isl_multi_aff *ma;
2418
2419         ma = extract_isl_multi_aff_from_basic_map(bmap);
2420         return isl_pw_multi_aff_alloc(domain, ma);
2421 }
2422
2423 /* Try and create an isl_pw_multi_aff that is equivalent to the given isl_map.
2424  * This obivously only works if the input "map" is single-valued.
2425  * If so, we compute the lexicographic minimum of the image in the form
2426  * of an isl_pw_multi_aff.  Since the image is unique, it is equal
2427  * to its lexicographic minimum.
2428  * If the input is not single-valued, we produce an error.
2429  *
2430  * As a special case, we first check if all output dimensions are uniquely
2431  * defined in terms of the parameters and input dimensions over the entire
2432  * domain.  If so, we extract the desired isl_pw_multi_aff directly
2433  * from the affine hull of "map" and its domain.
2434  */
2435 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_from_map(__isl_take isl_map *map)
2436 {
2437         int i;
2438         int sv;
2439         isl_pw_multi_aff *pma;
2440         isl_basic_map *hull;
2441
2442         if (!map)
2443                 return NULL;
2444
2445         hull = isl_map_affine_hull(isl_map_copy(map));
2446         sv = isl_basic_map_plain_is_single_valued(hull);
2447         if (sv >= 0 && sv)
2448                 return plain_pw_multi_aff_from_map(isl_map_domain(map), hull);
2449         isl_basic_map_free(hull);
2450         if (sv < 0)
2451                 goto error;
2452
2453         sv = isl_map_is_single_valued(map);
2454         if (sv < 0)
2455                 goto error;
2456         if (!sv)
2457                 isl_die(isl_map_get_ctx(map), isl_error_invalid,
2458                         "map is not single-valued", goto error);
2459         map = isl_map_make_disjoint(map);
2460         if (!map)
2461                 return NULL;
2462
2463         pma = isl_pw_multi_aff_empty(isl_map_get_space(map));
2464
2465         for (i = 0; i < map->n; ++i) {
2466                 isl_pw_multi_aff *pma_i;
2467                 isl_basic_map *bmap;
2468                 bmap = isl_basic_map_copy(map->p[i]);
2469                 pma_i = isl_basic_map_lexmin_pw_multi_aff(bmap);
2470                 pma = isl_pw_multi_aff_add_disjoint(pma, pma_i);
2471         }
2472
2473         isl_map_free(map);
2474         return pma;
2475 error:
2476         isl_map_free(map);
2477         return NULL;
2478 }
2479
2480 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_from_set(__isl_take isl_set *set)
2481 {
2482         return isl_pw_multi_aff_from_map(set);
2483 }
2484
2485 /* Return the piecewise affine expression "set ? 1 : 0".
2486  */
2487 __isl_give isl_pw_aff *isl_set_indicator_function(__isl_take isl_set *set)
2488 {
2489         isl_pw_aff *pa;
2490         isl_space *space = isl_set_get_space(set);
2491         isl_local_space *ls = isl_local_space_from_space(space);
2492         isl_aff *zero = isl_aff_zero_on_domain(isl_local_space_copy(ls));
2493         isl_aff *one = isl_aff_zero_on_domain(ls);
2494
2495         one = isl_aff_add_constant_si(one, 1);
2496         pa = isl_pw_aff_alloc(isl_set_copy(set), one);
2497         set = isl_set_complement(set);
2498         pa = isl_pw_aff_add_disjoint(pa, isl_pw_aff_alloc(set, zero));
2499
2500         return pa;
2501 }
2502
2503 /* Plug in "subs" for dimension "type", "pos" of "aff".
2504  *
2505  * Let i be the dimension to replace and let "subs" be of the form
2506  *
2507  *      f/d
2508  *
2509  * and "aff" of the form
2510  *
2511  *      (a i + g)/m
2512  *
2513  * The result is
2514  *
2515  *      floor((a f + d g')/(m d))
2516  *
2517  * where g' is the result of plugging in "subs" in each of the integer
2518  * divisions in g.
2519  */
2520 __isl_give isl_aff *isl_aff_substitute(__isl_take isl_aff *aff,
2521         enum isl_dim_type type, unsigned pos, __isl_keep isl_aff *subs)
2522 {
2523         isl_ctx *ctx;
2524         isl_int v;
2525
2526         aff = isl_aff_cow(aff);
2527         if (!aff || !subs)
2528                 return isl_aff_free(aff);
2529
2530         ctx = isl_aff_get_ctx(aff);
2531         if (!isl_space_is_equal(aff->ls->dim, subs->ls->dim))
2532                 isl_die(ctx, isl_error_invalid,
2533                         "spaces don't match", return isl_aff_free(aff));
2534         if (isl_local_space_dim(subs->ls, isl_dim_div) != 0)
2535                 isl_die(ctx, isl_error_unsupported,
2536                         "cannot handle divs yet", return isl_aff_free(aff));
2537
2538         aff->ls = isl_local_space_substitute(aff->ls, type, pos, subs);
2539         if (!aff->ls)
2540                 return isl_aff_free(aff);
2541
2542         aff->v = isl_vec_cow(aff->v);
2543         if (!aff->v)
2544                 return isl_aff_free(aff);
2545
2546         pos += isl_local_space_offset(aff->ls, type);
2547
2548         isl_int_init(v);
2549         isl_int_set(v, aff->v->el[1 + pos]);
2550         isl_int_set_si(aff->v->el[1 + pos], 0);
2551         isl_seq_combine(aff->v->el + 1, subs->v->el[0], aff->v->el + 1,
2552                         v, subs->v->el + 1, subs->v->size - 1);
2553         isl_int_mul(aff->v->el[0], aff->v->el[0], subs->v->el[0]);
2554         isl_int_clear(v);
2555
2556         return aff;
2557 }
2558
2559 /* Plug in "subs" for dimension "type", "pos" in each of the affine
2560  * expressions in "maff".
2561  */
2562 __isl_give isl_multi_aff *isl_multi_aff_substitute(
2563         __isl_take isl_multi_aff *maff, enum isl_dim_type type, unsigned pos,
2564         __isl_keep isl_aff *subs)
2565 {
2566         int i;
2567
2568         maff = isl_multi_aff_cow(maff);
2569         if (!maff || !subs)
2570                 return isl_multi_aff_free(maff);
2571
2572         if (type == isl_dim_in)
2573                 type = isl_dim_set;
2574
2575         for (i = 0; i < maff->n; ++i) {
2576                 maff->p[i] = isl_aff_substitute(maff->p[i], type, pos, subs);
2577                 if (!maff->p[i])
2578                         return isl_multi_aff_free(maff);
2579         }
2580
2581         return maff;
2582 }
2583
2584 /* Plug in "subs" for dimension "type", "pos" of "pma".
2585  *
2586  * pma is of the form
2587  *
2588  *      A_i(v) -> M_i(v)
2589  *
2590  * while subs is of the form
2591  *
2592  *      v' = B_j(v) -> S_j
2593  *
2594  * Each pair i,j such that C_ij = A_i \cap B_i is non-empty
2595  * has a contribution in the result, in particular
2596  *
2597  *      C_ij(S_j) -> M_i(S_j)
2598  *
2599  * Note that plugging in S_j in C_ij may also result in an empty set
2600  * and this contribution should simply be discarded.
2601  */
2602 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_substitute(
2603         __isl_take isl_pw_multi_aff *pma, enum isl_dim_type type, unsigned pos,
2604         __isl_keep isl_pw_aff *subs)
2605 {
2606         int i, j, n;
2607         isl_pw_multi_aff *res;
2608
2609         if (!pma || !subs)
2610                 return isl_pw_multi_aff_free(pma);
2611
2612         n = pma->n * subs->n;
2613         res = isl_pw_multi_aff_alloc_size(isl_space_copy(pma->dim), n);
2614
2615         for (i = 0; i < pma->n; ++i) {
2616                 for (j = 0; j < subs->n; ++j) {
2617                         isl_set *common;
2618                         isl_multi_aff *res_ij;
2619                         common = isl_set_intersect(
2620                                         isl_set_copy(pma->p[i].set),
2621                                         isl_set_copy(subs->p[j].set));
2622                         common = isl_set_substitute(common,
2623                                         type, pos, subs->p[j].aff);
2624                         if (isl_set_plain_is_empty(common)) {
2625                                 isl_set_free(common);
2626                                 continue;
2627                         }
2628
2629                         res_ij = isl_multi_aff_substitute(
2630                                         isl_multi_aff_copy(pma->p[i].maff),
2631                                         type, pos, subs->p[j].aff);
2632
2633                         res = isl_pw_multi_aff_add_piece(res, common, res_ij);
2634                 }
2635         }
2636
2637         isl_pw_multi_aff_free(pma);
2638         return res;
2639 }
2640
2641 /* Extend the local space of "dst" to include the divs
2642  * in the local space of "src".
2643  */
2644 __isl_give isl_aff *isl_aff_align_divs(__isl_take isl_aff *dst,
2645         __isl_keep isl_aff *src)
2646 {
2647         isl_ctx *ctx;
2648         int *exp1 = NULL;
2649         int *exp2 = NULL;
2650         isl_mat *div;
2651
2652         if (!src || !dst)
2653                 return isl_aff_free(dst);
2654
2655         ctx = isl_aff_get_ctx(src);
2656         if (!isl_space_is_equal(src->ls->dim, dst->ls->dim))
2657                 isl_die(ctx, isl_error_invalid,
2658                         "spaces don't match", goto error);
2659
2660         if (src->ls->div->n_row == 0)
2661                 return dst;
2662
2663         exp1 = isl_alloc_array(ctx, int, src->ls->div->n_row);
2664         exp2 = isl_alloc_array(ctx, int, dst->ls->div->n_row);
2665         if (!exp1 || !exp2)
2666                 goto error;
2667
2668         div = isl_merge_divs(src->ls->div, dst->ls->div, exp1, exp2);
2669         dst = isl_aff_expand_divs(dst, div, exp2);
2670         free(exp1);
2671         free(exp2);
2672
2673         return dst;
2674 error:
2675         free(exp1);
2676         free(exp2);
2677         return isl_aff_free(dst);
2678 }
2679
2680 /* Adjust the local spaces of the affine expressions in "maff"
2681  * such that they all have the save divs.
2682  */
2683 __isl_give isl_multi_aff *isl_multi_aff_align_divs(
2684         __isl_take isl_multi_aff *maff)
2685 {
2686         int i;
2687
2688         if (!maff)
2689                 return NULL;
2690         if (maff->n == 0)
2691                 return maff;
2692         maff = isl_multi_aff_cow(maff);
2693         if (!maff)
2694                 return NULL;
2695
2696         for (i = 1; i < maff->n; ++i)
2697                 maff->p[0] = isl_aff_align_divs(maff->p[0], maff->p[i]);
2698         for (i = 1; i < maff->n; ++i) {
2699                 maff->p[i] = isl_aff_align_divs(maff->p[i], maff->p[0]);
2700                 if (!maff->p[i])
2701                         return isl_multi_aff_free(maff);
2702         }
2703
2704         return maff;
2705 }
2706
2707 __isl_give isl_aff *isl_aff_lift(__isl_take isl_aff *aff)
2708 {
2709         aff = isl_aff_cow(aff);
2710         if (!aff)
2711                 return NULL;
2712
2713         aff->ls = isl_local_space_lift(aff->ls);
2714         if (!aff->ls)
2715                 return isl_aff_free(aff);
2716
2717         return aff;
2718 }
2719
2720 /* Lift "maff" to a space with extra dimensions such that the result
2721  * has no more existentially quantified variables.
2722  * If "ls" is not NULL, then *ls is assigned the local space that lies
2723  * at the basis of the lifting applied to "maff".
2724  */
2725 __isl_give isl_multi_aff *isl_multi_aff_lift(__isl_take isl_multi_aff *maff,
2726         __isl_give isl_local_space **ls)
2727 {
2728         int i;
2729         isl_space *space;
2730         unsigned n_div;
2731
2732         if (ls)
2733                 *ls = NULL;
2734
2735         if (!maff)
2736                 return NULL;
2737
2738         if (maff->n == 0) {
2739                 if (ls) {
2740                         isl_space *space = isl_multi_aff_get_domain_space(maff);
2741                         *ls = isl_local_space_from_space(space);
2742                         if (!*ls)
2743                                 return isl_multi_aff_free(maff);
2744                 }
2745                 return maff;
2746         }
2747
2748         maff = isl_multi_aff_cow(maff);
2749         maff = isl_multi_aff_align_divs(maff);
2750         if (!maff)
2751                 return NULL;
2752
2753         n_div = isl_aff_dim(maff->p[0], isl_dim_div);
2754         space = isl_multi_aff_get_space(maff);
2755         space = isl_space_lift(isl_space_domain(space), n_div);
2756         space = isl_space_extend_domain_with_range(space,
2757                                                 isl_multi_aff_get_space(maff));
2758         if (!space)
2759                 return isl_multi_aff_free(maff);
2760         isl_space_free(maff->space);
2761         maff->space = space;
2762
2763         if (ls) {
2764                 *ls = isl_aff_get_domain_local_space(maff->p[0]);
2765                 if (!*ls)
2766                         return isl_multi_aff_free(maff);
2767         }
2768
2769         for (i = 0; i < maff->n; ++i) {
2770                 maff->p[i] = isl_aff_lift(maff->p[i]);
2771                 if (!maff->p[i])
2772                         goto error;
2773         }
2774
2775         return maff;
2776 error:
2777         if (ls)
2778                 isl_local_space_free(*ls);
2779         return isl_multi_aff_free(maff);
2780 }
2781
2782
2783 /* Extract an isl_pw_aff corresponding to output dimension "pos" of "pma".
2784  */
2785 __isl_give isl_pw_aff *isl_pw_multi_aff_get_pw_aff(
2786         __isl_keep isl_pw_multi_aff *pma, int pos)
2787 {
2788         int i;
2789         int n_out;
2790         isl_space *space;
2791         isl_pw_aff *pa;
2792
2793         if (!pma)
2794                 return NULL;
2795
2796         n_out = isl_pw_multi_aff_dim(pma, isl_dim_out);
2797         if (pos < 0 || pos >= n_out)
2798                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
2799                         "index out of bounds", return NULL);
2800
2801         space = isl_pw_multi_aff_get_space(pma);
2802         space = isl_space_drop_dims(space, isl_dim_out,
2803                                     pos + 1, n_out - pos - 1);
2804         space = isl_space_drop_dims(space, isl_dim_out, 0, pos);
2805
2806         pa = isl_pw_aff_alloc_size(space, pma->n);
2807         for (i = 0; i < pma->n; ++i) {
2808                 isl_aff *aff;
2809                 aff = isl_multi_aff_get_aff(pma->p[i].maff, pos);
2810                 pa = isl_pw_aff_add_piece(pa, isl_set_copy(pma->p[i].set), aff);
2811         }
2812
2813         return pa;
2814 }