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