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