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