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  */
1392 __isl_give isl_basic_set *isl_aff_zero_basic_set(__isl_take isl_aff *aff)
1393 {
1394         isl_constraint *ineq;
1395         isl_basic_set *bset;
1396
1397         ineq = isl_equality_from_aff(aff);
1398
1399         bset = isl_basic_set_from_constraint(ineq);
1400         bset = isl_basic_set_simplify(bset);
1401         return bset;
1402 }
1403
1404 /* Return a basic set containing those elements in the shared space
1405  * of aff1 and aff2 where aff1 is greater than or equal to aff2.
1406  */
1407 __isl_give isl_basic_set *isl_aff_ge_basic_set(__isl_take isl_aff *aff1,
1408         __isl_take isl_aff *aff2)
1409 {
1410         aff1 = isl_aff_sub(aff1, aff2);
1411
1412         return isl_aff_nonneg_basic_set(aff1);
1413 }
1414
1415 /* Return a basic set containing those elements in the shared space
1416  * of aff1 and aff2 where aff1 is smaller than or equal to aff2.
1417  */
1418 __isl_give isl_basic_set *isl_aff_le_basic_set(__isl_take isl_aff *aff1,
1419         __isl_take isl_aff *aff2)
1420 {
1421         return isl_aff_ge_basic_set(aff2, aff1);
1422 }
1423
1424 __isl_give isl_aff *isl_aff_add_on_domain(__isl_keep isl_set *dom,
1425         __isl_take isl_aff *aff1, __isl_take isl_aff *aff2)
1426 {
1427         aff1 = isl_aff_add(aff1, aff2);
1428         aff1 = isl_aff_gist(aff1, isl_set_copy(dom));
1429         return aff1;
1430 }
1431
1432 int isl_aff_is_empty(__isl_keep isl_aff *aff)
1433 {
1434         if (!aff)
1435                 return -1;
1436
1437         return 0;
1438 }
1439
1440 /* Check whether the given affine expression has non-zero coefficient
1441  * for any dimension in the given range or if any of these dimensions
1442  * appear with non-zero coefficients in any of the integer divisions
1443  * involved in the affine expression.
1444  */
1445 int isl_aff_involves_dims(__isl_keep isl_aff *aff,
1446         enum isl_dim_type type, unsigned first, unsigned n)
1447 {
1448         int i;
1449         isl_ctx *ctx;
1450         int *active = NULL;
1451         int involves = 0;
1452
1453         if (!aff)
1454                 return -1;
1455         if (n == 0)
1456                 return 0;
1457
1458         ctx = isl_aff_get_ctx(aff);
1459         if (first + n > isl_aff_dim(aff, type))
1460                 isl_die(ctx, isl_error_invalid,
1461                         "range out of bounds", return -1);
1462
1463         active = isl_local_space_get_active(aff->ls, aff->v->el + 2);
1464         if (!active)
1465                 goto error;
1466
1467         first += isl_local_space_offset(aff->ls, type) - 1;
1468         for (i = 0; i < n; ++i)
1469                 if (active[first + i]) {
1470                         involves = 1;
1471                         break;
1472                 }
1473
1474         free(active);
1475
1476         return involves;
1477 error:
1478         free(active);
1479         return -1;
1480 }
1481
1482 __isl_give isl_aff *isl_aff_drop_dims(__isl_take isl_aff *aff,
1483         enum isl_dim_type type, unsigned first, unsigned n)
1484 {
1485         isl_ctx *ctx;
1486
1487         if (!aff)
1488                 return NULL;
1489         if (type == isl_dim_out)
1490                 isl_die(aff->v->ctx, isl_error_invalid,
1491                         "cannot drop output/set dimension",
1492                         return isl_aff_free(aff));
1493         if (type == isl_dim_in)
1494                 type = isl_dim_set;
1495         if (n == 0 && !isl_local_space_is_named_or_nested(aff->ls, type))
1496                 return aff;
1497
1498         ctx = isl_aff_get_ctx(aff);
1499         if (first + n > isl_local_space_dim(aff->ls, type))
1500                 isl_die(ctx, isl_error_invalid, "range out of bounds",
1501                         return isl_aff_free(aff));
1502
1503         aff = isl_aff_cow(aff);
1504         if (!aff)
1505                 return NULL;
1506
1507         aff->ls = isl_local_space_drop_dims(aff->ls, type, first, n);
1508         if (!aff->ls)
1509                 return isl_aff_free(aff);
1510
1511         first += 1 + isl_local_space_offset(aff->ls, type);
1512         aff->v = isl_vec_drop_els(aff->v, first, n);
1513         if (!aff->v)
1514                 return isl_aff_free(aff);
1515
1516         return aff;
1517 }
1518
1519 /* Project the domain of the affine expression onto its parameter space.
1520  * The affine expression may not involve any of the domain dimensions.
1521  */
1522 __isl_give isl_aff *isl_aff_project_domain_on_params(__isl_take isl_aff *aff)
1523 {
1524         isl_space *space;
1525         unsigned n;
1526         int involves;
1527
1528         n = isl_aff_dim(aff, isl_dim_in);
1529         involves = isl_aff_involves_dims(aff, isl_dim_in, 0, n);
1530         if (involves < 0)
1531                 return isl_aff_free(aff);
1532         if (involves)
1533                 isl_die(isl_aff_get_ctx(aff), isl_error_invalid,
1534                     "affine expression involves some of the domain dimensions",
1535                     return isl_aff_free(aff));
1536         aff = isl_aff_drop_dims(aff, isl_dim_in, 0, n);
1537         space = isl_aff_get_domain_space(aff);
1538         space = isl_space_params(space);
1539         aff = isl_aff_reset_domain_space(aff, space);
1540         return aff;
1541 }
1542
1543 __isl_give isl_aff *isl_aff_insert_dims(__isl_take isl_aff *aff,
1544         enum isl_dim_type type, unsigned first, unsigned n)
1545 {
1546         isl_ctx *ctx;
1547
1548         if (!aff)
1549                 return NULL;
1550         if (type == isl_dim_out)
1551                 isl_die(aff->v->ctx, isl_error_invalid,
1552                         "cannot insert output/set dimensions",
1553                         return isl_aff_free(aff));
1554         if (type == isl_dim_in)
1555                 type = isl_dim_set;
1556         if (n == 0 && !isl_local_space_is_named_or_nested(aff->ls, type))
1557                 return aff;
1558
1559         ctx = isl_aff_get_ctx(aff);
1560         if (first > isl_local_space_dim(aff->ls, type))
1561                 isl_die(ctx, isl_error_invalid, "position out of bounds",
1562                         return isl_aff_free(aff));
1563
1564         aff = isl_aff_cow(aff);
1565         if (!aff)
1566                 return NULL;
1567
1568         aff->ls = isl_local_space_insert_dims(aff->ls, type, first, n);
1569         if (!aff->ls)
1570                 return isl_aff_free(aff);
1571
1572         first += 1 + isl_local_space_offset(aff->ls, type);
1573         aff->v = isl_vec_insert_zero_els(aff->v, first, n);
1574         if (!aff->v)
1575                 return isl_aff_free(aff);
1576
1577         return aff;
1578 }
1579
1580 __isl_give isl_aff *isl_aff_add_dims(__isl_take isl_aff *aff,
1581         enum isl_dim_type type, unsigned n)
1582 {
1583         unsigned pos;
1584
1585         pos = isl_aff_dim(aff, type);
1586
1587         return isl_aff_insert_dims(aff, type, pos, n);
1588 }
1589
1590 __isl_give isl_pw_aff *isl_pw_aff_add_dims(__isl_take isl_pw_aff *pwaff,
1591         enum isl_dim_type type, unsigned n)
1592 {
1593         unsigned pos;
1594
1595         pos = isl_pw_aff_dim(pwaff, type);
1596
1597         return isl_pw_aff_insert_dims(pwaff, type, pos, n);
1598 }
1599
1600 __isl_give isl_pw_aff *isl_pw_aff_from_aff(__isl_take isl_aff *aff)
1601 {
1602         isl_set *dom = isl_set_universe(isl_aff_get_domain_space(aff));
1603         return isl_pw_aff_alloc(dom, aff);
1604 }
1605
1606 #undef PW
1607 #define PW isl_pw_aff
1608 #undef EL
1609 #define EL isl_aff
1610 #undef EL_IS_ZERO
1611 #define EL_IS_ZERO is_empty
1612 #undef ZERO
1613 #define ZERO empty
1614 #undef IS_ZERO
1615 #define IS_ZERO is_empty
1616 #undef FIELD
1617 #define FIELD aff
1618 #undef DEFAULT_IS_ZERO
1619 #define DEFAULT_IS_ZERO 0
1620
1621 #define NO_EVAL
1622 #define NO_OPT
1623 #define NO_MOVE_DIMS
1624 #define NO_LIFT
1625 #define NO_MORPH
1626
1627 #include <isl_pw_templ.c>
1628
1629 static __isl_give isl_set *align_params_pw_pw_set_and(
1630         __isl_take isl_pw_aff *pwaff1, __isl_take isl_pw_aff *pwaff2,
1631         __isl_give isl_set *(*fn)(__isl_take isl_pw_aff *pwaff1,
1632                                     __isl_take isl_pw_aff *pwaff2))
1633 {
1634         if (!pwaff1 || !pwaff2)
1635                 goto error;
1636         if (isl_space_match(pwaff1->dim, isl_dim_param,
1637                           pwaff2->dim, isl_dim_param))
1638                 return fn(pwaff1, pwaff2);
1639         if (!isl_space_has_named_params(pwaff1->dim) ||
1640             !isl_space_has_named_params(pwaff2->dim))
1641                 isl_die(isl_pw_aff_get_ctx(pwaff1), isl_error_invalid,
1642                         "unaligned unnamed parameters", goto error);
1643         pwaff1 = isl_pw_aff_align_params(pwaff1, isl_pw_aff_get_space(pwaff2));
1644         pwaff2 = isl_pw_aff_align_params(pwaff2, isl_pw_aff_get_space(pwaff1));
1645         return fn(pwaff1, pwaff2);
1646 error:
1647         isl_pw_aff_free(pwaff1);
1648         isl_pw_aff_free(pwaff2);
1649         return NULL;
1650 }
1651
1652 /* Compute a piecewise quasi-affine expression with a domain that
1653  * is the union of those of pwaff1 and pwaff2 and such that on each
1654  * cell, the quasi-affine expression is the better (according to cmp)
1655  * of those of pwaff1 and pwaff2.  If only one of pwaff1 or pwaff2
1656  * is defined on a given cell, then the associated expression
1657  * is the defined one.
1658  */
1659 static __isl_give isl_pw_aff *pw_aff_union_opt(__isl_take isl_pw_aff *pwaff1,
1660         __isl_take isl_pw_aff *pwaff2,
1661         __isl_give isl_basic_set *(*cmp)(__isl_take isl_aff *aff1,
1662                                         __isl_take isl_aff *aff2))
1663 {
1664         int i, j, n;
1665         isl_pw_aff *res;
1666         isl_ctx *ctx;
1667         isl_set *set;
1668
1669         if (!pwaff1 || !pwaff2)
1670                 goto error;
1671
1672         ctx = isl_space_get_ctx(pwaff1->dim);
1673         if (!isl_space_is_equal(pwaff1->dim, pwaff2->dim))
1674                 isl_die(ctx, isl_error_invalid,
1675                         "arguments should live in same space", goto error);
1676
1677         if (isl_pw_aff_is_empty(pwaff1)) {
1678                 isl_pw_aff_free(pwaff1);
1679                 return pwaff2;
1680         }
1681
1682         if (isl_pw_aff_is_empty(pwaff2)) {
1683                 isl_pw_aff_free(pwaff2);
1684                 return pwaff1;
1685         }
1686
1687         n = 2 * (pwaff1->n + 1) * (pwaff2->n + 1);
1688         res = isl_pw_aff_alloc_size(isl_space_copy(pwaff1->dim), n);
1689
1690         for (i = 0; i < pwaff1->n; ++i) {
1691                 set = isl_set_copy(pwaff1->p[i].set);
1692                 for (j = 0; j < pwaff2->n; ++j) {
1693                         struct isl_set *common;
1694                         isl_set *better;
1695
1696                         common = isl_set_intersect(
1697                                         isl_set_copy(pwaff1->p[i].set),
1698                                         isl_set_copy(pwaff2->p[j].set));
1699                         better = isl_set_from_basic_set(cmp(
1700                                         isl_aff_copy(pwaff2->p[j].aff),
1701                                         isl_aff_copy(pwaff1->p[i].aff)));
1702                         better = isl_set_intersect(common, better);
1703                         if (isl_set_plain_is_empty(better)) {
1704                                 isl_set_free(better);
1705                                 continue;
1706                         }
1707                         set = isl_set_subtract(set, isl_set_copy(better));
1708
1709                         res = isl_pw_aff_add_piece(res, better,
1710                                                 isl_aff_copy(pwaff2->p[j].aff));
1711                 }
1712                 res = isl_pw_aff_add_piece(res, set,
1713                                                 isl_aff_copy(pwaff1->p[i].aff));
1714         }
1715
1716         for (j = 0; j < pwaff2->n; ++j) {
1717                 set = isl_set_copy(pwaff2->p[j].set);
1718                 for (i = 0; i < pwaff1->n; ++i)
1719                         set = isl_set_subtract(set,
1720                                         isl_set_copy(pwaff1->p[i].set));
1721                 res = isl_pw_aff_add_piece(res, set,
1722                                                 isl_aff_copy(pwaff2->p[j].aff));
1723         }
1724
1725         isl_pw_aff_free(pwaff1);
1726         isl_pw_aff_free(pwaff2);
1727
1728         return res;
1729 error:
1730         isl_pw_aff_free(pwaff1);
1731         isl_pw_aff_free(pwaff2);
1732         return NULL;
1733 }
1734
1735 /* Compute a piecewise quasi-affine expression with a domain that
1736  * is the union of those of pwaff1 and pwaff2 and such that on each
1737  * cell, the quasi-affine expression is the maximum of those of pwaff1
1738  * and pwaff2.  If only one of pwaff1 or pwaff2 is defined on a given
1739  * cell, then the associated expression is the defined one.
1740  */
1741 static __isl_give isl_pw_aff *pw_aff_union_max(__isl_take isl_pw_aff *pwaff1,
1742         __isl_take isl_pw_aff *pwaff2)
1743 {
1744         return pw_aff_union_opt(pwaff1, pwaff2, &isl_aff_ge_basic_set);
1745 }
1746
1747 __isl_give isl_pw_aff *isl_pw_aff_union_max(__isl_take isl_pw_aff *pwaff1,
1748         __isl_take isl_pw_aff *pwaff2)
1749 {
1750         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2,
1751                                                         &pw_aff_union_max);
1752 }
1753
1754 /* Compute a piecewise quasi-affine expression with a domain that
1755  * is the union of those of pwaff1 and pwaff2 and such that on each
1756  * cell, the quasi-affine expression is the minimum of those of pwaff1
1757  * and pwaff2.  If only one of pwaff1 or pwaff2 is defined on a given
1758  * cell, then the associated expression is the defined one.
1759  */
1760 static __isl_give isl_pw_aff *pw_aff_union_min(__isl_take isl_pw_aff *pwaff1,
1761         __isl_take isl_pw_aff *pwaff2)
1762 {
1763         return pw_aff_union_opt(pwaff1, pwaff2, &isl_aff_le_basic_set);
1764 }
1765
1766 __isl_give isl_pw_aff *isl_pw_aff_union_min(__isl_take isl_pw_aff *pwaff1,
1767         __isl_take isl_pw_aff *pwaff2)
1768 {
1769         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2,
1770                                                         &pw_aff_union_min);
1771 }
1772
1773 __isl_give isl_pw_aff *isl_pw_aff_union_opt(__isl_take isl_pw_aff *pwaff1,
1774         __isl_take isl_pw_aff *pwaff2, int max)
1775 {
1776         if (max)
1777                 return isl_pw_aff_union_max(pwaff1, pwaff2);
1778         else
1779                 return isl_pw_aff_union_min(pwaff1, pwaff2);
1780 }
1781
1782 /* Construct a map with as domain the domain of pwaff and
1783  * one-dimensional range corresponding to the affine expressions.
1784  */
1785 static __isl_give isl_map *map_from_pw_aff(__isl_take isl_pw_aff *pwaff)
1786 {
1787         int i;
1788         isl_space *dim;
1789         isl_map *map;
1790
1791         if (!pwaff)
1792                 return NULL;
1793
1794         dim = isl_pw_aff_get_space(pwaff);
1795         map = isl_map_empty(dim);
1796
1797         for (i = 0; i < pwaff->n; ++i) {
1798                 isl_basic_map *bmap;
1799                 isl_map *map_i;
1800
1801                 bmap = isl_basic_map_from_aff(isl_aff_copy(pwaff->p[i].aff));
1802                 map_i = isl_map_from_basic_map(bmap);
1803                 map_i = isl_map_intersect_domain(map_i,
1804                                                 isl_set_copy(pwaff->p[i].set));
1805                 map = isl_map_union_disjoint(map, map_i);
1806         }
1807
1808         isl_pw_aff_free(pwaff);
1809
1810         return map;
1811 }
1812
1813 /* Construct a map with as domain the domain of pwaff and
1814  * one-dimensional range corresponding to the affine expressions.
1815  */
1816 __isl_give isl_map *isl_map_from_pw_aff(__isl_take isl_pw_aff *pwaff)
1817 {
1818         if (!pwaff)
1819                 return NULL;
1820         if (isl_space_is_set(pwaff->dim))
1821                 isl_die(isl_pw_aff_get_ctx(pwaff), isl_error_invalid,
1822                         "space of input is not a map",
1823                         return isl_pw_aff_free(pwaff));
1824         return map_from_pw_aff(pwaff);
1825 }
1826
1827 /* Construct a one-dimensional set with as parameter domain
1828  * the domain of pwaff and the single set dimension
1829  * corresponding to the affine expressions.
1830  */
1831 __isl_give isl_set *isl_set_from_pw_aff(__isl_take isl_pw_aff *pwaff)
1832 {
1833         if (!pwaff)
1834                 return NULL;
1835         if (!isl_space_is_set(pwaff->dim))
1836                 isl_die(isl_pw_aff_get_ctx(pwaff), isl_error_invalid,
1837                         "space of input is not a set",
1838                         return isl_pw_aff_free(pwaff));
1839         return map_from_pw_aff(pwaff);
1840 }
1841
1842 /* Return a set containing those elements in the domain
1843  * of pwaff where it is non-negative.
1844  */
1845 __isl_give isl_set *isl_pw_aff_nonneg_set(__isl_take isl_pw_aff *pwaff)
1846 {
1847         int i;
1848         isl_set *set;
1849
1850         if (!pwaff)
1851                 return NULL;
1852
1853         set = isl_set_empty(isl_pw_aff_get_domain_space(pwaff));
1854
1855         for (i = 0; i < pwaff->n; ++i) {
1856                 isl_basic_set *bset;
1857                 isl_set *set_i;
1858
1859                 bset = isl_aff_nonneg_basic_set(isl_aff_copy(pwaff->p[i].aff));
1860                 set_i = isl_set_from_basic_set(bset);
1861                 set_i = isl_set_intersect(set_i, isl_set_copy(pwaff->p[i].set));
1862                 set = isl_set_union_disjoint(set, set_i);
1863         }
1864
1865         isl_pw_aff_free(pwaff);
1866
1867         return set;
1868 }
1869
1870 /* Return a set containing those elements in the domain
1871  * of pwaff where it is zero (if complement is 0) or not zero
1872  * (if complement is 1).
1873  */
1874 static __isl_give isl_set *pw_aff_zero_set(__isl_take isl_pw_aff *pwaff,
1875         int complement)
1876 {
1877         int i;
1878         isl_set *set;
1879
1880         if (!pwaff)
1881                 return NULL;
1882
1883         set = isl_set_empty(isl_pw_aff_get_domain_space(pwaff));
1884
1885         for (i = 0; i < pwaff->n; ++i) {
1886                 isl_basic_set *bset;
1887                 isl_set *set_i, *zero;
1888
1889                 bset = isl_aff_zero_basic_set(isl_aff_copy(pwaff->p[i].aff));
1890                 zero = isl_set_from_basic_set(bset);
1891                 set_i = isl_set_copy(pwaff->p[i].set);
1892                 if (complement)
1893                         set_i = isl_set_subtract(set_i, zero);
1894                 else
1895                         set_i = isl_set_intersect(set_i, zero);
1896                 set = isl_set_union_disjoint(set, set_i);
1897         }
1898
1899         isl_pw_aff_free(pwaff);
1900
1901         return set;
1902 }
1903
1904 /* Return a set containing those elements in the domain
1905  * of pwaff where it is zero.
1906  */
1907 __isl_give isl_set *isl_pw_aff_zero_set(__isl_take isl_pw_aff *pwaff)
1908 {
1909         return pw_aff_zero_set(pwaff, 0);
1910 }
1911
1912 /* Return a set containing those elements in the domain
1913  * of pwaff where it is not zero.
1914  */
1915 __isl_give isl_set *isl_pw_aff_non_zero_set(__isl_take isl_pw_aff *pwaff)
1916 {
1917         return pw_aff_zero_set(pwaff, 1);
1918 }
1919
1920 /* Return a set containing those elements in the shared domain
1921  * of pwaff1 and pwaff2 where pwaff1 is greater than (or equal) to pwaff2.
1922  *
1923  * We compute the difference on the shared domain and then construct
1924  * the set of values where this difference is non-negative.
1925  * If strict is set, we first subtract 1 from the difference.
1926  * If equal is set, we only return the elements where pwaff1 and pwaff2
1927  * are equal.
1928  */
1929 static __isl_give isl_set *pw_aff_gte_set(__isl_take isl_pw_aff *pwaff1,
1930         __isl_take isl_pw_aff *pwaff2, int strict, int equal)
1931 {
1932         isl_set *set1, *set2;
1933
1934         set1 = isl_pw_aff_domain(isl_pw_aff_copy(pwaff1));
1935         set2 = isl_pw_aff_domain(isl_pw_aff_copy(pwaff2));
1936         set1 = isl_set_intersect(set1, set2);
1937         pwaff1 = isl_pw_aff_intersect_domain(pwaff1, isl_set_copy(set1));
1938         pwaff2 = isl_pw_aff_intersect_domain(pwaff2, isl_set_copy(set1));
1939         pwaff1 = isl_pw_aff_add(pwaff1, isl_pw_aff_neg(pwaff2));
1940
1941         if (strict) {
1942                 isl_space *dim = isl_set_get_space(set1);
1943                 isl_aff *aff;
1944                 aff = isl_aff_zero_on_domain(isl_local_space_from_space(dim));
1945                 aff = isl_aff_add_constant_si(aff, -1);
1946                 pwaff1 = isl_pw_aff_add(pwaff1, isl_pw_aff_alloc(set1, aff));
1947         } else
1948                 isl_set_free(set1);
1949
1950         if (equal)
1951                 return isl_pw_aff_zero_set(pwaff1);
1952         return isl_pw_aff_nonneg_set(pwaff1);
1953 }
1954
1955 /* Return a set containing those elements in the shared domain
1956  * of pwaff1 and pwaff2 where pwaff1 is equal to pwaff2.
1957  */
1958 static __isl_give isl_set *pw_aff_eq_set(__isl_take isl_pw_aff *pwaff1,
1959         __isl_take isl_pw_aff *pwaff2)
1960 {
1961         return pw_aff_gte_set(pwaff1, pwaff2, 0, 1);
1962 }
1963
1964 __isl_give isl_set *isl_pw_aff_eq_set(__isl_take isl_pw_aff *pwaff1,
1965         __isl_take isl_pw_aff *pwaff2)
1966 {
1967         return align_params_pw_pw_set_and(pwaff1, pwaff2, &pw_aff_eq_set);
1968 }
1969
1970 /* Return a set containing those elements in the shared domain
1971  * of pwaff1 and pwaff2 where pwaff1 is greater than or equal to pwaff2.
1972  */
1973 static __isl_give isl_set *pw_aff_ge_set(__isl_take isl_pw_aff *pwaff1,
1974         __isl_take isl_pw_aff *pwaff2)
1975 {
1976         return pw_aff_gte_set(pwaff1, pwaff2, 0, 0);
1977 }
1978
1979 __isl_give isl_set *isl_pw_aff_ge_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_ge_set);
1983 }
1984
1985 /* Return a set containing those elements in the shared domain
1986  * of pwaff1 and pwaff2 where pwaff1 is strictly greater than pwaff2.
1987  */
1988 static __isl_give isl_set *pw_aff_gt_set(__isl_take isl_pw_aff *pwaff1,
1989         __isl_take isl_pw_aff *pwaff2)
1990 {
1991         return pw_aff_gte_set(pwaff1, pwaff2, 1, 0);
1992 }
1993
1994 __isl_give isl_set *isl_pw_aff_gt_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_gt_set);
1998 }
1999
2000 __isl_give isl_set *isl_pw_aff_le_set(__isl_take isl_pw_aff *pwaff1,
2001         __isl_take isl_pw_aff *pwaff2)
2002 {
2003         return isl_pw_aff_ge_set(pwaff2, pwaff1);
2004 }
2005
2006 __isl_give isl_set *isl_pw_aff_lt_set(__isl_take isl_pw_aff *pwaff1,
2007         __isl_take isl_pw_aff *pwaff2)
2008 {
2009         return isl_pw_aff_gt_set(pwaff2, pwaff1);
2010 }
2011
2012 /* Return a set containing those elements in the shared domain
2013  * of the elements of list1 and list2 where each element in list1
2014  * has the relation specified by "fn" with each element in list2.
2015  */
2016 static __isl_give isl_set *pw_aff_list_set(__isl_take isl_pw_aff_list *list1,
2017         __isl_take isl_pw_aff_list *list2,
2018         __isl_give isl_set *(*fn)(__isl_take isl_pw_aff *pwaff1,
2019                                     __isl_take isl_pw_aff *pwaff2))
2020 {
2021         int i, j;
2022         isl_ctx *ctx;
2023         isl_set *set;
2024
2025         if (!list1 || !list2)
2026                 goto error;
2027
2028         ctx = isl_pw_aff_list_get_ctx(list1);
2029         if (list1->n < 1 || list2->n < 1)
2030                 isl_die(ctx, isl_error_invalid,
2031                         "list should contain at least one element", goto error);
2032
2033         set = isl_set_universe(isl_pw_aff_get_domain_space(list1->p[0]));
2034         for (i = 0; i < list1->n; ++i)
2035                 for (j = 0; j < list2->n; ++j) {
2036                         isl_set *set_ij;
2037
2038                         set_ij = fn(isl_pw_aff_copy(list1->p[i]),
2039                                     isl_pw_aff_copy(list2->p[j]));
2040                         set = isl_set_intersect(set, set_ij);
2041                 }
2042
2043         isl_pw_aff_list_free(list1);
2044         isl_pw_aff_list_free(list2);
2045         return set;
2046 error:
2047         isl_pw_aff_list_free(list1);
2048         isl_pw_aff_list_free(list2);
2049         return NULL;
2050 }
2051
2052 /* Return a set containing those elements in the shared domain
2053  * of the elements of list1 and list2 where each element in list1
2054  * is equal to each element in list2.
2055  */
2056 __isl_give isl_set *isl_pw_aff_list_eq_set(__isl_take isl_pw_aff_list *list1,
2057         __isl_take isl_pw_aff_list *list2)
2058 {
2059         return pw_aff_list_set(list1, list2, &isl_pw_aff_eq_set);
2060 }
2061
2062 __isl_give isl_set *isl_pw_aff_list_ne_set(__isl_take isl_pw_aff_list *list1,
2063         __isl_take isl_pw_aff_list *list2)
2064 {
2065         return pw_aff_list_set(list1, list2, &isl_pw_aff_ne_set);
2066 }
2067
2068 /* Return a set containing those elements in the shared domain
2069  * of the elements of list1 and list2 where each element in list1
2070  * is less than or equal to each element in list2.
2071  */
2072 __isl_give isl_set *isl_pw_aff_list_le_set(__isl_take isl_pw_aff_list *list1,
2073         __isl_take isl_pw_aff_list *list2)
2074 {
2075         return pw_aff_list_set(list1, list2, &isl_pw_aff_le_set);
2076 }
2077
2078 __isl_give isl_set *isl_pw_aff_list_lt_set(__isl_take isl_pw_aff_list *list1,
2079         __isl_take isl_pw_aff_list *list2)
2080 {
2081         return pw_aff_list_set(list1, list2, &isl_pw_aff_lt_set);
2082 }
2083
2084 __isl_give isl_set *isl_pw_aff_list_ge_set(__isl_take isl_pw_aff_list *list1,
2085         __isl_take isl_pw_aff_list *list2)
2086 {
2087         return pw_aff_list_set(list1, list2, &isl_pw_aff_ge_set);
2088 }
2089
2090 __isl_give isl_set *isl_pw_aff_list_gt_set(__isl_take isl_pw_aff_list *list1,
2091         __isl_take isl_pw_aff_list *list2)
2092 {
2093         return pw_aff_list_set(list1, list2, &isl_pw_aff_gt_set);
2094 }
2095
2096
2097 /* Return a set containing those elements in the shared domain
2098  * of pwaff1 and pwaff2 where pwaff1 is not equal to pwaff2.
2099  */
2100 static __isl_give isl_set *pw_aff_ne_set(__isl_take isl_pw_aff *pwaff1,
2101         __isl_take isl_pw_aff *pwaff2)
2102 {
2103         isl_set *set_lt, *set_gt;
2104
2105         set_lt = isl_pw_aff_lt_set(isl_pw_aff_copy(pwaff1),
2106                                    isl_pw_aff_copy(pwaff2));
2107         set_gt = isl_pw_aff_gt_set(pwaff1, pwaff2);
2108         return isl_set_union_disjoint(set_lt, set_gt);
2109 }
2110
2111 __isl_give isl_set *isl_pw_aff_ne_set(__isl_take isl_pw_aff *pwaff1,
2112         __isl_take isl_pw_aff *pwaff2)
2113 {
2114         return align_params_pw_pw_set_and(pwaff1, pwaff2, &pw_aff_ne_set);
2115 }
2116
2117 __isl_give isl_pw_aff *isl_pw_aff_scale_down(__isl_take isl_pw_aff *pwaff,
2118         isl_int v)
2119 {
2120         int i;
2121
2122         if (isl_int_is_one(v))
2123                 return pwaff;
2124         if (!isl_int_is_pos(v))
2125                 isl_die(isl_pw_aff_get_ctx(pwaff), isl_error_invalid,
2126                         "factor needs to be positive",
2127                         return isl_pw_aff_free(pwaff));
2128         pwaff = isl_pw_aff_cow(pwaff);
2129         if (!pwaff)
2130                 return NULL;
2131         if (pwaff->n == 0)
2132                 return pwaff;
2133
2134         for (i = 0; i < pwaff->n; ++i) {
2135                 pwaff->p[i].aff = isl_aff_scale_down(pwaff->p[i].aff, v);
2136                 if (!pwaff->p[i].aff)
2137                         return isl_pw_aff_free(pwaff);
2138         }
2139
2140         return pwaff;
2141 }
2142
2143 __isl_give isl_pw_aff *isl_pw_aff_floor(__isl_take isl_pw_aff *pwaff)
2144 {
2145         int i;
2146
2147         pwaff = isl_pw_aff_cow(pwaff);
2148         if (!pwaff)
2149                 return NULL;
2150         if (pwaff->n == 0)
2151                 return pwaff;
2152
2153         for (i = 0; i < pwaff->n; ++i) {
2154                 pwaff->p[i].aff = isl_aff_floor(pwaff->p[i].aff);
2155                 if (!pwaff->p[i].aff)
2156                         return isl_pw_aff_free(pwaff);
2157         }
2158
2159         return pwaff;
2160 }
2161
2162 __isl_give isl_pw_aff *isl_pw_aff_ceil(__isl_take isl_pw_aff *pwaff)
2163 {
2164         int i;
2165
2166         pwaff = isl_pw_aff_cow(pwaff);
2167         if (!pwaff)
2168                 return NULL;
2169         if (pwaff->n == 0)
2170                 return pwaff;
2171
2172         for (i = 0; i < pwaff->n; ++i) {
2173                 pwaff->p[i].aff = isl_aff_ceil(pwaff->p[i].aff);
2174                 if (!pwaff->p[i].aff)
2175                         return isl_pw_aff_free(pwaff);
2176         }
2177
2178         return pwaff;
2179 }
2180
2181 /* Assuming that "cond1" and "cond2" are disjoint,
2182  * return an affine expression that is equal to pwaff1 on cond1
2183  * and to pwaff2 on cond2.
2184  */
2185 static __isl_give isl_pw_aff *isl_pw_aff_select(
2186         __isl_take isl_set *cond1, __isl_take isl_pw_aff *pwaff1,
2187         __isl_take isl_set *cond2, __isl_take isl_pw_aff *pwaff2)
2188 {
2189         pwaff1 = isl_pw_aff_intersect_domain(pwaff1, cond1);
2190         pwaff2 = isl_pw_aff_intersect_domain(pwaff2, cond2);
2191
2192         return isl_pw_aff_add_disjoint(pwaff1, pwaff2);
2193 }
2194
2195 /* Return an affine expression that is equal to pwaff_true for elements
2196  * where "cond" is non-zero and to pwaff_false for elements where "cond"
2197  * is zero.
2198  * That is, return cond ? pwaff_true : pwaff_false;
2199  */
2200 __isl_give isl_pw_aff *isl_pw_aff_cond(__isl_take isl_pw_aff *cond,
2201         __isl_take isl_pw_aff *pwaff_true, __isl_take isl_pw_aff *pwaff_false)
2202 {
2203         isl_set *cond_true, *cond_false;
2204
2205         cond_true = isl_pw_aff_non_zero_set(isl_pw_aff_copy(cond));
2206         cond_false = isl_pw_aff_zero_set(cond);
2207         return isl_pw_aff_select(cond_true, pwaff_true,
2208                                  cond_false, pwaff_false);
2209 }
2210
2211 int isl_aff_is_cst(__isl_keep isl_aff *aff)
2212 {
2213         if (!aff)
2214                 return -1;
2215
2216         return isl_seq_first_non_zero(aff->v->el + 2, aff->v->size - 2) == -1;
2217 }
2218
2219 /* Check whether pwaff is a piecewise constant.
2220  */
2221 int isl_pw_aff_is_cst(__isl_keep isl_pw_aff *pwaff)
2222 {
2223         int i;
2224
2225         if (!pwaff)
2226                 return -1;
2227
2228         for (i = 0; i < pwaff->n; ++i) {
2229                 int is_cst = isl_aff_is_cst(pwaff->p[i].aff);
2230                 if (is_cst < 0 || !is_cst)
2231                         return is_cst;
2232         }
2233
2234         return 1;
2235 }
2236
2237 __isl_give isl_aff *isl_aff_mul(__isl_take isl_aff *aff1,
2238         __isl_take isl_aff *aff2)
2239 {
2240         if (!isl_aff_is_cst(aff2) && isl_aff_is_cst(aff1))
2241                 return isl_aff_mul(aff2, aff1);
2242
2243         if (!isl_aff_is_cst(aff2))
2244                 isl_die(isl_aff_get_ctx(aff1), isl_error_invalid,
2245                         "at least one affine expression should be constant",
2246                         goto error);
2247
2248         aff1 = isl_aff_cow(aff1);
2249         if (!aff1 || !aff2)
2250                 goto error;
2251
2252         aff1 = isl_aff_scale(aff1, aff2->v->el[1]);
2253         aff1 = isl_aff_scale_down(aff1, aff2->v->el[0]);
2254
2255         isl_aff_free(aff2);
2256         return aff1;
2257 error:
2258         isl_aff_free(aff1);
2259         isl_aff_free(aff2);
2260         return NULL;
2261 }
2262
2263 /* Divide "aff1" by "aff2", assuming "aff2" is a piecewise constant.
2264  */
2265 __isl_give isl_aff *isl_aff_div(__isl_take isl_aff *aff1,
2266         __isl_take isl_aff *aff2)
2267 {
2268         int is_cst;
2269         int neg;
2270
2271         is_cst = isl_aff_is_cst(aff2);
2272         if (is_cst < 0)
2273                 goto error;
2274         if (!is_cst)
2275                 isl_die(isl_aff_get_ctx(aff2), isl_error_invalid,
2276                         "second argument should be a constant", goto error);
2277
2278         if (!aff2)
2279                 goto error;
2280
2281         neg = isl_int_is_neg(aff2->v->el[1]);
2282         if (neg) {
2283                 isl_int_neg(aff2->v->el[0], aff2->v->el[0]);
2284                 isl_int_neg(aff2->v->el[1], aff2->v->el[1]);
2285         }
2286
2287         aff1 = isl_aff_scale(aff1, aff2->v->el[0]);
2288         aff1 = isl_aff_scale_down(aff1, aff2->v->el[1]);
2289
2290         if (neg) {
2291                 isl_int_neg(aff2->v->el[0], aff2->v->el[0]);
2292                 isl_int_neg(aff2->v->el[1], aff2->v->el[1]);
2293         }
2294
2295         isl_aff_free(aff2);
2296         return aff1;
2297 error:
2298         isl_aff_free(aff1);
2299         isl_aff_free(aff2);
2300         return NULL;
2301 }
2302
2303 static __isl_give isl_pw_aff *pw_aff_add(__isl_take isl_pw_aff *pwaff1,
2304         __isl_take isl_pw_aff *pwaff2)
2305 {
2306         return isl_pw_aff_on_shared_domain(pwaff1, pwaff2, &isl_aff_add);
2307 }
2308
2309 __isl_give isl_pw_aff *isl_pw_aff_add(__isl_take isl_pw_aff *pwaff1,
2310         __isl_take isl_pw_aff *pwaff2)
2311 {
2312         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2, &pw_aff_add);
2313 }
2314
2315 __isl_give isl_pw_aff *isl_pw_aff_union_add(__isl_take isl_pw_aff *pwaff1,
2316         __isl_take isl_pw_aff *pwaff2)
2317 {
2318         return isl_pw_aff_union_add_(pwaff1, pwaff2);
2319 }
2320
2321 static __isl_give isl_pw_aff *pw_aff_mul(__isl_take isl_pw_aff *pwaff1,
2322         __isl_take isl_pw_aff *pwaff2)
2323 {
2324         return isl_pw_aff_on_shared_domain(pwaff1, pwaff2, &isl_aff_mul);
2325 }
2326
2327 __isl_give isl_pw_aff *isl_pw_aff_mul(__isl_take isl_pw_aff *pwaff1,
2328         __isl_take isl_pw_aff *pwaff2)
2329 {
2330         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2, &pw_aff_mul);
2331 }
2332
2333 static __isl_give isl_pw_aff *pw_aff_div(__isl_take isl_pw_aff *pa1,
2334         __isl_take isl_pw_aff *pa2)
2335 {
2336         return isl_pw_aff_on_shared_domain(pa1, pa2, &isl_aff_div);
2337 }
2338
2339 /* Divide "pa1" by "pa2", assuming "pa2" is a piecewise constant.
2340  */
2341 __isl_give isl_pw_aff *isl_pw_aff_div(__isl_take isl_pw_aff *pa1,
2342         __isl_take isl_pw_aff *pa2)
2343 {
2344         int is_cst;
2345
2346         is_cst = isl_pw_aff_is_cst(pa2);
2347         if (is_cst < 0)
2348                 goto error;
2349         if (!is_cst)
2350                 isl_die(isl_pw_aff_get_ctx(pa2), isl_error_invalid,
2351                         "second argument should be a piecewise constant",
2352                         goto error);
2353         return isl_pw_aff_align_params_pw_pw_and(pa1, pa2, &pw_aff_div);
2354 error:
2355         isl_pw_aff_free(pa1);
2356         isl_pw_aff_free(pa2);
2357         return NULL;
2358 }
2359
2360 /* Compute the quotient of the integer division of "pa1" by "pa2"
2361  * with rounding towards zero.
2362  * "pa2" is assumed to be a piecewise constant.
2363  *
2364  * In particular, return
2365  *
2366  *      pa1 >= 0 ? floor(pa1/pa2) : ceil(pa1/pa2)
2367  *
2368  */
2369 __isl_give isl_pw_aff *isl_pw_aff_tdiv_q(__isl_take isl_pw_aff *pa1,
2370         __isl_take isl_pw_aff *pa2)
2371 {
2372         int is_cst;
2373         isl_set *cond;
2374         isl_pw_aff *f, *c;
2375
2376         is_cst = isl_pw_aff_is_cst(pa2);
2377         if (is_cst < 0)
2378                 goto error;
2379         if (!is_cst)
2380                 isl_die(isl_pw_aff_get_ctx(pa2), isl_error_invalid,
2381                         "second argument should be a piecewise constant",
2382                         goto error);
2383
2384         pa1 = isl_pw_aff_div(pa1, pa2);
2385
2386         cond = isl_pw_aff_nonneg_set(isl_pw_aff_copy(pa1));
2387         f = isl_pw_aff_floor(isl_pw_aff_copy(pa1));
2388         c = isl_pw_aff_ceil(pa1);
2389         return isl_pw_aff_cond(isl_set_indicator_function(cond), f, c);
2390 error:
2391         isl_pw_aff_free(pa1);
2392         isl_pw_aff_free(pa2);
2393         return NULL;
2394 }
2395
2396 /* Compute the remainder of the integer division of "pa1" by "pa2"
2397  * with rounding towards zero.
2398  * "pa2" is assumed to be a piecewise constant.
2399  *
2400  * In particular, return
2401  *
2402  *      pa1 - pa2 * (pa1 >= 0 ? floor(pa1/pa2) : ceil(pa1/pa2))
2403  *
2404  */
2405 __isl_give isl_pw_aff *isl_pw_aff_tdiv_r(__isl_take isl_pw_aff *pa1,
2406         __isl_take isl_pw_aff *pa2)
2407 {
2408         int is_cst;
2409         isl_pw_aff *res;
2410
2411         is_cst = isl_pw_aff_is_cst(pa2);
2412         if (is_cst < 0)
2413                 goto error;
2414         if (!is_cst)
2415                 isl_die(isl_pw_aff_get_ctx(pa2), isl_error_invalid,
2416                         "second argument should be a piecewise constant",
2417                         goto error);
2418         res = isl_pw_aff_tdiv_q(isl_pw_aff_copy(pa1), isl_pw_aff_copy(pa2));
2419         res = isl_pw_aff_mul(pa2, res);
2420         res = isl_pw_aff_sub(pa1, res);
2421         return res;
2422 error:
2423         isl_pw_aff_free(pa1);
2424         isl_pw_aff_free(pa2);
2425         return NULL;
2426 }
2427
2428 static __isl_give isl_pw_aff *pw_aff_min(__isl_take isl_pw_aff *pwaff1,
2429         __isl_take isl_pw_aff *pwaff2)
2430 {
2431         isl_set *le;
2432         isl_set *dom;
2433
2434         dom = isl_set_intersect(isl_pw_aff_domain(isl_pw_aff_copy(pwaff1)),
2435                                 isl_pw_aff_domain(isl_pw_aff_copy(pwaff2)));
2436         le = isl_pw_aff_le_set(isl_pw_aff_copy(pwaff1),
2437                                 isl_pw_aff_copy(pwaff2));
2438         dom = isl_set_subtract(dom, isl_set_copy(le));
2439         return isl_pw_aff_select(le, pwaff1, dom, pwaff2);
2440 }
2441
2442 __isl_give isl_pw_aff *isl_pw_aff_min(__isl_take isl_pw_aff *pwaff1,
2443         __isl_take isl_pw_aff *pwaff2)
2444 {
2445         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2, &pw_aff_min);
2446 }
2447
2448 static __isl_give isl_pw_aff *pw_aff_max(__isl_take isl_pw_aff *pwaff1,
2449         __isl_take isl_pw_aff *pwaff2)
2450 {
2451         isl_set *ge;
2452         isl_set *dom;
2453
2454         dom = isl_set_intersect(isl_pw_aff_domain(isl_pw_aff_copy(pwaff1)),
2455                                 isl_pw_aff_domain(isl_pw_aff_copy(pwaff2)));
2456         ge = isl_pw_aff_ge_set(isl_pw_aff_copy(pwaff1),
2457                                 isl_pw_aff_copy(pwaff2));
2458         dom = isl_set_subtract(dom, isl_set_copy(ge));
2459         return isl_pw_aff_select(ge, pwaff1, dom, pwaff2);
2460 }
2461
2462 __isl_give isl_pw_aff *isl_pw_aff_max(__isl_take isl_pw_aff *pwaff1,
2463         __isl_take isl_pw_aff *pwaff2)
2464 {
2465         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2, &pw_aff_max);
2466 }
2467
2468 static __isl_give isl_pw_aff *pw_aff_list_reduce(
2469         __isl_take isl_pw_aff_list *list,
2470         __isl_give isl_pw_aff *(*fn)(__isl_take isl_pw_aff *pwaff1,
2471                                         __isl_take isl_pw_aff *pwaff2))
2472 {
2473         int i;
2474         isl_ctx *ctx;
2475         isl_pw_aff *res;
2476
2477         if (!list)
2478                 return NULL;
2479
2480         ctx = isl_pw_aff_list_get_ctx(list);
2481         if (list->n < 1)
2482                 isl_die(ctx, isl_error_invalid,
2483                         "list should contain at least one element",
2484                         return isl_pw_aff_list_free(list));
2485
2486         res = isl_pw_aff_copy(list->p[0]);
2487         for (i = 1; i < list->n; ++i)
2488                 res = fn(res, isl_pw_aff_copy(list->p[i]));
2489
2490         isl_pw_aff_list_free(list);
2491         return res;
2492 }
2493
2494 /* Return an isl_pw_aff that maps each element in the intersection of the
2495  * domains of the elements of list to the minimal corresponding affine
2496  * expression.
2497  */
2498 __isl_give isl_pw_aff *isl_pw_aff_list_min(__isl_take isl_pw_aff_list *list)
2499 {
2500         return pw_aff_list_reduce(list, &isl_pw_aff_min);
2501 }
2502
2503 /* Return an isl_pw_aff that maps each element in the intersection of the
2504  * domains of the elements of list to the maximal corresponding affine
2505  * expression.
2506  */
2507 __isl_give isl_pw_aff *isl_pw_aff_list_max(__isl_take isl_pw_aff_list *list)
2508 {
2509         return pw_aff_list_reduce(list, &isl_pw_aff_max);
2510 }
2511
2512 #undef BASE
2513 #define BASE aff
2514
2515 #include <isl_multi_templ.c>
2516
2517 /* Create an isl_pw_multi_aff with the given isl_multi_aff on a universe
2518  * domain.
2519  */
2520 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_from_multi_aff(
2521         __isl_take isl_multi_aff *ma)
2522 {
2523         isl_set *dom = isl_set_universe(isl_multi_aff_get_domain_space(ma));
2524         return isl_pw_multi_aff_alloc(dom, ma);
2525 }
2526
2527 /* Create a piecewise multi-affine expression in the given space that maps each
2528  * input dimension to the corresponding output dimension.
2529  */
2530 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_identity(
2531         __isl_take isl_space *space)
2532 {
2533         return isl_pw_multi_aff_from_multi_aff(isl_multi_aff_identity(space));
2534 }
2535
2536 __isl_give isl_multi_aff *isl_multi_aff_add(__isl_take isl_multi_aff *maff1,
2537         __isl_take isl_multi_aff *maff2)
2538 {
2539         int i;
2540         isl_ctx *ctx;
2541
2542         maff1 = isl_multi_aff_cow(maff1);
2543         if (!maff1 || !maff2)
2544                 goto error;
2545
2546         ctx = isl_multi_aff_get_ctx(maff1);
2547         if (!isl_space_is_equal(maff1->space, maff2->space))
2548                 isl_die(ctx, isl_error_invalid,
2549                         "spaces don't match", goto error);
2550
2551         for (i = 0; i < maff1->n; ++i) {
2552                 maff1->p[i] = isl_aff_add(maff1->p[i],
2553                                             isl_aff_copy(maff2->p[i]));
2554                 if (!maff1->p[i])
2555                         goto error;
2556         }
2557
2558         isl_multi_aff_free(maff2);
2559         return maff1;
2560 error:
2561         isl_multi_aff_free(maff1);
2562         isl_multi_aff_free(maff2);
2563         return NULL;
2564 }
2565
2566 /* Given two multi-affine expressions A -> B and C -> D,
2567  * construct a multi-affine expression [A -> C] -> [B -> D].
2568  */
2569 __isl_give isl_multi_aff *isl_multi_aff_product(
2570         __isl_take isl_multi_aff *ma1, __isl_take isl_multi_aff *ma2)
2571 {
2572         int i;
2573         isl_aff *aff;
2574         isl_space *space;
2575         isl_multi_aff *res;
2576         int in1, in2, out1, out2;
2577
2578         in1 = isl_multi_aff_dim(ma1, isl_dim_in);
2579         in2 = isl_multi_aff_dim(ma2, isl_dim_in);
2580         out1 = isl_multi_aff_dim(ma1, isl_dim_out);
2581         out2 = isl_multi_aff_dim(ma2, isl_dim_out);
2582         space = isl_space_product(isl_multi_aff_get_space(ma1),
2583                                   isl_multi_aff_get_space(ma2));
2584         res = isl_multi_aff_alloc(isl_space_copy(space));
2585         space = isl_space_domain(space);
2586
2587         for (i = 0; i < out1; ++i) {
2588                 aff = isl_multi_aff_get_aff(ma1, i);
2589                 aff = isl_aff_insert_dims(aff, isl_dim_in, in1, in2);
2590                 aff = isl_aff_reset_domain_space(aff, isl_space_copy(space));
2591                 res = isl_multi_aff_set_aff(res, i, aff);
2592         }
2593
2594         for (i = 0; i < out2; ++i) {
2595                 aff = isl_multi_aff_get_aff(ma2, i);
2596                 aff = isl_aff_insert_dims(aff, isl_dim_in, 0, in1);
2597                 aff = isl_aff_reset_domain_space(aff, isl_space_copy(space));
2598                 res = isl_multi_aff_set_aff(res, out1 + i, aff);
2599         }
2600
2601         isl_space_free(space);
2602         isl_multi_aff_free(ma1);
2603         isl_multi_aff_free(ma2);
2604         return res;
2605 }
2606
2607 /* Exploit the equalities in "eq" to simplify the affine expressions.
2608  */
2609 static __isl_give isl_multi_aff *isl_multi_aff_substitute_equalities(
2610         __isl_take isl_multi_aff *maff, __isl_take isl_basic_set *eq)
2611 {
2612         int i;
2613
2614         maff = isl_multi_aff_cow(maff);
2615         if (!maff || !eq)
2616                 goto error;
2617
2618         for (i = 0; i < maff->n; ++i) {
2619                 maff->p[i] = isl_aff_substitute_equalities(maff->p[i],
2620                                                     isl_basic_set_copy(eq));
2621                 if (!maff->p[i])
2622                         goto error;
2623         }
2624
2625         isl_basic_set_free(eq);
2626         return maff;
2627 error:
2628         isl_basic_set_free(eq);
2629         isl_multi_aff_free(maff);
2630         return NULL;
2631 }
2632
2633 __isl_give isl_multi_aff *isl_multi_aff_scale(__isl_take isl_multi_aff *maff,
2634         isl_int f)
2635 {
2636         int i;
2637
2638         maff = isl_multi_aff_cow(maff);
2639         if (!maff)
2640                 return NULL;
2641
2642         for (i = 0; i < maff->n; ++i) {
2643                 maff->p[i] = isl_aff_scale(maff->p[i], f);
2644                 if (!maff->p[i])
2645                         return isl_multi_aff_free(maff);
2646         }
2647
2648         return maff;
2649 }
2650
2651 __isl_give isl_multi_aff *isl_multi_aff_add_on_domain(__isl_keep isl_set *dom,
2652         __isl_take isl_multi_aff *maff1, __isl_take isl_multi_aff *maff2)
2653 {
2654         maff1 = isl_multi_aff_add(maff1, maff2);
2655         maff1 = isl_multi_aff_gist(maff1, isl_set_copy(dom));
2656         return maff1;
2657 }
2658
2659 int isl_multi_aff_is_empty(__isl_keep isl_multi_aff *maff)
2660 {
2661         if (!maff)
2662                 return -1;
2663
2664         return 0;
2665 }
2666
2667 int isl_multi_aff_plain_is_equal(__isl_keep isl_multi_aff *maff1,
2668         __isl_keep isl_multi_aff *maff2)
2669 {
2670         int i;
2671         int equal;
2672
2673         if (!maff1 || !maff2)
2674                 return -1;
2675         if (maff1->n != maff2->n)
2676                 return 0;
2677         equal = isl_space_is_equal(maff1->space, maff2->space);
2678         if (equal < 0 || !equal)
2679                 return equal;
2680
2681         for (i = 0; i < maff1->n; ++i) {
2682                 equal = isl_aff_plain_is_equal(maff1->p[i], maff2->p[i]);
2683                 if (equal < 0 || !equal)
2684                         return equal;
2685         }
2686
2687         return 1;
2688 }
2689
2690 /* Return the set of domain elements where "ma1" is lexicographically
2691  * smaller than or equal to "ma2".
2692  */
2693 __isl_give isl_set *isl_multi_aff_lex_le_set(__isl_take isl_multi_aff *ma1,
2694         __isl_take isl_multi_aff *ma2)
2695 {
2696         return isl_multi_aff_lex_ge_set(ma2, ma1);
2697 }
2698
2699 /* Return the set of domain elements where "ma1" is lexicographically
2700  * greater than or equal to "ma2".
2701  */
2702 __isl_give isl_set *isl_multi_aff_lex_ge_set(__isl_take isl_multi_aff *ma1,
2703         __isl_take isl_multi_aff *ma2)
2704 {
2705         isl_space *space;
2706         isl_map *map1, *map2;
2707         isl_map *map, *ge;
2708
2709         map1 = isl_map_from_multi_aff(ma1);
2710         map2 = isl_map_from_multi_aff(ma2);
2711         map = isl_map_range_product(map1, map2);
2712         space = isl_space_range(isl_map_get_space(map));
2713         space = isl_space_domain(isl_space_unwrap(space));
2714         ge = isl_map_lex_ge(space);
2715         map = isl_map_intersect_range(map, isl_map_wrap(ge));
2716
2717         return isl_map_domain(map);
2718 }
2719
2720 #undef PW
2721 #define PW isl_pw_multi_aff
2722 #undef EL
2723 #define EL isl_multi_aff
2724 #undef EL_IS_ZERO
2725 #define EL_IS_ZERO is_empty
2726 #undef ZERO
2727 #define ZERO empty
2728 #undef IS_ZERO
2729 #define IS_ZERO is_empty
2730 #undef FIELD
2731 #define FIELD maff
2732 #undef DEFAULT_IS_ZERO
2733 #define DEFAULT_IS_ZERO 0
2734
2735 #define NO_NEG
2736 #define NO_EVAL
2737 #define NO_OPT
2738 #define NO_INVOLVES_DIMS
2739 #define NO_MOVE_DIMS
2740 #define NO_INSERT_DIMS
2741 #define NO_LIFT
2742 #define NO_MORPH
2743
2744 #include <isl_pw_templ.c>
2745
2746 #undef UNION
2747 #define UNION isl_union_pw_multi_aff
2748 #undef PART
2749 #define PART isl_pw_multi_aff
2750 #undef PARTS
2751 #define PARTS pw_multi_aff
2752 #define ALIGN_DOMAIN
2753
2754 #define NO_EVAL
2755
2756 #include <isl_union_templ.c>
2757
2758 /* Given a function "cmp" that returns the set of elements where
2759  * "ma1" is "better" than "ma2", return the intersection of this
2760  * set with "dom1" and "dom2".
2761  */
2762 static __isl_give isl_set *shared_and_better(__isl_keep isl_set *dom1,
2763         __isl_keep isl_set *dom2, __isl_keep isl_multi_aff *ma1,
2764         __isl_keep isl_multi_aff *ma2,
2765         __isl_give isl_set *(*cmp)(__isl_take isl_multi_aff *ma1,
2766                                     __isl_take isl_multi_aff *ma2))
2767 {
2768         isl_set *common;
2769         isl_set *better;
2770         int is_empty;
2771
2772         common = isl_set_intersect(isl_set_copy(dom1), isl_set_copy(dom2));
2773         is_empty = isl_set_plain_is_empty(common);
2774         if (is_empty >= 0 && is_empty)
2775                 return common;
2776         if (is_empty < 0)
2777                 return isl_set_free(common);
2778         better = cmp(isl_multi_aff_copy(ma1), isl_multi_aff_copy(ma2));
2779         better = isl_set_intersect(common, better);
2780
2781         return better;
2782 }
2783
2784 /* Given a function "cmp" that returns the set of elements where
2785  * "ma1" is "better" than "ma2", return a piecewise multi affine
2786  * expression defined on the union of the definition domains
2787  * of "pma1" and "pma2" that maps to the "best" of "pma1" and
2788  * "pma2" on each cell.  If only one of the two input functions
2789  * is defined on a given cell, then it is considered the best.
2790  */
2791 static __isl_give isl_pw_multi_aff *pw_multi_aff_union_opt(
2792         __isl_take isl_pw_multi_aff *pma1,
2793         __isl_take isl_pw_multi_aff *pma2,
2794         __isl_give isl_set *(*cmp)(__isl_take isl_multi_aff *ma1,
2795                                     __isl_take isl_multi_aff *ma2))
2796 {
2797         int i, j, n;
2798         isl_pw_multi_aff *res = NULL;
2799         isl_ctx *ctx;
2800         isl_set *set = NULL;
2801
2802         if (!pma1 || !pma2)
2803                 goto error;
2804
2805         ctx = isl_space_get_ctx(pma1->dim);
2806         if (!isl_space_is_equal(pma1->dim, pma2->dim))
2807                 isl_die(ctx, isl_error_invalid,
2808                         "arguments should live in the same space", goto error);
2809
2810         if (isl_pw_multi_aff_is_empty(pma1)) {
2811                 isl_pw_multi_aff_free(pma1);
2812                 return pma2;
2813         }
2814
2815         if (isl_pw_multi_aff_is_empty(pma2)) {
2816                 isl_pw_multi_aff_free(pma2);
2817                 return pma1;
2818         }
2819
2820         n = 2 * (pma1->n + 1) * (pma2->n + 1);
2821         res = isl_pw_multi_aff_alloc_size(isl_space_copy(pma1->dim), n);
2822
2823         for (i = 0; i < pma1->n; ++i) {
2824                 set = isl_set_copy(pma1->p[i].set);
2825                 for (j = 0; j < pma2->n; ++j) {
2826                         isl_set *better;
2827                         int is_empty;
2828
2829                         better = shared_and_better(pma2->p[j].set,
2830                                         pma1->p[i].set, pma2->p[j].maff,
2831                                         pma1->p[i].maff, cmp);
2832                         is_empty = isl_set_plain_is_empty(better);
2833                         if (is_empty < 0 || is_empty) {
2834                                 isl_set_free(better);
2835                                 if (is_empty < 0)
2836                                         goto error;
2837                                 continue;
2838                         }
2839                         set = isl_set_subtract(set, isl_set_copy(better));
2840
2841                         res = isl_pw_multi_aff_add_piece(res, better,
2842                                         isl_multi_aff_copy(pma2->p[j].maff));
2843                 }
2844                 res = isl_pw_multi_aff_add_piece(res, set,
2845                                         isl_multi_aff_copy(pma1->p[i].maff));
2846         }
2847
2848         for (j = 0; j < pma2->n; ++j) {
2849                 set = isl_set_copy(pma2->p[j].set);
2850                 for (i = 0; i < pma1->n; ++i)
2851                         set = isl_set_subtract(set,
2852                                         isl_set_copy(pma1->p[i].set));
2853                 res = isl_pw_multi_aff_add_piece(res, set,
2854                                         isl_multi_aff_copy(pma2->p[j].maff));
2855         }
2856
2857         isl_pw_multi_aff_free(pma1);
2858         isl_pw_multi_aff_free(pma2);
2859
2860         return res;
2861 error:
2862         isl_pw_multi_aff_free(pma1);
2863         isl_pw_multi_aff_free(pma2);
2864         isl_set_free(set);
2865         return isl_pw_multi_aff_free(res);
2866 }
2867
2868 static __isl_give isl_pw_multi_aff *pw_multi_aff_union_lexmax(
2869         __isl_take isl_pw_multi_aff *pma1,
2870         __isl_take isl_pw_multi_aff *pma2)
2871 {
2872         return pw_multi_aff_union_opt(pma1, pma2, &isl_multi_aff_lex_ge_set);
2873 }
2874
2875 /* Given two piecewise multi affine expressions, return a piecewise
2876  * multi-affine expression defined on the union of the definition domains
2877  * of the inputs that is equal to the lexicographic maximum of the two
2878  * inputs on each cell.  If only one of the two inputs is defined on
2879  * a given cell, then it is considered to be the maximum.
2880  */
2881 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_union_lexmax(
2882         __isl_take isl_pw_multi_aff *pma1,
2883         __isl_take isl_pw_multi_aff *pma2)
2884 {
2885         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
2886                                                     &pw_multi_aff_union_lexmax);
2887 }
2888
2889 static __isl_give isl_pw_multi_aff *pw_multi_aff_union_lexmin(
2890         __isl_take isl_pw_multi_aff *pma1,
2891         __isl_take isl_pw_multi_aff *pma2)
2892 {
2893         return pw_multi_aff_union_opt(pma1, pma2, &isl_multi_aff_lex_le_set);
2894 }
2895
2896 /* Given two piecewise multi affine expressions, return a piecewise
2897  * multi-affine expression defined on the union of the definition domains
2898  * of the inputs that is equal to the lexicographic minimum of the two
2899  * inputs on each cell.  If only one of the two inputs is defined on
2900  * a given cell, then it is considered to be the minimum.
2901  */
2902 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_union_lexmin(
2903         __isl_take isl_pw_multi_aff *pma1,
2904         __isl_take isl_pw_multi_aff *pma2)
2905 {
2906         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
2907                                                     &pw_multi_aff_union_lexmin);
2908 }
2909
2910 static __isl_give isl_pw_multi_aff *pw_multi_aff_add(
2911         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
2912 {
2913         return isl_pw_multi_aff_on_shared_domain(pma1, pma2,
2914                                                 &isl_multi_aff_add);
2915 }
2916
2917 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_add(
2918         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
2919 {
2920         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
2921                                                 &pw_multi_aff_add);
2922 }
2923
2924 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_union_add(
2925         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
2926 {
2927         return isl_pw_multi_aff_union_add_(pma1, pma2);
2928 }
2929
2930 /* Given two piecewise multi-affine expressions A -> B and C -> D,
2931  * construct a piecewise multi-affine expression [A -> C] -> [B -> D].
2932  */
2933 static __isl_give isl_pw_multi_aff *pw_multi_aff_product(
2934         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
2935 {
2936         int i, j, n;
2937         isl_space *space;
2938         isl_pw_multi_aff *res;
2939
2940         if (!pma1 || !pma2)
2941                 goto error;
2942
2943         n = pma1->n * pma2->n;
2944         space = isl_space_product(isl_space_copy(pma1->dim),
2945                                   isl_space_copy(pma2->dim));
2946         res = isl_pw_multi_aff_alloc_size(space, n);
2947
2948         for (i = 0; i < pma1->n; ++i) {
2949                 for (j = 0; j < pma2->n; ++j) {
2950                         isl_set *domain;
2951                         isl_multi_aff *ma;
2952
2953                         domain = isl_set_product(isl_set_copy(pma1->p[i].set),
2954                                                  isl_set_copy(pma2->p[j].set));
2955                         ma = isl_multi_aff_product(
2956                                         isl_multi_aff_copy(pma1->p[i].maff),
2957                                         isl_multi_aff_copy(pma2->p[i].maff));
2958                         res = isl_pw_multi_aff_add_piece(res, domain, ma);
2959                 }
2960         }
2961
2962         isl_pw_multi_aff_free(pma1);
2963         isl_pw_multi_aff_free(pma2);
2964         return res;
2965 error:
2966         isl_pw_multi_aff_free(pma1);
2967         isl_pw_multi_aff_free(pma2);
2968         return NULL;
2969 }
2970
2971 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_product(
2972         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
2973 {
2974         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
2975                                                 &pw_multi_aff_product);
2976 }
2977
2978 /* Construct a map mapping the domain of the piecewise multi-affine expression
2979  * to its range, with each dimension in the range equated to the
2980  * corresponding affine expression on its cell.
2981  */
2982 __isl_give isl_map *isl_map_from_pw_multi_aff(__isl_take isl_pw_multi_aff *pma)
2983 {
2984         int i;
2985         isl_map *map;
2986
2987         if (!pma)
2988                 return NULL;
2989
2990         map = isl_map_empty(isl_pw_multi_aff_get_space(pma));
2991
2992         for (i = 0; i < pma->n; ++i) {
2993                 isl_multi_aff *maff;
2994                 isl_basic_map *bmap;
2995                 isl_map *map_i;
2996
2997                 maff = isl_multi_aff_copy(pma->p[i].maff);
2998                 bmap = isl_basic_map_from_multi_aff(maff);
2999                 map_i = isl_map_from_basic_map(bmap);
3000                 map_i = isl_map_intersect_domain(map_i,
3001                                                 isl_set_copy(pma->p[i].set));
3002                 map = isl_map_union_disjoint(map, map_i);
3003         }
3004
3005         isl_pw_multi_aff_free(pma);
3006         return map;
3007 }
3008
3009 __isl_give isl_set *isl_set_from_pw_multi_aff(__isl_take isl_pw_multi_aff *pma)
3010 {
3011         if (!pma)
3012                 return NULL;
3013
3014         if (!isl_space_is_set(pma->dim))
3015                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
3016                         "isl_pw_multi_aff cannot be converted into an isl_set",
3017                         return isl_pw_multi_aff_free(pma));
3018
3019         return isl_map_from_pw_multi_aff(pma);
3020 }
3021
3022 /* Given a basic map with a single output dimension that is defined
3023  * in terms of the parameters and input dimensions using an equality,
3024  * extract an isl_aff that expresses the output dimension in terms
3025  * of the parameters and input dimensions.
3026  *
3027  * Since some applications expect the result of isl_pw_multi_aff_from_map
3028  * to only contain integer affine expressions, we compute the floor
3029  * of the expression before returning.
3030  *
3031  * This function shares some similarities with
3032  * isl_basic_map_has_defining_equality and isl_constraint_get_bound.
3033  */
3034 static __isl_give isl_aff *extract_isl_aff_from_basic_map(
3035         __isl_take isl_basic_map *bmap)
3036 {
3037         int i;
3038         unsigned offset;
3039         unsigned total;
3040         isl_local_space *ls;
3041         isl_aff *aff;
3042
3043         if (!bmap)
3044                 return NULL;
3045         if (isl_basic_map_dim(bmap, isl_dim_out) != 1)
3046                 isl_die(isl_basic_map_get_ctx(bmap), isl_error_invalid,
3047                         "basic map should have a single output dimension",
3048                         goto error);
3049         offset = isl_basic_map_offset(bmap, isl_dim_out);
3050         total = isl_basic_map_total_dim(bmap);
3051         for (i = 0; i < bmap->n_eq; ++i) {
3052                 if (isl_int_is_zero(bmap->eq[i][offset]))
3053                         continue;
3054                 if (isl_seq_first_non_zero(bmap->eq[i] + offset + 1,
3055                                            1 + total - (offset + 1)) != -1)
3056                         continue;
3057                 break;
3058         }
3059         if (i >= bmap->n_eq)
3060                 isl_die(isl_basic_map_get_ctx(bmap), isl_error_invalid,
3061                         "unable to find suitable equality", goto error);
3062         ls = isl_basic_map_get_local_space(bmap);
3063         aff = isl_aff_alloc(isl_local_space_domain(ls));
3064         if (!aff)
3065                 goto error;
3066         if (isl_int_is_neg(bmap->eq[i][offset]))
3067                 isl_seq_cpy(aff->v->el + 1, bmap->eq[i], offset);
3068         else
3069                 isl_seq_neg(aff->v->el + 1, bmap->eq[i], offset);
3070         isl_seq_clr(aff->v->el + 1 + offset, aff->v->size - (1 + offset));
3071         isl_int_abs(aff->v->el[0], bmap->eq[i][offset]);
3072         isl_basic_map_free(bmap);
3073
3074         aff = isl_aff_remove_unused_divs(aff);
3075         aff = isl_aff_floor(aff);
3076         return aff;
3077 error:
3078         isl_basic_map_free(bmap);
3079         return NULL;
3080 }
3081
3082 /* Given a basic map where each output dimension is defined
3083  * in terms of the parameters and input dimensions using an equality,
3084  * extract an isl_multi_aff that expresses the output dimensions in terms
3085  * of the parameters and input dimensions.
3086  */
3087 static __isl_give isl_multi_aff *extract_isl_multi_aff_from_basic_map(
3088         __isl_take isl_basic_map *bmap)
3089 {
3090         int i;
3091         unsigned n_out;
3092         isl_multi_aff *ma;
3093
3094         if (!bmap)
3095                 return NULL;
3096
3097         ma = isl_multi_aff_alloc(isl_basic_map_get_space(bmap));
3098         n_out = isl_basic_map_dim(bmap, isl_dim_out);
3099
3100         for (i = 0; i < n_out; ++i) {
3101                 isl_basic_map *bmap_i;
3102                 isl_aff *aff;
3103
3104                 bmap_i = isl_basic_map_copy(bmap);
3105                 bmap_i = isl_basic_map_project_out(bmap_i, isl_dim_out,
3106                                                         i + 1, n_out - (1 + i));
3107                 bmap_i = isl_basic_map_project_out(bmap_i, isl_dim_out, 0, i);
3108                 aff = extract_isl_aff_from_basic_map(bmap_i);
3109                 ma = isl_multi_aff_set_aff(ma, i, aff);
3110         }
3111
3112         isl_basic_map_free(bmap);
3113
3114         return ma;
3115 }
3116
3117 /* Create an isl_pw_multi_aff that is equivalent to
3118  * isl_map_intersect_domain(isl_map_from_basic_map(bmap), domain).
3119  * The given basic map is such that each output dimension is defined
3120  * in terms of the parameters and input dimensions using an equality.
3121  */
3122 static __isl_give isl_pw_multi_aff *plain_pw_multi_aff_from_map(
3123         __isl_take isl_set *domain, __isl_take isl_basic_map *bmap)
3124 {
3125         isl_multi_aff *ma;
3126
3127         ma = extract_isl_multi_aff_from_basic_map(bmap);
3128         return isl_pw_multi_aff_alloc(domain, ma);
3129 }
3130
3131 /* Try and create an isl_pw_multi_aff that is equivalent to the given isl_map.
3132  * This obviously only works if the input "map" is single-valued.
3133  * If so, we compute the lexicographic minimum of the image in the form
3134  * of an isl_pw_multi_aff.  Since the image is unique, it is equal
3135  * to its lexicographic minimum.
3136  * If the input is not single-valued, we produce an error.
3137  *
3138  * As a special case, we first check if all output dimensions are uniquely
3139  * defined in terms of the parameters and input dimensions over the entire
3140  * domain.  If so, we extract the desired isl_pw_multi_aff directly
3141  * from the affine hull of "map" and its domain.
3142  */
3143 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_from_map(__isl_take isl_map *map)
3144 {
3145         int i;
3146         int sv;
3147         isl_pw_multi_aff *pma;
3148         isl_basic_map *hull;
3149
3150         if (!map)
3151                 return NULL;
3152
3153         hull = isl_map_affine_hull(isl_map_copy(map));
3154         sv = isl_basic_map_plain_is_single_valued(hull);
3155         if (sv >= 0 && sv)
3156                 return plain_pw_multi_aff_from_map(isl_map_domain(map), hull);
3157         isl_basic_map_free(hull);
3158         if (sv < 0)
3159                 goto error;
3160
3161         sv = isl_map_is_single_valued(map);
3162         if (sv < 0)
3163                 goto error;
3164         if (!sv)
3165                 isl_die(isl_map_get_ctx(map), isl_error_invalid,
3166                         "map is not single-valued", goto error);
3167         map = isl_map_make_disjoint(map);
3168         if (!map)
3169                 return NULL;
3170
3171         pma = isl_pw_multi_aff_empty(isl_map_get_space(map));
3172
3173         for (i = 0; i < map->n; ++i) {
3174                 isl_pw_multi_aff *pma_i;
3175                 isl_basic_map *bmap;
3176                 bmap = isl_basic_map_copy(map->p[i]);
3177                 pma_i = isl_basic_map_lexmin_pw_multi_aff(bmap);
3178                 pma = isl_pw_multi_aff_add_disjoint(pma, pma_i);
3179         }
3180
3181         isl_map_free(map);
3182         return pma;
3183 error:
3184         isl_map_free(map);
3185         return NULL;
3186 }
3187
3188 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_from_set(__isl_take isl_set *set)
3189 {
3190         return isl_pw_multi_aff_from_map(set);
3191 }
3192
3193 /* Return the piecewise affine expression "set ? 1 : 0".
3194  */
3195 __isl_give isl_pw_aff *isl_set_indicator_function(__isl_take isl_set *set)
3196 {
3197         isl_pw_aff *pa;
3198         isl_space *space = isl_set_get_space(set);
3199         isl_local_space *ls = isl_local_space_from_space(space);
3200         isl_aff *zero = isl_aff_zero_on_domain(isl_local_space_copy(ls));
3201         isl_aff *one = isl_aff_zero_on_domain(ls);
3202
3203         one = isl_aff_add_constant_si(one, 1);
3204         pa = isl_pw_aff_alloc(isl_set_copy(set), one);
3205         set = isl_set_complement(set);
3206         pa = isl_pw_aff_add_disjoint(pa, isl_pw_aff_alloc(set, zero));
3207
3208         return pa;
3209 }
3210
3211 /* Plug in "subs" for dimension "type", "pos" of "aff".
3212  *
3213  * Let i be the dimension to replace and let "subs" be of the form
3214  *
3215  *      f/d
3216  *
3217  * and "aff" of the form
3218  *
3219  *      (a i + g)/m
3220  *
3221  * The result is
3222  *
3223  *      (a f + d g')/(m d)
3224  *
3225  * where g' is the result of plugging in "subs" in each of the integer
3226  * divisions in g.
3227  */
3228 __isl_give isl_aff *isl_aff_substitute(__isl_take isl_aff *aff,
3229         enum isl_dim_type type, unsigned pos, __isl_keep isl_aff *subs)
3230 {
3231         isl_ctx *ctx;
3232         isl_int v;
3233
3234         aff = isl_aff_cow(aff);
3235         if (!aff || !subs)
3236                 return isl_aff_free(aff);
3237
3238         ctx = isl_aff_get_ctx(aff);
3239         if (!isl_space_is_equal(aff->ls->dim, subs->ls->dim))
3240                 isl_die(ctx, isl_error_invalid,
3241                         "spaces don't match", return isl_aff_free(aff));
3242         if (isl_local_space_dim(subs->ls, isl_dim_div) != 0)
3243                 isl_die(ctx, isl_error_unsupported,
3244                         "cannot handle divs yet", return isl_aff_free(aff));
3245
3246         aff->ls = isl_local_space_substitute(aff->ls, type, pos, subs);
3247         if (!aff->ls)
3248                 return isl_aff_free(aff);
3249
3250         aff->v = isl_vec_cow(aff->v);
3251         if (!aff->v)
3252                 return isl_aff_free(aff);
3253
3254         pos += isl_local_space_offset(aff->ls, type);
3255
3256         isl_int_init(v);
3257         isl_seq_substitute(aff->v->el, pos, subs->v->el,
3258                             aff->v->size, subs->v->size, v);
3259         isl_int_clear(v);
3260
3261         return aff;
3262 }
3263
3264 /* Plug in "subs" for dimension "type", "pos" in each of the affine
3265  * expressions in "maff".
3266  */
3267 __isl_give isl_multi_aff *isl_multi_aff_substitute(
3268         __isl_take isl_multi_aff *maff, enum isl_dim_type type, unsigned pos,
3269         __isl_keep isl_aff *subs)
3270 {
3271         int i;
3272
3273         maff = isl_multi_aff_cow(maff);
3274         if (!maff || !subs)
3275                 return isl_multi_aff_free(maff);
3276
3277         if (type == isl_dim_in)
3278                 type = isl_dim_set;
3279
3280         for (i = 0; i < maff->n; ++i) {
3281                 maff->p[i] = isl_aff_substitute(maff->p[i], type, pos, subs);
3282                 if (!maff->p[i])
3283                         return isl_multi_aff_free(maff);
3284         }
3285
3286         return maff;
3287 }
3288
3289 /* Plug in "subs" for dimension "type", "pos" of "pma".
3290  *
3291  * pma is of the form
3292  *
3293  *      A_i(v) -> M_i(v)
3294  *
3295  * while subs is of the form
3296  *
3297  *      v' = B_j(v) -> S_j
3298  *
3299  * Each pair i,j such that C_ij = A_i \cap B_i is non-empty
3300  * has a contribution in the result, in particular
3301  *
3302  *      C_ij(S_j) -> M_i(S_j)
3303  *
3304  * Note that plugging in S_j in C_ij may also result in an empty set
3305  * and this contribution should simply be discarded.
3306  */
3307 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_substitute(
3308         __isl_take isl_pw_multi_aff *pma, enum isl_dim_type type, unsigned pos,
3309         __isl_keep isl_pw_aff *subs)
3310 {
3311         int i, j, n;
3312         isl_pw_multi_aff *res;
3313
3314         if (!pma || !subs)
3315                 return isl_pw_multi_aff_free(pma);
3316
3317         n = pma->n * subs->n;
3318         res = isl_pw_multi_aff_alloc_size(isl_space_copy(pma->dim), n);
3319
3320         for (i = 0; i < pma->n; ++i) {
3321                 for (j = 0; j < subs->n; ++j) {
3322                         isl_set *common;
3323                         isl_multi_aff *res_ij;
3324                         common = isl_set_intersect(
3325                                         isl_set_copy(pma->p[i].set),
3326                                         isl_set_copy(subs->p[j].set));
3327                         common = isl_set_substitute(common,
3328                                         type, pos, subs->p[j].aff);
3329                         if (isl_set_plain_is_empty(common)) {
3330                                 isl_set_free(common);
3331                                 continue;
3332                         }
3333
3334                         res_ij = isl_multi_aff_substitute(
3335                                         isl_multi_aff_copy(pma->p[i].maff),
3336                                         type, pos, subs->p[j].aff);
3337
3338                         res = isl_pw_multi_aff_add_piece(res, common, res_ij);
3339                 }
3340         }
3341
3342         isl_pw_multi_aff_free(pma);
3343         return res;
3344 }
3345
3346 /* Extend the local space of "dst" to include the divs
3347  * in the local space of "src".
3348  */
3349 __isl_give isl_aff *isl_aff_align_divs(__isl_take isl_aff *dst,
3350         __isl_keep isl_aff *src)
3351 {
3352         isl_ctx *ctx;
3353         int *exp1 = NULL;
3354         int *exp2 = NULL;
3355         isl_mat *div;
3356
3357         if (!src || !dst)
3358                 return isl_aff_free(dst);
3359
3360         ctx = isl_aff_get_ctx(src);
3361         if (!isl_space_is_equal(src->ls->dim, dst->ls->dim))
3362                 isl_die(ctx, isl_error_invalid,
3363                         "spaces don't match", goto error);
3364
3365         if (src->ls->div->n_row == 0)
3366                 return dst;
3367
3368         exp1 = isl_alloc_array(ctx, int, src->ls->div->n_row);
3369         exp2 = isl_alloc_array(ctx, int, dst->ls->div->n_row);
3370         if (!exp1 || !exp2)
3371                 goto error;
3372
3373         div = isl_merge_divs(src->ls->div, dst->ls->div, exp1, exp2);
3374         dst = isl_aff_expand_divs(dst, div, exp2);
3375         free(exp1);
3376         free(exp2);
3377
3378         return dst;
3379 error:
3380         free(exp1);
3381         free(exp2);
3382         return isl_aff_free(dst);
3383 }
3384
3385 /* Adjust the local spaces of the affine expressions in "maff"
3386  * such that they all have the save divs.
3387  */
3388 __isl_give isl_multi_aff *isl_multi_aff_align_divs(
3389         __isl_take isl_multi_aff *maff)
3390 {
3391         int i;
3392
3393         if (!maff)
3394                 return NULL;
3395         if (maff->n == 0)
3396                 return maff;
3397         maff = isl_multi_aff_cow(maff);
3398         if (!maff)
3399                 return NULL;
3400
3401         for (i = 1; i < maff->n; ++i)
3402                 maff->p[0] = isl_aff_align_divs(maff->p[0], maff->p[i]);
3403         for (i = 1; i < maff->n; ++i) {
3404                 maff->p[i] = isl_aff_align_divs(maff->p[i], maff->p[0]);
3405                 if (!maff->p[i])
3406                         return isl_multi_aff_free(maff);
3407         }
3408
3409         return maff;
3410 }
3411
3412 __isl_give isl_aff *isl_aff_lift(__isl_take isl_aff *aff)
3413 {
3414         aff = isl_aff_cow(aff);
3415         if (!aff)
3416                 return NULL;
3417
3418         aff->ls = isl_local_space_lift(aff->ls);
3419         if (!aff->ls)
3420                 return isl_aff_free(aff);
3421
3422         return aff;
3423 }
3424
3425 /* Lift "maff" to a space with extra dimensions such that the result
3426  * has no more existentially quantified variables.
3427  * If "ls" is not NULL, then *ls is assigned the local space that lies
3428  * at the basis of the lifting applied to "maff".
3429  */
3430 __isl_give isl_multi_aff *isl_multi_aff_lift(__isl_take isl_multi_aff *maff,
3431         __isl_give isl_local_space **ls)
3432 {
3433         int i;
3434         isl_space *space;
3435         unsigned n_div;
3436
3437         if (ls)
3438                 *ls = NULL;
3439
3440         if (!maff)
3441                 return NULL;
3442
3443         if (maff->n == 0) {
3444                 if (ls) {
3445                         isl_space *space = isl_multi_aff_get_domain_space(maff);
3446                         *ls = isl_local_space_from_space(space);
3447                         if (!*ls)
3448                                 return isl_multi_aff_free(maff);
3449                 }
3450                 return maff;
3451         }
3452
3453         maff = isl_multi_aff_cow(maff);
3454         maff = isl_multi_aff_align_divs(maff);
3455         if (!maff)
3456                 return NULL;
3457
3458         n_div = isl_aff_dim(maff->p[0], isl_dim_div);
3459         space = isl_multi_aff_get_space(maff);
3460         space = isl_space_lift(isl_space_domain(space), n_div);
3461         space = isl_space_extend_domain_with_range(space,
3462                                                 isl_multi_aff_get_space(maff));
3463         if (!space)
3464                 return isl_multi_aff_free(maff);
3465         isl_space_free(maff->space);
3466         maff->space = space;
3467
3468         if (ls) {
3469                 *ls = isl_aff_get_domain_local_space(maff->p[0]);
3470                 if (!*ls)
3471                         return isl_multi_aff_free(maff);
3472         }
3473
3474         for (i = 0; i < maff->n; ++i) {
3475                 maff->p[i] = isl_aff_lift(maff->p[i]);
3476                 if (!maff->p[i])
3477                         goto error;
3478         }
3479
3480         return maff;
3481 error:
3482         if (ls)
3483                 isl_local_space_free(*ls);
3484         return isl_multi_aff_free(maff);
3485 }
3486
3487
3488 /* Extract an isl_pw_aff corresponding to output dimension "pos" of "pma".
3489  */
3490 __isl_give isl_pw_aff *isl_pw_multi_aff_get_pw_aff(
3491         __isl_keep isl_pw_multi_aff *pma, int pos)
3492 {
3493         int i;
3494         int n_out;
3495         isl_space *space;
3496         isl_pw_aff *pa;
3497
3498         if (!pma)
3499                 return NULL;
3500
3501         n_out = isl_pw_multi_aff_dim(pma, isl_dim_out);
3502         if (pos < 0 || pos >= n_out)
3503                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
3504                         "index out of bounds", return NULL);
3505
3506         space = isl_pw_multi_aff_get_space(pma);
3507         space = isl_space_drop_dims(space, isl_dim_out,
3508                                     pos + 1, n_out - pos - 1);
3509         space = isl_space_drop_dims(space, isl_dim_out, 0, pos);
3510
3511         pa = isl_pw_aff_alloc_size(space, pma->n);
3512         for (i = 0; i < pma->n; ++i) {
3513                 isl_aff *aff;
3514                 aff = isl_multi_aff_get_aff(pma->p[i].maff, pos);
3515                 pa = isl_pw_aff_add_piece(pa, isl_set_copy(pma->p[i].set), aff);
3516         }
3517
3518         return pa;
3519 }
3520
3521 /* Return an isl_pw_multi_aff with the given "set" as domain and
3522  * an unnamed zero-dimensional range.
3523  */
3524 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_from_domain(
3525         __isl_take isl_set *set)
3526 {
3527         isl_multi_aff *ma;
3528         isl_space *space;
3529
3530         space = isl_set_get_space(set);
3531         space = isl_space_from_domain(space);
3532         ma = isl_multi_aff_zero(space);
3533         return isl_pw_multi_aff_alloc(set, ma);
3534 }
3535
3536 /* Add an isl_pw_multi_aff with the given "set" as domain and
3537  * an unnamed zero-dimensional range to *user.
3538  */
3539 static int add_pw_multi_aff_from_domain(__isl_take isl_set *set, void *user)
3540 {
3541         isl_union_pw_multi_aff **upma = user;
3542         isl_pw_multi_aff *pma;
3543
3544         pma = isl_pw_multi_aff_from_domain(set);
3545         *upma = isl_union_pw_multi_aff_add_pw_multi_aff(*upma, pma);
3546
3547         return 0;
3548 }
3549
3550 /* Return an isl_union_pw_multi_aff with the given "uset" as domain and
3551  * an unnamed zero-dimensional range.
3552  */
3553 __isl_give isl_union_pw_multi_aff *isl_union_pw_multi_aff_from_domain(
3554         __isl_take isl_union_set *uset)
3555 {
3556         isl_space *space;
3557         isl_union_pw_multi_aff *upma;
3558
3559         if (!uset)
3560                 return NULL;
3561
3562         space = isl_union_set_get_space(uset);
3563         upma = isl_union_pw_multi_aff_empty(space);
3564
3565         if (isl_union_set_foreach_set(uset,
3566                                     &add_pw_multi_aff_from_domain, &upma) < 0)
3567                 goto error;
3568
3569         isl_union_set_free(uset);
3570         return upma;
3571 error:
3572         isl_union_set_free(uset);
3573         isl_union_pw_multi_aff_free(upma);
3574         return NULL;
3575 }
3576
3577 /* Convert "pma" to an isl_map and add it to *umap.
3578  */
3579 static int map_from_pw_multi_aff(__isl_take isl_pw_multi_aff *pma, void *user)
3580 {
3581         isl_union_map **umap = user;
3582         isl_map *map;
3583
3584         map = isl_map_from_pw_multi_aff(pma);
3585         *umap = isl_union_map_add_map(*umap, map);
3586
3587         return 0;
3588 }
3589
3590 /* Construct a union map mapping the domain of the union
3591  * piecewise multi-affine expression to its range, with each dimension
3592  * in the range equated to the corresponding affine expression on its cell.
3593  */
3594 __isl_give isl_union_map *isl_union_map_from_union_pw_multi_aff(
3595         __isl_take isl_union_pw_multi_aff *upma)
3596 {
3597         isl_space *space;
3598         isl_union_map *umap;
3599
3600         if (!upma)
3601                 return NULL;
3602
3603         space = isl_union_pw_multi_aff_get_space(upma);
3604         umap = isl_union_map_empty(space);
3605
3606         if (isl_union_pw_multi_aff_foreach_pw_multi_aff(upma,
3607                                         &map_from_pw_multi_aff, &umap) < 0)
3608                 goto error;
3609
3610         isl_union_pw_multi_aff_free(upma);
3611         return umap;
3612 error:
3613         isl_union_pw_multi_aff_free(upma);
3614         isl_union_map_free(umap);
3615         return NULL;
3616 }
3617
3618 /* Local data for bin_entry and the callback "fn".
3619  */
3620 struct isl_union_pw_multi_aff_bin_data {
3621         isl_union_pw_multi_aff *upma2;
3622         isl_union_pw_multi_aff *res;
3623         isl_pw_multi_aff *pma;
3624         int (*fn)(void **entry, void *user);
3625 };
3626
3627 /* Given an isl_pw_multi_aff from upma1, store it in data->pma
3628  * and call data->fn for each isl_pw_multi_aff in data->upma2.
3629  */
3630 static int bin_entry(void **entry, void *user)
3631 {
3632         struct isl_union_pw_multi_aff_bin_data *data = user;
3633         isl_pw_multi_aff *pma = *entry;
3634
3635         data->pma = pma;
3636         if (isl_hash_table_foreach(data->upma2->dim->ctx, &data->upma2->table,
3637                                    data->fn, data) < 0)
3638                 return -1;
3639
3640         return 0;
3641 }
3642
3643 /* Call "fn" on each pair of isl_pw_multi_affs in "upma1" and "upma2".
3644  * The isl_pw_multi_aff from upma1 is stored in data->pma (where data is
3645  * passed as user field) and the isl_pw_multi_aff from upma2 is available
3646  * as *entry.  The callback should adjust data->res if desired.
3647  */
3648 static __isl_give isl_union_pw_multi_aff *bin_op(
3649         __isl_take isl_union_pw_multi_aff *upma1,
3650         __isl_take isl_union_pw_multi_aff *upma2,
3651         int (*fn)(void **entry, void *user))
3652 {
3653         isl_space *space;
3654         struct isl_union_pw_multi_aff_bin_data data = { NULL, NULL, NULL, fn };
3655
3656         space = isl_union_pw_multi_aff_get_space(upma2);
3657         upma1 = isl_union_pw_multi_aff_align_params(upma1, space);
3658         space = isl_union_pw_multi_aff_get_space(upma1);
3659         upma2 = isl_union_pw_multi_aff_align_params(upma2, space);
3660
3661         if (!upma1 || !upma2)
3662                 goto error;
3663
3664         data.upma2 = upma2;
3665         data.res = isl_union_pw_multi_aff_alloc(isl_space_copy(upma1->dim),
3666                                        upma1->table.n);
3667         if (isl_hash_table_foreach(upma1->dim->ctx, &upma1->table,
3668                                    &bin_entry, &data) < 0)
3669                 goto error;
3670
3671         isl_union_pw_multi_aff_free(upma1);
3672         isl_union_pw_multi_aff_free(upma2);
3673         return data.res;
3674 error:
3675         isl_union_pw_multi_aff_free(upma1);
3676         isl_union_pw_multi_aff_free(upma2);
3677         isl_union_pw_multi_aff_free(data.res);
3678         return NULL;
3679 }
3680
3681 /* Given two isl_multi_affs A -> B and C -> D,
3682  * construct an isl_multi_aff (A * C) -> (B, D).
3683  */
3684 __isl_give isl_multi_aff *isl_multi_aff_flat_range_product(
3685         __isl_take isl_multi_aff *ma1, __isl_take isl_multi_aff *ma2)
3686 {
3687         int i, n1, n2;
3688         isl_aff *aff;
3689         isl_space *space;
3690         isl_multi_aff *res;
3691
3692         if (!ma1 || !ma2)
3693                 goto error;
3694
3695         space = isl_space_range_product(isl_multi_aff_get_space(ma1),
3696                                         isl_multi_aff_get_space(ma2));
3697         space = isl_space_flatten_range(space);
3698         res = isl_multi_aff_alloc(space);
3699
3700         n1 = isl_multi_aff_dim(ma1, isl_dim_out);
3701         n2 = isl_multi_aff_dim(ma2, isl_dim_out);
3702
3703         for (i = 0; i < n1; ++i) {
3704                 aff = isl_multi_aff_get_aff(ma1, i);
3705                 res = isl_multi_aff_set_aff(res, i, aff);
3706         }
3707
3708         for (i = 0; i < n2; ++i) {
3709                 aff = isl_multi_aff_get_aff(ma2, i);
3710                 res = isl_multi_aff_set_aff(res, n1 + i, aff);
3711         }
3712
3713         isl_multi_aff_free(ma1);
3714         isl_multi_aff_free(ma2);
3715         return res;
3716 error:
3717         isl_multi_aff_free(ma1);
3718         isl_multi_aff_free(ma2);
3719         return NULL;
3720 }
3721
3722 /* Given two aligned isl_pw_multi_affs A -> B and C -> D,
3723  * construct an isl_pw_multi_aff (A * C) -> (B, D).
3724  */
3725 static __isl_give isl_pw_multi_aff *pw_multi_aff_flat_range_product(
3726         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
3727 {
3728         isl_space *space;
3729
3730         space = isl_space_range_product(isl_pw_multi_aff_get_space(pma1),
3731                                         isl_pw_multi_aff_get_space(pma2));
3732         space = isl_space_flatten_range(space);
3733         return isl_pw_multi_aff_on_shared_domain_in(pma1, pma2, space,
3734                                             &isl_multi_aff_flat_range_product);
3735 }
3736
3737 /* Given two isl_pw_multi_affs A -> B and C -> D,
3738  * construct an isl_pw_multi_aff (A * C) -> (B, D).
3739  */
3740 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_flat_range_product(
3741         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
3742 {
3743         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
3744                                             &pw_multi_aff_flat_range_product);
3745 }
3746
3747 /* If data->pma and *entry have the same domain space, then compute
3748  * their flat range product and the result to data->res.
3749  */
3750 static int flat_range_product_entry(void **entry, void *user)
3751 {
3752         struct isl_union_pw_multi_aff_bin_data *data = user;
3753         isl_pw_multi_aff *pma2 = *entry;
3754
3755         if (!isl_space_tuple_match(data->pma->dim, isl_dim_in,
3756                                  pma2->dim, isl_dim_in))
3757                 return 0;
3758
3759         pma2 = isl_pw_multi_aff_flat_range_product(
3760                                         isl_pw_multi_aff_copy(data->pma),
3761                                         isl_pw_multi_aff_copy(pma2));
3762
3763         data->res = isl_union_pw_multi_aff_add_pw_multi_aff(data->res, pma2);
3764
3765         return 0;
3766 }
3767
3768 /* Given two isl_union_pw_multi_affs A -> B and C -> D,
3769  * construct an isl_union_pw_multi_aff (A * C) -> (B, D).
3770  */
3771 __isl_give isl_union_pw_multi_aff *isl_union_pw_multi_aff_flat_range_product(
3772         __isl_take isl_union_pw_multi_aff *upma1,
3773         __isl_take isl_union_pw_multi_aff *upma2)
3774 {
3775         return bin_op(upma1, upma2, &flat_range_product_entry);
3776 }
3777
3778 /* Replace the affine expressions at position "pos" in "pma" by "pa".
3779  * The parameters are assumed to have been aligned.
3780  *
3781  * The implementation essentially performs an isl_pw_*_on_shared_domain,
3782  * except that it works on two different isl_pw_* types.
3783  */
3784 static __isl_give isl_pw_multi_aff *pw_multi_aff_set_pw_aff(
3785         __isl_take isl_pw_multi_aff *pma, unsigned pos,
3786         __isl_take isl_pw_aff *pa)
3787 {
3788         int i, j, n;
3789         isl_pw_multi_aff *res = NULL;
3790
3791         if (!pma || !pa)
3792                 goto error;
3793
3794         if (!isl_space_tuple_match(pma->dim, isl_dim_in, pa->dim, isl_dim_in))
3795                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
3796                         "domains don't match", goto error);
3797         if (pos >= isl_pw_multi_aff_dim(pma, isl_dim_out))
3798                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
3799                         "index out of bounds", goto error);
3800
3801         n = pma->n * pa->n;
3802         res = isl_pw_multi_aff_alloc_size(isl_pw_multi_aff_get_space(pma), n);
3803
3804         for (i = 0; i < pma->n; ++i) {
3805                 for (j = 0; j < pa->n; ++j) {
3806                         isl_set *common;
3807                         isl_multi_aff *res_ij;
3808                         int empty;
3809
3810                         common = isl_set_intersect(isl_set_copy(pma->p[i].set),
3811                                                    isl_set_copy(pa->p[j].set));
3812                         empty = isl_set_plain_is_empty(common);
3813                         if (empty < 0 || empty) {
3814                                 isl_set_free(common);
3815                                 if (empty < 0)
3816                                         goto error;
3817                                 continue;
3818                         }
3819
3820                         res_ij = isl_multi_aff_set_aff(
3821                                         isl_multi_aff_copy(pma->p[i].maff), pos,
3822                                         isl_aff_copy(pa->p[j].aff));
3823                         res_ij = isl_multi_aff_gist(res_ij,
3824                                         isl_set_copy(common));
3825
3826                         res = isl_pw_multi_aff_add_piece(res, common, res_ij);
3827                 }
3828         }
3829
3830         isl_pw_multi_aff_free(pma);
3831         isl_pw_aff_free(pa);
3832         return res;
3833 error:
3834         isl_pw_multi_aff_free(pma);
3835         isl_pw_aff_free(pa);
3836         return isl_pw_multi_aff_free(res);
3837 }
3838
3839 /* Replace the affine expressions at position "pos" in "pma" by "pa".
3840  */
3841 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_set_pw_aff(
3842         __isl_take isl_pw_multi_aff *pma, unsigned pos,
3843         __isl_take isl_pw_aff *pa)
3844 {
3845         if (!pma || !pa)
3846                 goto error;
3847         if (isl_space_match(pma->dim, isl_dim_param, pa->dim, isl_dim_param))
3848                 return pw_multi_aff_set_pw_aff(pma, pos, pa);
3849         if (!isl_space_has_named_params(pma->dim) ||
3850             !isl_space_has_named_params(pa->dim))
3851                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
3852                         "unaligned unnamed parameters", goto error);
3853         pma = isl_pw_multi_aff_align_params(pma, isl_pw_aff_get_space(pa));
3854         pa = isl_pw_aff_align_params(pa, isl_pw_multi_aff_get_space(pma));
3855         return pw_multi_aff_set_pw_aff(pma, pos, pa);
3856 error:
3857         isl_pw_multi_aff_free(pma);
3858         isl_pw_aff_free(pa);
3859         return NULL;
3860 }