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