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