isl_aff.c

   1 /*
   2  * Copyright 2011      INRIA Saclay
   3  * Copyright 2011      Sven Verdoolaege
   4  * Copyright 2012      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_list_private.h>
  23 #include <isl/constraint.h>
  24 #include <isl/seq.h>
  25 #include <isl/set.h>
  26 #include <isl_config.h>
  27
  28 __isl_give isl_aff *isl_aff_alloc_vec(__isl_take isl_local_space *ls,
  29         __isl_take isl_vec *v)
  30 {
  31         isl_aff *aff;
  32
  33         if (!ls || !v)
  34                 goto error;
  35
  36         aff = isl_calloc_type(v->ctx, struct isl_aff);
  37         if (!aff)
  38                 goto error;
  39
  40         aff->ref = 1;
  41         aff->ls = ls;
  42         aff->v = v;
  43
  44         return aff;
  45 error:
  46         isl_local_space_free(ls);
  47         isl_vec_free(v);
  48         return NULL;
  49 }
  50
  51 __isl_give isl_aff *isl_aff_alloc(__isl_take isl_local_space *ls)
  52 {
  53         isl_ctx *ctx;
  54         isl_vec *v;
  55         unsigned total;
  56
  57         if (!ls)
  58                 return NULL;
  59
  60         ctx = isl_local_space_get_ctx(ls);
  61         if (!isl_local_space_divs_known(ls))
  62                 isl_die(ctx, isl_error_invalid, "local space has unknown divs",
  63                         goto error);
  64         if (!isl_local_space_is_set(ls))
  65                 isl_die(ctx, isl_error_invalid,
  66                         "domain of affine expression should be a set",
  67                         goto error);
  68
  69         total = isl_local_space_dim(ls, isl_dim_all);
  70         v = isl_vec_alloc(ctx, 1 + 1 + total);
  71         return isl_aff_alloc_vec(ls, v);
  72 error:
  73         isl_local_space_free(ls);
  74         return NULL;
  75 }
  76
  77 __isl_give isl_aff *isl_aff_zero_on_domain(__isl_take isl_local_space *ls)
  78 {
  79         isl_aff *aff;
  80
  81         aff = isl_aff_alloc(ls);
  82         if (!aff)
  83                 return NULL;
  84
  85         isl_int_set_si(aff->v->el[0], 1);
  86         isl_seq_clr(aff->v->el + 1, aff->v->size - 1);
  87
  88         return aff;
  89 }
  90
  91 /* Return an affine expression that is equal to the specified dimension
  92  * in "ls".
  93  */
  94 __isl_give isl_aff *isl_aff_var_on_domain(__isl_take isl_local_space *ls,
  95         enum isl_dim_type type, unsigned pos)
  96 {
  97         isl_space *space;
  98         isl_aff *aff;
  99
 100         if (!ls)
 101                 return NULL;
 102
 103         space = isl_local_space_get_space(ls);
 104         if (!space)
 105                 goto error;
 106         if (isl_space_is_map(space))
 107                 isl_die(isl_space_get_ctx(space), isl_error_invalid,
 108                         "expecting (parameter) set space", goto error);
 109         if (pos >= isl_local_space_dim(ls, type))
 110                 isl_die(isl_space_get_ctx(space), isl_error_invalid,
 111                         "position out of bounds", goto error);
 112
 113         isl_space_free(space);
 114         aff = isl_aff_alloc(ls);
 115         if (!aff)
 116                 return NULL;
 117
 118         pos += isl_local_space_offset(aff->ls, type);
 119
 120         isl_int_set_si(aff->v->el[0], 1);
 121         isl_seq_clr(aff->v->el + 1, aff->v->size - 1);
 122         isl_int_set_si(aff->v->el[1 + pos], 1);
 123
 124         return aff;
 125 error:
 126         isl_local_space_free(ls);
 127         isl_space_free(space);
 128         return NULL;
 129 }
 130
 131 __isl_give isl_aff *isl_aff_copy(__isl_keep isl_aff *aff)
 132 {
 133         if (!aff)
 134                 return NULL;
 135
 136         aff->ref++;
 137         return aff;
 138 }
 139
 140 __isl_give isl_aff *isl_aff_dup(__isl_keep isl_aff *aff)
 141 {
 142         if (!aff)
 143                 return NULL;
 144
 145         return isl_aff_alloc_vec(isl_local_space_copy(aff->ls),
 146                                  isl_vec_copy(aff->v));
 147 }
 148
 149 __isl_give isl_aff *isl_aff_cow(__isl_take isl_aff *aff)
 150 {
 151         if (!aff)
 152                 return NULL;
 153
 154         if (aff->ref == 1)
 155                 return aff;
 156         aff->ref--;
 157         return isl_aff_dup(aff);
 158 }
 159
 160 void *isl_aff_free(__isl_take isl_aff *aff)
 161 {
 162         if (!aff)
 163                 return NULL;
 164
 165         if (--aff->ref > 0)
 166                 return NULL;
 167
 168         isl_local_space_free(aff->ls);
 169         isl_vec_free(aff->v);
 170
 171         free(aff);
 172
 173         return NULL;
 174 }
 175
 176 isl_ctx *isl_aff_get_ctx(__isl_keep isl_aff *aff)
 177 {
 178         return aff ? isl_local_space_get_ctx(aff->ls) : NULL;
 179 }
 180
 181 /* Externally, an isl_aff has a map space, but internally, the
 182  * ls field corresponds to the domain of that space.
 183  */
 184 int isl_aff_dim(__isl_keep isl_aff *aff, enum isl_dim_type type)
 185 {
 186         if (!aff)
 187                 return 0;
 188         if (type == isl_dim_out)
 189                 return 1;
 190         if (type == isl_dim_in)
 191                 type = isl_dim_set;
 192         return isl_local_space_dim(aff->ls, type);
 193 }
 194
 195 __isl_give isl_space *isl_aff_get_domain_space(__isl_keep isl_aff *aff)
 196 {
 197         return aff ? isl_local_space_get_space(aff->ls) : NULL;
 198 }
 199
 200 __isl_give isl_space *isl_aff_get_space(__isl_keep isl_aff *aff)
 201 {
 202         isl_space *space;
 203         if (!aff)
 204                 return NULL;
 205         space = isl_local_space_get_space(aff->ls);
 206         space = isl_space_from_domain(space);
 207         space = isl_space_add_dims(space, isl_dim_out, 1);
 208         return space;
 209 }
 210
 211 __isl_give isl_local_space *isl_aff_get_domain_local_space(
 212         __isl_keep isl_aff *aff)
 213 {
 214         return aff ? isl_local_space_copy(aff->ls) : NULL;
 215 }
 216
 217 __isl_give isl_local_space *isl_aff_get_local_space(__isl_keep isl_aff *aff)
 218 {
 219         isl_local_space *ls;
 220         if (!aff)
 221                 return NULL;
 222         ls = isl_local_space_copy(aff->ls);
 223         ls = isl_local_space_from_domain(ls);
 224         ls = isl_local_space_add_dims(ls, isl_dim_out, 1);
 225         return ls;
 226 }
 227
 228 /* Externally, an isl_aff has a map space, but internally, the
 229  * ls field corresponds to the domain of that space.
 230  */
 231 const char *isl_aff_get_dim_name(__isl_keep isl_aff *aff,
 232         enum isl_dim_type type, unsigned pos)
 233 {
 234         if (!aff)
 235                 return NULL;
 236         if (type == isl_dim_out)
 237                 return NULL;
 238         if (type == isl_dim_in)
 239                 type = isl_dim_set;
 240         return isl_local_space_get_dim_name(aff->ls, type, pos);
 241 }
 242
 243 __isl_give isl_aff *isl_aff_reset_domain_space(__isl_take isl_aff *aff,
 244         __isl_take isl_space *dim)
 245 {
 246         aff = isl_aff_cow(aff);
 247         if (!aff || !dim)
 248                 goto error;
 249
 250         aff->ls = isl_local_space_reset_space(aff->ls, dim);
 251         if (!aff->ls)
 252                 return isl_aff_free(aff);
 253
 254         return aff;
 255 error:
 256         isl_aff_free(aff);
 257         isl_space_free(dim);
 258         return NULL;
 259 }
 260
 261 /* Reset the space of "aff".  This function is called from isl_pw_templ.c
 262  * and doesn't know if the space of an element object is represented
 263  * directly or through its domain.  It therefore passes along both.
 264  */
 265 __isl_give isl_aff *isl_aff_reset_space_and_domain(__isl_take isl_aff *aff,
 266         __isl_take isl_space *space, __isl_take isl_space *domain)
 267 {
 268         isl_space_free(space);
 269         return isl_aff_reset_domain_space(aff, domain);
 270 }
 271
 272 /* Reorder the coefficients of the affine expression based
 273  * on the given reodering.
 274  * The reordering r is assumed to have been extended with the local
 275  * variables.
 276  */
 277 static __isl_give isl_vec *vec_reorder(__isl_take isl_vec *vec,
 278         __isl_take isl_reordering *r, int n_div)
 279 {
 280         isl_vec *res;
 281         int i;
 282
 283         if (!vec || !r)
 284                 goto error;
 285
 286         res = isl_vec_alloc(vec->ctx,
 287                             2 + isl_space_dim(r->dim, isl_dim_all) + n_div);
 288         isl_seq_cpy(res->el, vec->el, 2);
 289         isl_seq_clr(res->el + 2, res->size - 2);
 290         for (i = 0; i < r->len; ++i)
 291                 isl_int_set(res->el[2 + r->pos[i]], vec->el[2 + i]);
 292
 293         isl_reordering_free(r);
 294         isl_vec_free(vec);
 295         return res;
 296 error:
 297         isl_vec_free(vec);
 298         isl_reordering_free(r);
 299         return NULL;
 300 }
 301
 302 /* Reorder the dimensions of the domain of "aff" according
 303  * to the given reordering.
 304  */
 305 __isl_give isl_aff *isl_aff_realign_domain(__isl_take isl_aff *aff,
 306         __isl_take isl_reordering *r)
 307 {
 308         aff = isl_aff_cow(aff);
 309         if (!aff)
 310                 goto error;
 311
 312         r = isl_reordering_extend(r, aff->ls->div->n_row);
 313         aff->v = vec_reorder(aff->v, isl_reordering_copy(r),
 314                                 aff->ls->div->n_row);
 315         aff->ls = isl_local_space_realign(aff->ls, r);
 316
 317         if (!aff->v || !aff->ls)
 318                 return isl_aff_free(aff);
 319
 320         return aff;
 321 error:
 322         isl_aff_free(aff);
 323         isl_reordering_free(r);
 324         return NULL;
 325 }
 326
 327 __isl_give isl_aff *isl_aff_align_params(__isl_take isl_aff *aff,
 328         __isl_take isl_space *model)
 329 {
 330         if (!aff || !model)
 331                 goto error;
 332
 333         if (!isl_space_match(aff->ls->dim, isl_dim_param,
 334                              model, isl_dim_param)) {
 335                 isl_reordering *exp;
 336
 337                 model = isl_space_drop_dims(model, isl_dim_in,
 338                                         0, isl_space_dim(model, isl_dim_in));
 339                 model = isl_space_drop_dims(model, isl_dim_out,
 340                                         0, isl_space_dim(model, isl_dim_out));
 341                 exp = isl_parameter_alignment_reordering(aff->ls->dim, model);
 342                 exp = isl_reordering_extend_space(exp,
 343                                         isl_aff_get_domain_space(aff));
 344                 aff = isl_aff_realign_domain(aff, exp);
 345         }
 346
 347         isl_space_free(model);
 348         return aff;
 349 error:
 350         isl_space_free(model);
 351         isl_aff_free(aff);
 352         return NULL;
 353 }
 354
 355 int isl_aff_plain_is_zero(__isl_keep isl_aff *aff)
 356 {
 357         if (!aff)
 358                 return -1;
 359
 360         return isl_seq_first_non_zero(aff->v->el + 1, aff->v->size - 1) < 0;
 361 }
 362
 363 int isl_aff_plain_is_equal(__isl_keep isl_aff *aff1, __isl_keep isl_aff *aff2)
 364 {
 365         int equal;
 366
 367         if (!aff1 || !aff2)
 368                 return -1;
 369
 370         equal = isl_local_space_is_equal(aff1->ls, aff2->ls);
 371         if (equal < 0 || !equal)
 372                 return equal;
 373
 374         return isl_vec_is_equal(aff1->v, aff2->v);
 375 }
 376
 377 int isl_aff_get_denominator(__isl_keep isl_aff *aff, isl_int *v)
 378 {
 379         if (!aff)
 380                 return -1;
 381         isl_int_set(*v, aff->v->el[0]);
 382         return 0;
 383 }
 384
 385 int isl_aff_get_constant(__isl_keep isl_aff *aff, isl_int *v)
 386 {
 387         if (!aff)
 388                 return -1;
 389         isl_int_set(*v, aff->v->el[1]);
 390         return 0;
 391 }
 392
 393 int isl_aff_get_coefficient(__isl_keep isl_aff *aff,
 394         enum isl_dim_type type, int pos, isl_int *v)
 395 {
 396         if (!aff)
 397                 return -1;
 398
 399         if (type == isl_dim_out)
 400                 isl_die(aff->v->ctx, isl_error_invalid,
 401                         "output/set dimension does not have a coefficient",
 402                         return -1);
 403         if (type == isl_dim_in)
 404                 type = isl_dim_set;
 405
 406         if (pos >= isl_local_space_dim(aff->ls, type))
 407                 isl_die(aff->v->ctx, isl_error_invalid,
 408                         "position out of bounds", return -1);
 409
 410         pos += isl_local_space_offset(aff->ls, type);
 411         isl_int_set(*v, aff->v->el[1 + pos]);
 412
 413         return 0;
 414 }
 415
 416 __isl_give isl_aff *isl_aff_set_denominator(__isl_take isl_aff *aff, isl_int v)
 417 {
 418         aff = isl_aff_cow(aff);
 419         if (!aff)
 420                 return NULL;
 421
 422         aff->v = isl_vec_cow(aff->v);
 423         if (!aff->v)
 424                 return isl_aff_free(aff);
 425
 426         isl_int_set(aff->v->el[0], v);
 427
 428         return aff;
 429 }
 430
 431 __isl_give isl_aff *isl_aff_set_constant(__isl_take isl_aff *aff, isl_int v)
 432 {
 433         aff = isl_aff_cow(aff);
 434         if (!aff)
 435                 return NULL;
 436
 437         aff->v = isl_vec_cow(aff->v);
 438         if (!aff->v)
 439                 return isl_aff_free(aff);
 440
 441         isl_int_set(aff->v->el[1], v);
 442
 443         return aff;
 444 }
 445
 446 __isl_give isl_aff *isl_aff_add_constant(__isl_take isl_aff *aff, isl_int v)
 447 {
 448         if (isl_int_is_zero(v))
 449                 return aff;
 450
 451         aff = isl_aff_cow(aff);
 452         if (!aff)
 453                 return NULL;
 454
 455         aff->v = isl_vec_cow(aff->v);
 456         if (!aff->v)
 457                 return isl_aff_free(aff);
 458
 459         isl_int_addmul(aff->v->el[1], aff->v->el[0], v);
 460
 461         return aff;
 462 }
 463
 464 __isl_give isl_aff *isl_aff_add_constant_si(__isl_take isl_aff *aff, int v)
 465 {
 466         isl_int t;
 467
 468         isl_int_init(t);
 469         isl_int_set_si(t, v);
 470         aff = isl_aff_add_constant(aff, t);
 471         isl_int_clear(t);
 472
 473         return aff;
 474 }
 475
 476 /* Add "v" to the numerator of the constant term of "aff".
 477  */
 478 __isl_give isl_aff *isl_aff_add_constant_num(__isl_take isl_aff *aff, isl_int v)
 479 {
 480         if (isl_int_is_zero(v))
 481                 return aff;
 482
 483         aff = isl_aff_cow(aff);
 484         if (!aff)
 485                 return NULL;
 486
 487         aff->v = isl_vec_cow(aff->v);
 488         if (!aff->v)
 489                 return isl_aff_free(aff);
 490
 491         isl_int_add(aff->v->el[1], aff->v->el[1], v);
 492
 493         return aff;
 494 }
 495
 496 /* Add "v" to the numerator of the constant term of "aff".
 497  */
 498 __isl_give isl_aff *isl_aff_add_constant_num_si(__isl_take isl_aff *aff, int v)
 499 {
 500         isl_int t;
 501
 502         if (v == 0)
 503                 return aff;
 504
 505         isl_int_init(t);
 506         isl_int_set_si(t, v);
 507         aff = isl_aff_add_constant_num(aff, t);
 508         isl_int_clear(t);
 509
 510         return aff;
 511 }
 512
 513 __isl_give isl_aff *isl_aff_set_constant_si(__isl_take isl_aff *aff, int v)
 514 {
 515         aff = isl_aff_cow(aff);
 516         if (!aff)
 517                 return NULL;
 518
 519         aff->v = isl_vec_cow(aff->v);
 520         if (!aff->v)
 521                 return isl_aff_free(aff);
 522
 523         isl_int_set_si(aff->v->el[1], v);
 524
 525         return aff;
 526 }
 527
 528 __isl_give isl_aff *isl_aff_set_coefficient(__isl_take isl_aff *aff,
 529         enum isl_dim_type type, int pos, isl_int v)
 530 {
 531         if (!aff)
 532                 return NULL;
 533
 534         if (type == isl_dim_out)
 535                 isl_die(aff->v->ctx, isl_error_invalid,
 536                         "output/set dimension does not have a coefficient",
 537                         return isl_aff_free(aff));
 538         if (type == isl_dim_in)
 539                 type = isl_dim_set;
 540
 541         if (pos >= isl_local_space_dim(aff->ls, type))
 542                 isl_die(aff->v->ctx, isl_error_invalid,
 543                         "position out of bounds", return isl_aff_free(aff));
 544
 545         aff = isl_aff_cow(aff);
 546         if (!aff)
 547                 return NULL;
 548
 549         aff->v = isl_vec_cow(aff->v);
 550         if (!aff->v)
 551                 return isl_aff_free(aff);
 552
 553         pos += isl_local_space_offset(aff->ls, type);
 554         isl_int_set(aff->v->el[1 + pos], v);
 555
 556         return aff;
 557 }
 558
 559 __isl_give isl_aff *isl_aff_set_coefficient_si(__isl_take isl_aff *aff,
 560         enum isl_dim_type type, int pos, int v)
 561 {
 562         if (!aff)
 563                 return NULL;
 564
 565         if (type == isl_dim_out)
 566                 isl_die(aff->v->ctx, isl_error_invalid,
 567                         "output/set dimension does not have a coefficient",
 568                         return isl_aff_free(aff));
 569         if (type == isl_dim_in)
 570                 type = isl_dim_set;
 571
 572         if (pos >= isl_local_space_dim(aff->ls, type))
 573                 isl_die(aff->v->ctx, isl_error_invalid,
 574                         "position out of bounds", return isl_aff_free(aff));
 575
 576         aff = isl_aff_cow(aff);
 577         if (!aff)
 578                 return NULL;
 579
 580         aff->v = isl_vec_cow(aff->v);
 581         if (!aff->v)
 582                 return isl_aff_free(aff);
 583
 584         pos += isl_local_space_offset(aff->ls, type);
 585         isl_int_set_si(aff->v->el[1 + pos], v);
 586
 587         return aff;
 588 }
 589
 590 __isl_give isl_aff *isl_aff_add_coefficient(__isl_take isl_aff *aff,
 591         enum isl_dim_type type, int pos, isl_int v)
 592 {
 593         if (!aff)
 594                 return NULL;
 595
 596         if (type == isl_dim_out)
 597                 isl_die(aff->v->ctx, isl_error_invalid,
 598                         "output/set dimension does not have a coefficient",
 599                         return isl_aff_free(aff));
 600         if (type == isl_dim_in)
 601                 type = isl_dim_set;
 602
 603         if (pos >= isl_local_space_dim(aff->ls, type))
 604                 isl_die(aff->v->ctx, isl_error_invalid,
 605                         "position out of bounds", return isl_aff_free(aff));
 606
 607         aff = isl_aff_cow(aff);
 608         if (!aff)
 609                 return NULL;
 610
 611         aff->v = isl_vec_cow(aff->v);
 612         if (!aff->v)
 613                 return isl_aff_free(aff);
 614
 615         pos += isl_local_space_offset(aff->ls, type);
 616         isl_int_addmul(aff->v->el[1 + pos], aff->v->el[0], v);
 617
 618         return aff;
 619 }
 620
 621 __isl_give isl_aff *isl_aff_add_coefficient_si(__isl_take isl_aff *aff,
 622         enum isl_dim_type type, int pos, int v)
 623 {
 624         isl_int t;
 625
 626         isl_int_init(t);
 627         isl_int_set_si(t, v);
 628         aff = isl_aff_add_coefficient(aff, type, pos, t);
 629         isl_int_clear(t);
 630
 631         return aff;
 632 }
 633
 634 __isl_give isl_aff *isl_aff_get_div(__isl_keep isl_aff *aff, int pos)
 635 {
 636         if (!aff)
 637                 return NULL;
 638
 639         return isl_local_space_get_div(aff->ls, pos);
 640 }
 641
 642 __isl_give isl_aff *isl_aff_neg(__isl_take isl_aff *aff)
 643 {
 644         aff = isl_aff_cow(aff);
 645         if (!aff)
 646                 return NULL;
 647         aff->v = isl_vec_cow(aff->v);
 648         if (!aff->v)
 649                 return isl_aff_free(aff);
 650
 651         isl_seq_neg(aff->v->el + 1, aff->v->el + 1, aff->v->size - 1);
 652
 653         return aff;
 654 }
 655
 656 /* Remove divs from the local space that do not appear in the affine
 657  * expression.
 658  * We currently only remove divs at the end.
 659  * Some intermediate divs may also not appear directly in the affine
 660  * expression, but we would also need to check that no other divs are
 661  * defined in terms of them.
 662  */
 663 __isl_give isl_aff *isl_aff_remove_unused_divs( __isl_take isl_aff *aff)
 664 {
 665         int pos;
 666         int off;
 667         int n;
 668
 669         if (!aff)
 670                 return NULL;
 671
 672         n = isl_local_space_dim(aff->ls, isl_dim_div);
 673         off = isl_local_space_offset(aff->ls, isl_dim_div);
 674
 675         pos = isl_seq_last_non_zero(aff->v->el + 1 + off, n) + 1;
 676         if (pos == n)
 677                 return aff;
 678
 679         aff = isl_aff_cow(aff);
 680         if (!aff)
 681                 return NULL;
 682
 683         aff->ls = isl_local_space_drop_dims(aff->ls, isl_dim_div, pos, n - pos);
 684         aff->v = isl_vec_drop_els(aff->v, 1 + off + pos, n - pos);
 685         if (!aff->ls || !aff->v)
 686                 return isl_aff_free(aff);
 687
 688         return aff;
 689 }
 690
 691 /* Given two affine expressions "p" of length p_len (including the
 692  * denominator and the constant term) and "subs" of length subs_len,
 693  * plug in "subs" for the variable at position "pos".
 694  * The variables of "subs" and "p" are assumed to match up to subs_len,
 695  * but "p" may have additional variables.
 696  * "v" is an initialized isl_int that can be used internally.
 697  *
 698  * In particular, if "p" represents the expression
 699  *
 700  *      (a i + g)/m
 701  *
 702  * with i the variable at position "pos" and "subs" represents the expression
 703  *
 704  *      f/d
 705  *
 706  * then the result represents the expression
 707  *
 708  *      (a f + d g)/(m d)
 709  *
 710  */
 711 void isl_seq_substitute(isl_int *p, int pos, isl_int *subs,
 712         int p_len, int subs_len, isl_int v)
 713 {
 714         isl_int_set(v, p[1 + pos]);
 715         isl_int_set_si(p[1 + pos], 0);
 716         isl_seq_combine(p + 1, subs[0], p + 1, v, subs + 1, subs_len - 1);
 717         isl_seq_scale(p + subs_len, p + subs_len, subs[0], p_len - subs_len);
 718         isl_int_mul(p[0], p[0], subs[0]);
 719 }
 720
 721 /* Look for any divs in the aff->ls with a denominator equal to one
 722  * and plug them into the affine expression and any subsequent divs
 723  * that may reference the div.
 724  */
 725 static __isl_give isl_aff *plug_in_integral_divs(__isl_take isl_aff *aff)
 726 {
 727         int i, n;
 728         int len;
 729         isl_int v;
 730         isl_vec *vec;
 731         isl_local_space *ls;
 732         unsigned pos;
 733
 734         if (!aff)
 735                 return NULL;
 736
 737         n = isl_local_space_dim(aff->ls, isl_dim_div);
 738         len = aff->v->size;
 739         for (i = 0; i < n; ++i) {
 740                 if (!isl_int_is_one(aff->ls->div->row[i][0]))
 741                         continue;
 742                 ls = isl_local_space_copy(aff->ls);
 743                 ls = isl_local_space_substitute_seq(ls, isl_dim_div, i,
 744                                             aff->ls->div->row[i], len, i + 1);
 745                 vec = isl_vec_copy(aff->v);
 746                 vec = isl_vec_cow(vec);
 747                 if (!ls || !vec)
 748                         goto error;
 749
 750                 isl_int_init(v);
 751
 752                 pos = isl_local_space_offset(aff->ls, isl_dim_div) + i;
 753                 isl_seq_substitute(vec->el, pos, aff->ls->div->row[i],
 754                                         len, len, v);
 755
 756                 isl_int_clear(v);
 757
 758                 isl_vec_free(aff->v);
 759                 aff->v = vec;
 760                 isl_local_space_free(aff->ls);
 761                 aff->ls = ls;
 762         }
 763
 764         return aff;
 765 error:
 766         isl_vec_free(vec);
 767         isl_local_space_free(ls);
 768         return isl_aff_free(aff);
 769 }
 770
 771 /* Swap divs "a" and "b" in "aff", which is assumed to be non-NULL.
 772  *
 773  * Even though this function is only called on isl_affs with a single
 774  * reference, we are careful to only change aff->v and aff->ls together.
 775  */
 776 static __isl_give isl_aff *swap_div(__isl_take isl_aff *aff, int a, int b)
 777 {
 778         unsigned off = isl_local_space_offset(aff->ls, isl_dim_div);
 779         isl_local_space *ls;
 780         isl_vec *v;
 781
 782         ls = isl_local_space_copy(aff->ls);
 783         ls = isl_local_space_swap_div(ls, a, b);
 784         v = isl_vec_copy(aff->v);
 785         v = isl_vec_cow(v);
 786         if (!ls || !v)
 787                 goto error;
 788
 789         isl_int_swap(v->el[1 + off + a], v->el[1 + off + b]);
 790         isl_vec_free(aff->v);
 791         aff->v = v;
 792         isl_local_space_free(aff->ls);
 793         aff->ls = ls;
 794
 795         return aff;
 796 error:
 797         isl_vec_free(v);
 798         isl_local_space_free(ls);
 799         return isl_aff_free(aff);
 800 }
 801
 802 /* Merge divs "a" and "b" in "aff", which is assumed to be non-NULL.
 803  *
 804  * We currently do not actually remove div "b", but simply add its
 805  * coefficient to that of "a" and then zero it out.
 806  */
 807 static __isl_give isl_aff *merge_divs(__isl_take isl_aff *aff, int a, int b)
 808 {
 809         unsigned off = isl_local_space_offset(aff->ls, isl_dim_div);
 810
 811         if (isl_int_is_zero(aff->v->el[1 + off + b]))
 812                 return aff;
 813
 814         aff->v = isl_vec_cow(aff->v);
 815         if (!aff->v)
 816                 return isl_aff_free(aff);
 817
 818         isl_int_add(aff->v->el[1 + off + a],
 819                     aff->v->el[1 + off + a], aff->v->el[1 + off + b]);
 820         isl_int_set_si(aff->v->el[1 + off + b], 0);
 821
 822         return aff;
 823 }
 824
 825 /* Sort the divs in the local space of "aff" according to
 826  * the comparison function "cmp_row" in isl_local_space.c,
 827  * combining the coefficients of identical divs.
 828  *
 829  * Reordering divs does not change the semantics of "aff",
 830  * so there is no need to call isl_aff_cow.
 831  * Moreover, this function is currently only called on isl_affs
 832  * with a single reference.
 833  */
 834 static __isl_give isl_aff *sort_divs(__isl_take isl_aff *aff)
 835 {
 836         int i, j, n;
 837         unsigned off;
 838
 839         if (!aff)
 840                 return NULL;
 841
 842         off = isl_local_space_offset(aff->ls, isl_dim_div);
 843         n = isl_aff_dim(aff, isl_dim_div);
 844         for (i = 1; i < n; ++i) {
 845                 for (j = i - 1; j >= 0; --j) {
 846                         int cmp = isl_mat_cmp_div(aff->ls->div, j, j + 1);
 847                         if (cmp < 0)
 848                                 break;
 849                         if (cmp == 0)
 850                                 aff = merge_divs(aff, j, j + 1);
 851                         else
 852                                 aff = swap_div(aff, j, j + 1);
 853                         if (!aff)
 854                                 return NULL;
 855                 }
 856         }
 857
 858         return aff;
 859 }
 860
 861 /* Normalize the representation of "aff".
 862  *
 863  * This function should only be called of "new" isl_affs, i.e.,
 864  * with only a single reference.  We therefore do not need to
 865  * worry about affecting other instances.
 866  */
 867 __isl_give isl_aff *isl_aff_normalize(__isl_take isl_aff *aff)
 868 {
 869         if (!aff)
 870                 return NULL;
 871         aff->v = isl_vec_normalize(aff->v);
 872         if (!aff->v)
 873                 return isl_aff_free(aff);
 874         aff = plug_in_integral_divs(aff);
 875         aff = sort_divs(aff);
 876         aff = isl_aff_remove_unused_divs(aff);
 877         return aff;
 878 }
 879
 880 /* Given f, return floor(f).
 881  * If f is an integer expression, then just return f.
 882  * If f is a constant, then return the constant floor(f).
 883  * Otherwise, if f = g/m, write g = q m + r,
 884  * create a new div d = [r/m] and return the expression q + d.
 885  * The coefficients in r are taken to lie between -m/2 and m/2.
 886  */
 887 __isl_give isl_aff *isl_aff_floor(__isl_take isl_aff *aff)
 888 {
 889         int i;
 890         int size;
 891         isl_ctx *ctx;
 892         isl_vec *div;
 893
 894         if (!aff)
 895                 return NULL;
 896
 897         if (isl_int_is_one(aff->v->el[0]))
 898                 return aff;
 899
 900         aff = isl_aff_cow(aff);
 901         if (!aff)
 902                 return NULL;
 903
 904         aff->v = isl_vec_cow(aff->v);
 905         if (!aff->v)
 906                 return isl_aff_free(aff);
 907
 908         if (isl_aff_is_cst(aff)) {
 909                 isl_int_fdiv_q(aff->v->el[1], aff->v->el[1], aff->v->el[0]);
 910                 isl_int_set_si(aff->v->el[0], 1);
 911                 return aff;
 912         }
 913
 914         div = isl_vec_copy(aff->v);
 915         div = isl_vec_cow(div);
 916         if (!div)
 917                 return isl_aff_free(aff);
 918
 919         ctx = isl_aff_get_ctx(aff);
 920         isl_int_fdiv_q(aff->v->el[0], aff->v->el[0], ctx->two);
 921         for (i = 1; i < aff->v->size; ++i) {
 922                 isl_int_fdiv_r(div->el[i], div->el[i], div->el[0]);
 923                 isl_int_fdiv_q(aff->v->el[i], aff->v->el[i], div->el[0]);
 924                 if (isl_int_gt(div->el[i], aff->v->el[0])) {
 925                         isl_int_sub(div->el[i], div->el[i], div->el[0]);
 926                         isl_int_add_ui(aff->v->el[i], aff->v->el[i], 1);
 927                 }
 928         }
 929
 930         aff->ls = isl_local_space_add_div(aff->ls, div);
 931         if (!aff->ls)
 932                 return isl_aff_free(aff);
 933
 934         size = aff->v->size;
 935         aff->v = isl_vec_extend(aff->v, size + 1);
 936         if (!aff->v)
 937                 return isl_aff_free(aff);
 938         isl_int_set_si(aff->v->el[0], 1);
 939         isl_int_set_si(aff->v->el[size], 1);
 940
 941         return aff;
 942 }
 943
 944 /* Compute
 945  *
 946  *      aff mod m = aff - m * floor(aff/m)
 947  */
 948 __isl_give isl_aff *isl_aff_mod(__isl_take isl_aff *aff, isl_int m)
 949 {
 950         isl_aff *res;
 951
 952         res = isl_aff_copy(aff);
 953         aff = isl_aff_scale_down(aff, m);
 954         aff = isl_aff_floor(aff);
 955         aff = isl_aff_scale(aff, m);
 956         res = isl_aff_sub(res, aff);
 957
 958         return res;
 959 }
 960
 961 /* Compute
 962  *
 963  *      pwaff mod m = pwaff - m * floor(pwaff/m)
 964  */
 965 __isl_give isl_pw_aff *isl_pw_aff_mod(__isl_take isl_pw_aff *pwaff, isl_int m)
 966 {
 967         isl_pw_aff *res;
 968
 969         res = isl_pw_aff_copy(pwaff);
 970         pwaff = isl_pw_aff_scale_down(pwaff, m);
 971         pwaff = isl_pw_aff_floor(pwaff);
 972         pwaff = isl_pw_aff_scale(pwaff, m);
 973         res = isl_pw_aff_sub(res, pwaff);
 974
 975         return res;
 976 }
 977
 978 /* Given f, return ceil(f).
 979  * If f is an integer expression, then just return f.
 980  * Otherwise, create a new div d = [-f] and return the expression -d.
 981  */
 982 __isl_give isl_aff *isl_aff_ceil(__isl_take isl_aff *aff)
 983 {
 984         if (!aff)
 985                 return NULL;
 986
 987         if (isl_int_is_one(aff->v->el[0]))
 988                 return aff;
 989
 990         aff = isl_aff_neg(aff);
 991         aff = isl_aff_floor(aff);
 992         aff = isl_aff_neg(aff);
 993
 994         return aff;
 995 }
 996
 997 /* Apply the expansion computed by isl_merge_divs.
 998  * The expansion itself is given by "exp" while the resulting
 999  * list of divs is given by "div".
1000  */
1001 __isl_give isl_aff *isl_aff_expand_divs( __isl_take isl_aff *aff,
1002         __isl_take isl_mat *div, int *exp)
1003 {
1004         int i, j;
1005         int old_n_div;
1006         int new_n_div;
1007         int offset;
1008
1009         aff = isl_aff_cow(aff);
1010         if (!aff || !div)
1011                 goto error;
1012
1013         old_n_div = isl_local_space_dim(aff->ls, isl_dim_div);
1014         new_n_div = isl_mat_rows(div);
1015         if (new_n_div < old_n_div)
1016                 isl_die(isl_mat_get_ctx(div), isl_error_invalid,
1017                         "not an expansion", goto error);
1018
1019         aff->v = isl_vec_extend(aff->v, aff->v->size + new_n_div - old_n_div);
1020         if (!aff->v)
1021                 goto error;
1022
1023         offset = 1 + isl_local_space_offset(aff->ls, isl_dim_div);
1024         j = old_n_div - 1;
1025         for (i = new_n_div - 1; i >= 0; --i) {
1026                 if (j >= 0 && exp[j] == i) {
1027                         if (i != j)
1028                                 isl_int_swap(aff->v->el[offset + i],
1029                                              aff->v->el[offset + j]);
1030                         j--;
1031                 } else
1032                         isl_int_set_si(aff->v->el[offset + i], 0);
1033         }
1034
1035         aff->ls = isl_local_space_replace_divs(aff->ls, isl_mat_copy(div));
1036         if (!aff->ls)
1037                 goto error;
1038         isl_mat_free(div);
1039         return aff;
1040 error:
1041         isl_aff_free(aff);
1042         isl_mat_free(div);
1043         return NULL;
1044 }
1045
1046 /* Add two affine expressions that live in the same local space.
1047  */
1048 static __isl_give isl_aff *add_expanded(__isl_take isl_aff *aff1,
1049         __isl_take isl_aff *aff2)
1050 {
1051         isl_int gcd, f;
1052
1053         aff1 = isl_aff_cow(aff1);
1054         if (!aff1 || !aff2)
1055                 goto error;
1056
1057         aff1->v = isl_vec_cow(aff1->v);
1058         if (!aff1->v)
1059                 goto error;
1060
1061         isl_int_init(gcd);
1062         isl_int_init(f);
1063         isl_int_gcd(gcd, aff1->v->el[0], aff2->v->el[0]);
1064         isl_int_divexact(f, aff2->v->el[0], gcd);
1065         isl_seq_scale(aff1->v->el + 1, aff1->v->el + 1, f, aff1->v->size - 1);
1066         isl_int_divexact(f, aff1->v->el[0], gcd);
1067         isl_seq_addmul(aff1->v->el + 1, f, aff2->v->el + 1, aff1->v->size - 1);
1068         isl_int_divexact(f, aff2->v->el[0], gcd);
1069         isl_int_mul(aff1->v->el[0], aff1->v->el[0], f);
1070         isl_int_clear(f);
1071         isl_int_clear(gcd);
1072
1073         isl_aff_free(aff2);
1074         return aff1;
1075 error:
1076         isl_aff_free(aff1);
1077         isl_aff_free(aff2);
1078         return NULL;
1079 }
1080
1081 __isl_give isl_aff *isl_aff_add(__isl_take isl_aff *aff1,
1082         __isl_take isl_aff *aff2)
1083 {
1084         isl_ctx *ctx;
1085         int *exp1 = NULL;
1086         int *exp2 = NULL;
1087         isl_mat *div;
1088
1089         if (!aff1 || !aff2)
1090                 goto error;
1091
1092         ctx = isl_aff_get_ctx(aff1);
1093         if (!isl_space_is_equal(aff1->ls->dim, aff2->ls->dim))
1094                 isl_die(ctx, isl_error_invalid,
1095                         "spaces don't match", goto error);
1096
1097         if (aff1->ls->div->n_row == 0 && aff2->ls->div->n_row == 0)
1098                 return add_expanded(aff1, aff2);
1099
1100         exp1 = isl_alloc_array(ctx, int, aff1->ls->div->n_row);
1101         exp2 = isl_alloc_array(ctx, int, aff2->ls->div->n_row);
1102         if (!exp1 || !exp2)
1103                 goto error;
1104
1105         div = isl_merge_divs(aff1->ls->div, aff2->ls->div, exp1, exp2);
1106         aff1 = isl_aff_expand_divs(aff1, isl_mat_copy(div), exp1);
1107         aff2 = isl_aff_expand_divs(aff2, div, exp2);
1108         free(exp1);
1109         free(exp2);
1110
1111         return add_expanded(aff1, aff2);
1112 error:
1113         free(exp1);
1114         free(exp2);
1115         isl_aff_free(aff1);
1116         isl_aff_free(aff2);
1117         return NULL;
1118 }
1119
1120 __isl_give isl_aff *isl_aff_sub(__isl_take isl_aff *aff1,
1121         __isl_take isl_aff *aff2)
1122 {
1123         return isl_aff_add(aff1, isl_aff_neg(aff2));
1124 }
1125
1126 __isl_give isl_aff *isl_aff_scale(__isl_take isl_aff *aff, isl_int f)
1127 {
1128         isl_int gcd;
1129
1130         if (isl_int_is_one(f))
1131                 return aff;
1132
1133         aff = isl_aff_cow(aff);
1134         if (!aff)
1135                 return NULL;
1136         aff->v = isl_vec_cow(aff->v);
1137         if (!aff->v)
1138                 return isl_aff_free(aff);
1139
1140         isl_int_init(gcd);
1141         isl_int_gcd(gcd, aff->v->el[0], f);
1142         isl_int_divexact(aff->v->el[0], aff->v->el[0], gcd);
1143         isl_int_divexact(gcd, f, gcd);
1144         isl_seq_scale(aff->v->el + 1, aff->v->el + 1, gcd, aff->v->size - 1);
1145         isl_int_clear(gcd);
1146
1147         return aff;
1148 }
1149
1150 __isl_give isl_aff *isl_aff_scale_down(__isl_take isl_aff *aff, isl_int f)
1151 {
1152         isl_int gcd;
1153
1154         if (isl_int_is_one(f))
1155                 return aff;
1156
1157         aff = isl_aff_cow(aff);
1158         if (!aff)
1159                 return NULL;
1160
1161         if (isl_int_is_zero(f))
1162                 isl_die(isl_aff_get_ctx(aff), isl_error_invalid,
1163                         "cannot scale down by zero", return isl_aff_free(aff));
1164
1165         aff->v = isl_vec_cow(aff->v);
1166         if (!aff->v)
1167                 return isl_aff_free(aff);
1168
1169         isl_int_init(gcd);
1170         isl_seq_gcd(aff->v->el + 1, aff->v->size - 1, &gcd);
1171         isl_int_gcd(gcd, gcd, f);
1172         isl_seq_scale_down(aff->v->el + 1, aff->v->el + 1, gcd, aff->v->size - 1);
1173         isl_int_divexact(gcd, f, gcd);
1174         isl_int_mul(aff->v->el[0], aff->v->el[0], gcd);
1175         isl_int_clear(gcd);
1176
1177         return aff;
1178 }
1179
1180 __isl_give isl_aff *isl_aff_scale_down_ui(__isl_take isl_aff *aff, unsigned f)
1181 {
1182         isl_int v;
1183
1184         if (f == 1)
1185                 return aff;
1186
1187         isl_int_init(v);
1188         isl_int_set_ui(v, f);
1189         aff = isl_aff_scale_down(aff, v);
1190         isl_int_clear(v);
1191
1192         return aff;
1193 }
1194
1195 __isl_give isl_aff *isl_aff_set_dim_name(__isl_take isl_aff *aff,
1196         enum isl_dim_type type, unsigned pos, const char *s)
1197 {
1198         aff = isl_aff_cow(aff);
1199         if (!aff)
1200                 return NULL;
1201         if (type == isl_dim_out)
1202                 isl_die(aff->v->ctx, isl_error_invalid,
1203                         "cannot set name of output/set dimension",
1204                         return isl_aff_free(aff));
1205         if (type == isl_dim_in)
1206                 type = isl_dim_set;
1207         aff->ls = isl_local_space_set_dim_name(aff->ls, type, pos, s);
1208         if (!aff->ls)
1209                 return isl_aff_free(aff);
1210
1211         return aff;
1212 }
1213
1214 __isl_give isl_aff *isl_aff_set_dim_id(__isl_take isl_aff *aff,
1215         enum isl_dim_type type, unsigned pos, __isl_take isl_id *id)
1216 {
1217         aff = isl_aff_cow(aff);
1218         if (!aff)
1219                 return isl_id_free(id);
1220         if (type == isl_dim_out)
1221                 isl_die(aff->v->ctx, isl_error_invalid,
1222                         "cannot set name of output/set dimension",
1223                         goto error);
1224         if (type == isl_dim_in)
1225                 type = isl_dim_set;
1226         aff->ls = isl_local_space_set_dim_id(aff->ls, type, pos, id);
1227         if (!aff->ls)
1228                 return isl_aff_free(aff);
1229
1230         return aff;
1231 error:
1232         isl_id_free(id);
1233         isl_aff_free(aff);
1234         return NULL;
1235 }
1236
1237 /* Exploit the equalities in "eq" to simplify the affine expression
1238  * and the expressions of the integer divisions in the local space.
1239  * The integer divisions in this local space are assumed to appear
1240  * as regular dimensions in "eq".
1241  */
1242 static __isl_give isl_aff *isl_aff_substitute_equalities_lifted(
1243         __isl_take isl_aff *aff, __isl_take isl_basic_set *eq)
1244 {
1245         int i, j;
1246         unsigned total;
1247         unsigned n_div;
1248
1249         if (!eq)
1250                 goto error;
1251         if (eq->n_eq == 0) {
1252                 isl_basic_set_free(eq);
1253                 return aff;
1254         }
1255
1256         aff = isl_aff_cow(aff);
1257         if (!aff)
1258                 goto error;
1259
1260         aff->ls = isl_local_space_substitute_equalities(aff->ls,
1261                                                         isl_basic_set_copy(eq));
1262         if (!aff->ls)
1263                 goto error;
1264
1265         total = 1 + isl_space_dim(eq->dim, isl_dim_all);
1266         n_div = eq->n_div;
1267         for (i = 0; i < eq->n_eq; ++i) {
1268                 j = isl_seq_last_non_zero(eq->eq[i], total + n_div);
1269                 if (j < 0 || j == 0 || j >= total)
1270                         continue;
1271
1272                 isl_seq_elim(aff->v->el + 1, eq->eq[i], j, total,
1273                                 &aff->v->el[0]);
1274         }
1275
1276         isl_basic_set_free(eq);
1277         aff = isl_aff_normalize(aff);
1278         return aff;
1279 error:
1280         isl_basic_set_free(eq);
1281         isl_aff_free(aff);
1282         return NULL;
1283 }
1284
1285 /* Exploit the equalities in "eq" to simplify the affine expression
1286  * and the expressions of the integer divisions in the local space.
1287  */
1288 static __isl_give isl_aff *isl_aff_substitute_equalities(
1289         __isl_take isl_aff *aff, __isl_take isl_basic_set *eq)
1290 {
1291         int n_div;
1292
1293         if (!aff || !eq)
1294                 goto error;
1295         n_div = isl_local_space_dim(aff->ls, isl_dim_div);
1296         if (n_div > 0)
1297                 eq = isl_basic_set_add(eq, isl_dim_set, n_div);
1298         return isl_aff_substitute_equalities_lifted(aff, eq);
1299 error:
1300         isl_basic_set_free(eq);
1301         isl_aff_free(aff);
1302         return NULL;
1303 }
1304
1305 /* Look for equalities among the variables shared by context and aff
1306  * and the integer divisions of aff, if any.
1307  * The equalities are then used to eliminate coefficients and/or integer
1308  * divisions from aff.
1309  */
1310 __isl_give isl_aff *isl_aff_gist(__isl_take isl_aff *aff,
1311         __isl_take isl_set *context)
1312 {
1313         isl_basic_set *hull;
1314         int n_div;
1315
1316         if (!aff)
1317                 goto error;
1318         n_div = isl_local_space_dim(aff->ls, isl_dim_div);
1319         if (n_div > 0) {
1320                 isl_basic_set *bset;
1321                 isl_local_space *ls;
1322                 context = isl_set_add_dims(context, isl_dim_set, n_div);
1323                 ls = isl_aff_get_domain_local_space(aff);
1324                 bset = isl_basic_set_from_local_space(ls);
1325                 bset = isl_basic_set_lift(bset);
1326                 bset = isl_basic_set_flatten(bset);
1327                 context = isl_set_intersect(context,
1328                                             isl_set_from_basic_set(bset));
1329         }
1330
1331         hull = isl_set_affine_hull(context);
1332         return isl_aff_substitute_equalities_lifted(aff, hull);
1333 error:
1334         isl_aff_free(aff);
1335         isl_set_free(context);
1336         return NULL;
1337 }
1338
1339 __isl_give isl_aff *isl_aff_gist_params(__isl_take isl_aff *aff,
1340         __isl_take isl_set *context)
1341 {
1342         isl_set *dom_context = isl_set_universe(isl_aff_get_domain_space(aff));
1343         dom_context = isl_set_intersect_params(dom_context, context);
1344         return isl_aff_gist(aff, dom_context);
1345 }
1346
1347 /* Return a basic set containing those elements in the space
1348  * of aff where it is non-negative.
1349  */
1350 __isl_give isl_basic_set *isl_aff_nonneg_basic_set(__isl_take isl_aff *aff)
1351 {
1352         isl_constraint *ineq;
1353         isl_basic_set *bset;
1354
1355         ineq = isl_inequality_from_aff(aff);
1356
1357         bset = isl_basic_set_from_constraint(ineq);
1358         bset = isl_basic_set_simplify(bset);
1359         return bset;
1360 }
1361
1362 /* Return a basic set containing those elements in the domain space
1363  * of aff where it is negative.
1364  */
1365 __isl_give isl_basic_set *isl_aff_neg_basic_set(__isl_take isl_aff *aff)
1366 {
1367         aff = isl_aff_neg(aff);
1368         aff = isl_aff_add_constant_num_si(aff, -1);
1369         return isl_aff_nonneg_basic_set(aff);
1370 }
1371
1372 /* Return a basic set containing those elements in the space
1373  * of aff where it is zero.
1374  */
1375 __isl_give isl_basic_set *isl_aff_zero_basic_set(__isl_take isl_aff *aff)
1376 {
1377         isl_constraint *ineq;
1378         isl_basic_set *bset;
1379
1380         ineq = isl_equality_from_aff(aff);
1381
1382         bset = isl_basic_set_from_constraint(ineq);
1383         bset = isl_basic_set_simplify(bset);
1384         return bset;
1385 }
1386
1387 /* Return a basic set containing those elements in the shared space
1388  * of aff1 and aff2 where aff1 is greater than or equal to aff2.
1389  */
1390 __isl_give isl_basic_set *isl_aff_ge_basic_set(__isl_take isl_aff *aff1,
1391         __isl_take isl_aff *aff2)
1392 {
1393         aff1 = isl_aff_sub(aff1, aff2);
1394
1395         return isl_aff_nonneg_basic_set(aff1);
1396 }
1397
1398 /* Return a basic set containing those elements in the shared space
1399  * of aff1 and aff2 where aff1 is smaller than or equal to aff2.
1400  */
1401 __isl_give isl_basic_set *isl_aff_le_basic_set(__isl_take isl_aff *aff1,
1402         __isl_take isl_aff *aff2)
1403 {
1404         return isl_aff_ge_basic_set(aff2, aff1);
1405 }
1406
1407 __isl_give isl_aff *isl_aff_add_on_domain(__isl_keep isl_set *dom,
1408         __isl_take isl_aff *aff1, __isl_take isl_aff *aff2)
1409 {
1410         aff1 = isl_aff_add(aff1, aff2);
1411         aff1 = isl_aff_gist(aff1, isl_set_copy(dom));
1412         return aff1;
1413 }
1414
1415 int isl_aff_is_empty(__isl_keep isl_aff *aff)
1416 {
1417         if (!aff)
1418                 return -1;
1419
1420         return 0;
1421 }
1422
1423 /* Check whether the given affine expression has non-zero coefficient
1424  * for any dimension in the given range or if any of these dimensions
1425  * appear with non-zero coefficients in any of the integer divisions
1426  * involved in the affine expression.
1427  */
1428 int isl_aff_involves_dims(__isl_keep isl_aff *aff,
1429         enum isl_dim_type type, unsigned first, unsigned n)
1430 {
1431         int i;
1432         isl_ctx *ctx;
1433         int *active = NULL;
1434         int involves = 0;
1435
1436         if (!aff)
1437                 return -1;
1438         if (n == 0)
1439                 return 0;
1440
1441         ctx = isl_aff_get_ctx(aff);
1442         if (first + n > isl_aff_dim(aff, type))
1443                 isl_die(ctx, isl_error_invalid,
1444                         "range out of bounds", return -1);
1445
1446         active = isl_local_space_get_active(aff->ls, aff->v->el + 2);
1447         if (!active)
1448                 goto error;
1449
1450         first += isl_local_space_offset(aff->ls, type) - 1;
1451         for (i = 0; i < n; ++i)
1452                 if (active[first + i]) {
1453                         involves = 1;
1454                         break;
1455                 }
1456
1457         free(active);
1458
1459         return involves;
1460 error:
1461         free(active);
1462         return -1;
1463 }
1464
1465 __isl_give isl_aff *isl_aff_drop_dims(__isl_take isl_aff *aff,
1466         enum isl_dim_type type, unsigned first, unsigned n)
1467 {
1468         isl_ctx *ctx;
1469
1470         if (!aff)
1471                 return NULL;
1472         if (type == isl_dim_out)
1473                 isl_die(aff->v->ctx, isl_error_invalid,
1474                         "cannot drop output/set dimension",
1475                         return isl_aff_free(aff));
1476         if (type == isl_dim_in)
1477                 type = isl_dim_set;
1478         if (n == 0 && !isl_local_space_is_named_or_nested(aff->ls, type))
1479                 return aff;
1480
1481         ctx = isl_aff_get_ctx(aff);
1482         if (first + n > isl_local_space_dim(aff->ls, type))
1483                 isl_die(ctx, isl_error_invalid, "range out of bounds",
1484                         return isl_aff_free(aff));
1485
1486         aff = isl_aff_cow(aff);
1487         if (!aff)
1488                 return NULL;
1489
1490         aff->ls = isl_local_space_drop_dims(aff->ls, type, first, n);
1491         if (!aff->ls)
1492                 return isl_aff_free(aff);
1493
1494         first += 1 + isl_local_space_offset(aff->ls, type);
1495         aff->v = isl_vec_drop_els(aff->v, first, n);
1496         if (!aff->v)
1497                 return isl_aff_free(aff);
1498
1499         return aff;
1500 }
1501
1502 /* Project the domain of the affine expression onto its parameter space.
1503  * The affine expression may not involve any of the domain dimensions.
1504  */
1505 __isl_give isl_aff *isl_aff_project_domain_on_params(__isl_take isl_aff *aff)
1506 {
1507         isl_space *space;
1508         unsigned n;
1509         int involves;
1510
1511         n = isl_aff_dim(aff, isl_dim_in);
1512         involves = isl_aff_involves_dims(aff, isl_dim_in, 0, n);
1513         if (involves < 0)
1514                 return isl_aff_free(aff);
1515         if (involves)
1516                 isl_die(isl_aff_get_ctx(aff), isl_error_invalid,
1517                     "affine expression involves some of the domain dimensions",
1518                     return isl_aff_free(aff));
1519         aff = isl_aff_drop_dims(aff, isl_dim_in, 0, n);
1520         space = isl_aff_get_domain_space(aff);
1521         space = isl_space_params(space);
1522         aff = isl_aff_reset_domain_space(aff, space);
1523         return aff;
1524 }
1525
1526 __isl_give isl_aff *isl_aff_insert_dims(__isl_take isl_aff *aff,
1527         enum isl_dim_type type, unsigned first, unsigned n)
1528 {
1529         isl_ctx *ctx;
1530
1531         if (!aff)
1532                 return NULL;
1533         if (type == isl_dim_out)
1534                 isl_die(aff->v->ctx, isl_error_invalid,
1535                         "cannot insert output/set dimensions",
1536                         return isl_aff_free(aff));
1537         if (type == isl_dim_in)
1538                 type = isl_dim_set;
1539         if (n == 0 && !isl_local_space_is_named_or_nested(aff->ls, type))
1540                 return aff;
1541
1542         ctx = isl_aff_get_ctx(aff);
1543         if (first > isl_local_space_dim(aff->ls, type))
1544                 isl_die(ctx, isl_error_invalid, "position out of bounds",
1545                         return isl_aff_free(aff));
1546
1547         aff = isl_aff_cow(aff);
1548         if (!aff)
1549                 return NULL;
1550
1551         aff->ls = isl_local_space_insert_dims(aff->ls, type, first, n);
1552         if (!aff->ls)
1553                 return isl_aff_free(aff);
1554
1555         first += 1 + isl_local_space_offset(aff->ls, type);
1556         aff->v = isl_vec_insert_zero_els(aff->v, first, n);
1557         if (!aff->v)
1558                 return isl_aff_free(aff);
1559
1560         return aff;
1561 }
1562
1563 __isl_give isl_aff *isl_aff_add_dims(__isl_take isl_aff *aff,
1564         enum isl_dim_type type, unsigned n)
1565 {
1566         unsigned pos;
1567
1568         pos = isl_aff_dim(aff, type);
1569
1570         return isl_aff_insert_dims(aff, type, pos, n);
1571 }
1572
1573 __isl_give isl_pw_aff *isl_pw_aff_add_dims(__isl_take isl_pw_aff *pwaff,
1574         enum isl_dim_type type, unsigned n)
1575 {
1576         unsigned pos;
1577
1578         pos = isl_pw_aff_dim(pwaff, type);
1579
1580         return isl_pw_aff_insert_dims(pwaff, type, pos, n);
1581 }
1582
1583 __isl_give isl_pw_aff *isl_pw_aff_from_aff(__isl_take isl_aff *aff)
1584 {
1585         isl_set *dom = isl_set_universe(isl_aff_get_domain_space(aff));
1586         return isl_pw_aff_alloc(dom, aff);
1587 }
1588
1589 #undef PW
1590 #define PW isl_pw_aff
1591 #undef EL
1592 #define EL isl_aff
1593 #undef EL_IS_ZERO
1594 #define EL_IS_ZERO is_empty
1595 #undef ZERO
1596 #define ZERO empty
1597 #undef IS_ZERO
1598 #define IS_ZERO is_empty
1599 #undef FIELD
1600 #define FIELD aff
1601 #undef DEFAULT_IS_ZERO
1602 #define DEFAULT_IS_ZERO 0
1603
1604 #define NO_EVAL
1605 #define NO_OPT
1606 #define NO_MOVE_DIMS
1607 #define NO_LIFT
1608 #define NO_MORPH
1609
1610 #include <isl_pw_templ.c>
1611
1612 static __isl_give isl_set *align_params_pw_pw_set_and(
1613         __isl_take isl_pw_aff *pwaff1, __isl_take isl_pw_aff *pwaff2,
1614         __isl_give isl_set *(*fn)(__isl_take isl_pw_aff *pwaff1,
1615                                     __isl_take isl_pw_aff *pwaff2))
1616 {
1617         if (!pwaff1 || !pwaff2)
1618                 goto error;
1619         if (isl_space_match(pwaff1->dim, isl_dim_param,
1620                           pwaff2->dim, isl_dim_param))
1621                 return fn(pwaff1, pwaff2);
1622         if (!isl_space_has_named_params(pwaff1->dim) ||
1623             !isl_space_has_named_params(pwaff2->dim))
1624                 isl_die(isl_pw_aff_get_ctx(pwaff1), isl_error_invalid,
1625                         "unaligned unnamed parameters", goto error);
1626         pwaff1 = isl_pw_aff_align_params(pwaff1, isl_pw_aff_get_space(pwaff2));
1627         pwaff2 = isl_pw_aff_align_params(pwaff2, isl_pw_aff_get_space(pwaff1));
1628         return fn(pwaff1, pwaff2);
1629 error:
1630         isl_pw_aff_free(pwaff1);
1631         isl_pw_aff_free(pwaff2);
1632         return NULL;
1633 }
1634
1635 /* Compute a piecewise quasi-affine expression with a domain that
1636  * is the union of those of pwaff1 and pwaff2 and such that on each
1637  * cell, the quasi-affine expression is the better (according to cmp)
1638  * of those of pwaff1 and pwaff2.  If only one of pwaff1 or pwaff2
1639  * is defined on a given cell, then the associated expression
1640  * is the defined one.
1641  */
1642 static __isl_give isl_pw_aff *pw_aff_union_opt(__isl_take isl_pw_aff *pwaff1,
1643         __isl_take isl_pw_aff *pwaff2,
1644         __isl_give isl_basic_set *(*cmp)(__isl_take isl_aff *aff1,
1645                                         __isl_take isl_aff *aff2))
1646 {
1647         int i, j, n;
1648         isl_pw_aff *res;
1649         isl_ctx *ctx;
1650         isl_set *set;
1651
1652         if (!pwaff1 || !pwaff2)
1653                 goto error;
1654
1655         ctx = isl_space_get_ctx(pwaff1->dim);
1656         if (!isl_space_is_equal(pwaff1->dim, pwaff2->dim))
1657                 isl_die(ctx, isl_error_invalid,
1658                         "arguments should live in same space", goto error);
1659
1660         if (isl_pw_aff_is_empty(pwaff1)) {
1661                 isl_pw_aff_free(pwaff1);
1662                 return pwaff2;
1663         }
1664
1665         if (isl_pw_aff_is_empty(pwaff2)) {
1666                 isl_pw_aff_free(pwaff2);
1667                 return pwaff1;
1668         }
1669
1670         n = 2 * (pwaff1->n + 1) * (pwaff2->n + 1);
1671         res = isl_pw_aff_alloc_size(isl_space_copy(pwaff1->dim), n);
1672
1673         for (i = 0; i < pwaff1->n; ++i) {
1674                 set = isl_set_copy(pwaff1->p[i].set);
1675                 for (j = 0; j < pwaff2->n; ++j) {
1676                         struct isl_set *common;
1677                         isl_set *better;
1678
1679                         common = isl_set_intersect(
1680                                         isl_set_copy(pwaff1->p[i].set),
1681                                         isl_set_copy(pwaff2->p[j].set));
1682                         better = isl_set_from_basic_set(cmp(
1683                                         isl_aff_copy(pwaff2->p[j].aff),
1684                                         isl_aff_copy(pwaff1->p[i].aff)));
1685                         better = isl_set_intersect(common, better);
1686                         if (isl_set_plain_is_empty(better)) {
1687                                 isl_set_free(better);
1688                                 continue;
1689                         }
1690                         set = isl_set_subtract(set, isl_set_copy(better));
1691
1692                         res = isl_pw_aff_add_piece(res, better,
1693                                                 isl_aff_copy(pwaff2->p[j].aff));
1694                 }
1695                 res = isl_pw_aff_add_piece(res, set,
1696                                                 isl_aff_copy(pwaff1->p[i].aff));
1697         }
1698
1699         for (j = 0; j < pwaff2->n; ++j) {
1700                 set = isl_set_copy(pwaff2->p[j].set);
1701                 for (i = 0; i < pwaff1->n; ++i)
1702                         set = isl_set_subtract(set,
1703                                         isl_set_copy(pwaff1->p[i].set));
1704                 res = isl_pw_aff_add_piece(res, set,
1705                                                 isl_aff_copy(pwaff2->p[j].aff));
1706         }
1707
1708         isl_pw_aff_free(pwaff1);
1709         isl_pw_aff_free(pwaff2);
1710
1711         return res;
1712 error:
1713         isl_pw_aff_free(pwaff1);
1714         isl_pw_aff_free(pwaff2);
1715         return NULL;
1716 }
1717
1718 /* Compute a piecewise quasi-affine expression with a domain that
1719  * is the union of those of pwaff1 and pwaff2 and such that on each
1720  * cell, the quasi-affine expression is the maximum of those of pwaff1
1721  * and pwaff2.  If only one of pwaff1 or pwaff2 is defined on a given
1722  * cell, then the associated expression is the defined one.
1723  */
1724 static __isl_give isl_pw_aff *pw_aff_union_max(__isl_take isl_pw_aff *pwaff1,
1725         __isl_take isl_pw_aff *pwaff2)
1726 {
1727         return pw_aff_union_opt(pwaff1, pwaff2, &isl_aff_ge_basic_set);
1728 }
1729
1730 __isl_give isl_pw_aff *isl_pw_aff_union_max(__isl_take isl_pw_aff *pwaff1,
1731         __isl_take isl_pw_aff *pwaff2)
1732 {
1733         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2,
1734                                                         &pw_aff_union_max);
1735 }
1736
1737 /* Compute a piecewise quasi-affine expression with a domain that
1738  * is the union of those of pwaff1 and pwaff2 and such that on each
1739  * cell, the quasi-affine expression is the minimum of those of pwaff1
1740  * and pwaff2.  If only one of pwaff1 or pwaff2 is defined on a given
1741  * cell, then the associated expression is the defined one.
1742  */
1743 static __isl_give isl_pw_aff *pw_aff_union_min(__isl_take isl_pw_aff *pwaff1,
1744         __isl_take isl_pw_aff *pwaff2)
1745 {
1746         return pw_aff_union_opt(pwaff1, pwaff2, &isl_aff_le_basic_set);
1747 }
1748
1749 __isl_give isl_pw_aff *isl_pw_aff_union_min(__isl_take isl_pw_aff *pwaff1,
1750         __isl_take isl_pw_aff *pwaff2)
1751 {
1752         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2,
1753                                                         &pw_aff_union_min);
1754 }
1755
1756 __isl_give isl_pw_aff *isl_pw_aff_union_opt(__isl_take isl_pw_aff *pwaff1,
1757         __isl_take isl_pw_aff *pwaff2, int max)
1758 {
1759         if (max)
1760                 return isl_pw_aff_union_max(pwaff1, pwaff2);
1761         else
1762                 return isl_pw_aff_union_min(pwaff1, pwaff2);
1763 }
1764
1765 /* Construct a map with as domain the domain of pwaff and
1766  * one-dimensional range corresponding to the affine expressions.
1767  */
1768 static __isl_give isl_map *map_from_pw_aff(__isl_take isl_pw_aff *pwaff)
1769 {
1770         int i;
1771         isl_space *dim;
1772         isl_map *map;
1773
1774         if (!pwaff)
1775                 return NULL;
1776
1777         dim = isl_pw_aff_get_space(pwaff);
1778         map = isl_map_empty(dim);
1779
1780         for (i = 0; i < pwaff->n; ++i) {
1781                 isl_basic_map *bmap;
1782                 isl_map *map_i;
1783
1784                 bmap = isl_basic_map_from_aff(isl_aff_copy(pwaff->p[i].aff));
1785                 map_i = isl_map_from_basic_map(bmap);
1786                 map_i = isl_map_intersect_domain(map_i,
1787                                                 isl_set_copy(pwaff->p[i].set));
1788                 map = isl_map_union_disjoint(map, map_i);
1789         }
1790
1791         isl_pw_aff_free(pwaff);
1792
1793         return map;
1794 }
1795
1796 /* Construct a map with as domain the domain of pwaff and
1797  * one-dimensional range corresponding to the affine expressions.
1798  */
1799 __isl_give isl_map *isl_map_from_pw_aff(__isl_take isl_pw_aff *pwaff)
1800 {
1801         if (!pwaff)
1802                 return NULL;
1803         if (isl_space_is_set(pwaff->dim))
1804                 isl_die(isl_pw_aff_get_ctx(pwaff), isl_error_invalid,
1805                         "space of input is not a map",
1806                         return isl_pw_aff_free(pwaff));
1807         return map_from_pw_aff(pwaff);
1808 }
1809
1810 /* Construct a one-dimensional set with as parameter domain
1811  * the domain of pwaff and the single set dimension
1812  * corresponding to the affine expressions.
1813  */
1814 __isl_give isl_set *isl_set_from_pw_aff(__isl_take isl_pw_aff *pwaff)
1815 {
1816         if (!pwaff)
1817                 return NULL;
1818         if (!isl_space_is_set(pwaff->dim))
1819                 isl_die(isl_pw_aff_get_ctx(pwaff), isl_error_invalid,
1820                         "space of input is not a set",
1821                         return isl_pw_aff_free(pwaff));
1822         return map_from_pw_aff(pwaff);
1823 }
1824
1825 /* Return a set containing those elements in the domain
1826  * of pwaff where it is non-negative.
1827  */
1828 __isl_give isl_set *isl_pw_aff_nonneg_set(__isl_take isl_pw_aff *pwaff)
1829 {
1830         int i;
1831         isl_set *set;
1832
1833         if (!pwaff)
1834                 return NULL;
1835
1836         set = isl_set_empty(isl_pw_aff_get_domain_space(pwaff));
1837
1838         for (i = 0; i < pwaff->n; ++i) {
1839                 isl_basic_set *bset;
1840                 isl_set *set_i;
1841
1842                 bset = isl_aff_nonneg_basic_set(isl_aff_copy(pwaff->p[i].aff));
1843                 set_i = isl_set_from_basic_set(bset);
1844                 set_i = isl_set_intersect(set_i, isl_set_copy(pwaff->p[i].set));
1845                 set = isl_set_union_disjoint(set, set_i);
1846         }
1847
1848         isl_pw_aff_free(pwaff);
1849
1850         return set;
1851 }
1852
1853 /* Return a set containing those elements in the domain
1854  * of pwaff where it is zero (if complement is 0) or not zero
1855  * (if complement is 1).
1856  */
1857 static __isl_give isl_set *pw_aff_zero_set(__isl_take isl_pw_aff *pwaff,
1858         int complement)
1859 {
1860         int i;
1861         isl_set *set;
1862
1863         if (!pwaff)
1864                 return NULL;
1865
1866         set = isl_set_empty(isl_pw_aff_get_domain_space(pwaff));
1867
1868         for (i = 0; i < pwaff->n; ++i) {
1869                 isl_basic_set *bset;
1870                 isl_set *set_i, *zero;
1871
1872                 bset = isl_aff_zero_basic_set(isl_aff_copy(pwaff->p[i].aff));
1873                 zero = isl_set_from_basic_set(bset);
1874                 set_i = isl_set_copy(pwaff->p[i].set);
1875                 if (complement)
1876                         set_i = isl_set_subtract(set_i, zero);
1877                 else
1878                         set_i = isl_set_intersect(set_i, zero);
1879                 set = isl_set_union_disjoint(set, set_i);
1880         }
1881
1882         isl_pw_aff_free(pwaff);
1883
1884         return set;
1885 }
1886
1887 /* Return a set containing those elements in the domain
1888  * of pwaff where it is zero.
1889  */
1890 __isl_give isl_set *isl_pw_aff_zero_set(__isl_take isl_pw_aff *pwaff)
1891 {
1892         return pw_aff_zero_set(pwaff, 0);
1893 }
1894
1895 /* Return a set containing those elements in the domain
1896  * of pwaff where it is not zero.
1897  */
1898 __isl_give isl_set *isl_pw_aff_non_zero_set(__isl_take isl_pw_aff *pwaff)
1899 {
1900         return pw_aff_zero_set(pwaff, 1);
1901 }
1902
1903 /* Return a set containing those elements in the shared domain
1904  * of pwaff1 and pwaff2 where pwaff1 is greater than (or equal) to pwaff2.
1905  *
1906  * We compute the difference on the shared domain and then construct
1907  * the set of values where this difference is non-negative.
1908  * If strict is set, we first subtract 1 from the difference.
1909  * If equal is set, we only return the elements where pwaff1 and pwaff2
1910  * are equal.
1911  */
1912 static __isl_give isl_set *pw_aff_gte_set(__isl_take isl_pw_aff *pwaff1,
1913         __isl_take isl_pw_aff *pwaff2, int strict, int equal)
1914 {
1915         isl_set *set1, *set2;
1916
1917         set1 = isl_pw_aff_domain(isl_pw_aff_copy(pwaff1));
1918         set2 = isl_pw_aff_domain(isl_pw_aff_copy(pwaff2));
1919         set1 = isl_set_intersect(set1, set2);
1920         pwaff1 = isl_pw_aff_intersect_domain(pwaff1, isl_set_copy(set1));
1921         pwaff2 = isl_pw_aff_intersect_domain(pwaff2, isl_set_copy(set1));
1922         pwaff1 = isl_pw_aff_add(pwaff1, isl_pw_aff_neg(pwaff2));
1923
1924         if (strict) {
1925                 isl_space *dim = isl_set_get_space(set1);
1926                 isl_aff *aff;
1927                 aff = isl_aff_zero_on_domain(isl_local_space_from_space(dim));
1928                 aff = isl_aff_add_constant_si(aff, -1);
1929                 pwaff1 = isl_pw_aff_add(pwaff1, isl_pw_aff_alloc(set1, aff));
1930         } else
1931                 isl_set_free(set1);
1932
1933         if (equal)
1934                 return isl_pw_aff_zero_set(pwaff1);
1935         return isl_pw_aff_nonneg_set(pwaff1);
1936 }
1937
1938 /* Return a set containing those elements in the shared domain
1939  * of pwaff1 and pwaff2 where pwaff1 is equal to pwaff2.
1940  */
1941 static __isl_give isl_set *pw_aff_eq_set(__isl_take isl_pw_aff *pwaff1,
1942         __isl_take isl_pw_aff *pwaff2)
1943 {
1944         return pw_aff_gte_set(pwaff1, pwaff2, 0, 1);
1945 }
1946
1947 __isl_give isl_set *isl_pw_aff_eq_set(__isl_take isl_pw_aff *pwaff1,
1948         __isl_take isl_pw_aff *pwaff2)
1949 {
1950         return align_params_pw_pw_set_and(pwaff1, pwaff2, &pw_aff_eq_set);
1951 }
1952
1953 /* Return a set containing those elements in the shared domain
1954  * of pwaff1 and pwaff2 where pwaff1 is greater than or equal to pwaff2.
1955  */
1956 static __isl_give isl_set *pw_aff_ge_set(__isl_take isl_pw_aff *pwaff1,
1957         __isl_take isl_pw_aff *pwaff2)
1958 {
1959         return pw_aff_gte_set(pwaff1, pwaff2, 0, 0);
1960 }
1961
1962 __isl_give isl_set *isl_pw_aff_ge_set(__isl_take isl_pw_aff *pwaff1,
1963         __isl_take isl_pw_aff *pwaff2)
1964 {
1965         return align_params_pw_pw_set_and(pwaff1, pwaff2, &pw_aff_ge_set);
1966 }
1967
1968 /* Return a set containing those elements in the shared domain
1969  * of pwaff1 and pwaff2 where pwaff1 is strictly greater than pwaff2.
1970  */
1971 static __isl_give isl_set *pw_aff_gt_set(__isl_take isl_pw_aff *pwaff1,
1972         __isl_take isl_pw_aff *pwaff2)
1973 {
1974         return pw_aff_gte_set(pwaff1, pwaff2, 1, 0);
1975 }
1976
1977 __isl_give isl_set *isl_pw_aff_gt_set(__isl_take isl_pw_aff *pwaff1,
1978         __isl_take isl_pw_aff *pwaff2)
1979 {
1980         return align_params_pw_pw_set_and(pwaff1, pwaff2, &pw_aff_gt_set);
1981 }
1982
1983 __isl_give isl_set *isl_pw_aff_le_set(__isl_take isl_pw_aff *pwaff1,
1984         __isl_take isl_pw_aff *pwaff2)
1985 {
1986         return isl_pw_aff_ge_set(pwaff2, pwaff1);
1987 }
1988
1989 __isl_give isl_set *isl_pw_aff_lt_set(__isl_take isl_pw_aff *pwaff1,
1990         __isl_take isl_pw_aff *pwaff2)
1991 {
1992         return isl_pw_aff_gt_set(pwaff2, pwaff1);
1993 }
1994
1995 /* Return a set containing those elements in the shared domain
1996  * of the elements of list1 and list2 where each element in list1
1997  * has the relation specified by "fn" with each element in list2.
1998  */
1999 static __isl_give isl_set *pw_aff_list_set(__isl_take isl_pw_aff_list *list1,
2000         __isl_take isl_pw_aff_list *list2,
2001         __isl_give isl_set *(*fn)(__isl_take isl_pw_aff *pwaff1,
2002                                     __isl_take isl_pw_aff *pwaff2))
2003 {
2004         int i, j;
2005         isl_ctx *ctx;
2006         isl_set *set;
2007
2008         if (!list1 || !list2)
2009                 goto error;
2010
2011         ctx = isl_pw_aff_list_get_ctx(list1);
2012         if (list1->n < 1 || list2->n < 1)
2013                 isl_die(ctx, isl_error_invalid,
2014                         "list should contain at least one element", goto error);
2015
2016         set = isl_set_universe(isl_pw_aff_get_domain_space(list1->p[0]));
2017         for (i = 0; i < list1->n; ++i)
2018                 for (j = 0; j < list2->n; ++j) {
2019                         isl_set *set_ij;
2020
2021                         set_ij = fn(isl_pw_aff_copy(list1->p[i]),
2022                                     isl_pw_aff_copy(list2->p[j]));
2023                         set = isl_set_intersect(set, set_ij);
2024                 }
2025
2026         isl_pw_aff_list_free(list1);
2027         isl_pw_aff_list_free(list2);
2028         return set;
2029 error:
2030         isl_pw_aff_list_free(list1);
2031         isl_pw_aff_list_free(list2);
2032         return NULL;
2033 }
2034
2035 /* Return a set containing those elements in the shared domain
2036  * of the elements of list1 and list2 where each element in list1
2037  * is equal to each element in list2.
2038  */
2039 __isl_give isl_set *isl_pw_aff_list_eq_set(__isl_take isl_pw_aff_list *list1,
2040         __isl_take isl_pw_aff_list *list2)
2041 {
2042         return pw_aff_list_set(list1, list2, &isl_pw_aff_eq_set);
2043 }
2044
2045 __isl_give isl_set *isl_pw_aff_list_ne_set(__isl_take isl_pw_aff_list *list1,
2046         __isl_take isl_pw_aff_list *list2)
2047 {
2048         return pw_aff_list_set(list1, list2, &isl_pw_aff_ne_set);
2049 }
2050
2051 /* Return a set containing those elements in the shared domain
2052  * of the elements of list1 and list2 where each element in list1
2053  * is less than or equal to each element in list2.
2054  */
2055 __isl_give isl_set *isl_pw_aff_list_le_set(__isl_take isl_pw_aff_list *list1,
2056         __isl_take isl_pw_aff_list *list2)
2057 {
2058         return pw_aff_list_set(list1, list2, &isl_pw_aff_le_set);
2059 }
2060
2061 __isl_give isl_set *isl_pw_aff_list_lt_set(__isl_take isl_pw_aff_list *list1,
2062         __isl_take isl_pw_aff_list *list2)
2063 {
2064         return pw_aff_list_set(list1, list2, &isl_pw_aff_lt_set);
2065 }
2066
2067 __isl_give isl_set *isl_pw_aff_list_ge_set(__isl_take isl_pw_aff_list *list1,
2068         __isl_take isl_pw_aff_list *list2)
2069 {
2070         return pw_aff_list_set(list1, list2, &isl_pw_aff_ge_set);
2071 }
2072
2073 __isl_give isl_set *isl_pw_aff_list_gt_set(__isl_take isl_pw_aff_list *list1,
2074         __isl_take isl_pw_aff_list *list2)
2075 {
2076         return pw_aff_list_set(list1, list2, &isl_pw_aff_gt_set);
2077 }
2078
2079
2080 /* Return a set containing those elements in the shared domain
2081  * of pwaff1 and pwaff2 where pwaff1 is not equal to pwaff2.
2082  */
2083 static __isl_give isl_set *pw_aff_ne_set(__isl_take isl_pw_aff *pwaff1,
2084         __isl_take isl_pw_aff *pwaff2)
2085 {
2086         isl_set *set_lt, *set_gt;
2087
2088         set_lt = isl_pw_aff_lt_set(isl_pw_aff_copy(pwaff1),
2089                                    isl_pw_aff_copy(pwaff2));
2090         set_gt = isl_pw_aff_gt_set(pwaff1, pwaff2);
2091         return isl_set_union_disjoint(set_lt, set_gt);
2092 }
2093
2094 __isl_give isl_set *isl_pw_aff_ne_set(__isl_take isl_pw_aff *pwaff1,
2095         __isl_take isl_pw_aff *pwaff2)
2096 {
2097         return align_params_pw_pw_set_and(pwaff1, pwaff2, &pw_aff_ne_set);
2098 }
2099
2100 __isl_give isl_pw_aff *isl_pw_aff_scale_down(__isl_take isl_pw_aff *pwaff,
2101         isl_int v)
2102 {
2103         int i;
2104
2105         if (isl_int_is_one(v))
2106                 return pwaff;
2107         if (!isl_int_is_pos(v))
2108                 isl_die(isl_pw_aff_get_ctx(pwaff), isl_error_invalid,
2109                         "factor needs to be positive",
2110                         return isl_pw_aff_free(pwaff));
2111         pwaff = isl_pw_aff_cow(pwaff);
2112         if (!pwaff)
2113                 return NULL;
2114         if (pwaff->n == 0)
2115                 return pwaff;
2116
2117         for (i = 0; i < pwaff->n; ++i) {
2118                 pwaff->p[i].aff = isl_aff_scale_down(pwaff->p[i].aff, v);
2119                 if (!pwaff->p[i].aff)
2120                         return isl_pw_aff_free(pwaff);
2121         }
2122
2123         return pwaff;
2124 }
2125
2126 __isl_give isl_pw_aff *isl_pw_aff_floor(__isl_take isl_pw_aff *pwaff)
2127 {
2128         int i;
2129
2130         pwaff = isl_pw_aff_cow(pwaff);
2131         if (!pwaff)
2132                 return NULL;
2133         if (pwaff->n == 0)
2134                 return pwaff;
2135
2136         for (i = 0; i < pwaff->n; ++i) {
2137                 pwaff->p[i].aff = isl_aff_floor(pwaff->p[i].aff);
2138                 if (!pwaff->p[i].aff)
2139                         return isl_pw_aff_free(pwaff);
2140         }
2141
2142         return pwaff;
2143 }
2144
2145 __isl_give isl_pw_aff *isl_pw_aff_ceil(__isl_take isl_pw_aff *pwaff)
2146 {
2147         int i;
2148
2149         pwaff = isl_pw_aff_cow(pwaff);
2150         if (!pwaff)
2151                 return NULL;
2152         if (pwaff->n == 0)
2153                 return pwaff;
2154
2155         for (i = 0; i < pwaff->n; ++i) {
2156                 pwaff->p[i].aff = isl_aff_ceil(pwaff->p[i].aff);
2157                 if (!pwaff->p[i].aff)
2158                         return isl_pw_aff_free(pwaff);
2159         }
2160
2161         return pwaff;
2162 }
2163
2164 /* Assuming that "cond1" and "cond2" are disjoint,
2165  * return an affine expression that is equal to pwaff1 on cond1
2166  * and to pwaff2 on cond2.
2167  */
2168 static __isl_give isl_pw_aff *isl_pw_aff_select(
2169         __isl_take isl_set *cond1, __isl_take isl_pw_aff *pwaff1,
2170         __isl_take isl_set *cond2, __isl_take isl_pw_aff *pwaff2)
2171 {
2172         pwaff1 = isl_pw_aff_intersect_domain(pwaff1, cond1);
2173         pwaff2 = isl_pw_aff_intersect_domain(pwaff2, cond2);
2174
2175         return isl_pw_aff_add_disjoint(pwaff1, pwaff2);
2176 }
2177
2178 /* Return an affine expression that is equal to pwaff_true for elements
2179  * where "cond" is non-zero and to pwaff_false for elements where "cond"
2180  * is zero.
2181  * That is, return cond ? pwaff_true : pwaff_false;
2182  */
2183 __isl_give isl_pw_aff *isl_pw_aff_cond(__isl_take isl_pw_aff *cond,
2184         __isl_take isl_pw_aff *pwaff_true, __isl_take isl_pw_aff *pwaff_false)
2185 {
2186         isl_set *cond_true, *cond_false;
2187
2188         cond_true = isl_pw_aff_non_zero_set(isl_pw_aff_copy(cond));
2189         cond_false = isl_pw_aff_zero_set(cond);
2190         return isl_pw_aff_select(cond_true, pwaff_true,
2191                                  cond_false, pwaff_false);
2192 }
2193
2194 int isl_aff_is_cst(__isl_keep isl_aff *aff)
2195 {
2196         if (!aff)
2197                 return -1;
2198
2199         return isl_seq_first_non_zero(aff->v->el + 2, aff->v->size - 2) == -1;
2200 }
2201
2202 /* Check whether pwaff is a piecewise constant.
2203  */
2204 int isl_pw_aff_is_cst(__isl_keep isl_pw_aff *pwaff)
2205 {
2206         int i;
2207
2208         if (!pwaff)
2209                 return -1;
2210
2211         for (i = 0; i < pwaff->n; ++i) {
2212                 int is_cst = isl_aff_is_cst(pwaff->p[i].aff);
2213                 if (is_cst < 0 || !is_cst)
2214                         return is_cst;
2215         }
2216
2217         return 1;
2218 }
2219
2220 __isl_give isl_aff *isl_aff_mul(__isl_take isl_aff *aff1,
2221         __isl_take isl_aff *aff2)
2222 {
2223         if (!isl_aff_is_cst(aff2) && isl_aff_is_cst(aff1))
2224                 return isl_aff_mul(aff2, aff1);
2225
2226         if (!isl_aff_is_cst(aff2))
2227                 isl_die(isl_aff_get_ctx(aff1), isl_error_invalid,
2228                         "at least one affine expression should be constant",
2229                         goto error);
2230
2231         aff1 = isl_aff_cow(aff1);
2232         if (!aff1 || !aff2)
2233                 goto error;
2234
2235         aff1 = isl_aff_scale(aff1, aff2->v->el[1]);
2236         aff1 = isl_aff_scale_down(aff1, aff2->v->el[0]);
2237
2238         isl_aff_free(aff2);
2239         return aff1;
2240 error:
2241         isl_aff_free(aff1);
2242         isl_aff_free(aff2);
2243         return NULL;
2244 }
2245
2246 /* Divide "aff1" by "aff2", assuming "aff2" is a piecewise constant.
2247  */
2248 __isl_give isl_aff *isl_aff_div(__isl_take isl_aff *aff1,
2249         __isl_take isl_aff *aff2)
2250 {
2251         int is_cst;
2252         int neg;
2253
2254         is_cst = isl_aff_is_cst(aff2);
2255         if (is_cst < 0)
2256                 goto error;
2257         if (!is_cst)
2258                 isl_die(isl_aff_get_ctx(aff2), isl_error_invalid,
2259                         "second argument should be a constant", goto error);
2260
2261         if (!aff2)
2262                 goto error;
2263
2264         neg = isl_int_is_neg(aff2->v->el[1]);
2265         if (neg) {
2266                 isl_int_neg(aff2->v->el[0], aff2->v->el[0]);
2267                 isl_int_neg(aff2->v->el[1], aff2->v->el[1]);
2268         }
2269
2270         aff1 = isl_aff_scale(aff1, aff2->v->el[0]);
2271         aff1 = isl_aff_scale_down(aff1, aff2->v->el[1]);
2272
2273         if (neg) {
2274                 isl_int_neg(aff2->v->el[0], aff2->v->el[0]);
2275                 isl_int_neg(aff2->v->el[1], aff2->v->el[1]);
2276         }
2277
2278         isl_aff_free(aff2);
2279         return aff1;
2280 error:
2281         isl_aff_free(aff1);
2282         isl_aff_free(aff2);
2283         return NULL;
2284 }
2285
2286 static __isl_give isl_pw_aff *pw_aff_add(__isl_take isl_pw_aff *pwaff1,
2287         __isl_take isl_pw_aff *pwaff2)
2288 {
2289         return isl_pw_aff_on_shared_domain(pwaff1, pwaff2, &isl_aff_add);
2290 }
2291
2292 __isl_give isl_pw_aff *isl_pw_aff_add(__isl_take isl_pw_aff *pwaff1,
2293         __isl_take isl_pw_aff *pwaff2)
2294 {
2295         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2, &pw_aff_add);
2296 }
2297
2298 __isl_give isl_pw_aff *isl_pw_aff_union_add(__isl_take isl_pw_aff *pwaff1,
2299         __isl_take isl_pw_aff *pwaff2)
2300 {
2301         return isl_pw_aff_union_add_(pwaff1, pwaff2);
2302 }
2303
2304 static __isl_give isl_pw_aff *pw_aff_mul(__isl_take isl_pw_aff *pwaff1,
2305         __isl_take isl_pw_aff *pwaff2)
2306 {
2307         return isl_pw_aff_on_shared_domain(pwaff1, pwaff2, &isl_aff_mul);
2308 }
2309
2310 __isl_give isl_pw_aff *isl_pw_aff_mul(__isl_take isl_pw_aff *pwaff1,
2311         __isl_take isl_pw_aff *pwaff2)
2312 {
2313         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2, &pw_aff_mul);
2314 }
2315
2316 static __isl_give isl_pw_aff *pw_aff_div(__isl_take isl_pw_aff *pa1,
2317         __isl_take isl_pw_aff *pa2)
2318 {
2319         return isl_pw_aff_on_shared_domain(pa1, pa2, &isl_aff_div);
2320 }
2321
2322 /* Divide "pa1" by "pa2", assuming "pa2" is a piecewise constant.
2323  */
2324 __isl_give isl_pw_aff *isl_pw_aff_div(__isl_take isl_pw_aff *pa1,
2325         __isl_take isl_pw_aff *pa2)
2326 {
2327         int is_cst;
2328
2329         is_cst = isl_pw_aff_is_cst(pa2);
2330         if (is_cst < 0)
2331                 goto error;
2332         if (!is_cst)
2333                 isl_die(isl_pw_aff_get_ctx(pa2), isl_error_invalid,
2334                         "second argument should be a piecewise constant",
2335                         goto error);
2336         return isl_pw_aff_align_params_pw_pw_and(pa1, pa2, &pw_aff_div);
2337 error:
2338         isl_pw_aff_free(pa1);
2339         isl_pw_aff_free(pa2);
2340         return NULL;
2341 }
2342
2343 /* Compute the quotient of the integer division of "pa1" by "pa2"
2344  * with rounding towards zero.
2345  * "pa2" is assumed to be a piecewise constant.
2346  *
2347  * In particular, return
2348  *
2349  *      pa1 >= 0 ? floor(pa1/pa2) : ceil(pa1/pa2)
2350  *
2351  */
2352 __isl_give isl_pw_aff *isl_pw_aff_tdiv_q(__isl_take isl_pw_aff *pa1,
2353         __isl_take isl_pw_aff *pa2)
2354 {
2355         int is_cst;
2356         isl_set *cond;
2357         isl_pw_aff *f, *c;
2358
2359         is_cst = isl_pw_aff_is_cst(pa2);
2360         if (is_cst < 0)
2361                 goto error;
2362         if (!is_cst)
2363                 isl_die(isl_pw_aff_get_ctx(pa2), isl_error_invalid,
2364                         "second argument should be a piecewise constant",
2365                         goto error);
2366
2367         pa1 = isl_pw_aff_div(pa1, pa2);
2368
2369         cond = isl_pw_aff_nonneg_set(isl_pw_aff_copy(pa1));
2370         f = isl_pw_aff_floor(isl_pw_aff_copy(pa1));
2371         c = isl_pw_aff_ceil(pa1);
2372         return isl_pw_aff_cond(isl_set_indicator_function(cond), f, c);
2373 error:
2374         isl_pw_aff_free(pa1);
2375         isl_pw_aff_free(pa2);
2376         return NULL;
2377 }
2378
2379 /* Compute the remainder of the integer division of "pa1" by "pa2"
2380  * with rounding towards zero.
2381  * "pa2" is assumed to be a piecewise constant.
2382  *
2383  * In particular, return
2384  *
2385  *      pa1 - pa2 * (pa1 >= 0 ? floor(pa1/pa2) : ceil(pa1/pa2))
2386  *
2387  */
2388 __isl_give isl_pw_aff *isl_pw_aff_tdiv_r(__isl_take isl_pw_aff *pa1,
2389         __isl_take isl_pw_aff *pa2)
2390 {
2391         int is_cst;
2392         isl_pw_aff *res;
2393
2394         is_cst = isl_pw_aff_is_cst(pa2);
2395         if (is_cst < 0)
2396                 goto error;
2397         if (!is_cst)
2398                 isl_die(isl_pw_aff_get_ctx(pa2), isl_error_invalid,
2399                         "second argument should be a piecewise constant",
2400                         goto error);
2401         res = isl_pw_aff_tdiv_q(isl_pw_aff_copy(pa1), isl_pw_aff_copy(pa2));
2402         res = isl_pw_aff_mul(pa2, res);
2403         res = isl_pw_aff_sub(pa1, res);
2404         return res;
2405 error:
2406         isl_pw_aff_free(pa1);
2407         isl_pw_aff_free(pa2);
2408         return NULL;
2409 }
2410
2411 static __isl_give isl_pw_aff *pw_aff_min(__isl_take isl_pw_aff *pwaff1,
2412         __isl_take isl_pw_aff *pwaff2)
2413 {
2414         isl_set *le;
2415         isl_set *dom;
2416
2417         dom = isl_set_intersect(isl_pw_aff_domain(isl_pw_aff_copy(pwaff1)),
2418                                 isl_pw_aff_domain(isl_pw_aff_copy(pwaff2)));
2419         le = isl_pw_aff_le_set(isl_pw_aff_copy(pwaff1),
2420                                 isl_pw_aff_copy(pwaff2));
2421         dom = isl_set_subtract(dom, isl_set_copy(le));
2422         return isl_pw_aff_select(le, pwaff1, dom, pwaff2);
2423 }
2424
2425 __isl_give isl_pw_aff *isl_pw_aff_min(__isl_take isl_pw_aff *pwaff1,
2426         __isl_take isl_pw_aff *pwaff2)
2427 {
2428         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2, &pw_aff_min);
2429 }
2430
2431 static __isl_give isl_pw_aff *pw_aff_max(__isl_take isl_pw_aff *pwaff1,
2432         __isl_take isl_pw_aff *pwaff2)
2433 {
2434         isl_set *ge;
2435         isl_set *dom;
2436
2437         dom = isl_set_intersect(isl_pw_aff_domain(isl_pw_aff_copy(pwaff1)),
2438                                 isl_pw_aff_domain(isl_pw_aff_copy(pwaff2)));
2439         ge = isl_pw_aff_ge_set(isl_pw_aff_copy(pwaff1),
2440                                 isl_pw_aff_copy(pwaff2));
2441         dom = isl_set_subtract(dom, isl_set_copy(ge));
2442         return isl_pw_aff_select(ge, pwaff1, dom, pwaff2);
2443 }
2444
2445 __isl_give isl_pw_aff *isl_pw_aff_max(__isl_take isl_pw_aff *pwaff1,
2446         __isl_take isl_pw_aff *pwaff2)
2447 {
2448         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2, &pw_aff_max);
2449 }
2450
2451 static __isl_give isl_pw_aff *pw_aff_list_reduce(
2452         __isl_take isl_pw_aff_list *list,
2453         __isl_give isl_pw_aff *(*fn)(__isl_take isl_pw_aff *pwaff1,
2454                                         __isl_take isl_pw_aff *pwaff2))
2455 {
2456         int i;
2457         isl_ctx *ctx;
2458         isl_pw_aff *res;
2459
2460         if (!list)
2461                 return NULL;
2462
2463         ctx = isl_pw_aff_list_get_ctx(list);
2464         if (list->n < 1)
2465                 isl_die(ctx, isl_error_invalid,
2466                         "list should contain at least one element",
2467                         return isl_pw_aff_list_free(list));
2468
2469         res = isl_pw_aff_copy(list->p[0]);
2470         for (i = 1; i < list->n; ++i)
2471                 res = fn(res, isl_pw_aff_copy(list->p[i]));
2472
2473         isl_pw_aff_list_free(list);
2474         return res;
2475 }
2476
2477 /* Return an isl_pw_aff that maps each element in the intersection of the
2478  * domains of the elements of list to the minimal corresponding affine
2479  * expression.
2480  */
2481 __isl_give isl_pw_aff *isl_pw_aff_list_min(__isl_take isl_pw_aff_list *list)
2482 {
2483         return pw_aff_list_reduce(list, &isl_pw_aff_min);
2484 }
2485
2486 /* Return an isl_pw_aff that maps each element in the intersection of the
2487  * domains of the elements of list to the maximal corresponding affine
2488  * expression.
2489  */
2490 __isl_give isl_pw_aff *isl_pw_aff_list_max(__isl_take isl_pw_aff_list *list)
2491 {
2492         return pw_aff_list_reduce(list, &isl_pw_aff_max);
2493 }
2494
2495 #undef BASE
2496 #define BASE aff
2497
2498 #include <isl_multi_templ.c>
2499
2500 /* Construct an isl_multi_aff in the given space with value zero in
2501  * each of the output dimensions.
2502  */
2503 __isl_give isl_multi_aff *isl_multi_aff_zero(__isl_take isl_space *space)
2504 {
2505         int n;
2506         isl_multi_aff *ma;
2507
2508         if (!space)
2509                 return NULL;
2510
2511         n = isl_space_dim(space , isl_dim_out);
2512         ma = isl_multi_aff_alloc(isl_space_copy(space));
2513
2514         if (!n)
2515                 isl_space_free(space);
2516         else {
2517                 int i;
2518                 isl_local_space *ls;
2519                 isl_aff *aff;
2520
2521                 space = isl_space_domain(space);
2522                 ls = isl_local_space_from_space(space);
2523                 aff = isl_aff_zero_on_domain(ls);
2524
2525                 for (i = 0; i < n; ++i)
2526                         ma = isl_multi_aff_set_aff(ma, i, isl_aff_copy(aff));
2527
2528                 isl_aff_free(aff);
2529         }
2530
2531         return ma;
2532 }
2533
2534 /* Create an isl_multi_aff in the given space that maps each
2535  * input dimension to the corresponding output dimension.
2536  */
2537 __isl_give isl_multi_aff *isl_multi_aff_identity(__isl_take isl_space *space)
2538 {
2539         int n;
2540         isl_multi_aff *ma;
2541
2542         if (!space)
2543                 return NULL;
2544
2545         if (isl_space_is_set(space))
2546                 isl_die(isl_space_get_ctx(space), isl_error_invalid,
2547                         "expecting map space", goto error);
2548
2549         n = isl_space_dim(space, isl_dim_out);
2550         if (n != isl_space_dim(space, isl_dim_in))
2551                 isl_die(isl_space_get_ctx(space), isl_error_invalid,
2552                         "number of input and output dimensions needs to be "
2553                         "the same", goto error);
2554
2555         ma = isl_multi_aff_alloc(isl_space_copy(space));
2556
2557         if (!n)
2558                 isl_space_free(space);
2559         else {
2560                 int i;
2561                 isl_local_space *ls;
2562                 isl_aff *aff;
2563
2564                 space = isl_space_domain(space);
2565                 ls = isl_local_space_from_space(space);
2566                 aff = isl_aff_zero_on_domain(ls);
2567
2568                 for (i = 0; i < n; ++i) {
2569                         isl_aff *aff_i;
2570                         aff_i = isl_aff_copy(aff);
2571                         aff_i = isl_aff_add_coefficient_si(aff_i,
2572                                                             isl_dim_in, i, 1);
2573                         ma = isl_multi_aff_set_aff(ma, i, aff_i);
2574                 }
2575
2576                 isl_aff_free(aff);
2577         }
2578
2579         return ma;
2580 error:
2581         isl_space_free(space);
2582         return NULL;
2583 }
2584
2585 /* Create an isl_pw_multi_aff with the given isl_multi_aff on a universe
2586  * domain.
2587  */
2588 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_from_multi_aff(
2589         __isl_take isl_multi_aff *ma)
2590 {
2591         isl_set *dom = isl_set_universe(isl_multi_aff_get_domain_space(ma));
2592         return isl_pw_multi_aff_alloc(dom, ma);
2593 }
2594
2595 /* Create a piecewise multi-affine expression in the given space that maps each
2596  * input dimension to the corresponding output dimension.
2597  */
2598 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_identity(
2599         __isl_take isl_space *space)
2600 {
2601         return isl_pw_multi_aff_from_multi_aff(isl_multi_aff_identity(space));
2602 }
2603
2604 __isl_give isl_multi_aff *isl_multi_aff_add(__isl_take isl_multi_aff *maff1,
2605         __isl_take isl_multi_aff *maff2)
2606 {
2607         int i;
2608         isl_ctx *ctx;
2609
2610         maff1 = isl_multi_aff_cow(maff1);
2611         if (!maff1 || !maff2)
2612                 goto error;
2613
2614         ctx = isl_multi_aff_get_ctx(maff1);
2615         if (!isl_space_is_equal(maff1->space, maff2->space))
2616                 isl_die(ctx, isl_error_invalid,
2617                         "spaces don't match", goto error);
2618
2619         for (i = 0; i < maff1->n; ++i) {
2620                 maff1->p[i] = isl_aff_add(maff1->p[i],
2621                                             isl_aff_copy(maff2->p[i]));
2622                 if (!maff1->p[i])
2623                         goto error;
2624         }
2625
2626         isl_multi_aff_free(maff2);
2627         return maff1;
2628 error:
2629         isl_multi_aff_free(maff1);
2630         isl_multi_aff_free(maff2);
2631         return NULL;
2632 }
2633
2634 /* Given two multi-affine expressions A -> B and C -> D,
2635  * construct a multi-affine expression [A -> C] -> [B -> D].
2636  */
2637 __isl_give isl_multi_aff *isl_multi_aff_product(
2638         __isl_take isl_multi_aff *ma1, __isl_take isl_multi_aff *ma2)
2639 {
2640         int i;
2641         isl_aff *aff;
2642         isl_space *space;
2643         isl_multi_aff *res;
2644         int in1, in2, out1, out2;
2645
2646         in1 = isl_multi_aff_dim(ma1, isl_dim_in);
2647         in2 = isl_multi_aff_dim(ma2, isl_dim_in);
2648         out1 = isl_multi_aff_dim(ma1, isl_dim_out);
2649         out2 = isl_multi_aff_dim(ma2, isl_dim_out);
2650         space = isl_space_product(isl_multi_aff_get_space(ma1),
2651                                   isl_multi_aff_get_space(ma2));
2652         res = isl_multi_aff_alloc(isl_space_copy(space));
2653         space = isl_space_domain(space);
2654
2655         for (i = 0; i < out1; ++i) {
2656                 aff = isl_multi_aff_get_aff(ma1, i);
2657                 aff = isl_aff_insert_dims(aff, isl_dim_in, in1, in2);
2658                 aff = isl_aff_reset_domain_space(aff, isl_space_copy(space));
2659                 res = isl_multi_aff_set_aff(res, i, aff);
2660         }
2661
2662         for (i = 0; i < out2; ++i) {
2663                 aff = isl_multi_aff_get_aff(ma2, i);
2664                 aff = isl_aff_insert_dims(aff, isl_dim_in, 0, in1);
2665                 aff = isl_aff_reset_domain_space(aff, isl_space_copy(space));
2666                 res = isl_multi_aff_set_aff(res, out1 + i, aff);
2667         }
2668
2669         isl_space_free(space);
2670         isl_multi_aff_free(ma1);
2671         isl_multi_aff_free(ma2);
2672         return res;
2673 }
2674
2675 /* Exploit the equalities in "eq" to simplify the affine expressions.
2676  */
2677 static __isl_give isl_multi_aff *isl_multi_aff_substitute_equalities(
2678         __isl_take isl_multi_aff *maff, __isl_take isl_basic_set *eq)
2679 {
2680         int i;
2681
2682         maff = isl_multi_aff_cow(maff);
2683         if (!maff || !eq)
2684                 goto error;
2685
2686         for (i = 0; i < maff->n; ++i) {
2687                 maff->p[i] = isl_aff_substitute_equalities(maff->p[i],
2688                                                     isl_basic_set_copy(eq));
2689                 if (!maff->p[i])
2690                         goto error;
2691         }
2692
2693         isl_basic_set_free(eq);
2694         return maff;
2695 error:
2696         isl_basic_set_free(eq);
2697         isl_multi_aff_free(maff);
2698         return NULL;
2699 }
2700
2701 __isl_give isl_multi_aff *isl_multi_aff_scale(__isl_take isl_multi_aff *maff,
2702         isl_int f)
2703 {
2704         int i;
2705
2706         maff = isl_multi_aff_cow(maff);
2707         if (!maff)
2708                 return NULL;
2709
2710         for (i = 0; i < maff->n; ++i) {
2711                 maff->p[i] = isl_aff_scale(maff->p[i], f);
2712                 if (!maff->p[i])
2713                         return isl_multi_aff_free(maff);
2714         }
2715
2716         return maff;
2717 }
2718
2719 __isl_give isl_multi_aff *isl_multi_aff_add_on_domain(__isl_keep isl_set *dom,
2720         __isl_take isl_multi_aff *maff1, __isl_take isl_multi_aff *maff2)
2721 {
2722         maff1 = isl_multi_aff_add(maff1, maff2);
2723         maff1 = isl_multi_aff_gist(maff1, isl_set_copy(dom));
2724         return maff1;
2725 }
2726
2727 int isl_multi_aff_is_empty(__isl_keep isl_multi_aff *maff)
2728 {
2729         if (!maff)
2730                 return -1;
2731
2732         return 0;
2733 }
2734
2735 int isl_multi_aff_plain_is_equal(__isl_keep isl_multi_aff *maff1,
2736         __isl_keep isl_multi_aff *maff2)
2737 {
2738         int i;
2739         int equal;
2740
2741         if (!maff1 || !maff2)
2742                 return -1;
2743         if (maff1->n != maff2->n)
2744                 return 0;
2745         equal = isl_space_is_equal(maff1->space, maff2->space);
2746         if (equal < 0 || !equal)
2747                 return equal;
2748
2749         for (i = 0; i < maff1->n; ++i) {
2750                 equal = isl_aff_plain_is_equal(maff1->p[i], maff2->p[i]);
2751                 if (equal < 0 || !equal)
2752                         return equal;
2753         }
2754
2755         return 1;
2756 }
2757
2758 __isl_give isl_multi_aff *isl_multi_aff_set_dim_name(
2759         __isl_take isl_multi_aff *maff,
2760         enum isl_dim_type type, unsigned pos, const char *s)
2761 {
2762         int i;
2763
2764         maff = isl_multi_aff_cow(maff);
2765         if (!maff)
2766                 return NULL;
2767
2768         maff->space = isl_space_set_dim_name(maff->space, type, pos, s);
2769         if (!maff->space)
2770                 return isl_multi_aff_free(maff);
2771
2772         if (type == isl_dim_out)
2773                 return maff;
2774         for (i = 0; i < maff->n; ++i) {
2775                 maff->p[i] = isl_aff_set_dim_name(maff->p[i], type, pos, s);
2776                 if (!maff->p[i])
2777                         return isl_multi_aff_free(maff);
2778         }
2779
2780         return maff;
2781 }
2782
2783 __isl_give isl_multi_aff *isl_multi_aff_drop_dims(__isl_take isl_multi_aff *maff,
2784         enum isl_dim_type type, unsigned first, unsigned n)
2785 {
2786         int i;
2787
2788         maff = isl_multi_aff_cow(maff);
2789         if (!maff)
2790                 return NULL;
2791
2792         maff->space = isl_space_drop_dims(maff->space, type, first, n);
2793         if (!maff->space)
2794                 return isl_multi_aff_free(maff);
2795
2796         if (type == isl_dim_out) {
2797                 for (i = 0; i < n; ++i)
2798                         isl_aff_free(maff->p[first + i]);
2799                 for (i = first; i + n < maff->n; ++i)
2800                         maff->p[i] = maff->p[i + n];
2801                 maff->n -= n;
2802                 return maff;
2803         }
2804
2805         for (i = 0; i < maff->n; ++i) {
2806                 maff->p[i] = isl_aff_drop_dims(maff->p[i], type, first, n);
2807                 if (!maff->p[i])
2808                         return isl_multi_aff_free(maff);
2809         }
2810
2811         return maff;
2812 }
2813
2814 /* Return the set of domain elements where "ma1" is lexicographically
2815  * smaller than or equal to "ma2".
2816  */
2817 __isl_give isl_set *isl_multi_aff_lex_le_set(__isl_take isl_multi_aff *ma1,
2818         __isl_take isl_multi_aff *ma2)
2819 {
2820         return isl_multi_aff_lex_ge_set(ma2, ma1);
2821 }
2822
2823 /* Return the set of domain elements where "ma1" is lexicographically
2824  * greater than or equal to "ma2".
2825  */
2826 __isl_give isl_set *isl_multi_aff_lex_ge_set(__isl_take isl_multi_aff *ma1,
2827         __isl_take isl_multi_aff *ma2)
2828 {
2829         isl_space *space;
2830         isl_map *map1, *map2;
2831         isl_map *map, *ge;
2832
2833         map1 = isl_map_from_multi_aff(ma1);
2834         map2 = isl_map_from_multi_aff(ma2);
2835         map = isl_map_range_product(map1, map2);
2836         space = isl_space_range(isl_map_get_space(map));
2837         space = isl_space_domain(isl_space_unwrap(space));
2838         ge = isl_map_lex_ge(space);
2839         map = isl_map_intersect_range(map, isl_map_wrap(ge));
2840
2841         return isl_map_domain(map);
2842 }
2843
2844 #undef PW
2845 #define PW isl_pw_multi_aff
2846 #undef EL
2847 #define EL isl_multi_aff
2848 #undef EL_IS_ZERO
2849 #define EL_IS_ZERO is_empty
2850 #undef ZERO
2851 #define ZERO empty
2852 #undef IS_ZERO
2853 #define IS_ZERO is_empty
2854 #undef FIELD
2855 #define FIELD maff
2856 #undef DEFAULT_IS_ZERO
2857 #define DEFAULT_IS_ZERO 0
2858
2859 #define NO_NEG
2860 #define NO_EVAL
2861 #define NO_OPT
2862 #define NO_INVOLVES_DIMS
2863 #define NO_MOVE_DIMS
2864 #define NO_INSERT_DIMS
2865 #define NO_LIFT
2866 #define NO_MORPH
2867
2868 #include <isl_pw_templ.c>
2869
2870 #undef UNION
2871 #define UNION isl_union_pw_multi_aff
2872 #undef PART
2873 #define PART isl_pw_multi_aff
2874 #undef PARTS
2875 #define PARTS pw_multi_aff
2876 #define ALIGN_DOMAIN
2877
2878 #define NO_EVAL
2879
2880 #include <isl_union_templ.c>
2881
2882 /* Given a function "cmp" that returns the set of elements where
2883  * "ma1" is "better" than "ma2", return the intersection of this
2884  * set with "dom1" and "dom2".
2885  */
2886 static __isl_give isl_set *shared_and_better(__isl_keep isl_set *dom1,
2887         __isl_keep isl_set *dom2, __isl_keep isl_multi_aff *ma1,
2888         __isl_keep isl_multi_aff *ma2,
2889         __isl_give isl_set *(*cmp)(__isl_take isl_multi_aff *ma1,
2890                                     __isl_take isl_multi_aff *ma2))
2891 {
2892         isl_set *common;
2893         isl_set *better;
2894         int is_empty;
2895
2896         common = isl_set_intersect(isl_set_copy(dom1), isl_set_copy(dom2));
2897         is_empty = isl_set_plain_is_empty(common);
2898         if (is_empty >= 0 && is_empty)
2899                 return common;
2900         if (is_empty < 0)
2901                 return isl_set_free(common);
2902         better = cmp(isl_multi_aff_copy(ma1), isl_multi_aff_copy(ma2));
2903         better = isl_set_intersect(common, better);
2904
2905         return better;
2906 }
2907
2908 /* Given a function "cmp" that returns the set of elements where
2909  * "ma1" is "better" than "ma2", return a piecewise multi affine
2910  * expression defined on the union of the definition domains
2911  * of "pma1" and "pma2" that maps to the "best" of "pma1" and
2912  * "pma2" on each cell.  If only one of the two input functions
2913  * is defined on a given cell, then it is considered the best.
2914  */
2915 static __isl_give isl_pw_multi_aff *pw_multi_aff_union_opt(
2916         __isl_take isl_pw_multi_aff *pma1,
2917         __isl_take isl_pw_multi_aff *pma2,
2918         __isl_give isl_set *(*cmp)(__isl_take isl_multi_aff *ma1,
2919                                     __isl_take isl_multi_aff *ma2))
2920 {
2921         int i, j, n;
2922         isl_pw_multi_aff *res = NULL;
2923         isl_ctx *ctx;
2924         isl_set *set = NULL;
2925
2926         if (!pma1 || !pma2)
2927                 goto error;
2928
2929         ctx = isl_space_get_ctx(pma1->dim);
2930         if (!isl_space_is_equal(pma1->dim, pma2->dim))
2931                 isl_die(ctx, isl_error_invalid,
2932                         "arguments should live in the same space", goto error);
2933
2934         if (isl_pw_multi_aff_is_empty(pma1)) {
2935                 isl_pw_multi_aff_free(pma1);
2936                 return pma2;
2937         }
2938
2939         if (isl_pw_multi_aff_is_empty(pma2)) {
2940                 isl_pw_multi_aff_free(pma2);
2941                 return pma1;
2942         }
2943
2944         n = 2 * (pma1->n + 1) * (pma2->n + 1);
2945         res = isl_pw_multi_aff_alloc_size(isl_space_copy(pma1->dim), n);
2946
2947         for (i = 0; i < pma1->n; ++i) {
2948                 set = isl_set_copy(pma1->p[i].set);
2949                 for (j = 0; j < pma2->n; ++j) {
2950                         isl_set *better;
2951                         int is_empty;
2952
2953                         better = shared_and_better(pma2->p[j].set,
2954                                         pma1->p[i].set, pma2->p[j].maff,
2955                                         pma1->p[i].maff, cmp);
2956                         is_empty = isl_set_plain_is_empty(better);
2957                         if (is_empty < 0 || is_empty) {
2958                                 isl_set_free(better);
2959                                 if (is_empty < 0)
2960                                         goto error;
2961                                 continue;
2962                         }
2963                         set = isl_set_subtract(set, isl_set_copy(better));
2964
2965                         res = isl_pw_multi_aff_add_piece(res, better,
2966                                         isl_multi_aff_copy(pma2->p[j].maff));
2967                 }
2968                 res = isl_pw_multi_aff_add_piece(res, set,
2969                                         isl_multi_aff_copy(pma1->p[i].maff));
2970         }
2971
2972         for (j = 0; j < pma2->n; ++j) {
2973                 set = isl_set_copy(pma2->p[j].set);
2974                 for (i = 0; i < pma1->n; ++i)
2975                         set = isl_set_subtract(set,
2976                                         isl_set_copy(pma1->p[i].set));
2977                 res = isl_pw_multi_aff_add_piece(res, set,
2978                                         isl_multi_aff_copy(pma2->p[j].maff));
2979         }
2980
2981         isl_pw_multi_aff_free(pma1);
2982         isl_pw_multi_aff_free(pma2);
2983
2984         return res;
2985 error:
2986         isl_pw_multi_aff_free(pma1);
2987         isl_pw_multi_aff_free(pma2);
2988         isl_set_free(set);
2989         return isl_pw_multi_aff_free(res);
2990 }
2991
2992 static __isl_give isl_pw_multi_aff *pw_multi_aff_union_lexmax(
2993         __isl_take isl_pw_multi_aff *pma1,
2994         __isl_take isl_pw_multi_aff *pma2)
2995 {
2996         return pw_multi_aff_union_opt(pma1, pma2, &isl_multi_aff_lex_ge_set);
2997 }
2998
2999 /* Given two piecewise multi affine expressions, return a piecewise
3000  * multi-affine expression defined on the union of the definition domains
3001  * of the inputs that is equal to the lexicographic maximum of the two
3002  * inputs on each cell.  If only one of the two inputs is defined on
3003  * a given cell, then it is considered to be the maximum.
3004  */
3005 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_union_lexmax(
3006         __isl_take isl_pw_multi_aff *pma1,
3007         __isl_take isl_pw_multi_aff *pma2)
3008 {
3009         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
3010                                                     &pw_multi_aff_union_lexmax);
3011 }
3012
3013 static __isl_give isl_pw_multi_aff *pw_multi_aff_union_lexmin(
3014         __isl_take isl_pw_multi_aff *pma1,
3015         __isl_take isl_pw_multi_aff *pma2)
3016 {
3017         return pw_multi_aff_union_opt(pma1, pma2, &isl_multi_aff_lex_le_set);
3018 }
3019
3020 /* Given two piecewise multi affine expressions, return a piecewise
3021  * multi-affine expression defined on the union of the definition domains
3022  * of the inputs that is equal to the lexicographic minimum of the two
3023  * inputs on each cell.  If only one of the two inputs is defined on
3024  * a given cell, then it is considered to be the minimum.
3025  */
3026 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_union_lexmin(
3027         __isl_take isl_pw_multi_aff *pma1,
3028         __isl_take isl_pw_multi_aff *pma2)
3029 {
3030         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
3031                                                     &pw_multi_aff_union_lexmin);
3032 }
3033
3034 static __isl_give isl_pw_multi_aff *pw_multi_aff_add(
3035         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
3036 {
3037         return isl_pw_multi_aff_on_shared_domain(pma1, pma2,
3038                                                 &isl_multi_aff_add);
3039 }
3040
3041 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_add(
3042         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
3043 {
3044         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
3045                                                 &pw_multi_aff_add);
3046 }
3047
3048 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_union_add(
3049         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
3050 {
3051         return isl_pw_multi_aff_union_add_(pma1, pma2);
3052 }
3053
3054 /* Given two piecewise multi-affine expressions A -> B and C -> D,
3055  * construct a piecewise multi-affine expression [A -> C] -> [B -> D].
3056  */
3057 static __isl_give isl_pw_multi_aff *pw_multi_aff_product(
3058         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
3059 {
3060         int i, j, n;
3061         isl_space *space;
3062         isl_pw_multi_aff *res;
3063
3064         if (!pma1 || !pma2)
3065                 goto error;
3066
3067         n = pma1->n * pma2->n;
3068         space = isl_space_product(isl_space_copy(pma1->dim),
3069                                   isl_space_copy(pma2->dim));
3070         res = isl_pw_multi_aff_alloc_size(space, n);
3071
3072         for (i = 0; i < pma1->n; ++i) {
3073                 for (j = 0; j < pma2->n; ++j) {
3074                         isl_set *domain;
3075                         isl_multi_aff *ma;
3076
3077                         domain = isl_set_product(isl_set_copy(pma1->p[i].set),
3078                                                  isl_set_copy(pma2->p[j].set));
3079                         ma = isl_multi_aff_product(
3080                                         isl_multi_aff_copy(pma1->p[i].maff),
3081                                         isl_multi_aff_copy(pma2->p[i].maff));
3082                         res = isl_pw_multi_aff_add_piece(res, domain, ma);
3083                 }
3084         }
3085
3086         isl_pw_multi_aff_free(pma1);
3087         isl_pw_multi_aff_free(pma2);
3088         return res;
3089 error:
3090         isl_pw_multi_aff_free(pma1);
3091         isl_pw_multi_aff_free(pma2);
3092         return NULL;
3093 }
3094
3095 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_product(
3096         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
3097 {
3098         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
3099                                                 &pw_multi_aff_product);
3100 }
3101
3102 /* Construct a map mapping the domain of the piecewise multi-affine expression
3103  * to its range, with each dimension in the range equated to the
3104  * corresponding affine expression on its cell.
3105  */
3106 __isl_give isl_map *isl_map_from_pw_multi_aff(__isl_take isl_pw_multi_aff *pma)
3107 {
3108         int i;
3109         isl_map *map;
3110
3111         if (!pma)
3112                 return NULL;
3113
3114         map = isl_map_empty(isl_pw_multi_aff_get_space(pma));
3115
3116         for (i = 0; i < pma->n; ++i) {
3117                 isl_multi_aff *maff;
3118                 isl_basic_map *bmap;
3119                 isl_map *map_i;
3120
3121                 maff = isl_multi_aff_copy(pma->p[i].maff);
3122                 bmap = isl_basic_map_from_multi_aff(maff);
3123                 map_i = isl_map_from_basic_map(bmap);
3124                 map_i = isl_map_intersect_domain(map_i,
3125                                                 isl_set_copy(pma->p[i].set));
3126                 map = isl_map_union_disjoint(map, map_i);
3127         }
3128
3129         isl_pw_multi_aff_free(pma);
3130         return map;
3131 }
3132
3133 __isl_give isl_set *isl_set_from_pw_multi_aff(__isl_take isl_pw_multi_aff *pma)
3134 {
3135         if (!pma)
3136                 return NULL;
3137
3138         if (!isl_space_is_set(pma->dim))
3139                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
3140                         "isl_pw_multi_aff cannot be converted into an isl_set",
3141                         return isl_pw_multi_aff_free(pma));
3142
3143         return isl_map_from_pw_multi_aff(pma);
3144 }
3145
3146 /* Given a basic map with a single output dimension that is defined
3147  * in terms of the parameters and input dimensions using an equality,
3148  * extract an isl_aff that expresses the output dimension in terms
3149  * of the parameters and input dimensions.
3150  *
3151  * Since some applications expect the result of isl_pw_multi_aff_from_map
3152  * to only contain integer affine expressions, we compute the floor
3153  * of the expression before returning.
3154  *
3155  * This function shares some similarities with
3156  * isl_basic_map_has_defining_equality and isl_constraint_get_bound.
3157  */
3158 static __isl_give isl_aff *extract_isl_aff_from_basic_map(
3159         __isl_take isl_basic_map *bmap)
3160 {
3161         int i;
3162         unsigned offset;
3163         unsigned total;
3164         isl_local_space *ls;
3165         isl_aff *aff;
3166
3167         if (!bmap)
3168                 return NULL;
3169         if (isl_basic_map_dim(bmap, isl_dim_out) != 1)
3170                 isl_die(isl_basic_map_get_ctx(bmap), isl_error_invalid,
3171                         "basic map should have a single output dimension",
3172                         goto error);
3173         offset = isl_basic_map_offset(bmap, isl_dim_out);
3174         total = isl_basic_map_total_dim(bmap);
3175         for (i = 0; i < bmap->n_eq; ++i) {
3176                 if (isl_int_is_zero(bmap->eq[i][offset]))
3177                         continue;
3178                 if (isl_seq_first_non_zero(bmap->eq[i] + offset + 1,
3179                                            1 + total - (offset + 1)) != -1)
3180                         continue;
3181                 break;
3182         }
3183         if (i >= bmap->n_eq)
3184                 isl_die(isl_basic_map_get_ctx(bmap), isl_error_invalid,
3185                         "unable to find suitable equality", goto error);
3186         ls = isl_basic_map_get_local_space(bmap);
3187         aff = isl_aff_alloc(isl_local_space_domain(ls));
3188         if (!aff)
3189                 goto error;
3190         if (isl_int_is_neg(bmap->eq[i][offset]))
3191                 isl_seq_cpy(aff->v->el + 1, bmap->eq[i], offset);
3192         else
3193                 isl_seq_neg(aff->v->el + 1, bmap->eq[i], offset);
3194         isl_seq_clr(aff->v->el + 1 + offset, aff->v->size - (1 + offset));
3195         isl_int_abs(aff->v->el[0], bmap->eq[i][offset]);
3196         isl_basic_map_free(bmap);
3197
3198         aff = isl_aff_remove_unused_divs(aff);
3199         aff = isl_aff_floor(aff);
3200         return aff;
3201 error:
3202         isl_basic_map_free(bmap);
3203         return NULL;
3204 }
3205
3206 /* Given a basic map where each output dimension is defined
3207  * in terms of the parameters and input dimensions using an equality,
3208  * extract an isl_multi_aff that expresses the output dimensions in terms
3209  * of the parameters and input dimensions.
3210  */
3211 static __isl_give isl_multi_aff *extract_isl_multi_aff_from_basic_map(
3212         __isl_take isl_basic_map *bmap)
3213 {
3214         int i;
3215         unsigned n_out;
3216         isl_multi_aff *ma;
3217
3218         if (!bmap)
3219                 return NULL;
3220
3221         ma = isl_multi_aff_alloc(isl_basic_map_get_space(bmap));
3222         n_out = isl_basic_map_dim(bmap, isl_dim_out);
3223
3224         for (i = 0; i < n_out; ++i) {
3225                 isl_basic_map *bmap_i;
3226                 isl_aff *aff;
3227
3228                 bmap_i = isl_basic_map_copy(bmap);
3229                 bmap_i = isl_basic_map_project_out(bmap_i, isl_dim_out,
3230                                                         i + 1, n_out - (1 + i));
3231                 bmap_i = isl_basic_map_project_out(bmap_i, isl_dim_out, 0, i);
3232                 aff = extract_isl_aff_from_basic_map(bmap_i);
3233                 ma = isl_multi_aff_set_aff(ma, i, aff);
3234         }
3235
3236         isl_basic_map_free(bmap);
3237
3238         return ma;
3239 }
3240
3241 /* Create an isl_pw_multi_aff that is equivalent to
3242  * isl_map_intersect_domain(isl_map_from_basic_map(bmap), domain).
3243  * The given basic map is such that each output dimension is defined
3244  * in terms of the parameters and input dimensions using an equality.
3245  */
3246 static __isl_give isl_pw_multi_aff *plain_pw_multi_aff_from_map(
3247         __isl_take isl_set *domain, __isl_take isl_basic_map *bmap)
3248 {
3249         isl_multi_aff *ma;
3250
3251         ma = extract_isl_multi_aff_from_basic_map(bmap);
3252         return isl_pw_multi_aff_alloc(domain, ma);
3253 }
3254
3255 /* Try and create an isl_pw_multi_aff that is equivalent to the given isl_map.
3256  * This obviously only works if the input "map" is single-valued.
3257  * If so, we compute the lexicographic minimum of the image in the form
3258  * of an isl_pw_multi_aff.  Since the image is unique, it is equal
3259  * to its lexicographic minimum.
3260  * If the input is not single-valued, we produce an error.
3261  *
3262  * As a special case, we first check if all output dimensions are uniquely
3263  * defined in terms of the parameters and input dimensions over the entire
3264  * domain.  If so, we extract the desired isl_pw_multi_aff directly
3265  * from the affine hull of "map" and its domain.
3266  */
3267 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_from_map(__isl_take isl_map *map)
3268 {
3269         int i;
3270         int sv;
3271         isl_pw_multi_aff *pma;
3272         isl_basic_map *hull;
3273
3274         if (!map)
3275                 return NULL;
3276
3277         hull = isl_map_affine_hull(isl_map_copy(map));
3278         sv = isl_basic_map_plain_is_single_valued(hull);
3279         if (sv >= 0 && sv)
3280                 return plain_pw_multi_aff_from_map(isl_map_domain(map), hull);
3281         isl_basic_map_free(hull);
3282         if (sv < 0)
3283                 goto error;
3284
3285         sv = isl_map_is_single_valued(map);
3286         if (sv < 0)
3287                 goto error;
3288         if (!sv)
3289                 isl_die(isl_map_get_ctx(map), isl_error_invalid,
3290                         "map is not single-valued", goto error);
3291         map = isl_map_make_disjoint(map);
3292         if (!map)
3293                 return NULL;
3294
3295         pma = isl_pw_multi_aff_empty(isl_map_get_space(map));
3296
3297         for (i = 0; i < map->n; ++i) {
3298                 isl_pw_multi_aff *pma_i;
3299                 isl_basic_map *bmap;
3300                 bmap = isl_basic_map_copy(map->p[i]);
3301                 pma_i = isl_basic_map_lexmin_pw_multi_aff(bmap);
3302                 pma = isl_pw_multi_aff_add_disjoint(pma, pma_i);
3303         }
3304
3305         isl_map_free(map);
3306         return pma;
3307 error:
3308         isl_map_free(map);
3309         return NULL;
3310 }
3311
3312 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_from_set(__isl_take isl_set *set)
3313 {
3314         return isl_pw_multi_aff_from_map(set);
3315 }
3316
3317 /* Return the piecewise affine expression "set ? 1 : 0".
3318  */
3319 __isl_give isl_pw_aff *isl_set_indicator_function(__isl_take isl_set *set)
3320 {
3321         isl_pw_aff *pa;
3322         isl_space *space = isl_set_get_space(set);
3323         isl_local_space *ls = isl_local_space_from_space(space);
3324         isl_aff *zero = isl_aff_zero_on_domain(isl_local_space_copy(ls));
3325         isl_aff *one = isl_aff_zero_on_domain(ls);
3326
3327         one = isl_aff_add_constant_si(one, 1);
3328         pa = isl_pw_aff_alloc(isl_set_copy(set), one);
3329         set = isl_set_complement(set);
3330         pa = isl_pw_aff_add_disjoint(pa, isl_pw_aff_alloc(set, zero));
3331
3332         return pa;
3333 }
3334
3335 /* Plug in "subs" for dimension "type", "pos" of "aff".
3336  *
3337  * Let i be the dimension to replace and let "subs" be of the form
3338  *
3339  *      f/d
3340  *
3341  * and "aff" of the form
3342  *
3343  *      (a i + g)/m
3344  *
3345  * The result is
3346  *
3347  *      (a f + d g')/(m d)
3348  *
3349  * where g' is the result of plugging in "subs" in each of the integer
3350  * divisions in g.
3351  */
3352 __isl_give isl_aff *isl_aff_substitute(__isl_take isl_aff *aff,
3353         enum isl_dim_type type, unsigned pos, __isl_keep isl_aff *subs)
3354 {
3355         isl_ctx *ctx;
3356         isl_int v;
3357
3358         aff = isl_aff_cow(aff);
3359         if (!aff || !subs)
3360                 return isl_aff_free(aff);
3361
3362         ctx = isl_aff_get_ctx(aff);
3363         if (!isl_space_is_equal(aff->ls->dim, subs->ls->dim))
3364                 isl_die(ctx, isl_error_invalid,
3365                         "spaces don't match", return isl_aff_free(aff));
3366         if (isl_local_space_dim(subs->ls, isl_dim_div) != 0)
3367                 isl_die(ctx, isl_error_unsupported,
3368                         "cannot handle divs yet", return isl_aff_free(aff));
3369
3370         aff->ls = isl_local_space_substitute(aff->ls, type, pos, subs);
3371         if (!aff->ls)
3372                 return isl_aff_free(aff);
3373
3374         aff->v = isl_vec_cow(aff->v);
3375         if (!aff->v)
3376                 return isl_aff_free(aff);
3377
3378         pos += isl_local_space_offset(aff->ls, type);
3379
3380         isl_int_init(v);
3381         isl_seq_substitute(aff->v->el, pos, subs->v->el,
3382                             aff->v->size, subs->v->size, v);
3383         isl_int_clear(v);
3384
3385         return aff;
3386 }
3387
3388 /* Plug in "subs" for dimension "type", "pos" in each of the affine
3389  * expressions in "maff".
3390  */
3391 __isl_give isl_multi_aff *isl_multi_aff_substitute(
3392         __isl_take isl_multi_aff *maff, enum isl_dim_type type, unsigned pos,
3393         __isl_keep isl_aff *subs)
3394 {
3395         int i;
3396
3397         maff = isl_multi_aff_cow(maff);
3398         if (!maff || !subs)
3399                 return isl_multi_aff_free(maff);
3400
3401         if (type == isl_dim_in)
3402                 type = isl_dim_set;
3403
3404         for (i = 0; i < maff->n; ++i) {
3405                 maff->p[i] = isl_aff_substitute(maff->p[i], type, pos, subs);
3406                 if (!maff->p[i])
3407                         return isl_multi_aff_free(maff);
3408         }
3409
3410         return maff;
3411 }
3412
3413 /* Plug in "subs" for dimension "type", "pos" of "pma".
3414  *
3415  * pma is of the form
3416  *
3417  *      A_i(v) -> M_i(v)
3418  *
3419  * while subs is of the form
3420  *
3421  *      v' = B_j(v) -> S_j
3422  *
3423  * Each pair i,j such that C_ij = A_i \cap B_i is non-empty
3424  * has a contribution in the result, in particular
3425  *
3426  *      C_ij(S_j) -> M_i(S_j)
3427  *
3428  * Note that plugging in S_j in C_ij may also result in an empty set
3429  * and this contribution should simply be discarded.
3430  */
3431 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_substitute(
3432         __isl_take isl_pw_multi_aff *pma, enum isl_dim_type type, unsigned pos,
3433         __isl_keep isl_pw_aff *subs)
3434 {
3435         int i, j, n;
3436         isl_pw_multi_aff *res;
3437
3438         if (!pma || !subs)
3439                 return isl_pw_multi_aff_free(pma);
3440
3441         n = pma->n * subs->n;
3442         res = isl_pw_multi_aff_alloc_size(isl_space_copy(pma->dim), n);
3443
3444         for (i = 0; i < pma->n; ++i) {
3445                 for (j = 0; j < subs->n; ++j) {
3446                         isl_set *common;
3447                         isl_multi_aff *res_ij;
3448                         common = isl_set_intersect(
3449                                         isl_set_copy(pma->p[i].set),
3450                                         isl_set_copy(subs->p[j].set));
3451                         common = isl_set_substitute(common,
3452                                         type, pos, subs->p[j].aff);
3453                         if (isl_set_plain_is_empty(common)) {
3454                                 isl_set_free(common);
3455                                 continue;
3456                         }
3457
3458                         res_ij = isl_multi_aff_substitute(
3459                                         isl_multi_aff_copy(pma->p[i].maff),
3460                                         type, pos, subs->p[j].aff);
3461
3462                         res = isl_pw_multi_aff_add_piece(res, common, res_ij);
3463                 }
3464         }
3465
3466         isl_pw_multi_aff_free(pma);
3467         return res;
3468 }
3469
3470 /* Extend the local space of "dst" to include the divs
3471  * in the local space of "src".
3472  */
3473 __isl_give isl_aff *isl_aff_align_divs(__isl_take isl_aff *dst,
3474         __isl_keep isl_aff *src)
3475 {
3476         isl_ctx *ctx;
3477         int *exp1 = NULL;
3478         int *exp2 = NULL;
3479         isl_mat *div;
3480
3481         if (!src || !dst)
3482                 return isl_aff_free(dst);
3483
3484         ctx = isl_aff_get_ctx(src);
3485         if (!isl_space_is_equal(src->ls->dim, dst->ls->dim))
3486                 isl_die(ctx, isl_error_invalid,
3487                         "spaces don't match", goto error);
3488
3489         if (src->ls->div->n_row == 0)
3490                 return dst;
3491
3492         exp1 = isl_alloc_array(ctx, int, src->ls->div->n_row);
3493         exp2 = isl_alloc_array(ctx, int, dst->ls->div->n_row);
3494         if (!exp1 || !exp2)
3495                 goto error;
3496
3497         div = isl_merge_divs(src->ls->div, dst->ls->div, exp1, exp2);
3498         dst = isl_aff_expand_divs(dst, div, exp2);
3499         free(exp1);
3500         free(exp2);
3501
3502         return dst;
3503 error:
3504         free(exp1);
3505         free(exp2);
3506         return isl_aff_free(dst);
3507 }
3508
3509 /* Adjust the local spaces of the affine expressions in "maff"
3510  * such that they all have the save divs.
3511  */
3512 __isl_give isl_multi_aff *isl_multi_aff_align_divs(
3513         __isl_take isl_multi_aff *maff)
3514 {
3515         int i;
3516
3517         if (!maff)
3518                 return NULL;
3519         if (maff->n == 0)
3520                 return maff;
3521         maff = isl_multi_aff_cow(maff);
3522         if (!maff)
3523                 return NULL;
3524
3525         for (i = 1; i < maff->n; ++i)
3526                 maff->p[0] = isl_aff_align_divs(maff->p[0], maff->p[i]);
3527         for (i = 1; i < maff->n; ++i) {
3528                 maff->p[i] = isl_aff_align_divs(maff->p[i], maff->p[0]);
3529                 if (!maff->p[i])
3530                         return isl_multi_aff_free(maff);
3531         }
3532
3533         return maff;
3534 }
3535
3536 __isl_give isl_aff *isl_aff_lift(__isl_take isl_aff *aff)
3537 {
3538         aff = isl_aff_cow(aff);
3539         if (!aff)
3540                 return NULL;
3541
3542         aff->ls = isl_local_space_lift(aff->ls);
3543         if (!aff->ls)
3544                 return isl_aff_free(aff);
3545
3546         return aff;
3547 }
3548
3549 /* Lift "maff" to a space with extra dimensions such that the result
3550  * has no more existentially quantified variables.
3551  * If "ls" is not NULL, then *ls is assigned the local space that lies
3552  * at the basis of the lifting applied to "maff".
3553  */
3554 __isl_give isl_multi_aff *isl_multi_aff_lift(__isl_take isl_multi_aff *maff,
3555         __isl_give isl_local_space **ls)
3556 {
3557         int i;
3558         isl_space *space;
3559         unsigned n_div;
3560
3561         if (ls)
3562                 *ls = NULL;
3563
3564         if (!maff)
3565                 return NULL;
3566
3567         if (maff->n == 0) {
3568                 if (ls) {
3569                         isl_space *space = isl_multi_aff_get_domain_space(maff);
3570                         *ls = isl_local_space_from_space(space);
3571                         if (!*ls)
3572                                 return isl_multi_aff_free(maff);
3573                 }
3574                 return maff;
3575         }
3576
3577         maff = isl_multi_aff_cow(maff);
3578         maff = isl_multi_aff_align_divs(maff);
3579         if (!maff)
3580                 return NULL;
3581
3582         n_div = isl_aff_dim(maff->p[0], isl_dim_div);
3583         space = isl_multi_aff_get_space(maff);
3584         space = isl_space_lift(isl_space_domain(space), n_div);
3585         space = isl_space_extend_domain_with_range(space,
3586                                                 isl_multi_aff_get_space(maff));
3587         if (!space)
3588                 return isl_multi_aff_free(maff);
3589         isl_space_free(maff->space);
3590         maff->space = space;
3591
3592         if (ls) {
3593                 *ls = isl_aff_get_domain_local_space(maff->p[0]);
3594                 if (!*ls)
3595                         return isl_multi_aff_free(maff);
3596         }
3597
3598         for (i = 0; i < maff->n; ++i) {
3599                 maff->p[i] = isl_aff_lift(maff->p[i]);
3600                 if (!maff->p[i])
3601                         goto error;
3602         }
3603
3604         return maff;
3605 error:
3606         if (ls)
3607                 isl_local_space_free(*ls);
3608         return isl_multi_aff_free(maff);
3609 }
3610
3611
3612 /* Extract an isl_pw_aff corresponding to output dimension "pos" of "pma".
3613  */
3614 __isl_give isl_pw_aff *isl_pw_multi_aff_get_pw_aff(
3615         __isl_keep isl_pw_multi_aff *pma, int pos)
3616 {
3617         int i;
3618         int n_out;
3619         isl_space *space;
3620         isl_pw_aff *pa;
3621
3622         if (!pma)
3623                 return NULL;
3624
3625         n_out = isl_pw_multi_aff_dim(pma, isl_dim_out);
3626         if (pos < 0 || pos >= n_out)
3627                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
3628                         "index out of bounds", return NULL);
3629
3630         space = isl_pw_multi_aff_get_space(pma);
3631         space = isl_space_drop_dims(space, isl_dim_out,
3632                                     pos + 1, n_out - pos - 1);
3633         space = isl_space_drop_dims(space, isl_dim_out, 0, pos);
3634
3635         pa = isl_pw_aff_alloc_size(space, pma->n);
3636         for (i = 0; i < pma->n; ++i) {
3637                 isl_aff *aff;
3638                 aff = isl_multi_aff_get_aff(pma->p[i].maff, pos);
3639                 pa = isl_pw_aff_add_piece(pa, isl_set_copy(pma->p[i].set), aff);
3640         }
3641
3642         return pa;
3643 }
3644
3645 /* Return an isl_pw_multi_aff with the given "set" as domain and
3646  * an unnamed zero-dimensional range.
3647  */
3648 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_from_domain(
3649         __isl_take isl_set *set)
3650 {
3651         isl_multi_aff *ma;
3652         isl_space *space;
3653
3654         space = isl_set_get_space(set);
3655         space = isl_space_from_domain(space);
3656         ma = isl_multi_aff_zero(space);
3657         return isl_pw_multi_aff_alloc(set, ma);
3658 }
3659
3660 /* Add an isl_pw_multi_aff with the given "set" as domain and
3661  * an unnamed zero-dimensional range to *user.
3662  */
3663 static int add_pw_multi_aff_from_domain(__isl_take isl_set *set, void *user)
3664 {
3665         isl_union_pw_multi_aff **upma = user;
3666         isl_pw_multi_aff *pma;
3667
3668         pma = isl_pw_multi_aff_from_domain(set);
3669         *upma = isl_union_pw_multi_aff_add_pw_multi_aff(*upma, pma);
3670
3671         return 0;
3672 }
3673
3674 /* Return an isl_union_pw_multi_aff with the given "uset" as domain and
3675  * an unnamed zero-dimensional range.
3676  */
3677 __isl_give isl_union_pw_multi_aff *isl_union_pw_multi_aff_from_domain(
3678         __isl_take isl_union_set *uset)
3679 {
3680         isl_space *space;
3681         isl_union_pw_multi_aff *upma;
3682
3683         if (!uset)
3684                 return NULL;
3685
3686         space = isl_union_set_get_space(uset);
3687         upma = isl_union_pw_multi_aff_empty(space);
3688
3689         if (isl_union_set_foreach_set(uset,
3690                                     &add_pw_multi_aff_from_domain, &upma) < 0)
3691                 goto error;
3692
3693         isl_union_set_free(uset);
3694         return upma;
3695 error:
3696         isl_union_set_free(uset);
3697         isl_union_pw_multi_aff_free(upma);
3698         return NULL;
3699 }
3700
3701 /* Convert "pma" to an isl_map and add it to *umap.
3702  */
3703 static int map_from_pw_multi_aff(__isl_take isl_pw_multi_aff *pma, void *user)
3704 {
3705         isl_union_map **umap = user;
3706         isl_map *map;
3707
3708         map = isl_map_from_pw_multi_aff(pma);
3709         *umap = isl_union_map_add_map(*umap, map);
3710
3711         return 0;
3712 }
3713
3714 /* Construct a union map mapping the domain of the union
3715  * piecewise multi-affine expression to its range, with each dimension
3716  * in the range equated to the corresponding affine expression on its cell.
3717  */
3718 __isl_give isl_union_map *isl_union_map_from_union_pw_multi_aff(
3719         __isl_take isl_union_pw_multi_aff *upma)
3720 {
3721         isl_space *space;
3722         isl_union_map *umap;
3723
3724         if (!upma)
3725                 return NULL;
3726
3727         space = isl_union_pw_multi_aff_get_space(upma);
3728         umap = isl_union_map_empty(space);
3729
3730         if (isl_union_pw_multi_aff_foreach_pw_multi_aff(upma,
3731                                         &map_from_pw_multi_aff, &umap) < 0)
3732                 goto error;
3733
3734         isl_union_pw_multi_aff_free(upma);
3735         return umap;
3736 error:
3737         isl_union_pw_multi_aff_free(upma);
3738         isl_union_map_free(umap);
3739         return NULL;
3740 }
3741
3742 /* Local data for bin_entry and the callback "fn".
3743  */
3744 struct isl_union_pw_multi_aff_bin_data {
3745         isl_union_pw_multi_aff *upma2;
3746         isl_union_pw_multi_aff *res;
3747         isl_pw_multi_aff *pma;
3748         int (*fn)(void **entry, void *user);
3749 };
3750
3751 /* Given an isl_pw_multi_aff from upma1, store it in data->pma
3752  * and call data->fn for each isl_pw_multi_aff in data->upma2.
3753  */
3754 static int bin_entry(void **entry, void *user)
3755 {
3756         struct isl_union_pw_multi_aff_bin_data *data = user;
3757         isl_pw_multi_aff *pma = *entry;
3758
3759         data->pma = pma;
3760         if (isl_hash_table_foreach(data->upma2->dim->ctx, &data->upma2->table,
3761                                    data->fn, data) < 0)
3762                 return -1;
3763
3764         return 0;
3765 }
3766
3767 /* Call "fn" on each pair of isl_pw_multi_affs in "upma1" and "upma2".
3768  * The isl_pw_multi_aff from upma1 is stored in data->pma (where data is
3769  * passed as user field) and the isl_pw_multi_aff from upma2 is available
3770  * as *entry.  The callback should adjust data->res if desired.
3771  */
3772 static __isl_give isl_union_pw_multi_aff *bin_op(
3773         __isl_take isl_union_pw_multi_aff *upma1,
3774         __isl_take isl_union_pw_multi_aff *upma2,
3775         int (*fn)(void **entry, void *user))
3776 {
3777         isl_space *space;
3778         struct isl_union_pw_multi_aff_bin_data data = { NULL, NULL, NULL, fn };
3779
3780         space = isl_union_pw_multi_aff_get_space(upma2);
3781         upma1 = isl_union_pw_multi_aff_align_params(upma1, space);
3782         space = isl_union_pw_multi_aff_get_space(upma1);
3783         upma2 = isl_union_pw_multi_aff_align_params(upma2, space);
3784
3785         if (!upma1 || !upma2)
3786                 goto error;
3787
3788         data.upma2 = upma2;
3789         data.res = isl_union_pw_multi_aff_alloc(isl_space_copy(upma1->dim),
3790                                        upma1->table.n);
3791         if (isl_hash_table_foreach(upma1->dim->ctx, &upma1->table,
3792                                    &bin_entry, &data) < 0)
3793                 goto error;
3794
3795         isl_union_pw_multi_aff_free(upma1);
3796         isl_union_pw_multi_aff_free(upma2);
3797         return data.res;
3798 error:
3799         isl_union_pw_multi_aff_free(upma1);
3800         isl_union_pw_multi_aff_free(upma2);
3801         isl_union_pw_multi_aff_free(data.res);
3802         return NULL;
3803 }
3804
3805 /* Given two isl_multi_affs A -> B and C -> D,
3806  * construct an isl_multi_aff (A * C) -> (B, D).
3807  */
3808 __isl_give isl_multi_aff *isl_multi_aff_flat_range_product(
3809         __isl_take isl_multi_aff *ma1, __isl_take isl_multi_aff *ma2)
3810 {
3811         int i, n1, n2;
3812         isl_aff *aff;
3813         isl_space *space;
3814         isl_multi_aff *res;
3815
3816         if (!ma1 || !ma2)
3817                 goto error;
3818
3819         space = isl_space_range_product(isl_multi_aff_get_space(ma1),
3820                                         isl_multi_aff_get_space(ma2));
3821         space = isl_space_flatten_range(space);
3822         res = isl_multi_aff_alloc(space);
3823
3824         n1 = isl_multi_aff_dim(ma1, isl_dim_out);
3825         n2 = isl_multi_aff_dim(ma2, isl_dim_out);
3826
3827         for (i = 0; i < n1; ++i) {
3828                 aff = isl_multi_aff_get_aff(ma1, i);
3829                 res = isl_multi_aff_set_aff(res, i, aff);
3830         }
3831
3832         for (i = 0; i < n2; ++i) {
3833                 aff = isl_multi_aff_get_aff(ma2, i);
3834                 res = isl_multi_aff_set_aff(res, n1 + i, aff);
3835         }
3836
3837         isl_multi_aff_free(ma1);
3838         isl_multi_aff_free(ma2);
3839         return res;
3840 error:
3841         isl_multi_aff_free(ma1);
3842         isl_multi_aff_free(ma2);
3843         return NULL;
3844 }
3845
3846 /* Given two aligned isl_pw_multi_affs A -> B and C -> D,
3847  * construct an isl_pw_multi_aff (A * C) -> (B, D).
3848  */
3849 static __isl_give isl_pw_multi_aff *pw_multi_aff_flat_range_product(
3850         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
3851 {
3852         isl_space *space;
3853
3854         space = isl_space_range_product(isl_pw_multi_aff_get_space(pma1),
3855                                         isl_pw_multi_aff_get_space(pma2));
3856         space = isl_space_flatten_range(space);
3857         return isl_pw_multi_aff_on_shared_domain_in(pma1, pma2, space,
3858                                             &isl_multi_aff_flat_range_product);
3859 }
3860
3861 /* Given two isl_pw_multi_affs A -> B and C -> D,
3862  * construct an isl_pw_multi_aff (A * C) -> (B, D).
3863  */
3864 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_flat_range_product(
3865         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
3866 {
3867         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
3868                                             &pw_multi_aff_flat_range_product);
3869 }
3870
3871 /* If data->pma and *entry have the same domain space, then compute
3872  * their flat range product and the result to data->res.
3873  */
3874 static int flat_range_product_entry(void **entry, void *user)
3875 {
3876         struct isl_union_pw_multi_aff_bin_data *data = user;
3877         isl_pw_multi_aff *pma2 = *entry;
3878
3879         if (!isl_space_tuple_match(data->pma->dim, isl_dim_in,
3880                                  pma2->dim, isl_dim_in))
3881                 return 0;
3882
3883         pma2 = isl_pw_multi_aff_flat_range_product(
3884                                         isl_pw_multi_aff_copy(data->pma),
3885                                         isl_pw_multi_aff_copy(pma2));
3886
3887         data->res = isl_union_pw_multi_aff_add_pw_multi_aff(data->res, pma2);
3888
3889         return 0;
3890 }
3891
3892 /* Given two isl_union_pw_multi_affs A -> B and C -> D,
3893  * construct an isl_union_pw_multi_aff (A * C) -> (B, D).
3894  */
3895 __isl_give isl_union_pw_multi_aff *isl_union_pw_multi_aff_flat_range_product(
3896         __isl_take isl_union_pw_multi_aff *upma1,
3897         __isl_take isl_union_pw_multi_aff *upma2)
3898 {
3899         return bin_op(upma1, upma2, &flat_range_product_entry);
3900 }
3901
3902 /* Replace the affine expressions at position "pos" in "pma" by "pa".
3903  * The parameters are assumed to have been aligned.
3904  *
3905  * The implementation essentially performs an isl_pw_*_on_shared_domain,
3906  * except that it works on two different isl_pw_* types.
3907  */
3908 static __isl_give isl_pw_multi_aff *pw_multi_aff_set_pw_aff(
3909         __isl_take isl_pw_multi_aff *pma, unsigned pos,
3910         __isl_take isl_pw_aff *pa)
3911 {
3912         int i, j, n;
3913         isl_pw_multi_aff *res = NULL;
3914
3915         if (!pma || !pa)
3916                 goto error;
3917
3918         if (!isl_space_tuple_match(pma->dim, isl_dim_in, pa->dim, isl_dim_in))
3919                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
3920                         "domains don't match", goto error);
3921         if (pos >= isl_pw_multi_aff_dim(pma, isl_dim_out))
3922                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
3923                         "index out of bounds", goto error);
3924
3925         n = pma->n * pa->n;
3926         res = isl_pw_multi_aff_alloc_size(isl_pw_multi_aff_get_space(pma), n);
3927
3928         for (i = 0; i < pma->n; ++i) {
3929                 for (j = 0; j < pa->n; ++j) {
3930                         isl_set *common;
3931                         isl_multi_aff *res_ij;
3932                         int empty;
3933
3934                         common = isl_set_intersect(isl_set_copy(pma->p[i].set),
3935                                                    isl_set_copy(pa->p[j].set));
3936                         empty = isl_set_plain_is_empty(common);
3937                         if (empty < 0 || empty) {
3938                                 isl_set_free(common);
3939                                 if (empty < 0)
3940                                         goto error;
3941                                 continue;
3942                         }
3943
3944                         res_ij = isl_multi_aff_set_aff(
3945                                         isl_multi_aff_copy(pma->p[i].maff), pos,
3946                                         isl_aff_copy(pa->p[j].aff));
3947                         res_ij = isl_multi_aff_gist(res_ij,
3948                                         isl_set_copy(common));
3949
3950                         res = isl_pw_multi_aff_add_piece(res, common, res_ij);
3951                 }
3952         }
3953
3954         isl_pw_multi_aff_free(pma);
3955         isl_pw_aff_free(pa);
3956         return res;
3957 error:
3958         isl_pw_multi_aff_free(pma);
3959         isl_pw_aff_free(pa);
3960         return isl_pw_multi_aff_free(res);
3961 }
3962
3963 /* Replace the affine expressions at position "pos" in "pma" by "pa".
3964  */
3965 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_set_pw_aff(
3966         __isl_take isl_pw_multi_aff *pma, unsigned pos,
3967         __isl_take isl_pw_aff *pa)
3968 {
3969         if (!pma || !pa)
3970                 goto error;
3971         if (isl_space_match(pma->dim, isl_dim_param, pa->dim, isl_dim_param))
3972                 return pw_multi_aff_set_pw_aff(pma, pos, pa);
3973         if (!isl_space_has_named_params(pma->dim) ||
3974             !isl_space_has_named_params(pa->dim))
3975                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
3976                         "unaligned unnamed parameters", goto error);
3977         pma = isl_pw_multi_aff_align_params(pma, isl_pw_aff_get_space(pa));
3978         pa = isl_pw_aff_align_params(pa, isl_pw_multi_aff_get_space(pma));
3979         return pw_multi_aff_set_pw_aff(pma, pos, pa);
3980 error:
3981         isl_pw_multi_aff_free(pma);
3982         isl_pw_aff_free(pa);
3983         return NULL;
3984 }