isl_aff.c

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