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 MIT 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 #define ISL_DIM_H
  16 #include <isl_map_private.h>
  17 #include <isl_union_map_private.h>
  18 #include <isl_aff_private.h>
  19 #include <isl_space_private.h>
  20 #include <isl_local_space_private.h>
  21 #include <isl_mat_private.h>
  22 #include <isl_list_private.h>
  23 #include <isl/constraint.h>
  24 #include <isl/seq.h>
  25 #include <isl/set.h>
  26 #include <isl_config.h>
  27
  28 __isl_give isl_aff *isl_aff_alloc_vec(__isl_take isl_local_space *ls,
  29         __isl_take isl_vec *v)
  30 {
  31         isl_aff *aff;
  32
  33         if (!ls || !v)
  34                 goto error;
  35
  36         aff = isl_calloc_type(v->ctx, struct isl_aff);
  37         if (!aff)
  38                 goto error;
  39
  40         aff->ref = 1;
  41         aff->ls = ls;
  42         aff->v = v;
  43
  44         return aff;
  45 error:
  46         isl_local_space_free(ls);
  47         isl_vec_free(v);
  48         return NULL;
  49 }
  50
  51 __isl_give isl_aff *isl_aff_alloc(__isl_take isl_local_space *ls)
  52 {
  53         isl_ctx *ctx;
  54         isl_vec *v;
  55         unsigned total;
  56
  57         if (!ls)
  58                 return NULL;
  59
  60         ctx = isl_local_space_get_ctx(ls);
  61         if (!isl_local_space_divs_known(ls))
  62                 isl_die(ctx, isl_error_invalid, "local space has unknown divs",
  63                         goto error);
  64         if (!isl_local_space_is_set(ls))
  65                 isl_die(ctx, isl_error_invalid,
  66                         "domain of affine expression should be a set",
  67                         goto error);
  68
  69         total = isl_local_space_dim(ls, isl_dim_all);
  70         v = isl_vec_alloc(ctx, 1 + 1 + total);
  71         return isl_aff_alloc_vec(ls, v);
  72 error:
  73         isl_local_space_free(ls);
  74         return NULL;
  75 }
  76
  77 __isl_give isl_aff *isl_aff_zero_on_domain(__isl_take isl_local_space *ls)
  78 {
  79         isl_aff *aff;
  80
  81         aff = isl_aff_alloc(ls);
  82         if (!aff)
  83                 return NULL;
  84
  85         isl_int_set_si(aff->v->el[0], 1);
  86         isl_seq_clr(aff->v->el + 1, aff->v->size - 1);
  87
  88         return aff;
  89 }
  90
  91 __isl_give isl_aff *isl_aff_copy(__isl_keep isl_aff *aff)
  92 {
  93         if (!aff)
  94                 return NULL;
  95
  96         aff->ref++;
  97         return aff;
  98 }
  99
 100 __isl_give isl_aff *isl_aff_dup(__isl_keep isl_aff *aff)
 101 {
 102         if (!aff)
 103                 return NULL;
 104
 105         return isl_aff_alloc_vec(isl_local_space_copy(aff->ls),
 106                                  isl_vec_copy(aff->v));
 107 }
 108
 109 __isl_give isl_aff *isl_aff_cow(__isl_take isl_aff *aff)
 110 {
 111         if (!aff)
 112                 return NULL;
 113
 114         if (aff->ref == 1)
 115                 return aff;
 116         aff->ref--;
 117         return isl_aff_dup(aff);
 118 }
 119
 120 void *isl_aff_free(__isl_take isl_aff *aff)
 121 {
 122         if (!aff)
 123                 return NULL;
 124
 125         if (--aff->ref > 0)
 126                 return NULL;
 127
 128         isl_local_space_free(aff->ls);
 129         isl_vec_free(aff->v);
 130
 131         free(aff);
 132
 133         return NULL;
 134 }
 135
 136 isl_ctx *isl_aff_get_ctx(__isl_keep isl_aff *aff)
 137 {
 138         return aff ? isl_local_space_get_ctx(aff->ls) : NULL;
 139 }
 140
 141 /* Externally, an isl_aff has a map space, but internally, the
 142  * ls field corresponds to the domain of that space.
 143  */
 144 int isl_aff_dim(__isl_keep isl_aff *aff, enum isl_dim_type type)
 145 {
 146         if (!aff)
 147                 return 0;
 148         if (type == isl_dim_out)
 149                 return 1;
 150         if (type == isl_dim_in)
 151                 type = isl_dim_set;
 152         return isl_local_space_dim(aff->ls, type);
 153 }
 154
 155 __isl_give isl_space *isl_aff_get_domain_space(__isl_keep isl_aff *aff)
 156 {
 157         return aff ? isl_local_space_get_space(aff->ls) : NULL;
 158 }
 159
 160 __isl_give isl_space *isl_aff_get_space(__isl_keep isl_aff *aff)
 161 {
 162         isl_space *space;
 163         if (!aff)
 164                 return NULL;
 165         space = isl_local_space_get_space(aff->ls);
 166         space = isl_space_from_domain(space);
 167         space = isl_space_add_dims(space, isl_dim_out, 1);
 168         return space;
 169 }
 170
 171 __isl_give isl_local_space *isl_aff_get_domain_local_space(
 172         __isl_keep isl_aff *aff)
 173 {
 174         return aff ? isl_local_space_copy(aff->ls) : NULL;
 175 }
 176
 177 __isl_give isl_local_space *isl_aff_get_local_space(__isl_keep isl_aff *aff)
 178 {
 179         isl_local_space *ls;
 180         if (!aff)
 181                 return NULL;
 182         ls = isl_local_space_copy(aff->ls);
 183         ls = isl_local_space_from_domain(ls);
 184         ls = isl_local_space_add_dims(ls, isl_dim_out, 1);
 185         return ls;
 186 }
 187
 188 /* Externally, an isl_aff has a map space, but internally, the
 189  * ls field corresponds to the domain of that space.
 190  */
 191 const char *isl_aff_get_dim_name(__isl_keep isl_aff *aff,
 192         enum isl_dim_type type, unsigned pos)
 193 {
 194         if (!aff)
 195                 return NULL;
 196         if (type == isl_dim_out)
 197                 return NULL;
 198         if (type == isl_dim_in)
 199                 type = isl_dim_set;
 200         return isl_local_space_get_dim_name(aff->ls, type, pos);
 201 }
 202
 203 __isl_give isl_aff *isl_aff_reset_domain_space(__isl_take isl_aff *aff,
 204         __isl_take isl_space *dim)
 205 {
 206         aff = isl_aff_cow(aff);
 207         if (!aff || !dim)
 208                 goto error;
 209
 210         aff->ls = isl_local_space_reset_space(aff->ls, dim);
 211         if (!aff->ls)
 212                 return isl_aff_free(aff);
 213
 214         return aff;
 215 error:
 216         isl_aff_free(aff);
 217         isl_space_free(dim);
 218         return NULL;
 219 }
 220
 221 /* Reset the space of "aff".  This function is called from isl_pw_templ.c
 222  * and doesn't know if the space of an element object is represented
 223  * directly or through its domain.  It therefore passes along both.
 224  */
 225 __isl_give isl_aff *isl_aff_reset_space_and_domain(__isl_take isl_aff *aff,
 226         __isl_take isl_space *space, __isl_take isl_space *domain)
 227 {
 228         isl_space_free(space);
 229         return isl_aff_reset_domain_space(aff, domain);
 230 }
 231
 232 /* Reorder the coefficients of the affine expression based
 233  * on the given reodering.
 234  * The reordering r is assumed to have been extended with the local
 235  * variables.
 236  */
 237 static __isl_give isl_vec *vec_reorder(__isl_take isl_vec *vec,
 238         __isl_take isl_reordering *r, int n_div)
 239 {
 240         isl_vec *res;
 241         int i;
 242
 243         if (!vec || !r)
 244                 goto error;
 245
 246         res = isl_vec_alloc(vec->ctx,
 247                             2 + isl_space_dim(r->dim, isl_dim_all) + n_div);
 248         isl_seq_cpy(res->el, vec->el, 2);
 249         isl_seq_clr(res->el + 2, res->size - 2);
 250         for (i = 0; i < r->len; ++i)
 251                 isl_int_set(res->el[2 + r->pos[i]], vec->el[2 + i]);
 252
 253         isl_reordering_free(r);
 254         isl_vec_free(vec);
 255         return res;
 256 error:
 257         isl_vec_free(vec);
 258         isl_reordering_free(r);
 259         return NULL;
 260 }
 261
 262 /* Reorder the dimensions of the domain of "aff" according
 263  * to the given reordering.
 264  */
 265 __isl_give isl_aff *isl_aff_realign_domain(__isl_take isl_aff *aff,
 266         __isl_take isl_reordering *r)
 267 {
 268         aff = isl_aff_cow(aff);
 269         if (!aff)
 270                 goto error;
 271
 272         r = isl_reordering_extend(r, aff->ls->div->n_row);
 273         aff->v = vec_reorder(aff->v, isl_reordering_copy(r),
 274                                 aff->ls->div->n_row);
 275         aff->ls = isl_local_space_realign(aff->ls, r);
 276
 277         if (!aff->v || !aff->ls)
 278                 return isl_aff_free(aff);
 279
 280         return aff;
 281 error:
 282         isl_aff_free(aff);
 283         isl_reordering_free(r);
 284         return NULL;
 285 }
 286
 287 __isl_give isl_aff *isl_aff_align_params(__isl_take isl_aff *aff,
 288         __isl_take isl_space *model)
 289 {
 290         if (!aff || !model)
 291                 goto error;
 292
 293         if (!isl_space_match(aff->ls->dim, isl_dim_param,
 294                              model, isl_dim_param)) {
 295                 isl_reordering *exp;
 296
 297                 model = isl_space_drop_dims(model, isl_dim_in,
 298                                         0, isl_space_dim(model, isl_dim_in));
 299                 model = isl_space_drop_dims(model, isl_dim_out,
 300                                         0, isl_space_dim(model, isl_dim_out));
 301                 exp = isl_parameter_alignment_reordering(aff->ls->dim, model);
 302                 exp = isl_reordering_extend_space(exp,
 303                                         isl_aff_get_domain_space(aff));
 304                 aff = isl_aff_realign_domain(aff, exp);
 305         }
 306
 307         isl_space_free(model);
 308         return aff;
 309 error:
 310         isl_space_free(model);
 311         isl_aff_free(aff);
 312         return NULL;
 313 }
 314
 315 int isl_aff_plain_is_zero(__isl_keep isl_aff *aff)
 316 {
 317         if (!aff)
 318                 return -1;
 319
 320         return isl_seq_first_non_zero(aff->v->el + 1, aff->v->size - 1) < 0;
 321 }
 322
 323 int isl_aff_plain_is_equal(__isl_keep isl_aff *aff1, __isl_keep isl_aff *aff2)
 324 {
 325         int equal;
 326
 327         if (!aff1 || !aff2)
 328                 return -1;
 329
 330         equal = isl_local_space_is_equal(aff1->ls, aff2->ls);
 331         if (equal < 0 || !equal)
 332                 return equal;
 333
 334         return isl_vec_is_equal(aff1->v, aff2->v);
 335 }
 336
 337 int isl_aff_get_denominator(__isl_keep isl_aff *aff, isl_int *v)
 338 {
 339         if (!aff)
 340                 return -1;
 341         isl_int_set(*v, aff->v->el[0]);
 342         return 0;
 343 }
 344
 345 int isl_aff_get_constant(__isl_keep isl_aff *aff, isl_int *v)
 346 {
 347         if (!aff)
 348                 return -1;
 349         isl_int_set(*v, aff->v->el[1]);
 350         return 0;
 351 }
 352
 353 int isl_aff_get_coefficient(__isl_keep isl_aff *aff,
 354         enum isl_dim_type type, int pos, isl_int *v)
 355 {
 356         if (!aff)
 357                 return -1;
 358
 359         if (type == isl_dim_out)
 360                 isl_die(aff->v->ctx, isl_error_invalid,
 361                         "output/set dimension does not have a coefficient",
 362                         return -1);
 363         if (type == isl_dim_in)
 364                 type = isl_dim_set;
 365
 366         if (pos >= isl_local_space_dim(aff->ls, type))
 367                 isl_die(aff->v->ctx, isl_error_invalid,
 368                         "position out of bounds", return -1);
 369
 370         pos += isl_local_space_offset(aff->ls, type);
 371         isl_int_set(*v, aff->v->el[1 + pos]);
 372
 373         return 0;
 374 }
 375
 376 __isl_give isl_aff *isl_aff_set_denominator(__isl_take isl_aff *aff, isl_int v)
 377 {
 378         aff = isl_aff_cow(aff);
 379         if (!aff)
 380                 return NULL;
 381
 382         aff->v = isl_vec_cow(aff->v);
 383         if (!aff->v)
 384                 return isl_aff_free(aff);
 385
 386         isl_int_set(aff->v->el[0], v);
 387
 388         return aff;
 389 }
 390
 391 __isl_give isl_aff *isl_aff_set_constant(__isl_take isl_aff *aff, isl_int v)
 392 {
 393         aff = isl_aff_cow(aff);
 394         if (!aff)
 395                 return NULL;
 396
 397         aff->v = isl_vec_cow(aff->v);
 398         if (!aff->v)
 399                 return isl_aff_free(aff);
 400
 401         isl_int_set(aff->v->el[1], v);
 402
 403         return aff;
 404 }
 405
 406 __isl_give isl_aff *isl_aff_add_constant(__isl_take isl_aff *aff, isl_int v)
 407 {
 408         if (isl_int_is_zero(v))
 409                 return aff;
 410
 411         aff = isl_aff_cow(aff);
 412         if (!aff)
 413                 return NULL;
 414
 415         aff->v = isl_vec_cow(aff->v);
 416         if (!aff->v)
 417                 return isl_aff_free(aff);
 418
 419         isl_int_addmul(aff->v->el[1], aff->v->el[0], v);
 420
 421         return aff;
 422 }
 423
 424 __isl_give isl_aff *isl_aff_add_constant_si(__isl_take isl_aff *aff, int v)
 425 {
 426         isl_int t;
 427
 428         isl_int_init(t);
 429         isl_int_set_si(t, v);
 430         aff = isl_aff_add_constant(aff, t);
 431         isl_int_clear(t);
 432
 433         return aff;
 434 }
 435
 436 /* Add "v" to the numerator of the constant term of "aff".
 437  */
 438 __isl_give isl_aff *isl_aff_add_constant_num(__isl_take isl_aff *aff, isl_int v)
 439 {
 440         if (isl_int_is_zero(v))
 441                 return aff;
 442
 443         aff = isl_aff_cow(aff);
 444         if (!aff)
 445                 return NULL;
 446
 447         aff->v = isl_vec_cow(aff->v);
 448         if (!aff->v)
 449                 return isl_aff_free(aff);
 450
 451         isl_int_add(aff->v->el[1], aff->v->el[1], v);
 452
 453         return aff;
 454 }
 455
 456 /* Add "v" to the numerator of the constant term of "aff".
 457  */
 458 __isl_give isl_aff *isl_aff_add_constant_num_si(__isl_take isl_aff *aff, int v)
 459 {
 460         isl_int t;
 461
 462         if (v == 0)
 463                 return aff;
 464
 465         isl_int_init(t);
 466         isl_int_set_si(t, v);
 467         aff = isl_aff_add_constant_num(aff, t);
 468         isl_int_clear(t);
 469
 470         return aff;
 471 }
 472
 473 __isl_give isl_aff *isl_aff_set_constant_si(__isl_take isl_aff *aff, int v)
 474 {
 475         aff = isl_aff_cow(aff);
 476         if (!aff)
 477                 return NULL;
 478
 479         aff->v = isl_vec_cow(aff->v);
 480         if (!aff->v)
 481                 return isl_aff_free(aff);
 482
 483         isl_int_set_si(aff->v->el[1], v);
 484
 485         return aff;
 486 }
 487
 488 __isl_give isl_aff *isl_aff_set_coefficient(__isl_take isl_aff *aff,
 489         enum isl_dim_type type, int pos, isl_int v)
 490 {
 491         if (!aff)
 492                 return NULL;
 493
 494         if (type == isl_dim_out)
 495                 isl_die(aff->v->ctx, isl_error_invalid,
 496                         "output/set dimension does not have a coefficient",
 497                         return isl_aff_free(aff));
 498         if (type == isl_dim_in)
 499                 type = isl_dim_set;
 500
 501         if (pos >= isl_local_space_dim(aff->ls, type))
 502                 isl_die(aff->v->ctx, isl_error_invalid,
 503                         "position out of bounds", return isl_aff_free(aff));
 504
 505         aff = isl_aff_cow(aff);
 506         if (!aff)
 507                 return NULL;
 508
 509         aff->v = isl_vec_cow(aff->v);
 510         if (!aff->v)
 511                 return isl_aff_free(aff);
 512
 513         pos += isl_local_space_offset(aff->ls, type);
 514         isl_int_set(aff->v->el[1 + pos], v);
 515
 516         return aff;
 517 }
 518
 519 __isl_give isl_aff *isl_aff_set_coefficient_si(__isl_take isl_aff *aff,
 520         enum isl_dim_type type, int pos, int v)
 521 {
 522         if (!aff)
 523                 return NULL;
 524
 525         if (type == isl_dim_out)
 526                 isl_die(aff->v->ctx, isl_error_invalid,
 527                         "output/set dimension does not have a coefficient",
 528                         return isl_aff_free(aff));
 529         if (type == isl_dim_in)
 530                 type = isl_dim_set;
 531
 532         if (pos >= isl_local_space_dim(aff->ls, type))
 533                 isl_die(aff->v->ctx, isl_error_invalid,
 534                         "position out of bounds", return isl_aff_free(aff));
 535
 536         aff = isl_aff_cow(aff);
 537         if (!aff)
 538                 return NULL;
 539
 540         aff->v = isl_vec_cow(aff->v);
 541         if (!aff->v)
 542                 return isl_aff_free(aff);
 543
 544         pos += isl_local_space_offset(aff->ls, type);
 545         isl_int_set_si(aff->v->el[1 + pos], v);
 546
 547         return aff;
 548 }
 549
 550 __isl_give isl_aff *isl_aff_add_coefficient(__isl_take isl_aff *aff,
 551         enum isl_dim_type type, int pos, isl_int v)
 552 {
 553         if (!aff)
 554                 return NULL;
 555
 556         if (type == isl_dim_out)
 557                 isl_die(aff->v->ctx, isl_error_invalid,
 558                         "output/set dimension does not have a coefficient",
 559                         return isl_aff_free(aff));
 560         if (type == isl_dim_in)
 561                 type = isl_dim_set;
 562
 563         if (pos >= isl_local_space_dim(aff->ls, type))
 564                 isl_die(aff->v->ctx, isl_error_invalid,
 565                         "position out of bounds", return isl_aff_free(aff));
 566
 567         aff = isl_aff_cow(aff);
 568         if (!aff)
 569                 return NULL;
 570
 571         aff->v = isl_vec_cow(aff->v);
 572         if (!aff->v)
 573                 return isl_aff_free(aff);
 574
 575         pos += isl_local_space_offset(aff->ls, type);
 576         isl_int_addmul(aff->v->el[1 + pos], aff->v->el[0], v);
 577
 578         return aff;
 579 }
 580
 581 __isl_give isl_aff *isl_aff_add_coefficient_si(__isl_take isl_aff *aff,
 582         enum isl_dim_type type, int pos, int v)
 583 {
 584         isl_int t;
 585
 586         isl_int_init(t);
 587         isl_int_set_si(t, v);
 588         aff = isl_aff_add_coefficient(aff, type, pos, t);
 589         isl_int_clear(t);
 590
 591         return aff;
 592 }
 593
 594 __isl_give isl_aff *isl_aff_get_div(__isl_keep isl_aff *aff, int pos)
 595 {
 596         if (!aff)
 597                 return NULL;
 598
 599         return isl_local_space_get_div(aff->ls, pos);
 600 }
 601
 602 __isl_give isl_aff *isl_aff_neg(__isl_take isl_aff *aff)
 603 {
 604         aff = isl_aff_cow(aff);
 605         if (!aff)
 606                 return NULL;
 607         aff->v = isl_vec_cow(aff->v);
 608         if (!aff->v)
 609                 return isl_aff_free(aff);
 610
 611         isl_seq_neg(aff->v->el + 1, aff->v->el + 1, aff->v->size - 1);
 612
 613         return aff;
 614 }
 615
 616 /* Remove divs from the local space that do not appear in the affine
 617  * expression.
 618  * We currently only remove divs at the end.
 619  * Some intermediate divs may also not appear directly in the affine
 620  * expression, but we would also need to check that no other divs are
 621  * defined in terms of them.
 622  */
 623 __isl_give isl_aff *isl_aff_remove_unused_divs( __isl_take isl_aff *aff)
 624 {
 625         int pos;
 626         int off;
 627         int n;
 628
 629         if (!aff)
 630                 return NULL;
 631
 632         n = isl_local_space_dim(aff->ls, isl_dim_div);
 633         off = isl_local_space_offset(aff->ls, isl_dim_div);
 634
 635         pos = isl_seq_last_non_zero(aff->v->el + 1 + off, n) + 1;
 636         if (pos == n)
 637                 return aff;
 638
 639         aff = isl_aff_cow(aff);
 640         if (!aff)
 641                 return NULL;
 642
 643         aff->ls = isl_local_space_drop_dims(aff->ls, isl_dim_div, pos, n - pos);
 644         aff->v = isl_vec_drop_els(aff->v, 1 + off + pos, n - pos);
 645         if (!aff->ls || !aff->v)
 646                 return isl_aff_free(aff);
 647
 648         return aff;
 649 }
 650
 651 /* Given two affine expressions "p" of length p_len (including the
 652  * denominator and the constant term) and "subs" of length subs_len,
 653  * plug in "subs" for the variable at position "pos".
 654  * The variables of "subs" and "p" are assumed to match up to subs_len,
 655  * but "p" may have additional variables.
 656  * "v" is an initialized isl_int that can be used internally.
 657  *
 658  * In particular, if "p" represents the expression
 659  *
 660  *      (a i + g)/m
 661  *
 662  * with i the variable at position "pos" and "subs" represents the expression
 663  *
 664  *      f/d
 665  *
 666  * then the result represents the expression
 667  *
 668  *      (a f + d g)/(m d)
 669  *
 670  */
 671 void isl_seq_substitute(isl_int *p, int pos, isl_int *subs,
 672         int p_len, int subs_len, isl_int v)
 673 {
 674         isl_int_set(v, p[1 + pos]);
 675         isl_int_set_si(p[1 + pos], 0);
 676         isl_seq_combine(p + 1, subs[0], p + 1, v, subs + 1, subs_len - 1);
 677         isl_seq_scale(p + subs_len, p + subs_len, subs[0], p_len - subs_len);
 678         isl_int_mul(p[0], p[0], subs[0]);
 679 }
 680
 681 /* Look for any divs in the aff->ls with a denominator equal to one
 682  * and plug them into the affine expression and any subsequent divs
 683  * that may reference the div.
 684  */
 685 static __isl_give isl_aff *plug_in_integral_divs(__isl_take isl_aff *aff)
 686 {
 687         int i, n;
 688         int len;
 689         isl_int v;
 690         isl_vec *vec;
 691         isl_local_space *ls;
 692         unsigned pos;
 693
 694         if (!aff)
 695                 return NULL;
 696
 697         n = isl_local_space_dim(aff->ls, isl_dim_div);
 698         len = aff->v->size;
 699         for (i = 0; i < n; ++i) {
 700                 if (!isl_int_is_one(aff->ls->div->row[i][0]))
 701                         continue;
 702                 ls = isl_local_space_copy(aff->ls);
 703                 ls = isl_local_space_substitute_seq(ls, isl_dim_div, i,
 704                                             aff->ls->div->row[i], len, i + 1);
 705                 vec = isl_vec_copy(aff->v);
 706                 vec = isl_vec_cow(vec);
 707                 if (!ls || !vec)
 708                         goto error;
 709
 710                 isl_int_init(v);
 711
 712                 pos = isl_local_space_offset(aff->ls, isl_dim_div) + i;
 713                 isl_seq_substitute(vec->el, pos, aff->ls->div->row[i],
 714                                         len, len, v);
 715
 716                 isl_int_clear(v);
 717
 718                 isl_vec_free(aff->v);
 719                 aff->v = vec;
 720                 isl_local_space_free(aff->ls);
 721                 aff->ls = ls;
 722         }
 723
 724         return aff;
 725 error:
 726         isl_vec_free(vec);
 727         isl_local_space_free(ls);
 728         return isl_aff_free(aff);
 729 }
 730
 731 /* Swap divs "a" and "b" in "aff", which is assumed to be non-NULL.
 732  *
 733  * Even though this function is only called on isl_affs with a single
 734  * reference, we are careful to only change aff->v and aff->ls together.
 735  */
 736 static __isl_give isl_aff *swap_div(__isl_take isl_aff *aff, int a, int b)
 737 {
 738         unsigned off = isl_local_space_offset(aff->ls, isl_dim_div);
 739         isl_local_space *ls;
 740         isl_vec *v;
 741
 742         ls = isl_local_space_copy(aff->ls);
 743         ls = isl_local_space_swap_div(ls, a, b);
 744         v = isl_vec_copy(aff->v);
 745         v = isl_vec_cow(v);
 746         if (!ls || !v)
 747                 goto error;
 748
 749         isl_int_swap(v->el[1 + off + a], v->el[1 + off + b]);
 750         isl_vec_free(aff->v);
 751         aff->v = v;
 752         isl_local_space_free(aff->ls);
 753         aff->ls = ls;
 754
 755         return aff;
 756 error:
 757         isl_vec_free(v);
 758         isl_local_space_free(ls);
 759         return isl_aff_free(aff);
 760 }
 761
 762 /* Sort the divs in the local space of "aff" according to
 763  * the comparison function "cmp_row" in isl_local_space.c
 764  *
 765  * Reordering divs does not change the semantics of "aff",
 766  * so there is no need to call isl_aff_cow.
 767  * Moreover, this function is currently only called on isl_affs
 768  * with a single reference.
 769  */
 770 static __isl_give isl_aff *sort_divs(__isl_take isl_aff *aff)
 771 {
 772         int i, j, n;
 773
 774         if (!aff)
 775                 return NULL;
 776
 777         n = isl_aff_dim(aff, isl_dim_div);
 778         for (i = 1; i < n; ++i) {
 779                 for (j = i - 1; j >= 0; --j) {
 780                         if (isl_mat_cmp_div(aff->ls->div, j, j + 1) <= 0)
 781                                 break;
 782                         aff = swap_div(aff, j, j + 1);
 783                         if (!aff)
 784                                 return NULL;
 785                 }
 786         }
 787
 788         return aff;
 789 }
 790
 791 /* Normalize the representation of "aff".
 792  *
 793  * This function should only be called of "new" isl_affs, i.e.,
 794  * with only a single reference.  We therefore do not need to
 795  * worry about affecting other instances.
 796  */
 797 __isl_give isl_aff *isl_aff_normalize(__isl_take isl_aff *aff)
 798 {
 799         if (!aff)
 800                 return NULL;
 801         aff->v = isl_vec_normalize(aff->v);
 802         if (!aff->v)
 803                 return isl_aff_free(aff);
 804         aff = plug_in_integral_divs(aff);
 805         aff = sort_divs(aff);
 806         aff = isl_aff_remove_unused_divs(aff);
 807         return aff;
 808 }
 809
 810 /* Given f, return floor(f).
 811  * If f is an integer expression, then just return f.
 812  * If f is a constant, then return the constant floor(f).
 813  * Otherwise, if f = g/m, write g = q m + r,
 814  * create a new div d = [r/m] and return the expression q + d.
 815  * The coefficients in r are taken to lie between -m/2 and m/2.
 816  */
 817 __isl_give isl_aff *isl_aff_floor(__isl_take isl_aff *aff)
 818 {
 819         int i;
 820         int size;
 821         isl_ctx *ctx;
 822         isl_vec *div;
 823
 824         if (!aff)
 825                 return NULL;
 826
 827         if (isl_int_is_one(aff->v->el[0]))
 828                 return aff;
 829
 830         aff = isl_aff_cow(aff);
 831         if (!aff)
 832                 return NULL;
 833
 834         aff->v = isl_vec_cow(aff->v);
 835         if (!aff->v)
 836                 return isl_aff_free(aff);
 837
 838         if (isl_aff_is_cst(aff)) {
 839                 isl_int_fdiv_q(aff->v->el[1], aff->v->el[1], aff->v->el[0]);
 840                 isl_int_set_si(aff->v->el[0], 1);
 841                 return aff;
 842         }
 843
 844         div = isl_vec_copy(aff->v);
 845         div = isl_vec_cow(div);
 846         if (!div)
 847                 return isl_aff_free(aff);
 848
 849         ctx = isl_aff_get_ctx(aff);
 850         isl_int_fdiv_q(aff->v->el[0], aff->v->el[0], ctx->two);
 851         for (i = 1; i < aff->v->size; ++i) {
 852                 isl_int_fdiv_r(div->el[i], div->el[i], div->el[0]);
 853                 isl_int_fdiv_q(aff->v->el[i], aff->v->el[i], div->el[0]);
 854                 if (isl_int_gt(div->el[i], aff->v->el[0])) {
 855                         isl_int_sub(div->el[i], div->el[i], div->el[0]);
 856                         isl_int_add_ui(aff->v->el[i], aff->v->el[i], 1);
 857                 }
 858         }
 859
 860         aff->ls = isl_local_space_add_div(aff->ls, div);
 861         if (!aff->ls)
 862                 return isl_aff_free(aff);
 863
 864         size = aff->v->size;
 865         aff->v = isl_vec_extend(aff->v, size + 1);
 866         if (!aff->v)
 867                 return isl_aff_free(aff);
 868         isl_int_set_si(aff->v->el[0], 1);
 869         isl_int_set_si(aff->v->el[size], 1);
 870
 871         return aff;
 872 }
 873
 874 /* Compute
 875  *
 876  *      aff mod m = aff - m * floor(aff/m)
 877  */
 878 __isl_give isl_aff *isl_aff_mod(__isl_take isl_aff *aff, isl_int m)
 879 {
 880         isl_aff *res;
 881
 882         res = isl_aff_copy(aff);
 883         aff = isl_aff_scale_down(aff, m);
 884         aff = isl_aff_floor(aff);
 885         aff = isl_aff_scale(aff, m);
 886         res = isl_aff_sub(res, aff);
 887
 888         return res;
 889 }
 890
 891 /* Compute
 892  *
 893  *      pwaff mod m = pwaff - m * floor(pwaff/m)
 894  */
 895 __isl_give isl_pw_aff *isl_pw_aff_mod(__isl_take isl_pw_aff *pwaff, isl_int m)
 896 {
 897         isl_pw_aff *res;
 898
 899         res = isl_pw_aff_copy(pwaff);
 900         pwaff = isl_pw_aff_scale_down(pwaff, m);
 901         pwaff = isl_pw_aff_floor(pwaff);
 902         pwaff = isl_pw_aff_scale(pwaff, m);
 903         res = isl_pw_aff_sub(res, pwaff);
 904
 905         return res;
 906 }
 907
 908 /* Given f, return ceil(f).
 909  * If f is an integer expression, then just return f.
 910  * Otherwise, create a new div d = [-f] and return the expression -d.
 911  */
 912 __isl_give isl_aff *isl_aff_ceil(__isl_take isl_aff *aff)
 913 {
 914         if (!aff)
 915                 return NULL;
 916
 917         if (isl_int_is_one(aff->v->el[0]))
 918                 return aff;
 919
 920         aff = isl_aff_neg(aff);
 921         aff = isl_aff_floor(aff);
 922         aff = isl_aff_neg(aff);
 923
 924         return aff;
 925 }
 926
 927 /* Apply the expansion computed by isl_merge_divs.
 928  * The expansion itself is given by "exp" while the resulting
 929  * list of divs is given by "div".
 930  */
 931 __isl_give isl_aff *isl_aff_expand_divs( __isl_take isl_aff *aff,
 932         __isl_take isl_mat *div, int *exp)
 933 {
 934         int i, j;
 935         int old_n_div;
 936         int new_n_div;
 937         int offset;
 938
 939         aff = isl_aff_cow(aff);
 940         if (!aff || !div)
 941                 goto error;
 942
 943         old_n_div = isl_local_space_dim(aff->ls, isl_dim_div);
 944         new_n_div = isl_mat_rows(div);
 945         if (new_n_div < old_n_div)
 946                 isl_die(isl_mat_get_ctx(div), isl_error_invalid,
 947                         "not an expansion", goto error);
 948
 949         aff->v = isl_vec_extend(aff->v, aff->v->size + new_n_div - old_n_div);
 950         if (!aff->v)
 951                 goto error;
 952
 953         offset = 1 + isl_local_space_offset(aff->ls, isl_dim_div);
 954         j = old_n_div - 1;
 955         for (i = new_n_div - 1; i >= 0; --i) {
 956                 if (j >= 0 && exp[j] == i) {
 957                         if (i != j)
 958                                 isl_int_swap(aff->v->el[offset + i],
 959                                              aff->v->el[offset + j]);
 960                         j--;
 961                 } else
 962                         isl_int_set_si(aff->v->el[offset + i], 0);
 963         }
 964
 965         aff->ls = isl_local_space_replace_divs(aff->ls, isl_mat_copy(div));
 966         if (!aff->ls)
 967                 goto error;
 968         isl_mat_free(div);
 969         return aff;
 970 error:
 971         isl_aff_free(aff);
 972         isl_mat_free(div);
 973         return NULL;
 974 }
 975
 976 /* Add two affine expressions that live in the same local space.
 977  */
 978 static __isl_give isl_aff *add_expanded(__isl_take isl_aff *aff1,
 979         __isl_take isl_aff *aff2)
 980 {
 981         isl_int gcd, f;
 982
 983         aff1 = isl_aff_cow(aff1);
 984         if (!aff1 || !aff2)
 985                 goto error;
 986
 987         aff1->v = isl_vec_cow(aff1->v);
 988         if (!aff1->v)
 989                 goto error;
 990
 991         isl_int_init(gcd);
 992         isl_int_init(f);
 993         isl_int_gcd(gcd, aff1->v->el[0], aff2->v->el[0]);
 994         isl_int_divexact(f, aff2->v->el[0], gcd);
 995         isl_seq_scale(aff1->v->el + 1, aff1->v->el + 1, f, aff1->v->size - 1);
 996         isl_int_divexact(f, aff1->v->el[0], gcd);
 997         isl_seq_addmul(aff1->v->el + 1, f, aff2->v->el + 1, aff1->v->size - 1);
 998         isl_int_divexact(f, aff2->v->el[0], gcd);
 999         isl_int_mul(aff1->v->el[0], aff1->v->el[0], f);
1000         isl_int_clear(f);
1001         isl_int_clear(gcd);
1002
1003         isl_aff_free(aff2);
1004         return aff1;
1005 error:
1006         isl_aff_free(aff1);
1007         isl_aff_free(aff2);
1008         return NULL;
1009 }
1010
1011 __isl_give isl_aff *isl_aff_add(__isl_take isl_aff *aff1,
1012         __isl_take isl_aff *aff2)
1013 {
1014         isl_ctx *ctx;
1015         int *exp1 = NULL;
1016         int *exp2 = NULL;
1017         isl_mat *div;
1018
1019         if (!aff1 || !aff2)
1020                 goto error;
1021
1022         ctx = isl_aff_get_ctx(aff1);
1023         if (!isl_space_is_equal(aff1->ls->dim, aff2->ls->dim))
1024                 isl_die(ctx, isl_error_invalid,
1025                         "spaces don't match", goto error);
1026
1027         if (aff1->ls->div->n_row == 0 && aff2->ls->div->n_row == 0)
1028                 return add_expanded(aff1, aff2);
1029
1030         exp1 = isl_alloc_array(ctx, int, aff1->ls->div->n_row);
1031         exp2 = isl_alloc_array(ctx, int, aff2->ls->div->n_row);
1032         if (!exp1 || !exp2)
1033                 goto error;
1034
1035         div = isl_merge_divs(aff1->ls->div, aff2->ls->div, exp1, exp2);
1036         aff1 = isl_aff_expand_divs(aff1, isl_mat_copy(div), exp1);
1037         aff2 = isl_aff_expand_divs(aff2, div, exp2);
1038         free(exp1);
1039         free(exp2);
1040
1041         return add_expanded(aff1, aff2);
1042 error:
1043         free(exp1);
1044         free(exp2);
1045         isl_aff_free(aff1);
1046         isl_aff_free(aff2);
1047         return NULL;
1048 }
1049
1050 __isl_give isl_aff *isl_aff_sub(__isl_take isl_aff *aff1,
1051         __isl_take isl_aff *aff2)
1052 {
1053         return isl_aff_add(aff1, isl_aff_neg(aff2));
1054 }
1055
1056 __isl_give isl_aff *isl_aff_scale(__isl_take isl_aff *aff, isl_int f)
1057 {
1058         isl_int gcd;
1059
1060         if (isl_int_is_one(f))
1061                 return aff;
1062
1063         aff = isl_aff_cow(aff);
1064         if (!aff)
1065                 return NULL;
1066         aff->v = isl_vec_cow(aff->v);
1067         if (!aff->v)
1068                 return isl_aff_free(aff);
1069
1070         isl_int_init(gcd);
1071         isl_int_gcd(gcd, aff->v->el[0], f);
1072         isl_int_divexact(aff->v->el[0], aff->v->el[0], gcd);
1073         isl_int_divexact(gcd, f, gcd);
1074         isl_seq_scale(aff->v->el + 1, aff->v->el + 1, gcd, aff->v->size - 1);
1075         isl_int_clear(gcd);
1076
1077         return aff;
1078 }
1079
1080 __isl_give isl_aff *isl_aff_scale_down(__isl_take isl_aff *aff, isl_int f)
1081 {
1082         isl_int gcd;
1083
1084         if (isl_int_is_one(f))
1085                 return aff;
1086
1087         aff = isl_aff_cow(aff);
1088         if (!aff)
1089                 return NULL;
1090
1091         if (isl_int_is_zero(f))
1092                 isl_die(isl_aff_get_ctx(aff), isl_error_invalid,
1093                         "cannot scale down by zero", return isl_aff_free(aff));
1094
1095         aff->v = isl_vec_cow(aff->v);
1096         if (!aff->v)
1097                 return isl_aff_free(aff);
1098
1099         isl_int_init(gcd);
1100         isl_seq_gcd(aff->v->el + 1, aff->v->size - 1, &gcd);
1101         isl_int_gcd(gcd, gcd, f);
1102         isl_seq_scale_down(aff->v->el + 1, aff->v->el + 1, gcd, aff->v->size - 1);
1103         isl_int_divexact(gcd, f, gcd);
1104         isl_int_mul(aff->v->el[0], aff->v->el[0], gcd);
1105         isl_int_clear(gcd);
1106
1107         return aff;
1108 }
1109
1110 __isl_give isl_aff *isl_aff_scale_down_ui(__isl_take isl_aff *aff, unsigned f)
1111 {
1112         isl_int v;
1113
1114         if (f == 1)
1115                 return aff;
1116
1117         isl_int_init(v);
1118         isl_int_set_ui(v, f);
1119         aff = isl_aff_scale_down(aff, v);
1120         isl_int_clear(v);
1121
1122         return aff;
1123 }
1124
1125 __isl_give isl_aff *isl_aff_set_dim_name(__isl_take isl_aff *aff,
1126         enum isl_dim_type type, unsigned pos, const char *s)
1127 {
1128         aff = isl_aff_cow(aff);
1129         if (!aff)
1130                 return NULL;
1131         if (type == isl_dim_out)
1132                 isl_die(aff->v->ctx, isl_error_invalid,
1133                         "cannot set name of output/set dimension",
1134                         return isl_aff_free(aff));
1135         if (type == isl_dim_in)
1136                 type = isl_dim_set;
1137         aff->ls = isl_local_space_set_dim_name(aff->ls, type, pos, s);
1138         if (!aff->ls)
1139                 return isl_aff_free(aff);
1140
1141         return aff;
1142 }
1143
1144 __isl_give isl_aff *isl_aff_set_dim_id(__isl_take isl_aff *aff,
1145         enum isl_dim_type type, unsigned pos, __isl_take isl_id *id)
1146 {
1147         aff = isl_aff_cow(aff);
1148         if (!aff)
1149                 return isl_id_free(id);
1150         if (type == isl_dim_out)
1151                 isl_die(aff->v->ctx, isl_error_invalid,
1152                         "cannot set name of output/set dimension",
1153                         goto error);
1154         if (type == isl_dim_in)
1155                 type = isl_dim_set;
1156         aff->ls = isl_local_space_set_dim_id(aff->ls, type, pos, id);
1157         if (!aff->ls)
1158                 return isl_aff_free(aff);
1159
1160         return aff;
1161 error:
1162         isl_id_free(id);
1163         isl_aff_free(aff);
1164         return NULL;
1165 }
1166
1167 /* Exploit the equalities in "eq" to simplify the affine expression
1168  * and the expressions of the integer divisions in the local space.
1169  * The integer divisions in this local space are assumed to appear
1170  * as regular dimensions in "eq".
1171  */
1172 static __isl_give isl_aff *isl_aff_substitute_equalities_lifted(
1173         __isl_take isl_aff *aff, __isl_take isl_basic_set *eq)
1174 {
1175         int i, j;
1176         unsigned total;
1177         unsigned n_div;
1178
1179         if (!eq)
1180                 goto error;
1181         if (eq->n_eq == 0) {
1182                 isl_basic_set_free(eq);
1183                 return aff;
1184         }
1185
1186         aff = isl_aff_cow(aff);
1187         if (!aff)
1188                 goto error;
1189
1190         aff->ls = isl_local_space_substitute_equalities(aff->ls,
1191                                                         isl_basic_set_copy(eq));
1192         if (!aff->ls)
1193                 goto error;
1194
1195         total = 1 + isl_space_dim(eq->dim, isl_dim_all);
1196         n_div = eq->n_div;
1197         for (i = 0; i < eq->n_eq; ++i) {
1198                 j = isl_seq_last_non_zero(eq->eq[i], total + n_div);
1199                 if (j < 0 || j == 0 || j >= total)
1200                         continue;
1201
1202                 isl_seq_elim(aff->v->el + 1, eq->eq[i], j, total,
1203                                 &aff->v->el[0]);
1204         }
1205
1206         isl_basic_set_free(eq);
1207         aff = isl_aff_normalize(aff);
1208         return aff;
1209 error:
1210         isl_basic_set_free(eq);
1211         isl_aff_free(aff);
1212         return NULL;
1213 }
1214
1215 /* Exploit the equalities in "eq" to simplify the affine expression
1216  * and the expressions of the integer divisions in the local space.
1217  */
1218 static __isl_give isl_aff *isl_aff_substitute_equalities(
1219         __isl_take isl_aff *aff, __isl_take isl_basic_set *eq)
1220 {
1221         int n_div;
1222
1223         if (!aff || !eq)
1224                 goto error;
1225         n_div = isl_local_space_dim(aff->ls, isl_dim_div);
1226         if (n_div > 0)
1227                 eq = isl_basic_set_add(eq, isl_dim_set, n_div);
1228         return isl_aff_substitute_equalities_lifted(aff, eq);
1229 error:
1230         isl_basic_set_free(eq);
1231         isl_aff_free(aff);
1232         return NULL;
1233 }
1234
1235 /* Look for equalities among the variables shared by context and aff
1236  * and the integer divisions of aff, if any.
1237  * The equalities are then used to eliminate coefficients and/or integer
1238  * divisions from aff.
1239  */
1240 __isl_give isl_aff *isl_aff_gist(__isl_take isl_aff *aff,
1241         __isl_take isl_set *context)
1242 {
1243         isl_basic_set *hull;
1244         int n_div;
1245
1246         if (!aff)
1247                 goto error;
1248         n_div = isl_local_space_dim(aff->ls, isl_dim_div);
1249         if (n_div > 0) {
1250                 isl_basic_set *bset;
1251                 isl_local_space *ls;
1252                 context = isl_set_add_dims(context, isl_dim_set, n_div);
1253                 ls = isl_aff_get_domain_local_space(aff);
1254                 bset = isl_basic_set_from_local_space(ls);
1255                 bset = isl_basic_set_lift(bset);
1256                 bset = isl_basic_set_flatten(bset);
1257                 context = isl_set_intersect(context,
1258                                             isl_set_from_basic_set(bset));
1259         }
1260
1261         hull = isl_set_affine_hull(context);
1262         return isl_aff_substitute_equalities_lifted(aff, hull);
1263 error:
1264         isl_aff_free(aff);
1265         isl_set_free(context);
1266         return NULL;
1267 }
1268
1269 __isl_give isl_aff *isl_aff_gist_params(__isl_take isl_aff *aff,
1270         __isl_take isl_set *context)
1271 {
1272         isl_set *dom_context = isl_set_universe(isl_aff_get_domain_space(aff));
1273         dom_context = isl_set_intersect_params(dom_context, context);
1274         return isl_aff_gist(aff, dom_context);
1275 }
1276
1277 /* Return a basic set containing those elements in the space
1278  * of aff where it is non-negative.
1279  */
1280 __isl_give isl_basic_set *isl_aff_nonneg_basic_set(__isl_take isl_aff *aff)
1281 {
1282         isl_constraint *ineq;
1283         isl_basic_set *bset;
1284
1285         ineq = isl_inequality_from_aff(aff);
1286
1287         bset = isl_basic_set_from_constraint(ineq);
1288         bset = isl_basic_set_simplify(bset);
1289         return bset;
1290 }
1291
1292 /* Return a basic set containing those elements in the domain space
1293  * of aff where it is negative.
1294  */
1295 __isl_give isl_basic_set *isl_aff_neg_basic_set(__isl_take isl_aff *aff)
1296 {
1297         aff = isl_aff_neg(aff);
1298         aff = isl_aff_add_constant_num_si(aff, -1);
1299         return isl_aff_nonneg_basic_set(aff);
1300 }
1301
1302 /* Return a basic set containing those elements in the space
1303  * of aff where it is zero.
1304  */
1305 __isl_give isl_basic_set *isl_aff_zero_basic_set(__isl_take isl_aff *aff)
1306 {
1307         isl_constraint *ineq;
1308         isl_basic_set *bset;
1309
1310         ineq = isl_equality_from_aff(aff);
1311
1312         bset = isl_basic_set_from_constraint(ineq);
1313         bset = isl_basic_set_simplify(bset);
1314         return bset;
1315 }
1316
1317 /* Return a basic set containing those elements in the shared space
1318  * of aff1 and aff2 where aff1 is greater than or equal to aff2.
1319  */
1320 __isl_give isl_basic_set *isl_aff_ge_basic_set(__isl_take isl_aff *aff1,
1321         __isl_take isl_aff *aff2)
1322 {
1323         aff1 = isl_aff_sub(aff1, aff2);
1324
1325         return isl_aff_nonneg_basic_set(aff1);
1326 }
1327
1328 /* Return a basic set containing those elements in the shared space
1329  * of aff1 and aff2 where aff1 is smaller than or equal to aff2.
1330  */
1331 __isl_give isl_basic_set *isl_aff_le_basic_set(__isl_take isl_aff *aff1,
1332         __isl_take isl_aff *aff2)
1333 {
1334         return isl_aff_ge_basic_set(aff2, aff1);
1335 }
1336
1337 __isl_give isl_aff *isl_aff_add_on_domain(__isl_keep isl_set *dom,
1338         __isl_take isl_aff *aff1, __isl_take isl_aff *aff2)
1339 {
1340         aff1 = isl_aff_add(aff1, aff2);
1341         aff1 = isl_aff_gist(aff1, isl_set_copy(dom));
1342         return aff1;
1343 }
1344
1345 int isl_aff_is_empty(__isl_keep isl_aff *aff)
1346 {
1347         if (!aff)
1348                 return -1;
1349
1350         return 0;
1351 }
1352
1353 /* Check whether the given affine expression has non-zero coefficient
1354  * for any dimension in the given range or if any of these dimensions
1355  * appear with non-zero coefficients in any of the integer divisions
1356  * involved in the affine expression.
1357  */
1358 int isl_aff_involves_dims(__isl_keep isl_aff *aff,
1359         enum isl_dim_type type, unsigned first, unsigned n)
1360 {
1361         int i;
1362         isl_ctx *ctx;
1363         int *active = NULL;
1364         int involves = 0;
1365
1366         if (!aff)
1367                 return -1;
1368         if (n == 0)
1369                 return 0;
1370
1371         ctx = isl_aff_get_ctx(aff);
1372         if (first + n > isl_aff_dim(aff, type))
1373                 isl_die(ctx, isl_error_invalid,
1374                         "range out of bounds", return -1);
1375
1376         active = isl_local_space_get_active(aff->ls, aff->v->el + 2);
1377         if (!active)
1378                 goto error;
1379
1380         first += isl_local_space_offset(aff->ls, type) - 1;
1381         for (i = 0; i < n; ++i)
1382                 if (active[first + i]) {
1383                         involves = 1;
1384                         break;
1385                 }
1386
1387         free(active);
1388
1389         return involves;
1390 error:
1391         free(active);
1392         return -1;
1393 }
1394
1395 __isl_give isl_aff *isl_aff_drop_dims(__isl_take isl_aff *aff,
1396         enum isl_dim_type type, unsigned first, unsigned n)
1397 {
1398         isl_ctx *ctx;
1399
1400         if (!aff)
1401                 return NULL;
1402         if (type == isl_dim_out)
1403                 isl_die(aff->v->ctx, isl_error_invalid,
1404                         "cannot drop output/set dimension",
1405                         return isl_aff_free(aff));
1406         if (type == isl_dim_in)
1407                 type = isl_dim_set;
1408         if (n == 0 && !isl_local_space_is_named_or_nested(aff->ls, type))
1409                 return aff;
1410
1411         ctx = isl_aff_get_ctx(aff);
1412         if (first + n > isl_local_space_dim(aff->ls, type))
1413                 isl_die(ctx, isl_error_invalid, "range out of bounds",
1414                         return isl_aff_free(aff));
1415
1416         aff = isl_aff_cow(aff);
1417         if (!aff)
1418                 return NULL;
1419
1420         aff->ls = isl_local_space_drop_dims(aff->ls, type, first, n);
1421         if (!aff->ls)
1422                 return isl_aff_free(aff);
1423
1424         first += 1 + isl_local_space_offset(aff->ls, type);
1425         aff->v = isl_vec_drop_els(aff->v, first, n);
1426         if (!aff->v)
1427                 return isl_aff_free(aff);
1428
1429         return aff;
1430 }
1431
1432 /* Project the domain of the affine expression onto its parameter space.
1433  * The affine expression may not involve any of the domain dimensions.
1434  */
1435 __isl_give isl_aff *isl_aff_project_domain_on_params(__isl_take isl_aff *aff)
1436 {
1437         isl_space *space;
1438         unsigned n;
1439         int involves;
1440
1441         n = isl_aff_dim(aff, isl_dim_in);
1442         involves = isl_aff_involves_dims(aff, isl_dim_in, 0, n);
1443         if (involves < 0)
1444                 return isl_aff_free(aff);
1445         if (involves)
1446                 isl_die(isl_aff_get_ctx(aff), isl_error_invalid,
1447                     "affine expression involves some of the domain dimensions",
1448                     return isl_aff_free(aff));
1449         aff = isl_aff_drop_dims(aff, isl_dim_in, 0, n);
1450         space = isl_aff_get_domain_space(aff);
1451         space = isl_space_params(space);
1452         aff = isl_aff_reset_domain_space(aff, space);
1453         return aff;
1454 }
1455
1456 __isl_give isl_aff *isl_aff_insert_dims(__isl_take isl_aff *aff,
1457         enum isl_dim_type type, unsigned first, unsigned n)
1458 {
1459         isl_ctx *ctx;
1460
1461         if (!aff)
1462                 return NULL;
1463         if (type == isl_dim_out)
1464                 isl_die(aff->v->ctx, isl_error_invalid,
1465                         "cannot insert output/set dimensions",
1466                         return isl_aff_free(aff));
1467         if (type == isl_dim_in)
1468                 type = isl_dim_set;
1469         if (n == 0 && !isl_local_space_is_named_or_nested(aff->ls, type))
1470                 return aff;
1471
1472         ctx = isl_aff_get_ctx(aff);
1473         if (first > isl_local_space_dim(aff->ls, type))
1474                 isl_die(ctx, isl_error_invalid, "position out of bounds",
1475                         return isl_aff_free(aff));
1476
1477         aff = isl_aff_cow(aff);
1478         if (!aff)
1479                 return NULL;
1480
1481         aff->ls = isl_local_space_insert_dims(aff->ls, type, first, n);
1482         if (!aff->ls)
1483                 return isl_aff_free(aff);
1484
1485         first += 1 + isl_local_space_offset(aff->ls, type);
1486         aff->v = isl_vec_insert_zero_els(aff->v, first, n);
1487         if (!aff->v)
1488                 return isl_aff_free(aff);
1489
1490         return aff;
1491 }
1492
1493 __isl_give isl_aff *isl_aff_add_dims(__isl_take isl_aff *aff,
1494         enum isl_dim_type type, unsigned n)
1495 {
1496         unsigned pos;
1497
1498         pos = isl_aff_dim(aff, type);
1499
1500         return isl_aff_insert_dims(aff, type, pos, n);
1501 }
1502
1503 __isl_give isl_pw_aff *isl_pw_aff_add_dims(__isl_take isl_pw_aff *pwaff,
1504         enum isl_dim_type type, unsigned n)
1505 {
1506         unsigned pos;
1507
1508         pos = isl_pw_aff_dim(pwaff, type);
1509
1510         return isl_pw_aff_insert_dims(pwaff, type, pos, n);
1511 }
1512
1513 __isl_give isl_pw_aff *isl_pw_aff_from_aff(__isl_take isl_aff *aff)
1514 {
1515         isl_set *dom = isl_set_universe(isl_aff_get_domain_space(aff));
1516         return isl_pw_aff_alloc(dom, aff);
1517 }
1518
1519 #undef PW
1520 #define PW isl_pw_aff
1521 #undef EL
1522 #define EL isl_aff
1523 #undef EL_IS_ZERO
1524 #define EL_IS_ZERO is_empty
1525 #undef ZERO
1526 #define ZERO empty
1527 #undef IS_ZERO
1528 #define IS_ZERO is_empty
1529 #undef FIELD
1530 #define FIELD aff
1531 #undef DEFAULT_IS_ZERO
1532 #define DEFAULT_IS_ZERO 0
1533
1534 #define NO_EVAL
1535 #define NO_OPT
1536 #define NO_MOVE_DIMS
1537 #define NO_LIFT
1538 #define NO_MORPH
1539
1540 #include <isl_pw_templ.c>
1541
1542 static __isl_give isl_set *align_params_pw_pw_set_and(
1543         __isl_take isl_pw_aff *pwaff1, __isl_take isl_pw_aff *pwaff2,
1544         __isl_give isl_set *(*fn)(__isl_take isl_pw_aff *pwaff1,
1545                                     __isl_take isl_pw_aff *pwaff2))
1546 {
1547         if (!pwaff1 || !pwaff2)
1548                 goto error;
1549         if (isl_space_match(pwaff1->dim, isl_dim_param,
1550                           pwaff2->dim, isl_dim_param))
1551                 return fn(pwaff1, pwaff2);
1552         if (!isl_space_has_named_params(pwaff1->dim) ||
1553             !isl_space_has_named_params(pwaff2->dim))
1554                 isl_die(isl_pw_aff_get_ctx(pwaff1), isl_error_invalid,
1555                         "unaligned unnamed parameters", goto error);
1556         pwaff1 = isl_pw_aff_align_params(pwaff1, isl_pw_aff_get_space(pwaff2));
1557         pwaff2 = isl_pw_aff_align_params(pwaff2, isl_pw_aff_get_space(pwaff1));
1558         return fn(pwaff1, pwaff2);
1559 error:
1560         isl_pw_aff_free(pwaff1);
1561         isl_pw_aff_free(pwaff2);
1562         return NULL;
1563 }
1564
1565 /* Compute a piecewise quasi-affine expression with a domain that
1566  * is the union of those of pwaff1 and pwaff2 and such that on each
1567  * cell, the quasi-affine expression is the better (according to cmp)
1568  * of those of pwaff1 and pwaff2.  If only one of pwaff1 or pwaff2
1569  * is defined on a given cell, then the associated expression
1570  * is the defined one.
1571  */
1572 static __isl_give isl_pw_aff *pw_aff_union_opt(__isl_take isl_pw_aff *pwaff1,
1573         __isl_take isl_pw_aff *pwaff2,
1574         __isl_give isl_basic_set *(*cmp)(__isl_take isl_aff *aff1,
1575                                         __isl_take isl_aff *aff2))
1576 {
1577         int i, j, n;
1578         isl_pw_aff *res;
1579         isl_ctx *ctx;
1580         isl_set *set;
1581
1582         if (!pwaff1 || !pwaff2)
1583                 goto error;
1584
1585         ctx = isl_space_get_ctx(pwaff1->dim);
1586         if (!isl_space_is_equal(pwaff1->dim, pwaff2->dim))
1587                 isl_die(ctx, isl_error_invalid,
1588                         "arguments should live in same space", goto error);
1589
1590         if (isl_pw_aff_is_empty(pwaff1)) {
1591                 isl_pw_aff_free(pwaff1);
1592                 return pwaff2;
1593         }
1594
1595         if (isl_pw_aff_is_empty(pwaff2)) {
1596                 isl_pw_aff_free(pwaff2);
1597                 return pwaff1;
1598         }
1599
1600         n = 2 * (pwaff1->n + 1) * (pwaff2->n + 1);
1601         res = isl_pw_aff_alloc_size(isl_space_copy(pwaff1->dim), n);
1602
1603         for (i = 0; i < pwaff1->n; ++i) {
1604                 set = isl_set_copy(pwaff1->p[i].set);
1605                 for (j = 0; j < pwaff2->n; ++j) {
1606                         struct isl_set *common;
1607                         isl_set *better;
1608
1609                         common = isl_set_intersect(
1610                                         isl_set_copy(pwaff1->p[i].set),
1611                                         isl_set_copy(pwaff2->p[j].set));
1612                         better = isl_set_from_basic_set(cmp(
1613                                         isl_aff_copy(pwaff2->p[j].aff),
1614                                         isl_aff_copy(pwaff1->p[i].aff)));
1615                         better = isl_set_intersect(common, better);
1616                         if (isl_set_plain_is_empty(better)) {
1617                                 isl_set_free(better);
1618                                 continue;
1619                         }
1620                         set = isl_set_subtract(set, isl_set_copy(better));
1621
1622                         res = isl_pw_aff_add_piece(res, better,
1623                                                 isl_aff_copy(pwaff2->p[j].aff));
1624                 }
1625                 res = isl_pw_aff_add_piece(res, set,
1626                                                 isl_aff_copy(pwaff1->p[i].aff));
1627         }
1628
1629         for (j = 0; j < pwaff2->n; ++j) {
1630                 set = isl_set_copy(pwaff2->p[j].set);
1631                 for (i = 0; i < pwaff1->n; ++i)
1632                         set = isl_set_subtract(set,
1633                                         isl_set_copy(pwaff1->p[i].set));
1634                 res = isl_pw_aff_add_piece(res, set,
1635                                                 isl_aff_copy(pwaff2->p[j].aff));
1636         }
1637
1638         isl_pw_aff_free(pwaff1);
1639         isl_pw_aff_free(pwaff2);
1640
1641         return res;
1642 error:
1643         isl_pw_aff_free(pwaff1);
1644         isl_pw_aff_free(pwaff2);
1645         return NULL;
1646 }
1647
1648 /* Compute a piecewise quasi-affine expression with a domain that
1649  * is the union of those of pwaff1 and pwaff2 and such that on each
1650  * cell, the quasi-affine expression is the maximum of those of pwaff1
1651  * and pwaff2.  If only one of pwaff1 or pwaff2 is defined on a given
1652  * cell, then the associated expression is the defined one.
1653  */
1654 static __isl_give isl_pw_aff *pw_aff_union_max(__isl_take isl_pw_aff *pwaff1,
1655         __isl_take isl_pw_aff *pwaff2)
1656 {
1657         return pw_aff_union_opt(pwaff1, pwaff2, &isl_aff_ge_basic_set);
1658 }
1659
1660 __isl_give isl_pw_aff *isl_pw_aff_union_max(__isl_take isl_pw_aff *pwaff1,
1661         __isl_take isl_pw_aff *pwaff2)
1662 {
1663         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2,
1664                                                         &pw_aff_union_max);
1665 }
1666
1667 /* Compute a piecewise quasi-affine expression with a domain that
1668  * is the union of those of pwaff1 and pwaff2 and such that on each
1669  * cell, the quasi-affine expression is the minimum of those of pwaff1
1670  * and pwaff2.  If only one of pwaff1 or pwaff2 is defined on a given
1671  * cell, then the associated expression is the defined one.
1672  */
1673 static __isl_give isl_pw_aff *pw_aff_union_min(__isl_take isl_pw_aff *pwaff1,
1674         __isl_take isl_pw_aff *pwaff2)
1675 {
1676         return pw_aff_union_opt(pwaff1, pwaff2, &isl_aff_le_basic_set);
1677 }
1678
1679 __isl_give isl_pw_aff *isl_pw_aff_union_min(__isl_take isl_pw_aff *pwaff1,
1680         __isl_take isl_pw_aff *pwaff2)
1681 {
1682         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2,
1683                                                         &pw_aff_union_min);
1684 }
1685
1686 __isl_give isl_pw_aff *isl_pw_aff_union_opt(__isl_take isl_pw_aff *pwaff1,
1687         __isl_take isl_pw_aff *pwaff2, int max)
1688 {
1689         if (max)
1690                 return isl_pw_aff_union_max(pwaff1, pwaff2);
1691         else
1692                 return isl_pw_aff_union_min(pwaff1, pwaff2);
1693 }
1694
1695 /* Construct a map with as domain the domain of pwaff and
1696  * one-dimensional range corresponding to the affine expressions.
1697  */
1698 static __isl_give isl_map *map_from_pw_aff(__isl_take isl_pw_aff *pwaff)
1699 {
1700         int i;
1701         isl_space *dim;
1702         isl_map *map;
1703
1704         if (!pwaff)
1705                 return NULL;
1706
1707         dim = isl_pw_aff_get_space(pwaff);
1708         map = isl_map_empty(dim);
1709
1710         for (i = 0; i < pwaff->n; ++i) {
1711                 isl_basic_map *bmap;
1712                 isl_map *map_i;
1713
1714                 bmap = isl_basic_map_from_aff(isl_aff_copy(pwaff->p[i].aff));
1715                 map_i = isl_map_from_basic_map(bmap);
1716                 map_i = isl_map_intersect_domain(map_i,
1717                                                 isl_set_copy(pwaff->p[i].set));
1718                 map = isl_map_union_disjoint(map, map_i);
1719         }
1720
1721         isl_pw_aff_free(pwaff);
1722
1723         return map;
1724 }
1725
1726 /* Construct a map with as domain the domain of pwaff and
1727  * one-dimensional range corresponding to the affine expressions.
1728  */
1729 __isl_give isl_map *isl_map_from_pw_aff(__isl_take isl_pw_aff *pwaff)
1730 {
1731         if (!pwaff)
1732                 return NULL;
1733         if (isl_space_is_set(pwaff->dim))
1734                 isl_die(isl_pw_aff_get_ctx(pwaff), isl_error_invalid,
1735                         "space of input is not a map",
1736                         return isl_pw_aff_free(pwaff));
1737         return map_from_pw_aff(pwaff);
1738 }
1739
1740 /* Construct a one-dimensional set with as parameter domain
1741  * the domain of pwaff and the single set dimension
1742  * corresponding to the affine expressions.
1743  */
1744 __isl_give isl_set *isl_set_from_pw_aff(__isl_take isl_pw_aff *pwaff)
1745 {
1746         if (!pwaff)
1747                 return NULL;
1748         if (!isl_space_is_set(pwaff->dim))
1749                 isl_die(isl_pw_aff_get_ctx(pwaff), isl_error_invalid,
1750                         "space of input is not a set",
1751                         return isl_pw_aff_free(pwaff));
1752         return map_from_pw_aff(pwaff);
1753 }
1754
1755 /* Return a set containing those elements in the domain
1756  * of pwaff where it is non-negative.
1757  */
1758 __isl_give isl_set *isl_pw_aff_nonneg_set(__isl_take isl_pw_aff *pwaff)
1759 {
1760         int i;
1761         isl_set *set;
1762
1763         if (!pwaff)
1764                 return NULL;
1765
1766         set = isl_set_empty(isl_pw_aff_get_domain_space(pwaff));
1767
1768         for (i = 0; i < pwaff->n; ++i) {
1769                 isl_basic_set *bset;
1770                 isl_set *set_i;
1771
1772                 bset = isl_aff_nonneg_basic_set(isl_aff_copy(pwaff->p[i].aff));
1773                 set_i = isl_set_from_basic_set(bset);
1774                 set_i = isl_set_intersect(set_i, isl_set_copy(pwaff->p[i].set));
1775                 set = isl_set_union_disjoint(set, set_i);
1776         }
1777
1778         isl_pw_aff_free(pwaff);
1779
1780         return set;
1781 }
1782
1783 /* Return a set containing those elements in the domain
1784  * of pwaff where it is zero (if complement is 0) or not zero
1785  * (if complement is 1).
1786  */
1787 static __isl_give isl_set *pw_aff_zero_set(__isl_take isl_pw_aff *pwaff,
1788         int complement)
1789 {
1790         int i;
1791         isl_set *set;
1792
1793         if (!pwaff)
1794                 return NULL;
1795
1796         set = isl_set_empty(isl_pw_aff_get_domain_space(pwaff));
1797
1798         for (i = 0; i < pwaff->n; ++i) {
1799                 isl_basic_set *bset;
1800                 isl_set *set_i, *zero;
1801
1802                 bset = isl_aff_zero_basic_set(isl_aff_copy(pwaff->p[i].aff));
1803                 zero = isl_set_from_basic_set(bset);
1804                 set_i = isl_set_copy(pwaff->p[i].set);
1805                 if (complement)
1806                         set_i = isl_set_subtract(set_i, zero);
1807                 else
1808                         set_i = isl_set_intersect(set_i, zero);
1809                 set = isl_set_union_disjoint(set, set_i);
1810         }
1811
1812         isl_pw_aff_free(pwaff);
1813
1814         return set;
1815 }
1816
1817 /* Return a set containing those elements in the domain
1818  * of pwaff where it is zero.
1819  */
1820 __isl_give isl_set *isl_pw_aff_zero_set(__isl_take isl_pw_aff *pwaff)
1821 {
1822         return pw_aff_zero_set(pwaff, 0);
1823 }
1824
1825 /* Return a set containing those elements in the domain
1826  * of pwaff where it is not zero.
1827  */
1828 __isl_give isl_set *isl_pw_aff_non_zero_set(__isl_take isl_pw_aff *pwaff)
1829 {
1830         return pw_aff_zero_set(pwaff, 1);
1831 }
1832
1833 /* Return a set containing those elements in the shared domain
1834  * of pwaff1 and pwaff2 where pwaff1 is greater than (or equal) to pwaff2.
1835  *
1836  * We compute the difference on the shared domain and then construct
1837  * the set of values where this difference is non-negative.
1838  * If strict is set, we first subtract 1 from the difference.
1839  * If equal is set, we only return the elements where pwaff1 and pwaff2
1840  * are equal.
1841  */
1842 static __isl_give isl_set *pw_aff_gte_set(__isl_take isl_pw_aff *pwaff1,
1843         __isl_take isl_pw_aff *pwaff2, int strict, int equal)
1844 {
1845         isl_set *set1, *set2;
1846
1847         set1 = isl_pw_aff_domain(isl_pw_aff_copy(pwaff1));
1848         set2 = isl_pw_aff_domain(isl_pw_aff_copy(pwaff2));
1849         set1 = isl_set_intersect(set1, set2);
1850         pwaff1 = isl_pw_aff_intersect_domain(pwaff1, isl_set_copy(set1));
1851         pwaff2 = isl_pw_aff_intersect_domain(pwaff2, isl_set_copy(set1));
1852         pwaff1 = isl_pw_aff_add(pwaff1, isl_pw_aff_neg(pwaff2));
1853
1854         if (strict) {
1855                 isl_space *dim = isl_set_get_space(set1);
1856                 isl_aff *aff;
1857                 aff = isl_aff_zero_on_domain(isl_local_space_from_space(dim));
1858                 aff = isl_aff_add_constant_si(aff, -1);
1859                 pwaff1 = isl_pw_aff_add(pwaff1, isl_pw_aff_alloc(set1, aff));
1860         } else
1861                 isl_set_free(set1);
1862
1863         if (equal)
1864                 return isl_pw_aff_zero_set(pwaff1);
1865         return isl_pw_aff_nonneg_set(pwaff1);
1866 }
1867
1868 /* Return a set containing those elements in the shared domain
1869  * of pwaff1 and pwaff2 where pwaff1 is equal to pwaff2.
1870  */
1871 static __isl_give isl_set *pw_aff_eq_set(__isl_take isl_pw_aff *pwaff1,
1872         __isl_take isl_pw_aff *pwaff2)
1873 {
1874         return pw_aff_gte_set(pwaff1, pwaff2, 0, 1);
1875 }
1876
1877 __isl_give isl_set *isl_pw_aff_eq_set(__isl_take isl_pw_aff *pwaff1,
1878         __isl_take isl_pw_aff *pwaff2)
1879 {
1880         return align_params_pw_pw_set_and(pwaff1, pwaff2, &pw_aff_eq_set);
1881 }
1882
1883 /* Return a set containing those elements in the shared domain
1884  * of pwaff1 and pwaff2 where pwaff1 is greater than or equal to pwaff2.
1885  */
1886 static __isl_give isl_set *pw_aff_ge_set(__isl_take isl_pw_aff *pwaff1,
1887         __isl_take isl_pw_aff *pwaff2)
1888 {
1889         return pw_aff_gte_set(pwaff1, pwaff2, 0, 0);
1890 }
1891
1892 __isl_give isl_set *isl_pw_aff_ge_set(__isl_take isl_pw_aff *pwaff1,
1893         __isl_take isl_pw_aff *pwaff2)
1894 {
1895         return align_params_pw_pw_set_and(pwaff1, pwaff2, &pw_aff_ge_set);
1896 }
1897
1898 /* Return a set containing those elements in the shared domain
1899  * of pwaff1 and pwaff2 where pwaff1 is strictly greater than pwaff2.
1900  */
1901 static __isl_give isl_set *pw_aff_gt_set(__isl_take isl_pw_aff *pwaff1,
1902         __isl_take isl_pw_aff *pwaff2)
1903 {
1904         return pw_aff_gte_set(pwaff1, pwaff2, 1, 0);
1905 }
1906
1907 __isl_give isl_set *isl_pw_aff_gt_set(__isl_take isl_pw_aff *pwaff1,
1908         __isl_take isl_pw_aff *pwaff2)
1909 {
1910         return align_params_pw_pw_set_and(pwaff1, pwaff2, &pw_aff_gt_set);
1911 }
1912
1913 __isl_give isl_set *isl_pw_aff_le_set(__isl_take isl_pw_aff *pwaff1,
1914         __isl_take isl_pw_aff *pwaff2)
1915 {
1916         return isl_pw_aff_ge_set(pwaff2, pwaff1);
1917 }
1918
1919 __isl_give isl_set *isl_pw_aff_lt_set(__isl_take isl_pw_aff *pwaff1,
1920         __isl_take isl_pw_aff *pwaff2)
1921 {
1922         return isl_pw_aff_gt_set(pwaff2, pwaff1);
1923 }
1924
1925 /* Return a set containing those elements in the shared domain
1926  * of the elements of list1 and list2 where each element in list1
1927  * has the relation specified by "fn" with each element in list2.
1928  */
1929 static __isl_give isl_set *pw_aff_list_set(__isl_take isl_pw_aff_list *list1,
1930         __isl_take isl_pw_aff_list *list2,
1931         __isl_give isl_set *(*fn)(__isl_take isl_pw_aff *pwaff1,
1932                                     __isl_take isl_pw_aff *pwaff2))
1933 {
1934         int i, j;
1935         isl_ctx *ctx;
1936         isl_set *set;
1937
1938         if (!list1 || !list2)
1939                 goto error;
1940
1941         ctx = isl_pw_aff_list_get_ctx(list1);
1942         if (list1->n < 1 || list2->n < 1)
1943                 isl_die(ctx, isl_error_invalid,
1944                         "list should contain at least one element", goto error);
1945
1946         set = isl_set_universe(isl_pw_aff_get_domain_space(list1->p[0]));
1947         for (i = 0; i < list1->n; ++i)
1948                 for (j = 0; j < list2->n; ++j) {
1949                         isl_set *set_ij;
1950
1951                         set_ij = fn(isl_pw_aff_copy(list1->p[i]),
1952                                     isl_pw_aff_copy(list2->p[j]));
1953                         set = isl_set_intersect(set, set_ij);
1954                 }
1955
1956         isl_pw_aff_list_free(list1);
1957         isl_pw_aff_list_free(list2);
1958         return set;
1959 error:
1960         isl_pw_aff_list_free(list1);
1961         isl_pw_aff_list_free(list2);
1962         return NULL;
1963 }
1964
1965 /* Return a set containing those elements in the shared domain
1966  * of the elements of list1 and list2 where each element in list1
1967  * is equal to each element in list2.
1968  */
1969 __isl_give isl_set *isl_pw_aff_list_eq_set(__isl_take isl_pw_aff_list *list1,
1970         __isl_take isl_pw_aff_list *list2)
1971 {
1972         return pw_aff_list_set(list1, list2, &isl_pw_aff_eq_set);
1973 }
1974
1975 __isl_give isl_set *isl_pw_aff_list_ne_set(__isl_take isl_pw_aff_list *list1,
1976         __isl_take isl_pw_aff_list *list2)
1977 {
1978         return pw_aff_list_set(list1, list2, &isl_pw_aff_ne_set);
1979 }
1980
1981 /* Return a set containing those elements in the shared domain
1982  * of the elements of list1 and list2 where each element in list1
1983  * is less than or equal to each element in list2.
1984  */
1985 __isl_give isl_set *isl_pw_aff_list_le_set(__isl_take isl_pw_aff_list *list1,
1986         __isl_take isl_pw_aff_list *list2)
1987 {
1988         return pw_aff_list_set(list1, list2, &isl_pw_aff_le_set);
1989 }
1990
1991 __isl_give isl_set *isl_pw_aff_list_lt_set(__isl_take isl_pw_aff_list *list1,
1992         __isl_take isl_pw_aff_list *list2)
1993 {
1994         return pw_aff_list_set(list1, list2, &isl_pw_aff_lt_set);
1995 }
1996
1997 __isl_give isl_set *isl_pw_aff_list_ge_set(__isl_take isl_pw_aff_list *list1,
1998         __isl_take isl_pw_aff_list *list2)
1999 {
2000         return pw_aff_list_set(list1, list2, &isl_pw_aff_ge_set);
2001 }
2002
2003 __isl_give isl_set *isl_pw_aff_list_gt_set(__isl_take isl_pw_aff_list *list1,
2004         __isl_take isl_pw_aff_list *list2)
2005 {
2006         return pw_aff_list_set(list1, list2, &isl_pw_aff_gt_set);
2007 }
2008
2009
2010 /* Return a set containing those elements in the shared domain
2011  * of pwaff1 and pwaff2 where pwaff1 is not equal to pwaff2.
2012  */
2013 static __isl_give isl_set *pw_aff_ne_set(__isl_take isl_pw_aff *pwaff1,
2014         __isl_take isl_pw_aff *pwaff2)
2015 {
2016         isl_set *set_lt, *set_gt;
2017
2018         set_lt = isl_pw_aff_lt_set(isl_pw_aff_copy(pwaff1),
2019                                    isl_pw_aff_copy(pwaff2));
2020         set_gt = isl_pw_aff_gt_set(pwaff1, pwaff2);
2021         return isl_set_union_disjoint(set_lt, set_gt);
2022 }
2023
2024 __isl_give isl_set *isl_pw_aff_ne_set(__isl_take isl_pw_aff *pwaff1,
2025         __isl_take isl_pw_aff *pwaff2)
2026 {
2027         return align_params_pw_pw_set_and(pwaff1, pwaff2, &pw_aff_ne_set);
2028 }
2029
2030 __isl_give isl_pw_aff *isl_pw_aff_scale_down(__isl_take isl_pw_aff *pwaff,
2031         isl_int v)
2032 {
2033         int i;
2034
2035         if (isl_int_is_one(v))
2036                 return pwaff;
2037         if (!isl_int_is_pos(v))
2038                 isl_die(isl_pw_aff_get_ctx(pwaff), isl_error_invalid,
2039                         "factor needs to be positive",
2040                         return isl_pw_aff_free(pwaff));
2041         pwaff = isl_pw_aff_cow(pwaff);
2042         if (!pwaff)
2043                 return NULL;
2044         if (pwaff->n == 0)
2045                 return pwaff;
2046
2047         for (i = 0; i < pwaff->n; ++i) {
2048                 pwaff->p[i].aff = isl_aff_scale_down(pwaff->p[i].aff, v);
2049                 if (!pwaff->p[i].aff)
2050                         return isl_pw_aff_free(pwaff);
2051         }
2052
2053         return pwaff;
2054 }
2055
2056 __isl_give isl_pw_aff *isl_pw_aff_floor(__isl_take isl_pw_aff *pwaff)
2057 {
2058         int i;
2059
2060         pwaff = isl_pw_aff_cow(pwaff);
2061         if (!pwaff)
2062                 return NULL;
2063         if (pwaff->n == 0)
2064                 return pwaff;
2065
2066         for (i = 0; i < pwaff->n; ++i) {
2067                 pwaff->p[i].aff = isl_aff_floor(pwaff->p[i].aff);
2068                 if (!pwaff->p[i].aff)
2069                         return isl_pw_aff_free(pwaff);
2070         }
2071
2072         return pwaff;
2073 }
2074
2075 __isl_give isl_pw_aff *isl_pw_aff_ceil(__isl_take isl_pw_aff *pwaff)
2076 {
2077         int i;
2078
2079         pwaff = isl_pw_aff_cow(pwaff);
2080         if (!pwaff)
2081                 return NULL;
2082         if (pwaff->n == 0)
2083                 return pwaff;
2084
2085         for (i = 0; i < pwaff->n; ++i) {
2086                 pwaff->p[i].aff = isl_aff_ceil(pwaff->p[i].aff);
2087                 if (!pwaff->p[i].aff)
2088                         return isl_pw_aff_free(pwaff);
2089         }
2090
2091         return pwaff;
2092 }
2093
2094 /* Assuming that "cond1" and "cond2" are disjoint,
2095  * return an affine expression that is equal to pwaff1 on cond1
2096  * and to pwaff2 on cond2.
2097  */
2098 static __isl_give isl_pw_aff *isl_pw_aff_select(
2099         __isl_take isl_set *cond1, __isl_take isl_pw_aff *pwaff1,
2100         __isl_take isl_set *cond2, __isl_take isl_pw_aff *pwaff2)
2101 {
2102         pwaff1 = isl_pw_aff_intersect_domain(pwaff1, cond1);
2103         pwaff2 = isl_pw_aff_intersect_domain(pwaff2, cond2);
2104
2105         return isl_pw_aff_add_disjoint(pwaff1, pwaff2);
2106 }
2107
2108 /* Return an affine expression that is equal to pwaff_true for elements
2109  * where "cond" is non-zero and to pwaff_false for elements where "cond"
2110  * is zero.
2111  * That is, return cond ? pwaff_true : pwaff_false;
2112  */
2113 __isl_give isl_pw_aff *isl_pw_aff_cond(__isl_take isl_pw_aff *cond,
2114         __isl_take isl_pw_aff *pwaff_true, __isl_take isl_pw_aff *pwaff_false)
2115 {
2116         isl_set *cond_true, *cond_false;
2117
2118         cond_true = isl_pw_aff_non_zero_set(isl_pw_aff_copy(cond));
2119         cond_false = isl_pw_aff_zero_set(cond);
2120         return isl_pw_aff_select(cond_true, pwaff_true,
2121                                  cond_false, pwaff_false);
2122 }
2123
2124 int isl_aff_is_cst(__isl_keep isl_aff *aff)
2125 {
2126         if (!aff)
2127                 return -1;
2128
2129         return isl_seq_first_non_zero(aff->v->el + 2, aff->v->size - 2) == -1;
2130 }
2131
2132 /* Check whether pwaff is a piecewise constant.
2133  */
2134 int isl_pw_aff_is_cst(__isl_keep isl_pw_aff *pwaff)
2135 {
2136         int i;
2137
2138         if (!pwaff)
2139                 return -1;
2140
2141         for (i = 0; i < pwaff->n; ++i) {
2142                 int is_cst = isl_aff_is_cst(pwaff->p[i].aff);
2143                 if (is_cst < 0 || !is_cst)
2144                         return is_cst;
2145         }
2146
2147         return 1;
2148 }
2149
2150 __isl_give isl_aff *isl_aff_mul(__isl_take isl_aff *aff1,
2151         __isl_take isl_aff *aff2)
2152 {
2153         if (!isl_aff_is_cst(aff2) && isl_aff_is_cst(aff1))
2154                 return isl_aff_mul(aff2, aff1);
2155
2156         if (!isl_aff_is_cst(aff2))
2157                 isl_die(isl_aff_get_ctx(aff1), isl_error_invalid,
2158                         "at least one affine expression should be constant",
2159                         goto error);
2160
2161         aff1 = isl_aff_cow(aff1);
2162         if (!aff1 || !aff2)
2163                 goto error;
2164
2165         aff1 = isl_aff_scale(aff1, aff2->v->el[1]);
2166         aff1 = isl_aff_scale_down(aff1, aff2->v->el[0]);
2167
2168         isl_aff_free(aff2);
2169         return aff1;
2170 error:
2171         isl_aff_free(aff1);
2172         isl_aff_free(aff2);
2173         return NULL;
2174 }
2175
2176 static __isl_give isl_pw_aff *pw_aff_add(__isl_take isl_pw_aff *pwaff1,
2177         __isl_take isl_pw_aff *pwaff2)
2178 {
2179         return isl_pw_aff_on_shared_domain(pwaff1, pwaff2, &isl_aff_add);
2180 }
2181
2182 __isl_give isl_pw_aff *isl_pw_aff_add(__isl_take isl_pw_aff *pwaff1,
2183         __isl_take isl_pw_aff *pwaff2)
2184 {
2185         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2, &pw_aff_add);
2186 }
2187
2188 __isl_give isl_pw_aff *isl_pw_aff_union_add(__isl_take isl_pw_aff *pwaff1,
2189         __isl_take isl_pw_aff *pwaff2)
2190 {
2191         return isl_pw_aff_union_add_(pwaff1, pwaff2);
2192 }
2193
2194 static __isl_give isl_pw_aff *pw_aff_mul(__isl_take isl_pw_aff *pwaff1,
2195         __isl_take isl_pw_aff *pwaff2)
2196 {
2197         return isl_pw_aff_on_shared_domain(pwaff1, pwaff2, &isl_aff_mul);
2198 }
2199
2200 __isl_give isl_pw_aff *isl_pw_aff_mul(__isl_take isl_pw_aff *pwaff1,
2201         __isl_take isl_pw_aff *pwaff2)
2202 {
2203         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2, &pw_aff_mul);
2204 }
2205
2206 static __isl_give isl_pw_aff *pw_aff_min(__isl_take isl_pw_aff *pwaff1,
2207         __isl_take isl_pw_aff *pwaff2)
2208 {
2209         isl_set *le;
2210         isl_set *dom;
2211
2212         dom = isl_set_intersect(isl_pw_aff_domain(isl_pw_aff_copy(pwaff1)),
2213                                 isl_pw_aff_domain(isl_pw_aff_copy(pwaff2)));
2214         le = isl_pw_aff_le_set(isl_pw_aff_copy(pwaff1),
2215                                 isl_pw_aff_copy(pwaff2));
2216         dom = isl_set_subtract(dom, isl_set_copy(le));
2217         return isl_pw_aff_select(le, pwaff1, dom, pwaff2);
2218 }
2219
2220 __isl_give isl_pw_aff *isl_pw_aff_min(__isl_take isl_pw_aff *pwaff1,
2221         __isl_take isl_pw_aff *pwaff2)
2222 {
2223         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2, &pw_aff_min);
2224 }
2225
2226 static __isl_give isl_pw_aff *pw_aff_max(__isl_take isl_pw_aff *pwaff1,
2227         __isl_take isl_pw_aff *pwaff2)
2228 {
2229         isl_set *ge;
2230         isl_set *dom;
2231
2232         dom = isl_set_intersect(isl_pw_aff_domain(isl_pw_aff_copy(pwaff1)),
2233                                 isl_pw_aff_domain(isl_pw_aff_copy(pwaff2)));
2234         ge = isl_pw_aff_ge_set(isl_pw_aff_copy(pwaff1),
2235                                 isl_pw_aff_copy(pwaff2));
2236         dom = isl_set_subtract(dom, isl_set_copy(ge));
2237         return isl_pw_aff_select(ge, pwaff1, dom, pwaff2);
2238 }
2239
2240 __isl_give isl_pw_aff *isl_pw_aff_max(__isl_take isl_pw_aff *pwaff1,
2241         __isl_take isl_pw_aff *pwaff2)
2242 {
2243         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2, &pw_aff_max);
2244 }
2245
2246 static __isl_give isl_pw_aff *pw_aff_list_reduce(
2247         __isl_take isl_pw_aff_list *list,
2248         __isl_give isl_pw_aff *(*fn)(__isl_take isl_pw_aff *pwaff1,
2249                                         __isl_take isl_pw_aff *pwaff2))
2250 {
2251         int i;
2252         isl_ctx *ctx;
2253         isl_pw_aff *res;
2254
2255         if (!list)
2256                 return NULL;
2257
2258         ctx = isl_pw_aff_list_get_ctx(list);
2259         if (list->n < 1)
2260                 isl_die(ctx, isl_error_invalid,
2261                         "list should contain at least one element",
2262                         return isl_pw_aff_list_free(list));
2263
2264         res = isl_pw_aff_copy(list->p[0]);
2265         for (i = 1; i < list->n; ++i)
2266                 res = fn(res, isl_pw_aff_copy(list->p[i]));
2267
2268         isl_pw_aff_list_free(list);
2269         return res;
2270 }
2271
2272 /* Return an isl_pw_aff that maps each element in the intersection of the
2273  * domains of the elements of list to the minimal corresponding affine
2274  * expression.
2275  */
2276 __isl_give isl_pw_aff *isl_pw_aff_list_min(__isl_take isl_pw_aff_list *list)
2277 {
2278         return pw_aff_list_reduce(list, &isl_pw_aff_min);
2279 }
2280
2281 /* Return an isl_pw_aff that maps each element in the intersection of the
2282  * domains of the elements of list to the maximal corresponding affine
2283  * expression.
2284  */
2285 __isl_give isl_pw_aff *isl_pw_aff_list_max(__isl_take isl_pw_aff_list *list)
2286 {
2287         return pw_aff_list_reduce(list, &isl_pw_aff_max);
2288 }
2289
2290 #undef BASE
2291 #define BASE aff
2292
2293 #include <isl_multi_templ.c>
2294
2295 /* Construct an isl_multi_aff in the given space with value zero in
2296  * each of the output dimensions.
2297  */
2298 __isl_give isl_multi_aff *isl_multi_aff_zero(__isl_take isl_space *space)
2299 {
2300         int n;
2301         isl_multi_aff *ma;
2302
2303         if (!space)
2304                 return NULL;
2305
2306         n = isl_space_dim(space , isl_dim_out);
2307         ma = isl_multi_aff_alloc(isl_space_copy(space));
2308
2309         if (!n)
2310                 isl_space_free(space);
2311         else {
2312                 int i;
2313                 isl_local_space *ls;
2314                 isl_aff *aff;
2315
2316                 space = isl_space_domain(space);
2317                 ls = isl_local_space_from_space(space);
2318                 aff = isl_aff_zero_on_domain(ls);
2319
2320                 for (i = 0; i < n; ++i)
2321                         ma = isl_multi_aff_set_aff(ma, i, isl_aff_copy(aff));
2322
2323                 isl_aff_free(aff);
2324         }
2325
2326         return ma;
2327 }
2328
2329 /* Create an isl_multi_aff in the given space that maps each
2330  * input dimension to the corresponding output dimension.
2331  */
2332 __isl_give isl_multi_aff *isl_multi_aff_identity(__isl_take isl_space *space)
2333 {
2334         int n;
2335         isl_multi_aff *ma;
2336
2337         if (!space)
2338                 return NULL;
2339
2340         if (isl_space_is_set(space))
2341                 isl_die(isl_space_get_ctx(space), isl_error_invalid,
2342                         "expecting map space", goto error);
2343
2344         n = isl_space_dim(space, isl_dim_out);
2345         if (n != isl_space_dim(space, isl_dim_in))
2346                 isl_die(isl_space_get_ctx(space), isl_error_invalid,
2347                         "number of input and output dimensions needs to be "
2348                         "the same", goto error);
2349
2350         ma = isl_multi_aff_alloc(isl_space_copy(space));
2351
2352         if (!n)
2353                 isl_space_free(space);
2354         else {
2355                 int i;
2356                 isl_local_space *ls;
2357                 isl_aff *aff;
2358
2359                 space = isl_space_domain(space);
2360                 ls = isl_local_space_from_space(space);
2361                 aff = isl_aff_zero_on_domain(ls);
2362
2363                 for (i = 0; i < n; ++i) {
2364                         isl_aff *aff_i;
2365                         aff_i = isl_aff_copy(aff);
2366                         aff_i = isl_aff_add_coefficient_si(aff_i,
2367                                                             isl_dim_in, i, 1);
2368                         ma = isl_multi_aff_set_aff(ma, i, aff_i);
2369                 }
2370
2371                 isl_aff_free(aff);
2372         }
2373
2374         return ma;
2375 error:
2376         isl_space_free(space);
2377         return NULL;
2378 }
2379
2380 /* Create an isl_pw_multi_aff with the given isl_multi_aff on a universe
2381  * domain.
2382  */
2383 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_from_multi_aff(
2384         __isl_take isl_multi_aff *ma)
2385 {
2386         isl_set *dom = isl_set_universe(isl_multi_aff_get_domain_space(ma));
2387         return isl_pw_multi_aff_alloc(dom, ma);
2388 }
2389
2390 /* Create a piecewise multi-affine expression in the given space that maps each
2391  * input dimension to the corresponding output dimension.
2392  */
2393 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_identity(
2394         __isl_take isl_space *space)
2395 {
2396         return isl_pw_multi_aff_from_multi_aff(isl_multi_aff_identity(space));
2397 }
2398
2399 __isl_give isl_multi_aff *isl_multi_aff_add(__isl_take isl_multi_aff *maff1,
2400         __isl_take isl_multi_aff *maff2)
2401 {
2402         int i;
2403         isl_ctx *ctx;
2404
2405         maff1 = isl_multi_aff_cow(maff1);
2406         if (!maff1 || !maff2)
2407                 goto error;
2408
2409         ctx = isl_multi_aff_get_ctx(maff1);
2410         if (!isl_space_is_equal(maff1->space, maff2->space))
2411                 isl_die(ctx, isl_error_invalid,
2412                         "spaces don't match", goto error);
2413
2414         for (i = 0; i < maff1->n; ++i) {
2415                 maff1->p[i] = isl_aff_add(maff1->p[i],
2416                                             isl_aff_copy(maff2->p[i]));
2417                 if (!maff1->p[i])
2418                         goto error;
2419         }
2420
2421         isl_multi_aff_free(maff2);
2422         return maff1;
2423 error:
2424         isl_multi_aff_free(maff1);
2425         isl_multi_aff_free(maff2);
2426         return NULL;
2427 }
2428
2429 /* Given two multi-affine expressions A -> B and C -> D,
2430  * construct a multi-affine expression [A -> C] -> [B -> D].
2431  */
2432 __isl_give isl_multi_aff *isl_multi_aff_product(
2433         __isl_take isl_multi_aff *ma1, __isl_take isl_multi_aff *ma2)
2434 {
2435         int i;
2436         isl_aff *aff;
2437         isl_space *space;
2438         isl_multi_aff *res;
2439         int in1, in2, out1, out2;
2440
2441         in1 = isl_multi_aff_dim(ma1, isl_dim_in);
2442         in2 = isl_multi_aff_dim(ma2, isl_dim_in);
2443         out1 = isl_multi_aff_dim(ma1, isl_dim_out);
2444         out2 = isl_multi_aff_dim(ma2, isl_dim_out);
2445         space = isl_space_product(isl_multi_aff_get_space(ma1),
2446                                   isl_multi_aff_get_space(ma2));
2447         res = isl_multi_aff_alloc(isl_space_copy(space));
2448         space = isl_space_domain(space);
2449
2450         for (i = 0; i < out1; ++i) {
2451                 aff = isl_multi_aff_get_aff(ma1, i);
2452                 aff = isl_aff_insert_dims(aff, isl_dim_in, in1, in2);
2453                 aff = isl_aff_reset_domain_space(aff, isl_space_copy(space));
2454                 res = isl_multi_aff_set_aff(res, i, aff);
2455         }
2456
2457         for (i = 0; i < out2; ++i) {
2458                 aff = isl_multi_aff_get_aff(ma2, i);
2459                 aff = isl_aff_insert_dims(aff, isl_dim_in, 0, in1);
2460                 aff = isl_aff_reset_domain_space(aff, isl_space_copy(space));
2461                 res = isl_multi_aff_set_aff(res, out1 + i, aff);
2462         }
2463
2464         isl_space_free(space);
2465         isl_multi_aff_free(ma1);
2466         isl_multi_aff_free(ma2);
2467         return res;
2468 }
2469
2470 /* Exploit the equalities in "eq" to simplify the affine expressions.
2471  */
2472 static __isl_give isl_multi_aff *isl_multi_aff_substitute_equalities(
2473         __isl_take isl_multi_aff *maff, __isl_take isl_basic_set *eq)
2474 {
2475         int i;
2476
2477         maff = isl_multi_aff_cow(maff);
2478         if (!maff || !eq)
2479                 goto error;
2480
2481         for (i = 0; i < maff->n; ++i) {
2482                 maff->p[i] = isl_aff_substitute_equalities(maff->p[i],
2483                                                     isl_basic_set_copy(eq));
2484                 if (!maff->p[i])
2485                         goto error;
2486         }
2487
2488         isl_basic_set_free(eq);
2489         return maff;
2490 error:
2491         isl_basic_set_free(eq);
2492         isl_multi_aff_free(maff);
2493         return NULL;
2494 }
2495
2496 __isl_give isl_multi_aff *isl_multi_aff_scale(__isl_take isl_multi_aff *maff,
2497         isl_int f)
2498 {
2499         int i;
2500
2501         maff = isl_multi_aff_cow(maff);
2502         if (!maff)
2503                 return NULL;
2504
2505         for (i = 0; i < maff->n; ++i) {
2506                 maff->p[i] = isl_aff_scale(maff->p[i], f);
2507                 if (!maff->p[i])
2508                         return isl_multi_aff_free(maff);
2509         }
2510
2511         return maff;
2512 }
2513
2514 __isl_give isl_multi_aff *isl_multi_aff_add_on_domain(__isl_keep isl_set *dom,
2515         __isl_take isl_multi_aff *maff1, __isl_take isl_multi_aff *maff2)
2516 {
2517         maff1 = isl_multi_aff_add(maff1, maff2);
2518         maff1 = isl_multi_aff_gist(maff1, isl_set_copy(dom));
2519         return maff1;
2520 }
2521
2522 int isl_multi_aff_is_empty(__isl_keep isl_multi_aff *maff)
2523 {
2524         if (!maff)
2525                 return -1;
2526
2527         return 0;
2528 }
2529
2530 int isl_multi_aff_plain_is_equal(__isl_keep isl_multi_aff *maff1,
2531         __isl_keep isl_multi_aff *maff2)
2532 {
2533         int i;
2534         int equal;
2535
2536         if (!maff1 || !maff2)
2537                 return -1;
2538         if (maff1->n != maff2->n)
2539                 return 0;
2540         equal = isl_space_is_equal(maff1->space, maff2->space);
2541         if (equal < 0 || !equal)
2542                 return equal;
2543
2544         for (i = 0; i < maff1->n; ++i) {
2545                 equal = isl_aff_plain_is_equal(maff1->p[i], maff2->p[i]);
2546                 if (equal < 0 || !equal)
2547                         return equal;
2548         }
2549
2550         return 1;
2551 }
2552
2553 __isl_give isl_multi_aff *isl_multi_aff_set_dim_name(
2554         __isl_take isl_multi_aff *maff,
2555         enum isl_dim_type type, unsigned pos, const char *s)
2556 {
2557         int i;
2558
2559         maff = isl_multi_aff_cow(maff);
2560         if (!maff)
2561                 return NULL;
2562
2563         maff->space = isl_space_set_dim_name(maff->space, type, pos, s);
2564         if (!maff->space)
2565                 return isl_multi_aff_free(maff);
2566
2567         if (type == isl_dim_out)
2568                 return maff;
2569         for (i = 0; i < maff->n; ++i) {
2570                 maff->p[i] = isl_aff_set_dim_name(maff->p[i], type, pos, s);
2571                 if (!maff->p[i])
2572                         return isl_multi_aff_free(maff);
2573         }
2574
2575         return maff;
2576 }
2577
2578 __isl_give isl_multi_aff *isl_multi_aff_drop_dims(__isl_take isl_multi_aff *maff,
2579         enum isl_dim_type type, unsigned first, unsigned n)
2580 {
2581         int i;
2582
2583         maff = isl_multi_aff_cow(maff);
2584         if (!maff)
2585                 return NULL;
2586
2587         maff->space = isl_space_drop_dims(maff->space, type, first, n);
2588         if (!maff->space)
2589                 return isl_multi_aff_free(maff);
2590
2591         if (type == isl_dim_out) {
2592                 for (i = 0; i < n; ++i)
2593                         isl_aff_free(maff->p[first + i]);
2594                 for (i = first; i + n < maff->n; ++i)
2595                         maff->p[i] = maff->p[i + n];
2596                 maff->n -= n;
2597                 return maff;
2598         }
2599
2600         for (i = 0; i < maff->n; ++i) {
2601                 maff->p[i] = isl_aff_drop_dims(maff->p[i], type, first, n);
2602                 if (!maff->p[i])
2603                         return isl_multi_aff_free(maff);
2604         }
2605
2606         return maff;
2607 }
2608
2609 /* Return the set of domain elements where "ma1" is lexicographically
2610  * smaller than or equal to "ma2".
2611  */
2612 __isl_give isl_set *isl_multi_aff_lex_le_set(__isl_take isl_multi_aff *ma1,
2613         __isl_take isl_multi_aff *ma2)
2614 {
2615         return isl_multi_aff_lex_ge_set(ma2, ma1);
2616 }
2617
2618 /* Return the set of domain elements where "ma1" is lexicographically
2619  * greater than or equal to "ma2".
2620  */
2621 __isl_give isl_set *isl_multi_aff_lex_ge_set(__isl_take isl_multi_aff *ma1,
2622         __isl_take isl_multi_aff *ma2)
2623 {
2624         isl_space *space;
2625         isl_map *map1, *map2;
2626         isl_map *map, *ge;
2627
2628         map1 = isl_map_from_multi_aff(ma1);
2629         map2 = isl_map_from_multi_aff(ma2);
2630         map = isl_map_range_product(map1, map2);
2631         space = isl_space_range(isl_map_get_space(map));
2632         space = isl_space_domain(isl_space_unwrap(space));
2633         ge = isl_map_lex_ge(space);
2634         map = isl_map_intersect_range(map, isl_map_wrap(ge));
2635
2636         return isl_map_domain(map);
2637 }
2638
2639 #undef PW
2640 #define PW isl_pw_multi_aff
2641 #undef EL
2642 #define EL isl_multi_aff
2643 #undef EL_IS_ZERO
2644 #define EL_IS_ZERO is_empty
2645 #undef ZERO
2646 #define ZERO empty
2647 #undef IS_ZERO
2648 #define IS_ZERO is_empty
2649 #undef FIELD
2650 #define FIELD maff
2651 #undef DEFAULT_IS_ZERO
2652 #define DEFAULT_IS_ZERO 0
2653
2654 #define NO_NEG
2655 #define NO_EVAL
2656 #define NO_OPT
2657 #define NO_INVOLVES_DIMS
2658 #define NO_MOVE_DIMS
2659 #define NO_INSERT_DIMS
2660 #define NO_LIFT
2661 #define NO_MORPH
2662
2663 #include <isl_pw_templ.c>
2664
2665 #undef UNION
2666 #define UNION isl_union_pw_multi_aff
2667 #undef PART
2668 #define PART isl_pw_multi_aff
2669 #undef PARTS
2670 #define PARTS pw_multi_aff
2671 #define ALIGN_DOMAIN
2672
2673 #define NO_EVAL
2674
2675 #include <isl_union_templ.c>
2676
2677 /* Given a function "cmp" that returns the set of elements where
2678  * "ma1" is "better" than "ma2", return the intersection of this
2679  * set with "dom1" and "dom2".
2680  */
2681 static __isl_give isl_set *shared_and_better(__isl_keep isl_set *dom1,
2682         __isl_keep isl_set *dom2, __isl_keep isl_multi_aff *ma1,
2683         __isl_keep isl_multi_aff *ma2,
2684         __isl_give isl_set *(*cmp)(__isl_take isl_multi_aff *ma1,
2685                                     __isl_take isl_multi_aff *ma2))
2686 {
2687         isl_set *common;
2688         isl_set *better;
2689         int is_empty;
2690
2691         common = isl_set_intersect(isl_set_copy(dom1), isl_set_copy(dom2));
2692         is_empty = isl_set_plain_is_empty(common);
2693         if (is_empty >= 0 && is_empty)
2694                 return common;
2695         if (is_empty < 0)
2696                 return isl_set_free(common);
2697         better = cmp(isl_multi_aff_copy(ma1), isl_multi_aff_copy(ma2));
2698         better = isl_set_intersect(common, better);
2699
2700         return better;
2701 }
2702
2703 /* Given a function "cmp" that returns the set of elements where
2704  * "ma1" is "better" than "ma2", return a piecewise multi affine
2705  * expression defined on the union of the definition domains
2706  * of "pma1" and "pma2" that maps to the "best" of "pma1" and
2707  * "pma2" on each cell.  If only one of the two input functions
2708  * is defined on a given cell, then it is considered the best.
2709  */
2710 static __isl_give isl_pw_multi_aff *pw_multi_aff_union_opt(
2711         __isl_take isl_pw_multi_aff *pma1,
2712         __isl_take isl_pw_multi_aff *pma2,
2713         __isl_give isl_set *(*cmp)(__isl_take isl_multi_aff *ma1,
2714                                     __isl_take isl_multi_aff *ma2))
2715 {
2716         int i, j, n;
2717         isl_pw_multi_aff *res = NULL;
2718         isl_ctx *ctx;
2719         isl_set *set = NULL;
2720
2721         if (!pma1 || !pma2)
2722                 goto error;
2723
2724         ctx = isl_space_get_ctx(pma1->dim);
2725         if (!isl_space_is_equal(pma1->dim, pma2->dim))
2726                 isl_die(ctx, isl_error_invalid,
2727                         "arguments should live in the same space", goto error);
2728
2729         if (isl_pw_multi_aff_is_empty(pma1)) {
2730                 isl_pw_multi_aff_free(pma1);
2731                 return pma2;
2732         }
2733
2734         if (isl_pw_multi_aff_is_empty(pma2)) {
2735                 isl_pw_multi_aff_free(pma2);
2736                 return pma1;
2737         }
2738
2739         n = 2 * (pma1->n + 1) * (pma2->n + 1);
2740         res = isl_pw_multi_aff_alloc_size(isl_space_copy(pma1->dim), n);
2741
2742         for (i = 0; i < pma1->n; ++i) {
2743                 set = isl_set_copy(pma1->p[i].set);
2744                 for (j = 0; j < pma2->n; ++j) {
2745                         isl_set *better;
2746                         int is_empty;
2747
2748                         better = shared_and_better(pma2->p[j].set,
2749                                         pma1->p[i].set, pma2->p[j].maff,
2750                                         pma1->p[i].maff, cmp);
2751                         is_empty = isl_set_plain_is_empty(better);
2752                         if (is_empty < 0 || is_empty) {
2753                                 isl_set_free(better);
2754                                 if (is_empty < 0)
2755                                         goto error;
2756                                 continue;
2757                         }
2758                         set = isl_set_subtract(set, isl_set_copy(better));
2759
2760                         res = isl_pw_multi_aff_add_piece(res, better,
2761                                         isl_multi_aff_copy(pma2->p[j].maff));
2762                 }
2763                 res = isl_pw_multi_aff_add_piece(res, set,
2764                                         isl_multi_aff_copy(pma1->p[i].maff));
2765         }
2766
2767         for (j = 0; j < pma2->n; ++j) {
2768                 set = isl_set_copy(pma2->p[j].set);
2769                 for (i = 0; i < pma1->n; ++i)
2770                         set = isl_set_subtract(set,
2771                                         isl_set_copy(pma1->p[i].set));
2772                 res = isl_pw_multi_aff_add_piece(res, set,
2773                                         isl_multi_aff_copy(pma2->p[j].maff));
2774         }
2775
2776         isl_pw_multi_aff_free(pma1);
2777         isl_pw_multi_aff_free(pma2);
2778
2779         return res;
2780 error:
2781         isl_pw_multi_aff_free(pma1);
2782         isl_pw_multi_aff_free(pma2);
2783         isl_set_free(set);
2784         return isl_pw_multi_aff_free(res);
2785 }
2786
2787 static __isl_give isl_pw_multi_aff *pw_multi_aff_union_lexmax(
2788         __isl_take isl_pw_multi_aff *pma1,
2789         __isl_take isl_pw_multi_aff *pma2)
2790 {
2791         return pw_multi_aff_union_opt(pma1, pma2, &isl_multi_aff_lex_ge_set);
2792 }
2793
2794 /* Given two piecewise multi affine expressions, return a piecewise
2795  * multi-affine expression defined on the union of the definition domains
2796  * of the inputs that is equal to the lexicographic maximum of the two
2797  * inputs on each cell.  If only one of the two inputs is defined on
2798  * a given cell, then it is considered to be the maximum.
2799  */
2800 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_union_lexmax(
2801         __isl_take isl_pw_multi_aff *pma1,
2802         __isl_take isl_pw_multi_aff *pma2)
2803 {
2804         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
2805                                                     &pw_multi_aff_union_lexmax);
2806 }
2807
2808 static __isl_give isl_pw_multi_aff *pw_multi_aff_union_lexmin(
2809         __isl_take isl_pw_multi_aff *pma1,
2810         __isl_take isl_pw_multi_aff *pma2)
2811 {
2812         return pw_multi_aff_union_opt(pma1, pma2, &isl_multi_aff_lex_le_set);
2813 }
2814
2815 /* Given two piecewise multi affine expressions, return a piecewise
2816  * multi-affine expression defined on the union of the definition domains
2817  * of the inputs that is equal to the lexicographic minimum of the two
2818  * inputs on each cell.  If only one of the two inputs is defined on
2819  * a given cell, then it is considered to be the minimum.
2820  */
2821 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_union_lexmin(
2822         __isl_take isl_pw_multi_aff *pma1,
2823         __isl_take isl_pw_multi_aff *pma2)
2824 {
2825         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
2826                                                     &pw_multi_aff_union_lexmin);
2827 }
2828
2829 static __isl_give isl_pw_multi_aff *pw_multi_aff_add(
2830         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
2831 {
2832         return isl_pw_multi_aff_on_shared_domain(pma1, pma2,
2833                                                 &isl_multi_aff_add);
2834 }
2835
2836 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_add(
2837         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
2838 {
2839         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
2840                                                 &pw_multi_aff_add);
2841 }
2842
2843 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_union_add(
2844         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
2845 {
2846         return isl_pw_multi_aff_union_add_(pma1, pma2);
2847 }
2848
2849 /* Given two piecewise multi-affine expressions A -> B and C -> D,
2850  * construct a piecewise multi-affine expression [A -> C] -> [B -> D].
2851  */
2852 static __isl_give isl_pw_multi_aff *pw_multi_aff_product(
2853         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
2854 {
2855         int i, j, n;
2856         isl_space *space;
2857         isl_pw_multi_aff *res;
2858
2859         if (!pma1 || !pma2)
2860                 goto error;
2861
2862         n = pma1->n * pma2->n;
2863         space = isl_space_product(isl_space_copy(pma1->dim),
2864                                   isl_space_copy(pma2->dim));
2865         res = isl_pw_multi_aff_alloc_size(space, n);
2866
2867         for (i = 0; i < pma1->n; ++i) {
2868                 for (j = 0; j < pma2->n; ++j) {
2869                         isl_set *domain;
2870                         isl_multi_aff *ma;
2871
2872                         domain = isl_set_product(isl_set_copy(pma1->p[i].set),
2873                                                  isl_set_copy(pma2->p[j].set));
2874                         ma = isl_multi_aff_product(
2875                                         isl_multi_aff_copy(pma1->p[i].maff),
2876                                         isl_multi_aff_copy(pma2->p[i].maff));
2877                         res = isl_pw_multi_aff_add_piece(res, domain, ma);
2878                 }
2879         }
2880
2881         isl_pw_multi_aff_free(pma1);
2882         isl_pw_multi_aff_free(pma2);
2883         return res;
2884 error:
2885         isl_pw_multi_aff_free(pma1);
2886         isl_pw_multi_aff_free(pma2);
2887         return NULL;
2888 }
2889
2890 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_product(
2891         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
2892 {
2893         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
2894                                                 &pw_multi_aff_product);
2895 }
2896
2897 /* Construct a map mapping the domain of the piecewise multi-affine expression
2898  * to its range, with each dimension in the range equated to the
2899  * corresponding affine expression on its cell.
2900  */
2901 __isl_give isl_map *isl_map_from_pw_multi_aff(__isl_take isl_pw_multi_aff *pma)
2902 {
2903         int i;
2904         isl_map *map;
2905
2906         if (!pma)
2907                 return NULL;
2908
2909         map = isl_map_empty(isl_pw_multi_aff_get_space(pma));
2910
2911         for (i = 0; i < pma->n; ++i) {
2912                 isl_multi_aff *maff;
2913                 isl_basic_map *bmap;
2914                 isl_map *map_i;
2915
2916                 maff = isl_multi_aff_copy(pma->p[i].maff);
2917                 bmap = isl_basic_map_from_multi_aff(maff);
2918                 map_i = isl_map_from_basic_map(bmap);
2919                 map_i = isl_map_intersect_domain(map_i,
2920                                                 isl_set_copy(pma->p[i].set));
2921                 map = isl_map_union_disjoint(map, map_i);
2922         }
2923
2924         isl_pw_multi_aff_free(pma);
2925         return map;
2926 }
2927
2928 __isl_give isl_set *isl_set_from_pw_multi_aff(__isl_take isl_pw_multi_aff *pma)
2929 {
2930         if (!isl_space_is_set(pma->dim))
2931                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
2932                         "isl_pw_multi_aff cannot be converted into an isl_set",
2933                         return isl_pw_multi_aff_free(pma));
2934
2935         return isl_map_from_pw_multi_aff(pma);
2936 }
2937
2938 /* Given a basic map with a single output dimension that is defined
2939  * in terms of the parameters and input dimensions using an equality,
2940  * extract an isl_aff that expresses the output dimension in terms
2941  * of the parameters and input dimensions.
2942  *
2943  * Since some applications expect the result of isl_pw_multi_aff_from_map
2944  * to only contain integer affine expressions, we compute the floor
2945  * of the expression before returning.
2946  *
2947  * This function shares some similarities with
2948  * isl_basic_map_has_defining_equality and isl_constraint_get_bound.
2949  */
2950 static __isl_give isl_aff *extract_isl_aff_from_basic_map(
2951         __isl_take isl_basic_map *bmap)
2952 {
2953         int i;
2954         unsigned offset;
2955         unsigned total;
2956         isl_local_space *ls;
2957         isl_aff *aff;
2958
2959         if (!bmap)
2960                 return NULL;
2961         if (isl_basic_map_dim(bmap, isl_dim_out) != 1)
2962                 isl_die(isl_basic_map_get_ctx(bmap), isl_error_invalid,
2963                         "basic map should have a single output dimension",
2964                         goto error);
2965         offset = isl_basic_map_offset(bmap, isl_dim_out);
2966         total = isl_basic_map_total_dim(bmap);
2967         for (i = 0; i < bmap->n_eq; ++i) {
2968                 if (isl_int_is_zero(bmap->eq[i][offset]))
2969                         continue;
2970                 if (isl_seq_first_non_zero(bmap->eq[i] + offset + 1,
2971                                            1 + total - (offset + 1)) != -1)
2972                         continue;
2973                 break;
2974         }
2975         if (i >= bmap->n_eq)
2976                 isl_die(isl_basic_map_get_ctx(bmap), isl_error_invalid,
2977                         "unable to find suitable equality", goto error);
2978         ls = isl_basic_map_get_local_space(bmap);
2979         aff = isl_aff_alloc(isl_local_space_domain(ls));
2980         if (!aff)
2981                 goto error;
2982         if (isl_int_is_neg(bmap->eq[i][offset]))
2983                 isl_seq_cpy(aff->v->el + 1, bmap->eq[i], offset);
2984         else
2985                 isl_seq_neg(aff->v->el + 1, bmap->eq[i], offset);
2986         isl_seq_clr(aff->v->el + 1 + offset, aff->v->size - (1 + offset));
2987         isl_int_abs(aff->v->el[0], bmap->eq[i][offset]);
2988         isl_basic_map_free(bmap);
2989
2990         aff = isl_aff_remove_unused_divs(aff);
2991         aff = isl_aff_floor(aff);
2992         return aff;
2993 error:
2994         isl_basic_map_free(bmap);
2995         return NULL;
2996 }
2997
2998 /* Given a basic map where each output dimension is defined
2999  * in terms of the parameters and input dimensions using an equality,
3000  * extract an isl_multi_aff that expresses the output dimensions in terms
3001  * of the parameters and input dimensions.
3002  */
3003 static __isl_give isl_multi_aff *extract_isl_multi_aff_from_basic_map(
3004         __isl_take isl_basic_map *bmap)
3005 {
3006         int i;
3007         unsigned n_out;
3008         isl_multi_aff *ma;
3009
3010         if (!bmap)
3011                 return NULL;
3012
3013         ma = isl_multi_aff_alloc(isl_basic_map_get_space(bmap));
3014         n_out = isl_basic_map_dim(bmap, isl_dim_out);
3015
3016         for (i = 0; i < n_out; ++i) {
3017                 isl_basic_map *bmap_i;
3018                 isl_aff *aff;
3019
3020                 bmap_i = isl_basic_map_copy(bmap);
3021                 bmap_i = isl_basic_map_project_out(bmap_i, isl_dim_out,
3022                                                         i + 1, n_out - (1 + i));
3023                 bmap_i = isl_basic_map_project_out(bmap_i, isl_dim_out, 0, i);
3024                 aff = extract_isl_aff_from_basic_map(bmap_i);
3025                 ma = isl_multi_aff_set_aff(ma, i, aff);
3026         }
3027
3028         isl_basic_map_free(bmap);
3029
3030         return ma;
3031 }
3032
3033 /* Create an isl_pw_multi_aff that is equivalent to
3034  * isl_map_intersect_domain(isl_map_from_basic_map(bmap), domain).
3035  * The given basic map is such that each output dimension is defined
3036  * in terms of the parameters and input dimensions using an equality.
3037  */
3038 static __isl_give isl_pw_multi_aff *plain_pw_multi_aff_from_map(
3039         __isl_take isl_set *domain, __isl_take isl_basic_map *bmap)
3040 {
3041         isl_multi_aff *ma;
3042
3043         ma = extract_isl_multi_aff_from_basic_map(bmap);
3044         return isl_pw_multi_aff_alloc(domain, ma);
3045 }
3046
3047 /* Try and create an isl_pw_multi_aff that is equivalent to the given isl_map.
3048  * This obivously only works if the input "map" is single-valued.
3049  * If so, we compute the lexicographic minimum of the image in the form
3050  * of an isl_pw_multi_aff.  Since the image is unique, it is equal
3051  * to its lexicographic minimum.
3052  * If the input is not single-valued, we produce an error.
3053  *
3054  * As a special case, we first check if all output dimensions are uniquely
3055  * defined in terms of the parameters and input dimensions over the entire
3056  * domain.  If so, we extract the desired isl_pw_multi_aff directly
3057  * from the affine hull of "map" and its domain.
3058  */
3059 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_from_map(__isl_take isl_map *map)
3060 {
3061         int i;
3062         int sv;
3063         isl_pw_multi_aff *pma;
3064         isl_basic_map *hull;
3065
3066         if (!map)
3067                 return NULL;
3068
3069         hull = isl_map_affine_hull(isl_map_copy(map));
3070         sv = isl_basic_map_plain_is_single_valued(hull);
3071         if (sv >= 0 && sv)
3072                 return plain_pw_multi_aff_from_map(isl_map_domain(map), hull);
3073         isl_basic_map_free(hull);
3074         if (sv < 0)
3075                 goto error;
3076
3077         sv = isl_map_is_single_valued(map);
3078         if (sv < 0)
3079                 goto error;
3080         if (!sv)
3081                 isl_die(isl_map_get_ctx(map), isl_error_invalid,
3082                         "map is not single-valued", goto error);
3083         map = isl_map_make_disjoint(map);
3084         if (!map)
3085                 return NULL;
3086
3087         pma = isl_pw_multi_aff_empty(isl_map_get_space(map));
3088
3089         for (i = 0; i < map->n; ++i) {
3090                 isl_pw_multi_aff *pma_i;
3091                 isl_basic_map *bmap;
3092                 bmap = isl_basic_map_copy(map->p[i]);
3093                 pma_i = isl_basic_map_lexmin_pw_multi_aff(bmap);
3094                 pma = isl_pw_multi_aff_add_disjoint(pma, pma_i);
3095         }
3096
3097         isl_map_free(map);
3098         return pma;
3099 error:
3100         isl_map_free(map);
3101         return NULL;
3102 }
3103
3104 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_from_set(__isl_take isl_set *set)
3105 {
3106         return isl_pw_multi_aff_from_map(set);
3107 }
3108
3109 /* Return the piecewise affine expression "set ? 1 : 0".
3110  */
3111 __isl_give isl_pw_aff *isl_set_indicator_function(__isl_take isl_set *set)
3112 {
3113         isl_pw_aff *pa;
3114         isl_space *space = isl_set_get_space(set);
3115         isl_local_space *ls = isl_local_space_from_space(space);
3116         isl_aff *zero = isl_aff_zero_on_domain(isl_local_space_copy(ls));
3117         isl_aff *one = isl_aff_zero_on_domain(ls);
3118
3119         one = isl_aff_add_constant_si(one, 1);
3120         pa = isl_pw_aff_alloc(isl_set_copy(set), one);
3121         set = isl_set_complement(set);
3122         pa = isl_pw_aff_add_disjoint(pa, isl_pw_aff_alloc(set, zero));
3123
3124         return pa;
3125 }
3126
3127 /* Plug in "subs" for dimension "type", "pos" of "aff".
3128  *
3129  * Let i be the dimension to replace and let "subs" be of the form
3130  *
3131  *      f/d
3132  *
3133  * and "aff" of the form
3134  *
3135  *      (a i + g)/m
3136  *
3137  * The result is
3138  *
3139  *      (a f + d g')/(m d)
3140  *
3141  * where g' is the result of plugging in "subs" in each of the integer
3142  * divisions in g.
3143  */
3144 __isl_give isl_aff *isl_aff_substitute(__isl_take isl_aff *aff,
3145         enum isl_dim_type type, unsigned pos, __isl_keep isl_aff *subs)
3146 {
3147         isl_ctx *ctx;
3148         isl_int v;
3149
3150         aff = isl_aff_cow(aff);
3151         if (!aff || !subs)
3152                 return isl_aff_free(aff);
3153
3154         ctx = isl_aff_get_ctx(aff);
3155         if (!isl_space_is_equal(aff->ls->dim, subs->ls->dim))
3156                 isl_die(ctx, isl_error_invalid,
3157                         "spaces don't match", return isl_aff_free(aff));
3158         if (isl_local_space_dim(subs->ls, isl_dim_div) != 0)
3159                 isl_die(ctx, isl_error_unsupported,
3160                         "cannot handle divs yet", return isl_aff_free(aff));
3161
3162         aff->ls = isl_local_space_substitute(aff->ls, type, pos, subs);
3163         if (!aff->ls)
3164                 return isl_aff_free(aff);
3165
3166         aff->v = isl_vec_cow(aff->v);
3167         if (!aff->v)
3168                 return isl_aff_free(aff);
3169
3170         pos += isl_local_space_offset(aff->ls, type);
3171
3172         isl_int_init(v);
3173         isl_seq_substitute(aff->v->el, pos, subs->v->el,
3174                             aff->v->size, subs->v->size, v);
3175         isl_int_clear(v);
3176
3177         return aff;
3178 }
3179
3180 /* Plug in "subs" for dimension "type", "pos" in each of the affine
3181  * expressions in "maff".
3182  */
3183 __isl_give isl_multi_aff *isl_multi_aff_substitute(
3184         __isl_take isl_multi_aff *maff, enum isl_dim_type type, unsigned pos,
3185         __isl_keep isl_aff *subs)
3186 {
3187         int i;
3188
3189         maff = isl_multi_aff_cow(maff);
3190         if (!maff || !subs)
3191                 return isl_multi_aff_free(maff);
3192
3193         if (type == isl_dim_in)
3194                 type = isl_dim_set;
3195
3196         for (i = 0; i < maff->n; ++i) {
3197                 maff->p[i] = isl_aff_substitute(maff->p[i], type, pos, subs);
3198                 if (!maff->p[i])
3199                         return isl_multi_aff_free(maff);
3200         }
3201
3202         return maff;
3203 }
3204
3205 /* Plug in "subs" for dimension "type", "pos" of "pma".
3206  *
3207  * pma is of the form
3208  *
3209  *      A_i(v) -> M_i(v)
3210  *
3211  * while subs is of the form
3212  *
3213  *      v' = B_j(v) -> S_j
3214  *
3215  * Each pair i,j such that C_ij = A_i \cap B_i is non-empty
3216  * has a contribution in the result, in particular
3217  *
3218  *      C_ij(S_j) -> M_i(S_j)
3219  *
3220  * Note that plugging in S_j in C_ij may also result in an empty set
3221  * and this contribution should simply be discarded.
3222  */
3223 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_substitute(
3224         __isl_take isl_pw_multi_aff *pma, enum isl_dim_type type, unsigned pos,
3225         __isl_keep isl_pw_aff *subs)
3226 {
3227         int i, j, n;
3228         isl_pw_multi_aff *res;
3229
3230         if (!pma || !subs)
3231                 return isl_pw_multi_aff_free(pma);
3232
3233         n = pma->n * subs->n;
3234         res = isl_pw_multi_aff_alloc_size(isl_space_copy(pma->dim), n);
3235
3236         for (i = 0; i < pma->n; ++i) {
3237                 for (j = 0; j < subs->n; ++j) {
3238                         isl_set *common;
3239                         isl_multi_aff *res_ij;
3240                         common = isl_set_intersect(
3241                                         isl_set_copy(pma->p[i].set),
3242                                         isl_set_copy(subs->p[j].set));
3243                         common = isl_set_substitute(common,
3244                                         type, pos, subs->p[j].aff);
3245                         if (isl_set_plain_is_empty(common)) {
3246                                 isl_set_free(common);
3247                                 continue;
3248                         }
3249
3250                         res_ij = isl_multi_aff_substitute(
3251                                         isl_multi_aff_copy(pma->p[i].maff),
3252                                         type, pos, subs->p[j].aff);
3253
3254                         res = isl_pw_multi_aff_add_piece(res, common, res_ij);
3255                 }
3256         }
3257
3258         isl_pw_multi_aff_free(pma);
3259         return res;
3260 }
3261
3262 /* Extend the local space of "dst" to include the divs
3263  * in the local space of "src".
3264  */
3265 __isl_give isl_aff *isl_aff_align_divs(__isl_take isl_aff *dst,
3266         __isl_keep isl_aff *src)
3267 {
3268         isl_ctx *ctx;
3269         int *exp1 = NULL;
3270         int *exp2 = NULL;
3271         isl_mat *div;
3272
3273         if (!src || !dst)
3274                 return isl_aff_free(dst);
3275
3276         ctx = isl_aff_get_ctx(src);
3277         if (!isl_space_is_equal(src->ls->dim, dst->ls->dim))
3278                 isl_die(ctx, isl_error_invalid,
3279                         "spaces don't match", goto error);
3280
3281         if (src->ls->div->n_row == 0)
3282                 return dst;
3283
3284         exp1 = isl_alloc_array(ctx, int, src->ls->div->n_row);
3285         exp2 = isl_alloc_array(ctx, int, dst->ls->div->n_row);
3286         if (!exp1 || !exp2)
3287                 goto error;
3288
3289         div = isl_merge_divs(src->ls->div, dst->ls->div, exp1, exp2);
3290         dst = isl_aff_expand_divs(dst, div, exp2);
3291         free(exp1);
3292         free(exp2);
3293
3294         return dst;
3295 error:
3296         free(exp1);
3297         free(exp2);
3298         return isl_aff_free(dst);
3299 }
3300
3301 /* Adjust the local spaces of the affine expressions in "maff"
3302  * such that they all have the save divs.
3303  */
3304 __isl_give isl_multi_aff *isl_multi_aff_align_divs(
3305         __isl_take isl_multi_aff *maff)
3306 {
3307         int i;
3308
3309         if (!maff)
3310                 return NULL;
3311         if (maff->n == 0)
3312                 return maff;
3313         maff = isl_multi_aff_cow(maff);
3314         if (!maff)
3315                 return NULL;
3316
3317         for (i = 1; i < maff->n; ++i)
3318                 maff->p[0] = isl_aff_align_divs(maff->p[0], maff->p[i]);
3319         for (i = 1; i < maff->n; ++i) {
3320                 maff->p[i] = isl_aff_align_divs(maff->p[i], maff->p[0]);
3321                 if (!maff->p[i])
3322                         return isl_multi_aff_free(maff);
3323         }
3324
3325         return maff;
3326 }
3327
3328 __isl_give isl_aff *isl_aff_lift(__isl_take isl_aff *aff)
3329 {
3330         aff = isl_aff_cow(aff);
3331         if (!aff)
3332                 return NULL;
3333
3334         aff->ls = isl_local_space_lift(aff->ls);
3335         if (!aff->ls)
3336                 return isl_aff_free(aff);
3337
3338         return aff;
3339 }
3340
3341 /* Lift "maff" to a space with extra dimensions such that the result
3342  * has no more existentially quantified variables.
3343  * If "ls" is not NULL, then *ls is assigned the local space that lies
3344  * at the basis of the lifting applied to "maff".
3345  */
3346 __isl_give isl_multi_aff *isl_multi_aff_lift(__isl_take isl_multi_aff *maff,
3347         __isl_give isl_local_space **ls)
3348 {
3349         int i;
3350         isl_space *space;
3351         unsigned n_div;
3352
3353         if (ls)
3354                 *ls = NULL;
3355
3356         if (!maff)
3357                 return NULL;
3358
3359         if (maff->n == 0) {
3360                 if (ls) {
3361                         isl_space *space = isl_multi_aff_get_domain_space(maff);
3362                         *ls = isl_local_space_from_space(space);
3363                         if (!*ls)
3364                                 return isl_multi_aff_free(maff);
3365                 }
3366                 return maff;
3367         }
3368
3369         maff = isl_multi_aff_cow(maff);
3370         maff = isl_multi_aff_align_divs(maff);
3371         if (!maff)
3372                 return NULL;
3373
3374         n_div = isl_aff_dim(maff->p[0], isl_dim_div);
3375         space = isl_multi_aff_get_space(maff);
3376         space = isl_space_lift(isl_space_domain(space), n_div);
3377         space = isl_space_extend_domain_with_range(space,
3378                                                 isl_multi_aff_get_space(maff));
3379         if (!space)
3380                 return isl_multi_aff_free(maff);
3381         isl_space_free(maff->space);
3382         maff->space = space;
3383
3384         if (ls) {
3385                 *ls = isl_aff_get_domain_local_space(maff->p[0]);
3386                 if (!*ls)
3387                         return isl_multi_aff_free(maff);
3388         }
3389
3390         for (i = 0; i < maff->n; ++i) {
3391                 maff->p[i] = isl_aff_lift(maff->p[i]);
3392                 if (!maff->p[i])
3393                         goto error;
3394         }
3395
3396         return maff;
3397 error:
3398         if (ls)
3399                 isl_local_space_free(*ls);
3400         return isl_multi_aff_free(maff);
3401 }
3402
3403
3404 /* Extract an isl_pw_aff corresponding to output dimension "pos" of "pma".
3405  */
3406 __isl_give isl_pw_aff *isl_pw_multi_aff_get_pw_aff(
3407         __isl_keep isl_pw_multi_aff *pma, int pos)
3408 {
3409         int i;
3410         int n_out;
3411         isl_space *space;
3412         isl_pw_aff *pa;
3413
3414         if (!pma)
3415                 return NULL;
3416
3417         n_out = isl_pw_multi_aff_dim(pma, isl_dim_out);
3418         if (pos < 0 || pos >= n_out)
3419                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
3420                         "index out of bounds", return NULL);
3421
3422         space = isl_pw_multi_aff_get_space(pma);
3423         space = isl_space_drop_dims(space, isl_dim_out,
3424                                     pos + 1, n_out - pos - 1);
3425         space = isl_space_drop_dims(space, isl_dim_out, 0, pos);
3426
3427         pa = isl_pw_aff_alloc_size(space, pma->n);
3428         for (i = 0; i < pma->n; ++i) {
3429                 isl_aff *aff;
3430                 aff = isl_multi_aff_get_aff(pma->p[i].maff, pos);
3431                 pa = isl_pw_aff_add_piece(pa, isl_set_copy(pma->p[i].set), aff);
3432         }
3433
3434         return pa;
3435 }
3436
3437 /* Return an isl_pw_multi_aff with the given "set" as domain and
3438  * an unnamed zero-dimensional range.
3439  */
3440 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_from_domain(
3441         __isl_take isl_set *set)
3442 {
3443         isl_multi_aff *ma;
3444         isl_space *space;
3445
3446         space = isl_set_get_space(set);
3447         space = isl_space_from_domain(space);
3448         ma = isl_multi_aff_zero(space);
3449         return isl_pw_multi_aff_alloc(set, ma);
3450 }
3451
3452 /* Add an isl_pw_multi_aff with the given "set" as domain and
3453  * an unnamed zero-dimensional range to *user.
3454  */
3455 static int add_pw_multi_aff_from_domain(__isl_take isl_set *set, void *user)
3456 {
3457         isl_union_pw_multi_aff **upma = user;
3458         isl_pw_multi_aff *pma;
3459
3460         pma = isl_pw_multi_aff_from_domain(set);
3461         *upma = isl_union_pw_multi_aff_add_pw_multi_aff(*upma, pma);
3462
3463         return 0;
3464 }
3465
3466 /* Return an isl_union_pw_multi_aff with the given "uset" as domain and
3467  * an unnamed zero-dimensional range.
3468  */
3469 __isl_give isl_union_pw_multi_aff *isl_union_pw_multi_aff_from_domain(
3470         __isl_take isl_union_set *uset)
3471 {
3472         isl_space *space;
3473         isl_union_pw_multi_aff *upma;
3474
3475         if (!uset)
3476                 return NULL;
3477
3478         space = isl_union_set_get_space(uset);
3479         upma = isl_union_pw_multi_aff_empty(space);
3480
3481         if (isl_union_set_foreach_set(uset,
3482                                     &add_pw_multi_aff_from_domain, &upma) < 0)
3483                 goto error;
3484
3485         isl_union_set_free(uset);
3486         return upma;
3487 error:
3488         isl_union_set_free(uset);
3489         isl_union_pw_multi_aff_free(upma);
3490         return NULL;
3491 }
3492
3493 /* Convert "pma" to an isl_map and add it to *umap.
3494  */
3495 static int map_from_pw_multi_aff(__isl_take isl_pw_multi_aff *pma, void *user)
3496 {
3497         isl_union_map **umap = user;
3498         isl_map *map;
3499
3500         map = isl_map_from_pw_multi_aff(pma);
3501         *umap = isl_union_map_add_map(*umap, map);
3502
3503         return 0;
3504 }
3505
3506 /* Construct a union map mapping the domain of the union
3507  * piecewise multi-affine expression to its range, with each dimension
3508  * in the range equated to the corresponding affine expression on its cell.
3509  */
3510 __isl_give isl_union_map *isl_union_map_from_union_pw_multi_aff(
3511         __isl_take isl_union_pw_multi_aff *upma)
3512 {
3513         isl_space *space;
3514         isl_union_map *umap;
3515
3516         if (!upma)
3517                 return NULL;
3518
3519         space = isl_union_pw_multi_aff_get_space(upma);
3520         umap = isl_union_map_empty(space);
3521
3522         if (isl_union_pw_multi_aff_foreach_pw_multi_aff(upma,
3523                                         &map_from_pw_multi_aff, &umap) < 0)
3524                 goto error;
3525
3526         isl_union_pw_multi_aff_free(upma);
3527         return umap;
3528 error:
3529         isl_union_pw_multi_aff_free(upma);
3530         isl_union_map_free(umap);
3531         return NULL;
3532 }
3533
3534 /* Local data for bin_entry and the callback "fn".
3535  */
3536 struct isl_union_pw_multi_aff_bin_data {
3537         isl_union_pw_multi_aff *upma2;
3538         isl_union_pw_multi_aff *res;
3539         isl_pw_multi_aff *pma;
3540         int (*fn)(void **entry, void *user);
3541 };
3542
3543 /* Given an isl_pw_multi_aff from upma1, store it in data->pma
3544  * and call data->fn for each isl_pw_multi_aff in data->upma2.
3545  */
3546 static int bin_entry(void **entry, void *user)
3547 {
3548         struct isl_union_pw_multi_aff_bin_data *data = user;
3549         isl_pw_multi_aff *pma = *entry;
3550
3551         data->pma = pma;
3552         if (isl_hash_table_foreach(data->upma2->dim->ctx, &data->upma2->table,
3553                                    data->fn, data) < 0)
3554                 return -1;
3555
3556         return 0;
3557 }
3558
3559 /* Call "fn" on each pair of isl_pw_multi_affs in "upma1" and "upma2".
3560  * The isl_pw_multi_aff from upma1 is stored in data->pma (where data is
3561  * passed as user field) and the isl_pw_multi_aff from upma2 is available
3562  * as *entry.  The callback should adjust data->res if desired.
3563  */
3564 static __isl_give isl_union_pw_multi_aff *bin_op(
3565         __isl_take isl_union_pw_multi_aff *upma1,
3566         __isl_take isl_union_pw_multi_aff *upma2,
3567         int (*fn)(void **entry, void *user))
3568 {
3569         isl_space *space;
3570         struct isl_union_pw_multi_aff_bin_data data = { NULL, NULL, NULL, fn };
3571
3572         space = isl_union_pw_multi_aff_get_space(upma2);
3573         upma1 = isl_union_pw_multi_aff_align_params(upma1, space);
3574         space = isl_union_pw_multi_aff_get_space(upma1);
3575         upma2 = isl_union_pw_multi_aff_align_params(upma2, space);
3576
3577         if (!upma1 || !upma2)
3578                 goto error;
3579
3580         data.upma2 = upma2;
3581         data.res = isl_union_pw_multi_aff_alloc(isl_space_copy(upma1->dim),
3582                                        upma1->table.n);
3583         if (isl_hash_table_foreach(upma1->dim->ctx, &upma1->table,
3584                                    &bin_entry, &data) < 0)
3585                 goto error;
3586
3587         isl_union_pw_multi_aff_free(upma1);
3588         isl_union_pw_multi_aff_free(upma2);
3589         return data.res;
3590 error:
3591         isl_union_pw_multi_aff_free(upma1);
3592         isl_union_pw_multi_aff_free(upma2);
3593         isl_union_pw_multi_aff_free(data.res);
3594         return NULL;
3595 }
3596
3597 /* Given two isl_multi_affs A -> B and C -> D,
3598  * construct an isl_multi_aff (A * C) -> (B, D).
3599  */
3600 __isl_give isl_multi_aff *isl_multi_aff_flat_range_product(
3601         __isl_take isl_multi_aff *ma1, __isl_take isl_multi_aff *ma2)
3602 {
3603         int i, n1, n2;
3604         isl_aff *aff;
3605         isl_space *space;
3606         isl_multi_aff *res;
3607
3608         if (!ma1 || !ma2)
3609                 goto error;
3610
3611         space = isl_space_range_product(isl_multi_aff_get_space(ma1),
3612                                         isl_multi_aff_get_space(ma2));
3613         space = isl_space_flatten_range(space);
3614         res = isl_multi_aff_alloc(space);
3615
3616         n1 = isl_multi_aff_dim(ma1, isl_dim_out);
3617         n2 = isl_multi_aff_dim(ma2, isl_dim_out);
3618
3619         for (i = 0; i < n1; ++i) {
3620                 aff = isl_multi_aff_get_aff(ma1, i);
3621                 res = isl_multi_aff_set_aff(res, i, aff);
3622         }
3623
3624         for (i = 0; i < n2; ++i) {
3625                 aff = isl_multi_aff_get_aff(ma2, i);
3626                 res = isl_multi_aff_set_aff(res, n1 + i, aff);
3627         }
3628
3629         isl_multi_aff_free(ma1);
3630         isl_multi_aff_free(ma2);
3631         return res;
3632 error:
3633         isl_multi_aff_free(ma1);
3634         isl_multi_aff_free(ma2);
3635         return NULL;
3636 }
3637
3638 /* Given two aligned isl_pw_multi_affs A -> B and C -> D,
3639  * construct an isl_pw_multi_aff (A * C) -> (B, D).
3640  */
3641 static __isl_give isl_pw_multi_aff *pw_multi_aff_flat_range_product(
3642         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
3643 {
3644         isl_space *space;
3645
3646         space = isl_space_range_product(isl_pw_multi_aff_get_space(pma1),
3647                                         isl_pw_multi_aff_get_space(pma2));
3648         space = isl_space_flatten_range(space);
3649         return isl_pw_multi_aff_on_shared_domain_in(pma1, pma2, space,
3650                                             &isl_multi_aff_flat_range_product);
3651 }
3652
3653 /* Given two isl_pw_multi_affs A -> B and C -> D,
3654  * construct an isl_pw_multi_aff (A * C) -> (B, D).
3655  */
3656 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_flat_range_product(
3657         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
3658 {
3659         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
3660                                             &pw_multi_aff_flat_range_product);
3661 }
3662
3663 /* If data->pma and *entry have the same domain space, then compute
3664  * their flat range product and the result to data->res.
3665  */
3666 static int flat_range_product_entry(void **entry, void *user)
3667 {
3668         struct isl_union_pw_multi_aff_bin_data *data = user;
3669         isl_pw_multi_aff *pma2 = *entry;
3670
3671         if (!isl_space_tuple_match(data->pma->dim, isl_dim_in,
3672                                  pma2->dim, isl_dim_in))
3673                 return 0;
3674
3675         pma2 = isl_pw_multi_aff_flat_range_product(
3676                                         isl_pw_multi_aff_copy(data->pma),
3677                                         isl_pw_multi_aff_copy(pma2));
3678
3679         data->res = isl_union_pw_multi_aff_add_pw_multi_aff(data->res, pma2);
3680
3681         return 0;
3682 }
3683
3684 /* Given two isl_union_pw_multi_affs A -> B and C -> D,
3685  * construct an isl_union_pw_multi_aff (A * C) -> (B, D).
3686  */
3687 __isl_give isl_union_pw_multi_aff *isl_union_pw_multi_aff_flat_range_product(
3688         __isl_take isl_union_pw_multi_aff *upma1,
3689         __isl_take isl_union_pw_multi_aff *upma2)
3690 {
3691         return bin_op(upma1, upma2, &flat_range_product_entry);
3692 }
3693
3694 /* Replace the affine expressions at position "pos" in "pma" by "pa".
3695  * The parameters are assumed to have been aligned.
3696  *
3697  * The implementation essentially performs an isl_pw_*_on_shared_domain,
3698  * except that it works on two different isl_pw_* types.
3699  */
3700 static __isl_give isl_pw_multi_aff *pw_multi_aff_set_pw_aff(
3701         __isl_take isl_pw_multi_aff *pma, unsigned pos,
3702         __isl_take isl_pw_aff *pa)
3703 {
3704         int i, j, n;
3705         isl_pw_multi_aff *res = NULL;
3706
3707         if (!pma || !pa)
3708                 goto error;
3709
3710         if (!isl_space_tuple_match(pma->dim, isl_dim_in, pa->dim, isl_dim_in))
3711                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
3712                         "domains don't match", goto error);
3713         if (pos >= isl_pw_multi_aff_dim(pma, isl_dim_out))
3714                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
3715                         "index out of bounds", goto error);
3716
3717         n = pma->n * pa->n;
3718         res = isl_pw_multi_aff_alloc_size(isl_pw_multi_aff_get_space(pma), n);
3719
3720         for (i = 0; i < pma->n; ++i) {
3721                 for (j = 0; j < pa->n; ++j) {
3722                         isl_set *common;
3723                         isl_multi_aff *res_ij;
3724                         int empty;
3725
3726                         common = isl_set_intersect(isl_set_copy(pma->p[i].set),
3727                                                    isl_set_copy(pa->p[j].set));
3728                         empty = isl_set_plain_is_empty(common);
3729                         if (empty < 0 || empty) {
3730                                 isl_set_free(common);
3731                                 if (empty < 0)
3732                                         goto error;
3733                                 continue;
3734                         }
3735
3736                         res_ij = isl_multi_aff_set_aff(
3737                                         isl_multi_aff_copy(pma->p[i].maff), pos,
3738                                         isl_aff_copy(pa->p[j].aff));
3739                         res_ij = isl_multi_aff_gist(res_ij,
3740                                         isl_set_copy(common));
3741
3742                         res = isl_pw_multi_aff_add_piece(res, common, res_ij);
3743                 }
3744         }
3745
3746         isl_pw_multi_aff_free(pma);
3747         isl_pw_aff_free(pa);
3748         return res;
3749 error:
3750         isl_pw_multi_aff_free(pma);
3751         isl_pw_aff_free(pa);
3752         return isl_pw_multi_aff_free(res);
3753 }
3754
3755 /* Replace the affine expressions at position "pos" in "pma" by "pa".
3756  */
3757 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_set_pw_aff(
3758         __isl_take isl_pw_multi_aff *pma, unsigned pos,
3759         __isl_take isl_pw_aff *pa)
3760 {
3761         if (!pma || !pa)
3762                 goto error;
3763         if (isl_space_match(pma->dim, isl_dim_param, pa->dim, isl_dim_param))
3764                 return pw_multi_aff_set_pw_aff(pma, pos, pa);
3765         if (!isl_space_has_named_params(pma->dim) ||
3766             !isl_space_has_named_params(pa->dim))
3767                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
3768                         "unaligned unnamed parameters", goto error);
3769         pma = isl_pw_multi_aff_align_params(pma, isl_pw_aff_get_space(pa));
3770         pa = isl_pw_aff_align_params(pa, isl_pw_multi_aff_get_space(pma));
3771         return pw_multi_aff_set_pw_aff(pma, pos, pa);
3772 error:
3773         isl_pw_multi_aff_free(pma);
3774         isl_pw_aff_free(pa);
3775         return NULL;
3776 }