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