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         aff->v = isl_vec_cow(aff->v);
1280         if (!aff->ls || !aff->v)
1281                 goto error;
1282
1283         total = 1 + isl_space_dim(eq->dim, isl_dim_all);
1284         n_div = eq->n_div;
1285         for (i = 0; i < eq->n_eq; ++i) {
1286                 j = isl_seq_last_non_zero(eq->eq[i], total + n_div);
1287                 if (j < 0 || j == 0 || j >= total)
1288                         continue;
1289
1290                 isl_seq_elim(aff->v->el + 1, eq->eq[i], j, total,
1291                                 &aff->v->el[0]);
1292         }
1293
1294         isl_basic_set_free(eq);
1295         aff = isl_aff_normalize(aff);
1296         return aff;
1297 error:
1298         isl_basic_set_free(eq);
1299         isl_aff_free(aff);
1300         return NULL;
1301 }
1302
1303 /* Exploit the equalities in "eq" to simplify the affine expression
1304  * and the expressions of the integer divisions in the local space.
1305  */
1306 static __isl_give isl_aff *isl_aff_substitute_equalities(
1307         __isl_take isl_aff *aff, __isl_take isl_basic_set *eq)
1308 {
1309         int n_div;
1310
1311         if (!aff || !eq)
1312                 goto error;
1313         n_div = isl_local_space_dim(aff->ls, isl_dim_div);
1314         if (n_div > 0)
1315                 eq = isl_basic_set_add(eq, isl_dim_set, n_div);
1316         return isl_aff_substitute_equalities_lifted(aff, eq);
1317 error:
1318         isl_basic_set_free(eq);
1319         isl_aff_free(aff);
1320         return NULL;
1321 }
1322
1323 /* Look for equalities among the variables shared by context and aff
1324  * and the integer divisions of aff, if any.
1325  * The equalities are then used to eliminate coefficients and/or integer
1326  * divisions from aff.
1327  */
1328 __isl_give isl_aff *isl_aff_gist(__isl_take isl_aff *aff,
1329         __isl_take isl_set *context)
1330 {
1331         isl_basic_set *hull;
1332         int n_div;
1333
1334         if (!aff)
1335                 goto error;
1336         n_div = isl_local_space_dim(aff->ls, isl_dim_div);
1337         if (n_div > 0) {
1338                 isl_basic_set *bset;
1339                 isl_local_space *ls;
1340                 context = isl_set_add_dims(context, isl_dim_set, n_div);
1341                 ls = isl_aff_get_domain_local_space(aff);
1342                 bset = isl_basic_set_from_local_space(ls);
1343                 bset = isl_basic_set_lift(bset);
1344                 bset = isl_basic_set_flatten(bset);
1345                 context = isl_set_intersect(context,
1346                                             isl_set_from_basic_set(bset));
1347         }
1348
1349         hull = isl_set_affine_hull(context);
1350         return isl_aff_substitute_equalities_lifted(aff, hull);
1351 error:
1352         isl_aff_free(aff);
1353         isl_set_free(context);
1354         return NULL;
1355 }
1356
1357 __isl_give isl_aff *isl_aff_gist_params(__isl_take isl_aff *aff,
1358         __isl_take isl_set *context)
1359 {
1360         isl_set *dom_context = isl_set_universe(isl_aff_get_domain_space(aff));
1361         dom_context = isl_set_intersect_params(dom_context, context);
1362         return isl_aff_gist(aff, dom_context);
1363 }
1364
1365 /* Return a basic set containing those elements in the space
1366  * of aff where it is non-negative.
1367  * If "rational" is set, then return a rational basic set.
1368  */
1369 static __isl_give isl_basic_set *aff_nonneg_basic_set(
1370         __isl_take isl_aff *aff, int rational)
1371 {
1372         isl_constraint *ineq;
1373         isl_basic_set *bset;
1374
1375         ineq = isl_inequality_from_aff(aff);
1376
1377         bset = isl_basic_set_from_constraint(ineq);
1378         if (rational)
1379                 bset = isl_basic_set_set_rational(bset);
1380         bset = isl_basic_set_simplify(bset);
1381         return bset;
1382 }
1383
1384 /* Return a basic set containing those elements in the space
1385  * of aff where it is non-negative.
1386  */
1387 __isl_give isl_basic_set *isl_aff_nonneg_basic_set(__isl_take isl_aff *aff)
1388 {
1389         return aff_nonneg_basic_set(aff, 0);
1390 }
1391
1392 /* Return a basic set containing those elements in the domain space
1393  * of aff where it is negative.
1394  */
1395 __isl_give isl_basic_set *isl_aff_neg_basic_set(__isl_take isl_aff *aff)
1396 {
1397         aff = isl_aff_neg(aff);
1398         aff = isl_aff_add_constant_num_si(aff, -1);
1399         return isl_aff_nonneg_basic_set(aff);
1400 }
1401
1402 /* Return a basic set containing those elements in the space
1403  * of aff where it is zero.
1404  * If "rational" is set, then return a rational basic set.
1405  */
1406 static __isl_give isl_basic_set *aff_zero_basic_set(__isl_take isl_aff *aff,
1407         int rational)
1408 {
1409         isl_constraint *ineq;
1410         isl_basic_set *bset;
1411
1412         ineq = isl_equality_from_aff(aff);
1413
1414         bset = isl_basic_set_from_constraint(ineq);
1415         if (rational)
1416                 bset = isl_basic_set_set_rational(bset);
1417         bset = isl_basic_set_simplify(bset);
1418         return bset;
1419 }
1420
1421 /* Return a basic set containing those elements in the space
1422  * of aff where it is zero.
1423  */
1424 __isl_give isl_basic_set *isl_aff_zero_basic_set(__isl_take isl_aff *aff)
1425 {
1426         return aff_zero_basic_set(aff, 0);
1427 }
1428
1429 /* Return a basic set containing those elements in the shared space
1430  * of aff1 and aff2 where aff1 is greater than or equal to aff2.
1431  */
1432 __isl_give isl_basic_set *isl_aff_ge_basic_set(__isl_take isl_aff *aff1,
1433         __isl_take isl_aff *aff2)
1434 {
1435         aff1 = isl_aff_sub(aff1, aff2);
1436
1437         return isl_aff_nonneg_basic_set(aff1);
1438 }
1439
1440 /* Return a basic set containing those elements in the shared space
1441  * of aff1 and aff2 where aff1 is smaller than or equal to aff2.
1442  */
1443 __isl_give isl_basic_set *isl_aff_le_basic_set(__isl_take isl_aff *aff1,
1444         __isl_take isl_aff *aff2)
1445 {
1446         return isl_aff_ge_basic_set(aff2, aff1);
1447 }
1448
1449 __isl_give isl_aff *isl_aff_add_on_domain(__isl_keep isl_set *dom,
1450         __isl_take isl_aff *aff1, __isl_take isl_aff *aff2)
1451 {
1452         aff1 = isl_aff_add(aff1, aff2);
1453         aff1 = isl_aff_gist(aff1, isl_set_copy(dom));
1454         return aff1;
1455 }
1456
1457 int isl_aff_is_empty(__isl_keep isl_aff *aff)
1458 {
1459         if (!aff)
1460                 return -1;
1461
1462         return 0;
1463 }
1464
1465 /* Check whether the given affine expression has non-zero coefficient
1466  * for any dimension in the given range or if any of these dimensions
1467  * appear with non-zero coefficients in any of the integer divisions
1468  * involved in the affine expression.
1469  */
1470 int isl_aff_involves_dims(__isl_keep isl_aff *aff,
1471         enum isl_dim_type type, unsigned first, unsigned n)
1472 {
1473         int i;
1474         isl_ctx *ctx;
1475         int *active = NULL;
1476         int involves = 0;
1477
1478         if (!aff)
1479                 return -1;
1480         if (n == 0)
1481                 return 0;
1482
1483         ctx = isl_aff_get_ctx(aff);
1484         if (first + n > isl_aff_dim(aff, type))
1485                 isl_die(ctx, isl_error_invalid,
1486                         "range out of bounds", return -1);
1487
1488         active = isl_local_space_get_active(aff->ls, aff->v->el + 2);
1489         if (!active)
1490                 goto error;
1491
1492         first += isl_local_space_offset(aff->ls, type) - 1;
1493         for (i = 0; i < n; ++i)
1494                 if (active[first + i]) {
1495                         involves = 1;
1496                         break;
1497                 }
1498
1499         free(active);
1500
1501         return involves;
1502 error:
1503         free(active);
1504         return -1;
1505 }
1506
1507 __isl_give isl_aff *isl_aff_drop_dims(__isl_take isl_aff *aff,
1508         enum isl_dim_type type, unsigned first, unsigned n)
1509 {
1510         isl_ctx *ctx;
1511
1512         if (!aff)
1513                 return NULL;
1514         if (type == isl_dim_out)
1515                 isl_die(aff->v->ctx, isl_error_invalid,
1516                         "cannot drop output/set dimension",
1517                         return isl_aff_free(aff));
1518         if (type == isl_dim_in)
1519                 type = isl_dim_set;
1520         if (n == 0 && !isl_local_space_is_named_or_nested(aff->ls, type))
1521                 return aff;
1522
1523         ctx = isl_aff_get_ctx(aff);
1524         if (first + n > isl_local_space_dim(aff->ls, type))
1525                 isl_die(ctx, isl_error_invalid, "range out of bounds",
1526                         return isl_aff_free(aff));
1527
1528         aff = isl_aff_cow(aff);
1529         if (!aff)
1530                 return NULL;
1531
1532         aff->ls = isl_local_space_drop_dims(aff->ls, type, first, n);
1533         if (!aff->ls)
1534                 return isl_aff_free(aff);
1535
1536         first += 1 + isl_local_space_offset(aff->ls, type);
1537         aff->v = isl_vec_drop_els(aff->v, first, n);
1538         if (!aff->v)
1539                 return isl_aff_free(aff);
1540
1541         return aff;
1542 }
1543
1544 /* Project the domain of the affine expression onto its parameter space.
1545  * The affine expression may not involve any of the domain dimensions.
1546  */
1547 __isl_give isl_aff *isl_aff_project_domain_on_params(__isl_take isl_aff *aff)
1548 {
1549         isl_space *space;
1550         unsigned n;
1551         int involves;
1552
1553         n = isl_aff_dim(aff, isl_dim_in);
1554         involves = isl_aff_involves_dims(aff, isl_dim_in, 0, n);
1555         if (involves < 0)
1556                 return isl_aff_free(aff);
1557         if (involves)
1558                 isl_die(isl_aff_get_ctx(aff), isl_error_invalid,
1559                     "affine expression involves some of the domain dimensions",
1560                     return isl_aff_free(aff));
1561         aff = isl_aff_drop_dims(aff, isl_dim_in, 0, n);
1562         space = isl_aff_get_domain_space(aff);
1563         space = isl_space_params(space);
1564         aff = isl_aff_reset_domain_space(aff, space);
1565         return aff;
1566 }
1567
1568 __isl_give isl_aff *isl_aff_insert_dims(__isl_take isl_aff *aff,
1569         enum isl_dim_type type, unsigned first, unsigned n)
1570 {
1571         isl_ctx *ctx;
1572
1573         if (!aff)
1574                 return NULL;
1575         if (type == isl_dim_out)
1576                 isl_die(aff->v->ctx, isl_error_invalid,
1577                         "cannot insert output/set dimensions",
1578                         return isl_aff_free(aff));
1579         if (type == isl_dim_in)
1580                 type = isl_dim_set;
1581         if (n == 0 && !isl_local_space_is_named_or_nested(aff->ls, type))
1582                 return aff;
1583
1584         ctx = isl_aff_get_ctx(aff);
1585         if (first > isl_local_space_dim(aff->ls, type))
1586                 isl_die(ctx, isl_error_invalid, "position out of bounds",
1587                         return isl_aff_free(aff));
1588
1589         aff = isl_aff_cow(aff);
1590         if (!aff)
1591                 return NULL;
1592
1593         aff->ls = isl_local_space_insert_dims(aff->ls, type, first, n);
1594         if (!aff->ls)
1595                 return isl_aff_free(aff);
1596
1597         first += 1 + isl_local_space_offset(aff->ls, type);
1598         aff->v = isl_vec_insert_zero_els(aff->v, first, n);
1599         if (!aff->v)
1600                 return isl_aff_free(aff);
1601
1602         return aff;
1603 }
1604
1605 __isl_give isl_aff *isl_aff_add_dims(__isl_take isl_aff *aff,
1606         enum isl_dim_type type, unsigned n)
1607 {
1608         unsigned pos;
1609
1610         pos = isl_aff_dim(aff, type);
1611
1612         return isl_aff_insert_dims(aff, type, pos, n);
1613 }
1614
1615 __isl_give isl_pw_aff *isl_pw_aff_add_dims(__isl_take isl_pw_aff *pwaff,
1616         enum isl_dim_type type, unsigned n)
1617 {
1618         unsigned pos;
1619
1620         pos = isl_pw_aff_dim(pwaff, type);
1621
1622         return isl_pw_aff_insert_dims(pwaff, type, pos, n);
1623 }
1624
1625 __isl_give isl_pw_aff *isl_pw_aff_from_aff(__isl_take isl_aff *aff)
1626 {
1627         isl_set *dom = isl_set_universe(isl_aff_get_domain_space(aff));
1628         return isl_pw_aff_alloc(dom, aff);
1629 }
1630
1631 #undef PW
1632 #define PW isl_pw_aff
1633 #undef EL
1634 #define EL isl_aff
1635 #undef EL_IS_ZERO
1636 #define EL_IS_ZERO is_empty
1637 #undef ZERO
1638 #define ZERO empty
1639 #undef IS_ZERO
1640 #define IS_ZERO is_empty
1641 #undef FIELD
1642 #define FIELD aff
1643 #undef DEFAULT_IS_ZERO
1644 #define DEFAULT_IS_ZERO 0
1645
1646 #define NO_EVAL
1647 #define NO_OPT
1648 #define NO_MOVE_DIMS
1649 #define NO_LIFT
1650 #define NO_MORPH
1651
1652 #include <isl_pw_templ.c>
1653
1654 static __isl_give isl_set *align_params_pw_pw_set_and(
1655         __isl_take isl_pw_aff *pwaff1, __isl_take isl_pw_aff *pwaff2,
1656         __isl_give isl_set *(*fn)(__isl_take isl_pw_aff *pwaff1,
1657                                     __isl_take isl_pw_aff *pwaff2))
1658 {
1659         if (!pwaff1 || !pwaff2)
1660                 goto error;
1661         if (isl_space_match(pwaff1->dim, isl_dim_param,
1662                           pwaff2->dim, isl_dim_param))
1663                 return fn(pwaff1, pwaff2);
1664         if (!isl_space_has_named_params(pwaff1->dim) ||
1665             !isl_space_has_named_params(pwaff2->dim))
1666                 isl_die(isl_pw_aff_get_ctx(pwaff1), isl_error_invalid,
1667                         "unaligned unnamed parameters", goto error);
1668         pwaff1 = isl_pw_aff_align_params(pwaff1, isl_pw_aff_get_space(pwaff2));
1669         pwaff2 = isl_pw_aff_align_params(pwaff2, isl_pw_aff_get_space(pwaff1));
1670         return fn(pwaff1, pwaff2);
1671 error:
1672         isl_pw_aff_free(pwaff1);
1673         isl_pw_aff_free(pwaff2);
1674         return NULL;
1675 }
1676
1677 /* Compute a piecewise quasi-affine expression with a domain that
1678  * is the union of those of pwaff1 and pwaff2 and such that on each
1679  * cell, the quasi-affine expression is the better (according to cmp)
1680  * of those of pwaff1 and pwaff2.  If only one of pwaff1 or pwaff2
1681  * is defined on a given cell, then the associated expression
1682  * is the defined one.
1683  */
1684 static __isl_give isl_pw_aff *pw_aff_union_opt(__isl_take isl_pw_aff *pwaff1,
1685         __isl_take isl_pw_aff *pwaff2,
1686         __isl_give isl_basic_set *(*cmp)(__isl_take isl_aff *aff1,
1687                                         __isl_take isl_aff *aff2))
1688 {
1689         int i, j, n;
1690         isl_pw_aff *res;
1691         isl_ctx *ctx;
1692         isl_set *set;
1693
1694         if (!pwaff1 || !pwaff2)
1695                 goto error;
1696
1697         ctx = isl_space_get_ctx(pwaff1->dim);
1698         if (!isl_space_is_equal(pwaff1->dim, pwaff2->dim))
1699                 isl_die(ctx, isl_error_invalid,
1700                         "arguments should live in same space", goto error);
1701
1702         if (isl_pw_aff_is_empty(pwaff1)) {
1703                 isl_pw_aff_free(pwaff1);
1704                 return pwaff2;
1705         }
1706
1707         if (isl_pw_aff_is_empty(pwaff2)) {
1708                 isl_pw_aff_free(pwaff2);
1709                 return pwaff1;
1710         }
1711
1712         n = 2 * (pwaff1->n + 1) * (pwaff2->n + 1);
1713         res = isl_pw_aff_alloc_size(isl_space_copy(pwaff1->dim), n);
1714
1715         for (i = 0; i < pwaff1->n; ++i) {
1716                 set = isl_set_copy(pwaff1->p[i].set);
1717                 for (j = 0; j < pwaff2->n; ++j) {
1718                         struct isl_set *common;
1719                         isl_set *better;
1720
1721                         common = isl_set_intersect(
1722                                         isl_set_copy(pwaff1->p[i].set),
1723                                         isl_set_copy(pwaff2->p[j].set));
1724                         better = isl_set_from_basic_set(cmp(
1725                                         isl_aff_copy(pwaff2->p[j].aff),
1726                                         isl_aff_copy(pwaff1->p[i].aff)));
1727                         better = isl_set_intersect(common, better);
1728                         if (isl_set_plain_is_empty(better)) {
1729                                 isl_set_free(better);
1730                                 continue;
1731                         }
1732                         set = isl_set_subtract(set, isl_set_copy(better));
1733
1734                         res = isl_pw_aff_add_piece(res, better,
1735                                                 isl_aff_copy(pwaff2->p[j].aff));
1736                 }
1737                 res = isl_pw_aff_add_piece(res, set,
1738                                                 isl_aff_copy(pwaff1->p[i].aff));
1739         }
1740
1741         for (j = 0; j < pwaff2->n; ++j) {
1742                 set = isl_set_copy(pwaff2->p[j].set);
1743                 for (i = 0; i < pwaff1->n; ++i)
1744                         set = isl_set_subtract(set,
1745                                         isl_set_copy(pwaff1->p[i].set));
1746                 res = isl_pw_aff_add_piece(res, set,
1747                                                 isl_aff_copy(pwaff2->p[j].aff));
1748         }
1749
1750         isl_pw_aff_free(pwaff1);
1751         isl_pw_aff_free(pwaff2);
1752
1753         return res;
1754 error:
1755         isl_pw_aff_free(pwaff1);
1756         isl_pw_aff_free(pwaff2);
1757         return NULL;
1758 }
1759
1760 /* Compute a piecewise quasi-affine expression with a domain that
1761  * is the union of those of pwaff1 and pwaff2 and such that on each
1762  * cell, the quasi-affine expression is the maximum of those of pwaff1
1763  * and pwaff2.  If only one of pwaff1 or pwaff2 is defined on a given
1764  * cell, then the associated expression is the defined one.
1765  */
1766 static __isl_give isl_pw_aff *pw_aff_union_max(__isl_take isl_pw_aff *pwaff1,
1767         __isl_take isl_pw_aff *pwaff2)
1768 {
1769         return pw_aff_union_opt(pwaff1, pwaff2, &isl_aff_ge_basic_set);
1770 }
1771
1772 __isl_give isl_pw_aff *isl_pw_aff_union_max(__isl_take isl_pw_aff *pwaff1,
1773         __isl_take isl_pw_aff *pwaff2)
1774 {
1775         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2,
1776                                                         &pw_aff_union_max);
1777 }
1778
1779 /* Compute a piecewise quasi-affine expression with a domain that
1780  * is the union of those of pwaff1 and pwaff2 and such that on each
1781  * cell, the quasi-affine expression is the minimum of those of pwaff1
1782  * and pwaff2.  If only one of pwaff1 or pwaff2 is defined on a given
1783  * cell, then the associated expression is the defined one.
1784  */
1785 static __isl_give isl_pw_aff *pw_aff_union_min(__isl_take isl_pw_aff *pwaff1,
1786         __isl_take isl_pw_aff *pwaff2)
1787 {
1788         return pw_aff_union_opt(pwaff1, pwaff2, &isl_aff_le_basic_set);
1789 }
1790
1791 __isl_give isl_pw_aff *isl_pw_aff_union_min(__isl_take isl_pw_aff *pwaff1,
1792         __isl_take isl_pw_aff *pwaff2)
1793 {
1794         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2,
1795                                                         &pw_aff_union_min);
1796 }
1797
1798 __isl_give isl_pw_aff *isl_pw_aff_union_opt(__isl_take isl_pw_aff *pwaff1,
1799         __isl_take isl_pw_aff *pwaff2, int max)
1800 {
1801         if (max)
1802                 return isl_pw_aff_union_max(pwaff1, pwaff2);
1803         else
1804                 return isl_pw_aff_union_min(pwaff1, pwaff2);
1805 }
1806
1807 /* Construct a map with as domain the domain of pwaff and
1808  * one-dimensional range corresponding to the affine expressions.
1809  */
1810 static __isl_give isl_map *map_from_pw_aff(__isl_take isl_pw_aff *pwaff)
1811 {
1812         int i;
1813         isl_space *dim;
1814         isl_map *map;
1815
1816         if (!pwaff)
1817                 return NULL;
1818
1819         dim = isl_pw_aff_get_space(pwaff);
1820         map = isl_map_empty(dim);
1821
1822         for (i = 0; i < pwaff->n; ++i) {
1823                 isl_basic_map *bmap;
1824                 isl_map *map_i;
1825
1826                 bmap = isl_basic_map_from_aff(isl_aff_copy(pwaff->p[i].aff));
1827                 map_i = isl_map_from_basic_map(bmap);
1828                 map_i = isl_map_intersect_domain(map_i,
1829                                                 isl_set_copy(pwaff->p[i].set));
1830                 map = isl_map_union_disjoint(map, map_i);
1831         }
1832
1833         isl_pw_aff_free(pwaff);
1834
1835         return map;
1836 }
1837
1838 /* Construct a map with as domain the domain of pwaff and
1839  * one-dimensional range corresponding to the affine expressions.
1840  */
1841 __isl_give isl_map *isl_map_from_pw_aff(__isl_take isl_pw_aff *pwaff)
1842 {
1843         if (!pwaff)
1844                 return NULL;
1845         if (isl_space_is_set(pwaff->dim))
1846                 isl_die(isl_pw_aff_get_ctx(pwaff), isl_error_invalid,
1847                         "space of input is not a map",
1848                         return isl_pw_aff_free(pwaff));
1849         return map_from_pw_aff(pwaff);
1850 }
1851
1852 /* Construct a one-dimensional set with as parameter domain
1853  * the domain of pwaff and the single set dimension
1854  * corresponding to the affine expressions.
1855  */
1856 __isl_give isl_set *isl_set_from_pw_aff(__isl_take isl_pw_aff *pwaff)
1857 {
1858         if (!pwaff)
1859                 return NULL;
1860         if (!isl_space_is_set(pwaff->dim))
1861                 isl_die(isl_pw_aff_get_ctx(pwaff), isl_error_invalid,
1862                         "space of input is not a set",
1863                         return isl_pw_aff_free(pwaff));
1864         return map_from_pw_aff(pwaff);
1865 }
1866
1867 /* Return a set containing those elements in the domain
1868  * of pwaff where it is non-negative.
1869  */
1870 __isl_give isl_set *isl_pw_aff_nonneg_set(__isl_take isl_pw_aff *pwaff)
1871 {
1872         int i;
1873         isl_set *set;
1874
1875         if (!pwaff)
1876                 return NULL;
1877
1878         set = isl_set_empty(isl_pw_aff_get_domain_space(pwaff));
1879
1880         for (i = 0; i < pwaff->n; ++i) {
1881                 isl_basic_set *bset;
1882                 isl_set *set_i;
1883                 int rational;
1884
1885                 rational = isl_set_has_rational(pwaff->p[i].set);
1886                 bset = aff_nonneg_basic_set(isl_aff_copy(pwaff->p[i].aff),
1887                                                 rational);
1888                 set_i = isl_set_from_basic_set(bset);
1889                 set_i = isl_set_intersect(set_i, isl_set_copy(pwaff->p[i].set));
1890                 set = isl_set_union_disjoint(set, set_i);
1891         }
1892
1893         isl_pw_aff_free(pwaff);
1894
1895         return set;
1896 }
1897
1898 /* Return a set containing those elements in the domain
1899  * of pwaff where it is zero (if complement is 0) or not zero
1900  * (if complement is 1).
1901  */
1902 static __isl_give isl_set *pw_aff_zero_set(__isl_take isl_pw_aff *pwaff,
1903         int complement)
1904 {
1905         int i;
1906         isl_set *set;
1907
1908         if (!pwaff)
1909                 return NULL;
1910
1911         set = isl_set_empty(isl_pw_aff_get_domain_space(pwaff));
1912
1913         for (i = 0; i < pwaff->n; ++i) {
1914                 isl_basic_set *bset;
1915                 isl_set *set_i, *zero;
1916                 int rational;
1917
1918                 rational = isl_set_has_rational(pwaff->p[i].set);
1919                 bset = aff_zero_basic_set(isl_aff_copy(pwaff->p[i].aff),
1920                                                 rational);
1921                 zero = isl_set_from_basic_set(bset);
1922                 set_i = isl_set_copy(pwaff->p[i].set);
1923                 if (complement)
1924                         set_i = isl_set_subtract(set_i, zero);
1925                 else
1926                         set_i = isl_set_intersect(set_i, zero);
1927                 set = isl_set_union_disjoint(set, set_i);
1928         }
1929
1930         isl_pw_aff_free(pwaff);
1931
1932         return set;
1933 }
1934
1935 /* Return a set containing those elements in the domain
1936  * of pwaff where it is zero.
1937  */
1938 __isl_give isl_set *isl_pw_aff_zero_set(__isl_take isl_pw_aff *pwaff)
1939 {
1940         return pw_aff_zero_set(pwaff, 0);
1941 }
1942
1943 /* Return a set containing those elements in the domain
1944  * of pwaff where it is not zero.
1945  */
1946 __isl_give isl_set *isl_pw_aff_non_zero_set(__isl_take isl_pw_aff *pwaff)
1947 {
1948         return pw_aff_zero_set(pwaff, 1);
1949 }
1950
1951 /* Return a set containing those elements in the shared domain
1952  * of pwaff1 and pwaff2 where pwaff1 is greater than (or equal) to pwaff2.
1953  *
1954  * We compute the difference on the shared domain and then construct
1955  * the set of values where this difference is non-negative.
1956  * If strict is set, we first subtract 1 from the difference.
1957  * If equal is set, we only return the elements where pwaff1 and pwaff2
1958  * are equal.
1959  */
1960 static __isl_give isl_set *pw_aff_gte_set(__isl_take isl_pw_aff *pwaff1,
1961         __isl_take isl_pw_aff *pwaff2, int strict, int equal)
1962 {
1963         isl_set *set1, *set2;
1964
1965         set1 = isl_pw_aff_domain(isl_pw_aff_copy(pwaff1));
1966         set2 = isl_pw_aff_domain(isl_pw_aff_copy(pwaff2));
1967         set1 = isl_set_intersect(set1, set2);
1968         pwaff1 = isl_pw_aff_intersect_domain(pwaff1, isl_set_copy(set1));
1969         pwaff2 = isl_pw_aff_intersect_domain(pwaff2, isl_set_copy(set1));
1970         pwaff1 = isl_pw_aff_add(pwaff1, isl_pw_aff_neg(pwaff2));
1971
1972         if (strict) {
1973                 isl_space *dim = isl_set_get_space(set1);
1974                 isl_aff *aff;
1975                 aff = isl_aff_zero_on_domain(isl_local_space_from_space(dim));
1976                 aff = isl_aff_add_constant_si(aff, -1);
1977                 pwaff1 = isl_pw_aff_add(pwaff1, isl_pw_aff_alloc(set1, aff));
1978         } else
1979                 isl_set_free(set1);
1980
1981         if (equal)
1982                 return isl_pw_aff_zero_set(pwaff1);
1983         return isl_pw_aff_nonneg_set(pwaff1);
1984 }
1985
1986 /* Return a set containing those elements in the shared domain
1987  * of pwaff1 and pwaff2 where pwaff1 is equal to pwaff2.
1988  */
1989 static __isl_give isl_set *pw_aff_eq_set(__isl_take isl_pw_aff *pwaff1,
1990         __isl_take isl_pw_aff *pwaff2)
1991 {
1992         return pw_aff_gte_set(pwaff1, pwaff2, 0, 1);
1993 }
1994
1995 __isl_give isl_set *isl_pw_aff_eq_set(__isl_take isl_pw_aff *pwaff1,
1996         __isl_take isl_pw_aff *pwaff2)
1997 {
1998         return align_params_pw_pw_set_and(pwaff1, pwaff2, &pw_aff_eq_set);
1999 }
2000
2001 /* Return a set containing those elements in the shared domain
2002  * of pwaff1 and pwaff2 where pwaff1 is greater than or equal to pwaff2.
2003  */
2004 static __isl_give isl_set *pw_aff_ge_set(__isl_take isl_pw_aff *pwaff1,
2005         __isl_take isl_pw_aff *pwaff2)
2006 {
2007         return pw_aff_gte_set(pwaff1, pwaff2, 0, 0);
2008 }
2009
2010 __isl_give isl_set *isl_pw_aff_ge_set(__isl_take isl_pw_aff *pwaff1,
2011         __isl_take isl_pw_aff *pwaff2)
2012 {
2013         return align_params_pw_pw_set_and(pwaff1, pwaff2, &pw_aff_ge_set);
2014 }
2015
2016 /* Return a set containing those elements in the shared domain
2017  * of pwaff1 and pwaff2 where pwaff1 is strictly greater than pwaff2.
2018  */
2019 static __isl_give isl_set *pw_aff_gt_set(__isl_take isl_pw_aff *pwaff1,
2020         __isl_take isl_pw_aff *pwaff2)
2021 {
2022         return pw_aff_gte_set(pwaff1, pwaff2, 1, 0);
2023 }
2024
2025 __isl_give isl_set *isl_pw_aff_gt_set(__isl_take isl_pw_aff *pwaff1,
2026         __isl_take isl_pw_aff *pwaff2)
2027 {
2028         return align_params_pw_pw_set_and(pwaff1, pwaff2, &pw_aff_gt_set);
2029 }
2030
2031 __isl_give isl_set *isl_pw_aff_le_set(__isl_take isl_pw_aff *pwaff1,
2032         __isl_take isl_pw_aff *pwaff2)
2033 {
2034         return isl_pw_aff_ge_set(pwaff2, pwaff1);
2035 }
2036
2037 __isl_give isl_set *isl_pw_aff_lt_set(__isl_take isl_pw_aff *pwaff1,
2038         __isl_take isl_pw_aff *pwaff2)
2039 {
2040         return isl_pw_aff_gt_set(pwaff2, pwaff1);
2041 }
2042
2043 /* Return a set containing those elements in the shared domain
2044  * of the elements of list1 and list2 where each element in list1
2045  * has the relation specified by "fn" with each element in list2.
2046  */
2047 static __isl_give isl_set *pw_aff_list_set(__isl_take isl_pw_aff_list *list1,
2048         __isl_take isl_pw_aff_list *list2,
2049         __isl_give isl_set *(*fn)(__isl_take isl_pw_aff *pwaff1,
2050                                     __isl_take isl_pw_aff *pwaff2))
2051 {
2052         int i, j;
2053         isl_ctx *ctx;
2054         isl_set *set;
2055
2056         if (!list1 || !list2)
2057                 goto error;
2058
2059         ctx = isl_pw_aff_list_get_ctx(list1);
2060         if (list1->n < 1 || list2->n < 1)
2061                 isl_die(ctx, isl_error_invalid,
2062                         "list should contain at least one element", goto error);
2063
2064         set = isl_set_universe(isl_pw_aff_get_domain_space(list1->p[0]));
2065         for (i = 0; i < list1->n; ++i)
2066                 for (j = 0; j < list2->n; ++j) {
2067                         isl_set *set_ij;
2068
2069                         set_ij = fn(isl_pw_aff_copy(list1->p[i]),
2070                                     isl_pw_aff_copy(list2->p[j]));
2071                         set = isl_set_intersect(set, set_ij);
2072                 }
2073
2074         isl_pw_aff_list_free(list1);
2075         isl_pw_aff_list_free(list2);
2076         return set;
2077 error:
2078         isl_pw_aff_list_free(list1);
2079         isl_pw_aff_list_free(list2);
2080         return NULL;
2081 }
2082
2083 /* Return a set containing those elements in the shared domain
2084  * of the elements of list1 and list2 where each element in list1
2085  * is equal to each element in list2.
2086  */
2087 __isl_give isl_set *isl_pw_aff_list_eq_set(__isl_take isl_pw_aff_list *list1,
2088         __isl_take isl_pw_aff_list *list2)
2089 {
2090         return pw_aff_list_set(list1, list2, &isl_pw_aff_eq_set);
2091 }
2092
2093 __isl_give isl_set *isl_pw_aff_list_ne_set(__isl_take isl_pw_aff_list *list1,
2094         __isl_take isl_pw_aff_list *list2)
2095 {
2096         return pw_aff_list_set(list1, list2, &isl_pw_aff_ne_set);
2097 }
2098
2099 /* Return a set containing those elements in the shared domain
2100  * of the elements of list1 and list2 where each element in list1
2101  * is less than or equal to each element in list2.
2102  */
2103 __isl_give isl_set *isl_pw_aff_list_le_set(__isl_take isl_pw_aff_list *list1,
2104         __isl_take isl_pw_aff_list *list2)
2105 {
2106         return pw_aff_list_set(list1, list2, &isl_pw_aff_le_set);
2107 }
2108
2109 __isl_give isl_set *isl_pw_aff_list_lt_set(__isl_take isl_pw_aff_list *list1,
2110         __isl_take isl_pw_aff_list *list2)
2111 {
2112         return pw_aff_list_set(list1, list2, &isl_pw_aff_lt_set);
2113 }
2114
2115 __isl_give isl_set *isl_pw_aff_list_ge_set(__isl_take isl_pw_aff_list *list1,
2116         __isl_take isl_pw_aff_list *list2)
2117 {
2118         return pw_aff_list_set(list1, list2, &isl_pw_aff_ge_set);
2119 }
2120
2121 __isl_give isl_set *isl_pw_aff_list_gt_set(__isl_take isl_pw_aff_list *list1,
2122         __isl_take isl_pw_aff_list *list2)
2123 {
2124         return pw_aff_list_set(list1, list2, &isl_pw_aff_gt_set);
2125 }
2126
2127
2128 /* Return a set containing those elements in the shared domain
2129  * of pwaff1 and pwaff2 where pwaff1 is not equal to pwaff2.
2130  */
2131 static __isl_give isl_set *pw_aff_ne_set(__isl_take isl_pw_aff *pwaff1,
2132         __isl_take isl_pw_aff *pwaff2)
2133 {
2134         isl_set *set_lt, *set_gt;
2135
2136         set_lt = isl_pw_aff_lt_set(isl_pw_aff_copy(pwaff1),
2137                                    isl_pw_aff_copy(pwaff2));
2138         set_gt = isl_pw_aff_gt_set(pwaff1, pwaff2);
2139         return isl_set_union_disjoint(set_lt, set_gt);
2140 }
2141
2142 __isl_give isl_set *isl_pw_aff_ne_set(__isl_take isl_pw_aff *pwaff1,
2143         __isl_take isl_pw_aff *pwaff2)
2144 {
2145         return align_params_pw_pw_set_and(pwaff1, pwaff2, &pw_aff_ne_set);
2146 }
2147
2148 __isl_give isl_pw_aff *isl_pw_aff_scale_down(__isl_take isl_pw_aff *pwaff,
2149         isl_int v)
2150 {
2151         int i;
2152
2153         if (isl_int_is_one(v))
2154                 return pwaff;
2155         if (!isl_int_is_pos(v))
2156                 isl_die(isl_pw_aff_get_ctx(pwaff), isl_error_invalid,
2157                         "factor needs to be positive",
2158                         return isl_pw_aff_free(pwaff));
2159         pwaff = isl_pw_aff_cow(pwaff);
2160         if (!pwaff)
2161                 return NULL;
2162         if (pwaff->n == 0)
2163                 return pwaff;
2164
2165         for (i = 0; i < pwaff->n; ++i) {
2166                 pwaff->p[i].aff = isl_aff_scale_down(pwaff->p[i].aff, v);
2167                 if (!pwaff->p[i].aff)
2168                         return isl_pw_aff_free(pwaff);
2169         }
2170
2171         return pwaff;
2172 }
2173
2174 __isl_give isl_pw_aff *isl_pw_aff_floor(__isl_take isl_pw_aff *pwaff)
2175 {
2176         int i;
2177
2178         pwaff = isl_pw_aff_cow(pwaff);
2179         if (!pwaff)
2180                 return NULL;
2181         if (pwaff->n == 0)
2182                 return pwaff;
2183
2184         for (i = 0; i < pwaff->n; ++i) {
2185                 pwaff->p[i].aff = isl_aff_floor(pwaff->p[i].aff);
2186                 if (!pwaff->p[i].aff)
2187                         return isl_pw_aff_free(pwaff);
2188         }
2189
2190         return pwaff;
2191 }
2192
2193 __isl_give isl_pw_aff *isl_pw_aff_ceil(__isl_take isl_pw_aff *pwaff)
2194 {
2195         int i;
2196
2197         pwaff = isl_pw_aff_cow(pwaff);
2198         if (!pwaff)
2199                 return NULL;
2200         if (pwaff->n == 0)
2201                 return pwaff;
2202
2203         for (i = 0; i < pwaff->n; ++i) {
2204                 pwaff->p[i].aff = isl_aff_ceil(pwaff->p[i].aff);
2205                 if (!pwaff->p[i].aff)
2206                         return isl_pw_aff_free(pwaff);
2207         }
2208
2209         return pwaff;
2210 }
2211
2212 /* Assuming that "cond1" and "cond2" are disjoint,
2213  * return an affine expression that is equal to pwaff1 on cond1
2214  * and to pwaff2 on cond2.
2215  */
2216 static __isl_give isl_pw_aff *isl_pw_aff_select(
2217         __isl_take isl_set *cond1, __isl_take isl_pw_aff *pwaff1,
2218         __isl_take isl_set *cond2, __isl_take isl_pw_aff *pwaff2)
2219 {
2220         pwaff1 = isl_pw_aff_intersect_domain(pwaff1, cond1);
2221         pwaff2 = isl_pw_aff_intersect_domain(pwaff2, cond2);
2222
2223         return isl_pw_aff_add_disjoint(pwaff1, pwaff2);
2224 }
2225
2226 /* Return an affine expression that is equal to pwaff_true for elements
2227  * where "cond" is non-zero and to pwaff_false for elements where "cond"
2228  * is zero.
2229  * That is, return cond ? pwaff_true : pwaff_false;
2230  */
2231 __isl_give isl_pw_aff *isl_pw_aff_cond(__isl_take isl_pw_aff *cond,
2232         __isl_take isl_pw_aff *pwaff_true, __isl_take isl_pw_aff *pwaff_false)
2233 {
2234         isl_set *cond_true, *cond_false;
2235
2236         cond_true = isl_pw_aff_non_zero_set(isl_pw_aff_copy(cond));
2237         cond_false = isl_pw_aff_zero_set(cond);
2238         return isl_pw_aff_select(cond_true, pwaff_true,
2239                                  cond_false, pwaff_false);
2240 }
2241
2242 int isl_aff_is_cst(__isl_keep isl_aff *aff)
2243 {
2244         if (!aff)
2245                 return -1;
2246
2247         return isl_seq_first_non_zero(aff->v->el + 2, aff->v->size - 2) == -1;
2248 }
2249
2250 /* Check whether pwaff is a piecewise constant.
2251  */
2252 int isl_pw_aff_is_cst(__isl_keep isl_pw_aff *pwaff)
2253 {
2254         int i;
2255
2256         if (!pwaff)
2257                 return -1;
2258
2259         for (i = 0; i < pwaff->n; ++i) {
2260                 int is_cst = isl_aff_is_cst(pwaff->p[i].aff);
2261                 if (is_cst < 0 || !is_cst)
2262                         return is_cst;
2263         }
2264
2265         return 1;
2266 }
2267
2268 __isl_give isl_aff *isl_aff_mul(__isl_take isl_aff *aff1,
2269         __isl_take isl_aff *aff2)
2270 {
2271         if (!isl_aff_is_cst(aff2) && isl_aff_is_cst(aff1))
2272                 return isl_aff_mul(aff2, aff1);
2273
2274         if (!isl_aff_is_cst(aff2))
2275                 isl_die(isl_aff_get_ctx(aff1), isl_error_invalid,
2276                         "at least one affine expression should be constant",
2277                         goto error);
2278
2279         aff1 = isl_aff_cow(aff1);
2280         if (!aff1 || !aff2)
2281                 goto error;
2282
2283         aff1 = isl_aff_scale(aff1, aff2->v->el[1]);
2284         aff1 = isl_aff_scale_down(aff1, aff2->v->el[0]);
2285
2286         isl_aff_free(aff2);
2287         return aff1;
2288 error:
2289         isl_aff_free(aff1);
2290         isl_aff_free(aff2);
2291         return NULL;
2292 }
2293
2294 /* Divide "aff1" by "aff2", assuming "aff2" is a piecewise constant.
2295  */
2296 __isl_give isl_aff *isl_aff_div(__isl_take isl_aff *aff1,
2297         __isl_take isl_aff *aff2)
2298 {
2299         int is_cst;
2300         int neg;
2301
2302         is_cst = isl_aff_is_cst(aff2);
2303         if (is_cst < 0)
2304                 goto error;
2305         if (!is_cst)
2306                 isl_die(isl_aff_get_ctx(aff2), isl_error_invalid,
2307                         "second argument should be a constant", goto error);
2308
2309         if (!aff2)
2310                 goto error;
2311
2312         neg = isl_int_is_neg(aff2->v->el[1]);
2313         if (neg) {
2314                 isl_int_neg(aff2->v->el[0], aff2->v->el[0]);
2315                 isl_int_neg(aff2->v->el[1], aff2->v->el[1]);
2316         }
2317
2318         aff1 = isl_aff_scale(aff1, aff2->v->el[0]);
2319         aff1 = isl_aff_scale_down(aff1, aff2->v->el[1]);
2320
2321         if (neg) {
2322                 isl_int_neg(aff2->v->el[0], aff2->v->el[0]);
2323                 isl_int_neg(aff2->v->el[1], aff2->v->el[1]);
2324         }
2325
2326         isl_aff_free(aff2);
2327         return aff1;
2328 error:
2329         isl_aff_free(aff1);
2330         isl_aff_free(aff2);
2331         return NULL;
2332 }
2333
2334 static __isl_give isl_pw_aff *pw_aff_add(__isl_take isl_pw_aff *pwaff1,
2335         __isl_take isl_pw_aff *pwaff2)
2336 {
2337         return isl_pw_aff_on_shared_domain(pwaff1, pwaff2, &isl_aff_add);
2338 }
2339
2340 __isl_give isl_pw_aff *isl_pw_aff_add(__isl_take isl_pw_aff *pwaff1,
2341         __isl_take isl_pw_aff *pwaff2)
2342 {
2343         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2, &pw_aff_add);
2344 }
2345
2346 __isl_give isl_pw_aff *isl_pw_aff_union_add(__isl_take isl_pw_aff *pwaff1,
2347         __isl_take isl_pw_aff *pwaff2)
2348 {
2349         return isl_pw_aff_union_add_(pwaff1, pwaff2);
2350 }
2351
2352 static __isl_give isl_pw_aff *pw_aff_mul(__isl_take isl_pw_aff *pwaff1,
2353         __isl_take isl_pw_aff *pwaff2)
2354 {
2355         return isl_pw_aff_on_shared_domain(pwaff1, pwaff2, &isl_aff_mul);
2356 }
2357
2358 __isl_give isl_pw_aff *isl_pw_aff_mul(__isl_take isl_pw_aff *pwaff1,
2359         __isl_take isl_pw_aff *pwaff2)
2360 {
2361         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2, &pw_aff_mul);
2362 }
2363
2364 static __isl_give isl_pw_aff *pw_aff_div(__isl_take isl_pw_aff *pa1,
2365         __isl_take isl_pw_aff *pa2)
2366 {
2367         return isl_pw_aff_on_shared_domain(pa1, pa2, &isl_aff_div);
2368 }
2369
2370 /* Divide "pa1" by "pa2", assuming "pa2" is a piecewise constant.
2371  */
2372 __isl_give isl_pw_aff *isl_pw_aff_div(__isl_take isl_pw_aff *pa1,
2373         __isl_take isl_pw_aff *pa2)
2374 {
2375         int is_cst;
2376
2377         is_cst = isl_pw_aff_is_cst(pa2);
2378         if (is_cst < 0)
2379                 goto error;
2380         if (!is_cst)
2381                 isl_die(isl_pw_aff_get_ctx(pa2), isl_error_invalid,
2382                         "second argument should be a piecewise constant",
2383                         goto error);
2384         return isl_pw_aff_align_params_pw_pw_and(pa1, pa2, &pw_aff_div);
2385 error:
2386         isl_pw_aff_free(pa1);
2387         isl_pw_aff_free(pa2);
2388         return NULL;
2389 }
2390
2391 /* Compute the quotient of the integer division of "pa1" by "pa2"
2392  * with rounding towards zero.
2393  * "pa2" is assumed to be a piecewise constant.
2394  *
2395  * In particular, return
2396  *
2397  *      pa1 >= 0 ? floor(pa1/pa2) : ceil(pa1/pa2)
2398  *
2399  */
2400 __isl_give isl_pw_aff *isl_pw_aff_tdiv_q(__isl_take isl_pw_aff *pa1,
2401         __isl_take isl_pw_aff *pa2)
2402 {
2403         int is_cst;
2404         isl_set *cond;
2405         isl_pw_aff *f, *c;
2406
2407         is_cst = isl_pw_aff_is_cst(pa2);
2408         if (is_cst < 0)
2409                 goto error;
2410         if (!is_cst)
2411                 isl_die(isl_pw_aff_get_ctx(pa2), isl_error_invalid,
2412                         "second argument should be a piecewise constant",
2413                         goto error);
2414
2415         pa1 = isl_pw_aff_div(pa1, pa2);
2416
2417         cond = isl_pw_aff_nonneg_set(isl_pw_aff_copy(pa1));
2418         f = isl_pw_aff_floor(isl_pw_aff_copy(pa1));
2419         c = isl_pw_aff_ceil(pa1);
2420         return isl_pw_aff_cond(isl_set_indicator_function(cond), f, c);
2421 error:
2422         isl_pw_aff_free(pa1);
2423         isl_pw_aff_free(pa2);
2424         return NULL;
2425 }
2426
2427 /* Compute the remainder of the integer division of "pa1" by "pa2"
2428  * with rounding towards zero.
2429  * "pa2" is assumed to be a piecewise constant.
2430  *
2431  * In particular, return
2432  *
2433  *      pa1 - pa2 * (pa1 >= 0 ? floor(pa1/pa2) : ceil(pa1/pa2))
2434  *
2435  */
2436 __isl_give isl_pw_aff *isl_pw_aff_tdiv_r(__isl_take isl_pw_aff *pa1,
2437         __isl_take isl_pw_aff *pa2)
2438 {
2439         int is_cst;
2440         isl_pw_aff *res;
2441
2442         is_cst = isl_pw_aff_is_cst(pa2);
2443         if (is_cst < 0)
2444                 goto error;
2445         if (!is_cst)
2446                 isl_die(isl_pw_aff_get_ctx(pa2), isl_error_invalid,
2447                         "second argument should be a piecewise constant",
2448                         goto error);
2449         res = isl_pw_aff_tdiv_q(isl_pw_aff_copy(pa1), isl_pw_aff_copy(pa2));
2450         res = isl_pw_aff_mul(pa2, res);
2451         res = isl_pw_aff_sub(pa1, res);
2452         return res;
2453 error:
2454         isl_pw_aff_free(pa1);
2455         isl_pw_aff_free(pa2);
2456         return NULL;
2457 }
2458
2459 static __isl_give isl_pw_aff *pw_aff_min(__isl_take isl_pw_aff *pwaff1,
2460         __isl_take isl_pw_aff *pwaff2)
2461 {
2462         isl_set *le;
2463         isl_set *dom;
2464
2465         dom = isl_set_intersect(isl_pw_aff_domain(isl_pw_aff_copy(pwaff1)),
2466                                 isl_pw_aff_domain(isl_pw_aff_copy(pwaff2)));
2467         le = isl_pw_aff_le_set(isl_pw_aff_copy(pwaff1),
2468                                 isl_pw_aff_copy(pwaff2));
2469         dom = isl_set_subtract(dom, isl_set_copy(le));
2470         return isl_pw_aff_select(le, pwaff1, dom, pwaff2);
2471 }
2472
2473 __isl_give isl_pw_aff *isl_pw_aff_min(__isl_take isl_pw_aff *pwaff1,
2474         __isl_take isl_pw_aff *pwaff2)
2475 {
2476         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2, &pw_aff_min);
2477 }
2478
2479 static __isl_give isl_pw_aff *pw_aff_max(__isl_take isl_pw_aff *pwaff1,
2480         __isl_take isl_pw_aff *pwaff2)
2481 {
2482         isl_set *ge;
2483         isl_set *dom;
2484
2485         dom = isl_set_intersect(isl_pw_aff_domain(isl_pw_aff_copy(pwaff1)),
2486                                 isl_pw_aff_domain(isl_pw_aff_copy(pwaff2)));
2487         ge = isl_pw_aff_ge_set(isl_pw_aff_copy(pwaff1),
2488                                 isl_pw_aff_copy(pwaff2));
2489         dom = isl_set_subtract(dom, isl_set_copy(ge));
2490         return isl_pw_aff_select(ge, pwaff1, dom, pwaff2);
2491 }
2492
2493 __isl_give isl_pw_aff *isl_pw_aff_max(__isl_take isl_pw_aff *pwaff1,
2494         __isl_take isl_pw_aff *pwaff2)
2495 {
2496         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2, &pw_aff_max);
2497 }
2498
2499 static __isl_give isl_pw_aff *pw_aff_list_reduce(
2500         __isl_take isl_pw_aff_list *list,
2501         __isl_give isl_pw_aff *(*fn)(__isl_take isl_pw_aff *pwaff1,
2502                                         __isl_take isl_pw_aff *pwaff2))
2503 {
2504         int i;
2505         isl_ctx *ctx;
2506         isl_pw_aff *res;
2507
2508         if (!list)
2509                 return NULL;
2510
2511         ctx = isl_pw_aff_list_get_ctx(list);
2512         if (list->n < 1)
2513                 isl_die(ctx, isl_error_invalid,
2514                         "list should contain at least one element",
2515                         return isl_pw_aff_list_free(list));
2516
2517         res = isl_pw_aff_copy(list->p[0]);
2518         for (i = 1; i < list->n; ++i)
2519                 res = fn(res, isl_pw_aff_copy(list->p[i]));
2520
2521         isl_pw_aff_list_free(list);
2522         return res;
2523 }
2524
2525 /* Return an isl_pw_aff that maps each element in the intersection of the
2526  * domains of the elements of list to the minimal corresponding affine
2527  * expression.
2528  */
2529 __isl_give isl_pw_aff *isl_pw_aff_list_min(__isl_take isl_pw_aff_list *list)
2530 {
2531         return pw_aff_list_reduce(list, &isl_pw_aff_min);
2532 }
2533
2534 /* Return an isl_pw_aff that maps each element in the intersection of the
2535  * domains of the elements of list to the maximal corresponding affine
2536  * expression.
2537  */
2538 __isl_give isl_pw_aff *isl_pw_aff_list_max(__isl_take isl_pw_aff_list *list)
2539 {
2540         return pw_aff_list_reduce(list, &isl_pw_aff_max);
2541 }
2542
2543 /* Mark the domains of "pwaff" as rational.
2544  */
2545 __isl_give isl_pw_aff *isl_pw_aff_set_rational(__isl_take isl_pw_aff *pwaff)
2546 {
2547         int i;
2548
2549         pwaff = isl_pw_aff_cow(pwaff);
2550         if (!pwaff)
2551                 return NULL;
2552         if (pwaff->n == 0)
2553                 return pwaff;
2554
2555         for (i = 0; i < pwaff->n; ++i) {
2556                 pwaff->p[i].set = isl_set_set_rational(pwaff->p[i].set);
2557                 if (!pwaff->p[i].set)
2558                         return isl_pw_aff_free(pwaff);
2559         }
2560
2561         return pwaff;
2562 }
2563
2564 /* Mark the domains of the elements of "list" as rational.
2565  */
2566 __isl_give isl_pw_aff_list *isl_pw_aff_list_set_rational(
2567         __isl_take isl_pw_aff_list *list)
2568 {
2569         int i;
2570
2571         if (!list)
2572                 return NULL;
2573         if (list->n == 0)
2574                 return list;
2575
2576         for (i = 0; i < list->n; ++i) {
2577                 isl_pw_aff *pa;
2578
2579                 pa = isl_pw_aff_list_get_pw_aff(list, i);
2580                 pa = isl_pw_aff_set_rational(pa);
2581                 list = isl_pw_aff_list_set_pw_aff(list, i, pa);
2582         }
2583
2584         return list;
2585 }
2586
2587 #undef BASE
2588 #define BASE aff
2589
2590 #include <isl_multi_templ.c>
2591
2592 /* Create an isl_pw_multi_aff with the given isl_multi_aff on a universe
2593  * domain.
2594  */
2595 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_from_multi_aff(
2596         __isl_take isl_multi_aff *ma)
2597 {
2598         isl_set *dom = isl_set_universe(isl_multi_aff_get_domain_space(ma));
2599         return isl_pw_multi_aff_alloc(dom, ma);
2600 }
2601
2602 /* Create a piecewise multi-affine expression in the given space that maps each
2603  * input dimension to the corresponding output dimension.
2604  */
2605 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_identity(
2606         __isl_take isl_space *space)
2607 {
2608         return isl_pw_multi_aff_from_multi_aff(isl_multi_aff_identity(space));
2609 }
2610
2611 __isl_give isl_multi_aff *isl_multi_aff_add(__isl_take isl_multi_aff *maff1,
2612         __isl_take isl_multi_aff *maff2)
2613 {
2614         int i;
2615         isl_ctx *ctx;
2616
2617         maff1 = isl_multi_aff_cow(maff1);
2618         if (!maff1 || !maff2)
2619                 goto error;
2620
2621         ctx = isl_multi_aff_get_ctx(maff1);
2622         if (!isl_space_is_equal(maff1->space, maff2->space))
2623                 isl_die(ctx, isl_error_invalid,
2624                         "spaces don't match", goto error);
2625
2626         for (i = 0; i < maff1->n; ++i) {
2627                 maff1->p[i] = isl_aff_add(maff1->p[i],
2628                                             isl_aff_copy(maff2->p[i]));
2629                 if (!maff1->p[i])
2630                         goto error;
2631         }
2632
2633         isl_multi_aff_free(maff2);
2634         return maff1;
2635 error:
2636         isl_multi_aff_free(maff1);
2637         isl_multi_aff_free(maff2);
2638         return NULL;
2639 }
2640
2641 /* Given two multi-affine expressions A -> B and C -> D,
2642  * construct a multi-affine expression [A -> C] -> [B -> D].
2643  */
2644 __isl_give isl_multi_aff *isl_multi_aff_product(
2645         __isl_take isl_multi_aff *ma1, __isl_take isl_multi_aff *ma2)
2646 {
2647         int i;
2648         isl_aff *aff;
2649         isl_space *space;
2650         isl_multi_aff *res;
2651         int in1, in2, out1, out2;
2652
2653         in1 = isl_multi_aff_dim(ma1, isl_dim_in);
2654         in2 = isl_multi_aff_dim(ma2, isl_dim_in);
2655         out1 = isl_multi_aff_dim(ma1, isl_dim_out);
2656         out2 = isl_multi_aff_dim(ma2, isl_dim_out);
2657         space = isl_space_product(isl_multi_aff_get_space(ma1),
2658                                   isl_multi_aff_get_space(ma2));
2659         res = isl_multi_aff_alloc(isl_space_copy(space));
2660         space = isl_space_domain(space);
2661
2662         for (i = 0; i < out1; ++i) {
2663                 aff = isl_multi_aff_get_aff(ma1, i);
2664                 aff = isl_aff_insert_dims(aff, isl_dim_in, in1, in2);
2665                 aff = isl_aff_reset_domain_space(aff, isl_space_copy(space));
2666                 res = isl_multi_aff_set_aff(res, i, aff);
2667         }
2668
2669         for (i = 0; i < out2; ++i) {
2670                 aff = isl_multi_aff_get_aff(ma2, i);
2671                 aff = isl_aff_insert_dims(aff, isl_dim_in, 0, in1);
2672                 aff = isl_aff_reset_domain_space(aff, isl_space_copy(space));
2673                 res = isl_multi_aff_set_aff(res, out1 + i, aff);
2674         }
2675
2676         isl_space_free(space);
2677         isl_multi_aff_free(ma1);
2678         isl_multi_aff_free(ma2);
2679         return res;
2680 }
2681
2682 /* Exploit the equalities in "eq" to simplify the affine expressions.
2683  */
2684 static __isl_give isl_multi_aff *isl_multi_aff_substitute_equalities(
2685         __isl_take isl_multi_aff *maff, __isl_take isl_basic_set *eq)
2686 {
2687         int i;
2688
2689         maff = isl_multi_aff_cow(maff);
2690         if (!maff || !eq)
2691                 goto error;
2692
2693         for (i = 0; i < maff->n; ++i) {
2694                 maff->p[i] = isl_aff_substitute_equalities(maff->p[i],
2695                                                     isl_basic_set_copy(eq));
2696                 if (!maff->p[i])
2697                         goto error;
2698         }
2699
2700         isl_basic_set_free(eq);
2701         return maff;
2702 error:
2703         isl_basic_set_free(eq);
2704         isl_multi_aff_free(maff);
2705         return NULL;
2706 }
2707
2708 __isl_give isl_multi_aff *isl_multi_aff_scale(__isl_take isl_multi_aff *maff,
2709         isl_int f)
2710 {
2711         int i;
2712
2713         maff = isl_multi_aff_cow(maff);
2714         if (!maff)
2715                 return NULL;
2716
2717         for (i = 0; i < maff->n; ++i) {
2718                 maff->p[i] = isl_aff_scale(maff->p[i], f);
2719                 if (!maff->p[i])
2720                         return isl_multi_aff_free(maff);
2721         }
2722
2723         return maff;
2724 }
2725
2726 __isl_give isl_multi_aff *isl_multi_aff_add_on_domain(__isl_keep isl_set *dom,
2727         __isl_take isl_multi_aff *maff1, __isl_take isl_multi_aff *maff2)
2728 {
2729         maff1 = isl_multi_aff_add(maff1, maff2);
2730         maff1 = isl_multi_aff_gist(maff1, isl_set_copy(dom));
2731         return maff1;
2732 }
2733
2734 int isl_multi_aff_is_empty(__isl_keep isl_multi_aff *maff)
2735 {
2736         if (!maff)
2737                 return -1;
2738
2739         return 0;
2740 }
2741
2742 int isl_multi_aff_plain_is_equal(__isl_keep isl_multi_aff *maff1,
2743         __isl_keep isl_multi_aff *maff2)
2744 {
2745         int i;
2746         int equal;
2747
2748         if (!maff1 || !maff2)
2749                 return -1;
2750         if (maff1->n != maff2->n)
2751                 return 0;
2752         equal = isl_space_is_equal(maff1->space, maff2->space);
2753         if (equal < 0 || !equal)
2754                 return equal;
2755
2756         for (i = 0; i < maff1->n; ++i) {
2757                 equal = isl_aff_plain_is_equal(maff1->p[i], maff2->p[i]);
2758                 if (equal < 0 || !equal)
2759                         return equal;
2760         }
2761
2762         return 1;
2763 }
2764
2765 /* Return the set of domain elements where "ma1" is lexicographically
2766  * smaller than or equal to "ma2".
2767  */
2768 __isl_give isl_set *isl_multi_aff_lex_le_set(__isl_take isl_multi_aff *ma1,
2769         __isl_take isl_multi_aff *ma2)
2770 {
2771         return isl_multi_aff_lex_ge_set(ma2, ma1);
2772 }
2773
2774 /* Return the set of domain elements where "ma1" is lexicographically
2775  * greater than or equal to "ma2".
2776  */
2777 __isl_give isl_set *isl_multi_aff_lex_ge_set(__isl_take isl_multi_aff *ma1,
2778         __isl_take isl_multi_aff *ma2)
2779 {
2780         isl_space *space;
2781         isl_map *map1, *map2;
2782         isl_map *map, *ge;
2783
2784         map1 = isl_map_from_multi_aff(ma1);
2785         map2 = isl_map_from_multi_aff(ma2);
2786         map = isl_map_range_product(map1, map2);
2787         space = isl_space_range(isl_map_get_space(map));
2788         space = isl_space_domain(isl_space_unwrap(space));
2789         ge = isl_map_lex_ge(space);
2790         map = isl_map_intersect_range(map, isl_map_wrap(ge));
2791
2792         return isl_map_domain(map);
2793 }
2794
2795 #undef PW
2796 #define PW isl_pw_multi_aff
2797 #undef EL
2798 #define EL isl_multi_aff
2799 #undef EL_IS_ZERO
2800 #define EL_IS_ZERO is_empty
2801 #undef ZERO
2802 #define ZERO empty
2803 #undef IS_ZERO
2804 #define IS_ZERO is_empty
2805 #undef FIELD
2806 #define FIELD maff
2807 #undef DEFAULT_IS_ZERO
2808 #define DEFAULT_IS_ZERO 0
2809
2810 #define NO_NEG
2811 #define NO_EVAL
2812 #define NO_OPT
2813 #define NO_INVOLVES_DIMS
2814 #define NO_MOVE_DIMS
2815 #define NO_INSERT_DIMS
2816 #define NO_LIFT
2817 #define NO_MORPH
2818
2819 #include <isl_pw_templ.c>
2820
2821 #undef UNION
2822 #define UNION isl_union_pw_multi_aff
2823 #undef PART
2824 #define PART isl_pw_multi_aff
2825 #undef PARTS
2826 #define PARTS pw_multi_aff
2827 #define ALIGN_DOMAIN
2828
2829 #define NO_EVAL
2830
2831 #include <isl_union_templ.c>
2832
2833 /* Given a function "cmp" that returns the set of elements where
2834  * "ma1" is "better" than "ma2", return the intersection of this
2835  * set with "dom1" and "dom2".
2836  */
2837 static __isl_give isl_set *shared_and_better(__isl_keep isl_set *dom1,
2838         __isl_keep isl_set *dom2, __isl_keep isl_multi_aff *ma1,
2839         __isl_keep isl_multi_aff *ma2,
2840         __isl_give isl_set *(*cmp)(__isl_take isl_multi_aff *ma1,
2841                                     __isl_take isl_multi_aff *ma2))
2842 {
2843         isl_set *common;
2844         isl_set *better;
2845         int is_empty;
2846
2847         common = isl_set_intersect(isl_set_copy(dom1), isl_set_copy(dom2));
2848         is_empty = isl_set_plain_is_empty(common);
2849         if (is_empty >= 0 && is_empty)
2850                 return common;
2851         if (is_empty < 0)
2852                 return isl_set_free(common);
2853         better = cmp(isl_multi_aff_copy(ma1), isl_multi_aff_copy(ma2));
2854         better = isl_set_intersect(common, better);
2855
2856         return better;
2857 }
2858
2859 /* Given a function "cmp" that returns the set of elements where
2860  * "ma1" is "better" than "ma2", return a piecewise multi affine
2861  * expression defined on the union of the definition domains
2862  * of "pma1" and "pma2" that maps to the "best" of "pma1" and
2863  * "pma2" on each cell.  If only one of the two input functions
2864  * is defined on a given cell, then it is considered the best.
2865  */
2866 static __isl_give isl_pw_multi_aff *pw_multi_aff_union_opt(
2867         __isl_take isl_pw_multi_aff *pma1,
2868         __isl_take isl_pw_multi_aff *pma2,
2869         __isl_give isl_set *(*cmp)(__isl_take isl_multi_aff *ma1,
2870                                     __isl_take isl_multi_aff *ma2))
2871 {
2872         int i, j, n;
2873         isl_pw_multi_aff *res = NULL;
2874         isl_ctx *ctx;
2875         isl_set *set = NULL;
2876
2877         if (!pma1 || !pma2)
2878                 goto error;
2879
2880         ctx = isl_space_get_ctx(pma1->dim);
2881         if (!isl_space_is_equal(pma1->dim, pma2->dim))
2882                 isl_die(ctx, isl_error_invalid,
2883                         "arguments should live in the same space", goto error);
2884
2885         if (isl_pw_multi_aff_is_empty(pma1)) {
2886                 isl_pw_multi_aff_free(pma1);
2887                 return pma2;
2888         }
2889
2890         if (isl_pw_multi_aff_is_empty(pma2)) {
2891                 isl_pw_multi_aff_free(pma2);
2892                 return pma1;
2893         }
2894
2895         n = 2 * (pma1->n + 1) * (pma2->n + 1);
2896         res = isl_pw_multi_aff_alloc_size(isl_space_copy(pma1->dim), n);
2897
2898         for (i = 0; i < pma1->n; ++i) {
2899                 set = isl_set_copy(pma1->p[i].set);
2900                 for (j = 0; j < pma2->n; ++j) {
2901                         isl_set *better;
2902                         int is_empty;
2903
2904                         better = shared_and_better(pma2->p[j].set,
2905                                         pma1->p[i].set, pma2->p[j].maff,
2906                                         pma1->p[i].maff, cmp);
2907                         is_empty = isl_set_plain_is_empty(better);
2908                         if (is_empty < 0 || is_empty) {
2909                                 isl_set_free(better);
2910                                 if (is_empty < 0)
2911                                         goto error;
2912                                 continue;
2913                         }
2914                         set = isl_set_subtract(set, isl_set_copy(better));
2915
2916                         res = isl_pw_multi_aff_add_piece(res, better,
2917                                         isl_multi_aff_copy(pma2->p[j].maff));
2918                 }
2919                 res = isl_pw_multi_aff_add_piece(res, set,
2920                                         isl_multi_aff_copy(pma1->p[i].maff));
2921         }
2922
2923         for (j = 0; j < pma2->n; ++j) {
2924                 set = isl_set_copy(pma2->p[j].set);
2925                 for (i = 0; i < pma1->n; ++i)
2926                         set = isl_set_subtract(set,
2927                                         isl_set_copy(pma1->p[i].set));
2928                 res = isl_pw_multi_aff_add_piece(res, set,
2929                                         isl_multi_aff_copy(pma2->p[j].maff));
2930         }
2931
2932         isl_pw_multi_aff_free(pma1);
2933         isl_pw_multi_aff_free(pma2);
2934
2935         return res;
2936 error:
2937         isl_pw_multi_aff_free(pma1);
2938         isl_pw_multi_aff_free(pma2);
2939         isl_set_free(set);
2940         return isl_pw_multi_aff_free(res);
2941 }
2942
2943 static __isl_give isl_pw_multi_aff *pw_multi_aff_union_lexmax(
2944         __isl_take isl_pw_multi_aff *pma1,
2945         __isl_take isl_pw_multi_aff *pma2)
2946 {
2947         return pw_multi_aff_union_opt(pma1, pma2, &isl_multi_aff_lex_ge_set);
2948 }
2949
2950 /* Given two piecewise multi affine expressions, return a piecewise
2951  * multi-affine expression defined on the union of the definition domains
2952  * of the inputs that is equal to the lexicographic maximum of the two
2953  * inputs on each cell.  If only one of the two inputs is defined on
2954  * a given cell, then it is considered to be the maximum.
2955  */
2956 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_union_lexmax(
2957         __isl_take isl_pw_multi_aff *pma1,
2958         __isl_take isl_pw_multi_aff *pma2)
2959 {
2960         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
2961                                                     &pw_multi_aff_union_lexmax);
2962 }
2963
2964 static __isl_give isl_pw_multi_aff *pw_multi_aff_union_lexmin(
2965         __isl_take isl_pw_multi_aff *pma1,
2966         __isl_take isl_pw_multi_aff *pma2)
2967 {
2968         return pw_multi_aff_union_opt(pma1, pma2, &isl_multi_aff_lex_le_set);
2969 }
2970
2971 /* Given two piecewise multi affine expressions, return a piecewise
2972  * multi-affine expression defined on the union of the definition domains
2973  * of the inputs that is equal to the lexicographic minimum of the two
2974  * inputs on each cell.  If only one of the two inputs is defined on
2975  * a given cell, then it is considered to be the minimum.
2976  */
2977 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_union_lexmin(
2978         __isl_take isl_pw_multi_aff *pma1,
2979         __isl_take isl_pw_multi_aff *pma2)
2980 {
2981         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
2982                                                     &pw_multi_aff_union_lexmin);
2983 }
2984
2985 static __isl_give isl_pw_multi_aff *pw_multi_aff_add(
2986         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
2987 {
2988         return isl_pw_multi_aff_on_shared_domain(pma1, pma2,
2989                                                 &isl_multi_aff_add);
2990 }
2991
2992 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_add(
2993         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
2994 {
2995         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
2996                                                 &pw_multi_aff_add);
2997 }
2998
2999 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_union_add(
3000         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
3001 {
3002         return isl_pw_multi_aff_union_add_(pma1, pma2);
3003 }
3004
3005 /* Given two piecewise multi-affine expressions A -> B and C -> D,
3006  * construct a piecewise multi-affine expression [A -> C] -> [B -> D].
3007  */
3008 static __isl_give isl_pw_multi_aff *pw_multi_aff_product(
3009         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
3010 {
3011         int i, j, n;
3012         isl_space *space;
3013         isl_pw_multi_aff *res;
3014
3015         if (!pma1 || !pma2)
3016                 goto error;
3017
3018         n = pma1->n * pma2->n;
3019         space = isl_space_product(isl_space_copy(pma1->dim),
3020                                   isl_space_copy(pma2->dim));
3021         res = isl_pw_multi_aff_alloc_size(space, n);
3022
3023         for (i = 0; i < pma1->n; ++i) {
3024                 for (j = 0; j < pma2->n; ++j) {
3025                         isl_set *domain;
3026                         isl_multi_aff *ma;
3027
3028                         domain = isl_set_product(isl_set_copy(pma1->p[i].set),
3029                                                  isl_set_copy(pma2->p[j].set));
3030                         ma = isl_multi_aff_product(
3031                                         isl_multi_aff_copy(pma1->p[i].maff),
3032                                         isl_multi_aff_copy(pma2->p[i].maff));
3033                         res = isl_pw_multi_aff_add_piece(res, domain, ma);
3034                 }
3035         }
3036
3037         isl_pw_multi_aff_free(pma1);
3038         isl_pw_multi_aff_free(pma2);
3039         return res;
3040 error:
3041         isl_pw_multi_aff_free(pma1);
3042         isl_pw_multi_aff_free(pma2);
3043         return NULL;
3044 }
3045
3046 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_product(
3047         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
3048 {
3049         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
3050                                                 &pw_multi_aff_product);
3051 }
3052
3053 /* Construct a map mapping the domain of the piecewise multi-affine expression
3054  * to its range, with each dimension in the range equated to the
3055  * corresponding affine expression on its cell.
3056  */
3057 __isl_give isl_map *isl_map_from_pw_multi_aff(__isl_take isl_pw_multi_aff *pma)
3058 {
3059         int i;
3060         isl_map *map;
3061
3062         if (!pma)
3063                 return NULL;
3064
3065         map = isl_map_empty(isl_pw_multi_aff_get_space(pma));
3066
3067         for (i = 0; i < pma->n; ++i) {
3068                 isl_multi_aff *maff;
3069                 isl_basic_map *bmap;
3070                 isl_map *map_i;
3071
3072                 maff = isl_multi_aff_copy(pma->p[i].maff);
3073                 bmap = isl_basic_map_from_multi_aff(maff);
3074                 map_i = isl_map_from_basic_map(bmap);
3075                 map_i = isl_map_intersect_domain(map_i,
3076                                                 isl_set_copy(pma->p[i].set));
3077                 map = isl_map_union_disjoint(map, map_i);
3078         }
3079
3080         isl_pw_multi_aff_free(pma);
3081         return map;
3082 }
3083
3084 __isl_give isl_set *isl_set_from_pw_multi_aff(__isl_take isl_pw_multi_aff *pma)
3085 {
3086         if (!pma)
3087                 return NULL;
3088
3089         if (!isl_space_is_set(pma->dim))
3090                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
3091                         "isl_pw_multi_aff cannot be converted into an isl_set",
3092                         return isl_pw_multi_aff_free(pma));
3093
3094         return isl_map_from_pw_multi_aff(pma);
3095 }
3096
3097 /* Given a basic map with a single output dimension that is defined
3098  * in terms of the parameters and input dimensions using an equality,
3099  * extract an isl_aff that expresses the output dimension in terms
3100  * of the parameters and input dimensions.
3101  *
3102  * Since some applications expect the result of isl_pw_multi_aff_from_map
3103  * to only contain integer affine expressions, we compute the floor
3104  * of the expression before returning.
3105  *
3106  * This function shares some similarities with
3107  * isl_basic_map_has_defining_equality and isl_constraint_get_bound.
3108  */
3109 static __isl_give isl_aff *extract_isl_aff_from_basic_map(
3110         __isl_take isl_basic_map *bmap)
3111 {
3112         int i;
3113         unsigned offset;
3114         unsigned total;
3115         isl_local_space *ls;
3116         isl_aff *aff;
3117
3118         if (!bmap)
3119                 return NULL;
3120         if (isl_basic_map_dim(bmap, isl_dim_out) != 1)
3121                 isl_die(isl_basic_map_get_ctx(bmap), isl_error_invalid,
3122                         "basic map should have a single output dimension",
3123                         goto error);
3124         offset = isl_basic_map_offset(bmap, isl_dim_out);
3125         total = isl_basic_map_total_dim(bmap);
3126         for (i = 0; i < bmap->n_eq; ++i) {
3127                 if (isl_int_is_zero(bmap->eq[i][offset]))
3128                         continue;
3129                 if (isl_seq_first_non_zero(bmap->eq[i] + offset + 1,
3130                                            1 + total - (offset + 1)) != -1)
3131                         continue;
3132                 break;
3133         }
3134         if (i >= bmap->n_eq)
3135                 isl_die(isl_basic_map_get_ctx(bmap), isl_error_invalid,
3136                         "unable to find suitable equality", goto error);
3137         ls = isl_basic_map_get_local_space(bmap);
3138         aff = isl_aff_alloc(isl_local_space_domain(ls));
3139         if (!aff)
3140                 goto error;
3141         if (isl_int_is_neg(bmap->eq[i][offset]))
3142                 isl_seq_cpy(aff->v->el + 1, bmap->eq[i], offset);
3143         else
3144                 isl_seq_neg(aff->v->el + 1, bmap->eq[i], offset);
3145         isl_seq_clr(aff->v->el + 1 + offset, aff->v->size - (1 + offset));
3146         isl_int_abs(aff->v->el[0], bmap->eq[i][offset]);
3147         isl_basic_map_free(bmap);
3148
3149         aff = isl_aff_remove_unused_divs(aff);
3150         aff = isl_aff_floor(aff);
3151         return aff;
3152 error:
3153         isl_basic_map_free(bmap);
3154         return NULL;
3155 }
3156
3157 /* Given a basic map where each output dimension is defined
3158  * in terms of the parameters and input dimensions using an equality,
3159  * extract an isl_multi_aff that expresses the output dimensions in terms
3160  * of the parameters and input dimensions.
3161  */
3162 static __isl_give isl_multi_aff *extract_isl_multi_aff_from_basic_map(
3163         __isl_take isl_basic_map *bmap)
3164 {
3165         int i;
3166         unsigned n_out;
3167         isl_multi_aff *ma;
3168
3169         if (!bmap)
3170                 return NULL;
3171
3172         ma = isl_multi_aff_alloc(isl_basic_map_get_space(bmap));
3173         n_out = isl_basic_map_dim(bmap, isl_dim_out);
3174
3175         for (i = 0; i < n_out; ++i) {
3176                 isl_basic_map *bmap_i;
3177                 isl_aff *aff;
3178
3179                 bmap_i = isl_basic_map_copy(bmap);
3180                 bmap_i = isl_basic_map_project_out(bmap_i, isl_dim_out,
3181                                                         i + 1, n_out - (1 + i));
3182                 bmap_i = isl_basic_map_project_out(bmap_i, isl_dim_out, 0, i);
3183                 aff = extract_isl_aff_from_basic_map(bmap_i);
3184                 ma = isl_multi_aff_set_aff(ma, i, aff);
3185         }
3186
3187         isl_basic_map_free(bmap);
3188
3189         return ma;
3190 }
3191
3192 /* Create an isl_pw_multi_aff that is equivalent to
3193  * isl_map_intersect_domain(isl_map_from_basic_map(bmap), domain).
3194  * The given basic map is such that each output dimension is defined
3195  * in terms of the parameters and input dimensions using an equality.
3196  */
3197 static __isl_give isl_pw_multi_aff *plain_pw_multi_aff_from_map(
3198         __isl_take isl_set *domain, __isl_take isl_basic_map *bmap)
3199 {
3200         isl_multi_aff *ma;
3201
3202         ma = extract_isl_multi_aff_from_basic_map(bmap);
3203         return isl_pw_multi_aff_alloc(domain, ma);
3204 }
3205
3206 /* Try and create an isl_pw_multi_aff that is equivalent to the given isl_map.
3207  * This obviously only works if the input "map" is single-valued.
3208  * If so, we compute the lexicographic minimum of the image in the form
3209  * of an isl_pw_multi_aff.  Since the image is unique, it is equal
3210  * to its lexicographic minimum.
3211  * If the input is not single-valued, we produce an error.
3212  */
3213 static __isl_give isl_pw_multi_aff *pw_multi_aff_from_map_base(
3214         __isl_take isl_map *map)
3215 {
3216         int i;
3217         int sv;
3218         isl_pw_multi_aff *pma;
3219
3220         sv = isl_map_is_single_valued(map);
3221         if (sv < 0)
3222                 goto error;
3223         if (!sv)
3224                 isl_die(isl_map_get_ctx(map), isl_error_invalid,
3225                         "map is not single-valued", goto error);
3226         map = isl_map_make_disjoint(map);
3227         if (!map)
3228                 return NULL;
3229
3230         pma = isl_pw_multi_aff_empty(isl_map_get_space(map));
3231
3232         for (i = 0; i < map->n; ++i) {
3233                 isl_pw_multi_aff *pma_i;
3234                 isl_basic_map *bmap;
3235                 bmap = isl_basic_map_copy(map->p[i]);
3236                 pma_i = isl_basic_map_lexmin_pw_multi_aff(bmap);
3237                 pma = isl_pw_multi_aff_add_disjoint(pma, pma_i);
3238         }
3239
3240         isl_map_free(map);
3241         return pma;
3242 error:
3243         isl_map_free(map);
3244         return NULL;
3245 }
3246
3247 /* Try and create an isl_pw_multi_aff that is equivalent to the given isl_map,
3248  * taking into account that the output dimension at position "d"
3249  * can be represented as
3250  *
3251  *      x = floor((e(...) + c1) / m)
3252  *
3253  * given that constraint "i" is of the form
3254  *
3255  *      e(...) + c1 - m x >= 0
3256  *
3257  *
3258  * Let "map" be of the form
3259  *
3260  *      A -> B
3261  *
3262  * We construct a mapping
3263  *
3264  *      A -> [A -> x = floor(...)]
3265  *
3266  * apply that to the map, obtaining
3267  *
3268  *      [A -> x = floor(...)] -> B
3269  *
3270  * and equate dimension "d" to x.
3271  * We then compute a isl_pw_multi_aff representation of the resulting map
3272  * and plug in the mapping above.
3273  */
3274 static __isl_give isl_pw_multi_aff *pw_multi_aff_from_map_div(
3275         __isl_take isl_map *map, __isl_take isl_basic_map *hull, int d, int i)
3276 {
3277         isl_ctx *ctx;
3278         isl_space *space;
3279         isl_local_space *ls;
3280         isl_multi_aff *ma;
3281         isl_aff *aff;
3282         isl_vec *v;
3283         isl_map *insert;
3284         int offset;
3285         int n;
3286         int n_in;
3287         isl_pw_multi_aff *pma;
3288         int is_set;
3289
3290         is_set = isl_map_is_set(map);
3291
3292         offset = isl_basic_map_offset(hull, isl_dim_out);
3293         ctx = isl_map_get_ctx(map);
3294         space = isl_space_domain(isl_map_get_space(map));
3295         n_in = isl_space_dim(space, isl_dim_set);
3296         n = isl_space_dim(space, isl_dim_all);
3297
3298         v = isl_vec_alloc(ctx, 1 + 1 + n);
3299         if (v) {
3300                 isl_int_neg(v->el[0], hull->ineq[i][offset + d]);
3301                 isl_seq_cpy(v->el + 1, hull->ineq[i], 1 + n);
3302         }
3303         isl_basic_map_free(hull);
3304
3305         ls = isl_local_space_from_space(isl_space_copy(space));
3306         aff = isl_aff_alloc_vec(ls, v);
3307         aff = isl_aff_floor(aff);
3308         if (is_set) {
3309                 isl_space_free(space);
3310                 ma = isl_multi_aff_from_aff(aff);
3311         } else {
3312                 ma = isl_multi_aff_identity(isl_space_map_from_set(space));
3313                 ma = isl_multi_aff_range_product(ma,
3314                                                 isl_multi_aff_from_aff(aff));
3315         }
3316
3317         insert = isl_map_from_multi_aff(isl_multi_aff_copy(ma));
3318         map = isl_map_apply_domain(map, insert);
3319         map = isl_map_equate(map, isl_dim_in, n_in, isl_dim_out, d);
3320         pma = isl_pw_multi_aff_from_map(map);
3321         pma = isl_pw_multi_aff_pullback_multi_aff(pma, ma);
3322
3323         return pma;
3324 }
3325
3326 /* Is constraint "c" of the form
3327  *
3328  *      e(...) + c1 - m x >= 0
3329  *
3330  * or
3331  *
3332  *      -e(...) + c2 + m x >= 0
3333  *
3334  * where m > 1 and e only depends on parameters and input dimemnsions?
3335  *
3336  * "offset" is the offset of the output dimensions
3337  * "pos" is the position of output dimension x.
3338  */
3339 static int is_potential_div_constraint(isl_int *c, int offset, int d, int total)
3340 {
3341         if (isl_int_is_zero(c[offset + d]))
3342                 return 0;
3343         if (isl_int_is_one(c[offset + d]))
3344                 return 0;
3345         if (isl_int_is_negone(c[offset + d]))
3346                 return 0;
3347         if (isl_seq_first_non_zero(c + offset, d) != -1)
3348                 return 0;
3349         if (isl_seq_first_non_zero(c + offset + d + 1,
3350                                     total - (offset + d + 1)) != -1)
3351                 return 0;
3352         return 1;
3353 }
3354
3355 /* Try and create an isl_pw_multi_aff that is equivalent to the given isl_map.
3356  *
3357  * As a special case, we first check if there is any pair of constraints,
3358  * shared by all the basic maps in "map" that force a given dimension
3359  * to be equal to the floor of some affine combination of the input dimensions.
3360  *
3361  * In particular, if we can find two constraints
3362  *
3363  *      e(...) + c1 - m x >= 0          i.e.,           m x <= e(...) + c1
3364  *
3365  * and
3366  *
3367  *      -e(...) + c2 + m x >= 0         i.e.,           m x >= e(...) - c2
3368  *
3369  * where m > 1 and e only depends on parameters and input dimemnsions,
3370  * and such that
3371  *
3372  *      c1 + c2 < m                     i.e.,           -c2 >= c1 - (m - 1)
3373  *
3374  * then we know that we can take
3375  *
3376  *      x = floor((e(...) + c1) / m)
3377  *
3378  * without having to perform any computation.
3379  *
3380  * Note that we know that
3381  *
3382  *      c1 + c2 >= 1
3383  *
3384  * If c1 + c2 were 0, then we would have detected an equality during
3385  * simplification.  If c1 + c2 were negative, then we would have detected
3386  * a contradiction.
3387  */
3388 static __isl_give isl_pw_multi_aff *pw_multi_aff_from_map_check_div(
3389         __isl_take isl_map *map)
3390 {
3391         int d, dim;
3392         int i, j, n;
3393         int offset, total;
3394         isl_int sum;
3395         isl_basic_map *hull;
3396
3397         hull = isl_map_unshifted_simple_hull(isl_map_copy(map));
3398         if (!hull)
3399                 goto error;
3400
3401         isl_int_init(sum);
3402         dim = isl_map_dim(map, isl_dim_out);
3403         offset = isl_basic_map_offset(hull, isl_dim_out);
3404         total = 1 + isl_basic_map_total_dim(hull);
3405         n = hull->n_ineq;
3406         for (d = 0; d < dim; ++d) {
3407                 for (i = 0; i < n; ++i) {
3408                         if (!is_potential_div_constraint(hull->ineq[i],
3409                                                         offset, d, total))
3410                                 continue;
3411                         for (j = i + 1; j < n; ++j) {
3412                                 if (!isl_seq_is_neg(hull->ineq[i] + 1,
3413                                                 hull->ineq[j] + 1, total - 1))
3414                                         continue;
3415                                 isl_int_add(sum, hull->ineq[i][0],
3416                                                 hull->ineq[j][0]);
3417                                 if (isl_int_abs_lt(sum,
3418                                                     hull->ineq[i][offset + d]))
3419                                         break;
3420
3421                         }
3422                         if (j >= n)
3423                                 continue;
3424                         isl_int_clear(sum);
3425                         if (isl_int_is_pos(hull->ineq[j][offset + d]))
3426                                 j = i;
3427                         return pw_multi_aff_from_map_div(map, hull, d, j);
3428                 }
3429         }
3430         isl_int_clear(sum);
3431         isl_basic_map_free(hull);
3432         return pw_multi_aff_from_map_base(map);
3433 error:
3434         isl_map_free(map);
3435         isl_basic_map_free(hull);
3436         return NULL;
3437 }
3438
3439 /* Given an affine expression
3440  *
3441  *      [A -> B] -> f(A,B)
3442  *
3443  * construct an isl_multi_aff
3444  *
3445  *      [A -> B] -> B'
3446  *
3447  * such that dimension "d" in B' is set to "aff" and the remaining
3448  * dimensions are set equal to the corresponding dimensions in B.
3449  * "n_in" is the dimension of the space A.
3450  * "n_out" is the dimension of the space B.
3451  *
3452  * If "is_set" is set, then the affine expression is of the form
3453  *
3454  *      [B] -> f(B)
3455  *
3456  * and we construct an isl_multi_aff
3457  *
3458  *      B -> B'
3459  */
3460 static __isl_give isl_multi_aff *range_map(__isl_take isl_aff *aff, int d,
3461         unsigned n_in, unsigned n_out, int is_set)
3462 {
3463         int i;
3464         isl_multi_aff *ma;
3465         isl_space *space, *space2;
3466         isl_local_space *ls;
3467
3468         space = isl_aff_get_domain_space(aff);
3469         ls = isl_local_space_from_space(isl_space_copy(space));
3470         space2 = isl_space_copy(space);
3471         if (!is_set)
3472                 space2 = isl_space_range(isl_space_unwrap(space2));
3473         space = isl_space_map_from_domain_and_range(space, space2);
3474         ma = isl_multi_aff_alloc(space);
3475         ma = isl_multi_aff_set_aff(ma, d, aff);
3476
3477         for (i = 0; i < n_out; ++i) {
3478                 if (i == d)
3479                         continue;
3480                 aff = isl_aff_var_on_domain(isl_local_space_copy(ls),
3481                                                 isl_dim_set, n_in + i);
3482                 ma = isl_multi_aff_set_aff(ma, i, aff);
3483         }
3484
3485         isl_local_space_free(ls);
3486
3487         return ma;
3488 }
3489
3490 /* Try and create an isl_pw_multi_aff that is equivalent to the given isl_map,
3491  * taking into account that the dimension at position "d" can be written as
3492  *
3493  *      x = m a + f(..)                                         (1)
3494  *
3495  * where m is equal to "gcd".
3496  * "i" is the index of the equality in "hull" that defines f(..).
3497  * In particular, the equality is of the form
3498  *
3499  *      f(..) - x + m g(existentials) = 0
3500  *
3501  * or
3502  *
3503  *      -f(..) + x + m g(existentials) = 0
3504  *
3505  * We basically plug (1) into "map", resulting in a map with "a"
3506  * in the range instead of "x".  The corresponding isl_pw_multi_aff
3507  * defining "a" is then plugged back into (1) to obtain a definition fro "x".
3508  *
3509  * Specifically, given the input map
3510  *
3511  *      A -> B
3512  *
3513  * We first wrap it into a set
3514  *
3515  *      [A -> B]
3516  *
3517  * and define (1) on top of the corresponding space, resulting in "aff".
3518  * We use this to create an isl_multi_aff that maps the output position "d"
3519  * from "a" to "x", leaving all other (intput and output) dimensions unchanged.
3520  * We plug this into the wrapped map, unwrap the result and compute the
3521  * corresponding isl_pw_multi_aff.
3522  * The result is an expression
3523  *
3524  *      A -> T(A)
3525  *
3526  * We adjust that to
3527  *
3528  *      A -> [A -> T(A)]
3529  *
3530  * so that we can plug that into "aff", after extending the latter to
3531  * a mapping
3532  *
3533  *      [A -> B] -> B'
3534  *
3535  *
3536  * If "map" is actually a set, then there is no "A" space, meaning
3537  * that we do not need to perform any wrapping, and that the result
3538  * of the recursive call is of the form
3539  *
3540  *      [T]
3541  *
3542  * which is plugged into a mapping of the form
3543  *
3544  *      B -> B'
3545  */
3546 static __isl_give isl_pw_multi_aff *pw_multi_aff_from_map_stride(
3547         __isl_take isl_map *map, __isl_take isl_basic_map *hull, int d, int i,
3548         isl_int gcd)
3549 {
3550         isl_set *set;
3551         isl_space *space;
3552         isl_local_space *ls;
3553         isl_aff *aff;
3554         isl_multi_aff *ma;
3555         isl_pw_multi_aff *pma, *id;
3556         unsigned n_in;
3557         unsigned o_out;
3558         unsigned n_out;
3559         int is_set;
3560
3561         is_set = isl_map_is_set(map);
3562
3563         n_in = isl_basic_map_dim(hull, isl_dim_in);
3564         n_out = isl_basic_map_dim(hull, isl_dim_out);
3565         o_out = isl_basic_map_offset(hull, isl_dim_out);
3566
3567         if (is_set)
3568                 set = map;
3569         else
3570                 set = isl_map_wrap(map);
3571         space = isl_space_map_from_set(isl_set_get_space(set));
3572         ma = isl_multi_aff_identity(space);
3573         ls = isl_local_space_from_space(isl_set_get_space(set));
3574         aff = isl_aff_alloc(ls);
3575         if (aff) {
3576                 isl_int_set_si(aff->v->el[0], 1);
3577                 if (isl_int_is_one(hull->eq[i][o_out + d]))
3578                         isl_seq_neg(aff->v->el + 1, hull->eq[i],
3579                                     aff->v->size - 1);
3580                 else
3581                         isl_seq_cpy(aff->v->el + 1, hull->eq[i],
3582                                     aff->v->size - 1);
3583                 isl_int_set(aff->v->el[1 + o_out + d], gcd);
3584         }
3585         ma = isl_multi_aff_set_aff(ma, n_in + d, isl_aff_copy(aff));
3586         set = isl_set_preimage_multi_aff(set, ma);
3587
3588         ma = range_map(aff, d, n_in, n_out, is_set);
3589
3590         if (is_set)
3591                 map = set;
3592         else
3593                 map = isl_set_unwrap(set);
3594         pma = isl_pw_multi_aff_from_map(set);
3595
3596         if (!is_set) {
3597                 space = isl_pw_multi_aff_get_domain_space(pma);
3598                 space = isl_space_map_from_set(space);
3599                 id = isl_pw_multi_aff_identity(space);
3600                 pma = isl_pw_multi_aff_range_product(id, pma);
3601         }
3602         id = isl_pw_multi_aff_from_multi_aff(ma);
3603         pma = isl_pw_multi_aff_pullback_pw_multi_aff(id, pma);
3604
3605         isl_basic_map_free(hull);
3606         return pma;
3607 }
3608
3609 /* Try and create an isl_pw_multi_aff that is equivalent to the given isl_map.
3610  *
3611  * As a special case, we first check if all output dimensions are uniquely
3612  * defined in terms of the parameters and input dimensions over the entire
3613  * domain.  If so, we extract the desired isl_pw_multi_aff directly
3614  * from the affine hull of "map" and its domain.
3615  *
3616  * Otherwise, we check if any of the output dimensions is "strided".
3617  * That is, we check if can be written as
3618  *
3619  *      x = m a + f(..)
3620  *
3621  * with m greater than 1, a some combination of existentiall quantified
3622  * variables and f and expression in the parameters and input dimensions.
3623  * If so, we remove the stride in pw_multi_aff_from_map_stride.
3624  *
3625  * Otherwise, we continue with pw_multi_aff_from_map_check_div for a further
3626  * special case.
3627  */
3628 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_from_map(__isl_take isl_map *map)
3629 {
3630         int i, j;
3631         int sv;
3632         isl_basic_map *hull;
3633         unsigned n_out;
3634         unsigned o_out;
3635         unsigned n_div;
3636         unsigned o_div;
3637         isl_int gcd;
3638
3639         if (!map)
3640                 return NULL;
3641
3642         hull = isl_map_affine_hull(isl_map_copy(map));
3643         sv = isl_basic_map_plain_is_single_valued(hull);
3644         if (sv >= 0 && sv)
3645                 return plain_pw_multi_aff_from_map(isl_map_domain(map), hull);
3646         if (sv < 0)
3647                 hull = isl_basic_map_free(hull);
3648         if (!hull)
3649                 goto error;
3650
3651         n_div = isl_basic_map_dim(hull, isl_dim_div);
3652         o_div = isl_basic_map_offset(hull, isl_dim_div);
3653
3654         if (n_div == 0) {
3655                 isl_basic_map_free(hull);
3656                 return pw_multi_aff_from_map_check_div(map);
3657         }
3658
3659         isl_int_init(gcd);
3660
3661         n_out = isl_basic_map_dim(hull, isl_dim_out);
3662         o_out = isl_basic_map_offset(hull, isl_dim_out);
3663
3664         for (i = 0; i < n_out; ++i) {
3665                 for (j = 0; j < hull->n_eq; ++j) {
3666                         isl_int *eq = hull->eq[j];
3667                         isl_pw_multi_aff *res;
3668
3669                         if (!isl_int_is_one(eq[o_out + i]) &&
3670                             !isl_int_is_negone(eq[o_out + i]))
3671                                 continue;
3672                         if (isl_seq_first_non_zero(eq + o_out, i) != -1)
3673                                 continue;
3674                         if (isl_seq_first_non_zero(eq + o_out + i + 1,
3675                                                     n_out - (i + 1)) != -1)
3676                                 continue;
3677                         isl_seq_gcd(eq + o_div, n_div, &gcd);
3678                         if (isl_int_is_zero(gcd))
3679                                 continue;
3680                         if (isl_int_is_one(gcd))
3681                                 continue;
3682
3683                         res = pw_multi_aff_from_map_stride(map, hull,
3684                                                                 i, j, gcd);
3685                         isl_int_clear(gcd);
3686                         return res;
3687                 }
3688         }
3689
3690         isl_int_clear(gcd);
3691         isl_basic_map_free(hull);
3692         return pw_multi_aff_from_map_check_div(map);
3693 error:
3694         isl_map_free(map);
3695         return NULL;
3696 }
3697
3698 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_from_set(__isl_take isl_set *set)
3699 {
3700         return isl_pw_multi_aff_from_map(set);
3701 }
3702
3703 /* Return the piecewise affine expression "set ? 1 : 0".
3704  */
3705 __isl_give isl_pw_aff *isl_set_indicator_function(__isl_take isl_set *set)
3706 {
3707         isl_pw_aff *pa;
3708         isl_space *space = isl_set_get_space(set);
3709         isl_local_space *ls = isl_local_space_from_space(space);
3710         isl_aff *zero = isl_aff_zero_on_domain(isl_local_space_copy(ls));
3711         isl_aff *one = isl_aff_zero_on_domain(ls);
3712
3713         one = isl_aff_add_constant_si(one, 1);
3714         pa = isl_pw_aff_alloc(isl_set_copy(set), one);
3715         set = isl_set_complement(set);
3716         pa = isl_pw_aff_add_disjoint(pa, isl_pw_aff_alloc(set, zero));
3717
3718         return pa;
3719 }
3720
3721 /* Plug in "subs" for dimension "type", "pos" of "aff".
3722  *
3723  * Let i be the dimension to replace and let "subs" be of the form
3724  *
3725  *      f/d
3726  *
3727  * and "aff" of the form
3728  *
3729  *      (a i + g)/m
3730  *
3731  * The result is
3732  *
3733  *      (a f + d g')/(m d)
3734  *
3735  * where g' is the result of plugging in "subs" in each of the integer
3736  * divisions in g.
3737  */
3738 __isl_give isl_aff *isl_aff_substitute(__isl_take isl_aff *aff,
3739         enum isl_dim_type type, unsigned pos, __isl_keep isl_aff *subs)
3740 {
3741         isl_ctx *ctx;
3742         isl_int v;
3743
3744         aff = isl_aff_cow(aff);
3745         if (!aff || !subs)
3746                 return isl_aff_free(aff);
3747
3748         ctx = isl_aff_get_ctx(aff);
3749         if (!isl_space_is_equal(aff->ls->dim, subs->ls->dim))
3750                 isl_die(ctx, isl_error_invalid,
3751                         "spaces don't match", return isl_aff_free(aff));
3752         if (isl_local_space_dim(subs->ls, isl_dim_div) != 0)
3753                 isl_die(ctx, isl_error_unsupported,
3754                         "cannot handle divs yet", return isl_aff_free(aff));
3755
3756         aff->ls = isl_local_space_substitute(aff->ls, type, pos, subs);
3757         if (!aff->ls)
3758                 return isl_aff_free(aff);
3759
3760         aff->v = isl_vec_cow(aff->v);
3761         if (!aff->v)
3762                 return isl_aff_free(aff);
3763
3764         pos += isl_local_space_offset(aff->ls, type);
3765
3766         isl_int_init(v);
3767         isl_seq_substitute(aff->v->el, pos, subs->v->el,
3768                             aff->v->size, subs->v->size, v);
3769         isl_int_clear(v);
3770
3771         return aff;
3772 }
3773
3774 /* Plug in "subs" for dimension "type", "pos" in each of the affine
3775  * expressions in "maff".
3776  */
3777 __isl_give isl_multi_aff *isl_multi_aff_substitute(
3778         __isl_take isl_multi_aff *maff, enum isl_dim_type type, unsigned pos,
3779         __isl_keep isl_aff *subs)
3780 {
3781         int i;
3782
3783         maff = isl_multi_aff_cow(maff);
3784         if (!maff || !subs)
3785                 return isl_multi_aff_free(maff);
3786
3787         if (type == isl_dim_in)
3788                 type = isl_dim_set;
3789
3790         for (i = 0; i < maff->n; ++i) {
3791                 maff->p[i] = isl_aff_substitute(maff->p[i], type, pos, subs);
3792                 if (!maff->p[i])
3793                         return isl_multi_aff_free(maff);
3794         }
3795
3796         return maff;
3797 }
3798
3799 /* Plug in "subs" for dimension "type", "pos" of "pma".
3800  *
3801  * pma is of the form
3802  *
3803  *      A_i(v) -> M_i(v)
3804  *
3805  * while subs is of the form
3806  *
3807  *      v' = B_j(v) -> S_j
3808  *
3809  * Each pair i,j such that C_ij = A_i \cap B_i is non-empty
3810  * has a contribution in the result, in particular
3811  *
3812  *      C_ij(S_j) -> M_i(S_j)
3813  *
3814  * Note that plugging in S_j in C_ij may also result in an empty set
3815  * and this contribution should simply be discarded.
3816  */
3817 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_substitute(
3818         __isl_take isl_pw_multi_aff *pma, enum isl_dim_type type, unsigned pos,
3819         __isl_keep isl_pw_aff *subs)
3820 {
3821         int i, j, n;
3822         isl_pw_multi_aff *res;
3823
3824         if (!pma || !subs)
3825                 return isl_pw_multi_aff_free(pma);
3826
3827         n = pma->n * subs->n;
3828         res = isl_pw_multi_aff_alloc_size(isl_space_copy(pma->dim), n);
3829
3830         for (i = 0; i < pma->n; ++i) {
3831                 for (j = 0; j < subs->n; ++j) {
3832                         isl_set *common;
3833                         isl_multi_aff *res_ij;
3834                         int empty;
3835
3836                         common = isl_set_intersect(
3837                                         isl_set_copy(pma->p[i].set),
3838                                         isl_set_copy(subs->p[j].set));
3839                         common = isl_set_substitute(common,
3840                                         type, pos, subs->p[j].aff);
3841                         empty = isl_set_plain_is_empty(common);
3842                         if (empty < 0 || empty) {
3843                                 isl_set_free(common);
3844                                 if (empty < 0)
3845                                         goto error;
3846                                 continue;
3847                         }
3848
3849                         res_ij = isl_multi_aff_substitute(
3850                                         isl_multi_aff_copy(pma->p[i].maff),
3851                                         type, pos, subs->p[j].aff);
3852
3853                         res = isl_pw_multi_aff_add_piece(res, common, res_ij);
3854                 }
3855         }
3856
3857         isl_pw_multi_aff_free(pma);
3858         return res;
3859 error:
3860         isl_pw_multi_aff_free(pma);
3861         isl_pw_multi_aff_free(res);
3862         return NULL;
3863 }
3864
3865 /* Compute the preimage of the affine expression "src" under "ma"
3866  * and put the result in "dst".  If "has_denom" is set (to one),
3867  * then "src" and "dst" have an extra initial denominator.
3868  * "n_div_ma" is the number of existentials in "ma"
3869  * "n_div_bset" is the number of existentials in "src"
3870  * The resulting "dst" (which is assumed to have been allocated by
3871  * the caller) contains coefficients for both sets of existentials,
3872  * first those in "ma" and then those in "src".
3873  * f, c1, c2 and g are temporary objects that have been initialized
3874  * by the caller.
3875  *
3876  * Let src represent the expression
3877  *
3878  *      (a(p) + b x + c(divs))/d
3879  *
3880  * and let ma represent the expressions
3881  *
3882  *      x_i = (r_i(p) + s_i(y) + t_i(divs'))/m_i
3883  *
3884  * We start out with the following expression for dst:
3885  *
3886  *      (a(p) + 0 y + 0 divs' + f \sum_i b_i x_i + c(divs))/d
3887  *
3888  * with the multiplication factor f initially equal to 1.
3889  * For each x_i that we substitute, we multiply the numerator
3890  * (and denominator) of dst by c_1 = m_i and add the numerator
3891  * of the x_i expression multiplied by c_2 = f b_i,
3892  * after removing the common factors of c_1 and c_2.
3893  * The multiplication factor f also needs to be multiplied by c_1
3894  * for the next x_j, j > i.
3895  */
3896 void isl_seq_preimage(isl_int *dst, isl_int *src,
3897         __isl_keep isl_multi_aff *ma, int n_div_ma, int n_div_bset,
3898         isl_int f, isl_int c1, isl_int c2, isl_int g, int has_denom)
3899 {
3900         int i;
3901         int n_param, n_in, n_out;
3902         int o_div_bset;
3903
3904         n_param = isl_multi_aff_dim(ma, isl_dim_param);
3905         n_in = isl_multi_aff_dim(ma, isl_dim_in);
3906         n_out = isl_multi_aff_dim(ma, isl_dim_out);
3907
3908         o_div_bset = has_denom + 1 + n_param + n_in + n_div_ma;
3909
3910         isl_seq_cpy(dst, src, has_denom + 1 + n_param);
3911         isl_seq_clr(dst + has_denom + 1 + n_param, n_in + n_div_ma);
3912         isl_seq_cpy(dst + o_div_bset,
3913                     src + has_denom + 1 + n_param + n_out, n_div_bset);
3914
3915         isl_int_set_si(f, 1);
3916
3917         for (i = 0; i < n_out; ++i) {
3918                 if (isl_int_is_zero(src[has_denom + 1 + n_param + i]))
3919                         continue;
3920                 isl_int_set(c1, ma->p[i]->v->el[0]);
3921                 isl_int_mul(c2, f, src[has_denom + 1 + n_param + i]);
3922                 isl_int_gcd(g, c1, c2);
3923                 isl_int_divexact(c1, c1, g);
3924                 isl_int_divexact(c2, c2, g);
3925
3926                 isl_int_mul(f, f, c1);
3927                 isl_seq_combine(dst + has_denom, c1, dst + has_denom,
3928                                 c2, ma->p[i]->v->el + 1, ma->p[i]->v->size - 1);
3929                 isl_seq_scale(dst + o_div_bset,
3930                                 dst + o_div_bset, c1, n_div_bset);
3931                 if (has_denom)
3932                         isl_int_mul(dst[0], dst[0], c1);
3933         }
3934 }
3935
3936 /* Compute the pullback of "aff" by the function represented by "ma".
3937  * In other words, plug in "ma" in "aff".  The result is an affine expression
3938  * defined over the domain space of "ma".
3939  *
3940  * If "aff" is represented by
3941  *
3942  *      (a(p) + b x + c(divs))/d
3943  *
3944  * and ma is represented by
3945  *
3946  *      x = D(p) + F(y) + G(divs')
3947  *
3948  * then the result is
3949  *
3950  *      (a(p) + b D(p) + b F(y) + b G(divs') + c(divs))/d
3951  *
3952  * The divs in the local space of the input are similarly adjusted
3953  * through a call to isl_local_space_preimage_multi_aff.
3954  */
3955 __isl_give isl_aff *isl_aff_pullback_multi_aff(__isl_take isl_aff *aff,
3956         __isl_take isl_multi_aff *ma)
3957 {
3958         isl_aff *res = NULL;
3959         isl_local_space *ls;
3960         int n_div_aff, n_div_ma;
3961         isl_int f, c1, c2, g;
3962
3963         ma = isl_multi_aff_align_divs(ma);
3964         if (!aff || !ma)
3965                 goto error;
3966
3967         n_div_aff = isl_aff_dim(aff, isl_dim_div);
3968         n_div_ma = ma->n ? isl_aff_dim(ma->p[0], isl_dim_div) : 0;
3969
3970         ls = isl_aff_get_domain_local_space(aff);
3971         ls = isl_local_space_preimage_multi_aff(ls, isl_multi_aff_copy(ma));
3972         res = isl_aff_alloc(ls);
3973         if (!res)
3974                 goto error;
3975
3976         isl_int_init(f);
3977         isl_int_init(c1);
3978         isl_int_init(c2);
3979         isl_int_init(g);
3980
3981         isl_seq_preimage(res->v->el, aff->v->el, ma, n_div_ma, n_div_aff,
3982                         f, c1, c2, g, 1);
3983
3984         isl_int_clear(f);
3985         isl_int_clear(c1);
3986         isl_int_clear(c2);
3987         isl_int_clear(g);
3988
3989         isl_aff_free(aff);
3990         isl_multi_aff_free(ma);
3991         res = isl_aff_normalize(res);
3992         return res;
3993 error:
3994         isl_aff_free(aff);
3995         isl_multi_aff_free(ma);
3996         isl_aff_free(res);
3997         return NULL;
3998 }
3999
4000 /* Compute the pullback of "ma1" by the function represented by "ma2".
4001  * In other words, plug in "ma2" in "ma1".
4002  */
4003 __isl_give isl_multi_aff *isl_multi_aff_pullback_multi_aff(
4004         __isl_take isl_multi_aff *ma1, __isl_take isl_multi_aff *ma2)
4005 {
4006         int i;
4007         isl_space *space = NULL;
4008
4009         ma2 = isl_multi_aff_align_divs(ma2);
4010         ma1 = isl_multi_aff_cow(ma1);
4011         if (!ma1 || !ma2)
4012                 goto error;
4013
4014         space = isl_space_join(isl_multi_aff_get_space(ma2),
4015                                 isl_multi_aff_get_space(ma1));
4016
4017         for (i = 0; i < ma1->n; ++i) {
4018                 ma1->p[i] = isl_aff_pullback_multi_aff(ma1->p[i],
4019                                                     isl_multi_aff_copy(ma2));
4020                 if (!ma1->p[i])
4021                         goto error;
4022         }
4023
4024         ma1 = isl_multi_aff_reset_space(ma1, space);
4025         isl_multi_aff_free(ma2);
4026         return ma1;
4027 error:
4028         isl_space_free(space);
4029         isl_multi_aff_free(ma2);
4030         isl_multi_aff_free(ma1);
4031         return NULL;
4032 }
4033
4034 /* Extend the local space of "dst" to include the divs
4035  * in the local space of "src".
4036  */
4037 __isl_give isl_aff *isl_aff_align_divs(__isl_take isl_aff *dst,
4038         __isl_keep isl_aff *src)
4039 {
4040         isl_ctx *ctx;
4041         int *exp1 = NULL;
4042         int *exp2 = NULL;
4043         isl_mat *div;
4044
4045         if (!src || !dst)
4046                 return isl_aff_free(dst);
4047
4048         ctx = isl_aff_get_ctx(src);
4049         if (!isl_space_is_equal(src->ls->dim, dst->ls->dim))
4050                 isl_die(ctx, isl_error_invalid,
4051                         "spaces don't match", goto error);
4052
4053         if (src->ls->div->n_row == 0)
4054                 return dst;
4055
4056         exp1 = isl_alloc_array(ctx, int, src->ls->div->n_row);
4057         exp2 = isl_alloc_array(ctx, int, dst->ls->div->n_row);
4058         if (!exp1 || !exp2)
4059                 goto error;
4060
4061         div = isl_merge_divs(src->ls->div, dst->ls->div, exp1, exp2);
4062         dst = isl_aff_expand_divs(dst, div, exp2);
4063         free(exp1);
4064         free(exp2);
4065
4066         return dst;
4067 error:
4068         free(exp1);
4069         free(exp2);
4070         return isl_aff_free(dst);
4071 }
4072
4073 /* Adjust the local spaces of the affine expressions in "maff"
4074  * such that they all have the save divs.
4075  */
4076 __isl_give isl_multi_aff *isl_multi_aff_align_divs(
4077         __isl_take isl_multi_aff *maff)
4078 {
4079         int i;
4080
4081         if (!maff)
4082                 return NULL;
4083         if (maff->n == 0)
4084                 return maff;
4085         maff = isl_multi_aff_cow(maff);
4086         if (!maff)
4087                 return NULL;
4088
4089         for (i = 1; i < maff->n; ++i)
4090                 maff->p[0] = isl_aff_align_divs(maff->p[0], maff->p[i]);
4091         for (i = 1; i < maff->n; ++i) {
4092                 maff->p[i] = isl_aff_align_divs(maff->p[i], maff->p[0]);
4093                 if (!maff->p[i])
4094                         return isl_multi_aff_free(maff);
4095         }
4096
4097         return maff;
4098 }
4099
4100 __isl_give isl_aff *isl_aff_lift(__isl_take isl_aff *aff)
4101 {
4102         aff = isl_aff_cow(aff);
4103         if (!aff)
4104                 return NULL;
4105
4106         aff->ls = isl_local_space_lift(aff->ls);
4107         if (!aff->ls)
4108                 return isl_aff_free(aff);
4109
4110         return aff;
4111 }
4112
4113 /* Lift "maff" to a space with extra dimensions such that the result
4114  * has no more existentially quantified variables.
4115  * If "ls" is not NULL, then *ls is assigned the local space that lies
4116  * at the basis of the lifting applied to "maff".
4117  */
4118 __isl_give isl_multi_aff *isl_multi_aff_lift(__isl_take isl_multi_aff *maff,
4119         __isl_give isl_local_space **ls)
4120 {
4121         int i;
4122         isl_space *space;
4123         unsigned n_div;
4124
4125         if (ls)
4126                 *ls = NULL;
4127
4128         if (!maff)
4129                 return NULL;
4130
4131         if (maff->n == 0) {
4132                 if (ls) {
4133                         isl_space *space = isl_multi_aff_get_domain_space(maff);
4134                         *ls = isl_local_space_from_space(space);
4135                         if (!*ls)
4136                                 return isl_multi_aff_free(maff);
4137                 }
4138                 return maff;
4139         }
4140
4141         maff = isl_multi_aff_cow(maff);
4142         maff = isl_multi_aff_align_divs(maff);
4143         if (!maff)
4144                 return NULL;
4145
4146         n_div = isl_aff_dim(maff->p[0], isl_dim_div);
4147         space = isl_multi_aff_get_space(maff);
4148         space = isl_space_lift(isl_space_domain(space), n_div);
4149         space = isl_space_extend_domain_with_range(space,
4150                                                 isl_multi_aff_get_space(maff));
4151         if (!space)
4152                 return isl_multi_aff_free(maff);
4153         isl_space_free(maff->space);
4154         maff->space = space;
4155
4156         if (ls) {
4157                 *ls = isl_aff_get_domain_local_space(maff->p[0]);
4158                 if (!*ls)
4159                         return isl_multi_aff_free(maff);
4160         }
4161
4162         for (i = 0; i < maff->n; ++i) {
4163                 maff->p[i] = isl_aff_lift(maff->p[i]);
4164                 if (!maff->p[i])
4165                         goto error;
4166         }
4167
4168         return maff;
4169 error:
4170         if (ls)
4171                 isl_local_space_free(*ls);
4172         return isl_multi_aff_free(maff);
4173 }
4174
4175
4176 /* Extract an isl_pw_aff corresponding to output dimension "pos" of "pma".
4177  */
4178 __isl_give isl_pw_aff *isl_pw_multi_aff_get_pw_aff(
4179         __isl_keep isl_pw_multi_aff *pma, int pos)
4180 {
4181         int i;
4182         int n_out;
4183         isl_space *space;
4184         isl_pw_aff *pa;
4185
4186         if (!pma)
4187                 return NULL;
4188
4189         n_out = isl_pw_multi_aff_dim(pma, isl_dim_out);
4190         if (pos < 0 || pos >= n_out)
4191                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
4192                         "index out of bounds", return NULL);
4193
4194         space = isl_pw_multi_aff_get_space(pma);
4195         space = isl_space_drop_dims(space, isl_dim_out,
4196                                     pos + 1, n_out - pos - 1);
4197         space = isl_space_drop_dims(space, isl_dim_out, 0, pos);
4198
4199         pa = isl_pw_aff_alloc_size(space, pma->n);
4200         for (i = 0; i < pma->n; ++i) {
4201                 isl_aff *aff;
4202                 aff = isl_multi_aff_get_aff(pma->p[i].maff, pos);
4203                 pa = isl_pw_aff_add_piece(pa, isl_set_copy(pma->p[i].set), aff);
4204         }
4205
4206         return pa;
4207 }
4208
4209 /* Return an isl_pw_multi_aff with the given "set" as domain and
4210  * an unnamed zero-dimensional range.
4211  */
4212 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_from_domain(
4213         __isl_take isl_set *set)
4214 {
4215         isl_multi_aff *ma;
4216         isl_space *space;
4217
4218         space = isl_set_get_space(set);
4219         space = isl_space_from_domain(space);
4220         ma = isl_multi_aff_zero(space);
4221         return isl_pw_multi_aff_alloc(set, ma);
4222 }
4223
4224 /* Add an isl_pw_multi_aff with the given "set" as domain and
4225  * an unnamed zero-dimensional range to *user.
4226  */
4227 static int add_pw_multi_aff_from_domain(__isl_take isl_set *set, void *user)
4228 {
4229         isl_union_pw_multi_aff **upma = user;
4230         isl_pw_multi_aff *pma;
4231
4232         pma = isl_pw_multi_aff_from_domain(set);
4233         *upma = isl_union_pw_multi_aff_add_pw_multi_aff(*upma, pma);
4234
4235         return 0;
4236 }
4237
4238 /* Return an isl_union_pw_multi_aff with the given "uset" as domain and
4239  * an unnamed zero-dimensional range.
4240  */
4241 __isl_give isl_union_pw_multi_aff *isl_union_pw_multi_aff_from_domain(
4242         __isl_take isl_union_set *uset)
4243 {
4244         isl_space *space;
4245         isl_union_pw_multi_aff *upma;
4246
4247         if (!uset)
4248                 return NULL;
4249
4250         space = isl_union_set_get_space(uset);
4251         upma = isl_union_pw_multi_aff_empty(space);
4252
4253         if (isl_union_set_foreach_set(uset,
4254                                     &add_pw_multi_aff_from_domain, &upma) < 0)
4255                 goto error;
4256
4257         isl_union_set_free(uset);
4258         return upma;
4259 error:
4260         isl_union_set_free(uset);
4261         isl_union_pw_multi_aff_free(upma);
4262         return NULL;
4263 }
4264
4265 /* Convert "pma" to an isl_map and add it to *umap.
4266  */
4267 static int map_from_pw_multi_aff(__isl_take isl_pw_multi_aff *pma, void *user)
4268 {
4269         isl_union_map **umap = user;
4270         isl_map *map;
4271
4272         map = isl_map_from_pw_multi_aff(pma);
4273         *umap = isl_union_map_add_map(*umap, map);
4274
4275         return 0;
4276 }
4277
4278 /* Construct a union map mapping the domain of the union
4279  * piecewise multi-affine expression to its range, with each dimension
4280  * in the range equated to the corresponding affine expression on its cell.
4281  */
4282 __isl_give isl_union_map *isl_union_map_from_union_pw_multi_aff(
4283         __isl_take isl_union_pw_multi_aff *upma)
4284 {
4285         isl_space *space;
4286         isl_union_map *umap;
4287
4288         if (!upma)
4289                 return NULL;
4290
4291         space = isl_union_pw_multi_aff_get_space(upma);
4292         umap = isl_union_map_empty(space);
4293
4294         if (isl_union_pw_multi_aff_foreach_pw_multi_aff(upma,
4295                                         &map_from_pw_multi_aff, &umap) < 0)
4296                 goto error;
4297
4298         isl_union_pw_multi_aff_free(upma);
4299         return umap;
4300 error:
4301         isl_union_pw_multi_aff_free(upma);
4302         isl_union_map_free(umap);
4303         return NULL;
4304 }
4305
4306 /* Local data for bin_entry and the callback "fn".
4307  */
4308 struct isl_union_pw_multi_aff_bin_data {
4309         isl_union_pw_multi_aff *upma2;
4310         isl_union_pw_multi_aff *res;
4311         isl_pw_multi_aff *pma;
4312         int (*fn)(void **entry, void *user);
4313 };
4314
4315 /* Given an isl_pw_multi_aff from upma1, store it in data->pma
4316  * and call data->fn for each isl_pw_multi_aff in data->upma2.
4317  */
4318 static int bin_entry(void **entry, void *user)
4319 {
4320         struct isl_union_pw_multi_aff_bin_data *data = user;
4321         isl_pw_multi_aff *pma = *entry;
4322
4323         data->pma = pma;
4324         if (isl_hash_table_foreach(data->upma2->dim->ctx, &data->upma2->table,
4325                                    data->fn, data) < 0)
4326                 return -1;
4327
4328         return 0;
4329 }
4330
4331 /* Call "fn" on each pair of isl_pw_multi_affs in "upma1" and "upma2".
4332  * The isl_pw_multi_aff from upma1 is stored in data->pma (where data is
4333  * passed as user field) and the isl_pw_multi_aff from upma2 is available
4334  * as *entry.  The callback should adjust data->res if desired.
4335  */
4336 static __isl_give isl_union_pw_multi_aff *bin_op(
4337         __isl_take isl_union_pw_multi_aff *upma1,
4338         __isl_take isl_union_pw_multi_aff *upma2,
4339         int (*fn)(void **entry, void *user))
4340 {
4341         isl_space *space;
4342         struct isl_union_pw_multi_aff_bin_data data = { NULL, NULL, NULL, fn };
4343
4344         space = isl_union_pw_multi_aff_get_space(upma2);
4345         upma1 = isl_union_pw_multi_aff_align_params(upma1, space);
4346         space = isl_union_pw_multi_aff_get_space(upma1);
4347         upma2 = isl_union_pw_multi_aff_align_params(upma2, space);
4348
4349         if (!upma1 || !upma2)
4350                 goto error;
4351
4352         data.upma2 = upma2;
4353         data.res = isl_union_pw_multi_aff_alloc(isl_space_copy(upma1->dim),
4354                                        upma1->table.n);
4355         if (isl_hash_table_foreach(upma1->dim->ctx, &upma1->table,
4356                                    &bin_entry, &data) < 0)
4357                 goto error;
4358
4359         isl_union_pw_multi_aff_free(upma1);
4360         isl_union_pw_multi_aff_free(upma2);
4361         return data.res;
4362 error:
4363         isl_union_pw_multi_aff_free(upma1);
4364         isl_union_pw_multi_aff_free(upma2);
4365         isl_union_pw_multi_aff_free(data.res);
4366         return NULL;
4367 }
4368
4369 /* Given two aligned isl_pw_multi_affs A -> B and C -> D,
4370  * construct an isl_pw_multi_aff (A * C) -> [B -> D].
4371  */
4372 static __isl_give isl_pw_multi_aff *pw_multi_aff_range_product(
4373         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
4374 {
4375         isl_space *space;
4376
4377         space = isl_space_range_product(isl_pw_multi_aff_get_space(pma1),
4378                                         isl_pw_multi_aff_get_space(pma2));
4379         return isl_pw_multi_aff_on_shared_domain_in(pma1, pma2, space,
4380                                             &isl_multi_aff_range_product);
4381 }
4382
4383 /* Given two isl_pw_multi_affs A -> B and C -> D,
4384  * construct an isl_pw_multi_aff (A * C) -> [B -> D].
4385  */
4386 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_range_product(
4387         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
4388 {
4389         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
4390                                             &pw_multi_aff_range_product);
4391 }
4392
4393 /* Given two aligned isl_pw_multi_affs A -> B and C -> D,
4394  * construct an isl_pw_multi_aff (A * C) -> (B, D).
4395  */
4396 static __isl_give isl_pw_multi_aff *pw_multi_aff_flat_range_product(
4397         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
4398 {
4399         isl_space *space;
4400
4401         space = isl_space_range_product(isl_pw_multi_aff_get_space(pma1),
4402                                         isl_pw_multi_aff_get_space(pma2));
4403         space = isl_space_flatten_range(space);
4404         return isl_pw_multi_aff_on_shared_domain_in(pma1, pma2, space,
4405                                             &isl_multi_aff_flat_range_product);
4406 }
4407
4408 /* Given two isl_pw_multi_affs A -> B and C -> D,
4409  * construct an isl_pw_multi_aff (A * C) -> (B, D).
4410  */
4411 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_flat_range_product(
4412         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
4413 {
4414         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
4415                                             &pw_multi_aff_flat_range_product);
4416 }
4417
4418 /* If data->pma and *entry have the same domain space, then compute
4419  * their flat range product and the result to data->res.
4420  */
4421 static int flat_range_product_entry(void **entry, void *user)
4422 {
4423         struct isl_union_pw_multi_aff_bin_data *data = user;
4424         isl_pw_multi_aff *pma2 = *entry;
4425
4426         if (!isl_space_tuple_match(data->pma->dim, isl_dim_in,
4427                                  pma2->dim, isl_dim_in))
4428                 return 0;
4429
4430         pma2 = isl_pw_multi_aff_flat_range_product(
4431                                         isl_pw_multi_aff_copy(data->pma),
4432                                         isl_pw_multi_aff_copy(pma2));
4433
4434         data->res = isl_union_pw_multi_aff_add_pw_multi_aff(data->res, pma2);
4435
4436         return 0;
4437 }
4438
4439 /* Given two isl_union_pw_multi_affs A -> B and C -> D,
4440  * construct an isl_union_pw_multi_aff (A * C) -> (B, D).
4441  */
4442 __isl_give isl_union_pw_multi_aff *isl_union_pw_multi_aff_flat_range_product(
4443         __isl_take isl_union_pw_multi_aff *upma1,
4444         __isl_take isl_union_pw_multi_aff *upma2)
4445 {
4446         return bin_op(upma1, upma2, &flat_range_product_entry);
4447 }
4448
4449 /* Replace the affine expressions at position "pos" in "pma" by "pa".
4450  * The parameters are assumed to have been aligned.
4451  *
4452  * The implementation essentially performs an isl_pw_*_on_shared_domain,
4453  * except that it works on two different isl_pw_* types.
4454  */
4455 static __isl_give isl_pw_multi_aff *pw_multi_aff_set_pw_aff(
4456         __isl_take isl_pw_multi_aff *pma, unsigned pos,
4457         __isl_take isl_pw_aff *pa)
4458 {
4459         int i, j, n;
4460         isl_pw_multi_aff *res = NULL;
4461
4462         if (!pma || !pa)
4463                 goto error;
4464
4465         if (!isl_space_tuple_match(pma->dim, isl_dim_in, pa->dim, isl_dim_in))
4466                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
4467                         "domains don't match", goto error);
4468         if (pos >= isl_pw_multi_aff_dim(pma, isl_dim_out))
4469                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
4470                         "index out of bounds", goto error);
4471
4472         n = pma->n * pa->n;
4473         res = isl_pw_multi_aff_alloc_size(isl_pw_multi_aff_get_space(pma), n);
4474
4475         for (i = 0; i < pma->n; ++i) {
4476                 for (j = 0; j < pa->n; ++j) {
4477                         isl_set *common;
4478                         isl_multi_aff *res_ij;
4479                         int empty;
4480
4481                         common = isl_set_intersect(isl_set_copy(pma->p[i].set),
4482                                                    isl_set_copy(pa->p[j].set));
4483                         empty = isl_set_plain_is_empty(common);
4484                         if (empty < 0 || empty) {
4485                                 isl_set_free(common);
4486                                 if (empty < 0)
4487                                         goto error;
4488                                 continue;
4489                         }
4490
4491                         res_ij = isl_multi_aff_set_aff(
4492                                         isl_multi_aff_copy(pma->p[i].maff), pos,
4493                                         isl_aff_copy(pa->p[j].aff));
4494                         res_ij = isl_multi_aff_gist(res_ij,
4495                                         isl_set_copy(common));
4496
4497                         res = isl_pw_multi_aff_add_piece(res, common, res_ij);
4498                 }
4499         }
4500
4501         isl_pw_multi_aff_free(pma);
4502         isl_pw_aff_free(pa);
4503         return res;
4504 error:
4505         isl_pw_multi_aff_free(pma);
4506         isl_pw_aff_free(pa);
4507         return isl_pw_multi_aff_free(res);
4508 }
4509
4510 /* Replace the affine expressions at position "pos" in "pma" by "pa".
4511  */
4512 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_set_pw_aff(
4513         __isl_take isl_pw_multi_aff *pma, unsigned pos,
4514         __isl_take isl_pw_aff *pa)
4515 {
4516         if (!pma || !pa)
4517                 goto error;
4518         if (isl_space_match(pma->dim, isl_dim_param, pa->dim, isl_dim_param))
4519                 return pw_multi_aff_set_pw_aff(pma, pos, pa);
4520         if (!isl_space_has_named_params(pma->dim) ||
4521             !isl_space_has_named_params(pa->dim))
4522                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
4523                         "unaligned unnamed parameters", goto error);
4524         pma = isl_pw_multi_aff_align_params(pma, isl_pw_aff_get_space(pa));
4525         pa = isl_pw_aff_align_params(pa, isl_pw_multi_aff_get_space(pma));
4526         return pw_multi_aff_set_pw_aff(pma, pos, pa);
4527 error:
4528         isl_pw_multi_aff_free(pma);
4529         isl_pw_aff_free(pa);
4530         return NULL;
4531 }
4532
4533 #undef BASE
4534 #define BASE pw_aff
4535
4536 #include <isl_multi_templ.c>