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