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 /* Check that the domain space of "aff" matches "space".
2663  *
2664  * Return 0 on success and -1 on error.
2665  */
2666 int isl_aff_check_match_domain_space(__isl_keep isl_aff *aff,
2667         __isl_keep isl_space *space)
2668 {
2669         isl_space *aff_space;
2670         int match;
2671
2672         if (!aff || !space)
2673                 return -1;
2674
2675         aff_space = isl_aff_get_domain_space(aff);
2676
2677         match = isl_space_match(space, isl_dim_param, aff_space, isl_dim_param);
2678         if (match < 0)
2679                 goto error;
2680         if (!match)
2681                 isl_die(isl_aff_get_ctx(aff), isl_error_invalid,
2682                         "parameters don't match", goto error);
2683         match = isl_space_tuple_match(space, isl_dim_in,
2684                                         aff_space, isl_dim_set);
2685         if (match < 0)
2686                 goto error;
2687         if (!match)
2688                 isl_die(isl_aff_get_ctx(aff), isl_error_invalid,
2689                         "domains don't match", goto error);
2690         isl_space_free(aff_space);
2691         return 0;
2692 error:
2693         isl_space_free(aff_space);
2694         return -1;
2695 }
2696
2697 #undef BASE
2698 #define BASE aff
2699
2700 #include <isl_multi_templ.c>
2701
2702 /* Create an isl_pw_multi_aff with the given isl_multi_aff on a universe
2703  * domain.
2704  */
2705 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_from_multi_aff(
2706         __isl_take isl_multi_aff *ma)
2707 {
2708         isl_set *dom = isl_set_universe(isl_multi_aff_get_domain_space(ma));
2709         return isl_pw_multi_aff_alloc(dom, ma);
2710 }
2711
2712 /* Create a piecewise multi-affine expression in the given space that maps each
2713  * input dimension to the corresponding output dimension.
2714  */
2715 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_identity(
2716         __isl_take isl_space *space)
2717 {
2718         return isl_pw_multi_aff_from_multi_aff(isl_multi_aff_identity(space));
2719 }
2720
2721 __isl_give isl_multi_aff *isl_multi_aff_add(__isl_take isl_multi_aff *maff1,
2722         __isl_take isl_multi_aff *maff2)
2723 {
2724         return isl_multi_aff_bin_op(maff1, maff2, &isl_aff_add);
2725 }
2726
2727 /* Subtract "ma2" from "ma1" and return the result.
2728  */
2729 __isl_give isl_multi_aff *isl_multi_aff_sub(__isl_take isl_multi_aff *ma1,
2730         __isl_take isl_multi_aff *ma2)
2731 {
2732         return isl_multi_aff_bin_op(ma1, ma2, &isl_aff_sub);
2733 }
2734
2735 /* Given two multi-affine expressions A -> B and C -> D,
2736  * construct a multi-affine expression [A -> C] -> [B -> D].
2737  */
2738 __isl_give isl_multi_aff *isl_multi_aff_product(
2739         __isl_take isl_multi_aff *ma1, __isl_take isl_multi_aff *ma2)
2740 {
2741         int i;
2742         isl_aff *aff;
2743         isl_space *space;
2744         isl_multi_aff *res;
2745         int in1, in2, out1, out2;
2746
2747         in1 = isl_multi_aff_dim(ma1, isl_dim_in);
2748         in2 = isl_multi_aff_dim(ma2, isl_dim_in);
2749         out1 = isl_multi_aff_dim(ma1, isl_dim_out);
2750         out2 = isl_multi_aff_dim(ma2, isl_dim_out);
2751         space = isl_space_product(isl_multi_aff_get_space(ma1),
2752                                   isl_multi_aff_get_space(ma2));
2753         res = isl_multi_aff_alloc(isl_space_copy(space));
2754         space = isl_space_domain(space);
2755
2756         for (i = 0; i < out1; ++i) {
2757                 aff = isl_multi_aff_get_aff(ma1, i);
2758                 aff = isl_aff_insert_dims(aff, isl_dim_in, in1, in2);
2759                 aff = isl_aff_reset_domain_space(aff, isl_space_copy(space));
2760                 res = isl_multi_aff_set_aff(res, i, aff);
2761         }
2762
2763         for (i = 0; i < out2; ++i) {
2764                 aff = isl_multi_aff_get_aff(ma2, i);
2765                 aff = isl_aff_insert_dims(aff, isl_dim_in, 0, in1);
2766                 aff = isl_aff_reset_domain_space(aff, isl_space_copy(space));
2767                 res = isl_multi_aff_set_aff(res, out1 + i, aff);
2768         }
2769
2770         isl_space_free(space);
2771         isl_multi_aff_free(ma1);
2772         isl_multi_aff_free(ma2);
2773         return res;
2774 }
2775
2776 /* Exploit the equalities in "eq" to simplify the affine expressions.
2777  */
2778 static __isl_give isl_multi_aff *isl_multi_aff_substitute_equalities(
2779         __isl_take isl_multi_aff *maff, __isl_take isl_basic_set *eq)
2780 {
2781         int i;
2782
2783         maff = isl_multi_aff_cow(maff);
2784         if (!maff || !eq)
2785                 goto error;
2786
2787         for (i = 0; i < maff->n; ++i) {
2788                 maff->p[i] = isl_aff_substitute_equalities(maff->p[i],
2789                                                     isl_basic_set_copy(eq));
2790                 if (!maff->p[i])
2791                         goto error;
2792         }
2793
2794         isl_basic_set_free(eq);
2795         return maff;
2796 error:
2797         isl_basic_set_free(eq);
2798         isl_multi_aff_free(maff);
2799         return NULL;
2800 }
2801
2802 __isl_give isl_multi_aff *isl_multi_aff_scale(__isl_take isl_multi_aff *maff,
2803         isl_int f)
2804 {
2805         int i;
2806
2807         maff = isl_multi_aff_cow(maff);
2808         if (!maff)
2809                 return NULL;
2810
2811         for (i = 0; i < maff->n; ++i) {
2812                 maff->p[i] = isl_aff_scale(maff->p[i], f);
2813                 if (!maff->p[i])
2814                         return isl_multi_aff_free(maff);
2815         }
2816
2817         return maff;
2818 }
2819
2820 __isl_give isl_multi_aff *isl_multi_aff_add_on_domain(__isl_keep isl_set *dom,
2821         __isl_take isl_multi_aff *maff1, __isl_take isl_multi_aff *maff2)
2822 {
2823         maff1 = isl_multi_aff_add(maff1, maff2);
2824         maff1 = isl_multi_aff_gist(maff1, isl_set_copy(dom));
2825         return maff1;
2826 }
2827
2828 int isl_multi_aff_is_empty(__isl_keep isl_multi_aff *maff)
2829 {
2830         if (!maff)
2831                 return -1;
2832
2833         return 0;
2834 }
2835
2836 int isl_multi_aff_plain_is_equal(__isl_keep isl_multi_aff *maff1,
2837         __isl_keep isl_multi_aff *maff2)
2838 {
2839         int i;
2840         int equal;
2841
2842         if (!maff1 || !maff2)
2843                 return -1;
2844         if (maff1->n != maff2->n)
2845                 return 0;
2846         equal = isl_space_is_equal(maff1->space, maff2->space);
2847         if (equal < 0 || !equal)
2848                 return equal;
2849
2850         for (i = 0; i < maff1->n; ++i) {
2851                 equal = isl_aff_plain_is_equal(maff1->p[i], maff2->p[i]);
2852                 if (equal < 0 || !equal)
2853                         return equal;
2854         }
2855
2856         return 1;
2857 }
2858
2859 /* Return the set of domain elements where "ma1" is lexicographically
2860  * smaller than or equal to "ma2".
2861  */
2862 __isl_give isl_set *isl_multi_aff_lex_le_set(__isl_take isl_multi_aff *ma1,
2863         __isl_take isl_multi_aff *ma2)
2864 {
2865         return isl_multi_aff_lex_ge_set(ma2, ma1);
2866 }
2867
2868 /* Return the set of domain elements where "ma1" is lexicographically
2869  * greater than or equal to "ma2".
2870  */
2871 __isl_give isl_set *isl_multi_aff_lex_ge_set(__isl_take isl_multi_aff *ma1,
2872         __isl_take isl_multi_aff *ma2)
2873 {
2874         isl_space *space;
2875         isl_map *map1, *map2;
2876         isl_map *map, *ge;
2877
2878         map1 = isl_map_from_multi_aff(ma1);
2879         map2 = isl_map_from_multi_aff(ma2);
2880         map = isl_map_range_product(map1, map2);
2881         space = isl_space_range(isl_map_get_space(map));
2882         space = isl_space_domain(isl_space_unwrap(space));
2883         ge = isl_map_lex_ge(space);
2884         map = isl_map_intersect_range(map, isl_map_wrap(ge));
2885
2886         return isl_map_domain(map);
2887 }
2888
2889 #undef PW
2890 #define PW isl_pw_multi_aff
2891 #undef EL
2892 #define EL isl_multi_aff
2893 #undef EL_IS_ZERO
2894 #define EL_IS_ZERO is_empty
2895 #undef ZERO
2896 #define ZERO empty
2897 #undef IS_ZERO
2898 #define IS_ZERO is_empty
2899 #undef FIELD
2900 #define FIELD maff
2901 #undef DEFAULT_IS_ZERO
2902 #define DEFAULT_IS_ZERO 0
2903
2904 #define NO_NEG
2905 #define NO_EVAL
2906 #define NO_OPT
2907 #define NO_INVOLVES_DIMS
2908 #define NO_MOVE_DIMS
2909 #define NO_INSERT_DIMS
2910 #define NO_LIFT
2911 #define NO_MORPH
2912
2913 #include <isl_pw_templ.c>
2914
2915 #undef UNION
2916 #define UNION isl_union_pw_multi_aff
2917 #undef PART
2918 #define PART isl_pw_multi_aff
2919 #undef PARTS
2920 #define PARTS pw_multi_aff
2921 #define ALIGN_DOMAIN
2922
2923 #define NO_EVAL
2924
2925 #include <isl_union_templ.c>
2926
2927 /* Given a function "cmp" that returns the set of elements where
2928  * "ma1" is "better" than "ma2", return the intersection of this
2929  * set with "dom1" and "dom2".
2930  */
2931 static __isl_give isl_set *shared_and_better(__isl_keep isl_set *dom1,
2932         __isl_keep isl_set *dom2, __isl_keep isl_multi_aff *ma1,
2933         __isl_keep isl_multi_aff *ma2,
2934         __isl_give isl_set *(*cmp)(__isl_take isl_multi_aff *ma1,
2935                                     __isl_take isl_multi_aff *ma2))
2936 {
2937         isl_set *common;
2938         isl_set *better;
2939         int is_empty;
2940
2941         common = isl_set_intersect(isl_set_copy(dom1), isl_set_copy(dom2));
2942         is_empty = isl_set_plain_is_empty(common);
2943         if (is_empty >= 0 && is_empty)
2944                 return common;
2945         if (is_empty < 0)
2946                 return isl_set_free(common);
2947         better = cmp(isl_multi_aff_copy(ma1), isl_multi_aff_copy(ma2));
2948         better = isl_set_intersect(common, better);
2949
2950         return better;
2951 }
2952
2953 /* Given a function "cmp" that returns the set of elements where
2954  * "ma1" is "better" than "ma2", return a piecewise multi affine
2955  * expression defined on the union of the definition domains
2956  * of "pma1" and "pma2" that maps to the "best" of "pma1" and
2957  * "pma2" on each cell.  If only one of the two input functions
2958  * is defined on a given cell, then it is considered the best.
2959  */
2960 static __isl_give isl_pw_multi_aff *pw_multi_aff_union_opt(
2961         __isl_take isl_pw_multi_aff *pma1,
2962         __isl_take isl_pw_multi_aff *pma2,
2963         __isl_give isl_set *(*cmp)(__isl_take isl_multi_aff *ma1,
2964                                     __isl_take isl_multi_aff *ma2))
2965 {
2966         int i, j, n;
2967         isl_pw_multi_aff *res = NULL;
2968         isl_ctx *ctx;
2969         isl_set *set = NULL;
2970
2971         if (!pma1 || !pma2)
2972                 goto error;
2973
2974         ctx = isl_space_get_ctx(pma1->dim);
2975         if (!isl_space_is_equal(pma1->dim, pma2->dim))
2976                 isl_die(ctx, isl_error_invalid,
2977                         "arguments should live in the same space", goto error);
2978
2979         if (isl_pw_multi_aff_is_empty(pma1)) {
2980                 isl_pw_multi_aff_free(pma1);
2981                 return pma2;
2982         }
2983
2984         if (isl_pw_multi_aff_is_empty(pma2)) {
2985                 isl_pw_multi_aff_free(pma2);
2986                 return pma1;
2987         }
2988
2989         n = 2 * (pma1->n + 1) * (pma2->n + 1);
2990         res = isl_pw_multi_aff_alloc_size(isl_space_copy(pma1->dim), n);
2991
2992         for (i = 0; i < pma1->n; ++i) {
2993                 set = isl_set_copy(pma1->p[i].set);
2994                 for (j = 0; j < pma2->n; ++j) {
2995                         isl_set *better;
2996                         int is_empty;
2997
2998                         better = shared_and_better(pma2->p[j].set,
2999                                         pma1->p[i].set, pma2->p[j].maff,
3000                                         pma1->p[i].maff, cmp);
3001                         is_empty = isl_set_plain_is_empty(better);
3002                         if (is_empty < 0 || is_empty) {
3003                                 isl_set_free(better);
3004                                 if (is_empty < 0)
3005                                         goto error;
3006                                 continue;
3007                         }
3008                         set = isl_set_subtract(set, isl_set_copy(better));
3009
3010                         res = isl_pw_multi_aff_add_piece(res, better,
3011                                         isl_multi_aff_copy(pma2->p[j].maff));
3012                 }
3013                 res = isl_pw_multi_aff_add_piece(res, set,
3014                                         isl_multi_aff_copy(pma1->p[i].maff));
3015         }
3016
3017         for (j = 0; j < pma2->n; ++j) {
3018                 set = isl_set_copy(pma2->p[j].set);
3019                 for (i = 0; i < pma1->n; ++i)
3020                         set = isl_set_subtract(set,
3021                                         isl_set_copy(pma1->p[i].set));
3022                 res = isl_pw_multi_aff_add_piece(res, set,
3023                                         isl_multi_aff_copy(pma2->p[j].maff));
3024         }
3025
3026         isl_pw_multi_aff_free(pma1);
3027         isl_pw_multi_aff_free(pma2);
3028
3029         return res;
3030 error:
3031         isl_pw_multi_aff_free(pma1);
3032         isl_pw_multi_aff_free(pma2);
3033         isl_set_free(set);
3034         return isl_pw_multi_aff_free(res);
3035 }
3036
3037 static __isl_give isl_pw_multi_aff *pw_multi_aff_union_lexmax(
3038         __isl_take isl_pw_multi_aff *pma1,
3039         __isl_take isl_pw_multi_aff *pma2)
3040 {
3041         return pw_multi_aff_union_opt(pma1, pma2, &isl_multi_aff_lex_ge_set);
3042 }
3043
3044 /* Given two piecewise multi affine expressions, return a piecewise
3045  * multi-affine expression defined on the union of the definition domains
3046  * of the inputs that is equal to the lexicographic maximum of the two
3047  * inputs on each cell.  If only one of the two inputs is defined on
3048  * a given cell, then it is considered to be the maximum.
3049  */
3050 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_union_lexmax(
3051         __isl_take isl_pw_multi_aff *pma1,
3052         __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_union_lexmax);
3056 }
3057
3058 static __isl_give isl_pw_multi_aff *pw_multi_aff_union_lexmin(
3059         __isl_take isl_pw_multi_aff *pma1,
3060         __isl_take isl_pw_multi_aff *pma2)
3061 {
3062         return pw_multi_aff_union_opt(pma1, pma2, &isl_multi_aff_lex_le_set);
3063 }
3064
3065 /* Given two piecewise multi affine expressions, return a piecewise
3066  * multi-affine expression defined on the union of the definition domains
3067  * of the inputs that is equal to the lexicographic minimum of the two
3068  * inputs on each cell.  If only one of the two inputs is defined on
3069  * a given cell, then it is considered to be the minimum.
3070  */
3071 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_union_lexmin(
3072         __isl_take isl_pw_multi_aff *pma1,
3073         __isl_take isl_pw_multi_aff *pma2)
3074 {
3075         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
3076                                                     &pw_multi_aff_union_lexmin);
3077 }
3078
3079 static __isl_give isl_pw_multi_aff *pw_multi_aff_add(
3080         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
3081 {
3082         return isl_pw_multi_aff_on_shared_domain(pma1, pma2,
3083                                                 &isl_multi_aff_add);
3084 }
3085
3086 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_add(
3087         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
3088 {
3089         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
3090                                                 &pw_multi_aff_add);
3091 }
3092
3093 static __isl_give isl_pw_multi_aff *pw_multi_aff_sub(
3094         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
3095 {
3096         return isl_pw_multi_aff_on_shared_domain(pma1, pma2,
3097                                                 &isl_multi_aff_sub);
3098 }
3099
3100 /* Subtract "pma2" from "pma1" and return the result.
3101  */
3102 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_sub(
3103         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
3104 {
3105         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
3106                                                 &pw_multi_aff_sub);
3107 }
3108
3109 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_union_add(
3110         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
3111 {
3112         return isl_pw_multi_aff_union_add_(pma1, pma2);
3113 }
3114
3115 /* Given two piecewise multi-affine expressions A -> B and C -> D,
3116  * construct a piecewise multi-affine expression [A -> C] -> [B -> D].
3117  */
3118 static __isl_give isl_pw_multi_aff *pw_multi_aff_product(
3119         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
3120 {
3121         int i, j, n;
3122         isl_space *space;
3123         isl_pw_multi_aff *res;
3124
3125         if (!pma1 || !pma2)
3126                 goto error;
3127
3128         n = pma1->n * pma2->n;
3129         space = isl_space_product(isl_space_copy(pma1->dim),
3130                                   isl_space_copy(pma2->dim));
3131         res = isl_pw_multi_aff_alloc_size(space, n);
3132
3133         for (i = 0; i < pma1->n; ++i) {
3134                 for (j = 0; j < pma2->n; ++j) {
3135                         isl_set *domain;
3136                         isl_multi_aff *ma;
3137
3138                         domain = isl_set_product(isl_set_copy(pma1->p[i].set),
3139                                                  isl_set_copy(pma2->p[j].set));
3140                         ma = isl_multi_aff_product(
3141                                         isl_multi_aff_copy(pma1->p[i].maff),
3142                                         isl_multi_aff_copy(pma2->p[i].maff));
3143                         res = isl_pw_multi_aff_add_piece(res, domain, ma);
3144                 }
3145         }
3146
3147         isl_pw_multi_aff_free(pma1);
3148         isl_pw_multi_aff_free(pma2);
3149         return res;
3150 error:
3151         isl_pw_multi_aff_free(pma1);
3152         isl_pw_multi_aff_free(pma2);
3153         return NULL;
3154 }
3155
3156 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_product(
3157         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
3158 {
3159         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
3160                                                 &pw_multi_aff_product);
3161 }
3162
3163 /* Construct a map mapping the domain of the piecewise multi-affine expression
3164  * to its range, with each dimension in the range equated to the
3165  * corresponding affine expression on its cell.
3166  */
3167 __isl_give isl_map *isl_map_from_pw_multi_aff(__isl_take isl_pw_multi_aff *pma)
3168 {
3169         int i;
3170         isl_map *map;
3171
3172         if (!pma)
3173                 return NULL;
3174
3175         map = isl_map_empty(isl_pw_multi_aff_get_space(pma));
3176
3177         for (i = 0; i < pma->n; ++i) {
3178                 isl_multi_aff *maff;
3179                 isl_basic_map *bmap;
3180                 isl_map *map_i;
3181
3182                 maff = isl_multi_aff_copy(pma->p[i].maff);
3183                 bmap = isl_basic_map_from_multi_aff(maff);
3184                 map_i = isl_map_from_basic_map(bmap);
3185                 map_i = isl_map_intersect_domain(map_i,
3186                                                 isl_set_copy(pma->p[i].set));
3187                 map = isl_map_union_disjoint(map, map_i);
3188         }
3189
3190         isl_pw_multi_aff_free(pma);
3191         return map;
3192 }
3193
3194 __isl_give isl_set *isl_set_from_pw_multi_aff(__isl_take isl_pw_multi_aff *pma)
3195 {
3196         if (!pma)
3197                 return NULL;
3198
3199         if (!isl_space_is_set(pma->dim))
3200                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
3201                         "isl_pw_multi_aff cannot be converted into an isl_set",
3202                         return isl_pw_multi_aff_free(pma));
3203
3204         return isl_map_from_pw_multi_aff(pma);
3205 }
3206
3207 /* Given a basic map with a single output dimension that is defined
3208  * in terms of the parameters and input dimensions using an equality,
3209  * extract an isl_aff that expresses the output dimension in terms
3210  * of the parameters and input dimensions.
3211  *
3212  * Since some applications expect the result of isl_pw_multi_aff_from_map
3213  * to only contain integer affine expressions, we compute the floor
3214  * of the expression before returning.
3215  *
3216  * This function shares some similarities with
3217  * isl_basic_map_has_defining_equality and isl_constraint_get_bound.
3218  */
3219 static __isl_give isl_aff *extract_isl_aff_from_basic_map(
3220         __isl_take isl_basic_map *bmap)
3221 {
3222         int i;
3223         unsigned offset;
3224         unsigned total;
3225         isl_local_space *ls;
3226         isl_aff *aff;
3227
3228         if (!bmap)
3229                 return NULL;
3230         if (isl_basic_map_dim(bmap, isl_dim_out) != 1)
3231                 isl_die(isl_basic_map_get_ctx(bmap), isl_error_invalid,
3232                         "basic map should have a single output dimension",
3233                         goto error);
3234         offset = isl_basic_map_offset(bmap, isl_dim_out);
3235         total = isl_basic_map_total_dim(bmap);
3236         for (i = 0; i < bmap->n_eq; ++i) {
3237                 if (isl_int_is_zero(bmap->eq[i][offset]))
3238                         continue;
3239                 if (isl_seq_first_non_zero(bmap->eq[i] + offset + 1,
3240                                            1 + total - (offset + 1)) != -1)
3241                         continue;
3242                 break;
3243         }
3244         if (i >= bmap->n_eq)
3245                 isl_die(isl_basic_map_get_ctx(bmap), isl_error_invalid,
3246                         "unable to find suitable equality", goto error);
3247         ls = isl_basic_map_get_local_space(bmap);
3248         aff = isl_aff_alloc(isl_local_space_domain(ls));
3249         if (!aff)
3250                 goto error;
3251         if (isl_int_is_neg(bmap->eq[i][offset]))
3252                 isl_seq_cpy(aff->v->el + 1, bmap->eq[i], offset);
3253         else
3254                 isl_seq_neg(aff->v->el + 1, bmap->eq[i], offset);
3255         isl_seq_clr(aff->v->el + 1 + offset, aff->v->size - (1 + offset));
3256         isl_int_abs(aff->v->el[0], bmap->eq[i][offset]);
3257         isl_basic_map_free(bmap);
3258
3259         aff = isl_aff_remove_unused_divs(aff);
3260         aff = isl_aff_floor(aff);
3261         return aff;
3262 error:
3263         isl_basic_map_free(bmap);
3264         return NULL;
3265 }
3266
3267 /* Given a basic map where each output dimension is defined
3268  * in terms of the parameters and input dimensions using an equality,
3269  * extract an isl_multi_aff that expresses the output dimensions in terms
3270  * of the parameters and input dimensions.
3271  */
3272 static __isl_give isl_multi_aff *extract_isl_multi_aff_from_basic_map(
3273         __isl_take isl_basic_map *bmap)
3274 {
3275         int i;
3276         unsigned n_out;
3277         isl_multi_aff *ma;
3278
3279         if (!bmap)
3280                 return NULL;
3281
3282         ma = isl_multi_aff_alloc(isl_basic_map_get_space(bmap));
3283         n_out = isl_basic_map_dim(bmap, isl_dim_out);
3284
3285         for (i = 0; i < n_out; ++i) {
3286                 isl_basic_map *bmap_i;
3287                 isl_aff *aff;
3288
3289                 bmap_i = isl_basic_map_copy(bmap);
3290                 bmap_i = isl_basic_map_project_out(bmap_i, isl_dim_out,
3291                                                         i + 1, n_out - (1 + i));
3292                 bmap_i = isl_basic_map_project_out(bmap_i, isl_dim_out, 0, i);
3293                 aff = extract_isl_aff_from_basic_map(bmap_i);
3294                 ma = isl_multi_aff_set_aff(ma, i, aff);
3295         }
3296
3297         isl_basic_map_free(bmap);
3298
3299         return ma;
3300 }
3301
3302 /* Create an isl_pw_multi_aff that is equivalent to
3303  * isl_map_intersect_domain(isl_map_from_basic_map(bmap), domain).
3304  * The given basic map is such that each output dimension is defined
3305  * in terms of the parameters and input dimensions using an equality.
3306  */
3307 static __isl_give isl_pw_multi_aff *plain_pw_multi_aff_from_map(
3308         __isl_take isl_set *domain, __isl_take isl_basic_map *bmap)
3309 {
3310         isl_multi_aff *ma;
3311
3312         ma = extract_isl_multi_aff_from_basic_map(bmap);
3313         return isl_pw_multi_aff_alloc(domain, ma);
3314 }
3315
3316 /* Try and create an isl_pw_multi_aff that is equivalent to the given isl_map.
3317  * This obviously only works if the input "map" is single-valued.
3318  * If so, we compute the lexicographic minimum of the image in the form
3319  * of an isl_pw_multi_aff.  Since the image is unique, it is equal
3320  * to its lexicographic minimum.
3321  * If the input is not single-valued, we produce an error.
3322  */
3323 static __isl_give isl_pw_multi_aff *pw_multi_aff_from_map_base(
3324         __isl_take isl_map *map)
3325 {
3326         int i;
3327         int sv;
3328         isl_pw_multi_aff *pma;
3329
3330         sv = isl_map_is_single_valued(map);
3331         if (sv < 0)
3332                 goto error;
3333         if (!sv)
3334                 isl_die(isl_map_get_ctx(map), isl_error_invalid,
3335                         "map is not single-valued", goto error);
3336         map = isl_map_make_disjoint(map);
3337         if (!map)
3338                 return NULL;
3339
3340         pma = isl_pw_multi_aff_empty(isl_map_get_space(map));
3341
3342         for (i = 0; i < map->n; ++i) {
3343                 isl_pw_multi_aff *pma_i;
3344                 isl_basic_map *bmap;
3345                 bmap = isl_basic_map_copy(map->p[i]);
3346                 pma_i = isl_basic_map_lexmin_pw_multi_aff(bmap);
3347                 pma = isl_pw_multi_aff_add_disjoint(pma, pma_i);
3348         }
3349
3350         isl_map_free(map);
3351         return pma;
3352 error:
3353         isl_map_free(map);
3354         return NULL;
3355 }
3356
3357 /* Try and create an isl_pw_multi_aff that is equivalent to the given isl_map,
3358  * taking into account that the output dimension at position "d"
3359  * can be represented as
3360  *
3361  *      x = floor((e(...) + c1) / m)
3362  *
3363  * given that constraint "i" is of the form
3364  *
3365  *      e(...) + c1 - m x >= 0
3366  *
3367  *
3368  * Let "map" be of the form
3369  *
3370  *      A -> B
3371  *
3372  * We construct a mapping
3373  *
3374  *      A -> [A -> x = floor(...)]
3375  *
3376  * apply that to the map, obtaining
3377  *
3378  *      [A -> x = floor(...)] -> B
3379  *
3380  * and equate dimension "d" to x.
3381  * We then compute a isl_pw_multi_aff representation of the resulting map
3382  * and plug in the mapping above.
3383  */
3384 static __isl_give isl_pw_multi_aff *pw_multi_aff_from_map_div(
3385         __isl_take isl_map *map, __isl_take isl_basic_map *hull, int d, int i)
3386 {
3387         isl_ctx *ctx;
3388         isl_space *space;
3389         isl_local_space *ls;
3390         isl_multi_aff *ma;
3391         isl_aff *aff;
3392         isl_vec *v;
3393         isl_map *insert;
3394         int offset;
3395         int n;
3396         int n_in;
3397         isl_pw_multi_aff *pma;
3398         int is_set;
3399
3400         is_set = isl_map_is_set(map);
3401
3402         offset = isl_basic_map_offset(hull, isl_dim_out);
3403         ctx = isl_map_get_ctx(map);
3404         space = isl_space_domain(isl_map_get_space(map));
3405         n_in = isl_space_dim(space, isl_dim_set);
3406         n = isl_space_dim(space, isl_dim_all);
3407
3408         v = isl_vec_alloc(ctx, 1 + 1 + n);
3409         if (v) {
3410                 isl_int_neg(v->el[0], hull->ineq[i][offset + d]);
3411                 isl_seq_cpy(v->el + 1, hull->ineq[i], 1 + n);
3412         }
3413         isl_basic_map_free(hull);
3414
3415         ls = isl_local_space_from_space(isl_space_copy(space));
3416         aff = isl_aff_alloc_vec(ls, v);
3417         aff = isl_aff_floor(aff);
3418         if (is_set) {
3419                 isl_space_free(space);
3420                 ma = isl_multi_aff_from_aff(aff);
3421         } else {
3422                 ma = isl_multi_aff_identity(isl_space_map_from_set(space));
3423                 ma = isl_multi_aff_range_product(ma,
3424                                                 isl_multi_aff_from_aff(aff));
3425         }
3426
3427         insert = isl_map_from_multi_aff(isl_multi_aff_copy(ma));
3428         map = isl_map_apply_domain(map, insert);
3429         map = isl_map_equate(map, isl_dim_in, n_in, isl_dim_out, d);
3430         pma = isl_pw_multi_aff_from_map(map);
3431         pma = isl_pw_multi_aff_pullback_multi_aff(pma, ma);
3432
3433         return pma;
3434 }
3435
3436 /* Is constraint "c" of the form
3437  *
3438  *      e(...) + c1 - m x >= 0
3439  *
3440  * or
3441  *
3442  *      -e(...) + c2 + m x >= 0
3443  *
3444  * where m > 1 and e only depends on parameters and input dimemnsions?
3445  *
3446  * "offset" is the offset of the output dimensions
3447  * "pos" is the position of output dimension x.
3448  */
3449 static int is_potential_div_constraint(isl_int *c, int offset, int d, int total)
3450 {
3451         if (isl_int_is_zero(c[offset + d]))
3452                 return 0;
3453         if (isl_int_is_one(c[offset + d]))
3454                 return 0;
3455         if (isl_int_is_negone(c[offset + d]))
3456                 return 0;
3457         if (isl_seq_first_non_zero(c + offset, d) != -1)
3458                 return 0;
3459         if (isl_seq_first_non_zero(c + offset + d + 1,
3460                                     total - (offset + d + 1)) != -1)
3461                 return 0;
3462         return 1;
3463 }
3464
3465 /* Try and create an isl_pw_multi_aff that is equivalent to the given isl_map.
3466  *
3467  * As a special case, we first check if there is any pair of constraints,
3468  * shared by all the basic maps in "map" that force a given dimension
3469  * to be equal to the floor of some affine combination of the input dimensions.
3470  *
3471  * In particular, if we can find two constraints
3472  *
3473  *      e(...) + c1 - m x >= 0          i.e.,           m x <= e(...) + c1
3474  *
3475  * and
3476  *
3477  *      -e(...) + c2 + m x >= 0         i.e.,           m x >= e(...) - c2
3478  *
3479  * where m > 1 and e only depends on parameters and input dimemnsions,
3480  * and such that
3481  *
3482  *      c1 + c2 < m                     i.e.,           -c2 >= c1 - (m - 1)
3483  *
3484  * then we know that we can take
3485  *
3486  *      x = floor((e(...) + c1) / m)
3487  *
3488  * without having to perform any computation.
3489  *
3490  * Note that we know that
3491  *
3492  *      c1 + c2 >= 1
3493  *
3494  * If c1 + c2 were 0, then we would have detected an equality during
3495  * simplification.  If c1 + c2 were negative, then we would have detected
3496  * a contradiction.
3497  */
3498 static __isl_give isl_pw_multi_aff *pw_multi_aff_from_map_check_div(
3499         __isl_take isl_map *map)
3500 {
3501         int d, dim;
3502         int i, j, n;
3503         int offset, total;
3504         isl_int sum;
3505         isl_basic_map *hull;
3506
3507         hull = isl_map_unshifted_simple_hull(isl_map_copy(map));
3508         if (!hull)
3509                 goto error;
3510
3511         isl_int_init(sum);
3512         dim = isl_map_dim(map, isl_dim_out);
3513         offset = isl_basic_map_offset(hull, isl_dim_out);
3514         total = 1 + isl_basic_map_total_dim(hull);
3515         n = hull->n_ineq;
3516         for (d = 0; d < dim; ++d) {
3517                 for (i = 0; i < n; ++i) {
3518                         if (!is_potential_div_constraint(hull->ineq[i],
3519                                                         offset, d, total))
3520                                 continue;
3521                         for (j = i + 1; j < n; ++j) {
3522                                 if (!isl_seq_is_neg(hull->ineq[i] + 1,
3523                                                 hull->ineq[j] + 1, total - 1))
3524                                         continue;
3525                                 isl_int_add(sum, hull->ineq[i][0],
3526                                                 hull->ineq[j][0]);
3527                                 if (isl_int_abs_lt(sum,
3528                                                     hull->ineq[i][offset + d]))
3529                                         break;
3530
3531                         }
3532                         if (j >= n)
3533                                 continue;
3534                         isl_int_clear(sum);
3535                         if (isl_int_is_pos(hull->ineq[j][offset + d]))
3536                                 j = i;
3537                         return pw_multi_aff_from_map_div(map, hull, d, j);
3538                 }
3539         }
3540         isl_int_clear(sum);
3541         isl_basic_map_free(hull);
3542         return pw_multi_aff_from_map_base(map);
3543 error:
3544         isl_map_free(map);
3545         isl_basic_map_free(hull);
3546         return NULL;
3547 }
3548
3549 /* Given an affine expression
3550  *
3551  *      [A -> B] -> f(A,B)
3552  *
3553  * construct an isl_multi_aff
3554  *
3555  *      [A -> B] -> B'
3556  *
3557  * such that dimension "d" in B' is set to "aff" and the remaining
3558  * dimensions are set equal to the corresponding dimensions in B.
3559  * "n_in" is the dimension of the space A.
3560  * "n_out" is the dimension of the space B.
3561  *
3562  * If "is_set" is set, then the affine expression is of the form
3563  *
3564  *      [B] -> f(B)
3565  *
3566  * and we construct an isl_multi_aff
3567  *
3568  *      B -> B'
3569  */
3570 static __isl_give isl_multi_aff *range_map(__isl_take isl_aff *aff, int d,
3571         unsigned n_in, unsigned n_out, int is_set)
3572 {
3573         int i;
3574         isl_multi_aff *ma;
3575         isl_space *space, *space2;
3576         isl_local_space *ls;
3577
3578         space = isl_aff_get_domain_space(aff);
3579         ls = isl_local_space_from_space(isl_space_copy(space));
3580         space2 = isl_space_copy(space);
3581         if (!is_set)
3582                 space2 = isl_space_range(isl_space_unwrap(space2));
3583         space = isl_space_map_from_domain_and_range(space, space2);
3584         ma = isl_multi_aff_alloc(space);
3585         ma = isl_multi_aff_set_aff(ma, d, aff);
3586
3587         for (i = 0; i < n_out; ++i) {
3588                 if (i == d)
3589                         continue;
3590                 aff = isl_aff_var_on_domain(isl_local_space_copy(ls),
3591                                                 isl_dim_set, n_in + i);
3592                 ma = isl_multi_aff_set_aff(ma, i, aff);
3593         }
3594
3595         isl_local_space_free(ls);
3596
3597         return ma;
3598 }
3599
3600 /* Try and create an isl_pw_multi_aff that is equivalent to the given isl_map,
3601  * taking into account that the dimension at position "d" can be written as
3602  *
3603  *      x = m a + f(..)                                         (1)
3604  *
3605  * where m is equal to "gcd".
3606  * "i" is the index of the equality in "hull" that defines f(..).
3607  * In particular, the equality is of the form
3608  *
3609  *      f(..) - x + m g(existentials) = 0
3610  *
3611  * or
3612  *
3613  *      -f(..) + x + m g(existentials) = 0
3614  *
3615  * We basically plug (1) into "map", resulting in a map with "a"
3616  * in the range instead of "x".  The corresponding isl_pw_multi_aff
3617  * defining "a" is then plugged back into (1) to obtain a definition fro "x".
3618  *
3619  * Specifically, given the input map
3620  *
3621  *      A -> B
3622  *
3623  * We first wrap it into a set
3624  *
3625  *      [A -> B]
3626  *
3627  * and define (1) on top of the corresponding space, resulting in "aff".
3628  * We use this to create an isl_multi_aff that maps the output position "d"
3629  * from "a" to "x", leaving all other (intput and output) dimensions unchanged.
3630  * We plug this into the wrapped map, unwrap the result and compute the
3631  * corresponding isl_pw_multi_aff.
3632  * The result is an expression
3633  *
3634  *      A -> T(A)
3635  *
3636  * We adjust that to
3637  *
3638  *      A -> [A -> T(A)]
3639  *
3640  * so that we can plug that into "aff", after extending the latter to
3641  * a mapping
3642  *
3643  *      [A -> B] -> B'
3644  *
3645  *
3646  * If "map" is actually a set, then there is no "A" space, meaning
3647  * that we do not need to perform any wrapping, and that the result
3648  * of the recursive call is of the form
3649  *
3650  *      [T]
3651  *
3652  * which is plugged into a mapping of the form
3653  *
3654  *      B -> B'
3655  */
3656 static __isl_give isl_pw_multi_aff *pw_multi_aff_from_map_stride(
3657         __isl_take isl_map *map, __isl_take isl_basic_map *hull, int d, int i,
3658         isl_int gcd)
3659 {
3660         isl_set *set;
3661         isl_space *space;
3662         isl_local_space *ls;
3663         isl_aff *aff;
3664         isl_multi_aff *ma;
3665         isl_pw_multi_aff *pma, *id;
3666         unsigned n_in;
3667         unsigned o_out;
3668         unsigned n_out;
3669         int is_set;
3670
3671         is_set = isl_map_is_set(map);
3672
3673         n_in = isl_basic_map_dim(hull, isl_dim_in);
3674         n_out = isl_basic_map_dim(hull, isl_dim_out);
3675         o_out = isl_basic_map_offset(hull, isl_dim_out);
3676
3677         if (is_set)
3678                 set = map;
3679         else
3680                 set = isl_map_wrap(map);
3681         space = isl_space_map_from_set(isl_set_get_space(set));
3682         ma = isl_multi_aff_identity(space);
3683         ls = isl_local_space_from_space(isl_set_get_space(set));
3684         aff = isl_aff_alloc(ls);
3685         if (aff) {
3686                 isl_int_set_si(aff->v->el[0], 1);
3687                 if (isl_int_is_one(hull->eq[i][o_out + d]))
3688                         isl_seq_neg(aff->v->el + 1, hull->eq[i],
3689                                     aff->v->size - 1);
3690                 else
3691                         isl_seq_cpy(aff->v->el + 1, hull->eq[i],
3692                                     aff->v->size - 1);
3693                 isl_int_set(aff->v->el[1 + o_out + d], gcd);
3694         }
3695         ma = isl_multi_aff_set_aff(ma, n_in + d, isl_aff_copy(aff));
3696         set = isl_set_preimage_multi_aff(set, ma);
3697
3698         ma = range_map(aff, d, n_in, n_out, is_set);
3699
3700         if (is_set)
3701                 map = set;
3702         else
3703                 map = isl_set_unwrap(set);
3704         pma = isl_pw_multi_aff_from_map(set);
3705
3706         if (!is_set) {
3707                 space = isl_pw_multi_aff_get_domain_space(pma);
3708                 space = isl_space_map_from_set(space);
3709                 id = isl_pw_multi_aff_identity(space);
3710                 pma = isl_pw_multi_aff_range_product(id, pma);
3711         }
3712         id = isl_pw_multi_aff_from_multi_aff(ma);
3713         pma = isl_pw_multi_aff_pullback_pw_multi_aff(id, pma);
3714
3715         isl_basic_map_free(hull);
3716         return pma;
3717 }
3718
3719 /* Try and create an isl_pw_multi_aff that is equivalent to the given isl_map.
3720  *
3721  * As a special case, we first check if all output dimensions are uniquely
3722  * defined in terms of the parameters and input dimensions over the entire
3723  * domain.  If so, we extract the desired isl_pw_multi_aff directly
3724  * from the affine hull of "map" and its domain.
3725  *
3726  * Otherwise, we check if any of the output dimensions is "strided".
3727  * That is, we check if can be written as
3728  *
3729  *      x = m a + f(..)
3730  *
3731  * with m greater than 1, a some combination of existentiall quantified
3732  * variables and f and expression in the parameters and input dimensions.
3733  * If so, we remove the stride in pw_multi_aff_from_map_stride.
3734  *
3735  * Otherwise, we continue with pw_multi_aff_from_map_check_div for a further
3736  * special case.
3737  */
3738 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_from_map(__isl_take isl_map *map)
3739 {
3740         int i, j;
3741         int sv;
3742         isl_basic_map *hull;
3743         unsigned n_out;
3744         unsigned o_out;
3745         unsigned n_div;
3746         unsigned o_div;
3747         isl_int gcd;
3748
3749         if (!map)
3750                 return NULL;
3751
3752         hull = isl_map_affine_hull(isl_map_copy(map));
3753         sv = isl_basic_map_plain_is_single_valued(hull);
3754         if (sv >= 0 && sv)
3755                 return plain_pw_multi_aff_from_map(isl_map_domain(map), hull);
3756         if (sv < 0)
3757                 hull = isl_basic_map_free(hull);
3758         if (!hull)
3759                 goto error;
3760
3761         n_div = isl_basic_map_dim(hull, isl_dim_div);
3762         o_div = isl_basic_map_offset(hull, isl_dim_div);
3763
3764         if (n_div == 0) {
3765                 isl_basic_map_free(hull);
3766                 return pw_multi_aff_from_map_check_div(map);
3767         }
3768
3769         isl_int_init(gcd);
3770
3771         n_out = isl_basic_map_dim(hull, isl_dim_out);
3772         o_out = isl_basic_map_offset(hull, isl_dim_out);
3773
3774         for (i = 0; i < n_out; ++i) {
3775                 for (j = 0; j < hull->n_eq; ++j) {
3776                         isl_int *eq = hull->eq[j];
3777                         isl_pw_multi_aff *res;
3778
3779                         if (!isl_int_is_one(eq[o_out + i]) &&
3780                             !isl_int_is_negone(eq[o_out + i]))
3781                                 continue;
3782                         if (isl_seq_first_non_zero(eq + o_out, i) != -1)
3783                                 continue;
3784                         if (isl_seq_first_non_zero(eq + o_out + i + 1,
3785                                                     n_out - (i + 1)) != -1)
3786                                 continue;
3787                         isl_seq_gcd(eq + o_div, n_div, &gcd);
3788                         if (isl_int_is_zero(gcd))
3789                                 continue;
3790                         if (isl_int_is_one(gcd))
3791                                 continue;
3792
3793                         res = pw_multi_aff_from_map_stride(map, hull,
3794                                                                 i, j, gcd);
3795                         isl_int_clear(gcd);
3796                         return res;
3797                 }
3798         }
3799
3800         isl_int_clear(gcd);
3801         isl_basic_map_free(hull);
3802         return pw_multi_aff_from_map_check_div(map);
3803 error:
3804         isl_map_free(map);
3805         return NULL;
3806 }
3807
3808 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_from_set(__isl_take isl_set *set)
3809 {
3810         return isl_pw_multi_aff_from_map(set);
3811 }
3812
3813 /* Convert "map" into an isl_pw_multi_aff (if possible) and
3814  * add it to *user.
3815  */
3816 static int pw_multi_aff_from_map(__isl_take isl_map *map, void *user)
3817 {
3818         isl_union_pw_multi_aff **upma = user;
3819         isl_pw_multi_aff *pma;
3820
3821         pma = isl_pw_multi_aff_from_map(map);
3822         *upma = isl_union_pw_multi_aff_add_pw_multi_aff(*upma, pma);
3823
3824         return *upma ? 0 : -1;
3825 }
3826
3827 /* Try and create an isl_union_pw_multi_aff that is equivalent
3828  * to the given isl_union_map.
3829  * The isl_union_map is required to be single-valued in each space.
3830  * Otherwise, an error is produced.
3831  */
3832 __isl_give isl_union_pw_multi_aff *isl_union_pw_multi_aff_from_union_map(
3833         __isl_take isl_union_map *umap)
3834 {
3835         isl_space *space;
3836         isl_union_pw_multi_aff *upma;
3837
3838         space = isl_union_map_get_space(umap);
3839         upma = isl_union_pw_multi_aff_empty(space);
3840         if (isl_union_map_foreach_map(umap, &pw_multi_aff_from_map, &upma) < 0)
3841                 upma = isl_union_pw_multi_aff_free(upma);
3842         isl_union_map_free(umap);
3843
3844         return upma;
3845 }
3846
3847 /* Try and create an isl_union_pw_multi_aff that is equivalent
3848  * to the given isl_union_set.
3849  * The isl_union_set is required to be a singleton in each space.
3850  * Otherwise, an error is produced.
3851  */
3852 __isl_give isl_union_pw_multi_aff *isl_union_pw_multi_aff_from_union_set(
3853         __isl_take isl_union_set *uset)
3854 {
3855         return isl_union_pw_multi_aff_from_union_map(uset);
3856 }
3857
3858 /* Return the piecewise affine expression "set ? 1 : 0".
3859  */
3860 __isl_give isl_pw_aff *isl_set_indicator_function(__isl_take isl_set *set)
3861 {
3862         isl_pw_aff *pa;
3863         isl_space *space = isl_set_get_space(set);
3864         isl_local_space *ls = isl_local_space_from_space(space);
3865         isl_aff *zero = isl_aff_zero_on_domain(isl_local_space_copy(ls));
3866         isl_aff *one = isl_aff_zero_on_domain(ls);
3867
3868         one = isl_aff_add_constant_si(one, 1);
3869         pa = isl_pw_aff_alloc(isl_set_copy(set), one);
3870         set = isl_set_complement(set);
3871         pa = isl_pw_aff_add_disjoint(pa, isl_pw_aff_alloc(set, zero));
3872
3873         return pa;
3874 }
3875
3876 /* Plug in "subs" for dimension "type", "pos" of "aff".
3877  *
3878  * Let i be the dimension to replace and let "subs" be of the form
3879  *
3880  *      f/d
3881  *
3882  * and "aff" of the form
3883  *
3884  *      (a i + g)/m
3885  *
3886  * The result is
3887  *
3888  *      (a f + d g')/(m d)
3889  *
3890  * where g' is the result of plugging in "subs" in each of the integer
3891  * divisions in g.
3892  */
3893 __isl_give isl_aff *isl_aff_substitute(__isl_take isl_aff *aff,
3894         enum isl_dim_type type, unsigned pos, __isl_keep isl_aff *subs)
3895 {
3896         isl_ctx *ctx;
3897         isl_int v;
3898
3899         aff = isl_aff_cow(aff);
3900         if (!aff || !subs)
3901                 return isl_aff_free(aff);
3902
3903         ctx = isl_aff_get_ctx(aff);
3904         if (!isl_space_is_equal(aff->ls->dim, subs->ls->dim))
3905                 isl_die(ctx, isl_error_invalid,
3906                         "spaces don't match", return isl_aff_free(aff));
3907         if (isl_local_space_dim(subs->ls, isl_dim_div) != 0)
3908                 isl_die(ctx, isl_error_unsupported,
3909                         "cannot handle divs yet", return isl_aff_free(aff));
3910
3911         aff->ls = isl_local_space_substitute(aff->ls, type, pos, subs);
3912         if (!aff->ls)
3913                 return isl_aff_free(aff);
3914
3915         aff->v = isl_vec_cow(aff->v);
3916         if (!aff->v)
3917                 return isl_aff_free(aff);
3918
3919         pos += isl_local_space_offset(aff->ls, type);
3920
3921         isl_int_init(v);
3922         isl_seq_substitute(aff->v->el, pos, subs->v->el,
3923                             aff->v->size, subs->v->size, v);
3924         isl_int_clear(v);
3925
3926         return aff;
3927 }
3928
3929 /* Plug in "subs" for dimension "type", "pos" in each of the affine
3930  * expressions in "maff".
3931  */
3932 __isl_give isl_multi_aff *isl_multi_aff_substitute(
3933         __isl_take isl_multi_aff *maff, enum isl_dim_type type, unsigned pos,
3934         __isl_keep isl_aff *subs)
3935 {
3936         int i;
3937
3938         maff = isl_multi_aff_cow(maff);
3939         if (!maff || !subs)
3940                 return isl_multi_aff_free(maff);
3941
3942         if (type == isl_dim_in)
3943                 type = isl_dim_set;
3944
3945         for (i = 0; i < maff->n; ++i) {
3946                 maff->p[i] = isl_aff_substitute(maff->p[i], type, pos, subs);
3947                 if (!maff->p[i])
3948                         return isl_multi_aff_free(maff);
3949         }
3950
3951         return maff;
3952 }
3953
3954 /* Plug in "subs" for dimension "type", "pos" of "pma".
3955  *
3956  * pma is of the form
3957  *
3958  *      A_i(v) -> M_i(v)
3959  *
3960  * while subs is of the form
3961  *
3962  *      v' = B_j(v) -> S_j
3963  *
3964  * Each pair i,j such that C_ij = A_i \cap B_i is non-empty
3965  * has a contribution in the result, in particular
3966  *
3967  *      C_ij(S_j) -> M_i(S_j)
3968  *
3969  * Note that plugging in S_j in C_ij may also result in an empty set
3970  * and this contribution should simply be discarded.
3971  */
3972 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_substitute(
3973         __isl_take isl_pw_multi_aff *pma, enum isl_dim_type type, unsigned pos,
3974         __isl_keep isl_pw_aff *subs)
3975 {
3976         int i, j, n;
3977         isl_pw_multi_aff *res;
3978
3979         if (!pma || !subs)
3980                 return isl_pw_multi_aff_free(pma);
3981
3982         n = pma->n * subs->n;
3983         res = isl_pw_multi_aff_alloc_size(isl_space_copy(pma->dim), n);
3984
3985         for (i = 0; i < pma->n; ++i) {
3986                 for (j = 0; j < subs->n; ++j) {
3987                         isl_set *common;
3988                         isl_multi_aff *res_ij;
3989                         int empty;
3990
3991                         common = isl_set_intersect(
3992                                         isl_set_copy(pma->p[i].set),
3993                                         isl_set_copy(subs->p[j].set));
3994                         common = isl_set_substitute(common,
3995                                         type, pos, subs->p[j].aff);
3996                         empty = isl_set_plain_is_empty(common);
3997                         if (empty < 0 || empty) {
3998                                 isl_set_free(common);
3999                                 if (empty < 0)
4000                                         goto error;
4001                                 continue;
4002                         }
4003
4004                         res_ij = isl_multi_aff_substitute(
4005                                         isl_multi_aff_copy(pma->p[i].maff),
4006                                         type, pos, subs->p[j].aff);
4007
4008                         res = isl_pw_multi_aff_add_piece(res, common, res_ij);
4009                 }
4010         }
4011
4012         isl_pw_multi_aff_free(pma);
4013         return res;
4014 error:
4015         isl_pw_multi_aff_free(pma);
4016         isl_pw_multi_aff_free(res);
4017         return NULL;
4018 }
4019
4020 /* Compute the preimage of a range of dimensions in the affine expression "src"
4021  * under "ma" and put the result in "dst".  The number of dimensions in "src"
4022  * that precede the range is given by "n_before".  The number of dimensions
4023  * in the range is given by the number of output dimensions of "ma".
4024  * The number of dimensions that follow the range is given by "n_after".
4025  * If "has_denom" is set (to one),
4026  * then "src" and "dst" have an extra initial denominator.
4027  * "n_div_ma" is the number of existentials in "ma"
4028  * "n_div_bset" is the number of existentials in "src"
4029  * The resulting "dst" (which is assumed to have been allocated by
4030  * the caller) contains coefficients for both sets of existentials,
4031  * first those in "ma" and then those in "src".
4032  * f, c1, c2 and g are temporary objects that have been initialized
4033  * by the caller.
4034  *
4035  * Let src represent the expression
4036  *
4037  *      (a(p) + f_u u + b v + f_w w + c(divs))/d
4038  *
4039  * and let ma represent the expressions
4040  *
4041  *      v_i = (r_i(p) + s_i(y) + t_i(divs'))/m_i
4042  *
4043  * We start out with the following expression for dst:
4044  *
4045  *      (a(p) + f_u u + 0 y + f_w w + 0 divs' + c(divs) + f \sum_i b_i v_i)/d
4046  *
4047  * with the multiplication factor f initially equal to 1
4048  * and f \sum_i b_i v_i kept separately.
4049  * For each x_i that we substitute, we multiply the numerator
4050  * (and denominator) of dst by c_1 = m_i and add the numerator
4051  * of the x_i expression multiplied by c_2 = f b_i,
4052  * after removing the common factors of c_1 and c_2.
4053  * The multiplication factor f also needs to be multiplied by c_1
4054  * for the next x_j, j > i.
4055  */
4056 void isl_seq_preimage(isl_int *dst, isl_int *src,
4057         __isl_keep isl_multi_aff *ma, int n_before, int n_after,
4058         int n_div_ma, int n_div_bmap,
4059         isl_int f, isl_int c1, isl_int c2, isl_int g, int has_denom)
4060 {
4061         int i;
4062         int n_param, n_in, n_out;
4063         int o_dst, o_src;
4064
4065         n_param = isl_multi_aff_dim(ma, isl_dim_param);
4066         n_in = isl_multi_aff_dim(ma, isl_dim_in);
4067         n_out = isl_multi_aff_dim(ma, isl_dim_out);
4068
4069         isl_seq_cpy(dst, src, has_denom + 1 + n_param + n_before);
4070         o_dst = o_src = has_denom + 1 + n_param + n_before;
4071         isl_seq_clr(dst + o_dst, n_in);
4072         o_dst += n_in;
4073         o_src += n_out;
4074         isl_seq_cpy(dst + o_dst, src + o_src, n_after);
4075         o_dst += n_after;
4076         o_src += n_after;
4077         isl_seq_clr(dst + o_dst, n_div_ma);
4078         o_dst += n_div_ma;
4079         isl_seq_cpy(dst + o_dst, src + o_src, n_div_bmap);
4080
4081         isl_int_set_si(f, 1);
4082
4083         for (i = 0; i < n_out; ++i) {
4084                 int offset = has_denom + 1 + n_param + n_before + i;
4085
4086                 if (isl_int_is_zero(src[offset]))
4087                         continue;
4088                 isl_int_set(c1, ma->p[i]->v->el[0]);
4089                 isl_int_mul(c2, f, src[offset]);
4090                 isl_int_gcd(g, c1, c2);
4091                 isl_int_divexact(c1, c1, g);
4092                 isl_int_divexact(c2, c2, g);
4093
4094                 isl_int_mul(f, f, c1);
4095                 o_dst = has_denom;
4096                 o_src = 1;
4097                 isl_seq_combine(dst + o_dst, c1, dst + o_dst,
4098                                 c2, ma->p[i]->v->el + o_src, 1 + n_param);
4099                 o_dst += 1 + n_param;
4100                 o_src += 1 + n_param;
4101                 isl_seq_scale(dst + o_dst, dst + o_dst, c1, n_before);
4102                 o_dst += n_before;
4103                 isl_seq_combine(dst + o_dst, c1, dst + o_dst,
4104                                 c2, ma->p[i]->v->el + o_src, n_in);
4105                 o_dst += n_in;
4106                 o_src += n_in;
4107                 isl_seq_scale(dst + o_dst, dst + o_dst, c1, n_after);
4108                 o_dst += n_after;
4109                 isl_seq_combine(dst + o_dst, c1, dst + o_dst,
4110                                 c2, ma->p[i]->v->el + o_src, n_div_ma);
4111                 o_dst += n_div_ma;
4112                 o_src += n_div_ma;
4113                 isl_seq_scale(dst + o_dst, dst + o_dst, c1, n_div_bmap);
4114                 if (has_denom)
4115                         isl_int_mul(dst[0], dst[0], c1);
4116         }
4117 }
4118
4119 /* Compute the pullback of "aff" by the function represented by "ma".
4120  * In other words, plug in "ma" in "aff".  The result is an affine expression
4121  * defined over the domain space of "ma".
4122  *
4123  * If "aff" is represented by
4124  *
4125  *      (a(p) + b x + c(divs))/d
4126  *
4127  * and ma is represented by
4128  *
4129  *      x = D(p) + F(y) + G(divs')
4130  *
4131  * then the result is
4132  *
4133  *      (a(p) + b D(p) + b F(y) + b G(divs') + c(divs))/d
4134  *
4135  * The divs in the local space of the input are similarly adjusted
4136  * through a call to isl_local_space_preimage_multi_aff.
4137  */
4138 __isl_give isl_aff *isl_aff_pullback_multi_aff(__isl_take isl_aff *aff,
4139         __isl_take isl_multi_aff *ma)
4140 {
4141         isl_aff *res = NULL;
4142         isl_local_space *ls;
4143         int n_div_aff, n_div_ma;
4144         isl_int f, c1, c2, g;
4145
4146         ma = isl_multi_aff_align_divs(ma);
4147         if (!aff || !ma)
4148                 goto error;
4149
4150         n_div_aff = isl_aff_dim(aff, isl_dim_div);
4151         n_div_ma = ma->n ? isl_aff_dim(ma->p[0], isl_dim_div) : 0;
4152
4153         ls = isl_aff_get_domain_local_space(aff);
4154         ls = isl_local_space_preimage_multi_aff(ls, isl_multi_aff_copy(ma));
4155         res = isl_aff_alloc(ls);
4156         if (!res)
4157                 goto error;
4158
4159         isl_int_init(f);
4160         isl_int_init(c1);
4161         isl_int_init(c2);
4162         isl_int_init(g);
4163
4164         isl_seq_preimage(res->v->el, aff->v->el, ma, 0, 0, n_div_ma, n_div_aff,
4165                         f, c1, c2, g, 1);
4166
4167         isl_int_clear(f);
4168         isl_int_clear(c1);
4169         isl_int_clear(c2);
4170         isl_int_clear(g);
4171
4172         isl_aff_free(aff);
4173         isl_multi_aff_free(ma);
4174         res = isl_aff_normalize(res);
4175         return res;
4176 error:
4177         isl_aff_free(aff);
4178         isl_multi_aff_free(ma);
4179         isl_aff_free(res);
4180         return NULL;
4181 }
4182
4183 /* Compute the pullback of "ma1" by the function represented by "ma2".
4184  * In other words, plug in "ma2" in "ma1".
4185  */
4186 __isl_give isl_multi_aff *isl_multi_aff_pullback_multi_aff(
4187         __isl_take isl_multi_aff *ma1, __isl_take isl_multi_aff *ma2)
4188 {
4189         int i;
4190         isl_space *space = NULL;
4191
4192         ma2 = isl_multi_aff_align_divs(ma2);
4193         ma1 = isl_multi_aff_cow(ma1);
4194         if (!ma1 || !ma2)
4195                 goto error;
4196
4197         space = isl_space_join(isl_multi_aff_get_space(ma2),
4198                                 isl_multi_aff_get_space(ma1));
4199
4200         for (i = 0; i < ma1->n; ++i) {
4201                 ma1->p[i] = isl_aff_pullback_multi_aff(ma1->p[i],
4202                                                     isl_multi_aff_copy(ma2));
4203                 if (!ma1->p[i])
4204                         goto error;
4205         }
4206
4207         ma1 = isl_multi_aff_reset_space(ma1, space);
4208         isl_multi_aff_free(ma2);
4209         return ma1;
4210 error:
4211         isl_space_free(space);
4212         isl_multi_aff_free(ma2);
4213         isl_multi_aff_free(ma1);
4214         return NULL;
4215 }
4216
4217 /* Extend the local space of "dst" to include the divs
4218  * in the local space of "src".
4219  */
4220 __isl_give isl_aff *isl_aff_align_divs(__isl_take isl_aff *dst,
4221         __isl_keep isl_aff *src)
4222 {
4223         isl_ctx *ctx;
4224         int *exp1 = NULL;
4225         int *exp2 = NULL;
4226         isl_mat *div;
4227
4228         if (!src || !dst)
4229                 return isl_aff_free(dst);
4230
4231         ctx = isl_aff_get_ctx(src);
4232         if (!isl_space_is_equal(src->ls->dim, dst->ls->dim))
4233                 isl_die(ctx, isl_error_invalid,
4234                         "spaces don't match", goto error);
4235
4236         if (src->ls->div->n_row == 0)
4237                 return dst;
4238
4239         exp1 = isl_alloc_array(ctx, int, src->ls->div->n_row);
4240         exp2 = isl_alloc_array(ctx, int, dst->ls->div->n_row);
4241         if (!exp1 || !exp2)
4242                 goto error;
4243
4244         div = isl_merge_divs(src->ls->div, dst->ls->div, exp1, exp2);
4245         dst = isl_aff_expand_divs(dst, div, exp2);
4246         free(exp1);
4247         free(exp2);
4248
4249         return dst;
4250 error:
4251         free(exp1);
4252         free(exp2);
4253         return isl_aff_free(dst);
4254 }
4255
4256 /* Adjust the local spaces of the affine expressions in "maff"
4257  * such that they all have the save divs.
4258  */
4259 __isl_give isl_multi_aff *isl_multi_aff_align_divs(
4260         __isl_take isl_multi_aff *maff)
4261 {
4262         int i;
4263
4264         if (!maff)
4265                 return NULL;
4266         if (maff->n == 0)
4267                 return maff;
4268         maff = isl_multi_aff_cow(maff);
4269         if (!maff)
4270                 return NULL;
4271
4272         for (i = 1; i < maff->n; ++i)
4273                 maff->p[0] = isl_aff_align_divs(maff->p[0], maff->p[i]);
4274         for (i = 1; i < maff->n; ++i) {
4275                 maff->p[i] = isl_aff_align_divs(maff->p[i], maff->p[0]);
4276                 if (!maff->p[i])
4277                         return isl_multi_aff_free(maff);
4278         }
4279
4280         return maff;
4281 }
4282
4283 __isl_give isl_aff *isl_aff_lift(__isl_take isl_aff *aff)
4284 {
4285         aff = isl_aff_cow(aff);
4286         if (!aff)
4287                 return NULL;
4288
4289         aff->ls = isl_local_space_lift(aff->ls);
4290         if (!aff->ls)
4291                 return isl_aff_free(aff);
4292
4293         return aff;
4294 }
4295
4296 /* Lift "maff" to a space with extra dimensions such that the result
4297  * has no more existentially quantified variables.
4298  * If "ls" is not NULL, then *ls is assigned the local space that lies
4299  * at the basis of the lifting applied to "maff".
4300  */
4301 __isl_give isl_multi_aff *isl_multi_aff_lift(__isl_take isl_multi_aff *maff,
4302         __isl_give isl_local_space **ls)
4303 {
4304         int i;
4305         isl_space *space;
4306         unsigned n_div;
4307
4308         if (ls)
4309                 *ls = NULL;
4310
4311         if (!maff)
4312                 return NULL;
4313
4314         if (maff->n == 0) {
4315                 if (ls) {
4316                         isl_space *space = isl_multi_aff_get_domain_space(maff);
4317                         *ls = isl_local_space_from_space(space);
4318                         if (!*ls)
4319                                 return isl_multi_aff_free(maff);
4320                 }
4321                 return maff;
4322         }
4323
4324         maff = isl_multi_aff_cow(maff);
4325         maff = isl_multi_aff_align_divs(maff);
4326         if (!maff)
4327                 return NULL;
4328
4329         n_div = isl_aff_dim(maff->p[0], isl_dim_div);
4330         space = isl_multi_aff_get_space(maff);
4331         space = isl_space_lift(isl_space_domain(space), n_div);
4332         space = isl_space_extend_domain_with_range(space,
4333                                                 isl_multi_aff_get_space(maff));
4334         if (!space)
4335                 return isl_multi_aff_free(maff);
4336         isl_space_free(maff->space);
4337         maff->space = space;
4338
4339         if (ls) {
4340                 *ls = isl_aff_get_domain_local_space(maff->p[0]);
4341                 if (!*ls)
4342                         return isl_multi_aff_free(maff);
4343         }
4344
4345         for (i = 0; i < maff->n; ++i) {
4346                 maff->p[i] = isl_aff_lift(maff->p[i]);
4347                 if (!maff->p[i])
4348                         goto error;
4349         }
4350
4351         return maff;
4352 error:
4353         if (ls)
4354                 isl_local_space_free(*ls);
4355         return isl_multi_aff_free(maff);
4356 }
4357
4358
4359 /* Extract an isl_pw_aff corresponding to output dimension "pos" of "pma".
4360  */
4361 __isl_give isl_pw_aff *isl_pw_multi_aff_get_pw_aff(
4362         __isl_keep isl_pw_multi_aff *pma, int pos)
4363 {
4364         int i;
4365         int n_out;
4366         isl_space *space;
4367         isl_pw_aff *pa;
4368
4369         if (!pma)
4370                 return NULL;
4371
4372         n_out = isl_pw_multi_aff_dim(pma, isl_dim_out);
4373         if (pos < 0 || pos >= n_out)
4374                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
4375                         "index out of bounds", return NULL);
4376
4377         space = isl_pw_multi_aff_get_space(pma);
4378         space = isl_space_drop_dims(space, isl_dim_out,
4379                                     pos + 1, n_out - pos - 1);
4380         space = isl_space_drop_dims(space, isl_dim_out, 0, pos);
4381
4382         pa = isl_pw_aff_alloc_size(space, pma->n);
4383         for (i = 0; i < pma->n; ++i) {
4384                 isl_aff *aff;
4385                 aff = isl_multi_aff_get_aff(pma->p[i].maff, pos);
4386                 pa = isl_pw_aff_add_piece(pa, isl_set_copy(pma->p[i].set), aff);
4387         }
4388
4389         return pa;
4390 }
4391
4392 /* Return an isl_pw_multi_aff with the given "set" as domain and
4393  * an unnamed zero-dimensional range.
4394  */
4395 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_from_domain(
4396         __isl_take isl_set *set)
4397 {
4398         isl_multi_aff *ma;
4399         isl_space *space;
4400
4401         space = isl_set_get_space(set);
4402         space = isl_space_from_domain(space);
4403         ma = isl_multi_aff_zero(space);
4404         return isl_pw_multi_aff_alloc(set, ma);
4405 }
4406
4407 /* Add an isl_pw_multi_aff with the given "set" as domain and
4408  * an unnamed zero-dimensional range to *user.
4409  */
4410 static int add_pw_multi_aff_from_domain(__isl_take isl_set *set, void *user)
4411 {
4412         isl_union_pw_multi_aff **upma = user;
4413         isl_pw_multi_aff *pma;
4414
4415         pma = isl_pw_multi_aff_from_domain(set);
4416         *upma = isl_union_pw_multi_aff_add_pw_multi_aff(*upma, pma);
4417
4418         return 0;
4419 }
4420
4421 /* Return an isl_union_pw_multi_aff with the given "uset" as domain and
4422  * an unnamed zero-dimensional range.
4423  */
4424 __isl_give isl_union_pw_multi_aff *isl_union_pw_multi_aff_from_domain(
4425         __isl_take isl_union_set *uset)
4426 {
4427         isl_space *space;
4428         isl_union_pw_multi_aff *upma;
4429
4430         if (!uset)
4431                 return NULL;
4432
4433         space = isl_union_set_get_space(uset);
4434         upma = isl_union_pw_multi_aff_empty(space);
4435
4436         if (isl_union_set_foreach_set(uset,
4437                                     &add_pw_multi_aff_from_domain, &upma) < 0)
4438                 goto error;
4439
4440         isl_union_set_free(uset);
4441         return upma;
4442 error:
4443         isl_union_set_free(uset);
4444         isl_union_pw_multi_aff_free(upma);
4445         return NULL;
4446 }
4447
4448 /* Convert "pma" to an isl_map and add it to *umap.
4449  */
4450 static int map_from_pw_multi_aff(__isl_take isl_pw_multi_aff *pma, void *user)
4451 {
4452         isl_union_map **umap = user;
4453         isl_map *map;
4454
4455         map = isl_map_from_pw_multi_aff(pma);
4456         *umap = isl_union_map_add_map(*umap, map);
4457
4458         return 0;
4459 }
4460
4461 /* Construct a union map mapping the domain of the union
4462  * piecewise multi-affine expression to its range, with each dimension
4463  * in the range equated to the corresponding affine expression on its cell.
4464  */
4465 __isl_give isl_union_map *isl_union_map_from_union_pw_multi_aff(
4466         __isl_take isl_union_pw_multi_aff *upma)
4467 {
4468         isl_space *space;
4469         isl_union_map *umap;
4470
4471         if (!upma)
4472                 return NULL;
4473
4474         space = isl_union_pw_multi_aff_get_space(upma);
4475         umap = isl_union_map_empty(space);
4476
4477         if (isl_union_pw_multi_aff_foreach_pw_multi_aff(upma,
4478                                         &map_from_pw_multi_aff, &umap) < 0)
4479                 goto error;
4480
4481         isl_union_pw_multi_aff_free(upma);
4482         return umap;
4483 error:
4484         isl_union_pw_multi_aff_free(upma);
4485         isl_union_map_free(umap);
4486         return NULL;
4487 }
4488
4489 /* Local data for bin_entry and the callback "fn".
4490  */
4491 struct isl_union_pw_multi_aff_bin_data {
4492         isl_union_pw_multi_aff *upma2;
4493         isl_union_pw_multi_aff *res;
4494         isl_pw_multi_aff *pma;
4495         int (*fn)(void **entry, void *user);
4496 };
4497
4498 /* Given an isl_pw_multi_aff from upma1, store it in data->pma
4499  * and call data->fn for each isl_pw_multi_aff in data->upma2.
4500  */
4501 static int bin_entry(void **entry, void *user)
4502 {
4503         struct isl_union_pw_multi_aff_bin_data *data = user;
4504         isl_pw_multi_aff *pma = *entry;
4505
4506         data->pma = pma;
4507         if (isl_hash_table_foreach(data->upma2->dim->ctx, &data->upma2->table,
4508                                    data->fn, data) < 0)
4509                 return -1;
4510
4511         return 0;
4512 }
4513
4514 /* Call "fn" on each pair of isl_pw_multi_affs in "upma1" and "upma2".
4515  * The isl_pw_multi_aff from upma1 is stored in data->pma (where data is
4516  * passed as user field) and the isl_pw_multi_aff from upma2 is available
4517  * as *entry.  The callback should adjust data->res if desired.
4518  */
4519 static __isl_give isl_union_pw_multi_aff *bin_op(
4520         __isl_take isl_union_pw_multi_aff *upma1,
4521         __isl_take isl_union_pw_multi_aff *upma2,
4522         int (*fn)(void **entry, void *user))
4523 {
4524         isl_space *space;
4525         struct isl_union_pw_multi_aff_bin_data data = { NULL, NULL, NULL, fn };
4526
4527         space = isl_union_pw_multi_aff_get_space(upma2);
4528         upma1 = isl_union_pw_multi_aff_align_params(upma1, space);
4529         space = isl_union_pw_multi_aff_get_space(upma1);
4530         upma2 = isl_union_pw_multi_aff_align_params(upma2, space);
4531
4532         if (!upma1 || !upma2)
4533                 goto error;
4534
4535         data.upma2 = upma2;
4536         data.res = isl_union_pw_multi_aff_alloc(isl_space_copy(upma1->dim),
4537                                        upma1->table.n);
4538         if (isl_hash_table_foreach(upma1->dim->ctx, &upma1->table,
4539                                    &bin_entry, &data) < 0)
4540                 goto error;
4541
4542         isl_union_pw_multi_aff_free(upma1);
4543         isl_union_pw_multi_aff_free(upma2);
4544         return data.res;
4545 error:
4546         isl_union_pw_multi_aff_free(upma1);
4547         isl_union_pw_multi_aff_free(upma2);
4548         isl_union_pw_multi_aff_free(data.res);
4549         return NULL;
4550 }
4551
4552 /* Given two aligned isl_pw_multi_affs A -> B and C -> D,
4553  * construct an isl_pw_multi_aff (A * C) -> [B -> D].
4554  */
4555 static __isl_give isl_pw_multi_aff *pw_multi_aff_range_product(
4556         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
4557 {
4558         isl_space *space;
4559
4560         space = isl_space_range_product(isl_pw_multi_aff_get_space(pma1),
4561                                         isl_pw_multi_aff_get_space(pma2));
4562         return isl_pw_multi_aff_on_shared_domain_in(pma1, pma2, space,
4563                                             &isl_multi_aff_range_product);
4564 }
4565
4566 /* Given two isl_pw_multi_affs A -> B and C -> D,
4567  * construct an isl_pw_multi_aff (A * C) -> [B -> D].
4568  */
4569 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_range_product(
4570         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
4571 {
4572         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
4573                                             &pw_multi_aff_range_product);
4574 }
4575
4576 /* Given two aligned isl_pw_multi_affs A -> B and C -> D,
4577  * construct an isl_pw_multi_aff (A * C) -> (B, D).
4578  */
4579 static __isl_give isl_pw_multi_aff *pw_multi_aff_flat_range_product(
4580         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
4581 {
4582         isl_space *space;
4583
4584         space = isl_space_range_product(isl_pw_multi_aff_get_space(pma1),
4585                                         isl_pw_multi_aff_get_space(pma2));
4586         space = isl_space_flatten_range(space);
4587         return isl_pw_multi_aff_on_shared_domain_in(pma1, pma2, space,
4588                                             &isl_multi_aff_flat_range_product);
4589 }
4590
4591 /* Given two isl_pw_multi_affs A -> B and C -> D,
4592  * construct an isl_pw_multi_aff (A * C) -> (B, D).
4593  */
4594 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_flat_range_product(
4595         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
4596 {
4597         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
4598                                             &pw_multi_aff_flat_range_product);
4599 }
4600
4601 /* If data->pma and *entry have the same domain space, then compute
4602  * their flat range product and the result to data->res.
4603  */
4604 static int flat_range_product_entry(void **entry, void *user)
4605 {
4606         struct isl_union_pw_multi_aff_bin_data *data = user;
4607         isl_pw_multi_aff *pma2 = *entry;
4608
4609         if (!isl_space_tuple_match(data->pma->dim, isl_dim_in,
4610                                  pma2->dim, isl_dim_in))
4611                 return 0;
4612
4613         pma2 = isl_pw_multi_aff_flat_range_product(
4614                                         isl_pw_multi_aff_copy(data->pma),
4615                                         isl_pw_multi_aff_copy(pma2));
4616
4617         data->res = isl_union_pw_multi_aff_add_pw_multi_aff(data->res, pma2);
4618
4619         return 0;
4620 }
4621
4622 /* Given two isl_union_pw_multi_affs A -> B and C -> D,
4623  * construct an isl_union_pw_multi_aff (A * C) -> (B, D).
4624  */
4625 __isl_give isl_union_pw_multi_aff *isl_union_pw_multi_aff_flat_range_product(
4626         __isl_take isl_union_pw_multi_aff *upma1,
4627         __isl_take isl_union_pw_multi_aff *upma2)
4628 {
4629         return bin_op(upma1, upma2, &flat_range_product_entry);
4630 }
4631
4632 /* Replace the affine expressions at position "pos" in "pma" by "pa".
4633  * The parameters are assumed to have been aligned.
4634  *
4635  * The implementation essentially performs an isl_pw_*_on_shared_domain,
4636  * except that it works on two different isl_pw_* types.
4637  */
4638 static __isl_give isl_pw_multi_aff *pw_multi_aff_set_pw_aff(
4639         __isl_take isl_pw_multi_aff *pma, unsigned pos,
4640         __isl_take isl_pw_aff *pa)
4641 {
4642         int i, j, n;
4643         isl_pw_multi_aff *res = NULL;
4644
4645         if (!pma || !pa)
4646                 goto error;
4647
4648         if (!isl_space_tuple_match(pma->dim, isl_dim_in, pa->dim, isl_dim_in))
4649                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
4650                         "domains don't match", goto error);
4651         if (pos >= isl_pw_multi_aff_dim(pma, isl_dim_out))
4652                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
4653                         "index out of bounds", goto error);
4654
4655         n = pma->n * pa->n;
4656         res = isl_pw_multi_aff_alloc_size(isl_pw_multi_aff_get_space(pma), n);
4657
4658         for (i = 0; i < pma->n; ++i) {
4659                 for (j = 0; j < pa->n; ++j) {
4660                         isl_set *common;
4661                         isl_multi_aff *res_ij;
4662                         int empty;
4663
4664                         common = isl_set_intersect(isl_set_copy(pma->p[i].set),
4665                                                    isl_set_copy(pa->p[j].set));
4666                         empty = isl_set_plain_is_empty(common);
4667                         if (empty < 0 || empty) {
4668                                 isl_set_free(common);
4669                                 if (empty < 0)
4670                                         goto error;
4671                                 continue;
4672                         }
4673
4674                         res_ij = isl_multi_aff_set_aff(
4675                                         isl_multi_aff_copy(pma->p[i].maff), pos,
4676                                         isl_aff_copy(pa->p[j].aff));
4677                         res_ij = isl_multi_aff_gist(res_ij,
4678                                         isl_set_copy(common));
4679
4680                         res = isl_pw_multi_aff_add_piece(res, common, res_ij);
4681                 }
4682         }
4683
4684         isl_pw_multi_aff_free(pma);
4685         isl_pw_aff_free(pa);
4686         return res;
4687 error:
4688         isl_pw_multi_aff_free(pma);
4689         isl_pw_aff_free(pa);
4690         return isl_pw_multi_aff_free(res);
4691 }
4692
4693 /* Replace the affine expressions at position "pos" in "pma" by "pa".
4694  */
4695 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_set_pw_aff(
4696         __isl_take isl_pw_multi_aff *pma, unsigned pos,
4697         __isl_take isl_pw_aff *pa)
4698 {
4699         if (!pma || !pa)
4700                 goto error;
4701         if (isl_space_match(pma->dim, isl_dim_param, pa->dim, isl_dim_param))
4702                 return pw_multi_aff_set_pw_aff(pma, pos, pa);
4703         if (!isl_space_has_named_params(pma->dim) ||
4704             !isl_space_has_named_params(pa->dim))
4705                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
4706                         "unaligned unnamed parameters", goto error);
4707         pma = isl_pw_multi_aff_align_params(pma, isl_pw_aff_get_space(pa));
4708         pa = isl_pw_aff_align_params(pa, isl_pw_multi_aff_get_space(pma));
4709         return pw_multi_aff_set_pw_aff(pma, pos, pa);
4710 error:
4711         isl_pw_multi_aff_free(pma);
4712         isl_pw_aff_free(pa);
4713         return NULL;
4714 }
4715
4716 /* Check that the domain space of "pa" matches "space".
4717  *
4718  * Return 0 on success and -1 on error.
4719  */
4720 int isl_pw_aff_check_match_domain_space(__isl_keep isl_pw_aff *pa,
4721         __isl_keep isl_space *space)
4722 {
4723         isl_space *pa_space;
4724         int match;
4725
4726         if (!pa || !space)
4727                 return -1;
4728
4729         pa_space = isl_pw_aff_get_space(pa);
4730
4731         match = isl_space_match(space, isl_dim_param, pa_space, isl_dim_param);
4732         if (match < 0)
4733                 goto error;
4734         if (!match)
4735                 isl_die(isl_pw_aff_get_ctx(pa), isl_error_invalid,
4736                         "parameters don't match", goto error);
4737         match = isl_space_tuple_match(space, isl_dim_in, pa_space, isl_dim_in);
4738         if (match < 0)
4739                 goto error;
4740         if (!match)
4741                 isl_die(isl_pw_aff_get_ctx(pa), isl_error_invalid,
4742                         "domains don't match", goto error);
4743         isl_space_free(pa_space);
4744         return 0;
4745 error:
4746         isl_space_free(pa_space);
4747         return -1;
4748 }
4749
4750 #undef BASE
4751 #define BASE pw_aff
4752
4753 #include <isl_multi_templ.c>
4754
4755 /* Scale the first elements of "ma" by the corresponding elements of "vec".
4756  */
4757 __isl_give isl_multi_aff *isl_multi_aff_scale_vec(__isl_take isl_multi_aff *ma,
4758         __isl_take isl_vec *vec)
4759 {
4760         int i, n;
4761         isl_int v;
4762
4763         if (!ma || !vec)
4764                 goto error;
4765
4766         n = isl_multi_aff_dim(ma, isl_dim_out);
4767         if (isl_vec_size(vec) < n)
4768                 n = isl_vec_size(vec);
4769
4770         isl_int_init(v);
4771         for (i = 0; i < n; ++i) {
4772                 isl_aff *aff;
4773
4774                 isl_vec_get_element(vec, i, &v);
4775
4776                 aff = isl_multi_aff_get_aff(ma, i);
4777                 aff = isl_aff_scale(aff, v);
4778                 ma = isl_multi_aff_set_aff(ma, i, aff);
4779         }
4780         isl_int_clear(v);
4781
4782         isl_vec_free(vec);
4783         return ma;
4784 error:
4785         isl_vec_free(vec);
4786         isl_multi_aff_free(ma);
4787         return NULL;
4788 }
4789
4790 /* Scale the first elements of "pma" by the corresponding elements of "vec".
4791  */
4792 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_scale_vec(
4793         __isl_take isl_pw_multi_aff *pma, __isl_take isl_vec *v)
4794 {
4795         int i;
4796
4797         pma = isl_pw_multi_aff_cow(pma);
4798         if (!pma || !v)
4799                 goto error;
4800
4801         for (i = 0; i < pma->n; ++i) {
4802                 pma->p[i].maff = isl_multi_aff_scale_vec(pma->p[i].maff,
4803                                                         isl_vec_copy(v));
4804                 if (!pma->p[i].maff)
4805                         goto error;
4806         }
4807
4808         isl_vec_free(v);
4809         return pma;
4810 error:
4811         isl_vec_free(v);
4812         isl_pw_multi_aff_free(pma);
4813         return NULL;
4814 }
4815
4816 /* This function is called for each entry of an isl_union_pw_multi_aff.
4817  * Replace the entry by the result of applying isl_pw_multi_aff_scale_vec
4818  * to the original entry with the isl_vec in "user" as extra argument.
4819  */
4820 static int union_pw_multi_aff_scale_vec_entry(void **entry, void *user)
4821 {
4822         isl_pw_multi_aff **pma = (isl_pw_multi_aff **) entry;
4823         isl_vec *v = user;
4824
4825         *pma = isl_pw_multi_aff_scale_vec(*pma, isl_vec_copy(v));
4826         if (!*pma)
4827                 return -1;
4828
4829         return 0;
4830 }
4831
4832 /* Scale the first elements of "upma" by the corresponding elements of "vec".
4833  */
4834 __isl_give isl_union_pw_multi_aff *isl_union_pw_multi_aff_scale_vec(
4835         __isl_take isl_union_pw_multi_aff *upma, __isl_take isl_vec *v)
4836 {
4837         upma = isl_union_pw_multi_aff_cow(upma);
4838         if (!upma || !v)
4839                 goto error;
4840
4841         if (isl_hash_table_foreach(upma->dim->ctx, &upma->table,
4842                                    &union_pw_multi_aff_scale_vec_entry, v) < 0)
4843                 goto error;
4844
4845         isl_vec_free(v);
4846         return upma;
4847 error:
4848         isl_vec_free(v);
4849         isl_union_pw_multi_aff_free(upma);
4850         return NULL;
4851 }