isl_aff.c

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