isl_aff.c

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