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 /* Merge divs "a" and "b" in "aff", which is assumed to be non-NULL.
 763  *
 764  * We currently do not actually remove div "b", but simply add its
 765  * coefficient to that of "a" and then zero it out.
 766  */
 767 static __isl_give isl_aff *merge_divs(__isl_take isl_aff *aff, int a, int b)
 768 {
 769         unsigned off = isl_local_space_offset(aff->ls, isl_dim_div);
 770
 771         if (isl_int_is_zero(aff->v->el[1 + off + b]))
 772                 return aff;
 773
 774         aff->v = isl_vec_cow(aff->v);
 775         if (!aff->v)
 776                 return isl_aff_free(aff);
 777
 778         isl_int_add(aff->v->el[1 + off + a],
 779                     aff->v->el[1 + off + a], aff->v->el[1 + off + b]);
 780         isl_int_set_si(aff->v->el[1 + off + b], 0);
 781
 782         return aff;
 783 }
 784
 785 /* Sort the divs in the local space of "aff" according to
 786  * the comparison function "cmp_row" in isl_local_space.c,
 787  * combining the coefficients of identical divs.
 788  *
 789  * Reordering divs does not change the semantics of "aff",
 790  * so there is no need to call isl_aff_cow.
 791  * Moreover, this function is currently only called on isl_affs
 792  * with a single reference.
 793  */
 794 static __isl_give isl_aff *sort_divs(__isl_take isl_aff *aff)
 795 {
 796         int i, j, n;
 797         unsigned off;
 798
 799         if (!aff)
 800                 return NULL;
 801
 802         off = isl_local_space_offset(aff->ls, isl_dim_div);
 803         n = isl_aff_dim(aff, isl_dim_div);
 804         for (i = 1; i < n; ++i) {
 805                 for (j = i - 1; j >= 0; --j) {
 806                         int cmp = isl_mat_cmp_div(aff->ls->div, j, j + 1);
 807                         if (cmp < 0)
 808                                 break;
 809                         if (cmp == 0)
 810                                 aff = merge_divs(aff, j, j + 1);
 811                         else
 812                                 aff = swap_div(aff, j, j + 1);
 813                         if (!aff)
 814                                 return NULL;
 815                 }
 816         }
 817
 818         return aff;
 819 }
 820
 821 /* Normalize the representation of "aff".
 822  *
 823  * This function should only be called of "new" isl_affs, i.e.,
 824  * with only a single reference.  We therefore do not need to
 825  * worry about affecting other instances.
 826  */
 827 __isl_give isl_aff *isl_aff_normalize(__isl_take isl_aff *aff)
 828 {
 829         if (!aff)
 830                 return NULL;
 831         aff->v = isl_vec_normalize(aff->v);
 832         if (!aff->v)
 833                 return isl_aff_free(aff);
 834         aff = plug_in_integral_divs(aff);
 835         aff = sort_divs(aff);
 836         aff = isl_aff_remove_unused_divs(aff);
 837         return aff;
 838 }
 839
 840 /* Given f, return floor(f).
 841  * If f is an integer expression, then just return f.
 842  * If f is a constant, then return the constant floor(f).
 843  * Otherwise, if f = g/m, write g = q m + r,
 844  * create a new div d = [r/m] and return the expression q + d.
 845  * The coefficients in r are taken to lie between -m/2 and m/2.
 846  */
 847 __isl_give isl_aff *isl_aff_floor(__isl_take isl_aff *aff)
 848 {
 849         int i;
 850         int size;
 851         isl_ctx *ctx;
 852         isl_vec *div;
 853
 854         if (!aff)
 855                 return NULL;
 856
 857         if (isl_int_is_one(aff->v->el[0]))
 858                 return aff;
 859
 860         aff = isl_aff_cow(aff);
 861         if (!aff)
 862                 return NULL;
 863
 864         aff->v = isl_vec_cow(aff->v);
 865         if (!aff->v)
 866                 return isl_aff_free(aff);
 867
 868         if (isl_aff_is_cst(aff)) {
 869                 isl_int_fdiv_q(aff->v->el[1], aff->v->el[1], aff->v->el[0]);
 870                 isl_int_set_si(aff->v->el[0], 1);
 871                 return aff;
 872         }
 873
 874         div = isl_vec_copy(aff->v);
 875         div = isl_vec_cow(div);
 876         if (!div)
 877                 return isl_aff_free(aff);
 878
 879         ctx = isl_aff_get_ctx(aff);
 880         isl_int_fdiv_q(aff->v->el[0], aff->v->el[0], ctx->two);
 881         for (i = 1; i < aff->v->size; ++i) {
 882                 isl_int_fdiv_r(div->el[i], div->el[i], div->el[0]);
 883                 isl_int_fdiv_q(aff->v->el[i], aff->v->el[i], div->el[0]);
 884                 if (isl_int_gt(div->el[i], aff->v->el[0])) {
 885                         isl_int_sub(div->el[i], div->el[i], div->el[0]);
 886                         isl_int_add_ui(aff->v->el[i], aff->v->el[i], 1);
 887                 }
 888         }
 889
 890         aff->ls = isl_local_space_add_div(aff->ls, div);
 891         if (!aff->ls)
 892                 return isl_aff_free(aff);
 893
 894         size = aff->v->size;
 895         aff->v = isl_vec_extend(aff->v, size + 1);
 896         if (!aff->v)
 897                 return isl_aff_free(aff);
 898         isl_int_set_si(aff->v->el[0], 1);
 899         isl_int_set_si(aff->v->el[size], 1);
 900
 901         return aff;
 902 }
 903
 904 /* Compute
 905  *
 906  *      aff mod m = aff - m * floor(aff/m)
 907  */
 908 __isl_give isl_aff *isl_aff_mod(__isl_take isl_aff *aff, isl_int m)
 909 {
 910         isl_aff *res;
 911
 912         res = isl_aff_copy(aff);
 913         aff = isl_aff_scale_down(aff, m);
 914         aff = isl_aff_floor(aff);
 915         aff = isl_aff_scale(aff, m);
 916         res = isl_aff_sub(res, aff);
 917
 918         return res;
 919 }
 920
 921 /* Compute
 922  *
 923  *      pwaff mod m = pwaff - m * floor(pwaff/m)
 924  */
 925 __isl_give isl_pw_aff *isl_pw_aff_mod(__isl_take isl_pw_aff *pwaff, isl_int m)
 926 {
 927         isl_pw_aff *res;
 928
 929         res = isl_pw_aff_copy(pwaff);
 930         pwaff = isl_pw_aff_scale_down(pwaff, m);
 931         pwaff = isl_pw_aff_floor(pwaff);
 932         pwaff = isl_pw_aff_scale(pwaff, m);
 933         res = isl_pw_aff_sub(res, pwaff);
 934
 935         return res;
 936 }
 937
 938 /* Given f, return ceil(f).
 939  * If f is an integer expression, then just return f.
 940  * Otherwise, create a new div d = [-f] and return the expression -d.
 941  */
 942 __isl_give isl_aff *isl_aff_ceil(__isl_take isl_aff *aff)
 943 {
 944         if (!aff)
 945                 return NULL;
 946
 947         if (isl_int_is_one(aff->v->el[0]))
 948                 return aff;
 949
 950         aff = isl_aff_neg(aff);
 951         aff = isl_aff_floor(aff);
 952         aff = isl_aff_neg(aff);
 953
 954         return aff;
 955 }
 956
 957 /* Apply the expansion computed by isl_merge_divs.
 958  * The expansion itself is given by "exp" while the resulting
 959  * list of divs is given by "div".
 960  */
 961 __isl_give isl_aff *isl_aff_expand_divs( __isl_take isl_aff *aff,
 962         __isl_take isl_mat *div, int *exp)
 963 {
 964         int i, j;
 965         int old_n_div;
 966         int new_n_div;
 967         int offset;
 968
 969         aff = isl_aff_cow(aff);
 970         if (!aff || !div)
 971                 goto error;
 972
 973         old_n_div = isl_local_space_dim(aff->ls, isl_dim_div);
 974         new_n_div = isl_mat_rows(div);
 975         if (new_n_div < old_n_div)
 976                 isl_die(isl_mat_get_ctx(div), isl_error_invalid,
 977                         "not an expansion", goto error);
 978
 979         aff->v = isl_vec_extend(aff->v, aff->v->size + new_n_div - old_n_div);
 980         if (!aff->v)
 981                 goto error;
 982
 983         offset = 1 + isl_local_space_offset(aff->ls, isl_dim_div);
 984         j = old_n_div - 1;
 985         for (i = new_n_div - 1; i >= 0; --i) {
 986                 if (j >= 0 && exp[j] == i) {
 987                         if (i != j)
 988                                 isl_int_swap(aff->v->el[offset + i],
 989                                              aff->v->el[offset + j]);
 990                         j--;
 991                 } else
 992                         isl_int_set_si(aff->v->el[offset + i], 0);
 993         }
 994
 995         aff->ls = isl_local_space_replace_divs(aff->ls, isl_mat_copy(div));
 996         if (!aff->ls)
 997                 goto error;
 998         isl_mat_free(div);
 999         return aff;
1000 error:
1001         isl_aff_free(aff);
1002         isl_mat_free(div);
1003         return NULL;
1004 }
1005
1006 /* Add two affine expressions that live in the same local space.
1007  */
1008 static __isl_give isl_aff *add_expanded(__isl_take isl_aff *aff1,
1009         __isl_take isl_aff *aff2)
1010 {
1011         isl_int gcd, f;
1012
1013         aff1 = isl_aff_cow(aff1);
1014         if (!aff1 || !aff2)
1015                 goto error;
1016
1017         aff1->v = isl_vec_cow(aff1->v);
1018         if (!aff1->v)
1019                 goto error;
1020
1021         isl_int_init(gcd);
1022         isl_int_init(f);
1023         isl_int_gcd(gcd, aff1->v->el[0], aff2->v->el[0]);
1024         isl_int_divexact(f, aff2->v->el[0], gcd);
1025         isl_seq_scale(aff1->v->el + 1, aff1->v->el + 1, f, aff1->v->size - 1);
1026         isl_int_divexact(f, aff1->v->el[0], gcd);
1027         isl_seq_addmul(aff1->v->el + 1, f, aff2->v->el + 1, aff1->v->size - 1);
1028         isl_int_divexact(f, aff2->v->el[0], gcd);
1029         isl_int_mul(aff1->v->el[0], aff1->v->el[0], f);
1030         isl_int_clear(f);
1031         isl_int_clear(gcd);
1032
1033         isl_aff_free(aff2);
1034         return aff1;
1035 error:
1036         isl_aff_free(aff1);
1037         isl_aff_free(aff2);
1038         return NULL;
1039 }
1040
1041 __isl_give isl_aff *isl_aff_add(__isl_take isl_aff *aff1,
1042         __isl_take isl_aff *aff2)
1043 {
1044         isl_ctx *ctx;
1045         int *exp1 = NULL;
1046         int *exp2 = NULL;
1047         isl_mat *div;
1048
1049         if (!aff1 || !aff2)
1050                 goto error;
1051
1052         ctx = isl_aff_get_ctx(aff1);
1053         if (!isl_space_is_equal(aff1->ls->dim, aff2->ls->dim))
1054                 isl_die(ctx, isl_error_invalid,
1055                         "spaces don't match", goto error);
1056
1057         if (aff1->ls->div->n_row == 0 && aff2->ls->div->n_row == 0)
1058                 return add_expanded(aff1, aff2);
1059
1060         exp1 = isl_alloc_array(ctx, int, aff1->ls->div->n_row);
1061         exp2 = isl_alloc_array(ctx, int, aff2->ls->div->n_row);
1062         if (!exp1 || !exp2)
1063                 goto error;
1064
1065         div = isl_merge_divs(aff1->ls->div, aff2->ls->div, exp1, exp2);
1066         aff1 = isl_aff_expand_divs(aff1, isl_mat_copy(div), exp1);
1067         aff2 = isl_aff_expand_divs(aff2, div, exp2);
1068         free(exp1);
1069         free(exp2);
1070
1071         return add_expanded(aff1, aff2);
1072 error:
1073         free(exp1);
1074         free(exp2);
1075         isl_aff_free(aff1);
1076         isl_aff_free(aff2);
1077         return NULL;
1078 }
1079
1080 __isl_give isl_aff *isl_aff_sub(__isl_take isl_aff *aff1,
1081         __isl_take isl_aff *aff2)
1082 {
1083         return isl_aff_add(aff1, isl_aff_neg(aff2));
1084 }
1085
1086 __isl_give isl_aff *isl_aff_scale(__isl_take isl_aff *aff, isl_int f)
1087 {
1088         isl_int gcd;
1089
1090         if (isl_int_is_one(f))
1091                 return aff;
1092
1093         aff = isl_aff_cow(aff);
1094         if (!aff)
1095                 return NULL;
1096         aff->v = isl_vec_cow(aff->v);
1097         if (!aff->v)
1098                 return isl_aff_free(aff);
1099
1100         isl_int_init(gcd);
1101         isl_int_gcd(gcd, aff->v->el[0], f);
1102         isl_int_divexact(aff->v->el[0], aff->v->el[0], gcd);
1103         isl_int_divexact(gcd, f, gcd);
1104         isl_seq_scale(aff->v->el + 1, aff->v->el + 1, gcd, aff->v->size - 1);
1105         isl_int_clear(gcd);
1106
1107         return aff;
1108 }
1109
1110 __isl_give isl_aff *isl_aff_scale_down(__isl_take isl_aff *aff, isl_int f)
1111 {
1112         isl_int gcd;
1113
1114         if (isl_int_is_one(f))
1115                 return aff;
1116
1117         aff = isl_aff_cow(aff);
1118         if (!aff)
1119                 return NULL;
1120
1121         if (isl_int_is_zero(f))
1122                 isl_die(isl_aff_get_ctx(aff), isl_error_invalid,
1123                         "cannot scale down by zero", return isl_aff_free(aff));
1124
1125         aff->v = isl_vec_cow(aff->v);
1126         if (!aff->v)
1127                 return isl_aff_free(aff);
1128
1129         isl_int_init(gcd);
1130         isl_seq_gcd(aff->v->el + 1, aff->v->size - 1, &gcd);
1131         isl_int_gcd(gcd, gcd, f);
1132         isl_seq_scale_down(aff->v->el + 1, aff->v->el + 1, gcd, aff->v->size - 1);
1133         isl_int_divexact(gcd, f, gcd);
1134         isl_int_mul(aff->v->el[0], aff->v->el[0], gcd);
1135         isl_int_clear(gcd);
1136
1137         return aff;
1138 }
1139
1140 __isl_give isl_aff *isl_aff_scale_down_ui(__isl_take isl_aff *aff, unsigned f)
1141 {
1142         isl_int v;
1143
1144         if (f == 1)
1145                 return aff;
1146
1147         isl_int_init(v);
1148         isl_int_set_ui(v, f);
1149         aff = isl_aff_scale_down(aff, v);
1150         isl_int_clear(v);
1151
1152         return aff;
1153 }
1154
1155 __isl_give isl_aff *isl_aff_set_dim_name(__isl_take isl_aff *aff,
1156         enum isl_dim_type type, unsigned pos, const char *s)
1157 {
1158         aff = isl_aff_cow(aff);
1159         if (!aff)
1160                 return NULL;
1161         if (type == isl_dim_out)
1162                 isl_die(aff->v->ctx, isl_error_invalid,
1163                         "cannot set name of output/set dimension",
1164                         return isl_aff_free(aff));
1165         if (type == isl_dim_in)
1166                 type = isl_dim_set;
1167         aff->ls = isl_local_space_set_dim_name(aff->ls, type, pos, s);
1168         if (!aff->ls)
1169                 return isl_aff_free(aff);
1170
1171         return aff;
1172 }
1173
1174 __isl_give isl_aff *isl_aff_set_dim_id(__isl_take isl_aff *aff,
1175         enum isl_dim_type type, unsigned pos, __isl_take isl_id *id)
1176 {
1177         aff = isl_aff_cow(aff);
1178         if (!aff)
1179                 return isl_id_free(id);
1180         if (type == isl_dim_out)
1181                 isl_die(aff->v->ctx, isl_error_invalid,
1182                         "cannot set name of output/set dimension",
1183                         goto error);
1184         if (type == isl_dim_in)
1185                 type = isl_dim_set;
1186         aff->ls = isl_local_space_set_dim_id(aff->ls, type, pos, id);
1187         if (!aff->ls)
1188                 return isl_aff_free(aff);
1189
1190         return aff;
1191 error:
1192         isl_id_free(id);
1193         isl_aff_free(aff);
1194         return NULL;
1195 }
1196
1197 /* Exploit the equalities in "eq" to simplify the affine expression
1198  * and the expressions of the integer divisions in the local space.
1199  * The integer divisions in this local space are assumed to appear
1200  * as regular dimensions in "eq".
1201  */
1202 static __isl_give isl_aff *isl_aff_substitute_equalities_lifted(
1203         __isl_take isl_aff *aff, __isl_take isl_basic_set *eq)
1204 {
1205         int i, j;
1206         unsigned total;
1207         unsigned n_div;
1208
1209         if (!eq)
1210                 goto error;
1211         if (eq->n_eq == 0) {
1212                 isl_basic_set_free(eq);
1213                 return aff;
1214         }
1215
1216         aff = isl_aff_cow(aff);
1217         if (!aff)
1218                 goto error;
1219
1220         aff->ls = isl_local_space_substitute_equalities(aff->ls,
1221                                                         isl_basic_set_copy(eq));
1222         if (!aff->ls)
1223                 goto error;
1224
1225         total = 1 + isl_space_dim(eq->dim, isl_dim_all);
1226         n_div = eq->n_div;
1227         for (i = 0; i < eq->n_eq; ++i) {
1228                 j = isl_seq_last_non_zero(eq->eq[i], total + n_div);
1229                 if (j < 0 || j == 0 || j >= total)
1230                         continue;
1231
1232                 isl_seq_elim(aff->v->el + 1, eq->eq[i], j, total,
1233                                 &aff->v->el[0]);
1234         }
1235
1236         isl_basic_set_free(eq);
1237         aff = isl_aff_normalize(aff);
1238         return aff;
1239 error:
1240         isl_basic_set_free(eq);
1241         isl_aff_free(aff);
1242         return NULL;
1243 }
1244
1245 /* Exploit the equalities in "eq" to simplify the affine expression
1246  * and the expressions of the integer divisions in the local space.
1247  */
1248 static __isl_give isl_aff *isl_aff_substitute_equalities(
1249         __isl_take isl_aff *aff, __isl_take isl_basic_set *eq)
1250 {
1251         int n_div;
1252
1253         if (!aff || !eq)
1254                 goto error;
1255         n_div = isl_local_space_dim(aff->ls, isl_dim_div);
1256         if (n_div > 0)
1257                 eq = isl_basic_set_add(eq, isl_dim_set, n_div);
1258         return isl_aff_substitute_equalities_lifted(aff, eq);
1259 error:
1260         isl_basic_set_free(eq);
1261         isl_aff_free(aff);
1262         return NULL;
1263 }
1264
1265 /* Look for equalities among the variables shared by context and aff
1266  * and the integer divisions of aff, if any.
1267  * The equalities are then used to eliminate coefficients and/or integer
1268  * divisions from aff.
1269  */
1270 __isl_give isl_aff *isl_aff_gist(__isl_take isl_aff *aff,
1271         __isl_take isl_set *context)
1272 {
1273         isl_basic_set *hull;
1274         int n_div;
1275
1276         if (!aff)
1277                 goto error;
1278         n_div = isl_local_space_dim(aff->ls, isl_dim_div);
1279         if (n_div > 0) {
1280                 isl_basic_set *bset;
1281                 isl_local_space *ls;
1282                 context = isl_set_add_dims(context, isl_dim_set, n_div);
1283                 ls = isl_aff_get_domain_local_space(aff);
1284                 bset = isl_basic_set_from_local_space(ls);
1285                 bset = isl_basic_set_lift(bset);
1286                 bset = isl_basic_set_flatten(bset);
1287                 context = isl_set_intersect(context,
1288                                             isl_set_from_basic_set(bset));
1289         }
1290
1291         hull = isl_set_affine_hull(context);
1292         return isl_aff_substitute_equalities_lifted(aff, hull);
1293 error:
1294         isl_aff_free(aff);
1295         isl_set_free(context);
1296         return NULL;
1297 }
1298
1299 __isl_give isl_aff *isl_aff_gist_params(__isl_take isl_aff *aff,
1300         __isl_take isl_set *context)
1301 {
1302         isl_set *dom_context = isl_set_universe(isl_aff_get_domain_space(aff));
1303         dom_context = isl_set_intersect_params(dom_context, context);
1304         return isl_aff_gist(aff, dom_context);
1305 }
1306
1307 /* Return a basic set containing those elements in the space
1308  * of aff where it is non-negative.
1309  */
1310 __isl_give isl_basic_set *isl_aff_nonneg_basic_set(__isl_take isl_aff *aff)
1311 {
1312         isl_constraint *ineq;
1313         isl_basic_set *bset;
1314
1315         ineq = isl_inequality_from_aff(aff);
1316
1317         bset = isl_basic_set_from_constraint(ineq);
1318         bset = isl_basic_set_simplify(bset);
1319         return bset;
1320 }
1321
1322 /* Return a basic set containing those elements in the domain space
1323  * of aff where it is negative.
1324  */
1325 __isl_give isl_basic_set *isl_aff_neg_basic_set(__isl_take isl_aff *aff)
1326 {
1327         aff = isl_aff_neg(aff);
1328         aff = isl_aff_add_constant_num_si(aff, -1);
1329         return isl_aff_nonneg_basic_set(aff);
1330 }
1331
1332 /* Return a basic set containing those elements in the space
1333  * of aff where it is zero.
1334  */
1335 __isl_give isl_basic_set *isl_aff_zero_basic_set(__isl_take isl_aff *aff)
1336 {
1337         isl_constraint *ineq;
1338         isl_basic_set *bset;
1339
1340         ineq = isl_equality_from_aff(aff);
1341
1342         bset = isl_basic_set_from_constraint(ineq);
1343         bset = isl_basic_set_simplify(bset);
1344         return bset;
1345 }
1346
1347 /* Return a basic set containing those elements in the shared space
1348  * of aff1 and aff2 where aff1 is greater than or equal to aff2.
1349  */
1350 __isl_give isl_basic_set *isl_aff_ge_basic_set(__isl_take isl_aff *aff1,
1351         __isl_take isl_aff *aff2)
1352 {
1353         aff1 = isl_aff_sub(aff1, aff2);
1354
1355         return isl_aff_nonneg_basic_set(aff1);
1356 }
1357
1358 /* Return a basic set containing those elements in the shared space
1359  * of aff1 and aff2 where aff1 is smaller than or equal to aff2.
1360  */
1361 __isl_give isl_basic_set *isl_aff_le_basic_set(__isl_take isl_aff *aff1,
1362         __isl_take isl_aff *aff2)
1363 {
1364         return isl_aff_ge_basic_set(aff2, aff1);
1365 }
1366
1367 __isl_give isl_aff *isl_aff_add_on_domain(__isl_keep isl_set *dom,
1368         __isl_take isl_aff *aff1, __isl_take isl_aff *aff2)
1369 {
1370         aff1 = isl_aff_add(aff1, aff2);
1371         aff1 = isl_aff_gist(aff1, isl_set_copy(dom));
1372         return aff1;
1373 }
1374
1375 int isl_aff_is_empty(__isl_keep isl_aff *aff)
1376 {
1377         if (!aff)
1378                 return -1;
1379
1380         return 0;
1381 }
1382
1383 /* Check whether the given affine expression has non-zero coefficient
1384  * for any dimension in the given range or if any of these dimensions
1385  * appear with non-zero coefficients in any of the integer divisions
1386  * involved in the affine expression.
1387  */
1388 int isl_aff_involves_dims(__isl_keep isl_aff *aff,
1389         enum isl_dim_type type, unsigned first, unsigned n)
1390 {
1391         int i;
1392         isl_ctx *ctx;
1393         int *active = NULL;
1394         int involves = 0;
1395
1396         if (!aff)
1397                 return -1;
1398         if (n == 0)
1399                 return 0;
1400
1401         ctx = isl_aff_get_ctx(aff);
1402         if (first + n > isl_aff_dim(aff, type))
1403                 isl_die(ctx, isl_error_invalid,
1404                         "range out of bounds", return -1);
1405
1406         active = isl_local_space_get_active(aff->ls, aff->v->el + 2);
1407         if (!active)
1408                 goto error;
1409
1410         first += isl_local_space_offset(aff->ls, type) - 1;
1411         for (i = 0; i < n; ++i)
1412                 if (active[first + i]) {
1413                         involves = 1;
1414                         break;
1415                 }
1416
1417         free(active);
1418
1419         return involves;
1420 error:
1421         free(active);
1422         return -1;
1423 }
1424
1425 __isl_give isl_aff *isl_aff_drop_dims(__isl_take isl_aff *aff,
1426         enum isl_dim_type type, unsigned first, unsigned n)
1427 {
1428         isl_ctx *ctx;
1429
1430         if (!aff)
1431                 return NULL;
1432         if (type == isl_dim_out)
1433                 isl_die(aff->v->ctx, isl_error_invalid,
1434                         "cannot drop output/set dimension",
1435                         return isl_aff_free(aff));
1436         if (type == isl_dim_in)
1437                 type = isl_dim_set;
1438         if (n == 0 && !isl_local_space_is_named_or_nested(aff->ls, type))
1439                 return aff;
1440
1441         ctx = isl_aff_get_ctx(aff);
1442         if (first + n > isl_local_space_dim(aff->ls, type))
1443                 isl_die(ctx, isl_error_invalid, "range out of bounds",
1444                         return isl_aff_free(aff));
1445
1446         aff = isl_aff_cow(aff);
1447         if (!aff)
1448                 return NULL;
1449
1450         aff->ls = isl_local_space_drop_dims(aff->ls, type, first, n);
1451         if (!aff->ls)
1452                 return isl_aff_free(aff);
1453
1454         first += 1 + isl_local_space_offset(aff->ls, type);
1455         aff->v = isl_vec_drop_els(aff->v, first, n);
1456         if (!aff->v)
1457                 return isl_aff_free(aff);
1458
1459         return aff;
1460 }
1461
1462 /* Project the domain of the affine expression onto its parameter space.
1463  * The affine expression may not involve any of the domain dimensions.
1464  */
1465 __isl_give isl_aff *isl_aff_project_domain_on_params(__isl_take isl_aff *aff)
1466 {
1467         isl_space *space;
1468         unsigned n;
1469         int involves;
1470
1471         n = isl_aff_dim(aff, isl_dim_in);
1472         involves = isl_aff_involves_dims(aff, isl_dim_in, 0, n);
1473         if (involves < 0)
1474                 return isl_aff_free(aff);
1475         if (involves)
1476                 isl_die(isl_aff_get_ctx(aff), isl_error_invalid,
1477                     "affine expression involves some of the domain dimensions",
1478                     return isl_aff_free(aff));
1479         aff = isl_aff_drop_dims(aff, isl_dim_in, 0, n);
1480         space = isl_aff_get_domain_space(aff);
1481         space = isl_space_params(space);
1482         aff = isl_aff_reset_domain_space(aff, space);
1483         return aff;
1484 }
1485
1486 __isl_give isl_aff *isl_aff_insert_dims(__isl_take isl_aff *aff,
1487         enum isl_dim_type type, unsigned first, unsigned n)
1488 {
1489         isl_ctx *ctx;
1490
1491         if (!aff)
1492                 return NULL;
1493         if (type == isl_dim_out)
1494                 isl_die(aff->v->ctx, isl_error_invalid,
1495                         "cannot insert output/set dimensions",
1496                         return isl_aff_free(aff));
1497         if (type == isl_dim_in)
1498                 type = isl_dim_set;
1499         if (n == 0 && !isl_local_space_is_named_or_nested(aff->ls, type))
1500                 return aff;
1501
1502         ctx = isl_aff_get_ctx(aff);
1503         if (first > isl_local_space_dim(aff->ls, type))
1504                 isl_die(ctx, isl_error_invalid, "position out of bounds",
1505                         return isl_aff_free(aff));
1506
1507         aff = isl_aff_cow(aff);
1508         if (!aff)
1509                 return NULL;
1510
1511         aff->ls = isl_local_space_insert_dims(aff->ls, type, first, n);
1512         if (!aff->ls)
1513                 return isl_aff_free(aff);
1514
1515         first += 1 + isl_local_space_offset(aff->ls, type);
1516         aff->v = isl_vec_insert_zero_els(aff->v, first, n);
1517         if (!aff->v)
1518                 return isl_aff_free(aff);
1519
1520         return aff;
1521 }
1522
1523 __isl_give isl_aff *isl_aff_add_dims(__isl_take isl_aff *aff,
1524         enum isl_dim_type type, unsigned n)
1525 {
1526         unsigned pos;
1527
1528         pos = isl_aff_dim(aff, type);
1529
1530         return isl_aff_insert_dims(aff, type, pos, n);
1531 }
1532
1533 __isl_give isl_pw_aff *isl_pw_aff_add_dims(__isl_take isl_pw_aff *pwaff,
1534         enum isl_dim_type type, unsigned n)
1535 {
1536         unsigned pos;
1537
1538         pos = isl_pw_aff_dim(pwaff, type);
1539
1540         return isl_pw_aff_insert_dims(pwaff, type, pos, n);
1541 }
1542
1543 __isl_give isl_pw_aff *isl_pw_aff_from_aff(__isl_take isl_aff *aff)
1544 {
1545         isl_set *dom = isl_set_universe(isl_aff_get_domain_space(aff));
1546         return isl_pw_aff_alloc(dom, aff);
1547 }
1548
1549 #undef PW
1550 #define PW isl_pw_aff
1551 #undef EL
1552 #define EL isl_aff
1553 #undef EL_IS_ZERO
1554 #define EL_IS_ZERO is_empty
1555 #undef ZERO
1556 #define ZERO empty
1557 #undef IS_ZERO
1558 #define IS_ZERO is_empty
1559 #undef FIELD
1560 #define FIELD aff
1561 #undef DEFAULT_IS_ZERO
1562 #define DEFAULT_IS_ZERO 0
1563
1564 #define NO_EVAL
1565 #define NO_OPT
1566 #define NO_MOVE_DIMS
1567 #define NO_LIFT
1568 #define NO_MORPH
1569
1570 #include <isl_pw_templ.c>
1571
1572 static __isl_give isl_set *align_params_pw_pw_set_and(
1573         __isl_take isl_pw_aff *pwaff1, __isl_take isl_pw_aff *pwaff2,
1574         __isl_give isl_set *(*fn)(__isl_take isl_pw_aff *pwaff1,
1575                                     __isl_take isl_pw_aff *pwaff2))
1576 {
1577         if (!pwaff1 || !pwaff2)
1578                 goto error;
1579         if (isl_space_match(pwaff1->dim, isl_dim_param,
1580                           pwaff2->dim, isl_dim_param))
1581                 return fn(pwaff1, pwaff2);
1582         if (!isl_space_has_named_params(pwaff1->dim) ||
1583             !isl_space_has_named_params(pwaff2->dim))
1584                 isl_die(isl_pw_aff_get_ctx(pwaff1), isl_error_invalid,
1585                         "unaligned unnamed parameters", goto error);
1586         pwaff1 = isl_pw_aff_align_params(pwaff1, isl_pw_aff_get_space(pwaff2));
1587         pwaff2 = isl_pw_aff_align_params(pwaff2, isl_pw_aff_get_space(pwaff1));
1588         return fn(pwaff1, pwaff2);
1589 error:
1590         isl_pw_aff_free(pwaff1);
1591         isl_pw_aff_free(pwaff2);
1592         return NULL;
1593 }
1594
1595 /* Compute a piecewise quasi-affine expression with a domain that
1596  * is the union of those of pwaff1 and pwaff2 and such that on each
1597  * cell, the quasi-affine expression is the better (according to cmp)
1598  * of those of pwaff1 and pwaff2.  If only one of pwaff1 or pwaff2
1599  * is defined on a given cell, then the associated expression
1600  * is the defined one.
1601  */
1602 static __isl_give isl_pw_aff *pw_aff_union_opt(__isl_take isl_pw_aff *pwaff1,
1603         __isl_take isl_pw_aff *pwaff2,
1604         __isl_give isl_basic_set *(*cmp)(__isl_take isl_aff *aff1,
1605                                         __isl_take isl_aff *aff2))
1606 {
1607         int i, j, n;
1608         isl_pw_aff *res;
1609         isl_ctx *ctx;
1610         isl_set *set;
1611
1612         if (!pwaff1 || !pwaff2)
1613                 goto error;
1614
1615         ctx = isl_space_get_ctx(pwaff1->dim);
1616         if (!isl_space_is_equal(pwaff1->dim, pwaff2->dim))
1617                 isl_die(ctx, isl_error_invalid,
1618                         "arguments should live in same space", goto error);
1619
1620         if (isl_pw_aff_is_empty(pwaff1)) {
1621                 isl_pw_aff_free(pwaff1);
1622                 return pwaff2;
1623         }
1624
1625         if (isl_pw_aff_is_empty(pwaff2)) {
1626                 isl_pw_aff_free(pwaff2);
1627                 return pwaff1;
1628         }
1629
1630         n = 2 * (pwaff1->n + 1) * (pwaff2->n + 1);
1631         res = isl_pw_aff_alloc_size(isl_space_copy(pwaff1->dim), n);
1632
1633         for (i = 0; i < pwaff1->n; ++i) {
1634                 set = isl_set_copy(pwaff1->p[i].set);
1635                 for (j = 0; j < pwaff2->n; ++j) {
1636                         struct isl_set *common;
1637                         isl_set *better;
1638
1639                         common = isl_set_intersect(
1640                                         isl_set_copy(pwaff1->p[i].set),
1641                                         isl_set_copy(pwaff2->p[j].set));
1642                         better = isl_set_from_basic_set(cmp(
1643                                         isl_aff_copy(pwaff2->p[j].aff),
1644                                         isl_aff_copy(pwaff1->p[i].aff)));
1645                         better = isl_set_intersect(common, better);
1646                         if (isl_set_plain_is_empty(better)) {
1647                                 isl_set_free(better);
1648                                 continue;
1649                         }
1650                         set = isl_set_subtract(set, isl_set_copy(better));
1651
1652                         res = isl_pw_aff_add_piece(res, better,
1653                                                 isl_aff_copy(pwaff2->p[j].aff));
1654                 }
1655                 res = isl_pw_aff_add_piece(res, set,
1656                                                 isl_aff_copy(pwaff1->p[i].aff));
1657         }
1658
1659         for (j = 0; j < pwaff2->n; ++j) {
1660                 set = isl_set_copy(pwaff2->p[j].set);
1661                 for (i = 0; i < pwaff1->n; ++i)
1662                         set = isl_set_subtract(set,
1663                                         isl_set_copy(pwaff1->p[i].set));
1664                 res = isl_pw_aff_add_piece(res, set,
1665                                                 isl_aff_copy(pwaff2->p[j].aff));
1666         }
1667
1668         isl_pw_aff_free(pwaff1);
1669         isl_pw_aff_free(pwaff2);
1670
1671         return res;
1672 error:
1673         isl_pw_aff_free(pwaff1);
1674         isl_pw_aff_free(pwaff2);
1675         return NULL;
1676 }
1677
1678 /* Compute a piecewise quasi-affine expression with a domain that
1679  * is the union of those of pwaff1 and pwaff2 and such that on each
1680  * cell, the quasi-affine expression is the maximum of those of pwaff1
1681  * and pwaff2.  If only one of pwaff1 or pwaff2 is defined on a given
1682  * cell, then the associated expression is the defined one.
1683  */
1684 static __isl_give isl_pw_aff *pw_aff_union_max(__isl_take isl_pw_aff *pwaff1,
1685         __isl_take isl_pw_aff *pwaff2)
1686 {
1687         return pw_aff_union_opt(pwaff1, pwaff2, &isl_aff_ge_basic_set);
1688 }
1689
1690 __isl_give isl_pw_aff *isl_pw_aff_union_max(__isl_take isl_pw_aff *pwaff1,
1691         __isl_take isl_pw_aff *pwaff2)
1692 {
1693         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2,
1694                                                         &pw_aff_union_max);
1695 }
1696
1697 /* Compute a piecewise quasi-affine expression with a domain that
1698  * is the union of those of pwaff1 and pwaff2 and such that on each
1699  * cell, the quasi-affine expression is the minimum of those of pwaff1
1700  * and pwaff2.  If only one of pwaff1 or pwaff2 is defined on a given
1701  * cell, then the associated expression is the defined one.
1702  */
1703 static __isl_give isl_pw_aff *pw_aff_union_min(__isl_take isl_pw_aff *pwaff1,
1704         __isl_take isl_pw_aff *pwaff2)
1705 {
1706         return pw_aff_union_opt(pwaff1, pwaff2, &isl_aff_le_basic_set);
1707 }
1708
1709 __isl_give isl_pw_aff *isl_pw_aff_union_min(__isl_take isl_pw_aff *pwaff1,
1710         __isl_take isl_pw_aff *pwaff2)
1711 {
1712         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2,
1713                                                         &pw_aff_union_min);
1714 }
1715
1716 __isl_give isl_pw_aff *isl_pw_aff_union_opt(__isl_take isl_pw_aff *pwaff1,
1717         __isl_take isl_pw_aff *pwaff2, int max)
1718 {
1719         if (max)
1720                 return isl_pw_aff_union_max(pwaff1, pwaff2);
1721         else
1722                 return isl_pw_aff_union_min(pwaff1, pwaff2);
1723 }
1724
1725 /* Construct a map with as domain the domain of pwaff and
1726  * one-dimensional range corresponding to the affine expressions.
1727  */
1728 static __isl_give isl_map *map_from_pw_aff(__isl_take isl_pw_aff *pwaff)
1729 {
1730         int i;
1731         isl_space *dim;
1732         isl_map *map;
1733
1734         if (!pwaff)
1735                 return NULL;
1736
1737         dim = isl_pw_aff_get_space(pwaff);
1738         map = isl_map_empty(dim);
1739
1740         for (i = 0; i < pwaff->n; ++i) {
1741                 isl_basic_map *bmap;
1742                 isl_map *map_i;
1743
1744                 bmap = isl_basic_map_from_aff(isl_aff_copy(pwaff->p[i].aff));
1745                 map_i = isl_map_from_basic_map(bmap);
1746                 map_i = isl_map_intersect_domain(map_i,
1747                                                 isl_set_copy(pwaff->p[i].set));
1748                 map = isl_map_union_disjoint(map, map_i);
1749         }
1750
1751         isl_pw_aff_free(pwaff);
1752
1753         return map;
1754 }
1755
1756 /* Construct a map with as domain the domain of pwaff and
1757  * one-dimensional range corresponding to the affine expressions.
1758  */
1759 __isl_give isl_map *isl_map_from_pw_aff(__isl_take isl_pw_aff *pwaff)
1760 {
1761         if (!pwaff)
1762                 return NULL;
1763         if (isl_space_is_set(pwaff->dim))
1764                 isl_die(isl_pw_aff_get_ctx(pwaff), isl_error_invalid,
1765                         "space of input is not a map",
1766                         return isl_pw_aff_free(pwaff));
1767         return map_from_pw_aff(pwaff);
1768 }
1769
1770 /* Construct a one-dimensional set with as parameter domain
1771  * the domain of pwaff and the single set dimension
1772  * corresponding to the affine expressions.
1773  */
1774 __isl_give isl_set *isl_set_from_pw_aff(__isl_take isl_pw_aff *pwaff)
1775 {
1776         if (!pwaff)
1777                 return NULL;
1778         if (!isl_space_is_set(pwaff->dim))
1779                 isl_die(isl_pw_aff_get_ctx(pwaff), isl_error_invalid,
1780                         "space of input is not a set",
1781                         return isl_pw_aff_free(pwaff));
1782         return map_from_pw_aff(pwaff);
1783 }
1784
1785 /* Return a set containing those elements in the domain
1786  * of pwaff where it is non-negative.
1787  */
1788 __isl_give isl_set *isl_pw_aff_nonneg_set(__isl_take isl_pw_aff *pwaff)
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;
1801
1802                 bset = isl_aff_nonneg_basic_set(isl_aff_copy(pwaff->p[i].aff));
1803                 set_i = isl_set_from_basic_set(bset);
1804                 set_i = isl_set_intersect(set_i, isl_set_copy(pwaff->p[i].set));
1805                 set = isl_set_union_disjoint(set, set_i);
1806         }
1807
1808         isl_pw_aff_free(pwaff);
1809
1810         return set;
1811 }
1812
1813 /* Return a set containing those elements in the domain
1814  * of pwaff where it is zero (if complement is 0) or not zero
1815  * (if complement is 1).
1816  */
1817 static __isl_give isl_set *pw_aff_zero_set(__isl_take isl_pw_aff *pwaff,
1818         int complement)
1819 {
1820         int i;
1821         isl_set *set;
1822
1823         if (!pwaff)
1824                 return NULL;
1825
1826         set = isl_set_empty(isl_pw_aff_get_domain_space(pwaff));
1827
1828         for (i = 0; i < pwaff->n; ++i) {
1829                 isl_basic_set *bset;
1830                 isl_set *set_i, *zero;
1831
1832                 bset = isl_aff_zero_basic_set(isl_aff_copy(pwaff->p[i].aff));
1833                 zero = isl_set_from_basic_set(bset);
1834                 set_i = isl_set_copy(pwaff->p[i].set);
1835                 if (complement)
1836                         set_i = isl_set_subtract(set_i, zero);
1837                 else
1838                         set_i = isl_set_intersect(set_i, zero);
1839                 set = isl_set_union_disjoint(set, set_i);
1840         }
1841
1842         isl_pw_aff_free(pwaff);
1843
1844         return set;
1845 }
1846
1847 /* Return a set containing those elements in the domain
1848  * of pwaff where it is zero.
1849  */
1850 __isl_give isl_set *isl_pw_aff_zero_set(__isl_take isl_pw_aff *pwaff)
1851 {
1852         return pw_aff_zero_set(pwaff, 0);
1853 }
1854
1855 /* Return a set containing those elements in the domain
1856  * of pwaff where it is not zero.
1857  */
1858 __isl_give isl_set *isl_pw_aff_non_zero_set(__isl_take isl_pw_aff *pwaff)
1859 {
1860         return pw_aff_zero_set(pwaff, 1);
1861 }
1862
1863 /* Return a set containing those elements in the shared domain
1864  * of pwaff1 and pwaff2 where pwaff1 is greater than (or equal) to pwaff2.
1865  *
1866  * We compute the difference on the shared domain and then construct
1867  * the set of values where this difference is non-negative.
1868  * If strict is set, we first subtract 1 from the difference.
1869  * If equal is set, we only return the elements where pwaff1 and pwaff2
1870  * are equal.
1871  */
1872 static __isl_give isl_set *pw_aff_gte_set(__isl_take isl_pw_aff *pwaff1,
1873         __isl_take isl_pw_aff *pwaff2, int strict, int equal)
1874 {
1875         isl_set *set1, *set2;
1876
1877         set1 = isl_pw_aff_domain(isl_pw_aff_copy(pwaff1));
1878         set2 = isl_pw_aff_domain(isl_pw_aff_copy(pwaff2));
1879         set1 = isl_set_intersect(set1, set2);
1880         pwaff1 = isl_pw_aff_intersect_domain(pwaff1, isl_set_copy(set1));
1881         pwaff2 = isl_pw_aff_intersect_domain(pwaff2, isl_set_copy(set1));
1882         pwaff1 = isl_pw_aff_add(pwaff1, isl_pw_aff_neg(pwaff2));
1883
1884         if (strict) {
1885                 isl_space *dim = isl_set_get_space(set1);
1886                 isl_aff *aff;
1887                 aff = isl_aff_zero_on_domain(isl_local_space_from_space(dim));
1888                 aff = isl_aff_add_constant_si(aff, -1);
1889                 pwaff1 = isl_pw_aff_add(pwaff1, isl_pw_aff_alloc(set1, aff));
1890         } else
1891                 isl_set_free(set1);
1892
1893         if (equal)
1894                 return isl_pw_aff_zero_set(pwaff1);
1895         return isl_pw_aff_nonneg_set(pwaff1);
1896 }
1897
1898 /* Return a set containing those elements in the shared domain
1899  * of pwaff1 and pwaff2 where pwaff1 is equal to pwaff2.
1900  */
1901 static __isl_give isl_set *pw_aff_eq_set(__isl_take isl_pw_aff *pwaff1,
1902         __isl_take isl_pw_aff *pwaff2)
1903 {
1904         return pw_aff_gte_set(pwaff1, pwaff2, 0, 1);
1905 }
1906
1907 __isl_give isl_set *isl_pw_aff_eq_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_eq_set);
1911 }
1912
1913 /* Return a set containing those elements in the shared domain
1914  * of pwaff1 and pwaff2 where pwaff1 is greater than or equal to pwaff2.
1915  */
1916 static __isl_give isl_set *pw_aff_ge_set(__isl_take isl_pw_aff *pwaff1,
1917         __isl_take isl_pw_aff *pwaff2)
1918 {
1919         return pw_aff_gte_set(pwaff1, pwaff2, 0, 0);
1920 }
1921
1922 __isl_give isl_set *isl_pw_aff_ge_set(__isl_take isl_pw_aff *pwaff1,
1923         __isl_take isl_pw_aff *pwaff2)
1924 {
1925         return align_params_pw_pw_set_and(pwaff1, pwaff2, &pw_aff_ge_set);
1926 }
1927
1928 /* Return a set containing those elements in the shared domain
1929  * of pwaff1 and pwaff2 where pwaff1 is strictly greater than pwaff2.
1930  */
1931 static __isl_give isl_set *pw_aff_gt_set(__isl_take isl_pw_aff *pwaff1,
1932         __isl_take isl_pw_aff *pwaff2)
1933 {
1934         return pw_aff_gte_set(pwaff1, pwaff2, 1, 0);
1935 }
1936
1937 __isl_give isl_set *isl_pw_aff_gt_set(__isl_take isl_pw_aff *pwaff1,
1938         __isl_take isl_pw_aff *pwaff2)
1939 {
1940         return align_params_pw_pw_set_and(pwaff1, pwaff2, &pw_aff_gt_set);
1941 }
1942
1943 __isl_give isl_set *isl_pw_aff_le_set(__isl_take isl_pw_aff *pwaff1,
1944         __isl_take isl_pw_aff *pwaff2)
1945 {
1946         return isl_pw_aff_ge_set(pwaff2, pwaff1);
1947 }
1948
1949 __isl_give isl_set *isl_pw_aff_lt_set(__isl_take isl_pw_aff *pwaff1,
1950         __isl_take isl_pw_aff *pwaff2)
1951 {
1952         return isl_pw_aff_gt_set(pwaff2, pwaff1);
1953 }
1954
1955 /* Return a set containing those elements in the shared domain
1956  * of the elements of list1 and list2 where each element in list1
1957  * has the relation specified by "fn" with each element in list2.
1958  */
1959 static __isl_give isl_set *pw_aff_list_set(__isl_take isl_pw_aff_list *list1,
1960         __isl_take isl_pw_aff_list *list2,
1961         __isl_give isl_set *(*fn)(__isl_take isl_pw_aff *pwaff1,
1962                                     __isl_take isl_pw_aff *pwaff2))
1963 {
1964         int i, j;
1965         isl_ctx *ctx;
1966         isl_set *set;
1967
1968         if (!list1 || !list2)
1969                 goto error;
1970
1971         ctx = isl_pw_aff_list_get_ctx(list1);
1972         if (list1->n < 1 || list2->n < 1)
1973                 isl_die(ctx, isl_error_invalid,
1974                         "list should contain at least one element", goto error);
1975
1976         set = isl_set_universe(isl_pw_aff_get_domain_space(list1->p[0]));
1977         for (i = 0; i < list1->n; ++i)
1978                 for (j = 0; j < list2->n; ++j) {
1979                         isl_set *set_ij;
1980
1981                         set_ij = fn(isl_pw_aff_copy(list1->p[i]),
1982                                     isl_pw_aff_copy(list2->p[j]));
1983                         set = isl_set_intersect(set, set_ij);
1984                 }
1985
1986         isl_pw_aff_list_free(list1);
1987         isl_pw_aff_list_free(list2);
1988         return set;
1989 error:
1990         isl_pw_aff_list_free(list1);
1991         isl_pw_aff_list_free(list2);
1992         return NULL;
1993 }
1994
1995 /* Return a set containing those elements in the shared domain
1996  * of the elements of list1 and list2 where each element in list1
1997  * is equal to each element in list2.
1998  */
1999 __isl_give isl_set *isl_pw_aff_list_eq_set(__isl_take isl_pw_aff_list *list1,
2000         __isl_take isl_pw_aff_list *list2)
2001 {
2002         return pw_aff_list_set(list1, list2, &isl_pw_aff_eq_set);
2003 }
2004
2005 __isl_give isl_set *isl_pw_aff_list_ne_set(__isl_take isl_pw_aff_list *list1,
2006         __isl_take isl_pw_aff_list *list2)
2007 {
2008         return pw_aff_list_set(list1, list2, &isl_pw_aff_ne_set);
2009 }
2010
2011 /* Return a set containing those elements in the shared domain
2012  * of the elements of list1 and list2 where each element in list1
2013  * is less than or equal to each element in list2.
2014  */
2015 __isl_give isl_set *isl_pw_aff_list_le_set(__isl_take isl_pw_aff_list *list1,
2016         __isl_take isl_pw_aff_list *list2)
2017 {
2018         return pw_aff_list_set(list1, list2, &isl_pw_aff_le_set);
2019 }
2020
2021 __isl_give isl_set *isl_pw_aff_list_lt_set(__isl_take isl_pw_aff_list *list1,
2022         __isl_take isl_pw_aff_list *list2)
2023 {
2024         return pw_aff_list_set(list1, list2, &isl_pw_aff_lt_set);
2025 }
2026
2027 __isl_give isl_set *isl_pw_aff_list_ge_set(__isl_take isl_pw_aff_list *list1,
2028         __isl_take isl_pw_aff_list *list2)
2029 {
2030         return pw_aff_list_set(list1, list2, &isl_pw_aff_ge_set);
2031 }
2032
2033 __isl_give isl_set *isl_pw_aff_list_gt_set(__isl_take isl_pw_aff_list *list1,
2034         __isl_take isl_pw_aff_list *list2)
2035 {
2036         return pw_aff_list_set(list1, list2, &isl_pw_aff_gt_set);
2037 }
2038
2039
2040 /* Return a set containing those elements in the shared domain
2041  * of pwaff1 and pwaff2 where pwaff1 is not equal to pwaff2.
2042  */
2043 static __isl_give isl_set *pw_aff_ne_set(__isl_take isl_pw_aff *pwaff1,
2044         __isl_take isl_pw_aff *pwaff2)
2045 {
2046         isl_set *set_lt, *set_gt;
2047
2048         set_lt = isl_pw_aff_lt_set(isl_pw_aff_copy(pwaff1),
2049                                    isl_pw_aff_copy(pwaff2));
2050         set_gt = isl_pw_aff_gt_set(pwaff1, pwaff2);
2051         return isl_set_union_disjoint(set_lt, set_gt);
2052 }
2053
2054 __isl_give isl_set *isl_pw_aff_ne_set(__isl_take isl_pw_aff *pwaff1,
2055         __isl_take isl_pw_aff *pwaff2)
2056 {
2057         return align_params_pw_pw_set_and(pwaff1, pwaff2, &pw_aff_ne_set);
2058 }
2059
2060 __isl_give isl_pw_aff *isl_pw_aff_scale_down(__isl_take isl_pw_aff *pwaff,
2061         isl_int v)
2062 {
2063         int i;
2064
2065         if (isl_int_is_one(v))
2066                 return pwaff;
2067         if (!isl_int_is_pos(v))
2068                 isl_die(isl_pw_aff_get_ctx(pwaff), isl_error_invalid,
2069                         "factor needs to be positive",
2070                         return isl_pw_aff_free(pwaff));
2071         pwaff = isl_pw_aff_cow(pwaff);
2072         if (!pwaff)
2073                 return NULL;
2074         if (pwaff->n == 0)
2075                 return pwaff;
2076
2077         for (i = 0; i < pwaff->n; ++i) {
2078                 pwaff->p[i].aff = isl_aff_scale_down(pwaff->p[i].aff, v);
2079                 if (!pwaff->p[i].aff)
2080                         return isl_pw_aff_free(pwaff);
2081         }
2082
2083         return pwaff;
2084 }
2085
2086 __isl_give isl_pw_aff *isl_pw_aff_floor(__isl_take isl_pw_aff *pwaff)
2087 {
2088         int i;
2089
2090         pwaff = isl_pw_aff_cow(pwaff);
2091         if (!pwaff)
2092                 return NULL;
2093         if (pwaff->n == 0)
2094                 return pwaff;
2095
2096         for (i = 0; i < pwaff->n; ++i) {
2097                 pwaff->p[i].aff = isl_aff_floor(pwaff->p[i].aff);
2098                 if (!pwaff->p[i].aff)
2099                         return isl_pw_aff_free(pwaff);
2100         }
2101
2102         return pwaff;
2103 }
2104
2105 __isl_give isl_pw_aff *isl_pw_aff_ceil(__isl_take isl_pw_aff *pwaff)
2106 {
2107         int i;
2108
2109         pwaff = isl_pw_aff_cow(pwaff);
2110         if (!pwaff)
2111                 return NULL;
2112         if (pwaff->n == 0)
2113                 return pwaff;
2114
2115         for (i = 0; i < pwaff->n; ++i) {
2116                 pwaff->p[i].aff = isl_aff_ceil(pwaff->p[i].aff);
2117                 if (!pwaff->p[i].aff)
2118                         return isl_pw_aff_free(pwaff);
2119         }
2120
2121         return pwaff;
2122 }
2123
2124 /* Assuming that "cond1" and "cond2" are disjoint,
2125  * return an affine expression that is equal to pwaff1 on cond1
2126  * and to pwaff2 on cond2.
2127  */
2128 static __isl_give isl_pw_aff *isl_pw_aff_select(
2129         __isl_take isl_set *cond1, __isl_take isl_pw_aff *pwaff1,
2130         __isl_take isl_set *cond2, __isl_take isl_pw_aff *pwaff2)
2131 {
2132         pwaff1 = isl_pw_aff_intersect_domain(pwaff1, cond1);
2133         pwaff2 = isl_pw_aff_intersect_domain(pwaff2, cond2);
2134
2135         return isl_pw_aff_add_disjoint(pwaff1, pwaff2);
2136 }
2137
2138 /* Return an affine expression that is equal to pwaff_true for elements
2139  * where "cond" is non-zero and to pwaff_false for elements where "cond"
2140  * is zero.
2141  * That is, return cond ? pwaff_true : pwaff_false;
2142  */
2143 __isl_give isl_pw_aff *isl_pw_aff_cond(__isl_take isl_pw_aff *cond,
2144         __isl_take isl_pw_aff *pwaff_true, __isl_take isl_pw_aff *pwaff_false)
2145 {
2146         isl_set *cond_true, *cond_false;
2147
2148         cond_true = isl_pw_aff_non_zero_set(isl_pw_aff_copy(cond));
2149         cond_false = isl_pw_aff_zero_set(cond);
2150         return isl_pw_aff_select(cond_true, pwaff_true,
2151                                  cond_false, pwaff_false);
2152 }
2153
2154 int isl_aff_is_cst(__isl_keep isl_aff *aff)
2155 {
2156         if (!aff)
2157                 return -1;
2158
2159         return isl_seq_first_non_zero(aff->v->el + 2, aff->v->size - 2) == -1;
2160 }
2161
2162 /* Check whether pwaff is a piecewise constant.
2163  */
2164 int isl_pw_aff_is_cst(__isl_keep isl_pw_aff *pwaff)
2165 {
2166         int i;
2167
2168         if (!pwaff)
2169                 return -1;
2170
2171         for (i = 0; i < pwaff->n; ++i) {
2172                 int is_cst = isl_aff_is_cst(pwaff->p[i].aff);
2173                 if (is_cst < 0 || !is_cst)
2174                         return is_cst;
2175         }
2176
2177         return 1;
2178 }
2179
2180 __isl_give isl_aff *isl_aff_mul(__isl_take isl_aff *aff1,
2181         __isl_take isl_aff *aff2)
2182 {
2183         if (!isl_aff_is_cst(aff2) && isl_aff_is_cst(aff1))
2184                 return isl_aff_mul(aff2, aff1);
2185
2186         if (!isl_aff_is_cst(aff2))
2187                 isl_die(isl_aff_get_ctx(aff1), isl_error_invalid,
2188                         "at least one affine expression should be constant",
2189                         goto error);
2190
2191         aff1 = isl_aff_cow(aff1);
2192         if (!aff1 || !aff2)
2193                 goto error;
2194
2195         aff1 = isl_aff_scale(aff1, aff2->v->el[1]);
2196         aff1 = isl_aff_scale_down(aff1, aff2->v->el[0]);
2197
2198         isl_aff_free(aff2);
2199         return aff1;
2200 error:
2201         isl_aff_free(aff1);
2202         isl_aff_free(aff2);
2203         return NULL;
2204 }
2205
2206 /* Divide "aff1" by "aff2", assuming "aff2" is a piecewise constant.
2207  */
2208 __isl_give isl_aff *isl_aff_div(__isl_take isl_aff *aff1,
2209         __isl_take isl_aff *aff2)
2210 {
2211         int is_cst;
2212         int neg;
2213
2214         is_cst = isl_aff_is_cst(aff2);
2215         if (is_cst < 0)
2216                 goto error;
2217         if (!is_cst)
2218                 isl_die(isl_aff_get_ctx(aff2), isl_error_invalid,
2219                         "second argument should be a constant", goto error);
2220
2221         if (!aff2)
2222                 goto error;
2223
2224         neg = isl_int_is_neg(aff2->v->el[1]);
2225         if (neg) {
2226                 isl_int_neg(aff2->v->el[0], aff2->v->el[0]);
2227                 isl_int_neg(aff2->v->el[1], aff2->v->el[1]);
2228         }
2229
2230         aff1 = isl_aff_scale(aff1, aff2->v->el[0]);
2231         aff1 = isl_aff_scale_down(aff1, aff2->v->el[1]);
2232
2233         if (neg) {
2234                 isl_int_neg(aff2->v->el[0], aff2->v->el[0]);
2235                 isl_int_neg(aff2->v->el[1], aff2->v->el[1]);
2236         }
2237
2238         isl_aff_free(aff2);
2239         return aff1;
2240 error:
2241         isl_aff_free(aff1);
2242         isl_aff_free(aff2);
2243         return NULL;
2244 }
2245
2246 static __isl_give isl_pw_aff *pw_aff_add(__isl_take isl_pw_aff *pwaff1,
2247         __isl_take isl_pw_aff *pwaff2)
2248 {
2249         return isl_pw_aff_on_shared_domain(pwaff1, pwaff2, &isl_aff_add);
2250 }
2251
2252 __isl_give isl_pw_aff *isl_pw_aff_add(__isl_take isl_pw_aff *pwaff1,
2253         __isl_take isl_pw_aff *pwaff2)
2254 {
2255         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2, &pw_aff_add);
2256 }
2257
2258 __isl_give isl_pw_aff *isl_pw_aff_union_add(__isl_take isl_pw_aff *pwaff1,
2259         __isl_take isl_pw_aff *pwaff2)
2260 {
2261         return isl_pw_aff_union_add_(pwaff1, pwaff2);
2262 }
2263
2264 static __isl_give isl_pw_aff *pw_aff_mul(__isl_take isl_pw_aff *pwaff1,
2265         __isl_take isl_pw_aff *pwaff2)
2266 {
2267         return isl_pw_aff_on_shared_domain(pwaff1, pwaff2, &isl_aff_mul);
2268 }
2269
2270 __isl_give isl_pw_aff *isl_pw_aff_mul(__isl_take isl_pw_aff *pwaff1,
2271         __isl_take isl_pw_aff *pwaff2)
2272 {
2273         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2, &pw_aff_mul);
2274 }
2275
2276 static __isl_give isl_pw_aff *pw_aff_div(__isl_take isl_pw_aff *pa1,
2277         __isl_take isl_pw_aff *pa2)
2278 {
2279         return isl_pw_aff_on_shared_domain(pa1, pa2, &isl_aff_div);
2280 }
2281
2282 /* Divide "pa1" by "pa2", assuming "pa2" is a piecewise constant.
2283  */
2284 __isl_give isl_pw_aff *isl_pw_aff_div(__isl_take isl_pw_aff *pa1,
2285         __isl_take isl_pw_aff *pa2)
2286 {
2287         int is_cst;
2288
2289         is_cst = isl_pw_aff_is_cst(pa2);
2290         if (is_cst < 0)
2291                 goto error;
2292         if (!is_cst)
2293                 isl_die(isl_pw_aff_get_ctx(pa2), isl_error_invalid,
2294                         "second argument should be a piecewise constant",
2295                         goto error);
2296         return isl_pw_aff_align_params_pw_pw_and(pa1, pa2, &pw_aff_div);
2297 error:
2298         isl_pw_aff_free(pa1);
2299         isl_pw_aff_free(pa2);
2300         return NULL;
2301 }
2302
2303 /* Compute the quotient of the integer division of "pa1" by "pa2"
2304  * with rounding towards zero.
2305  * "pa2" is assumed to be a piecewise constant.
2306  *
2307  * In particular, return
2308  *
2309  *      pa1 >= 0 ? floor(pa1/pa2) : ceil(pa1/pa2)
2310  *
2311  */
2312 __isl_give isl_pw_aff *isl_pw_aff_tdiv_q(__isl_take isl_pw_aff *pa1,
2313         __isl_take isl_pw_aff *pa2)
2314 {
2315         int is_cst;
2316         isl_set *cond;
2317         isl_pw_aff *f, *c;
2318
2319         is_cst = isl_pw_aff_is_cst(pa2);
2320         if (is_cst < 0)
2321                 goto error;
2322         if (!is_cst)
2323                 isl_die(isl_pw_aff_get_ctx(pa2), isl_error_invalid,
2324                         "second argument should be a piecewise constant",
2325                         goto error);
2326
2327         pa1 = isl_pw_aff_div(pa1, pa2);
2328
2329         cond = isl_pw_aff_nonneg_set(isl_pw_aff_copy(pa1));
2330         f = isl_pw_aff_floor(isl_pw_aff_copy(pa1));
2331         c = isl_pw_aff_ceil(pa1);
2332         return isl_pw_aff_cond(isl_set_indicator_function(cond), f, c);
2333 error:
2334         isl_pw_aff_free(pa1);
2335         isl_pw_aff_free(pa2);
2336         return NULL;
2337 }
2338
2339 /* Compute the remainder of the integer division of "pa1" by "pa2"
2340  * with rounding towards zero.
2341  * "pa2" is assumed to be a piecewise constant.
2342  *
2343  * In particular, return
2344  *
2345  *      pa1 - pa2 * (pa1 >= 0 ? floor(pa1/pa2) : ceil(pa1/pa2))
2346  *
2347  */
2348 __isl_give isl_pw_aff *isl_pw_aff_tdiv_r(__isl_take isl_pw_aff *pa1,
2349         __isl_take isl_pw_aff *pa2)
2350 {
2351         int is_cst;
2352         isl_pw_aff *res;
2353
2354         is_cst = isl_pw_aff_is_cst(pa2);
2355         if (is_cst < 0)
2356                 goto error;
2357         if (!is_cst)
2358                 isl_die(isl_pw_aff_get_ctx(pa2), isl_error_invalid,
2359                         "second argument should be a piecewise constant",
2360                         goto error);
2361         res = isl_pw_aff_tdiv_q(isl_pw_aff_copy(pa1), isl_pw_aff_copy(pa2));
2362         res = isl_pw_aff_mul(pa2, res);
2363         res = isl_pw_aff_sub(pa1, res);
2364         return res;
2365 error:
2366         isl_pw_aff_free(pa1);
2367         isl_pw_aff_free(pa2);
2368         return NULL;
2369 }
2370
2371 static __isl_give isl_pw_aff *pw_aff_min(__isl_take isl_pw_aff *pwaff1,
2372         __isl_take isl_pw_aff *pwaff2)
2373 {
2374         isl_set *le;
2375         isl_set *dom;
2376
2377         dom = isl_set_intersect(isl_pw_aff_domain(isl_pw_aff_copy(pwaff1)),
2378                                 isl_pw_aff_domain(isl_pw_aff_copy(pwaff2)));
2379         le = isl_pw_aff_le_set(isl_pw_aff_copy(pwaff1),
2380                                 isl_pw_aff_copy(pwaff2));
2381         dom = isl_set_subtract(dom, isl_set_copy(le));
2382         return isl_pw_aff_select(le, pwaff1, dom, pwaff2);
2383 }
2384
2385 __isl_give isl_pw_aff *isl_pw_aff_min(__isl_take isl_pw_aff *pwaff1,
2386         __isl_take isl_pw_aff *pwaff2)
2387 {
2388         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2, &pw_aff_min);
2389 }
2390
2391 static __isl_give isl_pw_aff *pw_aff_max(__isl_take isl_pw_aff *pwaff1,
2392         __isl_take isl_pw_aff *pwaff2)
2393 {
2394         isl_set *ge;
2395         isl_set *dom;
2396
2397         dom = isl_set_intersect(isl_pw_aff_domain(isl_pw_aff_copy(pwaff1)),
2398                                 isl_pw_aff_domain(isl_pw_aff_copy(pwaff2)));
2399         ge = isl_pw_aff_ge_set(isl_pw_aff_copy(pwaff1),
2400                                 isl_pw_aff_copy(pwaff2));
2401         dom = isl_set_subtract(dom, isl_set_copy(ge));
2402         return isl_pw_aff_select(ge, pwaff1, dom, pwaff2);
2403 }
2404
2405 __isl_give isl_pw_aff *isl_pw_aff_max(__isl_take isl_pw_aff *pwaff1,
2406         __isl_take isl_pw_aff *pwaff2)
2407 {
2408         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2, &pw_aff_max);
2409 }
2410
2411 static __isl_give isl_pw_aff *pw_aff_list_reduce(
2412         __isl_take isl_pw_aff_list *list,
2413         __isl_give isl_pw_aff *(*fn)(__isl_take isl_pw_aff *pwaff1,
2414                                         __isl_take isl_pw_aff *pwaff2))
2415 {
2416         int i;
2417         isl_ctx *ctx;
2418         isl_pw_aff *res;
2419
2420         if (!list)
2421                 return NULL;
2422
2423         ctx = isl_pw_aff_list_get_ctx(list);
2424         if (list->n < 1)
2425                 isl_die(ctx, isl_error_invalid,
2426                         "list should contain at least one element",
2427                         return isl_pw_aff_list_free(list));
2428
2429         res = isl_pw_aff_copy(list->p[0]);
2430         for (i = 1; i < list->n; ++i)
2431                 res = fn(res, isl_pw_aff_copy(list->p[i]));
2432
2433         isl_pw_aff_list_free(list);
2434         return res;
2435 }
2436
2437 /* Return an isl_pw_aff that maps each element in the intersection of the
2438  * domains of the elements of list to the minimal corresponding affine
2439  * expression.
2440  */
2441 __isl_give isl_pw_aff *isl_pw_aff_list_min(__isl_take isl_pw_aff_list *list)
2442 {
2443         return pw_aff_list_reduce(list, &isl_pw_aff_min);
2444 }
2445
2446 /* Return an isl_pw_aff that maps each element in the intersection of the
2447  * domains of the elements of list to the maximal corresponding affine
2448  * expression.
2449  */
2450 __isl_give isl_pw_aff *isl_pw_aff_list_max(__isl_take isl_pw_aff_list *list)
2451 {
2452         return pw_aff_list_reduce(list, &isl_pw_aff_max);
2453 }
2454
2455 #undef BASE
2456 #define BASE aff
2457
2458 #include <isl_multi_templ.c>
2459
2460 /* Construct an isl_multi_aff in the given space with value zero in
2461  * each of the output dimensions.
2462  */
2463 __isl_give isl_multi_aff *isl_multi_aff_zero(__isl_take isl_space *space)
2464 {
2465         int n;
2466         isl_multi_aff *ma;
2467
2468         if (!space)
2469                 return NULL;
2470
2471         n = isl_space_dim(space , isl_dim_out);
2472         ma = isl_multi_aff_alloc(isl_space_copy(space));
2473
2474         if (!n)
2475                 isl_space_free(space);
2476         else {
2477                 int i;
2478                 isl_local_space *ls;
2479                 isl_aff *aff;
2480
2481                 space = isl_space_domain(space);
2482                 ls = isl_local_space_from_space(space);
2483                 aff = isl_aff_zero_on_domain(ls);
2484
2485                 for (i = 0; i < n; ++i)
2486                         ma = isl_multi_aff_set_aff(ma, i, isl_aff_copy(aff));
2487
2488                 isl_aff_free(aff);
2489         }
2490
2491         return ma;
2492 }
2493
2494 /* Create an isl_multi_aff in the given space that maps each
2495  * input dimension to the corresponding output dimension.
2496  */
2497 __isl_give isl_multi_aff *isl_multi_aff_identity(__isl_take isl_space *space)
2498 {
2499         int n;
2500         isl_multi_aff *ma;
2501
2502         if (!space)
2503                 return NULL;
2504
2505         if (isl_space_is_set(space))
2506                 isl_die(isl_space_get_ctx(space), isl_error_invalid,
2507                         "expecting map space", goto error);
2508
2509         n = isl_space_dim(space, isl_dim_out);
2510         if (n != isl_space_dim(space, isl_dim_in))
2511                 isl_die(isl_space_get_ctx(space), isl_error_invalid,
2512                         "number of input and output dimensions needs to be "
2513                         "the same", goto error);
2514
2515         ma = isl_multi_aff_alloc(isl_space_copy(space));
2516
2517         if (!n)
2518                 isl_space_free(space);
2519         else {
2520                 int i;
2521                 isl_local_space *ls;
2522                 isl_aff *aff;
2523
2524                 space = isl_space_domain(space);
2525                 ls = isl_local_space_from_space(space);
2526                 aff = isl_aff_zero_on_domain(ls);
2527
2528                 for (i = 0; i < n; ++i) {
2529                         isl_aff *aff_i;
2530                         aff_i = isl_aff_copy(aff);
2531                         aff_i = isl_aff_add_coefficient_si(aff_i,
2532                                                             isl_dim_in, i, 1);
2533                         ma = isl_multi_aff_set_aff(ma, i, aff_i);
2534                 }
2535
2536                 isl_aff_free(aff);
2537         }
2538
2539         return ma;
2540 error:
2541         isl_space_free(space);
2542         return NULL;
2543 }
2544
2545 /* Create an isl_pw_multi_aff with the given isl_multi_aff on a universe
2546  * domain.
2547  */
2548 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_from_multi_aff(
2549         __isl_take isl_multi_aff *ma)
2550 {
2551         isl_set *dom = isl_set_universe(isl_multi_aff_get_domain_space(ma));
2552         return isl_pw_multi_aff_alloc(dom, ma);
2553 }
2554
2555 /* Create a piecewise multi-affine expression in the given space that maps each
2556  * input dimension to the corresponding output dimension.
2557  */
2558 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_identity(
2559         __isl_take isl_space *space)
2560 {
2561         return isl_pw_multi_aff_from_multi_aff(isl_multi_aff_identity(space));
2562 }
2563
2564 __isl_give isl_multi_aff *isl_multi_aff_add(__isl_take isl_multi_aff *maff1,
2565         __isl_take isl_multi_aff *maff2)
2566 {
2567         int i;
2568         isl_ctx *ctx;
2569
2570         maff1 = isl_multi_aff_cow(maff1);
2571         if (!maff1 || !maff2)
2572                 goto error;
2573
2574         ctx = isl_multi_aff_get_ctx(maff1);
2575         if (!isl_space_is_equal(maff1->space, maff2->space))
2576                 isl_die(ctx, isl_error_invalid,
2577                         "spaces don't match", goto error);
2578
2579         for (i = 0; i < maff1->n; ++i) {
2580                 maff1->p[i] = isl_aff_add(maff1->p[i],
2581                                             isl_aff_copy(maff2->p[i]));
2582                 if (!maff1->p[i])
2583                         goto error;
2584         }
2585
2586         isl_multi_aff_free(maff2);
2587         return maff1;
2588 error:
2589         isl_multi_aff_free(maff1);
2590         isl_multi_aff_free(maff2);
2591         return NULL;
2592 }
2593
2594 /* Given two multi-affine expressions A -> B and C -> D,
2595  * construct a multi-affine expression [A -> C] -> [B -> D].
2596  */
2597 __isl_give isl_multi_aff *isl_multi_aff_product(
2598         __isl_take isl_multi_aff *ma1, __isl_take isl_multi_aff *ma2)
2599 {
2600         int i;
2601         isl_aff *aff;
2602         isl_space *space;
2603         isl_multi_aff *res;
2604         int in1, in2, out1, out2;
2605
2606         in1 = isl_multi_aff_dim(ma1, isl_dim_in);
2607         in2 = isl_multi_aff_dim(ma2, isl_dim_in);
2608         out1 = isl_multi_aff_dim(ma1, isl_dim_out);
2609         out2 = isl_multi_aff_dim(ma2, isl_dim_out);
2610         space = isl_space_product(isl_multi_aff_get_space(ma1),
2611                                   isl_multi_aff_get_space(ma2));
2612         res = isl_multi_aff_alloc(isl_space_copy(space));
2613         space = isl_space_domain(space);
2614
2615         for (i = 0; i < out1; ++i) {
2616                 aff = isl_multi_aff_get_aff(ma1, i);
2617                 aff = isl_aff_insert_dims(aff, isl_dim_in, in1, in2);
2618                 aff = isl_aff_reset_domain_space(aff, isl_space_copy(space));
2619                 res = isl_multi_aff_set_aff(res, i, aff);
2620         }
2621
2622         for (i = 0; i < out2; ++i) {
2623                 aff = isl_multi_aff_get_aff(ma2, i);
2624                 aff = isl_aff_insert_dims(aff, isl_dim_in, 0, in1);
2625                 aff = isl_aff_reset_domain_space(aff, isl_space_copy(space));
2626                 res = isl_multi_aff_set_aff(res, out1 + i, aff);
2627         }
2628
2629         isl_space_free(space);
2630         isl_multi_aff_free(ma1);
2631         isl_multi_aff_free(ma2);
2632         return res;
2633 }
2634
2635 /* Exploit the equalities in "eq" to simplify the affine expressions.
2636  */
2637 static __isl_give isl_multi_aff *isl_multi_aff_substitute_equalities(
2638         __isl_take isl_multi_aff *maff, __isl_take isl_basic_set *eq)
2639 {
2640         int i;
2641
2642         maff = isl_multi_aff_cow(maff);
2643         if (!maff || !eq)
2644                 goto error;
2645
2646         for (i = 0; i < maff->n; ++i) {
2647                 maff->p[i] = isl_aff_substitute_equalities(maff->p[i],
2648                                                     isl_basic_set_copy(eq));
2649                 if (!maff->p[i])
2650                         goto error;
2651         }
2652
2653         isl_basic_set_free(eq);
2654         return maff;
2655 error:
2656         isl_basic_set_free(eq);
2657         isl_multi_aff_free(maff);
2658         return NULL;
2659 }
2660
2661 __isl_give isl_multi_aff *isl_multi_aff_scale(__isl_take isl_multi_aff *maff,
2662         isl_int f)
2663 {
2664         int i;
2665
2666         maff = isl_multi_aff_cow(maff);
2667         if (!maff)
2668                 return NULL;
2669
2670         for (i = 0; i < maff->n; ++i) {
2671                 maff->p[i] = isl_aff_scale(maff->p[i], f);
2672                 if (!maff->p[i])
2673                         return isl_multi_aff_free(maff);
2674         }
2675
2676         return maff;
2677 }
2678
2679 __isl_give isl_multi_aff *isl_multi_aff_add_on_domain(__isl_keep isl_set *dom,
2680         __isl_take isl_multi_aff *maff1, __isl_take isl_multi_aff *maff2)
2681 {
2682         maff1 = isl_multi_aff_add(maff1, maff2);
2683         maff1 = isl_multi_aff_gist(maff1, isl_set_copy(dom));
2684         return maff1;
2685 }
2686
2687 int isl_multi_aff_is_empty(__isl_keep isl_multi_aff *maff)
2688 {
2689         if (!maff)
2690                 return -1;
2691
2692         return 0;
2693 }
2694
2695 int isl_multi_aff_plain_is_equal(__isl_keep isl_multi_aff *maff1,
2696         __isl_keep isl_multi_aff *maff2)
2697 {
2698         int i;
2699         int equal;
2700
2701         if (!maff1 || !maff2)
2702                 return -1;
2703         if (maff1->n != maff2->n)
2704                 return 0;
2705         equal = isl_space_is_equal(maff1->space, maff2->space);
2706         if (equal < 0 || !equal)
2707                 return equal;
2708
2709         for (i = 0; i < maff1->n; ++i) {
2710                 equal = isl_aff_plain_is_equal(maff1->p[i], maff2->p[i]);
2711                 if (equal < 0 || !equal)
2712                         return equal;
2713         }
2714
2715         return 1;
2716 }
2717
2718 __isl_give isl_multi_aff *isl_multi_aff_set_dim_name(
2719         __isl_take isl_multi_aff *maff,
2720         enum isl_dim_type type, unsigned pos, const char *s)
2721 {
2722         int i;
2723
2724         maff = isl_multi_aff_cow(maff);
2725         if (!maff)
2726                 return NULL;
2727
2728         maff->space = isl_space_set_dim_name(maff->space, type, pos, s);
2729         if (!maff->space)
2730                 return isl_multi_aff_free(maff);
2731
2732         if (type == isl_dim_out)
2733                 return maff;
2734         for (i = 0; i < maff->n; ++i) {
2735                 maff->p[i] = isl_aff_set_dim_name(maff->p[i], type, pos, s);
2736                 if (!maff->p[i])
2737                         return isl_multi_aff_free(maff);
2738         }
2739
2740         return maff;
2741 }
2742
2743 __isl_give isl_multi_aff *isl_multi_aff_drop_dims(__isl_take isl_multi_aff *maff,
2744         enum isl_dim_type type, unsigned first, unsigned n)
2745 {
2746         int i;
2747
2748         maff = isl_multi_aff_cow(maff);
2749         if (!maff)
2750                 return NULL;
2751
2752         maff->space = isl_space_drop_dims(maff->space, type, first, n);
2753         if (!maff->space)
2754                 return isl_multi_aff_free(maff);
2755
2756         if (type == isl_dim_out) {
2757                 for (i = 0; i < n; ++i)
2758                         isl_aff_free(maff->p[first + i]);
2759                 for (i = first; i + n < maff->n; ++i)
2760                         maff->p[i] = maff->p[i + n];
2761                 maff->n -= n;
2762                 return maff;
2763         }
2764
2765         for (i = 0; i < maff->n; ++i) {
2766                 maff->p[i] = isl_aff_drop_dims(maff->p[i], type, first, n);
2767                 if (!maff->p[i])
2768                         return isl_multi_aff_free(maff);
2769         }
2770
2771         return maff;
2772 }
2773
2774 /* Return the set of domain elements where "ma1" is lexicographically
2775  * smaller than or equal to "ma2".
2776  */
2777 __isl_give isl_set *isl_multi_aff_lex_le_set(__isl_take isl_multi_aff *ma1,
2778         __isl_take isl_multi_aff *ma2)
2779 {
2780         return isl_multi_aff_lex_ge_set(ma2, ma1);
2781 }
2782
2783 /* Return the set of domain elements where "ma1" is lexicographically
2784  * greater than or equal to "ma2".
2785  */
2786 __isl_give isl_set *isl_multi_aff_lex_ge_set(__isl_take isl_multi_aff *ma1,
2787         __isl_take isl_multi_aff *ma2)
2788 {
2789         isl_space *space;
2790         isl_map *map1, *map2;
2791         isl_map *map, *ge;
2792
2793         map1 = isl_map_from_multi_aff(ma1);
2794         map2 = isl_map_from_multi_aff(ma2);
2795         map = isl_map_range_product(map1, map2);
2796         space = isl_space_range(isl_map_get_space(map));
2797         space = isl_space_domain(isl_space_unwrap(space));
2798         ge = isl_map_lex_ge(space);
2799         map = isl_map_intersect_range(map, isl_map_wrap(ge));
2800
2801         return isl_map_domain(map);
2802 }
2803
2804 #undef PW
2805 #define PW isl_pw_multi_aff
2806 #undef EL
2807 #define EL isl_multi_aff
2808 #undef EL_IS_ZERO
2809 #define EL_IS_ZERO is_empty
2810 #undef ZERO
2811 #define ZERO empty
2812 #undef IS_ZERO
2813 #define IS_ZERO is_empty
2814 #undef FIELD
2815 #define FIELD maff
2816 #undef DEFAULT_IS_ZERO
2817 #define DEFAULT_IS_ZERO 0
2818
2819 #define NO_NEG
2820 #define NO_EVAL
2821 #define NO_OPT
2822 #define NO_INVOLVES_DIMS
2823 #define NO_MOVE_DIMS
2824 #define NO_INSERT_DIMS
2825 #define NO_LIFT
2826 #define NO_MORPH
2827
2828 #include <isl_pw_templ.c>
2829
2830 #undef UNION
2831 #define UNION isl_union_pw_multi_aff
2832 #undef PART
2833 #define PART isl_pw_multi_aff
2834 #undef PARTS
2835 #define PARTS pw_multi_aff
2836 #define ALIGN_DOMAIN
2837
2838 #define NO_EVAL
2839
2840 #include <isl_union_templ.c>
2841
2842 /* Given a function "cmp" that returns the set of elements where
2843  * "ma1" is "better" than "ma2", return the intersection of this
2844  * set with "dom1" and "dom2".
2845  */
2846 static __isl_give isl_set *shared_and_better(__isl_keep isl_set *dom1,
2847         __isl_keep isl_set *dom2, __isl_keep isl_multi_aff *ma1,
2848         __isl_keep isl_multi_aff *ma2,
2849         __isl_give isl_set *(*cmp)(__isl_take isl_multi_aff *ma1,
2850                                     __isl_take isl_multi_aff *ma2))
2851 {
2852         isl_set *common;
2853         isl_set *better;
2854         int is_empty;
2855
2856         common = isl_set_intersect(isl_set_copy(dom1), isl_set_copy(dom2));
2857         is_empty = isl_set_plain_is_empty(common);
2858         if (is_empty >= 0 && is_empty)
2859                 return common;
2860         if (is_empty < 0)
2861                 return isl_set_free(common);
2862         better = cmp(isl_multi_aff_copy(ma1), isl_multi_aff_copy(ma2));
2863         better = isl_set_intersect(common, better);
2864
2865         return better;
2866 }
2867
2868 /* Given a function "cmp" that returns the set of elements where
2869  * "ma1" is "better" than "ma2", return a piecewise multi affine
2870  * expression defined on the union of the definition domains
2871  * of "pma1" and "pma2" that maps to the "best" of "pma1" and
2872  * "pma2" on each cell.  If only one of the two input functions
2873  * is defined on a given cell, then it is considered the best.
2874  */
2875 static __isl_give isl_pw_multi_aff *pw_multi_aff_union_opt(
2876         __isl_take isl_pw_multi_aff *pma1,
2877         __isl_take isl_pw_multi_aff *pma2,
2878         __isl_give isl_set *(*cmp)(__isl_take isl_multi_aff *ma1,
2879                                     __isl_take isl_multi_aff *ma2))
2880 {
2881         int i, j, n;
2882         isl_pw_multi_aff *res = NULL;
2883         isl_ctx *ctx;
2884         isl_set *set = NULL;
2885
2886         if (!pma1 || !pma2)
2887                 goto error;
2888
2889         ctx = isl_space_get_ctx(pma1->dim);
2890         if (!isl_space_is_equal(pma1->dim, pma2->dim))
2891                 isl_die(ctx, isl_error_invalid,
2892                         "arguments should live in the same space", goto error);
2893
2894         if (isl_pw_multi_aff_is_empty(pma1)) {
2895                 isl_pw_multi_aff_free(pma1);
2896                 return pma2;
2897         }
2898
2899         if (isl_pw_multi_aff_is_empty(pma2)) {
2900                 isl_pw_multi_aff_free(pma2);
2901                 return pma1;
2902         }
2903
2904         n = 2 * (pma1->n + 1) * (pma2->n + 1);
2905         res = isl_pw_multi_aff_alloc_size(isl_space_copy(pma1->dim), n);
2906
2907         for (i = 0; i < pma1->n; ++i) {
2908                 set = isl_set_copy(pma1->p[i].set);
2909                 for (j = 0; j < pma2->n; ++j) {
2910                         isl_set *better;
2911                         int is_empty;
2912
2913                         better = shared_and_better(pma2->p[j].set,
2914                                         pma1->p[i].set, pma2->p[j].maff,
2915                                         pma1->p[i].maff, cmp);
2916                         is_empty = isl_set_plain_is_empty(better);
2917                         if (is_empty < 0 || is_empty) {
2918                                 isl_set_free(better);
2919                                 if (is_empty < 0)
2920                                         goto error;
2921                                 continue;
2922                         }
2923                         set = isl_set_subtract(set, isl_set_copy(better));
2924
2925                         res = isl_pw_multi_aff_add_piece(res, better,
2926                                         isl_multi_aff_copy(pma2->p[j].maff));
2927                 }
2928                 res = isl_pw_multi_aff_add_piece(res, set,
2929                                         isl_multi_aff_copy(pma1->p[i].maff));
2930         }
2931
2932         for (j = 0; j < pma2->n; ++j) {
2933                 set = isl_set_copy(pma2->p[j].set);
2934                 for (i = 0; i < pma1->n; ++i)
2935                         set = isl_set_subtract(set,
2936                                         isl_set_copy(pma1->p[i].set));
2937                 res = isl_pw_multi_aff_add_piece(res, set,
2938                                         isl_multi_aff_copy(pma2->p[j].maff));
2939         }
2940
2941         isl_pw_multi_aff_free(pma1);
2942         isl_pw_multi_aff_free(pma2);
2943
2944         return res;
2945 error:
2946         isl_pw_multi_aff_free(pma1);
2947         isl_pw_multi_aff_free(pma2);
2948         isl_set_free(set);
2949         return isl_pw_multi_aff_free(res);
2950 }
2951
2952 static __isl_give isl_pw_multi_aff *pw_multi_aff_union_lexmax(
2953         __isl_take isl_pw_multi_aff *pma1,
2954         __isl_take isl_pw_multi_aff *pma2)
2955 {
2956         return pw_multi_aff_union_opt(pma1, pma2, &isl_multi_aff_lex_ge_set);
2957 }
2958
2959 /* Given two piecewise multi affine expressions, return a piecewise
2960  * multi-affine expression defined on the union of the definition domains
2961  * of the inputs that is equal to the lexicographic maximum of the two
2962  * inputs on each cell.  If only one of the two inputs is defined on
2963  * a given cell, then it is considered to be the maximum.
2964  */
2965 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_union_lexmax(
2966         __isl_take isl_pw_multi_aff *pma1,
2967         __isl_take isl_pw_multi_aff *pma2)
2968 {
2969         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
2970                                                     &pw_multi_aff_union_lexmax);
2971 }
2972
2973 static __isl_give isl_pw_multi_aff *pw_multi_aff_union_lexmin(
2974         __isl_take isl_pw_multi_aff *pma1,
2975         __isl_take isl_pw_multi_aff *pma2)
2976 {
2977         return pw_multi_aff_union_opt(pma1, pma2, &isl_multi_aff_lex_le_set);
2978 }
2979
2980 /* Given two piecewise multi affine expressions, return a piecewise
2981  * multi-affine expression defined on the union of the definition domains
2982  * of the inputs that is equal to the lexicographic minimum of the two
2983  * inputs on each cell.  If only one of the two inputs is defined on
2984  * a given cell, then it is considered to be the minimum.
2985  */
2986 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_union_lexmin(
2987         __isl_take isl_pw_multi_aff *pma1,
2988         __isl_take isl_pw_multi_aff *pma2)
2989 {
2990         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
2991                                                     &pw_multi_aff_union_lexmin);
2992 }
2993
2994 static __isl_give isl_pw_multi_aff *pw_multi_aff_add(
2995         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
2996 {
2997         return isl_pw_multi_aff_on_shared_domain(pma1, pma2,
2998                                                 &isl_multi_aff_add);
2999 }
3000
3001 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_add(
3002         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
3003 {
3004         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
3005                                                 &pw_multi_aff_add);
3006 }
3007
3008 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_union_add(
3009         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
3010 {
3011         return isl_pw_multi_aff_union_add_(pma1, pma2);
3012 }
3013
3014 /* Given two piecewise multi-affine expressions A -> B and C -> D,
3015  * construct a piecewise multi-affine expression [A -> C] -> [B -> D].
3016  */
3017 static __isl_give isl_pw_multi_aff *pw_multi_aff_product(
3018         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
3019 {
3020         int i, j, n;
3021         isl_space *space;
3022         isl_pw_multi_aff *res;
3023
3024         if (!pma1 || !pma2)
3025                 goto error;
3026
3027         n = pma1->n * pma2->n;
3028         space = isl_space_product(isl_space_copy(pma1->dim),
3029                                   isl_space_copy(pma2->dim));
3030         res = isl_pw_multi_aff_alloc_size(space, n);
3031
3032         for (i = 0; i < pma1->n; ++i) {
3033                 for (j = 0; j < pma2->n; ++j) {
3034                         isl_set *domain;
3035                         isl_multi_aff *ma;
3036
3037                         domain = isl_set_product(isl_set_copy(pma1->p[i].set),
3038                                                  isl_set_copy(pma2->p[j].set));
3039                         ma = isl_multi_aff_product(
3040                                         isl_multi_aff_copy(pma1->p[i].maff),
3041                                         isl_multi_aff_copy(pma2->p[i].maff));
3042                         res = isl_pw_multi_aff_add_piece(res, domain, ma);
3043                 }
3044         }
3045
3046         isl_pw_multi_aff_free(pma1);
3047         isl_pw_multi_aff_free(pma2);
3048         return res;
3049 error:
3050         isl_pw_multi_aff_free(pma1);
3051         isl_pw_multi_aff_free(pma2);
3052         return NULL;
3053 }
3054
3055 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_product(
3056         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
3057 {
3058         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
3059                                                 &pw_multi_aff_product);
3060 }
3061
3062 /* Construct a map mapping the domain of the piecewise multi-affine expression
3063  * to its range, with each dimension in the range equated to the
3064  * corresponding affine expression on its cell.
3065  */
3066 __isl_give isl_map *isl_map_from_pw_multi_aff(__isl_take isl_pw_multi_aff *pma)
3067 {
3068         int i;
3069         isl_map *map;
3070
3071         if (!pma)
3072                 return NULL;
3073
3074         map = isl_map_empty(isl_pw_multi_aff_get_space(pma));
3075
3076         for (i = 0; i < pma->n; ++i) {
3077                 isl_multi_aff *maff;
3078                 isl_basic_map *bmap;
3079                 isl_map *map_i;
3080
3081                 maff = isl_multi_aff_copy(pma->p[i].maff);
3082                 bmap = isl_basic_map_from_multi_aff(maff);
3083                 map_i = isl_map_from_basic_map(bmap);
3084                 map_i = isl_map_intersect_domain(map_i,
3085                                                 isl_set_copy(pma->p[i].set));
3086                 map = isl_map_union_disjoint(map, map_i);
3087         }
3088
3089         isl_pw_multi_aff_free(pma);
3090         return map;
3091 }
3092
3093 __isl_give isl_set *isl_set_from_pw_multi_aff(__isl_take isl_pw_multi_aff *pma)
3094 {
3095         if (!pma)
3096                 return NULL;
3097
3098         if (!isl_space_is_set(pma->dim))
3099                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
3100                         "isl_pw_multi_aff cannot be converted into an isl_set",
3101                         return isl_pw_multi_aff_free(pma));
3102
3103         return isl_map_from_pw_multi_aff(pma);
3104 }
3105
3106 /* Given a basic map with a single output dimension that is defined
3107  * in terms of the parameters and input dimensions using an equality,
3108  * extract an isl_aff that expresses the output dimension in terms
3109  * of the parameters and input dimensions.
3110  *
3111  * Since some applications expect the result of isl_pw_multi_aff_from_map
3112  * to only contain integer affine expressions, we compute the floor
3113  * of the expression before returning.
3114  *
3115  * This function shares some similarities with
3116  * isl_basic_map_has_defining_equality and isl_constraint_get_bound.
3117  */
3118 static __isl_give isl_aff *extract_isl_aff_from_basic_map(
3119         __isl_take isl_basic_map *bmap)
3120 {
3121         int i;
3122         unsigned offset;
3123         unsigned total;
3124         isl_local_space *ls;
3125         isl_aff *aff;
3126
3127         if (!bmap)
3128                 return NULL;
3129         if (isl_basic_map_dim(bmap, isl_dim_out) != 1)
3130                 isl_die(isl_basic_map_get_ctx(bmap), isl_error_invalid,
3131                         "basic map should have a single output dimension",
3132                         goto error);
3133         offset = isl_basic_map_offset(bmap, isl_dim_out);
3134         total = isl_basic_map_total_dim(bmap);
3135         for (i = 0; i < bmap->n_eq; ++i) {
3136                 if (isl_int_is_zero(bmap->eq[i][offset]))
3137                         continue;
3138                 if (isl_seq_first_non_zero(bmap->eq[i] + offset + 1,
3139                                            1 + total - (offset + 1)) != -1)
3140                         continue;
3141                 break;
3142         }
3143         if (i >= bmap->n_eq)
3144                 isl_die(isl_basic_map_get_ctx(bmap), isl_error_invalid,
3145                         "unable to find suitable equality", goto error);
3146         ls = isl_basic_map_get_local_space(bmap);
3147         aff = isl_aff_alloc(isl_local_space_domain(ls));
3148         if (!aff)
3149                 goto error;
3150         if (isl_int_is_neg(bmap->eq[i][offset]))
3151                 isl_seq_cpy(aff->v->el + 1, bmap->eq[i], offset);
3152         else
3153                 isl_seq_neg(aff->v->el + 1, bmap->eq[i], offset);
3154         isl_seq_clr(aff->v->el + 1 + offset, aff->v->size - (1 + offset));
3155         isl_int_abs(aff->v->el[0], bmap->eq[i][offset]);
3156         isl_basic_map_free(bmap);
3157
3158         aff = isl_aff_remove_unused_divs(aff);
3159         aff = isl_aff_floor(aff);
3160         return aff;
3161 error:
3162         isl_basic_map_free(bmap);
3163         return NULL;
3164 }
3165
3166 /* Given a basic map where each output dimension is defined
3167  * in terms of the parameters and input dimensions using an equality,
3168  * extract an isl_multi_aff that expresses the output dimensions in terms
3169  * of the parameters and input dimensions.
3170  */
3171 static __isl_give isl_multi_aff *extract_isl_multi_aff_from_basic_map(
3172         __isl_take isl_basic_map *bmap)
3173 {
3174         int i;
3175         unsigned n_out;
3176         isl_multi_aff *ma;
3177
3178         if (!bmap)
3179                 return NULL;
3180
3181         ma = isl_multi_aff_alloc(isl_basic_map_get_space(bmap));
3182         n_out = isl_basic_map_dim(bmap, isl_dim_out);
3183
3184         for (i = 0; i < n_out; ++i) {
3185                 isl_basic_map *bmap_i;
3186                 isl_aff *aff;
3187
3188                 bmap_i = isl_basic_map_copy(bmap);
3189                 bmap_i = isl_basic_map_project_out(bmap_i, isl_dim_out,
3190                                                         i + 1, n_out - (1 + i));
3191                 bmap_i = isl_basic_map_project_out(bmap_i, isl_dim_out, 0, i);
3192                 aff = extract_isl_aff_from_basic_map(bmap_i);
3193                 ma = isl_multi_aff_set_aff(ma, i, aff);
3194         }
3195
3196         isl_basic_map_free(bmap);
3197
3198         return ma;
3199 }
3200
3201 /* Create an isl_pw_multi_aff that is equivalent to
3202  * isl_map_intersect_domain(isl_map_from_basic_map(bmap), domain).
3203  * The given basic map is such that each output dimension is defined
3204  * in terms of the parameters and input dimensions using an equality.
3205  */
3206 static __isl_give isl_pw_multi_aff *plain_pw_multi_aff_from_map(
3207         __isl_take isl_set *domain, __isl_take isl_basic_map *bmap)
3208 {
3209         isl_multi_aff *ma;
3210
3211         ma = extract_isl_multi_aff_from_basic_map(bmap);
3212         return isl_pw_multi_aff_alloc(domain, ma);
3213 }
3214
3215 /* Try and create an isl_pw_multi_aff that is equivalent to the given isl_map.
3216  * This obviously only works if the input "map" is single-valued.
3217  * If so, we compute the lexicographic minimum of the image in the form
3218  * of an isl_pw_multi_aff.  Since the image is unique, it is equal
3219  * to its lexicographic minimum.
3220  * If the input is not single-valued, we produce an error.
3221  *
3222  * As a special case, we first check if all output dimensions are uniquely
3223  * defined in terms of the parameters and input dimensions over the entire
3224  * domain.  If so, we extract the desired isl_pw_multi_aff directly
3225  * from the affine hull of "map" and its domain.
3226  */
3227 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_from_map(__isl_take isl_map *map)
3228 {
3229         int i;
3230         int sv;
3231         isl_pw_multi_aff *pma;
3232         isl_basic_map *hull;
3233
3234         if (!map)
3235                 return NULL;
3236
3237         hull = isl_map_affine_hull(isl_map_copy(map));
3238         sv = isl_basic_map_plain_is_single_valued(hull);
3239         if (sv >= 0 && sv)
3240                 return plain_pw_multi_aff_from_map(isl_map_domain(map), hull);
3241         isl_basic_map_free(hull);
3242         if (sv < 0)
3243                 goto error;
3244
3245         sv = isl_map_is_single_valued(map);
3246         if (sv < 0)
3247                 goto error;
3248         if (!sv)
3249                 isl_die(isl_map_get_ctx(map), isl_error_invalid,
3250                         "map is not single-valued", goto error);
3251         map = isl_map_make_disjoint(map);
3252         if (!map)
3253                 return NULL;
3254
3255         pma = isl_pw_multi_aff_empty(isl_map_get_space(map));
3256
3257         for (i = 0; i < map->n; ++i) {
3258                 isl_pw_multi_aff *pma_i;
3259                 isl_basic_map *bmap;
3260                 bmap = isl_basic_map_copy(map->p[i]);
3261                 pma_i = isl_basic_map_lexmin_pw_multi_aff(bmap);
3262                 pma = isl_pw_multi_aff_add_disjoint(pma, pma_i);
3263         }
3264
3265         isl_map_free(map);
3266         return pma;
3267 error:
3268         isl_map_free(map);
3269         return NULL;
3270 }
3271
3272 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_from_set(__isl_take isl_set *set)
3273 {
3274         return isl_pw_multi_aff_from_map(set);
3275 }
3276
3277 /* Return the piecewise affine expression "set ? 1 : 0".
3278  */
3279 __isl_give isl_pw_aff *isl_set_indicator_function(__isl_take isl_set *set)
3280 {
3281         isl_pw_aff *pa;
3282         isl_space *space = isl_set_get_space(set);
3283         isl_local_space *ls = isl_local_space_from_space(space);
3284         isl_aff *zero = isl_aff_zero_on_domain(isl_local_space_copy(ls));
3285         isl_aff *one = isl_aff_zero_on_domain(ls);
3286
3287         one = isl_aff_add_constant_si(one, 1);
3288         pa = isl_pw_aff_alloc(isl_set_copy(set), one);
3289         set = isl_set_complement(set);
3290         pa = isl_pw_aff_add_disjoint(pa, isl_pw_aff_alloc(set, zero));
3291
3292         return pa;
3293 }
3294
3295 /* Plug in "subs" for dimension "type", "pos" of "aff".
3296  *
3297  * Let i be the dimension to replace and let "subs" be of the form
3298  *
3299  *      f/d
3300  *
3301  * and "aff" of the form
3302  *
3303  *      (a i + g)/m
3304  *
3305  * The result is
3306  *
3307  *      (a f + d g')/(m d)
3308  *
3309  * where g' is the result of plugging in "subs" in each of the integer
3310  * divisions in g.
3311  */
3312 __isl_give isl_aff *isl_aff_substitute(__isl_take isl_aff *aff,
3313         enum isl_dim_type type, unsigned pos, __isl_keep isl_aff *subs)
3314 {
3315         isl_ctx *ctx;
3316         isl_int v;
3317
3318         aff = isl_aff_cow(aff);
3319         if (!aff || !subs)
3320                 return isl_aff_free(aff);
3321
3322         ctx = isl_aff_get_ctx(aff);
3323         if (!isl_space_is_equal(aff->ls->dim, subs->ls->dim))
3324                 isl_die(ctx, isl_error_invalid,
3325                         "spaces don't match", return isl_aff_free(aff));
3326         if (isl_local_space_dim(subs->ls, isl_dim_div) != 0)
3327                 isl_die(ctx, isl_error_unsupported,
3328                         "cannot handle divs yet", return isl_aff_free(aff));
3329
3330         aff->ls = isl_local_space_substitute(aff->ls, type, pos, subs);
3331         if (!aff->ls)
3332                 return isl_aff_free(aff);
3333
3334         aff->v = isl_vec_cow(aff->v);
3335         if (!aff->v)
3336                 return isl_aff_free(aff);
3337
3338         pos += isl_local_space_offset(aff->ls, type);
3339
3340         isl_int_init(v);
3341         isl_seq_substitute(aff->v->el, pos, subs->v->el,
3342                             aff->v->size, subs->v->size, v);
3343         isl_int_clear(v);
3344
3345         return aff;
3346 }
3347
3348 /* Plug in "subs" for dimension "type", "pos" in each of the affine
3349  * expressions in "maff".
3350  */
3351 __isl_give isl_multi_aff *isl_multi_aff_substitute(
3352         __isl_take isl_multi_aff *maff, enum isl_dim_type type, unsigned pos,
3353         __isl_keep isl_aff *subs)
3354 {
3355         int i;
3356
3357         maff = isl_multi_aff_cow(maff);
3358         if (!maff || !subs)
3359                 return isl_multi_aff_free(maff);
3360
3361         if (type == isl_dim_in)
3362                 type = isl_dim_set;
3363
3364         for (i = 0; i < maff->n; ++i) {
3365                 maff->p[i] = isl_aff_substitute(maff->p[i], type, pos, subs);
3366                 if (!maff->p[i])
3367                         return isl_multi_aff_free(maff);
3368         }
3369
3370         return maff;
3371 }
3372
3373 /* Plug in "subs" for dimension "type", "pos" of "pma".
3374  *
3375  * pma is of the form
3376  *
3377  *      A_i(v) -> M_i(v)
3378  *
3379  * while subs is of the form
3380  *
3381  *      v' = B_j(v) -> S_j
3382  *
3383  * Each pair i,j such that C_ij = A_i \cap B_i is non-empty
3384  * has a contribution in the result, in particular
3385  *
3386  *      C_ij(S_j) -> M_i(S_j)
3387  *
3388  * Note that plugging in S_j in C_ij may also result in an empty set
3389  * and this contribution should simply be discarded.
3390  */
3391 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_substitute(
3392         __isl_take isl_pw_multi_aff *pma, enum isl_dim_type type, unsigned pos,
3393         __isl_keep isl_pw_aff *subs)
3394 {
3395         int i, j, n;
3396         isl_pw_multi_aff *res;
3397
3398         if (!pma || !subs)
3399                 return isl_pw_multi_aff_free(pma);
3400
3401         n = pma->n * subs->n;
3402         res = isl_pw_multi_aff_alloc_size(isl_space_copy(pma->dim), n);
3403
3404         for (i = 0; i < pma->n; ++i) {
3405                 for (j = 0; j < subs->n; ++j) {
3406                         isl_set *common;
3407                         isl_multi_aff *res_ij;
3408                         common = isl_set_intersect(
3409                                         isl_set_copy(pma->p[i].set),
3410                                         isl_set_copy(subs->p[j].set));
3411                         common = isl_set_substitute(common,
3412                                         type, pos, subs->p[j].aff);
3413                         if (isl_set_plain_is_empty(common)) {
3414                                 isl_set_free(common);
3415                                 continue;
3416                         }
3417
3418                         res_ij = isl_multi_aff_substitute(
3419                                         isl_multi_aff_copy(pma->p[i].maff),
3420                                         type, pos, subs->p[j].aff);
3421
3422                         res = isl_pw_multi_aff_add_piece(res, common, res_ij);
3423                 }
3424         }
3425
3426         isl_pw_multi_aff_free(pma);
3427         return res;
3428 }
3429
3430 /* Extend the local space of "dst" to include the divs
3431  * in the local space of "src".
3432  */
3433 __isl_give isl_aff *isl_aff_align_divs(__isl_take isl_aff *dst,
3434         __isl_keep isl_aff *src)
3435 {
3436         isl_ctx *ctx;
3437         int *exp1 = NULL;
3438         int *exp2 = NULL;
3439         isl_mat *div;
3440
3441         if (!src || !dst)
3442                 return isl_aff_free(dst);
3443
3444         ctx = isl_aff_get_ctx(src);
3445         if (!isl_space_is_equal(src->ls->dim, dst->ls->dim))
3446                 isl_die(ctx, isl_error_invalid,
3447                         "spaces don't match", goto error);
3448
3449         if (src->ls->div->n_row == 0)
3450                 return dst;
3451
3452         exp1 = isl_alloc_array(ctx, int, src->ls->div->n_row);
3453         exp2 = isl_alloc_array(ctx, int, dst->ls->div->n_row);
3454         if (!exp1 || !exp2)
3455                 goto error;
3456
3457         div = isl_merge_divs(src->ls->div, dst->ls->div, exp1, exp2);
3458         dst = isl_aff_expand_divs(dst, div, exp2);
3459         free(exp1);
3460         free(exp2);
3461
3462         return dst;
3463 error:
3464         free(exp1);
3465         free(exp2);
3466         return isl_aff_free(dst);
3467 }
3468
3469 /* Adjust the local spaces of the affine expressions in "maff"
3470  * such that they all have the save divs.
3471  */
3472 __isl_give isl_multi_aff *isl_multi_aff_align_divs(
3473         __isl_take isl_multi_aff *maff)
3474 {
3475         int i;
3476
3477         if (!maff)
3478                 return NULL;
3479         if (maff->n == 0)
3480                 return maff;
3481         maff = isl_multi_aff_cow(maff);
3482         if (!maff)
3483                 return NULL;
3484
3485         for (i = 1; i < maff->n; ++i)
3486                 maff->p[0] = isl_aff_align_divs(maff->p[0], maff->p[i]);
3487         for (i = 1; i < maff->n; ++i) {
3488                 maff->p[i] = isl_aff_align_divs(maff->p[i], maff->p[0]);
3489                 if (!maff->p[i])
3490                         return isl_multi_aff_free(maff);
3491         }
3492
3493         return maff;
3494 }
3495
3496 __isl_give isl_aff *isl_aff_lift(__isl_take isl_aff *aff)
3497 {
3498         aff = isl_aff_cow(aff);
3499         if (!aff)
3500                 return NULL;
3501
3502         aff->ls = isl_local_space_lift(aff->ls);
3503         if (!aff->ls)
3504                 return isl_aff_free(aff);
3505
3506         return aff;
3507 }
3508
3509 /* Lift "maff" to a space with extra dimensions such that the result
3510  * has no more existentially quantified variables.
3511  * If "ls" is not NULL, then *ls is assigned the local space that lies
3512  * at the basis of the lifting applied to "maff".
3513  */
3514 __isl_give isl_multi_aff *isl_multi_aff_lift(__isl_take isl_multi_aff *maff,
3515         __isl_give isl_local_space **ls)
3516 {
3517         int i;
3518         isl_space *space;
3519         unsigned n_div;
3520
3521         if (ls)
3522                 *ls = NULL;
3523
3524         if (!maff)
3525                 return NULL;
3526
3527         if (maff->n == 0) {
3528                 if (ls) {
3529                         isl_space *space = isl_multi_aff_get_domain_space(maff);
3530                         *ls = isl_local_space_from_space(space);
3531                         if (!*ls)
3532                                 return isl_multi_aff_free(maff);
3533                 }
3534                 return maff;
3535         }
3536
3537         maff = isl_multi_aff_cow(maff);
3538         maff = isl_multi_aff_align_divs(maff);
3539         if (!maff)
3540                 return NULL;
3541
3542         n_div = isl_aff_dim(maff->p[0], isl_dim_div);
3543         space = isl_multi_aff_get_space(maff);
3544         space = isl_space_lift(isl_space_domain(space), n_div);
3545         space = isl_space_extend_domain_with_range(space,
3546                                                 isl_multi_aff_get_space(maff));
3547         if (!space)
3548                 return isl_multi_aff_free(maff);
3549         isl_space_free(maff->space);
3550         maff->space = space;
3551
3552         if (ls) {
3553                 *ls = isl_aff_get_domain_local_space(maff->p[0]);
3554                 if (!*ls)
3555                         return isl_multi_aff_free(maff);
3556         }
3557
3558         for (i = 0; i < maff->n; ++i) {
3559                 maff->p[i] = isl_aff_lift(maff->p[i]);
3560                 if (!maff->p[i])
3561                         goto error;
3562         }
3563
3564         return maff;
3565 error:
3566         if (ls)
3567                 isl_local_space_free(*ls);
3568         return isl_multi_aff_free(maff);
3569 }
3570
3571
3572 /* Extract an isl_pw_aff corresponding to output dimension "pos" of "pma".
3573  */
3574 __isl_give isl_pw_aff *isl_pw_multi_aff_get_pw_aff(
3575         __isl_keep isl_pw_multi_aff *pma, int pos)
3576 {
3577         int i;
3578         int n_out;
3579         isl_space *space;
3580         isl_pw_aff *pa;
3581
3582         if (!pma)
3583                 return NULL;
3584
3585         n_out = isl_pw_multi_aff_dim(pma, isl_dim_out);
3586         if (pos < 0 || pos >= n_out)
3587                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
3588                         "index out of bounds", return NULL);
3589
3590         space = isl_pw_multi_aff_get_space(pma);
3591         space = isl_space_drop_dims(space, isl_dim_out,
3592                                     pos + 1, n_out - pos - 1);
3593         space = isl_space_drop_dims(space, isl_dim_out, 0, pos);
3594
3595         pa = isl_pw_aff_alloc_size(space, pma->n);
3596         for (i = 0; i < pma->n; ++i) {
3597                 isl_aff *aff;
3598                 aff = isl_multi_aff_get_aff(pma->p[i].maff, pos);
3599                 pa = isl_pw_aff_add_piece(pa, isl_set_copy(pma->p[i].set), aff);
3600         }
3601
3602         return pa;
3603 }
3604
3605 /* Return an isl_pw_multi_aff with the given "set" as domain and
3606  * an unnamed zero-dimensional range.
3607  */
3608 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_from_domain(
3609         __isl_take isl_set *set)
3610 {
3611         isl_multi_aff *ma;
3612         isl_space *space;
3613
3614         space = isl_set_get_space(set);
3615         space = isl_space_from_domain(space);
3616         ma = isl_multi_aff_zero(space);
3617         return isl_pw_multi_aff_alloc(set, ma);
3618 }
3619
3620 /* Add an isl_pw_multi_aff with the given "set" as domain and
3621  * an unnamed zero-dimensional range to *user.
3622  */
3623 static int add_pw_multi_aff_from_domain(__isl_take isl_set *set, void *user)
3624 {
3625         isl_union_pw_multi_aff **upma = user;
3626         isl_pw_multi_aff *pma;
3627
3628         pma = isl_pw_multi_aff_from_domain(set);
3629         *upma = isl_union_pw_multi_aff_add_pw_multi_aff(*upma, pma);
3630
3631         return 0;
3632 }
3633
3634 /* Return an isl_union_pw_multi_aff with the given "uset" as domain and
3635  * an unnamed zero-dimensional range.
3636  */
3637 __isl_give isl_union_pw_multi_aff *isl_union_pw_multi_aff_from_domain(
3638         __isl_take isl_union_set *uset)
3639 {
3640         isl_space *space;
3641         isl_union_pw_multi_aff *upma;
3642
3643         if (!uset)
3644                 return NULL;
3645
3646         space = isl_union_set_get_space(uset);
3647         upma = isl_union_pw_multi_aff_empty(space);
3648
3649         if (isl_union_set_foreach_set(uset,
3650                                     &add_pw_multi_aff_from_domain, &upma) < 0)
3651                 goto error;
3652
3653         isl_union_set_free(uset);
3654         return upma;
3655 error:
3656         isl_union_set_free(uset);
3657         isl_union_pw_multi_aff_free(upma);
3658         return NULL;
3659 }
3660
3661 /* Convert "pma" to an isl_map and add it to *umap.
3662  */
3663 static int map_from_pw_multi_aff(__isl_take isl_pw_multi_aff *pma, void *user)
3664 {
3665         isl_union_map **umap = user;
3666         isl_map *map;
3667
3668         map = isl_map_from_pw_multi_aff(pma);
3669         *umap = isl_union_map_add_map(*umap, map);
3670
3671         return 0;
3672 }
3673
3674 /* Construct a union map mapping the domain of the union
3675  * piecewise multi-affine expression to its range, with each dimension
3676  * in the range equated to the corresponding affine expression on its cell.
3677  */
3678 __isl_give isl_union_map *isl_union_map_from_union_pw_multi_aff(
3679         __isl_take isl_union_pw_multi_aff *upma)
3680 {
3681         isl_space *space;
3682         isl_union_map *umap;
3683
3684         if (!upma)
3685                 return NULL;
3686
3687         space = isl_union_pw_multi_aff_get_space(upma);
3688         umap = isl_union_map_empty(space);
3689
3690         if (isl_union_pw_multi_aff_foreach_pw_multi_aff(upma,
3691                                         &map_from_pw_multi_aff, &umap) < 0)
3692                 goto error;
3693
3694         isl_union_pw_multi_aff_free(upma);
3695         return umap;
3696 error:
3697         isl_union_pw_multi_aff_free(upma);
3698         isl_union_map_free(umap);
3699         return NULL;
3700 }
3701
3702 /* Local data for bin_entry and the callback "fn".
3703  */
3704 struct isl_union_pw_multi_aff_bin_data {
3705         isl_union_pw_multi_aff *upma2;
3706         isl_union_pw_multi_aff *res;
3707         isl_pw_multi_aff *pma;
3708         int (*fn)(void **entry, void *user);
3709 };
3710
3711 /* Given an isl_pw_multi_aff from upma1, store it in data->pma
3712  * and call data->fn for each isl_pw_multi_aff in data->upma2.
3713  */
3714 static int bin_entry(void **entry, void *user)
3715 {
3716         struct isl_union_pw_multi_aff_bin_data *data = user;
3717         isl_pw_multi_aff *pma = *entry;
3718
3719         data->pma = pma;
3720         if (isl_hash_table_foreach(data->upma2->dim->ctx, &data->upma2->table,
3721                                    data->fn, data) < 0)
3722                 return -1;
3723
3724         return 0;
3725 }
3726
3727 /* Call "fn" on each pair of isl_pw_multi_affs in "upma1" and "upma2".
3728  * The isl_pw_multi_aff from upma1 is stored in data->pma (where data is
3729  * passed as user field) and the isl_pw_multi_aff from upma2 is available
3730  * as *entry.  The callback should adjust data->res if desired.
3731  */
3732 static __isl_give isl_union_pw_multi_aff *bin_op(
3733         __isl_take isl_union_pw_multi_aff *upma1,
3734         __isl_take isl_union_pw_multi_aff *upma2,
3735         int (*fn)(void **entry, void *user))
3736 {
3737         isl_space *space;
3738         struct isl_union_pw_multi_aff_bin_data data = { NULL, NULL, NULL, fn };
3739
3740         space = isl_union_pw_multi_aff_get_space(upma2);
3741         upma1 = isl_union_pw_multi_aff_align_params(upma1, space);
3742         space = isl_union_pw_multi_aff_get_space(upma1);
3743         upma2 = isl_union_pw_multi_aff_align_params(upma2, space);
3744
3745         if (!upma1 || !upma2)
3746                 goto error;
3747
3748         data.upma2 = upma2;
3749         data.res = isl_union_pw_multi_aff_alloc(isl_space_copy(upma1->dim),
3750                                        upma1->table.n);
3751         if (isl_hash_table_foreach(upma1->dim->ctx, &upma1->table,
3752                                    &bin_entry, &data) < 0)
3753                 goto error;
3754
3755         isl_union_pw_multi_aff_free(upma1);
3756         isl_union_pw_multi_aff_free(upma2);
3757         return data.res;
3758 error:
3759         isl_union_pw_multi_aff_free(upma1);
3760         isl_union_pw_multi_aff_free(upma2);
3761         isl_union_pw_multi_aff_free(data.res);
3762         return NULL;
3763 }
3764
3765 /* Given two isl_multi_affs A -> B and C -> D,
3766  * construct an isl_multi_aff (A * C) -> (B, D).
3767  */
3768 __isl_give isl_multi_aff *isl_multi_aff_flat_range_product(
3769         __isl_take isl_multi_aff *ma1, __isl_take isl_multi_aff *ma2)
3770 {
3771         int i, n1, n2;
3772         isl_aff *aff;
3773         isl_space *space;
3774         isl_multi_aff *res;
3775
3776         if (!ma1 || !ma2)
3777                 goto error;
3778
3779         space = isl_space_range_product(isl_multi_aff_get_space(ma1),
3780                                         isl_multi_aff_get_space(ma2));
3781         space = isl_space_flatten_range(space);
3782         res = isl_multi_aff_alloc(space);
3783
3784         n1 = isl_multi_aff_dim(ma1, isl_dim_out);
3785         n2 = isl_multi_aff_dim(ma2, isl_dim_out);
3786
3787         for (i = 0; i < n1; ++i) {
3788                 aff = isl_multi_aff_get_aff(ma1, i);
3789                 res = isl_multi_aff_set_aff(res, i, aff);
3790         }
3791
3792         for (i = 0; i < n2; ++i) {
3793                 aff = isl_multi_aff_get_aff(ma2, i);
3794                 res = isl_multi_aff_set_aff(res, n1 + i, aff);
3795         }
3796
3797         isl_multi_aff_free(ma1);
3798         isl_multi_aff_free(ma2);
3799         return res;
3800 error:
3801         isl_multi_aff_free(ma1);
3802         isl_multi_aff_free(ma2);
3803         return NULL;
3804 }
3805
3806 /* Given two aligned isl_pw_multi_affs A -> B and C -> D,
3807  * construct an isl_pw_multi_aff (A * C) -> (B, D).
3808  */
3809 static __isl_give isl_pw_multi_aff *pw_multi_aff_flat_range_product(
3810         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
3811 {
3812         isl_space *space;
3813
3814         space = isl_space_range_product(isl_pw_multi_aff_get_space(pma1),
3815                                         isl_pw_multi_aff_get_space(pma2));
3816         space = isl_space_flatten_range(space);
3817         return isl_pw_multi_aff_on_shared_domain_in(pma1, pma2, space,
3818                                             &isl_multi_aff_flat_range_product);
3819 }
3820
3821 /* Given two isl_pw_multi_affs A -> B and C -> D,
3822  * construct an isl_pw_multi_aff (A * C) -> (B, D).
3823  */
3824 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_flat_range_product(
3825         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
3826 {
3827         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
3828                                             &pw_multi_aff_flat_range_product);
3829 }
3830
3831 /* If data->pma and *entry have the same domain space, then compute
3832  * their flat range product and the result to data->res.
3833  */
3834 static int flat_range_product_entry(void **entry, void *user)
3835 {
3836         struct isl_union_pw_multi_aff_bin_data *data = user;
3837         isl_pw_multi_aff *pma2 = *entry;
3838
3839         if (!isl_space_tuple_match(data->pma->dim, isl_dim_in,
3840                                  pma2->dim, isl_dim_in))
3841                 return 0;
3842
3843         pma2 = isl_pw_multi_aff_flat_range_product(
3844                                         isl_pw_multi_aff_copy(data->pma),
3845                                         isl_pw_multi_aff_copy(pma2));
3846
3847         data->res = isl_union_pw_multi_aff_add_pw_multi_aff(data->res, pma2);
3848
3849         return 0;
3850 }
3851
3852 /* Given two isl_union_pw_multi_affs A -> B and C -> D,
3853  * construct an isl_union_pw_multi_aff (A * C) -> (B, D).
3854  */
3855 __isl_give isl_union_pw_multi_aff *isl_union_pw_multi_aff_flat_range_product(
3856         __isl_take isl_union_pw_multi_aff *upma1,
3857         __isl_take isl_union_pw_multi_aff *upma2)
3858 {
3859         return bin_op(upma1, upma2, &flat_range_product_entry);
3860 }
3861
3862 /* Replace the affine expressions at position "pos" in "pma" by "pa".
3863  * The parameters are assumed to have been aligned.
3864  *
3865  * The implementation essentially performs an isl_pw_*_on_shared_domain,
3866  * except that it works on two different isl_pw_* types.
3867  */
3868 static __isl_give isl_pw_multi_aff *pw_multi_aff_set_pw_aff(
3869         __isl_take isl_pw_multi_aff *pma, unsigned pos,
3870         __isl_take isl_pw_aff *pa)
3871 {
3872         int i, j, n;
3873         isl_pw_multi_aff *res = NULL;
3874
3875         if (!pma || !pa)
3876                 goto error;
3877
3878         if (!isl_space_tuple_match(pma->dim, isl_dim_in, pa->dim, isl_dim_in))
3879                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
3880                         "domains don't match", goto error);
3881         if (pos >= isl_pw_multi_aff_dim(pma, isl_dim_out))
3882                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
3883                         "index out of bounds", goto error);
3884
3885         n = pma->n * pa->n;
3886         res = isl_pw_multi_aff_alloc_size(isl_pw_multi_aff_get_space(pma), n);
3887
3888         for (i = 0; i < pma->n; ++i) {
3889                 for (j = 0; j < pa->n; ++j) {
3890                         isl_set *common;
3891                         isl_multi_aff *res_ij;
3892                         int empty;
3893
3894                         common = isl_set_intersect(isl_set_copy(pma->p[i].set),
3895                                                    isl_set_copy(pa->p[j].set));
3896                         empty = isl_set_plain_is_empty(common);
3897                         if (empty < 0 || empty) {
3898                                 isl_set_free(common);
3899                                 if (empty < 0)
3900                                         goto error;
3901                                 continue;
3902                         }
3903
3904                         res_ij = isl_multi_aff_set_aff(
3905                                         isl_multi_aff_copy(pma->p[i].maff), pos,
3906                                         isl_aff_copy(pa->p[j].aff));
3907                         res_ij = isl_multi_aff_gist(res_ij,
3908                                         isl_set_copy(common));
3909
3910                         res = isl_pw_multi_aff_add_piece(res, common, res_ij);
3911                 }
3912         }
3913
3914         isl_pw_multi_aff_free(pma);
3915         isl_pw_aff_free(pa);
3916         return res;
3917 error:
3918         isl_pw_multi_aff_free(pma);
3919         isl_pw_aff_free(pa);
3920         return isl_pw_multi_aff_free(res);
3921 }
3922
3923 /* Replace the affine expressions at position "pos" in "pma" by "pa".
3924  */
3925 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_set_pw_aff(
3926         __isl_take isl_pw_multi_aff *pma, unsigned pos,
3927         __isl_take isl_pw_aff *pa)
3928 {
3929         if (!pma || !pa)
3930                 goto error;
3931         if (isl_space_match(pma->dim, isl_dim_param, pa->dim, isl_dim_param))
3932                 return pw_multi_aff_set_pw_aff(pma, pos, pa);
3933         if (!isl_space_has_named_params(pma->dim) ||
3934             !isl_space_has_named_params(pa->dim))
3935                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
3936                         "unaligned unnamed parameters", goto error);
3937         pma = isl_pw_multi_aff_align_params(pma, isl_pw_aff_get_space(pa));
3938         pa = isl_pw_aff_align_params(pa, isl_pw_multi_aff_get_space(pma));
3939         return pw_multi_aff_set_pw_aff(pma, pos, pa);
3940 error:
3941         isl_pw_multi_aff_free(pma);
3942         isl_pw_aff_free(pa);
3943         return NULL;
3944 }