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